- Add of ER model + tests

- Add of benchmark scripts that compares perf wrt pygmalion.
=> As expected it is bad in terms of reading cost (because of row by row
matrix read) but it is also not good in terms of simulation cost.
hcat seems to perform badly, should investigate.

bench/pygmalion/access_val_er.jl

+
+using BenchmarkTools
+using pygmalion
+using MarkovProcesses
+
+println("Pygmalion:")
+
+str_m = "enzymatic_reaction"
+str_d = "abc_er"
+pygmalion.load_model(str_m)
+str_oml = "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)
+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 random_trajectory_value_pyg(so::SystemObservation)
+    n_states = get_number_of_observations(so, "P")
+    return to_vec(so, "P", rand(1:n_states))
+end
+# Bench
+@timev random_trajectory_value_pyg(so)
+b1_pyg = @benchmark random_trajectory_value_pyg($so)
+@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
+
+println("MarkovProcesses:")
+
+MarkovProcesses.load_model("ER")
+ER.time_bound = 10.0
+σ = MarkovProcesses.simulate(ER)
+function random_trajectory_value(σ::AbstractTrajectory)
+    n_states = get_states_number(σ)
+    return σ["P"][rand(1:n_states)]
+end
+# Bench
+@timev random_trajectory_value(σ) 
+b1 = @benchmark random_trajectory_value($σ)
+@show minimum(b1), mean(b1), maximum(b1)
+

bench/pygmalion/access_var_values.jl

+
+using BenchmarkTools
+using pygmalion
+using MarkovProcesses
+
+println("Pygmalion:")
+
+str_m = "enzymatic_reaction"
+str_d = "abc_er"
+pygmalion.load_model(str_m)
+str_oml = "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)
+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)
+    vals = to_vec(so, "P")
+    n_states = length(vals)
+    res = 0.0
+    for i = 1:n_states
+        res += vals[i]
+    end
+    return res
+end
+function read_trajectory_v2(so::SystemObservation)
+    idx_P = om_findfirst("P", so.oml)
+    n_states = length(so.otll[idx_P]) 
+    res = 0.0
+    for i = 1:n_states
+        res += so.otll[idx_P][i][2][1]
+    end
+    return res
+end
+# Bench
+@timev read_trajectory_v1(so)
+b1_pyg = @benchmark read_trajectory_v1($so) 
+@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
+@timev read_trajectory_v2(so)
+b2_pyg = @benchmark read_trajectory_v2($so) 
+@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
+
+println("MarkovProcesses:")
+
+MarkovProcesses.load_model("ER")
+ER.time_bound = 10.0
+σ = MarkovProcesses.simulate(ER)
+function read_trajectory(σ::AbstractTrajectory)
+    n_states = get_states_number(σ)
+    res = 0
+    for i = 1:n_states
+        res += (σ["P"])[i]
+    end
+    return res
+end
+# Bench
+@timev read_trajectory(σ) 
+b1 = @benchmark read_trajectory($σ)
+@show minimum(b1), mean(b1), maximum(b1)
+

bench/pygmalion/long_sim_er.jl

+
+using BenchmarkTools
+
+println("Pygmalion:")
+
+using pygmalion
+str_m = "enzymatic_reaction"
+str_d = "abc_er"
+pygmalion.load_model(str_m)
+str_oml = "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)
+u = Control(10.0)
+tml = 1:400
+g_all = create_observation_function([ObserverModel(str_oml, tml)]) 
+b1_pyg = @benchmark pygmalion.simulate($f, $g_all, $x0, $u, $p_true; on = nothing, full_timeline = true)
+b2_pyg = @benchmark pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
+
+@timev pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
+@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
+@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
+
+println("MarkovProcesses:")
+
+using MarkovProcesses
+MarkovProcesses.load_model("ER")
+ER.time_bound = 10.0
+b1 = @benchmark MarkovProcesses.simulate($ER)
+b2 = @benchmark MarkovProcesses.simulate(ER)
+
+@timev MarkovProcesses.simulate(ER)
+@show minimum(b1), mean(b1), maximum(b1)
+@show minimum(b2), mean(b2), maximum(b2)
+

core/MarkovProcesses.jl

@@ -8,7 +8,7 @@ export load_model
 include("model.jl")
 
 export Observations, AbstractTrajectory
-export +,-,δ,get_obs_variables
+export +,-,δ,get_obs_variables,get_states_number
 include("observations.jl")
 
 end

core/model.jl

@@ -38,7 +38,7 @@ end
 function simulate(m::ContinuousTimeModel)
     full_values = zeros(m.d, 0)
     times = zeros(0)
