Commit 356fd36a by Victor

### updated files

parent 51bb3591
Pipeline #101414 canceled with stage
 ... ... @@ -6,16 +6,11 @@ makedocs(sitename="EvoId.jl", authors = "Victor Boussange", pages = [ "Home" => "index.md", "Manual" => ["manual/agent.md", "manual/space.md", "manual/world.md", "manual/run_world.md"], "Examples" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"examples"),join=true))], "Manual" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"manual"),join=true))], # "Mathematics" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"mathematics"),join=true))], "Library" => Any[ "Public" => "lib/public.md", # "Internals" => map( # s -> "lib/internals/$(s)", # sort(readdir(joinpath(@__DIR__, "src/lib/internals"))) # ), ], "Developping" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"dev"),join=true))], # "contributing.md", ],) ... ...  # CFM algorithm Champagnat Ferriere Méléard algorithm described in [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632). This algorithm similar to Gillepsie algorithms, excepts that it runs much faster! Indeed, at every time step, only the fitness of the individual picked at random is evaluated. Thus this algorithm is of order K times more efficient. 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 EvoId,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,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, b, d, 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 = [EvoId] Pages = ["algo/CFM.jl"]   # Gillsepie algorithm ## Mathematical foundations - The original article by Gillepsie: [**A general method for numerically simulating the stochastic time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub) ### Update Rates  b_i  and  d_i  represent respcetively birth and death rates of agents i. The total rate is given by the sum of all individual rates math R(t) = \left[ \sum_i b_i(t) + d_i(t) \right]  A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate R. @autodocs Modules = [EvoId] Pages = ["algo/Gillepsie.jl"]   # Wright Fisher algorithm ## Foundations The Wright Fisher process is an individual based model where the number of agents is constant through time. It is helpful to visualize it through marbles in jar: ![alt text](https://upload.wikimedia.org/wikipedia/commons/0/0b/Random_sampling_genetic_drift.svg) At each time step, N agents are picked up from previous generation to reproduce. Their number of offspring is proportional to their fitness, calculated as usual with **birth and death rates**. It takes thus **only one time step to go trough one generation**. Thus it is more suit- able for numerical simulations. In practice, the Moran and Wright–Fisher models give qualitatively similar results, but genetic drift runs twice as fast in the Moran model. From this perspective we can easily get that branches are less stable than in the Gillepsie scenario, for as as time goes to infinity the probability of going extinct is intuitively bigger than 0. ## Initial conditions One need to construct the world as an array of agents, which will be the ancestors of the following julia agents = [Agent( [2e-1] .* randn(1)) for i in 1:K0]  The function is then called p["tend"] -1 times. ## Scenarios You have several options available concerning the resource implemented and competition: -  mode="std" is the standard mode -  mode="grad2D" corresponds to a an ecological gradient >We are not sure whether this corresponds to the following two images -  mode="mountain" corresponds to a scenario where a mountain arises (with an ecological gradient) -  mode="split" corresponds to a scenario where the resource is splitted in two -  mode="graph" this guy is probably not working ## Parallelism You can run your script in parallel, which makes sense for large populations. To do so: julia using Distributed;addprocs() @everywhere using EvoId  Parallelism only works with Wright Fisher model. @autodocs Modules = [EvoId] Pages = ["algo/WF.jl"]   ... ... @@ -41,20 +41,18 @@ end # default initialiser """ Agent(s; ancestors=false,rates=true) Agent(s, pos;ancestors=false,rates=true) # Arguments *s is the underlying space * pos is the initial agent position. If not provided, initialises agent with 0 values everywhere # Keyword arguments * rates. Set rates=true when agents fitness needs to be updated at each time step. This is required for the Gillepsie algorithm, but not for CFM algorithm - ancestors. Set ancestors=true when you want to store ancestors traits. Agent(s; ancestors=false,rates=true) Agent(s, pos;ancestors=false,rates=true) Returns an Agent living on the underlying space s with initial position pos. If pos not provided, initialises agent with 0 values everywhere # Keyword arguments * rates. Set rates=true when agents fitness needs to be updated at each time step. This is required for the Gillepsie algorithm, but not for CFM algorithm - ancestors. Set ancestors=true when you want to store ancestors traits. """ function Agent(s::S; ancestors=false,rates=true) where {S <: AbstractSpacesTuple} T,pos = _initpos(s) ... ... @@ -123,14 +121,14 @@ Base.getindex(a::Agent,i) = a.x_history[end][i] """ get_x(a::Agent) Returns trait i of the agent Returns trait i of the agent. """ get_x(a::Agent) = a.x_history[end] @deprecate get_x(a) a[:] """$(SIGNATURES) Returns geotrait of agent a at time t Returns geotrait of agent a at time t. """ function get_geo(a::Agent{A,R,T,U,V},t::Number) where {A<:Ancestors{true},R,T,U,V} tarray = vcat(a.t_history[2:end],convert(U,t)) ... ...
 ... ... @@ -7,7 +7,7 @@ module EvoId include("Space.jl") include("Agent.jl") include("world.jl") include("World.jl") include("Sim.jl") include("metrics.jl") include("plot.jl") ... ...
File moved
 ... ... @@ -5,6 +5,42 @@ # In this file lie the function for Gillepsie algorithm """ $(TYPEDEF) Champagnat Ferriere Méléard algorithm described in [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632). This algorithm is similar to Gillepsie, excepts that it runs faster for higher number of individuals. Indeed, at every time step, only the fitness of the individual picked at random is evaluated. 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). # Example julia using EvoId,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,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, b, d, 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. Algorithm described in article DOI : 10.1016/j.tpb.2005.10.004 """ struct CFM <: AbstractAlg end export CFM ... ... @@ -22,8 +58,6 @@ end """$(SIGNATURES) Updating rule for CFM setting. algorithm described in article DOI : 10.1016/j.tpb.2005.10.004 """ function updateWorld!(w::World{A,S,T},c::CFM,b,d) where {A,S,T} # total update world ... ...
 ... ... @@ -2,6 +2,17 @@ """ $(TYPEDEF) Gillespie algorithm. Denote by b_i and d_i the birth and death rates of agent i. The total rate is given by the sum of all individual rates R(t) = \\left[ \\sum_i b_i(t) + d_i(t) \\right] A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate R. # The original article by Gillsepie: [**A general method for numerically simulating the stochastic time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub) """ struct Gillepsie <: AbstractAlg end export Gillepsie ... ...  # In this file lie the function for WF algorithm """$(TYPEDEF) Wright Fisher algorithm. In the Wright Fisher process the number of agents is constant through time. It is helpful to visualize it through marbles in jar ![alt text](https://upload.wikimedia.org/wikipedia/commons/0/0b/Random_sampling_genetic_drift.svg) At each time step, N agents are picked up from previous generation to reproduce. Their number of offspring is proportional to their fitness, calculated as usual with birth and death rates. It takes thus only one time step to go trough one generation. Thus it is more suitable for numerical simulations than CFM or Gillespie. """ struct WF <: AbstractAlg end export WF """ function updateWorld_WF!(world,newworld,C,p,update_rates!,t) If p["reflected"]=true, we reflect only first trait corresponding to geographic position """ function updateWorld_WF!(world,newworld,p,update_rates!,t) @debug "updating rates" update_rates!(world,p,t); ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!