version 0.1.8, with new documentation, old competition and birth function...

version 0.1.8, with new documentation, old competition and birth function deleted because now passed by the used.

Project.toml

 name = "ABMEv"
 uuid = "837ac870-fb52-4b0c-9a0e-030f2f36f5ed"
 authors = ["Victor Boussange "]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

README.md

@@ -17,13 +17,26 @@ Pkg.add("ABMEv#no_C_matrix")
 using ABMEv
 ```
 
+## The `Agent` structure
+This package is an Agent Based Model, where the atomic structure is the `Agent`, which is defined by a history of traits, a death rate and a birth rate.
+```julia
+mutable struct Agent{T,U}
+    # history of traits for geotraits
+    x_history::Array{U}
+    # 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
-- ```reflected=>true``` if ```true``` then reflection occurs on the first trait -which should stand for geographic position- in the domain $` [-1,1] `$
-- ```"alpha" => α``` is the competition function
--```"K" => K``` is the birth rate
-- ```"tend" => 1.5``` is the time to end simulation
+- ```"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
@@ -49,22 +62,31 @@ Two type of simulation algorithm can be used
 ### Gillepsie algorithm
 
 ```julia
-a = 0;
-sigma_K = .9;
-sigma_a = .7;
-K0 = 1000;
-K(X) = gaussian(X[1],0.,sigma_K)
-α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
-p_default = Dict(
-        "alpha" => α,
-        "K" => K,
-        "D" => [1e-2],
-        "mu" => [.1],
-        "tend" => 1.5,
-        "NMax" => Int(10000))
-na_init = K0
-world0 = new_world_G(na_init,p_default,spread = [.01 .01], offset = [-.5 -.5])
-worldall,p_default["tspan"] = runWorld_store_G(p_default,world0)
+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)
 ```
 As of now, no mode is implemented.
 #### Specific parameters
@@ -112,18 +134,45 @@ Parallelism only works with Wright Fisher model.
 You can access properties of the agent using the following functions
 - `get_xarray(world::Array{Agent{T}},trait::Int) where T`
 
-Returns trait of every agents of world in the form of an array
+Returns trait of every agents of `world` in the form of an array which dimensions corresponds to the input.
+Particularly suited for an array `world` corresponding to one time step, or a time series from a **Wright Fisher simulation**.
+>TODO: implement the possibility of getting geotrait
+
+- `get_xarray(world::Array{T,1},geotrait=false) where {T <: Agent}`
+
+Returns every traits of every agents of world in the form of an array
+If geotrait = true, then a last trait dimension is added, corresponding to geotrait.
+
+- `get_xhist(world::Vector{Agent},geotrait = false)`
+
+Returns the trait history of every agents of world in the form of an 3 dimensional array,
+with
+- - first dimension as the agent index
+- - second as time index
+- - third as trait index
+If geotrait = true, then a last trait dimension is added, corresponding to geotrait.
+Note that because number of ancestors are different between agents, we return an array which size corresponds to the minimum of agents ancestors,
+and return the last generations, dropping the youngest ones
+
+- `world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}`
+
+Converts the array of agent world to a datafram, where each column corresponds to a trait of the
+agent, and an extra column captures fitness.
+Each row corresponds to an agent
 
+### Other accessors
 > TODO: describe the following accessors
 ```julia
-get_x(a::Agent,i::Number) = a.x_history[Int(i):Int(i),end]
+get_x(a::Agent,i::Number) = a.x_history[Int(i),end]
 get_x(a::Agent) = a.x_history[:,end]
-get_xhist(a::Agent,i::Number) = a.x_history[Int(i):Int(i),:]
+get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
 get_xhist(a::Agent) = a.x_history
 get_geo(a::Agent) = sum(get_xhist(a,1))
 get_d(a::Agent) = a.d
 get_b(a::Agent) = a.b
 get_fitness(a::Agent) = a.b - a.d
+get_dim(a::Agent) = size(a.x_history,1)
+get_nancestors(a::Agent) = size(a.x_history,2)
 ```
 
 ## Plotting
@@ -136,12 +185,13 @@ using Plots;pyplot()
 Plots.plot(worldall,p_default,what=["x"],trait=2)
 ```
 You can specify what you want to plot in the array ```what```:
-- ```x``` plots the component specified by ```trait=2```
-- ```geo``` plots geotrait, computed from first component
-- ```3dgeo``` plots a 3d diagram with x axis as geotrait and y axis as the second component
-- ```3d``` plots a 3d diagram with first and second component as x and y axis
-- ```var``` plots the variance of the  component specified by ```trait=2```
-- ```vargeo``` plots the variance of the geotrait
+- ```"x"``` returns a scatter plot `(xaxis = time, yaxis = trait value)` where trait component is specified by ```trait=2```
+- `"xs"` only works for agent type `MixedAgent`, because it needs a discrete geographical space. It returns a scatter plot `(xaxis = geographical component, yaxis = trait value)`
+- ```"geo"``` returns a scatter plot where trait is the *geotrait* computed from first component
+- ```"3dgeo"``` plots a 3d diagram with x axis as geotrait and y axis as the second component
+- ```"3d"``` plots a 3d diagram with first and second component as x and y axis
+- ```"var"``` plots the variance of the  component specified by ```trait=2``` :question: with respect to time?
+- ```"vargeo"``` plots the variance of the geotrait 
 
 ## Developping the code
 I recommend to first clone your branch in the directory you like best, and then to 

examples/gillepsie_simple.jl

