diff --git a/bench/pygmalion/multiple_long_sim_er.jl b/bench/pygmalion/multiple_long_sim_er.jl
index 6729aba714497f3668cacdb509bf07d14800709f..23244dc0cc06ccb13f35fcbe9b949ad03b913689 100644
--- a/bench/pygmalion/multiple_long_sim_er.jl
+++ b/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)
 
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 5fb0ec4540c0c6e4ca0ac7eff053e5d72f4e2f54..85eb0f67deab563d5ee0533be257099f67e50b44 100644
--- a/core/MarkovProcesses.jl
+++ b/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
 
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index 632b3ffdafee62b9136cd6b12806cf5687aed079..9c40acbbaeca5660e8c50245584be07c29f11c7d 100644
--- a/core/_tests_simulate.jl
+++ b/core/_tests_simulate.jl
@@ -1,21 +1,28 @@
 
+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
+
+
diff --git a/core/common.jl b/core/common.jl
index bbe0914fd88844f07af51e5b0ab245e40ec42aa8..8f7b015d80612439b84a0d2751a891a12c240240 100644
--- a/core/common.jl
+++ b/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
diff --git a/core/lha.jl b/core/lha.jl
index c47b7b6a939c6bb2cb4e02bdebbf52a325515843..2445b23d076bbcb53096adf10f213ba725ee0cd0 100644
--- a/core/lha.jl
+++ b/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
diff --git a/core/trajectory.jl b/core/trajectory.jl
index cb8901babed4c938fb920864634443458b171c14..a02735cf0fbe9a59b3372065558734c307638ebd 100644
--- a/core/trajectory.jl
+++ b/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
diff --git a/tests/dist_lp/dist_case_2.jl b/tests/dist_lp/dist_case_2.jl
index be7c86d801b4bc071fe0d73187f29b2d191ba8b6..c66ce3b1dfc4dd66b67339fb312cc24f7eac78be 100644
--- a/tests/dist_lp/dist_case_2.jl
+++ b/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)
diff --git a/tests/dist_lp/dist_case_3.jl b/tests/dist_lp/dist_case_3.jl
index b27177121813e71a7453e8778311ef317c1d57f6..8254be090ae5715242245eaa3b1ecc906bcf0ff9 100644
--- a/tests/dist_lp/dist_case_3.jl
+++ b/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)
diff --git a/tests/dist_lp/dist_l1_case_1.jl b/tests/dist_lp/dist_l1_case_1.jl
index a9152e0436815d868d1fcce0674f10647c64f3a6..262c8c77c0140bbf8e23d71540f73fdbf4a412ab 100644
--- a/tests/dist_lp/dist_l1_case_1.jl
+++ b/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
diff --git a/tests/unit/check_trajectory_consistency.jl b/tests/unit/check_trajectory_consistency.jl
index 4ae5da9bd52213f41d913f47fe2b5f447ea92b80..96ff4b44641acdbc0f7f2736d753aa3648fda946 100644
--- a/tests/unit/check_trajectory_consistency.jl
+++ b/tests/unit/check_trajectory_consistency.jl
@@ -1,13 +1,29 @@
 
 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)