Commit eab9cde4 authored by Victor's avatar Victor
Browse files

Introducing new abstract type to distinguish between different types of agents...

Introducing new abstract type to distinguish between different types of agents : continuous, discrete or mixed.
parent b6004ca0
...@@ -15,7 +15,7 @@ module ABMEv ...@@ -15,7 +15,7 @@ module ABMEv
@reexport using Distributions @reexport using Distributions
export update_rates! export update_rates!
export Agent,get_fitness,get_x,get_dim,get_nancestors,get_xarray,get_xhist, export Mixed,Agent,get_fitness,get_x,get_dim,get_nancestors,get_xarray,get_xhist,
get_geo,get_b,get_d,increment_x!,get_inc_reflected, get_geo,get_b,get_d,increment_x!,get_inc_reflected,
split_move,split_merge_move,KK,tin,new_world_G split_move,split_merge_move,KK,tin,new_world_G
export copy,runWorld_store_WF,runWorld_store_G #,runWorld_G!,runWorld_WF!, export copy,runWorld_store_WF,runWorld_store_G #,runWorld_G!,runWorld_WF!,
......
mutable struct Agent{T} abstract type StdAgent end
abstract type MixedAgent end
mutable struct Agent{T,U}
# history of traits for geotraits # history of traits for geotraits
x_history::Array{T} x_history::Array{U}
# death rate # death rate
d::Float64 d::Float64
#birth rate #birth rate
...@@ -8,8 +11,12 @@ mutable struct Agent{T} ...@@ -8,8 +11,12 @@ mutable struct Agent{T}
end end
# Constructors # Constructors
# This constructor should be used when one wants to impose the type of the agent (e.g. Mixed)
Agent{T}(xhist::Array{U}) where T = Agent{T,U}(reshape(xhist,:,1),0.,1.)
# This constructor is used by default
Agent(xhist::Array{U}) where U <: Number = Agent{StdAgent}(xhist)
Agent(xhist::Array{T}) where T <: Number = Agent{T}(reshape(xhist,:,1),0.,1.)
Agent() = Agent(Float64[],0.,1.) Agent() = Agent(Float64[],0.,1.)
import Base.copy import Base.copy
copy(a::Agent) = Agent(a.x_history,a.d,a.b) copy(a::Agent) = Agent(a.x_history,a.d,a.b)
...@@ -40,15 +47,15 @@ get_fitness(a::Agent) = a.b - a.d ...@@ -40,15 +47,15 @@ get_fitness(a::Agent) = a.b - a.d
get_dim(a::Agent) = size(a.x_history,1) get_dim(a::Agent) = size(a.x_history,1)
get_nancestors(a::Agent) = size(a.x_history,2) get_nancestors(a::Agent) = size(a.x_history,2)
""" """
get_xarray(world::Array{Agent{T}},trait::Int) where T get_xarray(world::Array{Agent},trait::Int)
Returns trait of every agents of world in the form of an array which dimensions corresponds to the input. 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 a timeseries. Particularly suited for an array world corresponding to a timeseries.
""" """
get_xarray(world::Array{Agent{T}},trait::Int) where T = reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2)) get_xarray(world::Array{Agent},trait::Int) = reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2))
""" """
get_xhist(world::Vector{Agent{T}},geotrait = false) get_xhist(world::Vector{Agent},geotrait = false)
Returns the trait history of every agents of world in the form of an 3 dimensional array, Returns the trait history of every agents of world in the form of an 3 dimensional array,
with with
- first dimension as the agent index - first dimension as the agent index
...@@ -58,7 +65,7 @@ If geotrait = true, then a last trait dimension is added, corresponding to geotr ...@@ -58,7 +65,7 @@ If geotrait = true, then a last trait dimension is added, corresponding to geotr
Note that because number of ancestors are different between agents, we return an array which size corresponds to the minimum of agents ancestors, 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 and return the last generations, dropping the youngest ones
""" """
function get_xhist(world::Vector{Agent{T}},geotrait = false) where T function get_xhist(world::Vector{Agent},geotrait = false)
hist = minimum(get_nancestors.(world)) hist = minimum(get_nancestors.(world))
ntraits = get_dim(first(world)); ntraits = get_dim(first(world));
xhist = zeros(length(world), hist, ntraits + geotrait); xhist = zeros(length(world), hist, ntraits + geotrait);
...@@ -73,11 +80,11 @@ end ...@@ -73,11 +80,11 @@ end
""" """
increment_x!(a::Agent{Float64},p::Dict) increment_x!(a::Agent{StdAgent,U},p::Dict)
This function increments current position by inc and updates xhist, This function increments agent by random numbers specified in p
ONLY FOR CONTINUOUS DOMAINS ONLY FOR CONTINUOUS DOMAINS
""" """
function increment_x!(a::Agent{Float64},p::Dict) function increment_x!(a::Agent{StdAgent,U},p::Dict)
tdim = length(p["D"]) tdim = length(p["D"])
reflected = haskey(p,"reflected") ? p["reflected"] : false reflected = haskey(p,"reflected") ? p["reflected"] : false
if reflected if reflected
...@@ -93,23 +100,30 @@ function increment_x!(a::Agent{Float64},p::Dict) ...@@ -93,23 +100,30 @@ function increment_x!(a::Agent{Float64},p::Dict)
end end
""" """
function increment_x!(a::Agent{Int64},inc::Array{Float64}) function increment_x!(a::Agent{MixedAgent,U},p::Dict)
This function increments current position by inc and updates xhist, This function increments first trait of agent with integer values, that are automatically reflected between 1 and p["nodes"].
ONLY FOR GRAPH TYPE DOMAINS ONLY FOR GRAPH TYPE DOMAINS
""" """
function increment_x!(a::Agent{Int64},p::Dict) # function increment_x!(a::Agent{Mixed},p::Dict)
# we have to add 1 otherwise it stays on the same edge # tdim = length(p["D"])
inc = randomwalk(p["g"],get_x(a,1),Int(p["D"][1]) + 1) # inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn(),1,p["nodes"])]
# the last coef of inc corresponds to last edge reached # if tdim > 1
a.x_history = hcat(a.x_history, reshape(inc[end:end],:,1)); # inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
end # end
# else
# # inc = yes no mutation * mutation
""" # inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
get_inc_reflected(x::Float64,inc::Float64,s=-1,e=1) # end
# a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
# end
# end
"""
get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D) Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
""" """
function get_inc_reflected(x::Float64,inc::Float64,s=-1,e=1) function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
if x + inc < s if x + inc < s
inc = 2 * ( s - x ) - inc inc = 2 * ( s - x ) - inc
elseif x + inc > e elseif x + inc > e
...@@ -122,43 +136,43 @@ end ...@@ -122,43 +136,43 @@ end
# need to make sure that this is working correctly # need to make sure that this is working correctly
""" """
α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64}) α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Generalised gaussian competition kernel Generalised gaussian competition kernel
""" """
function α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64}) 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)) return exp( -.5* sum(sum((a1 .- a2).^n_alpha,dims=2)./ sigma_a[:].^n_alpha))
end end
""" """
α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64}) α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Default Gaussian competition kernel Default Gaussian competition kernel
""" """
α(a1::Array{Float64},a2::Array{Float64},sigma_a::Array{Float64}) = α(a1,a2,2,sigma_a) α(a1::Array{Number},a2::Array{Number},sigma_a::Array{Number}) = α(a1,a2,2,sigma_a)
""" """
K(x::Array{Float64},K0::Float64,μ::Array{Float64},sigma_K::Array{Float64}) K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
Gaussian resource kernel Gaussian resource kernel
""" """
function K(x::Array{Float64},K0::Float64,μ::Array{Float64},sigma_K::Array{Float64}) 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*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
return K0 * pdf(MvNormal(μ,sigma_K),x) return K0 * pdf(MvNormal(μ,sigma_K),x)
end end
""" """
K(x::Array{Float64},K0::Float64,sigma_K::Array{Float64}) K(x::Array{Number},K0::Number,sigma_K::Array{Number})
Gaussian resource kernel with mean 0. Gaussian resource kernel with mean 0.
""" """
function K(x::Array{Float64},K0::Float64,sigma_K::Array{Float64}) 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*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
return K0 * pdf(MvNormal(sigma_K),x) return K0 * pdf(MvNormal(sigma_K),x)
end end
KK(x::Array{Float64},K0::Float64,n_K::Float64,sigma_K::Array{Float64},μ1::Float64,μ2::Float64) = K(x,K0/2,μ1,sigma_K) + K(x,K0/2,μ2,sigma_K) 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::Float64,a::Float64,b::Float64) function tin(t::Number,a::Number,b::Number)
if t in [a,b) returns 1. else returns 0 if t in [a,b) returns 1. else returns 0
""" """
function tin(t::Float64,a::Float64,b::Float64) function tin(t::Number,a::Number,b::Number)
return t>=a && t<b ? 1. : 0. return t>=a && t<b ? 1. : 0.
end end
......
...@@ -63,7 +63,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second ...@@ -63,7 +63,7 @@ If trait > 0, returns the covariance matrix, with first row geotrait and second
# Arguments # Arguments
""" """
function covgeo(world::Array{U,1},trait = 0) where U <: Union{Missing,Agent{T}} where T function covgeo(world::Array{U,1},trait = 0) where U <: Union{Missing,Agent}
world = collect(skipmissing(world)) world = collect(skipmissing(world))
xarray = get_geo.(world) xarray = get_geo.(world)
if trait > 0 if trait > 0
...@@ -79,7 +79,7 @@ Returns a matrix H where H_ij = hamming(a_i,a_j). ...@@ -79,7 +79,7 @@ Returns a matrix H where H_ij = hamming(a_i,a_j).
The hamming distance is taken through the whole history The hamming distance is taken through the whole history
and functional space of the agents. and functional space of the agents.
""" """
function hamming(world::Array{Agent{T},1}) where T <: Int function hamming(world::Array{Agent,1}) where T <: Int
N = size(world,1) N = size(world,1)
H = zeros(N,N) H = zeros(N,N)
for (i,a) in enumerate(world) for (i,a) in enumerate(world)
......
using RecipesBase using RecipesBase
# function designed for Wright Fisher process # function designed for Wright Fisher process
@recipe function plot(world::Array{U},p;what=["x","H"],trait = 1) where U <: Union{Missing,Agent{T}} where T @recipe function plot(world::Array{U},p;what=["x","H"],trait = 1) where U <: Union{Missing,Agent}
# We reduce time interval if it is too big # We reduce time interval if it is too big
if size(world,2)*size(world,1) > 1e6 && size(world,2) >= 200 if size(world,2)*size(world,1) > 1e6 && size(world,2) >= 200
p = copy(p) p = copy(p)
......
""" """
generalised_gaussian(x::Float64,mu::Float64,sigma::Float64,epsilon::Float64) generalised_gaussian(x::Number,mu::Number,sigma::Number,epsilon::Number)
""" """
function generalised_gaussian(x::Float64,mu::Float64,sigma::Float64,epsilon::Float64) function generalised_gaussian(x::Number,mu::Number,sigma::Number,epsilon::Number)
return exp( -.5 * ((x-mu) / sigma)^epsilon) return exp( -.5 * ((x-mu) / sigma)^epsilon)
end end
""" """
gaussian(x::Float64,mu::Float64,sigma::Float64) = generalised_gaussian(x,mu,sigma,2) gaussian(x::Number,mu::Number,sigma::Number) = generalised_gaussian(x,mu,sigma,2)
""" """
gaussian(x::Float64,mu::Float64,sigma::Float64) = generalised_gaussian(x,mu,sigma,2.) gaussian(x::Number,mu::Number,sigma::Number) = generalised_gaussian(x,mu,sigma,2.)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment