{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ChiSquare Test:\n",
    "This Python macro illustrates the use of ChiSquare as a goodness-of-fit measure,\n",
    "how the ChiSquare distribution comes about, how to use the ChiSquare value (along with\n",
    "number of degrees of freedom, Ndof) as a test, and that it actually works!\n",
    "At the same time, it is also an illustration of how to do a linear fit analytically,\n",
    "that is without running a minimisation algorithm (in our case called Minuit).\n",
    "\n",
    "The program generates a certain number of datasets, each consisting of 9 points along\n",
    "a line. These are then fitted (both analytically and numerically), and the result and\n",
    "the Chi2 of the fit is recorded along with the probability of the fit.\n",
    "\n",
    "## References:\n",
    "- Barlow: Chapter 6\n",
    "- Cowan: Chapter 2.7, Chapter 7\n",
    "- Bevington: Chapter 6\n",
    "\n",
    "## Author(s), contact(s), and dates:\n",
    "- Author: Troels C. Petersen (NBI)\n",
    "- Email:  petersen@nbi.dk\n",
    "- Date:   8th of November 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "from probfit import BinnedLH, Chi2Regression, Extended # Helper tool for fitting\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../External_Functions')\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initial Question/Task:\n",
    "Make sure you've read the relevant references and that you understand not only what the ChiSquare is, but also that it follows the ChiSquare distribution, and that the probability of obtaining such a ChiSquare or worse (given that the hypothesis/PDF is the correct one!) can be calculated from it.\n",
    "\n",
    "Executive Summary: The ChiSquare method is used both for getting the best fit, and for testing the quality/probability of this fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "r = np.random                         # Random generator\n",
    "r.seed(1)                            # Set a random seed (but a fixed one - more on that later.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nexp = 1000\n",
    "Npoints = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters used in this problem (constant, linear coef., quadratic coef., and uncertainty in y) for generating a line (parabola) with N points and given uncertainty on points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha0 = 3.6\n",
    "alpha1 = 0.09\n",
    "alpha2 = 0.01\n",
    "sigmay = 0.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating and fitting data:\n",
    "In the following, we generate data points (a simple line) and fit these both **analytically** and **numerically**.\n",
    "The linear fit (and only this) can be done analytically, as discussed in Barlow chapter 6.2 and outlined here: http://www.nbi.dk/~petersen/Teaching/Stat2018/StraightLineFit.pdf .\n",
    "The numerical fit of the line (and any other function) is done iteratively by Minuit. The code is slightly shorter, but a lot slower, as the code will typically test many (50-5000) possible combination of fit parameters!\n",
    "\n",
    "The whole thing is contained in a loop, allowing us to repeat the generation and fitting, so that we can see the result on more than just a case-by-case. For this reason, we save the result of the fit (in pre-defined arrays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Fit: a0= 3.635+-0.218  a1=0.078+-0.039  p=0.0500\n",
      "  Fit: a0= 3.615+-0.218  a1=0.070+-0.039  p=0.2537\n",
      "  Fit: a0= 3.700+-0.218  a1=0.085+-0.039  p=0.6989\n",
      "  Fit: a0= 3.474+-0.218  a1=0.089+-0.039  p=0.9727\n",
      "  Fit: a0= 3.581+-0.218  a1=0.103+-0.039  p=0.2892\n",
      "  Fit: a0= 3.818+-0.218  a1=0.052+-0.039  p=0.5569\n",
      "  Fit: a0= 3.703+-0.218  a1=0.098+-0.039  p=0.8759\n",
      "  Fit: a0= 3.859+-0.218  a1=0.049+-0.039  p=0.1152\n",
      "  Fit: a0= 3.601+-0.218  a1=0.094+-0.039  p=0.4949\n",
      "  Fit: a0= 3.551+-0.218  a1=0.111+-0.039  p=0.9796\n"
     ]
    }
   ],
   "source": [
    "# Arrays for storing fit results:\n",
    "array_alpha0 = np.zeros(Nexp)\n",
    "array_alpha1 = np.zeros(Nexp)\n",
    "array_Chi2 = np.zeros(Nexp)\n",
    "array_Prob = np.zeros(Nexp)\n",
    "\n",
    "# Loop, repeating the data generation and fit:\n",
    "for iexp in range(Nexp) : \n",
    "\n",
    "    # Generating data by filling values into (x,y) and associated uncertainties (here chosen to be 0 for x):\n",
    "    x = np.arange(Npoints)+1\n",
    "    ex = np.zeros_like(x)\n",
    "    y = alpha0 + alpha1*x + r.normal(0, sigmay, Npoints)   # + alpha2*x**2\n",
    "    ey = sigmay*np.ones_like(x)\n",
    "\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Numerical fit:\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Define a fit function:\n",
    "    def fit_function(x, alpha0, alpha1):\n",
    "        return alpha0 + alpha1*x\n",
    "\n",
    "    # Now we define a ChiSquare to be minimised (using probfit), where we set various settings and starting parameters:\n",
    "    chi2_object = Chi2Regression(fit_function, x, y, ey) \n",
    "    minuit = Minuit(chi2_object, pedantic=False, alpha0=3.0, alpha1=0.0, print_level=0)  \n",
    "    minuit.migrad();                                                  # Perform the actual fit\n",
    "    minuit_output = [minuit.get_fmin(), minuit.get_param_states()]    # Save the output parameters in case needed\n",
    "    \n",
    "    alpha0_fit = minuit.values['alpha0']\n",
    "    alpha1_fit = minuit.values['alpha1']\n",
    "    sigma_alpha0_fit = minuit.errors['alpha0']\n",
    "    sigma_alpha1_fit = minuit.errors['alpha1']\n",
    "    \n",
    "    Nvar = 2                     # Number of variables (alpha0 and alpha1)\n",
    "    Ndof_fit = Npoints - Nvar    # Number of degrees of freedom = Number of data points - Number of variables\n",
    "    \n",
    "    # In Minuit, you can just ask the fit function for the Chi2:\n",
    "    Chi2_fit = minuit.fval                          # The chi2 value\n",
    "    Prob_fit = stats.chi2.sf(Chi2_fit, Ndof_fit)    # The chi2 probability given N degrees of freedom\n",
    "    \n",
    "    # Fill the arrays with fit results (to produce plots of these at the end):\n",
    "    array_alpha0[iexp] = alpha0_fit\n",
    "    array_alpha1[iexp] = alpha1_fit\n",
    "    array_Chi2[iexp] = Chi2_fit\n",
    "    array_Prob[iexp] = Prob_fit\n",
    "    \n",
    "    # Let us see what the fit gives for the first couple of data sets:\n",
    "    if (iexp < 10) :\n",
    "        print(f\"  Fit: a0={alpha0_fit:6.3f}+-{sigma_alpha0_fit:5.3f}  a1={alpha1_fit:5.3f}+-{sigma_alpha1_fit:5.3f}  p={Prob_fit:6.4f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit the data and plot the result:\n",
    "Below we plot the latest fit (i.e. the last one from above), so that we can actually see what is going on (and potentially wrong!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.errorbar(x, y, ey, fmt='ro', ecolor='k', elinewidth=1, capsize=2, capthick=1)\n",
    "ax.plot(x, fit_function(x, *minuit.args), '-r')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XtcVVX6P/DPOiiHyxEBAVGuGoYJXpCj5phJjt90tEzTsnTyWqSlo4Wm5aTY5FRmv7HLWKZlTVpjF6exmbIpL+MtFRC8kVdERgVFEZUQFHh+f2w5F4E4h4uHA5/363VecPbaZ+1nbzf6uHj2WkpEQEREREREGp2jAyAiIiIiakiYIBMRERERWWCCTERERERkgQkyEREREZEFJshERERERBaYIBMRERERWWCCTERERERkgQkyEREREZEFJshERERERBaaOerAfn5+Eh4e7qjDExEREVETk5KScl5E/Kvbz2EJcnh4OJKTkx11eCIiIiJqYpRSJ23ZjyUWREREREQWmCATEREREVlggkxEREREZIEJMhERERGRBSbIREREREQWmCATEREREVlggkxEREREZIEJMhERERGRBSbIREREREQWmCATEREREVlggkxEREREZIEJMhERERGRBSbIREREREQWmCATEREREVlggkxEREREZIEJMhERERGRBSbIREREBABITEyEUsr0SkxMdHRIRA6hRMQhBzYajZKcnOyQYxMREVHVlFJwVH5AVJ+UUikiYqxuv0Y5gvzHP/7R9L/f6OjoGvWRmJhY488SERERUUXO8lsKmxNkpZSLUipVKfWvStpClVKbbrTvU0oNrtsw7TNnzhxkZ2cjISHBkWHg+vXrmDp1Knx9feHl5YUJEybgl19+qdNjJCYmomPHjvDw8EDr1q0xZswYZGdn29XHwYMHMXDgQHh7e8PPzw8jRozA//73P6tjWN7MSimMHDnSqo9169YhJiYGHh4eCA4OxrPPPotr167VyTlWZc6cObjrrrvqtM9Lly5hzJgxaNGiBfz8/JCQkIDS0lKbP79v3z4MHz4cbdq0gaenJ2JjY7F27VqrfbZu3YqhQ4eidevWUEph8+bNFfp55ZVX0L59e7i7u6Njx45YsWJFbU+NiIjI4RITE02/nRAR50+QAUwH8HMVbX8E8LmIxAB4BMDS2gZWGwaDAYGBgTAYDI4MA/PmzcMXX3yBL7/8EuvXr8fmzZsxY8aMOj1G+/btsXTpUhw8eBDr16/HmTNn8NBDD9nVx9ChQ+Hh4YGdO3fixx9/RHZ2Nh577DGrfSIjI5GdnW16WSZsx44dw8iRIzF8+HAcPHgQH330ET755BP8+c9/rpNzrMylS5fw7rvvYs6cOXXa7+TJk5GSkoINGzZg1apVWLlyJV599VWbP5+Wlobw8HB8/vnn2L9/Px555BE89NBD+O9//2va5/Lly4iOjsaiRYsq7eNvf/sbXnrpJSxevBg///wznnnmGcTHx2PTpk21Pj8iIiKygYhU+wIQDGADgP4A/lVJ+zIAs2983xvAjur6jI2NlZr67rvvpG/fvuLt7S3u7u7Sv39/SUtLq7Df/PnzJSoqqtI+AMj8+fOla9euotfrZdCgQXL27NkKn12yZIkEBASIv7+/LFmyxNReUlIikyZNknbt2omrq6uEhITI/PnzpbS01NTu6+srb775pukzq1atEr1eLwUFBTU+9+qsW7dOlFJy9epVm/bPzc0VALJ582bTtpUrV4rBYDC9/7XrKCLyxRdfiIuLi+ncRUTGjRsn9913n93xnzhxQrTb8te9+uqrEh0dLWVlZXYfoyrnz58XnU4n//znP03bXn75ZQkODq5Vv126dJEZM2ZU2F5+7Tdt2mS1/emnn5a4uDirbWFhYbJ48eJaxUFEZCtb/h4mqg1H3WMAksWG3NfWEeQlAJ4DUFZFeyKA3yulTgH4FsC0ynZSSsUrpZKVUsm5ubk2HrqinJwcjBs3Dtu2bUNqairatm2LIUOG4Pr163b1884772DhwoXYtWsXzp49i2nTrMPOyMjA0aNHsXXrVkyZMgUJCQk4deoUAKC0tBQuLi5YuXIlDh06hHfffRdLlizBsmXLTJ/Ny8tD3759Tf3dfffdKC4uxv79+2t87r8mLy8Pf/vb3xAVFQU3NzebPtOqVStERUXhyy+/xLVr11BQUICvv/4aQ4YMsdrvxIkTCAoKQvv27fHEE0/gwoULprYePXqgefPmWLNmDUQEWVlZ2LJlS4U+6kpxcTGWLFmC5557DkqpOut3z549KCsrq/BndurUKeTk5NSoTxHBpUuX4OPjY/Nn7r77buzduxcHDx4EAGzatAnnz5/HgAEDahQDERER2am6DBrAfQCW3vg+DpWPID8LIEHMI8jpAHS/1m9tRpBvlpGRIQAkJSXFant1I8gzZ840vf/3v/8tLi4ucvHiRdNnPT09paioSEREioqKRCkl69atqzKOsWPHyv333y8iItu3bxcAkpWVJSNHjpSePXtKYWGhAPjVPmrim2++EU9PTwEgvXr1kuzsbLs+f/r0aenZs6fodDpRSsm9994rhYWFpvbvvvtO1q5dK/v27ZN//OMf0qFDB+nTp4/V6O22bdskMDBQmjVrJgBk7ty5NToXW0aQ33vvPQkLC5Pr16/X6BhVWb16tbi4uEhZWZkYjUZ5+OGHJT09XQDIvn37atTne++9Jy1btpTTp09XaKtqBFlEZOnSpeLq6irNmjUTT09P+frrr2t0fCKimqju72Gi2nLUPYY6HEHuA2CoUioTwN8B9FdKrbppn0kAPr+RcP8EwA2An93Zuo2OHz+O0aNHo3379mjRogU6d+4MACgoKLCrH8tZKqKiolBaWorMzEzTtuDgYOj1egCAXq+Hh4cH8vLyTO3vvfcejEYj/P39YTAY8Nlnn1UaQ5s2bRAaGmpXbPa45557kJqaih9//BEuLi6YOnWqzZ8VEUydOhVBQUHYuXMntmzZgoKCAowfP960z6BBgzB8+HB07twZw4YNw5o1a7B9+3aUT9N39uxZTJo0CZMnT0ZycjLWrl2LlStXYulS20rRs7KyYDAYYDAYEBUVBQCm9waDAVlZWaZ9y8rKsHjxYiQkJKBZs2YV+po8ebLVZ+1V/hBiSEgI2rZta/fnLW3ZsgUzZ87Ep59+aldfW7duxcKFC/HBBx9gz549eOWVVzBu3DgkJSXVKh4iIiKykS1ZdPkLVY8gfwdg/I3v7wBwBjfmWK7qVZsR5DvuuEPi4uJk06ZNcuTIEdm2bVulI3HVjSB/9NFHpveZmZkCQPbs2VPlZz09PWXlypUiIvL3v/9dXF1d5e2335Z9+/bJ0aNH5eGHH5Z+/fqJiMiRI0es+hMRycrKEgDy008/1fjcq1N+jF27dtm0/6ZNm0QpJZcuXTJt27dvnwCQI0eOVPqZ0tJScXFxkc8//1xEtGt185/nW2+9JW3btrUphuvXr8vRo0fl6NGjsnnzZgFgen/06FGrkeI1a9aIn5+f/PLLL5X2dfbsWavP2uM///mPAJC8vDzTti1btggAu0fld+3aJS1btpRVq1ZVuU9VI8j9+vWThIQEq20PPvigjB492q4YiIhqChxBpnrmqHsMNo4gVxyCs5FS6qUbB1kHIAHAcqXUMwDkRrJcLzOMX7hwAT///DOWLl2KuLg4AMDhw4cr3ddgMODq1atV9lVe41n+vU6nQ3h4uE1xbNu2Db1797Yarc3IyICnpycAbXYJHx8fbN26FTExMQC0EUW9Xm8a8a4PLi4uAFBhOrmioiLk5OTA29sb3t7epu35+fkAAJ3O/MuE8pHZqq7dsWPHUFpaarpW+fn5Vp8v7+PXrv3N+0ZERFgdu/z9zV577TVMmzYNHh4elbYHBAQgICDApuPeLCYmBjqdzjQNG6D9mQUFBSEwMNBq3/z8fOTn5yMwMLBCvXdqaip+97vfYfHixRgzZozdcdT2ehIREVEt2ZJF18erpiPIpaWl4ufnJ2PHjpVjx47Jjz/+KF27dq10JG7jxo2ilJJvvvlGsrOzrUYdAUirVq3k22+/lbS0NImJiZGRI0ea2qsbQX7rrbfEy8tLfvjhBzl8+LDMnj1bDAaDaQRZRGT27NkSEBAgGzZskO3bt0t4eLg8/vjjNTrvyly5ckWmTZsmmzZtkhMnTsjOnTvl//7v/yQ4OFiuXLlite+mTZtMM3dYysnJkZYtW8ro0aMlPT1d0tLS5N5775WIiAi5du2aiIhMnz5dNmzYIBkZGfLf//5XjEaj9O7d2zRrxXfffScA5PXXX5fjx4/Lxo0bpV27djJ27Fi7z+nXapD/85//iKenp1y4cMHufm01atQoiYyMlF27dsn69evFx8dHXn755Qr7zZ8/v9J7bv/+/dKqVSt58cUXJTs72/SyjPnKlSuSmppq+jNZvny5pKammvaZPXu2+Pj4yD/+8Q/JyMiQTz75RFxdXeXDDz+st/MmIrJU1d/DRHXFUfcYbBxBdroEWURLfKOjo0Wv10u3bt1k3bp1VT7sNGvWLPHx8REA8pe//MXyAsmLL74oUVFRotfrZeDAgZKTk2Nqry5BLi4ulkmTJom3t7d4e3vLtGnT5Mknn7RKkIuLi+Xpp58Wb29vadGihYwbN65Op3i7evWqPPjggxIUFCSurq4SEBAgw4YNk/T09Ar7VpUgi4js2LFD+vXrJ15eXuLr6ytDhw61Kq8YMWKE+Pv7S/PmzSU4OFgmTZok586ds+pj5cqVEh0dLe7u7hIUFCRPP/20XL58uc7OVUSkf//+lU6XVpfy8/Pl0UcfFU9PT/H19ZVnn31WSkpKKuxXVYJcvv3ml+V9Uf5ncfPL8t6aM2eOhIWFiZubm9x+++1WUwwSEdU3JshU3xp6gqykfiohqmU0GqX8IS9HUErhiy++qLAiHDVMSUlJ6NOnDzIyMhAcHOzocIiIGjWlFByVH1DT4Kh7TCmVIiLG6vazZyU9IofJy8vDihUrmBwTERFRvavxQ3pEt9LAgQMdHQIRERE1EU02QeavjoiIiIioMiyxICIiIiKy0CgT5D/+8Y+mFdEsV8uzR2JiYo0/S0RERETOq1EmyHPmzEF2djYSEhIcGkf5ghOtW7eGUgqbN2+ul+MsWLAAgYGB8PT0xLBhw3Du3Dm7Pn/p0iWMGTMGLVq0gJ+fHxISElBaWmq1z44dO9C3b194eHigVatWVktRA9riFvHx8QgICIC7uzu6deuGI0eO1PbUftWcOXNw11131WmftlyL6qSnp6Nv375wd3dH+/bt8fHHH1u1T548Gbfddhvc3NwQFBSEp556CpcvXza1X7t2DS+88ALCwsLg5uaG2NhYbNmypU7Oj4iIiKrXKBNkg8GAwMBAGAwGh8Zx+fJlREdHY9GiRfV2jGXLlmHRokVYunQptm7ditOnT9u9etvkyZORkpKCDRs2YNWqVVi5ciVeffVVU/vBgwcxYMAAxMTEYOfOndixYweGDx9uahcRDBs2DLt27cJnn32G9PR0vP7666ZVBevDpUuX8O6772LOnDl12m9116I6xcXFGDx4MAICApCUlIQZM2Zg4sSJ2LZtm2mfrl274uOPP8ahQ4fwxRdfYOvWrXjqqadM7YsWLcLy5cuxbNkyHDhwAP369cPgwYNx5syZOj1XIiIiqoItkyXXx6s2C4V899130rdvX/H29hZ3d3fp37+/pKWlVdivssU+yuHGohldu3YVvV4vgwYNkrNnz1b47JIlSyQgIED8/f2tFmsoKSmRSZMmSbt27cTV1VVCQkJk/vz5ptXlLOXm5la5kEltde3aVZ555hnT+23btgkAOXr0qE2fP3/+vOh0OvnnP/9p2vbyyy9LcHCw6f24ceOkf//+VfaxadMmadasmWRkZNTgDKz92kp6ll599VWJjo6WsrKyWh+znC3Xojr/+Mc/xMXFxWrlvAEDBsiYMWOq/Mxbb70lbdq0Mb3v2bOnzJo1y/S+fPXI119/3eY4qHG6eSGayhb+IaoLtvw9TFQbjrrHYONCIU45gpyTk4Nx48Zh27ZtSE1NRdu2bTFkyBBcv37drn7eeecdLFy4ELt27cLZs2cxbdo0q/aMjAwcPXoUW7duxZQpU5CQkIBTp04BAEpLS+Hi4oKVK1fi0KFDePfdd7FkyRIsW7aszs6zOsXFxdi/fz/69u1r2ta7d280b94cti7CsmfPHpSVlVn1cffdd+PUqVPIyckBAGzevBl33nkn7r//fgQEBOA3v/kNNm3aZNp/8+bNiIqKwpo1axAcHIzIyEjMnz8fJSUldXSm1oqLi7FkyRI899xzUErVWb+2XIvqJCUlISoqCr6+vlZ9VPXncebMGXz++eeIjY01bbt27Rr0er3pvU6ng6urK1JTU+09JWpkEhMTTTPwiAgSExMdGxARUSPllAny+PHjMWnSJERFRSEyMhIvvfQSTp8+jf3799vVz4QJEzBkyBB07doVL7/8Mr766ivk5+eb2nU6Hd544w3cfvvteOGFF1BWVmZKUlxdXbFs2TL069cP7dq1w5AhQ/DAAw/gu+++q9Nz/TUXLlxAWVkZ/Pz8MHPmTLRv3x5FRUXw8fFBbm6uTX3k5ubCxcUF3t7e6NGjB0aNGgU/Pz9TGwBkZ2fjr3/9K+666y58//33uPPOOzF48GDTfxays7Nx4sQJ/Pjjj1i7di3+/Oc/46233sKSJUvq5bw/+ugj6PV6PProo3Xary3XwpY+/Pz88L///Q/+/v5455134OfnV+HzS5cuhYeHB4KCgtCyZUt8+umnprZ+/frh008/RUZGBkpLS7F06VKcO3fO5hiIiIiodpwyQT5+/DhGjx6N9u3bo0WLFujcuTMAoKCgwK5+LGepiIqKQmlpKTIzM03bgoODTSN5er0eHh4eyMvLM7W/9957MBqN8Pf3h8FgwGeffWZ3DHXF398foaGhaNbM/qmty2f8CAkJQdu2bSu0l5WVITY2FrNnz0ZMTAwWL14MHx8ffPnll6b2y5cvY+XKlejZsydGjBiB+Ph4rFq1yqbjZ2VlwWAwwGAwICoqCgBM7w0GA7KysqxiWbx4MRISEio918mTJ1t9tq6vha1cXV0RFhZmNZJsacyYMUhLS8M333yDjIwMLFiwwNSWmJiIqKgoREREQK/X46uvvsKgQYPqdLSciIiIquaUC4Xcf//9aN26NT788EMEBQXh3LlzuOuuu1BWVlbrvst/fQmg0gSsvH3NmjWYPn063njjDfTr1w/u7u6YO3cuzp49W+sYbNWqVSvodDqcP38es2fPxuzZs1FWVoaLFy/C39/fpj78/f1RUlKCixcvYu3atQC02TfK2wDAz88Pt99+u+kzOp0OYWFhOH36tKndzc0NISEhpn3at29vaq9O27ZtkZaWBgA4ffo04uLiTO/L28t9+eWXyM/Px6RJkyrt66WXXsLMmTNtOu7NbLkWtvSxa9cutG7d2lRW8ac//anC51u2bImWLVvi9ttvR8uWLdGvXz/MnDkTgYGB8Pb2xrp161BYWIjLly8jMDAQffr0QWRkZI3Oi4iIiOzjdCPIFy5cwM8//4z58+cjLi4OHTp0wMWLFyvd12Aw4OrVq1X2dfDgQavvdTodwsPDbYpj27Zt6N27N6ZOnYrOnTsjIiICGRkZdp1Lben1ekRHR5uSOAD46aefcP36dRiNRqt9i4qKkJmZaVVCAgAxMTHQ6XRWfWzZsgVBQUEIDAw07WN5biKCU6dOISgoyNReVFRkNctCVlaWqb06zZo1Q0REBCIiIhAWFgYApvcRERFW/1F57bXXMG3aNHh4eFTaV0BAgNVn7WHLtSiXn5+PzMxMFBUVWW03Go04ePCg1W8atmzZUuHPw5KLiwtEpMK96uHhgcDAQJw8eRI7d+5EXFycXedDRERENWTLk3z18arpLBblT/SPHTtWjh07Jj/++KN07dq10lkiNm7cKEop+eabbyQ7O1t++eUXy6cYpVWrVvLtt99KWlqaxMTEyMiRI03tlc2A4enpKStXrhQRbeYBLy8v+eGHH+Tw4cMye/ZsMRgM0q9fP9P+V65ckdTUVNm0aZMAkOXLl0tqaqrVDAe19e6774qHh4d89dVXkpKSIkajUQYMGFBhv/IYKnvqfdSoURIZGSm7du2S9evXi4+Pj7z88sum9rVr14pSSpYuXSpHjx6VefPmiZubm5w8eVJERAoKCqRNmzZy//33S3p6unz//ffi5+cnixYtsvt8fm0Wi//85z/i6elZp9fvZtVdi3LlswncfM8VFRVJWFiYPPjgg3LgwAF56623RKfTyZYtW0RE5NixY/L888/L9u3bJTMzUzZt2iTdu3eXbt26mfo4e/asfPzxx3Ls2DH56aefpEePHnLHHXdIUVFRvZ03OZeqfkaI6grvMapvjrrHYOMsFk6XIItoiW90dLTo9Xrp1q2brFu3rspp1GbNmiU+Pj4CQP7yl79YXiB58cUXJSoqSvR6vQwcOFBycnJM7dUlyMXFxTJp0iTx9vYWb29vmTZtmjz55JNWCXJ5Unrzq7yPujJ//nwJCAgQd3d3eeCBB6zO4+ZYKkuQ8/Pz5dFHHxVPT0/x9fWVZ599VkpKSqz2eeuttyQsLEzc3d2le/fu8sMPP1i1p6WlyV133SV6vV5CQkJk7ty5cv369To9z/79+8uMGTPqtM+b2XItRKpOkEVEDhw4IH369BG9Xi/h4eFWf96nT5+WwYMHS2BgoLi6ukrbtm1l7NixcurUKdM+Z8+ele7du4ubm5t4e3vLI488IqdPn66P0yUnxeSF6hvvMapvDT1BVmJRc3srGY1GsXUqsvqglMIXX3yBkSNHOiwGsl1SUhL69OmDjIwMBAcHOzocIodSSsFRf3dT08B7jOqbo+4xpVSKiFRd93iD09UgU9OUl5eHFStWMDkmIiKieueUs1hQ0zNw4EBHh0BERERNRJNNkPmrIyIiIiKqDEssiIiIiIgsMEEmIiIiIrLQKBPkzZs3QymF8+fPV7tvYmKi1ZLTRERERNS0OW2CfOzYMTz44IPw9vaGwWBAnz59sGfPHrv7mTlzJv773/9W2L58+XL06NEDXl5e8Pf3x0MPPYTMzMw6iJycQflqhLb+R8vSK6+8gt69e8PNza3SlRmPHDmCIUOGwMfHB97e3vj9739vtRrktWvX8MILLyAsLAxubm6IjY3Fli1bantKREREZCOnTJCzs7Pxm9/8BhcvXsS//vUvJCUl4ZFHHkF2drbdfRkMBrRq1arC9m3btuHxxx/Hjh078OOPPyIvLw+DBg1CSUlJXZwCNXALFiyAp6dnjT5bWFiIhx9+GGPHjq3QVlJSgqFDh0Kn02H79u34/vvvsX//fkycONG0z6JFi7B8+XIsW7YMBw4cQL9+/TB48GCrpbyJiIioHtmymkh9vGqzkt60adPEy8tLCgoKKm0vXzVuzZo1EhkZKS1atJBx48ZZrez2yiuvmFa2u3nFvMrs2bNHAEhaWlqN4ybnsH37dunSpYt8//33AkByc3Nr1M/rr78uYWFhVtvS09MFgBw8eNC0bd26daKUkrNnz4qISM+ePWXWrFmm9vLl1V9//fUaxUGND7jKGdUz3mNU3xx1j8HGlfSccgR5/fr1+N3vflftCN+KFSvw5ZdfYtWqVVi1ahW+/PJLU9u0adOQnZ2NhIQEm46Zn58PAPD29q554NTgFRQUYPz48Xj//ffh6upa5/1fu3YNAKDX603b3NzcICLYu3evaR/Ldp1OB1dXV6SmptZ5PERERFSRUybIWVlZNq2oNm/ePERHR2Po0KHo2bMnkpKSTG2enp4IDAyEwWCotp/S0lLMmzcPI0eORFhYWK1ip4Zt+vTpGDp0KHr16lUv/UdGRqJ169ZYvHgxioqKkJeXhzfeeANKKeTm5gIA+vXrh08//RQZGRkoLS3F0qVLce7cOVM7ERER1S+nTJBtFRERYfre19cXeXl5Nepn2rRpuHDhApYvX15XoVEDtG7dOmzduhV/+tOf6u0Ybm5u+Oyzz7B+/Xp4enoiKCgIcXFxEBEopQBoM6tERUUhIiICer0eX331FQYNGmRqJyIiovrllAlyaGgoTp8+Xe1+zZpZLxQoNVg9b/bs2Vi/fj1++OEHllc0chs3bkRGRgZ8fHzg5uaGe++9FwAQHByMJUuW1Nlx7rnnHpw4cQLZ2dnIzc3FiBEjAABBQUEAtDKedevWoaCgAKdOncKGDRuQl5dnaiciIqL65ZQJ8sCBA/Hdd9+hsLCwXo+TmJiIzz77DBs3bmRy0gS88MILOHDgANLS0pCWloYVK1YA0ObVvnlGivz8fGRmZqKoqKjGxwsICIDBYMBnn30Gg8EAo9Fo1e7h4YHAwECcPHkSO3fuRFxcXI2PRURERLZzygR5zpw5aN68Oe6//37s2LEDhw4dwl//+lf8+9//trmPnJwc5OTkoKCgACUlJab35Q9Rvfrqq3jzzTfx6aefws3NzdR+9erV+jotcrCAgAB07NjR9AoNDQWgler4+vpa7btkyRK0a9cOO3furNBPVlYW0tLScObMGVy7ds2UcJdbv349tmzZgpMnT+Kjjz7Ca6+9hlmzZsHDwwMAcO7cOfztb3/D8ePHsXPnTjz00EOIjIzEqFGj6vHsiYiIqFyz6ndpeIKCgrBt2zbMnj0bv/vd71BSUoLOnTtj6dKlNvfRpk2bSt9v2rQJcXFxeO+995Cfn4++ffta7bdy5UqMHz++1udAjde8efPw8ccfm97HxMQAMJf4XLhwAU888QTOnj2LkJAQJCYmYubMmVZ9vPnmm3jyySfh5uaGQYMG4Y033rCa2YKIiIjqj6pJXW5dMBqNkpyc7JBjExE5M6VUjZ6pILIV7zGqb466x5RSKSJirG4/pyyxICIiIiKqL0yQiYiIiIgsMEEmIiIiIrLABJmIiIiIyEKTSZDj4uIwderUCttTU1PRrVs3NG/enCuVEREREVHTSZCrMmfOHISEhOD48ePIzs52dDhEiIuLg1Kqwuvpp592dGhERERNgs0JslLKRSmVqpT6VxXtDyul0pVSB5VSn9ZdiPXr+PGFqIq8AAAgAElEQVTjuPvuuxEaGorAwEBHh0OEtWvXIjs72/Tau3cvdDqdaUlqIiIiql/2jCBPB/BzZQ1KqQ4AngfQR0SiAMyog9iqNH78eNx3332YM2cOvL29ERISgi+++MJqn48//hghISHw9PTElClTUFZWZtUeHh4OpRSOHz+O5557zjRKR+Rovr6+CAwMNL3Wr1+PsLAw3HPPPY4OjYiIqEmwKUFWSgUDGAJgRRW7PAHgryJyEQBE5FzdhFe1TZs2oUWLFkhKSsKgQYMQHx9vWgY6PT0dEydOREJCAvbs2QMA2LFjh9Xnk5KSkJ2djeDgYMybN880WkfU0KxYsQLjx4/nf+CIiIhuEVtHkJcAeA5AWRXttwO4XSm1XSm1Uyk1qLKdlFLxSqlkpVRybm5uDcI1a9u2LebOnYsOHTpg1qxZyM/Px7FjxwAAH3zwAYxGI2bMmIHIyEj85S9/gcFgsPq8v78/AgMD4eLighYtWphG64gaks2bN+P48eNc3pyIiOgWqjZBVkrdB+CciKT8ym7NAHQAEAfgUQDLlVLeN+8kIu+LiFFEjP7+/jUMWXPbbbeZvvf19QUA5OXlAQCOHTuGTp06mdrd3NwQERFRq+MROcLy5csxYMAAhIaGOjoUIiKiJsOWEeQ+AIYqpTIB/B1Af6XUqpv2OQVgnYhcF5ETAI5AS5jrTbNmzSps47rx1Jjk5eXhq6++wsSJEx0dChERUZNSbYIsIs+LSLCIhAN4BMBGEfn9Tbt9DW30GEopP2glFxl1G6rtIiIikJ6ebnpfVFRkKr8gchaffPIJPD09MWzYMEeHQkRE1KTUeB5kpdRLSqmhN95+D+CCUiodwCYAs0TkQl0EWBMTJ05EcnIy3nzzTRw+fBgJCQkoKChwVDhENbJ8+XKMHj0aer3e0aEQERE1KXYlyCKyWUTuu/H9PBFZd+N7EZFnRaSTiHQWkb/XR7C2ioqKwgcffIDFixeje/fuKCsrw29+8xtHhkRkl59++gkHDx5keQUREZEDKEfV7RqNRklOTnbIsYmInJlSis9cUL3iPUb1zVH3mFIqRUSM1e3X5JeaJiIiIiKyxASZiIiIiMgCE2QiIiIiIgtMkImIiIiILDBBrkJiYiKio6MdHQYRERER3WJOmSCPHz8eSikopaDX6xEZGYlXXnkFpaWljg6NnNylS5cwZswYtGjRAn5+fkhISLDrvkpMTDTdm5avqKgo0z5xcXEV2t95551fbVdK4emnn67TcyUiIqLKVVyv2UkMGDAAn3zyCYqKirBu3TrMmDEDOp0Os2fPdnRo5MQmT56M1NRUbNiwAXl5eRg9ejR8fX0xd+5cmz4/c+ZMTJ482Wpbv379MGLECKttEyZMwJ///GfTey8vL9P3a9euxbVr10zvz507h5iYmAp9EBERUf1wyhFkANDr9QgMDER4eDj+8Ic/4Le//S2+/vprU3tmZiaUUvjXv/6FwYMHw8PDA4GBgdiwYQMAICsrC/fddx88PT3h4+ODCRMmVLra3nPPPQcvLy8EBQVhxYoVt+z86Na7cOECPv/8cyxatAg9e/bEoEGDkJCQgPfee8/mPgwGAwIDA02vrKwsHD16FBMmTLDar/x+LH95eHiY2nx9fa3a1q9fj7CwMNxzzz11dq5ERERUNadNkG/m7u6O69evV9g+c+ZMDBs2DPv27cNHH30ET09PAMDo0aORn5+PHTt2YN26ddiyZQtmzZpl9dlDhw4hOzsbu3fvxuzZsxEfH48DBw7ckvOhW2/Pnj0oKytD3759TdvuvvtunDp1Cjk5OTXqc/ny5YiLi0O7du2stq9ZswZ+fn7o3LkzXnnlFZSUlFTZx4oVK0xlRURE9Wr1aiA8HKUAEB6uvSdqgpw+QS4rK8O3336L9evXY8CAARXaH374YcTHxyMiIgKDBg3CnXfeiQMHDmD79u3461//iq5du6Jv375YsGABPvzwQ6tfbet0Orz99tvo2LEj/vCHP8BoNOLDDz+8ladHt1Bubi5cXFzg7e2NHj16YNSoUfDz8zO12augoAB///vfMWnSJKvtjz32GNasWYONGzfiqaeewiuvvIJ58+ZV2sfmzZtx/PhxjB8/3u7jO4Oba7YTExMdHRJR07V6NRAfD5w8qSUHJ09q75kkUxPktAny+vXrYTAY4ObmhuHDh2PcuHGV/uNqORpY7tixY9DpdFazVHTt2hXXrl1DVlaWaVtQUBC8vb1N76OionD8+PG6PRFqUMoTtZCQELRt27ZWfX322WdwcXHBgw8+aLV90qRJ6N+/P7p06YIpU6ZgwYIFeOuttypdcnP58uUYMGAAQkNDaxVLQ5WYmGg6bxFhgkzkSHPnAoWF1tsKC7XtRHXFSX5L4bQP6d199914//33odfr0bZtW7i4uFS6n4+PT50el7/mbrz8/f1RUlKCixcvYu3atQCArVu3mtrstXz5cjz66KNwd3f/1f1iY2Pxyy+/4Pz581bHycvLw1dffYWPP/7Y7mMTEdnNYoDIpu1E9ir/LUVhofVvKQBgzBhHRlaB044ge3h4ICIiAiEhIVUmx1W57bbbUFZWZlVPvHfvXjRv3txqpO706dPIz883vT948CAiIiJqHzw1SDExMdDpdKakGAC2bNmCoKAgBAYGWu2bn5+PzMxMFBUVVdrX3r17kZSUhIkTJ1Z73EOHDsHT09NUzlHuk08+gaenJ4YNG1aDsyEislNVv6lqpL/BIgdwot9SOG2CXBudO3dG7969MXXqVOzduxdbt27F/PnzMWHCBLi6upr2Kysrw/Tp03H48GG8/fbbSE5OtinhIefk5+eHhx56CM899xx2796N77//Hm+88QamTJlSYd8lS5agXbt22LlzZ6V9LV++HNHR0ejRo4fV9mPHjmHevHnYvXs3Tpw4gc8//xwvvvgipk6dWuG3E8uXL8fo0aOh1+vr7iSJiKqycCFgMaMOAO39woWOiYcaHyf6LYXTlljU1qeffoqnnnoKvXv3hqurKx544AG8/vrrVvt07NgRrVu3htFohJeXFz744AN06tTJQRHTrbBs2TJMmTIF/fv3h16vx4QJEzBnzhy7+rh69SpWr15d6YN3er0eP/74I95++20UFhYiPDwcM2bMwMyZM632++mnn3Dw4EF88skntTofIiKblf+Ke+5clJ08CV1YmJYcN7BffZMTCw3Vyioq297AqMoeDLoVjEajJCcnO+TYRNSwKKUqfUiRKsfrRfWN9xjVC4saZBMPD+D992/Zf8SUUikiYqxuvyZZYkFERERUW5yq0k5jxmjJcFgYygAgLOyWJsf24AgyETkcR6vsw+tF9Y33mH14veznqGvGEWQiIiIiohpggkxEREREZIEJMhERERGRBSbIREREREQWmCATEREREVlggkxEREREZIEJMhERERGRBSbIREREREQWmCATEREREVlggkxEREREZKGZowMgIiIioibgwgVg925g926EOjqWajBBJiJyFqtXA3PnohQAwsOBhQuBMWMcHBQRUSWKi4G0NGDXLu21ezdw7JjWphS6Oza6ajFBJiJyBqtXA/HxQGGhVht38qT2HmCSTESOJQIcPaolweUJcVoacP261t62LdCrF/D449rX2Fh87eXl2JiroUTEIQc2Go2SnJx8S4+ZmJiIBQsWmN7Pnz8fiYmJtzQGIqpIKQVH/V3kNMLDtaT4ZmFhQGbmrY6GGjn+TNqnyV2v8+fNo8LlXy9e1No8PYEePYCePbVkuFcvICioQheOumZKqRQRMVa7X1NKkMs1uRuZqIHjz6QNdDptlOZmSgFlZbc+HmrU+DNpn0Z9vYqKgNRU64Q4I0Nr0+mA6GgtCS5PiDt1Alxcqu22oSfILLEgInIGoaGVjyCHNvRHXYjIaZSVAUeOWJdK7N0LlJRo7cHBWhI8ebKWEMfGAgaDY2OuJ0yQiYicwcKFphpkEw8PbTsRUU2cO2c9MpyUBOTna20Gg1YqMXOmeYS4bVvHxnsLMUEmInIG5Q/izZ2LspMnoQsL4ywWRGS7q1eBPXusE+Ly5xdcXIDOnYGHHzbXDXfsaFOpRGPFGmQicjj+TNqH14vqG+8x+zS461VWBhw+bC6T2LUL2L/fXCoRGmpOhHv2BLp31x6uu4VYg0xERERE9Scnx7puOCkJuHxZa/Py0kolnnvOnBAHBjo2XifABJmIiIjIWRQWAikp1qUSWVlam4sL0KULMHq0eYQ4MlKbbYLsYnOCrJRyAZAM4LSI3FfFPiMAfAmgh4g4pn6CiIiIqDEoLQUOHbJejW7/fm07oM2P3rs3MH26lgzHxGgP71Kt2TOCPB3AzwAqXfpEKdXixj676iAuIiIioqblzBnrUonkZODKFa2tZUutPOL557WvPXsCrVs7Nt5GzKYEWSkVDGAIgIUAnq1itz8BeA3ArLoJjYiIiKiRKijQSiUsE+JTp7S2Zs2Arl2Bxx4z1w3ffjtLJW4hW0eQlwB4DkCLyhqVUt0BhIjIv5VSVSbISql4APEAEMrJ7YmIiKgpKC0F0tOt64YPHDCvgtm+PXDXXea64W7dAHd3x8bcxFWbICul7gNwTkRSlFJxlbTrAPw/AOOr60tE3gfwPqBN82ZvsEREREQN3qlTFUslfvlFa/Px0UaEH3jAPDrs7+/YeKkCW0aQ+wAYqpQaDMANgJdSapWI/P5GewsA0QA2K6UAIBDAOqXUUD6oR0RERI3alSuIA4DXXjMnxGfOaG3Nm2ujwRMmaIlwr15Ahw6Ali9RA2bXQiE3RpBnVjWLxY19Nt/Y51eTYy4UQkTl+DNpH14vqm+8x6pQUgIcPGhdKpGebi6VuO02c5lEr15aHbGbm2NjbqAa7UIhSqmXACSLyLqa9kFERETUIIlopRKWq9GlpGjzEAOAr682KjxiBAYvWIBvc3MBPz/Hxkx1hktNE5HD8WfSPrxeVN+a5D12+bJWK2yZEOfkaG2urtocw+VlEr16aaPFN0olmuT1qqVGO4JMRERE5JRKSrQFNywfpPv5Z23UGNDqhAcMMD9E17UroNc7Nma6pZggExERUeMloi3FbLkaXUoKcPWq1t6qlZYIjxqlfe3RQyufoCaNCTIRERE1HpcuAUlJ1gnx2bNam14PdO8OxMebSyXateOsElQBE2QiIiJyTteva6USlnXDhw6Z2yMjgYEDzaUSXbpo9cRE1WCCTERERA2fCJCZaV03vGcPUFSktfv7a4nwmDHaV6NRW5SDqAaYIBMREVHDk5+vJcOWCXFurtbm5gbExgJTpphLJcLCWCpBdYYJMhERETnWtWvAvn3WC3AcPmxuv+MOYMgQ8zRrnTtrq9QR1RMmyERERHTriAAnTljXDaemAsXFWnvr1loSPHaslhD36AG0bOnYmKnJYYJMRERE9Scvr+KsEufPa23u7lqpxNSp5lKJkBCWSpDDMUEmIiKiulFcDOzda103fPSo1qYU0KkTMHSouVQiOhpoxlSEGh7elURERGQ/EeD4ceuR4dRUrZ4YANq00ZLgiRO1hNhoBLy8HBszkY2YIBMREVH1LlywHhnevVsrnwAADw8tAZ4+3TzncHAwSyXIaTFBJiIiIiuugPVDdLt2aaPFAKDTAVFRwPDh5rrhTp1YKkGNCu9mIiKipkxEqxO2GB2+AgB33qm1BwVpSXB8vDYyHBsLtGjhyIiJ6h0TZCIioqYkN9ecDJcvxHHxotZmMABGI/4fgDlr12oJcVCQQ8MlcgQlIg45sNFolOTkZIccWykFR503EVXEn0n78HqRzYqKtAfnLEslTpzQ2nQ6bRaJ8jKJXr20BTlcXHiP2YnXy36OumZKqRQRMVa3H0eQiYiIGoOyMuDIEesH6fbuBUpKtPaQEC0JLl+eOTYW8PR0bMxEDRQTZCIiImd07pz1jBK7dwOXLmltLVpoK9DNmmWec7hNG8fGS+REmCATERE1dFevAnv2WCfEmZlam4sL0Lkz8Mgj5inWOnbUthNRjTBBJiIiakjKyoBDh6xLJfbtA0pLtfawMC0JLl+euXt3bR5iIqozTJCJiIgcKSfHPCq8axeQlARcvqy1eXlpyfDs2ebR4cBAx8ZL1AQwQSYiIrpVCguBlBTrUomsLK2tWTOgSxdgzBhz3XBkpDbbBBHdUkyQiYioUUpMTMSCBQtM7+fPn4/ExMRbF0BpqVYqYTnF2oED5lKJ8HCgd29gxgwtIe7eHXB3v3XxEVGVOA8yETkcfybtw+tln1t2vc6csS6VSE4GrlzR2lq2NI8Kl5dKBATUf0w1xHvMPrxe9uM8yERERI1NQYG5VKI8IT51Smtr1gzo1g147DFzQtyhA0sliJwIE2QiIqJfU1oKpKdbl0ocPKjNNgEA7dsDffuaR4hjYgA3N8fGTES1wgSZiIjI0qlT1lOsJScDv/yitfn4aInw8OHa1549AX9/x8ZLRHWOCTIRETVdV65oCbBlqcSZM1pb8+ZaqcSECea64Q4dAKUcGzMR1TsmyERE1DSUlGilEZZTrB08CJQ/KBQRAcTFmeuGu3UD9HqHhkxEjsEEmYiIGh8R4H//A3bvxiIAuPtu7aG6wkKt3ddXS4JHjDCPDrdq5ciIiagBYYJMRETO7/JlbQU6y9rhnBwAwB8A4No14PHHzQ/S3XYbSyWIqEpMkImIyLmUlAD791vXDf/8s7lUokMHYMAA08iwV69eKN6507ExE5FTYYJMRI6zejUwdy5KAW1VsYULtWV2icqJaEsxW9YNp6QAV69q7X5+2qjwqFFaQtyjh1Y+YeGaA8ImIufGBJmIHGP1aiA+HigshA4ATp7U3gNMkpuy/HzzrBLlCfHZs1qbXq8txxwfb36Qrl07lkoQUZ3jUtNEdSQxMRELFiwwvZ8/fz4SExMdF1BDFx6uJcU3CwsDMjNvdTROpdH8HXb9OrBvn3Xd8KFD5vbISPMDdL16AV26AK6udh+m0VyvW4jXzD68XvZr6EtNM0EmqmO8v2yk05lrRi0pZV6hjCrllPeYiPYfH8u64T17gKIird3f3zwq3LOnVirh41Mnh3bK6+VgvGb24fWyX0NPkFliQUSOERpa+QhyaOitj4Xq3sWL2qwSlglxbq7W5uamlUpMmWJOisPCWCpBRA0GE2QicoyFC001yCYeHtp2ci7XrmmlEuVlErt2AUeOmNvvuAMYPNicDHfurK1SRw2OZamYUoqlYtRkscSCqI7x/rLDjVksyk6ehC4sjLNY2Mih95gIkJFhXTecmgoUF2vtrVtb1w336AG0bOmYWG/gzyTVN95j9mvoJRa6WxEMEVGlxowBMjPhAmj1qUyOG568PGD9emDBAmDIECAgQFuSefRo4P33gWbNgKlTgTVrtD/D7Gzgn/8E5s7V5iJ2cHJMVK9WrwbCw81TVa5e7eCAGr7ExESoG+VUSqkG+xsKm0eQlVIuAJIBnBaR+25qexbA4wBKAOQCmCgilRQXmnEEmRor3l/24zWzT71dr+JiYO9e6ynWjh4tP6hWKlFeJtGrFxAV5RSlEry/qF5YTFVp4uGh/ceR/9lvsOp8FosbSbARgFclCfI9AHaJSKFSagqAOBEZ9Wv9MUGmxor3l/14zexTJ9dLBDh+3LpuOC1NqycGgMBA62TYaAS8vGofvAPw/qJ6wakqnVKdzmKhlAoGMATAQgDP3twuIpss3u4E8Hsb4yQiolvh/HltRLi8dnj3bq18AtBGvYxG4A9/MCfEwcGcVYLo12Rl2bednIqts1gsAfAcgBY27DsJwHc1joiIiGqnqEgbDbacYu34ca1NKa00Yvhw84N0UVFaLTER2Y5TVTZq1f6NqJS6D8A5EUlRSsVVs+/voZVh9KuiPR5APACE8gYiIqo9Ea1O2LJuOC1NW6UOANq21ZLgJ57QvsbGAi1sGesgol/FqSobNVuGDPoAGKqUGgzADYCXUmqViFiVUSilBgCYC6CfiBRX1pGIvA/gfUCrQa5V5ERETVFuLgYDwLx5WkKclKQtygEAnp5aqcQzz5hLJYKCHBktUeNV/iAep6pslOyaB/nGCPLMSh7SiwHwJYBBInLUlr74kB41Vry/7MdrVoWrV7U5hi3nHD5xQmvT6YDoaHOZRK9eQKdOgIuLY2NugHh/UX3jPeY86n2paaXUSwCSRWQdgNcBGAB8cWNuuywRGVrTvomImpyyMm31Ocu64b17gZISrT04WEuCJ09Gv9mz8d9LlwCDwbExExE1UlxJj6iO8f6yX5O8ZufOWdcN794NXLqktRkM2gp0livStW1r+miTvF61wOtF9Y33mPOo9xFkIiKyUWGhViphOedw+dPvLi5aqcSoUeZSiY4dWSpBRORATJCJiOpSWRlw6JB13fC+fUBpqdYeGqolwVOnal+7d9ceriMiogaDCTIRUW3k5FjXDSclAZcva21eXlqpxOzZ5lKJwEDHxktERNVigkxEZKvCQiAlxTohLl81y8UF6NIFGD3aXDvcsaM22wQRETkVJshERJUpLdVKJSwfpNu/31wqER4O9O4NTJ+uJcQxMdoiAURE5PSaVoK8ejUwdy5KAe0fN07oTUTlsrOtH6JLTgauXNHaWrbURoSff1772rMn0Lq1Y+MlIqJ603QS5NWrTUtC6gDtCfL4eK2NSTJR01JQoJVKWD5Id+qU1tasGdC1K/DYY+ZSidtvZ6kEEVET0nTmQQ4PN0+rZCksDMjMvHVxUKPH+TDtV6/XrLQUSE+3rhs+cECbbQIA2re3Xo0uJgZwc6ufWOoI7zH78HpRfeM95jw4D/LNyh+ksXU7ETmn06et64aTk7URYwDw8dGS4QceMI8O+/s7Nl4iImpwmk6CHBpa+QhyaOitj4WI6saVK+ZZJcoT4tOntbbmzYFu3YDx483JcIcOgFIODbk2EhMTsWDBAgDaiNX8+fORmJjo2KCIiBqhplNiYVGDbOLhAbz/PmuQqU7xV232s+malZQABw9a1w2np5tLJW67zVwm0auXVkfcwEsl6NbgzyTVN95jzoMlFjcrT4LnzkXZyZPQhYVxFguihkpEe2jOsm44Odn8H1xfXy0JHjHCPDrcqpVjYyYiokaj6YwgW+D/9Kg+8f6yn5dSuLxhg3VCnJ2tNbq6ag/OlSfCvXppo8VOXCpBtxZ/Jqm+8R5zHhxBJqKGqaREm0XCYs7hfAD47W+19g4dtO/LSyW6dAH0ekdGTERETQwTZCKqPyLaTDGWdcMpKcDVq1q7nx/QsycS09Px0vr1QI8eWvkEERGRA7HEgqiONen769IlICnJulTi7FmtTa8Hune3LpVo1w5QqmlfM6p3vL+ovvEecx4ssSCi+nX9OrB/v/UUa4cOaaPGABAZCQwcaE6Iu3TR6omJiIgaOCbIRFQ9EW3FSctSiT17gKIird3fX0uER4/WvhqN2qIcRERETogJMhFVlJ+vJcPlCfHu3cC5c1qbmxsQGwtMmWJ+kC4sjLNKEBFRo8EEmaipu3YN2LfPulTi8GFz+x13AIMHm+uGO3fWVqkjIiJqpJggEzUlIsCJE1ZTrCE1FSgu1tpbt9aS4LFjtYS4Rw+gZUvHxkxERHSLMUEmaszy8iqWSpw/r7W5u2ulElOnmh+kCw1lqQQRETV5TJCJGoviYmDvXusp1o4e1dqUAjp1Au6/31w3HB0NNONfAURERDfjv45EzkgEOH7cum44NVWrJwaANm20JHjCBPOsEl5ejo2ZiIjISTBBJnIGFy5YT7G2e7dWPgEAHh5aAjx9uvlBuuBglkoQERHVEBNkooamqAhIS7NOiI8f19qUAqKigOHDzaUSnTqxVIKIiKgO8V9VIkcS0eqELUeG09K0VeoAoG1bLQl+4gnta2ws0KKFY2MmIiJq5JggE91KubnWI8NJScDFi1qbp6c2rdozz5hHh4OCHBsvERFRE8QEmai+XL2qPThnmRCfOKG16XTaLBIjR5rrhjt1AlxcHBszERERMUEmqhNlZcCRI8CuXXgH0B6a27sXKCnR2oODtSR4yhQtIY6NBQwGR0ZMREREVWCCTFQT585Zr0aXlARcugQAGAtoU6rNnGlegKNtW4eGS0RERLZjgkxUncJCrVTCMiE+eVJrc3EBOncGRo0y1Q17R0ejdONGx8ZMRERENcYEmchSWRlw6JB13fC+fUBpqdYeGqolwtOmaSPD3btrD9dZduGAsImIiKjuMEGmpi0nx3pp5qQk4PJlrc3LS5tVYvZsc6lEYKBj4yUiIqJ6xwSZmo7CQiAlxTohzsrS2po1A7p0AUaPNk+xFhmpzTZBRERETQoTZGqcSku1UgnLBTj27zeXSoSHA717a8sz9+qllUq4uzs0ZCIiImoYmCBT43DmjHXdcHIycOWK1taypVYe8fzzWjLcowfQurVj4yUiIqIGiwkyOZ+CAq1UwjIhPnVKa2vWDOjWDXjsMXOpRIcOLJUgIiIimzFBpoattBRIT7euGz5wQJttAgDatwf69jWvRhcTA7i5OTZmIiIicmpMkKlhOX3aer7h5GTgl1+0Nh8fLREeNkz72rMn4O/v2HiJiIio0WGCTI5z5Yp5Vony15kzWlvz5lqpxIQJ5lKJiAhAKcfGTERERI2ezQmyUsoFQDKA0yJy301tegB/AxAL4AKAUSKSWYdxkrMrKQEOHrSuG05PN5dKREQAcXHmZLhbN0Cvd2jIRERE1DTZM4I8HcDPALwqaZsE4KKIRCilHgHwGoBRdRAfOSMR7aE5yynWkpO1eYgBwNdXS4JHjjSXSrRq5diYiYiIiG6wKUFWSgUDGAJgIYBnK9nlAQCJN77/EsA7SiklIlIXQVIDd/mylgBbJsTZ2Vqbq6v24Nzjj5tXo7vtNpZKEBERUYNl6wjyEgDPAWhRRXsQgP8BgIiUKKUuAWgF4LzlTkqpeADxABAaGlqTeMnRSkq0BTcsSyV+/lkbNQaA29EsZ1oAABF1SURBVG8Hfvtbc6lE165akkxEdKutXg3MnYtSQFscaOFCYMwYBwdFRM6g2gRZKXUfgHMikqKUiqvNwUTkfQDvA4DRaOTockMnoi3FbJkMp6QAV69q7X5+WhI8apR5AQ5fX8fGTEQEaMlxfDxQWAgdAJw8qb0HmCQTUbVsGUHuA2CoUmowADcAXkqpVSLye4t9TgMIAXBKKdUMQEtoD+uRM7l0CUhKMpVJFGzYAMONKdaKAOQGByPkySfNpRLt2rFUgogaprlzzc89lCss1LYzQSaialSbIIvI8wCeB4AbI8gzb0qOAWAdgHEAfgIwEsBG1h83cNeva6USllOsHTpkbu/YEYYRI4BevWB8+mkkFxcjhKUSROQssrLs205EZKHG8yArpV4CkCwi6wB8AOATpdQxAHkAHqmj+KguiACZmdalEnv2AEVFWntAgDYqPGaMuVTC29v08ZSnn2YdMRE5l9BQrayisu1ERNWwK0EWkc0ANt/4fp7F9iIAD9VlYFQL+fnmZLj8a26u1ubmBsTGAk89ZV6eOSyMpRJE1LgsXGiqQTbx8NC2ExFVgyvpObtr14B9+6ynWDt8WGtTCujYERgyxFw33LmztkodEVFjVl5nPHcuyk6ehC4sjLNYEJHNlKNKhY1GoyQnJzvk2EopOGWJtAiQkWFdKpGaChQXa+2tW5unV+vVCzAagZYta31Yp71eDsLrZT9eM6pPvL+ovvEecx5KqRQRMVa3H0eQG7K8POtSid27gfM3ppZ2d9cS4KlTzQlxSAhLJYiIiIhqiQlyQ1FcDOzda10qcfSo1qYU0KkTMHSouW44Ohpoxj8+IiIiorrGDMsRRIBjx6xLJdLStHpiAGjTRkuCJ07UEmKjEfDycmzMRERERE0EE+Rb4fx5c4lE+ehwXp7W5uGhJcDTp5sfpAsOZqkEERERkYMwQa5rRUXaaLBlqcTx41qbTgdERQHDh5vrhjt1YqlEY7F6NTB3LkoBIDycT8wTERE5KWZmtVFWppVKWK5Gt3evtkodAAQFaUnwE09oX2NjgRYtHBsz1Y/Vq01zruoAbYGC+HitjUkyERGRU+E0b/bIzbVefGP3bm1RDgDw9NRWoCsvk+jVS0uQGwFOX2OD8PDKV+0KC9NWMaRfxXuM6hPvL6pvvMecB6d5q62rV7U5hi0T4hMntDadTptF4qGHzAlxp06Ai4tjYybHycqybzsRERE1WEyQAa1U4sgR67rhvXuBkhKtPSRES4KnTDGXSnh6OjZmalhCQysfQQ4NvfWxOJHExEQsWLAAgDYCM3/+fCQmJjo2KCIiavKaZIlFa6Vw9p//NI8MJyUBly5pjS1aaKUS5WUSvXpp0641YfzVkQ0sapBNPDyA999nDTKRg/HvMKpvvMecB0ssKnPkCHDvvTgLAA88oJVEdO4MPPKIOSHu2JGlEmS/8iR47lyUnTwJXVgYZ7EgIiJyUk0rQQ4JAXr1QsLJk3hjyxage3eWSlDdGTMGGDMGLkpB+GAeERGR09I5OoBbyt0dWLMG/w8A+vZlckxEREREFTStBJmIiIiIqBpMkImIiIiILDBBJiIiIqqBxMREKKUAaDNZcJrKxqNJTvPG6Vjsw+tlH14vooaFP5P/v737jbHsrOsA/v2xu6bdCtTQwVSWUl8AwTSxxUlFq4S0lhQkYDQxkGKiMVYTQgqiRqzBKYYXJIb4Rk0KVaoukMof01TANrGoGLt1trSF/glRaGsB3VWpuK5CWX6+uLfJwzq7M9PO7Lmz+/kkN3vvmWfu/ebJ7Mz3nvucc4AnbfQ0b/YgAwDAQEEGAICBggwAAAMFGQAABgoyAAAMFGQAABgoyAAAMFCQAQBgoCADAMBAQQYAgIGCDAAAAwUZAAAGCjIAAAzOqIK8srKSqkqSVFVWVlamDQQAwMKp7p7khZeXl3t1dXWS12ZzqipT/ZzsROYLFov/k8CTqupgdy+vN+6M2oMMAADrUZABAGCgIAMAwEBBBgCAgYIMAAADBRkAAAbrFuSqOquq7qqqe6vq/qq6fo0xF1TVHVX1maq6r6pevT1xAQBge21kD/LXk1ze3d+f5OIkV1XVy44b85tJbu7uS5K8Psnvb21MAAA4NXavN6BnZ1c/Mn+4Z347/ozrneRZ8/vPTvLlrQoIAACn0obWIFfVrqq6J8mhJLd394HjhqwkeWNVPZbk40nevKUpAQDgFNlQQe7uY919cZJ9SS6tqouOG/KGJO/v7n1JXp3kT6rq/z13VV1TVatVtXr48OGnmx0AALbcps5i0d2PJ7kjyVXHfennk9w8H/P3Sc5Kct4a339Ddy939/LS0tJTSwwAANtoI2exWKqqc+f3z05yZZKHjhv2aJIr5mNekllBtosYAIAdZ92D9JKcn+SmqtqVWaG+ubtvrap3Jlnt7luSvC3Je6vqrZkdsPez84P7AABgR9nIWSzuS3LJGtvfMdx/IMllWxsNAABOPVfSAwCAgYIMAAADBRkAAAYKMgAADBRkAAAYKMgAADBQkAEAYKAgAwDAQEEGAICBggwAAAMFGQAABgoyAAAMFGQAABgoyAAAMFCQAQBgoCADAMBAQQYAgIGCDAAAAwWZE9u/P7nwwhxLkgsvnD0GADjN7Z46AAtq//7kmmuSo0dn76IeeWT2OEmuvnrKZAAA28oeZNZ23XXJ0aPfvu3o0dl2AIDTmILM2h59dHPbAQBOEwoya7vggs1tBwA4TSjIrO1d70r27v32bXv3zrYD7AArKyupqiRJVWVlZWXaQMCOUd09yQsvLy/36urqJK/NBu3fn1x3Xb71yCN5xgteMCvHDtBbV1Vlqv9XAMCJVdXB7l5eb5yzWHBiV1+dXH11dlWlH3546jQAAKeEJRYAADBQkAEAYKAgAwDAQEEGAICBggwAAAMFGQAABgoyAAAMFGQAABgoyLBFXNYWAE4PLjXNulw6GQA4HWz0UtP2IAMAwEBBBgCAgYIMAAADBRkAAAYKMgAADNYtyFV1VlXdVVX3VtX9VXX9Ccb9dFU9MB/zga2PCgAA22/3BsZ8Pcnl3X2kqvYk+XRVfaK773xyQFW9MMnbk1zW3V+tquduU14AANhW6xbknp0A98j84Z757fiT4v5Ckt/r7q/Ov+fQVoYEAIBTZUNrkKtqV1Xdk+RQktu7+8BxQ16U5EVV9XdVdWdVXXWC57mmqlaravXw4cNPLzkAAGyDDRXk7j7W3Rcn2Zfk0qq66Lghu5O8MMkrkrwhyXur6tw1nueG7l7u7uWlpaWnlxwAALbBps5i0d2PJ7kjyfF7iB9Lckt3P9HdX0zy+cwKMwAA7CgbOYvF0pN7g6vq7CRXJnnouGF/ntne41TVeZktufjCliYFAIBTYCNnsTg/yU1VtSuzQn1zd99aVe9MstrdtyT5yySvrKoHkhxL8qvd/e/blhoAALZJzU5SceotLy/36urqJK/N5lRVpvo5AQDYKlV1sLuX1xvnSnoAADBQkAEAYKAgAwDAQEEGAICBggwAAAMFGQAABgoyAAAMFGQAABgoyAAAMFCQAQBgoCADAMBAQQYAgIGCDAAAAwUZAAAGCjIAAAwUZAAAGCjIAAAwUJABAGCgIAMAwEBBBgCAgYIMAAADBRkAAAYKMgAADBRkAAAYKMgAADBQkAEAYKAgAwDAQEEGAICBggwAAAMFGQAABgoyAAAMFGQAABgoyAAAMFCQAQBgoCADAMBAQQYAgIGCDAAAAwUZAAAGCjIAAAwUZAAAGKxbkKvqrKq6q6rurar7q+r6k4z9qarqqlre2pgAAHBq7N7AmK8nuby7j1TVniSfrqpPdPed46CqemaSa5Mc2IacAABwSqy7B7lnjswf7pnfeo2hv53k3Un+d+viAQDAqbWhNchVtauq7klyKMnt3X3guK+/NMnzu/sv1nmea6pqtapWDx8+/JRDAwDAdtlQQe7uY919cZJ9SS6tqoue/FpVPSPJe5K8bQPPc0N3L3f38tLS0lPNDAAA22ZTZ7Ho7seT3JHkqmHzM5NclORTVfVwkpclucWBejvfyspKqipJUlVZWVmZNhAAwClQ3WstJx4GVC0leaK7H6+qs5PcluTd3X3rCcZ/KsmvdPfqyZ53eXm5V1dPOgQAALZMVR3s7nV34m5kD/L5Se6oqvuS/ENma5Bvrap3VtVrn25QAABYJOue5q2770tyyRrb33GC8a94+rEAAGAarqQHAAADBRkAAAYKMgAADBRkAAAYKMgAADBQkAEAYKAgAwDAQEEGAICBggwAAAMFGQAABgoyAAAMFGQAABgoyAAAMFCQAQBgUN09zQtXHU7yyCQvnpyX5N8meu2dyHxtjvnaPHO2OeZrc8zX5pmzzTFfmzfVnL2gu5fWGzRZQZ5SVa129/LUOXYK87U55mvzzNnmmK/NMV+bZ842x3xt3qLPmSUWAAAwUJABAGBwphbkG6YOsMOYr80xX5tnzjbHfG2O+do8c7Y55mvzFnrOzsg1yAAAcCJn6h5kAABYk4IMAACDM6ogV9UfVtWhqvrc1Fl2gqp6flXdUVUPVNX9VXXt1JkWWVWdVVV3VdW98/m6fupMO0FV7aqqz1TVrVNn2Qmq6uGq+mxV3VNVq1PnWXRVdW5VfbiqHqqqB6vqh6bOtKiq6sXzn6snb1+rqrdMnWvRVdVb57/zP1dVH6yqs6bOtMiq6tr5XN2/yD9fZ9Qa5Kp6eZIjSf64uy+aOs+iq6rzk5zf3XdX1TOTHEzyE939wMTRFlJVVZJzuvtIVe1J8ukk13b3nRNHW2hV9ctJlpM8q7tfM3WeRVdVDydZ7m4XJdiAqropyd929/uq6juS7O3ux6fOteiqaleSLyX5we6e6qJeC6+qnpfZ7/rv6+7/qaqbk3y8u98/bbLFVFUXJflQkkuTfCPJJ5P8Unf/46TB1nBG7UHu7r9J8h9T59gpuvsr3X33/P5/JXkwyfOmTbW4eubI/OGe+e3MeQf6FFTVviQ/nuR9U2fh9FNVz07y8iQ3Jkl3f0M53rArkvyTcrwhu5OcXVW7k+xN8uWJ8yyylyQ50N1Hu/ubSf46yU9OnGlNZ1RB5qmrqguTXJLkwLRJFtt8ucA9SQ4lub27zdfJ/W6SX0vyramD7CCd5LaqOlhV10wdZsF9b5LDSf5ovoznfVV1ztShdojXJ/ng1CEWXXd/KcnvJHk0yVeS/Gd33zZtqoX2uSQ/WlXPqaq9SV6d5PkTZ1qTgsy6quo7k3wkyVu6+2tT51lk3X2suy9Osi/JpfOPk1hDVb0myaHuPjh1lh3mR7r7pUleleRN86VjrG13kpcm+YPuviTJfyf59WkjLb75UpTXJvmzqbMsuqr6riSvy+zN2PckOaeq3jhtqsXV3Q8meXeS2zJbXnFPkmOThjoBBZmTmq+l/UiS/d390anz7BTzj3HvSHLV1FkW2GVJXjtfU/uhJJdX1Z9OG2nxzfdYpbsPJflYZmv5WNtjSR4bPsn5cGaFmZN7VZK7u/tfpw6yA/xYki929+HufiLJR5P88MSZFlp339jdP9DdL0/y1SSfnzrTWhRkTmh+0NmNSR7s7vdMnWfRVdVSVZ07v392kiuTPDRtqsXV3W/v7n3dfWFmH+f+VXfb83ISVXXO/IDZzJcKvDKzjyxZQ3f/S5J/rqoXzzddkcRBxut7Qyyv2KhHk7ysqvbO/2ZekdnxOpxAVT13/u8Fma0//sC0ida2e+oAp1JVfTDJK5KcV1WPJfmt7r5x2lQL7bIkP5Pks/N1tUnyG9398QkzLbLzk9w0P/r7GUlu7m6nLmMrfXeSj83+Dmd3kg909yenjbTw3pxk/3zZwBeS/NzEeRba/I3XlUl+ceosO0F3H6iqDye5O8k3k3wmC34J5QXwkap6TpInkrxpUQ+cPaNO8wYAAOuxxAIAAAYKMgAADBRkAAAYKMgAADBQkAEAYKAgAwDAQEEGAIDB/wHWjdjZE1aR+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = {'alpha0':   [alpha0_fit, sigma_alpha0_fit],\n",
    "     'alpha1':   [alpha1_fit, sigma_alpha1_fit],\n",
    "     'Chi2':     Chi2_fit,\n",
    "     'ndf':      Ndof_fit,\n",
    "     'Prob':     Prob_fit,\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.95, text, ax, fontsize=14)\n",
    "fig.tight_layout()\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (save_plots) :\n",
    "    fig.savefig(\"FitToLine.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Many experiments\n",
    "In case we have more than one \"experiment\" (i.e. fit), we would like to see how the fit results and the Chi2 values distribute themselves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAALICAYAAADseNpmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XuUVNWd8P3vT1FpabDBIA8inUYQGcR4ST8i4TWixugbMt7GeI9iMCSTy2TULEXfGE0yySIZHx0ZM1EmJJqLAWUwGkmeeEFNNIYoGkRUFBgItCgGbaANKMJ+/+iipxuqL0B11emu72etWl119jl1fvt07XPqV2effSKlhCRJkiRJKq09Sh2AJEmSJEkyQZckSZIkKRNM0CVJkiRJygATdEmSJEmSMsAEXZIkSZKkDDBBlyRJkiQpA0zQpU4SEYsiYlyp4yiliDgzIlZGRENEHLWTy94QET8r9LySJOXjcdvjtpQFJujSLoiI5RHxse2mTYiIJ7a9TikdllJ6rJ33qYmIFBE9OinUUrsR+FJKqTKl9Fypg2lNRFwQESsi4p2I+GVE9Ct1TJKkwvG43WGZP25HxMCIuD8iXsv9L2pKHZNUSCboUjeWgS8QHwQWlTiGNkXEYcDtwKeBAcDfgP8oaVCSpLLkcbtDtgL/F/iHUgcidQYTdKmTNP+1PiKOiYhnImJ9RLwRETflZvtd7m99rjvZmIjYIyK+ljujuyYifhIR+zV734tzZWsj4rrt1nNDRMyKiJ9FxHpgQm7dT0VEfUSsjohbI2LvZu+XIuILEfFqRGyIiG9FxNCI+EMu3rubz79dHfPGGhH7REQDsCewICKWtrL8LbmudOsjYn5EHNfKfNvOWEzK/WK+OiK+ut1se+fWvyHXTbG22fKTI2JpruzFiDiz2XIXAr9KKf0updQAXAecFRG988UiSeqePG53jeN2SumNlNJ/AE/n/UdKXZwJulQctwC3pJT6AEOBu3PTP5r7W5XrTvYUMCH3OAE4GKgEbgWIiJE0nt29EBgI7AcM2m5dpwOzgCrg58AW4HLgA8AY4CTgC9stcwrwYeBY4CpgGnARMBgYBZzfSr3yxppSejelVJmb54iU0tBWln8aOBLoB9wF3BMRPVuZl9x6DgE+DlwdLbsrngbMyNX7fnLbLGcpcByN2+sbwM8iYmCu7DBgwbYZU0pLgfeA4W3EIUnq3jxu55eF47bUrZmgS7vul7lft+sjop62u0VvBoZFxAdSSg0ppT+2Me+FwE0ppWW5M7rXAOdFY7e3s2k82/tESuk94OtA2m75p1JKv0wpbU0pbUwpzU8p/TGl9H5KaTmN3bmP326Z76WU1qeUFgEvAA/m1r8O+A3Q2kAxbcXarpTSz1JKa3Ox/R9gH+DQNhb5RkrpnZTSQuDHtPwC8kRK6dcppS3AT4Ejmq3nnpTSa7ltMhN4FTgmV1wJrNtuPesAz6BLUvficbt7HLelbs0EXdp1Z6SUqrY92PHX7eYm0nhG9uWIeDoiPtnGvAcCK5q9XgH0oPH66AOBldsKUkp/A9Zut/zK5i8iYnhEPBARr+e6z32Hxl/lm3uj2fONeV5Xkl9bsbYrIr4aES9FxLrcl6X98sTWXPO6rcitf5vXmz3/G9Bz2xeOXPfCPzf7Ujaq2XoagD7bracPsKEjdZAkdRket7vHcVvq1kzQpSJIKb2aUjofOAD4LjArInqx46/oAK/ROEjLNtXA+zQefFcDB20riIgKYP/tV7fd6x8ALwOH5LrqXQvErtemw7G2KXfd2lXAOUDf3Jelde3ENni7db3WgfV8EPhP4EvA/rn1vNBsPYto9qt9RBxM4xmBV9p7b0lS9+Rxe0cZOm5L3ZoJulQEEXFRRPRPKW0F6nOTtwJv5v4e3Gz2XwCXR8SQiKik8ZfzmSml92m8Ru3vI+IjuQFgbqD9A1ZvYD3QEBEjgH8sVL3aibU9vWn8UvAm0CMivs6OZ7K3d11E7BuNI69fCszswHq2faF6EyAiLqXxl/htfk7jNj0u9+Xrm8DslJJn0CWpTHncbjWuLBy3yV33vk/u5T7tXAcvdSkm6FJxnAosyo2QegtwXu46s78B3waezHXjOhb4EY3XYv0O+G9gE/BlgNy1Zl+mcVCV1TR2z14DvNvGur8KXEBjl+3/pGMHx45qNdYO+C2Nt0l5hcZub5vYrptfHo8DS4BHgBtTSg+2t5KU0ovA/wGeovEMweHAk83KFwGfpzFRX0PjF5C2uj1Kkro/j9s7ysRxO2cjjdsSGnsbbOxYFaTsi5Ty9dSR1BXkfv2up7Eb3H+XOp7OEhE1NH6R2KuDv/JLkpQ5Hrcltccz6FIXExF/n+su1gu4EVgILC9tVJIkKR+P25J2hgm61PWcTuMgK6/ReG/R85JdYSRJyiqP25I6zC7ukiRJkiRlgGfQJUmSJEnKgB6lDmB7H/jAB1JNTU2pw5AkKdPmz5//15RS/1LH4XFbkqT2dfS4nbkEvaamhmeeeabUYUiSlGkRsaLUMYDHbUmSOqKjx227uEuSJEmSlAEm6JIkSZIkZYAJuiRJkiRJGWCCLkmSJElSBpigS5IkSZKUASbokiRJkiRlgAm6JEmSJEkZYIIuSZIkSVIGmKBLkiRJkpQBJuiSJEmSJGWACbokSWUoIqoiYlZEvBwRL0XEmIjoFxEPRcSrub99Sx2nJEnlxARdkqTydAvwf1NKI4AjgJeAycAjKaVDgEdyryVJUpGYoEuSVGYiYj/go8B0gJTSeymleuB04M7cbHcCZ5QmQkmSylOPUgcgqXTGTplLXf3GFtMGVVXw5OQTSxSRpCIZArwJ/DgijgDmA18BBqSUVufmeR0YkG/hiJgETAKorq7u/GilMpXvOL2Nx2upezJBl8pYXf1Glk8Z32JazeQ5JYpGUhH1AI4GvpxSmhcRt7Bdd/aUUoqIlG/hlNI0YBpAbW1t3nkk7b58x+ltPF5L3ZNd3CUpwz7/+c/zrW99q9RhqPtZBaxKKc3LvZ5FY8L+RkQMBMj9XVOi+DLPtilJ6gwm6JLUjgkTJrD33ntTWVnZ9Lj++us7vPy4ceP44Q9/uEvrvu2227juuut2adli+pd/+RcOO+ww9thjD+64444WZatWrWLcuHHsu+++HH300bzwwgtNZZs3b2bixIn07t2b6upq7r777hbLTp06lQEDBtC3b1+uueaaYlSlLKSUXgdWRsShuUknAS8C9wOX5KZdAtxXgvA6zLbZPtumJHUtJuiS1AFXXXUVDQ0NTY9vfOMbpQ4pUw4++GCmTp3Khz/84R3KJk2axKhRo1i7di3nnnsu5557blPZzTffzMKFC1m5ciU/+clP+MxnPsPKlSsBmDdvHjfccANz587lhRdeYMaMGTskCdotXwZ+HhHPA0cC3wGmACdHxKvAx3KvM8222TbbpiR1LSbokrQb7rjjDsaMGcOECRPo3bs3Y8eO5a9//SsA3/nOd6isrOT3v/89X/rSl6isrNzhS3JNTQ233norH/7wh+nVqxennXYaAA888ACVlZXstddefO1rX9thvbNnz2bUqFH07duXT3ziE6xevbqp7IEHHuDQQw+ld+/eHHLIIfz2t7/txC3Q6IILLuCkk05in332aTF9/fr1PPjgg0yePJmKigouv/xyVqxYwfPPPw/APffcwz/90z9RVVXFuHHjGDNmDPfeey8As2bN4qyzzuKwww5j0KBBXHbZZcyYMaPT61IuUkp/TinVppQ+lFI6I6X0dkppbUrppJTSISmlj6WU3ip1nLvKttnItilJXYsJuiTtpqeffpoLL7yQNWvW8P777/OjH/0IgGuvvZaGhgaOO+44br31VhoaGpg/f/4Oy99+++385Cc/Yd26dVx77bUAfPKTn6ShoYELL7xwh/n/9Kc/MXHiRKZPn86bb77JUUcdxaRJk5rKP/vZz/LNb36TDRs28OCDDzJo0KBOqnn7lixZQs+ePamsrOS4445j2bJlDB06lMWLFwOwePFiRowYwUUXXcSMGTMYOXLkDmW33HILV155ZYsyqSNsm62zbUpSNpmgS1IH3HjjjVRVVTU9Zs+e3VQ2fPhwTj75ZCoqKjjppJN45ZVXduq9J02axGGHHUaPHj049thj251/+vTpTJgwgdGjR9OjRw+++tWvMmfOHN59910A9thjD5YuXcr69esZMmQIo0aN2rnKFtA777xDZWUlGzZsYNGiRbz99tv07t2bhoaGFuULFy6krq4ub9nSpUt59dVXW5RJ29g2d41tU5KyyQRdkjrgq1/9KvX19U2Ps846q6msX79+Tc/33ntvNm3atFPvfcghh+zU/CtXruT2229vSkiGDBnC3nvv3dSV9p577uGpp56iurqaY489tsXAT8XWq1cvGhoaGDx4MG+99RZjxoxhw4YNVFZWtihfsGABV155Zd6yqVOncv/997cok7axbe4a26YkZZMJuiR1sj32aHtX26NHj516v8GDB3Pddde1SEo2bdpETU0NAB/5yEf41a9+xRtvvMGIESPyXidbLMOGDWPjxo2sWrUKgPfee4+lS5dy6KGNg4cPHz6cl19+uWn+F198sUNlUiHYNm2bkpQ1JuiS1Mn+1//6XwU9UzZhwgRuu+02nn32WVJKrFmzhpkzZwKwdetWfvazn7Fhw4am5KNPnz4FW3drNm/ezKZNm0gpNT3funUrffr04ZRTTmHKlCls2rSJm2++merqag4//HAAzjnnHKZOnUp9fT2PP/44Tz31FGeeeSYAn/rUp5g9ezaLFi2irq6O6dOntxhlWtpdtk3bpiRljQm6JHXA9773vRb3Wr7ssss6vOyVV17JQw89xIEHHsi4ceM6tMzHP/5xKisr+fnPf9607gkTJgAwZswYbrrpJi699FL2228/amtrWwxw9bOf/YzBgwfzgQ98gNdee41vf/vbO1PVXfLZz36WiooK/vCHPzBp0iQqKir43e9+BzQOtLVw4UL69u3LjBkzmDlzJhEBwOWXX86oUaMYPHgwF110EdOnT2fw4MEAjB49muuvv54TTjiBUaNGcc4555gEaAe2zbbZNiWpa4mUUqljaKG2tjY988wzpQ5DKgs1k+ewfMr4dqdJyp6ImJ9Sqi11HB63pc7T1jHZ47XUtXT0uO0ZdEmSJEmSMsAEXZIkSZKkDNi54UklSZIkZdrYKXOpq9+Yt2xQVQVPTj6xyBFJ6igTdEmSJKkbqavf2Oa165Kyyy7ukiRJkiRlQEET9IjYMyKei4gHcq+HRMS8iFgSETMjYu9Crk+SOttjjz1GRHD11Vc3TTvmmGOabkVULh599FEOP/xw+vTpw8CBA7niiivYsmVLh5bduHEjn/vc5+jXrx9VVVV84QtfaCrbvHkzEydOpHfv3lRXV3P33Xe3WPbhhx/miCOOoFevXgwdOpQFCxYUtF7qumyb/2Pq1KkMGDCAvn37cs0113R4uVWrVjFu3Dj23Xdfjj766Bb3hP/rX//KmWeeSd++fTnwwAO54YYbmsree+89vvjFLzJgwAD2339/rrrqKrJ2VyBJ6qoKfQb9K8BLzV5/F7g5pTQMeBuYWOD1SVKn69+/Pw8//DAAf/nLX1i3bl2JIyq+ESNG8Otf/5p169bxyiuvsGjRIr7//e93aNnLL7+cJUuWsGjRItauXdsiQb/55ptZuHAhK1eu5Cc/+Qmf+cxnWLlyJQDLly/nrLPO4tprr2XdunU88cQTDBw4sFPqp67Jtgnz5s3jhhtuYO7cubzwwgvMmDFjhx+6WjNp0iRGjRrF2rVrOffcc1vcy/xrX/samzdv5rXXXuOPf/wj06ZNY86cxq7RP/jBD5g3bx4vv/wyzz33HDNnzuTnP/95p9RPkspNwRL0iDgIGA/8MPc6gBOBWblZ7gTOKNT6JKlY9t57b0aOHMlzzz3HrFmzOOuss1qUz549m1GjRtG3b18+8YlPsHr16qay8847jwEDBtCrVy/GjBnD888/31RWU1PD9ddfz7Bhw+jXrx+33XZb0eq0swYOHMjgwYOJCN5//322bt3KSy+91O5yGzdu5M477+Tf//3fGThwIHvuuSejRo1qKr/nnnv4p3/6J6qqqhg3bhxjxozh3nvvBeCOO+7g1FNP5dxzz6VHjx4MHDiQAw44oNPqqK7HtklTvQ877DAGDRrEZZddxowZM9pdbv369Tz44INMnjyZiooKLr/8clasWNG0HVasWMEnP/lJKioqqK6uZsyYMbz44osAPP7441xwwQX07duX6upqzj//fGbPnt2p9ZSkclHIM+j/BlwFbM293h+oTym9n3u9ChiUb8GImBQRz0TEM2+++WYBQ5JUCGOnzKVm8pwdHmOnzC11aEVz9tlnM2vWLH75y19y+umnN03/05/+xMSJE5k+fTpvvvkmRx11FJMmTWoqP+qoo1i4cCHr16/nox/9KJdcckmL9503bx7PP/8806dP56qrruL9998nq/7yl7+w33770a9fP+bNm8eECRPaXeaVV14hIvjlL3/JgAEDGDlyZFMCDrB48WJGjBjBRRddxIwZMxg5ciSLFy8GYMGCBfTr14/Ro0czYMAALrjggrI8Q6q2lXvb3NaGbrnlFq688soWbagtS5YsoWfPnlRWVnLcccexbNkyhg4d2rTsF77wBR588EHeeecdli9fzvz58zn55JMB8nZnX7JkSWErJkllqiAJekR8EliTUpq/K8unlKallGpTSrX9+/cvREiSCmjbaLDbP1q7hUt3dMopp/Bf//Vf7LPPPnzgAx9omj59+nQmTJjA6NGj6dGjB1/96leZM2cO7777LgBXX301BxxwAHvuuScTJkxocZYO4DOf+Qz77rsv48ePZ8OGDbz++utFrdfOqK6uZt26dSxdupQrrriCQw45pN1l1q9fz3vvvceyZctYsWIF3//+9/n0pz/ddCbznXfeobKykoULF1JXV0fv3r1paGhoWvbuu+/mtttuY9myZaxfv56vf/3rnVpHdT3l3ja3taGlS5fy6quvtmhDHVluw4YNLFq0iLfffrvFskcddRTr1q1jv/32Y8iQIVx22WUceeSRAIwbN4677rqLtWvXsnz5cu655x7+9re/dWo9JalcFOoM+ljgtIhYDsygsWv7LUBVRGy7ldtBQF2B1idJRdWzZ09OO+00Lr300hbTV65cye23305VVRVVVVUMGTKEvffem9WrV7NlyxauueYahg4dSlVVFaNHj2br1q0tBlfr168f0NhVF2DTpk3Fq9QuOvjggxkxYgRf/OIX25133333ZcuWLVx55ZX07NmTE044geHDh/PHP/4RgF69etHQ0MCCBQu48sor2bBhA5WVlU3LnnzyyRx11FH06tWLSZMm8dhjj3Vm1dQFlXvb3NaGpk6dyv3339+iDXVkucGDB/PWW28xZsyYFsued955HHHEEbzzzjusWLGCu+66q+na9n/8x3/k2GOP5fDDD2f8+PGcccYZ9OnTp1PrKUnloiAJekrpmpTSQSmlGuA8YG5K6ULgUeDs3GyXAPcVYn2SVArf+973uOCCC1pMGzx4MNdddx319fVNj02bNlFTU8Ndd93Fvffey6OPPkp9fT1PPPEEkL97aFeTUmLhwoXtznfwwQcTETuMrL1tGwwfPpyXX365afqLL77IoYceCsDQoUNbXU5qrpzbZlttqC3Dhg1j48aNrFq1CmgcmX3p0qVNyz777LNMmDCBffbZh+rqasaPH88jjzwCNP5oceutt/Laa6+xaNEiIoIjjjiiE2onSeWns++DfjVwRUQsofGa9OmdvD5JKqoJEyZw22238eyzz5JSYs2aNcycOROADRs2UFFRQd++fdmwYQPf+c53ShztrvvFL37B888/z9atW/nLX/7CjTfeyJgxY1rMM3v2bIYNG0Zd3f90lurbty/HH388N910E5s3b+b3v/89r7zyCsceeywA55xzDlOnTqW+vp7HH3+cp556ijPPPBOAM888k1//+tcsXLiQTZs2MX36dE444YTiVVpdWrm0zU996lPMnj2bRYsWUVdXx/Tp01uMxg7522afPn045ZRTmDJlCps2beLmm2+murqaww8/HGjs4v7Tn/6UzZs38/rrr/Pb3/6WkSNHAo3d45csWcLWrVt5+OGH+eEPf9ji+n7tqLWxXMptPBdJ7St4gp5Seiyl9Mnc82UppWNSSsNSSp9KKb1b6PVJUimNGTOGm266iUsvvZT99tuP2tpa5s9vHI7j4osvprq6mgMPPJAjjjiiKSntitauXcvZZ59Nnz59+N//+3/zoQ99iBtvvLHFPOvXr2fp0qVs3ry5xfQf/ehHLF68mKqqKiZOnMhPf/pTDjzwQKDxFmyjRo1i8ODBXHTRRUyfPp3BgwcDcPzxx/P1r3+dk08+mQMPPJCKigq++c1vFqfC6vLKpW2OHj2a66+/nhNOOIFRo0Zxzjnn7JCgt9Y2b7/9dhYuXEjfvn2ZMWMGM2fObOq1cscdd/DnP/+Z/v37c+SRR3L88cc33SLxb3/7G+PHj6eyspLPf/7z/Md//McOP9ippdbGcim38VwktS+y1p2rtrY2PfPMM6UOQyoLNZPnsHzK+F2a1tZ0SZ0vIuanlGpLHYfHbal9bR0vs1QmqfN09Ljdo70ZJGl3jZ0yd4czBIOqKnhy8oklikiSJEnKHhN0SZ1uW9e+5momzylRNJIkSVI2maBLkiRJuylfb7FtBlVVFHx9g6oqWv2xuzPWJ6k4TNAlSZKk3ZSvt1hn8jIxqXvq7NusSZIkSZKkDjBBl6R21NTUUFFRQWVlJUOHDuXWW29tKosIRo8e3fT6qquuIiJ47LHHALj//vs5/PDD6d27Nx/84Af59re/XezwC2bq1KkMGDCAvn37cs0113Romc2bN3PxxRczcOBA9ttvP0444QQWLVrUVL548WJOPfVUqqqqqKmpabHs1q1bOffccznooIOICJYvX17A2qg7sG022pW2CbBq1SrGjRvHvvvuy9FHH80LL7zQVHbzzTczZMgQ+vTpw7Bhw/jxj3/cVPb5z3+eysrKpsc+++zTdP90SdLuMUGXpA741a9+RUNDA7/4xS+4+uqrefzxx5vK6uvrWblyJQCPPPII/fv3B2DZsmVcfPHFfP/732f9+vX8/ve/Z9iwYSWJf3fNmzePG264gblz5/LCCy8wY8YM7r777naX27JlC8OGDePpp5/m7bff5rTTTuOMM85oKu/RowfnnXce//qv/5p3+Y985CPMmjWrYPVQ92Pb3LW2CTBp0iRGjRrF2rVrOffcc1vcP338+PHMnz+f9evX88gjj/C1r32NBQsWAHDbbbfR0NDQ9Dj77LP5h3/4h06pnySVG69Bl9RCvkFnHGzmfxxzzDEcdthhPP300xx//PEAnHXWWcyaNYvjjz+ekSNH8sYbbwDw7LPPMmTIED760Y8CUF1dTXV1dcli3x2zZs3irLPO4rDDDgPgsssuY8aMGZxzzjltLtezZ0++/vWvN72+9NJLueKKK3jzzTfp378/Q4cOZejQoTz88MM7LLvHHnvwla98hffff7+wlVG3ZNvcuba5fv16HnzwQZYvX05FRQWXX3453/rWt3j++ef50Ic+xPDhw5vmfe+999i6dSsvv/wyRxxxRIv3WbduHffee2+Ls++SpF3nGXRJLTw5+USWTxnf4uFANP9j/vz5LFq0iL/7u79rmnbGGWdw3333MWvWrBZnkQ4//HAWLVrE9ddfz4IFC0gplSLkgli8eDEjRozglltu4corr2TkyJEsXrx4p9/nqaee4sADD2T//ffvhChVzmybO9c2lyxZQs+ePamsrOS4445j2bJlDB06tMWyd911F7169WL48OH06dOHj3/84zu8z4wZMzj66KM5+OCDC1ovSSpXJuiS1AFnnHEG++67Lx/96EeZMmUK48f/z0i9+++/P/vssw+zZs3ilFNOaZp+6KGH8utf/5o//OEPjB49miFDhnDfffeVIvzd9s4771BZWcnSpUt59dVX6d27Nw0NDTv1HuvWreMrX/kKN910E3vs4eFHhWHb3LW2uW25DRs2sGjRIt5+++0dlr3gggtoaGjgiSee4LLLLmO//fbb4X1+/OMfc/HFFxe0TpJUzvyGJEkd8Mtf/pINGzbwz//8z/zhD3/YofzSSy/ltNNOo6Ki5eUAH/vYx3jooYd4++23ueqqqzj//PNZu3ZtscIumF69etHQ0MDUqVO5//772bBhA5WVlR1e/t133+XMM8/c4TpXaXfZNnetbW5bbvDgwbz11luMGTMm77IRwdixY1m6dCn/+Z//2aLspZdeYsGCBe12p5ckdZwJuiR10J577skNN9zAU0891TQS9DbnnXceN954Y6vLVlRU8IUvfIGePXuybNmyTo608IYPH87LL7/c9PrFF1/k0EMP7dCyW7Zs4fzzz2fYsGF861vf6qwQVcZsmzvfNocNG8bGjRtZtWoV0Hid+dKlS1tdNqXEwoULW0z78Y9/zN///d9TVVW1GzWQJDVngi5JO2Gvvfbiy1/+Mt/85jfbnXfBggVMnz6d+vp6tm7dyi9+8Qs2b97cYvClruJTn/oUs2fPZtGiRdTV1TF9+vQdzoTPnj2bYcOGUVdX12L65z73OSKCH/zgBzu8b0qJTZs2sXnz5qbn7733XlP5u+++y7vvvtv0fNOmTZ1QO3UHts2da5t9+vThlFNOYcqUKWzatImbb76Z6urqptulTZ06lVWrVpFS4qmnnmLGjBmMGTOmafktW7bw05/+1O7tBbBtcNZ8DwdplcqPo7hL0k6aOHFi09m6tlRWVnL33Xdz1VVXNX35nz17dt7rOLNu9OjRXH/99Zxwwgls3ryZz33uczskAevXr2fp0qVs3ry5adqKFSuYPn06FRUVLer9m9/8huOOO44VK1YwZMiQpukVFRUcf/zxTWdBDz30UFasWAHAiBEjALr0gF7qXLbNjrdNgNtvv52LLrqIvn37MmLECGbOnElEAI0/Ynz3u9+lvr6eAw44gKuvvpoLL7ywadnf/OY3bNmyhVNPPbXzK9nNORCrpOYia190amtr0zPPPFN2v4ogAAAgAElEQVTqMKSyUDN5DsunjG9/xlaMnTKXuvqNLaYNqqrY4ctGvvXs7rqlchcR81NKtaWOw+O21KirHNe6SpxSd9PR47Zn0CXtsny/+m9/D3VJkiRJHeM16JIkSZIkZYBn0CVJkqScfJdvbZPvMi5JKiQTdEmSylBELAc2AFuA91NKtRHRD5gJ1ADLgXNSSm+XKkapFOrqN7Z6jbaXcUnqbHZxlySpfJ2QUjqy2aA1k4FHUkqHAI/kXkuSpCLxDLrUhXV0FPW25pWkZk4HxuWe3wk8BlxdqmAkSSo3JuhSF5avG15r3e/a6rInqSwl4MGISMDtKaVpwICU0upc+evAgJJFJ0lSGTJBlySpPP0/KaW6iDgAeCgiXm5emFJKueR9BxExCZgEUF1d3fmRSpJUJrwGXZKkMpRSqsv9XQPcCxwDvBERAwFyf9e0suy0lFJtSqm2f//+xQpZkqRuzwRdkqQyExG9IqL3tufAx4EXgPuBS3KzXQLcV5oIJUkqT3ZxlySp/AwA7o0IaPwucFdK6f9GxNPA3RExEVgBnFPCGCVJKjsm6JIklZmU0jLgiDzT1wInFT8iSZIEJuhSyeW7/Vk+rd0+TZIkSVL3YIIulVhHb3/W2u3TJEmSJHUPDhInSZIkSVIGmKBLkiRJkpQBJuiSJEmSJGWACbokSZIkSRlggi5JkiRJUgY4irvUzQyqqsg74vugqooSRCNJkiSpo0zQpW7Ge6VLkiRJXZNd3CVJkiRJygATdEmSJEmSMsAEXZIkSZKkDDBBlyRJkiQpA0zQJUmSJEnKAEdxl7qIfLdP89ZpkiRJUvdhgi51Ed4+TZIkSere7OIuSZIkSVIGmKBLkiRJkpQBJuiSJEmSJGWACbokSZIkSRlggi5JkiRJUgY4irukksh327ht0x2xXpIkSeXIBF1SSbSWhOdL2iVJkqRyYBd3SZIkSZIywDPokiRJUplo7RKzbWVeZiaVlgm6JEmSVCbaSsC9zEwqPbu4S5IkSZKUASbokiRJkiRlgF3cJXVJY6fMpa5+Y4tpXjsnSZKkrswEXVKXVFe/keVTxreY5rVzkiRJ6srs4i5JkiRJUgYULEGPiJ4R8aeIWBARiyLiG7npQyJiXkQsiYiZEbF3odYpSZIkSVJ3Ucgz6O8CJ6aUjgCOBE6NiGOB7wI3p5SGAW8DEwu4TkmSJEmSuoWCJeipUUPu5V65RwJOBGblpt8JnFGodUqSJEmS1F0U9Br0iNgzIv4MrAEeApYC9Sml93OzrAIGFXKdkiRJkiR1BwVN0FNKW1JKRwIHAccAIzqyXERMiohnIuKZN998s5AhSZIkSZLUJXTKKO4ppXrgUWAMUBUR227ndhBQl2f+aSml2pRSbf/+/TsjJEmSJEmSMq2Qo7j3j4iq3PMK4GTgJRoT9bNzs10C3FeodUqSJEmS1F30aH+WDhsI3BkRe9KY+N+dUnogIl4EZkTEvwDPAdMLuE5JkiRJkrqFgiXoKaXngaPyTF9G4/XokiRJkiSpFZ1yDbokSZIkSdo5JuiSJJWp3O1Rn4uIB3Kvh0TEvIhYEhEzI2LvUscoSVI5KeQ16JLEoKoKaibP2WHa7i7/5OQTCxKfpBa+QuOArn1yr78L3JxSmhERtwETgR+UKjhJksqNCbqkgtrdRDrf8tsn7JJ2X0QcBIwHvg1cEREBnAhckJvlTuAGTNAlSSoaE3RJksrTvwFXAb1zr/cH6lNK7+derwIG5VswIiYBkwCqq6s7OUwpO/L18mpeVq7GTplLXf3GvGX2gpN2jgm6JEllJiI+CaxJKc2PiHE7u3xKaRowDaC2tjYVODwps0w086ur38jyKePzltkLTto5JuiSJJWfscBpEfEJoCeN16DfAlRFRI/cWfSDgLoSxihJUtlxFHdJkspMSumalNJBKaUa4DxgbkrpQuBR4OzcbJcA95UoREmSypIJuiRJ2uZqGgeMW0LjNenTSxyPJEllxS7ukiSVsZTSY8BjuefLgGNKGY8kSeXMM+iSJEmSJGWACbokSZIkSRlggi5JkiRJUgaYoEuSJEmSlAEOEidJkqRuZ+yUudTVb8xbNqiqgicnn1jkiCSpfSbokjJvUFUFNZPn7DBNkqTW1NVvZPmU8XnLtj+mSFJWmKBLyjzPckiSJKkceA26JEmSJEkZYIIuSZIkSVIGmKBLkiRJkpQBJuiSJEmSJGWACbokSZIkSRlggi5JkiRJUgaYoEuSJEmSlAHeB12SJEllZVBVBTWT57RaVq7a2y5PTj6xyBFJ5ccEXZIkSWXFRDO/trZLa4m7pMKyi7skSZIkSRlggi5JkiRJUgaYoEuSJEmSlAFegy5JkqQuaeyUudTVb8xbVs6DvUnqukzQJUmS1CXV1W9k+ZTxpQ5DkgrGLu6SJEmSJGWACbokSZIkSRlgF3dJkiRJmdLe+ALey17dlQm6JEmSpExpa3yBmslzihyNVDx2cZckSZIkKQM8gy6p2xhUVbHDr+p2g5MkSVJXYYIuqdvIl4jbDU6SJEldhQm6JEmSpDbl66XWvExSYZigS5IkSWqTl4tJxeEgcZIkSZIkZYAJuiRJkiRJGWCCLkmSJElSBngNusrK2Clzqavf2GJaZ9yGq1jrUWHl+7/tDP/HkiRJ2h0m6CordfUbWT5lfItpnXEbrmKtR4WV7/+2M/wfS5IkaXfYxV2SJEmSpAwwQZckSZIkKQNM0CVJkiRJygCvQZckSVKmtTaI56CqihJEI0mdxwRdkiRJmba7g3hKUldhgi4VyaCqiryjfPvrf2m0dis8qRxERE/gd8A+NH4XmJVSuj4ihgAzgP2B+cCnU0rvlS5SSZLKiwm6VCTeHztbPBujMvcucGJKqSEi9gKeiIjfAFcAN6eUZkTEbcBE4AelDFSSpHLiIHGSJJWZ1Kgh93Kv3CMBJwKzctPvBM4oQXiSJJUtz6BLklSGImJPGruxDwO+DywF6lNK7+dmWQUMamXZScAkgOrq6s4PVlKX1dolftvK7GEotWSCLklSGUopbQGOjIgq4F5gxE4sOw2YBlBbW5s6J0JJ3UFbCXhribtUzuziLklSGUsp1QOPAmOAqojY9uP9QUBdyQKTJKkMmaBLklRmIqJ/7sw5EVEBnAy8RGOifnZutkuA+0oToSRJ5cku7pJUZK3d4s3r8FREA4E7c9eh7wHcnVJ6ICJeBGZExL8AzwHTSxmkJEnlxgRdkoos3y3evA5PxZRSeh44Ks/0ZcAxxY9IkiSBXdwlSZIkScqEgiXoETE4Ih6NiBcjYlFEfCU3vV9EPBQRr+b+9i3UOiVJkiRJ6i4KeQb9feDKlNJI4FjgixExEpgMPJJSOgR4JPdakiRJkiQ1U7AEPaW0OqX0bO75BhpHgx0EnA7cmZvtTuCMQq1TkiRJkqTuolOuQY+IGhoHn5kHDEgprc4VvQ4M6Ix1SpIkSZLUlRV8FPeIqAT+C/jnlNL6iGgqSymliEh5lpkETAKorq4udEhSp2rtllmSJEmStDMKmqBHxF40Juc/TynNzk1+IyIGppRWR8RAYM32y6WUpgHTAGpra3dI4KUsy3fLLEmSJEnaWYUcxT2A6cBLKaWbmhXdD1ySe34JcF+h1ilJkiRJUndRyDPoY4FPAwsj4s+5adcCU4C7I2IisAI4p4DrlCRJkiSpWyhYgp5SegKIVopPKtR6JEmSJEnqjjplFHdJkiRJkrRzCj6Ku9TVDKqqoGbynB2mPTn5xBJFJEmSJKkcmaCr7OVLxLdP2CVJkiSps9nFXZIkSZKkDPAMuiRJkkpu7JS51NVvzFs2qKqiyNFIUmmYoEuSJKnk6uo3snzK+FKHIUklZRd3SZIkSZIywARdkiRJkqQMMEGXJEmSJCkDTNAlSZIkScoAE3RJkiRJkjLABF2SJEmSpAwwQZckSZIkKQNM0CVJkiRJygATdEmSJEmSMqBHqQOQupKxU+ZSV7+xxbRBVRUlikYdMaiqgprJc/JOlyRJkrLEBF3aCXX1G1k+ZXypw9BOeHLyiaUOQZIkSeoQu7hLkiRJkpQBJuiSJEmSJGWACbokSZIkSRlggi5JkiRJUgaYoEuSJEmSlAEm6JIkSZIkZYAJuiRJkiRJGWCCLkmSJElSBpigS5IkSZKUASbokiRJkiRlgAm6JEmSJEkZYIIuSZIkSVIG9Ch1AFIWDaqqoGbynLzTpdbk+9wMqqrgycknligiSZIkdSUm6FIeJlTaFfk+N/l+6JFKLSIGAz8BBgAJmJZSuiUi+gEzgRpgOXBOSuntUsUpSVK5sYu7JEnl533gypTSSOBY4IsRMRKYDDySUjoEeCT3WpIkFYkJuiRJZSaltDql9Gzu+QbgJWAQcDpwZ262O4EzShOhJEnlyQRdkqQyFhE1wFHAPGBASml1ruh1GrvAS5KkIvEadEmSylREVAL/BfxzSml9RDSVpZRSRKRWlpsETAKorq4uRqjqJsZOmUtd/ca8ZQ7EKkkm6JIklaWI2IvG5PznKaXZuclvRMTAlNLqiBgIrMm3bEppGjANoLa2Nm8SL+VTV7+R5VPGlzoMScosE3RJ6kSt3XpNKqVoPFU+HXgppXRTs6L7gUuAKbm/95UgPEmSypYJuiR1Im/Zp4waC3waWBgRf85Nu5bGxPzuiJgIrADOKVF8kiSVJRN0SZLKTErpCSBaKT6pmLGo+/E6c0nadSbokiRJKhivM5ekXedt1iRJkiRJygATdEmSJEmSMsAEXZIkSZKkDPAadEmSJElFl+9WpM3LdkV7gxR6dxVlnQm6JEmSpKLrjGS5rUEKW/sxQMoSu7hLkiRJkpQBJuiSJEmSJGWACbokSZIkSRngNeiSJEnaKe0NxCVJ2jUm6JIkSdopbQ3EJUnadSbokpQB+W414+1gJEmSyosJuiRlQL5E3NvBSJIklRcHiZMkSZIkKQNM0CVJkiRJygATdEmSJEmSMsAEXZIkSZKkDHCQOHULbd2PtTnvzSpJkiQpq0zQ1S14P1ZJkiRJXZ1d3CVJkiRJygATdEmSJEmSMqBgCXpE/Cgi1kTEC82m9YuIhyLi1dzfvoVanyRJkiRJ3Ukhz6DfAZy63bTJwCMppUOAR3KvJUmSJEnSdgqWoKeUfge8td3k04E7c8/vBM4o1PokSZIkSepOOnsU9wEppdW5568DA/LNFBGTgEkA1dXVnRySJHVt+W4rOKiqgicnn1iiiCRJKp5BVRXUTJ7TapnUlRXtNmsppRQRqZWyacA0gNra2rzzSJIa5butYGtfVCRJ6m78QVrdWWeP4v5GRAwEyP1d08nrkyRJkiSpS+rsBP1+4JLc80uA+zp5fZIkSZIkdUmFvM3aL4CngEMjYlVETASmACdHxKvAx3KvJUmSJEnSdgp2DXpK6fxWik4q1DokSZIkqdDyDcC6jQOxqpiKNkicJEmSJGVRvgFYt3EgVhVTZ1+DLkmSJEmSOsAEXZIkSZKkDDBBlyRJkiQpA7wGXZIkSVK3N6iqotXryQdVVezycg4gp0IyQZckSZLU7e1qIt3Wcg4gp0Kzi7skSZIkSRlggi5JkiRJUgbYxV1dytgpc6mr37jD9LauG5IkSZKkrsAEXV1KXf1Glk8ZX+owJEmSJKng7OIuSZIkSVIGmKBLklSGIuJHEbEmIl5oNq1fRDwUEa/m/vYtZYySJJUbE3RJksrTHcCp202bDDySUjoEeCT3WpIkFYkJuiRJZSil9Dvgre0mnw7cmXt+J3BGUYOSJKnMOUicJEnaZkBKaXXu+evAgHwzRcQkYBJAdXV1kUJTZ2jt7ijQeIeUJyefWOSIJKm8maCrIPId4D2wS8UzqKqCmslzOjSf7VIdkVJKEZFaKZsGTAOora3NO4+6hrbujtKRfYokqbBM0FUQ+Q7wHtil4ulo0m27VDveiIiBKaXVETEQWFPqgCRJKidegy5Jkra5H7gk9/wS4L4SxiJJUtkxQZckqQxFxC+Ap4BDI2JVREwEpgAnR8SrwMdyryVJUpHYxV2SpDKUUjq/laKTihqIJElq4hl0SZIkSZIywARdkiRJkqQMsIu7Mqu1W7dJkiRJUndkgq7MauverJIkSZLU3djFXZIkSZKkDDBBlyRJkiQpA0zQJUmSJEnKABN0SZIkSZIywARdkiRJkqQMcBR3dZpBVRXUTJ6Td/qTk08sQURS19JWG5IkSVL3Y4KuTtNaEp4v4ZC0I3/IkiRJKi92cZckSZIkKQNM0CVJkiRJygATdEmSJEmSMsBr0CVJkrSD1gaq3FYmSSo8E3RJkiTtwIEqJan4TNBVdG39Ir/9fJKKY+yUudTVb2wxzVsiSpIkFZcJuorOL/xS9tTVb2T5lPEtpnlLREmSpOJykDhJkiRJkjLAM+iSJEldQL5LUbbxkhQpe2yz2hUm6JIkSV1AvktRtvGSFCl7bLPaFXZxlyRJkiQpA0zQJUmSJEnKALu4dyOdcZskb70kdS+t3eZwd25r6H5CkqTC8dr18maC3o10xm2SvPWS1L10xkHd/YQkSYXjtevlzS7ukiRJkiRlgAm6JEmSJEkZYIIuSZIkSVIGeA26JElSRrQ3OFRrWhsAsj27M0CkpLbbnu1Lu8IEXTst347IHZCk5toaLd7RZ6XWtTU4VFtsV1Jp2PZUaCbo2mnuiCS1p7X9hKPPSpIktc5r0CVJkiRJygATdEmSJEmSMsAu7pIkSUW0qwPBSZK6PxN0SZKkItrVgeAkSd2fXdwlSZIkScoAz6AXQb6ubLt7q6HW3nN7u3pf1LbeU1J56Oj+Y3f3E52xj5QkSeqKTNCLIF9Xtt291VBHu8f5BVfSrirW/qMz9pGSJEldkQm6JEkqa+0N2rYrP1Y5EJyktrTVS62tfcTu9o5t7T3z7ec6Y9+o9pmgS5KkstZWr7Rd/SLsQHCS2rKryW1nJMWt7ec6Y9+o9hVlkLiIODUiFkfEkoiYXIx1SpKkXeNxW5Kk0uj0BD0i9gS+D/y/wEjg/IgY2dnrlSRJO8/jtiRJpVOMM+jHAEtSSstSSu8BM4DTi7BeSZK08zxuS5JUIpFS6twVRJwNnJpSuiz3+tPA6JTSl5rNMwmYlHt5KLC4wGF8APhrgd+z1LpbnaxPtnW3+kD3q5P1ybbOqM8HU0r9C/yeHrezze2yI7dJfm6X/Nwu+blddlSy43YmBolLKU0DpnXW+0fEMyml2s56/1LobnWyPtnW3eoD3a9O1ifbult9PG6XhttlR26T/Nwu+bld8nO77KiU26QYXdzrgMHNXh+UmyZJkrLH47YkSSVSjAT9aeCQiBgSEXsD5wH3F2G9kiRp53ncliSpRDq9i3tK6f2I+BLwW2BP4EcppUWdvd7tdFo3vBLqbnWyPtnW3eoD3a9O1ifbukx9PG5nmttlR26T/Nwu+bld8nO77Khk26TTB4mTJEmSJEntK0YXd0mSJEmS1A4TdEmSJEmSMqDbJOgRMTgiHo2IFyNiUUR8Jc88+0XEryJiQW6eS0sRa0dERM+I+FOzWL+RZ559ImJmRCyJiHkRUVP8SDuug3W6Ivc/fD4iHomID5Yi1o7oSH2azfsPEZEiIrO3sOhofSLinGbt7K5ix9lRHfy8Vef2G8/lPnOfKEWsOyMi9szF+0Cesi61T9imnTp1mX3CNm3Vp9k8md8ndLaI6BcRD0XEq7m/fVuZ75LcPK9GxCXNpj8WEYsj4s+5xwHFi76wIuLUXF2WRMTkPOWttu2IuCY3fXFEnFLMuDvbrm6XiKiJiI3NPhu3FTv2ztSB7fLRiHg2It6PiLO3K8vbnrq63dwmW5p9VrrVgJgd2C6tHmO762cFdnu7dP7nJaXULR7AQODo3PPewCvAyO3muRb4bu55f+AtYO9Sx95KfQKozD3fC5gHHLvdPF8Abss9Pw+YWeq4C1CnE4B9c8//Mct16kh9cmW9gd8BfwRqSx33bv5/DgGeA/rmXh9Q6rh3sz7TgH/MPR8JLC913B2o1xXAXcADecq61D6hg3XqMvuEjtQnV94l9glF2E7fAybnnk/ednzebp5+wLLc376559v2P491h+1H40B8S4GDgb2BBXm+v+Rt27n91gJgH2BI7n32LHWdMrBdaoAXSl2HEm6XGuBDwE+As5tNb7U9deXH7myTXFlDqetQwu2S9xjbXT8ru7tdivV56TZn0FNKq1NKz+aebwBeAgZtPxvQOyICqKQxQX+/qIF2UGrUkHu5V+6x/Yh+pwN35p7PAk7K1S2TOlKnlNKjKaW/5V7+kcb772ZSB/9HAN8CvgtsKlZsu6KD9fks8P2U0tu5ZdYUMcSd0sH6JKBP7vl+wGtFCm+XRMRBwHjgh63M0qX2CdB+nbrSPgE69D+CLrJPKILmn9c7gTPyzHMK8FBK6a3cfuch4NQixVcsxwBLUkrLUkrvATNo3DbNtda2TwdmpJTeTSn9N7Ak937dwe5sl+6s3e2SUlqeUnoe2Lrdst21Pe3ONunOOrJdWjvGdtfPCuzedimKbpOgN5fr4nQUjWfMmrsV+Dsav4QvBL6SUspsQ811k/wzsIbGRrJ9fQYBK6HxtjjAOmD/4ka5czpQp+YmAr8pTmS7pr36RMTRwOCU0pySBLiTOvD/GQ4Mj4gnI+KPEZHpnXUH6nMDcFFErAJ+DXy5yCHurH8DrqL1Lxhdbp9A+3VqLvP7BNqpT1fbJ3SyASml1bnnrwMD8szT9JnOWUXLH99/nOtmeF0XTszaq2OLebZr2x1Ztqvane0CMCQaLzV5PCKO6+xgi2h3/ufd9fOyu/XqGRHP5L7X5PuhsKva2e3S/BjbXT8rsHvbBYrweel2CXpEVAL/BfxzSmn9dsWnAH8GDgSOBG6NiD5kVEppS0rpSBp/tTkmIkaVOqbd1dE6RcRFQC3wr8WMb2e1VZ+I2AO4CbiyVPHtrA78f3rQ2M19HHA+8J8RUVXcKDuuA/U5H7gjpXQQ8Angp7n/W+ZExCeBNSml+aWOpVB2pk5dYZ/QXn264j5hd0XEwxHxQp7H9mcrEvl7ILXlwpTS4cBxucenCxS2ur7VQHVK6Shyl5xk+fueSu6DKaVa4ALg3yJiaKkDKraucIwthVa2S6d/XjL5RXRXRcReNCbnP08pzc4zy6XA7FzX1yXAfwMjihnjrkgp1QOPsmPXkjpgMEBE9KCxi+7a4ka3a9qoExHxMeD/+//Zu+/wqMq8jePfxwRCIEICIj0bmiAERIhIdJFQLCuK4ILSFoIgoqAuNkDp+ipYEFikKE2wgLIgza6Aq2aRJpEoSBEIQQGpoRPyvH+cSTaBQAKZ5Mwk9+e65pqZU+8zM5mT35znPAdoY609ld/ZLscFtudKIBJYbozZDjQBFhk/6BTqIu/PLmCRtfaMpznlrzgFu0+7yPb0BD7wTBMHFAOuyt90OXYz0MbzWZoDtDDGvHPONP72nZCTbfKn74TstsdvvxMul7W2lbU2MovbQmCPMaYCgOc+q1Nm0j/THpU9w7DWpt0n45zz769Nuy+4jVlNc87fdk7m9VeX/bp4mvzvB/D8YLYVpwVYQZCb97ygfl5ytV0Zvku24fRtcb03w7koR6/LBfaxBfWzArl7XfLl81JgCnRP07ZpwC/W2jEXmGwn0NIzfTmgFk6nBz7HGFM27cikMSYYuBXYeM5ki4C0XhXbA197jkL4pJxskzHmemAKzh+Dz57fDNlvj7X2sLX2KmtthLU2AuccljbW2tWuBM5GDj9zH+EcPccYcxXOPzz+/DeU8TvhWpwCfV9+5swpa+0ga21lz2epI87fe9dzJvOr74ScbJM/fSdktz3+9p2QDzJ+XrsDC7OY5jPgNmNMmHF6eb8N+MwYE+j5Dkr7cf4uYEM+ZM4Lq4CaxpiqxpiiOJ+dc3sGvtDf9iKgo3F6M6+K84PpD/mUO69d9uvi+f4PADDGVMN5XXxyX3UZcvK6XEiWf095lDM/XfZr4nktgjyPr8L5ofXnPEuav7J9XS6yjy2onxXIxeuSb58X6wO96XnjBvwVp3lcPE4z9h9xmqz2Afp4pqkIfI5z/vkGoKvbuS+yPfVxesuO92Qd6hk+0vNhAaeY+BCnU5gfgGpu5/bCNn0J7MnwHi5yO3dutuec6Zfjwz0O5/D9MThNdH/2/B11dDt3LrenDvAdTg+ePwK3uZ07h9sWg6eHcH/+TsjhNvnNd0JOtuecaXz6OyEfXqMywFfAZs/7XNozPAqYmmG6Bzyf6S1AD8+wEsAaz993AjAOP+69HOf/lV9xjvQ+d+7n5mJ/2zhHeLYCm4C/ub0tvvC6AH/3fC5+BNYCd7u9Lfn8utyA0+LtGE5Li4QM857391QQbpf7mgA34fw/s95z39Ptbcnn1+WC+9iC+lnJzeuSX58X41mZiIiIiIiIiLiowDRxFxEREREREfFnKtBFREREREREfIAKdBEREREREREfoAJdRERERERExAeoQBcRERERERHxASrQRURERERERHyACnQRERERERERH6ACXURERERERMQHqEAXERERERER8QEq0EVERERERER8gAp0ERERERERER+gAl1ERERERETEB6hAF3GBMSbBGBPjdg43GWPaGWMSjTFHjTHXX8b8M40xL1xk/FFjTLXcpRQREdF+G3K/377EdV10Hy9SkKlAF/EyY8x2Y0yrc4bFGmO+TXtura1rrV2ezXIijDHWGBOYR1Hd9irQz1obYq1dd+5I43jMGLPBGHPMGLPLGPOhMaZeThbuWe42z7Ke9iwn2RjzmzHmaS9vi4iI+Cntt3Msu/229eyvjxpjkowxY4wxAS7kFPFrKtBFCikf+AfiL0DCRcaPAx4HHgNKA9cAHwGtL2NdBugGhAF3AP2MMR0vYzkiIiKu8IP9NsB11toQoCXQGXjw3Al8YDtEfJoKdBEXZPy13hjT2Biz2hhzxBizxxgzxjPZN577Q55fozJ09coAACAASURBVKONMVcYYwYbY3YYY/YaY2YZY0plWG43z7j9xpgh56xnuDFmnjHmHWPMESDWs+44Y8whY8zvxpgJxpiiGZZnjTGPGGM2e44+P2+MqW6M+d6T94OM05+zjVlmNcYEGWOOAgHAemPM1izmrQn0BTpZa7+21p6y1h631r5rrR2VYdIwY8xST7aVxpjq52SvAWCtfdlau9Zam2Kt3QQsBG6+1PdNREQKJ+23L77fPpe1diPwHyAyw+s3wBgTDxwzxgQaY641xiz3bEuCMabNOYu5yhjzhWc7Vhhj/pKDt0rE76lAF3HfOGCctbYkUB34wDP8Fs99qKc5WRwQ67k1B6oBIcAEAGNMHWAi0AWoAJQCKp2zrnuAeUAo8C5wFugPXAVE4/zi/cg589wONAKaAM8AbwJdgSo4O95OF9iuLLN6iu0QzzTXWWurZzFvS2CXtfaHCyw7TUdgBM6R8S3A/2UzPcYYAzQl+6MAIiIiWdF+OxuebWsKZGwK3wmnFVwoTsu2xcDnwNXAo8C7xphaGabvAjzv2dYfPdsvUuCpQBfJGx95fhE+ZIw5hLMDvpAzQA1jzFXW2qPW2v9eZNouwBhr7TZr7VFgENDROM3F2gOLrbXfWmtPA0MBe878cdbaj6y1qdbaE9baNdba/3qOLG8HpgDNzpnnZWvtEWttArAB+Nyz/sPAJ8CFOoq5WNbslAF+z8F0C6y1P1hrU3B23A1yMM9wnO++GTmYVkRECgftt3O3306z1hhzEKf4nkrmfe14a22itfYEzo8HIcAoa+1pa+3XwBIy/3iw1Fr7jbX2FPAcEG2MqXIJWUT8kgp0kbzR1lobmnbj/F+3M+qJc371RmPMKmPMXReZtiKwI8PzHUAgUM4zLjFthLX2OLD/nPkTMz4xxlxjjFlijPnD03zuRZxfqjPak+HxiSyeh5C1i2XNzn6cownZ+SPD4+MXyQKAMaYfzrnorT07fBEREdB+O7usOdXQWhtmra1urR1srU29wLZUBBLPGb+DzC0IMr42R4EDnvlECjQV6CIus9ZuttZ2wmniNRqYZ4wpwfm/ogPsxumkJU04kIKz8/0dqJw2whgTjHMkOtPqznk+CdgI1PQ01XsWp9mZN1wsa3a+AiobY6K8lAVjzAPAQKCltXaXt5YrIiKFi/bbly3jtuwGqhhjMtYi4UBShufpR8uNMSE4Hcbu9lIWEZ+lAl3EZcaYrsaYsp5fkQ95BqcC+zz3Ga/l/T7Q3xhT1bOzehGY62niPQ+42xhzk6cDmOFkv9O+EjgCHDXG1AYe9tZ2ZZP1oqy1m3GaF75vjIkxxhQ1xhQzxnQ0xgy81CDGmC6e9d+aduk1ERGRy6H9tlesxGn59owxpohxrjF/NzAnwzR3GmP+6nltngf+a61NPH9RIgWLCnQR990BJHh6SB0HdPScZ3Ycp9Oz7zznxDUBpgOzcXqK/Q04idOxCp5zzR7F2bn9DhwF9gIXa8r9FM5lUJKBt4C5XtyuC2bNocdwOtJ5A+cfoK1AO5zz2i7VCzhHJVYZp2fdo8aYyZexHBEREe23c8lzzv3dwN+AP3F+lO/m6f09zXvAMJym7Y1wOroTKfCMtVm1xhERf+f59fsQTjO439zOIyIiIhem/baIgI6gixQoxpi7jTHFPefCvQr8BGx3N5WIiIhkRfttETmXCnSRguUenA5UdgM1cZrdqZmMiIiIb9J+W0QyURN3ERERERERER+gI+giIiIiIiIiPiDQ7QDnuuqqq2xERITbMURERHzamjVr/rTWlnU7h/bbIiIi2cvpftvnCvSIiAhWr17tdgwRERGfZozZ4XYG0H5bREQkJ3K631YTdxEREREREREfoAJdRERERERExAeoQBcRERERERHxASrQRURERERERHyACnQRERERERERH6ACXURERERERMQHqEAXERERERER8QEq0EVERERERER8gAp0ERERERERER+gAl1ERERERETEB6hAFxEREREREfEBKtBFREREREREfECg2wHEx50+DQsXwpIlsHYt7N4NxkClStCwIbRu7dyCg91OKiIiIiIi4tdUoEvWzp6FqVNhxAj4/XcoUwZuugluuQWshZ07ncJ95kwoWxYefxz694fixd1OLiIiIpInbh71NUmHTnh1mZVCg/luYAuvLlNE/JcKdDlfUhLcfz989x00bQrTp8Ott0JAQObpUlJg+XJ4/XUYPBjeegsmTYK//c2V2CIFUZ8+fahUqRJDhgxxO4qISKGXdOgE20e19uoyIwYu9eryRMS/5egcdGPMHcaYTcaYLcaYgVmMv8UYs9YYk2KMaZ9heANjTJwxJsEYE2+Mud+b4SUP/Pe/cP318OOPMGsWrFgBd9xxfnEOEBgIrVrB0qXOdCVKwJ13wjPPwJkz+Z9dJI/ExsZStGhRQkJC0m/Dhg3L8fwxMTFMnTr1stY9efJkvyjOX3jhBerWrcsVV1zBzJkzM43btWsXMTExFC9enIYNG7Jhw4b0cWfOnKFnz55ceeWVhIeH88EHH2Sad/z48ZQrV46wsDAGDRqUH5siIiIi4ppsC3RjTADwBvA3oA7QyRhT55zJdgKxwHvnDD8OdLPW1gXuAMYaY0JzG1ryyH/+4xwpL1kSVq2Cf/zDOd88J265BdasgT594JVX4O674dixvM0rko+eeeYZjh49mn4bMWKE25F8SrVq1Rg/fjyNGjU6b1zv3r2JjIxk//793H///dx///9+q3399df56aefSExMZNasWTzwwAMkJiYCsHLlSoYPH87XX3/Nhg0bmDNnznkFvIiIiEhBkpMj6I2BLdbabdba08Ac4J6ME1hrt1tr44HUc4b/aq3d7Hm8G9gLlPVKcvGudeucpumVKsE338C11176MooVc5q4v/UWfPEFtGwJ+/d7P6uID5k5cybR0dHExsZy5ZVXcvPNN/Pnn38C8OKLLxISEsJ//vMf+vXrR0hIyHkFbEREBBMmTKBRo0aUKFGCNm3aALBkyRJCQkIoUqQIgwcPPm+98+fPJzIykrCwMO68805+//339HFLliyhVq1aXHnlldSsWZPPPvssD18BR+fOnWnZsiVBQUGZhh85coTPP/+cgQMHEhwcTP/+/dmxYwfx8fEAfPjhhzz22GOEhoYSExNDdHQ0CxYsAGDevHnce++91K1bl0qVKtGrVy/mzJmT59siIiIi4pacFOiVgMQMz3d5hl0SY0xjoCiwNYtxvY0xq40xq/ft23epi5bc+v13aNMGwsLg66+hYsXcLa9XL/j3v51m8rffDkeOeCeniI9atWoVXbp0Ye/evaSkpDB9+nQAnn32WY4ePUrTpk2ZMGECR48eZc2aNefNP2XKFGbNmsXhw4d59tlnAbjrrrs4evQoXbp0OW/6H374gZ49ezJt2jT27dvH9ddfT+/evdPHP/jgg4wcOZLk5GQ+//xzKlW65K9sr9myZQvFihUjJCSEpk2bsm3bNqpXr86mTZsA2LRpE7Vr16Zr167MmTOHOnXqnDdu3LhxPPnkk5nGiYiIiBRE+XIddGNMBWA20MNam3rueGvtm9baKGttVNmyOsCer1JSoEMHOHAAFi/OfXGepm1bmD8f1q93mruf8G6PpyL57dVXXyU0NDT9Nn/+/PRx11xzDbfeeivBwcG0bNmSX3/99ZKW3bt3b+rWrUtgYCBNmjTJdvpp06YRGxvLjTfeSGBgIE899RRLly7l1KlTAFxxxRVs3bqVI0eOULVqVSIjIy9tY73o2LFjhISEkJycTEJCAgcPHuTKK6/k6NGjmcb/9NNPJCUlZTlu69atbN68OdM4ERERkYIoJwV6ElAlw/PKnmE5YowpCSwFnrPW/vfS4kmee+klp7f2t96CBg28u+w774TZs51z2//xD+fybCJ+6qmnnuLQoUPpt3vvvTd9XOnSpdMfFy1alJMnT17SsmvWrHlJ0ycmJjJlypT0HwuqVq1K0aJF05u5f/jhh8TFxREeHk6TJk0ydcqW30qUKMHRo0epUqUKBw4cIDo6muTkZEJCQjKNX79+PU8++WSW48aPH8+iRYsyjRMREREpiHJSoK8CahpjqhpjigIdgUU5Wbhn+gXALGvtvMuPKXli1SrnOuedOzu3vNCxI7z8stPk/f/+L2/WIeLjrrji4l+1gYGXdsXLKlWqMGTIkEw/GJw8eZKIiAgAbrrpJhYvXsyePXuoXbt2luew55caNWpw4sQJdu3aBcDp06fZunUrtWrVApzWBxs3bkyf/ueff87ROBEREZGCKNsC3VqbAvQDPgN+AT6w1iYYY0YaY9oAGGNuMMbsAjoAU4wxCZ7Z7wNuAWKNMT96bl4+TCuXJSUFHnoIypWDN97I23U9+SR07QpDhsDChXm7LhEfVL58ea8exY6NjWXy5MmsXbsWay179+5l7ty5AKSmpvLOO++QnJyc/sNAyZIlvbbuCzlz5gwnT57EWpv+ODU1lZIlS3L77bczatQoTp48yeuvv054eDj16tUD4L777mP8+PEcOnSIFStWEBcXR7t27QDo0KED8+fPJyEhgaSkJKZNm5apB3gRERGRgiZH56Bbaz+21l5jra1urf0/z7Ch1tpFnserrLWVrbUlrLVlPJdVw1r7jrW2iLW2QYbbj3m3OZJjEyc6PbePGweheXzlO2PgzTfhhhugWzf47be8XZ9IHnj55ZczXQe9V69eOZ73ySef5IsvvqBixYrExMTkaJ7bbruNkJAQ3n333fR1x8bGAhAdHc2YMWPo0aMHpUqVIioqKlPnc++88w5VqlThqquuYvfu3fxfPrReefDBBwkODub777+nd+/eBAcH88033wBOJ3g//fQTYWFhzJkzh7lz52I8l3Ds378/kZGRVKlSha5duzJt2jSqVHHOqrrxxhsZNmwYzZs3JzIykvvuu08FuoiIiBRoxvrYecFRUVF29erVbsco2PbuhZo1IToaPvkk59c6z63t253z3OvUcS7ldonNekVE5H+MMWustVFu59B+WwqTiIFL2T6qtc8vU0R8T0732/nSi7v4mBdfhGPHYOzY/CvOASIiYMoUiItzzn0XERERERGRdCrQC5sdO2DSJIiNhdq183/9998PPXo4HcZ9913+r19ERERERMRHqUAvbEaMcI6aDxvmXobx4yE8HHr1gku8HJWIiIiIiEhBpQK9MNmyBd5+G/r2hSpVsp8+r4SEOE3dN26E5593L4eIiIiIiIgPUYFemLz2GhQpAk8/7XYSuP126N4dRo92epMXEREREREp5FSgFxZ79sCMGc5lzsqXdzuNY8wYuOoq6NkTzp51O41IlpYvX44xhgEDBqQPa9y4cfplwgqLPn36ZLrMXFBQUPq1zLOzcOFCoqOjCQoKSr9UXEZjx44lPDyckJAQGjZsSGpqKgCnT5+mb9++lCtXjjJlyvDMM8/ga1ceEREREfEmFeiFxYQJcPo0PPmk20n+p3Rp5zrs69Y510kX8VFly5blyy+/BGDnzp0cPnzY5UT5b/LkyRw9ejT91r59e/7+97/naN5SpUrx9NNP07Nnz/PGzZkzh9dee40FCxaQnJzMzJkz03/8mDRpEitXrmTjxo2sW7eOuXPn8u6773p1u0RERER8iS5EXRgcOwZvvAFt20KtWm6nAeDmUV+TdOgE2BK8H16P2k8MoPnGUA4Fl8y3DJVCg/luYIt8W5/4r6JFi1KnTh3WrVvHsmXLuPfeexk1alT6+Pnz5zN06FCSkpKIjo5m2rRpVKhQAYCOHTuybNkyjh49Sv369ZkyZQr169cHICIigu7du/Puu+9y4MABXnzxRfr06ePKNl6Kw4cPs2DBAjZs2JCj6WNiYgBYu3Ytx48fzzRu8uTJDBo0iEaNGgGkvzYAK1asoHPnzoSFhREWFkanTp2YP38+Xbt29c6GFHLGmFBgKhAJWOABYBMwF4gAtgP3WWsPuhRRRESk0FGBXhi8/z4cPAhPPOF2knRJh06wfVRr58k/IqBBA348tQLGTc63DBEDl+bbusT/tW/fnnnz5vGf//yHl19+Ob1A/+GHH+jZsyeffvopjRo1YtiwYfTu3ZvFixcDcP311zN+/HjKlCnDs88+S/fu3VmXod+FlStXEh8fz2effUb37t3p1asXgYG+/dU8Z84cGjZsSLVq1XK9rPXr1xMTE0P16tU5c+YMPXr0YMSIEQBZNmffsmVLrtcp6cYBn1pr2xtjigLFgWeBr6y1o4wxA4GBwICLLURERES8R03cC4PJkyEyEm6+2e0kWYuMhEcfdZq5r1njdhqRLN1+++38+9//JigoiKuuuip9+LRp04iNjeXGG28kMDCQp556iqVLl3Lq1CkABgwYwNVXX01AQACxsbHEx8dnWu4DDzxA8eLFad26NcnJyfzxxx/5ul2XY8aMGXTr1s0ryzpy5Agff/wxcXFxfPvtt7z99tssWLAAcI68v/fee+zfv5/t27fz4YcfnncEXi6PMaYUcAswDcBae9paewi4B3jbM9nbQFt3EoqIiBROvn2YRnJv9Wqn6J0wAYz5X9PyXMiTpuHDh8N77zlH+Zcvd67VLuJDihUrRps2bWjQoEGm4YmJiSxfvpwZM2akDytatCi///47VapUYfDgwXzwwQfs37+f1NRUUlNTOXv2LAEBAQCULl06fR6AkydP5tMWXZ5ffvmF9evXc99993llecWLF6dHjx5cffXVANx7770sX76cdu3a8fDDD7N582bq1atHWFgYbdu2ZdmyZV5Zr1AV2AfMMMZcB6wBHgfKWWt/90zzB1DOpXwiIiKFkgr0gm7yZCheHDznbGZqWn6Z8qRpeKlSTpH+yCOwZAncfbf31yGSSy+//DKQuZl1lSpVGDJkCIMGDTpv+tmzZ7NgwQKWLVtGeHg48fHxXHfddX7dE/mMGTO4++67CQ0N9cryqlevfl6P+GmvT9GiRZkwYQITJkwA4KmnnuK6667zynqFQKAh8Ki1dqUxZhxOc/Z01lprjMnyw2qM6Q30BggPD8/rrCIiIoWGmrgXZIcPO+efd+7sFMC+rlcvuOYaeOYZSElxO41IjsTGxjJ58mTWrl2LtZa9e/cyd+5cAJKTkwkODiYsLIzk5GRefPFFl9PmztmzZ5k9e/YFm7fPnz+fGjVqkJSUdN58J0+e5OzZs+mPUzx/4+3atWPq1Kns37+f3bt3s2DBApo3bw7AsWPH2LJlC6mpqXz55ZdMnTqV3r175+1GFh67gF3W2pWe5/NwCvY9xpgKAJ77vVnNbK1901obZa2NKlu2bL4EFhERKQxUoBdk8+bB8eNO4esPihSB0aNh40aYPt3tNCI5Eh0dzZgxY+jRowelSpUiKiqKNZ6+FLp160Z4eDgVK1bkuuuuo0mTJi6nzZ1PPvmEs2fPcscdd2Q5/siRI2zdupUzZ85kGj579myCg4MZNWoU77zzDsHBwbzwwguAc45+ZGQk1apVo1GjRnTr1o127doBcPz4cVq3bk1ISAh9+vRh4sSJREdH5+1GFhLW2j+ARGNM2qU9WgI/A4uA7p5h3YGFLsQTEREptIyvNbWMioqyq1evdjtGwRATA7t3w6ZN6ed0Rwxc6pUm7nm2DGuhaVPYuhU2b4aQkFyt57IyiIj4AWPMGmttVC7mb4BzmbWiwDagB84P9x8A4cAOnMusHbjYcrTflsIkL/530P8jIoVDTvfbOoJeUO3YAStWQLdu/tXhmjHwyivwxx8wZozbaURECixr7Y+eZur1rbVtrbUHrbX7rbUtrbU1rbWtsivORURExLtUoBdU777r3Hs6h/Mr0dHQtq1ToB886HYaERERERGRfKECvSCyFmbPdpqKR0S4nebyjBjhdHL32mtuJxEREREREckXKtALoh9/dDpa88ej52nq14f77oNx4+DPP91OIyIiIiIikudUoBdE8+ZBQADce6/bSXJn2DA4dsw5J11ERERERKSAC3Q7gHiZtU6BHhMDV13ldprcqVPHuYb7hAnwxBNQrpzbiaSQioiIYM+ePQQEBFCuXDn69+9Pv379ADDG0LhxY1audC4n/cwzz/DKK6+wbNkyYmJiWLRoEc899xzbt2+ndOnS9O7dm+eee87Nzbksy5Yt47HHHmPHjh2UKFGCTp068corrxAQEHDR+Xbu3EmdOnUyDTt27Bjz5s3j73//O0lJSTz44IN89913hIaG8tJLL9G5c2cATp8+Tf/+/Zk3bx4pKSn07NmT0aNHY/yp40sRcc3No74m6dAJry6zUmiwV5cnInIuFegFTUIC/Por/POfbifxjqFD4f33neujq1d3cdHixYtp1aoVP/zwA82bN6devXo0a9YMgEOHDpGYmEiVKlX46quvKFu2LADbtm2jW7duLFq0iKZNm5KYmEhcXJybm3HZateuzccff0zlypU5evQo7du354033uCxxx676Hzh4eEcPXo0/fnmzZtp2LBh+rXU+/XrR5UqVfjzzz/59ttv+dvf/kZ0dDRVq1Zl0qRJrFy5ko0bN5KcnEzTpk2pX78+Xf359B0RyTdJh07o8mUi4nfUxL2g+fe/nUuVtWvndhLvuOYa51JxkybBnj1upxGhcePG1K1bl1WrVqUPu/fee5k3bx5r166lTp06FC1aFIC1a9dStWpVbrnlFowxhIeHc//997sVPVcqVKhAlSpVMMaQkpJCamoqv/zyyyUvZ8aMGbRr144SJUoAsGLFCh566CGKFClC8+bNadCgAYsXL04f17lzZ8LCwggPD6dTp07Mnz/fq9slIiIi4ktUoBc08+Y5vbeXL+92Eu8ZNAhOn4bXX3c7iQhr1qwhISGBa6+9Nn1Y27ZtWbhwYXqz7TT16tUjISGBYcOGsX79eqy1bkT2mp07d1KqVClKly7NypUriY2NvaT5U1NTmTVrFt26dUsfltVrsmXLlmzHiYiIiBREauJekGzaBBs2OD2f56FKocFEDFya62Xk2DXXQIcOMHEiDBgAYWG5WrfI5Wjbti2pqakYYxg1ahStW/+v2WSZMmUICgpi3rx5DBkyJH14rVq1+Pjjjxk9ejSjR4+mfPnyjBs3jnvuuceNTci18PBwDh8+zLZt25g1axY1a9a8pPk///xzjDG0aNEifVhMTAxTpkxh/PjxfPvtt6xbt47IyMj0cbNnz6Z79+4kJyfz4Ycf6vzzQiKvzh3+bmCL7CeUHNF7JCKSN1SgFyRpTT/zuPd2V3aegwbB3LlOh3EZCiCR/PLRRx/RvHlzhg4dyvfff8+jjz6aaXyPHj1YvXo1wcGZf3xq1aoVrVq14sSJE8yYMYNOnTqRmJhImTJl8jO+V1WrVo3atWvTt29f3n///RzPN2PGDLp06cIVV/yv8da//vUvHnroISpXrkyDBg1o06YNJUuWBODhhx9m8+bN1KtXj7CwMNq2bcuyZcu8vj3ie/Li3OHc/rAsmek9EhHJG2riXpAsWQKNGkHlym4n8b7rroO77oKxYyFDh1Mi+SkgIIDhw4cTFxfH8uXLM43r2LEjr7766gXnDQ4O5pFHHqFYsWJs27Ytj5PmPWstP/30U46nP3jwIAsXLszUvB2gcuXKLF26lH379vHFF1/w22+/cd111wFQtGhRJkyYwO7du0lISMAYkz5OREREpCBSgV5Q7N8P//0vtC7AvZU++ywcOABvvul2EinEihQpwqOPPsrIkSOznXb9+vVMmzaNQ4cOkZqayvvvv8+ZM2e45ppr8iGpd73//vvEx8eTmprKzp07efXVV4mOjs40zfz586lRowZJSUnnzf/ee+8RGRl53iXXEhMTOXDgAKdOnWLs2LEkJibSztPJ5bFjx9iyZQupqal8+eWXTJ06ld69e+fdRoqIiIi4TAV6QfHpp5CaWrAL9OhoaN4cXn0VTp1yO40UYj179mTVqlXZXjItJCSEDz74gOrVqxMaGsprr73G/PnzKVWqVD4l9Z79+/fTvn17SpYsyQ033ED9+vXPazFw5MgRtm7dypkzZ86bf8aMGecdPQfYtGkTkZGRhIWFMXfuXD755JP0Ju7Hjx+ndevWhISE0KdPHyZOnHjejwIiIiIiBYnOQS8oliyBq6+GqCi3k+StZ5+FW2+FmTPhoYfcTiOFxPbt2zM9Dw0NJTk5Gci6p/Fdu3alP/7ss8/yNFt+6devH/369bvoNLGxsRfs2X316tVZDm/VqhW7d+/OclzZsmXZtGnTJeUUERER8Wc6gl4QpKQ4R9DvvBOuKOBvacuW0LgxjB7tbLeIiIiIiEgBUcCruUIiLg4OHSrYzdvTGAMDB8Jvv8GCBW6nERERERER8RoV6AXBkiUQGOg0/S4M2rSBGjWcc9GzaF4sIiIiIiLij3QOekHwySfQtCn4YcdTlyUgAJ54Ah55BL77Dv76V7cTiYiIiIjk2M2jvibp0AmvLrNSaDDfDWzh1WVK/lOB7u/++AN++glGjXI7Sf7q3h2GDHGOoqtAFxERERE/knToBNtHeff01IiBS726PHGHmrj7uy+/dO4LS/P2NMWLO0fQFy2CX391O42IiIiIiEiu6Qi6nzm3OcyrS2fSstiVNJybhP3g92znrxQanJfx8lffvvDyy/D66zBpkttpREREREREckUFup/J1BzGWnjnIbj7Dn4bfbe7wdxQrhx06+ZcE33kSChb1u1EIiIiIiIily1HBbox5g5gHBAATLXWjjpn/C3AWKA+0NFaOy/DuO7AYM/TF6y1b3sjuAAbN0JSUuFr3p7RE0/AW2/BxIkwbJjbaUREREQuSaXQYK+fO6zOwkT8V7YFujEmAHgDuBXYBawyxiyy1v6cYbKdQCzw1DnzlgaGAVGABdZ45j3onfiFXNr5561auZvDTbVrw113wYQJ8MwzEFyAmvCLiIhIgZcXhbQ6CxPxXznpJK4xsMVau81aexqYA9yTcQJrbiBQVgAAIABJREFU7XZrbTyQes68twNfWGsPeIryL4A7vJBbAL74AqpXh6pV3U7irqeegj//hFmz3E4iIiIiIiJy2XJSoFcCEjM83+UZlhM5mtcY09sYs9oYs3rfvn05XHQhd+YMLF9euI+ep7nlFmjUCMaOdc7LFxERERER8UM+cZk1a+2b1tooa21UWXX0lTNr1kByMrRs6XYS9xkDjz/unJP/xRdupxEREREREbksOSnQk4AqGZ5X9gzLidzMKxezYoVz36yZuzl8xX33Ob26jxvndhIREREREZHLkpMCfRVQ0xhT1RhTFOgILMrh8j8DbjPGhBljwoDbPMMkt775Bq69Fq6+2u0kviEoCB5+GD7+GDZvdjuNiIiIiIjIJcu2QLfWpgD9cArrX4APrLUJxpiRxpg2AMaYG4wxu4AOwBRjTIJn3gPA8zhF/ipgpGeY5MbZs/Dtt8651/I/Dz0ERYrAv/7ldhIREREREZFLlqProFtrPwY+PmfY0AyPV+E0X89q3unA9FxklHOtXw9HjqhAP1f58tCxI8yYAc8/D6VKuZ1IREREREQkx3JUoIuP+eYb514F+vkefxxmz3aK9H/+0+00IiLiIyqFBnv92tCVQoPz5BrWIiJSeKlA90crVkC1alA5y0YLhVujRnDzzU4z90cfhYAAtxOJiIgPyItC2tsFv4iIiE9cZk1yzthU+M9/1Hv7xTz2GGzb5nQYJyIiIiIi4idUoPuZmn/uhP371bz9Ytq1c1oX6JJrIiIiIiLiR1Sg+5nGiQnOAx1Bv7AiRaBvX/jqK0hIcDuNiIiIiIhIjugcdD/TJHGDc3Q4IsLtKL7twQdhxAgYPx6mTHE7jYiIiGtuHvU1SYdOeHWZlUKDvbo8ERFxqED3J9YStSsB7r4djHE7jW8rUwa6dnV6dH/pJShd2u1EIiI+xRizHUgGzgIp1tooY0xpYC4QAWwH7rPWHnQro3hH0qETbB/V2u0YIiKSA2ri7k927qT80QNOL+WSvccegxMnYNo0t5OIiPiq5tbaBtbaKM/zgcBX1tqawFee5yIiIpJPVKD7k7g45z462t0c/qJePedc/UmT4OxZt9OIiPiDe4C3PY/fBtq6mEVERKTQURN3f/L99xwvEkTx+vXdTuI/+vaF++6DTz6Bu+5yO42IiC+xwOfGGAtMsda+CZSz1v7uGf8HUC6rGY0xvYHeAOHh4fmR1SdVCg32+rXQK4UG58k120VExD+oQPcncXGsr3AN0YF623KsbVuoWBHeeEMFuohIZn+11iYZY64GvjDGbMw40lprPcX7eTzF/JsAUVFRWU5TGORFIe3tgl9ERPyLmrj7i+PH4ccfWVuxtttJ/EuRIvDQQ/Dpp7B5s9tpRER8hrU2yXO/F1gANAb2GGMqAHju97qXUEREpPBRge4vVq+GlBTWVlKBfskefBACA51z0UVEBGNMCWPMlWmPgduADcAioLtnsu7AQncSioiIFE4q0P2Fp4M4HUG/DBUqQPv2MGMGHDvmdhoREV9QDvjWGLMe+AFYaq39FBgF3GqM2Qy08jwXERGRfKKTmf1FXBzUrMnB4qXcTuKf+vaFOXPgvfecI+oiIoWYtXYbcF0Ww/cDLfM/kYgURjeP+pqkQyfcjpEtdd7o2/Lic+Tme64C3R9YC99/D3/7m9tJ/NfNN0P9+k5ncb16gTFuJxIREREp1JIOnWD7qNZux8iWOm/0bXnxOXLzPVcTd3+wbRvs2wc33eR2Ev9lDPTrB+vXOz92iIiIiIiI+BgV6P4graCMjnY3h7/r3BlKlYIJE9xOIiIiIiIich41cfcHK1dCSAjUrQskup3Gf5UoAT16OAX6mDFupxEREZEMKoUGe7VZaaXQYK8tS0S8I6/OFy9IVKD7gx9+gEaNICDA7ST+75FHYOxYeOstoJHbaURERMRDnXCJFHz+0u+Am9TE3dedPu2cN33DDW4nKRhq1oTbb4cpUwg8m+J2GhERERERkXQq0H1dfLxTpKtA956+fWH3bm7d/F+3k4iIiIiIiKRTge7rVq1y7lWge8+dd0JEBN3XLnE7iYiIiIiISDqdg+7rVq2CMmUgIsLtJAVHQAA8/DBNBgyADRsgMtLtRCIiIiIiPkUdurlDBbqvW7XKOXpujNtJCpYHHuDUs4MJmjgRJk50O42IiIiIiE9Rh27uUBN3X3bsGPz8s5q354WrrmLRtc1g1iw4fNjtNCIiIiIiIjqC7tPWrYPUVBXoeeTTZvfSYcOXDGs/gLcb3X1Zy6gUGqzLwoiIiIiIiFeoQPdl6iAuT02b8DCsfpsRicsZ8fmkyzqNIGLg0jxIJiIiIiIihZGauPuyVaugcmUoX97tJAVX376waRN89ZXbSUREREREpJBTge7L0jqIk7zToQOULQsTJridRERERERECjkV6L7q4EHYskUFel4rVgx69YLFi2HHDrfTiIiIiIhIIaYC3VetXu3cq0DPe336OPeTJ7ubQ0RERERECjUV6L4qrYO4qCh3cxQG4eHQpg1MnQonT7qdRkRERERECikV6L5q7VqoXh1CQ91OUjj06wd//gkffOB2EhERERERKaRUoPuqdevg+uvdTlF4tGgBtWvDG2+4nURERERERAopFei+6NAh2LZNBXp+Msa55NoPPzg3ERERERGRfBbodgDJwo8/Ovcq0PNXt24waJBzFL1xY7fTXJaZM2fSo0cPSpUqxW+//UZYWFj6uJSUFIoUKcKwYcMYPnz4JS83NTWVBx54wMuJ/2f48OGMGDEi2+lmzJhBTEwMVatW5a233qJXr14XnT4iIoKYmBhmzpyZq3zGmMt67fLS9u3bmTlzJt26daNatWp5uq7U1FSeeOIJ5s6dy549e2jTpg1jx46latWqzJgxg9jYWODin5W0z6e1Nk+zivizSqHBRAxc6vVlikjB5+3vD313uEMFui9at865b9jQ3RyFTcmSTpE+bRq8+qpzfXQ/dfjwYUaPHs2oUaO8sryZM2eSkpKSpwV6r169uOOOO9KfL126lBdeeIEPP/yQypUrpw+vXr06x44dy/FyFyxYQMmSJXOdLy4uLlMOX7B9+3ZGjBjBX//61zwv0OfNm8e4ceN47bXXiI6OpkyZMlSoUIG4uDiqV6+ePl1+fFZECrLvBrZwO4KI+Cl9fxQMKtB90bp1UKEClCvndpLCp29fmDjRKdIHDnQ7zWW77bbb+Ne//kX//v0p5yefo8qVK2cqgDdu3AhAgwYNqFGjRqZpL6VAv95LLVGaNGnileX4q19++QWAf/7zn1xxxf/Ojirsr4uIiIiIN+kcdF+kDuLcU6cONG8OkybB2bNup7lsgwcPBuCFF17IdtoffviBVq1aERISQokSJWjZsiU/ZDgPPyYmhhUrVvDdd99hjMEYQ0xMTPr43377jS5dulC2bFmCgoJo0KABCxYs8Po2ZeXs2bMMHTqUChUqEBoayt13382uXbsyTRMREZHe/Brgjz/+oHv37lSsWJGgoCAqVKjAXXfdxd69ey+6LmNMpubtv/76K+3atePqq6+mWLFihIeH06FDB1JSUi64jO3bt2OMYcqUKdnmPnPmDIMHDyYiIoKiRYsSERHB4MGDOXPmDADLly+nefPmANx6663p783y5csvuP7PP/+cO++8kwoVKlC8eHEiIyN57bXXOJvNZz0iIiJ92wMCAjDGMHPmzPTtSTt9ILvPioiIiIhcXI4KdGPMHcaYTcaYLcaY8w4rGmOCjDFzPeNXGmMiPMOLGGPeNsb8ZIz5xRgzyLvxC6ATJ+CXX1Sgu6lfP9i5E5YscTvJZatQoQL9+vXjzTffZMeOHRecLj4+nmbNmnHw4EFmzpzJrFmzOHLkCM2aNWP9+vUATJw4keuvv5769esTFxdHXFwcEydOBCAxMZEbb7yR9evX8/rrr7No0SIaNmzI3//+dxYtWpS+nrRCztvnb7/00kts2bKF6dOnM27cOOLi4ujatetF5/nHP/5BXFwcr7zyCl988QXjx4+ncuXKHD9+/JLW3bp1a5KSkpg0aRKfffYZo0aNIigoiNTUVK/k7t69O6NGjaJbt24sWbKE2NhYRo8eTffu3QFo2LAhb3iuOjB+/Pj096bhRU6N2bZtGy1btmT69OksXbqU7t27M3z4cJ577rmL5l2wYEH6jxxp62nduvV5013ssyIiIiIi2cu2ibsxJgB4A7gV2AWsMsYsstb+nGGynsBBa20NY0xHYDRwP9ABCLLW1jPGFAd+Nsa8b63d7u0NKTA2bHCO3KpAd0+bNlC5stNZ3D33uJ3msg0YMIApU6YwYsQIpk+fnuU0I0eOJCgoiK+++orQ0FDAORobERHBiBEjmD9/PnXq1KFkyZKkpKSc15x5+PDhWGtZsWIFZcqUAeD2228nMTGRoUOH0qZNG8A5+hwQEJCpabQ3RERE8N5776U/37dvH08//TS7d++mYsWKWc4TFxfHiy++SJcuXdKHdejQ4ZLW++eff7JlyxYWLlyYvo0AnTt39kruDRs28P7772fqlO62224jMDCQIUOGMHDgQOrXr0+dOnUAuPbaa3PU1LxPnz7pj621NG3alNOnT/Pqq6/y4osvXvD9uf7666lUqRKQuUn7uacaXOyzIiIi+ScvOhvMC+qETOR8OTkHvTGwxVq7DcAYMwe4B8hYoN8DDPc8ngdMMMYYwAIljDGBQDBwGjjinegFlDqIc19gIPTpA4MHw8aNzvXR/VDp0qV58sknGTFiBAMGDMjUkVeab775hrvuuiu9OAcoWbIkbdq0YfHixdmu49NPP+XOO++kVKlSmZp233777Tz99NMcOXKEkiVL8pe//OWiTb8v15133pnpeb169QDYuXPnBQv0G264gVdeeQVrLS1atCAyMhLn6yrnypQpQ7Vq1Rg4cCB79uwhJiaGmjVrei33N998A3DeUfWuXbsyZMgQVqxYQf369S8pM8Dvv//O8OHD+fTTT9m9e3em92Tv3r2UL1/+kpcpIiK+R52FifivnBzOqgQkZni+yzMsy2mstSnAYaAMTrF+DPgd2Am8aq09cO4KjDG9jTGrjTGr9+3bd8kbUaCsWwehoRAR4XaSwu3BB6FoUafDOD/Wv39/SpcuzdChQ7Mcf+DAASpUqHDe8PLly3Pw4MFsl793715mzZpFkSJFMt2efvppAPbv35+7DchG6dKlMz0PCgoC4OTJkxecZ+7cubRp04aXX36Z+vXrU6lSJUaOHJmjpulpjDF88cUXREVFMWjQIK655hqqVavGpEmTvJL7wAHna/Lc9yatgE4bfylSU1Np06YNS5YsYfDgwXz99desWrUqvXn7xV4zEREREckfed1JXGPgLFARqAo8aYw571pA1to3rbVR1tqosn58aSuvWLsWGjSASzyiJ1529dXQoQO8/TYkJ7ud5rKFhIQwaNAgPvzwQ3788cfzxpcuXZo//vjjvOF//PFHpmuoX0iZMmVo3749q1atyvJ2oaPYbrr66qt54403SEpKYuPGjcTGxjJs2DCmTJlyScupVq0as2bNYt++faxbt44WLVrwyCOP8Mknn+Q6Y1oBf+57k/b83AI/J7Zu3crq1asZPXo0Dz74IE2bNiUqKoqAgIBc582p2NhYXQNdRERE5CJyUqAnAVUyPK/sGZblNJ7m7KWA/UBn4FNr7Rlr7V7gOyAqt6ELrJQUiI/X+ee+ol8/OHIE3nnH7SS58sgjj1CpUqX0nt0zatasGR9//DHJGX6ESE5OZvHixZl63w4KCuLEiRPnzX/HHXcQHx9P3bp1iYqKOu+WdmTYV9WqVYsXX3yRsLAwNmzYcFnLMMbQoEEDxowZA3DZy8nolltuAWDOnDmZhr/77rsA6e9N2uub1XtzrrRO8IoUKZI+7MyZM+nL9JYLfVZEREREJHs5OQd9FVDTGFMVpxDviFN4Z7QI6A7EAe2Br6211hizE2gBzDbGlACaAGO9Fb7A2bQJTp5Uge4rbrzR6QvgjTecc9L9tFVDUFAQQ4cOpXfv3ueNGzJkCEuWLKFly5YMGDAAYwyjR4/m+PHjmZrF16lTh4kTJzJ37lyqV6/OlVdeSa1atRg5ciSNGzfmlltuoV+/fkRERHDw4EE2bNjAtm3b0jun27FjB9WrV2fo0KEXbG6fHw4fPkyrVq3o0qULtWvXpkiRIixcuJCDBw9y22235Xg58fHxPP7449x///3UqFGDs2fPMnPmTAIDA2nRIvfn/UVGRtKpUyeGDx9OSkoKN910E3FxcTz//PN06tQp/Zz1a665hsDAQKZPn07p0qUJCgqiVq1aXHnllect89prr+Uvf/kLzz33HAEBARQpUoTXX38911nPdaHPCjidEo4cOZJTp07l65F7EREREX+RbYFurU0xxvQDPgMCgOnW2gRjzEhgtbV2ETANpwjfAhzAKeLB6f19hjEmATDADGttfF5sSIGQ1kGcCnTfYIxzFP2BB2DFCvDj6zn36NGDV155hc2bN2caXr9+fZYvX85zzz1H9+7dsdbSpEkTVqxYwXXXXZc+3YABA9i0aRO9evXi6NGjNGvWjOXLlxMeHs7q1asZPnw4zz77LPv27aNMmTJERkamXw4MnB7Dz549e0nneeeFYsWK0bBhQ9566y127NjBFVdcQa1atXj33Xe55xJ67C9fvjzh4eGMGTOGXbt2UaxYMerVq8eSJUto1KiRV7LOnDmTatWqMX36dF544QUqVqzIgAEDGDZsWPo0ZcqUYcKECYwePZpmzZpx9uxZli1bluW1x4sWLcpHH31Ev3796NatG6VLl+aBBx4gPDycBx980CuZ4cKfFXDOgz979qyauYuIiIhcgPG1f5SioqLs6tWr3Y7hjiefdDolS052ehLPQsTApWwfdf71hyWPnDjhXHKtRQv48MPzRuv9EBG3GGPWWGtdP20sL/bb+m4VEV+SF99J+p7zbXnx/uR0v53XncTJpVi3DurXv2BxLi4IDoaePWHBAti1y+00IiIiIiJSgKlA9xXWwo8/Oj24i295+GFITYVL7OVbRERERETkUqhA9xVJSXDwIGQ471d8RNWq0Lo1vPkmnDrldhoREa8xxgQYY9YZY5Z4nlc1xqw0xmwxxsw1xhR1O6OIiEhhogLdV8R7+s6rX9/dHJK1fv1g717497/dTiIi4k2PA79keD4aeN1aWwM4CPR0JZWIiEghpQLdV6QV6J7LJ4mPufVWqFHDueSaiEgBYIypDLQGpnqeG5xLo87zTPI20NaddCIiIoWTCnRfER8Pf/kLlCrldhLJyhVXQN++8P33/7scnoiIfxsLPAOkXf+wDHDIWpvieb4LqORGMBERkcJKBbqviI9X83ZfFxsLxYvrKLqI+D1jzF3AXmvtmsucv7cxZrUxZvW+ffu8nE5ERKTwUoHuC06dgo0bVaD7utBQ6NoV3n0XDhxwO42ISG7cDLQxxmwH5uA0bR8HhBpj0q71WRlIympma+2b1tooa21U2bJl8yOviIhIoaAC3Rf88gucPasC3R/07QsnT8KMGW4nERG5bNbaQdbaytbaCKAj8LW1tguwDGjvmaw7sNCliCIiIoWSCnRfoB7c/Uf9+nDLLTBhgvOjiohIwTIAeMIYswXnnPRpLucREREpVAKzn0TyXHw8FCvm9BIuvu+xx6B9e1i8GCjidhoRkVyx1i4HlnsebwMau5lHRESkMNMRdF8QHw9160Kgfi/xC/fc4/S4P3as20lERERERKQAUYHuC9SDu38JDIRHH4UVK6izZ5vbaUREREREpIBQge62PXucmwp0/9KzJ5QoQY/Vi9xOIiIiIiIiBYQKdLf99JNzrwLdv4SGQmwsbX5Z7vzAIiIiIiIikksq0N2W1oN7vXru5pBL9+ijBJ1NgcmT3U4iIiIiIiIFgAp0t8XHQ4UKULas20nkUtWqxdfVomDSJDh1yu00IiIiIiLi51Sgu00dxPm16VH3OE3c5851O4qIiIiIiPg5FehuSkmBhAQV6H7s24gGUKeOc8k1a92OIyIiIiIifkwX3nbTr7/C6dMq0P2ZMfD44/DQQ/Dtt9C0qduJRERERMSLKoUGEzFwqdeXKZIVFehuSusgTgW6f+vaFQYNco6iq0AXERERKVC+G9jC7QhSiKiJu5vi4yEwEGrXdjuJ5Ebx4tC7N3z0EWzf7nYaERERERHxUyrQ3bRhA9SqBUWLup1EcqtvX6e5+4QJbicRERERERE/pQLdTQkJULeu2ynEGypXhg4dYOpUSE52O42IiIiIiPghFehuOX4cfvtNBXpB8vjjcPgwvP2220lERERERMQPqZO4fHTzqK9JOnQCgLp/bGGptTy89gSfXEKvkOrx0Yc1aQI33gjjxsHDD0NAgNuJRERERETEj6hAz0dJh06wfVRr58ns2fA2THqpG1x7rbvBxHuefBLuuw8WLYJ27dxOIyIiIiIifkRN3N2SkABFikCNGm4nEW9q1w6qVoVXX3U7iYiIiIiI+BkV6G5JSHB6cC9SxO0k4k2BgdC/P3z/vXMTERERERHJIRXoblEP7gXXAw9AWJiOoouIiIiIyCVRge6GY8fUg3tBVqIEPPIIfPQRbN7sdhoREREREfETKtDd8Msvzr0K9IKrXz/n9IUxY9xOIiIiIiIifkIFuhsSEpx7FegFV/ny0K0bzJwJ+/a5nUZERERERPyACnQ3JCRA0aJQvbrbSSQvPfEEnDwJEye6nURERERERPyACnQ3/Pyz04N7oC5DX6Bdey3cdRdMmAAnTridRkREREREfJwKdDeoB/fC4+mn4c8/4e233U4iIiIiIiI+TgV6fvv/9u4+TuZ6///448VaNisrpGw2cpmLw5ZcdE5FnS4Jucgqp7ZIUudQfFHfijollZzy1eVpuz5CovakfrrQNclVZBaFQ4iSXK0O1vb+/fEZ25KL2d2Z+czuPO+329xm5nMx+5y3WZ957ef9eb9zc2HtWhXo8eKcc+Css7zB4vLz/U4jIiIiIiIxTAV6tGkE9/hiBsOGedOt/fvffqcREREREZEYpgI92jSCe/zp3h3q1YOHHgLn/E4jIiIiIiIxKqQC3cwuMbOVZrbKzEYeZn1FM5sSXD/PzOoWWvcHM5trZgEz+9rMKoUvfikUCEDFihrBPZ4kJMDQoTB3Lnzyid9pREREREQkRh2zQDez8sDjwKVAU6CPmTU9ZLN+wDbnXAPgH8CDwX0TgFeAgc65ZkAHIC9s6UujQACaNIHy5f1OItF0/fVw4onwwAN+JxERERERkRgVyhn0NsAq59wa59w+YDLQ9ZBtugIHhqmeBlxgZgZcBCx1zi0BcM5tdc7F90hZgQA0PfTvG1LmJSV586LPmgULF/qdRkREREREYlAoE3GnAusLPd8AtD3SNs65/Wa2A6gONAKcmc0CagKTnXMPHfoDzGwAMAAgLS2tqO+h1Ki89xf47jsYMMDvKBImqSlJ1B05M6Rtk/eexpyKlfn06sHc3O32g17j85HnRyqiiIiIiIiUEqEU6CV9/T8BZwG/AB+Y2ULn3AeFN3LOPQM8A9C6desyO4pWg63Bv3NogLgyo8iFdeUldBozhk6Z9b1LHSDkAl9ERERERMq2ULq4bwTqFHp+SnDZYbcJXndeFdiKd7b9E+fcT865X4C3gTNKGrq0avTTd94DFejxa/BgqFQJHnzQ7yQiIiIiIhJjQinQ5wMNzayemSUCGUD2IdtkA9cGH/cEZjvnHDALaGFmxwUL9/OAnPBEL30a/vSdV5yddprfUcQvNWvCDTfAK694lzuIiIiIiIgEHbNAd87tB27BK7aXA1OdcwEzu9fMugQ3ywKqm9kq4DZgZHDfbcB4vCL/K2CRcy5u+/M2+uk7jeAu3pRrAOPG+ZtDRERERERiSkjXoDvn3sbrnl542d2FHu8Beh1h31fwplqLew1/+g7OvsjvGOK3tDT4y1/g2Wfhzjv9TiMiIiIiIjEilC7uEg47d5K6a4uuPxfPiBGwZw889pjfSUQkDplZJTP70syWmFnAzO4JLq9nZvPMbJWZTQle2iYiIiJRogI9WnKCl95rDnQBaNwYevSAxx+nyt7dfqcRkfizFzjfOdcSaAVcYmbtgAeBfzjnGgDbgH4+ZhQREYk7KtCjJRDw7nUGXQ644w7YsYPMBYeOuSgiElnOkxt8WiF4c8D5wLTg8heBbj7EExERiVsq0KMlJ4c9CYlQr57fSSRWpKdDly70W/Am7NjhdxoRiTNmVt7MvgJ+BN4DVgPbg4PDgjdVauoR9h1gZgvMbMGWLVuiE1hERCQOqECPlkCAVdXraAR3OdioUaTsyYX/+z+/k4hInHHO5TvnWgGnAG2AJkXY9xnnXGvnXOuaNWtGLKOIiEi8UYEeLYEA39RI8zuFxJozzuC9Bm1g/HjYudPvNCISh5xz24EPgfZAipkdmOHlFGCjb8FERETikAr0aNixAzZs4FsV6HIYj/3xKti2TWfRRSRqzKymmaUEHycBFwLL8Qr1nsHNrgXe9CehiIhIfFKBHg3BEdx1Bl0OZ9lJDaBzZ3jkEZ1FF5FoORn40MyWAvOB95xzbwEjgNvMbBVQHcjyMaOIiEjcSTj2JlJiwRHcv6lxqs9BJGaNGgVnnQUTJ3qju4uIRJBzbimQfpjla/CuRxcREREf6Ax6NAQCkJTEhqon+p1EYlXr1tCpk86ii4iIiIjEMRXo0RAIwOmn40zNLUcxahT8/LN3Fl1EREREROKOKsZoyMmBZs38TiGx7qyzfjuLvmuX32lERERERCTKVKBH2vbtsHGjCnQJzYGz6BrRXUREREQk7qhAj7TgCO4q0CUkZ53ljej+8MPe1GsiIiIiIhI3VKBHWnAEdxXoErL77/d6Xowb53cSERERERGJIhXokRYIwHHHwamaYk1C9Ic/QJ8+8OijsHmz32lERERERCTpl3JlAAAgAElEQVRKVKBHWnAEd8qpqaUI7rkH9u6FMWP8TiIiIiIiIlGiqjHSAgF1b5eia9gQrr8ennoK1q3zO42IiIiIiESBCvRI2rYNNm1SgS7Fc/fdXs+Le+7xO4mIiIiIiESBCvRI0gBxUhKnnAI33wwvvggrVvidRkREREREIkwFeiRpijUpqZEjvUEG777b7yQiIiIiIhJhKtAjKRCAypUhLc3vJFJa1awJt90Gr70Gixb5nUZERERERCJIBXokBQLQtKlGcJeSue02qF4dRowA5/xOIyIiIiIiEaLKMZI0gruEQ9WqcNdd8P77MGuW32lERERERCRCVKBHys8/w+bN3hl0kZK66SaoXx/+538gP9/vNCIiIiIiEgEq0CNFI7hLOCUmwtixsGwZvPCC32lERERERCQCVKBHigp0CbcePaB9e6+7++7dfqcREREREZEwU4EeKYEAJCdrBHcJHzMYNw42bYJHHvE7jYiIiIiIhJkK9EjJyfGuPzfzO4mUJWef7Z1Jf+ghb4wDEREREREpM1SgR4pGcJdIGTsW9u6FUaP8TiIiIiIiImGkAj0Stm6FH35QgS6R0aABDBoEzz7rDRonIiIiIiJlggr0SNAAcRJpd9/tzY8+ZAg453caEREREREJAxXokXCgQNcc6BIp1avDvffCBx/AjBl+pxERERERkTBI8DtAmRQIQJUqUKeO30mkLBs4EJ5+mu/730zHz2FvhYrFfqnUlCQ+H3l+GMOJiIiIiEhRqUCPhEBAI7hL5CUkwIQJ1D7/fFZWy4E77yz2S9UdOTOMwUREREREpDjUxT0SNIK7REvHjrzd6GwYMwbWr/c7jYiIiIiIlIDOoIfbli3eTQW6hCg1JalEZ7DP7HYzl03oB8OHw6uvhjGZiIiIiIhEkwr0cMvJ8e5VoEuIwnLtd6WV3qBxN90E555b8tcTEREREZGoUxf3cDswL7UKdImmESO8QQn/9jfYv9/vNCIiIiIiUgwq0MMtEIDjj4fUVL+TSDw57jgYPx6WLIGJE/1OIyIiIiIixRBSgW5ml5jZSjNbZWYjD7O+oplNCa6fZ2Z1D1mfZma5ZjYsPLFjWCAAzZtrBHeJvh494NJL4a67YMMGv9OIiIiIiEgRHbNAN7PywOPApUBToI+ZNT1ks37ANudcA+AfwIOHrB8PvFPyuDHOOa+Lu7q3ix/MvLPn+/fDkCF+pxERERERkSIK5Qx6G2CVc26Nc24fMBnoesg2XYEXg4+nAReYeaeQzawb8B8gEJ7IMeyHH+Dnn70z6CJ+OO007wz666/DTM1tLiIiIiJSmoRSoKcChSdY3hBcdthtnHP7gR1AdTNLBkYA9xztB5jZADNbYGYLtmzZEmr22BMI/g1CZ9DFT8OGwemnwy23wC+/+J1GRERERERCFOlB4kYD/3DO5R5tI+fcM8651s651jVr1oxwpAjSCO4SCxIT4cknYe1auO8+v9OIiIiIiEiIQinQNwJ1Cj0/JbjssNuYWQJQFdgKtAUeMrO1wBDgDjO7pYSZY1cgANWrQ61afieReHfeeXDttfDww5CT43caEREREREJQSgF+nygoZnVM7NEIAPIPmSbbODa4OOewGznOcc5V9c5Vxd4FBjjnCu7c0AFAt7Zc43gLrHg4Ye9Kf9uuAHy8/1OIyIxxMzqmNmHZpZjZgEzGxxcfoKZvWdm3wbvq/mdVUREJJ4cs0APXlN+CzALWA5Mdc4FzOxeM+sS3CwL75rzVcBtwO+mYivzNIK7xJqaNeHRR2HOHHj8cb/TiEhs2Q8Mdc41BdoBNwdnaBkJfOCcawh8QDwez0VERHyUEMpGzrm3gbcPWXZ3ocd7gF7HeI3RxchXemzcCDt3agR3iS19+8Krr8Ltt0Pnzt4o7yIS95xzm4BNwce7zGw53oCvXYEOwc1eBD7CG+xVREREoiDSg8TFD43gLrHIDJ5+GsqX97q6O+d3IhGJMWZWF0gH5gG1gsU7wGbgsIOqlJnZV0RERGKMCvRw0QjuEqvq1IFx42D2bHj2Wb/TiEgMCU6H+jowxDm3s/A655wDDvtXvTIz+4qIiEiMUYEeLoGAN3p7jRp+JxH5vRtugPPPh6FDYf16v9OISAwwswp4xfm/nHPTg4t/MLOTg+tPBn70K5+IiEg8CukadIE/jp3Nxu3/PeL6N/7f5+yuVIurR8484japKUmRiCZybGbwz39CixYwcCC89ZZmGxCJY2ZmeAO8LnfOjS+06sCsLGOD92/6EE9ERCRuqUAP0cbt/2Xt2E6HX/nrrzCxN/Trd+RtRPx22mnwwAMweDC88AJcd53fiUTEP38E/gJ8bWZfBZfdgVeYTzWzfsA64Eqf8omIiMQlFejh8N13sHu3rj+X2HfLLTBjhlekd+gA9er5nUhEfOCc+ww4UjeaC6KZRURERH6ja9DDQSO4S2lRrpx39twMrrkG8vP9TiQiIiIiIkEq0MNBI7hLaXLqqTBxInz2mTe6u4iIiIiIxAQV6OEQCEBqKqSk+J1EJDR9+0LPnnDXXfDVV8feXkREREREIk4FejgEAjp7LqWLGTz1lDctYN++VNy/z+9EIiIiIiJxTwV6SeXnQ04ONG/udxKRoqleHZ5/HgIBRnz0gt9pRERERETingr0kvrPf2DPHp1Bl9Lp4ovhb3/j+oXZkJ3tdxoRERERkbimAr2kDozgrjPoUlo99BDLatX35kVfv97vNCIiIiIicUsFekkdGMG9aVN/c4gUV8WK3NJlOOzbB336wP79ficSEREREYlLKtBLatkyb9qq5GS/k4gU29oTUuHpp+Hzz2HUKL/jiIiIiIjEJRXoJfX11/CHP/idQqTkrroKrr8eHngA3n/f7zQiIiIiInFHBXpJ7N0LK1aoQJeyY8IEaNIErr4aNm70O42IiIiISFxRgV4Sy5d706ypQJeyonJleO012L0bevXyrksXEREREZGoUIFeEkuXevcq0KUsadYMnnsO5s6FW2/1O42IiIiISNxQgV4SS5dCxYrQoIHfSUTC68orYdgweOIJeOklv9OIiIiIiMQFFeglsXSpd7YxIcHvJCLh98AD0LEj3HgjLF7sdxoRERERkTJPBXpJaAR3KcsSEmDyZKhRA7p3h61b/U4kIiIiIlKmqUAvrh9/hM2bVaBL2XbiifD66/D999CzpwaNExERERGJIBXoxfX11969CnQp69q08QaN++gjuPlmcM7vRCIiIiIiZZIuni4ujeAu8eTqq71pBe+/H04/HW67ze9EIiIiIiJljgr04lq6FE46CWrW9DuJSHTcey+sWOGN7t6oEXTu7HciEREREZEyRV3ci2vpUmjRwu8UItFTrhy8+CKkp0OfPr/1IhERERERkbBQgV4c+/dDIKDu7RJ/KleG7Gw4/njo1Ak2bPA7kYiIiIhImaECvTi+/Rb27lWBLvEpNRXefht27oRLLoFt2/xOJCIiIiJSJqhALw6N4C7xrmVLmDEDvvkGunaFPXv8TiQiIlIsAwcO5O9//7vfMUREAA0SVzxLl0L58t5o1iJlQGpKEnVHzizyfp0vvZWJ2Q/xTssLuLnrCH4tV75EGT4feX6x9xcRkdIrMzOTSZMmkZiYWLBs6NCh3HPPPSHt36FDB/r27Uv//v2L/LOfeuqpIu/jh/vuu49XX32V5cuX89xzz5GZmVmwbsOGDfTt25cvv/ySJk2a8NJLL9G8eXMA8vLyGDhwIFOnTqVatWqMGzeOK6+8smDfCRMmcP/997Nv3z4GDhzIAw88EO23JiKFqEAvjq++giZNoGJFv5OIhEXxC+NO8FhtLh0yhDW7Z8HEiWBWrFcqzh8IRESk7Bg+fDj33Xef3zFi1mmnncaECRMYOXLk79YNGDCA5s2b88477zBhwgR69+5NIBAA4B//+Adff/0169ev56uvvqJz5860b9+eOnXqMG/ePEaPHs2nn35KSkoKf/rTn0hPTz+ogBeR6FIX9+JYvNgbyVpEYPBgGD4cnngC7rzT7zQiIlLGvPDCC7Rv357MzEyqVKnCH//4R3766ScAxowZQ3JyMp9++im33HILycnJnHnmmQftX7duXSZOnMiZZ55J5cqV6dKlCwBvvfUWycnJVKhQgTsPc/yaPn06zZs3p1q1alx22WVs2rSpYN1bb71F48aNqVKlCg0bNmTWrFkRbAHPVVddxQUXXEDFQ04Q7dy5k3fffZeRI0eSlJTErbfeyrp161ganG3ltdde429/+xspKSl06NCB9u3bM2PGDACmTZtG9+7dadasGampqfTv35/JkydH/L2IyJGpQC+qH36A779XgS5S2NixcOONMGYM3H+/32lERKSMmT9/PldffTU//vgj+/fv57nnngPgjjvuIDc3l3POOYeJEyeSm5vLwoULf7f/008/zUsvvcSOHTu44447AOjcuTO5ublcffXVv9v+yy+/pF+/fmRlZbFlyxbS09MZMGBAwfobbriBe++9l127dvHuu++SmpoaoXd+bKtWraJSpUokJydzzjnnsGbNGurXr8/KlSsBWLlyJU2aNKFv375MnjyZpk2b/m7dY489xtChQw9aJyL+UIFeVIsXe/dnnOFvDpFYYuadQf/LX7yz6OPH+51IRERKmXHjxpGSklJwmz59esG6Ro0aceGFF5KUlMQFF1zAN998U6TXHjBgAM2aNSMhIYF27dodc/usrCwyMzNp27YtCQkJDBs2jJkzZ7J3714AypUrx+rVq9m5cyf16tUruN7bD7t37yY5OZldu3YRCATYtm0bVapUITc396D1X3/9NRs3bjzsutWrV/Ptt98etE5E/KECvagOFOitWvmbQyTWlCsHzz0HPXvC0KHw5JN+JxIRkVJk2LBhbN++veDWvXv3gnUnnHBCwePExET2FHH2kIYNGxZp+/Xr1/P0008X/LGgXr16JCYmFnRzf+2115g7dy5paWm0a9eOZcuWFen1w6ly5crk5uZSp04dfv75Z9q3b8+uXbtITk4+aP2SJUsYOnToYddNmDCB7Ozsg9aJiD9UoBfVokVw2mmQkuJ3EpHYk5AA//oXdOoEgwZBVpbfiUREJA6UK3f0r7QJCUUbF7lOnTrcddddB/3BYM+ePdStWxeAs88+m3//+9/88MMPNGnS5LDXsEdLgwYN+O9//8uGDRsA2LdvH6tXr6Zx48aA1/tgxYoVBdvn5OSEtE5E/KFR3ItKA8SJHF1iIkybBt26Qf/+sG8f3HTTMXcr7lRvh76GpmoTEYk/J510UljPYmdmZpKRkcHFF19Meno6W7Zs4cMPP6R37978+uuvTJo0ia5du1KpUiUAjj/++LD97CPJy8sjPz8f5xx5eXns2bOHxMREjj/+eC6++GLGjh3LuHHjeOyxx0hLS6NFixYAXHnllUyYMIHOnTuzZMkS5s6dW3ANf69evbj00ku59dZbSUlJISsrS9OsifhMBXpR7NgBq1fD9df7nUQktlWqBG+8Ab16eWfS9+6FIUOOuks4CmtN1SYiUno99NBDPProowXPMzIyePbZZ0Pad+jQoVxzzTXUrl2bRo0a8dFHHx1zn4suuog5c+awd+9ezIxHH32Unj17FowaP378eK677jr+85//kJKSQkZGBr179wbglVde4ZZbbsE5R9u2bcmKQo+xG264gRdffBGAOXPmMGDAAD788EM6dOjA008/Td++falWrRpNmjRhypQpWHDa01tvvZUVK1ZQp06dgiK8Tp06ALRt25ZRo0bRsWNH8vLyuPHGGwveo4j4w5xzx97I7BLgMaA88Kxzbuwh6ysCLwFnAluB3s65tWZ2ITAWSAT2Af/jnJt9tJ/VunVrt2DBguK8l4iqO3Imay9Nhg4d4O234dJL/Y4kEvv27YM+fWD6dG+k9xEjIvrj6o6cydqxnSL6M0RihZktdM619jtHJI7b+l0WERE/ReI4FOpx+5jXoJtZeeBx4FKgKdDHzJoeslk/YJtzrgHwD+DB4PKfgMudcy2Aa4GXQ38LMWjRIu9eXdxFQpOYCFOmQEYGjBwJo0dDCH8UFBERERGJR6F0cW8DrHLOrQEws8lAVyCn0DZdgdHBx9OAiWZmzrnFhbYJAElmVtE5t7fEyf2weDGcfDKcdJLfSURKj4QEeOUVr9v7PffAli0wYQKUL+93MhERERGRmBLKKO6pwPpCzzcElx12G+fcfmAHUP2QbXoAiw5XnJvZADNbYGYLtmzZEmr26NMAcSLFU768NwXb8OHefOkZGVDEKXJEJLzM7Dkz+9HMlhVadoKZvWdm3wbvq/mZUUREJN5EZZo1M2uG1+39xsOtd84945xr7ZxrXbNmzWhEKrKKeXth+XI44wy/o4iUTmbw4IPwyCPeKO+XXuoNvCgifnkBuOSQZSOBD5xzDYEPgs9FREQkSkIp0DcCdQo9PyW47LDbmFkCUBVvsDjM7BRgBnCNc251SQP7pcmWtZCfrzPoIiV1221el/fPPoPzzoPgvK0iEl3OuU+Anw9Z3BV4Mfj4RaBbVEOJiIjEuVAK9PlAQzOrZ2aJQAaQfcg22XiDwAH0BGY755yZpQAzgZHOuc/DFdoPLTav8h7oDLpIyV19Nbz1FqxZA23aQAzO3CASp2o55zYFH28Gah1uo1JzaZqUCh999BFmxohCM320adOmYJqweDJhwgRq1apFtWrVuP3224u8/6effoqZHTQ93U8//cQVV1xBtWrVqF27NqNHjy5YN3DgQJKTkwtuFStWLJg/XUT8ccwCPXhN+S3ALGA5MNU5FzCze82sS3CzLKC6ma0CbuO3LnG3AA2Au83sq+DtxLC/iyhotekbOPFEOPVUv6OIlA0XXwxz5ngjvZ97rtftXURihvPmYT3stAul4dI0KV1q1qzJ+++/D8B3333Hjji8BGrevHmMHj2a2bNns2zZMiZPnszUqVND3n///v2MHDmSJk2aHLT8zjvvJC8vj++//54vvviCZ555hpkzZwLw1FNPkZubW3Dr2bMnPXr0COv7EpGiCekadOfc2865Rs65+s65+4PL7nbOZQcf73HO9XLONXDOtTkw4rtz7j7nXGXnXKtCtx8j93Yip9X3K70zfXH411yRiGneHL780rt0pFcvuP9+TcMm4q8fzOxkgOB9qTxmS+mTmJhI06ZNWbx4MdOmTaN79+4HrZ8+fTrNmzenWrVqXHbZZWzatKlgXUZGBrVq1aJy5cq0b9+epUuXFqyrW7cuo0aNokGDBpxwwgk89dRTUXtPRXXgfTdr1ozU1FT69+/P5MmTQ97///7v/+jUqRO1ah3c8WXdunV07tyZpKQk0tLSaN++PTk5Ob/bf8eOHcyYMYNrrrmmxO9FRIovKoPElXo7dtDg5w3Qtq3fSUTKnhNPhA8+gL594c47ve7vu3f7nUokXhW+ZO1a4E0fs0ic6dmzJ9OmTeONN96ga9euBcu//PJL+vXrR1ZWFlu2bCE9PZ0BAwYUrE9PT+frr79m586dnHvuuVx77bUHve68efNYunQpWVlZDB8+nP3790ftPRXFypUradKkCY899hhDhw6ladOmrFy5MqR9N2/ezPPPP89tt932u3WDBg3i3XffZffu3axdu5aFCxdy4YUX/m67yZMnc8YZZ3DaaaeV+L2ISPGpQA/Fgetj27TxN4dIWVWpErz0EowZA5MnQ7t28M03fqcSKdPM7FVgLtDYzDaYWT9gLHChmX0L/Dn4XCQqLr74Yl5//XUqVqxIjRo1CpZnZWWRmZlJ27ZtSUhIYNiwYcycOZO9e72Ze0eMGMGJJ55I+fLlyczMPOgMOsD111/PcccdR6dOndi1axebN2+O6vsK1e7du0lOTmb16tV8++23VKlShdzc3JD2HTZsGHfccQeVKlX63br09HR27NhB1apVqVevHv3796dVq1a/2+7555/X2XORGKACPRTz5nn3Z53lbw6RsswMbr8dZs2CTZu837cZM/xOJVJmOef6OOdOds5VcM6d4pzLcs5tdc5d4Jxr6Jz7s3Pu0FHeRSKmUqVKdOnSheuuu+6g5evXr+fpp58mJSWFlJQU6tWrR2JiIps2bSI/P5/bb7+d+vXrk5KSQtu2bfn111/Jz88v2P+EE04AvG70AHv27InemyqCypUrk5uby4QJE8jOzmbXrl0kJycfc7/PP/+cNWvWkJGRcdj1GRkZtGzZkt27d7Nu3TomTZr0u2vbly9fzpIlS7jyyivD8l5EpPhUoIfiyy9ZfUIqVKvmdxKRsu/CC2HRImjSBLp3h+HDIUa7I4qISHg99NBDXHXVVQctq1OnDnfddRfbt28vuO3Zs4e6desyadIkZsyYwYcffsj27dv57LPPAHClcDyTRo0asWLFioLnOTk5NG7c+Jj7zZ8/n7lz52JmmBkff/wxN9xwA0OGDAFg0aJFZGZmUrFiRdLS0ujUqRMffPDBQa/x/PPPc/nll5OSkhLeNyUiRaYC/Vicg3nzWHJyI7+TiMSPtDT45BO46SZ4+GFvvvS1a/1OJSIiPsjMzOSpp55i0aJFOOf48ccfmTJlCgC7du0iKSmJatWqsWvXLsaMGeNz2uLr1asX06dPJxAIsHHjRrKysujdu/dB20yfPp0GDRqwcePGgmVDhgzBOVdwO++88/jnP//Jo48+Cnhd3F9++WXy8vLYvHkzs2bNomnTpgX75+fn8/LLL6t7u0iMUIF+LBs3wubNKtBFoq1iRXjiCZg0CZYtg5YtvcciIhJX2rdvz/jx47nuuuuoWrUqrVu3ZuHChQBcc801pKWlUbt2bVq2bEm7du18Tlt8bdu2ZdSoUXTs2JHmzZtz5ZVX/q5A37lzJ6tXryYvLy/k133hhRf46quvqFmzJq1ateK8885j0KBBBevfeecd8vPzueSSS8L2XkSk+CzWugC1bt3aLTgwKFsseP116NmTrn95hDdf+v3ImCISBWvXeqO8f/65d//443D88b/brO7Imawd2yn6+UR8YGYLnXOt/c4RieO2fpdFRMRPkTgOhXrc1hn0Y/nyS6hQgeUnasoJEd/UrQsffQSjR3tn0Vu18p6LiIiIiJQhKtCPZe5cSE9nX0IFv5OIxLeEBBg1Cj79FMqVg44d4eabYdcuv5OJiIiIiISFCvSj2bvXO4N+zjl+JxGRA84+G5YuhVtvhSefhBYt4L33/E4lIiIiIlJiKtCPZv58r0j/05/8TiIihR13HIwfD599BpUqwUUXQb9+VPtlh9/JRERERESKTQX60QTn0lSBLhKjzj4bFi+GESPgpZeY/c+B8Mwz8OuvficTERERESkyFehH8+mncPrpUKOG30lE5EiSkmDsWPjqK1bWPBVuvBHatYPgFDwiIhL76tatS1JSEsnJydSvX5+JEycWrDMz2rZtW/B8+PDhmBkfBQcLzc7OpkWLFlSpUoVTTz2V+++/P9rxw2bChAnUqlWLatWqcfvtt4e0T15eHtdccw0nn3wyVatWpWPHjgQCgYL1P/30E1dccQXVqlWjdu3ajB49umDdvn37uPnmm6lVqxbVq1dn+PDhxNoMTyLxRgX6kfz6qzelk86ei5QOzZqR0ecBeOUV+O47OOssGDgQfvjB72QiIhKCf//73+Tm5vLqq68yYsQIPv7444J127dvZ/369QB88MEH1KxZE4A1a9ZwzTXX8Pjjj7Nz504+/fRTGjRo4Ev+kpo3bx6jR49m9uzZLFu2jMmTJzN16tRj7pefn0+DBg2YP38+27Zto0uXLnTr1q1g/Z133kleXh7ff/89X3zxBc888wwzZ84E4Mknn2TevHmsWLGCxYsXM2XKFP71r39F7D2KyLGpQD+SZctgxw4NECdSmpjB1VfDypXwt79BVhY0aAD33Qe//OJ3OhERCUGbNm1o1qwZ8+fPL1jWvXt3pk2bxqJFi2jatCmJiYkALFq0iHr16nHuuediZqSlpdG7d2+/opfItGnT6N69O82aNSM1NZX+/fszefLkY+5XqVIl7r77bk455RTKlSvHddddx6pVq9iyZQsA69ato3PnziQlJZGWlkb79u3JyckB4OOPP+aqq66iWrVqpKWl0adPH6ZPnx7R9ykiR6cC/Uh0/blI6VW1Kjz6KAQCcOGFcNdd0LAhPP885Of7nU5ERI5i4cKFBAIBTj/99IJl3bp1480332TatGn06NGjYHmLFi0IBAKMGjWKJUuWlOru2StXrqRJkyY89thjDB06lKZNm7Jy5coiv87cuXOpXbs21atXB2DQoEG8++677N69m7Vr17Jw4UIuvPBCgMO216pVq0r2RkSkRFSgH8knn0BqKtSt63cSESmuRo1g+nRvPIk6deD666FlS5g2TQPJiYjEmG7dunHcccdx7rnnMnbsWDp16lSwrnr16lSsWJFp06Zx8cUXFyxv3Lgxb7/9NnPmzKFt27bUq1ePN99804/4JbZ7926Sk5NZvXo13377LVWqVCE3N7dIr7Fjxw4GDx7M+PHjKVfO+5qfnp7Ojh07qFq1KvXq1aN///60atUKgA4dOjBp0iS2bt3K2rVree211/hFPc5EfKUC/XB+/RVmz4bzz/e6zIpI6fanP8HcuTBlCuzfD716QXq6V7yrUBcRiQlvvPEGu3btYsiQIcyZM+d366+77jq6dOlCUlLSQcv//Oc/895777Ft2zaGDx9Onz592Lp1a7Rih03lypXJzc1lwoQJZGdns2vXLpKTk0Pef+/evVxxxRX07t37oG7+GRkZtGzZkt27d7Nu3TomTZpUcG37TTfdRLt27WjRogWdOnWiW7duHH/88WF/byISOhXoh/P117BlC/z5z34nEZFwMYMrr/S6vb/yCuzZAz16wBlnwOuvq+u7iEgMKF++PKNHj2bu3LkFo7QfkJGRwbhx4464b1JSEoMGDaJSpUqsWbMmwknDr1GjRqxYsaLgeU5ODo0bN0R/DvwAAA89SURBVA5p3/z8fPr06UODBg34+9//ftC6RYsWkZmZScWKFUlLS6NTp0588MEHACQmJjJx4kS+//57AoEAZkbLli3D96ZEpMgS/A4Qk957z7u/4AJ/c4hI+JUv7w0k17s3TJ4M994LPXt616jfdhsdt5zKf3aX7Kx6akoSn488P0yBRUTiS4UKFfjrX//KvffeS4cOHY667ZIlS1iwYAE9evTg+OOPZ8qUKeTl5dGoUaPohA2jXr16cemll3LrrbeSkpJCVlYWDzzwwEHbTJ8+neHDh/Pxxx+TmppasPzGG2/EzHjyySd/97rp6em8/PLLjBkzhq1btzJr1iz69esHeN3qN23axGmnncbs2bN59tlneeeddyL7RkXkqFSgH87773vznxf6j09EypiEBOjbFzIyvK7uDz8MN93Ea8dVpcaI22DQIKhRo1gvXXfkzDCHFRGJL/369Ss4k340ycnJTJ06leHDhxcU5tOnT6dq1apRSho+bdu2ZdSoUXTs2JG8vDxuvPHG341Iv3PnTlavXk1eXl7BsnXr1pGVlUVSUtJB7/udd97hnHPO4YUXXmDQoEHUrFmTSpUq0atXLwYNGgTAL7/8QqdOnVi/fj21a9fmiSeeoH379tF5wyJyWCrQD7V3rzdA3A03+J1ERKIhIcHr+t6rF3z8MUv6D+eCUaNg7Fi46iqvUD/jDL9TioiUaWvXrj3oeUpKCrt27QIOP9L4hg0bCh7PmjUrotmiafDgwQwePPiI6zMzM8nMzDxo2amnnnrU0esbNmzIewd6hx6iZs2axRopXkQiR9egH2rOHPjvf3X9uUi8MYMOHejXcxQsW+adXX/1VTjzTGjXDl580fu/QUREREQkQlSgH2rmTEhMhGNc8yQiZVizZvDMM7BxI0yYADt2QGYmnHIKDBkCixdDKZ5rV0RERERikwr0wpyDN9+Ejh2hShW/04iI31JS4K9/hZwc+PBDr2fNk096Xd5btoRHHoHNm/1OKSIiIiJlhAr0wlauhFWroEsXv5OISCwJdn9nyhTYtAmeeAKOOw6GDfMGk7zsMnjpJdi+3e+kIiIiIlKKaZC4wrKzvfvLL/c3h4jErhNOgJtu8m4rVniF+SuvwLXXQoUKcNFF3JDUlD8MyWVnpeRi/YhwTNP2x7Gz2bi9ZNfMa7o4ERERkehSgV5Ydjakp0OdOn4nEZHSoEkTGDMG7rsPvvwSpk2D117jf7+byf9WqOB1ib/8cujUCdLSQn7ZcEzTtnH7f1k7tlOJXkPTxYmIiIhEl7q4H7BpkzeCu7q3i0hRlSvnjfQ+bhysXQvz5nmDya1c6U3Tduqp0KIFjBwJn34K+/f7nVhEREREYpAK9AMmT/YGicvI8DuJiJRmZtCmDTz0kDemRU6OV7jXrOkNKnfuud7jXr28AedWrtSI8CIiIiICqIv7b/71L29k5iZN/E4iImWFGZx+uncbOtSbru2997zpHN9/3+sSD1C7Npx/vjeDxPm65ltEREQkXqlAB+8M1sKF3tktEZFIqVoVevb0bs7B6tUwe7Y3hdu773qDzQFzk6vD6o7Qvr13O+MMqFjR5/AiIiIiEmkq0AGysqB8eXVvF5HoMYMGDbzbgAFewR6cb/3LJ6bRdcGC386wJybCmWd6xXrr1t5glg0bev9viYiIiEiZERcF+tGmG6qYt5cvnniKufXbMmjCYmDxYbdLTUmKYEIRCYfUlKQSjzzu2++6GTRrBs2aMXhDPbqO7eQNXjl37m+3xx+HvXu97StXhpYtvWI9Pd07y960aVjPtIerPTVVm4iIiEho4qJAP+p0Qy+8AHt2cdlT97G2Y8eo5hKR8CpzheDJJ0P37t4NIC/PO8u+eLF3W7QIXnzRK9zBO6PeoAE0bcrQzYkwaYdXtDduDElF/8NDONpTU7WJiIiIhC4uCvQj+vVX77rzpk2hQwe/04iIHF2FCt5Z85YtITPTW/brr9617IsWwbJlXgGfk8NNK7+BuVO8bcygXj1vEMz69b0ivn5971avnq5vFxEREYkR8V2gT53qfaGdNMn7AisiEgOK17U8GWgHDdtBQzg1OYGPrziloGAnEIBvvoFPPoHc3N92M4M6dX4r3OvW9eZtP/VUaNQITjwxjO9MRERERI4mfgv0vDwYPdq75rN3b7/TiIgUCGtX/WbNDn7uHGzZ4p11X73am6v9wOM334Qff/xt28GD4dFHw5dFRERERI4qfgv0ceO86dWys6FcOb/TiIhEh5l3VvzEE71R4Q/1yy/w3Xewbp03P7uIiIiIRE18Fujz58M990CvXnD55X6nERGJHccd512r3qSJ30lERERE4k5IBbqZXQI8BpQHnnXOjT1kfUXgJeBMYCvQ2zm3NrjudqAfkA/8zTk3K2zpi2P5cujWDU466beRj0VEJCLCMVVbODKUuRH+I+xYx30RERGJjGMW6GZWHngcuBDYAMw3s2znXE6hzfoB25xzDcwsA3gQ6G1mTYEMoBlQG3jfzBo55/LD/UaOyjnYuhVefRXuussbsXjWLKhZM6oxRETiTSwUxn7/gaC0CfG4LyIiIhEQysXXbYBVzrk1zrl9wGSg6yHbdAVeDD6eBlxgZhZcPtk5t9c59x9gVfD1oue551gxvgfUqAF//Su0agVz5kDz5lGNISIiUkqEctwXERGRCAili3sqsL7Q8w1A2yNt45zbb2Y7gOrB5V8csm/qoT/AzAYAA4JPc81sZUjpQ1cD+AmAjz7yphOSkvitPSUc1J7ho7YMrzLTnvag3wki0panhvn1DjjmcT8ax217sGx89mJEmfldjgFqy/BSe4aX2jN8InEcCum4HRODxDnnngGeidTrm9kC51zrSL1+vFF7hpfaM3zUluGl9gyfstaWOm6XLmrP8FFbhpfaM7zUnuHjZ1uG0sV9I1Cn0PNTgssOu42ZJQBV8QaLC2VfERERiR06douIiPgklAJ9PtDQzOqZWSLeoG/Zh2yTDVwbfNwTmO2cc8HlGWZW0czqAQ2BL8MTXURERCIglOO+iIiIRMAxu7gHrym/BZiFN93Kc865gJndCyxwzmUDWcDLZrYK+BnvYE5wu6lADrAfuDnqI7h7ItYNL06pPcNL7Rk+asvwUnuGT6lpyyMd96Mco9S0Vymh9gwftWV4qT3DS+0ZPr61pXknukVERERERETET6F0cRcRERERERGRCFOBLiIiIiIiIhIDynyBbmaXmNlKM1tlZiP9zlPamdlaM/vazL4yswV+5ylNzOw5M/vRzJYVWnaCmb1nZt8G76v5mbE0OUJ7jjazjcHP51dmdpmfGUsLM6tjZh+aWY6ZBcxscHC5Pp/FcJT21OfzEMc6RgcHmZ0SXD/PzOpGP2XpEEJb3hb8TC41sw/MLKT5eONVqN8fzayHmTkz09RWRxFKe5rZlYX+35wU7YylRQi/62nBY9Di4O973B9rjuRw3yUPWW9mNiHY1kvN7Iyo5CrL16CbWXngG+BCYAPeyLR9nHM5vgYrxcxsLdDaOfeT31lKGzM7F8gFXnLONQ8uewj42Tk3NvifbDXn3Ag/c5YWR2jP0UCuc26cn9lKGzM7GTjZObfIzKoAC4FuQCb6fBbZUdrzSvT5LBDKMdrMBgF/cM4NNLMM4ArnXG9fAsewENuyIzDPOfeLmd0EdFBbHl6o3x+Dv98zgUTgFuecTlwcRoifz4bAVOB859w2MzvROfejL4FjWIht+Qyw2Dn3pJk1Bd52ztX1I2+sO9x3yUPWXwb8FbgMaAs85pxrG+lcZf0MehtglXNujXNuHzAZ6OpzJolTzrlP8GY5KKwr8GLw8Yt4X+IlBEdoTykG59wm59yi4ONdwHIgFX0+i+Uo7SkHC+UYXfgzOA24wMwsihlLi2O2pXPuQ+fcL8GnX+DNby+HF+r3x78DDwJ7ohmuFAqlPW8AHnfObQNQcX5EobSlA44PPq4KfB/FfKVKCN8lu+IV78459wWQEvwjfESV9QI9FVhf6PkG9CWppBzwrpktNLMBfocpA2o55zYFH28GavkZpoy4JdgN6Tl1yS66YBfidGAe+nyW2CHtCfp8FhbKMbpgG+fcfmAHUD0q6UqXon7f6Qe8E9FEpdsx2zPY1bWOc25mNIOVUqF8PhsBjczsczP7wswuiVq60iWUthwN9DWzDcDbeGeApXh8qSXLeoEu4fcn59wZwKXAzcGuIRIGzrvepOxecxIdTwL1gVbAJuARf+OULmaWDLwODHHO7Sy8Tp/PojtMe+rzKb4zs75Aa+Bhv7OUVmZWDhgPDPU7SxmSADQEOgB9gH+aWYqviUqvPsALzrlT8Lpmvxz8zEopUdb/sTYCdQo9PyW4TIrJObcxeP8jMAOvq40U3w8HusoE79WlqwSccz845/Kdc78C/0Sfz5CZWQW8YvJfzrnpwcX6fBbT4dpTn8/fCeUYXbCNmSXgddfcGpV0pUtI33fM7M/A/wJdnHN7o5StNDpWe1YBmgMfBcfmaQdka6C4Iwrl87kByHbO5Tnn/oN3nXXDKOUrTUJpy3541/PjnJsLVAJqRCVd2eNLLVnWC/T5QEMzq2dmiUAGkO1zplLLzCoHB0TBzCoDFwGHHfVQQpYNXBt8fC3wpo9ZSr1Drgu6An0+QxK8pjcLWO6cG19olT6fxXCk9tTn83dCOUYX/gz2BGa7sjy6bfEdsy3NLB14Gq841x/bju6o7emc2+Gcq+GcqxscfOsLvHbVIHGHF8rv+ht4Z88xsxp4Xd7XRDNkKRFKW34HXABgZqfjFehbopqy7MgGrgmO5t4O2FHo0r+ISYj0D/CTc26/md0CzALKA8855wI+xyrNagEzguPzJACTnHP/z99IpYeZvYp38KkRvC5oFDAWmGpm/YB1eKM8SwiO0J4dzKwVXlfstcCNvgUsXf4I/AX42sy+Ci67A30+i+tI7dlHn8/fHOkYbWb3Agucc9l4f+h42cxW4Q3kk+Ff4tgVYls+DCQDrwWP498557r4FjqGhdieEqIQ23MWcJGZ5QD5wP8459Rb5hAhtuVQvEsEbsU73mTqD5uHd4TvkhUAnHNP4V3DfxmwCvgFuC4qufTvJSIiIiIiIuK/st7FXURERERERKRUUIEuIiIiIiIiEgNUoIuIiIiIiIjEABXoIiIiIiIiIjFABbqIiIiIiIhIDFCBLiIiIiIiIhIDVKCLiIiIiIiIxID/Dx0ut8Z9oei3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if (Nexp > 1) :\n",
    "    \n",
    "    fig2, ax2 = plt.subplots(nrows=2, ncols=2, figsize=(14,10))\n",
    "    ax2 = ax2.flatten()\n",
    "\n",
    "    # Distribution of alpha0 values:\n",
    "    ax2[0].hist(array_alpha0, bins=80, range=(2.8, 4.4), histtype='step')\n",
    "    ax2[0].set_title('Histogram of alpha0')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_alpha0), array_alpha0.mean(), array_alpha0.std(ddof=1))\n",
    "    ax2[0].text(0.62, 0.95, string, family='monospace', transform=ax2[0].transAxes, fontsize=12, verticalalignment='top')\n",
    "\n",
    "    # Distribution of alpha1 values:\n",
    "    ax2[1].hist(array_alpha1, bins=60, range=(-0.05, 0.25), histtype='step')\n",
    "    ax2[1].set_title('Histogram of alpha1')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_alpha1), array_alpha1.mean(), array_alpha1.std(ddof=1))\n",
    "    ax2[1].text(0.05, 0.95, string, family='monospace', transform=ax2[1].transAxes, fontsize=12, verticalalignment='top')\n",
    "\n",
    "    # Distribution (normed) of Chi2 values, along with Chi2-distribution (for Ndof=7) drawn on top:\n",
    "    ax2[2].hist(array_Chi2, bins=25, range=(0.0, 25.0), histtype='step', density=True)\n",
    "    ax2[2].set_title('Histogram of Chi2')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_Chi2), array_Chi2.mean(), array_Chi2.std(ddof=1))\n",
    "    ax2[2].text(0.62, 0.95, string, family='monospace', transform=ax2[2].transAxes, fontsize=12, verticalalignment='top')\n",
    "    x_axis = np.linspace(0, 20, 1000)\n",
    "    y_axis = stats.chi2.pdf(x_axis, df=7)\n",
    "    ax2[2].plot(x_axis, y_axis, 'r-')\n",
    "    ax2[2].text(0.45, 0.65, \"Note: This is not a fit!\", transform=ax2[2].transAxes, fontsize=16, verticalalignment='top')\n",
    "\n",
    "    # Distribution of Chi2 fit probabilities:\n",
    "    ax2[3].hist(array_Prob, bins=20, range=(0.0, 1.0), histtype='step')\n",
    "    ax2[3].set_title('Histogram of Prob')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMS {:>10.3f}\".format(len(array_Prob), array_Prob.mean(), array_Prob.std(ddof=1))\n",
    "    ax2[3].text(0.3, 0.2, string, family='monospace', transform=ax2[3].transAxes, fontsize=12, verticalalignment='top')\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    if (save_plots) :\n",
    "        fig2.savefig('FitResultDistributions.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remark on $\\chi^2$ rule of thumb:\n",
    "A general rule is that the \"reduced chisquare\" (defined as $chi^2$/Ndof) should be around 1. That is true, but the rule definitely has its limitations, since \"around\" depends on the size of Ndof! This is examplified below, considering two cases with the same reduced chisquare of 1.5, but with **very** different conclusions.\n",
    "\n",
    "The reason for this can be found in the mean and width of the $\\chi^2$ distribution (with Ndof = k), which are k and sqrt(2k). From the central limit theorem, the $\\chi^2$ should approach a Gaussian for \"large\" values of k (k > 50 is probably enough), but the **relative** width (i.e. width/mean) becomes smaller and smaller as k grows, and a reduced chisquare of e.g. 1.5 becomes more and more unlikely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  Considering the relation between Chi2/Ndof ratio and Prob(Chi2,Ndof) vs. Ndof:\n",
      "  Low number of degrees of freedom:   Prob(  3.0,  2) = 0.2231302\n",
      "  High number of degrees of freedom:  Prob(300.0,200) = 0.0000059\n"
     ]
    }
   ],
   "source": [
    "# Numbers to illustrate the point, that using the ratio Chi2/Ndof (\"reduced chisquare\") is only a rule of thumb:\n",
    "print(\"\\n  Considering the relation between Chi2/Ndof ratio and Prob(Chi2,Ndof) vs. Ndof:\")\n",
    "print(\"  Low number of degrees of freedom:   Prob(  3.0,  2) = %9.7f\"%(stats.chi2.sf(  3.0,   2)))\n",
    "print(\"  High number of degrees of freedom:  Prob(300.0,200) = %9.7f\"%(stats.chi2.sf(300.0, 200)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------------------------------------------- \n",
    "\n",
    "\n",
    "\n",
    "# Questions:\n",
    "\n",
    " 1) Run the code such that it does exactly one fit (initial setting), and take a look at the line fitted.\n",
    "    Does this look reasonable, and are you satisfied with the fit probability? What is the chance that the input for the data could\n",
    "    actually be from a flat distribution, i.e. without the slope? [There are actually two simple/quick ways of testing this].\n",
    "\n",
    " 2) Does the fit reproduce the input numbers well (include the uncertainties in your\n",
    "    answer)?\n",
    "\n",
    " 3) Now increase the number of experiments to e.g. 1000, and rerun the macro. Figure\n",
    "    out what you see in the window split in 2-by-2, and go through each of these to\n",
    "    see, if you understand every feature of the distributions shown, and if you are\n",
    "    happy with them! Is this what you expected? Are they also as wide as you would expect?\n",
    "    Possibly try to make some fits of these output distributions to check, if you're satisfied!\n",
    "\n",
    "    This somehow makes this the \"long\" question without any leading questions.\n",
    "    But you should try to be statistically minded enough and have an idea of what to\n",
    "    look for, at least to some degree. We will discuss this further in class.\n",
    "\n",
    " 4) Investigate if the distributions of probabilities is flat, or if it has some\n",
    "    slope to it. Do you understand why the distributions of probabilities is flat?\n",
    "    (This is typically conceptually hard to begin with, but don't worry!)\n",
    "\n",
    " 5) Find the line where the random points for the line are generated, and add a\n",
    "    quadratic term in x with the coefficient 0.01. Run the program again with Nexp=1,\n",
    "    and start by looking at the single fit in the graph. Can you see this change?\n",
    "    Now run 1000 experiments again, and see what has changed in the distributions\n",
    "    in the 2-by-2 window when including such a term. Are these changes as you would\n",
    "    have expected them?\n",
    "\n",
    "    NOTE: The reason why the quadratic term is more clear in the latter case is NOT\n",
    "    that it is a better method. It is simply because there is 1000 times more statistics!\n",
    "\n",
    "\n",
    "### Advanced questions:\n",
    "\n",
    " 1) Change the coefficient from question 5) to -0.02. Of course the linear fit does\n",
    "    not do very well. What changes are needed for the fit to be good again? Make\n",
    "    these changes, and see that the condition in question 4) is again met.\n"
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}