Welcome to fwdpp’s documentation!

Introduction

Note

The documentation is currently in transition from doxygen to Sphinx. This change requires a lot of re-writing. The rewriting is good because we’ll get stuff up to date. The down side is that it will take a while!

fwdpp is a C++14 template library for implementing forward-time simulations of population- and quantitative- genetic models.

Some of you may have landed here because you want to run simulations. We encourange you to check out fwdpy11, which is developed using fwdpp

Citation

If you use fwdpp in your research, or a package developed with it, then please cite the following paper:

@Article{,
  author =   {K. R. Thornton},
  title =    {A C++ Template Library for Efficient Forward-Time Population Genetic Simulation of Large Populations},
  journal =          {Genetics},
  year =     {2014},
  OPTkey =   {},
  volume =   {198},
  OPTnumber =        {},
  pages =    {157-166},
  OPTmonth =         {},
  OPTnote =          {},
  annote =   {doi:/10.1534/genetics.114.165019}
}

License and Distribution

License

fwdpp is distributed under the terms of the GNU Public License, version 3 or later. This is also referred to as GPL3+.

The license file is available with the source code via the files LICENSE and COPYING.

Distribution

fwdpp can be obtained via GitHub.

Build status

Main branch:

https://circleci.com/gh/molpopgen/fwdpp.svg?style=svg

Development branch:

https://circleci.com/gh/molpopgen/fwdpp.svg?style=svg

Compiling, installing, etc.

Obtaining the source

The source code can be obtained from GitHub. For example,

git clone https://github.com/molpopgen/fwdpp
cd fwdpp
git submodule update --init --recursive

Note

It is important to update the submodules! Some of the examples and several of the unit tests make use of tskit, which is one of the submodules.

Warning

Be careful if you download the source code from the “releases” section of GitHub. Do not download the files that GitHub automatically generates! They do not include the submodule.

Dependencies

  1. A compiler that supports C++14 or later.

  2. GNU autotools

  3. The GNU Scientific Library

  4. The boost libraries. The boost-test and boost-program-options are used for the test suite and for some examples programs, respectively.

All of these dependencies are available via package managers. On Debian/Ubuntu-like systems, you may use apt to install them. On other operating systems, conda and/or brew will do the job.

Building examples and running tests

Note

The test suite will not be compiled if the boost test library is not found on your system.

To build the examples and the test suite:

autoreconf --install
./configure
make

The compilation can take a long time. To enable parallel compilation, you may execute make -j 6 to use six threads, for example.

To run the tests:

make check

To install the library:

# This may require sudo
make install

Changing the compilation flags

The default compilation flags are -DNDEBUG -O2 -g. The -DNDEBUG specifies compiling in “release mode”. To compile in “debug” mode:

./configure --enable-debug=yes

Note

Compiling in debug mode makes simulations run a lot slower! The test suite is always compiled in debug mode.

To change where the library will be installed:

To change the optimization levels, etc.:

CXXFLAGS="-O3 -W -Wall -Wconversion" ./configure

The above example enables more aggressive optimizations and sets more warning flags.

The default C++ mode is C++14. To enable C++17:

./configure --enable-cpp17=yes

If you have dependencies installed in non-standard locations, then you may need to provide those locations to the configure script. When doing so, keep in mind that CPPFLAGS is used to change where the compiler looks for headers, as it is the variable affecting the preprocessor. The CXXFLAGS variable is for the compiler and not the preprocessor!

Introduction

This section defines the lowest-level types required by a simulation.

These types can be included all at once via File forward_types.hpp.

Mutations

Most types of simulations must have mutations in order to do anything “interesting”. fwdpp allows you to define your own mutation type.

A valid mutation type publicly inherits from fwdpp::mutation_base. This class defines a very minimal interface. It has no concept of an “effect size” for a mutation, etc.. To add the relevant concepts for your model, create a derived class.

For example:

#pragma once

#include <tuple>
#include <fwdpp/type_traits.hpp>
#include <fwdpp/fundamental_types/mutation_base.hpp>

/*
 * The simplest mutation type, adding just a selection coefficient
 * and dominance to the interface.
*/
struct mutation : public fwdpp::mutation_base
{
    // selection coefficient
    double s;
    // dominance coefficient
    double h;
    mutation(const double &position, const double &sel_coeff,
             const double &dominance = 0.5) noexcept
        : mutation_base(position, (sel_coeff == 0)), s(sel_coeff), h(dominance)
    {
    }

    bool
    operator==(const mutation &rhs) const
    {
        return std::tie(this->s, this->h) == std::tie(rhs.s, rhs.h) && is_equal(rhs);
    }
};

static_assert(fwdpp::traits::is_mutation_v<mutation>,
              "Mutation is not a valid mutation type!");

Containers of mutations

In a simulation, mutation objects must be stored in random-access containers of objects. You must not use, for example, smart pointers to fwdpp::mutation_base.

In other words, the value_type of your container must equal a type derived from the mutation base class.

For example:

using mutation_container = std::vector<mutation>;

By using such a container, each mutation is represented only once in a simulation. To represent individual genotypes, we need a method of tracking which mutations are present in genomes. See here for details.

Haploid genomes

This section builds on the information found in Mutations.

A haploid genome is the information inherited from a single parent. For example, a diploid individual contains two haploid genomes. In the case of obligate outcrossers, one genome may be labelled “maternal” and the other “paternal”.

