new patch correcting an error in the calculation of full_gradients and full_transition_matrix

Author: Socrats

CHANGELOG.md

@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/) and this project adheres
 to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
 
+## [0.1.11.patch2] - 16-11-2022
+
+### Fixed
+
+- Fixed error in `calculate_full_transition_matrix`. The previous implementation tried to calculate the fitness for
+  strategies with 0 counts in the population. This would lead to an error when calculating the probability of group
+  forming. Now, we calculate directly the probability of a strategy with 0 counts increasing by 1 individual in the
+  population, i.e., prob = (n_decreasing_strategy / population_size) * (mu / (1 - nb_strategies)).
+
 ## [0.1.11.patch1] - 10-11-2022
 
 ### Fixed

cpp/src/version.h

@@ -18,6 +18,6 @@
 #ifndef EGTTOOLS_VERSION_H
 #define EGTTOOLS_VERSION_H
 
-#define EGTTOOLS_VERSION 0.1.11.1
+#define EGTTOOLS_VERSION 0.1.11.2
 
 #endif//EGTTOOLS_VERSION_H

pybind11

@@ -461,10 +461,17 @@ def full_gradient_selection(self, population_state: np.ndarray, beta: float) ->
         probability_selecting_strategy_first = population_state / self.pop_size
         probability_selecting_strategy_second = population_state / self.pop_size
         probabilities = np.outer(probability_selecting_strategy_first, probability_selecting_strategy_second)
-        fitness = np.asarray([[self.full_fitness(i, j, population_state) for i in
-                               range(len(population_state))] for j in range(len(population_state))])
-        fitness[np.isnan(fitness)] = 0
-        return (probabilities * np.tanh((beta / 2) * fitness)).sum(axis=0) * (1 - self.mu) + self.mu
+        fitness = np.zeros(shape=(self.nb_strategies, self.nb_strategies))
+        for j in range(self.nb_strategies):
+            if population_state[j] == 0:
+                continue
+            for i in range(self.nb_strategies):
+                if population_state[i] == 0:
+                    continue
+                fitness[j, i] = self.full_fitness(i, j, population_state)
+
+        return (probabilities * np.tanh((beta / 2) * fitness)).sum(axis=0) * (1 - self.mu) + (
+                    self.mu / (self.nb_strategies - 1)) * probability_selecting_strategy_second
 
     def full_gradient_selection_without_mutation(self, population_state: np.ndarray, beta: float) -> np.ndarray:
         """
@@ -488,9 +495,14 @@ def full_gradient_selection_without_mutation(self, population_state: np.ndarray,
         probability_selecting_strategy_first = population_state / self.pop_size
         probability_selecting_strategy_second = population_state / self.pop_size
         probabilities = np.outer(probability_selecting_strategy_first, probability_selecting_strategy_second)
-        fitness = np.asarray([[self.full_fitness(i, j, population_state) for i in
-                               range(len(population_state))] for j in range(len(population_state))])
-        fitness[np.isnan(fitness)] = 0
+        fitness = np.zeros(shape=(self.nb_strategies, self.nb_strategies))
+        for j in range(self.nb_strategies):
+            if population_state[j] == 0:
+                continue
+            for i in range(self.nb_strategies):
+                if population_state[i] == 0:
+                    continue
+                fitness[j, i] = self.full_fitness(i, j, population_state)
 
         return (probabilities * np.tanh((beta / 2) * fitness)).sum(axis=0)
 
@@ -588,17 +600,21 @@ def calculate_full_transition_matrix(self, beta: float, *args: Optional[list]) -
                     increase = np.where(np.array(permutation) == 1)[0][0]
                     decrease = np.where(np.array(permutation) == -1)[0][0]
 
-                    # now we calculate the transition probability
-                    fitness_diff = self.full_fitness(decrease, increase, current_state)
-                    # Probability that the individual that will die is selected and that the individual that
-                    # will be imitated is selected times the probability of imitation
-                    prob = not_mu * (current_state[increase] / (self.pop_size - 1))
-                    prob *= StochDynamics.fermi(beta, fitness_diff)
-                    # The probability that there will not be a mutation event times the probability of the transition
-                    # plus the probability that if there is a mutation event, the dying strategy is selected
-                    # times the probability that it mutates into the increasing strategy
-                    prob = (current_state[decrease] / self.pop_size) * (prob + mutation_probability)
-                    total_prob += prob
+                    if current_state[increase] == 0:
+                        prob = (current_state[decrease] / self.pop_size) * mutation_probability
+                    else:
+                        # now we calculate the transition probability
+                        fitness_diff = self.full_fitness(decrease, increase, current_state)
+                        # Probability that the individual that will die is selected and that the individual that
+                        # will be imitated is selected times the probability of imitation
+                        prob = not_mu * (current_state[increase] / (self.pop_size - 1))
+                        prob *= StochDynamics.fermi(beta, fitness_diff)
+                        # The probability that there will not be a mutation event times the probability
+                        # of the transition
+                        # plus the probability that if there is a mutation event, the dying strategy is selected
+                        # times the probability that it mutates into the increasing strategy
+                        prob = (current_state[decrease] / self.pop_size) * (prob + mutation_probability)
+                        total_prob += prob
 
                     transitions[i, new_state_index] = prob