Repository: AureumChaos/LEAP

Issue:
Add stochastic universal sampling (SUS)
Implement as:

```python3
@curry
@listiter_op
def sus_selection(population: List, n=None, shuffle=True, key=lambda x: x.fitness) -> Iterator
```

This will take a `population` of individuals and randomly find `n` evenly spaced points from which to sample the `population`.

Once all of the `n` points have been sampled you can choose to `shuffle` and resample or just resample by cycling through the existing points, allowing for an endless selection method.

A different selection criteria can be used by changing `key`.

![Statistically_Uniform](https://user-images.githubusercontent.com/33433321/122237720-8ca0ef80-ce8d-11eb-8a88-9cbd2c9cc4bc.png)

This is an example with `n=4` where individuals A, B, C, and F are selected. If `shuffle=True` here then a new 4 points would be generated. If `shuffle=False` here then selection would cycle: `A -> B -> C -> F -> A -> B ...`

Discussion:

Patch:
diff --git a/leap_ec/ops.py b/leap_ec/ops.py
index e58a4d4..e4f4767 100644
--- a/leap_ec/ops.py
+++ b/leap_ec/ops.py
@@ -543,29 +543,122 @@ def proportional_selection(population: List, offset=0, exponent: int = 1,
         >>> pop = Individual.evaluate_population(pop)
         >>> selected = proportional_selection(pop)
     """
-    population_total = 0.0
-    values = np.zeros(len(population), dtype=float)
-
-    # find minimum value to use as offset for negative values
-    if offset == 'pop-min':
-        offset = -min([key(ind) for ind in population])
-
-    # compute values and population total
-    for idx, ind in enumerate(population):
-        ind_val = (key(ind) + offset)**exponent
-        if ind_val < 0:
-            raise ValueError(('negative value found after applying offset. '
-                              f' Offending Individual: {ind}'))
-        population_total += ind_val
-        values[idx] = ind_val
+    # scale and shift to account for possible negative values
+    values = compute_population_values(population, offset=offset,
+                                       exponent=exponent, key=key)
+    assert(len(values) == len(population))
+    values = np.array(values)
+
+    # throw error on negative values since the algorithm does not
+    # work otherwise
+    if (values < 0.0).any():
+        raise ValueError('negative value found after applying offset.')
+
+    population_total = np.sum(values)
     proportions = values / population_total
 
-    # select individuals
     while True:
         choices = random.choices(population, weights=proportions)
         yield choices[0]
 
 
+##############################
+# Function sus_selection
+##############################
+@curry
+@listiter_op
+def sus_selection(population: List, n=None, shuffle: bool = True,
+                  offset=0, exponent: int = 1,
+                  key=lambda x: x.fitness) -> Iterator:
+    """ Returns an individual from a population in proportion to their
+        fitness or another given metric using the stochastic universal
+        sampling algorithm.
+
+        To deal with negative fitness values use `offset='pop-min'` or set a
+        custom offset. A `ValueError` is thrown if the result of adding
+        `offset` to a fitness value results in a negative number. The value
+        of an individual is calculated as follows
+
+        `value = (fitness + offset)^exponent`
+
+        :param population: the population to select from.
+            Should be a list, not an iterator.
+        :param n: the number of evenly spaced points to use in the algorithm.
+            Default is None which uses `len(population)`.
+        :param bool shuffle: if True, `n` points are resampled after one full
+            pass over them. If False, selection repeats over the same `n`
+            points. Defaults to True.
+        :param offset: the offset from zero. If negative fitness values are
+            possible and the minimum is unknown use `offset='pop-min'` for
+            an adaptive offset. Defaults to 0.
+        :param int exponent: the power to which fitness values are raised to.
+            This can be tuned to increase or decrease selection pressure by
+            creating larger or smaller differences between fitness values in
+            the population. Defaults to 1.
+        :param key: a function that computes the metric used to compare
+            individuals. Defaults to fitness.
+        :return: a random individual based on the proportion of the given
+            metric in the population.
+
+        >>> from leap_ec import Individual
+        >>> from leap_ec.binary_rep.problems import MaxOnes
+        >>> from leap_ec.ops import sus_selection
+
+        >>> pop = [Individual([0, 0, 0], problem=MaxOnes()),
+        ...        Individual([0, 0, 1], problem=MaxOnes())]
+        >>> pop = Individual.evaluate_population(pop)
+        >>> selected = sus_selection(pop)
+    """
+    # determine number of points to sample if not specified
+    if n is None:
+        n = len(population)
+
+    # check for non-positive number of points
+    if n <= 0:
+        raise ValueError(f'cannot sample {n} number of points')
+
+    # scale and shift to account for possible negative values
+    values = compute_population_values(population, offset=offset,
+                                       exponent=exponent, key=key)
+    assert(len(values) == len(population))
+    values = np.array(values)
+
+    # throw error on negative values since the algorithm does not
+    # work otherwise
+    if (values < 0.0).any():
+        raise ValueError('negative value found after applying offset.')
+
+    population_total = np.sum(values)
+    even_spacing = population_total / n
+    random_start = np.random.uniform(low=0.0, high=even_spacing)
+    selection_points = [random_start + i*even_spacing for i in range(0, n)]
+    selection_idx = 0
+    population_idx = 0
+    running_sum = 0.0
+    while True:
+        # check if all points have been selected
+        if selection_idx == len(selection_points):
+            # reset to allow for continuous selection
+            if shuffle:
+                random_start = np.random.uniform(low=0.0, high=even_spacing)
+                selection_points = [random_start + i*even_spacing
+                                    for i in range(0, n)]
+            selection_idx = 0
+            running_sum = 0.0
+            population_idx = 0
+
+        current_point = selection_points[selection_idx]
+        # continue until the running sum is greater than the point
+        while running_sum < current_point:
+            running_sum += values[population_idx]
+            population_idx += 1
+        selection_idx += 1
+
+        # yield the individual that caused the running_sum
+        # to move past the current_point
+        yield population[population_idx-1]
+
+
 ##############################
 # Function truncation_selection
 ##############################
@@ -1089,3 +1182,29 @@ def compute_expected_probability(expected_num_mutations: float,
     :return: the corresponding probability of mutation
     """
     return 1.0 / len(individual_genome) * expected_num_mutations
+
+
+##############################
+# function compute_population_values
+##############################
+def compute_population_values(population: List, offset=0, exponent: int = 1,
+                     key=lambda x: x.fitness) -> List:
+    """ Returns a list of values where the zero-point of the population is
+        shifted and the values are scaled by exponentiation.
+
+        :param population: the population to compute values from.
+        :param offset: the offset from zero. Specifying `offset='pop-min'`
+            will use the population's minimum value as the new zero-point.
+            Defaults to 0.
+        :param int exponent: the power to which values are raised to.
+            Defaults to 1.
+        :param key: a function that computes a metric based
+            on an `Individual`.
+        :return: a list of values that have been shifted by `offset` and
+            scaled by `exponent` corresponding to each individual in the
+            population.
+    """
+    values = [key(ind) for ind in population]
+    if offset == 'pop-min':
+        offset = -min(values)
+    return [(val + offset)**exponent for val in values]


Test patch:
diff --git a/tests/test_selection.py b/tests/test_selection.py
index a658b48..c6bb843 100644
--- a/tests/test_selection.py
+++ b/tests/test_selection.py
@@ -14,11 +14,175 @@ from leap_ec.real_rep.problems import SpheroidProblem
 
 
 ##############################
-# Tests for fitness_proportional_selection()
+# Tests for sus_selection()
+##############################
+def test_sus_selection1():
+    ''' Test of a deterministic case of stochastic universal sampling '''
+    # Make a population where sus_selection has an obvious
+    # reproducible choice
+    pop = [Individual([0, 0, 0], problem=MaxOnes()),
+           Individual([1, 1, 1], problem=MaxOnes())]
+
+    pop = Individual.evaluate_population(pop)
+    # This selection operator will always choose the [1, 1, 1] individual
+    # since [0, 0, 0] has zero fitness
+    selector = ops.sus_selection(pop)
+
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+    # run one more time to test shuffle
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+
+@pytest.mark.stochastic
+def test_sus_selection_shuffle():
+    ''' Test of a stochastic case of SUS selection '''
+    # Make a population where sus_selection has an obvious
+    # reproducible choice
+    # Proportions here should be 1/4 and 3/4, respectively
+    pop = [Individual([0, 1, 0], problem=MaxOnes()),
+           Individual([1, 1, 1], problem=MaxOnes())]
+
+    # Assign a unique identifier to each individual
+    pop[0].id = 0
+    pop[1].id = 1
+
+    # We first need to evaluate all the individuals so that
+    # selection has fitnesses to compare
+    pop = Individual.evaluate_population(pop)
+    selected = ops.sus_selection(pop)
+
+    N = 1000
+    p_thresh = 0.1
+    observed_dist = statistical_helpers.collect_distribution(
+        lambda: next(selected).id, samples=N)
+    expected_dist = {pop[0].id: 0.25*N, pop[1].id: 0.75*N}
+    print(f"Observed: {observed_dist}")
+    print(f"Expected: {expected_dist}")
+    assert(statistical_helpers.stochastic_equals(expected_dist,
+                                                 observed_dist, p=p_thresh))
+
+
+def test_sus_selection_offset():
+    ''' Test of SUS selection with a non-default offset '''
+    pop = [Individual([0, 0, 0], problem=MaxOnes()),
+           Individual([1, 1, 1], problem=MaxOnes())]
+
+    # evaluate population and negate fitness of second individual
+    pop = Individual.evaluate_population(pop)
+    pop[1].fitness = -pop[1].fitness
+
+    # now we try to evaluate normally (this should throw a ValueError)
+    # due to the negative fitness
+    with pytest.raises(ValueError):
+        selector = ops.sus_selection(pop)
+        selected = next(selector)
+    # it should work by setting the offset to +3
+    # this adds 3 to each fitness value, making the second
+    # individual's fitness 0.
+    selector = ops.sus_selection(pop, offset=3)
+
+    # we expect the first individual to always be selected
+    # since the new zero point is now -3.
+    selected = next(selector)
+    assert(selected.genome == [0, 0, 0])
+
+    selected = next(selector)
+    assert(selected.genome == [0, 0, 0])
+
+
+def test_sus_selection_pop_min():
+    ''' Test of SUS selection with pop-min offset '''
+    # Create a population of positive fitness individuals
+    # scaling the fitness by the population minimum makes it so the
+    # least fit member never gets selected.
+    pop = [Individual([0, 1, 0], problem=MaxOnes()),
+           Individual([1, 1, 1], problem=MaxOnes())]
+
+    pop = Individual.evaluate_population(pop)
+
+    selector = ops.sus_selection(pop, offset='pop-min')
+
+    # we expect that the second individual is always selected
+    # since the new zero point will be at the minimum fitness
+    # of the population
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+
+def test_sus_selection_custom_key():
+    ''' Test of SUS selection with custom evaluation '''
+    pop = [Individual([0, 0, 0], problem=MaxOnes()),
+           Individual([1, 1, 1], problem=MaxOnes())]
+
+    def custom_key(individual):
+        ''' Returns fitness based on MaxZeros '''
+        return individual.genome.count(0)
+
+    pop = Individual.evaluate_population(pop)
+    selector = ops.sus_selection(pop, key=custom_key)
+
+    # we expect the first individual to always be selected
+    # since its genome is all 0s
+    selected = next(selector)
+    assert(selected.genome == [0, 0, 0])
+
+    selected = next(selector)
+    assert(selected.genome == [0, 0, 0])
+
+
+def test_sus_selection_num_points():
+    ''' Test of SUS selection with varying `n` random points '''
+    # the second individual should always be selected
+    pop = [Individual([0, 0, 0], problem=MaxOnes()),
+           Individual([1, 1, 1], problem=MaxOnes())]
+
+    pop = Individual.evaluate_population(pop)
+    # with negative points
+    with pytest.raises(ValueError):
+        selector = ops.sus_selection(pop, n=-1)
+        selected = next(selector)
+
+    # with n = None (default)
+    selector = ops.sus_selection(pop, n=None)
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+    # with n less than len(population)
+    selector = ops.sus_selection(pop, n=1)
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+    # with n greater than len(population)
+    selector = ops.sus_selection(pop, n=3)
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+    selected = next(selector)
+    assert(selected.genome == [1, 1, 1])
+
+
+##############################
+# Tests for proportional_selection()
 ##############################
 def test_proportional_selection1():
     ''' Test of a deterministic case of proportional selection '''
-    # Make a population where fitness_proportional_selection has an obvious
+    # Make a population where proportional_selection has an obvious
     # reproducible choice
     pop = [Individual([0, 0, 0], problem=MaxOnes()),
            Individual([1, 1, 1], problem=MaxOnes())]