week08d-us-population-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Modeling exponential growth\n",
    "\n",
    "The goal of this exerice is to find the parameters of an exponential fit of US population growth with time\n",
    "\n",
    "$$ y = a_1 \\exp(a_2t) $$\n",
    "\n",
    "where $y$ is population and $t$ is time in years.\n",
    "\n",
    "Data for this exercise is located in a text file at [data/uspopulation](data/uspopulation)\n",
    "\n",
    "This exercise is based on Chapter 38 of Shapiro, B. E. (2016) Scientific Computation: Python Hacking for Math Junkies, Sherwood Forest Publishing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, read in the data and put into variables that are easy to remember and type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/uspopulation/us-population.txt',header=2)\n",
    "yr = df['YR']\n",
    "pop = df['USPOP']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this example is focused on the US, create a new variable where the year is referenced to the founding of the US."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = yr - 1776\n",
    "logpop = np.log(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'population')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEGCAYAAABYV4NmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbrElEQVR4nO3df5BV5Z3n8fcnYNxeoyKKLjYwkEjY0TUl8QbNspNyK7tAUjOBMGZDdicy0Skylm50Z5YamFStRmtXHSe642R1QlbXH5v4YwwiU8YhrOikdkfRRoiIDEsnotKwyIgoqTBG8Lt/nOeWp9vbze2+fe659/bnVXXq3vu955z7vacbvn2e5znPUURgZmY22j5UdgJmZtaZXGDMzKwQLjBmZlYIFxgzMyuEC4yZmRVifNkJtIrTTjstpk+fXnYaZmZtZdOmTX8fEZNqvecCk0yfPp2enp6y0zAzayuSXhnsPTeRmZlZIVxgzMysEC4wZmZWCBcYMzMrRGEFRtJUSU9K2i5pm6SrUvxaSX2StqTl87ltVkrqlbRD0vxc/HxJW9N7t0lSih8v6cEU3yhpem6bpZJ2pmVpUd/TzMxqK3IU2RHgDyPieUknApskrU/v3RoRf5pfWdLZwBLgHOBM4H9J+nhEHAXuAJYBzwA/AhYAjwOXAW9GxFmSlgA3AV+WNBG4BqgAkT57bUS8WeD3NTNrK2s293Hzuh3sOXiYMyd0sXz+LBbN7h61/Rd2BhMReyPi+fT8ELAdGCrzhcADEfFORLwM9AJzJE0GToqIpyOb+vleYFFum3vS84eBz6azm/nA+og4kIrKerKiZGZmZMVl5eqt9B08TAB9Bw+zcvVW1mzuG7XPaEofTGq6mg1sTKErJb0g6S5Jp6RYN/BabrPdKdadng+M99smIo4AbwGnDrGvgXktk9QjqWf//v0j/4JmZm3m5nU7OPzu0X6xw+8e5eZ1O0btMwovMJI+AvwQuDoi3iZr7voYcB6wF/h2ddUam8cQ8ZFu834gYlVEVCKiMmlSzQtRzcw60p6Dh4cVH4lCC4yk48iKy/cjYjVAROyLiKMR8R7wPWBOWn03MDW3+RRgT4pPqRHvt42k8cDJwIEh9mVmZsCZE7qGFR+JIkeRCbgT2B4Rt+Tik3OrfRF4MT1fCyxJI8NmADOBZyNiL3BI0oVpn5cAj+a2qY4QuxjYkPpp1gHzJJ2SmuDmpZiZmQHL58+i67hx/WJdx41j+fxZo/YZRY4imwt8FdgqaUuK/THwFUnnkTVZ7QK+DhAR2yQ9BLxENgLtijSCDOBy4G6gi2z02OMpfidwn6ResjOXJWlfByRdDzyX1rsuIg4U8i3NzNpQdbRYkaPIlP3Bb5VKJTzZpZnZ8EjaFBGVWu/5Sn4zMyuEC4yZmRXCBcbMzArhAmNmZoVwgTEzs0K4wJiZWSFcYMzMrBAuMGZmVggXGDMzK4QLjJmZFaLIucjMzKxgRd+VshEuMGZmbap6V8rqjcOqd6UEWqLIuInMzKxNNeOulI1wgTEza1PNuCtlI1xgzMzaVDPuStkIFxgzszbVjLtSNsKd/GZmbaoZd6VshAuMmVkbWzS7u2UKykBuIjMzs0K4wJiZWSFcYMzMrBAuMGZmVggXGDMzK4QLjJmZFcIFxszMCuECY2ZmhXCBMTOzQrjAmJlZIVxgzMysEJ6LzMysZK182+NGuMCYmZWo1W973IjCmsgkTZX0pKTtkrZJuirFJ0paL2lnejwlt81KSb2Sdkian4ufL2lreu82SUrx4yU9mOIbJU3PbbM0fcZOSUuL+p5mZo1o9dseN6LIPpgjwB9GxK8DFwJXSDobWAE8EREzgSfSa9J7S4BzgAXA7ZKqd9K5A1gGzEzLghS/DHgzIs4CbgVuSvuaCFwDXADMAa7JFzIzs1bR6rc9bkRhBSYi9kbE8+n5IWA70A0sBO5Jq90DLErPFwIPRMQ7EfEy0AvMkTQZOCkino6IAO4dsE11Xw8Dn01nN/OB9RFxICLeBNbzflEyM2sZrX7b40Y0ZRRZarqaDWwEzoiIvZAVIeD0tFo38Fpus90p1p2eD4z32yYijgBvAacOsS8zs5bS6rc9bkThnfySPgL8ELg6It5O3Sc1V60RiyHiI90mn9sysqY3pk2bNlheZmaFafXbHjei0AIj6Tiy4vL9iFidwvskTY6Ivan56/UU3w1MzW0+BdiT4lNqxPPb7JY0HjgZOJDiFw3Y5qmB+UXEKmAVQKVS+UABMjNrhla+7XEjihxFJuBOYHtE3JJ7ay1QHdW1FHg0F1+SRobNIOvMfzY1ox2SdGHa5yUDtqnu62JgQ+qnWQfMk3RK6tyfl2JmZtYkRZ7BzAW+CmyVtCXF/hi4EXhI0mXAq8CXACJim6SHgJfIRqBdERHVsXuXA3cDXcDjaYGsgN0nqZfszGVJ2tcBSdcDz6X1rouIAwV9TzMzq0HZH/xWqVSip6en7DTMzNqKpE0RUan1nuciMzOzQrjAmJlZIVxgzMysEC4wZmZWCBcYMzMrhAuMmZkVwgXGzMwK4QJjZmaF8B0tzcwa1Km3PG6UC4yZWQM6+ZbHjXITmZlZAzr5lseNcoExM2tAJ9/yuFEuMGZmDejkWx43ygXGzKwBnXzL40a5k9/MrAGdfMvjRrnAmJk1qFNvedwoN5GZmVkhXGDMzKwQLjBmZlYIFxgzMyuEC4yZmRXCBcbMzArhAmNmZoVwgTEzs0K4wJiZWSFcYMzMrBAuMGZmVggXGDMzK4QLjJmZFcIFxszMClH3dP2SxgFn5LeJiFeLSMrMzNpfXQVG0r8HrgH2Ae+lcACfKCgvM7OmWrO5zzcNG2X1NpFdBcyKiHMi4ty0DFlcJN0l6XVJL+Zi10rqk7QlLZ/PvbdSUq+kHZLm5+LnS9qa3rtNklL8eEkPpvhGSdNz2yyVtDMtS+v8jmY2Rq3Z3MfK1VvpO3iYAPoOHmbl6q2s2dxXdmptrd4C8xrw1jD3fTewoEb81og4Ly0/ApB0NrAEOCdtc3tqkgO4A1gGzExLdZ+XAW9GxFnArcBNaV8Tyc62LgDmANdIOmWYuZvZGHLzuh0cfvdov9jhd49y87odJWXUGertg/k58JSkx4B3qsGIuGWwDSLiJ/mzimNYCDwQEe8AL0vqBeZI2gWcFBFPA0i6F1gEPJ62uTZt/zDwnXR2Mx9YHxEH0jbryYrS/XXmYmZjzJ6Dh4cVt/rUewbzKrAe+DBwYm4ZiSslvZCa0KpnFt1kZ0lVu1OsOz0fGO+3TUQcITvDOnWIfX2ApGWSeiT17N+/f4Rfx8za3ZkTuoYVt/rUdQYTEd8CkHRi9jJ+McLPuwO4nmyAwPXAt4FLAdX62CHijHCb/sGIVcAqgEqlUnMdM+t8y+fPYuXqrf2aybqOG8fy+bNKzKr91XUGI+mfSdoMvAhsk7RJ0jnD/bCI2BcRRyPiPeB7ZH0kkJ1lTM2tOgXYk+JTasT7bSNpPHAycGCIfZmZ1bRodjc3LD6X7gldCOie0MUNi8/1KLIG1dsHswr4g4h4EkDSRWQF4p8P58MkTY6IvenlF8kKFsBa4AeSbgHOJOvMfzYijko6JOlCYCNwCfDnuW2WAk8DFwMbIiIkrQP+S675bR6wcjh5mtnYs2h2twvKKKu3wJxQLS4AEfGUpBOG2kDS/cBFwGmSdpON7LpI0nlkTVa7gK+n/W2T9BDwEnAEuCIiqueql5ONSOsi69x/PMXvBO5LAwIOkI1CIyIOSLoeeC6td121w9/MzJpHEcfuepD0CPA8cF8K/Q5QiYhFxaXWXJVKJXp6espOw8ysrUjaFBGVWu/VO4rsUmASsBp4JD3/2uikZ2ZmnajeUWRvAt8oOBczM+sgQxYYSf81Iq6W9FfUGOobEV8oLDMzM2trxzqDqfa5/GnRiZiZWWcZssBExKb09LyI+LP8e5KuAv6mqMTMzKy91dvJX2tG4t8dxTzMzBq2ZnMfc2/cwIwVjzH3xg2eDblkx+qD+Qrwb4EZktbm3joReKPIxMzMhqM65X51upfqlPuAL6AsybH6YP4W2AucRjZvWNUh4IWikjIzG66hptx3gSnHsfpgXgFeAT7dnHTMzEbGU+63nnonu7xQ0nOSfiHpV5KOSnq76OTMzOrlKfdbT72d/N8BvgLsJJsT7Pd4f9JJM7PSLZ8/i67jxvWLecr9ctU72SUR0StpXJqE8n9I+tsC8zIzG5ZqP8vN63aw5+BhzpzQxfL5s9z/UqJ6C8wvJX0Y2CLpT8g6/oecTdnMrNk85X5rqbfAfBUYB1wJ/AeyG3r9dlFJmdnYtWZzn89COkS9k12+kp4eBr5VXDpmNpb5WpbOcqwLLbcyyP3sASLiE6OekZmNWb6WpbMc6wzmN5uShZkZvpal09RzoaWZWVOcOaGLvhrFxNeytKd6L7Q8JOnttPyDL7Q0syL4WpbOUm8n/4n515IWAXOKSMjMxi5fy9JZ6r7QMi8i1khaMdrJmJn5WpbOUVeBkbQ49/JDQIUhRpeZmZnVewbzW7nnR4BdwMJRz8bMzDpGvX0wXys6ETMz6yz1jiL7qKS/krRf0uuSHpX00aKTMzOz9lXvdP0/AB4CJgNnAn8J3F9UUmZm1v7q7YNRRNyXe/0/JV1ZREJm1v48YaVB/QXmyTQs+QGy0WNfBh6TNBEgIg4UlJ+ZtRlPWGlV9RaYL6fHrw+IX0pWcNwfY2aAJ6y099U7imxG0YmYWWfwhJVWVe+FlscBlwOfSaGngO9GxLsF5WVmbcoTVlpVvaPI7gDOB25Py/kpNihJd6UhzS/mYhMlrZe0Mz2ekntvpaReSTskzc/Fz5e0Nb13mySl+PGSHkzxjZKm57ZZmj5jp6SldX5HMxsFnrDSquotMJ+KiKURsSEtXwM+dYxt7gYWDIitAJ6IiJnAE+k1ks4GlgDnpG1ul1T9Db0DWAbMTEt1n5cBb0bEWcCtwE1pXxOBa4ALyCbkvCZfyMysWItmd3PD4nPpntCFgO4JXdyw+Fz3v4xB9XbyH5X0sYj4GWQXXgJHh9ogIn6SP6tIFgIXpef3kDW1/VGKPxAR7wAvS+oF5kjaBZwUEU+nz70XWAQ8nra5Nu3rYeA76exmPrC+OrJN0nqyouTrdsyaxBNWGtRfYJaTDVX+eXo9HRjJ9DFnRMRegIjYK+n0FO8GnsmttzvF3k3PB8ar27yW9nVE0lvAqfl4jW36kbSM7OyIadOmjeDrmJnZYOptIvs/wHeB99LyXeDpUcxDNWIxRHyk2/QPRqyKiEpEVCZNmlRXomZmVp96C8y9wAzg+rTMAO4bcova9kmaDJAeX0/x3cDU3HpTgD0pPqVGvN82ksYDJwMHhtiXmZk1Ub0FZlZE/F5EPJmWZcDHR/B5a4HqqK6lwKO5+JI0MmwGWWf+s6k57ZCkC1P/yiUDtqnu62JgQ0QEsA6YJ+mU1Lk/L8XMzKyJ6u2D2Szpwoh4BkDSBWTNZoOSdD9Zh/5pknaTjey6EXhI0mXAq8CXACJim6SHgJfI7jdzRURUBxFcTjYirYusc//xFL8TuC8NCDhANgqNiDgg6XrgubTedZ7Kxmx4PJeYjQZlf/QfYyVpOzCLrCgATAO2k/XHRER8orAMm6RSqURPT0/ZaZiVbuBcYpBdx+KhxlaLpE0RUan1Xr1nMAOvZzGzDuW5xGy01DsX2StFJ2JmrcFzidloqbeT38zGiMHmDPNcYjZcLjBm1o/nErPRUm8fjJmNEdV+Fo8is0a5wJjZB3guMRsNbiIzM7NCuMCYmVkhXGDMzKwQLjBmZlYIFxgzMyuER5GZdShPWGllc4Ex60ADJ6zsO3iYlau3ArjIWNO4icysAw01YaVZs7jAmHUgT1hprcAFxqwDecJKawUuMGYdyBNWWitwJ79ZB/KEldYKXGDMOpQnrLSyuYnMzMwK4QJjZmaFcBOZWQvz1fjWzlxgzFqUr8a3ducmMrMW5avxrd25wJi1KF+Nb+3OBcasRflqfGt3LjBmLcpX41u7cye/WYvy1fjW7lxgzFqYr8a3duYmMjMzK4QLjJmZFcIFxszMClFKgZG0S9JWSVsk9aTYREnrJe1Mj6fk1l8pqVfSDknzc/Hz0356Jd0mSSl+vKQHU3yjpOlN/5JmyZrNfcy9cQMzVjzG3Bs3sGZzX9kpmTVFmWcw/zIizouISnq9AngiImYCT6TXSDobWAKcAywAbpdUHbt5B7AMmJmWBSl+GfBmRJwF3Arc1ITvY/YB1ele+g4eJnh/uhcXGRsLWqmJbCFwT3p+D7AoF38gIt6JiJeBXmCOpMnASRHxdEQEcO+Abar7ehj4bPXsxqyZPN2LjWVlFZgAfixpk6RlKXZGROwFSI+np3g38Fpu290p1p2eD4z32yYijgBvAacOTELSMkk9knr2798/Kl/MLM/TvdhYVlaBmRsRnwQ+B1wh6TNDrFvrzCOGiA+1Tf9AxKqIqEREZdKkScfK2WzYPN2LjWWlFJiI2JMeXwceAeYA+1KzF+nx9bT6bmBqbvMpwJ4Un1Ij3m8bSeOBk4EDRXwXs6F4uhcby5peYCSdIOnE6nNgHvAisBZYmlZbCjyanq8FlqSRYTPIOvOfTc1ohyRdmPpXLhmwTXVfFwMbUj+NWVMtmt3NDYvPpXtCFwK6J3Rxw+JzfXW+jQllTBVzBvBI6nMfD/wgIv5a0nPAQ5IuA14FvgQQEdskPQS8BBwBroiIaq/p5cDdQBfweFoA7gTuk9RLduaypBlfzKwWT/diY5X8h32mUqlET09P2WmYmbUVSZtyl5v048kuzY5hzeY+z2hsNgIuMGZDqF4oWb2WpXqhJOAiY3YMrXShpVnL8YWSZiPnAmM2BF8oaTZyLjBmQ/CFkmYj5wJjNgRfKGk2cu7ktzFhpCPBqut4FJnZ8LnAWMdrdCSYL5Q0Gxk3kVnH80gws3K4wFjH80gws3K4wFjH80gws3K4wFjH80gws3K4k9/ahkeCmbUXFxhrCx4JZtZ+3ERmbcEjwczajwuMtQWPBDNrPy4w1hY8Esys/bjAWFOt2dzH3Bs3MGPFY8y9cQNrNvfVtZ1Hgpm1H3fyW9M00lHvkWBm7ccFxppmqI56jwQz6zxuIrOmcUe92djiM5gxaKQXLDbqzAld9NUoJu6oN+tMPoMZY6r9IH0HDxO83w9Sb2d7I9xRbza2uMCMMWVesLhodjc3LD6X7gldCOie0MUNi891v4pZh3IT2RhTdj+IO+rNxg4XmDY10n4U94OYWbO4iawNNdKP4n4QM2sWn8GUaKRnIY1cT+ILFs2sWVxgGjTSItHIVe2N9qO4H8TMmsFNZA1opKmqkdFcnvjRzNqBC0wDGikSjZyFuB/FzNpBRxcYSQsk7ZDUK2nFaO+/kSLRyFmIrycxs3bQsX0wksYB/w3418Bu4DlJayPipdH6jEaG/C6fP6tfHwwM7yzE/Shm1uo6+QxmDtAbET+PiF8BDwALR/MDGmmq8lmImXW6jj2DAbqB13KvdwMX5FeQtAxYBjBt2rRhf0CjQ359FmJmnayTC4xqxKLfi4hVwCqASqUSNdY/JhcJM7PaOrmJbDcwNfd6CrCnpFzMzMacTi4wzwEzJc2Q9GFgCbC25JzMzMaMjm0ii4gjkq4E1gHjgLsiYlvJaZmZjRkdW2AAIuJHwI/KzsPMbCzq5CYyMzMrkSJGNHiq40jaD7xSYgqnAX9f4ucPxnkNX6vm5ryGx3nV59ciYlKtN1xgWoSknoiolJ3HQM5r+Fo1N+c1PM6rcW4iMzOzQrjAmJlZIVxgWseqshMYhPMavlbNzXkNj/NqkPtgzMysED6DMTOzQrjAmJlZIVxgSiBpqqQnJW2XtE3SVSl+raQ+SVvS8vkSctslaWv6/J4UmyhpvaSd6fGUJuc0K3dMtkh6W9LVZRwvSXdJel3Si7nYoMdH0sp0R9UdkuY3Oa+bJf2dpBckPSJpQopPl3Q4d9z+osl5DfpzK/l4PZjLaZekLSnezOM12P8Npf+OjUhEeGnyAkwGPpmenwj8X+Bs4FrgP5ac2y7gtAGxPwFWpOcrgJtKzG8c8P+AXyvjeAGfAT4JvHis45N+pj8FjgdmAD8DxjUxr3nA+PT8plxe0/PrlXC8av7cyj5eA97/NvCfSjheg/3fUPrv2EgWn8GUICL2RsTz6fkhYDvZDdJa1ULgnvT8HmBReanwWeBnEVHKrAsR8RPgwIDwYMdnIfBARLwTES8DvWR3Wm1KXhHx44g4kl4+Q3bLiqYa5HgNptTjVSVJwL8B7i/is4cyxP8Npf+OjYQLTMkkTQdmAxtT6MrUpHFXs5uikgB+LGlTuuMnwBkRsReyfwDA6SXkVbWE/v/wyz5eMPjxqXVX1bL+kLgUeDz3eoakzZL+RtJvlJBPrZ9bqxyv3wD2RcTOXKzpx2vA/w3t8Dv2AS4wJZL0EeCHwNUR8TZwB/Ax4DxgL9lperPNjYhPAp8DrpD0mRJyqEnZfX2+APxlCrXC8RrKMe+q2pQkpG8CR4Dvp9BeYFpEzAb+APiBpJOamNJgP7eWOF7AV+j/R0zTj1eN/xsGXbVGrGWuPXGBKYmk48h+gb4fEasBImJfRByNiPeA71HCqW5E7EmPrwOPpBz2SZqc8p4MvN7svJLPAc9HxL6UY+nHKxns+JR+V1VJS4HfBP5dpEb71JzyRnq+iazd/uPNymmIn1srHK/xwGLgwWqs2cer1v8NtPDv2FBcYEqQ2njvBLZHxC25+OTcal8EXhy4bcF5nSDpxOpzsk7iF8nuBLo0rbYUeLSZeeX0+8uy7OOVM9jxWQsskXS8pBnATODZZiUlaQHwR8AXIuKXufgkSePS84+mvH7exLwG+7mVerySfwX8XUTsrgaaebwG+7+BFv0dO6ayRxmMxQX4F2SnsS8AW9LyeeA+YGuKrwUmNzmvj5KNSPkpsA34ZoqfCjwB7EyPE0s4Zv8YeAM4ORdr+vEiK3B7gXfJ/nq8bKjjA3yT7C/eHcDnmpxXL1n7fPV37C/Sur+dfr4/BZ4HfqvJeQ36cyvzeKX43cDvD1i3mcdrsP8bSv8dG8niqWLMzKwQbiIzM7NCuMCYmVkhXGDMzKwQLjBmZlYIFxgzMyuEC4yZmRXCBcasg1QvCDRrBS4wZiWRdH31fh/p9X+W9A1JyyU9lyaD/Fbu/TVpEtJtuYlIkfQLSddJ2gh8uslfw2xQLjBm5bmTNP2HpA+RzRS9j2y6jzlkk0Gen5tw9NKIOB+oAN+QdGqKn0B2v5ILIuJ/NzF/syGNLzsBs7EqInZJekPSbOAMYDPwKbI54Dan1T5CVnB+QlZUvpjiU1P8DeAo2eSIZi3FBcasXP8d+F3gnwB3kd1Q7YaI+G5+JUkXkU3E+OmI+KWkp4B/lN7+h4g42qR8zermJjKzcj0CLCA7c1mXlkvT/UCQ1C3pdOBk4M1UXP4pcGFZCZvVy2cwZiWKiF9JehI4mM5Cfizp14Gns5nb+QXwO8BfA78v6QWyWXOfKStns3p5NmWzEqXO/eeBL0X/W/SatT03kZmVRNLZZPdsecLFxTqRz2DMzKwQPoMxM7NCuMCYmVkhXGDMzKwQLjBmZlYIFxgzMyvE/wdAvNFOWRTjqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(t,pop,'o')\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('population')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'ln(population)')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUgUlEQVR4nO3dfbBcdX3H8fcXiBrlIfKkEIgBq6ko1sCV+gTVUg04CBEfGduhg1OqtYPWaQYoM0WrDtKMdvpklRYKOohohcjUUWAQoXZEmxAeghhQgZaQkgimYImU4Ld/7Lmyuezee+69e87Z3fN+zezc3d/u3v1ydvnk3N/5ne9GZiJJao9dmi5AklQvg1+SWsbgl6SWMfglqWUMfklqmd2aLqCMfffdN5cuXdp0GZI0UtatW/fTzNxv6vhIBP/SpUtZu3Zt02VI0kiJiPt6jTvVI0ktY/BLUssY/JLUMga/JLWMwS9JLTMSq3okqW3WrN/E6qs38sC27Ry4aCGrVixj5fLFA/ndBr8kDZk16zdx9hW3s/2JJwHYtG07Z19xO8BAwt+pHkkaMquv3vir0J+0/YknWX31xoH8foNfkobMA9u2z2p8tgx+SRoyBy5aOKvx2TL4JakCa9Zv4rWf/BaHnPV1XvvJb7Fm/abSz121YhkLF+y609jCBbuyasWygdTmwV1JGrD5HpydfIyreiRpREx3cLZseK9cvnhgQT+VUz2SNGBVH5ydL4Nfkgas6oOz82XwS1Ifcz1AW/XB2flyjl+SepjPAdqqD87Ol8EvST3M9wBtlQdn58upHknqYdgP0M6HwS9JPQz7Adr5MPglqYdhP0A7H87xSxprc+1rP+wHaOfD4Jc0tgbROmEcgn4qp3okja2q+9qPKoNf0tga55U581FZ8EfERRGxJSI2dI2tjogfRsRtEXFlRCyq6vUlaZxX5sxHlXv8FwPHTRm7FnhZZr4cuAs4u8LXl9Ry47wyZz4qC/7MvBF4eMrYNZm5o7h5E3BQVa8vSSuXL+a8kw9n8aKFBLB40ULOO/nwsTxgOxtNruo5Dbi8350RcTpwOsCSJUvqqknSmBnXlTnz0UjwR8Q5wA7g0n6PycwLgAsAJiYmsqbSJA2hua7FV2+1B39EnAqcABybmQa6pGnNdy2+nq7W5ZwRcRxwJnBiZj5W52tLGk2uxR+8KpdzXgZ8F1gWEfdHxHuBvwP2AK6NiFsi4rNVvb6k8eBa/MGrbKonM0/pMXxhVa8naTwduGghm3qEfNvX4s+HZ+5KGmquxR88m7RJqtx8VuWMc5fMphj8kio1iFU5rsUfLKd6JFXKVTnDx+CXVClX5Qwfp3oklTLXeXpX5Qwf9/glzWhynn7Ttu0kT83Tr1m/acbnuipn+Bj8kmY0n3l6O2QOH6d6JM1ovvP0rsoZLu7xS5qR32Q1Xgx+STNynn68ONUjaUaePTteDH6pRebbOsGgHw8Gv9QSfqGJJjnHL7WErRM0yeCXWsLWCZpk8Est4ZJMTTL4pZZwSaYmeXBXagmXZGqSwS+NGJdkar4MfmmEuCRTg+AcvzRCXJKpQTD4pRHikkwNgsEvjRCXZGoQDH5phLgkU4PgwV1phLgkU4Ng8EsjxiWZmi+neiSpZdzjlxown5OwpPky+KWaeRKWmuZUj1QzT8JS0wx+qWaehKWmGfxSzTwJS00z+KWaeRKWmubBXalmnoSlphn8UgM8CUtNcqpHklrGPX5pDjwBS6PM4JdmyROwNOoqm+qJiIsiYktEbOgae0dE3BERv4yIiapeW6qSJ2Bp1FU5x38xcNyUsQ3AycCNFb6uVClPwNKoKx38EfHciHhpRBwaETM+LzNvBB6eMnZnZrpbpJHmCVgaddMGeETsFRF/FhG3AzcBnwO+DNwXEV+JiDdUVVhEnB4RayNi7datW6t6GWnWPAFLo26mg7v/AnweODozt3XfERFHAr8XEYdm5oWDLiwzLwAuAJiYmMhB/35prjwBS6Nu2uDPzDdOc986YN3AK5JGgCdgaZSVXs4ZEYuBF3Q/p5jHlySNkFLBHxHnA+8CfgBMrmNLplmdExGXAa8H9o2I+4Fz6Rzs/VtgP+DrEXFLZq6Yc/WSpFkru8e/EliWmY+X/cWZeUqfu64s+zskSYNXNvh/AiwASge/NOxsu6C2Khv8jwG3RMR1dIV/Zp5RSVVSxWy7oDYrG/xXFRdpLEzXdsHg17grFfyZeUlEPAN4cTG0MTOfqK4sqVq2XVCblWrZEBGvB+4G/h74DHBXRBxTXVlStWy7oDYr26vnU8CbMvO3MvMYYAXwV9WVJVXLtgtqs7Jz/Au6m6tl5l0RsaCimqTK2XZBbVY2+NdGxIXAF4rb78F2DRpxtl1QW5UN/vcDHwDOAILOGbufqaooSVJ1yq7qeRz4dHGRJI2waYM/Ir6cme8s+vE/rTVyZr68ssokSZWYaY//g8XPE6ouRJJUj5n68W8urv5RZp7ZfV/RsfPMpz9Lqo/9dqTZK7uOv9cXshw/yEKk2Zrst7Np23aSp/rtrFm/qenSpKE203fuvr+Y318WEbd1Xe4BbqunRKm36frtSOpvpjn+LwLfAM4DzuoafzQzH66sKqkE++1IczPtHn9m/k9m3puZp2TmfcB2Oqt7do+IJbVUKPVhvx1pbso2aXtLRNwN3APcANxL5y8BqTH225HmpuzB3Y8DrwLuysxDgGOBf6+sKqmElcsXc97Jh7N40UICWLxoIeedfLireqQZlG3Z8ERmPhQRu0TELpl5fbGcU2qU/Xak2Ssb/NsiYnc6PXoujYgtwI7qypIkVaXsVM9JdA7s/gnwTeDHwFuqKkqSVJ2yTdr+t+vmJRXVopby7FupXjM1aXuUHs3Z6LRmzszcs5Kq1BqTZ99Onog1efYtYPhLFZlpHf8emblnj8sehr4GwbNvpfqVmurpd7JWZv7nYMtR23j2rVS/sqt6vt51/VnAIcBG4KUDr0itcuCihWzqEfKefStVp9Sqnsw8vOvyIuAo4DvVlqY28OxbqX5l9/h3kpk3R8QrB12M2mfyAK6reqT6lJ3j/3DXzV2AI4CtlVSk1vHsW6leZff49+i6voPOnP9XB1+OJKlqZU/g+ihAROzZuZmPVlqVJKkyZdsyTxTfxHUbcHtE3BoRR1ZbmiSpCmWnei6i84Xr/wYQEa8D/hl4eVWFSZKqUbZJ26OToQ+Qmd8BnO6RpBFUdo//+xHxOeAyOr173gV8OyKOgM7yzorqkyQNWNngf0Xx89wp46+h8w/Bbw+qII0eu2tKo6Xsqp43VF2IRpPdNaXRU3ZVz14R8emIWFtcPhURe1VdnIaf3TWl0VP24O5FdA7mvrO4PEJnVU9fEXFRRGyJiA1dY3tHxLURcXfx87lzLVzDwe6a0ugpG/wvzMxzM/MnxeWjwKEzPOdi4LgpY2cB1xWN3q4rbmuE9euiaXdNaXiVDf7txdp9ACLitXS+g7evzLwReHjK8Ek89dWNlwArS76+hpTdNaXRU3ZVz/uBS4p5/aAT6KfO4fWel5mbATJzc0Ts3++BEXE6cDrAkiU9vwdGQ8DumtLoicxeX6nb58GdXj1k5iMlH78U+NfMfFlxe1tmLuq6/2eZOeM8/8TERK5du7Z0nZIkiIh1mTkxdbzsqp59IuJvgG8D10fEX0fEPnOo48GIOKD4nQcAW+bwOyRJ81B2jv9LdPrvvw14e3H98jm83lU8NUV0KvC1OfwOSdI8lA3+vTPzY5l5T3H5OLBouidExGXAd4FlEXF/RLwX+CTwxoi4G3hjcVuSVKOyB3evj4h3A18ubr+dnb+A/Wky85Q+dx1b8jUlSRUou8f/h8AXgceLy5eAD0fEoxFR6kCvJGk4lO3Vs8fMj5IkjYJp9/iL5ZjT3R8RcdBAK5IkVWqmPf7VEbELndU36+is5nkW8GvAG+jM158L3F9lkaqerZWl9pg2+DPzHRFxGPAe4DTgADqtGu6kc3D3E5n5i8qrVKVsrSy1y4xz/Jn5A+CcGmpRQ6ZrrWzwS+On7HJOIuI1wNLu52Tm5yuoSTWztbLULqWCPyK+ALwQuAWY3DVMwOAfAwcuWsimHiFva2VpPJXd458ADsvZdHTTyFi1YtlOc/xga2VpnJUN/g3A84HNFdaihthaWWqXssG/L/CDiPg+nTN3AcjMEyupSrVbuXyxQS+1RNng/0iVRUiS6lO2ZcMNVRciSarHtMEfEY/SWb3ztLuAzMw9K6lKklSZmc7ctTmbJI2Zsm2ZJUljovSZu6qHzdIkVc3gHyI2S5NUB6d6hsh0zdIkaVAM/iFiszRJdTD4h0i/pmg2S5M0SAb/EFm1YhkLF+y605jN0iQNmgd3h4jN0iTVweAfMjZLk1Q1p3okqWUMfklqGYNfklrG4JekljH4JallDH5JahmDX5JaxuCXpJYx+CWpZQx+SWoZg1+SWsbgl6SWMfglqWUMfklqmUaCPyI+GBEbIuKOiPhQEzVIUlvVHvwR8TLgD4CjgN8AToiIF9VdhyS1VRN7/C8BbsrMxzJzB3AD8NYG6pCkVmoi+DcAx0TEPhHxbODNwMFTHxQRp0fE2ohYu3Xr1tqLlKRxVXvwZ+adwPnAtcA3gVuBHT0ed0FmTmTmxH777VdzlZI0vho5uJuZF2bmEZl5DPAwcHcTdUhSGzXyZesRsX9mbomIJcDJwKubqKMqa9ZvYvXVG3lg23YOXLSQVSuW+QXqkoZGI8EPfDUi9gGeAD6QmT9rqI6BW7N+E2dfcTvbn3gSgE3btnP2FbcDGP6ShkIjwZ+ZRzfxunVYffXGX4X+pO1PPMnqqzca/JKGgmfuDtgD27bPalyS6mbwD9iBixbOalyS6mbwD9iqFctYuGDXncYWLtiVVSuWNVSRJO2sqYO7Y2tyHt9VPZKGlcFfgZXLFxv0koaWUz2S1DIGvyS1jMEvSS1j8EtSyxj8ktQyBr8ktYzBL0kt4zr+HmyrLGmcGfxT2FZZ0rhzqmeK6doqS9I4MPinsK2ypHFn8E9hW2VJ487gn8K2ypLGnQd3p7CtsqRxZ/D3YFtlSePMqR5JahmDX5JaxuCXpJYx+CWpZQx+SWoZg1+SWmZsl3PaYVOSehvL4LfDpiT1N5ZTPXbYlKT+xjL47bApSf2NZfDbYVOS+hvL4LfDpiT1N5YHd+2wKUn9jWXwgx02JamfsZzqkST1Z/BLUssY/JLUMga/JLWMwS9JLROZ2XQNM4qIrcB9DZawL/DTBl+/H+uaHeuaHeuanWGs6wWZud/UwZEI/qZFxNrMnGi6jqmsa3asa3asa3aGta5enOqRpJYx+CWpZQz+ci5ouoA+rGt2rGt2rGt2hrWup3GOX5Jaxj1+SWoZg1+SWsbg7xIRB0fE9RFxZ0TcEREfLMY/EhGbIuKW4vLmBmq7NyJuL15/bTG2d0RcGxF3Fz+fW3NNy7q2yS0R8UhEfKip7RURF0XElojY0DXWdxtFxNkR8aOI2BgRK2qua3VE/DAibouIKyNiUTG+NCK2d227z9ZcV9/3ruHtdXlXTfdGxC3FeC3ba5psaPzzNSeZ6aW4AAcARxTX9wDuAg4DPgL8acO13QvsO2XsL4GziutnAec3WN+uwH8DL2hqewHHAEcAG2baRsX7eivwTOAQ4MfArjXW9SZgt+L6+V11Le1+XAPbq+d71/T2mnL/p4A/r3N7TZMNjX++5nJxj79LZm7OzJuL648CdwLD3NT/JOCS4volwMrmSuFY4MeZ2dgZ1pl5I/DwlOF+2+gk4EuZ+Xhm3gP8CDiqrroy85rM3FHcvAk4qIrXnm1d02h0e02KiADeCVxWxWtPU1O/bGj88zUXBn8fEbEUWA58rxj64+LP8ovqnlIpJHBNRKyLiNOLsedl5mbofDCB/Ruoa9K72fl/xqa316R+22gx8F9dj7uf5v6RPw34RtftQyJifUTcEBFHN1BPr/duWLbX0cCDmXl311it22tKNozC5+tpDP4eImJ34KvAhzLzEeAfgBcCrwA20/lTs26vzcwjgOOBD0TEMQ3U0FNEPAM4EfhKMTQM22sm0WOs9rXNEXEOsAO4tBjaDCzJzOXAh4EvRsSeNZbU770biu0FnMLOOxi1bq8e2dD3oT3GhmbtvME/RUQsoPPGXpqZVwBk5oOZ+WRm/hL4Rxr4ky0zHyh+bgGuLGp4MCIOKOo+ANhSd12F44GbM/PBosbGt1eXftvofuDgrscdBDxQZ2ERcSpwAvCeLCaGi6mBh4rr6+jMDb+4rpqmee+GYXvtBpwMXD45Vuf26pUNDPHnazoGf5di/vBC4M7M/HTX+AFdD3srsGHqcyuu6zkRscfkdToHBjcAVwGnFg87FfhanXV12WkvrOntNUW/bXQV8O6IeGZEHAK8CPh+XUVFxHHAmcCJmflY1/h+EbFrcf3Qoq6f1FhXv/eu0e1V+B3gh5l5/+RAXdurXzYwpJ+vGTV9dHmYLsDr6Pw5dhtwS3F5M/AF4PZi/CrggJrrOpTOCoFbgTuAc4rxfYDrgLuLn3s3sM2eDTwE7NU11sj2ovOPz2bgCTp7XO+dbhsB59DZQ9wIHF9zXT+iMwc8+Tn7bPHYtxXv8a3AzcBbaq6r73vX5PYqxi8G3jflsbVsr2myofHP11wutmyQpJZxqkeSWsbgl6SWMfglqWUMfklqGYNfklrG4JekljH4pZpMnmgkNc3gl3qIiI9N9lwvbn8iIs6IiFUR8R9FE7OPdt2/pmigd0dXEz0i4ucR8RcR8T3g1TX/Z0g9GfxSbxdSnIofEbvQ6T76IJ1T74+i08TsyK5meadl5pHABHBGROxTjD+HTr/438zM79RYv9TXbk0XIA2jzLw3Ih6KiOXA84D1wCvp9ElaXzxsdzr/ENxIJ+zfWowfXIw/BDxJp7GXNDQMfqm/fwJ+H3g+cBGdL5s5LzM/1/2giHg9nQZir87MxyLi28Czirt/kZlP1lSvVIpTPVJ/VwLH0dnTv7q4nFb0ZCciFkfE/sBewM+K0P914FVNFSyV4R6/1Edm/l9EXA9sK/bar4mIlwDf7XTp5efA7wLfBN4XEbfR6cR4U1M1S2XYnVPqozioezPwjtz5q/6kkeZUj9RDRBxGp2f+dYa+xo17/JLUMu7xS1LLGPyS1DIGvyS1jMEvSS1j8EtSy/w/KE+ssYn+AfIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(t,logpop,'o')\n",
    "plt.xlabel('year')\n",
    "plt.ylabel('ln(population)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercises:\n",
    "\n",
    "* Estimate the parameters $a_1$ and $a_2$ is to performing a linear least-squares fit to the log-transformed data shown above. Note that log-transforming both sides of the first equation gives $\\ln(y) = \\ln(a_1) + a_2 t$\n",
    "* Plot the original data (not transformed) and the exponential fit using these parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the cost function\n",
    "\n",
    "An extremely inefficient way of optimizing the parameters would be to try every possible combination of parameters within some range. This is called the \"brute force\" method because it relies on computer power rather than intelligence.\n",
    "\n",
    "Even though it is inefficient, it is a good way to visualize the cost function in two dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-7-1a1c6e573cc9>:18: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.\n",
      "  plt.pcolor(alla2,alla1,np.log(rms),cmap='viridis')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$a_1$')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAACI/ElEQVR4nO29e/R9SVUf+Nnnfr9NCIq8TdPdCUbRqBhROoSZZI0aYiAkEcxAbCYjuEJWq8FRV8wKoJnRcUJGJ76CRhQj8hjkMaiBGBlFjTJOoLF1kLfaCkpLh04DYksiv+/3nD1/VO2qXbuqzuPec+/3VXut8/udetc55373e+8iZkaDBg0aNLh60J31Bho0aNCgwdlAIwANGjRocEWhEYAGDRo0uKLQCECDBg0aXFFoBKBBgwYNrig0AtCgQYMGVxQaAWjQoEGDKwqNADRo0KDBFYVGAC4BENH7ieivrzTX/05E37jGXBcNiOitRPS5Z72PXeCsv99leIdXCRoBaBCAiB4K4BkAftiX309E/5WI/piI/hMRvYSIPsmMeT8RXSOih5j6txERE9EjfPmvEtF/JKKPEdFHiOj/JaK/dKBHmwvfBeDbz3oT24L9fivN+VVE9MUj7X9ARI9WVRf6HV41aASggYavAvAzzPxfVd3fYeZPAvBoAF8A4HmFce8D8HQpENHnAbivKt8fwE8D+H4ADwJwA4D/FcAndt0wER3NqZs5x+sBfAkRXb/rvs4Ivgr599sKiOiriejLY5FuVWWpfAiAhwF4j6q+6O/wSkEjAJcMiOizieiXiOgPiehdRPRlqu0Liej/I6J7iej/IqJXE9E/V8P/JoBfLs3LzP8JwM/CEQILL4fjPAWeCeBlqvyZfo5XMnPPzP+VmX+Omd8+8hwPJ6KfIKL/TETvI6KvV23vJ6LnENHbAXyciI4qdWPvIuvPzH8C4NcA/I3avgr7fBgRvZ6IPkREf0RE/84TvLExxWcjok/30tEXqn73CAfu9/w8Ino3EX2UiH6MiP6Umjr5ftvsTcGLAXw6gG8E8C8ADABep+b+DAAfgMMhHyaiD2/7DhucHTQCcImAiI4B/DsAPwfHmf1PAF5BRJ9FRNcB+CkAL4Hjwl8J4MvNFJ8H4Dcrc98Ih2DuKDS/BcD9PcLdAPgKAP+nav8tAD0RvZSI/iYRPXDiOTr/HL8BJy08HsA3EtETVLenA/hbAB7AzKe2DgDV3sXEHO8B8PlqLz9IRD84st37w0k2fxbAIwA8BMBXb/NszPw7AJ7j9/mnAfwYgJcw8y+pKf4+gCfAIefPBPDPVJv9fov2VgBW//eqDGa+A8A/AfBaZv4kZn5w7R02OL/QCMDlgscB+CQA38HM15j5F+FUL0/3bUcAXsDMJ8z8kwDeasY/AMC9pu7fEtG9cNze3QC+tbK2SAFfCuC9AP5AGpj5jwD8VTgE8iMA/rPnTD+1MtdfAvBQZv52/xy/68fdovq8gJk/YNQdum7sXYzNca9/D7L3f8TM/6iyTzDzHcz8Rmb+BDN/BMAbATyQiD6FnEH0j4noUXOfjZl/BMBvA7gNwPUAvsUs+QN+zx8B8HzzPA+A+n61vQEAEf03RPRmIvplInqlZx40/AM41d73+T1cB+DJps/nA3hb4bUk77DB+YVGAC4XPBzAB5h5UHW/B8dpPhzAH3Ca//sDZvxHAXyyqXsKM38ygC8G8BfguMgSvBzA/wCnh36ZbWTm9zDzVzHzjQAe5ffzfZW5/hyAh3vVzR8S0R8C+GYAmmDYvdu6sXcxNscnA/jDyr4yIKKneYP23X6fz4WTeP4LnHTxWjNkzrP9CNw7+n5mtnYSveffg3tOgeT7jexNxv41Zv4iAL8Lg9yZ+Yc9k+CL/MPM/G/NXh4NJ8lYWPQOG5wdNAJwueCDAG7yagaBPwvHjd8F4AYiItV2kxn/dnh9vQVm/mU49dF3Vdp/D45jfBKAnyz1UX3f6+d6VKXLBwC8j5kfoK5PZuYn6WlKU6v7sXcxNsdno4zUMiCivwbgO+H05A+HI453A3ibl7L+c2HY6LOR87L6PgA/CuDbiOhBZrz+Zn/WP6dA+H5jewMAZv6gknxO4XT8GTCzVUHJs3dw3+9thWGz32GDs4VGAC4X3Abg4wD+KREde+Ph3wHwKgBvhtPjfp03kD4ZwGPN+J8B8EUj838fgC+l1O1Pw7PguMqP60oi+gtE9E3ejgAiuglOdfGWyjxvBfBH3kh7XyLaENGjaJnb6Ni7KAIR3QfAY+BUJXPg8+EQ+nvhVCsvhrM3vHtkzNSz/SsAv8bM/xDAvwfwQ2b8s4noRk8YvhnAq1Wb/n6z9kZEnwZn2/npmc8scF9/JThki3fY4AyhEYBLBMx8DcCXwf1B3wPgBwE8g5nf69v+LhyS/kMA/yPcH71WMbwMwJOI6L4ogOdoXwbgf660/w4z315ouhfAXwZwGxF9HA7xvxPAN1Xm6eGQ9aPhpIp7APwbAJ9SefTSHNV3MTLsywD8EjMHrpqIfoiILBIWeAWAYwD/Ce5d/jaAd/u1a/uqPpsnyk8E8DW++z8G8IVE9PfVFD8OZ9j+XX9pLy79/Sb35j2CXgrgK8f2XHmOj8MRp3cT0Z2qKXuHDc4vUDsS8uoCEd0G4IeY+cdU3b8AcDczf9+ZbeyMwL+PZzHzO1ec8yUAvmuNOYno/QD+ITP//EifWd+PXOzD6wB8tzeQrwL7eIcN9geNAFwhIKIvgnMTvAfOnfCHAPx5Zr7rTDd2SYGIfgaO0/89AD/MzC/Zcb73Y4IALJjrKwF8L5wkBgAvZOZXjwxpcAlhUcTkNuADVd4E4D5+vdcy87d6Hear4fyT3w/g7zHzR/2Y58GpKnoAX8/MP+vrHwNnPLwvnL7zG5iZvd7xZXC6xw8D+Apmfv++n+0CwmcBeA2ce+TvAHhqQ/77A2O0PlfAzC+H89xqcIVh7xKA9zq5HzP/sfc1/hUA3wCnj/4IM38HET0XwAOZ+TlE9DlwQUqPhfNe+HkAn8nMPRG91Y99CxwBeAEzv4GI/hGAv8jMX0NEtwD4cmb+ir0+WIMGDRpccNi7EZgd/LEvHvuL4fyOX+rrXwrgKf7+yQBe5YNX3gcXefpYcrlF7s/Mb/a+7C8zY2Su1wJ4vHF3bNCgQYMGBg7iBeRd3d4G54f8Rma+DcCnivrB//8w3/0GpMEud/q6G/y9rU/GsAtH/xiAB+/lYRo0aNDgksDebQBAcH17NBE9AMBPmdB4CyXOnUfqx8akExPdCuBWAKDj6x5z3UMfVg4FWgCTYsYu8y8Zq/puvSeujLf9bb9Ce2zjvA+bNUr7IdNMaRuLgEdpf7b9Sn3sGnotM962V/sUy/ZBkT2H7UuVvuR3TWYMqb4k86h+pF4umfmSrZDqVxzH6nWU14j/p3PKvrLxxbU5Xz/bj957+hydmTt7VgIIlH0GXaPv5P7X3v6Je5j5odgSnvAl9+MPf6Sf1ffX3v6Jn2XmJ2671hI4CAEQYOY/JKJfgvN1/hARXc/Md3n1zt2+251Iox1vhIt2vNPf23o95k7v3vYpAD5SWP9FAF4EAH/qhpv4z331PwYA0ARCHKvLxi5p50K/me1k+1kCwPU+xKadY10oDwAp5O3KfvwgdRzKNEDdsyr7+56TssxNPZtn5gRTcUfgDYE7VT5ydQAwbADeEAb/S+YNYdjIPTBsCKzKvAEGkXs7N17kYN4A3CH27zyy9+2uL6u9uDbuIkXkjqNc3UUCwB2DOl8nH4EA8mXqGEQAde7FdX5s51961w0gAja+nYix6Tggzw0xOuKk/YiG0H7UDeh8HeDmPerifQfGUdfHNhrQ+b4bYhxRH/Zy7O83Mjf12GBI2jdq/DH12IBNecBG9gLGMZ2GtTYYknKHAdf5cgfGhgZsIGud+rn9e8CAY9Xu2oBjKROwAeHY/8Y6EDYgdP5DHdMm3G+oQ+c/5ub63/497AD3fKTHbT9743RHAMfX/04t3crqsHcVEBE91HP+8AEqfx0uOvH1cGmD4f+XVLOvB3ALEd3HRyk+EsBbvZroXiJ6nNfvP8OMkbmeCuAXed/W7TOGhAsFimw/U6HfXBiZf87ai+YGHOIPBEdxe8z+UgSIY/+EcM0BFo5PzaXawv+Fed1aFNqKxF21F19+mNeKFaZoxkqZmcL9sODjSt8BhAFxfG0OWz9wV7wv9kV9bG3fAzoMCh316NCPlEM953Xl+dNcFz0z+gKKGMoZMVYARs/DrOuQcAgJ4HoALyWXJrgD8Bpm/mkiejOA1xDRswD8PoCnAQAzv4uIXgMXsn4K4NlehQQAX4voBvoGfwEub8rLiegOOM5fZ408LMhvbyH5YfIIRY+n5fMk+9iRBIqahWq0VDh1084JhnVXYOqJsv7lxQXbk/vLFXUDAzwwqKO0m6Id2b2lK7X/7b2sN7Zdvd4YPmYCMyNzTZC5w6v079xPqptZ/z7sNL5NkN+GGAMocMJjMDChIw6IuqPe13vul3LVxcAUJADAIegOsV/PXeD4XTsl6LtHh80MZNszoVPP7J7J9AEF1Y+uA5A9f89OCojzeQlLvdgh/GAGdLQOj8wATjFPBXRI2DsBYHfoxxcU6j8Mlwu9NOb5cKlubf3tKCQQY3cIxdN23qyFGUiUFWJa3L4ASQcCUYKFyJ4pHSL3U/Mnz2LWDIxYiTMe3QtlBIbUWpn1p4AFNZIu3hfawtw1CcDuQ/UPCL/6UIX5JxC/PBrZ9aoIn0DEUVIgTpAy+x+MbRdVib63UGvTxAVwyBlAQNCxXcopAYnlwY9PCY7MuTHPPBhi0XOXESW7duhbJA6cIHy3hqiMbP06HDmjLHGcNRzUBtCgQYMGVxWGJVzagaAlg8OOuvJDg1avVKD6M5saSzP6IG9f893lenWxAURbgNbRiz0gGKK1vaCgoyd1hfWMlKDbZVzVISDo/SmMr4Lpmz6m1fmX27XefhuGUo9P7jPON5Z7s7c+65uiEauXt/2tLt+uXdP163l7UDZvbexQeE/WJhDX2Q+SZjjJY851SGgEYCnMQZD7HL/LugvA4qiiPbNGOIkUMfH3HQEdhTFMFF05S1NoJMyIf7EWiRuiYBEwMZdVPZwagjOkD2RIP665gGPwfcuEgRKjru4fiqrdEgkNgkRrRuIxg3GJINTmj2OmDbtjYPtbBD9wl60BlIhHqU+JMFQIRsUQvA9j8ACedR0SGgE4rzDj7ynDQxXispMnUGX+oifQHohb5NRZuaWyY+uGiPwzgqGkgJRgmAVKnL7uVyMOqs2uF6lcYYixS1T3ghTh54LRuERQ8/oJ7SMfShCvJQgWKe9CMKprW2mixNFbiaHiCVTzBsrrKraQFZExI3oeTV2HhGYD2BcsNMwCFUOvMoaOGoJre8DyfeT7yg216Tq0nT6iuqCaS+bWHklMiUePM+4KcYjcNoc22zeWs+9Uu5fyto8phMF+QLsXKXojb+4lFMuhz5YgCLIzBuK8X2p0tR5AAoKsxRsoGnm90dc/RGf6W28gawi2nkB6btlXnDvdV1/whLKeQK4f+/GWuKxnBD5pNoALCAsk/a3a1+Karcpmontw1aR4P1ezEfqMSBy1a87e08mM3l8CzIwaKEzl7QDVGIExiSBIG8ZGkO0pzpNKDZSOt9JEMkcuIbAnbIVwiNBeg7JEQKPtFmrSApDbAcbaxmIISuUSB2/9/m2cQOw3vnapj6sr69v3poJhR3jmXIeERgAULFKVnIUef9f156ppinr98SEByU/9orRdgDD/FxiQuueKBVNKNHJB5RMRuhpX2k4F4ZcIxA4Md2ZwTh/PUnAabddqHyEUWjVSQvg1m8AcgmDL+zAEl9RSc4jFLLVR/KmYsTnsQxUjvMuc65DQCMBZwYoEpKiPn+ozBQWuPSmb+r1CFidgyyhz5YW2ImeuxwI5kjbl1BicE5/UO0jdZ88FZBTZrGWNxHONwnOgah+Y4Ohj/bhNICMgM2wINWlgrCzj5nkJFeoqyH5dQzCFPU1dh4RGAM4Z1JDsmUocFQJT9wKC5/JJ3ft6q8ytLTkgzTGkkeuQI/OiOqiC+DNCIJCpdNL2SbXOHPAvrYhzJlQ+yVZN36JEUDD6TqV/iNJCXW0z5QGUIezF0sA0WqpJDHNUQjVDcEn9s5ZnDiNKIVPXIaEZgVeEzMZH2B5RjKwBVFQR8lu3bSvsY5YhuLoxOxfCXsOfZyWtBOCIQEi4xmLIjfp2MIIRmJgC0eCOQQrZEjN07sg4Vm8OSTurR8pwSdiLGqt+BHmbfuB8vdQIjFlRwTp1BKu1xVCrI4Ft+1RUsDYSl9qSVBDGUDtwHtGbju9CRLAbn0YFi11B5rBRx3pcKSo4MSJXDMF6/tg3jxJeAxjAtXPIb5+/HZ0x1DnaPc2t11gKE/uaYlQDB6/+D/cTkPQZ28Ok7WA8HsBC5OpFr8/BHbQcJGYXRPZiwqMHQjFiBOa0fWsoqIVEvZOqfFICVfPzXyMoLNRn+n2v92dKDL5WZVEz1Lr557qBLpcOagFhaTK5mkpof1x/Nq9/31PXIaERgF1g1291VmqdOQSNCv1mjBNDsFYPzfIEqk7oMKBWBwEIFrMS0i8ZhKsxAYb7zyKAgSLRAFI7QNVAHB7c3deJR4m9zwlE6V6DVo3oPrOCvSoEobRGOq6MRqb08nMzfM41BM/zEsr3WdP0r2kIZuBc2gCaCqhBgwYN9gwMKhK7s4bzt6MGW8FcT6BtJMyqh6KSCBbNHQzD22zGSwTKtTNL5aCkBldGqhKqlHXfWnnMlVP6V6OOC48yaohGrhbKuk+ohGr5fuZEBWuvnSnf/rStrEKqjT30+QA1T6B9p4VYSwVERC8moruJ6J2Ftn9CRExEsw6VudoEYA/SVoYIt7QflJDsTvtdW12lDZZjdpPEE6jQce6+tCooQeLRDoBBiIPuw0YFxKMqI7dW2h62WrEnzEL4ltiUXpoZn+GjkfadXUIL6qK8T92Ncw1PoDJCp9HyHE8gqUvtFe6ygVf7CgRbWQX0ErhTFRMgopsAfCnc+Sqz4GoTgKWwkiG4+re6YP59GJSnuPhZBltB+GN7s64cHcYN0BkiVBUZQUAZMVd17+nWSgg/09tPSAA2O+joutpGkDyi1f+X2y3XWJcIqNg+J5fP3KhgaySeziWUJ3ybEwcwFXRWWjvOVarjrNyvbAhmEE74aNY1ORfzm1A48hbA9wL4p5iUPSM0G0AB5Hczy9NDfk/7YRzmAY2vz6jgZFLtnNYnY+zcfr3kPS18D2HsvO5qXr9YcPt0SC24aiqEyz5vUCj7f3TAV5aYTY+3zyN99XOzQtZTRkPF+RMz8rNtlMuq9T31ZXtiWDK9mk9OA5tzIEw8Ecz11fWu7ObTeXZc2YE97Wsong6myshPB9P9S6d+WYi5gMp5hvTcG1g30dz9U5Q82WljKxqCFxh4H0JEt6vyi/x55lUgoi8D8AfM/Bu0QLXaCMBZw4oEJCNcY3NPEI2t+8o+1AlhVo1Fti8RqLZAoZrUEZHMCpkC3u8/voA0hoAy3/zMV18/a+V/TWyyvan5JI4gPGgNacuc9u/W1Nnkb7os9/bUrjFIksEpBD8dF1CaK0Xy1hffIvVSgjm3b1e2ieTsaWICOfHoAE6PsRw7MjKrZ0a3jW1qApipKKlU4B5mvnluZyL60wC+BcDfWLqvpgLagx1grXUyW4L8v+OeVzUEA9P2gNJc5jmmzwco5AEK+n+lalExAXVXUD1XxQ4A5Kok1O0Ac1RMtn/VDmD3gYJKKNTXBY+p9NBjYFU3aVsX/h9L7ZAZbmeka5hU98xID12ae2l66L3EAfh3P3VtAZ8O4NMA/AYRvR/AjQB+nYj+zNTAJgEshZkce4Xhm25fIBGMqqpKXPuu0oYdr8rVvXgOv/oyyGw02BBinY5AJmYXlRyifUUCUPqsEQ6fRhAmOF2akXLxCUefbjGtC+vF5w7r1v6+2au1DMcfVF3SrN79lEQgah8AMSrYtAMjKiDFJY+ni14WFZxz/3lU8Ng5wT0TQGkqaZseurb23PTQru96UcGMMvFaZW7mdwB4mJQ9EbiZme+ZGtskgBFYxCkfSpLYcv0pnK+5clb/F5nUuYbgUreSkXjpu1NSgOXqNTdP0k+WMZz1WECYtRWUPIT0HCEobPYz0KhL6aGyg05lCk3WUBzqPo6JTNM/p9zwXLfPOcdELk0PvY40QEECmromZyJ6JYA3A/gsIrqTiJ617a6aBHBeYKGefbU1BcbWrkkTwilrBn/ucwiLKKGZ8ruf4XItUoBb3LPlgyor//9SbqBEh68QfTAIa0mnhPzDPgrShJkvNShHjr5I9BiIrrKc/BdtHu4m4fR9eecDYioc/1Dg6MUw7PpLfZrPx+r1Z+0h4/a7YAeQNebYAbJDZozOP0QEkzVu54bgNYABnPA6MzPz0yfaHzF3riYBXCSY+YeUMXBFLn752lYiSMqZTh9ZWgiemqO014APOWVvBXkH/X9+aPxU4FctJiBD+rIVw/XHhzBrzAWRGriQHbQgGewjO6itt3PUsoOm/eZlB51KD237x73OswOkB9mXuek5dgAJClvbDVSkmKnrkNAkAGA73fgMTjfT85sxU3YC6QNUdOslDnwO7CptTIxf9Fyh7DlZcOJFZFS6vjMndAHitwqkSHhwqmLrtZNx6FYiQJyLdX8ssAPI/+oDFj2OknG5hKBdTEvZQackAq3nL2YH9fdLs4NOHRNZOsYxtqVcvXULlT5urd2yg9Y8gJZkB10L5ibEOyQ0AtCgQYMGewbG/ozAu8D529E5gzMzBG8zlzW8WrUMZjL+FeOvhsQQbFQ6xT3pA2Jse23eOcBI1EOkVUUDp221qOHinKqPrjftwO6G4No+SobgRAtWMQS7Nl+3wLVw7fTQtb25cu5aWuKQ10gPXYsInpseeg1gUHhfU9choUkA5wl2Vc3suHawi3GsA5S2YkLtk45VRtlsLaMjKrmC1hYTzYEP7iLAGYDl76ZnUBf131pfz55YkNKPhwdEwahbsgcYZF0zBNu2oDpSKp2qIRjI34FRGdVcQMcjhVOV0DKVT12lo6OC3TjnxinIWQy0mVrHGJFdXWrQlXFjhuC+YBiWtYObaChvHxW8CzBjVpqHQ8PeJQAiuomI/gMRvYeI3kVE3+Drv42I/oCI3uavJ6kxzyOiO4joN4noCar+MUT0Dt/2AvIxz0R0HyJ6ta+/jYgese/nOjOwXH4F5hqCeeZ8xbVrY8fWsiBYJZubgDFpQJ+gnQSGIePww/GSmtsOxCEfkwV+GW4/bLGE+M34SQjzjxiCNR1gc2CMaq4ZgnWbzQ5aPBO4ysXPzw465lZ6ltlBexTOIa4877oHw88LAtsyEGxrOIQK6BTANzHzZwN4HIBnE9Hn+LbvZeZH++tnAMC33QLgc+Ey3v0gEQkxfiGAWwE80l+SEe9ZAD7KzJ8BlxDpO7fa6dJ3PxMRZ14uW3zjIvJUc23j1bMtFDQT1XnFEygpU/pekvH2F1mIEE6OppSU0HJgTJIqWryCUgKRIGiL9DUxCOvlRKHkBTQWJZwlh6uofcaSw41lB82Tx1H4fyw5nL6fGxegPYOmDnkXmJMddN5hMTRa3jU7qIU1zwReKw5gTdj7asx8FzP/ur+/F8B7ANwwMuTJAF7FzJ9g5vcBuAPAY4noegD3Z+Y3s0v28jIAT1FjXurvXwvg8SIdrPIMNe61BAsR/Ojcaz3BPpiKkoShEPrcd8ZmTHHu0vKC5AdA6/2JEf6a3b0JDvNEQhODMVfR1NeyciX70i9g+jmk75gdIEHuvr+2A9iTv7QL6DZQDBSr2AGAelDY0uygUpfMtWV20LmpILKANK4Hhe0K59EN9KCredXMFwC4zVd9HRG93R9w8EBfdwOAD6hhd/q6G/y9rU/GMPMpgI8BePA+nmHvsKWEkM0xAtWftl+b1f/6PptnSXrocJH6X6l/1Jzu/wVEN1HjMGhwVw1hW44/yQ2kt15C9JYgwHD19vGTdfR8EwSCgVycso9d5vpLbfqYyCmO38Ks4yRH1DrbRAWPzS3z6zWqfv9zYggq32GdY2AcMOYdBjPnQJg14WAEgIg+CcBPAPhGZv4jOHXOpwN4NIC7AHy3dC0M55H6sTF2D7cS0e1EdHv/8Y8ve4ALCIvUT3N+dxqRV8YmXPzcOYGEIGiiUpWItdfPEJG45vBJiMGYmqeA4K2KxyL4zEaQ7Gtk7uJz6Cv18snak8c3yD99LQFqEkFJ5z/G8WvQdoApjr60Xmn+ycRyFQSfp3nYNl1ExSuIyyeFbQNXVgIgomM45P8KZv5JAGDmDzFzz8wDgB8B8Fjf/U4AN6nhNwL4oK+/sVCfjCGiIwCfgsKBCcz8Ima+mZlv3tzvfpXNLn24eWPmGGUXzVG4X8Q11+bdcsyYCkfsAPaw+ADavaRAEDKJIZk8w4wRQSs1kLiEOukAiUvoKIFQU5eMyFYiqBp/hetXc40ThlxC4DECgXkSQUoQ1FGTBZ3/mIQwni20LDFMZQd1faakgWm7wJLsoHkOo/2AOxBmM+s6JBzCC4gA/CiA9zDz96j661W3Lwcg51u+HsAt3rPn0+CMvW9l5rsA3EtEj/NzPgPA69SYZ/r7pwL4ReaVyLaCRYhyIVIdtQOsIRWW5th17glpAKg/l7UTJIRB/ypnIP3E99+rgQTZj+v3uWrMrdoCCiogawROYgJMW/4sci1PDpfawssIOiMYM6FkBB7zUhlLDpch6Ky8zA7QM2XEJPMSWpAKomgH4HVjAhwf0s26DgmHcEz9KwC+EsA7iOhtvu6bATydiB4N927eD+CrAYCZ30VErwHwbjgPomczszjvfi3ceZj3BfAGfwGOwLyciO6A4/xv2esTHQII4xzijuNrOjUZq9vlPujbLGdKlHrmICJzG1MQuHm1OSa1ltzPJUzBv9ENIOJw75K/Rb9/mxoi+OmHsvej1+odjUBLCLqmAiruNa7N4WHLnUOaihWTw8kpYaGfiQOw6aA1lJLD6eRxMr8ry5jo9z/nlDA9VpC7jJG5LX88FRPg5k2Tw9VOCQPqMQFrwIITwQ4GeycAzPwrKP85/8zImOcDeH6h/nYAjyrU/wmAp+2wzRTM39y5hZn7zPBMjTjMfe5SP1WXrCcI37KhWtVjdC1sTwgrBYaVpACp9vfEBB7gJAEAPMBw8yoSFwrpqnKyNS9QaBfSRAVjf+Vm7uJ7zeaylNO06/dkEHxCtE3fmAvIV1DkhjfECXEIBMEgeA36BLA8W2iK3PUe3Hq+XEDq5RxDsd0GhQE5Qq8Fk7m1bRBa+myl3EFrADMdnLufA+cvNK1BgwYNLiEc2sd/Dpy/HV1UmCHdTRmC59gYqoZeSvsshW3HlIzBYT8l9X1XeQ8duUt5ArkLSdn1t4tyvDSXr907lS4+XL2MjW1VNU9Jl2+5d+Rr5C8gvrTJaGFlQ4iPSklQWNDxc9pu+8u9QCkorBY0VjPqjqWLHgskmxMUNta/ZghOYg4q9opSnqCSIXhtYzCrPZ2nSOAmASyEmoRehF31+GvPVVPx7NveoIhBUIkoQkGmb4LfiZwuP+QAAhLp3uYGShLgE3jgoB4iQZxKx5+rgBiSbJrUVK4DMoSceQHZVxEoWJoOmljZAWRxplSVleACt/ckx4/R7WU6f44pmmW+qHaZ/uB5LiCvIoI9RMWVe9UfSPX6c46JBKKKxh4CX0sPDerQKXRtD4Ox6aFLa6d7jc9VOypyG2AQTobDevjMgUYAzjOsgZhRn4NNNztWtyf6ZZQNwYBHwGrfI3ZONTZOTjLGnhA29w8xTZcZ7AC6HGwCHSKXjdwonEgF6lnYtCdIndLxu3y/YJMA8m85ZQ/gmcnh/MbPQ3I4ID8bID/ti7Y2BAPLk8Ot6QnU0kFfNDiENFZRlSyeYwKm1E/Sh+fup9YnqGzMmjodNNkyZeOT8pysHiK3ixpIjL8qNxANcNlCtZrG9xOX0TkxAVZyyNQ9BcSfu5VSOrbSPjano3VK5cNGG2X6l5LD6fq5yeFKuYAOmRyuFhC2bXK4fK56Wutt4bxGAjcJYEsocrYzOPYpFdJo+xRnrbjEOZx3be5FY5BLBLkX0Iy5rdonSAQc6uL/6WTaBZWRqlFoANhz/OISat06a26hlsMfsw0kZf0NgzCyxQvWhMB8zJzjT+nklEQgXD+AeFKYdiU1bqE1qQDQbqRWjRPnB9JTvCxnH6UB4dCj+siVUw4+2AF4MF5CRtVUkBisp8/AUQ+pVWNrnw9csl2cNZy/HV02WIPDX2MPY+W11hiZVyQCJhTPCQ5g0kNrQ3CcS6SHwoI62dsAYEgjg3WkcJHLZyEUbOoKnDqAjCDA9EXartvS+Uj9TxXiTbBGXvsCrdG3mgoC6RigHvm7VnK42LecSsLVncPkcCtogZijobodCNNgPljd7x5AMa3FtaVdM/WW6w9zSVDYQmkkcvxIAsNcpTMEAwAGEyMQJuA4PiBiRwgkKDxIACIRSEoICSIzEkApEIwQ27XuX6qDHQBxTBjLCz8jk9o70pdipA3xDAoBYba7lwAsVz4XtOE3DxxDqJdyzQ4Qy6TO982DwqbsADASQ81oXDooZtugsF2AQTg9h0bgJgFMwdIfwUyOf45Ofs4cYZ5EhTKjv/Qr9F1kByj1NQxq2kbpOHNMZFBny15Lv9CxX62yAYSUEL27MKS5gHK9O5e5/WT+eNl+Gdev+xTGF3MD1da1Y9Xj1nID2fY8V1Bsm0oOZ2E8F1AX/h9N8Kb09NYFcu5ZA1PZPpccEpO7q67LjYsUMnUdEhoB2AH2Ka1V516gUsoQ/hTsoK6ya2XEaWpuayQO85AiFIpYjJwYpg+E0XEBWg1UUgXl6iDOEX3y0KaO03qN+ClR80xLRFluoAx5UxGh2z4CFrmV2pYkhwPqRmMLuRHXqn20+mg6JqCUG2h0flCWG6h+JrAdu05MACO+v2YEvoqgdSjnYbwpV9VAu65jIEgX8rcoCysVkNchqA5KRRQQv19E4gJKGc8GAtEQBjOnuYFYJ4sDQD0luYIS4y9ylVBJZaRdN1m/iwpHr9VHjkDIZOX+1dxAXs+mPYJKbqBJKgi/yAapEbiUG2jMTbSUG6g3faxap5QbyK3VF9JBdOF91HIDxdxBqZF4Y9B3tg+eFxOwO7RUEA0aNGhwZeHQUb5z4PyRpPMIO6hGLjQo1Q0X7msqJi749c9VR2U2ADLzzYGg4vEsuXflSO0A6lyAxCYQpYKaekj/HEoxAelekEgBuarJ6s7s5dozAaegkqodElNqL8UE1Lx8puIELNTSR8/R60t/UdEsjQmY4zVUshuUVEK1k8K2geYFdEkhqCbsj2WGyibzjrFqGds+NY/WDoysP7XukrXD+DDQ1LGZK+jzuVyWtfV9B+gUzxii/j/7c9EO8ez6UhfVH8TII4MV8tbIPKh3KiqfWf/r3weP/yQyD6Ls2XIVUVT5cGoS8WtHnT4n5VJMwJhnUEkFlKqIrJdPl6R1yD2AtHqolCo6jeS1XkFatTMUsoHasugUg8rIxCzU97YOQm5eQA3WlSIWzLXYoLzDPu1amSF4CqyXEKAMwapd9882wbkUwAwaBucRJBIAw7F5OpkcIzkxrHSoTNVLyHLvKBiQheO3UoWVCPx91QCN3AgskkLN66cWCZy/Pgr/Z2sU5qp5DWmOdo4h2HoWpdLC+BnAlnuucdNnFQkMRInoPCWDawTgIsGuv40JJFzEM7rdqH9mE5ACUQhXLSjMZAfNoEP49VZVRAli5zRIjJUaSCP2pE9dDRRUQQXiEB67prbhQtsUsEpkVyM6yaNTsVxC6rUD40O7qqtnA80Jgi1nwVjmFC9gm8ygU+3Tx066jKB56og1iQDjfHoBNQKwBJZ+mxn9M934XE65MM/SdRMYGzNFOArEIEnfPLbuxN60HSBEBetIYJ0iQu4Hdak4AHtIfGg3toBq5O9IXUDodizS+vRB9eVVM5U18peU12+TG6gYJTyh869xqXNzA03p9XfNDVTytpk6QjL2s2uvZws4j0dCNgKwEuyLcGcEQsMUsUjUKAsXthz5gnGlcjJXgsQLZdOfCZHjJ8Kiv5FEAhhSgqDVQIkqCAa5z0wNUULahuOP4wgZ8Sjun4IqqGQ0jucDWI7fvoaUi59zYPxY+2hyt9CnbvTNg67GD4y30sPkmcC+ferAeEkTYfe/Oic+k/ufsy4RvZiI7iaid6q6/42I3k5EbyOinyOih8/ZViMAh4YtOfzVxsscY+UdwUoEU9KJVgO5Ppow1NagKGVUfsVk1T4mQljUQBliV5JAkdNXCF7UQJrj12UrHYxx9Ll3UPm5wlyZbg3Ji7TEYTQSWM11qNxAJT39rrmB5MD4dC8zPIMKyHfN6FwncK5mA3gJgCeaun/JzH+RmR8N4KcB/C9zJmpeQEtBvs9KYuF5BEYF71Lazqo63FvOs5IbKJQpHSNtYX0fFCZ/m0SMxGtIPIRGooKT3EDDEDyKaODkzGAMFF1CfTt10UuILfLm9H/r6UNsfiYFlVCaaTQ+V9HjiFVuIPs7NB9t7Lxg2y6eL0mQGEdvGDkvOPH6wXq5gQZ0ENesDR0mN1BE6tO5gdYABnA6rMNvM/ObiOgRpu6PVPF+mImhGgFYEZhyBDiXYGRjNVatzV2YA0B0CS1h6Fr/kb2yIOQ5P6nSHGr9/Dl9heb+wXFvMiZIE64/zYkNEJfQcKiMpySCROUsAOGQTWRwUNMIQeic62RGEPQzF9Q/4V4/OwunvZCTSPbGsEgfNQJBIhFwUnZ980jgOcni6pHCZZfQISDpcrpoICJ3G+lbSzUNWu4SWjowfmMIDrB2JPB47MQaQETPB/AMAB8D8CVzxjQV0FnBmr+FmXNN2hNKdVvu066VrT02NyHYBLQdINgCkkAz6efVQYXzgnVuoKD/90Fh0Ruo7BFUUgONqYTGbABjrqOZykirgfx9GSe7F1lzCY2voa5qqbmEzjlPWJfHAsX6kfFTuYEGYwuYowYanR/zcgM576C1VECLbAAPIaLb1XXrrDWYv4WZbwLwCgBfN2dMkwC2hS2Yt3O1/sT4qhpI2pWahjEinZTWMRIBBLGrhUn3LWwk/O169U91rwPHuQdyZwfr9NBeDeQmZZcaWiFvUQMBUS1jJYCw9ogEkL0ejfhRR8L5GEHInAb/mfWCuk5z+bas7gemIFUtzQ2k+8T7aAcIqqJQLucGcm3l4yJrB8QcIjfQWrDAx/8eZr55h6V+HMC/B/CtUx2bBNCgQYMG+wbebxwAET1SFb8MwHvnjGsSwMqQ6dUFZnDsU3aAJXtIuENrE9h2XcttFtrCslr9QPGYRpm7+p7UnsK0KjuoswdwVPVoYzDgDbxmUjkAhvzhMHLEU8egTTQCy2Ex0QgsevX4zMQMFiNyp+wAqPxfUQ1Zg3HsE58rMxJb3KB0/vFl+SZRGwURo9CupAR7vKQYggGnoxdDMKAkghFdeWLkZXNMoxwSXzgU3nLg9kjHPFMomSMh0zQSpQPjQTBG5NwQHFvWEfMZ6xmBieiVAL4YTlV0Jxyn/yQi+iy4v5TfA/A1c+a60gRgjmF1r7AQwY/ud8Fci5+bJtQ8tf1A7UmVk/WDvp/Tsh8YVESCFwnBEOzKFD2BBOLfbwoe+UtuIDB7hO/nGpwKKPHZH9TWBqd60mob7f6ZRerq96URPxTeDkaOFOEHDyLV7uwAHMdn5o5IrGz655geKVX7ANEIW8sNJOog2z6WG0jqpVxG3JzdA/vJDSR6vuChZFVGolOkFOmvZbgVG8AqczE/vVD9o9vM1VRAu8CedIWL1t9lD3a8KWtcti2wnVPup/YeED8l5cQQbKOAa4fEhDgAVq6g7iIfHCYGYh0ZHLl2DgjaGoETAzFSIpkYlBk5AbVzLYSU4FhqAPWykQWLjQWObZMaYlYqCDb5epSBVdpq5wVvmxtorZiANUDHZYxdh4RGAA4JM74tU/J3uzuSX7DuovlI7bVwn8+vUkMYgpBcJjdQtv/g7ePX6pCknOAxYmByAzlPIESvIM/xw0cG27TREVnHyGBLBIoePxWJYIwgTKaGsC+5sF6WGsJuh2MfwEYCx6nnpoYYi/QdJRAjqR567tKgMjO2lBtoaYroHlTOS7QyerySyeCI6CYi+g9E9B4iehcRfYOvfxARvZGIftv//0A15nlEdAcR/SYRPUHVP4aI3uHbXkDebYGI7kNEr/b1t9kgiTHYF8GdPe9CBJ8RiMpcxT5T7bbvkrVra9n+av2pPSb9JDlc1p880RmJCk7SQSg3UMkQqpG9zg0kyNyUawi/ivQrYwA9jlLiUQP/QlJX0ZzTL3L4jKRf2p5LB7umhhjj2ktIP0H8I6khBk7z9w+cJ3OrpYaw5VJqiKS8lhso79cIvC0cQgI4BfBNzPzZAB4H4NlE9DkAngvgF5j5kQB+wZfh224B8Llw4c4/SESiqnshgFsBPNJfEg79LAAfZebPAPC9AL7zAM91fmBtCaGGmLeEjFktrWc59kpuoChlmORwKleQjg3IN6OJweAMwiZGIJQLqqCMMFTUQNavP0P4SMdluraMoFDex/SvqoEUEre+/Wy4/n2khqj5/Uv/JMfPiOTgystTQ4yVi2NWQvoWrqQKiJnvYuZf9/f3AngPgBsAPBnAS323lwJ4ir9/MoBXMfMnmPl9AO4A8Fgiuh7A/Zn5zezCGV9mxshcrwXweJEO9g4LOfi5/YtI8oygincCMi7cV/pm5QKCn8oN5CQCpMnhaplBkweJUoC7EM8I6GNbQPyFQLGUi49qoJpkoF9iKWjMEoeamqemBtJIW+YsEZRtTgtLCEIBMU0hfNs+mjtISwcHShE9LxAsTxO9HRD6oZt1HRIOuppXzXwBgNsAfCoz3wU4IgHgYb7bDQA+oIbd6etu8Pe2PhnDzKdwodAPLqx/q0TX9R//eKhfpNpYAIvmXLj+6NwGqWbzl9oL5RqxqiJ506fUb2mK6GxPXeTuOZEQYiRwog6yaiNlBA4G4FD2ht6QLTTl4Oekii4j6vJlDcZJ5O8MmJIQ2BAI66UUuc60LPcCGRde+fh1CaHMtZfSH5eIh6y55HhJi+Clf57zf5o4rAGs9njVVEAAACL6JAA/AeAbTeKirGuhjkfqx8akFcwvYuabmfnmzf3uN7XliwW7/m7mSBxLpR01rkaMkjpCxu1HHX96RWmAkjFcUQfFswK0CojDWcGkzg0WI689N9gZiY2NoGAbCNsrEQd53Ak10JzUEBkENVA5NURSHFG1aDWQVQUdOjVEjy43ChvOPRlrDoAZCimhp1JE1+beCdTPbuo6JByEABDRMRzyfwUz/6Sv/pBX68D/f7evvxPATWr4jQA+6OtvLNQnY4joCMCnAPjI+k8yAocg3IdlDhIY/V0axLyTdCJtHdK8/4mKR43vZF3fbhF/CYIKyLiBDoNqc+qdxEZgjofUaqBJjyCYsiEAuVFYrviiqt5AiouPL9Fc0Bw/pV0UUi+5fUL10W2zvX4KBMGWp04LG1PZCGefGoHTtSaNxH7+fZ4NcFW9gAguSOE9zPw9qun1AJ7p758J4HWq/hbv2fNpcMbet3o10b1E9Dg/5zPMGJnrqQB+kXk5Ld2XGqg474y6qf2Mqq7GuOza+rV1S8ym5tJr62tuPRlLaT8Y7t6sk3H+ABJDr7cDcIHr5xIx4MFfkRCIQTjxCuo59wIyRuDZHkElTt8QiKgGogmKWwArIUi1UfMk/XXRjCmNndT5Z2oaTSC6rH8sd8X7Yt+C7WCs/xyX0JIaqHRa2C4ghPa8GYEPEQn8VwB8JYB3ENHbfN03A/gOAK8homcB+H0ATwMAZn4XEb0GwLvhPIiezcwS1ve1cIch3BfAG/wFOALzciK6A47zv2XPz9SgQYMGC+Dw+v05sHcCwMy/gjqf+PjKmOcDeH6h/nYAjyrU/wk8AWmwHJhS/XQA+Wq2jQp1UyCcPpDlKcrWn0gNQR05zl3P7Tn8JCU0EcgyciH1J9wcNIB632nDwQ4AALSR2ACvggn2AN/uD4sJWQpMbiBCvC+ljSBWEo/h3OHf0+zUENlfmOPcdY4fpypieTWJKsimhrCpJrZJDRHGwh0eM5UaQtbWiLKjs0kNIetArbMrDEP2kc4crnQuoNVhKWLcAZEWx4217Tq3qh77GTOlj1UlLKz6a2RTeydK9cPEoZznBkJUbIpKKPyVK+TpKtwGev9H7k8LI++KJ6qgkDxug6hqAQLSjvl5KEHaGqlHBK7K+jkZOcJP1EXuGeV95nmHHHIP79t+S/+Oo3onEgdRT5BC8JogEFFMEe03N5ULKCMCTAmy7cDRGEs6zXOeJrpDr4hCmkguHhAj5fQEMNCADZThl7qQIlrmTxLBKaS/phJIR1yfJ2ipIAzs6xvNntfozEf1/FNzq7mSPiWbwNScFTtAYgsozWX7hLE7poYgNXfJ46d6qZ+82ADE/cKnhkhcMuTUMOMRpG0CgRCIp9BcO4C6apHBi8DaEDJbQKl/+jq0S2jSVfVbmhqiFBlcsgOMRQkn5cS4mx4QY72GgHl2gBKsbQc4j26gTQK4jLCFZFFVA9Xm20Z6KY3zZcVc+htVYQgECOGQFmHiA/eqMoOSiQjO/rQGBtDH/sEY7CWAnp0aSDjhnkEbGj0z2EoAYU1WaiApq+dleS8FNVAyt34RHF8cMRcYgYj0KVkITloIAhEFNZB7LekZwaLySVVE8Qdj24sSAEpnCKdcvF67pAYCSgfGpGf8RknCz50dCNMn7eAhmR8YcjXQSnBoF8850CSAtWGEq96l/xyOfDWYkjgmhmuOv8TQWIkgk06KkoRw8HGMlkBCaoiO1No6IAypOkhgYKDv/RWzhMLmCPIuoKQlBsX5O4nASwwV7n9RZHBxjrj3ycjgSW6/4BKqyyNc9zapIbL+mou3HLv19NHumwXJYaxcSg0xJ0Po2oFgQPrOr5IX0IWDjBtdcd7Zcy7ksOdw8Mlz6fkra+VcuTQU1k4q8vWTKTQiksNi1NyBQMiALt0rGSMxa82OPyISQIr4AWElzeYVy+7TQojOn3WQGADa+EPhRQLo4Q6KF6JiJABr9LXIPkP++jHZI9qlbGOYz4yldL4kREJ+D+qVaHvBwEgkAlBqFJbDYhIJoMDxC8QDYaS/IgiVw2Ki9JHq6TfwBENJE9rIu5G5jcRQsxlsjHF5LR55zfMA1oQmAewLzvpbb7F+pmefmm+ptFObSyOfrC1y9EXuv5AaohYTUD03QAd/9UOQCGxgGEnq6JA5FOkVooUBGyNASLn2fUQGlxkA9+JKkcFW518LBLMGzJgraDwy2EoIVZvBiM4/0+fvnCdoPBWEZAhdPSCLz6cE0AjACOzjWyRIVsOMtYpqoDMiNHPUQFE9M96elIEM4Seg1EDZPIQk8Is7AnddvLQhWESGxODbBykgpIawgWFJIBg7KSBJF61UQxU1UJEQGAKwSmSw/kgFwrI0MlgjxiWRwXNSRUwdFqODs9ZSAwnyt/PXxu0Mxe9VuA4ITQV0nkHrBGbAqOpKzRXURfJ3xJV2lMulfYW1Q0VhfbVUUQ2k5s6Igfwt+glCsldvJA7rd+pe6/8Bpw84LVEjv3Y/AJ2/AFDfA32XxwUM6mFYjffqoaiOUvfCbZc4fVVOPonmvL2L5+yfQyAG6r3CzEemvyasnI7Z9czgmhoo9vf1XoTrlHtn6czgzM1TqZCCGggAeIhqIN8+dWZwrgZaB86jG2gjAPuEhQh8L+tjhz3Y8bX5tnnOkblzgqMqlMQQitmh8Z5L7Th4AoV5Nh3Q+7KMS3wce0BS8vbsjcHluAB3oDwCQZAzgxO9vuyrM4SvQgCYKu1Q9cD8wLDwfsuBYW4ozwoMs3EA254ZDIpEAUDRI8gGhu37zOAsUMz31B5Bu8J59AJqBGACMmS0z3lnINIil78rodlyvMXd2ZzIhIvsmZlSxJhJJ6W9KYKg1T/SJEZZGpxdgBQBSMrZA3GUAgBnAzgdwBuP8MUoLOUhjQTWqiHA70MQdoE4ZAje3GuOf+vAMMvl6/UVh38eAsMA7xFEeWBYGsjVKwIh880LDBNpwgaG6XIWGLYCyDs9b9AIQIMGDRrsGxjgc5gK4mobgc/f98hBcbdzocpoqLmqhuhKe9EAXVt7bL+EjGt347aMDB45LczdE3jjr2AT6NyVeQNJhtDeXX08OEZ7BAVOX4zCKkdQdoiMv4J3kKmrXbnXD5nyyDuGnp/C2OQ78ZzDYnQ7ZWPDaysYdbeJDLaeQXMig0uHxo9FBs89NH4fcQDNCHwVYVc9/Fp7WFNFVFTL+P+XruNxUxJroHThsV5VGAIBQhYZ7IZQtAMA3omcYjDYZhPzAAn0fVxL4gJOvUroaHBqH0lExxTPCJDNqD9iGmKMQK7CMYjcImD7Po1KKMypXwSXXpyah0xksFEnha7G6Gsjg0EOkUqXtSKD3X2eJO6QkcG5LWCtZHCHd/GcA1dbAgBmcdeB21wZprjwRWOnOPbKHIm3TWWuYn9dP7aI4firEkJN+qj2p8j9q70la2UxApRy/jYyOLiD+mjg01PgtI+RwSowTOIBAmc/wLuElg6TR8K9T3H/1TGw7XHvk5HBKIzX1Tw/MjjGAciVcvi6j77PJQAqSg/OLTM9LEaDjQy2Rz8mfQvcvY1C1rEANZfQVaBJAFcY1ubCD7CXGjOZzIN8rmDYLbRlc5s+a0UGAwguoSIdsLiBesshnSpCEDbgCQAQU0NINtDTAXTUOVUQBOlzcPUkpgRZawNxLVV0kjdIlQWhp8ibQoX0DUZg0y4eQXF+3yYEU6uFgheQWpugkD8nZTIc/8AE4unIYFnPSgQ1l1DdNhYZLGsC8yKDgVxisEbj0Gctxk8RyvMEjQAsACMp7xdmIOlcTSIN2823aD+mHLawzToKFwUPIKUKSb2CKFUD2fTQBZdQ7jzSJIpqIMC5hJ52KQEAIiYcXFQwnXq1wFEH7hndqVd1bHxyOEGcA4P6eAaBiwnwfQdXb1NFa84+Qei6X6H/bJdChcTdeyz8SESa0CocRuIFJOkgwtqkkS6nUawFAgEg9xIycQLSxyHjQZUVcQn93FKC3DViX/usgNWgEYBzCmOIc88wyWVrWIhcR+ceQ7IGCZfmne2GKhyn6ZosVZgr2fvUcxMCYsvOCugo+OkHF9BOYehOSQDCBotdgChKAYBXA/VgExeg4wC87iS+KM1VLwkM83WJJJVwkRwRt19ASxBhLotzuHAv3L/i8E1zouohsnEArlNwE/UiICmENxYYVnMJjeVcpz9UbQBeWlA5/QfFyXfo81xCTkz07fWzAnaGM8AvU9BsAIeEs2YAtCplAVQZl9p8Xg0zm+EhJHNJTh9rF0gun+8nW0eNY0LwBnL30RuIyd07jyDyl/lz8BKAswWcgk69HUA8g/ohBob5w+OL5wYU7AGC3BMvIc7LAUrEQc8V7uUFuXuybepFZfln/FibRsLq+KcOkZ9qL6eOsB5BZS8gV5/PFcrWXrDGmcFreQQVfgPNBnCeYAaHvYhjXwBF9dJMjn8RRz4HauNnSEolpnP2/MifRfT/iZqrcGQkq3J2VoDYAMQjKOgLPPff+XhPr19mlkCwzieG8+2nPXDUgU69BLDp8sAwQfhQyBwIHkFBO+GJWJAQxHtJSQwHOTEMqp+q0zYC0flrlRAre4Lo/LcNDLMeQSAkgWF6rEgdfZgLmY5fI2xnA5CX7XT8rj/8+N5JE3q8OStgFcgI8PmARgAuGsxAwhaqtgutTagRoEp7kegU9pUYhH0zaRUPDO7xBYaoE4RVTvcSyh6JBtW2JQaEEF9ANhK469LUEMYgzDw4NYYyClPPYEkFoRLEAQgG4agSUuqhTnHicAiaBvNcCcJWl7wnRVAktXNiBJ4CjYRGUkVDkH0shvV1f/GGtWcIB9dNZSMoIXxLENx9nEPKogYCooE3GIOLaSG0Tj+3Aej+jnjZ9lQNtBacx1QQTQW0JWxNzDXy2kf/uXOuOb4037b7rsyt1T8lSNrVJWogdHBuoSYwLKqBPDEgAlEHIuftw0oFJG6h1Pvr1LiGhmMiJXsooorGB4rlKhskdbUrPmjaX9enc5L6v5AqmuXl1lNFhyKnfea4gdrAMdteU+PkKp5O1Xcoq4/S7J4S7NWrsTq7qO3vXEN1ttHOjzmf2UCJ6MVEdDcRvVPV/Usiei8RvZ2IfoqIHjBnS40AWDhDKW0RUcn/Zpe1V9YP/XT/CaRr1y3+hi1Stvd2r6q+rOP33HqhnOj8u/ysgNB3Q2qc1Jn0X+GMAIkN6KPO30cHUx/tAVbvr08Ti4QAGcKfixxWjwvI5vBVxj7Aqk7a4ytK8+dvc2KY+P3X4gI06LgA8du3c8fytI4/LVNWtnXbglMRTl8z4CUAnmjq3gjgUcz8FwH8FoDnzZmoEYAGDRo02DcsIPKTUzG/CcBHTN3PMfOpL74FwI1zttVsAGcFhJmK2wOsP3MvowZwpcu29VwbU5rbzLF1YBgAJh/g1MW5RCoAAHQdSFQ/gOP+T3vtZA6wVwMBoNNT4GgT4wJOO9BmAB25/uIRlMQBKBfRkBEU6l7p+KG4c7EPVAPDgnqHQ/tWgWGQdZTotmVgGKuHWXpkpIDEAdjAsEMeGRmNwGlg2G5QEpf3Bv8AwKvndGwEoAQ1ZGagalzdAapIdgaSzsbOfI6tYGI/dgtL5w7jOS0Hj6BAAHyFWii8B68GCoFhG4+Q1ZnBvKEYGSzEQNxBOwrunQCcGui0B46EAGxAxwOo915BpxzODQaUugdKdaNSRWuCYPX6XCpLP8T7zOtnChQSj78T8zGFuGwZGEbsabL6QdrAsOR8AGDCQLwsMGztswJWg/l/hw8hottV+UXM/KI5A4noWwCcAnjFnP6NAJwl7CoFLETwkxx8ZBZHA8NmeQSZNRIOXzVL12QpMxfDnBhm8ZVIARAEH5GsTn9AfiyH/AEdiIaI8L0RWMpEHRgxEIz73kkBp85OQMdHPj1EdBvtegb7SOGSR5CWDhAZ1fhcKoo5UQkIQUBa1id8pS96RmBYNla9c47t4hGUeG6ptbcJDBN3T6AUCbxbYNigfhBbB4YBkLMCVoP5f+v3MPPNS6cnomcC+NsAHs88z+eoEYAdYRSp7jjv3uMCdpUQSuNLay0ldAmSyteyKqBS3EBoU2UW5C+5fAjgrgN8ugZ0HWizATrPAQoWChKAjwkIKiAfF+DjBOh0cKkhkrgAf+9TRASJoPO+9ZaDNyofSwCY0rJ+r1mcQHghfgFWJ4ZZYgBKiAUpX92pE8OIKeYH8otvc2KYvtcqHo3QNx6hd8nBL0MqEVCUAqybZznbaGy3EoM1Gu8E+5DEPRDREwE8B8AXMfN/mTtu66cjoufM7FdyWfo2IvoDInqbv56k2p5HRHcQ0W8S0RNU/WOI6B2+7QXkFb1EdB8ierWvv42IHrHtM+WbX22m/cLCfVZVkQZZZvNr5FqYs+StU/zNk2pX/496BNm9k76oUKbYXyJ9xRtIuYNiQ/Hc4A7RHVRUQSZRHPe9JwR9iA6Gv0jOCwjZQW1mUI6cvIkQXuoNZD2C8rGCnGvzKo8gKyUYwmLdQBmAdfNMzgcAZeOmI4FTDyDd12bv1P31wfG9WUPm64MLqXf7BJlyl7qNepfQtbx/APhvvo4XEBG9EsCbAXwWEd1JRM8C8AMAPhnAGz1O/aE525otARDRa3QRwKMBfOeMoS/xm3uZqf9eZv4us8bnALgFwOcCeDiAnyeiz2TmHsALAdwKZ+H+GTg3qDcAeBaAjzLzZxDRLX5PXzH3udaCre0Bu3Lha8BSDh0Tz1uTAlBQ8dTmDgWpE4TGyfyBQMgAbVj1xGAyMCzo/AXxb2JZ5Qpi7/+fSAA+LsCVO3SbLh4huVESwCAIwO9bAthMZLB+brbtrN5dhrxdBy1RqFiu0B4Jw4xMoWJgtmuTIHdOyhIY5h4tzjc3MEz6pvaBTkkBMTAMnKp4dGBYnhxut8Cw1WClv3Fmfnqh+ke3mWvJ0/0RM/89fz0NwM/PGVRyWRqBJwN4FTN/gpnfB+AOAI8lousB3J+Z3+x1Wy8D8BQ15qX+/rUAHi/SwSpwhlJAkQGpceElznlJe2X9hOOuwcRcO/3uC5IFU2FvWVyAkSoIwMZ7ANUCw7rOqYSUBOBUQtEwzH0friRH0OkwGRcQzwuQ+igRaKkgxhAgbS+ohTSHX5MCqh/Fv6BISKy4lc6TxQb4exsYBiD49A9mjIwDyhJAuFA5XayS2ydIBML17xgYJsFhlx2WPOHzTflbdlz763zU2ouJ6IG+7gYAH1B97vR1N/h7W5+M8X6wHwPw4NKCRHQrEd1ORLf3H//4jttfGc6Q0GQwtRfhFC1yxjThCX0Mop5UCYV6Ssp2L+yZwxjw5RC9CwiTdk8ICAbxk1cJeYR/tMmPjmROooL1JYQgInkOQWOhzqpuBOlX1DtTKqD0ntKy/gwV9U7yUaoqICQIHJxOYduLBMEi99kI39chDQzb5cAY3b92YIxea7VAMPvtK9chYZIAENGDiehrAXwxET2WiO4LAMw8l6svwQsBfDqcGukuAN8tyxX68kj92Ji8kvlFzHwzM9+8ud/95u+2xnXb+feAxEtIFsDs/ZT08kXYde927gozudU6JY5f5q29H9VmJYFQt6G0XqeE2HQuLmCzAcklUsAwuCsjBINPEcHxSjh+LwX0GvGzU2sViECm5x+QI4oa4YApGwLhDMKqzXykUqZQjdBtu+X0ywRhvF3f20yhut1mB5X/S6kitBQwGMSvbQB6PYHVzwW2P8jadUCY84Q/BeChAP4FgH8J4GNE9N5dFmXmDzFzzy7t4o8AeKxvuhPATarrjQA+6OtvLNQnY4joCMCnYL7K6XzB2kh4AkZ/bwbJZvOX2ivlSfWRQdKZRFCYK1cByUXVMhOiTl8ZhEUFhI1SB9kcQeIiCsTUEGIPEFXQiVcHhYPkGV3vrkgMENQ/iXrIqIAyYmBTShfSSUQpIEX02iBcliZSI26xD3KkH5ozFREyrr6E8KtGYN0XaXuJMxdjsDYIa9Dz9QrxA+VD5EsG4Z1Bvtmc64Aw58k+mZm/HcCHmPmLADwdwI/tsqjX6Qt8OQDxEHo9gFu8Z8+nAXgkgLcy810A7iWix3n9/jMAvE6Neaa/fyqAX5zrA9ugQYMGh4LzqAKa4wX0J/7/TxDRfZn5J4jolzHPA0hclr4YLrrtTgDfCqdOejQcXXw/gK8GAGZ+l/c2ejdcNNuzvQcQAHwtnEfRfeG8f97g638UwMuJ6A44zv+WOftym0NFWbR9f2E2DhYXgHl7OkhcQGnu0lqk9Hdz1tIMHcN57gRXH6lTfeVePILkm8jxkF06l/j7c+c9gnRgmKiCAPd/Lyk+EQzBdHLi2o98YNipDwzbuLgAd8/jcQFyXoBw24MXOHSgWEFdw5SWrYF4Ki7AtXES+JVEAIe5OTGD6EBBez6AnB0QJQXnASSBYOIRtDGeP0AeCSx1aWBXGiiWpIeWVBEqp39yRgANMT0E4k9BVEVjHkE7wzlkS+cQgO8iogfB5ZZ4MRH9R0QD7CQsdVli5ucjNziDmW8H8KhC/Z8AeNrc/Zx72AcSnoAicTF7SfrI/CN7zQihRh61ddS0yTA1P1NaJ+qiMIWofWRxVZbjIpMDYyQ9BFy/oPoBnBfQ0VFw+xSPIJYjI0M8gLiBnoJOj2Kk8JFKEyFxAUIQBqRI2OYKEgRv9fMjZavPn3VgTJgrunU65K6+glYn+Q6lA2MSN9CRA2OYCS4QmuQzFNU+QDlvkA4UA4bswJg0xcR4egiJEg6zMUHCtTcY1rUDXEQCwMw/4W+/h4i+EsDnwbleXg7YRgrA9JgqUt0BRhH1mmttMV+C8GdKAfDVe4kLkDXmxAWIBLAhlx5i4xG8EIONlgg2wIkkXQT45DRICKSDxADQaQc6cmt1pwN4sxmNC9AHyHjGNnnuvcQF+PtqoriwlumvJpPvs+uBMcAYwpd24fqRlCMnbyQCj/zzuIIoEYT0EH7vyTnDKyWDOwv1zhxYlAqCmV++r400WBkM0p1SA40SLKNmyZD8CMEozauFh+oY3dnsJSEaggvs3gDYRHFBWpD2jTNYkk8F4VRAzi3UtW8cMTjyfyZH3uArBOF0AHgAexUQXTsKXkNhvEgLAzkJwNMO6sh5A3VCEAga34R7TSgpfgOr4tF1gEL4hssfxUKBQ2clBfi5NPUI0kLcnJ5VEsVpjyGVosm3xYqxRHHhPhAPl/ohpo7oQnoI3yORAsI9RRWPzf6Z5wqSH9Pah8KP/erPBlouIGA7tcsMLnknKWAbyeSMpQBgQkoJndL6jCBYwqSqgvon4B6TLrokDQRmz22OdY4fJQEQsUO8IVkcgdUZwCEwzCN4Pj0F+8AvKZOSAND3IK8OkgjhqALiIAW4slEBTaSLhpcGEoKA2O6Qf+TaySDpnDgg/TuwUoB+576dArVJibbEBZBC6hlBAIJNABTvc4Qf66Ws00N3YOf14/t3SKN9e4/8Q24g5wqW2AzS3EBdkgJ6TRXQapmlV4RGAC4YLDIyG8JW08tX20f2QBZhTKydzWG6LQJS4400olVAgYsVTlPapN3fB1vnhhzSt+cFHHmO/ujIIfeNZ+M3G2A4DecF4OQkSAEAQMdHYDEIHw3O8BsQPoW4AMBJAtpIHA6Ql71J6gh58MEQW2/PsDYCq/NXOFs9OAc1kK6OL1lHC4s0oGwGXuUTuzskTmoiy+UHKcB/KJszaPxA+fS8ACAi8WAErkoEUuvLqJwXAECnlFgFLroK6NLD2lw0FiJsC3vYz+L1gcj9je1F9bVSwBwvJMUMhiV1V7sV3YERFyA9iV2HkBOETaQmPDjVDAcundLzAo427gwApeLhk1PI4SI4PY1SAACcnIKORF1k8gR1jK6D4vgpxAK4Mly6ClZlQvAKYnk2i/CRlpN00cLZAygahI2Kh+zHUP9rAiIqotBsvIACQSBBuuXzAfR9nh7at3vqP6gfTJACkEsQHZwUoM8LsLmCAAQCkhKDYb1cQIztcMCeoRGAXeAMEXRR1VLjus0+d7IHlPqMvYcSEraahW3eo+A15NJIYDwNx2/dQtGl6iPuohTgyhW3UJEIPDEQLp8HjlIAABwfOa8gAOzPDaBeSwDRCEw9QBuo8wIYPFB0vZS9KxUPGJmbaHj2Acm7tl5FRYOwIOyiW+hEuuhEReQISDwgZtwttAMw+AftNHKXelNXShe9llsolEpInxy2CjQCcAFgKTKaUHUIbG0PmDn/XqHEhk/0tZLPlPpJIDEAq2UziaBCVDLCJFy/LG7KTkUcsEdSDhHBwsUP3i3UnwhGRz14c+IiVgB/YAy5FBEA6OQUOPL3xxt3hOSxdws9ZWDDIE9LYhSwR5pil1Q2Aa2uCvp//R6VyihISYHLFySvXzSH1zLtFmrb41ohHXR475zYAGycgHULxRjHX5IIjJcQCFjqFqqPhCydGob46teDRgAarAlV9VKJiO1JCgjrjxGHuVJAbXwJFJVgxAlI6UFY+mmJQBlWnUeQDw5DbAsJ57xHEKkDY3C0iRLAxhiFxR6gvYKO3Z8YnRw5N1A5LWzjUkSgF2nDSwEK4VPnpADA1Wu3z3BvkbImuIX2bLx00UZfrriFKq8gN59/FkrnAymkr8o2XbSWEJaeHxxjBcbTRffEsG6hpXTR1i3UrdUlv7NdoamALgpsw3XPkBwO6hW0NpTY8Bl9t7UHlKQAfa+YzTiP75ARpkAAKEoBvsxd5IrJnx/MgSuneHAMEHIDUXAL9QheOP6+d15Bvsw+P5Drewo63qA78QnljrxBuBeCgKAGkrVJooEBlyxOCBgUErZHSmo1jZIIXF8VA2FUPEWvoIqKxxIbeAOwVRFFw64KDvN71ZHASVI5uAAxeazSITFCFIB4nKR2C03dRKkuBWDELdR9hHWlgEYAGjRo0OAKAqNJAA32A1WDcKkOsX51t9DauoW5s7Wkeqmkk+05zRPEuo9wySoXkKiBAMQ8QTKXcP+D19sfsbuXuICjjbcJeEPvySnAJ9Et9Nq1KC0cHwU1EOBsALRhdKJe6imogYBoAwhBav5AebFJhjxBlivXkcIlvb28GsXxB4NwzS1UpAGj0lFuP0ENFOajOJ/o/zOJQH2nRGePPE9QneN3H3RQkk0aJ5B69ohBeCwuIA0Ku+KpIK40bIOMZiBLYEtu4DyogYBcDzPRd+9uoWz6W/VTQgDU4r4cbAAc1UCAUwU5xOsRmaSGOJLUD94gLHr+06N4WhgAnJyCr11zbdcdBzUQAHQnXVQDAe7QGJWXSNJCxEPkEdVAQHQR1Qi+RBAUUk6NwnosBTVQ0qZVPCqtRKDSCXGJ71WIgzYCS3RwmF8haWJFk3leniDrBmrjBHRcgMQC6HEyvwSJBRVTIU/QatAIQIN9wXlwCwUUwq31t212LewgBfgJ2BiEWffTNoE52UK1DWDoohQAuKRvIgUAIQ1EcAvt++gRdO0kSAEAgkGYfFxA17GXAlyZukK2UKJgJC5mC9XvzbiBRptBqveXtkVuoSIZaGu76g9wQpSFgKRnAuQEwdXPOz849xpCaAerMwFISQEAYNxCY54g7Va6vg2A0FRAFxN2VEnUYB9uoYu8gnaBEhs+0fdQbqFJfyt9BAKgFldlMQjbbKHRE8dkCz3aRCkAAI6OQBsVGcyDUwsB4GvXghQAIBiE+Ug4fgpSgJR1cjinDkJA+GIQzgy3I+mjrddQLU/QPLdQTV3UHL7DWLZQSQsRVELqGzljsZEmCm6iAnU30BjgJVKArC1EJo5XKqGaQXhX4HWFibWgEYA5MBOpHxR2Rep7kAKSfmPvbK4UUBtfAkUlIsdPWMUt1GYL9Rk+nT0guoXSdccuJbQnAHw6BGLAp6dOCtBpIqxbqJcCAEeMRAoAEE4RC9JNp6QAIHL8+r0ZlVCC1BXWFYS/2C00TE1RCpB2cErsFQHRtgAgTRwn9oBSniA3hoLaR8pj2UL70N+XPYLvg8GinC3Uzd1UQA22hRkIeie30KXz7ksKmDO3ZhbN3vbmFurnStJESOdAAChKAb485hYq9oDELVSkAMBJAl4KAOCQv7fasrcHkMQQnByDjrrULfSUAEkVsa1bKMf25MX4NjLldGzk4Oe5hWrqoteiKAX4ZusWWsoWKt8oSAG+YixbqNRFFdJK2UIBwASI7QyNAFxwWBuJ7gq7SiZm/E5eQRoXWLVLSSKw/ZG2h25L37mRBKwRmJGWgxQQBmDUK2g4InSZV5Aki9sEKQCAiwvQ9gBvBwAAunaSJIvrTgZwR+i8BLCzV5C8N5sqQiN1Yw/Y2isISgpQ/aWi5BU0ZQ8QBE3+x1GSAoCyV5D2Agr2gNA/9QrqXT7wzCsIMPaAFaDZAK4azJQCgPV/HIuMwlNja8+xJQEadRstzalwS4mOlCQCXQjGSvJ3CQFQm/HlMa8gYgp6ex7IpYkQFdBw7E4Hu+7YlU+MPeBU5wk6CVIAgGKyOJECAAQPITIIn9SzJG6fViXky2RtBEAMEDMc/CKvoCTNBJKPPMcryNoDal5BYg9I3UBTryCt85/jFZQEilXsAatAIwDnCGi6y+i4QyDsubCrZLJw/OhehdPUhE3PX1mrRggThF7bi3Qs7CVIAYW9RQnAqDKk2Hm9tfYK0sniNh2oZ5UraJMkgMN1x5k9QIzCODkBfeKayizag466JFlc56UAAKBTpF5BlCeLowLHX00WlyB0d+3kFURIVD7beAUBKj6ApD71CqICwt/VK0ikALd+2R6wM+j3fY7g6hKAQ8GeCMYc2MUraLYUMHO+0b3VpACoeoVbEoIgSCv+fadDWI+lKAVIByVdgABslLqBXLK4kDiGEVJGA44Y8FEX3UIlWdyxkwDCGcFAsAfwqcoTdHQEnPi+RxtnE/CnjfGGQ8poADFXkBiBhZs3yeI0x8+6LMSATP/QVkgTsdQriPX4+FGEHgTbivEKKqWJsF5BQsckTUTtDOEpryD4MDDrFSTQgVC0B6wAzQvonMFBue5tkO4+97Pj+DlSQNJPI5tSe21eCgxmvoaZMtyPqY8swic1lVcHseL4Q8poGdtBxQX4w2PEs2fDwNEQVULHx+H84O64x3DtWvQKunbNGYyPPUE42gDHG5AcQH9E7gAZLzB0nbcJqCMkJWW0bD3RWAhHb+MElDrLSge184YZiFKAv98lWZz1CtLSQy1KeEmyuNRI7KOEhaCQxBI4qJ0hLHOtmRO02QAaNGjQ4KpCIwDnD86U614bKuqm2W6hZvwir6BdVV12vF1LqqnA5VfmKqmLmFKOviQRRJ2/53RFoeyPXIwu5P7MYDnG8agDhg2o92qdfgBd5zn+vndeQX0sQ4LDACc1nBwFo3B3MvhMpGJgZnQ91Gll8ZK9J27rntvXOv/SGcIAgmoo6OQH37jULdS868RoDDOe8myhMpeogfzSiVvonGyh415C3k00Jn0uq4EAaIPwztBsAA3mqoGA9cXFXbyCFoE8o5p7yi1029gAaYJqpsK9Vj2QGsT+PqqAAOsVxJ1SfXQAbeLJWOyPk2Svt0fPLh3EkUdPR0fRI+j0FOgHd2oYPEE4PXWuoQDo6MgfICMuph2o90ZmwBuBCZ0cPmPOEAalZwgHA7gu6/eq4hsiQvddvepI/xaDGghA1S00KvmNzYBBSG0GDFLNqZE3IQgonx5WSxYHLtkH6sniwF04MwAA7Olha9kAhLc4b7BquuuLCkyKM1oKS7/szP5b70fW2AXM+GwvU+2Fvsk7pry92t+uVRozBRQvNy+5iNqkrnB16v8ODrGpMncUcgUFg7D/n4865/PvL7ruGHR8BNrEM4VdbMA1Zw84OXESwUmvrgGdv6hndw3u6k5TKcDVK6nAM/EkundG3i5lVv8LMlft4vWk50ovcvEI4cXlfcTVU+f+j6YC8u3xfmByFygZC7g2KQ/mknYBnTxO5hu4Q88UPIMGJvRw1+Drpa1Ht+q5wLOuCSCiFxPR3UT0TlX3NCJ6FxENRHTz3C01CeAKwVXwCrLaKCsRzPEKmswVJAisA6ByBeGIwEMH9BIYxjFK+PgY8GogN3hID4+5duI8iLzRmI7caWPkYw66nsCnHFNDdOROFPN7GwiJBBCOkJRHE4LKsV2rg7QRWL51csDLoOMjKEgBbi6nJrNeQcl4xP6iMpI377pHjp7URm0uIGJafIawSAFAPENYyoPn/jsTJ+DaVKTdCrCiF9BLAPwAgJepuncC+LsAfnjJRI0AKGheQduPH1Vdqbn25hWk282yVa8gTstaUgDg4wM4nVAhupAr6CiKLbwxuYIG/yfW9842IHEAfR+PkATi2QESKXy0AfnYAADLcwXJ+cbyLNZNlFCNEnZjKVQIMUjUR/qjdeO5giDIXnkRsRCF0D3+IPR5wgMTiAidyiM0IBJhR8MVsvd6vhJR0O3iNdTb/qTyBjEhvOA1YJe/Rz0N85uI6BGm7j0AgpvwXNi7CqgirjyIiN5IRL/t/3+ganseEd1BRL9JRE9Q9Y8honf4theQf1Iiug8RvdrX32ZfzFI4U9XL2qD/4BUUn7FWp+ozVVlh/kTNU3kf2RwlsOPtXvwlKpxR1RoV+oU6yvuJiseP4Y2/OoLYBYLqp4uXHB8pKiA+3ih10DFwfAy67jrQddf5nEEb8MDOLuDVQfiEv05OQaen6E56fw0ufbS/OlEJBdUQxi/h8PUlL7GmEhr8eQVK1RODwdR4087644BiZlCtMjIqIX0NSVmQvqhvRIpQqiRLAIBMJRTHd0Zt5CUB7lRbvF8NPAGdcx0SDmEDeAmAJ5q65wL4BWZ+JIBf8GUQ0ecAuAXA5/oxP0gk6bDwQgC3Anikv2TOZwH4KDN/BoDvBfCde3uStWHG72sn+8TEvBmMIdGxsWNjZhAE227nn7NeYjMwCJ/NfTKv6Z8TCUr7dJTYCKCJwMZJAOHqOvDREfjoyHn5XHccLrruOmcT6HyCOX92AJ+cgE9OgE9cA11T9oDTHtQrW8Cp8wqicLH739sIYAhAUPPIVSIOApo4cOyj50gIgr+nrC3ep0ge0Hp/EbIYMO0jBEFfyG0BmmBkieOApL1X/QEoe0BqK9gZikSzcAEPIaLb1XXrOhvIYe8qoJK4AuDJAL7Y378UwC8BeI6vfxUzfwLA+4joDgCPJaL3A7g/M78ZAIjoZQCeAuANfsy3+bleC+AHiIiYtQy6cM8iDW4zQ5RkV+1/0dRBq6WKsP2RtjOwc8I40ZPMPkBG1C6cHiAzbAjdgBgodtRFT5lhA+qHqOLpnU0A2i207wHvFYSjI+CaSjV91AXXUADgznkIhUjhDSVuocEGYCOFlZ7f2gMCERA1uBBkJmCIthL2nlI6dURiwGSTKyh52dJXewGptkqUcPASIhrNFVSLEp7KFRTKI1HCu8KCv997mHm2IXcXOCsbwKcy810AwMx3EdHDfP0NAN6i+t3p6078va2XMR/wc50S0ccAPBjAPbtu8kyR7tpQ2E+V0JX2bupKunnXUJh7i3cxOl7/TXJax4VmMvfJENZjPWLT63iuXy8WjcCeddWpIjZKJcEEFgw8dMB1G5cwDnAxAv2AoBe/di2oglzZGYUpEIQO1HUgv1a3SY3C3an/IMqAPRCh89MN8izyaEbiWmQUHgDCiFHY2wPiQTpuX2kmUo4pNRSxEPdQ4d/EBVQnkhPXUPkWuxiFJUo4zR3kM4Uq99A1oKWCmIYSubV/07p+bEw+uROlbgWAowc8cJv97Q9mIMmdCNLSeWfuB5i5pwoWDmtPtdf26nFPgvS1tID0h5L8mHwh+SFZCSAJrnKLV1NFbKKKwrV30Ylks3GpIo59334AhiEi/L4HhtNYPj11rqEiAVzbOAKwtlFYnlNJB0HqSiQEChWTRmGQK1kCod57ySjsRnLRKJwQhBGjcEwG5/tz4fyAglHY5g7SbauAlnJ2BCJ6JZwG5SFEdCeAbwXwEQDfD+ChAP49Eb2NmZ9Qn8XBWRGADxHR9Z77vx7A3b7+TgA3qX43Avigr7+xUK/H3ElERwA+Be5lZMDMLwLwIgD4UzfeNOtzHFQKWBvplubHwj0tHT8lKRTmSvpYLG3JvPRH2p5wtjOfz/59s58gKCcET1hEqSKFBUG5stNNBIlB5w066sBOPHBT9UcAD6BhiGMHhj5ABpuTeN7w0WlygAy83UFyA3XecB0IALlModkBMvpZtDpIt3niYKOKpUNIHa1+jLMihcWddvBSQeivPKv8S9efcOCI4OV91w6QGYiqUoB7bZwgdlH/xIAvHTQW1UGrwEoEgJmfXmn6qaVznVUg2OsBPNPfPxPA61T9Ld6z59PgjL1v9eqie4nocd775xlmjMz1VAC/uIv+v0GDBg3WBgLOpRfQ3iWAirjyHQBeQ0TPAvD7AJ4GAMz8LiJ6DYB3AzgF8GxmlmQcXwvnUXRfOOPvG3z9jwJ4uTcYfwTOi2hVOFOuuwJrSyZVNRDyvrPsAbX2pVLR1PjSHilWWQ4/VVTUVEKKs9VrkHD3nosfAHuADKAMmIk6yNsDJHX0dRsQHztVEOAkgWEIqSEAgE/UATISQbzx2UI35OwAkhvoNGYMBRyHnaSK6CmTABKpqSAdhFQRoi5SKiH/SPKaXJOWhLiQKiK8mDxVRHgv3eDmC9+1fHZALVXE1AliOotoLKepIuzZAavBOWRLD+EFVBNXHl/p/3wAzy/U3w7gUYX6P4EnIPuGg6iD9kQw5sAio/DU2LHnKGDhKa8gO39tvZI9oEQEkvtU+5AYgRMVkNgDRJWxsYvDIb2NqD40AfD2ADkA5gguz48kg/P2gIA0T04BZrAcL3lyAmw60EZsAATqOnQVr6DEfRXO2ElDfBlFe4AmDloFZDQhYg8IhACUqHSsVxANznArRt+o+0f4P/xMxB4Q1nbI39NJbDpkXkGdIhDbeAXpA2Uye8Cq6aDPHwU4b0bgBgLbIN2V5p9lFF7IyY/uVc01ZRQu9kfannC2M/ZYMxoHe4CSKDJOuVMDOkH4HMqDeO2Ie6gYcZmda6gghb4H8X1CsjgMzkOIT7xEsHGpIXAUU0V0XReikLsT8obg+OxyjCTgkDn1ce+jXkCeOMT3gOhGCgR7QHaAjNg3PMIPAzoOnLx7dreBxKtIGYxdv/jS9SeM9gCRLigcI+nm3s0rqGoP2BUYzQvoMsC+kO4u/beWTHaVNgrjZ6uDRtYeU/lUx1s2X61fE0pk6nGJIG6G5J/A1KebC3EB0mETm/nIE4fAqXZAYhQ+dly1pJLWyB9wKSQ2yivoaAN0HTpPULgj7woqmyHwKYL7pDMKc/Aa0gSgI5VLSD+jQvijB8gg/WZBHTR1gIwYlfU7F3VQ51VjIv0pg3JJHcSaIADYNlVE8UD5teD8CQCNAGwLZ4Z09wEFQrNPm0Bx7NR7Ue2WQDCQqYIAzFYHZQSB0vZMAkiCq+JmnOZAIT6fPlr2Mhx1CIzu4CslcdwxO2lgEC5ZVEEn4WH45BR0JEeEnQDKLbTrHNvOXgcU1EEhLsBdwtFW1UGqr5UQFJPs66WDUweJ7UTUQdYrSKQRdkojhJcRXi4ix8/xvTmdTkTw/eBUQa7MwTVU9rYJkhgWnx+gg8JSddDucGgD7xxoBGAHOE82gZ32MjInUCEEO6w1SiD8/c42Ad9XnVRYlBYyCYDzdh0pnHgGCpK0h8oHnFFPHDcwoeMOOJZ9e4wp/YfBp3RQbqI6UrjrnBQgRmAfJNadeMTp1UEB/3dwcQLepUJLAJkRWNRBqj0QPwA2qhiiDjIcvkb40AifTaRwSR2kOH4txoVuSuUDqMNl/G/Lnh8gSwcCoKj8PHXQCtAIQIOtYQbSPahNoDQecY4po/DebQJ6LelmpYXKc2R0IuAKYxMwHDEITu2jNhvUFeY0MWy6xBNGYgSCTWBwncNM1665/32gmDUKd92pIwbBQO1fhJIAuON4qPxQUAdpyUcTBLkXWuTVQWFvIiHosqIHtUjh8CMYOKa4kH/Ce+iCZ5BMpe0J/dBh0w1p6ggw+sSrKP4ohDCkRl9FABLjx8o2gEYALh9cWZtAbS+m/sxsAqjXVWhS1SagjcTb2QR89VFURYTgpZpNQAWFuQ6DcxkNRuMBOO1jemkvAQiCdy6h6kQxgrERKHVQ7+rDcwqytwRAvUNn8PbFwbwXXw6G3VqksC+TPnXNv6twxkIwICPMJZ5B0nkqUrhXNgVRB5U8gyD16gewVi4gApoR+DLDwWwCM/uvrZ66aDYBwBAG4OxsAqEzV20CDq8xgvM9s4sTCARhAOMkxA3IkZK45nX+nTuNrBM/UElRHdxjUIwUhjxfH58z4/h9WXP8MFy/6y9cPBvKueRMYYqRwoj9RBPWeVVa4qHKCHEEXTeMRgqfosMRhirXfzp0OOr2YwNIiNw5gatNAJZy4BNwEJvAGcIuNoExldA+bAKy5s42gUK7qCvcSVhIVUAa6SdU0NVXD5NBB9Z/jYO7tErIPepJKKMfXL4gIKqDxAjsA8XC8bbUuW0np53F+zmJ47KyyomUE4OU35czjF13llpfjplB5yaO06eJ2TOFEyMxymcK68CzBOlzXTrYFZoK6DzCHogAsIM0sEQSWIp0V9rLLGlgoWSzuk3At2U2AaVtmGUT0HtT08gEmU3Aqk5qeYP8X17KyXrDMAAcIxqGgSgJKBUQ8xCQqnMT3SQEgLou4PiuI6CPKqngFqruQ+K4Pu4JCIx9BJEQVJmTOAGVTwnxm7GuKNkEdJsgePE+Ui6lwwCXCM93H4YOnUfg8dxgWazDBkPiFbTYJrAGiNRzzqARgAYNGjQ4ADQbwHmFhdzqHDiYIXZG333sZfachjvf2ihcUxnJeIrdR/MGGa2MdvooLVlUCan5GBTcG+0hLKUNJFkvdaAYI3oGyXOGZDxaipEOJ45TD7l+nIto8IfvgI4IIP8n7k8y64IRANAeQu5ycw1e4Z8YhfUmChJALV2EvJc5RmH32OUYAdZiiDJAW6OwqH+sUVj2JjECUfKaYRReC5oEcM5hD+qgg+j9ZqqDDmKfKPXdlgiMrF21CRTawlQWgxvMX2om0y0hCCR1asHk/ABkh8lo7xW9WtR0SF0X6vzAPJcMncY4ATlbQHIHXXP2ANI2AWWgThB8h52NwhmapPCPR+qIgWM1o7CfbN9GYXjXUSAahU99/6NuwKl/67p+Z2DExHznCBoBsLAyEdgalu5jJhEAdpAGzLhFRmHDXi+KEyhh4Rn7tNJEkA4Utko8g6hOBPR92Irpzx3F3ATSZ1CdiXPpQPCxP0xmkMheH40VThTz0kDYmyATj/DdYO8aCgB0CtpEAhCigq1RWO0tMQqD1ClZ8fmluy2z2lySN8jX6wRvNaOwa+Nxo/DgRmijMAb1XPA0UQSdmUZhGygG2JiA3aEZgS8KWF3ADnDwOAFMj1lbJXQQF9EKFq4ahf29nb/mGZQsQaZsthYRvqg3VH+NMwgmn7B9uNh5AIcMlq7cJVoVLQnovcUBQ7wAf6JYl6iEMqOwpBryBuCpvEGWICQqIfW/GMPj+zRGYUQi4DtAG3m1UZjFoDwzUGwAhXOEgTxQTIzCU4FigIkJWAMaAbi6cLA4ARlzYJXQYmlgJhGottWQfEFaSBB6oRwQlzDBdr24VDa9bUcXc/DPcxHVSL9ABAIilGFxc4lKQSSBMHgATk8jAQgqIP8nTwT29oDulP07SSnfgEgQNGSBYqEhvZc9O0banLTG6pupPEKRPPiSpIVQthMdHyeBYlElxMEzyJW9CkhR0rFAMW3MSNxDdwRCkwAuHizlwCfgoHECM4kAsK50UpUGlu5lhIsfRfIjOpvA0FXKATGRGm5VRmZ6rQuPdUqdoSWARKlvqaAfIYjOxwdE90Wfmf44GoZllgQGhbB6yQEN0IlLDBRUQhQ33sMbh2WyE/NgXiUU6IMc0WQ3YghAljwustn+6f176tj00SohIBiGAZVPyP8XIoZ9WVRNI8njEoTv27REUDtLYCfQLr3nCBoBmII9EAFgXaQ72h/TY9Y2Vu8SJ1Dby6jNoCYx1KQF1MuJSogUEVDPUVP5RK7V3wyG680kArNxxICmDpyeJwBFBApAHbn5RBKQtBGDx9YnPkWCJ06DigoOxmEhfIRgGJb30IFjjn15DHngApM8qj0h0fT7vqDEyKtzJDl1D4VzEpyHkJIuJPPoEJG2GIZDGSouwMcM9OYEshCuwbQfIzCaBHBxYSYiXQKXQiU0spfZKqHC3sa8gOaqhMLaY9KC7Y90LpEGEgRf2mtcOpUGnIXRFTqfMtlKBDIgJJFj30QoCQ8DxEicQ4pwTx3W026ipzE9qCYGYhwOKqEMc/sc+6pePIU0BGaZUq/QwnTw8lEopqmlYxu86ohNJLDMmhACv4lhUO+n8/tINDxp4FjyGAUPobWgxQFcdMhk/d3g0qiERuaepRKqEIGwl7lEwNzPjhpeohLSe6wsnREBwHuqqGMUE1WQWdSXa0TA3ZelgbA2eWlARxKf9olRGNekq3MBjQRBzQEgeAipOAF3cEoKWuPDSN8bCIl7rL4JuYPg0f2MtBEB4ZOxCcB5CAVVmvcQCgKAf4li7BX1kCYCWdqINYChw5PPDTQCsC0sRcAV2FklhAX7mNl/H9LAzioh0z6pEvL3s4LGbH+oufQWqU4E9L3dEqCkARU0luYZK1OmqBJykJoqC0RAcfgAQvI48OAuSSd9ShBlvkP8RwrhBx9KszdFnIjSPVH6nkogcW25DSF+pIDuw3tKj5fM3ETBCIfMh80wwhvrBmCgmPoHYiZxFZtuyNxEtXoIK6qA1tQgrAWNADRo0KDBAaDZAC4brCQFAOdTHbS2dxCwvTposU1A3Y/ZBOzQzCag91GyCeTLbW8TGHkZs2wCFLn6sD9SKT61TUCrg6hPvIIgM6p0DDZOoDPSiYrFyuMETNnGCdi+IgW4vhWbQFAfRSlAylEKkOcYwnsfOom4du0SJ9AH0crFCUgbVnIDdZs7fxSgEYBdYakaZgTOIxEAdlAHIV/jYDYBWdvq+Ed0Nlu5iZpHLdCfdEtr2gSIfMIIRwios3FnFFVCp73XO6U2AelHFAlCiRgkcQKQc3jjHm3Q2hy0mT5h+twJup9jE5ARnb8LKiGvDlKHzA+gJFVErwy/fXAjQhowtiswmhH4UsPKNgFgXT38Lv33ETS2mk3A9ykSqwJLvhc3UbO/sGf1R1+VBsZsAqQeqC/YBIYoARBxYph1Xj59QnjolJRNgMOD0Gnv2k/iDvVWmIBNUpO+uEGQsnrAIhGo2QaSpvQHs7VNIEsO1Pm9Dui6+B6tW2gpangNcJ+ySQCXG1YiAgJrq2F26X9mkcMFQrBV5LC636ubqH5GQDmYI1N7WJVQ9g3YvoBIjZxbZjl9hDyc6+MR5QkBistPjpdkBp0OAFQ0sYoT2EAFiwVIX5yOYrbpK7aTBnIiAGyhEsoCAbpABNzeKHELdYRAwor1U6wATQK4ArCiSgjYjxqm2n9fKiG9n21UQoXxU5HDtbZ9q4T0craOVYBSpvphiieM2cFxhrRMMV1DB/JZOpVEYDfTIW6WCBQyiTIwcDg/GNQDJykC70L+B5lIT+5enN5Lkq0BBdxX+Q1p3l533EkllPjPOmlg8DsK0oBRCcmTrxUIDDQJ4GrBeZEGluzlPKuEKuOnpAFghspnrkrINIc6pRLKJAYzTmkXMhVRMY/QpJHYwQAGUXro/ACV/I3IMfjKTZTFBNBzohJCzyClsA7xAZA0FHZjuUqoQ3TNjAbrebBXlVCg4FElpE+IxNAp20a0DewMjPy3fA7gTAkAEb0fwL1wEuYpM99MRA8C8GoAjwDwfgB/j5k/6vs/D8CzfP+vZ+af9fWPAfASAPcF8DMAvoH5HJDbpVz4BBzMSHyGKiGggNCBHO9NEQE1ZlTlU5g63Vg+Pgw31CA5bxg5sSgKG8rl3nGqelO5Wigc20jOgEw+Y6fl9gdye5FD4ocgEUQJIN73DumrYyd5GECiEaLeKEI2AAYQ11VC+pxe0cB05hjM4m+MS23qpSueH9hCJaTmgeyu4CUEpPECuwOjeQGV4UuY+R5Vfi6AX2Dm7yCi5/ryc4jocwDcAuBzATwcwM8T0Wcycw/ghQBuBfAWOALwRABvOORDTMJKEsHBjMQziNe+9rKNkXjUCGznHZkrkwj0vpCOL6p81HqLXUY7JU0McOkP1IRMUJ4k6cOJV0705PH1imCAEJC6i/4V8QAgGsBBJYQkfwKdKiIjq3MXomc36JL0DICc0OX2N7DxEtp4hlhJCHFeH0msg8YyYqAJ40KVUCFy2Hk5cXhGkQbIbm5HaAfCzIMnA/hif/9SAL8E4Dm+/lXM/AkA7yOiOwA81ksR92fmNwMAEb0MwFNw3ggAsBe1ELBnQiD9Z0gEa+5lWyNxRpQMISsieIOV144bkCWSrci8GuHrLYm9YCxuQLtrSv5+Qfy+uVNRtwPF8wGoG2LOH7EN9N5+0PcpJdPSAKI6RxRCPZDEJBMThuOUvA1Qh66Y95Soh9iULQPuCY3OIwQmFbOQq4TAUR1FHZJkcmFSSZKnl+xWxP+M1dxAiejFAP42gLuZ+VG+rqo5GYMVTdxbAQP4OSL6NSK61dd9KjPfBQD+/4f5+hsAfECNvdPX3eDvbf304pm8fwAgda0EWz/Hkn3M7Lf2Xorz1eosx277UayXKxurVThz2tXFdg+kcLys2Xkkrsboen0hqadin6R9o8obAm8Iw5G/jgnDxv9/TODjDsN1G3UdgY834cJmk6qJmIHeXXTSozsZQKfu2lwb0OnrlNGdwP2v7ql3HH136i6SS9dJH3+RD19Irh6gnpz6y9e5QDeKajF9DQTyF0tf/+J5IAx9Bx4oXINcfaxbBcTOMnVNw0vgtBwaRHPySAC/4MuTcNYSwF9h5g8S0cMAvJGI3jvSt/QVLGOl6/MJHJG5FQCOHvBA13EXNcausKJEsLM6ZiW10D72UpyvthczftQt1LZbvUxh6nRj+fjqD1JJAzKvfQQjbDjffJFsxD6gddj2OXS0rkeEFFQbXjrwnPyg+nYE6HTQTiIY4olgHjGJjYB7gGgI3CNn7PqQUthj90Qy3jkdmVxC2Qvz/29M2d4blRAzgTYcJp08X4BQkAhkagreQqvAWn/rzG8iokeY6prmZBTOVAJg5g/6/+8G8FMAHgvgQ0R0PQD4/+/23e8EcJMafiOAD/r6Gwv1pfVexMw3M/PNm/vdb81HadCgQYNRIE9Epy4ADyGi29V169TcqGtORuHMJAAiuh+Ajpnv9fd/A8C3A3g9gGcC+A7//+v8kNcD+HEi+h44I/AjAbyVmXsiupeIHgfgNgDPAPD92+xpJ8+WbWEmV70EdsrmOXcfM6WGtfcyO5WEGj/HLlBr03MXo5ALS49sJTTWMotqu0Ry3rAYiMPklKoLQkfF2YpdAN4rqFfNg5JIqPMGZFfsiMBd54PD4NxhxryEAHTX1B64c/uXZ/XPMhxrEUh5CSnNB2+8q74++MyWFdtKjMTg7LoTYsya2ADiRyStaemQegkNaqqOsSqPPN8L6B5mvnm9hetwliqgTwXwU94odQTgx5n5/yaiXwXwGiJ6FoDfB/A0AGDmdxHRawC8G06Qfbb3AAKAr0V0A30DdjAAN5XQgn0sUAnt3WVUIc6x8WcRN6C3F+q9SojVWFIInoEQKxAqPCIF4I3EFPc8eBWO0liwVQkRRaOwsleIyyhtPJbtBhcbIHsbyHkJ6TxCiZeQH6YejrjDIBZrVlcoa5WQihng6CUUllLIn3yZ1FYSyslIVEJgAm8mvIQUsXIH40DBOn+QxCrQbj/wISK6npnvMpqTUTgzAsDMvwvg8wv1Hwbw+MqY5wN4fqH+dgCPWn2PTSKYv48aAlZ7CF33KRGMYV+D6MNeKu0lbt9uYcwAtY3LaEkA0fdJ3MCIRBBdRlOJILiFemOqrN11BA6RwF4i8Fwx9YOzL3g3UUccNNvOUVoA0DGl3jRMIO6UROCQfzhvV92LNBA8Jje5hNAxgpto8rJ8f2d/MF9JnRim3xt5Jj96FVHYOA9KclgD9hsHUNOcjMJZG4HPPTSJYOE+JgiB7OVgQWQlbF2TCMYIgZ5vpkQwtjSZ/dm8QqAUx5a56BkSQeDiyxKB3IMQ1CYDAdRR4Fg7ArhjdD4ewPndDzlH2wsSBbqTIe79qANzDLaKrk+yVyUNDN5ArN5LUAvJc2/Me9igohbyHVTZPSOHj8vMELWQrBf21TEK52NuDysRACJ6JZzB9yFEdCeAb4VD/JnmZAoaAVgAZ0IMmkSwaC+LJAJft7WNoLCFGrqo0JhcIlAM+14kAmsjAICeHNL3xIQIiY2Auw5dzxgk/ULPQK/TMQypK0/v43C9VNDBI2ClEiKOyw8JQvfSgahsOKqF5Lm1Wkie25b1fI5gsJpgQiLQRHmtv3PGakEFzPz0SlNRczIGjQBsCRddPXTwiOIJiWCnvSyRCOxezPidbASFqdONqTlM/6pEoDrOkggUnhuzEYjbKIAsrUQgAB5fdz2h7zhk0KQNg045Rrb25JLJDdoAEc86wKlzCw3H+DJh4C4SnCNX59o88heXUwYwpDYCUQsBIvmYspEQiBEPld/MkwjcXlZj2sOznTdoBGAHuGzqobNOLSF7WVsiCN1mseHzJQI7NLLlce2xtBIzthI5U9GYqG+Vcrrqkm3odnJIOhIEtZiPKNYJ1iQASyajDhg2HjGeMqhzRABwqhIaOHgcYeAE2cUUCJH6dANigJUEZsnaR3Ht4cglyBtCu1MTiQ2AN1FKcHtZSSLwbev9XbMkFzpX0AjASrA1Et0FVlYPbU3QtlEP7VMiqOylqh6qse0ziUENyY9NnbX7eWzZ9rMqoKTdvlctETCcDSAgdU6ymBJFLjyoh7QLqUTcwiFZOkU0Eg8ciIKUIRlGZXOMYDNwB5MNUU+vpIHhWNaqSAQ+klcCt1iIgZQ999/pspII5FC0IBF0ADacSwQAEvfRXYGxrjixEjQCsDKcmVSwB2Kw1f7nSgULJAJgvb3sGlWcjZsiEiNTl9qLSF2tyXpttT/RXYc1xIXUqkIE53JKDLTJgMghd1EB0cAuqlgIwODbPdstxCCWOUoFftM0qM0z66OL0SXEwKWwqEkEvCHQkbIRDJ4obMLUkSiUyjb82HsNCfFyBIFkY+nH3RXOnwDQCMA+4cyJwY5r7osT36bvPvZSJC62b2HsmApo1MWUkq75eD03V7ai2rXRWIhB5s1SUxHpvQ7kCX7UH+nUFUSOGAQE74lBlAjISwGmHPozeOAoYRiCQAOAE3ffDR1YzR0kgl6XHREAvDpoiMRJ1EOB1gw+bkBJBAMjEDcnDUSJIFEPdbRqHFizAVxhuOgqIssI7U1FNIN4rU0MZhuNpySKsXlN25iraDJ3xMlZ2RqNE7dP4fgrBCC5l3mEqxZiEDU0gesHEIhBMBoPpn1wgVe5jUCpVjpO3nW47wcQUyYR0FFcixVCp40L9BoEk3n1EGsVkJIAMDj1UCh7o3EgEIOTMuS5Y5bRFaARgAYNGjS4gsBQwQznBxoBODCcuVoIWM1GEKbeRhq4yDaCjFU3eyntszBVrU+pn3D/ZMqJjUBz/AUbQdFrSN8n6iWKKiDewkag3EhpIBcbsIkqH8qMxOrphygddAMSjp6OgGGIqR1oA/ARVJmdmke4+I3vIyqiLp2PpSyqHqMu4tUCwbh5ATVIYSe1yi5QQWDbwlbEYFsbQaX/mdkICnVjRCSZqkRMRqbWexKEnfj9q3JiEJaJrI0A8V6ri2hIbQTMqQ0gqoSkP3nEyaGs3UhpcEg4IOmeQRtlNB44qJ+I042TdykVlRAP5Iy04ua5YXBPMa3FkVPhCLHhjXMlDemevUpo8P15kxKAhBgoO8gq0FRADcbgzKWDFSWDrW0Ec/YxQTzWtldU5yvte4oYFAzBpaVLfUbbKUXaAdlHu2s6idG/pzEDkSiEvgyE84hZiIDo+BGPsIQQh2j0dcRAIfwNBTsBACcdSBv7fkI8QgSYLw/snHN8/64nLwEgrMUbRwgATwx6jkdMbihIBYBD+EMfbQJVYrArNBVQgyVwJtLBipLBwY3GI/1XiW9Qcy82Gle2RqaxaDQm067mtl5GY0bjQCBKz1+SDiZUROAc4Wtfe1JupuRyQSgJIKqJgFw6AFNOEAQGRxSE42dPDMJejjyCVxy+SAVurQkVkTvf0t1vImHYHTguco6gEYALAmdOEFZUEwFbqorm7ONQ0oGfvzhfjVUfkw4K3ceWtn1qNgNtI0gIhFEBJZNw/r+WCMK9tHdONRPaxUtIq1KUqyZ5PXuIBdAIf0ilg4wg6BTZMjdHSYEZYK9icntxSewkMrg7csi/8yoih+RJSQgR6dMAsDrzYGdoKqAGa8HBCcJ5sBss3cdC6WDrvWhkOqUu2kY6qA1TCL1oFB5pt26giXRQ6JtEKZcIQofEhsBdJAiOAETDrosBoKCGIU5VOKQifzOCwPCH08jcxp10YL+eb/ZGXolC5p7QdQhuo92GMGxyFREALylgHWgqoAb7hItMENbmynfpu9VeRubN5ssWLDfUCAIbQrIzQdALakRqCIIKC5gmCHKW8aDaZxAEwEsHhiCkEkLs48rk7AFasuG4eSEgQZffO6lAYhi4Y0cEAkFA6DtsHIFYDZoXUINDwZkThLCR5VOVou/3ojI6FEFQc08ShArMJQhzJYAx+0GYbwFBsP+LVCDlxMuIU5URlArHEQNEjv8o1kk79aQIABuJAc74HNJQmL31Lgp5oyWEHuhOY1lLB7xaIJh+gecHrjYBKIjwlxV2Qqq7wEqSwmo2hLE9zOi79j5mexiNbyuocbL2MYJgOhftAtJuylYC0M9CdsyUyoihIn8FibMpm3ZtJO7jAzmCoPtLWe1/4EAgxJ015qWLB98MG14vDoDRJIBzCytyrxcJanmu9koYliDkEdhp7ytKCatKKyNbybpajt0WLUGwen0YYqOJEZvxkcEO7QzznBrBSx8Ow119gvDjhsncO2KRIvzYzqmEYNs9srcG5lTdpPZu3qFTPbnKzQDwmn8MTQK4YLASsrpocHBpYcX3vNXelzAAYwxhRc0zex8j88/aiiEKmVqnRBRk8gJhyKSAiXY27SXVkIxLBBAjHVgEnRIISqSLIC2EMqX9g0uqkih6rX4yEoI+x4BX/t03AtCgQYMGVxCYwX0/3e/A0AjAUpjBAV5GGEuLvhfpYMX3vPXel6oGJySZnd7hTHVRUSpQnUu6/8ROoPsUuP5QXZAOrAooNFtO2qqQOGWOCQV1kRprVUZQQWd5u6SeHlcZhbKVDtb8bTc30EsOU/ai8/f9V4ErRxyWrD+z71ZqoxJRqDk2UGGuKeJQIhaFcoL49U1FhaRVQXq8tGniwGYuIQ7V9kyFRIlLamyPD52qgBCIiVM3rfgDbiqgKw5XkECMIVhgDwRixXe81d6XrD+j75rSS5H7V32LaS1K+1hAHGTeTGpQ5UQakPVq0ger+WSwkR4yI68pi9QQ5s4kCrlX+Yp2BebmBdRgAub81i4ZkZhCssDKRGLFd3zWBGLONpdIMVMEojZnkbE1KqNsHwbBc6Wt5HFUMk5biaBmiCZdhykCQOv+vTUJoMHOsIQhOX+/t61gDpEQWIVYrPSOt9r3SgxntocZ76UqAdTmDBUTc0mXyh4sXrTSgB5ac9+0EkNNHVUar5F+SZpYC7hJAA0OCtsikwtMOJYgXWAFgrH0HdeQ4BbfatGQGYRgDmJPphwhWpPPU5k7kwZM3xIhGZNIqm21PrX7XYEZ6BsBaHARYA0u9IIQkW0Qr4VFRGTX9TQy3GauUiDYjL6lPWyL5GXu0vqzn8kSBZl2juQxQUD29ttt6aD3B0T0RAD/Cu7Mn3/DzN9xxlu62rCyKiODc0Rg1iAiNah67qwA1X1vSxh0F57uV1xm7netSVJ6/W2n2sNvi4FwiM15gktBAIhoA+BfA/hSAHcC+FUiej0zv/tsd9Zgb7BvArMGrPD3vk/iUoUV1Gh72/dcA/1I20FyYFlwJ9ecwcLjcCkIAIDHAriDmX8XAIjoVQCeDKARgAZnBxeBSK0AZ0KkdoCz4sObBLA/uAHAB1T5TgB/2XYiolsB3OqLn/id5/zjdx5gb4eEhwC456w3sQdoz3Vx4DI+EwB81i6D78VHf/bnh9c8ZGb3g72/y0IAZqUTY+YXAXgRABDR7cx88743dki4jM8EtOe6SHAZnwlwz7XLeGZ+4lp7WRPWOvP+rOFOADep8o0APnhGe2nQoEGDCwGXhQD8KoBHEtGnEdF1AG4B8Poz3lODBg0anGu4FCogZj4loq8D8LNwbqAvZuZ3TQx70f53dnC4jM8EtOe6SHAZnwm4pM9FfA7zUzRo0KBBg/3DZVEBNWjQoEGDhdAIQIMGDRpcUbgUBICInkhEv0lEdxDRcwvtREQv8O1vJ6IvnBpLRA8iojcS0W/7/x94qOeZ2ptq3+a5vo2I/oCI3uavJx3qecb2pdrHnunFRHQ3Eb3TjLno36r2XBfyWxHRTUT0H4joPUT0LiL6BjXmwn6riec602+1NTDzhb7gjL6/A+DPA7gOwG8A+BzT50kA3gAXL/A4ALdNjQXwfwB4rr9/LoDvvCTP9W0A/slF+1a+7b8D8IUA3mnGXNhvNfFcF/JbAbgewBf6+08G8FuX5O9q7LnO7Fvtcl0GCSCkgWDmawAkDYSGJwN4GTt4C4AHENH1E2OfDOCl/v6lAJ6y5+ewsK/nOkvY5ZnAzG8C8JHCvBf5W40911nC1s/EzHcx868DADPfC+A9cNH6MuZCfquJ57qQcBkIQCkNhP0otT5jYz+Vme8CAP//w1bc8xzY13MBwNd50fbFBxbBd3mmMbjI32oKLvS3IqJHAPgCALf5qkvxrQrPBZzdt9oaLgMBmJMGotZnVgqJM4J9PdcLAXw6gEcDuAvAd2+5v21gl2c6z7Cv57rQ34qIPgnATwD4Rmb+oxX3tgvs67nO8lttDZeBAMxJA1HrMzb2QyKi+//vXnHPc2Avz8XMH2LmnpkHAD8CJxIfCnZ5pjG4yN+qChf5WxHRMRySfAUz/6Tqc6G/Ve25zvhbbQ2XgQDMSQPxegDP8Nb9xwH4mBc/x8a+HsAz/f0zAbxu3w9iYC/PJX98Hr4cwCEzou7yTGNwkb9VFS7qtyIiAvCjAN7DzN9TGHMhv9XYc53xt9oeztoKvcYFZ7X/LTjr/rf4uq8B8DX+nuAOjPkdAO8AcPPYWF//YAC/AOC3/f8PuiTP9XLf9+1wP/TrL9AzvRJOvD6B49KedUm+Ve25LuS3AvBX4VQmbwfwNn896aJ/q4nnOtNvte3VUkE0aNCgwRWFy6ACatCgQYMGW0AjAA0aNGhwRaERgAYNGjS4otAIQIMGDRpcUWgEoEGDBg2uKDQC0KBBgwZXFBoBaHClgYieSkRvIaLfIKJfIaKHnvWeGjQ4FLQ4gAZXGojowcz8YX//rQDuYeZ/fcbbatDgINAIQIMrDUT0TQC+AsB9APwZAN8M4MMA/hZcpsp/zcw/d3Y7bNBgf3B01hto0OCsgIieAZe0668x8x8T0ZsAvItdDvh/61P6fheARgAaXEpoNoAGVxk+D8B/9Mj/vwfw38LlcxH4Z3A5YRo0uJTQCECDqwwvBfD1RPT/APhMAL/LzB/3WSC/E8Ab2J8A1aDBZYRmA2jQwAARfT1cquJfBfA2Zv6hM95SgwZ7gUYAGjRo0OCKQlMBNWjQoMEVhUYAGjRo0OCKQiMADRo0aHBFoRGABg0aNLii0AhAgwYNGlxRaASgQYMGDa4oNALQoEGDBlcUGgFo0KBBgysKjQA0aNCgwRWF/x9zr4a2yfhJcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the (square root of) the cost function \n",
    "# as a function of fit paramters\n",
    "alla1 = np.linspace(0,3e4,200) \n",
    "alla2 = np.linspace(0,0.025,200)\n",
    "\n",
    "# pre-allocate error matrix\n",
    "rms = np.nan*np.ones([len(alla1),len(alla2)])\n",
    "\n",
    "# nested for loop\n",
    "# go through each combination of parameter pairs in that range\n",
    "# 200 x 200 = 40000 values\n",
    "for ii,a1 in enumerate(alla1):\n",
    "    for jj,a2 in enumerate(alla2):\n",
    "        popfit = a1*np.exp(a2*t)\n",
    "        rms[ii,jj] = np.sqrt(np.mean((popfit-pop)**2))\n",
    "\n",
    "plt.figure()\n",
    "plt.pcolor(alla2,alla1,np.log(rms),cmap='viridis')\n",
    "#plt.plot(slope,np.exp(intercept),'ro')\n",
    "plt.colorbar()\n",
    "plt.title('log(RMS error): $a_1\\exp(a_2*t)$')\n",
    "plt.xlabel('$a_2$')\n",
    "plt.ylabel('$a_1$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-linear optimization\n",
    "\n",
    "The first step in this process is to define the function that you want to fit to data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exponential_growth(t,a1,a2):\n",
    "    '''inputs:\n",
    "    t - vector of times\n",
    "    a1 - initial population at time 0\n",
    "    a2 - specific growth rate\n",
    "    \n",
    "    returns: modeled population based on exponential growth\n",
    "    '''\n",
    "    # insert code here\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function exponential_growth in module __main__:\n",
      "\n",
      "exponential_growth(x, a, b)\n",
      "    computes f = a*exp(b*x)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(exponential_growth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10.         27.18281828]\n"
     ]
    }
   ],
   "source": [
    "fvalues = exponential_growth(np.array([0,1]),10,1)\n",
    "print(fvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scipy.optimize` library can be used to find the unknown coefficients. The `curve_fit` function is especially good at find the minima of cost functions. It needs to have 1) the function being fit (`exponential_growth1`) and 2) the data values (`t`,`pop`). It also really helps to have a good initial guess for the parameter values (`p0`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "popt,pcov = curve_fit(exponential_growth,t,pop,p0=[]) # insert first guess for p0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first variable in the output (`popt`) is the two coefficients that give the optimal fit to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.17678552e+04, 1.45465803e-02])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimating error of coefficients\n",
    "\n",
    "The second variable in the output (`pcov`) is the error covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.59526049e+06, -7.10599907e-01],\n",
       "       [-7.10599907e-01,  3.26096426e-07]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diagonal elements of the covariance matrix can be thought of as the variance that there would be in the coefficient if you resampled many times and repeated the fit procedure. This could be estimated from the data using subsets of the data.\n",
    "\n",
    "The square root of the diagonals are called the _standard errors_ of the parameter estimates (e.g. $\\sigma_a$ for parameter estimate $\\hat{a}$). If you resampled many times, and the errors are normally distributed, then about 68% of the parameter estimates would fall in the interval $\\hat{a}-\\sigma_a < a < \\hat{a}+\\sigma_a$. Confidence intervals can be computed using the standard error and the t-statistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.26303622e+03 5.71048533e-04]\n"
     ]
    }
   ],
   "source": [
    "err = np.sqrt(np.diag(pcov))\n",
    "print(err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diagonals of the error covariance matrix decribe correlations that occur in estimates of the different parameters. For example, if you resample many times, does parameter $a$ usually increase when parameter $b$ increases? If so, this means that the parameters covary when resampling, and this would result in a positive error covariance term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.pcolor(alla2,alla1,np.log(rms),cmap='viridis')\n",
    "plt.colorbar()\n",
    "#plt.plot(slope,np.exp(intercept),'ro')\n",
    "#plt.plot(popt[1],popt[0],'rs')\n",
    "plt.title('log(RMS error): $a_1\\exp(a_2*t)$')\n",
    "plt.xlabel('$a_2$')\n",
    "plt.ylabel('$a_1$')\n",
    "\n",
    "plt.plot([popt[1]-err[1],popt[1]+err[1]],\n",
    "         [popt[0],popt[0]],'r-')\n",
    "plt.plot([popt[1],popt[1]],\n",
    "         [popt[0]-err[0],popt[0]+err[0]],'r-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![images/cost-function-with-errors.png](images/cost-function-with-errors.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the parameter estimate on the map of the cost function shows that the nonlinear optimization procedure (square) succeeds in finding the minimum of the cost function, thus minimizing the sum of squared differences between the data values and corresponding model values.\n",
    "\n",
    "The standard errors can be visualized in this map of the cost function. The height and width of the red error bars are determined by the standard errors of the coefficients. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Plot the parameters obtained from the linear regression of the log-transformed data on the cost function map.\n",
    "* Use the optimization techniques described above to fit an exponential curve to the data from phytoplankton growth experiments\n",
    "* Write a function that can be used to fit an exponential curve to any data set, and return the standard error of the coefficients."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}