Simplify gaussian stabilizing selection Python API.

* Deprecate the current 4 classes. * Add fwdpy11.GaussianStabilizingSelection to replace * Update tests * Update manual Closes #463
molpopgen · Aug 19, 2023 · bced910 · bced910
1 parent 401ed3b
commit bced910
Show file tree

Hide file tree

Showing 30 changed files with 268 additions and 111 deletions.
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
@@ -54,6 +54,7 @@ set(GENETIC_VALUE_TO_FITNESS_SOURCES
     genetic_value_to_fitness/GeneticValueIsFitness.cc
     genetic_value_to_fitness/GeneticValueIsTrait.cc
     genetic_value_to_fitness/GSSmo.cc
+    genetic_value_to_fitness/GaussianStabilizingSelection.cc
     genetic_value_to_fitness/MultivariateGSSmo.cc
     genetic_value_to_fitness/Optimum.cc
     genetic_value_to_fitness/PleiotropicOptima.cc)

diff --git a/cpp/genetic_value_to_fitness/GaussianStabilizingSelection.cc b/cpp/genetic_value_to_fitness/GaussianStabilizingSelection.cc
@@ -0,0 +1,19 @@
+#include <tuple>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <fwdpy11/genetic_value_to_fitness/GaussianStabilizingSelection.hpp>
+
+namespace py = pybind11;
+
+void
+init_GaussianStabilizingSelection(py::module& m)
+{
+    py::class_<fwdpy11::GaussianStabilizingSelection, fwdpy11::GeneticValueIsTrait>(
+        m, "_ll_GaussianStabilizingSelection")
+        .def(py::init([](const fwdpy11::GSSmo& single_trait) {
+            return fwdpy11::GaussianStabilizingSelection(single_trait);
+        }))
+        .def(py::init([](const fwdpy11::MultivariateGSSmo& pleiotropy) {
+            return fwdpy11::GaussianStabilizingSelection(pleiotropy);
+        }));
+}
diff --git a/cpp/genetic_value_to_fitness/init.cc b/cpp/genetic_value_to_fitness/init.cc
@@ -11,6 +11,7 @@ void init_GeneticValueToFitnessMap(py::module&);
 void init_GeneticValueIsTrait(py::module&);
 void init_GeneticValueIsFitness(py::module&);
 void init_GSSmo(py::module&);
+void init_GaussianStabilizingSelection(py::module&);
 
 // Multivariate classes
 void init_MultivariateGSSmo(py::module&);
@@ -25,6 +26,7 @@ initialize_genetic_value_to_fitness(py::module& m)
     init_GeneticValueIsTrait(m);
     init_GeneticValueIsFitness(m);
     init_GSSmo(m);
+    init_GaussianStabilizingSelection(m);
 
     init_MultivariateGSSmo(m);
 }
diff --git a/doc/long_vignettes/gss_vignette.md b/doc/long_vignettes/gss_vignette.md
@@ -38,7 +38,7 @@ class Recorder(object):
 ```
 
 Now, we will set up and simulate the model.
-We use {class}`fwdpy11.GSSmo` to specify when/how the optimum value shifts.
+We use {class}`fwdpy11.GaussianStabilizingSelection` to specify when/how the optimum value shifts.
 
 We will simulate the population for $10N$ generations around an optimum of zero.
 Then, we shift the optimum to 1 and evolve another 200 generations.
@@ -51,7 +51,7 @@ pop = fwdpy11.DiploidPopulation(500, 1.0)
 
 rng = fwdpy11.GSLrng(54321)
 
-GSSmo = fwdpy11.GSSmo(
+gssmo = fwdpy11.GaussianStabilizingSelection.single_trait(
     [
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),
@@ -62,7 +62,7 @@ rho = 1000.
 
 p = {
     "nregions": [],
-    "gvalue": fwdpy11.Additive(2.0, GSSmo),
+    "gvalue": fwdpy11.Additive(2.0, gssmo),
     "sregions": [fwdpy11.GaussianS(0, 1., 1, 0.1)],
     "recregions": [fwdpy11.PoissonInterval(0, 1., rho / float(4 * pop.N))],
     "rates": (0.0, 1e-3, None),

diff --git a/doc/long_vignettes/tskitconvert_vignette.md b/doc/long_vignettes/tskitconvert_vignette.md
@@ -138,11 +138,13 @@ assert "model_params" not in ts.metadata
 Here is an example of Gaussian stabilizing selection around an optimum trait value of `0.0`:
 
 ```{code-cell} python
