Building genetic maps

#include <vector>
#include <memory>
#include <fwdpp/GSLrng_t.hpp>
#include <fwdpp/recbinder.hpp>
#include <fwdpp/genetic_map/poisson_interval.hpp>
#include <fwdpp/genetic_map/binomial_point.hpp>
#include <fwdpp/genetic_map/genetic_map.hpp>

/* Below, we have two methods to create a fwdpp::genetic_map.
 * Both create the same map.  Each map has 3 "regions":
 * 1. A region that generates a Poisson number of breakpoints on the continuous interval [0, 10)
 * 2. A region that generates a breakpoint at position 10 with probability 0.5.
 * 3. A region that generates a Poisson number of breakpoints on the continuous interval [10, 20)
 *
 * This model is therefore two regions separated by 50 centiMorgans.
 *
 * methods 1 and 2 are the more efficient of the three, requiring fewer temporary objects.
 * method 3 is useful for situations where objects must be accessed directly.  For example,
 * when constructing Python interfaces via pybind11, one cannot access a unique_ptr.
 * It is unlikely that any method would be a measurable performance bottleneck.
 */

fwdpp::genetic_map
method1()
{
    std::vector<std::unique_ptr<fwdpp::genetic_map_unit>> callbacks;
    callbacks.emplace_back(std::make_unique<fwdpp::poisson_interval>(0, 10, 1e-3));
    callbacks.emplace_back(std::make_unique<fwdpp::binomial_point>(10., 0.5));
    callbacks.emplace_back(std::make_unique<fwdpp::poisson_interval>(10., 20, 1e-3));
    return fwdpp::genetic_map(std::move(callbacks));
}

fwdpp::genetic_map
method2()
{
    fwdpp::genetic_map gmap;
    gmap.add_callback(std::make_unique<fwdpp::poisson_interval>(0, 10, 1e-3));
    gmap.add_callback(std::make_unique<fwdpp::binomial_point>(10., 0.5));
    gmap.add_callback(std::make_unique<fwdpp::poisson_interval>(10., 20, 1e-3));
    return gmap;
}

fwdpp::genetic_map
method3()
{
    fwdpp::genetic_map gmap;
    gmap.add_callback(fwdpp::poisson_interval(0, 10, 1e-3));
    gmap.add_callback(fwdpp::binomial_point(10., 0.5));
    gmap.add_callback(fwdpp::poisson_interval(10., 20, 1e-3));
    return gmap;
}

int
main(int /*argc*/, char** /*argv*/)
{
    auto map_method1 = method1();
    auto map_method2 = method2();
    auto map_method3 = method3();

    fwdpp::GSLrng_mt rng(42);

    // A simulation expects a genetic map to be equivalent to
    // std::function<std::vector<double>(void)>.
    // Our variables map_method1 and map_method2 are not.
    // The function fwdpp::recbinder rebinds them so that
    // they conform to API requirements.

    auto bound_map1 = fwdpp::recbinder(std::cref(map_method1), rng.get());
    auto bound_map2 = fwdpp::recbinder(std::cref(map_method2), rng.get());
    auto bound_map3 = fwdpp::recbinder(std::cref(map_method3), rng.get());

    // Generate breakpoints from our model.
    // You'd use "auto" here, but we write the exact return
    // type for documentation purposes.
    std::vector<double> breakpoints = bound_map1();
    breakpoints = bound_map2();
    breakpoints = bound_map3();
}
group genetic_map_types

These types define a flexible API for building genetics maps for simulations.

fwdpp::genetic_map_unit defines an interface class that is the basis for a genetic map.

Multiple genetic map units may be collected in a fwdpp::genetic_map.

The library provides several subclasses of fwdpp::genetic_map_unit that cover many common use cases.

Typedefs

typedef binomial_interval_t<double> binomial_interval

A fwdpp::binomial_interval_t using double to represent breakpoints.

typedef binomial_point_t<double> binomial_point

A fwdpp::binomial_point_t using double to represent breakpoints.

typedef fixed_number_crossovers_t<double> fixed_number_crossovers

A fwdpp::fixed_number_crossovers_t using double to represent breakpoints.

typedef poisson_interval_t<double> poisson_interval

A fwdpp::poisson_interval_t using double to represent breakpoints.

typedef poisson_point_t<double> poisson_point

A poisson_point_t using a double to represent breakpoints.

template<typename T>
struct fwdpp::binomial_interval_t : public fwdpp::genetic_map_unit
#include <fwdpp/genetic_map/binomial_interval.hpp>

Generate a breakpoint in an interval with a given probability.

Public Functions

binomial_interval_t(T b, T e, double p)

Note

The interval is half-open on [b, e).

Parameters
  • b: Beginning of interval

  • e: End of interval

  • p: Probability of a breakpoint

void operator()(const gsl_rng *r, std::vector<double> &breakpoints) const final

