Adding time dependency to birth and death coefficient

.gitlab-ci.yml

@@ -16,13 +16,10 @@
 
 # Name a test and select an appropriate image.
 # images comes from Docker hub
-test:0.7:
-  image: julia:0.7
+test:1.4:
+  image: julia:1.4
   <<: *test_definition
 
-test:1.0:
-  image: julia:1.0
-  <<: *test_definition
 
 ## Below commented because gitlab.ethz.ch has disabled Pages feature
 # pages:

Project.toml

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

src/ABMEv.jl

@@ -38,5 +38,5 @@ module ABMEv
         truncvar,get_xhist_mat
     export update_afterbirth_std!,update_afterdeath_std!
     export generalised_gaussian,gaussian,ma,geomsmooth,arithsmooth,eth_grad_std,
-        DiversityFunction,geomsmooth2D,arithsmooth2D,interpolate_df,groupby
+        DiversityFunction,geomsmooth2D,arithsmooth2D,interpolate_df,groupby,numargs
 end

src/ABMEv_Agent.jl

@@ -200,20 +200,3 @@ function increment_x!(a::AbstractAgent{A,R},s::AbstractSpacesTuple,p::Dict,t::T)
     a.x_history[1] = _get_xinc(a,s,p,t)
     return a
 end
-
-"""
-    function tin(t::Number,a::Number,b::Number)
-if t in [a,b) returns 1. else returns 0
-"""
-
-function tin(t::Number,a::Number,b::Number)
-    return t>=a && t<b ? 1. : 0.
-end
-
-function split_move(t)
-    return .0 + 1/100*(t-20.)*tin(t,20.,120.) + tin(t,120.,Inf64)
-end
-
-function split_merge_move(t)
-    return .0 + 1/30*(t-10.)*tin(t,10.,40.) + tin(t,40.,70.) + (1- 1/30*(t-70.))*tin(t,70.,100.)
-end

src/ABMEv_runworld.jl

-function  update_rates_graph!(world,C,p::Dict,t::Float64)
-    for e in edges(p["g"])
-        # agents on e
-        aidx_e = findall(a -> get_x(a,1)==e,world)
-        na = length(aidx_e)
-        for i in aidx_e
-            world[i].d = na^2
-            world[i].b = 1
-        end
-    end
-end
-
 """
 $(SIGNATURES)
 Run `w` with algorithm `alg`, until `tend` is reached.
@@ -23,6 +11,7 @@ first corresponding to initial conditions and last corresponding to world in the
 function run!(w::World{A,S,T},alg::L,tend::Number;
                 dt_saving=nothing,
                 cb=(names = String[],agg =nothing)) where {A,S,T,L<:AbstractAlg}
+    _check_timedep(w.p)
     n=size(w);
     NMax = maxsize(w)
     t = .0
@@ -55,3 +44,20 @@ function run!(w::World{A,S,T},alg::L,tend::Number;
     @info "simulation stopped at t=$(t), after $(i) generations"
     return sim
 end
+
+"""
+function _correct_timedep!(p::Dict)
+
+checks time dependency of birth and death functions,
+and overloads the function if not provided
+"""
+
+function _check_timedep(p::Dict)
+    @unpack d,b = p
+    if numargs(p["b"]) < 2
+        throw(ArgumentError("Birth function needs `X` and `t` arguments"))
+    end
+    if numargs(p["d"]) < 3
+        throw(ArgumentError("Death function needs `X`, `Y` and `t` arguments"))
+    end
+end

src/ABMEv_utils.jl

@@ -5,6 +5,16 @@ function generalised_gaussian(x::Number,mu::Number,sigma::Number,epsilon::Number
         return exp( -.5 * ((x-mu) / sigma)^epsilon)
 end
 
+"""
+    function tin(t::Number,a::Number,b::Number)
+if t in [a,b) returns 1. else returns 0
+"""
+
+function tin(t::Number,a::Number,b::Number)
+    return a <= t < b ? 1. : 0.
+end
+
+
 """
         gaussian(x::Number,mu::Number,sigma::Number) = generalised_gaussian(x,mu,sigma,2)
 """
@@ -105,3 +115,12 @@ groupby(f, list::Array{T}) where {T}= begin
   end
   groups
 end
+
+"""
+$(SIGNATURES)
+returns the number of arguments of function `f`
+"""
+
+function numargs(f)
+    return maximum([length(m.sig.parameters) - 1 for m in methods(f)] )
+end

src/algo/ABMEv_CFM.jl

@@ -38,8 +38,8 @@ function updateWorld!(w::World{A,S,T},c::CFM) where {A,S,T}
     x = get_x(w[i])
     W = rand()
     if dt > 0.
-        deathprob = (sum(d.(get_x.(alive),Ref(x))) - d(x,x)) / (Cbar*(n+1))
-        birthprob = b(x) / (Cbar*(n+1))
+        deathprob = (sum(d.(get_x.(alive),Ref(x)),w.t) - d(x,x,w.t)) / (Cbar*(n+1))
+        birthprob = b(x,w.t) / (Cbar*(n+1))
         if W <= deathprob
             updateDeathEvent!(w,c,i)
         elseif W <= deathprob + birthprob

src/algo/ABMEv_Gillepsie.jl

@@ -21,21 +21,21 @@ function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int)
     x_offspring = get_x(offspring)
     _agents = agents(w)
     for a in _agents
-        a.d += d(get_x(a),x_offspring)
+        a.d += d(get_x(a),x_offspring,w.t)
     end
     # Now updating new agent
-    offspring.d = sum(d.(get_x.(_agents),Ref(x_offspring))) #- d(x_offspring,x_offspring)
-    offspring.b = b(x_offspring)
+    offspring.d = sum(d.(get_x.(_agents),Ref(x_offspring),w.t)) #- d(x_offspring,x_offspring)
+    offspring.b = b(x_offspring,w.t)
     addAgent!(w,offspring)
 end
 
 function updateDeathEvent!(world::World,::Gillepsie,i_event::Int)
-    @unpack d,b = parameters(world)
+    @unpack d = parameters(world)
     x_death = get_x(world[i_event])
     # updating death rate only the two columns corresponding to agent idx
     removeAgent!(world,i_event)
     for a in agents(world)
-        a.d -= d(get_x(a),x_death)
+        a.d -= d(get_x(a),x_death,world.t)
     end
 end
 
@@ -55,7 +55,7 @@ function  update_rates!(w::World,::Gillepsie)
     # 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])
+            C = d(traits[i],traits[j],w.t)
             D[i] += C
             D[j] += C
         end
@@ -63,7 +63,7 @@ function  update_rates!(w::World,::Gillepsie)
     # Here we can do  it in parallel as well
     for (i,a) in enumerate(_agents)
         a.d = D[i]
-        a.b = b(traits[i])
+        a.b = b(traits[i],w.t)
     end
 end
 

test/cfm.jl

@@ -10,8 +10,8 @@ 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

test/gillepsie.jl

@@ -10,8 +10,8 @@ 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

test/metrics_hamming.jl

@@ -10,8 +10,8 @@ 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

test/runtests.jl

@@ -9,4 +9,5 @@ using ABMEv, Test, JLD2,Random
     include("simulation.jl")
     include("space_agent.jl")
     include("world.jl")
+    include("utils.jl")
 end

test/simulation.jl

@@ -4,8 +4,8 @@ using Revise,ABMEv
 using UnPack
 myspace = (GraphSpace(SimpleGraph(10,10)),RealSpace{1,Float64}())
 myagents = [Agent(myspace,ancestors=true,rates=true) for i in 1:10]
-d(X,Y) = gaussian(X[1],Y[1],1)
-b(X,Y) = gaussian(X[1],0,1)
+d(X,Y,t) = gaussian(X[1],Y[1],1)
+b(X,Y,t) = gaussian(X[1],0,1)
 D = (Int64(1),Float64(1.))
 mu = [1,1]
 NMax = 100

test/utils.jl

+# using Test
+# using Revise
+# using ABMEv
+# using UnPack
+#
+# b(X) = gaussian(X[1],0.,sigma_K)
+# d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+# p = Dict{String,Any}();@pack! p = d,b
+# X = [.5]; Y = [.6];t = 0.
+# _correct_timedep(p)
+# @testset "Utils" begin
+#     @test numargs(ABMEv.d) == 3
+#     @test numargs(ABMEv.b) == 2
+#     ABMEv.b(X,t) == b(X)
+# end