aggreg_last_dt : added option to save cb at specified time steps

src/ABMEv_Sim.jl

@@ -53,6 +53,15 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
     return nothing
 end
 
+function add_entry_cb_only!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
+    push!(s.agentarray,[Agent(w.space)]) # pushing NA agents
+    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_xnt(s::Simulation;trait = 1)
     return [getindex.(wa,trait) for wa in s.agentarray],[fill(t,size(s[j])) for (j,t) in enumerate(s.tspan)]

src/ABMEv_runworld.jl

@@ -4,25 +4,35 @@
 Run `w` with algorithm `alg`, until `tend` is reached. User needs to provide `b` the birth function,
 which takes as arguments `X,t`, and provide `d` the death function, with arguments `X,Y,t`.
 Returns a `Simulation` type.
-- `worldall` stores the world every `p["dt_saving"]` time steps.
+- if `dt_saving` specified, world is saved every time steps.
 If `dt_saving` not specified, `sim` contains an array of two elements,
 first corresponding to initial conditions and last corresponding to world in the last time step.
+- if `t_saving_cb::Array{Float64}` specified, callbacks are computed at each steps time specified in the array.
+This functionality is as of now only compatible with `dt_saving` not specified.
 - `cb` correspond to callbacks function. Look at the documentation for more information
 - the run stops if the number of agents reaches`p["NMax"]`.
 """
 function run!(w::World{A,S,T},alg::L,tend::Number,b,d;
                 dt_saving=nothing,
+                t_saving_cb=nothing,
                 cb=(names = String[],agg =nothing)) where {A,S,T,L<:AbstractAlg}
+    # argument check
     _check_timedep(b,d)
+    (!isnothing(dt_saving) && !isnothing(dt_saving)) ? ArgumentError("For now, can not specify both `dt_saving` and `t_saving_cb`") : nothing
+    isnothing(dt_saving) ? dt_saving =  tend + 1. : nothing
+    isnothing(t_saving_cb) ? t_saving_cb =  [tend + 1.] : nothing
+
+    # var init
     n=size(w);
     NMax = maxsize(w)
     t = .0
     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!(w,alg,b,d)
     end
+
+    # start
     while t<tend
         if dt < 0
             throw("We obtained negative time step dt = $dt at event $i")
@@ -37,6 +47,11 @@ function run!(w::World{A,S,T},alg::L,tend::Number,b,d;
             @info "saving world @ t = $(t)/ $(tend)"
             add_entry!(sim,w)
         end
+        if  t >= first(t_saving_cb)
+            @info "saving callback only @ t = $(t)/ $(tend)"
+            add_entry_cb_only!(sim,w)
+            popfirst!(t_saving_cb)
+        end
         dt = updateWorld!(w,alg,b,d)
         t +=  dt
         i += 1

test/simulation.jl

@@ -33,14 +33,15 @@ 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
+b(X,t) = gaussian(X[1],0.,sigma_K)
+d(X,Y,t) = 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,mu,NMax
 
+## testing cb
 myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
 w0 = World(myagents,myspace,p,0.)
 w1 = copy(w0)
@@ -53,6 +54,20 @@ eltype(cb.agg)
 @test typeof(get_world(sim,get_size(sim))) <: World
 @test typeof(sim[2]) <: Vector
 
+##testing t_saving_cb
+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,b,d,cb=cb,t_saving_cb = collect(1.0:0.1:1.5))
+@test typeof(sim["gamma_div"]) <: Vector
+@test get_size(sim) == length(sim["gamma_div"])
+@test length(sim[1]) == 1000
+@test length(sim[end]) > 1
+@test length(sim[2]) == 1
+
+
 ## testing plot
 using Plots
 plot(get_tspan(sim),sim["gamma_div"])