diff --git a/bench/pygmalion/multiple_long_sim_er.jl b/bench/pygmalion/multiple_long_sim_er.jl
index 3e27d6a2c97d25635e8b3af45877141a1a326220..6729aba714497f3668cacdb509bf07d14800709f 100644
--- a/bench/pygmalion/multiple_long_sim_er.jl
+++ b/bench/pygmalion/multiple_long_sim_er.jl
@@ -31,6 +31,13 @@ MarkovProcesses.load_model("ER")
 ER.time_bound = 20.0
 @timev MarkovProcesses.simulate(ER)
 
+println("Default buffer size=10")
+b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
+b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
+@show minimum(b1), mean(b1), maximum(b1)
+@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
 @show minimum(b1), mean(b1), maximum(b1)
diff --git a/bench/pygmalion/read_trajectory_er.jl b/bench/pygmalion/read_trajectory_er.jl
index 36ec8cce500b156dc3b06d3eaeb19345e476d499..f40df1dc3273505f9a7c44a4545aad37fefef132 100644
--- a/bench/pygmalion/read_trajectory_er.jl
+++ b/bench/pygmalion/read_trajectory_er.jl
@@ -16,7 +16,7 @@ u = Control(10.0)
 tml = 1:400
 g_all = create_observation_function([ObserverModel(str_oml, tml)]) 
 so = pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
-function read_trajectory_v1(so::SystemObservation)
+function read_trajectory_pyg_v1(so::SystemObservation)
     vals = to_vec(so, "P")
     n_states = length(vals)
     res = 0.0
@@ -25,7 +25,7 @@ function read_trajectory_v1(so::SystemObservation)
     end
     return res
 end
-function read_trajectory_v2(so::SystemObservation)
+function read_trajectory_pyg_v2(so::SystemObservation)
     idx_P = om_findfirst("P", so.oml)
     n_states = length(so.otll[idx_P]) 
     res = 0.0
@@ -35,11 +35,11 @@ function read_trajectory_v2(so::SystemObservation)
     return res
 end
 # Bench
-@timev read_trajectory_v1(so)
-b1_pyg = @benchmark read_trajectory_v1($so) 
+@timev read_trajectory_pyg_v1(so)
+b1_pyg = @benchmark read_trajectory_pyg_v1($so) 
 @show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
-@timev read_trajectory_v2(so)
-b2_pyg = @benchmark read_trajectory_v2($so) 
+@timev read_trajectory_pyg_v2(so)
+b2_pyg = @benchmark read_trajectory_pyg_v2($so) 
 @show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
 
 println("MarkovProcesses:")
@@ -47,16 +47,28 @@ println("MarkovProcesses:")
 MarkovProcesses.load_model("ER")
 ER.time_bound = 10.0
 σ = MarkovProcesses.simulate(ER)
-function read_trajectory(σ::AbstractTrajectory)
+function read_trajectory_v1(σ::AbstractTrajectory)
     n_states = length_states(σ)
     res = 0
     for i = 1:n_states
-        res += (σ["P"])[i]
+        res += σ["P",i]
+    end
+    return res
+end
+function read_trajectory_v2(σ::AbstractTrajectory)
+    n_states = length_states(σ)
+    vals = σ["P"]
+    res = 0
+    for i = 1:n_states
+        res += vals[i]
     end
     return res
 end
 # Bench
-@timev read_trajectory(σ) 
-b1 = @benchmark read_trajectory($σ)
+@timev read_trajectory_v1(σ) 
+b1 = @benchmark read_trajectory_v1($σ)
 @show minimum(b1), mean(b1), maximum(b1)
+@timev read_trajectory_v2(σ) 
+b2 = @benchmark read_trajectory_v2($σ)
+@show minimum(b2), mean(b2), maximum(b2)
 
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 029af35579835f397f3a304c30bf4371c47ee1f2..e9d26143b64737cb6bfd92bb7230dffc295c5ba6 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -1,6 +1,6 @@
 module MarkovProcesses
 
