{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "be083add",
   "metadata": {},
   "source": [
    "## Epidemiological SEIR model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43460fe1",
   "metadata": {},
   "source": [
    "In compartmental modeling in epidemiology, SEIR (Susceptible, Exposed, Infectious, Recovered) is a simplified set of equations to model how an infectious disease spreads through a population. \n",
    "See for example [the Wikipedia article](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology) for more information.\n",
    "\n",
    "The form we consider here, the model consists of a system of four non-linear differential equations:\n",
    "\n",
    "\\begin{align*}\n",
    "\\tfrac{\\mathrm{d}S}{\\mathrm{d}t} &= - \\beta IS \\tag{% Susceptible} \\\\\n",
    "\\tfrac{\\mathrm{d}E}{\\mathrm{d}t} &= \\beta IS - \\alpha E  \\tag{% Exposed} \\\\\n",
    "\\tfrac{\\mathrm{d}I}{\\mathrm{d}t} &= -\\gamma I + \\alpha E \\tag{% Infectious} \\\\\n",
    "\\tfrac{\\mathrm{d}R}{\\mathrm{d}t} &= \\gamma I \\tag{% Recovered}\n",
    "\\end{align*}\n",
    "\n",
    "where $S(t)$, $E(t)$, $I(t)$ and $R(t)$ are stochastic processes varying in time.\n",
    "The model has three model parameters: $\\alpha$, $\\beta$ and $\\gamma$, which determine how fast the disease spreads through the population and are different for every infectious disease, so they have to be estimated.\n",
    "\n",
    "We can implement the relationship of these ordinary equations in terms of Python code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cc163d18",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.139337Z",
     "iopub.status.busy": "2021-05-18T10:57:22.138984Z",
     "iopub.status.idle": "2021-05-18T10:57:22.146777Z",
     "shell.execute_reply": "2021-05-18T10:57:22.146425Z"
    }
   },
   "outputs": [],
   "source": [
    "def ode_seir(variables, coordinates, parameters):\n",
    "    var_s, var_e, var_i, var_r = variables\n",
    "    alpha, beta, gamma = parameters\n",
    "    \n",
    "    delta_s = -beta*var_i*var_s\n",
    "    delta_e = beta*var_i*var_s-alpha*var_e\n",
    "    delta_i = -gamma*var_i+alpha*var_e\n",
    "    delta_r = gamma*var_i\n",
    "    \n",
    "    return delta_s, delta_e, delta_i, delta_r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e01fae2",
   "metadata": {},
   "source": [
    "### Initial condition\n",
    "\n",
    "The initial condition $(S(0), E(0), I(0), R(0)) = (1-\\delta, \\delta, 0, 0)$ for some small $\\delta$. Note that the state $(1,0,0,0)$ implies that nobody has been exposed, so we must assume $\\delta>0$ to let the model to actually model spread of the disease. Or in terms of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "10c2e69e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.148836Z",
     "iopub.status.busy": "2021-05-18T10:57:22.148512Z",
     "iopub.status.idle": "2021-05-18T10:57:22.155762Z",
     "shell.execute_reply": "2021-05-18T10:57:22.155410Z"
    }
   },
   "outputs": [],
   "source": [
    "def initial_condition(delta):\n",
    "    return 1-delta, delta, 0, 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394a8221",
   "metadata": {},
   "source": [
    "### Model parameters\n",
    "\n",
    "The model parameters $\\alpha$, $\\beta$ and $\\gamma$ are assumed to have a value, but are in all practical applications unknown. Because of this, it make more sense to assume that the parameters are inherently uncertain and can only be described through a probability distribution. For this example, we will assume that all parameters are uniformly distributed with \n",
    "\n",
    "\\begin{align*}\n",
    "\\alpha &\\sim \\mathcal{U}(0.15, 0.25) & \\beta &\\sim \\mathcal{U}(0.95, 1.05) & \\gamma &\\sim \\mathcal{U}(0.45, 0.55)\n",
    "\\end{align*}\n",
    "\n",
    "Or using `chaospy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6e4aabe7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.158026Z",
     "iopub.status.busy": "2021-05-18T10:57:22.157708Z",
     "iopub.status.idle": "2021-05-18T10:57:22.165495Z",
     "shell.execute_reply": "2021-05-18T10:57:22.165142Z"
    }
   },
   "outputs": [],
   "source": [
    "import chaospy\n",
    "\n",
    "alpha = chaospy.Uniform(0.15, 0.25)\n",
    "beta = chaospy.Uniform(0.95, 1.05)\n",
    "gamma = chaospy.Uniform(0.45, 0.55)\n",
    "distribution = chaospy.J(alpha, beta, gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2015813",
   "metadata": {},
   "source": [
    "### Deterministic model\n",
    "\n",
    "To have a base line of how this model works, we will first assume the uncertain parameters have some fixed value.\n",
    "For example the expected value of the uncertain parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "800f8c00",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.167462Z",
     "iopub.status.busy": "2021-05-18T10:57:22.167148Z",
     "iopub.status.idle": "2021-05-18T10:57:22.202134Z",
     "shell.execute_reply": "2021-05-18T10:57:22.201797Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.2, 1. , 0.5])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parameters = chaospy.E(distribution)\n",
    "parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397797f9",
   "metadata": {},
   "source": [
    "We then solve the SEIR model on the time interval $[0, 200]$ using $1000$ steps using `scipy.integrate`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b4f20637",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.205180Z",
     "iopub.status.busy": "2021-05-18T10:57:22.204840Z",
     "iopub.status.idle": "2021-05-18T10:57:22.215243Z",
     "shell.execute_reply": "2021-05-18T10:57:22.214919Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "time_span = numpy.linspace(0, 200, 1000)\n",
    "responses = odeint(ode_seir, initial_condition(delta=1e-4), time_span, args=(parameters,))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcb665a4",
   "metadata": {},
   "source": [
    "We then use `matplotlib` to plot the four processes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2faa948e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.218544Z",
     "iopub.status.busy": "2021-05-18T10:57:22.218198Z",
     "iopub.status.idle": "2021-05-18T10:57:22.351116Z",
     "shell.execute_reply": "2021-05-18T10:57:22.350772Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f477d8c0130>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAERCAYAAACXT3dwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABLaklEQVR4nO3dd5xU1dnA8d+dmZ3Z3thll945dERARJSuYifGGEtsSYzdWBKNJpbYYuU1mtgSE1usMSpqLDQRCb2ItENdYIFddtnep9z3jzu7DLDLzi4zO1ue74f9zNw6D3dnn3vuueeeY5imiRBCiI7BFukAhBBCtBxJ+kII0YFI0hdCiA5Ekr4QQnQgkvSFEKIDkaQvhBAdiCR90aYYhtHLMIxPDMP41jCMeYZhfGMYxk3+ZcMNw1hqGIZpGMYy/7Lan83+da4wDGOzYRjVhmF845/3jGEYOYZh5PrXXWYYxlbDMK6JwP/v9/5YHgxi3fuDXVeIWoa00xdtiWEYC4APTNN8wT89GXjeNM3h/unewE5ggGma2wK2+8Y0zcn+91cDj5im2T1g+WuAwzTNn/mnLwPeAAYF7qcl+GPJMk3zwVCuKwRISV+0PScB39ROmKb5DfCvILa7q4mfMxuwAyc0cTshWjVJ+qKt2QXcZRhGXO0M0zQfb2hlwzB6G4bxmmmay5v4OQ7/a3Y9+3T6q4FMwzBuMgzjS3910ETDMO40DGOxYRhLDMNID9hmrGEYC/3VUgsNwxgbsKyfYRiL/Nv8E4g54vNGB2w3zzCMQU38vwhRR5K+aGtuAc4D9hqG8U/DMCY1sN6//HX27zbzcy4CHjdNc+mRC0zTrKmtKgJcpmnOAP6KdcWx2DTNCUAe8AsAwzCSgC+AB0zTnAj8AfjCMIxk/z7eBv5rmuZ4/7LTaz/Lv+2XwIOmaU4CZgGfGIYhf7uiWeSLI9oU0zTnAT2B3wC9gQWGYbxcz6qX+xPzJU3Y/en+EnwW8HvgtSC2meN/XQ/EBpwk1gF9/e/PBUr8VVGYprkIKATONwyjF1aV1Vv+ZXuB7wL2fy5QZprmfP/yz4FMYFwT/l9C1JGkL9oc0zTLTdP8u2maU4ApwC8Nw+jbwLpZpmleHeSu5/hPFCcANcBjQWxT6n/1BLyvnXb633fHKvkHyvPP7+Kfzg9YVhDwvjuQGtgSyb9tpyBiE+IokvRFm2IYxouB06ZpLgQOAkmNbNdQNdBRTNMswirpX+BvDXS89gDpR8xLx7pfsD9gulZgQt8DZJumObn2BzgR+DoEcYkOSJK+aGumG4ZxUu2EP5n7gM2NbPfHJn7Of4BtwM1N3K4+nwEJhmFMBDAMYwKQAsw2TXMXsBy4wr+sGzDpiG3Tam/8+m9gL6CRk5wQDXE0vooQrcqTwFOGYfiwmlT6gJmmaVb6W7XM8q/3vGEY5QHbpYP1cBbwOyC9tu2+YRjPADOsxcYbpmleaZqmzzCMJ/z7GWSa5rmBQRiGUVvSftcwjKuAZ4FM//2Fz4GrgWjDMO4wTXOWYRgzgGf8N2BN4Cz/FQXAZcDrhmGch/WMwVzgasMw9pmm+YphGGf7tzUAA+uGcJ5hGPf7464yDGOPaZqvHu/BFe2fPJwlhBAdiFTvCCFEByJJXwghOhBJ+kII0YFI0hdCiA5Ekr4QQnQgrb7JZl5eabObFyUlxVBcXBnKcEJC4moaiatpJK6maa9xpacnGPXNb9clfatZc+sjcTWNxNU0ElfTdLS42nXSF0IIcThJ+kII0YFI0hdCiA5Ekr4QQnQgkvSFEKIDCUuTTaVUJvAIMFJrPbae5TasASpKsUY/elVrfdSwdEIIIUIrXO30TwU+wRqBqD4XA4la698ppVKBpUqpwVprbyg+vNrj48tNucTFuaisqMFmGBgG1g8GNn9LqEPzDWxQ997wL8MAmwE2DOw2gyi7gcNm4LDbrFebQVTte7tBlM2Go3Ydm9Fqm4IJITqusCR9rfW/lVKTj7HKOfhH/tFaFyilqoChWOOKHreKGg+PfL01FLtqNpsBMVF2Yp126zXKTkyUjRinnaRYF04DkmKiSImNIiUmimT/a2psFOnxLuw2OWEIIUIvUk/kdubw8URL/POOkpQU0+QSsz3ayY9P7AaAz2dimuAzTUysV0zwmWBiWq8mmKYZsA5gmofW8YHH58PtNete3V4fniOmA5d7fSblNV7Ka5p+8eKwGWQmRdM9OYZuKTH0TIllYEY8AzMS6J4cgy1MJwS73UZycmxY9n08JK6mkbiapqPFFamkfwBICJhO9M87SnMfQ/7dlH4kJ8dSVFTRrO2Pl8fro9Lto8LtpdL/U1FjvRpRDnILyimq9FBY6aao0k1RhZvCSjf55TUcLK8hu7CS7MJKaxylALFRdvqlxTGyWyKjuicxsmsiSTFRIYk5ksfrWCSuppG4mqa9xpWenlDv/BZL+kqpOCBWa52HNZzcROBNf51+NLChpWJpCQ67jQS7jYToow9xY7/Mao+PnJIq9pdUsb+kml0FlWzPL2drfjkHy2v4YX8JP+wv4a2V2QAMTI9jUv9OTOqXxsDOcXIvQQjRoHC13pmENdBzF6XUH4BnsMYMHQ5cD7wPjFJKPQD0BK4M1U3c9sDlsNErNZZeqUdf2hVW1LD5QBlrs4tZs7eEDftL2JJXzpa8cv62ZDddk6I5d2gG5w3NIDMxOgLRCyFas1Y/Ru7x9LLZXi/bAlV7fKzaU8TCbQdZuP0gB8trAGv07Al9U7n6pB6M7JbU4nGFksTVNBJX07TXuBrqZbPVd60sjs3lsHFKn1RO6ZPK3dP7s2JXEZ+sz+Gbbfl8t6OA73YUMKp7EtdP6MWJ3ZMjHa4QIsIk6bcjNsNgXO8UxvVOobCihndX7+X9tftYk13Mde+tY9rANG6Z2IduSTGRDlUIESHSDUM7lRLr5IZT+/DpteP41fheuBw25m3J56evreK91XutpqtCiA5Hkn47F+9ycO0pvfjw52M5Q6VT7fHx9ILt3PTBOnJLqyMdnhCihUnS7yAyElw8eu5gnjx/CCkxUazcU8yVb61mTXZxpEMTQrQgSfodzJQBabx79WjG9kymoMLNDR+s46N1+yMdlhCihUjS74BSY5089+PhXD66O16fyWNztvLPZbtp7c13hRDHT5J+B+WwGdw2uS93T+uPAbzwXRZPz9kiiV+Idk6Sfgd30QldefTcwdhtBq8s2sk/lu2OdEhCiDCSpC84XaXz8NmDsBnw0uJdvLd6b6RDEkKEiSR9AViJ/5ELhgHwzILtLN5REOGIhBDhIElf1PnJ6O78anwvTOD3n29i58HW1x+JEOL4SNIXh/nF+J5MG5hGeY2X33yygfIaT6RDEkKEkCR9cRibYfDADMWA9Dh2F1by5LxtkQ5JCBFCkvTFUWKi7Dx2zmBcDhv/3XiALzblRjokIUSISNIX9erdKZY7p/QD4Im529hfUhXhiIQQoSBJXzRo5vBMJvfvRHmNlyfmbpMHt4RoByTpiwYZhsHd0/oT77KzeGcBc3RepEMSQhwnSfrimNLiXdw6sS8AT8/fTnGlO8IRCSGOhyR90agLhmcyqnsShZVu/rZkV6TDEUIcB0n6olE2w+C3U/thM+Dfa/fJQ1tCtGGS9EVQBqTHc8HwTLwm/HnhjkiHI4RoJkn6ImjXT+hNnNO6qbs0S/rmEaItkqQvgpYa6+SacT0BeHHxLmnCKUQbJElfNMnFo7qSGhvFxpxSFklPnEK0OZL0RZPERNm56qQeALy8OAuflPaFaFMk6Ysmu3BEF9LjnWzJK+ebrfmRDkcI0QSS9EWTRUfZufokq27/teV7pG5fiDZEkr5olvOHZZASE8Wm3DJW7imKdDhCiCBJ0hfNEh1l5+JRXQF4c0V2hKMRQgRLkr5ototO6Eq0w8aSrEK25pVFOhwhRBAc4dqxUmo6cCFwADC11n88Ynkf4GlgBXAC8LbWena44hGhlxwTxQXDM3lvzT7eXJHNQ2cPinRIQohGhKWkr5SKBV4CbtdaPwiMUEpNO2K1u4DvtNaPA08Az4QjFhFel43ujt2Ar3UeeWXVkQ5HCNGIcJX0xwO7tNa1WWAxcA4wL2CdXCDd/z4dWBWmWEQYdU2KZlL/NOZvzefjdTlce0qvSIckWoDPNPH5TLym9d7rM/3zwGuaAfOs5aYJRV6T4uJK6tp6mdS9N7HWOXyZNc8/WTevdvrw9c1D+zJrlwfXqiyhpJrS0uaNDNfUhmtNWT2+qIrOLjtJMVFN+5BGhCvpdwZKA6ZL/PMCzQI+UkrNAk4CHq5vR0lJMRiG0awg7HYbycmxzdo2nNpbXNec2sdK+utzuO1MRZQ9tBeQ7e14hYtpmpRVe9lbXEV+WQ3FlW7Kq72U13gor/ZQUeOtey2r8VBZ46Xa46PG48Pt9VHj9VHjManx1L4/NN/jNa1k7rNepZVuy3j1ytFM7JIU0n2GK+kfABICphP98wK9Bvxda/2OUiod2KqU6qu1PuzZ/uLiymYHkZwcS1FR6+sGuL3FpVKi6ZMay86CCmav2sO0gemNb9QCcYVbS8ZV5baSeW5pNbml1eT4Xw+UVlNY4aao0vrx+FouG9sMsNsMbIaB3TCw2bBeDQObzcBuUPfeZoDNZsP0+Q4rxBlA7aSBgf+fNW1Y8wLLfNb6xlHb4p9/+La1ax2bw2HD4/HVbddUTd0k2M9wOOzYPN5mf8fS0xPqnR+upL8E6KWUcvmreCYALyilUgGP1roE6AHs969fCPiQ1kRtkmEYXHRCF56av50P1u4LedLvSDxeH9vzK9iYW0pWQYX1c7CC/SXVQVUNxEbZSY1zkuCykxjtINbpINZpJzbKfvir005MlB2n3SDKbsNptxFlN3A6bIdNW+8NHDbbYUneZtDkK3A5eTdNuOIKS9LXWlcopW4AnlNK5QHrtNbzlFJPAgXA48DtwG1KqVOAPsC9Wmt5pr+NOntIBn9dlMWqPcXsOFhO305xkQ6pTSipcrNyTzGr9xSxMaeULXnlVPtLnYHsNoNuSdFkJrjICPxJdNEp1klyTBRJMVG4HLZWm8RE6xC2Jpta6znAnCPm3RXw/jvgu3B9vmhZ8S4HZw3pzIff7+ejdTncOaVfpENqlUzTZPOBMhZszWdpViGbc8uOKsH3TIlhcEY8/dPi6J0aS+/UWLonR+MI8b0S0TGFLemLjmfm8Ew+/H4/X2zM5daJfUJ+Q7ct25ZXzucbc5m/JY99JYeatkbZDUZ0TWRMj2SGd01kcEY8idGhba0hRCBJ+iJkVOd4BqTHsTWvnEXbDzK1g9ftV3t8fL35AB+ty+GH/SV18zvFOZk6II1T+6YyqnsSMVH2CEYpOhpJ+iJkDMPgvGGZzFqwndnrczts0i+v8fCf7/fz1spsCircAMQ57cwY3JkZgzozolsitmY2QxbieEnSFyF11qDOPLdwB0uyCjhQWk3nBFekQ2ox1R4f767eyxsr9lBS5QGsq5+LR3XldJUuJXrRKkjSFyGVHBvFxH6dmL81n/9uzOVq/5i67Zlpmszfms9zC3fU1def0C2Ra8b1ZHzvlGY/XChEOEjSFyF3/rBM5m/N59MNuVx1Uo92nfTyyqp5bM5WvvOPF9wvLZbbJvXl5N6pEY5MiPpJ0hchN653CunxTnYXVrJuXwkju4X2MfLWwDRNPlqzl4c/30RptYd4l52bT+vDBcO74LC135OcaPsk6YuQc9gMzhzUmbdWZvPlpgPtLulXub08MW8bn23IBeDUvqncM31Ah7p/IdouaUgtwmLGIKt/vblb8vF4j37CtK3KLqrk5++s5bMNuURH2bjvjIHMmjlUEr5oM6SkL8JiYOe4uk7Ylu0qYkLftl/HvSGnlDs+Wk9BhZueKTG8cNmJZERLixzRtkhJX4SFYRicOdhqp//l5iM7WG17vttxkOvf+56CCjfjeiXz+uWjUJn192IoRGsmSV+EzZn+Kp6F2/KpdHsjHE3zLdiaz28+2UiVx8c5QzN49kfDiHfJRbJomyTpi7DpnhzD8C4JVLp9fLvtYKTDaZaF2/K557NNeH0mPxvTnQfOHCgdn4k2Tb69IqxmDLZK+22xiud/Owv43aeHEv6tE/u062cORMcgSV+E1bSB6dgNWJJVSJG/H5q2YHNuKb/7dCMen8llo7tJwhfthiR9EVad4pyM7ZWC12eyYFvbGCMnp6SK2z7aQKXbx4zBnbltUl9J+KLdaHLSV0p1Ckcgov063d/b5lydF+FIGlfp9nL7Rxs4WF7D6B5J3HfGQEn4ol1ptAmCUioeOJ1DA52fB/wknEGJ9mVS/048Ntdg1Z4iiircJMe2zkFCTNPkT3O2si2/nJ4pMTx5/hCcDrkYFu1LMN/oz4BpWOPY9gHa/lM2okUlxUQxtmcyXpNWXcXzn3X7+WLTAaIdNp44f4iMYCXapWAaG+/UWt9cO6GU6hPGeEQ7NX1gGkuzCpm3JY8fjegS6XCOsjm3lGcWbAfg3jMG0D9NBnYX7VMwJf0spdTpSqleSqmewFXhDkq0P5P6p2E3YOXuolbXiqfK7eX+/2rcXpMfj+zCWYMzIh2SEGETTNK/HrgXeA14HbginAGJ9ik5JspqxdMKq3heXJzFzoIKeqXEcNukvpEOR4iwCibp36O1nlL7A1wX7qBE+zR9YBoA87a0nlY8q/YU8c6qvdgN+OPZg4iWIQ1FO9donb7W+jWl1HRgJLBWaz03/GGJ9mhS/zT+NGdrXRVPpFvxVLm9PPL1Fkzg5yf3ZKh0oCY6gEZL+kqp+4A7gF7Ab/zTQjRZckwUY3taVTzftIIqnteW7yG7qIp+abH8vAOM5SsEBFe949Ran621vlVrfRYQG+6gRPs1ra6KJ7JJP+tgBa8v3wPAPdMHSCdqosMI5pt+5LBH7WcYJNHiJg+wWvGs2F1IUWVkWvGYpskT87bi8ZlcMCyz3Q3nKMSxBNNO36OUmg3sAPoBy8IbkmjPaqt4lu4q5Jut+cyMQJv9+VvzWbmnmKRoBzdPlMdORMfSaElfa/0w8BcgG3hOa/1I2KMS7Vokq3jcXh/Pf7sTgOsn9CY5Rp66FR1LUMP/aK2/Br4GUEr9VGv9XlijEu3a5P5pPD53a10VT0sm3g/W7mNvcRV9UmMjcpUhGub1eigszMPjqWnRz83NNTBNs0U/MxhNiSsxsROxsfFBrdtg0ldKva21vkwptROo/WQDSAQk6YtmS46NYkzPZJbtKuLbbQc5f3hmi3xuSZWbV5fuBuDWSX1w2KT3zNaksDCP6OhY4uIyW7RnU7vdhtfb+m5VBhtXTU01RUX5QSf9Y1XvPOh/naW17uv/6QP8Iag9C3EM02q7W27BB7VeW7aHkioPY3omM6GP9BvY2ng8NcTFJUpX1k0UFeXE5/MEvX6DJX2t9Rb/27obt0qpEUB1MDv2P9B1IXAAMLXWfzxiuQHc4p/sDSRrrX8edOSiTZvcvxNPzN3K8t1FFFe6SQpzFc/B8hreX7sPQEbBasXk99J0TT1mwTTZnFH7Rmu9DlCNbaCUigVeAm7XWj8IjFBKTTtitZ8BRVrr57TWdwDPBhu0aPtSYp2c2CMZr89k4fbwD5r+xoo9VHt8TOrXicEZ8uStaD2ysnZyzjmns3v3rsPmP/TQAzz55GMh/7xj1elfBVwN9FJKTfbPNoCqIPY7Htilta69KlgMnAPMC1jncuBLpdStQCbw9yZFLtq86QPTWLG7iPlb8jl/WPjq9fPLa/jw+/0AXDu+V9g+R4TG2Ge+Ddm+Vtw5sd755eXlPP74Iyg1CLe7hmXLlvLSS6+G7HMbMnPm2Xz88X8BWLJkMWVlpZx++gz69u131LpTp05j0aLQHYtax2q98zHwDfAr4BX/PC+wP4j9dgZKA6ZL/PMC9QIStdYPKaUGYp0ABmutvYErJSXFNPuSz263kZzc+h4glrgsF4zuwRPztrF8dyE2VxSJDVTxHG9cLyzZRbXHx+mDOzNOHfk1bD75PTZNY3Hl5hrYQ/xkdEP727t3N/v2ZfPggw/jdDoZP34CTzzxCHa7neuuu5GbbrqOG264mSFDhvKnPz3C8OEjWL/+B+666x5+/eub6N69BwMHKlasWM5TT/0f+/bt5cMP36dfvwGsWrWCBx54mPXr1zF37tf07t2XlSuXccMNt5Cbm8vbb7/B2LHjmD9/Li6XkxkzzsYwDD777BOczig2btzIrFnPAWCzWcfkH//4G3a7nZqaGnw+H9ddd+Nh/x/DMIL+nR+rTr8YKAZ+HzhfKXUisLqR/R7g0PCKYLX4OXDEOiX47xdorbcopRKBHkBW4ErFxZWNfFTDkpNjKSqqaPb24SJxWezAiT2SWbm7iNmr93Du0PpL+8cTV35ZNW/7u1u4akz3kP7/5PfYNI3FZZomXq+vwdJ5czTU+mXgwMHMnHkRt912M5WVFUyaNJVTT53E4sWLSEhIYtiwEfh8JiUlpWzZohk3bjw33fRrUlI6MWzYCE44YRRnnHEWlZVVvPfeOyxduphBg4bg8XiIjo5h48YN3Hffvbzzzr9JTk5hzJiT6NdvABkZGVx22ZUATJ48lcWLF+H1+jBNkylTpjN48BDuvfe3zJs3h9jY2LoYXnjheW6//bc4HFEsXfq/o/5fpmkedWzT0+uvxgxmjNzuwK1AGlb1znBgTCObLcGqFnL5q3gmAC8opVIBj9a6BKuqp6//MxKxckBOY/GI9mX6wDRW7i5i3pb8BpP+8Xhn9V6qPT4m9++E6hxckzbR/hUUHGTYsOGcd94FgMl5583gnnvuw+u1WsFUVlqFzZSUFP7+99dZtmwp11zzMz766PPD9hPYjv6UU05j3LiT2bRpIykpKfWuV1tr4fE03Nqmvrb5DoeDSy65HICMjOP7Ownm4azHgQ+xbuj+G7i4sQ201hVKqRuA55RSecA6rfU8pdSTQIF/n08ATyql7sXq3uEqrXUw9wtEOzK5fxpPztvG0qxCSqs8JEQH9bxgUMqqPXV1+Vef1CNk+xVtX1VVFX/5y7OceOIYqqoq6devPyNGjOTVV1/hgw/eZefOHSxYMI/k5BQ+/PB9Bg0azJgxY4mJiQFg2bKlZGdns27dWp5++llOPXUir7zyIllZO9i9exe33/5bHnzwUV588Xm6dOlGSUkJAwYMZOTIUbz44l+IiYkhO3sPWm9i1aoV7NixnU8//ZiFC+dTUVHB5MnTePzxR9i8eSOVlZVce+0N/N//PUWnTmkkJBxfQwSjsSe+lFK/0Vo/rZS6W2v9hFLqPn/XDC0iL6+02Y/KtdXL3EiJVFzXv/89q/YU8+AMxTlDjx6qsLlxvbUymz8v3MGo7km88tORoQg1JHGFW1uNKydnF5mZLX+jvakPZz355GNMmHAaEyacFsaomhZXfccuPT2h3puhwdw1Ga2U6gWkK6V+BkwJKgohglT7oFYoR9TyeH28syobgCvGdA/ZfkXHlp+fzw8/fM+CBfMaX7mVCuZa+lmsm7IvAk8Bz4czINHxTBmQxlPztrF0VyFl1R7iXcdfxfO1zuNAWQ19UmOZ0FeevhWhkZaWxptvtu1eaIIZLjGwK+ULlVInhzEe0QGlxTkZ1T2J1dnFfLv9IGcPObqKpylM0+StlVYp/2djumOTpzyFqHOsh7P+0cCiETTeekeIJpk2MJ3V2cXM25J/3El/5Z4ituaV0ynOyYzBoWuXL0R7cKySvg94s575V4QpFtGBTR3Qiafnb2NpVsFxV/F8sNZqsfPjEV1wOmQYRCECHesv6zatddmRM5VSO8IYj+ig0uJdnNAtkTV7S1i04yBnDW5eaT+npIpvt+Vjtxn8aETLdNksRFPNnz+X556bVdclQ0s61hO5ZQBKqSuPWHQe8JNwBiU6pmkD01mzt4T5W/KbnfQ/WrcfrwnTB6SRFu8KcYSiJbRE3ztr1qzixhuv5dJLryAtrRP5+Qfp3Llz3QNQ4TZ16nSee25Wi3zWkYK5hr4Gqw8egJ5A8B03C9EEUwem8cyC7fxvZwHlNR7inE2r4qnx+Pj4B+uh7p+MklGxRMNGjRpNRkYmM2deSJ8+fXjqqceJjY3lhBMG8+mnX3H33Xdwxx13sWHDev71r9e54oprWLNmFTNn/phTT53IQw/dz8CBip07d3D66TPo378/TzzxKMOGjWDDhh94/PFn+Oc//47dbsfttvrLufbaG7j//ntJT+9Mp06dIvZ/D+av6lda6621E0qpG4+1shDNlR7vYmS3RNbuLeG77QWc2cSbsPO35lNQ4aZfWiyjuiWFKUoRbqHse6cxb7/9Jp06dSIxMZHzz/8RFRUV3HLL9TzxxCyUGkTnzhksXLiAK664mnPPPZ9LL/0xd931e9zuGi677AqKigq59NKLeOWVf7J16xZOPvkUbr75NsrLy3nppb9w++2/xel0snTp/9B6M1lZO3jooccoKyvl/fffabH/Z6Bgkn61Uqqn/30iMBl4IWwRiQ5t2sB01u4tYe6WvCYn/Q/8g6RcfEJXGYxDBOWyy66gT58+1NS4yc3NpUePnrhcLrZv34ZSgw5b91i9FxzZR8+HH84+qr+c1vKdDCbpLwR2YnW2VgL8LawRiQ5tygCrimdJVmGTqni25pWxbl8JcU47M5p5P0B0HGvWrCI3N4e33nqdlJQUKisriYpysnr1Sp566lkuuOAsDMNg6NBh5OUd4J133mLFimX87nf3MXHiZBYtWsg777xFVtYOHnzwUXbt2nVYHz1xcfFH9ZczZco0evfuy3PPzSImJobc3Bw2bFjP0KHDWvT/HkzfO2dqrb9qoXiOIn3vtJzWEte1765l7d4S/niW4uwhGUHF9fT8bby3Zh8XjezC3dMHtEicreV4HamtxtUa+97ZvXsXTz31J55//qUWjiqyfe/MV0rdoJT6i1LqRqWUM6gohGims/zVOl9sPHIIhvrVeHx8ucla94Lh0kxThM6CBfPYsWM7WVk7Ix1KyAST9P+B1Yf+DqyncRt6UleIkJg2MB2HzWD57kLyy6obXf+bbfkUV3kYmB7HIBn/VoTQVVf9nM8/n0Pv3n0iHUrIBFNhekBrfWfthFLq2fCFIwQkxURxat9Uvtl2kK8259G/e8ox15+93mqmKaV8IRoXTEm/8IjpbACl1FmhD0cIy1n+/ne+2HTsKp59xVUs31WE025IPztCBCGYkv4lSqlfYo1d2wsoVEqdg/Wg1tFDuAsRAhP6pJLgcqAPlLE1t5R0l73e9T5dn4OJ1eonMbr+gdWFEIcEk/S/Ap6rZ748pCXCxuWwMW1gGh//kMPsdfv5xdijB0Lx+kw+3ZALwPnDpGpHhE5ZWRkPPHAvXbt2Y9SoE5k69fRIhxQyjVbv+Ovzi7AGRi/UWu/y/9wd7uBEx3bWEKu6Zvb3+/DV07R4+e5Cckur6ZoUzZieyS0cnWjL1qxZxfjxJ/K3v71Y7/L163+gc+cM7rzz7qAT/syZZ9e9X7JkMXPmfBmSWEOt0ZK+Uup84K9YdfvJSqmbtNafhj0y0eGd0C2JzAQX+4qrWJNdzOgeyYctn/2DVco/b2iGDJTSjrREh2u1fe+ceuokJk0az4UXXkx+fh5Op5N7772fuXO/QutNzJ8/l+7de/Cf/7xP3779Wb16Jffd9xDr1//A/Plf07t3X1auXM71199c97DX2LHjmDdvDk6nk9NPn8GsWU+SmtqJgoKDDBo0hH79+nP11ZexZMlqvv32G2bNepKPP/4vr776MmBQWFjACSeM4swzw3PbNJjqnTOAflrrGqVUNNbwiZL0RdjZDOvm7GvL9/DZhtzDkn5plYdvt+djAOfWM5i6EMGIi4tj2LARTJkyjaFDh3H++TMAmDJlGk6nk6lTp3PdddcwaNAQfD4fMTGxaL2JBx/8Pe+882+Sk1MYO3Yc/fsPICMjk5/97Kq67RcvXsS2bVtZu3Y1b7zxLj6fj4kTx/HddyvIyLCqIydOnMysWU8CsGvXLlwuF+eeez59+4bvdmkwSX+X1roGQGtdpZTaHbZohDjCecMyeW35HubqPO6c0q9ucJV5W/Ko8ZqM6ZlMZmJ0hKMUodSSHa7VSkiwnu+w2eqv8T7llNMYN+5kNm3aSErK4U2Ia3s1qO1bx+OpvyPiwL53HI4oPB4PNTU1dfNuvvk2SkuL+fvfX6ZLl67cfvtvmv8fOoZgkn4/pdQdWA9n9cNqwSNEi+iZEsO4Pqks21nA15sPcOHIrgD819+U82xppimaobbvnffee5sdO7azaNFChgwZSm5uDqtXr2TBgnlovYkNG9bzm9/cwyuvvEhW1g52797F7bf/lgcffJQXX3yeLl26UVJSwoABAxk5chQvvvgXYmJiyM7eg9ab+PnPf8WoUaN5/fV/UFhYwL33PgDA2Wefy6xZT9K//4C6Pnhmz/4Pffr0IzOzC0oNDtv/PZi+d+KBe7Gexl0LPF7fiFrhIn3vtJzWGte3u4q489/rGJwRzxs/O5F9xVVc8PfluBw2vrz+5OMaWvF4tNbj1Vbjao1970RSxPre8Sf4p4AHgWdaMuELAXDmkAwSox1syi1DHyjjq81WKX9Sv04RS/hCtFWNJn2l1BXABqw+dzbWM3yiEGHlirLXdcL28br9/Hej1Wrn7CFyA1eIpgqmG4bzgd5a6xFAH2BmWCMSoh61/ep8vjGXrIJKUmKiGNcrObJBCdEGBZP01wS23gFWACiluoUzMCECDUiPZ0hmApVuq47zjEHpOOzBfH2FEIGCqRAdppR6CKv1Tl+gs7+K5zzgJ+EMTohAPx7RhY05pcChPveFEE0TTFGpC+DFaqrpBfZjVfOkhjEuIY6SFHOojFLl8UYwEiHC77XXXuWWW64P+X6DKenfqrX+4ciZSqmhIY9GiGOYo/Pq3n+wdj+jexy7n30hGrJmzSpuvPFaLr30CpKSElm48Bsee+xJMjO7RDq0OlOnTmf58mUh32+jSb++hO+fv+FY2ymlpgMXAgcAU2v9xwbWuxx4C0iQ5qCiIeU1Hr7ZdhAAmwHfbM0np6RKnsZth1qy752ZMy+kT58+/u4S1rBy5csoNYjt27dy2mmTGTVqNI8++iAjRpzA6tUrufrqX1JSUsLcuV/Rt28/tN7Mtddezy9+cQVXXfULzjrrHG655Xruued+li1bgt1ux+2uwefzceqpk7jxxmu58MKfsHPndsaPn0B29h66dOnKgQO59OzZi/POm8mdd97KkCFDcbvdITsOgcJyJ0wpFQu8BNyutX4QGKGUmlbPeoOBIeGIQbQvC7bmU+3xMapbIqerdLwmfLB2X6TDEm3cp59+wv3334vNZsfpdLJ+/Tq8Xi/p6Z1ZtWoFH3zwLj169OSnP72Mu+/+PZ07Z/Dgg7/npptu5dJLf4bbXYPWm7nzzrvJytqBYRhMnDiZXr1689JLf8HpdBIXF8+GDesZPHgIw4YNZ+rU6TzzzHOMGjWaL774DIC0tHRWrVrJggVziYuL41e/upHp08PTnXODJX2l1KvAr4Aof6udphiP1WdP7QCni4FzgHkB+48F7gKuw3riV4gG/dc/SPpZQzIYkB7HV5vz+PD7/Vwzrqc8oNXOtGTfO+eddwE9evTgvPNmMHjwELp27cYll1yOz+dj0aKFhw2IbhgGPl/9T8hOn34mzz77DPHxiVx++RUAOBwOLrnkcoC6DtbA6uenth+exMTkunUWLJiHzxf+e1XH+mvZqLX2KqV+DzxUO1MpdaXW+o1G9tsZKA2YLvHPC/Qo8JC/984Gd5SUFHNYR0VNYbfbSE6Obda24SRxNU1eWQ0r9xQRZTe4cGxPkmKi6vrj+WxzHtdPiswAbq31eLXVuHJzDewt1Ax39epV5Obm8tlnn3DLLbfxhz88wJ//PIvq6irefPOfFBcXM23a6Zx88ngefvgB3n33X2zatIE77vgtDz/8GC+99Bd69+5LdHQ006efjmEYXHzxJWzatJH09HQArrvuRp599mnS0tJISEggJ2cfO3fuYPbsj/j1r+9AKcX48eP5299exGazMWCAYsqUqXz22WxeeeUFiouL2LlzB3v27Gp0YHbDMIL+nTfY945S6l2gikN97gAYwHCt9Zhj7dRflXOv1nqaf/oOoLvW+g7/dA/gYWCzf5M/AQ8A/9Varwzcl/S903Jaa1wfbsjl8S81Uwak8eT5Vm3gsqxCbv7wB1Jjo/jklycRHVX/cIrh1FqPV1uNS/reOVy4+t45Vkn/GmAUcC3wesD8K4KIYQnQSynl8lfxTABeUEqlAh6t9R7g6tqVlVJ/AmbJjVxRn9nf7wc4bODzk3olMzgjnk25Zcxen8vFo7pGKjwh2pQGr6W01pVa6/9hNdlcWPsD3NbYTrXWFcANwHNKqUeAdVrrecDvCBhbVymVrpT6g3/yLnnKVxxp58EKNu4vId5lZ0KfQ4+GGIbB1Sf1AOCtlXvwtMKSmhCtUTB3wLoqpf7BoWqeXwK6sY201nOAOUfMu+uI6TzgEf+PEEf5cpPVudq0Aem4HIeXUSYPSKN3agxZBZV8sj6HH4+U0n5bZ5pms+/hdVSNdY9/pGDumtwF3I71RO5vsUrrQoSdaZp8udl6IGtGPd0u2AyD6yf0BuDvS3ZT5ZandNsyh8NJeXlJk5NYR+d212CzBd+CLZg1N2utl/vfL1VKTWlWZEI00bp9JewrriIj0cWJPZLqXWfqgLS6uv331+zjSn+Vj2h7UlLSKSzMo6ysqEU/1zCMVnmiCT4ug8TE4HvFCSbpD1RKjQJ2Yg2X2D/ovQtxHL70D4l47vAu2Bq45DcMgxtP7c0tH67n9RV7+NGILiRES7v9tshud5CW1vLdILTV1k7NFUz1zjPAX4Fs4M/AkyGPQogjeLw+5m7JB+D8Rurqx/VKYXSPJEqqPLy2fHdLhCdEmxVM3zubgVNaIBYh6izdVUhRpZs+nWIZnJlAcXFlg+sahsEtE/tyzb/W8PaqvZw3LJPeqa3v4SQhWgMZhUK0Sl/UdrswuHNQrTmGZiZw/vBMPD6Tp+dva5V1tEK0BpL0RatTXuNh4XarR80zBwU/WMpNp/YmweVg2a6iuh45hRCHk6QvWp2F2w5S7fExsmsiXZOC7zo5JdZZ14TzmQXbKav2hClCIdquoJK+UmqIUupzpdTXSqkfhzso0bF94W+1U1/b/MZcOLILgzPiyS2t5rlvd4Q6NCHavAaTvlKqU8DkpcCPtdZnAMPDHpXosA6W17B8VyF2m8F0ld7k7R02g/tnKBw2g4/W5bA0qyAMUQrRdh2rpP+yUupS//sa4GKl1NlYY+YKERZzdB4+E07pnUJyTFSz9tE/LY5fnWL1OPjI11sprgzPCERCtEXH6nDtIiBRKfUv4G3/ugOA+1soNtEBfXkcVTuBrhjbg6GZCeSWVvPQV1ukNY8Qfses09dav4zV384fgVit9Z+11rktEpnocHYVVLAhp5TYKDsT+3VqfINjcNgMHj13EAkuB99uP8jbq/aGKEoh2rZj1elfpJSagzXW7aNAvlLqQ3WsYa6EOA6fb/T3qDkwLSSDonRLiuH+MwcC8PyinazaU3Tc+xSirTtWSf8ErfXpwMXAZVrr97EGVLmjRSITHYrXZ/L5BivpnzM0I2T7nTwgjZ+N6Y7XZ3LX7I1kFbS+PlaEaEnHSvrp/hu3FwFuAK11gdb6uhaJTHQoK/cUcaCshq5J0YzqXn+Pms1182l9mNivEyVVHm7/aD1FFXJjV3Rcx0r692HduLUDj7dMOKKjqi3lnzsko8EeNZvLbjN4+OxBqM7xZBdVcet/fqC0Sh7cEh1Tgx2uaa0PYPWqKURYlVV7mL/V6lHz7KHH12qnIbFOO//3o6Fc++73bMot45YPf+AvFw0n3iXdMIuORbphEBE3b0se1R4fJ3ZPoltSTNg+Jz3exUsXj6BroosNOaXc8uEPFEkbftHBSNIXEVdXtRPCG7gNyUyM5sWLR5KZ4GL9/lJ++c5asosa7rZZiPZGkr6IqOyiStbsLSHaYWPqwLQW+cyuSdH8/dITGJAex67CSn7xzlq+31vcIp8tRKRJ0hcRVVvKnzYwjThny9WvZyS4eOWnIxnXK5mCCjfXvb+Ot1Zmy5O7ot2TpC8ixusz+TQMbfODFe9y8OyPhnHZ6G54fSZ/XriDOz7eQF5ZdYvHIkRLkaQvImZJVgG5pdX0SI5mdI/kiMTgsNu4fXI/nr5gCAkuB9/tKODi11by8br9+KTUL9ohSfoiYv7z/X4AZg7vEvK2+U01qX8a7141mtP6plJW7eXROVu5+l9rpOsG0e5I0hcRkVNSxeKdBThsBucOa/mqnfp0TnDxzMyhPHrOINLjnWzKLeP699dx+0frWbevJNLhCRES8mSKiIhP1+fiM2HqgDRSY52RDqeOYRicMagzp/XrxL9WZvPmimy+21HAdzsKGNUtkctGd+fUvqk47FJeEm2TJH3R4jw+k49/sKp2LhyZGeFo6hcTZeeX43tx4cguvLt6Lx+s3ceavSWs2buR1NgozhmSwaUn9yLNacOIcNWUEE0hSV+0uP/tLOBAWU1Eb+AGKzXWyY2n9uHKsT345IccPv5hP1kFlby5Mps3V2bTMyWGyf07MbFfJ4Z2ScRhkxOAaN0k6YsW99E6q5T/oxGRv4EbrHiXg8vHdOey0d1Yt6+ETzfk8u32g+wurOSNFdm8sSKb2Cg7I7slMrpHMiO7JjKwczyxzuMfF0CIUJKkL1rU7sJKFu8owGk3WqTbhVAzDIOR3ZIY2S2J+IRoFm7M4ZttB/nfzgJ2F1ayJKuQJVmF1rpAr9QYVOd4BqbH0zMlhh4pMXRPjsHlkHsCIjIk6YsW9f6avZhYY+CmtKIbuM3hsNsY3SOZ0T2SuXNKP/LKqlm9p5hV2UVs2F/K9oMVZBVUklVQyVeb8+q2M7CeCO6WHE16vIvO8U7S4l2kxzlJj3eSGuskIdpBgsuBXaqLRIiFLekrpaYDFwIHAFNr/ccjlt8NZAL7gTHA/VrrzeGKR0ReWbWHT9dbT+BecmK3CEcTeunxLs4c3Jkz/YO613h8bD9YzubcMrbnl7O7sJI9RZXsL64ip7SanNLGn/yNd9lJjI4i0eUgIdpBnNNOTJSd6Cgb0Q47MVE2oqPs1o/DRkyUnaSEaGqq3DjsBlE2A4fdwGGzEWU3cNgMomw2/zwDh92GwzAwDLAFvNoM66rGgLr3ta+ibQtL0ldKxWKNrTtUa13tH1t3mtZ6XsBq8cAdWmtTKfVT4CngvHDEI1qH2etzqHB7GdMjiQHp8ZEOJ+ycDhuDMxIYnJFw2HyP18fe4ipySqrJK68mr6zG/2O9L6p0U1LloazaQ1m1l7JqL/si9H+oT+BJwFZ3YrBOGIYB1hzrfS0Da5sj+zYKPIkYh80/+nODWbehU1LttvVtZ7MZ+Hyt7+lrm83gvjMGhryxQ7hK+uOBXVrr2qLMYuAcoC7pa63vC1jfBpSFKRbRCnh9Ju+t3gvAJSd2j3A0keWw2+iVGkuv1Nhjruf1mZRVeyit9lBS5aGkyk1FjZcqj49Kt5cqt//V46MqYBq7jYoqNx6faf14fXh8Jm6vicfnw+M1cQfM9/pMTMBnmpim9eozwTzy1R+XzwRMEy9waK4Ih2qPL+T7DFfS7wyUBkyX+OcdRSnlBK4CbqpveVJSTLMvKe12G8nJx/7DioSOGNecjbnsK6mmZ2os557YvUl11R3xeNXq1Ixt7HYbXm/okwUcOgEcfmI4/ORQu17dNgFxeQKS2GGni3rWD5xtBswNvFiob90jlxzax+H/j1rhPF7Hw263kRwdRUyIW4CFK+kfAAKvaRP98w7jT/gvAr/XWm+vb0fFxc0f4CI5OZaioopmbx8uHS0u0zR54ZttAPxkZBdKS5r2O+1ox+t4tba4ak/vicmxFFW2krgCyhzJidGt6njVqo2rupmhpacn1Ds/XO3GlgC9lFIu//QE4HOlVKpSKhHq6v1fBmZprVcppX4cplhEhK3cU8T6/aUkx0RxwfDW+QSuEB1FWEr6WusKpdQNwHNKqTxgndZ6nlLqSaAAeBx4CxgG9FFKAcQBH4YjHhFZ/1i2B4BLT+xGTJQ8rCREJIWtyabWeg4w54h5dwW8vzBcny1ajx/2lbBydxFxTjs/OaFrpMMRosOTxwJFWP1z2W4AfnJCVxKi5VlAISJNkr4Imw05pSzaUYDLYePS0e3vYSwh2iJJ+iJs/rpoJ2A9fdua+swXoiOTpC/CYtmuQlbsLiLB5eDKsR37YSwhWhNJ+iLkTNOsK+VfMbY7idFREY5ICFFLkr4IuXlb8tmUW0anOGe77FhNiLZMkr4IqSq3lz8v3AHAteN7Srt8IVoZSfoipF5fvoec0moGpscxc3iXSIcjhDiCJH0RMnuLK3ljhfX07W+n9pcBQIRohSTpi5AwTZMn5m6jxmsyY3BnTuieFOmQhBD1kKQvQuLzjbksySokMdrBryf2iXQ4QogGyHPx4rjll1Uza4F18/aOyf1Ii3c1ssXRSqs87DhYzq6CSoqr3FS5fbgcNpJiHAztmUqGyy7dOAgRAvJXJI6LaZo8OmcrpdUeTumTwtlD6h0rp15ZBRV8uekAS7IK2ZRT2ugYTP3T4jilTwpTBqQxNDNBxmsVohkk6Yvj8s7qvXy3o4AEl4N7pg9oNBH7TJOF2w7yr5XZfL+vpG5+lN2gf1ocPVNiSItz4YqyUe32UVhZw+6iKrYeKGNbfjnb8st5Y0U2fTrFcsGwTC4Ynkm8S77GQgRL/lpEs23MKeX5b60nb+8/cyCZidENrmuaJgu25vPy/3ax46A1FFBslJ3pKo2pA9M5sXtSg236k5NjOZBfxtq9xXy3o4CvNh9g58EKnl24g78v3cVFI7tyyYnd6BQn/fsI0RhJ+qJZCitquOfTjXh8Jhef0JXJA9IaXDe7qJIn521jSVYhABkJLq4c251zh2YSG+T4n06HjZN6pXBSrxRundiH73YU8M7qvazOLua15Xt4d/VeLhvTnSvGdJeSvxDHIH8doslqPD7umr2RfSXVDMlM4NZJfetdz+318dbKbF5duptqj48El4MbTu3NzOGZRNmb33DMYbcxeUAakweksW5fCa8t282iHQX8Y+lu/vP9fn5xck9+PLLLcX2GEO2VJH3RJKZp8tjcrazdW0LneCfPXDAEl+Po5Lo6u4jH52xjZ4FVlXPW4M7cNrlvyLtYHtE1kVk/Gsb3e4t5/tudfL+vhGcWbOe9NXu5+bQ+TB2QJjd8hQggSV8EzTRNnl24g8835OJy2Hhm5tCjmmcWVbh57tsdfLohF4CeKTHcPa0/J/VKCWtsI7sl8bdLRvLt9gL+umgnOwsq+N2nmxjeJZHbJvdlRNfEsH6+EG2FJH0RtL8t2cXbq/bisBk8ef4QBmUk1C0zTZPPNuTy54U7KK7yEGU3uOaknlx5Uo96rwTCwTAMJvXvxIS+qcz+YT8v/28XP+wv4RfvrGX6wDRuOq0P3ZNjWiQWIVorSfqiUaZp8vy3O3lzZTY2Ax49ZxCn9EmtW77zYAV/mruVNdnFAIztmczd0/rTKzU2IvE6bAYXjuzKGYM68+aKPfxr1V7mbsnnm20HuXhUV34+ridJMdLHv+iYJOmLY/J4ffxp7lZmr8/FbjP44wzF1IHpAFS6vfxj6W7eWpmNx2eSGhvFbZP7MmNQ51ZRjx7vcnDDqX24cGRXXlycxX835PL2qr18tiGXX5zck4tGdsXZQlchQrQWkvRFg/LLa7j3042s2VuCy2HjifOHMKFPKqZp8s22g8xasJ2c0moAfjQik5tP69MqR8nKSHDx4AzFpaO68ey3O1i5u4j/+2YH76/Zx82n9WHaQLnZKzoOSfqiXmuzi7n3803kldWQHu/kyfOHMKxLIlkHK5j1zfa6NveDOsdz17T+DG8DN0pVRjwvXDScxTsLeG6hdbP3ns82MSQzgatP6sGk/p2wSfIX7ZwkfXGYKreXlxbv4u1V2ZjAqG6JPHbeELw+k0e+2sKnG3LwmdS1ub9wRJeQ9ptf6faSX1ZDabWH0moPVW4v8XHRVFRU47AbJEZHkRTtIDkmijinvckldMMwOLVvJ07ufehm78acUu6avZHeqTFcObYHZw3ujEPa+It2yjDNxrq5iqy8vNJmB5icHEtRUUUowwmJ1hiXaZqsPVDOI59vYndhJTYDrjqpBzOHZ/Lemn18+P1+qj0+7AZcMLwL103odVxt7n2myY6DFazbW8y6fSXsKqxkX3EVBRXuoPcR77LTIzmGbkkx9EyJZkB6PKpzPN2So4MusVe6vcz+IYe3VmbXVVWlxTm5YHgmM4dnHrNridb4ewSJq6naa1zp6Qn1/hFI0o+A1hbXhpxSXvoui6W7rCqb3qkx/Gp8L5ZkFfLFpgN4fNavYPrANK6f0LtZrXKq3F425JSybl8J3+8tYd2+EkqrPUetF2U3SI9zkhgdRXy0gxiHDUeUHbfbS43HR0mVh+IqN4UVbqo8vno/K85pZ0B6HIMyEhicEc+QzAR6psQc80Tg8fr4anMer6/Yw05/30A2Ayb0SeWC4ZmM75161E3f1vZ7rCVxNU17jUuSfivSGuLymSYrdhfx+vI9rNhdBFhVNqf0SSG3tJq1e60eMG0GTB2QzlUndT+sXX5jDpbX8P2+Er73l+Q355bVnTxqZSS4GNk1kZHdEhmQHk/XpGjS451HJef6jpdpmhRVutlTVEV2USVZBRVsOVCOPlBGfnnNUfHEOe11J4DBGQkMyUygS6LrqOoh0zRZnV3Mf77fz/yt+XUxxzntTO7fidMHdeaknslE2W2t4vdYH4mradprXJL0W5FIxpVTUsXnG3P5dH0ue4urAHA5bHRPjia3tJqyai8A0Q4b5wzN4PLR3emRcuwHmnymSVZBBd/vLeH7fSWs21vMnqKqw9axGVZ/+CO7JdUl+mNVnQRq6vE6WF7DlrwyNuWUsSm3lI05pRwoO/pEkBwTddiJYGhm/GFPGBdU1PDZ+ly+2nyALXnldfPjnHbG9kxm6pAMRna2TlatiXzvm6a9xiVJvxVpybhM02RLXjmLth9k0Y4CNuaU1i1zOWz4TBO399AhHpwRz8wRXThDpTfYW2VxpZv1OaVs2F/CD/utpFpSdXhVTUyUjWFdEusS/LAuic3u/TIUxyu/rJqNuWVszCn1nwjKKKo8+v5BeryTIf4rgcGZ8QzOSCA5JopdBRXM0XnM3ZLH9vzDY+maFM2Irol1P/3S4nBEcFB4+d43TXuNS5J+KxLOuEqrPGzLL/fXnRfz/b6SoxLykQakxzFlQBozT+xOuutQV8den0l2USXb88vZnl/B9oPlbM0rZ3dh5VH7SI93+hN8Ul11TagSXziOl2ma5JRWszHHOgFszC1lU04p5TXeo9btmhTNwHRrgJfuyTHEOe3klFSh8ytYvP3gUdu4HDb6doqlX1ocA9Lj6OcfHKZzvCukLZ0a0hG/98ejvcbVUNIPW5NNpdR04ELgAGBqrf94xPJo4GlgLzAAeFxrvSVc8bQnZdUeckqq2V9Sxd7iSrbnV7Atv5w9hZUUN5LgwarWGNMjiVHdkxjYOR6Xw8b+kmq+3ZbP9pxS9hVXWT8lVVTXc7PU5bAxqHM8w7okMqxLAsO6JJCRcHT9eGtmGAZdEqPpkhjNNP8Txj7TZE9hpf8EYF0VbD5QVnc8juRy2OiS6KKP0/ozqnB7OVheQ3GVh025ZWzKLTtsfYfNoGtSNF2ToumWFE1GgotOsU46xTtJi3PSKc5JSkxUi5wYRMcVlqSvlIoFXgKGaq2rlVIfKqWmaa3nBax2G7Bba/2kUmo48CpwWjjiiTTTNPGaVgsRj8/EV17N3sJKSqo8lFa5Ka7yUFbjoazaQ1m1l/IaL+U1HqulSqWb0ioPZTVeKt1eqjw+vL7gL35sQHKsg6RoJ9FRVuuTSreX5buLmLslv9HtMxNc9EuLo19arP81jn6dYttlO3abYdArNZZeqbGcNTgDAI/PJOtgBTsOlrOnqNK6cVxYyZ6iSgoq3GQVHH3VUx/Dv6/dhZX1XikFrhfjtBMXZSPW6SDOZSfe6SDeZSfB5SAhOoqYKBuxTjtxUQ5inXZinTZiohy4HDaio2yke00qyqpx2AwcdgO7Yb06bDYcNkNOKh1cuEr644FdWutq//Ri4BwgMOmfA9wLoLX+QSk1UimVqLUu4ThtOVDKVf9aC3BosO2j3xxzIO6jljWwcuuuHAMfUFDhoaDi6CsAuwHJsU5SY6PISHDRp3M8adEOuiZG15VIO/ooVA6bQf/0OPqnxx21zB7tRO8p5EBZNXll1RwoqyGvrJq8shryymoornRTUuWhwu0N+ntiAhU1XipqvFAe/DMLzWEAhmG9YhgYRyyjdlndvEMzaucfurgzjpp35KnFsBkEVicb9ax1PBeLzdrUsE72Pn9crel0aDMMbp3Yl3OGZoR0v+H6i+4MlAZMl/jnBbPOYUk/KSmm6U9dFlYe1TywvbAbVlt2V13Jzo7LYcPlsOF02Ih22ImPdhDvcpBwxGu8y0FKrJO0eCed4l2kxERhCyj12e02vN76275Hkt3fPLK1sdttjBmQ3uh61R4fJZVuiipqKKp0U1ThpqzaOhlU1HipqPZYr24vZVUeSqqsq7tKt5dqj/V8Qo3Hh9tr4vFZV4ten4nPBJ/PbHbBwwTqcnBQ9/ba599Ua1bk9ob8ux+upH8ACGzUneif19R1KC4O7vI5UJrLzlmD0nG6HLjdXmyGVQ6xGVZdroF1Fq09l9TOtxlg85dmbP4SAAHzjbr5YPdfJtsMG1F2G067QZTdwOmwEWWzEeUwcNkMnA67Nc9us7YxICUplvLy6npLWrWx1haqDH8cLod1aR6yenO3hxL34aX/9npDK1yaElcUkO6yWzfKk0PbxNPrM/H4TH81okliYgwFhRX4TP+JwbROEmbAum6vjxqvD7fHeu/2mfhMax+125iAzwcmJqYPfFjLTNPENK3zhJfa99Zr7WeamIddXPtMiIlxUllZXTvrsPNM7fu6c9ARJ5hGz0n1LD9qH0cu98+IiYmistId9NV9S4mOieKUHknN/u6np9f/XE24kv4SoJdSyuWv4pkAvKCUSgU8/iqcz7GqgRb56/S/D0XVDkBKrJOHzhncupNF+6sSFxFiP6KePiE6Cm8rHC+gVf89dqC4wpJ6tNYVwA3Ac0qpR4B1/pu4vwNu9K/2Z6wTwx+AO4FfhCMWIYQQh4TtLp3Weg4w54h5dwW8rwRuCtfnCyGEOJpUMgghRAciSV8IIToQSfpCCNGBSNIXQogORJK+EEJ0IK2+l00hhBChIyV9IYToQCTpCyFEByJJXwghOpB22W9uYwO4tGAc/YBHgNVAd+Cg1vohpdSDwOSAVR/1P8HckrEtBWpHBvFqraf5+0Z6HNiBNbDNvVrr3BaOqzdWF9x7/LMSgXVAFi18zJRSmVi/v5Fa67H+eQ0O/qOU+hkwCvAC27XWL7dgXHcDmcB+YAxwv9Z6s39ZFtbxA9irtb68BeO6GrieQ9+1V7XWb/qXRfJ4vQr0C1htODBaa53VEsfrGLmhwb9BpdRvsf4eUoCvtdazm/PZ7S7pBzmAS0tJBd7VWn/ij22jUupzAK315AjEE+hLrfWDR8x7DJirtX5fKXUeVnK7ooXjKgWu01rPBfCfIOcC0yNwzE4FPgFOCJh3G/UM/qOU6g78BhiltTaVUiuUUvO11ltbKK544A7/Z/8UeAo4z7/stXp+1+FQX1wAl2itswJntILj9bXW+j1/LIlYx6g2xpY4Xg3lhmup529QKTUOmKK1Plsp5QA2KaUWaq2Lm/rB7S7pE9wALi1Ca73iiFk2oBxAKfV7oBqwA8/7O6lrScP9pcMYYIXW+nOs4/Sof/li4PUWjgmt9UGsJI9SygWM0Vo/qJSa3tLHTGv9b6XU5CNm1zv4D3AmsEprXdscbglwFhDyJFZfXFrr+wImbUDgWI2nKaXuwurK/Aut9f9CHVNDcfndrJTKAWKBv2itC4j88XovYPLnwD8CpsN+vI6RGxr6GzwX6xihtfYopTYBk4Aml/bbY51+MAO4tDil1I+Ar/yX3B8Az2qtn8aK9fkIhPSE1voJ4GHgXqXURA4/diVAir9UESmXAu/637eGYwYNf79axfdOKeUErgL+EDD7Hq31k8CfgH8opfq3YEgLsb5rTwMrsX6P0HqOlw3rBPR5wOwWPV5H5IaG/gZDdrzaY9IPanCWlqSUmgJMAW4H0Fpv0FqX+xfPB6a2dExa6+X+Vy+wyB9f4LFLBAq11o2PtB4+PwHeg9ZxzPwa+n5F/HvnT/gvAr/XWm+vnR/wu64A1mKNb9EitNY7tdZ5/sn5wCSllJ1WcLz8zgc+D7jiaNHjdWRuoOG/wZAdr/aY9OsGcPFPT+Dws3iLUkqdg1WS+DWQqZQar5R6KmCVAcD2ejcOX0yDlFKB4xfUxlA7sA1E/rhNBpZord3+6YgeswB1x+iIwX++AkYrpWpHMxkPfNFSQfnvZb0MzNJar1JK/dg/f5pSakbAqv1pwWOnlPpTwNXiACDLX9CI6PEKcBXwWu1ESx6v+nIDDf8NBn7vooDBwLfN+dx2+USuUup04CIgD3BHsPXOaKzL25X+WXHAXwGFVb95AKvVwP21LUBaKK6uwF+ANVglhijgDiAZeALYhdWy4Xct3XonIMZ3gFu01vn+6T/RwsdMKTUJuBKYgVWCfsa/6GmsVjL9gceOaL0zBqs1ypYwtkapL65/AcOAff7V4rTWY/0npgeBVUBXYJ/W+rEWjOtX/rh2Yv3e/qy1XupfP2LHS2tdqZQ6Abhca/3bgHVb5HgdIzfMpoG/QX/rnRT/zxfNbb3TLpO+EEKI+rXH6h0hhBANkKQvhBAdiCR9IYToQCTpCyFEByJJX4jj4G9z3uY/Q3Qc7bEbBtHBKKUWAcuATlgd7f3NvygN66nKP2qtLwnD556C1Zb6qSPmPwmcFMK+gm5QSm3QWi8I0f5EByYlfdEe/ENr/RusduGFWuvf+KcXaq01VncOIaWUSsZqrz+rnsUvhPjj/go8oJRKD/F+RQckJX3R5mmt/9nQfKXUrVgPnvVWSl2DVfJ/BhiBdSXwT6ynIgcA52qtS5RSQ4G7gR+AQVjdOO84YvcXAcv9T5eilDoR+COwHHDXrqSU6ot1Yvgf1sNJzwDZwNuAC6sX0+7Ac1i9Tk7BelgJwKm1/oO/F8rFWCev55p1kITwk5K+aNe01s8FvP8nsBlYrbW+AqvHzgSt9S+wnk4+3b/q34GXtNZPAW9y6EncQEOx+tSv9TLwiNb6YQ49ZQlQAzzk78Dr/7D6xckHbgRStda7/fv5UGv9DdYTrB9orR8C/huwn73+zxTiuEhJX3REtX2pFAW8L+RQh1YjgDP8PY/GcHg3xbVcQGBndEM51C1w4FWBG7hEKXUWVpcX6QBa621KqSxlDfgziUMl+EuBx5RSGf55/wvYT0zT/ptCHE1K+kIc7XvgP1rrx7EGlqmv47k9WANh1NoIDPS/7xsw/3dAmdb6UawBVwL9GasaKTagJ8oErfWPgJlYVwa1UoHdTf+vCHE4SfqiXVBKxWBVjSQppX4eMP9G/7xL/R3x9QKu9tfBj8AalWggMBE4z1/C/gVwq1Lqd8CTWHXwR/oYGBcwfT3WzdaHgelYPb2eA3wInK6UehRrVKReSqlpAP7RwboAbwTs52r/YDG3cviYASf79yXEcZEO14RoJqXUA8Da2iHvmritC6u+/69a6xsbWfd0rKHy7m1epEIcIiV9IZrJ32X3zmZu/k+s6p33GlsRyJWEL0JFSvpCCNGBSElfCCE6EEn6QgjRgUjSF0KIDkSSvhBCdCCS9IUQogORpC+EEB3I/wMcdobcYxPxlQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot\n",
    "\n",
    "labels = ['Susceptible', 'Exposed', 'Infectious', 'Recovered']\n",
    "for response, label in zip(responses.T, labels):\n",
    "    pyplot.plot(time_span, response, label=label)\n",
    "\n",
    "pyplot.title('SEIR model')\n",
    "pyplot.xlabel('Time (days)')\n",
    "pyplot.ylabel('% of population')\n",
    "pyplot.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b06905a0",
   "metadata": {},
   "source": [
    "### Stochastic model\n",
    "\n",
    "We now have our deterministic base line model, and can observe that it works. \n",
    "Let us change the assumption to assume that the parameters are random, and model it using polynomial chaos expansion (PCE).\n",
    "\n",
    "We start by generating a PCE bases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "623b3315",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.353795Z",
     "iopub.status.busy": "2021-05-18T10:57:22.353457Z",
     "iopub.status.idle": "2021-05-18T10:57:22.416346Z",
     "shell.execute_reply": "2021-05-18T10:57:22.416617Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q2-0.5, q1-1.0, q0-0.2, q2**2-q2+0.24917])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polynomial_order = 3\n",
    "polynomial_expansion = chaospy.generate_expansion(\n",
    "    polynomial_order, distribution)\n",
    "polynomial_expansion[:5].round(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09a1167a",
   "metadata": {},
   "source": [
    "Generate our quadrature nodes and weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fbc87293",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.418762Z",
     "iopub.status.busy": "2021-05-18T10:57:22.418429Z",
     "iopub.status.idle": "2021-05-18T10:57:22.671882Z",
     "shell.execute_reply": "2021-05-18T10:57:22.672183Z"
    }
   },
   "outputs": [],
   "source": [
    "quadrature_order = 8\n",
    "abscissas, weights = chaospy.generate_quadrature(\n",
    "    quadrature_order, distribution, rule=\"gaussian\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76b585e4",
   "metadata": {},
   "source": [
    "We wrap the deterministic model solution into a function of the model parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "484f6711",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.674412Z",
     "iopub.status.busy": "2021-05-18T10:57:22.674100Z",
     "iopub.status.idle": "2021-05-18T10:57:22.681548Z",
     "shell.execute_reply": "2021-05-18T10:57:22.681830Z"
    }
   },
   "outputs": [],
   "source": [
    "def model_solver(parameters, delta=1e-4):\n",
    "    return odeint(ode_seir, initial_condition(delta), time_span, args=(parameters,))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76985eed",
   "metadata": {},
   "source": [
    "Now we're going to evaluate the model taking parameters from the quadrature. To reduce the computational load, we use multiprocessing to increase the computational speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bdc6e583",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:22.684049Z",
     "iopub.status.busy": "2021-05-18T10:57:22.683724Z",
     "iopub.status.idle": "2021-05-18T10:57:23.015641Z",
     "shell.execute_reply": "2021-05-18T10:57:23.015948Z"
    }
   },
   "outputs": [],
   "source": [
    "from multiprocessing import Pool\n",
    "\n",
    "with Pool(4) as pool:\n",
    "    evaluations = pool.map(model_solver, abscissas.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c4c6f9f",
   "metadata": {},
   "source": [
    "And finally we're calculating the PCE Fourier coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3a3c4a38",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:23.018270Z",
     "iopub.status.busy": "2021-05-18T10:57:23.017954Z",
     "iopub.status.idle": "2021-05-18T10:57:23.939558Z",
     "shell.execute_reply": "2021-05-18T10:57:23.939235Z"
    }
   },
   "outputs": [],
   "source": [
    "model_approx = chaospy.fit_quadrature(\n",
    "    polynomial_expansion, abscissas, weights, evaluations)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6d180ca",
   "metadata": {},
   "source": [
    "With a model approximation we can calculate the mean and the standard deviations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b807a488",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:23.942443Z",
     "iopub.status.busy": "2021-05-18T10:57:23.941758Z",
     "iopub.status.idle": "2021-05-18T10:57:24.175127Z",
     "shell.execute_reply": "2021-05-18T10:57:24.174793Z"
    }
   },
   "outputs": [],
   "source": [
    "expected = chaospy.E(model_approx, distribution)\n",
    "std = chaospy.Std(model_approx, distribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44dfedb4",
   "metadata": {},
   "source": [
    "Finally we can plot the data with uncertainty intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4aae4e70",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:24.177938Z",
     "iopub.status.busy": "2021-05-18T10:57:24.177619Z",
     "iopub.status.idle": "2021-05-18T10:57:24.278937Z",
     "shell.execute_reply": "2021-05-18T10:57:24.279164Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f477d931730>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAERCAYAAACXT3dwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABj50lEQVR4nO2dd3wkR5X4v909eTSjUd7V5uTaXa93ne21cbZJxj4OOA4DPgwcHPkHhiOYA0w2Bnyc4cgccGQ4ksk457y2N9euN+8qZ00O3b8/ukcaaUfSjFYjaaX6fj6t6a6urn7T03pd/erVe5plWSgUCoVifqDPtAAKhUKhmD6U0lcoFIp5hFL6CoVCMY9QSl+hUCjmEUrpKxQKxTxCKX2FQqGYRyilrygZTdOWaZr2e03THtA07W5N0+7TNO2dFTjPxzVNa9M07eapbnvUeZZrmrZzkseOeS00TTtN07THNE2zNE173NmXX3Y7da7XNG23pmkpTdPuc8q+7Hzvdqfu45qm7dU07Y1T9qVL/34fLfU3mK7fSzFFWJalFrWUtAD3Au8o2L4U2Faw/QPg5ik615S1VdDmQeDSUWWRCl2L5YAFrB513H0F6zcAR4t87x8XbL8WyI5uZ5p+75J/g0r8XmqpzKJ6+opyOBe4L79hWdZ9wE9mSpipwLKsvkkeOtlr8cEyz3MHYACnl3mcQlEUpfQV5XAI+KCmacF8gWVZtwBomvb/gBcDNzimiTc75edomna/Ywa5X9O0c/LHaprWqGnar5x9j2qa9p+apvkLzlevadqPNE3brmnaDwuOW6Rp2q81TbtT07SHC80KmqY1aZr2F03T7tU07SFN0z7klH8fWAB8xZHvLMcsY2mattyp49I07RZN0x5xZPqlpmkry70WxXBMST+wLOuJiS7yKFzO59EibXqc72JpmvZOTdP+6piDLtY07f3OtXlU07SGgmPG+z1WaZr2oHPM9wH/qPOdVXDc3ZqmrS3zuyhmAzP9qqGWk2cBrgC6gT7g+8Alo/b/gIJXfKAa6MIxqQAXOdsRZ/vvwCecdQ+wBVhe0NYzgBfwOefd7OxbA7yk4Dz3Alc467cCH3LWg8BDBfUOcrx5xyo4503AnYDhbH8NuGGS12K50/bj2G8EjwE/GFXnBiY27/wr8PkJfhcLuNFZfy9wBDjf2b4D+HCJv8fjwEec9UVAT/73dI7tBC53tq8GJKAX++3VMnsX1dNXlIxlWXcDS4EPYCu1ezVN+9Y4h7wMGLBs0weWZT0I9ALXapq2CLgKW2FiWVYaeBO2Yslzr2VZKcuyksAeYIVTfgS4wumR3wesA85y9vUAL9E07VTLsmLAC8v4im8EfmRZVs7Z/hxwf7GKZVyL11mWdSnwmjLkuMrpwR8EPoqtUCfiTudzOxCwLOsxZ3srkH9bGe/3WIZtsvqxs+8Y8FBB+y8DopZl3ePs/xP2m9N5ZXwvxSxAKX1FWViWFbMs67uWZV0GXAb86zgmkMWMVOI424udhcL9lmU96yjqPAMF6ynstwGADwMXY/fuLwX+CgScfV8Efg38QtO0Z7F7pKUyQl7LslosyzowVuVyroVlWQcty7qhRDnudL7X6UAa++EzEYPOZ7ZgPb+dv27j/R4Lne2ugn09BeuLgdpCTyTn2LoSZFPMIpTSV5SMpmnfKNy2LOt+bBNH9RiHHAEaRpU1YNunjxRs59tfqWlaTQminAs8YFlWwtl2F+xrtCzrq5ZlbcDuhf9I07RVJbR5nLyaptXl7f2jmcS1yB93SYmyYNmDzB8F/mEsOcpkvN+jtWA7T6FCP4Jtiro0vwBnYpvoFCcRSukryuFKTdPOzW84CswEdjtFg0BA07Sgpmk/Af4IhDRNu9ipfyFQA9xhWVYLtkniBmefF/glw73S8XgeOEfTNN0ZSH1Bwb7Pa5p2urP+OHZPWRsl32XOwPNofgBcr2ma4WzfAmwaQ4aJrsVYfHKC/aP5Dfb3fVeZxxVjvN/jEPAEcL2zbxFwyahj6/MDv851v5cJHnKKWchMDyqo5eRZgLdg27jvBR7AHqA8v2D/Zmyl9wTwWqfsLKfeA86x5xTUbwR+5ex7GHiVU34j0IY98PqPwMexB0x3A5dj25LvBZ7FdpO816n7WmxzzgPAPcDTwHsKzvcuYCf2oOqpwN3Yg6CPYQ9curEV/aPY9uzPTuZaAGuBPztt/wX4v4Jlh1Pneuf7pHB894EvO9+7HfjfgnO9EYgCfywix98LvoNwrkkS+BZwrXNd2hge6B3v91jlfO/HgJ9hP3AOAm8ddez9zvEvc8o/XvB7vXmm71O1jL9ozo+mUCgUinmAMu8oFArFPEIpfYVCoZhHKKWvUCgU8wil9BUKhWIeoZS+QqFQzCNcE1eZWTo7ByftXlRd7ae/PzFxxWlGyVUeSq7yUHKVx1yVq6EhpBUrn9M9fU0r+p1nHCVXeSi5ykPJVR7zTa45rfQVCoVCMRKl9BUKhWIeoZS+QqFQzCOU0lcoFIp5hFL6CoVCMY+oiMumEGIB8Blgk5TynCL7dezEEIPYWYe+J6V8bHQ9hUKhUEwtlfLTfwHwe+zMP8V4NRCWUn5YCFELPCaEWCelzI1RvyxSWZO/7monGPSSiKfRNQ1NA00Dnfy6hgboBeuaU08vst+l24vb0Ed9OvtGlRv67HQDUygU85uKKH0p5f8JIS4dp8rVOBl3pJQ9QogkdnzzrVNx/ng6y2f+vncqmpo0HkPD7zYIeAwC+U+PQcDjIlLlwa9r1AY81ATc1Abc1AQ81AXcNFR51QNDoVBUjJmakdvIyDyeA07ZcVRX+8uepGD4PLzyzEUAmKaFZYFpWVjYn1hgWmBh2Z+WnUxmuA5gWcN1TMiaJpmcRSZnks1ZpHMmWXN4O5MzncVeT+cs0rks/clsWbK7dI2GkJcFYR/NER9LawOsWxBCLAhRF/TiNuy3Cbcx9cMxhqETiQQmrjjNKLnKQ8lVHvNNrplS+h1AqGA77JQdx2SnIX/4slVEIgH6+uKTOv5EsCyLVNYknskRT+dIOJ/5bdMwaO2O0ZvI0BNP0x1L0x3L0BVL05fI0NqfpLU/yTNHRrZb43eztNbP6vogq+uDLKr24nEZeF06HkPH49LxuXS8Lh2f2y7Xy3hgztT1mgglV3koucpjrsrV0BAqWj5tSl8IEQQCUspO4E/AxcCPHJu+D9gxXbJUGk3T8LkNfG6D2iIP6vF+zGQmR+tAisO9cZ7virG3M8a+rhhH+5L0JjL0Hsvw3LEBAMI+F6ctDLNpUZgVdYHjFLwGeF36kJkp6DGo8roIeIyyHgYKhWLuUCnvnUuwc4AuFEL8B3buzxuA04C3YSfAPkMI8QlgKfAvUzWIe7LjcxusqAuwoi7AJavrh8oHU1meO9rPE4d72dUeZV9XnIFklocP9PDwgR5CXhfnLouweXkNNQE7t7gFJLMmyaxJbyIz1JauQcjrotrnJuJ3E/a71ENAoZgnzPocuScSZXOuvrb1xTMc7YvzXOsgW48N8FzLAN2xNGD37jcsDHGlaGBxxF9Se4auURtws7o5gpHNzroHwFz9HSuFkqs85qpcY0XZnPWhlRXHEwm4iQSqWV4XZFNzmJfG0hzsSfDwgR62HhtgW+sg21oHOW1hiKtPbaKhyjtueznTojOaJn6sn0wyzYKQj+ZqHx6XmrunUMw1lNI/ianyutiwMExPPE3AY5uFBjZkuHdvN48c6GFb6yA726NctrqOK0VDSR4/6ZzF4b4Ex/qTLKz2sjTix1UBTyGFQjEzqP/mOUBtwMOZiyM0h32EfW7+4bQFfPSFazhnaYScaXHXni6+dM8+jvSV7gmVsyyO9iV58kgfbQPJCkqvUCimE6X05wiGrrG6Icj6BSEMXSPsc/OaMxfxrotWsCDkpSuW5qv3H+C+vV2UM46TyVns6Yyxo3WATM6s4DdQKBTTgVL6c4z6oIfTF4XxOfb4FXUB3nvpSi5cUUvOsvjDjnZ+/NTRshV4dzzDliP9DJY52UyhUMwulNKfgwQ9LjYtqiboMQBwGzqv2LSQN563BK9L59ljA3zjoYNEU+Up8FTOZGvLAF2Op5BCoTj5UEp/juJ16WxsDg8pfoANC8O8+6IV1PjdHOpN8PWHDjKQzIzTyvHkLItdbYPKzq9QnKQopT+HcRs6py0M43cP/8wLq3285xLbzt8+mOLrDx2kP1Ge4reAvZ0x2gdTUyyxQqGoNErpz3E8Lp0NC8O4jeF5GmGfm7e/YDkLw146o2m+/eghEunyJkRbwJ6O6NCkMIVCcXKglP48wO82WNcUonB6XpXXxdsuXE5TyEvbQIrvP3647MFdC9jdHi17bEChUMwcSunPEyJ+NyvqRkZ/q/K6+NfNSwn7XOzrjvO/jx4qy50TbBv/jtZB0lnlzqlQnAwopT+PWBzxUxdwjyirDXh4y+ZleF06Tx7q5YF93WW3m8qZ7GofLPuBoVAoph+l9OcZpzRW4TFGxmFqrvbxmjPspDN/3NHOvq5Y2e32J7Mc6p1c7gOFQjF9KKU/z3AbOqvrg8eVb1wU5kXrmzAt+NGTRydlpz/SmyjbE0ihUEwvSunPQ+qrvNQHPceVX7upmZV1AQZTWX71bEvZ5hoL2NMZJWcqM49CMVtRSn+esqo+eFwCdkPXuO6sRfhcOttbB3nicF/Z7SYyJgd6Zl9scoVCYaOU/jzF69JZVnN8kpXagId/3LgQgN9va6M3Xr65prU/WfZMX4VCMT0opT+Paa72jZitm+esJdWctjBEKmvy262tZbebn7GrvHkUitmHUvrzGF3TWFEkc7umafzjxoV4XTo72gbZ1jJQdtuxdI5j/So+j0Ix21BKf55TX+Ul7D0+gVq1381L1zcC8NutrSQz5eetP9ybUDH4FYpZhlL6CpYX6e0DXLCilqU1fvqTWf6+u7PsdrOmxUE1qKtQzCqU0lcQCbip9h3f29c1jVduWogGPLS/Z1Jx9NsGUsTSKjaPQjFbUEpfAcCymuK9/cURP2cvjZCzLP64va3sdi3gYLfq7SsUswWl9BWA09v3u4vue8m6RjyGxrbWwUmFaOiOZ9RMXYVilqCUvmKIZWPY9qv9bi5dUw/AH7a3TcoVU9n2FYrZgVL6iiEaQl4CbqPovktX1xP2uTjSl2R762DZbfcns/TEVcIVhWKmUUpfMYLFEV/Rcq9L54pT7N7+X3d3YE6it3+oR0XhVChmGqX0FSNoDHlHpFYs5PxlNUT8btoGUjx3rPwJW4OprEqvqFDMMMf76SnmNbqmsTDk43Df8b1yl6Fzpajn/55t5W+7O9jYHD4uaNtEHO5NUFckwudcx7IsUlmTdM4knbPI5EwyOYusaZLNWZiWRc60yFl23ZxlYVlgOp8WgAWWvYY1ou2R56qqihGNTm/S+lLe+6qCMaKxVGmVxzzP1If2CAajxGLTe71KYbNwUdzYemIopa84jgVhL0f6EkX/vc5dWsO9e7rojKZ55mg/Zy+NlNX2YCpLbzxNTWBuKn7LsohnckRTOWLpLPF0jkQmRzJjVkBdFcd+mMy+uEdZ05qVYbdNy15mHRWKXaWUvuI4fG6D2qCnqCnG0DWuWtvIz7cc4649nZy5pBpdK7+3P5eUfiydpTeeoS+RoT+ZnZWKTaHIUzGlL4S4EngF0AFYUspPjtq/AvgS8CRwOvBTKeUdlZJHUR7NYe+Y9vczF1fzt10ddEbTbG8dZGNzuKy2+5NZ+hOZMecFnAwkMjm6O6Psb+0nkVHxhRQnDxUZyBVCBIBvAu+TUt4MbBRCXDGq2geBh6SUtwBfAL5cCVkUk6Mm4MHnKn57GLrGpWvqALhnT9ek/PaPFhkzmO1YlkVXLM1zx/p58nAfB7vjSuErTjoq5b2zGTgkpcyPjjwMXD2qTjvQ4Kw3AE9XSBbFJFkQLu6+CbZtP+gxONKX4PlJztI9mWLydAymePpIPzvbBulPnjxyKxSjqZR5pxEonMEz4JQVchvwWyHEbcC5wKeLNVRd7Ucr02acxzB0IpHis0xnkpNFLn/QS1c6V3Q8KQhcsbaRO7a2cv++Hk5fXlf2+fpzsKiE6zCT16svnmZPR5TBZBbN4yLoGf6XMXSNYNA7I3KNRyXksiyLTM4imc2RypgkMzlSWXPICynjeCFlC7ySMo6XUta0P03LGvJUMi0wzWEvpZxlYZrD+3LmyHr2+rAnk+W4Mw1t25tYhWXOjXvc9lCZ5Rxz/PbI716wPs5w/HgvvCO9rcZpo2BdYwffuf4sLl7TMGb9yVAppd8BhAq2w05ZIT8Aviul/JkQogHYK4RYKaXsKazU3z95M0AkEqCvb/ZN/z+Z5PKaJt1jpEw8Z3E1f9vZzq62QXYf62NJ5Pj0i+OxP56izq3jHcOMNJ5clSabM9nXHad9cGxXvmDQOytd/SaSK5016U9mGEhmiaVtL6NYKkc8v57OEUvbSj2v5NM5c3Z6uMxhLCAaTU363m9oCBUtr5TSfxRYJoTwOiaeC4GvCyFqgayUcgBYAuRz8fUCJmqy2KyjKeQdU+kHPAbnL6/h/ue7uXdvF/9yzpKy2jYtaOlPsqJudr319CUyyPYoqZM0AUwqm6OlP0lnNEVXLE13LD00eN6XyJKYREIcAJeu4XXpIxa3oePSNdyGjqFruHUNl6Hh0nXchjZcpuv4/W6ymSy6pmFoGrpmZ2kzNA1dt+eI6Fr+01nXtRHlGoAGGgxZAOx1nH12ncJ6YB/PUD2tYB9UBb3E4+mhdij4HEYrsjZW3WJHjdzQjmvl+PYuOKURd25yv9V4VETpSynjQoi3A7cLITqBrVLKu4UQtwI9wC3A+4D3CiEuAFYAN0kpuyohj2Ly1AY9uHRtTL/vi1fV8eC+bra1DNAbz1ATKM8jp3UgydIaf9mTvCrF4d4Eh3ri0+ZTfyLkTIv2wRQt/UmO9Sdp6U/SEU0xMMGYg6FpVPtdhH22uSroNexPj0HQYxDw2Nt+t47HpeNz6Xhdxgn/RrP2zcjrguzUK9cTxdA1NHPq/y8q5rIppbwTuHNU2QcL1h8CHqrU+RVTg65pNFZ5aRkonu824nezsTnMs8cGePhADy87tams9rOmRdtgkkXV5ZmGppqcaSE7opNKFDNdDCazHOiOc6AnxoHuOC0DqaJzAgxdoy7gpr7KS0PQQ13QQ8TvptrvotrvJugxyp5boZg7qMlZiglpDHnGVPpg9/afPTbAYwd7uUo0TGijH01Lf5LmsG/SA/YnSiZnsqN1kIHU7PLKyeRM9nXF2NUeZXd78QdSfdBDc7WPRdU+mqt9NIW8LK6vIpGYvQ8vxcyilL5iQsI+N363PqZP+rLaAEtr/BzuTbDlSB+bV9SW1X4iY9ITz8xITJ501mRb6wCx9Ox4vU9mcmxvHeS5Y/3s7YqRyQ335D2GxrLaACtqA6yoD7A04sdXJBS2PktMZYrZiVL6ipJoqPJyuHdsT6qLV9Xx46eO8uD+Hs5fXlN2r/1Yf3LalX46a7K1ZYD4JAc2p4qcabGzbZCnj/Sxqz06YvxkcbWPtU1VrG0KzaqxD8U0UKE3X6X0FSXRFBpf6W9sDlPtc9E+mGJPZwzRWFVW+30Je7JWoR98Jcnm7B7+TCr8/kSGxw728tih3qHBVw1YVRfg9MXVbFgYIuw7eUNVGJqGodtjDLqmUe134zFNxyPneG8draAsv5730tELvHPy64VeOHmvnNHeNyPKCmQr9PypjgQY6I8fV2dEQ8cdP/b3HmtXuR2hSNBDX9/UmxyV0leUhN9tUOU1iKaKK0lD17hwZS1/3tnBA/u6y1b6AMf6kpwyiePKxbQsdrQNzphJp20gyd17unj2WP+Q73tDlYfzltVwxuJqIrMwJpHb0PAY+tDidmm4HbdMl+O2WbgYunackput81PyrqfzBaX0FSXTWOUlmhr7n/b85TX8fXcnu9ujdEZTNFSVNyu0M5pmZZ2Jq8L/gHs6ojMSSuFoX4K7ZCfbnHSTuma/IV2woobV9cEZG8gG2wff7zYIuA28bttN0+c28Llst03l7TN3UEpfUTL1QQ/7u8dW+kGPizMXV/PE4T4ePdjLtRsWlNV+zrJoG0yxuMyZveVwpDdBR3R6PVv6Ehn+vLOdp4/0A7aCPXdZhMvW1FM7zSGmNexJdVVD/vm2T76nTI8rxcmLUvqKkvG5Dap9rnF7yResqOWJw308caiPl6xrLPu1uaU/yaLqyrhv9sbTHOyZPvNCOmtyz94u7nu+i0zOsk1gK2q5bE3dtNnqDU0j7LP988NeF1U+Fy41GDyvUUpfURb1VZ5xlf6SGj9LIj6O9CV59tgA55SZWSuZtWP91E+xJ08qayI7otM20/b5zhi/fLZlKCfBpuYwV5/aNC0eSlVeg1q3nxq/m7DfpUwzihEopa8oi/qgl31d4/eWL1hRyy+eaeGRAz1lK32A1v7klCt92RElnau8yk9mctyxvY3HD/UBdurJV21qrnh8oaDHoLHKS32Vh4UNoVk5YKqYHSilrygLr0uf0MRz+qJq7tjexuHeBEf7EmXb6HsTGeLpHAHP1KSFPtqXoC9RPGjcVHKkL8GPnzxKVyxtp5U8pZ7LTqnHpVfGXu427BAZTSEvVV71r6woDXWnKMpmIhOPx6VzztIID+zr4ZEDPbz6jEVln6NlIMnq+uCJiAlAPJ2ruB3fsiwe3NfDH3e0k7MsFoa9vP7sxeMmoTkRqrwGzWEfjSGvMt0oykYN2SvKpr6EBB2bl9uhGJ452k9iEv7w7YOpMSN7loplWezpiFY0DnwmZ/Ljp47y++1t5CyLC1fU8v8uWVkRhR/xuzltYYgzF0dYEPYpha+YFErpK8rG69IJT2BOaAx5WdMQJJ2zeOpIX9nnyJkWndETC8PbMpCsaBC1vnia/37wAM8eG8Dr0rnh3CW8YtPCKZ/oU+1zsbE5zMbmMDXT7OKpmHsopa+YFPVVEyufC5zAa48e7J1U8vSW/rEje05EKmtysKdyydeP9SX4/F8lR/qS1AbcvOfiFZzWHJ7Sc/jdOusXhNi0aHbO0lWcnCilr5gUpXjXnLogRNiJxzORx08xYuncpAdgD3THisaanwoOdMf5+kMH6U9kWFkXmHJzjq7Bsho/Zy2JTLkXk0JRttIXQpSfAVsx5/A5sXjGw9A1zl9WA8AjB3rGrTsWrePE8R+LvkSmYrNuZUeUbz9ykGTW5IwlEf7tgmVT6jlT7XNx1pIIy2oDymavqAgT3q1CiCrgKoYTnV8D/FMlhVKcHNQHx4/FA3Y8nrv2dLKtdYCBZKbsmajdsTSpMiJhWpbFvq5YWecolR2tg/zwiSPkLItzlka44cIVJKcoWYmuwfLaQEVDUCgUUFpP/4/AFdh5bFcA5WXIUMxZ6oITK/Bqv5tTF4QwLXjCmbBUDvnk6aXSPpiqSPTMPR1RfvikrfAvWlnLq89onrLY9gG3wemLqpXCV0wLpbyXHpBSviu/IYRYUUF5FCcR+eTZY2XUyrN5RS3bWgd57GAvl59SX7bZ4lhfgkidf8J4PDnT4lAFBm8PdMf5/uOHyZm2S+Y/nLZgymIDNVR5WNNQpeLhKKaNUnr6B4UQVwkhlgkhlgJvqLRQipOHUnz21zQEqQu46U1kkO3Rss+RypolJSxv6U+Syo3/ACqXY30JvvvoIdI5i7OXRHj5xqlR+BqwojbAuqaQUviKaaUUpf824CbgB8APgesrKZDi5KIUE4+uaZzvuG9OdkB3IhNPNmdypG9qe/n9iQzffewwyazJxuYwrz6jeUoGVw1NY92CEEtqlDlHMf2UYt75iJTyB/kNIcSVlRNHcbIR9rnxGvqEPexzl0b4664OdrVH6Ymny44j35/MjptO8Wh/8oRn8BaSyub43mOHGUhmWVkX4HVnLZoSG77b0NiwIEzIpyKgKGaGCXv6UsofCCGuFEK8XwhxhZTyrukQTHHyUFtCb7/Ka88qtYDHD/ZO6jyt/cVn6GZzJsdOYCLXaEzL4idPHRtK1v6Gc5dMSTYvr0vn9EXVSuErZpQJ72QhxMeAG4FlwAecbYViiFInEF2w3PbZf/xQ36QmTrVHi8fjOdqfnNKJWH/d1cGOtkH8bp1/PX/plPjh+906m5rD+N1TEzlUoZgspdzNHinlS/MbQojPV1AexUlItd+NoWsTKt4VdQEWhLy0DabY3jrApkXVZZ0nZ1q0DyZZVD1sC8+a1gmFaxjNzrZB7t7ThQa84dwlNIbKy/NbDL9bZ2NzNV6VklAxCyjlLhxtrJ1a9wjFSY+uadQGJjbxaJrG5hV2b//RSZp4Riv41im05XfH0vz06aMAvHR9I2saqk64TZ9L57SFYaXwFbOGUnr6WSHEHcB+YBXweGVFUpyM1AU8dJYQ+uCsJRH+tKOdvZ0xOqMpGqrK60knMia98TQ1AQ+mZXGsf2o8drI5k/998giJjMn6BSEuXVN/wm16DI3TmsP4lElHMYsoZSD308DXgKPA7VLKz1RcKsVJR23ATSnOLX63wemLbbPOowdOrLffPpiashSIf9jRzlEnYuZ1Zy46YddMQ9fYsFDZ8BWzj5JGqKSUfwf+DiCE+Gcp5S8qKpXipMNl6FT77AlYE3HBilqeONTHk4f7eMn6xrLjz/fEMyTSWY5OkV++7Ijy0P4edA3+5ZwlJ5ymUddgfVOVSmFYJrlclt7eTrLZygTLG4v2dm1Sob8rTTlyhcN1BAKlmSPHvCuFED+VUr5WCHEAyJ9ZA8KAUvqK46gLekpS+ksifpZEfBzpS/LcsQHOLjN5ugXsbI9OGP6hFOLpLD/fcgyAF61tnJIJU6vrgyrZySTo7e3E5wsQDE5dmItSMAyd3BTP5J4KSpUrnU7R19d14kofuNn5vE1K+dV8oRDi7aU07EziegXQAVhSyk+O2q8B73Y2lwMRKeWbSpJaMSupC3p4vsQIl5uX13Lk2RYeOdhTttIH2N0eZVG1D/0EJ0z9+rlWBpJZltf6uWwK7PiLqn0Vy40718lm09Ou8OcCbrcH0yw9Q9yY79VSyj3O6tDArRBiIzBhDjshRAD4JvA+KeXNwEYhxBWjqr0e6JNS3i6lvBH4SslSK2YlXpc+YYz9PKcvrsbn0jnUkyjb5TKZyTGQzE46wUqeLUf7efbYAB5D57qzFp/wjNsav5uVdYETamO+oxR++ZR7zUoxpr44vyKl3AqIEo7ZDBySUuYfEA8DV4+q8zqgVgjxHiHE54DyI3EpZh11JZo1vC59qIf/6MHy4vF0O8HXeuKTV/rRVJbfbm0F4B9OazrhDFU+l866piqltBRlc/DgAa6++ioOHz40ovxTn/oEt976uSk/33g2/TcANwDLhBCXOsUaUEq3rBEYLNgecMoKWQaEpZSfEkKcAvxVCLFOSjkiGHp19cQhdcfCMHQikdnX85rLchk+D13p0uyjl69r4qH9PTx9pJ9Xn7N0TNdGQ9cIOtE8szmTFEl8TigD09AJlZmYBeBXz7UST+cQTSEuXz85k0JeLl2Ds5fVTEqOSnCy3l/t7RqGoXPmrfdN2Tm3fPDSouWxWIzPf/7TCLGWVCrFk08+zre+9T9Tdt6xuPbaF3PHHX8F4JFHHiIajfLCF76YVatWYxg6RoFTw+WXX8HDDz84omwsNE0r+Tcfz6b/O+A+4K3At52yHNBaQrsdDGfaAnvwt2NUnQEc05GUco8QIgwsAQ4WVuo/AT/sSCRAX1/5uVkrzVyXK5fKkMxOrPir3Tor6wLs747zkOxg84ri+XmCQS+xmP3S2BVNkSgw6xzpjLKstjwFt6cjymMHenDpGq84rYl4fHLeInm5VtcHySUz9CVPzNw0VZys95dlWVM+oDpWewcOHODIkSN89KM34/f7OPfczXz2s59E1w3e+ta38+53v423ve1drFu3ni984bNs2LCRHTu28YEPfIT3ve+dLF68hDVrBE899QS33nobLS0t/OY3v2TlytVs2fIUH/vYp9i+fRv33PN3li9fyVNPPcHb3vYu2tra+OEPv88555zHnXf+HY/HwxVXvBDLsvjtb3+Dx+Nm166dfOlL/wWAadrX5Pvf/y6GYZDJpDFNk7e8ZeTQqmVZx13bhoYQxRhT6Usp+4F+4KOF5UKIM4EtE1zrR7HfELyOiedC4OtCiFogK6UcAO4GVjpthgEDaJugXcVJQF3QU3IAtM3La9jfHeeRA72cv7xm3B63ZVn0jLLjR1M5UtkcXldpYwnprMn/PdsCwAvXNlBf5uSw0dQHPTRXq4HbqeTJ919c8XOsW7eef/zHV/H+97+HRCLOJZdczkUXXcrDDz9IJFLDaadtAuw3gr1793D++Rfwrne9l/r6ek47bROnn34GL3zhS0ilUvzqV7/gscceZu3a9Zimid8fQMpd3HzzR/nZz/6PSKSGc845j9Wr19DUtIDXv95OSXLZZXZPPs8VV1zFunXruemmf+e+++4mEAgMyfDNb36N973v3/F4PDz22CMn9N1LyZG7GHgPUI9t3jkNOHu8Y6SUccfL53YhRCewVUp5txDiVqAHuAX4AnCrEOIm7Jm+b5BSTl0QFcWMUY7S39gc5nfb2mgZSHK4NzFurz2azpHJHu+33B1N01xiqsG79nTSHc+wIOzl0tUn5q3jdemcUqMU/slIT083GzacxjXX/ANgcc01L+YjH/kYuZztBZNI2BaGmpoavvvdH/L444/xxje+nt/+9k8j2in0o7/ggos477zz2bVrJzU1NUXr5Ts12ezY3jbFfPNdLhevec3rAGhqWlDmtx3VVgl1bgF+jT2g+3/Aq0tpWEp5J3DnqLIPFqz3A/9WsqSKk4ZqnwuXrpUUE8dl6Jy7LMK9e7t55EDvuEq/dwwzTF8yS5NpYujj2z47oynu29sNwKtPP/Ect+sXhtEzpbvKKWYPyWSSr33tK5x55tkkkwlWrVrNxo2b+N73vs2vfvVzDhzYz7333k0kUsOvf/1L1q5dx9lnn4Pfb3cuHn/8MY4ePcrWrc/ypS99hRe84GK+/e1vcPDgfg4fPsT73vfv3HzzZ/nGN77KwoWLGBgYYM2aU9i06Qy+8Y2v4ff7OXr0CFLu4umnn2T//n384Q+/4/777yEej3PppVdwyy2fYffunSQSCd7ylrfzn//5Rerq6gmFipttSkWbaMaXEOIDUsovCSE+JKX8ghDiY05ohmmhs3Nw0lPlTlbb5kwxlXLJ9ijt0Qm9ewHoiqX5/J17cekan3jxKQRGJUoJBr30DSTY2xljrNu1scpDwwQRMb/36CF2tkc5Z2mE15y5qCTZxmJRtY+zVjfM+d9xKplIrra2QyxYsGwaJbIpd3LWrbd+jgsvvIgLL7yoglKVJ1exa9fQECraqynFZfMsIcQyoEEI8XrgspKkUMxrSkmjmKc+6EE0VpE1LZ463F+0Tl8iM6bCB+iOp8edsr67fZCd7VG8Lp2Xrh/tSFYeAbfBCuWPPy/p6upi27bnuPfeu2dalElTinnnK9ieON8Avgh8ddzaCgVQE/Cga1Bq1OPNy2uQHVEeOdjDRatqRwzoWpZF7wQ++TnTfjAUC3+QNU1+t832EbhKNBA+AddKDTilMTgluXIVJx/19fX86EcndxSaCZW+lLIwlPIrhBDnV1AexRzB0DVq/G66S5xAtX5BiLDPRWc0zb6uOKsbgkP7oqksmRKiaXbF0tT43TBKIT+8v4fOaJqGKg8XrSruFloqzdW+E3poKBQzzXiTs8aaqbCRCbx3FAqwvXhKVfqGrnH+shr+Ljt55EDPCKXfEyutjXTWYiCZJewfVsrxdI47ZScA125YgGuCwd7x8Ll0lpc5J0ChmG2M19M3gR8VKb++QrIo5hh1QQ9aZ4xSR+LPX17DXXs62dY6QF8iQ8TvJmdaDJQx6akrlh6h9O/Z20kiY7KmIci6phPLhLWmIXjCHj8KxUwzntJ/r5TyuHg4Qoj9FZRHMYdwGzphn4v+ZGlujdV+Nxubwzx7bICH9vfwslObnAHc0h24EhmTWCpL0OuiN57hwX12XJ+r1zedUFychiqPCpesmDLuuecubr/9Nn73uz9P+7nHm5EbBRBC/MuoXdcA/1RJoRRzh7qgp2SlD3DxqjqePTbAYwd7uUo02JE0y0yy0hVLE/S6+PvuDrKmxaZF4ROKk2/oGivrghNXVEwJ53z5gSlra6zZvc888zTveMdbuO6666mvr6Orq5vGxsahCVCV5vLLr+T222+blnONphTvnTdix+ABWAqo2SiKkqkPetjfXbrP+LLaAMtq/BzqTfDYwR7qg158ZSr9aCrHoe44Tx7uQ9fgJetOzEVzWY1fJTafY5xxxlk0NS3g5S9/BStWrOCLX7yFQCDA6aev4w9/+Bsf+tCN3HjjB9mxYzs/+ckPuf76N/LMM0/z8pe/khe84GI+9amPc8opggMH9nPVVS9m9erVI2L03HLLl4vGy/n4x2+ioaGRurq6GfvupSj9t0op9+Y3hBDvqKA8ijmGz21Q5TWIpnITV3a4eHUdP3ryKA/u7+HlGyY35fyPO9uxsF1By02+XkjQY7BIxdaZVqYj9k6en/70R9TV1REOh7n22n8kHo/z7ne/jS984TaEWEtjYxP3338v119/Ay972bVcd90r+eAHP0omk+a1r72evr5errvuVXz7298fEaOnWLwcKXdz8OB+PvWpzxGNDvLLX/5s2r5nIaUo/ZQQYqmzHgYuBb5eMYkUc466gIdoqvRoqactDBPx2zb5o31J1vjLc5HsiKbY3x3HY2hcJRrKFXcEq+qCKkb+HOa1r72eFStWkE5naG9vZ8mSpXi9Xvbtex4h1o6oO97Y0ugYPb/+9R3HxcuZLfdRKUr/fuAA9ryUAeA7FZVIMeeor/JwqLd0pW/oGucsreFO2cn2tgHWLCwv1sizR+1ZvZsWVZ+QT31d0EMkoHzy5yLPPPM07e1t/PjHP6SmpoZEIoHb7WHLlqf44he/wj/8w0vQNI1TT91AZ2cHP/vZj3nyycf58Ic/xsUXX8qDD97Pz372Yw4e3M/NN3+WQ4cOjYjREwxWHRcv57LLrmD58pXcfvtt+P1+2tvb2LFjO6eeumFav3spsXdeJKX82zTJcxwq9s70UUm5njrcRzxTuonn+c4o33n0MFnT4tVnLaKqxNDJHdEUf9zRjkvX+OfTmzm1OVRy2OVCdA3OWhLBP0ZiF5ifv+OJcDLG3jl8+BBf/OLn+epXvznNUs1s7J17hBBvF0J8TQjxDiGE8ltTlE056QhNyyKTs1jjTNDaeqx4PJ5i5Hv56xeE8LoNOqOTS5DSXO0bV+Er5gf33ns3+/fv4+DBAzMtypRRitL/H+wY+vuxZ+NWPqeYYs5RX1W60h9MZjEtOHWBbdbZ0x4llp7YaawjmuJofxKXrrHBObY/kSVZxhsGgEvXWFpifH7F3OYNb3gTf/rTnSxfvmKmRZkySrHpd0gp35/fEEJ8pXLiKOYqVV4XPpdeUhrFPic7VtjnZnmtn4M9Cba3DnLesppxjxvq5TeFRuTb7YimWVqGn/7SGj+uMt1EFYqThVLu7N5R20cBhBAvmXpxFHOZUnr7WdMklh7umW9srgZAdkTH7bGP6OWPGvgdTGZJpEvr7ftcukp/qJjTlNLTf40Q4l+xE5YvA3qFEFdjT9RaVUHZFHOMhqCXo33jp1EcSGRHxM2vD3pYUuPnSG+CXe1RzlhcXfS4Z48V7+Xn6YimSkqgvqw2oMImK+Y0pSj9vwG3FylXk7QUZRHyTWziKRay4cylEY70JtjZPsiGhSHco0wvPfE0R/uSGLrGqWO4d0ZTuaGYPGMR9Bg0ljH2oJi7RKNRPvGJm2huXsQZZ5zJ5ZdfNdMiTRkTmncce34fdmL0XinlIWf5UKWFU8w9xjPxZHIm8SJmmOZqH41VHlJZE9lxXAxAtrUMACAaguN63LQNjp++cXltYNZMoFFUlmeeeZrNm8/kO9/5RtH927dvo7Gxife//0MlK/yXv/ylQ+uPPvowd9751ymRdaqZsKcvhLgW+G9s235ECPFOKeUfKi6ZYk4ynolnrMBsmqaxsTnMXXu62N46yLqm0FCI48Fklv3dcTQNNiwMj3vuZMakP5GhusgM37DXRV0ZbqWKyjEdAdfysXde8IJLuOSSzbziFa+mq6sTj8fDTTd9nLvu+htS7uKee+5i8eIl/OY3v2TlytVs2fIUH/vYp9i+fRv33PN3li9fyVNPPcHb3vauocle55xzHnfffScej4errnoxt912K7W1dfT0dLN27XpWrVrNDTe8lkcf3cIDD9zHbbfdyu9+92e+971vARq9vT2cfvoZvOhFlRk2LcW880JglZQyLYTwYadPVEpfMSlCPhd+t04ic7yJpz8xdtz8JRE/NX43vYkMezqjrGuyzTjb2wawgFV1AarGMd3kaR9MEfa5juvRL6tVLprzkWAwyIYNG7nssis49dQNXHvtiwG47LIr8Hg8XH75lfzbv72RtWvXY5omfn8AKXdx880f5Wc/+z8ikRrOOec8Vq9eQ1PTAl7/+jcMHf/www/y/PN7efbZLfzv//4c0zS5+OLzeOihJ2lqsmNKXXzxpdx2260AHDp0CK/Xy8tedi0rV1ZuuLQUpX9ISpkGkFImhRCHKyaNYl7QEPRyuG9kWIZU1iRZ5EGQR9M0Tl9Uzb3Pd/HcsQHWNFSRyZns6YgBsHGCXn6eTM6iO5YZYWaq9rlUrPxZxHQGXMsTCtmdCH2MzGoXXHAR5513Prt27aSmZqTrcD6qQb4jkc2O/caax+Vyk81mSaeHJw++613vZXCwn+9+91ssXNjM+973gcl/oXEoRemvEkLciD05axW2B49CMWnqqzzHKf1SsmMtr/VTG3DTE88gOwZJZkxylmW/BZShtDtjKSJ+15AvfilePYq5RT72zi9+8VP279/Hgw/ez/r1p9Le3saWLU9x7713I+UuduzYzgc+8BG+/e1vcPDgfg4fPsT73vfv3HzzZ/nGN77KwoWLGBgYYM2aU9i06Qy+8Y2v4ff7OXr0CFLu4k1veitnnHEWP/zh/9Db28NNN30CgJe+9GXcdtutrF69ZigGzx13/IYVK1axYMFChFhXse9eSuydKuAm7Nm4zwK3FMuoVSlU7J3pYzrlevpI3wh//Oc7Y6TG8Orx+VwkHXv/4d44d+3pwufSMC1I5yyuXt9EU6i88MkRv4tFET/VPhebFhV3A52wDfU7lsXJGHtnJpmx2DuOgv8icDPw5elU+Iq5S2GM+2QmN6bCH82SiJ/6oIdk1iKds2gKectW+AB9iSzxdFb18hXzjgmVvhDiemAHdsydnUXSJyoUZVPoD19OOkXbtj9sv19/AsnOB1NZqn2lWDgVirlDKWEYrgWWSyk3AiuAl1dUIsW8wOc2hhRuKfb8QtIFbwU98clF0QQIely0DIw/Q1ihmGuUovSfKfTeAZ4EEEIsqqRgirlPY8hLMpMjnS192MayLHa2D1sYt7eVFoFzNH6PTpXXxcGeRMmmJYViLlDKu+0GIcSnsL13VgKNjonnGuCfKimcYm7TEPQwmCpPYXdG03TF0nhdOk1VHg73JdlytJ+LVpaXaDo/ppAzLZ7vig2FcVYo5jql9PQXAjlsV80c0Ipt5qmtoFyKeYDL0JnAeew4drQNAiAaqzh3WQ26Bns7Y3THSjfz+Nw6oYKJXN2xNJ3R8UM0KBTTzQ9+8D3e/e63TXm7pfT03yOl3Da6UAhx6pRLo5hXRFPZsrJTRVNZDvbYIRfWNVYR9LpY1xRiR9sgTxzu5cVrG0uKnVMsi9fzXTEifvdxwdwUc5Nnnnmad7zjLVx33fVUV4e5//77+NznbmXBgoUzLdoQl19+JU888fiUtzuh0i+m8J3yHeMdJ4S4EngF0AFYUspPjlHvdcCPgZByB51fdEbThLwGLh1KMavvao9iAStqA0PRMk9fFGZvZ5TWgRRH+hIsrRnfBdPj0ggX8djJ5Cz2dsZYr8w8M850xt55+ctfwYoVK5xwCc/w1FPfQoi17Nu3l4suupQzzjiLz372ZjZuPJ0tW57ihhv+lYGBAe6662+sXLkKKXfzlre8jTe/+Xre8IY385KXXM273/02PvKRj/P4449iGAaZTBrTNHnBCy7hHe94C694xT9x4MA+Nm++kKNHj7BwYTMdHe0sXbqMa655Oe9//3tYv/5UMpnyHBxKpSLdGiFEAPgm8D4p5c3ARiHEFUXqrQPWV0IGxeynM5oCTSMSOD4A2miyueEIm4X2d6/LGIqx/9ihXrITTGapD3rGfBvoiqVpU94884o//OH3fPzjN6HrBh6Ph+3bt5LL5WhoaOTpp5/kV7/6OUuWLOWf//m1fOhDH6WxsYmbb/4o73zne7juuteTyaSRcjfvf/+HOHhwP5qmcfHFl7Js2XK++c2v4fF4CAar2LFjO+vWrWfDhtO4/PIr+fKXb+eMM87iL3/5IwD19Q08/fRT3HvvXQSDQd761ndw5ZWVCec8Zk9fCPE94K2A2/HaKYfN2DF78obSh4GrgbsL2g8AHwT+DXvGr2IeEU1lh+LqR/xuuqLj92qe74qRzpk0VHlorBo5GWtdU4i9nTF64hmeaxngrCWRom24DK1ohM1C9nXHifjdRROxKKaH6Yy9c801/8CSJUu45poXs27depqbF/Ga17wO0zR58MH7RyRE1zQN0yzeqbjyyhfxla98maqqMK973fUAuFwuXvOa1wEMBVgDO85PvuMRDkeG6tx7792YZnn5nCfDeOadnVLKnBDio8Cn8oVCiH+RUv7vBO02AoMF2wNOWSGfBT7lRO8cs6Hqav+kY5wbhk4kMvtmXCq5oKsjSjBoK+8gUJsyiY/heqkBu5xe/ulLIviKmGcuPaWB3zzbwrbWAdYvCheNxbOw2k+oauLZu0fjWc5cUoWuj3/fqd+xPCaSq71dw5imMZUtW56mvb2dP/7x97z73e/lP/7jE/zXf91GKpXkRz/6Pv39/VxxxVWcf/5mPv3pT/Dzn/+EXbt2cOON/86nP/05vvnNr7F8+Up8Ph9XXnkVmqbx6le/hl27dtLQ0ADAv/3bO/jKV75EfX09oVCItrYWDhzYzx13/Jb/9/9uRAjB5s2b+c53voGu66xZI7jsssv54x/v4Nvf/jr9/X0cOLCfI0cOTZiYXdO0kn/zMWPvCCF+DiQZjrkD9v/faVLKs8dr1DHl3CSlvMLZvhFYLKW80dleAnwa2O0c8nngE8CfpZRPFbalYu9MH9Mp1xOHekdk0OqLpznWX9yDpiOe5o/b2gi4DV59evOYyvih/d3s6YyxMOw9blDX0GFNQ9VQHP6JWBzxsbIuOG4d9TuWh4q9Ux6Vir0zXk//jcAZwFuAHxaUX1+CDI8Cy4QQXsfEcyHwdSFELZCVUh4BbshXFkJ8HrhNDeTODwaT2eNSJoZ9btoGUxS7x7c6+W/XNY3f+z57SYRDvQlaB1I83xVjTcNwiIaagKdkhQ9wtC9J2Ocu6umjUJzMjPkuJaVMSCkfwXbZvD+/AO+dqFEpZRx4O3C7EOIzwFYp5d3AhynIrSuEaBBC/Iez+UE1y3d+0FXEp17Xi9vb+xMZDvckMDQN0Th+nB2f2+DcpREAHj/UOzRTV9OgroTB4tHs6YiSyFTexqpQTCel+Ok3CyH+h2Ezz78CcqKDpJR3AneOKvvgqO1O4DPOopgnjDURqjbgpic2ckA3PxlrVX2gpMHV1fVBDvbEOdKX5OEDPVx1SgM1Ac9Q7PxyyJoWO1oHOX1xNa4y3hIUk8eyLJWnuEwmCo8/mlL+Ez4IvA97Ru6/Y/fWFYpJUcy0k8frMgh4hm/JVNbk+S47M1apYRI0TePCFbV4DI2jfUme74qdUO7beCbH7vbBsv+xFOXjcnmIxQbUtS6TTCaNrpceLbaUmrullE84648JIS6blGQKBXbWqvGoC3iIp20P4T2dUbKmxeIyM2MFPC7OX17LA/u6efxQLxetqsPrmrzi74ln2NcVZ3XD+AO7ihOjpqaB3t5OotG+aT2vpmmz8kFTulwa4XDpUXFKUfqnCCHOAA5gp0tcXXLrCsUouqLjx8gJ+1y4DY1U1mSXY9rZuLj8zFar6gIc6olzqDfBj586yjsvWo5rjPynpdAykMTr0llSoxKoVwrDcFFfP/1hEE5Wb6fJUsp/wZeB/waOAv8F3DrlUijmBeOZdobQNGoDHg73Joimc4R9LpbVlq9oNU3jSlFPjd/N4d4Ef97ZMUmphznQE1czdhUnPaXE3tkNXDANsijmOBOZdvLUBFzsdHr565tCkx7YWxzx8/pzFvPfDx7g/ue7WV0fPOHYOns7Yxi6NiLdo0JxMqFCCiqmjc4JTDt5WgdStA2mcBsaayZpR3e7NEJeF8trA7xkfRMAP3366AmHULaA3e1RulQoZsVJilL6imlhIJkpOUPVg/u6ARANVZMOdVwXGA6sdunqOk5dECKRMfmfxw6fsO+9hR3xU5l6FCcjSukrpoVSe/kDyQzPHB1AA85bFpnUuQwdagomY+maxmvPWsSCsJeOaJofP3UU8wS9NSxgZ+sALf1K8StOLkpS+kKI9UKIPwkh/i6EeGWlhVLMPSby2snzyIFecpbFhoUhVk3StFMT8KCPGgfwuQ3edN5SAh6D3e1Rfr+t7YTd9CzLjv55sHv2eX4oFGMxptIXQhQmHb0OeKWU8oXAaRWXSjGn6E9kSJUQOCqTM3n0YA+A41tvFE14Mh6aZs/sLUZd0MMN5y7B0DUe2t/D3Xu6ymp7LA73JdjVPkjOnH2+3grFaMbr6X9LCHGds54GXi2EeCl2zlyFomQ6S8xf+8zRfqKpHM3VPlbW2WFiG6rKm1Rl+/mPfVuvqg/yurMWoQF/2dXBowd6ymp/LDqjaZ5r6SepYvUoZjnjBVx7FRAWQvwE+KlTdw3w8WmSTTEHsCyrJE8Xy7J4cL/Ty19ZOzQI63MbRCZIfFJIKVExNy2q5pWb7L7Lr59r5anDfSW3Px7RVI4tR/uLBpRTKGYL49r0pZTfwo6380kgIKX8Lyll+7RIppgT9CeypHMTmz32d8dp6U9S5R1Of5inKeyjFE/9oNcoOePV5hW1vGRdIxbw8y3HePxgb0nHTUTWtNjZNsjznTFl7lHMSsaz6b9KCHEndq7bzwJdQohfi/HSXCkUo+go0Z/9/udtN80LltceZ57xug0i/olt++UGVrtSNPDS9bbi/+WzLTy8f2pMPWCHbdhytI/+RGWSWysUk2W8nv7pUsqrgFcDr5VS/hI7ocqN0yKZ4qTHtKySTB1dsTQ72wYxdI0LVhQPHNUQ8jJedGOvSyfkLW/QF+CKUxq4doM9ees3W1v5y872KQu+lciYPNcywJ6OKJlZmJlJMT8Z77+kwRm4rQUyAFLKHuxE5grFhPTGM2RLMHE8tK8bCzhzcTWhMbx13IZOfdBDxxiun3XB8pOk5LlkdT0el85vnmvlrj1ddMczvOaM5knF4C9G22CKrliapTV+mqt9x7mTKhTTyXh39cewB24N4JbpEUcxlygl5EEik+MJZyD14lV149atC3pwG8crTJdRPOtWOWxeXsubzl+K16XzzNF+vv7QQfqm0DSTNS32d8d56nAf7YOpWRnKVzE/GLOnL6XswI6qqVCUTc606I5NrDSfONRLKmuyuj5Ic7Vv3Lq6rrEg7ONIb2JEeW3APSW953VNId510Qq++9ghDvUmuO3efbzu7MUTpmksh2TWRHZEOdyrszjipynkVT1/xbSiwjAoKkJ3LE1ugt5szrR4cJ89eDpRLz9P2Oci5B320NG1kSEXTpTmah83XrqKUxqDxNI5vvPIIf64o23KbfKJjMnezhiPH+rlYHdc+fcrpg2l9BUVoRSvnR2tg/QmMtQFPaxbUHpvemHYRz4fSrXffULJUYpR5XXxls3LeOHaBgDu3dvNf963n8Oj3jCmgkzO4nBfgicP97GtZYCOwZRy9VRUlPLdHRSKCcjkTHrjE5t2Hthvu2letLK2LBOH26XTVOWlbTB1Qvlvx0PXNF60thHRWMXPtxyjfTDFVx/YzwtW1vFC0YDfU9p8gFKxgN5Eht5EBkPTqAm4qQ96qA24p2xAWaEApfQVFaAjmmKivuqR3gQHuuP4XDrnTiKaZm3QDqrmdVVWIS6vDfD+y1bxl10dPPB8Nw/s62bLkT5esr6Ry9YtqMg5c46ra1csjYZt0qr2u4n43YR9LjUGoDghlNJXTDkdgxP75j/gxMw/b3kNXtfkes2bV9TwfFeMTAkzfk8Et6Fz7YYFnLm4mt9va2N/d5xfPdvKQwd6uWJNPZsWhSumiC2gP5mlP5nlcG8CXYOQ10WV10XI6yLkc+Fz6ZPOLqaYfyilr5hS4ukcg6nsuHV642mePdaPrsELVhafjDURYa+Lhiovhq6xvXVwUm2Uy+KIn3e8YDnPHhvgTzvaae1P8uOnjnKn9HL5mnpOXxSuuCnGtIYfAnkMTSPgMezFbZAxDNLpLD6XgTHejDbFvEQpfcWU0jE48QDuA/u6MS04Y3E1tYHJ2eQX19jJ0msDHpZG/Bzum/pB1mJomsYZi6s5rTnE1vYYf97aSvtgip9tOcYfdrSxeXktm5fXnPC8gXLIWRaDqezQw7YjlSPm5CP2GBpel4HXpeNx6XgN+9NjaLgN3Vk0ZTKaRyilr5gyLMuifQKlH09neexgHwCXrS7NTXM0frdOXYGb5rJaP7F0lu4SBo+nCpeuc9HqejY2BdlypJ8H9/fQ0p/kTtnJXbKTNQ1BzloSYcPCUMlB4CpBOmeRzmWZ6FlsaBouXcNlOJ+6hqFruHQdXbf3G06ZrmkYmj1vQtc0dA3n017XRn+CMj/NIpTSV0wZfSUkS3n4QC/pnMkpjUEWRfyTOs+SiH+EEtE0DdEUYmuLHY9/OnHpOucuq+GcpREOdMd5aH8P29sG2dMZY09nDLehsbaxinULQqxrqiLsm743gHLIWRa5nEWlLp8GVFVFicfTzkPAeRhgr9gPB/s3zf+0mlOW/6ULnxv537/oPooUFpaP2hWKZRiMJseN5KpNFOf1BJ5pYx263l8ZzzSl9BVTxkS9/EzO5CFnAPfyNfWTOofX0GkMeY8rd+kaGxaEefZYP8kSE7BPJZqmsbI+yMr6IPF0ludaBthypJ/93XG2tQ6yzRl3WBLxs6YhyMq6AMvrAvhn8C1gOrGwxyOOn4Mw83MS4pZGbBbmQFiVyVGJu0MpfcWUkM2ZE0bUfPJwH9F0jsURH6vrJ5f/dlFk7IBlHpfOxuYwz7UMkJoBxZ8n4HE5tv1aeuMZdrUPsrNtkL2dMY70JTjSl+CevXYPb2HYx7JaP4siPhZV+1kY9o6b+UuhOFGU0ldMCe3RFONNJDUti/v22jlpL19TPykbr8uJvTMePrfBpuYwW1sGZqTHP5qagJsLVtRywYpa0lmTfd0x9nfF2d8d40hvkpYBe8mja9BY5aUx5KWhykN90ENDlb0e9BjKNq44YZTSV0wJbQPjm3a2tgzQHbdDLpzWHJ7UOZrDPlwluCD63AabFlWzvXWAWHr2xLTxuHTWNYVY1xQCbHPXoZ4ER/sSHOtPcqw/ScdgijZnGY3b0Kj22ZO0qv0u+9PnJuR1EfTaLptBj4GvQrZgxdygYkpfCHEl8AqgA7CklJ8ctf9DwAKgFTgb+LiUcnel5FFUjsFkdlzlaloWd8lOAC5dXTcp90BD01gUGb+XX4jXpbNpUTWyI0r3LLTXgj3pa3VDkNUNw6audNakfTBFZzRFZzRNZyxNl7OezJpDM3UnwufSCXoM/B7bXdNr6Hjdznp+MXS8bh2PoTseO/qQ947bcLb1Qq8e3fbgcTx3lFfOyUlFlL4QIoCdZvFUKWXKSbN4hZTy7oJqVcCNUkpLCPHPwBeBayohj6KytBaYJ4qxo3WQ1oEU1X4X5y6NTOocCyZh63bpGqcuCHG4N8GhnvgsGDKcGI9LZ0mNnyU1x3s2JTM5+pNZ+hIZ+hIZ+hMZ+hNZomn7oRtP54ilssQzOZJZ0zZvVdiNdbS7pq5p6PrxZXl3T6xhzx1GuXMOl2nDdRj24smvF5aNaMfZUegdNBGapmEYOrlxvM6KtqOVUKfIuSasU7D+9z1dvOHsRaxpmLrQ3lC5nv5m4JCUMv+O+jBwNTCk9KWUHyuorwPRCsmiqCDZnEnnGNmswPbd/7vTy798Tf2kZqzqmj0bdrIsrfET8bvY0xEjfhKHMPa57cTvTUW8lwrxBzx09yWIpbMkMiapbI5U1hxjyZHJWWRNi0zOJGva69mcRdY0C9aH91uW/faW98gxLYvZ4IUz9+jnpWsbWNMwta1WSuk3AoVz4wecsuMQQniANwDvLLa/uto/6VdIw9CJRAKTOraSzCW5jvTG8Y0zq/bZI3209Cep9ru5fP2CSXmmLKkN0jRJb588EWBxY5hDPXEO9cSnJHyxoWsEg+Mr4JnA0DUaagJMsa44DtMafgDkTAvTsjBNZ9uyME1r6KGQsyw0GIqTZDkPjXzKhTG3sf9Y2OfK/2qWNXp7uE5+/2jGchY1NG0494NVvM7IQmviOuMfMsZxI0uW1wc5c1U9kQmcF8qlUkq/AwgVbIedshE4Cv8bwEellPuKNdTfP/np9ZFIgL6++KSPrxRzSS55pG/M3rNlWfzhuRbAnn2bTmYo17qua7BkVd2UXa8al0ag1s+hnjgdE3gcTUQw6B0KdzCbmEm5NOz8qkZ+QxteUderPC4UjRjZ3KTv/YaGUNHySjkEPwosE0Lku0EXAn8SQtQKIcIwZPf/FnCblPJpIcQrKySLokL0xtPjmkt2tUc51p8k5HVx/vKaSZ1jQciHd4onMHldOqc0VnH2kgiLI6V5BCkUc4WK9PSllHEhxNuB24UQncBWKeXdQohbgR7sROs/BjYAK4QQAEHg15WQR1EZWvrHHsA1LYu/7LRf7i5fUz8ps46uUXRAc6rwuQ1W1gVZXhugO5amI5qmL56ZMM2jQnEyUzGXTSnlncCdo8o+WLD+ikqdW1F5EpncuAHOnj3aT8tAkmq/i80rJtnLD/sqniQFbO8SewKUl5xp0Z+0vWMGHM8Ylb1QMZdQk7MUk+JY39i9/Kxp8tdddi//RWsbJ9XLNzSNpSfgsTNZDF2jNuAZCvlsWRbxjO0OGU87rpCZHOmcqcxCipMSpfQVZZPJmeMGV3vsYC/d8QxNIS9nL4lM6hyLqn14pqGXPxGaphH0uAh6jv9XiUQC9PTGhlwac6bjtWLZXitWwWeeQm+T4XMUrE8gy3h18u2Ew34GBoo7QJT6mKrEy814cs0ks1WuKp+bRHTqXYyV0leUTUt/cky7dyqb407HL/8l6xonlbnJpWssLmP27Uyiaxoel8ZsCnwQCfvwmjMfd2g0Sq7y8Lp0KvEomvmulOKkImdaIwKEjeaePV1EUzmW1fjZsLC4y9hELK3xVzztoEIxX1H/WYqyaB1IjpmIvCuW5t7n7Xj51562YFKT6rwunebqk6OXr1CcjCilrygZ07I4Ok4u2ju2tZEzLc5aUs3y2snNOF5RG1D5WhWKCqKUvqJkWvqTpMfo5cuOKDvaBvG6dK5e3zSp9sNeV9GsWAqFYupQSl9REjlz7F5+1jT53dZWAK4UDVT7J5cHdtUJxtdRKBQTo5S+oiSOjdPLv3tPFx3RNPVBDxevrJ1U+wtCXkI+5UymUFQapfQVE5LJmRwZo5ffNpDkbmmnQfyn05sn5XXj0jVW1M2+qKMKxVxEKX3FhBzqSRQNRWxaFr94poWcZXH+8poRGaDKYWVdQCUDVyimCfWfphiXWDo7ZmasB57v5nBvgrDPxctOndzgbbXPNWGyc4VCMXUoI6piXJ7vjBWdkn+0L8GfnSiar9rUjH+s8MeWRSpnksiYpLMm6ZxphyswLdDglIYqnj7Sh6bZ8XYKU+u5nM9aE+LRJC5Dx61reFw6bievq0KhKA+l9BVj0j6Yoj+ZPa48lTX5yVNHyVkWF6yo4dRRM29zpsVgKuskTM8yVvrRBWGvnY5vnKTqAN1ps2iSC0PX8Ll0fG6DgNsg4DEIOotK2K1QFEcpfUVRMjmT/d2xovvu2NZGRzRNU8jLNacuGCpPpHP0xNMMJCcORxzwGNQFTyxiTc60iKVzxNI5ugvKDU2jymtQ7XNT7XcR9rknFQNIoZiLKKWvKMrzXbGi4RYeO9jDY4d6cekarz97MR6XTiyVpSOaIp4uLWiVoVM0oFo+16qm2YHMJjszN2dZ9Cez9ltKn52Mpdrnpjbopi7gwTfFmbgUipMJpfQVx9ExmKIzenw22wPdcX7zXBsAr9y0kNqAm0M9caKp0sO/5r2AHj3QS0c0RVcsTU/MTruYzJgjxg/chkbI66La7ybkddEU8tIU8rIg5KUx5C25925a0JvI0JvIsI84VV6D+qCXhirP2GMRCsUcRSl9xQiSmRzPdx1v1umNp/nhE0fIWRYvWFnDkoif/V3xCeOuW06v+0hvgiN9CTqjKcaY4wXYPvuWZffWMzmLnniGniIZurwunaU1fpbXBlhZH2BlbaDkOQLRVI5oKs7Bnjghr4uGKg+NVd5ZEb9foag0SukrhjBNi13tUbKjDPLRVJZvPXKIwVSW5bV+1jZWFVXEhcTT9sPj+c4ofQWDwRr27NultX6awz7qqzzUBTxUeQ28LmNE7z2VzTGYypFBo6UnRvtgivbBFC39SXriGfZ2xtjbGQMJHkNjdX2QtU0h1jZVlTxeMJjKMpjKcqA7Tk3ATVPIS13Qo4K+KeYsSukrhpDtgwymRnrrJDM5vvPoITqjaeqCbi5cUYdpFVeIlmXRNpBie9sAR/uSQ28BXpfOsho/5yyNsK4phN9TmknF67IfBMGgl4XBkfF8BpIZDvYkONgdZ09nlNaBFDvbo+xsjwLQXO1jY3OYjc1hmkoI4mbB0FuFS9dorLJNSSo0hGKuoe5oBWD73bcnR9rmk5kc33vsMEf7koS8Lq46pbFoonLLsjjcm2BrywCdMXssQNNgWcTPmoYgy2sDrG4ITums27DPzcZmNxubwwD0JzLs7oiyuz2K7IjS0p+kpT/JX3d10BTysrE5zKbmMAvC3gndObNOopiWgSRBj0FTyKvMP4o5g1L6CjqjKfZ3xwkGh3vEsXSWbz98iKP9SQJugxetbSQwqodumhb7umNsax2kL2Gbe7wunVMXhBCNVfjdtrlmea1/QoXvNXT8HgOvodsmHs0e9E1nTVweg1hx79Ehqv1uzltWw3nLasjkTPZ2xtjaMsD21kHaB1PcKTu5U3bSUOVhk/MG0Fztm/ABEEvn2N8d50B3nNqgh6YqD7XK/KM4iVFKf57THUsjO6IjyrqiKb7z6GG6YmlCXhcvXts4wsyRzZns6YyxrXWAmDOxKugx2LAwzCkFPfq8wi/mIuk2NGoDHmoDbqp97nF70ZFIgK4eH4PJDH2JLD3x9NB5i+E2dNYvCLF+QYicafF8V4ytxwbY1jpAZzTNXXu6uGtPF/VBD5sW2Q+ARRM8ACznWnXH0rgNjYYqL6d4JxdCWqGYSTRrjATXs4XOzsFJCxiJBOjri0+lOFPCbJGrK5pid0d0aCJVMOjlyX1d/Pipo6SyJhG/mxevbSDgsRV+Kmuyq32QnW2DJLO2T361z8XG5jAr64IjBmHdhsayWj9e17DC1zWoDXhYEPZS43eXPGu22PVKZHJ0RtN0DKaIZ0pzGc05byZbjw2wrWWAaMGDoy7osU1Ai8IsLuENAOzrZaYzNFbZLqSzxf1zttxfo1FylceJytXQECp6EyulPwPMBrla+pPs6xqOq5MzLf62p5N7dndiAUsifi5ZVYfHpRNPZ9nRNsju9igZ5wlR7yjJZTX+4xSk362zpGbYpOMxNBaGfSwM+yZlF5/oeg0kM7QNpOiMpYtGAy2GaVns74rzXEs/W1sGRsw1qAu42bgozMbmapZExn4ABIPeEeEh8u6f9cGZnQA2G+6vYii5ykMp/UkwV3/MEyGv7FoKImce6Y3zsy3HaB+0B2E3Noc5a3E1fYks29sG2NcVG3obaA7bXjELxxgQrQm4aajykMlZ6Jq9HfS4yJoWmZxJ1rQwLWuoPQ0n0JoOLl3HbWi4dR23S8Nj6PhcOg11VSRjKTyGNm4PPGtadEZTtA4ky5owVvgA2NYy0oOpNuBmw8IQa5tCx4WAHq30Cwl5XdQHPdTPwAQwdd+Xx1yVSyn9WcRMyZXI5NjdHh1SarF0hj9u7+CpI32YFgS9Bi9YUYuuaWxvHeBI3/CDYXmNn9OawzRUHe/+mM7lyJkWIa8LXddwaRp1VV7C47g7ZnImg8ksyaxJMpOzP7MmlmVReEvqukZ10IOZzeFz6VR5XdQFPNQG3fjdLjwuzXbtNHQ8Ln3owTCYtENCd0bT5Mq4x03L4kB3nOeODbC1ZWDEA2D0XICljaExlX4hAbcdZ6g24Cbsc1U8GJy678tjrso1ltJXA7nzAMuyONaftJOhWBaJdI5793bx6MHeIXv4yroADSEvjx/qG/LEMTSNNQ1BNiwMEfYND1rmLJNkxiSRyZHKmvjdBhGfG7/HRX3QPTQGkM6atA0mHffJFD3xNP2JLP3JzLgDsaWgAVVeF1Vegyqvi2qfm5qAvTRW2eMGC6rsiVaxdJa+RIZU1sRt2G8TY3nf6JrGqvogq+qDvHzjAg72xNnVZruCtgwkR8wFaAh5WVHjZ1V9kJX1AWoDxSeExTM54n32jGSXrhHxO7L63SoOkGLaUT39GWA65eqNp9nfHSeWztEbT3H/8z08ebhvaCA27HMR8hi0DQ73iP1uHdFYZU+kcpRSOpcjkTFJpHNkTBPL0gh4DGoCLuoDHtyGTmcsbSv4AVvRd0XTY4Zp0DUI+VwE3PYELJ9bx+vSMRxlrDl/cqZFzoJ4OkvaeRuIpXIlDd5qzverCXioC9q9bJ/LjsMfciJw+lwGbkPHY+h4HJPSWO6lhXMB9nREh65hnhq/m5X1AZbXBlgc8dMc9k4YGsLn0qn2u6n22TGGpsIUpO778pircinzzixiOuTqiac50pugI5piW8sATxzuGxErx+vSyZkmhXprcY2fU+qDLI340XSGevOJTG5ogFQDsiaYlp0YpSeepmUgRbxIz13XoDHkpTnso7naR33QQ8RvK9sqr6tob1vXbNk8ho7LSZRSU+0nGk2ioWFhjwekMzl6Ehl642m64xm6orY7ZW/cDqzWE8/Qn8iMGxvI0CDodRFylirnM+xzURtwE/K5hmTJJ27Jm49ypkVvOsf2I33s645zoDtGImOOal9jYdjL4ho/iyP2QHZTlXfcGckew34ghYdkMsrOOzyf7/vJMFflUuadeUAqa9I+mGTL0X62twywoy3Kge5Y0QBnKUfb1/jdrKgLsLIuSHXIQ380RVc8TTKTJZkxiaZzxNO24o+mcvQnM0Vj5fvdBs3Vwwq+udrHgtDYPV2fSyfoMQh4XPjdOn63gd9tFPXuiUQC9Hkn7gHbbwLDbyQDqSytzptHVzRNj/OA6Iml6YlnGExlGUjaSzE8hjb0IMg/FMI+F5GAm/qgm9qwnw0LQ5yxuBqXrtEdT3OoxzbjHO1L0jGY4mh/kqP9I9NNhn2uEdFC64Ie6gJuIgE3oA/NBxi+tjpBj8tOEOO1P30uXSWKUUwKpfRPYrKmxZ6OQba2DLDTsTsf7U8UjYOfRwMaqjwsivhYHPbhcun0xTI83xUl1W4xmLDt7YOp7HGB1/LH1wc9NFf7WOQo94VhHxF/8QFKQ9McZWXb3vOKvhKpDj0uu0demHJ3Y3MYy7Jss1A6R8J5iMUzOXoTGbqiKXriGbpj9kOhw3lj6EtkSI8T5RPsRDB+l/PA8tgPraDHoC7gYUVtgIDbIJHN0RvP2GGko3b7+QfN3s6R04w1IOx3BqoDbiJ++20j7HU+ncVt6OgaQw/KgNs2j/nc9sNgtr+9K2aWiil9IcSVwCuADsCSUn5y1H4f8CXgGLAGuEVKuadS8pxMWJadbrA3nqE7nqG1P0GbE13yWF+StsEUvfFMyZOSQl5bOXhdOrqmkc6ZHOxOsLMtOtTjL0bAbbAw7GVB2MfCsNfuvYd9RePvaNhKKOCxlyqnVzobeqSapg0pSIIj96WzwyasZMZ+U7DNWvb1bx1I0TqYGno49CUyTmjmrP3wSOeA8SOOajgmK5dOxO/GbWhogIk9uzmbs4a8l/oTWfoTWfZ3j91e3ovJftjoQ6ki89+xJuTFC0OD3IVmq5DXRcjnJuDWyzYbKeYGFVH6QogA8E3gVCllSgjxayHEFVLKuwuqvRc4LKW8VQhxGvA94KKpOP9gIs13HjuM2+MimbTtupbjG24Blgkm9obpuAiaMKIeOOUWWEPrzjFOvZzpZHuyLHI559OyMJ3BR/vTbi9nWkNtm9jKJmea5EyGjjELZJxKBlN2iOJiGBrOQKeHpmofEa/BwrCPBWEvIe9w793QNNyGNqS8PIZu9y5dtqLxufWTMh5N/u2g2n98SAXTskhlTVIZk1TOdB4QdkauvniG/ozJka4YPYk0g8lhM1giY79JJDI5UhmTjDms1KcCu63jk9yUgwZOAnpGJKM3tOGk9PnE9Pa6jsuw13VNG5pbYQxt2w9X3WlDH8p+ZtfRNA2/10U2kxsq1/LHOKP2dsgl+4E4dCvZYZjsus5Gwa6hNU136jFcj6HjRh47+jb1+z0kEmnnzCOvUbGN8e7y8f4FNG30GUbtH7XtD3g5oynIooh/nKPKp1I9/c3AISll3on5YeBqoFDpXw3cBCCl3CaE2CSECEspB0705M93xfnZlpYTbeakRdPsXnre3BDxu4eWGr9tO87Hjq8PeobMBZFIgOhgEkNjSAnkJ0ydjAr9RNEL3xCKUDjQlsnZ5qN4Kkcia5uRklmTVNZ+GAwks07sftt8Fk+bTrYw2+01kzPJ5CzSzufIbXt9KjsDFrZ50H4OKXPQbOVtFyzjzZuXTWmblVL6jcBgwfaAU1ZKnRFKv7r6+Gn+E7EwkUUf68msjfM0L9we53hteBUAHW2o1zLUK8n3ZqCgZ+P0igwNt273nN2GgbtggpHXpeN1es5Bj/1q7nPbvVGvyxjqmXqdxePYlENeF1U+2/Yb9BiTMqkYhk4ufHzu2pnGMHQikcBMi3Eco+VqKONYy7Ls2cmm8/aXf8uz7Dc++61zWB1b5vDDIJUzyWbtz0zWIu3Mdk5nTdtcZFkkMzmyWZOc85aZzz9sv3WapHP2XAv7GIuMaT9c8uvZnC1fNjeqjZxlv/lCgbzDb8zW0HrBG3R+H/nvaK+Tf/POv1kPXZzjH0PWyD8j9lujCwprlfjmPG6dEp+Jk390jv1AX9xQNeX3fqWUfgcQKtgOO2Xl1qG/P1H2yRf4XTx+48Vz1hVrTCyTbCLNJC4ZMHdd1ypFJeXKmyoKcQF+QwPDgHHcPufj9ToRZrtck5WtoSFUtLxSIzmPAsuEEPk5+xcCfxJC1Aohwk7Zn7DNQDg2/eemwrSjUCgUirGpiNKXUsaBtwO3CyE+A2x1BnE/DLzDqfZf2A+G/wDeD7y5ErIoFAqFYpiKuWxKKe8E7hxV9sGC9QTwzkqdX6FQKBTHoxx1FQqFYh6hlL5CoVDMI5TSVygUinmEUvoKhUIxj1BKX6FQKOYRsz6evkKhUCimDtXTVygUinmEUvoKhUIxj1BKX6FQKOYRczJz1kQJXKZRjlXAZ4AtwGKgW0r5KSHEzcClBVU/68xgnk7ZHgPyefxyUsorhBC1wC3AfuzENjdJKdunWa7l2CG4jzhFYWArcJBpvmZCiAXYv98mKeU5TtmYyX+EEK8HzgBywD4p5bemUa4PAQuAVuBs4ONSyt3OvoPY1w/gmJTyddMo1w3A2xi+174npfyRs28mr9f3gFUF1U4DzpJSHpyO6zWObhjzf1AI8e/Y/w81wN+llHdM5txzTumXmMBluqgFfi6l/L0j204hxJ8ApJSXzoA8hfxVSnnzqLLPAXdJKX8phLgGW7ldP81yDQL/JqW8C8B5QN4FXDkD1+wFwO+B0wvK3kuR5D9CiMXAB4AzpJSWEOJJIcQ9Usq90yRXFXCjc+5/Br4IXOPs+0GR37oSFJML4DVSyoOFBbPgev1dSvkLR5Yw9jXKyzgd12ss3fAWivwPCiHOAy6TUr5UCOECdgkh7pdS9pd74jmn9Cktgcu0IKV8clSRDsQAhBAfBVKAAXzVCVI3nZzm9A79wJNSyj9hX6fPOvsfBn44zTIhpezGVvI4UVrPllLeLIS4crqvmZTy/4QQl44qLpr8B3gR8LSUMu8O9yjwEmDKlVgxuaSUHyvY1IFowfZFQogPYocy/4uU8pGplmksuRzeJYRoAwLA16SUPcz89fpFweabgP8p2K749RpHN4z1P/gy7GuElDIrhNgFXAKU3dufizb9UhK4TDtCiH8E/ua8cv8K+IqU8kvYsn51BkT6gpTyC8CngZuEEBcz8toNADVOr2KmuA74ubM+G64ZjH1/zYr7TgjhAd4A/EdB8UeklLcCnwf+RwixehpFuh/7XvsS8BT27wiz53rp2A+gPxUUT+v1GqUbxvofnLLrNReVfknJWaYTIcRlwGXA+wCklDuklDFn9z3A5dMtk5TyCeczBzzoyFd47cJAr5QyO92yFfBPwC9gdlwzh7Hurxm/7xyF/w3go1LKffnygt86DjyLnd9iWpBSHpBSdjqb9wCXCCEMZsH1crgW+FPBG8e0Xq/RuoGx/wen7HrNRaVfNIHLTAkjhLgauyfx/4AFQojNQogvFlRZA+wrenDlZForhCjMX5CXYSixDTN/3S4FHpVSZpztGb1mBYyV/OdvwFlCiHzCq83AX6ZLKGcs61vAbVLKp4UQr3TKrxBCvLig6mqm8doJIT5f8La4BjjodDRm9HoV8AbgB/mN6bxexXQDY/8PFt53bmAd8MBkzjsnZ+QKIa4CXgV0ApkZ9N45C/v19imnKAj8NyCw7Zsd2F4DH897gEyTXM3A14BnsHsMbuBGIAJ8ATiE7dnw4en23imQ8WfAu6WUXc7255nmayaEuAT4F+DF2D3oLzu7voTtJbMa+Nwo752zsb1R9lTQG6WYXD8BNgAtTrWglPIc58F0M/A00Ay0SCk/N41yvdWR6wD27/ZfUsrHnPozdr2klAkhxOnA66SU/15Qd1qu1zi64Q7G+B90vHdqnOUvk/XemZNKX6FQKBTFmYvmHYVCoVCMgVL6CoVCMY9QSl+hUCjmEUrpKxQKxTxCKX2F4gRwfM5P+nMo5g9zMQyDYp4hhHgQeByoww609x1nVz32rMpPSilfU4HzXoDtS/3FUeW3AudOYaygtwshdkgp752i9hTzGNXTV8wF/kdK+QFsv/BeKeUHnO37pZQSO5zDlCKEiGD7699WZPfXp/h0/w18QgjRMMXtKuYhqqevOOmRUn5/rHIhxHuwJ54tF0K8Ebvn/2VgI/abwPexZ0WuAV4mpRwQQpwKfAjYBqzFDuO8f1TzrwKecGaXIoQ4E/gk8ASQyVcSQqzEfjA8gj056cvAUeCngBc7iuli4HbsqJOXYU9WAvBIKf/DiUL5MPbD6/ZJXSSFwkH19BVzGinl7QXr3wd2A1uklNdjR+wMSSnfjD07+Sqn6neBb0opvwj8iOGZuIWcih1TP8+3gM9IKT/N8CxLgDTwKSeA139ix8XpAt4B1EopDzvt/FpKeR/2DNZfSSk/Bfy5oJ1jzjkVihNC9fQV85F8LJW+gvVehgNabQRe6EQe9TMyTHEeL1AYjO5UhsMCF74VZIDXCCFegh3yogFASvm8EOKgsBP+XMJwD/464HNCiCan7JGCdvzlfU2F4nhUT1+hOJ7ngN9IKW/BTixTLPDcEexEGHl2Aqc46ysLyj8MRKWUn8VOuFLIf2GbkQIFkShDUsp/BF6O/WaQpxY4XP5XUShGopS+Yk4ghPBjm0aqhRBvKih/h1N2nROIbxlwg2OD34idlegU4GLgGqeH/WbgPUKIDwO3YtvgR/M74LyC7bdhD7Z+GrgSO9Lr1cCvgauEEJ/Fzoq0TAhxBYCTHWwh8L8F7dzgJIt5DyNzBpzvtKVQnBAq4JpCMUmEEJ8Ans2nvCvzWC+2vf+/pZTvmKDuVdip8m6anKQKxTCqp69QTBInZPeBSR7+fWzzzi8mqgi0K4WvmCpUT1+hUCjmEaqnr1AoFPMIpfQVCoViHqGUvkKhUMwjlNJXKBSKeYRS+gqFQjGPUEpfoVAo5hH/HzR7KB8eqK8JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for mu, sigma, label in zip(expected.T, std.T, labels):\n",
    "    pyplot.fill_between(\n",
    "        time_span, mu-sigma, mu+sigma, alpha=0.3)\n",
    "    pyplot.plot(time_span, mu, label=label)\n",
    "    \n",
    "pyplot.xlabel(\"Time (days)\")\n",
    "pyplot.ylabel(\"% of population\")\n",
    "pyplot.title('Stochastic SEIR model')\n",
    "pyplot.legend()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}