Change of Trajectory.values type from Matrix to Vector of Vector in

order to replace vcat by push!.

The benchmarks are at least equal, and better when the buffer size is
small. However, no significative difference of performance with adequate
buffer size.

My opinion the gain can be seen when simulations are going to be very
long with more complicated models than we have implemented for now.

All tests passed.

bench/pygmalion/multiple_long_sim_er.jl

@@ -38,8 +38,8 @@ b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
 @show minimum(b2), mean(b2), maximum(b2)
 println("Buffer size 100")
 ER.buffer_size = 100
-b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
-b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
+b1_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
+b2_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
 @show minimum(b1), mean(b1), maximum(b1)
 @show minimum(b2), mean(b2), maximum(b2)
 

core/MarkovProcesses.jl

@@ -22,7 +22,7 @@ export load_automaton, get_index, get_value, length_var, isaccepted
 
 # Model related methods
 export simulate, set_param!, get_param, set_observed_var!
-export set_time_bound!
+export set_time_bound!, getproperty
 export isbounded, isaccepted, check_consistency
 export load_model, get_module_path
 

core/_tests_simulate.jl

 
+struct OldTrajectory <: AbstractTrajectory 
+    m::ContinuousTimeModel
+    values::Matrix{Int}
+    times::Vector{Float64}
+    transitions::Vector{Union{String,Nothing}}
+end 
+
 # File for benchmarking simulation and memory access of the package.
 
 # Trajectories
 
-_get_values_col(σ::AbstractTrajectory, var::String) = 
+_get_values_col(σ::OldTrajectory, var::String) = 
 @view σ.values[(σ.m)._map_obs_var_idx[var],:] 
-_get_values_row(σ::AbstractTrajectory, var::String) = 
+_get_values_row(σ::OldTrajectory, var::String) = 
 @view σ.values[:,(σ.m)._map_obs_var_idx[var]] 
 
-_get_state_col(σ::AbstractTrajectory, idx::Int) = 
+_get_state_col(σ::OldTrajectory, idx::Int) = 
 @view σ.values[:,idx]
-_get_state_row(σ::AbstractTrajectory, idx::Int) = 
+_get_state_row(σ::OldTrajectory, idx::Int) = 
 @view σ.values[idx,:]
 
-_get_value_col(σ::AbstractTrajectory, var::String, idx::Int) = 
+_get_value_col(σ::OldTrajectory, var::String, idx::Int) = 
 σ.values[(σ.m)._map_obs_var_idx[var],idx] 
-_get_value_row(σ::AbstractTrajectory, var::String, idx::Int) = 
+_get_value_row(σ::OldTrajectory, var::String, idx::Int) = 
 σ.values[idx,(σ.m)._map_obs_var_idx[var]] 
 
 # Model
@@ -51,7 +58,7 @@ function _simulate_col(m::ContinuousTimeModel)
             transitions[end] = nothing
         end
     end
-    return Trajectory(m, values, times, transitions)
+    return OldTrajectory(m, values, times, transitions)
 end
 
 function _simulate_row(m::ContinuousTimeModel)
@@ -85,7 +92,7 @@ function _simulate_row(m::ContinuousTimeModel)
             transitions[end] = nothing
         end
     end
-    return Trajectory(m, values, times, transitions)
+    return OldTrajectory(m, values, times, transitions)
 end
 
 
@@ -126,7 +133,7 @@ function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
             transitions[end] = nothing
         end
     end
-    return Trajectory(m, values, times, transitions)
+    return OldTrajectory(m, values, times, transitions)
 end
 
 function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
@@ -166,7 +173,7 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
             transitions[end] = nothing
         end
     end
-    return Trajectory(m, values, times, transitions)
+    return OldTrajectory(m, values, times, transitions)
 end
 
 function _simulate_without_view(m::ContinuousTimeModel)
@@ -211,6 +218,59 @@ function _simulate_without_view(m::ContinuousTimeModel)
         end
     end
     values = full_values[:,m._g_idx] 
-    return Trajectory(m, values, times, transitions)
+    return OldTrajectory(m, values, times, transitions)
 end
 
+# With trajectory values in Matrix type
+function _simulate_27d56(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 = 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{String}(undef, m.buffer_size)
+    isabsorbing::Bool = m.isabsorbing(m.p,xn)
+    # use sizehint! ?
+    while !isabsorbing && tn <= m.time_bound
+        i = 0
+        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound
+            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, :)
+            isabsorbing = m.isabsorbing(m.p,xn)
+        end
+        full_values = vcat(full_values, view(mat_x, 1:i, :))
+        append!(times, view(l_t, 1:i))
+        append!(transitions,  view(l_tr, 1:i))
+        n += i
+    end
+    if isbounded(m)
+        #=
+        if times[end] > m.time_bound
+            full_values[end,:] = full_values[end-1,:]
+            times[end] = m.time_bound
+            transitions[end] = nothing
+        else
+        end
+        =#
+        full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
+        push!(times, m.time_bound)
+        push!(transitions, nothing)
+    end
+    values = view(full_values, :, m._g_idx)
+    return OldTrajectory(m, values, times, transitions)
+end
+
+

