Commit 618d7910 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Try to optimize the allocation by rewriting simulate function. There's a

slight improvement in allocation and time in benchmark but too small.
parent 35ac8882
......@@ -29,7 +29,7 @@ println("MarkovProcesses:")
using MarkovProcesses
MarkovProcesses.load_model("ER")
ER.time_bound = 20.0
ER.estim_min_states = 7000
ER.estim_min_states = 9000
set_param!(ER, "k1", 0.2)
set_param!(ER, "k2", 40.0)
@timev MarkovProcesses.simulate(ER)
......
......@@ -108,7 +108,7 @@ function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
# at time n+1
mat_x = zeros(Int, m.d, buffer_size)
l_t = zeros(Float64, buffer_size)
l_tr = Vector{String}(undef, buffer_size)
l_tr = Vector{Union{String,Nothing}}(undef, buffer_size)
isabsorbing = m.isabsorbing(m.p,xn)::Bool
while !isabsorbing && (tn <= m.time_bound)
i = 0
......@@ -148,7 +148,7 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
# at time n+1
mat_x = zeros(Int, buffer_size, m.d)
l_t = zeros(Float64, buffer_size)
l_tr = Vector{String}(undef, buffer_size)
l_tr = Vector{Union{String,Nothing}}(undef, buffer_size)
isabsorbing = m.isabsorbing(m.p,xn)::Bool
while !isabsorbing && (tn <= m.time_bound)
i = 0
......@@ -189,7 +189,7 @@ function _simulate_without_view(m::ContinuousTimeModel)
# 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)
l_tr = Vector{Union{String,Nothing}}(undef, m.buffer_size)
isabsorbing = m.isabsorbing(m.p,xn)::Bool
while !isabsorbing && (tn <= m.time_bound)
i = 0
......@@ -235,7 +235,7 @@ function _simulate_27d56(m::ContinuousTimeModel)
# 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)
l_tr = Vector{Union{String,Nothing}}(undef, m.buffer_size)
isabsorbing::Bool = m.isabsorbing(m.p,xn)
# use sizehint! ?
while !isabsorbing && tn <= m.time_bound
......@@ -274,3 +274,126 @@ function _simulate_27d56(m::ContinuousTimeModel)
end
function _simulate_d7458(m::ContinuousTimeModel)
# trajectory fields
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
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{Union{String,Nothing}}(undef, m.buffer_size)
isabsorbing::Bool = m.isabsorbing(m.p,xn)
end_idx = -1
# use sizehint! ?
while !isabsorbing && tn <= m.time_bound
for i = 1:m.buffer_size
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
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))
end
append!(times, view(l_t, 1:end_idx))
append!(transitions, view(l_tr, 1:end_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)
end
values = full_values[m._g_idx]
return Trajectory(m, values, times, transitions)
end
function _simulate_last(product::SynchronizedModel)
# trajectory fields
m = product.m
A = product.automaton
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 = Union{String,Nothing}[nothing]
reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
S0 = init_state(A, m.x0, m.t0)
# values at time n
n = 0
xn = reshaped_x0
tn = m.t0
Sn = copy(S0)
# 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)
isacceptedLHA::Bool = isaccepted(Sn)
Snplus1 = copy(Sn)
while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
i = 0
while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
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, :)
tr_n = l_tr[i]
next_state!(Snplus1, A, xn, tn, tr_n, Sn)
Sn = Snplus1
isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1)
end
for k = 1:m.d
append!(full_values[k], view(mat_x, 1:i, k))
end
append!(times, view(l_t, 1:i))
append!(transitions, view(l_tr, 1:i))
n += i
end
# When the trajectory is accepted, we should not add an end value
if isbounded(m) && !isaccepted(Sn)
@assert times[end] < m.time_bound
for k = 1:m.d
push!(full_values[k], full_values[k][end])
end
push!(times, m.time_bound)
push!(transitions, nothing)
end
values = full_values[m._g_idx]
return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
end
......@@ -16,16 +16,21 @@ function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{
push!(times, time_bound)
push!(transitions, nothing)
end
"""
`simulate(m)`
# Simulation
Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a
Linear Hybrid Automaton.
"""
function simulate(m::ContinuousTimeModel)
# First alloc
full_values = Vector{Vector{Int}}(undef, m.d)
for i = 1:m.d full_values[i] = zeros(Int, m.estim_min_states) end
for i = eachindex(full_values) 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
for i = eachindex(full_values) full_values[i][1] = m.x0[i] end
times[1] = m.t0
transitions[1] = nothing
# Values at time n
......@@ -34,6 +39,7 @@ function simulate(m::ContinuousTimeModel)
tn = m.t0
# at time n+1
isabsorbing::Bool = m.isabsorbing(m.p,xn)
# If x0 is absorbing
if isabsorbing
_resize_trajectory!(full_values, times, transitions, 1)
values = full_values[m._g_idx]
......@@ -55,7 +61,7 @@ function simulate(m::ContinuousTimeModel)
n += 1
xn = view(mat_x, 1, :)
# Updating value
for k = 1:m.d full_values[k][n] = xn[k] end
for k = eachindex(full_values) full_values[k][n] = xn[k] end
times[n] = tn
transitions[n] = l_tr[1]
isabsorbing = m.isabsorbing(m.p,xn)
......@@ -63,7 +69,7 @@ function simulate(m::ContinuousTimeModel)
break
end
end
# If simulation ended
# If simulation ended before the estimation of states
if n < m.estim_min_states
_resize_trajectory!(full_values, times, transitions, n)
values = full_values[m._g_idx]
......@@ -76,7 +82,14 @@ function simulate(m::ContinuousTimeModel)
mat_x = zeros(Int, m.buffer_size, m.d)
l_t = zeros(Float64, m.buffer_size)
l_tr = Vector{Transition}(undef, m.buffer_size)
# Alloc buffer values
tmp_full_values = Vector{Vector{Int}}(undef, m.d)
for i = eachindex(tmp_full_values) tmp_full_values[i] = zeros(Int, 0) end
tmp_times = zeros(Float64, 0)
tmp_transitions = Vector{Transition}(undef, 0)
tmp_idx = 0
while !isabsorbing && tn <= m.time_bound
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+m.buffer_size)
i = 0
while i < m.buffer_size
i += 1
......@@ -92,15 +105,22 @@ function simulate(m::ContinuousTimeModel)
break
end
end
_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]
rng_tmp = (tmp_idx+1):(tmp_idx+m.buffer_size)
for k = eachindex(tmp_full_values) tmp_full_values[k][rng_tmp] = view(mat_x, :, k) end
tmp_times[rng_tmp] = l_t
tmp_transitions[rng_tmp] = l_tr
if i < m.buffer_size
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+i)
end
tmp_idx += i
n += i
end
# Push the temporary values
for k = eachindex(full_values) append!(full_values[k], tmp_full_values[k]) end
append!(times, tmp_times)
append!(transitions, tmp_transitions)
values = full_values[m._g_idx]
if isbounded(m)
# Add last value: the convention is the last transition is nothing,
......@@ -109,38 +129,89 @@ function simulate(m::ContinuousTimeModel)
end
return Trajectory(m, values, times, transitions)
end
function simulate(product::SynchronizedModel)
# trajectory fields
m = product.m
A = product.automaton
# 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 = Union{String,Nothing}[nothing]
reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
for i = eachindex(full_values) 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 = eachindex(full_values) full_values[i][1] = m.x0[i] end
times[1] = m.t0
transitions[1] = nothing
S0 = init_state(A, m.x0, m.t0)
# values at time n
n = 0
xn = reshaped_x0
tn = m.t0
# Values at time n
n = 1
xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
tn = m.t0
Sn = copy(S0)
# 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)
isacceptedLHA::Bool = isaccepted(Sn)
if isabsorbing || isacceptedLHA
_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 SynchronizedTrajectory(Sn, product, 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]
Snplus1 = copy(Sn)
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, :)
tr_n = l_tr[1]
next_state!(Snplus1, A, xn, tn, tr_n, Sn)
Sn = Snplus1
# Updating value
for k = eachindex(full_values) full_values[k][n] = xn[k] end
times[n] = tn
transitions[n] = l_tr[1]
isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA
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 SynchronizedTrajectory(Snplus1, product, 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{Transition}(undef, m.buffer_size)
# Alloc buffer values
tmp_full_values = Vector{Vector{Int}}(undef, m.d)
for i = eachindex(tmp_full_values) tmp_full_values[i] = zeros(Int, 0) end
tmp_times = zeros(Float64, 0)
tmp_transitions = Vector{Transition}(undef, 0)
tmp_idx = 0
while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+m.buffer_size)
i = 0
while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
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
i -= 1
break
end
xn = view(mat_x, i, :)
......@@ -149,27 +220,36 @@ function simulate(product::SynchronizedModel)
Sn = Snplus1
isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA
break
end
end
for k = 1:m.d
append!(full_values[k], view(mat_x, 1:i, k))
# Update values
rng_tmp = (tmp_idx+1):(tmp_idx+m.buffer_size)
for k = eachindex(tmp_full_values) tmp_full_values[k][rng_tmp] = view(mat_x, :, k) end
tmp_times[rng_tmp] = l_t
tmp_transitions[rng_tmp] = l_tr
if i < m.buffer_size
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+i)
end
append!(times, view(l_t, 1:i))
append!(transitions, view(l_tr, 1:i))
tmp_idx += i
n += i
end
# When the trajectory is accepted, we should not add an end value
# Push the temporary values
for k = eachindex(full_values) append!(full_values[k], tmp_full_values[k]) end
append!(times, tmp_times)
append!(transitions, tmp_transitions)
values = full_values[m._g_idx]
if isbounded(m) && !isaccepted(Sn)
@assert times[end] < m.time_bound
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 SynchronizedTrajectory(Snplus1, product, values, times, transitions)
end
function simulate(m::ContinuousTimeModel, n::Int)
end
......
......@@ -10,7 +10,7 @@ p = [1.0, 1.0, 1.0]
x0 = [100, 100, 0, 0]
t0 = 0.0
function ER_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
function ER_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{Union{Nothing,String}},
xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
......
......@@ -9,7 +9,7 @@ l_tr = ["R1","R2"]
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},
function SIR_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{Union{Nothing,String}},
xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
......
......@@ -12,6 +12,16 @@ sync_SIR = A_F * SIR
function test_last_state(A::LHA, m::SynchronizedModel)
σ = simulate(m)
test = (get_state_from_time(σ, (σ.S).time)[1] == (σ.S)["n"]) && ((σ.S)["d"] == 0)
if !test
@show σ.S
@show times(σ)[end]
@show σ[end]
@show times(σ)[end-1]
@show σ[end-1]
@show get_state_from_time(σ, (σ.S).time)[1]
@show (σ.S)["n"], (σ.S)["d"]
error("Ouch")
end
return test
end
......
......@@ -23,5 +23,6 @@ using Test
@test include("unit/create_models.jl")
@test include("unit/create_automata.jl")
@test include("unit/model_prior.jl")
@test include("unit/absorbing_x0.jl")
end
using MarkovProcesses
load_model("SIR")
new_x0 = [95, 0, 0]
SIR.x0 = new_x0
σ = simulate(SIR)
check_consistency(σ)
return length_states(σ) == 1 && length(times(σ)) == 1 && length(transitions(σ)) == 1
......@@ -8,8 +8,8 @@ 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, "length(σ.values[1]), length(σ.times), length(σ.transitions)")
println(io, "$(length(σ.values[1])), l$(length(σ.times)), $(length(σ.transitions))")
println(io, "isbounded(m), isbounded(σ)")
println(io, "$(isbounded(m)), $(isbounded(σ))")
println(io, σ.values)
......
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