Commit 79b0ee60 authored by Victor's avatar Victor
Browse files

updated doc

parent 356fd36a
Pipeline #101415 canceled with stage
......@@ -9,7 +9,7 @@ alt="" width="400"></img> </div>
EvoId.jl (for **Evo**lutionary **I**n**d**ividual-based models) is a package aimed at simulating the eco-evolutionary dynamics of a population in a multidimensional space, at the individual level. The dynamics is specified under the framework of [stochastic models for structured populations](https://arxiv.org/abs/1506.04165).
Individuals are characterised by **a set of traits** in some **combination of evolutionary spaces**. An evolutionary space can represent for example a geographical landscape, a trait space, or genetic structure. Individuals give birth at a rate given by the birth function `b`, and die at a rate given by the death function `d`. When an individual give birth, its offspring can move on the underlying evolutionary spaces. The movement can capture whether migration or mutation processes, and is characterised by a probability `m` and movement range `D`.
Individuals are characterised by **a set of traits** in some **combination of evolutionary spaces**. An evolutionary space can represent for example a geographical landscape, a trait space, or genetic structure. Spaces can be of any dimensions, discrete or continuous, bounded or unbounded. They can equally consist of graphs. Individuals give birth at a rate given by the birth function `b`, and die at a rate given by the death function `d`. When an individual give birth, its offspring can move on the underlying evolutionary spaces. The movement can capture whether migration or mutation processes, and is characterised by a probability `m` and movement range `D`.
The user can provide **any birth and death functions**, which should depend on the system state and the individuals' trait. Together with the **movement rate and movement range**, this defines the dynamics of the system.
......@@ -17,10 +17,10 @@ EvoId.jl provides a **numerical laboratory** for eco-evolutionary dynamics, supp
- flexible types for **individuals**, which can
- evolve over any combination of space,
- store ancestors trait,
- [store ancestors trait](https://vboussange.github.io/EvoId.jl/dev/examples/gradient.html#lineages),
- flexible types for **evolutionary spaces**, that can consist of multidimensional **discrete or continuous domains**, as well as **graphs**,
- the possibility to use **callback functions** to save the state of the system at any time step
- several **algorithms** for the simulations (Gillespie, Wright Fisher, etc...),
- several **algorithms** for the simulations ([Gillespie](https://en.wikipedia.org/wiki/Gillespie_algorithm),[Wright-Fisher](https://en.wikipedia.org/wiki/Moran_process), etc...),
- **utility functions** to analyse simulation results.
## Installation
......@@ -34,7 +34,7 @@ Pkg.add("https://github.com/vboussange/EvoId.jl")
This will download latest version from git repo and download all dependencies.
## Getting started
Check out the tutorial prodived below. You can also dive into the [documentation](https://vboussange.github.io/EvoId.jl/dev) if you want to use the advanced features of EvoId.jl.
Check out the tutorial prodived below. You can also look at the `example` folder, or dive into the [documentation](https://vboussange.github.io/EvoId.jl/dev) if you want to use the advanced features of EvoId.jl.
## Related papers
- [Topology and habitat assortativity drive neutral and adaptive diversification in spatial graphs](https://www.biorxiv.org/content/10.1101/2021.07.06.451404v2), Boussange et al. 2021.
......@@ -54,7 +54,6 @@ We provide here a tutorial that sums up the 5 steps necessary to launch a simula
Let's import EvoId.jl, and LightGraphs.jl
```julia
using EvoId
using LightGraphs
```
### 1. Define the evolutionary spaces
......@@ -69,7 +68,7 @@ evolspace = (landscape,traitspace)
```
### 2. Define birth and death function
Birth and death functions depend on agents position in the combination of spaces defined above, i.e. position on the graph and the adaptive trait.
Birth and death functions depend on individuals position in the combination of spaces defined above, i.e. position on the graph and the adaptive trait.
We decide that each vertex selects for an optimal trait value $`\theta_i \in \{-1,1\}`$.
```julia
......
......@@ -9,7 +9,11 @@ makedocs(sitename="EvoId.jl",
"Manual" => ["manual/agent.md",
"manual/space.md",
"manual/world.md",
"manual/run_world.md"],
"manual/run_world.md",
"manual/plots.md",
"manual/diversity.md",
"manual/utils.md",
],
"Examples" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"examples"),join=true))],
"Developping" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"dev"),join=true))],
# "contributing.md",
......
docs/src/assets/logo.png

85.9 KB | W: | H:

docs/src/assets/logo.png

419 KB | W: | H:

docs/src/assets/logo.png
docs/src/assets/logo.png
docs/src/assets/logo.png
docs/src/assets/logo.png
  • 2-up
  • Swipe
  • Onion skin
# Developping the code
I recommend to first clone your branch in the directory you like best, and then to
To develop, you can
I recommend to first clone your branch in the directory you like best, and then to develop, you can
```julia
using Pkg
Pkg.dev("path_to_EvoId_dir")
......@@ -8,13 +7,11 @@ Pkg.dev("path_to_EvoId_dir")
You can also do the same trick with directly the gitlab address, cf [Pkg.jl](https://docs.julialang.org/en/v1/stdlib/Pkg/index.html)
## Future directions
- Try to improve parallelism with the help of Threads.@Threads and @inbounds (cf [tutorial](https://juliagpu.gitlab.io/CUDA.jl/tutorials/introduction/#Introduction-1) )
- Make use of CUDA.jl to accelerate the simulations wih the use of GPU.
## Todo
- We do not need to have the `Rates{}` parameter for `Agent` type.
- Fix `WF` algorithm
- Implement Moran process
- extend the birth function to the form `b(X,Y,t)` for coherence with death function
- make mutation and disperal range features of agents, so that they can also evolve.
- Simplify composite type `Rates{}` parameter for `Agents`.
```julia
abstract type AbstractAgent{A<:Ancestors,R<:Rates} end
```
- extend the birth function to the form `b(X,Y,t)` for coherence with death function
- make mutation and disperal range features of agents, so that they can also evolve.
# EvoId.jl: Evolutionary Individual Based Model
# EvoId.jl
EvoId.jl (for **Evo**lutionary **I**n**d**ividual-based model) is a package aimed at simulating the evolutionary dynamics of a population in a multidimensional space. The population is modelled at the individual level. Individuals experience four elementary events : birth, death, mutation and migration.
- EvoId.jl hence falls in the realm of *Agent Based Model*.
EvoId.jl provides a numerical laboratory for evolutionary dynamics, supplying
- flexible types for individuals, which can
- evolve over any combination of space
- store ancestors trait,
- flexible types for evolutionary spaces, comprising multidimensional discrete and continuous sets, as well as graphs,
- the possibility for the user to provide any birth and death functions,
- several algorithms for the simulations,
- utility functions to analyse simulation results.
## Features
Agents consist of sets of traits in some combination of evolutionary spaces. An evolutionary space can represent for example a geographical landscape, a trait space, or genetic structure. Spaces can be of any dimensions, discrete or continuous, bounded or unbounded. They can equally consist of graphs.
Vector spaces are used to define birth and death processes, as well as mutation processes.
### Specificities
- [EvoId.jl allows to keep track of agents' trait lineages](@ref lineages)
- [EvoId.jl enables to run evolutionary dynamics on graphs!](@ref genetic_structure)
## Getting started
```@repl
using EvoId
```
## Tutorial
We strongly advise to have a look at the tutorial section. All the scripts of the examples can be found [here](https://gitlab.ethz.ch/bvictor/EvoId/-/tree/master/examples).
```@contents
Pages = [
"examples/delta_competition_example.md",
"examples/changing_environment.md",
"examples/sympatric_speciation.md",
"examples/gradient.md",
"examples/genetic_structure.md",
""
]
Depth = 2
```
## How it works
General workflow to launch any simulation is the following
- [Define the combination of vector spaces you are interested in.](manual/space.md)
- Define birth and death function, that depend on agents position in the space
- Define mutation function
- [Define initial population state and time](manual/world)
- [Run the simulation according to some updating algorithm](manual/run_world.md)
- [Obtain a summary of the population state](manual/callbacks.md)
### Available algorithms
As of now, three types of simulation algorithm can be used:
- [Gillepsie](manual/gillepsie.md)
- [Wright-Fisher](manual/wright_fisher.md)
- [CFM](CFM.md)
## References
- [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632)
- [Champagnat and Ferriere second article - 2008](https://www.tandfonline.com/doi/full/10.1080/15326340802437710)
## Similar packages
[Agents.jl](https://juliadynamics.github.io/Agents.jl/) This package is oriented towards general ABM modelling, and thus is not as efficient and easy to deploy as EvoId.jl for simulating stochastic models of structured populations.
EvoId.jl (for **Evo**lutionary **I**n**d**ividual-based models) is a package aimed at simulating the eco-evolutionary dynamics of a population in a multidimensional space, at the individual level. The dynamics is specified under the framework of [stochastic models for structured populations](https://arxiv.org/abs/1506.04165).
\ No newline at end of file
# Agent properties
!!! note "Specificities"
The type `Agent` has two important composite types
- `Ancestors{bool}` : when `bool = true`, the ancestors traits are stored,
- `Rates{bool}` : when `bool = true`, the rates `d` and `b` of agents are updated at each time step. This is needed in e.g. Gillepsie Algorithm
# Agent
```@autodocs
Modules = [EvoId]
......
# Callbacks
Callbacks can be used to extract properties of the world at each `dt_saving` time steps of your simulation.
## Constructing the callbacks
A callback has to be of the form
```julia
cb = (names = String[], agg = Function[])
```
It is a tuple, with first value corresponding to the names of the aggregate properties of the world. The second correspond to the aggregation functions.
We provide here an example on how to extract the ``\gamma`` diversity of a simulation biological population. ``\gamma`` diversity can be calculated as the variance of the trait distribution of the population.
Here is how we write the function
```julia
cb = (names = ["gamma_div"], agg = Function[w -> var((get_x(w,1)))])
```
Here is how we use it
```julia
myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
b(X) = gaussian(X[1],0.,sigma_K)
d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
D = (1e-2,)
mu = [.1]
NMax = 10000
tend = 1.5
p = Dict{String,Any}();@pack! p = D,mu,NMax
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
w1 = copy(w0)
@time sim = run!(w1,Gillepsie(),tend,b,d,cb=cb,dt_saving = .1)
```
## Accessing the callbacks
You can easily access the properties, using `sim` as you would for a usual `Dictionary`.
```julia
using Plots
plot(get_tspan(sim),sim["gamma_div"])
```
# Diversity measures
# Metrics
```@autodocs
Modules = [EvoId]
......
......@@ -5,6 +5,7 @@ run!
```
For now three algorithms
```@docs
Gillepsie
CFM
......
# Simulation
A `Simulation` object is returned by the function `run!`. It is a container for snapshots of the world at every `dt_saving` time steps. It renders post processing easier, through dedicated methods to obtain time series of quantities.
!!! note "Default behaviour"
If `df_saving` is not provided, initial and last time steps will be saved.
```@autodocs
Modules = [EvoId]
Pages = ["Sim.jl"]
```
# World
```@autodocs
Modules = [EvoId]
Pages = ["world.jl"]
Pages = ["World.jl"]
```
......@@ -7,6 +7,12 @@ export AbstractAgentM
"""
$(TYPEDEF)
!!! note "Specificities"
The type `Agent` has two important composite types
- `Ancestors{bool}` : when `bool = true`, the ancestors traits are stored,
- `Rates{bool}` : when `bool = true`, the rates `d` and `b` of agents are updated at each time step. This is needed in e.g. `Gillepsie` algorithm.
"""
mutable struct Agent{A<:Ancestors,R<:Rates,T<:Tuple,U,V} <: AbstractAgent{A,R}
# history of traits for geotraits
......@@ -49,10 +55,9 @@ with initial position `pos`. If `pos` not provided,
initialises agent with 0 values everywhere
# Keyword arguments
* `rates`. Set `rates=true` when agents fitness
needs to be updated at each time step. This
* `rates`. Set `rates=true` when agents fitness needs to be updated at each time step. This
is required for the Gillepsie algorithm, but not for CFM algorithm
- `ancestors`. Set `ancestors=true` when you want to store ancestors traits.
* `ancestors`. Set `ancestors=true` when you want to store ancestors traits.
"""
function Agent(s::S; ancestors=false,rates=true) where {S <: AbstractSpacesTuple}
T,pos = _initpos(s)
......
module EvoId
using Distributions,LinearAlgebra,Reexport,StatsBase
using LightGraphs
using UnPack
using DocStringExtensions
using Arpack
......@@ -18,7 +17,7 @@ module EvoId
include("runworld.jl")
@reexport using Distributions, DataFrames
@reexport using Distributions, DataFrames, LightGraphs
export
......
......@@ -13,7 +13,9 @@ end
"""
$(SIGNATURES)
`Simulation` object, output by the method `run!`
`Simulation` object, output by the method `run!`.
See `run!` documentation for more insights.
"""
function Simulation(w0::World{A,S,T};cb = nothing) where {A,S,T}
tspan = zeros(1)
......
......@@ -22,13 +22,15 @@ abstract type AbstractStatSpace{Dim,T,I} <: AbstractSpace{Dim,T,I} end
"""
$(TYPEDEF)
Creates a Graph Space.
# Example
```julia
using LightGraphs
g = star_graph(7)
GraphSpace(g)
```
Creates a Graph Space.
# Example
```julia
using LightGraphs
g = star_graph(7)
GraphSpace(g)
```
"""
struct GraphSpace{T} <: AbstractStatSpace{1,T,IsFinite{true}}
......@@ -39,14 +41,16 @@ abstract type AbstractSegment{T<:Number} <: AbstractStatSpace{1,T,IsFinite{true
"""
$(TYPEDEF)
Creates a segment space, where individuals are reflected at both ends.
# Arguments
* `s` start of the segment
* `e` end of the segment
# Example
```julia
ContinuousSegment(1., 2.)
```
Creates a segment space, where individuals are reflected at both ends.
# Arguments
* `s` start of the segment
* `e` end of the segment
# Example
```julia
ContinuousSegment(1., 2.)
```
"""
struct ContinuousSegment{T<:AbstractFloat} <: AbstractSegment{T}
s::T
......@@ -55,14 +59,18 @@ end
"""
$(TYPEDEF)
Creates a discrete segement space, where individuals are reflected at both ends.
# Arguments
* `s` start of the segment
* `e` end of the segment
# Example
```julia
DiscreteSegment(1, 2)
```
Creates a discrete segement space, where individuals are reflected at both ends.
# Arguments
* `s` start of the segment
* `e` end of the segment
# Example
```julia
DiscreteSegment(1, 2)
```
"""
struct DiscreteSegment{T<:Integer} <: AbstractSegment{T}
s::T
......@@ -71,10 +79,12 @@ end
"""
$(TYPEDEF)
Creates a real space.
# Arguments
* `N` dimension of the space
* `T` type of the underlying traits.
Creates a real space.
# Arguments
* `N` dimension of the space
* `T` type of the underlying traits.
"""
struct RealSpace{N,T} <: AbstractStatSpace{N,T,IsFinite{false}} end
RealSpace(N) = RealSpace{N,Float64}()
......@@ -89,6 +99,7 @@ struct NaturalSpace{N,T} <: AbstractStatSpace{N,T,IsFinite{false}} end
"""
$(SIGNATURES)
Returns increment corresponding to space `s`
"""
get_inc(x,D,s::AbstractStatSpace,t) = get_inc(x,D,s) # this is defined to skip representation of t for following specialised methods
......@@ -145,14 +156,16 @@ end
abstract type AbstractDynSpace{Dim,T<:Number} <: AbstractSpace{Dim,T,IsFinite{true}} end
"""
$(TYPEDEF)
A dynamic graph space.
# Arguments
* `g` the underlying graph
* `f` a function that takes as argument time,
and returns the index of the graph to pick at time `t` from array `g`
# Example
`DynGraphSpace(g,f)`
A dynamic graph space.
# Arguments
* `g` the underlying graph
* `f` a function that takes as argument time,
and returns the index of the graph to pick at time `t` from array `g`
# Example
`DynGraphSpace(g,f)`
"""
struct DynGraphSpace{T<:Number} <: AbstractDynSpace{1,T}
g::Vector{AbstractGraph{T}}
......@@ -165,6 +178,7 @@ function DynGraphSpace(g::Array{A},f) where A <: AbstractGraph
"""
$SIGNATURES
Returns the graph correseponding to `d::DynGraphSpace` at time `t`
"""
get_graph(d::DynGraphSpace,t) = d.g[d.f(t)]
......
......@@ -9,41 +9,44 @@ end
"""
$(SIGNATURES)
Constructs a world.
# Arguments
* `w` a vector of agents,
* `s` a tuple of evolutionary spaces,
* `p` a dictionary that contains
- `"mu"`, a vector of same length as `s`,
that contains the mutation rate. If `mu[i]`
has dimension greater than 1,
then mutations happen independently at each dimension
of `s[i]`.
- `"sigma"`, a vector of same length as `s`,
that contains the dispersal ranges. Only `nothing` is
supported for `GraphSpace`, equivalent to a random walk of length 1.
- `"NMax"` the maximum number of individuals allowed during the simulation
# Examples
```julia
nodes = 7
g = star_graph(nodes)
landscape = GraphSpace(g)
θ = [rand([-1,1]) for i in 1:nodes]
traitspace = RealSpace(1)
evolspace = (landscape,traitspace)
D = [nothing,5e-2]
mu = [1f-1,1f-1]
p = Dict("NMax" => 2000,
"D" => D,
"mu" => mu)
myagents = [Agent(evolspace,[rand(1:nodes),randn() * D[2]]) for i in 1:K]
Constructs a world.
# Arguments
* `w` a vector of agents,
* `s` a tuple of evolutionary spaces,
* `p` a dictionary that contains
- `"mu"`, a vector of same length as `s`,
that contains the mutation rate. If `mu[i]`
has dimension greater than 1,
then mutations happen independently at each dimension
of `s[i]`.
- `"sigma"`, a vector of same length as `s`,
that contains the dispersal ranges. Only `nothing` is
supported for `GraphSpace`, equivalent to a random walk of length 1.
- `"NMax"` the maximum number of individuals allowed during the simulation
w0 = World(myagents,evolspace,p)
```
# Examples
```julia
nodes = 7
g = star_graph(nodes)
landscape = GraphSpace(g)
θ = [rand([-1,1]) for i in 1:nodes]
traitspace = RealSpace(1)
evolspace = (landscape,traitspace)
D = [nothing,5e-2]
mu = [1f-1,1f-1]
p = Dict("NMax" => 2000,
"D" => D,
"mu" => mu)
myagents = [Agent(evolspace,[rand(1:nodes),randn() * D[2]]) for i in 1:K]
w0 = World(myagents,evolspace,p)
```
"""
function World(w::Vector{A},s::S,p::Dict;t::T=0.) where {A<:AbstractAgent,S<:AbstractSpacesTuple,T}
......@@ -104,6 +107,7 @@ Base.copy(w::W) where {W<:World} = W(copy.(w.agents),w.space,w.p,copy(w.t))
## Accessors
"""
$(SIGNATURES)
Get x of world without geotrait.
"""
Base.getindex(w::World,i) = w.agents[i]
......@@ -132,6 +136,7 @@ get_geo(w::World) = map(a-> get_geo(a,time(w)), agents(w))
"""
$(SIGNATURES)
Returns trait of every agents of world in the form of an array which dimensions corresponds to the input.
If `trait = 0` , we return the geotrait.
!!! warning "Warning"
......@@ -151,6 +156,7 @@ end
"""
$(SIGNATURES)
Returns every traits of every agents of `world` in the form **of a one dimensional array** (in contrast to `get_x`).
If `geotrait=true` the geotrait is also added to the set of trait, in the last column.
If you do not want to specify `t` (only useful for geotrait), it is also possible to use `get_xarray(world::Array{T,1}) where {T <: Agent}`.
......@@ -167,7 +173,7 @@ end
@deprecate get_xarray(world,geotrait=false) get_x(world,Colon())
"""
function give_birth(mum_idx::Int,w::World)
give_birth(mum_idx::Int,w::World)
Copies agent within index `mum_idx`, and increment it by dx.
Return new agent (offspring).
"""
......
......@@ -6,12 +6,14 @@ Gillespie algorithm.
Denote by ``b_i`` and ``d_i`` the birth and death rates of
agent ``i``. The total rate is given by the sum of all individual rates
``R(t) = \\left[ \\sum_i b_i(t) + d_i(t) \\right]``
A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate ``R``.
``R(t) = \\left[ \\sum_i b_i(t) + d_i(t) \\right]``.
A particular event, birth or death,
is chosen at random with a probability equal to the rate of this event divided by the total rate ``R``.
# The original article by Gillsepie:
[**A general method for numerically simulating the stochastic time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub)
[**A general method for numerically simulating the stochastic
time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub)
"""
struct Gillepsie <: AbstractAlg end
......
......@@ -6,7 +6,7 @@ hamming(k::Array, l::Array) = sum(k .!= l)
"""
H_discrete(s)
Interconnectedness measure as in Nordbotten 2018 for discrete setup
Interconnectedness measure as in Nordbotten 2018 for discrete setup.
"""
function H_discrete(s)
N = length(s)
......@@ -22,9 +22,7 @@ end
"""
findclusters(v::Vector,allextrema =true)
Returns a tuple with the cluster mean and its associated weight
# Arguments
Returns a tuple with the cluster mean and its associated weight.
"""
function findclusters(v::Vector,allextrema =true)
# this function fits a histogram to some distribution and extracts its local maxima
......@@ -44,7 +42,7 @@ end
import Statistics:var,mean
# TODO: rename this to gamma diversity
"""
function var(world::World;trait=1)
var(world::World;trait=1)
Return the variance of the `world`'s `trait` distribution.
If trait = 0, returns the variance of the geotrait,
knowing that by default it is associated with position trait 1.
......@@ -65,7 +63,7 @@ function var(world::World;trait=1)
end
"""
function mean(world::World;trait=1)
mean(world::World;trait=1)
Returns the mean of the `world`'s `trait` distribution.
If trait = 0, returns the variance of the geotrait,
"""
......@@ -98,7 +96,7 @@ function covgeo(world::World,trait = 0)
end
"""
function hamming(world::Array{Agent,1})
hamming(world::Array{Agent,1})
Returns a matrix H where H_ij = hamming(a_i,a_j).
The hamming distance is taken through the whole history
and functional space of the agents.
......@@ -164,13 +162,14 @@ end
"""
$(SIGNATURES)
get_xhist_mat(agentarray, trait=1, time = 0) where {A<:AbstractAgent}
returns `xhist,ttot`, where `xhist` is a matrix with dimension `lenght(world)` x `length(thist)+1`,
which consists of geographical history of ancestors at every time step.
If `time` is specified and is higher that the highest time step in world,
then the last column of xhist corresponds to actual position of agents
"""
function get_xhist_mat(agentarray::Vector{A},trait=1,time = 0) where {A<:AbstractAgent}
function get_xhist_mat(agentarray::Vector{A}, trait=1, time = 0) where {A<:AbstractAgent}
thist = vcat(get_thist.(agentarray)...)
ttot = sort!(unique(thist))
xhist = zeros(length(agentarray),length(ttot))
......
"""
function run!(w::World{A,S,T},alg::L,tend::Number,b,d;dt_saving=nothing,cb=(names = String[],agg =nothing))
run!(w, alg, tend, b, d; dt_saving=nothing, cb=(names = String[],agg =nothing))
Run `w` with algorithm `alg`, until `tend` is reached. User needs to provide `b` the birth function,
which takes as arguments `X,t`, and provide `d` the death function, with arguments `X,Y,t`.
Returns a `Simulation` type.
The run stops if the number of agents reaches`p["NMax"]`, where `p` is the parameter dictionary in the world `w`.