Repository: popsim-consortium/stdpopsim (language: python)

Issue:
HomSap/Zigzag_1S14 population sizes don't match Schiffels & Durbin (2014)
**Short description of the problem:**
The population sizes in the `Zigzag_1S14` demographic model swing between 1431 and 14312. However, the values in [Schiffels & Durbin (2014)](https://doi.org/10.1038/ng.3015) are five times higher.

**How to reproduce the problem:**
Compare the zigzag parameters in [stdpopsim's catalog](https://stdpopsim.readthedocs.io/en/latest/catalog.html#sec_catalog_homsap_models_zigzag_1s14) with Figure 2 in Schiffels & Durbin (2014). Also see the `ms` command in the supplementary material, section 7, pg. 25.

This appears to have been deliberate, according to #105, but is not documented, and doesn't seem to have come up during the model QC process (#498).

**Suggested fixes:**
The original `ms` command given in the supplement includes a `-eN 0 5` command, which has not been reproduced in the stdpopsim model. From the ms manual:
`-eN t x` `Set all subpop’s to size x ∗ N0 and growth rates to zero.`

We should at the very least document the difference between the stdpopsim ZigZag model, and that in the original publication.

CC: @Chris1221, @apragsdale.

[EDIT: [translation of the ms command, as stolen from the QC implementation](https://gist.github.com/grahamgower/7a2c0115493547e5ddc02887030abee9)]


Discussion / hints:
Yep, that looks like a mistake on our part. We should probably fix those values in our implementation to follow the original publication, instead of just documenting it. Maybe the next major release could emit a warning the first time someone uses the zigzag model, that the previous version of stdpopsim had sizes that were 5x smaller than intended?

I'm sorry for missing that in the QC. I'll be very happy the day that I never have to translate an ms command again.
Happy day indeed! It'll come: https://github.com/popsim-consortium/demes-python/pull/102

(more quickly, if someone would like to pick this PR up...)
I agree, that does look like a mistake in our implementation. My apologies for missing that initially, my ms-fu is pretty weak and I completely didn't see the scaling factor. Comparing it to the literature values would have been a good sanity check that I also missed... I agree with @apragsdale's warning suggestion, it's important that people know that it was incorrect!  
Also, (I think) it looks like this change may need to be propagated to the values in the Zigzag model description in the demes-docs gallery: https://popsim-consortium.github.io/demes-docs/main/gallery.html#zigzag-model
> Also, (I think) it looks like this change may need to be propagated to the values in the Zigzag model description in the demes-docs gallery: https://popsim-consortium.github.io/demes-docs/main/gallery.html#zigzag-model

Yep, thanks @Chris1221, I've opened an issues in Demes for this too: https://github.com/popsim-consortium/demes-python/issues/179.

> Also, (I think) it looks like this change may need to be propagated to the values in the Zigzag model description in the demes-docs gallery: https://popsim-consortium.github.io/demes-docs/main/gallery.html#zigzag-model

Yes, for sure. I think I used our implementation from stdpopsim to translate over to demes. Great example of how easily these errors can propagate.. similar to the OOA model issue from earlier this year. We'll need fix it there as well.
Great to catch this issue, but I don't think it's that serious in the scheme of things since it's a toy model. It's just another example of how easy it is to make mistakes with today's toolchain. So, no worries @Chris1221 and @apragsdale!

In terms of compatibility, I'm not sure what to do. I guess emitting a warning when it's used is probably the best we can do, which we remove after (say) 6 months.

Fix patch:
diff --git a/stdpopsim/engines.py b/stdpopsim/engines.py
index 6c2e7be93..a99c6d06c 100644
--- a/stdpopsim/engines.py
+++ b/stdpopsim/engines.py
@@ -104,15 +104,6 @@ def get_version(self):
         """
         raise NotImplementedError()
 
-    def _warn_zigzag(self, demographic_model):
-        if demographic_model.id == "Zigzag_1S14":
-            warnings.warn(
-                "In stdpopsim <= 0.1.2, the Zigzag_1S14 model produced population "
-                "sizes 5x lower than those in Schiffels & Durbin (2014). "
-                "The population sizes have now been corrected. For details see: "
-                "https://github.com/popsim-consortium/stdpopsim/issues/745"
-            )
-
     def _warn_mutation_rate_mismatch(self, contig, demographic_model):
         if demographic_model.mutation_rate is not None and not math.isclose(
             demographic_model.mutation_rate, contig.mutation_rate
@@ -271,8 +262,6 @@ def simulate(
                 "Samples must be a dict of the form {population_name:num_samples}."
             )
 
-        # TODO: remove this after a release or two. See #745.
-        self._warn_zigzag(demographic_model)
         self._warn_mutation_rate_mismatch(contig, demographic_model)
         self._warn_recombination_rate_mismatch(contig, demographic_model)
 
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index f8aa7129d..8e9568d50 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -1625,8 +1625,6 @@ def simulate(
         # by msprime later
         slim_rate_map = get_slim_mutation_rate_map(contig)
 
-        # TODO: remove this after a release or two. See #745.
-        self._warn_zigzag(demographic_model)
         self._warn_mutation_rate_mismatch(contig, demographic_model)
         self._warn_recombination_rate_mismatch(contig, demographic_model)
 


Test patch:
diff --git a/tests/test_models.py b/tests/test_models.py
index 291be6fff..31ee902f5 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -369,17 +369,6 @@ def test_population_order(self):
             assert sum(allow_sample_status[num_sampling:]) == 0
 
 
-class TestZigZagWarning:
-    def test_zigzag_produces_warning(self):
-        species = stdpopsim.get_species("HomSap")
-        model = species.get_demographic_model("Zigzag_1S14")
-        contig = species.get_contig("chr22")
-        samples = {"generic": 5}
-        for engine in stdpopsim.all_engines():
-            with pytest.warns(UserWarning, match="Zigzag_1S14 model.*sizes 5x lower"):
-                engine.simulate(model, contig, samples, dry_run=True)
-
-
 class TestMutationRates:
     def test_mutation_rate_warning(self):
         species = stdpopsim.get_species("HomSap")