Commit 8e33b260 authored by Victor's avatar Victor

big changes on passing functions

parent e0f9c937
Pipeline #77939 failed with stage
in 19 minutes and 5 seconds
......@@ -8,17 +8,17 @@ for K0 in [10,50,100,500,1000,5000]
myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
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
tend = 2
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
myagents = [Agent(myspace,(0,),ancestors=false,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
@show K0
@time sim = run!(w0,Gillepsie(),tend)
@time sim = run!(w0,Gillepsie(),tend,b,d)
end
println("--------------------------------")
......@@ -28,18 +28,18 @@ for K0 in [10,50,100,500,1000,5000]
myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
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
tend = 2
Cbar=2
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar
myagents = [Agent(myspace,(0,),ancestors=true) for i in 1:K0]
dm = d([0],[0],0.);bm = 1.
p = Dict{String,Any}();@pack! p = dm,bm,D,mu,NMax,Cbar
myagents = [Agent(myspace,(0,),ancestors=false) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
@show K0
@time sim = run!(w0,CFM(),tend)
@time sim = run!(w0,CFM(),tend,b,d)
end
println("--------------------------------")
......
......@@ -10,10 +10,10 @@ D = (1e-2,)
mu = [.1]
NMax = 2000
tend = 1500
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
p = Dict{String,Any}();@pack! p = D,mu,NMax
myagents = [Agent(myspace,(1e-2 * randn(),),rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
@time sim = run!(w0,Gillepsie(),tend,dt_saving = 10)
@time sim = run!(w0,Gillepsie(),tend,dt_saving = 10,b,d)
using JLD2
@save joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim
......
......@@ -7,15 +7,15 @@ myspace = (RealSpace{1,Float64}(),)
K0 = 1000
b(X,t) = 1.
d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b)
Cbar = b([0],0.) + d([0],[0],0.)
D = (1e-2,)
mu = [.1]
NMax = 2000
tend = 1500
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar
dm = d([0],[0],0.);bm = 1.
p = Dict{String,Any}();@pack! p = dm,bm,D,mu,NMax
myagents = [Agent(myspace,(1e-2 * randn(),)) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
@time sim = run!(w0,CFM(),tend,dt_saving = 4)
@time sim = run!(w0,CFM(),tend,dt_saving = 20,b,d)
using JLD2
@save joinpath(@__DIR__,"sim_sympatric_speciation_CFM.jld2") sim
......
......@@ -96,9 +96,11 @@ end
import Base:copy,show
Base.copy(a::A) where {A<:AbstractAgent} = A(copy(a.x_history),copy(a.t_history),copy(a.d),copy(a.b))
# this function only copies the trait history and time (x,t), and set birth and death rates to 0.
copyxt(a::Agent{A,R,T,U,V}) where {A,R,T,U,V<:Number} = Agent{A,R,T,U,V}(copy(a.x_history),copy(a.t_history),zero(V),zero(V))
copyxt(a::Agent{A,R,T,U,Nothing}) where {A,R,T,U} = Agent{A,R,T,U,Nothing}(copy(a.x_history),copy(a.t_history),nothing,nothing)
# this has to be overloaded for Base.copy(a::Agent) to work properly
Base.copy(m::Missing) = missing
Base.copy(n::Nothing) = nothing
function Base.show(io::IO, a::Agent{A,R,T,U,V}) where {A,R,T,U,V}
......
"""
function run!(w::World{A,S,T},alg::L,tend::Number;dt_saving=nothing,cb=(names = String[],agg =nothing))
function run!(w::World{A,S,T},alg::L,tend::Number,b,d;dt_saving=nothing,cb=(names = String[],agg =nothing))
Run `w` with algorithm `alg`, until `tend` is reached.
Run `w` with algorithm `alg`, until `tend` is reached. User needs to provide `b` the birth function,
which takes as arguments `X,t`, and provide `d` the death function, with arguments `X,Y,t`.
Returns a `Simulation` type.
- `worldall` stores the world every `p["dt_saving"]` time steps.
If `dt_saving` not specified, `sim` contains an array of two elements,
......@@ -9,10 +10,10 @@ first corresponding to initial conditions and last corresponding to world in the
- `cb` correspond to callbacks function. Look at the documentation for more information
- the run stops if the number of agents reaches`p["NMax"]`.
"""
function run!(w::World{A,S,T},alg::L,tend::Number;
function run!(w::World{A,S,T},alg::L,tend::Number,b,d;
dt_saving=nothing,
cb=(names = String[],agg =nothing)) where {A,S,T,L<:AbstractAlg}
_check_timedep(w.p)
_check_timedep(b,d)
n=size(w);
NMax = maxsize(w)
t = .0
......@@ -20,7 +21,7 @@ function run!(w::World{A,S,T},alg::L,tend::Number;
isnothing(dt_saving) ? dt_saving = tend + 1. : nothing
sim = Simulation(w,cb=cb)
if A <: AbstractAgent{AA,Rates{true}} where {AA}
update_rates!(w,alg)
update_rates!(w,alg,b,d)
end
while t<tend
if dt < 0
......@@ -36,7 +37,7 @@ function run!(w::World{A,S,T},alg::L,tend::Number;
@info "saving world @ t = $(t)/ $(tend)"
add_entry!(sim,w)
end
dt = updateWorld!(w,alg)
dt = updateWorld!(w,alg,b,d)
t += dt
i += 1
end
......@@ -52,12 +53,11 @@ end
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
function _check_timedep(b,d)
if numargs(b) < 2
throw(ArgumentError("Birth function needs `X` and `t` arguments"))
end
if numargs(p["d"]) < 3
if numargs(d) < 3
throw(ArgumentError("Death function needs `X`, `Y` and `t` arguments"))
end
end
......@@ -25,9 +25,9 @@ Updating rule for CFM setting.
algorithm described in article
DOI : 10.1016/j.tpb.2005.10.004
"""
function updateWorld!(w::World{A,S,T},c::CFM) where {A,S,T}
function updateWorld!(w::World{A,S,T},c::CFM,b,d) where {A,S,T}
# total update world
@unpack d,b,bm,dm = parameters(w)
@unpack bm,dm = parameters(w)
alive = agents(w)
# Total rate of events
n = size(w)
......
# In this file lie the function for Gillepsie algorithm
"""
$(TYPEDEF)
"""
......@@ -9,28 +10,28 @@ export Gillepsie
Used for Gillepsie setting
"""
function give_birth(mum_idx::Int,w::World)
new_a = copy(w[mum_idx])
new_a = copyxt(w[mum_idx])
increment_x!(new_a,space(w),parameters(w),time(w))
return new_a
end
function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int)
function updateBirthEvent!(w::World,::Gillepsie,mum_idx::Int,b,d)
# updating competition only the two columns corresponding to agent idx
@unpack d,b = parameters(w)
offspring = give_birth(mum_idx,w)
x_offspring = get_x(offspring)
_agents = agents(w)
for a in _agents
a.d += d(get_x(a),x_offspring,w.t)
# Warning : symmetric competition
tmp = d(get_x(a),x_offspring,w.t)
a.d += tmp
offspring.d += tmp
end
# Now updating new agent
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 = parameters(world)
function updateDeathEvent!(world::World,::Gillepsie,i_event::Int,d)
x_death = get_x(world[i_event])
# updating death rate only the two columns corresponding to agent idx
removeAgent!(world,i_event)
......@@ -40,13 +41,12 @@ function updateDeathEvent!(world::World,::Gillepsie,i_event::Int)
end
"""
update_rates_std!(world,p::Dict,t::Float64)
update_rates!(w::World,::Gillepsie,b,d)
This standard updates takes
- competition kernels of the form α(x,y) and
- carrying capacity of the form K(x)
"""
function update_rates!(w::World,::Gillepsie)
@unpack b,d = parameters(w)
function update_rates!(w::World,::Gillepsie,b,d)
_agents = agents(w)
traits = get_x.(_agents)
# traits = get_xhist.(world)
......@@ -68,13 +68,11 @@ function update_rates!(w::World,::Gillepsie)
end
"""
function updateWorld_G!(world,t,p)
function updateWorld!(w::World{A,S,T},g::G,b,d)
Updating rule for gillepsie setting.
Returning dt drawn from an exponential distribution with parameter the total rates of events.
# Args
t is for now not used but might be used for dynamic landscape
Returning `dt` drawn from an exponential distribution with parameter the total rates of events.
"""
function updateWorld!(w::World{A,S,T},g::G) where {A,S,T,G <: Gillepsie}
function updateWorld!(w::World{A,S,T},g::G,b,d) where {A,S,T,G <: Gillepsie}
# total update world
alive = agents(w)
events_weights = ProbabilityWeights(vcat(get_d.(alive),get_b.(alive)))
......@@ -89,11 +87,11 @@ function updateWorld!(w::World{A,S,T},g::G) where {A,S,T,G <: Gillepsie}
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)
updateDeathEvent!(w,g,i_event,d)
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)
updateBirthEvent!(w,g,i_event-n,b,d)
end
return dt
else
......
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