diff --git a/docs/src/assets/tutorials/optimal_trait.png b/docs/src/assets/tutorials/optimal_trait.png new file mode 100644 index 0000000000000000000000000000000000000000..a8a32973c130495ae28ca6f56be6876a7d15fcd1 Binary files /dev/null and b/docs/src/assets/tutorials/optimal_trait.png differ diff --git a/docs/src/assets/tutorials/time_varying_pop.png b/docs/src/assets/tutorials/time_varying_pop.png new file mode 100644 index 0000000000000000000000000000000000000000..b20f5234e81ee06ddc6943d823a3e7ba4556c5c1 Binary files /dev/null and b/docs/src/assets/tutorials/time_varying_pop.png differ diff --git a/docs/src/examples/changing_environment.md b/docs/src/examples/changing_environment.md new file mode 100644 index 0000000000000000000000000000000000000000..487b4f9192f84cc68cddcd0b7411ef2d27a25d2a --- /dev/null +++ b/docs/src/examples/changing_environment.md @@ -0,0 +1,36 @@ +# Changing environments + +In this tutorial, we define a birth function that is time dependent. This can be related to changing environment, where the optimal adaptive trait changes because of underlying resource variability, e.g. related to climate. + +## Defining the variation +julia + ω = 2* π / 150 # angular frequency + optimal_trait(t) = sin(ω * t) + tend = 300 + Plots.plot(1:tend,optimal_trait,label = "Optimal trait",xlabel = "time") + +![](../assets/tutorials/optimal_trait.png) + +## Running +optimal_trait function is fed into the birth function, that we define as gaussian. + +julia + myspace = (RealSpace{1,Float64}(),) + K0 = 1000 # We will have in total 1000 individuals + b(X,t) = gaussian(X[1],optimal_trait(t),1) + d(X,Y,t) = 1/K0 + D = (5e-2,) + mu = [1.] + NMax = 2000 + p = Dict{String,Any}();@pack! p = d,b,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.) + + +## Plotting + +julia +Plots.plot(sim) + +![](../assets/tutorials/time_varying_pop.png) diff --git a/docs/src/examples/delta_competition_example.md b/docs/src/examples/delta_competition_example.md index 3bdb01732d058d9d6eab1feb3b7d66298ce518b2..1a91b1ec220ffe04ff0eed9261b1a63a04b9a812 100644 --- a/docs/src/examples/delta_competition_example.md +++ b/docs/src/examples/delta_competition_example.md @@ -96,7 +96,7 @@ Plots.plot(get_tspan(sim),wsize, xlabel ="time", grid = false)  -![](tutorials/delta_comp_wsize.png) +![](../assets/tutorials/delta_comp_wsize.png) !!! notes "Callbacks" @@ -113,4 +113,4 @@ Plots.plot(sim, grid = false, markersize = 10)  -![](tutorials/delta_comp_pos.png) +![delta_comp_pos](../assets/tutorials/delta_comp_pos.png) diff --git a/docs/src/examples/genetic_structure.md b/docs/src/examples/genetic_structure.md index 4b17a458132a08063a4600a42c2571144c3ab184..b6e18ee1026d9bcdb8548571bb4ec4394833ef25 100644 --- a/docs/src/examples/genetic_structure.md +++ b/docs/src/examples/genetic_structure.md @@ -6,32 +6,32 @@ The genotype space is inspired from the article [The architecture of an empirica ## Defining the space julia -##### Genotype space##### -dim_neutr = 1000 -magicprop = 523728 / 32896 -g = SimpleGraph{Int16}(dim_neutr,round(Int16,dim_neutr * magicprop)) -initnode = argmax(eigenvector_centrality(g)) # This is the central node the we will use to instantiate the populations -myspace = (DiscreteSegment(Int8(1),Int8(nodes)),GraphSpace(g)) # union of vector spaces + ##### Genotype space##### + dim_neutr = 1000 + magicprop = 523728 / 32896 + g = SimpleGraph{Int16}(dim_neutr,round(Int16,dim_neutr * magicprop)) + initnode = argmax(eigenvector_centrality(g)) # This is the central node the we will use to instantiate the populations + myspace = (DiscreteSegment(Int8(1),Int8(nodes)),GraphSpace(g)) # union of vector spaces  ## Defining birth, death processes and mutation julia -K0 = 1000 -sigma_K = 1.; -sigma_a = .8; -mu = [1.,1.] -b(X,t) = 1 / nodes -d(X,Y,t) = (X[1] ≈ Y[1]) / K0 -D = (5e-1,1.4825) + K0 = 1000 + mu = [1.,1.] + b(X,t) = 1 / nodes + d(X,Y,t) = (X[1] ≈ Y[1]) / K0 + D = (5e-1,1.4825)  ## Final set up julia -NMax = 2000 -# tend = 1.5 -tend = 3000 -p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu -myagents = [Agent(myspace,(rand(Int8(1):Int8(nodes)),initnode),ancestors=true,rates=true) for i in 1:round(K0/nodes)] + NMax = 2000 + # tend = 1.5 + tend = 3000 + p_default = Dict{String,Any}();@pack! p_default = d,b,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.)  + +## Plotting diff --git a/docs/src/examples/sympatric_speciation.md b/docs/src/examples/sympatric_speciation.md index 9cb28f8aea7277bfeadd8d75b6c1fd1c89a438fe..27f4ca1c0e22d5b02a3e4e225b7af5f5476d3bb6 100644 --- a/docs/src/examples/sympatric_speciation.md +++ b/docs/src/examples/sympatric_speciation.md @@ -1,3 +1,20 @@ # Modelling Sympatric speciation This script aims at reproducing the 1999 article of Doebeli [On The Origin of Species By Sympatric Speciation](http://www.nature.com/articles/22521). + +In this article, birth and death functions are defined as gaussian, with respective variance \sigma_b and \sigma_d. It is shown that when \sigma_d < \sigma_b, speciation in the trait space occurs. This is what we reproduce here. + +## Running the world +![]() + +## Plotting lineages +A cool feature of ABMEv.jl is its ability to track agents ancestors traits (cf [Agent section](../manual/agent.md)) + +On can plot it, to get an idea of the coalescence of the population. + +![]() + +Beautiful, isn't it? + +!!! tip "Making sense of trait histories" + Some metrics are available (cf [Metrics section](../manual/metrics.md)) that summarize the divergence in trait value (or geographical position) through time). diff --git a/docs/src/index.md b/docs/src/index.md index 0da73013631b3dd52f0e8882e0bc1025ede600a6..4731a5ed9a6817c87e2be76ac2f2fe09954e6633 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -17,7 +17,17 @@ using ABMEv ## Tutorial We strongly advise to have a look at the tutorial section. - +@contents +Pages = [ + "examples/delta_competition_example.md", + "examples/changing_environment.md", + "examples/sympatric_speciation.md", + "examples/gradient_establishment.md", + "examples/genetic_structure.md", + "" + ] +Depth = 2 + ## How it works There general workflow to launch any simulation is the following diff --git a/docs/src/manual/run_world.md b/docs/src/manual/run_world.md index a963792d3446472ccceb0a6185968131f73d777d..55c70b5077c529daa9027a73e1de0eb2a68ccb1a 100644 --- a/docs/src/manual/run_world.md +++ b/docs/src/manual/run_world.md @@ -1,5 +1,8 @@ # Run the World +- NMax Maximum number of individuals that can be attained. If attained, then the programm stops. + + For now three algorithms @docs Gillepsie diff --git a/docs/src/manual/simulation.md b/docs/src/manual/simulation.md index 414114f7a8945ee844c883d51be34dccbc902d38..96f35b3e38974f26e3692ede29ea28abb0923122 100644 --- a/docs/src/manual/simulation.md +++ b/docs/src/manual/simulation.md @@ -1,11 +1,13 @@ # Simulation -A Simulation object is return by the function run! +A Simulation object is returned by the function run!. It is a container for snapshots of the world at every dt_saving time steps. It renders post processing easier, through dedicated methods to obtain time series of quantities. + + +!!! note "Default behaviour" + If df_saving is not provided, initial and last time steps will be saved. + + -- dt_saving = 10. -will allow to save the world every 10. time steps. If not specified, the algorithm will return first and last time step world. -- NMax Maximum number of individuals that can be attained. If attained, then the programm stops. -# @autodocs Modules = [ABMEv] Pages = ["ABMEv_Sim.jl"] diff --git a/examples/genetic_structure.jl b/examples/genetic_structure.jl new file mode 100644 index 0000000000000000000000000000000000000000..de3b3547987de707c63f31e78fbcadb21f6432f0 --- /dev/null +++ b/examples/genetic_structure.jl @@ -0,0 +1,27 @@ +using ABMEv,UnPack + +##### Genotype space##### +dim_neutr = 1000 +magicprop = 523728 / 32896 +g = SimpleGraph{Int16}(dim_neutr,round(Int16,dim_neutr * magicprop)) +initnode = argmax(eigenvector_centrality(g)) # This is the central node the we will use to instantiate the populations +myspace = (DiscreteSegment(Int8(1),Int8(nodes)),GraphSpace(g)) # union of vector spaces + +K0 = 1000 +mu = [1.,1.] +b(X,t) = 1 / nodes +d(X,Y,t) = (X[1] ≈ Y[1]) / K0 +D = (5e-1,1.4825) + +NMax = 2000 +# tend = 1.5 +tend = 3000 +p_default = Dict{String,Any}();@pack! p_default = d,b,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.) + +@time sim = run!(w0,Gillepsie(),tend,dt_saving=3.) + +using GraphPlot +# This is to plot with consistency +locs_x, locs_y = spring_layout(g) diff --git a/examples/sympatric_speciation.jl b/examples/sympatric_speciation.jl new file mode 100644 index 0000000000000000000000000000000000000000..560c801cd80fa45b1d6ba781559ba39e1ac927df --- /dev/null +++ b/examples/sympatric_speciation.jl @@ -0,0 +1,35 @@ +using ABMEv,UnPack,Plots + +myspace = (RealSpace{1,Float64}(),) +σ_b = .9; +σ_d = .7; +b(X,t) = 1. +d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b) +D = (1e-2,) +mu = [.1] +NMax = 2000 +tend = 1000 +p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax +myagents = [Agent(myspace,(1e-2 * randn(),),ancestors=true,rates=true) for i in 1:K0] +w0 = World(myagents,myspace,p,0.) +@time sim = run!(w0,Gillepsie(),tend,dt_saving = 4) + +Plots.plot(sim, ylabel = "Adaptive trait") +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")) + +using JLD2 +@save joinpath(@__DIR__,"sim_sympatric_speciation.jld2") sim diff --git a/examples/time_varying.jl b/examples/time_varying.jl new file mode 100644 index 0000000000000000000000000000000000000000..fab485de837963e112d45a1376a7f78270673577 --- /dev/null +++ b/examples/time_varying.jl @@ -0,0 +1,25 @@ +using ABMEv,UnPack + + +ω = 2* π / 150 # angular frequency +optimal_trait(t) = sin(ω * t) +tend = 300 +Plots.plot(1:tend,optimal_trait,label = "Optimal trait",xlabel = "time") +# savefig(joinpath(@__DIR__, "optimal_trait.png")) + +myspace = (RealSpace{1,Float64}(),) +sigma_K = 1.; +K0 = 1000 # We will have in total 1000 individuals +b(X,t) = gaussian(X[1],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,b,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.) + +using Plots +Plots.plot(sim, ylabel = "Adaptive trait") +savefig(joinpath(@__DIR__, "time_varying_pop.png"))