From 1ebf54929e515201289d587b862f1ec248e8856d Mon Sep 17 00:00:00 2001
From: awb9691 <abby.broscius@gmail.com>
Date: Mon, 19 Jun 2023 11:01:17 -0400
Subject: [PATCH 1/7] add comments to DirectedSphereExclusion, predict_radius

---
 DiverseSelector/methods/partition.py | 16 ++++++++++++----
 DiverseSelector/methods/utils.py     | 20 ++++++++++++++++----
 2 files changed, 28 insertions(+), 8 deletions(-)

diff --git a/DiverseSelector/methods/partition.py b/DiverseSelector/methods/partition.py
index 6456e111..89122da9 100644
--- a/DiverseSelector/methods/partition.py
+++ b/DiverseSelector/methods/partition.py
@@ -48,11 +48,15 @@ class DirectedSphereExclusion(SelectionBase):
     Starting point is chosen as the reference point and not included in the selected molecules. The
     distance of each point is calculated to the reference point and the points are then sorted based
     on the ascending order of distances. The points are then evaluated in their sorted order, and
-    are selected if their distance to all the other selected points is at least r away. Euclidian
+    are selected if their distance to all the other selected points is at least r away. Euclidean
     distance is used by default and the r value is automatically generated if not passed to satisfy
     the number of molecules requested.
 
-    Adapted from https://doi.org/10.1021/ci025554v
+    Notes
+    -----
+    Gobbi, A., and Lee, M.-L. (2002). DISE: directed sphere exclusion.
+    Journal of Chemical Information and Computer Sciences,
+    43(1), 317–323. https://doi.org/10.1021/ci025554v
     """
 
     def __init__(self, r=None, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42):
@@ -106,25 +110,29 @@ def algorithm(self, arr, uplimit):
         count = 0
         candidates = np.delete(np.arange(0, len(arr)), self.starting_idx)
         distances = []
+        # calculate distance from reference point to all data points
         for idx in candidates:
             ref_point = arr[self.starting_idx]
             data_point = arr[idx]
             distance = spatial.distance.minkowski(ref_point, data_point, p=self.p)
             distances.append((distance, idx))
+        # order points by distance from reference
         distances.sort()
         order = [idx for dist, idx in distances]
-
+        # initialize search variables
         kdtree = spatial.KDTree(arr)
         bv = bitarray.bitarray(len(arr))
         bv[:] = 0
         bv[self.starting_idx] = 1
-
+        # select points
         for idx in order:
             if not bv[idx]:
                 selected.append(idx)
                 count += 1
+                # finished selecting # of points required, return
                 if count > uplimit:
                     return selected
+                # eliminate points within radius from consideration
                 elim = kdtree.query_ball_point(arr[idx], self.r, eps=self.eps, p=self.p, workers=-1)
                 for index in elim:
                     bv[index] = 1
diff --git a/DiverseSelector/methods/utils.py b/DiverseSelector/methods/utils.py
index 7a5b1753..895ecf14 100644
--- a/DiverseSelector/methods/utils.py
+++ b/DiverseSelector/methods/utils.py
@@ -31,8 +31,7 @@
 
 
 def predict_radius(obj, arr, num_selected, cluster_ids=None):
-    """
-    Algorithm that uses sphere_exclusion for selecting points from cluster.
+    """Algorithm that uses sphere_exclusion for selecting points from cluster.
 
     Parameters
     ----------
@@ -50,43 +49,56 @@ def predict_radius(obj, arr, num_selected, cluster_ids=None):
     selected: list
         List of ids of selected molecules
     """
+    # set the limits on # of selected points according to the tolerance percentage
     error = num_selected * obj.tolerance / 100
     lowlimit = num_selected - error
     uplimit = num_selected + error
-
+    # sort into clusters if data is labelled
     if cluster_ids is not None:
         arr = arr[cluster_ids]
 
     original_r = None
+    # use specified radius if passed in
     if obj.r is not None:
         original_r = obj.r
+        # run the selection algorithm
         result = obj.algorithm(arr, uplimit)
-    # Use numpy.optimize.bisect instead
+    # calculate radius from largest feature values
     else:
         rg = max(np.ptp(arr, axis=0)) / num_selected * 3
         obj.r = rg
         result = obj.algorithm(arr, uplimit)
+
+    # correct number of points chosen, return selected
     if len(result) == num_selected:
         return result
 
+    # incorrect number of points chosen
     low = obj.r if len(result) > num_selected else 0
     high = obj.r if low == 0 else None
     bounds = [low, high]
     count = 0
     while (len(result) < lowlimit or len(result) > uplimit) and count < 10:
+        # too many points selected, make radius larger
         if bounds[1] is None:
             rg = bounds[0] * 2
+        # too few points selected, make radius smaller
         else:
             rg = (bounds[0] + bounds[1]) / 2
         obj.r = rg
+        # re-run selection with new radius
         result = obj.algorithm(arr, uplimit)
+
+        # adjust upper/lower bounds of radius size to fine tune
         if len(result) > num_selected:
             bounds[0] = rg
         else:
             bounds[1] = rg
         count += 1
+    # cannot find radius that produces desired number of selected points
     if count >= 10:
         print(f"Optimal radius finder failed to converge, selected {len(result)} molecules instead "
               f"of requested {num_selected}.")
+    # undo any changes to radius
     obj.r = original_r
     return result

From bb2415210fa9f41e297b0098548b7a2798594606 Mon Sep 17 00:00:00 2001
From: AWBroscius <abby.broscius@gmail.com>
Date: Fri, 23 Jun 2023 12:37:11 -0400
Subject: [PATCH 2/7] optimize Directed_Spheres, add comments, tests

---
 DiverseSelector/methods/partition.py          | 72 ++++++++++--------
 .../methods/tests/test_partition.py           | 73 +++++++++++++------
 2 files changed, 93 insertions(+), 52 deletions(-)

diff --git a/DiverseSelector/methods/partition.py b/DiverseSelector/methods/partition.py
index 89122da9..31a8d98b 100644
--- a/DiverseSelector/methods/partition.py
+++ b/DiverseSelector/methods/partition.py
@@ -26,6 +26,8 @@
 import math
 
 import bitarray
+import scipy.spatial
+
 from DiverseSelector.methods.base import SelectionBase
 from DiverseSelector.diversity import compute_diversity
 from DiverseSelector.methods.utils import predict_radius
@@ -45,7 +47,8 @@
 class DirectedSphereExclusion(SelectionBase):
     """Selecting points using Directed Sphere Exclusion algorithm.
 
-    Starting point is chosen as the reference point and not included in the selected molecules. The
+    Starting point is chosen as the reference point
+    and not included in the selected molecules. The
     distance of each point is calculated to the reference point and the points are then sorted based
     on the ascending order of distances. The points are then evaluated in their sorted order, and
     are selected if their distance to all the other selected points is at least r away. Euclidean
@@ -59,7 +62,7 @@ class DirectedSphereExclusion(SelectionBase):
     43(1), 317–323. https://doi.org/10.1021/ci025554v
     """
 
-    def __init__(self, r=None, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42):
+    def __init__(self, r=None, tolerance=0.05, eps=1e-8, p=2, start_id=0, random_seed=42):
         """
         Initializing class.
 
@@ -69,7 +72,7 @@ def __init__(self, r=None, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42
             Initial guess of radius for directed sphere exclusion algorithm, no points within r
             distance to an already selected point can be selected.
         tolerance: float
-            Percentage error of number of molecules actually selected from number of molecules
+            Percentage error of number of points actually selected from number of points
             requested.
         eps: float
             Approximate nearest neighbor search for eliminating close points. Branches of the tree
@@ -90,74 +93,85 @@ def __init__(self, r=None, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42
         self.starting_idx = start_id
         self.random_seed = random_seed
 
-    def algorithm(self, arr, uplimit):
+    def algorithm(self, x, uplimit):
         """
