Merge pull request #22 from mcveanlab/port-tennyson

Quickly port Tennessen et al model.
popsim-consortium · Mar 5, 2019 · 20fb835 · 20fb835
2 parents 79d26b4 + b3d45ea
commit 20fb835
Show file tree

Hide file tree

Showing 2 changed files with 64 additions and 1 deletion.
diff --git a/stdpopsim/homo_sapiens.py b/stdpopsim/homo_sapiens.py
@@ -106,7 +106,7 @@ class GutenkunstThreePopOutOfAfrica(models.Model):
     """
 
     def __init__(self):
-
+        super().__init__()
         # First we set out the maximum likelihood values of the various parameters
         # given in Table 1.
         N_A = 7300
@@ -161,3 +161,46 @@ def __init__(self):
             msprime.PopulationParametersChange(
                 time=T_AF, initial_size=N_A, population_id=0)
         ]
+
+
+class TennessenEuropean(models.Model):
+    """
+    The model is derived from the Tennesen et al.
+    `analysis <https://doi.org/10.1126/science.1219240>`_  of the jSFS from
+    European Americans and African Americans.
+
+    .. todo:: document this model, including the original publications
+        and clear information about what the different population indexes
+        mean.
+
+    """
+    def __init__(self):
+        super().__init__()
+        # Population sizes
+        N_A = 7310
+        N_AF = 14474
+        N_B = 1861
+        N_EU1 = 9475
+        # Times
+        generation_time = 25
+        T_AF = 148000 / generation_time
+        T_B = 51000 / generation_time
+        T_EU0 = 23000 / generation_time
+        T_EU1 = 5115 / generation_time
+        # Rates and Present Ne
+        r_EU0 = 0.00307
+        r_EU1 = 0.0195
+        N_EU = N_EU1 / math.exp(-r_EU1 * T_EU1)
+        self.population_configurations = [
+            msprime.PopulationConfiguration(initial_size=N_EU, growth_rate=r_EU1)
+        ]
+        self.demographic_events = [
+            msprime.PopulationParametersChange(
+                time=T_EU1, initial_size=N_EU1, growth_rate=r_EU0, population_id=0),
+            msprime.PopulationParametersChange(
+                time=T_EU0, initial_size=N_B, growth_rate=0, population_id=0),
+            msprime.PopulationParametersChange(
+                time=T_B, initial_size=N_AF, growth_rate=0, population_id=0),
+            msprime.PopulationParametersChange(
+                time=T_AF, initial_size=N_A, growth_rate=0, population_id=0)
+        ]
diff --git a/tests/test_homo_sapiens.py b/tests/test_homo_sapiens.py
@@ -72,3 +72,23 @@ def test_debug_runs(self):
         model.debug(output)
         s = output.getvalue()
         self.assertGreater(len(s), 0)
+
+
+class TestTennessenEuropean(unittest.TestCase):
+    """
+    Basic tests for the TennessenEuropean model.
+    """
+
+    def test_simulation_runs(self):
+        model = homo_sapiens.TennessenEuropean()
+        ts = msprime.simulate(
+            samples=[msprime.Sample(0, 0) for _ in range(10)],
+            **model.asdict())
+        self.assertEqual(ts.num_populations, 1)
+
+    def test_debug_runs(self):
+        model = homo_sapiens.TennessenEuropean()
+        output = io.StringIO()
+        model.debug(output)
+        s = output.getvalue()
+        self.assertGreater(len(s), 0)