Skip to content

Commit

Permalink
Enforce initial population deme sizes match the graph's expectations.
Browse files Browse the repository at this point in the history
  • Loading branch information
molpopgen committed Feb 25, 2023
1 parent 1397769 commit f81c1e3
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 14 deletions.
60 changes: 60 additions & 0 deletions cpptests/test_evolvets.cc
@@ -1,7 +1,9 @@
#include <boost/test/tools/old/interface.hpp>
#include <boost/test/unit_test.hpp>

#include <core/evolve_discrete_demes/evolvets.hpp>
#include <core/demes/forward_graph.hpp>
#include <cstdint>
#include <fwdpy11/evolvets/SampleRecorder.hpp>
#include <fwdpy11/genetic_values/dgvalue_pointer_vector.hpp>
#include <fwdpy11/regions/MutationRegions.hpp>
Expand All @@ -12,6 +14,7 @@
#include <fwdpy11/types/DiploidPopulation.hpp>

#include "forward_demes_graph_fixtures.hpp"
#include "fwdpy11/discrete_demography/exceptions.hpp"

BOOST_AUTO_TEST_SUITE(test_evolvets)

Expand Down Expand Up @@ -159,4 +162,61 @@ BOOST_FIXTURE_TEST_CASE(test_simlen_longer_than_model_length, common_setup)
BOOST_REQUIRE_EQUAL(pop.generation, 10);
}

BOOST_FIXTURE_TEST_CASE(test_invalid_deme_metadata, common_setup)
{
auto model = SingleDemeModel();
fwdpy11_core::ForwardDemesGraph forward_demes_graph(model.yaml, 10);

// Make some of the individuals in deme 1 but all must be in deme 0
// b/c that is the model
for (std::size_t i = 0; i < pop.diploid_metadata.size(); ++i)
{
if (i % 2 == 0)
{
pop.diploid_metadata[i].deme = 1;
}
}

BOOST_CHECK_THROW(
{
evolve_with_tree_sequences_refactor(
rng, pop, recorder, 10, forward_demes_graph, 10, 0., 0., mregions,
recregions, gvalue_ptrs, sample_recorder_callback, stopping_criterion,
post_simplification_recorder, options);
},
fwdpy11::discrete_demography::DemographyError);
// pop hasn't evolved!
BOOST_REQUIRE_EQUAL(pop.generation, 0);
}

BOOST_FIXTURE_TEST_CASE(test_initial_pop_size_invalid, common_setup)
{
// The demes model specifies N = 100.
// Here, we will start with a different N.
// This is a hard error!
// We test values ~100 because TDD found cases of memory errors
// when N >> the correct N but things can pass when N =~ the
// correct N. Gotta love UB :).
std::vector<std::uint32_t> initial_popsizes{50, 99, 101, 200};
for (auto initial_n : initial_popsizes)
{
// reset the fixture
pop = fwdpy11::DiploidPopulation(initial_n, 10.0);
auto model = SingleDemeModel();
fwdpy11_core::ForwardDemesGraph forward_demes_graph(model.yaml, 10);

BOOST_CHECK_THROW(
{
evolve_with_tree_sequences_refactor(
rng, pop, recorder, 10, forward_demes_graph, 10, 0., 0.,
mregions, recregions, gvalue_ptrs, sample_recorder_callback,
stopping_criterion, post_simplification_recorder, options);
},
// TODO: is this the type that we want?
fwdpy11::discrete_demography::DemographyError);
}
// pop hasn't evolved!
BOOST_REQUIRE_EQUAL(pop.generation, 0);
}

BOOST_AUTO_TEST_SUITE_END()
55 changes: 41 additions & 14 deletions lib/evolve_discrete_demes/evolvets.cc
@@ -1,4 +1,5 @@
#include <algorithm>
#include <cstdint>
#include <fwdpp/simparams.hpp>
#include <fwdpp/ts/simplify_tables.hpp>
#include <fwdpp/ts/simplify_tables_output.hpp>
Expand All @@ -14,6 +15,7 @@

#include "fwdpy11/discrete_demography/exceptions.hpp"
#include "fwdpy11/discrete_demography/simulation/multideme_fitness_bookmark.hpp"
#include "fwdpy11/types/Diploid.hpp"
#include "util.hpp"
#include "diploid_pop_fitness.hpp"
#include "index_and_count_mutations.hpp"
Expand All @@ -25,6 +27,7 @@
#include "runtime_checks.hpp"

#include <core/evolve_discrete_demes/evolvets.hpp>
#include <sstream>

namespace ddemog = fwdpy11::discrete_demography;

Expand Down Expand Up @@ -663,6 +666,43 @@ evolve_with_tree_sequences(
demography.set_model_state(std::move(current_demographic_state));
}

void
check_initial_deme_sizes(const std::vector<fwdpy11::DiploidMetadata> &metadata,
const fwdpy11_core::ForwardDemesGraph &demography)
{
auto parental_deme_sizes = demography.parental_deme_sizes();
auto num_extant_parental_demes
= std::count_if(std::begin(parental_deme_sizes), std::end(parental_deme_sizes),
[](auto a) { return a > 0.0; });
if (num_extant_parental_demes == 0)
{
throw fwdpy11::discrete_demography::DemographyError(
"all parental deme sizes are zero");
}
std::vector<std::uint32_t> current_deme_sizes(demography.number_of_demes(), 0);
for (const auto &md : metadata)
{
if (md.deme >= demography.number_of_demes() || md.deme < 0)
{
throw ddemog::DemographyError("individual has invalid deme");
}
current_deme_sizes[md.deme] += 1;
}
std::size_t i = 0;
for (auto psize : parental_deme_sizes)
{
auto intsize = static_cast<std::uint32_t>(psize);
if (intsize != current_deme_sizes[i])
{
std::ostringstream o;
o << "initial size of deme " << i << " is " << current_deme_sizes[i]
<< " but the required size is " << intsize;
throw ddemog::DemographyError(o.str());
}
++i;
}
}

void
evolve_with_tree_sequences_refactor(
const fwdpy11::GSLrng_t &rng, fwdpy11::DiploidPopulation &pop,
Expand Down Expand Up @@ -772,20 +812,7 @@ evolve_with_tree_sequences_refactor(
{
throw std::runtime_error("maxdemes must be > 0");
}
// NOTE: the else block may be unnecessary but we will keep
// it for now.
else
{
auto parental_deme_sizes = demography.parental_deme_sizes();
auto num_extant_parental_demes = std::count_if(
std::begin(parental_deme_sizes), std::end(parental_deme_sizes),
[](auto a) { return a > 0.0; });
if (num_extant_parental_demes == 0)
{
throw fwdpy11::discrete_demography::DemographyError(
"all parental deme sizes are zero");
}
}
check_initial_deme_sizes(pop.diploid_metadata, demography);

// NOTE: should this be != instead of > ??
if (gvalue_pointers.genetic_values.size()
Expand Down

0 comments on commit f81c1e3

Please sign in to comment.