+optimum = fwdpy11.Optimum(optimum=0.0, VS=10.0, when=0)
+gss = fwdpy11.GaussianStabilizingSelection.single_trait([optimum])
 pdict = {
     'nregions': [],
     'sregions': [fwdpy11.GaussianS(0, 1, 1, 0.10)],
     'recregions': [fwdpy11.PoissonInterval(0, 1, 1e-3)],
-    'gvalue': fwdpy11.Additive(2., fwdpy11.GSS(optimum=0.0, VS=10.0)),
+    'gvalue': fwdpy11.Additive(2., gss),
     'rates': (0.0, 1e-3, None),
     'simlen': 10 * pop.N,
     'prune_selected': False,
@@ -165,16 +167,18 @@ Imagine, for example, that we simulate the above model until equilibrium, and th
 The details are omitted here:
 
 1. This simulation requires two calls to {func}`fwdpy11.evolvets`
-2. It could have been done instead with a single call to {func}`fwdpy11.evolvets` and using
-   {class}`fwdpy11.GSSmo` and one `ModelParams` object instead of {class}`fwdpy11.GSS`
-   and two parameters objects.
+2. It could have been done instead with a single call to {func}`fwdpy11.evolvets` 
+   by passing a list of {class}`fwdpy11.Optimum` to {class}`fwdpy11.GaussianStabilizingSelection`.
+   This approach would have required one `ModelParams` instance.
 :::
 
 ```{code-cell} python
 import copy
 pdict2 = copy.deepcopy(pdict)
 
-pdict2['gvalue'] = fwdpy11.Additive(2., fwdpy11.GSS(optimum=5.0, VS=10.0)),
+optimum = fwdpy11.Optimum(optimum=5.0, VS=10.0, when=0)
+gss = fwdpy11.GaussianStabilizingSelection.single_trait([optimum])
+pdict2['gvalue'] = fwdpy11.Additive(2., gss),
 pdict2['simlen'] = 100
 
 params2 = fwdpy11.ModelParams(**pdict2)

diff --git a/doc/misc/changelog.md b/doc/misc/changelog.md
@@ -3,6 +3,22 @@
 Major changes are listed below.  Each release likely contains fiddling with back-end code,
 updates to latest `fwdpp` version, etc.
 
+## Next release
+
+Deprecations
+
+The following are deprecated:
+
+* `fwdpy11.GSS`
+* `fwdpy11.GSSmo`
+* `fwdpy11.MutivariateGSS`
+* `fwdpy11.MutivariateGSSmo`
+
+Their functionality is replaced with {class}`fwdpy11.GaussianStabilizingSelection`.
+
+PR {pr}`1166`.
+Issue {issue}{`463`}.
+
 ## 0.21.0
 
 Breaking changes
@@ -643,12 +659,12 @@ Changes to `fwdpy11.conditional_models`:
 
 Bug fixes
 
-* Fixed bug in updating {class}`fwdpy11.MultivariateGSSmo`.
+* Fixed bug in updating `fwdpy11.MultivariateGSSmo`.
   PR {pr}`867`.
 
 Back end changes:
 
-* {class}`fwdpy11.GSSmo` and {class}`fwdpy11.MultivariateGSSmo` now handle cases where population start time is greater than zero.
+* `fwdpy11.GSSmo` and `fwdpy11.MultivariateGSSmo` now handle cases where population start time is greater than zero.
   PR {pr}`867`.
 
 Build system:
@@ -1260,7 +1276,7 @@ Maintenance release and one new feature:
 
 * Allow the first generation of a simulation to be preserved. PR {pr}`470`
   See {ref}`recapitation`.
-* Parameterizing classes like {class}`fwdpy11.GSSmo` is now more Pythonic,
+* Parameterizing classes like `fwdpy11.GSSmo` is now more Pythonic,
   and some existing `init` methods are deprecated in favor of the
   new approach. PR {pr}`461`.
 
@@ -1479,7 +1495,7 @@ Miscellaneous changes:
 The following bugs are fixed:
 
 * Mutations were not being recycled properly during simulations with tree sequences, resulting in excessive memory consumption. PR {pr}`317`
-* Several interface issues with {class}`fwdpy11.MultivariateGSSmo` are fixed. PR {pr}`313`
+* Several interface issues with `fwdpy11.MultivariateGSSmo` are fixed. PR {pr}`313`
 * Fix a bug that could lead to fixations with tree sequences not "pruning" selected fixations when that behavior is desired. {issue}`287`, fixed in PR {pr}`289`
 * A memory safety issue was fixed in the implementation of {attr}`fwdpy11.TreeIterator.samples_below`. PR {pr}`300`.  {issue}`299`
 

diff --git a/doc/pages/advancedtopics.md b/doc/pages/advancedtopics.md
@@ -560,7 +560,8 @@ class MyGeneticValueToFitness(fwdpy11.GeneticValueIsTrait):
 ```
 
 Two complete examples come from the test suite.  The first example reimplements
-{class}`fwdpy11.GSS` in Python.  The second example implements a model
+{class}`fwdpy11.GaussianStabilizingSelection` in Python for a single trait.
+The second example implements a model
 where the optimum changes to a new value each generation. (Note that
 the second model requires that {func}`numpy.random.seed` be called elsewhere
 so that results are reproducible.)

diff --git a/doc/pages/gvalue_to_fitness.md b/doc/pages/gvalue_to_fitness.md
@@ -13,16 +13,6 @@
 .. autoclass:: fwdpy11.GeneticValueIsFitness
 ```
 
-```{eval-rst}
-.. autoclass:: fwdpy11.GSS
-    :members: asdict, fromdict, asblack
-```
-
-```{eval-rst}
-.. autoclass:: fwdpy11.GSSmo
-    :members: asdict, fromdict, asblack
-```
-
 ```{eval-rst}
 .. autoclass:: fwdpy11.Optimum
     :members: asdict, fromdict, asblack
@@ -33,14 +23,10 @@
     :members: asdict, fromdict, asblack
 ```
 
-```{eval-rst}
-.. autoclass:: fwdpy11.MultivariateGSS
-    :members: asdict, fromdict, asblack
-```
 
 ```{eval-rst}
-.. autoclass:: fwdpy11.MultivariateGSSmo
-    :members: asdict, fromdict, asblack
+.. autoclass:: fwdpy11.GaussianStabilizingSelection
+    :members:
 ```
 
 
diff --git a/doc/short_vignettes/evolvets_vignette.md b/doc/short_vignettes/evolvets_vignette.md
@@ -35,7 +35,7 @@ pop = fwdpy11.DiploidPopulation(500, 1.0)
 
 rng = fwdpy11.GSLrng(54321)
 
-GSSmo = fwdpy11.GSSmo(
+GSSmo = fwdpy11.GaussianStabilizingSelection.single_trait(
     [
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),

diff --git a/doc/short_vignettes/gvalues_traits.md b/doc/short_vignettes/gvalues_traits.md
@@ -43,23 +43,24 @@ To specify an additive effects model of a trait under Gaussian stabilizing selec
 ```{code-cell} python
 import fwdpy11
 
+gss = fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=0.0, VS=1.0, when=0)])
 gvalue = fwdpy11.Additive(
-    scaling=2.0, gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=1.0)
+    scaling=2.0, gvalue_to_fitness=gss,
 )
 ```
 
 Here, we are using a second parameter to initialize a "genetic value to fitness" map stored in an instance of {class}`fwdpy11.Additive`.
 ({class}`fwdpy11.Multiplicative` also supports such maps.)
-See {class}`fwdpy11.GSS` for details.
+See {class}`fwdpy11.GaussianStabilizingSelection` for details.
 
 We can also add Gaussian noise to an individual's trait value:
 
 ```{code-cell} python
 import numpy as np
