From 5b7fc78fc7daa9739e326196bd7622bcd35a0283 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Mon, 23 Nov 2020 22:27:35 +0100
Subject: [PATCH] Change of Trajectory.values type from Matrix to Vector of
Vector in order to replace vcat by push!.
The benchmarks are at least equal, and better when the buffer size is
small. However, no significative difference of performance with adequate
buffer size.
My opinion the gain can be seen when simulations are going to be very
long with more complicated models than we have implemented for now.
All tests passed.
---
bench/pygmalion/multiple_long_sim_er.jl | 4 +-
core/MarkovProcesses.jl | 2 +-
core/_tests_simulate.jl | 82 +++++++++++++++++++---
core/common.jl | 4 +-
core/lha.jl | 9 ++-
core/trajectory.jl | 24 ++++---
tests/dist_lp/dist_case_2.jl | 12 ++--
tests/dist_lp/dist_case_3.jl | 12 ++--
tests/dist_lp/dist_l1_case_1.jl | 12 ++--
tests/unit/check_trajectory_consistency.jl | 20 +++++-
10 files changed, 135 insertions(+), 46 deletions(-)
diff --git a/bench/pygmalion/multiple_long_sim_er.jl b/bench/pygmalion/multiple_long_sim_er.jl
index 6729aba..23244dc 100644
--- a/bench/pygmalion/multiple_long_sim_er.jl
+++ b/bench/pygmalion/multiple_long_sim_er.jl
@@ -38,8 +38,8 @@ b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
@show minimum(b2), mean(b2), maximum(b2)
println("Buffer size 100")
ER.buffer_size = 100
-b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
-b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
+b1_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
+b2_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
@show minimum(b1), mean(b1), maximum(b1)
@show minimum(b2), mean(b2), maximum(b2)
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 5fb0ec4..85eb0f6 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -22,7 +22,7 @@ export load_automaton, get_index, get_value, length_var, isaccepted
# Model related methods
export simulate, set_param!, get_param, set_observed_var!
-export set_time_bound!
+export set_time_bound!, getproperty
export isbounded, isaccepted, check_consistency
export load_model, get_module_path
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index 632b3ff..9c40acb 100644
--- a/core/_tests_simulate.jl
+++ b/core/_tests_simulate.jl
@@ -1,21 +1,28 @@
+struct OldTrajectory <: AbstractTrajectory
+ m::ContinuousTimeModel
+ values::Matrix{Int}
+ times::Vector{Float64}
+ transitions::Vector{Union{String,Nothing}}
+end
+
# File for benchmarking simulation and memory access of the package.
# Trajectories
-_get_values_col(σ::AbstractTrajectory, var::String) =
+_get_values_col(σ::OldTrajectory, var::String) =
@view σ.values[(σ.m)._map_obs_var_idx[var],:]
-_get_values_row(σ::AbstractTrajectory, var::String) =
+_get_values_row(σ::OldTrajectory, var::String) =
@view σ.values[:,(σ.m)._map_obs_var_idx[var]]
-_get_state_col(σ::AbstractTrajectory, idx::Int) =
+_get_state_col(σ::OldTrajectory, idx::Int) =
@view σ.values[:,idx]
-_get_state_row(σ::AbstractTrajectory, idx::Int) =
+_get_state_row(σ::OldTrajectory, idx::Int) =
@view σ.values[idx,:]
-_get_value_col(σ::AbstractTrajectory, var::String, idx::Int) =
+_get_value_col(σ::OldTrajectory, var::String, idx::Int) =
σ.values[(σ.m)._map_obs_var_idx[var],idx]
-_get_value_row(σ::AbstractTrajectory, var::String, idx::Int) =
+_get_value_row(σ::OldTrajectory, var::String, idx::Int) =
σ.values[idx,(σ.m)._map_obs_var_idx[var]]
# Model
@@ -51,7 +58,7 @@ function _simulate_col(m::ContinuousTimeModel)
transitions[end] = nothing
end
end
- return Trajectory(m, values, times, transitions)
+ return OldTrajectory(m, values, times, transitions)
end
function _simulate_row(m::ContinuousTimeModel)
@@ -85,7 +92,7 @@ function _simulate_row(m::ContinuousTimeModel)
transitions[end] = nothing
end
end
- return Trajectory(m, values, times, transitions)
+ return OldTrajectory(m, values, times, transitions)
end
@@ -126,7 +133,7 @@ function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
transitions[end] = nothing
end
end
- return Trajectory(m, values, times, transitions)
+ return OldTrajectory(m, values, times, transitions)
end
function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
@@ -166,7 +173,7 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
transitions[end] = nothing
end
end
- return Trajectory(m, values, times, transitions)
+ return OldTrajectory(m, values, times, transitions)
end
function _simulate_without_view(m::ContinuousTimeModel)
@@ -211,6 +218,59 @@ function _simulate_without_view(m::ContinuousTimeModel)
end
end
values = full_values[:,m._g_idx]
- return Trajectory(m, values, times, transitions)
+ return OldTrajectory(m, values, times, transitions)
end
+# With trajectory values in Matrix type
+function _simulate_27d56(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 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
+ 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)
+ isabsorbing::Bool = m.isabsorbing(m.p,xn)
+ # use sizehint! ?
+ while !isabsorbing &&Â tn <= m.time_bound
+ i = 0
+ while i < m.buffer_size && !isabsorbing &&Â tn <= m.time_bound
+ i += 1
+ m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
+ tn = l_t[i]
+ if tn > m.time_bound
+ i -= 1 #Â 0 is an ok value, 1:0 is allowed
+ break
+ end
+ xn = view(mat_x, i, :)
+ isabsorbing = m.isabsorbing(m.p,xn)
+ 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
+ end
+ if isbounded(m)
+ #=
+ if times[end] > m.time_bound
+ full_values[end,:] = full_values[end-1,:]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ else
+ end
+ =#
+ full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
+ push!(times, m.time_bound)
+ push!(transitions, nothing)
+ end
+ values = view(full_values, :, m._g_idx)
+ return OldTrajectory(m, values, times, transitions)
+end
+
+
diff --git a/core/common.jl b/core/common.jl
index bbe0914..8f7b015 100644
--- a/core/common.jl
+++ b/core/common.jl
@@ -26,7 +26,7 @@ end
struct Trajectory <: AbstractTrajectory
m::ContinuousTimeModel
- values::Matrix{Int}
+ values::Vector{Vector{Int}}
times::Vector{Float64}
transitions::Vector{Transition}
end
@@ -65,7 +65,7 @@ end
struct SynchronizedTrajectory <: AbstractTrajectory
S::StateLHA
m::SynchronizedModel
- values::Matrix{Int}
+ values::Vector{Vector{Int}}
times::Vector{Float64}
transitions::Vector{Transition}
end
diff --git a/core/lha.jl b/core/lha.jl
index c47b7b6..2445b23 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -24,7 +24,7 @@ end
isaccepted(S::StateLHA) = (S.loc in (S.A).l_loc_final)
# Methods for synchronize / read the trajectory
-function init_state(A::LHA, x0::SubArray{Int,1}, t0::Float64)
+function init_state(A::LHA, x0::Vector{Int}, t0::Float64)
S0 = StateLHA(A, "", zeros(length_var(A)), t0)
for loc in A.l_loc_init
if A.Λ[loc](A,S0)
@@ -119,12 +119,17 @@ function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
A_new = LHA(A, (σ.m)._map_obs_var_idx)
l_t = times(σ)
l_tr = transitions(σ)
+ mat_x = zeros(Int, length_states(σ), σ.m.d)
+ for (i,var) in enumerate(σ.m.g)
+ mat_x[:,i] = σ[var]
+ end
Sn = init_state(A_new, σ[1], l_t[1])
Snplus1 = copy(Sn)
if verbose println("Init: ") end
if verbose @show Sn end
for n in 1:length_states(σ)
- next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn; verbose = verbose)
+ xn = view(mat_x, n, :)
+ next_state!(Snplus1, A_new, xn, l_t[n], l_tr[n], Sn; verbose = verbose)
if Snplus1.loc in A_new.l_loc_final
break
end
diff --git a/core/trajectory.jl b/core/trajectory.jl
index cb8901b..a02735c 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -64,7 +64,7 @@ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String;
end
# Distance function. Vectorized version
-function dist_lp(x_obs::SubArray{Int,1}, t_x::Vector{Float64}, y_obs::SubArray{Int,1}, t_y::Vector{Float64};
+function dist_lp(x_obs::Vector{Int}, t_x::Vector{Float64}, y_obs::Vector{Int}, t_y::Vector{Float64};
verbose::Bool = false, p::Int = 1)
current_y_obs = y_obs[1]
current_t_y = t_y[2]
@@ -150,23 +150,31 @@ function _riemann_sum(f::Function, t_begin::Real, t_end::Real, step::Float64)
end
function check_consistency(σ::AbstractTrajectory)
- @assert length(σ.times) == length(σ.transitions) == size(σ.values)[1]
- @assert length_obs_var(σ) == length(σ.m.g) == size(σ.values)[2]
+ test_length_var = true
+ for i = 1:(σ.m).dobs
+ test_length_i = (length(σ.values[1]) == length(σ.values[i]))
+ test_length_var = test_length_var && test_length_i
+ end
+ @assert begin (length(σ.times) == length(σ.transitions)) &&
+ (length(σ.times) == length(σ.values[1])) &&
+ test_length_var
+ end
+ @assert length_obs_var(σ) == σ.m.dobs
return true
end
-issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, σ[end])
+issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, view(reshape(σ[end], 1, m.d), 1, :))
# Properties of the trajectory
length_states(σ::AbstractTrajectory) = length(σ.times)
-length_obs_var(σ::AbstractTrajectory) = size(σ.values)[2]
+length_obs_var(σ::AbstractTrajectory) = length(σ.values)
get_obs_var(σ::AbstractTrajectory) = (σ.m).g
isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing
isaccepted(σ::SynchronizedTrajectory) = isaccepted(σ.S)
# 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]]
+get_var_values(σ::AbstractTrajectory, var::String) = σ.values[(σ.m)._map_obs_var_idx[var]]
+get_state(σ::AbstractTrajectory, idx::Int) = [σ.values[i][idx] for i = 1:length(σ.values)] # /!\ Creates an array
+get_value(σ::AbstractTrajectory, var::String, idx::Int) = get_var_values(σ, var)[idx]
# Operation @
function get_state_from_time(σ::AbstractTrajectory, t::Float64)
@assert t >= 0.0
diff --git a/tests/dist_lp/dist_case_2.jl b/tests/dist_lp/dist_case_2.jl
index be7c86d..c66ce3b 100644
--- a/tests/dist_lp/dist_case_2.jl
+++ b/tests/dist_lp/dist_case_2.jl
@@ -11,21 +11,21 @@ for p = 1:2
y_obs = [1, 2, 5, 3, 3]
t_y = [0.0, 0.5, 0.7, 1.4, 3.0]
- values = zeros(length(y_obs), 1)
- values[:,1] = y_obs
+ values = [zeros(length(y_obs))]
+ values[1] = y_obs
l_tr = Vector{Nothing}(nothing, length(y_obs))
σ2 = Trajectory(SIR, values, t_y, l_tr)
x_obs = [4, 2, 3, 1, 1]
t_x = [0.0, 0.6, 0.7, 0.8, 3.0]
- values = zeros(length(x_obs), 1)
- values[:,1] = x_obs
+ values = [zeros(length(x_obs))]
+ values[1] = x_obs
l_tr = Vector{Nothing}(nothing, length(x_obs))
σ1 = Trajectory(SIR, values, t_x, l_tr)
- test_1 = isapprox(dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p), res; atol = 1E-10)
+ test_1 = isapprox(dist_lp(x_obs, t_x, y_obs, 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(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; 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 = 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 b271771..8254be0 100644
--- a/tests/dist_lp/dist_case_3.jl
+++ b/tests/dist_lp/dist_case_3.jl
@@ -7,15 +7,15 @@ test_all = true
for p = 1:2
x_obs =[5, 6, 5, 4, 3, 2, 1, 1]
t_x = [0.0, 3.10807, 4.29827, 4.40704, 5.67024, 7.1299, 11.2763, 20.0]
- values = zeros(length(x_obs), 1)
- values[:,1] = x_obs
+ values = [zeros(length(x_obs))]
+ values[1] = x_obs
l_tr = Vector{Nothing}(nothing, length(x_obs))
σ1 = Trajectory(SIR, values, t_x, l_tr)
y_obs =[5, 4, 5, 4, 3, 4, 3, 2, 3, 2, 1, 2, 3, 4, 3, 4, 4]
t_y = [0.0, 0.334082, 1.21012, 1.40991, 1.58866, 2.45879, 2.94545, 4.66746, 5.44723, 5.88066, 7.25626, 11.4036, 13.8373, 17.1363, 17.8193, 18.7613, 20.0]
- values = zeros(length(y_obs), 1)
- values[:,1] = y_obs
+ values = [zeros(length(y_obs))]
+ values[1] = y_obs
l_tr = Vector{Nothing}(nothing, length(y_obs))
σ2 = Trajectory(SIR, values, t_y, l_tr)
@@ -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(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p)
+ res1 = dist_lp(x_obs, t_x, y_obs, t_y; p=p)
res2 = dist_lp(σ1,σ2; p=p)
- res1_bis = dist_lp(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=p)
+ res1_bis = dist_lp(y_obs, t_y, x_obs, 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 a9152e0..262c8c7 100644
--- a/tests/dist_lp/dist_l1_case_1.jl
+++ b/tests/dist_lp/dist_l1_case_1.jl
@@ -8,20 +8,20 @@ load_model("SIR")
x_obs = [2, 4, 3, 3]
t_x = [0.0, 0.5, 1.2, 2.2]
t_y_bis = [0.0, 0.5, 1.2, 2.2]
-values = zeros(length(x_obs), 1)
-values[:,1] = x_obs
+values = [zeros(length(x_obs))]
+values[1] = x_obs
l_tr = Vector{Nothing}(nothing, length(x_obs))
σ1 = Trajectory(SIR, values, t_x, l_tr)
y_obs = [6, 6]
t_y = [0.0, 2.2]
t_x_bis = [0.0, 2.2]
-values = zeros(length(y_obs), 1)
-values[:,1] = y_obs
+values = [zeros(length(y_obs))]
+values[1] = y_obs
l_tr = Vector{Nothing}(nothing, length(y_obs))
σ2 = Trajectory(SIR, values, t_y, l_tr)
-test_1 = dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=1) == 6.4
+test_1 = dist_lp(x_obs, t_x, y_obs, 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(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=1) == 6.4
+test_1_bis = dist_lp(y_obs, t_y, x_obs, 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
diff --git a/tests/unit/check_trajectory_consistency.jl b/tests/unit/check_trajectory_consistency.jl
index 4ae5da9..96ff4b4 100644
--- a/tests/unit/check_trajectory_consistency.jl
+++ b/tests/unit/check_trajectory_consistency.jl
@@ -1,13 +1,29 @@
using MarkovProcesses
+
load_model("SIR")
-SIR.time_bound = 100.0
load_model("ER")
-ER.time_bound = 10.0
test_all = true
nb_sim = 1000
+for i = 1:nb_sim
+ σ = simulate(SIR)
+ σ2 = simulate(ER)
+ try
+ global test_all = test_all && MarkovProcesses.check_consistency(σ) && MarkovProcesses.check_consistency(σ2)
+ catch err
+ @show err
+ @show length(σ.values), length(σ.times), length(σ.transitions)
+ println()
+ @show σ.values
+ @show σ.times
+ @show σ.transitions
+ throw(err)
+ end
+end
+SIR.time_bound = 100.0
+ER.time_bound = 10.0
for i = 1:nb_sim
σ = simulate(SIR)
σ2 = simulate(ER)
--
GitLab