Commit 68ea8cfa by Victor

parent f80c1a7b
Pipeline #78102 passed with stage
in 20 minutes and 47 seconds
 ... ... @@ -22,10 +22,10 @@ In this tutorial, we define a birth function that is time dependent. This can b D = (5e-2,) mu = [1.] NMax = 2000 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax 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.) @time sim = run!(w0,Gillepsie(),tend,dt_saving=3.) @time sim = run!(w0, Gillepsie(), tend, b, d, dt_saving=3.)  ## Plotting ... ...
 ... ... @@ -74,12 +74,12 @@ We keep track of individuals' ancestors by setting ancestors=true. Because we using UnPack# useful macro @pack! NMax = 2000 tend = 300. p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax p = Dict{String,Any}();@pack! p = D,mu,NMax myagents = [Agent(myspace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,Gillepsie(),tend) @time sim = run!(w0,Gillepsie(), b, d, tend)  This is the simplest run you can do. Now time to more interesting things This is the simplest run you can do. Now is time for more interesting things ## Analysis ... ...
 # Dynamic graph in ABMEv.jl, individuals can evolve over graphs which connectivity can change over time. To define such a dynamic graph, one needs to use the constructor DynGraphSpace. @docs DynGraphSpace  Here is an example julia using UnPack,ABMEv,LightGraphs nodes = 10 g1 = LightGraphs.grid(Int8[9,1]) g2 = SimpleGraph(Int8(9)) g = [g1,g2] # array of graphs # This is the function to implement DynGraphSpace. # Note that it returns indices function update_g(t) T = 100 if sin(t/T*2*π) > 0 1 # use graph g1 else 2 # use graph g2 end end mydynamicgraph = DynGraphSpace(g,update_g) wholespace = (mydynamicgraph,) # Definition of birth and death rate K0 = 1000 # We will have in total 1000 individuals b(X,t) = 1 / nodes d(X,Y,t) = (X[1] ≈ Y[1]) / K0 # Mutation / dispersal parameters mu = [1.] D = (1.5,) # maximum size, tend NMax = 2000 tend = 300. # wrapping up all the parameters 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.) @time sim = run!(w0,Gillepsie(),tend,b,d) 
 # [Modelling agents with a genetic structure](@id genetic_structure) In this example, we show how to model a population which evolve on a linear geographic space, and is defined by a genotype graph. Any two connected node in the genotype graph should be thought of as two neighbour genomes, i.e that are very similar in their alleles / nucleotide. The genotype space is inspired from the article [The architecture of an empirical genotype-phenotype map](http://doi.wiley.com/10.1111/evo.13487) In this example, we show how to model a population evolving over a linear geographic space. The population is also defined by a genotype, which corresponds to the node of a **genotype graph**. The concept of genotype graph is better explained in [The architecture of an empirical genotype-phenotype map](http://doi.wiley.com/10.1111/evo.13487). Any two connected nodes in the genotype graph should be thought of as two related genomes, i.e that are very similar in their alleles / nucleotide. As such, if the genotype of an offspring differs from its parent, it will be located in the neighbourhood of it parent's genotype node. ## Defining the space julia ... ... @@ -29,7 +27,7 @@ The genotype space is inspired from the article [The architecture of an empirica NMax = 2000 # tend = 1.5 tend = 3000 p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu p_default = Dict{String,Any}();@pack! p_default = NMax,mu myagents = [Agent(myspace,(rand(Int8(1):Int8(nodes)),initnode),ancestors=true,rates=true) for i in 1:round(K0/nodes)] w0 = World(myagents,myspace,p_default,0.)  ... ...
 ... ... @@ -29,13 +29,13 @@ NMax = 2000 D = (5e-1,5e-2) # tend = 1.5 tend = 1500 p = Dict{String,Any}();@pack! p = d,b,NMax,mu,D 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.) s = run!(w0,Gillepsie(),tend,dt_saving=5); s = run!(w0,Gillepsie(),tend,b,d,dt_saving=5); Plots.plot(s, ylabel = "Adaptive trait",trait=2)  ![](../assets/tutorials/gradient_adaptive_trait.png) ... ...
 ... ... @@ -17,10 +17,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,b,d,dt_saving = 10) Plots.plot(sim, ylabel = "Adaptive trait")  ![](../assets/tutorials/sympatric_speciation.png)
 ... ... @@ -20,7 +20,7 @@ using ABMEv  ## Tutorial We strongly advise to have a look at the tutorial section. We strongly advise to have a look at the tutorial section. All the scripts of the examples can be found [here](https://gitlab.ethz.ch/bvictor/abmev/-/tree/master/examples). @contents Pages = [ "examples/delta_competition_example.md", ... ... @@ -52,5 +52,5 @@ As of now, three types of simulation algorithm can be used: - [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632) - [Champagnat and Ferriere second article - 2008](https://www.tandfonline.com/doi/full/10.1080/15326340802437710) ## Similar packages: [Agents.jl](https://juliadynamics.github.io/Agents.jl/) This package is oriented towards general ABM modelling, and miss many of the specificities of ABMEv.jl , such as vector spaces operations, update algorithm, and analysis tools. ## Similar packages [Agents.jl](https://juliadynamics.github.io/Agents.jl/) This package is oriented towards general ABM modelling, and thus is not as efficient and easy to deploy as ABMEv.jl for simulating stochastic models of structured populations.
 ... ... @@ -20,10 +20,10 @@ D = (1e-2,) mu = [.1] NMax = 2000 tend = 1500 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar p = Dict{String,Any}();@pack! p = D,mu,NMax,Cbar 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, b, d, dt_saving = 4)  !!! warning "Development" ... ...
 ... ... @@ -20,7 +20,7 @@ end The type Agent has two important composite types - Ancestors{bool} : when bool = true, the ancestors traits are stored, - Rates{bool} : when bool = true, the rates d and b of agents are updated at each time step. This is need in e.g. Gillepsie Algorithm - Rates{bool} : when bool = true, the rates d and b of agents are updated at each time step. This is needed in e.g. Gillepsie Algorithm @autodocs Modules = [ABMEv] ... ...
 ... ... @@ -29,12 +29,12 @@ D = (1e-2,) mu = [.1] NMax = 10000 tend = 1.5 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax 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.) w1 = copy(w0) @time sim = run!(w1,Gillepsie(),tend,cb=cb,dt_saving = .1) @time sim = run!(w1,Gillepsie(),tend,b,d,cb=cb,dt_saving = .1)  ## Accessing the callbacks ... ...
 ... ... @@ -11,25 +11,8 @@ math R(t) = \left[ \sum_i b_i(t) + d_i(t) \right]  A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate R A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate R. ### Time steps An event is exponentiallly distributed in time, with parameter \lambda = U(t). This makes events memoryless, meaning that the probability of having a birth or death event is always the same, no matter when (P(X > s_t | X > t) = P(X > s) . @autodocs Modules = [ABMEv] ... ...
 # Taken from gillepsie algorithm documentation !!! tip "Inversion method" Let B(t) = \sum_i b_i(t) and  D(t) = \sum_i d_i(t). Let T_b, T_d the time for a birth or death event to occur. Then we have P(T_b < T_d) = \frac{B(t)}{B(t) + D(t)} (competing exponentials). Let U be an \mathcal{U}_{(0,1)}-distributed random variable and F \colon \R \to [0,1] be a distribution function. Then we have math P(I_F(U) \leq x ) = P(U \leq F(x)) = F(x)  Thanks to the ***inversion method*** we get the incremental time step dt, exponentially distributed with parameter \lambda = R(t), as math dt(\omega) = -\frac{\log(U(\omega))}{R(t)} \iff X(\omega) = \exp(-U(t)dt(\omega))  ### Time steps An event is exponentiallly distributed in time, with parameter \lambda = U(t). This makes events memoryless, meaning that the probability of having a birth or death event is always the same, no matter when (P(X > s_t | X > t) = P(X > s) .
 ... ... @@ -56,11 +56,11 @@ for K0 in [10,50,100,500,1000,5000] NMax = 10000 tend = 2 Cbar=2 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar 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) @time sim = run!(w0,CFM(),tend,b,d) end println("--------------------------------") ... ... @@ -78,11 +78,11 @@ for tend in [1,10,50,100] mu = [.1] NMax = 10000 Cbar=2 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar 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) @time sim = run!(w0,CFM(),tend,b,d) end println("--------------------------------") ... ... @@ -100,11 +100,11 @@ for tend in [1,10,50,100] mu = [.1] NMax = 10000 Cbar=2 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar 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) @time sim = run!(w0,CFM(),tend,b,d) end println("--------------------------------") ... ... @@ -121,9 +121,9 @@ for tend in [1,10,50,100] mu = [.1] NMax = 10000 # Cbar=2 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax #,Cbar 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) @time sim = run!(w0,Gillepsie(),tend,b,d) end
 using Revise,ABMEv,UnPack myspace = (RealSpace{1,Float16}(),) # CFM algorithm myspace = (RealSpace{1,Float64}(),) sigma_K = .9; sigma_a = .7; K0 = 5000 ... ... @@ -10,16 +11,15 @@ D = (Float16(1e-1),) mu = [1.] NMax = 7000 tend = 50 Cbar= b([0],0.) + d([0],[0],0.) p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar bm= b([0],0.); dm = d([0],[0],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,dt_saving=1.) @time mysim = run!(w0,CFM(),tend,b,d,dt_saving=1.) x,t = get_xnt(mysim,trait=1) using Plots Plots.scatter(t,x,color =:blue,labels="" ) Plots.plot(mysim) ########### GILLEPSIE ... ... @@ -33,11 +33,9 @@ D = (1e-1,) mu = [1.] NMax = 10000 tend = 100 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax #,Cbar 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.) @time mysim = run!(w0,Gillepsie(),tend,dt_saving=1.) x,t = get_xnt(sim,trait=1) using Plots Plots.scatter(t,x,color =:blue,labels="" ) Plots.plot(mysim)
 using Revise,ABMEv # Space as a discrete segement using UnPack,ABMEv nodes = 10 mysegment = DiscreteSegment(1,nodes) using ABMEv, LightGraphs # other possibility: defining the space as a graph nodes = 10 g = grid([nodes,1]) g = LightGraphs.grid([nodes,1]) # cf LightGraphs.jl for options to generate a graph mysegmentgraph = GraphSpace(g) wholespace = (mysegmentgraph,) using GraphPlot # plotting the graph gplot(g, collect(1:nodes), collect(1:nodes)) # Definition of birth and death rate K0 = 1000 # We will have in total 1000 individuals b(X,t) = 1 / nodes d(X,Y,t) = (X[1] ≈ Y[1]) / K0 # Mutation / dispersal parameters mu = [1.] D = (1.5,) using UnPack # maximum size, tend NMax = 2000 tend = 300. p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax # wrapping up all the parameters 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.) @time sim = run!(w0,Gillepsie(),tend) @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 @time sim = run!(w0,Gillepsie(),tend,dt_saving=2.) @time sim = run!(w0,Gillepsie(),tend,dt_saving=2.,b,d) wsize = [length(w) for w in sim[:]] using Plots Plots.plot(get_tspan(sim),wsize, ... ...
 using LightGraphs using Test using Revise,ABMEv mysegmentgraph = GraphSpace(grid([10,1])) mysegment = DiscreteSegment(1,10) myplot = Plots.plot(grid = false) for D in [1e-3,1e-2,1e-1,1e0,1e1] graphinc = [get_inc(5,D,mysegmentgraph) for i in 1:10000] seginc = [get_inc(5,D,mysegment) for i in 1:10000] using StatsPlots StatsPlots.density!(graphinc,label="Dispersal distribution for graph segment") StatsPlots.density!(seginc,label="Dispersal distribution for regular segment",linestyle=:dash) end myplot
 # Dynamic graph using UnPack,ABMEv,LightGraphs nodes = 10 g1 = LightGraphs.grid(Int8[9,1]) g2 = SimpleGraph(Int8(9)) g = [g1,g2] # array of graphs # This is the function to implement DynGraphSpace. # Note that it returns indices function update_g(t) T = 100 if sin(t/T*2*π) > 0 1 # use graph g1 else 2 # use graph g2 end end mydynamicgraph = DynGraphSpace(g,update_g) wholespace = (mydynamicgraph,) # Definition of birth and death rate K0 = 1000 # We will have in total 1000 individuals b(X,t) = 1 / nodes d(X,Y,t) = (X[1] ≈ Y[1]) / K0 # Mutation / dispersal parameters mu = [1.] D = (1.5,) # maximum size, tend NMax = 2000 tend = 300. # wrapping up all the parameters 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.) @time sim = run!(w0,Gillepsie(),tend,b,d)
 ... ... @@ -17,11 +17,11 @@ D = (3e-1,5e-1) NMax = 2000 # tend = 1.5 tend = 200 p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu,D 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.) @time sim = run!(w0,Gillepsie(),tend,dt_saving=3.) @time sim = run!(w0,Gillepsie(),tend,b,d,dt_saving=3.) using GraphPlot,StatsBase ... ...
 using ABMEv,UnPack cd(@__DIR__) using Dates using JLD2 using Random;Random.seed!(0) a = 1; D = [1e0;a]; mu = [1.;.001]; sigma_K = 1.; sigma_a = .8; nodes = 10 K0 = 1000 K(X) = gaussian(X[2],X[1] * a,sigma_K) / nodes alpha(X,Y) = (X[1] ≈ Y[1]) * gaussian(X[2],Y[2],sigma_a) / K0 NMax = 5000 tend = 1. reflected = true; dt_saving = 15. p = Dict{String,Any}() @pack! p = alpha,K,D,NMax,tend,reflected,nodes, dt_saving,mu agent0 = [Agent{MixedAgent}( Float32[rand(1:p["nodes"]), 1e-2* randn() + (p["nodes"] + 1 ) * a/2 ] ) for i in 1:K0] world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - K0))) worldall,p["tspan"] = runWorld_store_G(p,world0); # @save "worldall_discrete_a=$(a)_mu2=$(mu[2])_D2=$(D[2])_D1=$(D[1]).jld2" using Plots Plots.plot(worldall,p,trait=2)
 ... ... @@ -17,13 +17,13 @@ NMax = 2000 D = (5e-1,5e-2) # tend = 1.5 tend = 1500 p = Dict{String,Any}();@pack! p = d,b,NMax,mu,D 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.) s = run!(w0,Gillepsie(),tend,dt_saving=5); 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")) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!