From b8c2329378acd97c25ac7007be8cea9d4220f720 Mon Sep 17 00:00:00 2001
From: ckyriazis <CKyriazis33@gmail.com>
Date: Wed, 25 Oct 2023 12:21:58 -0700
Subject: [PATCH] Vaquita QC

---
 stdpopsim/qc/PhoSin.py   | 45 +++++++++++++++++++
 stdpopsim/qc/__init__.py |  1 +
 tests/test_PhoSin.py     | 96 +++++++++++++++++++---------------------
 3 files changed, 92 insertions(+), 50 deletions(-)
 create mode 100644 stdpopsim/qc/PhoSin.py

diff --git a/stdpopsim/qc/PhoSin.py b/stdpopsim/qc/PhoSin.py
new file mode 100644
index 000000000..28c0e25b8
--- /dev/null
+++ b/stdpopsim/qc/PhoSin.py
@@ -0,0 +1,45 @@
+import msprime
+
+import stdpopsim
+
+
+_species = stdpopsim.get_species("PhoSin")
+
+
+def Robinson2022_TwoEpoch():
+    id = "QC-Vaquita2Epoch_1R22"
+
+    populations = [
+        stdpopsim.Population(id="Vaquita", description=""),
+    ]
+
+    # Time of second epoch
+    T_2 = 2162
+    # population sizes
+    N_ANC = 4485
+    N_2 = 2807
+
+    return stdpopsim.DemographicModel(
+        id=id,
+        description=id,
+        long_description=id,
+        generation_time=_species.generation_time,
+        mutation_rate=5.83e-9,
+        populations=populations,
+        population_configurations=[
+            msprime.PopulationConfiguration(
+                initial_size=N_2, metadata=populations[0].asdict()
+            ),
+        ],
+        demographic_events=[
+            msprime.PopulationParametersChange(
+                time=T_2, initial_size=N_ANC, population_id=0
+            ),
+        ],
+        population_id_map=[{"Vaquita": 0}] * 2,
+    )
+
+
+_species.get_demographic_model("Vaquita2Epoch_1R22").register_qc(
+    Robinson2022_TwoEpoch()
+)
diff --git a/stdpopsim/qc/__init__.py b/stdpopsim/qc/__init__.py
index 20b4013b2..f9a76cb9f 100644
--- a/stdpopsim/qc/__init__.py
+++ b/stdpopsim/qc/__init__.py
@@ -9,3 +9,4 @@
 from . import PanTro  # NOQA
 from . import OrySat  # NOQA
 from . import MusMus  # NOQA
+from . import PhoSin  # NOQA
diff --git a/tests/test_PhoSin.py b/tests/test_PhoSin.py
index be93fe6ff..24d0864f7 100644
--- a/tests/test_PhoSin.py
+++ b/tests/test_PhoSin.py
@@ -21,45 +21,42 @@ def test_common_name(self):
     # independently referring to the citations provided in the
     # species definition, filling in the appropriate values
     # and deleting the pytest "skip" annotations.
-    @pytest.mark.skip("Population size QC not done yet")
     def test_qc_population_size(self):
-        assert self.species.population_size == -1
+        assert self.species.population_size == 3500
 
-    @pytest.mark.skip("Generation time QC not done yet")
     def test_qc_generation_time(self):
-        assert self.species.generation_time == -1
+        assert self.species.generation_time == 11.9
 
 
 class TestGenomeData(test_species.GenomeTestBase):
 
     genome = stdpopsim.get_species("PhoSin").genome
 
-    @pytest.mark.skip("Recombination rate QC not done yet")
     @pytest.mark.parametrize(
         ["name", "rate"],
         {
-            "1": -1,
-            "2": -1,
-            "3": -1,
-            "4": -1,
-            "5": -1,
-            "6": -1,
-            "7": -1,
-            "8": -1,
-            "9": -1,
-            "10": -1,
-            "11": -1,
-            "12": -1,
-            "13": -1,
-            "14": -1,
-            "15": -1,
-            "16": -1,
-            "17": -1,
-            "18": -1,
-            "19": -1,
-            "20": -1,
-            "21": -1,
-            "X": -1,
+            "1": 1e-8,
+            "2": 1e-8,
+            "3": 1e-8,
+            "4": 1e-8,
+            "5": 1e-8,
+            "6": 1e-8,
+            "7": 1e-8,
+            "8": 1e-8,
+            "9": 1e-8,
+            "10": 1e-8,
+            "11": 1e-8,
+            "12": 1e-8,
+            "13": 1e-8,
+            "14": 1e-8,
+            "15": 1e-8,
+            "16": 1e-8,
+            "17": 1e-8,
+            "18": 1e-8,
+            "19": 1e-8,
+            "20": 1e-8,
+            "21": 1e-8,
+            "X": 1e-8,
         }.items(),
     )
     def test_recombination_rate(self, name, rate):
@@ -67,32 +64,31 @@ def test_recombination_rate(self, name, rate):
             self.genome.get_chromosome(name).recombination_rate
         )
 
-    @pytest.mark.skip("Mutation rate QC not done yet")
     @pytest.mark.parametrize(
         ["name", "rate"],
         {
-            "1": -1,
-            "2": -1,
-            "3": -1,
-            "4": -1,
-            "5": -1,
-            "6": -1,
-            "7": -1,
-            "8": -1,
-            "9": -1,
-            "10": -1,
-            "11": -1,
-            "12": -1,
-            "13": -1,
-            "14": -1,
-            "15": -1,
-            "16": -1,
-            "17": -1,
-            "18": -1,
-            "19": -1,
-            "20": -1,
-            "21": -1,
-            "X": -1,
+            "1": 5.83e-9,
+            "2": 5.83e-9,
+            "3": 5.83e-9,
+            "4": 5.83e-9,
+            "5": 5.83e-9,
+            "6": 5.83e-9,
+            "7": 5.83e-9,
+            "8": 5.83e-9,
+            "9": 5.83e-9,
+            "10": 5.83e-9,
+            "11": 5.83e-9,
+            "12": 5.83e-9,
+            "13": 5.83e-9,
+            "14": 5.83e-9,
+            "15": 5.83e-9,
+            "16": 5.83e-9,
+            "17": 5.83e-9,
+            "18": 5.83e-9,
+            "19": 5.83e-9,
+            "20": 5.83e-9,
+            "21": 5.83e-9,
+            "X": 5.83e-9,
         }.items(),
     )
     def test_mutation_rate(self, name, rate):