runworld.jl 4.27 KB
Newer Older
1
"""
Victor's avatar
Victor committed
2
    run!(w, alg, tend, b, d; dt_saving=nothing, cb=(names = String[],agg =nothing))
Victor's avatar
Victor committed
3

Victor's avatar
Victor committed
4
5
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`.
Victor's avatar
Victor committed
6
7
8
9
10
11
12
The run stops if the number of agents reaches`p["NMax"]`, where `p` is the parameter dictionary in the world `w`.


Returns a `Simulation` object. 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.

# Keyword arguments
13
- if `dt_saving` specified, world is saved every time steps.
14
If `dt_saving` not specified, `sim` contains an array of two elements,
Victor's avatar
Victor committed
15
first corresponding to initial conditions and last corresponding to world in the last time step.
Victor's avatar
Victor committed
16
- if `t_saving_cb` specified, callbacks are computed at each steps time specified in the array.
17
This functionality is as of now only compatible with `dt_saving` not specified.
Victor's avatar
Victor committed
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
- `cb` correspond to callbacks function. 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)))])
```

# Example

```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)
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"])
```
68
"""
Victor's avatar
Victor committed
69
function run!(w::World{A,S,T},alg::L,tend::Number,b,d;
70
                dt_saving=nothing,
71
                t_saving_cb=nothing,
Victor's avatar
Victor committed
72
                cb=nothing) where {A,S,T,L<:AbstractAlg}
73
    # argument check
Victor's avatar
Victor committed
74
    _check_timedep(b,d)
75
76
77
78
79
    (!isnothing(dt_saving) && !isnothing(dt_saving)) ? ArgumentError("For now, can not specify both `dt_saving` and `t_saving_cb`") : nothing
    isnothing(dt_saving) ? dt_saving =  tend + 1. : nothing
    isnothing(t_saving_cb) ? t_saving_cb =  [tend + 1.] : nothing

    # var init
Victor's avatar
Victor committed
80
    n=size(w);
Victor's avatar
Victor committed
81
    NMax = maxsize(w)
82
    t = .0
Victor's avatar
Victor committed
83
84
85
    i = 1;j=1;dt = 0.
    sim = Simulation(w,cb=cb)
    if A <: AbstractAgent{AA,Rates{true}} where {AA}
Victor's avatar
Victor committed
86
        update_rates!(w,alg,b,d)
Victor's avatar
Victor committed
87
    end
88
89

    # start
Victor's avatar
Victor committed
90
    while t<tend
91
92
        if dt < 0
            throw("We obtained negative time step dt = $dt at event $i")
Victor's avatar
Victor committed
93
        elseif size(w) == NMax
Victor's avatar
Victor committed
94
            @info "We have reached the maximum number of individuals allowed"
95
            break
Victor's avatar
Victor committed
96
        elseif size(w) == 0
Victor's avatar
Victor committed
97
            @info "All individuals have died :("
98
99
            break
        end
Victor's avatar
Victor committed
100
        if  t - get_tend(sim) >= dt_saving
Victor's avatar
Victor committed
101
            # @info "saving world @ t = $(t)/ $(tend)"
Victor's avatar
Victor committed
102
            add_entry!(sim,w,cb)
103
        end
104
        if  t >= first(t_saving_cb)
Victor's avatar
Victor committed
105
            # @info "saving callback only @ t = $(t)/ $(tend)"
Victor's avatar
Victor committed
106
            add_entry_cb_only!(sim,w,cb)
107
108
            popfirst!(t_saving_cb)
        end
Victor's avatar
Victor committed
109
        dt = updateWorld!(w,alg,b,d)
Victor's avatar
Victor committed
110
        t +=  dt
111
112
        i += 1
    end
Victor's avatar
Victor committed
113
    # Saving last time step
Victor's avatar
Victor committed
114
    add_entry!(sim,w,cb)
vboussange's avatar
vboussange committed
115
    @info "simulation stopped at t=$(t), after $(i) steps"
Victor's avatar
Victor committed
116
    return sim
117
end
118
119

"""
Victor's avatar
Victor committed
120
    _check_timedep(b,d)
121

Victor's avatar
Victor committed
122
Checks number of arguments of functions,
Victor's avatar
Victor committed
123
and throws error if problem
124
"""
Victor's avatar
Victor committed
125
126
function _check_timedep(b,d)
    if numargs(b) < 2
127
128
        throw(ArgumentError("Birth function needs `X` and `t` arguments"))
    end
Victor's avatar
Victor committed
129
    if numargs(d) < 3
130
131
132
        throw(ArgumentError("Death function needs `X`, `Y` and `t` arguments"))
    end
end