core/common.jl

@@ -26,7 +26,7 @@ end
 
 struct Trajectory <: AbstractTrajectory
     m::ContinuousTimeModel
-    values::Matrix{Int}
+    values::Vector{Vector{Int}}
     times::Vector{Float64}
     transitions::Vector{Transition}
 end
@@ -65,7 +65,7 @@ end
 struct SynchronizedTrajectory <: AbstractTrajectory
     S::StateLHA
     m::SynchronizedModel
-    values::Matrix{Int}
+    values::Vector{Vector{Int}}
     times::Vector{Float64}
     transitions::Vector{Transition}
 end

core/lha.jl

@@ -24,7 +24,7 @@ end
 isaccepted(S::StateLHA) = (S.loc in (S.A).l_loc_final)
 
 # Methods for synchronize / read the trajectory
-function init_state(A::LHA, x0::SubArray{Int,1}, t0::Float64)
+function init_state(A::LHA, x0::Vector{Int}, t0::Float64)
     S0 = StateLHA(A, "", zeros(length_var(A)), t0)
     for loc in A.l_loc_init
         if A.Λ[loc](A,S0) 
@@ -119,12 +119,17 @@ function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
     A_new = LHA(A, (σ.m)._map_obs_var_idx)
     l_t = times(σ)
     l_tr = transitions(σ)
+    mat_x = zeros(Int, length_states(σ), σ.m.d)
+    for (i,var) in enumerate(σ.m.g)
+        mat_x[:,i] = σ[var]
+    end
     Sn = init_state(A_new, σ[1], l_t[1])
     Snplus1 = copy(Sn)
     if verbose println("Init: ") end
     if verbose @show Sn end
     for n in 1:length_states(σ)
-        next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn; verbose = verbose)
+        xn = view(mat_x, n, :)
+        next_state!(Snplus1, A_new, xn, l_t[n], l_tr[n], Sn; verbose = verbose)
         if Snplus1.loc in A_new.l_loc_final 
             break 
         end

core/trajectory.jl

@@ -64,7 +64,7 @@ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String;
 end
 
 # Distance function. Vectorized version
-function dist_lp(x_obs::SubArray{Int,1}, t_x::Vector{Float64}, y_obs::SubArray{Int,1}, t_y::Vector{Float64}; 
+function dist_lp(x_obs::Vector{Int}, t_x::Vector{Float64}, y_obs::Vector{Int}, t_y::Vector{Float64}; 
                  verbose::Bool = false, p::Int = 1)
     current_y_obs = y_obs[1]
     current_t_y = t_y[2]
@@ -150,23 +150,31 @@ function _riemann_sum(f::Function, t_begin::Real, t_end::Real, step::Float64)
 end
 
 function check_consistency(σ::AbstractTrajectory)
-    @assert length(σ.times) == length(σ.transitions) == size(σ.values)[1]
-    @assert length_obs_var(σ) == length(σ.m.g) == size(σ.values)[2]
+    test_length_var = true
+    for i = 1:(σ.m).dobs 
+        test_length_i = (length(σ.values[1]) == length(σ.values[i]))
+        test_length_var = test_length_var && test_length_i 
+    end
+    @assert begin (length(σ.times) == length(σ.transitions)) && 
+                  (length(σ.times) == length(σ.values[1])) &&
+                  test_length_var
+            end
+    @assert length_obs_var(σ) == σ.m.dobs
     return true
 end
-issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, σ[end])
+issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, view(reshape(σ[end], 1, m.d), 1, :))
 
 # Properties of the trajectory
 length_states(σ::AbstractTrajectory) = length(σ.times)
-length_obs_var(σ::AbstractTrajectory) = size(σ.values)[2]
+length_obs_var(σ::AbstractTrajectory) = length(σ.values)
 get_obs_var(σ::AbstractTrajectory) = (σ.m).g
 isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing 
 isaccepted(σ::SynchronizedTrajectory) = isaccepted(σ.S)
 
 # Access to trajectory values
-get_var_values(σ::AbstractTrajectory, var::String) = view(σ.values, :, (σ.m)._map_obs_var_idx[var])
-get_state(σ::AbstractTrajectory, idx::Int) = view(σ.values, idx, :)
-get_value(σ::AbstractTrajectory, var::String, idx::Int) = σ.values[idx,(σ.m)._map_obs_var_idx[var]] 
+get_var_values(σ::AbstractTrajectory, var::String) = σ.values[(σ.m)._map_obs_var_idx[var]]
+get_state(σ::AbstractTrajectory, idx::Int) = [σ.values[i][idx] for i = 1:length(σ.values)] # /!\ Creates an array
+get_value(σ::AbstractTrajectory, var::String, idx::Int) = get_var_values(σ, var)[idx]
 # Operation @
 function get_state_from_time(σ::AbstractTrajectory, t::Float64)
     @assert t >= 0.0

