diff --git a/bench/array_memory_order/access_trajectory_sir.jl b/bench/array_memory_order/access_trajectory_sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6189405f69e86a353a9509d7025033f628545a41
--- /dev/null
+++ b/bench/array_memory_order/access_trajectory_sir.jl
@@ -0,0 +1,49 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+bound_time = 200.0
+l_var = ["S", "I", "R"]
+
+println("Col buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+function access_trajectory_col(m::Model)
+    res = 0
+    n_sim = 100
+    for k = 1:n_sim
+        σ = _simulate_col_buffer(m)
+        t = _get_values_col(σ, "I")
+        res += t[end-1]
+    end
+    return res
+end
+# Bench
+@timev access_trajectory_col(SIR_col_buffer) 
+b1_col = @benchmark access_trajectory_col($SIR_col_buffer)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Row buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+function access_trajectory_row(m::Model)
+    res = 0
+    n_sim = 100
+    for k = 1:n_sim
+        σ = _simulate_row_buffer(m)
+        t = _get_values_row(σ, "I")
+        res += t[end-1]
+    end
+    return res
+end
+# Bench
+@timev access_trajectory_row(SIR_row_buffer) 
+b1_row = @benchmark access_trajectory_row($SIR_row_buffer)
+@show minimum(b1_row), mean(b1_row), maximum(b1_row)
+
diff --git a/bench/array_memory_order/multiple_sim_sir.jl b/bench/array_memory_order/multiple_sim_sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..2a9efbdbe07d7eb1a436a187d27f34f314b700e1
--- /dev/null
+++ b/bench/array_memory_order/multiple_sim_sir.jl
@@ -0,0 +1,48 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+l_var = ["S", "I", "R"]
+bound_time = 200.0
+nbr_sim = 10000
+
+load_model("_bench_perf_test/SIR_col")
+SIR_col.time_bound = bound_time
+set_observed_var!(SIR_col, l_var)
+
+println("Col")
+b1_col = @benchmark for i = 1:$(nbr_sim) _simulate_col($SIR_col) end
+@timev _simulate_col(SIR_col)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+
+println("Col + buffer:")
+b1_col_buffer = @benchmark for i = 1:$(nbr_sim) _simulate_col_buffer($SIR_col_buffer) end
+@timev _simulate_col_buffer(SIR_col_buffer)
+@show minimum(b1_col_buffer), mean(b1_col_buffer), maximum(b1_col_buffer)
+
+println("Col + buffer_10:")
+b1_col_buffer_10 = @benchmark for i = 1:$(nbr_sim) _simulate_col_buffer($SIR_col_buffer; buffer_size = 10) end
+@timev _simulate_col_buffer(SIR_col_buffer; buffer_size = 10)
+@show minimum(b1_col_buffer_10), mean(b1_col_buffer_10), maximum(b1_col_buffer_10)
+
+load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+
+println("Row + buffer:")
+b1_row_buffer = @benchmark for i = 1:$(nbr_sim) _simulate_row_buffer($SIR_row_buffer) end
+@timev _simulate_row_buffer(SIR_row_buffer)
+@show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
+
+println("Row + buffer_10:")
+b1_row_buffer_10 = @benchmark for i = 1:$(nbr_sim) _simulate_row_buffer($SIR_row_buffer; buffer_size = 10) end
+@timev _simulate_row_buffer(SIR_row_buffer; buffer_size = 10)
+@show minimum(b1_row_buffer_10), mean(b1_row_buffer_10), maximum(b1_row_buffer_10)
+
diff --git a/bench/array_memory_order/read_random_state_trajectory.jl b/bench/array_memory_order/read_random_state_trajectory.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c8efbcfd7fed8687aaa8c0e7d62a584d851c9dca
--- /dev/null
+++ b/bench/array_memory_order/read_random_state_trajectory.jl
@@ -0,0 +1,51 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+l_var = ["S", "I", "R"]
+bound_time = 200.0
+
+println("Col buffer:")
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+function random_trajectory_value_col(m::ContinuousTimeModel)
+    σ = _simulate_col_buffer(m)
+    n_states = get_states_number(σ)
+    nb_rand = 1000
+    res = 0
+    for i = 1:nb_rand
+        a = _get_state_col(σ, rand(1:n_states))
+        res += a[2]
+    end
+    return res
+end
+# Bench
+@timev random_trajectory_value_col(SIR_col_buffer) 
+b1_col_buffer = @benchmark random_trajectory_value_col($SIR_col_buffer)
+@show minimum(b1_col_buffer), mean(b1_col_buffer), maximum(b1_col_buffer)
+
+println("Row buffer:")
+
+load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+function random_trajectory_value_row(m::ContinuousTimeModel)
+    σ = _simulate_row_buffer(m)
+    n_states = get_states_number(σ)
+    nb_rand = 1000
+    res = 0
+    for i = 1:nb_rand
+        a = _get_state_row(σ, rand(1:n_states))
+        res += a[2]
+    end
+    return res
+end
+# Bench
+@timev random_trajectory_value_row(SIR_row_buffer) 
+b1_row_buffer = @benchmark random_trajectory_value_row($SIR_row_buffer)
+@show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
+
diff --git a/bench/array_memory_order/read_trajectory_sir.jl b/bench/array_memory_order/read_trajectory_sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..2a32ccfba46c3fcbfbe7af05d6f2a626fd6c1568
--- /dev/null
+++ b/bench/array_memory_order/read_trajectory_sir.jl
@@ -0,0 +1,53 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+bound_time = 200.0
+l_var = ["S", "I", "R"]
+
+println("Col buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+function read_trajectory_col(m::Model)
+    res = 0
+    σ = _simulate_col_buffer(m)
+    n_states = get_states_number(σ)
+    n_read = 100000
+    for k = 1:n_read
+        for i = 1:n_states
+            res += _get_value_col(σ, "I", i)
+        end
+    end
+    return res
+end
+# Bench
+@timev read_trajectory_col(SIR_col_buffer) 
+b1_col = @benchmark read_trajectory_col($SIR_col_buffer)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Row buffer:")
+
+MarkovProcesses.load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+function read_trajectory_row(m::Model)
+    res = 0
+    σ = _simulate_row_buffer(m)
+    n_states = get_states_number(σ)
+    n_read = 100000
+    for k = 1:n_read
+        for i = 1:n_states
+            res += _get_value_row(σ, "I", i)
+        end
+    end
+    return res
+end
+# Bench
+@timev read_trajectory_row(SIR_row_buffer) 
+b1_row = @benchmark read_trajectory_row($SIR_row_buffer)
+@show minimum(b1_row), mean(b1_row), maximum(b1_row)
+
diff --git a/bench/array_memory_order/sim_sir.jl b/bench/array_memory_order/sim_sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..84a5d14803ddd57355bcd01556a77db08b52258a
--- /dev/null
+++ b/bench/array_memory_order/sim_sir.jl
@@ -0,0 +1,48 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+l_var = ["S", "I", "R"]
+bound_time = 200.0
+
+load_model("_bench_perf_test/SIR_col")
+SIR_col.time_bound = bound_time
+set_observed_var!(SIR_col, l_var)
+
+println("Col")
+b1_col = @benchmark _simulate_col($SIR_col)
+@timev _simulate_col(SIR_col)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = bound_time
+set_observed_var!(SIR_col_buffer, l_var)
+
+println("Col + buffer:")
+b1_col_buffer = @benchmark _simulate_col_buffer($SIR_col_buffer)
+@timev _simulate_col_buffer(SIR_col_buffer)
+@show minimum(b1_col_buffer), mean(b1_col_buffer), maximum(b1_col_buffer)
+
+println("Col + buffer_10:")
+b1_col_buffer_10 = @benchmark _simulate_col_buffer($SIR_col_buffer; buffer_size = 10)
+@timev _simulate_col_buffer(SIR_col_buffer; buffer_size = 10)
+@show minimum(b1_col_buffer_10), mean(b1_col_buffer_10), maximum(b1_col_buffer_10)
+
+load_model("_bench_perf_test/SIR_row_buffer")
+SIR_row_buffer.time_bound = bound_time
+set_observed_var!(SIR_row_buffer, l_var)
+
+println("Row + buffer:")
+b1_row_buffer = @benchmark _simulate_row_buffer($SIR_row_buffer)
+@timev _simulate_row_buffer(SIR_row_buffer)
+@show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
+
+println("Row + buffer_10:")
+b1_row_buffer_10 = @benchmark _simulate_row_buffer($SIR_row_buffer; buffer_size = 10)
+@timev _simulate_row_buffer(SIR_row_buffer; buffer_size = 10)
+@show minimum(b1_row_buffer_10), mean(b1_row_buffer_10), maximum(b1_row_buffer_10)
+
+
diff --git a/bench/pygmalion/read_random_state_trajectory.jl b/bench/pygmalion/read_random_state_trajectory.jl
index 05c444c5d43a69ceec90a45d03b3898824be0c10..9877751d9c8468cc1c187582d9ba71aab1fb2e79 100644
--- a/bench/pygmalion/read_random_state_trajectory.jl
+++ b/bench/pygmalion/read_random_state_trajectory.jl
@@ -8,7 +8,8 @@ println("Pygmalion:")
 str_m = "enzymatic_reaction"
 str_d = "abc_er"
 pygmalion.load_model(str_m)
-str_oml = "P,R,time"
+l_var = ["E", "S", "ES", "P"]
+str_oml = "E,S,ES,P,R,time"
 ll_om = split(str_oml, ",")
 x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
 p_true = Parameters(1.0, 1.0, 1.0)
@@ -20,6 +21,10 @@ function random_trajectory_value_pyg(so::SystemObservation)
     n_states = get_number_of_observations(so, "P")
     return to_vec(so, "P", rand(1:n_states))
 end
+function random_trajectory_state_pyg(so::SystemObservation)
+    n_states = get_number_of_observations(so, "P")
+    return to_vec(so, so.oml, rand(1:n_states))
+end
 # Bench
 @timev random_trajectory_value_pyg(so)
 b1_pyg = @benchmark random_trajectory_value_pyg($so)
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index 06e01d179268aead933b13e5fc3db11d1fb52a97..1562f41c8d348a1746c1b86b081db888c6f556e3 100644
--- a/core/_tests_simulate.jl
+++ b/core/_tests_simulate.jl
@@ -4,14 +4,14 @@
 # Trajectories
 
 _get_values_col(σ::AbstractTrajectory, var::String) = 
-σ.values[(σ.m)._map_obs_var_idx[var],:] 
+@view σ.values[(σ.m)._map_obs_var_idx[var],:] 
 _get_values_row(σ::AbstractTrajectory, var::String) = 
-σ.values[:,(σ.m)._map_obs_var_idx[var]] 
+@view σ.values[:,(σ.m)._map_obs_var_idx[var]] 
 
 _get_state_col(σ::AbstractTrajectory, idx::Int) = 
-σ.values[:,idx]
+@view σ.values[:,idx]
 _get_state_row(σ::AbstractTrajectory, idx::Int) = 
-σ.values[idx,:]
+@view σ.values[idx,:]
 
 _get_value_col(σ::AbstractTrajectory, var::String, idx::Int) = 
 σ.values[(σ.m)._map_obs_var_idx[var],idx] 
diff --git a/core/model.jl b/core/model.jl
index d37ad6d6fe8425399956b223ff4ed95b6af88e2e..a1b584c8ebe28da198a624cf7a27aae8c1cb0d39 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -32,6 +32,14 @@ function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_tra
         _g_idx[i] = map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::String => idx in state space )
         _map_obs_var_idx[g[i]] = i
     end
+  
+    if length(methods(f!)) >= 2
+        @warn "You have possibly redefined a function Model.f! used in a previously instantiated model."
+    end
+    if length(methods(is_absorbing)) >= 2
+        @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)
 end
 
