Skip to content

Commit

Permalink
Further prep for upcoming refactor. (#1067)
Browse files Browse the repository at this point in the history 
* Annotate existing code for stuff that will
  change.
* Duplicate another internal function and annotate
  it.
  • Loading branch information
@molpopgen
molpopgen committed Feb 22, 2023
1 parent 02ab084 commit 932b98e
Show file tree
Hide file tree
Showing 2 changed files with 147 additions and 4 deletions.
98 changes: 98 additions & 0 deletions fwdpy11/headers/fwdpy11/evolvets/evolve_generation_ts.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,5 +173,103 @@ namespace fwdpy11
throw std::runtime_error("error in book-keeping offspring nodes");
}
}
template <typename rng_t, typename poptype, typename genetic_param_holder>
void
evolve_generation_ts_refactor(
const rng_t& rng, poptype& pop, genetic_param_holder& genetics,
const fwdpy11::discrete_demography::DiscreteDemographyState&
current_demographic_state,
const fwdpy11::discrete_demography::multideme_fitness_lookups<std::uint32_t>&
fitness_lookup,
const fwdpy11::discrete_demography::migration_lookup& miglookup,
const fwdpp::uint_t generation, fwdpp::ts::edge_buffer& new_edge_buffer,
std::vector<fwdpy11::DiploidGenotype>& offspring,
std::vector<fwdpy11::DiploidMetadata>& offspring_metadata,
std::int32_t next_index)
{
fwdpp::debug::all_haploid_genomes_extant(pop);

genetics.haploid_genome_recycling_bin
= fwdpp::make_haploid_genome_queue(pop.haploid_genomes);

fwdpp::zero_out_haploid_genomes(pop);

offspring.clear();
offspring_metadata.clear();

// Generate the offspring
auto next_index_local = next_index;

// NOTE: this variable changes
using maxdeme_type = typename std::remove_const<
decltype(current_demographic_state.maxdemes)>::type;
// NOTE: this loop changes -- we can probably just iterate
// over the offspring deme sizes iterator
for (maxdeme_type deme = 0; deme < current_demographic_state.maxdemes; ++deme)
{
auto next_N_deme = current_demographic_state.current_deme_parameters
.next_deme_sizes.get()[deme];
for (decltype(next_N_deme) ind = 0; ind < next_N_deme; ++ind)
{
// NOTE: this API changes
// Get the parents
auto pdata = fwdpy11::discrete_demography::pick_parents(
rng, deme, miglookup,
current_demographic_state.current_deme_parameters
.current_deme_sizes,
current_demographic_state.current_deme_parameters
.selfing_rates,
current_demographic_state.fitness_bookmark, fitness_lookup);
fwdpy11::DiploidGenotype dip{
std::numeric_limits<std::size_t>::max(),
std::numeric_limits<std::size_t>::max()};
auto offspring_data = generate_offspring(
rng, std::make_pair(pdata.parent1, pdata.parent2), pop, dip,
genetics);
auto p1id = parent_nodes_from_metadata(
pdata.parent1, pop.diploid_metadata,
offspring_data.first.swapped);
auto p2id = parent_nodes_from_metadata(
pdata.parent2, pop.diploid_metadata,
offspring_data.second.swapped);
fwdpp::ts::table_index_t offspring_node_1
= fwdpp::ts::record_diploid_offspring(
offspring_data.first.breakpoints, p1id, deme, generation,
*pop.tables, new_edge_buffer);
fwdpp::ts::record_mutations_infinite_sites(
offspring_node_1, pop.mutations,
offspring_data.first.mutation_keys, *pop.tables);
fwdpp::ts::table_index_t offspring_node_2
= fwdpp::ts::record_diploid_offspring(
offspring_data.second.breakpoints, p2id, deme,
generation, *pop.tables, new_edge_buffer);
fwdpp::ts::record_mutations_infinite_sites(
offspring_node_2, pop.mutations,
offspring_data.second.mutation_keys, *pop.tables);

// Add metadata for the offspring
offspring_metadata.emplace_back(fwdpy11::DiploidMetadata{
0.0,
0.0,
1.,
{0, 0, 0},
offspring_metadata.size(),
{pdata.parent1, pdata.parent2},
deme,
0,
{offspring_node_1, offspring_node_2}});
offspring.emplace_back(std::move(dip));

next_index_local = offspring_node_2;
}
}
assert(static_cast<std::size_t>(next_index_local)
== pop.tables->num_nodes() - 1);
if (next_index_local
!= static_cast<decltype(next_index_local)>(pop.tables->num_nodes() - 1))
{
throw std::runtime_error("error in book-keeping offspring nodes");
}
}
} // namespace fwdpy11
#endif
53 changes: 49 additions & 4 deletions lib/evolve_discrete_demes/evolvets.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -730,6 +730,10 @@ evolve_with_tree_sequences_refactor(
}
}

// TODO: here, we need to simply set the model
// to the value of pop.generation, which means
// that pop.generation is a parental generation
// (which has always been true for fwdpy11).
if (pop.generation == 0)
{
// If we have an already-initialized model state,
Expand All @@ -748,22 +752,30 @@ evolve_with_tree_sequences_refactor(
}
}

// TODO: this goes away
auto current_demographic_state = demography.get_model_state();

// TODO: this is replaced by calling update_state
// on the ForwardGraph
current_demographic_state.initialize(pop);