tests/dist_lp/dist_case_2.jl

@@ -11,21 +11,21 @@ for p = 1:2
 
     y_obs = [1, 2, 5, 3, 3]
     t_y = [0.0, 0.5, 0.7, 1.4, 3.0]
-    values = zeros(length(y_obs), 1)
-    values[:,1] = y_obs
+    values = [zeros(length(y_obs))]
+    values[1] = y_obs
     l_tr = Vector{Nothing}(nothing, length(y_obs))
     σ2 = Trajectory(SIR, values, t_y, l_tr)
 
     x_obs = [4, 2, 3, 1, 1]
     t_x = [0.0, 0.6, 0.7, 0.8, 3.0]
-    values = zeros(length(x_obs), 1)
-    values[:,1] = x_obs
+    values = [zeros(length(x_obs))]
+    values[1] = x_obs
     l_tr = Vector{Nothing}(nothing, length(x_obs))
     σ1 = Trajectory(SIR, values, t_x, l_tr)
 
-    test_1 = isapprox(dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p), res; atol = 1E-10) 
+    test_1 = isapprox(dist_lp(x_obs, t_x, y_obs, 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(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; 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 = 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)

tests/dist_lp/dist_case_3.jl

@@ -7,15 +7,15 @@ test_all = true
 for p = 1:2
     x_obs =[5, 6, 5, 4, 3, 2, 1, 1]
     t_x = [0.0, 3.10807, 4.29827, 4.40704, 5.67024, 7.1299, 11.2763, 20.0]
-    values = zeros(length(x_obs), 1)
-    values[:,1] = x_obs
+    values = [zeros(length(x_obs))]
+    values[1] = x_obs
     l_tr = Vector{Nothing}(nothing, length(x_obs))
     σ1 = Trajectory(SIR, values, t_x, l_tr)
 
     y_obs =[5, 4, 5, 4, 3, 4, 3, 2, 3, 2, 1, 2, 3, 4, 3, 4, 4]
     t_y = [0.0, 0.334082, 1.21012, 1.40991, 1.58866, 2.45879, 2.94545, 4.66746, 5.44723, 5.88066, 7.25626, 11.4036, 13.8373, 17.1363, 17.8193, 18.7613, 20.0]
-    values = zeros(length(y_obs), 1)
-    values[:,1] = y_obs
+    values = [zeros(length(y_obs))]
+    values[1] = y_obs
     l_tr = Vector{Nothing}(nothing, length(y_obs))
     σ2 = Trajectory(SIR, values, t_y, l_tr)
 
@@ -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(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p)
+    res1 = dist_lp(x_obs, t_x, y_obs, t_y; p=p)
     res2 = dist_lp(σ1,σ2; p=p)
-    res1_bis = dist_lp(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=p)
+    res1_bis = dist_lp(y_obs, t_y, x_obs, 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)

tests/dist_lp/dist_l1_case_1.jl

@@ -8,20 +8,20 @@ load_model("SIR")
 x_obs = [2, 4, 3, 3]
 t_x = [0.0, 0.5, 1.2, 2.2]
 t_y_bis = [0.0, 0.5, 1.2, 2.2]
-values = zeros(length(x_obs), 1)
-values[:,1] = x_obs
+values = [zeros(length(x_obs))]
+values[1] = x_obs
 l_tr = Vector{Nothing}(nothing, length(x_obs))
 σ1 = Trajectory(SIR, values, t_x, l_tr)
 
 y_obs = [6, 6]
 t_y = [0.0, 2.2]
 t_x_bis = [0.0, 2.2]
-values = zeros(length(y_obs), 1)
-values[:,1] = y_obs
+values = [zeros(length(y_obs))]
+values[1] = y_obs
 l_tr = Vector{Nothing}(nothing, length(y_obs))
 σ2 = Trajectory(SIR, values, t_y, l_tr)
 
-test_1 = dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=1) == 6.4
+test_1 = dist_lp(x_obs, t_x, y_obs, 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(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=1) == 6.4 
+test_1_bis = dist_lp(y_obs, t_y, x_obs, 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

tests/unit/check_trajectory_consistency.jl

 
 using MarkovProcesses
+
 load_model("SIR")
-SIR.time_bound = 100.0
 load_model("ER")
-ER.time_bound = 10.0
 
 test_all = true
 nb_sim = 1000
 
+for i = 1:nb_sim
+    σ = simulate(SIR)
+    σ2 = simulate(ER)
+    try
+        global test_all = test_all && MarkovProcesses.check_consistency(σ) && MarkovProcesses.check_consistency(σ2)
+    catch err
+        @show err
+        @show length(σ.values), length(σ.times), length(σ.transitions)
+        println()
+        @show σ.values
+        @show σ.times
+        @show σ.transitions 
+        throw(err)
+    end
+end
+SIR.time_bound = 100.0
+ER.time_bound = 10.0
 for i = 1:nb_sim
     σ = simulate(SIR)
     σ2 = simulate(ER)