diff --git a/core/model.jl b/core/model.jl
index 8e8fa2eb5852cb662c4e1a171c0555a0186c0432..ed108270f69e28dcc1db055e5c97fda659acda8b 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -1,8 +1,67 @@
+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
diff --git a/core/observations.jl b/core/observations.jl
index 931f6bbd21f9900f5098e609e53db7fb88dbdfd4..86fca14a7de6e8f6382d3d263f5cc6710a899170 100644
--- a/core/observations.jl
+++ b/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
diff --git a/models/sir.jl b/models/sir.jl
index 9f960757f430a30baf0a43d57c712c49ae5f42c4..de72297c19c29f392692fcee202c303e6cff0dbd 100644
--- a/models/sir.jl
+++ b/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
diff --git a/tests/README b/tests/README
new file mode 100644
index 0000000000000000000000000000000000000000..b34713ceca6793f5bb4d746927678333e22f62fa
--- /dev/null
+++ b/tests/README
@@ -0,0 +1,7 @@
+
+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
+
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index ac173f135ed083dc2ddefa1e624bf5cc519ddc5c..9aee0767b25e17afcdba9fa691dbb46470330a2a 100644
--- a/tests/run_unit.jl
+++ b/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
diff --git a/tests/simulation/sim_sir.jl b/tests/simulation/sim_sir.jl
index 9601a670c94b606a72b739ce969030eb863935b6..cfad9a6b579514ff7c358870b8ea3c5675f36991 100644
--- a/tests/simulation/sim_sir.jl
+++ b/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()
diff --git a/tests/unit/simulate_sir.jl b/tests/unit/simulate_sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..a5b0c7560e5b9d1c4d84db949b3ae282c159dfe8
--- /dev/null
+++ b/tests/unit/simulate_sir.jl
@@ -0,0 +1,16 @@
+
+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
+