-
+gss = fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=0.0, VS=2.0 / 3.0, when=0)])
 gvalue = fwdpy11.Additive(
     scaling=2.0,
-    gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=2.0 / 3.0),
+    gvalue_to_fitness=gss,
     noise=fwdpy11.GaussianNoise(mean=0.0, sd=np.sqrt(1.0 / 3.0)),
     )
 ```
@@ -69,24 +70,15 @@ The last example requires some explanation:
 * We want `VS = 1.0`.  We can decompose `VS = VW + VE`, where `VW` and `VE` are the additive contributions of genetic and environmental effects.
 * Here, the environmental effect is a Gaussian with mean zero and variance `1/3`.
   The class is parameterized with the standard deviation, however, so we need to pass on the square root.
-* We then set `VS = 1 - 1/3 = 2/3` when initializing {class}`fwdpy11.GSS`.
+* We then set `VS = 1 - 1/3 = 2/3` when initializing {class}`fwdpy11.Optimum`.
 
 Yes, this is a nomenclature issue!
-The `VS` argument to {class}`fwdpy11.GSS` really should be called `VW` and we'll fix that in a future version and hopefully not break people's code.
+The `VS` argument to {class}`fwdpy11.Optimum` really should be called `VW` and we'll fix that in a future version and hopefully not break people's code.
 
 In general, there's a good bit of subtlety to properly modeling quantitative traits.
 The machinery described here was used in {cite}`Thornton2019-nu`. {cite}`Burger2000-ul` is an excellent technical reference on the topic.
 {cite}`Walsh2018-ux` also thoroughly covers a lot of relevant material.
 
-:::{note}
-
-Under the hood, the `GSS` and `GSSmo` classes aren't that different.
-Their multivariate analogs are rather similar, too.
-Thus, we envision a future with one single `fwdpy11.GaussianStabilizingSelection` class to handle all cases.
-The types discussed here would remain as simple Python wrappers so that we don't break existing simulations.
-
-:::
-
 For an example of another approach to modeling phenotypes often attributed to {cite}`Eyre-Walker2010-rs`, see {ref}`here <eyre-walker>`.
 
 :::{todo}
@@ -98,11 +90,11 @@ Write (and refer to) an advanced section on pleiotropic models.
 ### A sudden optimum shift
 
 The previous example set up a model where the optimum is stable for the entire simulation.
-We can parameterize a shifting optimum using {class}`fwdpy11.GSSmo`.
+We can parameterize a shifting optimum using {class}`fwdpy11.GaussianStabilizingSelection`.
 For example, to shift the optimum from `0.0` to `1.0` at generation `100`:
 
 ```{code-cell} python
