From c5a2f29f496783b3dd37b2e0ca0d02bb63fdc5bc Mon Sep 17 00:00:00 2001 From: "Kevin R. Thornton" Date: Wed, 16 Aug 2023 08:28:33 -0700 Subject: [PATCH] Simplify gaussian stabilizing selection Python API. * Deprecate the current 4 classes. * Add fwdpy11.GaussianStabilizingSelection to replace * Update tests * Update manual --- cpp/CMakeLists.txt | 1 + .../GaussianStabilizingSelection.cc | 19 ++++ cpp/genetic_value_to_fitness/init.cc | 2 + doc/long_vignettes/gss_vignette.md | 6 +- doc/long_vignettes/tskitconvert_vignette.md | 14 ++- doc/misc/changelog.md | 24 ++++- doc/pages/advancedtopics.md | 3 +- doc/pages/gvalue_to_fitness.md | 18 +--- doc/short_vignettes/evolvets_vignette.md | 2 +- doc/short_vignettes/gvalues_traits.md | 30 +++--- ...iate_gaussian_effect_sizes_across_demes.md | 3 +- .../neutralmutsafter_vignette.md | 2 +- doc/short_vignettes/workingexample_trait.md | 2 +- doc/technical/genetic_values.md | 2 +- fwdpy11/__init__.py | 1 + fwdpy11/genetic_values.py | 93 ++++++++++++++++--- .../GaussianStabilizingSelection.hpp | 44 +++++++++ tests/demes-spec | 2 +- tests/test_add_mutation.py | 2 +- tests/test_dgvalue_pointer_vector.py | 3 +- tests/test_genetic_value_to_fitness.py | 36 +++++-- tests/test_genetic_values.py | 28 +++--- tests/test_inherited_noise.py | 3 +- tests/test_multivariate_effects.py | 2 +- tests/test_pickling.py | 3 +- tests/test_python_genetic_values.py | 9 +- tests/test_record_genetic_value_matrix.py | 4 +- tests/test_recreate_pop_from_own_treeseq.py | 2 +- tests/test_tree_sequences.py | 17 ++-- ...uences_different_des_in_different_demes.py | 2 +- 30 files changed, 268 insertions(+), 111 deletions(-) create mode 100644 cpp/genetic_value_to_fitness/GaussianStabilizingSelection.cc create mode 100644 fwdpy11/headers/fwdpy11/genetic_value_to_fitness/GaussianStabilizingSelection.hpp diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index aae4e7fb30..317541158a 100644 --- 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 new file mode 100644 index 0000000000..ec5552f4aa --- /dev/null +++ b/cpp/genetic_value_to_fitness/GaussianStabilizingSelection.cc @@ -0,0 +1,19 @@ +#include +#include +#include +#include + +namespace py = pybind11; + +void +init_GaussianStabilizingSelection(py::module& m) +{ + py::class_( + 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 index 5b25895e37..25da5e33ea 100644 --- 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 index 8232942306..77c4f352d3 100644 --- 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 index b8f9a8a6c4..2ac36e2b50 100644 --- 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 index 3e1585a5c5..4610a3aaa7 100644 --- 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 index fa69ab3303..af31872398 100644 --- 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 index 6484b37fd1..93f18f7a21 100644 --- 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 index 6eb26d80a5..bb07c69857 100644 --- 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 index 467209a7bb..b9901f5acc 100644 --- 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 `. :::{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 index 666a3dce0a..a1d4ef9164 100644 --- 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 index eec0c64d81..6735cd6df2 100644 --- 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 index b0f07c4ca1..35d48b950a 100644 --- 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 index 5ef4a768ef..fcc55ed650 100644 --- 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 index da9063c63e..ef95c87697 100644 --- a/fwdpy11/__init__.py +++ b/fwdpy11/__init__.py @@ -48,6 +48,7 @@ from .genetic_values import ( # NOQA PleiotropicOptima, Optimum, + GaussianStabilizingSelection, GSS, GSSmo, MultivariateGSS, diff --git a/fwdpy11/genetic_values.py b/fwdpy11/genetic_values.py index 18451aa90b..c5584743f0 100644 --- a/fwdpy11/genetic_values.py +++ b/fwdpy11/genetic_values.py @@ -17,13 +17,17 @@ # along with fwdpy11. If not, see . # +import enum import typing +import warnings import attr import numpy as np from ._fwdpy11 import (GeneticValueIsTrait, GeneticValueNoise, _ll_Additive, - _ll_GaussianNoise, _ll_GBR, _ll_GSSmo, + _ll_GaussianNoise, _ll_GBR, + _ll_GaussianStabilizingSelection, + _ll_GSSmo, _ll_Multiplicative, _ll_MultivariateGSSmo, _ll_NoNoise, _ll_Optimum, _ll_PleiotropicOptima, _ll_StrictAdditiveMultivariateEffects) @@ -51,13 +55,6 @@ class Optimum(_ll_Optimum): :param when: The time when the optimum shifts :type when: int or None - .. note:: - - When used to model a stable optimum (e.g., - :class:`fwdpy11.GSS`), the `when` parameter is omitted. - The `when` parameter is used for moving optima - (:class:`fwdpy11.GSSmo`). - .. versionadded:: 0.7.1 .. versionchanged:: 0.8.0 @@ -98,13 +95,6 @@ class PleiotropicOptima(_ll_PleiotropicOptima): :param when: The time when the optimum shifts :type when: int or None - .. note:: - - When used to model stable optima (e.g., - :class:`fwdpy11.MultivariateGSS`), the `when` parameter is omitted. - The `when` parameter is used for moving optima - (:class:`fwdpy11.MultivariateGSSmo`). - .. versionadded:: 0.7.1 .. versionchanged:: 0.8.0 @@ -137,6 +127,71 @@ def __eq__(self, other): return optima_equal and VS_equal and when_equal +@attr.s(auto_attribs=True, frozen=True, repr_ns="fwdpy11") +class GaussianStabilizingSelection(_ll_GaussianStabilizingSelection): + """ + Define a mapping of phenotype-to-fitness according to a + Gaussian stabilizing selection model. + + Instances of this trait must be constructed by one of the + various class methods available. + """ + is_single_trait: bool + optima: list + + def __attrs_post_init__(self): + if self.optima != sorted(self.optima, key=lambda x: x.when): + raise ValueError("optima must be sorted by time from past to present") + if self.is_single_trait is True: + inner = _ll_GSSmo(self.optima) + else: + inner = _ll_MultivariateGSSmo(self.optima) + super(GaussianStabilizingSelection, self).__init__(inner) + + def __getstate__(self): + return self.asdict() + + def __setstate__(self, d): + self.__dict__.update(**d) + self.__attrs_post_init__() + + @classmethod + def single_trait(cls, optima: typing.List[Optimum]): + """ + Stabilizing selection on a single trait + + :param optima: The optimum values. + Multiple values specify a moving optimum. + :type optima: List[fwdpy11.Optimum] + + """ + new_optima = [] + for o in optima: + if o.when is None: + new_optima.append(Optimum(o.optimum, o.VS, 0)) + else: + new_optima.append(o) + return cls(is_single_trait=True, optima=new_optima) + + @classmethod + def pleiotropy(cls, optima: typing.List[PleiotropicOptima]): + """ + Stabilizing selection with pleiotropy. + + :param optima: The optimum values. + Multiple values specify a moving optimum. + :type optima: List[fwdpy11.PleiotropicOptima] + """ + return cls(is_single_trait=False, optima=optima) + + def asdict(self): + return {"is_single_trait": self.is_single_trait, "optima": self.optima} + + @classmethod + def fromdict(cls, data): + return cls(**data) + + @attr_add_asblack @attr_class_to_from_dict_no_recurse @attr.s(auto_attribs=True, frozen=True, repr_ns="fwdpy11") @@ -172,6 +227,8 @@ class GSS(_ll_GSSmo): VS: typing.Optional[float] = None def __attrs_post_init__(self): + warnings.warn("use GaussianStabilizingSelection instead", + DeprecationWarning) if self.VS is None: super(GSS, self).__init__( [Optimum(optimum=self.optimum.optimum, @@ -237,6 +294,8 @@ def validate_optima(self, attribute, value): raise ValueError("Optimum.when is None") def __attrs_post_init__(self): + warnings.warn("use GaussianStabilizingSelection instead", + DeprecationWarning) super(GSSmo, self).__init__(self.optima) @@ -286,6 +345,8 @@ class MultivariateGSS(_ll_MultivariateGSSmo): VS: typing.Optional[float] = None def __attrs_post_init__(self): + warnings.warn("use GaussianStabilizingSelection instead", + DeprecationWarning) if self.VS is None: super(MultivariateGSS, self).__init__([self.optima]) else: @@ -342,6 +403,8 @@ def validate_optima(self, attribute, value): raise ValueError("PleiotropicOptima.when is None") def __attrs_post_init__(self): + warnings.warn("use GaussianStabilizingSelection instead", + DeprecationWarning) super(MultivariateGSSmo, self).__init__(self.optima) diff --git a/fwdpy11/headers/fwdpy11/genetic_value_to_fitness/GaussianStabilizingSelection.hpp b/fwdpy11/headers/fwdpy11/genetic_value_to_fitness/GaussianStabilizingSelection.hpp new file mode 100644 index 0000000000..893bb0d9ea --- /dev/null +++ b/fwdpy11/headers/fwdpy11/genetic_value_to_fitness/GaussianStabilizingSelection.hpp @@ -0,0 +1,44 @@ +#pragma once + +#include + +#include "GeneticValueIsTrait.hpp" +#include "GSSmo.hpp" +#include "MultivariateGSSmo.hpp" + +namespace fwdpy11 +{ + class GaussianStabilizingSelection : public GeneticValueIsTrait + { + private: + std::shared_ptr pimpl; + + public: + explicit GaussianStabilizingSelection(const GSSmo &input) + : GeneticValueIsTrait(input.total_dim), pimpl(input.clone()) + { + } + + explicit GaussianStabilizingSelection(const MultivariateGSSmo &input) + : GeneticValueIsTrait(input.total_dim), pimpl(input.clone()) + { + } + + double + operator()(const DiploidGeneticValueToFitnessData data) const final { + return this->pimpl->operator()(data); + } + + void + update(const DiploidPopulation &pop) final + { + return this->pimpl->update(pop); + } + + std::shared_ptr + clone() const final + { + return this->pimpl->clone(); + } + }; +} diff --git a/tests/demes-spec b/tests/demes-spec index 402bbbca16..44803a99f2 160000 --- a/tests/demes-spec +++ b/tests/demes-spec @@ -1 +1 @@ -Subproject commit 402bbbca16ba9d0da520ab6b17089dc804708aae +Subproject commit 44803a99f28209560948acab1a6b4106ec2c0283 diff --git a/tests/test_add_mutation.py b/tests/test_add_mutation.py index 9989b32c03..2879dc8354 100644 --- a/tests/test_add_mutation.py +++ b/tests/test_add_mutation.py @@ -92,7 +92,7 @@ def set_up_quant_trait_model(): # mu = theta/(4*N) r = rho / (4 * N) Opt = fwdpy11.Optimum - GSSmo = fwdpy11.GSSmo( + GSSmo = fwdpy11.GaussianStabilizingSelection.single_trait( [Opt(when=0, optimum=0.0, VS=1.0), Opt(when=N, optimum=1.0, VS=1.0)] ) a = fwdpy11.Additive(2.0, GSSmo) diff --git a/tests/test_dgvalue_pointer_vector.py b/tests/test_dgvalue_pointer_vector.py index 7fe3b7b327..c3480c150c 100644 --- a/tests/test_dgvalue_pointer_vector.py +++ b/tests/test_dgvalue_pointer_vector.py @@ -31,7 +31,8 @@ def test_construct_from_empty_list(self): def test_invalid_dimensionality(self): a = fwdpy11.Additive(2.0) - mvgss = fwdpy11.MultivariateGSS(np.zeros(2), 10.0) + mvgss = fwdpy11.GaussianStabilizingSelection.pleiotropy( + [fwdpy11.PleiotropicOptima(np.zeros(2), 10.0, when=0)]) mv = fwdpy11.StrictAdditiveMultivariateEffects(2, 0, mvgss) from fwdpy11._fwdpy11 import _dgvalue_pointer_vector diff --git a/tests/test_genetic_value_to_fitness.py b/tests/test_genetic_value_to_fitness.py index 455cf80c48..011605958e 100644 --- a/tests/test_genetic_value_to_fitness.py +++ b/tests/test_genetic_value_to_fitness.py @@ -43,15 +43,16 @@ class testGSS(unittest.TestCase): def setUp(self): self.VS = 1.0 self.opt = 0.0 - self.g = fwdpy11.GSS(VS=self.VS, optimum=self.opt) + self.g = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=self.VS, optimum=self.opt, when=None)]) def testPickle(self): import pickle p = pickle.dumps(self.g, -1) up = pickle.loads(p) - self.assertEqual(up.VS, self.VS) - self.assertEqual(up.optimum, self.opt) + self.assertEqual(up.optima[0].VS, self.VS) + self.assertEqual(up.optima[0].optimum, self.opt) def test_mapping(self): self.assertEqual(self.g.maps_to_fitness, False) @@ -66,7 +67,7 @@ def setUp(self): Opt(when=0, optimum=0.0, VS=1.0), Opt(when=100, optimum=1.0, VS=1.0), ] - self.g = fwdpy11.GSSmo(self.optima) + self.g = fwdpy11.GaussianStabilizingSelection.single_trait(self.optima) def testPickle(self): import pickle @@ -92,8 +93,9 @@ def setUp(self): timepoints = [0, 100] po = [] for i, t in enumerate(timepoints): - po.append(fwdpy11.PleiotropicOptima(when=t, optima=optima[i, :], VS=1.0)) - self.mvgssmo = fwdpy11.MultivariateGSSmo(po) + po.append(fwdpy11.PleiotropicOptima( + when=t, optima=optima[i, :], VS=1.0)) + self.mvgssmo = fwdpy11.GaussianStabilizingSelection.pleiotropy(po) def test_pickle(self): import pickle @@ -107,5 +109,27 @@ def test_mapping(self): self.assertEqual(self.mvgssmo.maps_to_trait_value, True) +def test_gaussian_stabilizing_selection_single_trait(): + gss = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(0.0, 10.0, None)]) + d = gss.asdict() + gss2 = gss.fromdict(d) + assert gss == gss2 + s = str(d) + gss2 = gss.fromdict(eval(s)) + assert gss == gss2 + + +def test_gaussian_stabilizing_selection_pleiotropy(): + gss = fwdpy11.GaussianStabilizingSelection.pleiotropy( + [fwdpy11.PleiotropicOptima([0, 0, 0, 0], 1.0, 0)]) + d = gss.asdict() + gss2 = gss.fromdict(d) + assert gss == gss2 + s = str(d) + gss2 = gss.fromdict(eval(s)) + assert gss == gss2 + + if __name__ == "__main__": unittest.main() diff --git a/tests/test_genetic_values.py b/tests/test_genetic_values.py index 38189c19e2..05925479aa 100644 --- a/tests/test_genetic_values.py +++ b/tests/test_genetic_values.py @@ -8,9 +8,10 @@ class testAdditive(unittest.TestCase): def setUp(self): GN = fwdpy11.GaussianNoise self.w = fwdpy11.Additive(2.0) - self.t = fwdpy11.Additive(2.0, fwdpy11.GSS(0.0, 1.0)) + self.t = fwdpy11.Additive(2.0, fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(0.0, 1.0, 0)])) self.tn = fwdpy11.Additive( - 1.0, fwdpy11.GSS(0.0, 1.0), GN(mean=0.1, sd=2.0)) + 1.0, fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(0.0, 1.0, 0)]), GN(mean=0.1, sd=2.0)) def testScaling(self): self.assertEqual(self.w.scaling, 2.0) @@ -74,13 +75,15 @@ def testPickleTraitWithNoiseToFile(self): class testMultiplicative(unittest.TestCase): - @classmethod + @ classmethod def setUp(self): GN = fwdpy11.GaussianNoise self.w = fwdpy11.Multiplicative(2.0) - self.t = fwdpy11.Multiplicative(2.0, fwdpy11.GSS(0.0, 1.0)) + self.t = fwdpy11.Multiplicative( + 2.0, fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(0.0, 1.0)])) self.tn = fwdpy11.Multiplicative( - 1.0, fwdpy11.GSS(0.0, 1.0), GN(mean=0.1, sd=2.0) + 1.0, fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(0.0, 1.0, 0)]), GN(mean=0.1, sd=2.0) ) def testScaling(self): @@ -131,9 +134,10 @@ def testPickleTraitWithNoise(self): class testGBR(unittest.TestCase): - @classmethod + @ classmethod def setUp(self): - self.gss = fwdpy11.GSS(0.0, 1.0) + self.gss = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(0.0, 1.0, 0)]) self.gnoise = fwdpy11.GaussianNoise(mean=0.0, sd=1.0) self.nonoise = fwdpy11.NoNoise() @@ -160,15 +164,5 @@ def testPicklingGSSGaussianNoise(self): self.assertEqual(type(self.gss), type(up.gvalue_to_fitness)) -class testGSS(unittest.TestCase): - @classmethod - def setUp(self): - self.x = fwdpy11.GSS(0.0, 1.0) - - def testProperties(self): - self.assertEqual(self.x.optimum, 0.0) - self.assertEqual(self.x.VS, 1.0) - - if __name__ == "__main__": unittest.main() diff --git a/tests/test_inherited_noise.py b/tests/test_inherited_noise.py index 45854a3482..ea66e44aa2 100644 --- a/tests/test_inherited_noise.py +++ b/tests/test_inherited_noise.py @@ -41,7 +41,8 @@ def setUp(self): [100] * 3, 1), "simlen": 3, "gvalue": fwdpy11.Additive( - 2.0, fwdpy11.GSS(optimum=0.0, VS=1.0), InheritedNoise() + 2.0, fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(optimum=0.0, VS=1.0, when=0)]), InheritedNoise() ), } diff --git a/tests/test_multivariate_effects.py b/tests/test_multivariate_effects.py index 12c77d7272..07da2cfd1e 100644 --- a/tests/test_multivariate_effects.py +++ b/tests/test_multivariate_effects.py @@ -17,7 +17,7 @@ def set_up_quant_trait_model(): po = [] for i, j in enumerate(timepoints): po.append(PO(when=int(j), optima=optima[i, :], VS=1.0)) - GSSmo = fwdpy11.MultivariateGSSmo(po) + GSSmo = fwdpy11.GaussianStabilizingSelection.pleiotropy(po) cmat = np.identity(ntraits) np.fill_diagonal(cmat, 0.1) a = fwdpy11.StrictAdditiveMultivariateEffects(ntraits, 0, GSSmo) diff --git a/tests/test_pickling.py b/tests/test_pickling.py index b480eabf11..0abddb1e86 100644 --- a/tests/test_pickling.py +++ b/tests/test_pickling.py @@ -63,7 +63,8 @@ def setUp(self): self.mu = self.theta / (4 * self.N) self.r = self.rho / (4 * self.N) - self.GSS = fwdpy11.GSS(VS=1.0, optimum=0.0) + self.GSS = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=1.0, optimum=0.0, when=0)]) a = fwdpy11.Additive(2.0, self.GSS) demography = fwdpy11.ForwardDemesGraph.tubes([self.N], burnin=100, diff --git a/tests/test_python_genetic_values.py b/tests/test_python_genetic_values.py index 6bc7252bc7..266d416b7d 100644 --- a/tests/test_python_genetic_values.py +++ b/tests/test_python_genetic_values.py @@ -93,8 +93,8 @@ def test_run(self): # ndemes=2, scaling=2, gvalue_to_fitness= ndemes=2, scaling=2, - gvalue_to_fitness=fwdpy11.GSS( - fwdpy11.Optimum(optimum=0.0, VS=1.0)), + gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait([ + fwdpy11.Optimum(optimum=0.0, VS=1.0)]), ) params = fwdpy11.ModelParams(**pdict) pop2 = fwdpy11.DiploidPopulation([1000, 1000], 1.0) @@ -124,7 +124,8 @@ def test_run(self): import pynoise - GSS = fwdpy11.GSS(optimum=0.0, VS=1.0) + GSS = fwdpy11.fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(optimum=0.0, VS=1.0, when=0)]) Noise = pynoise.PyNoise() pdict = build_model(GSS, Noise, 5, [1000, 1000]) params = fwdpy11.ModelParams(**pdict) @@ -202,7 +203,7 @@ def test_run(self): md = np.array(pop.diploid_metadata) pdict["gvalue"] = fwdpy11.Additive( - 2.0, fwdpy11.GSS(optimum=0.0, VS=1.0)) + 2.0, fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=0.0, VS=1.0, when=0)])) params = fwdpy11.ModelParams(**pdict) pop2 = fwdpy11.DiploidPopulation(1000, 1.0) rng = fwdpy11.GSLrng(1010) diff --git a/tests/test_record_genetic_value_matrix.py b/tests/test_record_genetic_value_matrix.py index b2d160f7d9..76dcaa0f1d 100644 --- a/tests/test_record_genetic_value_matrix.py +++ b/tests/test_record_genetic_value_matrix.py @@ -38,7 +38,7 @@ def set_up_quant_trait_model(): rho = 1.0 r = rho / (4 * N) Opt = fwdpy11.Optimum - GSSmo = fwdpy11.GSSmo( + GSSmo = fwdpy11.GaussianStabilizingSelection.single_trait( [Opt(when=0, optimum=0.0, VS=1.0), Opt( when=N, optimum=1.0, VS=1.0)] ) @@ -72,7 +72,7 @@ def set_up_two_trait_quant_trait_model(): when=N, optima=np.array([np.sqrt(2.0), 0]), VS=2.0 ), ] - GSSmo = fwdpy11.MultivariateGSSmo(optima) + GSSmo = fwdpy11.GaussianStabilizingSelection.pleiotropy(optima) a = fwdpy11.StrictAdditiveMultivariateEffects(2, 0, GSSmo) vcov = np.identity(2) np.fill_diagonal(vcov, 0.25) diff --git a/tests/test_recreate_pop_from_own_treeseq.py b/tests/test_recreate_pop_from_own_treeseq.py index 2e52946a56..4710399e6c 100644 --- a/tests/test_recreate_pop_from_own_treeseq.py +++ b/tests/test_recreate_pop_from_own_treeseq.py @@ -69,7 +69,7 @@ def pop_rng_params(): "recregions": [fwdpy11.BinomialInterval(0, L, 1)], "rates": (0.0, mu, None), "gvalue": fwdpy11.Additive( - scaling=2, gvalue_to_fitness=fwdpy11.GSS(optimum=optimum, VS=VS) + scaling=2, gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=optimum, VS=VS, when=0)]) ), "simlen": simlen, "prune_selected": False, diff --git a/tests/test_tree_sequences.py b/tests/test_tree_sequences.py index 10c1237a41..f271fe5859 100644 --- a/tests/test_tree_sequences.py +++ b/tests/test_tree_sequences.py @@ -59,7 +59,7 @@ def set_up_quant_trait_model(simlen=1.0): # mu = theta/(4*N) r = rho / (4 * N) Opt = fwdpy11.Optimum - GSSmo = fwdpy11.GSSmo( + GSSmo = fwdpy11.GaussianStabilizingSelection.single_trait( [Opt(when=0, optimum=0.0, VS=1.0), Opt(when=N, optimum=1.0, VS=1.0)] ) a = fwdpy11.Additive(2.0, GSSmo) @@ -998,7 +998,8 @@ def setUpClass(self): self.nreps = 500 self.mu = self.theta / (4 * self.N) self.r = self.rho / (4 * self.N) - self.GSS = fwdpy11.GSS(VS=1.0, optimum=0.0) + self.GSS = fwdpy11.fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=1.0, optimum=0.0, when=0)]) a = fwdpy11.Additive(2.0, self.GSS) self.p = { "nregions": [], @@ -1031,7 +1032,8 @@ def testQtraitSim(self): rho = 1.0 r = rho / (4 * N) - GSS = fwdpy11.GSS(VS=1.0, optimum=1.0) + GSS = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=1.0, optimum=1.0, when=0)]) a = fwdpy11.Additive(2.0, GSS) p = { "nregions": [], @@ -1115,7 +1117,8 @@ def testQtraitSim(self): rho = 1.0 r = rho / (4 * N) - GSS = fwdpy11.GSS(VS=1.0, optimum=1.0) + GSS = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=1.0, optimum=1.0, when=0)]) a = fwdpy11.Additive(2.0, GSS) p = { "nregions": [], @@ -1208,7 +1211,8 @@ def setUpClass(self): self.mu = self.theta / (4 * self.N) self.r = self.rho / (4 * self.N) - self.GSS = fwdpy11.GSS(VS=1.0, optimum=0.0) + self.GSS = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=1.0, optimum=0.0, when=0)]) a = fwdpy11.Additive(2.0, self.GSS) self.p = { "nregions": [], @@ -1353,7 +1357,8 @@ def __call__(self, pop): self.rho = 1.0 self.r = self.rho / (4 * self.N) - self.GSS = fwdpy11.GSS(VS=1.0, optimum=0.0) + self.GSS = fwdpy11.GaussianStabilizingSelection.single_trait( + [fwdpy11.Optimum(VS=1.0, optimum=0.0, when=0)]) a = fwdpy11.Additive(2.0, self.GSS) self.p = { "nregions": [], diff --git a/tests/test_tree_sequences_different_des_in_different_demes.py b/tests/test_tree_sequences_different_des_in_different_demes.py index 207cd0eb44..2767ff656b 100644 --- a/tests/test_tree_sequences_different_des_in_different_demes.py +++ b/tests/test_tree_sequences_different_des_in_different_demes.py @@ -135,7 +135,7 @@ def setUpClass(self): "demography": demography, "simlen": 100, "gvalue": fwdpy11.Additive( - ndemes=2, scaling=2, gvalue_to_fitness=fwdpy11.GSS(optimum=0.0, VS=1.0) + ndemes=2, scaling=2, gvalue_to_fitness=fwdpy11.GaussianStabilizingSelection.single_trait([fwdpy11.Optimum(optimum=0.0, VS=1.0, when=0)]) ), }