-        Directed sphere exclusion algorithm logic.
+        Directed sphere exclusion algorithm.
+
+        Given a reference point, sorts all points by distance to the reference point.
+        Then using a KDTree, the closest points are selected and a sphere
+        is built around the point. All points within the sphere are excluded
+        from the search. This process iterates until the number of selected
+        points is greater than `uplimit`, or the algorithm runs out of points
+        to select from.
 
         Parameters
         ----------
-        arr: np.ndarray
-            Coordinate array of points.
+        x: np.ndarray
+            Feature matrix.
         uplimit: int
             Maximum number of points to select.
 
         Returns
         -------
         selected: list
-            List of ids of selected molecules
+            List of ids of selected points.
         """
         selected = []
         count = 0
-        candidates = np.delete(np.arange(0, len(arr)), self.starting_idx)
-        distances = []
         # calculate distance from reference point to all data points
-        for idx in candidates:
-            ref_point = arr[self.starting_idx]
-            data_point = arr[idx]
-            distance = spatial.distance.minkowski(ref_point, data_point, p=self.p)
-            distances.append((distance, idx))
+        ref_point = x[self.starting_idx]
+        distances = scipy.spatial.minkowski_distance(ref_point, x, p=self.p)
         # order points by distance from reference
-        distances.sort()
-        order = [idx for dist, idx in distances]
-        # initialize search variables
-        kdtree = spatial.KDTree(arr)
-        bv = bitarray.bitarray(len(arr))
+        order = np.argsort(distances)
+        # Construct KDTree to make it easier to search neighbors
+        kdtree = spatial.KDTree(x)
+        # bv tracks viable candidates
+        bv = bitarray.bitarray(len(x))
         bv[:] = 0
         bv[self.starting_idx] = 1
-        # select points
+        # select points based on closest to reference point
         for idx in order:
+            # If it isn't already part of any hyperspheres
             if not bv[idx]:
+                # Then select it to be part of it
                 selected.append(idx)
                 count += 1
                 # finished selecting # of points required, return
                 if count > uplimit:
                     return selected
-                # eliminate points within radius from consideration
-                elim = kdtree.query_ball_point(arr[idx], self.r, eps=self.eps, p=self.p, workers=-1)
+                # find all points now within radius of newly selected point
+                elim = kdtree.query_ball_point(x[idx], self.r, eps=self.eps, p=self.p, workers=-1)
+                # turn 'on' bits in bv to make for loop skip indices of eliminated points
+                #   eliminate points from selection
                 for index in elim:
                     bv[index] = 1
 
         return selected
 
-    def select_from_cluster(self, arr, num_selected, cluster_ids=None):
+    def select_from_cluster(self, x, num_selected, cluster_ids=None):
         """
         Algorithm that uses sphere_exclusion for selecting points from cluster.
 
         Parameters
         ----------
-        arr: np.ndarray
-            Coordinate array of points
+        x: np.ndarray
+            Feature points.
         num_selected: int
-            Number of molecules that need to be selected.
+            Number of points that need to be selected.
         cluster_ids: np.ndarray
-            Indices of molecules that form a cluster
+            Indices of points that form a cluster
 
         Returns
         -------
         selected: list
             List of ids of selected molecules
         """
-        return predict_radius(self, arr, num_selected, cluster_ids)
+        if x.shape[0] < num_selected:
+            raise RuntimeError(
+                f"The number of selected points {num_selected} is greater than the number of points"
+                f"provided {x.shape[0]}."
+            )
+        return predict_radius(self, x, num_selected, cluster_ids)
 
 
 class GridPartitioning(SelectionBase):
diff --git a/DiverseSelector/methods/tests/test_partition.py b/DiverseSelector/methods/tests/test_partition.py
index ced8f90a..a489b736 100644
--- a/DiverseSelector/methods/tests/test_partition.py
+++ b/DiverseSelector/methods/tests/test_partition.py
@@ -22,36 +22,63 @@
 
 """Test Partition-Based Selection Methods."""
 
+import numpy as np
 from DiverseSelector.methods.partition import DirectedSphereExclusion, GridPartitioning, Medoid
 from DiverseSelector.methods.tests.common import generate_synthetic_data
-from numpy.testing import assert_equal
+from numpy.testing import assert_equal, assert_raises
 
 
+def test_directed_sphere_num_selected_error():
+    """Test DirectedSphereExclusion error when too many points requested."""
+    x = np.array([[1, 9]]*100)
+    selector = DirectedSphereExclusion()
+    assert_raises(RuntimeError, selector.select, x, num_selected=105)
 
-def test_directedsphereexclusion():
-    """Testing DirectedSphereExclusion class."""
-    coords, _, _ = generate_synthetic_data(n_samples=100,
-                                           n_features=2,
-                                           n_clusters=1,
-                                           pairwise_dist=True,
-                                           metric="euclidean",
-                                           random_state=42)
 
-    coords_cluster, class_labels_cluster, _ = generate_synthetic_data(n_samples=100,
-                                                                      n_features=2,
-                                                                      n_clusters=3,
-                                                                      pairwise_dist=True,
-                                                                      metric="euclidean",
-                                                                      random_state=42)
-    selector = DirectedSphereExclusion()
-    selected_ids = selector.select(arr=coords_cluster, size=12, labels=class_labels_cluster)
-    # make sure all the selected indices are the same with expectation
-    assert_equal(selected_ids, [95, 14, 88, 84, 76, 68, 93, 50, 29, 19, 54])
+def test_directed_sphere_same_number_of_pts():
+    """Test DirectSphereExclusion with `num_selected` = number of points in dataset."""
+    # (0,0) as the reference point
+    x = np.array([[0,0],[0,1],[0,2],[0,3]])
+    selector = DirectedSphereExclusion(r=1, tolerance=0)
+    selected = selector.select(arr=x, num_selected=3)
+    expected = [1,2,3]
+    assert_equal(selected, expected)
+    assert_equal(selector.r, 0.5)
 
-    selector = DirectedSphereExclusion()
-    selected_ids = selector.select(arr=coords, size=12)
-    # make sure all the selected indices are the same with expectation
-    assert_equal(selected_ids, [17, 92, 64, 6, 12, 76, 10, 87, 73, 66, 11, 57])
+
+def test_directed_sphere_exclusion_select_more_number_of_pts():
+    """Test DirectSphereExclusion on points on the line with `num_selected` < number of points in dataset."""
+    # (0,0) as the reference point
+    x = np.array([[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [0, 6]])
+    selector = DirectedSphereExclusion(r=0.5, tolerance=0)
+    selected = selector.select(arr=x, num_selected=3)
+    expected = [1, 3, 5]
+    assert_equal(selected, expected)
+    assert_equal(selector.r, 1.0)
+
+
+def test_directed_sphere_exclusion_on_line_with_():
+    """Test Direct Sphere Exclusion on points on line with smaller distribution than the radius."""
+    # (0,0) as the reference point
+    x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
+                  [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5], [0, 6]])
+    selector = DirectedSphereExclusion(r=0.5, tolerance=0)
+    selected = selector.select(arr=x, num_selected=3)
+    expected = [1, 5, 9]
+    assert_equal(selected, expected)
+    assert_equal(selector.r, 1.0)
+
+
+def test_directed_sphere_on_line_with_larger_radius():
+    """Test Direct Sphere Exclusion on points on the line with a too large radius size."""
+    # (0,0) as the reference point
+    x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
+                  [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5]])
+    selector = DirectedSphereExclusion(r=2.0, tolerance=0)
+    selected = selector.select(arr=x, num_selected=3)
+    expected = [1, 5, 9]
+    assert_equal(selected, expected)
+    assert_equal(selector.r, 1.0)
 
 
 def test_gridpartitioning():

From 7250e73c67a248a1266c910ea78a48237a02dc56 Mon Sep 17 00:00:00 2001
From: AWBroscius <abby.broscius@gmail.com>
Date: Fri, 23 Jun 2023 12:38:11 -0400
Subject: [PATCH 3/7] rename predict_radius, add comments, docs

---
 DiverseSelector/methods/dissimilarity.py |  4 +-
 DiverseSelector/methods/utils.py         | 78 ++++++++++++++----------
 2 files changed, 48 insertions(+), 34 deletions(-)

diff --git a/DiverseSelector/methods/dissimilarity.py b/DiverseSelector/methods/dissimilarity.py
index b339241a..67a70592 100644
--- a/DiverseSelector/methods/dissimilarity.py
+++ b/DiverseSelector/methods/dissimilarity.py
@@ -23,7 +23,7 @@
 """Module for Dissimilarity-Based Selection Methods."""
 
 from DiverseSelector.methods.base import SelectionBase
-from DiverseSelector.methods.utils import predict_radius
+from DiverseSelector.methods.utils import optimize_radius
 import numpy as np
 from scipy import spatial
 
@@ -272,6 +272,6 @@ def select_from_cluster(self, arr, num_selected, cluster_ids=None):
         selected: list
             List of ids of selected molecules
         """
-        return predict_radius(self, arr, num_selected, cluster_ids)
+        return optimize_radius(self, arr, num_selected, cluster_ids)
 
 
diff --git a/DiverseSelector/methods/utils.py b/DiverseSelector/methods/utils.py
index 895ecf14..a68beae5 100644
--- a/DiverseSelector/methods/utils.py
+++ b/DiverseSelector/methods/utils.py
@@ -21,84 +21,98 @@
 # --
 
 """Module for Selection Utilities."""
+import warnings
 
 import numpy as np
 
 
 __all__ = [
-    "predict_radius",
+    "optimize_radius",
 ]
 
 
-def predict_radius(obj, arr, num_selected, cluster_ids=None):
-    """Algorithm that uses sphere_exclusion for selecting points from cluster.
+def optimize_radius(obj, x, num_selected, cluster_ids=None):
+    """Algorithm that uses sphere exclusion for selecting points from cluster.
+
+    Iteratively searches for the optimal radius to obtain the correct number
+    of selected points. If the radius cannot converge to return `num_selected` points,
+    the function returns the closest number of points to `num_selected` as possible.
 
     Parameters
     ----------
     obj: object
-        Instance of dissimilarity selection class
-    arr: np.ndarray
-        Coordinate array of points
+        Instance of 'DirectedSphereExclusion' or 'OptiSim' selection class.
+    x: np.ndarray
+        Feature array.
     num_selected: int
-        Number of molecules that need to be selected.
+        Number of points that need to be selected.
     cluster_ids: np.ndarray
-        Indices of molecules that form a cluster
+        Indices of points that form a cluster.
 
     Returns
     -------
     selected: list
-        List of ids of selected molecules
+        List of ids of selected points.
     """
+    if x.shape[0] < num_selected:
+        raise RuntimeError(
+            f"The number of selected points {num_selected} is greater than the number of points"
+            f"provided {x.shape[0]}."
+        )
     # set the limits on # of selected points according to the tolerance percentage
-    error = num_selected * obj.tolerance / 100
-    lowlimit = num_selected - error
-    uplimit = num_selected + error
+    error = num_selected * obj.tolerance
+    lowlimit = round(num_selected - error)
+    uplimit = round(num_selected + error)
     # sort into clusters if data is labelled
     if cluster_ids is not None:
-        arr = arr[cluster_ids]
+        x = x[cluster_ids]
 
-    original_r = None
     # use specified radius if passed in
     if obj.r is not None:
-        original_r = obj.r
         # run the selection algorithm
-        result = obj.algorithm(arr, uplimit)
+        result = obj.algorithm(x, uplimit)
     # calculate radius from largest feature values
     else:
-        rg = max(np.ptp(arr, axis=0)) / num_selected * 3
+        rg = max(np.ptp(x, axis=0)) / num_selected * 3
         obj.r = rg
-        result = obj.algorithm(arr, uplimit)
+        result = obj.algorithm(x, uplimit)
 
     # correct number of points chosen, return selected
     if len(result) == num_selected:
         return result
 
-    # incorrect number of points chosen
-    low = obj.r if len(result) > num_selected else 0
-    high = obj.r if low == 0 else None
-    bounds = [low, high]
-    count = 0
-    while (len(result) < lowlimit or len(result) > uplimit) and count < 10:
+    # if incorrect number of points chosen, then optimize radius
+    if len(result) > num_selected:
+        # Too many points, so the radius should be bigger.
+        bounds = [obj.r, np.inf]
+    else:
+        bounds = [0, obj.r]
+    niter = 0
+    print(f"Number of results {len(result)}, {lowlimit}, {uplimit}")
+    while (len(result) < lowlimit or len(result) > uplimit) and niter < 10:
         # too many points selected, make radius larger
-        if bounds[1] is None:
+        if bounds[1] == np.inf:
             rg = bounds[0] * 2
         # too few points selected, make radius smaller
         else:
             rg = (bounds[0] + bounds[1]) / 2
         obj.r = rg
         # re-run selection with new radius
-        result = obj.algorithm(arr, uplimit)
+        result = obj.algorithm(x, uplimit)
 
         # adjust upper/lower bounds of radius size to fine tune
         if len(result) > num_selected:
             bounds[0] = rg
         else:
             bounds[1] = rg
-        count += 1
+        print(f"Radius ", rg)
+        niter += 1
     # cannot find radius that produces desired number of selected points
-    if count >= 10:
-        print(f"Optimal radius finder failed to converge, selected {len(result)} molecules instead "
-              f"of requested {num_selected}.")
-    # undo any changes to radius
-    obj.r = original_r
+    if niter >= 10:
+        warnings.warn(
+            f"Optimal radius finder failed to converge, selected {len(result)} molecules instead "
+              f"of requested {num_selected}."
+        )
+
+    print(f"radius after optimization: ", obj.r)
     return result

From d8ef2cc0f61820f62c74e70066b9c22a75a9d4ae Mon Sep 17 00:00:00 2001
From: AWBroscius <abby.broscius@gmail.com>
Date: Fri, 23 Jun 2023 16:23:05 -0400
Subject: [PATCH 4/7] add n_iter as user specified argument

---
 DiverseSelector/methods/dissimilarity.py      |  5 +++-
 DiverseSelector/methods/partition.py          | 24 ++++++++++--------
 .../methods/tests/test_partition.py           |  8 +++---
 DiverseSelector/methods/utils.py              | 25 +++++++++----------
 4 files changed, 33 insertions(+), 29 deletions(-)

diff --git a/DiverseSelector/methods/dissimilarity.py b/DiverseSelector/methods/dissimilarity.py
index 67a70592..e1c160f4 100644
--- a/DiverseSelector/methods/dissimilarity.py
+++ b/DiverseSelector/methods/dissimilarity.py
@@ -171,7 +171,7 @@ class OptiSim(SelectionBase):
     Adapted from  https://doi.org/10.1021/ci970282v
     """
 
-    def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42):
+    def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42, n_iter=10):
         """
         Initializing class.
 
@@ -197,6 +197,8 @@ def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_s
             Index for the first point to be selected.
         random_seed: int
             Seed for random selection of points be evaluated.
+        n_iter: int
+            Number of iterations to execute when optimizing the size of exclusion radius. Default is 10.
         """
         self.r = r
         self.k = k
@@ -205,6 +207,7 @@ def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_s
         self.p = p
         self.start_id = start_id
         self.random_seed = random_seed
+        self.n_iter = n_iter
 
     def algorithm(self, arr, uplimit) -> list:
         """
diff --git a/DiverseSelector/methods/partition.py b/DiverseSelector/methods/partition.py
index 31a8d98b..d990d631 100644
--- a/DiverseSelector/methods/partition.py
+++ b/DiverseSelector/methods/partition.py
@@ -62,13 +62,13 @@ class DirectedSphereExclusion(SelectionBase):
     43(1), 317–323. https://doi.org/10.1021/ci025554v
     """
 
