added a function get_xnt(s::Simulation). Deleted old files for simulating gillepsie alg

docs/src/index.md

@@ -2,28 +2,22 @@
 
 This is a suite for simulating an Agent Based Model that captures the evolutionary dynamics of a population in a multidimensional space.
 
+## How it works
+This library helps you studying the evolution of a population of agents that evolve into some multidimensional space.
+
+- Define a space
+- Define birth and death function, that depend on agents traits, population state and time
+- Define mutation function
+- Initialise the world and run the simulation according to some updating algorithm
+- Obtain a summary of the population state
+
+
+
 ## Getting started
 ```julia
 using ABMEv
 ```
 
-## The `Agent` structure
-This package is an Agent Based Model, where the atomic structure is the `Agent`. 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
-mutable struct Agent{T,U}
-    # history of traits for geotraits
-    x_history::Array{U}
-    # birth time of ancestors
-    t_history::Array{U,1}
-    # death rate
-    d::Float64
-    #birth rate
-    b::Float64
-end
-```
-For more understanding of the composite types `T,U`, check the wiki: Gillepsie / Agent Type.
 ## Parameters of the simulation
 Parameters are stored in the parameter dictionary `p`
 ### General parameters

docs/src/manual/agent.md

 # Agent properties
 
+## The `Agent` structure
+The atomic structure is the `Agent`. 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
+mutable struct Agent{A<:Ancestors,R<:Rates,T<:Tuple,U,V} <: AbstractAgent{A,R}
+    # history of traits for geotraits
+    x_history::Array{T,1}
+    # birth time of ancestors
+    t_history::Array{U,1}
+    # death rate
+    d::V
+    #birth rate
+    b::V
+end
+```
+!!! 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/mathematics/pde.md

@@ -35,3 +35,15 @@ Then we get the **large population limit without mutation scaling**
 ## Mutations
 Assuming a small mutational variance $`\max \sigma^2_i << 1`$ and a mutation rate $`U`$, the mutational effects can be approximated by an elliptic operator $`\sum_{i=1}^{n} (\mu_i^2/x)\partial_{ii}`$ with $`\mu_i = \sigma_i\sqrt{U}`$
 > :warning: check that with Burger
+
+In other words (from Champagnat, Ferriere and Meleard 2006), we have
+
+```math
+\partial_t u(t,x) =  u(t,x)\big(b(t,x) - d(x,u(t,x)) \big) + \frac{1}{2} \Delta (\sigma^2 r \mu u)(t,x)
+```
+where 
+```math
+d(x,u(t,x)) = (u \ast c)(x)
+```
+and ``c`` is the competition function.
+For us, ``x \in ``
\ No newline at end of file

examples/benchmark.jl

+using ABMEv,UnPack
+
+println("--------------------------------")
+println("""TESTING GILLEPSIE with popoulation""")
+println("--------------------------------")
+
+for K0 in [10,50,100,500,1000,5000]
+    myspace = (RealSpace{1,Float64}(),)
+    sigma_K = .9;
+    sigma_a = .7;
+    b(X) = gaussian(X[1],0.,sigma_K)
+    d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+    D = (1e-2,)
+    mu = [.1]
+    NMax = 10000
+    tend = 2
+    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.)
+    @show K0
+    @time sim = run!(w0,Gillepsie(),tend)
+end
+
+println("--------------------------------")
+println("""TESTING CFM with population""")
+println("--------------------------------")
+for K0 in [10,50,100,500,1000,5000]
+    myspace = (RealSpace{1,Float64}(),)
+    sigma_K = .9;
+    sigma_a = .7;
+    b(X) = gaussian(X[1],0.,sigma_K)
+    d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+    D = (1e-2,)
+    mu = [.1]
+    NMax = 10000
+    tend = 2
+    Cbar=2
+    p = Dict{String,Any}();@pack! p = d,b,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)
+end
+
+println("--------------------------------")
+println("""TESTING CFM with population without ancestors""")
+println("--------------------------------")
+for K0 in [10,50,100,500,1000,5000]
+    myspace = (RealSpace{1,Float64}(),)
+    sigma_K = .9;
+    sigma_a = .7;
+    b(X) = gaussian(X[1],0.,sigma_K)
+    d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+    D = (1e-2,)
+    mu = [.1]
+    NMax = 10000
+    tend = 2
+    Cbar=2
+    p = Dict{String,Any}();@pack! p = d,b,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)
+end
+
+println("--------------------------------")
+println("""TESTING CFM with time""")
+println("--------------------------------")
+
+for tend in [1,10,50,100]
+    K0 = 1000
+    myspace = (RealSpace{1,Float64}(),)
+    sigma_K = .9;
+    sigma_a = .7;
+    b(X) = gaussian(X[1],0.,sigma_K)
+    d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+    D = (1e-2,)
+    mu = [.1]
+    NMax = 10000
+    Cbar=2
+    p = Dict{String,Any}();@pack! p = d,b,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)
+end
+
+println("--------------------------------")
+println("""TESTING CFM with time without Ancestors""")
+println("--------------------------------")
+
+for tend in [1,10,50,100]
+    K0 = 1000
+    myspace = (RealSpace{1,Float64}(),)
+    sigma_K = .9;
+    sigma_a = .7;
+    b(X) = gaussian(X[1],0.,sigma_K)
+    d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+    D = (1e-2,)
+    mu = [.1]
+    NMax = 10000
+    Cbar=2
+    p = Dict{String,Any}();@pack! p = d,b,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)
+end
+
+println("--------------------------------")
+println("""TESTING Gillepsie with time""")
+println("--------------------------------")
+for tend in [1,10,50,100]
+    K0 = 1000
+    myspace = (RealSpace{1,Float64}(),)
+    sigma_K = .9;
+    sigma_a = .7;
+    b(X) = gaussian(X[1],0.,sigma_K)
+    d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+    D = (1e-2,)
+    mu = [.1]
+    NMax = 10000
+    # Cbar=2
+    p = Dict{String,Any}();@pack! p = d,b,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)
+end

examples/cfm.jl

+using Revise,ABMEv,UnPack
+
+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
+D = (Float16(1e-1),)
+mu = [1.]
+NMax = 7000
+tend = 50
+Cbar= b([0]) + d([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]
+w0 = World(myagents,myspace,p,0.)
+@time mysim = run!(w0,CFM(),tend,dt_saving=1.)
+
+x,t = get_xnt(mysim,trait=1)
+using Plots
+Plots.scatter(t,x,color =:blue,labels="" )
+
+
+########### GILLEPSIE
+K0 = 1000
+myspace = (RealSpace{1,Float64}(),)
+sigma_K = .9;
+sigma_a = .7;
+b(X) = gaussian(X[1],0.,sigma_K)
+d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
+D = (1e-1,)
+mu = [1.]
+NMax = 10000
+tend = 100
+p = Dict{String,Any}();@pack! p = d,b,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,dt_saving=1.)
+x,t = get_xnt(sim,trait=1)
+using Plots
+Plots.scatter(t,x,color =:blue,labels="" )

examples/gillepsie_simple.jl → examples/scaling.jl

@@ -5,13 +5,13 @@ using ABMEv
 sigma_K = 1.; #bandwith of resource
 sigma_a = 1.2; #bandwith of competition
 K0 = 1000 #carrying capacity
-K(X) = gaussian(X[1],0,Float32(sigma_K)) #birth
-alpha(X,Y) = gaussian(X[1],Y[1],Float32(sigma_a)) / Float32(K0) #competition
+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 = 100.
+tend = 2.
 using UnPack
 p = Dict{String,Any}()
 @pack! p = K,alpha,D,mu,NMax,dt_saving,tend

src/ABMEv.jl

@@ -31,7 +31,7 @@ module ABMEv
     export World,parameters,time,space,agents,size,maxsize,addAgent!,removeAgent!
     export run!,give_birth,updateWorld!,update_clock!,updateBirthEvent!,
         updateDeathEvent!#,runWorld_G!,runWorld_WF!,
-    export Simulation,add_entry!,get_tend,get_size,get_tspan,get_world
+    export Simulation,add_entry!,get_tend,get_size,get_tspan,get_world,get_xnt
     export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div,
         get_dist_hist,get_pairwise_average_isolation,
         get_local_pairwise_average_isolation,

src/ABMEv_Sim.jl

@@ -54,19 +54,8 @@ function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F}
 end
 
 #TODO : code it
-function get_x(s::Simulation,i)
-    i = get_size(s)
-    j = size(s.agentarray,2)
-    if i == j
-        # we double the siwe of agent array
-        s.agentarray = hcat(s.agentarray,Array{Missing}(missing,maxsize(w),j))
-    end
-    s.agentarray[1:size(w),i+1] .= copy.(agents(w))
-    push!(s.tspan,w.t)
-    if !(F==Nothing)
-        push!(s.df_agg,Dict(s.cb.names .=> [f(w) for f in s.cb.agg]))
-    end
-    return nothing
+function get_xnt(s::Simulation;trait = i)
+    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))
 # @deprecate get_x(agentarray::Array{T},t::Number,trait::Integer)