diff --git a/docs/src/assets/tutorials/delta_comp_pos.png b/docs/src/assets/tutorials/delta_comp_pos.png new file mode 100644 index 0000000000000000000000000000000000000000..4cd928ee3d0cba5792cc9f8166948c350359722b Binary files /dev/null and b/docs/src/assets/tutorials/delta_comp_pos.png differ diff --git a/docs/src/assets/tutorials/delta_comp_wsize.png b/docs/src/assets/tutorials/delta_comp_wsize.png new file mode 100644 index 0000000000000000000000000000000000000000..628ef495d452225b72d2df45c8ebc74997167af1 Binary files /dev/null and b/docs/src/assets/tutorials/delta_comp_wsize.png differ diff --git a/docs/src/examples/delta_competition_example.md b/docs/src/examples/delta_competition_example.md index ff55926af05c5bc776ac56cb1c349df515e7e4e2..3bdb01732d058d9d6eab1feb3b7d66298ce518b2 100644 --- a/docs/src/examples/delta_competition_example.md +++ b/docs/src/examples/delta_competition_example.md @@ -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. !!! 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. @@ -15,6 +15,7 @@ There are two ways of implementing a linear landscape. The first one uses a Dis using ABMEv nodes = 10 mysegment = DiscreteSegment(1,nodes) +wholespace = (mysegment,)  ### grid @@ -23,12 +24,15 @@ using ABMEv, LightGraphs nodes = 10 g = grid([nodes,1]) 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. julia -using GraphPlots -gplot(g,locs_x = 1:nodes,locs_y=1:nodes) +using GraphPlot +gplot(g, collect(1:nodes), collect(1:nodes))  ## Defining competition processes @@ -49,6 +53,9 @@ d(X,Y,t) = (X[1] ≈ Y[1]) / K0  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 We assume that anytime an offspring is born, it is given a chance to move (\mu = 1). julia @@ -58,31 +65,52 @@ 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. !!! 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 +using UnPack# useful macro @pack! NMax = 2000 -tend = 2 +tend = 300. 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] w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,Gillepsie(),tend)  +This is the simplest run you can do. Now time to more interesting things + ## Analysis + +### Size of the world 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 +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 -@time sim = run!(w0,Gillepsie(),tend,dt_saving=1) -Plots.plot() -@animate for (i,t) in tspan(sim) - Plots. -end +Plots.plot(get_tspan(sim),wsize, + label = "", + ylabel = "Metapopulation size", + xlabel ="time", + 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. -Agents values are accessed by get_x(w::World). +!!! notes "Callbacks" + 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" - This will be proposed as a blackbox with PlotRecipes.jl in the near future +### Position through time + +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) diff --git a/examples/delta_competition.jl b/examples/delta_competition.jl new file mode 100644 index 0000000000000000000000000000000000000000..6cf66c44f9a248056a10bfe482e64a787a4ea33c --- /dev/null +++ b/examples/delta_competition.jl @@ -0,0 +1,47 @@ +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")) diff --git a/src/ABMEv_Sim.jl b/src/ABMEv_Sim.jl index 7ae053ad272a77f845bf9faf083732f9f3566853..51684a6b30a710038192747ff1fc99304068537c 100644 --- a/src/ABMEv_Sim.jl +++ b/src/ABMEv_Sim.jl @@ -54,7 +54,7 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F} end #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)] 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)) diff --git a/src/ABMEv_plot.jl b/src/ABMEv_plot.jl index ac365e4b6d45a3b52a659ad1e0d5c810091ced98..05bd65ceab07d06a61fb5f3b4eb054e1d7221e29 100644 --- a/src/ABMEv_plot.jl +++ b/src/ABMEv_plot.jl @@ -2,101 +2,48 @@ using RecipesBase using Colors 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 -- what = ["x","H"]: the plots you want to obtain -- trait = 1: the trait that will plotted regarding what you asked. trait = 0 will plot the geotrait -- tplot = 0 used when calling xs, as it plots a snapshot of the world at a particular time -It should correspond to an integer, as it indexes the column to plot +- 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 -# Options available -- "x" -- "xs"` +!!! warning "To be implemented" + We might want to get a 3 dimensional scatter plot + with time, trait1 and trait2 as axis """ -@recipe function plot(sim::Simulation;what=["x","H"],trait = 1,tplot = 0) - world = sim.agentarray - tot_dim = size(world,2)*size(world,1) - tspan = sim.tspan - p = sim.p - # We reduce time interval if it is too big - # 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 +@recipe function plot(sim::Simulation;trait = 1,time = nothing) + # world = sim.agentarray + # tspan = sim.tspan + # p = sim.p + if length(trait) == 1 + xarray,tspan = get_xnt(sim, trait=trait) d_i = [] - for i in 1:length(tspan) - x = get_x.(clean_world(world[:,i]),tspan[i],trait)[:] - append!(d_i,pdf(kde(x),x)) + for x in xarray + push!(d_i,pdf(kde(x),x)) 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 - 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 - markercolor := eth_grad_small[d_i ./ maximum(d_i)] + markercolor := eth_grad_std[vcat(d_i...)] markerstrokewidth := 0 - seriesalpha :=1. + seriesalpha := .2 # xlabel := "time" # ylabel := "trait value" label := "" + xlabel := "time" grid := false # markersize := 2.3/1000*size(world_sm,1) - tspan_ar[:],xarray[:] + vcat(tspan...),vcat(xarray...) end end - # 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 "gs" in what - _world = clean_world(world[:, tplot > 0 ? tplot : length(tspan) ]) - y = get_x(_world,tend,2)[:] - x = get_x(_world,tend,0)[:] + if length(trait) == 2 + isnothing(time) ? w = get_world(sim[end]) : w = get_world(sim[time]) + y = get_x(w,trait[1]); + x=get_x(w,trait[2]) X = hcat(x,y) d = kde(X) # by density @@ -119,100 +66,9 @@ It should correspond to an integer, as it indexes the column to plot x,y 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 + if length(trait) == 2 + throw(ArgumentError("Plot for three traits not yet implemented")) 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 @@ -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]) # 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.]) + + +########## 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