ABMEv_Agent.jl 7.62 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
    # birth time of ancestors
8
    t_history::Array{Float64,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_x(a::Agent) = a.x_history[:,end]
42
43
44
45
46
47
48
49
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]
50
get_t(a::Agent) = a.t_history[end]
51
52
get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
get_xhist(a::Agent) = a.x_history
53
get_thist(a::Agent) = a.t_history
54
55
56
get_d(a::Agent) = a.d
get_b(a::Agent) = a.b
get_fitness(a::Agent) = a.b - a.d
57
58
get_dim(a::Agent) = size(a.x_history,1)
get_nancestors(a::Agent) = size(a.x_history,2)
59
60

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"))
61
"""
62
    get_x(world::Array{T},t::Number,trait::Integer) where {T <: Agent}
63
Returns trait of every agents of world in the form of an array which dimensions corresponds to the input.
64
If trait = 0 , we return the geotrait.
65

66
"""
67
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
"""
    get_xarray(world::Array{Agent,1})
Returns every traits of every agents of world in the form of an array
"""
73
74
75
76
77
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}
78
79
    xarray = hcat(get_x.(world)...)
    if geotrait
80
        xarray = vcat( xarray, get_geo.(world,t)')
81
    end
Victor's avatar
Victor committed
82
    return xarray
83
84
end

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
# """
#     get_xhist(world::Vector{Agent},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{T}) 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]';
#     end
#     return xhist
# end
105

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# 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


121
function world2df(world::Array{T,1},geotrait=false) where {T <: Agent}
122
123
124
125
126
127
128
129
130
131
132
    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

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


152

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

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

190
191
192
193


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

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

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