From c5f43847b3b29308ecbef1616bb21f468540126f Mon Sep 17 00:00:00 2001 From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr> Date: Sat, 14 Nov 2020 15:32:02 +0100 Subject: [PATCH] First commit with. Details the organization of the code. Two folders: - core/ which contains the essential files for the packages - tests/ which contains tests and benchmarks of other packages / methods For now on: - we wrote minimal tests in tests/ (a simulation of SIR model) that should run for a first minimal version of the package - we described in core/ the general structure of essential types, how they articulate each others and the minimum methods we should implement for a first version that runs the tests. --- core/MarkovProcesses.jl | 14 +++++++++++ core/model.jl | 10 ++++++++ core/observations.jl | 26 +++++++++++++++++++++ tests/load_module.jl | 3 +++ tests/models/sir.jl | 52 +++++++++++++++++++++++++++++++++++++++++ tests/sim_sir.jl | 12 ++++++++++ 6 files changed, 117 insertions(+) create mode 100644 core/MarkovProcesses.jl create mode 100644 core/model.jl create mode 100644 core/observations.jl create mode 100644 tests/load_module.jl create mode 100644 tests/models/sir.jl create mode 100644 tests/sim_sir.jl diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl new file mode 100644 index 0000000..e1f47fe --- /dev/null +++ b/core/MarkovProcesses.jl @@ -0,0 +1,14 @@ +module MarkovProcesses + +import Base: +, -, getfield, getindex + +export Model, ContinuousTimeModel, DiscreteTimeModel +export simulate, set_param!, get_param +include("model.jl") + +export Observations, AbstractTrajectory +export +,-,δ,get_obs_variables +include("observations.jl") + +end + diff --git a/core/model.jl b/core/model.jl new file mode 100644 index 0000000..4c093a8 --- /dev/null +++ b/core/model.jl @@ -0,0 +1,10 @@ + +abstract type Model end +abstract type ContinuousTimeModel <: Model end +abstract type DiscreteTimeModel <: Model 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 +function get_param(m::Model)::AbstractVector{Real} end + diff --git a/core/observations.jl b/core/observations.jl new file mode 100644 index 0000000..931f6bb --- /dev/null +++ b/core/observations.jl @@ -0,0 +1,26 @@ + +abstract type AbstractTrajectory end +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} +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 + diff --git a/tests/load_module.jl b/tests/load_module.jl new file mode 100644 index 0000000..548f0bc --- /dev/null +++ b/tests/load_module.jl @@ -0,0 +1,3 @@ + +using MarkovProcesses + diff --git a/tests/models/sir.jl b/tests/models/sir.jl new file mode 100644 index 0000000..c4db3e0 --- /dev/null +++ b/tests/models/sir.jl @@ -0,0 +1,52 @@ + +import StaticArrays: SVector, SMatrix, @SMatrix + +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"] +p = Parameters(0.0012, 0.05) +x0 = State(95, 5, 0) + +function f(xn::State, p::Parameters, tn::Real) + 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 0.0; + 1.0 -1.0; + 0.0 1.0] + + 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 = l_nu[:,reaction] + xnplus1 = State(xn[1]+nu[1], xn[2]+nu[2], xn[3]+nu[3]) + tnplus1 = tn + tau + transition = "R$(reaction)" + + return xnplus1 +end + +g = SVector(1) + +# Gamma should be constructed automatically in the case of + +#m = Model(d,dobs,k,dict,l_name,param,p,x0,f,g) + diff --git a/tests/sim_sir.jl b/tests/sim_sir.jl new file mode 100644 index 0000000..2feb900 --- /dev/null +++ b/tests/sim_sir.jl @@ -0,0 +1,12 @@ + +using MarkovProcesses +using PyPlot + +include("models/sir.jl") + +σ = simulate(SIR) +plt.figure() +plot(σ["S,times"], σ["S,values"]) +plt.savefig("sim_sir.png") +plt.close() + -- GitLab