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 5b94cf5c authored by Victor's avatar Victor
Browse files

added many modifications

parent eea8f623
......@@ -140,34 +140,31 @@ function get_xarray(world::Array{T,1},t::Number,geotrait::Bool=false) where {T <
return xarray
end
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[:g] = get_geo.(world)
## 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}
@unpack D,mu = p
_x = get_x(a)
inc = similar(_x)
for (i,s) in enumerate(ss)
if rand() < mu[i]
inc[i] = get_inc.(_x[i],D[i],s)
end
end
return dfw
push!(a.t_history,t)
a.x_history = push!(a.x_history, _x + inc);
end
"""
world2df(world::Array{T,1},t::Number,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},t::Number,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[:g] = get_geo.(world,t)
end
return dfw
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];
end
"""
......
# this used to be world
mutable struct Simulation{A<:AbstractAgentM, S<:AbstractSpacesTuple,T<:Number,N<:Number}
worldarray::Array{AbstractAgentM,2}
tspan::Vector{T}
aggregates::Vector{Dict{String,Number}}
p::Dict{String,Any}
end
function Simulation()
tspan = zeros(1);
worldall = reshape(copy.(world0),N,1);
worldall = hcat(worldall,Array{Missing}(missing,N,1))
end
get_tend(s::Simulation) = s.tspan[end]
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)
end
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[:g] = get_geo.(world)
end
return dfw
end
"""
world2df(world::Array{T,1},t::Number,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},t::Number,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[:g] = get_geo.(world,t)
end
return dfw
end
......@@ -6,52 +6,60 @@ when space is finite
"""
abstract type IsFinite{T} end
Base.isfinite(::Type{IsFinite{T}}) where {T} = T
#ife stands for is finite
abstract type AbstractSpace{Dim,T,I} end
AbstractSpacesTuple = Tuple{Vararg{AbstractSpace}}
Base.ndims(x::AbstractSpace{Dim,T,I}) where {Dim,T,I} = Dim
Base.isfinite(x::AbstractSpace{Dim,T,I}) where {Dim,T,F,I} = isfinite(I) #not sure we need this
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
SpaceType=Union{Nothing, AbstractSpace} # not sure what is this used for
abstract type AbstractDiscreteSpace{Dim,T,I} <: AbstractSpace{Dim,T,I} end
abstract type AbstractGraphSpace{T} <: AbstractDiscreteSpace{1,T,IsFinite{true}} end
struct GraphSpace{T} <: AbstractGraphSpace{T}
struct GraphSpace{T} <: AbstractSpace{1,T,IsFinite{true}}
g::AbstractGraph{T}
end
struct DiscreteSegment{T} <: AbstractDiscreteSpace{1,T,IsFinite{true}}
abstract type AbstractSegment{T} <: AbstractSpace{1,T,IsFinite{true}} end
struct ContinuousSegment{T} <: AbstractSegment{T}
s::T
e::T
end
struct Real1DSpace{T} <: AbstractSpace{1,T,IsFinite{false}} end
struct DiscreteSegment{T} <: AbstractSegment{T}
s::T
e::T
end
struct Real1DSpace{T} <: AbstractSpace{1,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
"""
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::Agent{StdAgent,U},t,p::Dict) where U
tdim = length(p["D"])
reflected = haskey(p,"reflected") ? p["reflected"] : false
if reflected
inc = [get_inc_reflected(get_x(a,1),p["D"][1] *randn())]
if tdim > 1
inc = vcat(inc,(rand(tdim-1) < p["mu"][2:end]) .* p["D"][2:end] .* randn(tdim-1))
end
function get_inc(D,s::AbstractSpace{Dim,T,I}) where {Dim,T,I<:IsFinite{false}}
if Dim > 1
return D[:] .* rand(T,Dim)
else
# inc = yes no mutation * mutation
inc = (rand(tdim) < vec(p["mu"])) .* vec(p["D"][:]) .* randn(tdim)
return D * rand(T)
end
a.x_history = hcat(a.x_history, get_x(a) + reshape(inc,:,1));
push!(a.t_history,t)
end
end
get_inc(x,D,s::AbstractSpace{Dim,T,I}) where {Dim,T,I<:IsFinite{false}} = get_inc(D,s)
function get_inc(x,D,s::ContinuousSegment{T}) where {T}
inc = D * rand(T)
return reflect1D(x,inc,s)
end
function get_inc(x,D,s::DiscreteSegment{T}) where {T}
inc = D * rand(T)
return round(reflect1D(x,inc,s))
end
function get_inc(x,D,s::GraphSpace{T}) where {T}
niter = round(rand(T) < D)
inc = D * rand(T)
return round(reflect1D(x,inc,s))
end
"""
function increment_x!(a::Agent{MixedAgent,U},t::U,p::Dict) where U
......@@ -73,13 +81,13 @@ end
function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
Here we increment the trajectory of trait 1 such that it follows a reflected brownian motion (1D)
"""
function get_inc_reflected(x::Number,inc::Number,s=-1,e=1)
function reflect1D(x::Number,inc::Number,s<:AbstractSegment)
if x + inc < s
inc = 2 * ( s - x ) - inc
inc = 2 * ( s.s - x ) - inc
elseif x + inc > e
inc = 2 * ( e - x ) - inc
inc = 2 * ( s.e - x ) - inc
else
return inc
end
get_inc_reflected(x,inc,s,e)
reflect1D(x,inc,s)
end
struct AgentBasedModel{A, S, F, P}
agents::A
space::S
scheduler::F
properties::P
end
function AgentBasedModel(a::Vector{A},s::S) where {A<:AbstractAgentM, S<:AbstractSpacesTuple}
AgentBasedModel(a,s,nothing, nothing)
end
function AgentBasedModel(N <:Int, Max <: Int ,s::S) where {A<:AbstractAgentM, S<:AbstractSpacesTuple}
# To be implemented
end
# for now we consider that competition is local within an array
abstract type AbstractAlg end
function solve(m::Model;dt_saving=nothing,callbacks=nothing)
end
"""
update_rates_std!(world,p::Dict,t::Float64)
......@@ -6,24 +10,25 @@ This standard updates takes
- competition kernels of the form α(x,y) and
- carrying capacity of the form K(x)
"""
function update_rates_std!(world,p::Dict,t::Float64)
α = p["alpha"];K=p["K"];
traits = get_x.(world)
function update_rates_std!(w::World)
@unpack b,k = parameters(w)
_agents = agents(w)
traits = get_x.(_agents)
# traits = get_xhist.(world)
N = length(traits)
D = zeros(N)
n = size(world)
D = zeros(n)
# Here you should do a shared array to compute in parallel
for i in 1:(N-1)
for j in i+1:N
C = α(traits[i],traits[j])
for i in 1:(n-1)
for j in i+1:n
C = d(traits[i],traits[j])
D[i] += C
D[j] += C
end
end
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
for (i,a) in _agents
a.d = D[i]
a.b = K(traits[i])
a.b = b(traits[i])
end
end
......@@ -39,151 +44,35 @@ function update_rates_graph!(world,C,p::Dict,t::Float64)
end
end
function update_rates_2D!(world,C,p::Dict,t::Float64)
# Competition matrix
traits = get_x.(world)
# traits = get_xhist.(world)
N = length(traits)
# C = SharedArray{Float64}((N,N))
# Here you should do a shared array to compute in parallel
for i in 1:(N-1)
C[i,i] = 1.
for j in i+1:N
# be careful, for xhist what is return is an array hence no need of []
C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
C[j,i] = C[i,j]
end
end
C[N,N] = 1.
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
# we only update death rate
# a.d = sum(C[i,:]) / K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1])
a.d = sum(C[i,:])
# /!| not ideal to assign at every time step the birth rate that is constant
a.b = KK(traits[i][1:1],
p["K0"],
p["n_K"],
p["sigma_K"][1:1],
split_merge_move(t),
-split_merge_move(t))/K(traits[i][2:2],
p["K0"],
p["n_K"],
p["sigma_K"][2:2])
end
end
function update_rates_grad2D!(world,C,p::Dict,t::Float64)
# Competition matrix
traits = get_x.(world)
# traits = get_xhist.(world)
N = length(traits)
# C = SharedArray{Float64}((N,N))
# Here you should do a shared array to compute in parallel
for i in 1:(N-1)
C[i,i] = 1.
for j in i+1:N
# be careful, for xhist what is return is an array hence no need of []
C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
C[j,i] = C[i,j]
end
end
C[N,N] = 1.
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
# we only update death rate
a.d = sum(C[i,:])
a.b = K(traits[i][2:2],p["K0"],p["n_K"],p["sigma_K"], μ = p["a"]*traits[i][1])
end
end
function update_rates_mountain!(world,C,p::Dict,t::Float64)
# Competition matrix
traits = get_x.(world)
# traits = get_xhist.(world)
N = length(traits)
# C = SharedArray{Float64}((N,N))
# Here you should do a shared array to compute in parallel
for i in 1:(N-1)
C[i,i] = 1.
for j in i+1:N
# be careful, for xhist what is return is an array hence no need of []
C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
C[j,i] = C[i,j]
end
end
C[N,N] = 1.
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
# we only update death rate
a.d = sum(C[i,:])
a.b = KK(traits[i][1:1],
p["K0"],
p["n_K"],
p["sigma_K"][1:1],
split_move(t),
-split_move(t))/K(traits[i][2:2],
p["K0"],
p["n_K"],
p["sigma_K"][2:2],
μ=(tin(traits[i][1],-1.,0.)*(traits[i][1]+1) - tin(traits[i][1],0.,1.)*(traits[i][1]-1) ) * split_move(t))
end
end
function update_rates_std_split!(world,C,p::Dict,t::Float64)
# Competition matrix
traits = get_x.(world)
# traits = get_xhist.(world)
N = length(traits)
# C = SharedArray{Float64}((N,N))
# Here you should do a shared array to compute in parallel
for i in 1:(N-1)
C[i,i] = 1.
for j in i+1:N
C[i,j] = α(traits[i],traits[j],p["n_alpha"],p["sigma_a"])
C[j,i] = C[i,j]
end
end
C[N,N] = 1.
# Here we can do it in parallel as well
for (i,a) in enumerate(world)
# we only update death rate
# this could be imptrove since \alpha is symmetric, by using a symmetric matrix
# a.d = sum(C[i,:]) / KK(traits[i][:,end],p["K0"],p["n_K"],p["sigma_K"],split_merge_move(t),-split_merge_move(t))
a.d = sum(C[i,:])
# /!| not ideal to assign at every time step the birth rate that is constant
a.b = KK(traits[i][:,end],p["K0"],p["n_K"],p["sigma_K"],split_move(t),-split_move(t))
end
end
"""
function runWorld_store_WF(p,world0;mode="std")
Wright Fisher process. Returns an array worldall with every agents.
"""
function runWorld_store_WF(p,world0;mode="std")
worldall = repeat(world0,inner = (1,length(1:Int(p["tend"]))))
N=length(world0);
newworld = copy.(world0)
if mode == "std"
update_rates! = update_rates_std!
elseif mode == "2D"
update_rates! = update_rates_2D!
elseif mode == "grad2D"
update_rates! = update_rates_grad2D!
elseif mode == "mountain"
update_rates! = update_rates_mountain!
elseif mode == "split"
update_rates! = update_rates_std_split!
elseif mode == "graph"
update_rates! = update_rates_graph!
else
error("Mode $mode is not recognized")
end
for i in 1:(Int(p["tend"])-1)
# we have to take views, otherwise it does not affect worldall
world = @view worldall[:,i];newworld = @view worldall[:,i+1];
updateWorld_WF!(world,newworld,p,update_rates!,Float64(i))
end
return worldall,collect(0:Int(p["tend"]-1))
# worldall = repeat(world0,inner = (1,length(1:Int(p["tend"]))))
# N=length(world0);
# newworld = copy.(world0)
# if mode == "std"
# update_rates! = update_rates_std!
# elseif mode == "2D"
# update_rates! = update_rates_2D!
# elseif mode == "grad2D"
# update_rates! = update_rates_grad2D!
# elseif mode == "mountain"
# update_rates! = update_rates_mountain!
# elseif mode == "split"
# update_rates! = update_rates_std_split!
# elseif mode == "graph"
# update_rates! = update_rates_graph!
# else
# error("Mode $mode is not recognized")
# end
# for i in 1:(Int(p["tend"])-1)
# # we have to take views, otherwise it does not affect worldall
# world = @view worldall[:,i];newworld = @view worldall[:,i+1];
# updateWorld_WF!(world,newworld,p,update_rates!,Float64(i))
# end
# return worldall,collect(0:Int(p["tend"]-1))
end
"""
function runWorld_store_G(p,world0)
......@@ -194,48 +83,36 @@ 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 runWorld_store_G(p,world0)
# we store the value of world every 100 changes by default
function run(w::World{AbstractAgent{A,R},S,T},g::G;dt_saving=nothing,callbacks=nothing) where {G<:Gillepsie,A,R,S,T}
@unpack tend = parameters(w)
n=size(w);
NMax = sizemax(w)
t = .0
i = 1;j=1;dt = 1.
N=length(world0);
tspanarray = zeros(1);
dt_saving = haskey(p,"dt_saving") ? p["dt_saving"] : p["tend"] + 1.
# Array that stores the agents
worldall = reshape(copy.(world0),N,1);
# We add a second row for final time step
worldall = hcat(worldall,Array{Missing}(missing,N,1))
# we instantiate C as the biggest size it can take
update_rates_std!(skipmissing(world0),p,0.)
while t<p["tend"]
isnothing(dt_saving) ? dt_saving = p["tend"] + 1. : nothing
sim = Simulation(w,callbacks=callbacks)
if R <: Rates{true}
update_rates!(w)
end
while t<tend
if dt < 0
throw("We obtained negative time step dt = $dt at event $i")
elseif count(ismissing,world0) == p["NMax"]
elseif size(world) == NMax
@info "All individuals have died :("
break
elseif count(ismissing,world0) == 0
elseif size(world) == 0
@info "We have reached the maximum number of individuals allowed"
break
end
if t - tspanarray[end] >= dt_saving
@info "saving world @ t = $(t)/ $(p["tend"])"
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)
if t - get_tend(sim) >= dt_saving
@info "saving world @ t = $(t)/ $(tend)"
add_entry!(s,w)
end
dt = updateWorld_G!(world0,p,update_rates_std!,t)
t += dt
dt = updateWorld!!(w,g)
i += 1
end
# Saving last time step
j+=1
worldall[1:Int(N - count(ismissing,world0)),j] .= copy.(collect(skipmissing(world0)));
push!(tspanarray,t)
add_entry!(sim,w)
@info "simulation stopped at t=$(tspanarray[end]), after $(i) generations"
return worldall[:,1:j],tspanarray
return sim
end
# this used to be the worldalive
mutable struct World{A<:AbstractAgent, S<:AbstractSpacesTuple,T<:Number}
agents::Vector{AbstractAgentM}
space::S
parameters::Dict{String,Any}
t::T
end
parameters(world::World) = world.parameters
time(w::World) = w.time
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)
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))
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
return nothing
end
update_clock!(w::World{A,S,T},dt) where {A,S,T} = begin
w.t = convert(T,sum(w.t + dt))
return nothing
end
# In this file lie the function for Gillepsie algorithm
struct Gillepsie <: AbstractAlg end
"""
function give_birth(a::Agent,t,p::Dict)
Used for Gillepsie setting
"""
function give_birth(a::Agent,t,p::Dict)
new_a = copy(a)
increment_x!(new_a,t,p)
function give_birth(mum_idx::Int,w::World{AbstractAgent{A,R},S,T}) where {A <: Ancestors{true}}
new_a = copy(w[mum_idx])
increment_x!(new_a,space(w),parameters(w),time(w))
return new_a
end
function update_afterbirth_std!(world,idx_offspring,p::Dict) where T
function updateBirthEvent!(world::World,::Gillepsie,mum_idx::Int)
# updating competition only the two columns corresponding to agent idx
@unpack d,b = p
x_offspring = get_x(world[idx_offspring])
for a in skipmissing(world)
@unpack d,b = parameters(world)
offspring = give_birth(w,mum_idx)
x_offspring = get_x(offspring)
_agents = agents(world)
for a in _agents
a.d += d(get_x(a),x_offspring)
end
# Now updating new agent
world[idx_offspring].d = sum(d.(get_x.(skipmissing(world)),Ref(x_offspring))) - d(x_offspring,x_offspring)
world[idx_offspring].b = b(x_offspring)
offspring.d = sum(d.(get_x.(_agents),Ref(x_offspring))) - d(x_offspring,x_offspring)