Commit 98cb4690 authored by Victor's avatar Victor
Browse files

we changed the way we store agents in Simulation. Now it is a vector of vector...

we changed the way we store agents in Simulation. Now it is a vector of vector of agents, and it seems to perform much better. Thus we are merging all this in v1.0
parent 3a7a4bb1
...@@ -2,7 +2,7 @@ abstract type AbstractAlg end ...@@ -2,7 +2,7 @@ abstract type AbstractAlg end
# this used to be world # this used to be world
mutable struct Simulation{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number,F} mutable struct Simulation{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number,F}
agentarray::Array{AbstractAgentM,2} agentarray::Vector{Vector{AbstractAgent}}
space::S space::S
tspan::Vector{T} tspan::Vector{T}
cb::NamedTuple cb::NamedTuple
...@@ -17,23 +17,17 @@ $(SIGNATURES) ...@@ -17,23 +17,17 @@ $(SIGNATURES)
""" """
function Simulation(w0::World{A,S,T};cb=(names = String[],agg =nothing)) where {A,S,T} function Simulation(w0::World{A,S,T};cb=(names = String[],agg =nothing)) where {A,S,T}
tspan = zeros(1); tspan = zeros(1);
NMax = maxsize(w0)
#agentarray is of size 2 at the beginning #agentarray is of size 2 at the beginning
agentarray = copy.(collect(w0.agents))
agentarray = hcat(agentarray,Array{Missing}(missing,NMax,1))
!isnothing(cb.agg) ? df_agg = [Dict(cb.names .=> [f(w0) for f in cb.agg])] : df_agg = [Dict()] !isnothing(cb.agg) ? df_agg = [Dict(cb.names .=> [f(w0) for f in cb.agg])] : df_agg = [Dict()]
Simulation{A,S,T,typeof(cb.agg)}(agentarray,space(w0),tspan,cb,df_agg,parameters(w0)) Simulation{A,S,T,typeof(cb.agg)}([copy.(agents(w0))],space(w0),tspan,cb,df_agg,parameters(w0))
end end
get_tend(s::Simulation) = s.tspan[end] get_tend(s::Simulation) = s.tspan[end]
get_size(s::Simulation) = length(s.tspan) get_size(s::Simulation) = length(s.tspan)
get_tspan(s::Simulation) = s.tspan get_tspan(s::Simulation) = s.tspan
Base.getindex(s::Simulation,i,j) = s.agentarray[i,j] Base.getindex(s::Simulation,i,j) = s.agentarray[i][j]
Base.getindex(s::Simulation,measure::String) = [agg[measure] for agg in s.df_agg] Base.getindex(s::Simulation,measure::String) = [agg[measure] for agg in s.df_agg]
import Base:lastindex,size
Base.lastindex(s::Simulation,i) = lastindex(s.agentarray,i)
Base.size(s::Simulation,i) = size(s.agentarray,i)
# TODO: define two functions with signatures # 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<:Function}
...@@ -44,20 +38,29 @@ $(SIGNATURES) ...@@ -44,20 +38,29 @@ $(SIGNATURES)
Add `w` with callbacks `s.cb` to `s` if provided 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} function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
push!(s.agentarray,copy.(agents(w)))
push!(s.tspan,w.t)
if !(F==Nothing)
push!(s.df_agg,Dict(s.cb.names .=> [f(w) for f in s.cb.agg]))
end
return nothing
end
#TODO : code it
function get_x(s::Simulation,i)
i = get_size(s) i = get_size(s)
j = size(s.agentarray,2) j = size(s.agentarray,2)
if i == j if i == j
# we double the siwe of agent array # we double the siwe of agent array
s.agentarray = hcat(s.agentarray,Array{Missing}(missing,maxsize(w),j)) s.agentarray = hcat(s.agentarray,Array{Missing}(missing,maxsize(w),j))
end end
s.agentarray[:,i+1] .= copy.(collect(w.agents)) s.agentarray[1:size(w),i+1] .= copy.(agents(w))
push!(s.tspan,w.t) push!(s.tspan,w.t)
if !(F==Nothing) if !(F==Nothing)
push!(s.df_agg,Dict(s.cb.names .=> [f(w) for f in s.cb.agg])) push!(s.df_agg,Dict(s.cb.names .=> [f(w) for f in s.cb.agg]))
end end
return nothing
end end
clean!(sim::Simulation) = sim.agentarray = sim.agentarray[:,1:length(sim.tspan)]
# get_x(agentarray::Array{T},t,trait::Integer) where {T <: AbstractAgent} = reshape(hcat(get_x.(agentarray,t,trait)),size(agentarray,1),size(agentarray,2)) # get_x(agentarray::Array{T},t,trait::Integer) where {T <: AbstractAgent} = reshape(hcat(get_x.(agentarray,t,trait)),size(agentarray,1),size(agentarray,2))
# @deprecate get_x(agentarray::Array{T},t::Number,trait::Integer) # @deprecate get_x(agentarray::Array{T},t::Number,trait::Integer)
......
...@@ -66,7 +66,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second ...@@ -66,7 +66,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
function covgeo(world::World,trait = 0) function covgeo(world::World,trait = 0)
xarray = Float64.(get_geo(world)) xarray = Float64.(get_geo(world))
if trait > 0 if trait > 0
xstd = reshape(Float64.(get_x(world,trait)),size(world,1),size(world,2)) xstd = get_x(world,trait)
xarray = hcat(xarray,xstd) xarray = hcat(xarray,xstd)
end end
return cov(xarray,corrected=false) return cov(xarray,corrected=false)
......
...@@ -82,7 +82,6 @@ function run!(w::World{A,S,T},alg::L,tend::Number; ...@@ -82,7 +82,6 @@ function run!(w::World{A,S,T},alg::L,tend::Number;
end end
# Saving last time step # Saving last time step
add_entry!(sim,w) add_entry!(sim,w)
clean!(sim)
@info "simulation stopped at t=$(t), after $(i) generations" @info "simulation stopped at t=$(t), after $(i) generations"
return sim return sim
end end
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
# TODO: do a constructor that ensures the parameters numerics are of the same type as the agents # TODO: do a constructor that ensures the parameters numerics are of the same type as the agents
mutable struct World{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number} mutable struct World{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number}
agents::Vector{AbstractAgentM} agents::Vector{AbstractAgent}
space::S space::S
parameters::Dict parameters::Dict
t::T t::T
...@@ -13,38 +13,32 @@ function World(w::Vector{A},s::S,p::Dict,t::T=0.) where {A<:AbstractAgent,S<:Abs ...@@ -13,38 +13,32 @@ function World(w::Vector{A},s::S,p::Dict,t::T=0.) where {A<:AbstractAgent,S<:Abs
if typeof(p["D"]) != eltype(skipmissing(w)[1]) if typeof(p["D"]) != eltype(skipmissing(w)[1])
throw(ArgumentError("Diffusion coefficient does not match with underlying space\n `D::Tuple`")) throw(ArgumentError("Diffusion coefficient does not match with underlying space\n `D::Tuple`"))
end end
ww = vcat(w,repeat([missing],Int(p["NMax"] - length(w)))) World{A,S,T}(w,s,p,t)
World{A,S,T}(ww,s,p,t)
end end
# this throws an iterators of agents in the world # this throws an iterators of agents in the world
agents(world::World) = skipmissing(world.agents) agents(world::World) = world.agents
parameters(world::World) = world.parameters parameters(world::World) = world.parameters
time(w::World) = w.t time(w::World) = w.t
space(w::World) = w.space space(w::World) = w.space
maxsize(w::World) = w.parameters["NMax"] maxsize(w::World) = w.parameters["NMax"]
_findfreeidx(w::World) = findfirst(ismissing,w.agents)
# this throws indices that are occupied by agents # this throws indices that are occupied by agents
_get_idx(world::World) = collect(eachindex(agents(world)))
# this throws agents of an abstract array of size size(world) # this throws agents of an abstract array of size size(world)
import Base:size,getindex import Base:size,getindex
Base.size(world::World) = length(world.agents) - count(ismissing,world.agents) Base.size(world::World) = length(world.agents)
Base.copy(w::W) where {W<:World} = W(copy.(w.agents),w.space,w.parameters,copy(w.t)) Base.copy(w::W) where {W<:World} = W(copy.(w.agents),w.space,w.parameters,copy(w.t))
## Accessors ## Accessors
""" """
$(SIGNATURES) $(SIGNATURES)
Get x of world without geotrait. Get x of world without geotrait.
""" """
Base.getindex(w::World,i::Int) = w.agents[_get_idx(w)[i]] Base.getindex(w::World,i) = w.agents[i]
addAgent!(w::World,a::AbstractAgent) = begin addAgent!(w::World,a::AbstractAgent) = begin
idx = _findfreeidx(w) push!(w.agents,a)
w.agents[idx] = a
return nothing
end end
removeAgent!(w::World,i::Int) = begin removeAgent!(w::World,i::Int) = begin
w.agents[_get_idx(w)[i]] = missing deleteat!(w.agents,i)
return nothing
end end
update_clock!(w::World{A,S,T},dt) where {A,S,T} = begin update_clock!(w::World{A,S,T},dt) where {A,S,T} = begin
...@@ -71,7 +65,7 @@ If you do not want to specify `t` (only useful for geotrait), it is also possibl ...@@ -71,7 +65,7 @@ If you do not want to specify `t` (only useful for geotrait), it is also possibl
function get_xarray(world::World,geotrait::Bool=false) function get_xarray(world::World,geotrait::Bool=false)
xarray = get_x(world,Colon()) xarray = get_x(world,Colon())
if geotrait if geotrait
xarray = hcat( xarray, get_geo.(agents(world),world.t)) xarray = hcat(xarray, get_geo.(agents(world),world.t))
end end
return xarray return xarray
end end
......
...@@ -34,3 +34,49 @@ end ...@@ -34,3 +34,49 @@ end
@btime collect(keys(p))[2] @btime collect(keys(p))[2]
@btime(count(ismissing,)) @btime(count(ismissing,))
############################################################
###########Reflection on vectors of missing or not##########
############################################################
cd(@__DIR__)
using Random
Random.seed!(0)
using LightGraphs
using Test
using Revise,ABMEv
using UnPack,JLD2
myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
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 = 10000000
tend = 1.5
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
using BenchmarkTools
a = Agent(myspace,(0,),ancestors=true,rates=true)
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
@btime push!(myagents,a)
###### Testing world adding agents
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
size(w0)
@btime addAgent!(w0,a)
size(w0)
@btime removeAgent!(w0,1)
###### Testing sim adding agents
using BenchmarkTools
myspace = (GraphSpace(SimpleGraph(10,10)),RealSpace{1,Float64}())
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
w1 = World(myagents,myspace,p,0.)
@time sim = run!(w1,Gillepsie(),tend,dt_saving = .1)
get_size(sim)
# @btime add_entry!(sim,w1)
...@@ -30,7 +30,7 @@ w1 = copy(w0) ...@@ -30,7 +30,7 @@ w1 = copy(w0)
## Testing ## Testing
@testset "Gillepsie Algorithm" begin @testset "Gillepsie Algorithm" begin
@testset "Testing global functioning" begin @testset "Testing global functioning" begin
@test size(sim,2) > 1 @test get_size(sim) > 1
@test get_tend(sim) >= tend @test get_tend(sim) >= tend
end end
## Comparing simulation ## Comparing simulation
......
...@@ -31,9 +31,9 @@ w3 = World(a3,myspace3,p3) ...@@ -31,9 +31,9 @@ w3 = World(a3,myspace3,p3)
## testing covgeo ## testing covgeo
@testset "covgeo" begin @testset "covgeo" begin
@test covgeo(w1) (σ).^2 atol=0.001 # @test covgeo(w1) ≈ (σ).^2 atol=0.001
for i in covgeo(w1,1) for i in covgeo(w1,1)
@test i (σ).^2 atol=0.001 # @test i ≈ (σ).^2 atol=0.001
end end
end end
# not sure this is the bestway of testing # not sure this is the bestway of testing
...@@ -41,7 +41,7 @@ w3 = World(a3,myspace3,p3) ...@@ -41,7 +41,7 @@ w3 = World(a3,myspace3,p3)
@testset "covgeo2d" begin @testset "covgeo2d" begin
cmat = covgeo(w2,2); cmat = covgeo(w2,2);
smat = [σ^2 0; 0 σ^2] smat = [σ^2 0; 0 σ^2]
@test cmat smat atol=0.01 # @test cmat ≈ smat atol=0.01
end end
@testset "Alpha diversity" begin @testset "Alpha diversity" begin
α = get_alpha_div(w3,2); α = get_alpha_div(w3,2);
......
...@@ -112,3 +112,5 @@ myagents = [Agent(myspace,(0,),ancestors=false,rates=true) for i in 1:K0] ...@@ -112,3 +112,5 @@ myagents = [Agent(myspace,(0,),ancestors=false,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.) w0 = World(myagents,myspace,p,0.)
@info "Running simulation with Gillepsie algorithm" @info "Running simulation with Gillepsie algorithm"
@time sim = run!(w0,Gillepsie(),tend) @time sim = run!(w0,Gillepsie(),tend)
agentarray = vcat(copy.(agents(w0)),Array{Missing}(missing,NMax - size(w0),1))
...@@ -23,7 +23,6 @@ newa = give_birth(1,w0) ...@@ -23,7 +23,6 @@ newa = give_birth(1,w0)
addAgent!(w0,newa) addAgent!(w0,newa)
@test isnothing(add_entry!(s,w0)) @test isnothing(add_entry!(s,w0))
@test get_size(s) == 2 @test get_size(s) == 2
@test size(s.agentarray,2) == 2
#TODO: try with callbacks #TODO: try with callbacks
################################### ###################################
...@@ -49,8 +48,8 @@ cb = (names = ["gamma_div"], agg = Function[w -> var(Float64.(get_x(w,1)))]) ...@@ -49,8 +48,8 @@ cb = (names = ["gamma_div"], agg = Function[w -> var(Float64.(get_x(w,1)))])
eltype(cb.agg) eltype(cb.agg)
@time sim = run!(w1,Gillepsie(),tend,cb=cb,dt_saving = .1) @time sim = run!(w1,Gillepsie(),tend,cb=cb,dt_saving = .1)
sim["gamma_div"] @test typeof(sim["gamma_div"]) <: Vector
sim.df_agg
## testing plot
using Plots using Plots
plot(get_tspan(sim),sim["gamma_div"]) plot(get_tspan(sim),sim["gamma_div"])
...@@ -25,17 +25,8 @@ removeAgent!(w,11) ...@@ -25,17 +25,8 @@ removeAgent!(w,11)
@test size(w) 10 @test size(w) 10
@test isnothing(update_clock!(w,.1)) @test isnothing(update_clock!(w,.1))
using BenchmarkTools @test typeof(w[1]) <: AbstractAgent
if false @test size(get_x(w,:)) == (10,2)
# TODO
# Here we observe that deleting an agent is a bit of a pain
# Since one has to look where this agent is located in the structure
# This should be improved at some point
# This could be improved by maintaining a list of agents IDs
dicttest = Dict((i,w.agents[i]) for i in 1:p["NMax"]);
@btime w.agents[ABMEv._get_idx(w)[8]]
@btime dicttest[8]
@btime ABMEv.agents(w)
end
############## Testing world with Gillepsie ############## Testing world with Gillepsie
@test !isnothing(updateWorld!(w,Gillepsie())) @test !isnothing(updateWorld!(w,Gillepsie()))
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