-    def __init__(self, r=None, tolerance=0.05, eps=1e-8, p=2, start_id=0, random_seed=42):
+    def __init__(self, r_0=None, tolerance=0.05, eps=1e-8, p=2, start_id=0, random_seed=42, n_iter=10):
         """
         Initializing class.
 
         Parameters
         ----------
-        r: float
+        r_0: float
             Initial guess of radius for directed sphere exclusion algorithm, no points within r
             distance to an already selected point can be selected.
         tolerance: float
@@ -82,16 +82,19 @@ def __init__(self, r=None, tolerance=0.05, eps=1e-8, p=2, start_id=0, random_see
             Which Minkowski p-norm to use. Should be in the range [1, inf]. A finite large p may
             cause a ValueError if overflow can occur.
         start_id: int
-            Index for the first point to be selected.
+            Index for the first point to be selected. Default is 0.
         random_seed: int
             Seed for random selection of points be evaluated.
+        n_iter: int
+            Number of iterations to execute when optimizing the size of exclusion radius. Default is 10.
         """
-        self.r = r
+        self.r = r_0
         self.tolerance = tolerance
         self.eps = eps
         self.p = p
         self.starting_idx = start_id
         self.random_seed = random_seed
+        self.n_iter = n_iter
 
     def algorithm(self, x, uplimit):
         """
@@ -116,8 +119,7 @@ def algorithm(self, x, uplimit):
         selected: list
             List of ids of selected points.
         """
-        selected = []
-        count = 0
+
         # calculate distance from reference point to all data points
         ref_point = x[self.starting_idx]
         distances = scipy.spatial.minkowski_distance(ref_point, x, p=self.p)
@@ -130,14 +132,14 @@ def algorithm(self, x, uplimit):
         bv[:] = 0
         bv[self.starting_idx] = 1
         # select points based on closest to reference point
+        selected = []
         for idx in order:
-            # If it isn't already part of any hyperspheres
-            if not bv[idx]:
-                # Then select it to be part of it
+            # If point isn't already part of any hyperspheres
+            if bv[idx] == 0:
+                # Then add point to selection
                 selected.append(idx)
-                count += 1
                 # finished selecting # of points required, return
-                if count > uplimit:
+                if len(selected) > uplimit:
                     return selected
                 # find all points now within radius of newly selected point
                 elim = kdtree.query_ball_point(x[idx], self.r, eps=self.eps, p=self.p, workers=-1)
diff --git a/DiverseSelector/methods/tests/test_partition.py b/DiverseSelector/methods/tests/test_partition.py
index a489b736..f8c8b572 100644
--- a/DiverseSelector/methods/tests/test_partition.py
+++ b/DiverseSelector/methods/tests/test_partition.py
@@ -39,7 +39,7 @@ def test_directed_sphere_same_number_of_pts():
     """Test DirectSphereExclusion with `num_selected` = number of points in dataset."""
     # (0,0) as the reference point
     x = np.array([[0,0],[0,1],[0,2],[0,3]])
-    selector = DirectedSphereExclusion(r=1, tolerance=0)
+    selector = DirectedSphereExclusion(r_0=1, tolerance=0)
     selected = selector.select(arr=x, num_selected=3)
     expected = [1,2,3]
     assert_equal(selected, expected)
@@ -50,7 +50,7 @@ def test_directed_sphere_exclusion_select_more_number_of_pts():
     """Test DirectSphereExclusion on points on the line with `num_selected` < number of points in dataset."""
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [0, 6]])
-    selector = DirectedSphereExclusion(r=0.5, tolerance=0)
+    selector = DirectedSphereExclusion(r_0=0.5, tolerance=0)
     selected = selector.select(arr=x, num_selected=3)
     expected = [1, 3, 5]
     assert_equal(selected, expected)
@@ -62,7 +62,7 @@ def test_directed_sphere_exclusion_on_line_with_():
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
                   [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5], [0, 6]])
-    selector = DirectedSphereExclusion(r=0.5, tolerance=0)
+    selector = DirectedSphereExclusion(r_0=0.5, tolerance=0)
     selected = selector.select(arr=x, num_selected=3)
     expected = [1, 5, 9]
     assert_equal(selected, expected)
@@ -74,7 +74,7 @@ def test_directed_sphere_on_line_with_larger_radius():
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
                   [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5]])
-    selector = DirectedSphereExclusion(r=2.0, tolerance=0)
+    selector = DirectedSphereExclusion(r_0=2.0, tolerance=0)
     selected = selector.select(arr=x, num_selected=3)
     expected = [1, 5, 9]
     assert_equal(selected, expected)
diff --git a/DiverseSelector/methods/utils.py b/DiverseSelector/methods/utils.py
index a68beae5..d1579e49 100644
--- a/DiverseSelector/methods/utils.py
+++ b/DiverseSelector/methods/utils.py
@@ -63,11 +63,8 @@ def optimize_radius(obj, x, num_selected, cluster_ids=None):
     error = num_selected * obj.tolerance
     lowlimit = round(num_selected - error)
     uplimit = round(num_selected + error)
-    # sort into clusters if data is labelled
-    if cluster_ids is not None:
-        x = x[cluster_ids]
 
-    # use specified radius if passed in
+    # use initial radius if passed in
     if obj.r is not None:
         # run the selection algorithm
         result = obj.algorithm(x, uplimit)
@@ -83,13 +80,15 @@ def optimize_radius(obj, x, num_selected, cluster_ids=None):
 
     # if incorrect number of points chosen, then optimize radius
     if len(result) > num_selected:
-        # Too many points, so the radius should be bigger.
+        # Too many points selected, make radius bigger.
         bounds = [obj.r, np.inf]
     else:
+        # Too few points selected, make radius smaller
         bounds = [0, obj.r]
-    niter = 0
+    n_iter = 0
     print(f"Number of results {len(result)}, {lowlimit}, {uplimit}")
-    while (len(result) < lowlimit or len(result) > uplimit) and niter < 10:
+
+    while (len(result) < lowlimit or len(result) > uplimit) and n_iter < obj.n_iter:
         # too many points selected, make radius larger
         if bounds[1] == np.inf:
             rg = bounds[0] * 2
@@ -102,17 +101,17 @@ def optimize_radius(obj, x, num_selected, cluster_ids=None):
 
         # adjust upper/lower bounds of radius size to fine tune
         if len(result) > num_selected:
+            # Too many points selected, raise lower bound
             bounds[0] = rg
         else:
+            # Too few points selected, lower upper bound
             bounds[1] = rg
-        print(f"Radius ", rg)
-        niter += 1
+        n_iter += 1
     # cannot find radius that produces desired number of selected points
-    if niter >= 10:
+    if n_iter >= obj.n_iter:
         warnings.warn(
-            f"Optimal radius finder failed to converge, selected {len(result)} molecules instead "
-              f"of requested {num_selected}."
+            f"Optimal radius finder failed to converge, selected {len(result)} points instead "
+            f"of requested {num_selected}."
         )
 
-    print(f"radius after optimization: ", obj.r)
     return result

From 59e2af6fa3f6a5671915a3fd6582e709c8fc1ffb Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Mon, 10 Jul 2023 14:05:02 -0400
Subject: [PATCH 5/7] Fix failed tests, wrong imports, and black updates

The `test_gridpartitioning` is commented, because the GridPartitioning
is not comptable with the new API.
---
 DiverseSelector/diversity.py                  | 38 +++++++----
 DiverseSelector/methods/partition.py          |  4 +-
 .../methods/tests/test_partition.py           | 66 +++++++++----------
 3 files changed, 59 insertions(+), 49 deletions(-)

diff --git a/DiverseSelector/diversity.py b/DiverseSelector/diversity.py
index e34b7b08..e0dbdff3 100644
--- a/DiverseSelector/diversity.py
+++ b/DiverseSelector/diversity.py
@@ -23,11 +23,11 @@
 
 """Molecule dataset diversity calculation module."""
 