diff --git a/core/observations.jl b/core/observations.jl
index 5967ea533a2601d3c6690cbb3e769303eb8f057a..e82b222085b280a0ae6f910c30424f8a91caff8d 100644
--- a/core/observations.jl
+++ b/core/observations.jl
@@ -4,9 +4,9 @@ ContinuousObservations = AbstractVector{AbstractTrajectory}
 
 struct Trajectory <: AbstractTrajectory
     m::ContinuousTimeModel
-    values::AbstractMatrix{Float64}
-    times::AbstractVector{Float64}
-    transitions::AbstractVector{Union{String,Nothing}}
+    values::Matrix{Int}
+    times::Vector{Float64}
+    transitions::Vector{Union{String,Nothing}}
 end
 
 function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
diff --git a/models/SIR.jl b/models/SIR.jl
index 7416631b94d471155178b23f2195376b4173ba29..d448cc2ac141b678786bc7b162a0622b74dbcbcb 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -9,7 +9,7 @@ l_tr = ["R1","R2"]
 p = [0.0012, 0.05]
 x0 = [95, 5, 0]
 t0 = 0.0
-function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
+function SIR_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
            xn::Vector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
@@ -41,9 +41,9 @@ function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
     tnplus1[1] = tn + tau
     tr[1] = "R$(reaction)"
 end
