slidernexttofigure_and_header_sir.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4a9ea413460942d7a26dae8d7793d550",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HTML(value='\\n<style>\\n.container {\\n  position: relative;\\n  text-align: center;\\n  color: white;\\n}\\n\\n.bott…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, clear_output\n",
    "from ipywidgets import interact, interactive,fixed, IntSlider, FloatSlider, HBox, Layout, Output, VBox, HTML,HTMLMath,Box,Text,link,FloatText\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "from numpy import *\n",
    "\n",
    "HTML(\"\"\"\n",
    "<style>\n",
    ".container {\n",
    "  position: relative;\n",
    "  text-align: center;\n",
    "  color: white;\n",
    "}\n",
    "\n",
    ".bottom-left {\n",
    "  position: absolute;\n",
    "  bottom: 8px;\n",
    "  left: 16px;\n",
    "}\n",
    "\n",
    ".top-left {\n",
    "  position: absolute;\n",
    "  top: 8px;\n",
    "  left: 16px;\n",
    "}\n",
    "\n",
    ".top-right {\n",
    "  position: absolute;\n",
    "  top: 8px;\n",
    "  right: 16px;\n",
    "}\n",
    "\n",
    ".bottom-right {\n",
    "  position: absolute;\n",
    "  bottom: 8px;\n",
    "  right: 16px;\n",
    "}\n",
    "\n",
    ".centered {\n",
    "  position: absolute;\n",
    "  top: 50%;\n",
    "  left: 50%;\n",
    "  transform: translate(-50%, -50%);\n",
    "}\n",
    "\n",
    "</style>\n",
    "    \"\"\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a1601298bc4843b28d7a314d500b12ca",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value='\\n<div class=\"container\">\\n<img src=\"https://raw.githubusercontent.com/bigfooted/co…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "out = Output(layout={'width': '60%'})\n",
    "\n",
    "N = 1000000  # initial population size\n",
    "I = 10       # initial infected\n",
    "R = 0        # initial removed\n",
    "\n",
    "def update_plot(a,b): \n",
    "    \n",
    "    with out: \n",
    "        clear_output(wait=True) # prevents distortions\n",
    "        \n",
    "        fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2, figsize=(16, 8))\n",
    "\n",
    "        def dX_dt(X, t=0):\n",
    "            \"\"\" Return the susceptible and infected population. \"\"\"\n",
    "            return array([ -   a*X[0]*X[1]/N ,              # susceptible\n",
    "                               a*X[0]*X[1]/N - b*X[1],      # infected\n",
    "                                               b*X[1]])     # removed\n",
    "        #t = linspace(0, 90,  180)              # time\n",
    "        t=np.arange(0,90,0.5)\n",
    "        X0 = array([N, I,R])                     # initials conditions: 10 susceptible and 0 infected\n",
    "        X, infodict = integrate.odeint(dX_dt, X0, t, full_output=True)\n",
    "        susceptible, infected,removed = X.T        \n",
    "\n",
    "        scat, = ax1.plot(t, susceptible,'r-',label='susceptible',linewidth=4)\n",
    "        scat, = ax1.plot(t, infected,'b-',label='infected',linewidth=4)\n",
    "        scat, = ax1.plot(t, removed,'k-',label='removed',linewidth=4)\n",
    "        ax1.legend(loc='best')\n",
    "        ax1.grid()\n",
    "        ax1.set_title('susceptible and infected part of the population',fontsize=18)\n",
    "        ax1.set_xlabel('time [days]', fontsize=15)\n",
    "        ax1.set_ylabel('population', fontsize=15)\n",
    "        ax1.set_title('susceptible and infected',fontsize=18)\n",
    "\n",
    "        # we are on this line:\n",
    "        # V = delta*x -gamma*ln(x) + beta*y - alpha*ln(y)\n",
    "        #our_V = c*X0[0] -d*log(X0[0]) + b*X0[1] - a*log(X0[1])\n",
    "        # ---------------------------------------------------------\n",
    "        values  = linspace(0.33, 2, 6)                          # position of the orbits around our initial conditions\n",
    "        \n",
    "                \n",
    "        vcolors = plt.cm.autumn_r(linspace(0.3, 1., len(values)))  # colors for each trajectory\n",
    "        #X_f0 = array([     0. ,  0.])\n",
    "        #X_f1 = array([ d/c, a/b])  # stable point (equilibrium)\n",
    "        #equil = \"equilibrium point = (\" + str(X_f1[0]) + \",\" + str(X_f1[1]) + \")\"\n",
    "    \n",
    "        #------------------------------------------------------------------\n",
    "        # plot trajectories\n",
    "        for v, col in zip(values, vcolors):\n",
    "        \n",
    "            #X0 = v * X_f1                               # starting point\n",
    "            X_0 = v*X0                               # starting point\n",
    "            X = integrate.odeint( dX_dt, X_0, t)         # we don't need infodict here\n",
    "            ax2.plot( X[:,0], X[:,1], lw=5.0, color=col, label='X0=(%.f, %.f)' % ( X_0[0], X_0[1]) )\n",
    "\n",
    "        # location of the equilibrium point     \n",
    "        #ax2.plot(X_f1[0],X_f1[1],'o')    \n",
    "        \n",
    "        #-------------------------------------------------------\n",
    "        # define a grid and compute direction at each point\n",
    "        ymax = plt.ylim(ymin=0)[1]                        # get axis limits\n",
    "        xmax = plt.xlim(xmin=0)[1]\n",
    "        nb_points   = 20\n",
    "\n",
    "        x = linspace(0, xmax, nb_points)\n",
    "        y = linspace(0, ymax, nb_points)\n",
    "\n",
    "        X1 , Y1  = meshgrid(x, y)                       # create a grid\n",
    "        DX1, DY1, DZ1 = dX_dt([X1, Y1])                 # compute growth rate on the grid (for S,I,R)\n",
    "        M = (hypot(DX1, DY1))                           # Norm of the growth rate \n",
    "        M[ M == 0] = 1.                                 # Avoid zero division errors \n",
    "        DX1 /= M                                        # Normalize each arrows\n",
    "        DY1 /= M\n",
    "\n",
    "        ax2.set_title('Trajectories and direction fields',fontsize=18)\n",
    "        ax2.set_xlabel('Number of Susceptible', fontsize=15)\n",
    "        ax2.set_ylabel('Number of Infected',fontsize=15)\n",
    "        Q=ax2.quiver(X1, Y1, DX1, DY1, M, pivot='mid', cmap=plt.cm.jet)\n",
    "    \n",
    "\n",
    "\n",
    "        \n",
    "        #fig.tight_layout()\n",
    "        plt.show()\n",
    "        fig.savefig('susceptible_and_infected.png')\n",
    "\n",
    "        \n",
    "        \n",
    "headertextbox = HTML(value=\"\"\"\n",
    "<div class=\"container\">\n",
    "<img src=\"https://raw.githubusercontent.com/bigfooted/covid19-seir/master/coronavirus_cropped.jpg\" alt=\"corona\" style=\"width:100%;margin:-8px 0px\">\n",
    "   <div class=\"centered\">\n",
    "          <font size=\"10\">\n",
    "          SIR model of a virus outbreak\n",
    "              </font> \n",
    "     </div>\n",
    "</div>\n",
    "    \"\"\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#textbox_R0 = Text(\n",
    "#value=R0text\n",
    "#    ,layout=Layout(width='200px')\n",
    "#)\n",
    "\n",
    "\n",
    "\n",
    "textbox = HTMLMath(\n",
    "value=\"\"\"<body><font size=\"4\"> Classic SIR model: Kermack-McKendrick: <br> \\\n",
    "        $$\\\\dot S = -\\\\beta I S/N$$<br> \\\n",
    "        $$\\\\dot I = +\\\\beta I S/N - \\\\gamma I $$ <br> \\\n",
    "        $$\\\\dot R = \\\\gamma I$$ <br> \\\n",
    "        $S(t)$ : number of Susceptible persons (persons that are not infected or immune) <br> \\\n",
    "        $I(t)$ : number of Infected persons <br> \\\n",
    "        $R(t)$ : number of Removed persons (immune or dead) <br> \\\n",
    "        $N = S+I+R$ : total population <br> \\\n",
    "        $\\\\beta=1/T_l$ : infection rate = inverse of latent period <br> \\\n",
    "        $\\\\gamma=1/T_r$ : recovery rate = inverse of recovery time.  <br> \\\n",
    "        A virus can spread when the number of infections increases: $\\\\dot I>0$. This means that S=N, <br><br>\\\n",
    "        $$R_0 = \\\\beta / \\\\gamma > 1$$ \\\n",
    "        </font></body> \n",
    "       \"\"\"\n",
    "    ,layout=Layout(width='30%')\n",
    ")\n",
    "\n",
    "\n",
    "# to do: an output box displaying the location of the equilibrium plot\n",
    "\n",
    "\n",
    "layout=Layout(border='0px solid black', width='400px',height='50px')\n",
    "\n",
    "style = {'description_width': '150px','width':'500px'}\n",
    "sliderInfectionrate = FloatSlider(min=0.2, max=2.0, step=0.1, value=0.4, description=\"infection rate:\",orientation='horizontal',style=style,layout=layout) \n",
    "sliderRecoveryrate = FloatSlider(min=0.02, max=0.5, step=0.01, value=0.2, description=\"recovery rate:\",orientation='horizontal',style=style,layout=layout) \n",
    "#l = link((sliderInfectionrate, 'value'), (textbox_R0, 'value'))\n",
    "\n",
    "r0 = Output(\n",
    "    layout={'width': '90%'}\n",
    ")\n",
    "def update_R0(a,b):\n",
    "    with r0:\n",
    "        clear_output(wait=True)\n",
    "        # update the text\n",
    "        print(\"      Basic reproduction number R0 = {a:2.2f}\".format(a=(sliderInfectionrate.value/sliderRecoveryrate.value)))\n",
    "        \n",
    "\n",
    "\n",
    "sliders = interactive(update_plot,a=sliderInfectionrate,b=sliderRecoveryrate)\n",
    "textr0 = interactive(update_R0,a=sliderInfectionrate,b=sliderRecoveryrate)\n",
    "\n",
    "update_plot(a=0.4,b=0.2) # to show the initial plot\n",
    "update_R0(a=0.4,b=0.2)   # to show the initial R0\n",
    "display(VBox([headertextbox,HBox([out, VBox([sliderInfectionrate,sliderRecoveryrate, r0],layout=Layout(width='30%')), textbox])]))\n"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}