adding a function clean to remove unused rows in sim, and some minor changes

src/ABMEv.jl

@@ -11,7 +11,7 @@ module ABMEv
     include("ABMEv_world.jl")
     include("ABMEv_Sim.jl")
     include("ABMEv_metrics.jl")
-    # include("ABMEv_plot.jl")
+    include("ABMEv_plot.jl")
     include("ABMEv_utils.jl")
     include("algo/ABMEv_WF.jl")
     include("algo/ABMEv_Gillepsie.jl")

src/ABMEv_Agent.jl

@@ -141,7 +141,7 @@ $(SIGNATURES)
 Returns trait `i` of the agent.
 Geotrait corresponds to dimension `i=0`.
 """
-get_x(a::Agent,t::Number,i::Integer) = i > 0 ? a.x_history[end][Int(i)] : get_geo(a,t)
+get_x(a::AbstractAgent,t::Number,i::Integer) = i > 0 ? a.x_history[end][Int(i)] : get_geo(a,t)
 
 """
 $(SIGNATURES)

src/ABMEv_Sim.jl

@@ -57,6 +57,11 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
     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))
+# @deprecate get_x(agentarray::Array{T},t::Number,trait::Integer)
+
+
 #TODO: code it
 # function world2df(world::Array{T,1},geotrait=false) where {T <: Agent}
 #     xx = get_xarray(world)

src/ABMEv_plot.jl

@@ -15,29 +15,32 @@ It should correspond to an integer, as it indexes the column to plot
 - `"xs"`
 """
 
-@recipe function plot(world::Array{U},p;what=["x","H"],trait = 1,tplot = 0) where U <: Union{Missing,Agent}
+@recipe function plot(sim::Simulation;what=["x","H"],trait = 1,tplot = 0)
+    world = sim.agentarray
     tot_dim = size(world,2)*size(world,1)
+    tspan = sim.tspan
+    p = sim.p
     # We reduce time interval if it is too big
     # if tot_dim > 1e6 && size(world,2) >= 200
     #     p = copy(p)
     #     idx_reduced = floor.(Int,range(1,size(world,2),length = 200))
-    #     p["tspan" ] = p["tspan"][idx_reduced]
+    #     p["tspan" ] = tspan[idx_reduced]
     #     world = world[:,idx_reduced]
     # end
     # second condition is to make sure that the world corresponds to all the time steps
     # If not, then we can not get "x"
-    if count(ismissing,world) > 0 && length(p["tspan"]) == size(world,2)
-        tspan_ar = vcat([p["tspan"][i]*ones(Int(p["NMax"] - count(ismissing,world[:,i]))) for i in 1:length(p["tspan"]) ]...);
+    if count(ismissing,world) > 0 && length(tspan) == size(world,2)
+        tspan_ar = vcat([tspan[i]*ones(Int(p["NMax"] - count(ismissing,world[:,i]))) for i in 1:length(tspan) ]...);
     else
-        tspan_ar = repeat(p["tspan"],inner = size(world,1))
+        tspan_ar = repeat(tspan,inner = size(world,1))
     end
     # tspan = Float64.(tspan)
-    tend = p["tspan"][tplot > 0 ? tplot : length(p["tspan"])]
+    tend = tspan[tplot > 0 ? tplot : length(tspan)]
     world_sm = clean_world(world)
     if "x" in what
         d_i = []
-        for i in 1:size(world,2)
-            x = get_x(clean_world(world[:,i]),p["tspan"][i],trait)[:]
+        for i in 1:length(tspan)
+            x = get_x.(clean_world(world[:,i]),tspan[i],trait)[:]
             append!(d_i,pdf(kde(x),x))
         end
         @series begin
@@ -59,40 +62,39 @@ It should correspond to an integer, as it indexes the column to plot
     end
     # we use this for discrete agents
     # world should be a one dimensional vector, corresponding to one time step only
-    if "xs" in what
-        d_i = []; xt_array = []; x1_array = []
-        world_df_all = world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true)
-        world_df_g = groupby(world_df_all,:x1)
-        for world_df in world_df_g
-            if trait == 0
-                x = Float64.(world_df.g)
-            else
-                # fitness occupies first spot
-                x = world_df[:,trait+1] ;
-            end
-            x1 =  world_df.x1;
-            append!(d_i,pdf(kde(x),x))
-            append!(xt_array,x)
-            append!(x1_array,x1)
-        end
-        # TODO: we stopped here
-        @series begin
-            seriestype := :scatter
-            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
-            # markercolor := :blue
-            markerstrokewidth := 0
-            # seriesalpha := 1.
-            xaxis := "geographical position"
-            xticks :=  sort!(unique(world_df_all.x1))
-            yaxis := "trait value"
-            label := ""
-            grid := false
-            # marker := (:rect,20,1.)
-            x1_array[:],xt_array[:]
-        end
-    end
+    # if "xs" in what
+    #     d_i = []; xt_array = []; x1_array = []
+    #     world_df_all = world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true)
+    #     world_df_g = groupby(world_df_all,:x1)
+    #     for world_df in world_df_g
+    #         if trait == 0
+    #             x = Float64.(world_df.g)
+    #         else
+    #             # fitness occupies first spot
+    #             x = world_df[:,trait+1] ;
+    #         end
+    #         x1 =  world_df.x1;
+    #         append!(d_i,pdf(kde(x),x))
+    #         append!(xt_array,x)
+    #         append!(x1_array,x1)
+    #     end
+    #     @series begin
+    #         seriestype := :scatter
+    #         markercolor := eth_grad_small[d_i ./ maximum(d_i)]
+    #         # markercolor := :blue
+    #         markerstrokewidth := 0
+    #         # seriesalpha := 1.
+    #         xaxis := "geographical position"
+    #         xticks :=  sort!(unique(world_df_all.x1))
+    #         yaxis := "trait value"
+    #         label := ""
+    #         grid := false
+    #         # marker := (:rect,20,1.)
+    #         x1_array[:],xt_array[:]
+    #     end
+    # end
     if "gs" in what
-        _world = clean_world(world[:, tplot > 0 ? tplot : length(p["tspan"]) ])
+        _world = clean_world(world[:, tplot > 0 ? tplot : length(tspan) ])
         y = get_x(_world,tend,2)[:]
         x = get_x(_world,tend,0)[:]
         X = hcat(x,y)
@@ -121,8 +123,8 @@ It should correspond to an integer, as it indexes the column to plot
         d_i = []
         for i in 1:size(world,2)
             _world = clean_world(world[:,i])
-            x = get_x(_world,p["tspan"][i],2)[:]
-            y = get_x(_world,p["tspan"][i],0)[:]
+            x = get_x(_world,tspan[i],2)[:]
+            y = get_x(_world,tspan[i],0)[:]
             X = hcat(x,y)
             # d = kde(X)
             # di_temp = diag(pdf(d,X[:,1],X[:,2]))
@@ -149,8 +151,8 @@ It should correspond to an integer, as it indexes the column to plot
     if "3d" in what
         d_i = []
         for i in 1:size(world,2)
-            x = get_x(clean_world(world[:,i]),p["tspan"][i],1)[:]
-            y = get_x(clean_world(world[:,i]),p["tspan"][i],2)[:]
+            x = get_x(clean_world(world[:,i]),tspan[i],1)[:]
+            y = get_x(clean_world(world[:,i]),tspan[i],2)[:]
             X = hcat(x,y)
             d = kde(X)
             di_temp = diag(pdf(d,X[:,1],X[:,2]))
@@ -188,7 +190,7 @@ It should correspond to an integer, as it indexes the column to plot
             label := "Variance"
             xlabel := "Time"
             ylabel := "Variance"
-            p["tspan"],var(world_sm,trait=trait)[:]
+            tspan,var(world_sm,trait=trait)[:]
         end
     end
     if "vargeo" in what
@@ -198,7 +200,7 @@ It should correspond to an integer, as it indexes the column to plot
             label := "Variance of geotrait"
             xlabel := "Time"
             ylabel := "Variance"
-            p["tspan"],i->first(covgeo(world_sm[:,Int(i)]))
+            tspan,i->first(covgeo(world_sm[:,Int(i)]))
         end
     end
     # if "density_t" in what
@@ -208,7 +210,7 @@ It should correspond to an integer, as it indexes the column to plot
     #         label := "Variance of geotrait"
     #         xlabel := "Time"
     #         ylabel := "Variance"
-    #         p["tspan"],i->first(covgeo(world_sm[:,Int(i)]))
+    #         tspan,i->first(covgeo(world_sm[:,Int(i)]))
     #     end
     # end
 end

src/ABMEv_runworld.jl

@@ -82,6 +82,7 @@ function run!(w::World{A,S,T},alg::L,tend::Number;
     end
     # Saving last time step
     add_entry!(sim,w)
+    clean!(sim)
     @info "simulation stopped at t=$(t), after $(i) generations"
     return sim
 end

src/algo/ABMEv_CFM.jl

@@ -15,8 +15,8 @@ function updateBirthEvent!(w::World,::CFM,mum_idx::Int)
     addAgent!(w,offspring)
 end
 
-function updateDeathEvent!(world::World,::CFM,i_event::Int)
-    removeAgent!(world,i_event)
+function updateDeathEvent!(world::World,::CFM,i::Int)
+    removeAgent!(world,i)
 end
 
 """
@@ -27,34 +27,26 @@ DOI : 10.1016/j.tpb.2005.10.004
 """
 function updateWorld!(w::World{A,S,T},c::CFM) where {A,S,T}
     # total update world
+    @unpack Cbar,d,b = parameters(w)
     alive = agents(w)
-    events_weights = ProbabilityWeights(vcat(get_d.(alive),get_b.(alive)))
     # Total rate of events
-    ∑ = sum(events_weights)
-    dt = - log(rand(T))/∑
+    n = size(w)
+    ∑ = (n + 1)
+    dt = - log(rand(T))/(Cbar * ∑)
     update_clock!(w,dt)
+    i = rand(1:n)
+    x = get_x(w[i])
+    W = rand()
     if dt > 0.
-        i_event = sample(events_weights)
-        # This is ok since size is a multiple of 2
-        n = size(w)
-        if i_event <= n
-            # DEATH EVENT
-            # In this case i_event is also the index of the individual to die in the world_alive
-            updateDeathEvent!(w,g,i_event)
-        else
-            # birth event
-            # i_event - n is also the index of the individual to give birth in the world_alive
-            updateBirthEvent!(w,g,i_event-n)
+        deathprob = (sum(d.(get_x.(alive),Ref(x))) - d(x,x)) / (Cbar*(n+1))
+        birthprob = b(x) / (Cbar*(n+1))
+        if W <= deathprob
+            updateDeathEvent!(w,c,i)
+        elseif W <= deathprob + birthprob
+            updateBirthEvent!(w,c,i)
         end
         return dt
     else
         return -1
     end
 end
-
-"""
-    clean_world(world::Array{T}) where {T <: Agent}
-Get rid of missing value in `world`
-"""
-#TODO : put some type specs here
-clean_world(world) = collect(skipmissing(world))

test/cfm.jl

+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 = 10000
+tend = 100
+Cbar = 2
+p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar
+
+myagents = [Agent(myspace,(0,),ancestors=false,rates=false) for i in 1:K0]
+w0 = World(myagents,myspace,p,0.)
+w1 = copy(w0)
+@info "Running simulation with CFM algorithm"
+@time sim = run!(w1,CFM(),tend,dt_saving=1.)
+
+ABMEv.clean!(sim)
+using Plots
+Plots.plot(sim)

test/profiling_ABMEv.jl

@@ -79,3 +79,36 @@ tspan=zeros(1)
 # using BenchmarkTools
 # update_afterbirth_std!(skipmissing(world0),1,p)
 # @btime update_afterdeath_std!(skipmissing(world0),[.8],p_default)
+
+################################
+#####Testing Ancestors##########
+################################
+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 = 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.)
+@info "Running simulation with Gillepsie algorithm"
+@time sim = run!(w0,Gillepsie(),tend)
+
+myagents = [Agent(myspace,(0,),ancestors=false,rates=true) for i in 1:K0]
+w0 = World(myagents,myspace,p,0.)
+@info "Running simulation with Gillepsie algorithm"
+@time sim = run!(w0,Gillepsie(),tend)