first commit for package ABMEv.jl

src/ABMEv.jl

+__precompile__()
+
+module ABMEv
+    using Distributions,LinearAlgebra,Reexport,StatsBase
+    using SharedArrays,Distributed,LightGraphs
+
+    include("ABMEv_Agent.jl")
+    include("ABMEv_WF.jl")
+    include("ABMEv_Gillepsie.jl")
+    include("ABMEv_runworld.jl")
+    include("ABMEv_metrics.jl")
+    include("ABMEv_plot.jl")
+
+
+    @reexport using Distributions,Plots
+    export update_rates!
+    export Agent,get_fitness,get_x,get_xarray,get_xhist,
+        get_geo,get_b,get_d,increment_x!,get_inc_reflected,
+        split_move,split_merge_move,KK,tin
+    export copy,runWorld_store_WF,runWorld_store_G #,runWorld_G!,runWorld_WF!,
+    export H_discrete,findclusters,var,covgeo,hamming
+end

src/ABMEv_Agent.jl

+mutable struct Agent{T}
+    # history of traits for geotraits
+    x_history::Array{T}
+    # death rate
+    d::Float64
+    #birth rate
+    b::Float64
+end
+
+# Constructors
+
+Agent(xhist::Array) = Agent(reshape(xhist,:,1),0.,1.)
+Agent() = Agent([],0,1)
+import Base.copy
+copy(a::Agent) = Agent(a.x_history,a.d,a.b)
+copy(m::Missing) = missing
+
+
+# returns trait i of the agent
+get_x(a::Agent,i::Number) = a.x_history[Int(i):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) = 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_xarray(world::Array{Agent{T}},trait::Int) where T
+Returns trait of every agents of world in the form of an array
+"""
+get_xarray(world::Array{Agent{T}},trait::Int) where T = reshape(hcat(get_x.(world,trait)...),size(world,1),size(world,2))
+
+"""
+    function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
+This function increments current position by inc and updates xhist,
+    ONLY FOR CONTINUOUS DOMAINS
+"""
+function increment_x!(a::Agent{Float64},p::Dict;reflected=false)
+    tdim = length(p["D"])
+    if reflected
+        inc = [get_inc_reflected(get_x(a)[1],p["D"][1] *randn())]
+        if  tdim > 1
+            inc = vcat(inc,rand.(Binomial.(1,p["mu"][2:end])) .* p["D"][2:end] .* randn(tdim-1))
+        end
+    else
+        inc = rand.(Binomial.(1,p["mu"][:])) .* p["D"][:] .* randn(length(tdim))
+    end
+    a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
+ end
+
+ """
+     function increment_x!(a::Agent{Int64},inc::Array{Float64})
+ This function increments current position by inc and updates xhist,
+     ONLY FOR GRAPH TYPE DOMAINS
+ """
+ function increment_x!(a::Agent{Int64},p::Dict;reflected=false)
+     # we have to add 1 otherwise it stays on the same edge
+     inc = randomwalk(p["g"],get_x(a,1)[1],Int(p["D"][1]) + 1)
+     # the last coef of inc corresponds to last edge reached
+     a.x_history = hcat(a.x_history, reshape(inc[end:end],:,1));
+  end
+
+
+"""
+    Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
+    Careful though, we do not implement reflections
+"""
+function get_inc_reflected(x::Float64,inc::Float64,s=-1,e=1)
+    if x + inc < s
+        return 2 * ( s - x ) - inc
+    elseif  x + inc > e
+        return 2 * ( e - x ) - inc
+    else
+        return inc
+    end
+end
+
+# need to make sure that this is working correctly
+"""
+    α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64})
+Gaussian competition kernel
+"""
+function α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64})
+        return exp( - sum(sum((a1 .- a2).^n_alpha,dims=2)./ sigma_a[:].^n_alpha))
+end
+
+"""
+    function K(x::Array{Float64},K0::Float64,n_K::Float64,sigma_K::Array{Float64};μ::Float64=.0)
+Gaussian resource kernel
+"""
+function K(x::Array{Float64},K0::Float64,n_K::Float64,sigma_K::Array{Float64};μ::Float64=.0)
+    return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
+end
+
+KK(x::Array{Float64},K0::Float64,n_K::Float64,sigma_K::Array{Float64},μ1::Float64,μ2::Float64) = K(x,K0,n_K,sigma_K,μ=μ1) + K(x,K0,n_K,sigma_K,μ=μ2)
+
+"""
+    function tin(t::Float64,a::Float64,b::Float64)
+if t in [a,b) returns 1. else returns 0
+"""
+
+function tin(t::Float64,a::Float64,b::Float64)
+    return t>=a && t<b ? 1. : 0.
+end
+
+function split_move(t)
+    return .0 + 1/100*(t-20.)*tin(t,20.,120.) + tin(t,120.,Inf64)
+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

src/ABMEv_Gillepsie.jl

+# In this file lie the function for Gillepsie algorithm
+
+
+"""
+    give_birth(a::Agent,p)
+Used for Gillepsie setting
+"""
+function give_birth(a::Agent,p::Dict,reflected)
+    new_a = copy(a)
+    increment_x!(new_a,p,reflected=reflected)
+    # ! carefull you also need to change get_xhist
+    return new_a
+end
+
+function update_afterbirth_std!(world,C,idx::Int,p::Dict) where T
+    traits = get_x.(world)
+    # updating competition only the two columns corresponding to agent idx
+    for i in 1:length(traits)
+        C[i,idx] = α(traits[i],traits[idx],p["n_alpha"],p["sigma_a"])
+        C[idx,i] = C[i,idx]
+    end
+    # updating death rate only the two columns corresponding to agent idx
+    for (i,a) in enumerate(world)
+        a.d += C[i,idx] / p["K0"]
+    end
+    # Now updating new agent
+    world[idx].d = sum(C[idx,:]) / p["K0"]
+    world[idx].b = K(traits[idx][:,end],1,p["n_K"],p["sigma_K"])
+end
+
+function update_afterdeath_std!(world,C,idx::Int,p::Dict) where T
+    traits = get_x.(world)
+    # updating death rate only the two columns corresponding to agent idx
+    for (i,a) in enumerate(world)
+        a.d -= C[i,idx] / p["K0"]
+    end
+    # updating competition only the two columns corresponding to agent idx
+    for i in 1:length(traits)
+        C[i,idx] = .0
+        C[idx,i] = .0
+    end
+end
+
+"""
+    function updateWorld_G!(world,C,tspan,p)
+Updating rule for gillepsie setting.
+Returning dt drawn from an exponential distribution with parameter the total rates of events.
+"""
+function updateWorld_G!(world,C,p,update_rates!,tspan,reflected)
+    # total update world
+    world_alive = skipmissing(world)
+    idx_world = collect(eachindex(world_alive))
+    # Total rate of events
+    U = sum(get_d.(world_alive) .+ get_b.(world_alive))
+    dt = - log(rand(Float64))/U
+    events_weights = ProbabilityWeights(vcat(get_d.(world_alive),get_b.(world_alive)))
+    i_event = sample(events_weights)
+    # This is ok since length is a multiple of 2
+    I = Int(length(events_weights) / 2)
+    if i_event <= I
+        # DEATH EVENT
+        idx_offspring = idx_world[i_event]
+        world[idx_offspring] = missing
+        update_afterdeath_std!(world_alive,C,idx_offspring,p)
+    else
+        # birth event
+        idx_offspring = findfirst(ismissing,world)
+        world[idx_offspring] = give_birth(world[idx_world[i_event-I]],p,reflected)
+        update_afterbirth_std!(world_alive,C,idx_offspring,p)
+
+    end
+    return dt
+end

src/ABMEv_WF.jl

+# In this file lie the function for WF algorithm
+    """
+        function updateWorld_WF!(world,newworld,C,p,update_rates!,t,reflected)
+            If reflected=true, we reflect only first trait corresponding to geographic position
+    """
+    function updateWorld_WF!(world,newworld,C,p,update_rates!,t,reflected)
+        @debug "updating rates"
+        update_rates!(world,C,p,t);
+        # normalise to make sure that we have a probability vector
+        fitness = get_fitness.(world)
+        # we need to substract the minimum, otherwise fitness can be negative
+        prob = normalize(fitness .- minimum(fitness) .+ eps(1.),1)
+        offsprings = rand(Multinomial(length(world),prob))
+        # this guy can be done in parallel
+        @debug "filling new offspring row"
+        for (parentID,nboffspring) in collect(enumerate(offsprings))
+            for j in 1:nboffspring
+                newworld[sum(offsprings[1:(parentID-1)]) + j] = copy(world[parentID])
+            end
+        end
+        # we introduce randomness here
+        @debug "incrementing offsprings traits"
+        for w in newworld
+            increment_x!(w,p,reflected=reflected)
+        end
+    end

src/ABMEv_metrics.jl

+import StatsBase:Histogram,fit,step,normalize
+
+d1(x,y) = norm(x .- y,1)
+d2(x,y) = norm(x .- y,2)
+hamming(k::Array, l::Array) = sum(k .!= l)
+
+"""
+    H_discrete(s)
+Interconnectedness measure as in Nordbotten 2018 for discrete setup
+"""
+function H_discrete(s)
+    N = length(s)
+    H = 0
+    for x in s
+        for y in s
+            H += d2(x,y)
+        end
+    end
+    return H / N^2
+end
+
+
+"""
+    findclusters(v::Vector,allextrema =true)
+Returns a tuple with the cluster mean and its associated weight
+# Arguments
+
+"""
+function findclusters(v::Vector,allextrema =true)
+    # this function fits a histogram to some distribution and extracts its local maxima
+    h = fit(Histogram,v)
+    s = step(h.edges...)
+    x = normalize(h.weights,1)
+    dx = x[2:end] .- x[1:end-1]
+    if allextrema
+        sdx = dx[2:end] .* dx[1:end-1]
+        idx = findall(i -> i < 0, sdx)
+    else
+        idx = findall(i -> dx[i] > 0 && dx[i+1] < 0, 1:length(dx))
+    end
+    return collect(h.edges...)[idx .+ 1] .+ s/2, x[idx .+ 1]
+end
+
+import Statistics.var
+"""
+    var(world::Array{Agent};trait=:)
+If trait = 0, returns the variance of the geotrait,
+knowing that by default it is associated with position trait 1.
+If trait > 0, returns the covariance matrix, with first row geotrait and second row trait.
+# Arguments
+
+"""
+function var(world::Array{Agent};trait=1)
+    xarray = get_xarray(world,trait)
+    return var(xarray,dims=1)
+end
+"""
+    covgeo(world::Array{Agent,1},trait = 0)
+If trait = 0, returns the variance of the geotrait,
+knowing that by default it is associated with position trait 1.
+If trait > 0, returns the covariance matrix, with first row geotrait and second row trait
+# Arguments
+
+"""
+function covgeo(world::Array{Agent{T},1},trait = 0) where T
+    xarray = get_geo.(world)
+    if trait > 0
+        xstd = reshape(hcat(get_x.(world,trait)...),size(world,1),size(world,2))
+        xarray = hcat(xarray,xstd)
+    end
+    return cov(xarray)
+end
+
+"""
+    function hamming(world::Array{Agent,1})
+Returns a matrix H where H_ij = hamming(a_i,a_j).
+The hamming distance is taken through the whole history
+and functional space of the agents.
+"""
+function hamming(world::Array{Agent{T},1}) where T <: Int
+    N = size(world,1)
+    H = zeros(N,N)
+    for (i,a) in enumerate(world)
+            for (j,b) in enumerate(world)
+                    H[i,j] = hamming(get_xhist(a),get_xhist(b))
+            end
+    end
+    return H
+end

src/ABMEv_plot.jl

+using Plots
+# function designed for Wright Fisher process
+@recipe function plot(world::Array{Agent{T}},p;what=["x","H"],trait = 1) where T
+    # here  we can take xhist because every agent is updated at the same time
+    # world should correspond to a one dimensional array
+    N = size(world,1)
+    tspan = collect(1:Int(p["tend"]))
+    if "x" in what
+        @series begin
+        xarray = get_xarray(world,trait)
+            seriestype := :scatter
+            markercolor := "blue"
+            markerstrokewidth := 0
+            alpha :=.1
+            xlabel := "time"
+            ylabel := "trait value"
+            label := ""
+            markersize := 2.3/1000*size(world,1)
+            repeat(tspan,inner = N),xarray[:]
+        end
+    end
+    if "geo" in what
+        @series begin
+        xarray = get_geo.(world)
+            seriestype := :scatter
+            markercolor := "blue"
+            markerstrokewidth := 0
+            alpha :=.1
+            xlabel := "time"
+            ylabel := "trait value"
+            label := ""
+            markersize := 2.3/1000*size(world,1)
+            repeat(tspan,inner = N),xarray[:]
+        end
+    end
+    if "3dgeo" in what
+        @series begin
+        xarray = get_geo.(world)
+        yarray = get_xarray(world,2)
+            seriestype := :scatter3d
+            markercolor := "blue"
+            markerstrokewidth := 0
+            alpha :=.1
+            xlabel := "time"
+            ylabel := "geotrait"
+            zlabel := "trait value"
+            label := ""
+            markersize := 2.3/1000*size(world,1)
+            repeat(tspan,inner = N),xarray[:],yarray[:]
+        end
+    end
+    if "3d" in what
+        @series begin
+        xarray = get_xarray(world,1)
+        yarray = get_xarray(world,2)
+            seriestype := :scatter3d
+            markercolor := "blue"
+            markerstrokewidth := 0
+            alpha :=.1
+            xlabel := "time"
+            ylabel := "position"
+            zlabel := "trait value"
+            label := ""
+            markersize := 2.3/1000*size(world,1)
+            repeat(tspan,inner = N),xarray[:],yarray[:]
+        end
+    end
+    # if "H" in what
+    #     @series begin
+    #         x = get_x.(world,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,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[:,Int(i)]))
+        end
+    end
+end
+
+@recipe function plot(world::Array{Union{Agent{T},Missing}},p;what=["x","H"],trait = 1) where T
+    # here world should correspond to a multidimensional array
+    xarray = vcat(get_x.(skipmissing(world),trait)...)
+    N = size(xarray,1)
+    tspan = p["tspan"]
+    tspanarray = vcat([tspan[i]*ones(Int(p["NMax"] - count(ismissing,world[:,i]))) for i in 1:length(tspan) ]...)
+    if "x" in what
+        @series begin
+            seriestype := :scatter
+            markercolor := "blue"
+            markerstrokewidth := 0
+            alpha :=.1
+            xlabel := "time"
+            ylabel := "trait value"
+            label := ""
+            markersize := 2.3/1000*size(world,1)
+            tspanarray,xarray
+        end
+    end
+    if "H" in what
+        @series begin
+            x = get_x.(world,trait)
+            linewidth := 2
+            seriestype := :line
+            label := "Interconnectedness"
+            tspan,N^2 / 2 .* [H_discrete(x[:,i]) for i in tspan]
+        end
+    end
+end

src/ABMEv_runworld.jl

+# for now we consider that competition is local within an array
+
+function  update_rates_std!(world,C,p::Dict,t::Float64)
+    # Competition matrix
+    traits = get_x.(world)
+    # traits = get_xhist.(world)
+    N = length(traits)
+    # C = SharedArray{Float64}((N,N))
+    # Here you should do a shared array to compute in parallel
+    @sync @distributed for i in 1:(N-1)
+        for j in i+1:N
+            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
+            C[j,i] = C[i,j]
+        end
+    end
+    # Here we can do  it in parallel as well
+    for (i,a) in enumerate(world)
+        # we only update death rate
+        # this could be imptrove since \alpha is symmetric, by using a symmetric matrix
+        a.d = sum(C[i,:]) / p["K0"]
+        # /!| not ideal to assign at every time step the birth rate that is constant
+        a.b = K(traits[i][:,end],1,p["n_K"],p["sigma_K"])
+    end
+end
+
+function  update_rates_graph!(world,C,p::Dict,t::Float64)
+    for e in edges(p["g"])
+        # agents on e
+        aidx_e = findall(a -> get_x(a,1)==[e],world)
+        na = length(aidx_e)
+        for i in aidx_e
+            world[i].d = na^2
+            world[i].b = 1
+        end
+    end
+end
+
+function  update_rates_2D!(world,C,p::Dict,t::Float64)
+    # Competition matrix
+    traits = get_x.(world)
+    # traits = get_xhist.(world)
+    N = length(traits)
+    # C = SharedArray{Float64}((N,N))
+    # Here you should do a shared array to compute in parallel
+    @sync @distributed for i in 1:(N-1)
+        C[i,i] = 1.
+        for j in i+1:N
+            # be careful, for xhist what is return is an array hence no need of []
+            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
+            C[j,i] = C[i,j]
+        end
+    end
+    C[N,N] = 1.
+    # Here we can do  it in parallel as well
+    for (i,a) in enumerate(world)
+        # we only update death rate
+        # a.d = sum(C[i,:]) / K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1])
+        a.d = sum(C[i,:])
+        # /!| not ideal to assign at every time step the birth rate that is constant
+        a.b = KK(traits[i][1:1],
+                    p["K0"],
+                    p["n_K"],
+                    p["sigma_K"][1:1],
+                    split_merge_move(t),
+                    -split_merge_move(t))/K(traits[i][2:2],
+                    p["K0"],
+                    p["n_K"],
+                    p["sigma_K"][2:2])
+    end
+end
+
+function  update_rates_grad2D!(world,C,p::Dict,t::Float64)
+    # Competition matrix
+    traits = get_x.(world)
+    # traits = get_xhist.(world)
+    N = length(traits)
+    # C = SharedArray{Float64}((N,N))
+    # Here you should do a shared array to compute in parallel
+    @sync @distributed for i in 1:(N-1)
+        C[i,i] = 1.
+        for j in i+1:N
+            # be careful, for xhist what is return is an array hence no need of []
+            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
+            C[j,i] = C[i,j]
+        end
+    end
+    C[N,N] = 1.
+    # Here we can do  it in parallel as well
+    for (i,a) in enumerate(world)
+        # we only update death rate
+        a.d = sum(C[i,:])
+        a.b = K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1])
+    end
+end
+
+function  update_rates_mountain!(world,C,p::Dict,t::Float64)
+    # Competition matrix
+    traits = get_x.(world)
+    # traits = get_xhist.(world)
+    N = length(traits)
+    # C = SharedArray{Float64}((N,N))
+    # Here you should do a shared array to compute in parallel
+    @sync @distributed for i in 1:(N-1)
+        C[i,i] = 1.
+        for j in i+1:N
+            # be careful, for xhist what is return is an array hence no need of []
+            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
+            C[j,i] = C[i,j]
+        end
+    end
+    C[N,N] = 1.
+    # Here we can do  it in parallel as well
+    for (i,a) in enumerate(world)
+        # we only update death rate
+        a.d = sum(C[i,:])
+        a.b = KK(traits[i][1:1],
+                p["K0"],
+                p["n_K"],
+                p["sigma_K"][1:1],
+                split_move(t),
+                -split_move(t))/K(traits[i][2:2],
+                p["K0"],
+                p["n_K"],
+                p["sigma_K"][2:2],
+                μ=(tin(traits[i][1],-1.,0.)*(traits[i][1]+1) - tin(traits[i][1],0.,1.)*(traits[i][1]-1) ) * split_move(t))
+    end
+end
+
+function  update_rates_std_split!(world,C,p::Dict,t::Float64)
+    # Competition matrix
+    traits = get_x.(world)
+    # traits = get_xhist.(world)
+    N = length(traits)
+    # C = SharedArray{Float64}((N,N))
+    # Here you should do a shared array to compute in parallel
+    @sync @distributed for i in 1:(N-1)
+        C[i,i] = 1.
+        for j in i+1:N
+            C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
+            C[j,i] = C[i,j]
+        end
+    end
+    C[N,N] = 1.
+    # Here we can do  it in parallel as well
+    for (i,a) in enumerate(world)
+        # we only update death rate
+        # this could be imptrove since \alpha is symmetric, by using a symmetric matrix