compartmental_models/SEIR Simulator with Parameter Estimation in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation & Parameter Estimation of SEIR Model\n",
    "\n",
    "SEIR is a type of compartmental models which are used in modelling of infectious disease using differential equations. These types of models divide the population into groups or compartments and the dynamics of these groups are expressed with the help of a system of differential equations.\n",
    "\n",
    "These system of equations are parametrized which capture the mechanistic nature of the disease. For simulation, you select values of these parameters and the resulting curves simulate the behaviour by solving the set of equations. Finally the results are plotted in a graph to visually understand the effect of the parameters.\n",
    "\n",
    "## SEIR Model\n",
    "\n",
    "For completeness the SEIR model is produced below:\n",
    "\n",
    "<img src=\"images/seir_model.png\">\n",
    "\n",
    "$\\displaystyle \\frac{dS}{dt} = -\\frac{\\beta S I}{N}$<br><br>\n",
    "$\\displaystyle \\frac{dE}{dt} = \\frac{\\beta S I}{N} - \\sigma E$<br><br>\n",
    "$\\displaystyle \\frac{dI}{dt} = \\sigma E -  \\gamma I$<br><br>\n",
    "$\\displaystyle \\frac{dR}{dt} = \\gamma I$<br><br>\n",
    "$N = S + E + I + R$<br><br>\n",
    "Where,<br><br>\n",
    "$\\beta$ is infection rate or the rate of spread<br><br>\n",
    "$\\sigma$ is the incubation rate or the rate of latent individuals becoming infectious (average duration of incubation is $1/\\sigma$)<br><br>\n",
    "$\\gamma$ is the recovery rate or mortality rate. If the duration of indection is D then $\\gamma$ = 1/D\n",
    "\n",
    "#### A supplementary blogpost to this notebook can be found at [Estimating Parameters of Compartmental Models from Observed Data](https://dilbertai.com/2020/04/27/estimating-parameters-of-compartmental-models-from-observed-data/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.integrate import odeint\n",
    "import plotly.graph_objects as go\n",
    "import plotly.io as pio\n",
    "import requests\n",
    "from lmfit import minimize, Parameters, Parameter, report_fit\n",
    "pio.renderers.default = \"notebook\"\n",
    "%matplotlib inline\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jupyter Specifics\n",
    "from IPython.display import HTML\n",
    "from ipywidgets.widgets import interact, IntSlider, FloatSlider, Layout, ToggleButton\n",
    "\n",
    "style = {'description_width': '100px'}\n",
    "slider_layout = Layout(width='99%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_model(z, t, beta, sigma, gamma):\n",
    "    \"\"\"\n",
    "    Reference https://www.idmod.org/docs/hiv/model-seir.html\n",
    "    \"\"\"\n",
    "    S, E, I, R = z\n",
    "    N = S + E + I + R\n",
    "    dSdt = -beta*S*I/N\n",
    "    dEdt = beta*S*I/N - sigma*E\n",
    "    dIdt = sigma*E - gamma*I\n",
    "    dRdt = gamma*I\n",
    "    return [dSdt, dEdt, dIdt, dRdt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_solver(t, initial_conditions, params):\n",
    "    initE, initI, initR, initN = initial_conditions\n",
    "    beta, sigma, gamma = params['beta'].value, params['sigma'].value, params['gamma'].value\n",
    "    initS = initN - (initE + initI + initR)\n",
    "    res = odeint(ode_model, [initS, initE, initI, initR], t, args=(beta, sigma, gamma))\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Request Success? True\n"
     ]
    }
   ],
   "source": [
    "response = requests.get('https://api.rootnet.in/covid19-in/stats/history')\n",
    "print('Request Success? {}'.format(response.status_code == 200))\n",
    "covid_history = response.json()['data']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "keys = ['day', 'total', 'confirmedCasesIndian', 'confirmedCasesForeign', 'confirmedButLocationUnidentified',\n",
    "        'discharged', 'deaths']\n",
    "df_covid_history = pd.DataFrame([[d.get('day'), \n",
    "                                  d['summary'].get('total'), \n",
    "                                  d['summary'].get('confirmedCasesIndian'), \n",
    "                                  d['summary'].get('confirmedCasesForeign'),\n",
    "                                  d['summary'].get('confirmedButLocationUnidentified'),\n",
    "                                  d['summary'].get('discharged'), \n",
    "                                  d['summary'].get('deaths')] \n",
    "                                 for d in covid_history],\n",
    "                    columns=keys)\n",
    "df_covid_history = df_covid_history.sort_values(by='day')\n",
    "df_covid_history['infected'] = df_covid_history['total'] - df_covid_history['discharged'] - df_covid_history['deaths']\n",
    "df_covid_history['total_recovered_or_dead'] = df_covid_history['discharged'] + df_covid_history['deaths']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ref: https://www.medrxiv.org/content/10.1101/2020.04.01.20049825v1.full.pdf\n",
    "initN = 1380000000\n",
    "# S0 = 966000000\n",
    "initE = 1000\n",
    "initI = 47\n",
    "initR = 0\n",
    "sigma = 1/5.2\n",
    "gamma = 1/2.9\n",
    "R0 = 4\n",
    "beta = R0 * gamma\n",
    "days = 112\n",
    "\n",
    "params = Parameters()\n",
    "params.add('beta', value=beta, min=0, max=10)\n",
    "params.add('sigma', value=sigma, min=0, max=10)\n",
    "params.add('gamma', value=gamma, min=0, max=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(initE, initI, initR, initN, beta, sigma, gamma, days, param_fitting):\n",
    "    initial_conditions = [initE, initI, initR, initN]\n",
    "    params['beta'].value, params['sigma'].value,params['gamma'].value = [beta, sigma, gamma]\n",
    "    tspan = np.arange(0, days, 1)\n",
    "    sol = ode_solver(tspan, initial_conditions, params)\n",
    "    S, E, I, R = sol[:, 0], sol[:, 1], sol[:, 2], sol[:, 3]\n",
    "    \n",
    "    # Create traces\n",
    "    fig = go.Figure()\n",
    "    if not param_fitting:\n",
    "        fig.add_trace(go.Scatter(x=tspan, y=S, mode='lines+markers', name='Susceptible'))\n",
    "        fig.add_trace(go.Scatter(x=tspan, y=E, mode='lines+markers', name='Exposed'))\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=I, mode='lines+markers', name='Infected'))\n",
    "    fig.add_trace(go.Scatter(x=tspan, y=R, mode='lines+markers',name='Recovered'))\n",
    "    if param_fitting:\n",
    "        fig.add_trace(go.Scatter(x=tspan, y=df_covid_history.infected, mode='lines+markers',\\\n",
    "                             name='Infections Observed', line = dict(dash='dash')))\n",
    "        fig.add_trace(go.Scatter(x=tspan, y=df_covid_history.total_recovered_or_dead, mode='lines+markers',\\\n",
    "                             name='Recovered/Deceased Observed', line = dict(dash='dash')))\n",
    "    \n",
    "    if days <= 30:\n",
    "        step = 1\n",
    "    elif days <= 90:\n",
    "        step = 7\n",
    "    else:\n",
    "        step = 30\n",
    "    \n",
    "    # Edit the layout\n",
    "    fig.update_layout(title='Simulation of SEIR Model',\n",
    "                       xaxis_title='Day',\n",
    "                       yaxis_title='Counts',\n",
    "                       title_x=0.5,\n",
    "                      width=900, height=600\n",
    "                     )\n",
    "    fig.update_xaxes(tickangle=-90, tickformat = None, tickmode='array', tickvals=np.arange(0, days + 1, step))\n",
    "    if not os.path.exists(\"images\"):\n",
    "        os.mkdir(\"images\")\n",
    "    fig.write_image(\"images/seir_simulation.png\")\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a787662501c544789daa0c0b3a4b3c80",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1000, description='initE', layout=Layout(width='99%'), max=100000, style…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(main, initE=IntSlider(min=0, max=100000, step=1, value=initE, description='initE', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               initI=IntSlider(min=0, max=100000, step=10, value=initI, description='initI', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               initR=IntSlider(min=0, max=100000, step=10, value=initR, description='initR', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               initN=IntSlider(min=0, max=1380000000, step=1000, value=initN, description='initN', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               beta=FloatSlider(min=0, max=4, step=0.01, value=beta, description='Infection rate', style=style, \\\n",
    "                                layout=slider_layout),\n",
    "               sigma=FloatSlider(min=0, max=4, step=0.01, value=sigma, description='Incubation rate', style=style, \\\n",
    "                                 layout=slider_layout),\n",
    "               gamma=FloatSlider(min=0, max=4, step=0.01, value=gamma, description='Recovery rate', style=style, \\\n",
    "                                 layout=slider_layout),\n",
    "               days=IntSlider(min=0, max=600, step=7, value=days, description='Days', style=style, \\\n",
    "                              layout=slider_layout),\n",
    "               param_fitting=ToggleButton(value=False, description='Fitting Mode', disabled=False, button_style='', \\\n",
    "             tooltip='Click to show fewer plots', icon='check-circle')\n",
    "        );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error(params, initial_conditions, tspan, data):\n",
    "    sol = ode_solver(tspan, initial_conditions, params)\n",
    "    return (sol[:, 2:4] - data).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_conditions = [initE, initI, initR, initN]\n",
    "beta = 1.08\n",
    "sigma = 0.02\n",
    "gamma = 0.02\n",
    "params['beta'].value = beta\n",
    "params['sigma'].value = sigma\n",
    "params['gamma'].value = gamma\n",
    "days = 45\n",
    "tspan = np.arange(0, days, 1)\n",
    "data = df_covid_history.loc[0:(days-1), ['infected', 'total_recovered_or_dead']].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><th> name </th><th> value </th><th> initial value </th><th> min </th><th> max </th><th> vary </th></tr><tr><td> beta </td><td>  1.08000000 </td><td> 1.3793103448275863 </td><td>  0.00000000 </td><td>  10.0000000 </td><td> True </td></tr><tr><td> sigma </td><td>  0.02000000 </td><td> 0.1923076923076923 </td><td>  0.00000000 </td><td>  10.0000000 </td><td> True </td></tr><tr><td> gamma </td><td>  0.02000000 </td><td> 0.3448275862068966 </td><td>  0.00000000 </td><td>  10.0000000 </td><td> True </td></tr></table>"
      ],
      "text/plain": [
       "Parameters([('beta', <Parameter 'beta', value=1.08, bounds=[0:10]>),\n",
       "            ('sigma', <Parameter 'sigma', value=0.02, bounds=[0:10]>),\n",
       "            ('gamma', <Parameter 'gamma', value=0.02, bounds=[0:10]>)])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit model and find predicted values\n",
    "result = minimize(error, params, args=(initial_conditions, tspan, data), method='leastsq')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><th> name </th><th> value </th><th> standard error </th><th> relative error </th><th> initial value </th><th> min </th><th> max </th><th> vary </th></tr><tr><td> beta </td><td>  0.25240521 </td><td>  0.02426587 </td><td> (9.61%) </td><td> 1.08 </td><td>  0.00000000 </td><td>  10.0000000 </td><td> True </td></tr><tr><td> sigma </td><td>  0.07915451 </td><td>  0.00865739 </td><td> (10.94%) </td><td> 0.02 </td><td>  0.00000000 </td><td>  10.0000000 </td><td> True </td></tr><tr><td> gamma </td><td>  0.02201161 </td><td>  0.00129203 </td><td> (5.87%) </td><td> 0.02 </td><td>  0.00000000 </td><td>  10.0000000 </td><td> True </td></tr></table>"
      ],
      "text/plain": [
       "Parameters([('beta',\n",
       "             <Parameter 'beta', value=0.25240520565167124 +/- 0.0243, bounds=[0:10]>),\n",
       "            ('sigma',\n",
       "             <Parameter 'sigma', value=0.07915450950986003 +/- 0.00866, bounds=[0:10]>),\n",
       "            ('gamma',\n",
       "             <Parameter 'gamma', value=0.022011605792732203 +/- 0.00129, bounds=[0:10]>)])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Fit Statistics]]\n",
      "    # fitting method   = leastsq\n",
      "    # function evals   = 64\n",
      "    # data points      = 90\n",
      "    # variables        = 3\n",
      "    chi-square         = 25843295.5\n",
      "    reduced chi-square = 297049.373\n",
      "    Akaike info crit   = 1137.09769\n",
      "    Bayesian info crit = 1144.59712\n",
      "[[Variables]]\n",
      "    beta:   0.25240521 +/- 0.02426587 (9.61%) (init = 1.08)\n",
      "    sigma:  0.07915451 +/- 0.00865739 (10.94%) (init = 0.02)\n",
      "    gamma:  0.02201161 +/- 0.00129203 (5.87%) (init = 0.02)\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(beta, sigma)  = -0.991\n",
      "    C(beta, gamma)  =  0.437\n",
      "    C(sigma, gamma) = -0.327\n"
     ]
    }
   ],
   "source": [
    "# display fitted statistics\n",
    "report_fit(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "        \n",
       "        \n",
       "            <div id=\"8dabc243-3547-47fc-b346-22f5a83e4c41\" class=\"plotly-graph-div\" style=\"height:600px; width:1000px;\"></div>\n",
       "            <script type=\"text/javascript\">\n",
       "                require([\"plotly\"], function(Plotly) {\n",
       "                    window.PLOTLYENV=window.PLOTLYENV || {};\n",
       "                    \n",
       "                if (document.getElementById(\"8dabc243-3547-47fc-b346-22f5a83e4c41\")) {\n",
       "                    Plotly.newPlot(\n",
       "                        '8dabc243-3547-47fc-b346-22f5a83e4c41',\n",
       "                        [{\"line\": {\"dash\": \"dot\"}, \"mode\": \"markers\", \"name\": \"Observed Infections\", \"type\": \"scatter\", \"x\": [0, 1, 2, 3, 4, 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], \"y\": [47, 60, 73, 70, 72, 95, 99, 120, 134, 149, 196, 256, 329, 424, 469, 553, 634, 640, 810, 901, 1118, 1238, 1649, 1860, 2322, 2784, 3219, 3851, 4312, 4714, 5218, 6039, 6634, 7409, 8048, 9272, 10197, 10824, 11616, 12289, 13295, 14255, 15122, 15859, 16689]}, {\"line\": {\"dash\": \"dot\"}, \"mode\": \"markers\", \"name\": \"Observed Recovered/Deceased\", \"type\": \"scatter\", \"x\": [0, 1, 2, 3, 4, 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], \"y\": [0, 0, 0, 12, 12, 15, 15, 17, 17, 24, 27, 27, 31, 43, 50, 53, 60, 84, 99, 123, 133, 159, 185, 209, 225, 288, 358, 430, 477, 560, 647, 722, 895, 1038, 1304, 1543, 1736, 1935, 2219, 2503, 2821, 3401, 3863, 4612, 5011]}, {\"mode\": \"lines+markers\", \"name\": \"Fitted Infections\", \"type\": \"scatter\", \"x\": [0, 1, 2, 3, 4, 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], \"y\": [47.0, 121.93945599912564, 191.78665420458944, 258.2448665323269, 322.7946114892799, 386.7408035519474, 451.25231885706546, 517.3953867147102, 586.1619882130099, 658.4942314388144, 735.3055166080162, 817.4991632130607, 905.9850715342029, 1001.6948955335492, 1105.5961330646087, 1218.7054833713573, 1342.1017733786975, 1476.9387208887902, 1624.4577682189138, 1786.001206670372, 1963.0257897031563, 2157.1170236610033, 2370.004320661652, 2603.577192918343, 2859.9026694397153, 3141.2441209246235, 3450.0816830673652, 3789.1344799591507, 4161.384860741957, 4570.1048757543, 5018.8852368285125, 5511.667026114684, 6052.776439046636, 6646.962872967847, 7299.440710038885, 8015.935120587657, 8802.732392127722, 9666.735075937857, 10615.522645805648, 11657.417941202446, 12801.560232392845, 14057.985398098288, 15437.713894108158, 16952.84751221901, 18616.675394705802]}, {\"mode\": \"lines+markers\", \"name\": \"Fitted Recovered/Deceased\", \"type\": \"scatter\", \"x\": [0, 1, 2, 3, 4, 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], \"y\": [0.0, 1.870368893590941, 5.330871494492761, 10.288621864948826, 16.68567978708348, 24.494682056501006, 33.71542489863358, 44.37224745082952, 56.512093004353055, 70.20314957284408, 85.53399010944248, 102.61314988787299, 121.56909110905595, 142.5505170287964, 165.7270075382864, 191.28995618656197, 219.4537960485481, 250.45750762764698, 284.566408430343, 322.07422757749725, 363.30547435651334, 408.61811365039114, 458.4065649870086, 513.1050462817893, 573.1912872214667, 639.1906411253515, 711.6806281883144, 791.295948272033, 878.7340039987671, 974.760981958488, 1080.218543424491, 1196.0311822480396, 1323.2143136053699, 1462.8831637593842, 1616.2625402634815, 1784.697556624641, 1969.6654304405097, 2172.78841975851, 2395.8480622648176, 2640.8007779352192, 2909.795033756946, 3205.190184086683, 3529.577145951644, 3885.8011437747773, 4276.986632913569]}],\n",
       "                        {\"height\": 600, \"template\": {\"data\": {\"bar\": [{\"error_x\": {\"color\": \"#2a3f5f\"}, \"error_y\": {\"color\": \"#2a3f5f\"}, \"marker\": {\"line\": {\"color\": \"#E5ECF6\", \"width\": 0.5}}, \"type\": \"bar\"}], \"barpolar\": [{\"marker\": {\"line\": {\"color\": \"#E5ECF6\", \"width\": 0.5}}, \"type\": \"barpolar\"}], \"carpet\": [{\"aaxis\": {\"endlinecolor\": \"#2a3f5f\", \"gridcolor\": \"white\", \"linecolor\": \"white\", \"minorgridcolor\": \"white\", \"startlinecolor\": \"#2a3f5f\"}, \"baxis\": {\"endlinecolor\": \"#2a3f5f\", \"gridcolor\": \"white\", \"linecolor\": \"white\", \"minorgridcolor\": \"white\", \"startlinecolor\": \"#2a3f5f\"}, \"type\": \"carpet\"}], \"choropleth\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"type\": \"choropleth\"}], \"contour\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"colorscale\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"type\": \"contour\"}], \"contourcarpet\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"type\": \"contourcarpet\"}], \"heatmap\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"colorscale\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"type\": \"heatmap\"}], \"heatmapgl\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"colorscale\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"type\": \"heatmapgl\"}], \"histogram\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"histogram\"}], \"histogram2d\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"colorscale\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"type\": \"histogram2d\"}], \"histogram2dcontour\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"colorscale\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"type\": \"histogram2dcontour\"}], \"mesh3d\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"type\": \"mesh3d\"}], \"parcoords\": [{\"line\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"parcoords\"}], \"pie\": [{\"automargin\": true, \"type\": \"pie\"}], \"scatter\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scatter\"}], \"scatter3d\": [{\"line\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scatter3d\"}], \"scattercarpet\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scattercarpet\"}], \"scattergeo\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scattergeo\"}], \"scattergl\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scattergl\"}], \"scattermapbox\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scattermapbox\"}], \"scatterpolar\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scatterpolar\"}], \"scatterpolargl\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scatterpolargl\"}], \"scatterternary\": [{\"marker\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"type\": \"scatterternary\"}], \"surface\": [{\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}, \"colorscale\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"type\": \"surface\"}], \"table\": [{\"cells\": {\"fill\": {\"color\": \"#EBF0F8\"}, \"line\": {\"color\": \"white\"}}, \"header\": {\"fill\": {\"color\": \"#C8D4E3\"}, \"line\": {\"color\": \"white\"}}, \"type\": \"table\"}]}, \"layout\": {\"annotationdefaults\": {\"arrowcolor\": \"#2a3f5f\", \"arrowhead\": 0, \"arrowwidth\": 1}, \"coloraxis\": {\"colorbar\": {\"outlinewidth\": 0, \"ticks\": \"\"}}, \"colorscale\": {\"diverging\": [[0, \"#8e0152\"], [0.1, \"#c51b7d\"], [0.2, \"#de77ae\"], [0.3, \"#f1b6da\"], [0.4, \"#fde0ef\"], [0.5, \"#f7f7f7\"], [0.6, \"#e6f5d0\"], [0.7, \"#b8e186\"], [0.8, \"#7fbc41\"], [0.9, \"#4d9221\"], [1, \"#276419\"]], \"sequential\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]], \"sequentialminus\": [[0.0, \"#0d0887\"], [0.1111111111111111, \"#46039f\"], [0.2222222222222222, \"#7201a8\"], [0.3333333333333333, \"#9c179e\"], [0.4444444444444444, \"#bd3786\"], [0.5555555555555556, \"#d8576b\"], [0.6666666666666666, \"#ed7953\"], [0.7777777777777778, \"#fb9f3a\"], [0.8888888888888888, \"#fdca26\"], [1.0, \"#f0f921\"]]}, \"colorway\": [\"#636efa\", \"#EF553B\", \"#00cc96\", \"#ab63fa\", \"#FFA15A\", \"#19d3f3\", \"#FF6692\", \"#B6E880\", \"#FF97FF\", \"#FECB52\"], \"font\": {\"color\": \"#2a3f5f\"}, \"geo\": {\"bgcolor\": \"white\", \"lakecolor\": \"white\", \"landcolor\": \"#E5ECF6\", \"showlakes\": true, \"showland\": true, \"subunitcolor\": \"white\"}, \"hoverlabel\": {\"align\": \"left\"}, \"hovermode\": \"closest\", \"mapbox\": {\"style\": \"light\"}, \"paper_bgcolor\": \"white\", \"plot_bgcolor\": \"#E5ECF6\", \"polar\": {\"angularaxis\": {\"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\"}, \"bgcolor\": \"#E5ECF6\", \"radialaxis\": {\"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\"}}, \"scene\": {\"xaxis\": {\"backgroundcolor\": \"#E5ECF6\", \"gridcolor\": \"white\", \"gridwidth\": 2, \"linecolor\": \"white\", \"showbackground\": true, \"ticks\": \"\", \"zerolinecolor\": \"white\"}, \"yaxis\": {\"backgroundcolor\": \"#E5ECF6\", \"gridcolor\": \"white\", \"gridwidth\": 2, \"linecolor\": \"white\", \"showbackground\": true, \"ticks\": \"\", \"zerolinecolor\": \"white\"}, \"zaxis\": {\"backgroundcolor\": \"#E5ECF6\", \"gridcolor\": \"white\", \"gridwidth\": 2, \"linecolor\": \"white\", \"showbackground\": true, \"ticks\": \"\", \"zerolinecolor\": \"white\"}}, \"shapedefaults\": {\"line\": {\"color\": \"#2a3f5f\"}}, \"ternary\": {\"aaxis\": {\"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\"}, \"baxis\": {\"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\"}, \"bgcolor\": \"#E5ECF6\", \"caxis\": {\"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\"}}, \"title\": {\"x\": 0.05}, \"xaxis\": {\"automargin\": true, \"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\", \"title\": {\"standoff\": 15}, \"zerolinecolor\": \"white\", \"zerolinewidth\": 2}, \"yaxis\": {\"automargin\": true, \"gridcolor\": \"white\", \"linecolor\": \"white\", \"ticks\": \"\", \"title\": {\"standoff\": 15}, \"zerolinecolor\": \"white\", \"zerolinewidth\": 2}}}, \"title\": {\"text\": \"Observed vs Fitted\", \"x\": 0.5}, \"width\": 1000, \"xaxis\": {\"title\": {\"text\": \"Day\"}}, \"yaxis\": {\"title\": {\"text\": \"Counts\"}}},\n",
       "                        {\"responsive\": true}\n",
       "                    ).then(function(){\n",
       "                            \n",
       "var gd = document.getElementById('8dabc243-3547-47fc-b346-22f5a83e4c41');\n",
       "var x = new MutationObserver(function (mutations, observer) {{\n",
       "        var display = window.getComputedStyle(gd).display;\n",
       "        if (!display || display === 'none') {{\n",
       "            console.log([gd, 'removed!']);\n",
       "            Plotly.purge(gd);\n",
       "            observer.disconnect();\n",
       "        }}\n",
       "}});\n",
       "\n",
       "// Listen for the removal of the full notebook cells\n",
       "var notebookContainer = gd.closest('#notebook-container');\n",
       "if (notebookContainer) {{\n",
       "    x.observe(notebookContainer, {childList: true});\n",
       "}}\n",
       "\n",
       "// Listen for the clearing of the current output cell\n",
       "var outputEl = gd.closest('.output');\n",
       "if (outputEl) {{\n",
       "    x.observe(outputEl, {childList: true});\n",
       "}}\n",
       "\n",
       "                        })\n",
       "                };\n",
       "                });\n",
       "            </script>\n",
       "        </div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "final = data + result.residual.reshape(data.shape)\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(x=tspan, y=data[:, 0], mode='markers', name='Observed Infections', line = dict(dash='dot')))\n",
    "fig.add_trace(go.Scatter(x=tspan, y=data[:, 1], mode='markers', name='Observed Recovered/Deceased', line = dict(dash='dot')))\n",
    "fig.add_trace(go.Scatter(x=tspan, y=final[:, 0], mode='lines+markers', name='Fitted Infections'))\n",
    "fig.add_trace(go.Scatter(x=tspan, y=final[:, 1], mode='lines+markers', name='Fitted Recovered/Deceased'))\n",
    "fig.update_layout(title='Observed vs Fitted',\n",
    "                       xaxis_title='Day',\n",
    "                       yaxis_title='Counts',\n",
    "                       title_x=0.5,\n",
    "                      width=1000, height=600\n",
    "                     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(48, 2)\n"
     ]
    }
   ],
   "source": [
    "observed_IR = df_covid_history.loc[:, ['infected', 'total_recovered_or_dead']].values\n",
    "print(observed_IR.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "tspan_fit_pred = np.arange(0, observed_IR.shape[0], 1)\n",
    "params['beta'].value = result.params['beta'].value\n",
    "params['sigma'].value = result.params['sigma'].value\n",
    "params['gamma'].value = result.params['gamma'].value\n",
    "fitted_predicted = ode_solver(tspan_fit_pred, initial_conditions, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(48, 2)\n"
     ]
    }
   ],
   "source": [
    "fitted_predicted_IR = fitted_predicted[:, 2:4]\n",
    "print(fitted_predicted_IR.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted MAE\n",
      "Infected:  587.1564520730407\n",
      "Recovered/Deceased:  220.23711453849532\n",
      "\n",
      "Fitted RMSE\n",
      "Infected:  702.6622222460597\n",
      "Recovered/Deceased:  283.8331485287845\n"
     ]
    }
   ],
   "source": [
    "print(\"Fitted MAE\")\n",
    "print('Infected: ', np.mean(np.abs(fitted_predicted_IR[:days, 0] - observed_IR[:days, 0])))\n",
    "print('Recovered/Deceased: ', np.mean(np.abs(fitted_predicted_IR[:days, 1] - observed_IR[:days, 1])))\n",
    "\n",
    "print(\"\\nFitted RMSE\")\n",
    "print('Infected: ', np.sqrt(np.mean((fitted_predicted_IR[:days, 0] - observed_IR[:days, 0])**2)))\n",
    "print('Recovered/Deceased: ', np.sqrt(np.mean((fitted_predicted_IR[:days, 1] - observed_IR[:days, 1])**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Predicted MAE\n",
      "Infected:  3500.8548519859637\n",
      "Recovered/Deceased:  894.9307531861135\n",
      "\n",
      "Predicted RMSE\n",
      "Infected:  3590.026376405217\n",
      "Recovered/Deceased:  901.1249714019484\n"
     ]
    }
   ],
   "source": [
    "print(\"Predicted MAE\")\n",
    "print('Infected: ', np.mean(np.abs(fitted_predicted_IR[days:observed_IR.shape[0], 0] - observed_IR[days:, 0])))\n",
    "print('Recovered/Deceased: ', np.mean(np.abs(fitted_predicted_IR[days:observed_IR.shape[0], 1] - observed_IR[days:, 1])))\n",
    "\n",
    "print(\"\\nPredicted RMSE\")\n",
    "print('Infected: ', np.sqrt(np.mean((fitted_predicted_IR[days:observed_IR.shape[0], 0] - observed_IR[days:, 0])**2)))\n",
    "print('Recovered/Deceased: ', np.sqrt(np.mean((fitted_predicted_IR[days:observed_IR.shape[0], 1] - observed_IR[days:, 1])**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6c4e20b45f724cec829aa0bba02c0ca2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=1000, description='initE', layout=Layout(width='99%'), max=100000, style…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact(main, initE=IntSlider(min=0, max=100000, step=1, value=initE, description='initE', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               initI=IntSlider(min=0, max=100000, step=10, value=initI, description='initI', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               initR=IntSlider(min=0, max=100000, step=10, value=initR, description='initR', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               initN=IntSlider(min=0, max=1380000000, step=1000, value=initN, description='initN', style=style, \\\n",
    "                               layout=slider_layout),\n",
    "               beta=FloatSlider(min=0, max=4, step=0.01, value=result.params['beta'].value, description='Infection rate', style=style, \\\n",
    "                                layout=slider_layout),\n",
    "               sigma=FloatSlider(min=0, max=4, step=0.01, value=result.params['sigma'].value, description='Incubation rate', style=style, \\\n",
    "                                 layout=slider_layout),\n",
    "               gamma=FloatSlider(min=0, max=4, step=0.01, value=result.params['gamma'].value, description='Recovery rate', style=style, \\\n",
    "                                 layout=slider_layout),\n",
    "               days=IntSlider(min=0, max=600, step=7, value=240, description='Days', style=style, \\\n",
    "                              layout=slider_layout),\n",
    "               param_fitting=ToggleButton(value=False, description='Fitting Mode', disabled=True, button_style='', \\\n",
    "             tooltip='Click to show fewer plots', icon='check-circle')\n",
    "        );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References:**<br>\n",
    "1. SEIR and SEIRS Model https://www.idmod.org/docs/hiv/model-seir.html<br>\n",
    "2. Compartmental models in epidemiology https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model<br>\n",
    "3. Solve Differential Equations in Python https://www.youtube.com/watch?v=VV3BnroVjZo<br>\n",
    "4. Computational Statistics in Python https://people.duke.edu/~ccc14/sta-663/CalibratingODEs.html<br>\n",
    "5. Ordinary Differential Equations (ODE) with Python and Jupyter https://elc.github.io/posts/ordinary-differential-equations-with-python/<br>\n",
    "6. SEIRS+ Model https://github.com/ryansmcgee/seirsplus<br>\n",
    "7. Stack Overflow https://stackoverflow.com/questions/40753159/why-is-scipy-minimize-ignoring-my-constraints<br>\n",
    "8. Lotka–Volterra equations https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations<br>\n",
    "9. SEIR and Regression Model based COVID-19 outbreak predictions in India https://www.medrxiv.org/content/10.1101/2020.04.01.20049825v1.full.pdf<br>\n",
    "\n",
    "A simulator built with RShiny which provides many more parameters https://alhill.shinyapps.io/COVID19seir/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}