diff --git a/bench/array_memory_order/Notes b/bench/array_memory_order/Notes
new file mode 100644
index 0000000000000000000000000000000000000000..e6916bd05e24d41e0aa82f41c69c47d7ba7ce5cb
--- /dev/null
+++ b/bench/array_memory_order/Notes
@@ -0,0 +1,2 @@
+julia script.jl SIR|ER to run benchmarks
+
diff --git a/bench/array_memory_order/access_trajectory.jl b/bench/array_memory_order/access_trajectory.jl
new file mode 100644
index 0000000000000000000000000000000000000000..39be2cd79ebcba33b979bb4d67f34966f227a63f
--- /dev/null
+++ b/bench/array_memory_order/access_trajectory.jl
@@ -0,0 +1,71 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+if ARGS[1] == "SIR"
+    bound_time = 200.0
+    l_var = ["S", "I", "R"]
+
+    load_model("_bench_perf_test/SIR_col_buffer")
+    SIR_col_buffer.time_bound = bound_time
+    set_observed_var!(SIR_col_buffer, l_var)
+    load_model("_bench_perf_test/SIR_row_buffer")
+    SIR_row_buffer.time_bound = bound_time
+    set_observed_var!(SIR_row_buffer, l_var)
+
+    model_col_buffer = SIR_col_buffer
+    model_row_buffer = SIR_row_buffer
+elseif ARGS[1] == "ER"
+    l_var = ["E","S","ES","P"]
+    bound_time = 20.0
+
+    load_model("_bench_perf_test/ER_col_buffer")
+    ER_col_buffer.time_bound = bound_time
+    set_observed_var!(ER_col_buffer, l_var)
+    load_model("_bench_perf_test/ER_row_buffer")
+    ER_row_buffer.time_bound = bound_time
+    set_observed_var!(ER_row_buffer, l_var)
+
+    model_col_buffer = ER_col_buffer
+    model_row_buffer = ER_row_buffer
+else
+    error("Unavailable model")
+end
+
+println("Col buffer:")
+
+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(model_col_buffer) 
+b1_col = @benchmark access_trajectory_col($model_col_buffer)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Row buffer:")
+
+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(model_row_buffer) 
+b1_row = @benchmark access_trajectory_row($model_row_buffer)
+@show minimum(b1_row), mean(b1_row), maximum(b1_row)
+
diff --git a/bench/array_memory_order/access_trajectory_sir.jl b/bench/array_memory_order/access_trajectory_sir.jl
deleted file mode 100644
index 6189405f69e86a353a9509d7025033f628545a41..0000000000000000000000000000000000000000
--- a/bench/array_memory_order/access_trajectory_sir.jl
+++ /dev/null
@@ -1,49 +0,0 @@
-
-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.jl b/bench/array_memory_order/multiple_sim.jl
new file mode 100644
index 0000000000000000000000000000000000000000..1fb5fc6e5aed941800c9f713afb2f7abe0316e2d
--- /dev/null
+++ b/bench/array_memory_order/multiple_sim.jl
@@ -0,0 +1,78 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+if ARGS[1] == "SIR"
+    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)
+    load_model("_bench_perf_test/SIR_col_buffer")
+    SIR_col_buffer.time_bound = bound_time
+    set_observed_var!(SIR_col_buffer, l_var)
+    load_model("_bench_perf_test/SIR_row_buffer")
+    SIR_row_buffer.time_bound = bound_time
+    set_observed_var!(SIR_row_buffer, l_var)
+
+    model_col = SIR_col
+    model_col_buffer = SIR_col_buffer
+    model_row_buffer = SIR_row_buffer
+elseif ARGS[1] == "ER"
+    l_var = ["E","S","ES","P"]
+    bound_time = 20.0
+    nbr_sim = 10000
+
+    load_model("_bench_perf_test/ER_col")
+    ER_col.time_bound = bound_time
+    set_observed_var!(ER_col, l_var)
+    load_model("_bench_perf_test/ER_col_buffer")
+    ER_col_buffer.time_bound = bound_time
+    set_observed_var!(ER_col_buffer, l_var)
+    load_model("_bench_perf_test/ER_row_buffer")
+    ER_row_buffer.time_bound = bound_time
+    set_observed_var!(ER_row_buffer, l_var)
+
+    model_col = ER_col
+    model_col_buffer = ER_col_buffer
+    model_row_buffer = ER_row_buffer
+else
+    error("Unavailable model")
+end
+nbr_sim = 10000
+
+println("Col")
+b1_col = @benchmark for i = 1:$(nbr_sim) _simulate_col($model_col) end
+@timev _simulate_col(model_col)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Col + buffer:")
+b1_col_buffer = @benchmark for i = 1:$(nbr_sim) _simulate_col_buffer($model_col_buffer) end
+@timev _simulate_col_buffer(model_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($model_col_buffer; buffer_size = 10) end
+@timev _simulate_col_buffer(model_col_buffer; buffer_size = 10)
+@show minimum(b1_col_buffer_10), mean(b1_col_buffer_10), maximum(b1_col_buffer_10)
+
+
+println("Row + buffer:")
+b1_row_buffer = @benchmark for i = 1:$(nbr_sim) _simulate_row_buffer($model_row_buffer) end
+@timev _simulate_row_buffer(model_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($model_row_buffer; buffer_size = 10) end
+@timev _simulate_row_buffer(model_row_buffer; buffer_size = 10)
+@show minimum(b1_row_buffer_10), mean(b1_row_buffer_10), maximum(b1_row_buffer_10)
+
+println("Row + buffer_100:")
+b1_row_buffer_100 = @benchmark for i = 1:$(nbr_sim) _simulate_row_buffer($model_row_buffer; buffer_size = 100) end
+@timev _simulate_row_buffer(model_row_buffer; buffer_size = 100)
+@show minimum(b1_row_buffer_100), mean(b1_row_buffer_100), maximum(b1_row_buffer_100)
+
diff --git a/bench/array_memory_order/multiple_sim_sir.jl b/bench/array_memory_order/multiple_sim_sir.jl
deleted file mode 100644
index 2a9efbdbe07d7eb1a436a187d27f34f314b700e1..0000000000000000000000000000000000000000
--- a/bench/array_memory_order/multiple_sim_sir.jl
+++ /dev/null
@@ -1,48 +0,0 @@
-
-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
index c8efbcfd7fed8687aaa8c0e7d62a584d851c9dca..c7b869c32075721e93767539e20dafd1524312a2 100644
--- a/bench/array_memory_order/read_random_state_trajectory.jl
+++ b/bench/array_memory_order/read_random_state_trajectory.jl
@@ -4,14 +4,39 @@ import BenchmarkTools: mean
 using MarkovProcesses
 include(get_module_path() * "/core/_tests_simulate.jl")
 
-l_var = ["S", "I", "R"]
-bound_time = 200.0
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+if ARGS[1] == "SIR"
+    bound_time = 200.0
+    l_var = ["S", "I", "R"]
+
+    load_model("_bench_perf_test/SIR_col_buffer")
+    SIR_col_buffer.time_bound = bound_time
+    set_observed_var!(SIR_col_buffer, l_var)
+    load_model("_bench_perf_test/SIR_row_buffer")
+    SIR_row_buffer.time_bound = bound_time
+    set_observed_var!(SIR_row_buffer, l_var)
+
+    model_col_buffer = SIR_col_buffer
+    model_row_buffer = SIR_row_buffer
+elseif ARGS[1] == "ER"
+    l_var = ["E","S","ES","P"]
+    bound_time = 20.0
+
+    load_model("_bench_perf_test/ER_col_buffer")
+    ER_col_buffer.time_bound = bound_time
+    set_observed_var!(ER_col_buffer, l_var)
+    load_model("_bench_perf_test/ER_row_buffer")
+    ER_row_buffer.time_bound = bound_time
+    set_observed_var!(ER_row_buffer, l_var)
+
+    model_col_buffer = ER_col_buffer
+    model_row_buffer = ER_row_buffer
+else
+    error("Unavailable model")
+end
 
 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(σ)
@@ -24,15 +49,12 @@ function random_trajectory_value_col(m::ContinuousTimeModel)
     return res
 end
 # Bench
-@timev random_trajectory_value_col(SIR_col_buffer) 
-b1_col_buffer = @benchmark random_trajectory_value_col($SIR_col_buffer)
+@timev random_trajectory_value_col(model_col_buffer) 
+b1_col_buffer = @benchmark random_trajectory_value_col($model_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(σ)
@@ -45,7 +67,7 @@ function random_trajectory_value_row(m::ContinuousTimeModel)
     return res
 end
 # Bench
-@timev random_trajectory_value_row(SIR_row_buffer) 
-b1_row_buffer = @benchmark random_trajectory_value_row($SIR_row_buffer)
+@timev random_trajectory_value_row(model_row_buffer) 
+b1_row_buffer = @benchmark random_trajectory_value_row($model_row_buffer)
 @show minimum(b1_row_buffer), mean(b1_row_buffer), maximum(b1_row_buffer)
 
diff --git a/bench/array_memory_order/read_trajectory.jl b/bench/array_memory_order/read_trajectory.jl
new file mode 100644
index 0000000000000000000000000000000000000000..65d9e4738f592f9dd83beaaeb5f4dfd7454cc87f
--- /dev/null
+++ b/bench/array_memory_order/read_trajectory.jl
@@ -0,0 +1,75 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+if ARGS[1] == "SIR"
+    bound_time = 200.0
+    l_var = ["S", "I", "R"]
+
+    load_model("_bench_perf_test/SIR_col_buffer")
+    SIR_col_buffer.time_bound = bound_time
+    set_observed_var!(SIR_col_buffer, l_var)
+    load_model("_bench_perf_test/SIR_row_buffer")
+    SIR_row_buffer.time_bound = bound_time
+    set_observed_var!(SIR_row_buffer, l_var)
+
+    model_col_buffer = SIR_col_buffer
+    model_row_buffer = SIR_row_buffer
+elseif ARGS[1] == "ER"
+    l_var = ["E","S","ES","P"]
+    bound_time = 20.0
+
+    load_model("_bench_perf_test/ER_col_buffer")
+    ER_col_buffer.time_bound = bound_time
+    set_observed_var!(ER_col_buffer, l_var)
+    load_model("_bench_perf_test/ER_row_buffer")
+    ER_row_buffer.time_bound = bound_time
+    set_observed_var!(ER_row_buffer, l_var)
+
+    model_col_buffer = ER_col_buffer
+    model_row_buffer = ER_row_buffer
+else
+    error("Unavailable model")
+end
+
+println("Col buffer:")
+
+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(model_col_buffer) 
+b1_col = @benchmark read_trajectory_col($model_col_buffer)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Row buffer:")
+
+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(model_row_buffer) 
+b1_row = @benchmark read_trajectory_row($model_row_buffer)
+@show minimum(b1_row), mean(b1_row), maximum(b1_row)
+
diff --git a/bench/array_memory_order/read_trajectory_sir.jl b/bench/array_memory_order/read_trajectory_sir.jl
deleted file mode 100644
index 2a32ccfba46c3fcbfbe7af05d6f2a626fd6c1568..0000000000000000000000000000000000000000
--- a/bench/array_memory_order/read_trajectory_sir.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-
-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.jl b/bench/array_memory_order/sim.jl
new file mode 100644
index 0000000000000000000000000000000000000000..79fbc4fc8b377b17fbc8bcb5837d751e41be9ed2
--- /dev/null
+++ b/bench/array_memory_order/sim.jl
@@ -0,0 +1,73 @@
+
+using BenchmarkTools
+import BenchmarkTools: mean
+BenchmarkTools.DEFAULT_PARAMETERS.samples = 20000
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+
+if ARGS[1] == "SIR"
+    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)
+    load_model("_bench_perf_test/SIR_col_buffer")
+    SIR_col_buffer.time_bound = bound_time
+    set_observed_var!(SIR_col_buffer, l_var)
+    load_model("_bench_perf_test/SIR_row_buffer")
+    SIR_row_buffer.time_bound = bound_time
+    set_observed_var!(SIR_row_buffer, l_var)
+
+    model_col = SIR_col
+    model_col_buffer = SIR_col_buffer
+    model_row_buffer = SIR_row_buffer
+elseif ARGS[1] == "ER"
+    l_var = ["E","S","ES","P"]
+    bound_time = 20.0
+    nbr_sim = 10000
+
+    load_model("_bench_perf_test/ER_col")
+    ER_col.time_bound = bound_time
+    set_observed_var!(ER_col, l_var)
+    load_model("_bench_perf_test/ER_col_buffer")
+    ER_col_buffer.time_bound = bound_time
+    set_observed_var!(ER_col_buffer, l_var)
+    load_model("_bench_perf_test/ER_row_buffer")
+    ER_row_buffer.time_bound = bound_time
+    set_observed_var!(ER_row_buffer, l_var)
+
+    model_col = ER_col
+    model_col_buffer = ER_col_buffer
+    model_row_buffer = ER_row_buffer
+else
+    error("Unavailable model")
+end
+nbr_sim = 10000
+
+println("Col")
+b1_col = @benchmark _simulate_col($model_col)
+@timev _simulate_col(model_col)
+@show minimum(b1_col), mean(b1_col), maximum(b1_col)
+
+println("Col + buffer:")
+b1_col_buffer = @benchmark _simulate_col_buffer($model_col_buffer)
+@timev _simulate_col_buffer(model_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($model_col_buffer; buffer_size = 10)
+@timev _simulate_col_buffer(model_col_buffer; buffer_size = 10)
+@show minimum(b1_col_buffer_10), mean(b1_col_buffer_10), maximum(b1_col_buffer_10)
+
+println("Row + buffer:")
+b1_row_buffer = @benchmark _simulate_row_buffer($model_row_buffer)
+@timev _simulate_row_buffer(model_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($model_row_buffer; buffer_size = 10)
+@timev _simulate_row_buffer(model_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/sim_sir.jl b/bench/array_memory_order/sim_sir.jl
deleted file mode 100644
index 84a5d14803ddd57355bcd01556a77db08b52258a..0000000000000000000000000000000000000000
--- a/bench/array_memory_order/sim_sir.jl
+++ /dev/null
@@ -1,48 +0,0 @@
-
-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/models/_bench_perf_test/ER_col.jl b/models/_bench_perf_test/ER_col.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f237a78448398b3e6882acc26b6e428f62395a42
--- /dev/null
+++ b/models/_bench_perf_test/ER_col.jl
@@ -0,0 +1,51 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+d=4
+k=3
+dict_var = Dict("E" => 1, "S" => 2, "ES" => 3, "P" => 4)
+dict_p = Dict("k1" => 1, "k2" => 2, "k3" => 3)
+l_tr = ["R1","R2","R3"]
+p = [1.0, 1.0, 1.0]
+x0 = [100, 100, 0, 0]
+t0 = 0.0
+
+function ER_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{String}, 
+           xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
+    a1 = p[1] * xn[1] * xn[2]
+    a2 = p[2] * xn[3]
+    a3 = p[3] * xn[3]
+    l_a = SVector(a1, a2, a3)
+    asum = sum(l_a)
+    l_nu = @SMatrix [-1 1 1;
+                     -1 1 0;
+                     1 -1 -1;
+                     0 0 1]
+    u1, u2 = rand(), rand()
+    tau = - log(u1) / asum
+    b_inf = 0.0
+    b_sup = a1
+    reaction = 0
+    for i = 1:3
+        if b_inf < asum*u2 < b_sup
+            reaction = i
+            break
+        end
+        b_inf += l_a[i]
+        b_sup += l_a[i+1]
+    end
+ 
+    nu = @view l_nu[:,reaction]
+    for i = 1:4
+        xnplus1[i] = xn[i]+nu[i]
+    end
+    tnplus1[1] = tn + tau
+    tr[1] = "R$(reaction)"
+end
+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)
+export ER_col
+
diff --git a/models/_bench_perf_test/ER_col_buffer.jl b/models/_bench_perf_test/ER_col_buffer.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c1c640aeba20ebd08534e45f0af5431a78c1b074
--- /dev/null
+++ b/models/_bench_perf_test/ER_col_buffer.jl
@@ -0,0 +1,50 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+d=4
+k=3
+dict_var = Dict("E" => 1, "S" => 2, "ES" => 3, "P" => 4)
+dict_p = Dict("k1" => 1, "k2" => 2, "k3" => 3)
+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,
+           xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
+    a1 = p[1] * xn[1] * xn[2]
+    a2 = p[2] * xn[3]
+    a3 = p[3] * xn[3]
+    l_a = SVector(a1, a2, a3)
+    asum = sum(l_a)
+    l_nu = @SMatrix [-1 1 1;
+                     -1 1 0;
+                     1 -1 -1;
+                     0 0 1]
+    u1, u2 = rand(), rand()
+    tau = - log(u1) / asum
+    b_inf = 0.0
+    b_sup = a1
+    reaction = 0
+    for i = 1:3
+        if b_inf < asum*u2 < b_sup
+            reaction = i
+            break
+        end
+        b_inf += l_a[i]
+        b_sup += l_a[i+1]
+    end
+ 
+    nu = @view l_nu[:,reaction] # macro for avoiding a copy
+    for i = 1:4
+        mat_x[i,idx] = xn[i]+nu[i]
+    end
+    l_t[idx] = tn + tau
+    l_tr[idx] = "R$(reaction)"
+end
+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)
+export ER_col_buffer
+
diff --git a/models/_bench_perf_test/ER_row_buffer.jl b/models/_bench_perf_test/ER_row_buffer.jl
new file mode 100644
index 0000000000000000000000000000000000000000..bc6104dd4c65cdbebfad1ff17cb09588941766ca
--- /dev/null
+++ b/models/_bench_perf_test/ER_row_buffer.jl
@@ -0,0 +1,50 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+d=4
+k=3
+dict_var = Dict("E" => 1, "S" => 2, "ES" => 3, "P" => 4)
+dict_p = Dict("k1" => 1, "k2" => 2, "k3" => 3)
+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,
+           xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
+    a1 = p[1] * xn[1] * xn[2]
+    a2 = p[2] * xn[3]
+    a3 = p[3] * xn[3]
+    l_a = SVector(a1, a2, a3)
+    asum = sum(l_a)
+    l_nu = @SMatrix [-1 1 1;
+                     -1 1 0;
+                     1 -1 -1;
+                     0 0 1]
+    u1, u2 = rand(), rand()
+    tau = - log(u1) / asum
+    b_inf = 0.0
+    b_sup = a1
+    reaction = 0
+    for i = 1:3
+        if b_inf < asum*u2 < b_sup
+            reaction = i
+            break
+        end
+        b_inf += l_a[i]
+        b_sup += l_a[i+1]
+    end
+ 
+    nu = @view l_nu[:,reaction] # macro for avoiding a copy
+    for i = 1:4
+        mat_x[idx,i] = xn[i]+nu[i]
+    end
+    l_t[idx] = tn + tau
+    l_tr[idx] = "R$(reaction)"
+end
+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)
+export ER_row_buffer
+
diff --git a/tests/run_simulation.jl b/tests/run_simulation.jl
index 0336e5019d9a3c92d3dfa8c7894b56e199d9583e..3bc43c9fdeec0f49782c1156a2681ecfe2acff0e 100644
--- a/tests/run_simulation.jl
+++ b/tests/run_simulation.jl
@@ -11,5 +11,6 @@ if !isdir(str_dir_pics) mkdir(str_dir_pics) end
     @test include("simulation/sim_sir_col_buffer_bounded.jl")
     @test include("simulation/sim_sir_row_buffer_bounded.jl")
     @test include("simulation/sim_er.jl")
+    @test include("simulation/sim_er_row_buffer_bounded.jl")
 end
 
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 443aef4f0a5bb31da78dc23b3cfe22583217cb07..9934547b322a2158be3476c4f8c75832cbc74fdc 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -3,6 +3,7 @@ using Test
 
 @testset "Unit tests" begin
     @test include("unit/load_model.jl")
+    @test include("unit/load_model_bench.jl")
     @test include("unit/load_module.jl")
     @test include("unit/simulate_sir.jl")
     @test include("unit/simulate_sir_bounded.jl")
diff --git a/tests/simulation/sim_er_row_buffer_bounded.jl b/tests/simulation/sim_er_row_buffer_bounded.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6da83769f8beed6be145263dcb1354556451614f
--- /dev/null
+++ b/tests/simulation/sim_er_row_buffer_bounded.jl
@@ -0,0 +1,16 @@
+
+using MarkovProcesses 
+include(get_module_path() * "/core/_tests_simulate.jl")
+using PyPlot
+
+load_model("_bench_perf_test/ER_row_buffer")
+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.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er_row_buffer_bounded.png")
+plt.close()
+
+return true
+
diff --git a/tests/unit/load_model_bench.jl b/tests/unit/load_model_bench.jl
new file mode 100644
index 0000000000000000000000000000000000000000..223e8fe6d350fc24664ae15a68c9c4cd9ba1a142
--- /dev/null
+++ b/tests/unit/load_model_bench.jl
@@ -0,0 +1,12 @@
+
+using MarkovProcesses
+
+load_model("_bench_perf_test/SIR_col")
+load_model("_bench_perf_test/SIR_col_buffer")
+load_model("_bench_perf_test/SIR_row_buffer")
+load_model("_bench_perf_test/ER_col")
+load_model("_bench_perf_test/ER_col_buffer")
+load_model("_bench_perf_test/ER_row_buffer")
+
+return true
+