From 0d36339e0c4509d4fb00854edfea3e81c7476cd8 Mon Sep 17 00:00:00 2001 From: "Kevin R. Thornton" Date: Sun, 7 Jul 2024 12:52:09 -0700 Subject: [PATCH] Fix handling of "super lethal" mutations. (#1319) * for multiplicative fitness, return 0 as soon as w <= 0. * extract C++ back end to locations that enable C++ testing in cpptests/ --- cpp/genetic_values/Additive.cc | 82 +----- cpp/genetic_values/Multiplicative.cc | 88 +------ .../cpp_neutral_benchmark.cc | 78 +----- cpptests/CMakeLists.txt | 1 + cpptests/test_evolvets.cc | 82 +----- cpptests/test_site_dependent_genetic_value.cc | 211 +++++++++++++++ .../genetic_values/DiploidAdditive.hpp | 96 +++++++ .../genetic_values/DiploidMultiplicative.hpp | 96 +++++++ .../fwdpp_wrappers/fwdpp_genetic_value.hpp | 16 +- .../site_dependent_genetic_value.hpp | 249 ++++++++++++++++++ 10 files changed, 687 insertions(+), 312 deletions(-) create mode 100644 cpptests/test_site_dependent_genetic_value.cc create mode 100644 fwdpy11/headers/fwdpy11/genetic_values/DiploidAdditive.hpp create mode 100644 fwdpy11/headers/fwdpy11/genetic_values/DiploidMultiplicative.hpp create mode 100644 fwdpy11/headers/fwdpy11/genetic_values/site_dependent_genetic_value.hpp diff --git a/cpp/genetic_values/Additive.cc b/cpp/genetic_values/Additive.cc index ffe78a6655..d544f8c817 100644 --- a/cpp/genetic_values/Additive.cc +++ b/cpp/genetic_values/Additive.cc @@ -21,97 +21,25 @@ #include #include -#include +#include #include #include -namespace -{ - struct single_deme_additive_het - { - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d += m.s * m.h; - } - }; - - struct multi_deme_additive_het - { - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d += m.esizes[deme] * m.heffects[deme]; - } - }; - - struct single_deme_additive_hom - { - double scaling; - single_deme_additive_hom(double s) : scaling(s) - { - } - - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d += scaling * m.s; - } - }; - - struct multi_deme_additive_hom - { - double scaling; - multi_deme_additive_hom(double s) : scaling(s) - { - } - - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d += scaling * m.esizes[deme]; - } - }; - - struct final_additive_trait - { - inline double - operator()(double d) const - { - return d; - } - }; - - struct final_additive_fitness - { - inline double - operator()(double d) const - { - return std::max(0.0, 1. + d); - } - }; - - using DiploidAdditive = fwdpy11::stateless_site_dependent_genetic_value_wrapper< - single_deme_additive_het, single_deme_additive_hom, multi_deme_additive_het, - multi_deme_additive_hom, 0>; -} - namespace py = pybind11; void init_Additive(py::module& m) { - py::class_(m, "_ll_Additive") + py::class_(m, "_ll_Additive") .def(py::init([](double scaling, const fwdpy11::GeneticValueIsTrait* gvalue_to_fitness, const fwdpy11::GeneticValueNoise* noise, std::size_t ndemes) { if (gvalue_to_fitness != nullptr) { - return DiploidAdditive(ndemes, scaling, final_additive_trait(), - gvalue_to_fitness, noise); + return fwdpy11::additive_trait_model(ndemes, scaling, + gvalue_to_fitness, noise); } - return DiploidAdditive(ndemes, scaling, final_additive_fitness(), - gvalue_to_fitness, noise); + return fwdpy11::additive_fitness_model(ndemes, scaling, noise); }), py::arg("scaling"), py::arg("gvalue_to_fitness"), py::arg("noise"), py::arg("ndemes")); diff --git a/cpp/genetic_values/Multiplicative.cc b/cpp/genetic_values/Multiplicative.cc index 40aca63c93..e8cdf382a7 100644 --- a/cpp/genetic_values/Multiplicative.cc +++ b/cpp/genetic_values/Multiplicative.cc @@ -17,104 +17,28 @@ // along with fwdpy11. If not, see . // #include -#include #include -#include +#include #include #include namespace py = pybind11; -namespace -{ - struct single_deme_multiplicative_het - { - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + m.s * m.h); - } - }; - - struct multi_deme_multiplicative_het - { - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + m.esizes[deme] * m.heffects[deme]); - } - }; - - struct single_deme_multiplicative_hom - { - double scaling; - single_deme_multiplicative_hom(double s) : scaling(s) - { - } - - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + scaling * m.s); - } - }; - - struct multi_deme_multiplicative_hom - { - double scaling; - multi_deme_multiplicative_hom(double s) : scaling(s) - { - } - - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + scaling * m.esizes[deme]); - } - }; - - struct final_multiplicative_trait - { - inline double - operator()(double d) const - { - return d - 1.0; - } - }; - - struct final_multiplicative_fitness - { - inline double - operator()(double d) const - { - return std::max(0.0, d); - } - }; - - using DiploidMultiplicative - = fwdpy11::stateless_site_dependent_genetic_value_wrapper< - single_deme_multiplicative_het, single_deme_multiplicative_hom, - multi_deme_multiplicative_het, multi_deme_multiplicative_hom, 1>; -} - void init_Multiplicative(py::module& m) { - py::class_(m, - "_ll_Multiplicative") + py::class_( + m, "_ll_Multiplicative") .def(py::init([](double scaling, const fwdpy11::GeneticValueIsTrait* gvalue_to_fitness, const fwdpy11::GeneticValueNoise* noise, std::size_t ndemes) { if (gvalue_to_fitness != nullptr) { - return DiploidMultiplicative(ndemes, scaling, - final_multiplicative_trait(), - gvalue_to_fitness, noise); + return fwdpy11::multiplicative_trait_model( + ndemes, scaling, gvalue_to_fitness, noise); } - return DiploidMultiplicative(ndemes, scaling, - final_multiplicative_fitness(), - gvalue_to_fitness, noise); + return fwdpy11::multiplicative_fitness_model(ndemes, scaling, noise); }), py::arg("scaling"), py::arg("gvalue_to_fitness"), py::arg("noise"), py::arg("ndemes")); diff --git a/cpp_neutral_benchmark/cpp_neutral_benchmark.cc b/cpp_neutral_benchmark/cpp_neutral_benchmark.cc index 5c21098123..202f9b9186 100644 --- a/cpp_neutral_benchmark/cpp_neutral_benchmark.cc +++ b/cpp_neutral_benchmark/cpp_neutral_benchmark.cc @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -21,73 +22,7 @@ namespace po = boost::program_options; // the Python library C++ code. They are needed // to define how Multiplicative works. -struct single_deme_multiplicative_het -{ - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + m.s * m.h); - } -}; - -struct multi_deme_multiplicative_het -{ - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + m.esizes[deme] * m.heffects[deme]); - } -}; - -struct single_deme_multiplicative_hom -{ - double scaling; - single_deme_multiplicative_hom(double s) : scaling(s) - { - } - - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + scaling * m.s); - } -}; - -struct multi_deme_multiplicative_hom -{ - double scaling; - multi_deme_multiplicative_hom(double s) : scaling(s) - { - } - - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d *= (1. + scaling * m.esizes[deme]); - } -}; - -struct final_multiplicative_trait -{ - inline double - operator()(double d) const - { - return d - 1.0; - } -}; - -struct final_multiplicative_fitness -{ - inline double - operator()(double d) const - { - return std::max(0.0, d); - } -}; - -using DiploidMultiplicative = fwdpy11::stateless_site_dependent_genetic_value_wrapper< - single_deme_multiplicative_het, single_deme_multiplicative_hom, - multi_deme_multiplicative_het, multi_deme_multiplicative_hom, 1>; +using fwdpy11::DiploidMultiplicative; struct RecordNothing // copied from internal code @@ -175,8 +110,9 @@ simulate(const command_line_options& options) fwdpy11_core::ForwardDemesGraph forward_demes_graph(o.str(), options.nsteps); - DiploidMultiplicative fitness(1, 2., final_multiplicative_fitness(), nullptr, - nullptr); + DiploidMultiplicative fitness( + 1, 2., fwdpy11::final_multiplicative_fitness(), + [](const double) { return false; }, nullptr, nullptr); fwdpy11::dgvalue_pointer_vector_ gvalue_pointers(fitness); fwdpy11::no_ancient_samples no_ancient_samples{}; @@ -199,9 +135,7 @@ simulate(const command_line_options& options) evolve_with_tree_sequences( rng, pop, sample_recorder, options.simplification_interval, forward_demes_graph, - options.nsteps, 0., 0., mmodel, genetic_map, gvalue_pointers, - [](const fwdpy11::DiploidPopulation& /*pop*/, fwdpy11::SampleRecorder& /*sr*/) { - }, + options.nsteps, 0., 0., mmodel, genetic_map, gvalue_pointers, no_ancient_samples, stopping_criteron, record_nothing, tsoptions); std::cout << pop.generation << ' ' << pop.N << ' ' << pop.tables->edges.size() << ' ' << pop.tables->nodes.size() << '\n'; diff --git a/cpptests/CMakeLists.txt b/cpptests/CMakeLists.txt index 7970dda2d7..b65e1ceb4b 100644 --- a/cpptests/CMakeLists.txt +++ b/cpptests/CMakeLists.txt @@ -8,6 +8,7 @@ SET(CPPTEST_SOURCES forward_demes_graph_fixtures.cc test_evolvets.cc test_core_genetic_map_regions.cc + test_site_dependent_genetic_value.cc ) add_executable(fwdpy11_cpp_tests ${CPPTEST_SOURCES}) diff --git a/cpptests/test_evolvets.cc b/cpptests/test_evolvets.cc index 1f2340c176..dfe3ea095a 100644 --- a/cpptests/test_evolvets.cc +++ b/cpptests/test_evolvets.cc @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -36,79 +37,7 @@ * correct */ -namespace -{ - // NOTE: all of this is a copy/paste - // from code internal to the pybind11 module. - - struct single_deme_additive_het - { - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d += m.s * m.h; - } - }; - - struct multi_deme_additive_het - { - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d += m.esizes[deme] * m.heffects[deme]; - } - }; - - struct single_deme_additive_hom - { - double scaling; - single_deme_additive_hom(double s) : scaling(s) - { - } - - inline void - operator()(double& d, const fwdpy11::Mutation& m) const - { - d += scaling * m.s; - } - }; - - struct multi_deme_additive_hom - { - double scaling; - multi_deme_additive_hom(double s) : scaling(s) - { - } - - inline void - operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const - { - d += scaling * m.esizes[deme]; - } - }; - - struct final_additive_trait - { - inline double - operator()(double d) const - { - return d; - } - }; - - struct final_additive_fitness - { - inline double - operator()(double d) const - { - return std::max(0.0, 1. + d); - } - }; - - using DiploidAdditive = fwdpy11::stateless_site_dependent_genetic_value_wrapper< - single_deme_additive_het, single_deme_additive_hom, multi_deme_additive_het, - multi_deme_additive_hom, 0>; -} +using fwdpy11::DiploidAdditive; template struct common_setup { @@ -128,7 +57,12 @@ template struct common_setup common_setup() : rng{42}, pop{100, 10.0}, mregions{{}, {}}, recregions{0., {}}, - additive{1, 2.0, final_additive_fitness(), nullptr, nullptr}, + additive{1, + 2.0, + fwdpy11::final_additive_fitness(), + [](const auto) { return false; }, + nullptr, + nullptr}, gvalue_ptrs{additive}, recorder{}, post_simplification_recorder{[](const fwdpy11::DiploidPopulation&) {}}, sample_recorder_callback{ diff --git a/cpptests/test_site_dependent_genetic_value.cc b/cpptests/test_site_dependent_genetic_value.cc new file mode 100644 index 0000000000..3e673fc8db --- /dev/null +++ b/cpptests/test_site_dependent_genetic_value.cc @@ -0,0 +1,211 @@ +#include "fwdpp/fundamental_types/haploid_genome.hpp" +#include "fwdpp/fundamental_types/mutation_base.hpp" +#include "fwdpy11/genetic_value_data/genetic_value_data.hpp" +#include "fwdpy11/genetic_values/site_dependent_genetic_value.hpp" +#include "fwdpy11/types/Diploid.hpp" +#include "fwdpy11/types/DiploidPopulation.hpp" +#include +#include +#include +#include +#include + +struct Mutation : public fwdpp::mutation_base +{ + double s; + double h; + Mutation(double pos, double s, double h) + : fwdpp::mutation_base(pos, 0, false), s(s), h(h) + { + } +}; + +struct Pop +{ + std::vector mutations; + std::vector genomes; + std::vector metadata; + + Pop() : mutations{}, genomes{}, metadata{} + { + } +}; + +struct Fwdpy11Pop +{ + fwdpy11::DiploidPopulation pop; + std::vector offspring_metadata; + + Fwdpy11Pop() : pop{1, 100.}, offspring_metadata() + { + } +}; + +BOOST_AUTO_TEST_SUITE(test_low_level_gvalue_calculation) + +BOOST_FIXTURE_TEST_CASE(test_multiplicative_fitness_0, Pop) +{ + mutations.emplace_back(1.0, -0.55, 1.0); + genomes.emplace_back(1); + genomes.emplace_back(1); + genomes[0].smutations.emplace_back(0); + BOOST_REQUIRE_EQUAL(genomes.size(), 2); + fwdpy11::site_dependent_genetic_value sdgv{[](const double) { return false; }}; + auto gv = sdgv( + genomes[0].smutations.begin(), genomes[0].smutations.end(), + genomes[1].smutations.begin(), genomes[1].smutations.end(), mutations, + [](double &w, const auto &m) { w *= (1. + 2.0 * m.s); }, + [](double &w, const auto &m) { w *= (1. + m.s * m.h); }, + [](const double w) { return std::max(w, 0.0); }, 1.0); + BOOST_REQUIRE_EQUAL(gv, 1.0 - 0.55); +} + +BOOST_FIXTURE_TEST_CASE(test_multiplicative_fitness_1, Pop) +{ + mutations.emplace_back(1.0, -1.1, 1.0); + genomes.emplace_back(1); + genomes.emplace_back(1); + genomes[0].smutations.emplace_back(0); + BOOST_REQUIRE_EQUAL(genomes.size(), 2); + fwdpy11::site_dependent_genetic_value sdgv{[](const double) { return false; }}; + auto gv = sdgv( + genomes[0].smutations.begin(), genomes[0].smutations.end(), + genomes[1].smutations.begin(), genomes[1].smutations.end(), mutations, + [](double &w, const auto &m) { w *= (1. + 2.0 * m.s); }, + [](double &w, const auto &m) { w *= (1. + m.s * m.h); }, + [](const double w) { return std::max(w, 0.0); }, 1.0); + BOOST_REQUIRE_EQUAL(gv, 0.0); +} + +BOOST_FIXTURE_TEST_CASE(test_multiplicative_fitness_2, Pop) +{ + mutations.emplace_back(1.0, -1.1, 1.0); + mutations.emplace_back(1.0, -1.1, 1.0); + genomes.emplace_back(1); + genomes.emplace_back(1); + genomes[0].smutations.emplace_back(0); + genomes[0].smutations.emplace_back(1); + BOOST_REQUIRE_EQUAL(genomes.size(), 2); + // Gives the previous behavior where two "more than lethal" + // variants give w > 0.0 when multiplied together + fwdpy11::site_dependent_genetic_value sdgv{[](const double) { return false; }}; + auto gv = sdgv( + genomes[0].smutations.begin(), genomes[0].smutations.end(), + genomes[1].smutations.begin(), genomes[1].smutations.end(), mutations, + [](double &w, const auto &m) { w *= (1. + 2.0 * m.s); }, + [](double &w, const auto &m) { w *= (1. + m.s * m.h); }, + fwdpy11::final_multiplicative_fitness{}, 1.0); + BOOST_REQUIRE(gv > 0.0); + + // Redo with a new closure that will return a fitness of 0.0 + // when we first see w <= 0.0 + sdgv + = fwdpy11::site_dependent_genetic_value{[](const double w) { return w <= 0.0; }}; + gv = sdgv( + genomes[0].smutations.begin(), genomes[0].smutations.end(), + genomes[1].smutations.begin(), genomes[1].smutations.end(), mutations, + [](double &w, const auto &m) { w *= (1. + 2.0 * m.s); }, + [](double &w, const auto &m) { w *= (1. + m.s * m.h); }, + fwdpy11::final_multiplicative_fitness{}, 1.0); + BOOST_REQUIRE_EQUAL(gv, 0.0); +} + +BOOST_FIXTURE_TEST_CASE(test_additive_fitness, Pop) +{ + mutations.emplace_back(1.0, -1.1, 1.0); + mutations.emplace_back(1.0, -1.1, 1.0); + genomes.emplace_back(1); + genomes.emplace_back(1); + genomes[0].smutations.emplace_back(0); + genomes[0].smutations.emplace_back(1); + BOOST_REQUIRE_EQUAL(genomes.size(), 2); + // Gives the previous behavior where two "more than lethal" + // variants give w > 0.0 when multiplied together + fwdpy11::site_dependent_genetic_value sdgv{[](const double) { return false; }}; + auto gv = sdgv( + genomes[0].smutations.begin(), genomes[0].smutations.end(), + genomes[1].smutations.begin(), genomes[1].smutations.end(), mutations, + [](double &w, const auto &m) { w += (2.0 * m.s); }, + [](double &w, const auto &m) { w += (m.s * m.h); }, + fwdpy11::final_additive_fitness(), 0.0); + BOOST_REQUIRE_EQUAL(gv, 0.0); + + // Additive models DIFFER!! + // If we use the same callback as for the multiplicative case, + // we incorrectly get a final fitness > 0. + sdgv + = fwdpy11::site_dependent_genetic_value{[](const double w) { return w <= 0.0; }}; + gv = sdgv( + genomes[0].smutations.begin(), genomes[0].smutations.end(), + genomes[1].smutations.begin(), genomes[1].smutations.end(), mutations, + [](double &w, const auto &m) { w += (m.s); }, + [](double &w, const auto &m) { w += (m.s * m.h); }, + fwdpy11::final_additive_fitness(), 0.0); + BOOST_REQUIRE(gv > 0.0); +} + +BOOST_AUTO_TEST_SUITE_END() + +// NOTE: the tests below rely on knowledge of how +// code that is not really part of the public API +// works and are thus fragile. + +BOOST_AUTO_TEST_SUITE(test_diploid_multiplicative) + +BOOST_FIXTURE_TEST_CASE(test_multiplicative_fitness_0, Fwdpy11Pop) +{ + pop.mutations.emplace_back(false, 0.1, -0.55, 1.0, 0); + pop.haploid_genomes[0].smutations.emplace_back(0); + pop.haploid_genomes.emplace_back(1); + BOOST_REQUIRE_EQUAL(pop.haploid_genomes.size(), 2); + auto m = fwdpy11::multiplicative_fitness_model(1, 2, nullptr); + offspring_metadata.emplace_back( + fwdpy11::DiploidMetadata{0., 0., 1., {0., 0., 0.}, 0, {0, 0}, 0, 0, {0, 0}}); + fwdpy11::GSLrng_t rng(42); + pop.diploids[0] = {0, 1}; + fwdpy11::DiploidGeneticValueData data{rng, + pop, + pop.diploid_metadata[0], + pop.diploid_metadata[0], + 0, + offspring_metadata[0]}; + m(data); + BOOST_REQUIRE_CLOSE(offspring_metadata[0].g, 1.0 - 0.55, 1e-3); + BOOST_REQUIRE_CLOSE(offspring_metadata[0].w, 1.0 - 0.55, 1e-3); +} + +BOOST_FIXTURE_TEST_CASE(test_multiplicative_fitness_1, Fwdpy11Pop) +{ + pop.mutations.emplace_back(false, 0.1, -1.55, 1.0, 0); + pop.mutations.emplace_back(false, 0.2, -1.55, 1.0, 0); + pop.haploid_genomes[0].smutations.emplace_back(0); + pop.haploid_genomes[0].smutations.emplace_back(1); + pop.haploid_genomes.emplace_back(1); + BOOST_REQUIRE_EQUAL(pop.haploid_genomes.size(), 2); + // This lambda will NOT exit when w first is <= 0.0 + fwdpy11::DiploidMultiplicative m( + 1, 2., fwdpy11::final_multiplicative_fitness(), + [](const double) { return false; }, nullptr, nullptr); + offspring_metadata.emplace_back( + fwdpy11::DiploidMetadata{0., 0., 1., {0., 0., 0.}, 0, {0, 0}, 0, 0, {0, 0}}); + fwdpy11::GSLrng_t rng(42); + pop.diploids[0] = {0, 1}; + BOOST_CHECK_EQUAL(pop.haploid_genomes[0].smutations.size(), 2); + fwdpy11::DiploidGeneticValueData data{rng, + pop, + pop.diploid_metadata[0], + pop.diploid_metadata[0], + 0, + offspring_metadata[0]}; + m(data); + BOOST_REQUIRE(offspring_metadata[0].g > 0.0); + BOOST_REQUIRE(offspring_metadata[0].w > 0.0); + + // we now pass in this correct lambda expression + m = fwdpy11::multiplicative_fitness_model(1, 2., nullptr); + m(data); + BOOST_REQUIRE_EQUAL(offspring_metadata[0].g, 0.); + BOOST_REQUIRE_EQUAL(offspring_metadata[0].w, 0.); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/fwdpy11/headers/fwdpy11/genetic_values/DiploidAdditive.hpp b/fwdpy11/headers/fwdpy11/genetic_values/DiploidAdditive.hpp new file mode 100644 index 0000000000..134a44859d --- /dev/null +++ b/fwdpy11/headers/fwdpy11/genetic_values/DiploidAdditive.hpp @@ -0,0 +1,96 @@ +#pragma once + +#include +#include +#include "fwdpp_wrappers/fwdpp_genetic_value.hpp" + +namespace fwdpy11 +{ + struct single_deme_additive_het + { + inline void + operator()(double& d, const fwdpy11::Mutation& m) const + { + d += m.s * m.h; + } + }; + + struct multi_deme_additive_het + { + inline void + operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const + { + d += m.esizes[deme] * m.heffects[deme]; + } + }; + + struct single_deme_additive_hom + { + double scaling; + single_deme_additive_hom(double s) : scaling(s) + { + } + + inline void + operator()(double& d, const fwdpy11::Mutation& m) const + { + d += scaling * m.s; + } + }; + + struct multi_deme_additive_hom + { + double scaling; + multi_deme_additive_hom(double s) : scaling(s) + { + } + + inline void + operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const + { + d += scaling * m.esizes[deme]; + } + }; + + struct final_additive_trait + { + inline double + operator()(double d) const + { + return d; + } + }; + + struct final_additive_fitness + { + inline double + operator()(double d) const + { + return std::max(0.0, 1. + d); + } + }; + + using DiploidAdditive = fwdpy11::stateless_site_dependent_genetic_value_wrapper< + single_deme_additive_het, single_deme_additive_hom, multi_deme_additive_het, + multi_deme_additive_hom, 0>; + + inline DiploidAdditive + additive_fitness_model(std::size_t ndemes, double scaling, + const GeneticValueNoise* noise) + { + return DiploidAdditive( + ndemes, scaling, final_additive_fitness(), + [](const double) { return false; }, nullptr, noise); + } + + inline DiploidAdditive + additive_trait_model(std::size_t ndemes, double scaling, + const GeneticValueIsTrait* gvalue_to_fitness, + const GeneticValueNoise* noise) + { + return DiploidAdditive( + ndemes, scaling, final_additive_trait(), [](const double) { return false; }, + gvalue_to_fitness, noise); + } + +} diff --git a/fwdpy11/headers/fwdpy11/genetic_values/DiploidMultiplicative.hpp b/fwdpy11/headers/fwdpy11/genetic_values/DiploidMultiplicative.hpp new file mode 100644 index 0000000000..c0bb4416f3 --- /dev/null +++ b/fwdpy11/headers/fwdpy11/genetic_values/DiploidMultiplicative.hpp @@ -0,0 +1,96 @@ +#pragma once + +#include +#include +#include "fwdpp_wrappers/fwdpp_genetic_value.hpp" + +namespace fwdpy11 +{ + struct single_deme_multiplicative_het + { + inline void + operator()(double& d, const fwdpy11::Mutation& m) const + { + d *= (1. + m.s * m.h); + } + }; + + struct multi_deme_multiplicative_het + { + inline void + operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const + { + d *= (1. + m.esizes[deme] * m.heffects[deme]); + } + }; + + struct single_deme_multiplicative_hom + { + double scaling; + single_deme_multiplicative_hom(double s) : scaling(s) + { + } + + inline void + operator()(double& d, const fwdpy11::Mutation& m) const + { + d *= (1. + scaling * m.s); + } + }; + + struct multi_deme_multiplicative_hom + { + double scaling; + multi_deme_multiplicative_hom(double s) : scaling(s) + { + } + + inline void + operator()(const std::size_t deme, double& d, const fwdpy11::Mutation& m) const + { + d *= (1. + scaling * m.esizes[deme]); + } + }; + + struct final_multiplicative_trait + { + inline double + operator()(double d) const + { + return d - 1.0; + } + }; + + struct final_multiplicative_fitness + { + inline double + operator()(double d) const + { + return std::max(0.0, d); + } + }; + + using DiploidMultiplicative + = fwdpy11::stateless_site_dependent_genetic_value_wrapper< + single_deme_multiplicative_het, single_deme_multiplicative_hom, + multi_deme_multiplicative_het, multi_deme_multiplicative_hom, 1>; + + inline DiploidMultiplicative + multiplicative_fitness_model(std::size_t ndemes, double scaling, + const GeneticValueNoise* noise) + { + return DiploidMultiplicative( + ndemes, scaling, final_multiplicative_fitness(), + [](const double w) { return w <= 0.0; }, nullptr, noise); + } + + inline DiploidMultiplicative + multiplicative_trait_model(std::size_t ndemes, double scaling, + const GeneticValueIsTrait* gvalue_to_fitness, + const GeneticValueNoise* noise) + { + return DiploidMultiplicative( + ndemes, scaling, final_multiplicative_trait(), + [](const double) { return false; }, gvalue_to_fitness, noise); + } +} diff --git a/fwdpy11/headers/fwdpy11/genetic_values/fwdpp_wrappers/fwdpp_genetic_value.hpp b/fwdpy11/headers/fwdpy11/genetic_values/fwdpp_wrappers/fwdpp_genetic_value.hpp index fde032eaa9..66e0cfe07c 100644 --- a/fwdpy11/headers/fwdpy11/genetic_values/fwdpp_wrappers/fwdpp_genetic_value.hpp +++ b/fwdpy11/headers/fwdpy11/genetic_values/fwdpp_wrappers/fwdpp_genetic_value.hpp @@ -1,10 +1,11 @@ #ifndef FWDPY11_GENETIC_VALUES_WRAPPERS_FWDPP__GVALUE_HPP__ #define FWDPY11_GENETIC_VALUES_WRAPPERS_FWDPP__GVALUE_HPP__ +#include #include #include #include "../DiploidGeneticValue.hpp" -#include +#include #include namespace fwdpy11 @@ -26,7 +27,7 @@ namespace fwdpy11 } inline double - operator()(const fwdpp::site_dependent_genetic_value& gv, + operator()(const fwdpy11::site_dependent_genetic_value& gv, const std::size_t diploid_index, const DiploidMetadata& /*metadata*/, const DiploidPopulation& pop) const @@ -47,7 +48,7 @@ namespace fwdpy11 } inline double - operator()(const fwdpp::site_dependent_genetic_value& gv, + operator()(const fwdpy11::site_dependent_genetic_value& gv, const std::size_t diploid_index, const DiploidMetadata& metadata, const DiploidPopulation& pop) const { @@ -75,7 +76,7 @@ namespace fwdpy11 }; using callback_type = std::function; using make_return_value_t = std::function; @@ -90,7 +91,7 @@ namespace fwdpy11 return multi_deme_callback(aa_scaling); } - fwdpp::site_dependent_genetic_value gv; + fwdpy11::site_dependent_genetic_value gv; double aa_scaling; make_return_value_t make_return_value; callback_type callback; @@ -99,8 +100,9 @@ namespace fwdpy11 public: stateless_site_dependent_genetic_value_wrapper( std::size_t ndim, double scaling, make_return_value_t mrv, - const GeneticValueToFitnessMap* gv2w_, const GeneticValueNoise* noise_) - : DiploidGeneticValue{ndim, gv2w_, noise_}, gv{}, aa_scaling(scaling), + std::function clamp, const GeneticValueToFitnessMap* gv2w_, + const GeneticValueNoise* noise_) + : DiploidGeneticValue{ndim, gv2w_, noise_}, gv{clamp}, aa_scaling(scaling), make_return_value(std::move(mrv)), callback(init_callback(ndim, aa_scaling)), isfitness(gv2w->isfitness) { diff --git a/fwdpy11/headers/fwdpy11/genetic_values/site_dependent_genetic_value.hpp b/fwdpy11/headers/fwdpy11/genetic_values/site_dependent_genetic_value.hpp new file mode 100644 index 0000000000..2c23f09684 --- /dev/null +++ b/fwdpy11/headers/fwdpy11/genetic_values/site_dependent_genetic_value.hpp @@ -0,0 +1,249 @@ +#pragma once + +#include +#include + +namespace fwdpy11 +{ + struct site_dependent_genetic_value + { + std::function clamp; + //! The return value type + using result_type = double; + + template + inline result_type + operator()(iterator_t first1, iterator_t last1, iterator_t first2, + iterator_t last2, const MutationContainerType &mutations, + const updating_policy_hom &fpol_hom, + const updating_policy_het &fpol_het, + const make_return_value &rv_function, + const double starting_value) const noexcept + /*! + Range-based call operator. Calculates genetic values over ranges of + mutation keys first1/last1 + and first2/last2, which are iterators derived from the 'smutations' + of two haploid_genomes in a diploid. + + \param first1 Iterator to first mutation derived from haploid_genome 1 + \param last1 Iterator to one past the last mutation derived from + haploid_genome 1 + \param first2 Iterator to first mutation derived from haploid_genome 2 + \param last2 Iterator to one past the last mutation derived from + haploid_genome 2 + \param mutations The container of mutations for the simulation + \param fpol_hom Policy that updates genetic value for a homozygous + mutation. + \param fpol_het Policy that updates genetic value for a heterozygous + mutation. + \param make_return_value Policy generated the final return value. + Must be equivalent to std::function + \param starting_value The initial genetic value. + + \returns Fitness (double) + */ + { + static_assert(fwdpp::traits::is_mutation< + typename MutationContainerType::value_type>::value, + "MutationContainerType::value_type must be a mutation type"); + static_assert( + std::is_convertible< + updating_policy_hom, + std::function>::value, + "decltype(fpol_hom) must be convertible to " + "std::function>::value, + "decltype(fpol_het) must be convertible to " + "std::function + inline result_type + operator()(iterator_t first1, iterator_t last1, iterator_t first2, + iterator_t last2, const MutationContainerType &mutations, + const updating_policy_hom &fpol_hom, + const updating_policy_het &fpol_het, + const double starting_value) const noexcept + { + return this->operator()( + first1, last1, first2, last2, mutations, fpol_hom, fpol_het, + [](double d) { return d; }, starting_value); + } + + template + inline result_type + operator()(const HaploidGenomeType &g1, const HaploidGenomeType &g2, + const MutationContainerType &mutations, + const updating_policy_hom &fpol_hom, + const updating_policy_het &fpol_het, + const make_return_value &rv_function, + const double starting_value) const noexcept + /*! + Calculates genetic value for a diploid whose genotype + across sites is given by haploid_genomes g1 and g2. + + \param g1 A haploid_genome + \param g2 A haploid_genome + \param mutations The container of mutations for the simulation + \param fpol_hom Policy that updates genetic value for a homozygous + mutation. + \param fpol_het Policy that updates genetic value for a heterozygous + mutation. + \param starting_value The initial genetic value. + + \returns Fitness (double) + */ + { + static_assert(fwdpp::traits::is_haploid_genome::value, + "HaploidGenomeType::value_type must be a haploid_genome " + "type"); + return this->operator()(g1.smutations.cbegin(), g1.smutations.cend(), + g2.smutations.cbegin(), g2.smutations.cend(), + mutations, fpol_hom, fpol_het, rv_function, + starting_value); + } + + template + inline result_type + operator()(const HaploidGenomeType &g1, const HaploidGenomeType &g2, + const MutationContainerType &mutations, + const updating_policy_hom &fpol_hom, + const updating_policy_het &fpol_het, + const double starting_value) const noexcept + { + return this->operator()( + g1, g2, mutations, fpol_hom, fpol_het, [](double d) { return d; }, + starting_value); + } + + template + inline result_type + operator()(const DiploidType &dip, const GenomeContainerType &haploid_genomes, + const MutationContainerType &mutations, + const updating_policy_hom &fpol_hom, + const updating_policy_het &fpol_het, + const make_return_value &rv_function, + const double starting_value) const noexcept + /*! + Calculates genetic value for a diploid type. + + \param dip A diploid + \param haploid_genomes The container of haploid_genomes for the simulation. + \param mutations The container of mutations for the simulation + \param fpol_hom Policy that updates genetic value for a homozygous + mutation. + \param fpol_het Policy that updates genetic value for a heterozygous + mutation. + \param starting_value The initial genetic value. + + \returns Fitness (double) + */ + { + return this->operator()(haploid_genomes[dip.first], + haploid_genomes[dip.second], mutations, fpol_hom, + fpol_het, rv_function, starting_value); + } + + template + inline result_type + operator()(const DiploidType &dip, const GenomeContainerType &haploid_genomes, + const MutationContainerType &mutations, + const updating_policy_hom &fpol_hom, + const updating_policy_het &fpol_het, + const double starting_value) const noexcept + { + return this->operator()( + dip, haploid_genomes, mutations, fpol_hom, fpol_het, + [](double d) { return d; }, starting_value); + } + }; +}