To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 20e4134a authored by Victor's avatar Victor
Browse files

Many changes, we are almost send a full run. We need now to modify constructor...

Many changes, we are almost send a full run. We need now to modify constructor of agent, to distinguish between AgentA and AgentNA.
Then we need to continue unit testing with simulation, maybe adding a callback function.
After that, we are good to go.
parent b2e2a80b
......@@ -24,11 +24,14 @@ module ABMEv
export GraphSpace,ContinuousSegment,DiscreteSegment,RealSpace,
AbstractSpacesTuple,get_inc
export update_rates!
export MixedAgent,StdAgent,Agent,get_fitness,get_x,get_t,get_dim,
get_nancestors,get_xarray,get_xhist,
export AbstractAgent,Agent,get_fitness,get_x,get_t,get_dim,
nancestors,get_xarray,get_xhist,
get_thist,get_geo,get_b,get_d,increment_x!,get_inc_reflected,world2df,
split_move,split_merge_move,tin,new_world_G
export copy,runWorld_store_WF,runWorld_store_G,clean_world #,runWorld_G!,runWorld_WF!,
export World,parameters,time,space,agents,size,maxsize,addAgent!,removeAgent!
export runWorld!,give_birth,updateWorld!,update_clock!,updateBirthEvent!,
updateDeathEvent!#,runWorld_G!,runWorld_WF!,
export Simulation,add_entry!
export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div,
get_dist_hist,get_pairwise_average_isolation,
get_local_pairwise_average_isolation,
......
......@@ -5,13 +5,27 @@ hasrates(::Type{Rates{T}}) where {T} = T # not sure we need it
abstract type AbstractAgent{A<:Ancestors,R<:Rates} end # tc for time contingency, fit for fitness coeff
AbstractAgentM = Union{Missing,AbstractAgent}
export AbstractAgentM
"""
$(TYPEDEF)
"""
mutable struct Agent{A<:Ancestors,R<:Rates,T,U,V} <: AbstractAgent{A,R}
abstract type Agent{A<:Ancestors,R<:Rates,T<:Tuple,U,V} <: AbstractAgent{A,R}
mutable struct AgentA{R<:Rates,T<:Tuple,U,V} <: Agent{Ancestors{true},R,T,U,V}
# history of traits for geotraits
x_history::Array{T,1}
# birth time of ancestors
t_history::Array{U,1}
# death rate
d::V
#birth rate
b::V
end
struct AgentNA{R<:Rates,T<:Tuple,U,V} <: Agent{Ancestors{false},R,T,U,V}
# history of traits for geotraits
pos::Array{T,1}
x_history::Array{T,1}
# birth time of ancestors
t_history::Array{U,1}
# death rate
......@@ -20,6 +34,8 @@ mutable struct Agent{A<:Ancestors,R<:Rates,T,U,V} <: AbstractAgent{A,R}
b::V
end
eltype(a::Agent{A,R,T,U,V}) where {A,R,T,U,V} = T
# infers position type and zeros
function initpos(s::S) where {S<:AbstractSpacesTuple}
Eltype = eltype.(s)
......@@ -27,18 +43,23 @@ function initpos(s::S) where {S<:AbstractSpacesTuple}
pos = tuple()
for i in 1:length(Eltype)
if Dims[i] > 1
pos = (pos...,tuple(zeros(Eltype[i],Dims[i])...))
pos = (pos...,Eltype[i](ones(Dims[i])))
else
pos = (pos...,zero(Eltype[i]))
pos = (pos...,one(Eltype[i]))
end
end
Tuple{Eltype...},pos
end
# default initialiser
"""
$(SIGNATURES)
Initialises agent with 0 values everywhere
"""
function Agent(s::S;ancestors=false,rates=false) where {S <: AbstractSpacesTuple}
T,pos = initpos(s)
t = ancestors ? [Float64(.0)] : [nothing]
U = ancestors ? Float64 : Nothing
t = zeros(Float64,1)
U = Float64
d = rates ? Float64(.0) : nothing
b = d
V = rates ? Float64 : Nothing
......@@ -46,25 +67,36 @@ function Agent(s::S;ancestors=false,rates=false) where {S <: AbstractSpacesTupl
end
# here pos is provided
"""
$(SIGNATURES)
Initialises agent with `pos` provided
"""
function Agent(s::S, pos::P;ancestors=false,rates=false) where {P,S <: AbstractSpacesTuple}
T = Tuple{eltype.(s)...}
if !(T <: P)
throw(ArgumentError("Position provided does not match with underlying space"))
T = eltype.(s)
for (i,p) in enumerate(pos)
if typeof(p) !== T[i]
try
p = convert(T[i],p)
catch e
throw(ArgumentError("Position provided does not match with underlying space"))
end
end
end
t = ancestors ? Float64(.0) : [nothing]
U = ancestors ? Float64 : Nothing
t = zeros(Float64,1)
U = Float64
d = rates ? Float64(.0) : nothing
b = d
V = rates ? Float64 : Nothing
Agent{Ancestors{ancestors},Rates{rates},T,U,V}([pos],t,d,b)
@show pos, T
Agent{Ancestors{ancestors},Rates{rates},Tuple{T...},U,V}([pos],t,d,b)
end
# TODO : implement pretty print
import Base.copy
copy(a::Agent{A,R,T,U,V}) where {A,R,T,U,V} = Agent{A,R,T,U,V}(copy(a.pos),copy(a.t_history),copy(a.d),copy(a.b))
copy(m::Missing) = missing
copy(n::Nothing) = nothing
Base.copy(a::Agent{A,R,T,U,V}) where {A,R,T,U,V} = Agent{A,R,T,U,V}(copy(a.x_history),copy(a.t_history),copy(a.d),copy(a.b))
Base.copy(m::Missing) = missing
Base.copy(n::Nothing) = nothing
#####################
###Agent accessors###
......@@ -99,14 +131,14 @@ get_x(a::AbstractAgent,i::Integer) = a.x_history[end][Int(i)]
Get time when agent born.
"""
get_t(a::AbstractAgent) = a.t_history[end]
get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
get_xhist(a::Agent,i::Number) = [a.x_history[t][Int(i)] for t in 1:length(a.xhistory)]
get_xhist(a::Agent) = a.x_history
get_thist(a::Agent) = a.t_history
get_d(a::Agent) = a.d
get_b(a::Agent) = a.b
get_fitness(a::Agent) = a.b - a.d
get_dim(a::Agent) = size(a.x_history,1)
get_nancestors(a::Agent) = size(a.x_history,2)
ndims(a::Agent) = size(a.x_history[end],1)
nancestors(a::Agent) = length(a.x_history)
#####################
###World accessors###
......@@ -144,31 +176,39 @@ function get_xarray(world::Array{T,1},t::Number,geotrait::Bool=false) where {T <
return xarray
end
## Modifiers
"""
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
"""
function increment_x!(a::AbstractAgent{A,R},s::AbstractSpacesTuple,p::Dict{String,Any},t::T) where {A<:Ancestors{true},R,T}
import Base.zero
Base.zero(t::Tuple{Vararg{Union{Number,Tuple{Vararg{Number}}}}}) = [zero.(e) for e in t]
import Base.(+)
(+)(t1::Tuple{Vararg{T,N}},t2::Tuple{Vararg{T,N}}) where {T<:Number,N}= tuple([t1[i] + t2[i] for i in 1:length(t1)]...)
function _get_xinc(a::AbstractAgent,s::AbstractSpacesTuple,p::Dict,t::Number)
@unpack D,mu = p
_x = get_x(a)
inc = similar(_x)
for (i,s) in enumerate(ss)
inc = zero(_x)
for (i,ss) in enumerate(s)
if rand() < mu[i]
inc[i] = get_inc.(_x[i],D[i],s)
inc[i] = get_inc(_x[i],D[i],ss)
end
end
tuple((_x .+ inc)...)
end
## Modifiers
"""
$(SIGNATURES)
This function increments agent by random numbers specified in p
ONLY FOR CONTINUOUS DOMAINS
"""
function increment_x!(a::AbstractAgent{A,R},s::AbstractSpacesTuple,p::Dict,t::T) where {A<:Ancestors{true},R,T}
push!(a.t_history,t)
a.x_history = push!(a.x_history, _x + inc);
a.x_history = push!(a.x_history,_get_xinc(a,s,p,t))
return a
end
function increment_x!(a::AbstractAgent{A,R},s::AbstractSpacesTuple,p::Dict{String,Any},t::T) where {A<:Ancestors{false},R,T}
@unpack D,mu = p
_x = get_x(a)
inc = get_inc.(_x,D,s)
a.t_history = [t]
a.x_history = [_x + inc];
function increment_x!(a::AbstractAgent{A,R},s::AbstractSpacesTuple,p::Dict,t::T) where {A<:Ancestors{false},R,T}
a.t_history = [ t ]
a.x_history = [_get_xinc(a,s,p,t)]
return a
end
"""
......
abstract type AbstractAlg end
# this used to be world
mutable struct Simulation{A<:AbstractAgentM, S<:AbstractSpacesTuple,T<:Number,N<:Number}
worldarray::Array{AbstractAgentM,2}
mutable struct Simulation{A<:AbstractAgentM, S<:AbstractSpacesTuple,T<:Number,F}
agentarray::Vector{AbstractAgentM}
space::S
tspan::Vector{T}
aggregates::Vector{Dict{String,Number}}
cb::F
df_agg
p::Dict{String,Any}
end
function Simulation()
tspan = zeros(1);
worldall = reshape(copy.(world0),N,1);
worldall = hcat(worldall,Array{Missing}(missing,N,1))
# callbacks has to be of the form (names" => String[],"aggregates" => Function)
"""
$(SIGNATURES)
"""
function Simulation(w0::World{A,S,T};cb=(names = String[],agg =nothing)) where {A,S,T}
tspan = zeros(1);
NMax = maxsize(w0)
#agentarray is of size 2 at the beginning
agentarray = copy.(collect(w0.agents))
agentarray = hcat(agentarray,Array{Missing}(missing,NMax,1))
!isnothing(cb.agg) ? df_agg = [Dict(cb.names .=> cb.agg(world))] : df_agg = nothing
Simulation{A,S,T,cb.agg}([agentarray],space(w0),tspan,cb.agg,df_agg,parameters(w0))
end
get_tend(s::Simulation) = s.tspan[end]
get_size(s::Simulation) = length(s.tspan)
function add_entry!(s::Simulation,w::World)
j+=1;sw = size(worldall,2);
# we use <= because we save at the end of the wile loop
if sw <= j
# we double the size of worldall
worldall = hcat(worldall,Array{Missing}(missing,N,sw))
end
worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
push!(tspanarray,t)
"""
$(SIGNATURES)
Add `w` with callbacks `s.cb` to `s`
"""
function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F<:Function}
push!(s.agentarray,copy.(collect(w.agents)))
push!(s.tspan,w.t)
push!(s.df_agg,Dict(s.cb.names .=> s.cb.agg(world)))
end
"""
$(SIGNATURES)
Add `w` to `s`
"""
function add_entry!(s::Simulation{A,S,T,F},w::World) where {A,S,T,F<:Nothing}
push!(s.agentarray,copy.(collect(w.agents)))
push!(s.tspan,w.t)
end
function world2df(world::Array{T,1},geotrait=false) where {T <: Agent}
......
......@@ -9,9 +9,12 @@ abstract type IsFinite{T} end
#ife stands for is finite
abstract type AbstractSpace{Dim,T,I} end
AbstractSpacesTuple = Tuple{Vararg{AbstractSpace}}
import Base:ndims,isfinite,eltype
Base.ndims(x::AbstractSpace{Dim,T,I}) where {Dim,T,I} = Dim
Base.isfinite(x::AbstractSpace{Dim,T,IsFinite{t}}) where {Dim,T,t} = t #not sure we need this
Base.eltype(::AbstractSpace{Dim,T,I}) where {Dim,T,I} = T
Base.eltype(::AbstractSpace{Dim,T,I}) where {Dim,T,I} = Dim > 1 ? Tuple{Vararg{T,Dim}} : T
Base.ndims(ss::AbstractSpacesTuple) = length(ss)
Base.eltype(ss::AbstractSpacesTuple) where {Dim,T,I} = Tuple{eltype.(ss)...}
SpaceType=Union{Nothing, AbstractSpace} # not sure what is this used for
......@@ -22,12 +25,12 @@ struct GraphSpace{T} <: AbstractSpace{1,T,IsFinite{true}}
g::AbstractGraph{T}
end
abstract type AbstractSegment{T} <: AbstractSpace{1,T,IsFinite{true}} end
abstract type AbstractSegment{T<:Number} <: AbstractSpace{1,T,IsFinite{true}} end
"""
$(TYPEDEF)
"""
struct ContinuousSegment{T} <: AbstractSegment{T}
struct ContinuousSegment{T<:AbstractFloat} <: AbstractSegment{T}
s::T
e::T
end
......@@ -35,40 +38,47 @@ end
"""
$(TYPEDEF)
"""
struct DiscreteSegment{T} <: AbstractSegment{T}
struct DiscreteSegment{T<:Integer} <: AbstractSegment{T}
s::T
e::T
end
"""
$(TYPEDEF)
A real pace with dimension N and type T
"""
struct RealSpace{N,T} <: AbstractSpace{N,T,IsFinite{false}} end
# TODO: find a way to put a type on get_inc
# TODO: there is probably a better way of dealing with get_inc
# TODO:
"""
$(SIGNATURES)
function get_inc(D,s::AbstractSpace{Dim,T,I}) where {Dim,T,I<:IsFinite{false}}
Returns increment
"""
# TODO: find a way to put a type on D in get_inc
function get_inc(D,s::AbstractSpace{Dim,T,I}) where {Dim,T<:AbstractFloat,I<:IsFinite{false}}
if Dim > 1
tuple(T.(return D[:] .* rand(Dim)))
return Tuple(D .* randn(T,Dim))
else
return T(D * rand())
return D * randn(T)
end
end
get_inc(x,D,s::AbstractSpace{Dim,T,I}) where {Dim,T,I<:IsFinite{false}} = get_inc(D,s)
#TODO: there is probably a better way of dealing with those two functions
function get_inc(x,D,s::ContinuousSegment{T}) where {T}
inc = D * rand()
inc = D * randn(T)
return _reflect1D(x,inc,s)
end
function get_inc(x,D,s::DiscreteSegment{T}) where {T}
inc = D * rand()
inc = D * randn()
return round(T,_reflect1D(x,inc,s))
end
function get_inc(x,D,s::GraphSpace{T}) where {T}
niter = round(Int,D*rand())
niter = round(Int,abs(D*randn()))
if niter > 0
return last(randomwalk(s.g,x,niter)) - x
else
......
......@@ -78,7 +78,10 @@ If `p["dt_saving"]` not specified, it returns an array with two columns,
first corresponding to initial conditions and last corresponding to world in the last time step.
>:warning: if you choose `nagents = 1` then nothing is saved until the end of the simulation.
"""
function run(w::World{AbstractAgent{A,R},S,T},g::G;dt_saving=nothing,callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
# function run(w::World{AbstractAgent{A,R},S,T},g::G;dt_saving=nothing,callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
function run(w::World{AbstractAgent{A,R},S,T},g::G,tend::Number;
dt_saving=nothing,
callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
@unpack tend = parameters(w)
n=size(w);
NMax = sizemax(w)
......@@ -103,7 +106,7 @@ function run(w::World{AbstractAgent{A,R},S,T},g::G;dt_saving=nothing,callbacks=n
@info "saving world @ t = $(t)/ $(tend)"
add_entry!(s,w)
end
dt = updateWorld!!(w,g)
dt = updateWorld!(w,g)
i += 1
end
# Saving last time step
......
# this used to be the worldalive
# TODO: do a constructor that ensures the parameters numerics are of the same type as the agents
mutable struct World{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number}
agents::Vector{AbstractAgentM}
space::S
parameters::Dict{String,Any}
parameters::Dict
t::T
end
#constructor
function World(w::Vector{A},s::S,p::Dict,t::T=0.) where {A<:AbstractAgent,S<:AbstractSpacesTuple,T}
if typeof(p["D"]) != eltype(skipmissing(w)[1])
throw(ArgumentError("Diffusion coefficient does not match with underlying space"))
end
ww = vcat(w,repeat([missing],Int(p["NMax"] - length(w))))
World{A,S,T}(ww,s,p,t)
end
parameters(world::World) = world.parameters
time(w::World) = w.time
time(w::World) = w.t
space(w::World) = w.space
# this throws indices that are occupied by agents
_get_idx(world::World) = collect(eachindex(skipmissing(world.agents)))
# this throws agents of an abstract array of size size(world)
import Base:size,getindex
Base.size(world::World) = length(world.agents) - count(ismissing,world.agents)
Base.getindex(w::World,i::Int) = w.agents[_get_idx(w)[i]]
# this throws an iterators of agents in the world
agents(world) = skipmissing(world.agents)
size(world) = size(world.agents) - count(ismissing,world.agents)
maxsize(w::World) = length(w.agents)
_findfreeidx(w::World) = findfirst(findfirst(ismissing,w.agents))
agents(world::World) = skipmissing(world.agents)
maxsize(w::World) = w.parameters["NMax"]
_findfreeidx(w::World) = findfirst(ismissing,w.agents)
addAgent!(w::World,a::AbstractAgent) = begin
idx = _findfreeidx(w)
w.agents[idx] = a
return nothing
end
removeAgent!(w::World,i::Int) = begin
world[i] = missing
w.agents[_get_idx(w)[i]] = missing
return nothing
end
......
# In this file lie the function for Gillepsie algorithm
"""
$(TYPEDEF)
"""
struct Gillepsie <: AbstractAlg end
export Gillepsie
"""
function give_birth(a::Agent,t,p::Dict)
Used for Gillepsie setting
"""
function give_birth(mum_idx::Int,w::World{AbstractAgent{A,R},S,T}) where {A <: Ancestors{true},R,S,T}
function give_birth(mum_idx::Int,w::World)
new_a = copy(w[mum_idx])
increment_x!(new_a,space(w),parameters(w),time(w))
return new_a
......@@ -44,27 +46,27 @@ Returning dt drawn from an exponential distribution with parameter the total rat
# Args
t is for now not used but might be used for dynamic landscape
"""
function updateWorld!(w::World,g::G) where {G <: Gillepsie}
function updateWorld!(w::World{A,S,T},g::G) where {A,S,T,G <: Gillepsie}
# total update world
agents = agents(world)
events_weights = ProbabilityWeights(vcat(get_d.(agents),get_b.(agents)))
alive = agents(w)
events_weights = ProbabilityWeights(vcat(get_d.(alive),get_b.(alive)))
# Total rate of events
U = sum(events_weights)
dt = - log(rand(T))/U
update_clock!(world,dt)
= sum(events_weights)
dt = - log(rand(T))/
update_clock!(w,dt)
if dt > 0.
i_event = sample(events_weights)
# This is ok since size is a multiple of 2
n = size(world)
n = size(w)
if i_event <= n
# DEATH EVENT
# In this case i_event is also the index of the individual to die in the world_alive
updateDeathEvent!(world,G,i_event)
updateDeathEvent!(w,g,i_event)
else
# birth event
# i_event - n is also the index of the individual to give birth in the world_alive
mum = world[i_event-n]
update_afterbirth_std!(world,idx_offspring)
mum = w[i_event-n]
updateBirthEvent!(w,g,idx_offspring)
end
return dt
else
......@@ -78,17 +80,3 @@ Get rid of missing value in `world`
"""
#TODO : put some type specs here
clean_world(world) = collect(skipmissing(world))
"""
function new_world_G(nagents::Int,p::Dict; spread = 1., offset = 0.)
Returns an array `world0` of size `p["NMax"]`, with `nagents` ancestors agents.
Those agents have traits distributed according to the normal distribution with mean `offset` and standard deviation `spread`.
The dimension of the domain is determined by the size of the array `spread`.
"""
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
using LightGraphs
using Test
using Revise,ABMEv
using UnPack
myspace = (GraphSpace(SimpleGraph(10,10)),RealSpace{1,Float64}())
myagents = [Agent(myspace,ancestors=true,rates=true) for i in 1:10]
d(X,Y) = gaussian(X[1],Y[1],1)
b(X,Y) = gaussian(X[1],0,1)
D = (Int64(1),Float64(1.))
mu = [1,1]
NMax = 100
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
w = World(myagents,myspace,p,0.)
s = Simulation(w)
......@@ -5,38 +5,59 @@ using Revise,ABMEv
##### SPACES #####
mysegment = DiscreteSegment(1,10)
mygraph = GraphSpace(SimpleGraph(10,10))
myline = RealSpace{2,Float64}()
mycontinuoussegment = ContinuousSegment(-1,1)
myspace = (mysegment,mygraph,myline,mycontinuoussegment)
myspace2 = (mysegment,mygraph,myline)
real2d = RealSpace{2,Float64}()
myline = RealSpace{1,Float16}()
mycontinuoussegment = ContinuousSegment(-1.,1.)
myspace = (mysegment,mygraph,myline)
myspace2 = (mysegment,mycontinuoussegment,real2d)
# checking essential properties of spaces
@test isfinite(mysegment) true
@test isfinite(mygraph) true
@test isfinite(myline) false
@test ndims(myline) == 2
@test ndims(real2d) 2
@test isfinite(mycontinuoussegment) true
@test typeof(myspace) <: AbstractSpacesTuple
@test eltype(myspace2) == Tuple{Int64,Float64,Tuple{Float64,Float64}}
# still does not work
@test get_inc([0.],myline) (0.,0.)
@test get_inc([0.,0.],myline) (0.,0.)
# increment on infinite spaces
@test get_inc(0.,myline) (0.)
@test !(get_inc(1.,myline) 0.)
@test !(get_inc(1,1.,myline) 0.)
@test typeof(get_inc([1.,0.],real2d)) == Tuple{Float64,Float64}
@test typeof(get_inc([1.,0.],[1.,0.],real2d)) == Tuple{Float64,Float64}
@test typeof(get_inc([1.],real2d)) == Tuple{Float64,Float64}
@test typeof(get_inc(1.,real2d)) == Tuple{Float64,Float64}
get_inc([1.],real2d)
ABMEv.initpos(myspace2)
# increment on finite spaces
# checking if reflection works
@test mysegment.s < get_inc(5.,100.,mysegment) + 5. < mysegment.e
@test mycontinuoussegment.s < get_inc(0.,100.,mycontinuoussegment) < mycontinuoussegment.e
#checking if graph works
@test get_inc(1,10,mygraph) + 1 vertices(mygraph.g)
@test prod([get_inc(1,10,mygraph) + 1 vertices(mygraph.g) for i in 1:30])
##### AGENTS #######
a1 = Agent(myspace)
a1bis = Agent(myspace2)
a2 = Agent(myspace,ancestors = true)
a3 = Agent(myspace,(1,1,1.))
a4 = Agent(myspace,(1,1,1.),rates=true)
typeof(a1) <: AbstractAgent
a4 = Agent(myspace2,(1,1,(1.,1)),rates=true)
a5 = Agent(myspace2,ancestors=true)
agents = [a1,a2,a3,a4]
agentsm = [a1,a2,a3,a4,missing]
# basic test
@test typeof(a1) <: AbstractAgent
@test eltype(a1) == eltype(myspace)
@test eltype(a5) == eltype(myspace2)
@test typeof(a1) <: AbstractAgent