Skip to content

Commit

Permalink
Bug fix: fix handling of fixation pruning
Browse files Browse the repository at this point in the history 
Reorganize back end to properly handle
pruning in the generations when simplification
happens.
  • Loading branch information
@molpopgen
molpopgen committed Jul 18, 2024
1 parent 63c01e8 commit 973e017
Show file tree
Hide file tree
Showing 3 changed files with 246 additions and 101 deletions.
195 changes: 149 additions & 46 deletions cpptests/test_fixation_pruning_during_simulation.cc
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,11 @@
#include <algorithm>
#include <boost/test/tools/fpc_op.hpp>
#include <boost/test/tools/old/interface.hpp>
#include <boost/test/unit_test.hpp>

#include <boost/test/unit_test_suite.hpp>
#include <core/evolve_discrete_demes/evolvets.hpp>
#include <core/demes/forward_graph.hpp>
#include <cstdint>
#include <cwchar>
#include <exception>
#include <fwdpy11/evolvets/SampleRecorder.hpp>
#include <fwdpy11/genetic_values/dgvalue_pointer_vector.hpp>
#include <fwdpy11/regions/MutationRegions.hpp>
Expand All @@ -22,32 +19,92 @@
#include <memory>
#include <stdexcept>

#include "forward_demes_graph_fixtures.hpp"
#include "fwdpp/ts/simplification/simplification.hpp"
#include "fwdpy11/discrete_demography/exceptions.hpp"
#include "core/genetic_maps/regions.hpp"
#include "fwdpy11/genetic_values/DiploidMultiplicative.hpp"
#include "fwdpy11/mutation_dominance/MutationDominance.hpp"
#include "fwdpy11/regions/Sregion.hpp"

const char *demographic_model = R"(
description: single deme model
time_units: generations
demes:
- name: A
epochs:
- start_size: 1000
)";

void
sample_recorder_callback(const fwdpy11::DiploidPopulation &pop,
const fwdpy11::SampleRecorder &)
{
const auto validate_fixations = [](const fwdpy11::DiploidPopulation &pop,
std::size_t genome) {
for (const auto m : pop.haploid_genomes[genome].smutations)
{
if (pop.mcounts[m] == 2 * pop.N)
{
auto found = false;
for (std::size_t i = 0; i < pop.fixations.size(); ++i)
{
if (pop.fixations[i].g == pop.mutations[m].g
&& pop.fixations[i].pos == pop.mutations[m].pos
&& pop.fixations[i].s == pop.mutations[m].s)
{
found = true;
if (pop.fixation_times[i] != pop.generation)
{
throw std::runtime_error(
"fixation found in "
"genome "
"that did not occur "
"this "
"generation");
}
}
}
if (found == false)
{
throw std::runtime_error("fixation not found");
}
}
}
};
for (const auto diploid : pop.diploids)
{
validate_fixations(pop, diploid.first);
validate_fixations(pop, diploid.second);
}
}

fwdpy11::MutationRegions
strong_positive_selection()
{
auto nweights = std::vector<double>();
auto sweights = std::vector<double>(1.);
std::vector<std::unique_ptr<fwdpy11::Sregion>> nregions, sregions;

sregions.emplace_back(fwdpy11::ExpS(fwdpy11::Region(0., 10., 1.0, true, 0), 2.0,
sregions.emplace_back(fwdpy11::ExpS(fwdpy11::Region(0., 1e6, 1.0, true, 0), 1.0,
0.025, fwdpy11::fixed_dominance(1.))
.clone());
return fwdpy11::MutationRegions::create(0., nweights, sweights, nregions, sregions);
}

fwdpy11::GeneralizedGeneticMap
poisson_map()
{
auto p = fwdpy11_core::PoissonInterval(0., 1e6, 0.01, false);
std::vector<std::unique_ptr<fwdpy11::PoissonCrossoverGenerator>> pr;
pr.emplace_back(p.ll_clone());
std::vector<std::unique_ptr<fwdpy11::NonPoissonCrossoverGenerator>> nr;
return fwdpy11::GeneralizedGeneticMap(std::move(pr), std::move(nr));
}

struct common_setup
{
fwdpy11::GSLrng_t rng;
fwdpy11::DiploidPopulation pop;
fwdpy11::MutationRegions mregions;
fwdpy11::RecombinationRegions recregions;
fwdpy11::GeneralizedGeneticMap recregions;
fwdpy11::DiploidMultiplicative multiplicative;
fwdpy11::dgvalue_pointer_vector_ gvalue_ptrs;
fwdpy11::SampleRecorder recorder;
Expand All @@ -58,15 +115,16 @@ struct common_setup
evolve_with_tree_sequences_options options;

common_setup()
: rng{42}, pop{100, 10.0}, mregions{strong_positive_selection()},
recregions{0., {}}, multiplicative{1,
2.0,
fwdpy11::final_multiplicative_fitness(),
[](const auto) { return false; },
nullptr,
nullptr},
: rng{54321}, pop{1000, 1e6}, mregions{strong_positive_selection()},
recregions(poisson_map()),
multiplicative{1,
2.0,
fwdpy11::final_multiplicative_fitness(),
[](const auto) { return false; },
nullptr,
nullptr},
gvalue_ptrs{multiplicative}, recorder{},
forward_demes_graph(SingleDemeModel().yaml, 10000),
forward_demes_graph(std::string(demographic_model), 10000),
post_simplification_recorder{[](const fwdpy11::DiploidPopulation &) {}},
stopping_criterion{[](const fwdpy11::DiploidPopulation &, const bool) -> bool {
return false;
Expand All @@ -80,42 +138,87 @@ BOOST_AUTO_TEST_SUITE(test_remove_fixations_multiplicative)

BOOST_FIXTURE_TEST_CASE(test_remove_fixations, common_setup)
{
const auto sample_recorder_callback = [](const fwdpy11::DiploidPopulation &pop,
fwdpy11::SampleRecorder &) {
for (const auto diploid : pop.diploids)
{
for (const auto m : pop.haploid_genomes[diploid.first].smutations)
{
if (pop.mcounts[m] == 2 * pop.N)
{
for (std::size_t i = 0; i < pop.fixations.size(); ++i)
{
if (pop.fixations[i].g == pop.mutations[m].g
&& pop.fixations[i].pos
== pop.mutations[m].pos
&& pop.fixations[i].s == pop.mutations[m].s)
{
if (pop.fixation_times[i]
!= pop.generation)
{
throw std::runtime_error(
"fixation found in genome "
"that did not occur this "
"generation");
}
}
}
}
}
}
};
BOOST_REQUIRE_EQUAL(mregions.weights.size(), mregions.regions.size());
BOOST_REQUIRE_EQUAL(mregions.weights.size(), 1);
options.track_mutation_counts_during_sim = true;
BOOST_REQUIRE_EQUAL(options.preserve_selected_fixations, false);
BOOST_REQUIRE_EQUAL(options.track_mutation_counts_during_sim, true);
evolve_with_tree_sequences(rng, pop, recorder, 10, forward_demes_graph, 1000, 0.,
0.02, mregions, recregions, gvalue_ptrs,
1e-3, mregions, recregions, gvalue_ptrs,
sample_recorder_callback, stopping_criterion,
post_simplification_recorder, options);
BOOST_REQUIRE_EQUAL(pop.generation, 1000);
BOOST_REQUIRE_EQUAL(pop.fixations.size(), pop.fixation_times.size());
BOOST_REQUIRE(!pop.fixations.empty());
}

BOOST_FIXTURE_TEST_CASE(test_retain_fixations, common_setup)
{
BOOST_REQUIRE_EQUAL(mregions.weights.size(), mregions.regions.size());
BOOST_REQUIRE_EQUAL(mregions.weights.size(), 1);
options.track_mutation_counts_during_sim = true;
options.preserve_selected_fixations = true;
BOOST_REQUIRE_EQUAL(options.preserve_selected_fixations, true);
BOOST_REQUIRE_EQUAL(options.track_mutation_counts_during_sim, true);
// Disable: the recorder in the fixture will raise an exception
// when preserving selected fixations.
auto sample_recorder_callback
= [](const fwdpy11::DiploidPopulation &, const fwdpy11::SampleRecorder &) {};
evolve_with_tree_sequences(rng, pop, recorder, 10, forward_demes_graph, 1000, 0.,
1e-3, mregions, recregions, gvalue_ptrs,
sample_recorder_callback, stopping_criterion,
post_simplification_recorder, options);
BOOST_REQUIRE_EQUAL(pop.generation, 1000);
BOOST_REQUIRE_EQUAL(pop.fixations.size(), pop.fixation_times.size());
BOOST_REQUIRE(!pop.fixations.empty());
BOOST_REQUIRE(pop.fixations.size() > 1);
for (const auto &g : pop.haploid_genomes)
{
if (g.n > 0)
{
auto found = 0;
for (auto m : g.smutations)
{
if (pop.mcounts[m] == 2 * pop.N)
{
found += 1;
}
}
BOOST_REQUIRE_EQUAL(found, pop.fixations.size());
}
}
}

BOOST_FIXTURE_TEST_CASE(test_remove_fixations_simplify_each_generation, common_setup)
{
BOOST_REQUIRE_EQUAL(pop.N, 1000);
BOOST_REQUIRE_EQUAL(mregions.weights.size(), mregions.regions.size());
BOOST_REQUIRE_EQUAL(mregions.weights.size(), 1);
options.track_mutation_counts_during_sim = true;
BOOST_REQUIRE_EQUAL(options.preserve_selected_fixations, false);
BOOST_REQUIRE_EQUAL(options.track_mutation_counts_during_sim, true);
options.suppress_edge_table_indexing = !options.suppress_edge_table_indexing;
BOOST_REQUIRE_EQUAL(options.suppress_edge_table_indexing, true);
evolve_with_tree_sequences(rng, pop, recorder, 1, forward_demes_graph, 1000, 0.,
1e-3, mregions, recregions, gvalue_ptrs,
sample_recorder_callback, stopping_criterion,
post_simplification_recorder, options);
BOOST_REQUIRE_EQUAL(pop.generation, 1000);
BOOST_REQUIRE_EQUAL(pop.fixations.size(), pop.fixation_times.size());
BOOST_REQUIRE(!pop.fixations.empty());
}

BOOST_FIXTURE_TEST_CASE(test_remove_fixations_simplify_each_generation_index_tables_after_simplify, common_setup)
{
BOOST_REQUIRE_EQUAL(pop.N, 1000);
BOOST_REQUIRE_EQUAL(mregions.weights.size(), mregions.regions.size());
BOOST_REQUIRE_EQUAL(mregions.weights.size(), 1);
options.track_mutation_counts_during_sim = true;
BOOST_REQUIRE_EQUAL(options.preserve_selected_fixations, false);
BOOST_REQUIRE_EQUAL(options.track_mutation_counts_during_sim, true);
BOOST_REQUIRE_EQUAL(options.suppress_edge_table_indexing, false);
evolve_with_tree_sequences(rng, pop, recorder, 1, forward_demes_graph, 1000, 0.,
1e-3, mregions, recregions, gvalue_ptrs,
sample_recorder_callback, stopping_criterion,
post_simplification_recorder, options);
BOOST_REQUIRE_EQUAL(pop.generation, 1000);
Expand Down
Loading

0 comments on commit 973e017

Please sign in to comment.