ABMEv_Agent.jl 6.37 KB
Newer Older
1
2
3
4
abstract type StdAgent end
abstract type MixedAgent end

mutable struct Agent{T,U}
5
    # history of traits for geotraits
6
    x_history::Array{U}
7
8
9
10
11
12
13
    # death rate
    d::Float64
    #birth rate
    b::Float64
end

# Constructors
14
15
16
17
18
# 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)
19

20
Agent() = Agent(Float64[],0.,1.)
21
22
23
24
import Base.copy
copy(a::Agent) = Agent(a.x_history,a.d,a.b)
copy(m::Missing) = missing

25
26
27
28
29
30
31
32
33
34
35
36
"""
    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
37
38

# returns trait i of the agent
39
get_x(a::Agent,i::Number) = a.x_history[Int(i),end]
40
get_x(a::Agent) = a.x_history[:,end]
41
get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
42
43
44
45
46
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
47
48
get_dim(a::Agent) = size(a.x_history,1)
get_nancestors(a::Agent) = size(a.x_history,2)
49
"""
50
    get_xarray(world::Array{Agent},trait::Int)
51
52
53
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.

54
"""
55
get_xarray(world::Array{Agent},trait::Int) = reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2))
56

57
"""
58
    get_xhist(world::Vector{Agent},geotrait = false)
59
60
61
62
63
64
65
66
67
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
"""
68
function get_xhist(world::Vector{Agent},geotrait = false)
69
    hist = minimum(get_nancestors.(world))
Victor's avatar
Victor committed
70
71
72
    ntraits = get_dim(first(world));
    xhist = zeros(length(world), hist, ntraits + geotrait);
    for (i,a) in enumerate(world)
73
74
75
76
77
78
79
80
81
        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


82
"""
83
84
increment_x!(a::Agent{StdAgent,U},p::Dict)
    This function increments agent by random numbers specified in p
85
86
    ONLY FOR CONTINUOUS DOMAINS
"""
87
function increment_x!(a::Agent{StdAgent,U},p::Dict)
88
    tdim = length(p["D"])
Victor's avatar
Victor committed
89
    reflected = haskey(p,"reflected") ? p["reflected"] : false
90
    if reflected
Victor Boussange's avatar
Victor Boussange committed
91
        inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
92
        if  tdim > 1
Victor Boussange's avatar
Victor Boussange committed
93
            inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
94
95
        end
    else
Victor Boussange's avatar
Victor Boussange committed
96
97
        # inc = yes no mutation * mutation
        inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
98
99
100
101
102
    end
    a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
 end

 """
103
104
     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"].
105
106
     ONLY FOR GRAPH TYPE DOMAINS
 """
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#  function increment_x!(a::Agent{Mixed},p::Dict)
#      tdim = length(p["D"])
#          inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn(),1,p["nodes"])]
#          if  tdim > 1
#              inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
#          end
#      else
#          # inc = yes no mutation * mutation
#          inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
#      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)
124
125
    Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
"""
126
function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
127
    if x + inc < s
128
        inc = 2 * ( s - x ) - inc
129
    elseif  x + inc > e
130
        inc = 2 * ( e - x ) - inc
131
132
133
    else
        return inc
    end
134
    get_inc_reflected(x,inc,s,e)
135
136
137
138
end

# need to make sure that this is working correctly
"""
139
    α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Victor Boussange's avatar
Victor Boussange committed
140
Generalised gaussian competition kernel
141
"""
142
function α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Victor Boussange's avatar
Victor Boussange committed
143
        return exp( -.5* sum(sum((a1 .- a2).^n_alpha,dims=2)./ sigma_a[:].^n_alpha))
144
end
Victor Boussange's avatar
Victor Boussange committed
145
"""
146
    α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Victor Boussange's avatar
Victor Boussange committed
147
148
Default Gaussian competition kernel
"""
149
α(a1::Array{Number},a2::Array{Number},sigma_a::Array{Number}) = α(a1,a2,2,sigma_a)
150
151

"""
152
    K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
153
154
Gaussian resource kernel
"""
155
function K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
156
157
158
159
    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
    return K0 * pdf(MvNormal(μ,sigma_K),x)
end
"""
160
    K(x::Array{Number},K0::Number,sigma_K::Array{Number})
161
162
Gaussian resource kernel with mean 0.
"""
163
function K(x::Array{Number},K0::Number,sigma_K::Array{Number})
164
165
    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
    return K0 * pdf(MvNormal(sigma_K),x)
166
167
end

168
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)
169
170

"""
171
    function tin(t::Number,a::Number,b::Number)
172
173
174
if t in [a,b) returns 1. else returns 0
"""

175
function tin(t::Number,a::Number,b::Number)
176
177
178
179
180
181
182
183
184
185
    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