Commit 07e3153e authored by Victor's avatar Victor
Browse files

version 0.1.2

reflected is now passed as an argument in p
an other argument "dt_saving" is introduced in p for Gillepsie, which allows one to save results every dt_saving step instead of fixed number. if not specified, only last time step is saved.
parent 41745feb
name = "ABMEv"
uuid = "837ac870-fb52-4b0c-9a0e-030f2f36f5ed"
authors = ["Victor Boussange "]
version = "0.1.0"
version = "0.1.2"
[deps]
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
......
......@@ -45,12 +45,13 @@ Returns trait of every agents of world in the form of an array
get_xarray(world::Array{Agent{T}},trait::Int) where T = reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2))
"""
increment_x!(a::Agent{Float64},p::Dict;reflected=false)
increment_x!(a::Agent{Float64},p::Dict)
This function increments current position by inc and updates xhist,
ONLY FOR CONTINUOUS DOMAINS
"""
function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
function increment_x!(a::Agent{Float64},p::Dict)
tdim = length(p["D"])
reflected = haskey(p,"reflected") ? p["reflected"] : false
if reflected
inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
if tdim > 1
......@@ -68,7 +69,7 @@ function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
This function increments current position by inc and updates xhist,
ONLY FOR GRAPH TYPE DOMAINS
"""
function increment_x!(a::Agent{Int64},p::Dict;reflected=false)
function increment_x!(a::Agent{Int64},p::Dict)
# we have to add 1 otherwise it stays on the same edge
inc = randomwalk(p["g"],get_x(a,1),Int(p["D"][1]) + 1)
# the last coef of inc corresponds to last edge reached
......
......@@ -5,9 +5,9 @@
give_birth(a::Agent,p)
Used for Gillepsie setting
"""
function give_birth(a::Agent,p::Dict,reflected)
function give_birth(a::Agent,p::Dict)
new_a = copy(a)
increment_x!(new_a,p,reflected=reflected)
increment_x!(new_a,p)
return new_a
end
......@@ -32,13 +32,13 @@ function update_afterdeath_std!(world,x_death,p::Dict) where T
end
"""
function updateWorld_G!(world,C,tspan,p)
function updateWorld_G!(world,tspan,p)
Updating rule for gillepsie setting.
Returning dt drawn from an exponential distribution with parameter the total rates of events.
# Args
tspan is for now not used but might be used for dynamic landscape
"""
function updateWorld_G!(world,p,update_rates!,tspan,reflected)
function updateWorld_G!(world,p,update_rates!,tspan)
# total update world
world_alive = skipmissing(world)
idx_world = collect(eachindex(world_alive))
......@@ -62,7 +62,7 @@ function updateWorld_G!(world,p,update_rates!,tspan,reflected)
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,reflected)
world[idx_offspring] = give_birth(mum,p)
update_afterbirth_std!(world,idx_offspring,p)
end
return dt
......
# In this file lie the function for WF algorithm
"""
function updateWorld_WF!(world,newworld,C,p,update_rates!,t,reflected)
If reflected=true, we reflect only first trait corresponding to geographic position
function updateWorld_WF!(world,newworld,C,p,update_rates!,t)
If p["reflected"]=true, we reflect only first trait corresponding to geographic position
"""
function updateWorld_WF!(world,newworld,p,update_rates!,t,reflected)
function updateWorld_WF!(world,newworld,p,update_rates!,t)
@debug "updating rates"
update_rates!(world,p,t);
# normalise to make sure that we have a probability vector
......@@ -21,6 +21,6 @@
# we introduce randomness here
@debug "incrementing offsprings traits"
for w in newworld
increment_x!(w,p,reflected=reflected)
increment_x!(w,p)
end
end
......@@ -156,10 +156,10 @@ function update_rates_std_split!(world,C,p::Dict,t::Float64)
end
end
"""
function runWorld_store_WF(p,world0;init = ([.0],),reflected=false,mode="std")
function runWorld_store_WF(p,world0;mode="std")
Wright Fisher process. Returns an array worldall with every agents.
"""
function runWorld_store_WF(p,world0;init = ([.0],),reflected=false,mode="std")
function runWorld_store_WF(p,world0;mode="std")
worldall = repeat(world0,inner = (1,length(1:Int(p["tend"]))))
N=length(world0);
newworld = copy.(world0)
......@@ -181,24 +181,24 @@ function runWorld_store_WF(p,world0;init = ([.0],),reflected=false,mode="std")
for i in 1:(Int(p["tend"])-1)
# we have to take views, otherwise it does not affect worldall
world = @view worldall[:,i];newworld = @view worldall[:,i+1];
updateWorld_WF!(world,newworld,p,update_rates!,Float64(i),reflected)
updateWorld_WF!(world,newworld,p,update_rates!,Float64(i))
end
return worldall,collect(0:Int(p["tend"]-1))
end
"""
function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
function runWorld_store_G(p,world0)
Gillepsie process. Returns a tuple worldall,tspanarray
If specified p["dt_saving"] determines the time step to save the simulation. If not only last time step is saved
"""
function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
function runWorld_store_G(p,world0)
# we store the value of world every 100 changes by default
tspan = zeros(1)
i = 1;j=0;dt = 1.
N=length(world0);
tspanarray = [];
ninit = Int(length(world0) - count(ismissing,world0))
# Array that stores the agents.
# We decide that we store agents ninit events where ninit is the initial population
# length of worldall should be changed at some point
tspanarray = zeros(1);
dt_saving = haskey(p,"dt_saving") ? p["dt_saving"] : p["tend"] + 1.
# Array that stores the agents
worldall = reshape(copy.(world0),N,1)
# we instantiate C as the biggest size it can take
update_rates_std!(skipmissing(world0),p,0.)
......@@ -212,8 +212,7 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
@info "We have reached the maximum number of individuals allowed"
break
end
# we save every ninit times
if mod(i,ninit) == 1
if tspan[end] - tspanarray[end] >= dt_saving
@info "saving world @ t = $(tspan[i])/ $(p["tend"])"
j+=1;sw = size(worldall,2);
# we use <= because we save at the end of the wile loop
......@@ -224,7 +223,7 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
push!(tspanarray,tspan[i])
end
dt = updateWorld_G!(world0,p,update_rates_std!,tspan,reflected)
dt = updateWorld_G!(world0,p,update_rates_std!,tspan)
push!(tspan, tspan[end] + dt)
i += 1
end
......
......@@ -18,7 +18,7 @@ p = Dict(
"tend" => 10.)
na_init = K0
agents = [Agent( [1e-2] .* randn(1) .- .5) for i in 1:K0]
@time worldall_test,p["tspan"] = runWorld_store_WF(p,agents,reflected=false);
@time worldall_test,p["tspan"] = runWorld_store_WF(p,agents);
## load to compare simulation
# @save "wrightfisher_test.jld2" worldall p
......
Supports Markdown
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