Commit af372d32 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

new notebook examples: creation of a ContinuousTimeModel, ABC-SMC with ContinuousTimeModel

parent 2f893aa0
%% Cell type:markdown id: tags:
# 1. About simulation of models and trajectories
The goal of this notebook is to present the basics of the simulation in our package.
Let's get familiar with the package. First, we shoud load it.
%% Cell type:code id: tags:
``` julia
using MarkovProcesses
```
%%%% Output: stream
┌ Info: Precompiling MarkovProcesses [top-level]
└ @ Base loading.jl:1260
%% Cell type:markdown id: tags:
## Models
In this package, we are focused on Continuous-Time Markov Chains models that can be described by Chemical Reaction Networks.
Let's simulate our first model. We consider the famous SIR epidemiology model described by two reactions:
$$
Infection: S + I \xrightarrow{k_i} 2I \\
Recovery: I \xrightarrow{k_r} R
$$
In the MarkovProcesses package, models are objects that can be instantiatied and their types all derived from the abstract type `Model`. Let's load the SIR model.
In the MarkovProcesses package, models are objects and their types all derived from the abstract type `Model`. Let's load the SIR model.
%% Cell type:code id: tags:
``` julia
load_model("SIR")
```
%% Cell type:markdown id: tags:
This function searchs a model file called SIR.jl in the models/ directory (at the root of the package). This files creates the necessary resources to create a `ContinuousTimeModel` and store it in a variable called SIR.
This function searchs a model file called SIR.jl in the models/ directory (at the root of the package). This file creates the necessary resources to create a `ContinuousTimeModel` and store it in a variable called SIR.
%% Cell type:code id: tags:
``` julia
println(SIR)
```
%%%% Output: stream
SIR SSA pkg model (ContinuousTimeModel)
- 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.02, 0.07]
x0 = [95, 5, 0]
t0 = 0.0
time bound = Inf
%% Cell type:markdown id: tags:
This variable contains a parameter vector and an initial point. It is ready to use.
%% Cell type:code id: tags:
``` julia
@show SIR.p
@show SIR.x0
```
%%%% Output: stream
SIR.p = [0.0012, 0.05]
SIR.x0 = [95, 5, 0]
%%%% Output: execute_result
3-element Array{Int64,1}:
95
5
0
%% Cell type:markdown id: tags:
You can change the parameters with the function `set_param!`
%% Cell type:code id: tags:
``` julia
set_param!(SIR, :ki, 0.015)
@show SIR.p
set_param!(SIR, [0.02, 0.07])
@show SIR.p
```
%%%% Output: stream
SIR.p = [0.015, 0.05]
SIR.p = [0.02, 0.07]
%%%% Output: execute_result
2-element Array{Float64,1}:
0.02
0.07
%% Cell type:markdown id: tags:
## Trajectories
The simulation of the model is done by the function `simulate`.
%% Cell type:code id: tags:
``` julia
σ = simulate(SIR)
```
%%%% Output: execute_result
Trajectory
- Model name: SIR SSA pkg
- Variable trajectories:
* I: [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 46, 47, 48, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 79, 78, 79, 80, 81, 82, 83, 82, 81, 82, 83, 84, 85, 86, 87, 86, 87, 86, 87, 88, 87, 86, 85, 84, 85, 86, 87, 86, 85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
- times = [0.0, 0.24577855621552194, 0.2478149161565589, 0.2853150202535163, 0.4573465669300486, 0.5156802507122017, 0.5488713853151262, 0.5694106885211497, 0.5788872503885709, 0.6245912291868114, 0.6566121315415033, 0.7931283435353416, 0.83270020199368, 0.8386998075002131, 0.8417042952608857, 0.8427551376647552, 0.8701999554468036, 0.8951256292071661, 0.9077288218029982, 0.9376622319229073, 0.9827229649909012, 0.9953032576765457, 0.9986177826032556, 1.0153848291967649, 1.0447914188625567, 1.0492014122265003, 1.070524783103734, 1.0726269740951375, 1.072967138802752, 1.0943818500991553, 1.116489485084964, 1.1285735695980215, 1.1723375316129172, 1.213845977572097, 1.2286687178948017, 1.2295093483811073, 1.2312638191767724, 1.2325709282468362, 1.244481593691647, 1.2624870979965666, 1.2834145939861266, 1.285764796090055, 1.310972419324214, 1.324617274813682, 1.3451928372043, 1.390603165691919, 1.4116004938074227, 1.4505176830738031, 1.4547213579879592, 1.4624171730488276, 1.4632236357632131, 1.4702808868992607, 1.4742187030030112, 1.4754940369410445, 1.4755237349411963, 1.4797055276946836, 1.4822125736745166, 1.5234583191077604, 1.5472145358013738, 1.5643622279323768, 1.587883934525096, 1.6060487480766195, 1.6719097587058893, 1.692707314337466, 1.697088843328341, 1.7073061626597648, 1.710582536000933, 1.7131112181403199, 1.7241992997835416, 1.7634968816756837, 1.78550991157585, 1.8672533491313623, 1.868775294555033, 1.8780308579480804, 1.9271595351999966, 1.9433721804249033, 2.0115414023682696, 2.0307236004972524, 2.0331919769914935, 2.111600134727418, 2.129193253114897, 2.13941932331246, 2.1802017796212882, 2.203441732414025, 2.21058492861084, 2.2181127678095187, 2.2382137733157794, 2.269453289860572, 2.273242781280825, 2.3235086535004847, 2.354182479039521, 2.3687838068308014, 2.394257920516965, 2.4510845536452797, 2.5223982640674647, 2.5811289105807145, 2.6021582222564925, 2.7226222986870483, 2.7910197312036473, 2.9425350009980282, 3.02055077757811, 3.039594333770209, 3.133377977116126, 3.1817075376142987, 3.4147235863411423, 3.539782357009317, 3.587084899226989, 3.620447160024518, 3.6274441686043883, 3.8364757218525587, 3.8391673517442046, 3.915668148051359, 4.013608159662592, 4.153708574532225, 4.181307293785593, 4.383741529594886, 4.755130254057704, 5.4036065556424635, 5.410621801836842, 5.4330236196083534, 5.7776665099328985, 5.802258390336277, 5.953384202813328, 5.9957703036539165, 6.006366484766268, 6.033378984272586, 6.1441075669762135, 6.157705333270836, 6.212908365173471, 6.252648412582357, 6.392052855060355, 6.4025853260551, 6.518930517440784, 6.592596845941982, 6.7164258553402245, 7.132577546602882, 7.1635240737101045, 7.786486929175749, 7.801298892611052, 7.880452010019803, 8.023739641544415, 8.194073072136348, 8.227854219758548, 8.389105633271651, 8.997315201732246, 9.49806435329397, 9.684927150070815, 9.757433953477268, 10.384896097233423, 10.412553909276932, 10.482340218974503, 10.72193174545675, 10.729699137569275, 10.934608965950016, 11.11917441898455, 11.20446838592195, 11.446224333854014, 11.803183235069083, 12.468580993540257, 12.681442546818868, 13.245705889763494, 13.256101393479387, 14.089151139720176, 14.155401738255222, 14.177565343853095, 14.357063891191089, 14.894584527602456, 16.417373253382813, 16.43496115585706, 16.676589190572557, 17.44669462094279, 17.466179493638776, 18.621869021895897, 19.519031289839436, 20.33698130395561, 21.012241021922648, 22.94396170910104, 24.779124484090048, 25.904276533056056, 26.291353569941716, 26.757984872221343, 27.572970488225828, 28.090717842386944, 28.3914382857189, 30.188892511626793, 31.083610922646052, 31.275510714328856, 33.987557235729874, 36.0759352242401, 43.47952103086246, 43.64008186023434, 45.106699045418154, 46.59516722084342, 49.974791679941056, 65.89654404976704, 114.54035862296453]
- transitions = Union{Nothing, Symbol}[nothing, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R2, :R1, :R1, :R2, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R2, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R1, :R2, :R2, :R1, :R1, :R1, :R1, :R1, :R2, :R2, :R1, :R1, :R1, :R1, :R1, :R1, :R2, :R1, :R2, :R1, :R1, :R2, :R2, :R2, :R2, :R1, :R1, :R1, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2, :R2]
%% Cell type:markdown id: tags:
`simulate` returns a path 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 notebook). It is easy to access the values of a trajectory:
%% Cell type:code id: tags:
``` 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)
```
%%%% Output: stream
σ[3] = [7]
length_states(σ) = 196
σ[:I, 4] = 8
σ.I[4] = 8
get_state_from_time(σ, 2.3) = [83]
%%%% Output: execute_result
1-element Array{Int64,1}:
83
%% Cell type:markdown id: tags:
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.
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.
%% Cell type:code id: tags:
``` julia
@show SIR.map_var_idx
@show SIR.g
@show size(σ.values), length(σ[:I]) # Only one column which corresponds to the I variable
```
%%%% Output: stream
SIR.map_var_idx = Dict(:I => 2,:R => 3,:S => 1)
SIR.g = [:I]
(size(σ.values), length(σ[:I])) = ((1,), 196)
%%%% Output: execute_result
((1,), 196)
%% Cell type:markdown id: tags:
The SIR model is by default unbounded, i.e. each trajectory is simulated until it reaches an absorbing state.
%% Cell type:code id: tags:
``` julia
@show isbounded(SIR)
@show isbounded(σ)
```
%%%% Output: stream
isbounded(SIR) = false
isbounded(σ) = false
%%%% Output: execute_result
false
%% Cell type:markdown id: tags:
For example, computations of Lp integral distances needs bounded trajectories.
%% Cell type:code id: tags:
``` julia
dist_lp(simulate(SIR),simulate(SIR); p=2)
```
%%%% Output: stream
┌ Warning: Lp distance computed on unbounded trajectories. Result should be wrong
└ @ MarkovProcesses /home/moud/MarkovProcesses/markovprocesses.jl/core/trajectory.jl:61
%%%% Output: execute_result
31.26273873394255
%% Cell type:markdown id: tags:
But this feature can be changed.
%% Cell type:code id: tags:
``` julia
set_time_bound!(SIR, 120.0)
```
%%%% Output: execute_result
120.0
%% Cell type:code id: tags:
``` julia
@show dist_lp(simulate(SIR),simulate(SIR))
@show isbounded(SIR)
@show isbounded(simulate(SIR))
```
%%%% Output: stream
dist_lp(simulate(SIR), simulate(SIR)) = 164.33150971327115
isbounded(SIR) = true
isbounded(simulate(SIR)) = true
%%%% Output: execute_result
true
%% Cell type:markdown id: tags:
## Buffer size
A useful feature is that we can change the preallocation of the memory size during the simulation in order to gain computational performance.
%% Cell type:code id: tags:
``` julia
load_model("ER")
@show ER.buffer_size
@timev σ = simulate(ER)
@show length_states(σ)
```
%%%% Output: stream
ER.buffer_size = 10
0.036003 seconds (54.02 k allocations: 2.972 MiB)
elapsed time (ns): 36002718
bytes allocated: 3115887
pool allocs: 53998
non-pool GC allocs:18
length_states(σ) = 385
%%%% Output: execute_result
385
%% Cell type:markdown id: tags:
By default, the buffer size is 10. It means that a matrix of size 10 is allocated anf filled. When it is full, another matrix of size 10 is allocated. But if you know that your simulations will have a certain number of states, you can change the buffer size. It can make a big difference in terms of performance.
%% Cell type:code id: tags:
``` julia
ER.buffer_size = 100
@timev σ2 = simulate(ER)
@show length_states(σ2)
```
%%%% Output: stream
0.000154 seconds (411 allocations: 59.938 KiB)
elapsed time (ns): 153854
bytes allocated: 61376
pool allocs: 399
non-pool GC allocs:12
length_states(σ2) = 379
%%%% Output: execute_result
379
%% Cell type:markdown id: tags:
## Plots
%% Cell type:code id: tags:
``` julia
```
......
%% Cell type:markdown id: tags:
# Create a model
In this package, probabilistic models are represented by `ContinuousTimeModel` objects. There are several ways to create these kind of objects.
%% Cell type:code id: tags:
``` julia
using MarkovProcesses
```
%% Cell type:markdown id: tags:
## 1. With load_model()
A bunch of models are already available in 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`.
%% Cell type:code id: tags:
``` julia
# Poisson process
load_model("poisson")
poisson
```
%% Cell type:code id: tags:
``` julia
# Repressilator model
load_model("repressilator")
repressilator
```
%% Cell type:markdown id: tags:
## 2. With @network_model
A useful feature is avaiable for CTMCs induced by 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.
%% Cell type:code id: tags:
``` julia
easy_sir = @network_model begin
Infection: (S + I => 2I, ki*I*S)
Recovery: (I => R, kr*I)
end "My awesome SIR"
```
%% Cell type:markdown id: tags:
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.
%% Cell type:code id: tags:
``` julia
set_param!(easy_sir, [0.0012, 0.05])
set_x0!(easy_sir, [95, 5, 0])
using Plots
σ = simulate(easy_sir)
plot(times(σ), σ.I, linetype = :steppost)
```
%% Cell type:markdown id: tags:
## 3. Manually (advanced)
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)
```
First, one has to specify the dimensions of the state space and the parameter space.
%% Cell type:code id: tags:
``` julia
dim_state_sir, dim_params_sir = 3, 2
```
%% Cell type:markdown id: tags:
`map_var_idx` is a dictionary that maps each model variable (represented by a `Symbol`) to an index in the state space.
%% Cell type:code id: tags:
``` julia
map_var_idx_sir = Dict(:S => 1, :I => 2, :R => 3)
```
%% Cell type:markdown id: tags:
`map_var_params` is the equivalent of `map_var_idx` for parameters.
%% Cell type:code id: tags:
``` julia
map_params_idx_sir = Dict(:ki => 1, :kr => 2)
```
%% Cell type:markdown id: tags:
`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.
%% Cell type:code id: tags:
``` julia
transitions_sir = [:Infection, :Recovery]
p_sir = [0.0012, 0.05]
x0_sir = [95, 5, 0]
t0_sir = 0.0
```
%% Cell type:markdown id: tags:
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:
%% Cell type:code id: tags:
``` 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
```
%% Cell type:markdown id: tags:
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:
%% Cell type:code id: tags:
``` julia
@everywhere sir_isabsorbing(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
```
%% Cell type:markdown id: tags:
Finally one sets the observed variables and the model can be created.
%% Cell type:code id: tags:
``` julia
g_sir = [:I]
tryhard_sir = ContinuousTimeModel(dim_state_sir, dim_params_sir,
map_var_idx_sir, map_params_idx_sir,
transitions_sir, p_sir, x0_sir, t0_sir,
getfield(Main, :sir_f!), getfield(Main, :sir_isabsorbing); g = g_sir, name = "Handmade SIR")
```
%% Cell type:code id: tags:
``` julia
using Plots
σ = simulate(tryhard_sir)
plot(times(σ), σ.I, linetype = :steppost)
```
%% Cell type:code id: tags:
``` julia
```
%% Cell type:markdown id: tags:
# ABC-SMC with ContinuousTimeModel
The package allows the run of ABC-SMC algorithm with models defined by the package.
%% Cell type:code id: tags:
``` julia
using MarkovProcesses
using Distributions
load_model("SIR")
set_time_bound!(SIR, 100.0)
parametric_SIR = ParametricModel(SIR, (:ki, Uniform(0.001, 0.0015)))
vec_observations = [simulate(SIR)]
function dist_obs(vec_sim::Vector{Trajectory}, vec_observations::Vector{Trajectory})
return dist_lp(vec_sim[1], vec_observations[1]; p=2)
end
epsilon = 0.5 * dist_obs([simulate(SIR)], vec_observations)
@show epsilon
res_abc = abc_smc(parametric_SIR, vec_observations, dist_obs, nbr_particles = 100, tolerance = epsilon)
```
%% Cell type:code id: tags:
``` julia
using Plots
histogram(res_abc.mat_p_end, weights = res_abc.weights_end)
```
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment