Skip to content

Commit

Permalink
molpopgen published a site update
Browse files Browse the repository at this point in the history
  • Loading branch information
molpopgen committed Jan 19, 2024
1 parent 0f3bcce commit 40d98fd
Show file tree
Hide file tree
Showing 94 changed files with 2,693 additions and 555 deletions.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file not shown.
181 changes: 181 additions & 0 deletions _sources/long_vignettes/demes_event_timings.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,181 @@
---
jupytext:
formats: md:myst
text_representation:
extension: .md
format_name: myst
kernelspec:
display_name: Python 3
language: python
name: python3
---

(demes_event_timings)=

# Timing of events in demographic models

This document provides several examples showing *when* events happen in demographic models.

We will define very simple setups that allow very simple assertions about their outcomes.
While most simulations will involve more complex setups, the advantage of simplicity here
is that we can easily validate the results because we can easily reason about the models.

These examples come from the `fwdpy11` test suite.

The key concept here is that the [demes specification](https://popsim-consortium.github.io/demes-spec-docs/main/introduction.html)
defines an epoch as a half open interval on `[present, past)`.
Our reasoning about the outcomes of data that we record during simulation involves breaking our models up into half open intervals and asking questions about them.

Why do we spend time going over this?
Imagine that you want to know the details about the ancestry of a given generation in your simulation.
For example, the ancestry of the generation born during a pulse.
One way to do that would be to preserve both the pulse generation and its parent generation as ancient samples.
(See [here](recorders_vignette).)
To do such sampling you need to know when to sample!

## Some setup

During some of our simulations, we will record the current time point
and the part of the metadata containing individual parents.

We can use a simple data class to record this information:


```{literalinclude} ../../tests/test_demes_event_timings.py
:language: python
:lines: 8-20
```

There will be no genetics happening during the simulation.
All that happens is that we apply our recorder.
We will always "burn in" the demographic model with 100
generations of matings.

```{literalinclude} ../../tests/test_demes_event_timings.py
:language: python
:lines: 24-42
```

```{code-cell} python
---
tags: ['remove-input']
---
import demes
import demesdraw
```

## Intervals when demes exist

In the following graph, demes completely replace one another over time:

```{literalinclude} ../../tests/demes_event_examples/deme_existence.yaml
:language: yaml
```

The following graphic depicts the model:

```{code-cell} python
---
tags: ['remove-input']
---
graph = demes.load("../../tests/demes_event_examples/deme_existence.yaml")
demesdraw.tubes(graph);
```

The intervals of a deme's existence are from $[i, j)$ generations ago, where $j > i$.
Therefore `deme2` exists from $[0, 25)$ generations ago, etc..

During simulation, a node's birth time and deme are recorded.
Therefore, we can export the simulation to a tree sequence and check that the node
times match up with their expected deme labels:

```{literalinclude} ../../tests/test_demes_event_timings.py
:language: python
:lines: 83-93
```

## A single pulse event

The following graph defines a two deme model with a single pulse event 10 generations ago.


```{literalinclude} ../../tests/demes_event_examples/single_pulse.yaml
:language: yaml
```

The following graphic depicts the model:


```{code-cell} python
---
tags: ['remove-input']
---
graph = demes.load("../../tests/demes_event_examples/single_pulse.yaml")
demesdraw.tubes(graph);
```

Given this model:

* The time interval $[0, 10)$ generations ago is the *post* pulse interval
(thinking forwards in time).
During this interval, all individuals in `deme0` have parents from that same deme.
Likewise for `deme1`.
* All individuals born in the pulse generation (10 generations ago) have all of their
parents drawn from `deme0`.
The reason for this is that there are no events affecting the ancestry of `deme0` and
the proportion of `deme1` ancestry that comes from `deme0` is 1.0, which is all of it.
* The time interval $[11, \infty)$ generations ago is the *pre* pulse interval
(thinking forwards in time).
During this interval, all individuals in `deme0` have parents from that same deme.
Likewise for `deme1`.

Using the reasoning about the pre and post epochs with respect to the timing of
the pulse, we can make the following assertions:

```{literalinclude} ../../tests/test_demes_event_timings.py
:language: python
:lines: 45-61
```

Let's work through the logic behind assertions in detail:

* The size of each deme is 100 diploids.
* The metadata are ordered by deme id.
Thus, the first 100 entries are the indexes parents of individuals in `deme0`, etc..
* Due to the metadata ordering, an index $< 100$ (or, $[0, 100)$) indicates
a parent from `deme0`.
An index in $[100, 200)$ indicates a parent from `deme1`.


## A multi-generation "burst" of migration

This model is a slight twist on the previous.
Rather than a single generation pulse, there are several consecutive generations of migration
in one direction:

```{literalinclude} ../../tests/demes_event_examples/burst_of_migration.yaml
:language: yaml
```
The following graphic depicts the model:

```{code-cell} python
---
tags: ['remove-input']
---
graph = demes.load("../../tests/demes_event_examples/burst_of_migration.yaml")
demesdraw.tubes(graph);
```

In this model, migration occurs during the interval $[5, 10)$ generations ago.
Because the migration is one way (`deme0` to `deme1`) and the rate is 1.0, the parents
of all offspring born in this interval have their parents sampled from `deme0`.

Outside of this interval, the parents of any individual are drawn from that
individual's deme.

The following code asserts that these intervals are what we expect:

```{literalinclude} ../../tests/test_demes_event_timings.py
:language: python
:lines: 64-80
```
39 changes: 2 additions & 37 deletions _sources/long_vignettes/tskitconvert_vignette.md
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,7 @@ pdict = {
'rates': (0.0, 1e-3, None),
'simlen': 10 * pop.N,
'prune_selected': False,
'demography': fwdpy11.ForwardDemesGraph.tubes([pop.N], 1)
}
params = fwdpy11.ModelParams(**pdict)
Expand All @@ -160,47 +161,11 @@ ts = pop.dump_tables_to_tskit(model_params=params)
assert fwdpy11.ModelParams(**eval(ts.metadata["model_params"])) == params
```

Some simulation designs make use of multiple `ModelParams` objects.
Imagine, for example, that we simulate the above model until equilibrium, and then simulate another 100 generations while adapting to a new optimum:

:::{note}
The details are omitted here:

1. This simulation requires two calls to {func}`fwdpy11.evolvets`
2. It could have been done instead with a single call to {func}`fwdpy11.evolvets`
by passing a list of {class}`fwdpy11.Optimum` to {class}`fwdpy11.GaussianStabilizingSelection`.
This approach would have required one `ModelParams` instance.
:::

```{code-cell} python
import copy
pdict2 = copy.deepcopy(pdict)
optimum = fwdpy11.Optimum(optimum=5.0, VS=10.0, when=0)
gss = fwdpy11.GaussianStabilizingSelection.single_trait([optimum])
pdict2['gvalue'] = fwdpy11.Additive(2., gss)
pdict2['simlen'] = 100
params2 = fwdpy11.ModelParams(**pdict2)
```

We can pass both objects to `tskit` via a {class}`dict`:

```{code-cell} python
ts = pop.dump_tables_to_tskit(model_params={'burnin': params, 'adaptation': params2})
type(ts.metadata["model_params"])
```

```{code-cell} python
assert fwdpy11.ModelParams(**eval(ts.metadata["model_params"]["burnin"])) == params
assert fwdpy11.ModelParams(**eval(ts.metadata["model_params"]["adaptation"])) == params2
```

### Including a `demes` graph

You may include a {class}`demes.Graph` in the top-level metadata.
If you also include a {class}`fwdpy11.ModelParams` (see above), including the `demes` graph gives redundant information.
However, including the graph may be useful if downstream analysis will involve other tools compatible with the `deme` specification.
However, including the graph may be useful if downstream analysis will involve other tools compatible with the `demes` specification.
With the graph as metadata, you can extract it and reconstruct the original `YAML` file, or send it to another Python package that understands it.

The following hidden code block defines a function to return a {class}`demes.Graph` from `YAML` input stored in a string literal.
Expand Down
17 changes: 17 additions & 0 deletions _sources/misc/changelog.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,23 @@
Major changes are listed below. Each release likely contains fiddling with back-end code,
updates to latest `fwdpp` version, etc.

## 0.22.0

Bug fixes:

* Fixed incorrect handling of models with end times > 0 in a demes Graph.
When exporting to tskit, alive node times were treated as 0.0, causing
potential problems for some use cases.
Issue {issue}`1253`
PR {issue}`1255`

Deprecations

* Deprecate the `demes_graph` argument to {func}`fwdpy11.DiploidPopulation.dump_tables_to_tskit`
PR {issue}`1265`
* Deprecate using a dict for the the `model_params` argument to {func}`fwdpy11.DiploidPopulation.dump_tables_to_tskit`
PR {issue}`1265`

## 0.21.6

Behavior changes:
Expand Down
9 changes: 9 additions & 0 deletions _sources/pages/citation.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
# Citation

If you use this software for research, please cite the following publications:

* {cite}`Thornton2019-nu`
* {cite}`Thornton2014-hx`

This software was developed for the first paper.
The second paper describes a key part of this software's back end.
27 changes: 21 additions & 6 deletions _sources/pages/tablefs.md
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ We have to define some conventions:
matrix.

For our examples, we will initialize a {class}`fwdpy11.DiploidPopulation` from
a {class}`tskit.TreeSequence` generated by {func}`msprime.simulate`.
a {class}`tskit.TreeSequence` generated by {func}`msprime.sim_ancestry`.

:::{note}

Expand All @@ -55,17 +55,32 @@ for neutral and non-neutral mutations. See

```{code-cell} python
import fwdpy11
import demes
import msprime
import numpy as np
yaml = """
time_units: generations
demes:
- name: deme0
epochs:
- start_size: 250
- name: deme1
epochs:
- start_size: 250
migrations:
- {demes: [deme0, deme1], rate: 0.1}
"""
graph = demes.loads(yaml)
demography = msprime.Demography.from_demes(graph)
rng = fwdpy11.GSLrng(4321678)
config = [msprime.PopulationConfiguration(500), msprime.PopulationConfiguration(500)]
ts = msprime.simulate(
population_configurations=config,
Ne=500.0,
ts = msprime.sim_ancestry(
samples = {0: 250, 1: 250},
random_seed=777,
migration_matrix=np.array([0, 0.1, 0.1, 0]).reshape(2, 2),
recombination_rate=0.25,
demography = demography,
)
pop = fwdpy11.DiploidPopulation.create_from_tskit(ts)
Expand Down
2 changes: 1 addition & 1 deletion _sources/pages/tstypes.md
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,7 @@ any further.

:::{note}

The trees generated by `msprime.simulate` are bifurcating, which is a consequence of simulating from the Kingman
The trees generated by `msprime.sim_ancestry` are bifurcating, which is a consequence of simulating from the Kingman
coalescent. In forward-time simulations, it is not uncommon to have more than two descendants of a node. When that
happens, left_child and right_child refer to the left-most and right-most children, respectively. Thus, to "walk"
along the descendants of a node, you proceed to left_child, and then march along right_sib until a value of -1 is
Expand Down
5 changes: 5 additions & 0 deletions _sources/pages/types.md
Original file line number Diff line number Diff line change
Expand Up @@ -262,6 +262,11 @@ that is handled on the C++ side.

# Types related to discrete demographic events

```{eval-rst}
.. autoclass:: fwdpy11.ForwardTimeInterval
:members:
```

```{eval-rst}
.. autoclass:: fwdpy11.ForwardDemesGraph
:members:
Expand Down
5 changes: 5 additions & 0 deletions _sources/short_vignettes/demes_vignette.md
Original file line number Diff line number Diff line change
Expand Up @@ -147,3 +147,8 @@ for deme in demography.demes_at_final_generation:

It is probably useful to read more about how ``fwdpy11`` handles ``demes``
models {ref}`here <demes_details>`.

# More information

See [here](demes_event_timings) for concrete examples about reasoning through
when demographic events occur.
6 changes: 6 additions & 0 deletions examples/gss_divergent_optima.html
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,7 @@
<li class="toctree-l1"><a class="reference internal" href="../pages/installation.html">Installation</a></li>
<li class="toctree-l1"><a class="reference internal" href="../pages/deploymenttools.html">Deployment tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../pages/invokingpython.html">Invoking Python</a></li>
<li class="toctree-l1"><a class="reference internal" href="../pages/citation.html">Citation</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Short vignettes - operations on populations</span></p>
<ul class="nav bd-sidenav">
Expand Down Expand Up @@ -186,6 +187,7 @@
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="../short_vignettes/demography_vignettes_intro.html">Introduction</a></li>
<li class="toctree-l1"><a class="reference internal" href="../short_vignettes/demes_vignette.html">Using <code class="docutils literal notranslate"><span class="pre">demes</span></code> to specify demographic models.</a></li>

</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Short vignettes - varying fitness conditions across demes</span></p>
<ul class="nav bd-sidenav">
Expand All @@ -202,6 +204,10 @@
<li class="toctree-l1"><a class="reference internal" href="../short_vignettes/selective_sweep.html">Selective sweeps</a></li>
<li class="toctree-l1"><a class="reference internal" href="../short_vignettes/incomplete_sweep.html">Incomplete selective sweeps</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Longer vignettes - demographic models</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="../long_vignettes/demes_event_timings.html">Timing of events in demographic models</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Longer vignettes - exporting data to tskit</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="../long_vignettes/tskitconvert_vignette.html">Converting data to tskit format</a></li>
Expand Down
Loading

0 comments on commit 40d98fd

Please sign in to comment.