Reorganisation of the tests directory:

- a script to run the different categories
- now the execution of simulation/ scripts is checked
- several minor changes and add of tests

New directory : models/_bench_perf_tests.  Models for performance benchmarking purposes.

core/_tests_simulate.jl

+
+# 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
+

models/_bench_perf_test/SIR_col.jl

+
+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
+

models/_bench_perf_test/SIR_col_buffer.jl

+
+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
+

models/_bench_perf_test/SIR_row_buffer.jl

+
+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
+

tests/run_all.jl

+
+include("run_unit.jl")
+include("run_simulation.jl")
+

tests/run_simulation.jl

+
+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
+

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
+

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
+

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
+

tests/simulation/sim_sir_col_buffer_bounded.jl

+
+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
+