+import warnings
 from typing import List
 
 import numpy as np
 from scipy.spatial.distance import euclidean
-import warnings
 
 __all__ = [
     "compute_diversity",
@@ -42,7 +42,7 @@
 
 def compute_diversity(
     features: np.array,
-    div_type: str = "total_diversity_volume",
+    div_type: str = "hypersphere_overlap_of_subset",
 ) -> float:
     """Compute diversity metrics.
 
@@ -53,7 +53,7 @@ def compute_diversity(
     div_type : str, optional
         Method of calculation diversity for a given molecule set, which
         includes "entropy", "logdet", "shannon_entropy", "wdud",
-        gini_coefficient" and "total_diversity_volume". Default is "total_diversity_volume".
+        gini_coefficient" and "hypersphere_overlap_of_subset".
     mols : List[rdkit.Chem.rdchem.Mol], optional
         List of RDKit molecule objects. This is only needed when using the
         "explicit_diversity_index" method. Default=None.
@@ -68,7 +68,7 @@ def compute_diversity(
         "logdet": logdet,
         "shannon_entropy": shannon_entropy,
         "wdud": wdud,
-        "total_diversity_volume": total_diversity_volume,
+        "hypersphere_overlap_of_subset": hypersphere_overlap_of_subset,
         "gini_coefficient": gini_coefficient,
     }
 
@@ -201,7 +201,10 @@ def shannon_entropy(x: np.ndarray) -> float:
         p_i = np.count_nonzero(x[:, i]) / num_mols
         # sum all non-zero terms
         if p_i == 0:
-            raise ValueError(f"Feature {i} has value 0 for all molecules. Remove extraneous feature from data set.")
+            raise ValueError(
+                f"Feature {i} has value 0 for all molecules."
+                "Remove extraneous feature from data set."
+            )
         h_x += (-1 * p_i) * np.log10(p_i)
     return h_x
 
@@ -246,7 +249,9 @@ def wdud(x: np.ndarray) -> float:
     min_x = np.min(x, axis=0)
     # Normalization of each feature to [0, 1]
     if np.any(np.abs(max_x - min_x) < 1e-30):
-        raise ValueError(f"One of the features is redundant and causes normalization to fail.")
+        raise ValueError(
+            f"One of the features is redundant and causes normalization to fail."
+        )
     x_norm = (x - min_x) / (max_x - min_x)
     ans = []  # store the Wasserstein distance for each feature
     for i in range(0, num_features):
@@ -254,7 +259,7 @@ def wdud(x: np.ndarray) -> float:
         y = np.sort(x_norm[:, i])
         # Round to the sixth decimal place and count number of unique elements
         #    to construct an accurate cumulative discrete distribution func \sum_{x <= y_{i + 1}} 1/k
-        y, counts = np.unique(np.round(x_norm[:,i], decimals=6), return_counts=True)
+        y, counts = np.unique(np.round(x_norm[:, i], decimals=6), return_counts=True)
         p = 0
         # Ignore 0 and because v_min= 0
         for j in range(1, len(counts)):
@@ -264,7 +269,7 @@ def wdud(x: np.ndarray) -> float:
             # Make a grid from yi1 to yi
             grid = np.linspace(yi1, yi, num=1000, endpoint=True)
             # Evaluate the integrand  |x - \sum_{x <= y_{i + 1}} 1/k|
-            p += counts[j-1]
+            p += counts[j - 1]
             integrand = np.abs(grid - p / num_mols)
             # Integrate using np.trapz
             wdu += np.trapz(y=integrand, x=grid)
@@ -311,23 +316,28 @@ def hypersphere_overlap_of_subset(lib: np.ndarray, x: np.array) -> float:
     min_x = np.min(lib, axis=0)
     # Normalization of each feature to [0, 1]
     if np.any(np.abs(max_x - min_x) < 1e-30):
-        raise ValueError(f"One of the features is redundant and causes normalization to fail.")
+        raise ValueError(
+            f"One of the features is redundant and causes normalization to fail."
+        )
     x_norm = (x - min_x) / (max_x - min_x)
     # r_o = hypersphere radius
     r_o = d * np.sqrt(1 / k)
     if r_o > 0.5:
-        warnings.warn(f"The number of molecules should be much larger"
-                      " than the number of features.")
+        warnings.warn(
+            f"The number of molecules should be much larger"
+            " than the number of features."
+        )
     g_s = 0
     edge = 0
-    lam = (d - 1.0) / d   # Lambda parameter controls edge penalty
+    # lambda parameter controls edge penalty
+    lam = (d - 1.0) / d
     # calculate overlap volume
     for i in range(0, (k - 1)):
         for j in range((i + 1), k):
             dist = np.linalg.norm(x_norm[i] - x_norm[j])
             # Overlap penalty
             if dist <= (2 * r_o):
-                with np.errstate(divide='ignore'):
+                with np.errstate(divide="ignore"):
                     # min(100) ignores the inf case with divide by zero
                     g_s += min(100, 2 * (r_o / dist) - 1)
         # Edge penalty: lambda (1 - \sum^d_j e_{ij} / (dr_0)
@@ -342,7 +352,7 @@ def hypersphere_overlap_of_subset(lib: np.ndarray, x: np.array) -> float:
             if dist > r_o:
                 dist = r_o
             edge_pen += dist
-        edge_pen /= (d * r_o)
+        edge_pen /= d * r_o
         # print("Should be positive value only", (1.0 - edge_pen))
         edge_pen = lam * (1.0 - edge_pen)
         edge += edge_pen
diff --git a/DiverseSelector/methods/partition.py b/DiverseSelector/methods/partition.py
index d990d631..2e3ddfd3 100644
--- a/DiverseSelector/methods/partition.py
+++ b/DiverseSelector/methods/partition.py
@@ -30,7 +30,7 @@
 
 from DiverseSelector.methods.base import SelectionBase
 from DiverseSelector.diversity import compute_diversity
-from DiverseSelector.methods.utils import predict_radius
+from DiverseSelector.methods.utils import optimize_radius
 import numpy as np
 from scipy import spatial
 from sklearn.decomposition import PCA
@@ -173,7 +173,7 @@ def select_from_cluster(self, x, num_selected, cluster_ids=None):
                 f"The number of selected points {num_selected} is greater than the number of points"
                 f"provided {x.shape[0]}."
             )
-        return predict_radius(self, x, num_selected, cluster_ids)
+        return optimize_radius(self, x, num_selected, cluster_ids)
 
 
 class GridPartitioning(SelectionBase):
diff --git a/DiverseSelector/methods/tests/test_partition.py b/DiverseSelector/methods/tests/test_partition.py
index f8c8b572..c27df3f5 100644
--- a/DiverseSelector/methods/tests/test_partition.py
+++ b/DiverseSelector/methods/tests/test_partition.py
@@ -28,30 +28,30 @@
 from numpy.testing import assert_equal, assert_raises
 
 
-def test_directed_sphere_num_selected_error():
+def test_directed_sphere_size_error():
     """Test DirectedSphereExclusion error when too many points requested."""
-    x = np.array([[1, 9]]*100)
+    x = np.array([[1, 9]] * 100)
     selector = DirectedSphereExclusion()
-    assert_raises(RuntimeError, selector.select, x, num_selected=105)
+    assert_raises(ValueError, selector.select, x, size=105)
 
 
 def test_directed_sphere_same_number_of_pts():
-    """Test DirectSphereExclusion with `num_selected` = number of points in dataset."""
+    """Test DirectSphereExclusion with `size` = number of points in dataset."""
     # (0,0) as the reference point
     x = np.array([[0,0],[0,1],[0,2],[0,3]])
     selector = DirectedSphereExclusion(r_0=1, tolerance=0)
-    selected = selector.select(arr=x, num_selected=3)
+    selected = selector.select(arr=x, size=3)
     expected = [1,2,3]
     assert_equal(selected, expected)
     assert_equal(selector.r, 0.5)
 
 
 def test_directed_sphere_exclusion_select_more_number_of_pts():
-    """Test DirectSphereExclusion on points on the line with `num_selected` < number of points in dataset."""
+    """Test DirectSphereExclusion on points on the line with `size` < number of points in dataset."""
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [0, 6]])
     selector = DirectedSphereExclusion(r_0=0.5, tolerance=0)
-    selected = selector.select(arr=x, num_selected=3)
+    selected = selector.select(arr=x, size=3)
     expected = [1, 3, 5]
     assert_equal(selected, expected)
     assert_equal(selector.r, 1.0)
@@ -63,7 +63,7 @@ def test_directed_sphere_exclusion_on_line_with_():
     x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
                   [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5], [0, 6]])
     selector = DirectedSphereExclusion(r_0=0.5, tolerance=0)
-    selected = selector.select(arr=x, num_selected=3)
+    selected = selector.select(arr=x, size=3)
     expected = [1, 5, 9]
     assert_equal(selected, expected)
     assert_equal(selector.r, 1.0)
@@ -75,36 +75,36 @@ def test_directed_sphere_on_line_with_larger_radius():
     x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
                   [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5]])
     selector = DirectedSphereExclusion(r_0=2.0, tolerance=0)
-    selected = selector.select(arr=x, num_selected=3)
+    selected = selector.select(arr=x, size=3)
     expected = [1, 5, 9]
     assert_equal(selected, expected)
     assert_equal(selector.r, 1.0)
 
 
-def test_gridpartitioning():
-    """Testing DirectedSphereExclusion class."""
-    coords, _, _ = generate_synthetic_data(n_samples=100,
-                                           n_features=2,
-                                           n_clusters=1,
-                                           pairwise_dist=True,
-                                           metric="euclidean",
-                                           random_state=42)
-
-    coords_cluster, class_labels_cluster, _ = generate_synthetic_data(n_samples=100,
-                                                                      n_features=2,
-                                                                      n_clusters=3,
-                                                                      pairwise_dist=True,
-                                                                      metric="euclidean",
-                                                                      random_state=42)
-    selector = GridPartitioning(cells=3)
-    selected_ids = selector.select(arr=coords_cluster, size=12, labels=class_labels_cluster)
-    # make sure all the selected indices are the same with expectation
-    assert_equal(selected_ids, [2, 25, 84, 56, 8, 70, 58, 78, 4, 46, 65, 29])
-
-    selector = GridPartitioning(cells=3)
-    selected_ids = selector.select(arr=coords, size=12)
-    # make sure all the selected indices are the same with expectation
-    assert_equal(selected_ids, [7, 55, 70, 57, 29, 91, 9, 65, 28, 11, 54, 88])
+# def test_gridpartitioning():
+#     """Testing DirectedSphereExclusion class."""
+#     coords, _, _ = generate_synthetic_data(n_samples=100,
+#                                            n_features=2,
+#                                            n_clusters=1,
+#                                            pairwise_dist=True,
+#                                            metric="euclidean",
+#                                            random_state=42)
+
+#     coords_cluster, class_labels_cluster, _ = generate_synthetic_data(n_samples=100,
+#                                                                       n_features=2,
+#                                                                       n_clusters=3,
+#                                                                       pairwise_dist=True,
+#                                                                       metric="euclidean",
+#                                                                       random_state=42)
+#     selector = GridPartitioning(cells=3)
+#     selected_ids = selector.select(arr=coords_cluster, size=12, labels=class_labels_cluster)
+#     # make sure all the selected indices are the same with expectation
+#     assert_equal(selected_ids, [2, 25, 84, 56, 8, 70, 58, 78, 4, 46, 65, 29])
+
+#     selector = GridPartitioning(cells=3)
+#     selected_ids = selector.select(arr=coords, size=12)
+#     # make sure all the selected indices are the same with expectation
+#     assert_equal(selected_ids, [7, 55, 70, 57, 29, 91, 9, 65, 28, 11, 54, 88])
 
 
 def test_medoid():

From d5c1fbc6576edc0c59c15716b58a15f43b841d32 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Fri, 14 Jul 2023 09:06:08 -0400
Subject: [PATCH 6/7] Finalize DirectedSphereExclusion implementation

---
 DiverseSelector/methods/dissimilarity.py      |   6 +-
 DiverseSelector/methods/partition.py          | 165 +++++++++---------
 .../methods/tests/test_partition.py           |   8 +-
 DiverseSelector/methods/utils.py              |  98 +++++------
 4 files changed, 136 insertions(+), 141 deletions(-)

diff --git a/DiverseSelector/methods/dissimilarity.py b/DiverseSelector/methods/dissimilarity.py
index e1c160f4..f7caa5b1 100644
--- a/DiverseSelector/methods/dissimilarity.py
+++ b/DiverseSelector/methods/dissimilarity.py
@@ -171,7 +171,7 @@ class OptiSim(SelectionBase):
     Adapted from  https://doi.org/10.1021/ci970282v
     """
 
-    def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_seed=42, n_iter=10):
+    def __init__(self, r=None, k=10, tol=5.0, eps=0, p=2, start_id=0, random_seed=42, n_iter=10):
         """
         Initializing class.
 
@@ -183,7 +183,7 @@ def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_s
         k: int
             Amount of points to add to subsample before selecting one of the points with the
             greatest minimum distance to the previously selected points.
-        tolerance: float
+        tol: float
             Percentage error of number of molecules actually selected from number of molecules
             requested.
         eps: float
@@ -202,7 +202,7 @@ def __init__(self, r=None, k=10, tolerance=5.0, eps=0, p=2, start_id=0, random_s
         """
         self.r = r
         self.k = k
