Week3_2-MoreSIR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #ddd;\n",
       "    color: #000000;\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 2.2em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0, 80, 120);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 1.9em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(200,100,0);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(94,127,192);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "        font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "HTML{ASCIIString}(\"<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\\n<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\\n<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\\n<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\\n<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\\n\\n<style>\\n\\n@font-face {\\n    font-family: \\\"Computer Modern\\\";\\n    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\\n}\\n\\n#notebook_panel { /* main background */\\n    background: #ddd;\\n    color: #000000;\\n}\\n\\n\\n\\n/* Formatting for header cells */\\n.text_cell_render h1 {\\n    font-family: 'Philosopher', sans-serif;\\n    font-weight: 400;\\n    font-size: 2.2em;\\n    line-height: 100%;\\n    color: rgb(0, 80, 120);\\n    margin-bottom: 0.1em;\\n    margin-top: 0.1em;\\n    display: block;\\n}\\t\\n.text_cell_render h2 {\\n    font-family: 'Philosopher', serif;\\n    font-weight: 400;\\n    font-size: 1.9em;\\n    line-height: 100%;\\n    color: rgb(200,100,0);\\n    margin-bottom: 0.1em;\\n    margin-top: 0.1em;\\n    display: block;\\n}\\t\\n\\n.text_cell_render h3 {\\n    font-family: 'Philosopher', serif;\\n    margin-top:12px;\\n    margin-bottom: 3px;\\n    font-style: italic;\\n    color: rgb(94,127,192);\\n}\\n\\n.text_cell_render h4 {\\n    font-family: 'Philosopher', serif;\\n}\\n\\n.text_cell_render h5 {\\n    font-family: 'Alegreya Sans', sans-serif;\\n    font-weight: 300;\\n    font-size: 16pt;\\n    color: grey;\\n    font-style: italic;\\n    margin-bottom: .1em;\\n    margin-top: 0.1em;\\n    display: block;\\n}\\n\\n.text_cell_render h6 {\\n    font-family: 'PT Mono', sans-serif;\\n    font-weight: 300;\\n    font-size: 10pt;\\n    color: grey;\\n    margin-bottom: 1px;\\n    margin-top: 1px;\\n}\\n\\n.CodeMirror{\\n        font-family: \\\"PT Mono\\\";\\n        font-size: 100%;\\n}\\n\\n</style>\\n\\n\")"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setting up a custom stylesheet in IJulia\n",
    "file = open(\"style.css\") # A .css file in the same folder as this notebook file\n",
    "styl = readall(file) # Read the file\n",
    "HTML(\"$styl\") # Output as HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> More on SIR models </h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>In this lecture</h2>\n",
    "\n",
    "- [Outcome](#Outcome)\n",
    "- [The threshold prediction of our SIR model](#The-threshold-prediction-of-our-SIR-model)\n",
    "- [The duration of infectiousness implied by our model](#The-duration-of-infectiousness-implied-by-our-model)\n",
    "- [A phase plane: plots of $I$ vs $S$](#A-phase-plane:-plots-of-$I$-vs-$S$)\n",
    "- [Putting it all together: estimates of $\\gamma$ and $\\lambda$](#Putting-it-all-together:-estimates-of-$\\gamma$-and-$\\lambda$)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to the top](#In-this-lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Outcome</h2>\n",
    "\n",
    "After this lecture, you will be able to\n",
    "- Explain the threshold phenomenon in SIR models\n",
    "- Use three-week duration of infectiousness to estimate $\\gamma$ in our EVD model\n",
    "- Plot the course of an epidemic in the $SI$ phase plane\n",
    "- Use country populations plus the threshold phenomenonn to estimate $\\lambda$ in our EVD model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to the top](#In-this-lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>The threshold prediction of our SIR model</h2>\n",
    "\n",
    "Recall the three equations of our SIR model:\n",
    "\n",
    "1 ...     $ S(t_{i+1}) = S(t_i) - \\lambda S(t_i)I(t_i)\\,dt$\n",
    "\n",
    "2 ...     $I(t_{i+1}) = I(t_i) + \\lambda S(t_i)I(t_i)\\,dt - \\gamma I(t_i)\\,dt$\n",
    "\n",
    "3 ...     $R(t_{i+1}) = R(t_i) + \\gamma I(t_i)\\,dt$\n",
    "\n",
    "And recall that $S(t)$ is a symbol for the size of the susceptible population (not yet infected) at time $t$, likewise $I(t)$ represents the infecteds, and $R(t)$ the removeds.\n",
    "\n",
    "The parameter $\\lambda$ models the likelihood that when an infected person meets a susceptible person, the susceptible person becomes infected.\n",
    "\n",
    "The parameter $\\gamma$ models the rate at which infected people become removed.\n",
    "\n",
    "Our task in this lecture is to propose estimates for the parameter values. To do so, we will use the *threshold phenomenon* of our model: for given values of $\\lambda$ and $\\gamma$, there is no epidemic unless $S(0)$ is big enough. \n",
    "\n",
    "The reasoning is simple: for an epidemic to happen, the infection must spread. But that means more people must be getting infected than are recovering. That is, the gain minus the loss of infecteds must be positive. In symbols, we must have\n",
    "\n",
    "$\\lambda SI - \\gamma I > 0$\n",
    "\n",
    "This can only happen if $\\lambda S - \\gamma > 0$. In other words, the number of infecteds can increase only when $S > \\gamma/\\lambda$.\n",
    "\n",
    "This implies that if the initial population is below the threshold, indicated by $S(0) < \\gamma/\\lambda$, then the number of infecteds, even if initially large, will decrease. That is, in a small enough population, any initial infection simply dies out. On the other hand, if $S(0)>\\gamma/\\lambda$, the number of infecteds will increase, at least for a while.\n",
    "\n",
    "It is of course obvious that for our model to apply to the West African epidemic, our parameter values must be such that the susceptible population at the beginning of 2014 was larger than $\\gamma/\\lambda$. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to the top](#In-this-lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>The duration of infectiousness implied by our model</h2>\n",
    "\n",
    "We note that the loss rate of infecteds is given by $\\gamma I$. This means that every day, the fraction of infected people becoming removed is $\\gamma$.\n",
    "\n",
    "For example, consider $\\gamma = 1/10$. Then every day a tenth of ill people recover. If we assume that every person is ill for 10 days, and as many people are getting ill as recovering, then a tenth of all ill people recover every day.\n",
    "\n",
    "In other words, it makes sense to thing of the duration of one person's illness as $1/\\gamma$. Recall that the duration of EVD is approximately three weeks. This means we may estimate that $1/\\gamma\\approx 21$. We will therefore set $\\gamma=0.05$ as a reasonable estimate.\n",
    "\n",
    "It is crude to think that they are equally infectious every day of these three weeks, but that is what we are doing. Our assumption was that the only difference between people is whether they are susceptible, infected or removed.\n",
    "\n",
    "Likewise, it is crude to think that all infected people are ill in the same way and for the same duration, but again this is what we do. In fact, without a lot more data (especially data on duration of illness) and a considerably more complicated model, this is all we can do. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to the top](#In-this-lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>A phase plane: plots of $I$ vs $S$</h2>\n",
    "\n",
    "A good way to illustrate the threshold phenomenon is via plots of $I$ versus $S$. That is, we run the model as we did in the previous lecture, but now we plot the phase diagram. (Well, part of it---the full phase diagram gives all three dependent variables.)\n",
    "\n",
    "We re-use the code from the previous lecture (just cut and paste):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "updateSIR (generic function with 1 method)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function updateSIR(popnvector)       # exactly the same function as before\n",
    "    susceptibles = popnvector[1];\n",
    "    infecteds    = popnvector[2]; \n",
    "    removeds     = popnvector[3];\n",
    "    newS = susceptibles - lambda*susceptibles*infecteds*dt\n",
    "    newI = infecteds + lambda*susceptibles*infecteds*dt - gam*infecteds*dt  \n",
    "    newR = removeds + gam*infecteds*dt\n",
    "    return [newS newI newR] \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots  # and of course we need to load Plots and choose to use the PyPlot backend\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# likewise, a run of the model uses exactly the same code ... but we'll play a bit with the values that determine a run\n",
    "\n",
    "gam = 1/20.       # recovery rate parameter  (ditto)\n",
    "dt = 0.5         # length of time step in days\n",
    "tfinal = 610.;    # respecting community values: lowercase only in the names \n",
    "s0 = 22000000.     # initial susceptibles, note that we use the  type Float64 from the start\n",
    "i0 = 4.          # initial infecteds; set this to 1. to  mimic an epidemic with an index case\n",
    "r0 = 0.          # not always the case, of course\n",
    "lambda = gam/(0.1*s0)\n",
    "\n",
    "# initialise \n",
    "nsteps = round(Int64, tfinal/dt)    # note the use of round() with type Int64 to ensure that nsteps is an integer\n",
    "resultvals = Array(Float64, nsteps+1, 3)  #initialise array of type Float64 to hold results\n",
    "timevec = Array(Float64, nsteps+1)        # ... ditto for time values\n",
    "resultvals[1,:] = [s0, i0, r0]  # ... and assign them to the first row\n",
    "timevec[1] = 0.                 # also Float64, of course.\n",
    "\n",
    "for step  = 1:nsteps\n",
    "    resultvals[step+1, :] = updateSIR(resultvals[step, :])  # NB! pay careful attention to the rows being used\n",
    "    timevec[step+1] = timevec[step] + dt\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVWX+B/DPc84VFAQ3zBBBcM0lDNAcFSU3dEyvIqlYmss05qRW06TtqZWVTP1qUifHickpFXdxzz0VzQ1cEzVzwcxUXFFcuOc8vz8IRgQRkHvPvfd83q8XL+Nw7j3PgU+XL8/z3OcRUkoJIioTtWrVgpeXF44cOWJ0U/KsW7cOnTt3xgcffIA333zTYdd1xu8FlVxkZCS2b9+O7Oxso5tC5FIUoxtA5MxOnjwJRVHu+aGqKq5evZp3vhACQgi7tGXdunVQFAUffvihXZ6/rNnr+1CUt99+G4qiYOvWraV63MKFC+3Usge3dOlSdOvWDTVq1ICHhweqV6+O0NBQPPfcc1i2bJndrlsWmU5ISICiKJg1a1YZtYrI+VmMbgCRK6hXrx4GDBhQ4LgQAuXLl8/7fNOmTYYUFpSjtMWAPQvjsvDOO+9gwoQJqFixInr06IHatWsjOzsbBw4cwOzZs3Hs2DF0797d6GYWyZm/v0T2wAKLqBjq1auHd999977nhYSE2K0NHM2/P3f8Hh07dgwfffQRQkJC8MMPP+Chhx7K9/WbN29i586dBrWueNzx50J0PxwiJCpDtWrVQoMGDfIdu3PYKiEhAeHh4fD29kZ0dDSAnF8+06ZNw+OPP45q1arBy8sLQUFB6NmzJ5KTkwHk9GBER0dDCJH3fLlDlL/++mup27t//3706dMHNWrUQPny5VG3bl387W9/w6VLl8rk/MLEx8dDVVV07doVWVlZRZ57+fJlfPzxx4iKikLNmjXh6emJWrVqYfDgwThx4kS+c9u2bZs3fBoZGZn3Pbr75/GgoqKiYLFY7vl9f+GFF6AoCjZu3Jh3bN68eYiKikKNGjVQoUIFBAQEoEuXLli8ePF9r7d9+3bouo7Y2NgCxRUAlC9fHm3bti1W2+8cZt68eTOeeOIJ+Pr6omrVqujbty+OHz9erOcBAJvNhk8++QTNmjWDl5cXKleujI4dO2LFihX5zhs4cCCGDRsGABgwYEDez8XDw6PY1yJyRQ7pwbp+/Tri4+OxY8cO7NixA5cuXcL06dPx7LPPlur52rdvn+/F607lypXDrVu3HqS5RKVW2DBI7vDThx9+iI0bN8JqtaJr1655v2BeffVVfPbZZ6hfvz4GDBgAb29v/Prrr9i8eTPWr1+PyMhIdOjQAadOncI333yDDh06oF27dnnP7evrW6q2bty4Ed26dYPNZkPfvn0RFBSELVu24LPPPsPy5cuxbds2VK5cudTnF+aVV17B559/jqeffhrTp0+HxVL0S9CPP/6I8ePHo0OHDoiNjYW3tzcOHjyIGTNmYOXKlUhNTUVAQAAA4E9/+hMURUFycjKGDh2KoKAgAEDVqlVL9f25l2effRbJycmYNWsWXn311Xxfy87Oxty5c1G7dm1ERUUBACZNmoSXXnoJAQEBiI2NRdWqVXHmzBns2LEDixcvRs+ePYu8XrVq1QCgTN8skJycjPfeew/dunXDiy++iAMHDmD+/PnYvHkztm/fnve9uxcpJWJiYrB8+XI88sgjGDVqFDIzMzFnzhx0794dkyZNwogRIwAAsbGxyMzMxNKlS9G7d2+EhoYCAFRVLbP7IXJK0gFOnDghhRAyODhYdujQQSqKIv/73/+W+vnWrl0rZ86cme9j2rRpUgghe/ToUYYtJ7PLzW79+vXluHHjCnxs27Yt3/m1atWS9evXz3fs7bfflkIIWalSJZmWllbgGpUqVZK1a9eWt27dKvC1S5cu5f332rVrpRBCTpgwoUT3UNjjNE2TwcHBUlVVuWHDhnznv/LKK1IIIYcPH17q86XM/73Izs6WzzzzjFQURf71r38tdtuvXLkiL1++XOg9qaoqX3jhhXzH3377bakoityyZUuxr3Hn4xYsWFCsNlWoUEE2a9aswNcWLVokhRDyrbfeyjvWrFkz6eXlJS9evFjg/MKO3S0zM1PWqlVLCiGk1WqVM2fOlD/99NN9H1eY3CwoiiK//vrrfF+bMmWKFELI3r175zseGRkpy5Url+9YQkKCFELI6OhoabPZ8o6fPHlSVqtWTXp6esr09PS841999ZVUFEXOnDmzVO0mckUOKbBu374tz549K6WUcteuXVII8UAFVmFmzJghhRBy9uzZZfq8ZG65BZaiKIV+/OMf/8h3flEF1uuvv17oNSpVqiTr168vs7Ozi2xLWRZYGzZskEII2atXrwLnX716VVapUkX6+PhITdNKdb6U//teXL9+XXbp0kUqiiI//vjjErW9KI0bN5YNGjTId8wRBZaUUvbt21cqiiL379+f73jv3r2loijy4MGDeceaNWsmK1WqJK9cuVKiNt0pNTVVNmnSRCqKIoUQUgghq1SpIq1Wq0xKSir28+RmoUmTJgW+pmmarFu3rrRYLPkK+8IKrHbt2klFUeSePXsKPM97771X4GfNAovMyCFzsMqVK1fo3IHCrFy5Eu3atUPFihXh6+uL7t274+DBg/d93MyZM1GxYkVYrdYHbS5RAV26dIGmaQU+XnzxxWI9XgiBFi1aFPq1uLg4/Pzzz2jatCnGjh2L77//Hjdv3izL5hewe/duCCHyhrHu5OPjg/DwcFy/fh0//fRTqc7Pdf36dXTo0AHr169HQkICXnvttRK3df369ejZsydq1qwJDw+PvDk8aWlpDzT/7EEMHDgQUkp8++23eceuXLmCFStWIDw8HI0aNco7HhcXh8zMTDRp0gSvvfYavvvuO2RmZpboemFhYThw4ACSk5PxwQcfICYmBhaLBcuWLUNMTAyGDBlSoueLjIwscExRFLRu3Rq6rmPfvn1FPn7Pnj3w8fFBs2bNCnytffv2kFJiz549JWoTkbtxqknu3377Lbp37w4fHx/Ex8fj3XffRVpaGtq2bYv09PR7Pi4jIwNr165FTEwMKlSo4MAWExVfjRo1Cj0+ZcoUTJw4ERaLBR988AE6dOiAatWqYejQobh48aJd2pK7dte92uTv75/vvJKef+d19u3bh2rVquXNGyuJ2bNno3Pnzti0aRPatWuHV155BWPHjsW4ceMQGBho2HzLrl27ws/PD4mJiXnH5syZg1u3bmHgwIH5zn399dfx73//G/7+/vj000/RrVs3+Pn5ITY2tsjXtcK0atUKb775JhYsWIBz585h/vz5qFixIr755hssX7682M9zr59j7vErV64U+fjMzMwSZ4HIbJymwLp+/TpeeuklDBs2DMuXL8cLL7yAV199Fdu2bYOu60Uurjh79mxomoZnnnnGgS0mKpl7rQOkqipeffVVHDhwAKdPn8asWbPQpk0bTJ8+HYMGDbJLW3x9fSGlxNmzZwv9+m+//ZZ3XmnOz+Xv74+FCxfi8uXLiIqKwrFjx0rUzrFjx8Lb2xu7d+/G7Nmz8fHHH2Ps2LF49913DX0XmsViQVxcHE6fPo0NGzYAyPkDMff43YYOHYodO3bg3LlzWLRoEXr27IlFixY9cI97TExMXi/q+vXri/24e/0cc49XqlSpyMf7+PiUOAtEZuM0BdaaNWtw5coVxMXF4cKFC3kfQgi0bNky70WsMLNmzUL16tXRqVMnB7aYqOw9/PDDiIuLw6pVqxASEoJVq1blbVGS+64rTdMe+DphYWEAgO+//77A165du4aUlBR4e3ujfv36pTr/Tl27dkVSUhIuXLiAqKgo/Pzzz8Vu57Fjx9CkSRMEBwfnO/7LL78UWKYBKNvv0f3kDhPOmDEDJ0+exNatW9G5c+cip0NUrVoVPXv2xNy5c9GuXTvs37+/REsjFKZixYolXmcqd/mPO+m6jq1bt0JRlLx3+t1LWFgYMjMzCx0GzH2tfuyxx/KOOfLnQuQsnKbA+umnnyClRPv27VG9evW8j4ceeghr1qzB+fPnC33c8ePHsW3bNsTFxUFRnOZ2iIrl5s2b2L59e4HjmZmZuH79OsqVK5f3yyl3uYFTp0498HXbtWuH4OBgLFu2rMCSJ+PHj8fly5fz1iwqzfl369KlC5YsWZJXZB09erRY7axduzaOHDmCCxcu5B27efMmhg8fXugv66pVq0JKWSbfo/tp0aIFGjRogAULFmDatGkAUGB4EEChS8pkZ2fnDf96enoWeZ3t27dj5syZuH37doGvnTt3DgkJCRBCFDqv6l7S0tLwn//8J9+xL7/8EseOHUPPnj3vu9zGoEGDIKXE66+/nu/nkJ6ejs8//xweHh54+umn84478udC5CycZiV3XdchhMCMGTMKHdu/13o5M2fOhBAi3//MRK7i+vXraNWqFRo2bIjw8HAEBQUhMzMTy5Ytw/nz5/HGG2/kFS2NGzdGjRo1MHPmTKiqioCAAAgh8PLLL8Pb27tE11UUBdOnT0e3bt0QHR2db12rTZs2oWHDhpgwYUKpzy9Mp06dsHTpUlitVkRFRWHDhg33XQR01KhReOWVV9CsWTPExsYiOzsbq1evRrly5fDoo48iLS0t3/nt27eHEAJjxozBnj17UKlSJVStWhV/+ctf7vs9KWkvEJBTUL3zzjv45JNPULFiRfTq1avAOd27d0e1atXQsmVLBAcH4/bt21i1ahUOHz6M/v37o2bNmkVe45dffsHAgQMxYsQItGvXDg0bNoSqqjh58iSWLVuGrKws9OzZEzExMcVud9euXTFixAgsXboUjRs3xv79+7Fs2TI8/PDD+L//+7/7Pn7w4MFYuHAhVqxYgdDQUHTv3h1Xr17F3LlzcfnyZXzxxRcIDAzMO79169bw9PTEp59+inPnzqF69epQFAVvvPFGsdtM5HIc/bbFey3TMG/ePCmEkGvWrCnR8zVu3LjA2+KJysqJEyekoiiyW7duxTq/Vq1aJVo64Pbt2zI+Pl526dJFBgYGyvLly0t/f3/Zvn17OW/evALnb9u2TUZFRUlfX9+8pSJOnz5dZJvWrl0rFUWRH374YYGv7du3T/bp00dWr15denp6yjp16si//e1v91yfqSTnF/a9kFLK9evXS29vbxkQECCPHDlSZNullPLLL7+UTZs2lV5eXjIgIEAOHz5cXrhwQUZGRkoPD48C53/99dcyNDRUVqhQQSqKUqzXh5Iu05DrxIkTUlVVqSiKHDx4cKHn/POf/5Q9e/aUISEh0svLS1avXl22atVKfvXVV/mWtbiXzMxMOXPmTPnss8/KRx99VFatWlV6eHjIGjVqyOjoaPnNN98Uu713LtmxadMmGRUVJX18fGSVKlVk37595fHjxws85l7fZ5vNJj/55JO873WlSpVkx44d5cqVKwu99rJly+Tjjz8uvb29paIohT4nkTsRUjp2k6iUlBS0aNGiwErumZmZCAwMRHh4OFavXl2gxyojIwN+fn75ju3Zswfh4eEYO3Ysxo4d65D2ExG5qnXr1qFz58744IMP8OabbxrdHCK35rAhwilTpuDy5cs4ffo0AGDJkiV54/EvvvgifHx88OWXX+LZZ59FeHg44uLiUL16daSnp2P58uWIjIzEF198ke85Z8yYweFBIiIicjoO68EKCQm555ovx48fz9v7atOmTfj444+xbds23Lp1CwEBAWjbti1GjhyZ904mIGe+RFBQEPz9/bFjxw5H3AIRkUtjDxaR4zisB6u4b0Vu165dsRYkFELwHSlERCWUu/k4EdmXw+dgEREREbk7LhxFREREVMbsNkSYkZGBVatWITg4mPsDEhERkdu4ceMGTpw4gS5duhRY4SCX3QqsVatWYcCAAfZ6eiIiIiJDzZgx4577INutwMrdO2zGjBlo1KiRvS5TIjdsEjvOS2w7K7H9nMTJaznHG1cBWj4k0KK6QI0KAj4WoGI5wNNy/4mg17Mlzt0Ezt+QOHsDOHdD4vzNnH9PZwE/XwWkBPzKA2F+QISfQISfghCfe2/+624GDx6M6dOnG90MMhAzQMwAuVMG0tLSMGDAgAL7pN7JbgVW7rBgo0aNEB4ebq/LFEuWTWLKjzomHtBx4RYQ4gN0aaOgc4BAh5oCVcvbr9C5clti61mJjWdyPj45L2G7BFQvD7TzF4h6WCDKX0HTqoDipgVXcHCw4RkgYzEDxAyQO2agqClQTrMXoT3c0iSmpemYsEfHhZvAnxoq+FuogvqVHFfIVPIQ+GOgwB9/35brWrbED3cUXK9u13Fb11HVE2j7sECUv0Cv2gpCfN2z2CIiIjIDtyywsnWJ6Uck3k/VcDoLGFhP4N1wFXWcoGipWE6gcy2BzrVyPr9hk9h27n8F15s7dbyyTUebGgID6gn0raPYtYeNiIiIyp7bFViHL0vErLEh7TIQV1dgbLiKRyo7b4FSwSLQvqZA+5o5n1/Pllh8UmLGUR0jt+p48QcdTwYKDKiv4MlAgfLFmBfmbH777Tejm0AGYwaIGTCH9PR0ZGRkFPq1n3/+GampqQ5uUdnw8/PL23GmuNyqwFqRrqP/eg0B3sDu3hY8Vs31ihHvcgJP1xN4up6Cs1kSc47pmHFU4qm1Gip5AH1Ccoqttg8Ll5mzpaqq0U0ggzEDxAy4v/T0dDRq1AhZWVn3PCciIsKBLSo7Xl5eSEtLK1GR5RYFlpQSE/fqeHOnju5BAjPaq/D1cI3ioyg1vARebKrixaY5PXMzj+qYcVTHV4c1BFUEnqmnYEA9BY2rOPe9/vDDD0Y3gQzGDBAz4P4yMjKQlZXlVKsHlIXcdwxmZGSYq8DKskkM3ahhzjGJd8IUjItQXKZnpyQaVhZ4r7mK8REKtp6VmHFUYmqajo/26AirBgyop6B/PQX+Xu5370RE5DqcYfUAZ+DSW+VcuS3RbqmGpekS8zqqeK+56pbF1Z2EEGjzsIIvI1X89owFSZ1V1PUVeHOXjlqzbHjyOxvWndbBLSaJiIiM47IFVraeMy/p56sSW6wWPFXHZW+l1DxUgZ7BCuZ1suC3Zyz4V6SKX7MkOq3QEL7IhplHdWTrLLSIiIjudPLkSbRv3x6VK1e2W2+bS1YlUkq8kKxh4xmJRZ1Vl5zMXtYqewo894iC1BgL1nZT8XAFgQEbNNSZbcOn+zRcvW1coTVkyBDDrk3OgRkgZoCcia+vLyZMmIDExES7XcMlC6y/79Px1WGJf7dV8URNl7wFuxFCoGOAgpV/tGBfrAWdAgTe2KkjcJYNo7drOHXN8YVWdHS0w69JzoUZIGaAjPDpp5/i+eefz/v8ypUrqF69OoQQaN26Nby8vOx2bZeb5D7/mI7Xduh4O0zBoAYsroryaFWBr6MsmNBcYtKPOqam6fh8v45+dQX+9qiKMD/H9Pz179/fIdch58UMEDNAAHDsqsTl26V7bGUPlHjB8Oeeew4NGzbE3//+d/j6+uLrr79Gr169ULly5dI1ogRcqsA6cFFi4Pca+tcVeC+CxVVx1fQW+OhxFW8+puA/R3R8tl/HzKM2dKwp8Gqogi61hGk2niYiImNk3JSoP9eG0k4NVgXw2wAL/Eqwu0mlSpXw1FNP4T//+Q9efvllfPnll5g7d27pGlBCLlNg6VLiz5s1hPgA/2mnsiAoBR8PgZeaqhjRWMHC4xJ/36fjj99paFoF+Fuoiv51BTxVfl+JiKjs+ZUX+Kmv5YF6sEpSXOUaNWoUrFYrHnnkETz00ENo1qxZ6RpQQi5TYE1N07HtnMTmHqpLbhfjTCyKQN+6An3qCGz+TeKTfTqGbNTw5k7gxSYKnm+koIpn2X2Pk5OTERkZWWbPR66HGSBmgICSD/GVhYYNG6JOnToYNmwYPvnkk3xfk1LabVkjlxhnO31d4vUdOp5/REHkwy7RZJcghEA7fwVLuliQ1seC7kEC41JzJsS/8oOGczfKJnTx8fFl8jzkupgBYgbISH/+85+haRpiY2MBADdu3EBgYCD69euXtwXOW2+9VabXdIkerFFbNXhbgI8fZ3FlL49UFpjW1oL3IyQmH9TxxQEd/z6s45VHFfztUeWBth6aPXt2GbaUXBEzQMwAGWnDhg144YUX8vbErFChAk6dOmXXazp9xZJ0QseiExJftFZRuQyHrahwNbwE3m+u4licBcMbKYjfq+etpXXTVroeLXu+DZZcAzNAzAAZ4cyZM2jUqBF2796Nl19+2aHXVgDgpZdeQkhICBRFwb59++77oMGDB0NRFFy9etWujdN0iTHbNXStJfBUCIsrR6pWXuDvLVUc7WdBbIiC13boqD/Xhq8O6bBxdXgiInIB/v7+SEtLQ3JyMry9vR16bQUA+vTpgy1btiA4OPi+D1i0aBE8PDwc8i6++cclfroKvN9c4bsGDRLgLfCvtioOPmVBmxoCf96socl8G+Yd06Fzv0MiIqJCKQAQGRmJmjVr3ncm/dmzZ/HRRx/hs88+s/tmwrqUmLBHQ3SAQPPqTj+S6fYaVBaY3dGC1BgL6voK9F2noUWSDatO3X9j6dGjRzuoleSsmAFiBshsSjTJfdiwYfj73//ukG62ZekS+y8Ck7uzuHImYX4CK7pasOmMjjd26uj6nYYof4GPWihoVaPwn1VQUJCDW0nOhhkgZsA80tLSjG5CmSrt/RS7wEpISEDt2rURFRVVqguVhJQSE3briKyRs4wAOZ92/gqSewgsT5d4a5eG1ks0WGvrmNBcRdOq+YdzR40aZVAryVkwA8QMuD8/Pz94eXlhwIABRjelzHl5ecHPz69Ejyl29bJhwwYsXrwYderUQUhICAAgNDQUe/fuLfJx3bp1g9VqzffRqlUrJCUl5Ttv9erVsFqtAIB1v0rsOC/xdpiCESNGICEhId+5qampsFqtyMjIyHd87NixmDhxYr5j6enpsFqtOHToUL7jkyZNKtBlnZWVBavViuTk5HzHExMTC90Jvl+/fkXex53c8T6EEOheW0FC8H6EpX6GfRk2hC6wYeAGG45dlS5zH7lc/efB++B98D54H0beR1BQENLS0jB58mS0bdsWKSkp+T769OmDd955J9+xGTNmoG3btli7dm2+43/+858xatSofMeWLVuGtm3bYv78+fmOjx49GgMHDsx3LDk5GW3btsVXX32V7/iECRPQo0ePAm3r3LkzPvnkk3zHcu9jwoQJaNWqFUaOHAmr1YqQkBDExcUV+F7cTcg7JtCEhIRg8eLFCA0Nve8DFUXBlStX4OPjU+jXU1NTERERgZSUFISHh9/3+e4Us9qG45kSu3tbOLndhdzWJBIO63gvVceFW8CwRxS8HabgYS/+DImIyH0Up8ZRAGD48OEIDAzE6dOn0aVLFzRo0ABATiU8bdq0Qh8ohLDLRPdzNySWpUv8qSHfOehqPFSBvzRW8XOcBe9HKJh5VEfdOTY8/905XL7Fdxya2d1/OZP5MANktgwoADB16lScOnUKt2/fxpkzZ3DkyBEAwPjx4zFs2LBCH6hpGnx9fcu8QbOO6hACeLoe5165Ki+LwGuPqTgeZ8HLTRUkHK+AkNk2/N8+DdlcQ8uUxowZY3QTyGDMAJktA05XxUw/osNaW6BaKXbMJudS2VNgQgsV2zpdQVxdBaN36Gi2wIY1v+hGN40cbPLkyUY3gQzGDJDZMuBUBdbuDIm9F4HB9Z2qWfSAmtevhS8jVaTGWOBXXiB6pYbea2w4kcneLLPgW/SJGSCzZcCpKpnpR3TUqAB0DWTvlTtqVk1gY3cVs9qr2HFeotE8G8alaLhRyj0OiYiInJXTFFiaLjH7mI5n6imwKCyw3JUQAv3rKTjUx4K/Pqrgoz06Gs2zYcHx+68IT0RE5CqcpsDaeV7i3A0gJpjFlbu5e00WAKhYTuDDFioOPGVB0yoCT63V0HmFhoOXWGS5o8IyQObCDJDZMuA0BdbSdIlqnkCrh1hguZusrKx7fq1+JYFlXS1Y1kXFyWsSzRbY8MoPGq7cZqHlTorKAJkDM0Bmy0C+hUbLUkkXGg1dkI3Hqgp8075E2yOSG7mlSXy2X8cHu3V4lwM+bqFiUAMBheuhERGREyn2QqNGO5mZs7Fzj9pO0RwyiKcq8PpjKg71saBjTYGhm3L2ONx5nss6EBGRa3GKimZpuo5yCtClFnsqCKhVUWBWBws2dVdx0ybRMknDc5tsOHeDw4ZEROQanKLAWpYuEeUv4OvBAssd3b0ZaXG19VewK8aCyW0ULDwh0WCuDf84wNXgXVFpM0Dugxkgs2XA8AIrW5fY/JtEdACLK3c1dOjQUj/Wogi80FjFkb4WxNVR8NcfdIQttGHDrxw2dCUPkgFyD8wAmS0DhhdYqRkSWTYgyp8FlrsaN27cAz+HX3mBqW1V7IqxoJKHQIflGvqutSH9GnuzXEFZZIBcGzNAZsuA4QXWpjMS3hYgzI8FlrsqzrtIi/1cfgLJPVR8+4SK5LMSj8y14f1UDTe5GrxTK8sMkGtiBshsGTC8wNp4RqJ1DYFyXL2dikkIgQH1FRzuY8HIJgre362jyXwb1p7msCERETkHQwssTZdIPivR7mEWV1RyPh4C8S1V7I+1IKiiQOcVGgZ/b8OFm+zNIiIiYxlaYO2/BFy5DbTj/Cu3lpCQYNfnb1hZYP2TKr5qq2LxyZxNpBOPcm9DZ2LvDJDzYwbIbBkwtMDa8lvO+lePV2eB5c5SU1Ptfg0hBP70iIK0PhY84S/w9AYN3VdpnATvJByRAXJuzACZLQOGFli7MiRCqwqUt7DAcmdTpkxx2LUe9hKY28mCJdEq9l2UaDwvZ+0sjWtnGcqRGSDnxAyQ2TJgaIGVkiERwXcPkh30qK3gx6csGNwgZ+2s1ks07LvAIouIiBzDsALrhk3i4CWwwCK78fUQmNxGRbJVxbVsiYhFNry1k0s6EBGR/RlWYO29IKFJFlhkf61rKEjtbcE74Qo+2acjdIENG89wSQciIrIfwwqslAyJcgrQtKpRLSBHsVqtRjcBnqrAu+Eq9vS2oEYFgSeWafjzJhsu3WJvliM4QwbIWMwAmS0DhhZYTavk/OIj9zZy5Eijm5AMVvNOAAAgAElEQVSnURWBjT1UfNlGwdxjOUs6zD/GJR3szZkyQMZgBshsGTCswNp3EQirxuLKDKKjo41uQj6KEBjeWMXBPha0ekigzzoNvdZo+IVLOtiNs2WAHI8ZILNlwJACS5cSaZclmlRhgUXGCfAWWBRtwYJOKnael2g834Z/HtSgszeLiIgekCEFVvo1IMsGNGaBRU6gd4iCg09Z8HRdBSO26Gi7VMPBSyyyiIio9AwpsHJ/ebHAMoekpCSjm3BflT0FprZVsbG7igs3JR5baMO4FA23NBZaZcEVMkD2xQyQ2TJgTIF1WaJiOSDQ24irk6MlJiYa3YRia+evYE9vC15rpmDCbh1hC23Y8huXdHhQrpQBsg9mgMyWAcN6sBpVFhCCPVhmMGfOHKObUCLlLQLvN1exu7cFlTwEIpdqeCFZw5Xb7M0qLVfLAJU9ZoDMlgGDerCARpWNuDJR8TWtKpDcQ8Wk1gq+Paqj8TwblpxkbxYREd2fIQXWkSsSDSux94qcn6oIjGyi4uBTFoRVE+i5WsOg7224zAVKiYioCA4vsC7dkrh0C6jrywKLXEdgRYGlXVRMj1KRdEKi6QIbVp1ibxYRERXO4QXWz1dz/vKv6+voK5NRhgwZYnQTyoQQAoMaKDjwlAWNKwt0/U7D85s1ZHJu1n25Swao9JgBMlsGDCiwcv5lD5Z5uNvqvYEVBVb9UcXUSAUzj+ZsHv39r+zNKoq7ZYBKjhkgs2XAkB6sKp5AFU8WWGbRv39/o5tQ5oQQeL6Rin2xFgRVFGi/XMNLWzVk2dibVRh3zACVDDNAZsuAwwuso1cl6rH3itxEHV+BDd1VfPYHBdMO6XhsgQ1bz7I3i4jI7BxeYB3LBOr4OPqqRPajCIGXH1Wxp7cF1coLtF2q4bXtGm6yN4uIyLQcXmClX5MI9mEPlpkkJycb3QSHaFhZYHMPFROaK/j8gI6IRTaknGeRBZgnA3RvzACZLQMOLbB0KfHLdW6RYzbx8fFGN8FhLIrA64+pSImxwFMFWi62YWyKhtsm39PQTBmgwjEDZLYMOLTAOnsDyNaBoIrswTKT2bNnG90Eh2taVWB7LwveDlPw4W4df1hsw/6L5i2yzJgByo8ZILNlQAGAl156CSEhIVAUBfv27Sv0xAMHDiAqKgqNGzdGaGgonnvuOdy6datEF0u/lvMLJpAFlql4eXkZ3QRDlFMExkWo2NbTgts6ELHIho/2aLDp5iu0zJoB+h9mgMyWAQUA+vTpgy1btiA4OPieJ5YvXx5TpkzBwYMHsXfvXly7dg0TJ04s0cVOXcv5N4hDhGQiEdUFUmIseOVRBW/v0hG5VMPhy+YrsoiIzEQBgMjISNSsWRNS3vtFv169emjatCmAnDWAWrRogRMnTpToYqeuS3hbgCqepW8wkSvyVAU+flxFcg8VF29KPLbQhs/2a9CL+H+OiIhcV6nmYF2/fh1fffUVevXqVaLH/XIdCPDOKdDIPEaPHm10E5xGqxoK9sRa8HwjBa9s09F+mYZjV92/yGIGiBkgs2WgxAVWdnY24uLi0LVrV1it1hI99rcsiYcrsLgym6CgIKOb4FS8LAKft1Kx4UkV6dckQhfYMPWgVmQPsqtjBogZILNloEQFls1mQ79+/RAQEIDPPvusWI/p1q0brFYrrFYrVm/bi6MpyWjVqhWSkpLynbd69epCC7YRI0YgISEh37HU1FRYrVZkZGTkOz527NgC88LS09NhtVpx6NChfMcnTZpUoJrOysqC1WotsFZHYmJioZtU9uvXj/dRjPsYNWqUW9zHncriPs5snIM2a0dhQD0Ff9mio8tKDaeuSZe7j+L8PHIz4Or3kYv3UfL7GDVqlFvcB+AePw8j7uPQoUMueR+JiYl5dYzVakVISAji4uIKPMfdhLzjz+aQkBAsXrwYoaGhBU7UNA19+/ZF1apV8e9///u+T5yamoqIiAikpKQgPDwcAPDo/Gy0r6ngi9bqfR9PZCarTun402YNmbeBf7RWMai+4FA6EZGTKqzGuZsCAMOHD0dgYCBOnz6NLl26oEGDBgByKshp06YBAObMmYOkpCTs2rULYWFhCA8Pz9czURxnbwA1KjzILRG5py6BCg7EWtArWGDIRg3W1RrOZLnvkCERkbuzAMDUqVML/eL48ePz/vvpp5/G008/XeoL2XSJjJtADc7BMp1Dhw7hkUceMboZTq+yp8B/n7Cgd7COYckams634Z9tVPSr6/AdrcocM0DMAJktAw575T5/E5BgD5YZjRkzxugmuJSewQp+fMqCjjUF4tZreHaDDVdvu3ZvFjNAzACZLQMOK7DO3sj5lwWW+UyePNnoJrgcv/ICczqq+OYJFUknc9bN2npWN7pZpcYMEDNAZsuA4wqs3+eTPOzFIUKzMdtbc8uKEAID6yvY29sCfy+Btks1jE1xza12mAFiBshsGXBYgfXb7z1YD5V31BWJ3EOIr8DG7irGhiuYsDtnq52fTbA4KRGRK3PgEKFEJQ+gvIU9WEQlZVEE3g3P2Wrn/I2cIcPpR3S3XpyUiMiVOXQOFudfmVNJNwWne/tDDQV7elvwVEjOcg5912m4eNP5iyxmgJgBMlsGHNqDxSUazCkrK8voJrgVHw+Br6MsmNtRxbpfJUIX2rD+tHNPgGcGiBkgs2Ug30ruZenuVU47r7Chkgcwv5PFHpcjMqVfrkkM2qhhw68Sr4YqeL+5Ak+Vf8gQEdlTsVdyd4SLtySqefKFn6gs1aoosKabiviWCj4/oOMPi21Iu+T8Q4ZERO7OYQXWpVtAFU9HXY3IPBQh8Gqoiu09LbilAeGLbPjnQY0T4ImIDOTYAsvDUVcjZ3L3LulkH2F+ArtiLBjaUMGILTp6rNJw7oZzFFnMADEDZLYMOKTA0qXEldtAFQ4RmtLQoUONboJpeFkEprRRsTRaxY7zEo8usGFFuvET4JkBYgbIbBlwSIF15XbOPoQcIjSncePGGd0E0+leW8H+WAua+wk8uUrDyC0abtiM681iBogZILNlwCEF1qVbOf+ywDKne73DguyrhpfAsi4qJrdWkHBYR/NFNuy5YEyRxQwQM0Bmy4BjCywPDhESOZIQAiOaqEiJsaCcAjyeZMMn+zTonABPRGRXjimwbue8mFdlDxaRIRpXEdjey4KXmioYvV1H9AoNp6+zyCIishcOEZLdJSQkGN0EAuCpCvy9pYo13VSkXc6ZAL/guGMmwDMDxAyQ2TLgsAJLAPDlMg2mlJqaanQT6A6dAhTsi7WgQ02Bp9Zq+NNGG65l27c3ixkgZoDMlgEHFVgSlT1zFkQk85kyZYrRTaC7VCsvMK+jioR2KuYck3hsoQ3bz9mvN4sZIGaAzJYBB83B4iKjRM5GCIGhDRXs6W2BX3mBNks0vJ+qwaZzbhYR0YNySIF1+RZQmfOviJxSvUoCm3uoePMxBeNSdTyxTMOJTBZZREQPwiEFVma2hG85Dg8SOatyisB7zVVs6q7il+s5Q4aOmgBPROSOHFRgARXLOeJK5IysVqvRTaBiavNwzpBh54CcCfB/SS6bFeCZAWIGyGwZcFiB5cMCy7RGjhxpdBOoBCp7CsztqOJfkSqmH9HRcrENaZcerMhiBogZILNlwIEFFocIzSo6OtroJlAJCSEwrJGCnb0s0HQgYpENCYd0yFKuAM8MEDNAZsuAQwqsa9mSPVhELqhpVYGdMRYMqC/w3GYNT2/QcOU2J8ATEd0P52ARUZG8LALT2lowu4OKFekS4Qtt2HmeE+CJiIrCOVhkd0lJSUY3gcpAv7oKdv++ZlbrxRo+LcGm0cwAMQNktgzYvcCSUuIa52CZWmJiotFNoDJSxzdnzay/Pqrg1e06uq/ScO7G/YssZoCYATJbBuxeYN3UAAn2YJnZnDlzjG4ClSEPVSC+pYqVXVXsOp+zZtb600UPGTIDxAyQ2TJg9wLrenbOvyywiNxL10AFe2MtaFxZoNMKDW/v5DY7RES57F5gZdly/uUkdyL34+8lsOqPKj5oruDjvTnb7KRfY5FFROSwAotzsIjck6oIvBmmYmN3Fad+32ZnEbfZISKTs3uBdUPL+WvWmz1YpjVkyBCjm0AOkLvNzhP+Ar3Xahi5RcPN37fZYQaIGSCzZcDuBdYtLeffCqq9r0TOymyr95pZFU+BBZ1UTGmj4KvDOdvsHLosmQFiBsh0GbB/gfX7SEEFi72vRM6qf//+RjeBHEgIgRcaq9je04JbWs42O7ci+pV6mx1yD3wdILNlwP7LNPw+B4s9WETm0qyaQEqMBXF1BIZs1DDwew2Z3GaHiEzCcUOE7MEiMh3vcgIJURbMbK9iyUmJ8EU2pJxnkUVE7s/+PVi6RDkFsCh8F6FZJScnG90EMljQb1uRGmNBJQ+BVkts+Hy/xiFDk+HrAJktAwoAvPTSSwgJCYGiKNi3b989T162bBkaNWqEhg0b4qmnnsK1a9fue4FbGocHzS4+Pt7oJpDB4uPjUa+SwFarilFNFPx1mw7rag0ZN1lkmQVfB8hsGVAAoE+fPtiyZQuCg4PveeL169fx3HPPYcmSJTh8+DD8/f3x3nvv3fcCt2wcHjS72bNnG90EMlhuBjxUgU//oGJZFxXbzkk0W2DD979yzSwz4OsAmS0DCgBERkaiZs2aRXbZr1y5EuHh4ahfvz4A4IUXXijWxo03dfZgmZ2Xl5fRTSCD3Z2BJ4Ny1sxqUEmgw3INY1O4zY674+sAmS0DxZ6DlZ6ejtq1a+d9HhwcjN9++w26XvRfn7c0wIs9WER0lwBvgbXdVIyPUPDBbh0dlmv4hdvsEJGbcMgyDRUsnOBORAWpisA74Sq+f1LF8UyJZgttWJHOIUMicn3FLrCCgoJw4sSJvM+PHz8Of39/KErRT7FizTocP/wjrFZr3kerVq2QlJSU77zVq1fDarUWePyIESOQkJCQ71hqaiqsVisyMjLyHR87diwmTpyY71h6ejqsVisOHTqU7/ikSZMwevTofMeysrJgtVoLvNMhMTGx0CX++/Xrx/soxn3c2T5Xvo878T5Kdh+5z3Wv+5j9wSi8mjUHrR4SeHKVhjd2aNiR4nz3kcvVfx5G3Mfo0aPd4j4A9/h5GHEfYWFhLnkfiYmJ+WqYkJAQxMXFFXiOAuQdgoOD5d69e2VhMjMzZY0aNeThw4ellFKOHDlSjh49utBzpZQyJSVFApCd/rlddl6efc/zyP198cUXRjeBDFbcDGi6LifusUn137dl2yXZ8pdrup1bRo7C1wFypwzk1jgpKSn3PEdIKeXw4cOxfPlynD17FtWqVYOPjw+OHDmCsWPHIiAgAMOGDQOQs0zD6NGjoWkamjZtiv/+97/w8fEptHBLTU1FREQEIidtR9UG4VgczYlYRFQ8yb/piFuv4bYGzGivIrqW3WczEBEVW26Nk5KSgvDw8ELPsQDA1KlTC/3i+PHj833evXt3dO/evUSN4CR3IiqpyIcV7I4RGPi9hq4rNbwVJjEuXIHKBYuJyEXYf5I7FxololKoXkFgRVcV7zdX8OEeHZ1XaPgti+8yJCLXYP+9CHW+i9Ds7p6QSOZT2gwoQuCtMBXruqlIuyzx2EIb1p/muwxdEV8HyGwZsH+BZWMPltmNGTPG6CaQwR40A0/UzFmYtGkVgU4rNLyXqkHjwqQuha8DZLYM2L3Auq0DniywTG3y5MlGN4EMVhYZqOElsOqPKsaGKxiXoqPrdxrO3WCR5Sr4OkBmy4DdC6xsHfDgG4BMLSgoyOgmkMHKKgOqIjA2QsWabir2XcwZMtx0hkOGroCvA2S2DNi99LFJwIM9WERUhjoG/G8vw/bLNXy0R4NexF6qRESOZv8eLI09WERU9vy9cvYyfKOZgrd26njyOw0ZN1lkEZFzsH+BJVlgmd3d2x2Q+dgrAxZF4IMWKlZ2VbErQyJsoQ1bfuOQoTPi6wCZLQOOmYPFIUJTy8rKMroJZDB7Z6BLoILdMRYEVxSIWqbhk30aJIcMnQpfB8hsGRDSTq9CucvI463t+M/A5hjSkN1YRGRf2brEO7t0TNyro0eQwPQoFVXLcx0+IipbxdkqxyFVD3uwiMgRyikCHz+uYlkXFVvOSoQtsmH7OQ4ZEpHjOabAYucVETnQk0EKdve2oKaXQNulGj7fzyFDInIsFlhkdxkZGUY3gQxmRAaCKgps7K5iVBMFf92mI3athsu3WGQZha8DZLYMcIiQ7G7o0KFGN4EMZlQGPFSBT/+gIqmzig2/SoQvsmHXeQ4ZGoGvA2S2DLAHi+xu3LhxRjeBDGZ0BnoGK0iNscCvvECbJRqm/MghQ0czOgNkPLNlgAUW2d293mFB5uEMGQjxFdjcQ8XwRgpGbtXRb52Gq7dZZDmKM2SAjGW2DHCIkIhMw1MV+EdrFfM7qVj1i0TEIhv2XGCRRURlz0E9WFyHhoicR2yIgtTeFviUA/6w2IZ/pXHIkIjKFocIye4SEhKMbgIZzBkzUNdXYKvVgqENFAxP1jFgg4br2Syy7MUZM0COZbYMcIiQ7C41NdXoJpDBnDUD5S0C/4xUkdhBxeKTEo8n2XDoMosse3DWDJDjmC0D7MEiu5syZYrRTSCDOXsG4uoq2NnLAgmgRZIN845xKYey5uwZIPszWwbYg0VEBKBRFYEdvSzoHiTQd52Gl3/QcFtjbxYRlQ57sIiIflexnMCs9iomtVbwz4M6nlim4ZdrLLKIqORYYBER3UEIgZFNVGzqruLU9ZzV39ed5pAhEZWMQ0qfciywTM1qtRrdBDKYK2bgDzVyVn9/rJpA9EoNE3Zr0LmUQ6m5YgaobJktAyywyO5GjhxpdBPIYK6ageoVBFZ2VfF2mIK3d+mwrtJwiRtGl4qrZoDKjtky4JDSh+uMmlt0dLTRTSCDuXIGVEVgfISKFV1V/HBOInyhDSnnWWSVlCtngMqG2TJg9wJLETlzGoiIXNkfA/+3YXTrJTZMS9O5+jsR3ZPdCyyVw4NE5CZq+wgkW1X8qaGC55M1DNmoIcvGIouICrJ/gcXOK9NLSkoyuglkMHfKgKeas/r7t0+omHtMotViG366wiLrftwpA1Q6ZsuA/Qsse1+AnF5iYqLRTSCDuWMGBtRXsKOXBTc1oPkiGxYe51IORXHHDFDJmC0D9p+DxSFC05szZ47RTSCDuWsGmlYV2NnLguhaArFrNby6TUO2zt6swrhrBqj4zJYBDhESET0AXw+BuR1VfPYHBf84oKPjcg1nslhkEZkdCywiogckhMDLj6r4vruKn69KhC20YeMZDhkSmRkLLCKiMtLm4ZylHBpXFui4XEP8Xo1LORCZFAsssrshQ4YY3QQymJkyUMNLYHU3FWOaKXhth46YNRouc/V3U2WACme2DHAdLLI7s63eSwWZLQMWReDDFiqWRKvYeEaieZINey6Yu8gyWwaoILNlgMs0kN3179/f6CaQwcyagR61FaTEWOBbDmi12IavD5t3XpZZM0D/Y7YMOGSrHCIis6rjK7DVasGAegJDN2l4bpMNN7j6O5Hbyyuwjh49ijZt2qBhw4Zo2bIl0tLSCn3AxIkT0aRJE4SFhaF169bYuXNnkRfgECERmV15i8C/21nwn3YqZh6VaLPEhmNXWWQRubO88uf555/H8OHDcfjwYYwZMwaDBg0qcPLevXvx5ZdfYteuXdi9ezdGjBiBkSNHFnkBTnKn5ORko5tABmMGcgxpqOCHnhZczQYiFtmw9KR5hgyZATJbBhQAOH/+PFJSUvDMM88AAGJjY3Hq1CkcO3Ys38lCCNhsNmRmZgIALl++jMDAwCIvwAKL4uPjjW4CGYwZ+J/Hqgns6mVBlL+AdbWGN3Zo0Eyw+jszQGbLgAUATp06BX9/fyh37GsTFBSE9PR01KlTJ+9YaGgoXn75ZYSEhKBatWrw9PTEpk2birwACyyaPXu20U0ggzED+VX2FFjUWcXf9+l4Y6eOXRkSs9qrqF7BfV8wmQEyWwZKNEPqxIkTWLhwIY4dO4b09HS8/PLL6Nu3b5GPYYFFXl5eRjeBDMYMFCSEwJhmKtZ2U7H3gkTEIht2nnffIUNmgMyWAQUAAgMDcebMGej6//7nTk9PR1BQUL6TFyxYgNDQUNSoUQNAzqJhW7Zsgc1mu+cFfjywD1arNd9Hq1atkJSUlO+81atXw2q1Fnj8iBEjkJCQkO9YamoqrFYrMjIy8h0fO3YsJk6cmO9Yeno6rFYrDh06lO/4pEmTMHr06HzHsrKyYLVaC4wTJyYmFrpAWr9+/XgfvA/eB+/jge6jfc2cpRz8vQQiF9vQbOQXLnkfd3Llnwfvg/dx930kJibmq2FCQkIQFxdX4DnuJuTv+zh06NABgwYNwqBBgzB//nzEx8djx44d+U5etGgR3n33XWzbtg3e3t6YM2cOxo0bV+g7DlNTUxEREYEWn2/Hjpcev29DiIjM7JYm8fIPOqam6fhTQ4HJrVWUt3AIgMgZ5dY4KSkpCA8PL/ScvCHCqVOn4l//+hcaNmyI+Ph4TJ8+HUBOFTlt2jQAQExMDKxWK5o3b46wsDBMmjQJs2bNKrIRFi7TYHp3/2VB5sMM3J+nKvBlpIqvo1TMOCoRuVTDyUz3mfzODJDZMmDJ/Y8GDRpg69atBU4YP358vs8nTJiACRMmFPsCXGiU7h5qJvNhBopvcAMFoVUFYtfYELHIhsQOKjrXcv2/VJkBMlsGuNkz2d2oUaOMbgIZjBkomXA/gV0xFjSvLtBlpYYPd2vQpWv3ZjEDZLYMsMAiInJC1coLLO+i4u0wBW/t0tF7jYYrt127yCIyExZYREROSlUE3muuYkm0iu/PSLRYZMOBiyyyiFwBCyyyu7vfUkvmwww8mB61FezqZUF5C9BysQ1zfna99bKYATJbBuxeYHGSO40ZM8boJpDBmIEHV6+SwA9WC3rVFohbr+GVHzRku9AWO8wAmS0Ddi+wuEwDTZ482egmkMGYgbLhXU5gRnsVX7RSMOlHHR2Xa/gtyzWKLGaAzJYBu5c/7MAis701lwpiBsqOEAKjmqrY0F3FT1ckwhfZsPWs8w8ZMgNktgxwDhYRkQuKfFhBam8L6voIRC3VMPlHDdLFl3IgcicO6MFihUVEZA/+XgLru6sY0UTBqK06nv1eQ5aNRRaRM7B/gcX6yvTu3rCTzIcZsJ9yisDnrVTMaq9i4QmJVottOHrF+YosZoDMlgG+i5DsLisry+gmkMGYAfvrX0/Btp4WZNmA5kk2LDvpXPOymAEyWwaEtNOgfe5O030TdmDO0Bb2uAQREd3l8i2JZ7/XsDRd4t1wBe+GKVD5ly5RmcqtcVJSUhAeHl7oOezBIiJyI5U9BZKiVXzQXMH7qTq6r9Jw8abzDRkSuTsWWEREbkYRAm+FqVjZVcWO8xLNk2zYncEii8iRuAwo2V1GRobRTSCDMQPG6BKoICXGgiqeQOslNvz3iHHzspgBMlsG2INFdjd06FCjm0AGYwaME+wjkNzDgv51BQZv1PBCsoZbmuN7s5gBMlsG7F9g2fsC5PTGjRtndBPIYMyAsSpYBBLaqfhXpIqEwzqilmk4fd2xRRYzQGbLANfBIru71zssyDyYAeMJITCskYLNPVT8cl0iYpENm884bsiQGSCzZYBDhEREJvL4QwpSelnQsJJAh+XcYofIXrjZMxGRydTwElj75P+22Bm8UcMNbrFDVKY4REh2l5CQYHQTyGDMgPPJ3WLn2ydUzDsmEbnUhpOZ9iuymAEyWwY4REh2l5qaanQTyGDMgPMaUF/B1p4WXLwJRCyyYf1p+8zLYgbIbBngECHZ3ZQpU4xuAhmMGXBuj1UT2BVjQZifQOeVGj7ZV/bzspgBMlsGuEwDERGhWnmB77qqGB2qYPR2Hf3Xa7iezXlZRKXFIUIiIgIAqIrAx4+rmNtRxbJ0iVZLbPj5KossotJgBxMREeXTp46CbT0tuGEDmi+yYeUp47bYIXJV7MEiu7NarUY3gQzGDLieplUFdvayoE0NgSe/0zBhtwb9AeZlMQNktgywwCK7GzlypNFNIIMxA66psqfAki4q3glX8PYuHbFrNFy9Xboiixkgs2WA62CR3UVHRxvdBDIYM+C6FCEwPkLF4mgV63+VaLnYhkOXS15kMQNktgzwXYRERHRf1toKdvSyQAB4PMmGxSc4L4uoKFwHi4iIiqVhZYHtPS3oHCDQa42Gd3Zp0HS+y5CoMBwiJLtLSkoyuglkMGbAffh4CMzvpOLDFgom7NbRY7WGS7fuX2QxA2S2DHCIkOwuMTHR6CaQwZgB9yKEwBuPqVjRVcW2cxItkmzYf7HoIosZILNlgO8iJLubM2eO0U0ggzED7qlroIJdvSzwtgB/WGzD3J/vPS+LGSCzZYBDhEREVGp1fAW2Wi3oWVug33oNY7ZrsHFeFhEs9r4Ae7CIiNybdzmBme1VtKiuY/R2HakZErM7qvArz18AZF4OeBch/wcjInJ3Qgj89VEVa7qp2HtRovkiG1Iz2JNF5sVlGsjuhgwZYnQTyGDMgHm0r6kgJcYCv/ICbZbY8O1POfOymAEyWwY4REh2Z7bVe6kgZsBcgioKbO6h4oUtGp79XsOu8xIdOncxullkMLO9DuT1YB09ehRt2rRBw4YN0bJlS6SlpRX6gFOnTsFqteKRRx5B06ZNMWXKlCIvwEnu1L9/f6ObQAZjBsyngkXgP+1UTG6t4J8HdST4PoXzNzhkaGZmex3IK7Cef/55DB8+HIcPH8aYMWMwaNCgQh8QExODwYMH49ChQzhw4AD69u1bvAsQEZGpCCEwoomKdU+qSLss0TyJ87LIPBQAOH/+PFJSUvDMM88AAGJjY3Hq1CkcO3Ys38nr1q1D+fLl0bt377xj1atXL/IC7MEiIjK3dv4562U99Pu8rBk/cR9Dcn8KkDPs5+/vD0X5X39TUFAQ0kyM3xEAAB6mSURBVNPT85188OBB+Pn5oX///ggPD0dsbCyOHz9e5AVYYFFycrLRTSCDMQN0cs8WbOqhol8dgYHfa3jlB66XZTZmex0o0QiezWbDhg0bMHbsWKSmpiI6Ovq+Q4Ssryg+Pt7oJpDBmAGKj49HBYvA11Eq/tFKwRc/6uiyUkPGTRZZZmG21wEFAAIDA3HmzBno+v+6bdPT0xEUFJTv5KCgIISFheGRRx4BAAwcOBC7d++Gpmn3vED8xHhYrdZ8H61atSqw6ePq1athtVoLPH7EiBFISEjIdyw1NRVWqxUZGRn5jo8dOxYTJ07Mdyw9PR1WqxWHDh3Kd3zSpEkYPXp0vmNZWVmwWq0FquzExMRC317ar18/3kcx7mP27NlucR934n2U7D5yM+Dq95GL91Hy+5g9ezYmTZqEMWPG4MWmKtZ2U7HvokTEwmxEDXjRZe4DcI+fhxH34e/v75L3kZiYmK+GCQkJQVxcXIHnuJuQUkoA6NChAwYNGoRBgwZh/vz5iI+Px44dOwo0ODQ0FJs2bULNmjUxb948vPfee9i/f3+BJ05NTUVERAQ+TNqJN3o2v29DiIjIXE5mSsSsseHQZSChnYr+9fi2KHINuTVOSkoKwsPDCz0nbx2sqVOnYvDgwfjwww9RqVIlTJ8+HUBOFRkQEIBhw4bBy8sLU6dOxZNPPgkAqFSpUr7eicJwiJCIiApT20dgi9WCYZs1PL1BQ+oFiY9aKLBwAUVyA3kFVoMGDbB169YCJ4wfPz7f5506dcLu3bvt3zIiInJ7FSwC3zyhItwvZx/DvRckZndQUZX7GJKLY38s2d3dY+NkPswAFZWB3H0MV/1RRWpGznpZ+y5w8ru7MdvrgP33IuQfIaZ395slyHyYASpOBjoGKNgVY4FvOaDVEhvm/sz1styJ2V4H2INFdjdq1Cijm0AGYwaouBkI9hHY2tMCa22Bfus1vL5Dg8b1styC2V4H7L7ZMzuwiIioJLwsArPaq4jw0/HaDh17LkgkdlBRxZO/Uch1cIiQiIicjhACr4aq+K6rih3nJVok2fDjRfZkkevgECHZ3d2LwpH5MANU2gx0rpWzj6GXBWi52IaFxzkvy1WZ7XXA/j1Y9r4AOb0xY8YY3QQyGDNAD5KBOr4CP1gt6BYoELtWw9s7NeiSvVmuxmyvA+zBIrubPHmy0U0ggzED9KAZ8C4nMKejio9aKPhwjw7rKg2Xb7HIciVmex3gHCyyO7O9NZcKYgaoLDIghMDrj6lY3lXFlrMSjyfZkHaJRZarMNvrAHuwiIjIpfwxUMHOXhZ4qDnzshaf4Lwscj6cg0VERC6nXqWceVmdAwR6rdEwLoXzssi5cIiQ7G7ixIlGN4EMxgyQPTLg4yEwr5OK9yMUjE/VEbtGQ+ZtFlnOymyvAxwiJLvLysoyuglkMGaA7JUBRQi8Ha5icbSKdb9KtFpiw89XWWQ5I7O9Dggp7dOnmpqaioiICHy2bCdefrK5PS5BRESU5+AliZ6rbbhwC5jbUUWnAPYhkH3k1jgpKSkIDw8v9Bymj4iI3ELjKgI7elnQwk+gy0oNn+3XYKc+BKL7YoFFRERuo4qnwIquKv72qIJXtukYvFHDTRuLLHI8vouQ7C4jI8PoJpDBmAFyZAZURSC+pYpvn1Ax55hE1DINp6+zyDKa2V4HWGCR3Q0dOtToJpDBmAEyIgMD6itI7qHi1yyJ5ots+OEs18syktleB7hMA9nduHHjjG4CGYwZIKMy0Lx6zqKkdXwFnlim4evDLLKMYrbXAc7BIru71zssyDyYATIyAw97Cax/UsWgBgJDN2l4cauGbJ1Dho5mttcBi70vwA4sIiIymqcq8K9IFWHVdLy4VceBixJzO6nwK8/fUmQf7MEiIiJTEELgL41VrH1Sxf5LEi0W2bDvAnuyyD5YYJHdJSQkGN0EMhgzQM6UgSh/Bbt6WVDZE2i1xIb5xzgvyxGcKQOOwEnuZHepqalGN4EMxgyQs2Wgto9Acg8LugcJ9Fmn4Z1d3Cza3pwtA/bGOVhkd1OmTDG6CWQwZoCcMQPe5QRmd1DxWDUdb+3Use+ixLdPqPD14G8ue3DGDNgThwiJiMi0hBB44zEVS7uo+P5XiT8stuGnK+zJogfHIUIiIjK9/2/v7qOqKvM9gH+fvY+mpOJrSgqC8WJGKiiikaJYqBkH1CZpvDOo06g3ra6tgdudmdKZbq7BXmbN0lZZMdHLhI5MHhBHE1MTzERg0lsihkqHTA0sx0FQO/s89w9GEgEFYZ99ztnfz1r+4XEfzm/nd21/7f07zzMzQMH+JAs0CYyzOfBhFeeyqGO4kjsRERGA4b0F9idaMGGgwAMfanjpEDeLppvHR4SkO6vVanQJZDBmgDwlA71vEdgcryJ1pIJf7Xfi57s11HOz6E7hKRnoLGywSHfLli0zugQyGDNAnpQBVRH4wzgV709R8bcTEhM3a/i6lk1WR3lSBjoDZ7BId/Hx8UaXQAZjBsgTM/BIsIJCqwXf1kuMtTmw9zTnsjrCEzPQEZzBIiIiakVkf4HiWRaE+gpM2aLhjSNssqht+IiQiIjoOm7rLrDjARW/CFOwqEDDsr3cLJpujHewSHc2m83oEshgzAB5ega6qgKv3qvi1RgF68qcuP/vGqrr2WS1h6dnoL3YYJHusrKyjC6BDMYMkLdkYMkIFTtnqjj8vUSUzYHPuFl0m3lLBtqKjwhJdxs2bDC6BDIYM0DelIGJfgqKZ1nQrxtwT44Dfz3Guay28KYMtAUbLCIionYK6CFQkGBBUqDA3J0afnOAm0VTU/pv9sxnhERE5IV8LAJ/mdKwWfTTRQ2bRb83RYUvN4smcAaLiIjopgkhkDZKxZbpKgpOS0TbHDh6jney6KoGq6KiAjExMQgLC0N0dDTKysqu+8b58+dDURScP39e9yLJsy1YsMDoEshgzAB5ewZm+CsoSmx4KDQux4Gt3Cy6GW/PwLUaG6zFixdjyZIlKC8vR1paGlJSUlp906ZNm9C1a1eINjz/4yNCMtvqvdQcM0BmyEBob4H9SRbcO1Bg5jYNqw9ys+irmSEDVxNSSlldXY2QkBB89913UJSGnsvPzw979+7FsGHDmrzhzJkzSEhIwK5du9CzZ0+cO3cOvXr1avaDS0tLMWbMGLy1/QDm3z/WJSdDRERkNM0p8WyJE6s+c+KROwTenKTCx8K7Dd7kSo9TUlKCyMjIFo9RAKCqqgp+fn6NzRUABAQEwG63N3vDokWL8MILL+DWW2/VqWwiIiLPpSoCz0epWB+nwlYpMXGzA3ZuFm067Rpyz8jIwNChQxEbG9vm97BnJyIiM5p7h4JPEi04exGIsjlQyM2iTUUBAH9/f5w6dQpO549/+Xa7HQEBAU0O3rVrF3JycjBs2DAEBQUBAEaOHImDBw+2+gHLn3oKVqu1ya8JEyY0WzJ/+/btsFqtzd6/dOlSZGRkNHmttLQUVqsVNTU1TV5fsWIF0tPTm7xmt9thtVpx5MiRJq+vWbMGqampTV6rq6uD1WpFYWFhk9ezsrJaHM6bO3cuz6MN53F1HZ58HlfjebTvPK68x9PP4wqeR/vPo7Cw0CvOA2jf30fZ9vWYsPMp3NlbIG6LhnVlmkeeR2f8fcyePdsjzyMrK6tJDxMUFITk5ORmP6MZ+W9TpkyRmZmZUkopN27cKKOiouSNCCHk+fPnW/yzkpISCUC+nX/ghj+HvFtCQoLRJZDBmAEyewYua065tNAh8fpluaTAIS9rTqNLcjlvysCVHqekpKTVYxoXGn3ttdcwf/58rFq1Cr6+vsjMzGzsIgcPHoxFixY1a86EEPyGBN3Q+vXrjS6BDMYMkNkz0EURWBujYlQ/gaV7NZSdk8i+T0X/buYZpDFbBhobrNDQUHzyySfNDvjd737X6ps1TbvhB3CZBvLx8TG6BDIYM0DMQINfDlcw3BeYvUND1CYHcqdZcHdfc/xDabYMcCV3IiIiF5rop6A4yYJeXYEJOQ7YKjn87o242TMREZGLDe0psNdqwXR/gVn5Gv63lIuSehsX3MHiPSyzu/bbHWQ+zAAxA8316CLw16kqVkYqeKbEieSdGuoc3ttkmS0DvINFurt2uQ8yH2aAmIGWKUJgxRgV2fepyLNL3JvrQJWXLkpqtgwIqdM9ySvLyL+3oxjzpo7R4yOIiIi8xsGzEtbtDlzUgE33q7hnIO+BuKs2b5WjJz4gJCIiurFR/QQOJFkQ5iswOU/Dn8s5/O7J9G+w2GERERG1yW3dBXY8oGJ+qMAv9mhYvk+Dw+mdjwy9He9gke6u3daAzIcZIGag7bqqAuvuVbH2HgVrvnDigW0avr/k+U2W2TLAB7yku7S0NKNLIIMxA8QMtI8QAkvvUrF9hoqSGolomwNHznl2k2W2DLDBIt2tXbvW6BLIYMwAMQM3J26wgqIkC7ooQLTNga1VnjuXZbYMsMEi3Zntq7nUHDNAzMDNu6OXwL5EC2L9BGZu0/DiIc9clNRsGWCDRURE5OZ6dRWwxav4n9EKUvc7kbJbw0UvXpTUG7DBIiIi8gCKEHg+SsX7U1RsPCERm6fhmwtsstwVv0VIuktPTze6BDIYM0DMQOd5JFhBQYKKk3USUTYHDlR7xlyW2TLAO1iku7q6OqNLIIMxA8QMdK6xAxQcSLLAv4fAxM0a/lLh/k2W2TKg+1Y56z8qxtw4bpVDRETU2S46JBYXanjnS4n/HqXg+bEKVIXPjvTWlq1yLC6uiYiIiDpJN4tAZqyKUX2dSC1y4vPvJN6PU9GrK5sso/ERIRERkQcTQuCpkSq2TFNReEZifI4DFf/k8LvR9G+w2ESbXk1NjdElkMGYAWIG9DfdX8H+RAs0CYzLceCjk+41l2W2DPAOFulu4cKFRpdABmMGiBlwjbDeAp8mWjBugMC0rRrWfO4+i5KaLQNssEh3K1euNLoEMhgzQMyA6/S5RSBvmoonwxU8sc+JRQUaLmvGN1lmywAbLNJda9+wIPNgBogZcC2LIvDSeBVvxap450uJqVs0fFtvbJNltgywwSIiIvJS80MV7H5QxZfnGxYl/eys8XeyzIIruRMREXmxCQMVFCdZMKCbQEyuA9nH3Wv43VvxDhbpLiMjw+gSyGDMADEDxhrSQ2BPggrrUIGffKRhRYkGp4uH382WATZYpLvS0lKjSyCDMQPEDBjPxyLw/hQVz49V8FypEw9/pOHCD65rssyWAd23ytmwsxgPT+FWOURERO4ip9KJebs0hPgCufENexpS27VlqxzewSIiIjKZxEAFn1gt+P4SEGVzYN8ZzmV1NjZYREREJjSyn0BRkgUhvQQm52l490s2WZ2J3yIkIiIyqdu6C+yYqWJesMDPd2t4ukiD5uRSDp2Bd7BId1ar1egSyGDMADED7usWVSBjkoqXxyt44ZATs/I1/Oty5zdZZssAGyzS3bJly4wugQzGDBAz4N6EEFh+t4q8aSo+PiVxT64DJ853bpNltgywwSLdxcfHG10CGYwZIGbAM8zwV/BpogX1WsPw+55TnTeXZbYMsMEiIiKiRnf2EdifaMGofgJTt2h48wiH32+G/kPunHInIiLyKP26CWyboWLRnQp+WaDhv/ZpcHD4vV14B4t0Z7PZjC6BDMYMEDPgebooAq/EqHglRsHaL5yYuU3DuUs332SZLQNssEh3WVlZRpdABmMGiBnwXI+NULF9hooDNRLjcxw4eu7mmiyzZYANFuluw4YNRpdABmMGiBnwbHGDFRQlWiAEEJ3jQP7X7Z/LMlsG2GARERHRDQX7CnyaaMH42wRmbNOw9gsNOm1n7BUaG6yKigrExMQgLCwM0dHRKCsra3bw559/jtjYWIwYMQIjR47Eo48+ikuXLl33AzjjTkRE5B18uwrkTVPxZLiCxz9x4j8LnfiBw+8tamywFi9ejCVLlqC8vBxpaWlISUlpdnC3bt3wyiuv4PDhwzh48CBqa2uRnp7u0oKJiIjIOKoi8NJ4FRmTVPz5qBPxf9dw9iKbrGspAFBdXY2SkhLMmzcPADBnzhxUVVXh+PHjTQ4ODg5GeHg4gIZVX6OiolBZWenaisnjLFiwwOgSyGDMADED3mdhmIKdM1V8/r3EOJsDX3x3/SbLbBlQAKCqqgp+fn5QlB9HsgICAmC321t944ULF/Dmm28iKSlJ/yrJo5lt9V5qjhkgZsA73TtIwYEkC27tAkzIdWCLvfXhd7Nl4KaG3H/44QckJydj+vTpptu8kdrvkUceMboEMhgzQMyA9wrsKbA3wYIptwskfKjhxUMtD7+bLQMKAPj7++PUqVNwOn/sPO12OwICApq9weFwYO7cuRg8eDD++Mc/3vADFi9ZDKvV2uTXhAkTmi04tn379habtaVLlyIjI6PJa6WlpbBaraipqWny+ooVK5rNhNntdlitVhw5cqTJ62vWrEFqamqT1+rq6mC1WlFYWNjk9aysrBZvbc6dO5fnwfPgefA8eB48D9OfR8+uApvuV/H0aAWp+50IWJGPS1rTJssTzuNqV/4+srKymvQwQUFBSE5ObvYzriXkv9vMuLg4pKSkICUlBdnZ2Vi9ejWKioqaHKxpGh5++GH07dsXb7zxxnV/cGlpKcaMGYPs3cWYEzvmhoUQERGR5/tLhRO/2KNhTH+BD+5TMdDH+9YTuNLjlJSUIDIyssVjGh8Rvvbaa1i3bh3CwsKwevVqZGZmAmjoIl9//XUADYuE2Ww2FBcXIyIiApGRkXj88cf1PxPyaNf+HwSZDzNAzIB5zAtW8PGDKo6flxiX48DBsw13ssyWgcY7WJ2Nd7DoCqvVitzcXKPLIAMxA8QMmM/XtRJJ+RrKzkm8N1nFW08meU0G2nUHi0gv69evN7oEMhgzQMyA+QzpIbAnQcWDAQKzd2iI/O1fTbXyu8XoAsj7+fj4GF0CGYwZIGbAnHwsAuvjVIT3ceLZEqC8VsOfJ6nobvG+uaxr6X4Hy/v/ExIREVFrhBB4JlLFxqkqciolJm3W8M0F77+TxUeEREREpLuHhinYa7XgdL1ElM2BA9WtL0rqDdhgke6uXZ+EzIcZIGaAUlNTEdFf4ECSBQE9BCZt1pBV4b1NFhss0l1LC9aSuTADxAzQlQwM8hHYNVPFT4IEfrpLwzPFGpxeOPyu+zINf9tdjNlcpoGIiIiuIqXEC4eceLrIiaRAgXcmq+jRxTMmt91imQbP+E9FREREriSEQNooFTnxKvJPStyb68BX//KeO1l8REhERESGSRiqYJ/Vgn9eBsblOPDJGe+Yy2KDRbq7dmNOMh9mgJgBul4GwvsKFCVZMNxXYEqehsyjnt9kscEi3aWlpRldAhmMGSBmgG6UgQHdBfIfUPHzEIEFH2tI3a9Bc3ruI0Ou5E66W7t2rdElkMGYAWIGqC0Z6KoKvD5Rxd19nVj+qROHv5fIilPRq6vnTXTzDhbpjl/PJmaAmAFqawaEEHgiXMXW6Sr2npG4J9eB4+c9704Wv0VIREREbid+iIJPEy24qAHROQ4UnPKsuSzewSIiIiK3NLy3wP5EC8L7CEz9u4a3yj2nyWKDRbpLT083ugQyGDNAzADdbAb6dRPY/oCKBaEKFu7xnOF3DrmT7urq6owugQzGDBAzQB3JQBdF4LV7FdzVB1j+qRNHzkm8P0VFTzceftd9q5xNu4uRxK1yiIiIqBNsq3Ji7kcaAnoAm6dZENjT9U0Wt8ohIiIirzLdX8G+RAvqHMA4mwN7T7vnXBZnsIiIiMijjOgjsD/JghF9BOK2aHjbDVd+Z4NFuqupqTG6BDIYM0DMAHV2Bvp3E9g+Q8XPQgTmf6zh6SINTn2mnm4KGyzS3cKFC40ugQzGDBAzQHpkoKsq8MZEFS+PV/DCISdm52uo/cE9miw2WKS7lStXGl0CGYwZIGaA9MqAEALL71aRG69i5zcSMbkO2GuNb7L0b7A45W56rX3DgsyDGSBmgPTOwMwABfusFvzrByDK5sC+M8bOZfEOFhEREXmFu/o2rPwe2ktgcp6G9740rslig0VEREReY0B3gR0zVcwLFvjZbg2/PmDM8DsbLNJdRkaG0SWQwZgBYgbIlRm4RRXImKTihWgFf/jMiYd2uH74nQ0W6a60tNToEshgzAAxA+TqDAgh8KuRKnLiVeSflJi42YEqFw6/675Vju3jYiRO4lY5REREZIz/+04i4UMHLmpATryK6Ns6dn+JW+UQERGR6d3dV6AoyYI7egnE5mnIqtB/+J2PCImIiMjr3dZdYOdMFXOHCfx0l4ZnivUdfrfo9pOJiIiI3MgtqkBmrIq7+jjxdJETZeck3o5VcWuXzn/exjtYpDur1Wp0CWQwZoCYAXKXDAghkDZKxab7VWyrkpiU58DXOgy/s8Ei3S1btszoEshgzAAxA+RuGUgMVLDXakF1PTAux4Gibzt3LotD7qS7+Ph4o0sggzEDxAyQO2ZgVL+G4fehPRqG3zcc67wmi3ewiIiIyLQG+QjsmqliTpBA8k4NK0s6Z/idQ+5ERERkat0sAu9OVjGitxO/KW4Yfn8rVoWP5eafw/EOFunOZrMZXQIZjBkgZoDcPQNCCPw6QsXf7lORZ5eI3azh5IWbv5PFBot0l56ebnQJZDBmgJgB8pQMzA5SUJhgwel6iXE2B0qqb67JamywKioqEBMTg7CwMERHR6OsrKzFN+Tl5eHOO+9EWFgYHnroIdTW1l73AwSn3E1vwIABRpdABmMGiBkgT8pARP+G4fchtwpM3OzAxuPtH35vbLAWL16MJUuWoLy8HGlpaUhJSWl28IULF/Doo48iNzcX5eXl8PPzw+9///uOnQURERGRm/HzEdj9oIrEQIGHP9LwXKmG9mzfrABAdXU1SkpKMG/ePADAnDlzUFVVhePHjzc5eOvWrYiMjERISAgA4LHHHkNWVlZnnQsRERGR2+huEXh/iornxih4tsSJn+7SUO9oW5OlAEBVVRX8/PygKD+OZAUEBMButzc52G63Y+jQoY2/DwwMxOnTp+F06r9pIhEREZGrCSHw20gVG6eqyKmUiM3TUF1/4yZLt2Ua6uvrAQAV5UdQeisHscysqKgIpaWlRpdBBmIGiBkgT8/AMABvDJVY/qmG5G+OAPix12mJBQD8/f1x6tQpOJ3OxrtYdrsdAQEBTQ4OCAhAfn5+4+9PnDjR7M7XFZWVlQCA5b/8jw6dEHmHMWPGGF0CGYwZIGaAvC0DlZWViImJafHPhPz3xFZcXBxSUlKQkpKC7OxsrF69GkVFRU0Orq2tRXBwMPbs2YPQ0FA8/vjj6N69O1avXt3sB9fU1ODDDz9EYGAgunfvrsNpEREREblefX09KisrMW3aNPTv37/FYxobrKNHj2L+/Pk4e/YsfH19kZmZiREjRmDFihUYPHgwFi1aBKBhmYbU1FRomobw8HC8/fbb6Nmzp+vOioiIiMjNNTZYRERERNQ5uJI7ERERUSdjg0VERETUyTrcYOm1xQ55jrZk4KuvvoLFYkFkZCQiIiIQGRmJEydOGFAt6eHJJ59EUFAQFEXBoUOHWj2O1wHv1ZYM8Drg3S5duoRZs2Zh+PDhiIiIwLRp03Ds2LEWjzXFtUB2UFxcnHznnXeklFJmZ2fLqKioZsfU1tbKgQMHyqNHj0oppVy2bJlMTU3t6EeTm2hLBiorK2WfPn1cXRq5SEFBgTx58qQMCgqSBw8ebPEYXge8W1sywOuAd7t48aLcunVr4+/Xrl0rJ0+e3Ow4s1wLOtRgffvtt9LX11dqmtb42qBBg+SxY8eaHLdx40Y5Y8aMxt8fPnxYDhkypCMfTW6irRmorKyUvXv3dnV55GKBgYGt/uPK64A5XC8DvA6YS3FxsQwKCmr2ulmuBR16RMgtdqitGQCAuro6REVFYezYsXjuuefatWkmeT5eBwjgdcBM/vSnPyEpKanZ62a5FnDInVzi9ttvx8mTJ3HgwAHs2LEDBQUFeOmll4wui4hciNcB81i1ahWOHTuGVatWGV2KYTrUYF29xc4VrW2xc2XrHOD6W+yQZ2lrBrp06dK42m3v3r2xcOFCFBQUuLRWMhavA8TrgDm8+OKLsNls2LZtG7p169bsz81yLejQ2QwYMACRkZF49913AQDZ2dnw9/fHsGHDmhw3ffp0/OMf/8DRo0cBAK+++iqSk5M78tHkJtqagerqajgcDgAN3zT54IMPEBER4fJ6yTi8DhCvA97v5Zdfxvr165Gfn9/qLi+muRZ0dIirvLxcTpgwQYaGhsqoqCj5xRdfSCmlfPbZZ+W6desaj9u8ebMcPny4DAkJkbNmzZLnz5/v6EeTm2hLBj744AMZHh4uR48eLcPDw+UTTzwhL1++bGTZ1IkWL14shwwZIrt06SIHDRokQ0JCpJS8DphJWzLA64B3+/rrr6UQQgYHB8uIiAg5evRoOX78eCmlOa8F3CqHiIiIqJN51wNPIiIiIjfw/xuZavR42fl9AAAAAElFTkSuQmCC\" />"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# now we can plot the outcome  \n",
    "\n",
    "svals = resultvals[:,1];  # get the results\n",
    "ivals = resultvals[:,2];\n",
    "\n",
    "plot(svals, ivals, title = \"First look at I vs S plot\")        # and plot them\n",
    "\n",
    "# we first vary S(0); the initial threshold is (0.1)/(1/20000) = 2000\n",
    "# one can also vary gam and lambda and dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to the top](#In-this-lecture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Putting it all together: estimates of $\\gamma$ and $\\lambda$</h2>\n",
    "\n",
    "We already have the estimate $\\gamma = 0.05$, based on a 20-day infectious period per patient.\n",
    "\n",
    "We now need estimates of $\\lambda$ and of $S(0)$ (we will discuss $I(0)$ and $R(0)$ later).\n",
    "\n",
    "Since at the beginning of the epidemic, one could not safely say that any person in the three countries concerned would never be exposed to the Ebola virus, we should assume that the susceptible population was the sum of the populations of Guinea, Liberia and Sierra Leone. Using Wikipedia, this gives us $S(0)=22\\times10^6$, because in total there are approximately 22 million people in these countries.\n",
    "\n",
    "Now we use the threshold, which we write as $S^\\ast = \\gamma/\\lambda$, so that $\\lambda = \\gamma/S^\\ast$. Assuming that our $S(0)$ is far larger than $S^\\ast$, because the epidemic was so severe, let us set $S^\\ast = 0.1S(0)=2.2\\times10^6$. This gives the estimate $\\lambda = 0.05/2.2\\times10^6 \\approx 2.3\\times10^{-8}$. \n",
    "\n",
    "Let's go back to our phase plot and put in these estimates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to the top](#In-this-lecture)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.6",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}