Commit 7208248f authored by Victor Boussange's avatar Victor Boussange
Browse files

added test files, modifying functioning of update_rates. Looks pretty good right now

parent 7b466e4d
......@@ -9,8 +9,8 @@ end
# Constructors
Agent(xhist::Array) = Agent(reshape(xhist,:,1),0.,1.)
Agent() = Agent([],0,1)
Agent(xhist::Array{T}) where T <: Number = Agent{T}(reshape(xhist,:,1),0.,1.)
Agent() = Agent(Float64[],0.,1.)
import Base.copy
copy(a::Agent) = Agent(a.x_history,a.d,a.b)
copy(m::Missing) = missing
......@@ -29,9 +29,9 @@ function new_world_G(nagents::Int,p::Dict; spread = 1., offset = 0.)
end
# returns trait i of the agent
get_x(a::Agent,i::Number) = a.x_history[Int(i):Int(i),end]
get_x(a::Agent,i::Number) = a.x_history[Int(i),end]
get_x(a::Agent) = a.x_history[:,end]
get_xhist(a::Agent,i::Number) = a.x_history[Int(i):Int(i),:]
get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
get_xhist(a::Agent) = a.x_history
get_geo(a::Agent) = sum(get_xhist(a,1))
get_d(a::Agent) = a.d
......@@ -42,11 +42,11 @@ get_fitness(a::Agent) = a.b - a.d
get_xarray(world::Array{Agent{T}},trait::Int) where T
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))
get_xarray(world::Array{Agent{T}},trait::Int) where T = reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2))
"""
function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
This function increments current position by inc and updates xhist,
increment_x!(a::Agent{Float64},p::Dict;reflected=false)
This function increments current position by inc and updates xhist,
ONLY FOR CONTINUOUS DOMAINS
"""
function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
......@@ -70,13 +70,15 @@ function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
"""
function increment_x!(a::Agent{Int64},p::Dict;reflected=false)
# we have to add 1 otherwise it stays on the same edge
inc = randomwalk(p["g"],get_x(a,1)[1],Int(p["D"][1]) + 1)
inc = randomwalk(p["g"],get_x(a,1),Int(p["D"][1]) + 1)
# the last coef of inc corresponds to last edge reached
a.x_history = hcat(a.x_history, reshape(inc[end:end],:,1));
end
"""
get_inc_reflected(x::Float64,inc::Float64,s=-1,e=1)
Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
Careful though, we do not implement reflections
"""
......
......@@ -19,7 +19,7 @@ function update_afterbirth_std!(world,idx_offspring,p::Dict) where T
a.d += α(get_x(a),x_offspring)
end
# Now updating new agent
world[idx_offspring].d = sum(α.(get_x.(skipmissing(world)),Ref(x_offspring)))
world[idx_offspring].d = sum(α.(get_x.(skipmissing(world)),Ref(x_offspring))) - α(x_offspring,x_offspring)
world[idx_offspring].b = K(x_offspring)
end
......@@ -35,6 +35,8 @@ end
function updateWorld_G!(world,C,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)
# total update world
......
......@@ -3,9 +3,9 @@
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,reflected)
function updateWorld_WF!(world,newworld,p,update_rates!,t,reflected)
@debug "updating rates"
update_rates!(world,C,p,t);
update_rates!(world,p,t);
# normalise to make sure that we have a probability vector
fitness = get_fitness.(world)
# we need to substract the minimum, otherwise fitness can be negative
......
......@@ -50,7 +50,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
# Arguments
"""
function var(world::Array{Agent};trait=1)
function var(world::Array{Agent{T}};trait=1) where T
xarray = get_xarray(world,trait)
return var(xarray,dims=1)
end
......@@ -65,7 +65,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
function covgeo(world::Array{Agent{T},1},trait = 0) where T
xarray = get_geo.(world)
if trait > 0
xstd = reshape(hcat(get_x.(world,trait)...),size(world,1),size(world,2))
xstd = reshape(get_x.(world,trait),size(world,1),size(world,2))
xarray = hcat(xarray,xstd)
end
return cov(xarray)
......
# for now we consider that competition is local within an array
"""
update_rates_std!(world,C,p::Dict,t::Float64)
update_rates_std!(world,p::Dict,t::Float64)
This standard updates takes
- competition kernels of the form α(x,y) and
- carrying capacity of the form K(x)
......@@ -11,17 +11,18 @@ function update_rates_std!(world,p::Dict,t::Float64)
traits = get_x.(world)
# traits = get_xhist.(world)
N = length(traits)
# C = SharedArray{Float64}((N,N))
D = SharedArray{Float64}(N)
# Here you should do a shared array to compute in parallel
for i in 1:(N-1)
@sync @distributed for i in 1:(N-1)
for j in i+1:N
C = α(traits[i],traits[j])
world[i].d += C
world[j].d += C
D[i] += C
D[j] += C
end
end
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
a.d = D[i]
a.b = K(traits[i])
end
end
......@@ -29,7 +30,7 @@ end
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)
aidx_e = findall(a -> get_x(a,1)==e,world)
na = length(aidx_e)
for i in aidx_e
world[i].d = na^2
......@@ -161,7 +162,6 @@ Wright Fisher process. Returns an array worldall with every agents.
function runWorld_store_WF(p,world0;init = ([.0],),reflected=false,mode="std")
worldall = repeat(world0,inner = (1,length(1:Int(p["tend"]))))
N=length(world0);
C = SharedArray{Float64}((N,N));
newworld = copy.(world0)
if mode == "std"
update_rates! = update_rates_std!
......@@ -181,7 +181,7 @@ 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 ,C,p,update_rates!,Float64(i),reflected)
updateWorld_WF!(world,newworld,p,update_rates!,Float64(i),reflected)
end
return worldall,collect(0:Int(p["tend"]-1))
end
......@@ -216,7 +216,8 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
if mod(i,ninit) == 1
@info "saving world @ t = $(tspan[i])/ $(p["tend"])"
j+=1;sw = size(worldall,2);
if sw < j
# we use <= because we save at the end of the wile loop
if sw <= j
# we double the size of worldall
worldall = hcat(worldall,Array{Missing}(missing,N,sw))
end
......@@ -227,6 +228,9 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
push!(tspan, tspan[end] + dt)
i += 1
end
@info "simulation stopped at t=$(tspan[end])"
# Saving laste time step
worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
push!(tspanarray,tspan[i])
@info "simulation stopped at t=$(tspanarray[end])"
return worldall[:,1:j],tspanarray
end
using ABMEv
using Revise,ABMEv,Test
import ABMEv:update_rates_std!
a = 0;
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
K(X) = gaussian(X[1],0.,sigma_K)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/1000
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 1000.,
"tend" => 1.5,
"NMax" => Int(10000))
na_init = 200
world0 = new_world_G(na_init,p_default,spread = .01, offset = -.25)
@time worldall,p_default["tspan"] = runWorld_store_G(p_default,world0)
# ======================================================================
using JLD2
@save "gillepsie_test.jld2" worldall,p_default
using Plots
Plots.plot(worldall,p_default)
na_init = K0
world0 = new_world_G(na_init,p_default,spread = .01)
worldall,p_default["tspan"] = runWorld_store_G(p_default,world0)
world_alive = collect(skipmissing(world0))
@testset "testing global functioning" begin
@test size(worldall,2) > 1
@test p_default["tspan"][end] >= p_default["tend"]
end
@testset "testing update rates matrix" begin
bs_end = get_b.(world_alive);ds_end = get_d.(world_alive)
update_rates_std!(world_alive,p_default,0.);
bs_recalculated = get_b.(world_alive);ds_recalculated = get_d.(world_alive);
@testset "birth coefficient" begin
for i in 1:length(bs_end)
@test bs_end[i] bs_recalculated[i];
end
end
@testset "death coefficient" begin
for i in 1:length(bs_end)
@test ds_end[i] ds_recalculated[i];
end
end
end
using ABMEv,Revise,Test
K0 = 1000; σ = 1e-1
agents1 = [Agent( [σ] .* randn(1) .- .5) for i in 1:K0]
agents2 = [Agent( [σ σ] .* randn(2) .- .5) for i in 1:K0]
world
## testing variance
@testset "Testing metrics" begin
@test first(var(agents1)) (σ).^2 atol=0.001
@test first(var(agents2,trait=2)) (σ).^2 atol=0.001
## testing covgeo
@test covgeo(agents1) (σ).^2 atol=0.001
@testset "covgeo" begin
for i in covgeo(agents1,1)
@test i (σ).^2 atol=0.001
end
end
# not sure this is the bestway of testing
# there is a problem here but we do not know hor to
@testset "covgeo2d" begin
cmat = covgeo(agents2,2);
smat = [σ^2 σ^2; σ^2 σ^2]
for i in 1:length(cmat)
@test cmat[i] smat[i] atol=0.001
end
end
end
# TODO needs to test hamming distance
......@@ -20,5 +20,38 @@ C = SharedArray{Float64}((Int(p["K0"]),Int(p["K0"])))
@btime update_afterbirth_std!(skipmissing(world0),C,1,p)
## Testing get_inc_reflected
a = Agent
@btime get_inc_reflected(get_x(a,1),p["D"][1] *randn())
using ABMEv,BenchmarkTools
a = Agent(rand(1))
@btime get_inc_reflected(get_x(a)[1],.1 *randn())
@btime rand()
@btime get_x(a)[1]
@btime get_x(a,1)
@btime get_d(a)
@btime get_xhist(a)
@btime a.x_history[1,end]
@btime get_x(a)
@btime get_xhist(a)[:,end]
##
using ABMEv,BenchmarkTools
a = 0;
sigma_K = .9;
sigma_a = .7;
K(X) = gaussian(X[1],0.,sigma_K)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/1000
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 1000.,
"NMax" => Int(10000))
na_init = 1000
world0 = new_world_G(na_init,p_default,spread = .01, offset = -.25)
tspan=zeros(1)
import ABMEv:update_rates_std!,updateWorld_G!
@btime update_rates_std!(skipmissing(world0),p_default,0.)
@btime updateWorld_G!(world0,p_default,update_rates_std!,tspan,true)
@btime update_afterbirth_std!(skipmissing(world0),1,p_default)
@btime update_afterdeath_std!(skipmissing(world0),[].8],p_default)
cd(@__DIR__)
using Revise,ABMEv
a = 0;
sigma_K = .9;
sigma_a = 1.251;
K0 = 1000;
# K(X) = gaussian(X[1],0.,sigma_K)
K(X) = 1 - 0.125 * sum(X.^2)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 1000.)
na_init = K0
agents = [Agent( [1e-2] .* randn(1) .- .5) for i in 1:K0]
@time worldall,p_default["tspan"] = runWorld_store_WF(p_default,agents,reflected=false);
# ======================================================================
using JLD2
@save "wrightfisher_test.jld2" worldall p_default
using Plots
Plots.plot(worldall,p_default,what = ["var"])
var(agents)
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