Add tests of evolving demes models in multiple stages

cpptests/discrete_demography_roundtrips.cc

@@ -54,15 +54,20 @@ DiscreteDemography_roundtrip(
 // 1. Initial deme labels are set by the user. NOTE: validated by manager object
 //
 {
+    bool demographic_model_needs_updating = true;
     auto current_demographic_state
         = fwdpy11::discrete_demography::initialize_model_state(
-            pop.generation, pop.diploid_metadata, demography);
+            pop.generation, pop.diploid_metadata, demography,
+            &demographic_model_needs_updating);
     decltype(pop.diploid_metadata) offspring_metadata;
     offspring_metadata.reserve(pop.N);
     std::vector<MatingEventRecord> rv;
-    fwdpy11::discrete_demography::update_demography_manager(
-        rng, pop.generation, pop.diploid_metadata, demography,
-        *current_demographic_state);
+    if (demographic_model_needs_updating)
+        {
+            fwdpy11::discrete_demography::update_demography_manager(
+                rng, pop.generation, pop.diploid_metadata, demography,
+                *current_demographic_state);
+        }
     if (current_demographic_state->will_go_globally_extinct() == true)
         {
             std::ostringstream o;

doc/misc/changelog.md

@@ -11,6 +11,18 @@ Bug fixes
   This does not affect previous results because unsigned overflow doing the "right thing" ended up with final values being correct.
   {pr}`766`
   {issue}`765`
+* Fix a bug where stopping/restarting the evolution of demographic models at time points
+  where a deme goes extinct.
+  It is not possible that this bug affected results from earlier versions, as attempting to stop/start at these time points raised exceptions.
+  {issue}`775`
+  {pr}`774`
+
+Behavior changes
+
+* If a demographic model is evolved, pickled, unpickled, and then used to evolve,
+  it is now possible that exceptions will raise.
+  This change is due to the fix for {issue}`775` introduced in {pr}`774`.
+  See issue {issue}`777` for more background.
 
 New features
 

fwdpy11/headers/fwdpy11/discrete_demography/simulation/demographic_model_state.hpp

@@ -103,15 +103,19 @@ namespace fwdpy11
             }
         };
 
+        /// model_needs_update added in 0.16.0
+        /// in relation to GitHub issue 775
         template <typename METADATATYPE>
         inline demographic_model_state_pointer
         initialize_model_state(std::uint32_t generation,
                                const std::vector<METADATATYPE>& metadata,
-                               DiscreteDemography& demography)
+                               DiscreteDemography& demography, bool* model_needs_update)
         {
             // "steal" pointer from input
             auto rv = demography.get_model_state();
 
+            *model_needs_update = (rv == nullptr);
+
             if (rv == nullptr || generation == 0)
                 // If there is no state, then we need to make
                 // one.  If there is a state, but the generation
@@ -121,6 +125,7 @@ namespace fwdpy11
                 {
                     demography.update_event_times(generation);
                     rv.reset(new demographic_model_state(metadata, demography));
+                    *model_needs_update = true;
                 }
             return rv;
         }

fwdpy11/headers/fwdpy11/discrete_demography/simulation/pick_parents.hpp

@@ -21,6 +21,7 @@
 #define FWDPY11_DISCRETE_DEMOGRAPY_PICK_PARENTS_HPP
 
 #include <cstdint>
+#include <sstream>
 #include "../../rng.hpp"
 #include "../MigrationMatrix.hpp"
 #include "migration_lookup.hpp"
@@ -70,7 +71,9 @@ namespace fwdpy11
                 }
             if (miglookup.lookups[offspring_deme] == nullptr)
                 {
-                    throw DemographyError("parental deme lookup is NULL");
+                    std::ostringstream o;
+                    o <<"parental deme lookup is NULL for deme " << offspring_deme; 
+                    throw DemographyError(o.str());
                 }
             std::int32_t pdeme = static_cast<std::int32_t>(gsl_ran_discrete(
                 rng.get(), miglookup.lookups[offspring_deme].get()));

fwdpy11/src/evolve_population/evolvets.cc

@@ -222,8 +222,12 @@ evolve_with_tree_sequences(
         }
 
     // Set up discrete demography types. New in 0.6.0
-    auto current_demographic_state = ddemog::initialize_model_state(
-        pop.generation, pop.diploid_metadata, demography);
+    /// model_needs_update added in 0.16.0
+    /// in relation to GitHub issue 775
+    bool demographic_model_needs_update = true;
+    auto current_demographic_state
+        = ddemog::initialize_model_state(pop.generation, pop.diploid_metadata,
+                                         demography, &demographic_model_needs_update);
     if (gvalue_pointers.genetic_values.size()
         > static_cast<std::size_t>(current_demographic_state->maxdemes))
         {
@@ -285,8 +289,11 @@ evolve_with_tree_sequences(
                               record_genotype_matrix);
     pop.genetic_value_matrix.swap(new_diploid_gvalues);
     pop.diploid_metadata.swap(offspring_metadata);
-    ddemog::update_demography_manager(rng, pop.generation, pop.diploid_metadata,
-                                      demography, *current_demographic_state);
+    if (demographic_model_needs_update == true)
+        {
+            ddemog::update_demography_manager(rng, pop.generation, pop.diploid_metadata,
+                                              demography, *current_demographic_state);
+        }
     if (current_demographic_state->will_go_globally_extinct() == true)
         {
             std::ostringstream o;

tests/conftest.py

@@ -113,3 +113,53 @@ def small_split_model():
   ancestors: [ancestral]
 """
     return demes.loads(yaml)
+
+
+@pytest.fixture
+def two_deme_split_with_ancestral_size_change():
+    yaml = """
+description: Two deme model with migration and size changes.
+time_units: generations
+demes:
+- name: ancestral
+  description: ancestral deme, two epochs
+  epochs:
+  - {end_time: 20, start_size: 100}
+  - {end_time: 10, start_size: 200}
+- name: deme1
+  description: child 1
+  epochs:
+  - {start_size: 250, end_size: 500, end_time: 0}
+  ancestors: [ancestral]
+- name: deme2
+  description: child 2
+  epochs:
+  - {start_size: 50, end_size: 200, end_time: 0}
+  ancestors: [ancestral]
+migrations:
+- {demes: [deme1, deme2], rate: 1e-3}
+"""
+    return demes.loads(yaml)
+
+
+@pytest.fixture
+def start_demes_model_with_two_pops():
+    yaml = """
+description: after slicing original graph at time 50 (below slice)
+time_units: generations
+demes:
+- name: deme1
+  epochs:
+  - {start_size: 100, end_time: 50}
+  - {start_size: 100, end_size: 200, end_time: 0}
+- name: deme2
+  epochs:
+  - {start_size: 75, end_time: 50}
+  - {start_size: 75, end_time: 20}
+  - {start_size: 300, end_time: 0}
+migrations:
+- demes: [deme1, deme2]
+  rate: 1e-2
+  start_time: 50
+"""
+    return demes.loads(yaml)

tests/test_demes2fwdpy11.py

@@ -1,8 +1,10 @@
+import unittest
+from dataclasses import dataclass
+
 import demes
+import fwdpy11
 import numpy as np
-import unittest
 import pytest
-import fwdpy11
 
 
 def check_valid_demography(cls):
@@ -521,3 +523,135 @@ def test_pulse_migration_matrix(self):
         self.assertTrue(
             np.all(self.demog.model.set_migration_rates[1].migrates == [0, 1])
         )
+
+
+def test_split_model_population_size_history(two_deme_split_with_ancestral_size_change):
+    """
+    This is a detailed test of the complete size history
+    of a model that we run all the way through.
+    """
+    model = fwdpy11.discrete_demography.from_demes(
+        two_deme_split_with_ancestral_size_change, burnin=1
+    )
+
+    @dataclass
+    class DemeSizeAtTime:
+        when: int
+        size: int
+
+    class DemeSizes(object):
+        def __init__(self):
+            self.sizes = dict()
+
+        def __call__(self, pop, _):
+            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)
+    pop = fwdpy11.DiploidPopulation(100, 1.0)
+    rng = fwdpy11.GSLrng(90210)
+    recorder = DemeSizes()
+    fwdpy11.evolvets(rng, pop, params, 100, recorder=recorder)
+
+    # The ancestral deme exists until generation 110,
+    # and we only see offspring from birth time 1 on.
+    assert [i.when for i in recorder.sizes[0]] == [i for i in range(1, 111)]
+    # 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)
+        ]
+    # initial daughter deme sizes
+    assert recorder.sizes[1][0].size == 250
+    assert recorder.sizes[2][0].size == 50
+    # 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:
+            assert i.size == 100, f"{i}"
+        else:
+            assert i.size == 200, f"{i}"
+
+
+@pytest.mark.parametrize("when", [i for i in range(75, 120)])
+def test_evolve_population_in_two_stages(
+    when, two_deme_split_with_ancestral_size_change
+):
+    model = fwdpy11.discrete_demography.from_demes(
+        two_deme_split_with_ancestral_size_change, burnin=1
+    )
+    pdict = {
+        "gvalue": fwdpy11.Multiplicative(2.0),
+        "rates": (0, 0, 0),
+        "demography": model,
+        "simlen": when,
+    }
+    params = fwdpy11.ModelParams(**pdict)
+    pop = fwdpy11.DiploidPopulation(100, 1.0)
+    rng = fwdpy11.GSLrng(90210)
+    fwdpy11.evolvets(rng, pop, params, 100)
+
+    pdict["simlen"] = model.metadata["total_simulation_length"] - when
+    params = fwdpy11.ModelParams(**pdict)
+
+    fwdpy11.evolvets(rng, pop, params, 100, check_demographic_event_timings=False)
+
+    counts = np.unique(np.array(pop.diploid_metadata)["deme"], return_counts=True)
+    assert counts[1][0] == 500
+    assert counts[1][1] == 200
+
+
+# NOTE: update this test to have burnin=0 once GittHub issue 776
+# is fixed
+# NOTE: models like this are important to applications where
+# the ancestor is simulated with msprime, and the split part
+# happens in fwdpy11
+@pytest.mark.parametrize("burnin", [1])
+def test_evolve_demes_model_starting_with_two_pops_and_no_ancestry(
+    burnin,
+    start_demes_model_with_two_pops,
+):
+    model = fwdpy11.discrete_demography.from_demes(
+        start_demes_model_with_two_pops, burnin=burnin
+    )
+    pdict = {
+        "gvalue": fwdpy11.Multiplicative(2.0),
+        "rates": (0, 0, 0),
+        "demography": model,
+        "simlen": model.metadata["total_simulation_length"],
+    }
+    params = fwdpy11.ModelParams(**pdict)
+    print(model.metadata["initial_sizes"])
+    initial_sizes = [i for i in model.metadata["initial_sizes"].values()]
+    pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
+    for key, value in pop.deme_sizes(as_dict=True).items():
+        if key == 0:
+            assert value == 100
+        elif key == 1:
+            assert value == 75
+        else:
+            raise RuntimeError("unexpected key")
+    rng = fwdpy11.GSLrng(101)
+    fwdpy11.evolvets(rng, pop, params, 100)
+    for key, value in pop.deme_sizes(as_dict=True).items():
+        if key == 0:
+            assert value == 200
+        elif key == 1:
+            assert value == 300
+        else:
+            raise RuntimeError("unexpected key")

tests/test_discrete_demography_with_tree_sequences.py

@@ -20,9 +20,8 @@
 import unittest
 from collections import namedtuple
 
-import numpy as np
-
 import fwdpy11
+import numpy as np
 
 
 def validate_alive_node_metadata(pop):
@@ -59,8 +58,7 @@ class TestConstantPopulationSize(unittest.TestCase):
 
     @classmethod
     def setUp(self):
-        from test_demographic_models import setup_pop_rng
-        from test_demographic_models import setup_pdict
+        from test_demographic_models import setup_pdict, setup_pop_rng
 
         d = fwdpy11.DiscreteDemography()
         self.pop, self.rng = setup_pop_rng()
@@ -810,6 +808,9 @@ def test_evolve_in_two_steps_restart_with_two_demes(self):
         self.assertEqual(deme_sizes[1][0], self.N0t)
         self.assertEqual(deme_sizes[1][1], self.N1t)
 
+    @unittest.skip(
+        reason="Fixing issue 775 means that this no longer works because pickling DiscreteDemography does not fully round-trip the internal model state"
+    )
     def test_evolve_in_two_steps_restart_with_two_demes_and_pickle(self):
         """
         Tests the more complex case of restarting a sim