Repository: Happy-Algorithms-League/hal-cgp
Issue category: Performance Issue

Issue:
Issue warning when `mutation_rate` is high and `mutate` potentially takes long
#157 made sure that exactly the correct number of mutations occurs in a genome by checking that the mutated genome differs from the original one in `n_mutations` position. when using a high mutation rate and small number of primitives, this can lead to a long runtime for `mutate` as it struggles to find positions in the genome to mutate and alternative genes to insert. we should figure out exactly under which conditions this tends to occur and issue a user warning when we detect such a situation.

@HenrikMettler, since you're already familiar with that part of the code, maybe you want to look into this?


Hints:
btw, an alternative solution to this problem would be to have a maximum number iterations in the while loop (https://github.com/Happy-Algorithms-League/hal-cgp/blob/master/cgp/genome.py#L377) and return even if the desired number of differences was not reached while issuing a warning to the user
What about keeping track of the genes that were previously mutated and only allowing `gene_idx` to take different values? 
Perhaps I am missing something here, but won't this problem go away once we finally fix the current implementation to not apply a fixed number of mutations but rather mutate each gene with a fixed probability, as the variable names suggest, i.e. once #124 is closed?
Actually #124 (sorry I will do a PR soonish) only calculates a different `n_mutations` for every individual in every generation (from a binomial distribution). It does not apply a mutation with a probability to every single gene (this would be computationally more expensive I think). So the `while` loop in `Genome.mutate` still looks the same, just with different values for `n_mutations` 
I doubt that making a coin flip for each gene in each individual is slower. 
For instance, running this code:

```python
import numpy as np
import time

n_individuals = 10
n_genes = 15

p = 0.2

t0 = time.time()
for i in range(1000):
    np.random.binomial(n_genes, p, size=n_individuals)

print(f"Binomial took {time.time() - t0} seconds.")

t0 = time.time()
for i in range(1000):
    indx = np.random.rand(n_individuals, n_genes) < p

print(f"Bernoulli took {time.time() - t0} seconds.")
```

gives me:
```
Binomial took 0.0036504268646240234 seconds.
Bernoulli took 0.002948284149169922 seconds.
```

yes, it seems like you are right @mschmidt87, my intuition would've been with @HenrikMettler, oh well this is why one should always benchmark :flushed:

however, i don't think checking every gene individually would remove the fundamental problem that we would need to detect if a user wants to modify more genes than could actually be changed. two things can happen otherwise: i) we have an infinite loop ii) the effective mutation rate is smaller than the one requested by the user without them noticing.

one suggestion to resolve this (just discussed with @HenrikMettler) would be to determine the set of all "mutable genes", make sure that set is larger than `n_mutations`, and then pick `n_mutations` genes from that set to actually be mutated. a PR will follow
If I understand correctly, there are several types of genes that cannot be mutated:
- input region: no gene can be mutated
- output region: function gene and is_active_input gene cannot be mutated

Thus, the number of mutable genes in a genome is `len(genome) - n_immutable_genes`. This leads to a problem if the user understand `mutation_rate` as determining the number of mutated genes by doing `n_mutations = int(mutation_rate * len(self.dna))`. This is what we're doing now and how it is documented. If we do it like that, then we need to add additional checks etc.

However, if `mutation_rate` is understood as the probability for each mutable gene to be mutated, then these problems just go away. We can simply do:
```python
## is_mutable = array of length len(dna) with True/False for each gene
mutated_idx = np.where(np.logical_and(np.random.rand(len(dna)) < p, is_mutable))
dna[mutated_idx] = rng.choice(permissible_inputs, size=len(mutated_idx))
```

This is just a question of the exact definition/understanding of the mutation rate, and I would opt for the interpretation that is simpler and makes more sense, in my opinion.
I agree with your analysis and (of course short ;)) implementation. i would like to add that for some parameter configurations there might be only one option for a gene, which makes it questionable whether replacing this gene with itself should be counted as a mutation or not.

that being said, I think the two approaches (calculating the total number of mutations once as `n_mutations = int(mutation_rate * len(self.dna))` or mutating each gene with a certain probability) are identical as long as we count trying to mutate an gene that can not be mutated (as in your definition above) as a successful mutation. in both cases, the number of genes that are different between parent and offspring will be smaller than expected based on the total number of genes and the mutation rate.

