Commit 67c2f4a7 by Victor

parent 5614f4ca
Pipeline #75039 passed with stage
in 23 minutes and 41 seconds
 ... ... @@ -6,9 +6,9 @@ makedocs(sitename="ABMEv.jl", authors = "Victor Boussange", pages = [ "Home" => "index.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))], "Examples" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"examples"),join=true))], "Library" => Any[ "Public" => "lib/public.md", # "Internals" => map( ... ...
 # A first model of meta population In this script, we model the evolution of a population where agents are simply defined by their position on some landscape. We implement the simplest possible birth and death function. ## The landscape Let's start by a linear landscape. We define a discrete segment of length 9, with reflecting boundary conditions.
 # Modelling agents with a 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. The genotype space is inspired from the article [The architecture of an empirical genotype-phenotype map](http://doi.wiley.com/10.1111/evo.13487) ## Defining the space julia ##### Genotype space##### dim_neutr = 1000 magicprop = 523728 / 32896 g = SimpleGraph{Int16}(dim_neutr,round(Int16,dim_neutr * magicprop)) initnode = argmax(eigenvector_centrality(g)) # This is the central node the we will use to instantiate the populations myspace = (DiscreteSegment(Int8(1),Int8(nodes)),GraphSpace(g)) # union of vector spaces  ## Defining birth, death processes and mutation julia K0 = 1000 sigma_K = 1.; sigma_a = .8; mu = [1.,1.] b(X,t) = 1 / nodes d(X,Y,t) = (X[1] ≈ Y[1]) / K0 D = (5e-1,1.4825)  ## Final set up julia NMax = 2000 # tend = 1.5 tend = 3000 p_default = Dict{String,Any}();@pack! p_default = d,b,NMax,mu myagents = [Agent(myspace,(rand(Int8(1):Int8(nodes)),initnode),ancestors=true,rates=true) for i in 1:round(K0/nodes)] w0 = World(myagents,myspace,p_default,0.) 
 # Modelling Sympatric speciation This script aims at reproducing the 1999 article of Doebeli [On The Origin of Species By Sympatric Speciation](http://www.nature.com/articles/22521).
 # ABMEv.jl Documentation 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 This is a suite for simulating the evolutionary dynamics of a population in a multidimensional space. The population is modelled at the level of the individual - this suite hence falls in the realm of *Agent Based Modelling*. The purpose of this package is to provide a numerical laboratory for evolutionary dynamics, supplying - a flexible atomic structure of agents and underlying evolutionary space - algorithms and update rules for the simulations - 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. Vector spaces are used to define birth and death processes, as well as mutation processes. ## Getting started julia using ABMEv  ## Parameters of the simulation Parameters are stored in the parameter dictionary p ### General parameters - "reflected"=>true: if true then reflection occurs on the first trait -which should stand for geographic position. Depending on the agent type, reflections occurs in the domain  [-1,1]  or between nodes 1 and p["nodes"] - "alpha" => α: is the competition function - "K" => K: is the birth rate - "tend" => 1.5: is the time to end simulation :warning: Check how to define functions α and K in the algorithm section. ### Mutation If anisotropy in mutation, the following parameters should be declared as arrays where each entry corresponds to a dimension. - mu The probability of mutation. - D If mutation happens on the agent, the increment follows \mathcal{N}_{ 0, D} ### Birth #### Growth - K is the birth coefficient ( b(x) = K(x) ) ### Death #### Competition - Competition between agent with trait x and y is defined as α(x,y) - Death coefficient is defined as d(x^{(i)}) = \sum_j^{N(t)} \alpha(x^{(i)},x^{(j)}) ### Fitness Fitness is defined internally as b - d. > TODO b here is confounded with K. ## Launching simulation Two type of simulation algorithm can be used @contents Pages = ["manual/gillepsie.md", "manual/wright_fisher.md"] Depth = 5  ### Gillepsie algorithm julia @repl using ABMEv ## Defining parameters 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 D = [1e-2] #mutation range mu = [1.] #probability of mutation NMax = 2000 #number of individual dt_saving = 1.0 #time step saving tend = 1000. using UnPack p = Dict{String,Any}() @pack! p = K,alpha,D,mu,NMax,dt_saving,tend ## Initial conditions agent0 = [Agent(.1 .* randn(Float32,1)) for i in 1:K0] world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - K0))) ## launch simulation worldall,p["tspan"] = runWorld_store_G(p,world0); using Plots Plots.plot(worldall,p)  ## Tutorial We strongly advise to have a look at the tutorial section. ### Wright Fisher algorithm julia sigma_K = .9; sigma_a = 1.251; K0 = 1000; # K(X) = gaussian(X[1],0.,sigma_K) K(X) = 1 - 0.125 * sum(X.^2) α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0 # α(X,Y) = 0. p = Dict( "alpha" => α, "K" => K, "D" => [1e-2], "mu" => [.1], "tend" => 10.) na_init = K0 agents = [Agent( [1e-2] .* randn(1) .- .5) for i in 1:K0] @time worldall_test,p["tspan"] = runWorld_store_WF(p,agents,mode="std");  ## How it works There 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 - Define mutation function - [Define initial population state and time](manual/world) - [Run the simulation according to some updating algorithm](manual/run_world.md) - [Obtain a summary of the population state](manual/simulation.md) ### 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) ## References - [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632) - [Champagnat and Ferriere second article - 2008](https://www.tandfonline.com/doi/full/10.1080/15326340802437710) ## Similar packages: [Agents.jl](https://juliadynamics.github.io/Agents.jl/) It would be worth to have a look! It has been designed by Ali R. Vahdati, from UZH. [Agents.jl](https://juliadynamics.github.io/Agents.jl/) This package is oriented towards general ABM modelling, and miss many of the specificities of ABMEv.jl , such as vector spaces operations, update algorithm, and analysis tools.
 # 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). !!! 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.
 ... ... @@ -6,59 +6,30 @@ [**A general method for numerically simulating the stochastic time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub) ### Rates Each individual is assigned a birth  b_i  and death  d_i  rate. The total rate is given by the sum of all individual rates ### 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 > This has to be checked, we are still not hundred percent sure ### Time steps An event is exponentiallly distributed in time, with parameter \lambda = U(t). This makes events memoryless, meaning that the probability of having a birth or death event is always the same, no matter when (P(X > s_t | X > t) = P(X > s) . > Let B(t) = \sum_i b_i(t) and  D(t) = \sum_i d_i(t). Let T_b, T_d the time for a birth or death event to occur. Then we have P(T_b < T_d) = \frac{B(t)}{B(t) + D(t)} (competing exponentials). #### Inversion method Let U be an \mathcal{U}_{(0,1)}-distributed random variable and F \colon \R \to [0,1] be a distribution function. Then we have math P(I_F(U) \leq x ) = P(U \leq F(x)) = F(x)  Thanks to the ***inversion method*** we get the incremental time step dt, exponentially distributed with parameter \lambda = R(t), as math dt(\omega) = -\frac{\log(U(\omega))}{R(t)} \iff X(\omega) = \exp(-U(t)dt(\omega))  ## Initialize @docs new_world_G  ## Run @docs runWorld_store_G  ## Scenarios As of now, no mode is implemented. For further examples, check the folder examples in source code. ## Specific parameters - dt_saving = 10. will allow to save the world every 10. time steps. If not specified, the algorithm will return first and last time step world. - NMax Maximum number of individuals that can be attained. If attained, then the programm stops. ## Developping ### Efficiency The simulation are still very long. :flushed: How to improve it? - We think it would be more efficient if we found an other way of incrementing mutations !!! tip "Inversion method" Let B(t) = \sum_i b_i(t) and  D(t) = \sum_i d_i(t). Let T_b, T_d the time for a birth or death event to occur. Then we have P(T_b < T_d) = \frac{B(t)}{B(t) + D(t)} (competing exponentials). Let U be an \mathcal{U}_{(0,1)}-distributed random variable and F \colon \R \to [0,1] be a distribution function. Then we have math P(I_F(U) \leq x ) = P(U \leq F(x)) = F(x)  ### Parallelism > For now there is no parallelism implemented for one run > But we think we should rather set up a pmap or the macro @Threads to explore parameter space Thanks to the ***inversion method*** we get the incremental time step dt, exponentially distributed with parameter \lambda = R(t), as math dt(\omega) = -\frac{\log(U(\omega))}{R(t)} \iff X(\omega) = \exp(-U(t)dt(\omega))  @autodocs Modules = [ABMEv] ... ...
 # Simulation A Simulation object is return by the function run! - dt_saving = 10. will allow to save the world every 10. time steps. If not specified, the algorithm will return first and last time step world. - NMax Maximum number of individuals that can be attained. If attained, then the programm stops. # `@autodocs Modules = [ABMEv] Pages = ["ABMEv_Sim.jl"] ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!