ABMEv_Agent.jl 5.92 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
mutable struct Agent{T}
    # history of traits for geotraits
    x_history::Array{T}
    # death rate
    d::Float64
    #birth rate
    b::Float64
end

# Constructors

12
13
Agent(xhist::Array{T}) where T <: Number = Agent{T}(reshape(xhist,:,1),0.,1.)
Agent() = Agent(Float64[],0.,1.)
14
15
16
17
import Base.copy
copy(a::Agent) = Agent(a.x_history,a.d,a.b)
copy(m::Missing) = missing

18
19
20
21
22
23
24
25
26
27
28
29
"""
    function new_world_G(nagents::Int,p::Dict; spread = 1., offset = 0.)
Returns an array of type Array{Union{Missing,Agent}} initialised with normal distribution.
Only relevant for Gillepsie algorithm as of now.
"""
function new_world_G(nagents::Int,p::Dict; spread = 1., offset = 0.)
    typeof(spread) <: Array ? spread = spread[:] : nothing;
    typeof(offset) <: Array ? offset = offset[:] : nothing;
    agent0 = [Agent( spread  .* randn(length(spread)) .+ offset) for i in 1:nagents]
    world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - nagents)))
    return world0
end
30
31

# returns trait i of the agent
32
get_x(a::Agent,i::Number) = a.x_history[Int(i),end]
33
get_x(a::Agent) = a.x_history[:,end]
34
get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
35
36
37
38
39
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
40
41
get_dim(a::Agent) = size(a.x_history,1)
get_nancestors(a::Agent) = size(a.x_history,2)
42
43
"""
    get_xarray(world::Array{Agent{T}},trait::Int) where T
44
45
46
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.

47
"""
48
get_xarray(world::Array{Agent{T}},trait::Int) where T = reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2))
49

50
51
52
53
54
55
56
57
58
59
60
61
62
"""
    get_xhist(world::Vector{Agent{T}},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
"""
function get_xhist(world::Vector{Agent{T}},geotrait = false) where T
    hist = minimum(get_nancestors.(world))
Victor's avatar
Victor committed
63
64
65
    ntraits = get_dim(first(world));
    xhist = zeros(length(world), hist, ntraits + geotrait);
    for (i,a) in enumerate(world)
66
67
68
69
70
71
72
73
74
        xhist[i,:,1:end-geotrait] = get_xhist(a)[:,end-hist+1:end]';
        if geotrait
            xhist[i,:,ntraits+geotrait] = cumsum(get_xhist(a,1))[end-hist+1:end]
        end
    end
    return xhist
end


75
"""
Victor's avatar
Victor committed
76
increment_x!(a::Agent{Float64},p::Dict)
77
    This function increments current position by inc and updates xhist,
78
79
    ONLY FOR CONTINUOUS DOMAINS
"""
Victor's avatar
Victor committed
80
function increment_x!(a::Agent{Float64},p::Dict)
81
    tdim = length(p["D"])
Victor's avatar
Victor committed
82
    reflected = haskey(p,"reflected") ? p["reflected"] : false
83
    if reflected
Victor Boussange's avatar
Victor Boussange committed
84
        inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
85
        if  tdim > 1
Victor Boussange's avatar
Victor Boussange committed
86
            inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
87
88
        end
    else
Victor Boussange's avatar
Victor Boussange committed
89
90
        # inc = yes no mutation * mutation
        inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
91
92
93
94
95
96
97
98
99
    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
 """
Victor's avatar
Victor committed
100
 function increment_x!(a::Agent{Int64},p::Dict)
101
     # we have to add 1 otherwise it stays on the same edge
102
     inc = randomwalk(p["g"],get_x(a,1),Int(p["D"][1]) + 1)
103
104
105
106
107
108
     # the last coef of inc corresponds to last edge reached
     a.x_history = hcat(a.x_history, reshape(inc[end:end],:,1));
  end


"""
109
get_inc_reflected(x::Float64,inc::Float64,s=-1,e=1)
110
111
112
113
    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)
    if x + inc < s
114
        inc = 2 * ( s - x ) - inc
115
    elseif  x + inc > e
116
        inc = 2 * ( e - x ) - inc
117
118
119
    else
        return inc
    end
120
    get_inc_reflected(x,inc,s,e)
121
122
123
124
125
end

# need to make sure that this is working correctly
"""
    α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64})
Victor Boussange's avatar
Victor Boussange committed
126
Generalised gaussian competition kernel
127
128
"""
function α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64})
Victor Boussange's avatar
Victor Boussange committed
129
        return exp( -.5* sum(sum((a1 .- a2).^n_alpha,dims=2)./ sigma_a[:].^n_alpha))
130
end
Victor Boussange's avatar
Victor Boussange committed
131
132
133
134
135
"""
    α(a1::Array{Float64},a2::Array{Float64},n_alpha::Float64,sigma_a::Array{Float64})
Default Gaussian competition kernel
"""
α(a1::Array{Float64},a2::Array{Float64},sigma_a::Array{Float64}) = α(a1,a2,2,sigma_a)
136
137

"""
138
    K(x::Array{Float64},K0::Float64,μ::Array{Float64},sigma_K::Array{Float64})
139
140
Gaussian resource kernel
"""
141
142
143
144
145
146
147
148
149
150
151
function K(x::Array{Float64},K0::Float64,μ::Array{Float64},sigma_K::Array{Float64})
    # 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{Float64},K0::Float64,sigma_K::Array{Float64})
Gaussian resource kernel with mean 0.
"""
function K(x::Array{Float64},K0::Float64,sigma_K::Array{Float64})
    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
    return K0 * pdf(MvNormal(sigma_K),x)
152
153
end

154
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)
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171

"""
    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