after this discussion and considering the implementation i'm starting to like this version a bit better, however the issue remains that adopting this version offsprings can be generated _which are genetically identical to their parents_. as we rely on mutations to drive search this null operation seems wasteful to me, but from quickly scanning other implementations it seems quite common. this might be related to, as you said, this interpretation of the mutation rate being easier to explain/understand.

@HenrikMettler do you have any strong feelings either way? i mean besides having already working code for one ;)
I don't, given that computationally there seems to be no significant difference. (I vaguely remember reading somewhere in the literature that there is one, so I will quickly check again - but probably it is just my memory fooling me)

I would just argue that for the `permissible_inputs` we do: 
`permissible_inputs = self._permissible_inputs()`
`permissible_inputs.remove(dna[idx])`

Since a gene for which this list is empty (eg. a hidden gene, that can only have one possible input, in a specific graph structure) it should be that: `is_mutable = false`  I would do this before defining which elements should be mutated

Patch:
diff --git a/cgp/genome.py b/cgp/genome.py
index 3c18a87f..4e8f035d 100644
--- a/cgp/genome.py
+++ b/cgp/genome.py
@@ -107,7 +107,7 @@ def dna(self) -> List[int]:
     def dna(self, value: List[int]) -> None:
         self._validate_dna(value)
         self._dna = value
-        self._initialize_unkown_parameters()
+        self._initialize_unknown_parameters()
 
     @property
     def _n_hidden(self) -> int:
@@ -332,6 +332,9 @@ def _is_output_region(self, region_idx: int) -> bool:
     def _is_function_gene(self, gene_idx: int) -> bool:
         return (gene_idx % self._length_per_region) == 0
 
+    def _is_hidden_input_gene(self, gene_idx: int, region_idx: int) -> bool:
+        return self._is_hidden_region(region_idx) & ((gene_idx % self._length_per_region) != 0)
+
     def _is_active_input_gene(self, gene_idx: int) -> bool:
         input_index = gene_idx % self._length_per_region
         assert input_index > 0
@@ -346,87 +349,113 @@ def _is_active_input_gene(self, gene_idx: int) -> bool:
         else:
             assert False  # should never be reached
 
+    def _select_gene_indices_for_mutation(
+        self, mutation_rate: float, rng: np.random.RandomState
+    ) -> List[int]:
+        """Selects the gene indices for mutations
+
+        Parameters
+        ----------
+        mutation_rate : float
+            Probability of a gene to be mutated, between 0 (excluded) and 1 (included).
+        rng : numpy.random.RandomState
+            Random number generator instance to use for selecting the indices.
+
+        Returns
+        ----------
+        selected_gene_indices: np.ndarray
+            indices of the genes selected for mutation.
+        """
+
+        selected_gene_indices = np.nonzero(rng.rand(len(self.dna)) < mutation_rate)[0]
+        return selected_gene_indices
+
     def mutate(self, mutation_rate: float, rng: np.random.RandomState):
         """Mutate the genome.
 
         Parameters
         ----------
         mutation_rate : float
-            Proportion of genes to be mutated, between 0 and 1.
+            Probability of a gene to be mutated, between 0 (excluded) and 1 (included).
         rng : numpy.random.RandomState
             Random number generator instance to use for mutating.
 
+
         Returns
         ----------
         bool
             True if only inactive regions of the genome were mutated, False otherwise.
         """
 
-        def count_dna_differences(dna0: List[int], dna1: List[int]) -> int:
-            return len([1 for gene0, gene1 in zip(dna0, dna1) if gene0 != gene1])
-
-        n_mutations = int(mutation_rate * len(self.dna))
-        assert n_mutations > 0
-
         graph = CartesianGraph(self)
         active_regions = graph.determine_active_regions()
-
         dna = list(self._dna)
-
         only_silent_mutations = True
-        while count_dna_differences(self.dna, dna) < n_mutations:
 
-            gene_idx = rng.randint(0, self._n_genes)
-            region_idx = gene_idx // self._length_per_region
+        selected_gene_indices = self._select_gene_indices_for_mutation(mutation_rate, rng)
 
-            if self._is_input_region(region_idx):
-                continue  # nothing to do here
+        for (gene_idx, allele) in zip(selected_gene_indices, np.array(dna)[selected_gene_indices]):
 
