Gillepsie algorithm, as well as callbacks, working

src/ABMEv.jl

@@ -31,7 +31,7 @@ module ABMEv
     export World,parameters,time,space,agents,size,maxsize,addAgent!,removeAgent!
     export run!,give_birth,updateWorld!,update_clock!,updateBirthEvent!,
         updateDeathEvent!#,runWorld_G!,runWorld_WF!,
-    export Simulation,add_entry!,get_tend,get_size
+    export Simulation,add_entry!,get_tend,get_size,get_tspan
     # 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,

src/ABMEv_Agent.jl

@@ -102,9 +102,7 @@ Base.summary(A::AbstractAgent) = string(TYPE_COLOR,nameof(typeof(a)),NO_COLOR,"
 Returns trait i of the agent
 """
 
-Base.getindex(a::Agent,i::Int) = a.x_history[end][i]
-
-Base.getindex(a::Agent,I::Vararg{Int}) = a.x_history[end][I]
+Base.getindex(a::Agent,i) = a.x_history[end][i]
 
 get_x(a::Agent) = a.x_history[end]
 @deprecate get_x(a) a[:]
@@ -134,11 +132,11 @@ Get time when agent born.
 get_t(a::Agent) = a.t_history[end]
 
 # TODO: change this with getindex
-get_xhist(a::Agent,i::Number) = [a.x_history[t][Int(i)] for t in 1:length(a.xhistory)]
+get_xhist(a::AbstractAgent,i::Number) = [a.x_history[t][Int(i)] for t in 1:length(a.x_history)]
 
-get_xhist(a::Agent) = a.x_history
+get_xhist(a::AbstractAgent) = a.x_history
 
-get_thist(a::Agent) = a.t_history
+get_thist(a::AbstractAgent) = a.t_history
 
 get_d(a::AbstractAgent) = a.d
 
@@ -146,7 +144,9 @@ get_b(a::AbstractAgent) = a.b
 
 get_fitness(a::AbstractAgent) = a.b - a.d
 
-Base.ndims(a::Agent{A,R,T,U,V}) where {A,R,T<:(Tuple{Vararg{S,N}} where {S,N}),U,V} = N
+# TODO : we can surely extract N in Agent{A,R,Tuple{Vararg{S,N}},U,V}
+# Inspiration : where U <: Union{Missing,Agent{T}} where T
+Base.length(a::AbstractAgent) = length(a.x_history[end])
 
 nancestors(a::Agent) = length(a.x_history)
 

src/ABMEv_Sim.jl

@@ -5,8 +5,8 @@ mutable struct Simulation{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number,F}
     agentarray::Array{AbstractAgentM,2}
     space::S
     tspan::Vector{T}
-    cb::F
-    df_agg
+    cb::NamedTuple
+    df_agg::Vector{Dict}
     p::Dict{String,Any}
 end
 
@@ -21,13 +21,16 @@ function Simulation(w0::World{A,S,T};cb=(names = String[],agg =nothing)) where {
     #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 .=> cb.agg(world))] : df_agg = nothing
-    Simulation{A,S,T,typeof(cb.agg)}(agentarray,space(w0),tspan,cb.agg,df_agg,parameters(w0))
+    !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))
  end
 
 get_tend(s::Simulation) = s.tspan[end]
 get_size(s::Simulation) = length(s.tspan)
+get_tspan(s::Simulation) = s.tspan
 Base.getindex(s::Simulation,i,j) = s.agentarray[i,j]
+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)
@@ -49,8 +52,8 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
     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)))
+    if !(F==Nothing)
+        push!(s.df_agg,Dict(s.cb.names .=> [f(w) for f in s.cb.agg]))
     end
 end
 

src/ABMEv_metrics.jl

@@ -51,8 +51,8 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
 # Arguments
 
 """
-function var(world::Array{U,1};trait=1) where U <: Union{Missing,Agent{T}} where T
-    world = collect(skipmissing(world))
+function var(world::World;trait=1)
+    world = agents(world)
     xarray = get_x(world,trait)
     return var(xarray,dims=1,corrected=false)
 end

src/ABMEv_runworld.jl

-
-"""
-    update_rates_std!(world,p::Dict,t::Float64)
-This standard updates takes
-    - competition kernels of the form α(x,y) and
-    - carrying capacity of the form K(x)
-"""
-function  update_rates_std!(w::World)
-    @unpack b,d = parameters(w)
-    _agents = agents(w)
-    traits = get_x.(_agents)
-    # traits = get_xhist.(world)
-    n = size(w)
-    D = zeros(n)
-    # Here you should do a shared array to compute in parallel
-    for i in 1:(n-1)
-        for j in i+1:n
-            C = d(traits[i],traits[j])
-            D[i] += C
-            D[j] += C
-        end
-    end
-    # Here we can do  it in parallel as well
-    for (i,a) in enumerate(_agents)
-        a.d = D[i]
-        a.b = b(traits[i])
-    end
-end
-
 function  update_rates_graph!(world,C,p::Dict,t::Float64)
     for e in edges(p["g"])
         # agents on e
@@ -79,9 +50,9 @@ 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{A,S,T},g::G,tend::Number;
+function run!(w::World{A,S,T},alg::L,tend::Number;
                 dt_saving=nothing,
-                cb=(names = String[],agg =nothing)) where {A,S,T,G<:Gillepsie}
+                cb=(names = String[],agg =nothing)) where {A,S,T,L<:AbstractAlg}
     n=size(w);
     NMax = maxsize(w)
     t = .0
@@ -89,7 +60,7 @@ function run!(w::World{A,S,T},g::G,tend::Number;
     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)
+        update_rates!(w,alg)
     end
     while t<tend
         if dt < 0
@@ -103,9 +74,9 @@ function run!(w::World{A,S,T},g::G,tend::Number;
         end
         if  t - get_tend(sim) >= dt_saving
             @info "saving world @ t = $(t)/ $(tend)"
-            add_entry!(s,w)
+            add_entry!(sim,w)
         end
-        dt = updateWorld!(w,g)
+        dt = updateWorld!(w,alg)
         t +=  dt
         i += 1
     end

src/ABMEv_world.jl

@@ -17,20 +17,26 @@ function World(w::Vector{A},s::S,p::Dict,t::T=0.) where {A<:AbstractAgent,S<:Abs
     World{A,S,T}(ww,s,p,t)
 end
 
+# this throws an iterators of agents in the world
+agents(world::World) = skipmissing(world.agents)
 parameters(world::World) = world.parameters
 time(w::World) = w.t
 space(w::World) = w.space
+maxsize(w::World) = w.parameters["NMax"]
+_findfreeidx(w::World) = findfirst(ismissing,w.agents)
 # this throws indices that are occupied by agents
-_get_idx(world::World) = collect(eachindex(skipmissing(world.agents)))
+_get_idx(world::World) = collect(eachindex(agents(world)))
 # this throws agents of an abstract array of size size(world)
 import Base:size,getindex
 Base.size(world::World) = length(world.agents) - count(ismissing,world.agents)
+Base.copy(w::W) where {W<:World} = W(copy.(w.agents),w.space,w.parameters,copy(w.t))
+## Accessors
+"""
+$(SIGNATURES)
+Get x of world without geotrait.
+"""
 Base.getindex(w::World,i::Int) = w.agents[_get_idx(w)[i]]
 
-# this throws an iterators of agents in the world
-agents(world::World) = skipmissing(world.agents)
-maxsize(w::World) = w.parameters["NMax"]
-_findfreeidx(w::World) = findfirst(ismissing,w.agents)
 addAgent!(w::World,a::AbstractAgent) = begin
     idx = _findfreeidx(w)
     w.agents[idx] = a
@@ -46,38 +52,25 @@ update_clock!(w::World{A,S,T},dt) where {A,S,T} = begin
     return nothing
 end
 
+#TODO : code it
+"""
+$(SIGNATURES)
+Returns trait of every agents of world in the form of an array which dimensions corresponds to the input.
+If `trait = 0` , we return the geotrait.
+"""
+get_x(w::World,trait) = vcat(getindex.(agents(w),trait)...)
 
-## Accessors
 """
 $(SIGNATURES)
-Get x of world without geotrait.
+Returns every traits of every agents of `world` in the form **of a one dimensional array** (in contrast to `get_x`).
+If `geotrait=true` the geotrait is also added to the set of trait, in the last line.
+If you do not want to specify `t` (only useful for geotrait), it is also possible to use `get_xarray(world::Array{T,1}) where {T <: Agent}`.
 """
-Base.getindex(w::World,i::Integer) = getindex.(agents(w),i)
-#TODO : code it
-# """
-# $(SIGNATURES)
-# Returns trait of every agents of world in the form of an array which dimensions corresponds to the input.
-# 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
-# """
-# function get_xarray(world::Array{T,1}) where {T <: Agent}
-#     return hcat(get_x.(world)...)
-# end
-# """
-# $(SIGNATURES)
-# Returns every traits of every agents of `world` in the form **of a one dimensional array** (in contrast to `get_x`).
-# If `geotrait=true` the geotrait is also added to the set of trait, in the last line.
-# If you do not want to specify `t` (only useful for geotrait), it is also possible to use `get_xarray(world::Array{T,1}) where {T <: Agent}`.
-# """
-# function get_xarray(world::Array{T,1},t::Number,geotrait::Bool=false) where {T <: Agent}
-#     xarray = hcat(get_x.(world)...)
-#     if geotrait
-#         xarray = vcat( xarray, get_geo.(world,t)')
-#     end
-#     return xarray
-# end
+function get_xarray(world::World,geotrait::Bool=false)
+    xarray = hcat(get_x.(agents.(world))...)
+    if geotrait
+        xarray = vcat( xarray, get_geo.(agents(world),world.t)')
+    end
+    return xarray
+end
+@deprecate get_xarray(world,geotrait=false) get_x(world,Colon())

src/algo/ABMEv_Gillepsie.jl

@@ -24,7 +24,7 @@ function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int)
         a.d += d(get_x(a),x_offspring)
     end
     # Now updating new agent
-    offspring.d = sum(d.(get_x.(_agents),Ref(x_offspring))) - d(x_offspring,x_offspring)
+    offspring.d = sum(d.(get_x.(_agents),Ref(x_offspring))) #- d(x_offspring,x_offspring)
     offspring.b = b(x_offspring)
     addAgent!(w,offspring)
 end
@@ -39,6 +39,34 @@ function updateDeathEvent!(world::World,::Gillepsie,i_event::Int)
     end
 end
 
+"""
+    update_rates_std!(world,p::Dict,t::Float64)
+This standard updates takes
+    - competition kernels of the form α(x,y) and
+    - carrying capacity of the form K(x)
+"""
+function  update_rates!(w::World,::Gillepsie)
+    @unpack b,d = parameters(w)
+    _agents = agents(w)
+    traits = get_x.(_agents)
+    # traits = get_xhist.(world)
+    n = size(w)
+    D = zeros(n)
+    # Here you should do a shared array to compute in parallel
+    for i in 1:(n-1)
+        for j in i+1:n
+            C = d(traits[i],traits[j])
+            D[i] += C
+            D[j] += C
+        end
+    end
+    # Here we can do  it in parallel as well
+    for (i,a) in enumerate(_agents)
+        a.d = D[i]
+        a.b = b(traits[i])
+    end
+end
+
 """
     function updateWorld_G!(world,t,p)
 Updating rule for gillepsie setting.

test/benchmark_dict.jl

+p = Dict(1:1:1000 .=> 1:1:1000)
+delete!(p,999)
+a = Any[rand() for i in 1:1000]
+a[2:2:1000] .= missing
+get_idx(w) = collect(eachindex(skipmissing(w)))
+using BenchmarkTools
+@btime a[get_idx(a)[3]]
+@btime p[3]
+
+
+@btime [true for i in 1:1000]
+[true for i in 1:1000]
+
+@btime sort!(values(p))
+function findfreeidx(p)
+    i = 1
+    while haskey(p,i)
+        i+=1
+    end
+    i
+end
+
+function get_xp(a,j)
+    for (i,a) in enumerate(skipmissing(a))
+        if i==j
+            return a
+            break
+        end
+    end
+end
+@btime get_xp(a,998)
+
+@btime findfreeidx(p)
+@btime collect(keys(p))[2]
+
+@btime(count(ismissing,))

test/gillepsie.jl

@@ -15,18 +15,18 @@ d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
 D = (1e-2,)
 mu = [.1]
 NMax = 10000
-tend = 1.5
+tend = 100
 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
 
-myagents = [Agent(myspace,ancestors=true,rates=true) for i in 1:K0]
+myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
 w0 = World(myagents,myspace,p,0.)
 w1 = copy(w0)
+@info "Running simulation with Gillepsie algorithm"
+@time sim = run!(w1,Gillepsie(),tend)
+@test typeof(sim) <: Simulation
 
-sim = run!(w0,Gillepsie(),tend)
-
-world_alive_test = collect(skipmissing(sim[:,end]))
 # @save "gillepsie_test.jld2" world_alive
-@load "gillepsie_test.jld2" world_alive
+# @load "gillepsie_test.jld2" world_alive
 ## Testing
 @testset "Gillepsie Algorithm" begin
         @testset "Testing global functioning" begin
@@ -34,19 +34,35 @@ world_alive_test = collect(skipmissing(sim[:,end]))
                 @test get_tend(sim) >= tend
         end
         ## Comparing simulation
-        @testset "Matching new vs old results " begin
-                xarray = get_x(world_alive,1);xarray_test = get_x(world_alive_test,1);
-                @test xarray ≈ xarray_test
-        end
+        # @testset "Matching new vs old results " begin
+        #         xarray = get_x(w1,1);xarray_test = get_x(world_alive,1);
+        #         @test prod(xarray .≈ xarray_test)
+        # end
 
         @testset "Testing update rates matrix" begin
-                bs_end = get_b.(world_alive);ds_end = get_d.(world_alive)
-                update_rates_std!(world_alive,p_default,0.);
-                bs_recalculated = get_b.(world_alive);ds_recalculated = get_d.(world_alive);
-
-                @test bs_end ≈ bs_recalculated;
+                @testset "birth event" begin
+                        myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
+                        w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie())
+                        mum_idx = 1
+                        updateBirthEvent!(w0,Gillepsie(),1)
+                        bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1))
+                        update_rates!(w1,Gillepsie());
+                        bs_recalculated = get_b.(agents(w1));ds_recalculated = get_d.(agents(w1))
+                        @test prod(bs_end .≈ bs_recalculated)
+                        @test prod(ds_end .≈ ds_recalculated)
 
-                @test ds_end ≈ ds_recalculated;
+                end
 
+                @testset "death event" begin
+                        myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
+                        w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie())
+                        mum_idx = 1
+                        updateDeathEvent!(w0,Gillepsie(),1)
+                        bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1))
+                        update_rates!(w1,Gillepsie());
+                        bs_recalculated = get_b.(agents(w1));ds_recalculated = get_d.(agents(w1))
+                        @test prod(bs_end .≈ bs_recalculated)
+                        @test prod(ds_end .≈ ds_recalculated)
+                end
         end
 end

test/metrics.jl

-
+myspace1 = (RealSpace{1,Float64}(),)
+myspace2 = (RealSpace{2,Float64}(),)
+myspace3 = (DiscreteSegment(Int16(1),Int16(10)),RealSpace{1,Float64}())
 K0 = 1000; σ = 1e-1
-agents1 = [Agent( [σ]  .* randn(1) .- .5) for i in 1:K0]
-agents2 = [Agent( [σ σ]  .* randn(2) .- .5) for i in 1:K0]
-agentsd = [Agent{MixedAgent}( Float16[rand(1:10), 1e-1* randn() + 5.5] ) for i in 1:K0]
+w1 = [Agent(myspace1, (σ,)  .* randn() .- .5) for i in 1:K0]
+w2 = [Agent(myspace2,(tuple((σ, σ)  .* randn(2) .- .5...),)) for i in 1:K0]
+w3 = [Agent(myspace3, (rand(Int16.(1:10)), 1e-1* randn() + 5.5 )) for i in 1:K0]
 p = Dict("mu" => [1. 1.],"D" => [0. 0.], "nodes" =>10 )
 ## testing variance
 @testset "Testing metrics" begin
     @testset "var" begin
-        @test first(var(agents1)) ≈ (σ).^2 atol=0.001
-        @test first(var(agents2,trait=2)) ≈ (σ).^2 atol=0.001
+        @test first(var(w1)) ≈ (σ).^2 atol=0.001
+        @test first(var(w2,trait=2)) ≈ (σ).^2 atol=0.001
     end
 
     ## testing covgeo
     @testset "covgeo" begin
-        @test covgeo(agents1) ≈ (σ).^2 atol=0.001
-         for i in covgeo(agents1,1)
+        @test covgeo(w1) ≈ (σ).^2 atol=0.001
+         for i in covgeo(w1,1)
              @test i ≈ (σ).^2 atol=0.001
          end
      end
      # not sure this is the bestway of testing
      # there is a problem here
      @testset "covgeo2d" begin
-         cmat = covgeo(agents2,2);
+         cmat = covgeo(w2,2);
          smat = [σ^2 0; 0 σ^2]
          @test cmat ≈ smat atol=0.01
       end
       @testset "Alpha diversity" begin
-          α = get_alpha_div(agentsd,1.0,2);
+          α = get_alpha_div(w3,1.0,2);
           @test abs(α) < Inf
        end
        @testset "Beta diversity" begin
-           β = get_beta_div(agentsd,1.0,2);
+           β = get_beta_div(w3,1.0,2);
            @test abs(β) < Inf
        end
        @testset "Isolation by history - hamming distance" begin

test/simulation.jl

@@ -26,3 +26,31 @@ addAgent!(w0,newa)
 @test size(s.agentarray,2) == 2
 
 #TODO: try with callbacks
+###################################
+#######CALLBACKS###################
+###################################
+
+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 = 10000
+tend = 1.5
+p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
+
+myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
+w0 = World(myagents,myspace,p,0.)
+w1 = copy(w0)
+cb = (names = ["gamma_div"], agg = Function[w -> var(Float64.(get_x(w,1)))])
+eltype(cb.agg)
+@time sim = run!(w1,Gillepsie(),tend,cb=cb,dt_saving = .1)
+
+sim["gamma_div"]
+sim.df_agg
+
+using Plots
+plot(get_tspan(sim),sim["gamma_div"])