Commit 2e9a6662 authored by Victor Boussange's avatar Victor Boussange
Browse files

adding Project.toml, as well as some line to make sure time step is never negative

parent 4c26f995
name = "ABMEv"
uuid = "837ac870-fb52-4b0c-9a0e-030f2f36f5ed"
authors = ["Victor Boussange "]
version = "0.1.0"
......@@ -52,21 +52,25 @@ function updateWorld_G!(world,C,p,update_rates!,tspan,reflected)
# Total rate of events
U = sum(get_d.(world_alive) .+ get_b.(world_alive))
dt = - log(rand(Float64))/U
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
# DEATH EVENT
idx_offspring = idx_world[i_event]
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)
if dt > 0.
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
# DEATH EVENT
idx_offspring = idx_world[i_event]
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)
end
return dt
else
return -1.
end
return dt
end
......@@ -189,7 +189,7 @@ Gillepsie process. Returns a tuple worldall,tspanarray
function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
# we store the value of world every 100 changes by default
tspan = zeros(1)
i = 1;j=0;
i = 1;j=0;dt = 1.
N=length(world0);
tspanarray = [];
ninit = Int(length(world0) - count(ismissing,world0))
......@@ -200,7 +200,7 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
# 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"] && count(ismissing,world0) < p["NMax"] && count(ismissing,world0) > 0
while tspan[i]<p["tend"] && dt > 0. && count(ismissing,world0) < p["NMax"] && count(ismissing,world0) > 0
# we save every ninit times
if mod(i,ninit) == 1
@info "saving world @ t = $(tspan[i])/ $(p["tend"])"
......@@ -212,7 +212,8 @@ 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
push!(tspan, tspan[end] + updateWorld_G!(world0,C,p,update_rates_std!,tspan,reflected))
dt = updateWorld_G!(world0,C,p,update_rates_std!,tspan,reflected)
push!(tspan, tspan[end] + dt)
i += 1
end
return worldall[:,1:j],tspanarray
......
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