-            elif self._is_output_region(region_idx):
-                silent = self._mutate_output_region(dna, gene_idx, rng)
+            region_idx = gene_idx // self._length_per_region
 
-            elif self._is_hidden_region(region_idx):
-                silent = self._mutate_hidden_region(dna, gene_idx, active_regions, rng)
+            permissible_values = self._determine_alternative_permissible_values(
+                gene_idx, allele, region_idx
+            )
+            if len(permissible_values) > 0:
 
-            else:
-                assert False  # should never be reached
+                dna[gene_idx] = rng.choice(permissible_values)
+                silent = region_idx not in active_regions
 
-            only_silent_mutations = only_silent_mutations and silent
+                only_silent_mutations = only_silent_mutations and silent
 
         self.dna = dna
         return only_silent_mutations
 
-    def _mutate_output_region(
-        self, dna: List[int], gene_idx: int, rng: np.random.RandomState
-    ) -> bool:
-        assert self._is_gene_in_output_region(gene_idx)
+    def _determine_alternative_permissible_values(
+        self, gene_idx: int, allele: int, region_idx: int
+    ) -> List[int]:
 
-        if not self._is_function_gene(gene_idx) and self._is_active_input_gene(gene_idx):
-            permissible_inputs = self._permissible_inputs_for_output_region()
-            dna[gene_idx] = rng.choice(permissible_inputs)
-            return False
-        else:
-            return True
+        if self._is_input_region(region_idx):
+            return []  # genes in input regions have no alternative permissible values
 
-    def _mutate_hidden_region(
-        self, dna: List[int], gene_idx: int, active_regions: List[int], rng: np.random.RandomState
-    ) -> bool:
+        elif self._is_hidden_region(region_idx):
+            return self._determine_alternative_permissible_values_hidden_gene(
+                gene_idx, allele, region_idx
+            )
 
-        assert self._is_gene_in_hidden_region(gene_idx)
+        elif self._is_output_region(region_idx):
+            return self._determine_alternative_permissible_values_output_gene(gene_idx, allele)
 
-        region_idx = gene_idx // self._length_per_region
-        silent_mutation = region_idx not in active_regions
+        else:
+            assert False  # should never be reached
 
+    def _determine_alternative_permissible_values_hidden_gene(
+        self, gene_idx: int, allele: int, region_idx: int
+    ) -> List[int]:
         if self._is_function_gene(gene_idx):
-            dna[gene_idx] = self._primitives.sample_allele(rng)
-            return silent_mutation
+            permissible_values = list(np.arange(len(self._primitives._primitives)))
 
+        elif self._is_hidden_input_gene(gene_idx, region_idx):
+            permissible_values = self._permissible_inputs(region_idx)
+
+        else:
+            assert False
+        permissible_values.remove(allele)
+        return permissible_values
+
+    def _determine_alternative_permissible_values_output_gene(
+        self, gene_idx: int, allele: int
+    ) -> List[int]:
+        input_index = (
+            gene_idx % self._length_per_region
+        )  # assumes that the second gene in all output regions is the index of the input
+        if input_index == 1:
+            permissible_values = self._permissible_inputs_for_output_region()
+            permissible_values.remove(allele)
         else:
-            permissible_inputs = self._permissible_inputs(region_idx)
-            dna[gene_idx] = rng.choice(permissible_inputs)
+            permissible_values = []
 
-            silent_mutation = silent_mutation or (not self._is_active_input_gene(gene_idx))
-            return silent_mutation
+        return permissible_values
 
     @property
     def primitives(self) -> Primitives:
@@ -484,7 +513,7 @@ def update_parameters_from_torch_class(self, torch_cls: "torch.nn.Module") -> bo
 
         return any_parameter_updated
 
-    def _initialize_unkown_parameters(self) -> None:
+    def _initialize_unknown_parameters(self) -> None:
         for region_idx, region in self.iter_hidden_regions():
             node_id = region[0]
             node_type = self._primitives[node_id]
diff --git a/cgp/population.py b/cgp/population.py
index f2ccce23..114d70f9 100755
--- a/cgp/population.py
+++ b/cgp/population.py
@@ -25,8 +25,7 @@ def __init__(
         n_parents : int
             Number of parent individuals.
         mutation_rate : float
-            Rate of mutations determining the number of genes to be
-            mutated for offspring creation, between 0 and 1.
+            Probability of a gene to be mutated, between 0 (excluded) and 1 (included).
         seed : int
             Seed for internal random number generator.
         genome_params : dict
@@ -34,8 +33,8 @@ def __init__(
         """
         self.n_parents = n_parents  # number of individuals in parent population
 
-        if not (0.0 < mutation_rate and mutation_rate < 1.0):
-            raise ValueError("mutation rate needs to be in (0, 1)")
+        if not (0.0 < mutation_rate and mutation_rate <= 1.0):
+            raise ValueError("mutation rate needs to be in (0, 1]")
         self._mutation_rate = mutation_rate  # probability of mutation per gene
 
         self.seed = seed


Test patch:
diff --git a/test/test_ea_mu_plus_lambda.py b/test/test_ea_mu_plus_lambda.py
index cf0e1b2b..8f0521ea 100644
--- a/test/test_ea_mu_plus_lambda.py
+++ b/test/test_ea_mu_plus_lambda.py
@@ -196,6 +196,8 @@ def objective(individual):
         individual.fitness = float(individual.idx)
         return individual
 
+    population_params["mutation_rate"] = 1.0  # ensures every offspring has mutations
+
     pop = cgp.Population(**population_params, genome_params=genome_params)
     ea = cgp.ea.MuPlusLambda(**ea_params)
 
diff --git a/test/test_genome.py b/test/test_genome.py
index 55d1a892..17106b86 100644
--- a/test/test_genome.py
+++ b/test/test_genome.py
@@ -3,6 +3,7 @@
 
 import cgp
 from cgp.genome import ID_INPUT_NODE, ID_OUTPUT_NODE, ID_NON_CODING_GENE
+from cgp.cartesian_graph import CartesianGraph
 
 
 def test_check_dna_consistency():
@@ -346,53 +347,91 @@ def test_is_gene_in_output_region(rng_seed):
     assert not genome._is_gene_in_output_region(11)
 
 
-def test_mutate_hidden_region(rng_seed):
+def test_mutation_rate(rng_seed, mutation_rate):
+    n_inputs = 1
+    n_outputs = 1
+    n_columns = 4
+    n_rows = 3
+    genome = cgp.Genome(n_inputs, n_outputs, n_columns, n_rows, None, (cgp.Add, cgp.Sub))
     rng = np.random.RandomState(rng_seed)
-    genome = cgp.Genome(1, 1, 3, 1, None, (cgp.Add, cgp.ConstantFloat))
-    dna = [
+    genome.randomize(rng)
+
+    def count_n_immutable_genes(n_inputs, n_outputs, n_rows):
+        length_per_region = genome.primitives.max_arity + 1  # function gene + input gene addresses
+        n_immutable_genes = n_inputs * length_per_region  # none of the input genes are mutable
+        n_immutable_genes += n_outputs * (
+            length_per_region - 1
+        )  # only one gene per output can be mutated
+        if n_inputs == 1:
+            n_immutable_genes += (
+                n_rows * genome.primitives.max_arity
+            )  # input gene addresses in the first (hidden) column can't be mutated
+            # if only one input node exists
+        return n_immutable_genes
+
+    def count_mutations(dna0, dna1):
+        n_differences = 0
+        for (allele0, allele1) in zip(dna0, dna1):
+            if allele0 != allele1:
+                n_differences += 1
+        return n_differences
+
+    n_immutable_genes = count_n_immutable_genes(n_inputs, n_outputs, n_rows)
+    n_mutations_mean_expected = mutation_rate * (len(genome.dna) - n_immutable_genes)
+    n_mutations_std_expected = np.sqrt(
+        (len(genome.dna) - n_immutable_genes) * mutation_rate * (1 - mutation_rate)
+    )
+
+    n = 10000
+    n_mutations = []
+    for _ in range(n):
+        dna_old = genome.dna
+        genome.mutate(mutation_rate, rng)
+        n_mutations.append(count_mutations(dna_old, genome.dna))
+
+    assert np.mean(n_mutations) == pytest.approx(n_mutations_mean_expected, rel=0.04)
+    assert np.std(n_mutations) == pytest.approx(n_mutations_std_expected, rel=0.04)
+
+
+def test_only_silent_mutations(genome_params, mutation_rate, rng_seed):
+    genome = cgp.Genome(**genome_params)
+    rng = np.random.RandomState(rng_seed)
+    genome.randomize(rng)
+
+    only_silent_mutations = genome.mutate(mutation_rate=0, rng=rng)
+    assert only_silent_mutations is True
+
+    only_silent_mutations = genome.mutate(mutation_rate=1, rng=rng)
+    assert not only_silent_mutations
+
+    dna_fixed = [
+        ID_INPUT_NODE,
+        ID_NON_CODING_GENE,
+        ID_NON_CODING_GENE,
         ID_INPUT_NODE,
         ID_NON_CODING_GENE,
         ID_NON_CODING_GENE,
+        2,
         1,
         0,
         0,
-        1,
         0,
         0,
         0,
         0,
-        2,
-        ID_OUTPUT_NODE,
-        3,
-        ID_NON_CODING_GENE,
-    ]
-    genome.dna = list(dna)
-    active_regions = cgp.CartesianGraph(genome).determine_active_regions()
-
-    # mutating any gene in inactive region returns True
-    assert genome._mutate_hidden_region(list(dna), 3, active_regions, rng) is True
-    assert genome._mutate_hidden_region(list(dna), 4, active_regions, rng) is True
-    assert genome._mutate_hidden_region(list(dna), 5, active_regions, rng) is True
-
-    # mutating function gene in active region returns False
-    assert genome._mutate_hidden_region(list(dna), 6, active_regions, rng) is False
-    # mutating inactive genes in active region returns True
-    assert genome._mutate_hidden_region(list(dna), 7, active_regions, rng) is True
-    assert genome._mutate_hidden_region(list(dna), 8, active_regions, rng) is True
-
-    # mutating any gene in active region without silent genes returns False
-    assert genome._mutate_hidden_region(list(dna), 9, active_regions, rng) is False
-    assert genome._mutate_hidden_region(list(dna), 10, active_regions, rng) is False
-    assert genome._mutate_hidden_region(list(dna), 11, active_regions, rng) is False
-
-
-def test_mutate_output_region(rng_seed):
-    rng = np.random.RandomState(rng_seed)
-    genome = cgp.Genome(1, 1, 2, 1, None, (cgp.Add,))
-    dna = [
-        ID_INPUT_NODE,
-        ID_NON_CODING_GENE,
-        ID_NON_CODING_GENE,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
+        0,
         0,
         0,
         0,
@@ -403,32 +442,110 @@ def test_mutate_output_region(rng_seed):
         2,
         ID_NON_CODING_GENE,
     ]
+    genome.dna = dna_fixed
+    graph = CartesianGraph(genome)
+    active_regions = graph.determine_active_regions()
+    length_per_region = genome.primitives.max_arity + 1  # function gene + input gene addresses
+    gene_to_be_mutated_non_active = (
+        3 * length_per_region
+    )  # operator gene of 2nd hidden node, a silent node in this dna configuration
+
+    def select_gene_indices_silent(mutation_rate, rng):
+        selected_gene_indices = [gene_to_be_mutated_non_active]
+        return selected_gene_indices
+
+    genome._select_gene_indices_for_mutation = (
+        select_gene_indices_silent  # monkey patch the selection of indices to select a silent gene
+    )
+    genome.dna = dna_fixed
+    only_silent_mutations = genome.mutate(mutation_rate, rng)
+    assert only_silent_mutations is True
+
+    gene_to_be_mutated_active = (
+        active_regions[-1] * length_per_region
+    )  # operator gene of the 1st active hidden node,
+    # should always be mutable and result in a non-silent-mutation
+
+    def select_gene_indices_non_silent(mutation_rate, rng):
+        selected_gene_indices = [gene_to_be_mutated_active]
+        return selected_gene_indices
+
+    genome._select_gene_indices_for_mutation = select_gene_indices_non_silent
+    # monkey patch the selection of indices to select a non-silent gene
+    only_silent_mutations = genome.mutate(mutation_rate, rng)
+    assert not only_silent_mutations
+
+
+def test_permissible_values_hidden(rng_seed):
+    genome_params = {
+        "n_inputs": 2,
+        "n_outputs": 1,
+        "n_columns": 3,
+        "n_rows": 3,
+        "levels_back": 2,
+        "primitives": (cgp.Add, cgp.Sub, cgp.ConstantFloat),
+    }
+
+    genome = cgp.Genome(**genome_params)
+
+    # test function gene
+    gene_idx = 6  # function gene of first hidden region
+    allele = 0  #
+    region_idx = None  # can be any number, since not touched for function gene
+    permissible_func_gene_values = genome._determine_alternative_permissible_values_hidden_gene(
+        gene_idx, allele, region_idx
+    )
+    assert permissible_func_gene_values == [
+        1,
+        2,
+    ]  # function idx 1,2 (cgp.Sub, cgp.ConstantFloat)
+
+    # test input gene
+    gene_idx = 16  # first input gene in second column of hidden region
+    allele = 0
+    region_idx = gene_idx // (
+        genome.primitives.max_arity + 1
+    )  # function gene + input gene addresses
+    permissible_input_gene_values = genome._determine_alternative_permissible_values_hidden_gene(
+        gene_idx, allele, region_idx
+    )
 
-    assert genome._mutate_output_region(list(dna), 9, rng) is True
-    assert genome._mutate_output_region(list(dna), 10, rng) is False
-    assert genome._mutate_output_region(list(dna), 11, rng) is True
+    assert permissible_input_gene_values == [
+        1,
+        2,
+        3,
+        4,
+    ]  # gene indices of input, 1st hidden layer
 
 
-@pytest.mark.parametrize("mutation_rate", [0.02, 0.05, 0.2])
-def test_correct_number_of_mutations(mutation_rate, rng_seed):
+def test_permissible_values_output(rng_seed):
+    genome_params = {
+        "n_inputs": 2,
+        "n_outputs": 1,
+        "n_columns": 3,
+        "n_rows": 3,
+        "levels_back": 2,
+        "primitives": (cgp.Add, cgp.Sub, cgp.ConstantFloat),
+    }
+    genome = cgp.Genome(**genome_params)
 
-    n_inputs = 2
-    n_outputs = 1
-    n_columns = 10
-    n_rows = 2
-    levels_back = 5
-    primitives = (cgp.Add, cgp.Sub, cgp.Mul, cgp.Div, cgp.ConstantFloat)
+    gene_idx_function = 33
+    gene_idx_input0 = 34
+    gene_idx_input1 = 35
 
-    rng = np.random.RandomState(rng_seed)
-    genome = cgp.Genome(n_inputs, n_outputs, n_columns, n_rows, levels_back, primitives)
-    genome.randomize(rng)
+    allele = 5  # set current value for input0 -> should be excluded from permissible values
 
-    n_mutations = 0
-    genome_new = genome.clone()
-    genome_new.mutate(mutation_rate, rng)
-    for (gene_0, gene_1) in zip(genome.dna, genome_new.dna):
-        if gene_0 != gene_1:
-            n_mutations += 1
+    permissible_values = genome._determine_alternative_permissible_values_output_gene(
+        gene_idx_function, allele
+    )
+    assert permissible_values == []
 
-    n_mutations_expected = int(mutation_rate * len(genome.dna))
-    assert n_mutations == n_mutations_expected
+    permissible_values = genome._determine_alternative_permissible_values_output_gene(
+        gene_idx_input0, allele
+    )
+    assert permissible_values == [0, 1, 2, 3, 4, 6, 7, 8, 9, 10]
+
+    permissible_values = genome._determine_alternative_permissible_values_output_gene(
+        gene_idx_input1, allele
+    )
+    assert permissible_values == []
diff --git a/test/test_hl_api.py b/test/test_hl_api.py
index 9ea76dd1..50837686 100644
--- a/test/test_hl_api.py
+++ b/test/test_hl_api.py
@@ -65,6 +65,7 @@ def test_parallel_population(population_params, genome_params, ea_params):
         fitness_per_n_processes.append(
             _test_population(population_params, genome_params, ea_params)
         )
+
     assert fitness_per_n_processes[0] == pytest.approx(fitness_per_n_processes[1])
     assert fitness_per_n_processes[0] == pytest.approx(fitness_per_n_processes[2])