diff --git a/bench/pygmalion/multiple_long_sim_er.jl b/bench/pygmalion/multiple_long_sim_er.jl
index 18552e08e1a5f063f3bc75ba7980e580517e939f..f804a69355db94b91d6efb297b882587fdcff7fd 100644
--- a/bench/pygmalion/multiple_long_sim_er.jl
+++ b/bench/pygmalion/multiple_long_sim_er.jl
@@ -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)
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index 9c40acbbaeca5660e8c50245584be07c29f11c7d..cdc1ac7d6a58b55284b2b54a4d453756caf27e8c 100644
--- a/core/_tests_simulate.jl
+++ b/core/_tests_simulate.jl
@@ -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
+
+
diff --git a/core/model.jl b/core/model.jl
index a07e731490cb00858bc0584ffa23e3c0c835d269..7865a8b44e4b4cb2563131a8378285e196309cf4 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -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
diff --git a/models/_bench_perf_test/ER_col.jl b/models/_bench_perf_test/ER_col.jl
index e548336ac6a5a3da140e66ceabaf2703592d2e6b..366e368d379f7b7f58f2a38a965c0e1bdb50717c 100644
--- a/models/_bench_perf_test/ER_col.jl
+++ b/models/_bench_perf_test/ER_col.jl
@@ -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]
diff --git a/models/_bench_perf_test/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl
index 381a487f3fc6874a4fcd48012fd46fdd24ade348..cc522164b740861bf70391ca281e59c539accf34 100644
--- a/models/_bench_perf_test/SIR_col.jl
+++ b/models/_bench_perf_test/SIR_col.jl
@@ -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]
diff --git a/tests/automata/sync_simulate_last_state_F.jl b/tests/automata/sync_simulate_last_state_F.jl
index 7f65ae6f381759fc44e8c3149d0294742847c07e..5a1cf7adefaf35bd0fbbb767cf823b70aafa24ea 100644
--- a/tests/automata/sync_simulate_last_state_F.jl
+++ b/tests/automata/sync_simulate_last_state_F.jl
@@ -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
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 3202f53e9c18a73cbaa30f3594f5c5298d5026e4..a384dcb288e0879e736960ecead45c40be5dee84 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -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
diff --git a/tests/unit/absorbing_x0.jl b/tests/unit/absorbing_x0.jl
new file mode 100644
index 0000000000000000000000000000000000000000..222fe459c9032d6f55094430853b27ba819a8bb1
--- /dev/null
+++ b/tests/unit/absorbing_x0.jl
@@ -0,0 +1,11 @@
+
+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
+
diff --git a/tests/unit/check_trajectory_consistency.jl b/tests/unit/check_trajectory_consistency.jl
index 33e468615a795e744ac7cc0e9a7099c7185ed1d3..f60a07a0ae4b77a031ffe750aaa7d59e8573ce8c 100644
--- a/tests/unit/check_trajectory_consistency.jl
+++ b/tests/unit/check_trajectory_consistency.jl
@@ -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)