+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 = 100.
+using UnPack
+p = Dict{String,Any}()
+@pack! p = K,alpha,D,mu,NMax,dt_saving,tend
+
+## Initial conditions
+agent0 = [Agent(- Float32(1.5) .+ .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)

src/ABMEv.jl

@@ -17,7 +17,7 @@ module ABMEv
     export update_rates!
     export MixedAgent,StdAgent,Agent,get_fitness,get_x,get_dim,get_nancestors,get_xarray,get_xhist,
         get_geo,get_b,get_d,increment_x!,get_inc_reflected,world2df
-        split_move,split_merge_move,KK,tin,new_world_G
+        split_move,split_merge_move,tin,new_world_G
     export copy,runWorld_store_WF,runWorld_store_G #,runWorld_G!,runWorld_WF!,
     export H_discrete,findclusters,var,covgeo,hamming
     export update_afterbirth_std!,update_afterdeath_std!

src/ABMEv_Agent.jl

@@ -60,7 +60,7 @@ get_xarray(world::Array{T},trait::Int) where {T <: Agent}= reshape(hcat(get_x.(w
 Returns every traits of every agents of world in the form of an array
 If geotrait = true, then a last trait dimension is added, corresponding to geotrait.
 """
-function get_xarray(world::Array{T,1},geotrait=false) where {T <: Agent}
+function get_xarray(world::Array{T,1},geotrait::Bool=false) where {T <: Agent}
     xarray = hcat(get_x.(world)...)
     if geotrait
         xarray = vcat( xarray, get_geo.(world)')
@@ -91,6 +91,25 @@ function get_xhist(world::Vector{T},geotrait = false) where {T <: Agent}
     return xhist
 end
 
+"""
+    world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}
+Converts the array of agent world to a datafram, where each column corresponds to a trait of the
+agent, and an extra column captures fitness.
+Each row corresponds to an agent
+"""
+function world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}
+    xx = get_xarray(world)
+    dfw = DataFrame(:f => get_fitness.(world))
+    for i in 1:size(xx,1)
+        dfw[Meta.parse("x$i")] = xx[i,:]
+    end
+    if geotrait
+        dfw[:g] = get_geo.(world)
+    end
+    return dfw
+end
+
+
 
 """
 increment_x!(a::Agent{StdAgent,U},p::Dict)
@@ -115,7 +134,8 @@ function increment_x!(a::Agent{StdAgent,U},p::Dict) where U
  """
      function increment_x!(a::Agent{MixedAgent,U},p::Dict)
  This function increments first trait of agent with integer values, that are automatically reflected between 1 and p["nodes"].
-     ONLY FOR GRAPH TYPE DOMAINS
+Other traits are incremented as usual.
+TODO : make it work for a graph type landscape, where domain is not a line anymore.
  """
  function increment_x!(a::Agent{MixedAgent,U},p::Dict) where U
      tdim = length(p["D"])
@@ -143,39 +163,6 @@ function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
     get_inc_reflected(x,inc,s,e)
 end
 
-# need to make sure that this is working correctly
-"""
-    α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
-Generalised gaussian competition kernel
-"""
-function α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
-        return exp( -.5* sum(sum((a1 .- a2).^n_alpha,dims=2)./ sigma_a[:].^n_alpha))
-end
-"""
-    α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
-Default Gaussian competition kernel
-"""
-α(a1::Array{Number},a2::Array{Number},sigma_a::Array{Number}) = α(a1,a2,2,sigma_a)
-
-"""
-    K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
-Gaussian resource kernel
-"""
-function K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
-    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
-    return K0 * pdf(MvNormal(μ,sigma_K),x)
-end
-"""
-    K(x::Array{Number},K0::Number,sigma_K::Array{Number})
-Gaussian resource kernel with mean 0.
-"""
-function K(x::Array{Number},K0::Number,sigma_K::Array{Number})
-    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
-    return K0 * pdf(MvNormal(sigma_K),x)
-end
-
-KK(x::Array{Number},K0::Number,n_K::Number,sigma_K::Array{Number},μ1::Number,μ2::Number) = K(x,K0/2,μ1,sigma_K) + K(x,K0/2,μ2,sigma_K)
-
 """
     function tin(t::Number,a::Number,b::Number)
 if t in [a,b) returns 1. else returns 0
@@ -192,21 +179,3 @@ end
 function split_merge_move(t)
     return .0 + 1/30*(t-10.)*tin(t,10.,40.) + tin(t,40.,70.) + (1- 1/30*(t-70.))*tin(t,70.,100.)
 end
-
-"""
-    world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}
-Converts the array of agent world to a datafram, where each column corresponds to a trait of the
-agent, and an extra column captures fitness.
-Each row corresponds to an agent
-"""
-function world2df(world::Array{T,1}; geotrait = false) where {T <: Agent}
-    xx = get_xarray(world)
-    dfw = DataFrame(:f => get_fitness.(world))
-    for i in 1:size(xx,1)
-        dfw[Meta.parse("x$i")] = xx[i,:]
-    end
-    if geotrait
-        dfw[:geo] = get_geo.(world)
-    end
-    return dfw
-end

src/ABMEv_plot.jl

@@ -38,6 +38,31 @@ import KernelDensity:kde,pdf
             tspan_ar[:],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 = []
+        world_df_g = groupby(world2df(world_sm),:x1)
+        for world_df in world_df_g
+            x = world_df.x2
+            append!(d_i,pdf(kde(x),x))
+        end
+        # TODO: we stopped here
+        @series begin
+        xarray = get_xarray(world_sm,trait)
+            seriestype := :scatter
+            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
+            # markercolor := :blue
+            markerstrokewidth := 0
+            alpha :=1.
+            xlabel := "time"
+            ylabel := "trait value"
+            label := ""
+            grid := false
+            # markersize := 2.3/1000*size(world_sm,1)
+            tspan_ar[:],xarray[:]
+        end
+    end
     if "geo" in what
         @series begin
         xarray = get_geo.(world_sm)