ABMEv_Agent.jl 7.22 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
# This  constructor should be used when one wants to impose the type of the agent (e.g. Mixed)
15
Agent{T}(xhist::Array{U}) where {T,U} = Agent{T,U}(reshape(xhist,:,1),0.,1.)
16
17

# This constructor is used by default
18
Agent(xhist::Array{U}) where {U <: Number} = Agent{StdAgent}(xhist)
19

20
Agent() = Agent(Float64[],0.,1.)
21
import Base.copy
22
copy(a::Agent{T,U}) where {T,U} = Agent{T,U}(a.x_history,a.d,a.b)
23
24
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
Mainly works for WF-type world
52
53
54
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.

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

58
59
60
61
62
63
64
65
66
67
68
69
"""
    get_xarray(world::Array{Agent,1})
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}
    xarray = hcat(get_x.(world)...)
    if geotrait
        xarray = vcat( xarray, get_geo.(world)')
    end
end

70
"""
71
    get_xhist(world::Vector{Agent},geotrait = false)
72
73
74
75
76
77
78
79
80
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
"""
Victor's avatar
Victor committed
81
function get_xhist(world::Vector{T},geotrait = false) where {T <: Agent}
82
    hist = minimum(get_nancestors.(world))
Victor's avatar
Victor committed
83
84
85
    ntraits = get_dim(first(world));
    xhist = zeros(length(world), hist, ntraits + geotrait);
    for (i,a) in enumerate(world)
86
87
88
89
90
91
92
93
94
        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


95
"""
96
97
increment_x!(a::Agent{StdAgent,U},p::Dict)
    This function increments agent by random numbers specified in p
98
99
    ONLY FOR CONTINUOUS DOMAINS
"""
100
function increment_x!(a::Agent{StdAgent,U},p::Dict) where U
101
    tdim = length(p["D"])
Victor's avatar
Victor committed
102
    reflected = haskey(p,"reflected") ? p["reflected"] : false
103
    if reflected
Victor Boussange's avatar
Victor Boussange committed
104
        inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
105
        if  tdim > 1
Victor Boussange's avatar
Victor Boussange committed
106
            inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
107
108
        end
    else
Victor Boussange's avatar
Victor Boussange committed
109
110
        # inc = yes no mutation * mutation
        inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
111
112
113
114
115
    end
    a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
 end

 """
116
117
     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"].
118
119
     ONLY FOR GRAPH TYPE DOMAINS
 """
120
121
122
123
124
125
126
127
128
 function increment_x!(a::Agent{MixedAgent,U},p::Dict) where U
     tdim = length(p["D"])
     inc = [round(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
     a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
end

129
130
131
132


"""
get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
133
134
    Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
"""
135
function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
136
    if x + inc < s
137
        inc = 2 * ( s - x ) - inc
138
    elseif  x + inc > e
139
        inc = 2 * ( e - x ) - inc
140
141
142
    else
        return inc
    end
143
    get_inc_reflected(x,inc,s,e)
144
145
146
147
end

# need to make sure that this is working correctly
"""
148
    α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Victor Boussange's avatar
Victor Boussange committed
149
Generalised gaussian competition kernel
150
"""
151
function α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Victor Boussange's avatar
Victor Boussange committed
152
        return exp( -.5* sum(sum((a1 .- a2).^n_alpha,dims=2)./ sigma_a[:].^n_alpha))
153
end
Victor Boussange's avatar
Victor Boussange committed
154
"""
155
    α(a1::Array{Number},a2::Array{Number},n_alpha::Number,sigma_a::Array{Number})
Victor Boussange's avatar
Victor Boussange committed
156
157
Default Gaussian competition kernel
"""
158
α(a1::Array{Number},a2::Array{Number},sigma_a::Array{Number}) = α(a1,a2,2,sigma_a)
159
160

"""
161
    K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
162
163
Gaussian resource kernel
"""
164
function K(x::Array{Number},K0::Number,μ::Array{Number},sigma_K::Array{Number})
165
166
167
168
    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
    return K0 * pdf(MvNormal(μ,sigma_K),x)
end
"""
169
    K(x::Array{Number},K0::Number,sigma_K::Array{Number})
170
171
Gaussian resource kernel with mean 0.
"""
172
function K(x::Array{Number},K0::Number,sigma_K::Array{Number})
173
174
    # return K0*exp(-sum(sum((x .- μ).^n_K,dims=2)./sigma_K[:].^n_K))
    return K0 * pdf(MvNormal(sigma_K),x)
175
176
end

177
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)
178
179

"""
180
    function tin(t::Number,a::Number,b::Number)
181
182
183
if t in [a,b) returns 1. else returns 0
"""

184
function tin(t::Number,a::Number,b::Number)
185
186
187
188
189
190
191
192
193
194
    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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212

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