optimize_r.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "import gripql\n",
    "import pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "conn = gripql.Connection(\"http://192.168.1.147:8201\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "G  = conn.graph(\"covid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-04-17 09:42:48,293\t2 results received in 0 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<AttrDict({'gid': '41051', 'label': 'SummaryLocation', 'data': {'county': 'Multnomah County', 'projects': [], 'province_state': 'OR', 'submitter_id': '41051', 'type': 'summary_location', 'zipcodes': ['97292', '97205', '97292', '97206', '97292', '97216', '97292', '97223', '97292', '97266', '97292', '97239', '97292', '97227', '97292', '97203', '97292', '97208', '97292', '97282', '97292', '97133', '97292', '97010', '97292', '97233', '97292', '97220', '97292', '97236', '97292', '97034', '97292', '97210', '97292', '97212', '97292', '97214', '97292', '97218', '97292', '97290', '97292', '97294', '97292', '97254', '97292', '97293', '97292', '97009', '97292', '97296', '97292', '97250', '97292', '97232', '97292', '97231', '97292', '97221', '97292', '97258', '97292', '97030', '97292', '97019', '97292', '97211', '97292', '97213', '97292', '97215', '97292', '97014', '97292', '97283', '97292', '97240', '97292', '97251', '97292', '97080', '97292', '97060', '97292', '97204', '97292', '97222', '97292', '97207', '97292', '97252', '97292', '97230', '97292', '97024', '97292', '97124', '97292', '97253', '97292', '97056', '97292', '97217', '97292', '97280', '97292', '97228', '97292', '97219', '97292', '97229', '97292', '97035', '97292', '97209', '97292', '97202', '97292', '97201', '97292', '97225', '97292', '97286', '97292', '97238', '97292', '97242']}})>,\n",
       " <AttrDict({'gid': '41067', 'label': 'SummaryLocation', 'data': {'county': 'Washington County', 'projects': [], 'province_state': 'OR', 'submitter_id': '41067', 'type': 'summary_location', 'zipcodes': ['97078', '97133', '97078', '97008', '97078', '97225', '97078', '97075', '97078', '97124', '97078', '97035', '97078', '97229', '97078', '97070', '97078', '97006', '97078', '97116', '97078', '97007', '97078', '97132', '97078', '97231', '97078', '97298', '97078', '97062', '97078', '97113', '97078', '97005', '97078', '97291', '97078', '97210', '97078', '97056', '97078', '97223', '97078', '97119', '97078', '97125', '97078', '97221', '97078', '97219', '97078', '97224', '97078', '97117', '97078', '97106', '97078', '97140', '97078', '97003', '97078', '97123', '97078', '97109', '97078', '97144', '97078', '97281']}})>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G.query().V().hasLabel(\"SummaryLocation\").has(gripql.contains(\"zipcodes\", \"97229\")).execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-04-17 09:42:48,309\t23 results received in 0 seconds\n"
     ]
    }
   ],
   "source": [
    "q = G.query().V(\"41067\").out(\"summary_reports\").render([\"date\", \"confirmed\", \"deaths\", \"recovered\"])\n",
    "\n",
    "data = {}\n",
    "for row in q:\n",
    "    d = datetime.datetime.strptime(row[0], \"%Y-%m-%d %H:%M:%S\")\n",
    "    data[d] = {\"confirmed\":int(row[1]), \"deaths\" : int(row[2]), \"recovered\":int(row[3])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pandas.DataFrame(data).transpose().sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = pandas.Series( (df.index - df.index[0]).round(\"D\").astype('timedelta64[D]'), index=df.index, name=\"days\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddf = df.join(delta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>confirmed</th>\n",
       "      <th>deaths</th>\n",
       "      <th>recovered</th>\n",
       "      <th>days</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2020-03-23 23:19:34</th>\n",
       "      <td>69</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-24 23:37:31</th>\n",
       "      <td>76</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-25 23:33:19</th>\n",
       "      <td>96</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-26 23:48:35</th>\n",
       "      <td>104</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-27 22:14:55</th>\n",
       "      <td>122</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-28 23:05:00</th>\n",
       "      <td>140</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-29 23:08:00</th>\n",
       "      <td>154</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-30 22:52:00</th>\n",
       "      <td>168</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-03-31 23:43:56</th>\n",
       "      <td>186</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-01 21:58:49</th>\n",
       "      <td>189</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-02 23:25:00</th>\n",
       "      <td>211</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-03 22:46:37</th>\n",
       "      <td>228</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-04 23:34:00</th>\n",
       "      <td>228</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-05 23:06:45</th>\n",
       "      <td>259</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-06 23:22:00</th>\n",
       "      <td>259</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-07 23:04:49</th>\n",
       "      <td>271</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-08 22:51:58</th>\n",
       "      <td>282</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-09 23:02:37</th>\n",
       "      <td>311</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-10 22:54:07</th>\n",
       "      <td>311</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-11 22:45:33</th>\n",
       "      <td>320</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-12 23:18:00</th>\n",
       "      <td>348</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-13 23:07:54</th>\n",
       "      <td>360</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2020-04-14 23:33:31</th>\n",
       "      <td>365</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     confirmed  deaths  recovered  days\n",
       "2020-03-23 23:19:34         69       1          0     0\n",
       "2020-03-24 23:37:31         76       2          0     1\n",
       "2020-03-25 23:33:19         96       2          0     2\n",
       "2020-03-26 23:48:35        104       3          0     3\n",
       "2020-03-27 22:14:55        122       3          0     4\n",
       "2020-03-28 23:05:00        140       3          0     5\n",
       "2020-03-29 23:08:00        154       3          0     6\n",
       "2020-03-30 22:52:00        168       3          0     7\n",
       "2020-03-31 23:43:56        186       3          0     8\n",
       "2020-04-01 21:58:49        189       3          0     9\n",
       "2020-04-02 23:25:00        211       4          0    10\n",
       "2020-04-03 22:46:37        228       4          0    11\n",
       "2020-04-04 23:34:00        228       4          0    12\n",
       "2020-04-05 23:06:45        259       4          0    13\n",
       "2020-04-06 23:22:00        259       4          0    14\n",
       "2020-04-07 23:04:49        271       5          0    15\n",
       "2020-04-08 22:51:58        282       6          0    16\n",
       "2020-04-09 23:02:37        311       6          0    17\n",
       "2020-04-10 22:54:07        311       6          0    18\n",
       "2020-04-11 22:45:33        320       6          0    19\n",
       "2020-04-12 23:18:00        348       6          0    20\n",
       "2020-04-13 23:07:54        360       7          0    21\n",
       "2020-04-14 23:33:31        365       7          0    22"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ddf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-04-17 09:42:48,358\t25 results received in 0 seconds\n"
     ]
    }
   ],
   "source": [
    "q = G.query().V(\"41067\").out(\"census\").has(gripql.eq(\"gender\", None)).render([\"population\"])\n",
    "population = sum(list(a[0] for a in q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'vertex_labels': ['CensusReport', 'SummaryLocation', 'SummaryReport'],\n",
       " 'edge_labels': ['census', 'census_of', 'performed_at', 'summary_reports']}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G.listLabels()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe10cb0f400>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xW5f3/8deHkBBmGAkQQiAsWTINQ62KWLcWtVVxIThAxWpb68/RodYO22qtVWuLg1FRQUAL7lm1DiDsDQGCIQQIBAJJyLyv7x859hcZkn3u8X4+Hvcj577OOcknx5O3h+u+znXMOYeIiISXRn4XICIidU/hLiIShhTuIiJhSOEuIhKGFO4iImGosd8FAMTHx7uUlBS/yxARCSlLlizZ45xLONq6oAj3lJQU0tLS/C5DRCSkmNm2Y61Tt4yISBhSuIuIhCGFu4hIGFK4i4iEIYW7iEgYUriLiIQhhbuISBgKinHuIiKRpKQswOodeSzamsuApDhO7Rlf5z9D4S4iUs8OlZSzLHMfi7bmsmhrLku/3kdRaQCAW0f1ULiLiISCA0WlLMnYx8KtuSzaupdVWXmUljvMoF9iK64a3oXhKW0Z1q0t8S2a1EsNCncRkVrKKyzli817WLg1l8UZuazNPoBzEB1lDOzcmptO687wlLYM7dqGuKbRDVKTwl1EpIZKygJM/yKDv324iYPFZcRGN+Kkrm2486xeDO/WliHJbWgaE+VLbQp3EZEa+Gj9Lh5+Yx1b9xRwZu8EJp/Zk4GdWxPTODgGISrcRUSqIX13Pr99cy3/2ZBD9/jmTB0/jDP7tPe7rCMo3EVEqiDvUClPfriJaV9k0DQ6il9e2JdxJ6cEzZX64Y4b7mYWC3wKNPG2n+Oce8DMpgFnAHnepuOdc8vNzIAngAuAQq99aX0ULyJS38oDjtlpmTz67gZyC0sYOyyZu87pXW+jXOpKVa7ci4HRzrl8M4sG/mtmb3vr7nbOzTls+/OBXt5rBPCM91VEJKQs2prLQwvWsGbHAYaltGH6xcM5MSnO77Kq5Ljh7pxzQL73Ntp7ue/YZQwww9vvKzNrbWaJzrnsWlcrItIAsvYf4pG317NgxQ46xcXy5FVDuGhgIhUdE6GhSn3uZhYFLAF6Ak875xaa2a3A78zs18CHwL3OuWIgCcistPt2ry37sO85EZgI0KVLl9r+HiIitVZUWs4/PtnMPz7ZjHNw51m9uOWMHr4NZ6yNKoW7c64cGGxmrYHXzOxE4D5gJxADTAHuAX5T1R/snJvi7Udqaup3/UtARKTelQccE6Yu5sste7loYCL3XdCXpNZN/S6rxqr1Ma9zbj/wMXCecy7bVSgGpgLDvc2ygORKu3X22kREgtaTH23iyy17eeSyATx19dCQDnaoQribWYJ3xY6ZNQXOBtabWaLXZsAlwGpvl/nAOKswEshTf7uIBLMvN+/lbx9u4tIhSVw5LPn4O4SAqnTLJALTvX73RsBs59wbZvaRmSUABiwHbvG2f4uKYZDpVAyFnFD3ZYuI1I29+cXc+coyUto15+FLTgypD02/S1VGy6wEhhylffQxtnfA5NqXJiJSvwIBx12vrmD/oVKmThhGiybhc19ncN5aJSLSAJ777xb+syGHX13Yl/6dQmP8elUp3EUkIi39eh9/emcD55/YkWtHdvW7nDqncBeRiJN3qJQ7Xl5Gx7hYHvnhwLDpZ68sfDqYRESqwDnHvXNXsjOviFdvObnBHp7R0HTlLiIR5cWvtvH26p3cfW5vhnRp43c59UbhLiIRY82OPB5+cx2jeidw82nd/S6nXincRSQiFBSX8eOXltGmWTSPXT6IRo3Cr5+9MvW5i0hE+NXrq8nYW8BLN4+kXZDPxV4XdOUuImFvzpLtzFuWxR1n9WJk93Z+l9MgFO4iEtbSdx/kV6+vZmT3tvx4dC+/y2kwCncRCVtFpeXc/tIymsZE8cTYIUSFeT97ZepzF5Gw9Zs31rJ+50GmTRhGh1axfpfToHTlLiJh6c2V2by08GsmndGdUb3b+11Og1O4i0jYWZG5n3vnrmRIl9b8/JzefpfjC3XLiEhYKA843l+7kxc+z2DR1lzaNo/hb2OHEB0VmdewCncRCWl5h0qZtfhrpn+xjaz9h0hq3ZRfXNCXK4Ylh+28MVWhcBeRkLQ5J59pn2cwd+l2CkvKGd6tLb+6qC/f79uBxhF6tV6Zwl1EQkYg4Ph0Uw5TP8/gk405xEQ14geDOzH+lBROTAqvh23UlsJdRIJeQXEZ85ZuZ9oXGWzOKSChZRN++v0TuHpEFxJahv9UAjWhcBeRehcIOBas3EHOweJq77tjfxFzlmRyoKiMgZ3jePzKQVw4oBMxjdX18l0U7iJS7579bAt/eHt9jfaNamScd2JHbjg1haFd2oTlU5Pqg8JdROrV6qw8Hn1vA+f278Cjlw+q9v7RUY2IjY6qh8rCm8JdROpNYUkZd7y8jHbNm/DIZQNpGRu5QxMbmsJdROrNw2+sY+veAmbeOII2zWP8LieiHPcTCTOLNbNFZrbCzNaY2UNeezczW2hm6WY2y8xivPYm3vt0b31K/f4KIhKM3lm9k5cXfc3E07tzSs94v8uJOFX5uLkYGO2cGwQMBs4zs5HAH4HHnXM9gX3Ajd72NwL7vPbHve1EJILszCvi3nkrGZAUx11nR+bcLn47bri7Cvne22jv5YDRwByvfTpwibc8xnuPt/4s08fbIhEjEHDc9epyiksD/HXsYA1Z9EmVjrqZRZnZcmA38D6wGdjvnCvzNtkOJHnLSUAmgLc+D4iM51qJCM9+toXP0/fywMX96JHQwu9yIlaVwt05V+6cGwx0BoYDfWr7g81sopmlmVlaTk5Obb+diASBb4Y9nte/I1cOS/a7nIhWrX8vOef2Ax8DJwOtzeyb0TadgSxvOQtIBvDWxwF7j/K9pjjnUp1zqQkJCTUsX0SCxbeGPf5wgG428llVRsskmFlrb7kpcDawjoqQ/5G32fXAv73l+d57vPUfOedcXRYtIsHnm2GPf7lyEK2badij36oyzj0RmG5mUVT8z2C2c+4NM1sLvGJmvwWWAc972z8P/MvM0oFcYGw91C0iQeSbYY+TzujOKT007DEYHDfcnXMrgSFHad9CRf/74e1FwOV1Up2IBD0NewxOGqMkIjWmYY/BS/8lRKTGNOwxeCncRaRGNOwxuCncRaTaNOwx+GlWSBGptv/N9njTCA17DFK6cheRavnfsMfTe2jYYxDTlbuIVMnBolKe+2wrUz7dwoCkOH529gl+lyTfQeEuIt+pqLScF7/axtMfp7OvsJQLByTyy4v6athjkFO4i8hRlZUHmLcsi7++v5EdeUWc1iueu8/tzcDOrf0uTapA4S4i3+Kc4901u3j0vQ2k785nUOc4/nz5IE7V05RCisJdRP7ni817+OM7G1iRuZ8eCc35x7VDObd/Rw11DEEKdxFh1fY8/vTuej7btIfEuFj+9MOBXDY0icZR6lcPVQp3kQi2JSefx97fyJsrs2nTLJpfXtiXa0d2JTY6yu/SpJYU7iIRaG9+MY+9v5FZizNp0rgRd4zuyU2nd6dVbLTfpUkdUbiLRJDS8gAzvtzGXz/YyKGScq4d0YXbR/cioWUTv0uTOqZwF4kQ/9mwm4ffWMvmnAJO6xXPAxf3o2f7ln6XJfVE4S4S5rbk5PPbN9fx0frdpLRrxvPXpzK6T3uNgAlzCneRMHWgqJSnPkpn6udbadI4ivvO78P4U1No0lgflkYChbtImAkEHK8uyeTP725gb0EJl5/UmZ+f25v2LWP9Lk0akMJdJIykZeTy0IK1rMrK46SubXhh/DBNFxChFO4iYSA77xB/eGs981fsoGOrWJ4YO5gfDOqkfvUIpnAXCRJLtu3j6Y/TKQu4au3nnCMtYx8B57hjdE9uGdWDZjH60450OgNEgkDeoVJum7mEsnJHcttm1d7//BM78tOzT6jRvhKeFO4iQeDhN9ayJ7+EebeewqBk9ZFL7WlWIBGffbR+F3OWbGfS6d0V7FJnjhvuZpZsZh+b2VozW2Nmd3rtD5pZlpkt914XVNrnPjNLN7MNZnZuff4CIqEsr7CU++at4oQOLbjz+738LkfCSFW6ZcqAu5xzS82sJbDEzN731j3unHu08sZm1g8YC/QHOgEfmNkJzrnyuixcJBz8xuuOeXZcqm4ukjp13Ct351y2c26pt3wQWAckfccuY4BXnHPFzrmtQDowvC6KFQknH63fxdyl27nljO4aiy51rlp97maWAgwBFnpNt5vZSjN7wczaeG1JQGal3bZzlP8ZmNlEM0szs7ScnJxqFy4SyvIKS7l37ip6d2jJHWepO0bqXpXD3cxaAHOBnzjnDgDPAD2AwUA28Fh1frBzbopzLtU5l5qQkFCdXUVC3kNvrGFvQQmPXj5I3TFSL6oU7mYWTUWwz3TOzQNwzu1yzpU75wLAs/z/rpcsILnS7p29NhEBPli7i3lLs7htVA8GdI7zuxwJU1UZLWPA88A659xfKrUnVtrsUmC1tzwfGGtmTcysG9ALWFR3JYuErrzCUu5/bRV9Orbkx6PVHSP1pyqjZU4FrgNWmdlyr+1+4CozGww4IAOYBOCcW2Nms4G1VIy0mayRMiIVHlpQ0R3zwvhhxDTWbSZSf44b7s65/wJHm33ore/Y53fA72pRl0jY+WDtLuYty+KO0T05MUndMVK/dOkg0gD2F5Zwn9cdc7u6Y6QBaG4ZkQbw0IK17CsoYaq6Y6SB6CwTqWfvr93Fa8uyuO1MdcdIw1G4i9Sj/YUl3P/aKvomtuL2M3v6XY5EEHXLiNSjB+evYV9BCdMmqDtGGpbONpF68t6anby+fAeTz+xJ/07qjpGGpXAXqQf7Ckq4/7XV9E1sxWR1x4gP1C0j8h0CAUfAVe+ZpgAPLljD/sISpt+g7hjxh8Jd5BgWrNjBffNWkV9cVqP9f/L9XuqOEd8o3EUOEwg4Hnt/A09/vJmTurZh1AnVn7W0XYsmXJ7auR6qE6kahbtIJfnFZfzkleV8sG4XV6Ym8/AlJ6pbRUKSwl3E8/XeQm6asZjNOQU8eHE/rj8lhYpJUUVCj8JdBPhi8x5um7kU52D6hOF8r1e83yWJ1IrCXSKac44Xv9rGgwvW0i2+Oc+NSyUlvrnfZYnUmsJdIlZJWYAHF6zhpYVfM7pPe54YO5iWsdF+lyVSJxTuEpH25hdz28ylLNyay62jevDzc3oT1Uj96xI+FO4ScdZlH+DmGWnkHCzmibGDGTM4ye+SROqcwl0iyjurd/Kz2ctpGduY2ZNOZlBya79LEqkXCneJCM45nvwonb+8v5HBya2Zct1JtG8V63dZIvVG4S5hr7Q8wE9nLeeNldlcNiSJ3182gNjoKL/LEqlXCncJa+UBx12zV/DGymzuOa8Pt5zRXTcmSURQuEvYcs7xy9dXM3/FDu45rw+3jurhd0kiDUaTZkhYcs7xyNvreXnR19w2qoeCXSKOwl3C0tMfp/PPT7cw7uSu3H1ub7/LEWlwxw13M0s2s4/NbK2ZrTGzO732tmb2vplt8r628drNzP5mZulmttLMhtb3LyFS2dTPt/Loexu5bEgSD17cX33sEpGqcuVeBtzlnOsHjAQmm1k/4F7gQ+dcL+BD7z3A+UAv7zUReKbOqxY5hlfTMnlowVrO6deBP/1oII1016lEqOOGu3Mu2zm31Fs+CKwDkoAxwHRvs+nAJd7yGGCGq/AV0NrMEuu8cpHDvL0qm3vmruR7PeN58uohNI5Sr6NErmqd/WaWAgwBFgIdnHPZ3qqdQAdvOQnIrLTbdq/t8O810czSzCwtJyenmmWLfNsnG3O445VlFTcojTuJJo01jl0iW5XD3cxaAHOBnzjnDlRe55xzQLWeIuycm+KcS3XOpSYkVP8xZiLfWJyRy6R/pdGrfUumThhOsxiN8BWpUribWTQVwT7TOTfPa971TXeL93W3154FJFfavbPXJlLnVmflccPUxXRq3ZQZNw4nrqmm7BWBqo2WMeB5YJ1z7i+VVs0HrveWrwf+Xal9nDdqZiSQV6n7RqTOpO8+yLgXFtGqaTQv3jiC+BZN/C5JJGhU5d+vpwLXAavMbLnXdj/wCDDbzG4EtgFXeOveAi4A0oFCYEKdViwCZOYWcu1zi2hkxos3jaBT66Z+lyQSVI4b7s65/wLHGk921lG2d8DkWtYlcky7DhRxzXMLOVRazqxJI+mmx+KJHEFjxSSk7Cso4brnF7Inv5hpE4bRp2Mrv0sSCUoaViBBKb+4jIw9BWzZU0CG99q6t4D03fkUlwWYNmEYQ7q08btMkaClcBffFJaUkbGnkIy9BWz9JsT3FrB1TyF78ou/tW1iXCwp7Zpz0cBOXDokieHd2vpUtUhoULhLgyssKePuV1fy5qpvD6JKaNmEbu2ac1af9qTEN6dbfDNS4pvTtW1zmsbopiSR6lC4S4PKOVjMTdMXsyorj5tP68ag5NaktGtOSnxzWjTR6ShSV/TXJA0mffdBxk9dzN78Ev55XSpn9+tw/J1EpEYU7tIgvty8l0n/SiOmcRSzJo1kYOfWfpckEtYU7lLvXlu2nf83ZyVd2zVn6vhhJLdt5ndJImFP4S71xjnHUx+l89j7GxnZvS3/vDaVuGaa+0WkISjcpV6Ulgf4xWurmJ22ncuGJPHIDwcS01j3zIk0FIW71LkDRaVMnrmUzzbt4Y6zevHT7/fSo+5EGpjCXerUjv2HuGHaYtJ35/PnHw3k8tTk4+8kInVO4S51ZnVWHjdMW8yhknKmTRjO93rF+12SSMRSuEud+HjDbm6fuZS4ptHMufUUends6XdJIhFN4S61NnPhNn797zX0TWzJ89cPo0OrWL9LEol4CnepsfziMn7/1jpeWvg1Z/ZO4Kmrh9JcUwiIBAX9JUqNfLl5L3fPWUHW/kNMOr07d5/bm8ZRGuooEiwU7lIth0rK+dO765n6eQZd2zVj9qSTGZai6XdFgo3CXapsybZ9/PzVFWzdU8D1J3flnvP70CxGp5BIMNJfphxXUWk5j3+wkWc/3UJiXFNeumkEp/TUMEeRYKZwl++0cvt+7pq9gk2787lqeDL3X9CXlrGaH0Yk2Cnc5ahKygI8+dEm/v6fzSS0aMK0CcMY1bu932WJSBUp3OUI67IP8LPZK1iXfYDLhibxwMX9iWuqq3WRUKJwl/8pKw/wj08288SHm4hrGsOz4/S0JJFQpXAXoOIReHfNXsGK7XlcPKgTv/lBf9o0j/G7LBGpoePedWJmL5jZbjNbXantQTPLMrPl3uuCSuvuM7N0M9tgZufWV+FSN8oDjimfbuaCv/2Xr3MLefrqoTx51RAFu0iIq8qV+zTgKWDGYe2PO+cerdxgZv2AsUB/oBPwgZmd4Jwrr4NapY5t3VPAz19dwZJt+zinXwd+d+kAElo28bssEakDxw1359ynZpZSxe83BnjFOVcMbDWzdGA48GWNK5Q6Fwg4pn+ZwR/fWU9MVCMev3IQlwxO0gM1RMJIbfrcbzezcUAacJdzbh+QBHxVaZvtXtsRzGwiMBGgS5cutShDqiMzt5C756zgqy25jOqdwCOXDaRjnGZxFAk3NZ3p6RmgBzAYyAYeq+43cM5Ncc6lOudSExISaliGVJVzjpkLt3HeXz9lddYB/vjDAUwdP0zBLhKmanTl7pzb9c2ymT0LvOG9zQIqP1ets9cmPtqx/xD3zF3JZ5v2cGrPdvzxhwPp3KaZ32WJSD2qUbibWaJzLtt7eynwzUia+cBLZvYXKj5Q7QUsqnWVUiPOOeYs2c5vFqylLOB4eEx/rhnRlUaN1LcuEu6OG+5m9jIwCog3s+3AA8AoMxsMOCADmATgnFtjZrOBtUAZMFkjZfyx+0AR981bxYfrdzM8pS1/vnwgXds197ssEWkg5pzzuwZSU1NdWlqa32WEjfkrdvDrf6/mUEk5d5/bmxtO7aardZEwZGZLnHOpR1unO1TDSCDg+MPb63j2s60MTm7NY1cMokdCC7/LEhEfKNzDRFFpOT+dtZy3V+9k3Mld+fVF/fTYO5EIpnAPA3vzi7lpRhrLM/fzywv7cuP3uumGJJEIp3APcVty8pkwbTE784p45pqhnHdiot8liUgQULiHsMUZudw8I41GZrw8cSRDu7TxuyQRCRIK9xC1YMUO7pq9gs5tmjJ1wjANcxSRb1G4hxjnHP/4ZAt/fGc9w1LaMOW6VE3PKyJHULiHkLLyAL+ev4aXFn7NxYM68ecfDSQ2OsrvskQkCCncQ0R+cRmTZy7lk4053DaqBz8/p7duTBKRY1K4h4CdeUXcMG0xG3Yd5PeXDuDqEZoiWUS+m8I9yK3LPsAN0xZz4FApz1+fyqje7f0uSURCgMI9iH22KYdbX1xK8yZRzL7lZPp3ivO7JBEJEQr3IDV7cSb3v7aKnu1bMHXCMBLjmvpdkoiEEIV7kHHO8dh7G3nq43RO6xXP368ZSsvYaL/LEpEQo3APIsVl5dwzZyWvL9/BlanJ/PbSE4nW5F8iUgMK9yCRV1jKxH+lsXBrLnef25vbRvXQ5F8iUmMK9yCQmVvI+KmLyMw9xBNjBzNmcJLfJYlIiFO4+2x55n5umr6Y0nLHv24czoju7fwuSUTCgMLdR++u2cmdrywjoWUTXhk/nJ7t9dQkEakbCnefvPDfrTz85loGdm7N89enEt+iid8liUgYUbg3sPKA47dvrmXq5xmc278Df71yCE1jNPmXiNQthXsDOlRSzp2vLOO9tbu44dRu/OLCvkRp8i8RqQcK9waSc7DiOacrt+/ngYv7MeHUbn6XJCJhTOHeADbtOsgN0xeTc7CYf157Euf07+h3SSIS5o4b7mb2AnARsNs5d6LX1haYBaQAGcAVzrl9VnHXzRPABUAhMN45t7R+Sg9uJWUBPli3i1mLM/l0Uw7tmscwa+LJDEpu7XdpIhIBqnLlPg14CphRqe1e4EPn3CNmdq/3/h7gfKCX9xoBPON9jRgbdh5k1uJMXl+eRW5BCYlxsfz4zJ5cM7IrHVrF+l2eiESI44a7c+5TM0s5rHkMMMpbng78h4pwHwPMcM454Csza21mic657LoqOBgdLCplwYpsZqVlsiJzP9FRxtn9OnBFajKn9UrQh6Yi0uBq2ufeoVJg7wQ6eMtJQGal7bZ7bUeEu5lNBCYCdOkSek8Wcs6xOGMfsxZn8taqbA6VlnNChxb88sK+XDokiXYaty4iPqr1B6rOOWdmrgb7TQGmAKSmplZ7f7/sPlDE3KVZvJqWyZY9BbRo0phLhiRx5bBkBnWO02RfIhIUahruu77pbjGzRGC3154FJFfarrPXFvICAcdv31zH9C8zKA84hqe05bYze3LBgI40i9GgIxEJLjVNpfnA9cAj3td/V2q/3cxeoeKD1Lxw6G8vDzjumbuSOUu2M3ZYMhNP7073BM0DIyLBqypDIV+m4sPTeDPbDjxARajPNrMbgW3AFd7mb1ExDDKdiqGQE+qh5gZVWh7gZ7NXsGDFDn76/RO446ye6noRkaBXldEyVx1j1VlH2dYBk2tbVLAoLivnjpeX8e6aXdx7fh9uOaOH3yWJiFSJOouPoai0nFtfXMLHG3I0XYCIhByF+1EUlpRx84w0vti8l99fOoCrR4TeUE0RiWwK98McLCrlhmmLWbJtH49dPojLhnb2uyQRkWpTuFeSV1jKuKmLWJOVx9+uGsJFAzv5XZKISI0o3D25BSVc+9xC0nfn8/drhmrmRhEJaQp3YPfBIq59biHb9hYyZdxJjOrd3u+SRERqJeLDPTvvENc8u5CdB4qYOn4Yp/SM97skEZFai+hwz8wt5OrnvmJfQSkzbhhOakpbv0sSEakTERvuW/cUcPWzX1FYUs7Mm0boIRoiElYiLtyLy8r5dOMe7n9tFeUBx8s3j6Rfp1Z+lyUiUqciItxLygJ8vnkPb67M5t01OzlYVEZiXCwv3TSCXh1a+l2eiEidC9twLysP8MXmvby5Mpt31uwk71ApLWMbc06/jlw0KJFTe8QT07iR32WKiNSLsAr38oBj4Za9LPCu0HMLSmjRpDFn9+vAhQMSOe2EeJo0jvK7TBGRehfy4V4ecCzOyOXNldm8vTqbPfklNIuJ4qy+HbhoYCJnnJBAbLQCXUQiS0iH+0frd3Hv3FXsPlhMbHQjRvdpz0UDO3Fm7/Y0jVGgi0jkCulw79S6KUO6tObCgZ04q097mjcJ6V9HRKTOhHQa9unYin9el+p3GSIiQUfDRUREwpDCXUQkDCncRUTCkMJdRCQMKdxFRMKQwl1EJAwp3EVEwpDCXUQkDJlzzu8aMLMcYFsNd48H9tRhOeFCx+VIOiZH0jE5Uigdk67OuYSjrQiKcK8NM0tzzuk21cPouBxJx+RIOiZHCpdjom4ZEZEwpHAXEQlD4RDuU/wuIEjpuBxJx+RIOiZHCotjEvJ97iIicqRwuHIXEZHDKNxFRMJQSIe7mZ1nZhvMLN3M7vW7nmBgZhlmtsrMlptZmt/1+MXMXjCz3Wa2ulJbWzN738w2eV/b+FljQzvGMXnQzLK882W5mV3gZ40NzcySzexjM1trZmvM7E6vPeTPlZANdzOLAp4Gzgf6AVeZWT9/qwoaZzrnBofDWN1amAacd1jbvcCHzrlewIfe+0gyjSOPCcDj3vky2Dn3VgPX5Lcy4C7nXD9gJDDZy5GQP1dCNtyB4UC6c26Lc64EeAUY43NNEiScc58CuYc1jwGme8vTgUsatCifHeOYRDTnXLZzbqm3fBBYByQRBudKKId7EpBZ6f12ry3SOeA9M1tiZhP9LibIdHDOZXvLO4EOfhYTRG43s5Vet03IdT/UFTNLAYYACwmDcyWUw12O7nvOuaFUdFdNNrPT/S4oGLmKMcAaBwzPAD2AwUA28Ji/5fjDzFoAc4GfOOcOVF4XqudKKId7FpBc6X1nry2iOeeyvK+7gdeo6L6SCrvMLBHA+7rb53p855zb5Zwrd84FgGeJwPPFzKKpCPaZzrl5XnPInyuhHO6LgV5m1s3MYoCxwHyfa/KVmTU3s5bfLAPnAKu/e6+IMh+43na7B6IAAACzSURBVFu+Hvi3j7UEhW8CzHMpEXa+mJkBzwPrnHN/qbQq5M+VkL5D1Ru29VcgCnjBOfc7n0vylZl1p+JqHaAx8FKkHhMzexkYRcX0rbuAB4DXgdlAFyqmmL7CORcxHzAe45iMoqJLxgEZwKRKfc1hz8y+B3wGrAICXvP9VPS7h/S5EtLhLiIiRxfK3TIiInIMCncRkTCkcBcRCUMKdxGRMKRwFxEJQwp3EZEwpHAXEQlD/wfMyqYk1SVYGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ddf['days'], df['confirmed'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Based on model found at https://github.com/omerka-weizmann/2_day_workweek/blob/master/code.ipynb\n",
    "def SEIR_model(y,t,config):\n",
    "    \"\"\"\n",
    "    SEIR model\n",
    "    @y,t: - variables for the differential equations\n",
    "    @config: include - rates beta,gamma for differential equations\n",
    "    @rfunuc: a function that maps time to viral reproduction rate\n",
    "    \"\"\"\n",
    "    S,E,I,R = y\n",
    "    Tinc,Tinf = config[\"Tinc\"],config[\"Tinf\"]\n",
    "    Rt = config[\"Rt\"]\n",
    "    dydt = [-Rt/Tinf * (I*S),\n",
    "            Rt/Tinf * (I*S) - (1/Tinc)*E,\n",
    "            (1/Tinc)*E - (1/Tinf)*I,\n",
    "            (1/Tinf)*I]\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tinc = 3\n",
    "Tinf = 15\n",
    "startI = 0.00005\n",
    "\n",
    "tmax = ddf['days'].max()+1\n",
    "t = np.linspace(1,tmax,tmax)\n",
    "config = {'Rt' : 3.0, 'Tinc': Tinc, 'Tinf': Tinf, 'beta': 0.25, 'gamma': 0.25}\n",
    "SEIR_y0 = [1-startI,startI/2,startI/2,0]\n",
    "\n",
    "modelOutput = odeint(SEIR_model, SEIR_y0, t, args=(config,), atol=1e-12, rtol=1e-12)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelSums = pandas.DataFrame(modelOutput[:,[1,2,3]]).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe10c3d60a0>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xUVfrH8c+TRuihhBBSCCX0EiA0FUWwwwoqoKCIyoq97Koruv7U3dVV1+7qqlhWadIEQUVQELAiECC0UEJNQkghpBHSZs7vjzusUUFC2p3yvF+vvDJz5k7mYRy/OTn33HPEGINSSinv4md3AUoppWqehrtSSnkhDXellPJCGu5KKeWFNNyVUsoLBdhdAEDLli1NTEyM3WUopZRHSUhIyDbGhJ7qMbcI95iYGDZs2GB3GUop5VFE5ODpHtNhGaWU8kIa7kop5YU03JVSygtpuCullBfScFdKKS+k4a6UUl5Iw10ppbyQW8xzV0opX1JYUs6W1FwSU/LoGdGU82Jb1vhraLgrpVQtKnc42ZVRwOaUXBJTctmcksuezEJObqVxx9AOGu5KKeXOjDGkHjtBYmoumw/lkpiay9a0PIrLnAA0axBIXFQII3q2oXdUU3pHhtCsYVCt1KLhrpRS1VBc5mB+Qiqrd2aSmJpLdmEpAEEBfvRo04QJA9rSO6opfaKaEdW8PiJSJ3VpuCulVBWUljuZtyGFN1Ylk55XTPvQhlzQqRVx0SHERYbQJbwxgf72zVnRcFdKqbNQ7nCycGMar329h9RjJ+jXthkvju3N4A4t6qxXXhka7kopVQkOp2FJYhqvrtjDgaNF9IpsylOje3BBp1C3CvWTzhjuIhIMfAPUcx2/wBjzhIh8AFwA5LkOvckYs1msf+WrwBVAkat9Y20Ur5RStc3pNHy+NZ1XVuxmb9ZxuoY34Z0b47moayu3DPWTKtNzLwGGGWMKRSQQ+E5EvnA99pAxZsGvjr8ciHV9DQTedH1XSimPYYxh+fYMXlmxm51HCugU1og3r+/Lpd1b4+fnvqF+0hnD3RhjgELX3UDXl/mdp4wCpruet1ZEQkQk3BiTXu1qlVKqlhlj+HpnJi99tZvth/Np37Ihr14Xx8hebfD3gFA/qVJj7iLiDyQAHYE3jDE/icgdwNMi8jiwEphqjCkBIoCUCk9PdbWl/+pnTgGmAERHR1f336GUUtW2JTWXxxdvZ3NKLtHNG/DC2N6MjmtDgI2zXqqqUuFujHEAcSISAiwSkR7AI8ARIAiYBjwM/L2yL2yMmeZ6HvHx8b/3l4BSStW6/dnHmfjeOuoH+vPM1T0Z0y/S1qmM1XVWlRtjcoFVwGXGmHRjKQH+CwxwHZYGRFV4WqSrTSml3FJeURmTP1iPn8C82wYzfkC0Rwc7VCLcRSTU1WNHROoDFwM7RSTc1SbAaGCb6ylLgBvFMgjI0/F2pZS7KnM4uXN2AinHinh7YjzRLRrYXVKNqMywTDjwoWvc3Q+YZ4z5TES+FpFQQIDNwO2u45diTYNMxpoKeXPNl62UUtVnjOFvn27n++SjPD+mFwPaNbe7pBpTmdkyW4A+p2gfdprjDXBX9UtTSqna9eEPB5i59hC3XdCesfFRZ36CB/HsQSWllKqi1bsy+ftnO7i4WxgPX9rF7nJqnIa7Usrn7Mko4J7Zm+jcugmvXBvnERclnS0Nd6WUT8k5XsrkDzdQL9CfdyfF07Cedy6x5Z3/KqWUOoWScge3z0jgSH4xc6cMIiKkvt0l1RrtuSulfIIxhscWbWPdgRxeGNubPtHN7C6pVmm4K6V8wrRv9jE/IZV7h8dyZe82dpdT6zTclVJe76sdGTy7bCcjeoVz//BYu8upExruSimvtv1wHvfN2USviKa8OLa3V86MORUNd6WU18osKObWDzfQJDiQd26MJzjQ3+6S6ozOllFKeaXiMgdTpidwrKiM+bcPplWTYLtLqlMa7kopr2OM4S8LtrA5JZe3buhHj4imdpdU5zTclVJeo6i0nC+2HmHO+kOsP3CMv1zWmct6tLa7LFtouCulPJoxhk0puczfkMKniekUlpQT06IB/zeyG7ecG2N3ebbRcFdKeaSsghIWbUpl3oZUkjMLqR/oz4he4YyLj6J/TDOsrSZ8l4a7UspjlDucrNqVxbwNKazamUm509A3OoRnr+7JiF7hNA4OtLtEt6HhrpRye8mZhcxPSGHhxjSyCkpo2SiIyee1Y2x8JB1bNba7PLek4a6UqnV5J8p4ZmkS2YWlZ/3crMISElNy8fcTLuzcinHxkVzYpZXH73Fa2zTclVK1yuE03D9nE9/uyaZT2Nn3soMD/Xjk8i5c1TeCVo19a656dWi4K6Vq1Qtf7mLVriyeGt2DGwa1tbscn6F/1yilas2niYd5c/VeJgyM1mCvY2cMdxEJFpF1IpIoIttF5G+u9nYi8pOIJIvIXBEJcrXXc91Pdj0eU7v/BKWUO9p+OI+HFiTSP6YZT/6hu93l+JzK9NxLgGHGmN5AHHCZiAwCngNeNsZ0BI4Bk13HTwaOudpfdh2nlPIhRwtLmDI9gWYNgvjP9f0ICtBBgrp2xnfcWApddwNdXwYYBixwtX8IjHbdHuW6j+vx4eLrVxMo5UPKHE7unLWR7MIS3p7Yj9DG9ewuySdV6tepiPiLyGYgE/gK2AvkGmPKXYekAhGu2xFACoDr8TygRU0WrZRyX099toOf9ufw7DU96RUZYnc5PqtS4W6McRhj4oBIYADQpbovLCJTRGSDiGzIysqq7o9TSrmBeetT+PDHg9w6pB1X9Ym0uxz3d3gTlBXXyo8+q4EwY0wusAoYDISIyMmplJFAmut2GhAF4Hq8KXD0FD9rmjEm3hgTHxoaWsXylVLuYuOhYzz2yTaGxLbk4cuq3f/zfumJ8N8RsPzRWvnxlZktEyoiIa7b9YGLgSSskB/jOmwSsNh1e4nrPq7HvzbGmJosWinlXjLyi7l9RgKtmwbz7/F9CNCrR39fbgrMGgf1m8H5D9XKS1TmIqZw4EMR8cf6ZTDPGPOZiOwA5ojIU8Am4D3X8e8BM0QkGcgBrquFupVSbqK4zMGUGQkUlpQzY/JAQhoE2V2SeyvOg9njoOwETF4OTcJr5WXOGO7GmC1An1O078Maf/91ezEwtkaqU0q5NWMMj32yjUTXjkedW+siXr+rvBTmToTsPXDDx9Cqa629lC4/oJSqsg9+OMCChFTuGx7rszseVZox8Ol9sH8NjH4L2l9Qqy+nA2NKqSr5ITmbpz5P4pJuYdw3PNbuctzfmucgcTYMfRTixtf6y2m4K6XOWkpOEXfO3kj7lg156do4/Pz0OsXftXk2rH4G4q6HC/5SJy+p4a6UOitFpeXcOn0DTqfhnRvjaVRPR3d/177VsOQeaD8U/vAq1NEF+xruSqlKS887wT2zN7E7o4DXJ/QlpmVDu0tybxnbrROoLTvDuOngX3fbAOqvXKXUGSVnFvL2mr18sjkNp4HHR3bj/E568eHvyk+35rIHNYTr50Fw0zp9eQ13pdRpbTp0jLfW7OXLHRnUC/Dj+oFtmXxeO6KaN7C7NPdWUgCzx0JxLtz8BTSt+6UYNNyVUr9gjOGbPdm8uTqZtftyaFo/kHsu7Mikc2Jo0UhXeDwjRznMvwkydsCEeRDey5YyNNyVUgCUO5x8se0Ib67ey470fFo3CeaxEV0ZPyCahnrStHKMgaUPQPIK6+Rp7EW2laL/xZTyccVlDj7emMq0b/Zx8GgR7UMb8q8xvRgdF6GbbJyt716GhA/gvD9Dv5tsLUXDXSkflV9cxqy1h3jvu/1kF5bQO7Ipj9zQj0u6hem89arYugBW/g16jIFh/2d3NRruSvmawpJy/vvdft75dh/5xeUMiW3JHUPjGNy+BbppWhUd+B4+uQPanguj/wN+9v/Fo+GulI8oKi1n+o8HeXvNXo4VlXFRV2vZgJ6RdTtFz+ukJ8JH4yGkLVw7EwLc46SzhrtSXq64zMHMtQd5a81esgtLuaBTKH++uBO9o3QLvGrL3AkzroLgJjBxETRobndF/6PhrpSXKil3MHd9Cm+sSiYjv4RzOrTg7Ymd6NfWfQLIox3dC9NHgV8A3LgYQqLsrugXNNyV8jJlDicLElJ5/etk0nJP0D+mGa9c24fBHXSf+hqTm2IFu6MUbl4KLTrYXdFvaLgr5SXKHU4+2XyY11bu4VBOEXFRITx7TU/O69hST5TWpIIMK9iL82HSklrdcKM6NNyVciNFpeU4z3LHYWMMX+/M5NUVe9iXfZweEU14/6Z4LuzcSkO9phXlwIzRUHDEGmNvE2d3Rael4a6Um3h1xR5eXrG7ys/v0roxb93Qj0u7h2mo14biPOvk6dG9cP18iB5od0W/S8NdKTewamcmL6/YzcXdwhgQc/YnPKNbNODirnrxUa0pPW6t8JixDa6bXetb5NUEDXelbJaWe4I/zdtM1/Am/Ht8H4ID/e0uSVVUVgxzJkDqOhjzPnS61O6KKuWMl1GJSJSIrBKRHSKyXUTuc7U/KSJpIrLZ9XVFhec8IiLJIrJLRDzjnVDKBqXlTu6evZFyh+E/1/fVYHc3jjJrhcd9q2HUG9D9KrsrqrTK9NzLgQeMMRtFpDGQICJfuR572RjzQsWDRaQbcB3QHWgDrBCRTsYYR00WrpQ3+NeynWw6lMsbE/rSTnc1ci9OByycAru/gBEvQtwEuys6K2fsuRtj0o0xG123C4AkIOJ3njIKmGOMKTHG7AeSgQE1UaxS3mTZtiO8+91+bjonhhG9wu0uR1XkdMKSe2H7Qrj4H9D/j3ZXdNbOanUbEYkB+gA/uZruFpEtIvK+iDRztUUAKRWelsopfhmIyBQR2SAiG7Kyss66cKU82aGjRTy0INFaifGKLnaXoyoyBpZNhc0z4YKpcO69dldUJZUOdxFpBHwM3G+MyQfeBDoAcUA68OLZvLAxZpoxJt4YEx8aqnsxKt9RXObgztkJCPD6hL7UC9BxdrdhjLVs77q3YfDdMHSq3RVVWaXCXUQCsYJ9ljFmIYAxJsMY4zDGOIF3+HnoJQ2ouMhCpKtNKQU8/XkS29LyeXFcnO5F6k6MgTXPWRtuxN8ClzwFHny9QGVmywjwHpBkjHmpQnvFQcKrgG2u20uA60Sknoi0A2KBdTVXslKea0niYWasPciU89tzcbcwu8tRJxkDK56E1c9A7wlwxYseHexQudky5wITga0istnV9igwXkTiAAMcAG4DMMZsF5F5wA6smTZ36UwZpWBvViGPfLyFfm2b8dClne0uR53kdFpj7OvetnrsV7zoFpttVNcZw90Y8x1wql9hS3/nOU8DT1ejLqW8yolSB3fN2khQgB+vT+hDoL/nh4dXcDrg03th00xrjN3Dh2Iq0itUlaoDTy7Zzq6MAv57U3/Cm9a3uxwF1gVKi26HbQvggodh6CNeE+yg4a5UrVuQkMrcDSncM6wjQzu3srscBVBeAvNvhl2fw0VPwnl/sruiGqfhrlQt2nWkgMc+2crg9i24/6JOdpejAEqLYO4NsHclXP4vGHib3RXVCg13pWrJ8ZJy7pyVQKN6gbw6Pg5/XbHRfiUFMPs6OPg9XPk69J1od0W1RsNdqVpgjOHRRVvZn32cmX8cSKvGwXaXpE4cg5lj4PAmuOZd6DnG7opqlYa7UrXgo3UpLN58mAcv6cQ5HVraXY46nm3toJS1C8ZNh64j7a6o1mm4K3Ua29LyeG7ZTvKLy8/6uUmH8zm/Uyh3Du1YC5Wps5Kfbu15mnsQxn8EHS+yu6I6oeGu1K84nYb3v9/Pc8t2EtIgiG7hTc76Z4zoFc5jI7rqzkh2yz0EH14JhZlw/QJoN8TuiuqMhrtSFWQVlPDA/ES+2Z3Fxd3C+Nc1vWjWMMjuslRVHN1rBXtJAdy4GKL6211RndJwV8plze4sHpi3mYLicv4xugc3DIzWjaY9VcYOa4zdWQ43fQrhve2uqM5puCufV1Lu4Pllu3j3u/10DmvMrD8OonPrxnaXpapq3xprHntgA7jpc2jV1e6KbKHhrnzavqxC7p2ziW1p+dw4uC2PXtFV9zH1ZIlzYPHd0KIDXD8fQqLtrsg2Gu7KJxljWJCQyhNLthMU4Me0if24pHtru8tSVWUMfPMCrHoKYobAtTOhfojdVdlKw135nPziMv66aBufJh5mcPsWvHxtHK2b6kVGHstRBp//GTZOh57jYNTrEFDP7qpsp+GufErCwWPcN2cT6XnFPHRpZ26/oIMuC+DJSgpg3iRrnZghD8Kwx7xqZcfq0HBXPsHhNLy5OpmXV+yhTUgw828fTN/oZmd+onJf+ekwe6w1M+YPr0K/m+yuyK1ouCuvl19cxu0zEvhh71Gu7N2Gp67qQZPgQLvLUtWRsQNmjYXiXJgwF2Ivtrsit6PhrrxadmEJk95fx+6MAv41phdj+0Xq3HVP97+pjvXh5qU+OYe9MjTcldc6nHuCG979icN5J3h3Un8u6BRqd0mqunSqY6VpuCuvtC+rkInvrSP/RBkzJg+kf0xzu0tS1aFTHc/aGXfpFZEoEVklIjtEZLuI3Odqby4iX4nIHtf3Zq52EZHXRCRZRLaISN/a/kcoVdGOw/mMe/tHisscfDRlkAa7p3OUWZtYr3oKel0LNyzUYK+EymzBXg48YIzpBgwC7hKRbsBUYKUxJhZY6boPcDkQ6/qaArxZ41UrdRoJB3O4btqPBPr7Me/2wfSIaGp3Sao6ivPho+usOexDHoSr3oYAXcitMs44LGOMSQfSXbcLRCQJiABGAUNdh30IrAYedrVPN8YYYK2IhIhIuOvnKFVrvt2TxZTpCbRuGsyMyQOIbNbA7pJUdWTvgTkTrNUddarjWTurMXcRiQH6AD8BYRUC+wgQ5rodAaRUeFqqq+0X4S4iU7B69kRH60kRVT3LtqVz70eb6dCqEdNvGUBoY71C0aPtXAqLbgP/QGu5Xh9ah72mVGZYBgARaQR8DNxvjMmv+Jirl27O5oWNMdOMMfHGmPjQUJ3FoKpu/oYU7py1kR4RTZhz6yANdk/mdMLqZ2HOeGjeDqas0WCvokr13EUkECvYZxljFrqaM04Ot4hIOJDpak8Doio8PdLVplSNe/+7/fz9sx0MiW3J2xP70SBIJ4B5rOJ8q7e+ayn0Hg8jX7bmsqsqqcxsGQHeA5KMMS9VeGgJMMl1exKwuEL7ja5ZM4OAPB1vVzXNGMMrK3bz9892cFn31rw7KV6D3ZNl7YZ3hsHu5XDZczD6TQ32aqrM/w3nAhOBrSKy2dX2KPAsME9EJgMHgXGux5YCVwDJQBFwc41WrHye02l46vMk3v9+P2P6RfLs1T0J8K/0CKNyNzuXwsIp1kqOk5ZAzHl2V+QVKjNb5jvgdNdrDz/F8Qa4q5p1KXVK5Q4nUxduZUFCKjedE8PjI7vpJtSeyumENc/BmmehTR/rwqSmkXZX5TX071jllsodTtLzikk9doKUY0WkHjtBak4RO9Lz2XmkgPuGx3L/RbG6ToynKs6DhbfB7i+g9wQY+ZIOw9QwDXdlC2MMmQUlpOQUWeGdY4V4Ss4JUnOLSM8tptz58wQsEQhvEkxk8wY8c3VPxg/Q6bMeK2u3NX/92H64/HkYcKuuwV4LNNxVncvIL+aejzaxbn/OL9pDG9cjsll9+kQ148re9Yls1oCoZg2Ial6f8Kb1CQrQcXWPt/Nzq8ceGAw3LoGYc+2uyGtpuKs69UNyNvfO2URRqYOHL+tCl9aNiWpen4iQBtQP0o2pvZbT4Rpffw7a9IVrZ+j4ei3TcFd1wuk0vLlmLy9+uYv2oY2YM6UvHVs1trssVRfy02HhrXDgW4i7Hka8ZPXcVa3ScFe1LreolD/N3cyqXVmMimvDP6/qScN6+tHzCbuWwSd3QHkxjHrDCncdX68T+n+YqlWJKbncOWsjWQUl/GN0D24YGK0zXHxBeQl89QT89CaE9YQx70NoJ7ur8ika7qpWGGOYufYg//gsidDG9Zh/+2B6R+ka3D4hOxkW3AxHtsDA2+Giv+kwjA003FWNO15SzqOLtrJ482GGdWnFS+N6E9JA1+D2esZA4kfw+YPW1abj50Dny+2uymdpuKsalZxZwO0zN7Ivq5CHLu3MHRd00CtIfUFxPnz+AGydB23Pg2vegSZt7K7Kp2m4qxqzeHMajyzcSoMgf2ZOHsg5HVvaXZKqC2kJsGAy5B6EC/8KQx4AP53WajcNd1VtJeUOnvosiRlrD9I/phmvT+hLWBMdY/V6Tif8+Dqs/Bs0ag03LYW2g+2uSrlouKtqOZB9nPvmbCIxNY/bzm/Pg5d2JlBXaPR+hZmw6HbYuxK6jIQr/w0NdCNyd6LhrqrkRKmDN1cn89Y3+6jn78fbE/txaffWdpel6sLer60lBEryrQuS4m/RuetuSMNdnRVjDMu3Z/CPz3aQlnuC0XFteOSKrjoM4wtKCuGrx2HDexDaxdrbNKyb3VWp09BwV5W2L6uQJ5Zs59s92XRp3Zi5UwYxsH0Lu8tSdWHfGlhyN+SmwOC7rROnQQ3srkr9Dg13dUbHS8p5fVUy7367j+AAf574QzcmDmqrux/5goq99eYd4JZlED3I7qpUJWi4q9MyxvD51nSe/jyJ9LxixvSL5OHLuhDauJ7dpam6oL11j6bhrk5pT0YBTyzZzg97j9K9TRNen9CHfm11NoRP0N66V9BwV79QUFzGayv38N/vD9AgyJ9/jOrOhIFt8derTH2D9ta9hoa7AqwhmMWbD/PPpUlkFZZwXf8oHrykMy0a6RCMT9Deutc5Y7iLyPvASCDTGNPD1fYkcCuQ5TrsUWPMUtdjjwCTAQdwrzFmeS3UrWpQUno+TyzezroDOfSObMq0G+OJ0xUcfYf21r1SZXruHwCvA9N/1f6yMeaFig0i0g24DugOtAFWiEgnY4yjBmpVNSzvRBkvf7WbGWsP0iQ4gGeu7sm18VG60Jev0N66VztjuBtjvhGRmEr+vFHAHGNMCbBfRJKBAcCPVa5Q1Tin0/DxxlSeW7aTo8dLuX5gNA9e0lmX5fUVxlgbVS+bCnmp2lv3UtUZc79bRG4ENgAPGGOOARHA2grHpLrafkNEpgBTAKKjo6tRhjob29LyeHzxNjYeyqVvdAgf3DyAHhFN7S5L1ZWc/fDFw7BnObTqDrcsh+iBdlelakFVw/1N4B+AcX1/EbjlbH6AMWYaMA0gPj7eVLEOVUm5RaW88OUuZv10iOYNgnh+TC+u6RupQzC+orwEvn8Nvn0B/ALgkqdh4G3gH2h3ZaqWVCncjTEZJ2+LyDvAZ667aUBUhUMjXW3KJk6nYe6GFP61bCd5J8qYNDiGP13ciab19X9qn7H3a2t3pJy90G00XPpPaHrKP6iVF6lSuItIuDEm3XX3KmCb6/YSYLaIvIR1QjUWWFftKlWVbE7J5YnF20hMzWNATHP+Nqo7XcOb2F2Wqiv5h2H5X2H7QmjeHm5YCB2H212VqiOVmQr5ETAUaCkiqcATwFARicMaljkA3AZgjNkuIvOAHUA5cJfOlKl7RwtLeH75LuZuSKFlo3q8cm0co+LaILosq29wlMO6t2HVP8FRZp0sPede3aTax4gx9g93x8fHmw0bNthdhldYsSODB+YncryknJvPjeHe4bE0DtYhGJ9xaK21l2nGNoi9BC5/zuq1K68kIgnGmPhTPaZXqHqJcoeTF77czVtr9tIjogkvjYujU1hju8tSdeX4UVjxOGyaCU0i4dqZ1g5J+teaz9Jw9wKZBcXcM3sTP+3PYcLAaB4f2Y3gQN2g2Cc4ymHjB/D1U1BSAOfeDxf8BYIa2l2ZspmGu4dbu+8o93y0iYLiMl4a15ur+0baXZKqC8bA7uXw1f9B9m6IGQJXvACtuthdmXITGu4eyuk0vP3NPp5fvpOYFg2ZOXkgnVvrMIxPSE+ELx+D/d9Ai45w3WzofIUOwahf0HD3QHlFZTwwP5EVSRmM6BXOc9f0olE9/U/p9fLSrOGXxI+gfjO4/HmIv1kvRFKnpIngYbal5XHHrASO5BXz5B+6MemcGJ3i6O1KCuC7V+DHN8A44dx7YcgDEKzLRqjT03D3EMYYPlqXwpOfbqdlwyDm3jaYvtHN7C5L1SZHOWyabs1XP54FPcbA8MehWVu7K1MeQMPdAxSVlvPYom0s3JTGkNiWvHpdH5o31BUcvZYxsOcr62Rp1k6IHgzj50JkP7srUx5Ew93N7c0q5M6ZG9mdWcD9F8Vyz7BY3fLOmx3Zap0s3bfauvhI56urKtJwd2Ofb0nnLwsSCQrw48ObB3B+p1C7S1K15eheWPMcbJkH9UPgsucg/hYI0L/QVNVouLuh0nIn/1yaxAc/HKBPdAhvTOhLm5D6dpelakPOfvjmeWsGTECwdbL0vD9bAa9UNWi4u5nDuSe4a/ZGNh3K5eZzY3jk8q4EBfjZXZaqaccOWmurb5plTWUcdCecex80amV3ZcpLaLi7kTW7s7h/ziZKy528MaEvI3qF212Sqml5qfDti7BxhjWOPuBWOO9P0Li13ZUpL6Ph7gYcTsNrK/fw2td76NSqMf+5oS8dQhvZXZaqSfnp8N1LkPCBNRum3yRr+EU3zVC1RMPdZkcLS7h/7ma+3ZPN1X0jeHp0T+oH6aJfXqMgA75/Bda/B8YBfW6wLkAK0X2DVe3ScLdRwsFj3D17I0ePl/LM1T25rn+UXm3qLY5nW6G+7l1wlELceDj/IWgWY3dlykdouNvAGMN/vz/AP5cmER4SzMI7zqFHhF5K7hXyUuHH/1jDL+UnoNe1Vqi36GB3ZcrHaLjXsYLiMh7+eAtLtx7hoq5hvDi2N00b6MJPHi8zCb5/FbbOt8bUe46xQr1lrN2VKR+l4V6Hdh7J586ZGzmYU8TUy7tw2/ntdRjG0x380Rp+2b0MAhtA/1th8J06pq5sp+FeRz5OSOWvn2ylcXAgs/44kEHtW9hdkqoqpxN2f2H11FN+ggYtYOij1rTGBs3trk4poBLhLiLvAyOBTGNMD1dbc2AuEAMcAMYZY46J1Q19FbgCKAJuMsZsrJ3S3V+Zw8l3ydnM35DC0q1HGNS+OZSZoboAAA1USURBVK+N70OrxroLvUcqL4Wt8+D71yB7l9U7v+IFiLseghrYXZ1Sv1CZnvsHwOvA9AptU4GVxphnRWSq6/7DwOVArOtrIPCm67vPKHc4Wbsvh8+2HGbZ9iPkFpXRODiAe4d15N7hsQT469WmHqc43zpBuvY/UJAOYT3hmveg22jw1z9+lXs64yfTGPONiMT8qnkUMNR1+0NgNVa4jwKmG2MMsFZEQkQk3BiTXlMFuyOn07D+QA6fbUnni23pZBeW0jDIn4u7hTGyVxuGdGpJvQCdu+5x8tJg/bvWHPWSPGh3Pox6AzoM01UaldurarcjrEJgHwHCXLcjgJQKx6W62n4T7iIyBZgCEB3teSefjDFsSsnl08TDLN2aTkZ+CcGBfgzvEsbIXuFc2KUVwYEa6B7HGDj4PaybBkmfWTsfdbvSWvclQtdTV56j2n9TGmOMiJgqPG8aMA0gPj7+rJ9vB2MM29Ly+WzLYT7bkk5a7gmC/P0Y2jmUkb3bMLxLKxrqXqaeqfS4NY3xp2mQuR2CQ2DwXdD/j7rzkfJIVU2ijJPDLSISDmS62tOAqArHRbraPF7eiTLunr2Rb/dkE+AnDIltyZ8v7sTF3cNoEqzz1D1Wzn5r6GXTDCjOs8bTr/y3taWdniRVHqyq4b4EmAQ86/q+uEL73SIyB+tEap43jLcfOlrEzR+s41BOEY+N6Mo1fSNpptvceS6nE/Z9Devegd3Lwc8ful4JA6ZA9CAdT1deoTJTIT/COnnaUkRSgSewQn2eiEwGDgLjXIcvxZoGmYw1FfLmWqi5TiUczOHW6Qk4nIYZk3V+ukcrzofNs2H9O3A0GRqGWleRxt8MTdrYXZ1SNaoys2XGn+ah4ac41gB3Vbcod7F4cxoPLdhCm6bBvH9Tf9rrMrye6fBm2DgdtsyF0kKI7A9XvwPdRkFAPburU6pW6Nm/UzDG8NrKZF5esZsB7Zrz9g39dBjG05w4Blvmw6bp1qbT/vWgx9XW0EtEX7urU6rWabj/Skm5g6kfb2XRpjSu7hvBM1f31DnqnsLphAPfWr30pE/BUQKte1lXkfYcA/Wb2V2hUnVGw72CnOOl3DZjA+sPHOPBSzpx14UddWEvT5CXZo2lb5oBuQchuCn0vRH6ToTw3nZXp5QtNNxdkjMLmfzhetLzinl9Qh9G9tITbG6tvNRavGvjDNi70rrYqN35MOz/oOtICKxvd4VK2UrDHfghOZvbZyYQ6O/HnCmD6Butf767rYwdsHkWJM6Bomxo3Mbati7uemjezu7qlHIbPh/u89an8OiirbRr2ZD3b+pPVHO9cMXt5KbAtgWwdQFkbAO/QOh8uTX00mGYNU9dKfULPhvuTqfhX8t38daavQyJbckb1/fVK03dSVEObF9kBfqhH6y2yP5w+fPQ/SpoFGpvfUq5OZ8M97wTZUz9eAtfbDvChIHR/O3K7gTqUrz2Kz0Ou76w1nhJXgHOcmjZGYY9Zi0HoMMuSlWaz4T7/uzjrEzKYGVSJusP5OAwhsdGdGXyee10RoydHGWwd5UV6Ds/h7Lj0CQCBt0JPcdC6566HIBSVeC14V7mcLLhwDG+3mkF+r7s4wB0CmvEH4e0Z0TPcHpGNrW5Sh/ldMChtbB9oTX0UnTUWoWx11gr0KPPAT/9S0qp6vCqcM8tKmX1rixW7sxkza5M8ovLCfL3Y2D75kw6J4ZhXVrpCVO7OMpg/zeQtMTqoR/PgoD61onRnmOh43BdCkCpGuTR4W6MYW9WISuSMvk6KZMNB3NwGmjZKIhLu7dmeNcwzottSSNdY90eZSesIZekJbBrqbWkbmBD6HQJdP0DxF4C9RrbXaVSXsmjU29BQioPLdgCQLfwJtx1YUeGdw2jV0RT/Px0nNYWJQWw50vr8v/dX1pj6MFNofMV1rK6HS7UC4yUqgMeHe7ndwrlqdE9GNalFW1CNDBsU5RjzXJJ+hT2fm2t6dIwFHqNs3ro7c4Hf51mqlRd8uhwD2sSzA2DdAu0OmcMZO+BPcutzS4O/WhNW2wSCfG3WHuORg3Ui4uUspFHh7uqQ2XFcPA7a6hlz3I4dsBqb9UNzrnH6qG36avTFpVyExru6vTyD1s98z1fwr7VUFYEAcHQ7gIr0GMvhZCoM/4YpVTd03BXP3M6IC3BFejLrU0uAJpGQdwEK8zbDdETokp5AA13X5d7yOqVn/wqOgrib42ZX/SkFeituupwi1IeRsPd15w4Bvu//TnMc/Za7Y1aQ8eLrLnnHYfrrkVKeTgNd29XXgIp62DfKivMD2+yNrYIbAgx58GAW6H9UAjtor1zpbxItcJdRA4ABYADKDfGxItIc2AuEAMcAMYZY45Vr0xVaU6nteb5/jXW1aEHf4DyE9ZQS2Q8nP8XK8wj+kGAbvqtlLeqiZ77hcaY7Ar3pwIrjTHPishU1/2Ha+B11Kk4yiA9EQ5+bwX5oR+ty/zBWi63743WVaFtz4XgJvbWqpSqM7UxLDMKGOq6/SGwGg33mlN2wprRcvAHK9BT1llTFAFadIRuo6wgb3c+NNF9YJXyVdUNdwN8KSIGeNsYMw0IM8akux4/AoSd6okiMgWYAhAdHV3NMrxYSQGk/OQK8x+sYHeUAgJhPaDPRGh7DkQPhsanfKuVUj6ouuF+njEmTURaAV+JyM6KDxpjjCv4f8P1i2AaQHx8/CmP8TnGQM4+SN0AqeutryNbrBOg4g9t+sDA262eefRAndGilDqtaoW7MSbN9T1TRBYBA4AMEQk3xqSLSDiQWQN1eqfiPKsn/r8w3wAncqzHghpZYT7kQatnHtkf6jWyt16llMeocriLSEPAzxhT4Lp9CfB3YAkwCXjW9X1xTRTq8ZwOyEz6OcRT10P2bqyRLbGmInYZYc1oiexv3deFt5RSVVSdnnsYsMi1/2gAMNsYs0xE1gPzRGQycBAYV/0yPYyjHI7usWaxpCdC+hZrfnmZtdUfDVpYAd5zrBXmEX2tNc+VUqqGVDncjTH7gN6naD8KDK9OUR6lvAQyd1gBfjLMM7Zbc8vB2kqudQ9rbZaoAVaYN2unFwwppWqVXqF6NkoKreBOT4QjriDPTLLWMgeo1wRa97LWNA/vDeG9oEUs+OvbrJSqW5o6p1JaBNm7IHMnZCVZAZ65E/IO/XxMgxZWgJ9zkSvIe0NIDPj52Va2Ukqd5NvhXlZsjY1n7rSGVrJ2WkF+7ADWiU7AP8jqfUcNgH43QqvuVpA3aaNDK0opt+X94e50QF6KNX88Zx/k7Ieje61Qz9lnzSEH8AuA5h2s4O59nTVbpVU3aN5eh1WUUh7HO1LLUQ65B63g/l+I77OWsz12EJxlPx8bEGwFdmgX6H6VtVZ5aFfr0n1dSEsp5SU8O9x3fwnLHrY2nDh5UhOs5Wybt7d63l1GWrebt4cWHax1y3VcXCnl5Tw73Bu2sGandL/q5wBv3h4ahel4uFLKp3l2uEf0g3Ef2l2FUkq5HR2fUEopL6ThrpRSXkjDXSmlvJCGu1JKeSENd6WU8kIa7kop5YU03JVSygtpuCullBcSY+zfm1pEsrB2baqKlkB2DZbjLfR9+S19T35L35Pf8qT3pK0xJvRUD7hFuFeHiGwwxsTbXYe70fflt/Q9+S19T37LW94THZZRSikvpOGulFJeyBvCfZrdBbgpfV9+S9+T39L35Le84j3x+DF3pZRSv+UNPXellFK/ouGulFJeyKPDXUQuE5FdIpIsIlPtrscdiMgBEdkqIptFZIPd9dhFRN4XkUwR2VahrbmIfCUie1zfm9lZY107zXvypIikuT4vm0XkCjtrrGsiEiUiq0Rkh4hsF5H7XO0e/1nx2HAXEX/gDeByoBswXkS62VuV27jQGBPnDXN1q+ED4LJftU0FVhpjYoGVrvu+5AN++54AvOz6vMQZY5bWcU12KwceMMZ0AwYBd7lyxOM/Kx4b7sAAINkYs88YUwrMAUbZXJNyE8aYb4CcXzWPAk7uy/ghMLpOi7LZad4Tn2aMSTfGbHTdLgCSgAi84LPiyeEeAaRUuJ/qavN1BvhSRBJEZIrdxbiZMGNMuuv2ESDMzmLcyN0issU1bONxww81RURigD7AT3jBZ8WTw12d2nnGmL5Yw1V3icj5dhfkjow1B1jnAcObQAcgDkgHXrS3HHuISCPgY+B+Y0x+xcc89bPiyeGeBkRVuB/pavNpxpg01/dMYBHW8JWyZIhIOIDre6bN9djOGJNhjHEYY5zAO/jg50VEArGCfZYxZqGr2eM/K54c7uuBWBFpJyJBwHXAEptrspWINBSRxidvA5cA237/WT5lCTDJdXsSsNjGWtzCyQBzuQof+7yIiADvAUnGmJcqPOTxnxWPvkLVNW3rFcAfeN8Y87TNJdlKRNpj9dYBAoDZvvqeiMhHwFCs5VszgCeAT4B5QDTWEtPjjDE+c4LxNO/JUKwhGQMcAG6rMNbs9UTkPOBbYCvgdDU/ijXu7tGfFY8Od6WUUqfmycMySimlTkPDXSmlvJCGu1JKeSENd6WU8kIa7kop5YU03JVSygtpuCullBf6fyRgOvBT9BWHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ddf['days'], ddf['confirmed'])\n",
    "plt.plot(ddf['days'], modelSums[ddf['days']] * population )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0      39.492450\n",
       "1      43.174918\n",
       "2      59.201916\n",
       "3      62.643499\n",
       "4      75.525109\n",
       "5      87.843219\n",
       "6      95.574759\n",
       "7     102.682710\n",
       "8     113.119345\n",
       "9     107.828010\n",
       "10    120.743982\n",
       "11    127.794759\n",
       "12    116.899964\n",
       "13    135.970989\n",
       "14    122.910462\n",
       "15    120.611563\n",
       "16    115.957224\n",
       "17    127.819216\n",
       "18    109.057130\n",
       "19     97.517248\n",
       "20    103.031304\n",
       "21     90.415121\n",
       "22     68.467114\n",
       "dtype: float64"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ddf['confirmed'].values - (modelSums[ddf['days']] * population)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_delta(df, R=3.0, Tinc=3, Tinf=15, startI=0.00005, beta=0.25, gamma=0.25, Toffset=0):\n",
    "    tmax = df['days'].max()+1+Toffset\n",
    "    t = np.linspace(1,tmax,tmax)\n",
    "    config = {'Rt' : R, 'Tinc': Tinc, 'Tinf': Tinf, 'beta': 0.25, 'gamma': 0.25}\n",
    "    SEIR_y0 = [1-startI,startI/2,startI/2,0]\n",
    "    \n",
    "    print(config)\n",
    "\n",
    "    modelOutput = odeint(SEIR_model, SEIR_y0, t, args=(config,), atol=1e-12, rtol=1e-12)\n",
    "    # exposed + infected + recovered \n",
    "    modelSums = pandas.DataFrame(modelOutput[:,[1,2,3]]).sum(axis=1)\n",
    "    # compare to confirmed numbers\n",
    "    delta = np.sum(np.power(df['confirmed'].values - (modelSums[df['days']+Toffset] * population),2))\n",
    "    print(delta)\n",
    "    return delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Rt': 2.6, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "112200.77776798235\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "112200.77776798235"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calc_delta(ddf, R=2.6, Toffset=5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Rt': 3.0, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "55755.57898798663\n",
      "{'Rt': 3.00000001, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "55755.58004834627\n",
      "{'Rt': 1.0, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "780066.0980877717\n",
      "{'Rt': 1.00000001, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "780066.0948056707\n",
      "{'Rt': 2.8856461263772695, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "53461.390149150895\n",
      "{'Rt': 2.8856461363772694, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "53461.38955868635\n",
      "{'Rt': 2.9265480638899053, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52141.052988511125\n",
      "{'Rt': 2.9265480738899052, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52141.05294161367\n",
      "{'Rt': 2.930076969095161, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.18470333635\n",
      "{'Rt': 2.930076979095161, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.18470569595\n",
      "{'Rt': 2.92990792110805, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481163363\n",
      "{'Rt': 2.92990793110805, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481162914\n",
      "{'Rt': 2.9299082421201454, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481147264\n",
      "{'Rt': 2.9299082521201454, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.164811472016\n",
      "{'Rt': 2.929908294110805, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.1648114621\n",
      "{'Rt': 2.9299083041108047, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481145587\n",
      "{'Rt': 2.929908502073442, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481144943\n",
      "{'Rt': 2.929908512073442, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481144746\n",
      "{'Rt': 2.9299093339239914, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481200527\n",
      "{'Rt': 2.9299093439239914, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.164812016505\n",
      "{'Rt': 2.9299085689634192, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481145731\n",
      "{'Rt': 2.929908578963419, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.164811456736\n",
      "{'Rt': 2.929908511468016, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.164811448594\n",
      "{'Rt': 2.929908521468016, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481144978\n",
      "{'Rt': 2.929908511468016, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.164811448594\n",
      "{'Rt': 2.929908521468016, 'Tinc': 3, 'Tinf': 15, 'beta': 0.25, 'gamma': 0.25}\n",
      "52133.16481144978\n"
     ]
    }
   ],
   "source": [
    "out = minimize(lambda x: calc_delta(ddf, R=x[0], Toffset=5), (3), bounds=((1,5),), method=\"L-BFGS-B\", options={\"ftol\":1e-12})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe10c3c7610>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD6CAYAAABamQdMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xUVfrH8c9DChBaCIQYUghVmtRIsyO6SllUFKWjKFjY1RVdy7qrv3V1Xbu7lhUrIB1RERAUARsIBgi9hRJICEkgEGrazPn9ca9rRCCTMrkzk+f9es0rM+fem3mGYb5czpx7jhhjUEopFViqOV2AUkqpiqfhrpRSAUjDXSmlApCGu1JKBSANd6WUCkAa7kopFYA8DncRCRKRdSIy3378oYjsEZFk+9bJbhcR+beIpIjIBhHp4q3ilVJKnV1wKfa9H9gK1C3W9rAxZs4Z+10PtLRv3YG37J/n1LBhQ5OQkFCKUpRSSq1Zs+aQMSbybNs8CncRiQX6Ac8AD5aw+0BgsrGujvpRRMJFJNoYk3GuAxISEkhKSvKkFKWUUjYRST3XNk+7ZV4F/gy4z2h/xu56eUVEqtttMcD+Yvuk2W1nFjVWRJJEJCk7O9vDMpRSSnmixHAXkf5AljFmzRmbHgNaAxcDEcAjpXliY8xEY0yiMSYxMvKs/6tQSilVRp6cuV8C/F5E9gIzgN4i8pExJsNY8oEPgG72/ulAXLHjY+02pZRSlaTEcDfGPGaMiTXGJAC3AUuNMcNFJBqs0THADcAm+5B5wEh71EwPIPd8/e1KKaUqXmlGy5xpqohEAgIkA3fb7QuBvkAKcAq4vVwVKqWUKrVShbsxZjmw3L7f+xz7GOC+8hamlFKq7PQKVaWUCkAa7kop5ZTl/4L0MwciVozy9LkrpZQqq13LYPmzYFwQ07XCf72euSulVGUrzIMFEyCiGVxa0kX/ZaNn7kopVdl+eBVydsGITyCkhleeQs/clVKqMh3eBd+9BO1vhuZnHXRYITTclVKqshgDCx6E4Jrwu2e9+lTaLaOUUpVl4xzYvRz6vgh1orz6VHrmrpRSleH0EVj8mDUyJvEOrz+dnrkrpVRl+PppOHUYhn8M1YK8/nR65q6UUt6WlgRJ70P3uyG6Y6U8pYa7Ukp5k6sI5j8AdaLhqscr7Wm1W0Yppbxp9dtwcCMMngzV61Ta0+qZu1JKeUtuGix9BlpeC21+X6lPreGulFLesuhRMG7o+wKIVOpTa7grpZQ3bF8EWz+HK/4M9RMq/ek9DncRCRKRdSIy337cVERWiUiKiMwUkVC7vbr9OMXenuCd0pVSykcVnISFD0Nka+g53pESSnPmfj+wtdjjfwGvGGNaAEeAMXb7GOCI3f6KvZ9SSlUd3zwPufug38sQHOpICR6Fu4jEAv2Ad+3HAvQG5ti7TMJaJBtgoP0Ye/vV9v5KKRX4MrfAyteh03BIuMSxMjw9c38V+DPgth83AI4aY4rsx2lAjH0/BtgPYG/Ptff/FREZKyJJIpKUnZ1dxvKVUsqHuN3WxGDV68I1f3e0lBLDXUT6A1nGmApdC8oYM9EYk2iMSYyMjKzIX62UUs5Y8z7sW2kFe63fnNNWKk8uYroE+L2I9AVqAHWB14BwEQm2z85jgXR7/3QgDkgTkWCgHnC4witXSilfkrMHvvwbNLsKOg93upqSz9yNMY8ZY2KNMQnAbcBSY8wwYBlws73bKOAz+/48+zH29qXGGFOhVSullC9xu+Gz8daEYANfr/Qx7WdTnnHujwAPikgKVp/6e3b7e0ADu/1B4NHylaiUUj5u9URI/d5agKNerNPVAKWcW8YYsxxYbt/fDXQ7yz55wC0VUJtSSvm+Qymw5ClrigEf6I75mV6hqpRSZeV2wWf3WmPZB/zbJ7pjfqazQiqlVFmtfAP2r4IbJ0LdaKer+RU9c1dKqbLI3g5L/wEX9oMOg52u5jc03JVSqrRcRfDJ3RBaCwa86lPdMT/TbhmllCqtFa/BgbVw8/tQu5HT1ZyVnrkrpVRpZG6GZf+EtjdA+0FOV3NOGu5KKeUpV6HVHVOjHvR7yelqzku7ZZRSylPfvQQHN8CtH0Gthk5Xc1565q6UUp44kAzfvgAX3QJtBjhdTYk03JVSqiRF+fDpvRDWAK5/vsJ+rTGGgiJ3yTuWgXbLKKVUSb75F2RthiEzISyizL/mWF4hG/bnkrz/COv2HSV5/1FG90rgD1e3rMBiLRruSil1Pmlr4PtXoNMwuPA6jw9zuQ07Mo+TvP8o6/YdIXn/UXZmneDnOXJbNKrNVa0bcVFsPa+UreGulFLnkncMPh4DdaKtGR9LcODoaaat2kdSag4b03I5WeACoH5YCJ3iwunfoTGd4sLpGBdOvZohXi1dw10ppc7GGJj/JziaCqMXQM3wc+6afTyfN5enMPXHfbiNoW3jutzcNZZO8eF0jqtPkwZhVPZS0hruSil1NslTYdMcuOoJaNLrrLvknirk7W938cEPeylwubm5Syx/7NOSmPCalVzsb5UY7iJSA/gWqG7vP8cY86SIfAhcgbUANsBoY0yyWP88vQb0BU7Z7Wu9UbxSSnlF9nZY+DAkXAaXPfibzSfzi/jghz28/e1ujucVMaBjY/7UpyXNIms7UOzZeXLmng/0NsacEJEQ4HsR+cLe9rAxZs4Z+18PtLRv3YG37J9KKeX7Ck/DnDsgpCbc9I61dJ4tr9DF1FX7eHNZCodPFtCnTRQTrm1Fm+i6DhZ8diWGu73+6Qn7YYh9O9+aqAOByfZxP4pIuIhEG2Myyl2tUkp525dPQOYmGDr7f3O0F7rczE5K4z9Ld5KRm8clLRow4doL6RJf3+Fiz82jPncRCQLWAC2AN4wxq0TkHuAZEfkb8DXwqDEmH4gB9hc7PM1uyzjjd44FxgLEx8eX93UopVT5bZkHP70LPcdDq2sBSNqbw4TZ60k9fIou8eG8NLgjvZr79tQD4OEVqsYYlzGmExALdBOR9sBjQGvgYiACa8FsjxljJhpjEo0xiZGRkaUsWymlKtjRfTBvPDTuDFc/CUDq4ZPcOTkJgPdHJ/LxPb38ItihlNMPGGOOAsuA64wxGcaSD3zAL4tlpwNxxQ6LtduUUso3uQphzhhwu6052oNDOZ5XyJ2TrGCffEc3ereOqvThjOVRYriLSKSIhNv3awLXANtEJNpuE+AGYJN9yDxgpFh6ALna366U8mnLnoW01daqShHNcLkND8xIZvehk7w5tAtNGtRyusJS86TPPRqYZPe7VwNmGWPmi8hSEYkEBEgG7rb3X4g1DDIFayjk7RVftlJKVZBdS63pBbqMhItuBuCFxdv5elsWfx/Yjl4t/KMb5kyejJbZAHQ+S3vvc+xvgPvKX5pSSnnZiSyYOw4atoLr/gXAp+vS+e83uxjaPZ4RPZo4XGDZ6RWqSqmqye2GT8ZB/jEY+SmEhpG8/yh//ngD3ZtG8NSAdn7Vx34mDXelVNW04t9Wl0z/VyCqHQdz8xg7OYlGdarz1vCuhAb793IX/l29UkqVxf6fYOnT0HYgdL2dvEIX46YkcTK/iHdHJRJRK9TpCstNz9yVUlXLyUMw53ao0xgG/BsDPPLxBtan5TJxRFdaX+B7UwmUhYa7UqrqcBXC7NHWF6l3LIKa4by1PIXPkg/w0LWtuLbdBU5XWGE03JVSVceXf4W938EN/4WYLizZkskLi7fTv0M0913VwunqKpT2uSulqob1M2DVW9D9Hug0hB2Zx7l/xjraN67HCzd39OuRMWej4a6UCnwH1sHn91vzs1/7NEdOFnDnpCTCqgczcWRXaoYGlfw7/Ix2yyilAtuJbJgxHGpFktt/InN/TGPKylQO5uYxY1wPous5v2qSN2i4K6UCl6sQM3sU7pPZvBT7H959dT0FRW46xNbj7RFdfXo+9vLScFdKBaTMY3kcnHk/HdN/YELBvSxNDWfIxTEMvjiOdo3rOV2e12m4K6UCRpHLzbLt2cz8aR/hOz/mxeAZLKh1I1cOHM9z7S+gRkjg9a2fi4a7Uson5ZwsYNKKveQVujza/1SBi8WbD5J1PJ/La+3nuZD3OB1zKf1ufxeCql7UVb1XrJTyeYUuN3dPWcNPqTlU93COl2oi9GzWgBEXhXHFNw8hcgHBQyZXyWAHDXellA96ZsFWVu/N4bXbOjGwU4znB7oKYfJAOHUExiyGWg28V6SP03BXSvmUuWvT+HDFXu64pGnpgh1g8V8g9Qe46V2I7uidAv2EJ8vs1RCR1SKyXkQ2i8j/2e1NRWSViKSIyEwRCbXbq9uPU+ztCd59CUqpQLEpPZfH5m6kR7MIHuvbunQHr5sKq9+GnuOhwy3eKdCPeNKZlQ/0NsZ0BDoB19lro/4LeMUY0wI4Aoyx9x8DHLHbX7H3U0qp88o5WcC4KWuIqBXK60O7EBJUigvoU1fC/D9B0yugz/95r0g/UuKfnrGcsB+G2DcD9Abm2O2TsBbJBhhoP8befrUE2qQNSqkKVeRy84fpa8k+kc9/h3elYe3qnh98aCfMGALhcXDLh1X2C9QzefRPo4gEiUgykAV8BewCjhpjiuxd0oCfO8digP0A9vZc4DffaojIWBFJEpGk7Ozs8r0KpZRfe2Hxdn5IOcw/bmhPx7hwzw88kQ0fDQIJgmFzICzCe0X6GY/C3RjjMsZ0AmKBbkApO8PO+jsnGmMSjTGJkZGR5f11Sik/NX/DAd7+djfDe8QzODHO8wMLTsH0W6252YfOgoim3ivSD5VqVkhjzFFgGdATCBeRn///Ewuk2/fTgTgAe3s94HCFVKuUCijbDh7j4dkb6NqkPn/r387zA90u+PhOSF8LN78HsV29V6Sf8mS0TKSIhNv3awLXAFuxQv5me7dRwGf2/Xn2Y+ztS40xpiKLVkr5v9xThYybsoY6NYJ5a1gXzxekNgYWPQrbF8D1z0Prft4t1E958s1DNDBJRIKw/jGYZYyZLyJbgBki8g9gHfCevf97wBQRSQFygNu8ULdSyo+53Ib7Z67jwNHTzBjbg0Z1a3h+8Mo3YPVEa8hj97HeK9LPlRjuxpgNQOeztO/G6n8/sz0P0EGmSqlzenXJDpZvz+YfN7Sna5NSfAm6+VP48gloOxCuedp7BQYAXYlJKVWpFm8+yH+WpjA4MZZh3eM9P3DfKpg7FuK6wY1vQzWNr/PRPx2lVKXZcuAYE2atp2NsPf4+sL3n65YeSoHpt0G9WLhtOoQE5upJFUlH+yulvG7bwWO8uWwX8zccIKJWKG8N7+r53OonD8HUQSACw+dU6cnASkPDXSnlNWv3HeHNZSks2ZpFrdAg7rq8GWMubUqjOh5+gVpwCqbdCscPwqj5ENHMuwUHEA13pVSFMsbwQ8ph3liWwsrdhwkPC+HBa1oxqmcC9cJCPP9FbhfMvQvS18CtUyDuYu8VHYA03JVSFcLtNizZmskby3exfv9RoupW54l+bRjSLZ5a1UsZNcbAggdh23y47jloM8A7RQcwDXelVLkUudzM35DBm8tT2JF5gviIMJ698SIGdY2henAZ1iw1BhY/Dms+hEsfhB73VHjNVYGGu1KqTFxuw+frD/DKkh2kHj5Fq6javHZbJ/pdFE1waabrPdOyZ+HHN6H73XD13yqu4CpGw10pVSrGGBZvzuTlr7azI/MEbaLrMnFEV/q0iaJatXLO7v39K/Dt89B5BPzun9YIGVUmGu5KKY8YY/hu5yFe/HI7G9JyaRZZi9eHdqZv++jyhzrA6ndgyVPQfhAMeE0vUionDXelVIl+2pvDC4u3s3pPDjHhNXn+5g7c1DmmfN0vxa2bCgsfggv72leflqGvXv2KhrtS6pw2pefy4pfbWb49m8g61fn7wHbcenFc2b4oPZfNn8C88dDsKrj5AwgqxXBJdU4a7kpVEYdP5FPk9mz27ezj+by5PIWFGw8SHhbCo9e3ZlTPBGqGVvAZ9Y7F1rzscd3htqkQUorZIdV5abgrVQW8sSyFFxZvL9UxtUKDuP/qloy5rCl1a3jhbHr3NzBzBFxwEQydCaG1Kv45qjANd6UC3Ipd1peg17SN4qoLG3l0THCQ0KdNFBG1Qr1T1L5VMH0INGgOw+dCjXreeZ4qrMRwF5E4YDIQBRhgojHmNRF5CrgL+Hl168eNMQvtYx4DxgAu4I/GmMVeqF0pVYJDJ/K5f0YyzRrW4tVbO5X+SlFvOJAMU2+BOhfAiE91UWsv8eSdLgImGGPWikgdYI2IfGVve8UY82LxnUWkLdbqS+2AxsASEWlljHFVZOFKqfNzuw1/mpnMsdOFTL6jm28Ee9ZWmHKjdaY+ah7UiXK6ooBV4jgmY0yGMWatff841vqpMec5ZCAwwxiTb4zZA6RwlhWblFLe9dY3u/hu5yGeHNCONtF1nS4HDm6ED/tDUCiM/NSam115TakGqYpIAtaSe6vspvEiskFE3heR+nZbDLC/2GFpnOUfAxEZKyJJIpKUnZ195malVDn8tDeHl7/awYCOjRnSLc7pcqyZHT/sD8E14PaFVl+78iqPw11EagMfAw8YY44BbwHNgU5ABvBSaZ7YGDPRGJNojEmMjIwszaFKqfM4crKAP05fR2z9mjx7YylWO/KW1JUwaaDVFaPBXmk8CncRCcEK9qnGmLkAxphMY4zLGOMG3uGXrpd0oPipQqzdppTyMmMME2av5/CJAt4Y2oU63hjCWBq7l8NHN1lfnt6xCOo3cbaeKqTEcBfrn/33gK3GmJeLtUcX2+1GYJN9fx5wm4hUF5GmQEtgdcWVrJQ6l/e+38PSbVn8pV8b2sc4PLxwx5cwdTDUT7DO2Os2draeKsaTr88vAUYAG0Uk2W57HBgiIp2whkfuBcYBGGM2i8gsYAvWSJv7dKSMUt63bt8RnvtiG79rF8XIng6fIW+ZB3PugKi2MPwTXffUASWGuzHme+BsnXYLz3PMM8Az5ahLKVUKuacKGT9tHVF1a/D8oI7O9rNvmA2fjIOYLjBsDtQMd66WKkzn1FTKzxljeOTjDWQey+P1oZ1Lt05pRVs7xVr3NL4njPhEg91BGu5K+bnJK1NZtPkgj1zXms7x9Us+wFtWv2PN7tj8Khg2G6rXca4WpeGulD/blJ7LMwu20rt1I8Zc2tS5Qn749y/zsQ+ZAaFhztWiAJ04TCm/dTyvkPumrSWiVigv3tKxYlZDKi1j4Jt/wfJ/QtsbYNC7Oh+7j9BwV8pHbEg7yp/nbCD96GmP9i9yGQpcbqbf1cN7szeej6sIFk6ANR9CxyHw+9chSCPFV+g7oZTDjDG8/8NenvtiK5G1qzOoS6zH60L3at6Qbk0dmFWx4CTMvh12LoZL/wS9/6ZrnvoYDXelHHTkZAEPzV7P19uyuKZtFC/c3IHwMAfOwkvjRBZMGwwZ66Hvi9DtLqcrUmeh4a6UQ1bvyeGP09eRc7KApwa0ZVSvBOfngSnJoRSYOgiOZ8KtU6F1X6crUueg4a5UJXO5DW8sS+HVJTuIjwhj7r29nJ8qwBP7V8O0W0EERs+H2ESnK1LnoeGuVCXKOpbHAzOTWbHrMDd0asw/bryI2r6wiEZJts6Hj8dAnWgY/rHO7OgH/OBvlVKBYfn2LCbMWs+pAhfP39yBW7rG+n43DFgXJy182JpOYOgsqNXQ6YqUBzTclfKyQpebF7/cztvf7ObCqDq8MawzLRr5wdWbbjd8/RT88Jp1cdKg9/TiJD+i4a6UF2Udy2PcR2tYt+8oQ7vH87f+bakREuR0WSUryodP74VNcyBxDPR9Aar5Qd3qfzTclfKS9KOnGfbOj2Qdz+eNoV3o1yG65IN8wakcmDUS9n4HfZ6CSx7A44H3ymdouCvlBXsPnWTYu6s4llfIlDHd6NrEgQuNyiJzM8wYCscOwE3vQIfBTlekysiTlZjiRGSZiGwRkc0icr/dHiEiX4nITvtnfbtdROTfIpJiL57dxdsvQilfsiPzOLe8vZJTBUVMv6uH/wT7lnnw7jVQmAejF2qw+zlPrhcuAiYYY9oCPYD7RKQt8CjwtTGmJfC1/Rjgeqyl9VoCY7EW0laqStiUnsutb69EgFnjevrH+HW3G5Y+A7NGWCsnjV0OcRc7XZUqpxLD3RiTYYxZa98/DmwFYoCBwCR7t0nADfb9gcBkY/kRCD9jvVWlAtKa1CMMeedHwkKDmTWuJy2j/GBETN4xmDkMvn0eOg+H0Qugrn5cA0Gp+txFJAHoDKwCoowxGfamg0CUfT8G2F/ssDS7LQOlAtSKlEPcOTmJRnWqM/WuHsSE13S6pJIdSrH61w+nwPUvWHPE6BenAcPjcBeR2sDHwAPGmGPFL74wxhgRMaV5YhEZi9VtQ3x8fGkOVcqnLN2Wyd0frSWhQRgfjelOo7o1nC6pZDuXWAtYBwXDyM+g6WVOV6QqmEdzdIpICFawTzXGzLWbM3/ubrF/Ztnt6UBcscNj7bZfMcZMNMYkGmMSIyMjy1q/Uo5auDGDcVPW0CqqNjPG9vT9YDcGvn8Vpt4M9ePhrmUa7AHKk9EyArwHbDXGvFxs0zxglH1/FPBZsfaR9qiZHkBuse4bpQLGx2vSGD9tLR1iw5nm1IIZpVFwCj6+E5Y8Ce1uhDu+hPpNnK5KeYkn3TKXACOAjSKSbLc9DjwHzBKRMUAq8PO4qYVAXyAFOAXcXqEVK+UDPvoxlSc+3cQlLRrwzshEwkJ9/JKRI6kwczgc3KgXJlURJf6NNMZ8D5zrb8HVZ9nfAPeVsy6lHOdyGzJyT7M/5zT7c06x/8gp9uecIjXnFOv2HaV360a8OayL708nsGUefDbeuj90FrS61tl6VKXw8dMNpSrHxrRcvkvJ/lWQHzh6mkLXL+MEqglE16tJXERN7r2yOQ/0aUVosA8vLVeYB1/9FVZPhMZd4Ob3IaKp01WpSqLhrqo0YwzvfLebfy3ajsttiKgVSlxEGBfF1KPvRdHER4QRVz+M+IgwosNrEBLkw2Fe3OFdMHs0HNwAPcfD1U9CsI9/J6AqlIa7qrJO5Bfx8Oz1fLHpINe3v4Bnb7yI+r7+pagnNsyG+Q9AUAgMmQkXXud0RcoBGu6qSkrJOs64KWvYc+gkj/dtzV2XNfOPhTPOp+AUfPFnWDcF4nvCoHehXqzTVSmHaLirKmfBhgz+PGc9NUKC+OjO7vRqHgArC2Vthdm3Q/Y2uOwhuPIx6wIlVWXpu6+qjCKXm38t2sY73+2hc3w4bw7rQnQ9P5gm4HyMgXUfWcvgVa8NI+ZC895OV6V8gIa7qhKyj+czftpaVu3JYWTPJjzRr61vj3TxRP5xmP8gbJwFTS+Hm96FOlElH6eqBA13FfDWpOZw79S15J4u5OXBHbmpSwD0Q6euhE/vgaOpcNUTcNmDugye+hUNdxWwjDFMXpnK0/O30Di8JnPv6UbbxnWdLqt8CvNg6dOw8g1r6oDRC6BJL6erUj5Iw10FpCMnC3jq8818lnyAq1s34uXBnagXFuJ0WeWTvhY+uRsObbcWrb7m71Y/u1JnoeGuAorbbZi9Zj/PfbGNY3lFTLimFfdd1YJq1fx4mGNRAXz3Inz7ItSOguFzocVvZv5Q6lc03FXA2HLgGE98upG1+46S2KQ+T9/QnjbRft4Nk7kFPhlnXWnacQhc9xzUDHe6KuUHNNyV3zueV8jLX+1g0oq9hIeF8sLNHRjUJda/z9bdLljxb1j2LNSoB7dOhTb9na5K+RENd+W3jDHMW3+AZxZsJftEPkO7xfPw7y4kPMzPpxA4vMvqW09bDW1+D/1fgVoBcKGVqlQa7sovpWSd4G+fbWLFrsNcFFOPd0Ym0jHOz7sr3C746V346kkIrg6D3oP2g3TedVUmGu7Kr5wucPGfpTt557vd1AgJ4ukb2jO0WzxB/twFA5CxwZrsK30NtLgGfv8fqBvtdFXKj5UY7iLyPtAfyDLGtLfbngLuArLt3R43xiy0tz0GjAFcwB+NMYu9ULeqYowxfLUlk//7fAvpR09zU5cYHru+DZF1qjtdWvkUnITl/4SVb0JYhJ6tqwrjyZn7h8DrwOQz2l8xxrxYvEFE2gK3Ae2AxsASEWlljHFVQK2qitqfc4qn5m3m621ZtIqqzcyxPejerIHTZZXf9kWw8CHI3Q9dR1vL39Ws73BRKlB4sszetyKS4OHvGwjMMMbkA3tEJAXoBqwsc4WqysovcjHxm928viyF4GrCX/q2YfQlCf6zYMa5HDsAXzwCW+dBZBu4YzHE93C6KhVgytPnPl5ERgJJwARjzBEgBvix2D5pdttviMhYYCxAfHx8OcpQgejbHdk8OW8zew6dpN9F0TzRv43/z+D48xemXz8N7kK4+m/Q8w+6QpLyirKG+1vA04Cxf74E3FGaX2CMmQhMBEhMTDQl7K6qiIO5eTw9fwsLNmbQtGEtJt/RjctbRTpdVvllrIfP74cD66wpefu9BBHNnK5KBbAyhbsxJvPn+yLyDjDffpgOxBXbNdZuU+q8Cl1uJq3Yyytf7aDIbZhwTSvGXtGM6sF+PtNh3jFY/hysegvCGuoXpqrSlCncRSTaGJNhP7wR2GTfnwdME5GXsb5QbQmsLneVKqCt3pPDXz/dxPbM4/Ru3YinBrQjvkGY02WVj9tlLaKx9Gk4mQ1db7e/MPXzsfjKb3gyFHI6cCXQUETSgCeBK0WkE1a3zF5gHIAxZrOIzAK2AEXAfTpSRp1LfpGLZxZsZfLKVGLCazJxRFeuaRvl/2uZ7v0eFj0KBzdCXA8YOgtiujhdlapixBjnu7sTExNNUlKS02WoSpSRe5p7PlpL8v6jjLm0KROubUVYqJ9fU5ezB776K2z9HOrFwTX/B+1u0i4Y5TUissYYk3i2bX7+aVL+aEXKIf4wfR15hS7eGtaF6y/y8ysx847Bdy/Bj29CtWBrZaRe4yHEz0f3KL+m4a4qjTGGt7/dzfOLttEssjb/Hd6VFo38eLEJtwuSp1pDG09mWVPyXv2kThugfIKGu6oUx/MKeWj2ehZvzqRfh2ieH9SBWtX9+K/f3u9h0WPWPOtx3WHoDIjp6nRVSv2PH3+6lL/YkXmcu6esITXnFE/0a8OYS5v675emmZth6TOwfQHUjbnQRBkAAA7ASURBVNWhjcpnabgrr/p8/QH+PGcDtaoHM+3O7v47J0zOblj2T9g4G6rXsfrVe94HoX4+ZFMFLA135RWFLjf/XLiN93/YQ9cm9XlzWBei6tZwuqzSO3YAvnke1k2BaiFwyf3WLSzC6cqUOi8Nd1Xhso7lcd+0tfy09wijeyXweN82hAb72WRfJw/DD6/A6nesL0673g6XPwR1LnC6MqU8ouGuKtRPe3O4d+paTuQV8dptnRjY6azzxvmuvGPWkMYVr0PBCeh4G1z5KNRPcLoypUpFw11VCGMMH/ywl2cXbiW2fk2mjOlG6wvqOl2W5wpPw0/vWePVT+dAmwFw1V+gURunK1OqTDTcVbmdzC/i0bkb+Xz9Afq0ieKlwR2pVzPE6bI8U3ASkj6AFf+BEwetGRt7P6HDGpXf03BX5bI7+wR3f7SGlKwTPPy7C7nniuZU84f1TE8ftfrTf3zTOlNvejkMeheaXuZ0ZUpVCA13VWaLNx/koVnrCQ4SJt3Rjcta+sG86ycPwco3rEUz8o9Bq+vgsocg7mKnK1OqQmm4q1JzuQ0vfrmdt5bvokNsPd4a3pWYcB+fRyU33ep6WfMhFOVB24Fw2QSI7uB0ZUp5hYa7KpXDJ/K5f0Yy36ccYki3eJ4c0JYaIT68oEbOHvjhVVg3FYwbOtwKl/4JIls5XZlSXqXhrjyWvP8o9360hkMnC3h+UAcGXxxX8kFOyVhvdb9snGPN1NhlpHXxUf0mTlemVKXwZLGO94H+QJYxpr3dFgHMBBKwFusYbIw5ItaEIa8BfYFTwGhjzFrvlK4qS6HLzYyf9vP051toVLc6c+/pRfuYek6X9VtuF+xYBCvfhNTvIaQW9LgHeo7XmRpVlePJmfuHwOvA5GJtjwJfG2OeE5FH7cePANdjLa3XEuiOtZB294osWFUOt9uwem8O89Yf4IuNGRw5VcgVrSJ59dZO1K8V6nR5v5Z/3Op2WfVfOLLHWijj2n9A5xG6rJ2qskoMd2PMtyKScEbzQKyl9wAmAcuxwn0gMNlYyzv9KCLhZ6y3qnyYMYaN6bnMSz7A/A0ZHDyWR82QIPq0jeL3HRvTu3UjgnxpmOPRfbDqbVg7BfJzIbYb9HkSWg+AIO1xVFVbWT8BUcUC+yAQZd+PAfYX2y/NbvtNuIvIWGAsQHx8fBnLUBVhZ+Zx5q0/wOfrD7D38ClCgoQrWjXi8X5t6NOmke8tf7d/tdWfvvVz63G7G6DHvRB71tXGlKqSyv2pNcYYESn1QqzGmInARLDWUC1vHap0juUV8tGPqcxLPsC2g8epJtCzeQPuubI517WLpl6Yj11hWngatnxmXXiUngQ16llL2XUbC/Vina5OKZ9T1nDP/Lm7RUSigSy7PR0oPoQi1m5TPiTtyClu/+AndmadoHN8OE8OaEu/DtE0quODU/Jm77DGpidPhbyj0KAF9H3RWtKuuh8v0aeUl5U13OcBo4Dn7J+fFWsfLyIzsL5IzdX+dt+yIe0oYyYlkVfoYtqd3enVoqHTJf1WUb7V5bLmQ9j7nTWPepsBkHg7JFymqx4p5QFPhkJOx/rytKGIpAFPYoX6LBEZA6QCg+3dF2INg0zBGgp5uxdqVmW0ZEsmf5i+johaoUy7szsto+o4XdKv5ey2An3dR3DqMIQ3gT5PQafhUNsPpjZQyod4MlpmyDk2XX2WfQ1wX3mLUhVv8sq9PDVvM+1j6vHuqETf6YJxFcL2LyDpfdi9DCQILrzeOktv1huq+dkiH0r5CB8bBqEqmttteHbhVt79fg992kTx7yGdnB/9Ygwc3ADJ0601SU8dshabvuov1th0veBIqXLTcA9gpwtc/GlmMos2H2R0rwT+2r+ts+PUjx+0wjx5OmRthqBQa1bGTsOg5TVQzYfnqFHKz2i4B6hDJ/K5c1IS69OO8tf+bRlzaVNnCinMg+0LrEDf9bU1eVdMojXipf0gXWhaKS/RcA9Au7JPMPqD1WQfz+etYV25rn0lL+psjHWh0fppsOkT6+rRujFwyQPWEEadkVEpr9NwDzCrdh9m7JQ1BFcTpt/Vg87x9SvniY2BzE2waS5s/sSa4yUkzBrC2HGItdKRdrsoVWk03APEnkMnWbgxg9eW7CQ2oiYfju5GfIMw7z9x1lYrzDfNhcM7rdEuTS+Hyx+Gtr+H6j423FKpKkLD3U+53IZ1+47w1dZMlmzJZFf2SQAua9mQ/wzpTHiYF2duPJQCm+dagZ69FRBIuNSaXrftQKjlgxdGKVXFaLj7kZP5RXy3M5uvtmSxbHsWOScLCK4m9GjWgBE9mnB1myjiIrx0tp6zxzpD3zwXDm602uJ7wvUvWGfodSq5X18pdV4a7j4uI/c0X2/NYsnWTFakHKbA5aZujWB6t25En7ZRXN4qkro1vDDJlzHWakbbF8K2BVZ/OlgjXX73LLS9AerFVPzzKqUqhIa7DztVUMQVzy+nwOWmSYMwRva0zs4TE+oTEuSFKzddhZD6A2xbaIV67n5AIL4HXPO01eWiy9Qp5Rc03H1YWGgwL9zSgXaN69I8sjbijQmz8k9AyhIrzHcsgrxcCK4BzXvDFY9YFxnpvC5K+R0Ndx83sJMXuj6O7LUCfcdi2P0NuPKhZgRc2A9a94PmV0ForYp/XqVUpdFwrwoK86wFo1O+hp1fWUMWwZp18eI7oXVfiOuhS9MpFUD00xyoDu+yzs5TlsCe76DoNARVt4YsXjwGWvSxFr7QudGVCkga7oEi7xikrrDmb0lZYs2NDhDRDLqMtCbmanIJhFbChU1KKcdpuPurwjzYvwr2fAN7voX0tWBcEFwTml4G3e+BFldDg+ZOV6qUckC5wl1E9gLHARdQZIxJFJEIYCaQAOwFBhtjjpSvTIWrCA6sgz3LrTDft8r6IlSCIKYrXPon67L/uO4Q4iMLcSilHFMRZ+5XGWMOFXv8KPC1MeY5EXnUfvxIBTxP1eIqtBa0SF1phXnqCig4bm2Lusj6IrTZFdZVojXqOlurUsrneKNbZiDWmqsAk4DlaLiXLP8EpP0E+36EfSsgLQkKT1nbIppDh1usM/OEy6FWA2drVUr5vPKGuwG+FBEDvG2MmQhEGWMy7O0HgaizHSgiY4GxAPHx8eUsww+dyLKDfKV1y9hg9ZlLNYhqby03F9/DOjPXZeeUUqVU3nC/1BiTLiKNgK9EZFvxjcYYYwf/b9j/EEwESExMPOs+AaOoADI3Wl96piVZZ+g5u6xtwTWs+Voue9AK89iLoUY9Z+tVSvm9coW7MSbd/pklIp8A3YBMEYk2xmSISDSQVQF1+g9jrGGI6WusW1qS1XfuKrC2146ywrzrKIjvBdEdIdiL0/MqpaqkMoe7iNQCqhljjtv3rwX+DswDRgHP2T8/q4hCfZIxcOyAFd4HkiE9yQr00/bgoJAwaNwZuo+zAj2mK9SL1QuHlFJeV54z9yjgE3syq2BgmjFmkYj8BMwSkTFAKjC4/GX6ALcLDqdYc5lnrLcC/eBGOHXY3kGgURto3R9i7SCPbKOX9CulHFHm5DHG7AY6nqX9MHB1eYpyXP5xyN5uhffBjVaQZ27+ZfRKUKgV5Bf2hQs6QHQH60vQ6rWdrVsppWxV+7Qy7xgc2mGtA5q9zbplbYNjab/sU70eXHARdBllhfgFHSDyQgjywgIZSilVQQI/3N1uOJ4BR/ZYk2n9L8y3/zrEg2tAw5bQpJcV3o3aQKO2UD9B+8iVUn4nMMLdVWStGpSz27od2Wvf32OFelHeL/sWD/FGrSHSvtVPgGpBTr0CpZSqUP4d7ju+hEWPwNF94C76pT24JkQ0tWZEbHH1L/frN4XweA1xpVTA8+9wD2tgjRNvd6MV3BHNrCCvfQFU88Iao0op5Sf8O9xju8ItHzpdhVJK+Rw9vVVKqQCk4a6UUgFIw10ppQKQhrtSSgUgDXellApAGu5KKRWANNyVUioAabgrpVQAEmOcX+FORLKx5n4vi4bAoQosxxcF+msM9NcHgf8a9fU5o4kxJvJsG3wi3MtDRJKMMYlO1+FNgf4aA/31QeC/Rn19vke7ZZRSKgBpuCulVAAKhHCf6HQBlSDQX2Ogvz4I/Neor8/H+H2fu1JKqd8KhDN3pZRSZ9BwV0qpABQQ4S4iT4lIuogk27e+TtdUEUTkOhHZLiIpIvKo0/V4g4jsFZGN9vuW5HQ95SUi74tIlohsKtYWISJfichO+2d9J2ssr3O8xoD5DIpInIgsE5EtIrJZRO632/3qfQyIcLe9YozpZN8WOl1MeYlIEPAGcD3QFhgiIm2drcprrrLfN78aR3wOHwLXndH2KPC1MaYl8LX92J99yG9fIwTOZ7AImGCMaQv0AO6zP3t+9T4GUrgHmm5AijFmtzGmAJgBDHS4JlUCY8y3QM4ZzQOBSfb9ScANlVpUBTvHawwYxpgMY8xa+/5xYCsQg5+9j4EU7uNFZIP9X0af/u+Sh2KA/cUep9ltgcYAX4rIGhEZ63QxXhJljMmw7x8EopwsxosC7TOIiCQAnYFV+Nn76DfhLiJLRGTTWW4DgbeA5kAnIAN4ydFiVWlcaozpgtX9dJ+IXO50Qd5krLHHgTj+OOA+gyJSG/gYeMAYc6z4Nn94H4OdLsBTxpg+nuwnIu8A871cTmVIB+KKPY612wKKMSbd/pklIp9gdUd962xVFS5TRKKNMRkiEg1kOV1QRTPGZP58PxA+gyISghXsU40xc+1mv3of/ebM/XzsP+if3QhsOte+fuQnoKWINBWRUOA2YJ7DNVUoEaklInV+vg9cS2C8d2eaB4yy748CPnOwFq8IpM+giAjwHrDVGPNysU1+9T4GxBWqIjIF67+DBtgLjCvWN+a37OFkrwJBwPvGmGccLqlCiUgz4BP7YTAwzd9fo4hMB67EmiI2E3gS+BSYBcRjTW092Bjjt19InuM1XkmAfAZF5FLgO2Aj4LabH8fqd/eb9zEgwl0ppdSvBUS3jFJKqV/TcFdKqQCk4a6UUgFIw10ppQKQhrtSSgUgDXellApAGu5KKRWA/h8GgTMxCePdIAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Tinc = 3\n",
    "Tinf = 15\n",
    "startI = 0.00005\n",
    "Rt = 2.9\n",
    "Toffset = 5\n",
    "tmax = ddf['days'].max()+1+Toffset\n",
    "t = np.linspace(1,tmax,tmax)\n",
    "config = {'Rt' : Rt, 'Tinc': Tinc, 'Tinf': Tinf, 'beta': 0.25, 'gamma': 0.25}\n",
    "SEIR_y0 = [1-startI,startI/2,startI/2,0]\n",
    "\n",
    "\n",
    "modelOutput = odeint(SEIR_model, SEIR_y0, t, args=(config,), atol=1e-12, rtol=1e-12)\n",
    "modelSums = pandas.DataFrame(modelOutput[:,[1,2,3]]).sum(axis=1)\n",
    "\n",
    "plt.plot(ddf['days'], ddf['confirmed'])\n",
    "plt.plot(t-1-Toffset, modelSums[t-1] * population )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe10cad9850>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hUZfr/8fedQu8QWighEKUTIHQLoqKiLoiIoEgX3XVtq+vKumtbXdG1ru5XRUFAqgoIsqAIoqDUAKGEGkogIYQUII20yfP74wy/jQqmTs7Myf26rlyZOXNmco9jPjk85z7PI8YYlFJKOYuf3QUopZQqfxruSinlQBruSinlQBruSinlQBruSinlQAF2FwDQqFEjExISYncZSinlU7Zv355sjAm61GNeEe4hISFERkbaXYZSSvkUEYm93GM6LKOUUg6k4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg7kFX3uSilVmbgKDPtOpfHTkWS6Btelf7tG5f4zNNyVUsrDjDEcScpk45FkfopJZvPRVM5fyAPg9wPbargrpZSvSDh/gZ9iUtgYk8zGIymcTssGILhedW7q1IT+bRvRv21DGtep5pGfr+GulFLlJP7cBab/cIQNh5M5mpwJQIOaVejXtiED2jZiQLuGtGpQAxHxeC0a7kopVUb5rgJmb4rljdUHcRUY+rdtyD19WtG/bSPaN62Nn5/nw/yXNNyVUqoM9sSdZ+rS3eyNT+O6K4N4cWhnWjaoYXdZGu5KKVUaGTn5vLn6ELM2HqNhrar8554eDOnStEKGXIqjyHAXkWrAeqCqe/8vjDHPicgs4FrgvHvX8caYKLHe2TvAECDLvX2HJ4pXSik7fLsvkeeW7SUhLZt7+7Tizze1p271QLvL+pniHLnnAIOMMRkiEgj8KCKr3I/92RjzxS/2vwUIc3/1Ad53f1dKKZ92+nw2zy+P5uvo01zZpDbv3tODnq3r213WJRUZ7sYYA2S47wa6v8xvPGUoMMf9vM0iUk9EmhljEspcrVJK2cBVYJi7OZZ/fXOQPFcBT918JfdfHUqgv/de5F+sykTEX0SigDPAt8aYLe6HXhaR3SLylohUdW8LBk4Wenqce9svX3OKiESKSGRSUlIZ3oJSSnnO8eRMhr+/keeWR9O9VT2+ffxa/jCwnVcHOxQz3I0xLmNMONAC6C0inYGpQHugF9AA+EtJfrAxZroxJsIYExEUdMklAJVSylY5+S4enLud48mZvDMqnDkTe9Oqof2dMMVRoj89xphzwDrgZmNMgrHkAJ8Avd27xQMtCz2thXubUkr5lHfWHObA6XTeHNmNoeHBXtMJUxxFhruIBIlIPfft6sCNwAERaebeJsAwYK/7KcuBsWLpC5zX8XallK/ZHnuWD344wsiIFlzfoYnd5ZRYcbplmgGzRcQf64/BZ8aYFSLynYgEAQJEAQ+691+J1QYZg9UKOaH8y1ZKKc/Jys3nyc930axudf5+W0e7yymV4nTL7Aa6X2L7oMvsb4CHyl6aUkrZ49VVBziWnMn8+/tQu5p39a8Xl3ef7lVKqQr2U0wyszfFMmFACP3blv9UvBVFw10ppdzSsvP48+e7CA2qyV9ubm93OWWic8sopZTbi1/tIzE9h8W/70+1QH+7yykTPXJXSims+WK+2B7HHwa2JbxlPbvLKTMNd6VUpZeSkcPUJbvp1LwODw8Ks7uccqHDMkqpSs0Yw9++3EvahXzmTQ6nSoAzjnmd8S6UUqqUlkWdYtXe0/xp8BVc2bS23eWUGw13pVSldfp8Ns8u20vP1vW5/+pQu8spVxruSqlKyRjDU4t3k+cyvHFXN/xtWOfUkzTclVKV0rwtJ1h/KIm/DmlPSKOadpdT7jTclVKVTmxKJv9cuZ+rwxoxpm9ru8vxCO2WUUpVGjn5Lr7cGc9762Lw9xNeG9HVp6bxLQkNd6WU42Xk5LNw6wk+3nCM02nZdGpeh3+N6EazutXtLs1jNNyVUo6VmpnLrJ+OMXtTLOcv5NEvtCGvjejK1WGNHHvEfpGGu1LKceLPXeCj9UdZuO0E2XkFDO7YhAcHtqVHq/p2l1ZhNNyVUo5xODGdD344yrIoa2XPoeHBPHhtKGFNnHNxUnFpuCulvNLJ1CxeX32QzBxXsfbPyMlj89FUqgX6MaZvayZf3YYW9X1jMWtP0HBXSnkdV4Hh8UVRRJ9Ko00xe9D9/OCRQe0YP6ANDWpW8XCF3k/DXSnldWZvPE5k7FneHNmN4T1a2F2OT9KLmJRSXiU2JZPXvjnAoPaNuaN7sN3l+Kwiw11EqonIVhHZJSLRIvKCe3sbEdkiIjEiskhEqri3V3Xfj3E/HuLZt6CUcoqCAsNTX+wm0M+Pl+/o7Ph2RU8qzpF7DjDIGNMNCAduFpG+wKvAW8aYdsBZYJJ7/0nAWff2t9z7KaVUkeZtPcGWY6n87bYOjr7AqCIUGe7GkuG+G+j+MsAg4Av39tnAMPftoe77uB+/XvTPr1KqCHFns5jmnu9lZERLu8vxecUacxcRfxGJAs4A3wJHgHPGmHz3LnHAxcGxYOAkgPvx80DD8ixaKeUsxhimLtkDwLQ7nTvfS0UqVrgbY1zGmHCgBdAbaF/WHywiU0QkUkQik5KSyvpySikftmjbSTYcTmbqkA4E19PhmPJQom4ZY8w5YB3QD6gnIhdbKVsA8e7b8UBLAPfjdYGUS7zWdGNMhDEmIigoqJTlK6V8XcL5C7z83/30C23IPb1b2V2OYxSnWyZIROq5b1cHbgT2Y4X8CPdu44Bl7tvL3fdxP/6dMcaUZ9FKKWe4OByTX2B49c6u+DlsNSQ7FecipmbAbBHxx/pj8JkxZoWI7AMWishLwE5ghnv/GcCnIhIDpAKjPFC3UsoBFu+I5/uDSTx3e0daNay8UwV4QpHhbozZDXS/xPajWOPvv9yeDdxVLtUppRwrMS2bF7+KpldIfcb1C7G7HMfRK1SVUhXOGMMzS/eSk1/AayO66XCMB2i4K6Uq3PJdp1izP5EnB19Z7InBVMlouCulKlRSeg7PLY+me6t6TLyqjd3lOJaGu1KqQj27bC9ZOS7+NaIr/joc4zEa7kqpCmGMYc6m46zae5pHbwijXePKtzpSRdL53JVSHrfr5Dle+CqaHSfO0Te0AQ9cE2p3SY6n4a6U8pjEtGxe+/ogi3fE0ahWVV67sysjerbQ7pgKoOGulCp32XkuZvx4jP+siyHfZXjw2rY8dF1balcLLNkL5edCgC6ZVxoa7kqpcmOMYdXe0/xz5X7izl5gcMcmPHNrB1o3LEW7Y+xGWPIA3LMImnQs/2IdTsNdKVUuok+d58Wv9rHlWCrtm9Zm/uQ+9G/XqHQvlnYKPhsH1epAXV1qrzQ03JVSZZKckcMbqw+ycNtJ6lUP5KVhnRnVqyUB/qVsxsvPgc/GQm4mjPsKqtUt34IrCQ13pVSpFBQYFm47ybRV+8nKdTGhfxsevT6MujVKOK7+S18/DXHb4K7Z0LjMS0dUWhruSqkSO5SYzl+X7CEy9ix9Qxvw0rAutGtcq+wvvONTiJwJAx6DTsOK3l9dloa7UqrYsvNcvPvdYT784Si1qgXwrxFWa2O5LIsXvx3++wSEDoRBfy/761VyGu5KqWL58XAyz3y5h9iULIb3COaZIR1oWKtq+bx4RhIsGgu1msCIT8Bfo6ms9L+gUuo3JWfk8PJ/97N0ZzwhDWswb3IfBpS2C+ZSXPnwxQTISoaJ30CNBuX32pWYhrtS6pKMMXweGcc/V+0nMyefhwe146Hr2lEt0L98f9Ca5+D4Bhj2ATQPL9/XrsQ03JWqJC7kunAVcznj+LMX+PuyvWw9lkpE6/q8MrwLYU08MNHXni9g03vQewqEjy7/16/ENNyVqgQWbD3BM0v3UFCCperrVAvgleFduDuipWfmgkmMhuUPQ8u+MPjl8n/9Sk7DXSmHO30+m5dW7KNHq/rc1KlpsZ4T4C/c2rUZjWtX80xRF87Cwnuhah0YOUfnj/EADXelHO7FFdHkFxjeHBlOq4Y17C4HCgpgyRQ4HwcTVkLtJnZX5EhFXh8sIi1FZJ2I7BORaBF51L39eRGJF5Eo99eQQs+ZKiIxInJQRG7y5BtQSl3eugNnWLnnNI9cH+YdwQ7wwzQ4vBpumQYte9tdjWMV58g9H3jCGLNDRGoD20XkW/djbxljXi+8s4h0BEYBnYDmwBoRucIY4yrPwpVSv+1Crotnl++lXeNa3H+1lyyOcWAl/PAqhI+BiEl2V+NoRR65G2MSjDE73LfTgf3Ab03TNhRYaIzJMcYcA2IA/fOsVAV7b91hTqZe4KVhnakS4AUrap6KgsWToXl3uPV1KI+rWtVllegTF5EQoDuwxb3pjyKyW0Rmikh997Zg4GShp8VxiT8GIjJFRCJFJDIpKanEhSulLu9wYjrT1x/lzh4t6Bva0O5yrPH1+XdbFyiNXgSB1e2uyPGKHe4iUgtYDDxmjEkD3gfaAuFAAvBGSX6wMWa6MSbCGBMRFBRUkqcqpX6DMYZnvtxLjSoB/HWIF8yqmJ0G80ZCXhbc85meQK0gxQp3EQnECvZ5xpglAMaYRGOMyxhTAHzE/4Ze4oGWhZ7ewr1NKVUBFu+IZ+uxVKbe0r785n4pLVcefD4ekg/CyNm6olIFKk63jAAzgP3GmDcLbW9WaLc7gL3u28uBUSJSVUTaAGHA1vIrWSl1OWczc/nnyv30bF2fkREti36CJxkDK5+EI2vhtreg7SB766lkitMtMwC4D9gjIlHubX8FRotIOGCA48ADAMaYaBH5DNiH1WnzkHbKKFUxXv36AOcv5PHyHZ09c1VpSWz8N2yfBVf9CXqMtbeWSqjIcDfG/Ahc6v+Slb/xnJcBvZ5YqQoUeTyVhdtO8sA1obRvWsfeYqK/hG+fhU7DdW52m3hBf5RSqqzyXAU8s3QvwfWq8+gNYfYWc3IbLH3AmjNm2PvgpzFjB51+QCkHmPnjMQ4mpvPR2AhqVLHx1zr1GCwYBbWbwaj5EOihuWlUkfRPqlI+Lu5sFm+vOcyNHZtwY0cb2wyzUmHeXWBccO8XUNML+usrMT1yV8rHPb98n/X9d53sKyI/BxbdB+diYewyaNTOvloUoEfuSvm01dGnWbM/kcdvDCO4nk1XfRoDyx+B2B9h6P9B6/721KF+RsNdKR+VmZPP88ujad+0NhMGtLGvkO+nwe6FcN3foOtd9tWhfkaHZZTyEueycvn7smhOpGYVa/+0C3mcOp/Nu/d0J9DfpuO0LR9aU/iGj4FrnrSnBnVJGu5KeYG4s1mM/2QbJ1Ky6Nu24SUvLPmletUDGd8/hJ6tG3i8vkuKmg+rnoL2t8Ht7+gsj15Gw10pm+07lcb4T7ZyIc/FnEm9vWMWx6LsWw7LHoLQgTBiJvhrlHgb/USUstFPMck88Ol2alcL4IsH+3Nl09p2l1S0mLXwxUQIjrB62QNsnpxMXZKGu1I2+XJnPH/+YhehjWoxa2IvmtX1gTnOT2y2FrYOag/3fg5VatpdkboMDXelKpgxhg/XH2XaqgP0DW3Ah/dFULd6oN1lFS1hl3WRUt1guG8pVK9nd0XqN2i4K1WBXAWGf6zYx6yNx7m9W3Nev6srVQP87S6raEmH4NPhULUO3Pcl1NIFdrydhrtSFSQ7z8VjC6P4Ovo0U64J5emb29s/LW9xnDsBnw6zumHGLoN6Ns8Tr4pFw12pCnAuK5f750QSGXuWZ2/ryMSrbLzoqCTSE2HOUMjNgPErdVoBH6LhrpSHFe5hf290D27t2qzoJ3mDrFTriD090Tpib9rZ7opUCWi4K+VBR5MyGDV9M9m+1MMOkJMO80ZASozVFdOyl90VqRLScFfKQ1Iychj/yTZcBYYvft+fK5r4QA87QG4WLBgNp6Lg7k+tC5WUz9FwV8oDsvNcTJ4TSWJaNgum9PWhYM+E+XfD8R9h+HRof6vdFalS0nBXqpwVFBgeXxRF1MlzvH9vD3q0qm93ScWTkw7zRsLJzXDHh9B1pN0VqTIocio5EWkpIutEZJ+IRIvIo+7tDUTkWxE57P5e371dROTfIhIjIrtFpIen34RS3uSVVftZtfc0zwzpwM2dfeTk6YVz8OkdcHIL3DkDut1td0WqjIozT2g+8IQxpiPQF3hIRDoCTwNrjTFhwFr3fYBbgDD31xTg/XKvWikvNWfTcT7acIzx/UOY5CvtjlmpVrvjqSgYOQc6D7e7IlUOigx3Y0yCMWaH+3Y6sB8IBoYCs927zQaGuW8PBeYYy2agnoj4yOGLUqW3Zl8izy+P5oYOTfj7bR0RX5gCNzMFZv8OzuyDUfOgw212V6TKSYlm+BeREKA7sAVoYoxJcD90Gri4Mm8wcLLQ0+Lc2375WlNEJFJEIpOSkkpYtlLeZXfcOR5esJPOwXX59+hw/H3hytOMMzDrVkg5DKMXwhU32V2RKkfFDncRqQUsBh4zxqQVfswYYwBTkh9sjJlujIkwxkQEBek8Fcp3nUzNYuKsSBrWqsKMcb2oUcUH+hTSEqxgPxcL93wG7a63uyJVzooV7iISiBXs84wxS9ybEy8Ot7i/n3FvjwcKTz7Rwr1NKcc5fyGPCbO2kZvvYtaEXgTV9oG5zc/HwawhkHYKxiyG0Gvtrkh5QHG6ZQSYAew3xrxZ6KHlwDj37XHAskLbx7q7ZvoC5wsN3yjlGLn5BTz46XZiUzKZPjaCdo19oJf9bCx8cgtkJlvT9rbub3dFykOK8+/HAcB9wB4RiXJv+yswDfhMRCYBscDFptiVwBAgBsgCJpRrxUp5AWMMTy/ezaajKbx9d7hvTCuQcsQ6eZqbYc0VE6xdyk5WZLgbY36Ey67X+6uBOvf4+0NlrEspr/bWmsMs2RnPk4OvYFj3X/ULeJ+kg1a7oysXxn0FzbraXZHyMB8486OUvc5m5nIkKcP9lcmhxHS+P5jE3REteeg6H5gC98QWmD8S/KvAuBXQpKPdFakKoOGuFNYKSXFns6wAP5P5szBPzcz9//tVDfCjTaOajO8fwjO3dvD+Xvb9K2DxJKgTbJ08beAjF1apMtNwV5XemfRsxs/cxr6E/3X4NqpVhdCgWtzUqSltg2rStnEt2gXVonm96r7Rww6wbQasfBKad7faHWs2srsiVYE03FWlFnc2izEfb+FMeg4v/K4TnYPrENqoFvVrVrG7tNIzBta9DOv/BWE3wV2fQJWadlelKpiGu6q0jiRlMObjLWTm5PPppD70bO0jszf+FlcerHgMds6F7vfBbW+Dv/6aV0b6qatKKfrUecbO2IoILJzSj47N69hdUtnlZsLn4+Hwarj2LzBwqrWotaqUNNxVpbM9NpXxn2yjdtUA5k7uQ2hQLbtLKrvMZJh3FyREWUfrEXp5SWWn4a4qlQ2Hk5gyZztN61Zj7uQ+BNerbndJZZd6FObeaU0ncPc8aD/E7oqUF9BwV5XG13tP88iCnYQG1eTTSX18Yx6YosTvsHrYC/Kti5Na9ra7IuUlSjTlr1K+avH2OB6av4NOwXVYNKWfM4L98BqYdRsEVIeJqzXY1c9ouCvHm73xOE98vou+oQ2YO6kPdWsE2l1S2RgDG9+D+XdBg1CYtBqCrrC7KuVldFhGOZYxhv+si+H11Ye4sWMT3h3dnWqB/naXVTZ52bDicdg1HzrcDsM+gKoOOCGsyp2Gu3IkV4Hh1a8PMH39UYZ3D+a1EV0J8Pfxf6imn4ZFYyBuG1z7tNXu6Ofj70l5jIa7cpzEtGweXxTFxiMp3Ne3NS/8rhN+vjJlwOXE74CF90L2ObhrNnQaVvRzVKWm4a4cZd2BMzzx+S4u5Lp47c6u3BXRwvsn9yrKni9g2UNQMwgmfqPT9api0XBXjpCbX8BrXx/g4x+P0b5pbd67p7tvrIz0Wwpc8N0/4Me3oFV/GDkHaul6w6p4NNyVzzuenMnDC3ayJ/48Y/u15q9DOvj+idPsNFhyPxz6GnqMgyGvQ4APT2amKpyGu/JpX+6M55mlewjw9+ODMT25uXNTu0squ5QjsGA0pMRYod5rss4Ro0pMw135pMycfJ5dFs3iHXH0CqnP26O6O2MqgSPrrMm/RKwFrEOvtbsi5aM03JXPiT51nofn7+RYSiaPXB/GI4Pa+X6bY0EBbHgDvv8nNLoSRi/QVZNUmWi4K59hjGHWxuO8svIA9WsGMn9yX/q1bWh3WWWXkWSNrx9dB51HwO1vQ1UfPxmsbFfk4Y6IzBSRMyKyt9C250UkXkSi3F9DCj02VURiROSgiNzkqcJV5ZKamcv9cyJ54at9XB3WiFWPXuOMYD/+I3xwFcRuhNvfgTs/1mBX5aI4R+6zgPeAOb/Y/pYx5vXCG0SkIzAK6AQ0B9aIyBXGGFc51Koqqc1HU3hsYRSpmbk8e1tHJgwI8f3e9YIC+PENWPdPa36YMV9A0y52V6UcpMhwN8asF5GQYr7eUGChMSYHOCYiMUBvYFOpK1SVVr6rgH9/F8N73x2mdcOaLBnXn87Bde0uq+wyk61hmCPf6TCM8piyjLn/UUTGApHAE8aYs0AwsLnQPnHubb8iIlOAKQCtWrUqQxnKiU6du8BjC6PYejyV4T2CeXFoZ2pVdcApouM/weJJkJVqrZjUc7y2OSqPKG2LwftAWyAcSADeKOkLGGOmG2MijDERQUF61Z36n9XRp7nlnQ3sPXWeN0d2482R4b4f7Be7YWbfBoE14P611lJ4GuzKQ0r1G2OMSbx4W0Q+Ala478YDLQvt2sK9TakiZee5eGXlfmZviqVzcB3eHd2DNo1q2l1W2WUmw5IpcGQtdL7TOnGqwzDKw0oV7iLSzBiT4L57B3Cxk2Y5MF9E3sQ6oRoGbC1zlcrxYs5k8PCCnexPSGPigDb85ZYrqRrg41MIABz9Hpb+HrJSdBhGVagiw11EFgADgUYiEgc8BwwUkXDAAMeBBwCMMdEi8hmwD8gHHtJOGfVbjDF8vj2O55ZFU72KPzPHRzCofRO7yyq73CxY+wJs+QAahsE9i3Q2R1WhxBhjdw1ERESYyMhIu8tQFcxVYPjbl3tZsPUE/UIb8vaocJrUqWZ3WWUXvx2WPAAph6HP7+GG5yDQAVMjKK8jItuNMRGXeszHz1IpX5WbX8CfPotixe4E/jCwLU8MvhJ/X19Qw5UH61+H9f+C2k1h7DIIHWh3VaqS0nBXFe5Cros/zNvOuoNJPDOkA/dfE2p3SWWXdNA6aZoQBd1Gw83ToHo9u6tSlZiGu6pQ6dl5TJodybbjqbwyvAuje/v4NQ4FBbD1Q1jzvNXiOHIOdBxqd1VKabiripOamcv4T7ay71Qa/x7Vndu7Nbe7pLI5dxKW/QGOrYcrbobb/w21HXAyWDmChruqEIlp2Yz5eAsnUrOYPranb3fEGAO7FsKqp8AUwO/ehe73aYuj8ioa7srjTqRkce+MzaRm5DJ7Ym/6hvrwbI7nTsLKJ63l71r1g2Hv67zryitpuCuPOpyYzpgZW8jJL2D+/X3p1tJHTzIWuGDrdFj7D8DA4Jeh7+/BzwEXWilH0nBXHrMn7jxjZ24h0N+PRVP6cWVTH73kPmE3fPUInNoJ7W6EW9+A+q3trkqp36Thrjxiy9EUJs2OpF6NQOZN7kPrhj44R0xuJnw/DTb9B2o0gBEzodNwHVtXPkHDXZW7dQfP8OCn22nZoAZzJ/WhaV0fvOo0Zg2seBzOnYAeY+GGF6yAV8pHaLircrVi9ykeWxhF+2a1mT2hNw1rVbW7pJLJSIJvpsKez605YcavhJABdlelVIlpuKtys2jbCaYu2UPP1vWZMb4XdaoF2l1S8RkDO+fC6r9ZwzHXPg1X/wkCfOyPk1JuGu6qXHy84Sgv/Xc/114RxAdjelK9ig91kSTshq+fhtifrPbG29+BoCvtrkqpMtFwV2VijOHtNYd5Z+1hhnRpytt3d6dKQGkX+Kpgmcnw3UuwYzZUq2eFevex4Ocj9Sv1GzTcVakZY/jHiv3M/OkYd/VswSvDuxDg7wPB6MqDrR9ZnTC5GdD7ARj4F6he3+7KlCo3Gu6qVFwFhqlLdvNZZBwTB7Thb7d2wM8XpuyNWQNfT4XkQ9B2ENz0CjRub3dVSpU7DXdVYrn5BTy2aCcr95zm0evDeOyGMMTbe79TjsA3z8ChVdAgFEYvtCb78va6lSolDXdVIhdyXTw4dzs/HErib7d2YPLVXj4Xe3aatXjG5vchoBrc+CL0eVC7YJTjabirYkvLzmPyrEi2xaby6p1duLuXF8/FXuCCqPmw9kXIPAPhY+D6Z3VKXlVpaLirYknNzGXszC0cSEjn3dHdua2rl87FbgwcXGmFetIBaNEb7lkIwT3trkypClVkuIvITOA24IwxprN7WwNgERACHAdGGmPOijXw+g4wBMgCxhtjdnimdFUR8lwFrI5O5I1vDxJ/9gIfjY3guvaN7S7r0o5tgLUvQNw2aNgO7pptrYqk4+qqEirOkfss4D1gTqFtTwNrjTHTRORp9/2/ALcAYe6vPsD77u/Kx5xJy2b+1hMs2HqCxLQcWtSvzpyJvenjjXOxn4qyjtSPrIXaza3FM7rdA/76D1NVeRX5f78xZr2IhPxi81BgoPv2bOB7rHAfCswxxhhgs4jUE5FmxpiE8ipYeY4xhq3HUpmzOZZv9p4mv8BwzRVBvDysNde1b4y/t7U6phyxLkKKXmL1qA9+CXpNhsDqdlemlO1Ke2jTpFBgnwYunqUKBk4W2i/Ove1X4S4iU4ApAK1aefGJuUogMyefpTvjmbs5lgOn06lTLYBx/UMY07c1bRp54VS9aQnww6uwY47V9XLNn6H/w1Ctrt2VKeU1yvzvVmOMERFTiudNB6YDRERElPj5quxizmQwd3Msi7fHkZ6TT8dmdZg2vAtDw4O9c26YzBTY+G/Y8iEU5EOvSVaw1/LScwBK2ai04Z54cbhFRJoBZ9zb44GWhfZr4d6mvMx73x3m9dWHCPQXhnRpxth+renRqr53XoyUngib3oVtMyEvC7qOhIFTde1SpX5DacN9OTAOmOb+vqzQ9j+KyEKsE6nndbzd+3y84Sivrz7E7fCV0CoAAArqSURBVN2a8+xtHQmq7aUX9KSdgp/ege2zwJULnUfA1U/odAFKFUNxWiEXYJ08bSQiccBzWKH+mYhMAmKBke7dV2K1QcZgtUJO8EDNqgzmbo7lpf/uZ0iXprw1spt3TvR1NhZ+fAui5oEpgG6j4Ko/QcO2dlemlM8oTrfM6Ms8dP0l9jXAQ2UtSnnGkh1x/H3ZXga1b8zbd3f3vmBPOQIb3oTdC0H8oPsYGPCYLkatVCloI3AlsXJPAk9+vot+oQ35v3t7eNec62f2w4Y3YO9i8K8Cve6HAY9AHS+9ClYpH6DhXgl8dyCRRxbspHur+nw0NoJqgV7QCWMMnNgEm/4DB1ZAYE3o90erpVG7X5QqMw13h/spJpkH5+6gQ7M6fDKhFzWr2vyRu/Jg3zLY9B6c2gnVG8A1T0Hf30ONBvbWppSDaLg7WOTxVCbPjqRNw5rMmdjb3gWrL5yzlrPb8iGkxVtzv9z6JnQbDVVq2FeXUg6l4e5Qu+POMeGTbTSrW41PJ/emfs0q9hSSeswK9J2fWkvahVxthXrYYF2rVCkP0nB3oIOn0xk7cyt1qgcyd3IfGteuVrEFGAMnt1pDLwdWWJ0vnUdAvz9As24VW4tSlZSGu8McTcrg3o+3UDXAj/n396F5vQqcRCvvAuxdApEzIH67NdfLgEeh9xTtfFGqgmm4O0hsSib3frwFYwzzJveldcMKmvQr+TBEzrQuOso+D42ugFv+BeH3QNVaFVODUupnNNx9XHaei9X7ElmyI44Nh5OpWcWfhVP60a5xbc/+YFeeNeQSOROOrQe/AOhwO0RMgpCrdIEMpWym4e6DCgoMW4+nsmRHHKv2nCY9J5/mdavxwDWhjO7dipYNPNh9cu6k1fWyYw5kJELdVjDo79D9Pl2fVCkvouHuQ44kZbB0RzxLd8YTf+4CNav4c0uXZgzvEUzfNg3x89RiGgUuOPIdbJsBh7+xTphecRNETIR2N4CfF1wUpZT6GQ13L5eamcuK3adYvCOeXSfP4SdwVVgQT918JYM7NvXsvOtJByFqPuxeBOkJUDMIrnoceo6HerrAilLeTMPdi2Xk5NN/2lqy8wpo37Q2zwzpwNDw5jSu48HWxgtnrY6XqPkQHwnib/Wk3/IqXHELBNjUL6+UKhENdy9Wq2oAz93eiW4t6tGxeR3P/SBXPhxdZ3W7HFgJrhxo3AkGv2wtjKFzvSjlczTcvdzo3h4c/jhzAHbNh12LIOO0Nc9LxASrhbFpV+14UcqHabhXNmdjIXqJNb3u6T1WC2PYYCvQw27SYRelHELDvTJIOwXRX1qBHh9pbQuOgJtegS53Qa0ge+tTSpU7DXenykiC/cusk6OxGwFjDbXc8Dx0ugPqh9hbn1LKozTcnSQrFQ781zpCP/aDtf5ooyth4FToPBwahdldoVKqgmi4+7pzJ+HgSmsqgOM/gXFB/TZWP3rnO6FxRz0xqlQlpOHua4yBM/usI/QDKyBhl7U9qL01A2OH26F5dw10pSq5MoW7iBwH0gEXkG+MiRCRBsAiIAQ4Dow0xpwtW5mVXIELTm75X6CfPQ4ItOwNN74IV94KjdrZXaVSyouUx5H7dcaY5EL3nwbWGmOmicjT7vt/KYefU7lcOAtHv4fDa+DQKshKAf8qEDrQGnK54hadqEspdVmeGJYZCgx0354NfI+Ge9EKCuD0LohZYwV63DZr/LxaXasPvf2t1iRdVT08la9SyhHKGu4GWC0iBvjQGDMdaGKMSXA/fhq45OGliEwBpgC0alVJJ6HKSrVmW4xZY31lJlnbm4XD1X+CdjdCcE/w11MjSqmSKWtqXGWMiReRxsC3InKg8IPGGOMO/l9x/yGYDhAREXHJfRwnPxdO7bCGW2LWWEvRmQKoXh/aXg9hN0LbQTqXi1KqzMoU7saYePf3MyKyFOgNJIpIM2NMgog0A86UQ52+qcBldbMcW299ndgEeVmAQHAPuOYpK9Cbd9c50ZVS5arU4S4iNQE/Y0y6+/Zg4EVgOTAOmOb+vqw8CvUJBQWQtP9/YX78J8g5bz0W1B66j4GQq61l6Go0sLdWpZSjleXIvQmwVKx+6gBgvjHmaxHZBnwmIpOAWGBk2cv0Uq58OBMNJ7ZA7E9w/EfIcjcO1W8DnYZBm2usQNfOFqVUBSp1uBtjjgLdLrE9Bbi+LEV5rezzEBdp9Zyf2GyNmedmWI/VCbaGWC6Geb2W9taqlKrUtA3jcoyBc7FwcqsV5Ce3QGI0YED8oEkn6DYaWvWFln2gbgu9KlQp5TU03MEK8vQEOBUFp3ZCQpR1O9N9LrhKbWgRAQOftq4KDY6Aah5cGUkppcqo8oV74SBPcId54SAXP2smxXbXWz3mLftYR+nazaKU8iHODvfcTEg6CEkH4Mx+63vCLshItB4vHOTNwq2WxKadoUpNe+tWSqkycka452ZB8qGfh/iZ/XDuBNZFtFjzsjQMsy4SahYOzcOhaRcNcqWUI/l2uB/6BlY9Za0LejHE/QKtRSmCe0L4vdC4PQR1gAahehm/UqrS8O20q9nIOgrvNtq6SKjxxRAPtLsypZSylW+He3BPGDnb7iqUUsrr+NldgFJKqfKn4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg6k4a6UUg4kxti/NrWIJGGt2lQajYDkcizHGzn9PTr9/YHz36O+P3u0NsYEXeoBrwj3shCRSGNMhN11eJLT36PT3x84/z3q+/M+OiyjlFIOpOGulFIO5IRwn253ARXA6e/R6e8PnP8e9f15GZ8fc1dKKfVrTjhyV0op9Qsa7kop5UCOCHcReV5E4kUkyv01xO6ayoOI3CwiB0UkRkSetrseTxCR4yKyx/25RdpdT1mJyEwROSMiewttayAi34rIYff3+nbWWFaXeY+O+R0UkZYisk5E9olItIg86t7uU5+jI8Ld7S1jTLj7a6XdxZSViPgD/wFuAToCo0Wko71Vecx17s/Np/qIL2MWcPMvtj0NrDXGhAFr3fd92Sx+/R7BOb+D+cATxpiOQF/gIffvnk99jk4Kd6fpDcQYY44aY3KBhcBQm2tSRTDGrAdSf7F5KHBxPcjZwLAKLaqcXeY9OoYxJsEYs8N9Ox3YDwTjY5+jk8L9jyKy2/1PRq/+51IxBQMnC92Pc29zGgOsFpHtIjLF7mI8pIkxJsF9+zTQxM5iPMhpv4OISAjQHdiCj32OPhPuIrJGRPZe4mso8D7QFggHEoA3bC1WlcRVxpgeWMNPD4nINXYX5EnG6j12Yv+x434HRaQWsBh4zBiTVvgxX/gcA+wuoLiMMTcUZz8R+QhY4eFyKkI80LLQ/RbubY5ijIl3fz8jIkuxhqPW21tVuUsUkWbGmAQRaQacsbug8maMSbx42wm/gyISiBXs84wxS9ybfepz9Jkj99/i/g990R3A3svt60O2AWEi0kZEqgCjgOU211SuRKSmiNS+eBsYjDM+u19aDoxz3x4HLLOxFo9w0u+giAgwA9hvjHmz0EM+9Tk64gpVEfkU65+DBjgOPFBobMxnudvJ3gb8gZnGmJdtLqlciUgosNR9NwCY7+vvUUQWAAOxpohNBJ4DvgQ+A1phTW090hjjsyckL/MeB+KQ30ERuQrYAOwBCtyb/4o17u4zn6Mjwl0ppdTPOWJYRiml1M9puCullANpuCullANpuCullANpuCullANpuCullANpuCullAP9PzRq9FxZwXkMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ddf['days'], ddf['confirmed'])\n",
    "plt.plot(ddf['days'] - 5, modelSums[ddf['days']] * population )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[5e-05,\n",
       " 5.543238127036979e-05,\n",
       " 6.192841716842426e-05,\n",
       " 6.936424363812374e-05,\n",
       " 7.768964225547436e-05,\n",
       " 8.690200264991526e-05,\n",
       " 9.703053412724029e-05,\n",
       " 0.00010812679224723995,\n",
       " 0.00012025910851309146,\n",
       " 0.000133509441855809,\n",
       " 0.00014797174280281128,\n",
       " 0.00016375127504527492,\n",
       " 0.0001809645544544315,\n",
       " 0.00019973969976459313,\n",
       " 0.0002202170712016502,\n",
       " 0.00024255012410190428,\n",
       " 0.0002669064362341259,\n",
       " 0.0002934688861972096,\n",
       " 0.00032243697532357684,\n",
       " 0.00035402828982530714,\n",
       " 0.000388480107201827,\n",
       " 0.0004260511560935002,\n",
       " 0.000467023538400154,\n",
       " 0.0005117048254381321,\n",
       " 0.0005604303404759352,\n",
       " 0.0006135656568721837,\n",
       " 0.0006715093020862655,\n",
       " 0.0007346957133569787]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "modelSums.to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelSums = pandas.DataFrame(modelOutput[:,[1,2,3]]).sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1     0.000055\n",
       "2     0.000062\n",
       "10    0.000148\n",
       "dtype: float64"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "modelSums[[1,2,10]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DatetimeIndex(['2020-03-24 23:19:34', '2020-03-25 23:37:31',\n",
       "               '2020-03-26 23:33:19', '2020-03-27 23:48:35',\n",
       "               '2020-03-28 22:14:55', '2020-03-29 23:05:00',\n",
       "               '2020-03-30 23:08:00', '2020-03-31 22:52:00',\n",
       "               '2020-04-01 23:43:56', '2020-04-02 21:58:49',\n",
       "               '2020-04-03 23:25:00', '2020-04-04 22:46:37',\n",
       "               '2020-04-05 23:34:00', '2020-04-06 23:06:45',\n",
       "               '2020-04-07 23:22:00', '2020-04-08 23:04:49',\n",
       "               '2020-04-09 22:51:58', '2020-04-10 23:02:37',\n",
       "               '2020-04-11 22:54:07', '2020-04-12 22:45:33',\n",
       "               '2020-04-13 23:18:00', '2020-04-14 23:07:54',\n",
       "               '2020-04-15 23:33:31'],\n",
       "              dtype='datetime64[ns]', freq=None)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.index + datetime.timedelta(days=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Timestamp('2020-03-23 23:19:34')"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ddf.index[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DatetimeIndex(['2020-03-18 23:19:34', '2020-03-19 23:19:34',\n",
       "               '2020-03-20 23:19:34', '2020-03-21 23:19:34',\n",
       "               '2020-03-22 23:19:34', '2020-03-23 23:19:34',\n",
       "               '2020-03-24 23:19:34', '2020-03-25 23:19:34',\n",
       "               '2020-03-26 23:19:34', '2020-03-27 23:19:34',\n",
       "               '2020-03-28 23:19:34', '2020-03-29 23:19:34',\n",
       "               '2020-03-30 23:19:34', '2020-03-31 23:19:34',\n",
       "               '2020-04-01 23:19:34', '2020-04-02 23:19:34',\n",
       "               '2020-04-03 23:19:34', '2020-04-04 23:19:34',\n",
       "               '2020-04-05 23:19:34', '2020-04-06 23:19:34',\n",
       "               '2020-04-07 23:19:34', '2020-04-08 23:19:34',\n",
       "               '2020-04-09 23:19:34', '2020-04-10 23:19:34',\n",
       "               '2020-04-11 23:19:34', '2020-04-12 23:19:34',\n",
       "               '2020-04-13 23:19:34', '2020-04-14 23:19:34'],\n",
       "              dtype='datetime64[ns]', freq=None)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pandas.to_timedelta(modelSums.index - 5, unit=\"D\") + ddf.index[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}