The type fwdpp::haploid_genome represents a haploid genome in fwdpp. This type is itself a typedef of fwdpp::haploid_genome_base.

The key data members are:

The two classes of mutation keys are stored in different containers so that we can skip “neutral” mutations when calculating genetic values.

The type of the index containers is fwdpp::haploid_genome_base::mutation_container. The value_type of this container is fwdpp::uint_t.

Note

The field fwdpp::haploid_genome_base::n is not equivalent to the frequency of a genome. Rather, it is the number of occurrences of a specific genome object.

Containers of genomes

In a simulation, genome objects must be stored in random-access containers of objects.

In other words, the value_type of your container must equal a genome type;

For example:

using genome_container = std::vector<fwdpp::haploid_genome>;

Diploids

This section builds on the information in Haploid genomes.

The simplest representation of a diploid is a pair of integers referring to the simulation’s genome container. fwdpp treats std::pair with integer types as the simplest form of diploid. At compile time, the library ensures that first_type and second_type are both integer types and that both are the same type.

For example:

#include <utility>
#include <cstdint>
#include <fwdpp/type_traits.hpp>

using diploid_size_t = std::pair<std::size_t, std::size_t>;
using diploid_int32 = std::pair<std::int32_t, std::int32_t>;

int
main(int, char **)
{
    static_assert(fwdpp::traits::is_diploid_v<diploid_size_t>,
                  "diploid_size_t is not a diploid!");
    static_assert(!fwdpp::traits::is_custom_diploid_v<diploid_size_t>,
                  "diploid_size_t is a custom diploid!");
    static_assert(fwdpp::traits::is_diploid_v<diploid_int32>,
                  "diploid_int32 is not a diploid!");
    static_assert(!fwdpp::traits::is_custom_diploid_v<diploid_int32>,
                  "diploid_size_t is a custom diploid!");
}

In the above example, the static assertions tell us that our std::pair are diploids. They also tell us that they are not custom diploids..

Diploids are more than just pairs of genomes. You may want to record other meta data along with the genome indexes. To do so, you may define your own types that duck-type the std::pair interface, adding whatever features you would like. We refer to such types as custom diploid types:

struct dip_with_fitness
{
    using first_type = std::size_t;
    using second_type = std::size_t;
    std::size_t first, second;
    double fitness;

    dip_with_fitness() : first{}, second{}, fitness{}
    {
    }
    dip_with_fitness(first_type f, second_type s) : first{f}, second{s}, fitness{0}
    {
    }
    dip_with_fitness(first_type f, second_type s, double w)
        : first{f}, second{s}, fitness{w}
    {
    }
};

static_assert(fwdpp::traits::is_diploid_v<dip_with_fitness>,
              "Error: dip with fitness is not a valid diploid type");
static_assert(fwdpp::traits::is_custom_diploid_v<dip_with_fitness>,
              "Error: dip with fitness is not a valid diploid type");

Internally, fwdpp doesn’t care much about custom diploids versus std::pair. The support for custom diploids exists to give modeling flexibility. However, the library is flexible enough to let you take other paths as well:

  • You could represent meta data as a separate array of structures.

  • You could represent the meta data as a separate structure of arrays.

Containers of diploids

In a simulation, diploid objects must be stored in random-access containers of objects.

In other words, the value_type of your container must equal a diploid type;

For example:

using diploid_container = std::vector<my_diploid_type>;

Typedefs

The following typedefs are used throughout fwdpp:

Random number generators

The random number generator type is fwdpp::GSLrng_mt, which is a mersenne twister. The class constructor accepts a single, 32-bit unsigned integer. This class is effectively a unique_ptr with a custom deleter. As such, it defines a move-only type.

fwdpp::GSLrng_mt is defined as a template typedef of fwdpp::GSLrng_t. The function fwdpp::GSLrng_t::get() returns the underlying const gsl_rng *.

These types are defined in File GSLrng_t.hpp.

To use this type:

#include <fwdpp/GSLrng_t.hpp>

int main(int, char **)
{
    // Initialize with a seed.
    fwdpp::GSLrng_mt(42);
}

Note

Most or all of fwdpp is compatible with using a “bare” gsl_rng *. However, you give up exception safety by doing so.

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.

Change log

0.9.2

  • Discrete genetic map unit objects now check that they have valid lengths. #326.

0.9.1

  • Fixes a build issue on macOS. 8574a3e

  • Adds GitHub action to test macOS/C++14 #324

0.9.0

This is a pretty big release. In fact, it is too big to list all of the changes here. They are collected on GitHub under the 0.9.0 milestone.

Importantly, this release deprecates several features. Where possible, a C++14-style [[deprecated]] marker is used.

We plan to start releasing more often, so that individual releases aren’t so big.

In future releases, we plan to improve the organization of the library headers. We will also clean up the test suite to reflect the new layout. We will keep the existing headers and use the preprocessor to emit a warning at compile time.

Developers guide

Calculating code coverage

Code coverage requires lcov, which will be available via your favorite package manager.

The coverage is calculated from the test suite and is automated via a make target:

make -C testsuite coverage-local

To clean up after a coverage calculation:

make clean-local

The output will be in the directory fwdpp_coverage, which will be generated in the root directory of the source repository. To view the results, point your browser to fwdpp_coverage/index.html.

Note

Coverage calculations are a bit odd for header-only libraries like this one. The output only applies to library files that are included in the test suite. Thus, there may be files for which no coverage calculation is possible. We are working to ensure that all files do get covered.

Indices and tables