epi_uni_app.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Epidemic at a university\n",
    "\n",
    "This app will simulate an epidemic on a university campus (or more broadly, a fairly small, 'closed' population). Key assumptions of the model are:\n",
    "\n",
    "* There is a simple Susceptible -> Exposed -> Infectious -> Recovered framework\n",
    "* The exposed and infectious period both have mean length of 1 week, and immunity is long-lasting.\n",
    "* 1000 people live in households of 10. Most mixing is within the household, with a small amount of random mixing.\n",
    "* Everyone attends a daily 2-hour class with a fixed group of 100 students.\n",
    "\n",
    "The app will show the outcome in the default case when no interventions are applied. You can use the sliders below to choose a mix of different interventions:\n",
    "\n",
    "* **Testing** - There is weekly testing of a percentage of the population with immediate results. The tests are 95% accurate at finding cases (exposed and infectious) and have a 1% false positive rate for susceptible and recovered individuals.\n",
    "* **Isolation** - Positive cases and those contact traced (below) fully isolate for a chosen number of days.\n",
    "* **Tracing** - A chosen (maximum) number of close contacts are identified from the household and class of identified cases to quarantine for the isolation period.\n",
    "* **Mixing reduction** - Restrictions on mixing of individuals outside of households and classes can be enacted.\n",
    "* **Class size** - The size of the daily classes can be lowered (creating more classes).\n",
    "\n",
    "Once you have chosen your intervention, click 'Run Interact'. The model will run 10 times (each run can take 10-20 seconds) and then plot the output both with the default assumptions and your chosen interventions. It will also tell you the peak and total number infected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e799679683d94f9e8d7c58ef2fc7b7b4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>IntSlider</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "IntSlider(value=0, description='% Testing', step=5, style=SliderStyle(description_width='initial'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d7664fd8842e49658becf3af70acf90f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>IntSlider</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "IntSlider(value=1, description='Days isolating', max=10, min=1, style=SliderStyle(description_width='initial'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c485a40a517345adb5c31bdfb8a3460e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>IntSlider</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "IntSlider(value=0, description='Max no. contacts', max=10, style=SliderStyle(description_width='initial'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b689cdc338064beaa779472a9b02c860",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>IntSlider</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "IntSlider(value=0, description='% Mixing reduction', step=10, style=SliderStyle(description_width='initial'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fedadc630df44a88ac37055d30eb2e09",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>IntSlider</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "IntSlider(value=100, description='Class size', min=10, step=10, style=SliderStyle(description_width='initial'))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56169e7031f945acaf4be31e8b3ccefb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(Button(description='Run Interact', style=ButtonStyle()), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.runmodel>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### Import the necessary libraries\n",
    "#!jupyter nbextension enable --py widgetsnbextension --sys-prefix\n",
    "#!jupyter serverextension enable voila --sys-prefix\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact_manual\n",
    "\n",
    "style = {'description_width': 'initial'}\n",
    "test_w=widgets.IntSlider(min=0, max=100, step=5, value=0, description='% Testing',style=style)\n",
    "isolate_w=widgets.IntSlider(min=1, max=10, step=1, value=0, description='Days isolating',style=style)\n",
    "trace_w=widgets.IntSlider(min=0, max=10, step=1, value=0, description='Max no. contacts',style=style)\n",
    "class_w=widgets.IntSlider(min=10, max=100, step=10, value=100, description='Class size',style=style)\n",
    "mix_w=widgets.IntSlider(min=0, max=100, step=10, value=0, description='% Mixing reduction',style=style)\n",
    "\n",
    "### Global constants\n",
    "def runmodel():\n",
    "    \n",
    "    POPSIZE=1000 # Total individuals in population\n",
    "    NOHOUSES=100 # No. of households (mean household size = POPSIZE / NOHOUSES)\n",
    "\n",
    "    OMEGA=1/7 # Latent rate\n",
    "    GAMMA=1/7 # Recovery rate\n",
    "    WANE=0 # Waning immunity rate\n",
    "    ISOLATE=5\n",
    "    BETA=np.zeros(2)\n",
    "    BETA[0]=0.02 # Transmission coefficient at home\n",
    "    BETA[1]=0.02 # Transmission coefficient at class\n",
    "    #BETASS=0.02\n",
    "\n",
    "    REPS=10 # No. of replicate runs\n",
    "    I0=10 # Initial no. of infected individuals (mean)\n",
    "    MAX_TIME=150 # Max time\n",
    "    CLASS_START = 0.4 # Time in day when move to class \n",
    "    CLASS_END = 0.48333 # Time in day when move to house\n",
    "\n",
    "    statstore=np.zeros(4)\n",
    "\n",
    "    plt.rcParams.update({'font.size': 16})\n",
    "    \n",
    "    solstore=np.zeros((10,150))\n",
    "    solstore[0]=[0.0, 11.0, 11.0, 10.0, 9.0, 10.0, 8.0, 9.0, 13.0, 17.0, 16.0, 17.0, 23.0, 28.0, 29.0, 27.0, 31.0, 34.0, 36.0, 39.0, 39.0, 49.0, 42.0, 42.0, 47.0, 51.0, 59.0, 61.0, 69.0, 70.0, 73.0, 75.0, 87.0, 93.0, 110.0, 125.0, 133.0, 141.0, 151.0, 157.0, 166.0, 170.0, 187.0, 201.0, 194.0, 197.0, 205.0, 204.0, 194.0, 199.0, 193.0, 181.0, 175.0, 170.0, 176.0, 162.0, 151.0, 146.0, 139.0, 131.0, 123.0, 108.0, 99.0, 92.0, 85.0, 81.0, 78.0, 74.0, 74.0, 68.0, 50.0, 48.0, 42.0, 41.0, 37.0, 38.0, 34.0, 29.0, 32.0, 31.0, 28.0, 25.0, 23.0, 21.0, 20.0, 20.0, 18.0, 15.0, 12.0, 11.0, 9.0, 10.0, 9.0, 7.0, 7.0, 7.0, 6.0, 5.0, 5.0, 5.0, 5.0, 5.0, 4.0, 4.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 5.0, 4.0, 4.0, 4.0, 4.0, 4.0, 3.0, 3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[1]=[0.0, 11.0, 12.0, 15.0, 13.0, 12.0, 15.0, 14.0, 15.0, 15.0, 18.0, 17.0, 24.0, 24.0, 28.0, 31.0, 33.0, 39.0, 40.0, 39.0, 42.0, 43.0, 50.0, 56.0, 59.0, 71.0, 88.0, 94.0, 100.0, 110.0, 113.0, 125.0, 122.0, 133.0, 143.0, 152.0, 151.0, 154.0, 157.0, 157.0, 151.0, 161.0, 150.0, 159.0, 162.0, 158.0, 160.0, 160.0, 147.0, 140.0, 139.0, 135.0, 131.0, 125.0, 119.0, 115.0, 110.0, 107.0, 102.0, 93.0, 88.0, 84.0, 74.0, 76.0, 74.0, 68.0, 63.0, 59.0, 50.0, 50.0, 45.0, 42.0, 36.0, 37.0, 37.0, 35.0, 34.0, 33.0, 27.0, 24.0, 24.0, 18.0, 15.0, 12.0, 10.0, 9.0, 8.0, 8.0, 7.0, 6.0, 5.0, 4.0, 4.0, 3.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[2]=[0.0, 7.0, 6.0, 6.0, 8.0, 12.0, 12.0, 12.0, 14.0, 16.0, 17.0, 19.0, 20.0, 23.0, 20.0, 25.0, 26.0, 35.0, 44.0, 46.0, 45.0, 45.0, 56.0, 58.0, 62.0, 71.0, 71.0, 80.0, 88.0, 99.0, 104.0, 114.0, 113.0, 122.0, 128.0, 146.0, 154.0, 167.0, 173.0, 172.0, 163.0, 181.0, 173.0, 174.0, 175.0, 170.0, 174.0, 171.0, 165.0, 161.0, 158.0, 148.0, 140.0, 127.0, 125.0, 125.0, 130.0, 118.0, 113.0, 103.0, 99.0, 91.0, 84.0, 81.0, 75.0, 70.0, 66.0, 62.0, 60.0, 59.0, 53.0, 46.0, 45.0, 42.0, 42.0, 39.0, 34.0, 32.0, 26.0, 23.0, 23.0, 24.0, 20.0, 19.0, 17.0, 15.0, 14.0, 12.0, 9.0, 8.0, 6.0, 4.0, 4.0, 6.0, 6.0, 5.0, 4.0, 3.0, 3.0, 2.0, 2.0, 3.0, 4.0, 3.0, 3.0, 3.0, 4.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[3]=[0.0, 10.0, 11.0, 12.0, 11.0, 12.0, 11.0, 12.0, 12.0, 16.0, 16.0, 18.0, 22.0, 23.0, 27.0, 27.0, 29.0, 31.0, 38.0, 38.0, 40.0, 40.0, 39.0, 49.0, 57.0, 59.0, 64.0, 73.0, 85.0, 92.0, 92.0, 92.0, 100.0, 117.0, 123.0, 124.0, 125.0, 133.0, 155.0, 165.0, 181.0, 180.0, 179.0, 175.0, 170.0, 166.0, 167.0, 179.0, 174.0, 161.0, 160.0, 152.0, 148.0, 150.0, 143.0, 144.0, 131.0, 123.0, 118.0, 118.0, 108.0, 102.0, 92.0, 89.0, 88.0, 81.0, 76.0, 68.0, 63.0, 58.0, 54.0, 51.0, 50.0, 45.0, 45.0, 45.0, 46.0, 42.0, 37.0, 36.0, 33.0, 35.0, 31.0, 29.0, 26.0, 24.0, 24.0, 22.0, 21.0, 19.0, 18.0, 18.0, 14.0, 14.0, 14.0, 12.0, 12.0, 12.0, 11.0, 8.0, 7.0, 6.0, 5.0, 4.0, 4.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[4]=[0.0, 10.0, 8.0, 9.0, 7.0, 9.0, 9.0, 10.0, 12.0, 13.0, 17.0, 19.0, 21.0, 21.0, 24.0, 31.0, 31.0, 28.0, 34.0, 39.0, 40.0, 49.0, 52.0, 67.0, 78.0, 82.0, 91.0, 91.0, 109.0, 116.0, 122.0, 125.0, 137.0, 152.0, 161.0, 165.0, 172.0, 178.0, 185.0, 184.0, 192.0, 197.0, 191.0, 189.0, 186.0, 169.0, 159.0, 153.0, 152.0, 151.0, 145.0, 139.0, 143.0, 141.0, 130.0, 119.0, 108.0, 95.0, 84.0, 82.0, 81.0, 72.0, 68.0, 62.0, 60.0, 58.0, 45.0, 43.0, 43.0, 39.0, 33.0, 30.0, 30.0, 29.0, 24.0, 23.0, 18.0, 18.0, 17.0, 12.0, 11.0, 11.0, 11.0, 12.0, 10.0, 10.0, 9.0, 8.0, 7.0, 7.0, 6.0, 4.0, 4.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 3.0, 4.0, 4.0, 4.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[5]=[0.0, 9.0, 8.0, 11.0, 10.0, 13.0, 16.0, 14.0, 15.0, 12.0, 14.0, 17.0, 21.0, 24.0, 25.0, 26.0, 28.0, 26.0, 30.0, 39.0, 43.0, 45.0, 50.0, 50.0, 60.0, 70.0, 77.0, 85.0, 83.0, 98.0, 98.0, 104.0, 104.0, 114.0, 120.0, 136.0, 139.0, 129.0, 138.0, 145.0, 152.0, 152.0, 146.0, 142.0, 150.0, 142.0, 151.0, 150.0, 143.0, 145.0, 145.0, 140.0, 135.0, 131.0, 121.0, 123.0, 119.0, 105.0, 89.0, 86.0, 87.0, 83.0, 85.0, 80.0, 74.0, 67.0, 72.0, 69.0, 66.0, 56.0, 55.0, 46.0, 42.0, 42.0, 36.0, 34.0, 28.0, 27.0, 28.0, 28.0, 22.0, 20.0, 18.0, 17.0, 16.0, 16.0, 18.0, 16.0, 13.0, 15.0, 12.0, 13.0, 14.0, 11.0, 10.0, 8.0, 8.0, 7.0, 4.0, 5.0, 5.0, 4.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[6]=[0.0, 13.0, 13.0, 16.0, 16.0, 19.0, 21.0, 22.0, 27.0, 31.0, 36.0, 38.0, 40.0, 45.0, 48.0, 51.0, 54.0, 59.0, 67.0, 69.0, 85.0, 93.0, 90.0, 92.0, 104.0, 118.0, 123.0, 126.0, 129.0, 138.0, 144.0, 151.0, 160.0, 165.0, 165.0, 154.0, 156.0, 153.0, 151.0, 162.0, 158.0, 152.0, 137.0, 132.0, 130.0, 127.0, 126.0, 120.0, 118.0, 121.0, 120.0, 115.0, 104.0, 95.0, 88.0, 77.0, 72.0, 70.0, 70.0, 65.0, 64.0, 65.0, 63.0, 65.0, 61.0, 57.0, 58.0, 52.0, 49.0, 48.0, 46.0, 39.0, 35.0, 34.0, 29.0, 27.0, 28.0, 28.0, 25.0, 24.0, 23.0, 20.0, 18.0, 17.0, 15.0, 12.0, 10.0, 10.0, 8.0, 8.0, 7.0, 6.0, 6.0, 6.0, 7.0, 6.0, 4.0, 5.0, 5.0, 4.0, 4.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]    \n",
    "    solstore[7]=[0.0, 9.0, 11.0, 13.0, 13.0, 11.0, 14.0, 16.0, 16.0, 18.0, 16.0, 21.0, 26.0, 35.0, 33.0, 36.0, 39.0, 42.0, 45.0, 48.0, 51.0, 54.0, 57.0, 67.0, 74.0, 72.0, 80.0, 89.0, 101.0, 111.0, 112.0, 118.0, 135.0, 136.0, 145.0, 154.0, 158.0, 168.0, 172.0, 179.0, 174.0, 184.0, 188.0, 178.0, 164.0, 165.0, 154.0, 152.0, 152.0, 149.0, 143.0, 135.0, 134.0, 132.0, 128.0, 119.0, 110.0, 108.0, 102.0, 91.0, 84.0, 82.0, 79.0, 72.0, 65.0, 58.0, 57.0, 51.0, 46.0, 46.0, 40.0, 39.0, 39.0, 33.0, 29.0, 23.0, 25.0, 24.0, 22.0, 19.0, 16.0, 15.0, 15.0, 12.0, 13.0, 13.0, 11.0, 8.0, 8.0, 7.0, 6.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[8]=[0.0, 9.0, 10.0, 9.0, 9.0, 7.0, 6.0, 5.0, 7.0, 8.0, 9.0, 12.0, 16.0, 18.0, 16.0, 17.0, 17.0, 17.0, 21.0, 26.0, 29.0, 30.0, 34.0, 39.0, 36.0, 33.0, 34.0, 45.0, 49.0, 69.0, 71.0, 76.0, 80.0, 91.0, 98.0, 118.0, 130.0, 143.0, 136.0, 133.0, 137.0, 143.0, 155.0, 152.0, 158.0, 157.0, 158.0, 154.0, 159.0, 163.0, 169.0, 170.0, 172.0, 157.0, 164.0, 154.0, 150.0, 150.0, 144.0, 135.0, 128.0, 123.0, 113.0, 111.0, 102.0, 93.0, 90.0, 93.0, 87.0, 75.0, 69.0, 64.0, 58.0, 51.0, 47.0, 41.0, 37.0, 33.0, 29.0, 27.0, 29.0, 28.0, 25.0, 25.0, 20.0, 19.0, 17.0, 15.0, 14.0, 13.0, 12.0, 11.0, 12.0, 8.0, 7.0, 7.0, 6.0, 5.0, 4.0, 4.0, 4.0, 3.0, 2.0, 3.0, 2.0, 2.0, 3.0, 3.0, 3.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    solstore[9]=[0.0, 8.0, 9.0, 11.0, 10.0, 12.0, 12.0, 12.0, 11.0, 16.0, 18.0, 19.0, 19.0, 21.0, 20.0, 24.0, 29.0, 32.0, 34.0, 34.0, 40.0, 41.0, 41.0, 51.0, 55.0, 63.0, 68.0, 77.0, 85.0, 97.0, 97.0, 100.0, 107.0, 114.0, 127.0, 133.0, 145.0, 143.0, 153.0, 151.0, 162.0, 171.0, 180.0, 179.0, 170.0, 171.0, 171.0, 170.0, 170.0, 160.0, 164.0, 153.0, 151.0, 145.0, 134.0, 131.0, 132.0, 124.0, 115.0, 109.0, 92.0, 84.0, 81.0, 77.0, 67.0, 64.0, 61.0, 61.0, 54.0, 50.0, 46.0, 37.0, 35.0, 36.0, 38.0, 35.0, 31.0, 25.0, 21.0, 18.0, 19.0, 20.0, 16.0, 15.0, 15.0, 13.0, 12.0, 12.0, 13.0, 12.0, 10.0, 10.0, 9.0, 9.0, 7.0, 7.0, 7.0, 7.0, 7.0, 5.0, 6.0, 5.0, 5.0, 4.0, 5.0, 5.0, 5.0, 4.0, 4.0, 4.0, 4.0, 4.0, 3.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 3.0, 3.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
    "    \n",
    "    ### Functions\n",
    "\n",
    "    # Function to count numbers in each compartment\n",
    "    def count_type(type):\n",
    "        return np.count_nonzero(PERSON[:,3] == type)\n",
    "\n",
    "    # Function to check the current scale\n",
    "    def findscale(home):\n",
    "        # Check location\n",
    "        loc = NOHOUSES if home==1 else NOCLASSES\n",
    "        counts = np.zeros((loc, 4))\n",
    "        icounts = np.zeros((loc, 4))\n",
    "        icountss = np.zeros((loc,4))\n",
    "        for j in range(0, loc):\n",
    "            for c in range(1, 4):\n",
    "                counts[j, c] = np.count_nonzero(np.logical_and(PERSON[:,home]==j, PERSON[:,3]==c)) \n",
    "            for c in [0,2]:\n",
    "                icounts[j, c] = np.count_nonzero(np.logical_and(np.logical_and(PERSON[:,home]==j, PERSON[:,3]==c), PERSON[:,4]==0))\n",
    "               # icountss[j, c] = np.count_nonzero(np.logical_and(np.logical_and(PERSON[0:200,home]==j, PERSON[0:200,3]), PERSON[0:200,4]==0))\n",
    "        # Sum of all rates\n",
    "        scale = np.sum(OMEGA*counts[:,1] + GAMMA*counts[:,2] + WANE*counts[:,3]) + (1-MIX)*BETA[home-1]*np.sum(icounts[:,0]*icounts[:,2])+MIX*(BETA[home-1]*np.sum(icounts[:,2]))*np.sum(icounts[:,0])\n",
    "        return scale\n",
    "\n",
    "    # Main function\n",
    "\n",
    "    MIX=0.0001*(100-mix_w.value)\n",
    "    TRACING=(trace_w.value)\n",
    "    TESTNO=(10*test_w.value)  \n",
    "    CLASS=(class_w.value)\n",
    "    ISOLATE=(isolate_w.value)\n",
    "\n",
    "    classes=CLASS\n",
    "    NOCLASSES=int(np.ceil(POPSIZE/(classes))) \n",
    "\n",
    "    peaks=[] # For storing peak no. of infecteds\n",
    "    tots=[] # For storing total no. infected\n",
    "\n",
    "    for reps in range(0,REPS):\n",
    "        \n",
    "        if reps==0:\n",
    "            print(\"Running simulation 1\", end=\" \")\n",
    "        else:\n",
    "            print(\", %i\" %int(reps+1), end=\" \")\n",
    "\n",
    "        shuff1=np.arange(POPSIZE)\n",
    "        shuff2=np.arange(POPSIZE)\n",
    "\n",
    "        np.random.shuffle(shuff1)\n",
    "        np.random.shuffle(shuff2)\n",
    "\n",
    "        PERSON=np.zeros((POPSIZE,5))\n",
    "        # 0th col: ID - First 200 are superspreaders\n",
    "        # 1st col: house\n",
    "        # 2nd col: class\n",
    "        # 3rd col: SIR status\n",
    "        # 4th col: Isolation status\n",
    "        infs=np.sort(np.random.choice(POPSIZE,I0))\n",
    "        for i in range(0,POPSIZE):\n",
    "            PERSON[i][0] = i\n",
    "            if i in infs:\n",
    "                PERSON[i][3] = 2 # Initially everyone susceptible except I0 individuals\n",
    "\n",
    "        for index in range(NOHOUSES):\n",
    "            for i in shuff1[index::NOHOUSES]:\n",
    "                PERSON[i][1]=index\n",
    "\n",
    "        for index in range(NOCLASSES):\n",
    "            for i in shuff2[index::NOCLASSES]:\n",
    "                PERSON[i][2]=index \n",
    "\n",
    "        # Some local constants / lists\n",
    "        tsteps=[0]\n",
    "        infecteds=[count_type(2)]\n",
    "        susceptibles=[count_type(0)]\n",
    "        exposed=[count_type(1)]\n",
    "        current_t=0  \n",
    "        home=2 # Everyone starts at home\n",
    "        total_infections=0 # Since immunity wanes, count every infection event\n",
    "\n",
    "        # Main run\n",
    "        while current_t < MAX_TIME and np.any(PERSON[:,3] == 2):\n",
    "\n",
    "            # Find proposed time to next event\n",
    "            scale = findscale(home)\n",
    "            dt = -np.log(np.random.uniform()) / scale\n",
    "            proposed_t = tsteps[-1] + dt\n",
    "\n",
    "            if int(proposed_t) > int(current_t): # Has to come first or misses days\n",
    "                # If new day, change isolation status and test, then start again\n",
    "                # Update isolation status\n",
    "                for ppl in range(POPSIZE):\n",
    "                    if PERSON[ppl,4]>=1:\n",
    "                        PERSON[ppl,4]-=1\n",
    "\n",
    "                # Jump is only one new day\n",
    "                current_t = int(current_t)+1\n",
    "\n",
    "                # Weekly Test\n",
    "                if int(current_t)%7==1:\n",
    "\n",
    "                    # Only non-isolators test\n",
    "                    findx=np.where(PERSON[:,4]==0)\n",
    "                    # Adjust test number if too many isolating\n",
    "                    testnum = TESTNO if len(findx[0])>TESTNO else len(findx[0])\n",
    "                    # Randomise list\n",
    "                    np.random.shuffle(findx[0])\n",
    "                    # Isolate Exposed and Infected individuals for 10 days\n",
    "                    for tests in range(testnum):\n",
    "                        if PERSON[findx[0][tests],3]==1 or PERSON[findx[0][tests],3]==2:\n",
    "                            # False negative / ignore rate of 5%\n",
    "                            if np.random.uniform()<0.95:\n",
    "                                PERSON[findx[0][tests],4]=ISOLATE\n",
    "                                # contact trace\n",
    "                                num_traced = np.random.randint(TRACING+1)\n",
    "                                case_contacts=np.where((PERSON[:,1]==PERSON[findx[0][tests],1]) | (PERSON[:,2]==PERSON[findx[0][tests],2]))\n",
    "                                np.random.shuffle(case_contacts[0])\n",
    "                                #print(num_traced,case_contacts[0].size)\n",
    "                                for trace in range(num_traced):\n",
    "                                    if PERSON[case_contacts[0][trace],4]==0:\n",
    "                                        PERSON[case_contacts[0][trace],4]=ISOLATE\n",
    "                        else:\n",
    "                            # False positive rate of 1%\n",
    "                            if np.random.uniform()>0.99:\n",
    "                                PERSON[findx[0][tests],4]=ISOLATE\n",
    "                                # contact trace\n",
    "                                num_traced = np.random.randint(TRACING+1)\n",
    "                                case_contacts=np.where((PERSON[:,1]==PERSON[findx[0][tests],1]) | (PERSON[:,2]==PERSON[findx[0][tests],2]))\n",
    "                                np.random.shuffle(case_contacts[0])\n",
    "                                for trace in range(num_traced):\n",
    "                                    if PERSON[case_contacts[0][trace],4]==0:\n",
    "                                        PERSON[case_contacts[0][trace],4]=ISOLATE\n",
    "            elif home == 1 and proposed_t > int(proposed_t) + CLASS_START and proposed_t < int(proposed_t) + CLASS_END:\n",
    "                # If students are home and proposed time of next event is later \n",
    "                # than class starts, then no event occurs before class starts\n",
    "                current_t = int(proposed_t)+CLASS_START\n",
    "                home = 2\n",
    "            elif home == 2 and proposed_t > int(proposed_t) + CLASS_END:\n",
    "                # If students are in class and proposed time of next event is \n",
    "                # later than class ends, then no event occurs before class ends\n",
    "                current_t = int(proposed_t)+CLASS_END\n",
    "                home = 1 \n",
    "            else:\n",
    "                # Next event occurs before class starts/ends\n",
    "                current_t = proposed_t\n",
    "\n",
    "                # Find event\n",
    "                eventcheck = np.random.uniform()\n",
    "                # Need to find non-isolating infecteds, superspreaders and susceptibles\n",
    "                #countSS=np.count_nonzero(np.logical_and(PERSON[:,3]==2, PERSON[:,4]==0)) \n",
    "                countNSS=np.count_nonzero(np.logical_and(PERSON[:,3]==2, PERSON[:,4]==0)) \n",
    "                countsus=np.count_nonzero(np.logical_and(PERSON[:,3]==0, PERSON[:,4]==0))\n",
    "\n",
    "                if eventcheck < GAMMA*infecteds[-1]/scale: #Event is recovery  \n",
    "\n",
    "                    # If there are any infected people, randomly choose one to recover\n",
    "                    infected_indices = np.where(PERSON[:,3] == 2)\n",
    "                    if infected_indices:\n",
    "                        PERSON[np.random.choice(infected_indices[0]), 3] = 3\n",
    "\n",
    "                elif eventcheck < (GAMMA*infecteds[-1] + OMEGA*exposed[-1])/scale: # Event is latent->infected\n",
    "\n",
    "                    # If there are any latents, randomly choose one to become infected\n",
    "                    latent_indices = np.where(PERSON[:,3] == 1)\n",
    "                    if latent_indices:\n",
    "                        PERSON[np.random.choice(latent_indices[0]), 3] = 2                                                                                \n",
    "\n",
    "                elif eventcheck < (GAMMA*infecteds[-1] + OMEGA*exposed[-1]+WANE*count_type(3))/scale: \n",
    "                    # Event is waned immunity\n",
    "                    # If there are any immunes, randomly choose one to become susceptible\n",
    "                    immune_indices = np.where(PERSON[:,3] == 3)\n",
    "                    if immune_indices:\n",
    "                        PERSON[np.random.choice(immune_indices[0]), 3] = 0                                 \n",
    "\n",
    "                elif eventcheck <(GAMMA*infecteds[-1] + OMEGA*exposed[-1]+WANE*count_type(3) + MIX*(BETA[home-1]*countNSS)*countsus)/scale:\n",
    "\n",
    "                    # Event is mixed transmission\n",
    "\n",
    "                    sus_indices = np.where((PERSON[:,3] == 0) & (PERSON[:,4]==0))\n",
    "                    if sus_indices[0].size>0:\n",
    "                        PERSON[np.random.choice(sus_indices[0]), 3] = 1\n",
    "                        total_infections+=1\n",
    "\n",
    "                else: #Event is Transmission by househould contact\n",
    "                    findx=np.where((PERSON[:,3]==0) & (PERSON[:,4]==0)) # Find susceptible hosts who are not isolating\n",
    "                    np.random.shuffle(findx[0]) # randomly shuffle\n",
    "                    for tryx in findx[0]:\n",
    "                        loc=PERSON[tryx,home]\n",
    "                        contacts=np.where(np.logical_and(np.logical_and(PERSON[:,home]==loc, PERSON[:,3]==2), PERSON[:,4]==0))\n",
    "                        #contacts=np.where(PERSON[:,home]==loc)\n",
    "\n",
    "                        if contacts[0].size>0:\n",
    "                            PERSON[tryx,3]=1\n",
    "                            total_infections+=1\n",
    "                            break  \n",
    "\n",
    "            # Update lists\n",
    "            tsteps.append(current_t)\n",
    "            infecteds.append(count_type(2))\n",
    "            exposed.append(count_type(1))\n",
    "            susceptibles.append(count_type(0))\n",
    "\n",
    "            # Stop if infections has finished\n",
    "            if infecteds[-1]==0 & exposed[-1]==0:\n",
    "                tsteps.append(MAX_TIME)\n",
    "                infecteds.append(0)\n",
    "                exposed.append(0)\n",
    "                susceptibles.append(count_type(0))\n",
    "                break\n",
    "\n",
    "       # print(count_type(0), count_type(1), count_type(2), count_type(3))\n",
    "\n",
    "        # Find peak no. infected\n",
    "        peaks.append(max(infecteds)/POPSIZE)\n",
    "        tots.append(total_infections+I0)\n",
    "        infplot=[]\n",
    "        for i in range(len(infecteds)):\n",
    "            infplot.append(infecteds[i]/10)\n",
    "        if reps==0:\n",
    "            plt.plot(tsteps,infplot,'b:',label='Interventions')\n",
    "        else:\n",
    "            plt.plot(tsteps,infplot,'b:')\n",
    "    plt.plot(solstore[0]/10,'r:',label='Default')\n",
    "    for i in range(1,10):\n",
    "        plt.plot(solstore[i]/10,'r:')\n",
    "    plt.xlabel('Days')\n",
    "    plt.ylabel('% Infected')\n",
    "        \n",
    "    print(\"\\nThe mean peak infected was reduced from 18% to\", int(np.mean(peaks)*100), \"%.\")\n",
    "    print(\"The mean total infected was reduced from 96% to\", int(np.mean(tots)/10), \"%.\")\n",
    "\n",
    "display(test_w)\n",
    "display(isolate_w)\n",
    "display(trace_w)\n",
    "display(mix_w)\n",
    "display(class_w)\n",
    "interact_manual(runmodel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Caveats\n",
    "\n",
    "There are of course many important elements of a real epidemic that are not included in the model. These include:\n",
    "\n",
    "* Vaccination - Everyone in this model has the same vaccine status (which could be assumed to be no-one or everyone vaccinated). A more accurate model would account for variation in vaccine coverage.\n",
    "* Super-spreading - The model assumes everyone has the same potential to infect others. In many diseases 'super-spreading' plays a key role, where some individuals cause far more infections than others.\n",
    "* External contacts - The model assumes that the population is compltely closed, with no infection arriving in from other sources.\n",
    "* Testing accuracy - Testing accuracy in the model depends simply on your epidemic status, but a more realistic model would link this to infection growth dynamics in the body.\n",
    "\n",
    "There are many more simplifications and assumptions besides. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}