Fix the redefinition of the f! function of models (side effect of

previously loaded model). The package warns when this happen.

Several benchmarks for the data structure. Now we chose the row buffer
method according to results (mainly because the trajectory access with
far way faster, but strangly not that much differences in other cases).

bench/array_memory_order/access_trajectory_sir.jl

+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+bound_time = 200.0
+l_var = ["S", "I", "R"]
+
+println("Col buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+function access_trajectory_col(m::Model)
+    res = 0
+    n_sim = 100
+    for k = 1:n_sim
+        σ = _simulate_col_buffer(m)
+        t = _get_values_col(σ, "I")
+        res += t[end-1]
+    end
+    return res
+end
+# Bench
+@timev access_trajectory_col(SIR_col_buffer) 
+b1_col = @benchmark access_trajectory_col($SIR_col_buffer)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Row buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+function access_trajectory_row(m::Model)
+    res = 0
+    n_sim = 100
+    for k = 1:n_sim
+        σ = _simulate_row_buffer(m)
+        t = _get_values_row(σ, "I")
+        res += t[end-1]
+    end
+    return res
+end
+# Bench
+@timev access_trajectory_row(SIR_row_buffer) 
+b1_row = @benchmark access_trajectory_row($SIR_row_buffer)
+@show minimum(b1_row), mean(b1_row), maximum(b1_row)
+

bench/array_memory_order/multiple_sim_sir.jl

+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+l_var = ["S", "I", "R"]
+bound_time = 200.0
+nbr_sim = 10000
+
+load_model("_bench_perf_test/SIR_col")
+SIR_col.time_bound = bound_time
+set_observed_var!(SIR_col, l_var)
+
+println("Col")
+b1_col = @benchmark for i = 1:$(nbr_sim) _simulate_col($SIR_col) end
+@timev _simulate_col(SIR_col)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+
+println("Col + buffer:")
+b1_col_buffer = @benchmark for i = 1:$(nbr_sim) _simulate_col_buffer($SIR_col_buffer) end
+@timev _simulate_col_buffer(SIR_col_buffer)
+@show minimum(b1_col_buffer), mean(b1_col_buffer), maximum(b1_col_buffer)
+
+println("Col + buffer_10:")
+b1_col_buffer_10 = @benchmark for i = 1:$(nbr_sim) _simulate_col_buffer($SIR_col_buffer; buffer_size = 10) end
+@timev _simulate_col_buffer(SIR_col_buffer; buffer_size = 10)
+@show minimum(b1_col_buffer_10), mean(b1_col_buffer_10), maximum(b1_col_buffer_10)
+
+load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+
+println("Row + buffer:")
+b1_row_buffer = @benchmark for i = 1:$(nbr_sim) _simulate_row_buffer($SIR_row_buffer) end
+@timev _simulate_row_buffer(SIR_row_buffer)
+@show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
+
+println("Row + buffer_10:")
+b1_row_buffer_10 = @benchmark for i = 1:$(nbr_sim) _simulate_row_buffer($SIR_row_buffer; buffer_size = 10) end
+@timev _simulate_row_buffer(SIR_row_buffer; buffer_size = 10)
+@show minimum(b1_row_buffer_10), mean(b1_row_buffer_10), maximum(b1_row_buffer_10)
+

bench/array_memory_order/read_random_state_trajectory.jl

+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+l_var = ["S", "I", "R"]
+bound_time = 200.0
+
+println("Col buffer:")
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+function random_trajectory_value_col(m::ContinuousTimeModel)
+    σ = _simulate_col_buffer(m)
+    n_states = get_states_number(σ)
+    nb_rand = 1000
+    res = 0
+    for i = 1:nb_rand
+        a = _get_state_col(σ, rand(1:n_states))
+        res += a[2]
+    end
+    return res
+end
+# Bench
+@timev random_trajectory_value_col(SIR_col_buffer) 
+b1_col_buffer = @benchmark random_trajectory_value_col($SIR_col_buffer)
+@show minimum(b1_col_buffer), mean(b1_col_buffer), maximum(b1_col_buffer)
+
+println("Row buffer:")
+
+load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+function random_trajectory_value_row(m::ContinuousTimeModel)
+    σ = _simulate_row_buffer(m)
+    n_states = get_states_number(σ)
+    nb_rand = 1000
+    res = 0
+    for i = 1:nb_rand
+        a = _get_state_row(σ, rand(1:n_states))
+        res += a[2]
+    end
+    return res
+end
+# Bench
+@timev random_trajectory_value_row(SIR_row_buffer) 
+b1_row_buffer = @benchmark random_trajectory_value_row($SIR_row_buffer)
+@show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
+

bench/array_memory_order/read_trajectory_sir.jl

+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+bound_time = 200.0
+l_var = ["S", "I", "R"]
+
+println("Col buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+function read_trajectory_col(m::Model)
+    res = 0
+    σ = _simulate_col_buffer(m)
+    n_states = get_states_number(σ)
+    n_read = 100000
+    for k = 1:n_read
+        for i = 1:n_states
+            res += _get_value_col(σ, "I", i)
+        end
+    end
+    return res
+end
+# Bench
+@timev read_trajectory_col(SIR_col_buffer) 
+b1_col = @benchmark read_trajectory_col($SIR_col_buffer)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Row buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+function read_trajectory_row(m::Model)
+    res = 0
+    σ = _simulate_row_buffer(m)
+    n_states = get_states_number(σ)
+    n_read = 100000
+    for k = 1:n_read
+        for i = 1:n_states
+            res += _get_value_row(σ, "I", i)
+        end
+    end
+    return res
+end
+# Bench
+@timev read_trajectory_row(SIR_row_buffer) 
+b1_row = @benchmark read_trajectory_row($SIR_row_buffer)
+@show minimum(b1_row), mean(b1_row), maximum(b1_row)
+

bench/array_memory_order/sim_sir.jl

+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+l_var = ["S", "I", "R"]
+bound_time = 200.0
+
+load_model("_bench_perf_test/SIR_col")
+SIR_col.time_bound = bound_time
+set_observed_var!(SIR_col, l_var)
+
+println("Col")
+b1_col = @benchmark _simulate_col($SIR_col)
+@timev _simulate_col(SIR_col)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+
+println("Col + buffer:")
+b1_col_buffer = @benchmark _simulate_col_buffer($SIR_col_buffer)
+@timev _simulate_col_buffer(SIR_col_buffer)
+@show minimum(b1_col_buffer), mean(b1_col_buffer), maximum(b1_col_buffer)
+
+println("Col + buffer_10:")
+b1_col_buffer_10 = @benchmark _simulate_col_buffer($SIR_col_buffer; buffer_size = 10)
+@timev _simulate_col_buffer(SIR_col_buffer; buffer_size = 10)
+@show minimum(b1_col_buffer_10), mean(b1_col_buffer_10), maximum(b1_col_buffer_10)
+
+load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+
+println("Row + buffer:")
+b1_row_buffer = @benchmark _simulate_row_buffer($SIR_row_buffer)
+@timev _simulate_row_buffer(SIR_row_buffer)
+@show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
+
+println("Row + buffer_10:")
+b1_row_buffer_10 = @benchmark _simulate_row_buffer($SIR_row_buffer; buffer_size = 10)
+@timev _simulate_row_buffer(SIR_row_buffer; buffer_size = 10)
+@show minimum(b1_row_buffer_10), mean(b1_row_buffer_10), maximum(b1_row_buffer_10)
+
+

bench/pygmalion/read_random_state_trajectory.jl

@@ -8,7 +8,8 @@ println("Pygmalion:")
 str_m = "enzymatic_reaction"
 str_d = "abc_er"
 pygmalion.load_model(str_m)
-str_oml = "P,R,time"
+l_var = ["E", "S", "ES", "P"]
+str_oml = "E,S,ES,P,R,time"
 ll_om = split(str_oml, ",")
 x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
 p_true = Parameters(1.0, 1.0, 1.0)
@@ -20,6 +21,10 @@ function random_trajectory_value_pyg(so::SystemObservation)
     n_states = get_number_of_observations(so, "P")
     return to_vec(so, "P", rand(1:n_states))
 end
+function random_trajectory_state_pyg(so::SystemObservation)
+    n_states = get_number_of_observations(so, "P")
+    return to_vec(so, so.oml, rand(1:n_states))
+end
 # Bench
 @timev random_trajectory_value_pyg(so)
 b1_pyg = @benchmark random_trajectory_value_pyg($so)

core/_tests_simulate.jl

@@ -4,14 +4,14 @@
 # Trajectories
 
 _get_values_col(σ::AbstractTrajectory, var::String) = 
-σ.values[(σ.m)._map_obs_var_idx[var],:] 
+@view σ.values[(σ.m)._map_obs_var_idx[var],:] 
 _get_values_row(σ::AbstractTrajectory, var::String) = 
-σ.values[:,(σ.m)._map_obs_var_idx[var]] 
+@view σ.values[:,(σ.m)._map_obs_var_idx[var]] 
 
 _get_state_col(σ::AbstractTrajectory, idx::Int) = 
-σ.values[:,idx]
+@view σ.values[:,idx]
 _get_state_row(σ::AbstractTrajectory, idx::Int) = 
-σ.values[idx,:]
+@view σ.values[idx,:]
 
 _get_value_col(σ::AbstractTrajectory, var::String, idx::Int) = 
 σ.values[(σ.m)._map_obs_var_idx[var],idx] 

core/model.jl

@@ -32,6 +32,14 @@ function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_tra
         _g_idx[i] = map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::String => idx in state space )
         _map_obs_var_idx[g[i]] = i
     end
+  
+    if length(methods(f!)) >= 2
+        @warn "You have possibly redefined a function Model.f! used in a previously instantiated model."
+    end
+    if length(methods(is_absorbing)) >= 2
+        @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)
 end
 

core/observations.jl

@@ -4,9 +4,9 @@ ContinuousObservations = AbstractVector{AbstractTrajectory}
 
 struct Trajectory <: AbstractTrajectory
     m::ContinuousTimeModel
-    values::AbstractMatrix{Float64}
-    times::AbstractVector{Float64}
-    transitions::AbstractVector{Union{String,Nothing}}
+    values::Matrix{Int}
+    times::Vector{Float64}
+    transitions::Vector{Union{String,Nothing}}
 end
 
 function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end

models/SIR.jl

@@ -9,7 +9,7 @@ l_tr = ["R1","R2"]
 p = [0.0012, 0.05]
 x0 = [95, 5, 0]
 t0 = 0.0
-function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
+function SIR_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
            xn::Vector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
@@ -41,9 +41,9 @@ function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
     tnplus1[1] = tn + tau
     tr[1] = "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::Vector{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,f!,is_absorbing_sir; g=g)
+SIR = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_f!,is_absorbing_SIR; g=g)
 export SIR
 

models/_bench_perf_test/SIR_col.jl

@@ -9,7 +9,7 @@ l_tr = ["R1","R2"]
 p = [0.0012, 0.05]
 x0 = [95, 5, 0]
 t0 = 0.0
-function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
+function SIR_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
            xn::Vector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
