ABMEv_runworld.jl 8.16 KB
Newer Older
1
2
# for now we consider that competition is local within an array

Victor Boussange's avatar
Victor Boussange committed
3
"""
4
    update_rates_std!(world,p::Dict,t::Float64)
Victor Boussange's avatar
Victor Boussange committed
5
6
7
8
This standard updates takes
    - competition kernels of the form α(x,y) and
    - carrying capacity of the form K(x)
"""
9
function  update_rates_std!(world,p::Dict,t::Float64)
Victor Boussange's avatar
Victor Boussange committed
10
    α = p["alpha"];K=p["K"];
11
12
13
    traits = get_x.(world)
    # traits = get_xhist.(world)
    N = length(traits)
14
    D = SharedArray{Float64}(N)
15
16
17
    # Here you should do a shared array to compute in parallel
    @sync @distributed for i in 1:(N-1)
        for j in i+1:N
18
            C = α(traits[i],traits[j])
19
20
            D[i] += C
            D[j] += C
21
22
23
24
        end
    end
    # Here we can do  it in parallel as well
    for (i,a) in enumerate(world)
25
        a.d = D[i]
26
        a.b = K(traits[i])
27
28
29
30
31
32
    end
end

function  update_rates_graph!(world,C,p::Dict,t::Float64)
    for e in edges(p["g"])
        # agents on e
33
        aidx_e = findall(a -> get_x(a,1)==e,world)
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
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
        na = length(aidx_e)
        for i in aidx_e
            world[i].d = na^2
            world[i].b = 1
        end
    end
end

function  update_rates_2D!(world,C,p::Dict,t::Float64)
    # Competition matrix
    traits = get_x.(world)
    # traits = get_xhist.(world)
    N = length(traits)
    # C = SharedArray{Float64}((N,N))
    # Here you should do a shared array to compute in parallel
    @sync @distributed for i in 1:(N-1)
        C[i,i] = 1.
        for j in i+1:N
            # be careful, for xhist what is return is an array hence no need of []
            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
            C[j,i] = C[i,j]
        end
    end
    C[N,N] = 1.
    # Here we can do  it in parallel as well
    for (i,a) in enumerate(world)
        # we only update death rate
        # a.d = sum(C[i,:]) / K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1])
        a.d = sum(C[i,:])
        # /!| not ideal to assign at every time step the birth rate that is constant
        a.b = KK(traits[i][1:1],
                    p["K0"],
                    p["n_K"],
                    p["sigma_K"][1:1],
                    split_merge_move(t),
                    -split_merge_move(t))/K(traits[i][2:2],
                    p["K0"],
                    p["n_K"],
                    p["sigma_K"][2:2])
    end
end

function  update_rates_grad2D!(world,C,p::Dict,t::Float64)
    # Competition matrix
    traits = get_x.(world)
    # traits = get_xhist.(world)
    N = length(traits)
    # C = SharedArray{Float64}((N,N))
    # Here you should do a shared array to compute in parallel
    @sync @distributed for i in 1:(N-1)
        C[i,i] = 1.
        for j in i+1:N
            # be careful, for xhist what is return is an array hence no need of []
            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
            C[j,i] = C[i,j]
        end
    end
    C[N,N] = 1.
    # Here we can do  it in parallel as well
    for (i,a) in enumerate(world)
        # we only update death rate
        a.d = sum(C[i,:])
        a.b = K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1])
    end
end

function  update_rates_mountain!(world,C,p::Dict,t::Float64)
    # Competition matrix
    traits = get_x.(world)
    # traits = get_xhist.(world)
    N = length(traits)
    # C = SharedArray{Float64}((N,N))
    # Here you should do a shared array to compute in parallel
    @sync @distributed for i in 1:(N-1)
        C[i,i] = 1.
        for j in i+1:N
            # be careful, for xhist what is return is an array hence no need of []
            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
            C[j,i] = C[i,j]
        end
    end
    C[N,N] = 1.
    # Here we can do  it in parallel as well
    for (i,a) in enumerate(world)
        # we only update death rate
        a.d = sum(C[i,:])
        a.b = KK(traits[i][1:1],
                p["K0"],
                p["n_K"],
                p["sigma_K"][1:1],
                split_move(t),
                -split_move(t))/K(traits[i][2:2],
                p["K0"],
                p["n_K"],
                p["sigma_K"][2:2],
                μ=(tin(traits[i][1],-1.,0.)*(traits[i][1]+1) - tin(traits[i][1],0.,1.)*(traits[i][1]-1) ) * split_move(t))
    end
end

function  update_rates_std_split!(world,C,p::Dict,t::Float64)
    # Competition matrix
    traits = get_x.(world)
    # traits = get_xhist.(world)
    N = length(traits)
    # C = SharedArray{Float64}((N,N))
    # Here you should do a shared array to compute in parallel
    @sync @distributed for i in 1:(N-1)
        C[i,i] = 1.
        for j in i+1:N
            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
            C[j,i] = C[i,j]
        end
    end
    C[N,N] = 1.
    # Here we can do  it in parallel as well
    for (i,a) in enumerate(world)
        # we only update death rate
        # this could be imptrove since \alpha is symmetric, by using a symmetric matrix
        # a.d = sum(C[i,:]) / KK(traits[i][:,end],p["K0"],p["n_K"],p["sigma_K"],split_merge_move(t),-split_merge_move(t))
        a.d = sum(C[i,:])
        # /!| not ideal to assign at every time step the birth rate that is constant
        a.b =  KK(traits[i][:,end],p["K0"],p["n_K"],p["sigma_K"],split_move(t),-split_move(t))
    end
end
"""
Victor's avatar
Victor committed
159
    function runWorld_store_WF(p,world0;mode="std")
160
161
Wright Fisher process. Returns an array worldall with every agents.
"""
Victor's avatar
Victor committed
162
function runWorld_store_WF(p,world0;mode="std")
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
    worldall = repeat(world0,inner = (1,length(1:Int(p["tend"]))))
    N=length(world0);
    newworld = copy.(world0)
    if mode == "std"
        update_rates! = update_rates_std!
    elseif mode == "2D"
        update_rates! = update_rates_2D!
    elseif mode == "grad2D"
        update_rates! = update_rates_grad2D!
    elseif mode == "mountain"
        update_rates! = update_rates_mountain!
    elseif mode == "split"
        update_rates! = update_rates_std_split!
    elseif mode == "graph"
        update_rates! = update_rates_graph!
    else
        error("Mode $mode is not recognized")
    end
    for i in 1:(Int(p["tend"])-1)
        # we have to take views, otherwise it does not affect worldall
        world = @view worldall[:,i];newworld = @view worldall[:,i+1];
Victor's avatar
Victor committed
184
        updateWorld_WF!(world,newworld,p,update_rates!,Float64(i))
185
    end
186
    return worldall,collect(0:Int(p["tend"]-1))
187
188
end
"""
Victor's avatar
Victor committed
189
    function runWorld_store_G(p,world0)
190
Gillepsie process. Returns a tuple worldall,tspanarray
Victor's avatar
Victor committed
191
192
If specified p["dt_saving"] determines the time step to save the simulation. If not only last time step is saved

193
"""
Victor's avatar
Victor committed
194
function runWorld_store_G(p,world0)
195
196
    # we store the value of world every 100 changes by default
    tspan = zeros(1)
197
    i = 1;j=0;dt = 1.
198
    N=length(world0);
Victor's avatar
Victor committed
199
200
201
    tspanarray = zeros(1);
    dt_saving = haskey(p,"dt_saving") ? p["dt_saving"] : p["tend"] + 1.
    # Array that stores the agents
202
    worldall = reshape(copy.(world0),N,1)
Victor Boussange's avatar
Victor Boussange committed
203
    # we instantiate C as the biggest size it can take
204
205
206
207
208
209
210
211
212
213
214
    update_rates_std!(skipmissing(world0),p,0.)
    while tspan[i]<p["tend"]
        if dt < 0
            throw("We obtained negative time step dt = $dt at event $i")
        elseif count(ismissing,world0) == p["NMax"]
            @info "All individuals have died :("
            break
        elseif count(ismissing,world0) == 0
            @info "We have reached the maximum number of individuals allowed"
            break
        end
Victor's avatar
Victor committed
215
        if  tspan[end] - tspanarray[end] >= dt_saving
216
217
            @info "saving world @ t = $(tspan[i])/ $(p["tend"])"
            j+=1;sw = size(worldall,2);
218
219
            # we use <= because we save at the end of the wile loop
            if sw <= j
220
221
222
223
                # we double the size of worldall
                worldall = hcat(worldall,Array{Missing}(missing,N,sw))
            end
            worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
224
225
            push!(tspanarray,tspan[i])
        end
Victor's avatar
Victor committed
226
        dt = updateWorld_G!(world0,p,update_rates_std!,tspan)
227
        push!(tspan, tspan[end] + dt)
228
229
        i += 1
    end
230
231
232
233
    # Saving laste time step
    worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
    push!(tspanarray,tspan[i])
    @info "simulation stopped at t=$(tspanarray[end])"
234
235
    return worldall[:,1:j],tspanarray
end