-        self.tolerance = tolerance
+        self.tol = tol
         self.eps = eps
         self.p = p
         self.start_id = start_id
diff --git a/DiverseSelector/methods/partition.py b/DiverseSelector/methods/partition.py
index 2e3ddfd3..8744ac4f 100644
--- a/DiverseSelector/methods/partition.py
+++ b/DiverseSelector/methods/partition.py
@@ -45,15 +45,22 @@
 
 
 class DirectedSphereExclusion(SelectionBase):
-    """Selecting points using Directed Sphere Exclusion algorithm.
-
-    Starting point is chosen as the reference point
-    and not included in the selected molecules. The
-    distance of each point is calculated to the reference point and the points are then sorted based
-    on the ascending order of distances. The points are then evaluated in their sorted order, and
-    are selected if their distance to all the other selected points is at least r away. Euclidean
-    distance is used by default and the r value is automatically generated if not passed to satisfy
-    the number of molecules requested.
+    """Select samples using Directed Sphere Exclusion (DISE) algorithm.
+
+    In a nutshell, this algorithm iteratively excludes any sample within a given radius from
+    any already selected sample. The radius of the exclusion sphere is an adjustable parameter.
+    Compared to Sphere Exclusion algorithm, the Directed Sphere Exclusion algorithm achieves a
+    more evenly distributed subset selection by abandoning the random selection approach and
+    instead imposing a directed selection.
+
+    Reference sample is chosen based on the `ref_index`, which is excluded from the selected
+    subset. All samples are sorted (ascending order) based on their Minkowski p-norm distance
+    from the reference sample. Looping through sorted samples, the sample is selected if it is
+    not already excluded. If selected, all its neighboring samples within a sphere of radius r
+    (i.e., exclusion sphere) are excluded from being selected. When the selected number of points
+    is greater than specified subset `size`, the selection process terminates. The `r0` is used
+    as the initial radius of exclusion sphere, however, it is optimized to select the desired
+    number of samples.
 
     Notes
     -----
