Commit bc0a5c8a authored by Victor Boussange's avatar Victor Boussange
Browse files

worldall is not allocated at every time steps anymore

parent cb5c5990
using Distributed;addprocs()
@everywhere push!(LOAD_PATH,homedir()*"/ETHZ/projects/ABMEv.jl/src")
@everywhere using ABMEv,BenchmarkTools,SharedArrays
p= Dict("K0" => 1000.,
"D" => [1e-2 - 1e-3],
"mu" => [.1],
"sigma_a" => [5e-1],
"sigma_K" => [1e0],
"n_alpha" => 2.,
"n_K" => 2.,
"tend" => 5.,
"NMax" => Int(1e4))
na_init = 500
agent0 = [Agent([ .01 * randn() .- .5]) for i in 1:na_init]
world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - na_init)))
C = SharedArray{Float64}((Int(p["K0"]),Int(p["K0"])))
@btime update_afterbirth_std!(skipmissing(world0),C,1,p)
......@@ -16,7 +16,8 @@ module ABMEv
export update_rates!
export Agent,get_fitness,get_x,get_xarray,get_xhist,
get_geo,get_b,get_d,increment_x!,get_inc_reflected,
split_move,split_merge_move,KK,tin
split_move,split_merge_move,KK,tin,new_world_G
export copy,runWorld_store_WF,runWorld_store_G #,runWorld_G!,runWorld_WF!,
export H_discrete,findclusters,var,covgeo,hamming
export update_afterbirth_std!,update_afterdeath_std!
end
......@@ -15,6 +15,18 @@ import Base.copy
copy(a::Agent) = Agent(a.x_history,a.d,a.b)
copy(m::Missing) = missing
"""
function new_world_G(nagents::Int,p::Dict; spread = 1., offset = 0.)
Returns an array of type Array{Union{Missing,Agent}} initialised with normal distribution.
Only relevant for Gillepsie algorithm as of now.
"""
function new_world_G(nagents::Int,p::Dict; spread = 1., offset = 0.)
typeof(spread) <: Array ? spread = spread[:] : nothing;
typeof(offset) <: Array ? offset = offset[:] : nothing;
agent0 = [Agent( spread .* randn(length(spread)) .+ offset) for i in 1:nagents]
world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - nagents)))
return world0
end
# returns trait i of the agent
get_x(a::Agent,i::Number) = a.x_history[Int(i):Int(i),end]
......
......@@ -8,7 +8,6 @@ Used for Gillepsie setting
function give_birth(a::Agent,p::Dict,reflected)
new_a = copy(a)
increment_x!(new_a,p,reflected=reflected)
# ! carefull you also need to change get_xhist
return new_a
end
......@@ -25,7 +24,7 @@ function update_afterbirth_std!(world,C,idx::Int,p::Dict) where T
end
# Now updating new agent
world[idx].d = sum(C[idx,:]) / p["K0"]
world[idx].b = K(traits[idx][:,end],1,p["n_K"],p["sigma_K"])
world[idx].b = K(traits[idx][:,end],1.,p["n_K"],p["sigma_K"])
end
function update_afterdeath_std!(world,C,idx::Int,p::Dict) where T
......
......@@ -19,7 +19,7 @@ function update_rates_std!(world,C,p::Dict,t::Float64)
# this could be imptrove since \alpha is symmetric, by using a symmetric matrix
a.d = sum(C[i,:]) / p["K0"]
# /!| not ideal to assign at every time step the birth rate that is constant
a.b = K(traits[i][:,end],1,p["n_K"],p["sigma_K"])
a.b = K(traits[i][:,end],1.,p["n_K"],p["sigma_K"])
end
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=1;
N=length(world0);
tspanarray = [];
ninit = Int(length(world0) - count(ismissing,world0))
......@@ -197,15 +197,19 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
# We decide that we store agents ninit events where ninit is the initial population
# 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"] && 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"])"
j+=1;
# worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
# we are reallocating worldall everytime which is bad
worldall = hcat(worldall,copy.(world0))
@info "saving world @ t = $(tspan[i])/ $(p["tend"])"
j+=1;sw = size(worldall,2);
if sw < j
# 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)));
push!(tspanarray,tspan[i])
end
push!(tspan, tspan[end] + updateWorld_G!(world0,C,p,update_rates_std!,tspan,reflected))
......
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