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

fix of syntaxes in first notebook

parent 8fe47a29
%% 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.
%% 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.
%% Cell type:code id: tags:
``` julia
SIR
println(SIR)
```
%%%% Output: execute_result
%%%% Output: stream
ContinuousTimeModel(3, 2, Dict("S" => 1,"I" => 2,"R" => 3), Dict("I" => 1), Dict("kr" => 2,"ki" => 1), Union{Nothing, String}["R1", "R2"], [0.0012, 0.05], [95, 5, 0], 0.0, MarkovProcesses.SIR_f!, ["I"], [2], MarkovProcesses.isabsorbing_SIR, Inf, 10)
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)
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(ContinuousTimeModel(3, 2, Dict("S" => 1,"I" => 2,"R" => 3), Dict("I" => 1), Dict("kr" => 2,"ki" => 1), Union{Nothing, String}["R1", "R2"], [0.02, 0.07], [95, 5, 0], 0.0, MarkovProcesses.SIR_f!, ["I"], [2], MarkovProcesses.isabsorbing_SIR, Inf, 10), [5; 6; … ; 1; 0], [0.0, 0.0002090421844976379, 0.026738149867351305, 0.1546747797181483, 0.1670337962687422, 0.1864041143902683, 0.22499455696452758, 0.2609816573473689, 0.26211457357364876, 0.34981665245812765 … 28.156954777152635, 29.8789510639123, 31.349526365163364, 31.418132567003155, 43.3115664747499, 49.84826982651157, 50.582546017297425, 51.33856085518521, 55.25637949865968, 63.081159577787396], Union{Nothing, String}[nothing, "R1", "R1", "R1", "R1", "R1", "R1", "R1", "R1", "R1" … "R2", "R2", "R2", "R2", "R2", "R2", "R2", "R2", "R2", "R2"])
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 σ.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
get_state_from_time(σ, 2.3) = [76]
σ[:I, 4] = 8
σ.I[4] = 8
get_state_from_time(σ, 2.3) = [83]
%%%% Output: execute_result
1-element view(::Array{Int64,2}, 84, :) with eltype Int64:
76
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.
%% 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
@show size(σ.values), length(σ[:I]) # Only one column which corresponds to the I variable
```
%%%% Output: stream
SIR.map_var_idx = Dict("S" => 1,"I" => 2,"R" => 3)
SIR.g = ["I"]
(size(σ.values), length(σ["I"])) = ((196, 1), 196)
SIR.map_var_idx = Dict(:I => 2,:R => 3,:S => 1)
SIR.g = [:I]
(size(σ.values), length(σ[:I])) = ((1,), 196)
%%%% Output: execute_result
((196, 1), 196)
((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
20.942860320770663
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)) = 270.22871331339087
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.070949 seconds (130.66 k allocations: 7.141 MiB)
elapsed time (ns): 70948721
bytes allocated: 7488324
pool allocs: 130585
non-pool GC allocs:73
length_states(σ) = 407
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
407
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.000232 seconds (3.02 k allocations: 212.000 KiB)
elapsed time (ns): 231641
bytes allocated: 217088
pool allocs: 3004
non-pool GC allocs:11
length_states(σ2) = 425
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
425
379
%% Cell type:markdown id: tags:
## Plots
%% Cell type:code id: tags:
``` julia
```
......
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