### updated version with all files in folder `examples` running

parent 16f6fd9f
Pipeline #93156 failed with stage
in 14 minutes and 39 seconds
 ... ... @@ -16,8 +16,8 @@ The user can provide **any birth and death functions**, which should depend on t ABMEv.jl provides a **numerical laboratory** for eco-evolutionary dynamics, supplying - flexible types for **individuals**, which can - evolve over any combination of space, - flexible types for **individuals**, which can - evolve over any combination of space, - store ancestors trait, - flexible types for **evolutionary spaces**, that can consist of multidimensional **discrete or continuous domains**, as well as **graphs**, - the possibility to use **callback functions** to save the state of the system at any time step ... ... @@ -42,9 +42,9 @@ Check out the documentation if you want to use the advanced features of ABMEv.jl ----- ## A first tutorial We provide here a tutorial that sums up the 5 steps necessary to launch a simulation. For the sake of the tutorial, we propose to model a population that is structured over the vertices of a graph and characterised by a trait under selection. We provide here a tutorial that sums up the 5 steps necessary to launch a simulation. For the sake of the tutorial, we propose to model a population that is structured over the vertices of a graph and characterised by a trait under selection. ### 0. Import the relevant libraries ### 0. Import the relevant libraries Let's import ABMEv.jl, and LightGraphs.jl ```julia using ABMEv ... ... @@ -52,7 +52,7 @@ using LightGraphs.jl ``` ### 1. Define the evolutionary spaces We define the geographical space as star graph with 7 vertices (i.e. the abstraction of the landscape), and a continuous trait space. We define the geographical space as star graph with 7 vertices (i.e. the abstraction of the landscape), and a continuous trait space. ```julia nodes = 7 ... ... @@ -71,10 +71,10 @@ K = 1000 # carrying capacity θ = [rand([-1,1]) for i in 1:nodes] # optimal trait values # X is the geographical position # X corresponds to the adaptive traits b(X,t) = 1 - 0.5 * (θ[X] - X)^2 d(X,Y,t) = (X ≈ Y) / K b(X,t) = max(1 - 0.5 * (θ[X] - X)^2,0.) d(X,Y,t) = (X ≈ Y) / K ``` > :warning: birth and death functions should have the same number of > :warning: birth and death functions should have the same number of arguments as above. ### 3. Define how individuals move over the evolutionary space ... ... @@ -116,7 +116,7 @@ sim = run!(w0, tend, b, d, cb = cb, cb = cb, t_saving_cb = t_saving_cb) ``` ### Plotting ... ...
 using ABMEv,UnPack println("--------------------------------") println("""TESTING GILLEPSIE with popoulation""") println("--------------------------------") for K0 in [10,50,100,500,1000,5000] myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X,t) = gaussian(X,0.,sigma_K) d(X,Y,t) = gaussian(X,Y,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=false,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @show K0 @time sim = run!(w0,Gillepsie(),tend,b,d) end println("--------------------------------") println("""TESTING CFM with population""") println("--------------------------------") for K0 in [10,50,100,500,1000,5000] myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X,t) = gaussian(X,0.,sigma_K) d(X,Y,t) = gaussian(X,Y,sigma_a)/K0 D = (1e-2,) mu = [.1] NMax = 10000 tend = 2 dm = d(,,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,b,d) end println("--------------------------------") println("""TESTING CFM with population without ancestors""") println("--------------------------------") for K0 in [10,50,100,500,1000,5000] myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X) = gaussian(X,0.,sigma_K) d(X,Y) = gaussian(X,Y,sigma_a)/K0 D = (1e-2,) mu = [.1] NMax = 10000 tend = 2 Cbar=2 p = Dict{String,Any}();@pack! p = D,mu,NMax,Cbar myagents = [Agent(myspace,(0,),ancestors=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @show K0 @time sim = run!(w0,CFM(),tend,b,d) end println("--------------------------------") println("""TESTING CFM with time""") println("--------------------------------") for tend in [1,10,50,100] K0 = 1000 myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X) = gaussian(X,0.,sigma_K) d(X,Y) = gaussian(X,Y,sigma_a)/K0 D = (1e-2,) mu = [.1] NMax = 10000 Cbar=2 p = Dict{String,Any}();@pack! p = D,mu,NMax,Cbar myagents = [Agent(myspace,(0,),ancestors=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @show tend @time sim = run!(w0,CFM(),tend,b,d) end println("--------------------------------") println("""TESTING CFM with time without Ancestors""") println("--------------------------------") for tend in [1,10,50,100] K0 = 1000 myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X) = gaussian(X,0.,sigma_K) d(X,Y) = gaussian(X,Y,sigma_a)/K0 D = (1e-2,) mu = [.1] NMax = 10000 Cbar=2 p = Dict{String,Any}();@pack! p = D,mu,NMax,Cbar myagents = [Agent(myspace,(0,)) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @show tend @time sim = run!(w0,CFM(),tend,b,d) end println("--------------------------------") println("""TESTING Gillepsie with time""") println("--------------------------------") for tend in [1,10,50,100] K0 = 1000 myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X) = gaussian(X,0.,sigma_K) d(X,Y) = gaussian(X,Y,sigma_a)/K0 D = (1e-2,) mu = [.1] NMax = 10000 # Cbar=2 p = Dict{String,Any}();@pack! p = D,mu,NMax #,Cbar myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @show tend @time sim = run!(w0,Gillepsie(),tend,b,d) end
 using Revise,ABMEv,UnPack # CFM algorithm myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; K0 = 5000 b(X,t) = gaussian(X,0.,sigma_K) d(X,Y,t) = gaussian(X,Y,sigma_a)/K0 D = (Float16(1e-1),) mu = [1.] NMax = 7000 tend = 50 bm= b(,0.); dm = d(,,0.) p = Dict{String,Any}();@pack! p = D,mu,NMax,bm,dm # myagents = [Agent(myspace,(- Float16(0.5) .+ .01 .* randn(Float16),)) for i in 1:K0] myagents = [Agent(myspace,(0.,)) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @time mysim = run!(w0,CFM(),tend,b,d,dt_saving=1.) using Plots Plots.plot(mysim) ########### GILLEPSIE K0 = 1000 myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; b(X) = gaussian(X,0.,sigma_K) d(X,Y) = gaussian(X,Y,sigma_a)/K0 D = (1e-1,) mu = [1.] NMax = 10000 tend = 100 p = Dict{String,Any}();@pack! p = D,mu,NMax #,Cbar myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @time mysim = run!(w0,Gillepsie(),tend,b,d,dt_saving=1.) Plots.plot(mysim)
 ... ... @@ -28,12 +28,12 @@ p = Dict{String,Any}();@pack! p = D,mu,NMax # definining world 0 and running myagents = [Agent(wholespace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,wholespace,p,0.) w0 = World(myagents,wholespace,p) @time sim = run!(w0,Gillepsie(),tend,b,d) ### Plotting size of the world myagents = [Agent(wholespace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,wholespace,p,0.) # we need to reinitialise the world w0 = World(myagents,wholespace,p) # we need to reinitialise the world @time sim = run!(w0,Gillepsie(),tend,dt_saving=2.,b,d) wsize = [length(w) for w in sim[:]] using Plots ... ... @@ -42,7 +42,7 @@ Plots.plot(get_tspan(sim),wsize, ylabel = "Metapopulation size", xlabel ="time", grid = false) savefig(joinpath(@__DIR__, "delta_comp_wsize.png")) # savefig(joinpath(@__DIR__, "delta_comp_wsize.png")) ## plotting position through time using Plots ... ... @@ -51,4 +51,4 @@ Plots.plot(sim, ylabel = "Geographical position", grid = false, markersize = 10) savefig(joinpath(@__DIR__, "delta_comp_pos.png")) # savefig(joinpath(@__DIR__, "delta_comp_pos.png"))
 ... ... @@ -33,5 +33,5 @@ p = Dict{String,Any}();@pack! p = D,mu,NMax # definining world 0 and running myagents = [Agent(wholespace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,wholespace,p,0.) w0 = World(myagents,wholespace,p) @time sim = run!(w0,Gillepsie(),tend,b,d)
 ... ... @@ -19,7 +19,7 @@ NMax = 2000 tend = 200 p_default = Dict{String,Any}();@pack! p_default = NMax,mu,D myagents = [Agent(myspace,(Int8(5),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) @time sim = run!(w0,Gillepsie(),tend,b,d,dt_saving=3.) ... ...
 using Revise,ABMEv,Plots,UnPack using ABMEv,Plots,UnPack nodes = 9 ... ... @@ -14,18 +14,17 @@ a = 1 b(X,t) = gaussian(X,X * a,sigma_K) / nodes d(X,Y,t) = (X ≈ Y) * gaussian(X,Y,sigma_a) / K0 NMax = 2000 D = (5e-1,5e-2) D = [5e-1,5e-2] # tend = 1.5 tend = 1500 p = Dict{String,Any}();@pack! p = NMax,mu,D myagents = [Agent(myspace,(Int8(5),Float16(5) + Float16(5e-2) * randn(Float16),),ancestors=true,rates=true) for i in 1:round(K0/nodes)] w0 = World(myagents,myspace,p,0.) w0 = World(myagents,myspace,p) s = run!(w0,Gillepsie(),tend,b,d,dt_saving=5); Plots.plot(s, ylabel = "Adaptive trait",trait=2) savefig(joinpath(@__DIR__, "gradient_adaptive_trait.png")) using Printf anim = @animate for i in 1:get_size(s) ... ... @@ -52,6 +51,3 @@ aplot = Plots.plot(thist,shistall, ylabel = "Lineages adaptive trait", ylims = (-0.5,10.5) ) savefig(aplot,joinpath(@__DIR__, "gradient_lineages_adaptive_trait.png")) # Plotting lineages
 ... ... @@ -11,15 +11,12 @@ mu = [.1] NMax = 2000 tend = 1500 p = Dict{String,Any}();@pack! p = D,mu,NMax myagents = [Agent(myspace,(1e-2 * randn(Float64),),rates=true) for i in 1:K0] myagents = [Agent(myspace,(1e-2 * randn(Float64),)) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,Gillepsie(),tend,dt_saving = 10,b,d) @time sim = run!(w0,Gillepsie(),tend,b,d, dt_saving = 10) using JLD2 @load joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim Plots.plot(sim, ylabel = "Adaptive trait", ylims = (-1,1), markersize = 2.) savefig(joinpath(@__DIR__, "sympatric_speciation.png"))
 using Revise using ABMEv,UnPack,Plots myspace = (RealSpace{1,Float64}(),) ... ... @@ -13,15 +12,12 @@ NMax = 2000 tend = 1500 dm = d(,,0.);bm = 1. p = Dict{String,Any}();@pack! p = dm,bm,D,mu,NMax myagents = [Agent(myspace,(1e-2 * randn(Float64),)) for i in 1:K0] myagents = [Agent(myspace,(1e-2 * randn(Float64)),rates=false) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,CFM(),tend,dt_saving = 10,b,d) @time sim = run!(w0,CFM(),tend,b,d,dt_saving = 10) using JLD2 @save joinpath(@__DIR__,"sim_sympatric_speciation_CFM.jld2") sim Plots.plot(sim, ylabel = "Adaptive trait", ylims = (-1,1), markersize = 2.) savefig(joinpath(@__DIR__, "sympatric_speciation_CFM.png"))
 using ABMEv,UnPack ############################### # this example implements a birth rate that is time dependent ############################## using ABMEv,UnPack using Plots ω = 2* π / 150 # angular frequency optimal_trait(t) = sin(ω * t) ... ... @@ -9,17 +13,16 @@ Plots.plot(1:tend,optimal_trait,label = "Optimal trait",xlabel = "time") myspace = (RealSpace{1,Float64}(),) sigma_K = 1.; K0 = 1000 # We will have in total 1000 individuals K0 = 300 # We will have in total 1000 individuals b(X,t) = gaussian(X,optimal_trait(t),sigma_K) d(X,Y,t) = 1/K0 D = (5e-2,) mu = [1.] NMax = 2000 p = Dict{String,Any}();@pack! p = D,mu,NMax myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) myagents = [Agent(myspace,(0,)) for i in 1:K0] w0 = World(myagents,myspace,p) @time sim = run!(w0,Gillepsie(),tend, b, d, dt_saving=3.) using Plots Plots.plot(sim, ylabel = "Adaptive trait") savefig(joinpath(@__DIR__, "time_varying_pop.png")) # savefig(joinpath(@__DIR__, "time_varying_pop.png"))
 using ABMEv using LightGraphs using Plots nodes = 7 g = star_graph(nodes) landscape = GraphSpace(g) θ = [rand([-1,1]) for i in 1:nodes] traitspace = RealSpace(1) evolspace = (landscape,traitspace) K = 150 b(X,t) = max(1 - 0.5 * (θ[X] - X)^2,0) # b(X,t) =1. d(X,Y,t) = (X ≈ Y) / K D = [nothing,5e-2] mu = [1f-1,1f-1] tend = 200. t_saving_cb = collect(range(0.,tend,length=300)) cb(w) = Dict("N" => size(w)) p = Dict("NMax" => 2000, "D" => D, "mu" => mu ) myagents = [Agent(evolspace,[rand(1:nodes),randn() * D]) for i in 1:K] w0 = World(myagents,evolspace,p) println("Running simulation with callback") @time sim = run!(w0,Gillepsie(),tend,b,d,cb=cb,t_saving_cb = t_saving_cb) Plots.plot(sim.tspan, sim["N"]) println("Running simulation with `dt_saving`") w0 = World(myagents,evolspace,p) @time sim = run!(w0,Gillepsie(),tend,b,d,dt_saving=2.0) Plots.plot(sim,trait=2)
 ... ... @@ -97,6 +97,11 @@ function get_inc(x,D,s::DiscreteSegment{T}) where {T} return round(T,_reflect1D(x,inc,s)) end function get_inc(x,D::Nothing,s::DiscreteSegment{T}) where {T} inc = rand([one(T),-one(T)]) return round(T,_reflect1D(x,inc,s)) end # normal dispersal kernel that gets truncated function get_inc(x,D::Number,s::GraphSpace{T}) where {T} niter = round(T,abs(D*randn(Float32))) + 1 ... ...
 ... ... @@ -25,12 +25,26 @@ function World(w::Vector{A},s::S,p::Dict;t::T=0.) where {A<:AbstractAgent,S<:Abs length(D) == length(s) ? nothing : throw(ArgumentError("Length of parameter D should correspond to dimension of underlying space")) SS = eltype.(s) _SS = _get_types_dim(s) for (i,_S) in enumerate(SS) if _S <: Integer # handling for discrete space is different Di = D[i] (isnothing(Di) || typeof(Di) <: AbstractFloat) ? nothing : throw(ArgumentError("`D` at dimension \$i should be whether nothing or a float")) _SS[i] = typeof(Di) end end D2 = Union{_SS...}[] for (i,Di) in enumerate(D) try push!(D2, convert(_SS[i],Di)) catch e throw(ArgumentError("`D` at dimension \$i does not match with underlying space")) # making sure that we only modify D for continuous space # for discrete space, D should be whether Nothing or AbstractFloat if SS[i] <: AbstractFloat try push!(D2, convert(_SS[i],Di)) catch e throw(ArgumentError("`D` at dimension \$i does not match with underlying space")) end else push!(D2, Di) end end p["D"] = D2 ... ...
 p = Dict(1:1:1000 .=> 1:1:1000) delete!(p,999) a = Any[rand() for i in 1:1000] a[2:2:1000] .= missing get_idx(w) = collect(eachindex(skipmissing(w))) using BenchmarkTools @btime a[get_idx(a)] @btime p @btime [true for i in 1:1000] [true for i in 1:1000] @btime sort!(values(p)) function findfreeidx(p) i = 1 while haskey(p,i) i+=1 end i end function get_xp(a,j) for (i,a) in enumerate(skipmissing(a)) if i==j return a break end end end @btime get_xp(a,998) @btime findfreeidx(p) @btime collect(keys(p)) @btime(count(ismissing,)) ############################################################ ###########Reflection on vectors of missing or not########## ############################################################ cd(@__DIR__) using Random Random.seed!(0) using LightGraphs using Test using Revise,ABMEv using UnPack,JLD2 myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; K0 = 1000; b(X) = gaussian(X,0.,sigma_K) d(X,Y) = gaussian(X,Y,sigma_a)/K0 D = (1e-2,) mu = [.1] NMax = 10000000 tend = 1.5 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax using BenchmarkTools a = Agent(myspace,(0,),ancestors=true,rates=true) myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] @btime push!(myagents,a) ###### Testing world adding agents myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) size(w0) @btime addAgent!(w0,a) size(w0) @btime removeAgent!(w0,1) ###### Testing sim adding agents using BenchmarkTools myspace = (GraphSpace(SimpleGraph(10,10)),RealSpace{1,Float64}()) myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] w1 = World(myagents,myspace,p,0.) @time sim = run!(w1,Gillepsie(),tend,dt_saving = .1) get_size(sim) # @btime add_entry!(sim,w1)
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