To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

delta_competition_example.md 4.13 KB
Newer Older
Victor's avatar
Victor committed
1 2 3 4 5
# A first model of meta population

In this script, we model the evolution of a population where agents are simply defined by their position on some landscape. We implement the simplest possible birth and death function.

## The landscape
Victor's avatar
Victor committed
6 7 8
Let's start by a linear landscape. We define a discrete segment of length `9`, with reflecting boundary conditions. In fact, reflecting boundary conditions are implemented for any finite space.

!!! warning "1D Reflections"
Victor's avatar
Victor committed
9
    As of v1, only reflections on one dimensional space are implemented. We have a two dimensional reflection method that will be released in the future.
Victor's avatar
Victor committed
10 11 12 13 14 15 16 17

There are two ways of implementing a linear landscape. The first one uses a `DiscreteSegment` while the second relies on LightGraphs.jl library. Both are *almost* equivalent.

### DiscreteSegment
```julia
using ABMEv
nodes = 10
mysegment = DiscreteSegment(1,nodes)
Victor's avatar
Victor committed
18
wholespace = (mysegment,)
Victor's avatar
Victor committed
19 20 21 22 23 24 25 26
```

### grid
```julia
using ABMEv, LightGraphs
nodes = 10
g = grid([nodes,1])
mysegmentgraph = GraphSpace(g)
Victor's avatar
Victor committed
27
wholespace = (mysegmentgraph,)
Victor's avatar
Victor committed
28
```
Victor's avatar
Victor committed
29 30
!!! warning "Space tuple"
    Notice that the whole space should be a *tuple of spaces*. Even where there is only one sub vector space as here, you need to have brackets and comma around the unit vector space.
Victor's avatar
Victor committed
31 32 33
Here is how you can visualise the landscape.

```julia
Victor's avatar
Victor committed
34 35
using GraphPlot
gplot(g, collect(1:nodes), collect(1:nodes))
Victor's avatar
Victor committed
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
```

## Defining competition processes
We propose that any individual have a constant birth rate, and competes with all the individuals present in the same patch. Let ``i \in \N``,``x_i \in \{1,2,\dots,9\}``.
The competition pressure experience by individual ``i`` is such that

```math
d(x_i) = \sum_j \delta(x_i-x_j)
```
where ``\delta`` is the dirac function.

In this way, we recover a logistic growth function for subpopulation within a patch.

```julia
K0 = 1000 # We will have in total 1000 individuals
b(X,t) = 1 / nodes
d(X,Y,t) = (X[1]  Y[1]) / K0
```
At equilibrium, population size in each deme will reach `K0 / nodes`.

Victor's avatar
Victor committed
56 57 58
!!! warning "Time dependency"
    Even though time is not used, you have to write birth and death functions with time dependency.

Victor's avatar
Victor committed
59 60 61 62 63 64 65 66 67
## Dispersal
We assume that anytime an offspring is born, it is given a chance to move (`\mu = 1`).
```julia
mu = [1.]
D = (1.5,)
```
## Running the world
We initialise the world with initial population of size ``K_0 / 9`` located on patch `5`. We keep track of individuals' ancestors by setting `ancestors=true`. Because we wish to use `Gillepsie` algorithm, we need `rates=true` as agents' internal birth and death rates are updated at every time step.
!!! note "Warning"
Victor's avatar
Victor committed
68
    `rates` treatment is something we might implement in the library internals.
Victor's avatar
Victor committed
69 70

```julia
Victor's avatar
Victor committed
71
using UnPack# useful macro @pack!
Victor's avatar
Victor committed
72
NMax = 2000
Victor's avatar
Victor committed
73
tend = 300.
Victor's avatar
Victor committed
74 75 76 77 78
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
myagents = [Agent(myspace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes]
w0 = World(myagents,myspace,p,0.)
@time sim = run!(w0,Gillepsie(),tend)
```
Victor's avatar
Victor committed
79 80
This is the simplest run you can do. Now time to more interesting things

Victor's avatar
Victor committed
81
## Analysis
Victor's avatar
Victor committed
82 83

### Size of the world
Victor's avatar
Victor committed
84
Let's verify that the population's growth is logistic. We will plot the population size over time.
Victor's avatar
Victor committed
85
To do so, one need to define `dt_saving` < `tend` to save every `dt_saving` time steps of the world.
Victor's avatar
Victor committed
86 87

```julia
Victor's avatar
Victor committed
88 89 90 91
myagents = [Agent(wholespace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes]
w0 = World(myagents,wholespace,p,0.) # we need to reinitialise the world
@time sim = run!(w0,Gillepsie(),tend,dt_saving=2.)
wsize = [length(w) for w in sim[:]]
Victor's avatar
Victor committed
92
using Plots
Victor's avatar
Victor committed
93 94 95 96 97
Plots.plot(get_tspan(sim),wsize,
                label = "",
                ylabel = "Metapopulation size",
                xlabel ="time",
                grid = false)
Victor's avatar
Victor committed
98
```
Victor's avatar
Victor committed
99 100
![](tutorials/delta_comp_wsize.png)

Victor's avatar
Victor committed
101

Victor's avatar
Victor committed
102 103
!!! notes "Callbacks"
    Note that one could also use a callback function to obtain time series of size of the world computed at simulation time. See [Simulation page](../manual/simulation.md).
Victor's avatar
Victor committed
104

Victor's avatar
Victor committed
105 106 107 108 109 110 111 112 113 114 115 116
### Position through time

One might be tempted to plot the agents position for some time steps.

```julia
Plots.plot(sim,
        label = "",
        ylabel = "Geographical position",
        grid = false,
        markersize = 10)
```
![](tutorials/delta_comp_pos.png)