To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit c594f84d authored by Victor's avatar Victor
Browse files

run! not crashing.

now we need to unit test gillepsiee.jl, and especially the update rates
almost there!
parent 1c79e98f
......@@ -29,9 +29,9 @@ module ABMEv
get_thist,get_geo,get_b,get_d,increment_x!,get_inc_reflected,world2df,
split_move,split_merge_move,tin,new_world_G
export World,parameters,time,space,agents,size,maxsize,addAgent!,removeAgent!
export runWorld!,give_birth,updateWorld!,update_clock!,updateBirthEvent!,
export run!,give_birth,updateWorld!,update_clock!,updateBirthEvent!,
updateDeathEvent!#,runWorld_G!,runWorld_WF!,
export Simulation,add_entry!
export Simulation,add_entry!,get_tend,get_size
# export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div,
# get_dist_hist,get_pairwise_average_isolation,
# get_local_pairwise_average_isolation,
......
......@@ -104,7 +104,7 @@ Returns trait i of the agent
Base.getindex(a::Agent,i::Int) = a.x_history[end][i]
Base.getindex(a::Agent,I::Vararg{Int, 2}) = a.x_history[end][I]
Base.getindex(a::Agent,I::Vararg{Int}) = a.x_history[end][I]
get_x(a::Agent) = a.x_history[end]
@deprecate get_x(a) a[:]
......
......@@ -2,7 +2,7 @@ abstract type AbstractAlg end
# this used to be world
mutable struct Simulation{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number,F}
agentarray::Vector{AbstractAgentM}
agentarray::Array{AbstractAgentM,2}
space::S
tspan::Vector{T}
cb::F
......@@ -22,29 +22,36 @@ function Simulation(w0::World{A,S,T};cb=(names = String[],agg =nothing)) where {
agentarray = copy.(collect(w0.agents))
agentarray = hcat(agentarray,Array{Missing}(missing,NMax,1))
!isnothing(cb.agg) ? df_agg = [Dict(cb.names .=> cb.agg(world))] : df_agg = nothing
Simulation{A,S,T,cb.agg}([agentarray],space(w0),tspan,cb.agg,df_agg,parameters(w0))
Simulation{A,S,T,typeof(cb.agg)}(agentarray,space(w0),tspan,cb.agg,df_agg,parameters(w0))
end
get_tend(s::Simulation) = s.tspan[end]
get_size(s::Simulation) = length(s.tspan)
Base.getindex(s::Simulation,i,j) = s.agentarray[i,j]
import Base:lastindex,size
Base.lastindex(s::Simulation,i) = lastindex(s.agentarray,i)
Base.size(s::Simulation,i) = size(s.agentarray,i)
"""
$(SIGNATURES)
Add `w` with callbacks `s.cb` to `s`
"""
function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F<:Function}
push!(s.agentarray,copy.(collect(w.agents)))
push!(s.tspan,w.t)
push!(s.df_agg,Dict(s.cb.names .=> s.cb.agg(world)))
end
# TODO: define two functions with signatures
# function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F<:Function}
# function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F<:Nothing}
"""
$(SIGNATURES)
Add `w` to `s`
Add `w` with callbacks `s.cb` to `s` if provided
"""
function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F<:Nothing}
push!(s.agentarray,copy.(collect(w.agents)))
function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
i = get_size(s)
j = size(s.agentarray,2)
if i == j
# we double the siwe of agent array
s.agentarray = hcat(s.agentarray,Array{Missing}(missing,maxsize(w),j))
end
s.agentarray[:,i+1] .= copy.(collect(w.agents))
push!(s.tspan,w.t)
if F <: Function
push!(s.df_agg,Dict(s.cb.names .=> s.cb.agg(world)))
end
end
#TODO: code it
......
......@@ -6,11 +6,11 @@ This standard updates takes
- carrying capacity of the form K(x)
"""
function update_rates_std!(w::World)
@unpack b,k = parameters(w)
@unpack b,d = parameters(w)
_agents = agents(w)
traits = get_x.(_agents)
# traits = get_xhist.(world)
n = size(world)
n = size(w)
D = zeros(n)
# Here you should do a shared array to compute in parallel
for i in 1:(n-1)
......@@ -21,7 +21,7 @@ function update_rates_std!(w::World)
end
end
# Here we can do it in parallel as well
for (i,a) in _agents
for (i,a) in enumerate(_agents)
a.d = D[i]
a.b = b(traits[i])
end
......@@ -79,26 +79,25 @@ first corresponding to initial conditions and last corresponding to world in the
>:warning: if you choose `nagents = 1` then nothing is saved until the end of the simulation.
"""
# function run(w::World{AbstractAgent{A,R},S,T},g::G;dt_saving=nothing,callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
function run(w::World{AbstractAgent{A,R},S,T},g::G,tend::Number;
function run!(w::World{A,S,T},g::G,tend::Number;
dt_saving=nothing,
callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
@unpack tend = parameters(w)
cb=(names = String[],agg =nothing)) where {A,S,T,G<:Gillepsie}
n=size(w);
NMax = sizemax(w)
NMax = maxsize(w)
t = .0
i = 1;j=1;dt = 1.
isnothing(dt_saving) ? dt_saving = p["tend"] + 1. : nothing
sim = Simulation(w,callbacks=callbacks)
if R <: Rates{true}
update_rates!(w)
i = 1;j=1;dt = 0.
isnothing(dt_saving) ? dt_saving = tend + 1. : nothing
sim = Simulation(w,cb=cb)
if A <: AbstractAgent{AA,Rates{true}} where {AA}
update_rates_std!(w)
end
while t<tend
if dt < 0
throw("We obtained negative time step dt = $dt at event $i")
elseif size(world) == NMax
elseif size(w) == NMax
@info "All individuals have died :("
break
elseif size(world) == 0
elseif size(w) == 0
@info "We have reached the maximum number of individuals allowed"
break
end
......@@ -107,10 +106,11 @@ function run(w::World{AbstractAgent{A,R},S,T},g::G,tend::Number;
add_entry!(s,w)
end
dt = updateWorld!(w,g)
t += dt
i += 1
end
# Saving last time step
add_entry!(sim,w)
@info "simulation stopped at t=$(tspanarray[end]), after $(i) generations"
@info "simulation stopped at t=$(t), after $(i) generations"
return sim
end
......@@ -11,7 +11,7 @@ end
#constructor
function World(w::Vector{A},s::S,p::Dict,t::T=0.) where {A<:AbstractAgent,S<:AbstractSpacesTuple,T}
if typeof(p["D"]) != eltype(skipmissing(w)[1])
throw(ArgumentError("Diffusion coefficient does not match with underlying space"))
throw(ArgumentError("Diffusion coefficient does not match with underlying space\n `D::Tuple`"))
end
ww = vcat(w,repeat([missing],Int(p["NMax"] - length(w))))
World{A,S,T}(ww,s,p,t)
......@@ -60,7 +60,7 @@ Base.getindex(w::World,i::Integer) = getindex.(agents(w),i)
# If `trait = 0` , we return the geotrait.
# """
# get_x(w::World,t::Number,trait::Integer) = trait > 0 ? w[i] : reshape(hcat(get_geo.(w,t)),size(w,1),size(w,2))
#
#
# """
# $(SIGNATURES)
# Returns every traits of every agents of world in the form of an array
......
......@@ -14,12 +14,12 @@ function give_birth(mum_idx::Int,w::World)
return new_a
end
function updateBirthEvent!(world::World,::Gillepsie,mum_idx::Int)
function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int)
# updating competition only the two columns corresponding to agent idx
@unpack d,b = parameters(world)
offspring = give_birth(w,mum_idx)
@unpack d,b = parameters(w)
offspring = give_birth(mum_idx,w)
x_offspring = get_x(offspring)
_agents = agents(world)
_agents = agents(w)
for a in _agents
a.d += d(get_x(a),x_offspring)
end
......@@ -65,8 +65,7 @@ function updateWorld!(w::World{A,S,T},g::G) where {A,S,T,G <: Gillepsie}
else
# birth event
# i_event - n is also the index of the individual to give birth in the world_alive
mum = w[i_event-n]
updateBirthEvent!(w,g,idx_offspring)
updateBirthEvent!(w,g,i_event-n)
end
return dt
else
......
cd(@__DIR__)
using Random
Random.seed!(0)
import ABMEv:update_rates_std!
using LightGraphs
using Test
using Revise,ABMEv
using UnPack,JLD2
a = 0;
myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
K(X) = gaussian(X[1],0.,sigma_K)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 1.5,
"NMax" => Int(10000))
na_init = K0
world0 = new_world_G(na_init,p_default,spread = .01)
worldall,p_default["tspan"] = runWorld_store_G(p_default,world0)
world_alive_test = clean_world(worldall[:,end])
b(X) = gaussian(X[1],0.,sigma_K)
d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
D = (1e-2,)
mu = [.1]
NMax = 10000
tend = 1.5
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
myagents = [Agent(myspace,ancestors=true,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
w1 = copy(w0)
sim = run!(w0,Gillepsie(),tend)
world_alive_test = collect(skipmissing(sim[:,end]))
# @save "gillepsie_test.jld2" world_alive
@load "gillepsie_test.jld2" world_alive
## Testing
@testset "Gillepsie Algorithm" begin
@testset "Testing global functioning" begin
@test size(worldall,2) > 1
@test p_default["tspan"][end] >= p_default["tend"]
@test size(sim,2) > 1
@test get_tend(sim) >= tend
end
## Comparing simulation
@testset "Matching new vs old results " begin
......
......@@ -10,6 +10,19 @@ D = (Int64(1),Float64(1.))
mu = [1,1]
NMax = 100
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
w = World(myagents,myspace,p,0.)
w0 = World(myagents,myspace,p,0.)
s = Simulation(w)
# Basic test
s = Simulation(w0)
typeof(s)
@test get_tend(s) 0.
@test get_size(s) 1
# adding an entry to sim
newa = give_birth(1,w0)
addAgent!(w0,newa)
@test isnothing(add_entry!(s,w0))
@test get_size(s) == 2
@test size(s.agentarray,2) == 2
#TODO: try with callbacks
Markdown is supported
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