diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
new file mode 100644
index 0000000000000000000000000000000000000000..0f55f06e852ea5bd27bd5f5a86668cd9b0a06e71
--- /dev/null
+++ b/core/_tests_simulate.jl
@@ -0,0 +1,158 @@
+
+# File for benchmarking simulation and memory access of the package.
+
+_get_values_col(σ::AbstractTrajectory, var::String) =
+σ.values[(σ.m)._map_obs_var_idx[var],:]
+_get_values_row(σ::AbstractTrajectory, var::String) =
+σ.values[:,(σ.m)._map_obs_var_idx[var]]
+
+
+function _simulate_col(m::ContinuousTimeModel)
+ # trajectory fields
+ full_values = zeros(m.d, 0)
+ times = zeros(0)
+ transitions = Vector{Union{String,Nothing}}(undef,0)
+ # values at time n
+ n = 0
+ xn = m.x0
+ tn = m.t0
+ tr = [""]
+ # at time n+1
+ xnplus1 = zeros(Int, m.d)
+ tnplus1 = zeros(Float64, 1)
+ is_absorbing = (m.is_absorbing(m.p,xn))::Bool
+ while !is_absorbing && (tn <= m.time_bound)
+ m.f!(xnplus1, tnplus1, tr, xn, tn, m.p)
+ full_values = hcat(full_values, xnplus1)
+ push!(times, tnplus1[1])
+ push!(transitions, tr[1])
+ xn, tn = xnplus1, tnplus1[1]
+ n += 1
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ values = @view full_values[m._g_idx,:]
+ if is_bounded(m)
+ if times[end] > m.time_bound
+ values[:, end] = values[:,end-1]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ end
+ end
+ return Trajectory(m, values, times, transitions)
+end
+
+function _simulate_row(m::ContinuousTimeModel)
+ # trajectory fields
+ full_values = zeros(m.d, 0)
+ times = zeros(0)
+ transitions = Vector{Union{String,Nothing}}(undef,0)
+ # values at time n
+ n = 0
+ xn = m.x0
+ tn = m.t0
+ tr = [""]
+ # at time n+1
+ xnplus1 = zeros(Int, m.d)
+ tnplus1 = zeros(Float64, 1)
+ is_absorbing = (m.is_absorbing(m.p,xn))::Bool
+ while !is_absorbing && (tn <= m.time_bound)
+ m.f!(xnplus1, tnplus1, tr, xn, tn, m.p)
+ full_values = vcat(full_values, xnplus1)
+ push!(times, tnplus1[1])
+ push!(transitions, tr[1])
+ xn, tn = xnplus1, tnplus1[1]
+ n += 1
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ values = @view full_values[m._g_idx,:]
+ if is_bounded(m)
+ if times[end] > m.time_bound
+ values[:, end] = values[:,end-1]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ end
+ end
+ return Trajectory(m, values, times, transitions)
+end
+
+
+function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
+ # trajectory fields
+ full_values = zeros(m.d, 0)
+ times = zeros(0)
+ transitions = Vector{Union{String,Nothing}}(undef,0)
+ # values at time n
+ n = 0
+ xn = m.x0
+ tn = m.t0
+ # at time n+1
+ mat_x = zeros(Int, m.d, buffer_size)
+ l_t = zeros(Float64, buffer_size)
+ l_tr = Vector{String}(undef, buffer_size)
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ while !is_absorbing && (tn <= m.time_bound)
+ i = 0
+ while i < 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]
+ tn = l_t[i]
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ full_values = hcat(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
+ end
+ values = @view full_values[m._g_idx,:]
+ if is_bounded(m)
+ if times[end] > m.time_bound
+ values[:, end] = values[:,end-1]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ end
+ end
+ return Trajectory(m, values, times, transitions)
+end
+
+function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
+ # trajectory fields
+ full_values = zeros(0, m.d)
+ times = zeros(0)
+ transitions = Vector{Union{String,Nothing}}(undef,0)
+ # values at time n
+ n = 0
+ xn = m.x0
+ tn = m.t0
+ # at time n+1
+ mat_x = zeros(Int, buffer_size, m.d)
+ l_t = zeros(Float64, buffer_size)
+ l_tr = Vector{String}(undef, buffer_size)
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ while !is_absorbing && (tn <= m.time_bound)
+ i = 0
+ while i < 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]
+ tn = l_t[i]
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ 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
+ is_absorbing = m.is_absorbing(m.p,xn)::Bool
+ end
+ values = @view full_values[:,m._g_idx]
+ if is_bounded(m)
+ if times[end] > m.time_bound
+ values[end,:] = values[end-1,:]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ end
+ end
+ return Trajectory(m, values, times, transitions)
+end
+
diff --git a/models/_bench_perf_test/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl
new file mode 100644
index 0000000000000000000000000000000000000000..00cf2e48e53862d49bb00f4310e1cffe142151dc
--- /dev/null
+++ b/models/_bench_perf_test/SIR_col.jl
@@ -0,0 +1,49 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+d=3
+k=2
+dict_var = Dict("S" => 1, "I" => 2, "R" => 3)
+dict_p = Dict("ki" => 1, "kr" => 2)
+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},
+ xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+ a1 = p[1] * xn[1] * xn[2]
+ a2 = p[2] * xn[2]
+ l_a = SVector(a1, a2)
+ asum = sum(l_a)
+ # column-major order
+ l_nu = @SMatrix [-1 0;
+ 1 -1;
+ 0 1]
+
+ u1, u2 = rand(), rand()
+ tau = - log(u1) / asum
+ b_inf = 0.0
+ b_sup = a1
+ reaction = 0
+ for i = 1:2
+ 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
+ xnplus1[1] = xn[1]+nu[1]
+ xnplus1[2] = xn[2]+nu[2]
+ xnplus1[3] = xn[3]+nu[3]
+ 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
+g = ["I"]
+
+SIR_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_sir; 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
new file mode 100644
index 0000000000000000000000000000000000000000..d37828bdbb2eee8145513b1d2217afe11749294d
--- /dev/null
+++ b/models/_bench_perf_test/SIR_col_buffer.jl
@@ -0,0 +1,49 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+d=3
+k=2
+dict_var = Dict("S" => 1, "I" => 2, "R" => 3)
+dict_p = Dict("ki" => 1, "kr" => 2)
+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,
+ xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
+ a1 = p[1] * xn[1] * xn[2]
+ a2 = p[2] * xn[2]
+ l_a = SVector(a1, a2)
+ asum = sum(l_a)
+ # column-major order
+ l_nu = @SMatrix [-1 0;
+ 1 -1;
+ 0 1]
+
+ u1, u2 = rand(), rand()
+ tau = - log(u1) / asum
+ b_inf = 0.0
+ b_sup = a1
+ reaction = 0
+ for i = 1:2
+ 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:3
+ mat_x[i,idx] = xn[i]+nu[i]
+ end
+ 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
+g = ["I"]
+
+SIR_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_sir; 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
new file mode 100644
index 0000000000000000000000000000000000000000..366b68b9d120fc8368d32bacdd810af87b8ea3d8
--- /dev/null
+++ b/models/_bench_perf_test/SIR_row_buffer.jl
@@ -0,0 +1,49 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+d=3
+k=2
+dict_var = Dict("S" => 1, "I" => 2, "R" => 3)
+dict_p = Dict("ki" => 1, "kr" => 2)
+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,
+ xn::AbstractVector{Int}, tn::Float64, p::Vector{Float64})
+ a1 = p[1] * xn[1] * xn[2]
+ a2 = p[2] * xn[2]
+ l_a = SVector(a1, a2)
+ asum = sum(l_a)
+ # column-major order
+ l_nu = @SMatrix [-1 0;
+ 1 -1;
+ 0 1]
+
+ u1, u2 = rand(), rand()
+ tau = - log(u1) / asum
+ b_inf = 0.0
+ b_sup = a1
+ reaction = 0
+ for i = 1:2
+ 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:3
+ mat_x[idx,i] = xn[i]+nu[i]
+ end
+ 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
+g = ["I"]
+
+SIR_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f!,is_absorbing_sir; g=g)
+export SIR_row_buffer
+
diff --git a/tests/run_all.jl b/tests/run_all.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d8c11f700af560185a8f28890fe4a1e262b4b8d9
--- /dev/null
+++ b/tests/run_all.jl
@@ -0,0 +1,4 @@
+
+include("run_unit.jl")
+include("run_simulation.jl")
+
diff --git a/tests/run_simulation.jl b/tests/run_simulation.jl
new file mode 100644
index 0000000000000000000000000000000000000000..224820fe1b65b18a6d63c2dddd4fa5b7bc44fe03
--- /dev/null
+++ b/tests/run_simulation.jl
@@ -0,0 +1,14 @@
+
+using Test
+import MarkovProcesses: get_module_path
+
+str_dir_pics = get_module_path() * "/tests/simulation/res_pics"
+if !isdir(str_dir_pics) mkdir(str_dir_pics) end
+
+@testset "Simulation tests" begin
+ @test include("simulation/sim_sir.jl")
+ @test include("simulation/sim_sir_bounded.jl")
+ @test include("simulation/sim_sir_col_buffer_bounded.jl")
+ @test include("simulation/sim_er.jl")
+end
+
diff --git a/tests/simulation/sim_er.jl b/tests/simulation/sim_er.jl
index f46aaacbdf2bf477e026638f4834b02eb118387a..09d4575b3867d5c5121448031d15717e3a14ed05 100644
--- a/tests/simulation/sim_er.jl
+++ b/tests/simulation/sim_er.jl
@@ -7,6 +7,8 @@ load_model("ER")
σ = simulate(ER)
plt.figure()
plt.step(σ["times"], σ["P"], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig("sim_er.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er.png", dpi=480)
plt.close()
+return true
+
diff --git a/tests/simulation/sim_sir.jl b/tests/simulation/sim_sir.jl
index dca855fa53c4dc48d0fb2ab7cef2871407bd9e49..2a1f00fa21f22f798834fbfd7b4f8edc98ca358b 100644
--- a/tests/simulation/sim_sir.jl
+++ b/tests/simulation/sim_sir.jl
@@ -7,6 +7,8 @@ load_model("SIR")
σ = simulate(SIR)
plt.figure()
plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig("sim_sir.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir.png")
plt.close()
+return true
+
diff --git a/tests/simulation/sim_sir_bounded.jl b/tests/simulation/sim_sir_bounded.jl
index bbca0535d2090101c912ff99167c1d499a7e4f99..95d64343b1769a9d6bf7d2b9cde10d7f913dfacc 100644
--- a/tests/simulation/sim_sir_bounded.jl
+++ b/tests/simulation/sim_sir_bounded.jl
@@ -8,6 +8,8 @@ SIR.time_bound = 100.0
σ = simulate(SIR)
plt.figure()
plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig("sim_sir_bounded.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_bounded.png")
plt.close()
+return true
+
diff --git a/tests/simulation/sim_sir_col_buffer_bounded.jl b/tests/simulation/sim_sir_col_buffer_bounded.jl
new file mode 100644
index 0000000000000000000000000000000000000000..94335f9a0ab35e28bc81d90d2c602c9ad900b0fa
--- /dev/null
+++ b/tests/simulation/sim_sir_col_buffer_bounded.jl
@@ -0,0 +1,16 @@
+
+using MarkovProcesses
+include(get_module_path() * "/core/_tests_simulate.jl")
+using PyPlot
+
+load_model("_bench_perf_test/SIR_col_buffer")
+SIR_col_buffer.time_bound = 100.0
+
+σ = _simulate_col_buffer(SIR_col_buffer)
+plt.figure()
+plt.step(σ["times"], σ["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()
+
+return true
+