SIMC_SIR_Exercise_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "533bea75",
   "metadata": {},
   "source": [
    "# Exercise 3: Model Calibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a9539a4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load libraries\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Use scipy.optimize.minimize to mnimizes sum of squared errors function\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "384c6e0a",
   "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>Week_of_year</th>\n",
       "      <th>Deaths</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>39</td>\n",
       "      <td>31</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>40</td>\n",
       "      <td>18</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>41</td>\n",
       "      <td>12</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>42</td>\n",
       "      <td>9</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>43</td>\n",
       "      <td>13</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Week_of_year  Deaths  time\n",
       "0            39      31     0\n",
       "1            40      18     1\n",
       "2            41      12     2\n",
       "3            42       9     3\n",
       "4            43      13     4"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load plague Data\n",
    "plague_dat = pd.read_csv('plague_data.csv')\n",
    "plague_dat['time'] = plague_dat.index.values\n",
    "\n",
    "plague_dat.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d1a05b6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function for numerically integrating the SIR model\n",
    "def sir_model(state, time, parameters):\n",
    "    # Unpack Parameters\n",
    "    S, I, R = state\n",
    "    beta, gamma = parameters\n",
    "    \n",
    "    # Calculate total population size\n",
    "    N = S + I + R\n",
    "    \n",
    "    # Calculate force of infection\n",
    "    lambd = beta * I / N\n",
    "    \n",
    "    # Calculate derivatives\n",
    "    dS = -lambd * S\n",
    "    dI = lambd * S - gamma * I\n",
    "    dR = gamma * I\n",
    "    \n",
    "    # Return derivative\n",
    "    return [dS, dI, dR]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8bb2164",
   "metadata": {},
   "source": [
    "Define a cost function for quantifying the sum of squared errors - the difference between the data and our model:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6470fd5a",
   "metadata": {},
   "source": [
    "## Question: What is R0 for your fitted model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "75c1e9fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sir_sse(parameters, dat):\n",
    "    # Calculate model output\n",
    "    output = odeint(sir_model, initial_state_values, times, args=(parameters,))\n",
    "    output = pd.DataFrame(output, columns=['S', 'I', 'R'])\n",
    "    output['time'] = output.index + 10 # plus nine, because we are going to start on week 10\n",
    "    \n",
    "    # Merge with data\n",
    "    output = pd.merge(output, dat, on='time')\n",
    "    \n",
    "    # Calculate sum of squared errors (SSE)\n",
    "    output['sse'] = (output['I'] - output['Deaths'])**2\n",
    "    \n",
    "    sse = output['sse'].sum()\n",
    "    \n",
    "    return sse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6b155d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = {'beta': 1.0, 'gamma': 0.5}\n",
    "\n",
    "# Define times - start at week 10\n",
    "times = np.arange(10, 61, 1)\n",
    "\n",
    "# Initial state values\n",
    "initial_state_values = [54999, 1, 0]\n",
    "\n",
    "# Perform the minimization of Sum of Squared Errors\n",
    "fitval = minimize(sir_sse, x0=list(parameters.values()), args=(plague_dat,), method='BFGS')\n",
    "\n",
    "# Parameter values from the optimization:\n",
    "fitval_params = fitval.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "35cbb07b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.67546942, 2.20012301])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitval_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9e70505f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the model with the optimized parameters\n",
    "sir_fit = odeint(sir_model, initial_state_values, times, args=(fitval_params,))\n",
    "sir_fit = pd.DataFrame(sir_fit, columns=['S', 'I', 'R'])\n",
    "sir_fit['time'] = times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32103f46",
   "metadata": {},
   "source": [
    "Create a plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f5f0fbb4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvAAAAIoCAYAAAAPy984AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABtcklEQVR4nO3deXxcdb3/8denyXRNV7qwdNrSFgpYCEtJAaXUgizX6wp4uXhF/CkEtUEFf1FxAe5VkYiyVK9WrwhXyqLgBj9AEEhLWdKCNMhONxqgULo3XdJp8v39cWYySybJTDIz58zM+/l4zGP4nvM9Zz7TL8l85pvP+R5zziEiIiIiIsVhgN8BiIiIiIhI5pTAi4iIiIgUESXwIiIiIiJFRAm8iIiIiEgRUQIvIiIiIlJElMCLiIiIiBSRSr8DKCZjx451U6ZM8eW1d+7cybBhw3x5bUlPYxJMGpfg0ZgEk8YleDQmweTXuDz77LMbnXPj0u1TAp+FKVOm8Mwzz/jy2o2NjcydO9eX15b0NCbBpHEJHo1JMGlcgkdjEkx+jYuZvdHdPpXQiIiIiIgUESXwIiIiIiJFRAm8iIiIiEgRUQIvIiIiIlJElMCLiIiIiBQRJfAiIiIiIkVECbyIiIiISBFRAi8iIiIiUkR0IycRERGRPmhra2Pz5s3s2LGD9vb2fp9v5MiRvPzyyzmITHIpF+NSUVHB8OHDGTNmDIMGDep3TErgRURERLLU1tbGunXrGD16NFOmTCEUCmFm/Trnjh07GD58eI4ilFzp77g454hEImzfvp1169YxadKkfifxKqERERERydLmzZsZPXo0Y8eOZeDAgf1O3qV0mRkDBw5k7NixjB49ms2bN/f7nErgRURERLK0Y8cORowY4XcYUmRGjBjBjh07+n0eJfAiIiIiWWpvbycUCvkdhhSZUCiUk+sllMCLiIiI9IHKZiRbufp/Rgm8iIiIiEgRUQIvIiIiIlJElMCLiIiIiBQRJfAiIpKVlpYW6urqqKmpoa6ujpaWFr9DEhGfmJmuBfCBbuQkIiIZa2lpobq6mtbWViKRCCtWrGDRokU0NzcTDof9Dk9EpCxoBl5ERDLW0NDQmbwDRCIRWltbaWho8DkyEZHyoQReREQy1tTU1Jm8x0QiEZYtW+ZTRCIi5UcJvIiIZGz27Nldbl4TCoWoqanxKSIRkfKjGngRkVK3eTMsWgSRCITD3mPSJJgwASoqsjpVfX09ixYt6iyjCYVCVFVVUV9fn6fgRYpMPy7oHJ7DMDLmnB+vKv2kBF5EpJTt2wcf+AC8/HLXfZWVcNBB8YQ+ltyfeCIce2za04XDYZqbm2loaGDZsmXU1NRQX1+vC1hFRApICbyISClbujR98g5ecv/GG94j1Q9+AFdckfawcDjMggULchikiIhkQzXwIiKl7LHH+nbcD34Ae/fmNhYREckJJfAiIqWssTG5fdxxcPTRMGZMz8ft2gXLl+crKpHS5VyfHzu2b+/X8X16SFFSAi8iUqp274ann07edvfd8NxzsGkTtLbCK6/AQw/Bb34DJ5yQ3Hfx4sLFKiIiGVMNvIhIqXr66eQymClTvEfMsGEwY4b3AGhvT074Gxu7rYMXERH/aAZeRKRUpZbPfPCDPfefOze5/cQTqoMXEQkgJfAiIqUqNYFPTdBTTZ8OBx4Yb+/aBc88k+uoRESkn5TAi4iUonT176ec0vMxZl2TfNXBi0gPnHM4XQxbcErgRURKUWr9+8EHw+TJvR+XmuSnzuKLiIjvlMCLiJSi1PXfeyuf6a7f0qUQieQiIhERyREl8CIipSjb+veYQw6BAw6It1UHLyISOErgRURKza5d0NSUvK23+veYdHXwKqMREQkUJfAiIqWmr/XvManJvi5kFREJFCXwIiKlpq/lM931Vx28iEigKIEXESk12d7AKdWhh8L++8fbO3fCs8/2OywREckNJfAiIqWkP/XvMaqDFxEJNCXwIiKl5Kmnkuvfp06FSZOyP4/WgxcRCSwl8CIipaS/9e/dHffEE6qDFxEJCCXwIiKlJFcJ/IwZMGFCvN3aCv/4R1+jEhGRHFICLyJSKvpZ/97S0kJdXR01NTXUXXopu2pqkjuojEZEJBAq/Q5ARERy5Kmnkstcsqh/b2lpobq6mtbWViKRCCtWrGBgKMRPEjs1NsI3vpHLiEVEpA80Ay8iUir6UT7T0NDQmbwDRCIRHk68GBa89eD37etfjCIiAbF27VrMjAsvvNDvULKmBF5EpFT0I4FvamrqTN5j/rlvH5srE/5Qqzp4EUlhZkmPQYMGMW7cOI499li+8IUv8MADD9De3u5rfHP7ei1QgKmERkSkFKSrf8/iQ2v27NmsWLEiKYkPhUKsnTyZMStXxjs2NkJqbbyIlL0rr7wSgPb2drZu3cqLL77I7373O37zm98wa9YsFi1axKGHHupzlKVDCbyISCl48snk+vdp0yAczvjw+vp6Fi1a1FlGEwqFqKqq4uALL4TvfCfesbER6utzFraIlIarrrqqy7Z3332Xuro6/vCHP3DaaafxzDPPMH78+MIHV4JUQiMiUgr6uXxkOBymubmZ2tpaampqqK2tpbm5mdGf+ERyx8cfVx28iGRkwoQJ3HnnncydO5eWlhZ++MMfdumzefNmvvWtb3H44YczZMgQRo4cyamnnspDDz3Upe+2bdv48Y9/zLx585g4cSIDBw5k3LhxfPSjH+Xpp59O6nvLLbdgZgAsXrw4qcwn3ZeNtWvXct555zF27FgGDx7MrFmzuO+++7r027t3LzfddBPHHnsso0ePZujQoUyZMoWPfexj/P3vf+/jv1T2NAMvIlIKcrD+ezgcZsGCBckbnYNx4+C997x2ays89xwcf3yfwhSR8jJgwAC+853v0NjYyB133MH111/fmVi/8cYbzJ07l7Vr13LyySdz5plnsnPnTu677z7OPPNMFi5cyEUXXdR5rpdffplvf/vbzJkzhw9/+MOMHj2adevW8de//pUHHniAe++9lzPPPBOAo48+miuvvJKrr76ayZMnJ12omloT/8Ybb1BTU8PUqVP5zGc+w+bNm7nrrrs6k/JZs2Z19r3wwgu54447mDlzJhdccAFDhgzh7bffZunSpTz44IOcdtpp+fvHTOSc0yPDx3HHHef88thjj/n22pKexiSYynJcWludC4Wc89Jt77FuXe7Of+65yeduaMjq8LIckyKgcemfl156Kefn3L59e5dt69atc/Pnz3fHH3+8mz9/vluXy5/tHACcl052b8+ePa6ystIBbvXq1Z3bTznlFGdm7o477kjqv2XLFlddXe0GDx7s3nnnnc7tW7dude+9916X87e0tLgDDjjAHXbYYWnjO+WUU9LGtWbNms74r7rqqqR9Dz74oAPcWWed1TkuW7dudWbmjjvuOLdv374u59u4cWP3/wgJMv1/B3jGdZOTqoRGRKTYpa7/nmX9e69SbwalGzqJFETs/gwLFy5k+fLlLFy4kOrqalpaWvwOLSuDBg1iv/32A+C96F/zmpubWbx4MWeffTbnnXdeUv9Ro0Zx9dVXs2fPHu65557O7SNHjmTs2LFdzj9x4kTOOeccXnnlFdatW5d1fJMnT+Y7idf6AGeccQaTJk1i2bJlndvMDOccgwYNYsCAril07D0WgkpoRESKXQ7KZ3qUer5YHXylPkJE8ind/RlaW1tpaGjoWu4WcN6EMp3lM0899RTg1bWnq0mPJfovv/xy0vYnnniCG2+8kaeeeooNGzawN+V+FW+99RaTMryBXczRRx9NRUVFl+3hcLgzToARI0bwkY98hHvvvZejjz6as88+m5NPPpnZs2czdOjQrF6zv/TbV0Sk2OU7gT/iCBg7FjZu9No7dqgOXqQA0t2fIRKJJM0KF4M9e/awefNmAMaNGwfApk2bAHj44Yd5+OGHuz22tbW187//9Kc/cc455zB48GA+9KEPMW3aNIYNG8aAAQNobGxk8eLFtLW1ZR3fqFGj0m6vrKyko6Mjadtdd93Ftddey+233965dObgwYM555xzuO6665gwYULWr98XSuBFRIrZzp2Q+mGe6wTezCujSfhTNosXK4EXybPu7s9QU2T3Yli6dCn79u1jwoQJTJkyBfDKYQBuvPFGLr300ozO893vfpeBAwfyzDPPcPjhhyftq62tZfHixTmNO50hQ4Zw1VVXcdVVV9HS0sKSJUu45ZZbuO2221i7di2PP/543mMALSMpIlLcUtd/nz4dJk7M/eukfilQHbxI3tXX11NVVUUoFALovD9DfRHdi6Gjo4Mf/OAHAJx//vmd20844QSArBLelStXcsQRR3RJ3js6Oli6dGnaYwYMGJC3O8GGw2E+/elP87e//Y1DDjmEpUuXdv5lId+UwIuIFLN8l890d16tBy+Sd93dnyGcy4vU82jDhg2cd955NDY2MmnSJK644orOfbNmzeLkk0/mj3/8IzfffHPa4//5z3+yYcOGzvaUKVN4/fXXefvttzu3Oee4+uqreemll9KeY7/99svZRb/vvfceTal3vAZ27tzJjh07qKysZODAgTl5rd6ohEZEpJgVKoFPrYPfvh1WrICE9ZFFJPfS3p8hgGIXonZ0dLB161ZefPFFli5dyt69e6mpqWHRokVdVpC5/fbbmTdvHp///Oe56aabmD17NqNGjeLNN9/k+eef54UXXuCpp57qvHvr1772NS655BKOOeYYzj77bEKhEE888QQvvfRS58WlqU499VTuvPNOPvKRj3DcccdRWVnJnDlzmDNnTtbv8a233uKEE07g8MMP59hjjyUcDrN9+3buu+8+3nnnHS699FKGDx+e/T9eHyiBFxEpVunq31OXfMyVAQNgzhz44x/j2xYvVgIvIgBcffXVAAwcOJDhw4czefJkLrjgAs4++2xOP/30tMsuTpw4kWeffZYFCxZwzz33sGjRItrb29l///054ogjqKur48gjj+zsX1tby6BBg7jhhhu49dZbGTJkCCeffDK//e1vueeee9Im8DfeeCNmxiOPPML9999PR0cHV155ZZ8S+ClTpnD11VfT2NjIY489xsaNGxkzZgwzZszgRz/6UZflMPPJYsv6SO9mzZrlnnnmGV9eu7Gxscudw8RfGpNgKqtxefhhOP30eHv6dHj99fy93oIFkHix2b/+K6T5wExVVmNSRDQu/fPyyy93qcXurx07dhRsBlcyl+txyfT/HTN71jmXdpZENfAiIsWqUOUz3Z1/yRLI08VhIiLSPSXwIiLF6rHHktv5TuDf9z5IvNNgrA5eREQKSgm8iEgxam2F5cuTt+U7gY/VwSfScpIiIgWnBF5EpBg9/XTyMo6HHAIHHZT/1039klCAG6eIiEgyJfAiIsUodc3jk04qzOuqDl5ExHdK4EVEitFrryW3DzusMK87cyaMGRNvb9sGq1YV5rVFRARQAi8iUpxSl4s85JDCvO6AAZCwLjMAr75amNcWERFACbyISHFKnYE/9NDCvfaMGcntV14p3GuLBIjupSPZytX/M0rgRUSKTVsbrFuXvG3atMK9fmq5jmbgpQxVVFQQiUT8DkOKTCQSoaKiot/nUQIvIlJsVq+Gjo54OxyGoUML9/qpCbxm4KUMDR8+nO3bt/sdhhSZ7du35+SurkrgRUSKTWr5TKHq32NUQiPCmDFj2LJlCxs3bmTv3r0qp5FuOefYu3cvGzduZMuWLYxJXAigjypzEJeIiBRS6gWshax/B5g8GQYN8kp5ADZtgo0bYezYwsYh4qNBgwYxadIkNm/ezNq1a2nPwXKqe/bsYfDgwTmITnIpF+NSUVHB8OHDmTRpEoMGDep3TErgRUSKjd8z8BUV3mu+8EJ826uvKoGXsjNo0CAOOOAADjjggJycr7GxkWOOOSYn55LcCeK4BKaExsw+bGYPmdmbZrbbzFab2R/M7MRu+p9kZveb2WYz22Vmz5vZV82s2ysDzOyzZrbMzFrNbJuZNZrZv+bvXYmI5IFfS0gm0oWsIiK+CUQCb2bXAvcBxwIPAjcC/wA+BjxhZv+R0v9jwBJgDvAn4OfAQOB64M5uXuM64BbgAODXwG3AkcC9ZjY/529KRCRf/FxCMkYXsoqI+Mb3BN7M9ge+DrwLHOGc+4Jz7pvOuXOAMwAD/jOh/wi8BLwdmOuc+7xz7v8CRwNPAeeY2Xkpr3EScDmwCjjKOfc159yXgeOAzcB1ZjYlv+9URCQHWlvh7bfj7QED4OCDCx9H6oWsaWbgW1paqKur45VXXqGuro6WlpYCBSciUtp8T+CByXhxNDnnNiTucM49BuwAxiVsPifavtM590xC3z3Ad6LNL6a8xiXR5x8457YkHLMWb/Z+EPC5fr8TEZF8W7kyuX3wwTBwYOHj6GUGvqWlherqahYuXMjOnTtZuHAh1dXVSuJFRHIgCAn868BeoMbMkq6AMrM5wHDg7wmb50WfH0xzriXALuAkM0u8xLenYx5I6SMiElxBqH+HrmU7q1bB3r2dzYaGBlpbWztvdBOJRGhtbaWhoaGQUYqIlCTfE3jn3GbgG8AE4CUz+5WZXWNmvwceAh4GahMOif3dNqUIFJxz+4A1eKvrTAUws2HAQUCrc259mhBin4Y+FJGKiGTJ7yUkY0aMgAMPjLfb270bTEU1NTV1uUtlJBJh2bJlhYpQRKRkBWIZSefcDWa2FrgZuChh10rglpTSmpHR523dnC62fVQf+ycxs4uBiwEmTJhAY2NjN6fJr9bWVt9eW9LTmARTqY/LYUuWsH9C+3XgLZ/eb/WECYxOqMd/4e672fiBDwBQV1fHe++9h3OOiRMnct1112FmjBs3rqTHp5iU+s9KMdKYBFMQxyUQCbyZ1QM/BG4Cfga8AxwGXAMsMrOjnXP1mZ4u+pztLdHS9nfO/Qr4FcCsWbPc3LlzszxtbjQ2NuLXa0t6GpNgKvlx+fa3k5qHfPjDHOLX+z3hBHjuuc7mzFAIorHEauBbW1u55ppr+Na3vkVVVRXNzc2Ew2F/4pUkJf+zUoQ0JsEUxHHxvYTGzOYC1wJ/dc5d5pxb7Zzb5Zz7B/AJ4C3gcjObGj0kNmM+ssvJPCNS+vXWv7cZehGR4PD7Jk6JeriQNRwO09zcTG1tLcOGDaO2tlbJu4hIjviewAOxGyk9lrrDObcLWIYXZ+wWWLG1yroUfppZJXAwsA9YHT3HTrwvAVVmlu5WabFPvy419SIigbJlC2zcGG8PHAiTJvkXT+pSkikr0YTDYRYsWMBhhx3GggULlLyLiORIEBL42Gox47rZH9seW97g0ejzmWn6zgGGAk8659oStvd0zFkpfUREgin1AtZp06Ci25tP51+6u7G6bKsXRUQkW0FI4B+PPl9sZgcl7jCzs4D3A3uAJ6Ob7wY2AueZ2ayEvoOB70ebv0h5jV9Gn79tZqMTjpkCfBloA37b73ciIpJPQVlCMiYchiFD4u0tW+C99/yLR0SkTAThIta78dZ5Pw142cz+hHcR6+F45TUGfNM5twnAObfdzC6KHtdoZnfi3U31o3hLTN4N3JX4As65J83sp8BlwPNmdjcwEPg3YAxQF72pk4hIcAVlCcmYAQO8GJqb49tefRXGj/cvJhGRMuD7DLxzrgP4F+BrwEt4F65eDpwA3A+c4Zy7MeWYPwOn4N246WygDojgJejnOdf1b7jOucuBC/G+HFwMXAC8CHzEOfezPLw1EZE+aWlpoa6ujpqaGurq6uJ3Lw3SBawxvdyRVUREci8IM/A45yLADdFHpsc8gZf4Z/M6twK3ZnOMiEghJS6/GIlEWLFiBYsWLfJWcAnaDDz0eiGriIjknu8z8CIiEtfQ0NCZvIN399LW1lYarr22OGbgX301fT8REckZJfAiIgHS1NTUmbzHRCIRXn/ySdi+Pb5x6FA48MACR5eGZuBFRApOCbyISIDMnj2bUCiUtC0UCnHGwQcndzzkEDDDd6llPGvWQFtb+r4iIpITSuBFRAKkvr6eqqqqziQ+FApRVVXFZ086KbljEMpnAKqqvOUkYzo6YOVK/+IRESkDSuBFRAIkHA7T3NxMbW0tNTU11NbW0tzczJhNm5I7BuEC1hiV0YiIFFQgVqEREZG4cDjMggULkjcG8QLWmMMOg7//Pd7WhawiInmlGXgRkWIQxCUkYzQDLyJSUErgRUSCrqOjawIftBn4RJqBFxHJKyXwIiJB9/bbsHt3vD1qFIwd61s4XaS7G2vXG2KLiEiOKIEXEQm6dPXvQVhCMuagg2DYsHh7+3Z45x3/4hERKXFK4EVEgi7I9e/gfZlIrYNXGY2ISN4ogRcRCbog17/H6EJWEZGCUQIvIhJ0QV5CMkYXsoqIFIwSeBGRoAt6CQ2kv5BVRETyQgm8iEiQ7dsHq1YlbwviDLxKaERECkYJvIhIkK1bB5FIvD1+PIwc6V883UldGeeNN5KXvhQRkZxRAi8iEmTFUP8OMHQoTJoUbzvXtfRHRERyQgm8iEiQFUP9e4wuZBURKQgl8CIiQVYMS0jG6EJWEZGCUAIvIhJkqSU0QZ6B182cREQKQgm8iEiQaQZeRERSKIEXEQmqvXth7drkbdOn+xJKRtLNwDvnTywiIiVMCbyISFCtXg0dHfH2xIneai9BdcABMHx4vN3aCm+/7V88IiIlSgm8iEhQFcsSkjFmuqGTiEgBKIEXEQmqYlpCMkZLSYqI5J0SeBGRoCqmC1hjdCGriEjeKYEXEQmqYlpCMkYlNCIieacEXkQkqEphBl4lNCIiOacEXkQkiHbtgjffjLcHDICpU/2LJ1PTp3sXs8asW8eA3bv9i0dEpAQpgRcRCaKVK5PbU6bAwIG+hJKVwYPh4IOTNg196y2fghERKU1K4EVEgqjYlpBMlFJGM3TdOp8CEREpTUrgRUSCqBiXkIxJuZBVCbyISG4pgRcRCaISmoEf0tLiUyAiIqVJCbyISBBpBl5ERLqhBF5EJIiKcQnJmNQa+DffhI4On4IRESk9SuBFRIJm2zbYsCHeDoVg8uQu3VpaWqirq6Ompoa6ujpaglKqMn48jBrV2azYsyd5SUwREemXSr8DEBGRFKmz79OmQUVF0qaWlhaqq6tpbW0lEomwYsUKFi1aRHNzM+FwuIDBpmHmldE0NcW3vfIKTJrkX0wiIiVEM/AiIkGTegFrmvr3hoaGzuQdIBKJ0NraSkNDQyEi7J3uyCoikjdK4EVEgiaD+vempqbO5D0mEomwbNmyLn19KbVJuZCVV17J/2uKiJQJldCIiARNBktIzp49mxUrViQl8aFQiJqamqR+vpXaaAZeRCRvNAMvIhI0GSwhWV9fT1VVFaFQCPCS96qqKurr65P6+VZqoxl4EZG8UQIvIhIkzmVUQhMOh2lubqa2tpaamhpqa2vTzqpnU2qTU6kX3r71FuzYkd/XFBEpEyqhEREJko0bYevWeHvoUDjwwLRdw+EwCxYs6PF0mZba5NygQTB1avKXkVdfhVmz8vu6IiJlQDPwIiJBkjr7Pn06DOj7r+pMS23yIrX0Z+XK/L+miEgZUAIvIhIkGSwhmY1MS23yYvr05PaqVfl/TRGRMqASGhGRIElNctPUv2crk1KbvEhN4DUDLyKSE5qBFxEJktWrk9tTp/oTRy4ogRcRyQsl8CIiQbJmTXL74IP9iSMXpk1LbquERkQkJ5TAi4gESSkl8JMnJy8luX497NzpXzwiIiVCCbyISFDs3g3vvBNvDxgAhbjYNF8GDvSS+ESahRcR6Tcl8CIiQbF2bXI7HIbo8o9FSyvRiIjknBJ4EZGgKKXymZjUOnhdyCoi0m9K4EVEgqIUE3itRCMiknNK4EVEgkIJvIiIZEAJvIhIUJRiAq+lJEVEck4JvIhIUJRiAp96I6p166CtzZ9YRERKhBJ4EZGgKMUEfsgQ9owbF2871/V9iohIVpTAi4gEwdat3iNm0CDYf3+/osmp3QcdlLxBZTQiIv2iBF5EJAhSZ6WnTPFu5FQC9hx4YPIGXcgqItIvpfHpICJS7FavTm6XQvlMVJcZeCXwIiL9ogReRCQISrH+PUoJvIhIbimBFxEJgtQEPnX1liK2O7WERjXwIiL9ogReRCQISnkGPjWBX7MG9u3zJxgRkRKgBF5EpABaWlqoq6ujpqaGuro6WlpakjuUcALfPmwYjB8f37Bvn7cevIiI9Eml3wGIiJS6lpYWqquraW1tJRKJsGLFChYtWkRzczPhcNhbG33t2uSDSiiBB2D6dNiwId5etaqkyoRERApJM/AiInnW0NDQmbwDRCIRWltbaWho8Dq88w7s2RM/YORIGD3ah0jzaNq05LYuZBUR6TMl8CIiedbU1NSZvMdEIhGWLVvmNUq4fKbT9OnJbSXwIiJ9pgReRCTPZs+eTSgUStoWCoWoqanxGuWYwGslGhGRPlMCLyKSZ/X19VRVVXUm8aFQiKqqKurr670O5ZDAq4RGRCRnlMCLiORZOBymubmZ2tpaampqqK2tjV/ACuWRwKebge/o8CcWEZEip1VoREQKIBwOs2DBgvQ7yyGBHzMGRo2CrVu99p498PbbMHGin1GJiBQlzcCLiPitHBJ4M9XBi4jkiBJ4ERE/RSJdb2o0ZYovoeSd6uBFRHJCCbyIiJ9aWpJrwSdMgKFD/Ysnn7SUpIhITiiBFxHxUzmUz8SohEZEJCeUwIuI+KmcE3jNwIuI9IkSeBERP6Um8FOn+hNHIaSrgXfOn1hERIqYEngRET+V0wz8/vsn1/fv2AEbN/oXj4hIkVICLyLip3JK4NMtJakyGhGRrCmBFxHxUzkl8KClJEVEckAJvIiIX3btgnffjbcrKiAc9i+eQtAMvIhIvwUqgTezk83sHjNbb2Zt0eeHzOxf0vQ9yczuN7PNZrbLzJ43s6+aWUUP5/+smS0zs1Yz22ZmjWb2r/l9VyIi3Vi7NrkdDkNlpS+hFIyWkhQR6bfAJPBm9h1gCTAHeBD4CXAvMBqYm9L3Ywl9/wT8HBgIXA/c2c35rwNuAQ4Afg3cBhwJ3Gtm83P9fkREelVu5TOgGXgRkRwIxFSPmZ0L/Bfwd+CTzrkdKftDCf89Ai8BbwfmOueeiW7/LvAocI6ZneecuzPhmJOAy4FVwPHOuS3R7T8GngWuM7P7nHNr8/cuRURSlGMCrxp4EZF+830G3swGANcCu4DzU5N3AOdcJKF5DjAOuDOWvEf77AG+E21+MeUUl0SffxBL3qPHrMWbvR8EfK5/70REJEvlmMBPnAgDB8bbmzbB1q2+hSMiUox8T+CBk4CDgfuBLWb2YTP7hpl9xcxOTNN/XvT5wTT7luB9ETjJzAZleMwDKX1ERAqjHBP4ioquN6tSHbyISFaCkMAfH31+F/gHcB/wI+AG4EkzW2xm4xL6z4g+v5Z6IufcPmANXmnQVAAzGwYcBLQ659anef3Xo8+H9u9tiIhkafXq5HY5JPCgMhoRkX4KQg38+OjzJXjJ92lAEzAZ70LWM4A/EL+QdWT0eVs354ttH9XH/knM7GLgYoAJEybQ2NjYzWnyq7W11bfXlvQ0JsFUNOPiHB9YuTLpl/CT69eztxhiz1LqmEwfPJiJCftXP/ww6yZMKHhc5a5oflbKiMYkmII4LkFI4GPLPhpwjnOuOdp+0cw+gTfTfoqZneiceyqD81n02WUZR9r+zrlfAb8CmDVrlps7d26Wp82NxsZG/HptSU9jEkxFMy6bN8POnfH24MGc9MlPencrLTFdxuSFF+CeezqbUzs6mFoMY1ZiiuZnpYxoTIIpiOMShBKa2EWlqxOSdwCcc7uBv0WbNdHn2Iz5SNIbkdKvt/69zdCLiOReav37lCklmbynpaUkRUT6JQgJ/KvR563d7I8l+ENS+nepWTezSrwLYvcBqwGcczuBt4AqMzsgzfkPiT53qakXEcmbcryANUY18CIi/RKEBH4JXsJ9iJkNTLN/ZvR5bfT50ejzmWn6zgGGAk8659oStvd0zFkpfURE8q+cE/jJk73VaGLWr08uJxIRkR75nsA75zYCd+GVsnwvcZ+ZfQjvItZtxJeAvBvYCJxnZrMS+g4Gvh9t/iLlZX4Zff62mY1OOGYK8GWgDfhtDt6OiEhmUhP41KUVS9nAgV4Snyh1RR4REemW7wl81GXASrwEe4mZXWdmf8Bbo70duMg5txXAObcduAjv4tdGM/sfM2sAVgAn4iX4dyWe3Dn3JPBTYBrwvJldb2Y/B54BxgBf111YRaSgynkGHlQHLyLSD0FYhQbn3AYzm413J9VPACcAO4D/B1zjnHs6pf+fzewU4NvA2cBgvC8AlwE3Oee6rCjjnLvczJ4H5uMtC9mBt+78j51z9+XtzYmIpFPuCXxqHbxu5iQikrFAJPAAzrnNeAn4ZRn2fwL4lyxf41bg1uyjExHJoY4OWLs2eVu5JfApM/BLb7mFy37/e2bPnk19fT3hcNinwEREgi8oJTQiIuXjnXegLeE6+1GjvEc5SUng2156ieXLl7Nw4UKqq6tpaWnxKTARkeBTAi8iUmjlXj4DXUpopkUrHyORCK2trTQ0NPgRlYhIUVACLyJSaErgu6y6EwZi6whHIhGWLVtW8JBERIqFEngRkUJLXTKxHBP4IUNg4sTOZgUwJfrfoVCImpqadEeJiAhK4EVECk8z8J6UOvjpeMl7VVUV9fX1/sQkIlIElMCLiBSaEnhPSh38vEmTqK2tpbm5WavQiIj0IDDLSIqIlA0l8J6UGfjLP/YxuOkmn4IRESkemoEXESmkSATefDN525QpvoTiO92NVUSkT5TAi4gU0rp13o2cYvbf37ugsxyl3o1VCbyISEaUwIuIFJLKZ+JSE/i1a2HfPl9CEREpJkrgRUQKSQl83IgRMH58vB2JgO7AKiLSKyXwIiKFlJrAp9zQqOyoDl5EJGtK4EVECkkz8MlSy2hWrfInDhGRIqIEXkSkkJTAJ9MMvIhI1pTAi4gUkhL4ZErgRUSypgReRKRQdu6EDRvi7YoKmDjRv3iCQCU0IiJZUwIvIlIoa9cmtydNgsoyvyF26gz8qlXJ6+SLiEgXSuBFRApl9erkdrmXzwCMGQOjRsXbu3fD+vW+hSMiUgyUwIuIFIrq37syUx28iEiWlMCLiBSKEvj0UuvglcCLiPRICbyISKEogU/vkEOS26++6k8cIiJFQgm8iEihKIFPb8aM5PZrr/kTh4hIkVACLyJSCM4pge9OagKvGXgRkR4pgRcRKYTNm2HHjnh7yBCYMMG/eIIkNYFftQr27fMnFhGRIqAEXkSkEFJn36dM8VZgERgxAvbfP96ORLr+e4mISKecJvBmNtrMhuXynCIiJUHlMz1TGY2ISMayTuDN7FQzazCz0QnbxpvZYmAjsNnMfprLIEVEip4S+J4demhyWwm8iEi3+jIDXwd80jm3JWHbdcDJwEpgE/AVM/tUDuITESkNqQn81Kn+xBFUWolGRCRjfUngq4GlsYaZDQHOAR52zs0AZgAtwCU5iVBEpBSk3pxIM/DJVEIjIpKxviTw44G3E9qzgcHALQDOuR3AfXiJvIiIQNcEPvXmReVOCbyISMb6ksC3AUMS2icDDliSsG07MKYfcYmIlI62Nli3LnmbSmiSHXwwVFbG2++8A9u3+xePiEiA9SWBXwPMS2ifDbzunHsrYVsY74JWERFZuxY6OuLtiRNh6FDfwgmkykqYNi15m2bhRUTS6ksCfytwpJk1mdnjwJHA7Sl9jgX0m1dEBOD115Pb06f7E0fQ6UJWEZGMVPbepYtfACcA/wYYcC9wbWynmdUAhwN35CJAEZGil1r/rgQ+PdXBi4hkJOsE3jkXAc43s0u8ptuR0mU1cAywtv/hiYiUACXwmVECLyKSkb7MwAPgnEt7dZFzbiOqfxcRiVMCnxkl8CIiGelLDbyIiGRDS0hmJl0NfOLFvyIiAvQxgTezU8zsPjPbYGYRM2tP89iX62BFRIpOJOKtQpModbUV8YwdC6NGxdu7d8Obb/oWjohIUGVdQmNmHwb+DFQA6/BWm1GyLiKSztq10N4ebx9wAAwb5ls4gWbmzcI3NcW3vfYaTJrkX0wiIgHUlxr4q4AI8GHn3EO5DUdEpMSo/j07qQn8q6/Caaf5F4+ISAD1pYRmJnCXkncRkQyo/j07upBVRKRXfUngW4HNuQ5ERKQkaQY+O0rgRUR61ZcE/hHgxFwHIiJSknQX1uwogRcR6VVfEvhvANPM7DtmZrkOSESkpGgGPjvTpnkXs8asW+etRiMiIp16vYjVzG5Os/lF4Grg/5jZCmBrmj7OOff5fkUnIlLM9u2DNWuStymB79mQITB5cnzpTee8L0FHHulrWCIiQZLJKjQX9rBvSvSRjgOUwItI+Vq3zkviYyZMgOHD/YunWMyYkbx2/quvKoEXEUmQSQJ/cN6jEBEpRSqf6ZsZM+Bvf4u3VQcvIpKk1wTeOfdGIQIRESk5uoC1b3Qhq4hIj7K+iNXMvmdmc3rpc7KZfa/vYYmIlADNwPeNEngRkR71ZRWaq4C5vfSZA1zZh3OLiJQOJfB9c+ihye3XXvMuZhUREaBvCXwmKoGOPJ1bRKQ46C6sfXPQQTB0aLy9dSu8955v4YiIBE2+EvjjgI15OreISPC1t8Pq1cnbpk3zJ5ZiM2BA11l4ldGIiHTKZBUazOzRlE0XmtncNF0rgDAwGbijX5GJiBSzlhbYuzfeHjsWRo3yLZyiM2MGrFgRb7/6Kpx8sm/hiIgESUYJPMk1747u13/vADYBdwFf60dcIiLFTfXv/aMLWUVEupVRAu+c6yy1MbMO4Crn3H/mLSoRkWKnBL5/lMCLiHQr0xn4RJ8Dnst1ICIiJUUXsPZPupVoREQE6EMC75y7NR+BiIiUFM3A909qAr9qFUQiEAr5E4+ISID0ZQa+k5lNBA4CBqXb75xb0p/zi4gULd2FtX9GjIADDoD16732vn2wZk3XxF5EpAz1KYE3s9OB64HDeula0Zfzi4gUtY4Ob8Y4kRL47M2YEU/gwauDVwIvIpL9OvBmNhu4DxgF/AwwYAnwa+CVaPteQBe5ikh5eustaGuLt8eM8R6SHV3IKiKSVl9u5HQFsAc43jn3lei2x5xzlwAzgf8CTgPuzk2IIiJFRvXvuZGawOtCVhERoG8J/InAX51zb6eex3muBF4Grs5BfCIixUcJfG7obqwiImn1JYEfCaxLaO8FhqX0eQKY09egRESKmi5gzQ2V0IiIpNWXBH4DMDqlPS2lTwgY0tegRESKmmbgc2PKlORlI999F7Zt8y0cEZGg6EsC/xrJCfvTwIfM7FAAM9sfOBt4Pc2xIiIlpaWlhbq6Ompqaqirq6OlpUU3ccqVysquX340Cy8i0qdlJB8Evm9mY5xzm4EbgU8Cz5nZS8AhwHCgPndhiogET0tLC9XV1bS2thKJRFixYgWLbruNjXv2JM+OaAa+72bMgJdfjrdffRVqavyLR0QkAPoyA78Qr749AuCcewI4F1iDtwrNeuCLzrn/zVWQIiJB1NDQ0Jm8A0QiEYa3tjJgz554p5EjYb/9fIqwBGglGhGRLrKegXfObQeaUrb9CfhTroISESkGTU1Nncl7zJR9+5I7TZ8OZgWMqsRoJRoRkS76MgMvIiLA7NmzCSVeZAnMqEi5AbXKZ/pHK9GIiHTR5wTezMaZ2SVmdqOZ/U/K9hoz0yo0IlLS6uvrqaqq6kziQ6EQh6ck9LqAtZ9SE/jXX4eODn9iEREJiD4l8Gb2eWAt8HOgDvhcwu4JwFPA+f0NTkQkyMLhMM3NzdTW1lJTU0NtbS218+Yld9IMfP+MHQtjxsTbu3dDS4t/8YiIBEDWNfBm9iHgV8DzwJXAGcAlsf3OuRfM7EXg48BvchOmiEgwhcNhFixYEN9wzDHJHZTA99+MGfDUU/H2q6/C5Mn+xSMi4rO+zMB/A2+lmVOcc3/Fu5FTqueBI/oTmIhI0XFOd2HNB61EIyKSpC8J/CzgvuhqNN15E9i/byGJiBSpd9+FnTvj7aoqGD/ev3hKhVaiERFJ0pcEfiCws5c+o4D2PpxbRKR4pbsDq5aQ7D+tRCMikqQvCfxa4Lhe+swG9BtWRMpLagKv8pncUAIvIpKkLwn8X4CTzezcdDvN7HPAUcA9/QlMRKToKIHPj+nTYUDCx9W6dbBrl3/xiIj4rC8JfAOwDrjDzO4CTgQws/nR9q+A14EF3Z9CRKQE6QLW/Bg0CKZMSd6W+mVJRKSMZJ3AO+e2AKcAS4FzgdMBA26Ktp8ETnXO9VYnLyJSWtLVwEtuqIxGRKRT1uvAAzjn1gFzzewovBn4/YBtwNPOuWdzGJ+ISHFwTiU0+XToofDAA/G2EngRKWN9SuBjnHPP4635LiJS3jZuhO0Jq+sOHQr7azXdnNEMvIhIp77UwANgZpPNbJaZHWdmk3IZlJl9xsxc9PGFbvqcZGb3m9lmM9tlZs+b2VfNrKKH837WzJaZWauZbTOzRjP711zGLiJlKt3su5aQzB0l8CIinbJK4M1srJn91MzWA6uBJmAZsMbM3jazH5vZmP4EZGZhvAtgW3vo8zFgCTAH+BPwc7z16a8H7uzmmOuAW4ADgF8DtwFHAvea2fz+xCwiogtY8ywlgW/9xz+omz+flpYWnwISEfFPxgm8mR0CPAN8BZiAd6OmDcB70f/eH7gMeMbMpvYlGDMz4LfAJuCX3fQZgZeAtwNznXOfd879X+Bo4CngHDM7L+WYk4DLgVXAUc65rznnvoy3nv1m4Dozm9KXmEVEAF3Amm8HHkjHsGGdzar2dv60cCHV1dVK4kWk7GSUwJvZAGARMAlYDJwGVDnnDnDO7Q8Mx1uNZgkwBW92uy8uBeYBn6P7u72eA4wD7nTOPRPb6JzbA3wn2vxiyjGXRJ9/EF1FJ3bMWrzZ+0HR1xQR6RtdwJpfZrw5dGjSpqn79tHa2kpDQ4NPQYmI+CPTGfjTgVnA7/GWiHzUObc3ttM51+ac+zte8n03MNvMPpRNIGZ2OPAj4Ebn3JIeus6LPj+YZt8SYBdwkpkNyvCYB1L6iIhkTwl83r3c3p7UngFEIhGWLVvmT0AiIj7JNIE/G2gD6pxzrrtO0X3zgQjeTHlGzKwS+B3eDaKu6KV7rBDytTSvvw9Yg7e6ztTouYcBBwGtzrn1ac4XK1w9NNN4RUS6UAKfd+3TpiW1DwdCoRA1NTX+BCQi4pNMl5E8FnjCOfdebx2dcxvMbGn0mEx9DzgG+IBzbncvfUdGn7d1sz+2fVQf+ycxs4uBiwEmTJhAY2NjL+HlR2trq2+vLelpTILJj3Gp3LaND2zprM6jfeBAHn/tNd0tNCpXY3LQGWfA8uWd7XOmT2fAl7/M+973Pv0s9oF+hwWPxiSYgjgumSbwYbw7r2bqReDfM+loZjV4s+4/cc49lcVrdHvK6HO3fynoRtr+zrlfAb8CmDVrlps7d27fI+uHxsZG/HptSU9jEky+jEtTU1Kz4tBDmTtPVXkxORuTgw6C73+/sznmzTc5+5OfJDwppysZlw39DgsejUkwBXFcMi2hGQFszeK8W/EubO1RQunMa8B3Mzx3bMZ8ZDf7R6T0661/bzP0IiI9U/lMYUybBokr0ezZQ3hAn29nIiJStDL9zTcQb9nGTHVEj+lNFV7t+eHAnoSbNzngymifX0e33RBtx+7e0aVmPfqF4GBgH9469TjndgJvAVVmdkCaGGJrvXWpqReR8tXS0kJdXR01NTXU1dX1vFShEvjCGDAAjjoqeVtzsz+xiIj4KNMSGsi+JCUTbcBvutl3LF5d/FK8pD1WXvMo8GngTOCOlGPmAEOBJc65toTtjwKfiR7z25RjzkroIyJCS0sL1dXVtLa2EolEWLFiBYsWLaK5uZlwONz1ACXwhVNdDU8lVFs2N8OHP+xfPCIiPsjmb49XmVl7Jg+8i1J75Zzb7Zz7QroH8Ndot1uj2+6Ktu8GNgLnmdms2LnMbDAQK478RcpLxW4K9W0zG51wzBTgy3hfJFITexEpUw0NDZ3JO3hLFfa43njqXVh1E6f8Ofro5PaKFX5EISLiq2xm4K33LknyMWOPc267mV2El8g3mtmdeHdT/SjeEpN3A3elHPOkmf0U706xz5vZ3XglPv8GjMFbHnNtPuIVkeLT1NTUmbzH9LjeuGbgC6e6OrmtEhoRKUMZzcA75wb04VGRr6Cdc38GTsG7cdPZQB3e2vOXAeelW6veOXc5cCHwDt6ykBfgrZbzEefcz/IVq4gUn9mzZxMKhZK2dbve+JYtsGlTvD1oEEycmOcIy9iRR4IlzCe9/jrs7O7G3SIipSmwl+87565yzplz7n+62f+Ec+5fnHOjnXNDnHNHOueud851e7Gtc+5W59zxzrlhzrnhzrlTnHP35e9diEgxqq+vp6qqqjOJD4VCVFVVUV9f37XzqlXJ7alTvYstJT+GDUsuUXIO/vlP/+IREfGBPmVERFKEw2Gam5upra2lpqaG2tpaXcAaJCqjEZEyl00NvIhI2QiHwyxYsKD3jrqAtfCqq+EPf4i3dSGriJQZzcCLiPSHZuALL3UlGs3Ai0iZUQIvItIfr6XcA27aNH/iKCepJTTPPw8dHf7EIiLiAyXwIiJ91dEBL7yQvO2II/yJpZwcdBCMGRNv79zZ9WJiEZESpgReRKSv1q6F1tZ4e9QoL7mU/DJTGY2IlDUl8CIifZW6fOFRRyWvUS75o5VoRKSMKYEXEemr559Pbh95pD9xlKPUBF4r0YhIGVECLyLSV+lm4KUwVEIjImVMCbyISF9pBt4/hx8O0TvlAtDSAps3+xePiEgBKYEXEemL3bu73sRp5kx/YilHAwd6SXwizcKLSJlQAi8i0hcvv5y89vjBB8Pw4f7FU45URiMiZUoJvIhIX6SWz6j+vfC0Eo2IlCkl8CIifZF6Aavq3wtPK9GISJlSAi8i0he6gNV/qQn8Sy9BJOJPLCIiBaQEXkSkL7SEpP/Gjk2+8+3evfDKK/7FIyJSIErgRUSytWEDvPtuvD1oEEyf7l885UxlNCJShpTAi4hkK3X2/X3vg8pKf2Ipd7qQVUTKkBJ4EZFs6QLW4NBSkiJShpTAi4hkS0tIBke6EhrnfAlFRKRQlMCLiGRLM/DBMX06DBkSb2/cCOvX+xePiEgBKIEXEclGezu88ELyNiXw/qmo6PoXEF3IKiIlTgm8iEg2Vq2CPXvi7XHjYMIE/+IRXcgqImVHCbyISDbS3cDJzJ9YxKMEXkTKjBJ4EZFs6AZOwZO6Eo1KaESkxCmBFxHJhi5gDZ7UMXj9ddi1y59YREQKQAm8iEg2tIRk8AwfDtOmxdsdHV0vNBYRKSFK4EVEMtXaCqtXx9tmcMQR/sUjcSqjEZEyogReRCRTL76YfJOg6dNh6FD/4pE4XcgqImVECbyISKZ0AWtwKYEXkTKiBF5EJFPplpCUYEgtoWlu9mrhRURKkBJ4EZFMaQY+uMJhGDUq3m5thTVrfAtHRCSflMCLiGTCOS0hGWRmKqMRkbKhBF5EJBPr18OmTfH20KEwdap/8UhXWolGRMqEEngRkUykzr7PnAkD9Cs0UDQDLyJlQp8+IiKZ0A2cgk8JvIiUCSXwIiKZUP178B1xBFRWxttvvAFbtvgXj4hIniiBFxHJhJaQDL7Bg+Gww5K3pY6biEgJUAIvItKbSARefjl5mxL4YEopo7nugguoq6ujpaXFp4BERHJPCbyISG9efx327o23DzgAxo71Lx7pXspKNGPWrWPhwoVUV1criReRkqEEXkSkN7qAtXikzMBXA5FIhNbWVhoaGvyJSUQkx5TAi4j0RhewFo+UBH4mUImXxC9btsyXkEREck0JvIhIbzQDXzzGj2fb0KGdzUHADCAUClFTU+NbWCIiuaQEXkSkN5qBLyoDUxL14yoqqKqqor6+3qeIRERySwm8iEhPtm3z1hOPqaiAww/3Lx7p1ZDZs5PanznqKJqbmwmHwz5FJCKSW5W9dxERKWMvvJDcnjEDBg3yJxbJTMpKNKcNHw5K3kWkhGgGXkTKSktLC3V1ddTU1GS2PrjKZ4rPCSckt5uaoK3Nn1hERPJAM/AiUjZaWlqorq6mtbWVSCTCihUrWLRoUc/lFbqAtfhMnuzNuMe+nLW1wfLl8IEP+BuXiEiOaAZeRMpGQ0NDZ/IOGa4Prhn44mMGc+Ykb1uyxJ9YRETyQAm8iJSNpqamzuQ9psf1wZ3rmsBrBr44nHJKcnvxYn/iEBHJAyXwIlI2Zs+eTSgUStrW4/rgLS3eKjQxI0bApEl5jFByJnUG/oknYN8+f2IREckxJfAiUjbq6+upqqrqTOJDoVDP64On1r/PnOmVZ0jwHXoojB8fb+/cCc895188IiI5pAReRMpGOBymubmZ2tpaampqqK2t7fkCVpXPFC/VwYtICdMqNCJSVsLhMAsWLMissy5gLW5z5sDdd8fbS5bA5Zf7F4+ISI5oBl5EpDtaQrK4pc7AP/44dHT4E4uISA4pgRcRSaetDV59NXnbzJn+xCJ9M3MmjBoVb2/ZAi++6Fs4IiK5ogReRCSdV15JXrVk0qTkZFCCr6Ki682bVAcvIiVACbyISDqqfy8NupBVREqQEngRkXRU/14a0iXwzvkTi4hIjiiBFxFJRzPwpeHYY2HYsHj7nXdg5Ur/4hERyQEl8CIiqZyDFSuStymBL06hEJx0UvK2xYv9iUVEJEeUwIuIpFq1ypupjRk6FGbM8C8e6R/VwYtIiVECLyKSKnWG9qSTvJlcKU5K4EWkxCiBFxFJlZrgpSaAUlxqamDgwHj7jTe8h4hIkVICLyKSKjWBP+UUf+KQ3Bg8GGbPTt72+OP+xCIikgNK4EVEEq1bB2vXxtsDB3ozuFLcVEYjIiVECbyISKLUxG72bG8GV4qbEngRKSFK4EVEEqn+vTSdeCJUVMTbr74K777rXzwiIv2gBF5EJJHq30vT8OHeTZ0SqQ5eRIqUEngRkZh33vFmZmMqKryZWykNqV/GVEYjIkVKCbyISEzqjOxxx0FVlT+xSO6llkPpjqwiUqSUwIuIxKh8prR94ANgFm//85+webN/8YiI9JESeBGRGF3AWtpGj4Yjj4y3nYMnnvAvHhGRPlICLyIC3kzsP/8Zb5vB+9/vXzySH1pOUkRKgBJ4ERGApUu9GdmYo47yZmyltCiBF5ESoAReRARU/14uTj45uf3ss9Da6k8sIiJ9pAReRAS6rkii+vfStP/+cOih8XZ7Ozz1lH/xiIj0gRJ4EZEdO+Af/0jeljpTK6VDZTQiUuSUwIuIPPkkdHTE24cfDuPH+xeP5JcSeBEpckrgRUS0fGR5Sb2+oakJ9uzxJxYRkT5QAi8iovr38jJpEkyeHG+3tcHy5f7FIyKSJSXwIlLedu+GZcuStymBL32pY5z6JU5EJMB8T+DNbD8z+4KZ/cnMVprZbjPbZmZLzezzZpY2RjM7yczuN7PNZrbLzJ43s6+aWUUPr/VZM1tmZq3R12g0s3/N37sTkcBraoJIJN6eOhUmTvQvHikM1cGLSBHzPYEHzgV+DcwGmoAbgHuAmcD/AL83M0s8wMw+BiwB5gB/An4ODASuB+5M9yJmdh1wC3BA9PVuA44E7jWz+Tl+TyJSLFQ+U55Sx/nJJ5O/yImIBFgQEvjXgI8CE51zn3bOfcs593+Aw4AW4Gzgk7HOZjYCLwFvB+Y65z7vnPu/wNHAU8A5ZnZe4guY2UnA5cAq4Cjn3Necc18GjgM2A9eZ2ZT8vk0RCSRdwFqeDjkEJkyIt3fuhOee8y8eEZEs+J7AO+cedc7d65zrSNn+DvDLaHNuwq5zgHHAnc65ZxL67wG+E21+MeVlLok+/8A5tyXhmLV4s/eDgM/1752ISNHZu7frTXx0B9byYKYyGhEpWr4n8L2I/T1zX8K2edHnB9P0XwLsAk4ys0EZHvNASh8RKRfPPutdxBpz0EFw8MH+xSOFlZrANzb6EoaISLYq/Q6gO2ZWCVwQbSYm3jOiz6+lHuOc22dma4D3AVOBl81sGHAQ0OqcW5/mpV6PPh+aZh9mdjFwMcCECRNo9OkXfGtrq2+vLelpTIIpm3GZdPvtTE1ovztjBi9rNZKcC+rPyrAhQzg+od3+8MM8+cADtA8Z4ltMhRTUcSlnGpNgCuK4BDaBB36EdyHr/c65vyVsHxl93tbNcbHto/rYP4lz7lfArwBmzZrl5s6d21PMedPY2Ihfry3paUyCKatxaWhIak741KeYoDHNucD+rJxyCvzXf8EbbwBQsXcvJ2/bBmed5XNghRHYcSljGpNgCuK4BLKExswuxbvo9BXgM9keHn12WR6XbX8RKWbt7bB0afI2XcBaXszg3HOTt/3hD/7EIiKShcAl8Gb2ZeBG4CXgg865zSldYjPmI0lvREq/3vr3NkMvIqVoxQrYsSPeHjsWDjvMt3DEJ5/6VHL7/vuT/78QEQmgQCXwZvZV4GfAC3jJ+ztpur0afe5Ssx6tmz8Y76LX1QDOuZ3AW0CVmR2Q5nyHRJ+71NSLSAlLt3xk8i0npBzMmgVTpsTbe/bAfff5Fo6ISCYCk8Cb2TfwbsS0Ai9539BN10ejz2em2TcHGAo86Zxry/CYs1L6iEg5SE3gtXxkeVIZjYgUoUAk8Gb2XbyLVp8FTnXObeyh+93ARuA8M5uVcI7BwPejzV+kHBNbT/7bZjY64ZgpwJeBNuC3/XkPIlJEOjrg8ceTt6n+vXypjEZEiozvq9CY2WeB/8S7s+rjwKXW9c/Ya51ztwA457ab2UV4iXyjmd2JdzfVj+ItMXk3cFfiwc65J83sp8BlwPNmdjcwEPg3YAxQF72pk4iUg5degk2b4u2RI+HII/2LR/x13HFeGc3atV67rc0ro/n3f/czKhGRbgVhBj5215QK4KvAlWkeFyYe4Jz7M3AK3o2bzgbq8G76dBlwnnOuy4oyzrnLo+d5B29d9wuAF4GPOOd+ltN3JCIF19LSQl1dHa+88gp1dXW0tLR03zm1fObkk6GiIr8BSnCZdZ2F//3v/YlFRCQDvifwzrmrnHPWy2NumuOecM79i3NutHNuiHPuSOfc9c659h5e61bn3PHOuWHOueHOuVOcc7paSaTItbS0UF1dzcKFC9m5cycLFy6kurq6+yQ+3QWsUt5SE/gHHoDt2/2JRUSkF74n8CIi/dXQ0EBrayuRSASASCRCa2srDSk3agLAOSXw0tWxx8LBB8fbsTIaEZEAUgIvIkWvqampM3mPiUQiLFu2rGvnlSth/fp4e9gwL3mT8qYyGhEpIkrgRaTozZ49m1AolLQtFApRU1PTtXPq7PtJJ0HKsVKmUpeTfPBBldGISCApgReRoldfX09VVVVnEh8KhaiqqqK+vr5rZ5XPSHeOPRamTo2329rg3nv9i0dEpBtK4EWk6IXDYZqbm6mtrWXYsGHU1tbS3NxMOBxO7ugcLF6cvE0JvMSkK6PRTZ1EJICUwItISQiHwyxYsIDDDjuMBQsWdE3eAf75T3jjjXh70CBIV2Yj5Su1jEar0YhIACmBF5HysWhRcnvePBg82J9YJBBi9w+oqanx7h8wdixMmxbvsHcv/PWv/gUoIpKGEngRKQ8dHXDHHcnbzj/fn1gkEBLvH7B8+XLv/gFHH832M85I7qgyGhEJGCXwIlIeli6FxBs7DR0KH/+4b+GI/7q7f8AvNm1K7vjgg7Btmw8RioikpwReRMrD7bcntz/2Maiq8icWCYTu7h/wx9WrYfr0+Ma9e7UajYgEihJ4ESl9e/d2LYNQ+UzZ6/b+AbNnd72YVTd1EpEAUQIvIqXvb3+DzZvj7f32g9Q6Zyk7Pd4/IHU5yb/9TWU0IhIYSuBFpPSlrj5z7rm6+6ok3T+gpqYm+f4B1dVdy2i0Go2IBESl3wGIiOTVjh1dE69Pf9qfWCRwYvcP6CJ2U6cf/jC+7fe/h898pnDBiYh0QzPwIlLa/vxn2L073p40CU46ybdwpIikltE89BBs3epLKCIiiZTAi0hpSy2fOf98GKBffZKBo46CQw6Jt1VGIyIBoU8xESld774Lf/978jaVz0imYmU0ibQajYgEgBJ4ESldv/89tLfH20ceCTNn+hePFJ/U5SRVRiMiAaAEXkRKV+rNmzT7Ltk66ig49NB4OxKBv/zFv3hERFACLyKlatUqePrp5G3nnedPLFK8zHRTJxEJHCXwIlKa7rgjuX3yyTB5sj+xSHFLrYN/8EFYudKfWEREUAIvIqXIufSrz4j0xZFHeqU0MR0d0NDgXzwiUvaUwItI6VmxAl55Jd6urOxaBiGSKTP4xjeSt916K7z1lj/xiEjZUwIvIqUndfb9zDNhv/38iUVKw6c+BVOnxtt798JPf+pfPCJS1pTAi0igtbS0UFdXR01NDXV1dbS0tPR8QHt71/p3rT4j/VVZCfX1ydt++UvYtMmfeESkrCmBF5HAamlpobq6moULF7J8+XIWLlxIdXV1j0n8qOefh7ffjm8YNgw++tECRCsl77OfhQMOiLd37YKbbvIvHhEpW0rgRSSwGhoaaG1tJRKJABCJRGhtbaWhhwsIJ6TeefUTn4ChQ/MZppSLwYPh8suTty1YADt2+BOPiJQtJfAiElhNTU2dyXtMJBJh2bJl6Q9oa2Pc4sXJ21Q+I7l08cUwenS8vWULLFzoXzwiUpaUwItIYM2ePZtQKJS0LRQKUVNTk/6A+++ncufOeHvcODjttDxGKGVn+HC49NLkbT/5CezZ4088IlKWlMCLSGDV19dTVVXVmcSHQiGqqqqoT72YMOb225Pb//Zv3sWHIrlUV+ddWxHzzjtwyy2+hSMi5UcJvIgEVjgcprm5mdraWmpqaqitraW5uZlwONy187ZtcO+9ydt08ybJh/32g0suSd7W0AD79vkTj4iUHU1NiUighcNhFixY0HvHP/0J2tri7alT4YQT8heYlLfLLvMuYN2712uvWQN33aVrLkSkIDQDLyKlIfXmTeef791BUyQfDjwQLrwweds110BHhy/hiEh5UQIvIsVvzRp49NHkbSqfkRxKe0Ox+noYkPAx+uKLXcu4RETyQAm8iBS/730veebzmGPg8MP9i0dKSrc3FBs40LtQOtEPfwjO+ROoiJQNJfAiUtyef75r+czXvuZPLFKSeryh2De/mdx52TJ47DEfohSRcqIEXkSK27e+lTTj2Tp1qspnJKd6vKHYUUfBRz6SfMAPf1jA6ESkHCmBF5HitWQJ3H9/0qbVX/gCVFT4FJCUol5vKHbFFckHPPIINDUVKDoRKUdK4EWkODnXtXzhAx9gs5aOlBzr9YZiJ5wAc+cmH3TNNYUNUkTKihJ4EfFF2lU9snHvvfDUU8nbfvQjLR0pOZfRDcVSZ+H/8hdvVRoRkTzQjZxEpOBiq3rELgxcsWIFixYt6v4uq6na27smTB/5CLz//dDYmJeYpbz1ekOx006DWbPgmWfi2665Bm67Lf/BiUjZ0Qy8iBRcj6t6ZOJ3v0ue3TTThYPiL7OuXypvv927TkNEJMeUwItIwfW4qkdv9uyBK69M3nbBBTBzZg4jFOmDj30Mjjgi3nYOPvtZ2L7dv5hEpCQpgReRgut1VY+e/OIXsG5dvD1wIFx9dY4jFOmDAQPgppuSt61dC5dd5ks4IlK6lMCLSMH1uqpHd7Ztgx/8IHnbl74EkyfnKVKRLJ16Klx6afK23/wG/vpXf+IRkZKkBF5ECi6jVT3S+clPYNOmeHv48K51xyJ++9GP4LDDkrdddBG8954/8YhIydEqNCLii15X9Uj17rvw058mb/v612HcuNwGJtJfQ4Z4F1qfeCLs2+dt27ABLr4Y/vhHLXUqIv2mGXgRKQ7f/z7s3Blvjxun2mIJrlmz4LvfTd725z/Drbf6Eo6IlBYl8CISfKtXw8KFydu++12oqvInHpFMXHEFpF6Yfeml3oWtIiL9oAReRILvu9+FxGUnp0yB2lrfwhHJSGWlV0ozZEh8244dcOGF0NHhW1giUvyUwItIsK1Y4d0QJ9F//Ze3fKRI0B16KPz4x8nbFi+G66/3Jx4RKQlK4EUk2FJXmTnqKDj/fH9iEemLL30JTj89edsVV8ALL/gTj4gUPSXwIhJcv/sdPPBA8rZrrvFumCNSLMzg5pth9Oj4tr174T/+w3sWEcmSPgVFJJiefhq+8IXkbSefDGed5U88IhloaWmhrq6Ompoa6urqaGlp8XYcdBD8938nd25uhquuKniMIlL8lMCLSPC0tMDHP548Ozl4MNx4o9bQlsBqaWmhurqahQsXsnz5chYuXEh1dXU8iT/vPO+R6Npr4cknCx+siBQ1JfAiEiy7dnnJ+7vvJm//zW/gmGN8CUkkEw0NDbS2thKJrpgUiURobW2loaEh3unnP4cDD4y3Ozrg05+Gt94qcLQiUsyUwItIcDgHn/sc/OMfydu/9S1duCqB19TU1Jm8x0QiEZYtWxbfMGYM/Pa3yQeuXQunnOL95UlEJANK4EUkOL7/ffj975O3ffSj3naRgJs9ezahUChpWygUoib1Zk6nnw7z5ydvW7XKS+J1kycRyYASeBEJhj/+Eb73veRtM2fCbbdp1RkpCvX19VRVVXUm8aFQiKqqKurr67t2/slP4BOfSN62Zo2XxK9eXYBoRaSY6VNRRHKm2xU4etPcDJ/5TPK2sWPhr3+F4cNzH6hIHoTDYZqbm6mtraWmpoba2lqam5sJh8NdOw8cCHfdBeeem7x93ToviX/99cIELSJFqdLvAEQk+FpaWmhoaKCpqYnZs2dTX1/fJSmJrcARu4hvxYoVLFq0qPsEJmbDBq9MZteu+LbKSrj7bjj44Dy9I5H8CIfDLFiwILPOoZB3l+HKSrjjjvj2N9/0kvjHHoMZM/ITqIgUNc3Ai0iPel0aLyqjFThS7d0LZ5/tzTom+vnPvQRGpER1/rXqpJO4dPRodn7yk8kd1q/3fgZeesmfAEUk0JTAi0iPMk3MM1qBI5Fz3i3mly5N3j5/Plx8cc7iFwma1C/Fv/z1r5n86KPsTC2nefddmDsX/vlPX+IUkeBSAi8iPco0Mc94BY6Ym27y1nZPdOqpcP31/Y5ZJMjSfSnevnMn3xo/Hi66KLnze+/BBz/oXSciIhKlBF5EepRpYp7xChwdHXDNNXDZZcnbp0/3lpCs1KU5Utq6+1LctHw5/PKX8MUvJh+waRPMm9f1/ggiUraUwItIjzJNzDNagWPLFu8uq1dc4SXyMSNGeCvOjBlTgHck4q8evxQPGOBdA3LppckHbd4Mc+Z4+xJ/dkSkLCmBF5EeZbM0XmwFjqamJhYsWJDc57nn4Ljj4N57kw+qqPBW4Dj88Dy/E5Fg6PVLsRnccANcfnnygTt3eteIzJvn3fhJRMqW/lYtIr3Kamm8dG6+2btgta0teft++3nL6J1+ev8CFCkisS/FDQ0NLFu2jJqamq5Ls5rBj3/M9l27GPGLXySfYPFiOPJI+OEPoa7O+xIsImVFCbyI5M/u3d6M4c03d91XUwN/+ANMmlT4uER8lsmX4pY336T6jjs4t6KCH7e3MyJx5+7d8LWveT9DN9+s9eJFyoxKaAIutlbwK6+8kt2dLUX8tno1nHRS+uT9S1+CJUuUvIv0oKGhgdadO/lVezszgQfTdXrySTj6aPjxj2HfvsIG6LM+3/m5iJXje5b0lMAHWOJawTt37uz2BjoigXPvvV69+4oVyduHDoXbbvMuxBs0yJfQRIpF4mo1LcBZwIXAjtSSmT17oL7e+8L84ou9nteviaFcJp+Z3mCulJTje5buKYEPsD7d2VLETzt2wDe/CR/9KGzdmrzv0EOhqQk+/WlfQhMpNulWq7k9FOJHn/kMfOQjXQ9YvpzIUUfx64MO4htf+ELaxC4fE0OZJOa5Tj7L8fOxHN+zdE8JfIBlfWdLkSzlbEZs61b4r/+CKVPg2mu77j/7bFi+HGbO7E+4ImWlu9VqLvnP/4S//MX7a1bK0quhjg4uevttrvzNb/j79Omsf/TRpP3ZJIG5TMxznXyW4+djOb5n6Z4S+ADL+s6WIlnIyYzYpk3w3e/C5Mnwve95a1UnqqiAn/zEu9BuxIj05xCRtHpcwtXM+2vWiy/CJz/Z5dihwOf27uWAU0+Fs86Chx4C5zJOAnOdmOc6+SyWz8dMJ0kyKWsqlvcsBeKc0yPDx3HHHecKad26dW706NEuFAq56667zoVCITd69Gi3bt26gsYh6T322GN+h9Av8+fPd6FQyAGdj1Ao5ObPn9/7we++61x9vXNVVc5B+sf++zu3ZEn+30iKYh+XUqQxybOODvfNadPcu939LMYe73ufW/TBD7rhlZUOcNddd123P/eZ/n44/vjjk/rEHjU1NX06X6YSPx9j5wra52OmMWb6WV8M77lU+fU7DHjGdZOTagY+wBJnX4YNG9bjDXR0Zbpkq08zYm+/7S1dN2UKNDRAa2vXPiNHerPyL70EJ5+c26BFpCszWs86i/dVVvJ/gXXd9XvxRc5/7DFWt7fzgwEDqNq2rds7K2f6+yHTWeFM7+gMmX2eZXODuXzIJMZM/zqRaT+/37METHeZvR7+z8An6unbn76V58a6devc/Pnz3fHHH+/mz5/f679fkGcVM3kvGc+ItbU599BDzl10kXODBnU/uzdmjHPf/75zW7cW6F2mF+RxKVcak/xL/ByoBHd+RYVbXlHR84w8uDfGjXNbv/IV55qanGtv7zxfpr8fsvn8if1eqqmp6fb3Ur4+z7L9/d7buTKJMdO/TiT2i/1VJF0/8U8QZ+B9T4qL6RHUBD7Xf5rMRi5/KfqpLx8aQU1K+vJn2y79tmxxbtEi5/7t35wbMaLnJGD8eOcaGpzbscOfN5wiqONSzjQmhZE2QX7qKec+9SnnMkjm3fjxzn32s879/veu5YUXcpqYZyofn2d9+ZKRi8mPvvTrqaxJ/KMEvsgfQU3gM/2W75w/sxDFoC8fGn78QOd0Zt0lf/B+74IL3OYrr3Ru3jznKit7/7A/8EDnbrjBuZ07C/DOM6dkMXg0JgHwxhvOff3rzo0c2fvPNjhXWel2n3ii++NJJ7nzZ850l37pSwX53Z7N51mmcv3XhExjzHUNvPgniAm8auBLQKY1iNmsOpLL+r5s5bqeP5Pz+b08Vy6Xa8vovXR0wMqVhJ9+mgXDhtG0ezdX/+//Mvrqq+HRR3u+o+OkSfDf/w2rVsFXvuLdnElEgm3SJO9urW++CQsWwKxZPffft4/BTz3FJ558kkUvvMCNN99M+Jxz4ItfhF//Gp59Ftrach5mPlZayfT3e6afaZnGmGnNejbXu4l06i6z16N4ZuAz/Zbv1yxENnI9q5/rf5vYOefPn+9+8Ytf9PhXjEz/2pHrGFP7VYI7trLS/e+ppzr3la84d/LJvZfFpD7GjXPu//wf5/7yF68mPsA02xs8GpNgeuKee5y7+Wbnzj7bueHDs/udEJ2ld9XVzn3uc84tWODc4sXOrVvn3L59fY4pH3/ZzfWKOvn867N+VoIpiDPwvifFhXwAE4GbgbeBNmAtcAMwOpPjg5rAO5dZDWKul/vKtuwkX+UfuTifn8t95eTDpa3NuZUrnXv0UbfxJz9x1wwe7G4xc8vB7c72Qzn2OPxw5775TeeefLJfH8iFpg/A4NGYBFPSuLS1OffII85ddplzM2b07XdG7BEKOTd9unMf+pBzF1/s3DXXOHfnnd6Fshs2ONfR0WNcuaypj50vXxM5uYoxRj8rwaQE3t/kfRrwbvSH8s/Aj4BHo+1XgP16O0eQE/hM+DkL4VdtYV+uD+jpF3KmFxtl80HQbYzHH+/c9u3OrVnj3DPPuJ999KPuPyoq3Hxw14K7A9xTZm7r0KHOmfXvAxecGzDAuTlznPvJT5x77bUucRYLfQAGj8YkmHocl5UrnbvxRuc+/GHnJkzo/++XxMeQIc5Nnuzc8cd75//c57z7Slx3nXP/+7/OPfCAc88+683mb9uWtDpOX/m5Ak429LMSTEFM4CspH/8NjAcudc4tiG00s58CXwN+AFziU2wFUV9fz6JFizpr/Lpbh3f27NmsWLEiqWawp/q+hoYGli1bRk1NDfX19Wnr9nqqLVywoHM4Mn7tXJ8v9n4Sj02nu1rKZ5qaYNcu2LMH2tpYt2QJUyMRBgGDgSpgWCTCqPvvhyOP9NZPjz5+sns3b5gx1Dn2g87HuGefTbp76ZejjyTOea/bFyNHwjHHeI9Zs+CMM2C//fp2LhEpLdOmwaWXeg+A9evhH/+A557znv/xD3jjjb6de/du79hMjzeD4cO934cjRyY/xx7DhsGQId5j6NAuz+EhQ1hw8cVQVwcDB3p3id60yfvv6CObzzQRv5mX4Jc2M5sKrMIrmZnmnOtI2DccWA8YMN45t7O788yaNcs988wzeY42vcbGRubOndvv87S0tPT6yyl2sWRqot+fi2pqampYvnx52u1NTU1Zv3auz5epuro6jv35zznDOUYPHUr7rl0MhuB/Ez7wwHiyHntMmeJ9MJaYXP2sSO5oTIKp3+OyaROsWBFP6FeuhDVrvO3FKBTykvlBg6Cy0mtXVnZ9JG6vqIABA7zn3v67u4dZ53+/9fbbHBQOe9sSH7F+3T2g522x/87kOZttqft629adTPv61G/lypVM/9jH4EMfyuy8OWJmzzrn0l5xHvi8I0fmRZ8fSkzeAZxzO8zsCeB04ATgkUIHV0iZzDLnYxYi17P6+fgrQSbq6+tZsXAhB0YifZ/5zgczOOAAmDzZe0ya5D1Pm+Yl6+PH+x2hiJSa/faDU0/1Hom2b4e1a71kPt1jZ7fzZP6KRLyHj/Ed5NsrS0+mg/eXowIn8D0plwR+RvT5tW72v46XwB9KiSfwmcok0c9GpuU7mb52rs+XqXA4zJizzoK//jUn5+vVkCHeh2TqIzVZnzjRmzkSEfHbiBFw1FHeI5VzsGMHbNgQf7z3XnI79ti0yfsyENSEX8RH5ZLAj4w+b+tmf2z7qNQdZnYxcDHAhAkTaGxszHVsGWltbfXttXPlrrvu4p133mHnzp0MGzaM/fffn1WrVrFq1apAnC9Th7W1MSxlmzOjIxSiY+BAXPS5Y+DAzm3tQ4bQPniw95z6iG7vGDKESFUV+0aOJDJiBJERI+gYPLj3gNat8x4ClMbPSqnRmARTIMZl1Cjvccgh3Xax9nYqdu2icudOKnbupHLXLipaW6lM2DZg714q2toYsGdP/L+jj4rYtj17sH37GBCJxJ8jEQZE/1ukJ+vXr+dVv39eEpRLAt+bWBFUlwsCnHO/An4FXg28X3WcqiENkPe9D9raeOLZZ3n/vHkweDBWWUmFGRV+xyb6WQkgjUkwaVwSOOeVz+zd692gat8+7xGJxP87dVskAu3t3o3x2tt7/2/nvOd0j+i+1195hUOmT0+/fk+sX+ojFn9322L/nclzNttS9/W2rad/+yD3A958800mnn02BwTo56VcEvjYDPvIbvaPSOkn0r1x4wCIrFzprYwgIiLFzSy+Ik1VlW9hvNXYyCEBShLFs7KxkYkBG5cBfgdQIK9Gnw/tZn/sb3fd1ciLiIiIiARCuSTwj0WfTzezpPccXUby/cBu4OlCByYiIiIiko2ySOCdc6uAh4ApdL0XztXAMOB/e1oDXkREREQkCMqlBh7gS8CTwE1mdirwMjAb+CBe6cy3fYxNRERERCQjZTEDD52z8LOAW/AS98uBacBNwInOuSK9dZyIiIiIlJNymoHHOdcCfM7vOERERERE+qpsZuBFREREREqBEngRERERkSKiBF5EREREpIgogRcRERERKSJK4EVEREREiogSeBERERGRIqIEXkRERESkiCiBFxEREREpIkrgRURERESKiBJ4EREREZEiogReRERERKSIKIEXERERESkiSuBFRERERIqIOef8jqFomNl7wBs+vfxYYKNPry3paUyCSeMSPBqTYNK4BI/GJJj8GpfJzrlx6XYogS8SZvaMc26W33FInMYkmDQuwaMxCSaNS/BoTIIpiOOiEhoRERERkSKiBF5EREREpIgogS8ev/I7AOlCYxJMGpfg0ZgEk8YleDQmwRS4cVENvIiIiIhIEdEMvIiIiIhIEVECLyIiIiJSRJTAB5iZTTSzm83sbTNrM7O1ZnaDmY32O7ZSZmbnmNkCM3vczLabmTOz23o55iQzu9/MNpvZLjN73sy+amYVhYq7lJnZfmb2BTP7k5mtNLPdZrbNzJaa2efNLO3vMo1LfpnZtWb2iJm1RMdks5k9Z2ZXmtl+3RyjMSkwM/tM9PeYM7MvdNNH45JH0c9v183jnW6O0ZgUgJmdbGb3mNn6aK613sweMrN/SdM3MGOiGviAMrNpwJPAeOAvwCtADfBB4FXg/c65Tf5FWLrMbAVQDbQCbwKHAYucc//RTf+PAfcAe4C7gM3AR4AZwN3OuXMLEHZJM7NLgF8A64HHgHXABOCTwEi8f/9zXcIvNI1L/pnZXuAfwEvABmAYcAIwC3gbOME515LQX2NSYGYWBv4JVABVwEXOuf9J6aNxyTMzWwuMAm5Is7vVOXddSn+NSQGY2XeA/8K7SdN9eJ8xY4FjgMecc/UJfYM1Js45PQL4AP4GOKAuZftPo9t/6XeMpfrA+5J0CGDA3Oi/923d9B2Bl7i0AbMStg/G+wLmgPP8fk/F/gDm4f2iHJCyfX+8ZN4BZ2tcCj4ug7vZ/oPov/F/a0x8HR8D/g6sAn4c/Tf+QkofjUthxmItsDbDvhqTwozJudF/y4eB4Wn2h4I8JiqhCSAzmwqcjvcD//OU3VcCO4HPmNmwAodWFpxzjznnXnfRn85enAOMA+50zj2TcI49wHeizS/mIcyy4px71Dl3r3OuI2X7O8Avo825Cbs0LgUQ/fdM5/fR50MStmlMCu9SvC+/n8P73EhH4xI8GpM8i5ZdXgvsAs53zu1I7eOciyQ0AzcmlYV8McnYvOjzQ2kSlh1m9gRegn8C8Eihg5MksbF6MM2+JXi/HE4ys0HOubbChVVWYr9k9yVs07j46yPR5+cTtmlMCsjMDgd+BNzonFtiZvO66apxKZxBZvYfwCS8L1TPA0ucc+0p/TQm+XcScDBwN7DFzD4MzMQrj1nmnHsqpX/gxkQJfDDNiD6/1s3+1/ES+ENRAu+3bsfKObfPzNYA7wOmAi8XMrByYGaVwAXRZuIvVo1LAZnZ1/Hqq0fi1b9/AC85+VFCN41JgUR/Ln6HV152RS/dNS6Fsz/euCRaY2afc84tTtimMcm/46PP7+Jdx3Nk4k4zWwKc45x7L7opcGOiEppgGhl93tbN/tj2UfkPRXqhsfLXj/BmTe53zv0tYbvGpbC+jlfe91W85P1B4PSEDz/QmBTS9/AuwrvQObe7l74al8L4LXAqXhI/DC9hXAhMAR4ws+qEvhqT/Bsffb4EGAKcBgzH+zz5GzAH+ENC/8CNiRL44mTRZy0hFHwaqzwxs0uBy/FWaPpMtodHnzUuOeCc2985Z3jJySfxZqGeM7NjsziNxiQHzKwGb9b9J2nKAPp0yuizxqUfnHNXR6/ledc5t8s594Jz7hK8hSmGAFdlcTqNSf/Fln00vJn2R5xzrc65F4FP4K1Ad4qZnZjh+Qo+Jkrggyn2TW5kN/tHpPQT/2isfGBmXwZuxFu+8IPOuc0pXTQuPogmJ3/CK/HbD/jfhN0akzxLKJ15DfhuhodpXPwVuwh/TsI2jUn+bYk+r3bONSfuiP7VKvYX3Zroc+DGRAl8ML0afT60m/2xlR26q5GXwul2rKIfpgfjXVy5upBBlTIz+yrwM+AFvOQ93U1QNC4+cs69gffl6n1mNja6WWOSf1V4/76HA3sSbxaEV+IE8OvothuibY2LvzZEnxNXldOY5F/s33hrN/tjCf6QlP6BGRMl8MH0WPT59NQ7TJrZcOD9wG7g6UIHJl08Gn0+M82+OcBQ4EmtFJAbZvYN4HpgBV7yvqGbrhoX/x0YfY6tsKExyb824DfdPJ6L9lkabcfKazQu/oqVaCQmfhqT/FuCl3AfYmYD0+yfGX1eG30O3pgUctF5PTJ/oBs5BeJBZjdyeo8A3dyhVB94JQEOeAYY00tfjUv+x+MwYP802wcQv5HTExqTYDzwaqy7u5GTxiW///bvS/c7C5iMt6qcA67QmBR8XG6L/lt+P2X7h4AOvNn5UUEdE4sGIAFjZtPw/qcYD/wFb1mi2Xh3CX0NOMk5t8m/CEuXmX0c+Hi0uT9wBt7syOPRbRudc19P6X833vqxd+LdXvmjRG+vDHzK6QetX8zss8AteLO5C0hfZ7jWOXdLwjEfR+OSN9FSph/jzWStAjYBE4BT8C5ifQc41Tn3UsIxH0dj4gszuwqvjOYi59z/pOz7OBqXvIn+238T76/ra4AdwDTgw3gJ4P3AJ5xzexOO+Tgak7wys/HAE8B0vM/3ZXhfqj6Bl5Cf75z7Q0L/jxOkMfH7G5Ae3T+AMN7SU+uBvcAbeBfu9Tj7qEe//92vwvvh7e6xNs0x78f7JbwFr7zpn8DXgAq/308pPDIYEwc0alwKOiYz8e4UvQLYiPfn6G3A8uh4pf09pTHxbbxiP0Nf6Ga/xiV///anAHfgrZi1Fe/mc+8BD+Pdx8I0Jr6NzRi8yoY10TxrE96k6QlBHxPNwIuIiIiIFBFdxCoiIiIiUkSUwIuIiIiIFBEl8CIiIiIiRUQJvIiIiIhIEVECLyIiIiJSRJTAi4iIiIgUESXwIiIiIiJFRAm8iIiIiEgRUQIvIhIAZnahmTkzu9DvWLJhZmPMbLOZ/dzvWPrKzOZG/+2vyvF5z46e99RcnldERAm8iEiORZO2bB4X+h1zP/wnMAT4YWyDmd0bfV9npjvAzF6L7r+lm/3/Gd3/3bxEXDh/BP4B/NTM9HkrIjlT6XcAIiIl6Oo0274KjARuBLam7FsBrAGeBtbnMa6cMrNJQC3wW+fcWwm7HgH+FTgVeDDlmIOAQwAX3Z/OvITzFC3nnDOza4G7gPOA230OSURKhBJ4EZEcc85dlbotOss+ErjBObe2m0O35S+qvKjF+xy5JWX7o9HneXQVS9rvBs41s0Occ6/HdprZMKAGaAWW5TRaf/wF7wvbl1ACLyI5oj/piYgEQHc18Ga2NvqoMrPrzazFzHab2Qoz+3i0T6WZXWFmr5vZHjNbZWbze3itM8zsfjPbaGZt0f4/NrNRWcRrwOeAFufckym7/wm8BxxtZmNS9s0DdgHXJrQTnQyEgMXOuX0Jr3eYmd0Sff9tZvaumd1uZjO6iW+omX0r+u+008xazewpM/v3LN7jYDO7OzouP4+VwZjZBDO7zsxejZ57a/S/bzGzqYnncM61AX8G3m9mh2X62iIiPVECLyISfCHgYeBf8GZ0fwdMA+6JXiB5F94MbyPwP0AVsMDM/i31RGb2PbyyltnA/wNuAlYCXweeMLMRGcb0PuAA4InUHc45BzyG9xkzN2X3PGApXm34e3Qto4m1O8tnorX0/wA+DSzHK0N6BPgksMzMjk15j6Oir/FDoB24GbgVGAfcbmbf7+3NmdlovH/zTwLfcs592TnXYWZDo+/5cuAN4BfAb/C+tHwMOCLN6WL/Rqf19roiIplQCY2ISPAdiJfAzo3O6GJmvwOWAH8AVgEznXNbo/t+CrwCfBMvuSe6/YN49flPAf8S6x/ddyHw2+j+r2UQ0weiz890s/8R4FN4Cfsfo69xCBAG/jtaH74YmGtmFk36IaX+PZpI34E3az/HOfdSQszvA5rwvrQkJvE3AMcA33DONST0H4w3G36Fmd3tnFuRLnAzmww8AEwHLnDO3Zaw+1S8L083OOe+lnLcQGBQmlMujz7PAX6W7jVFRLKhGXgRkeLw1VjyDuCcexzvwtfReInq1oR9q/FmfY80s4qEc1wafb4osX/0mFvwLqb9dIbxTIo+d3fRbbo6+Hkp+x7DmxU/EjqT9aOBDXgz2gAXAKOAKxOT92jMLwK/Bo4xsyOi59gP+A/gmcTkPdp/D/ANwIDz0wVtZkfjfcE5CDgrJXlPtDt1g3Nur3NuR5q+70SfJ6XZJyKSNc3Ai4gE31bn3Ko0298GDgaeTbPvLaAC2D/63wAnAhG8i0fPTXPMQGCcme3nnNvUS0z7RZ+3pNvpnFtpZuuAw83sAOfcerwEfntCvI9Fn+cBz+OV2wwAHk2YkT8x+lzdzTrth0afDwdeAo7He9/dreseSuif6gPAZcAOvNn+5jR9FuP9e34zWrpzP96XpRXOufY0/QE2R5/HdrNfRCQrSuBFRIKvu9Vp9gE459Ltj10AGkrYth/e7/0re3m9KqC3BD42Az24hz6P4F3oOs/Mbgc+CCyJJbrOuZfN7F28spQbSFP/TvyLwkUZxJzY//joo7f+iY4BhgNP4pUgdeGc225mJ+CVGn0UOCO6a6OZ/TfwfedcJOWwIdHnLrP2IiJ9oRIaEZHysQ3Y4pyzXh5vZHCuDdHn/Xrok1hGcyReucxjKX0agTnRUp9067/HvpxU9xLzrSn9r++l/wfTxPszvItSzwD+amZD0vTBOfemc+7zwHhgJl5p0ibge9FHqti/0YY0+0REsqYEXkSkfDwNjI5e/Nlfz0efe1oaMZaIz6Nr/XvMY8AI4CN4ZS1rnHNrEvY/HX0+OcO4lgEdWfRP5JxzX8L7a8DpwP+LrkvfbWfn3IvOuQXAh6KbP56ma+zfaEUfYhIR6UIJvIhI+bg++vxrMzswdaeZDYuWh2TicbwlGrvtH617fwWYAvwfvFrw51O6xWbk/zP6nHr31d/i3QjpSjOrSRPzADObm/CaG4BFwCwz+66ZdSkVNbNpZnZwD3F/DbgGr+Tnb4lLa5rZTDObkuawCdHnXWn2xf6NUv/6ICLSJ6qBFxEpE865R8zsm3jJ6etmdj/eSjZVwGTgFLz108/M4FzbzOwRvGUgRzvn0l7MipeQH4ZXQvNH51xHynleM7O3o/tj/RP3bzKzc4A/AU9HX/NFvFn2SXgXue5Hci3+fOAQvC8FnzGzpcC7eMtxHo5XG//v0ffe3fu7wsz24NW6P2xmZ0bf42nAT80sVie/AZiItwZ8B/DjNKc7He9LSOpfH0RE+kQz8CIiZcQ5dy3eeuT/D3g/8FXgXLxlE38FfCeL0/033so15/XQJzEh724GOrbdkSbJdc49AhwVfb0pwCXAF/Dqzx9NfX3n3Ha8LyN1wEbgbLzVZT6It8LM1/Bu0tQj59x/AvVADfCImY0F/oZXYjMYL2m/HO/f82HgZOfc3YnnMLND8Wbgb3XOpZudFxHJmsVX6hIREclc9MLTfwJ7gWOcPlC6MLOf4P1F4PDo+vwiIv2mGXgREemT6HKQXweqgU/6HE7gmNkBwBeBBUreRSSXVAMvIiJ95py738y+Qs/rwZerKcC1wI0+xyEiJUYlNCIiIiIiRUQlNCIiIiIiRUQJvIiIiIhIEVECLyIiIiJSRJTAi4iIiIgUESXwIiIiIiJFRAm8iIiIiEgR+f/mwUcs9LKNNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x648 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (12,9))\n",
    "plt.plot(sir_fit['time'], sir_fit['I'], color='red', linewidth=4, label='I')\n",
    "plt.scatter(plague_dat['time'], plague_dat['Deaths'], color='black', s=30, label='Deaths')\n",
    "plt.xlabel(\"Time (Weeks)\",fontsize = 20)\n",
    "plt.ylabel(\"Deaths\",fontsize = 20)\n",
    "plt.xticks(fontsize = 20)\n",
    "plt.yticks(fontsize = 20)\n",
    "plt.legend(fontsize = 20)\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eb6b027",
   "metadata": {},
   "source": [
    "## Question: What is R0 for your fitted model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "dabe7c23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2160544687197459"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitval_params[0]/fitval_params[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a9283b3",
   "metadata": {},
   "source": [
    "## Question: What do you notice about your fitted model?\n",
    "Look at the plot above: what do you notice about it? Where does it do a good job of fitting the data? Where does it fail to fit the data?\n",
    "\n",
    "What might we change about the assumptions of our model to improve the model fit?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "796905a2",
   "metadata": {},
   "source": [
    "## Question: What is missing from our model?\n",
    "\n",
    "The SIR model is very simplistic - can you think of any aspects of disease transmission which we could add to our model in order to make the fit more realistic?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b68016e",
   "metadata": {},
   "source": [
    "\n",
    "## Optional questions to consider\n",
    "\n",
    "* If we vary the initial guess for the parameters, does the answer change in a meaningful way\n",
    "* What happens if we do not assume that the initial population size is known at the start? How should we change the calibration?\n",
    "* What happens if we assume that not all cases are reported on, and that the data represented here are systematically under-counting, but we do not know by how much? How should we change the calibration?\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}