Program Listing for File fitness_models.hpp¶
↰ Return to documentation for file (fwdpp/fitness_models.hpp)
#ifndef _FITNESS_MODELS_HPP_
#define _FITNESS_MODELS_HPP_
#include <fwdpp/forward_types.hpp>
#include <fwdpp/fwd_functional.hpp>
#include <fwdpp/type_traits.hpp>
#include <fwdpp/util/named_type.hpp>
#include <cmath>
#include <stdexcept>
#include <type_traits>
#include <algorithm>
#include <functional>
/*
Revision history
KRT June 9 2016:
1. Added range-based operator() to additive/multiplicative models.
1a. I have not (yet) made these new operator() called by the old ones taking
haploid_genomes as args.
2. Updated documentation for site_dependent_genetic_value
3. Changed typedef for result_type in additive/multiplicative models to refer
to value from site_dependent_genetic_value
KRT June 8 2016:
1. Added range-based operator() to site_dependent_genetic_value.
2. Changed site_dependent_genetic_value::operator() that took haploid_genomes as
arguments
to call the new range-based function.
3. Fixed site_dependent_genetic_value::operator() for custom diploids.
This had never been updated from the iterator-based library to the
index-based library,
but it was never caught b/c most calls go through the specific
multiplicative/additive models.
KRT April 24 2017
1. additive_diploid and multiplicative_diploid changed to
accept closures allowing them to calculate "fitness" or "trait value"
more naturally. This is in response to #49 on github. This change
does not break API compatibility. The default is still to calculate
"fitness".
*/
namespace fwdpp
{
struct no_selection
{
using result_type = double;
template <typename haploid_genome_type, typename mcont_t>
inline result_type
operator()(const haploid_genome_type &, const haploid_genome_type &,
const mcont_t &) const noexcept
{
static_assert(
traits::is_haploid_genome<haploid_genome_type>::value,
"haploid_genome_type::value_type must be a haploid_genome "
"type");
static_assert(
traits::is_mutation<typename mcont_t::value_type>::value,
"mcont_t::value_type must be a mutation type");
return 1.;
}
template <typename T>
inline result_type
operator()(const T &) const noexcept
{
return 1.;
}
};
struct site_dependent_genetic_value
{
using result_type = double;
template <typename iterator_t, typename mcont_t,
typename updating_policy_hom, typename updating_policy_het,
typename make_return_value>
inline result_type
operator()(iterator_t first1, iterator_t last1, iterator_t first2,
iterator_t last2, const mcont_t &mutations,
const updating_policy_hom &fpol_hom,
const updating_policy_het &fpol_het,
const make_return_value & rv_function,
const double starting_value) const noexcept
{
static_assert(
traits::is_mutation<typename mcont_t::value_type>::value,
"mcont_t::value_type must be a mutation type");
static_assert(
std::is_convertible<
updating_policy_hom,
std::function<void(double &, const typename mcont_t::
value_type &)>>::value,
"decltype(fpol_hom) must be convertible to "
"std::function<void(double &,const typename "
"mcont_t::value_type");
static_assert(
std::is_convertible<
updating_policy_het,
std::function<void(double &, const typename mcont_t::
value_type &)>>::value,
"decltype(fpol_het) must be convertible to "
"std::function<void(double &,const typename "
"mcont_t::value_type");
result_type w = starting_value;
if (first1 == last1 && first2 == last2)
return rv_function(w);
else if (first1 == last1)
{
for (; first2 != last2; ++first2)
fpol_het(w, mutations[*first2]);
return rv_function(w);
}
else if (first2 == last2)
{
for (; first1 != last1; ++first1)
fpol_het(w, mutations[*first1]);
return rv_function(w);
}
for (; first1 != last1; ++first1)
{
for (; first2 != last2 && *first1 != *first2
&& mutations[*first2].pos < mutations[*first1].pos;
++first2)
// All mutations in this range are Aa
{
fpol_het(w, mutations[*first2]);
}
if (first2 < last2
&& (*first1 == *first2
|| mutations[*first1].pos
== mutations[*first2].pos))
// mutation with index first1 is homozygous
{
fpol_hom(w, mutations[*first1]);
++first2; // increment so that we don't re-process
// this site as a het next time 'round
}
else // mutation first1 is heterozygous
{
fpol_het(w, mutations[*first1]);
}
}
for (; first2 != last2; ++first2)
{
fpol_het(w, mutations[*first2]);
}
return rv_function(w);
}
template <typename iterator_t, typename mcont_t,
typename updating_policy_hom, typename updating_policy_het>
inline result_type
operator()(iterator_t first1, iterator_t last1, iterator_t first2,
iterator_t last2, const mcont_t &mutations,
const updating_policy_hom &fpol_hom,
const updating_policy_het &fpol_het,
const double starting_value) const noexcept
{
return this->operator()(first1, last1, first2, last2, mutations, fpol_hom, fpol_het,
[](double d) { return d; }, starting_value);
}
template <typename haploid_genome_type, typename mcont_t,
typename updating_policy_hom, typename updating_policy_het,
typename make_return_value>
inline result_type
operator()(const haploid_genome_type &g1,
const haploid_genome_type &g2, const mcont_t &mutations,
const updating_policy_hom &fpol_hom,
const updating_policy_het &fpol_het,
const make_return_value & rv_function,
const double starting_value) const noexcept
{
static_assert(
traits::is_haploid_genome<haploid_genome_type>::value,
"haploid_genome_type::value_type must be a haploid_genome "
"type");
return this->operator()(
g1.smutations.cbegin(), g1.smutations.cend(),
g2.smutations.cbegin(), g2.smutations.cend(), mutations,
fpol_hom, fpol_het, rv_function, starting_value);
}
template <typename haploid_genome_type, typename mcont_t,
typename updating_policy_hom, typename updating_policy_het>
inline result_type
operator()(const haploid_genome_type &g1,
const haploid_genome_type &g2, const mcont_t &mutations,
const updating_policy_hom &fpol_hom,
const updating_policy_het &fpol_het,
const double starting_value) const noexcept
{
return this->operator()(g1,g2,mutations,fpol_hom, fpol_het,[](double d){return d;}, starting_value);
}
template <typename diploid, typename gcont_t, typename mcont_t,
typename updating_policy_hom, typename updating_policy_het,
typename make_return_value>
inline result_type
operator()(const diploid &dip, const gcont_t &haploid_genomes,
const mcont_t &mutations,
const updating_policy_hom &fpol_hom,
const updating_policy_het &fpol_het,
const make_return_value & rv_function,
const double starting_value) const noexcept
{
return this->operator()(haploid_genomes[dip.first],
haploid_genomes[dip.second], mutations,
fpol_hom, fpol_het, rv_function, starting_value);
}
template <typename diploid, typename gcont_t, typename mcont_t,
typename updating_policy_hom, typename updating_policy_het>
inline result_type
operator()(const diploid &dip, const gcont_t &haploid_genomes,
const mcont_t &mutations,
const updating_policy_hom &fpol_hom,
const updating_policy_het &fpol_het,
const double starting_value) const noexcept
{
return this->operator()(dip, haploid_genomes, mutations, fpol_hom,
fpol_het, [](double d) { return d; }, starting_value);
}
};
using site_dependent_fitness = site_dependent_genetic_value;
struct haplotype_dependent_trait_value
{
using result_type = double;
template <typename haploid_genome_type, typename mcont_t,
typename haplotype_policy, typename diploid_policy>
inline result_type
operator()(const haploid_genome_type &g1,
const haploid_genome_type &g2, const mcont_t &mutations,
const haplotype_policy &hpol,
const diploid_policy &dpol) const noexcept
{
static_assert(traits::is_haploid_genome_v<haploid_genome_type>,
"haploid_genome_type must be a haploid_genome type");
static_assert(
traits::is_mutation<typename mcont_t::value_type>::value,
"mcont_t::value_type must be a mutation type");
return dpol(hpol(g1, mutations), hpol(g2, mutations));
}
template <typename diploid_t, typename gcont_t, typename mcont_t,
typename haplotype_policy, typename diploid_policy>
inline result_type
operator()(const diploid_t &diploid, const gcont_t &haploid_genomes,
const mcont_t &mutations, const haplotype_policy &hpol,
const diploid_policy &dpol) const noexcept
{
static_assert(traits::is_diploid<diploid_t>::value,
"diploid_t must represent a diploid");
return this->operator()(haploid_genomes[diploid.first],
haploid_genomes[diploid.second], mutations,
hpol, dpol);
}
};
using haplotype_dependent_fitness = haplotype_dependent_trait_value;
struct genetic_value_is_trait
{
};
struct genetic_value_is_fitness
{
};
using trait = strong_types::named_type<double, genetic_value_is_trait>;
using fitness = strong_types::named_type<double, genetic_value_is_fitness>;
inline bool
assign_is_trait_value(const trait &)
{
return true;
}
inline bool
assign_is_trait_value(const fitness &)
{
return false;
}
inline bool
assign_is_fitness_value(const trait &)
{
return false;
}
inline bool
assign_is_fitness_value(const fitness &)
{
return true;
}
struct multiplicative_diploid
{
std::function<double(double)>
assign_f(trait &)
{
return [](const double d) { return d - 1.0; };
}
std::function<double(double)>
assign_f(fitness &)
{
return [](const double d) { return std::max(0.0, d); };
}
const double scaling;
const bool gvalue_is_trait;
const bool gvalue_is_fitness;
using result_type = site_dependent_genetic_value::result_type;
const std::function<double(double)> make_return_value;
multiplicative_diploid(trait t)
: scaling{ t.get() }, gvalue_is_trait(assign_is_trait_value(t)),
gvalue_is_fitness(assign_is_fitness_value(t)), make_return_value{
assign_f(t)
}
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
template<typename make_return_value_fxn>
multiplicative_diploid(fitness f, make_return_value_fxn&& make_rv)
: scaling{ f.get() }, gvalue_is_trait{ assign_is_trait_value(f) },
gvalue_is_fitness{ assign_is_fitness_value(f) },
make_return_value{ std::forward<make_return_value_fxn>(make_rv) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
multiplicative_diploid(fitness gvtype)
: scaling{ gvtype.get() },
gvalue_is_trait(assign_is_trait_value(gvtype)),
gvalue_is_fitness(assign_is_fitness_value(gvtype)),
make_return_value{ assign_f(gvtype) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
template<typename make_return_value_fxn>
multiplicative_diploid(trait t, make_return_value_fxn&& make_rv)
: scaling{ t.get() }, gvalue_is_trait{ assign_is_trait_value(t) },
gvalue_is_fitness{ assign_is_fitness_value(t) },
make_return_value{ std::forward<make_return_value_fxn>(make_rv) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
template <typename iterator_t, typename mcont_t>
inline result_type
operator()(iterator_t first1, iterator_t last1, iterator_t first2,
iterator_t last2, const mcont_t &mutations) const noexcept
{
using __mtype = typename mcont_t::value_type;
return site_dependent_genetic_value()(
first1, last1, first2, last2, mutations,
[this](double &value, const __mtype &mut) noexcept {
value *= (1. + scaling * mut.s);
},
[](double &value, const __mtype &mut) noexcept {
value *= (1. + mut.h * mut.s);
},
make_return_value,
1.);
}
template <typename haploid_genome_type, typename mcont_t>
inline double
operator()(const haploid_genome_type &g1,
const haploid_genome_type &g2,
const mcont_t &mutations) const noexcept
{
using __mtype = typename mcont_t::value_type;
return site_dependent_genetic_value()(
g1, g2, mutations,
[this](double &value, const __mtype &mut) noexcept {
value *= (1. + scaling * mut.s);
},
[](double &value, const __mtype &mut) noexcept {
value *= (1. + mut.h * mut.s);
},
make_return_value,
1.);
}
template <typename diploid, typename gcont_t, typename mcont_t>
inline result_type
operator()(const diploid &dip, const gcont_t &haploid_genomes,
const mcont_t &mutations) const noexcept
{
return this->operator()(haploid_genomes[dip.first],
haploid_genomes[dip.second], mutations);
}
};
struct additive_diploid
{
std::function<double(double)> assign_f(trait)
{
return [](const double d) { return d; };
}
std::function<double(double)> assign_f(fitness)
{
return [](const double d) { return std::max(0.0, 1.0 + d); };
}
const double scaling;
const bool gvalue_is_trait;
const bool gvalue_is_fitness;
const std::function<double(double)> make_return_value;
using result_type = site_dependent_genetic_value::result_type;
additive_diploid(fitness f)
: scaling{ f.get() }, gvalue_is_trait{ assign_is_trait_value(f) },
gvalue_is_fitness{ assign_is_fitness_value(f) },
make_return_value{ assign_f(f) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
template<typename make_return_value_fxn>
additive_diploid(fitness f, make_return_value_fxn&& make_rv)
: scaling{ f.get() }, gvalue_is_trait{ assign_is_trait_value(f) },
gvalue_is_fitness{ assign_is_fitness_value(f) },
make_return_value{ std::forward<make_return_value_fxn>(make_rv) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
additive_diploid(trait t)
: scaling{ t.get() }, gvalue_is_trait{ assign_is_trait_value(t) },
gvalue_is_fitness{ assign_is_fitness_value(t) },
make_return_value{ assign_f(t) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
template<typename make_return_value_fxn>
additive_diploid(trait t, make_return_value_fxn&& make_rv)
: scaling{ t.get() }, gvalue_is_trait{ assign_is_trait_value(t) },
gvalue_is_fitness{ assign_is_fitness_value(t) },
make_return_value{ std::forward<make_return_value_fxn>(make_rv) }
{
if (!std::isfinite(scaling))
{
throw std::invalid_argument(
"scaling parameter must be finite");
}
}
template <typename iterator_t, typename mcont_t>
inline result_type
operator()(iterator_t first1, iterator_t last1, iterator_t first2,
iterator_t last2, const mcont_t &mutations) const noexcept
{
using __mtype = typename mcont_t::value_type;
return site_dependent_genetic_value()(
first1, last1, first2, last2, mutations,
[this](double &value, const __mtype &mut) noexcept {
value += (scaling * mut.s);
},
[](double &value, const __mtype &mut) noexcept {
value += (mut.h * mut.s);
},
make_return_value,
0.);
}
template <typename haploid_genome_type, typename mcont_t>
inline result_type
operator()(const haploid_genome_type &g1,
const haploid_genome_type &g2,
const mcont_t &mutations) const noexcept
{
using __mtype = typename mcont_t::value_type;
return site_dependent_genetic_value()(
g1, g2, mutations,
[this](double &value, const __mtype &mut) noexcept {
value += (scaling * mut.s);
},
[](double &value, const __mtype &mut) noexcept {
value += (mut.h * mut.s);
},
make_return_value,
0.);
}
template <typename diploid, typename gcont_t, typename mcont_t>
inline result_type
operator()(const diploid &dip, const gcont_t &haploid_genomes,
const mcont_t &mutations) const noexcept
{
return this->operator()(haploid_genomes[dip.first],
haploid_genomes[dip.second], mutations);
}
};
} // namespace fwdpp
#endif /* _FITNESS_MODELS_HPP_ */