diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 23c1828641e4bf3a3d0f169afc5f910c7c747e24..cef5b9048a8811114d3c24b5ae5d3a8b8af5cfe7 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -9,7 +9,8 @@ export load_model, get_module_path
include("model.jl")
export Observations, AbstractTrajectory, Trajectory
-export +,-,δ,get_obs_variables,get_states_number
+export +,-,δ
+export get_obs_variables,get_states_number
include("observations.jl")
end
diff --git a/core/model.jl b/core/model.jl
index a1b584c8ebe28da198a624cf7a27aae8c1cb0d39..8aa7b2ee4ce1c35da59868397a5faf64dcef0cfa 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -20,11 +20,13 @@ mutable struct CTMC <: ContinuousTimeModel
_g_idx::Vector{Int} # of dimension dobs
is_absorbing::Function
time_bound::Float64
+ buffer_size::Int
end
function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::Vector{String},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
- f!::Function, is_absorbing::Function; g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf)
+ f!::Function, is_absorbing::Function;
+ g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10)
dobs = length(g)
_map_obs_var_idx = Dict()
_g_idx = Vector{Int}(undef, dobs)
@@ -40,41 +42,46 @@ function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_tra
@warn "You have possibly redefined a function Model.is_absorbing used in a previously instantiated model."
end
- return CTMC(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound)
+ return CTMC(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
end
function simulate(m::ContinuousTimeModel)
# trajectory fields
- full_values = zeros(m.d, 0)
+ full_values = zeros(0, m.d)
times = zeros(0)
transitions = Vector{Union{String,Nothing}}(undef,0)
# values at time n
n = 0
xn = m.x0
tn = m.t0
- tr = [""]
# at time n+1
- xnplus1 = zeros(Int, m.d)
- tnplus1 = zeros(Float64, 1)
- is_absorbing = (m.is_absorbing(m.p,xn))::Bool
+ mat_x = zeros(Int, m.buffer_size, m.d)
+ l_t = zeros(Float64, m.buffer_size)
+ l_tr = Vector{String}(undef, m.buffer_size)
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
while !is_absorbing &&Â (tn <= m.time_bound)
- m.f!(xnplus1, tnplus1, tr, xn, tn, m.p)
- full_values = hcat(full_values, xnplus1)
- push!(times, tnplus1[1])
- push!(transitions, tr[1])
- xn, tn = xnplus1, tnplus1[1]
- n += 1
+ i = 0
+ while i < m.buffer_size && !is_absorbing &&Â (tn <= m.time_bound)
+ i += 1
+ m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
+ xn = @view mat_x[i,:]
+ tn = l_t[i]
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ full_values = vcat(full_values, @view mat_x[1:i,:])
+ append!(times, @view l_t[1:i])
+ append!(transitions, @view l_tr[1:i])
+ n += i
is_absorbing = m.is_absorbing(m.p,xn)::Bool
end
- values = full_values[m._g_idx,:]
+ values = @view full_values[:,m._g_idx]
if is_bounded(m)
if times[end] > m.time_bound
- values[:, end] = values[:,end-1]
+ values[end,:] = values[end-1,:]
times[end] = m.time_bound
transitions[end] = nothing
end
end
-
return Trajectory(m, values, times, transitions)
end
diff --git a/core/observations.jl b/core/observations.jl
index e82b222085b280a0ae6f910c30424f8a91caff8d..99d2b55fff238285155d297277762f697bd6d464 100644
--- a/core/observations.jl
+++ b/core/observations.jl
@@ -12,21 +12,39 @@ end
function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
function -(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
function δ(σ1::AbstractTrajectory,t::Float64) end
-function get_obs_variables(σ::AbstractTrajectory) end
-get_values(σ::AbstractTrajectory, var::String) =
-σ.values[(σ.m)._map_obs_var_idx[var],:]
+# Properties of the trajectory
+get_states_number(σ::AbstractTrajectory) = length(σ.times)
+get_obs_variables(σ::AbstractTrajectory) = (σ.m).g
-get_states_number(σ::AbstractTrajectory) =
-length(σ.times)
+# Access to trajectory values
+get_var_values(σ::AbstractTrajectory, var::String) =
+@view σ.values[:,(σ.m)._map_obs_var_idx[var]]
+get_state(σ::AbstractTrajectory, idx::Int) = @view σ.values[idx,:]
+get_value(σ::AbstractTrajectory, var::String, idx::Int) =
+σ.values[idx,(σ.m)._map_obs_var_idx[var]]
-function getindex(σ::AbstractTrajectory, idx::String)
- if idx == "times"
+# Get var values ["I"]
+function getindex(σ::AbstractTrajectory, var::String)
+ if var == "times"
return σ.times
- elseif idx == "transitions"
+ elseif var == "transitions"
return σ.transitions
else
- return get_values(σ, idx)
+ return get_var_values(σ, var)
+ end
+end
+# Get i-th state [i]
+getindex(σ::AbstractTrajectory, idx::Int) = get_state(σ, i)
+# Get i-th value of var ["I", idx]
+function getindex(σ::AbstractTrajectory, var_idx::Tuple{String,Int})
+ var, idx = var_idx[1], var_idx[2]
+ if var == "times"
+ return σ.times[idx]
+ elseif var == "transitions"
+ return σ.transitions[idx]
+ else
+ return get_value(σ, var, idx)
end
end
diff --git a/models/ER.jl b/models/ER.jl
index eb0d4676a929a93bc51e6d2eeb73ed1afc3ad5d3..7fa9477ec33fe2a2ff82ac02b190d654aec530b4 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -9,9 +9,8 @@ l_tr = ["R1","R2","R3"]
p = [1.0, 1.0, 1.0]
x0 = [100, 100, 0, 0]
t0 = 0.0
-
-function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
- xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+ xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
a3 = p[3] * xn[3]
@@ -35,18 +34,17 @@ function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
b_sup += l_a[i+1]
end
- nu = @view l_nu[:,reaction]
- xnplus1[1] = xn[1]+nu[1]
- xnplus1[2] = xn[2]+nu[2]
- xnplus1[3] = xn[3]+nu[3]
- xnplus1[4] = xn[4]+nu[4]
- tnplus1[1] = tn + tau
- tr[1] = "R$(reaction)"
+ nu = @view l_nu[:,reaction] # macro for avoiding a copy
+ for i = 1:4
+ mat_x[idx,i] = xn[i]+nu[i]
+ end
+ l_t[idx] = tn + tau
+ l_tr[idx] = "R$(reaction)"
end
-is_absorbing_er(p::Vector{Float64},xn::Vector{Int}) =
+is_absorbing_ER(p::Vector{Float64},xn::AbstractVector{Int}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g = ["P"]
-ER = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_er; g=g)
+ER = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_f!,is_absorbing_ER; g=g)
export ER
diff --git a/models/SIR.jl b/models/SIR.jl
index d448cc2ac141b678786bc7b162a0622b74dbcbcb..6489084d275628e57fe74ac2347ae8ac21e0c829 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -9,8 +9,8 @@ l_tr = ["R1","R2"]
p = [0.0012, 0.05]
x0 = [95, 5, 0]
t0 = 0.0
-function SIR_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
- xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+ xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
l_a = SVector(a1, a2)
@@ -35,13 +35,13 @@ function SIR_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{Strin
end
nu = @view l_nu[:,reaction] # macro for avoiding a copy
- xnplus1[1] = xn[1]+nu[1]
- xnplus1[2] = xn[2]+nu[2]
- xnplus1[3] = xn[3]+nu[3]
- tnplus1[1] = tn + tau
- tr[1] = "R$(reaction)"
+ for i = 1:3
+ mat_x[idx,i] = xn[i]+nu[i]
+ end
+ l_t[idx] = tn + tau
+ l_tr[idx] = "R$(reaction)"
end
-is_absorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+is_absorbing_SIR(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g = ["I"]
SIR = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_f!,is_absorbing_SIR; g=g)