@@ -62,118 +69,106 @@ class DirectedSphereExclusion(SelectionBase):
     43(1), 317–323. https://doi.org/10.1021/ci025554v
     """
 
-    def __init__(self, r_0=None, tolerance=0.05, eps=1e-8, p=2, start_id=0, random_seed=42, n_iter=10):
-        """
-        Initializing class.
+    def __init__(self, r0=None, ref_index=0, p=2.0, eps=0.0, tol=0.05, n_iter=10, random_seed=42):
+        """Initialize class.
 
         Parameters
         ----------
-        r_0: float
-            Initial guess of radius for directed sphere exclusion algorithm, no points within r
-            distance to an already selected point can be selected.
-        tolerance: float
-            Percentage error of number of points actually selected from number of points
-            requested.
-        eps: float
-            Approximate nearest neighbor search for eliminating close points. Branches of the tree
-            are not explored if their nearest points are further than r / (1 + eps), and branches
-            are added in bulk if their furthest points are nearer than r * (1 + eps).
-        p: float
-            Which Minkowski p-norm to use. Should be in the range [1, inf]. A finite large p may
-            cause a ValueError if overflow can occur.
-        start_id: int
-            Index for the first point to be selected. Default is 0.
-        random_seed: int
+        r0: float, optional
+            Initial guess for radius of the exclusion sphere.
+        ref_index: int, optional
+            Index of the reference sample to start the selection algorithm from.
+            This sample is not included in the selected subset.
+        p: float, optional
+            Which Minkowski p-norm to use. The values of `p` should be within [1, inf].
+            A finite large p may cause a ValueError if overflow can occur.
+        eps: float, optional
+            Approximate nearest neighbor search used in `KDTree.query_ball_tree`.
+            Branches of the tree are not explored if their nearest points are further than
+            r/(1+eps), and branches are added in bulk if their furthest points are nearer than
+            r * (1+eps). eps has to be non-negative.
+        tol: float, optional
+            Percentage error of number of samples actually selected from number of samples requested.
+        n_iter: int, optional
+            Number of iterations for optimizing the radius of exclusion sphere.
+        random_seed: int, optional
             Seed for random selection of points be evaluated.
-        n_iter: int
-            Number of iterations to execute when optimizing the size of exclusion radius. Default is 10.
         """
-        self.r = r_0
-        self.tolerance = tolerance
-        self.eps = eps
+        self.r = r0
+        self.ref_index = ref_index
         self.p = p
-        self.starting_idx = start_id
-        self.random_seed = random_seed
+        self.eps = eps
+        self.tol = tol
         self.n_iter = n_iter
+        self.random_seed = random_seed
 
-    def algorithm(self, x, uplimit):
-        """
-        Directed sphere exclusion algorithm.
-
-        Given a reference point, sorts all points by distance to the reference point.
-        Then using a KDTree, the closest points are selected and a sphere
-        is built around the point. All points within the sphere are excluded
-        from the search. This process iterates until the number of selected
-        points is greater than `uplimit`, or the algorithm runs out of points
-        to select from.
+    def algorithm(self, X, max_size):
+        """Return selected samples based on directed sphere exclusion algorithm.
 
         Parameters
         ----------
-        x: np.ndarray
-            Feature matrix.
-        uplimit: int
-            Maximum number of points to select.
+        X: ndarray of shape (n_samples, n_features)
+           Feature matrix of `n_samples` samples in `n_features` dimensional space.
+        max_size: int
+            Maximum number of samples to select.
 
         Returns
         -------
         selected: list
-            List of ids of selected points.
+            List of indices of selected samples.
         """
 
-        # calculate distance from reference point to all data points
-        ref_point = x[self.starting_idx]
-        distances = scipy.spatial.minkowski_distance(ref_point, x, p=self.p)
-        # order points by distance from reference
-        order = np.argsort(distances)
-        # Construct KDTree to make it easier to search neighbors
-        kdtree = spatial.KDTree(x)
-        # bv tracks viable candidates
-        bv = bitarray.bitarray(len(x))
-        bv[:] = 0
-        bv[self.starting_idx] = 1
-        # select points based on closest to reference point
+        # calculate distance of all samples from reference sample; distance is a (n_samples,) array
+        distances = scipy.spatial.minkowski_distance(X[self.ref_index], X, p=self.p)
+        # get sorted index of samples based on their distance from reference (closest to farthest)
+        index_sorted = np.argsort(distances)
+        # construct KDTree for quick nearest-neighbor lookup
+        kdtree = spatial.KDTree(X)
+
+        # construct bitarray to track selected samples (1 means exclude)
+        bv = bitarray.bitarray(list(np.zeros(len(X), dtype=int)))
+        bv[self.ref_index] = 1
+
         selected = []
-        for idx in order:
-            # If point isn't already part of any hyperspheres
+        for idx in index_sorted:
+            # select sample if it is not already excluded from consideration
+            # indexing a single item of a bitarray will always return an integer
             if bv[idx] == 0:
-                # Then add point to selection
                 selected.append(idx)
-                # finished selecting # of points required, return
-                if len(selected) > uplimit:
+                # return indices of selected samples, if desired number is selected
+                if len(selected) > max_size:
                     return selected
-                # find all points now within radius of newly selected point
-                elim = kdtree.query_ball_point(x[idx], self.r, eps=self.eps, p=self.p, workers=-1)
-                # turn 'on' bits in bv to make for loop skip indices of eliminated points
-                #   eliminate points from selection
-                for index in elim:
+                # find index of all samples within radius of sample idx (this includes the sample index itself)
+                index_exclude = kdtree.query_ball_point(X[idx], self.r, eps=self.eps, p=self.p, workers=-1)
+                # exclude samples within radius r of sample idx (measure by Minkowski p-norm) from
+                # future consideration by setting their bitarray value to 1
+                for index in index_exclude:
                     bv[index] = 1
 
         return selected
 
-    def select_from_cluster(self, x, num_selected, cluster_ids=None):
-        """
-        Algorithm that uses sphere_exclusion for selecting points from cluster.
+    def select_from_cluster(self, X, size, cluster_ids=None):
+        """Return selected samples from a cluster based on directed sphere exclusion algorithm
 
         Parameters
         ----------
-        x: np.ndarray
-            Feature points.
-        num_selected: int
-            Number of points that need to be selected.
+        X: ndarray of shape (n_samples, n_features)
+           Feature matrix of `n_samples` samples in `n_features` dimensional space.
+        size: int
+            Number of samples to be selected.
         cluster_ids: np.ndarray
-            Indices of points that form a cluster
+            Indices of samples that form a cluster.
 
         Returns
         -------
         selected: list
-            List of ids of selected molecules
+            List of indices of selected samples.
         """
-        if x.shape[0] < num_selected:
+        if X.shape[0] < size:
             raise RuntimeError(
-                f"The number of selected points {num_selected} is greater than the number of points"
-                f"provided {x.shape[0]}."
+                f"Number of samples is less than the requested sample size: {X.shape[0]} < {size}."
             )
-        return optimize_radius(self, x, num_selected, cluster_ids)
+        return optimize_radius(self, X, size, cluster_ids)
 
 
 class GridPartitioning(SelectionBase):
diff --git a/DiverseSelector/methods/tests/test_partition.py b/DiverseSelector/methods/tests/test_partition.py
index c27df3f5..14513472 100644
--- a/DiverseSelector/methods/tests/test_partition.py
+++ b/DiverseSelector/methods/tests/test_partition.py
@@ -39,7 +39,7 @@ def test_directed_sphere_same_number_of_pts():
     """Test DirectSphereExclusion with `size` = number of points in dataset."""
     # (0,0) as the reference point
     x = np.array([[0,0],[0,1],[0,2],[0,3]])
-    selector = DirectedSphereExclusion(r_0=1, tolerance=0)
+    selector = DirectedSphereExclusion(r0=1, tol=0)
     selected = selector.select(arr=x, size=3)
     expected = [1,2,3]
     assert_equal(selected, expected)
