fix split model deme overlap

doc/misc/changelog.md

@@ -21,6 +21,11 @@ Bug fixes
   leading to runtime exceptions.
   PR {pr}`802`
   PR {pr}`803`
+* Fix error in `demes` models where "replacement" models had 1 generation of overlap between ancestral/derived demes.
+  Issue {issue}`814`
+  PR {pr}`815`
+  {user}`apragsdale`
+  {user}`molpopgen`
 
 Behavior changes
 

fwdpy11/_functions/import_demes.py

@@ -1,35 +1,23 @@
-from typing import (
-    Dict,
-    List,
-    Optional,
-    Union,
-)
-import numpy as np
-import attr
-import math
-import sys
 import copy
 import itertools
+import math
+import sys
+from typing import Dict, List, Optional, Union
 
+import attr
 import demes
-
-from ..discrete_demography import (
-    MassMigration,
-    MigrationMatrix,
-    SetDemeSize,
-    SetExponentialGrowth,
-    SetMigrationRates,
-    SetSelfingRate,
-    copy_individuals,
-    move_individuals,
-    DiscreteDemography,
-)
+import numpy as np
 
 from .. import class_decorators
-
-from ..demographic_models import DemographicModelDetails, DemographicModelCitation
-
 from .._demography import exponential_growth_rate
+from ..demographic_models import (DemographicModelCitation,
+                                  DemographicModelDetails)
+from ..discrete_demography import (DiscreteDemography, MassMigration,
+                                   MigrationMatrix, SetDemeSize,
+                                   SetExponentialGrowth, SetMigrationRates,
+                                   SetSelfingRate, copy_individuals,
+                                   move_individuals)
+
 
 # TODO: need type hints for dg
 def demography_from_demes(
@@ -451,7 +439,7 @@ def _process_epoch(
         # Handle size change functions
         events.set_deme_sizes.append(
             SetDemeSize(
-                when=when - 1,
+                when=when,
                 deme=idmap[deme_id],
                 new_size=e.start_size,
             )
@@ -575,8 +563,8 @@ def _process_migrations(
                         destination=idmap[m.dest],
                         rate_change=-m.rate,
                         from_deme_graph=True,
-                        )
                     )
+                )
             except AttributeError:
                 for source, dest in itertools.permutations(m.demes, 2):
                     events.migration_rate_changes.append(
@@ -601,7 +589,7 @@ def _process_pulses(
         when = burnin_generation + int(model_times.model_start_time - p.time)
         events.migration_rate_changes.append(
             _MigrationRateChange(
-                when=when - 1,
+                when=when,
                 source=idmap[p.source],
                 destination=idmap[p.dest],
                 rate_change=p.proportion,
@@ -610,7 +598,7 @@ def _process_pulses(
         )
         events.migration_rate_changes.append(
             _MigrationRateChange(
-                when=when,
+                when=when + 1,
                 source=idmap[p.source],
                 destination=idmap[p.dest],
                 rate_change=-p.proportion,
@@ -632,7 +620,7 @@ def _process_admixtures(
         for parent, proportion in zip(a.parents, a.proportions):
             events.migration_rate_changes.append(
                 _MigrationRateChange(
-                    when=when - 1,
+                    when=when,
                     source=idmap[parent],
                     destination=idmap[a.child],
                     rate_change=proportion,
@@ -641,7 +629,7 @@ def _process_admixtures(
             )
             events.migration_rate_changes.append(
                 _MigrationRateChange(
-                    when=when,
+                    when=when + 1,
                     source=idmap[parent],
                     destination=idmap[a.child],
                     rate_change=-proportion,
@@ -663,7 +651,7 @@ def _process_mergers(
         for parent, proportion in zip(m.parents, m.proportions):
             events.migration_rate_changes.append(
                 _MigrationRateChange(
-                    when=when - 1,
+                    when=when,
                     source=idmap[parent],
                     destination=idmap[m.child],
                     rate_change=proportion,
@@ -672,7 +660,7 @@ def _process_mergers(
             )
             events.migration_rate_changes.append(
                 _MigrationRateChange(
-                    when=when,
+                    when=when + 1,
                     source=idmap[parent],
                     destination=idmap[m.child],
                     rate_change=-proportion,
@@ -703,7 +691,7 @@ def _process_splits(
             # one generation of migration to move lineages from parent to children
             events.migration_rate_changes.append(
                 _MigrationRateChange(
-                    when=when - 1,
+                    when=when,
                     source=idmap[s.parent],
                     destination=idmap[c],
                     rate_change=1.0,
@@ -713,7 +701,7 @@ def _process_splits(
             # turn off that migration after one generation
             events.migration_rate_changes.append(
                 _MigrationRateChange(
-                    when=when,
+                    when=when + 1,
                     source=idmap[s.parent],
                     destination=idmap[c],
                     rate_change=-1.0,
@@ -742,7 +730,7 @@ def _process_branches(
         # turn on migration for one generation at "when"
         events.migration_rate_changes.append(
             _MigrationRateChange(
-                when=when - 1,
+                when=when,
                 source=idmap[b.parent],
                 destination=idmap[b.child],
                 rate_change=1.0,
@@ -752,7 +740,7 @@ def _process_branches(
         # end that migration after one generation
         events.migration_rate_changes.append(
             _MigrationRateChange(
-                when=when,
+                when=when + 1,
                 source=idmap[b.parent],
                 destination=idmap[b.child],
                 rate_change=-1.0,

tests/test_demes2fwdpy11.py

@@ -1,12 +1,11 @@
 import copy
+import typing
 import unittest
 from dataclasses import dataclass
 
 import demes
 import fwdpy11
 import numpy as np
-import unittest
-import typing
 import pytest
 
 
@@ -78,9 +77,9 @@ def setUpClass(self):
 
     def test_size_change_params(self):
         self.assertTrue(len(self.demog.model.set_deme_sizes) == 1)
-        self.assertTrue(
-            self.demog.model.set_deme_sizes[0].when
-            == self.demog.metadata["burnin_time"] - 1
+        self.assertEqual(
+            self.demog.model.set_deme_sizes[0].when,
+            self.demog.metadata["burnin_time"],
         )
         self.assertTrue(self.demog.model.set_deme_sizes[0].new_size == 2000)
 
@@ -105,7 +104,7 @@ def setUpClass(self):
     def test_conversion_to_generations(self):
         self.assertTrue(
             self.demog.model.set_deme_sizes[0].when
-            == self.demog.metadata["burnin_time"] - 1
+            == self.demog.metadata["burnin_time"]
         )
         self.assertTrue(
             self.demog.metadata["total_simulation_length"]
@@ -173,25 +172,30 @@ def setUpClass(self):
         self.demog = fwdpy11.discrete_demography.from_demes(self.g, 10)
 
     def test_size_changes(self):
-        self.assertTrue(len(self.demog.model.set_deme_sizes) == 3)
-        self.assertTrue(self.demog.model.set_deme_sizes[0].deme == 1)
-        self.assertTrue(self.demog.model.set_deme_sizes[0].new_size == 100)
-        self.assertTrue(
-            self.demog.model.set_deme_sizes[0].when
-            == self.demog.metadata["burnin_time"] - 1
-        )
-        self.assertTrue(self.demog.model.set_deme_sizes[1].deme == 2)
-        self.assertTrue(self.demog.model.set_deme_sizes[1].new_size == 100)
-        self.assertTrue(
-            self.demog.model.set_deme_sizes[1].when
-            == self.demog.metadata["burnin_time"] - 1
-        )
-        self.assertTrue(self.demog.model.set_deme_sizes[2].deme == 0)
-        self.assertTrue(self.demog.model.set_deme_sizes[2].new_size == 0)
-        self.assertTrue(
-            self.demog.model.set_deme_sizes[2].when
-            == self.demog.metadata["burnin_time"]
-        )
+        self.assertEqual(len(self.demog.model.set_deme_sizes), 3)
+
+        check_debugger_passes(self.demog)
+
+        # NOTE: commented these out as they test internal sort order and not
+        # the model validity
+        # self.assertTrue(self.demog.model.set_deme_sizes[0].deme == 1)
+        # self.assertTrue(self.demog.model.set_deme_sizes[0].new_size == 100)
+        # self.assertTrue(
+        #     self.demog.model.set_deme_sizes[0].when
+        #     == self.demog.metadata["burnin_time"]
+        # )
+        # self.assertTrue(self.demog.model.set_deme_sizes[1].deme == 2)
+        # self.assertTrue(self.demog.model.set_deme_sizes[1].new_size == 100)
+        # self.assertTrue(
+        #     self.demog.model.set_deme_sizes[1].when
+        #     == self.demog.metadata["burnin_time"]
+        # )
+        # self.assertTrue(self.demog.model.set_deme_sizes[2].deme == 0)
+        # self.assertTrue(self.demog.model.set_deme_sizes[2].new_size == 0)
+        # self.assertTrue(
+        #     self.demog.model.set_deme_sizes[2].when
+        #     == self.demog.metadata["burnin_time"]
+        # )
 
 
 @check_valid_demography
@@ -512,11 +516,11 @@ def test_pulse_migration_matrix(self):
         self.assertTrue(len(self.demog.model.set_migration_rates) == 2)
         self.assertTrue(
             self.demog.model.set_migration_rates[0].when
-            == self.demog.metadata["burnin_time"] - 1
+            == self.demog.metadata["burnin_time"]
         )
         self.assertTrue(
             self.demog.model.set_migration_rates[1].when
-            == self.demog.metadata["burnin_time"]
+            == self.demog.metadata["burnin_time"] + 1
         )
         self.assertTrue(self.demog.model.set_migration_rates[0].deme == 1)
         self.assertTrue(self.demog.model.set_migration_rates[1].deme == 1)
@@ -777,6 +781,7 @@ def test_yamls_with_migration(data):
     demog = fwdpy11.discrete_demography.from_demes(g, 1)
     check_debugger_passes(demog)
 
+
 def test_split_model_population_size_history(two_deme_split_with_ancestral_size_change):
     """
     This is a detailed test of the complete size history
@@ -822,19 +827,21 @@ def __call__(self, pop, _):
     # The daughter demes are seen from 110 till the end
     for deme in [1, 2]:
         assert [i.when for i in recorder.sizes[deme]] == [
-            i for i in range(110, model.metadata["total_simulation_length"] + 1)
+            i for i in range(111, model.metadata["total_simulation_length"] + 1)
         ]
     # initial daughter deme sizes
-    assert recorder.sizes[1][0].size == 250
-    assert recorder.sizes[2][0].size == 50
+    g1 = 2 ** (1 / 10) - 1
+    g2 = 4 ** (1 / 10) - 1
+    assert recorder.sizes[1][0].size == int(np.rint(250 * (1 + g1)))
+    assert recorder.sizes[2][0].size == int(np.rint(50 * (1 + g2)))
     # final daughter deme sizes
     assert recorder.sizes[1][-1].size == 500
     assert recorder.sizes[2][-1].size == 200
 
     # At generation 100, the ancestral pop size changed from 100
     # to 200
     for i in recorder.sizes[0]:
-        if i.when < 100:
+        if i.when < 101:
             assert i.size == 100, f"{i}"
         else:
             assert i.size == 200, f"{i}"
@@ -936,3 +943,70 @@ def test_evolve_demes_model_starting_with_two_pops_and_no_ancestry(
         else:
             raise RuntimeError("unexpected key")
 
+
+def test_four_deme_split_model():
+    """
+    Tests primarily for lack of temporal
+    overlap b/w ancestral/derived demes.
+    See GitHub issue #814
+    """
+    yaml = """time_units: generations
+demes:
+- name: A
+  epochs:
+    - {end_time: 45, start_size: 10}
+- name: B
+  ancestors: [A] 
+  epochs:
+    - {end_time: 30, start_size: 20}
+- name: C
+  ancestors: [B] 
+  epochs:
+    - {end_time: 15, start_size: 30}
+- name: D
+  ancestors: [C] 
+  epochs:
+    - {end_time: 0, start_size: 40}
+"""
+    burnin = 1
+    g = demes.loads(yaml)
+    model = fwdpy11.discrete_demography.from_demes(g, burnin=burnin)
+
+    @dataclass
+    class DemeSizeAtTime:
+        when: int
+        size: int
+
+    class DemeSizes(object):
+        def __init__(self):
+            self.sizes = dict()
+
+        def __call__(self, pop, _):
+            deme_sizes = pop.deme_sizes()
+            assert len(deme_sizes[0]) == 1
+            for key, value in pop.deme_sizes(as_dict=True).items():
+                if key not in self.sizes:
+                    self.sizes[key] = [DemeSizeAtTime(when=pop.generation, size=value)]
+                else:
+                    self.sizes[key].append(
+                        DemeSizeAtTime(when=pop.generation, size=value)
+                    )
+
+    pdict = {
+        "gvalue": fwdpy11.Multiplicative(2.0),
+        "rates": (0, 0, 0),
+        "demography": model,
+        "simlen": model.metadata["total_simulation_length"],
+    }
+    params = fwdpy11.ModelParams(**pdict)
+    initial_size = g["A"].epochs[0].start_size
+    pop = fwdpy11.DiploidPopulation(initial_size, 1.0)
+    rng = fwdpy11.GSLrng(90210)
+    recorder = DemeSizes()
+    fwdpy11.evolvets(rng, pop, params, 100, recorder=recorder)
+
+    for deme in recorder.sizes:
+        if deme > 0:
+            assert len(recorder.sizes[deme]) == 15
+        else:
+            assert len(recorder.sizes[deme]) == burnin * initial_size