Program Listing for File haploid_genome_cleaner.hpp

↰ Return to documentation for file (fwdpp/internal/haploid_genome_cleaner.hpp)

#ifndef __FWDPP_INTERNAL_HAPLOID_GENOME_CLEANER_HPP__
#define __FWDPP_INTERNAL_HAPLOID_GENOME_CLEANER_HPP__

#include <vector>
#include <algorithm>
#include <limits>
#include <utility>
#include <type_traits>
#include <fwdpp/forward_types.hpp>
#include <fwdpp/fwd_functional.hpp>

namespace fwdpp
{
    namespace fwdpp_internal
    {

        // First, we have a set of helper functions:

        template <typename GenomeContainerType_iterator>
        inline GenomeContainerType_iterator
        next_extant_haploid_genome(GenomeContainerType_iterator beg,
                                   GenomeContainerType_iterator end) noexcept
        {
            return std::find_if(
                beg, end,
                [](const typename GenomeContainerType_iterator::value_type &g) noexcept {
                    return g.n;
                });
        }

        struct find_fixation
        {
            template <typename mut_index_cont, typename mcounts_t>
            inline typename mut_index_cont::const_iterator
            operator()(const mut_index_cont &mc, const mcounts_t &mcounts,
                       const uint_t twoN) const noexcept
            {
                return std::find_if(mc.cbegin(), mc.cend(),
                                    [&mcounts, &twoN](const std::size_t &i) noexcept {
                                        return mcounts[i] == twoN;
                                    });
            }

            template <typename mut_index_cont, typename MutationContainerType,
                      typename mcounts_t, typename mutation_removal_policy>
            inline typename mut_index_cont::const_iterator
            operator()(const mut_index_cont &mc, const MutationContainerType &mutations,
                       const mcounts_t &mcounts, const uint_t twoN,
                       mutation_removal_policy &&mp) const noexcept
            {
                return std::find_if(
                    mc.cbegin(), mc.cend(),
                    [&mutations, &mcounts, &twoN, &mp](const std::size_t &i) noexcept {
                        return mcounts[i] == twoN && mp(mutations[i]);
                    });
            }
        };

        struct haploid_genome_cleaner_erase_remove_idiom_wrapper
        {
            template <typename mut_index_cont, typename mcounts_t>
            inline void
            operator()(mut_index_cont &mc, const mcounts_t &mcounts,
                       const typename mut_index_cont::value_type &first_fixation,
                       const uint_t twoN) const noexcept
            {
                /*
                  The first call to std::find relies on a de-referencing of the
                  return value of
                  find_fixation, passed to this function as first_fixation.
                  Because mc is sorted according to mutation position,
                  first_fixation is the fixation with the smallest position.
                  Using std::find like this avoids some calls to the lambda
                  expression, which experiences cache-misses because mcounts is
                  NOT sorted according to mutation position.
                */
                mc.erase(std::remove_if(
                             std::find(mc.begin(), mc.end(), first_fixation), mc.end(),
                             [&mcounts, &twoN](
                                 const typename mut_index_cont::value_type &i) noexcept {
                                 return mcounts[i] == twoN;
                             }),
                         mc.end());
            }

            template <typename mut_index_cont, typename MutationContainerType,
                      typename mcounts_t, typename mutation_removal_policy>
            inline void
            operator()(mut_index_cont &mc, const MutationContainerType &mutations,
                       const mcounts_t &mcounts,
                       const typename mut_index_cont::value_type &first_fixation,
                       const uint_t twoN, mutation_removal_policy &&mp) const noexcept
            {
                /*
                  The first call to std::find relies on a de-referencing of the
                  return value of
                  find_fixation, passed to this function as first_fixation.
                  Because mc is sorted according to mutation position,
                  first_fixation is the fixation with the smallest position.
                  Using std::find like this avoids some calls to the lambda
                  expression, which experiences cache-misses because mcounts is
                  NOT sorted according to mutation position.
                */
                mc.erase(std::remove_if(
                             std::find(mc.begin(), mc.end(), first_fixation), mc.end(),
                             [&mcounts, &mutations, &twoN, &mp](
                                 const typename mut_index_cont::value_type &i) noexcept {
                                 return mcounts[i] == twoN && mp(mutations[i]);
                             }),
                         mc.end());
            }

            template <typename mut_index_cont, typename mcounts_t>
            inline void
            operator()(mut_index_cont &mc, const mcounts_t &mcounts,
                       const uint_t twoN) const noexcept
            {
                /*
                  \note Added in 0.5.0 to address Issue #41
                */
                mc.erase(std::remove_if(
                             mc.begin(), mc.end(),
                             [&mcounts, &twoN](
                                 const typename mut_index_cont::value_type &i) noexcept {
                                 return mcounts[i] == twoN;
                             }),
                         mc.end());
            }

            template <typename mut_index_cont, typename MutationContainerType,
                      typename mcounts_t, typename mutation_removal_policy>
            inline void
            operator()(mut_index_cont &mc, const MutationContainerType &mutations,
                       const mcounts_t &mcounts, const uint_t twoN,
                       mutation_removal_policy &&mp) const noexcept
            {
                /*
                  \note Added in 0.5.0 to address Issue #41
                */
                mc.erase(std::remove_if(
                             mc.begin(), mc.end(),
                             [&mcounts, &mutations, &twoN, &mp](
                                 const typename mut_index_cont::value_type &i) noexcept {
                                 return mcounts[i] == twoN && mp(mutations[i]);
                             }),
                         mc.end());
            }
        };

        template <typename GenomeContainerType, typename MutationContainerType,
                  typename mutation_removal_policy>
        inline typename std::enable_if<
            std::is_same<mutation_removal_policy, fwdpp::remove_nothing>::value>::type
        haploid_genome_cleaner(GenomeContainerType &, const MutationContainerType &,
                               const std::vector<uint_t> &, const uint_t,
                               const mutation_removal_policy &)
        {
            return;
        }

        template <typename GenomeContainerType, typename fixation_finder,
                  typename idiom_wrapper>
        inline void
        haploid_genome_cleaner_details(GenomeContainerType &haploid_genomes,
                                       const fixation_finder &ff,
                                       const idiom_wrapper &iw) noexcept
        {
            auto extant_haploid_genome = next_extant_haploid_genome(
                haploid_genomes.begin(), haploid_genomes.end());
            const auto fixation_n = ff(extant_haploid_genome->mutations);
            bool neutral_fixations_exist
                = (fixation_n != extant_haploid_genome->mutations.cend());
            const auto fixation_s = ff(extant_haploid_genome->smutations);
            bool selected_fixations_exist
                = (fixation_s != extant_haploid_genome->smutations.cend());
            if (!neutral_fixations_exist && !selected_fixations_exist)
                return;

            const auto gend = haploid_genomes.end();
            // Assign values to avoid tons of de-referencing later
            const auto fixation_n_value
                = (fixation_n == extant_haploid_genome->mutations.cend())
                      ? typename decltype(fixation_n)::value_type()
                      : *fixation_n;
            const auto fixation_s_value
                = (fixation_s == extant_haploid_genome->smutations.cend())
                      ? typename decltype(fixation_s)::value_type()
                      : *fixation_s;
            while (extant_haploid_genome < gend)
                {
                    if (neutral_fixations_exist)
                        {
                            iw(extant_haploid_genome->mutations, fixation_n_value);
                        }
                    if (selected_fixations_exist)
                        {
                            iw(extant_haploid_genome->smutations, fixation_s_value);
                        }
                    extant_haploid_genome
                        = next_extant_haploid_genome(extant_haploid_genome + 1, gend);
                }
        }

        template <typename GenomeContainerType, typename fixation_finder>
        std::pair<bool, bool>
        fixation_finder_search_all(const GenomeContainerType &haploid_genomes,
                                   const fixation_finder &ff) noexcept
        {
            bool neutral = false, selected = false;
            for (const auto &g : haploid_genomes)
                {
                    if (!neutral)
                        {
                            auto itr = ff(g.mutations);
                            if (itr != g.mutations.cend())
                                {
                                    neutral = true;
                                }
                        }
                    if (!selected)
                        {
                            auto itr = ff(g.smutations);
                            if (itr != g.smutations.cend())
                                {
                                    selected = true;
                                }
                        }
                    if (neutral && selected)
                        return std::make_pair(neutral, selected);
                }
            return std::make_pair(neutral, selected);
        }

        template <typename GenomeContainerType, typename fixation_finder,
                  typename idiom_wrapper>
        inline void
        haploid_genome_cleaner_details(GenomeContainerType &haploid_genomes,
                                       const fixation_finder &ff,
                                       const idiom_wrapper &iw, std::true_type) noexcept
        {
            auto t = fixation_finder_search_all(haploid_genomes, ff);
            auto neutral_fixations_exist = t.first;
            auto selected_fixations_exist = t.second;
            if (!neutral_fixations_exist && !selected_fixations_exist)
                return;

            auto extant_haploid_genome = next_extant_haploid_genome(
                haploid_genomes.begin(), haploid_genomes.end());
            const auto gend = haploid_genomes.end();
            while (extant_haploid_genome < gend)
                {
                    if (neutral_fixations_exist)
                        {
                            iw(extant_haploid_genome->mutations);
                        }
                    if (selected_fixations_exist)
                        {
                            iw(extant_haploid_genome->smutations);
                        }
                    extant_haploid_genome
                        = next_extant_haploid_genome(extant_haploid_genome + 1, gend);
                }
        }

        /*
          Now, we have the two overloads of haploid_genome_cleaner.
          The std::cref(...) below are required, o/w copies will get bound,
          which is bad for performance.
        */

        template <typename GenomeContainerType, typename MutationContainerType,
                  typename mutation_removal_policy>
        inline typename std::enable_if<
            std::is_same<mutation_removal_policy, std::true_type>::value>::type
        haploid_genome_cleaner(GenomeContainerType &haploid_genomes,
                               const MutationContainerType &,
                               const std::vector<uint_t> &mcounts, const uint_t twoN,
                               const mutation_removal_policy &)
        {
            haploid_genome_cleaner_details(
                haploid_genomes,
                [&mcounts,
                 twoN](const typename GenomeContainerType::value_type::mutation_container
                           &mc) { return find_fixation()(mc, mcounts, twoN); },
                [&mcounts,
                 twoN](typename GenomeContainerType::value_type::mutation_container &mc,
                       const typename GenomeContainerType::value_type::
                           mutation_container::value_type v) {
                    return haploid_genome_cleaner_erase_remove_idiom_wrapper()(
                        mc, mcounts, v, twoN);
                });
        }

        template <typename GenomeContainerType, typename MutationContainerType,
                  typename mutation_removal_policy>
        inline typename std::enable_if<
            !std::is_same<mutation_removal_policy, std::true_type>::value
            && !std::is_same<mutation_removal_policy,
                             fwdpp::remove_nothing>::value>::type
        haploid_genome_cleaner(GenomeContainerType &haploid_genomes,
                               const MutationContainerType &mutations,
                               const std::vector<uint_t> &mcounts, const uint_t twoN,
                               const mutation_removal_policy &mp)
        {
            haploid_genome_cleaner_details(
                haploid_genomes,
                [&mutations, &mcounts, &mp,
                 twoN](const typename GenomeContainerType::value_type::mutation_container
                           &mc) {
                    return find_fixation()(mc, mutations, mcounts, twoN, mp);
                },
                [&mutations, &mcounts, &mp,
                 twoN](typename GenomeContainerType::value_type::mutation_container &mc,
                       typename GenomeContainerType::value_type::mutation_container::
                           value_type v) {
                    return haploid_genome_cleaner_erase_remove_idiom_wrapper()(
                        mc, mutations, mcounts, v, twoN, mp);
                });
        }

        template <typename GenomeContainerType, typename MutationContainerType,
                  typename mutation_removal_policy>
        inline typename std::enable_if<
            std::is_same<mutation_removal_policy, std::true_type>::value>::type
        haploid_genome_cleaner(GenomeContainerType &haploid_genomes,
                               const MutationContainerType &,
                               const std::vector<uint_t> &mcounts, const uint_t twoN,
                               const mutation_removal_policy &, std::true_type)
        {
            haploid_genome_cleaner_details(
                haploid_genomes,
                [&mcounts,
                 twoN](const typename GenomeContainerType::value_type::mutation_container
                           &mc) { return find_fixation()(mc, mcounts, twoN); },
                [&mcounts, twoN](
                    typename GenomeContainerType::value_type::mutation_container &mc) {
                    return haploid_genome_cleaner_erase_remove_idiom_wrapper()(
                        mc, mcounts, twoN);
                },
                std::true_type());
        }

        template <typename GenomeContainerType, typename MutationContainerType,
                  typename mutation_removal_policy>
        inline typename std::enable_if<
            !std::is_same<mutation_removal_policy, std::true_type>::value
            && !std::is_same<mutation_removal_policy,
                             fwdpp::remove_nothing>::value>::type
        haploid_genome_cleaner(GenomeContainerType &haploid_genomes,
                               const MutationContainerType &mutations,
                               const std::vector<uint_t> &mcounts, const uint_t twoN,
                               const mutation_removal_policy &mp, std::true_type)
        {
            haploid_genome_cleaner_details(
                haploid_genomes,
                [&mutations, &mcounts, &mp,
                 twoN](const typename GenomeContainerType::value_type::mutation_container
                           &mc) {
                    return find_fixation()(mc, mutations, mcounts, twoN, mp);
                },
                [&mutations, &mcounts, &mp, twoN](
                    typename GenomeContainerType::value_type::mutation_container &mc) {
                    return haploid_genome_cleaner_erase_remove_idiom_wrapper()(
                        mc, mutations, mcounts, twoN, mp);
                },
                std::true_type());
        }
    }
}

#endif