Commit bcbc425f by Victor

update docs tuts

Pipeline #75086 passed with stage
in 24 minutes and 10 seconds

143 KB

19.7 KB

 ... @@ -6,7 +6,7 @@ In this script, we model the evolution of a population where agents are simply d ... @@ -6,7 +6,7 @@ In this script, we model the evolution of a population where agents are simply d Let's start by a linear landscape. We define a discrete segment of length 9, with reflecting boundary conditions. In fact, reflecting boundary conditions are implemented for any finite space. Let's start by a linear landscape. We define a discrete segment of length 9, with reflecting boundary conditions. In fact, reflecting boundary conditions are implemented for any finite space. !!! warning "1D Reflections" !!! warning "1D Reflections" As of v1, only reflections on one dimensional space are implemented. We have a two dimensional reflection method that will be released in the future. As of v1, only reflections on one dimensional space are implemented. We have a two dimensional reflection method that will be released in the future. There are two ways of implementing a linear landscape. The first one uses a DiscreteSegment while the second relies on LightGraphs.jl library. Both are *almost* equivalent. There are two ways of implementing a linear landscape. The first one uses a DiscreteSegment while the second relies on LightGraphs.jl library. Both are *almost* equivalent. ... @@ -15,6 +15,7 @@ There are two ways of implementing a linear landscape. The first one uses a Dis ... @@ -15,6 +15,7 @@ There are two ways of implementing a linear landscape. The first one uses a Dis using ABMEv using ABMEv nodes = 10 nodes = 10 mysegment = DiscreteSegment(1,nodes) mysegment = DiscreteSegment(1,nodes) wholespace = (mysegment,)   ### grid ### grid ... @@ -23,12 +24,15 @@ using ABMEv, LightGraphs ... @@ -23,12 +24,15 @@ using ABMEv, LightGraphs nodes = 10 nodes = 10 g = grid([nodes,1]) g = grid([nodes,1]) mysegmentgraph = GraphSpace(g) mysegmentgraph = GraphSpace(g) wholespace = (mysegmentgraph,)   !!! warning "Space tuple" Notice that the whole space should be a *tuple of spaces*. Even where there is only one sub vector space as here, you need to have brackets and comma around the unit vector space. Here is how you can visualise the landscape. Here is how you can visualise the landscape. julia julia using GraphPlots using GraphPlot gplot(g,locs_x = 1:nodes,locs_y=1:nodes) gplot(g, collect(1:nodes), collect(1:nodes))   ## Defining competition processes ## Defining competition processes ... @@ -49,6 +53,9 @@ d(X,Y,t) = (X[1] ≈ Y[1]) / K0 ... @@ -49,6 +53,9 @@ d(X,Y,t) = (X[1] ≈ Y[1]) / K0   At equilibrium, population size in each deme will reach K0 / nodes. At equilibrium, population size in each deme will reach K0 / nodes. !!! warning "Time dependency" Even though time is not used, you have to write birth and death functions with time dependency. ## Dispersal ## Dispersal We assume that anytime an offspring is born, it is given a chance to move (\mu = 1). We assume that anytime an offspring is born, it is given a chance to move (\mu = 1). julia julia ... @@ -58,31 +65,52 @@ D = (1.5,) ... @@ -58,31 +65,52 @@ D = (1.5,) ## Running the world ## 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. 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" !!! note "Warning" rates treatment is something we might implement in the library internals. rates treatment is something we might implement in the library internals. julia julia using UnPack# useful macro @pack! NMax = 2000 NMax = 2000 tend = 2 tend = 300. p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax myagents = [Agent(myspace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] myagents = [Agent(myspace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,myspace,p,0.) w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,Gillepsie(),tend) @time sim = run!(w0,Gillepsie(),tend)   This is the simplest run you can do. Now time to more interesting things ## Analysis ## Analysis ### Size of the world Let's verify that the population's growth is logistic. We will plot the population size over time. Let's verify that the population's growth is logistic. We will plot the population size over time. To do so, one need to define dt_saving. To do so, one need to define dt_saving < tend to save every dt_saving time steps of the world. julia julia 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.) wsize = [length(w) for w in sim[:]] using Plots using Plots @time sim = run!(w0,Gillepsie(),tend,dt_saving=1) Plots.plot(get_tspan(sim),wsize, Plots.plot() label = "", @animate for (i,t) in tspan(sim) ylabel = "Metapopulation size", Plots. xlabel ="time", end grid = false)   ![](tutorials/delta_comp_wsize.png) One might be tempted to plot the world for some time steps. This requires some knowledge of the library Plots.jl. !!! notes "Callbacks" Agents values are accessed by get_x(w::World). Note that one could also use a callback function to obtain time series of size of the world computed at simulation time. See [Simulation page](../manual/simulation.md). !!! notes "Plotting" ### Position through time This will be proposed as a blackbox with PlotRecipes.jl in the near future One might be tempted to plot the agents position for some time steps. julia Plots.plot(sim, label = "", ylabel = "Geographical position", grid = false, markersize = 10)  ![](tutorials/delta_comp_pos.png)
 using Revise,ABMEv nodes = 10 mysegment = DiscreteSegment(1,nodes) using ABMEv, LightGraphs nodes = 10 g = grid([nodes,1]) mysegmentgraph = GraphSpace(g) wholespace = (mysegmentgraph,) using GraphPlot gplot(g, collect(1:nodes), collect(1:nodes)) K0 = 1000 # We will have in total 1000 individuals b(X,t) = 1 / nodes d(X,Y,t) = (X[1] ≈ Y[1]) / K0 mu = [1.] D = (1.5,) using UnPack NMax = 2000 tend = 300. p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax 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) ### 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.) wsize = [length(w) for w in sim[:]] using Plots Plots.plot(get_tspan(sim),wsize, label = "", ylabel = "Metapopulation size", xlabel ="time", grid = false) savefig(joinpath(@__DIR__, "delta_comp_wsize.png")) ## plotting position through time using Plots Plots.plot(sim, label = "", ylabel = "Geographical position", grid = false, markersize = 10) savefig(joinpath(@__DIR__, "delta_comp_pos.png"))
 ... @@ -54,7 +54,7 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F} ... @@ -54,7 +54,7 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F} end end #TODO : code it #TODO : code it function get_xnt(s::Simulation;trait = i) function get_xnt(s::Simulation;trait = 1) return [getindex.(wa,trait) for wa in s.agentarray],[fill(t,size(s[j])) for (j,t) in enumerate(s.tspan)] return [getindex.(wa,trait) for wa in s.agentarray],[fill(t,size(s[j])) for (j,t) in enumerate(s.tspan)] end end # get_x(agentarray::Array{T},t,trait::Integer) where {T <: AbstractAgent} = reshape(hcat(get_x.(agentarray,t,trait)),size(agentarray,1),size(agentarray,2)) # get_x(agentarray::Array{T},t,trait::Integer) where {T <: AbstractAgent} = reshape(hcat(get_x.(agentarray,t,trait)),size(agentarray,1),size(agentarray,2)) ... ...
 ... @@ -2,101 +2,48 @@ using RecipesBase ... @@ -2,101 +2,48 @@ using RecipesBase using Colors using Colors import KernelDensity:kde,pdf import KernelDensity:kde,pdf """ """ function plot(world::Array{U},p;what=["x","H"],trait = 1,tplot = false) where U <: Union{Missing,Agent} function plot(sim::Simulation;trait = 1) Plot recipe for ABMEv.jl # ARGS # ARGS - what = ["x","H"]: the plots you want to obtain - if length(trait) == 1 then we scatter plot trait along time - trait = 1: the trait that will plotted regarding what you asked. trait = 0 will plot the geotrait - if 2 <= length(trait) <= 3 then we project world of the - tplot = 0 used when calling xs, as it plots a snapshot of the world at a particular time last time step in the two (three) dimensional trait space define by trait It should correspond to an integer, as it indexes the column to plot # Options available !!! warning "To be implemented" - "x" We might want to get a 3 dimensional scatter plot - "xs" with time, trait1 and trait2 as axis """ """ @recipe function plot(sim::Simulation;what=["x","H"],trait = 1,tplot = 0) @recipe function plot(sim::Simulation;trait = 1,time = nothing) world = sim.agentarray # world = sim.agentarray tot_dim = size(world,2)*size(world,1) # tspan = sim.tspan tspan = sim.tspan # p = sim.p p = sim.p if length(trait) == 1 # We reduce time interval if it is too big xarray,tspan = get_xnt(sim, trait=trait) # if tot_dim > 1e6 && size(world,2) >= 200 # p = copy(p) # idx_reduced = floor.(Int,range(1,size(world,2),length = 200)) # p["tspan" ] = tspan[idx_reduced] # world = world[:,idx_reduced] # end # second condition is to make sure that the world corresponds to all the time steps # If not, then we can not get "x" if count(ismissing,world) > 0 && length(tspan) == size(world,2) tspan_ar = vcat([tspan[i]*ones(Int(p["NMax"] - count(ismissing,world[:,i]))) for i in 1:length(tspan) ]...); else tspan_ar = repeat(tspan,inner = size(world,1)) end # tspan = Float64.(tspan) tend = tspan[tplot > 0 ? tplot : length(tspan)] world_sm = clean_world(world) if "x" in what d_i = [] d_i = [] for i in 1:length(tspan) for x in xarray x = get_x.(clean_world(world[:,i]),tspan[i],trait)[:] push!(d_i,pdf(kde(x),x)) append!(d_i,pdf(kde(x),x)) end 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 ] @series begin @series begin if length(world_sm) !== length(tspan_ar[:]) throw(DimensionMismatch("You want to plot a world with missing data")) end xarray = get_x.(world_sm,tspan_ar[:],trait) seriestype := :scatter seriestype := :scatter markercolor := eth_grad_small[d_i ./ maximum(d_i)] markercolor := eth_grad_std[vcat(d_i...)] markerstrokewidth := 0 markerstrokewidth := 0 seriesalpha :=1. seriesalpha := .2 # xlabel := "time" # xlabel := "time" # ylabel := "trait value" # ylabel := "trait value" label := "" label := "" xlabel := "time" grid := false grid := false # markersize := 2.3/1000*size(world_sm,1) # markersize := 2.3/1000*size(world_sm,1) tspan_ar[:],xarray[:] vcat(tspan...),vcat(xarray...) end end end end # we use this for discrete agents if length(trait) == 2 # world should be a one dimensional vector, corresponding to one time step only isnothing(time) ? w = get_world(sim[end]) : w = get_world(sim[time]) # if "xs" in what y = get_x(w,trait[1]); # d_i = []; xt_array = []; x1_array = [] x=get_x(w,trait[2]) # world_df_all = world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true) # world_df_g = groupby(world_df_all,:x1) # for world_df in world_df_g # if trait == 0 # x = Float64.(world_df.g) # else # # fitness occupies first spot # x = world_df[:,trait+1] ; # end # x1 = world_df.x1; # append!(d_i,pdf(kde(x),x)) # append!(xt_array,x) # append!(x1_array,x1) # end # @series begin # seriestype := :scatter # markercolor := eth_grad_small[d_i ./ maximum(d_i)] # # markercolor := :blue # markerstrokewidth := 0 # # seriesalpha := 1. # xaxis := "geographical position" # xticks := sort!(unique(world_df_all.x1)) # yaxis := "trait value" # label := "" # grid := false # # marker := (:rect,20,1.) # x1_array[:],xt_array[:] # end # end if "gs" in what _world = clean_world(world[:, tplot > 0 ? tplot : length(tspan) ]) y = get_x(_world,tend,2)[:] x = get_x(_world,tend,0)[:] X = hcat(x,y) X = hcat(x,y) d = kde(X) d = kde(X) # by density # by density ... @@ -119,100 +66,9 @@ It should correspond to an integer, as it indexes the column to plot ... @@ -119,100 +66,9 @@ It should correspond to an integer, as it indexes the column to plot x,y x,y end end end end if "3dgeo" in what if length(trait) == 2 d_i = [] throw(ArgumentError("Plot for three traits not yet implemented")) for i in 1:size(world,2) _world = clean_world(world[:,i]) x = get_x(_world,tspan[i],2)[:] y = get_x(_world,tspan[i],0)[:] X = hcat(x,y) # d = kde(X) # di_temp = diag(pdf(d,X[:,1],X[:,2])) di_temp = y di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp)) # here we normalise with respect to maximum value at each time step append!(d_i,di_temp) end @series begin xarray = get_geo.(world_sm,tspan_ar) yarray = get_x(world_sm,2) seriestype := :scatter3d markercolor := eth_grad_std[d_i ./ 1.] markerstrokewidth := 0 seriesalpha :=.1 xlabel := "time" ylabel := "geotrait" zlabel := "trait value" label := "" # markersize := 2.3/1000*size(world_sm,1) tspan_ar,xarray[:],yarray[:] end end if "3d" in what d_i = [] for i in 1:size(world,2) x = get_x(clean_world(world[:,i]),tspan[i],1)[:] y = get_x(clean_world(world[:,i]),tspan[i],2)[:] X = hcat(x,y) d = kde(X) di_temp = diag(pdf(d,X[:,1],X[:,2])) di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp)) append!(d_i,di_temp) end @series begin xarray = get_x(world_sm,1) yarray = get_x(world_sm,2) seriestype := :scatter3d markercolor := eth_grad_small[d_i ./ maximum(d_i)] markerstrokewidth := 0 seriesalpha :=.1 xlabel := "time" ylabel := "position" zlabel := "trait value" label := "" # markersize := 2.3/1000*size(world_sm,1) tspan_ar,xarray[:],yarray[:] end end # if "H" in what # @series begin # x = get_x.(world_sm,trait) # linewidth := 2 # seriestype := :line # label := "Interconnectedness" # tspan,N^2 / 2 .* [H_discrete(x[:,i]) for i in tspan] # end # end if "var" in what @series begin linewidth := 2 seriestype := :line label := "Variance" xlabel := "Time" ylabel := "Variance" tspan,var(world_sm,trait=trait)[:] end end end if "vargeo" in what @series begin linewidth := 2 seriestype := :line label := "Variance of geotrait" xlabel := "Time" ylabel := "Variance" tspan,i->first(covgeo(world_sm[:,Int(i)])) end end # if "density_t" in what # @series begin # linewidth := 2 # seriestype := :plot3d # label := "Variance of geotrait" # xlabel := "Time" # ylabel := "Variance" # tspan,i->first(covgeo(world_sm[:,Int(i)])) # end # end end end ... @@ -221,3 +77,137 @@ import Plots:cgrad ... @@ -221,3 +77,137 @@ import Plots:cgrad const eth_grad_small = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,.1]) const eth_grad_small = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,.1]) # symmetry between red and blue # symmetry between red and blue const eth_grad_std = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,1.]) const eth_grad_std = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,1.]) ########## OLD PLOT RECIPES ############## # To get some inspiration # we use this for discrete agents # world should be a one dimensional vector, corresponding to one time step only # if "xs" in what # d_i = []; xt_array = []; x1_array = [] # world_df_all = world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true) # world_df_g = groupby(world_df_all,:x1) # for world_df in world_df_g # if trait == 0 # x = Float64.(world_df.g) # else # # fitness occupies first spot # x = world_df[:,trait+1] ; # end # x1 = world_df.x1; # append!(d_i,pdf(kde(x),x)) # append!(xt_array,x) # append!(x1_array,x1) # end # @series begin # seriestype := :scatter # markercolor := eth_grad_small[d_i ./ maximum(d_i)] # # markercolor := :blue # markerstrokewidth := 0 # # seriesalpha := 1. # xaxis := "geographical position" # xticks := sort!(unique(world_df_all.x1)) # yaxis := "trait value" # label := "" # grid := false # # marker := (:rect,20,1.) # x1_array[:],xt_array[:] # end # end ## # if "3dgeo" in what # d_i = [] # for i in 1:size(world,2) # _world = clean_world(world[:,i]) # x = get_x(_world,tspan[i],2)[:] # y = get_x(_world,tspan[i],0)[:] # X = hcat(x,y) # # d = kde(X) # # di_temp = diag(pdf(d,X[:,1],X[:,2])) # di_temp = y # di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp)) # # here we normalise with respect to maximum value at each time step # append!(d_i,di_temp) # end # @series begin # xarray = get_geo.(world_sm,tspan_ar) # yarray = get_x(world_sm,2) # seriestype := :scatter3d # markercolor := eth_grad_std[d_i ./ 1.] # markerstrokewidth := 0 # seriesalpha :=.1 # xlabel := "time" # ylabel := "geotrait" # zlabel := "trait value" # label := "" # # markersize := 2.3/1000*size(world_sm,1) # tspan_ar,xarray[:],yarray[:] # end # end # if "3d" in what # d_i = [] # for i in 1:size(world,2) # x = get_x(clean_world(world[:,i]),tspan[i],1)[:] # y = get_x(clean_world(world[:,i]),tspan[i],2)[:] # X = hcat(x,y) # d = kde(X) # di_temp = diag(pdf(d,X[:,1],X[:,2])) # di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp)) # append!(d_i,di_temp) # end # @series begin # xarray = get_x(world_sm,1) # yarray = get_x(world_sm,2) # seriestype := :scatter3d # markercolor := eth_grad_small[d_i ./ maximum(d_i)] # markerstrokewidth := 0 # seriesalpha :=.1 # xlabel := "time" # ylabel := "position" # zlabel := "trait value" # label := "" # # markersize := 2.3/1000*size(world_sm,1) # tspan_ar,xarray[:],yarray[:] # end # end # if "H" in what # @series begin # x = get_x.(world_sm,trait) # linewidth := 2 # seriestype := :line # label := "Interconnectedness" # tspan,N^2 / 2 .* [H_discrete(x[:,i]) for i in tspan] # end # end # if "var" in what # @series begin # linewidth := 2 # seriestype := :line # label := "Variance" # xlabel := "Time" # ylabel := "Variance" # tspan,var(world_sm,trait=trait)[:] # end # end # if "vargeo" in what # @series begin # linewidth := 2 # seriestype := :line # label := "Variance of geotrait" # xlabel := "Time" # ylabel := "Variance" # tspan,i->first(covgeo(world_sm[:,Int(i)])) # end # end # if "density_t" in what # @series begin # linewidth := 2 # seriestype := :plot3d # label := "Variance of geotrait" # xlabel := "Time" # ylabel := "Variance" # tspan,i->first(covgeo(world_sm[:,Int(i)])) # end # end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!