-import Base: +, -, getfield, getindex
+import Base: +, -, getfield, getindex, ImmutableDict
 
 export Model, ContinuousTimeModel, DiscreteTimeModel
 export simulate, set_param!, get_param, set_observed_var!
@@ -10,7 +10,7 @@ include("model.jl")
 
 export Observations, AbstractTrajectory, Trajectory
 export +, -, δ, dist_lp
-export get_obs_var, length_states, length_obs_var, is_bounded
+export get_obs_var, length_states, length_obs_var, is_bounded, times, transitions
 include("trajectory.jl")
 
 end
diff --git a/core/model.jl b/core/model.jl
index 541a4bbc4aa771180e9649deb8be99e6ad261ef7..e2b23b07b0b24c4bac76e3449a3cc51f3a47b75d 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -2,15 +2,14 @@
 import StaticArrays: SVector
 
 abstract type Model end 
-abstract type ContinuousTimeModel <: Model end 
 abstract type DiscreteTimeModel <: Model end 
 
-mutable struct CTMC <: ContinuousTimeModel
+mutable struct ContinuousTimeModel <: Model
     d::Int # state space dim
     k::Int # parameter space dim
-    map_var_idx::Dict # maps str to full state space
-    _map_obs_var_idx::Dict # maps str to observed state space
-    map_param_idx::Dict # maps str in parameter space
+    map_var_idx::Dict{String,Int} # maps str to full state space
+    _map_obs_var_idx::Dict{String,Int} # maps str to observed state space
+    map_param_idx::Dict{String,Int} # maps str in parameter space
     l_name_transitions::Vector{String}
     p::Vector{Float64}
     x0::Vector{Int}
@@ -23,7 +22,7 @@ mutable struct CTMC <: ContinuousTimeModel
     buffer_size::Int
 end
 
-function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::Vector{String}, 
+function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::Vector{String}, 
               p::Vector{Float64}, x0::Vector{Int}, t0::Float64, 
               f!::Function, is_absorbing::Function; 
               g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10)
@@ -42,7 +41,7 @@ function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_tra
         @warn "You have possibly redefined a function Model.is_absorbing used in a previously instantiated model."
     end
 
-    return CTMC(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
+    return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
 end
 
 function simulate(m::ContinuousTimeModel)
@@ -53,27 +52,27 @@ function simulate(m::ContinuousTimeModel)
     transitions = Union{String,Nothing}[nothing]
     # values at time n
     n = 0
-    xn = @view m.x0[:] # View for type stability
+    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)
-    is_absorbing = m.is_absorbing(m.p,xn)::Bool
+    is_absorbing::Bool = m.is_absorbing(m.p,xn)
     while !is_absorbing && (tn <= m.time_bound)
         i = 0
         while i < m.buffer_size && !is_absorbing && (tn <= m.time_bound)
             i += 1
             m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
-            xn = @view mat_x[i,:]
+            xn = view(mat_x, i, :)
             tn = l_t[i]
-            is_absorbing = m.is_absorbing(m.p,xn)::Bool
+            is_absorbing = m.is_absorbing(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])
+        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
-        is_absorbing = m.is_absorbing(m.p,xn)::Bool
+        is_absorbing = m.is_absorbing(m.p,xn)
     end
     if is_bounded(m)
         if times[end] > m.time_bound
@@ -86,7 +85,7 @@ function simulate(m::ContinuousTimeModel)
             push!(transitions, nothing)
         end
     end
-    values = @view full_values[:,m._g_idx] 
+    values = view(full_values, :, m._g_idx)
     return Trajectory(m, values, times, transitions)
 end
 
diff --git a/core/trajectory.jl b/core/trajectory.jl
index 5ef494095e510c410d60a1ea62c600a2bbaba2bf..6f8ce4ad82bd4a37f03a74d3ef35e54a6b54ef87 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -74,11 +74,11 @@ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String;
     if !is_bounded(σ1) || !is_bounded(σ2)
         @warn "Lp distance computed on unbounded trajectories. Result should be wrong"
     end
-    return dist_lp(σ1[var], σ1["times"], σ2[var], σ2["times"]; verbose = false, p = p)
+    return dist_lp(σ1[var], times(σ1), σ2[var], times(σ2); verbose = false, p = p)
 end
 
 # Distance function. Vectorized version
-function dist_lp(x_obs::SubArray{Int,1}, t_x::SubArray{Float64,1}, y_obs::SubArray{Int,1}, t_y::SubArray{Float64,1}; 
+function dist_lp(x_obs::SubArray{Int,1}, t_x::Vector{Float64}, y_obs::SubArray{Int,1}, t_y::Vector{Float64}; 
                  verbose::Bool = false, p::Int = 1)
     current_y_obs = y_obs[1]
     current_t_y = t_y[2]
@@ -179,32 +179,23 @@ is_bounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing
 
 # 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,:]
+    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]] 
-function δ(σ1::AbstractTrajectory,t::Float64) end
-
-# Get var values ["I"]
-function getindex(σ::AbstractTrajectory, var::String)
-    if var  == "times"
-        return @view σ.times[:]
-    elseif var == "transitions"
-        return @view σ.transitions[:] 
-    else
-        return get_var_values(σ, var)
-    end
-end
+    σ.values[idx,(σ.m)._map_obs_var_idx[var]] 
+
+
+δ(σ::AbstractTrajectory,t::Float64) = true
+
+states(σ::AbstractTrajectory) = σ.values
+times(σ::AbstractTrajectory) = σ.times
+transitions(σ::AbstractTrajectory) = σ.transitions
+
 # Get i-th state [i]
 getindex(σ::AbstractTrajectory, idx::Int) = get_state(σ, idx)
 # Get i-th value of var ["I", idx]
-function getindex(σ::AbstractTrajectory, var::String, idx::Int)
-    if var  == "times"
-        return σ.times[idx]
-    elseif var == "transitions"
-        return σ.transitions[idx]
-    else
-        return get_value(σ, var, idx)
-    end
-end
+getindex(σ::AbstractTrajectory, var::String, idx::Int) = get_value(σ, var, idx)
+# Get the path of a variable ["I"]
+getindex(σ::AbstractTrajectory, var::String) = get_var_values(σ, var)
 
diff --git a/models/ER.jl b/models/ER.jl
index 3c694b9345074a3a4510e6fc7ec38e4d12275631..243e0b8aa01b599851d47551462518c5ac6128b6 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -45,6 +45,6 @@ is_absorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) =
     (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
 g = ["P"]
 
-ER = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_f!,is_absorbing_ER; g=g)
+ER = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_f!,is_absorbing_ER; g=g)
 export ER
 
diff --git a/models/SIR.jl b/models/SIR.jl
index 215a4a87bc86347c46bd478bcb66ce4e7f9c4dd3..b3ab3aa31663a921638595a7af803779d82f2a81 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -34,7 +34,7 @@ function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String},
         b_sup += l_a[i+1]
     end
  
-    nu = @view l_nu[:,reaction] # macro for avoiding a copy
+    nu = view(l_nu, :, reaction) # macro for avoiding a copy
     for i = 1:3
         mat_x[idx,i] = xn[i]+nu[i]
     end
@@ -44,6 +44,6 @@ end
 is_absorbing_SIR(p::Vector{Float64}, xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_f!,is_absorbing_SIR; g=g)
+SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_f!,is_absorbing_SIR; g=g)
 export SIR
 
diff --git a/models/_bench_perf_test/ER_col.jl b/models/_bench_perf_test/ER_col.jl
index f237a78448398b3e6882acc26b6e428f62395a42..2546fd41d7434ea9ce918705d75c8c0f74969dab 100644
--- a/models/_bench_perf_test/ER_col.jl
+++ b/models/_bench_perf_test/ER_col.jl
@@ -46,6 +46,6 @@ is_absorbing_ER_col(p::Vector{Float64},xn::AbstractVector{Int}) =
     (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
 g = ["P"]
 
-ER_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,is_absorbing_ER_col; g=g)
+ER_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,is_absorbing_ER_col; g=g)
 export ER_col
 
diff --git a/models/_bench_perf_test/ER_col_buffer.jl b/models/_bench_perf_test/ER_col_buffer.jl
index c1c640aeba20ebd08534e45f0af5431a78c1b074..c7315c4e64b875868875bdb30273603fa6bc3e43 100644
--- a/models/_bench_perf_test/ER_col_buffer.jl
+++ b/models/_bench_perf_test/ER_col_buffer.jl
@@ -45,6 +45,6 @@ is_absorbing_ER_col_buffer(p::Vector{Float64},xn::AbstractVector{Int}) =
     (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
 g = ["P"]
 
-ER_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_buffer_f!,is_absorbing_ER_col_buffer; g=g)
+ER_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_buffer_f!,is_absorbing_ER_col_buffer; g=g)
 export ER_col_buffer
 
diff --git a/models/_bench_perf_test/ER_row_buffer.jl b/models/_bench_perf_test/ER_row_buffer.jl
index bc6104dd4c65cdbebfad1ff17cb09588941766ca..5c1489d1093cdd6f8dc4fd08bfb2ffc336f21c6f 100644
--- a/models/_bench_perf_test/ER_row_buffer.jl
+++ b/models/_bench_perf_test/ER_row_buffer.jl
@@ -45,6 +45,6 @@ is_absorbing_ER_row_buffer(p::Vector{Float64},xn::AbstractVector{Int}) =
     (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
 g = ["P"]
 
-ER_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_row_buffer_f!,is_absorbing_ER_row_buffer; g=g)
+ER_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_row_buffer_f!,is_absorbing_ER_row_buffer; g=g)
 export ER_row_buffer
 
diff --git a/models/_bench_perf_test/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl
index 09d2b2ac7758cf38894adccf1612db8d6974a5ae..3052b5e6b3a20238712d548d8b57c7b770fb35da 100644
--- a/models/_bench_perf_test/SIR_col.jl
+++ b/models/_bench_perf_test/SIR_col.jl
@@ -44,6 +44,6 @@ end
 is_absorbing_SIR_col(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,is_absorbing_SIR_col; g=g)
+SIR_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,is_absorbing_SIR_col; g=g)
 export SIR_col
 
diff --git a/models/_bench_perf_test/SIR_col_buffer.jl b/models/_bench_perf_test/SIR_col_buffer.jl
index 9f79b3e79eec3b7b40c498b20d786be57677b412..ee703931ed7a0b40708c2188b898ae9c3f431dba 100644
--- a/models/_bench_perf_test/SIR_col_buffer.jl
+++ b/models/_bench_perf_test/SIR_col_buffer.jl
@@ -44,6 +44,6 @@ end
 is_absorbing_SIR_col_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,is_absorbing_SIR_col_buffer; g=g)
+SIR_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,is_absorbing_SIR_col_buffer; g=g)
 export SIR_col_buffer
 
diff --git a/models/_bench_perf_test/SIR_row_buffer.jl b/models/_bench_perf_test/SIR_row_buffer.jl
index fcf6a86ea698314e924a4962af7a35b1ece809c6..4eaed57ebac6680c98ec4665abb7476990cf07ad 100644
--- a/models/_bench_perf_test/SIR_row_buffer.jl
+++ b/models/_bench_perf_test/SIR_row_buffer.jl
@@ -44,6 +44,6 @@ end
 is_absorbing_SIR_row_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g = ["I"]
 
-SIR_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,is_absorbing_SIR_row_buffer; g=g)
+SIR_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,is_absorbing_SIR_row_buffer; g=g)
 export SIR_row_buffer
 
diff --git a/tests/dist_lp/dist_case_2.jl b/tests/dist_lp/dist_case_2.jl
index ce29df428ad2ef34bff2392ccb96cf856fc5491e..be7c86d801b4bc071fe0d73187f29b2d191ba8b6 100644
--- a/tests/dist_lp/dist_case_2.jl
+++ b/tests/dist_lp/dist_case_2.jl
@@ -23,9 +23,9 @@ for p = 1:2
     l_tr = Vector{Nothing}(nothing, length(x_obs))
     σ1 = Trajectory(SIR, values, t_x, l_tr)
 
