Version 0.2. Changed the structure of agents, to better integrate the...

Version 0.2. Changed the structure of agents, to better integrate the geotrait. Changes have telescoped to accessor get_geo that now takes t as an argument, and all the functions that are using it.

Project.toml

 name = "ABMEv"
 uuid = "837ac870-fb52-4b0c-9a0e-030f2f36f5ed"
 authors = ["Victor Boussange "]
-version = "0.1.8"
+version = "0.2.1"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

examples/gillepsie_discrete_2d.jl

-
 using ABMEv,UnPack
 cd(@__DIR__)
 using Dates
@@ -16,7 +15,7 @@ K0 = 1000
 K(X) = gaussian(X[2],X[1] * a,sigma_K) / nodes
 alpha(X,Y) = (X[1] ≈ Y[1]) * gaussian(X[2],Y[2],sigma_a) / K0
 NMax = 5000
-tend = 3000.
+tend = 1.
 reflected = true;
 dt_saving = 15.
 

src/ABMEv_Agent.jl

@@ -4,6 +4,8 @@ abstract type MixedAgent end
 mutable struct Agent{T,U}
     # history of traits for geotraits
     x_history::Array{U}
+    # birth time of ancestors
+    t_history::Array{U,1}
     # death rate
     d::Float64
     #birth rate
@@ -12,14 +14,14 @@ end
 
 # Constructors
 # This  constructor should be used when one wants to impose the type of the agent (e.g. Mixed)
-Agent{T}(xhist::Array{U}) where {T,U} = Agent{T,U}(reshape(xhist,:,1),0.,1.)
+Agent{T}(xhist::Array{U}) where {T,U} = Agent{T,U}(reshape(xhist,:,1),[0.],0.,1.)
 
 # This constructor is used by default
 Agent(xhist::Array{U}) where {U <: Number} = Agent{StdAgent}(xhist)
 
-Agent() = Agent(Float64[],0.,1.)
+Agent() = Agent(Float64[],0.,0.,1.)
 import Base.copy
-copy(a::Agent{T,U}) where {T,U} = Agent{T,U}(a.x_history,a.d,a.b)
+copy(a::Agent{T,U}) where {T,U} = Agent{T,U}(copy(a.x_history),copy(a.t_history),copy(a.d),copy(a.b))
 copy(m::Missing) = missing
 
 """
@@ -39,8 +41,14 @@ end
 get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
 get_xhist(a::Agent) = a.x_history
 get_x(a::Agent) = a.x_history[:,end]
-get_geo(a::Agent) = sum(get_xhist(a,1))
-get_x(a::Agent,i::Number) = i > 0 ? a.x_history[Int(i),end] : get_geo(a)
+function get_geo(a::Agent{U,T},t::Number) where {U,T}
+    tarray = vcat(a.t_history[2:end],convert(T,t))
+    tarray .-= a.t_history
+    return sum(get_xhist(a,1) .* tarray)
+end
+# This method can acces geotrait, while the second not
+get_x(a::Agent,t::Number,i::Integer) = i > 0 ? a.x_history[Int(i),end] : get_geo(a,t)
+get_x(a::Agent,i::Integer) = a.x_history[Int(i),end]
 get_d(a::Agent) = a.d
 get_b(a::Agent) = a.b
 get_fitness(a::Agent) = a.b - a.d
@@ -48,23 +56,29 @@ get_dim(a::Agent) = size(a.x_history,1)
 get_nancestors(a::Agent) = size(a.x_history,2)
 """
     get_xarray(world::Array{Agent},trait::Int)
-If trait = 0 , we return the geotrait.
 Mainly works for WF-type world
 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.
 Particularly suited for an array world corresponding to a timeseries.
 
 """
-get_xarray(world::Array{T},trait::Int) where {T <: Agent} = trait > 0 ? reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2)) : reshape(hcat(get_geo.(world)),size(world,1),size(world,2))
+get_x(world::Array{T},trait::Integer) where {T <: Agent} = trait > 0 ? reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2)) : throw(ErrorException("Not the right method, need `t` as an argument"))
+
+get_x(world::Array{T},t::Number,trait::Integer) where {T <: Agent} = trait > 0 ? reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2)) : reshape(hcat(get_geo.(world,t)),size(world,1),size(world,2))
 
 """
     get_xarray(world::Array{Agent,1})
 Returns every traits of every agents of world in the form of an array
 If geotrait = true, then a last trait dimension is added, corresponding to geotrait.
 """
-function get_xarray(world::Array{T,1},geotrait::Bool=false) where {T <: Agent}
+function get_xarray(world::Array{T,1}) where {T <: Agent}
+    return hcat(get_x.(world)...)
+end
+
+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)')
+        xarray = vcat( xarray, get_geo.(world,t)')
     end
     return xarray
 end
@@ -80,33 +94,57 @@ If geotrait = true, then a last trait dimension is added, corresponding to geotr
 Note that because number of ancestors are different between agents, we return an array which size corresponds to the minimum of agents ancestors,
 and return the last generations, dropping the youngest ones
 """
-function get_xhist(world::Vector{T},geotrait = false) where {T <: Agent}
+function get_xhist(world::Vector{T}) where {T <: Agent}
     hist = minimum(get_nancestors.(world))
     ntraits = get_dim(first(world));
     xhist = zeros(length(world), hist, ntraits + geotrait);
     for (i,a) in enumerate(world)
         xhist[i,:,1:end-geotrait] = get_xhist(a)[:,end-hist+1:end]';
-        if geotrait
-            xhist[i,:,ntraits+geotrait] = cumsum(get_xhist(a,1))[end-hist+1:end]
-        end
     end
     return xhist
 end
 
+# TODO: This method broken, when one ask for the geotraits
+# function get_xhist(world::Vector{T},t::Number,geotrait = false) where {T <: Agent}
+#     hist = minimum(get_nancestors.(world))
+#     ntraits = get_dim(first(world));
+#     xhist = zeros(length(world), hist, ntraits + geotrait);
+#     for (i,a) in enumerate(world)
+#         xhist[i,:,1:end-geotrait] = get_xhist(a)[:,end-hist+1:end]';
+#         if geotrait
+#             xhist[i,:,ntraits+geotrait] = cumsum(get_xhist(a,1))[end-hist+1:end]
+#         end
+#     end
+#     return xhist
+# end
+
+
+function world2df(world::Array{T,1}) where {T <: Agent}
+    xx = get_xarray(world)
+    dfw = DataFrame(:f => get_fitness.(world))
+    for i in 1:size(xx,1)
+        dfw[Meta.parse("x$i")] = xx[i,:]
+    end
+    if geotrait
+        dfw[:g] = get_geo.(world)
+    end
+    return dfw
+end
+
 """
-    world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}
+    world2df(world::Array{T,1},t::Number,geotrait = false) where {T <: Agent}
 Converts the array of agent world to a datafram, where each column corresponds to a trait of the
 agent, and an extra column captures fitness.
 Each row corresponds to an agent
 """
-function world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}
+function world2df(world::Array{T,1},t::Number,geotrait = false) where {T <: Agent}
     xx = get_xarray(world)
     dfw = DataFrame(:f => get_fitness.(world))
     for i in 1:size(xx,1)
         dfw[Meta.parse("x$i")] = xx[i,:]
     end
     if geotrait
-        dfw[:g] = get_geo.(world)
+        dfw[:g] = get_geo.(world,t)
     end
     return dfw
 end
@@ -114,11 +152,11 @@ end
 
 
 """
-increment_x!(a::Agent{StdAgent,U},p::Dict)
+increment_x!(a::Agent{StdAgent,U},t::U,p::Dict) where U
     This function increments agent by random numbers specified in p
     ONLY FOR CONTINUOUS DOMAINS
 """
-function increment_x!(a::Agent{StdAgent,U},p::Dict) where U
+function increment_x!(a::Agent{StdAgent,U},t,p::Dict) where U
     tdim = length(p["D"])
     reflected = haskey(p,"reflected") ? p["reflected"] : false
     if reflected
@@ -131,21 +169,23 @@ function increment_x!(a::Agent{StdAgent,U},p::Dict) where U
         inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
     end
     a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
+    push!(a.t_history,t)
  end
 
  """
-     function increment_x!(a::Agent{MixedAgent,U},p::Dict)
+     function increment_x!(a::Agent{MixedAgent,U},t::U,p::Dict) where U
  This function increments first trait of agent with integer values, that are automatically reflected between 1 and p["nodes"].
 Other traits are incremented as usual.
 TODO : make it work for a graph type landscape, where domain is not a line anymore.
  """
- function increment_x!(a::Agent{MixedAgent,U},p::Dict) where U
+ function increment_x!(a::Agent{MixedAgent,U},t,p::Dict) where U
      tdim = length(p["D"])
      inc = [round(get_inc_reflected(get_x(a,1),p["D"][1] *randn(),1,p["nodes"]))]
      if  tdim > 1
          inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
      end
      a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
+     push!(a.t_history,t)
 end
 
 

src/ABMEv_Gillepsie.jl

@@ -5,9 +5,9 @@
     give_birth(a::Agent,p)
 Used for Gillepsie setting
 """
-function give_birth(a::Agent,p::Dict)
+function give_birth(a::Agent,t,p::Dict)
     new_a = copy(a)
-    increment_x!(new_a,p)
+    increment_x!(new_a,t,p)
     return new_a
 end
 
@@ -62,7 +62,7 @@ function updateWorld_G!(world,p,update_rates!,t)
             idx_offspring = findfirst(ismissing,world)
             # i_event - N is also the index of the individual to give birth in the world_alive
             mum = world[idx_world[i_event-N]]
-            world[idx_offspring] = give_birth(mum,p)
+            world[idx_offspring] = give_birth(mum,t+dt,p)
             update_afterbirth_std!(world,idx_offspring,p)
         end
         return dt

src/ABMEv_metrics.jl

@@ -53,7 +53,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
 """
 function var(world::Array{U,1};trait=1) where U <: Union{Missing,Agent{T}} where T
     world = collect(skipmissing(world))
-    xarray = get_xarray(world,trait)
+    xarray = get_x(world,trait)
     return var(xarray,dims=1)
 end
 """
@@ -64,9 +64,9 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
 # Arguments
 
 """
-function covgeo(world::Array{U,1},trait = 0) where U <: Union{Missing,Agent}
+function covgeo(world::Array{U,1},t::Number,trait = 0) where U <: Union{Missing,Agent}
     world = collect(skipmissing(world))
-    xarray = get_geo.(world)
+    xarray = get_geo.(world,t)
     if trait > 0
         xstd = reshape(get_x.(world,trait),size(world,1),size(world,2))
         xarray = hcat(xarray,xstd)

src/ABMEv_plot.jl

@@ -30,15 +30,16 @@ It should correspond to an integer, as it indexes the column to plot
         tspan_ar = repeat(p["tspan"],inner = size(world,1))
     end
     # tspan = Float64.(tspan)
-    world_sm = collect(skipmissing(world))
+    tend = p["tspan"][tplot > 0 ? tplot : size(world,2)]
+    world_sm = clean_world(world)
     if "x" in what
         d_i = []
         for i in 1:size(world,2)
-            x = get_x.(skipmissing(world[:,i]),trait)
+            x = get_x.(clean_world(world[:,i]),tend,trait)
             append!(d_i,pdf(kde(x),x))
         end
         @series begin
-        xarray = get_xarray(world_sm,trait)
+        xarray = get_xarray(world_sm,tend,trait)
             seriestype := :scatter
             markercolor := eth_grad_small[d_i ./ maximum(d_i)]
             # markercolor := :blue
@@ -56,7 +57,7 @@ It should correspond to an integer, as it indexes the column to plot
     # 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_g = groupby(world2df(collect(skipmissing(world[:, tplot > 0 ? tplot : size(world,2) ])),geotrait=true),:x1)
+        world_df_g = groupby(world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true),:x1)
         for world_df in world_df_g
             if trait == 0
                 x = world_df.g
@@ -86,8 +87,8 @@ It should correspond to an integer, as it indexes the column to plot
     end
     if "3dgeo" in what
         @series begin
-        xarray = get_geo.(world_sm)
-        yarray = get_xarray(world_sm,2)
+        xarray = get_geo.(world_sm,tspan_ar)
+        yarray = get_x(world_sm,2)
             seriestype := :scatter3d
             markercolor := "blue"
             markerstrokewidth := 0
@@ -103,7 +104,7 @@ It should correspond to an integer, as it indexes the column to plot
     if "3d" in what
         @series begin
         xarray = get_xarray(world_sm,1)
-        yarray = get_xarray(world_sm,2)
+        yarray = get_x(world_sm,2)
             seriestype := :scatter3d
             markercolor := "blue"
             markerstrokewidth := 0

test/gillepsie.jl

@@ -20,7 +20,7 @@ p_default = Dict(
 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 = collect(skipmissing(worldall[:,end]))
+world_alive_test = clean_world(worldall[:,end])
 # @save "gillepsie_test.jld2" world_alive
 @load "gillepsie_test.jld2" world_alive
 ## Testing
@@ -31,7 +31,7 @@ world_alive_test = collect(skipmissing(worldall[:,end]))
         end
         ## Comparing simulation
         @testset "Matching new vs old results " begin
-                xarray = get_xarray(world_alive,1);xarray_test = get_xarray(world_alive_test,1);
+                xarray = get_x(world_alive,1);xarray_test = get_x(world_alive_test,1);
                 @test xarray ≈ xarray_test
         end
 

test/parallel.jl

@@ -31,7 +31,7 @@ agents = [Agent( [1e-2]  .* randn(1) .- .5) for i in 1:K0]
 # using JLD2
 # @save "wrightfisher_test.jld2" worldall p
 # @load "wrightfisher_test.jld2" worldall
-# xarray = get_xarray(worldall,1); xarray_test = get_xarray(worldall_test,1)
+# xarray = get_x(worldall,1); xarray_test = get_x(worldall_test,1)
 
 
 # @test xarray ≈ xarray_test

test/wrightfisher.jl

@@ -23,7 +23,7 @@ worldall_test,p["tspan"] = runWorld_store_WF(p,agents);
 ## load to  compare simulation
 # @save "wrightfisher_test.jld2" worldall p
 @load "wrightfisher_test.jld2" worldall
-xarray = get_xarray(worldall,1); xarray_test = get_xarray(worldall_test,1)
+xarray = get_x(worldall,1); xarray_test = get_x(worldall_test,1)
 @testset "Wright Fisher Algorithm" begin
         @test xarray ≈ xarray_test
 end