From 06dca928f1b9e20d2a3b020a4e52ced099923a63 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Sun, 15 Nov 2020 17:43:10 +0100
Subject: [PATCH] - 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     | 41 ++++++++++++++++++
 bench/pygmalion/access_var_values.jl | 62 ++++++++++++++++++++++++++++
 bench/pygmalion/long_sim_er.jl       | 35 ++++++++++++++++
 core/MarkovProcesses.jl              |  2 +-
 core/model.jl                        | 11 ++++-
 core/observations.jl                 |  7 +++-
 models/ER.jl                         | 54 ++++++++++++++++++++++++
 models/SIR.jl                        |  6 +--
 tests/run_unit.jl                    |  2 +
 tests/simulation/sim_er.jl           | 12 ++++++
 tests/simulation/sim_sir_bounded.jl  | 13 ++++++
 tests/unit/simulate_er.jl            | 16 +++++++
 tests/unit/simulate_sir_bounded.jl   | 17 ++++++++
 13 files changed, 271 insertions(+), 7 deletions(-)
 create mode 100644 bench/pygmalion/access_val_er.jl
 create mode 100644 bench/pygmalion/access_var_values.jl
 create mode 100644 bench/pygmalion/long_sim_er.jl
 create mode 100644 models/ER.jl
 create mode 100644 tests/simulation/sim_er.jl
 create mode 100644 tests/simulation/sim_sir_bounded.jl
 create mode 100644 tests/unit/simulate_er.jl
 create mode 100644 tests/unit/simulate_sir_bounded.jl

diff --git a/bench/pygmalion/access_val_er.jl b/bench/pygmalion/access_val_er.jl
new file mode 100644
index 0000000..05c444c
--- /dev/null
+++ b/bench/pygmalion/access_val_er.jl
@@ -0,0 +1,41 @@
+
+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)
+
diff --git a/bench/pygmalion/access_var_values.jl b/bench/pygmalion/access_var_values.jl
new file mode 100644
index 0000000..a108da2
--- /dev/null
+++ b/bench/pygmalion/access_var_values.jl
@@ -0,0 +1,62 @@
+
+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)
+
diff --git a/bench/pygmalion/long_sim_er.jl b/bench/pygmalion/long_sim_er.jl
new file mode 100644
index 0000000..6a066ba
--- /dev/null
+++ b/bench/pygmalion/long_sim_er.jl
@@ -0,0 +1,35 @@
+
+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)
+
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index eec79d0..04ac4c7 100644
--- a/core/MarkovProcesses.jl
+++ b/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
diff --git a/core/model.jl b/core/model.jl
index ed10827..5e87907 100644
--- a/core/model.jl
+++ b/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
diff --git a/core/observations.jl b/core/observations.jl
index 86fca14..7a35b81 100644
--- a/core/observations.jl
+++ b/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"
diff --git a/models/ER.jl b/models/ER.jl
new file mode 100644
index 0000000..876b129
--- /dev/null
+++ b/models/ER.jl
@@ -0,0 +1,54 @@
+
+
+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
+
diff --git a/models/SIR.jl b/models/SIR.jl
index de72297..122889b 100644
--- a/models/SIR.jl
+++ b/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
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 9aee076..e85df99 100644
--- a/tests/run_unit.jl
+++ b/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
 
diff --git a/tests/simulation/sim_er.jl b/tests/simulation/sim_er.jl
new file mode 100644
index 0000000..f46aaac
--- /dev/null
+++ b/tests/simulation/sim_er.jl
@@ -0,0 +1,12 @@
+
+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()
+
diff --git a/tests/simulation/sim_sir_bounded.jl b/tests/simulation/sim_sir_bounded.jl
new file mode 100644
index 0000000..bbca053
--- /dev/null
+++ b/tests/simulation/sim_sir_bounded.jl
@@ -0,0 +1,13 @@
+
+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()
+
diff --git a/tests/unit/simulate_er.jl b/tests/unit/simulate_er.jl
new file mode 100644
index 0000000..32340d1
--- /dev/null
+++ b/tests/unit/simulate_er.jl
@@ -0,0 +1,16 @@
+
+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
+
diff --git a/tests/unit/simulate_sir_bounded.jl b/tests/unit/simulate_sir_bounded.jl
new file mode 100644
index 0000000..c0b68e6
--- /dev/null
+++ b/tests/unit/simulate_sir_bounded.jl
@@ -0,0 +1,17 @@
+
+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
+
-- 
GitLab