From 7aa17e74315eeec5ddee4e7fdf809291b9279542 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Tue, 17 Nov 2020 00:25:51 +0100
Subject: [PATCH] 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).
---
.../access_trajectory_sir.jl | 49 +++++++++++++++++
bench/array_memory_order/multiple_sim_sir.jl | 48 +++++++++++++++++
.../read_random_state_trajectory.jl | 51 ++++++++++++++++++
.../array_memory_order/read_trajectory_sir.jl | 53 +++++++++++++++++++
bench/array_memory_order/sim_sir.jl | 48 +++++++++++++++++
.../pygmalion/read_random_state_trajectory.jl | 7 ++-
core/_tests_simulate.jl | 8 +--
core/model.jl | 8 +++
core/observations.jl | 6 +--
models/SIR.jl | 6 +--
models/_bench_perf_test/SIR_col.jl | 6 +--
models/_bench_perf_test/SIR_col_buffer.jl | 6 +--
models/_bench_perf_test/SIR_row_buffer.jl | 6 +--
13 files changed, 282 insertions(+), 20 deletions(-)
create mode 100644 bench/array_memory_order/access_trajectory_sir.jl
create mode 100644 bench/array_memory_order/multiple_sim_sir.jl
create mode 100644 bench/array_memory_order/read_random_state_trajectory.jl
create mode 100644 bench/array_memory_order/read_trajectory_sir.jl
create mode 100644 bench/array_memory_order/sim_sir.jl
diff --git a/bench/array_memory_order/access_trajectory_sir.jl b/bench/array_memory_order/access_trajectory_sir.jl
new file mode 100644
index 0000000..6189405
--- /dev/null
+++ b/bench/array_memory_order/access_trajectory_sir.jl
@@ -0,0 +1,49 @@
+
+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)
+
diff --git a/bench/array_memory_order/multiple_sim_sir.jl b/bench/array_memory_order/multiple_sim_sir.jl
new file mode 100644
index 0000000..2a9efbd
--- /dev/null
+++ b/bench/array_memory_order/multiple_sim_sir.jl
@@ -0,0 +1,48 @@
+
+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)
+
diff --git a/bench/array_memory_order/read_random_state_trajectory.jl b/bench/array_memory_order/read_random_state_trajectory.jl
new file mode 100644
index 0000000..c8efbcf
--- /dev/null
+++ b/bench/array_memory_order/read_random_state_trajectory.jl
@@ -0,0 +1,51 @@
+
+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)
+
diff --git a/bench/array_memory_order/read_trajectory_sir.jl b/bench/array_memory_order/read_trajectory_sir.jl
new file mode 100644
index 0000000..2a32ccf
--- /dev/null
+++ b/bench/array_memory_order/read_trajectory_sir.jl
@@ -0,0 +1,53 @@
+
+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)
+
diff --git a/bench/array_memory_order/sim_sir.jl b/bench/array_memory_order/sim_sir.jl
new file mode 100644
index 0000000..84a5d14
--- /dev/null
+++ b/bench/array_memory_order/sim_sir.jl
@@ -0,0 +1,48 @@
+
+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)
+
+
diff --git a/bench/pygmalion/read_random_state_trajectory.jl b/bench/pygmalion/read_random_state_trajectory.jl
index 05c444c..9877751 100644
--- a/bench/pygmalion/read_random_state_trajectory.jl
+++ b/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)
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index 06e01d1..1562f41 100644
--- a/core/_tests_simulate.jl
+++ b/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]
diff --git a/core/model.jl b/core/model.jl
index d37ad6d..a1b584c 100644
--- a/core/model.jl
+++ b/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
diff --git a/core/observations.jl b/core/observations.jl
index 5967ea5..e82b222 100644
--- a/core/observations.jl
+++ b/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
diff --git a/models/SIR.jl b/models/SIR.jl
index 7416631..d448cc2 100644
--- a/models/SIR.jl
+++ b/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
diff --git a/models/_bench_perf_test/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl
index 00cf2e4..919e8fc 100644
--- a/models/_bench_perf_test/SIR_col.jl
+++ b/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
diff --git a/models/_bench_perf_test/SIR_col_buffer.jl b/models/_bench_perf_test/SIR_col_buffer.jl
index d37828b..9f79b3e 100644
--- a/models/_bench_perf_test/SIR_col_buffer.jl
+++ b/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
diff --git a/models/_bench_perf_test/SIR_row_buffer.jl b/models/_bench_perf_test/SIR_row_buffer.jl
index 366b68b..fcf6a86 100644
--- a/models/_bench_perf_test/SIR_row_buffer.jl
+++ b/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
--
GitLab