From 35ac88825629e5681aa736af8bcb99e99d8d2559 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Tue, 24 Nov 2020 16:17:17 +0100
Subject: [PATCH] Change of the simulate function with resize!. Gain of
performance is small.
Fix of a misbehaviour of isbounded() for trajectories.
---
core/model.jl | 118 +++++++++++++++------
core/trajectory.jl | 2 +-
models/ER.jl | 2 +-
models/SIR.jl | 2 +-
models/_bench_perf_test/ER_col_buffer.jl | 2 +-
models/_bench_perf_test/ER_row_buffer.jl | 2 +-
models/_bench_perf_test/SIR_col_buffer.jl | 2 +-
models/_bench_perf_test/SIR_row_buffer.jl | 2 +-
tests/unit/check_trajectory_consistency.jl | 81 +++++++++++---
9 files changed, 155 insertions(+), 58 deletions(-)
diff --git a/core/model.jl b/core/model.jl
index e449830..a07e731 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -1,64 +1,112 @@
load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl")
+function _resize_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
+ transitions::Vector{Transition}, size::Int)
+ for i = eachindex(values) resize!(values[i], size) end
+ resize!(times, size)
+ resize!(transitions, size)
+end
+
+
+function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
+ transitions::Vector{Transition}, time_bound::Float64)
+
+ for i = eachindex(values) push!(values[i], values[i][end]) end
+ push!(times, time_bound)
+ push!(transitions, nothing)
+end
+
# Simulation
function simulate(m::ContinuousTimeModel)
- # trajectory fields
+ # First alloc
full_values = Vector{Vector{Int}}(undef, m.d)
- for i = 1:m.d full_values[i] = Int[m.x0[i]] end
- for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
- times = Float64[m.t0]
- transitions = Transition[nothing]
- # values at time n
- n = 0
+ for i = 1:m.d full_values[i] = zeros(Int, m.estim_min_states) end
+ times = zeros(Float64, m.estim_min_states)
+ transitions = Vector{Transition}(undef, m.estim_min_states)
+ # Initial values
+ for i = 1:m.d full_values[i][1] = m.x0[i] end
+ times[1] = m.t0
+ transitions[1] = nothing
+ # Values at time n
+ n = 1
xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
tn = m.t0
# at time n+1
+ isabsorbing::Bool = m.isabsorbing(m.p,xn)
+ if isabsorbing
+ _resize_trajectory!(full_values, times, transitions, 1)
+ values = full_values[m._g_idx]
+ if isbounded(m)
+ _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
+ end
+ return Trajectory(m, values, times, transitions)
+ end
+ # First we fill the allocated array
+ mat_x = zeros(Int, 1, m.d)
+ l_t = Float64[0.0]
+ l_tr = Transition[nothing]
+ for i = 2:m.estim_min_states
+ m.f!(mat_x, l_t, l_tr, 1, xn, tn, m.p)
+ tn = l_t[1]
+ if tn > m.time_bound
+ break
+ end
+ n += 1
+ xn = view(mat_x, 1, :)
+ # Updating value
+ for k = 1:m.d full_values[k][n] = xn[k] end
+ times[n] = tn
+ transitions[n] = l_tr[1]
+ isabsorbing = m.isabsorbing(m.p,xn)
+ if isabsorbing
+ break
+ end
+ end
+ # If simulation ended
+ if n < m.estim_min_states
+ _resize_trajectory!(full_values, times, transitions, n)
+ values = full_values[m._g_idx]
+ if isbounded(m)
+ _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
+ end
+ return Trajectory(m, values, times, transitions)
+ end
+ # Otherwise, buffering system
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)
- end_idx = -1
- # use sizehint! ?
+ l_tr = Vector{Transition}(undef, m.buffer_size)
while !isabsorbing &&Â tn <= m.time_bound
- for i = 1:m.buffer_size
+ i = 0
+ while i < m.buffer_size
+ 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
- end_idx = i
+ i -= 1
break
end
xn = view(mat_x, i, :)
isabsorbing = m.isabsorbing(m.p,xn)
if isabsorbing
- end_idx = i
break
end
end
- if end_idx != -1
- break
- end
- for k = 1:m.d
- append!(full_values[k], view(mat_x, :, k))
- end
- append!(times, l_t)
- append!(transitions, l_tr)
- n += m.buffer_size
- end
- for k = 1:m.d
- append!(full_values[k], view(mat_x, 1:end_idx, k))
+ _resize_trajectory!(full_values, times, transitions, n+i)
+ # Update values
+ rng_traj = (n+1):(n+i)
+ rng_buffer = 1:i
+ for k = 1:m.d full_values[k][rng_traj] = view(mat_x, rng_buffer, k) end
+ times[rng_traj] = l_t[rng_buffer]
+ transitions[rng_traj] = l_tr[rng_buffer]
+ n += i
end
- append!(times, view(l_t, 1:end_idx))
- append!(transitions, view(l_tr, 1:end_idx))
+ values = full_values[m._g_idx]
if isbounded(m)
- for k = 1:m.d
- push!(full_values[k], full_values[k][end])
- end
- push!(times, m.time_bound)
- push!(transitions, nothing)
+ # Add last value: the convention is the last transition is nothing,
+ # the trajectory is bounded
+ _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end
- values = full_values[m._g_idx]
return Trajectory(m, values, times, transitions)
end
diff --git a/core/trajectory.jl b/core/trajectory.jl
index a02735c..9ceb4b0 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -168,7 +168,7 @@ issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, view(reshap
length_states(σ::AbstractTrajectory) = length(σ.times)
length_obs_var(σ::AbstractTrajectory) = length(σ.values)
get_obs_var(σ::AbstractTrajectory) = (σ.m).g
-isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing
+isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing && length_states(σ) >= 2
isaccepted(σ::SynchronizedTrajectory) = isaccepted(σ.S)
# Access to trajectory values
diff --git a/models/ER.jl b/models/ER.jl
index 3c4a3c4..a376d6d 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -9,7 +9,7 @@ l_tr_ER = ["R1","R2","R3"]
p_ER = [1.0, 1.0, 1.0]
x0_ER = [100, 100, 0, 0]
t0_ER = 0.0
-function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int,
xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
diff --git a/models/SIR.jl b/models/SIR.jl
index c9033d4..48e9354 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -9,7 +9,7 @@ l_tr_SIR = ["R1","R2"]
p_SIR = [0.0012, 0.05]
x0_SIR = [95, 5, 0]
t0_SIR = 0.0
-function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int,
xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
diff --git a/models/_bench_perf_test/ER_col_buffer.jl b/models/_bench_perf_test/ER_col_buffer.jl
index 4261271..134c218 100644
--- a/models/_bench_perf_test/ER_col_buffer.jl
+++ b/models/_bench_perf_test/ER_col_buffer.jl
@@ -9,7 +9,7 @@ l_tr = ["R1","R2","R3"]
p = [1.0, 1.0, 1.0]
x0 = [100, 100, 0, 0]
t0 = 0.0
-function ER_col_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+function ER_col_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int,
xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
diff --git a/models/_bench_perf_test/ER_row_buffer.jl b/models/_bench_perf_test/ER_row_buffer.jl
index 0fb918f..37fbd15 100644
--- a/models/_bench_perf_test/ER_row_buffer.jl
+++ b/models/_bench_perf_test/ER_row_buffer.jl
@@ -9,7 +9,7 @@ l_tr = ["R1","R2","R3"]
p = [1.0, 1.0, 1.0]
x0 = [100, 100, 0, 0]
t0 = 0.0
-function ER_row_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx::Int,
+function ER_row_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int,
xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
diff --git a/models/_bench_perf_test/SIR_col_buffer.jl b/models/_bench_perf_test/SIR_col_buffer.jl
index 771149d..8acc8d7 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 SIR_col_buffer_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{Union{Nothing,String}}, idx::Int,
xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
diff --git a/models/_bench_perf_test/SIR_row_buffer.jl b/models/_bench_perf_test/SIR_row_buffer.jl
index ee241e6..930f1c9 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 SIR_row_buffer_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{Union{Nothing,String}}, idx::Int,
xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
diff --git a/tests/unit/check_trajectory_consistency.jl b/tests/unit/check_trajectory_consistency.jl
index 96ff4b4..33e4686 100644
--- a/tests/unit/check_trajectory_consistency.jl
+++ b/tests/unit/check_trajectory_consistency.jl
@@ -7,38 +7,87 @@ load_model("ER")
test_all = true
nb_sim = 1000
+function show_traj(io::IO, σ::AbstractTrajectory, m::Model)
+ println(io, "length(σ.values), length(σ.times), length(σ.transitions)")
+ println(io, "$(length(σ.values)), $(length(σ.times)), $(length(σ.transitions))")
+ println(io, "isbounded(m), isbounded(σ)")
+ println(io, "$(isbounded(m)), $(isbounded(σ))")
+ println(io, σ.values)
+ println(io, "$(σ.values)")
+ println(io, σ.times)
+ println(io, "$(σ.times)")
+ println(io, "σ.transitions")
+ println(io, "$(σ.transitions)")
+end
+
+σ_dump = simulate(SIR)
+σ2_dump = simulate(ER)
+
for i = 1:nb_sim
σ = simulate(SIR)
σ2 = simulate(ER)
try
- global test_all = test_all && MarkovProcesses.check_consistency(σ) && MarkovProcesses.check_consistency(σ2)
+ global test_all = test_all && check_consistency(σ) && check_consistency(σ2) &&
+ !isbounded(σ) && !isbounded(σ2)
catch err
- @show err
- @show length(σ.values), length(σ.times), length(σ.transitions)
- println()
- @show σ.values
- @show σ.times
- @show σ.transitions
+ show_traj(stdout,σ, SIR)
+ show_traj(stdout,σ2, ER)
+ global σ_dump = σ
+ global σ2_dump = σ2
throw(err)
end
+ if !test_all
+ show_traj(stdout,σ, SIR)
+ show_traj(stdout,σ2, ER)
+ global σ_dump = σ
+ global σ2_dump = σ2
+ error("Ouch")
+ end
end
-SIR.time_bound = 100.0
-ER.time_bound = 10.0
+
+SIR.time_bound = 1.0
+ER.time_bound = 0.01
for i = 1:nb_sim
σ = simulate(SIR)
σ2 = simulate(ER)
try
- global test_all = test_all && MarkovProcesses.check_consistency(σ) && MarkovProcesses.check_consistency(σ2)
+ global test_all = test_all && check_consistency(σ) && check_consistency(σ2) &&
+ isbounded(σ) && isbounded(σ2)
catch err
- @show err
- @show length(σ.values), length(σ.times), length(σ.transitions)
- println()
- @show σ.values
- @show σ.times
- @show σ.transitions
+ show_traj(stdout, σ, SIR)
+ show_traj(stdout, σ2, ER)
+ global σ_dump = σ
+ global σ2_dump = σ2
throw(err)
end
+ if !test_all
+ show_traj(stdout, σ, SIR)
+ show_traj(stdout, σ2, ER)
+ global σ_dump = σ
+ global σ2_dump = σ2
+ error("Ouch")
+ end
end
+SIR.time_bound = Inf
+ER.time_bound = Inf
+new_x0_SIR = view(reshape([95, 0, 5], 1, :), 1, :)
+new_x0_ER = view(reshape([0, 0, 0, 100], 1, :), 1, :)
+SIR.x0 = new_x0_SIR
+ER.x0 = new_x0_ER
+σ = simulate(SIR)
+σ2 = simulate(ER)
+test_all = test_all && SIR.isabsorbing(SIR.p,new_x0_SIR) && ER.isabsorbing(ER.p,new_x0_ER) &&
+ length_states(σ) == 1 && length_states(σ2) == 1 &&
+ check_consistency(σ) && check_consistency(σ2)
+
+SIR.time_bound = 1.0
+ER.time_bound = 0.01
+σ = simulate(SIR)
+σ2 = simulate(ER)
+test_all = test_all && SIR.isabsorbing(SIR.p,new_x0_SIR) && ER.isabsorbing(ER.p,new_x0_ER) &&
+ length_states(σ) == 2 && length_states(σ2) == 2 &&
+ check_consistency(σ) && check_consistency(σ2)
+
return test_all
--
GitLab