-moving_optimum = fwdpy11.GSSmo(
+moving_optimum = fwdpy11.GaussianStabilizingSelection.single_trait(
     [
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=100, optimum=1.0, VS=1.0),
@@ -137,7 +129,7 @@ while last_time < 1000:
     if last_time < 1000:
         optima.append(fwdpy11.Optimum(when=last_time, optimum=last_optimum, VS=10.0))
 
-random_moving_optimum = fwdpy11.GSSmo(optima)
+random_moving_optimum = fwdpy11.GaussianStabilizingSelection.single_trait(optima)
 random_moving_optimum
 ```
 
@@ -147,4 +139,4 @@ Note the cast to `int` when updating the time.
 {class}`fwdpy11.Optimum` is very picky about its input.
 It requires `int` for `when` and will raise an exception if the {attr}`numpy.int64` from {func}`numpy.random.geometric` gets passed in.
 
-:::
+S
diff --git a/doc/short_vignettes/multivariate_gaussian_effect_sizes_across_demes.md b/doc/short_vignettes/multivariate_gaussian_effect_sizes_across_demes.md
@@ -78,7 +78,8 @@ pdict = {
     "demography": model,
     "simlen": model.metadata["total_simulation_length"],
     "gvalue": fwdpy11.Additive(
-        ndemes=3, scaling=2, gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=10.0)
+        ndemes=3, scaling=2,
+        gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=0.0, VS=10.0, when=0)])
     ),
 }
 ```

diff --git a/doc/short_vignettes/neutralmutsafter_vignette.md b/doc/short_vignettes/neutralmutsafter_vignette.md
@@ -26,7 +26,7 @@ pop = fwdpy11.DiploidPopulation(500, 1.0)
 
 rng = fwdpy11.GSLrng(54321)
 
-GSSmo = fwdpy11.GSSmo(
+GSSmo = fwdpy11.GaussianStabilizingSelection.single_trait(
     [
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=10 * pop.N - 200, optimum=1.0, VS=1.0),

diff --git a/doc/short_vignettes/workingexample_trait.md b/doc/short_vignettes/workingexample_trait.md
@@ -22,7 +22,7 @@ rho = 1000.0
 mu_neutral = 0.0
 mu_selected = 1e-3
 
-GSSmo = fwdpy11.GSSmo(
+GSSmo = fwdpy11.GaussianStabilizingSelection.single_trait(
     [
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=10 * N - 200, optimum=1.0, VS=1.0),

diff --git a/doc/technical/genetic_values.md b/doc/technical/genetic_values.md
@@ -162,7 +162,7 @@ gv.noise_function.update(pop)
 The definition of our classes and the way that `update` functions are applied
 implies that all three objects are updated *independently from one another*.
 This independence covers a large number of use cases.  For example, a moving
-trait optimum ({class}`fwdpy11.GSSmo`) only needs acces to
+trait optimum ({class}`fwdpy11.GaussianStabilizingSelection`) only needs acces to
 {attr}`fwdpy11.DiploidPopulation.generation` to know if it needs to update its
 internal state.
 

diff --git a/fwdpy11/__init__.py b/fwdpy11/__init__.py
@@ -48,6 +48,7 @@
 from .genetic_values import (  # NOQA
     PleiotropicOptima,
     Optimum,
+    GaussianStabilizingSelection,
     GSS,
     GSSmo,
     MultivariateGSS,