Skip to content
Snippets Groups Projects
Commit 7b9f1841 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud Committed by Mahmoud Bentriou
Browse files

Add of new documentation pages + only doc build for ci for now

parent c3f95163
No related branches found
No related tags found
No related merge requests found
stages:
- build-and-test
# - build-and-test
- deploy
test:latest:
......
# Model test related methods
```@autodocs
Modules = [MarkovProcesses]
Pages = ["src/model.jl"]
```
# Plots related methods
```@autodocs
Modules = [MarkovProcesses]
Pages = ["src/plots.jl"]
```
# Trajectory related methods
```@autodocs
Modules = [MarkovProcesses]
Pages = ["src/trajectory.jl"]
```
docs/src/assets/logo.png

2.98 KiB

docs/src/assets/sir_trajectory.png

132 KiB

# Create a model
The package offers different ways to create models based on CRNs.
## Load a pre-written model
A bunch of models are already writtne within the package. If `str_model::String` is the name of an implemented model, then `load_model(str_model)` creates a variable with name `str_model`.
```julia
load_model("poisson")
```
Available models are listed below.
## Define a Chemical Reaction Network
Let's consider the Chemical Reaction Network of the SIR model:
``
Infection: S + I \xrightarrow{k_i} 2I \\
Recovery: I \xrightarrow{k_r} R
``
The macro `@network_model` creates easily a CTMC stored in a `ContinuousTimeModel` variable based on this formalism.
```julia
julia> easy_sir = @network_model begin
Infection: (S + I => 2I, ki*I*S)
Recovery: (I => R, kr*I)
end "My awesome SIR"
My_awesome_SIRModel <: ContinuousTimeModel model
- variables :
* I (index = 2 in state space)
* R (index = 3 in state space)
* S (index = 1 in state space)
- parameters :
* ki (index = 1 in parameter space)
* kr (index = 2 in parameter space)
- transitions : Infection,Recovery
- observed variables :
* S (index = 1 in observed state space, index = 1 in state space)
* I (index = 2 in observed state space, index = 2 in state space)
* R (index = 3 in observed state space, index = 3 in state space)
p = [0.0, 0.0]
x0 = [0, 0, 0]
t0 = 0.0
time bound = Inf
```
In the first reaction, `ki*I*S` is the reaction rate of the reaction `Infection`. This model is almost ready to use, we have to set the initial state and the parameters.
```julia
julia> set_param!(easy_sir, [0.0012, 0.05])
set_x0!(easy_sir, [95, 5, 0])
σ = simulate(easy_sir)
load_plots()
plot(σ)
```
![Plot of a simulated SIR trajectory](assets/sir_trajectory.png)
## Manually (advanced)
This page is intented to advanced uses of the package, in order to use
### Based on an existing model
### From scratch
When the above cases don't fit your application one can create manually a `ContinuousTimeModel`. Let's take a look about the signature of the constructor method:
```julia
function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict{VariableModel,Int},
map_param_idx::Dict{ParameterModel,Int}, transitions::Vector{<:Transition},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, isabsorbing::Function; kwargs)
```
Let's construct an SIR model manually. First, one has to specify the dimensions of the state space and the parameter space.
```julia
dim_state_sir, dim_params_sir = 3, 2
```
`map_var_idx` is a dictionary that maps each model variable (represented by a `Symbol`) to an index in the state space.
```julia
map_var_idx_sir = Dict(:S => 1, :I => 2, :R => 3)
```
`map_var_params` is the equivalent of `map_var_idx` for parameters.
```julia
map_params_idx_sir = Dict(:ki => 1, :kr => 2)
```
`transitions` are the transitions/reactions of the model (vector of `Symbol`), `p`, `x0` and `t0` are respectively the parameters, the initial state and initial time of the model.
```julia
transitions_sir = [:Infection, :Recovery]
p_sir = [0.0012, 0.05]
x0_sir = [95, 5, 0]
t0_sir = 0.0
```
The two last arguments are functions, the first one, called `f!` must have the signature:
```julia
function f!(xnplus1::Vector{Int}, ptr_t::Vector{Float64}, ptr_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
```
It should return nothing. `xnplus1`, `ptr_t` and `ptr_tr` are vectors where the next values are stored. `ptr_t` is of length 1 and stores the next time value (`ptr_t[1] = tn + delta_t`) whereas `ptr_tr` stores the name of the next transition/reaction (`ptr_tr[1] = :Infection` for example). This function is implemented in the package as:
```julia
@everywhere function sir_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
@inbounds a1 = p[1] * xn[1] * xn[2]
@inbounds a2 = p[2] * xn[2]
l_a = (a1, a2)
asum = sum(l_a)
if asum == 0.0
copyto!(xnplus1, xn)
return nothing
end
nu_1 = (-1, 1, 0)
nu_2 = (0, -1, 1)
l_nu = (nu_1, nu_2)
l_str_R = (:Infection, :Recovery)
u1 = rand()
u2 = 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
@inbounds b_inf += l_a[i]
@inbounds b_sup += l_a[i+1]
end
nu = l_nu[reaction]
for i = 1:3
@inbounds xnplus1[i] = xn[i]+nu[i]
end
@inbounds l_t[1] = tn + tau
@inbounds l_tr[1] = l_str_R[reaction]
end
```
i
The second function called `isaborbing` must have the signature:
```julia
isabsorbing(p::Vector{Float64}, xn::Vector{Int})
```
This function checks if the state `xn` is an absorbing state according to the model parametrised by `p`. It has to return true or false.
For a CTMC, a state is an absorbing state if the total exit rate is zero. In the case of the SIR model;
```julia
@everywhere sir_isabsorbing(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
```
Finally one sets the observed variables and the model can be created. The following lines creates a new type `TryhardSIRModel <: ContinuousTimeModel`, and the core of simulation.
```julia
g_sir = [:I]
# Generates simulate method for the new model
@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:TryhardSIRModel))
@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:TryhardSIRModel))
@everywhere @eval $(MarkovProcesses.generate_code_simulation(:TryhardSIRModel, :sir_f!, :sir_isabsorbing))
tryhard_sir = TryhardSIRModel(dim_state_sir, dim_params_sir,
map_var_idx_sir, map_params_idx_sir,
transitions_sir, p_sir, x0_sir, t0_sir,
:sir_f!, :sir_isabsorbing; g = g_sir)
σ = simulate(tryhard_sir)
```
## List of pre-written models
- `load_model("poisson")`: Poisson process
- `load_model("ER")`: Michaelis-Menten kinetics (Enzymatic Reactions)
- `load_model("SIR")`: Susceptible-Infected-Removed
- `load_model("doping_3way_oscillator")`: Three-way oscillator with doping reactions
- `load_model("repressilator")`: A repressilator model
# Getting Started
## Installation
1. Launch Julia's REPL (for example by entering julia in your shell)
2. Enter Pkg's REPL by typing `]`
3. Enter
```julia
pkg> add https://gitlab-research.centralesupelec.fr/2017bentrioum/markovprocesses.jl/
```
## Context - Mathematical framework
In this package, we are focused on Continuous-Time Markov Chains (CTMC, also called Markov jump processes), that can be described by Chemical Reaction Networks. The future state only depends on the current state. It is defined by two properties:
- ``\forall t, s \in \mathbb{R}_{\geq 0}, \mathbb{P}(S_t | S_v, v \in [0,s]) = \mathbb{P}(S_t| S_s)`` (Memoryless/Markov property)
- ``\forall t,v \in \mathbb{R}_{\geq 0}, t > v, \mathbb{P}(S_t | S_v) = \mathbb{P}(S_{t-v} | S_0)`` (Time-homogeneity).
A Chemical Reaction Network (CRN) is a formalism that describes biological phenomena. An example is the Susceptible-Infected-Removed model (SIR):
``
R_1: S + I \xrightarrow{k_i} 2I
``
``
R_2: I \xrightarrow{k_r} R
``
This CRN has two reactions that models two phenomena: infection ($R_1$) or immunisation ($R_2$). Each reaction is parametrized by a kinetic rate ($k_i$ or $k_r$). The stochastic dynamics of a CRN are described by CTMCs.
## Models
In the package, models are objects and their types all derived from the abstract type `Model`. Let's load the SIR model, which is pre-written within the package.
```julia
julia> load_model("SIR")
create_SIR (generic function with 1 method)
julia> println(SIR)
SIRModel <: ContinuousTimeModel model
- variables :
* I (index = 2 in state space)
* R (index = 3 in state space)
* S (index = 1 in state space)
- parameters :
* ki (index = 1 in parameter space)
* kr (index = 2 in parameter space)
- transitions : R1,R2
- observed variables :
* I (index = 1 in observed state space, index = 2 in state space)
p = [0.0012, 0.05]
x0 = [95, 5, 0]
t0 = 0.0
time bound = Inf
```
`load_model` has created a variable called SIR, which contains all the information for simulating the SIR model described above. It also contains a parameter vector and an initial point.
```julia
julia> @show SIR.p
@show SIR.x0
SIR.p = [0.0012, 0.05]
SIR.x0 = [95, 5, 0]
```
You can change the parameters or the initial state with the functions `set_param!` and `set_x0!`.
```julia
julia> set_param!(SIR, :ki, 0.015)
julia> @show SIR.p
SIR.p = [0.015, 0.05]
julia> set_param!(SIR, [0.02, 0.07])
julia> @show SIR.p
SIR.p = [0.02, 0.07]
julia> set_x0!(SIR, :S, 93)
julia> @show SIR.x0
SIR.x0 = [93, 5, 0]
```
## Trajectories
The simulation of the model is done by the function simulate.
```julia
julia> σ = simulate(SIR)
```
`simulate` returns a trajectory, which type is derived from `AbstractTrajectory`. It can be either an object of type `Trajectory` for models `::ContinuousTimeModel` or `SynchronizedTrajectory` for models that includes an automaton (but this is the subject of another section). It is easy to access the values of a trajectory.
```julia
julia> @show σ[3] # the third state of the trajectory
@show length_states(σ) # number of states
@show σ[:I, 4] # Fourth value of the variable I
@show σ.I[4] # Fourth value of the variable I
@show get_state_from_time(σ, 2.3)
σ[3] = [7]
length_states(σ) = 196
σ[:I, 4] = 8
σ.I[4] = 8
get_state_from_time(σ, 2.3) = [79]
```
The SIR object includes an observation model symbolized by the vector `SIR.g`. Even if the variables of the model are `[:S, :I, :R]`, only I will be observed.
```julia
julia> @show SIR.map_var_idx
@show SIR.g
@show size(σ.values), length(σ[:I]) # Only one column which corresponds to the I variable
SIR.map_var_idx = Dict(:I => 2,:R => 3,:S => 1)
SIR.g = [:I]
(size(σ.values), length(σ[:I])) = ((1,), 196)
```
The SIR model is by default unbounded, i.e. each trajectory is simulated until it reaches an absorbing state.
```julia
julia> @show isbounded(SIR)
@show isbounded(σ)
isbounded(SIR) = false
isbounded(σ) = false
```
We can bound the SIR's trajectories until time 120 by running `set_time_bound!(SIR, 120.0)`.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment