Commit 9a5e454e authored by Victor's avatar Victor
Browse files

changes on CFM algo

parent 0e0d0521
Pipeline #76424 passed with stage
in 20 minutes and 17 seconds
......@@ -27,19 +27,19 @@ 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)
@unpack d,b,bm,dm = parameters(w)
alive = agents(w)
# Total rate of events
n = size(w)
= (n + 1)
dt = - log(rand(T))/(Cbar * )
Cbar = bm + dm*n
dt = rand(Exponential(Cbar)) / n
update_clock!(w,dt)
i = rand(1:n)
x = get_x(w[i])
W = rand()
if dt > 0.
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))
deathprob = (sum(d.(get_x.(alive),Ref(x),w.t)) .- d(x,x,w.t)) / Cbar
birthprob = b(x,w.t) / Cbar
if W <= deathprob
updateDeathEvent!(w,c,i)
elseif W <= deathprob + birthprob
......
......@@ -10,21 +10,22 @@ myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
b(X,t) = gaussian(X[1],0.,sigma_K)
d(X,Y,t) = gaussian(X[1],Y[1],sigma_a)/K0
b(X,t) = 1.
d(X,Y,t) = gaussian(X[1],Y[1],sigma_a)/gaussian(X[1],0.,sigma_K)/K0
D = (1e-2,)
mu = [.1]
NMax = 10000
tend = 100
Cbar = 2
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar
NMax = 2000
tend = 2000
# Cbar = b([0],0.)/K0 + d([0],[0],0.)
dm = d([0],[0],0.);bm = 1.
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,dm,bm
myagents = [Agent(myspace,(0,),ancestors=false,rates=false) for i in 1:K0]
myagents = [Agent(myspace,(-.01 + 1e-2 * randn(),),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.)
@time sim = run!(w1,CFM(),tend,dt_saving=10.)
ABMEv.clean!(sim)
# ABMEv.clean!(sim)
using Plots
Plots.plot(sim)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment