{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Create a model\n",
"\n",
"In this package, probabilistic models are represented by `ContinuousTimeModel` objects. There are several ways to create these kind of objects."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using MarkovProcesses"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. With load_model()\n",
"\n",
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Poisson process\n",
"load_model(\"poisson\")\n",
"poisson"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Repressilator model\n",
"load_model(\"repressilator\")\n",
"repressilator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. With @network_model\n",
"\n",
"A useful feature is avaiable for CTMCs induced by a Chemical Reaction Network. Let's consider the Chemical Reaction Network of the SIR model:\n",
"\n",
"$$\n",
"Infection: S + I \\xrightarrow{k_i} 2I \\\\\n",
"Recovery: I \\xrightarrow{k_r} R\n",
"$$\n",
"\n",
"The macro `@network_model` creates easily a CTMC stored in a `ContinuousTimeModel` variable based on this formalism."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"easy_sir = @network_model begin\n",
" Infection: (S + I => 2I, ki*I*S)\n",
" Recovery: (I => R, kr*I)\n",
"end \"My awesome SIR\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"set_param!(easy_sir, [0.0012, 0.05])\n",
"set_x0!(easy_sir, [95, 5, 0])\n",
"using Plots\n",
"σ = simulate(easy_sir)\n",
"plot(times(σ), σ.I, linetype = :steppost)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Manually (advanced)\n",
"\n",
"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:\n",
"\n",
"```julia\n",
"function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict{VariableModel,Int}, \n",
" map_param_idx::Dict{ParameterModel,Int}, transitions::Vector{<:Transition},\n",
" p::Vector{Float64}, x0::Vector{Int}, t0::Float64, \n",
" f!::Function, isabsorbing::Function; kwargs)\n",
"```\n",
"\n",
"First, one has to specify the dimensions of the state space and the parameter space."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dim_state_sir, dim_params_sir = 3, 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"map_var_idx_sir = Dict(:S => 1, :I => 2, :R => 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`map_var_params` is the equivalent of `map_var_idx` for parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"map_params_idx_sir = Dict(:ki => 1, :kr => 2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"transitions_sir = [:Infection, :Recovery]\n",
"p_sir = [0.0012, 0.05]\n",
"x0_sir = [95, 5, 0]\n",
"t0_sir = 0.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The two last arguments are functions, the first one, called `f!` must have the signature:\n",
"\n",
"```julia\n",
"function f!(xnplus1::Vector{Int}, ptr_t::Vector{Float64}, ptr_tr::Vector{Transition},\n",
" xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n",
"\n",
"```\n",
"\n",
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@everywhere function sir_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},\n",
" xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n",
" @inbounds a1 = p[1] * xn[1] * xn[2]\n",
" @inbounds a2 = p[2] * xn[2]\n",
" l_a = (a1, a2)\n",
" asum = sum(l_a)\n",
" if asum == 0.0\n",
" copyto!(xnplus1, xn)\n",
" return nothing\n",
" end\n",
" nu_1 = (-1, 1, 0)\n",
" nu_2 = (0, -1, 1)\n",
" l_nu = (nu_1, nu_2)\n",
" l_str_R = (:Infection, :Recovery)\n",
"\n",
" u1 = rand()\n",
" u2 = rand()\n",
" tau = - log(u1) / asum\n",
" b_inf = 0.0\n",
" b_sup = a1\n",
" reaction = 0\n",
" for i = 1:2 \n",
" if b_inf < asum*u2 < b_sup\n",
" reaction = i\n",
" break\n",
" end\n",
" @inbounds b_inf += l_a[i]\n",
" @inbounds b_sup += l_a[i+1]\n",
" end\n",
" \n",
" nu = l_nu[reaction]\n",
" for i = 1:3\n",
" @inbounds xnplus1[i] = xn[i]+nu[i]\n",
" end\n",
" @inbounds l_t[1] = tn + tau\n",
" @inbounds l_tr[1] = l_str_R[reaction]\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second function called `isaborbing` must have the signature:\n",
"\n",
"```julia\n",
"isabsorbing(p::Vector{Float64}, xn::Vector{Int})\n",
"```\n",
"\n",
"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.\n",
"\n",
"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",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@everywhere sir_isabsorbing(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally one sets the observed variables and the model can be created."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"g_sir = [:I]\n",
"\n",
"# Generates simulate method for the new model\n",
"@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:TryhardSIRModel))\n",
"@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:TryhardSIRModel))\n",
"@everywhere @eval $(MarkovProcesses.generate_code_simulation(:TryhardSIRModel, :sir_f!, :sir_isabsorbing))\n",
"\n",
"\n",
"tryhard_sir = TryhardSIRModel(dim_state_sir, dim_params_sir, \n",
" map_var_idx_sir, map_params_idx_sir, \n",
" transitions_sir, p_sir, x0_sir, t0_sir, \n",
" :sir_f!, :sir_isabsorbing; g = g_sir)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"σ = simulate(tryhard_sir)\n",
"plot(times(σ), σ.I, linetype = :steppost)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.3",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}