Skip to content
Snippets Groups Projects
Commit 66d44c7e authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

update of notebooks according to the new simulation methods

parent 920cd91a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
core/MarkovProcesses.jl
+ 1
− 1
View file @ 66d44c7e
......@@ -10,7 +10,7 @@ import Distributed: @everywhere, @distributed
import Distributions: Product, Uniform, Normal
import Distributions: Distribution, Univariate, Continuous, UnivariateDistribution,
MultivariateDistribution, product_distribution
import Distributions: insupport, pdf
import Distributions: insupport, isbounded, pdf
import FunctionWrappers: FunctionWrapper
import Random: rand, rand!
import StaticArrays: SVector, @SVector
......
examples/notebooks/1_Simulation.ipynb
+ 3
− 3
View file @ 66d44c7e
......@@ -271,15 +271,15 @@
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.4.2",
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.4"
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.2"
"version": "1.5.3"
}
},
"nbformat": 4,
......
examples/notebooks/2_Create_a_model.ipynb
+ 9
− 2
View file @ 66d44c7e
%% 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,
# 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,
getfield(Main, :sir_f!), getfield(Main, :sir_isabsorbing); g = g_sir, name = "Handmade SIR")
:sir_f!, :sir_isabsorbing; g = g_sir)
```
%% Cell type:code id: tags:
``` julia
using Plots
σ = simulate(tryhard_sir)
plot(times(σ), σ.I, linetype = :steppost)
```
%% Cell type:code id: tags:
``` julia
```
......
examples/notebooks/3_ABC-SMC_with_ContinuousTimeModel.ipynb
+ 1
− 0
View file @ 66d44c7e
%% 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(vec(res_abc.mat_p_end), weights = res_abc.weights)
```
%% Cell type:code id: tags:
``` julia
```
......
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