To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit e391c168 authored by Victor's avatar Victor
Browse files

updated CFM algorithm

parent e0cb73c8
Pipeline #75130 passed with stage
in 25 minutes and 11 seconds
...@@ -15,13 +15,33 @@ D = (5e-1,1.4825) ...@@ -15,13 +15,33 @@ D = (5e-1,1.4825)
NMax = 2000 NMax = 2000
# tend = 1.5 # tend = 1.5
tend = 3000 tend = 100
p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu,D
myagents = [Agent(myspace,(rand(Int8(1):Int8(nodes)),initnode),ancestors=true,rates=true) for i in 1:round(K0/nodes)] myagents = [Agent(myspace,(rand(Int8(1):Int8(nodes)),initnode),ancestors=true,rates=true) for i in 1:round(K0/nodes)]
w0 = World(myagents,myspace,p_default,0.) w0 = World(myagents,myspace,p_default,0.)
@time sim = run!(w0,Gillepsie(),tend,dt_saving=3.) @time sim = run!(w0,Gillepsie(),tend,dt_saving=3.)
using GraphPlot using GraphPlot,StatsBase
# Plotting genetic space
# This is to plot with consistency # This is to plot with consistency
locs_x, locs_y = spring_layout(g) locs_x, locs_y = spring_layout(g)
xarray,tspan = get_xnt(sim, trait=2)
# here we compute the number of individuals per genome, throught all the time steps
d_i = []
for x in xarray
_d_i = zeros(nv(g))
c = countmap(x)
_d_i[collect(keys(c))] .= values(c)
push!(d_i,_d_i)
end
# Here we normalize with respect to the max size of the whole world through time
d_i = [ (_d_i .- minimum(_d_i)) ./ (maximum(_d_i) .- minimum(_d_i)) for _d_i in d_i ]anim = @animate for t in 1:size(worldall,2)
Plots.plot(worldall,p_default,what=["gs"],trait=2,tplot=t,
title = @sprintf(" t = %1.2f",p_default["tspan"][t]),
ylims = (1,p_default["nodes"]),
xlims = (0,maximum(xgeo))
)
end
gif(anim,"gif_d2=$(p_default["D"][2])_d1=$(p_default["D"][1])_tend=$(p_default["tend"])_mu=1e0_1e0.gif",fps = 13)
...@@ -3,17 +3,21 @@ using ABMEv,UnPack,Plots ...@@ -3,17 +3,21 @@ using ABMEv,UnPack,Plots
myspace = (RealSpace{1,Float64}(),) myspace = (RealSpace{1,Float64}(),)
σ_b = .9; σ_b = .9;
σ_d = .7; σ_d = .7;
K0 = 1000
b(X,t) = 1. b(X,t) = 1.
d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b) d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b)
D = (1e-2,) D = (1e-2,)
mu = [.1] mu = [.1]
NMax = 2000 NMax = 2000
tend = 1000 tend = 1500
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
myagents = [Agent(myspace,(1e-2 * randn(),),ancestors=true,rates=true) for i in 1:K0] myagents = [Agent(myspace,(1e-2 * randn(),),ancestors=true,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.) w0 = World(myagents,myspace,p,0.)
@time sim = run!(w0,Gillepsie(),tend,dt_saving = 4) @time sim = run!(w0,Gillepsie(),tend,dt_saving = 4)
using JLD2
@save joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim
Plots.plot(sim, ylabel = "Adaptive trait") Plots.plot(sim, ylabel = "Adaptive trait")
savefig(joinpath(@__DIR__, "sympatric_speciation.png")) savefig(joinpath(@__DIR__, "sympatric_speciation.png"))
...@@ -28,8 +32,5 @@ xplot = Plots.plot(thist,xhistall, ...@@ -28,8 +32,5 @@ xplot = Plots.plot(thist,xhistall,
grid = false, grid = false,
xlabel = "time", xlabel = "time",
ylabel = "Historical adaptive trait" ylabel = "Historical adaptive trait"
) )a
savefig(joinpath(@__DIR__, "x_hist_sympatric_speciation.png")) savefig(joinpath(@__DIR__, "x_hist_sympatric_speciation.png"))
using JLD2
@save joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim
...@@ -38,7 +38,7 @@ function updateWorld!(w::World{A,S,T},c::CFM) where {A,S,T} ...@@ -38,7 +38,7 @@ function updateWorld!(w::World{A,S,T},c::CFM) where {A,S,T}
x = get_x(w[i]) x = get_x(w[i])
W = rand() W = rand()
if dt > 0. if dt > 0.
deathprob = (sum(d.(get_x.(alive),Ref(x)),w.t) - d(x,x,w.t)) / (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)) birthprob = b(x,w.t) / (Cbar*(n+1))
if W <= deathprob if W <= deathprob
updateDeathEvent!(w,c,i) updateDeathEvent!(w,c,i)
......
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