ABMEv_Agent.jl 7.65 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
    # birth time of ancestors
    t_history::Array{U,1}
9
10
11
12
13
14
15
    # death rate
    d::Float64
    #birth rate
    b::Float64
end

# Constructors
16
# This  constructor should be used when one wants to impose the type of the agent (e.g. Mixed)
17
Agent{T}(xhist::Array{U}) where {T,U} = Agent{T,U}(reshape(xhist,:,1),[0.],0.,1.)
18
19

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

22
Agent() = Agent(Float64[],0.,0.,1.)
23
import Base.copy
24
copy(a::Agent{T,U}) where {T,U} = Agent{T,U}(copy(a.x_history),copy(a.t_history),copy(a.d),copy(a.b))
25
26
copy(m::Missing) = missing

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

# returns trait i of the agent
41
get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
42
get_xhist(a::Agent) = a.x_history
43
get_x(a::Agent) = a.x_history[:,end]
44
45
46
47
48
49
50
51
function get_geo(a::Agent{U,T},t::Number) where {U,T}
    tarray = vcat(a.t_history[2:end],convert(T,t))
    tarray .-= a.t_history
    return sum(get_xhist(a,1) .* tarray)
end
# This method can acces geotrait, while the second not
get_x(a::Agent,t::Number,i::Integer) = i > 0 ? a.x_history[Int(i),end] : get_geo(a,t)
get_x(a::Agent,i::Integer) = a.x_history[Int(i),end]
52
53
54
get_d(a::Agent) = a.d
get_b(a::Agent) = a.b
get_fitness(a::Agent) = a.b - a.d
55
56
get_dim(a::Agent) = size(a.x_history,1)
get_nancestors(a::Agent) = size(a.x_history,2)
57
"""
58
    get_xarray(world::Array{Agent},trait::Int)
59
Mainly works for WF-type world
60
Returns trait of every agents of world in the form of an array which dimensions corresponds to the input.
61
If trait = 0 , we return the geotrait.
62
63
Particularly suited for an array world corresponding to a timeseries.

64
"""
65
66
67
get_x(world::Array{T},trait::Integer) where {T <: Agent} = trait > 0 ? reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2)) : throw(ErrorException("Not the right method, need `t` as an argument"))

get_x(world::Array{T},t::Number,trait::Integer) where {T <: Agent} = trait > 0 ? reshape(hcat(get_x.(world,trait)),size(world,1),size(world,2)) : reshape(hcat(get_geo.(world,t)),size(world,1),size(world,2))
68

69
70
71
72
73
"""
    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.
"""
74
75
76
77
78
function get_xarray(world::Array{T,1}) where {T <: Agent}
    return hcat(get_x.(world)...)
end

function get_xarray(world::Array{T,1},t::Number,geotrait::Bool=false) where {T <: Agent}
79
80
    xarray = hcat(get_x.(world)...)
    if geotrait
81
        xarray = vcat( xarray, get_geo.(world,t)')
82
    end
Victor's avatar
Victor committed
83
    return xarray
84
85
end

86
"""
87
    get_xhist(world::Vector{Agent},geotrait = false)
88
89
90
91
92
93
94
95
96
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
"""
97
function get_xhist(world::Vector{T}) where {T <: Agent}
98
    hist = minimum(get_nancestors.(world))
Victor's avatar
Victor committed
99
100
101
    ntraits = get_dim(first(world));
    xhist = zeros(length(world), hist, ntraits + geotrait);
    for (i,a) in enumerate(world)
102
103
104
105
106
        xhist[i,:,1:end-geotrait] = get_xhist(a)[:,end-hist+1:end]';
    end
    return xhist
end

107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# TODO: This method broken, when one ask for the geotraits
# function get_xhist(world::Vector{T},t::Number,geotrait = false) where {T <: Agent}
#     hist = minimum(get_nancestors.(world))
#     ntraits = get_dim(first(world));
#     xhist = zeros(length(world), hist, ntraits + geotrait);
#     for (i,a) in enumerate(world)
#         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


function world2df(world::Array{T,1}) 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[:g] = get_geo.(world)
    end
    return dfw
end

134
"""
135
    world2df(world::Array{T,1},t::Number,geotrait = false) where {T <: Agent}
136
137
138
139
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
"""
140
function world2df(world::Array{T,1},t::Number,geotrait = false) where {T <: Agent}
141
142
143
144
145
146
    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
147
        dfw[:g] = get_geo.(world,t)
148
149
150
151
152
    end
    return dfw
end


153

154
"""
155
increment_x!(a::Agent{StdAgent,U},t::U,p::Dict) where U
156
    This function increments agent by random numbers specified in p
157
158
    ONLY FOR CONTINUOUS DOMAINS
"""
159
function increment_x!(a::Agent{StdAgent,U},t,p::Dict) where U
160
    tdim = length(p["D"])
Victor's avatar
Victor committed
161
    reflected = haskey(p,"reflected") ? p["reflected"] : false
162
    if reflected
Victor Boussange's avatar
Victor Boussange committed
163
        inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
164
        if  tdim > 1
Victor Boussange's avatar
Victor Boussange committed
165
            inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
166
167
        end
    else
Victor Boussange's avatar
Victor Boussange committed
168
169
        # inc = yes no mutation * mutation
        inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
170
171
    end
    a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
172
    push!(a.t_history,t)
173
174
175
 end

 """
176
     function increment_x!(a::Agent{MixedAgent,U},t::U,p::Dict) where U
177
 This function increments first trait of agent with integer values, that are automatically reflected between 1 and p["nodes"].
178
179
Other traits are incremented as usual.
TODO : make it work for a graph type landscape, where domain is not a line anymore.
180
 """
181
 function increment_x!(a::Agent{MixedAgent,U},t,p::Dict) where U
182
183
184
185
186
187
     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));
188
     push!(a.t_history,t)
189
190
end

191
192
193
194


"""
get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
195
196
    Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
"""
197
function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
198
    if x + inc < s
199
        inc = 2 * ( s - x ) - inc
200
    elseif  x + inc > e
201
        inc = 2 * ( e - x ) - inc
202
203
204
    else
        return inc
    end
205
    get_inc_reflected(x,inc,s,e)
206
207
208
end

"""
209
    function tin(t::Number,a::Number,b::Number)
210
211
212
if t in [a,b) returns 1. else returns 0
"""

213
function tin(t::Number,a::Number,b::Number)
214
215
216
217
218
219
220
221
222
223
    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