Commit 8ec7f87f authored by Victor Boussange's avatar Victor Boussange
Browse files

Merge remote-tracking branch 'origin/no_C_matrix'

* origin/no_C_matrix:
  uploading jld2 wrightfisher unit test
  uploading gillepsie jld2 test
  modifying project.toml and unit test
  Adding many unit test files
  added test files, modifying functioning of update_rates. Looks pretty good right now
  no message
  significantly modifying ABMEv_Gillepse

# Conflicts:
#	src/ABMEv_runworld.jl
parents 1f927406 ab835c9c
......@@ -21,11 +21,29 @@ version = "3.5.0+2"
[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
[[BinaryProvider]]
deps = ["Libdl", "SHA"]
git-tree-sha1 = "5b08ed6036d9d3f0ee6369410b830f8873d4024c"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.8"
[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "0becdab7e6fbbcb7b88d8de5b72e5bb2f28239f3"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "0.5.8"
[[CodecZlib]]
deps = ["BinaryProvider", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "05916673a2627dd91b4969ff8ba6941bc85a960e"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.6.0"
[[CompilerSupportLibraries_jll]]
deps = ["Libdl", "Pkg"]
git-tree-sha1 = "ff8101d6736414bc93c0f8df77b1e4095ca988c3"
git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612"
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.3.2+0"
version = "0.3.3+0"
[[DataAPI]]
git-tree-sha1 = "674b67f344687a88310213ddfa8a2b3c76cc4252"
......@@ -56,6 +74,15 @@ git-tree-sha1 = "c4ed10355637fcb0725dc6a27060f74df24f13cd"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.23.2"
[[FileIO]]
deps = ["Pkg"]
git-tree-sha1 = "3d7cb2c4c850439f19c4d6d3fbe1dce6481cddb1"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.2.4"
[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "3eb5253af6186eada40de3df524a1c10f0c6bfa2"
......@@ -71,7 +98,20 @@ version = "0.1.2"
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
[[JLD2]]
deps = ["CodecZlib", "DataStructures", "FileIO", "Mmap", "Pkg", "Printf", "UUIDs"]
git-tree-sha1 = "d6cfa7c24e27d7eaa2290372739c8298257dae16"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.1.12"
[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "2eadbbde5534346cbb837c3a75b377cba477a06d"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.7.13"
[[LibGit2]]
deps = ["Printf"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
[[Libdl]]
......@@ -90,6 +130,12 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "1c41621653250b2824b6e664ac5bd805558aeff9"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "0.4.3"
[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "f7d2e3f654af75f01ec49be82c231c382214223a"
......@@ -134,7 +180,7 @@ uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.9.12"
[[Pkg]]
deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Test", "UUIDs"]
deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
[[Printf]]
......@@ -166,6 +212,12 @@ git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "0.2.0"
[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "UUIDs", "Unicode"]
git-tree-sha1 = "6cefbc0e3b62146e564a3cf209e7370a839883da"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "2.6.0"
[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "86c5647b565873641538d8f812c04e4c9dbeb370"
......@@ -243,6 +295,12 @@ uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.5"
[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
......
......@@ -6,9 +6,14 @@ version = "0.1.0"
[deps]
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
......@@ -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,22 +42,23 @@ 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)
tdim = length(p["D"])
if reflected
inc = [get_inc_reflected(get_x(a)[1],p["D"][1] *randn())]
inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
if tdim > 1
inc = vcat(inc,rand.(Binomial.(1,p["mu"][2:end])) .* p["D"][2:end] .* randn(tdim-1))
inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
end
else
inc = rand.(Binomial.(1,p["mu"][:])) .* p["D"][:] .* randn(length(tdim))
# inc = yes no mutation * mutation
inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
end
a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
end
......@@ -69,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
"""
......
......@@ -11,33 +11,23 @@ function give_birth(a::Agent,p::Dict,reflected)
return new_a
end
function update_afterbirth_std!(world,C,idx::Int,p::Dict) where T
traits = get_x.(world)
function update_afterbirth_std!(world,idx_offspring,p::Dict) where T
# updating competition only the two columns corresponding to agent idx
α = p["alpha"];K=p["K"]
for i in 1:length(traits)
C[i,idx] = α(traits[i],traits[idx])
C[idx,i] = C[i,idx]
end
# updating death rate only the two columns corresponding to agent idx
for (i,a) in enumerate(world)
a.d += C[i,idx]
α = p["alpha"];K=p["K"];
x_offspring = get_x(world[idx_offspring])
for a in skipmissing(world)
a.d += α(get_x(a),x_offspring)
end
# Now updating new agent
world[idx].d = sum(C[idx,:])
world[idx].b = K(traits[idx])
world[idx_offspring].d = sum(α.(get_x.(skipmissing(world)),Ref(x_offspring))) - α(x_offspring,x_offspring)
world[idx_offspring].b = K(x_offspring)
end
function update_afterdeath_std!(world,C,idx::Int,p::Dict) where T
traits = get_x.(world)
function update_afterdeath_std!(world,x_death,p::Dict) where T
α = p["alpha"]
# updating death rate only the two columns corresponding to agent idx
for (i,a) in enumerate(world)
a.d -= C[i,idx]
end
# updating competition only the two columns corresponding to agent idx
for i in 1:length(traits)
C[i,idx] = .0
C[idx,i] = .0
for a in skipmissing(world)
a.d -= α(get_x(a),x_death)
end
end
......@@ -45,8 +35,10 @@ 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,C,p,update_rates!,tspan,reflected)
function updateWorld_G!(world,p,update_rates!,tspan,reflected)
# total update world
world_alive = skipmissing(world)
idx_world = collect(eachindex(world_alive))
......@@ -57,18 +49,21 @@ function updateWorld_G!(world,C,p,update_rates!,tspan,reflected)
events_weights = ProbabilityWeights(vcat(get_d.(world_alive),get_b.(world_alive)))
i_event = sample(events_weights)
# This is ok since length is a multiple of 2
I = Int(length(events_weights) / 2)
if i_event <= I
N = length(world) - count(ismissing,world)
if i_event <= N
# DEATH EVENT
# In this case i_event is also the index of the individual to die in the world_alive
idx_offspring = idx_world[i_event]
x_death = get_x(world[idx_offspring])
update_afterdeath_std!(world,x_death,p)
world[idx_offspring] = missing
update_afterdeath_std!(world_alive,C,idx_offspring,p)
else
# birth event
idx_offspring = findfirst(ismissing,world)
world[idx_offspring] = give_birth(world[idx_world[i_event-I]],p,reflected)
update_afterbirth_std!(world_alive,C,idx_offspring,p)
# 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)
update_afterbirth_std!(world,idx_offspring,p)
end
return dt
else
......
......@@ -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)
"""
function update_rates_std!(world,C,p::Dict,t::Float64)
function update_rates_std!(world,p::Dict,t::Float64)
α = p["alpha"];K=p["K"];
# Competition matrix
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
@sync @distributed for i in 1:(N-1)
for j in i+1:N
C[i,j] = α(traits[i],traits[j])
C[j,i] = C[i,j]
C = α(traits[i],traits[j])
D[i] += C
D[j] += C
end
end
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
# we only update death rate
# this could be imptrove since \alpha is symmetric, by using a symmetric matrix
a.d = sum(C[i,:])
# /!| not ideal to assign at every time step the birth rate that is constant
a.d = D[i]
a.b = K(traits[i])
end
end
......@@ -33,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
......@@ -165,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!
......@@ -185,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
......@@ -205,9 +201,23 @@ function runWorld_store_G(p,world0;init = ([.0],),reflected=false)
# length of worldall should be changed at some point
worldall = reshape(copy.(world0),N,1)
# we instantiate C as the biggest size it can take
<<<<<<< HEAD
C = SharedArray{Float64}((N,N))
update_rates_std!(skipmissing(world0),C,p,0.)
while tspan[i]<p["tend"] && count(ismissing,world0) < p["NMax"] && count(ismissing,world0) > 0
=======
update_rates_std!(skipmissing(world0),p,0.)
while tspan[i]<p["tend"]
if dt < 0
throw("We obtained negative time step dt = $dt at event $i")
elseif count(ismissing,world0) == p["NMax"]
@info "All individuals have died :("
break
elseif count(ismissing,world0) == 0
@info "We have reached the maximum number of individuals allowed"
break
end
>>>>>>> origin/no_C_matrix
# we save every ninit times
if dt < 0
throw("We obtained negative time step dt = $dt at event $i")
......@@ -215,16 +225,21 @@ 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
worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
push!(tspanarray,tspan[i])
end
dt = updateWorld_G!(world0,C,p,update_rates_std!,tspan,reflected)
dt = updateWorld_G!(world0,p,update_rates_std!,tspan,reflected)
push!(tspan, tspan[end] + dt)
i += 1
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
cd(@__DIR__)
Random.seed!(0)
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)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 1.5,
"NMax" => Int(10000))
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_test = collect(skipmissing(worldall[:,end]))
# @save "gillepsie_test.jld2" world_alive
@load "gillepsie_test.jld2" world_alive
## Testing
@testset "Gillepsie Algorithm" begin
@testset "Testing global functioning" begin
@test size(worldall,2) > 1
@test p_default["tspan"][end] >= p_default["tend"]
end
## Comparing simulation
xarray = get_xarray(world_alive,1);xarray_test = get_xarray(world_alive_test,1);
@test xarray xarray_test
@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);
@test bs_end bs_recalculated;
@test ds_end ds_recalculated;
end
end
K0 = 1000; σ = 1e-1
agents1 = [Agent( [σ] .* randn(1) .- .5) for i in 1:K0]
agents2 = [Agent( [σ σ] .* randn(2) .- .5) for i in 1:K0]
## 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
@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
## Testing parallel implementation
using Random;Random.seed!(0)
cd(@__DIR__)
using Distributed;addprocs(exeflags="--project")
using Test
@everywhere begin
using ABMEv
sigma_a = 1.251;
K0 = 1000;
K(X) = 1 - 0.125 * sum(X.^2)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
end
p = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 10.)
na_init = K0
agents = [Agent( [1e-2] .* randn(1) .- .5) for i in 1:K0]
# @test nprocs() == 5
# time should be less than 1 sec on the second run
# using BenchmarkTools
@time (worldall_test,p["tspan"]) = runWorld_store_WF(p,agents,reflected=false)
@test size(worldall_test,2) == 10
## load to compare simulation
# using JLD2
# @save "wrightfisher_test.jld2" worldall p
# @load "wrightfisher_test.jld2" worldall
# xarray = get_xarray(worldall,1); xarray_test = get_xarray(worldall_test,1)
# @test xarray ≈ xarray_test
......@@ -2,6 +2,7 @@ using Distributed;addprocs()
@everywhere push!(LOAD_PATH,homedir()*"/ETHZ/projects/ABMEv.jl/src")
@everywhere using ABMEv,BenchmarkTools,SharedArrays
## Testing update_afterbirth_std!
p= Dict("K0" => 1000.,
"D" => [1e-2 - 1e-3],
"mu" => [.1],
......@@ -17,3 +18,64 @@ world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - na_init)))
C = SharedArray{Float64}((Int(p["K0"]),Int(p["K0"])))
@btime update_afterbirth_std!(skipmissing(world0),C,1,p)
## Testing get_inc_reflected
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)
## parallelisation Gillepsie
# For now Gillepsie can not really be parallelised
using Distributed;addprocs(exeflags="--project")
@everywhere begin
using ABMEv
sigma_a = 1.251;
K0 = 1000;
K(X) = 1 - 0.125 * sum(X.^2)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
end