Commit 41f47779 authored by Victor Boussange's avatar Victor Boussange
Browse files

significantly modifying ABMEv_Gillepse

parent 4325917f
......@@ -11,33 +11,23 @@ function give_birth(a::Agent,p::Dict,reflected)
return new_a
end
function update_afterbirth_std!(world,C,idx::Int,p::Dict) where T
traits = get_x.(world)
function update_afterbirth_std!(world,idx_offspring,p::Dict) where T
# updating competition only the two columns corresponding to agent idx
α = p["alpha"];K=p["K"]
for i in 1:length(traits)
C[i,idx] = α(traits[i],traits[idx])
C[idx,i] = C[i,idx]
end
# updating death rate only the two columns corresponding to agent idx
for (i,a) in enumerate(world)
a.d += C[i,idx]
α = p["alpha"];K=p["K"];
x_offspring = get_x(world[idx_offspring])
for a in skipmissing(world)
a.d += α(get_x(a),x_offspring)
end
# Now updating new agent
world[idx].d = sum(C[idx,:])
world[idx].b = K(traits[idx])
world[idx_offspring].d = sum(α.(get_x.(skipmissing(world)),Ref(x_offspring)))
world[idx_offspring].b = K(x_offspring)
end
function update_afterdeath_std!(world,C,idx::Int,p::Dict) where T
traits = get_x.(world)
function update_afterdeath_std!(world,x_death,p::Dict) where T
α = p["alpha"]
# updating death rate only the two columns corresponding to agent idx
for (i,a) in enumerate(world)
a.d -= C[i,idx]
end
# updating competition only the two columns corresponding to agent idx
for i in 1:length(traits)
C[i,idx] = .0
C[idx,i] = .0
for a in skipmissing(world)
a.d -= α(get_x(a),x_death)
end
end
......@@ -46,7 +36,7 @@ end
Updating rule for gillepsie setting.
Returning dt drawn from an exponential distribution with parameter the total rates of events.
"""
function updateWorld_G!(world,C,p,update_rates!,tspan,reflected)
function updateWorld_G!(world,p,update_rates!,tspan,reflected)
# total update world
world_alive = skipmissing(world)
idx_world = collect(eachindex(world_alive))
......@@ -57,18 +47,21 @@ function updateWorld_G!(world,C,p,update_rates!,tspan,reflected)
events_weights = ProbabilityWeights(vcat(get_d.(world_alive),get_b.(world_alive)))
i_event = sample(events_weights)
# This is ok since length is a multiple of 2
I = Int(length(events_weights) / 2)
if i_event <= I
N = length(world) - count(ismissing,world)
if i_event <= N
# DEATH EVENT
# In this case i_event is also the index of the individual to die in the world_alive
idx_offspring = idx_world[i_event]
x_death = get_x(world[idx_offspring])
update_afterdeath_std!(world,x_death,p)
world[idx_offspring] = missing
update_afterdeath_std!(world_alive,C,idx_offspring,p)
else
# birth event
idx_offspring = findfirst(ismissing,world)
world[idx_offspring] = give_birth(world[idx_world[i_event-I]],p,reflected)
update_afterbirth_std!(world_alive,C,idx_offspring,p)
# i_event - N is also the index of the individual to give birth in the world_alive
mum = world[idx_world[i_event-N]]
world[idx_offspring] = give_birth(mum,p,reflected)
update_afterbirth_std!(world,idx_offspring,p)
end
return dt
else
......
......@@ -6,26 +6,22 @@ This standard updates takes
- competition kernels of the form α(x,y) and
- carrying capacity of the form K(x)
"""
function update_rates_std!(world,C,p::Dict,t::Float64)
function update_rates_std!(world,p::Dict,t::Float64)
α = p["alpha"];K=p["K"];
# 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)
for i in 1:(N-1)
for j in i+1:N
C[i,j] = α(traits[i],traits[j])
C[j,i] = C[i,j]
C = α(traits[i],traits[j])
world[i].d += C
world[j].d += C
end
end
# 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,:])
# /!| not ideal to assign at every time step the birth rate that is constant
a.b = K(traits[i])
end
end
......@@ -205,9 +201,17 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
# length of worldall should be changed at some point
worldall = reshape(copy.(world0),N,1)
# we instantiate C as the biggest size it can take
C = SharedArray{Float64}((N,N))
update_rates_std!(skipmissing(world0),C,p,0.)
while tspan[i]<p["tend"] && dt > 0. && count(ismissing,world0) < p["NMax"] && count(ismissing,world0) > 0
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
# we save every ninit times
if mod(i,ninit) == 1
@info "saving world @ t = $(tspan[i])/ $(p["tend"])"
......@@ -219,9 +223,10 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
push!(tspanarray,tspan[i])
end
dt = updateWorld_G!(world0,C,p,update_rates_std!,tspan,reflected)
dt = updateWorld_G!(world0,p,update_rates_std!,tspan,reflected)
push!(tspan, tspan[end] + dt)
i += 1
end
@info "simulation stopped at t=$(tspan[end])"
return worldall[:,1:j],tspanarray
end
using ABMEv
a = 0;
sigma_K = .9;
sigma_a = .7;
K(X) = gaussian(X[1],0.,sigma_K)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/1000
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2], # we let mutation constant and =1e-3
"mu" => [.1], # we consider there is a probability of 1 to have mutation
"tend" => 150.,
"NMax" => Int(10000))
na_init = 150
world0 = new_world_G(na_init,p_default,spread = .01, offset = -.25)
@time worldall,p_default["tspan"] = runWorld_store_G(p_default,world0)
# ======================================================================
using JLD2
@save "gillepsie_test.jld2" worldall,p_default
using Plots
Plots.plot(worldall,p_default)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment