From 8020a182e26e50a6d8d046f2eed898514154576a Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Wed, 18 Nov 2020 15:14:12 +0100
Subject: [PATCH] Fix of several type unstabilities due to @view macro. Add of
bench for @view macro.
---
.../multiple_dist_lp_view.jl | 43 +++++++++++++++
bench/array_memory_order/multiple_sim_view.jl | 28 ++++++++++
core/_tests_simulate.jl | 53 +++++++++++++++++--
core/model.jl | 2 +-
core/trajectory.jl | 6 +--
models/ER.jl | 4 +-
models/SIR.jl | 4 +-
models/_bench_perf_test/SIR_col.jl | 4 +-
tests/dist_lp/dist_case_2.jl | 5 +-
tests/dist_lp/dist_case_3.jl | 4 +-
tests/dist_lp/dist_l1_case_1.jl | 4 +-
11 files changed, 137 insertions(+), 20 deletions(-)
create mode 100644 bench/array_memory_order/multiple_dist_lp_view.jl
create mode 100644 bench/array_memory_order/multiple_sim_view.jl
diff --git a/bench/array_memory_order/multiple_dist_lp_view.jl b/bench/array_memory_order/multiple_dist_lp_view.jl
new file mode 100644
index 0000000..49972a4
--- /dev/null
+++ b/bench/array_memory_order/multiple_dist_lp_view.jl
@@ -0,0 +1,43 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+str_model = ARGS[1]
+load_model(str_model)
+if str_model == "SIR"
+ model = SIR
+ SIR.time_bound = 150
+elseif str_model == "ER"
+ model = ER
+ ER.time_bound = 10.0
+else
+ error("Unrecognized model")
+end
+
+
+function compute_dist(m::Model)
+ nbr_sim = 10000
+ for i = 1:nbr_sim
+ σ1, σ2 = simulate(m), simulate(m)
+ dist_lp(σ1, σ2)
+ end
+end
+
+function compute_dist_without_view(m::Model)
+ nbr_sim = 10000
+ for i = 1:nbr_sim
+ σ1, σ2 = simulate(m), simulate(m)
+ dist_lp(σ1, σ2)
+ end
+end
+
+@timev simulate(model)
+b1 = @benchmark compute_dist($model)
+@show minimum(b1), mean(b1), maximum(b1)
+
+@timev _simulate_without_view(model)
+b1_without_view = @benchmark compute_dist_without_view($model)
+@show minimum(b1_without_view), mean(b1_without_view), maximum(b1_without_view)
+
diff --git a/bench/array_memory_order/multiple_sim_view.jl b/bench/array_memory_order/multiple_sim_view.jl
new file mode 100644
index 0000000..6b1af30
--- /dev/null
+++ b/bench/array_memory_order/multiple_sim_view.jl
@@ -0,0 +1,28 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+str_model = ARGS[1]
+load_model(str_model)
+if str_model == "SIR"
+ model = SIR
+ SIR.time_bound = 150
+elseif str_model == "ER"
+ model = ER
+ ER.time_bound = 10.0
+else
+ error("Unrecognized model")
+end
+
+nbr_sim = 10000
+
+@timev simulate(model)
+b1 = @benchmark for i = 1:$(nbr_sim) simulate($model) end
+@show minimum(b1), mean(b1), maximum(b1)
+
+@timev _simulate_without_view(model)
+b1_without_view = @benchmark for i = 1:$(nbr_sim) _simulate_without_view($model) end
+@show minimum(b1_without_view), mean(b1_without_view), maximum(b1_without_view)
+
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index 1562f41..ca97957 100644
--- a/core/_tests_simulate.jl
+++ b/core/_tests_simulate.jl
@@ -27,7 +27,7 @@ function _simulate_col(m::ContinuousTimeModel)
transitions = Vector{Union{String,Nothing}}(undef,0)
# values at time n
n = 0
- xn = m.x0
+ xn = @view m.x0[:]
tn = m.t0
tr = [""]
# at time n+1
@@ -61,7 +61,7 @@ function _simulate_row(m::ContinuousTimeModel)
transitions = Vector{Union{String,Nothing}}(undef,0)
# values at time n
n = 0
- xn = m.x0
+ xn = @view m.x0[:]
tn = m.t0
tr = [""]
# at time n+1
@@ -96,7 +96,7 @@ function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
transitions = Vector{Union{String,Nothing}}(undef,0)
# values at time n
n = 0
- xn = m.x0
+ xn = @view m.x0[:]
tn = m.t0
# at time n+1
mat_x = zeros(Int, m.d, buffer_size)
@@ -136,7 +136,7 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
transitions = Vector{Union{String,Nothing}}(undef,0)
# values at time n
n = 0
- xn = m.x0
+ xn = @view m.x0[:]
tn = m.t0
# at time n+1
mat_x = zeros(Int, buffer_size, m.d)
@@ -169,3 +169,48 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
return Trajectory(m, values, times, transitions)
end
+function _simulate_without_view(m::ContinuousTimeModel)
+ # trajectory fields
+ full_values = Matrix{Int}(undef, 1, m.d)
+ full_values[1,:] = m.x0
+ times = Float64[m.t0]
+ transitions = Union{String,Nothing}[nothing]
+ # values at time n
+ n = 0
+ xn = m.x0
+ tn = m.t0
+ # at time n+1
+ 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)
+ 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 = mat_x[i,:]
+ tn = l_t[i]
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ full_values = vcat(full_values, mat_x[1:i,:])
+ append!(times, l_t[1:i])
+ append!(transitions, l_tr[1:i])
+ n += i
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ if is_bounded(m)
+ if times[end] > m.time_bound
+ full_values[end,:] = full_values[end-1,:]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ else
+ full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
+ push!(times, m.time_bound)
+ push!(transitions, nothing)
+ end
+ end
+ values = full_values[:,m._g_idx]
+ return Trajectory(m, values, times, transitions)
+end
+
diff --git a/core/model.jl b/core/model.jl
index f351a5f..541a4bb 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -53,7 +53,7 @@ function simulate(m::ContinuousTimeModel)
transitions = Union{String,Nothing}[nothing]
# values at time n
n = 0
- xn = m.x0
+ xn = @view m.x0[:] # View for type stability
tn = m.t0
# at time n+1
mat_x = zeros(Int, m.buffer_size, m.d)
diff --git a/core/trajectory.jl b/core/trajectory.jl
index 5fdb3aa..5ef4940 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -78,7 +78,7 @@ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String;
end
# Distance function. Vectorized version
-function dist_lp(x_obs::AbstractVector{Int}, t_x::Vector{Float64}, y_obs::AbstractVector{Int}, t_y::Vector{Float64};
+function dist_lp(x_obs::SubArray{Int,1}, t_x::SubArray{Float64,1}, y_obs::SubArray{Int,1}, t_y::SubArray{Float64,1};
verbose::Bool = false, p::Int = 1)
current_y_obs = y_obs[1]
current_t_y = t_y[2]
@@ -188,9 +188,9 @@ function δ(σ1::AbstractTrajectory,t::Float64) end
# Get var values ["I"]
function getindex(σ::AbstractTrajectory, var::String)
if var == "times"
- return σ.times
+ return @view σ.times[:]
elseif var == "transitions"
- return σ.transitions
+ return @view σ.transitions[:]
else
return get_var_values(σ, var)
end
diff --git a/models/ER.jl b/models/ER.jl
index 7fa9477..3c694b9 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -10,7 +10,7 @@ p = [1.0, 1.0, 1.0]
x0 = [100, 100, 0, 0]
t0 = 0.0
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})
+ xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
a3 = p[3] * xn[3]
@@ -41,7 +41,7 @@ function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, i
l_t[idx] = tn + tau
l_tr[idx] = "R$(reaction)"
end
-is_absorbing_ER(p::Vector{Float64},xn::AbstractVector{Int}) =
+is_absorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g = ["P"]
diff --git a/models/SIR.jl b/models/SIR.jl
index 6489084..215a4a8 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -10,7 +10,7 @@ p = [0.0012, 0.05]
x0 = [95, 5, 0]
t0 = 0.0
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})
+ xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
l_a = SVector(a1, a2)
@@ -41,7 +41,7 @@ function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String},
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(p::Vector{Float64}, xn::SubArray{Int,1}) = (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)
diff --git a/models/_bench_perf_test/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl
index 919e8fc..09d2b2a 100644
--- a/models/_bench_perf_test/SIR_col.jl
+++ b/models/_bench_perf_test/SIR_col.jl
@@ -10,7 +10,7 @@ p = [0.0012, 0.05]
x0 = [95, 5, 0]
t0 = 0.0
function SIR_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
- xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+ 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)
@@ -41,7 +41,7 @@ function SIR_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{S
tnplus1[1] = tn + tau
tr[1] = "R$(reaction)"
end
-is_absorbing_SIR_col(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::AbstractVector{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,SIR_col_f!,is_absorbing_SIR_col; g=g)
diff --git a/tests/dist_lp/dist_case_2.jl b/tests/dist_lp/dist_case_2.jl
index 786c0da..ce29df4 100644
--- a/tests/dist_lp/dist_case_2.jl
+++ b/tests/dist_lp/dist_case_2.jl
@@ -4,6 +4,7 @@ import QuadGK: quadgk
load_model("SIR")
test_all = true
+p=2
for p = 1:2
res = (4-1)^p * 0.5 + (4-2)^p * 0.1 + 0 + (5-3)^p * 0.1 + (5-1)^p * (1.4 - 0.8) + (3-1)^p * (3.0-1.4)
res = res^(1/p)
@@ -22,9 +23,9 @@ for p = 1:2
l_tr = Vector{Nothing}(nothing, length(x_obs))
σ1 = Trajectory(SIR, values, t_x, l_tr)
- test_1 = isapprox(dist_lp(x_obs, t_x, y_obs, t_y; p=p), res; atol = 1E-10)
+ test_1 = isapprox(dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=p), res; atol = 1E-10)
test_1 = test_1 && isapprox(dist_lp(σ1,σ2; p=p), res; atol = 1E-10)
- test_1_bis = isapprox(dist_lp(y_obs, t_y, x_obs, t_x; p=p), res; atol = 1E-10)
+ test_1_bis = isapprox(dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=p), res; atol = 1E-10)
test_1_bis = test_1_bis && isapprox(dist_lp(σ2,σ1; p=p), res; atol = 1E-10)
f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
diff --git a/tests/dist_lp/dist_case_3.jl b/tests/dist_lp/dist_case_3.jl
index 149a517..c7d9558 100644
--- a/tests/dist_lp/dist_case_3.jl
+++ b/tests/dist_lp/dist_case_3.jl
@@ -26,9 +26,9 @@ for p = 1:2
int_riemann = MarkovProcesses._riemann_sum(diff_f, 0.0, 20.0, 1E-5)
int_riemann = int_riemann^(1/p)
- res1 = dist_lp(x_obs, t_x, y_obs, t_y; p=p)
+ res1 = dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=p)
res2 = dist_lp(σ1,σ2; p=p)
- res1_bis = dist_lp(y_obs, t_y, x_obs, t_x; p=p)
+ res1_bis = dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=p)
res2_bis = dist_lp(σ2,σ1; p=p)
test_1 = isapprox(res1, int_riemann; atol = 1E-3)
test_1 = test_1 && isapprox(res2, int_riemann; atol = 1E-3)
diff --git a/tests/dist_lp/dist_l1_case_1.jl b/tests/dist_lp/dist_l1_case_1.jl
index 08e9d7e..d02b7f3 100644
--- a/tests/dist_lp/dist_l1_case_1.jl
+++ b/tests/dist_lp/dist_l1_case_1.jl
@@ -21,7 +21,7 @@ values[:,1] = y_obs
l_tr = Vector{Nothing}(nothing, length(y_obs))
σ2 = Trajectory(SIR, values, t_y, l_tr)
-test_1 = dist_lp(x_obs, t_x, y_obs, t_y; p=1) == 6.4
+test_1 = dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=1) == 6.4
test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4
f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
@@ -33,7 +33,7 @@ test_2 = isapprox(6.4, int; atol=err)
# Case 1 bis : inverse of case 1
-test_1_bis = dist_lp(y_obs, t_y, x_obs, t_x; p=1) == 6.4
+test_1_bis = dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=1) == 6.4
test_1_bis = test_1_bis && dist_lp(σ2, σ1; p=1) == 6.4
return test_1 && test_1_bis && test_2
--
GitLab