-is_absorbing_sir(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+is_absorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (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,f!,is_absorbing_sir; g=g)
+SIR = CTMC(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/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl
index 00cf2e48e53862d49bb00f4310e1cffe142151dc..919e8fc73fb40ce38b58ec33a8547c892ed8fc4a 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 f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
+function SIR_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
            xn::Vector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
@@ -41,9 +41,9 @@ function f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String},
     tnplus1[1] = tn + tau
     tr[1] = "R$(reaction)"
 end
-is_absorbing_sir(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+is_absorbing_SIR_col(p::Vector{Float64}, xn::Vector{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,f!,is_absorbing_sir; g=g)
+SIR_col = CTMC(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 d37828bdbb2eee8145513b1d2217afe11749294d..9f79b3e79eec3b7b40c498b20d786be57677b412 100644
--- a/models/_bench_perf_test/SIR_col_buffer.jl
+++ b/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 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{String}, idx::Int,
            xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
@@ -41,9 +41,9 @@ function f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx:
     l_t[idx] = tn + tau
     l_tr[idx] = "R$(reaction)"
 end
-is_absorbing_sir(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+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,f!,is_absorbing_sir; g=g)
+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)
 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 366b68b9d120fc8368d32bacdd810af87b8ea3d8..fcf6a86ea698314e924a4962af7a35b1ece809c6 100644
--- a/models/_bench_perf_test/SIR_row_buffer.jl
+++ b/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 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{String}, idx::Int,
            xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
@@ -41,9 +41,9 @@ function f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, idx:
     l_t[idx] = tn + tau
     l_tr[idx] = "R$(reaction)"
 end
-is_absorbing_sir(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+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,f!,is_absorbing_sir; g=g)
+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)
 export SIR_row_buffer