Improve early detection of invalid optimum values. (#1299)

See #1217 for context.

* Add new custom exception class
* Detect setup errors when constructing ModelParams
* Update tests and manual pages.

This PR breaks Python API for "gentic value" objects.
Such object must now provide a "validation" function
that is run by ModelParams

Note: deprecated classes are not updated to be compatible
with the new validation API.

Author: molpopgen

doc/pages/mvdes.md

@@ -60,8 +60,9 @@ migrations:
 """
 
 graph = demes.loads(yaml)
-demography = fwdpy11.discrete_demography.from_demes(graph, 1)
+demography = fwdpy11.ForwardDemesGraph.from_demes(graph, 1)
 
+optima = [fwdpy11.Optimum(when=0, optimum=0.0, VS=10.0)]
 pdict = {
     "nregions": [],
     "recregions": [],
@@ -74,7 +75,8 @@ pdict = {
     "demography": demography,
     "simlen": 100,
     "gvalue": fwdpy11.Additive(
-        ndemes=2, scaling=2, gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=10.0)
+        ndemes=2, scaling=2,
+        gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait(optima)
     ),
     "prune_selected": False,
 }

doc/short_vignettes/multivariate_gaussian_effect_sizes_across_demes.md

@@ -59,7 +59,7 @@ migrations:
     rate: 0.10
 """
 g = demes.loads(yaml)
-model = fwdpy11.discrete_demography.from_demes(g, burnin=1)
+model = fwdpy11.ForwardDemesGraph.from_demes(g, burnin=1)
 demesdraw.tubes(g);
 ```
 
@@ -76,7 +76,7 @@ pdict = {
     ],
     "rates": (0, 2.5e-3, None),
     "demography": model,
-    "simlen": model.metadata["total_simulation_length"],
+    "simlen": model.final_generation,
     "gvalue": fwdpy11.Additive(
         ndemes=3, scaling=2,
         gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=0.0, VS=10.0, when=0)])
@@ -105,7 +105,7 @@ Let's evolve the model now:
 params = fwdpy11.ModelParams(**pdict)
 # TODO: update this once we have a function to pull the sizes
 # automatically from demes-derived models:
-initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
+initial_sizes = model.initial_sizes
 pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
 rng = fwdpy11.GSLrng(1010)
 fwdpy11.evolvets(rng, pop, params, 10)
@@ -138,7 +138,7 @@ pdict["sregions"] = [
 params = fwdpy11.ModelParams(**pdict)
 # TODO: update this once we have a function to pull the sizes
 # automatically from demes-derived models:
-initial_sizes = [v for v in model.metadata["initial_sizes"].values()]
+initial_sizes = model.initial_sizes
 pop = fwdpy11.DiploidPopulation(initial_sizes, 1.0)
 fwdpy11.evolvets(rng, pop, params, 10)
 for i in pop.tables.mutations:

doc/short_vignettes/mutationdominance_vignette.md

@@ -48,18 +48,18 @@ pop = fwdpy11.DiploidPopulation(500, 1.0)
 
 rng = fwdpy11.GSLrng(54321)
 
-GSSmo = fwdpy11.GSSmo(
-    [
-        fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
-        fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),
-    ]
-)
+optima = [
+    fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
+    fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),
+]
+
+GSS = fwdpy11.GaussianStabilizingSelection.single_trait(optima)
 
 rho = 1000.
 
 p = {
     "nregions": [],
-    "gvalue": fwdpy11.Additive(2.0, GSSmo),
+    "gvalue": fwdpy11.Additive(2.0, GSS),
     "sregions": [des],
     "recregions": [fwdpy11.PoissonInterval(0, 1., rho / float(4 * pop.N))],
     "rates": (0.0, 1e-3, None),

doc/short_vignettes/poptour_vignette.md

@@ -30,18 +30,18 @@ pop = fwdpy11.DiploidPopulation(500, 1.0)
 
 rng = fwdpy11.GSLrng(54321)
 
-GSSmo = fwdpy11.GSSmo(
-    [
-        fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
-        fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),
-    ]
-)
+optima = [
+    fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
+    fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),
+]
+
+GSS = fwdpy11.GaussianStabilizingSelection.single_trait(optima)
 
 rho = 1000.
 
 p = {
     "nregions": [],
-    "gvalue": fwdpy11.Additive(2.0, GSSmo),
+    "gvalue": fwdpy11.Additive(2.0, GSS),
     "sregions": [fwdpy11.GaussianS(0, 1., 1, 0.1)],
     "recregions": [fwdpy11.PoissonInterval(0, 1., rho / float(4 * pop.N))],
     "rates": (0.0, 1e-3, None),

doc/short_vignettes/recorders_vignette.md

@@ -91,12 +91,12 @@ pop = fwdpy11.DiploidPopulation(500, 1.0)
 
 rng = fwdpy11.GSLrng(54321)
 
-GSSmo = fwdpy11.GSSmo(
-    [
-        fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
-        fwdpy11.Optimum(when=pop.N - 200, optimum=1.0, VS=1.0),
-    ]
-)
+optima= [
+    fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
+    fwdpy11.Optimum(when=pop.N - 200, optimum=1.0, VS=1.0),
+]
+
+GSS = fwdpy11.GaussianStabilizingSelection.single_trait(optima)
 
 rho = 1000.
 
@@ -105,7 +105,7 @@ des = fwdpy11.GaussianS(beg=0, end=1, weight=1, sd=0.1,
 
 p = {
     "nregions": [],
-    "gvalue": fwdpy11.Additive(2.0, GSSmo),
+    "gvalue": fwdpy11.Additive(2.0, GSS),
     "sregions": [des],
     "recregions": [fwdpy11.PoissonInterval(0, 1., rho / float(4 * pop.N))],
     "rates": (0.0, 1e-3, None),
@@ -169,7 +169,7 @@ Thus, we need to rebuild the `gvalue` field.
 import pickle
 pop = fwdpy11.DiploidPopulation(500, 1.0)
 rng = fwdpy11.GSLrng(54321)
-p['gvalue'] = fwdpy11.Additive(2.0, GSSmo)
+p['gvalue'] = fwdpy11.Additive(2.0, GSS)
 params = fwdpy11.ModelParams(**p)
 recorder = RandomSamples(5)
 fwdpy11.evolvets(rng, pop, params, recorder=recorder, simplification_interval=100) 

doc/short_vignettes/workingexample_trait.md

@@ -37,7 +37,7 @@ p = {
     "rates": (mu_neutral, mu_selected, None),
     # Keep mutations at frequency 1 in the pop if they affect fitness.
     "prune_selected": False,
-    "demography": None,
+    "demography": fwdpy11.ForwardDemesGraph.tubes([N], burnin=1),
     "simlen": 10 * N,
 }
 params = fwdpy11.ModelParams(**p)

examples/plugin/test_plugin.py

@@ -13,7 +13,7 @@
 """
 
 graph = demes.loads(yaml)
-demography = fwdpy11.discrete_demography.from_demes(graph, 1)
+demography = fwdpy11.ForwardDemesGraph.from_demes(graph, 1)
 
 pdict = {
     "demography": demography,

fwdpy11/__init__.py

@@ -61,6 +61,7 @@
     StrictAdditiveMultivariateEffects,
     AdditivePleiotropy,
     PyDiploidGeneticValue,
+    TimingError,
 )
 
 from ._types import (  # NOQA

fwdpy11/_types/model_params.py

@@ -17,10 +17,12 @@
 # along with fwdpy11.  If not, see <http://www.gnu.org/licenses/>.
 #
 
+import abc
 import typing
 import warnings
 
 import attr
+import demes
 import fwdpy11
 import numpy as np
 from fwdpy11.class_decorators import (
@@ -96,6 +98,20 @@ def _convert_rates(value):
         return MutationAndRecombinationRates(**value)
 
 
+@typing.runtime_checkable
+class ValidateEventTimings(typing.Protocol):
+    @abc.abstractmethod
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        raise NotImplementedError
+
+
+def validate_timings(
+    gvalue: ValidateEventTimings, deme: int, demography: ForwardDemesGraph
+) -> None:
+    assert isinstance(gvalue, ValidateEventTimings)
+    return gvalue.validate_timings(deme, demography)
+
+
 @attr_add_asblack
 @attr_class_to_from_dict_no_recurse
 @attr.s(kw_only=True, frozen=True, slots=True, repr_ns="fwdpy11")
@@ -299,15 +315,17 @@ def rates_validator(self, attribute, value):
     @gvalue.validator
     def validate_gvalue(self, attribute, value):
         try:
-            for i in value:
+            for deme, i in enumerate(value):
                 attr.validators.instance_of(fwdpy11.DiploidGeneticValue)(
                     self, attribute, i
                 )
+                validate_timings(i, deme, self.demography)
         except TypeError:
             try:
                 attr.validators.instance_of(fwdpy11.DiploidGeneticValue)(
                     self, attribute, value
                 )
+                validate_timings(value, 0, self.demography)
             except TypeError:
                 raise
 

fwdpy11/genetic_values.py

@@ -22,10 +22,12 @@
 import warnings
 
 import attr
+import demes
 import numpy as np
 
 from deprecated import deprecated
 
+
 from ._fwdpy11 import (
     GeneticValueIsTrait,
     GeneticValueNoise,
@@ -49,6 +51,12 @@
     attr_class_to_from_dict_no_recurse,
 )
 
+from ._types.forward_demes_graph import ForwardDemesGraph
+
+
+class TimingError(Exception):
+    pass
+
 
 @attr_add_asblack
 @attr_class_pickle_with_super
@@ -204,6 +212,20 @@ def asdict(self):
     def fromdict(cls, data):
         return cls(**data)
 
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        graph = demography.demes_graph
+        deme_obj = graph.demes[deme]
+        if deme_obj.start_time == float("inf"):
+            start_time = demography.to_backwards_time(0)
+        else:
+            start_time = deme_obj.start_time
+        for i in self.optima[:1]:
+            btime = demography.to_backwards_time(i.when)
+            if btime < start_time:
+                msg = f"deme {deme_obj.name}, event {i} is at time {btime} in the past "
+                msg += f", but the deme starts at {int(start_time)}"
+                raise TimingError(msg)
+
 
 @attr_add_asblack
 @attr_class_to_from_dict_no_recurse
@@ -508,6 +530,11 @@ def __attrs_post_init__(self):
             self.scaling, self.gvalue_to_fitness, self.noise, self.ndemes
         )
 
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        if self.gvalue_to_fitness is None:
+            return
+        self.gvalue_to_fitness.validate_timings(deme, demography)
+
 
 @attr_add_asblack
 @attr_class_pickle_with_super
@@ -555,6 +582,11 @@ def __attrs_post_init__(self):
             self.scaling, self.gvalue_to_fitness, self.noise, self.ndemes
         )
 
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        if self.gvalue_to_fitness is None:
+            return
+        self.gvalue_to_fitness.validate_timings(deme, demography)
+
 
 @attr_add_asblack
 @attr_class_pickle_with_super
@@ -591,6 +623,11 @@ class GBR(_ll_GBR):
     def __attrs_post_init__(self):
         super(GBR, self).__init__(self.gvalue_to_fitness, self.noise)
 
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        if self.gvalue_to_fitness is None:
+            return
+        self.gvalue_to_fitness.validate_timings(deme, demography)
+
 
 @attr_add_asblack
 @attr_class_pickle_with_super
@@ -637,6 +674,11 @@ def __attrs_post_init__(self):
             self.ndimensions, self.focal_trait, self.gvalue_to_fitness, self.noise
         )
 
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        if self.gvalue_to_fitness is None:
+            return
+        self.gvalue_to_fitness.validate_timings(deme, demography)
+
 
 @deprecated(reason="Use AdditivePleiotropy instead.")
 class StrictAdditiveMultivariateEffects(AdditivePleiotropy):
@@ -652,3 +694,8 @@ def __init__(self, ndim: int, gvalue_to_fitness=None, noise=None):
             gvalue_to_fitness,
             noise,
         )
+
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        if self.gvalue_to_fitness is None:
+            return
+        self.gvalue_to_fitness.validate_timings(deme, demography)

tests/pygss.py

@@ -24,6 +24,8 @@
 import fwdpy11
 import fwdpy11.custom_genetic_value_decorators
 
+from fwdpy11._types import ForwardDemesGraph
+
 
 @fwdpy11.custom_genetic_value_decorators.default_update
 @fwdpy11.custom_genetic_value_decorators.genetic_value_is_trait_default_clone()
@@ -41,6 +43,9 @@ def __call__(self, data: fwdpy11.DiploidGeneticValueToFitnessData) -> float:
             / (2 * self.VS)
         )
 
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        pass
+
 
 @fwdpy11.custom_genetic_value_decorators.genetic_value_is_trait_default_clone()
 @attr.s()
@@ -65,3 +70,6 @@ def update(self, pop: fwdpy11.DiploidPopulation) -> None:
             1,
         )[0]
         self.optima.append((pop.generation, self.opt))
+
+    def validate_timings(self, deme: int, demography: ForwardDemesGraph) -> None:
+        pass