{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NOTEBOOK_HEADER-->\n",
    "*This notebook contains material from [CBE40455-2020](https://jckantor.github.io/CBE40455-2020);\n",
    "content is available [on Github](https://github.com/jckantor/CBE40455-2020.git).*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.1 Campus SIR Modeling](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html) | [Contents](toc.html) | [2.3 Campus Re-opening Model](https://jckantor.github.io/CBE40455-2020/02.03-Campus-reopening.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE40455-2020/blob/master/docs/02.02-Campus-SEIR-modeling.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7wQY1Qx7eOKX",
    "nbpages": {
     "level": 1,
     "link": "[2.2 Campus SEIR Modeling](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2-Campus-SEIR-Modeling)",
     "section": "2.2 Campus SEIR Modeling"
    }
   },
   "source": [
    "# 2.2 Campus SEIR Modeling\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nDB7fXWZeVef",
    "nbpages": {
     "level": 2,
     "link": "[2.2.1 Campus infection data](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.1-Campus-infection-data)",
     "section": "2.2.1 Campus infection data"
    }
   },
   "source": [
    "## 2.2.1 Campus infection data\n",
    "\n",
    "The following data consists of new infections reported since August 3, 2020, from diagnostic testing administered by the Wellness Center and University Health Services at the University of Notre Dame. The data is publically available on the [Notre Dame Covid-19 Dashboard](https://here.nd.edu/our-approach/dashboard/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1035,
     "status": "ok",
     "timestamp": 1597680824455,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "POjTcEnHSDwi",
    "nbpages": {
     "level": 2,
     "link": "[2.2.1 Campus infection data](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.1-Campus-infection-data)",
     "section": "2.2.1 Campus infection data"
    },
    "outputId": "3768ce83-7730-44c6-c8f0-ce301e462d64"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeff/opt/anaconda3/lib/python3.7/site-packages/pandas/plotting/_matplotlib/converter.py:103: FutureWarning: Using an implicitly registered datetime converter for a matplotlib plotting method. The converter was registered by pandas on import. Future versions of pandas will require you to explicitly register matplotlib converters.\n",
      "\n",
      "To register the converters:\n",
      "\t>>> from pandas.plotting import register_matplotlib_converters\n",
      "\t>>> register_matplotlib_converters()\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAEICAYAAACHwyd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZSElEQVR4nO3df5xldX3f8ddbFmRlkRUwU1yQJRFRZAuWEUk1dhYkoEahj4dBKNGlxWzbxMQ2GxV9aIyxrdiIir+abCV1o9aVohYiqYYHYTBWA7IKWX5IQFiEFRYlu8IQFFc//eOerZdxZu+dnTtzxjuv5+NxH3PPOd97vp97Z7/7nu85Z86kqpAkSfPrCW0XIEnSYmQAS5LUAgNYkqQWGMCSJLXAAJYkqQUGsCRJLTCApQUkyR8m+UTbdcxGkn+Z5J4kE0meO8D9/kqS2wa1P6ltBrB+7iXZkuTR5j/8+5N8LMmyFmt58RzteyxJJfnIpPVfTnLuXPTZ1e+9M3jJe4DXVdWyqvrGLPqtJM/YtVxVf1NVR+3p/qSFxgDWsHh5VS0DjgOeC7x5PjtPsmSeunoEeHWSlfPU3544HLi57SKkhc4A1lCpqvuBL9IJYgCSnJjkK0l2JLkxyVjXtvEk70pyXZKHklyW5MCu7a9IcnPz2vEkz+7atiXJm5L8HfBIkk8BTwf+opmNv7GP/o9Ick2Sh5NcCRzc4y3uAD4GvH26Bkn+TZJbk2xP8sUkhzfr35Hkg83zvZM8kuSPm+WlSX7Q/d53s//xJO9M8n+buv8qycFJnphkAtgLuDHJt5r2T0vymSTfTXJXkt/t2tdeSd6S5FvNvjYlOSzJl5omNzaf5asmz8STPLupZUfzPXpF17aPJflwkiua/V6b5JeabUnyviQPNN/zzUmO6fW+pYGrKh8+fq4fwBbgxc3zQ4HNwEXN8grgQeCldH7gPKVZfmqzfRzYChwD7Ad8BvhEs+2ZdGacpwB7A28E7gD26er3BuAwYOnkWvrs/6vAe4EnAi8CHt7V/xTvcwy4F/gnwEPAUc36LwPnNs9Pb2p8NrAEeCvwlWbbScDm5vk/B74FXNu17cbd9du1PN689pnA0mb5gq7tBTyjef4EYBPwB8A+wC8CdwKnNtvf0Hy/jgICHAscNHk/k+tovh93AG9p9ntS89nt+kw+1nzOJzSfwyeBjc22U5ualjd9Phs4pO1/xz4W38MZsIbF/07yMHAP8AA/nSH+BvCXVfWXVfWTqroSuJ5OIO7y8aq6qaoeAd4GnJlkL+BVwBVVdWVV/YjOuc2ldMJrlw9U1T1V9eg0dU3bf5KnA88D3lZVP6yqLwF/0euNVmeW/yfAH02x+d8B76qqW6tqJ/BfgOOaWfBXgSOTHEQn7C8GVjTny/8FcE2vvrv8j6r6++Z9X0LXEYdJnkfnh40/qqrHqupO4L8DZzXbXwu8tapuq44bq+rBPvo/EVhGJ/gfq6q/Bj4PnN3V5nNVdV3zOXyyq8YfAfsDzwLSfFb39f/WpcEwgDUszqiq/enMkp7FTw/lHg78enOYckeSHcALgUO6XntP1/O76cyuDgae1iwDUFU/adqumOa1U9ld/08DtjfB391/P94NnJrk2Cn6u6irr3+gM8tb0YTl9XTC9kV0AvcrwAuYeQDf3/X8H+mE4VQOB5426f2/BRhpth9GZzY9U08D7mm+J7vczeO/N1PW2IT1h4APAw8kWZ/kyXtQgzQrBrCGSlVdQ+fw43uaVffQmeEu73rsV1UXdL3ssK7nT6czQ/oe8B06AQJ0zh02bbd2dzm5hEnLu+v/PuApSfab1H8/7/NB4P3AO6fo799O6m9pVX2l2X4NncO1zwW+1iyfSudQ7ZcYvHuAuybVs39VvbRr+y/twX6/AxyWpPv/sKfz+O/NtKrqA1V1PHA0nUPpb9iDGqRZMYA1jN4PnNLMDj8BvDzJqc0FP/s2F/Mc2tX+N5IcneRJdA7rXlpVP6ZzaPVlSU5OsjewDvghnVnjdLbROc+5y7T9V9XddGak70iyT5IXAi+fwft8L53D4c/uWvcnwJuTPAcgyQFJfr1r+zXAa4BbquoxOudvX0snJL87g777dR3wcHOx2tLmMzgmyfOa7R8F3pnkyObiqH/aHCKHn/0su11LZ1b7xuaCsjE6n93GXgUleV6S5zff00eAHwA/6fEyaeAMYA2dJkj+HPiDqrqHzoVJbwG+S2fG9QYe/2//43RmzfcD+wK/2+znNjrncD9IZ0b8cjq/7vTYbrp/F/DW5nDr7/fR/78Cnk/nUPHbm7r7fZ8PAf8VOLBr3efoHJ7emOQh4CbgJV0v+wqd89i7Zru30AmguZj90vwg82t0zr/eRedz/ChwQNPkvXR+0PkrOheWXdzUB/CHwIbmszxz0n4fo/P9eEmzz48Ar6mqb/ZR1pPpnIfeTuew9YPAH+/ZO5T2XKomHzGTFo8k43SuOv5o27VIWlycAUuS1AIDWJKkFngIWpKkFjgDliSpBfN1A3kADj744Fq5cuW89ffII4+w33779W4oDTnHgtQx32Nh06ZN36uqp061bV4DeOXKlVx//fXz1t/4+DhjY2Pz1p+0UDkWpI75HgtJpr27nYegJUlqgQEsSVILDGBJklpgAEuS1AIDWJKkFhjAkiS1oK8ATvIfk9yc5KYkn2r+pNoRSa5NckeSTyfZZ66LlSRpWPQM4CQr6Px5ttGqOgbYCziLzp88e19VPYPOn/U6by4LlSRpmPR7CHoJsDTJEuBJwH3AScClzfYNwBmDL0+SpOHU805YVbU1yXuAbwOP0vnD2ZuAHVW1s2l2L7BiqtcnWQusBRgZGWF8fHwAZfdnYmJiXvuTFirHghaDzVu/37PNyFL44Ccvm3b7qhUHDLKk3eoZwEmeApwOHAHsAP4XcFq/HVTVemA9wOjoaM3nLcC8/Z7U4VjQYnDu+Vf0bLNu1U4u3Dx99G05Z2yAFe1eP4egXwzcVVXfraofAZ8FXgAsbw5JAxwKbJ2jGiVJGjr9BPC3gROTPClJgJOBW4CrgVc2bdYA08/pJUnS4/QM4Kq6ls7FVl8HNjevWQ+8Cfi9JHcABwEXz2GdkiQNlb7+HGFVvR14+6TVdwInDLwiSZIWAe+EJUlSCwxgSZJaYABLktQCA1iSpBYYwJIktcAAliSpBQawJEktMIAlSWqBASxJUgsMYEmSWmAAS5LUAgNYkqQWGMCSJLXAAJYkqQUGsCRJLTCAJUlqQc8ATnJUkhu6Hg8l+Q9JDkxyZZLbm69PmY+CJUkaBj0DuKpuq6rjquo44HjgH4HPAecDV1XVkcBVzbIkSerDTA9Bnwx8q6ruBk4HNjTrNwBnDLIwSZKG2UwD+CzgU83zkaq6r3l+PzAysKokSRpyqar+Gib7AN8BnlNV25LsqKrlXdu3V9XPnAdOshZYCzAyMnL8xo0bB1N5HyYmJli2bNm89SctVI4FLQabt36/Z5uRpbDt0em3r1pxwAArgtWrV2+qqtGpti2ZwX5eAny9qrY1y9uSHFJV9yU5BHhgqhdV1XpgPcDo6GiNjY3NoMvZGR8fZz77kxYqx4IWg3PPv6Jnm3WrdnLh5umjb8s5YwOsaPdmcgj6bH56+BngcmBN83wNcNmgipIkadj1FcBJ9gNOAT7btfoC4JQktwMvbpYlSVIf+joEXVWPAAdNWvcgnauiJUnSDHknLEmSWmAAS5LUAgNYkqQWGMCSJLXAAJYkqQUGsCRJLTCAJUlqgQEsSVILDGBJklpgAEuS1AIDWJKkFhjAkiS1wACWJKkFBrAkSS0wgCVJaoEBLElSCwxgSZJa0FcAJ1me5NIk30xya5JfTnJgkiuT3N58fcpcFytJ0rDodwZ8EfCFqnoWcCxwK3A+cFVVHQlc1SxLkqQ+9AzgJAcALwIuBqiqx6pqB3A6sKFptgE4Y66KlCRp2KSqdt8gOQ5YD9xCZ/a7CXg9sLWqljdtAmzftTzp9WuBtQAjIyPHb9y4caBvYHcmJiZYtmzZvPUnLVSOBS0Gm7d+v2ebkaWw7dHpt69accAAK4LVq1dvqqrRqbb1E8CjwN8CL6iqa5NcBDwE/E534CbZXlW7PQ88Ojpa119//YzfwJ4aHx9nbGxs3vqTFirHghaDledf0bPNulU7uXDzkmm3b7ngZYMsiSTTBnA/54DvBe6tqmub5UuBfwZsS3JI08EhwAODKFaSpMWgZwBX1f3APUmOaladTOdw9OXAmmbdGuCyOalQkqQhNP08/PF+B/hkkn2AO4F/TSe8L0lyHnA3cObclChJ0vDpK4Cr6gZgqmPYJw+2HEmSFgfvhCVJUgsMYEmSWmAAS5LUAgNYkqQWGMCSJLXAAJYkqQUGsCRJLTCAJUlqgQEsSVILDGBJklpgAEuS1AIDWJKkFhjAkiS1wACWJKkFBrAkSS0wgCVJaoEBLElSC5b00yjJFuBh4MfAzqoaTXIg8GlgJbAFOLOqts9NmZIkDZeZzIBXV9VxVTXaLJ8PXFVVRwJXNcuSJKkPszkEfTqwoXm+AThj9uVIkrQ4pKp6N0ruArYDBfxpVa1PsqOqljfbA2zftTzptWuBtQAjIyPHb9y4cZD179bExATLli2bt/6khcqxoMVg89bv92wzshS2PTr99lUrDhhgRbB69epNXUeOH6evc8DAC6tqa5JfAK5M8s3ujVVVSaZM8qpaD6wHGB0drbGxsf4rn6Xx8XHmsz9poXIsaDE49/wrerZZt2onF26ePvq2nDM2wIp2r69D0FW1tfn6APA54ARgW5JDAJqvD8xVkZIkDZueAZxkvyT773oO/CpwE3A5sKZptga4bK6KlCRp2PRzCHoE+FznNC9LgP9ZVV9I8jXgkiTnAXcDZ85dmZIkDZeeAVxVdwLHTrH+QeDkuShKkqRh552wJElqgQEsSVILDGBJklpgAEuS1AIDWJKkFhjAkiS1wACWJKkFBrAkSS0wgCVJaoEBLElSCwxgSZJaYABLktQCA1iSpBYYwJIktcAAliSpBQawJEktMIAlSWpB3wGcZK8k30jy+Wb5iCTXJrkjyaeT7DN3ZUqSNFxmMgN+PXBr1/K7gfdV1TOA7cB5gyxMkqRh1lcAJzkUeBnw0WY5wEnApU2TDcAZc1GgJEnDaEmf7d4PvBHYv1k+CNhRVTub5XuBFVO9MMlaYC3AyMgI4+Pje1zsTE1MTMxrf9JC5VjQYrBu1c6ebUaW7r7dfI6TngGc5NeAB6pqU5KxmXZQVeuB9QCjo6M1NjbjXeyx8fFx5rM/aaFyLGgxOPf8K3q2WbdqJxdunj76tpwzNsCKdq+fGfALgFckeSmwL/Bk4CJgeZIlzSz4UGDr3JUpSdJw6XkOuKreXFWHVtVK4Czgr6vqHOBq4JVNszXAZXNWpSRJQ2Y2vwf8JuD3ktxB55zwxYMpSZKk4dfvRVgAVNU4MN48vxM4YfAlSZI0/LwTliRJLTCAJUlqgQEsSVILDGBJklpgAEuS1AIDWJKkFhjAkiS1wACWJKkFBrAkSS0wgCVJaoEBLElSCwxgSZJaYABLktQCA1iSpBYYwJIktcAAliSpBQawJEkt6BnASfZNcl2SG5PcnOQdzfojklyb5I4kn06yz9yXK0nScOhnBvxD4KSqOhY4DjgtyYnAu4H3VdUzgO3AeXNXpiRJw6VnAFfHRLO4d/Mo4CTg0mb9BuCMOalQkqQhtKSfRkn2AjYBzwA+DHwL2FFVO5sm9wIrpnntWmAtwMjICOPj47MsuX8TExPz2p+0UDkWtBisW7WzZ5uRpbtvN5/jpK8ArqofA8clWQ58DnhWvx1U1XpgPcDo6GiNjY3tQZl7Znx8nPnsT1qoHAtaDM49/4qebdat2smFm6ePvi3njA2wot2b0VXQVbUDuBr4ZWB5kl3v4lBg64BrkyRpaPVzFfRTm5kvSZYCpwC30gniVzbN1gCXzVWRkiQNm34OQR8CbGjOAz8BuKSqPp/kFmBjkv8EfAO4eA7rlCRpqPQM4Kr6O+C5U6y/EzhhLoqSJGnYeScsSZJaYABLktQCA1iSpBYYwJIktcAAliSpBQawJEktMIAlSWqBASxJUgsMYEmSWmAAS5LUAgNYkqQWGMCSJLXAAJYkqQUGsCRJLTCAJUlqgQEsSVILDGBJklrQM4CTHJbk6iS3JLk5yeub9QcmuTLJ7c3Xp8x9uZIkDYd+ZsA7gXVVdTRwIvDbSY4GzgeuqqojgauaZUmS1IeeAVxV91XV15vnDwO3AiuA04ENTbMNwBlzVaQkScMmVdV/42Ql8CXgGODbVbW8WR9g+67lSa9ZC6wFGBkZOX7jxo2zr7pPExMTLFu2bN76kxYqx4IWg81bv9+zzchS2Pbo9NtXrThggBXB6tWrN1XV6FTb+g7gJMuAa4D/XFWfTbKjO3CTbK+q3Z4HHh0dreuvv34Gpc/O+Pg4Y2Nj89aftFA5FrQYrDz/ip5t1q3ayYWbl0y7fcsFLxtkSSSZNoD7ugo6yd7AZ4BPVtVnm9XbkhzSbD8EeGAQxUqStBj0cxV0gIuBW6vqvV2bLgfWNM/XAJcNvjxJkobT9PPwn3oB8Gpgc5IbmnVvAS4ALklyHnA3cObclChJ0vDpGcBV9WUg02w+ebDlSJK0OHgnLEmSWmAAS5LUAgNYkqQWGMCSJLXAAJYkqQUGsCRJLTCAJUlqQT834pAkacb6uTdzL4O+N/NC4gxYkqQWGMCSJLXAQ9CS1CIP0y5ezoAlSWqBASxJUgsMYEmSWmAAS5LUAgNYkqQWGMCSJLWgZwAn+bMkDyS5qWvdgUmuTHJ78/Upc1umJEnDpZ8Z8MeA0yatOx+4qqqOBK5qliVJUp96BnBVfQn4h0mrTwc2NM83AGcMuC5JkoZaqqp3o2Ql8PmqOqZZ3lFVy5vnAbbvWp7itWuBtQAjIyPHb9y4cTCV92FiYoJly5bNW3/SQuVY6N/mrd+f9T5WrThgwfY3nxbiZzmyFLY9Opj++rF69epNVTU61bZZ34qyqirJtCleVeuB9QCjo6M1NjY22y77Nj4+znz2Jy1UjoX+nTuIW0OeM7Zg+5tPC/GzXLdqJxdunj765vOz3NOroLclOQSg+frA4EqSJGn47WkAXw6saZ6vAS4bTDmSJC0O/fwa0qeArwJHJbk3yXnABcApSW4HXtwsS5KkPvU8B1xVZ0+z6eQB1yJJ0qLhnbAkSWqBASxJUgtm/WtIkjSXVg7iV1kueNkAKpEGyxmwJEktMIAlSWqBASxJUgsMYEmSWmAAS5LUAgNYkqQW+GtI0hDo9as661bt3O1fipnJr+n4a0HSYDgDliSpBQawJEkt8BC0JC0isz2F4OmDwXEGLElSCwxgSZJa4CFoaQ54pbCkXpwBS5LUglnNgJOcBlwE7AV8tKouGEhVasUwX5zhjFTSQrPHM+AkewEfBl4CHA2cneToQRUmSdIwm80h6BOAO6rqzqp6DNgInD6YsiRJGm6pqj17YfJK4LSqem2z/Grg+VX1uknt1gJrm8WjgNv2vNwZOxj43jz2Jy1UjgWpY77HwuFV9dSpNsz5VdBVtR5YP9f9TCXJ9VU12kbf0kLiWJA6FtJYmM0h6K3AYV3LhzbrJElSD7MJ4K8BRyY5Isk+wFnA5YMpS5Kk4bbHh6CrameS1wFfpPNrSH9WVTcPrLLBaOXQt7QAORakjgUzFvb4IixJkrTnvBOWJEktMIAlSWrBgg/gJGckqSTPGvB+35zkjiS3JTm1WbdvkuuS3Jjk5iTvGGSf0p6ai3GQ5KAkVyeZSPKhSduOT7K5GSMfSJJB9SvNxnyOhST7J7mh6/G9JO8fVL8LPoCBs4EvN18Horll5lnAc4DTgI80t9b8IXBSVR0LHAecluTEQfUrzcLAxwHwA+BtwO9Pse2/Ab8JHNk8Thtgv9JszNtYqKqHq+q4XQ/gbuCzg+p0QQdwkmXAC4Hz6ATmrvVjST7ftfyhJOc2z1+a5JtJNjU/uX9+8n7p3DJzY1X9sKruAu4ATqiOiabN3s3Dq9TUqrkaB1X1SFV9mc5/Pt39HQI8uar+tjpXaf45cMZcvDdpJuZ7LEzq+5nALwB/M6j3s6ADmE5QfqGq/h54MMnxu2ucZF/gT4GXVNXxwJS3/wJWAPd0Ld/brCPJXkluAB4Arqyqa2f5HqTZmqtxMJ0VdMbELv9/fEgtm++x0O0s4NM1wF8dWugBfDadP/JA87XXIYdnAXc2s1qAT820w6r6cXOo4VDghCTHzHQf0oDN+ziQFqg2x8JZs3z9z5jze0HvqSQHAicBq5IUnZt9VJI3ADt5/A8P+85w9z1vo1lVO5JcTefc100z3L80EHM8Dqazlc6Y2MXbzKp1LY2FXX0fCyypqk2D3O9CngG/Evh4VR1eVSur6jDgLuBX6JwIPzrJE5MsB05uXnMb8ItJVjbLr5pm35cDZzWvP4LORSbXJXlqsz+SLAVOAb45B+9N6tdcjoMpVdV9wENJTmyufn4NcNns34o0K/M+FrqczRwcSVqwM2A6b/jdk9Z9Bji7qv59kkvozEzvAr4BUFWPJvkt4AtJHqFzv+qfUVU3N6+/hc5PTr9dVT9uLj7Z0FwR/QTgkqqa6iIuab7M2TgASLIFeDKwT5IzgF+tqluA3wI+BiwF/k/zkNrU1lgAOBN46SDfDAzhrSiTLKuqieYn9w8Dt1fV+9quS5pPjgOpYyGPhYV8CHpP/WZzFfPNwAF0roCTFhvHgdSxYMfC0M2AJUn6eTCMM2BJkhY8A1iSpBYYwJIktcAAliSpBQawJEkt+H9GqDOuM9v5CQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "from scipy.integrate import solve_ivp\n",
    "from scipy.optimize import minimize\n",
    "from datetime import timedelta\n",
    "\n",
    "data = [\n",
    "    [\"2020-08-03\", 0],\n",
    "    [\"2020-08-04\", 0],\n",
    "    [\"2020-08-05\", 0],\n",
    "    [\"2020-08-06\", 1],\n",
    "    [\"2020-08-07\", 0],\n",
    "    [\"2020-08-08\", 1],\n",
    "    [\"2020-08-09\", 2],\n",
    "    [\"2020-08-10\", 4],\n",
    "    [\"2020-08-11\", 4],\n",
    "    [\"2020-08-12\", 7],\n",
    "    [\"2020-08-13\", 10],\n",
    "    [\"2020-08-14\", 14],\n",
    "    [\"2020-08-15\", 3],\n",
    "    [\"2020-08-16\", 15],\n",
    "    [\"2020-08-17\", 80],\n",
    "]\n",
    "\n",
    "df = pd.DataFrame(data, columns=[\"date\", \"new cases\"])\n",
    "df[\"date\"] = pd.to_datetime(df[\"date\"])\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,4))\n",
    "ax.bar(df[\"date\"], df[\"new cases\"], width=0.6)\n",
    "ax.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO))\n",
    "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %d\"))\n",
    "plt.title(\"Reported New Infections\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NCrQCtMfHqBG",
    "nbpages": {
     "level": 2,
     "link": "[2.2.2 Fitting an SEIR model to campus data](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.2-Fitting-an-SEIR-model-to-campus-data)",
     "section": "2.2.2 Fitting an SEIR model to campus data"
    }
   },
   "source": [
    "## 2.2.2 Fitting an SEIR model to campus data\n",
    "\n",
    "Because of  the limited amount of data available at the time this notebook was prepared, the model fitting has been limited to an SEIR model for infectious disease in a homogeneous population. In an SEIR model, the progression of an epidemic can be modeled by the rate processes shown in the following diagram.\n",
    "\n",
    "$$\\text{Susceptible}\n",
    "\\xrightarrow {\\frac{\\beta S I}{N}} \n",
    "\\text{Exposed} \n",
    "\\xrightarrow{\\alpha E} \n",
    "\\text{Infectious} \n",
    "\\xrightarrow{\\gamma I} \n",
    "\\text{Recovered} $$\n",
    "\n",
    "which yeild the following model for the population of the four compartments\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\frac{dS}{dt} &= -\\beta S \\frac{I}{N} \\\\\n",
    "\\frac{dE}{dt} &= \\beta S \\frac{I}{N} - \\alpha E \\\\\n",
    "\\frac{dI}{dt} &= \\alpha E - \\gamma I \\\\\n",
    "\\frac{dR}{dt} &= \\gamma I \\\\\n",
    "\\end{align*}$$\n",
    "\n",
    "The recovery rate is given by $\\gamma = 1/t_{recovery}$ where the average recovery time $t_{recovery}$ is estimated as 8 days. \n",
    "\n",
    "\n",
    "| Parameter | Description | Estimated Value | Source |\n",
    "| :-- | :-- | :-- | :-- |\n",
    "| $N$ | campus population | 15,000 | estimate |\n",
    "| $\\alpha$ | 1/average latency period | 1/(3.0 d) |\n",
    "| $\\gamma$ | 1/average recovery period | 1/(8.0 d) | literature |\n",
    "| $\\beta$ | infection rate constant | tbd | fitted to data | \n",
    "| $I_0$ | initial infectives on Aug 3, 2020 | tbd | fitted to data \n",
    "| $R_0$ | reproduction number | ${\\beta}/{\\gamma}$ |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 312
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1705,
     "status": "ok",
     "timestamp": 1597680825137,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "UtL_lSxiWk8H",
    "nbpages": {
     "level": 2,
     "link": "[2.2.2 Fitting an SEIR model to campus data](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.2-Fitting-an-SEIR-model-to-campus-data)",
     "section": "2.2.2 Fitting an SEIR model to campus data"
    },
    "outputId": "090b7006-6087-4f1c-93db-83ec1cc74f8a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R0 = 42.3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEWCAYAAAB/tMx4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhV1b3/8fc3yYEQZgERBYUyiCiTDc5DK2ql9qJtvVZvVWy12OLtYOmgt71VO9qr1ttbettS9JafotVSq3SglTrbVgUsg4IgUUSQIYBAIAQyfH9/7J1wEk6Sk+Sc7LOTz+t58uTsffbe50se+LCy9tprmbsjIiLxkxd1ASIi0joKcBGRmFKAi4jElAJcRCSmFOAiIjGlABcRiSkFuIhITCnARURiSgEuIhJTCnDJKWa23sy+YmYrzGy3mT1sZoVJ73/EzJaZ2S4z+7uZjQv3f8rMfp903Btm9puk7XfMbEIjn3lWeK1d4XHXhvsvNrN/mtmecP9tSecUmtkDZrYjPG+xmQ0M3+ttZvea2WYz22Rm3zWz/PC9EWb2bPhn225mD2f2JyidiQJcctHlwEXAMGAccC2AmU0E7gNuAPoBvwAWmFlX4FngbDPLM7OjgS7A6eF57wN6ACsafpCZHQcsBH4CDAAmAMvCt/cB1wB9gIuBz5nZpeF704DewJCwls8C+8P3fgVUASOAicCFwPXhe98BngD6AoPDzxVpFQW45KL/cfd33X0n8HuCUAWYDvzC3V9y92p3nwscAE5z9zeBsvDYc4C/AO+a2WjgXOB5d69J8Vn/BvzV3R9y90p33+HuywDc/Rl3X+nuNe6+AngovBZAJUFwjwhrWerue8JW+IeBL7n7PnffBtwDXJF03nHA0e5e4e4vZOynJp2OAlxy0Zak1+UErWcIgm9m2GWxy8x2EbSAjw7ffxb4AEGAPws8QxC454bbqQwBSlK9YWanmtnTZlZqZrsJWtn9w7fvJ/hP4tdm9q6Z/ZeZJcIaE8DmpBp/ARwZnvc1wICXzew1M/t0Wj8RkRQU4BIn7wDfc/c+SV9F7v5Q+H5tgJ8dvn6W5gP8HWB4I+89CCwAhrh7b+DnBOFL2Fq/3d3HAGcAHyHobnmH4LeC/kk19nL3E8Pztrj7Z9z9aIKuoP81sxGt/olIp6YAlzj5JfDZsGVsZtY9vNHYM3z/WeCDQDd33wg8T9CX3g/4ZyPXnAecb2aXm1mBmfVLutnZE9jp7hVmdgpBdwsAZvZBMxsb3pzcQ9A1UuPumwn6uO82s15hn/xwMzs3PO9fzWxweJn3AAdSde2INEsBLrHh7kuAzwCzCMJvHeENzvD9tcBeguDG3fcAbwJ/c/fqRq65gaDPeiawk+AG5vjw7RnAt82sDPgW8EjSqUcB8wnCezXBfx73h+9dQ3ATdVVY53xgUPjeJOAlM9tL0Lr/Yth/L9JipgUdRETiSS1wEZGYUoCLiMSUAlxEJKYU4CIiMVXQnh/Wv39/Hzp0aHt+pIhI7C1dunS7uw9ouL9dA3zo0KEsWbKkPT9SRCT2zOztVPvVhSIiElMKcBGRmFKAi4jEVLv2gadSWVnJxo0bqaioiLqUSBUWFjJ48GASiUTUpYhIJpSUwN13wwMPwN690KMHXHUVzJwJwxubP61l2vVR+uLiYm94E/Ott96iZ8+e9OvXDzNrt1pyibuzY8cOysrKGDZsWNTliEhbLVwIl10GlZXBV61EIviaPx+mTEn7cma21N2LG+5PqwvFzG4K5y5+1cweCpeTGmZmL5nZunDZqy5pV5OkoqKiU4c3gJnRr1+/Tv9biEiHUFIShHd5ef3whmC7vDx4vyTlNPQt0myAm9kxwBeAYnc/CcgnWF3kh8A97j6CYMa161pbRGcO71r6GYh0EHfffXhwN1RZCffc0+aPSvcmZgHQzcwKgCJgM3AewTSZAHOBSxs5N3NKSmDGDOjVC/Lygu8zZmTkfzIRkYx44IH0Avz++5s+Jg3NBri7bwLuAjYQBPduYCmwy92rwsM2AsekOt/MppvZEjNbUlpa2vpKFy6EceNgzhwoKwP34PucOcH+hQtbf+0kt912G3fddVej7z/22GOsWrUqI58lIh3Q3r2ZPa4J6XSh9AUuIVgh/GigO8EqJ2lx99nuXuzuxQMGHPYkaHrasU+pOQpwEWlSjx7NH9OS45qQThfK+cBb7l7q7pXAo8CZQJ+wSwVgMLCpzdU0Jst9St/73vcYNWoUZ511FmvWrAHgl7/8JZMmTWL8+PF8/OMfp7y8nL///e8sWLCAr371q0yYMIGSkpKUx4lIJ3bVVcFIk6YkEnD11W3+qHQCfANwmpkVWXCnbTLBUlFPA5eFx0wDHm9zNY3JYp/S0qVL+fWvf82yZcv405/+xOLFiwH42Mc+xuLFi1m+fDknnHAC9957L2eccQZTp07lzjvvZNmyZQwfPjzlcSLSic2cmV6A33RTmz8qnT7wlwhuVr4CrAzPmQ18Hfiyma0jWDQ2e8mVxT6l559/no9+9KMUFRXRq1cvpk6dCsCrr77K2WefzdixY5k3bx6vvfZayvPTPU5EOonhw4Nx3kVFhwd5IhHsnz8/Iw/zpDUKxd1vdffR7n6Su1/t7gfc/U13P8XdR7j7v7r7gTZX05h27FOqde211zJr1ixWrlzJrbfe2ugY7XSPE5FOZMoUWLECpk+vP2pu+vRgfwse4mlKPOZCyWKf0jnnnMNjjz3G/v37KSsr4/e//z0AZWVlDBo0iMrKSubNm1d3fM+ePSkrK6vbbuw4Eenkhg+HWbNg926org6+z5qVscfoIS4BnsU+pZNPPplPfOITjB8/nilTpjBp0iQAvvOd73Dqqady5plnMnr06Lrjr7jiCu68804mTpxISUlJo8eJiGRb5HOhrF69mhNOOKH5kzM8t0AuSvtnISKdSpvmQskJ7dSnJCKSSX9ft523tu+juibzjeXIp5Ntkdo+pVmzoq5ERKRZ1TXOp+cupqKyhsJEHi/9x/n07pa5KaPj0wIXEYmZDTvLqaisAaBH10RGwxsU4CIiWbNmy6ERa8cflblhzrUU4CIiWVIvwAf2yvj1FeAiIlmyZuueutejj+qZ8esrwDNs6NChbN++vc3HiEj8vV6vC0UBLiISCxWV1azfvg8AMxg5UH3gWbF+/XpGjx7Ntddey6hRo/jkJz/JX//6V84880xGjhzJyy+/zM6dO7n00ksZN24cp512GitWrABgx44dXHjhhZx44olcf/31JD8Y9cADD3DKKacwYcIEbrjhBqqrq6P6I4pIO1u3bS+1Q7+PPaKIoi6ZH7WdU+PAh978x6xde/0dFzf5/rp16/jNb37Dfffdx6RJk3jwwQd54YUXWLBgAd///vcZMmQIEydO5LHHHuOpp57immuuYdmyZdx+++2cddZZfOtb3+KPf/xj3XSyq1ev5uGHH+Zvf/sbiUSCGTNmMG/ePK655pqs/RlFJHfUv4GZ+e4TyLEAj9KwYcMYO3YsACeeeCKTJ0/GzBg7dizr16/n7bff5re//S0A5513Hjt27GDPnj0899xzPProowBcfPHF9O3bF4Ann3ySpUuX1s2tsn//fo488sgI/mQiEoU1Ww8FeDZuYIICvE7Xrl3rXufl5dVt5+XlUVVVRaK5ybQacHemTZvGD37wg4zWKSLxUP8GZuaHEEKOBXhz3RxROvvss5k3bx7/+Z//yTPPPEP//v3p1asX55xzDg8++CDf/OY3WbhwIe+99x4AkydP5pJLLuGmm27iyCOPZOfOnZSVlXHcccdF/CcRkfawNssP8UAaAW5mxwMPJ+16H/At4P+F+4cC64HL3f29zJeYG2677TY+/elPM27cOIqKipg7dy4At956K1deeSUnnngiZ5xxBsceeywAY8aM4bvf/S4XXnghNTU1JBIJfvrTnyrARTqB3eWVbNkTLO7SpSCPof26Z+VzWjSdrJnlEyxefCpwI7DT3e8ws5uBvu7+9abOb9N0sp2AfhYiHcNLb+7gE7NfBGDMoF786Ytnt+l6mZpOdjJQ4u5vA5cAc8P9c4FL21ShiEgH0R43MKHlAX4F8FD4eqC7bw5fbwEGpjrBzKab2RIzW1JaWtrKMkVE4iPbT2DWSjvAzawLMBX4TcP3POiHSdkX4+6z3b3Y3YsHDBiQ8trtuSpQrtLPQKTjSL6BOSoXAhyYArzi7lvD7a1mNggg/L6tNQUUFhayY8eOTh1g7s6OHTsoLCyMuhQRaSN3b7culJYMI7ySQ90nAAuAacAd4ffHW1PA4MGD2bhxI529e6WwsJDBgwdHXYaItNG7uysoq6gCoFdhAUf1yl7DLK0AN7PuwAXADUm77wAeMbPrgLeBy1tTQCKRYNiwYa05VUQk56zZkjyFbC/MLGuflVaAu/s+oF+DfTsIRqWIiEhozZa9da9HZekBnlqajVBEJIOSW+DZeoS+lgJcRCSDkocQZvMGJijARUQyprK6hpLSpC6ULE0jW0sBLiKSIeu376OyOhgSfXTvQnp3a9kspi2lABcRyZDX2+kBnloKcBGRDFnTTo/Q11KAi4hkSHvewAQFuIhIxqzZmjSEcGB2hxCCAlxEJCP2HajinZ37AcjPM4YfmZ1FHJIpwEVEMmBt0gRWw/p3p2tBftY/UwEuIpIB7X0DExTgIiIZUe8GZpYf4KmlABcRyYDkLhS1wEVEYkRdKCIiMVRadoAd+w4CUNQlnyF9i9rlcxXgIiJtlNz6HjmwJ3l52VvEIVlaAW5mfcxsvpm9bmarzex0MzvCzBaZ2Rvh977ZLlZEJBe9nrwKTzvdwIT0W+A/Bv7s7qOB8cBq4GbgSXcfCTwZbouIdDpR3MCENALczHoD5wD3Arj7QXffBVwCzA0Pmwtcmq0iRURyWRQ3MCG9FvgwoBT4PzP7p5nNCRc5Hujum8NjtgADU51sZtPNbImZLensK8+LSMdTU+Os3XpoEYdcC/AC4GTgZ+4+EdhHg+4Sd3fAU53s7rPdvdjdiwcMGNDWekVEcsqGneXsr6wGoH+PLvTv0bXdPjudAN8IbHT3l8Lt+QSBvtXMBgGE37dlp0QRkdy1JqL+b0gjwN19C/COmR0f7poMrAIWANPCfdOAx7NSoYhIDkvu/872GpgNFaR53OeBeWbWBXgT+BRB+D9iZtcBbwOXZ6dEEZHctaadF3FIllaAu/syoDjFW5MzW46ISLwkjwE//qjsL+KQTE9iioi0UkVlNet3lANgBqMG9mjXz1eAi4i0UknpXqprggF4xx5RRFGXdHulM0MBLiLSSlHewAQFuIhIq0V5AxMU4CIirfZ6RI/Q11KAi4i0UvIkVmqBi4jExO7ySjbvrgCgS34ex/Xr3u41KMBFRFoh+RH64Uf2IJHf/nGqABcRaYU1yYs4RNB9AgpwEZFWifoGJijARURaJapVeJIpwEVEWsjd67fAI3iIBxTgIiIttnl3BWUVVQD0LCxgUO/CSOpQgIuItFDDJzDNLJI6FOAiIi0U5So8ydKaOsvM1gNlQDVQ5e7FZnYE8DAwFFgPXO7u72WnTBGR3LEmB/q/oWUt8A+6+wR3r13Y4WbgSXcfCTxJg4WORUQ6qvpDCNt3EYdkbelCuQSYG76eC1za9nJERHJbZXUNJdv21m3HoQXuwBNmttTMpof7Brr75vD1FmBgqhPNbLqZLTGzJaWlpW0sV0QkWm/v2MfB6hoABvUupHdRIrJa0l0+4ix332RmRwKLzOz15Dfd3c3MU53o7rOB2QDFxcUpjxERiYtceAKzVlotcHffFH7fBvwOOAXYamaDAMLv27JVpIhIrsiVG5iQRoCbWXcz61n7GrgQeBVYAEwLD5sGPJ6tIkVEckUutcDT6UIZCPwuHKheADzo7n82s8XAI2Z2HfA2cHn2yhQRyQ1r4hTg7v4mMD7F/h3A5GwUJSKSi8oPVrFhZzkA+XnG8AE9Iq1HT2KKiKRp7dZDwweH9iuiMJEfYTUKcBGRtNVfxCG6B3hqKcBFRNKUSzcwQQEuIpK2XFjEIZkCXEQkTQ2nkY2aAlxEJA3b9x5g+96DAHRL5DOkb1HEFSnARUTSktz6HjWwB3l50SzikEwBLiKShly7gQkKcBGRtKzNkTnAkynARUTS8PrW3JnEqpYCXESkGTU1zhs5NoQQFOAiIs16571yyg9WA9CvexcG9OwacUUBBbiISDNyaQbCZApwEZFmKMBFRGIqF29gggJcRKRZsW+Bm1m+mf3TzP4Qbg8zs5fMbJ2ZPWxmXbJXpohINA5UVfPW9n1126Ni2gL/IrA6afuHwD3uPgJ4D7guk4WJiOSCkm37qK5xAI49oojuXdNZibJ9pBXgZjYYuBiYE24bcB4wPzxkLnBpNgoUEYnSmq2HFnHIpe4TSL8F/t/A14CacLsfsMvdq8LtjcAxqU40s+lmtsTMlpSWlrapWBGR9lZvDpQc6j6BNALczD4CbHP3pa35AHef7e7F7l48YMCA1lxCRCQyuXoDE9JYlR44E5hqZh8GCoFewI+BPmZWELbCBwObslemiEg01ubYIg7Jmm2Bu/st7j7Y3YcCVwBPufsngaeBy8LDpgGPZ61KEZEI7N5fybu7KwDokp/H0P7dI66ovraMA/868GUzW0fQJ35vZkoSEckNyWtgvm9AdxL5ufXoTIvGw7j7M8Az4es3gVMyX5KISG54PYe7T0BPYoqINGrNluQhhLmxiEMyBbiISCPWbtlb91otcBGRmHB3Xt+Suw/xgAJcRCSlLXsq2FMRPKvYs7CAQb0LI67ocApwEZEUGj6BGcwgklsU4CIiKazN4ScwaynARURSWJPjQwhBAS4iklJyF0ouzQGeTAEuItJAVXUN60qThxDm3hhwUICLiBxm/Y5yDlYFs2cf1auQ3kWJiCtKTQEuItJALk8hm0wBLiLSQPIj9Ll6AxMU4CIih4nDDUxQgIuIHGbNVnWhiIjETvnBKjbsLAcgP88YcWSPiCtqnAJcRCTJG1v34h68HtqviMJEfrQFNSGdRY0LzexlM1tuZq+Z2e3h/mFm9pKZrTOzh82sS/bLFRHJrriMQIH0WuAHgPPcfTwwAbjIzE4Dfgjc4+4jgPeA67JXpohI+6g/iVVuPsBTK51Fjd3dax9JSoRfDpwHzA/3zwUuzUqFIiLtaG1MbmBCmn3gZpZvZsuAbcAioATY5e5V4SEbgWMaOXe6mS0xsyWlpaWZqFlEJGtyfR3MZGkFuLtXu/sEYDDBQsaj0/0Ad5/t7sXuXjxgwIBWlikikn079h5g+94DABQm8jj2iKKIK2pai0ahuPsu4GngdKCPmdWuaj8Y2JTh2kRE2tWaBg/w5OXl3iIOydIZhTLAzPqEr7sBFwCrCYL8svCwacDj2SpSRKQ9NFyFJ9cVNH8Ig4C5ZpZPEPiPuPsfzGwV8Gsz+y7wT+DeLNYpIpJ1cbqBCWkEuLuvACam2P8mQX+4iEiHUP8GZm4PIQQ9iSkiAkBNjddrgY86Kncfoa+lABcRATa+t5/yg9UAHNG9CwN6dI24ouYpwEVEaDAD4cCemOX2CBRQgIuIAPUXcYjDDUxQgIuIAPF6ArOWAlxEOr3yg1W8+OaOuu1RCnARkXj41d/Xs33vQSBYhX7sMb0jrig9CnAR6dR276/k58+U1G1/8fyRJPLjEY3xqFJEJEt++dyb7KkIJlYd2q+Iy94/OOKK0qcAF5FOq7TsAPf97a267ZsuGBWb1jcowEWkE/vfZ9bVPbwz+qie/Mu4oyOuqGUU4CLSKW3atZ95L26o2/7Khcfn/PSxDSnARaRT+p+/vsHB6hoAJh7bh8knHBlxRS2nABeRTqekdC/zX9lYt/3VDx0fi0fnG1KAi0inc8+itVTXOABnjejPGcP7R1xR66SzIs8QM3vazFaZ2Wtm9sVw/xFmtsjM3gi/981+uSIibfPau7v5w4rNddtf+dDxEVbTNum0wKuAme4+BjgNuNHMxgA3A0+6+0jgyXBbRCSn3f3E2rrXF44ZyIQhfSKspm2aDXB33+zur4SvywjWwzwGuASYGx42F7g0W0WKiGTCkvU7eer1bQCYwcwL49v6hhb2gZvZUILl1V4CBrp77e8hW4CBjZwz3cyWmNmS0tLSNpQqItJ67s5//WVN3falE46JzbSxjUk7wM2sB/Bb4Evuvif5PXd3wFOd5+6z3b3Y3YsHDBjQpmJFRFrr+Te28/JbOwEoyDO+dP7IiCtqu7QC3MwSBOE9z90fDXdvNbNB4fuDgG3ZKVFEpG3cnTuTWt+fmDSE4/p1j7CizEhnFIoB9wKr3f1HSW8tAKaFr6cBj2e+PBGRtvvLa1tYuWk3AF0L8vj8efFvfQMUpHHMmcDVwEozWxbu+w/gDuARM7sOeBu4PDslioi0XnWNc1fSyJNpZwzlqN6FEVaUOc0GuLu/ADT2iNLkzJYjIpJZj/1zE+u27QWgR9cCPnvu8Igryhw9iSkiHdbBqhru+euh1vf1Zw/jiO5dIqwosxTgItJhPbx4Axvf2w9A36IE1501LOKKMksBLiId0v6D1fzPU+vqtmd8YAQ9CxOpDy4pgRkzoFcvyMsLvs+YEezPYQpwEemQ5v5jPaVlBwAY2KsrV59+XOoDFy6EceNgzhwoKwP34PucOcH+hQvbr+gWUoCLSOZF3KLdU1HJz5IWKv7C5JEUJvJT13nZZVBeDpWV9d+rrAz2X3ZZzrbEFeAiklk50KKd89yb7N4fBPKxRxRxefGQ1Afefffhwd1QZSXcc0+GK8wMBbiIZE4OtGi37z3AnBcOLVT85aYWKn7ggfQC/P77M1hh5ijARSRzcqBF+7NnSuoWKj5+YE/+ZXwTCxXv3ZveRdM9rp0pwEUkcyJu0b67az/3v/h23fbMC0eR39RCxT16pHfhdI9rZwpwEcmciFu0P3nqDQ5WBQsVjx/ShwvGpJzl+pCrroJEI0MLayUScPXVGaowsxTgIpI5EbZo39q+j0eWHFqo+GvpLFQ8c2Z6AX7TTRmoMPMU4CKSORG2aJMXKj5jeD/OHJHGQsXDh8P8+VBUdHjdiUSwf/784LgcpAAXkcyJqEW76t09LFj+bt12ixYqnjIFVqyA6dPrj1ufPj3YP2VKRmvNJAW4iGRORC3aHy06tFjD+ScM5ORj+7bsAsOHw6xZsHs3VFcH32fNytmWdy0FuEhnlo0nJtu5Rbv07ff46+rkhYpHZfT6ucyC5SzbR3FxsS9ZsqTdPk9EmrBwYfBQTWVl/aF/iUTwNX9+TncfQLBU2pW/fJEX3wzWurxkwtH8+IqJEVeVeWa21N2LG+5PZ0m1+8xsm5m9mrTvCDNbZGZvhN9b+PuKiEQqB56YzIS/rdtRF975ecZN53ee1jek14XyK+CiBvtuBp5095HAk+G2iMRFDjwx2VbBQsWv121fXjyEof3jv1BxSzQb4O7+HLCzwe5LgLnh67nApRmuS0SyKeZzgAA8sWoryzcGCxV3KcjjC5NHRFxR+2vtTcyB7r45fL0FaPRxJzObbmZLzGxJaWlpKz9ORDIq5nOAVNc4dz9xaOTJNacdx6De3SKsKBptHoXiwV3QRu+Euvtsdy929+IBAwa09eNEJBNiPgfIguWbWLs1+M+le5d8PveB3B7uly2tDfCtZjYIIPy+LXMliUjWxXgOkINVNdyz6I267evOfh/9enSNsKLotDbAFwDTwtfTgMczU46INCqTY7ZjOgfI1j0V3PjgK2zYWQ5An6IE15/dsRYqbol0hhE+BPwDON7MNprZdcAdwAVm9gZwfrgtItmS6VVuYjYHSE2N88CLb3P+3c+yaNXWuv2fO3c4vRpbqLgT0IM8IrmupCQI6fLyxo8pKgqecmxp4JaUBEMF778/uGHZo0fQbXLTTTkT3uu2lXHLoytZvP69evv/7dRj+fbUEylobLWdDqSxB3kKoihGRFqgJWO2Z81q2bVr5wBp6Xnt4EBVNT97poT/fbqEg9U1dfuH9e/O9z86ltOH94uwutygFrhIruvVK+guSee43buzX087WLJ+Jzc/upJ12w4NYyzIM2449318/rxGVpjvwNQCF4mrmI/Zbok9FZX8cOHrzHtpQ739E4b04Y6Pj2X0Ub0iqiw3KcBFcl2PHum1wHN0zHa6/vzqFm5d8Cpb9xyo29e9Sz5f/dDxXH360KbXtuykOn7vv0hjsjGVajauG+Mx2+nYuqeCG+5fwmcfWFovvCePPpJFXz6Xa88cpvBuhPrApXPK1lSq2bhuNkehRKimxnnw5Q38cOHrlB2oqtvfv0dXbps6hovHDmp+TctOorE+cAW4dD7ZCsRsBm0HmLs7WWNDA6+YNIRbppxA76LOO7Y7lVbPBy7S4WRrKtVsTtEa43Ubkx2oquaeRWuZ8uPn64X3sP7deegzp3HHx8c1Hd7Z6vaKKbXApfPJ1rC8TjjcryUWr9/JLSmGBn723OH8+3kjmh8a2MF+C2kJtcAlvjLd6srWsLxONNyvJfZUVPKN363kX3/+j3rhPWFIH/7whbP4yoeObz68O8gKQpmmAJfcluk5QCB7U6nGfIrWTKqpcUpK9/LQyxu44EfP1hvX3b1LPrdPPZHffu6M9Md1d4AVhLJBAZ7r4tbnl8l6s9XqytawvA4+3K8p28oqWLRqK3f9ZQ1X3/sSE779BJPvfpZbHl1Zb2jg+ScEQwOnndHCcd0dYAWhbFAfeC6LW59fpuudMSNoaTf1DzeRCG7ktWQujziOQskh+w5UsXLTbpa/s4tl7+xi+Tu7eHd3RZPn9O/RldunnsiHxx7VuqGBeXnBb1/pHFdd3fLr5zgNI4ybuIVBNurN5k3BOI0Dj1BldQ1rtpSxfGMQ1Mvf2c0b28qoSSM2jujehQlD+vD+4/py1anHtW1oYCe/QRzPm5hxeVIuG9fNdp9fpn8G2ag3mzcFszUsL8bD/dydDTvKWbD8Xb7zh1Vc9rO/M/a2v/CRn7zAN373Ko8s2ciaranDuzCRx6Shfbn+rGH85MqJPP+1D7L0m+dz37WTuPGDI9o+rrsTd081JXdb4HFrIWX6unFrfWaj3k7e6moLd6eisoZd+w+yq7ySXeWV7K59ve2VkiwAAAljSURBVL/Bdrhvy+79vFfezH/CQJ7BqIE9GT+4D+OH9GHCkD6MGtgju/Nyx+030gzLSheKmV0E/BjIB+a4e5Mr86Qd4HHro8zGdbPV55etn0E26s1WH3gOcHeqa5zKaudgdQ2VtV9VDbarazhY5fW2D1TVUFZRxe79lewqPxTKu8srDwX2/koOVtU0X0gajunTjQlD+jB+SG/GD+7DScf0pnvXCObB62DdUy2R8QA3s3xgLXABsBFYDFzp7qsaOyftAE/6h/v80Ak8MfK0w4/Jy4cxY+Ccc9Iv+rlnYdVqqGkiQDJwXSfFTZq8PDhxDJx9dnrXnHMvHDzY/HFdusD117Wg1udh1SqoOfSP2xveVMrLC34GZ51Vb3eTf1Puuw8OHvpH5SnvU1lQ77XBcqoN/+od9lexbA/+6O+g+tA8GbU/W7fwdUEBPnUq9OhZV2Pt32lPumbtpd390J/DoXbLHWrcqfHgmJrDtp2ammDfoWObOd6husY5WJUcyIcCuh1/+U1br8KCulb1+MF9GDekN0f2LIy6rENisIJQNmQjwE8HbnP3D4XbtwC4+w8aOyftAE/61fkXp3yMH3zw062qUaQz65KfR5+iRPDVrQu9ixL06RZuF3Whd7dD79Ued0yfbppAKgdlY0GHY4B3krY3Aqem+ODpwHSAY489Nr0rd7In1aRzMgtCtkt+HomCPBL5RqJ2Oz+PREGwfWhfuF2QR6/CAnrXBm8YxL2TgrhPty4UJvIUxh1c1juy3H02MBuCFnhaJyVNYH/m+mV8+4mfpT6usBDuviv9YmbOhIoDzR+Xgetaqg6HwsJgtEa6XlsFc34Z9Bkn9xvn5wdf138m6JZpiS/PhAMpxuw2/E2sWzf40Y8OO6zJOHjtNZg9u65eq71mbb033ICddFL96zW4YKrrH36M1R1sdcdY0utD5xh2+PkNjw238vOC9/LMyDPIs+DcvKR9VvteXvJ26uPNgrk+akM4kW8kCg4FtOa4ljZz91Z9AacDf0navgW4palz3v/+93taPvc590TCPYiV1F+JhPuNN6Z3vbhe19193brgvF693PPygu833hjsb41s1pqNekXEgSWeIlPb0gdeQHATczKwieAm5r+5+2uNnaNRKDkwxClOtYoIkIUHedy9Cvh34C/AauCRpsK7RYYPD4YEFRUdPng/kQj2z5/f8oCJ23WzIU61ikiT2jTy3t3/5O6j3H24u38vU0UB8XtSLk5P4MWpVhFpVO4+iSkiIkBc50IREZFGKcBFRGJKAS4iElPt2gduZqXA2608vT+wPYPlZFuc6lWt2ROneuNUK8Sr3rbWepy7D2i4s10DvC3MbEmqTvxcFad6VWv2xKneONUK8ao3W7WqC0VEJKYU4CIiMRWnAJ8ddQEtFKd6VWv2xKneONUK8ao3K7XGpg9cRETqi1MLXEREkijARURiKhYBbmYXmdkaM1tnZjdHXU9jzGyImT1tZqvM7DUz+2LUNTXHzPLN7J9m9oeoa2mOmfUxs/lm9rqZrQ6X9ctJZnZT+HfgVTN7yMxyaGFJMLP7zGybmb2atO8IM1tkZm+E3/tGWWOyRuq9M/y7sMLMfmdmfaKssVaqWpPem2lmbmb9M/FZOR/g4eLJPwWmAGOAK82shcvQtJsqYKa7jwFOA27M4VprfZFgOuA4+DHwZ3cfDYwnR+s2s2OALwDF7n4SkA9cEW1Vh/kVcFGDfTcDT7r7SODJcDtX/IrD610EnOTu4wjWJrilvYtqxK84vFbMbAhwIbAhUx+U8wEOnAKsc/c33f0g8GvgkohrSsndN7v7K+HrMoKAOSbaqhpnZoOBi4E5UdfSHDPrDZwD3Avg7gfdfVe0VTWpAOgWLnxSBLwbcT31uPtzwM4Guy8B5oav5wKXtmtRTUhVr7s/Ea5LAPAiMLjdC0uhkZ8twD3A1yDVeoutE4cAT7V4cs6GYi0zGwpMBF6KtpIm/TfBX6iaqAtJwzCgFPi/sMtnjpl1j7qoVNx9E3AXQUtrM7Db3Z+Itqq0DHT3zeHrLcDAKItpoU8DC6MuojFmdgmwyd2XZ/K6cQjw2DGzHsBvgS+5+56o60nFzD4CbHP3pVHXkqYC4GTgZ+4+EdhHbv2KXyfsO76E4D+do4HuZnZVtFW1TLgOYyzGGJvZNwi6L+dFXUsqZlYE/AfwrUxfOw4BvgkYkrQ9ONyXk8wsQRDe89z90ajracKZwFQzW0/QLXWemT0QbUlN2ghsdPfa32jmEwR6LjofeMvdS929EngUOCPimtKx1cwGAYTft0VcT7PM7FrgI8AnPXcfahlO8J/58vDf22DgFTM7qq0XjkOALwZGmtkwM+tCcDNoQcQ1pWRmRtBHu9rdfxR1PU1x91vcfbC7DyX4mT7l7jnbSnT3LcA7ZnZ8uGsysCrCkpqyATjNzIrCvxOTydEbrg0sAKaFr6cBj0dYS7PM7CKCLsCp7t7EKt3RcveV7n6kuw8N/71tBE4O/063Sc4HeFYXT868M4GrCVqzy8KvD0ddVAfyeWCema0AJgDfj7ielMLfEuYDrwArCf6d5dRj32b2EPAP4Hgz22hm1wF3ABeY2RsEv0XcEWWNyRqpdxbQE1gU/lv7eaRFhhqpNTuflbu/dYiISFNyvgUuIiKpKcBFRGJKAS4iElMKcBGRmFKAi4jElAJcOg0zu83MvhJ1HSKZogAXEYkpBbh0aGb2DTNba2YvAMeH+z5jZovNbLmZ/TZ8YrKnmb0VToWAmfWq3TazL4RzvK8ws19H+gcSSaIAlw7LzN5PME3ABODDwKTwrUfdfZK7184pfl04/e8zBNPrEp73aDiXyc3AxHDe6c+24x9BpEkKcOnIzgZ+5+7l4ayQtXPonGRmz5vZSuCTwInh/jnAp8LXnwL+L3y9guAR/qsIZr0TyQkKcOmMfgX8u7uPBW4HCgHc/W/AUDP7AJDv7rVLYl1MsCrUycDicJEGkcgpwKUjew641My6mVlP4F/C/T2BzWF/9ycbnPP/gAcJW99mlgcMcfenga8DvYEe7VG8SHM0mZV0aOFk/9MI5rbeQDBD4D6CaUhLCVZM6unu14bHHwW8BQxy911hyD9NENwGPODuOTNLn3RuCnCRJGZ2GXCJu18ddS0izVFfnkjIzH4CTCEYsSKS89QCFxGJKd3EFBGJKQW4iEhMKcBFRGJKAS4iElMKcBGRmPr/68+bsg44zhsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 15000         # estimated campus population\n",
    "gamma = 1/8.0     # recovery rate = 1 / average recovery time in days\n",
    "alpha = 1/3.0\n",
    "\n",
    "def model(t, y, beta):\n",
    "    S, E, I, R = y\n",
    "    dSdt = -beta*S*I/N\n",
    "    dEdt = beta*S*I/N - alpha*E\n",
    "    dIdt = alpha*E - gamma*I\n",
    "    dRdt = gamma*I\n",
    "    return np.array([dSdt, dEdt, dIdt, dRdt])\n",
    "\n",
    "def solve_model(t, params):\n",
    "    beta, I_initial = params\n",
    "    IC = [N - I_initial, I_initial, 0.0, 0.0]\n",
    "    soln = solve_ivp(lambda t, y: model(t, y, beta), np.array([t[0], t[-1]]), \n",
    "                IC, t_eval=t, atol=1e-6, rtol=1e-9)\n",
    "    S, E, I, R = soln.y\n",
    "    U = beta*S*I/N\n",
    "    return S, E, I, R, U\n",
    "\n",
    "def residuals(df, params):\n",
    "    S, E, I, R, U = solve_model(df.index, params)\n",
    "    return np.linalg.norm(df[\"new cases\"] - U)\n",
    "\n",
    "def fit_model(df, params_est=[0.5, 0.5]):\n",
    "    return minimize(lambda params: residuals(df, params), params_est, method=\"Nelder-Mead\").x\n",
    "\n",
    "def plot_data(df):\n",
    "    plt.plot(df.index, np.array(df[\"new cases\"]), \"r.\", ms=20, label=\"data\")\n",
    "    plt.xlabel(\"days\")\n",
    "    plt.title(\"new cases\")\n",
    "    plt.legend()\n",
    "\n",
    "def plot_model(t, params):\n",
    "    print(\"R0 =\", round(beta/gamma, 1))\n",
    "    S, E, I, R, U = solve_model(t, params)\n",
    "    plt.plot(t, U, lw=3, label=\"model\")\n",
    "    plt.xlabel(\"days\")\n",
    "    plt.title(\"new cases\")\n",
    "    plt.legend()\n",
    "\n",
    "plot_data(df)\n",
    "beta, I_initial = fit_model(df)\n",
    "plot_model(df.index, [beta, I_initial])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_bgKAAIBUrEj",
    "nbpages": {
     "level": 2,
     "link": "[2.2.3 Fitted parameter values](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.3-Fitted-parameter-values)",
     "section": "2.2.3 Fitted parameter values"
    }
   },
   "source": [
    "## 2.2.3 Fitted parameter values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 136
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 458,
     "status": "ok",
     "timestamp": 1597680829013,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "tWmFUnGsOEgn",
    "nbpages": {
     "level": 2,
     "link": "[2.2.3 Fitted parameter values](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.3-Fitted-parameter-values)",
     "section": "2.2.3 Fitted parameter values"
    },
    "outputId": "05c20701-888d-4fc4-aede-278883e7ffdd"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter              Value\n",
      "-----------  ---------------\n",
      "N            15000\n",
      "I0               2.30269e-05\n",
      "beta             5.28165\n",
      "gamma            0.125\n",
      "R0              42.2532\n"
     ]
    }
   ],
   "source": [
    "from tabulate import tabulate\n",
    "parameter_table = [\n",
    "                   [\"N\", 15000],\n",
    "                   [\"I0\", I_initial],\n",
    "                   [\"beta\", beta],\n",
    "                   [\"gamma\", gamma],\n",
    "                   [\"R0\", beta/gamma]\n",
    "]\n",
    "print(tabulate(parameter_table, headers=[\"Parameter\", \"Value\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XBWZxEBkHxBJ",
    "nbpages": {
     "level": 2,
     "link": "[2.2.4 Short term predictions of newly confirmed cases](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.4-Short-term-predictions-of-newly-confirmed-cases)",
     "section": "2.2.4 Short term predictions of newly confirmed cases"
    }
   },
   "source": [
    "## 2.2.4 Short term predictions of newly confirmed cases\n",
    "\n",
    "Using the fitted parameters, the following code presents a short term projection of newly diagnosed infections. Roughly speaking, the model projects a 50% increase per day in newly diagnosed cases as a result of testing sympotomatic individuals. \n",
    "\n",
    "The number of infected but asympotomatic individuals is unknown at this time, but can be expected to be a 2x multiple of this projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3407,
     "status": "ok",
     "timestamp": 1597681072898,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "hl001U_gRuNs",
    "nbpages": {
     "level": 2,
     "link": "[2.2.4 Short term predictions of newly confirmed cases](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.4-Short-term-predictions-of-newly-confirmed-cases)",
     "section": "2.2.4 Short term predictions of newly confirmed cases"
    },
    "outputId": "9b7e85cd-f4ae-4ae4-e340-745ac1e6f638"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAEICAYAAABcYjLsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwcdZ3/8dene64ckxsChABZCQJBAxIBOSN3gEDUXRYEwZNVwegunvktHnEV2XVR2CCKF+IFKBjCTQQCrMhCkDNcAZOQhBzkziSZo7u/vz++VdM1Pd0z3XNV9/T7+Xh0uq7u+nR1pvrd3/5WlTnnEBERERGpdom4CxARERERKQcKxiIiIiIiKBiLiIiIiAAKxiIiIiIigIKxiIiIiAigYCwiIiIiAigYi1QkM9vPzJyZ1Qzwer9pZr8ZoHUtMrNPDsS68qy7/XWa2T5m1mRmyR48zxwz+1nfV9gzZjbezB41s+1m9t8DuN7PmNm6YDuODe7/YaDWn1PLgP0fFpHKM6AfqiKSZWaXAR8F3gX83jn30VgLkrycc28Cw7tbzsymA79xzu0deex3+7G0nrgE2ACMcAN0EnszqwWuBo5yzj0XTO52e4qIxEHBWCQ+bwH/AZwGDIm5lkHLzGqcc6m46ygT+wIvDVQoDowHGoAlxSxsZknnXLp/SxIRyU9dKURi4py73Tk3H9jY3bJmljSz75vZBjP7O3BmzvyPmdnLwU/kfzezf4nMe9HMZkbGa4PnOazAuq4xs5Vmts3Mnjaz43IWqTOzm4J1LTGzaZHH7mVmt5nZ22a2zMxmR+YdYWZ/NbMtZrbGzOaZWV1k/ilm9oqZbTWzeYB1sT2+aWZ/NLNbgjr+ZmZTI/OXm9lXzOx5YIeZ1ZjZUWb2eLD+54IW3nD5SWb2SPBcC4FxkXkduq2Y2Rgz+6WZvWVmm81svpkNA+4F9gq6CTQF26LDz/ZmdnawzbYEXUUOyqn5i2b2fLANbjGzhmDeODO7K3jcJjN7zMzy7r/N7Ggzeyp4jqfM7Ohg+o3AxcCXg/pOzvPYIWb232a2Inj8/5rZkJ7WbmYHAK8Gi20xs4eC5Z2Z7R/WZWbXm9k9ZrYDeH/wfF8Knm+Hmf3cfDeQe4P36M9mNjqy/h69t3le/3QzW2Vml5vZ+uD/6cci8+vN/x2+ab5ryI8j2+cRM/tQMHxM8BrPDMZPMrNnC6wzab7LzRtBjU+b2cRgXsG/RfN/T4uDeevM7Ooit8dHze8jtpv/G72g0PYQqUrOOd100y3GG77V+MZulvk08AowERgDPAw4oCaYfybwDnyYPAHYCbwnmPdl4JbIc50DvNDFui4ExuJ/UbocWAs0BPO+CTQDZwBJ4ErgiWBeAnga+DpQB/wD8HfgtGD+4cBRwfPuB7wMfCGYNw7YDvwjUAv8K5ACPlmgxm8CbZHlvwgsA2qD+cuBZ4PtNQSYgP8CckZQ5ynB+G7B8n/F/9xfDxwf1PKbYN5+Odv6buAWYHSw7hOC6dOBVXnqDJ/nAGBHsO7a4H15HaiL1PwksFfwHr8MfDqYdyXw4+BxtcBxgOXZLmOAzcBHgu18fjA+Nph/I/AfXbz31wGLgu2VBI4Otklvau+w/YJpDtg/UtNW4JjgvWkInu8JfGvzBGA98DfgsGD+Q8A3gsf3+L3N8/qn4//fzQ1e5xn4v6XRwfwfAAuC19gI3AlcGcybC/xPMDwHeAO4KjLvmgLr/BLwAvBO/N/v1Mj71dXf4l+BjwTDw/FdVbrcHsAwYBvwzmDZPYEpce8DddOtnG6xF6CbbtV+o7hg/FAYNILxU3PDRs7y84HPB8N7BWFgRDD+R+DLJdS3GZgaDH8T+HNk3sHArmD4SODNnMd+Dfhlgef9AvCnYPgigoAdjBuwiq6DcXT5BLAGOC4YXw58PDL/K8Cvc57jfnwL6j5BGBoWmfc78gTjIEhkCIJSzvNNp+tgfAVwa07Nq4HpkZovjMz/T+DHwfBc4A6CMNnFe/UR4MmcaX8FPhoM30iBYBzUsyt8r3Pm9ab29u0XmZ8bjG/KWd9y4ILI+G3A9ZHxzwHze/veFngPd+XUuh7/hc7wXw7eEZn3PmBZMHwS8HwwfB/wSbJfGh8BPlhgna8C5/Tgb/FR4FvAuJxlutoew4AtwIeAIcWsUzfdqu2mrhQilWEvYGVkfEV0ppnNMLMngp/Zt+Bbi8YBOOfeAv4CfMjMRgEzgN8Gj1ti2Z/+jwumfdF8t4ytwXONpOPPz2sjwzuBBvPdDPbFdyXYEt7wLWfjg+c9IOgOsNbMtgHfjTxvh9fnnHM5rzef6PIZfJDeK9/8oLZ/yqntWHzQ3QvY7JzbEVm+w/aNmAhscs5t7qa2fPaKPm9Q80p8C18od9uGB6n9F76F9oHgZ/CvFrOOwIqcdRQyDt8a+0Yf116MfO/1usjwrjzj4fP31Xsb2ug69kkPX8tuwFDg6ch67gumg/8CcoCZjQcOBW4CJprZOOAIfJDNZyL5t3l3f4ufwLfkv2K+y8xZ3W2PYDv8M/4XqDVmdreZHdjN9hCpKjr4TqQyrMF/gIb2CQfMrB7fonYRcIdzrs3M5tOxj+6v8C1YNcBfnXOrAZxzU6IrCcLxl/GtX0uccxkz25zzXIWsxLeeTS4w/3rgGeB859x2M/sCvitEp9dnZpbzevOJLp8A9sYf0BiKHmC2Et+K9qncJzGzfYHRZjYsEqD2yXl89HnGmNko59yWnHndHdD2Fv4MJOF6w9e4upvH4Zzbjv8p/XIzOwR4yMyecs49mGcd++ZM2wcf4LqzAd9N5h3Acznzelx7kXpzMGBfvbfd2YAP5FPCv58o59xOM3sa+DzwonOu1cweB/4NeMM5t6GL+t8BvJhTe5d/i865pcD5wf/9DwJ/NLOxdLE9gsfdD9wf9I3+D+Cn+K45IoIOvhOJjfkDwhrwfTmT5g9WKvRl9VZgtpntbf6go2iLYR2+/+TbQMrMZuC7WkTNB96D/9C+qYuyGvE/Pb8N1JjZ14ERRb6kJ4Ht5g96GxIcVHSImb038tzbgKagleozkcfeDUwxsw8G22A2sEc36zs8svwXgBZ8v9R8fgPMNLPTgroazB9otbdzbgWwGPiWmdWZ2bHAzHxP4pxbgz/I7kdmNtr8gYzHB7PXAWPNbGSBGm4FzgwOxKrFB90W4PFuXidmdpaZ7R8E0q1AGt+lI9c9+FbLDwf/v/4Z393lru7WEbQC/wK42vyBg0kze1/wxavHtQ+APnlvuxNsn58CPzCz3QHMbIKZnRZZ7BHgsuAefH/t6Hg+PwO+bWaTzXt3EHC7/Fs0swvNbLegrvBLWqar7WH+AMZzzB8s2gI0kf//kUjVUjAWic+/41ugvoo/yGZXMC2fn+L7CT6HPwjp9nBG0Jo4Gx9eNgMfxh8gRGSZXfhW5UnRx+ZxP7518TX8T87NdN+lIVxHGjgL/zPyMnwL28/wP/+CP0Duw/j+zj/FH8AWPnYD8E/A9/AHCk3Gd//oyh34n4XDg80+6JxrK1DbSvxBh3PwQWMl/qCncB/4YXwf6U3AN+j6y8NH8Af+vYLvf/qFYB2vAL8H/h78hB3t1oFz7lX8+/w/+G0zE5jpnGvt5nWC3x5/xgeZvwI/cs49nOd1bsS/B5fjt+OXgbO6aK3M9UX8gWBP4bfFVUCil7X3qz5+b7vzFXyXlieC7kB/xh80F3oEH2gfLTCez9X4v90H8F8cf44/YLS7v8XTgSVm1gRcA5znnNvVzfZI4Fuw38JvjxPo+AVVpOqZ78onIoNd0OJ0gHPuwrhr6S0z+yb+4K2Kfy0iIlI+1MdYpAqY2Rj8wTofibsWERGRcqWuFCKDnJl9Cv9z6r3Oua5+0hUREalq6kohIiIiIoJajEVEREREgDLpYzxu3Di33377xbLuHTt2MGzYsFjWLSLVRfsbERko2t8U9vTTT29wzu2Wb15ZBOP99tuPxYsXx7LuRYsWMX369FjWLSLVRfsbERko2t8UZmYFr4CprhQiIiIiIigYi4iIiIgACsYiIiIiIoCCsYiIiIgIoGAsIiIiIgIoGIuIiIiIAArGIiIiIjKQfvITuOceaG6Ou5JOyuI8xiIiIiJSBZqb4fLLYccOGD4cXn0V9tor7qraqcVYRERERAbGQw/5UAyw557+VkYUjEVERERkYNxxR3b4nHPALL5a8lAwFhEREZH+l8nAggXZ8XPOia+WAhSMRURERKT/PfkkrF3rh3fbDd73vnjryUPBWERERET6X7QbxcyZkEzGV0sBCsYiIiIi0v9y+xeXIQVjEREREelfS5fCyy/74SFD4OST462nAAVjEREREelf0dbiU0+FoUPjq6ULCsYiIiIi0r/mz88Oz5oVXx3dUDAWERERkf6zfj08/rgfTiTgrLPiracLCsYiIiIi0n/uuguc88PHHAPjxsVbTxcUjEVERESk/1TA2ShCCsYiIiIi0j927IAHHsiOKxiLiIiISFVauBCam/3wlCmw//7x1tMNBWMRERER6R8V1I0CFIxFREREpD+kUnDnndlxBWMRERERqUqPPw4bN/rhPfeEadPiracICsYiIiIi0vdyu1Ekyj92dluhmU00s4fN7CUzW2Jmnw+mjzGzhWa2NLgfHUw3M7vWzF43s+fN7D39/SJEREREpIw4V3H9i6G4FuMUcLlz7mDgKOBSMzsY+CrwoHNuMvBgMA4wA5gc3C4Bru/zqkVERESkfL30Erzxhh9ubIT3vz/eeorUbTB2zq1xzv0tGN4OvAxMAM4BfhUs9isgvPD1OcBNznsCGGVme/Z55SIiIiJSnqKtxaefDvX18dVSgpI6e5jZfsBhwP8B451za4JZa4HxwfAEYGXkYauCaSIiIiJSDebPzw7PmlV4uTJTU+yCZjYcuA34gnNum5m1z3POOTNzpazYzC7Bd7Vg/PjxLFq0qJSH95mmpqbY1i0i1UX7GxEZKHHub+o2bODop54CIJNM8viIEaQqZN9XVDA2s1p8KP6tc+72YPI6M9vTObcm6CqxPpi+GpgYefjewbQOnHM3ADcATJs2zU2fPr1nr6CXFi1aRFzrFpHqov2NiAyUWPc3P/5x+2Bi+nSOPeuseOrogWLOSmHAz4GXnXNXR2YtAC4Ohi8G7ohMvyg4O8VRwNZIlwsRERERGcwq8GwUoWJajI8BPgK8YGbPBtPmAN8DbjWzTwArgHODefcAZwCvAzuBj/VpxSIiIiJSnrZtgwcfzI4PtmDsnPtfwArMPinP8g64tJd1iYiIiEilue8+aGvzw4cdBvvsE289JSr/S5CIiIiISGWo4G4UoGAsIiIiIn2hrQ3uvjs7rmAsIiIiIlXpkUdg61Y/vO++MHVqvPX0gIKxiIiIiPRebjcKK3SIWvlSMBYRERGR3nGu4vsXg4KxiIiIiPTWs8/CypV+eNQoOO64eOvpIQVjEREREemdaGvxmWdCbW18tfSCgrGIiIiI9M78+dnhWbPiq6OXFIxFREREpOeWL4fnnvPDdXVw2mmxltMbCsYiIiIi0nMLFmSHTzoJGhvjq6WXFIxFREREpOcGwdkoQgrGIiIiItIzmzb5C3uEZs6Mr5Y+oGAsIiIiIj1zzz2QTvvhI4+EvfaKt55eUjAWERERkZ4ZRN0oQMFYRERERHqiuRnuuy87rmAsIiIiIlXpoYegqckP778/HHRQvPX0AQVjERERESldtBvFrFlgFl8tfUTBWERERERKk8l0PH/xIOhGAQrGIiIiIlKqp56CtWv98G67wfveF289fUTBWERERERKE+1GcdZZkEzGV0sfUjAWERERkdLMn58dnjUrvjr6mIKxiIiIiBRv6VJ4+WU/PGQInHxyvPX0IQVjERERESletBvFqafC0KHx1dLHFIxFREREpHiD7Gp3UQrGIiIiIlKc9evhL3/xw4mEP/BuEFEwFhEREZHi3HUXOOeHjznGn6ptEFEwFhEREZHiDOJuFKBgLCIiIiLF2LkTFi7MjisYi4iIiEhVeuAB2LXLDx98MOy/f7z19AMFYxERERHpXrQbxSC6qEeUgrGIiIiIdC2d9gfehQZhNwpQMBYRERGR7jz+OGzY4If33BOmTYu3nn6iYCwiIiIiXYt2ozj7bH8O40FocL4qEREREekbzsH8+dnxQdqNAhSMRURERKQrL70Eb7zhh4cPhxNPjLeefqRgLCIiIiKFRbtRzJgB9fXx1dLPug3GZvYLM1tvZi9Gpn3TzFab2bPB7YzIvK+Z2etm9qqZndZfhYuIiIjIABjkV7uLKqbF+Ebg9DzTf+CcOzS43QNgZgcD5wFTgsf8yMySfVWsiIiIiAygt96CJ5/0w8kknHFG18tXuG6DsXPuUWBTkc93DnCzc67FObcMeB04ohf1iYiIiEhcFizIDk+fDqNHx1bKQKjpxWMvM7OLgMXA5c65zcAE4InIMquCaZ2Y2SXAJQDjx49n0aJFvSil55qammJbt4hUF+1vRGSg9NX+5l2//CVjg+GlBx/M6kG+DzPnXPcLme0H3OWcOyQYHw9sABzwbWBP59zHzWwe8IRz7jfBcj8H7nXO/bGr5582bZpbvHhxb15Hjy1atIjp06fHsm4RqS7a34jIQOmT/c22bbDbbtDa6seXL4d99+1tabEzs6edc3mvUNKjs1I459Y559LOuQzwU7LdJVYDEyOL7h1MExEREZFKct992VB86KGDIhR3p0fB2Mz2jIx+AAjPWLEAOM/M6s1sEjAZeLJ3JYqIiIjIgIuejWLWrPjqGEDd9jE2s98D04FxZrYK+AYw3cwOxXelWA78C4BzbomZ3Qq8BKSAS51z6f4pXURERET6RVsb3HNPdnyQn6Yt1G0wds6dn2fyz7tY/jvAd3pTlIiIiIjE6NFHYcsWP7zvvjB1arz1DBBd+U5EREREOop2ozj7bDCLr5YBpGAsIiIiIlnOwfz52fEq6UYBCsYiIiIiEvXss7BypR8eNQqOPz7eegaQgrGIiIiIZEW7UZx5JtTWxlfLAFMwFhEREZGsaDCuom4UoGAsIiIiIqHly31XCoC6Ojj99FjLGWgKxiIiIiLiLViQHT7pJGhsjK+WGCgYi4iIiIhXxd0oQMFYRERERAA2b4ZHHsmOz5wZXy0xUTAWEREREbj7bkin/fARR8Bee8VbTwwUjEVERESkYzeKWbPiqyNGCsYiIiIi1a6lBe67Lztehf2LQcFYRERERB56CJqa/PD++8NBB8VbT0wUjEVERESqXe7ZKMziqyVGCsYiIiIi1SyTqfrTtIUUjEVERESq2VNPwdq1fnjcODj66HjriZGCsYiIiEg1i7YWz5wJyWR8tcRMwVhERESkmqkbRTsFYxEREZFqtXQpvPSSHx4yBE45Jd56YqZgLCIiIlKtoq3Fp54KQ4fGV0sZUDAWERERqVbqRtGBgrGIiIhINVq/Hh5/3A8nEnDWWfHWUwYUjEVERESq0V13+XMYgz9F2267xVtPGVAwFhEREalG6kbRiYKxiIiISLXZuRMWLsyOKxgDCsYiIiIi1WfhQti1yw8ffDBMnhxvPWVCwVhERESk2qgbRV4KxiIiIiLVJJ2GO+/MjisYt1MwFhEREakmjz8OGzb44T33hPe+N956yoiCsYiIiEg1iXajOPtsfw5jARSMRURERKqHczB/fnZc3Sg6UDAWERERqRYvvQRvvOGHhw+HE0+Mt54yo2AsIiIiUi2i3ShmzID6+vhqKUMKxiIiIiLVQqdp65KCsYiIiEg1eOstePJJP5xMwhlnxFtPGeo2GJvZL8xsvZm9GJk2xswWmtnS4H50MN3M7Foze93Mnjez9/Rn8SIiIiJSpAULssMnnACjR8dXS5kqpsX4RuD0nGlfBR50zk0GHgzGAWYAk4PbJcD1fVOmiIiIiPSKulF0q9tg7Jx7FNiUM/kc4FfB8K+AWZHpNznvCWCUme3ZV8WKiIiISA9s3w4PPZQdVzDOq6aHjxvvnFsTDK8FxgfDE4CVkeVWBdPWkMPMLsG3KjN+/HgWLVrUw1J6p6mpKbZ1i0h10f5GRAZK7v5mt0WLmNLaCsD2/ffn6WXLYNmymKorXz0Nxu2cc87MXA8edwNwA8C0adPc9OnTe1tKjyxatIi41i0i1UX7GxEZKJ32Nz/7Wftg4wUXaF9UQE/PSrEu7CIR3K8Ppq8GJkaW2zuYJiIiIiJxaGuDu+/OjqsbRUE9DcYLgIuD4YuBOyLTLwrOTnEUsDXS5UJEREREBtqjj8KWLX54n33g0EPjraeMdduVwsx+D0wHxpnZKuAbwPeAW83sE8AK4Nxg8XuAM4DXgZ3Ax/qhZhEREREpVu7ZKMziq6XMdRuMnXPnF5h1Up5lHXBpb4sSERERkT7gnE7TVgJd+U5ERERksHr2WXjzTT88ahQcf3y89ZQ5BWMRERGRwSraWnzGGVBbG18tFUDBWERERGSwigbjWbMKLyeAgrGIiIjI4LRihe9KAVBXB6efHm89FaDXF/gQERERkTKRSsGSJYx8/nm4997s9BNPhMbG+OqqEArGIiIiIpVu0yaYNw+uvRZaWniXc9DcnJ1/8snx1VZBFIxFREREKtnSpf5sE1u2tIfhTgHvqqvg7LNh8uQBL6+SqI+xiIiISKXatAmOOw7WrevYQpxrwwYfnjdtGrjaKpCCsYiIiEilmjcPtm71F/LoinO+Rfm66wamrgqlYCwiIiJSiVIp36e4q5biqOZmuOYaSKf7t64KpmAsIiIiUomWLIGWltIe09oKL77YP/UMAgrGIiIiIpVo2zZIJkt7TCLhHyd5KRiLiIiIVKIRI0rvFpHJ+MdJXgrGIiIiIpVoyhSory/tMfX1cMgh/VPPIKBgLCIiIlKJampg9mxoaChu+YYGv3yp3S+qiIKxiIiISKW67DIYNQrMul7OzC936aUDU1eFUjAWERERqVRjxsCjj8LIkYWXaWiA8eP9cmPGDFxtFUjBWERERKSSrVkDO3dmxxMJUsOGQWMjjBsHc+b4U7vpctDd6nQpbRERERGpEM89BzNn+vMTA7zjHXDDDbzw4oscdsIJ/kA79SkumoKxiIiISCX6+9/htNOy5yXeYw9YuBAmTWJrIgFTp8ZbXwVSVwoRERGRSrN2LZxyCqxb58dHjoT774dJk+Ktq8IpGIuIiIhUkq1b4fTTfYsx+IPr7rwT3v3ueOsaBBSMRURERCpFczOcfbbvWwy+//Ctt8Jxx8Vb1yChYCwiIiJSCVIpOO88f9q10M9/7g++kz6hYCwiIiJS7pyDf/kXuOOO7LT/+i+4+OL4ahqEFIxFREREyt2cOfCLX2THv/Ql+OIX46tnkFIwFhERESlnV18N3/tedvxjH4OrroqvnkFMwVhERESkXN10E1x+eXb87LPhhhvALL6aBjEFYxEREZFydPfd8PGPZ8ePOw5uvhlqdH22/qJgLCIiIlJu/vIX+Kd/gnTaj7/73bBgAQwZEm9dg5yCsYiIiEg5eeEFOOss2LXLj0+aBPfdB6NGxVtXFVAwFhERESkXy5fDaafBli1+fPx4eOAB2HPPWMuqFgrGIiIiIuVg/Xo45RRYs8aPjxgB994L++8fb11VRMFYREREJG7btsGMGfD66368vt5fzOOww+Ktq8ooGIuIiIjEqbkZZs2Cv/3NjycS8Pvfw/TpsZZVjXp1vg8zWw5sB9JAyjk3zczGALcA+wHLgXOdc5t7V6aIiIjIIJROwwUXwMMPZ6f95CfwgQ/EV1MV64sW4/c75w51zk0Lxr8KPOicmww8GIyLiIiISJRz8NnPwu23Z6ddeSV88pPx1VTl+qMrxTnAr4LhXwGz+mEdIiIiIpXt61/3V7EL/eu/wle+El89gjnnev5gs2XAZsABP3HO3WBmW5xzo4L5BmwOx3MeewlwCcD48eMPv/nmm3tcR280NTUxfPjwWNYtItVF+xsRCU247TYmz5vXPr72lFN45atf9f2L+4D2N4W9//3vfzrS06GD3gbjCc651Wa2O7AQ+BywIBqEzWyzc250V88zbdo0t3jx4h7X0RuLFi1iujq3i8gA0P5GRAD47W/hwguz42ecAfPnQ21tn61C+5vCzKxgMO7V1xLn3Orgfj3wJ+AIYJ2Z7RmseE9gfW/WISIiIjJo3HcffPSj2fGjj4Y//KFPQ7H0XI+DsZkNM7PGcBg4FXgRWABcHCx2MXBHb4sUERERqXh//St86EOQSvnxKVPgzjth6NB465J2vTld23jgT74bMTXA75xz95nZU8CtZvYJYAVwbu/LFBEREalgS5bAmWfCzp1+fN994f77YcyYeOuSDnocjJ1zfwem5pm+ETipN0WJiIiIDBpvvgmnnQabg8s6jBsHDzwAEybEW5d0oivfiYiIiPSXDRvg1FNh9Wo/Pny472d8wAHx1iV5KRiLiIiI9Ift2/0ZJ1591Y/X1cEdd8Dhh8dblxSkYCwiIiLS11pa4IMfhKee8uNm/jRtJ54Yb13SJQVjERERkb6UTsPFF8Of/5yd9qMfwT/+Y3w1SVEUjEVERET6inMwezbcckt22ty58OlPx1eTFE3BWERERKSvzJ3rW4dDn/sc/Pu/x1ePlETBWERERKQv/OhH8M1vZsfPPx9++EPfv1gqgoKxiIiISG/deitcdll2/LTT4MYbIaGoVUn0bomIiIj0xsKFcOGFvn8xwJFHwm23+dOzSUVRMBYRERHpqaeegg98ANra/PhBB8Hdd8OwYfHWJT2iYCwiIiLSE6+8AjNmwI4dfnziRLj/fhg7Nt66pMcUjEVERERKtWqVv9Tzxo1+fOxYeOABH46lYikYi4iIiJRi40Yfileu9OPDhsE998CBB8Zbl/SagrGIiIhIsXbsgLPOgpdf9uO1tXD77XDEEfHWJX1CwVhERESkO+m07yoxYwY88YSfZgY33eRbj2VQqIm7ABEREZGy9fzzPvz+7newZk3HeddeC+edF09d0i8UjEVERESi1qzxQfimm3wwzmfu3I4X9JBBQcFYREREyoZz0NoK9fUDvOKdO2H+fB+GFy6ETKbzMrvvDh/+MFx8MRx66AAXKANBwVhEREQGlHPQ0gK7dvk8umtXxxvA8eFi1xoAABkdSURBVMf7Lrz9KpOBRYvg17+GP/4Rmpo6L9PQALNmwUc+AieeCK++Ctu2wXPPwZQpUKMoNZjo3RQREakGqRQsWeJD3YgR/R7qouE3vIUhuLk5f4NsVHMzDBnST8W99JIPw7/9bfaUa7lOOAEuugg+9CF/4N28eX68pQWSST+tvh5mz/ZdKsaM6adiZSApGIuIiAxmmzb5UHfttX0e6sJuD/lafXft6j78FlJX56+w3KfBeP16uPlm31Xi6afzL/POd/qW4QsugP3289OWLvXN11u2+LQe1dQEV14J118Pjz4Kkyf3YcESBwVjERGRwaqPQl205Tc3BPcm/A4Zkr0NHZodTiZ79pydNDfDggW+dfjee/0Xglxjx8L55/tA/N73duy/sWkTHHecD9XOFV7HunV+Oy9ZEmvLsXP+h4G2NtixI0k63YfbskooGIuIiPTWAHdTKEqJoa71mSXsbBiTt+U3X54sRm1t59Ab3vpt82Qy8Je/+JbhP/wBtm7tvExdHcyc6btGnH66H89n3jz/+ELbL+Sc//Jx3XVwxRW9fgnpdDbg5rsvNC/6Pi1d2sixx8LIkb0up6ooGIuIiPRUP3ZT6IlUypfR0gK1V85j2JatJIoIdelNW1g95zpWXFR6qAvDb76W3wH9brB0qW8Z/vWvYfny/Mscc4xvGT73XBg9uuvnS6X8+5rb0l5IczNccw3MmQPJZHvrbSnBNpzX01b4fC9BSqNgLCIi0hMD3Pe0rS0beqO3Xbv86nbs8GW0tUG6JcWsH19LoqW4UJdsbWbv265hxQVz8v72XlOTv9V36NCYG8Y3boRbbvFhOLwaXa53vMOH4Qsv9MPFWrIE19JCKSfGSO9qZcnvXmTbpKmxhFIz/3749ytFQtc3LpmCsYiIlL9y66rQx31P84XenTt94A1Db0uLXy685fv5PDRu9RISbS0lvSRra2X39S/C1KmdQnBtbUlP1b9aWuCee3xXibvv9hshR2bkaFLnfpj0+ReSOvxIUmnz3RPWZrsphPfR4ei04c9s4xCXLCkoZSxBevM2UhN79xITCb/Na2qy99HhQvOSyWwX6V27mrptFJfOFIxFRKR8lVlXhXYl9D11W7bQevV1bP/CFbS0+JAbDbzR0BuG3ba2nv2cbuYD0kjbBonSjrqqqUtw8N7b4ODS19tX0mn/ujsF1pQj9eTfSP/hdlIL7iG1bQdpkqR4Z3BfQyrZQPp9x5I69QzckUdl+w3/rfN6nPPryWT8cLje6HBL2whciZ2rLZMhNXRE+3ixwTZ3GbX0xkfBWESkvwWtnSOff973a4y7tTNXubXGhsrgNFnRoBreWnakmPCDa6kpsu+pNTfjfngNf9hjDq3pZLdZupCwFbG21me+6PCwYTB8uO/aUF8Pw14fQc2PSjxiLpPx738R0umOt0ym87SeLBNqD66r3yKz8CEyf36IzNq1ZEiQYSIZEjgSpEmQ2f+duGOPJ33k0WSGNfpwu7rr0FvMe2BuCkcn6qklz0U/CkgMqeeQ8w6hpt7/CfX7BUqkz5XBnk9EZJDKae18l3P+kzLu1s4C9ZVNa2xYWx+fJsu5ziE3egqysOtC2IobnpYslfLn6m1t9Y8Zt3oJn9nRUtIHaCLVSuOKF9kwYWre+clkNuiGYbe+3gfd4cN98B02zE8P54XDtbWdA5ibNAUa6mFH8aEuXVPPG3WHkH65c5jN7cIRhstoq2s0fEbH8y0TTssdJ5Vi9LLF1L7wDK2vLmPjW804EsCI4BYYPQYOPxymTYPx4/205uCG3x6JhN+u4bZNJrPTovOiwx2XqWHtubPZ53dXkmwt4ktQQwOJL8xmyHCdH62SKRiLSGWroNbO9qrK4aIAZdAa26VSuips3sLO/7yOjZ+9gp07Yft2fwvDbXi/c2c23Ia3fK2V3Rm3bRsZKy38uESCkbaN1CgfZhsafOgdOtQPNzT4MBb+nB4GtDCYhuF85878/WLD1xLOc66GA06azYF/upKatu5DXaqmgWeOnc3TjyTzBto+lU7793bTJn/bvJn6t1fynmXzOWzzQyRdCodhODIkeImDeJUDSDWMIHnUe0kcfyzJKQeSrEnkDbeJRHFdEQqF4+hw+jOX4e66HrdxHdbVhjCDUaPg0kv7bjtJLMrg00Okn5RrYKo05bodq6y1c7DVl8n4MNfa6ltlw1baHTtg57YUx3z/WuqK7arQ0oy75hq+v31Ol4E1GvKcy95yx6PBEHxQCltkzWBLZgQJV3rf09aGEe1Ximtpgc2bO64vd93RGnpi3bsvY9K915NIrevylG0ZM1qGjGLxkZdS5EksOrTIhkE0HLdUK4nNG0lufJvEpg0kNqwnsXE9yQ3rSby9lsTGt0m4FImgY8Rwmngvi6mhjSSd65xqL3DIsBW89tNHSR14SJdBNt9woWnFdXMYA48/WvhLJPhvNaNG+S+Ruix0xSuDT7eYqM9f75VrjeUcmHKV6zaE8t6Og6i1sy8vClC0HtaXyfj/Ck1N2RbZpqZsl4PwfseO7Hhzsx9uafHD4RkXWluzrZ25IXDS9iUcubOFApdcyCvR1krbMy+ybMTUTs8XDZ5hIApbFBMJ/ycX3kf77w4Z4v+7R7s51NRA7T5TyCysh7bSuimsGnUILs+1JvpLW+MYbvv8o3zo2uOp37mFmlTnUJeqbaBt6Cge+vqjjJkwpmDYDVuxa2qgtnUHNW+vIbH2LRKr3vL3a1eTXLOaxFurSGxcT4IMSd8DuMOt87Q0w2nCoOBp0ZIuTXLHFg7511PiubLc5Ml+vddd589T3Nqabc4P94eXXlo+nyvSK+b6/DeS0k2bNs0tXrx4YFaW82Gfco4a9fkbPDV2FZig4zf7OK9pX87bEMp7O27aBAcf3HVrJ/j0M378wH+QplKwxx7+/KpFcmPHkl69DpdIduh7Gf6MHh2Ojkenhy2vYRgNb+H0MJC27Uox+7t7MHRX8fVtrRnLjPesoy2TbO9yEIbO3BbXvvCeHY/xPytmMiJTfIpsSo7kikPv5OVxx7WH2zDgRvvj1tVlA3AYdEv5+T10xP1zee+fr6Q2T9jM1VbTwFMnz+HJ065oP51WuL7ocO64Wfb0W9FuFmGID2sP+89GlwmXSySgZtsmRvz6OkbeeA3W1grJBKR9qNv1qdm0fupSbGwQihPOt/KuepPEm8tJrFxB4s3l2JsrYEVw27Kl+A1VyB57wL77+m9Wr7xSXF+WhgZ/8YyB/BKZK52GF1/MNmYcckjZXnN50aJFTJ8+Pe4yypKZPe2cm5Z3XlUF43L+sK+E+qC8ayz3wBQq520IfbIdu/qZurvx7g7WGfbfcxk+78qiLlyQqWvg7U/NYd0lV3QKk4UOAOrqoKF8ATUcDq84NnL5c3zoh8dS11J8a+KumkaumP4Ybwyf2uFUVeGVs/IF5Og9FP+T+wG7nuOmZccyLFN8fTsSjVw06TFeG5L/wLGeiobAMJzW1MA7m59j3nPHMjRdfI0t9Y3cculjbJpYuMbcluJwPPc+38/v0UCaSED9jk0cc8kU6rZ03ffUmZEeN551Dy5pD5+FgnBX0/tEOg3PPOPD7fbt/sW8+WY28K5Y4cd37uzdepJJ2HtvH3yjt/328/cTJ/r9XA++RDJ2rO/iU6ZhtJwoGBfWVTAuk99sB0AZ9Kmr6PoqocZ583Bbt3Z9gAS0n1c09cPrSM+5Im+LV3Q438+xhablBrzoss6Bbd7EPjOOI7lxfeE6m5tx69aRet/xLLl1CakRYzo8R751Fltf9Dmi88PhTAYO/MM8Dty0lZoitmNq4xae+/h1PHPWFZ1Cbb7njx6tHgbAaMgLnyNsPAqvHBXOI5XiP35e/NW8Eq3NDPnpNXz7rTm0ZZKdjpjPrScahHODcb5+n7mtpgDv2bGNM9uSJXUDaMskePuNbbw+rIQH9dDwzDbSlBYq0iQYntkGdD6CP9pKGT0na9hq29CQvQ9v4RXTwrMrhMuGLaOJzBRqLq6H7cUH4+SQeo765CEk6zrWFd63938tEIbD4eJPrzUGpnbf99RGjaLm0UeZMLmIfWE6nT1KMHrEYL6jCHs63BeXY2togH32KRx899qruC5hS5b4b5OlaG31LbZT+/ZLmkio34KxmZ0OXAMkgZ85577XX+sqShn0+SsUmjIZSF49j9otxYW6zOYtbPvOdWyd3THU5ftAh/yBrdhwFQ2IE385j302byVZRI3pjVtY+rnrePXcKzq89q7CWzic22oXDUrhctFryadS/p9PX3UtQ0o4r2jb96/hqow/WKdQEMqtPfISO22vUO5R3NHxM5+ex96buw+d5vyR9q99/jrmv+uKvO9V+Nz56ihUc6HlwmFLp7jtL9fm7YuYT01bM++46xo+u2oOaZKdnjc6nu++VAfsWoK1lvZBmky1knn+Rd7o49bOQpoSI0hS2oFZSTI0JfypqHLDWfSgr+hwdFp4yw2EuaG1thb23j6C2pVpKOGgrqH1Gb763RHUTssG3Wjf09xW1mhNPVMDL8z2fcWL+ZtuaKDm32ZzwEGRwB/uODpc0izd+XQO4X2+acXMmzsX7rsPHnjAH00Y9j1NJv3pxA46CP7zP4sLr6WGxP4yYkTn0BsNvrvv3jcn6N22rfSW30TCP06kn/RLVwozSwKvAacAq4CngPOdcy/lW77fu1L04Oea7YmRfGrfe3FWQ8YZzvmdQJoEjvCD3cgE09vvMf9548Bh/ub8YQXO4R/bfoiBn5Z0KW5fdywjXfH9trbYKGbu9kT7Ediu/R/oOGgdpnV6t13Hejo/PpiTSfHApiMYVWKNJ476W96jxF1kffn/B3bc6Wb/m+bfGR+Qeombtn+AYewour4dDOOixj/xWs3AXOYp6VIs2noYo9zmoh+z2UYzfeSzJZ8aqqd6vh1vL7AdO7+7lvd/Z+4y4YDzASsYPazt//jh9k/S6LYXXd92a+RLI37Cc3XvjTx3psPzhn8F2WmufXr4K7ZFplkCEuba5xuOpGVIWoYa2rjxrdNK6h+7M9nIzw69ltqEo76mjTpLUZtIUZ9MUW9t1EXuGxJtNCRaqU+2MSTR6seTbdQkMkFNdPzmkTucycCdd/qWt2LV1sLJJ/vhrn4ayTe9u59TCk1PpeCNN4rre2rmm6Cj5zbrzekcBrNk0v+S11XwHTVqYGp57jk49ljfz7hYjY3w2GNqMS6CulIUFkdXiiOA151zfw8KuBk4B8gbjPtdD36uSWTaaFu2itd4Z4krcwWGCzuAV6mhhA8poNa1svv653tQX88cwKvU9qDGSZsXD0iNw1lNumDMyi+N0bh9JQl2o7v3LRqWLO+0zo+PVmM4DmAptZT2/7DONTN1y0O8zv7tz5toX0fHgGmdasnOz3ZRzHSqK1zWgIm8UOT/2iyH4x+2/41mWtqfJ9w2Fvkq6G/+XQq+XrYflZ69hlU4LXeeH96blSX/P6x3zUzfejuH8TA1tJLEh9cafxFZammjhjQ1tFEbjCdJUU8LNcHFZv1wdlpt8Pjwsb6+nhua3s7spz/Wi2foZ21tcO+9cVdRmHO+5bXSmfmreIQnOO6P4bpSOvn0sylTfF+aUoJxfb0/4E2kn/RXMJ4ArIyMrwKOjC5gZpcAlwCMHz+eRYsW9VMpMPL553mXcyW92DQJhpfQalaczm22BjSynXSJH6sZEoxiE/XBZX6svd23c2AK120d7jsPh8uHLWpRe7A2J1J1z2HszSq2MCoypePz54Y0i9RvHaZ1DH7WIXQ5xrCRWkrrO1dLigmsIhOcSD7R/vw+mPl1ZCKh07WH0ux0IsPha3LBc9ChxkksI1HK79f4n9hP437exQuR1x2e9sgBGZLt9WaC19ExTIavLTuciYRPPy1JmiQZGtlOA0WeyDQwhF18nmtJk2h/Hn+fe8ueoskiLbH9rY42zuOPA7Q2KScukcAlk53uCce7mpdnmU7zIstEnzNdX0+moYF0fT3phob24cyQIdl54bRgONPQQCbfJex6KjwNyebif6GKw74zZ7LP735HsohfMNJ1dbw5cyYrHntsACqrfE1NTf2arQar2A6+c87dANwAvitFvzb3jx5d8s6mvibNycenec+ojSQTYAnaQ0U4HoaRZHCKGzP8B785EklIWhBYzC9fk8iG02QCLPiZuPHtLQy9LQVtxdc3tDbFF87fyq7dl2MWtiRGfnY22qcbjoTldJiwIMAFExPtP1ln7yFbY93azQz7RZpSGuuG1af49qc3k5mwvP2o6vb1BeEI6DzPImG4QJ/KTtJ18LUkpXyXGTIsyXe+NwYSu4p8RHQL9iDWrRoLVztKaTRuqIcL/m2SP8J7IKTT8LXnS2p9Sw4byt5X/b/SD53vaQC4+27fp7OYg4hqamDGDDjzzOLW2Zfz162Dq64qfMBTeKLcr33Nd/XKfXyhjsZdzS9l2bVr4Rvf8O91W56dT22tb2m88kqYMKHj0Wq5t0LTe/KYQtMzGXj9db89R470Z04JTzKc26k6GDazEr/Oy4B797v93/O6dV0ffGBGcswYJn3/+0yK+5SlFUJdKXrIOdfnN+B9wP2R8a8BXyu0/OGHH+76VVubc2PHFur9lv82bpxzqVT/1lUp9VVKjd/6lnMNDcXV1tDg3Ny5A1ebc5WxDZ0r/+24caNze+zhnFnXtZn55TZuHNj6cmudO9e/742NrnXYMOcaG/37OnduvLXlqc+NHFle9Ul1eO01/7daaL/T0ODnv/Za3JVWlIcffjjuEsoWsNgVyKT99WvmU8BkM5tkZnXAecCCflpX92pq/EUTGhqKW76hwS8/UOdJLPf6oDJqvOwyf9BIMa16cVzTvhK2IZT/dhwzxp/jefz4wtuyocHPj/sSrWPG+LPbrFsHjz3Gi9/9rj9waO1aPz3ulq+c+rjzzvKqT6pDeGW5OXP8eYobG/2vAo2NMG6cn75kSbwXZZKq0W8X+DCzM4Af4k/X9gvn3HcKLTsgF/jYtMl39C/i55pYLv5Q7vVBZdRYCRfPKPdtCOW/HcFvywq7RKt+2hTpRgVdWa7caX9TmK58Fyr3D/tyrw8qo8ZyD0yVsA2h/LdjqII+SPVBJSIDRfubwhSMo3I+7NsyGWoTifL5sK+EMFIJNUJ5B6ZK2YZQ3tuxwuiDSkQGivY3hSkY5xN82D/zyCMcdsIJ5fdhXwlhpBJqLHfahlVFH1QiMlC0vyksjgt8lL9kEqZOZevmzeV5BZ2gvrJWCTWWO21DERGRsjFQ59gXERERESlrZdGVwszeBlbEtPpxwIaY1i0i1UX7GxEZKNrfFLavc263fDPKIhjHycwWF+pnIiLSl7S/EZGBov1Nz6grhYiIiIgICsYiIiIiIoCCMcANcRcgIlVD+xsRGSja3/RA1fcxFhEREREBtRiLiIiIiAAKxiIiIiIiwCAIxmY2y8ycmR3Yx8/7NTN73cxeNbPTgmkNZvakmT1nZkvM7Ft9uU4RKW/9sb8xs7Fm9rCZNZnZvJx5h5vZC8G+6Fozs75ar4iUt4Hc35hZo5k9G7ltMLMf9tV6K0nFB2PgfOB/g/s+YWYHA+cBU4DTgR+ZWRJoAU50zk0FDgVON7Oj+mq9IlL2+nx/AzQDVwBfzDPveuBTwOTgdnofrldEytuA7W+cc9udc4eGN/xF127vw/VWjIoOxmY2HDgW+AQ+yIbTp5vZXZHxeWb20WD4DDN7xcyeDlpg7sp9XuAc4GbnXItzbhnwOnCE85qCZWqDm45eFKkC/bW/cc7tcM79L/4DK7q+PYERzrknnD9K+iZgVn+8NhEpLwO9v8lZ9wHA7sBjffaCKkhFB2N8gL3POfcasNHMDu9qYTNrAH4CzHDOHQ7kvRwgMAFYGRlfFUzDzJJm9iywHljonPu/Xr4GEakM/bW/KWQCft8Tat8PicigN9D7m6jzgFtclZ62rNKD8fnAzcHwzXT/c8OBwN+DVmCA35e6QudcOviZYW/gCDM7pNTnEJGKNOD7GxGpWnHub87r5eMrWk3cBfSUmY0BTgTeZWYOSALOzL4EpOgY+htKfPrVwMTI+N7BtHbOuS1m9jC+z9+LJT6/iFSQft7fFLIav+8JddoPicjgE9P+Jlz3VKDGOfd0Xz5vJankFuN/BH7tnNvXObefc24isAw4Dt9p/GAzqzezUcBJwWNeBf7BzPYLxv+5wHMvAM4LHj8Jf9DLk2a2W/B8mNkQ4BTglX54bSJSXvpzf5OXc24NsM3MjgrORnERcEfvX4qIlLkB399EnE8VtxZDBbcY49+8q3Km3Qac75z7jJndim/JXQY8A+Cc22VmnwXuM7MdwFP5ntg5tyR4/Ev4b2eXOufSwcEwvwrOUJEAbnXO5Tt4T0QGl37b3wCY2XJgBFBnZrOAU51zLwGfBW4EhgD3BjcRGdzi2t8AnAuc0ZcvptJU3SWhzWy4c64paIG5DljqnPtB3HWJyOCj/Y2IDBTtb/pGJXel6KlPBWeVWAKMxB/FKSLSH7S/EZGBov1NH6i6FmMRERERkXyqscVYRERERKQTBWMRERERERSMRUREREQABWMREREREUDBWEREREQEgP8PKb7Z6hUF/yYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# prediction horizon (days ahead)\n",
    "H = 1\n",
    "\n",
    "# retrospective lag\n",
    "K = 6\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(12, 4))\n",
    "for k in range(0, K+1):\n",
    "    # use data up to k days ago\n",
    "    if k > 0:\n",
    "        beta, I_initial = fit_model(df[:-k])\n",
    "        P = max(df[:-k].index) + H\n",
    "        c = 'b'\n",
    "        a = 0.25\n",
    "    else:\n",
    "        beta, I_initial = fit_model(df)\n",
    "        P = max(df.index) + H\n",
    "        c = 'r'\n",
    "        a = 1.0\n",
    "\n",
    "    # simulation\n",
    "    t = np.linspace(0, P, P+1)\n",
    "    S, E, I, R, U  = solve_model(t, [beta, I_initial])\n",
    "\n",
    "    # plotting\n",
    "    dates = [df[\"date\"][0] + timedelta(days=t) for t in t]\n",
    "    ax.plot(dates, U, c, lw=3, alpha=a)\n",
    "\n",
    "ax.plot(df[\"date\"], df[\"new cases\"], \"r.\", ms=25, label=\"new infections (data)\")\n",
    "ax.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO))\n",
    "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %d\"))\n",
    "ax.grid(True)\n",
    "ax.set_title(f\"{H} day-ahead predictions of confirmed new cases\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "nkJN6e5mkj-I",
    "nbpages": {
     "level": 2,
     "link": "[2.2.4 Short term predictions of newly confirmed cases](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html#2.2.4-Short-term-predictions-of-newly-confirmed-cases)",
     "section": "2.2.4 Short term predictions of newly confirmed cases"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.1 Campus SIR Modeling](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html) | [Contents](toc.html) | [2.3 Campus Re-opening Model](https://jckantor.github.io/CBE40455-2020/02.03-Campus-reopening.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE40455-2020/blob/master/docs/02.02-Campus-SEIR-modeling.ipynb\"> <img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open in Google Colaboratory\"></a><p><a href=\"https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.ipynb\"> <img align=\"left\" src=\"https://img.shields.io/badge/Github-Download-blue.svg\" alt=\"Download\" title=\"Download Notebook\"></a>"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyMLHLafHBVRbuqsiaP/J00M",
   "collapsed_sections": [],
   "name": "Campus-SEIR-modeling.ipynb",
   "provenance": [
    {
     "file_id": "14v51fHgXJppd-iWZ_I8hWcQsqultykEX",
     "timestamp": 1597525384103
    }
   ]
  },
  "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
}