Notes on information theory

Katie-fisher-information.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fisher information matrix is the expectation value of the variance of the log-likelihood function, which is often just the expectation of the Hessian of the log-likelihood function (i.e. the second derivatives with respect to the model parameters). That's just a boring definition! For a data set $D$ and model parameters $\\theta$, \n",
    "$$\n",
    "\\mathcal{F}_{jk} = -\\mathbb{E}\\left[\\frac{\\partial^2}{\\partial\\theta_j \\, \\partial\\theta_k}\\ln\\mathcal{L}(D \\,|\\, \\theta)\\right]\n",
    "$$\n",
    "where $\\mathbb{E}\\left[x\\right]$ is the expectation value of $x$ evaluated at the true parameter values.\n",
    "Conceptually, it is a matrix that tells us how much information the data $D$ contain about a given parameter, or how much covariance there is in our knowledge of a combination of parameters. \n",
    "\n",
    "I think of it as a result of a Taylor expansion of the likelihood function around the maximum likelihood parameter values. It is related to the curvature of the likelihood function around its maximum (because second derivatives are involves). It is not a fundamental quantity that is universally interpretable or meaningful. For example, if a model is a bad representation of data, or if the likelihood surface is extremely structured or multi-modal (i.e. cases where a Taylor expansion to low order is a bad representation of a function), then the Fisher information does not really summarize any useful properties of the data or model. \n",
    "\n",
    "# Example: Fitting a line to data\n",
    "\n",
    "As an example of how to compute this thing, let's use the ever-loved example of fitting a straight line to data $y$ at positions $x$ with uncertainties $\\sigma_y$ only on the data $y$. Let's generate some fake data to work with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndata = 8\n",
    "rng = np.random.default_rng(42)\n",
    "true_pars = rng.normal(0, 2.0, size=2)\n",
    "\n",
    "x = np.sort(rng.uniform(0, 10, ndata))\n",
    "yerr = 10 ** rng.uniform(-1, 0, ndata)\n",
    "y = rng.normal(true_pars[0] * x + true_pars[1], yerr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0.5, 0, 'x'), Text(0, 0.5, 'y')]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3AAAANwCAYAAABu6ojXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABYlAAAWJQFJUiTwAAA9k0lEQVR4nO3deZTVdf348dew6yRDhgpOgmKDZO6jZirlEmkuKTpmBzWXcE2tXL4/JzhfwJzGXRQ1y0wrOynOodw1zONJEzUn93UUBWlwAccZmFgE7u8PD/NlAiZgZu7nvmcej3M638/cz+fOfdH3nvg8+WxFuVwuFwAAABS8HlkPAAAAwLoRcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAInolfUArL+tt946mpqaYptttsl6FAAAYD2888470b9//3j33Xc36P0CLkFNTU2xePHirMegA9XV1UVERFlZWcaT0B34vpEvvmvkk+8b+dSe71t79+MFXIJWHnmrra3NeBI6yqhRoyIiYvr06RlPQnfg+0a++K6RT75v5FN7vm/l5eXt+mzXwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACSiKJfL5bIegvVTXl4eERG1tbUZTwIAAKyP9u7LOwIHAACQCAEHAACQCAEHAACQCAGXscbGxjjmmGPi4YcfznoUAACgwAm4jF1//fXxySefZD0GAACQAAGXoSeeeCIee+yxrMcAAAASIeAy0tTUFNddd13WYwAAAAkRcBm5/vrrY6ONNsp6DAAAICECLgNPPvlkPP/88zF27NisRwEAABLSK+sBupumpqa49tpr49xzz41NNtkk63EAAICECLg8u/HGG2OnnXaKfffdN1544YUN/j11dXUxatSoDX7/9OnTN/i9AADQnbVnP7yuri7Kyso2+P0CLo9mzJgRtbW1ccstt2Q9CgAAkCABlycLFiyIyZMnx7nnnhv9+/dv9+8rKytzFA0AADLQnv3w8vLydn22m5jkyY033hg77rhjjBw5MutRAACARDkClwdPPfVUPPvss3HzzTdnPQoAAJAwR+A62cKFC2Py5Mlx9tlnx4ABA7IeBwAASJgjcJ3sueeei969e8c999wT99xzT6t1CxcubFm+88474y9/+UtERGy77bZx1lln5XVOAACg8Am4TjZy5Mi1Xvf2wgsvxAUXXBAREccee2wcdNBB+RwNAABIjIADAAC6vIkTJ65xOTUCDgAA6PImTZrUspxywLmJSYaWLFnSsrx06dIMJwEAAFLgCFye3XTTTVFXVxcREXPmzGl5/dZbb43HHnssIiIGDRoUF154YRbjAQAABUzA5dkZZ5yR9QgAAECinEIJAACQCAEHAACQCAEHAACQCAEHAAB0aY2NjW3+nBIBBwAAdFlVVVVRWlra6rXS0tKoqqrKaKL2cRdKAACgS6qqqorx48ev9npzc3PL6+PGjcv3WO3iCBwAANDlNDY2RnV1dZvbVFdXR1NTU54m6hgCDgAA6HJqamqiubm5zW2am5ujpqYmTxN1DAEHAAB0OXPnzu3Q7QqFgAMAALqcwYMHd+h2hULAAQAAXU5FRUUUFxe3uU1xcXFUVFTkaaKOIeAAAIAup6SkJCorK9vcprKyMvr375+niTqGxwgAAABd0spHBFRXV7e6oUlxcXFUVlYm9wiBiIiiXC6Xy3oI1k95eXlERNTW1mY8CQAAFL6mpqYoKSlp+bmxsTGzI2/t3Zd3CiUAANCl/WespXba5KoEHAAAQCIEHAAAQCIEHAAAQCIEHAAAQCIEHAAAQCIEHAAAQCI8yBsAAOjyJkyYkPUIHULAAQAAXd7EiROzHqFDOIUSAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEb2yHgAAALqriRMnrnEZ1qYol8vlsh6C9VNeXh4REbW1tRlPAgBAexQVFbUs2y3vHtq7L+8USgAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgET0ynoAAAC6jokTJ65xGegYRblcLpf1EKyf8vLyiIiora3NeBIAgNaKiopalu1m/nf+++p+2rsv7xRKAACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARPTKeoDu5Pnnn4+HH344Fi1aFAsWLIj6+vrYcsst48ADD4yDDjooevbsmfWIAABAARNweTJ58uT45JNP4sILL4zi4uKIiGhqaoqLL744rrnmmnjkkUeiqqoqNtpoo4wnBQAACpVTKPPgueeei/vvvz9GjhzZEm8REf37948LLrggevbsGS+99FLceOONGU4JAAAUOgGXB2+++WZERFx77bXx17/+tdW6QYMGxW677RYREY888kg0NzfnfT4AAPKvsbGxzZ9hTQRcHgwbNiwiIhYtWhS/+tWvVltfWloaERHLli2LOXPm5HU2AADyr6qqqmUfcKXS0tKoqqrKaCJS4Rq4PNhjjz3ixz/+cTz55JMxcuTI1davevOSPn365HM0AADyrKqqKsaPH7/a683NzS2vjxs3Lt9jkQgBlyeHHnpoHHrooWtc969//SsiIjbbbLMYMmRIPscCACCPGhsbo7q6us1tqqur45xzzon+/fvnaSpSIuAy1tDQELW1tVFUVBRnnnnmOj9KoK6uLkaNGrXBnzt9+vQNfi8AABumpqbmv97zoLm5OWpqauKUU07J01Ssr/bsh9fV1UVZWdkGv1/A5cnSpUvjZz/7WWyyySZx4oknxhZbbBHNzc1x9dVXxyabbBLnnHNO7LvvvlmPCQBAJ5o7d26Hbkf3I+DyZNmyZTF79uz44IMP4plnnolBgwbFnDlzYpNNNokLL7wwdt999/X6fWVlZY6iAQAkZvDgwR26Hdloz354eXl5uz67KJfL5dr1G9hgS5YsiTvuuCNuv/322HHHHeOiiy6KzTff/L++b+X/02trazt7RACA9VJUVNSybDdzdY2NjVFaWtrmaZTFxcVRX1/vGrguqr378h4jkKG+ffvGiSeeGAcffHC89NJLcd5550VDQ0PWYwEA0ElKSkqisrKyzW0qKyvFG2sl4ArAYYcdFhERH3zwQdx2223ZDgMAQKcaN25cXHLJJVFcXNzq9eLi4rjkkks8QoA2Cbg8eOutt+KUU06J0aNHx0MPPbTa+lUfHfDUU0/lczQAADIwbty4qK+vb/VafX29eOO/EnB5cO2118Z7770XCxcujD/96U9tbrt48eI8TQUA0LEaGxvb/JnW/vM0SadNsi4EXB7MnDkzIiJ69uwZX//611db//7777csDx06NG9zAQB0lKqqqigtLW31WmlpaVRVVWU0EXRNHiOQB5tuumlsvvnm8fOf/zz69u272vonnniiZfmoo47K52gAAO1WVVUV48ePX+315ubmltedGggdwxG4PDjssMPiww8/jBUrVqy2bt68eXH33XdHRMSRRx4Z++23X56nAwDYcI2NjVFdXd3mNtXV1dHU1JSniaBrcwQuD7773e/G8uXLY9y4cbHPPvvE4MGDo3fv3jFz5syYNm1aFBcXx1lnnRUHHHBA1qMCAKyXmpqaNp9pFvHZkbiampo45ZRT8jQVdF0CLg+KiopizJgxMXr06HjxxRdj7ty5sXjx4hgwYED8/Oc/j2233TbrEQEANsjcuXM7dDugbQIujzbaaKP46le/mvUYAAAdZvDgwR26HdA218ABALDBKioqVnsg9X8qLi6OioqKPE0EXZuAAwBgg5WUlERlZWWb21RWVnrGGXQQp1ACANAuKx8RUF1d3eqGJsXFxVFZWekRAtCBinK5XC7rIVg/5eXlERFRW1ub8SQAAP+nqakpSkpKWn5ubGxs15G3iRMnrnG5KykqKmpZtlvePbR3X17AJUjAAQCFqiODpDvETXf4M9Jae/flXQMHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAACQCAEHAEDBaWxsbPNn6K4EHAAABaWqqipKS0tbvVZaWhpVVVUZTQSFo1fWAwAAwEpVVVUxfvz41V5vbm5ueX3cuHH5HgsKhiNwAAAUhMbGxqiurm5zm+rq6mhqasrTRFB4BBwAAAWhpqYmmpub29ymubk5ampq8jQRFB4BBwBAQZg7d26HbgddkYADAKAgDB48uEO3g65IwAEAUBAqKiqiuLi4zW2Ki4ujoqIiTxNB4RFwAAAUhJKSkqisrGxzm8rKyujfv3+eJoLC4zECAAAUjJWPCKiurm51Q5Pi4uKorKz0CAG6PQEHAEBBGTduXJxzzjlRUlLS8lp9fX2XPPI2YcKErEcgMQIOAICC85+x1hXjLSJi4sSJWY9AYgQcAAAdxhEl6FwCDgCADuOIEnQud6EEAABIhIADAABIhIADAABIhIADAABIhIADAABIhIADAABIhIADAABIhIADAABIhIADAABIhIADAABIRK+sBwAAgDWZMGFC1iNAwRFwAAAUpIkTJ2Y9AhQcp1ACAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkQsABAAAkolfWAwAAEDFx4sQ1LgOsqiiXy+WyHoL1U15eHhERtbW1GU8CAHSUoqKilmW7Z9B1tXdf3imUAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAieiV9QDdyTvvvBMPPvhgzJ8/P+bNmxcffvhhDBw4MPbaa684/PDDo3///lmPCAAAFDABlyf33XdfzJkzJ0488cQoLi6OiIjm5ua49dZb47bbbotp06ZFZWVl7L777hlPCgAAFCqnUObBq6++Gq+//nqcccYZLfEWEVFcXBxnn3127LPPPtHU1BSTJk2KWbNmZTgpAABQyARcHkybNi1Gjhy51vVHHXVUREQsXrw47rjjjnyNBQAAJEbA5cGsWbPi6aefjiVLlqxx/bbbbtuy/OKLL+ZrLAAAIDECLg823njjuPfee+OYY46Jm2++ebX1vXr936WIDQ0N+RwNAABIiIDLg2984xsREbFo0aKYOnVqzJ07t9X6efPmtSyXlpbmdTYAACAd7kKZB0cccUQ0NjbGY489FltttVUMHDiw1fp//OMfLcujR49ep99ZV1cXo0aN2uCZpk+fvsHvBQCA7qw9++F1dXVRVla2we8XcHnQs2fPOPnkk+Pkk09ebd3ixYvjT3/6U0REHHzwwXHIIYfkezwAACARAi5DuVwurrvuuqivr4/Ro0fHmWeeuc7vLSsrcxQNAAAy0J798PLy8nZ9toDLyLJly+Kqq66KJ554Ii666KI48MADsx4JAAAocAIuA01NTXHxxRfHggUL4oYbboghQ4ZkPRIAAJAAd6HMs9mzZ8c555wTQ4YMiSlTpog3AABgnQm4PHrttddi/PjxMXbs2Dj33HOjT58+rdavWLEiHn300YymAwAACp1TKPNkzpw5UVVVFZMmTYptt912jdvU19fHnDlz8jwZAACQCkfg8uTKK6+Mb33rW2uNt4jWz4MDAAD4T47A5cHMmTPjlVdeiU033TSeeeaZKCoqarV+xYoVMXv27PjDH/4QRx11VEZTAgAAhU7A5UFtbW1ERDz++OPx+OOPZzwNAACQKqdQ5sEHH3yQ9QgAAEAXUJTL5XJZD8H6Wfn09pVH9gCA9K16iYXdM+i62rsv7wgcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAABAIgQcAEDGGhsb2/wZYCUBBwCQoaqqqigtLW31WmlpaVRVVWU0EVDIemU9AABAd1VVVRXjx49f7fXm5uaW18eNG5fvsYAC5ggcAEAGGhsbo7q6us1tqquro6mpKU8TASkQcAAAGaipqYnm5uY2t2lubo6ampo8TQSkQMABAGRg7ty5Hbod0D0IOACADAwePLhDtwO6BwEHAJCBioqKKC4ubnOb4uLiqKioyNNEQAoEHABABkpKSqKysrLNbSorK6N///55mghIgccIAABkZOUjAqqrq1vd0KS4uDgqKys9QgBYTVEul8tlPQTrp7y8PCIiamtrM54EAOgITU1NUVJS0vJzY2OjI2/QRbV3X94plAAAGfvPWBNvwNoIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIuAxNmTIlfve732U9BgAAkAgBl5EPP/wwHn300azHAAAAEtIr6wG6m1wuF6+88kpcffXVsXDhwqzHAQAKxIQJE7IeAUiAgMuTmpqamDFjRjQ2NkZzc3PMmzcv65EAgAIyceLErEcAEiDg8qSioiIqKioiIuL999+PE044IeOJAACA1Ag4AArGqkcgHI0AgNUJOAAKxqRJk1qWBRwArE7AJaquri5GjRq1we+fPn16B04DAADdR3v2w+vq6qKsrGyD3+8xAgAAAIlwBC5RZWVljqIBAEAG2rMfXl5e3q7PdgQOAAAgEQIOAAAgEQIOAAAgEZ0ecC+//HJnfwQAAEC30OkBt8cee3T2RwAAAHQLnR5wS5YsiXfeeaezPwYAAKDLy8s1cGeffXYsWbIkHx8FAADQZeUl4B566KHYbrvt4rLLLos333wzHx9Z0FaN2aVLl2Y4CQAAkJK8PMj7c5/7XHz5y1+O3//+91FZWRkjRoyII488Mo488sjYc8898zFC5mpqamLGjBkRETFv3ryW1+++++549dVXo6ioKHr16hWXXXZZViMCAAAFriiXy+U68wN69OgRL730UnzlK1+JiIj3338/7rvvvrj33nvjkUceic9//vPxne98J4488sg44IADolevvDRl0lY+vb22tjbjSQA6VlFRUctyJ//1BACZaO++fKefQtm3b9+WeIuIGDRoUIwdOzbuvvvumD9/ftx0002xYsWKOOWUU2KzzTaLMWPGxNSpU2PBggWdPRoAAEBSOj3gZs2atdZ1/fr1i8MOOyxuuummmDNnTkyfPj2GDx8el112WWy++ebx7W9/O375y1/G+++/39ljAgAAFLxOP4VyQz366KNx+umnx8yZM6OoqCj22GOPGD16dBxxxBGx3XbbZT1eppxCCXRVTqEEoKtr7758QQXcO++8E3feeWdMnTo1XnjhhZbXVx2xqKgohg8fHkceeWSMHj2629wEZVUCDuiqBBwAXV3BXwP3ve99r831s2fPjiuvvDL22GOP+NKXvhTjxo2LF154IXK5XMt/RowYERMmTIhXX3013nrrrTjttNPiqaeein322Sd22GGHqKmp6ew/BgAAQOY6/Qhcv379YvHixa1eq6+vj6lTp8add94ZzzzzTMvrq44ybNiwOPbYY+PYY4+NnXbaaY2/e968eXHzzTdHdXV1nHbaaXHllVd2zh+iwDgCB3RVjsAB0NW1d1++0+/Zv3Tp0rj99ttj+PDh8dRTT8Vdd90VM2bMaPmLedW/oIcMGRLf/e5349hjj235g7Vl4MCBUVlZGUcffXTsvffeMWLEiBg7dmyn/VkAAACylJeHrp144omtfl412gYPHhzHHHNMHHvssfG1r31tg37/8OHD48ILL4wpU6YIOAAAoMvK21OzV422zTffPI4++ug49thjY+TIka1OmVlXn376adx5553x8ccfxxFHHBEDBw6MN998syNHBgAAKCh5CbhcLhcbb7xxfO9734vvfe97ccABB0SPHu27f8oPf/jDuOWWWyIi4ic/+Un06NEjtthii44YFwAAoCDlJeCOOeaY+MUvfhGbbrpph/3OF198MXK5XBQVFUUul4vly5evdqomAABAV9LpjxEoKiqKm2++uUPjLSLijDPOiB49erScmnnqqafGxRdf3KGfAQAAUEg6/Qhcnz59on///h3+e0866aTYa6+94sUXX4ztttsudt555w7/DAAAgELS6QH3xz/+sdN+94gRI2LEiBGd9vsBAAAKSaefQnnkkUd29kcAAAB0C50ecAAAAHQMAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQdAQWhsbGzzZwBAwAFQAKqqqqK0tLTVa6WlpVFVVZXRRABQmHplPQAA3VtVVVWMHz9+tdebm5tbXh83bly+xwKAguQIHACZaWxsjOrq6ja3qa6ujqampjxNBACFTcABkJmamppobm5uc5vm5uaoqanJ00QAUNgEHACZmTt3boduBwBdnYADIDODBw/u0O0AoKsTcABkpqKiIoqLi9vcpri4OCoqKvI0EQAUNgEHQGZKSkqisrKyzW0qKyujf//+eZoIAAqbxwgAkKmVjwiorq5udUOT4uLiqKys9AgBAFhFUS6Xy2U9BOunvLw8IiJqa2szngSg4zQ1NUVJSUnLz42NjY68AdDltHdf3imUABSE/4w18QYAqxNwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAiRBwAAAAieiV9QAAsNKECROyHgEACpqAA6BgTJw4MesRAKCgOYUSAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEZ4Dl2evvPJK3HHHHdHQ0BA9e/aMiIj9998/Dj/88JafAQAA1kTA5dGf//znuOWWW+J///d/Y4899oiIiPnz58ekSZPi73//e1RVVUWfPn0ynhIAAChUTqHMkyeeeCJ+8YtfxHHHHdcSbxERX/jCF2LcuHHx4osvxqWXXhq5XC7DKQEAgEIm4PKgubk5rr766sjlcnHwwQevtn6LLbaIXXfdNR5//PH461//msGEAABACgRcHkybNi0WLFgQQ4YMiQEDBqxxm5122ikiIm6//fZYvnx5HqcDAABSIeDy4P7774+IiKFDh651m2222SYiIv71r3/F888/n4+xAACAxAi4TvbOO+/E/PnzIyJi0KBBa91u8ODBLcvPPvtsp88FAACkx10oO9lLL73Usry20yf/c93LL7/8X39vXV1djBo1aoPnmj59+ga/FwAAurP27IfX1dVFWVnZBr/fEbhOtvLoW0TERhtttNbtNt544zW+BwAAYCVH4DpZY2Njy3KvXmv/r3vVdZ988sl//b1lZWWOogEAQAbasx9eXl7ers92BK6T9e3bt2V5xYoVa91u1XX9+vXr1JkAAIA0CbhO9vnPf75l+dNPP13rdkuXLm1ZbutaOQAAoPsScJ3si1/8YstyU1PTWrdbdd1WW23VqTMBAABpEnCdbNddd42ePXtGRNs3J1l1XXvPiwUAALomAdfJiouLY+edd46Izx7SvTbvvfdeRET06NEj9tprr7zMBgAApEXA5cGYMWMi4rNnPqztRiavv/56RHz2TInNN988b7MBAADpEHB5sPPOO8dee+0Vzc3NrR7svdLy5cvjmWeeieLi4jj++OMzmBAAAEiBgMuTiy66KLbeeuv4/e9/H8uXL2+17oEHHoiGhob43//93xg0aFBGEwIAAIWuKJfL5bIeorv4+OOP4/LLL4+lS5fGoYceGhtvvHE899xzMWPGjDj33HNjjz32WKffs/ImJ7W1tZ05LgAA0MHauy8v4DLw5ptvxssvvxxLly6NIUOGxO677x59+vRZ5/cLOAAASFN79+V7deQwrJvhw4fH8OHDsx4DAABIjGvgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEiHgAAAAEtEr6wFgXUycOHGNywAA0J0U5XK5XNZDsH7Ky8sjIqK2tjbjSfKnqKioZdlXFgCAVLV3X94plAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcAAAAIkQcBl544034vjjj896DAAAICECLiN33HFH1iMAAACJEXB51tDQEFdeeWU88cQTWY8CAAAkplfWA3QHs2bNiuuuuy6WLl0an3zyScyfPz/rkQAAgAQJuDwYOnRoXHXVVS0/n3/++fHiiy9mOBEAAJAip1ACAAAkQsABAAAkwimUiaqrq4tRo0Zt8PunT5/egdMAAED30Z798Lq6uigrK9vg9zsCBwAAkAhH4BJVVlbmKBoAAGSgPfvh5eXl7fpsR+AAAAAS0W2PwF1++eWdegTL0TEAAKCjdduAO/300+P73/9+1mMAAACss24bcCUlJVFSUpL1GAAAAOvMNXAAAACJEHAAAACJ6LanUNJ+EydOXOMyAADQOQRcBpYuXdrq/6Zq0qRJLcsCDgAAOp+Ay4NZs2bFddddFxERy5YtizfffDMiIhoaGuLMM8+MjTfeOCIijj/++Nh1110zmxMAAChsAi4Phg4dGldddVXWYwAAAIlzExMAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDg2SGNjY5s/AwAAHU/Asd6qqqqitLS01WulpaVRVVWV0UQAANA99Mp6ANJSVVUV48ePX+315ubmltfHjRuX77EAAKBbcASOddbY2BjV1dVtblNdXR1NTU15mggAALoXAcc6q6mpiebm5ja3aW5ujpqamjxNBAAA3YuAY53NnTu3Q7dbV26YAgAAnxFwrLPBgwd36Hbrwg1TAADg/wg41llFRUUUFxe3uU1xcXFUVFR0yOetvGHKf562ufKGKSIOAIDuRsCxzkpKSqKysrLNbSorK6N///7t/iw3TAEAgNUJONbLuHHj4pJLLlntSFxxcXFccsklHfYIATdMAQCA1XkOHOtt3Lhxcc4550RJSUnLa/X19R1y5G2lrG6YAgAAhcwRODbIf8ZaR8ZbRDY3TAEAgEIn4ChI+b5hCgAApEDAUZDyecMUAABIhWvgKFgrb4hSXV3d6oYmxcXFUVlZ2WE3TAEAgFQU5XK5XNZDsH7Ky8sjIqK2tjbTOYqKilqWO/Nr1NTU1OqGKY2NjY68AQCQpPbuyzuFkoLX2TdMAQCAVAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARAg4AACARPTKegDSNWHChKxHAACAbkXAscEmTpyY9QgAANCtOIUSAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEb2yHqC7WLRoUdxzzz0xc+bMWLBgQcyZMyeKiopil112iYMPPji+/OUvZz0iAABQ4ARcHsydOzd++ctfxvHHHx/HHntsy+tPPPFEXHXVVfHggw/GEUccEWeddVYUFRVlOCkAAFDInELZyVasWBFXXHFFnHfeefGlL32p1bp99903LrzwwsjlcvHnP/85br/99oymBAAAUiDgOtmzzz4bAwYMiP79+69x/d577x2DBw+OiIiamppobm7O53gAAEBCBFwnmzVrVrz99ttRX1+/1m2GDRsWERH//ve/o66uLl+jAQAAiRFwnay4uDjq6+vjpJNOip/85Cfx0UcfrbZNr17/dyliQ0NDPscDAAASIuA62Z577hn9+vWLXC4XL7/8cjzwwAOrbTNv3ryW5dLS0nyOBwAAJMRdKDvZwIEDY9KkSXHrrbfGwoULY4899mi1vqmpKd58882IiNhhhx1i+PDh6/R76+rqYtSoURs81/Tp0zf4vQAA0J21Zz+8rq4uysrKNvj9Ai4Pdtttt9htt93WuO6uu+6KTz/9NAYNGhSVlZV5ngwAAEiJgMvQs88+G1OnTo2hQ4dGdXV1bLbZZuv83rKyMkfRAAAgA+3ZDy8vL2/XZwu4jDz55JNxySWXxP777x8/+tGPYqONNsp6JAAAoMB124C7/PLLO/UIVlu/e+rUqXH77bfH2WefHYccckinzQAAAHQt3TbgTj/99Pj+97+f189ctmxZTJ48Od54442YMmVKDB06NK+fDwAApK3bBlxJSUmUlJTk7fMWL14cEyZMiEGDBsX1118fffv2XW2bRx55JL75zW/mbSYAACAtngOXB7lcLqqqqmL77bePn/zkJ2uMt4iIZ555Js+TAQAAKRFwefDggw/Gu+++GyeccMJat2lsbGx5HhwAAMCadNtTKPPp3nvvjaVLl8bf//736Nev32rrFyxYEHfffXcsW7Ysg+kAAIBUCLhO1tjYGG+99VZERFx88cVtbrvFFlvkYyQAACBRTqHsZB988EHWIwAAAF2EI3CdbPjw4Z36vDkAAKD7cAQOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEQIOAAAgEb2yHgDWxYQJE7IeAQAAMifgSMLEiROzHgEAADLnFEoAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgoAKNGjYpRo0ZlPQbdhO8b+eK7Rj75vpFPWX7fBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiBBwAAEAiinK5XC7rIVg/m266aSxevDi+/OUvZz0KHaSuri4iIsrKyjKehO7A94188V0jn3zfyKf2fN9ee+216NevX3z88ccb9NkCLkFbb711NDU1xTbbbJP1KAAAwHp45513on///vHuu+9u0PsFHAAAQCJcAwcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJAIAQcAAJCIXlkPAN3VO++8Ew8++GDMnz8/5s2bFx9++GEMHDgw9tprrzj88MOjf//+WY9IF9fY2Bhjx46NsWPHxkEHHZT1OADAOhBwkIH77rsv5syZEyeeeGIUFxdHRERzc3Pceuutcdttt8W0adOisrIydt9994wnpSu7/vrr45NPPsl6DLqwmTNnxl133RX19fXRq1evWL58eeyyyy4xZsyY6NOnT9bj0UV88skncc8998Q//vGPiIjo3bt3NDc3x4477hiHH354DB06NOMJSdWUKVOipKQkvv/976/T9o2NjXHHHXfECy+8EH369ImlS5fG8OHD47jjjovNNtusw+ZyCiXk2auvvhqvv/56nHHGGS3xFhFRXFwcZ599duyzzz7R1NQUkyZNilmzZmU4KV3ZE088EY899ljWY9CF3XvvvVFZWRn7779/XHvttXHVVVfF1VdfHfPmzYtx48ZFLpfLekS6gCeffDJOOeWUWLx4cVx66aUxZcqUuPrqq+OKK66IFStWxJlnnhkPPPBA1mOSoA8//DAeffTRdd5+zpw5ccYZZ8SSJUti8uTJMXny5JgyZUoMHDgwTj/99Hj11Vc7bDYBB3k2bdq0GDly5FrXH3XUURERsXjx4rjjjjvyNRbdSFNTU1x33XVZj0EXNn369JgyZUpcdNFFseeee7a8nsvlYtGiRfH888/HG2+8keGEdAUvvfRS/OxnP4vvfOc7cdppp7X6R9H+/fvHueeeG/vtt19ce+21HbrzTNeWy+Xi5ZdfjosuuigWLly4Tu9pamqKcePGxec+97k4++yzW84w6NmzZxx//PFRVlYWEyZMiLlz53bIjAIO8mzWrFnx9NNPx5IlS9a4ftttt21ZfvHFF/M1Ft3I9ddfHxtttFHWY9BFNTQ0xI033hjl5eWx6667tlp30003xd/+9reIiJg/f34W49GF/OpXv4ply5bFYYcdttZtTjrppFixYkVMmzYtj5ORopqamjj//PPj1FNPjaqqqnjvvffW+b2/+tWvor6+Pr71rW9Fjx6r59VBBx0Un3zySVxzzTUdMquAgzzbeOON4957741jjjkmbr755tXW9+r1f5emNjQ05HM0uoEnn3wynn/++Rg7dmzWo9BF3XnnnbFw4cLYb7/9Vlu3YsWKiIgYMmRI7LLLLvkdjC5l0aJFLUdxN9lkk7Vut/nmm8eAAQNcksB/VVFREVdddVX8+te/Xq/Qqq+vj0ceeSQiInbaaac1brPy9eeee65D/nFewEGefeMb34iIz/7ymTp16mqH0+fNm9eyXFpamtfZ6Nqampri2muvjXPPPdddTukUy5cvb7lm5Ctf+cpq63/4wx/GLbfcEjfddFOr091gfa1YsaLlOsrXX3+9zW2XLFkSm266aT7Goht66KGHYvny5VFUVBRDhgxZ4zYDBw5s+YeG+++/v92fKeAgz4444ogYM2ZMbLnllvHVr341Bg4c2Gr9yrtoRUSMHj063+PRhd14442x0047xb777pv1KHRRr7/+ejQ0NETPnj1j8ODBq63v0aNHDBkyJHr37p3BdHQlxcXFMWjQoIiImDx5csyZM2eN282cOTMWLVoU3/zmN/M5Ht3Is88+GxERJSUlbV6esPJ/E//5z3+2+yZOHiMAedazZ884+eST4+STT15t3eLFi+NPf/pTREQcfPDBccghh+R7PLqoGTNmRG1tbdxyyy1Zj0IXtvKUtpKSkujZs2fMmjUr7rzzznj//fejV69eMXTo0DjqqKPWGHewvo4++ui44YYbYs6cOXHqqafGgQceGAcffHDssMMOEfHZUbpf//rXcfDBB8eoUaMynpauaNGiRfHWW29FRMSAAQPa3Hbl+k8++SRmz57drsdbCDgoELlcLq677rqor6+P0aNHx5lnnpn1SHQRCxYsiMmTJzt1kk737rvvRsRnz+GaMWNG3HzzzXHBBRfE9ttvH83NzXH55ZfHKaecEv/v//2/NV4jB+vjyCOPjPfeey/uueeeWLZsWTz88MPx8MMPx5ZbbhmjRo2Kt99+O/bee+82b3IC7fHxxx+3HE37bzcH23jjjVuW58+f366AcwolFIBly5bF5ZdfHo8//nhcdNFFcdZZZ0VRUVHWY9FF3HjjjbHjjju2+fgK6AhNTU0R8dkNmK644or46U9/Gttvv31EfHbK2/nnnx9FRUVx+eWXx7/+9a8sR6WLOOecc+Kyyy5rdU1lfX19/Pa3v40nn3wynn766Xj77bcznJCurLGxsWV51ZvQrcmq61d934YQcJCxpqamuOiii2LmzJlxww03xIEHHpj1SHQhTz31VDz77LNx9tlnZz0K3cCiRYsiImLp0qWxww47xJe+9KVW6/v37x877rhjfPrpp55zSYf4+9//HjfccENceOGF8fvf/z5OOumk2GqrrSLis1Mon3rqqTjrrLPi4YcfznhSuqK+ffu2LK+8y+7arLp+1fdtCAEHGZo9e3acc845MWTIkJgyZcpa714EG2LhwoUxefLkOPvss//rufnQEVY9hWj33Xdf4zYjRoyIiM/+cQHa45577omJEyfGD37wg9hnn31i0KBBcdxxx8VvfvObmDJlSstdn1esWBHXXHPNWm90Ahvq85//fMvyp59+2ua2S5cuXeP7NoRr4CAjr732WlRXV8epp566xlPbVqxYEY899lgccMABGUxHV/Dcc89F796945577ol77rmn1bqFCxe2LN95553xl7/8JSI+e5D8WWedldc56TpWvcZyiy22WOM2JSUlEfHZhfyLFi3yUHk2yIIFC+Lmm2+OLbbYIr72ta+ttn7EiBExfvz4+NrXvhaXXnppLF++PB577LE4/vjjM5iWrmrAgAGxySabxIIFC1pOIV+blet79OjR7sdECTjIwJw5c6KqqiomTZoU22677Rq3qa+v96+FtMvIkSPXet3bCy+8EBdccEFERBx77LFx0EEH5XM0uqgvfvGLLct9+vRZ4zarvt6jhxOB2DCvvPJKLF68OHbcccc2rxk/8MAD4+GHH47nnnsuPv744zxOSHfQo0eP2GWXXeLxxx+PhoaGyOVya/0+zp8/PyIiysrK2n1DMf/LCRm48sor41vf+tZa4y2i9fPgAFKw8847tyyv7V+jGxoaIuKzB9u29zoQuq8lS5ZERLQ8HLktw4YNi4iIL3zhC506E93TPvvsExGffSc/+uijNW7z6aefxvvvvx8REXvvvXe7P9MROMizmTNnxiuvvBKbbrppPPPMM6v9S82KFSti9uzZ8Yc//CGOOuqojKYEWH9lZWUxaNCgeP/99+PDDz9c4zYrzyxY02lvsK5WPktw5T8ItGVl7K3tukxoj/322y9+97vfRX19fbzxxhux+eabr7ZNXV1dLF++PDbZZJM44ogj2v2ZAg7yrLa2NiIiHn/88Xj88ccznobuauUOTUTrC6uhPXr06BHHHXdcXHXVVfHMM8/EMccc02r9smXL4plnnonevXvH0UcfndGUdAVlZWUxbNiweOWVV6KhoWGtN4VYsmRJPPnkk7HLLrvEdtttl+cp6Q569uwZY8eOjYsvvjieeuqpNV66sPKmTccff3yrR15sKAEHefbBBx9kPQLd1E033RR1dXUREa2ur7z11lvjsccei4iIQYMGxYUXXpjFeHQRo0aNihkzZsTTTz8ds2bNavWw2oceeigWLlwYZ511Vrsv4qd7KyoqivPPPz/+53/+J6qrq2PChAmr7RgvX748rr322ujTp09UVlZmNCkpWt9/5Bw5cmQcd9xxcdddd0VFRUVss802Les++uijuP/+++Oggw7qsDOrinIrHx8OANABli5dGldffXW88MILMWbMmBg0aFC88sorcf/998fJJ58chxxySNYj0kXMmjUrrr/++pg9e3Z8+9vfjhEjRkRRUVHU19fHo48+GkOHDo0zzjgjPve5z2U9KgWupqYmZsyYERER8+bNi/r6+oiI6NevX5SVlUVRUVH06tUrLrvssrX+jj/+8Y9x3333xRFHHBFbb711zJ07N6ZNmxZ77713jB07Nnr27Nkhswo4AKBTzJ49O55//vn497//HVtuuWXstttudqTpFHPnzo2XXnopPv744+jXr18MHDgwdtllF9838q6hoSH+8Y9/xLx582LAgAGx++67r/G6uPYQcAAAAInwGAEAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAAIBECDgAyLNJkyZF7969o6ioaLX/9O7dO2677baIiHjwwQejb9++q23Tq1evGDNmTLZ/CAAyUZTL5XJZDwEA3dG0adPimGOOiRUrVkRExNe//vV45JFHonfv3q22O++88+Kaa66JrbbaKn7729/G17/+9ejZs2cWIwOQMUfgACAjRx11VIwdO7bl54aGhtXiLSJi0KBB0a9fv3jggQdi//33F28A3ZgjcACQofnz58fw4cPj448/joiIO++8M7773e+2rF+8eHHsuOOOcdppp8WFF16Y1ZgAFAhH4AAgQ1/4wheiqqqq5ecLLrgg/v3vf7f8PGnSpBgwYECcd955WYwHQIERcACQsdNOOy3Ky8sjIuK9996Ln//85xER8fzzz8d1110Xv/nNb5w2CUBEOIUSAArCU089FXvvvXfkcrno27dvvPjiizFmzJg49NBDY9KkSVmPB0CBEHAAUCB+8IMfxG9+85uIiNh8881js802i3/+85/Rp0+fjCcDoFA4hRIACsSll14aAwYMiIiIDz/8MM444wzxBkArAg4ACsRmm20WP/vZz1p+vvLKK2PRokUZTgRAoRFwAFBANt5445blWbNmRXV1dYbTAFBoXAMHAAXigw8+iO233z7OPffcmDRpUssNTV599dUYNmxY1uMBUAAcgQOAAvGjH/0oDj744JgwYUKccsopERGxZMmS+NGPfpTxZAAUCkfgAKAA3H///XHSSSfFq6++GptttlnMmzcvhg8fHg0NDRERce+998Zhhx2W8ZQAZM0ROADI2MKFC+PMM8+Mq6++OjbbbLOIiBg4cGBccsklLdv8+Mc/jiVLlmQ1IgAFQsABQMZ++tOfxrBhw+KEE05o9foZZ5wRu+22W0REvP3223HFFVdkMR4ABcQplACQoenTp8dBBx0UU6dOjYqKitXW/+1vf4tvfOMbERGx0UYbxWuvvRZDhw7N95gAFAhH4AAgA8uWLYv77rsvjjnmmMjlcjF79uw1bjds2LAoKiqKiIhFixbFCSecEB999FE+RwWggDgCBwB5NmbMmJg6dWosX7681eubbrppzJ8/PyI+e6TA0KFDY+nSpbGmv6r79u0bf/vb32LPPffMy8wAFAYBBwAAkAinUAIAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACRCwAEAACTi/wP1wxK5kTZPQwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 440,
       "width": 440
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(x, y, yerr, marker=\"o\", ls=\"none\")\n",
    "ax.set(xlabel=\"x\", ylabel=\"y\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, our model for the data is:\n",
    "$$\n",
    "f(x \\,;\\, \\theta) = a \\, x + b \\\\\n",
    "\\theta = (a, b)\n",
    "$$\n",
    "and we will assume that the data are generated from some process such that the uncertainties $\\sigma_y$ are Gaussian (and I'll drop the subscript $y$ because I'm lazy). Our likelihood function for a single data point $(x_i, y_i, \\sigma_i)$ is then:\n",
    "$$\n",
    "L(y_i \\,;\\, \\theta) = \\mathcal{N}(y_i \\,|\\, f(x \\,;\\, \\theta), \\sigma_i)\n",
    "$$\n",
    "($x_i$ and $\\sigma_i$ aren't really data because we assume we know them perfectly, but they are associated with a given data point $y$ as metadata). \n",
    "$\\mathcal{N}(w \\,|\\, \\mu, \\sigma)$ is the Normal distribution over $w$ with mean $\\mu$ and standard deviation $\\sigma$.\n",
    "\n",
    "In this model, the data points are independent, so the total likelihood $\\mathcal{L}$ is just the product of the likelihoods for each individual data point:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathcal{L} &= \\prod_i L(y_i \\,;\\, \\theta)\\\\ \n",
    "&= \\prod_i \\mathcal{N}(y_i \\,|\\, f_i, \\sigma_i)\n",
    "\\end{align}\n",
    "$$\n",
    "where, as a reminder, $f_i$ is the model predicted value of $y$ at a given location $x_i$ -- you could think of it as the model's prediction for the true value of $y_i$.\n",
    "In the log, this product is just a sum:\n",
    "$$\n",
    "\\ln\\mathcal{L} = \\sum_i \\ln\\mathcal{N}(y_i \\,|\\, f_i, \\sigma_i)\n",
    "$$\n",
    "\n",
    "So now we need to take second derivatives of this thing. The derivative here can be moved inside of the sum, so  \n",
    "$$\n",
    "\\frac{\\partial^2\\ln\\mathcal{L}}{\\partial \\theta_j \\partial \\theta_k} = \\sum_i \\frac{\\partial^2}{\\partial \\theta_j \\partial \\theta_k}\\ln\\mathcal{N}(y_i \\,|\\, f_i, \\sigma_i)\n",
    "$$\n",
    "\n",
    "To see how this works in this one-dimensional case, let's expand the expression for the log-Normal distribution:\n",
    "$$\n",
    "\\ln\\mathcal{N}(y_i \\,|\\, f_i, \\sigma_i) = -\\frac{1}{2} \\left[ \\frac{(y_i - f_i)^2}{\\sigma_i^2} + \\ln(2\\pi\\,\\sigma_i^2) \\right]\n",
    "$$\n",
    "(this may seem like a weird way to write it but trust me). Now, we need to take derivatives of this expression with respect to the model parameters $\\theta$. One thing to note is that the second term $\\ln(2\\pi\\,\\sigma^2)$ does not depend on the parameters (it is only a function of the data, i.e. the uncertainties), so the derivatives of this term are zero. The first term is a quadratic function of $y_i$ and $f_i$, but $y_i$ is data and $f_i$ is a function of the model parameters, so the only terms that survive the derivative are those that depend on $f_i$ (i.e. the model predictions):\n",
    "$$\n",
    "\\frac{\\partial^2}{\\partial \\theta_j \\partial \\theta_k} \n",
    "    \\ln\\mathcal{N}(y_i \\,|\\, f_i, \\sigma_i) = \n",
    "        \\frac{1}{\\sigma_i^2} \\frac{\\mathrm{d} f_i}{\\mathrm{d}\\theta_j} \n",
    "            \\frac{\\mathrm{d} f_i}{\\mathrm{d}\\theta_k}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model $f$ only has two parameters, so we can write this matrix out explicitly, for one data point:\n",
    "$$\n",
    "\\begin{align}\n",
    "&= \\frac{1}{\\sigma_i^2} \\, \\begin{pmatrix}\n",
    "\\frac{\\mathrm{d} f_i}{\\mathrm{d}a} \\frac{\\mathrm{d} f_i}{\\mathrm{d}a} & \\frac{\\mathrm{d} f_i}{\\mathrm{d}a} \\frac{\\mathrm{d} f_i}{\\mathrm{d}b} \\\\\n",
    "\\frac{\\mathrm{d} f_i}{\\mathrm{d}b} \\frac{\\mathrm{d} f_i}{\\mathrm{d}a} & \\frac{\\mathrm{d} f_i}{\\mathrm{d}b} \\frac{\\mathrm{d} f_i}{\\mathrm{d}b}\n",
    "\\end{pmatrix} \\\\\n",
    "&= \\frac{1}{\\sigma_i^2} \\, \\begin{pmatrix}\n",
    "x_i^2 & x_i \\\\\n",
    "x_i & 1 \n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "And so the Fisher information matrix is just the sum of these matrices over all data points:\n",
    "$$\n",
    "\\mathcal{F} = \\sum_i \\frac{1}{\\sigma_i^2} \\, \\begin{pmatrix}\n",
    "x_i^2 & x_i \\\\\n",
    "x_i & 1\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an interesting (but simple) result: the Fisher information here does not even contain the model parameter values! That is because our model $f$ is _linear_ in the model parameters and our likelihood is Gaussian, so the curvature is a constant. Let's look at how to compute this given our toy data above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2, 8)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is the matrix component of the Fisher information from above\n",
    "M = np.array([[x**2, x], [x, np.ones_like(x)]])\n",
    "M.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = np.sum(1 / yerr**2 * M, axis=-1)\n",
    "F.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10822.60583843,  1226.64589812],\n",
       "       [ 1226.64589812,   157.35632654]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK so we've computed the Fisher information! ...So what? It turns out that the inverse of the Fisher information matrix is a useful matrix itself: it is the covariance matrix of the maximum likelihood parameter estimates. This is a useful property because it tells us the expected precision of our parameter estimates given the data we have. If the Fisher information is large, then the parameter uncertainties (i.e. the inverses) are small, so the parameters are well-constrained by the data. Let's look at the inverse of the Fisher information matrix for our toy data above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00079335, -0.00618439],\n",
       "       [-0.00618439,  0.05456446]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Finv = np.linalg.inv(F)\n",
    "Finv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we ignore the covariances (i.e. the off-diagonal terms), this tells us the expected errors we would get on the parameter values if we fit the data with a straight line. The diagonal terms are the variances of the parameter estimates, and the square roots of these are the expected standard deviations of the parameter estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.02816638, 0.23359036])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(np.diag(Finv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now do a maximum likelihood fit using least-squares algebra and compare the results to the Fisher information matrix. With least-squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Design matrix:\n",
    "A = np.stack((x, np.ones_like(x)), axis=-1)\n",
    "\n",
    "# Data covariance matrix:\n",
    "Cinv = np.diag(1 / yerr**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the standard matrix algebra for computing the maximum likelihood parameter estimates for a linear model (`mle_pars` below). In this case, the covariance matrix of the parameter estimates, $P$, is the inverse of the Fisher information matrix! (BTW, as an aside, it's generally bad to construct the data covariance matrix this way and to do the matrix inverse explicitly here — a better way is to use, e.g., `numpy.linalg.leastsq` or `numpy.linalg.solve`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.58499692, -1.87999557])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = np.linalg.inv(A.T @ Cinv @ A)\n",
    "mle_pars = P @ A.T @ Cinv @ y\n",
    "mle_pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.02816638, 0.23359036])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Note that these are the same values as sqrt(diag(Finv))!\n",
    "np.sqrt(np.diag(P))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More generally: Multi-dimensional and nonlinear (but still Gaussian) models\n",
    "\n",
    "The example above is a special case where we can solve for the Fisher information analytically and our model is linear. We are often in more complex situations in astronomy. One example is our project (where our likelihood function involves orbit integrations; here the model is nonlinear and the data are multi-dimensional). Another example is solving the Kepler problem by fitting a two-body orbit to, e.g., radial velocity data (model is nonlinear but data are still one dimensional). \n",
    "\n",
    "In general, our data may be a vector $\\boldsymbol{y}$ and our model may be a function $\\boldsymbol{f}(\\boldsymbol{x} \\,;\\, \\boldsymbol{\\theta})$ where $\\boldsymbol{x}$ is a vector of independent variables and $\\boldsymbol{\\theta}$ is a vector of model parameters. If the likelihood is non-Gaussian, there's not much we can do analytically. But if the likelihood is Gaussian, even with a nonlinear model, we can simplify the Fisher information matrix calculation by using the chain rule for derivatives.\n",
    "\n",
    "In this case, our likelihood for one data point $\\boldsymbol{y}_i$ (note the vector $y$) is:\n",
    "$$\n",
    "L(\\boldsymbol{y}_i \\,;\\, \\boldsymbol{\\theta}) = \\mathcal{N}(\\boldsymbol{y}_i \\,|\\, \\boldsymbol{f}_i, \\boldsymbol{\\Sigma}_i)\n",
    "$$\n",
    "where now $\\boldsymbol{\\Sigma}_i$ is the covariance matrix of the data point $\\boldsymbol{y}_i$, and $\\boldsymbol{f}_i$ is a vector of model-predicted \"true\" values for the data. We can expand the log-Normal expression in the same way as we did for the one-dimensional case above, but now it involves some matrix algebra:\n",
    "$$\n",
    "\\ln L(\\boldsymbol{y}_i \\,;\\, \\boldsymbol{\\theta}) = \n",
    "    -\\frac{1}{2} \\left[ (\\boldsymbol{y}_i - \\boldsymbol{f}_i)^\\top \\boldsymbol{\\Sigma}_i^{-1} (\\boldsymbol{y}_i - \\boldsymbol{f}_i) + \\ln\\det\\boldsymbol{\\Sigma}_i + N \\ln(2\\pi) \\right]\n",
    "$$\n",
    "where $\\det$ is the determinant of a matrix and $N$ is the number of data points. We again want derivatives of this expression with respect to model parameters $\\boldsymbol{\\theta}$, so once again we can ignore the second and third terms because they are independent of the model parameters. The first term is a quadratic form in the data and model predictions, so the only terms that survive the derivative are those that depend on the model predictions:\n",
    "$$\n",
    "\\frac{\\partial^2}{\\partial \\theta_j \\partial \\theta_k} \n",
    "    \\ln L(\\boldsymbol{y}_i \\,;\\, \\boldsymbol{\\theta}) = \n",
    "        \\frac{1}{2} \\left[ \\frac{\\partial \\boldsymbol{f}_i}{\\partial \\theta_j} \\boldsymbol{\\Sigma}_i^{-1} \\frac{\\partial \\boldsymbol{f}_i}{\\partial \\theta_k} \\right]\n",
    "$$\n",
    "This is the quantity you might see in, for example, [Bonaca & Hogg 2018](https://ui.adsabs.harvard.edu/abs/2018ApJ...867..101B/abstract) or the draft by Sophia Lilleengen of our upcoming paper on stream information theory.\n",
    "\n",
    "Let's now do an example of this for a two-dimensional model with some nonlinear model parameters. For example, a sinusoidal model of two variables with unknown amplitudes and frequency:\n",
    "$$\n",
    "\\boldsymbol{f}(x \\,;\\, \\boldsymbol{\\theta}) = \n",
    "    \\begin{pmatrix}\n",
    "    a \\, \\sin(2\\pi \\, \\nu \\, x) \\\\\n",
    "    b \\, \\cos(2\\pi \\, \\nu \\, x)\n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "We again generate some toy data to play with (and yes I'm being naughty by using the same variable names as above):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((16, 2, 2), (16, 2))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ndata = 16\n",
    "rng = np.random.default_rng(8675309)\n",
    "true_pars = np.abs(rng.normal(0, 2.0, size=3))\n",
    "\n",
    "x = np.sort(rng.uniform(0, 2 / true_pars[2], ndata))\n",
    "yerr = 10 ** rng.uniform(-2, -0.5, (ndata, 2)) * true_pars[:2]\n",
    "\n",
    "y = np.stack(\n",
    "    (\n",
    "        rng.normal(true_pars[0] * np.sin(2 * np.pi * true_pars[2] * x), yerr[:, 0]),\n",
    "        rng.normal(true_pars[1] * np.cos(2 * np.pi * true_pars[2] * x), yerr[:, 1]),\n",
    "    ),\n",
    "    axis=-1,\n",
    ")\n",
    "\n",
    "# The inverse covariance matrix of the data (which here is diagonal):\n",
    "Cinv = np.stack([np.diag(1 / yerr[i] ** 2) for i in range(ndata)], axis=0)\n",
    "Cinv.shape, y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABbAAAALgCAYAAABF8VHCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABYlAAAWJQFJUiTwAAC8NElEQVR4nOzde3RTZdr//08PlNJIDyBILRQcLAIiUCgOIIg4FOgoKFqHL4L+POAwKvWMEoMUGEIBlWeUmWFGH9ajoiNirTiKMCA4CoqckUMF63AQaAGBNoVwaEnz+4PVTErSkJY0SdP3ay2WO/s+7GuTtu5e3LnuMLvdbhcAAAAAAAAAAEEmPNABAAAAAAAAAADgDglsAAAAAAAAAEBQIoENAAAAAAAAAAhKJLABAAAAAAAAAEGJBDYAAAAAAAAAICiRwAYAAAAAAAAABCUS2AAAAAAAAACAoEQCGwAAAAAAAAAQlEhgAwAAAAAAAACCEglsAAAAAAAAAEBQIoENAAAAAAAAAAhKJLABAAAAAAAAAEEpMtABAO60a9dOpaWluuaaawIdCgAAgPbu3avY2Fjt27cv0KEghPEMDAAAgkmwPAOTwEZQKi0t1dmzZ2s1tqCgQJKUkpLiy5DgY7xP9QPvU/3A+1R/8F7VD+7ep9o+lwA1UdtnYH621B+8V/UD71P9wPtUf/Be1Q/B/AxMAhtBqXLVyaZNm2o8Nj09XZK0YsUKn8YE3+J9qh94n+oH3qf6g/eqfnD3PvXs2TNQ4aABqe0zMD9b6g/eq/qB96l+4H2qP3iv6odgfgamBjYAAAAAAAAAICiRwAYAAAAAAAAABCUS2AAAAAAAAACAoEQCGwAAAAAAAAAQlEhgAwAAAAAAAACCEglsAAAAAAAAAEBQigx0AA3Nzp07tXDhQhUXFysiIkKSNHDgQA0bNszx2hc+/PBD7d27V88//7zXY86cOaPc3FytW7dOkZGRKisrU3JyskaNGqW2bdt6Nce6deuUl5en06dPKywsTJGRkbr99tt166231vZWAAAAAAAAADRQJLD9aPHixZo/f74mT56sXr16SZKOHz+uqVOn6ptvvpHZbFZUVNRlX8dqteqf//ynbrjhBq/HlJSU6LnnnlPr1q01a9YsGQwGSdKyZcs0fvx4GY1G9e3b1+Mcb775ppYvXy6z2awOHTpIkg4ePCiTyaT169frhRdeUFhYWO1vDAAAAAAAAECDQgkRP1mzZo3mzZun0aNHO5LXktS8eXOZTCZt27ZNM2fOlN1uv6zr7N27Vy+++KIOHz7s9ZiysjJlZ2eruLhYEydOdCSvJWno0KG65ZZbNGPGDO3atavaOfLy8rRo0SI9/vjjjuS1JLVu3VoTJkzQypUr9cYbb9TupgAAAAAAAAA0SCSw/cBqtWrOnDmy2+0aOnSoS/tVV12l1NRUrV69WitXrqzx/KtWrdKzzz6rcePG6aWXXlJ+fn6Nxn/wwQfKz8/XwIEDFR0d7dKekZGhc+fOafbs2bLZbC7tRUVFeuONN3TFFVeof//+Lu1dunRR69atlZubqx07dtQoNgAAAAAAAAANFwlsP8jLy9PJkyeVnJys+Ph4t326du0qSXr33XfdJok9ufXWW/Xqq6/q73//u959990ajbVarcrLy5OkakuOXHfddYqOjtaBAwe0atUql/bKmDt37lxtHe/K+3v77bdrFB8AAAAAAACAhosEth8sWbJEkjxuhHjNNddIkg4dOqStW7f6IyxJ0ldffaVTp05Jktq1a+e2T0REhNq0aSNJ+uyzz6q0nTlzRl9++aXH8dJ/72/r1q06ePDgZUYNAAAAAAAAoCFgE8c6tnfvXh0/flyS1KpVq2r7JSYmOo43btyonj171nlskrRhwwbHsaf4WrVqpYKCAu3atUunTp3SFVdcIelCQrq8vPyS4y++v9atW18ytoKCAqWnp1+yHwAAgDdq+lzh3L+goEApKSm+DgkAAADAJbACu45t377dcVxd+ZCL2/xZJ7ryWk2aNFHjxo2r7VcZX0VFRZUa286xBuP9AQAAAHWlvLxcxcXFHvscPXrUT9EAAACEJlZg17HK1dfShSRxdWJiYtyOqUvnz5+XxWKR5Dk2qfr4nI/dbQBZyXl+b+8vJSVFK1as8Kov6h/e2/qB96l+4H2qP3ivAuty/v799ek4+MfcuXMVFxen+++//7Lm+fDDD7Vo0SINHjxYaWlpSkpKUnh4uIqLi1VQUKCvvvpKnTt31tixY30UuXv8bKk/eK/qB96n+oH3qf7gvaofgvl9IoFdxyoTxJIUGVn9X7dzW0lJSV2G5GCxWGS3212u745zu/M9OcfaqFGjasc7tzmPBwAAAPzp6NGjWrVqlUaMGHHZc50/f15Wq1Uff/yxPv744yptjRo10n333adRo0Zd9nUAAAAaMhLYdcy5LEdFRUW1/ZzbPK1k9iXn63iK7eJ253vydg6bzeZ2PAAAAOAPdrtdO3fu1Jw5cxybmPtCly5dVF5ern379qmiokItW7ZUz549dccddyg5Odln1wEAAGioSGDXsYSEBMdx5WaH7pSVlTmOPdWS9iWDwaBGjRqpvLzcY2xS1fic78k5Vk9zOLc5jwcAAADqUm5urtauXSuLxSKr1apjx475dP7U1NTLLkUCAACA6pHArmOtW7d2HJeWllbbz7mtTZs2dRqTszZt2mjPnj06efKk7Ha7wsLC3PZzjs/5npxj9fb+nMcDAAAAdSkzM1OZmZmSpMOHD+u+++4LcEQAAACoifBABxDqUlNTFRERIcnz5oXObf7cJKjyWhUVFR53UK+MLyEhQe3bt3cZ79zH0/iLxwAAAAAAAABAdViBXccMBoO6deumzZs369ChQ9X2O3DggCQpPDxcvXv39ld46tu3rz788ENJ0qFDh9SsWTO3/Q4ePChJ6tOnT5VV2u3atVNSUpIOHTrk1f1V/n0AAAAAoeKHH37Q0qVLVVJSoiNHjqi0tFSdOnXSsGHDlJqaGujwAAAA6jUS2H5w7733avPmzSooKFBFRYXCw10Xvu/atUuSlJ6erpYtW/otti5duqhr167atm2bdu3apRtuuMGlz4kTJ3T06FFFRERo5MiRLu2jRo3SK6+84rgHdyrb7r77br9tUgkAAADUtX//+986deqUxo8fr6ioKEnSnj17lJ2dreeff14jRozQY4895vV8BQUFSk9Pr3U8K1asqPVYAAAQei7nuaKgoEApKSk+jKZ2KCHiB926dVPv3r1ltVq1fft2l3abzab169fLYDBozJgxVdqKior0yCOPaPjw4XrvvffqJL6xY8cqMjJS69atc9v+3XffSZKGDx+uq6++2qV90KBBuvbaa7V//34VFRW5tFutVu3YsUMtWrTQXXfd5dvgAQAAgAAyGAx69NFHHclrSfrVr36lF154QZL08ccfa/HixQGKDgAAoP5jBbafTJw4UU899ZQWLFigLl26OOpiS9Lnn3+u4uJimc1mtWrVqsq4NWvWaN++fZKkRYsWafTo0R6vc+7cOcdxeXm5V7F16tRJTzzxhObMmaPNmzerR48ejrbTp09r0aJF6tGjh8aNG+d2fEREhKZNm6asrCy9/fbbmjhxYpX2999/X+Hh4TKbzTIYDF7FBAAAAAS7Ro0auZTYq9SlSxe1adNGBw4c0IIFC5SRkaHGjRtfcs6UlBRWUQMAAJ+5nOeKYNnHjgS2nxgMBs2aNUuzZ8/WhAkTdNtttykmJkZbtmzR2rVrNW3atCqJ40r9+vXT8uXLdeTIEd1zzz1u5161apWWLFkiSSotLXWc/+abb/Tkk08qMvLC25ydna3Y2Fi3c2RkZCgsLEyvvPKKhgwZoo4dO6q4uFh5eXnq3LmzsrKyqiTdL9aiRQu98sorysnJ0eTJk/Wb3/xGkZGR+vbbb/XDDz9o9uzZuuaaa7z++wIAAACC3YgRIzw+I3fs2FEHDhxQaWmptm3bpl69evkxOgAAgNBAAtuPmjVrppkzZ+rHH3/Ujh079Msvv6h79+4aO3ZslY8cOktMTNSbb77pcd5bb71Vt95662XHN3ToUPXv31/r16/Xnj17ZDAYNHnyZLVu3dqr8a1bt9af//xnbd++Xbt371ZFRYUGDhyoZ555xuODPQAAAFAfXWpvl4SEBMfxgQMHSGADAADUAgnsAOjQoYM6dOgQ6DDcMhgMGjhwYK3Hh4WFqWvXruratasPo4I/TZkyxe0xAAAAqjp79qy+//57JSYmKjk52aXdefP2srIyf4YGAABCXEPK35DABlDF1KlTHceh/gMQAACgts6dO6fx48dr//79ioiI0N///ne1bdu2Sp+SkhLHccuWLf0cIQAACGUNKX8TfukuAAAAAABnX3/9tfbv3y9Jstls+vnnn136nDhxwnGckpLit9gAAABCCQlsAAAAAKgh5+R0r169dOONN1ZpLy8vV35+viTphhtuUJs2bfwaHwAAQKgggQ0AAAAAFykqKtIjjzyi4cOH67333nNprywXkpaWJrPZrMaNG1dp37x5s06dOqVGjRrpD3/4g19iBgAACEXUwAYAAADQIJw7d85xfKlNFdesWaN9+/ZJkhYtWqTRo0dXae/Vq5dSUlJUVlamsLCwKm1lZWV65513FBkZqaeffjpoN3AH0DA1pI3fAISGMLvdbg90EMDFevbsKUnatGlTgCNpeJx/AePHAwAAF/BsUn/l5uZq7dq1kqRjx46psLBQkhQdHa2UlBSFhYUpMjJSs2bNqjKuqKhIkydP1pEjR/S73/1OY8aMcZm7uLhYL7/8sqKiotS/f3/FxcXp2LFjWrx4scLDw/X444/r+uuv9zpWvs4AXA5vE9P8zgeEBn98LwfLswkJbASlYPkGaYh4mAEAwBXPJvDk559/1o4dO2SxWJSQkKBf/epXtVp1zdcZgMvh7e9y/M4HhIaGlMCmhAgAAAAAXIbk5GQlJycHOgwADZjFYnF5HRcXF6BoAMC32MQRAAAAABoAi8Wi4cOH69Zbb9Xw4cNdEl4A6iez2aykpKQq55KSkmQ2m136ukt0A0CwI4ENAAAAACGuMsH16aef6ssvv9Snn35abYILQP1hNps1adIkWa3WKuetVqsmTZpU5Xu8JoluAAgmJLABAAAAIITVJMEFoH6wWCyaO3eupk2b5rFfTk6OSktL+TkAoF4jgQ0AAAAAIcpisSgnJ8djn8oEF4D6oXIl9RNPPKGysjKPfa1Wq9555x1+DgCo10hgAwAAAECIys3NdVlxeTGr1arc3Fw/RQTgclS3ktqTlStX8nMAQL1GAhsAAAAAQlRRUZFP+wEIHG8+UXE5+DkAIFiRwAYAAACAEJWYmOjTfgACx5tPVFzMYDDoN7/5jVd9+TkAIFiRwAYAAACAEJWZmSmDweCxj8FgUGZmpp8iAlBbtVkhbTQadd999/FzAEC9RgIbAAAAAEJUXFycjEajxz5Go1GxsbF+ighAbdVkhbTBYND06dNlMpn4OQCg3osMdAAAAAAAgLpjMpkkSTk5OVXKDxgMBhmNRkc7gOCWmZmpJ5980qsyIoWFhVUS0vwcAFCfsQIbAAAAAEKcyWRSYWFhlXOFhYUkrYB6xJuV1JXcrabm5wCA+ooV2AAAAADQAMTGxio7O7vKawD1i6eV1N6szL74+56fAwDqAxLYAAAAANBATJkyJdAhALhMJpNJWVlZys3NVVFRkRITE5WZmak5c+YEOjQAfmKxWFxex8XFBSiaukcCGwAAAAAAoB6JjY3VQw89VOWct/9A5fxJDAD1j9lsVk5OTpVzSUlJIV3PngQ2AAAAAABAA8EnMYD6y2w2a9KkSS7nrVar43woJrHZxBEAAAAAAAAAgpjFYnFZeX2xnJwclZaW+iki/yGBDcDBXQ0lAAAAAAAABFZubu4lN2u1Wq3Kzc31U0T+QwIbgKQLH0NJSkqqci4pKUlmszlAEQEAAAAAAECSioqKfNqvPqEGNoAGW0MJAAAAAACgPkhMTPRpv/qEFdhAA9eQaygBAAAAAADUB5mZmTIYDB77GAwGZWZm+iki/yGBDTRwDbmGEgAAAAAAQH0QFxcno9HosY/RaFRsbKyfIvIfSogADVxDrqEEAAAAAABQX1SWd83JyamyGNFgMMhoNIZs+VcS2EAD15BrKAEAAAAAANQnJpNJWVlZiouLc5wrLCwMyZXXlSghAjRwDbmGEgAAAAAAQH1zcbI6lJPXEglsoMFryDWUAAAAAAAAENwoIQKgwdZQAgAAAAAAQHAjgQ1AUsOsoQQAAAAAAIDgRgkRAA4NrYYSAAAAAAAAghsJbAAAAAAAAABAUCKBDQAAAAAAAAAISiSwAQAAAAAAAABBiQQ2AAAAAAAAACAokcAGAAAAAAAAAAQlEtgAAAAAAAAAgKAUGegAAABAYEyZMsXtMQAAAADAv/j9rHphdrvdHugggIv17NlTkrRp06YAR9LwhIWFOY758QCENr7fAe/xbAJ/4OsMAICGq6a/n/nj97lgeTahhAgAAAAAAAAAICiRwAYAAAAAAAAABCUS2AAAAAAAAACAoEQCGwAAAAAAAAAQlEhgAwAAAAAAAACCEglsAAAAAAAAAEBQigx0AAAAAAAAAAAA72VnZwc6BL8hgQ0AAAAAAAAA9ciUKVMCHYLfUEIEAAAAAAAAABCUSGADAAAAAAAAAIISCWwAAAAAAAAAQFAigQ0AAAAAAAAAAWKxWDy+buhIYAMAAAAAAABAAJjNZiUlJVU5l5SUJLPZHKCIgk9koAMAAAAAAH+aO3eu4uLidP/99/tkvn379ukf//iHDh48qKioKNlsNvXu3VuZmZlq3LixT64BAABCj9ls1qRJk1zOW61Wx3mTyeTvsIIOCWwAVWRnZwc6BAAAgDpz9OhRrVq1SiNGjPDJfKtXr9bs2bOVlZWlF198UZJ06tQpvfzyy8rKytLLL7+suLg4n1wLAACEDovFopycHI99cnJylJWVpdjYWD9FFZwoIQKgiilTpjj+AAAAhAq73a4dO3Zo4sSJOnXqlE/mzM/P18yZMzVo0CANHjzYcf6KK67QxIkT9csvvyg7O1tlZWU+uR4AAAgdubm5slqtHvtYrVbl5ub6KaLgxQps+ERJSYliYmIUFRXltt1ms6mkpETNmzf3c2QAAABoyHJzc7V27VpZLBZZrVYdO3bMJ/PabDbNnj1bZWVlysjIcGlv0qSJBgwYoCVLlujDDz/U6NGjfXJdAAAQGoqKinzaL5SRwPaznTt3auHChSouLlZERIQkaeDAgRo2bJjjdW3Vtvbes88+q23btnl1jddff12dOnVyOZ+Tk6OioiINGTJE3bp1U/PmzWWz2XT06FHl5+fryy+/1NixY9WnT59a3x8AAABQU5mZmcrMzJQkHT58WPfdd59P5l25cqUOHTqkmJgYtW/f3m2frl27asmSJfroo4905513ymAw+OTaAACg/ktMTPRpv1BGAtuPFi9erPnz52vy5Mnq1auXJOn48eOaOnWqvvnmG5nN5mpXMF+Kv2rvxcfHuz1//vx5FRUV6a233nI75rnnntOvf/3ry74+AAAAEAyWLFkiSWrTpk21C1GuueYaSdLJkyf11Vdf6be//a3f4gMAAMEtMzNTTz75pMcyIgaDwfEP8Q0ZCWw/WbNmjebNm6cHH3zQkbyWpObNm8tkMun+++/XzJkz9dJLLyksLKxGc1fW3hs8eLDb2nv33nuvsrOzNXv27GoT5J06ddIDDzygK6+80m17QUGB/vSnP+mqq66qNo6bbrpJe/fu1dGjRxUVFaXExET1799fw4YNa/DF5gEAABA6Tp06pV27dkmSWrVqVW0/57aNGzeSwAYAAA5xcXEyGo2aNGlStX2MRiM5NZHA9gur1ao5c+bIbrdr6NChLu1XXXWVUlNTtXr1aq1cuVKDBg3yem5f1N4LCwvThAkT1KZNm2qvc+7cOXXs2FHh4dXv+zlixAh169bN69gBAACA+mjnzp2qqKiQVP0nFKULz+LR0dE6e/astm/f7tXcBQUFSk9Pr3VsK1asqPVYAADgXyaTSdKF0rzOK7ENBoOMRqOj/XJcznNFQUGBUlJSLjuGy1V9NhI+k5eXp5MnTyo5ObnaB9yuXbtKkt59913ZbDav5/a29p4kffTRR24/ltC8eXOPyWtJWrdunfr16+d1XLg8U6ZMcfwBAABAcDl+/LjjuEmTJh77VrZbLJYaPecDAICGwWQyqbCwsMq5wsJCnySvQwUrsP2gsj5e27Ztq+1TWR/v0KFD2rp1q3r27FmjuS+n9l7fvn09XqO0tFTLly/XX/7yF69iwuWbOnWq45gkNgAAQHCxWCyO48hIz79SVbbb7XZZLBY1a9bMY/+UlBRWUQMA0MBcXCbEl2VDLue5wtv8ZF0jgV3H9u7d61ih4ak+nvOOohs3bvTqC8RXtfcGDBjg8TrvvPOO7rjjDjVt2vSSMa1atUrr16+XxWJRUVGRJCk1NVV33323WrdufcnxAAAAQLBr3Lix47iylEh1nNudxwEAAMA7JLDrmHOtO0/18ZzbduzY4dXcdVl7r9LmzZu1e/duPfroo5fsO2/ePN15552aOHGi49zXX3+tmTNn6l//+pcmTJiggQMHen1t6v8BAABfCoX6fwgOCQkJjuPy8nKPfcvKyiRJUVFRMhgMdRoXAABAKKIGdh3ztj5eTEyM2zG+mNu5vSa1906ePKnXXntNzz33XLXlSZz17NnTZZPKm2++WaNHj1Z5eblefvll/fjjj15dGwAAAAhWzvvHlJaWVtuvoqLCsQcNn0YEAACoHVZg1zFv6+M5t5WUlPh0bud2b2vvSdLrr7+uQYMGeazdXalRo0bV1tLOyMjQ22+/rfLycr311luaMWPGJeeTqP8HAAB8KxTq/yE4tG/fXvHx8SopKfG4+OTEiROOT0ympaX5KzwAAICQwgrsOuZtfTzntujoaJ/OfXG7N7X3vv32W61bt06ZmZlexWIymXT99de7bWvWrJlatGgh6UJJkjNnzng1JwAAABCMwsLC1KdPH0kXNmGvzsGDBx3Hl9o4HQAAAO6RwK5j3tbHq6yNJ3muZ12buZ3n96b23tmzZ/XnP/9Z3bp1u2RpkkqX2uCxcsW3zWZTYWGhV3MCAAAAwWrkyJGKiIjQ4cOHq3wy0lnlhuvdu3evdrEHAAAAPCOBXceca915qo/n3OZcU8+Tuqq99/HHH+uXX35xrJr2xvHjx/Xtt9+quLjYbXtYWJjj2DlZDwAAAASjoqIiPfLIIxo+fLjee+89l/akpCQNGzZMdrtd69atczvH2rVrFRkZqYcffriuwwUAAEFuypQpjj+oGWpg17HU1FRFRETIZrN5rI/n3OZtjcW6qL1XVlamDz/8UNKFetne2LNnj55++mmdPn1aV111ld5++22XTR+d63rXJDEOAKgbF68WtFgsiouLC1A0AOAf586dcxxfalHFmjVrtG/fPknSokWLNHr0aJc+48aN0759+7Rw4UL179+/yqcXN2zYoPz8fE2YMEEdO3b0zQ0AAIB6a+rUqY5jktg1QwK7jhkMBnXr1k2bN2/2WB/vwIEDkqTw8HD17t3bq7kra+8tXbrUZ7X3vvvuO508eVLShcS3N/Ly8nT69GlJ0pEjR2S1WhUbG1ulT+XK7ObNm6t58+ZezQsAqBtms1k5OTlVziUlJcloNMpkMgUoKgCoG7m5uVq7dq0k6dixY47zn3zyifLz8xUWFqbIyEjNmjWryrh+/fpp+fLlOnLkiO655x63c0dGRmrq1Kl67bXX9NRTT2nEiBFKSEjQDz/8oOXLl+v5559Xenp63d0cAABAA0AC2w/uvfdebd68WQUFBaqoqFB4uGvllsr6eOnp6WrZsqXXc48cOVLLly931N5zt3quJrX3tmzZ4jjev3+/VzFUJrrDwsI0ZswYl+T1rl27dPbsWUlSRkZGlXIiAAD/MpvNmjRpkst5q9XqOE8SG0AoyczM9HpjcmeJiYl68803L9kvJiZGRqNRP//8s7Zu3ao9e/aoXbt2evPNNy+59wwAAAAujRrYftCtWzf17t1bVqtV27dvd2m32Wxav369DAaDxowZU6XN37X3jh49WuXa3tSrbtu2rSTpwQcf1P333+/Svnr1akkXfgmozS8PAADfsFgsLiuvL5aTk+NxXwUAgHvJyckaPny4Ro0apVtuuYXkNQAAgI+QwPaTiRMnql27dlqwYIFsNluVts8//1zFxcWaPHmyWrVqVaWtsvbemTNntGjRIrdzjxs3Tt27d9fChQt15syZKm2Vtfeefvppr2rvxcTEOI4rKiq8SmLccccdMhgMbpPdRUVF+uyzz5SQkKDJkyfzIO8Fd3VpAcAXcnNzHZv6VsdqtSo3N9dPEQEAAAAAJCk7O9vxB1VRQsRPDAaDZs2apdmzZ2vChAm67bbbFBMToy1btmjt2rWaNm2aevTo4TLO37X3Bg0apK+++kp2u10xMTEu5UDcadWqlWbNmqVXX31Vp0+fVpcuXRQdHa19+/YpNzdXqampGj9+vK688kqvYmjIqEsLoC4VFRX5tB8AAAAAwDfY2LF6YXa73R7oIBqaH3/8UTt27FBZWZmSk5OVlpamqKgon8xdWXvParUqMTFRvXr1qvGq5x07dig/P1+pqalKSUnxelxFRYXy8/P1008/6ezZs2rZsqU6d+7ssqrcGz179pQkbdq0qcZj66vq6tJWmj59OklsAJdl/vz5Gjt2rFf9HnroIT9EBNQfDfHZBP7H1xkAAKHJYrEoPj7e8bqkpMTtPnbBJlieTUhgIygFyzeIv1gsFiUlJXn8aL/BYFBhYaFXq+IBwB1+1gC119CeTRAYfJ0BABB6Kj9t7/x7mMFgqBeftg+WZxNqYANBgLq0APwhLi5ORqPRYx+j0UjyGgAAAAB8oPLT9hfnfKxWqyZNmiSz2RygyOoXEthAEKAuLQB/MZlMmj59ukt5KYPBQKkiAAAAAPARi8Xiss/ZxXJyclRaWuqniOovEthAEEhMTPRpPwDwxGQyqbCwsMq5wsJCktcAAAAA4CN82t53SGADQSAzM/OSm20aDAZlZmb6KSIAoe7iMiGUDQEAAAAA3+HT9r5DAhsIAtSlBQAAAAAACB182t53SGADQYK6tAAAAAAAAKGBT9v7DglsIIhQlxYAAAAAAKD+49P2vhMZ6AAAVEVdWgAAAAAAgPqvckFiTk5OlQ0dDQaDjEYjCxa9RAIbAAAAAAAAAOqAyWRSVlaW4uLiHOcKCwtZsFgDlBABAAAAAAAAgDrCp+0vDyuwAQAAAAAAAOAyTJkyxe0xLh8JbAAAAAAAAAC4DFOnTnUck8D2LUqIAAAAAAAAAACCEglsAAAAAAAAAEBQooQIAAAAAAAAANSh7OzsQIdQb5HABgAAAAAAAIA6RF3s2qOECAAAAAAAAAAgKJHABgAAAAAAAAAEJRLYAAAAAAAAAICgRAIbAAAAAAAAABCU2MQRCELsTAsAAAAAAACQwAaCEjvTAgAAAAAA1A8Wi8XldVxcXICiCT2UEAEAAAAAAACAWjCbzUpKSqpyLikpSWazOUARhR5WYAMAAAAAAABADZnNZk2aNMnlvNVqdZw3mUz+DivksAIbAAAAAAAAAGrAYrEoJyfHY5+cnByVlpb6KaLQRQIbAAAAAAAAAGogNzdXVqvVYx+r1arc3Fw/RRS6SGADAAAAAAAAQA0UFRX5tB+qRwIbAAAAAAAAAGogMTHRp/1QPRLYAAAAAAAAAFADmZmZMhgMHvsYDAZlZmb6KaLQRQIbAAAAAAAAAGogLi5ORqPRYx+j0ajY2Fg/RRS6IgMdAAAAAAAAAADUNyaTSZKUk5NTZUNHg8Ego9HoaMflIYENAAAAAAAAALVgMpmUlZWluLg4x7nCwkJWXvsQJUQAAAAAAAAAoJYuTlaTvPYtEtgAAAAAAAAAgKBEAhsAAAAAAAAAEJSogQ0AQAOVnZ0d6BAAAAAAAPCIBDYAAA3UlClTAh0CAAAAAAAeUUIEAAAAAAAAABCUSGADAAAAAAAAAIISCWwAAAAAAAAAQFAigQ0AAAAAAAAACEoksAEAAAAAAAAAQYkENgAAAAAAAAAgKEUGOgAAAAAAAAAAqM+ys7MDHULIIoENAAAAAAAAAJdhypQpgQ4hZFFCBAAAAABqoby8XMXFxR77HD161E/RAAAAhCZWYAMAAAAIeTt37tTChQtVXFysiIgISdLAgQM1bNgwx+ua+vDDD7Vo0SINHjxYaWlpSkpKUnh4uIqLi1VQUKCvvvpKnTt31tixY315KwAAAA0KCWwAAAAAIW3x4sWaP3++Jk+erF69ekmSjh8/rqlTp+qbb76R2WxWVFRUjec9f/68rFarPv74Y3388cdV2ho1aqT77rtPo0aN8sk9AAAANFSUEAEAAAAQstasWaN58+Zp9OjRjuS1JDVv3lwmk0nbtm3TzJkzZbfbazV/ly5ddN1116lx48Zq1KiRkpKSNHz4cP3tb38jeQ0AAOADrMAGAAAAEJKsVqvmzJkju92uoUOHurRfddVVSk1N1erVq7Vy5UoNGjSoxtdITU3V/fff74twAQAA4AYJbDRYzrvDslMsAABA6MnLy9PJkyfVtm1bxcfHu+3TtWtXbdq0Se+++64GDhxY63rYAAAAqBsksNFgTZ061XFMAhsAACD0LFmyRJLUtm3bavtcc801kqRDhw5p69at6tmzp19iAwAAgHdIYAMAAAAIOXv37tXx48clSa1ataq2X2JiouN448aNtUpg//DDD1q6dKlKSkp05MgRlZaWqlOnTho2bJhSU1NrNFdBQYHS09NrHEOlFStW1HosAAAIPZfzXFFQUKCUlBQfRlM7JLABAAAAhJzt27c7jqsrH3Jx244dO2p8nX//+986deqUxo8fr6ioKEnSnj17lJ2dreeff14jRozQY489VuN5AQAAcAEJbAAAAAAhp3L1tSQ1adKk2n4xMTFux3jLYDDo0UcfVVhYmOPcr371K73wwgt6+umn9fHHH+vqq6/WnXfe6dV8KSkprKIGAAA+cznPFcFSWi080AEAAAAAgK9ZLBbHcWRk9et2nNtKSkpqdI1GjRqpT58+VZLXlbp06aI2bdpIkhYsWKBz587VaG4AAABcQAIbAAAAQMhp3Lix47iioqLafs5t0dHRNbrGiBEjdM8991Tb3rFjR0lSaWmptm3bVqO5AQAAcAEJbAAAAAAhJyEhwXFcXl5ebb+ysjLHsada2e5ER0erUaNGXsVw4MCBGs0NAACAC0hgAwAAAAg5rVu3dhyXlpZW28+5rbLkh7fOnj2rdevW6eeff3bbHh7+31+3nBPlAAAA8B6bOAIAAAAIOampqYqIiJDNZvO4OaNzW002Kjp37pzGjx+v/fv3KyIiQn//+9/Vtm3bKn2ca2q3bNnS++ABAADgwApsAAAAACHHYDCoW7dukqRDhw5V26+ytEd4eLh69+7t9fxff/219u/fL0my2WxuV2GfOHHCcZySkuL13AAAAPgvEtiQJB09etRje3Fxsc6fP++naAAAAIDLd++990qSCgoKqt3IcdeuXZKk9PT0Gq2Sdk5O9+rVSzfeeGOV9vLycuXn50uSbrjhhhqXJwEAAMAFlBDxs507d2rhwoUqLi5WRESEJGngwIEaNmyY43Vt7du3T//4xz908OBBRUVFyWazqXfv3srMzKyyC/vFjh07pjFjxigtLU233HKLUlJSFBMTo9OnT+vAgQPasGGD/vOf/2jmzJmKjY31GMOKFSu0dOlS2Ww2VVRUKCYmRpmZmerVq9dl3RsAAABQU926dVPv3r313Xffafv27Y4V2ZVsNpvWr18vg8GgMWPGVGkrKirS5MmTdeTIEY0cOVKjR4+u0l5ZLiQtLU1ms1lhYWFV2jdv3qxTp06pUaNG+sMf/lAHdwcAANAwkMD2o8WLF2v+/PmaPHmyI6F7/PhxTZ06Vd98843MZrOioqJqNffq1as1e/ZsZWVl6cUXX5QknTp1Si+//LKysrL08ssvKy4uzu3Y8+fPy263a8OGDdqwYYNLe48ePTR9+nSPyeuKigrl5OTop59+0vTp05WUlCRJ2r17t0wmkzIyMvTwww/X6t4AAA3TlClT3B4DQE1MnDhRTz31lBYsWKAuXbpUWTTy+eefq7i4WGazWa1ataoybs2aNdq3b58kadGiRS4J7F69eiklJUVlZWUuyeuysjK98847ioyM1NNPP60OHTrUzc0BAAA0AGF2u90e6CAagjVr1uiPf/yjHnzwQf2///f/qrQdOXJE999/v2666Sa99NJLLg/Al5Kfn68JEyZo8ODBevLJJ6u0nTlzRvfee6/atm2r2bNnu02QHz58WI899phSU1O1ZcsWWa1WxcfHq3379srIyFC/fv0uGdO8efOUl5en119/XZ06darS9sUXX2jWrFl6/PHHdeedd3p1T5Ub6GzatMmr/rXhfE98GwBA8OHnNIKJP55NUHdOnDih2bNnq6ysTLfddptiYmK0ZcsWrV27Vk888YTbTws6r8D+3e9+57JCW7pQZu/ll19WVFSU+vfvr7i4OB07dkyLFy9WeHi4Hn/8cV1//fVex8nXGQAACCbB8mzCCmw/sFqtmjNnjux2u4YOHerSftVVVyk1NVWrV6/WypUrNWjQIK/nttlsjofxjIwMl/YmTZpowIABWrJkiT788EOXlSOVYmJi9NJLL3l/U062bdumvLw8tWvXziV5LUkDBgzQ3Llz9fe//12//vWvlZiYWKvrAAAAALXRrFkzzZw5Uz/++KN27NihX375Rd27d9fYsWOr/QRkYmKi3nzzTY/zJiQkaMaMGfr555+1Y8cOFRQUKCEhQc888wyrrmuBT94AAAB3SGD7QV5enk6ePKm2bdsqPj7ebZ+uXbtq06ZNevfddzVw4ECv62GvXLlShw4dUkxMjNq3b1/t3EuWLNFHH32kO++8UwaDoba34tY777wj6cLmNO40atRIHTt21ObNm/Xee+/pueee8+n1AQAAAG906NChThLLycnJSk5O9vm8Dc3UqVMdxySwAQBApfBAB9AQLFmyRNJ/N3px55prrpEkHTp0SFu3bq3x3G3atKk26V0598mTJ/XVV195Pbc3Dh48qO+//16S1K5du2r7VcawatUqnTlzxqcxAAAAAAAAAAhNrMCuY3v37tXx48clyWVjGGfOZTU2btzoqDHjyalTp7Rr165Lzu3ctnHjRv32t79126+0tFR5eXk6ePCgjh07piNHjqhly5a6+eabdfvtt6tx48YuYzZu3Oj2OhervL/y8nJ9//336t27t+ebAwAAAAAAANDgkcCuY9u3b3ccV1c+5OK2HTt2eDX3zp07VVFRccm5mzRpoujoaJ09e7ZKPM5KS0s1c+ZMPfbYY2rdurWkCxtAzp07V3/729+0bNkyzZgxQy1atKgyrjb3t337dq8S2AUFBUpPT79kv+qsWLGi1mMBAEDouZznioKCAqWkpPgwGgAAAADeoIRIHatcfS1dSCRXJyYmxu0YX8zt3G6xWGSz2Vzaz5w5o3HjxjmS15Vjnn32WbVp00b79u3T1KlTHQnzmsZQm/sDAAAAAAAA0LCxAruOWSwWx3FkZPV/3c5tJSUlPp3bud1ut8tisahZs2aOtkaNGikpKcltje6IiAgNHjxY8+fP1+7du7V69WoNGDCgxjE4tzmP8SQlJYVV1AAAwGcu57nCm/JuAAAAAHyPFdh1zLlu9MWrl505t0VHR/t07ovbL65l3bx5c82dO7fasR07dnQcr1mzplYxeLo+AAAAAAAAALhDAruOJSQkOI7Ly8ur7VdWVuY49lRLujZzO88fFRUlg8Hg0t60aVOvrnPgwIFaxeB8f85jAAAAAAAAAKA6JLDrmHNN6dLS0mr7Obe1adPGq7md+3mau6KiQlar1SUeZ7t379aGDRtkt9td2sLD//tl4pyIvng+X98fAAAAAAAAgIaNBHYdS01NVUREhCTPmxc6t3lbY7F9+/aO1dqe5j5x4oSjhEdaWppL+4IFCzR+/Hi9+OKLeuutt1zanWtyt2jRokqb83y+vj8AAAAAAAAADRsJ7DpmMBjUrVs3SdKhQ4eq7VdZmiM8PFy9e/f2au6wsDD16dPnknMfPHjQcdy3b98qbTabTe+//77j9Z49e1zGnzhxwnHcoUOHKm3dunVzlCTxJobqNosEAAAAAAAAgIuRwPaDe++9V5JUUFBQ7UaHu3btkiSlp6erZcuWXs89cuRIRURE6PDhw7JYLB7n7t69u66//voqbaWlpY7a1fHx8RozZozL+O+//17ShYR5RkZGlbbo6Gjdddddki6UIalOZQyjR4/25rYAAAAAAAAAgAS2P3Tr1k29e/eW1WrV9u3bXdptNpvWr18vg8HgkkAuKirSI488ouHDh+u9995zGZuUlKRhw4bJbrdr3bp1bq+/du1aRUZG6uGHH3Zpi4+PV3x8vMLDw/WXv/xF1113XZX2iooKffPNN5Kku+66S1dffbXLHHfffbdatGihbdu26cyZMy7tBw8e1IEDB5SSkqJbb73VbYwAAAAAAAAAcDES2H4yceJEtWvXTgsWLJDNZqvS9vnnn6u4uFiTJ09Wq1atqrStWbNG+/bt05kzZ7Ro0SK3c48bN07du3fXwoULXRLIGzZsUH5+vp5++ml17NjRZWxYWJhGjRqliooKl7gk6aOPPtKJEyfUu3dvjR071u31DQaD/vjHP8put2vhwoUu7W+//bZatGihadOmOeqBAwAAAAAAAMClhNntdnugg2goTpw4odmzZ6usrEy33XabYmJitGXLFq1du1ZPPPGEevXq5TKmqKhIkydP1pEjR/S73/3ObYkPSTp9+rRee+017du3TyNGjFBCQoJ++OEHLV++XA8++KDS09M9xvbBBx/oiy++UEZGhpKSkmSz2bR27VqtXbtWmZmZGjlypMLCwjzOsWvXLs2aNUvXX3+9+vTpo/Pnz+uLL75QSUmJjEaj29Xb1anc6HHTpk1ej6kp5/vh2wAAgg8/pxFM/PFsAjT0rzN+7gMAEFyC5dmEBHYA/Pjjj9qxY4fKysqUnJystLQ0RUVF+WTun3/+WVu3bpXValViYqJ69erl2GTxUqxWqzZu3KijR48qIiJCSUlJSk1NrVFsNptNmzdv1p49exQeHq5OnTqpS5cuNb4PEtgA0LBZLBbFx8c7XpeUlCguLi5wAaHBC5aHd4S2hv51xvM5AADBJVieTSIDevUGqkOHDurQoUOdzJ2cnKzk5ORajTUYDBowYMBlXT8iIkK9evVyu5ocAABvmM1m5eTkVDmXlJQko9Eok8kUoKgAAAAAAIFAAhsAAAQNs9msSZMmuZy3Wq2O8ySxAQAAAKDhYBNHAAAQFCwWi8vK64vl5OSotLTUTxEBAAAAAAKNBDYAAAgKubm5slqtHvtYrVbl5ub6KSIAAAAAQKCRwAYAAEGhqKjIp/0AAAAAAPUfCWwAABAUEhMTfdoPAAAAAFD/kcAGAABBITMzUwaDwWMfg8GgzMxMP0UEAAAAAAg0EtgAACAoxMXFyWg0euxjNBoVGxvrp4gAAAAAAIEWGegAAAAAKplMJklSTk5OlQ0dDQaDjEajox0AAAAA0DCQwAYAAEHFZDIpKytLcXFxjnOFhYWsvAaAEGaxWFxeO/9/AAB8ZcqUKW6PAQQvSogAAICgc3GymuQ1AIQus9mspKSkKueSkpJkNpsDFBGAUDZ16lTHHwD1AyuwAQAAAAABYTabNWnSJJfzVqvVcZ7yUQAANGyswAYAAAAA+J3FYlFOTo7HPjk5OSotLfVTRAAAIBiRwAYAAAAA+F1ubm6VDXvdsVqtys3N9VNEAAAgGJHABgAAAAD4XVFRkU/7AQCA0EQCGwAAAADgd4mJiT7tBwAAQhMJbAAAAACA32VmZspgMHjsYzAYlJmZ6aeIAABAMCKBDQAAAADwu7i4OBmNRo99jEajYmNj/RQRAAAIRpGBDgAAAAAA0DCZTCZJUk5OTpUNHQ0Gg4xGo6MdAAA0XCSwAQAAAAABYzKZlJWVpbi4OMe5wsJCVl4DAABJlBABAAAAAATYxclqktcAAKASCWwAAAAAAAAAQFAigQ0AAAAAAAAACEoksAEAAAAAAAAAQYkENhoki8Xi8TUAAAAAAACAwCOBjQbHbDYrKSmpyrmkpCSZzeYARQQAAAAAAADAnchABwD4k9ls1qRJk1zOW61Wx3mTyeTvsAAAAAAAAAC4wQpsNBgWi0U5OTke++Tk5Ki0tNRPEQEAAAAAAADwhBXYaDByc3NltVo99rFarcrNzdVDDz3kp6gAANXJzs4OdAgAAAAAgAAjgY0Go6ioyKf9AAB1a8qUKYEOAQAAAAAQYJQQQYORmJjo034AAAAAAAAA6hYJbDQYmZmZMhgMHvsYDAZlZmb6KSIAAAAAAAAAnpDARoMRFxcno9HosY/RaFRsbKyfIgIAAAAAAADgCTWw0aCYTCZJUk5OTpUNHQ0Gg4xGo6MdAAAAAAAAQOCRwEaDYzKZlJWVpbi4OMe5wsJCVl4DAACEsJ07d2rhwoUqLi5WRESEJGngwIEaNmyY43Vt7du3T//4xz908OBBRUVFyWazqXfv3srMzFTjxo19ET4AAECDRQIbDdLFyWqS1wAAAKFr8eLFmj9/viZPnqxevXpJko4fP66pU6fqm2++kdlsVlRUVK3mXr16tWbPnq2srCy9+OKLkqRTp07p5ZdfVlZWll5++eUqCycAAABQM9TABgAAABCy1qxZo3nz5mn06NGO5LUkNW/eXCaTSdu2bdPMmTNlt9trPHd+fr5mzpypQYMGafDgwY7zV1xxhSZOnKhffvlF2dnZKisr88m9AAAANEQksAEAAACEJKvVqjlz5shut2vo0KEu7VdddZVSU1O1evVqrVy5skZz22w2zZ49W2VlZcrIyHBpb9KkiQYMGKCdO3fqww8/rPU9AAAANHQksAEAAACEpLy8PJ08eVLJycmKj49326dr166SpHfffVc2m83ruVeuXKlDhw4pJiZG7du39zj3Rx99VGUDcQAAAHiPBDYAAACAkLRkyRJJUtu2bavtc80110iSDh06pK1bt9Z47jZt2lS7CWTl3CdPntRXX33l9dwAAAD4LxLYAAAAAELO3r17dfz4cUlSq1atqu2XmJjoON64caNXc586dUq7du265NzObd7ODQAAgKoiAx0AAAAAAPja9u3bHcfVlQ+5uG3Hjh1ezb1z505VVFRccu4mTZooOjpaZ8+erRKPJwUFBUpPT/eqrzsrVqyo9VgAABB6Lue5oqCgQCkpKT6MpnZYgQ0AAAAg5FSuvpYuJJKrExMT43aML+Z2brdYLDWqsQ0AAIALWIENAAAAIORYLBbHcWRk9b/2OLeVlJT4dG7ndrvdLovFombNmnnsn5KSwipqAADgM5fzXNGzZ08fRlJ7rMAGAAAAEHIaN27sOK4s9+GOc1t0dLRP57643XkcAAAAvEMCGwAAAEDISUhIcByXl5dX26+srMxx7KmedW3mdp4/KipKBoPBq/kBAADwX5QQAQAAABByWrdu7TguLS2ttp9zW5s2bbya27mfp7krKipktVpd4oF72dnZgQ4BQIhzLgFV+TouLi5A0QDwFglsAAAAACEnNTVVERERstlsHjdndG7zts5j+/btFR8fr5KSEo9znzhxwlFCJC0tzcvIG64pU6YEOgQAIcxsNisnJ6fKuaSkJBmNRplMpgBFBcAblBABAAAAEHIMBoO6desmSTp06FC1/Q4cOCBJCg8PV+/evb2aOywsTH369Lnk3AcPHnQc9+3b16u5AQC+ZzabNWnSJMenYipZrVZNmjRJZrM5QJEB8AYJbAAAAAAh6d5775UkFRQUVLvZ4q5duyRJ6enpatmypddzjxw5UhERETp8+LDLR9Ivnrt79+66/vrraxI6AMBHLBaLy8rri+Xk5HgsCQUgsEhgAwAAAAhJ3bp1U+/evWW1WrV9+3aXdpvNpvXr18tgMGjMmDFV2oqKivTII49o+PDheu+991zGJiUladiwYbLb7Vq3bp3b669du1aRkZF6+OGHfXNDAIAay83NdVl5fTGr1arc3Fw/RQSgpkhgAwAAAAhZEydOVLt27bRgwQLZbLYqbZ9//rmKi4s1efJktWrVqkrbmjVrtG/fPp05c0aLFi1yO/e4cePUvXt3LVy4UGfOnKnStmHDBuXn5+vpp59Wx44dfXtTAACvFRUV+bQfAP9jE0cAAAAAIctgMGjWrFmaPXu2JkyYoNtuu00xMTHasmWL1q5dq2nTpqlHjx4u4/r166fly5fryJEjuueee9zOHRkZqalTp+q1117TU089pREjRighIUE//PCDli9frueff17p6el1fYsAAA8SExN92g+A/4XZ7XZ7oIMALla5A/ymTZvq7BphYWGOY74NAACAJ/54NkHd+/HHH7Vjxw6VlZUpOTlZaWlpioqK8sncP//8s7Zu3Sqr1arExET16tVLBoOhRnPwdQYAvmexWJSUlOSxjIjBYFBhYaFiY2P9GBkQ/ILl2YQV2AAAAAAahA4dOqhDhw51MndycrKSk5PrZG4AQO3FxcXJaDRq0qRJ1fYxGo0kr4EgRgIbAAAAAAAAIctkMkmScnJyqqzENhgMMhqNjnYAwYkENgAAAAAAAEKayWRSVlaW4uLiHOcoGwLUD+GBDgAAAAAAAACoaxcnq0leA/UDCWwAAAAAAAAAQFAigQ0AAAAAAAAACEoksAEAAAAAAAAAQYkENgAAAAAAAAAgKJHAhk+UlJSorKys2nabzabjx4/7MSIAAAAAAAAA9V1koANoaHbu3KmFCxequLhYERERkqSBAwdq2LBhjte1tW/fPv3jH//QwYMHFRUVJZvNpt69eyszM1ONGzf2ONZms+nzzz/Xzp07ZbVadfDgQZ0/f16dO3dWenq60tLSPI7PyclRUVGRhgwZom7duql58+ay2Ww6evSo8vPz9eWXX2rs2LHq06fPZd0jAAAAAAAAgIaDBLYfLV68WPPnz9fkyZPVq1cvSdLx48c1depUffPNNzKbzYqKiqrV3KtXr9bs2bOVlZWlF198UZJ06tQpvfzyy8rKytLLL7+suLg4t2NLS0s1Z84c3XnnnRo2bJjj/LZt2zRz5kytWrVK/fv3l9FoVKNGjdzOcf78eRUVFemtt95yaYuPj9dzzz2nX//617W6NwAAAAAAAAANEyVE/GTNmjWaN2+eRo8e7UheS1Lz5s1lMpkcyWK73V7jufPz8zVz5kwNGjRIgwcPdpy/4oorNHHiRP3yyy/Kzs6utsTHK6+8oocffljdu3evcr5r166aNm2awsLCtHr1ar322mse47jpppt09dVXKzIyUjExMWrfvr0eeOABzZ8/n+Q1AAAAAAAAgBpjBbYfWK1WzZkzR3a7XUOHDnVpv+qqq5SamqrVq1dr5cqVGjRokNdz22w2zZ49W2VlZcrIyHBpb9KkiQYMGKAlS5boww8/1OjRo6u079+/X8ePH1ebNm3czn/ttdeqW7du2rp1q1asWKFRo0YpKSnJbd8RI0aoW7duXscOAAAAAAAAAJ6wAtsP8vLydPLkSSUnJys+Pt5tn65du0qS3n33XdlsNq/nXrlypQ4dOuRY8exp7o8++khWq7VK2/79+1VYWKiffvqp2mtUzltRUaHt27d7HRsAAAAAAAAAXA4S2H6wZMkSSVLbtm2r7XPNNddIkg4dOqStW7fWeO42bdpUuwlk5dwnT57UV199VaXNYDDo1KlTevTRRzVu3Djt3bvXZbzzvCUlJV7HBgAAAAAAAACXgxIidWzv3r06fvy4JKlVq1bV9ktMTHQcb9y4UT179rzk3KdOndKuXbsuObdz28aNG/Xb3/7W8fr6669X8+bNdfz4ce3Zs0cffvihnn/++SrjK+OXVG35kEqrVq3S+vXrZbFYVFRUJElKTU3V3XffrdatW1/yngAAAAAAAACgEgnsOuZccqO68iEXt+3YscOruXfu3KmKiopLzt2kSRNFR0fr7NmzLiVAoqOjNX36dM2bN09Hjx5Vv379qrTbbDZt2rRJ0oUke+/evau9zrx583TnnXdq4sSJjnNff/21Zs6cqX/961+aMGGCBg4c6NW9SVJBQYHS09O97n+xFStW1HosAAAIPZfzXFFQUKCUlBQfRgMAAADAGySw65jz6uUmTZpU2y8mJsbtGF/MXdl+9uxZWSwW2Wy2KmVBrr32Wr366qtuxy1btkwlJSVq2rSpJk+erEaNGlV7jZ49e7psUnnzzTfrwIEDeuutt/Tyyy8rKSlJHTp08Ob2AAAAAAAAADRwJLDrmMVicRxHRlb/1+3c5m2daW/ndm632+2yWCxq1qzZJef/6aef9Pe//11XXnmlZsyY4ail7U6jRo3Ut29ft20ZGRl6++23VV5errfeekszZsy45LUlKSUlhVXUAADAZy7nucKb8m4AAAAAfI9NHOtY48aNHceV5T7ccW6Ljo726dwXtzuPq05+fr6ef/55derUSfPmzfOYvJYkk8mk66+/3m1bs2bN1KJFC0nS5s2bdebMmUteHwAAAAAAAABIYNexhIQEx3F5eXm1/crKyhzHnupZ12Zu5/mjoqJkMBg89l21apVeeOEF3X333crJyfEqnqZNm3psr1zxbbPZVFhYeMn5AAAAAAAAAIASInWsdevWjuPS0tJq+zm3tWnTxqu5nft5mruiokJWq9UlnovZ7Xa9/fbbWrFihXJyctSlSxev4pAu1OPevXu3OnXqVCWxXiksLMxx7JysBwAAAAAAAIDqkMCuY6mpqYqIiJDNZvO4OaNzm7c1Ftu3b6/4+HiVlJR4nPvEiROOEiJpaWlu+1RUVOiVV17RyZMnNW/ePMXGxrr0+eKLLzRo0CCX83v27NHTTz+t06dP66qrrtLbb79dZZNIqWpd78pyIgAAAAAAAADgCSVE6pjBYFC3bt0kSYcOHaq234EDByRJ4eHh6t27t1dzh4WFqU+fPpec++DBg47j6jZa/POf/6ywsDBNmzbNbfJaktavX+/2fF5enk6fPi1JOnLkiGO1t7Pi4mJJUvPmzdW8efNqYwUAAAAAAACASiSw/eDee++VJBUUFFS72eKuXbskSenp6WrZsqXXc48cOVIRERE6fPiwLBaLx7m7d+/udqPFzZs3a/ny5crKyqpS6sNZeXm5vv/+e7dtJ06ckHQhoX7fffe5JMB37dqls2fPSpIyMjKqvQYAAAAAAAAAOKOEiB9069ZNvXv31nfffaft27c7VmRXstlsWr9+vQwGg8aMGVOlraioSJMnT9aRI0c0cuRIjR49ukp7UlKShg0bpsWLF2vdunUaPHiwy/XXrl2ryMhIPfzww27j++c//ym73a5vv/3W7WaMp0+f1ooVKxyJ6ou1bdtWGzZs0IMPPqhRo0a5tK9evVqSlJiYqMzMTLdzAAAAAAAAAMDFSGD7ycSJE/XUU09pwYIF6tKlS5Ua0Z9//rmKi4tlNpvVqlWrKuPWrFmjffv2SZIWLVrkksCWpHHjxmnfvn1auHCh+vfvryZNmjjaNmzYoPz8fE2YMEEdO3Z0G9vmzZtVVlamnJycWt3bHXfcoaVLl7rdnLGoqEifffaZEhISNHnyZBkMhlpdAwAAAAAAAEDDQwLbTwwGg2bNmqXZs2drwoQJuu222xQTE6MtW7Zo7dq1mjZtmnr06OEyrl+/flq+fLmOHDmie+65x+3ckZGRmjp1ql577TU99dRTGjFihBISEvTDDz9o+fLlev7555Wenu52bGlpqc6cOXNZ99aqVSvNmjVLr776qk6fPq0uXbooOjpa+/btU25urlJTUzV+/HhdeeWVl3UdAAAAAAAAAA1LmN1utwc6iIbmxx9/1I4dO1RWVqbk5GSlpaUpKirKJ3P//PPP2rp1q6xWqxITE9WrVy+/rXquqKhQfn6+fvrpJ509e1YtW7ZU586dXVaVe6Nnz56SpE2bNvk6TAfnWtx8GwAAAE/88WwC8HUGAHWPXADgvWB5NmEFdgB06NBBHTp0qJO5k5OTlZycXCdzX0p4eLi6dOmiLl26BOT6AAAAAAAAAEJLeKADAAAAAAAAAADAHRLYAAAAAAAAAICgRAIbAAAAAAAAABCUSGADAAAAAAAAAIISCWwAAAAAAAAAQFAigQ0AAAAAAAAACEqRgQ4ACJTs7OxAhwAAAAAAAADAAxLYaLCmTJkS6BAAAAAAAAAAeEACGwAAAAAAAA0Cn8YG6h8S2ACAkOP8CQs+bQEAAACgEr8fAPUPCWwAQMiZOnWq45gHVAAAgPqBRQgAAHdIYAMAAAAAgIBjEQIAwJ3wQAcAAAAAAAAAAIA7JLABAAAAAAAAAEGJBDYAAAAAAAAAICiRwAYAAAAAAAAABCUS2AAAAAAAAACAoEQCGwAAAAAAAAAQlEhgAwAAAAAAAACCEglsAAAAAKiFo0ePemwvLi7W+fPn/RQNAABAaIoMdAAAAAAAUFcsFosWLlyo77//XlFRUSorK1OHDh00evRotWjRotbzHjt2TGPGjFFaWppuueUWpaSkKCYmRqdPn9aBAwe0YcMG/ec//9HMmTMVGxvrwzsCAABoWEhgAwAAAAhJBw8e1IQJE9SnTx/96U9/UlRUlGw2m95//32NGzdO06dPV+fOnWs19/nz52W327VhwwZt2LDBpb1Hjx6aPn06yWsAAIDLRAIbAAAAQMgpLS2VyWTSFVdcofHjxys8/EL1xIiICI0ZM0bbt29Xdna2Xn/9dSUmJtbqGk2bNlVqaqq2bNkiq9Wq+Ph4tW/fXhkZGerXr5/CwsJ8eUsAAAANEglsAAAAoAE5d+6cli9frj179shut6tz58669dZbFRnp+qvBsWPH9Mknn6hDhw667rrr1LJlywBEXDtvvPGGCgsL9fvf/96RvHY2ZMgQ5eTk6H/+5380e/bsWl0jJiZGL7300uWGCgAAAA/q7SaOGRkZgQ4BAAAAqFfmzZun1q1b684779QzzzyjZ599VhkZGWrdurX++te/uvRv3LixysvLtWjRIvXv318JCQnq3bt3ACKvmcLCQn3xxReSpK5du7rtU3l+y5Yt2rZtm99iAwAAQM3U2xXYP/74Y6BDAAAAAOqNl156STNmzJDdbndpO3r0qLKysrRkyRItWrRIBoNB0oUSGX/4wx8kSTt27FDXrl3d1nsONsuWLZPNZlNYWJiSk5Pd9rnyyivVtGlTnTx5UkuWLKk20Y2GYcqUKW6PAQBA4NXLBPbZs2d18ODBQIcBAAAA1AvfffedZsyYoWbNmunRRx/Vb37zG1199dU6c+aMduzYocWLF+uTTz7R0qVLdeutt+pf//qX4uPjq8zRqFGjwARfCxs3bpQkxcXFqUmTJtX2S0xM1MmTJ7V582bZ7fZa1awuLS1VXl6eDh48qGPHjunIkSNq2bKlbr75Zt1+++1q3Lhxre8D/jN16lTHMQlsAACCi18T2B988IHOnTt32fN88803On/+vA8iAgAAAELf3LlzdeONN2rp0qUuiemuXbvq3nvv1d69e/XUU0/p008/1eDBg/XFF18oNjY2MAFfhjNnzuinn36SJJd7vVhle0lJiX7++We1bdu2RtcqLS3VzJkz9dhjj6l169aO68+dO1d/+9vftGzZMs2YMUMtWrTwes6CggKlp6fXKA5nK1asqPVYAAAQei7nuaKgoEApKSk+jKZ2/JrA/vDDD5WXl8du3AAAAIAfffvtt1q1apXHhO4111yjTz75RO+++67Gjx+voUOH6l//+peaNm3qv0B94MSJE44yKZ5WX0sXNmGsdPz48RonsM+cOaNx48Y5kteV13z22We1a9cu7du3T1OnTtXrr7/udiNJAAAAXJpfE9i5ubmaPHmypk+friuvvNJRW6+mjh8/LqvV6uPoAAAAgNDUpEkTXXPNNV71HTNmjNLS0jR8+HANHjxYy5YtU1xcXB1H6DsWi8VxHBnp+dcd53bncd5o1KiRkpKS3Ca9IyIiNHjwYM2fP1+7d+/W6tWrNWDAAK/mTUlJYRU1AADwmct5rujZs6cPI6k9v9fAnjZtmtavX6+bbrpJL730Uq3m+OGHH9SlSxcfRwYAAACEJueVxt7o2LGj1q9fr8zMTN1888364osv6igy33OuOV1RUeGxr3N7TWtVN2/eXHPnzq22vWPHjo7jNWvWeJ3ABgAAQFUB2cRxwoQJmjVrVq3Hd+rUSc2aNfNhRAAAAEDoGjx4sP7v//5PDz74oCTpP//5j/75z39qxIgRateundsx8fHx+te//qWnnnpK/fr105///OfLjmP27Nl1urp4xYoVSkhIcLwuLy/32L+srMxx7DzOW57KqzjPd+DAgRrPDf+5ePW9xWKpV586AAAg1AUkgd2vXz899dRTlzVHTTZCAdDwOO8ez07yAICGbuLEierXr5/CwsL0wAMP6Le//a1++uknvfnmm8rPz692XEREhObOnas33nhDo0aNuuw4xo0bp/vvv/+y5/EkPj5eTZs21cmTJ1VaWuqxb2V7eHi4kpKSanyt3bt3q7S0VGlpaS77/DjXvHZOlCO4mM1m5eTkVDmXlJQko9Eok8kUoKgAAICzgCSwGzdurLFjx17WHJezghtA6Js6darjmAQ2AKChi42N1ZIlSzRy5Ejl5OTop59+kt1u15kzZ7wa//vf/14dO3bU3XffrRMnTtQ6jri4uDpf2RoeHq7u3btr9erVKi4ult1ur3YT+ePHj0u6UHc6Nja2RtdZsGCB3nnnHUnSvffe61jdXqmkpMRxzOKb4GQ2mzVp0iSX81ar1XGeJDYAAIEXsK2wn3zyycsaP2zYMB9FAgAAAIS+Nm3a6Ntvv9Xrr7+u559/Xo899pg++eQTr8fffPPN2rBhgx5++OE6jNI3brrpJknSuXPn9Msvv7jtU15ersOHD0uS+vbtW6P5bTab3n//fcfrPXv2uPRxTvR36NChRvOj7lksFpeV1xfLycm55Cp+AABQ9+p0BfYXX3yh3Nxcde7cWbfccou6du1al5cDAEnUMQQAwNNz+JAhQzRkyJBazduuXTu98cYbvgqzztxyyy165513VFhYqN27d6tly5YufQoKCmSz2dS0aVPdcccdNZq/tLTUUV87Pj5eY8aMcenz/fffS5LCwsKUkZFRi7tAXcrNzZXVavXYx2q1Kjc3Vw899JCfogIAAO7UaQL7nnvuqfIv1vHx8erfv78GDBigW265Rd27d6/243wAUBvUMQQAgOfwiIgIjR07VtOmTdN3332n/v37u/T57rvvJEljxoyRwWBwaZ8zZ46+/PJLXXvttfrjH/+oK664wtEWHx+v+Ph4lZaW6i9/+YtLgryiokLffPONJOmuu+7S1Vdf7cvbgw8UFRX5tB8AAKg7dVpCZOjQobLb7Y4/xcXF+vTTT/Xcc88pLS1NzZo107Bhw/TKK69o48aNqqioqMtwAIS4yjqGF6+mqaxjaDabAxQZAAD+xXO41L9/f40ePVr//ve/tXfv3iptv/zyi5YsWaIhQ4borrvuchlbWlqqpUuX6uzZs9qxY4e2bNlSpT0sLEyjRo1SRUWFbDaby/iPPvpIJ06cUO/evS977x/UjcTERJ/2AwAAdSfMbrfb62pym82mpUuX6ssvv9TXX3+trVu3ujzgOa/8uOKKK3TTTTc5VoakpaUpIiKirsJDEOvZs6ckadOmTQGOBPWFxWJRUlKSx4+CGgwGFRYW1niTJtQ/zv9vqcP/zQFoQOrbswnP4f/1/vvv67PPPtMdd9yhdu3aqaioSHl5eerbt6/Gjh1b7X16WoFd6YMPPtAXX3yhjIwMJSUlyWazae3atVq7dq0yMzM1cuTIGq10r29fZ/UZz47BiWc4AAguwfJsUqcJ7IudOnVK3377rb7++mt9/fXX2rBhg86dO1c1IKf/YcXExKhv3766++679f/+3//jwaEBCZZvENQf8+fP92qF0/z586lj2ADwyw8AX6vvzyYN/Tm8uLhYGzZs0LFjxxQfH6+0tDS3dbFrw2q1auPGjTp69KgiIiKUlJSk1NRURUVF1Xiu+v51Vt9UfnqvOtOnT6cEnZ/xDAcAwSVYnk38msC+2Llz57Ru3TrHg/TatWtd/gW88n9g0dHRev7552U0Gmv1MIj6JVi+QVB/TJ8+XS+99JJX/fhFJPTxyw8AXwu1ZxOew4NTqH2d1QeV+6c4f/0bDAb2TwkQnuEAILgEy7NJndbAvpTGjRvr5ptv1qRJk7R8+XKVlJRo7dq1uummmy4EFx7uqNt35swZTZs2TV27dtVPP/0UyLABBCHqGAIA4D2ew4ELTCaTCgsLq5wrLCwkeQ0AQBAJaAL7YhEREfr1r3+tf//737r//vuVlpampUuX6vHHH1ezZs1kt9v1448/ql+/ftq8eXOgwwUQRDIzM2UwGDz2MRgMyszM9FNEAADUHzyHoyG7uEROfS+ZAwBAqAmqBHaliIgIzZ8/X+Hh4dq7d6/mzp2rwsJCvfHGG0pKStLRo0c1bNgwHTt2LNChAggScXFxMhqNHvsYjUZ+IQEAwAOewwEAABBsgjKBLV14eDaZTPrf//1fSVJUVJTGjh2rgoIC/f73v1dRUZEefPDBAEcJIJiYTCZNnz7dZSW2wWCg9jXcmjJliuMPAOACnsMBAAAQTOp8E8ft27frr3/9qw4cOKDu3btr9OjR6tSpk1djjxw5ol/96lcuG8pI0iuvvKIXXnhBX375pW6++WZfh40AC5Yi8aifSktLFRcX53htsVhYed3AeLsBEBsFAfBWfXw24Tm8/qmPX2ehgmcC/3BeNOBuAQHvAwAEl2B5NqnTBPY333yjQYMGqaysrMr5G2+8UQ8++KBGjhxZJcl0sfXr12vo0KE6ceKE2/bf/e53iomJ0VtvveXLsBEEguUbBPUXD78NGwlsAL5W355NeA6vn+rb11ko4ZnAPy7198z7AADBJVieTeq0hMgLL7wgSbryyisdu5jb7XatX79ejz76qFq1aqU777xT8+bN0+7du2Wz2RxjN27cqLFjx+q6666rdv7x48dr9erVdXkLAAAAQL3DczgAAABCRWRdTv79999r27ZtSklJ0bFjx/TZZ59pwYIF+ve//y273a5z587p008/1aeffipJCg8PV4sWLWS1WnXq1ClJ0muvvVbt/F26dNGhQ4fq8hYAAACAeofncKBmsrOzAx0CAACoRp0msFu0aKGzZ89KurD644EHHtADDzygffv26a233tJ7772n//znP47+NptNhw8fdrweMmSIHn/8ccfrTz/9VP/7v/+rtLQ0Pf/88yovL3fZrA0AJH4JAQA0bDyHAzXDhs4AAASvOq2BPXXqVK1fv17//Oc/FRER4bbPli1b9Pnnn2vdunXat2+fzp8/r+TkZN199916+OGHFR7+3yonzZs3V0lJiSRpzJgxGjlypP7yl79oyZIldXULCJBgqbEDoH6iBjYAX6tvzyY8h9dP9e3rDKgpamADQP0SLM8mdZrAPn/+vO666y4dOXJE//M//6O+ffte1nxdunRRfn6+wsLC1KRJE3Xq1EmzZ8/WwIEDfRQxgkWwfIMAqJ9IYAPwtfr2bMJzeP1U377OgJoigQ0A9UuwPJvU6SaOkZGR+uSTT/TII4/orrvu0i233KJ58+ZVu5v5pfzf//2ffvWrX8lut6tRo0b6wx/+wEMzAAAAcBGewwEAABAq6rQGtnThX1DHjh2rkSNH6r333tO7774ri8WiiRMn1niuXr166aefftKxY8fUrFmzKh9rBAAAAPBfPIcDAAAgFNRpCRGgtoLlIwoA6idKiADwNZ5N4A98nSHUUUIEAOqXYHk2YekEAAAAAAAAACAokcAGAAAAAAABZbFYPL4GADRcJLABAAAAAEDAmM1mJSUlVTmXlJQks9kcoIgAAMGkzjdxBAAAAAAAcMdsNmvSpEku561Wq+O8yWTyd1gAgCDCCmwAAAAAAOB3FotFOTk5Hvvk5OSotLTUTxEBAIIRCWwAAAAAAOB3ubm5slqtHvtYrVbl5ub6KSIAQDAigQ0AAAAAAPyuqKjIp/0AAKGJGth+tnPnTi1cuFDFxcWKiIiQJA0cOFDDhg1zvK6tffv26R//+IcOHjyoqKgo2Ww29e7dW5mZmWrcuPElx585c0a5ublat26dIiMjVVZWpuTkZI0aNUpt27b1KoZ169YpLy9Pp0+fVlhYmCIjI3X77bfr1ltvvax7AwBvudvBPi4uLkDRAAAAoDqJiYk+7QcACE2swPajxYsXa+LEiRo+fLj+/Oc/67XXXtPkyZO1atUqTZw4UWVlZbWee/Xq1crKylJaWpr++te/6k9/+pNycnL0448/KisryyWhc7GSkhJlZWXpP//5j2bNmqU//elP+utf/6ru3btr/Pjx+vbbby8Zw5tvvqlXXnlFDz/8sObOnavXX39dzzzzjN5++23NnDlTdru91vcHAN6oyQ727hLdAAAA8J/MzEwZDAaPfQwGgzIzM/0UEQAgGJHA9pM1a9Zo3rx5Gj16tHr16uU437x5c5lMJm3btq3WSd78/HzNnDlTgwYN0uDBgx3nr7jiCk2cOFG//PKLsrOzq02Ql5WVKTs7W8XFxZo4cWKVB4ihQ4fqlltu0YwZM7Rr165qY8jLy9OiRYv0+OOPq0OHDo7zrVu31oQJE7Ry5Uq98cYbNb43APBW5Q72F9dRrNzB3jmJXZNENwAAAOpGXFycjEajxz5Go1GxsbF+iggAEIxIYPuB1WrVnDlzZLfbNXToUJf2q666SqmpqVq9erVWrlxZo7ltNptmz56tsrIyZWRkuLQ3adJEAwYM0M6dO/Xhhx+6neODDz5Qfn6+Bg4cqOjoaJf2jIwMnTt3TrNnz5bNZnNpLyoq0htvvKErrrhC/fv3d2nv0qWLWrdurdzcXO3YsaNG9wcA3qjJDvY1SXQDAACgbplMJk2fPt1lJbbBYND06dNlMpkCFBkAIFiQwPaDvLw8nTx5UsnJyYqPj3fbp2vXrpKkd999122SuDorV67UoUOHFBMTo/bt23uc+6OPPnKbsMnLy5Mk3XDDDW7HX3fddYqOjtaBAwe0atUql/bKmDt37lxtHe/KGN5++23vbgwAasDbHezfeecdrxPdAAAA8A+TyaTCwsIq5woLC0leAwAkkcD2iyVLlkiSx40Qr7nmGknSoUOHtHXr1hrP3aZNm2qTx5Vznzx5Ul999VWVtq+++kqnTp2SJLVr187t+IiICLVp00aS9Nlnn1VpO3PmjL788kuP451j2Lp1qw4ePOjhjgCg5rzdmX7lypVeJbpzc3N9ERYAAAC8dHGZEMqGAAAqkcCuY3v37tXx48clSa1ataq2n/Ouyhs3bvRq7lOnTjnqUnua27nt4rk3bNjgtl91c+zatcuR8JYuJKTLy8svOb429wcA3vL1zvTeJsQBAAAAAEDdigx0AKFu+/btjuPqyodc3OZtneidO3eqoqLiknM3adJE0dHROnv2bJV4nK/VpEkTNW7c+JLxVVRUKD8/XzfeeKNLrDW5vzvvvLPavpUKCgqUnp5+yX7VWbFiRa3HAqhfMjMz9eSTT3pcXW0wGPSb3/xGixcvvuR8vk6IAwgOl/NcUVBQoJSUFB9GAwAAAMAbrMCuY5Wrr6ULSeLqxMTEuB3ji7md2y0Wi6PG9vnz52WxWLwaX118zsfuNoB0F5+39wcA3vJ2B/v77rvPZYOgixkMBmVmZvoyPAAAAAAAUEuswK5jlQliSYqMrP6v27mtpKTEp3M7t9vtdlksFjVr1kwWi0V2u71G4y++rnOsjRo1qna8c5vzeE9SUlJYRQ3Aa5Wb/OTk5FRZiW0wGGQ0Gh3tRqNRkyZNqnYeo9FIzUUgRF3Oc0XPnj19GAkAAAAAb5HArmPOZTkqy32449zmaSVzbea+uL1ynPN1ajO+JnNUrvq+eDwA+JLJZFJWVpbi4uIc5woLC6skpL1NdAMAAAAAgMAjgV3HEhISHMeVmx26U1ZW5jj2VEu6NnM7zx8VFeX4+LzBYFCjRo1UXl7u9fiLr+scq6c5nNucxwOAr3mzg703iW4AAAAAABB41MCuY61bt3Ycl5aWVtvPua1NmzZeze3cz9PcFRUVjlWGzvE4z3Hy5ElHOZFLxec8h7cxVDceAALFm0Q3AAAAAAAILBLYdSw1NVURERGSPG9e6NzmbY3F9u3bO1ZAe5r7xIkTjvIeaWlpVdoqr1VRUaHi4uJLxpeQkKD27du7jdXX9wcAAAAAAACgYSOBXccMBoO6desmSTp06FC1/Q4cOCBJCg8PV+/evb2aOywsTH369Lnk3AcPHnQc9+3bt0qb82tv5ujTp4/CwsIc59u1a6ekpKRLjq+8P+e/DwAAAAAAAADwhAS2H9x7772SpIKCgmo3Oty1a5ckKT09XS1btvR67pEjRyoiIkKHDx+WxWLxOHf37t11/fXXV2nr0qWLunbtWqXfxU6cOKGjR48qIiJCI0eOdGkfNWqUx/HObXfffbfXm1QCAAAAAAAAaNhIYPtBt27d1Lt3b1mtVm3fvt2l3Wazaf369TIYDBozZkyVtqKiIj3yyCMaPny43nvvPZexSUlJGjZsmOx2u9atW+f2+mvXrlVkZKQefvhht+1jx45VZGRkteO/++47SdLw4cN19dVXu7QPGjRI1157rfbv36+ioiKXdqvVqh07dqhFixa666673F4DAAAAAAAAAC5GAttPJk6cqHbt2mnBggWy2WxV2j7//HMVFxdr8uTJatWqVZW2NWvWaN++fTpz5owWLVrkdu5x48ape/fuWrhwoc6cOVOlbcOGDcrPz9fTTz+tjh07uh3fqVMnPfHEE/r++++1efPmKm2nT5/WokWL1KNHD40bN87t+IiICE2bNk3NmjXT22+/7dL+/vvvKzw8XGazWQaDwe0cAAAAAAAAAHCxMLvdbg90EA3FiRMnNHv2bJWVlem2225TTEyMtmzZorVr1+qJJ55Qr169XMYUFRVp8uTJOnLkiH73u9+5rNCudPr0ab322mvat2+fRowYoYSEBP3www9avny5HnzwQaWnp18yvmXLlumdd97RkCFD1LFjRxUXFysvL0/XXnutsrKy1KRJE4/jDx48qJycHDVv3ly/+c1vFBkZqW+//VY//PCDXnjhBV133XXe/UXpvxs9btq0yesxAFDJuVa/p//NedsPAHg2gT/wdYZQd6lnL57NACC4BMuzCQnsAPjxxx+1Y8cOlZWVKTk5WWlpaYqKivLJ3D///LO2bt0qq9WqxMRE9erVq0arnq1Wq9avX6/Dhw/LYDCoR48eat26tdfj7Xa7tm/frt27d6uiokLt27dXamqqIiIianQfwfINAqB+IoENwNd4NoE/8HWGUEcCGwDql2B5NokM6NUbqA4dOqhDhw51MndycrKSk5NrPd5gMGjgwIG1Hh8WFqauXbs6NoYEAAAAAAAAgNqiBjYAAAAAAAAAICiRwAYAAAAAAAAABCUS2AAAAAAAAACAoEQCGwAAAAAAAAAQlEhgAwAAAAAAAACCEglsAAAAAAAAAEBQIoENAAAAAAAAAAhKkYEOAAAAAAAAAIA0ZcoUt8dAQ0YCGwDQYGVnZwc6BACAH+3evVt//OMf9e677/pszjNnzig3N1fr1q1TZGSkysrKlJycrFGjRqlt27Y+uw4AoGGYOnWq45gENnABCWwAQIPFAyEANCwLFy706XwlJSV67rnn1Lp1a82aNUsGg0GStGzZMo0fP15Go1F9+/b16TUBAAAaGmpgAwAAAAhpxcXFeuWVV7RmzRqfzVlWVqbs7GwVFxdr4sSJjuS1JA0dOlS33HKLZsyYoV27dvnsmgAAAA0RK7ABAAAAhJz9+/fr9ddfV1lZmUpKSnT8+HGfzv/BBx8oPz9fd9xxh6Kjo13aMzIytGzZMs2ePVtvvvmmIiIifHp9AACAhoIENgAAAICQ07ZtW7366quO188++6y2bdvmk7mtVqvy8vIkSTfccIPbPtddd52io6N14MABrVq1Sunp6T65NgAAQENDCREAAAAAqIGvvvpKp06dkiS1a9fObZ+IiAi1adNGkvTZZ5/5KzQAAICQQwIbAAAAAGpgw4YNjuNWrVpV26+ybdeuXY6ENwAAAGqGEiIAAAAAUAM7duyQJDVp0kSNGzeutl98fLwkqaKiQvn5+brxxhsvOXdBQcFllRtZsWJFrccCAIDQcznPFQUFBUpJSfFhNLXDCmwAAAAA8NL58+dlsVgkXUhgexITE+M49vUmkgAAAA0FK7ABAAAAwEsWi0V2u12SFBnp+dcp5/bKpPelpKSksIoaIeni7wGLxaK4uLgARQMADcflPFf07NnTh5HUHiuwAQAAAMBL0dHRjuOKigqPfZ3bPZUaAUKd2WxWUlJSlXNJSUkym80BiggAUJ+wAhsAAACAX8yePbtOVxf7Y+WywWBQo0aNVF5ervLyco99y8rKHMcJCQl1HRoQlMxmsyZNmuRy3mq1Os6bTCZ/hwUAqEdIYAMAAADwi3Hjxun+++8PdBiXrU2bNtqzZ49Onjwpu92usLAwt/1KS0sdx61bt/ZXeEDQsFgsysnJ8dgnJydHWVlZio2N9VNUAID6hgQ2AAAAAL+Ii4sLiZq3PXv21J49e1RRUaHi4mI1a9bMbb/KjRsTEhLUvn17f4YIBIXc3FxZrVaPfaxWq3Jzc/XQQw/5KSoAQH1DDWwAAAAAqIG+ffs6jg8dOlRtv4MHD0qS+vTpU+0qbSCUFRUV+bQfEOrcbXYKgAQ2AAAAANRIly5d1LVrV0nSrl273PY5ceKEjh49qoiICI0cOdKf4QFBIzEx0af9gFDGZqdA9UhgAwAAAMBF5syZo2HDhunpp5/WqVOnXNrHjh2ryMhIrVu3zu347777TpI0fPhwXX311XUaKxCsMjMzZTAYPPYxGAzKzMz0U0RAcKrc7PTikjuVm52SxEZDRwIbAAAAQMgrKyur8l9PSktLtXTpUp09e1Y7duzQli1bXPp06tRJTzzxhL7//ntt3ry5Stvp06e1aNEi9ejRQ+PGjfPNDQD1UFxcnIxGo8c+RqORDRzRoHm72anzxsBAQ8MmjgAAAABCzv79+/X6669Lks6fP68ff/xRklRcXKxHH31UMTExkqQxY8YoNTW1ytjY2FhlZGToyy+/1LXXXuvSXikjI0NhYWF65ZVXNGTIEHXs2FHFxcXKy8tT586dlZWVpYiIiDq8SyD4mUwmSRcScM6rSw0Gg4xGo6MdaKjY7BS4NBLYAAAAAEJO27Zt9eqrr9Z6/DPPPKNnnnnmkv2GDh2q/v37a/369dqzZ48MBoMmT56s1q1b1/raQKgxmUzKyspSXFyc41xhYSErrwGx2SngDRLYAAAAAHAZDAaDBg4cGOgwgKB2cbKa5DVwAZudApdGAhsAvDBlyhS3xwAAAAB8Izs7O9AhAH6XmZmpJ5980mMZETY7RUNHAhsAvDB16lTHMQlsAAAAwPd4zkZDVLnZ6aRJk6rtw2anaOhIYAMAAAAAAAABwmangGcksAEAAAAAAIAAYrNToHrhgQ4AAAAAAAAAaOjY7BRwjwQ2AAAAAAAAACAokcAGAAAAAAAAAAQlEtgAAAAAAAAAgKBEAhsAAAAAAAAAEJRIYAMAAAAAAAAAghIJbAAAAAAAAABAUCKBDQAAAAAAAAAISiSwAQAAAAAAAABBiQQ2AAAAAAAAACAoRQY6AAAAfC07OzvQIQAAAAAAAB8ggQ0ACDlTpkwJdAgAAAAAAMAHKCECAAAAAAAAAAhKJLABAAAAAAAAAEGJEiIAAAAAAABAEGA/H8AVCWwAAAAAAAAgCLCfD+CKEiIAAAAAAAAAgKBEAhsAAAAAAAAAEJRIYAMAAAAAAAAAghIJbAAAAAAAAABAUCKBDQAAAAAAAAAISiSwAQAAAAAAAABBiQQ2AAAAAAAAACAokcAGAAAAAAAAAAQlEtgAAAAAAAAAgKBEAhsALsFisXh8DQAAAAAAgLpBAhsAPDCbzUpKSqpyLikpSWazOUARAQAAAAAANByRgQ4AAIKV2WzWpEmTXM5brVbHeZPJ5O+wAAAAAAAAGgxWYMMn7Ha7jh075rHP0aNH/RQNcPksFotycnI89snJyVFpaamfIgIAAAAAAGh4WIHtRxaLRQsXLtT333+vqKgolZWVqUOHDho9erRatGhx2fOfOXNGubm5WrdunSIjI1VWVqbk5GSNGjVKbdu2veT4nTt3auXKlSopKdGRI0d0/PhxXX311erXr59++9vfKjo6utqxX375pV599VXdeuut6t27t9q2bauIiAiVlpZq7969WrNmjWJiYvTiiy9e9n0C/pCbmyur1eqxj9VqVW5urh566CE/RQUAAAAAANCwkMD2k4MHD2rChAnq06eP/vSnPykqKko2m03vv/++xo0bp+nTp6tz5861nr+kpETPPfecWrdurVmzZslgMEiSli1bpvHjx8toNKpv377Vjn/77bcVFRWlP/zhD4qKipIkFRcX689//rPmzZunvLw8TZkyRddee63b8eXl5SorK9OyZcu0bNmyKm1hYWG68847NW7cuFrfH+BvRUVFPu0HAAAAAACAmqOEiB+UlpbKZDLpiiuu0Pjx4x0J4oiICI0ZM0YpKSnKzs6udSKsrKxM2dnZKi4u1sSJEx3Ja0kaOnSobrnlFs2YMUO7du1yO/6rr76SJI0aNcoRmyQlJCToxRdfVIcOHXTkyBG9+OKLKikpqTaOX/3qV+ratauio6MVGRmpVq1aKT09XXPnztVjjz2miIiIWt0fEAiJiYk+7QcAAAAAAICaI4HtB2+88YYKCws1ePBghYe7/pUPGTJEJSUl+p//+Z9azf/BBx8oPz9fAwcOdFvmIyMjQ+fOndPs2bNls9lc2nNzczVgwAC3c0dEROjOO++UdGFF9uLFi6uNo3379nr11Vf16aefaunSpVqwYIGef/55XXfddbW6LyCQMjMzq/xjkDsGg0GZmZl+iggAAAAAAKDhIYFdxwoLC/XFF19Ikrp27eq2T+X5LVu2aNu2bTWa32q1Ki8vT5J0ww03uO1z3XXXKTo6WgcOHNCqVatc2vft26dvvvnGbXJbupCYrlTT+ID6Ki4uTkaj0WMfo9Go2NhYP0UET6ZMmeL4A6D+43saAAAAQCVqYNexZcuWyWazKSwsTMnJyW77XHnllWratKlOnjypJUuWVJvoduerr77SqVOnJEnt2rVz2yciIkJt2rRRQUGBPvvsM6Wnp1dpNxgMeuutt/TRRx9pzJgxuuuuu1zGV/JUQgQINSaTSZKUk5NTZUNHg8Ego9HoaEfgTZ061XFMwguo//ieBgAAAFCJBHYd27hxo6QLqzmbNGlSbb/ExESdPHlSmzdvlt1uV1hYmFfzb9iwwXHcqlWravu1atVKBQUF2rVrl06dOqUrrrjC0TZgwADl5eXp5MmT+vvf/67hw4crMvK/XxrHjx93HCclJXmMZ//+/frkk0904sQJHT16VMePH9e1116roUOHql+/fl7fFxAsTCaTsrKyFBcX5zhXWFjIymsAAAAAAAA/IIFdh86cOaOffvpJkhQfH++xb2V7SUmJfv75Z7Vt29ara+zYsUOS1KRJEzVu3PiS81dUVCg/P1833nijo+3BBx9UeXm5NmzYoG7dulVJXkvS+vXrHccjRoyo9hrbt29XRUWFsrKyHLWDDx8+rGnTpmnatGnq37+/XnzxRZf5q1NQUOCyWrwmVqxYUeuxgLOLk9UkrwGgfrqc54qCggKlpKT4MBoA8C2LxaLc3FwVFRUpMTFRmZmZVRZhAABQX1EDuw6dOHFCdrtdkjyuvpakmJgYx7HzimdPzp8/L4vFctnzR0dH64knntCCBQv03HPPVWk7ceKEli9fLkn6//6//089evSo9hoVFRV65plnqmx816pVK7300kuKjIzU6tWr9b//+79e3RsAAAAAwDtms1lJSUkaO3asXnrpJY0dO1ZJSUkym82BDg0AgMvGCuw6VJlclnTJVcfO7c7jLjV/ZYK8LuYvKyvTrFmzZLVa9fvf/1733HNPtX0bNWqktLQ0RUVFubQlJiYqNTVVGzZs0CeffKK7775bLVq0uOT1U1JSWEUNAAB85nKeK3r27OnDSADAd8xmsyZNmuRy3mq1Os6zdwsAoD5jBXYdci7pUVFR4bGvc7unUiDOoqOj62z+06dP66WXXtKePXs0Y8YMj8lrSbr55pv1hz/8odr26667TtKFVePr1q275PUBAAAAAJ5ZLBbl5OR47JOTk6PS0lI/RQQAgO812BXYs2fPrtPVvStWrFBCQoLjdXl5ucf+ZWVljmPncZ4YDAY1atRI5eXlPp3/yJEjeumllxQbG6t58+bpyiuvvGQskZGRHleBN2vWzHF84MCBS84HAAAAAPAsNzdXVqvVYx+r1arc3Fw99NBDfooKAADfarAJ7HHjxun++++v02vEx8eradOmOnny5CX/xbuyPTw8XElJSV5fo02bNtqzZ49Onjwpu92usLAwj/NLUuvWraudb+fOnZo2bZpuv/12jR49WuHh3i3St9ls2rx5s+Li4tShQweXdue4nJPpAAAAAIDaKSoq8mk/AACCUYNNYMfFxdX5jszh4eHq3r27Vq9ereLiYo8J5sqNFVNSUhQbG+v1NXr27Kk9e/aooqJCxcXFVVY6u5s/ISFB7du3d9tn7dq1+tvf/qYXX3xR3bp1c2kvKSnRTz/9pLS0tCrn7Xa7Jk6cqK1bt0qSZs6c6VInsqSkxHHsTf1rAAAAAIBniYmJPu0HAEAwogZ2HbvpppskSefOndMvv/zitk95ebkOHz4sSerbt2+N5nfuf+jQoWr7HTx4UJLUp08ft0n0HTt26G9/+5tefvllt8lrSfrpp58ciXBn27dvdySvJWnv3r0ufU6cOOE4TklJqTZOAAAAAIB3MjMzZTAYPPYxGAzKzMz0U0QAAPgeCew6dsstt+jqq6+WJO3evdttn4KCAtlsNjVt2lR33HFHjebv0qWLunbtKknatWuX2z4nTpzQ0aNHFRERoZEjR7q0l5eXa9asWbr33nvVsmXLaq+1ceNGt+eLi4sdxx06dFB6erpLn++//16S1KpVK5fV2QAAAEB9dPToUY/txcXFOn/+vJ+iQUMUFxcno9HosY/RaKzRp3wBAAg2JLDrWEREhMaOHStJ+u6779z2qTw/ZswYt/96PmfOHA0bNkxPP/20Tp065dI+duxYRUZGat26dR7nHz58uCOZ7mzDhg06fPiwdu/erQ0bNrj8Wbt2rRYsWKB//vOfbudPTk6WJLVr105/+tOfXEqz7N+/Xz///LMk6dFHH/W6rjYAAADgS7t379aYMWN8MtexY8c0ZswYvfjii1q+fLn27t2rI0eOaO/evfr666/16quvymQy6fTp0z65HlAdk8mk6dOnu/wuaTAYNH36dJlMpgBFBgCAbzTYGtj+1L9/f40ePVoffvihMjMzdc011zjafvnlFy1ZskRDhgzRXXfd5TK2tLRUS5culXShzMeWLVvUv3//Kn06deqkJ554QnPmzNHmzZvVo0cPR9vp06e1aNEi9ejRQ+PGjXMb36ZNmyRJn376qT799NMa398111yjvn37at++fS7J6YqKCs2fP19hYWF64IEHalwiBQAAAPCVhQsX+myu8+fPy263OxZ9XKxHjx6aPn06K1/hFyaTSVlZWbrvvvt08uRJNW3aVAsWLODrDwAQEkhg+8kDDzygxo0ba9KkSbrjjjvUrl07FRUVKS8vT4MHD3as0r5YbGysMjIy9OWXX+raa69Vamqq234ZGRkKCwvTK6+8oiFDhqhjx44qLi5WXl6eOnfurKysLEVERLgde6mPPnpj4sSJev311/XCCy8oPT1dzZo108mTJ/X555/rxIkTmjp1qvr06XPZ1wEAAABqqri4WPPnz9eaNWt01VVX+Wzepk2bKjU1VVu2bJHValV8fLzat2+vjIwM9evXr9oN3IG6EBsbq08++STQYQAA4HNhdrvdHuggGpLi4mJt2LBBx44dU3x8vNLS0jzWna4pq9Wq9evX6/DhwzIYDOrRo4dat27ts/kv5ejRo9qyZYuKi4t1xRVXqF27drr++utr/PBeWSe7cnU4EGjOX8P82Aw+vD9AaAnG72meTeqf/fv36/XXX1dZWZlKSkp0/PhxlZeX66qrrtK777572fMfPnxYzz33nE/mqsTXGUJdMP58BwBUL1ieTViB7WcJCQkaPHhwnc1vMBg0cODAOpv/Ulq2bKkhQ4YE7PoAAACAJLVt21avvvqq4/Wzzz6rbdu2BTAiANnZ2YEOAQBQD5HABgAAAAAAdW7KlCmBDgEAUA+RwAYAAACAWiotLVVeXp4OHjyoY8eO6ciRI2rZsqVuvvlm3X777WrcuHGgQwQAAKjXSGADAAAAQC2UlpZq5syZeuyxxxz7zpw5c0Zz587V3/72Ny1btkwzZsxQixYtvJ6zoKBA6enptY5pxYoVtR4LAABCz+U8VxQUFCglJcWH0dROeKADAAAAAID66MyZMxo3blyVTdObNGmiZ599Vm3atNG+ffs0depUVVRUBDBKAACA+o0V2AAAAABQQ40aNVJSUpLatm3r0hYREaHBgwdr/vz52r17t1avXq0BAwZ4NW9KSgqrqAEAgM9cznNFz549fRhJ7bECGwAAAABqqHnz5po7d2617R07dnQcr1mzxh8hAQAAhCRWYAMAAADwi9mzZ9fp6mJ/r1xu2rRptW0JCQmO4wMHDvgjHAAAgJBEAhsAAACAX4wbN073339/oMPwmd27d6u0tFRpaWkKCwur0hYe/t8Pu5aVlfk7NAAAgJBBAhsAUK9ZLBaX13FxcQGKBgDgSVxcXMj8jF6wYIHeeecdSdK9996rBx98sEp7SUmJ47hFixb+DA0AACCkUAMbAFBvmc1mJSUlVTmXlJQks9kcoIgAAA2BzWbT+++/73i9Z88elz4nTpxwHHfo0MEvcQEAAIQiVmADAOols9msSZMmuZy3Wq2O8yaTyd9hAQAagNLSUpWXl0uS4uPjNWbMGJc+33//vSQpLCxMGRkZfo0PAAAglLACGwBQ71gsFuXk5Hjsk5OTo9LSUj9FBAAINXPmzNGwYcP09NNP69SpU1Xa4uPjFR8fr/DwcP3lL3/RddddV6W9oqJC33zzjSTprrvu0tVXX+23uAEAAEINCWwAQL2Tm5srq9XqsY/ValVubq6fIgIABLvKjRS92VCxtLRUS5cu1dmzZ7Vjxw5t2bKlSntYWJhGjRqliooK2Ww2l/EfffSRTpw4od69e2vs2LG+uQEAAIAGihIiAOCF7OzsQIcAJ0VFRT7tBwAIPfv379frr78uSTp//rx+/PFHSVJxcbEeffRRxcTESJLGjBmj1NTUKmNjY2OVkZGhL7/8Utdee61Lu3RhZXV5ebkmT56sjIwMJSUlyWazae3atVq7dq0efvhhjRw5UmFhYXV8pwAAAKGNBDYAeGHKlCmBDgFOEhMTfdoPABB62rZtq1dffbXW45955hk988wzHvuMHDlSt99+uzZu3Kiff/5ZERER6tevn7KyshQVFVXrawMAAOC/SGADAOqdzMxMPfnkkx7LiBgMBmVmZvoxKgBAQ2QwGDRgwIBAhwEAABCyqIENAKh34uLiZDQaPfYxGo2KjY31U0QAAAAAAKAusAIbAFAvmUwmSVJOTk6VldgGg0FGo9HRDgAAAAAA6i8S2ACAestkMikrK0txcXGOc4WFhay8BgAAAAAgRFBCBABQr12crCZ5DQAAAABA6CCBDQAAAAAAAAAISiSwAQAAAAAAAABBiQQ2AAAAAAAAACAokcAGAAAA/v/27j0uyjL///ibo8gooKmBeErD1LQ81aZm5a4nWltNaf2lVlq27raSm1kr4TekJBKtLKvtuJu1u5myblmmZdvJLVNMU9M00jwgKCowwIgCw/z+8MEs4wwjh2G4B17Px8PH4+a+rvu6r/GSm48frvu6AAAAABgSCWwAAAAAAAAAgCGRwAYAAAAAAAAAGBIJbAAAAAAAAACAIQU2dgcAAACA+li4cKHLYwAAAAC+jwQ2AAAAfFpycrL9mAQ2AAAA0LSwhAgAAAAAAAAAwJBIYAMAAAAAAAAADIkENgAAAAAAAADAkEhgAwAAwGeZzWa3XwMAAADwbSSwAQAA4JNSUlIUHR3tcC46OlopKSmN1CMAAAAAnhbY2B0AAAAAaislJUULFixwOm+xWOznExMTvd0tAAAAAB5GAhsAAAA+xWw2KzU11W2d1NRUxcfHKywszEu9AoCmZ+HChS6PAQDwJhLYAAAA8Cnp6emyWCxu61gsFqWnp+vuu+/2Uq8AoOlJTk62H5PABgA0FtbABgAAgGHUZFPGnJycGrVV03oAAAAAjIsENgAAAAyhppsyRkVF1ai9mtYDAAAAYFwksAEAANDoKjdlvHBpkMpNGasmsePi4mQymdy2ZzKZFBcX1yB9BQAAAOA9JLABAADQqGq6KWNhYaEkKTw8XAkJCW7rJyQksIEjAAAA0ASQwAYAAECjqs2mjJUSExO1aNEip5nYJpNJixYtUmJiYoP0FQAAAIB3BTZ2BwAAANC81XVTxsTERMXHxys8PNx+Ljs7m5nXAAAAQBNCAhsA4POSkpIauwsA6qE+mzJemKwmeQ0AAAA0LSSwAQA+b+HChY3dBQD1EBcXpzlz5rhdRoRNGQEAAIDmiTWwAQAA0KjYlBEAAABAdZiBDQAAgEZXueliamqqw0xsk8mkhIQENmUEAAAAmikS2AAAADCEum7KyDr4AAAAQNNFAhsAAACGUZdNGVkHHwAAAGi6WAMbAAAAAAAAAGBIJLABAAAAAAAAAIZEAhsAAAAAAAAAYEgksAEAAAAAAAAAhkQCGwAAAAAAAABgSCSwAQAAAAAAAACGRAIbAAAAAAAAAGBIJLABAAAAAAAAAIZEAhsAAAAAAAAAYEgksAEAAAAAAAAAhkQCGwAAAAAAAABgSCSwAQAAAAAAAACGRAIbkqSCggKVlpZWW261WnX69Gkv9ggAAAAAAABAcxfY2B1oTsxms1auXKmdO3cqODhYpaWl6tmzp6ZOnar27dvXu/2SkhKlp6dry5YtCgwMVGlpqbp06aLbb79dXbt2dXttamqqcnJyNGbMGF199dW65JJLZLValZubq7179+qzzz7TzJkzNWTIELftHDp0SP/85z+VlZWl4OBgWa1WXXfddYqLi1OLFi3q/RkBAAAAAAAANB8ksL0kKytLDz30kIYMGaJly5bZk7tvv/22Zs2apUWLFqlPnz51br+goEDz5s1Tp06dtHjxYplMJknShg0bNHv2bCUkJGjo0KHVXl9eXq6cnBy98cYbTmURERGaN2+efvGLX7jtw6ZNm5SWlqb4+Hg98sgjkqTi4mItWbJE8fHxWrJkicLDw+v8GQEAAAAAAAA0Lywh4gWFhYVKTExUq1atNHv2bAUHB0uSAgICNG3aNMXExCgpKUk5OTl1ar+0tFRJSUnKz8/X/Pnz7clrSRo7dqxuuukmPfHEE9q3b5/bdoYNG6aOHTsqMDBQoaGh6tGjh6ZPn67XX3/9osnrvXv36sknn9TIkSM1evRo+/lWrVpp/vz5OnnypJKSktwuUwIAAAAAAAAAVZHA9oJXXnlF2dnZGj16tPz9nf/Kx4wZo4KCAj3zzDN1av+dd97R3r17NWLECIWEhDiVx8bG6ty5c0pLS5PVaq22nVtvvVUrVqzQ+vXr9d577+mll17S1KlTFRYW5vb+VqtVaWlpKi0tVWxsrFN5y5YtdeONN2rPnj1avXp17T8gAAAAAAAAgGaJJUQaWHZ2tj755BNJ0lVXXeWyTuX5HTt2aNeuXdXWc8VisWjNmjWSpH79+rmsc8UVVygkJERHjx7Vp59+qlGjRtXmI1zUf/7zHx07dsw+a9uVq666SuvWrdO//vUvTZgwwWGWOAAAANAQSkpKtHbtWh08eFBFRUXKysqSn5+f+vfvr7Fjx6p3794euUdd96EBAADAxTEDu4Ft2LBBVqtVfn5+6tKli8s67dq1U+vWrSVJ69atq1X7X3zxhYqLiyVJ3bp1c1knICBAnTt3liR98MEHtWq/Jir73LlzZwUEBLisc9lll0mSioqK9MUXX3i8DwAAAEBVOTk5Wrx4sQYNGqSEhAQ98cQTevPNN3Xvvffqyy+/1Jw5c/TCCy/IZrPV+R4FBQWKj4/XgQMHtHjxYi1btkwvvvii+vfvr9mzZ+vrr7/24CcCAABonpiB3cC2bdsmSQoPD1fLli2rrRcVFaWioiJt375dNptNfn5+NWo/IyPDfhwZGVltvcjISGVmZmrfvn0qLi5Wq1atXNb79NNPtXXrVpnNZvua3AMGDNCkSZPUqVMnp/rFxcX2tbUvdv9K27Zt08033+z+gwEAAAB1VFFRoSVLlmjhwoVOy+Fdf/318vf3V1JSkt59912FhYXpjjvuqPU9qu5D8/zzzzss5Td27Fjt2bNHTzzxhJYuXapevXrV+zMBAAA0VySwG1BJSYl++uknSVJERITbupXlBQUFOnLkSI1fN/z+++8lnV9nukWLFhdtv6KiQnv37tW1117rVOcvf/mLJkyYoPnz59vPffnll3ryySf10Ucf6aGHHtKIESMcrtmzZ48qKioc7uFKy5YtFRISorNnz2r37t01+myZmZn1Wu5k48aNdb4WAAA0PfWJKzIzMxUTE+PB3qAhbdu2TREREdXu5TJ06FBFRUUpJydH6enpmjhxYq2XuKvch2b8+PHV7kOzYcMGpaWl6dVXX632TUUAAAC4xxIiDSgvL8/+SqK72deSFBoaaj8+ffp0jdovLy+X2Wz2WPuDBg3S2LFjHc7dcMMNmjp1qsrKyrRkyRL9+OOPDuVV27pYHyrLzWaz280kAQAAgPo4fPiwDhw4oOzs7GrrdO/eXZJ05swZZWZm1qr9uuxDAwAAgLphBnYDqkwuS1JgoPu/6qrlVa+7WPuVCfL6th8UFKShQ4e6vDY2NlYrVqxQWVmZ3njjDT3xxBMu26ppH2w2m8xms9q2beu2fkxMDLOoAQCAx9Qnrhg0aJAHe4KGZjKZlJ2drenTp+vKK6/UI488ovbt2zvUqRq75ufn16r92uxDk5mZqQ8++MDjG6kDAAA0FySwG1DVJT0ql9moTtVyd0uBVFX1VcX6tp+YmGjfSPJCbdu2Vfv27ZWbm6vt27erpKTEPpu6oT8jAAAAUFvXXnutffm677//Xh9++KHuuusuhzqnTp2yH0dHR9eqfU/vQ1MVy+gBAABPagrL6DXbBHZaWlqDBncbN25UmzZt7F+XlZW5rV9aWmo/rnqdOyaTSUFBQSorK6t3+9Ulryu1bdtWubm5slqtys7OVo8ePZzaqmkfgoODa73GIAAAAFBT7dq1U3Jysv72t7+puLhY11xzjUN5YWGhfWm8vn37qmfPnrVq35P70AAAAMC9ZpvAnjVrlu68884GvUdERIRat26toqIiFRYWuq1bWe7v71+rGSCdO3fWwYMHVVRUJJvNJj8/P7ftS1KnTp2cyk+fPq39+/erd+/eLhPcVdutmgzv3Lmzy3tcqKKiQhaLpdr7AwAAAJ40cOBADRw40GXZ6tWrVVZWpsjISCUkJNSqXU/vQ3MhltEDAACe1BSW0Wu2Cezw8HCFh4c36D38/f3Vv39/bdq0Sfn5+W4TzJUBbUxMTLW7pbsyaNAgHTx4UBUVFcrPz692XenK9tu0aWOfPV3p4MGDeuCBB3TmzBldeumlWrFihdMu6QUFBfbjqusH9ujRQxERESooKHAblOfl5dmXEBk8eHCNPx8AAADgSdu2bdOqVavUtWtXpaamOq2NfTGe3IcGAAAAF+ff2B1o6oYNGyZJOnfunE6ePOmyTllZmY4fPy5J1W6kWJ2q9Y8dO1ZtvaysLEnSkCFDnJLoa9as0ZkzZyRJJ06csM+UrqpyY5tLLrlEl1xyif28n5+fhgwZUuP7X9hnAAAAwFu+/vprPfrooxoxYoSWL19e6+S15Nl9aAAAAHBxzXYGtrfcdNNNevPNN5Wdna39+/erQ4cOTnUyMzNltVrVunVrjR8/vlbt9+3bV1dddZV27dqlffv2qV+/fk518vLylJubq4CAAE2ePNlluXQ+GT1t2jSnGeD79u3T2bNnJUmxsbFOCfDJkyfr448/1vHjx2U2m13ObN+3b58kqX///rryyitr9RkBAADQNHhjH5rqrFq1Sn//+981e/Zs3XzzzXW+hyf3oQEAAMDFkcBuYAEBAZo5c6Yee+wxffPNNxo+fLhTnW+++UaSNG3aNJebGz799NP67LPPdPnll+vxxx932r185syZmjt3rrZs2aLbbrut2vZ/85vfqGPHjk7lXbt2VUZGhmbMmKHbb7/dqXzTpk2SpKioKMXFxTmVR0dH65ZbbtG7776rLVu2aPTo0U51Nm/erMDAQN1zzz1OZQAAAGgevLEPzYXKy8u1bNky7d+/X8uXL1fXrl3r3aan9qEBAADAxZHA9oLhw4dr6tSpWr16teLi4nTZZZfZy06ePKl169ZpzJgxmjhxotO1hYWFWr9+vaTzu53v2LHDKQneu3dv3X///Xr66ae1fft2h81qzpw5o1WrVmngwIGaNWuWy/6NHz9e69evd5ghUiknJ0cffPCB2rRpo0cffdRlgl06/5+RQ4cOaeXKlRo+fLjDhjYZGRnau3evHnroIfXq1cvN3xQAAACaMm/sQ1PV2bNnlZSUpMjISD3//PMul/H45JNPNHLkyFq164l9aAAAAFAzJLC9ZPr06WrRooUWLFig8ePHq1u3bsrJydGaNWs0evRozZw50+V1YWFhio2Ntc/AHjBggMt6lUt7LF26VGPGjFGvXr2Un5+vNWvWqE+fPoqPj3famLFSZGSkFi9erKeeekpnzpxR3759FRISokOHDik9PV0DBgzQ7Nmz1a5du2o/X2BgoJKTk/Xss8/qT3/6k2699Va1adNGP/zwgz7++GM9/PDDGjVqVO3/4gAAAIA6sNlsSklJUZ8+fXTXXXdVW2/r1q21TmAPHTpUq1evlnR+H5jqEtju9qEBAABAzfjZKrfQhlfk5+crIyNDp06dUkREhAYPHuxyXey6slgs2rp1q44fPy6TyaSBAwfW+HXFiooK7d27Vz/99JPOnj2rDh06qE+fPoqMjKxVH44cOaLvvvtOFotFUVFRuuaaa6qduV2dQYMGSZK+/fbbWl0HAAB8X9VEn1FCVWIT3/Phhx/q7bff1ooVK+Tv73rverPZrDlz5uiNN96odfsPPvigdu3apd/97ncul/HLy8vT5MmTFRAQoL/+9a8ul/K7EP/OYDRGfB4DALzHKLEJM7C9rE2bNi7XiPYUk8mkESNG1Olaf39/9e3bV3379q1XH7p06aIuXbrUqw0AAACgPt5//32Vlpbqq6++UkhIiFN5UVGR3nvvPZWXl7u8vqH3oQEAAEDNkMAGAAAA0KSYzWb99NNPkqTHHnvMbd1LL73U6Zw39qEBAABAzZDABgAAANCknDhxol7Xe2MfGgAAANQMa2DDkIyyxg4AAPA+I665SmwCd+qzD01V/DuD0RjxeQwA8B6jxCbMwAYAAACAeqjPPjQAAABwz/V23AAAAAAAAAAANDJmYANAI1i4cKHLYwAAAAAAAPwPCWwAaATJycn2YxLYAAAAAAAArrGECAAAAAAAAADAkEhgAwAAAAAAB2az2e3XAAB4CwlsAAAAAABgl5KSoujoaIdz0dHRSklJaaQeAQCaM9bABgAAAAAAks4nrxcsWOB03mKx2M8nJiZ6u1sAgGaMBDYAAAC8ruoGtmxmCwDGYDablZqa6rZOamqq4uPjFRYW5qVeAQCaOxLYAAAA8Lrk5GT7MQlsADCG9PR0WSwWt3UsFovS09N19913e6lXAIDmjjWwAQAAAACAcnJyPFoPAABPIIENAAAAAAAUFRXl0XoAAHgCCWwAAAAAAKC4uDiZTCa3dUwmk+Li4rzUIwAASGADAAAAAABJ4eHhSkhIcFsnISGBDRwBAF7FJo4AAAAAAECSlJiYKElKTU112NDRZDIpISHBXg4AgLeQwAYAAAAAAHaJiYmKj49XeHi4/Vx2djYzrwEAjYIENgAAAAwlKSmpsbsAAM3ehclqktcAgMZCAhsAAACGsnDhwsbuAgAAAACDYBNHAAAAAAAAAIAhkcAGAAAAAAAAABgSCWwAAAAAAAAAgCGRwAYAAAAAAAAAGBIJbAAAAAAAAACAIZHABgAAAAAAAAAYEglsAAAAeJXZbHb7NQAAAABUIoENAAAAr0lJSVF0dLTDuejoaKWkpDRSjwAAAAAYWWBjdwAAAADNQ0pKihYsWOB03mKx2M8nJiZ6u1sAAAAADIwZ2AAAAGhwZrNZqampbuukpqaqsLDQSz0CAAAA4AuYgQ0ABrVw4UKXxwDgi9LT02WxWNzWsVgsSk9P19133+2lXgEAAAAwOhLYAGBQycnJ9mMS2AB8XU5OjkfrAQAAAGgeWEIEAAAADS4qKsqj9QAAAAA0DySwAcDLzGaz268BoCmKi4uTyWRyW8dkMikuLs5LPQIAAADgC0hgA4AXpaSkKDo62uFcdHS0UlJSGqlHAOAd4eHhSkhIcFsnISFBYWFhXuoRAAAAAF/AGtgA4CUpKSlasGCB03mLxWI/n5iY6O1uAYDXVD7jUlNTHTZ0NJlMSkhI4BkIAAAAwAkzsAHAC8xms1JTU93WSU1NVWFhoZd6BACNIzExUdnZ2Q7nsrOzSV4DAAAAcIkENgB4QXp6usNsQ1csFovS09O91CMAaDwXLhPCsiEAAAAAqkMCGwC8ICcnx6P1AAAAAAAAmgMS2ADgBVFRUR6tBwAAAAAA0ByQwAYAL4iLi5PJZHJbx2QyKS4uTtL5NbOruvBrAAAAAACA5oAENgB4QXh4uBISEtzWSUhIUFhYmFJSUhQdHe1QFh0drZSUlIbsIgAAAAAAgOEENnYHAKC5SExMlCSlpqY6bOhoMpmUkJCgxMREpaSkaMGCBU7XWiwW+/nKdgAAAAAAAJo6ZmADgBclJiYqOzvb4Vx2drYSExNlNpuVmprq9vrU1FQVFhY2ZBcBAAAAAAAMgwQ2AHhZWFiYy6/T09MdZma7YrFYlJ6e3mB9AwAAAAAAMBIS2ABgEDk5OR6tBwAAAAAA4OtIYAOAQURFRXm0HgAAAAAAgK8jgQ0ABhEXFyeTyeS2jslkUlxcnJd6BAAAAAAA0LhIYAOAQYSHhyshIcFtnYSEBKc1tAEAAAAAAJqqwMbuAADgfxITEyVJqampDhs6mkwmJSQk2MsBAAAAAACaAxLYAGAwiYmJio+PV3h4uP1cdnY2M68BAAAAAECzwxIiAGBAFyarSV4DAGA8ubm5bsvz8/NVXl7upd4AAAA0TczABgAAgNclJSU1dhfQDJSUlGjt2rU6ePCgioqKlJWVJT8/P/Xv319jx45V796969z2qVOnNG3aNA0ePFg33XSTYmJiFBoaqjNnzujo0aPKyMjQgQMH9OSTT/KLaAAAgHoggQ0AAACvW7hwYWN3AU1cTk6OXn75ZU2bNk2TJ0+2n//vf/+rp556SuvXr9f48eN13333yc/Pr9btl5eXy2azKSMjQxkZGU7lAwcO1KJFi0heAwAA1BMJbAAAAABNSkVFhZYsWaKFCxc6JZCvv/56+fv7KykpSe+++67CwsJ0xx131Ok+rVu31oABA7Rjxw5ZLBZFRESoR48eio2N1fXXX1+nxDgAAAAckcAGAAAA0KRs27ZNERER1c5+Hjp0qKKiopSTk6P09HRNnDhRJpOp1vcJDQ3V//3f/9W3uwAAAHCDTRwBAAAANCmHDx/WgQMHlJ2dXW2d7t27S5LOnDmjzMxMb3UNAAAAtUQCGwAAAECTYjKZlJ2drenTp+uBBx7QyZMnneoEBv7vZdT8/Hxvdg8AAAC1wBIiAAAAAJqUa6+9ViEhITp79qy+//57ffjhh7rrrrsc6pw6dcp+HB0dXed7FRYWas2aNcrKytKpU6d04sQJdejQQTfccIPGjRunFi1a1LltAAAAkMAGAAAA0MS0a9dOycnJ+tvf/qbi4mJdc801DuWFhYX68ccfJUl9+/ZVz54963SfwsJCPfnkk7rvvvvUqVMnSVJJSYmWL1+ul156SRs2bNATTzyh9u3b17jNzMxMjRo1qk79kaSNGzfW+VoAAND01CeuyMzMVExMjAd7UzcksOERNptNp0+fVrt27aqtk5ubqw4dOnixV4BvS0pKauwuAADgswYOHKiBAwe6LFu9erXKysoUGRmphISEOt+jpKREs2bNsievJally5Z68MEHtW/fPh06dEjJycl67rnn5O/P6o0AAAB1QQLbi8xms1auXKmdO3cqODhYpaWl6tmzp6ZOnVqrWRnVKSkpUXp6urZs2aLAwECVlpaqS5cuuv3229W1a9dqr/voo4+0dOnSGt1j3LhxmjNnjtP5zz77TE899ZR++ctf6rrrrlPXrl0VEBCgwsJC/fzzz/rvf/+r0NBQPfLII3X+fEBzs3DhwsbuAgAATc62bdu0atUqde3aVampqXWOw4OCghQdHe0yzg4ICNDo0aP1+uuva//+/dq0aZNuvPHGGrUbExPDLGoAAOAx9YkrBg0a5MGe1B0JbC/JysrSQw89pCFDhmjZsmUKDg6W1WrV22+/rVmzZmnRokXq06dPndsvKCjQvHnz1KlTJy1evFgmk0mStGHDBs2ePVsJCQkaOnRovT9HeHi4y/NlZWUqLS3Vhg0btGHDBocyPz8/TZgwQbNmzar3/QEAgLGZzWbdcccdKi4uVqtWrfTWW29VGz8A3vb1119r0aJFGjFihObMmaOWLVvWua1LLrlEy5cvr7a8V69e9uP//ve/NU5gAwAAwBEJbC8oLCxUYmKiWrVqpdmzZ9tfHwwICNC0adO0e/duJSUl6bnnnlNUVFSt2y8tLVVSUpLy8/P1/PPPKyQkxF42duxY7dmzR0888YSWLl3qEEhXFRISoj/+8Y+64oorFBAQ4FReXl6u5ORkdevWrdp+dO/eXa1atdKPP/6o8vJytWvXTv369dP48eN1xRVX1PpzAQAA35KSkqLU1FRZLBb7uejoaCUkJCgxMbERewajSEtLa9DZxe7aXrVqlf7+979r9uzZuvnmmz1yv9atW1db1qZNG/vx0aNHPXI/AACA5ogEthe88sorys7O1u9+9zuXa9+NGTNGqampeuaZZ5SWllbr9t955x3t3btX48ePd0heV4qNjdWGDRuUlpamV1991WWCesKECRo7dqzb+4SEhKh3797Vlvfo0UMPP/xwrfsPAAB8X0pKihYsWOB03mKx2M+TxMasWbN05513evWe5eXlWrZsmfbv36/ly5e7XVqvtvbv36/CwkINHjxYfn5+DmVV4/7S0lKP3RMAAKC5IYHdwLKzs/XJJ59Ikq666iqXdSrP79ixQ7t27aq2nisWi0Vr1qyRJPXr189lnSuuuEIhISE6evSoPv30U5e7j/bv39/tfU6cOKHg4GBdeumlNe4bAABoHsxms1JTU93WSU1NVXx8vMLCwrzUKxhReHi4V5eUOXv2rJKSkhQZGannn39eLVq0cKrzySefaOTIkbVu+6233tKbb74pSZoyZYpmzJjhUF5QUGA/9sR+NwAAAM0VW2E3sA0bNshqtcrPz09dunRxWaddu3b21w/XrVtXq/a/+OILFRcXS1K1y3sEBASoc+fOkqQPPvjAqbxDhw6KiYlxe58VK1Zo3LhxteobAABoHtLT0x2WDXHFYrEoPT3dSz0CJJvNppSUFPXp00cPPPCAy+S1JG3durXWbVfuZVPp4MGDTnXy8vLsxz179qz1PQAAAHAeM7Ab2LZt2ySdn23ibpOYqKgoFRUVafv27bLZbE6vIFYnIyPDfhwZGVltvcjISGVmZmrfvn32TZUqDRgwwO099u/fr0OHDmnevHkX7c/hw4f13nvvKS8vT7m5uTp9+rQuv/xyjR07Vtdff32NPxcAAPAdOTk5Hq0HeML69et16NAhJScnV1vHbDbrxx9/rHXbhYWFKisrkyRFRERo2rRpTnV27twp6fyG5rGxsbW+BwAAAM4jgd2ASkpK9NNPP0k6H9i6U1leUFCgI0eO1Hhtvu+//16S1LJly2pnlVRtv6KiQnv37tW1115bo/bPnTunpUuX6uGHH3a5fndVu3fvVkVFheLj42UymSRJx48f12OPPabHHntMw4cP1yOPPKLAwJr9s8vMzHS53ElNNeQGQQAA4H9qugl1XTar9qT6xBWZmZkXfWMNxvL++++rtLRUX331lct9YoqKivTee++pvLzc5fVPP/20PvvsM11++eV6/PHHHSaAREREKCIiQoWFhXrhhRfUoUMHh2srKir01VdfSZImTpyojh07evCTAQAANC8ksBtQXl6ebDabJLmdfS1JoaGh9uPTp0/XKIFdXl4us9lcp/Zr6uWXX9awYcNq9B+2iooKzZ07V8HBwfZzkZGR+r//+z/dfffd2rRpk1577TX9/ve/r/H9AQCA8cXFxWnOnDlulxExmUyKi4vzYq/QnJnNZvtEkscee8xtXVd7vBQWFmr9+vWSzk8Y2bFjh4YPH24v9/Pz0+23366//OUvslqtTtf/61//Ul5enq677jrNnDmzPh8FAACg2SOB3YAqk8uSLjrruGp51esu1n5lgrwh2s/IyNB3332nl19++aJ1g4KCNHjwYIfkdaWoqCgNGDBAGRkZeu+99zRp0qQabWQTExPDLGoAAHxAeHi4EhIStGDBgmrrJCQkNPoGjvWJKwYNGuTBnqChnThxol7Xh4WFKTY21j4D29WSexMnTlRZWZkeffRRxcbGKjo6WlarVZs3b9bmzZt1zz33aPLkySyhBwAAUE8ksBtQ1SU9Kioq3NatWu5uKZCqqr4K6en2S0pKtGzZMt11110KCgq6aP0bbrhBQ4YMqbb8iiuuUEZGhsrLy7VlyxY2hAQAoIlJTEyUJKWmpjrMxDaZTEpISLCXA97Qs2fPek+EmDt3rubOneu2zuTJkzVu3Dht27ZNR44cUUBAgK6//nrFx8e7nNgBAACA2mu2Cey0tLQGnd27ceNGtWnTxv515SYv1SktLbUfV73OHZPJpKCgIJWVlXm8/RUrVig3N1e/+MUvatSXwMBAt7PA27Ztaz8+evRojdoEAAC+JTExUfHx8QoPD7efy87ObvSZ10BDMplMuvHGGxu7GwAAAE1Ws01gz5o1S3feeWeD3iMiIkKtW7dWUVGRCgsL3datLPf391d0dHSN79G5c2cdPHhQRUVFstls1b6iWPX+nTp1ctvmiRMntHbtWoWEhDj8B9Qdq9Wq7du3Kzw8XD179nQqr9qvqsl0oLlKSkpq7C4AQIMICwtzeMaRvAYA30XMCgAwgmabwA4PD69xcrau/P391b9/f23atEn5+fluE8yVGyvGxMTU6j96gwYN0sGDB1VRUaH8/HyHmc6u2m/Tpo169Ojhts133nlHZWVlNV6vz2azaf78+fruu+8kSU8++aTTOpEFBQX245qsfw00dQsXLmzsLgBAg+EZBwBNA89zAIAR+Dd2B5q6YcOGSZLOnTunkydPuqxTVlam48ePS5KGDh1aq/ar1j927Fi19bKysiRJQ4YMcZuYLisr0+effy7p/Ezp4uLii/Zh9+7d9uS1JP38889OdfLy8uzHMTExF20TAAAAAAAAAEhgN7CbbrpJHTt2lCTt37/fZZ3MzExZrVa1bt1a48ePr1X7ffv21VVXXSVJ2rdvn8s6eXl5ys3NVUBAgCZPnuy2vZ9++klFRUX2rw8fPnzRPuTn59uPe/bsqVGjRjnV2blzpyQpMjLSaXY2AAAAAAAAALhCAruBBQQEaObMmZKkb775xmWdyvPTpk2TyWRyKn/66ad1yy236IEHHnA5I3rmzJkKDAzUli1b3Lb/m9/8xp5Mr05ubq7D10eOHHFbX5K6dOkiSerWrZuWLVvmtDTL4cOH7e384Q9/kL8//+wAAAAAAAAAXByZRC8YPny4pk6dqs8//9xpeY2TJ09q3bp1GjNmjCZOnOh0bWFhodavX6+zZ8/q+++/144dO5zq9O7dW/fff7927typ7du3O5SdOXNGq1at0sCBAzVr1qyL9jU0NNTh66qzq6tz2WWXaejQoSotLXVKTldUVOj111+Xn5+fZsyYUeslUgAAAAAAAAA0X812E0dvmz59ulq0aKEFCxZo/Pjx6tatm3JycrRmzRqNHj3aPkv7QmFhYYqNjdVnn32myy+/XAMGDHBZLzY2Vn5+flq6dKnGjBmjXr16KT8/X2vWrFGfPn0UHx+vgICAi/bz6quv1qWXXqoTJ05IkqKiomr0+ebPn6/nnntOf/7znzVq1Ci1bdtWRUVF+vDDD5WXl6fk5GQNGTKkRm0BAAAAAAAAgCT52Ww2W2N3ojnJz89XRkaGTp06pYiICA0ePFgdOnTwWPsWi0Vbt27V8ePHZTKZNHDgQHXq1KnWffziiy/UunVrjRgxolZLfuTm5mrHjh3Kz89Xq1at1K1bN1155ZVuN450pXKd7G+//bZW1wEAADQEYhN4A//OAACAkRglNmEGtpe1adNGo0ePbrD2TSaTRowYUa822rRpowkTJtTp2g4dOmjMmDH1uj8AAAAAAAAASKyBDQAAAAAAAAAwKBLYAAAAAAAAAABDIoENAAAAAAAAADAkEtgAAAAAAAAAAEMigY0mZ9SoURo1alRjdwMXwTj5BsbJNzBOvoOx8g2ME3wN/2Z9B2PlGxgn38A4+Q7GyjcYeZxIYAMAAAAAAAAADIkENgAAAAAAAADAkEhgAwAAAAAAAAAMiQQ2AAAAAAAAAMCQSGADAAAAAAAAAAyJBDYAAAAAAAAAwJBIYAMAAAAAAAAADIkENgAAAAAAAADAkEhgAwAAAAAAAAAMiQQ2AAAAAAAAAMCQSGADAAAAAAAAAAyJBDYAAAAAAAAAwJBIYAMAAAAAAAAADIkENgAAAAAAAADAkPxsNputsTsBXKht27Y6e/asevfuXetrMzMzJUkxMTGe7hY8iHHyDYyTb2CcfAdj5RtcjdMPP/ygkJAQ5eXlNVa30AzUNQbm2eI7GCvfwDj5BsbJdzBWvsHIMTAJbBhSt27dVFhYqMsuu6yxuwIAAKCff/5ZYWFhOnToUGN3BU0YMTAAADASo8TAJLABAAAAAAAAAIbEGtgAAAAAAAAAAEMigQ0AAAAAAAAAMCQS2AAAAAAAAAAAQyKBDQAAAAAAAAAwJBLYAAAAAAAAAABDIoENAAAAAAAAADAkEtgAAAAAAAAAAEMKbOwOAJ6yZ88erVy5Uvn5+QoICJAkjRgxQrfccov9azSsn3/+WevXr9fp06d16tQp5ebmql27drruuut0yy23KCwszO31ZrNZK1eu1M6dOxUcHKzS0lL17NlTU6dOVfv27b30KZons9msmTNnaubMmRozZsxF6zJO3nPw4EGtXr1a2dnZCgwMlNVqVf/+/TVlyhQFBwdXe11ubq7+8Y9/KDMz0z5OAwYM0OTJky/6vYjaKSgo0Nq1a5WRkSFJCgoKksViUb9+/XTLLbeoa9eubq8/dOiQ/vnPfyorK0vBwcGyWq267rrrFBcXpxYtWnjjIzRJy5cvV3h4uO68884a1ffEs23jxo1av369rFarKioqFBoaqri4OF1zzTX1+SiAW8TAjYv417cRAxsXMbDxEQMbU5ONgW1AE/Dvf//bNm7cONvWrVvt506dOmWLj4+3zZs3z3bu3LlG7F3z8P7779v+8pe/2IqLi+3niouLbcuXL7eNHDnSNnHiRFtGRka11x89etT2//7f/7M9++yz9vEqLy+3vfXWW7Zbb73VtmfPngb/DM3ZokWLbCNHjrRt2LDBbT3GybvWrl1r++1vf2vbsmWL/ZzVarUtWbLENm/ePFtFRYXL63bv3m2bMGGC7R//+IetvLzcZrPZbOfOnbM988wztttvv9127Ngxr/S/Ofjqq69st956q+3ll192eP6ZzWbbs88+a4uNjbWtW7eu2uu//PJL27hx42wfffSR/VxRUZHt0Ucftd177722goKCBu1/U3XixAnbhAkTbCtWrKhR/fo+26xWq23RokW26dOn27Kysuzn9+3bZ5s0aZLttddeq/uHAdwgBm5cxL++jxjYmIiBjY8Y2JiacgxMAhs+b9OmTbbRo0fb3n77baey48eP20aPHm1LTk6u9occ6m/Pnj22JUuWVFuelJRkGzlypG3cuHG2Q4cOOZWbzWbbnXfeaZs5c6bNarU6lT/88MO2uLg4W3Z2tkf7jfM2bdpkGzly5EWDd8bJuz7++GPbqFGjbNu3b3c4X15ebnvsscdsI0eOtP3www9O1x07dsw2adIkW0JCglNZeXm5bcaMGbbp06fbzGZzg/W9udi1a5dt7Nixtr/97W/V1lm8eLFt9OjRLoO/PXv22G6++WbbsmXLnMrOnDljmzBhgm3OnDkkoGqhoqLCtnv3btuMGTNsI0eOrFHw7oln24svvmgbOXKkbe/evU5lGzdutI0cOdL273//u1afBbgYYuDGRfzr+4iBjYkY2PiIgY2nOcTArIENn2axWPT000/LZrNp7NixTuWXXnqpBgwYoE2bNuk///lPI/SweVizZo2GDx9ebfnEiRMlSWfPntXKlSudyl955RVlZ2dr9OjR8vd3fiyNGTNGBQUFeuaZZzzXaUiSCgsL9dxzz9WoLuPkPfn5+XrxxRc1aNAgDRgwwKHspZde0pdffilJOn36tNO1zzzzjMxms2JjY53KAgICNGrUKGVlZenVV19tmM43I6+88orKy8s1bty4autMnz5dFRUVWrNmjcN5q9WqtLQ0lZaWuhyrli1b6sYbb9SePXu0evVqj/e9qUlPT9eDDz6oe++9VykpKTp69GiNr63vs23Xrl1as2aNunXrpt69ezuV33jjjQoNDdXLL7+snJycmn8owA1i4MZH/OvbiIGNiRjYNxADG0dzioFJYMOnrVmzRkVFRerSpYsiIiJc1rnqqqskSX//+99ltVq92Lvm4/Dhw9qyZYvOnTvnsrxHjx724127djmUZWdn65NPPpH0v7G6UOX5HTt2OF2P+nn++efVsmXLi9ZjnLzrnXfeUXFxsW666SansoqKCklSly5d1L9/f4eynTt36rvvvpMk9evXz2XbleO0ceNGkmn1UFJSov3790uSWrduXW29Dh06KCIiQocPH3Y4/5///EfHjh1TaGiowzOyqsqx+te//iWLxeKhnjdNcXFxeuqpp/Taa6/VKoHgiWfbm2++Kan677mgoCD16tVL5eXl+sc//lHjvgHuEAM3PuJf30YMbEzEwMZHDGwszSkGJoENn7Zu3TpJcrs5wGWXXSZJOnbsmP2HGjwrNDRU77//vm677TaXv9EODPzffrH5+fkOZRs2bJDVapWfn5+6dOnisv127drZfzhWjjnq7+uvv9Z3332nmTNnXrQu4+Q9VqtVn376qSTpyiuvdCr/4x//qNdff10vvfSSTCaTQ9mHH34oSYqIiKg2oVH5TLRarVq/fr0He968VFRUyGazSZL27dvntu65c+fUtm1bh3OV3yOdO3eudpO1yrEqKirSF198Ud8uw4X6PtuysrK0c+dOSVK3bt2qvU/lWH766acqKSnxQM/R3BEDNz7iX99FDGxMxMC+gRi4afDFGJgENnzWzz//bH91KDIystp6UVFR9uNt27Y1eL+aoxtvvFHS+d/Grlq1yuk32qdOnbIfR0dHO5RVjkl4eLjbWRCV47h9+3b7D0zUXWFhoZ599lndf//9NdqNm3Hynn379ik/P18BAQEOz69K/v7+6tKli4KCghzOV1RU6Ntvv5Xk/pkYGhqq8PBwSTwT68NkMtn/npctW6asrCyX9Q4ePKiSkhKNHDnSfq64uNge8Lsbq6pljFXDqO+zreq41CQWKSsrswf7QF0RAxsD8a9vIgY2LmJg30AM3DT4YgxMAhs+a/fu3fbj6n7LemHZ999/34A9ar7Gjx+vKVOmqGPHjvrFL36hdu3aOZRnZGTYj2+99Vb7cUlJiX766SdJ7sewanlBQYGOHDnimY43Yy+++KKuuuoqXX/99Retyzh5V+UreeHh4QoICNDhw4eVlpamuXPn6uGHH9YLL7zg8rXHo0ePymw2S6r5OP3000/MBq2HSZMmSTo/A+Hee+/V0qVLHX7OVFRU6LXXXtPYsWM1atQo+/k9e/bYX4N1N1YtW7ZUSEiIJMefefAMTzzb6hKLMJaoL2JgYyD+9U3EwMZFDOw7iIF9m6/GwIEXrwIYU9WNG9z9xig0NNTlNfCcgIAAzZgxQzNmzHAqO3v2rP79739LksaOHaubb77ZXpaXl2f/Ld7F1qC7cBzdvTIL9zZv3qxvv/1Wr7/+eo3qM07edejQIUnn1wzbvHmzXn31Vc2bN099+vSRxWJRWlqa7r77bv35z392WB+w6vOtMuCrTuU42mw25efn12gNSDibMGGCjh49qrVr16q8vFwfffSRPvroI3Xs2FGjRo3SgQMHNHToUKcNbmr686uy/OzZszKbzbJardW+aona88SzjVgEjYF/d8ZA/Ot7iIGNjRjYdxAD+zZfjYGZgQ2fVflbVslxjbkLVS0rKChoyC7hAjabTc8995yys7N16623au7cuQ7lNR3DC8urXofaKSoq0rJly2r82qTEOHlbYWGhpPPrZS5ZskSPPPKI+vTpI+n8K3sPPvig/Pz8lJaWpmPHjtmvq/p8u/DVygtVLee5WD/x8fFavHixw1qM2dnZWrFihb7++mtt2bJFBw4ccLimLt9TNpuN7ykP88SzrS6xCOOI+iIGNjbiX2MiBjY+YmDfQgzsu3w1BiaBDZ/VokUL+3HlayiuVC272G9k4Tnl5eVKS0vTpk2bNH/+fN13333y8/NzqFPTMbywvOp1qJ0XX3xR/fr10/Dhw2t8DePkXZWvM5aWlqpv3766/PLLHcrDwsLUr18/lZWVaeXKlfbzVZ9vFxsnq9Xq8jrU3ldffaUXXnhBDz30kN566y1Nnz5dnTt3lnR+HL755hvdd999+uijj+zX8D1lDJ4Yh7rEIowj6osY2LiIf42LGNj4iIF9CzGw7/LVGJgENnxWmzZt7MdlZWXV1istLbUfX2x9H3hGYWGh5s+fr4MHD+qFF17Qr371K5f1ajqGkuM4Vr0ONffNN99o27Ztmj17dq2uY5y8q+orWIMHD3ZZp1evXpLOj2mlqs+3quPgStVxZJzqbu3atVq4cKHuueceDRs2TJGRkZo6dar++te/avny5fYNvioqKvTMM8/YN7mpy/dUcHCwwwwX1J8nnm11iUX4nkN9EQMbE/GvcRED+wZiYN9BDOzbfDUGJoENn9WpUyf7ceXrRq5ULav8jSAazpEjRxQfH68uXbpo+fLl6tKlS7V1IyIi1Lp1a0nux7Bqub+/v9NO7ri44uJiLVu2TLNnz671f2IZJ++q+lrrpZde6rJO5Q7qBQUF9tkqnTp1kr//+R/rRUVFbu9ROU6tW7e2t4XaKSoq0quvvqpLL71UQ4YMcSrv1auXFixYoPnz50s6P+Pn888/l+T4s8jd91RFRYUsFoskx5958AxPPNuIRdAY+HdnPMS/xkUM7DuIgX0DMbDv89UYmE0c4bMGDBiggIAAWa1Wt4vBVy0bNGiQN7rWbP3www9KTU3Vvffe6/L1vIqKCn3++ef65S9/Ken8Q7B///7atGmT8vPzZbPZnF6zrFQ5jjExMTVetw7/s2PHDgUFBWnt2rVau3atQ1lxcbH9+J133tHHH38sSerRo4fuu+8+xsnLqgYDwcHBLutUPV8ZsIeFhenyyy/Xjz/+6PaZWFFRoby8PEnnn6OV16N29uzZo7Nnz6pfv37Vfj9I0q9+9St99NFH2rFjh/3vvUePHoqIiFBBQYHbscrLy7O/dlfdTCTUnSeebYMHD9a7777rUMfd9RKxCOqPGNhYiH+NjRjYdxAD+wZiYN/nqzEwCWz4LJPJpKuvvlrbt2932MThQkePHpV0/pv0uuuu81b3mp2srCylpKQoOTlZPXr0cFknOzvb/vpQpWHDhmnTpk06d+6cTp48qQ4dOjhdV1ZWpuPHj0uShg4d6vnONwPDhw+vds2/nTt3at68eZKkyZMna8yYMU51GCfvufrqq+3H1f02Oz8/X5LUrl07h7XEhg4dqh9//FHHjx+vdrfuEydO2F/zGjZsmCe73qycO3dOkuyzF9zp3r27duzYoUsuuUSS5OfnpyFDhmj9+vVuf35VfV7yPdUw6vtsu/rqq2UymWSxWGo0ltHR0fbd24G6IgY2DuJf4yMG9h3EwL6BGLhp8MUYmF85wadNmTJFkpSZmVntwvH79u2TJI0aNcrlNyU8Y+nSpRo9enS1wbskZWRkOJ276aab1LFjR0nS/v37XV6XmZkpq9Wq1q1ba/z48Z7pMGqFcfKemJgYRUZGSpJyc3Nd1qkMBC58bW/8+PFq1aqVysrKnHb9rlT5TIyOjravT4fai4qKkvS//0i5UxnoV51BMnnyZAUEBOj48ePV7shdOVb9+/fXlVdeWd8uw4X6PttCQkI0ceJEt9dL/xvLqVOneqLbADGwQRD/Nn2MlfcQA/sGYuCmwRdjYBLY8GlXX321rrvuOlksFu3evdup3Gq1auvWrTKZTJo2bVoj9LB5OHjwoPbs2aNDhw5p69atysjIcPizZcsWrV69WitWrHC6NiAgQDNnzpTkuBlHVZXnp02bxgYOjYRx8h5/f3/7D/itW7c6lZeXl2vr1q0KCgrSpEmTHMpatWplv3bLli0u268cp5kzZ7qcnYKaiYmJUffu3bVnzx63Afy5c+f09ddfq3///rriiivs56Ojo3XLLbfIZrNVO1abN29WYGCg7rnnHo/3H+d54tk2adIktW/fXrt27bKvx1lVVlaWjh49qpiYGPsSAkB9EQM3PuLf5oGx8h5iYN9ADNw0+GIMTAIbPm/+/Pnq1q2b3nrrLVmtVoeyDz/8UPn5+Xr00Uftv82F53377beSpE2bNikxMVGPPPKIw58FCxbolVdesW/EcKHhw4dr6tSp+vzzz/Xzzz87lJ08eVLr1q3TmDFj7L/hg2dV/mZccr9zN+PkPaNGjdLQoUO1e/duHT582KFsw4YNKi4u1r333utyk6C4uDiNHDlSa9eu1alTpxzKDhw4oE2bNunOO+/U9ddf36Cfoanz8/PTgw8+qKCgIKWmprp8vlmtVj377LMKDg5WQkKCU/msWbPUv39/rVy50inoy8jI0N69e/XAAw+oV69eDfY5mqKaPtMq1ffZZjKZ9Pjjj8tms2nlypVO5StWrFD79u312GOP8R9meBQxcOMi/vV9xMDGQwxsfMTAxtXUY2A/m81mq3crQCPLy8tTWlqaSktL9etf/1qhoaHasWOHNm/erPvvv1/XXHNNY3exSXv++ef13nvv1ajuHXfcoTvvvNNl2dtvv60PPvhA48ePV7du3ZSTk6M1a9Zo6NCh/Kbcw1566SVlZmZKOv+b0cqNNVq3bq3LLrtMkhQZGamHHnrI6VrGyTtKS0v19NNPa+fOnZoyZYoiIyO1Z88erVu3TjNmzNDNN99c7bXl5eV65ZVXtGXLFk2aNEmRkZH6+eeftXbtWo0fP16//e1vvfhJmrbDhw/r+eef15EjRxQbG6tevXrJz89P2dnZ+vTTT9W1a1f9/ve/V6tWrVxef+bMGT377LM6dOiQbr31VrVp00Y//PCDPv74Y82YMUOjRo3y8ifyTenp6dq8ebMk6dSpU8rOzpZ0/vXGmJgY+fn5KTAwUIsXL662jfo+2/bt26fFixfryiuv1JAhQ1ReXq5PPvlEBQUFSkhIsL+mCXgSMXDjIf71TcTAxkcM7BuIgY2hOcXAJLDRpPz444/6/vvvVVpaqi5dumjw4MHV7mAMY8rPz1dGRoZOnTqliIgIDR48mHUbDYhx8p4jR47ou+++05kzZ9SxY0cNHDiw2kDwQrm5udq2bZsKCgrUvn17XXPNNYqIiGjYDjdTOTk52r17t/Ly8hQSEqJ27dqpf//+NR6rynG2WCyKiorSNddcw2vIjaC+zzar1art27fr4MGD8vf3V+/evdW3b98G7DFwHjGwbyOu8h2MlfcQA/sGYuCmwRdiYBLYAAAAAAAAAABDYg1sAAAAAAAAAIAhkcAGAAAAAAAAABgSCWwAAAAAAAAAgCGRwAYAAAAAAAAAGBIJbAAAAAAAAACAIZHABgAAAAAAAAAYEglsAAAAAAAAAIAhkcAGAAAAAAAAABgSCWwAAAAAAAAAgCGRwAYAAAAAAAAAGBIJbAAAAAAAAACAIZHABgAAAAAAAAAYEglsAAAAAAAAAIAhkcAGAAAAAAAAABgSCWwAAAAAAAAAgCGRwAYAAAAAAAAAGBIJbAAAAAAAAACAIZHABgAAAAAAAAAYEglsAAAAAAAAAIAhkcAGAAAAAAAAABgSCWwAAAAAAAAAgCGRwAYA+LTk5GQFBQXJz8/P6U9QUJDeeOMNSdL69evVokULpzqBgYGaMmVK434IAAAAoBaIgQE0J342m83W2J0AAKC+1qxZo9tuu00VFRWSpBtuuEGffPKJgoKCHOrNnTtXzzzzjDp37qwVK1bohhtuUEBAQGN0GQAAAKgXYmAAzQEzsAEATcLEiRM1c+ZM+9f5+flOgbskRUZGKiQkRB9++KFGjBhB4A4AAACfRQwMoDlgBjYAoMk4ffq0evbsqby8PEnSO++8o9/+9rf28rNnz6pfv3763e9+p4ceeqixugkAAAB4DDEwgKaOGdgAgCbjkksuUUpKiv3refPm6cyZM/avk5OTFRERoblz5zZG9wAAAACPIwYG0NQxAxsA0KRUVFTo2muv1bfffitJSkxM1KJFi/Tdd99p2LBh+uabb9SvX79G7iUAAADgOcTAAJoyEtgAgCbnm2++0dChQ2Wz2dSiRQvt2rVLU6ZM0a9//WslJyc3dvcAAAAAjyMGBtBUkcAGADRJ99xzj/76179Kkjp06KD27dtr+/btCg4ObuSeAQAAAA2DGBhAU8Qa2ACAJunJJ59URESEJCk3N1e///3vCdwBAADQpBEDA2iKSGADAJqk9u3b6/HHH7d/vXTpUpWUlDRijwAAAICGRQwMoCkigQ0AaLJCQ0Ptx4cPH1Zqamoj9gYAAABoeMTAAJoa1sAGADRJJ06cUJ8+fXT//fcrOTnZvpnN3r171b1798buHgAAAOBxxMAAmiJmYAMAmqQ5c+Zo7NixSkpK0t133y1JOnfunObMmdPIPQMAAAAaBjEwgKaIGdgAgCZn3bp1mj59uvbu3av27dvr1KlT6tmzp/Lz8yVJ77//vsaNG9fIvQQAAAA8hxgYQFPFDGwAQJNSXFysP/zhD3r66afVvn17SVK7du20aNEie50//elPOnfuXGN1EQAAAPAoYmAATRkJbABAk/LII4+oe/fuuuOOOxzO//73v9fAgQMlSQcOHNCSJUsao3sAAACAxxEDA2jKWEIEANBkbNy4UWPGjNGqVasUFxfnVP7ll1/qxhtvlCS1bNlSP/zwg7p27ertbgIAAAAeQwwMoKljBjYAwOeVl5frgw8+0G233SabzaYjR464rNe9e3f5+flJkkpKSnTHHXfo5MmT3uwqAAAA4BHEwACaC2ZgAwB82pQpU7Rq1SpZrVaH823bttXp06clSSdOnFDXrl1VWloqVz/2WrRooS+//FLXXnutV/oMAAAA1AcxMIDmhAQ2AAAAAAAAAMCQWEIEAAAAAAAAAGBIJLABAAAAAAAAAIZEAhsAAAAAAAAAYEgksAEAAAAAAAAAhkQCGwAAAAAAAABgSCSwAQAAAAAAAACGRAIbAAAAAAAAAGBIJLABAAAAAAAAAIZEAhsAAAAAAAAAYEgksAEAAAAAAAAAhkQCGwAAAAAAAABgSCSwAQAAAAAAAACGRAIbAAAAAAAAAGBIJLABAAAAAAAAAIZEAhsAAAAAAAAAYEgksAEAAAAAAAAAhkQCGwAAAAAAAABgSCSwAQAAAAAAAACGRAIbAAAAAAAAAGBIJLABAAAAAAAAAIZEAhsAAAAAAAAAYEgksAEAAAAAAAAAhvT/ATuLyT1Dw/HzAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 368,
       "width": 728
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
    "\n",
    "for i in range(2):\n",
    "    axes[i].errorbar(x, y[..., i], yerr[..., i], marker=\"o\", ls=\"none\")\n",
    "    axes[i].set(xlabel=\"x\", ylabel=f\"$y_{i+1}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's implement this model in JAX so we can make use of the auto-differentiation capabilities. We'll again compute the Fisher information matrix and compare it to the maximum likelihood parameter estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "jax.config.update(\"jax_enable_x64\", True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jax.jit\n",
    "def model_f(x, pars):\n",
    "    return jnp.stack(\n",
    "        (\n",
    "            pars[0] * jnp.sin(2 * jnp.pi * pars[2] * x),\n",
    "            pars[1] * jnp.cos(2 * jnp.pi * pars[2] * x),\n",
    "        ),\n",
    "        axis=-1,\n",
    "    )\n",
    "\n",
    "\n",
    "@jax.jit\n",
    "def ln_likelihood_one_star(pars, x, y, Cinv):\n",
    "    y_model = model_f(x, pars)\n",
    "    dy = y - y_model\n",
    "    return -0.5 * dy.T @ Cinv @ dy\n",
    "\n",
    "\n",
    "ln_likelihood_helper = jax.vmap(ln_likelihood_one_star, in_axes=(None, 0, 0, 0))\n",
    "\n",
    "\n",
    "@jax.jit\n",
    "def ln_likelihood(pars, x, y, Cinv):\n",
    "    return jnp.sum(ln_likelihood_helper(pars, x, y, Cinv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try evaluating our log-likelihood function at the true parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(-11.87854157, dtype=float64)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ln_likelihood(true_pars, x, y, Cinv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([[ 4.25782094e+05, -0.00000000e+00, -7.83052999e+06],\n",
       "       [-0.00000000e+00,  6.54218186e+03,  6.83218150e+05],\n",
       "       [-7.83052999e+06,  6.83218150e+05,  6.17018283e+08]],      dtype=float64)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = -jax.hessian(ln_likelihood)(true_pars, x, y, Cinv)\n",
    "F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And once again, we can turn this into a prediction for the uncertainty on our model parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.78625210e-03, 1.34166769e-02, 4.98967477e-05])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Finv = np.linalg.inv(F)\n",
    "Fisher_param_uncertainties = np.sqrt(np.diag(Finv))\n",
    "Fisher_param_uncertainties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, because our model is nonlinear, we can't use the least-squares expressions above to analytically compute the maximum likelihood parameter estimates and precisions. Instead, we will specify a prior and use MCMC to generate samples from the posterior distribution. We can then compare the covariance matrix of the samples to the inverse of the Fisher information matrix (or just the diagonal terms, the parameter uncertainties)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jax.jit\n",
    "def ln_prior(pars):\n",
    "    lp = 0.0\n",
    "\n",
    "    # Gaussian priors on the amplitudes with mean=0, stddev=10:\n",
    "    lp += jax.scipy.stats.norm.logpdf(pars[0], 0, 10)\n",
    "    lp += jax.scipy.stats.norm.logpdf(pars[1], 0, 10)\n",
    "\n",
    "    # Uniform prior on the frequency:\n",
    "    lp += jax.scipy.stats.uniform.logpdf(pars[2], 0, 1)\n",
    "    lp = jnp.where((pars[2] > 0) & (pars[2] < 1), lp, -jnp.inf)\n",
    "\n",
    "    return lp\n",
    "\n",
    "\n",
    "@jax.jit\n",
    "def ln_posterior(pars, x, y, Cinv):\n",
    "    return ln_likelihood(pars, x, y, Cinv) + ln_prior(pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import blackjax\n",
    "\n",
    "rng_key, warmup_key, sample_key = jax.random.split(jax.random.PRNGKey(42), 3)\n",
    "init_pars = true_pars\n",
    "func = lambda pars: ln_posterior(pars, x, y, Cinv)\n",
    "\n",
    "warmup = blackjax.window_adaptation(blackjax.nuts, func)\n",
    "(state, parameters), _ = warmup.run(warmup_key, init_pars, num_steps=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inference_loop(rng_key, kernel, initial_state, num_samples):\n",
    "    @jax.jit\n",
    "    def one_step(state, rng_key):\n",
    "        state, _ = kernel(rng_key, state)\n",
    "        return state, state\n",
    "\n",
    "    keys = jax.random.split(rng_key, num_samples)\n",
    "    _, states = jax.lax.scan(one_step, initial_state, keys)\n",
    "\n",
    "    return states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = blackjax.nuts(func, **parameters).step\n",
    "states = inference_loop(sample_key, kernel, state, 10_000)\n",
    "\n",
    "mcmc_samples = states.position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's look at the parameter uncertainties as estimated by MCMC:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([1.80400079e-03, 1.35137301e-02, 5.07659010e-05], dtype=float64)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(mcmc_samples, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare that to our Fisher approach above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.78625210e-03, 1.34166769e-02, 4.98967477e-05])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fisher_param_uncertainties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are very similar, as we expect!\n",
    "\n",
    "We can also look at the full covariance matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.25474432e-06, -4.52384312e-06,  4.79426144e-08],\n",
       "       [-4.52384312e-06,  1.82639165e-04, -2.74545006e-07],\n",
       "       [ 4.79426144e-08, -2.74545006e-07,  2.57743445e-09]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.cov(mcmc_samples.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.19069658e-06, -4.78172971e-06,  4.57876380e-08],\n",
       "       [-4.78172971e-06,  1.80007218e-04, -2.60004737e-07],\n",
       "       [ 4.57876380e-08, -2.60004737e-07,  2.48968543e-09]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Finv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py312",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}