-    test_1 = isapprox(dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=p), res; atol = 1E-10) 
+    test_1 = isapprox(dist_lp(@view(x_obs[:]), t_x, @view(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[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); 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 = 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 c7d95586c3681f60edacea733f2bdafc44c63fe8..b27177121813e71a7453e8778311ef317c1d57f6 100644
--- a/tests/dist_lp/dist_case_3.jl
+++ b/tests/dist_lp/dist_case_3.jl
@@ -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[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=p)
+    res1 = dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p)
     res2 = dist_lp(σ1,σ2; p=p)
-    res1_bis = dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=p)
+    res1_bis = dist_lp(@view(y_obs[:]), t_y, @view(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 d02b7f3e343e071810a7021feb0b40ee6ab1f35f..a9152e0436815d868d1fcce0674f10647c64f3a6 100644
--- a/tests/dist_lp/dist_l1_case_1.jl
+++ b/tests/dist_lp/dist_l1_case_1.jl
@@ -21,7 +21,7 @@ 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[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=1) == 6.4
+test_1 = dist_lp(@view(x_obs[:]), t_x, @view(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[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=1) == 6.4 
+test_1_bis = dist_lp(@view(y_obs[:]), t_y, @view(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/dist_lp/dist_sim_sir.jl b/tests/dist_lp/dist_sim_sir.jl
index 5236040bbf550a07eb4a903e3c108534f82fd34e..13e3889bd426132f205c3cdf2f416ed6fdab8a0b 100644
--- a/tests/dist_lp/dist_sim_sir.jl
+++ b/tests/dist_lp/dist_sim_sir.jl
@@ -14,8 +14,8 @@ for p = 1:2
         d = dist_lp(σ1, σ2, "I"; p = p)
         d2 = dist_lp(σ1, σ2; p = p)
 
-        f_x(t::Float64) = MarkovProcesses._f_step(σ1["I"], σ1["times"], t)
-        f_y(t::Float64) = MarkovProcesses._f_step(σ2["I"], σ2["times"], t)
+        f_x(t::Float64) = MarkovProcesses._f_step(σ1["I"], times(σ1), t)
+        f_y(t::Float64) = MarkovProcesses._f_step(σ2["I"], times(σ2), t)
         diff_f(t) = abs(f_x(t) - f_y(t))^p
         int_riemann = MarkovProcesses._riemann_sum(diff_f, SIR.t0, SIR.time_bound, 1E-3)
         res_int_riemann = int_riemann^(1/p)
diff --git a/tests/simulation/sim_er.jl b/tests/simulation/sim_er.jl
index 09d4575b3867d5c5121448031d15717e3a14ed05..21e18e22659b79879c4345c516a02a1776c4cd1a 100644
--- a/tests/simulation/sim_er.jl
+++ b/tests/simulation/sim_er.jl
@@ -6,7 +6,7 @@ load_model("ER")
 
 σ = simulate(ER)
 plt.figure()
-plt.step(σ["times"], σ["P"], "ro--", marker="x", where="post", linewidth=1.0)
+plt.step(times(σ), σ["P"], "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er.png", dpi=480)
 plt.close()
 
diff --git a/tests/simulation/sim_er_row_buffer_bounded.jl b/tests/simulation/sim_er_row_buffer_bounded.jl
index 6da83769f8beed6be145263dcb1354556451614f..d684462f6f2e6246dffaf69edbfae167324f4b91 100644
--- a/tests/simulation/sim_er_row_buffer_bounded.jl
+++ b/tests/simulation/sim_er_row_buffer_bounded.jl
@@ -8,7 +8,7 @@ ER_row_buffer.time_bound = 10.0
 
 σ = _simulate_row_buffer(ER_row_buffer)
 plt.figure()
-plt.step(σ["times"], _get_values_row(σ,"P"), "ro--", marker="x", where="post", linewidth=1.0)
+plt.step(times(σ), _get_values_row(σ,"P"), "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er_row_buffer_bounded.png")
 plt.close()
 
diff --git a/tests/simulation/sim_sir.jl b/tests/simulation/sim_sir.jl
index 2a1f00fa21f22f798834fbfd7b4f8edc98ca358b..74fd911d82c083e54d1d85450cd11da3a0f1ca6a 100644
--- a/tests/simulation/sim_sir.jl
+++ b/tests/simulation/sim_sir.jl
@@ -6,7 +6,7 @@ load_model("SIR")
 
 σ = simulate(SIR)
 plt.figure()
-plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
+plt.step(times(σ), σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir.png")
 plt.close()
 
diff --git a/tests/simulation/sim_sir_bounded.jl b/tests/simulation/sim_sir_bounded.jl
index 95d64343b1769a9d6bf7d2b9cde10d7f913dfacc..469e7844bc6bd1ac8969b0c694ecb76b052c9ba4 100644
--- a/tests/simulation/sim_sir_bounded.jl
+++ b/tests/simulation/sim_sir_bounded.jl
@@ -7,7 +7,7 @@ SIR.time_bound = 100.0
 
 σ = simulate(SIR)
 plt.figure()
-plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
+plt.step(times(σ), σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_bounded.png")
 plt.close()
 
diff --git a/tests/simulation/sim_sir_col_buffer_bounded.jl b/tests/simulation/sim_sir_col_buffer_bounded.jl
index 334349553d9cd9720522b7f003ad7114750195b8..524487b03ba84e6756e7a7379d938b56e55424d5 100644
--- a/tests/simulation/sim_sir_col_buffer_bounded.jl
+++ b/tests/simulation/sim_sir_col_buffer_bounded.jl
@@ -8,7 +8,7 @@ SIR_col_buffer.time_bound = 100.0
 
 σ = _simulate_col_buffer(SIR_col_buffer)
 plt.figure()
-plt.step(σ["times"], _get_values_col(σ,"I"), "ro--", marker="x", where="post", linewidth=1.0)
+plt.step(times(σ), _get_values_col(σ,"I"), "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_col_buffer_bounded.png")
 plt.close()
 
diff --git a/tests/simulation/sim_sir_row_buffer_bounded.jl b/tests/simulation/sim_sir_row_buffer_bounded.jl
index b6a38981e1d574ffa49a6bce8f28b986a5c86c0f..b49c589841824c71223ce8f1692fd81699db4468 100644
--- a/tests/simulation/sim_sir_row_buffer_bounded.jl
+++ b/tests/simulation/sim_sir_row_buffer_bounded.jl
@@ -8,7 +8,7 @@ SIR_row_buffer.time_bound = 100.0
 
 σ = _simulate_row_buffer(SIR_row_buffer)
 plt.figure()
-plt.step(σ["times"], _get_values_row(σ,"I"), "ro--", marker="x", where="post", linewidth=1.0)
+plt.step(times(σ), _get_values_row(σ,"I"), "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_row_buffer_bounded.png")
 plt.close()
 
diff --git a/tests/unit/simulate_sir_bounded.jl b/tests/unit/simulate_sir_bounded.jl
index 0a7a10f02b475d8c365bd82d37316d60e94cd3e4..c872c02612cae7273e842c3592e73ca069efe316 100644
--- a/tests/unit/simulate_sir_bounded.jl
+++ b/tests/unit/simulate_sir_bounded.jl
@@ -11,7 +11,7 @@ d2 = Dict("I" => 1)
 
 bool_test = SIR.g == ["I"] && SIR._g_idx == [2] && 
             SIR.map_var_idx == d1 && 
-            SIR._map_obs_var_idx == d2 && σ["times"][end] <= SIR.time_bound
+            SIR._map_obs_var_idx == d2 && times(σ)[end] <= SIR.time_bound
 
 return bool_test