updated doc version

.gitignore

 docs/build
 docs/*.png
 docs/uploads
+*.jld2

docs/src/examples/delta_competition_example.md

@@ -63,7 +63,8 @@ mu = [1.]
 D = (1.5,)
 ```
 ## Running the world
-We initialise the world with initial population of size ``K_0 / 9`` located on patch `5`. We keep track of individuals' ancestors by setting `ancestors=true`. Because we wish to use `Gillepsie` algorithm, we need `rates=true` as agents' internal birth and death rates are updated at every time step.
+We initialise the world with initial population of size ``K_0 / 9`` located on patch `5`. `NMax` corresponds to the maximum number of individuals that can be attained. If attained, then the programm stops.
+We keep track of individuals' ancestors by setting `ancestors=true`. Because we wish to use `Gillepsie` algorithm, we need `rates=true` as agents' internal birth and death rates are updated at every time step.
 !!! note "Warning"
     `rates` treatment is something we might implement in the library internals.
 

docs/src/examples/genetic_structure.md

-# Modelling agents with a genetic structure
+# [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.
 
@@ -35,3 +35,9 @@ The genotype space is inspired from the article [The architecture of an empirica
 ```
 
 ## Plotting
+Plotting animations of graphs is a bit involved. If you want to learn about that, consult  `src/examples/genetic_structure.jl`.
+### Geographical space
+![](../assets/tutorials/space_genetic_struct.gif)
+
+### Genetic space
+![](../assets/tutorials/animated_gen_space.gif)

docs/src/examples/gradient.md

@@ -42,7 +42,7 @@ Here we show how to plot a cool animated scatter plot of the trait space through
 ```
 ![](../assets/tutorials/gradient_2Dtrait.gif)
 
-## Plotting lineages
+## [Plotting lineages](@id 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.

docs/src/examples/sympatric_speciation.md

@@ -5,16 +5,22 @@ This script aims at reproducing the 1999 article of Doebeli [On The Origin of Sp
 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).
+We need to run the simulation for a significant amount of time to observe sympatric speciation. As such, we set `ancestors=false`. The rest is pretty standard
+```julia
+myspace = (RealSpace{1,Float64}(),)
+σ_b = .9;
+σ_d = .7;
+K0 = 1000
+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 = 1500
+p = Dict{String,Any}();@pack! p = d,b,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)
+Plots.plot(sim, ylabel = "Adaptive trait")
+```
+![](../assets/tutorials/sympatric_speciation.png)

docs/src/index.md

@@ -7,9 +7,13 @@ The purpose of this package is to provide a numerical laboratory for evolutionar
 - analysis tools to investigate the simulations.
 
 ## Features
-Agents consist of a set of traits in some combination of vector spaces. A vector space can represent for example a geographical landscape, or trait space. Spaces can be of any dimensions, discrete or continuous, bounded or unbounded. They can equally consist of graphs.
+Agents consist of sets of traits in some combination of vector spaces. A vector space can represent for example **a geographical landscape, a trait space, or genetic structure**. Spaces can be of any dimensions, discrete or continuous, bounded or unbounded. They can equally consist of graphs.
 Vector spaces are used to define birth and death processes, as well as mutation processes.
 
+### Specificities
+- [ABMEv.jl allows to keep track of agents' trait lineages](@ref genetic_structure)
+- [ABMEv.jl enables to run evolutionary dynamics on graphs!](@ref lineages)
+
 ## Getting started
 ```@repl
 using ABMEv
@@ -22,14 +26,14 @@ Pages = [
     "examples/delta_competition_example.md",
     "examples/changing_environment.md",
     "examples/sympatric_speciation.md",
-    "examples/gradient_establishment.md",
+    "examples/gradient.md",
     "examples/genetic_structure.md",
     ""
     ]
 Depth = 2
 ```
 ## How it works
-There general workflow to launch any simulation is the following
+General workflow to launch any simulation is the following
 
 - [Define the combination of vector spaces you are interested in.](manual/space.md)
 - Define birth and death function, that depend on agents position in the space
@@ -40,9 +44,9 @@ There general workflow to launch any simulation is the following
 
 ### Available algorithms
 As of now, three types of simulation algorithm can be used:
-        - [Gillepsie](manual/gillepsie.md)
-        - [Wright-Fisher](manual/wright_fisher.md)
-        - [CFM](CFM.md)
+- [Gillepsie](manual/gillepsie.md)
+- [Wright-Fisher](manual/wright_fisher.md)
+- [CFM](CFM.md)
 
 ## References
 - [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632)

docs/src/manual/CFM.md

@@ -6,5 +6,30 @@ Indeed, at every time step, only the fitness of the individual picked at random
 
 In order to use it, you need to feed to the dictionary parameters `p` a constant `Cbar<:Real` that is the upperbound of the maximum of the sum of the birth and death rates (cf article).
 
+## An example on how to use it
+```julia
+using ABMEv,UnPack,Plots
+myspace = (RealSpace{1,Float64}(),)
+σ_b = .9;
+σ_d = .7;
+K0 = 1000
+b(X,t) = 1.
+d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b)
+Cbar = b([0],0.) + d([0],[0],0.)
+D = (1e-2,)
+mu = [.1]
+NMax = 2000
+tend = 1500
+p = Dict{String,Any}();@pack! p = d,b,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)
+```
+
 !!! warning "Development"
     CFM gives an approximate time step. As of now, we do not manage to obtain qualitatively the same results as the Gillepsie algorithm.
+
+```@autodocs
+Modules = [ABMEv]
+Pages   = ["algo/ABMEv_CFM.jl"]
+```

docs/src/manual/agent.md

 # Agent properties
 
 ## The `Agent` structure
-The atomic structure is the `Agent`. It has four attributes
+`Agent` is the atomic structure of ABMEv.jl. It has four attributes
 -  the ancestors' history of traits, and the corresponding time where the traits have changed,
 - a death rate and a birth rate.
 ```julia
@@ -16,15 +16,12 @@ mutable struct Agent{A<:Ancestors,R<:Rates,T<:Tuple,U,V} <: AbstractAgent{A,R}
     b::V
 end
 ```
-!!! note Specificities
+!!! note "Specificities"
     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
 
-
-#
-
 ```@autodocs
 Modules = [ABMEv]
 Pages   = ["ABMEv_Agent.jl"]

docs/src/manual/plots.md

 # Plotting
 
-ABMEv comes with Plot recipes:
-`plot(world::Array{U},p;what=["x","H"],trait = 1) where U <: Union{Missing,Agent{T}} where T`.
-
-An example to use it:
+ABMEv comes with Plot recipes.
 ```julia
-using Plots;pyplot()
-Plots.plot(worldall,p_default,what=["x"],trait=2)
+function plot(sim::Simulation;trait = 1)
 ```
-You can specify what you want to plot in the array ```what```:
-- ```"x"``` returns a scatter plot `(xaxis = time, yaxis = trait)` . Geotrait corresponds to `trait=0`
-- `"xs"` only works for agent type `MixedAgent`, because it needs a discrete geographical space. It returns a scatter plot `(xaxis = geographical component, yaxis = trait value)`. Very similar to a `histogram2d` plot, with nicer look.
-- ```"gs"``` returns a scatter plot `(xaxis = geotrait, yaxis = trait value)`
-- ```"3dgeo"``` plots a 3d diagram with x axis as geotrait and y axis as the second component. :warning: this is probably weekend
-- ```"3d"``` plots a 3d diagram with first and second component as x and y axis
-- ```"var"``` plots the variance of the  component specified by ```trait=2``` :question: with respect to time?
-- ```"vargeo"``` plots the variance of the geotrait
 
-```@autodocs
-Modules = [ABMEv]
-Pages   = ["ABMEv_plotting.jl"]
-```
+## Arguments
+- if `length(trait) == 1` then we scatter plot `trait` along time
+- if `2 <= length(trait) <= 3` then we project world of the last time step in the two  (three) dimensional trait space define by `trait`
+
+!!! warning "To be implemented"
+        We might want to get a 3 dimensional scatter plot
+        with time, trait1 and trait2 as axis

docs/src/manual/run_world.md

 # Run the World
 
-- ```NMax``` Maximum number of individuals that can be attained. If attained, then the programm stops.
-
+```@docs
+run!
+```
 
 For now three algorithms
 ```@docs
@@ -10,7 +11,3 @@ CFM
 WF
 ```
 > Warning : `WF` not implemented yet
-```@autodocs
-Modules = [ABMEv]
-Pages   = ["ABMEv_runworld.jl"]
-```

docs/src/manual/space.md

 # Space
 
-# Utils
 ```@autodocs
 Modules = [ABMEv]
 Pages   = ["ABMEv_Space.jl"]

examples/cfm.jl

@@ -4,13 +4,13 @@ myspace = (RealSpace{1,Float16}(),)
 sigma_K = .9;
 sigma_a = .7;
 K0 = 5000
-b(X) = gaussian(X[1],0.,sigma_K)
-d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+b(X,t) = gaussian(X[1],0.,sigma_K)
+d(X,Y,t) = gaussian(X[1],Y[1],sigma_a)/K0
 D = (Float16(1e-1),)
 mu = [1.]
 NMax = 7000
 tend = 50
-Cbar= b([0]) + d([0],[0])
+Cbar= b([0],0.) + d([0],[0],0.)
 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax,Cbar
 # myagents = [Agent(myspace,(- Float16(0.5) .+ .01 .* randn(Float16),)) for i in 1:K0]
 myagents = [Agent(myspace,(0.,)) for i in 1:K0]

examples/genetic_structure.jl

-using ABMEv,UnPack
+using ABMEv,LightGraphs,UnPack
 
 ##### Genotype space#####
+nodes = 9
 dim_neutr = 1000
 magicprop = 523728 / 32896
 g = SimpleGraph{Int16}(dim_neutr,round(Int16,dim_neutr * magicprop))
@@ -11,20 +12,20 @@ K0 = 1000
 mu = [1.,1.]
 b(X,t) = 1 / nodes
 d(X,Y,t) = (X[1] ≈ Y[1]) / K0
-D = (5e-1,1.4825)
+D = (3e-1,5e-1)
 
 NMax = 2000
 # tend = 1.5
-tend = 100
+tend = 200
 p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu,D
-myagents = [Agent(myspace,(rand(Int8(1):Int8(nodes)),initnode),ancestors=true,rates=true) for i in 1:round(K0/nodes)]
+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.)
 
 using GraphPlot,StatsBase
 
-# Plotting genetic space
+## Plotting genetic space
 # This is to plot with consistency
 locs_x, locs_y = spring_layout(g)
 xarray,tspan = get_xnt(sim, trait=2)
@@ -37,11 +38,49 @@ for x in xarray
         push!(d_i,_d_i)
 end
 # Here we normalize with respect to the max size of the whole world through time
-d_i = [ (_d_i .- minimum(_d_i)) ./ (maximum(_d_i) .- minimum(_d_i)) for _d_i in d_i ]anim = @animate for t in 1:size(worldall,2)
-Plots.plot(worldall,p_default,what=["gs"],trait=2,tplot=t,
-title = @sprintf("               t = %1.2f",p_default["tspan"][t]),
-ylims = (1,p_default["nodes"]),
-xlims = (0,maximum(xgeo))
-)
-                end
-gif(anim,"gif_d2=$(p_default["D"][2])_d1=$(p_default["D"][1])_tend=$(p_default["tend"])_mu=1e0_1e0.gif",fps = 13)
+d_i = [ (_d_i .- minimum(_d_i)) ./ (maximum(_d_i) .- minimum(_d_i)) for _d_i in d_i ]
+using Printf,Cairo,Compose
+for i in 1:length(d_i)
+        gp =gplot(g,locs_x,locs_y,
+                nodefillc=eth_grad_std[d_i][i])
+        draw(PNG(joinpath(@__DIR__,"gen_struct/gif_genetic_structure_$i.png"), 16cm, 16cm), gp)
+
+        # Plots.plot(xarray,p_default,what=["gs"],trait=2,tplot=t,
+        # title = @sprintf("               t = %1.2f",p_default["tspan"][t]),
+        # ylims = (1,p_default["nodes"]),
+        # xlims = (0,maximum(xgeo))
+        # )
+end
+## Plotting spatial space
+using Plots
+# This is to plot with consistency
+xarray,tspan = get_xnt(sim, trait=1)
+# here we compute the number of individuals per genome, throught all the time steps
+d_i = []
+for x in xarray
+        _d_i = zeros(nv(g))
+        c = countmap(x)
+        _d_i[collect(keys(c))] .= values(c)
+        push!(d_i,_d_i)
+end
+# Here we normalize with respect to the max size of the whole world through time
+d_i = [ (_d_i .- minimum(_d_i)) ./ (maximum(_d_i) .- minimum(_d_i)) for _d_i in d_i ]
+anim  = @animate for i in 1:length(d_i)
+        Plots.scatter(collect(1:nodes),ones(nodes),
+                markercolor=ABMEv.eth_grad_small[d_i][i],
+                markersize = 40,
+                grid = false,
+                xaxis = false,
+                yaxis = false,
+                label = "",
+                markerstrokewidth = 0.,
+                ylims = (-8,10),
+                xlims = (0,10)
+                )
+        # Plots.plot(xarray,p_default,what=["gs"],trait=2,tplot=t,
+        # title = @sprintf("               t = %1.2f",p_default["tspan"][t]),
+        # ylims = (1,p_default["nodes"]),
+        # xlims = (0,maximum(xgeo))
+        # )
+end
+gif(anim,joinpath(@__DIR__, "space_genetic_struct.gif"),fps = 13)

examples/gradient.jl

+using Revise,ABMEv,Plots,UnPack
+
+
+nodes = 9
+dim_neutr = 30
+geospace = DiscreteSegment(Int8(1),Int8(nodes))
+adaptivespace = RealSpace{1,Float16}()
+myspace = (geospace,adaptivespace)
+sigma_K = 1.;
+sigma_a = .8;
+K0 = 1000;
+mu = [1.,1.]
+a = 1
+b(X,t) = gaussian(X[2],X[1] * a,sigma_K) / nodes
+d(X,Y,t) =  (X[1] ≈ Y[1]) * gaussian(X[2],Y[2],sigma_a) / K0
+NMax = 2000
+D = (5e-1,5e-2)
+# tend = 1.5
+tend = 1500
+p = Dict{String,Any}();@pack! p = d,b,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);
+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)
+    Plots.plot(s, ylabel = "Adaptive trait",
+        trait=(1,2),
+        time=i,
+        xaxis = "Position",
+        yaxis = "Adaptive trait",
+        title = @sprintf("               t = %1.2f",s.tspan[i]),
+        ylims = (0,11),
+        xlims = (1,11))
+end
+gif(anim,joinpath(@__DIR__, "gradient_2Dtrait.gif"),fps = 13)
+
+world = get_world(s,get_size(s))
+shistall = get_xhist.(world[:],2)
+thist = get_thist.(world[:])
+aplot = Plots.plot(thist,shistall,
+            linecolor = eth_grad_std[0.],
+            label = "",
+            # title = latexstring("\\sigma_\\mu=",@sprintf("%1.2f",world.p["D"][2]),", \\sigma_D=",@sprintf("%1.2f",world.p["D"][1])),
+            grid = false,
+            xlabel = "time",
+            ylabel = "Lineages adaptive trait",
+            ylims = (-0.5,10.5)
+            )
+savefig(aplot,joinpath(@__DIR__, "gradient_lineages_adaptive_trait.png"))
+
+# Plotting lineages

examples/sympatric_speciation_CFM.jl

+using Revise
+using ABMEv,UnPack,Plots
+
+myspace = (RealSpace{1,Float64}(),)
+σ_b = .9;
+σ_d = .7;
+K0 = 1000
+b(X,t) = 1.
+d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b)
+Cbar = b([0],0.) + d([0],[0],0.)
+D = (1e-2,)
+mu = [.1]
+NMax = 2000
+tend = 1500
+p = Dict{String,Any}();@pack! p = d,b,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)
+
+using JLD2
+@save joinpath(@__DIR__,"sim_sympatric_speciation_CFM.jld2") sim
+
+Plots.plot(sim, ylabel = "Adaptive trait")
+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"))

src/ABMEv_runworld.jl

 """
-$(SIGNATURES)
+    function run!(w::World{A,S,T},alg::L,tend::Number;dt_saving=nothing,cb=(names = String[],agg =nothing))
+
 Run `w` with algorithm `alg`, until `tend` is reached.
 Returns a `Simulation` type.
 - `worldall` stores the world every `p["dt_saving"]` time steps.
 If `dt_saving` not specified, `sim` contains an array of two elements,
 first corresponding to initial conditions and last corresponding to world in the last time step.
+- the run stops if the number of agents reaches`p["NMax"]`.
 >:warning: if you choose `nagents = 1` then nothing is saved until the end of the simulation.
 """
-# function run(w::World{AbstractAgent{A,R},S,T},g::G;dt_saving=nothing,callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
 function run!(w::World{A,S,T},alg::L,tend::Number;
                 dt_saving=nothing,
                 cb=(names = String[],agg =nothing)) where {A,S,T,L<:AbstractAlg}
@@ -46,12 +47,11 @@ function run!(w::World{A,S,T},alg::L,tend::Number;
 end
 
 """
-function _correct_timedep!(p::Dict)
+    function _correct_timedep!(p::Dict)
 
 checks time dependency of birth and death functions,
 and overloads the function if not provided
 """
-
 function _check_timedep(p::Dict)
     @unpack d,b = p
     if numargs(p["b"]) < 2