Program Listing for File add_mutation.hpp¶
↰ Return to documentation for file (fwdpp/sugar/add_mutation.hpp)
#ifndef FWDPP_SUGAR_ADD_MUTATION_HPP
#define FWDPP_SUGAR_ADD_MUTATION_HPP
#include <exception>
#include <type_traits>
#include <algorithm>
#include <stdexcept>
#include <vector>
#include <unordered_map>
#include <fwdpp/simfunctions/recycling.hpp>
#include <fwdpp/poptypes/tags.hpp>
#include <fwdpp/debug.hpp>
namespace fwdpp
{
namespace sugar
{
template <typename MutationContainerType, typename mcounts_t>
std::size_t
get_mut_index(MutationContainerType &mutations, mcounts_t &mcounts,
typename MutationContainerType::value_type &new_mutation)
{
// If the mutation already exists, let's not create it...
auto mitr
= std::find(std::begin(mutations), std::end(mutations), new_mutation);
if (mitr != mutations.end())
return std::size_t(std::distance(std::begin(mutations), mitr));
/*
If we are here, then new_mutation does not currently exist, so we
have to add it.
We do this using fwdpp's recycling method. But, we don't need to
create recycling bin,
as we're just looking to add one mutation. Instead, we can
simply search for an extinct
mutation.
*/
// find an extinct mutation, if one exists
auto extinct_mut = std::find(std::begin(mcounts), std::end(mcounts), 0);
std::size_t mindex = mcounts.size(); // this will be the correct
// value if we use the else
// block below
if (extinct_mut != std::end(mcounts)) // then we can recycle
{
auto dist = std::distance(std::begin(mcounts), extinct_mut);
mutations[dist]
= std::move(new_mutation); // move new mutation into place
mindex = std::size_t(dist); // update our value
}
else
{
// cannot recycle, so add it to end
mutations.emplace_back(std::move(new_mutation));
mcounts.push_back(0); // Add a place for this variant
}
return mindex;
}
template <typename poptype, typename map_t>
void
add_mutation_details(poptype &p, const std::vector<std::size_t> &mindexes,
const map_t &gams)
{
auto gam_recycling_bin = make_haploid_genome_queue(p.haploid_genomes);
// Function object for calls to upper bound
auto inserter
= [&p](const double &__value, const std::size_t __mut) noexcept {
return __value < p.mutations[__mut].pos;
};
for (const auto &gi : gams)
{
auto n = p.haploid_genomes[gi.first].mutations;
auto s = p.haploid_genomes[gi.first].smutations;
for (auto mindex : mindexes)
{
auto pos = p.mutations[mindex].pos;
if (p.mutations[mindex].neutral)
{
debug::validate_mutation_key_ranges(
p.mutations, n.begin(), n.end());
n.insert(std::upper_bound(n.begin(), n.end(), pos,
inserter),
mindex);
}
else
{
debug::validate_mutation_key_ranges(
p.mutations, s.begin(), s.end());
s.insert(std::upper_bound(s.begin(), s.end(), pos,
inserter),
mindex);
}
// update mutation count
p.mcounts[mindex] +=
typename poptype::mcount_t::value_type(gi.second.size());
}
// get new haploid_genome
auto new_haploid_genome_key = recycle_haploid_genome(
p.haploid_genomes, gam_recycling_bin, n, s);
// update haploid_genome count
p.haploid_genomes[gi.first].n
-= decltype(p.haploid_genomes[gi.first].n)(gi.second.size());
p.haploid_genomes[new_haploid_genome_key].n += decltype(
p.haploid_genomes[new_haploid_genome_key].n)(gi.second.size());
// This updates every diploid to have a key to
// p.haploid_genomes[new_haploid_genome_key]
for (auto i : gi.second)
{
*i = new_haploid_genome_key;
}
}
}
template <typename poptype,
typename
= std::enable_if<std::is_same<typename poptype::popmodel_t,
fwdpp::poptypes::DIPLOID_TAG>::value>>
std::unordered_map<std::size_t,
std::vector<typename poptype::diploid_type::first_type *>>
collect_haploid_genomes(poptype &p, const std::vector<std::size_t> &indlist,
const std::vector<short> &clist)
{
std::unordered_map<std::size_t,
std::vector<typename poptype::diploid_type::first_type *>>
gams;
for (std::size_t i = 0; i < indlist.size(); ++i)
{
if (clist[i] == 0 || clist[i] == 2)
{
gams[p.diploids[indlist[i]].first].push_back(
&p.diploids[indlist[i]].first);
}
if (clist[i] > 0)
{
gams[p.diploids[indlist[i]].second].push_back(
&p.diploids[indlist[i]].second);
}
}
return gams;
}
} // namespace sugar
template <typename poptype, class... Args>
std::size_t
add_mutation(poptype &p, const std::vector<std::size_t> &indlist,
const std::vector<short> &clist, Args &&... args)
{
static_assert(std::is_same<typename poptype::popmodel_t,
fwdpp::poptypes::DIPLOID_TAG>::value,
"poptype must be a single-locus object type");
// Before we go deep into creating objects, let's do some checks
for (const auto &i : indlist)
{
if (i >= p.diploids.size())
throw std::out_of_range(
"indlist contains elements > p.diploids.size()");
}
for (const auto &c : clist)
{
if (c < 0 || c > 2)
throw std::out_of_range("clist contains elements < 0 and/or > 2");
}
if (indlist.size() != clist.size())
throw std::invalid_argument("indlist and clist must be same length");
// create a new mutation
typename poptype::mutation_container::value_type new_mutant(
std::forward<Args>(args)...);
auto mindex = sugar::get_mut_index(p.mutations, p.mcounts, new_mutant);
auto gams = sugar::collect_haploid_genomes(p, indlist, clist);
sugar::add_mutation_details(p, {mindex}, gams);
return mindex;
}
template <typename poptype>
void
add_mutations(poptype &p, const std::vector<std::size_t> &indlist,
const std::vector<short> &clist,
const std::vector<std::size_t> &mutation_indexes)
{
static_assert(std::is_same<typename poptype::popmodel_t,
fwdpp::poptypes::DIPLOID_TAG>::value,
"poptype must be a diploid population ype");
// Before we go deep into creating objects, let's do some checks
for (const auto &i : indlist)
{
if (i >= p.diploids.size())
throw std::out_of_range(
"indlist contains elements > p.diploids.size()");
}
for (const auto &c : clist)
{
if (c < 0 || c > 2)
throw std::invalid_argument(
"clist contains elements < 0 and/or > 2");
}
if (indlist.size() != clist.size())
throw std::invalid_argument("indlist and clist must be same length");
for (auto mi : mutation_indexes)
{
if (mi >= p.mutations.size())
throw std::out_of_range("mutation key out of range");
}
if (p.mcounts.size() != p.mutations.size())
throw std::invalid_argument("p.mcounts.size() != p.mutations.size()");
auto gams = sugar::collect_haploid_genomes(p, indlist, clist);
sugar::add_mutation_details(p, mutation_indexes, gams);
}
} // namespace fwdpp
#endif