{
 "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.0 Modeling](https://jckantor.github.io/CBE40455-2020/02.00-Modeling.html) | [Contents](toc.html) | [2.2 Campus SEIR Modeling](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE40455-2020/blob/master/docs/02.01-Campus-SIR-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.01-Campus-SIR-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.1 Campus SIR Modeling](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1-Campus-SIR-Modeling)",
     "section": "2.1 Campus SIR Modeling"
    }
   },
   "source": [
    "# 2.1 Campus SIR Modeling\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nDB7fXWZeVef",
    "nbpages": {
     "level": 2,
     "link": "[2.1.1 Campus infection data](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.1-Campus-infection-data)",
     "section": "2.1.1 Campus infection data"
    }
   },
   "source": [
    "## 2.1.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": 527,
     "status": "ok",
     "timestamp": 1597709326245,
     "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.1.1 Campus infection data](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.1-Campus-infection-data)",
     "section": "2.1.1 Campus infection data"
    },
    "outputId": "07f05dd4-c7b2-4113-d7b3-9a436aa2c5ec"
   },
   "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": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAEICAYAAACK3Vc9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAWl0lEQVR4nO3de5SlVX3m8e8jiLS0gohTg01Lk4gXQo8aSmRGY6pBB4QozFoGYYxCBtMzE42uCV7QpTHGlREn8YJRk+lIAl6WLYNmIOKNRSyMY0S7o6QFJLZcAh0uakBogmLrb/44b+uhrOqqOqcuu6q+n7XOqvPe9t7vqd39nL3PW+9JVSFJktrwkMVugCRJ+hmDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLC0RSX4/yYcXux3DSPKfktySZGeSp81hub+S5Pq5Kk9aTAazlrUkNyW5vwuC25NckGT1IrblOfNU9liSSvL+Ceu/mOTM+aizr95bZ3HIHwOvqKrVVfW1IeqtJI/fvVxVf1tVTxy0PKklBrNWgudX1WrgqcDTgNcvZOVJ9l6gqu4DXpJk3QLVN4hDgWsWuxFSywxmrRhVdTvwWXoBDUCSY5J8KcndSa5OMta3bTzJ25J8Jck9SS5JcmDf9hckuaY7djzJk/u23ZTkdUn+AbgvyUeBxwF/3Y3eXzuD+g9LcmWSe5NcDhw0zSneDVwAvHmqHZL8lyTXJbkryWeTHNqtf0uSP+mePzTJfUn+qFteleQH/ee+h/LHk7w1yf/r2v25JAcleViSncBewNVJvt3t/9gkH0/ynSQ3JnllX1l7JXlDkm93ZW1NsjbJF7pdru5eyxdNHLkneXLXlru739EL+rZdkOR9SS7ryr0qyS9225LkXUnu7H7n25IcOd15S3Oqqnz4WLYP4CbgOd3zQ4BtwHnd8hrge8CJ9N6kPrdbfky3fRzYARwJ7Ad8HPhwt+0J9EaozwUeCrwW2A7s01fv14G1wKqJbZlh/X8HvBN4GPBs4N7d9U9ynmPArcC/Be4Bntit/yJwZvf85K6NTwb2Bt4IfKnbdiywrXv+H4BvA1f1bbt6T/X2LY93xz4BWNUtn9u3vYDHd88fAmwFfg/YB/gF4Abg+G77a7rf1xOBAE8BHj2xnInt6H4f24E3dOUe2712u1+TC7rX+ejudfgIsLnbdnzXpgO6Op8MHLzY/djHyno4YtZK8H+T3AvcAtzJz0aUvwF8qqo+VVU/qarLgS30gnK3D1XVN6rqPuBNwKlJ9gJeBFxWVZdX1Y/ofXa6il6o7faeqrqlqu6fol1T1p/kccDTgTdV1Q+r6gvAX093otWbFfgz4A8m2fzfgLdV1XVVtQv4n8BTu1Hz3wGHJ3k0vTcB5wNrus/jfxW4crq6+/xlVf1jd94X0TdDMcHT6b0J+YOqeqCqbgD+HDit2/4y4I1VdX31XF1V35tB/ccAq+m9IXigqv4G+CRwet8+f1VVX+leh4/0tfFHwCOAJwHpXqvbZn7q0vAMZq0Ep1TVI+iNqp7Ez6aEDwV+vZvuvDvJ3cCzgIP7jr2l7/nN9EZjBwGP7ZYBqKqfdPuumeLYyeyp/scCd3VvCPrrn4m3A8cnecok9Z3XV9e/0BsVrulCdAu9EH42vSD+EvBMZh/Mt/c9/1d6ITmZQ4HHTjj/NwAj3fa19Ebfs/VY4Jbud7LbzTz4dzNpG7sQfy/wPuDOJJuSPHKANkgDM5i1YlTVlfSmMf+4W3ULvRHxAX2P/arq3L7D1vY9fxy9EdV3gX+mFyxA77PJbt8d/VVObMKE5T3VfxvwqCT7Tah/Juf5PeDdwFsnqe+/TqhvVVV9qdt+Jb1p36cBX+2Wj6c35fsF5t4twI0T2vOIqjqxb/svDlDuPwNrk/T///Y4Hvy7mVJVvaeqjgKOoDcl/5oB2iANzGDWSvNu4LndaPLDwPOTHN9daLRvdxHRIX37/0aSI5I8nN708MVV9WN6U7QnJTkuyUOBs4Ef0htlTuUOep+j7jZl/VV1M70R7FuS7JPkWcDzZ3Ge76Q3rf7kvnV/Brw+yS8BJNk/ya/3bb8SeClwbVU9QO/z4ZfRC8/vzKLumfoKcG93kdyq7jU4MsnTu+0fAN6a5PDuoqx/1021w8+/lv2uojcKfm13IdsYvddu83QNSvL0JM/ofqf3AT8AfjLNYdKcMpi1onQB80Hg96rqFnoXRL0B+A69EdprePC/iw/RG2XfDuwLvLIr53p6nxH/Cb0R9PPp/VnWA3uo/m3AG7tp21fPoP7/DDyD3pTzm7t2z/Q87wH+F3Bg37q/ojfNvTnJPcA3gOf1HfYlep+T7x4dX0svmOZjtEz3BufX6H2+eyO91/EDwP7dLu+k9wboc/QuaDu/ax/A7wMXdq/lqRPKfYDe7+N5XZnvB15aVd+cQbMeSe9z7rvoTX9/D/ijwc5QGkyqJs6uSYLen/7Quwr6A4vdFkkrhyNmSZIaYjBLktQQp7IlSWqII2ZJkhqyUDfX36ODDjqo1q1btyB13Xfffey3337T7yhNwT6kYdh/BLB169bvVtVjJtvWRDCvW7eOLVu2LEhd4+PjjI2NLUhdWp7sQxqG/UcASaa8k59T2ZIkNcRgliSpIQazJEkNmTaYk/xF96Xh3+hbd2CSy5N8q/v5qG59krwnyfYk/5Dkl+ez8ZIkLTczGTFfAJwwYd05wBVVdThwRbcMvXvTHt49NgJ/OjfNlCRpZZg2mLsvaP+XCatPBi7snl8InNK3/oPdl5p/GTggycFIkqQZGfQz5pGquq17fjs/+2LzNTz4y+Fv5cFfTi5JkvZg6L9jrqpKMuv7eibZSG+6m5GREcbHx4dtyozs3LlzwerS8mQf0jDsP5rOoMF8R5KDq+q2bqr6zm79DmBt336HdOt+TlVtAjYBjI6O1kL9wb1/3K9h2Yc0DPuPpjNoMF8KnAGc2/28pG/9K5JspvcF79/vm/KWpBVj3TmXTbr+7PW7OHOKbTede9J8NklLxLTBnOSjwBhwUJJbgTfTC+SLkpwF3Ayc2u3+KeBEYDvwr8BvzkObJUlatqYN5qo6fYpNx02ybwEvH7ZRkiStVN75S5KkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1ZO/FboAkafGsO+eyWe1/07knzVNLtJsjZkmSGuKIWZKWCEe3K4MjZkmSGmIwS5LUEINZkqSGGMySJDVkqGBO8j+SXJPkG0k+mmTfJIcluSrJ9iQfS7LPXDVWkqTlbuBgTrIGeCUwWlVHAnsBpwFvB95VVY8H7gLOmouGSpK0Egw7lb03sCrJ3sDDgduAY4GLu+0XAqcMWYckSStGqmrwg5NXAX8I3A98DngV8OVutEyStcCnuxH1xGM3AhsBRkZGjtq8efPA7ZiNnTt3snr16gWpS8uTfUgzsW3H9yddP7IK7rh/8mPWr9l/oDKnMl1581Wmprdhw4atVTU62baBbzCS5FHAycBhwN3A/wFOmOnxVbUJ2AQwOjpaY2NjgzZlVsbHx1mourQ82Yc0E2dOcTOQs9fv4h3bJv+v96YXjw1U5lSmK2++ytRwhpnKfg5wY1V9p6p+BHwCeCZwQDe1DXAIsGPINkqStGIME8z/BByT5OFJAhwHXAt8Hnhht88ZwCXDNVGSpJVj4GCuqqvoXeT198C2rqxNwOuA302yHXg0cP4ctFOSpBVhqC+xqKo3A2+esPoG4OhhypUkaaXyzl+SJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1ZKhgTnJAkouTfDPJdUn+fZIDk1ye5Fvdz0fNVWMlSVruhh0xnwd8pqqeBDwFuA44B7iiqg4HruiWJUnSDAwczEn2B54NnA9QVQ9U1d3AycCF3W4XAqcM20hJklaKVNVgByZPBTYB19IbLW8FXgXsqKoDun0C3LV7ecLxG4GNACMjI0dt3rx5oHbM1s6dO1m9evWC1KXlyT6kmdi24/uTrh9ZBXfcP/kx69fsP1CZU5muvPkqU9PbsGHD1qoanWzbMME8CnwZeGZVXZXkPOAe4Hf6gzjJXVW1x8+ZR0dHa8uWLQO1Y7bGx8cZGxtbkLq0PNmHNBPrzrls0vVnr9/FO7btPem2m849aaAypzJdefNVpqaXZMpgHuYz5luBW6vqqm75YuCXgTuSHNxVfDBw5xB1SJK0ogwczFV1O3BLkid2q46jN619KXBGt+4M4JKhWihJ0goy+XzKzP0O8JEk+wA3AL9JL+wvSnIWcDNw6pB1SJK0YgwVzFX1dWCyOfLjhilXkqSVyjt/SZLUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkNMZglSWqIwSxJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1JChgznJXkm+luST3fJhSa5Ksj3Jx5LsM3wzJUlaGeZixPwq4Lq+5bcD76qqxwN3AWfNQR2SJK0IQwVzkkOAk4APdMsBjgUu7na5EDhlmDokSVpJUlWDH5xcDLwNeATwauBM4MvdaJkka4FPV9WRkxy7EdgIMDIyctTmzZsHbsds7Ny5k9WrVy9IXVqe7EOaiW07vj/p+pFVcMf9kx+zfs3+A5U5lenKm68yNb0NGzZsrarRybbtPWihSX4NuLOqtiYZm+3xVbUJ2AQwOjpaY2OzLmIg4+PjLFRdWp7sQ5qJM8+5bNL1Z6/fxTu2Tf5f700vHhuozKlMV958lanhDBzMwDOBFyQ5EdgXeCRwHnBAkr2rahdwCLBj+GZKkrQyDPwZc1W9vqoOqap1wGnA31TVi4HPAy/sdjsDuGToVkqStELMx98xvw743STbgUcD589DHZIkLUvDTGX/VFWNA+Pd8xuAo+eiXEmSVhrv/CVJUkMMZkmSGmIwS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BCDWZKkhhjMkiQ1xGCWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIaYjBLktQQg1mSpIYYzJIkNcRgliSpIQazJEkN2XuxGyBJ0p6sO+eyWe1/07knzVNLFoYjZkmSGmIwS5LUEINZkqSGGMySJDXEi78kSXNmthdqwdK/WGuuOWKWJKkhBrMkSQ0xmCVJaojBLElSQwYO5iRrk3w+ybVJrknyqm79gUkuT/Kt7uej5q65kiQtb8OMmHcBZ1fVEcAxwMuTHAGcA1xRVYcDV3TLkiRpBgYO5qq6rar+vnt+L3AdsAY4Gbiw2+1C4JRhGylJ0kqRqhq+kGQd8AXgSOCfquqAbn2Au3YvTzhmI7ARYGRk5KjNmzcP3Y6Z2LlzJ6tXr16QurQ82Yc0E9t2fH/S9SOr4I77Jz9m/Zr9BypzKtOVNx9lzra8+ShzJue92DZs2LC1qkYn2zZ0MCdZDVwJ/GFVfSLJ3f1BnOSuqtrj58yjo6O1ZcuWodoxU+Pj44yNjS1IXVqe7EOaialutHH2+l28Y9vk93aa7kYb8/EtS3Nd5nzcYGQ5frtUkimDeag7fyV5KPBx4CNV9Ylu9R1JDq6q25IcDNw5TB2SJM21lsN+mKuyA5wPXFdV7+zbdClwRvf8DOCSwZsnSdLKMsyI+ZnAS4BtSb7erXsDcC5wUZKzgJuBU4droiRJK8fAwVxVXwQyxebjBi1XkqSVzDt/SZLUEINZkqSG+H3MkkTbV+lqZXHELElSQwxmSZIaYjBLktQQg1mSpIZ48Zekn2rhPsczKVNazhwxS5LUEINZkqSGGMySJDXEYJYkqSEGsyRJDTGYJUlqiMEsSVJDDGZJkhpiMEuS1BDv/CUtYX5VobT8OGKWJKkhBrMkSQ0xmCVJaojBLElSQwxmSZIa4lXZ0iT2dLXz2et3ceYk273iWdJccMQsSVJDHDFLC8S/OZY0E46YJUlqiCNmSUuOsw9azhwxS5LUEINZkqSGGMySJDXEYJYkqSFe/KVlwYuBJC0X8zJiTnJCkuuTbE9yznzUIUnScjTnI+YkewHvA54L3Ap8NcmlVXXtXNelpcnRrSRNbT5GzEcD26vqhqp6ANgMnDwP9UiStOykqua2wOSFwAlV9bJu+SXAM6rqFRP22whs7BafCFw/pw2Z2kHAdxeoLi1P9iENw/4jgEOr6jGTbVi0i7+qahOwaaHrTbKlqkYXul4tH/YhDcP+o+nMx1T2DmBt3/Ih3TpJkjSN+QjmrwKHJzksyT7AacCl81CPJEnLzpxPZVfVriSvAD4L7AX8RVVdM9f1DGHBp8+17NiHNAz7j/Zozi/+kiRJg/OWnJIkNcRgliSpIUsumJOckqSSPGmOy319dwvR65Mc363bN8lXklyd5Jokb5nLOrU45qMPJXl0ks8n2ZnkvRO2HZVkW9e/3pMkc1WvFsdC9qEkj0jy9b7Hd5O8e67qVXuWXDADpwNf7H7OiSRH0Lt6/JeAE4D3d7cW/SFwbFU9BXgqcEKSY+aqXi2aOe9DwA+ANwGvnmTbnwK/BRzePU6Yw3q1OBasD1XVvVX11N0P4GbgE3NYrxqzpII5yWrgWcBZ9IJ09/qxJJ/sW35vkjO75ycm+WaSrd1o5ZMTy6V3y9DNVfXDqroR2A4cXT07u30e2j28Wm4Jm68+VFX3VdUX6f3n2l/fwcAjq+rL1bvS8oPAKfNxbloYC92HJtT9BODfAH87Zyek5iypYKYXoJ+pqn8EvpfkqD3tnGRf4H8Dz6uqo4BJb38GrAFu6Vu+tVtHkr2SfB24E7i8qq4a8hy0uOarD01lDb3+tNtP+5aWrIXuQ/1OAz5W/jnNsrbUgvl0el+KQfdzummkJwE3dKNggI/OtsKq+nE3fXQIcHSSI2dbhpqy4H1Iy85i9qHThjxeS8Ci3St7tpIcCBwLrE9S9G5eUkleA+ziwW8y9p1l8dPeRrSq7k7yeXqfD35jluWrAfPch6ayg15/2s1b1C5hi9SHdtf9FGDvqto6l+WqPUtpxPxC4ENVdWhVrauqtcCNwK/QuxjiiCQPS3IAcFx3zPXALyRZ1y2/aIqyLwVO644/jN4FOl9J8piuPJKsovcd09+ch3PTwpjPPjSpqroNuCfJMd3V2C8FLhn+VLRIFrwP9TkdR8srwpIZMdPrlG+fsO7jwOlV9d+TXERvJHsj8DWAqro/yW8Dn0lyH737eP+cqrqmO/5aeu96X15VP+4u3Lmwu0L7IcBFVTXZxWNaGuatDwEkuQl4JLBPklOA/1hV1wK/DVwArAI+3T20NC1WHwI4FThxLk9GbVr2t+RMsrqqdnajlfcB36qqdy12u7R02Ic0LPuQZmMpTWUP6re6q6qvAfand3WkNBv2IQ3LPqQZW/YjZkmSlpKVMGKWJGnJMJglSWqIwSxJUkMMZkmSGmIwS5LUkP8PObpP5Ew5P1sAAAAASUVORK5CYII=\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\", 0],\n",
    "    [\"2020-08-07\", 0],\n",
    "    [\"2020-08-08\", 1],\n",
    "    [\"2020-08-09\", 2],\n",
    "    [\"2020-08-10\", 6],\n",
    "    [\"2020-08-11\", 5],\n",
    "    [\"2020-08-12\", 9],\n",
    "    [\"2020-08-13\", 14],\n",
    "    [\"2020-08-14\", 14],\n",
    "    [\"2020-08-15\", 4],\n",
    "    [\"2020-08-16\", 16],\n",
    "    [\"2020-08-17\", 99],\n",
    "    [\"2020-08-18\", 84],\n",
    "    [\"2020-08-19\", 85],\n",
    "    [\"2020-08-20\", 24],\n",
    "    [\"2020-08-21\", 26],\n",
    "    [\"2020-08-22\", 19],\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.1.2 Fitting an SIR model to campus data](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.2-Fitting-an-SIR-model-to-campus-data)",
     "section": "2.1.2 Fitting an SIR model to campus data"
    }
   },
   "source": [
    "## 2.1.2 Fitting an SIR 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 a standard SIR model for infectious disease in a homogeneous population.\n",
    "\n",
    "$$\\text{Susceptible}\n",
    "\\xrightarrow {\\frac{\\beta S I}{N}} \n",
    "\\text{Infectious} \n",
    "\\xrightarrow{\\gamma I} \n",
    "\\text{Recovered} $$\n",
    "\n",
    "The resulting differential equations are:\n",
    "\n",
    "$$\\begin{align*}\n",
    "\\frac{dS}{dt} &= -\\beta S \\frac{I}{N} \\\\\n",
    "\\frac{dI}{dt} &= \\beta S \\frac{I}{N} - \\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",
    "| $\\gamma$ | average recovery rate | 1/8.0 $d^{-1}$ | 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 | tbd | calculated as ${\\beta}/{\\gamma}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 312
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 11775,
     "status": "ok",
     "timestamp": 1597709340495,
     "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.1.2 Fitting an SIR model to campus data](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.2-Fitting-an-SIR-model-to-campus-data)",
     "section": "2.1.2 Fitting an SIR model to campus data"
    },
    "outputId": "de914a60-eab6-40aa-88b6-90c44cdde82c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R0 = 2.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxV9Z3/8dcnJBBCEsK+QwAFBAFBUCuLVqsVbV1a29q6oLWlLTPz6zh0WmfsVDutXcZap639tVV0SgVXRKUVfq371lEEREAQSVgDgYQEQkJYQvL9/XFOYhJukpvcm7uc+34+HveRe89y74fD5c033/M932POOUREJFjS4l2AiIhEn8JdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXZKCme0ws++Y2XozqzCzJ8wss9H6z5jZOjM7ZGZ/N7NJ/vJbzOzPjbbbamZPNXq928zOauEzZ/rvdcjf7mZ/+RVm9p6ZHfaX39Von0wzW2xmZf5+75rZAH9dTzN7yMyKzWyPmf3YzLr4604zs9f8P9sBM3siukdQUo3CXZLJF4HLgJHAJOBmADObAjwMfAPoA/wBWG5m3YDXgFlmlmZmg4GuwCf8/UYB2cD65h9kZiOAlcBvgH7AWcA6f/UR4CYgD7gC+JaZXe2vmwv0BIb5tXwTOOqv+yNwEjgNmAJcCnzNX/cj4G9AL2Co/7kiHaZwl2Tya+fcXudcOfBnvMAFmAf8wTn3jnOu1jm3CDgOnOec2wZU+tvOBv4K7DWzccAFwBvOuboQn/UV4EXn3GPOuRrnXJlzbh2Ac+5V59wG51ydc2498Jj/XgA1eKF+ml/LGufcYb/1fjnwz865I865EuA+4LpG+40ABjvnjjnn3ozaUZOUpHCXZLKv0fNqvFY3eKG4wO8GOWRmh/BazoP99a8BF+KF+2vAq3hhfIH/OpRhQGGoFWZ2rpm9YmalZlaB1zrv669+BO8/kMfNbK+Z/ZeZZfg1ZgDFjWr8A9Df3++7gAGrzOwDM/tqWEdEpAUKdwmC3cDdzrm8Ro8s59xj/vr6cJ/lP3+NtsN9NzC6hXWPAsuBYc65nsDv8YIZv5X/Q+fceOB84DN4XTi78X6b6Nuoxlzn3AR/v33Oua875wbjdS/9XzM7rcNHRFKewl2C4EHgm36L2sysh3/SM8df/xrwSaC7c64IeAOv774P8F4L77kE+JSZfdHM0s2sT6MTrzlAuXPumJmdg9eFA4CZfdLMJvonSg/jdbfUOeeK8frU7zWzXP8cwGgzu8Df7wtmNtR/m4OAA0J1F4mEReEuSc85txr4OnA/XjAW4J9s9dd/BFThhTrOucPANuAt51xtC++5C6+PfAFQjncydbK/ej7wn2ZWCfwAeLLRrgOBpXjBvhnvP5ZH/HU34Z3Q3eTXuRQY5K+bDrxjZlV4vxV82z9fINIhppt1iIgEj1ruIiIBpHAXEQkghbuISAAp3EVEAig93gUA9O3b1+Xn58e7DBGRpLJmzZoDzrl+odYlRLjn5+ezevXqeJchIpJUzGxnS+vULSMiEkAKdxGRAFK4i4gEUJt97mb2MN7kRyXOuTP9Zb2BJ4B8YAfwRefcQTMz4Fd4l21XAzc759Z2pLCamhqKioo4duxYR3YPjMzMTIYOHUpGRka8SxH5WGEh3HsvLF4MVVWQnQ033AALFsDoluZbk1hqc/oBM5uNNy/HnxqF+3/hTZz0MzO7HejlnPuemV0O/BNeuJ8L/Mo5d25bRUybNs01P6G6fft2cnJy6NOnD97/GanHOUdZWRmVlZWMHDky3uWIeFauhGuvhZoa71EvI8N7LF0Kc+bEr74UYmZrnHPTQq1rs1vGOfc63sRJjV0FLPKfLwKubrT8T87zNpBnZoPogGPHjqV0sAOYGX369En5314kgRQWesFeXd002MF7XV3trS8MORW+xFBH+9wH+FOYgncDhQH+8yF481bXK/KXncLM5pnZajNbXVpaGvJDUjnY6+kYSEK5995TQ725mhq4777Y1CMtiviEqvP6ddo9taRz7gHn3DTn3LR+/UKOwQ9fYSHMnw+5uZCW5v2cP1+tB5FoW7w4vHB/5JHWt5FO19Fw31/f3eL/LPGX78G7PVm9of6yzrNyJUyaBAsXQmUlOOf9XLjQW75yZdQ+6q677uIXv/hFi+ufffZZNm3aFLXPE0k4VVXR3U46TUfDfTneXd7xfz7XaPlN/t1wzgMqGnXfRF+C9f8p3CXwsrPb3qY920mnaTPczewx4H+BsWZWZGa3Aj8DLjGzrcCn/NcAK/DucFOAd+uz+Z1Sdb0Y9P/dfffdjBkzhpkzZ7JlyxYAHnzwQaZPn87kyZP5/Oc/T3V1NX//+99Zvnw5//qv/8pZZ51FYWFhyO1EktoNN3gjYlqTkQE33hibeqRlzrm4P84++2zX3KZNm05ZdoqcHOe8jpjWH7m5bb9XCKtXr3ZnnnmmO3LkiKuoqHCjR49299xzjztw4EDDNnfccYf79a9/7Zxzbu7cue6pp55qWNfSdu0V1rEQiYWCAueyslr/95aV5W0nnQ5Y7VrI1eS+QrWT+//eeOMNrrnmGrKyssjNzeXKK68EYOPGjcyaNYuJEyeyZMkSPvjgg5D7h7udSNIYPdobx56VdWoLPiPDW750qS5kSgDJHe5x6v+7+eabuf/++9mwYQN33nlni+PQw91OJKnMmQPr18O8eU1HqM2b5y3XBUwJIbnDvZP7/2bPns2zzz7L0aNHqays5M9//jMAlZWVDBo0iJqaGpYsWdKwfU5ODpWVlQ2vW9pOJOmNHg333w8VFVBb6/28/3612BNIcof7ggXhhfttt3Xo7adOncqXvvQlJk+ezJw5c5g+fToAP/rRjzj33HOZMWMG48aNa9j+uuuu45577mHKlCkUFha2uJ2ISGdrc26ZWAg1t8zmzZs544wz2t45Bea5CPtYiEhKiWhumYSn/j8RkVMkxG32Ilbf/3f//fGuREQkISR/y11ERE6hcBcRCSCFu4hIACncRUQCSOEeI/n5+Rw4cCDibUREwqFwFxEJIIV7K3bs2MG4ceO4+eabGTNmDNdffz0vvvgiM2bM4PTTT2fVqlWUl5dz9dVXM2nSJM477zzWr18PQFlZGZdeeikTJkzga1/7Go0vFlu8eDHnnHMOZ511Ft/4xjeora2N1x9RRAIqKca559/+fKe9946fXdHq+oKCAp566ikefvhhpk+fzqOPPsqbb77J8uXL+clPfsKwYcOYMmUKzz77LC+//DI33XQT69at44c//CEzZ87kBz/4Ac8//zwPPfQQ4F1t+sQTT/DWW2+RkZHB/PnzWbJkCTfddFOn/RlFJPUkRbjH08iRI5k4cSIAEyZM4OKLL8bMmDhxIjt27GDnzp08/fTTAFx00UWUlZVx+PBhXn/9dZYtWwbAFVdcQa9evQB46aWXWLNmTcM8NUePHqV///5x+JOJSJAp3NvQrVu3hudpaWkNr9PS0jh58iQZbU1c1oxzjrlz5/LTn/40qnWKiDSWFOHeVtdJPM2aNYslS5bwH//xH7z66qv07duX3NxcZs+ezaOPPsr3v/99Vq5cycGDBwG4+OKLueqqq7jtttvo378/5eXlVFZWMmLEiDj/SUQkSJIi3BPZXXfdxVe/+lUmTZpEVlYWixYtAuDOO+/ky1/+MhMmTOD8889n+PDhAIwfP54f//jHXHrppdTV1ZGRkcFvf/tbhbuIRFXyT/mbAnQsRJopLIR774XFi73baGZnezfvWbAgpW4YEuwpf0UktaxcCZMmwcKFUFnp3Za7stJ7PWmSt14U7iISB4WFMH9+03swzJ/vLW9rv2uvherqpjfnAe91dbW3vq33SQEJHe6J0GUUbzoGEjiRtLzvvffUUG+upgbuuy+6NSehhA33zMxMysrKUjrcnHOUlZWRmZkZ71JEoiPSlvfixeGF+yOPRKfeJJawo2WGDh1KUVERpaWl8S4lrjIzMxk6dGi8yxCJjva0vEPdWa2qKrzPCXe7AEvY0TIiEkC5uV4XTDjbVVREf/+A0WgZEUkMkba8b7gB2roqPCMDbryxfXUFkMJdRGInOzuy7RYsCC/cb7utfXUFkMJdRGIn0pb36NGwdClkZZ36PhkZ3vKlS1PqQqaWKNxFJHai0fKeMwfWr4d585qOk583z1s+Z050a05SOqEqIrG1cqU33LGmpunImYwM77F0qQI6TJ12QtXMbjOzD8xso5k9ZmaZZjbSzN4xswIze8LMukbyGSISMGp5x0SHW+5mNgR4ExjvnDtqZk8CK4DLgWXOucfN7PfA+86537X2Xmq5i4i0X2cOhUwHuptZOpAFFAMXAUv99YuAqyP8DBERaacOh7tzbg/wC2AXXqhXAGuAQ865k/5mRcCQUPub2TwzW21mq1P9KlQRkWjrcLibWS/gKmAkMBjoAVwW7v7OuQecc9Occ9P69evX0TJERCSESLplPgVsd86VOudqgGXADCDP76YBGArsibBGERFpp0jCfRdwnpllmZkBFwObgFeAa/1t5gLPRVaiiIi0VyR97u/gnThdC2zw3+sB4HvAv5hZAdAHeCgKdYqISDtENOWvc+5O4M5mi7cB50TyviIiEhlNPyAiEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIAEUU7maWZ2ZLzexDM9tsZp8ws95m9oKZbfV/9opWsSIiEp5IW+6/Av6fc24cMBnYDNwOvOScOx14yX8tIiIx1OFwN7OewGzgIQDn3Ann3CHgKmCRv9ki4OpIixQRkfaJpOU+EigF/sfM3jOzhWbWAxjgnCv2t9kHDAi1s5nNM7PVZra6tLQ0gjJERKS5SMI9HZgK/M45NwU4QrMuGOecA1yonZ1zDzjnpjnnpvXr1y+CMkREpLlIwr0IKHLOveO/XooX9vvNbBCA/7MkshJFRKS9Ohzuzrl9wG4zG+svuhjYBCwH5vrL5gLPRVShiIi0W3qE+/8TsMTMugLbgFvw/sN40sxuBXYCX4zwM0REpJ0iCnfn3DpgWohVF0fyviIiEhldoSoiEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRiZO6OkdtXciL+COmcBcRiYO9h45y08OrePCNbZ3y/pFexCQiIu3gnOPptXv44fIPqDx+klXby7lwbD/GDcyN6uco3EVEYqS08jj//swGXti0v2FZTV0dq7aXK9xFRJLRyg3F3PHsRsqPnGhYNqJPFvd+YTLT8ntH/fMU7iIinaiiuoYfLN/Ic+v2Nll+43kjuH3OOHp065wYVriLiHSSV7eU8L2n17P/8PGGZYN6ZvLzz09i9pjOvY+Fwl1EJMqqjp/k7uc389iqXU2Wf27qEO787AR6ds/o9BoU7iIiUfTOtjK+s/R9dpcfbVjWN7srd18zkU9PGBizOhTuIiJRcKymll/8dQsPvbUd1+i6pMsmDOTua86kT3a3mNajcBcRidD6okP8y5PvU1BS1bAsNzOd/7zqTK46azBmFvOaFO4iIh104mQd979SwG9fKWgyjcAFY/rx889PYmDPzLjVpnAXEemAzcWH+c5T7/PB3sMNy7K6duH7V4zny+cMi0trvTGFu4hIOxw5fpJfvbSVh97c3qS1fk5+b37xhckM75MVx+o+pnAXEQnT3z7Yx13LP2BvxbGGZV3T0/jup8fy1RkjSUuLb2u9MYW7iEgbig5Wc9fyTby4eX+T5eeN6s2Pr57Iaf2z41RZyxTuIiItqKmt4+E3t/PfL27laE1tw/I+PbpyxxVncM2UIXHvW2+Jwl1EJITVO8q545mNbNlf2WT5l88ZxvcuG0deVtc4VRYehbuISCOHqk/ws5Uf8vi7u5ssHzcwh7uvOZOzR0R/BsfOoHAXEcG7icaytXu4e8XmJtPyds/owm2XnM4tM0aS0SV5bl6ncBeRlFdQUskdz2zkne3lTZZfMn4Ad105gSF53eNUWccp3EUkZR2rqeX+lwv4w+uF1NR+PGZ9cM9M7rpyApfGcKKvaFO4i0jKcc7x4uYSfvSXTewqr25Y3iXNuHXmSL598emddhONWEnu6kVE2mnjngp+/Pwm3t7WtAtm6vA87r5mImcMiu69TONF4S4iKaG44ij3/HULz7y3p8mUvD27Z3D7nHF8adqwhLrCNFIKdxEJtKrjJ/nDa4U8+MY2jtXUNSzvkmZcf+5wvn3x6TGfaz0WIg53M+sCrAb2OOc+Y2YjgceBPsAa4Ebn3InW3kNEJNpO1tbx5OoifvnCRxyoOt5k3afO6M/tc85IyGkDoiUaLfdvA5uB+o6qnwP3OeceN7PfA7cCv4vC54iIhOXVLSX8ZMVmPtpf1WT5hMG53HHFGZw/um+cKoudiMLdzIYCVwB3A/9i3iQLFwFf8TdZBNyFwl1EYuDDfYe5+/nNvLH1QJPlA3Mz+c6nx/K5KUMC1a/emkhb7v8NfBfI8V/3AQ455076r4uAIaF2NLN5wDyA4cOHR1iGiKSyksPH+OULH/Hk6t00mmKdrK5d+NYFo/narFF079olfgXGQYfD3cw+A5Q459aY2YXt3d859wDwAMC0adNcG5uLiJyi+sRJHnx9O394vZDqEx/P2phm8KXpw7jtkjH0z4nfre7iKZKW+wzgSjO7HMjE63P/FZBnZul+630osCfyMkVEPlZTW8fSNUX894sfsf9w05Ols8f0447Lz2DswJwW9k4NHQ5359y/Af8G4Lfcv+Ocu97MngKuxRsxMxd4Lgp1iohQU1vHsrVF/OblAooOHm2ybuyAHP79ijO4YEy/OFWXWDpjnPv3gMfN7MfAe8BDnfAZIpJCamrreGbtHn7zylZ2lzcN9X453VhwyRi+MG0YXVLkZGk4ohLuzrlXgVf959uAc6LxviKS2k7W1vHMe3u4/5UCdpZVN1nXKyuDebNHc9MnRiT9PDCdQUdERBLOydo6nlu3l9+8vJUdzUI9LyuDr88axdzz88lWqLdIR0ZEEkZtneO5dXv4zcsFbD9wpMm6nt0zmDdboR4uHSERibvaOsef39/Lr1/ayrZmoZ6bmc7XZ43i5hn55GRmxKnC5KNwF5G4qa1z/GW9F+qFpU1DPSczna/NHMUtM/PJVai3m8JdRGKupraOv6zfy29fKaSgpOn8LzmZ6dw6cyS3zBhJz+4K9Y5SuItIzBw+VsPjq3bx8Js72Hf4WJN1Od3SuWXmSG6dqVCPBoW7iHS6vYeO8j9vbeexVbupOn6yybrsbuncMiOfW2eOJC+ra5wqDB6Fu4h0mo17Klj4xjb+sr6Yk3VNp5Dqm92Nm88fwQ3njVCodwKFu4hElXOO1z4q5cE3tvFWQdkp60/rn83XZ43kqrOGkJmRWjM1xpLCXUSi4vjJWpav28vCN7azZX/lKevPG9WbebNHceGY/ikzp3o8KdxFJCIV1TUsWbWTP761g5LKpjM0dkkzLp84iK/PGsmkoXlxqjA1KdxFpEN2lh1h0d938vi7u5rMpQ7eTTKumz6cW2bkM6x3VpwqTG0KdxEJ28naOl76sITFb+885VZ2AP1zunHLjJF85Zzh9MzScMZ4UriLSJuKK47y+KrdPP7urlNujgHeXOpfnz2KKycPpmt6WhwqlOYU7iISUl2d442CAyx5eycvfVhCbbOhjGbwybH9uekTI7hgTD/MdJI0kSjcRaSJsqrjPLm6iMdW7WJXefUp6/tmd+O66cO47pxhDO2l/vREpXAXEZxzrNpezpJ3drFyYzE1tafes/780X24/twRXDJ+gLpekoDCXSSFVRyt4Zm1RSx5Zxdbm03gBd4c6teePZSvnDuc0f2y41BhgioshHvvhcWLoaoKsrPhhhtgwQIYPTre1QEKd5GUU1vn+N/CMpatLWLFxmKO1dSdss2U4Xlcf+4IPjNpkK4ibW7lSrj2Wqip8R4AlZWwcCEsWgRLl8KcOfGtEYW7SMrYsq+SZe8V8ex7e0KOeMnq2oWrpwzh+nOHM2FwzzhUmAQKC71grz71XERD2F97LaxfH/cWvMJdJMBKK4+z/P29LFtbxAd7D4fcZtzAHK4/bwRXnzVYdzpqy733ftxab0lNDdx3H9x/f2xqaoE5d+qJk1ibNm2aW716dbzLEAmEYzW1vLBpP8vWFvH61gOnDGEE6NOjK5+dPJjPTx3KmUNyU28YY0f7zHNzvS6YtuTmQkVF9OptgZmtcc5NC7lO4S6S/OrqHO/uKGfZ2j2s2FBMZbM50wG6pqdxyfgBfG7KEGaP6UdGlxQd8RKqzxwgI8N7tNZnnpYG4WRmWhrU1ra9XYRaC3d1y4gksW2lVTzz3h6Wrd3DnkNHQ25zTn5vrpk6hMsnDtIdjiLtM8/ODq/lnh3/kUUKd5Eks620ipUb97FyYzEb94TuR8/vk8Xnpg7lmilDNHFXY5H2md9wgzcqprX3yMiAG2+MrM4oULeMSBIoKKlkxYZ9rNhQzIf7Qrcce3bP4LOTB3HNlKFMHZ6Xev3o4Yi0z7ywECZNCt3yr5eVFbPRMuqWEUkyzjm27PcCfeWG4pAXGAFkdDEuHNufz08dwifH9adbusakt6oq9HEMe7vRo70++bb67BPgQiaFu0iCcM7xwd7DrNhQzMqN+9h+4EjI7bqmp3HBmH5cPnEgF58xgFwNXwxfNPrM58zxWub33QePPPLxaJsbb4TbbkuIYAeFu0hcOed4v6iClRuKWbGxmN3loU+KZmakcdG4/sw5cxCfHNef7G76p9sh0eozHz3a65OP81j21ugbIhJjx2pqeXtbGa98WMKLm0taHOXSo2sXLjpjAJefOZALxvYjq6v+uUZswQJvioC2wv2222JXUyfRt0UkBoorjvLKh6W8/OF+3ioo42hN6DHQOd3SuWT8AC47cyCzx/TTvC7RlkR95pHqcLib2TDgT8AAwAEPOOd+ZWa9gSeAfGAH8EXn3MHISxVJHrV1jnW7D/LyhyW8/GEpm4tDD1kEb5TLpeMHcPnEQZx/Wh+dFO1sidBnHoNZJTs8FNLMBgGDnHNrzSwHWANcDdwMlDvnfmZmtwO9nHPfa+29NBRSgqCiuobXtpbyyoclvLqlhIPVLf/qP6pvDz45rj8XjevPOSN7p+7Voqkokitkm+mUoZDOuWKg2H9eaWabgSHAVcCF/maLgFeBVsNdJBk559haUsVLm0t45cMS1uw6GHIeF/CGLJ47sg8X+YGe37dHjKuVhBDDWSWj0uduZvnAFOAdYIAf/AD78LptQu0zD5gHMHz48GiUIdLpiiuO8lZBGX8vOMBbhQdCTp1br19ONy4a259PjuvPzNP7aoSLxHRWyYivUDWzbOA14G7n3DIzO+Scy2u0/qBzrldr76FuGUlUFdU1/O+2A7xVUMZbhQfYVhp67Dl4N4yeNDSPi8Z6rfMJg3NJS9NVotJIlGeV7LQrVM0sA3gaWOKcW+Yv3m9mg5xzxX6/fEkknyESS8dqanl3R7nXOi88wIY9Fa1OApiTmc6s0/ty0bgBXDi2H32zu8WuWEk+kV4h2w6RjJYx4CFgs3Pul41WLQfmAj/zfz4XUYUinehkbR3r91R43SwFZazZeZATtafedq5e1/Q0puf34vzRfZlxWl/OHJxLejKeDE2Ce4AGUgxnlYyk5T4DuBHYYGbr/GX/jhfqT5rZrcBO4IuRlSgSPcdqallfVMG7O8pZvaOc1TsOhpz7vF6awcShecwY3YcZp/Xl7BG9kn/seZLcAzSQYjirpGaFlEA7eOQEa3Ye5N2dXpBvKKpotWUOcFr/7IYwP3dUn2DNgZ5gsxqmnCgff80KKSnBOUfRwaO8u6Ocd3ccZPWO8hZnU2xsUM9MZpzWlxmn9eH80X0ZkJsZg2rjJInuARpIMbxCVi13SVona+vYsr+S1TsOssrvZmltaGK90f16MD2/N9PyezNtRC9G9Mlq39zn0eivjlefd4LdAzRlFRZG5QpZ3UNVkp5zjt3lR3m/6BDv7z7E+0WH2LjncItztNRLTzPOHNKT6fm9GsK8TyQjWqJxdWEUr1BstwS7B6hERuEuSedA1XHWFx1i3e4K1vuB3trl/PWyu6UzdUQvpo/wwvysYXl07xqlE6DR6C+Nd5+3Wu6Boj53SWhHjp9kw576EK9g3e5DLU6D29zA3EzOzv84zMcNzOm8oYnR6K+Od593Et0DVCKjlrvE1KHqE2wqPsymvYfZVHyYD/YcZmtJJS1MydJEbmY6k4flMXloHpOG9mTysLzYnvyMRqs33i3neP/mIFGllrvEXH0f+abiioYg31xcGXaLvGt6GmcOzmXS0DzOGpbH5GF5jOidFd/L+aNxdWEMr1AMKYXmM091CneJ2LGaWgpKqhpCfNPew2wuPtzqxUGNmcGY/jlMHtazIczHDsxJvGlwo3F1YQyvUGxRIsxnLp1O4S5hq61z7Cw7wkf7qygoqeSj/VVs2VdJYWkVJ8PpVwG6dkljzMBsxg/K9R6DezJhcC49kmHGxGj0VydKn3cS3ANUIqM+dzlFqBD/aH8l2w4c4cTJ1q/ubKxXVgbjB9eHeC7jB/VkVL8eidciD1cQRstIoKjPXUI6WVvHrvJqtpZUsXW/F+JbS6ooLK1qV4gD5PfJOiXIB+R2a9/FQYkuGv3V6vOWGFG4B5xzjtLK42w7cITt/mNb6RG2H6hiV3k1NbXt+81tQG43xgzI4bT+2YwZkMOYAdmMHZibXDeiiOTq0Gj0V6vPW2JA3TIBcfhYDdtL/fBuCPIqtpce4ciJ9l9pOCC3G6f3z+H0AR+H+Gn9cuiZleSTaMXz6lCRKFO3TAA45zhQdYJd5dXsLq9ml//YWeYF+YGqEx163/oQb9wSP71/AEI8lBjev1Ik3hTuCeRYTW2T4K4P8t3lR9lVXt3mPCotyRkG4iYAAAtjSURBVMlMZ1TfHozql83Ivj0aHvl9eyRXd0qk4n11qEgMqVsmhqpPnGTvoWPsPXSUvYeOsufQUfYcPNoQ5CWVbc9o2JKuXdLI75vlB3c2o/r2YGQ/L8T79OgarBObHRXvq0NFokzdMjFQV+c4cOQ4ew4ebQjwPY1CfO+ho2FNfNWanG7pDO+TxfDe3mOY/3Nk3x4MzutOF92MuXXxvjpUJIYU7mE4WVvHgaoT7Dt8jH0Vx9h/+Bj7Dh9jf8UxiiuOsbfiKMWHjrV5h5+2dEkzhuR1bxLcHwd5d3p2zwhGCzxec5knwtWhIjGS0uHunKPq+EkvrCuOe4HtB3jj5weqjoc1sVVb0tOMQXmZDMnrzuC87g0/6wN8UM/M5LzZcnvE8/6diXJ1qEgMBLLP/fjJWg5UnaC08jillcc5UHW84Xlp5XFKG73u6EnKUHp2z2gI7SF5mQyuD/Fe3rK+2d1Su+skWldndrTlr6tDJWAC2ee+q6ya5zcUNw1v/2fF0cj6tkPp06MrA3IzGdgzkwG5mQzI7cbA3EwG9MxkaF53BuV1T62RJx0RjdEqkbT8dXWopJCkbbm/VXCA6xe+E/Fnd0tPawjsgaHCOzeT/rnd6JYepbv5pLJIR6tEs+Wf7FeHxuu8hSSUQN5m76P9lVx63+sh13VJM/r06Eq/nG7eI7vbx8+bvc7ulh6Mk5TJINL7d86fH16f+bx5wR6nrqtsxRfIcD98rIZfv7i1aWj7wd0rq2t8b+ogoUXactc4dZ03kCYC2eeem5nB9z8zPt5lSHtEOlpF49R1la2ELeDj7qRTFBZ6XSS5uV4XSm6u97qwsPX9Fizwwrs1GRle33co4Y4/D/I49cWLwwv3Rx6JTT2SsBTu0j4rV3rdAgsXel0kzn08WmXSJG99S+pHq2RlnRryGRne8tZGq9xwQ3j/OQR5nLp+e5EwKdxTUUdb3o1nVWzeeqyp8ZZfe23r71M/l/m8eU0/f948b3lrJwIjbfkHgX57kTAp3FNNJC3v9vT3tqb+/p0VFd6omIoK73VbJwAjbfkHgX57kTAlZ7h3tOWZ6vtH2vJOhP7eSFr+QaDfXiRczrm4P84++2wXthUrnMvKci4jwzmv3ek9MjK85StWaP+WfOtbp+7X/JGR4dw//EPo/c1a37f+kZbW+p9BIhPpd0gCA1jtWsjVTglr4DJgC1AA3N7W9mGHe0GB9+VtLViysrzttP+pcnLCC+fc3M7ZX6KnoMD7Tzg31/vPNDfXe93S370EUmvhHvVuGTPrAvwWmAOMB75sZtEZkB5pn2+q7x/pSAv19yaOjp63kJQR9StUzewTwF3OuU/7r/8NwDn305b2CfsK1Xhf4Zjq++vqSJGE0toVqp1xQnUIsLvR6yJ/WfOi5pnZajNbXVpaGt47R9ryTPX9I215a7SKSNKI22gZ59wDzrlpzrlp/fr1C2+nSMf4pvr+0RhpkeqjVUSSRGeE+x5gWKPXQ/1lkYu05Znq+0er5a3+XpHE19KZ1o4+8CYj2waMBLoC7wMTWttHo2VitH/j99FIC5GkRxyGQl4OfAQUAne0tb3GucdwfxEJjJiHe3sf7Qp35yJveab6/iISCK2Fe9LerENEJNXFeiikiIjEmcJdRCSAFO4iIgGUEH3uZlYK7Ozg7n2BA1EsJ9pUX2RUX+QSvUbV13EjnHMhrwJNiHCPhJmtbumEQiJQfZFRfZFL9BpVX+dQt4yISAAp3EVEAigI4f5AvAtog+qLjOqLXKLXqPo6QdL3uYuIyKmC0HIXEZFmFO4iIgGUNOFuZpeZ2RYzKzCz20Os72ZmT/jr3zGz/BjWNszMXjGzTWb2gZl9O8Q2F5pZhZmt8x8/iFV9/ufvMLMN/mefMpGPeX7tH7/1ZjY1hrWNbXRc1pnZYTP752bbxPz4mdnDZlZiZhsbLettZi+Y2Vb/Z68W9p3rb7PVzObGqLZ7zOxD/+/vGTPLa2HfVr8LnVzjXWa2p9Hf4+Ut7Nvqv/dOrO+JRrXtMLN1Lewbk2MYkZZmFEukB9AFb/rgUXw8R/z4ZtvMB37vP78OeCKG9Q0CpvrPc/CmO25e34XAX+J4DHcAfVtZfzmwEjDgPOCdOP5d78O7OCOuxw+YDUwFNjZa9l/A7f7z24Gfh9ivN949DXoDvfznvWJQ26VAuv/856FqC+e70Mk13gV8J4zvQKv/3jurvmbr7wV+EM9jGMkjWVru5wAFzrltzrkTwOPAVc22uQpY5D9fClxsZhaL4pxzxc65tf7zSmAzIe4bm+CuAv7kPG8DeWY2KA51XAwUOuc6esVy1DjnXgfKmy1u/D1bBFwdYtdPAy8458qdcweBF4DLOrs259zfnHMn/Zdv490FLW5aOH7hCOffe8Raq8/Pji8Cj0X7c2MlWcI9nJtuN2zjf8ErgD4xqa4RvztoCvBOiNWfMLP3zWylmU2IaWHggL+Z2RozmxdifVg3No+B62j5H1Q8j1+9Ac65Yv/5PmBAiG0S4Vh+Fe83sVDa+i50tn/0u44ebqFbKxGO3yxgv3Nuawvr430M25Qs4Z4UzCwbeBr4Z+fc4War1+J1NUwGfgM8G+PyZjrnpgJzgH8ws9kx/vw2mVlX4ErgqRCr4338TuG8388Tbiyxmd0BnASWtLBJPL8LvwNGA2cBxXhdH4noy7Teak/4f0/JEu7h3HS7YRszSwd6AmUxqc77zAy8YF/inFvWfL1z7rBzrsp/vgLIMLO+sarPObfH/1kCPIP3q29jnXdj8/DNAdY65/Y3XxHv49fI/vruKv9nSYht4nYszexm4DPA9f5/PqcI47vQaZxz+51ztc65OuDBFj47rt9FPz8+BzzR0jbxPIbhSpZwfxc43cxG+q2764DlzbZZDtSPSrgWeLmlL3e0+f1zDwGbnXO/bGGbgfXnAMzsHLxjH5P/fMysh5nl1D/HO/G2sdlmy4Gb/FEz5wEVjbofYqXF1lI8j18zjb9nc4HnQmzzV+BSM+vldztc6i/rVGZ2GfBd4ErnXHUL24TzXejMGhufx7mmhc8O5997Z/oU8KFzrijUyngfw7DF+4xuuA9C3HQb+E+8LzJAJt6v8wXAKmBUDGubiffr+Xpgnf+4HPgm8E1/m38EPsA78/82cH4M6xvlf+77fg31x69xfQb81j++G4BpMf777YEX1j0bLYvr8cP7j6YYqMHr970V7zzOS8BW4EWgt7/tNGBho32/6n8XC4BbYlRbAV5fdf13sH702GBgRWvfhRgev0f879d6vMAe1LxG//Up/95jUZ+//I/137tG28blGEby0PQDIiIBlCzdMiIi0g4KdxGRAFK4i4gEkMJdRCSAFO4iIgGkcJeU589U+J141yESTQp3EZEAUrhLSjKzO8zsIzN7ExjrL/u6mb3rT072tJllmVmOmW33p5fAzHLrX5vZ/zFvDv/1ZvZ4XP9AIs0o3CXlmNnZeJe0n4V3JeR0f9Uy59x0501OthnvisVK4FXgCn+b6/ztavDmc5/inJuEdzWtSMJQuEsqmgU845yrdt7snfXzlpxpZm+Y2QbgeqB+WuGFwC3+81uA//GfrweWmNkNeLMwiiQMhbvIx/4I/KNzbiLwQ7z5inDOvQXkm9mFQBfnXP0kUVfgzcczFXjXn01QJCEo3CUVvQ5cbWbd/dn9PusvzwGK/f7165vt8yfgUfxWu5mlAcOcc68A38ObYjo7FsWLhEMTh0lK8m9oMRdvPvZdeDcDOYI3ZW4p3p20cpxzN/vbDwS2481ieMj/D+AVvFA3YLFz7mex/nOItEThLhIGM7sWuMo5d2O8axEJh/oIRdpgZr/Bu0vU5fGuRSRcarmLiASQTqiKiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgA/X/5SHpQMQnrBQAAAABJRU5ErkJggg==\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",
    "\n",
    "def model(t, y, beta):\n",
    "    S, I, R = y\n",
    "    dSdt = -beta*S*I/N\n",
    "    dRdt = gamma*I\n",
    "    dIdt = -dSdt - dRdt\n",
    "    return np.array([dSdt, dIdt, dRdt])\n",
    "\n",
    "def solve_model(t, params):\n",
    "    beta, I_initial = params\n",
    "    IC = [N - I_initial, I_initial, 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, I, R = soln.y\n",
    "    U = beta*S*I/N\n",
    "    return S, I, R, U\n",
    "\n",
    "def residuals(df, params):\n",
    "    S, 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, 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.1.3 Fitted parameter values](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.3-Fitted-parameter-values)",
     "section": "2.1.3 Fitted parameter values"
    }
   },
   "source": [
    "## 2.1.3 Fitted parameter values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 136
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 571,
     "status": "ok",
     "timestamp": 1597709346279,
     "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.1.3 Fitted parameter values](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.3-Fitted-parameter-values)",
     "section": "2.1.3 Fitted parameter values"
    },
    "outputId": "cd8b5c3e-2dd2-4bd6-cbd8-9b58186cefc4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter           Value\n",
      "-----------  ------------\n",
      "N            15000\n",
      "I0              24.5549\n",
      "beta             0.244001\n",
      "gamma            0.125\n",
      "R0               1.95201\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.1.4 Short term predictions of newly confirmed cases](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4-Short-term-predictions-of-newly-confirmed-cases)",
     "section": "2.1.4 Short term predictions of newly confirmed cases"
    }
   },
   "source": [
    "## 2.1.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": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 2057,
     "status": "ok",
     "timestamp": 1597709362062,
     "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.1.4 Short term predictions of newly confirmed cases](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4-Short-term-predictions-of-newly-confirmed-cases)",
     "section": "2.1.4 Short term predictions of newly confirmed cases"
    },
    "outputId": "24b19ea3-e000-4bbc-bb24-e8699380f5ca"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAEICAYAAACUDtg6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxb533n+88PXACukkhK1G7tkrXYli07dhwrcpzEWZtOm+m0adOkzYzvzMTjtJNMkvreTDO96bi9M9NeZ5xJb7olue2N0yZtkjb7YtpJvMq2ZEuyrdXWQkmUuIj7AuK5f/wAASTBHSRI8Pt+vc4LOOcAOA8AQvriwe88j4UQEBERERGR6YvkuwEiIiIiIoVC4VpEREREJEcUrkVEREREckThWkREREQkRxSuRURERERyROFaRERERCRHFK5F5iEzW2dmwcyKZ/m4nzazv5mlYzWY2b+ejWNlOfbV52lma82s08yKpvA495vZX+S+hVNjZvVm9piZdZjZ/5jF4/47M7uYfB1rk5cbZuv4w9oya3/DIrIwzep/zCKSZmb3Ah8EdgFfCSF8MK8NkqxCCKeByvFuZ2b7gL8JIazOuO9/ncGmTcU9wGWgOszSJAdmVgL8CXBrCOFgcvO4r6eIyHylcC2SP43AZ4C7gbI8t6VgmVlxCCGe73bMEdcAR2YrWCfVAzHg8ERubGZFIYTBmW2SiMjMUVmISJ6EEP4hhPANoHm825pZkZn9dzO7bGYngXcO2/9bZvZS8uf+k2b2v2XsO2Rm785YL0k+zu5RjvWgmZ0xs3Yze9bM7hh2k1Iz+3LyWIfNbE/GfVea2dfN7JKZnTKz+zL23WJmT5hZm5mdN7OHzKw0Y/9bzOxlM7tiZg8BNsbr8Wkz+5qZfTXZjufM7PqM/a+a2SfM7AWgy8yKzexWM3s8efyDyZ7m1O3Xm9mjycf6IVCXsW9ICY6Z1ZjZX5tZo5m1mtk3zKwC+C6wMlny0Jl8LYaUIJjZLyRfs7Zk2cu1w9r8MTN7IfkafNXMYsl9dWb2z8n7tZjZT80s67/fZvZ6M3sm+RjPmNnrk9u/CHwA+HiyfW/Oct8yM/sfZvZa8v4/M7OyqbbdzLYAryRv1mZmP0nePpjZplS7zOzzZvYdM+sC7kw+3n9KPl6Xmf2leUnLd5Pv0Y/MbEnG8af03mZ5/vvM7KyZfdTMmpJ/p7+VsT9q/jk8bV7m8mcZr8+jZvbLyeu3J5/jO5Prd5nZgVGOWWRePnQi2cZnzWxNct+on0Xzz9P+5L6LZvYnE3w9Pmj+b0SH+Wf010d7PURkikIIWrRoyeOC915/cZzb/FvgZWANUAM8AgSgOLn/ncBGPJC+EegGbkzu+zjw1YzHeg/w4hjH+g2gFv9l66PABSCW3PdpoBd4B1AEPAA8mdwXAZ4F/jNQCmwATgJ3J/ffBNyafNx1wEvA7yT31QEdwHuBEuB3gTjwr0dp46eBgYzbfww4BZQk978KHEi+XmXAKvxLzDuS7XxLcn1p8vZP4KULUWBvsi1/k9y3bthr/W3gq8CS5LHfmNy+DzibpZ2px9kCdCWPXZJ8X44DpRltfhpYmXyPXwL+bXLfA8CfJe9XAtwBWJbXpQZoBd6ffJ1/Lblem9z/ReAzY7z3nwMakq9XEfD65GsynbYPef2S2wKwKaNNV4Dbk+9NLPl4T+K93quAJuA5YHdy/0+A30/ef8rvbZbnvw//u/uD5PN8B/5ZWpLc/6fAt5LPsQr4J+CB5L4/AP5n8vr9wAngjzP2PTjKMf8T8CKwFf/8Xp/xfo31WXwCeH/yeiVedjPm6wFUAO3A1uRtVwA78v1voBYthbbkvQFatCz0hYmF65+kwkpy/a3DA8uw238D+Ejy+spkoKhOrn8N+Pgk2tcKXJ+8/mngRxn7tgM9yeuvA04Pu+/vAX89yuP+DvCPyeu/STKkJ9cNOMvY4Trz9hHgPHBHcv1V4Lcz9n8C+H+HPcb38Z7ctclAVZGx7/8jS7hOhpEEybA17PH2MXa4/hTwd8PafA7Yl9Hm38jY/38Bf5a8/gfAN0kG0jHeq/cDTw/b9gTwweT1LzJKuE62pyf1Xg/bN522X339MvYPD9dfHna8V4Ffz1j/OvD5jPX/AHxjuu/tKO9hz7C2NuFfCg3/grExY99twKnk9buAF5LXvwf8a9JfPB8FfmmUY74CvGcKn8XHgP8C1A27zVivRwXQBvwyUDaRY2rRomXyi8pCROaHlcCZjPXXMnea2dvN7MlkyUAb3mtVBxBCaAR+DvyymS0G3g78bfJ+hy1dxnBHctvHzEtMriQfaxFDf0q/kHG9G4iZl0xcg5dFtKUWvAevPvm4W5KlDRfMrB34rxmPO+T5hRDCsOebTebtE3gYX5ltf7Jt/3JY296Ah+WVQGsIoSvj9kNe3wxrgJYQQus4bctmZebjJtt8Bu9pTBn+2qZO/PtveE/xD5I/6X9yIsdIem3YMUZTh/cKn8hx2yci23t9MeN6T5b11OPn6r1NaQ5Da/RTz2UpUA48m3Gc7yW3g3+J2WJm9cANwJeBNWZWB9yCh+Fs1pD9NR/vs/gh/BeFl83Lf9413uuRfB3+Ff5L2Hkz+7aZbRvn9RCRSdIJjSLzw3n8P+GUtakrZhbFe/Z+E/hmCGHAzL7B0JrlL+E9acXAEyGEcwAhhB2ZB0kG7I/jvXCHQwgJM2sd9lijOYP34m0eZf/ngeeBXwshdJjZ7+BlHSOen5nZsOebTebtI8Bq/CTRlMyT9s7gvXn/ZviDmNk1wBIzq8gIYWuH3T/zcWrMbHEIoW3YvvFOEmzER4ZJHTf1HM+Ncz9CCB14WcBHzWwn8BMzeyaE8OMsx7hm2La1eAgcz2W85GcjcHDYvim3fYKmc4Jlrt7b8VzGQ/2O1OcnUwih28yeBT4CHAoh9JvZ48B/BE6EEC6P0f6NwKFhbR/zsxhCOAb8WvJv/5eAr5lZLWO8Hsn7fR/4frJW/DPAn+NlRiKSI+q5FskT85PsYnhta5H5CWCjfeH9O+A+M1ttfiJXZs9lKV5PegmIm9nb8bKRTN8AbsT/4//yGM2qwn9GvwQUm9l/Bqon+JSeBjrMTyQsS56otdPMbs547HagM9lb9u8y7vttYIeZ/VLyNbgPWD7O8W7KuP3vAH14nW42fwO828zuTrYrZn7y2uoQwmvAfuC/mFmpmb0BeHe2BwkhnMdPXPxfZrbE/OTQvcndF4FaM1s0Shv+Dnhn8uS2Ejws9wGPj/M8MbN3mdmmZKi9Agzi5SnDfQfvPX1f8u/rX+GlO/883jGSvdF/BfyJ+cmYRWZ2W/LL25TbPgty8t6OJ/n6/Dnwp2a2DMDMVpnZ3Rk3exS4N3kJXr+euZ7NXwD/p5ltNnddMiSP+Vk0s98ws6XJdqW+6CXGej3MTwp9j/kJuH1AJ9n/jkRkGhSuRfLn/8B7wj6Jn7jUk9yWzZ/jdZMH8RO7/iG1I9mreR8egFqB9+EnXZFxmx68d3t95n2z+D7ey3kU//m8l/HLM1LHGATehf8kfgrv6fsL/Kds8JMO34fXf/85flJg6r6XgX8J/BF+8tVmvJRlLN/Ef+JOncD3SyGEgVHadgY/kfN+PKycwU8kS/0b+D68ZrwF+H3G/gLyfvxkypfxetzfSR7jZeArwMnkz/GZJSqEEF7B3+f/ib827wbeHULoH+d5gr8eP8LD0BPA/wohPJLleTbj78FH8dfx48C7xug1He5j+Ml1z+CvxR8DkWm2fUbl+L0dzyfw8pwnk6VNP8JPREx5FA/Fj42yns2f4J/dH+BfPv8SPwl3vM/i24DDZtYJPAj8agihZ5zXI4L3pDfir8cbGfolV0RywLy0UUQKXbLna0sI4Tfy3ZbpMrNP4yfEzfvnIiIihUU11yILgJnV4CdAvT/fbRERESlkKgsRKXBm9m/wn4a/G0IY6+dpERERmSaVhYiIiIiI5Ih6rkVEREREcmRO1FzX1dWFdevW5eXYXV1dVFRU5OXYIoVInymR3NJnSiS3nn322cshhKXj33Jq5kS4XrduHfv378/LsRsaGti3b19eji1SiPSZEsktfaZEcsvMxpupdVpUFiIiIiIikiMK1yIiIiIiOaJwLSIiIiKSIwrXIiIiIiI5onAtIiIiIpIj44ZrM/srM2sys0NZ9n3UzIKZ1SXXzcw+a2bHzewFM7txJhotIiIiIjIXTaTn+ovA24ZvNLM1wFuB0xmb3w5sTi73AJ+ffhNFREREZtbTT8OFC/luhRSCccN1COExoCXLrj8FPg5kzp/+HuDLwT0JLDazFTlpqYiIiMgMaGyEAwfgW9+CH/0IQhj/PiKjmVLNtZm9BzgXQjg4bNcq4EzG+tnkNhEREZE5JwR48sn0elERmOWvPTL/TXqGRjMrB+7HS0KmzMzuwUtHqK+vp6GhYToPN2WdnZ15O7ZIIdJnSiS39JmaWWfOlHHgwGIAIpHAmjWXaGgYzHOrZD6byvTnG4H1wEHzr3argefM7BbgHLAm47ark9tGCCF8AfgCwJ49e0K+pnbVtLIiuaXPlEhu6TM1c+Jx+OpXYetWX7/xRtizZ1t+GyXz3qTLQkIIL4YQloUQ1oUQ1uGlHzeGEC4A3wJ+MzlqyK3AlRDC+dw2WURERGT6XngBurr8elkZXH99ftsjhWEiQ/F9BXgC2GpmZ83sQ2Pc/DvASeA48OfAv89JK0VERERyqLvbT2JMuflmKCnJX3ukcIxbFhJC+LVx9q/LuB6AD0+/WSIiIiIzZ/9+LwsBqKlJl4aITJdmaBQREZEFpaUFXn45vX7rrRohRHJH4VpEREQWlMyh99asgdWr89cWKTwK1yIiIrJgnDkDZ8/6dTPvtRbJJYVrERERWRASCXjiifT6tm2wZEn+2iOFSeFaREREFoSXX4a2Nr9eUgJ79uS3PVKYFK5FRESk4PX3+wghKTfc4GNbi+SawrWIiIgUvAMHoLfXr1dWwq5d+W2PFC6FaxERESlonZ3w4ovp9ZtvhuJxZ/oQmRqFaxERESloTz8Ng4N+felS2LQpv+2RwqZwLSIiIgWrqQmOH0+va8IYmWkK1yIiIlKwMieMWbcOVqzIW1NkgVC4FhERkYJ06hRcuODXIxF43evy2x5ZGBSuRUREpOAkEvDUU+n17dth0aL8tUcWDoVrERERKTiHD0N7u18vLYWbbspve2ThULgWERGRgtLXB889l16/8UaIRvPXHllYFK5FRESkoDz3nAdsgOpq2LEjv+2RhUXhWkRERApGe7uXhKS87nVQVJS/9sjCo3AtIiIiBeOpp/xkRoDly2H9+vy2RxYehWsREREpCOfP+/B7Kbfemr+2yMKlcC0iIiLzXghDJ4zZtAmWLctfe2ThGjdcm9lfmVmTmR3K2PbfzOxlM3vBzP7RzBZn7Ps9MztuZq+Y2d0z1XARERGRlBMn4NIlv15UBLfckt/2yMI1kZ7rLwJvG7bth8DOEMJ1wFHg9wDMbDvwq8CO5H3+l5npNAIRERGZMfE4PP10en3XLqiszF97ZGEbN1yHEB4DWoZt+0EIIZ5cfRJYnbz+HuDhEEJfCOEUcBzQd0cRERGZMYcOQWenX4/F4IYb8tseWdiKc/AYvw18NXl9FR62U84mt41gZvcA9wDU19fT0NCQg6ZMXmdnZ96OLVKI9JkSyS19psbW1xfhJz9ZRjxuAOzadYXHH+/Oc6tkIZtWuDaz/x2IA3872fuGEL4AfAFgz549Yd++fdNpypQ1NDSQr2OLFCJ9pkRyS5+psf30p7Bxo19fvBje+16IaLgGyaMph2sz+yDwLuCuEEJIbj4HrMm42erkNhEREZGcam2Fl19Or996q4K15N+U/gTN7G3Ax4FfCCFk/vbyLeBXzSxqZuuBzcDT2R5DREREZDqefNKH4ANYtQrWrs1ve0RgAj3XZvYVYB9QZ2Zngd/HRweJAj80M4AnQwj/NoRw2Mz+DjiCl4t8OIQwOFONFxERkYXp7Fk4cya9rgljZK4YN1yHEH4ty+a/HOP2fwj84XQaJSIiIjKa4RPGbN0KtbX5a49IJlUmiYiIyLzyyivQkhwkuLgYbr45v+0RyaRwLSIiIvPGwADs359ev+EGKC/PX3tEhlO4FhERkXnj4EHoTg6lUF4O112X3/aIDKdwLSIiIvNCVxe88EJ6/ZZbvCxEZC5RuBYREZF54ZlnIB7367W1sHlzftsjko3CtYiIiMx5ly/D0aPp9dtuAx8NWGRuUbgWERGROS9z6L1rroGVK/PXFpGxKFyLiIjInPbaa9DY6NfN4HWvy297RMaicC0iIiJzViIxtNd6+3ZYvDh/7REZj8K1iIiIzFlHjsCVK369tBRuuim/7REZj8K1iIiIzEn9/fDcc+n13bshFstfe0QmQuFaRERE5qTnnoPeXr9eVQU7d+a3PSIToXAtIiIic05HBxw6lF6/5RYoKspfe0QmSuFaRERE5pynnvKTGQGWLYONG/PbHpGJUrgWERGROeXiRTh5Mr1+2235a4vIZClci4iIyJzyxBPp6xs2QH19/toiMlkK1yIiInNRPA4HD7LohRfg4EFfXwBOnICmJr8eiWjCGJl/ivPdABEREcnQ0gIPPQSf/Sz09bErBJ+WMBqF++6De++Fmpp8t3JGDA7C00+n13fu9FFCROYThWsREZG54tgx2LsX2tqujkF39T/qzk544AH4/Ofhscdg8+a8NXOmHDrko4SAj2d94435bY/IVIxbFmJmf2VmTWZ2KGNbjZn90MyOJS+XJLebmX3WzI6b2Qtmpo+FiIjIRLS0wB13+Nl8qcGdh+vt9f179/rtC0hzM+zfn16/6SafkVFkvplIzfUXgbcN2/ZJ4MchhM3Aj5PrAG8HNieXe4DP56aZIiIiBe6hh3ye7xDGvl0I3rP9uc/NTrtmQX8//PCHXhYCXvVy7bX5bZPIVI0brkMIjwHDvx6/B/hS8vqXgF/M2P7l4J4EFpvZilw1VkREpCDF415jPVqP9XC9vfDgg+k0Os89+ii0t/v1khJ485v9ZEaR+WiqNdf1IYTzyesXgNQgOauAMxm3O5vcdp5hzOwevHeb+vp6GhoaptiU6ens7MzbsUUKkT5TIpNXcfw4u7u7J/Wfcry7m+f/+q/p2rRpxto1G06dquDQoeqr6zfe2MqBAxP8kiEyB037hMYQQjCzcX7Dynq/LwBfANizZ0/Yt2/fdJsyJQ0NDeTr2CKFSJ8pkSkoKvIC456eCd+luLSUm7du9TrteaqpCY4fh61bfX3HDrj99vy2SWS6pvqjy8VUuUfyMjkiJeeANRm3W53cJiIiIqOprp58iUci4febp3p74Uc/Sk9xvnQp3HprftskkgtTDdffAj6QvP4B4JsZ238zOWrIrcCVjPIRERERyWbHDh/HejKiUR8Ieh4KAR55xEcXBO+0f/ObvQNfZL6byFB8XwGeALaa2Vkz+xDwR8BbzOwY8ObkOsB3gJPAceDPgX8/I60WEREpJMXFPkFMLDax28difvt5mkYPHIAzGWdo3XmnJouRwjFuzXUI4ddG2XVXltsG4MPTbZSIiMiCc++9PkHMxYtjD8dnBosXw4fn53+3jY1Dx7O+/nq45pr8tUck1zTQjYiIyFxQU+MzL9bXj96DHYv5/scem5dToHd3w49/nP7usHw53HxzftskkmsK1yIiInPF5s1w+DDcfz/U1kJVFQMVFV4zUVfn2w8fnpdTnycSHqxTA6LEYnDXXRrPWgrPtIfiExERkRyqqYFPfcqD9KFDHHr0UXa/8Y1+8uI8rbEGePZZOJ8xxMFdd0FFRf7aIzJTFK5FRETmoqIiuP56rrS2emHyPHb6NDz/fHp9zx5YtSp/7RGZSfoxRkRERGZMZ6cPu5eyejXs3p2/9ojMNIVrERERmRGJhE8U09fn6xUV8KY3+YAnIoVK4VpERERmxJNP+hTn4IH6rrsmPpS3yHylcC0iIiI5d/IkHDqUXn/d63zoPZFCpxMaRUREJKeuXIFHH02vr1sH1103gTvG4z7UYHs7VFf7tPDFhR1VEgk4eBC2bNHoKYWisP9iRUREZFbF4/DDH8LAgK9XVcG+fePcqaUFHnoIPvtZL9AuKoLBQYhGfZr3e++dl5PmjOfKFT/Zs6nJhyl8+9tVj14IFK5FREQkbZq9x48/7lkZfIKYt7wFSkvHuMOxY7B3L7S1QW/v0H2dnfDAAz4t/GOPzcvJc0bz0kvwxBP+cgOcPeulNBs35rddMn0K1yIiIpKT3uOjR+Hll9Prr3+9Tyw55jHvuMO7blNzog/X2wsXL3oAP3x43vdg9/R4yczp0+ltkYhPA79hQ/7aJbmjExpFREQWumPHvIf6gQegudl7jK9c8cvmZt++Y4ffbhQtLfCzn6XXN22C7dvHOe5DD/lxRgvWKSF4z/bnPjfx5zQHvfYa/P3fDw3WS5bAv/gXPk+QSkIKg8K1iIjIQpbqPb54cWRZRkpm73Gq5iPDwICPZ50qcVi82B9yTPG495KPdsxsbXjwQe9Nn2cGBryq5fvfH/p0d+2CX/olqK3NX9sk9xSuRUREFrIc9B4/9pjvAi/PfstboKRknOMePpyeXWai+vuHju83D1y8CF//+tBymYoKeOc74bbbvPpGCovCtYiIyEKVg97jI0fgxIn0Te64w0sdxtXePvlkGYn4/eaBRAL274dvfWtokzduhPe+F1atyl/bZGbphEYREZGFajq9x9dfz6VLPjpIyrZtkxjQo7p68iUeiYTfb45ra/Mh9i5dSm8rLYU3vMFr0aWwKVyLiIgsVNPoPe7r8zrrRMI319bC7bdP4nF27PCRSDo7J36faBR27pxUc2fb4cPw1FPp+nOAlSt9rO/Kyrw1S2aRykJEREQWqmn0Hjc0QEeHbyot9TrrSeX04mIf4i8Wm9jtYzG//RwtUu7uhu9+F37+83SwjkTg1lu9vlrBeuGYVrg2s981s8NmdsjMvmJmMTNbb2ZPmdlxM/uqmY01dLyIiIjkS6r3eDKiUV5I7OS119Kb3vjGKVZr3HuvDy0y3hh0Zn67D394CgeZeadOwde+BmfOpLfV1PhIINddpyH2Fpoph2szWwXcB+wJIewEioBfBf4Y+NMQwiagFfhQLhoqIiIiOTaF3uOO37qPp/ane4937YL166d4/JoaH2qkvn70NsRivv+xx+bcBDL9/dDQ4NO9Z54Tet11Pnb1HGuuzJLploUUA2VmVgyUA+eBNwFfS+7/EvCL0zyGiIiIzJRJ9B6HRYv57oYPXx21b9kyeN3rpnn8zZu9UPn++71wu6oKFi3yy7o633748Jyb+vzCBR9i7+jR9LbKSnjXu7wUZI5Wr8gssDDeuJZj3dnsI8AfAj3AD4CPAE8me60xszXAd5M928Pvew9wD0B9ff1NDz/88JTbMR2dnZ1UqhBKJGf0mRLJrdn4TJWdPcsNH/kIxZ2dFPX3j9g/WFpKvLKSv/7gFznGFgBKSxPs3XuZsrIcTuoyOEjFqVMUd3cTLy+na/36OZdSEwl45ZUqjh8f+p6sWtXDrl1XKCmZeq6S2XHnnXc+G0LYM1OPP+XRQsxsCfAeYD3QBvw98LaJ3j+E8AXgCwB79uwJ+/btm2pTpqWhoYF8HVukEOkzJZJbs/aZesc7fIKYBx/0eodIxJNkNErRfffx4u0fJnK8hq3Jm7/97bBmzcw3ay5pbfUh9oqKYGvyhSgt9bG9N27Mb9tk7pjOUHxvBk6FEC4BmNk/ALcDi82sOIQQB1YD56bfTBEREZlRNTXwqU95GcahQz5MX3U17NzJ2fNFPP2d9E13715YwTqE9BB7mYOrrFrlQ+xVVOStaTIHTSdcnwZuNbNyvCzkLmA/8AjwXuBh4APAN6fbSBEREZklRUVw/fVXV7u64Cc/Se9euRL2zNgP6nNPVxc8+iicPZveVlTkteY7dmgkEBlpyuE6hPCUmX0NeA6IA8/jZR7fBh42s88kt/1lLhoqIiIisyuRgB//OD0SRnk5vOlNCydQnjwJP/3p0Eksa2v9NZjQFO+yIE1rhsYQwu8Dvz9s80ngluk8roiIiOTfM8/4qBjggfquuzxgF7r+fp8M5tixodtvuMF77SOagk/GoOnPRUREZITXXoODB9PrN98MK1bkrz2z5fx5P2kxc1b2ykq4886F8fxl+hSuRUREZIiODg+YKWvXDinDLkiDg7B//9AvFABbtsDrX++jgohMhMK1iIiIXBWP+4yDqeGuKyt9RIxCrrO+cAF+9jNoaUlvi0Zh795pzD4pC5bCtYiIiAAeqL/3Pbh82dcjEXjzmyc+O/p809Xlw+sdPz50++rV/oViIdSXS+4pXIuIiAi9vfCd76SDNXg5xLJl+WvTTBkchBdfhOee8576lOLi9BB7C0EI0Nzss8xL7ihci4iILHDd3R6sM8sibrsNtm/PX5tmyunT8PjjPkdOpg0b4NZbvQym0PX3w9GjPjFOezu8732aCCeXFK5FREQWsK4u+Od/hitX0tv27oVt2/LXpplw5Qo88YSH60w1Nd5Dv3Jlfto1m9raPFAfPQoDA+ntR474aDCSGwrXIiIiC1R7O3z72z46CPhJi/v2webNeW1WTg0MwPPPwwsv+KQ4KaWlHiivvbawx60OAc6c8RntM2eZTCkt9RknJXcUrkVERBagtjbvse7u9vVIxCeJKaTRMY4d8xMWU88x5dprPVgX6oma4KUfr7ySLv0YbskSry3fvBlKSma/fYVM4VpERGSBaWnxYJ2a1ryoCN7yFh/PuhBcvuwzLF68OHR7fT3cfnthn8DX1ua91EePDj1ZM+Waa2DnTli1avbbtlAoXIuIiCwgTU3w3e9CX5+vFxfD295WGDXHvb3w9NPw8stDt5eX+8mKmzblp10zbSKlH9u2+Qmq1dWz376FRuFaRERkgbhwwYN16mS20lJ4+9u9R3c+SyT8pLz9+9OT34CXulx3HezeXZilDxMp/di500s/ipX4Zo1eahERkQXg7Fn4wQ/SpQKxGLzjHfO/RKKx0YfWyxxGELzE5bbbYNGi/LRrJrW2pkf9GF76YZYu/SiEXyPmI4VrERGRAkNYRswAACAASURBVPfaaz6leWq0jPJyD9Y1Nflt13R0dsKTT8LJk0O3V1f70HqFUj+eEoIPI3joEJw7N3J/NJou/aiqmv32SZrCtYiISAE7cQIeeSQdrCsr4Z3vnL89uvG4D6t34MDI2RVvvBF27SqsoeX6+tKlH6khEzPV1Hgv9aZNKv2YK/Q2iIiIFKijR+HRR73XE7xX913vmr+zEJ465RPBdHYO3b55s09bXl6en3bNBJV+zF8K1yIiIgXoyBH42c/S64sXe7CejwG0tdXrqoeXQ9TVeQnI8uX5aVeuheAlPIcPj136sWPH/P2CtBAoXIuIiBSYF17weuSU2lovBZlvk6b09/sIIIcPp3vfwZ/HzTd70DTLX/tyRaUfhUVvkYiISAF59llfUpYt8+H2otH8tWmyQvCw+fTT6YluwIP0jh1w003z6/mM5vx5L/s4cSJ76ce6dR6qV6zIS/NkiqYVrs1sMfAXwE4gAL8NvAJ8FVgHvAr8SgihdVqtFBERkXE99RQcPJheX7HCJ4iZT2M8X7zoJSCXLg3dvnKll4DM5xFOwOvFjx71JdvY1NGoT8++fbtKP+ar6fZcPwh8L4TwXjMrBcqB+4EfhxD+yMw+CXwS+MQ0jyMiIiKjCMGn+z5yJL1t9Wp461vnTxlBczM8//zIofUqK312xQ0b8tOuXIjH/WTMo0ez11KDl+7s2KHSj0Iw5bfPzBYBe4EPAoQQ+oF+M3sPsC95sy8BDShci4iIzIgQfESQo0fT29atg7vumh9D0l2+DM89B6++OnR7URFcfz3ccMP8DZsXLqTLPlKzYmYqLfUwvWWLl+8sKAMDXmA+33+KyMJC5hkCk7mj2Q3AF4AjwPXAs8BHgHMhhMXJ2xjQmlofdv97gHsA6uvrb3r44Yen1I7p6uzspFK/u4jkjD5TIrk11mcqkYDnn19MY2PZ1W0rV/awe3cbkchstXBqWltLOHasiosXRxZPr1jRy/bt7ZSXD+ahZdPT0xPh7Nlyzpwpo6sr+7eCpUv7WLOmm+XL+ygqmloOm49iFy6w5OmnqXnmGZY89xxNd97J0Y99bNbbceeddz4bQtgzU48/nXC9B3gSuD2E8JSZPQi0A/8hM0ybWWsIYclYj7Vnz56wf//+KbVjuhoaGti3b19eji1SiPSZEsmt0T5Tg4Pwox/50G0pW7fC3r1zewSNCxe8p/rs2ZH71q+H3bvn35Ts8bi/D6+8kv15gU/as3Wrj8ldUTG77cub7m5oaIDvfx++972hP68ArFnjL9ws/8Ga2YyG6+n80HIWOBtCeCq5/jW8vvqima0IIZw3sxVA03QbKSIiImnxOPzgB0OD3M6dcNttczdYNzZ6qG5sHLlv40YP1fOtQqCpyQP1iRM+bOBwJSX+3LZuhfr62W/frAvBxxP83vc8UP/0pz7O4GjM/MzVAquJmXK4DiFcMLMzZrY1hPAKcBdeInIE+ADwR8nLb+akpSIiIkJ/v2eXCxfS2264AW65JX9tGsvZsx6qM9sLnqs2bfJQvXhE8ejc1d0Nx455qG5ry36bVas8UK9bN3/rxSespcV/Qvne9/wb32hnbIIPUL5vnw9hc/fd/iLN1W+D0zDdt/w/AH+bHCnkJPBbQAT4OzP7EPAa8CvTPIaIiIjgnYDf+c7QYer27IEbb8xfm0Zz+rSH6qZhv19HIl4asXu3T8c+HwwOevXC0aNw5szQCW1Sqqv9xMQtWwp8CL3BQR+A/Pvf9+Xpp734fzQ7dniQvvtuuOMOKCsb/bYFYlrhOoRwAMhWs3LXdB5XREREhurpgW9/2zsKU269Fa67Ln9tGi41ffdzz/koIJkiEe+ovOEGqKrKT/sm69IlD9THj2evbigp8SECt24tnCnYszp3Ll03/aMf+Xz0o1m8GN7yFg/Tb32r11UvMIX+Y4WIiMi819XlwTqzDOENb/CJRuaCEHwc5+eeGxr+wYfU27bNh9WbDz26PT1e9nH06MjnkrJypfdQr18/vybombDeXq+XTvVOHzo0+m3NvCbp7ru93OPmmxdALczYFvazFxERmeO6u4v41rd8SGDwLPPGN3q4y7dEwid9ee65kfXHxcU+0+D110N5eX7aN1HxuJd7HD3q5SzZyj6qqtJlH/Ol533CQvAnn+qdbmjwbxmjSU39effd8OY3+ww4cpXCtYiIyBx15Qr8/Oe1XHONr0ci8KY35X+2wkTCSyWef97bmKm42Mtsr7tubpfXxuMepE+e9Mt4fORtiov9td6yxfNkQZ1719LiIfoHP/BAnTmm43ClpV4vneqd3rmzwF6M3FK4FhERmYMuXfLM09vr0ywWFXkp69q1+WtTIuEdnM8/n+5JTykp8cy1a5cPCjEXDQx4D/VYgRq8fnrrVg/WBVP20dYGjz0GjzziywsvZO+iT9m8Od07vW/fAhqce/oUrkVEROaQRAIOHPBSi9QgDMXFnnFWrcpPmwYHfei5Awegs3PovtJSD9Q7d0J05GSLeTcw4J2yp055sB4tUC9enO6lni+jmIypvR1+9rN0mH7++bFH9aishLvuSo/ske+fR+YxhWsREZE54soVz0GZw9eVlATe8Y78jEYRj8PLL3uo7u4eui8a9dKPHTs8YM8l/f3pko8zZ/zLQTZLlniG3LDBr89rXV1Dw/Szz47+xMF/Ctmzx+uM7r7bZyCaa2/kPKVwLSIiMgccOQJPPjm0Z7W+HlatusTy5dtmtS0DA/DSS3Dw4Mjz2mIxP0lx+/a5VTLR3+891KlAPVonbU2Nh+n16+d5oO7uhscf97rpRx7x8aZH65YHL9jfvRvuvNOXN7yhQLro5x6FaxERkTzq7oZHH/VAmBKJeKfi9dfDo4+O0fuYY21tHqpfeWXkdN7l5d6ea6+dOyOt9ffDq696oD57dvxAvWHD/JoNcojeXv/2leqZfuqp7HOup5j5G3bnnV4zvXfvPH7y88sc+XiIiIgsPCdP+nDCmROULFnieaiubnbakEh4QD1yBBobR+6vqPCJX7Zt80qCfOvrSwfqc+dGD9R1dd47vWEDLFo0q03Mjf5+D9CpnunHH88+k02mnTvTPdNvfKN/q5BZp3AtIiIyy/r7vTz2+PGh26+7zufgmI0Q29Hh9dQvv5x9SOPqau/43LIl/6G6tzcdqBsbxw7UqR7qeVfxMDAA+/ene6Z//vOxx5oG/xkhM0wvXTo7bZUxKVyLiIjMonPnvDOyqyu9rbLSf7lfuXJmjx2Cl58cOeIn/A1nBtdc4/XUq1bldyjjzEB97tzoo8YtXZquoZ5Xgbq72+ukH3/cf7746U+H/lFks3lzOkzv21fgc67PXwrXIiIisyAe9yw1fCbpLVvg9a+f2YEauru9jvqll0YOpQde+rFtmy/5HM64qys9ykdj4+iBetmydKCeN7MlNjZ6b/Tjj/vl88+PfQIi+JPMDNP5GotRJkXhWkREZIZduuS/9GdOER6L+aR369fP3HEbG72X+tVXs5dSrF7tvdRr1/pJlLMtkfBhB0+f9h715ubRb1tfnw7UlZWz18YpGRz0b1GZYfrVV8e/39q16TB95535nTFIpkzhWkREZIZkmxAGPDPt3esjcORaX5/PonjkyMipycFD/datXq6bjzKKnh4P0mfO+AgfY52jt3x5OlDP6QkCOzr85MOf/9yXJ58cOYVlNtu3+88Wt9/uNdPr1mla8QKgcC0iIjIDsk0IU1zsc3Vce23uj9fU5IH6xInsc4csX+5Zbv362T1BMQS4fNl7p0+f9l780UQisGKF131v2DAzXz5y4vTpdJD++c99KvGxZj8E/1Zzyy0epG+/3f8QNJpHQVK4FhERybHRJoS5887c9hYPDPiII0eOZC+pKCnxmu5rr53dHNfX573SZ854Du3tHf22FRXek79mjZcUz6WJaQB/Ew8eTAfpxx/3JzeeFSs8RKd6pm+4QTMgLhAK1yIiIjky3oQwufrFv6XFA/WxYx6wh6ur817qjRtnL6w2N6fD9MWLo5+MaOa96GvWeKiec523bW3+zSgVpp96auTc78OZwa5d6SB9++0q8VjAFK5FRERyYLQJYd70Jqitnf7jDw76MY4c8fA6XHGxh+nt22dnuOP+fj9hMlXuMVb+LCtLh+nVq+dQB248DocPwzPP+FAuTzzh66N9M0ipqIBbb033TN966zydqUZmgsK1iIjINPT1eQfnTE0I09lZxJNP+lB62U7+W7zYA/XmzRCNTu9Y42lrS4fpCxfGLjNetszD9Nq1/uUi7524iYS/Sc88k16ef378iVrAvxmkeqRf/3p/c+fKHPAy50z7L8PMioD9wLkQwrvMbD3wMFALPAu8P4TQP93jiIiIzDUzNSFMf7+P3HbsGDzyyDK2bh26PxLxExO3b/fS3pkSj6d7p8+cGXsAjGjUM2hqicVmrl3jCsHrojOD9P792YdPGa6oyGt4MsP0mjUz32YpGLn42vUR4CUgdYrGHwN/GkJ42Mz+DPgQ8PkcHEdERGROmIkJYfr7PcSeOOFBNluvcGWlB+qtW73UItdC8HruxkbPpo2N2UceSamrS5d7LFuWx97p5uZ0iH76ab/MVjuTzerVPorHzTf75S23zIOBtGUum1a4NrPVwDuBPwT+o5kZ8CbgfcmbfAn4NArXIiJSIHI5IUw8Dq+9lg7U2YJsakrya6/1IJvLABsCtLZ6iG5shPPnxx53uqTEs2hqdI+8DJXX0eEDh2f2Sp86NbH71tZ6iM5cNIW45Nh0e67/b+DjQGry0VqgLYSQGnzoLKC5OkVEZN7L1YQw8Xh6iu/Tp0efAbuuzk9QXLOmibvv3pr9RlMwPEyPNUwe+EmZqTC9fPksz+TY1+djSKd6o595xudwH++EQ/De55tuGhqkNYKHzIIph2szexfQFEJ41sz2TeH+9wD3ANTX19PQ0DDVpkxLZ2dn3o4tUoj0mZJC1NlZxIEDS2htTY9rV1QU2LmznVism6efHvv+g4Nw6VKMxsYYFy7EGBzMHvCqqwdYtaqXFSt6qKgYpLUVBgen95nq7Czm8uVSmpujXL5cSn//2Ok4Gk1QW9tHXV0/y5b1UVY2SE+Pz/p49OiUmzEuGxig4rXXqDx2jKqXX6bqlVeoPHGCyGjfPjIkSkro3LiRjm3baN+6lY5t2+hes2bo2aSvveaLyAyzMJFvf9nuaPYA8H4gDsTwmut/BO4GlocQ4mZ2G/DpEMLdYz3Wnj17wv79+6fUjulqaGhg3759eTm2SCHSZ0oKyeCgj8y2f//kJ4RJJLxu+cQJz3T9o5zav2SJ91Bv2OAjfww32c/UlStDe6bHG6I5FvOTL1NLtjbkXGurT8xy4EB6OXIk+6Ddw0UiXnie2SN93XVzaHw/mevM7NkQwp6Zevwp91yHEH4P+D2AZM/1x0IIv25mfw+8Fx8x5APAN3PQThERkVkzOOjVBwcODA2n400Ik0j4CCInT3oZ8GiBevFiD9MbN3q4no729nSYbmycWJhescKD9IoVMzyJSwg+7ElmiD5wwOthJmrjxqFBevdunXAoc9pMDNL4CeBhM/sM8DzwlzNwDBERkZyLx+Hll0eGahh9QphEwnuIT5zwQD3aCYHV1Z4TN26cXqDt6BgapjOHAcwmGk2H6ZUr/XnMSNlxX5/3PmeG6IMHJzb8Xcr69f7N5aabfNSOPXvm4BSOImPLSbgOITQADcnrJ4FbcvG4IiIis2GsUF1eDjfc4KN1pEp4Q/BAffKkL6OdFFhVle6hrqubWtt6eoo4ejQdpjs7x759aenQMF1TMwNhurl5ZFnHSy+NfnZmtkbu3OlB+oYbfLnuulmqSRGZWZpeSEREFqx43DPhwYOjh+pt23wyvhB8VsJUD/Vo5RcVFeka6mXLJteeRMJz68WL0NTkl/v3L+PMmdHvU1IyNEzndDbERCJ7WcdYDRqupiYdoFPLtm3ecJECpHAtIiILTipUHzgwcvbr8nIv6922zUNqU5OH6ZMnRy/BKC9P91BPZjKVrq50kG5q8jG0x5q0BTyTLl8+NEznZHi8y5f97M3U8uKL/q2jvX3ij7FxYzpAp3qlV6/W8HeyoChci4jIghGPe1nwwYMjQ3VFRToLNjbCj3/sl6OdlFhW5iXCGzd62B0vPw4Oen7N7JUer14afMi/1avTYbqubpphurV1aIg+fNinmmxqmvhjRKNe1pHZG33ddWMPnyKyQChci4hIwRsrVMdisGqV58XDh+HnPx/9cWKxdKBesWLsQN3RMTRINzdnn9J8uKoqH+pv2TK/fPHFC7zpTdsm9kQzXbkyMkQfPuzF4pNRW+td+Zk90lu3qqxDZBQK1yIiUrBSofrAgaEnHXZ3+/rixT7IxYkToz9GZaXPTrh+vfccZ+s1Hhjwko5UkG5qGhnisyku9hCdCtLLlnmPeKZxe6k7OvxJDg/RZ8+O34BMZWV+1uaOHb6keqZXrlRZh8gkKFyLiEjBicc9Xx486CE6Hoe2Nmhp8dC7dGl6Ku/hc6kVFXmeXL3aQ/XwASxC8E7hzCDd0jKxGbkXLx4apJcsmUSJR1eXF4oPL+eYzJjR4F30mSE6taxbN3RGQxGZEoVrEREpGAMD6Z7q5mYPva2tfk5eNOphedOmkYF28WLft2aNh+7ijP8d+/q8VzrzxMPRxrLOVFo6slc6Gp3Ak2huhldeubrsfOwxH6bk1VcnluAzG7B168gQvXGjQrTIDFK4FhGReW9gAJ59FlI5tLU1PZN2LOaBur4+HapLS73OOtU7nZrwr7PTqymam325fHn8caVTamrSQbq+HhYtGqOaoq8Pjh/3AH306JAwTUvLkJuOOzx2cTFs2eJlHJkhetOmod8SZO5J/cTS3u4ng+7YofesAOgdFBGReSk11fhjj8FTT43IpMRisHatB95IxEfZSPVOL13qpcqXL3tlRSpMjzYZzHCx2NAe6aVLPbAPEQKcPTcyQB896r3QEzm7MVNREWzePLInevPmLAeXOa2lBR56CD77Wf+iVVTkw8lEo3DffXDvvZqZch5TuBYRkXmjs9PnLzl1Cvbv98tUD3VKKlSvXQvXXOP102VlXrLc3AxPPunZZqKTCUYiI3ulh4w419EBLw7rfT561JeJjLU3XHm590Rv3QpbtnAkBLb/yq/4tgnVlcicduwY7N3rJwEM/zbX2QkPPACf/7x/a9y8OT9tlGlRuBYRkTnryhUv87h40S8vX/axp8+dGxmqy8pg1y7PoOXlHp6bm736YqKlyqWlPvJcba33dNfVeT12JBH33uZXXoF/GtYTPdmh7cDrRdatS4fo1LJli9erZBSFNzU0sH3XrskfQyZutsozWlrgjju8cH+0P8reXv+D37vX26Qe7HlH4VpEROaEzElWLlzwJdWx19fn2zNDdarnedEinx1x+XK/3UQHzygv9/CcCtK15T1UXz7p4/I9dcIvT55Mz3c+PM1PxJIlQ4Nz6vqmTd7FLvk12+UZDz3k3xjH+7YXgvdsf+5z8KlP5e74MisUrkVEJC9SgTkVpIdP/R2Pp8P2pUt++4EB71CMRLyDcdOmdE31WCN4LFoEtTWBuqJW6jpfpbb5KGWnjsGPkiH6xImp9UCDT6aycePQHuhUmK6r0xjRc9Vsl2fE4x7iJ1rY39sLDz4I99+v0V3mGYVrERGZFe3t6SB94YJnmkyJhGea8+e99OPSJc8jRUXeyRuL+age5eVeT7106cgh9SJhkJrEJWq7z1J75SR1l1+m9vwhSl495gG6o2PqT2DFiuy90OvWaYSH+SYf5RmHD09sDMdM/f1+xu3110/v2DKr9K+BiIjkXCLhvc6pIH3xYnrGwsFBv97d7UtXl9/28mXPM6WlHqRraoYOgrFkSXJkjqpeIi2XKTl0kdqu09RdOUFt00vUnTvIkjMvEBmcQvkGeIq/5hrvhd6wwS9Ty4YNPi+5FIZ8lGe0t0++BzoS8fvJvKJwLSIi05ZZ4pGabKW/30N0V1c6SHd3p0P2wID3TPf3e+aoqUlmj/iA3+lSJ5WDbSwbvMCG+FGWtx+jpull6lpeoY7LVNHBpAsuKiqGhubMIL12rZd4SGHLV3lGdfXQuqeJSCSGDU0j84HCtYiITFpmiUdjoy+pEN3V5QE6FaJTiouhLBYoi/TR19xFcUcXJd1dlHV1eT1IdxdlXc3U9jWyjZfYySHW8RqLaaOYSYSS5cuzh+eNG72WRDXQC1u+yjN27PATJSc6KxH47XfunPoxJS8UrkVEZEypyVYaG30kjrNnvWQ1FaaHdgCGZM1HF2X97VQPtFDe1US8rZP21kGudESweD8xeqmgm3K6KaeLxbSxk0Ns4xXquTh2j3RJidc5ZyvfWL/ee6dFRpOv8oziYh+B5IEHJtZrHov57WfqZEbNDjlj9CqKiAjgv0BfuACvveYTtZw96+sdHZ6X43EgJKCnGzq7PF0ne5wjnZ1U9FyiurOR6sEWquigiwq6qKSLcmL0soIuNtJNGT0UkaCIQa7hNTZzjDWcIUKy/rWkJD0DTOaybp0va9Zo9ASZunyWZ9x7r49AcvHi2PXeZj7A+oc/PP1jDqfZIWecwrWIyALT1+eTq5w+7SH63DkfoePShTj9V3qGFkgPX3p6gEAxcarpoJp2FtHGas5SSzMJIrSxmGZqWUIbtbSMOP6q6GU2r+ph3aZiStdvgnVvHhqiV6wYOQyISK7kszyjpsaH9httCEDwHuvFi/12uQ65mh1yVkw5XJvZGuDLQD0QgC+EEB40sxrgq8A64FXgV0IIrdNvqoiITFRfnw+G0N7uJxeeOdZD4yvtXDjZS+uFPrrb+kYG577sP1WXEKeCThZxhRVcYBVnWc1ZamhlCa0MEuEUGzjOJjqphLJySlLTHNbWQE0tdesq2XRjNZtur6d8Ta3qniV/8l2esXmzl2N87nN+omR/v3+ZTCTSvccf/nDug7Vmh5w10+m5jgMfDSE8Z2ZVwLNm9kPgg8CPQwh/ZGafBD4JfGL6TRURkUz9/XClLXDlXCftx5u4cqqFC8c7aXy1n+amOJ2tcbo6BunrikP/xE7gKqebZTSxgvOs4hxrOc0KLrCIK5STPEOxvp6u1Vs5XnMnB8p20Vy+ZkiQpqwM8JHrNm3yZcmSmXoVZpjqUgtTvsszamp8aL/77/cTJVN/Xzt3zlzJk2aHnDVT/hcihHAeOJ+83mFmLwGrgPcA+5I3+xLQgMK1iMik9Td30HniIh2nLtNxupXOs210nO+k42I37Zd6aW1O0Nkap3OghC4q6KSS+AT+WS8iQRk9LOVSMkSfZQ1nWMdpalaXw+rVyWUDrN4Lq1cTVq2mpWotjYnlvHquhMbG7I8djfp5hZs2+aAd85bqUgtbvsszUoqKZmeCGM0OOassjPcNZiIPYrYOeAzYCZwOISxObjegNbU+7D73APcA1NfX3/Twww9Pux1T0dnZSWVlZV6OLVKI9JkaRwiUtLdT2tyMXWxl4EIXA0099F/qoa+5n762QfquDNLbCQMDEeIU0UPZ1XE1eiinmzJ6KSMxxpgakWSALqOXCjpZHmliaXUny+u6qKuPU74yyuCyGvqWLr26DNTUEDL+I+3sLOby5VIuX47S3FxKf3/2OuhIJLB8eR+rVnWzbFnfvC+XLjt7lhs+8hGKOzsp6u8fsX+wtJR4ZSUHHnyQntWrZ7w9+kzNnOL2dlZ94xus/vrXsYEBQiSCJRIkSko498u/zLlf/EXiBTDOdMXx4+y+7z6Kh4+POYZ4WRnPf/azdG3aNIMty48777zz2RDCnpl6/GmHazOrBB4F/jCE8A9m1pYZps2sNYQw5g+Ce/bsCfv3759WO6aqoaGBffv25eXYIoVowX6m4nGvZTx//urSd6bJe5wb273HuamHzuY+OgbL6KCKftLTD/YSvRqiUwG6m/Iht0mJMMhqzlFOF4EiWiN1lFbHKFtcyqK6EtasgZUbyli6eTG1W+tYtG0FkWV1454k2N6eHrO6sdHLsDPZYJwljYcp7W2nP1ZN+c072LStmHXrhs6kOK+1tMD27WPXpYKXC9TXz0pd6oL9TM2mwcHZK8/Ih5/+FN79bi8LmahFi+Cf/snrtAuMmc1ouJ5W4ZiZlQBfB/42hPAPyc0XzWxFCOG8ma0AmqbbSBGRvOjv95CVWlJTDzY1ES5cpPdCG10XOui40EVHcz+doZwOquigik4qhwXjShJUZQ3QPZQzyMjgG4kYZeVFxCqLiVWVsqSij9e3f5fXnf1HisIAFBURIRCiUTo/eB8lv3svlWsnHvQ6O4eG6dEGT4h2tbDjkYfY+chnKR7sw4qLiIRBrBBLJFSXujDNVnlGvmh2yFk1ndFCDPhL4KUQwp9k7PoW8AHgj5KX35xWC0VEciUED07DgvLw64mLl+hu6qTrykByrObsS4LlwMjC4gGKswboXqIEDIpLoKyMSHmMWGUxi6pKiS2KEl1cRqymnGhdJbG6KkqqyygqNhYtgpVdx7j5Y3sp7mwj0t+bOpDr7WTJnz0AD489hFZ399AwPd6cGNEobEwc43X/dS/FHW3Y8NFECm3oLtWlSqHS7JCzajo917cD7wdeNLMDyW3346H678zsQ8BrwK9Mr4kiImPo7fVBm0fpYR5yvamJgf4EXcm5ATupvHqZDs3r6OHacQ87SIQ+ovRSlg7QsVq6y2oYiFVBWRlF5TFi1aVEF0VZsriMaG0FsdoKopWlxGJDSymiUR9RY/HioUtlJUTaWmD7HdAyuSG0enuHhum2trGfU0mJDzG9cqUvtdaC7bgDmhfI0F35mhZbZKble/jBBWY6o4X8DEY9m+auqT6uiCxgAwMelC9fnvjS1XX17r1ER+llrqeLDXRRkbWGOZsERj8xeivr6K2opbdsCb2xRfSWVtNbUkV/SQXFFVGiVVGiS7zHuaasiJUx/38pFvOwOlx1S1xvtQAAFwdJREFU9cgAvXix335UkyhVSLS2cfYTn+Ppuz9Fy8j5W4YoLvYRPVJhum54WfYfLLASiXxNiy0yG/I9/OACosE6RWRmDA76yWHDw/BY4TnLyTap0TKGjphRRg+r6WbL1fVuyhlkAsEoGoPqKkJlNQOVS+gpr6E3tpi+0mp6o9X0FlfSG6kgxMoorSolGotQWprsWU5eRqPe6zxaDisuzh6gq6unMETyJEsVIn29LPvKg7TedD9EhjYwEhkappctG+Mcx4VYIqG6VClkc2X4wQVA4VpExjY46P8Qt7Skl9bWoesZ2245c8Z7k1tbR+0dSWBZAvP6jJP90pfj9jQXFXntREWlh5zqKqiqhqoqBioW0RddRE+yt7mvqJy+hD9eIuFBNxWWy0phcTI0l5aOP/u2mR+2qmpkiK6oyOEEhFMoVYjE+1ly7hBt11zP0qWwalU6TE843C/EEgnVpUqhy9fskAuMwrXIQtHTM2YoHnV9gkM3JbCrNcgtlNNDzYignLrsZZQaCIt4Mq2sHHcJFZUMRCvpsxh9/UZfn38PSOX5wUH/PyMVnlO9ztlKNUZTXu7huarKc3vqelWVNzOSyJi9b7Aalk9/9r7UOZctLd7Jn3i0nd2JogkWs7hISYS9N7Sz5Bcm93yHWIglEqpLlYUgH7NDLjAK1yLzRSLh/wi2tWVfrlxJX08F5MygPNGf95MGKKaXGL3UXQ3Ew5fM7X1ERz5I+cSC8tWlrAwiEQYHvUOlr8+rE8B7gkPwJZGAxCAU93svc0WFnww42f8XotGhgTkzQFdWjpGTW1rgM9Ofva+/Px2im5vTb1XqOQPUXKlm9yRLFYotwbJN1TDVYA0Lt0RCdamyUBT68IN5pHAtMlvi8YmH42xLR8f4J5aNwnuVs4fi0cJyun7ZvAu3otzD8tXLjKW8/OrlqcuXWL9zl28blnZD8HMW+/s9m6SWq6G5y18mM+9xraiYRs8rHo4zw/PwAD2liU+OHRu9ZnGUoelC8Lc+FaBTYXoi1QetK3eQKI5C3yyXKizUEgnVpYrINClci4wnBK8hbm/Pvly5Mvr2zMA8mZAyWlOAAUroI5oswYhdvT58GdGrXFQ8MgxnBuWr2yqG7kv2Jo9ncNBDcX8iAVVVV0NzIuHL4KAvxcV+qFxMkR2NDm3y8NKNsrLpH2OIlhafrWys2ft6ewkXLxJ//V72f/EwFwdqRvRGj6e83DNbbS3U1haTuPc+woMPYLNZqrCQSyRUlyoi06BwLYVrcNADbUeHL6OF4LECcmpJJHLbtOQYyRNdeonRTyl9scWE8nIoK/cEVlbmS3k5lJcN3T48SJeWTuosOzPPV5FIuocZ0pkyMzQnEl5/XHvhMPGzh6gp76V15Q4oKiYS8ceYTClyJJIl6w/L/+Xl0y5vnrwJDolnIWBtbRT/P5+j6V2jD00XiXgHqIfodKAe8aXg4/fCl/JQqrCQSyRUlyoiU6RwLXNHf//QMJy6PtFtw/d3d89YUwPQT+mIpY/oiG3p3uUYfeVL6CtbRLysevRQfPWybGSInkZ3b+pkvlRYHvJ8kj3MIXgPayLh5RsDA9kfa8jjpqbG/slnicT7GMQoIjBYHOXwm+7j8J330leR7uFL1UiPFporKrwTNGejbUxRIpH+s+rogPaWOLv/x2cpnWDtenG8l10/eZAD77ifECmirGxogK6t9Tw6obc0X6UKKpFQXaqITJrCtUze4KCXSaSWzs6h69m2dXaml9GCcn//rDQ/VVqRLRyPGpZLq+gvW0R/rJr+WHVyBr4YxMo8XJSlLpPbyjK2pwJ0NJaTWoiSEg/Kmb3KqcsQyFqOEY/7eXeTHVltPNUXj/EL/30vpd1tFMeHBa++Tm78/gPc8MTnaf7Hx4jt2kxFxQz0NsczRuyorvZa4QkcJAT//jUkQLcP/fPMVHPmMLsn+QKWhH5+YcMhqu+4fvolKvkqVVCJhIjIpChcF6IQ/D/A7u6Ry/AQPFoQHmt9kqNO5EqcIgYouRqMs10fsS1W7aE4FYyjVfRHk8W4Ew3HsWS9cg5FIt6DG4tlL73IlOpRTgXlgQF/CyZTwztVqZcg8+VILRV9Lax66x1EOpqwUUoGIv29RJovUv8vk1NjF+cwgLW0eJnGGCN29FfWDAnMmdc7OiY3GEZpbzshMrlygKKSCPVl7ZCr2u98lSrMhRKJKX6JEhGZbfqXaTalfnPv6fElW/jN1ZLjGuGpiFNEnOKrgTdO8dihOBKjv2wRA7EqBso8CA9EfSzjgdIKQjQ5p3Q0mpxfOuq9waNdTmQmkGlITTaSmqkvswc5c0m9FcN7k/v6Jh/wpqu4ePSwnLmkpu8e8+X7g4egI09TYx87RrjDSxWsb+SIHfHPPED/n3yeb330MdrrN0/5MBUV6RFGllVVUxyZI0PT5atUIR/HncCXKPWai8hcYmGKQ3vl0p49e8L+/ftn/8DxOD/94Q+5Y8+edODNtvT2jr1/MrebA6E3UwK7GnwzL1PXs6+XMhCrIh6rZCBWSTxawUBp6rKCeEkZAyXl6dk7oqVQGs0IxRkBOXNbcfGsFdqWlIwfjDOlhorLDMipIeX6+/PztmZOkJJ6KTOvDw/LZWXTG9ZuiHjc59Fubp74fWpr/cS4MXo6E4n0x6i7O/tlvKmFu393O9GOJiJj/PuVMKO3qp6///ThITXfmVLjXGeOMpK6Xlk5rKlTec51dXDhgk6Am6qxhj2EofXem6f+JWqua2hoYN++ffluhkjBMLNnQwh7ZurxF27P9dq1cOYMd+S7HaMIcDXMDiZ7gCeyXL1tcRnxaEV6KS1nMFrGQHG5h+CScuKl5SRKo+ngOzwIj7ZeMrlRJ3KlqMjDYeqkvMwwnOphzbyembsyw3GqBjlVOZPv75epkJwZkEcLzJnXswblzJ/OI9WwaYZ+Op/C1Nihr5+mnxyiY/31V39gyfwRJ3V9PDf+80OUdF8ZM1gDREKgtKeNm578HGd/61PTH+d6IQ9Nlw8THPaQixc9gB8+rB5sEZkTFm64nkTgSAXdVHAdpGjq1yOlDJaWEY+WE49WEi8tz1jKiJeUM1gSY7AkNrTuYDJLScmc+A+9qOj/b+/uY+QqrzuOf8+87I5t/MIabAfsgEvtUtsIY1IClDWIiMQFGUxJW/uPpKhNUqm1/EcV1NaqI0Vpm1RN0uDaRbRR/wgRdaISUkJbUtTaxjQtocbYeBEQg0PtJd6AzfqVfZmZ0z/uzPru7OzuzO69c2d3fh9pNHPv3Jln7lqPfO7znOfc4FEZBFe+LisvxCurTKcI5xs3m3IOdeWjWnBcU5BcrwZNnZfT+QsnzpBLpakn6WagkOKFfz/DiUkMMFohz6r/3D5y8eQoMoN9rPqPh1n1+NZo+kQrl6ZrtBrLHsaSdiQiMgktG1w/77/KBVtOfyYHuVlBwJudQaFtBoVsLnidzVHItFPMtkNbFjLZi8FrWzYYwc1W7BvvmCYIemF46bVwxYnK1IjwiHD5OVyRovxd5efK1IlG5hNPVDo98euY8j/riGu18AjyjJgXX9V5x8Bi8WI6S7mCSPkR3jfa+wAdx+Zwb3+BegZ+rVhkIDd+/nEuV6WEd2l7ztEu2qmz5MnAQLAIL4pcYZWma4x8PrhQrPVKuq8vqGSyNaKLKBGRSWjZ4PrYF/6e7FuvceKNLhYtX8X7V6zEI64IUY0V8lz6ThdtfWcYyM3h5KKVFCwzFIyGH6MFvFD9deVNPipHgsnn6fhZF5kLQdunF6/EshfPuZwyEfc5R/G3zmQupohUPobeszxzjgXtpi6dAytW0jYzMyI4jnTNY4wjyOHa0+VH4d1TLLijk9TJ0St20NdHsaeH/l9Zy/e+1MX59skHfBO5Jbe3tdN+4yqumT28nHc4gB53EeWJM5CpM3hKpYKLnKioNF38JpB2FOlFlIjIJLRecF0Kfu7/q+2kBvvJu5ExJ59uZ/+tW9h/y2Yu5DqGjeyGA91yABse0Q2nNoTzfVOpi8HtjA9Osfr5HXzkR8GNNjyVxooFCpl2DnRu4ZW1mxmc3TGUTzyZwZfKVIvKm3yE2652k4+ohNtNFy62W8y28/a9W+i+fzM2v2P8ILnK/jFTvmsJcOfEEPiURpC9t3fkbarPnaP4F1+muOMRjn37OT5YvGxEoJzPjwyew/uqLZpc8/QOLu89PXpgXZJyJ3u+l+XP7uTAGHcMrEU2C+2XZHjzni1c++SXSQ/Wln/c/tAW1m+Y5KjinDn1XwHGUbGjGUrTTWdnztT/d4z6IkpEZIJaq1rIOCvPBzM5BmbM45+2PMfZRcuGRo0nO7I55o02gHwmx8DMeTz1+cmVDZtQ29kcgzPnsfuLz9G3ZBmZDNE8jv6E7MfWYqcbvMp/rAAX8PYcxbnzeO+J5+j/8LKhxY3h59Fej3VcqvcU921dQS6CCha1skKeTz20iNz52qtX9M2az2Nf7SGbSw/lfodzw6vtq3w91B9OnQrSXWrJP164MJoFZ6rY0RoOHoTbbht5J5+xzJ4N+/ZNy5FrVQsRiZaqhUSlhpXn2Xwf6XM9/MbfrI0k+IFg9Pber3aOGXRl8n2kzvZw39fW8sOvd1GY20E6HQSp5UWB4dfV9lV7nTlzistv78TGuMlHZrCPzJke7v7K8NX25dSUcP50+HXlvvJCt74+8JOnWHRnJ5wce5W/9/SQv2Utr+zqYuCSjmF52uHnavuqHZM+fYr7/3Tsv7X192Hv9jD7nrU8E9G/McCaH+6grdYKFhd6WbG7/hFks+Ej+POPd5Ep1jd13m4DfObmw9jqKZp/rIodrWHlyuBqrp7gur09mDkQEUlYbMG1ma0DHgbSwDfd/StxtVWTGleel8t33fCjnbz6wLahVJByLeTKPOjK15WVMK761g7a+2ovG7bmv3dy7MFtIxYHlvOhBweHLxwMv1+5b/muHczvPU2mhtX2+ZO9HPydnbx0z7ZJl6Zb8/QOFtaQqmDu2OleCtt3cmiSqQoAa56NP8Ctpu4KFvk+Vu99mMIfbSWbS4+bClPeNyI+3HcG2tJQQ/m6od+aTsHZePOPB4tFsuX6gnHkH6tix/SniygRmcJiSQsxszTwBnAXcBx4Edjk7q9WOz72tJAJTCWXp8/rvd1x2GSm7SfTbpJtT6V2B2bPZ893esi0p0eO+mfGmBGoeM6+epAZH78Na/QUdrNNnRcKcPgwB/bu5Ybbb483/1g3F5n+kkg7alJKCxGJ1lRNC7kJOOLubwGY2S7gPqBqcB27Caw8T+UHuLT7MKeWTDwIufSdLlL5xrcbRdujjdCPNmJffsz9af2pClkfYG3HYQZ++fqhutjlRZ2V2+Hn8OtMVxdtVl+7bQzw8SsiqC6QT2jxVbNNnZdujX36/ffjz3tVxY7pT2UPRWSKiiu4vhI4Fto+Dnw0fICZfQ74HMDChQvZs2dPTD8F5h46xHXudZ1s0Zz+dw/QO9sAr3Jb7Mp9PrKSyNmX8DpvZOgpZ2b+Bc5m+0N1psNtDW83vF0+LpWC+T9/GdL1zUpYpsjS+Xu5bPn79f3okDkcglR97RatyIUTezm9YOLtTuTfeLBY5PDevUEwOAmzjhzhhoGButrODwxw4PXXOT/J2odXrV/Phx9/nHS5APUYCm1t/N/69by9b9+k2hzPuXPnYu3Pw3R2wq23MuvoUTIXLpCfOZPzS5cGgf6hQ435DRKrzKOPcuX3v8/iJ57ABgfxVAorFilms3Q/8ADdGzaQ7+6G7u6kf2psGtqnRGTS4koL+SSwzt0/U9r+FPBRd99c7fjY00KSmj5Pctq+1c45yb91khUsmnDqXFPYEotS2lErlj1UnxKJVtxpIVHePiOsG1gS2l5c2peM8vR5PaKYPk+q3STbbrV24eLiq1yutuOjXHxVnjpfuHD09nO54H1NnctUVko7orMzeG6RwFpEpp64gusXgWVmttTM2oCNwFMxtTW+pIKfJIOuVjvnJP/WEFSwmDdvnLvbEE8Fi3L+8datMH9+MCI/d27wfNllwf6uLi3sExERaYDYbiJjZncD3yAoxfcP7v7nox3bkJvIJDV9nuS0faudc9IpEs1QwaIJps41hS0SLfUpkWhN1bQQ3P1f3X25u18zVmDdMElNnyc5bd9q55x0ikQzjCBr6lxERCRRsQXXTalK8DM4a1b8wU+SQVdSbbdau2UdHbBtWzB6vm8f/OAHwfOJE8F+5TyLiIhMa7GlhdSjIWkhlRp5w4sq7SYybZ9U263WbovTFLZItNSnRKI1VW8i0/waecOLKu0mIqm2W61dERERaVmtlRYiIiIiIhKjpkgLMbN3gbcTav4y4L2E2haZjtSnRKKlPiUSrV9y99lxfXlTpIW4++VJtW1m/xtn3o1Iq1GfEomW+pRItMws1oV+SgsREREREYmIgmsRERERkYgouIa/S/oHiEwz6lMi0VKfEolWrH2qKRY0ioiIiIhMBxq5FhERERGJiIJrEREREZGITPng2sw2mJmb2bURf++fmNkRM3vdzD5R2pczsx+b2UEz6zKzL0bZpkgziKNPmdl8M9ttZufMbEfFezea2Sul/rbdzCyqdkWaQSP7lJnNNrOXQ4/3zOwbUbUr0gxi6lN3mdn+0v9H+83szirHPGVmh8f7rikfXAObgOdLz5EwsxXARmAlsA74WzNLA/3Ane5+PbAaWGdmN0fVrkiTiLxPAX3ANuDzVd57BPgssKz0WBdhuyLNoGF9yt3Puvvq8oPgBm3fi7BdkWYQR596D1jv7tcBvw08Fn7TzH4dOFfLF03p4NrMLgFuA36XIBgu77/DzJ4Obe8wswdLr+82s9dKVyXbw8eF3Afscvd+dz8KHAFu8kD5D5stPbQiVKaNuPqUu5939+cJAoJwex8C5rj7/3iwuvpbwIY4zk0kCY3uUxVtLwcWAPsiOyGRhMXYpw64+zulzS5ghpm1h9r8Q+DPavmNUzq4JgiCn3H3N4CTZnbjWAebWQ54FPg1d78RGO3OkFcCx0Lbx0v7MLO0mb0M/Bx41t1fmOQ5iDSTuPrUaK4k6F9lQ31NZJpodJ8K2wh8x1UWTKaXRvSpB4CX3L2/tP0l4GvAhVp+4FQPrjcBu0qvdzH+9MC1wFul0WiAf6y3QXcvlKbaFgM3mdmqer9DpIk1vE+JTHNJ9qmNk/y8SDOKtU+Z2UrgL4HfK22vBq5x9ydr/YGZWg9sNmbWAdwJXGdmDqQBN7OHgDzDLxxydX59N7AktL24tG+Iu/ea2W6C/NBxk9tFml3MfWo03QT9q2xEXxOZqhLqU+W2rwcy7r4/yu8VSVLcfcrMFgNPAp929zdLu28BPmJmPyWImxeY2R53v2O075nKI9efBB5z96vc/Wp3XwIcBToJFnCsMLN2M5sHfKz0mdeBXzCzq0vbvzXKdz8FbCx9finBIqsfm9nlpe/DzGYAdwGvxXBuIkmIs09V5e4/A86Y2c2lKiGfBv558qci0hQa3qdCNqFRa5l+YutTpc/8C/DH7v5f5f3u/oi7X+HuVxPker8xVmANUzu43kRwdRH2BLDJ3Y8B3yUYUf4ucADA3T8Afh94xsz2A2eB05Vf7O5dpc+9CjwD/IG7F4APAbvN7BDwIkHOdbUFkSJTUWx9CqB01f914EEzO16qykPp898kWDj8JvBvEZ6TSJKS6lMAv4mCa5l+4uxTm4FfBL4QKmW5YCI/suVuf25ml7j7udIo2U7gJ+7+10n/LpGpSn1KJFrqUyLRanSfmsoj1xP12VK1jy5gLsEKUhGZOPUpkWipT4lEq6F9quVGrkVERERE4tKKI9ciIiIiIrFQcC0iIiIiEhEF1yIiIiIiEVFwLSIiIiISEQXXIiIiIiIR+X9WbLcUU0bnyAAAAABJRU5ErkJggg==\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 = 3\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.4\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, 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": "markdown",
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.1.4.1 Simple logarithmic extrapolation](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4.1-Simple-logarithmic-extrapolation)",
     "section": "2.1.4.1 Simple logarithmic extrapolation"
    }
   },
   "source": [
    "### 2.1.4.1 Simple logarithmic extrapolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.1.4.1 Simple logarithmic extrapolation](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4.1-Simple-logarithmic-extrapolation)",
     "section": "2.1.4.1 Simple logarithmic extrapolation"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.2009787  0.24331935]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'New cases')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAE/CAYAAACuKr76AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU1f3H8fdJMglZSFgDgpElgIAsQYNUFA2ISFQUMCIiW7RQpVqLuFb7q621tSqiFq0iOgOJghhXqlBRikvrBpKwbwFkC6CsCeskOb8/ZqCRzUCWm5n5vJ4nDzPn3pn75Z4EP56cOcdYaxERERERCWVhThcgIiIiIuI0hWIRERERCXkKxSIiIiIS8hSKRURERCTkKRSLiIiISMhTKBYRERGRkKdQLCIiIiIhT6FYRKQCjDHrjTHbjTGxZdp+aYyZ52BZIiJymhSKRUQqLhy4y+kiRETkzCkUi4hU3JPAPcaYOic6aIxpa4yZY4zZaYxZaYwZ5G9vYYzZbYwJ8z9/2Rizvczrsowxvz3JeyYZY942xvxgjNlhjJnob082xsz1t/1ojHmtbF3GmPuNMZuNMYX+Wi73t4cZYx4wxuT7XzvDGFPPf6yWMSbb377bGPOtMaZRJd07EZEaQaFYRKTi5gPzgHuOPeCfVjEHeB1IBAYDLxhj2ltr1wF7gS7+0y8Fiowx7fzPLwM+PcF7hgP/BL4HmgNNgelHDgN/BZoA7YAk4BH/684F7gC6WmtrA1cC6/2vuxPo779mE2AX8Lz/2Aggwf9e9YHbgAPluTEiIoFCoVhEpHL8H3CnMabhMe3XAOuttW5rbbG1diHwFnCD//inwGXGmMb+5zn+5y2AeCDvBNe6EF9wvddau89ae9Ba+wWAtXaNtXaOtfaQtfYH4Gl8QRegBIgC2htjXNba9dbafP+x24CHrLWbrLWH8AXpDGNMBODFF4ZbWWtLrLULrLV7z/hOiYjUQBFOFyAiEgystUuMMf8EHgCWlznUDOhmjNldpi0CyPI//hS4FtgEfIZvxHkYcBD43FpbeoLLJQHfW2uLjz3gn9bwLNADqI1v8GOXv8Y1/ukYjwDnGWP+Bdxtrd3ir/MdY0zZ65UAjfy1JgHT/VMxsvEFaG85bo2ISEDQSLGISOX5AzAK33SGIzYCn1pr65T5irPW3u4//im+AJvmf/wFcDEnmTpR5j3P8Y/iHusvgAU6WmvjgaH4plQAYK193Vp7Cb4QbIG/lXnP9GPqrGWt3Wyt9Vpr/2itbQ90xzf6Pfx0boyISE2nUCwiUkmstWuAN4DflGn+J9DGGDPMGOPyf3U9Mm/YWrsa3/zcofjC815gG3A9Jw/F3wAFwOPGmFj/B+Eu9h+rDRQBe4wxTYF7j7zIGHOuMaaXMSYK30j0AeDIyPCLwGPGmGb+cxsaY67zP+5pjOnon8u8F990ihONYIuIBCyFYhGRyvUn4OiaxdbaQqAPvg/YbQG24hudjSrzmk+BHdbajWWeG+C7E13AWlsC9ANaARvwTb240X/4j8D5wB7gA+DtMi+NAh4HfvTXkQg86D/2LPA+8JExphD4CujmP9YY31znvfimhnzK/6Z/iIgEBWOtdboGERERERFHaaRYREREREKeQrGIiIiIhDyFYhEREREJeQrFIiIiIhLyFIpFREREJOTViB3tGjRoYJs3b+50GUFn3759xMbG/vyJUqOo3wKX+i5wqe8Cl/ouMDnZbwsWLPjRWtvw2PYaEYqbN2/O/PnznS4j6MybN4+0tDSny5DTpH4LXOq7wKW+C1zqu8DkZL8ZY74/UbumT4iIiIhIyFMoFhEREZGQp1AsIiIiIiFPoVhEREREQp5CsYiIiIiEPIViEREREQl5CsUiIiIiEvIUikVEREQk5CkUi4iIiEjIUygWERERkZCnUCwiIiIiIU+hWERERERCnkKxiIiIiIQ8hWIRERERCXkKxSIiIiIS8qokFBtj+htjXjbGvGGM6VMV1xAREamQ/HwYMwbi47msVy+Ij/c9z893ujKRoGSt5buC77DWOl3KCZU7FBtjXjXGbDfGLDmmva8xZqUxZo0x5gEAa+271tpRwG3AjZVbsoiISAXNmgWdOsHkyVBYiLEWCgt9zzt18h0XkUqxtWgrT/7nSTr8owMXTLqALzd96XRJJ3Q6I8UeoG/ZBmNMOPA8kA60B24yxrQvc8rD/uMiIiI1Q34+ZGTA/v3g9f70mNfra8/I0IixSAVtKdxCv2n9OPvps7nv4/tIiErgpWteokNiB6dLO6GI8p5orf3MGNP8mOYLgTXW2rUAxpjpwHXGmOXA48Asa+13lVSriIhIxY0ff3wYPpbXCxMmwMSJ1VOTSJBYWLCQH/b/QJ/kPtSPrs/3u7/n3u73MiJlBG0btHW6vFMypzOvwx+K/2mt7eB/ngH0tdb+0v98GNANWAWMAL4Fcq21L57gvUYDowEaNWp0wfTp0yv0F5HjFRUVERcX53QZcprUb4FLfRcYLrn6aiL27//Z84pjYvjigw+qoSKpCP3cOW/34d18vP1jZm+dTf6+fJrFNMPT1QP45hEbY457jZP91rNnzwXW2tRj28s9Unw6rLXPAc/9zDmTgEkAqampNi0trSpKCWnz5s1D9zXwqN8Cl/ouQBw4UK7TIg4eVH8GAP3cOevxLx7n95//nuLSYro26crdl93N4A6DqRdd75Svq4n9VtFQvBlIKvP8bH+biIhIzRQX5/tQXXnOE5GfWLJ9Ce6Fbu648A5a1G1Bl8Zd+G233zIiZUSNnStcXhVdku1boLUxpoUxJhIYDLxf8bJERESqyNCh4HKd+hyXC4YNq556qkKZ5eYIC9Nyc1IhOw/s5Plvnid1Uiod/9GRv3/zd77e/DUAV7a6kif7PBnwgRhOb0m2acCXwLnGmE3GmFuttcXAHcC/gOXADGvt0qopVUREpBKMG1e+UDx2bOVet7qC6jHLzaHl5qQC9h3eR7NnmnHHrDsosSU82/dZtozbwuAOg50urdKdzuoTN52k/UPgw0qrSEREpColJ0NOjm/ZNa/3pytRuFy+r5wc33mVZdas4693JKhOmeK7Xnp6xa9Tdrm5Yx25dkYGLFpUuX8/CRorflyBJ9fD+t3rmZ4xndjIWCZcOYHUJqmkNE5xurwq5eg2z8aYfsaYSXv27HGyDBERCTXp6b5gOHo0xMdjjfGN3I4e7WuvjIB6RHWui3w6y82J+O05uIeX5r/ERa9cRLvn2/HUf5/iQPEBDpccBuCX5/8y6AMxOByKrbUzrbWjExISnCxDRERCUXKybx3iPXv4dO5c2LPH97yyR1CrM6hmZ5fvWllZFb+WBLSS0hK8Jb7vlal5U7ntg9soPFTIU1c8xaa7N/He4PeIDI90uMrq5WgoFhERCXrVGVSLiir3PAk6a3au4eG5D9Pi2RZMyZsCwNBOQ/l21Lcsvn0x47qPo3FcY4erdEaVrFMsIiIiftUZVLXcnJyAtRZ3rht3rpsvNnxBmAnjyuQraVm3JQB1o+uSGn3cXhYhRyPFIiIiVam8AbQygmooLDcn5VJqS1n+w3IAjDG8OP9Fftj3A49f/jgbfruBD2/+kF4tejlcZc2ikWIREZGqNHSob5WJU02hqKygOm6cbzWLn7tWZS83JzXGul3rmJo3lSl5UygoKqBgXAF1atXhw5s/pH50/RNuuSw+Wn1CRESkKlXnushHlpuLiTn+mi6Xr72yl5uTGmHBlgX0mtKLls+15I+f/pHkeslM7jeZWhG1AGgQ00CB+Gdo9QkREZGqVN1B9Zjl5o5uFFIVy82JY6y1fLHhC/K25gEQ44phw54NPNrzUdb/dj1zhs3h5k43Hw3F8vM0fUJERKSqHQmqEyb4VpkoKvLNIR42zDdCXNkjt0eWm5s4sXLfVxy3cc9GpuZNxZPnYc3ONdzc8WayB2bTrmE7Vt+5WqPBFaBQLCIiUh0UVKWCbnnvFjy5HiyWtOZpPNzjYa5vf/3R4wrEFaNQLCIiIlLDWGv5evPXzFg6gyeueIKIsAg6N+rM7y/9PSNSRhxdTk0qj0KxiIiISA2xpXALWXlZePI8rPhxBdER0QzvPJyUxinc9Yu7nC4vqCkUi4iIiNQAuVtzuWDSBZTaUi5OupjJ/SZzw3k3EB8V73RpIcHRUGyM6Qf0a9WqlZNliIiIiFQray3fFXyHJ9dDYmwiv7/s93Rq1Ik/9/wzGe0zaF2/tdMlhhxHQ7G1diYwMzU1dZSTdYiIiIhUh21F23ht8Wt4cj0s3r6YqPAoRl8wGoAwE8aDPR50uMLQpekTIiIiIlXIW+LFFe5bo/r+j+9nSt4UujXtxotXv8iNHW6kTq06DlcooFAsIiIiUiXytubhyfWQvTibOcPmkNI4hYd6PMR9F99H+4btnS5PjqFQLCIiIlJJ9h3exysLX8GT62Hh1oW4wlxc1/Y6wk04gOYK12AKxSIiIiIVUFxazKa9m2hepzkAD819iDb12/Bc3+cY0nEI9WPqO1uglItCsYiIiMgZWPbDMtwL3WQtyqJhbEMW3baI2MhYVt6xkia1mzhdnpwmhWIRERGR0/DPVf/k0c8e5ZvN3xARFsHVra9mZMrIo8cViAOT1ikWEREROYWS0hI+XvsxXc7qQmJsIrsP7uaA9wBP93mamzvdTGJsotMlSiUIc/Li1tqZ1trRCQkJTpYhIiIicpxVO1bxu09+R7NnmtH3tb5kL8oGYEjHIeTdlsfYi8YqEAcRTZ8QERERKcNb4uXyqZfz+YbPCTNhpLdK59m+z3JNm2sA3yYbEnwUikVERCSkldpS5q2fx4ItC7j34ntxhbvokNiBa9pcw7BOwzir9llOlyjVQKFYREREQtLaXWvx5HqYkjeFDXs2UC+6Hrd3vZ24yDheuPoFp8uTaqZQLCIiIiEnKy+L4e8Ox2C4IvkK/tb7b1x37nVEu6KdLk0colAsIiIiQc1ay+cbPseT6+GaNtcwsN1ALm95OX/u+WeGdx5OUkKS0yVKDaBQLCIiIkFpw54NTMmdgifPw9pda4mLjKNzo86Aby3hhy59yOEKpSZRKBYREZGgUVJaQnhYOADXvH4Ni7cvpmfznjxy2SMMbDeQ2MhYhyuUmkqbd4iIiEhAs9by5aYvcS90Mzt/Nit+vYLYyFhevOZFmtRuQvM6zZ0uUQKAo6HYWjsTmJmamjrKyTpEREQk8Gzft53XN7zObc/fxsodK4lxxXBD+xsoPFxIbGQs3ZO6O12iBBBNnxAREZGAcbD4IHsP7SUxNpFNezfx8rqX6XFOD+6/+H4y2mdQO6q20yVKgFIoFhERkRrNWsv8LfNx57qZtmQa/dv2x32dmy6NuzCt2zQG9x3sdIkSBBSKRUREpMZ6ecHLPPv1syz9YSm1ImoxsN1ARnQeAYAxhsa1GjtcoQQLhWIRERGpMQ6XHGb2mtlc0+YawkwYy35YRnxUPC9d8xI3nncjCbUSnC5RgpRCsYiIiDgud2su7oVuXlv8GjsO7GDeiHlc1vwynuzzJBFhiitS9fRdJiIiIo5Zv3s9/af3J29bHpHhkfRv25/MlEwuOecSAAViqTb6ThMREZFq4y3xMmvNLPZ79zO4w2Ca1m5K47jGjDp/FDd1vIl60fWcLlFClEKxiIiIVLkl25fgXugme3E22/dtJ7VJKoM7DMYV7mL20NlOlyeiUCwiIiJV6/459/PEf58gIiyCfm36kZmSSd9WfZ0uS+QnFIpFRESk0hSXFjMnfw6ePA+P9nyUNvXbcE2ba2hSuwlDOg6hYWxDp0sUOSFHQ7Exph/Qr1WrVk6WISIiIhW04scVeHI9TM2bSkFRAfWj6zOs0zDa1G9Dj2Y96NGsh9MlipySo6HYWjsTmJmamjrKyTpERETk9FlrMcZQeKiQlBdTKC4tJr11OiM7j6Tfuf2IDI90ukSRctP0CRERESm3ktIS5q6biyfPQ0FhAXNHzKV2VG3evOFNujbtSuM47TAngUmhWERERH7Wul3reGXhK0zNm8rGvRupU6sOQzoMwVvixRXuot+5/ZwuUaRCFIpFRETkhAoPFRJmwoiNjGXO2jn89Yu/0ie5D0/1eYprz72WWhG1nC5RpNIoFIuIiMhRpbaUz77/DHeum5xlOTx5xZOM6TqGIR2HcHXrq2ka39TpEkWqhEKxiIiIUGpL+fNnf8aT62Hd7nXUjqzNzR1vpntSdwDiIuOIi4xzuEqRqqNQLCIiEqL2e/fzzeZvSGueRpgJ45N1n9Cybkse7fkoA9oNIMYV43SJItVGoVhERCSEWGv578b/4s51M2PpDA4WH6RgXAH1Y+ozZ9gcLaMmIUuhWEREJER89v1n/PL9X7J652piXbEMOm8QI1NGUi+6HoACsYQ0hWIREZEgdcB7gHdXvEuzOs3ontSds+PPpkntJjzU4yGub3+95giLlKFQLCIiEkSstXyz+RvcuW6mL5nOnkN7uLXLrXRP6k7Lui2ZN3Ke0yWK1EgKxSIiIkHk6tevZtaaWURHRHN9++vJTMkkrXma02WJ1HgKxSIiIgHqUPEhZq6aSc6yHKYOmEpkeCQ3nncjA9oOYNB5g0ioleB0iSIBQ6FYREQkgFhrWbh1Ie6Fbl5f8jo7D+ykae2mrNm5hvYN2zMiZYTTJYoEJEdDsTGmH9CvVatWTpYhIiJS41lrMcYwf8t8Lpx8IVHhUfRv25/MlEx6t+xNeFi40yWKBDRHQ7G1diYwMzU1dZSTdYiIiNRE3hIvH6z+AHeum5Z1WjKh7wRSm6Tivs7NdedeR93ouk6XKBI0wpwuQERERH5qyfYl3P2vu2n6dFMGvDGArzd9Tf2Y+gAYYxiZMjK0A3F+PowZA/HxXNarF8TH+57n5ztdmQQwzSkWERGpAXYe2EndWnUxxvDc18/hyfVw7bnXkpmSyZWtriQiTP/JBmDWLMjIAK8XvF4MQGEhTJ4MU6ZATg6kpztdpQQgjRSLiIg4pLi0mA9WfUDGjAwaP9WYb7d8C8AfLvsDBeMKyBmUw9VtrlYgPiI/3xeI9+/3heKyvF5fe0aGRozljOinTEREpJrtOrCLv37xV7IWZbG1aCsNYxpyx4V3kBibCEDT+KYOV1hDjR9/fBg+ltcLEybAxInVU5MEDY0Ui4iIVIPdB3ezsGAhALUiauHJ9dCtaTfevfFdNt29iaevfJrmdZo7W+SZKjPHl7Cwqpvjm51dvlCclVW515WQoJFiERGRKlJSWsIn6z7BnevmneXv0KxOM1b8egXRrmi+/+33RLuinS6x4o6Z4wtU3RzfoqLKPU+kDIViERGRKpCVl8Xv5v6OTXs3US+6HqPOH0Vml8yjx4MiEJed43usIyE5IwMWLYLk5IpfLy7OF7jLc57IadL0CRERkUqw99BeJn83mU17NwEQGxlLp0admJExgy13b+HvV/2d8886H2OMw5VWotOZ41sZhg4Fl+vU57hcMGxY5VxPQopCsYiIyBkqtaXMXTeXYe8Mo/FTjRk1cxTvLH8HgAFtB/DBkA+44bwbiIqIcrjSKlLdc3zHjStfKB47tnKuJyFF0ydERETOwKHiQ5z3wnnk78onISqB4Z2Hk5mSyYVNLwQIrhHhk6nuOb7Jyb45ysfOYQZfGHa5fMcrY6qGhByFYhERkXIoOlzEW8veYsn2JTzZ50miIqIY2mko59Y/l/5t+wfHHOHT5cQc3/R03xzlCRMgKwtbWIipXds3ZWLsWAViOWOaPiEiInIS1lo++/4zbnnvFho/1ZiR743k/VXvs9/r+2DZI2mPcFPHm0IzEINzc3yTk33rEO/Zw6dz58KePb7nCsRSAQrFIiIiJ/Hydy9zmecy3lz2JoM7DObzzM9Z8esVxLhinC6tZtAcXwkimj4hIiIC7Pfu553l7+DJ8zCs0zCGdx7OgLYDiI6IZmC7gcRGxjpdYs2jOb4SRBSKRUQkZFlr+WrTVzy16ik+/+pz9h7aS/M6zTH4PiTXMLYhwzprea9TOmaOL0VFvjnEmuMrAUahWEREQs7eQ3uJj4rHGMOvP/w1y7cvZ1DHQWSmZHJps0sJM5pdeFqOzPGdONHpSkTOmH7qRUQkJBwsPsiMpTNIfy2ds58+m72H9gKQPTCbty56iyn9p5DWPE2BWH4qPx/GjIH4eAgL8/05ZoyvXYKKRopFRCSordu1jqf++xTTlkxj18FdJMUncVe3uyguLQagfcP2bI/Y7nCVUiPNmnX8fOnCQpg8GaZM8c2XTk93tkapNArFIiISdLYWbeWA9wAt6rZgn3cfr+a+ysB2AxnZeSS9WvQiPCzc6RKlpsvP9wXi/fuPP3YkJGdk+OZTa950UHD0d0TGmH7GmEl79uxxsgwREQkCh0sO8/byt+k3rR9nP302D//7YQA6JHZg+z3beW3ga1yRfIUCsZTP+PHl28J6woTqqUeqnKOh2Fo701o7OiEhwckyREQkwD322WM0fbop18+4ngVbFnBP93t4uMfDR4/XjqrtYHUSkLKzyxeKs7Kqpx6pcpo+ISIiAefH/T+SsyyHUeePIjwsnMMlh0lrnkZmSiZ9kvsQEab/vEkFFRVV7nlS4+lfDRERCQjeEi+z18zGnevmn6v+ibfUy3kNz6NHsx78secfnS5Pgk1cnO9DdeU5T4KCQrGIiNR4q3espoe7B9v2bSMxNpE7L7yTkSkj6dioo9OlSbAaOtS3ysSpplC4XL5NSiQoKBSLiEiNs/PATqYvmY7BcHvX22lZtyVXtb6K/m37k94qHVe4y+kSJdiNG+dbdu3nQvHYsdVXk1QphWIREakRSlavYs7zd+Pe/i/eTS7mcAT0PZTE7fX6EJ6czKvXvep0iRJKkpN96xAfu04x+MKwy+U7ruXYgoa27REREefNmsWd955Het0P+CSpmNvmw8IXYdZTW6FTJ98mCiLVLT3dtw7x6NE/3dFu9GhfuzbuCCoaKRYRkWq35+Ae3lj6Bp5cDy91fpiOGTfwy4RiLl8F/VZBZMmRM7VJgjgsORkmTvR9SVBTKBYRkWpRakuZu24u7lw3by9/m4PFB2nfsD07sl4Cr5fzC+D8gpO8+MgmCQomIlJFFIpFRKRKHfAeINoVzb7D+7hu+nVEhkeSmZJJZkomqU1SMQkJ5d8kQaFYRKqIQrGIiFS6wkOF5CzLwZ3rZp93HwtGL6B2VG3mDp9L58adqRVR638na5MEEakBFIpFRKTSfFfwHc99/Rw5y3LY591Hm/ptGNl5JMWlxUSERdDt7G7Hv0ibJIhIDaDVJ0REKiI/H8aMgfh4LuvVy/fJ9DFjfO1VfL2jn4SvyuuVw/rd69l9cDcAi7Yt4u3lbzOk4xD+e8t/WfHrFTzY48FTb7s8dKhveatT0SYJIlLFFIpFRM7UrFm+5cImT4bCQoy1vhHPyZOrZhmxY65HVV/vFPZ795OVl8XlUy+nxbMtcC90AzC4w2C23rOVSf0mcVHSRRhjfv7Nxo0rXyjWJgkiUoUUikVEzkR+vm+ZsP37j/+QmNfra8/IqLwR3Oq+3kmU2lJGzxxN46caM/zd4azfvZ4/pf2J69tfD0CtiFrEuGJO702PbJIQE3N8OHa5fO3aJEFEqphCsYjImRg/vnwrJkyYEJjXK2PT3k28seQNAMJMGNv2bWNgu4HMGzGP1Xeu5veX/Z5zEs6p2EW0SYKIOEwftBMRORPZ2dW7jFg1X+9g8UHeXfEu7lw3c/LnEB4WzhXJV1Avuh7v3vhu+aZFnC5tkiAiDlIoFhE5E9W9jFg1Xm/W6lkMeXsIuw/u5pyEc3j40ocZ0XkE9aLrAVRNIBYRcZhCsYjImajuZcSq8HoFhQVkL8qmc+PO9EnuQ4fEDlzd+moyUzLp2aInYUYz7UQk+OlfOhGRM1Hdy4hV8vUOFR8iZ1kO17x+DUkTkrjv4/v4KP8jAJISksgemM3lLS9XIBaRkKF/7UREzkR1LyNWydfrOaUnN7x5A7lbc7nv4vtYecdKnurzVCUUKiISmDR9QkTkTBxZRiwjw/cBt7IfgnO5fF+VuYxYBa63fd92Xlv0Gm+veJuPhn5EtCuaBy55gKjwKHq37E14WHjl1CgiEsA0UiwicqaOWUbMGlO1y4idxrJl3hIv7614j/7T+9P06abc/dHdHC45zObCzQBce+61XNnqSgViERE/jRSLiFREmWXEPp03j7S0tGq73okcKj5EVEQUC7cupP8b/WkU24jfdvstI1NGcl7ieVVbm4hIAFMoFhEJcDv272Dakmm4c92c3/h8Xr72Zbo26cpHQz8irXkarvCfmYssIiIKxSIigeqTtZ/w4oIXeX/l+xwuOUyXxl34xdm/AHxrCV+RfIXDFYqIBA6FYhGRALLyx5W0qd8GYwzvrHiHeevnMSZ1DCNTRtK5cWenyxMRCVgKxSIiNdyuA7t4Y+kbuHPdfLP5Gz4b+Rk9mvXg0Z6P8vSVTxMZHul0iSIiAU+hWESkhtq+bzt3zb6Ld5a/w6GSQ3RI7MD4PuNp17AdAHWj6zpcoYhI8FAoFhGpQVbvWM3GvRvp1aIXdWrV4buC7/jl+b9kZMpILjjrAowxTpcoIhKUFIpFRBy299Be3lz6Ju5cN//Z+B9a1WvFqjtWERkeyYpfr1AQFhGpBpUeio0xLYGHgARrbUZlv7+ISDB57uvneODjBzhQfIC2Ddryt95/Y2inoUeDsAKxiEj1KNeOdsaYV40x240xS45p72uMWWmMWWOMeQDAWrvWWntrVRQrIhLo1u1axyPzHmHtrrUAtKnfhuGdh/PVrV+xbMwy7rv4PprUbuJwlSIioae8I8UeYCIw9UiDMSYceB64AtgEfGuMed9au6yyixQRCWRFh4t4a9lbePI8zFs/D4OhWUIzWtZtSd9Wfenbqq/TJYqIhLxyhWJr7WfGmObHNF8IrLHWrgUwxkwHrgMUikVE/A4WH6TZM83YeWAnreq14s89/8zwzsNJSkhyujQRESnDWGvLd6IvFP/TWtvB/zwD6Gut/aX/+TCgG/AH4DF8I8iTrbV/Pcn7jQZGAzRq1OiC6dOnV+gvIscrKioiLi7O6TLkNKnfKqbW5s0kzZhBo48/JvzAAUqio0EkyqkAAB2pSURBVNnWuzcbBw3iYNOmVXrtoqIi9kXs46NtH7Fx/0Z+1+53ALy7+V2S45LpEN9Bc4RrKP3cBS71XWByst969uy5wFqbemx7pX/Qzlq7A7itHOdNAiYBpKam2rS0tMouJeTNmzcP3dfAo36rgFmzYPRo8Hp9X0DE/v00nTWLph9/DDk5kJ5e6Zc94D3AOyve4em8p/lu93dYLD2b96Tbxd2IdkWTRlqlX1Mql37uApf6LjDVxH6rSCjeDJT9/d/Z/jYRkeqXnw8ZGbB///HHjoTkjAxYtAiSkyt8OWstJbaEiLAIPLkexnw4hkZRjfi/y/6PEZ1H0KJuiwpfQ0REqk9FQvG3QGtjTAt8YXgwMKRSqhIROV3jxx8dHT4prxcmTICJE8/4MlsKt5CVl4Unz8Nvu/2WX6X+ips63kTbBm2x6y290nqd8XuLiIhzyrsk2zTgS+BcY8wmY8yt1tpi4A7gX8ByYIa1dmnVlSoicgrZ2eULxVlZp/3W1lpmLJ3BVa9dRdKEJB745AEaxDSgabxvjnKdWnXo2aInYaZc/6SKiEgNVN7VJ246SfuHwIdnenFjTD+gX6tWrc70LUREfIqKKvU8ay3rd6+nRd0WGGOY8NUENu3dxIOXPMiIziNoXb91BYoVEZGaxtFtnq21M4GZqampo5ysQ0SCQFwcFBaW77xT2Fa0jexF2XjyPKzesZqCcQXUja7L24PeJjE2kfCw8EoqWEREahL9rk9EgsPQoeBynfoclwuGDTvhoSXbl3DttGtp+nRT7plzD7GuWJ7t+yyR4ZEAnFX7LAViEZEg5uhIsYhIpRk3DqZMOfW8YpcLxo49+jRvax5hJoyOjToSGR7J/C3zGXfROEamjKRdw3bVULSIiNQUCsUiEhySk33rEGdk/GSdYsAXhl0uyMnhx7MSeP3r53Dnusndmsug8wbxRsYbtKnfho1jN2o0WEQkRGn6hIgEj/R03zrEo0dDfDyEhfn+HD0aFi3iN8yiyfgm3DX7LsJNOH9P/zsvXPXC0ZcrEIuIhC5HR4q1+oSIVLrkZN86xBMnsnT7UqYtmcYfLvsDrnAXLba34M4L72Rkykg6NurodKUiIlKDaPUJEQkquw7sYtqSaXhyPXy75VsiwiLo16Yf3c7uxtiLxv78G4iISEjSnGIRCRrLflhGl5e6cLjkMJ0adeLpPk9zc6ebSYxNdLo0ERGp4RSKRSRgrfxxJZ5cD3GRcTx06UO0bdCW+7rfx8B2A+lyVhenyxMRkQCiUCwiAWXPwT3MWDoDd66bLzd9SbgJZ2inoQCEmTAe7fWowxWKiEggUigWkRqv1JZiMBhjuOeje5i8cDLtGrTjid5PMLTTUM6qfZbTJYqISIBTKBaRGit/Zz6eXA9T8qbwzo3vcEGTCxjXfRyjLhhF1yZdMcY4XaKIiAQJLckmIjXKweKDTF8yHXeum8++/wyDoU9yHywWgLYN2jpcoYiIBCNHN++w1s601o5OSEhwsgwRcVipLWVL4Zajj38z6zcUFBbwl15/YcPYDcweOpvUJqkOVykiIsFM0ydExDHrd69nat5UPLkeol3RLLl9CTGuGHJvy6VFnRaaHiEiItVG2zyLhKL8fBgz5qdbIY8Z42uvBnPXzeXyqZfT4tkW/GHeH2hRtwUPXvLg0SkSLeu2VCAWEZFqpZFikVAzaxZkZIDX6/sCKCyEyZNhyhTIyYH09Eq9pLWWLzd9Sat6rUiMTWRb0TbW7VrHH9P+yIjOI2hWp1mlXk9EROR0KRSLhJL8fF8g3r//+GNHQnJGBixaBMnJFb7cpr2byMrLwpPnYdWOVTx++ePcf8n9DDpvEDd2uJEwo19WiYhIzaBQLBJKxo//3+jwyXi9MGECTJx4xpcpKS3h2unXMnvNbEptKZc2u5QHL3mQjPYZAISHhZ/xe4uIiFQFDdOIhJLs7PKF4qys03pbay3fbP6G575+DvCF3qa1m/JQj4dYc+caPh35KSNTRhIXGXemlYuIiFQprVMsEkqKiir1vILCArIXZePJ87Dsh2XEumIZ3nk4dWrVYVK/SRUoVEREpHppnWKRUBJXzpHacpyXsyyHpAlJ3PfxfSREJTDpmklsvnszdWrVqWCRIiIi1U9zikVCydChvlUmTjWFwuWCYcN+0mStZeHWhXhyPaQ1T2Ngu4Fccs4l3Nv9XkamjOTcBudWceEiIiJVS6FYJJSMG+dbdu3nQvHYsQD8sO8HXlv8Gu5cN4u2LSIyPJJGsY2gHTSOa8xfe/+1mgoXERGpWgrFIqEkOdm3DvGx6xSDLwy7XNg338T4l2Pr+1pfviv4jq5NuvL8Vc8zuMNg6kXXc6h4ERGRqqPVJ0RCTXq6bx3i0aN/sqPd4tuvZ5x7MK3X/Iaiw74P2k24cgKLb1/MN6O+YUzXMWceiB3eQU9EROTnKBSLhKLkZJg4kd3bvmfil8+S+mRrOtWbzt9XZtG5cWd2HdgFwKXNLqVDYoeKXWvWLOjUyTeXubAQrP3fDnqdOvmOi4iIOEyhWCTEFJcWs/PATgDW717PnbPupNSW8mzfZ9kybgtvDXqLpISkyrlY2R30jp3H7PX62jMyNGIsIiKO05xikRCx/IfleHI9ZC3K4orkK5jSfwopjVNYNmYZ7Rq2q5qLVtMOeiIiIhWlkWKRIPf64tf5xeRf0P6F9oz/cjxdm3ZlUPtBR49XWSCGKttBT0REpLJpRzuRIFNSWsK89fPo2aInYSaMBVsWsM+7j/F9xnNzx5tpFNeo+oqp5B30REREqop2tBMJEqt3rOahTx6i+bPN6Z3Vm7nr5gLwl8v/wqLbFnH3RXdXbyCGSt1BT0REpCpp+oRIgNu8dzM93D1oM7ENj//ncTomdmRGxgwuOecSAKIiojDGOFPc0KG+9Y9P5QQ76ImIiFQ3fdBOJMCU2lI+Xf8pOw7sIKN9Bo3iGhERFsHjlz/O0E5DaRrf1OkS/+c0d9ATERFxikKxSIBYt2sdU/Om4snzsH73eto3bM/17a4nIiyCf4/4t9PlnVg5dtAjJ8d3noiIiIM0fUIkAPzp0z/R8rmW/PHTP9KqXiteG/ga80fNd25axOk4yQ56jB7ta09Pd7pCERERjRSL1DTWWr7Y8AXuhW7uvfhe2jZoS8/mPQnvGc6wzsM4J+Ecp0s8ff4d9LQWcYDJz/etNZ2d7VshJC7ON0983DiN7otI0FEoFqkhNu7ZyNS8qfzj23+w+bPNxEXG0bdVX9o2aEuPZj3o0ayH0yVKKJk16/hpL0e2554yxTftRaP8IhJEFIpFHGStxRjDfu9+2j3fjn3efaQkpPBYn8e4vv31xEVqqTJxQNntuY91JCRnZPimv2jEWESChEKxSDWz1vL15q9xL3Szdvda5gybQ4wrBk9/D+efdT4b8jaQlpLmdJkSyrQ9t4iEIH3QTqSaFBQW8Lcv/kb7F9pz0SsXkb04mya1m3Co+BAAGe0zaFm3pcNViqDtuUUkJGmbZ5EqdLD4IKW2lBhXDB+s/oAHPnmAS865hMn9JnPDeTcQHxXvdIkix9P23CISgrTNs0gls9Yyf8t87vjwDpqMb8Ir370CwOAOg1l1xyo+z/ycW8+/VYFYai5tzy0iIUhzikUqibWWZ756hldzX2XJ9iVEhUcxsN1AUpukAhAXGUfr+q0drlKkHIYO9a0y8XM7EWp7bhEJIppTLFIBh0sO89WmrwAwxvDuyneJccXwj6v/wdZ7tvL69a9zUdJFDlcpcprGjfOF3lPR9twiEmQ0UixyBvK25uHJ9ZC9OJvdB3ez+e7NJMYmMuvmWcS4YpwuT6RitD23iIQgjRSLnIavN31Nl5e6kPJSCi/Mf4G05mm8N/g96kXXA1AgluCh7blFJMRopFjkFIpLi5m9Zjb1o+tzUdJFNIprRERYBH9P/zs3dbiJ+jH1nS5RpOpoe24RCSEKxSInsHT7UqbkTSFrURZbi7YypOMQLkq6iOZ1mvPtqG+dLk9EREQqmUKxyDEG5wzmjaVvEBEWwdWtryYzJZOrWl/ldFkiIiJShRSKJaSVlJYwZ+0c3lj6Bi9e/SJREVH0Se5Dt6bduLnTzSTGJjpdooiIiFQDhWIJSat2rMKT62Fq3lQ2F26mXnQ97up2FymNU7ilyy1OlyciIiLVTKFYQk7e1jxSXkohzISR3iqdZ/s+yzVtriEqIsrp0kRERMQhCsUS1EptKf9e9288eR4SYxIZf+V4OjXqxAtXvUD/tv05q/ZZTpcoIiIiNYBCsQSltbvW4sn1MCVvChv2bCAhKoHbUm8DfDvP3d71docrFBERkZpEoViCRtHhImJdsRhjeOI/TzBpwST6JPfhb73/xnXnXke0K9rpEkVERKSGcnRHO2NMP2PMpD179jhZhgQway2frv+UzPcyafxUY77a9BUAD/V4iA1jNzB76GwGdxisQCwiIiKn5OhIsbV2JjAzNTV1lJN1SOApPFTIM189gyfPw9pda6kdWZubOtx0dLvlpIQkhysUERGRQKLpExIw9nv3s27XOs5LPA9XuItnvn6Gzo0688hljzCw3UBiI2OdLlFEREQClEKx1GjWWr7c9CXuhW7eWPoGjeIaseqOVdSKqMXa36wloVaC0yWKiIhIEFAolhrr7eVv87tPfsfKHSuJccVwQ/sbyEzJPHpcgVhEREQqi0Kx1BgHiw/y3or36J7UnaSEJAyGxNhE7r/4fjLaZ1A7qrbTJYqIiEiQUigWR1lrmb9lPu5cN9OWTGP3wd080fsJ7r34Xga0G8CAdgOcLlFERERCgEKxOKa4tJiuL3cld2sutSJqMbDdQDJTMunVopfTpYmIiEiIUSiWanO45DD/XPVP5m+Zz18u/wsRYRFc1eoqxqSOYdB5gzRHWERERByjUCxVbmHBQty5bl5f/Do7Duygae2mPHjJg9SOqs1jlz/mdHkiIiIiCsVStabmTWXEuyOIDI+kf9v+ZKZkckXLKwgPC3e6NBEREZGjFIql0nhLvMxaMwt3rpsBbQcwvPNwrmp9Fc9f9TyDOww+utuciIiISE2jUCwVtnjbYjy5HrIXZ7N933YSYxO5MvlKABrENGBM1zEOVygiIiJyagrFckYOeA8Q7YoG4Jb3byF3ay792vQjMyWTvq364gp3OVyhiIiISPkpFEu5FZcW81H+R7hz3Xy89mPW37WehFoJvHLtK5wVdxYNYxs6XaKIiIjIGVEolp+1ee9mnvv6ObIWZVFQVED96PoM7zScQyWHAOjUqJPDFYqIiIhUTJjTBUjNtOfgHjbs2eB7fGgP478cT2qTVN4e9DZbxm3h2fRnSYxNdLjKIJKfD2PGQHw8l/XqBfHxvuf5+U5XJiIiEhI0UixHlZSWMHfdXNy5bt5Z8Q792/Zn2vXTaN+wPdvu2Ub9mPpOlxicZs2CjAzwesHrxQAUFsLkyTBlCuTkQHq601WKiIgENYViAeCZr57h6S+fZuPejdSpVYdbUm7hli63HD2uQFxF8vN9gXj//uOP+UMyGRmwaBEkJ1d/fSIiIiFC0ydCVOGhQqbmTaW4tBiAHft30CGxA29kvEHBuAKev/p5LmhygcNVhoDx433B91S8XpgwoXrqERERCVEaKQ4hpbaUz77/DHeum5xlOez37qdJ7Sb0btmbP/X8E8YYp0sMPdnZ5QvFWVkwcWL11CQiIhKCHA3Fxph+QL9WrVo5WUZIWL97PT2n9GT97vXER8UztONQRqaM5Bdn/wJAgdgpRUWVe56IiIicEUenT1hrZ1prRyckJDhZRlDad3gf/9r6LyYtmARAUnwSF519EdkDsikYV8BL/V7ioqSLFIadFhdXueeJiIjIGdH0iSBireU/G/+De6GbGctmUHS4iO77uzP6gtGEh4Xz+vWvO12iHGvoUN8qE6eaQuFywbBh1VeTiIhICFIoDiL3fHQPT3/1NLGuWAadN4jOtjO/ue43TpclpzJunG/ZtZ8LxWPHVl9NIiIiIUirTwSoA94DTFs8jT5ZfcjbmgfAkI5D8FznYes9W3n1ulfpXKezpkfUdMnJvnWIY2J84bcsl8vXnpOj5dhERESqmEaKA4i1lq83f40n18P0JdPZc2gPzRKasaVwC50bd+aCJhdoGbVAlJ7uW4d4wgTIysIWFmJq1/ZNmRg7VoFYRESkGigUBwBviRdXuIv93v30ntqbUltKRvsMRqaMJK15GmFGA/4BLznZt+TaxIl8Om8eaWlpTlckIiISUhSKa6hDxYeYuWom7lw3W4u2Mn/UfGIjY/lgyAd0OasL8VHxTpcoIiIiEjQUimuY5T8s54VvX+D1Ja+z88BOmtZuyojOI/CWeokMj+Sy5pc5XaKIiIhI0FEorgG279tOZHgkdWrVYUHBAl7+7mUGtBtAZkoml7e4nPCwcKdLFBEREQlqCsUO8ZZ4+WD1B7hz3Xy4+kMev/xxxnUfR0b7DK5ufTV1o+s6XaKIiIhIyFAormbWWu756B6yFmXxw/4faBzXmLt/cTf9zu0HQK2IWtSKqOVwlSIiIiKhRaG4Gvy4/0e+2PAF/dv2xxjDyh0rubTZpWSmZHJlqyuJCFM3iIiIiDhJaayKFJcWM3vNbNy5bmaunEmJLaFgXAGJsYm8f9P7WkZNREREpAZRKK4C/173b4a8PYStRVtpENOAX3f9NZldMkmMTQRQIBYRERGpYRSKK8GuA7uYvmQ6reu3pnfL3rSp34ZuTbsxMmUkV7W+isjwSKdLFBEREZFTUCg+QyWlJXy89mPcuW7eXfEuh0oO8asLfkXvlr1pGt+Udwe/63SJIiIiIlJOCsVnqO9rffl47cfUi67HqPNHkdklky6NuzhdloiIiIicAYXicth7aC8zls7gzWVv8s6N7xDjimFM6hh+dcGv6NemH1ERUU6XKCIiIiIVoFB8EqW2lH+v+zeePA9vLXuLA8UHaNegHet3r6d9w/YMaDfA6RJFREREpJIoFB+jpLSE8LBw8rbm0TurNwlRCYzoPIKRKSO5sOmFGGOcLlFEREREKplCMVB0uIicZTm4c920rteayddOJqVxCu8Nfo8rWl5BtCva6RJFREREpAqFdCj+atNXTFowiRlLZ7DPu49W9VoxoK1vWoQxhmvPvdbhCkVERESkOoR0KJ6+ZDpvLnuTwR0Gk5mSSfek7poeISIiIhKCQjoUP3zpwzzW6zFiI2OdLkVEREREHBTSobhBTAOnSxARERGRGiDM6QJERERERJymUCwiIiIiIU+hWERERERCnkKxiIiIiIQ8hWIRERERCXkKxSInk58PY8ZAfDyEhfn+HDPG1y4iIiJBpdJDsTEm1hgzxRjzsjHm5sp+f5FqMWsWdOoEkydDYSFY6/tz8mRf+6xZTlcoIiIilahcodgY86oxZrsxZskx7X2NMSuNMWuMMQ/4mwcCOdbaUYD2SZbAk58PGRmwfz94vT895vX62jMyNGIsIiISRMo7UuwB+pZtMMaEA88D6UB74CZjTHvgbGCj/7SSyilTpBqNH398GD6W1wsTJlRPPSIiIlLlyhWKrbWfATuPab4QWGOtXWutPQxMB64DNuELxuV+f5EaJTu7fKE4K6t66hEREZEqV5FtnpvyvxFh8IXhbsBzwERjzNXAzJO92BgzGhgN0KhRI+bNm1eBUuREioqKdF/PwGVFRZhynGcLC/m0Cu6v+i1wqe8Cl/oucKnvAlNN7LeKhOITstbuAzLLcd4kYBJAamqqTUtLq+xSQt68efPQfT0DcXG+D9X9DFO7dpXcX/Vb4FLfBS71XeBS3wWmmthvFZnesBlIKvP8bH+bSGAbOhRcrlOf43LBsGHVU4+IiIhUuYqE4m+B1saYFsaYSGAw8H7llCXioHHjyheKx46tnnpERESkypV3SbZpwJfAucaYTcaYW621xcAdwL+A5cAMa+3SqitVpJokJ0NODsTEHB+OXS5fe06O7zwREREJCuWaU2ytvekk7R8CH1ZqRSI1QXo6LFrkW3YtKwuKinxzjYcN840QKxCLiIgElUr/oN3pMMb0A/q1atXKyTJETiw5GSZO9H2JiIhIUHN0HWFr7Uxr7eiEhAQnyxARERGREKfNNUREREQk5CkUi4iIiEjIUygWERERkZCnUCwiIiIiIa9GrD4B7DXGrHayliDVAPjR6SLktKnfApf6LnCp7wKX+i4wOdlvzU7UaKy11V2IVBNjzHxrbarTdcjpUb8FLvVd4FLfBS71XWCqif2m6RMiIiIiEvIUikVEREQk5CkUB7dJThcgZ0T9FrjUd4FLfRe41HeBqcb1m+YUi4iIiEjI00ixiIiIiIQ8heJqZIxJMsb82xizzBiz1Bhzl7+9njFmjjFmtf/Puv72m40xi4wxi40x/zXGdC7zXn2NMSuNMWuMMQ+c4poj/O+72hgzokz7Tf73XWSMmW2MaXCS15/wOsYYjzFmnTEm1/+VUhn3qCYKsn7rZYz5zhizxBgzxRjj6LKMVS1A++5VY8x2Y8ySY9of9b821xjzkTGmSUXvT00WZH33Rpl/K9cbY3Iren9qqkDrt5PV6z92g7+t1BhTo1ZJqApB1nePGGM2l/m5u6pcN8Faq69q+gLOAs73P64NrALaA08AD/jbHwD+5n/cHajrf5wOfO1/HA7kAy2BSCAPaH+C69UD1vr/rOt/XBff+tTbgQb+854AHjnB6096HcADZDh9T9Vv5e83fP8TvBFo4z/vT8CtTt9f9d1x73EpcD6w5Jj2+DKPfwO86PT9Vd+Vr++OOWc88H9O31/126nr9T9vB5wLzANSnb636rvT6rtHgHtO9x5opLgaWWsLrLXf+R8XAsuBpsB1wBT/aVOA/v5z/mut3eVv/wo42//4QmCNtXattfYwMN3/Hse6Ephjrd3pf585QF/A+L9ijTEGiAe2nOD15b1OUAuifqsPHLbWrvKfNwe4/rRvSAAJwL7DWvsZsPME7XvLPI0FgvoDIcHUd0f4Xz8ImPbzdyAwBVq/naJerLXLrbUrz/hmBJhg6rszpVDsEGNMc6AL8DXQyFpb4D+0FWh0gpfcCszyP26Kb8TviE2c+BvhhOdZa73A7cBifN9o7YFXyvv6Ms8f8/9qY4IxJuoErw86Ad5vPwIRZX4NmAEkneD1QSlA+u7n/g6PGWM2AjcD/3e6rw9UwdB3fj2AbdbakNjBNdD67Zh6Q1qQ9N0d/ozy6pEpHz9HodgBxpg44C3gt8eM/mB94/72mPN74vuGu7+Sru/C9w3XBWgCLAIePM23eRBoC3TF96uPSqmtJgv0fvPXOBiYYIz5BigESiqjtpou0PvuCGvtQ9baJOA14I7KqK2mC5a+87uJIB4lLivQ+u1U9YaaIOm7fwDJQApQgG/a0s9SKK5m/s5+C3jNWvu2v3mbMeYs//Gz8M2lOXJ+J2AycJ21doe/eTM/HeE7G9hsjOlWZlL5tSc7D983CdbafP83+Aygu3/S+pHX33aK1x/5tYW11h4C3Ph+XRK0gqjfvrTW9rDWXgh8hm8OVlALsL4rr9cI8qkvEFx9Z3wfah0IvHGatyHgBFq/naTekBQsfWet3WatLbHWlgIvU96MYmvA5O5Q+cI3R2Yq8Mwx7U/y00nsT/gfnwOsAbofc34EvgnpLfjfJPbzTnC9esA6fBPX6/of18P3f14FQEP/eY8C40/w+pNeBzirzN/pGeBxp++v+q1c/Zbo/zMK+ATo5fT9Vd+dsO7mHP9Bu9ZlHt8J5Dh9f9V35es7f3tf4FOn76v6rXz1HnPOPELjg3ZB03f4M4r/8VhgernugdOdEEpfwCX4fu2wCMj1f12F7wNQnwCrgY+Bev7zJwO7ypw7v8x7XYVvlC8feOgU17zF/027Bsgs034bvknpi4CZQP2TvP6E1wHm4pvvswTIBuKcvr/qt3L125P+16/E96smx++v+u6410/D9x8EL745drf629/y/7wdeX1Tp++v+q58fec/5gFuc/q+qt/KV6//2AB/Px4CtgH/cvr+qu/K3XdZ+DLKIuB9yoTkU31pRzsRERERCXmaUywiIiIiIU+hWERERERCnkKxiIiIiIQ8hWIRERERCXkKxSIiIiIS8hSKRURERCTkKRSLiIiISMhTKBYRERGRkPf/WmKFdpQmzXAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(12, 5))\n",
    "ax.semilogy(df[\"date\"], df[\"new cases\"], \"r.\", ms=20)\n",
    "\n",
    "t_fit = np.array(df[6:].index)\n",
    "y_fit = np.log(np.array(df[6:][\"new cases\"]))\n",
    "p = np.polyfit(t_fit, y_fit, 1)\n",
    "\n",
    "t = np.linspace(min(df.index), max(t_fit) + 3)\n",
    "y = np.exp(p[0]*t + p[1])\n",
    "dates = [df[\"date\"][0] + timedelta(days=t) for t in t]\n",
    "\n",
    "print(p)\n",
    "ax.semilogy(dates, y, 'g--')\n",
    "ax.grid(True)\n",
    "ax.set_ylim(0.8, 150)\n",
    "ax.set_title('New cases')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 526,
     "status": "aborted",
     "timestamp": 1597709228633,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "nkJN6e5mkj-I",
    "nbpages": {
     "level": 3,
     "link": "[2.1.4.1 Simple logarithmic extrapolation](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4.1-Simple-logarithmic-extrapolation)",
     "section": "2.1.4.1 Simple logarithmic extrapolation"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.22259873535947472"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(p[0]) - 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.1.4.1 Simple logarithmic extrapolation](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4.1-Simple-logarithmic-extrapolation)",
     "section": "2.1.4.1 Simple logarithmic extrapolation"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpages": {
     "level": 3,
     "link": "[2.1.4.1 Simple logarithmic extrapolation](https://jckantor.github.io/CBE40455-2020/02.01-Campus-SIR-modeling.html#2.1.4.1-Simple-logarithmic-extrapolation)",
     "section": "2.1.4.1 Simple logarithmic extrapolation"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [2.0 Modeling](https://jckantor.github.io/CBE40455-2020/02.00-Modeling.html) | [Contents](toc.html) | [2.2 Campus SEIR Modeling](https://jckantor.github.io/CBE40455-2020/02.02-Campus-SEIR-modeling.html) ><p><a href=\"https://colab.research.google.com/github/jckantor/CBE40455-2020/blob/master/docs/02.01-Campus-SIR-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.01-Campus-SIR-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": "ABX9TyMlC91CUu1UroEDeoXT2chU",
   "collapsed_sections": [],
   "name": "Campus-SIR-modeling.ipynb",
   "provenance": []
  },
  "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
}