Program Listing for File mutate_tables.hpp

↰ Return to documentation for file (fwdpp/ts/mutate_tables.hpp)

#ifndef FWDPP_TS_MUTATE_TABLES_HPP
#define FWDPP_TS_MUTATE_TABLES_HPP

#include <gsl/gsl_randist.h>

#include <vector>
#include "definitions.hpp"
#include "mark_multiple_roots.hpp"
#include "mutation_tools.hpp"
#include "table_collection_functions.hpp"

namespace fwdpp
{
    namespace ts
    {
        template <typename TableCollectionType, typename Samples, typename rng,
                  typename mfunction>
        unsigned
        mutate_tables(const rng &r, const mfunction &make_mutation,
                      TableCollectionType &tables, Samples &&samples, const double mu)
        {
            unsigned nmuts = 0;
            if (!(mu > 0.0))
                {
                    return nmuts;
                }
            auto mr = mark_multiple_roots(tables, std::forward<Samples>(samples));
            const double L = tables.genome_length();
            for (auto &i : mr)
                {
                    auto dt = tables.nodes[i.first].time;
                    for (auto j : i.second)
                        {
                            double mean = dt * (j.second - j.first) * mu / L;
                            auto nm = gsl_ran_poisson(r.get(), mean);
                            nmuts += nm;
                            for (unsigned m = 0; m < nm; ++m)
                                {
                                    unsigned g = static_cast<unsigned>(
                                        gsl_ran_flat(r.get(), 1, dt + 1));
                                    new_variant_record r
                                        = make_mutation(j.first, j.second, g);
                                    auto newsite = tables.emplace_back_site(
                                        r.s.position, r.s.ancestral_state);
                                    tables.emplace_back_mutation(i.first, r.key, newsite,
                                                                 r.derived_state,
                                                                 r.neutral);
                                }
                        }
                }
            for (auto &e : tables.edges)
                {
                    auto ct = tables.nodes[e.child].time;
                    auto pt = tables.nodes[e.parent].time;
                    auto dt = ct - pt;
                    double mean = dt * (e.right - e.left) * mu / L;
                    auto nm = gsl_ran_poisson(r.get(), mean);
                    for (unsigned m = 0; m < nm; ++m)
                        {
                            unsigned g = static_cast<unsigned>(
                                gsl_ran_flat(r.get(), pt + 1, ct + 1));
                            new_variant_record r = make_mutation(e.left, e.right, g);
                            auto site = tables.emplace_back_site(r.s.position,
                                                                 r.s.ancestral_state);
                            tables.emplace_back_mutation(
                                e.child, r.key, site, r.derived_state, r.neutral);
                        }
                    nmuts += nm;
                }
            sort_mutation_table_and_rebuild_site_table(tables);
            return nmuts;
        }
    } // namespace ts
} // namespace fwdpp

#endif