// TODO: this changes to graph.number_of_demes() > 0
if (current_demographic_state.maxdemes <= 0)
{
throw std::runtime_error("maxdemes must be > 0");
}
if (gvalue_pointers.genetic_values.size()
// TODO: this changes to graph.number_of_demes()
> static_cast<std::size_t>(current_demographic_state.maxdemes))
{
throw std::invalid_argument(
"list of genetic values is longer than maxdemes");
}
// TODO: this changes to graph.number_of_demes()
if (static_cast<std::size_t>(current_demographic_state.maxdemes) > 1
&& gvalue_pointers.genetic_values.size() > 1
&& gvalue_pointers.genetic_values.size()
// TODO: this changes to graph.number_of_demes()
< static_cast<std::size_t>(current_demographic_state.maxdemes))
{
throw std::invalid_argument("too few genetic value objects");
Expand Down Expand Up @@ -794,6 +806,11 @@ evolve_with_tree_sequences_refactor(
auto genetics = fwdpp::make_genetic_parameters(gvalue_pointers.genetic_values,
std::move(bound_mmodel),
std::move(bound_rmodel));
// NOTE: this could be a bit tricky!
// A demes model gives ids according to order
// in the graph.
// We may need an API to allow a Python-side
// mapping of deme name -> genetic value object.
std::vector<std::size_t> deme_to_gvalue_map(current_demographic_state.maxdemes, 0);
if (genetics.gvalue.size() > 1)
{
Expand All @@ -811,6 +828,8 @@ evolve_with_tree_sequences_refactor(
std::vector<fwdpy11::DiploidMetadata> offspring_metadata(pop.diploid_metadata);
std::vector<fwdpy11::DiploidGenotype> offspring;
std::vector<double> new_diploid_gvalues;

// NOTE: this API is likely to break.
calculate_diploid_fitness(rng, pop, genetics.gvalue, deme_to_gvalue_map,
offspring_metadata, new_diploid_gvalues,
options.record_gvalue_matrix);
Expand All @@ -819,11 +838,16 @@ evolve_with_tree_sequences_refactor(

ddemog::multideme_fitness_lookups<std::uint32_t> fitness_lookup{
current_demographic_state.maxdemes};

// NOTE: this goes away and is replaced by ancestry
// proportions per offspring deme.
ddemog::migration_lookup miglookup{current_demographic_state.maxdemes,
current_demographic_state.M.empty()};

// TODO: this entire if/else may not be necessary?
if (pop.generation == 0)
{
// NOTE: the next 3 function calls all change
current_demographic_state.early(rng, pop.generation, pop.diploid_metadata);
current_demographic_state.late(rng, pop.generation, miglookup,
pop.diploid_metadata);
Expand All @@ -847,13 +871,19 @@ evolve_with_tree_sequences_refactor(
current_demographic_state.M);
}

// TODO: check that we can have all deme sizes be zero with a demes model?
// The specification says that a deme size has an exclusive minimum of zerol.
// So, no -- this cannot happen.
// Demes go extinct after their last epoch (forwards in time) has run out
// and the simulation is otherwise continuing.
if (current_demographic_state.will_go_globally_extinct() == true)
{
std::ostringstream o;
o << "extinction at time " << pop.generation;
throw ddemog::GlobalExtinction(o.str());
}

// TODO: Do we have sufficient test coverage through here?
if (!pop.mutations.empty())
{
// It is possible that pop already has a tree sequence
Expand Down Expand Up @@ -934,10 +964,10 @@ evolve_with_tree_sequences_refactor(
for (std::uint32_t gen = 0; gen < simlen && !stopping_criteron_met; ++gen)
{
++pop.generation;
fwdpy11::evolve_generation_ts(rng, pop, genetics, current_demographic_state,
fitness_lookup, miglookup, pop.generation,
*new_edge_buffer, offspring,
offspring_metadata, next_index);
fwdpy11::evolve_generation_ts_refactor(
rng, pop, genetics, current_demographic_state, fitness_lookup, miglookup,
pop.generation, *new_edge_buffer, offspring, offspring_metadata,
next_index);
// TODO: abstract out these steps into a "cleanup_pop" function
// NOTE: by swapping the diploids here, it is not possible
// for genetics.value to make use of parental genotype information.
Expand All @@ -950,6 +980,7 @@ evolve_with_tree_sequences_refactor(
// for a bit more context.
pop.diploids.swap(offspring);

// NOTE: this is a no-op
current_demographic_state.early(rng, pop.generation, offspring_metadata);

// NOTE: the two swaps of the metadata ensure
Expand Down Expand Up @@ -1067,14 +1098,26 @@ evolve_with_tree_sequences_refactor(
}
}

// NOTE: this is where the model state is getting updated
// We will replace this with updating the state of the
// ForwardGraph.
current_demographic_state.late(rng, pop.generation, miglookup,
pop.diploid_metadata);

// NOTE: API change needed
// Easy mode: overload this function to accept a new type.
fitness_lookup.update(current_demographic_state.fitness_bookmark);
// NOTE: new function needed:
// If a pop has nonzero ancestry from a deme, but that
// parental deme has no individuals, throw an exception.
ddemog::validate_parental_state(
pop.generation, fitness_lookup,
current_demographic_state.current_deme_parameters,
current_demographic_state.M);

// NOTE: this may be irrelevant -- see comment above about
// whether or not this is even possible.
// See above -- this block will likely just go away entirely.
if (current_demographic_state.will_go_globally_extinct() == true)
{
simplification(
Expand Down Expand Up @@ -1167,5 +1210,7 @@ evolve_with_tree_sequences_refactor(
<< __LINE__;
throw std::runtime_error(o.str());
}

// This is irrelevant -- yay!
demography.set_model_state(std::move(current_demographic_state));
}

0 comments on commit 932b98e

Please sign in to comment.