@@ -41,9 +41,9 @@ function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
     tnplus1[1] = tn + tau
     tr[1] = "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_col(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_sir; g=g)
+SIR_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,is_absorbing_SIR_col; g=g)
 export SIR_col
 

models/_bench_perf_test/SIR_col_buffer.jl

@@ -9,7 +9,7 @@ l_tr = ["R1","R2"]
 p = [0.0012, 0.05]
 x0 = [95, 5, 0]
 t0 = 0.0
-function f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+function SIR_col_buffer_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]
@@ -41,9 +41,9 @@ function f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx:
     l_t[idx] = tn + tau
     l_tr[idx] = "R$(reaction)"
 end
-is_absorbing_sir(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+is_absorbing_SIR_col_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_sir; g=g)
+SIR_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,is_absorbing_SIR_col_buffer; g=g)
 export SIR_col_buffer
 

models/_bench_perf_test/SIR_row_buffer.jl

@@ -9,7 +9,7 @@ l_tr = ["R1","R2"]
 p = [0.0012, 0.05]
 x0 = [95, 5, 0]
 t0 = 0.0
-function f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+function SIR_row_buffer_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]
@@ -41,9 +41,9 @@ function f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx:
     l_t[idx] = tn + tau
     l_tr[idx] = "R$(reaction)"
 end
-is_absorbing_sir(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+is_absorbing_SIR_row_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_sir; g=g)
+SIR_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,is_absorbing_SIR_row_buffer; g=g)
 export SIR_row_buffer