Program Listing for File recycling.hpp¶
↰ Return to documentation for file (fwdpp/ts/recycling.hpp)
#ifndef FWDPP_TS_RECYCLING_HPP
#define FWDPP_TS_RECYCLING_HPP
#include <algorithm>
#include <numeric>
#include <stdexcept>
#include <vector>
#include <limits>
#include <cstdint>
#include <deque>
#include <queue>
#include <fwdpp/forward_types.hpp>
#include <fwdpp/simfunctions/recycling.hpp>
namespace fwdpp
{
namespace ts
{
inline flagged_mutation_queue
make_mut_queue(
const std::vector<std::uint32_t> &mcounts,
const std::vector<std::uint32_t> &counts_from_preserved_nodes)
{
flagged_mutation_queue::value_type mutation_recycling_bin;
for (std::size_t i = 0; i < mcounts.size(); ++i)
{
if (mcounts[i] + counts_from_preserved_nodes[i] == 0)
{
mutation_recycling_bin.push(i);
}
}
return flagged_mutation_queue(std::move(mutation_recycling_bin));
}
inline flagged_mutation_queue
make_mut_queue(
const std::vector<std::size_t> &preserved_mutation_indexes,
const std::size_t num_mutations)
{
std::vector<std::size_t> mindexes(num_mutations);
std::iota(mindexes.begin(), mindexes.end(), 0);
for (auto i : preserved_mutation_indexes)
{
mindexes[i] = std::numeric_limits<std::size_t>::max();
}
flagged_mutation_queue::value_type rv;
for (auto i : mindexes)
{
if (i != std::numeric_limits<std::size_t>::max())
{
rv.push(i);
}
}
return flagged_mutation_queue(std::move(rv));
}
namespace detail
{
template <typename mcont_t, typename lookup_table>
void
process_mutation_index(mcont_t &mutations, lookup_table &lookup,
const std::size_t i)
{
auto itr = lookup.equal_range(mutations[i].pos);
mutations[i].pos = std::numeric_limits<double>::max();
while (itr.first != itr.second)
{
if (itr.first->second == i)
{
lookup.erase(itr.first);
break;
}
++itr.first;
}
}
template <typename mcont_t, typename mutation_count_container,
typename lookup_table>
inline void
process_fixations(mcont_t &mutations,
mutation_count_container &mcounts,
mcont_t & /*fixations*/,
std::vector<uint_t> & /*fixation_times*/,
lookup_table &mut_lookup,
const uint_t /*generation*/, const std::size_t i,
std::false_type, std::false_type)
// Remove all fixations.
// Do not record fixations
{
process_mutation_index(mutations, mut_lookup, i);
mcounts[i] = 0;
}
template <typename mcont_t, typename mutation_count_container,
typename lookup_table>
inline void
process_fixations(mcont_t &mutations,
mutation_count_container &mcounts,
mcont_t &fixations,
std::vector<uint_t> &fixation_times,
lookup_table &mut_lookup,
const uint_t generation, const std::size_t i,
std::false_type, std::true_type)
// Remove all fixations
// Record fixations
{
fixations.push_back(mutations[i]);
fixation_times.push_back(generation);
process_mutation_index(mutations, mut_lookup, i);
mcounts[i] = 0;
}
struct equal_range_comparison
{
template <typename mutation_type>
inline bool
operator()(const mutation_type &m, const double p) const
{
return m.pos < p;
}
template <typename mutation_type>
inline bool
operator()(const double p, const mutation_type &m) const
{
return p < m.pos;
}
};
template <typename mcont_t, typename mutation_count_container,
typename lookup_table>
inline void
process_fixations(mcont_t &mutations,
mutation_count_container &mcounts,
mcont_t & /*fixations*/,
std::vector<uint_t> & /*fixation_times*/,
lookup_table &mut_lookup,
const uint_t /*generation*/, const std::size_t i,
std::true_type, std::false_type)
// Only remove neutral fixations
// Do not record fixations
{
if (mutations[i].neutral)
{
process_mutation_index(mutations, mut_lookup, i);
mcounts[i] = 0;
}
}
template <typename mcont_t, typename mutation_count_container,
typename lookup_table>
inline void
process_fixations(mcont_t &mutations,
mutation_count_container &mcounts,
mcont_t &fixations,
std::vector<uint_t> &fixation_times,
lookup_table &mut_lookup,
const uint_t generation, const std::size_t i,
std::true_type, std::true_type)
// Only remove neutral fixations.
// Record fixations.
// The recording is in order of mutation position.
{
if (mutations[i].neutral)
{
// Record the fixation here,
// as we will recycle this variant
auto loc = std::lower_bound(
fixations.begin(), fixations.end(),
mutations[i].pos,
[](const typename mcont_t::value_type &m,
const double val) { return m.pos < val; });
auto d = std::distance(fixations.begin(), loc);
fixations.insert(loc, mutations[i]);
fixation_times.insert(fixation_times.begin() + d,
generation);
process_mutation_index(mutations, mut_lookup, i);
mcounts[i] = 0;
}
else
{
// The logic here is tougher. We are preserving
// selected mutations, and thus we only
// want to record their fixations once. Thus,
// we have to make sure that this fixation
// hasn't been recorded at all.
auto loc_range = std::equal_range(
fixations.begin(), fixations.end(),
mutations[i].pos,
detail::equal_range_comparison());
if (std::find(loc_range.first, loc_range.second,
mutations[i])
== loc_range.second)
{
auto d = std::distance(fixations.begin(),
loc_range.first);
fixations.insert(loc_range.first,
mutations[i]);
fixation_times.insert(
fixation_times.begin() + d, generation);
}
}
}
} // namespace detail
template <typename poptype, typename mutation_count_container,
typename preserve_selected_fixations,
typename record_fixations>
void
flag_mutations_for_recycling(
poptype &pop,
mutation_count_container &mcounts_from_preserved_nodes,
const uint_t twoN, const uint_t generation,
const preserve_selected_fixations preserve,
const record_fixations record)
{
for (std::size_t i = 0; i < pop.mcounts.size(); ++i)
{
if (mcounts_from_preserved_nodes[i] == 0)
{
if (pop.mcounts[i] > twoN)
{
throw std::runtime_error(
"mutation count out of range");
}
if (pop.mcounts[i] == twoN)
{
detail::process_fixations(
pop.mutations, pop.mcounts,
pop.fixations, pop.fixation_times,
pop.mut_lookup, generation, i,
preserve, record);
}
else if (pop.mcounts[i] == 0)
{
detail::process_mutation_index(
pop.mutations, pop.mut_lookup, i);
}
}
}
}
} // namespace ts
} // namespace fwdpp
#endif