-    transitions = Vector{String}(undef, 0)
+    transitions = Vector{Union{String,Nothing}}(undef, 0)
     xn = m.x0
     tn = m.t0 
     n = 0
@@ -51,6 +51,14 @@ function simulate(m::ContinuousTimeModel)
         n += 1
     end
     values = 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
 
@@ -62,6 +70,7 @@ function simulate(m::ContinuousTimeModel, n::Int)
     return obs
 end
 
+is_bounded(m::Model) = m.time_bound < Inf
 function check_consistency(m::Model) end
 function simulate(m::Model, n::Int; bound::Real = Inf)::AbstractObservations end
 function set_param!(m::Model, p::AbstractVector{Real})::Nothing end

core/observations.jl

@@ -6,7 +6,7 @@ struct Trajectory <: AbstractTrajectory
     m::ContinuousTimeModel
     values::AbstractMatrix{Real}
     times::AbstractVector{Real}
-    transitions::AbstractVector{Union{String,Missing}}
+    transitions::AbstractVector{Union{String,Nothing}}
 end
 
 function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
@@ -15,7 +15,10 @@ function δ(σ1::AbstractTrajectory,t::Real) end
 function get_obs_variables(σ::AbstractTrajectory) end
 
 get_values(σ::AbstractTrajectory, var::String) = 
-    σ.values[(σ.m)._map_obs_var_idx[var],:] 
+σ.values[(σ.m)._map_obs_var_idx[var],:] 
+
+get_states_number(σ::AbstractTrajectory) =
+length(σ.times)
 
 function getindex(σ::AbstractTrajectory, idx::String)
     if idx  == "times"

models/ER.jl

+
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+State = SVector{4, Int}
+Parameters = SVector{3, Real}
+
+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 = Parameters(1.0, 1.0, 1.0)
+x0 = State(100, 100, 0, 0)
+t0 = 0.0
+
+function f(xn::State, tn::Real, p::Parameters)
+    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 = l_nu[:,reaction]
+    xnplus1 = State(xn[1]+nu[1], xn[2]+nu[2], xn[3]+nu[3], xn[4]+nu[4])
+    tnplus1 = tn + tau
+    transition = "R$(reaction)"
+    
+    return xnplus1, tnplus1, transition
+end
+is_absorbing_sir(p::Parameters,xn::State) = 
+    (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) == 0.0
+g = SVector("P")
+
+ER = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f,is_absorbing_sir; g=g)
+export ER
+

models/SIR.jl

@@ -18,9 +18,9 @@ function f(xn::State, tn::Real, p::Parameters)
     l_a = SVector(a1, a2)
     asum = sum(l_a)
     # column-major order
-    l_nu = @SMatrix [-1.0 0.0;
-                     1.0 -1.0;
-                     0.0 1.0]
+    l_nu = @SMatrix [-1 0;
+                     1 -1;
+                     0 1]
     
     u1, u2 = rand(), rand()
     tau = - log(u1) / asum

tests/run_unit.jl

@@ -5,5 +5,7 @@ using Test
     @test include("unit/load_model.jl")
     @test include("unit/load_module.jl")
     @test include("unit/simulate_sir.jl")
+    @test include("unit/simulate_sir_bounded.jl")
+    @test include("unit/simulate_er.jl")
 end
 

tests/simulation/sim_er.jl

+
+using MarkovProcesses 
+using PyPlot
+
+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.close()
+

tests/simulation/sim_sir_bounded.jl

+
+using MarkovProcesses 
+using PyPlot
+
+load_model("SIR")
+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.close()
+

tests/unit/simulate_er.jl

+
+using MarkovProcesses
+
+load_model("ER")
+
+σ = simulate(ER)
+
+d1 = Dict("E" => 1, "S"=> 2, "ES" => 3, "P" => 4)
+d2 = Dict("P" => 1)
+
+bool_test = ER.g == ["P"] && ER._g_idx == [4] && 
+            ER.map_var_idx == d1 && 
+            ER._map_obs_var_idx == d2
+
+return bool_test
+

tests/unit/simulate_sir_bounded.jl

+
+using MarkovProcesses
+
+load_model("SIR")
+SIR.time_bound = 100.0
+
+σ = simulate(SIR)
+
+d1 = Dict("S" => 1, "I"=> 2, "R" => 3)
+d2 = Dict("I" => 1)
+
+bool_test = SIR.g == ["I"] && SIR._g_idx == [2] && 
+            SIR.map_var_idx == d1 && 
+            SIR._map_obs_var_idx == d2
+
+return bool_test
+