Generate breakpints

Parameters
  • r: Random number generator

  • breakpoints: container to fill with breakpoints

std::unique_ptr<genetic_map_unit> clone() const final

Clone object

Return

std::unique_ptr to genetic_map_unit.

Public Members

double beg

Beginning of interval (inclusive)

double end

End of interval (exclusive)

double prob

Probability of a breakpoint.

template<typename T>
struct fwdpp::binomial_point_t : public fwdpp::genetic_map_unit
#include <fwdpp/genetic_map//binomial_point.hpp>

Generate a breakpoint at a specific position with a given probability.

Public Functions

binomial_point_t(const T pos, const double pr)

Parameters
  • pos: The breakpoint position

  • pr: The probability of a breakpoint

void operator()(const gsl_rng *r, std::vector<double> &breakpoints) const final

Note

Future revisions may change the return type to void and allow for a reusable vector.

std::unique_ptr<genetic_map_unit> clone() const final

Clone object

Return

std::unique_ptr to genetic_map_unit.

Public Members

T position

Breakpoint position.

double prob

Breakpoint probability.

template<typename T>
struct fwdpp::fixed_number_crossovers_t : public fwdpp::genetic_map_unit
#include <fwdpp/genetic_map/fixed_number_crossovers.hpp>

Generate a specific number of breakpoints in an interval.

Public Functions

fixed_number_crossovers_t(T b, T e, int n)

Note

The interval is half-open on [b, e).

Parameters
  • b: Beginning of interval

  • e: End of interval

  • n: Number of braekpoints to generate

void operator()(const gsl_rng *r, std::vector<double> &breakpoints) const final

Note

Future revisions may change the return type to void and allow for a reusable vector.

std::unique_ptr<genetic_map_unit> clone() const final

Clone object

Return

std::unique_ptr to genetic_map_unit.

Public Members

const double beg

Beginning of range (inclusive)

const double end

End of range (exclusive)

const int nxovers

Number of crossovers.

class fwdpp::genetic_map
#include <fwdpp/genetic_map/genetic_map.hpp>

A genetic map is a container of fwdpp::genetic_map_unit objects.

Public Functions

genetic_map(std::vector<std::unique_ptr<genetic_map_unit>> &&c)

Constructor

Parameters

void add_callback(std::unique_ptr<genetic_map_unit> &&gu)

Add a new callback by moving the input.

void add_callback(const genetic_map_unit &gu)

Add a new callback by cloning the input.

std::vector<double> operator()(const gsl_rng *r) const

Note

Future revisions may change the return type to void and allow for a reusable vector.

std::size_t size() const

Return number of stored callbacks.

struct fwdpp::genetic_map_unit : public fwdpp::util::abstract_cloneable<genetic_map_unit>
#include <fwdpp/genetic_map/genetic_map_unit.hpp>

Interface class.

Subclassed by fwdpp::binomial_interval_t< T >, fwdpp::binomial_point_t< T >, fwdpp::fixed_number_crossovers_t< T >, fwdpp::poisson_interval_t< T >, fwdpp::poisson_point_t< T >

Public Functions

genetic_map_unit()

Constructor.

void operator()(const gsl_rng*, std::vector<double>&) const = 0

Note

Future revisions may change the return type to void and allow for a reusable vector.

template<typename T>
struct fwdpp::poisson_interval_t : public fwdpp::genetic_map_unit
#include <fwdpp/genetic_map/poisson_interval.hpp>

Generate a Poisson number of breakpoints in an interval.

Public Functions

poisson_interval_t(T b, T e, double m)

Note

The interval is half-open on [b, e).

Parameters
  • b: Beginning of interval

  • e: End of interval

  • m: Mean number of breakpoints

void operator()(const gsl_rng *r, std::vector<double> &breakpoints) const final

Note

Future revisions may change the return type to void and allow for a reusable vector.

std::unique_ptr<genetic_map_unit> clone() const final

Clone object

Return

std::unique_ptr to genetic_map_unit.

Public Members

const double beg

Beginning of interval (inclusive)

const double end

End of interval (exclusive)

const double mean

Mean number of breakpoints.

template<typename T>
struct fwdpp::poisson_point_t : public fwdpp::genetic_map_unit
#include <fwdpp/genetic_map/poisson_point.hpp>

Generate a breakpoint at a position if a Poisson deviate is odd.

Public Functions

poisson_point_t(const T pos, const double m)

Parameters
  • pos: Breakpoint position

  • m: The mean of a Poisson distribution

void operator()(const gsl_rng *r, std::vector<double> &breakpoints) const final

Note

Future revisions may change the return type to void and allow for a reusable vector.

std::unique_ptr<genetic_map_unit> clone() const final

Clone object

Return

std::unique_ptr to genetic_map_unit.

Public Members

const T position

Breakpoint position.

double mean

Mean number of crossovers at position.