Now the data structure is: add belong the rows with a buffer system that

can be set by the user. Changes of SIR and ER.

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

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
 

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
 

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
 

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)