@@ -50,7 +50,7 @@ def test_directed_sphere_exclusion_select_more_number_of_pts():
     """Test DirectSphereExclusion on points on the line with `size` < number of points in dataset."""
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5], [0, 6]])
-    selector = DirectedSphereExclusion(r_0=0.5, tolerance=0)
+    selector = DirectedSphereExclusion(r0=0.5, tol=0)
     selected = selector.select(arr=x, size=3)
     expected = [1, 3, 5]
     assert_equal(selected, expected)
@@ -62,7 +62,7 @@ def test_directed_sphere_exclusion_on_line_with_():
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
                   [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5], [0, 6]])
-    selector = DirectedSphereExclusion(r_0=0.5, tolerance=0)
+    selector = DirectedSphereExclusion(r0=0.5, tol=0)
     selected = selector.select(arr=x, size=3)
     expected = [1, 5, 9]
     assert_equal(selected, expected)
@@ -74,7 +74,7 @@ def test_directed_sphere_on_line_with_larger_radius():
     # (0,0) as the reference point
     x = np.array([[0, 0], [0, 1], [0, 1.1], [0, 1.2], [0, 2],
                   [0, 3], [0, 3.1], [0, 3.2], [0, 4], [0, 5]])
-    selector = DirectedSphereExclusion(r_0=2.0, tolerance=0)
+    selector = DirectedSphereExclusion(r0=2.0, tol=0)
     selected = selector.select(arr=x, size=3)
     expected = [1, 5, 9]
     assert_equal(selected, expected)
diff --git a/DiverseSelector/methods/utils.py b/DiverseSelector/methods/utils.py
index d1579e49..7c50d0f8 100644
--- a/DiverseSelector/methods/utils.py
+++ b/DiverseSelector/methods/utils.py
@@ -31,87 +31,87 @@
 ]
 
 
-def optimize_radius(obj, x, num_selected, cluster_ids=None):
+def optimize_radius(obj, X, size, cluster_ids=None):
     """Algorithm that uses sphere exclusion for selecting points from cluster.
 
     Iteratively searches for the optimal radius to obtain the correct number
-    of selected points. If the radius cannot converge to return `num_selected` points,
-    the function returns the closest number of points to `num_selected` as possible.
+    of selected samples. If the radius cannot converge to return `size` points,
+    the function returns the closest number of samples to `size` as possible.
 
     Parameters
     ----------
     obj: object
-        Instance of 'DirectedSphereExclusion' or 'OptiSim' selection class.
-    x: np.ndarray
-        Feature array.
-    num_selected: int
-        Number of points that need to be selected.
+        Instance of `DirectedSphereExclusion` or `OptiSim` selection class.
+    X: ndarray of shape (n_samples, n_features)
+        Feature matrix of `n_samples` samples in `n_features` dimensional space.
+    size: int
+        Number of sample points to select (i.e. size of the subset).
     cluster_ids: np.ndarray
         Indices of points that form a cluster.
 
     Returns
     -------
     selected: list
-        List of ids of selected points.
+        List of indices of selected samples.
     """
-    if x.shape[0] < num_selected:
+    if X.shape[0] < size:
         raise RuntimeError(
-            f"The number of selected points {num_selected} is greater than the number of points"
-            f"provided {x.shape[0]}."
+            f"Size of samples to be selected is greater than existing the number of samples; "
+            f"{size} > {X.shape[0]}."
         )
     # set the limits on # of selected points according to the tolerance percentage
-    error = num_selected * obj.tolerance
-    lowlimit = round(num_selected - error)
-    uplimit = round(num_selected + error)
+    error = size * obj.tol
+    lower_size = round(size - error)
+    upper_size = round(size + error)
 
-    # use initial radius if passed in
+    # select `upper_size` number of samples
     if obj.r is not None:
-        # run the selection algorithm
-        result = obj.algorithm(x, uplimit)
-    # calculate radius from largest feature values
+        # use initial sphere radius
+        selected = obj.algorithm(X, upper_size)
     else:
-        rg = max(np.ptp(x, axis=0)) / num_selected * 3
-        obj.r = rg
-        result = obj.algorithm(x, uplimit)
+        # calculate a sphere radius based on maximum of n_features range
+        # np.ptp returns range of values (maximum - minimum) along an axis
+        obj.r = max(np.ptp(X, axis=0)) / size * 3
+        selected = obj.algorithm(X, upper_size)
 
-    # correct number of points chosen, return selected
-    if len(result) == num_selected:
-        return result
+    # return selected if the correct number of samples chosen
+    if len(selected) == size:
+        return selected
 
-    # if incorrect number of points chosen, then optimize radius
-    if len(result) > num_selected:
-        # Too many points selected, make radius bigger.
+    # optimize radius to select the correct number of samples
+    # first, set a sensible range for optimizing r value within that range
+    if len(selected) > size:
+        # radius should become bigger, b/c too many samples were selected
         bounds = [obj.r, np.inf]
     else:
-        # Too few points selected, make radius smaller
+        # radius should become smaller, b/c too few samples were selected
         bounds = [0, obj.r]
-    n_iter = 0
-    print(f"Number of results {len(result)}, {lowlimit}, {uplimit}")
 
-    while (len(result) < lowlimit or len(result) > uplimit) and n_iter < obj.n_iter:
-        # too many points selected, make radius larger
+    n_iter = 0
+    while (len(selected) < lower_size or len(selected) > upper_size) and n_iter < obj.n_iter:
+        # change sphere radius based on the defined bound
         if bounds[1] == np.inf:
-            rg = bounds[0] * 2
-        # too few points selected, make radius smaller
+            # make sphere radius larger by a factor of 2
+            obj.r = bounds[0] * 2
         else:
-            rg = (bounds[0] + bounds[1]) / 2
-        obj.r = rg
-        # re-run selection with new radius
-        result = obj.algorithm(x, uplimit)
-
-        # adjust upper/lower bounds of radius size to fine tune
-        if len(result) > num_selected:
-            # Too many points selected, raise lower bound
-            bounds[0] = rg
+            # make sphere radius smaller by a factor of 1/2
+            obj.r = (bounds[0] + bounds[1]) / 2
+
+        # re-select samples with the new radius
+        selected = obj.algorithm(X, upper_size)
+
+        # adjust lower/upper bounds of radius range
+        if len(selected) > size:
+            bounds[0] = obj.r
         else:
-            # Too few points selected, lower upper bound
-            bounds[1] = rg
+            bounds[1] = obj.r
         n_iter += 1
+
     # cannot find radius that produces desired number of selected points
     if n_iter >= obj.n_iter:
         warnings.warn(
-            f"Optimal radius finder failed to converge, selected {len(result)} points instead "
-            f"of requested {num_selected}."
+            f"Optimal radius finder failed to converge, selected {len(selected)} points instead "
+            f"of requested {size}."
         )
 
-    return result
+    return selected

From f62cea9202df7d16fba3e62df534e26ddce2b6c7 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Fri, 14 Jul 2023 09:36:41 -0400
Subject: [PATCH 7/7] Fix long line

---
 DiverseSelector/methods/partition.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/DiverseSelector/methods/partition.py b/DiverseSelector/methods/partition.py
index 8744ac4f..ded636d8 100644
--- a/DiverseSelector/methods/partition.py
+++ b/DiverseSelector/methods/partition.py
@@ -139,7 +139,9 @@ def algorithm(self, X, max_size):
                 if len(selected) > max_size:
                     return selected
                 # find index of all samples within radius of sample idx (this includes the sample index itself)
-                index_exclude = kdtree.query_ball_point(X[idx], self.r, eps=self.eps, p=self.p, workers=-1)
+                index_exclude = kdtree.query_ball_point(
+                    X[idx], self.r, eps=self.eps, p=self.p, workers=-1
+                )
                 # exclude samples within radius r of sample idx (measure by Minkowski p-norm) from
                 # future consideration by setting their bitarray value to 1
                 for index in index_exclude: