Merge pull request #92 from andrewkern/generation_time_addition

added generation times as class variable to demog models

dmel_example.py

@@ -7,7 +7,8 @@
 import stdpopsim
 from stdpopsim import drosophila_melanogaster
 
-chrom = drosophila_melanogaster.genome.chromosomes["chr2R"]
+chr_str = "chr2L"
+chrom = drosophila_melanogaster.genome.chromosomes[chr_str]
 recomb_map = chrom.recombination_map()
 
 print("Testing Li and Stephan (2006) model")
@@ -28,6 +29,21 @@
     **model.asdict())
 print("simulated:", ts.num_trees, ts.num_sites)
 
+print("====================\n")
+print(f"Testing Li and Stephan (2006) model- with on {chr_str} map")
+model = drosophila_melanogaster.LiStephanTwoPopulation()
+model.debug()
+
+samples = [msprime.Sample(population=0, time=0)] * 10 + [msprime.Sample(population=1, time=0)] * 10
+
+
+ts = msprime.simulate(
+    samples=samples,
+    recombination_map=chrom.recombination_map(),
+    mutation_rate=chrom.default_mutation_rate,
+    **model.asdict())
+print("simulated:", ts.num_trees, ts.num_sites)
+
 print("====================\n")
 print("Testing Sheehan and Song (2016) model")
 model = drosophila_melanogaster.SheehanSongThreeEpoch()
@@ -47,7 +63,7 @@
 print("simulated:", ts.num_trees, ts.num_sites)
 
 print("====================\n")
-print("Testing Sheehan and Song (2016) model on chrom2R map.")
+print(f"Testing Sheehan and Song (2016) model on {chr_str} map.")
 print("This will take a while- perhaps 2 hours\n")
 model = drosophila_melanogaster.SheehanSongThreeEpoch()
 

stdpopsim/arabidopsis_thaliana.py

@@ -74,6 +74,8 @@ class Salome2012(genetic_maps.GeneticMap):
 #
 ###########################################################
 
+default_generation_time = 1.0
+
 
 class Durvasula2017MSMC(models.Model):
     """
@@ -102,6 +104,8 @@ def __init__(self):
         # set the last 2 entries equal
         # to the size at 30 (~1.6Mya)
         self.sizes[30:32] = self.sizes[30]
+        # generation time is 1 year
+        self.generation_time = default_generation_time
         self.demographic_events = []
         for idx, t in enumerate(self.times):
             self.demographic_events.append(

stdpopsim/drosophila_melanogaster.py

@@ -72,6 +72,8 @@ class Comeron2012_dm6(genetic_maps.GeneticMap):
 #
 ###########################################################
 
+default_generation_time = 0.1
+
 
 class SheehanSongThreeEpoch(models.Model):
     """
@@ -98,6 +100,7 @@ def __init__(self):
         # generation_time = 10 / year
         t_1 = t_1_coal * 4 * N_ref
         t_2 = (t_1_coal + t_2_coal) * 4 * N_ref
+        self.generation_time = default_generation_time
         # Single population in this model
         self.population_configurations = [
             msprime.PopulationConfiguration(initial_size=N_R),
@@ -131,6 +134,7 @@ def __init__(self):
         N_A0 = 8.603e06
         t_A0 = 600000  # assuming 10 generations / year
         N_A1 = N_A0 / 5.0
+        self.generation_time = default_generation_time
         # European Parameter values from "Demo History of Euro Population"
         N_E0 = 1.075e06
         N_E1 = 2200

stdpopsim/e_coli.py

@@ -36,6 +36,8 @@
 #
 ###########################################################
 
+# TODO add a generation time here
+# default_generation_time = -1
 
 class LapierreConstant(models.Model):
     """
@@ -44,6 +46,7 @@ class LapierreConstant(models.Model):
     """
 
     def __init__(self):
+        super().__init__()
 
         N_e = 1.8e8
         # Single population

stdpopsim/homo_sapiens.py

@@ -162,6 +162,10 @@ def chromosome_factory(name, genetic_map=None, length_multiplier=1):
 ###########################################################
 
 
+# species wide default generation time
+default_generation_time = 25
+
+
 class GutenkunstThreePopOutOfAfrica(models.Model):
     """
     Model Name:
@@ -208,10 +212,10 @@ def __init__(self):
         N_EU0 = 1000
         N_AS0 = 510
         # Times are provided in years, so we convert into generations.
-        generation_time = 25
-        T_AF = 220e3 / generation_time
-        T_B = 140e3 / generation_time
-        T_EU_AS = 21.2e3 / generation_time
+        self.generation_time = default_generation_time
+        T_AF = 220e3 / self.generation_time
+        T_B = 140e3 / self.generation_time
+        T_EU_AS = 21.2e3 / self.generation_time
         # We need to work out the starting (diploid) population sizes based on
         # the growth rates provided for these two populations
         r_EU = 0.004
@@ -297,11 +301,11 @@ class TennessenTwoPopOutOfAfrica(models.Model):
     def __init__(self):
         super().__init__()
 
-        generation_time = 25
-        T_AF = 148e3 / generation_time
-        T_OOA = 51e3 / generation_time
-        T_EU0 = 23e3 / generation_time
-        T_EG = 5115 / generation_time
+        self.generation_time = default_generation_time
+        T_AF = 148e3 / self.generation_time
+        T_OOA = 51e3 / self.generation_time
+        T_EU0 = 23e3 / self.generation_time
+        T_EG = 5115 / self.generation_time
 
         # Growth rates
         r_EU0 = 0.00307
@@ -392,9 +396,9 @@ class TennessenOnePopAfrica(models.Model):
     def __init__(self):
         super().__init__()
 
-        generation_time = 25
-        T_AF = 148e3 / generation_time
-        T_EG = 5115 / generation_time
+        self.generation_time = default_generation_time
+        T_AF = 148e3 / self.generation_time
+        T_EG = 5115 / self.generation_time
 
         # Growth rate
         r_AF = 0.0166

stdpopsim/models.py

@@ -124,6 +124,7 @@ def __init__(self):
         self.demographic_events = []
         # Defaults to a single population
         self.migration_matrix = [[0]]
+        self.generation_time = -1
 
     @property
     def name(self):

stdpopsim/pongo.py

@@ -11,13 +11,15 @@
 # TODO: how do we define the Orangutan genome here? Are they similar enough to
 # humans that we just use humans? What about the different species of pongo?
 
+# TODO: add a default generation time to the species
 
 class LockeEtAlPongoIM(models.Model):
     '''
     http://doi.org/10.1038/nature09687
     '''
 
     def __init__(self):
+        super().__init__()
 
         # Parameters from paper:
         # ancestral size, before split

tests/test_models.py

@@ -247,6 +247,10 @@ def test_filtering_outside_classes(self):
         for model in models.all_models():
             self.assertNotIsInstance(model, DummyModel)
 
+    def test_generation_times_non_empty(self):
+        self.assertGreater(len([model.generation_time for model in
+                                models.all_models()]), 0)
+
 
 class TestModelsEqual(unittest.TestCase):
     """
@@ -291,3 +295,12 @@ def test_migration_matrices(self):
         self.assertTrue(m1.equals(m2))
         # If we have higher tolerances we catch the differences
         self.assertFalse(m1.equals(m2, atol=1e-10, rtol=1e-9))
+
+
+class TestModelProperties(unittest.TestCase):
+    def test_model_generation_time(self):
+        self.assertTrue(models.Model().generation_time == -1)
+        known_models = models.all_models()
+        n = len(known_models)
+        for j in range(n):
+            self.assertTrue(known_models[j].generation_time > -2)