diff --git a/docs/src/examples/changing_environment.md b/docs/src/examples/changing_environment.md index 487b4f9192f84cc68cddcd0b7411ef2d27a25d2a..cabd1914b8995fa3445252bf0c5ba7c938e7aac9 100644 --- a/docs/src/examples/changing_environment.md +++ b/docs/src/examples/changing_environment.md @@ -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 diff --git a/docs/src/examples/delta_competition_example.md b/docs/src/examples/delta_competition_example.md index 5ab1ce5a31abaa7c16499e8da912587c3de14ffc..1152f5852ba07ae45ab4768dfbaa3bf3f0d77b70 100644 --- a/docs/src/examples/delta_competition_example.md +++ b/docs/src/examples/delta_competition_example.md @@ -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 diff --git a/docs/src/examples/dynamic_graphs.md b/docs/src/examples/dynamic_graphs.md new file mode 100644 index 0000000000000000000000000000000000000000..b2a29bef38148a859a8488d17d1419e318c033e8 --- /dev/null +++ b/docs/src/examples/dynamic_graphs.md @@ -0,0 +1,48 @@ +# 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) + diff --git a/docs/src/examples/genetic_structure.md b/docs/src/examples/genetic_structure.md index 6cb5b81ca1f2c472b4654159d27bd0ef8dab6fe4..0bcfbbce3fc9753598ff8f578698bcf617e7eaac 100644 --- a/docs/src/examples/genetic_structure.md +++ b/docs/src/examples/genetic_structure.md @@ -1,8 +1,6 @@ # [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.)  diff --git a/docs/src/examples/gradient.md b/docs/src/examples/gradient.md index c50fc6cf5e52968c1cc39a110394d5dd724f775f..7dd85a6818737925379d691f3a38eda916f735cc 100644 --- a/docs/src/examples/gradient.md +++ b/docs/src/examples/gradient.md @@ -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) diff --git a/docs/src/examples/sympatric_speciation.md b/docs/src/examples/sympatric_speciation.md index 1abcbb58b430832bd83504b1583c781697eba850..e93bcc978d107555e592d212e521f93330f87842 100644 --- a/docs/src/examples/sympatric_speciation.md +++ b/docs/src/examples/sympatric_speciation.md @@ -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) diff --git a/docs/src/index.md b/docs/src/index.md index b468e8738a1ded234fe313924907f38f38daf46a..446dfd50fdbbf978066a5946204be82186e03090 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -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. diff --git a/docs/src/manual/CFM.md b/docs/src/manual/CFM.md index 4753ba32da10ac78226bc870ab1ca4a88fccdd80..1256c21eb3e9e2c29143fc77b7cb6c1e7f581317 100644 --- a/docs/src/manual/CFM.md +++ b/docs/src/manual/CFM.md @@ -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" diff --git a/docs/src/manual/agent.md b/docs/src/manual/agent.md index 2309ae1c38d621c32ac29e0c3eaa28897a79b51b..a30ccc79c58f1092e9008a65b3e0119027098b35 100644 --- a/docs/src/manual/agent.md +++ b/docs/src/manual/agent.md @@ -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] diff --git a/docs/src/manual/callbacks.md b/docs/src/manual/callbacks.md index abb3d60dad584fc30128e86e514f27a3efe5e90d..2ee6b7a9db0d7e19a5c2e9c978849f788b79b57e 100644 --- a/docs/src/manual/callbacks.md +++ b/docs/src/manual/callbacks.md @@ -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 diff --git a/docs/src/manual/gillepsie.md b/docs/src/manual/gillepsie.md index 1ca0d5bb346cae1b166587550ebdc3924e667cc7..7930a6ef9d578cb198644114ad73b52a1779fe89 100644 --- a/docs/src/manual/gillepsie.md +++ b/docs/src/manual/gillepsie.md @@ -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] diff --git a/docs/src/notes/notes.md b/docs/src/notes/notes.md new file mode 100644 index 0000000000000000000000000000000000000000..1dc022d80dcd9ba83d1152c170929850331db658 --- /dev/null +++ b/docs/src/notes/notes.md @@ -0,0 +1,19 @@ +# 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) . diff --git a/examples/benchmark.jl b/examples/benchmark.jl index a310092fe57fc8ecaf2d2d11f042ecb2c21816bd..be205c9c6b7e39091eb2a4f61e30413f3fc6b0e4 100644 --- a/examples/benchmark.jl +++ b/examples/benchmark.jl @@ -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 diff --git a/examples/cfm.jl b/examples/cfm_vs_gillepsie.jl similarity index 62% rename from examples/cfm.jl rename to examples/cfm_vs_gillepsie.jl index 50220cc2fe6296edbc7b4d94882e6a9f7208351b..dfc5510a1f1b415fa0947b56bbbccbd7efa5524f 100644 --- a/examples/cfm.jl +++ b/examples/cfm_vs_gillepsie.jl @@ -1,6 +1,7 @@ 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) diff --git a/examples/delta_competition.jl b/examples/delta_competition.jl index 6cf66c44f9a248056a10bfe482e64a787a4ea33c..5258046d1cf256226a44c54b2c14fbe4e3396548 100644 --- a/examples/delta_competition.jl +++ b/examples/delta_competition.jl @@ -1,33 +1,40 @@ -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, diff --git a/examples/dispersal_distrib.jl b/examples/dispersal_distrib.jl deleted file mode 100644 index 1e6b3465f4fca0d67db4c6d2fa70bd0876f9121e..0000000000000000000000000000000000000000 --- a/examples/dispersal_distrib.jl +++ /dev/null @@ -1,15 +0,0 @@ -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 diff --git a/examples/dynamic_graphs.jl b/examples/dynamic_graphs.jl new file mode 100644 index 0000000000000000000000000000000000000000..15117eeec3555004f62bdf2eb53cb744e34db598 --- /dev/null +++ b/examples/dynamic_graphs.jl @@ -0,0 +1,37 @@ +# 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) diff --git a/examples/genetic_structure.jl b/examples/genetic_structure.jl index b99fc18065fb68eb4a6f72d497f0b9f8619359dc..aef728949d762fe29e6b582e931895d885a0973f 100644 --- a/examples/genetic_structure.jl +++ b/examples/genetic_structure.jl @@ -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 diff --git a/examples/gillepsie_discrete_2d.jl b/examples/gillepsie_discrete_2d.jl deleted file mode 100644 index c6d08ed7903f39ffe5dec26c35968ecd4a29f99f..0000000000000000000000000000000000000000 --- a/examples/gillepsie_discrete_2d.jl +++ /dev/null @@ -1,32 +0,0 @@ -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) diff --git a/examples/gradient.jl b/examples/gradient.jl index 4b57e211d4482aacaa009dd03e007f97f2a7131c..1648f43bfa746c0282ff7774c385c670d7458a5b 100644 --- a/examples/gradient.jl +++ b/examples/gradient.jl @@ -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")) diff --git a/examples/scaling.jl b/examples/scaling.jl deleted file mode 100644 index 84fbe5c7a38cbf1f1416e7b06513eb995ca59e81..0000000000000000000000000000000000000000 --- a/examples/scaling.jl +++ /dev/null @@ -1,27 +0,0 @@ -using ABMEv - -## Defining parameters - -sigma_K = 1.; #bandwith of resource -sigma_a = 1.2; #bandwith of competition -K0 = 1000 #carrying capacity -b(X) = gaussian(X[1],0,Float32(sigma_K)) #birth -d(X,Y) = gaussian(X[1],Y[1],Float32(sigma_a)) / Float32(K0) #competition -D = [1e-2] #mutation range -mu = [1.] #probability of mutation -NMax = 2000 #number of individual -dt_saving = 1.0 #time step saving -tend = 2. -using UnPack -p = Dict{String,Any}() -@pack! p = K,alpha,D,mu,NMax,dt_saving,tend - -## Initial conditions -agent0 = [Agent(- Float32(1.5) .+ .1 .* randn(Float32,1)) for i in 1:K0] -world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - K0))) - -## launch simulation -worldall,p["tspan"] = runWorld_store_G(p,world0); - -using Plots -Plots.plot(worldall,p) diff --git a/examples/sympatric_speciation.jl b/examples/sympatric_speciation.jl index 31bebb06f29a071d29e4bbc045e76760e6df2b05..a7345c75bd85fddfc41b26a93a9a099e795ae676 100644 --- a/examples/sympatric_speciation.jl +++ b/examples/sympatric_speciation.jl @@ -11,26 +11,15 @@ mu = [.1] NMax = 2000 tend = 1500 p = Dict{String,Any}();@pack! p = D,mu,NMax -myagents = [Agent(myspace,(1e-2 * randn(),),rates=true) for i in 1:K0] +myagents = [Agent(myspace,(1e-2 * randn(Float64),),rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,Gillepsie(),tend,dt_saving = 10,b,d) using JLD2 -@save joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim +@load joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim -Plots.plot(sim, ylabel = "Adaptive trait") +Plots.plot(sim, + ylabel = "Adaptive trait", + ylims = (-1,1), + markersize = 2.) savefig(joinpath(@__DIR__, "sympatric_speciation.png")) - -# plotting lineages -world = get_world(sim,get_size(sim)) -xhistall = get_xhist.(world[:],1) -thist = get_thist.(world[:]) -xplot = Plots.plot(thist,xhistall, - linecolor = eth_grad_std[0.], - label = "", - # title = latexstring("\\sigma_\\mu=",@sprintf("%1.2f",world.p["D"][2][1]),", \\sigma_D=",@sprintf("%1.2f",world.p["D"][1])), - grid = false, - xlabel = "time", - ylabel = "Historical adaptive trait" - ) -savefig(joinpath(@__DIR__, "x_hist_sympatric_speciation.png")) diff --git a/examples/sympatric_speciation_CFM.jl b/examples/sympatric_speciation_CFM.jl index 519cc83d194cb165d2db5c55ad4fac5dada59cf1..1e42f4778c20adf65564c1f03b8d8027ea7d565c 100644 --- a/examples/sympatric_speciation_CFM.jl +++ b/examples/sympatric_speciation_CFM.jl @@ -13,26 +13,15 @@ NMax = 2000 tend = 1500 dm = d([0],[0],0.);bm = 1. p = Dict{String,Any}();@pack! p = dm,bm,D,mu,NMax -myagents = [Agent(myspace,(1e-2 * randn(),)) 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,CFM(),tend,dt_saving = 20,b,d) +@time sim = run!(w0,CFM(),tend,dt_saving = 10,b,d) using JLD2 @save joinpath(@__DIR__,"sim_sympatric_speciation_CFM.jld2") sim -Plots.plot(sim, ylabel = "Adaptive trait") +Plots.plot(sim, + ylabel = "Adaptive trait", + ylims = (-1,1), + markersize = 2.) savefig(joinpath(@__DIR__, "sympatric_speciation_CFM.png")) - -# plotting lineages -world = get_world(sim,get_size(sim)) -xhistall = get_xhist.(world[:],1) -thist = get_thist.(world[:]) -xplot = Plots.plot(thist,xhistall, - linecolor = eth_grad_std[0.], - label = "", - # title = latexstring("\\sigma_\\mu=",@sprintf("%1.2f",world.p["D"][2][1]),", \\sigma_D=",@sprintf("%1.2f",world.p["D"][1])), - grid = false, - xlabel = "time", - ylabel = "Historical adaptive trait" - ) -savefig(joinpath(@__DIR__, "x_hist_sympatric_speciation_CFM.png")) diff --git a/examples/time_varying.jl b/examples/time_varying.jl index fab485de837963e112d45a1376a7f78270673577..d316b55a6fdbcaf0505807b0510ea32ffdab59b9 100644 --- a/examples/time_varying.jl +++ b/examples/time_varying.jl @@ -15,10 +15,10 @@ d(X,Y,t) = 1/K0 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.) using Plots Plots.plot(sim, ylabel = "Adaptive trait") diff --git a/src/ABMEv_Space.jl b/src/ABMEv_Space.jl index 62cb57e17cb5b5762aeb4f7647b60d31d6d31033..05c2f240e3c60530e746d361d3c2e858cdcdefc5 100644 --- a/src/ABMEv_Space.jl +++ b/src/ABMEv_Space.jl @@ -111,6 +111,9 @@ abstract type AbstractDynSpace{Dim,T<:Number} <: AbstractSpace{Dim,T,IsFinite{tr """ \$(TYPEDEF) A dynamic graph space. + +# Example +DynGraphSpace(g,f) Function f(t) takes as argument time, and returns the index of the graph to pick at time t from array g """ struct DynGraphSpace{T<:Number} <: AbstractDynSpace{1,T} diff --git a/test/dyn_land.jl b/test/dyn_land.jl index ce7fba14661be15b2d8d93d15ea7e5bc30eb38cf..5397717d2dcc87067aa7cb9b95849ed4b59cad8e 100644 --- a/test/dyn_land.jl +++ b/test/dyn_land.jl @@ -35,11 +35,11 @@ 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,(1,),ancestors=true,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) w1 = copy(w0) @info "Running simulation with Gillepsie algorithm" -sim = run!(w1,Gillepsie(),tend) +sim = run!(w1,Gillepsie(),tend,b,d) @test !(isnothing(sim)) diff --git a/test/gillepsie.jl b/test/gillepsie.jl index bceafe18a32d965d3f7d19af4a86e018d2a3f2df..5ffdfe0ba29293450b673466a41e38c3570b3fce 100644 --- a/test/gillepsie.jl +++ b/test/gillepsie.jl @@ -16,13 +16,13 @@ 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) @info "Running simulation with Gillepsie algorithm" -@time sim = run!(w1,Gillepsie(),tend) +@time sim = run!(w1,Gillepsie(),tend,b,d) @test typeof(sim) <: Simulation # @save "gillepsie_test.jld2" world_alive @@ -42,11 +42,11 @@ w1 = copy(w0) @testset "Testing update rates matrix" begin @testset "birth event" begin myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] - w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie()) + w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie(),b,d,) mum_idx = 1 - updateBirthEvent!(w0,Gillepsie(),1) + updateBirthEvent!(w0,Gillepsie(),1,b,d) bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1)) - update_rates!(w1,Gillepsie()); + update_rates!(w1,Gillepsie(),b,d); bs_recalculated = get_b.(agents(w1));ds_recalculated = get_d.(agents(w1)) @test prod(bs_end .≈ bs_recalculated) @test prod(ds_end .≈ ds_recalculated) @@ -55,11 +55,11 @@ w1 = copy(w0) @testset "death event" begin myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0] - w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie()) + w1 = World(myagents,myspace,p,0.);update_rates!(w1,Gillepsie(),b,d) mum_idx = 1 - updateDeathEvent!(w0,Gillepsie(),1) + updateDeathEvent!(w0,Gillepsie(),1,d) bs_end = get_b.(agents(w1));ds_end = get_d.(agents(w1)) - update_rates!(w1,Gillepsie()); + update_rates!(w1,Gillepsie(),b,d); bs_recalculated = get_b.(agents(w1));ds_recalculated = get_d.(agents(w1)) @test prod(bs_end .≈ bs_recalculated) @test prod(ds_end .≈ ds_recalculated) diff --git a/test/metrics_hamming.jl b/test/metrics_hamming.jl index 67ae4f02e7e98a9db532b3b80f6e07f39c49c868..3984a793e9a556e1e9b15b85b4995895a62fab20 100644 --- a/test/metrics_hamming.jl +++ b/test/metrics_hamming.jl @@ -16,12 +16,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.) @info "Running simulation with Gillepsie algorithm" -@time sim = run!(w0,Gillepsie(),tend) +@time sim = run!(w0,Gillepsie(),tend,b,d) @testset "Hamming distances" begin @test typeof(get_xhist_mat(agents(w0))[1] )<: Array diff --git a/test/profiling_ABMEv.jl b/test/profiling_ABMEv.jl index e075a7f6813ec8e30f2160d623d54fdbb912b808..7cc6905a3def8bab98f6271d23deea778a315c46 100644 --- a/test/profiling_ABMEv.jl +++ b/test/profiling_ABMEv.jl @@ -101,16 +101,16 @@ 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.) @info "Running simulation with Gillepsie algorithm" -@time sim = run!(w0,Gillepsie(),tend) +@time sim = run!(w0,Gillepsie(),tend,b,d) myagents = [Agent(myspace,(0,),ancestors=false,rates=true) for i in 1:K0] w0 = World(myagents,myspace,p,0.) @info "Running simulation with Gillepsie algorithm" -@time sim = run!(w0,Gillepsie(),tend) +@time sim = run!(w0,Gillepsie(),tend,b,d) agentarray = vcat(copy.(agents(w0)),Array{Missing}(missing,NMax - size(w0),1)) diff --git a/test/simulation.jl b/test/simulation.jl index 50e850266225a9d3cb0615d27b323f076808bedb..48e5ab222a4781a5f712f8a0f16d39670c55c8f8 100644 --- a/test/simulation.jl +++ b/test/simulation.jl @@ -39,14 +39,14 @@ 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) cb = (names = ["gamma_div"], agg = Function[w -> var(Float64.(get_x(w,1)))]) eltype(cb.agg) -@time sim = run!(w1,Gillepsie(),tend,cb=cb,dt_saving = .1) +@time sim = run!(w1,Gillepsie(),tend,b,d,cb=cb,dt_saving = .1) @test typeof(sim["gamma_div"]) <: Vector diff --git a/test/world.jl b/test/world.jl index f314a1df9469ac9bfb4f9c076e1dd419944b8d10..7fd180a04f4d320960e00081f3551f30883e10d1 100644 --- a/test/world.jl +++ b/test/world.jl @@ -29,4 +29,4 @@ removeAgent!(w,11) @test size(get_x(w,:)) == (10,2) ############## Testing world with Gillepsie -@test !isnothing(updateWorld!(w,Gillepsie())) +@test !isnothing(updateWorld!(w,Gillepsie(),b,d))