README.md 5.46 KB
Newer Older
bvictor's avatar
bvictor committed
1
# ABMEv.jl Documentation
bvictor's avatar
bvictor committed
2
3
This is a suite for simulating an Agent Based Model that captures the evolutionary dynamics of a population in a multidimensional space.

bvictor's avatar
bvictor committed
4
5
6
7
8
9
## Installation
```julia
using Pkg;
Pkg.add("https://gitlab.ethz.ch/bvictor/abmev.git")
```
This will download latest version from git repo and download all dependencies.
bvictor's avatar
bvictor committed
10
11
12
13
14
To check out from an other branch than master, one has to do the trick
```julia
using Pkg;
Pkg.add("ABMEv#no_C_matrix")
```
bvictor's avatar
bvictor committed
15
16
17
18
## Getting started
```julia
using ABMEv
```
bvictor's avatar
bvictor committed
19

bvictor's avatar
bvictor committed
20
### Gillepsie algorithm
bvictor's avatar
bvictor committed
21

bvictor's avatar
bvictor committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
```julia
p_default = Dict("K0" => 1000.,
        "D" => [1e-2 - 1e-3],
        "mu" => [1. .1],
        "sigma_a" => [5e-1; 7e-1],
        "sigma_K" => [9e-1; 9e-1],
        "n_alpha" => 2.,
        "n_K" => 2.,
        "tend" => 150.,
        "NMax" => Int(2e3))
na_init = 500
world0 = new_world_G(na_init,p_default,spread = [.01 .01], offset = [-.5 -.5])
worldall,p_default["tspan"] = runWorld_store_G(p_default,world0)
```
As of now, no mode is implemented.

### Wright Fisher algorithm
bvictor's avatar
bvictor committed
39
40
#### Specific parameters
- ```NMax``` Maximum number of individuals that can be attained. If attained, then the programm stops.
bvictor's avatar
bvictor committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
```julia
p_default = Dict("K0" => 2000.,
        "D" => [2e-2], # this is the dispersion coefficient, that should be taken as constant
        "mu" => [1.],
        "sigma_a" => [2e0],
        "sigma_K" => [1e0],
        "a" => .95,
        "n_alpha" => 2.,
        "n_K" => 2.,
        "tend" => 10.,
        "NMax" => Int(1e4))
na_init = p_default["K0"]
agent0 = [Agent([-.5 + .01  * randn()]) for i in 1:na_init]
world0 = vcat(agent0[:],repeat([missing],Int(p_default["NMax"] - na_init)))
(worldall , p_default["tspan"]) = runWorld_store_WF(p_default,agent0,mode="std")
```
You have several options available concerning the resource implemented and competition:
- ``` mode="std"``` is the standard mode
- ``` mode="grad2D"``` corresponds to a an ecological gradient
bvictor's avatar
bvictor committed
60
>We are not sure whether this corresponds to the following two images
bvictor's avatar
bvictor committed
61
62
63
64
- ``` mode="mountain"``` corresponds to a scenario where a mountain arises (with an ecological gradient)
- ``` mode="split"``` corresponds to a scenario where the resource is splitted in two
- ``` mode="graph"``` this guy is probably not working

bvictor's avatar
bvictor committed
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
### Mutation
If anisotropy in mutation, the following parameters should be declared as arrays where each entry corresponds to a dimension.
- ```mu``` The probability of mutation.
- ```D``` If mutation happens on the agent, the increment follows $\mathcal{N}_{ 0, D}$
### Birth
#### Growth
- Resource kernel for agent with trait $x$ is defined as 
```math
K_{\mu,\sigma}(x) = K_0 \mathcal{N}_{\mu,\sigma}(x)
```
with $\mu$ and $\sigma$ potentially vectors.
> We just modified this in ABMEv_Agent.jl so you should check if it works.
- Dirth coefficient is defined as $`b(x) = K(x)`$
### Death
#### Competition
- Competition between agent with trait ```x``` and ```y``` is defined as
```math 
\alpha(x,y) = \exp(-\sum_i^{N(t)} \frac{1}{\sigma_{\alpha_i}^{n_\alpha}}\sum_j^T (x_{i,j} - y_{i,j})^{n_\alpha})
```
- Death coefficient is defined as $`d(x^{(i)}) = \sum_j^{N(t)} \alpha(x^{(i)},x^{(j)})`$
> We are not sure if the sum includes $`x^{(i)}`$ or not.
### Fitness
Fitness is defined as ```b - d```.

## Parameter description
- ```K0``` Carrying capacity
- ```a``` only used for mode ```grad2D``` where growth rate is set such as $`\mu = a x_1`$
> We are not sure if this is OK or not? Check it 
[Grad2D kernel explained](https://gitlab.ethz.ch/bvictor/abmev/-/wikis/Grad2D)
bvictor's avatar
bvictor committed
94
95
96
97
98
99
100
101
102

### Parallelism
You can run your script in parallel, which makes sense for large populations. To do so:
```julia
using Distributed;addprocs()
@everywhere push!(LOAD_PATH,homedir()*"path_to_ABMEv")
@everywhere using ABMEv
```

bvictor's avatar
bvictor committed
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
## Properties of agents
You can access properties of the agent using the following functions
- `get_xarray(world::Array{Agent{T}},trait::Int) where T`

Returns trait of every agents of world in the form of an array

> TODO: describe the following accessors
```julia
get_x(a::Agent,i::Number) = a.x_history[Int(i):Int(i),end]
get_x(a::Agent) = a.x_history[:,end]
get_xhist(a::Agent,i::Number) = a.x_history[Int(i):Int(i),:]
get_xhist(a::Agent) = a.x_history
get_geo(a::Agent) = sum(get_xhist(a,1))
get_d(a::Agent) = a.d
get_b(a::Agent) = a.b
get_fitness(a::Agent) = a.b - a.d
```

## Plotting
bvictor's avatar
bvictor committed
122
123
124
125
ABMEv comes with Plot recipes:
`function plot(world::Array{U},p;what=["x","H"],trait = 1) where U <: Union{Missing,Agent{T}} where T`.

An example to use it: 
bvictor's avatar
bvictor committed
126
127
128
129
130
131
132
133
134
135
136
```julia
using Plots;pyplot()
Plots.plot(worldall,p_default,what=["x"],trait=2)
```
You can specify what you want to plot in the array ```what```:
- ```x``` plots the component specified by ```trait=2```
- ```geo``` plots geotrait, computed from first component
- ```3dgeo``` plots a 3d diagram with x axis as geotrait and y axis as the second component
- ```3d``` plots a 3d diagram with first and second component as x and y axis
- ```var``` plots the variance of the  component specified by ```trait=2```
- ```vargeo``` plots the variance of the geotrait
bvictor's avatar
bvictor committed
137
138
139
140
141
142
143
144
145
146

## Developping the code
I recommend to first clone your branch in the directory you like best, and then to 
To develop, you ca
```julia
using Pkg
Pkg.dev("path_to_ABMEv_dir")
```
You can also do the same trick with directly the gitlab address, cf [https://docs.julialang.org/en/v1/stdlib/Pkg/index.html](Pkg.jl)

bvictor's avatar
bvictor committed
147
148
149
## 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)