Change of the simulate function with resize!. Gain of performance is

small.

Fix of a misbehaviour of isbounded() for trajectories.

core/model.jl

 
 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
 

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

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]

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]

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]

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]

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]

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]

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