Commit 8020a182 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Fix of several type unstabilities due to @view macro.

Add of bench for @view macro.
parent 963ef13b
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)
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)
......@@ -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
......@@ -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)
......
......@@ -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
......
......@@ -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"]
......
......@@ -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)
......
......@@ -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)
......
......@@ -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)
......
......@@ -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)
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment