Commit that implement the essentials for the lauch of a simulation of a

SIR model.

- Evolution of the ContinuousTimeModel type
- simulate, easy access to Trajectory fields
- fix of several errors

core/model.jl

 
+import StaticArrays: SVector
+
 abstract type Model end 
 abstract type ContinuousTimeModel <: Model end 
 abstract type DiscreteTimeModel <: Model end 
 
+mutable struct CTMC <: ContinuousTimeModel
+    d::Int # state space dim
+    k::Int # parameter space dim
+    map_var_idx::Dict # maps str to full state space
+    _map_obs_var_idx::Dict # maps str to observed state space
+    map_param_idx::Dict # maps str in parameter space
+    l_name_transitions::AbstractVector{String}
+    p::AbstractVector{Real}
+    x0::AbstractVector{Int}
+    t0::Real
+    f::Function
+    g::AbstractVector{String} # of dimension dobs
+    _g_idx::AbstractVector{Int} # of dimension dobs
+    is_absorbing::Function
+    time_bound::Real
+end
+
+function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::AbstractVector{String}, 
+              p::AbstractVector{Real}, x0::AbstractVector{Int}, t0::Real, 
+              f::Function, is_absorbing::Function; g::AbstractVector{String} = keys(map_var_idx), time_bound::Real = Inf)
+    dobs = length(g)
+    _map_obs_var_idx = Dict()
+    _g_idx = Vector{Int}(undef, dobs)
+    for i = 1:dobs
+        _g_idx[i] = map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::String => idx in state space )
+        _map_obs_var_idx[g[i]] = i
+    end
+    return CTMC(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f, g, _g_idx, is_absorbing, time_bound)
+end
+
+function simulate(m::ContinuousTimeModel)
+    full_values = zeros(m.d, 0)
+    times = zeros(0)
+    transitions = Vector{String}(undef, 0)
+    xn = m.x0
+    tn = m.t0 
+    n = 0
+    while !m.is_absorbing(m.p,xn) && tn <= m.time_bound
+        xnplus1, tnplus1, tr = f(xn, tn, m.p)
+        full_values = hcat(full_values, xnplus1)
+        push!(times, tnplus1)
+        push!(transitions, tr)
+        xn, tn = xnplus1, tnplus1
+        n += 1
+    end
+    values = full_values[m._g_idx,:] 
+    return Trajectory(m, values, times, transitions)
+end
+
+function simulate(m::ContinuousTimeModel, n::Int)
+    obs = ContinuousObservations(undef, n)
+    for i = 1:n
+        obs[i] = simulate(m)
+    end
+    return obs
+end
+
 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

@@ -5,22 +5,25 @@ ContinuousObservations = AbstractVector{AbstractTrajectory}
 struct Trajectory <: AbstractTrajectory
     m::ContinuousTimeModel
     values::AbstractMatrix{Real}
-    times::AbstractMatrix{Real}
-    transitions::AbstractVector{Union{String,Missing,Nothing}}
-end
-
-struct ObservedTrajectory <: AbstractTrajectory
-    m::ContinuousTimeModel
-    values::AbstractMatrix{Real}
-    times::AbstractMatrix{Real}
+    times::AbstractVector{Real}
+    transitions::AbstractVector{Union{String,Missing}}
 end
 
 function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
 function -(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
 function δ(σ1::AbstractTrajectory,t::Real) end
-function get_obs_variables(σ::Trajectory) end
-function get_obs_variables(σ::ObservedTrajectory) end
-function get_values(σ::AbstractTrajectory, variable::String) end
-function get_times(σ::AbstractTrajectory, variable::String) end
-function getindex(σ::AbstractTrajectory, idx::String) end
+function get_obs_variables(σ::AbstractTrajectory) end
+
+get_values(σ::AbstractTrajectory, var::String) = 
+    σ.values[(σ.m)._map_obs_var_idx[var],:] 
+
+function getindex(σ::AbstractTrajectory, idx::String)
+    if idx  == "times"
+        return σ.times
+    elseif idx == "transitions"
+        return σ.transitions
+    else
+        return get_values(σ, idx)
+    end
+end
 

models/sir.jl

@@ -5,14 +5,14 @@ State = SVector{3, Int}
 Parameters = SVector{2, Real}
 
 d=3
-dobs=1
 k=2
-dict=Dict("S" => 1, "I" => 2, "R" => 3)
-l_name_param = ["ki", "kr"]
+dict_var = Dict("S" => 1, "I" => 2, "R" => 3)
+dict_p = Dict("ki" => 1, "kr" => 2)
+l_tr = ["R1","R2"]
 p = Parameters(0.0012, 0.05)
 x0 = State(95, 5, 0)
-
-function f(xn::State, p::Parameters, tn::Real)
+t0 = 0.0
+function f(xn::State, tn::Real, p::Parameters)
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
     l_a = SVector(a1, a2)
@@ -27,7 +27,7 @@ function f(xn::State, p::Parameters, tn::Real)
     b_inf = 0.0
     b_sup = a1
     reaction = 0
-    for i = 1:2
+    for i = 1:2 
         if b_inf < asum*u2 < b_sup
             reaction = i
             break
@@ -41,12 +41,11 @@ function f(xn::State, p::Parameters, tn::Real)
     tnplus1 = tn + tau
     transition = "R$(reaction)"
 
-    return xnplus1
+    return xnplus1, tnplus1, transition
 end
+is_absorbing_sir(p::Parameters,xn::State) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) == 0.0
+g = SVector("I")
 
-g = SVector(1)
-
-# Gamma should be constructed automatically in the case of 
-
-#SIR = Model(d,dobs,k,dict,l_name,param,p,x0,f,g)
+SIR = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f,is_absorbing_sir; g=g)
+export SIR
 

tests/README

+
+unit/: ensures essential commands execute well, a test can contain a simulate() method to check it is well launched
+simulation/: tries to check the consistency of the produced simulations, comparison with other packages,
+             quantify if the simulation is probable
+
+Lauch unit tests: julia run_unit.jl
+

tests/run_unit.jl

@@ -4,5 +4,6 @@ using Test
 @testset "Unit tests" begin
     @test include("unit/load_model.jl")
     @test include("unit/load_module.jl")
+    @test include("unit/simulate_sir.jl")
 end
 

tests/simulation/sim_sir.jl

@@ -6,7 +6,7 @@ load_model("sir")
 
 σ = simulate(SIR)
 plt.figure()
-plot(σ["S,times"], σ["S,values"])
+plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
 plt.savefig("sim_sir.png")
 plt.close()
 

tests/unit/simulate_sir.jl

+
+using MarkovProcesses
+
+load_model("sir")
+
+σ = 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
+