{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. About simulation of models and trajectories\n",
    "\n",
    "The goal of this notebook is to present the basics of the simulation in our package.\n",
    "\n",
    "Let's get familiar with the package. First, we shoud load it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using MarkovProcesses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models\n",
    "\n",
    "In this package, we are focused on Continuous-Time Markov Chains models that can be described by Chemical Reaction Networks.\n",
    "\n",
    "Let's simulate our first model. We consider the famous SIR epidemiology model described by two reactions:\n",
    "\n",
    "$$\n",
    "Infection: S + I \\xrightarrow{k_i} 2I \\\\\n",
    "Recovery: I \\xrightarrow{k_r} R\n",
    "$$\n",
    "\n",
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_model(\"SIR\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(SIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This variable contains a parameter vector and an initial point. It is ready to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@show SIR.p\n",
    "@show SIR.x0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can change the parameters with the function `set_param!`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_param!(SIR, :ki, 0.015)\n",
    "@show SIR.p\n",
    "set_param!(SIR, [0.02, 0.07])\n",
    "@show SIR.p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trajectories\n",
    "\n",
    "The simulation of the model is done by the function `simulate`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ = simulate(SIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@show σ[3] # the third state of the trajectory\n",
    "@show length_states(σ) # number of states\n",
    "@show σ[:I, 4] # Fourth value of the variable I\n",
    "@show σ.I[4] # Fourth value of the variable I\n",
    "@show get_state_from_time(σ, 2.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@show SIR.map_var_idx\n",
    "@show SIR.g\n",
    "@show size(σ.values), length(σ[:I]) # Only one column which corresponds to the I variable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SIR model is by default unbounded, i.e. each trajectory is simulated until it reaches an absorbing state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@show isbounded(SIR)\n",
    "@show isbounded(σ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, computations of Lp integral distances needs bounded trajectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_lp(simulate(SIR),simulate(SIR); p=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this feature can be changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_time_bound!(SIR, 120.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@show dist_lp(simulate(SIR),simulate(SIR))\n",
    "@show isbounded(SIR)\n",
    "@show isbounded(simulate(SIR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Buffer size\n",
    "\n",
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_model(\"ER\")\n",
    "@show ER.buffer_size\n",
    "@timev σ = simulate(ER)\n",
    "@show length_states(σ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ER.buffer_size = 100\n",
    "@timev σ2 = simulate(ER)\n",
    "@show length_states(σ2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plots"
   ]
  },
  {
   "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
}