Added mutation rates

Author: mtomasini

docs/source/overview.md

@@ -18,7 +18,7 @@ Currently, each individual's culture is represented by a set of {math}`N` featur
 
 ### Neutral model
 
-The Neutral model follows loosely the set up selected by other authors working on the spread of cultural traits (_e.g._ *Patterns in space and time: simulating cultural transmission in archaeology*, Marko Porčić). While in most works the Neutral Model assumes some sort of mutation rate in the traits, in the form of errors in the process of copying traits, we assume no copying error and instead use the simplest model where one focal individual at random in each subpopulation and each generation copies one trait at random from another individual in the same sub-population. Because of the absence of mutations, given enough time, this process is expected to lead to uniformity in an isolated subpopulation. 
+The Neutral model follows loosely the set up selected by other authors working on the spread of cultural traits (_e.g._ *Patterns in space and time: simulating cultural transmission in archaeology*, Marko Porčić). While in most works the Neutral Model assumes some sort of mutation rate in the traits, in the form of errors in the process of copying traits, we assume no copying error and instead use the simplest model where one focal individual at random in each subpopulation and each generation copies one trait at random from another individual in the same subpopulation. Because of the absence of mutations, given enough time, this process is expected to lead to uniformity in an isolated subpopulation. 
 
 ### Axelrod model
 

exploration.ipynb

@@ -14,10 +14,11 @@ class Individual():
         original_deme_id (int): Identifier of the deme where the individual originated.
         number_of_features (int): Number of cultural features of the individual.
         number_of_traits (int): Number of traits per feature of the individual.
+        mutation_rate (float): Probability of a random mutation to occur during cultural transmission.
         features (List[int]): List of features of the individual.
         number_of_changes (int): The number of times this individual has changed set of features following an interaction.
     """
-    def __init__(self, id: int, original_deme_id: int, number_of_features: int, number_of_traits: int, features: List = None):
+    def __init__(self, id: int, original_deme_id: int, number_of_features: int, number_of_traits: int, mutation_rate: float = 0.0, features: List = None):
         """
         Create a new individual with a random set of features.
 
@@ -26,12 +27,14 @@ def __init__(self, id: int, original_deme_id: int, number_of_features: int, numb
             original_deme_id (int): Identifier of the deme where the individual originated.
             number_of_features (int): Number of cultural features of the individual.
             number_of_traits (int): Number of traits per feature of the individual.
+            mutation_rate (float, optional): Probability of a random mutation to occur during cultural transmission.
             features (List, optional): Preset set of features of the individual. Default is None.
         """
         self.id = id
         self.original_deme_id = original_deme_id
         self.number_of_features = number_of_features
         self.number_of_traits = number_of_traits
+        self.mutation_rate = mutation_rate
 
         if features is None:
             # number of traits is +1 as the argument high is exclusive
@@ -43,6 +46,7 @@ def __init__(self, id: int, original_deme_id: int, number_of_features: int, numb
                 raise ValueError("The input number of features does not match the input set of features!")
 
         self.number_of_changes = 0
+        self.number_of_mutations = 0
 
 
     def axelrod_interaction(self, interacting_individual: "Individual") -> None:
@@ -55,21 +59,34 @@ def axelrod_interaction(self, interacting_individual: "Individual") -> None:
             interacting_individual (Individual): Individual with which the self individual interacts. Currently accepts only "axelrod_interaction".
         """
         probability_of_interaction = 1 - np.count_nonzero(self.features - interacting_individual.features)/self.number_of_features
-        random_number = np.random.rand()
-        if (random_number <= probability_of_interaction) and (probability_of_interaction < 1.0):
-            index = np.random.choice(np.nonzero(self.features - interacting_individual.features)[0])
-            self.features[index] = interacting_individual.features[index]
-            self.number_of_changes += 1
+        [interaction_random_number, mutation_random_number] = np.random.rand(2) # generates two random numbers
+        if (interaction_random_number <= probability_of_interaction) and (probability_of_interaction < 1.0):
+            index_to_copy = np.random.choice(np.nonzero(self.features - interacting_individual.features)[0])
+            if mutation_random_number <= self.mutation_rate:
+                # if mutation is occurring, just chose a random trait from possible traits
+                self.features[index_to_copy] = np.random.randint(low = 1, high = self.number_of_traits+1, size=1)
+                self.number_of_mutations += 1
+                self.number_of_changes += 1
+            else:
+                self.features[index_to_copy] = interacting_individual.features[index_to_copy]
+                self.number_of_changes += 1         
+
 
     def neutral_interaction(self, interacting_individual: "Individual") -> None:
         """
         Interaction following a neutral model, where replication of a trait is purely based on frequency in the population. The focal indivdual changes one 
         trait at random copying from the source individual.
         """
-        index = np.random.choice(range(0, self.number_of_features))
-        self.features[index] = interacting_individual.features[index]
-        self.number_of_changes += 1
-            
+        mutation_random_number = np.random.rand()
+        index_to_copy = np.random.choice(range(0, self.number_of_features))
+        if mutation_random_number <= self.mutation_rate:
+            self.features[index_to_copy] = np.random.randint(low = 1, high = self.number_of_trait+1, size=1)
+            self.number_of_mutations += 1
+            self.number_of_changes += 1
+        else:
+            self.features[index_to_copy] = interacting_individual.features[index_to_copy]
+            self.number_of_changes += 1
+
 
     def interact(self, interacting_individual: "Individual", interaction_function: str) -> None:
         """

metapypulation/individual.py

@@ -24,6 +24,7 @@ class Metapopulation():
         carrying_capacities (List[int] | int): A list of carrying capacities (one for each subpopulation) or an integer (same carrying capacity for each subpopulation).
         number_of_features (int): Total number of cultural features per individual.
         number_of_traits (int, optional): Number of different possible traits for each cultural feature.
+        mutation_rate (float, optional): Probability of a mutation to occur.
         min_trait (int, optional): Minimum value for a trait in each feature. 
         max_trait (int, optional): Maximum value for a trait in each feature. 
     """
@@ -33,6 +34,7 @@ def __init__(self, number_of_subpopulations: int,
                  carrying_capacities: List[int] | int = 100,
                  number_of_features: int = 5,
                  number_of_traits: int = 10,
+                 mutation_rate: float = 0.0,
                  min_trait: int = 1,
                  max_trait: int = 10
                  ):
@@ -45,6 +47,7 @@ def __init__(self, number_of_subpopulations: int,
             carrying_capacities (List[int] | int, optional): Either a list of carrying capacities (of which the `len()` is the same as `number_of_subpopulations`) or single integer determining the same carrying capacity for all subpopulations. Defaults to 100.
             number_of_features (int, optional): Total number of cultural features per individual. Defaults to 5.
             number_of_traits (int, optional): Number of different possible traits for each cultural feature. Defaults to 10.
+            mutation_rate (float, optional): Probability of a mutation to occur. Defaults to 0.0.
             min_trait (int, optional): Minimum value for a trait in each feature. Defaults to 1.
             max_trait (int, optional): Maximum value for a trait in each feature. Deafults to 10.
         """
@@ -55,6 +58,7 @@ def __init__(self, number_of_subpopulations: int,
         self.carrying_capacities = carrying_capacities
         self.number_of_features = number_of_features
         self.number_of_traits = number_of_traits
+        self.mutation_rate = mutation_rate
         self.min_trait = min_trait
         self.max_trait = max_trait
 
@@ -72,14 +76,14 @@ def populate(self) -> None:
                     for i in range(self.carrying_capacities[subpopulation.id]):
                         derived_number_of_traits = self.max_trait - self.min_trait + 1 # e.g. if smaller trait is 1 and largest is 10, there are 10 traits: 10 - 1 + 1 
                         set_of_features = np.random.randint(low = self.min_trait, high = self.max_trait + 1, size = self.number_of_features)
-                        new_individual = Individual(i, subpopulation.id, self.number_of_features, derived_number_of_traits, set_of_features)
+                        new_individual = Individual(i, subpopulation.id, self.number_of_features, derived_number_of_traits, self.mutation_rate, set_of_features)
                         subpopulation.add_individual(new_individual)
             case int():
                 for subpopulation in self.subpopulations:
                     for i in range(self.carrying_capacities):
                         derived_number_of_traits = self.max_trait - self.min_trait + 1 # e.g. if smaller trait is 1 and largest is 10, there are 10 traits: 10 - 1 + 1 
                         set_of_features = np.random.randint(low = self.min_trait, high = self.max_trait + 1, size = self.number_of_features)
-                        new_individual = Individual(i, subpopulation.id, self.number_of_features, derived_number_of_traits, set_of_features)
+                        new_individual = Individual(i, subpopulation.id, self.number_of_features, derived_number_of_traits, self.mutation_rate, set_of_features)
                         subpopulation.add_individual(new_individual)
 
 

metapypulation/metapopulation.py

@@ -142,53 +142,6 @@ def run_single_replicate(self, replicate_id: int) -> None:
 
             print(f"{t} generations ran in {total_time}.")
 
-            
-    # def run_replicate_with_burn_in(self, replicate_id: int) -> None:
-    #     """
-    #     Run one replicate of the simulation.
-
-    #     Args:
-    #         replicate_id (int): The number of the current replicate (for the output data columns).
-    #     """
-    #     metapopulation = Metapopulation(self.number_of_subpopulations, self.interaction_type, self.migration_matrix, self.carrying_capacities)
-    #     metapopulation.populate()
-        
-    #     set_counts = []
-    #     shannon = []
-    #     metapop_counts = []
-    #     metapop_shannon = []
-        
-    #     start_time = time.time()
-        
-    #     for t in range(self.generations + 1):
-    #         if self.verbose:
-    #             if t%self.verbose_timing == 0:
-    #                 print(f"Replicate {replicate_id}, gen {t}!")
-    #                 # TODO print other fun stuff
-                    
-    #         if t%self.measure_timing == 0:
-    #             set_counts.append(np.mean(metapopulation.traits_sets_per_subpopulation()))
-    #             shannon.append(np.mean(metapopulation.shannon_diversity_per_subpopulation()))
-    #             metapop_counts.append(metapopulation.metapopulation_test_sets())
-    #             metapop_shannon.append(metapopulation.metapopulation_shannon_diversity())
-            
-    #         if t > self.burn_in:
-    #             metapopulation.migrate()
-    
-    #         metapopulation.make_interact()
-        
-    #     self.subpop_set_counts = pd.concat([self.subpop_set_counts, pd.Series(set_counts, name=replicate_id)], axis=1)
-    #     self.subpop_shannon = pd.concat([self.subpop_shannon, pd.Series(shannon, name=replicate_id)], axis=1)
-    #     self.metapop_set_counts = pd.concat([self.metapop_set_counts, pd.Series(metapop_counts, name=replicate_id)], axis=1)
-    #     self.metapop_shannoneturn  = pd.concat([self.metapop_shannon, pd.Series(metapop_shannon, name=replicate_id)], axis=1)
-                             
-    #     if self.verbose:
-    #         end_time = time.time()
-    #         total_time = end_time - start_time
-    #         total_time = time.strftime("%H:%M:%S", time.gmtime(total_time))
-
-    #         print(f"{t} generations ran in {total_time}.")
-           
 
     def run_simulation(self) -> None:
         """

metapypulation/simulation.py

@@ -0,0 +1,49 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import time
+
+from metapypulation.simulation import Simulation
+
+total_population = 200
+replicates = 1
+#migrations = ['island'] # 'stepping_stone'
+interactions = ['neutral_interaction'] # 'axelrod_interaction'
+#number_of_subpopulations = [4]
+subpopulations = 4
+population_sizes = [[50, 50, 50, 50]]
+
+start_time = time.time()
+
+count = 0
+
+for interaction in interactions:
+    for rate_of_migration in [0.001, 0.1]:
+        # rate_of_migration = 0.001
+
+        carrying_capacity = int(np.ceil(total_population / subpopulations)) # [283, 39, 39, 39]# 
+        generations = 200000 #350000
+        burn_in = 0
+
+        # create one big simulation
+        simulation = Simulation(generations = generations, 
+                                number_of_subpopulations=subpopulations, 
+                                migration_matrix = migration, 
+                                interaction = interaction,
+                                carrying_capacities = carrying_capacity,
+                                replicates = replicates,
+                                output_path = f"./Outputs/{subpopulations}subpop_popsize{population_sizes[0]}_{migration}_{interaction}_{rate_of_migration}_noburnin",  #TAG2024/04-migration-rates/{subpopulations}subpop_{migration}_{interaction}_m{rate_of_migration}",
+                                burn_in = burn_in,
+                                migration_rate = rate_of_migration)
+
+        # modify individuals in simulation
+        
+
+        simulation.run_simulation()
+        count = count + 1
+    #f"./Outputs/SourceSink/pop{total_population}/{migrations}/{subpopulations}subpop_m1e-3_burnin_{carrying_capacity[0]}",
+                            
+end_time = time.time() - start_time
+hours = round(end_time//3600)
+minutes = round(end_time//60) - hours*60
+seconds = round(end_time) - hours*3600 - minutes*60
+print(f"Simulation of {count} sets of parameters, {replicates} replicates each, finished in {hours}h, {minutes}m and {seconds}s")