/
simulate.jl
138 lines (112 loc) · 4.12 KB
/
simulate.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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
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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
"""
**Main simulation loop**
simulate(parameters, biomass; start::Int64=0, stop::Int64=500, use::Symbol=:stiff)
This function takes two mandatory arguments:
- `p` is a `Dict` as returned by `make_parameters`
- `biomass` is an `Array{Float64, 1}` with the initial biomasses of every species
Internally, the function will check that the length of `biomass` matches
with the size of the network.
In addition, the function takes three optional arguments:
- `start` (defaults to 0), the initial time
- `stop` (defaults to 500), the final time
- `use` (defaults to `:stiff`), a hint to select the solver
The integration method is, by default, `:stiff`, and can be changed to
`:nonstiff`. This is because internally, this function used the
`DifferentialEquations` package to pick the most appropriate algorithm.
The `simulate` function returns a `Dict{Symbol, Any}`, with three
top-level keys:
- `:p`, the parameters that were given as input
- `:t`, the timesteps
- `:B`, an `Array{Float64, 2}` with the biomasses
The array of biomasses has one row for each timestep, and one column for
each species.
"""
function simulate(parameters, biomass; n_concentration::Vector{Float64}=rand(Float64, 2).*10, start::Int64=0, stop::Int64=500, use::Symbol=:nonstiff, cb_interp_points::Int64=100, extinction_threshold::Float64=1e-6)
@assert stop > start
@assert length(biomass) == size(parameters[:A],1)
@assert length(n_concentration) == 2
if parameters[:productivity] == :nutrients
biomass = vcat(biomass, n_concentration)
end
@assert use ∈ vec([:stiff :nonstiff])
alg = use == :stiff ? Rodas4(autodiff=false) : Tsit5()
S = size(parameters[:A], 1)
# Pre-allocate the timeseries matrix
tspan = (float(start), float(stop))
t_keep = collect(start:0.25:stop)
# Perform the actual integration
prob = ODEProblem(dBdt, biomass, tspan, parameters)
ϵ = []
function species_under_extinction_threshold(u, t, integrator)
workingbm = deepcopy(u)
sort!(ϵ)
deleteat!(workingbm, unique(ϵ))
cond = any(x -> x < extinction_threshold, workingbm) ? -0.0 : 1.0
return cond
end
function remove_species!(integrator)
u = integrator.u
idϵ = findall(x -> x < extinction_threshold, u)
for e in idϵ
if !(e ∈ ϵ)
u[e] = 0.0
append!(ϵ,e)
end
end
sort!(ϵ)
nothing
end
function remove_species_and_update!(integrator)
remove_species!(integrator)
if parameters[:productivity] == :nutrients
workingbm = deepcopy(integrator.u[1:end-2])
else
workingbm = deepcopy(integrator.u)
end
parameters = update_rewiring_parameters(parameters, workingbm, integrator.t)
end
function remove_target_and_update!(u, t, integrator)
remove_species!(integrator)
if parameters[:productivity] == :nutrients
workingbm = deepcopy(integrator.u[1:end-2])
else
workingbm = deepcopy(integrator.u)
end
parameters = BioEnergeticFoodWebs.update_rewiring_parameters(parameters, workingbm, integrator.t)
end
cb = species_under_extinction_threshold
affect_function = remove_species_and_update!
if parameters[:rewire_method] == :ADBM
if parameters[:adbm_trigger] == :interval
Δt = parameters[:adbm_interval]
cb1 = PeriodicCallback(affect_function, Δt)
else
cb1 = ContinuousCallback(cb, affect_function, interp_points = cb_interp_points)
end
else
cb1 = ContinuousCallback(cb, affect_function, interp_points = cb_interp_points)
end
is_any_extinct = any(biomass .< extinction_threshold)
if is_any_extinct
cb2 = FunctionCallingCallback(remove_target_and_update!, funcat = [1.0])
extinction_callback = CallbackSet(cb1, cb2)
else
extinction_callback = cb1
end
sol = solve(prob, alg, callback = extinction_callback, saveat=t_keep, dense=false, save_timeseries=false, force_dtmin=false)
B = hcat(sol.u...)'
if parameters[:productivity] == :nutrients
output = Dict{Symbol,Any}(
:p => parameters,
:t => sol.t,
:B => B[:,1:S],
:C => B[:,S+1:end]
)
else
output = Dict{Symbol,Any}(
:p => parameters,
:t => sol.t,
:B => B)
end
return output
end