{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "870343b9-07ae-491a-b743-650b69e33f98",
   "metadata": {},
   "source": [
    "# How to: Use `jackknify` for line finding. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cd6f83b-8236-46fe-b940-93cccfc82367",
   "metadata": {},
   "source": [
    "`jackknify` can be used in a variaty of ways. Here, we show how to use the jackknifed data sets for inference of a line detection. According to JWST data, In the ALMA data we have been working with, there should be [OIII] emission coming from a galaxy which some 14 Billion light years away. Let's see if we can find it, or if it is undistinguishable from noise.\n",
    "\n",
    "The line finding is done with a code named [Source EXtractor](https://www.astromatic.net/software/sextractor/). Source extractor is integrated in a Python interface through the [`Interferopy` package](https://interferopy.readthedocs.io/en/latest/). Sadly, the combination of `interferopy` and `casatasks` restricts the usage of Python version to Python=3.8. We ran the linefinding seperately, but you can you find the output catalogs also on the [google drive](https://drive.google.com/file/d/1FlQNwy7VtAk0zcFfdW5tMY2aGKivxlyR/view?usp=sharing) and the script we used to generate the outputs in the same folder as the tutorials. \n",
    "\n",
    "Let's first load in everything. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0fde5d3-b707-4110-ad06-cd64a52eef69",
   "metadata": {},
   "source": [
    "### Load in "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2561edc8-1e1b-4f78-ba4f-fd081fe267fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from astropy.io import ascii, fits\n",
    "from astropy import constants as c\n",
    "from astropy import units as u\n",
    "from astropy.modeling import models\n",
    "from astropy.coordinates import SkyCoord\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['text.usetex'] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "db2df258",
   "metadata": {},
   "outputs": [],
   "source": [
    "def data(path, positive = True):\n",
    "    if positive:\n",
    "        return ascii.read(path+'findlcumps_clumpsP_minSNR_0_cropped.cat')\n",
    "    if not positive:\n",
    "        return ascii.read(path+'findlcumps_clumpsN_minSNR_0_cropped.cat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fcefe7d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#noise\n",
    "galaxy = 'Glass-z13'\n",
    "\n",
    "outdir = '../../output/findclumps/'\n",
    "paths  = [outdir+x+'/' for x in os.listdir(outdir) if x.startswith('.') is False]\n",
    "paths_real = [s for s in paths if 'Jack' not in s]\n",
    "paths_jack = [s for s in paths if 'Jack' in s]\n",
    "\n",
    "paths_real = [s for s in paths_real if galaxy in s][0]\n",
    "paths_jack = [s for s in paths_jack if galaxy in s]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb1571cd-906d-4405-a3de-3b45763984ed",
   "metadata": {},
   "source": [
    "## Plot sampled probabillity functions\n",
    "\n",
    "As explained in [Vio & Andreani 2016](https://arxiv.org/abs/1602.02392), the underlying distribution that sets the likelihood of false detection is the distribution of peaks of a smoothed (close to) Gaussian random field. With the jackknifed measurment sets we effectively sample this distribution. This is more complete and more acurately than, for isntance, using the distribution of negative peak values as used in [Walter+2016](https://arxiv.org/abs/1607.06768). \n",
    "\n",
    "Let's first convert the dataframes to numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "41437f16",
   "metadata": {},
   "outputs": [],
   "source": [
    "RAs_real_P   = np.array(data(paths_real)['RA'])\n",
    "Decs_real_P  = np.array(data(paths_real)['DEC'])\n",
    "data_real_P  = np.array(data(paths_real)['SNR']) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b77084bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "RAs_jack  = np.empty(0)\n",
    "Decs_jack = np.empty(0)\n",
    "data_jack = np.empty(0)\n",
    "\n",
    "for idx, path in enumerate(paths_jack):\n",
    "    RAs_jack  = np.append(RAs_jack,  np.array(data(path)['RA']))\n",
    "    Decs_jack = np.append(Decs_jack, np.array(data(path)['DEC']))\n",
    "    data_jack = np.append(data_jack, np.array(data(path)['SNR']))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3641813b-ec8a-4394-bf24-29567590a2c4",
   "metadata": {},
   "source": [
    "### Detection inference\n",
    "\n",
    "Now let's compare the noise distribution --- drawn fromt the jackknifed data sets --- with the positive peak distribution of the real initial data set. Since the astronomical signal is positive, any excess in peaks distribution in the data can be considered real."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f9f69655-26e5-4e1a-ac6f-d01f2f723b98",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGGCAYAAAAjCmXBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABcYUlEQVR4nO29f2wj6Xnn+a0iRaqlllRi93h+2D0zKjnOppNzspSE29wZWfhMJf5jFwhwYneABdZADi0a2xfcJZiI6T+MwJs/dFQ2axxyvWuxFwYmRg5RkwEC2L5zQk6MMYzL5SQxdmKMszNW9Xg6Y3e7u6nSz+bPqvujVMUf4o+qYtXLKur5AAKaRbLqJZ/m8633fZ73eThVVVUQBEEQhE/ghz0AgiAIgrACCRdBEAThK0i4CIIgCF9BwkUQBEH4ChIugiAIwleQcBEEQRC+goSLIAiC8BUkXARBEISvIOEiCIIgfAUJF0G4RCKRwOzsLDiOw/z8PBKJBCRJMp5fXl5GMpkc4ggJwp+QcBGEw0iShPn5eUiShEwmg/39fWxubqJYLCKbzQ57eH3JZrNYXl4+d1yWZcTjcczOzmJ+fp5ElxgawWEPgCBGjXg8DlEUkcvljGOxWAyxWGyIo+pPMplEOp1GJBLp+PytW7ewvLyMe/fuYWdnx5hBZjIZxiMlLjo04yIIB8lmsygUCtjc3Bz2UCyTSqWwv7/fcSYlSRJkWcbq6ioEQUAsFsPm5iay2SxkWWY/WOJCQ8JFEA6ytbWFaDQKURQtvzebzWJhYcGIibUvKyaTSSNmtrCwgHw+3/N4p/NzHHfuL51O9x2bKIrnxHhxcREAsLOzY/mzEsQgkHARhINIkmQ4dKsUi0Xcu3cPqqpic3MT8XgchUIBAJDP55HNZvHgwQOoqopUKoVIJNL1eCdWVlagqqrxt7q6ClEUsbq6amp87WKsC5bdz0sQdqEYF0E4iCzLEATB1nubBSQWi0EUReTzeUSjUciyjGKxaJxbj5dls9mOx/uRz+eRTqexu7tra6yAtrS4trZm+/MShF1oxkUQDiKKojFL0onH48ayXKdsvWbS6TTi8TgWFhZaUudjsRgikYhxDn0ZsdvxfsTjcaytrSEajVr8hI33R6NRpFIpW+8niEEg4SIIB1leXkY+n29JWMhkMlBVFWtraz3fu7CwgEwmg0Qigd3d3RZREQQBe3t72NzchCAIiMfj2NjY6Hpcj5XpfxsbG8a54vE4IpGIbdHRsyZJtIhhQcJFEA6ytrYGURQt73GSJAmFQgG5XK7nct/q6ioymQw2NzextbXV9fju7m5LPEsXTT0mZjeFPR6PY3l5mUSLGCokXAThMJlMBvfv3zeSK2RZRqFQaFn6a0dPqNAz/PS0ep1sNouNjQ3IsgxZlpHL5SCKYtfjndA3EKdSKVtLhPF4HEtLS7hx44ZxPUqFJ4aCShCE4+zv76urq6tqNBpVAaiiKKorKytqLpczXhOLxdS1tTXj8dramgpAFQRBXV1dVWOxmLq5uamqqqru7u6qsVhMFQRBBaDGYjF1f3+/6/FO6Odv/9PHsLm5ee45URRVVVXVvb29ju8F0PKZCIIFnKqqKmOtJAiCIAjbjGQ6fDqdNrK7VlZWbG0GJQiCILzJyMW4JEnC3t4eYrEY1tbWqBAoQRDEiOF54SoUClhYWDh3XJIkbGxstASnAS1ran5+vuV1BEEQxOjg6aXCbDbbcUMnoGU46bv+JUnCrVu3kMlkzlUuoKwngiCI0cLTwrWystLxePssSi+NA2gbNUmsCIIgRhfPLxV2Ip/PnyskGolEUCgUsLi4iGfPnhnH7Za0IQiCILyJp2dc3eg2oyoWi4jFYtjZ2UE+n4ckSbh3717X85TLZZTLZeOxoigoFou4cuUKOI5zetgEQRBEF1RVxdHREV555RXwfO85lS+Fqxu6oJlt07C+vo4vfvGLLo6IIAiCsMLDhw/xsY99rOdrfClcgiCgWCy2HGtu7WCWO3fu4Hd+53eMxwcHB3j11Vfx8OFDTE9POzFUgiAIwgSHh4e4du0apqam+r7Wl8Kltw1vx2pDu3A4jHA4fO749PQ0CRdBEMQQMBOm8U1yRnNcq70Sht51lhraEQRBjD6ennHl83nkcjkAWjxqaWnJSJHPZDJIJpNYWlrC9va27TYNAHD37l3cvXsX9XrdkXETBEEQ7kFFdps4PDzEzMwMDg4OaKmQIAiCIVb8r2+WCgmCIAgCIOEiCIIgfAYJF0EQBOErSLigJWdcv34dS0tLwx4KQRAE0QdKzmiCkjMIgiCGAyVnEARBECMLCRdxjlKthJtfv4mbX7+JUq007OEQBEG0QMJFEMRAbGxsgOO4jp3K0+l0S0dyVmNxguXlZWxsbAx8nkKhwLzbRKFQQCKRwPz8PDiOw+zsLDiOM+zkxOcaJp6unMEKqpxBeIFnx2Ucl2tDHcPlcBBXLp+v39kPQRBQKBRQKBSG2gOvWx3Ti0QymcT8/DxSqRTW19eRzWaxt7cHQCuPl0qlkEwmkcvljMpEfoOEC8Dt27dx+/ZtIzh40Sg8LuCgfGA8rtQrOKocAQDefvg2QoHQuffMhGcQfZGadDrFs+My/pc/+zvsn1YHOk+pWkepqmB8jMf4WMDya2YnxvC//8Y/tyxeoihCFEXDIQ6LaDR6oZvHptNpQ5wAnLuREEXREPZ0Og1Zln1Z45WE64JTeFzA5775ua7Pv/HtN7o+9+Zn3yTxcojjcg37p1WEAzzCY/ZX8KfHg3h6XMHTozKuToVx9XLrTcfT4wqKx5WOz5WrCvZPqzgu12zNulKpFObn54c+67rI3Lhxo6UfoSRJRn3XZpaXl5FOp30pWgDFuEaefokW+kwrwv8sXgku4ZXgEl4OLGAMExjDBF4OLBjH9b8I/7Mt7yWcIzzGYyIUHOjv1cgEPjp7CfsnFRyVasbxo1IN+ycVfHT2El6NTJx73yCCCWh387FYDOvr6w59G9bJ5/Mt8aRsNouFhQVwHIf5+Xlks9lz79nY2DBiQQsLC8jn8x3PnU6nkUgkjMeJRAKzs7OYn59HOp02jsuyjOXl5b7n64Qem0omk5AkqeU5WZaxsbHRMz7VLkSSJHXcn5rL5bC2tmZ6XF6DZlwXnPcenwAAHu79KpTSR7WDXBUTrz0CALz7o18H1LGW9/DjH2Jy7r/ivccn+PSrLEdLmOXF6XEAwKODxs3Ko4MSXpoZN55zg1QqhYWFBUiSdK790DAoFou4d+8eotEo8vk8lpeXsbu7a8wIE4kEdnZ2kMlkIIoidnZ2Wloo6eiitbu7CwCIx+OQZRkPHjxAsVjEwsICFhcXEY1GEY/HUSwWsbe3h0gkglu3bpkaa6FQwM7ODlKpFHZ2drCwsIB79+5hZWXFWP5LpVKmZ0m68MVisXPXKRaLvo4FknCNEO2xKqB/vOqdZ+8AACYjOwgH/77xhPJxAMDEtfN3i+VaBQBwMuREAqI37eLltmgBjRhTKpXyhGNsXjaLxWIQRRH5fB7RaBSyLCOdTmNvb88Q2XYnD2izuEQigVwuh2g0CkmSkM1msb+/D0EQIAgCUqkUtra2IAgC8vl8yznv3LnTcabXzs7OjjHeWCyGt956CwsLC9jd3cX6+rrl1k2FQsEYH6AJmT62QdpAeQESrhGhX6wK6B2vGhsr49KlU1PXqpdUkGQR3UilUlheXjYSBIZNOp1GLpeDJEkty2/5fB6CIPScGeZyOSSTScRiMUPUCoUCAGBubq7ltYuLi4ZY2JltRiKRlsf6DYAuXlbZ2tpCLBZDIpGAJEmIRqO4c+dOy4ytUChgfX0dkiThzp07AIDt7W0jE7H5eX2ZVJZliKLYMXbGChIujEY6vD7TevXyq5gKTxnHFVXBg4MHAIC5mTnwXGscY/9gAu/+47/ElZe3Aci2r1+qlQzhfPOzb2I86O6dPdGfx4clY3kQaMy83J516TOb9fV1pnu4OrGwsIBIJGKIT6e9Zr3I5/NG+ng+nzfEKxqNdhQTMzMrK6ytrRnXtprwks/njaXGbkSjUWP5VH/dysoKrly5Yjx/8+ZN5HK5ltlrMplEsVhsOcYSEi6MVjr8VHgKVy9dNR7X1ToeHj0EAFy5dAUBrjX9uX46CbV6FXaQnpSQe0eLhVWUMg6fa6ncb/3jY4R461lpnRAmQlh6PdL/hUQLzaLVLFSsxCuVSiEejw911iVJEgqFArqVY9WXC3vF41ZXV7G2tgZBEBCPx/HgwQNEo1EUCoWOqeSiKPY9pxXS6bQhXisrK6bPqY+v09JnL/T4XiwW65kqf+fOHczNzZFwEc6wfzCB+umk8VhR66jVtFnWT59Ngm8TLvloElbheW2h8P/+wQ/xje8dage5GsZf1mZ9v/Xn3wDU8/+11Hr4vEhyVUy89mUAwOmPPn8uEQQAMp//ZRIvC3QTrfaYl5vitbKyAkEQWuJc+XzeELNIJGIsw+nOeGNjo8Uxr6ysoFAoQJIkCIKATCZjKW6mL72l02msrq4im82iUCjg5s2bADSRWV1dRTweRyaTQSQSQT6fx/b2tiG4+oxxdXUVmUwGyWQSm5ubxvs2NzchiiKy2SwkScLa2pqRoPHWW2+hWCyaTs4A0CIWejwqlUrhypUrxszITHKGvkxoNd39/v37RpJJLwRBQCQSGdrWBxKuEUHPDnzvvX/RyA40+DQA4EmP9wcC5qNWY6EyKgAufXSr4/OTr//nru998fh3cUlpxAYUVPAopInp3CtT4NFIHjl6XsP7xVPIp5WWc9CyZHe6iZYOS/FKpVJIJBItiQ+iKOLGjRsQBMFYptrb2zPSyfXlKv19W1tbmJ+fRywWOxcD6ocgCFhbWzPSy2/cuHHOmW9ubiKZTGJ5eRmSJCEWi3WdJW5ubmJ+ft4QrGQyaSw9Li4uGu976623EI/HMTs7a8SYzMw8V1ZWjE3Bz549M6pfANqSoSiK+MxnPgNRFPsmV+gJJWbZ2dkxYoGLi4um3zcsqK1JE35ua/Kl7/wFvrL3BeDglxEOtmYOlmoKVFVFOMiD71AzLRCoIRQ+tnS9Ko5RqVVQVVSM8RyCQQ5laMuGYbwErm2LYFk5Rj3wFNOlz+Ij4Z8xjiuo4UP8OQDgo/gfwTfdSx2c1PHg0QT+3af+BX7pWmMJt6KU8X/84LcBAP/zz33p3LKkH5cXf/TsBP/uTwuYHtf2VNmhn2j1e+1ppYbDUg3/6d9E8doV6zNxM7QnGszOzuLBgweGeOjCpe9VWl1dbYnxmHHYy8vLXZcHLwKyLCMej5uuYJJOp7G7u4vNzU0jeaV55pvNZpHL5c7NdmdnZ7G/v+/YuK34X5pxjRjB0AkujbdmB45Dc0p1RcV4KIgAP3jBT6V2CfVaGJeCPMLBAFQoqJzN6UKYOSdcHB/ACZ7icPybOMQ3O55TQttM7TJw+ePAl7f/LdS/aYo99lmWVOth3P+f/rXvxAvQqlfARs5mc7WMqfEgTiu9zzE1HkS5FsKH+89Rrim4ejl0dm330ZfDZFlGJBKBIAjG3i+dvb09xONx3L9/33CYenacF/aHeRlBEGyX3TL73W5sbBhZiMOAhOsCwAGYCGnO7KRSw+SA4lWu1VGuKQifiZYZxnAZk+p8yyxt7Ezwus3UFJRR4n6CiWt/0vW83ZYldx9dw9Lr/4OFTzVcLoeDmJ0Yw/5pFeW6NQHRaw9GLocQCvI4LJkTvlCQR+RyCKeVOn56VMb4WACzE2O4HHbXLdy/fx+RSATb29uGg11dXcXGxgay2ayxoTcWiyGZTBqZevPz810daz6fx+LiIjKZDJWbsoAkScZWgWw2ey4DUZIkbG1tGc8Xi0XjxmNYiRkALRUCaE2Hf/fdd329VBh8/klcHu8sSioaMy+74mVHtHqdIxTkcIwfAgAu4+PnZmpVHENF6zaFulJHmdPEbpx7GXzTe0plHsq4hN+c/wP89qd+3db4hoWfq8Obxe6epH7E43Fks1kj/kPi5T9oqdAio5QO34tBZ15OiBYA473lmgIV6Pm/cAyXWx7XFRXPKxVw4QB4jkO4bVmyoqpgs+DlPFcuh10VDa/gRkVyv1eCIKxBRXYvGLp4BXgOJ2ezLzM4JVo64WAA4SCPSg0I1eYxhU+cm221U1dUnFRqCPBcxyQTwtvoqe3NBWkJwg4047qAWJ15OS1aOs0zr+bHnWgWrYnQGDh8wrFxEGyIRqOOZqERFxeacV1QzM683BItHX3mVa4pKNc6l9xqFa0gaK5FEBcbEq4LTD/xclu0dHqJF4kWQRDtkHBdcLqJFyvR0ukkXiRaBEF0gmJcxLmYV5DnUVPYiZZOc8yrrgA1RSHRIgjiHDTjIgA0xAvQBCPIsxUtnXAwYAgnzsZEokUQRDMkXNA2IF+/fh1LS0vDHspQqTTFl2qKYjpV3knqimqIVvuYCIIgABIuANoG5HfeeQfb29vDHsrQaI5pTY2PWd7n5QTNMa2p8bG+2YaEN9jY2ADHcR2bNKbT6aE3kxyU5eVlJJNJJtcqFApIJBKYn58Hx3GYnZ0Fx3HG96sXH77oUIyL6JiI4WRtQzN0SsSwss+rH35ohVIsFXFSORnqGCZDk4iMWy9OLAgCCoXC0PozjQLJZNJoZbK+vo5sNou9vT0AWs1AvRNzLpezXUR3VCDhuuB0yx50ujBvL3plD9oVr+buzIC5Ds3DbIdSLBWR/HYSclkeyvV1hLCA1K+kLIuXKIoQRdFwrIQ10um0IU4Azt0AiKJoVMnXe3Y5XTbLT5BwXWD6pbyzEC8zKe9WxEvvzvyN7x3ia3/bVMyVq2LiNa3dy+0//buOnZaB4XVbPqmcQC7LCPNhhIPDqVdYrpUhl2WcVE5szbpSqRTm5+dp1mWDGzdutFRblyTpXKV2QFu2TKfTF1q0AIpxXVjM7tOyW9vQDFb2aZmpsAEAwZAmTi+/9ADitQ+Mv7mPPUQ4VEE4VMHcxx42Hf8As1f2cGl6D9zY05Zuy6VaCTe/fhM3v34TpVrJqY/dk3AwjEvBS0P5G1QwRVFELBbD+vq6Q9+GdeLxONLptBFby+fzxnOJRAKzs7OYn59vqZeYzWaxsLAAjuMwPz9vtFGxgh6bSiaTLX3FAK2o8MbGRs/4VLsQSZLUMVksl8thbW3N8vhGDZpxXUCsbi52Y+ZlZ3OxmZkXB+3Y4aWvdz3P08ttRV6ntB/CZQB7h9ewjJfMfASiA6lUymgKOYyGj7IsY3NzE7IsI5VKIRaLAdAETZZlPHjwwOj3tbi4iGg0imKxiHv37iEajRodlHd3d03PGguFAnZ2dpBKpbCzs4OFhQXcu3cPKysrxvJfKpUyPUvShU8fe/N1isXiuU7EFxESrguG3YoYTorXIBUx+omX3rCyvYdXe8NKgEelWkddVTEeDKBSDUAZl/C8dmjrMxEa0WgU0WgUqVRqaA5WkiQ8ePDAEAq9CeL+/j4EQYAgCEilUtja2kI0Gm1ZoovFYhBFEfl83rRw7ezsGOeIxWJ46623jL5j6+vrlluuFAoFY5z6+PP5PPb29qh9yxkkXBeIQcs4OSFeTpRxMiNe7ahQUMGTs+dn8LyiQFVUTJ19hpqP+3h5jVQqheXlZSPRgDWxWKxldlMoFAAAc3NzLa9bXFw0/p1Op41OwO1Lff2IRFrjgbpw222aubW1hVgshkQiAUmSEI1GcefOnXOfaX19HZIkIZFIANBmm6IodoyNjRokXBcEp2oPDiJeTtYeHCRV/rRSg6JwHcf+o+N38a0PvgUAqNQrOKocAQDefvg2QoHQuXPNhGcQfZESEZrRZy3r6+tD2cPVaYkyGo12FZGFhQVEIhEkk0nEYrGO+9Gssra2hmQyaWnmppPP542lxm5Eo1HcvHkTuVyuZcaYTCZRLBZbjo0iJFwXAKcL5toRLzcK5toRL0VVAUXF5dBYy5g5VfspvPX4q3jr8VfPve+Nb7/R9ZxvfvZNEq82UqkU4vH40GZdzUSjURQKhY4p5JIkoVAoQFWd3WifTqcN8VpZWTEd79PH2R7fMsudO3cwNzc38sJFWYUjjltV3q1kG7pZ5d1stiHAg6+I4MoiLodC54SWV6dw/MM38EpwCT9/5efx81d+Hj8X+TmMB8YxHhjHz0V+zjiu/716+VUAwEH5wMFPNBqsrKxAEISWOFc+n8fs7CzS6TSy2ayxFKazsbGBbDZr/AGaI89ms8jn88aSmFVEUcTq6iri8bhxvWw2i42NDWOZT88yzGazxtKiFWRZbvmce3t7RlLG8vJyy/O90JcJ7aa7C4KASCRi6zP4CZpxQatVePfuXdTro1VayO3WJGZmXixak/SbeanQlgfritpzdqhWr0IpvY76mFa9QlHrgPKh9jlOX4XKtZ63Vn0C4AMHP8lokUqlkEgkjNmGvoR448YNCIKAaDSK5eVl7O3tGcKhL4/p79va2sL8/Dxisdi5WJIVNjc3kUwmjWXAxcVFI9NvbW3NSGW/ceOGZeFYWVkxNgU/e/bMqH4BaEuGoijiM5/5DERR7JtcMYhAXyQ41ek5so85PDzEzMwMDg4OMD09PezhWOJL3/kLfGXvCwg+/yQuj3NM+2l1EwbW/bQ6fWazolUpX8bjB/+69SBXxcRrXwYAnP7o8+c2LfPjH2Jy7o/xW9fXsbr0rwYa+8PDh/idt39n6BuQy0oZ//Ff/kdcm77myjXaExZmZ2fx4MEDJJNJLC8vG8Kl73laXV1tiRWNcladLMuIx+OmK49ks1nkcrlz2Zuzs7PY3993Y4iuYsX/0oxrBGHdBLLTzAsA8yaQ7TOvUDBgSrQAIBQ+xotzX0O93vhJqKijGtKWAa9c+0tjj5hOuaZtVj4p1wYe+2RoEkJYgFyWUa6UBz6fXYSwgMnQpKvX0GNNsiwjEolAEARj75fO3t4e4vE47t+/bzhmfWlxGPvDWCAIwsDlsjY2NnDnzh2HRuRdSLhGjHJNgQL2TSDbxQvAUJpANouXLmBmMx9D4eMOR8+y4iaOzj1TL6kYXLI0IuMRpH4l5dsiu1a4f/8+IpEItre3DUe9urpqxLj0DcKxWAzJZNKId83Pz4+saFlFkiRsbW0Ze9SKxaJxQzDqiRkACdfIEhpCE0itojuP04oWKwwH+aE0gQwFA4ZoBXne9cr2ThEZj7guGl5Ad6zt6d6dShl5ISvRi5iJl40ylFU4YvCc5qRPKzWwDl7WFRWnlTp4jgPPcTit1Jk3o9RjWgCMTsrUz8tbmM2wI4hukHCNGGMBbWOtJiLsxKs5EWMyHMRk2J3CvL1oT8SYCJlNlSdYUCgUIElSS4FbgrADLRWOIAFeE6+TSg2nlZrrcaZu2YMsm1F2yx50shklMRjRaNSX2W6E96AZ14iii5fbM69eKe9utkRppl/Ku/lNygRB+AESrhHGbfEys0/LbfEyu0+LxIsgRgcSrhHHLfGysrnYLfEyK1o6JF4EMRqQcF0AnBYvOxUxnBYvq6KlQ+JFEP6HhOuC4JR4DVLGySnxsitaOk6JV6fxl2ol3Pz6Tdz8+k2UaiXb5yYIojskXBeIQcXLidqDg4rXoKKlM6h4lWt11BjvUesGiSVx0SDhumDYFS8nC+baFS+nREvHrnjptSCDPqnI4TYbGxvgOA6zs7OYnZ0Fx3GYn583CuWaZWFhgSqju8gofb+0j+sCYnWflxtV3q02o3RatHSs7vNqLmBcrXGO1Sr0O4IgtOzRKhQKuHXrFra3t02XJmpvT084yyh9vzTjgtaP6/r161haWhr2UJhhdublZmsSszMvt0RLx+zMi3XVfT+jtyCx0phxZWXFduffi0Q+n8f8/Lzl91n9fu1ehwU04wJw+/Zt3L592+gH41VKtRI+983PAdDaxY8Hxwc6X7+ZF4t+Wv1mXm6Llk6/mVcv0fr/PvwBvvQd7d81pYrHxzIA4I//n68hyLf28AKAj1yO4N/80q84/Am8hyiKEAQBOzs7iEajwx4OMULQjOuC023mxbIJZLeZFyvR0uk28+omWvzZr+f7pf8TX9n7Ar6y9wX8yYN/j2flH+NZ+cf4kwf/3jje/Pe/fe82/vS733b1s3gBSZIgyzIWFxeNY8lkEvPz85idnT0Xb1leXkYymWx5rR4zW1hYQD6f73m83/nj8Tg2NjaQSCQwOzuL+fn5lvd2o31chUIBHNf4v2jmvBsbG5ifn+845ub3NddxjMfjSKfTSKfTxjnj8TiWl5chSRI4jgPHcUbR4mw2i4WFBSPGqLeD6fY5eo2723USiQSWl5dbztv+fbCAhIs4J141RWHeBLJdvGqKwlS0dNrFq9dMayLEIVT+JILPG3+B578AKJcA5RICz3+h5bng808CJ78AAPjpcZHJ5xkW+Xze6Gisz7bi8TgKhQJyuRwePHiAYrF4zgk2vz+bzeLBgwdQVRWpVAqRSKTrcTPnl2UZyWQS8XgcDx48QDQadSRZod95E4kEtra2kMlksL+/j1QqZYhNPB6HJEl48OABcrkcksmksbQqyzI2NzeRSqWQSqUQi8WQyWSQyWQgiiJUVYWqqkbcqlgs4t69e1BVFZubm8b3YWfc3a6TSCSQz+dbKvxvbm6ea1HjNrRU6FG23y9CPq20HKsoZRw+rwIA3vrHxwjxjRbv0pPB0qBblw211iSsm0Dq4nVSrhk9vViKls75ZcPuMa2JUOvYVHDQ21FeHufAtX2DxyV+JBM6ZFluuevWnaDeY6tQKCCbzWJ/f99wtJlMBrOzs8jn8+diL7Iso1gsGq/Vn9ebJrYfN3v+aDRq/LvT7MEu3c4ryzLS6TT29vaMJpj66/QmkPqYBUFAKpXC1taWIfa6qJlJqmhuIBmLxSCKIvL5fM9lWqvfRzQahSiKSKfThm3v37+Pe/fu9R2fk5BweZDt94uIf/lvzj/BVTHx2ikA4Paf/h2gNuIn/PghJucAnq8BOB9XIQg3ac8qbGdnZ8eIeTWzuLiIXC53TrhisRgikQg4jkMsFkMikTCSCzodN3v+5mVLfabmBN3Om8/nIQhCx87N+mxobm6u67lisZilTMB0Oo1cLgdJkiBJku1x90KfQa6traFQKECWZZpxEYB8WgE39hQvRU5xKdQ4rqIGOaTNwl782ENwTeYrBz7AIYBg6BSA9QST5piW3sWYRUuUZvSYlqKqmAhpXYxZtERpp3l5UHtMLVEGxWrzSEEQsLe3ZzjieDyOVCqFtbW1jsetnNcN7J43Go1id3e36/OdBK8bCwsLiEQiSCaTiMViWFhY6PseO+NeXV1FMpmEJEnY2tpiLloAxbg8yd7hO7j88f+A48h/wpPLjb+nl9OoBZ6gFniCp5fTLc8dXvo6AICDdefanogR5HnmzSjbEzGCPM+kJUo77TEtqm3oDLFYzEjWaGZnZ6fnNpTV1VVkMhlsbm5ia2ur63G757dDsWg+PhmNRiHLcsfZTzQaNWYsgyJJkhHfc3tLgSAIiMViRrxxGJuaSbg8yPPaIQBAOX0dY9VrmFBfx4T6Oi6pr4JXQ+DVEC6prxrH9b9JdR5juGzpWt2yB1n18wK6Zw+y6uel0y0Rg8RrcPRYymc+8xlDYOLxOERR7HjHns1msbGxAVmWIcsycrkcRFHsetzq+a0giqKxrCdJUktmnpn3rq6uGkkYsiwjm80imUyee675c/c7p36ufD4PSZKMZT49K9HK/jkr19FJJBLIZDIoFotD2XtHwuVhAuoU6rUpqLXLCGEGIcyAQwAcAsbj5j+nRMu4PgPx6pfyzkq8+m0uNiteHHhM4ROYwifAWfh5XYR6g/psYGFhAXNzc4hEIl2XyURRRC6Xw9zcHGZnZyHLMu7du9f1uNXzWyGRSGBnZ8dIsU8kEpaW8DY3NxGLxbC8vIzZ2Vlsbm7i5s2bxnPRaBQLCwvGc/2EIBqNIhqNYm5uzlgmFQQBa2trRmq7/l0MsjTa6To6elzxxo0bts8/CJyqqt6oFOoB9A3IBwcHmJ6eHto4vvSdv8BX9r6A4PNPYiyoGA41FORwjB8CAC7j45YcYztW9mm5tafLyj4tN/d0WamIMWj1jOOSitqlv8dvzv8BfvtTv24cH2RzudMb0wnCDAsLC0aavhNY8b804/I4bixTWRUiN2ZeVoXIrZmXVSGiZUOCgBG3G1aJLhIuH9DsLJUBJ8h2Z09Oipfd2ZPT4mV39kTiRVxU9ESS9fX1ln1jrKF0eJ8QDgagQkEJ0MTLxkrZoEt+VqvKd2LQJT+rVeW7MeiSn9Wq8m4yHhzH1r/a6v9CghgQvexUvzR+tyHh8hHhYABlVXP+5Vod40HzE2an4lSDiJdTcapBxcupKu9eEi+CYMHKygq8kBZBwuUjOPCY4X7WcLwc6qacpdPJFXbEy+nkCrvi5XRrErvi9aPjd/GtD75lPK7UKziqHAEA3n74NkKBUMvrZ8IziL5IFdYJAiDh8iVWnKVbGYFWxMutjECr4uVWPy0r9lAU7bm3Hn8Vbz3+asfXvPHtNzoef/Ozb5J4EQRIuHyLGWfpdmsSM+LldmsSs+LldhNIs/aoVSdx+sM38Il/9jZmZ06N5xRVwYODBwCAuZk58FxjGfiofIQPjj/AQfnA8XEThB8h4fIxvZwlq35avcSLVT+tfuLFqnOxGXtwHAe1ehXTYy/g6qWTxvNqHQ+PHgIArly6ggDXNs5jEARxBgmXz+nkLFk2gQQ6ixfAtglkN/FiJVo6/ewxFuRxDEA+mmx5n6LWUatps6yfPpsE3yRch9VTEATRgIRrBGh2lnUFqCkK0yaQQKt4nZS1jlOKyrYJZLt4BXkeNYWdaOn0skdV1b6b995rq9zNVTHx2vcBALs/+FRby5oPMTm3jfcen+DTr7L5DAThZUi4RgTtzl5zkgCYN4EEcOacA0YTyIlQgHkTSF28jkpV1BQFQZ6taOl0s0cofIwX576Ger31p6eijmpIi2FdufaXLVX+yzWtlY1+Q0AQF52RrJyht6QetDqyn6grquEkAaAyhIoO2v6yxhjKNYVJS5R2mj97TVGYtURpppc9QuFjXJqQ2/4OwPN18HwdlyYOWp4Lhk46XYIgLiwjKVw7OzuO9LjxC80xlKnxsaGUI2pPxGDdz0unOaY1NT7GvJ8X4A17EMQoM1ThKhQKHbt0SpKEjY2Nlt47Vhi0nL+f6JSIwbqWXqfsQZb9vHTaEzFY9/MC7NvDbjsUgriIDC3Glc1mWxq0NROPx406WJIk4datW8hkMqyH6Hl6ZQ+yKkfUK+XdidqGZumWPehUbUMzeMEeBHERGJpwdetK2t7iWhRF5PN543E2m+3YBntlZcVScze/Yybl3W1naWafFgvx6pfyzkK8vGAPgrgoeC6rMJ/PG22odSKRCAqFAqLR6MBtuEcBK/u03HKWVjYXuyleZvdpuSleXrAHQVwkPCdc3eJZxWLR9Dny+XzLEmQ02rm+W7lcRrlcNh4fHh6avsawsLO52GlnaacihhviZXVzsRvi5QV7EMRFw3PC1Q0rCRqxWMxUZ8719XV88YtfHGBUbBmkIoZTznKQMk5OipfdihhOipcX7EEQFxHPpS8JgnBudlUsFl3JErxz5w4ODg6Mv4cPHzp+DadwoozToNmGTtQedCLbcNAyTk5kG3rBHgRxUfGccHWbKS0uLjp+rXA4jOnp6ZY/L+Jk7UG7ztLJgrmDiJdTtQcHES8v2IMgLjKeEK7mZcD2zEBJkrC4uHhh9mW140bBXKvO0o0q73bEy+mCuXbEywv2IIiLztBiXPl8HrlcDoAWa1paWjIyBjOZDJLJJJaWlrC9ve36Hq67d+/i7t27qNe95TQUVcVJpe5KwVyzMRY3W5NYiXm5VeXdSszLzar7dmJepVoJn/vm5wBoTSbHg+MOjoggvAunquowysl5ksPDQ8zMzODg4GCoy4Zf+s5f4Ct7X0D96BcwHlKHsnEXYNdPq58gsGhN0u+zsmoV0+mzHpdU1C79PX5z/g/w25/6deO1JFzEKGHF/3piqZDoDM+535qk2zIVK9ECei8bsuqn1WvZkGV/M1o2JIj+kHB5mLEAx6Q1SbuzZClaOp3Ei3UTyE7ixbopJ0DiRRD98M0+LsJdmmMsepyFZRNIoDXmdVSqno2LbT+t9piXPi7W/c2a7VGtAbTLiyAa0IwLWnLG9evXsbS0NOyhDJVQk0AEeZ55E0hAE4kg3/hvGRrCxlytontjDOEgz7wpJzCcz04QfoCEC8Dt27fxzjvvYHt7e9hDGRr68iAAo+X9MJapyrW60bkYZ2NinT2kLVfWwXMceI7DaaXOvBllsz14bhiySRDehYTLwyiMEj7bY1oToeHEWJpjWhOhwFCaUTbHtCbDQUyG2fbzAs7bYyxAwkUQzZBweZhqXXXdWXZLxGCdINApEYN1M8pOiRism1EOIzGGIPwGJWd4GJ7jXG1+2M9JsioE2yt7kFUzyl7Zg6yaUfazx4+O38W3PviW8bhSr+CocgQAePvh2wgFQi2vnwnPIPpi584IBOFnSLjg3coZYwEOHO+OeJm9s3dbvMykvLstXmZS3t0Wr1724FTtZ/rW46/ircdf7fj+N779Rsfjb372TRIvYuQg4YKWnHH79m1j57aXcMNZWl2Ocku8rOzTcku8rOzTcku8+tmDV6dw/MM38Il/9jZmZ06N44qq4MHBAwDA3MwceK6x8n9UPsIHxx/goHww8PgIwmuQcHkcp52l3RiK0+JlZ3Ox0+JlZ3PxsOyhVq9CKb2O+tiJcUxR64DyofZZTl+FyjW+x1r1CYAPbI+LILwMCZcPcMpZDhr4d0q8BqmI4ZR4DVIRg7U9AgEtLf699xbaBlLFxGvfBwDs/uBTgDpmPMWPf4jJuW289/gEn37V8tAIwtOQcPmEQZ2lU9lqg4qXE2WcBhUvJ8o4sbRHKHyMF+e+hnq99eeqoo5q6AAqgMkX/y+MBccQPDtPuVYBAJyUa1Y/GkF4HhIuH2HXWTqdYm1XvJysPWhXvJysPcjSHqHwcYfzKKhDSygav3yASg3gz77bekkFSRYxqtA+Lp9hdV+RW/uCrO7zcqNgrtV9Xm4UzPWrPQjCz5BwwX+1Cs06S7c3s5p1lm5WeTcrXm5WefeiPViXqCIIlpBwwZ+1Cvs5S1YVGPqJF4vWJP3Ei0VrEq/Zo0bCRYwwJFw+ppuzZF02qJt4seyn1U28WPbT8pI9gi5fgyCGCQmXz2l3ljVFGUqtu3bxYt0EEjgvXjVFYd4E0iv2oBqHxChDWYUjgO4sT8o1nFa0Gc8wCrSezzZk2wQSaM821FqTsG4CycoeHHhM4ROOnpMg/ADNuAiCIAhfQcI1AugxFEVVMREKMO8fpdO6PDic1OzWmFYAisq2nxfgHXsQxKhCwgX/pcM30x74D/I80/5ROu0xrWHsK2pPxAjyPPNmlF6xB4kkMcqQcMGf6fBA92w11s0PuyVisBSvbtmDLJtReskelA5PjDIkXD6lX4o1K2fZL3uQhXj1S3lnIV5eswelwxOjDAmXDzG7L8htZ2k25d1N8TK7T8tN8fKiPSgdnhhlSLh8htXNrG45S6v7tNwQL6ubi90QL7/agyD8DAmXj7BbgcFpZ2nXSTopXnYrYjgpXn63B0H4FRIunzBo2SCnnOWgTtIJ8Rq0jJMT4jUq9iAIP0KVM3yAU7XuBm1+6JSTHKQZpVO1BwdpRukne0hPSsi988jSuISJEJZej1h6D0GwhITL4zhdoNWus3T6zt6OeDldMNeOePnFHjyvtZH8xvcO8bW/3bU8rsznf5nEi/AsJFzQNiDfvXsX9br3GvC5UaDVqrN0aznKini5VeXdini5VeXdDXsEQ6eoAHj9hTomlWnTY/nJQQlPjyuQTysWPwVBsIOEC9oG5Nu3b+Pw8BAzMzPDHo5Bta4iAHeqipt1lm7HUMyIl9utScyIl9utSdyyx8xkABGMWxrL02MSLcLbUHKGh1FUd1th9EsQYBX475WwwaqfVq+EDVb9tLxiD4LwOiRcHmYswLm+kbSbs2TtJDuJF8smkEBn8WLdBNIr9hiUUq2Em1+/iZtfv4lSrTTs4RAjBi0VehieY1P9oH2ZKsjzqCnsnWTzsmFdAWqKwrQJJNC6bHhS1hIc3J75tuMVexCEV6EZFwGg4SwBTTCC/HCcpNZ2XnPUOBsT6+JFzS1RmluTsMQr9iAIL0LCRRhUmuJLNUUZSmuMuqIaotU+JlaoaCSL4Ozfw6i17gV7EIQXMS1cd+7ccXMcxJBpjqFMjY8Npflhc0xranxsKM0o22NarPt56XjBHgThVUwL18bGBg4PD88d/5mf+RlHB0Q0qNZVJs6yPfDPun8U0DkRg3Uzyk6JGCz7eel4wR4E4WVMC5eqdv6x7O3tOTYYohUWbee7ZauxdJa9sgdZiVev7EGW4uUFe7RDGYKE1xg4xsUxyny7iIwFOFedZb8UaxbO0kzKu9viZSblnYV4ecEelabYHkF4FUrO8DA8556zNLsvyE1naWWfllviZWWflpvi5QV7nFZq+MkBzagI72NauDiOG9nZ1d27d3H9+nUsLS0NeyjncMNZWt3M6oaztLO52GnxsrO5eFTtcVqpYe/JCUJBupclvI/pDciqquL111/vePzKlSstxziOw9OnTwceHCu8WqtQZ5AWHO3YrcAwaAuOZgapiDFIS5RmBqmIMWr20EVrfCyA6XCQahUSnse0cKVSKTx79szNsRA9cMJZDlo2yAln6UQZp0HFy4kyTqNij2bREq9O4tlJ2fI4CII1poXrd3/3d90cB2GCQZylU7XuBnGWTtYetCteTtYe9Ls92kWLdXUQgrAL1Sr0GXacpdMFWu04SzcK5loVLzcK5nrVHv2i106LVuFxAQflA+NxpV7BUeUIAPD2w7cRCoTOvWcmPIPoi9GBrktcTGwJ11//9V8jl8tBkiTIsoxoNIpf/dVfxac//Wmnx0d0wIqzdKuquBXxcrPKu1nxcrPKuxftEQqpQIfTK6jg++o6ylwdM2P/K8Sr046I1ue++bmuz7/x7Te6PvfmZ98k8SIsY0m43n//fSQSCeRyOYiiCFEUAQCZTAapVAqLi4vIZDJ47bXXXBks0cCMs3S7FYYZ8WLRmqSfeLFoTeI1e5Rq9Y7CdVqpoczVwXMc5hxaHtRnWq9efhVT4SkAgKIqeHDwAAAwNzMHnmudAh6Vj/DB8QctszSCMIsl4YrFYhBFEXt7e5ibm2t5TpIkJJNJxGIx7O7uYnrafLtwwh69nCWr/k29xItlP61u4sWyn5aX7HFc5VAHUK7WgTHtudNKDdKTE/AvcggHeQQctshUeApXL10FANTVOh4ePQQAXLl0BQGuw2c+dvTyxAXC9KaN3/u934Moivirv/qrc6IFAKIoIpPJ4PXXX8fv/d7vOTpIojud9hWxbjrYaV8R6yaQwPl9XqybQALesUdoTLvOo4MSTs+EVI9phWmvFuFzTM+4/vzP/xyZTKbv6zY3N/Frv/ZrAw2KsEbznf5RqQoAzJsOts+89HGx7qfVPPPSZ18sm0AC3rEHAISCPN57rE1tJsJBvH51HP+V2SgIwh1M33pJkoRf+qVf6vs6URQhSdIgYyJsEOA5BPmGOUNDaDqoVXRvjCEc5Jk3gQRaP3uQ54eS5u0FewDA7EQjm++l6TClvBMjgekZl5WKEl6sPjHqlGt1o1NuTVEGruhgB215TAv8A8BppY7JEMfUWerLgwCM76Jcs19hwy5esAcAPD4sYeJs2fD9Z6eYe2EMOJ+Z3gI39hTfL/6/CH6g/Y77pbb/w9N/cH7gBNED08JlpU7hqNY09CrtMZS6wjtSjsgK7TEtAI6UI7JCp5hWuTZ4eSireMEeiqKlw4eCPD4euQwAkJ6eQHpyArwCdDPHc/4BLn/8P+ArewA6dCzqldoe5GlbKMEG0//T9vf3z9Uk7IYsy3bHQ1ikU+DfyVp6ZuiWiOFULT0zdEvEcKq2oVm8Yg89Hf6l6XHjuxCvTmLvaQUHtbo2tg6DUHACAIjwP4uXZzX30C+1HdBESwgL7nwggmjDUq1Cwlv0ylZj5Sx7ZQ86WQi2F/2yB1mJl6fsEdDS4bm2hphzVyfx9wqHcq2OU7WGy6HO64ZK+WXUTye0f6t1QPlQO//pq1A7pLbzY3UgTC1RCDZQrUKfYibF2m1naSbl3W3xMpvy7rZ4ec0ewbEAnncZQ5jXtgxIPz3Bx18IGUu7AKDn1nz4T/8NHpY+qj3gqph47fsAgN0ffApQxzpe/7/7b7+FyDSJF+E+lhelv/vd72JnZwdLS0v4xV/8RTfGxJy7d+/i7t27qNfdaw3vJFb2BbnlLK3s03JLvKzu03JLvLxoj2qf14eDPEJjAew9OcH8C5OGeE1e0sY/89LbCCizAAAVdVRDWoWLK9f+ElxbSY5qeRoHj/97VKvDyZwkLh6WdiLeuHED0WgUa2triEaj+I3f+A23xsWU27dv45133sH29vawh9IXO5tZnW5+aGdzsdPND+1uLna6GaVf7QEAc1cnMX4mXnompk5o/BCXJuSzvwPwfB08X8eliYOm49rfWPhwwE9AENYwLVx/+Id/CEmSsLe3h2KxiPfeew87Ozv4oz/6IzfHRzQxSAUGp5zlIBUxnBKvQStiOCVefrdHgOcg9hAvgvAqpoUrnU7jv/yX/2KUexJFEV/+8pfxZ3/2Z64NjmjgRNmgQZ2lE2WcBhUvp8o4DSpeo2KPdvEqV619FyqA0tnSK0GwwnSMq1PljFgsRuWdGOBkrTu7MRYnaw/ajXk5XXvQbszLT/Y4QWsVGwU11M/SNor4HnhovbsiV1U8Oizhx6X3+25Q1tHtoaqDLnYShDVcqZxBOIcbBVqtOks3CuZaFS+3CuZaFS+/2ENPoPgJ942u55Hwn5sGAWC28VBVAj3XY5rtEQ7yZ7u/CIINrlTOIJyhrqhQXKoqbtZZulnl3ax4uV3l3ax4uVnl3Wl7jOEyJtV5qGhd+lOhoIxHAIAwXgLXpk4qgEoVKNfDCIZUU/aolMg3EGwZuHKGqqodjz979mywkRGoKSouh9yrKt7PWbJoTdJPvFi1JuknXixakzhtjzFcPndMhYIKngAAQpg5J1wAEBoDTtWGPZpqBUMF8JxxqxiCaIcqZ3iYIM+53jupm7Nk2U+rm3ix7qfVTbxY9tPyoj0mQryxbHhaqUFROBItYqhQ5QwPw8oxtDvLcJDHaaXOtJ/WeWcZQLmmML+zbxcv/d8s+2l50R7cuLa6AkXF5dAYiRYxVKicMwGg3VlqrUlYt+HQneVJWRsDwL4JJNBp5sW2CSTgLXsclxXUzzIHL5u0R4ALYOHFBZdHSFxUqIc3QRAE4StIuAgA7TGUABTVmXJEVtBjWoqqYiIUcKw8lFVaY1rOlYeygtfswXMcOGAo9iCIdki4PAwrB9Ee+A/yvKO19MzQnogR5HlHaxuapT0Rw+nahmbwoj04ADzHDe1mgiCaIeHyMDVFdd1ZdstWc7oQbC+6ZQ86XZi3H92yB1mKl5ftAbC1B0F0g4TLwwR5zlVn2S/FmoWz7Jfyzkq8+qW8sxAvsgdBmIOEy8MEzvZxueEsze4LctNZmt2n5bazNLtPy03x8rM9FKpVSDCGhMvjuOEsrW5mdcNZWt1c7JZ4Wd1cTPbQaLZH8543gmABCZcPcNJZ2q3A4KSztFsRw2nxslsRg+yhoduD6pgSrCHh8glOOMtBywY54SwHLePklHgNWsaJ7KHBARh3uSwZQbRD/+N8xCDO0qlad4M4S6dqDw4qXk7VHrwo9uDAYwqfwBQ+0bEoL0GwZiRLPmWzWQDA9vY2lpeXEYvFhjwi57DT/NDpAq12mh86XTDXbjNKpwvm+tUeABXMJfzLyN0+5fN5SJKElZUVJBIJJJPJYQ/Jcazc6btVVdzKnb5bVd6tzrzcqvLuN3sAgKKq1JqE8C1DFa5CoYCFhfOFOCVJwsbGBrLZLDY2NiDLsulzxmIxrK2tGedZXFx0ariewoyzdLsVhhln6XZrErPi5XZrEj/ZQ1FVqBhOAWOCcIKhLRVms1mIoohCoXDuuXg8jt3dXQCa+Ny6dQuZTMbyNTY3N0e6j1ivZSpW/Zt6LVOx6qfVb9mQVT8tv9hDHTsr3+RwNqB8NGn6tYfVU0evTVwshiZcKysrHY9LktTyWBRF5PN543E2mz33Gv18oigajzc2NnDnzp2WY6NIJ2fJsukg0NlZAmybQHYTL5ZNIAF/2EMvmOvY9QI1AMB775lvY8KPf4jJuW289/gEn37VwcEQFwLPJWfk83lEIpGWY5FIBIVCAdFotKvgtZ8jFoshGo0im82aeo+faXaWdQWoKQrTpoNAq7M8KWuOTFHZxlDaxSvI86gpbJtAAt63x3OHrxUKH+PFua+hXu/sTko1BepZhXlFUTEW5FFXqgBgjI0grOA54eoWzyoWi6beL0kS4vE4RFGELMuIxWIjL1yAfmevOUkAzJsOAjBacOhNIPXWJCzRxeuoVEVNURDk2TeBBLxsD3euFQofd31uHMBRSROqEM9jIhTAcUkFSRZhF88JVzfMJmiIooj9/X1Try2XyyiXy8bjw8NDO0PzBHVFNZwkAFRqdeYOW8X5lveBEM/cYVeakiNqioK6wjMXUK/a41KIfT5WJ3sQxCB47n+QIAjnZlfFYhGCIDh+rfX1dczMzBh/165dc/waLGiOoUyNjw2l+WF7Igbr/lE6zTGtqfGxoVQx97o9WNLNHlSYlxgEzwlXt83CbqS137lzBwcHB8bfw4cPHb+G23QK/LNuftgpe5Bl/yid9kSMYbTg8IM9WIlGL3tU6yRchH08sVQoy7Ixo2rPAtT3Yrkx4wqHwwiHw46flxW9stXsVHSwQ6+Ud7sVHezQLXvQboUNO/jFHsc428vFYWj2qFSpojxhn6EJVz6fRy6XA6At2S0tLRlJFJlMBslkEktLS9je3ra1h2vUMZNi7bazNLNPi4V49Ut5ZyFefrIHDw6Kqs2EJ0OhodhjLMBRcgZhG05VabH57t27uHv3Lur1Ot59910cHBxgenp6aOP50nf+Al/Z+wKCzz+Jy+Pn3YrVfUFu7GWyurnYrb1MVj6bWxuiyR4NzH6245KK2qW/x2/O/wF++1O/7tDVCT9zeHiImZkZU/7XczGuYXD79m2888472N7eHvZQ+mLH4TgdY7EjAG7EvKwKgBsxL7JHA9abvYmLCwmXjxjkLtkpZznIrMVJZ2nXSTopXmSPBiRaBEs8kZxB9MeJpZ1BYyxOLLU5EfMa1Ek6EfMiezQYxB7SkxJy7zwy/XphIoSl1yP9X0iMNCRcPsDJeIRdZ+lkfGgQZ+nUnf0g4kX2aGDXHjyvpWZ843uH+Nrf7loab+bzv0zidcEh4UJrcobXcCOIbtVZupHUYMdZOr0cZUe8yB4NBrFHMHSKCoDXX6hjUmkE4lUFeHRUQqWm4OWZcYSCjWjG0fMa3i+eQj6tWLoWMXpQjAveTc5QVPeqipuNsbjZmsRKjMWtGIqVmJebVd4vqj1mJgP4yNS48ffizDh+4ZUZCBMh7J9WcTkcNJ6bukT32YQGCZeHqdZVV6uK93OWLPppmXGWbgf+zYgXi9YkZI/GGMSrkxgfC2DvyQnzMlWE9yHh8jA8534rjG7OklUTSKC3s2SVrdZLvFj20yJ7NMZA4kV0g4TLw4wFnG341412Z8nSSep0cpasU6w7iRfrJpAA2aN5DM3iValRmShCg4QLWnLG9evXsbS0NOyhDI1mZ3lUqjJ1kjrNzvKoVB3KvqB28WItWjpkj8YYdPH6sVxidl3C25BwwbvJGawJNTmkIM++hxWgOaog3/hvGWLoJHW0iu6NMYSD7HuKAWSP5jG8NO3fYtiE85BwEQAaMRQARst7lv2jdMq1utG5GGdjYl1MU1seq4PnOPAch9NKnWk/L4Ds0cxppYb3n50iFBjG7QPhRUi4PAyrvkntMZSJENv+UTrNMZSJUGAozSibY1qT4SAmw2z7eQFkj2ZOKzXsPTnB+FgAL89cYnRVwuuQcHmYal113Vl2C/yzbn7YKfDPuhllp0QM1s0oyR4NmkVLvDoJjrwVcQb9V/AwPOeus+yXrcbKWfbKVmPlLHtlD7ISL7JHg3bRGkZ8j/AuJFzwblbhWIBzzVmaTbF221maSbF221maSXl3W7zIHg1ItIh+kHDB21mFbjhLq/uC3HKWVvYFueUsrezTcku8yB4NSLQIM5BweRynnaXdzaxOO0s7m1mddpZ2NheTPRo4bQ8SLcIsJFw+wClnOWgFBqec5SAVGJxyloNUxCB7NHDKHiRahBVIuHzCoM7SqbJBgzpLJ8oGDeosnSjjRPZoMKg9SLQIq5Bw+Qi7ztLpWnd2naWTte7sOksnaw+SPRrYtUe5WifRIixDwuUzrDpLtwq0WnWWbhRoteos3SiYS/ZoYEe8Hh2USLQIy5Bw+RCzztLtquJmnaWbVcXNOks3q7yTPRqYtYdy9h2FgjyJFmEZEi54dx9XL/o5S1atMPo5SxatMPo5SxatScgeDczYo3Q2tpemx0m0CMuQcMHb+7h60c1Zsu7f1M1Zsuzf1M1ZsuynRfZo0NcenPYdcCRahA2Cwx4AMRi6szw96x01EQqgXFOY92/SHWG5qdkf6/5NurM8qdRwWqkhHORxWqkz7adF9mjQyx7BsQCeMxkFMYrQjGsE0J1lcwsO1k0HgfY7ffZNB4H2O32tNQnrJpBkjwZesAcxepBwEQRBEL6ChGsE0GMoiqpiIhRg3j9KpzWGwr5/FNAe0wpAUdn2jwLIHs14wR7E6EHC5XPaA/9BnmfaP0qnPfDPun8UcD4RI8jzzJsfkj0a9LJHpcq+mzMxOpBw+Zhu2Wqsmx92y1Zj6Sy7ZQ+ybH5I9mjQ1x5n3b1VxrNQYjQg4fIp/VKsWTnLfinWLJxlv5R3FuJF9mhgxh7jZ2N7dFhivoRK+B8SLvhvA7LZfUFuO0uz+4LcdJZm92m5KV5kjwZm7cGffUeVmgLp6QmJF2EJEi74awOy1c2sbjlLq5tZ3XCWVjcXuyFeZI8GdjZ7vzQzjlK1TuJFWIKEy0fYrcDgtLO0W4HBSWdptyKGk+JF9mhg1x7hsQDmX5gk8SIsQcLlEwYtG+SUsxy0bJATznLQMk5OiBfZo8Gg9pgIBUm8CEuQcPkAp2rdDeosnap1N4izdKr24CDiRfZo4JQ9SLwIK5BweRynC7TadZZOF2i14yydLphrR7zIHg2ctgeJF2EWEi6P40ZVcavO0q2q4lacpVtV3q2Il1tV3skeDUi8CDOQcHmYal11rUCrWWfpdisMM87S7dYkZsTL7dYkZI8GJF5EP0i4PIyiultVvJ+zZNW/qZezZNVPq5d4seqnRfZoQOJF9IKEy8OMBTjXW2F0c5Ysmw4CnZ0lyyaQQGfxYt0EkuzRoF28VKX/e4iLATWS9DA8x6ZrUXvzwyDPo6aw79/U3PywrgA1RWHaBBJobX54Uq4BcH/m2w7Zo4EuXntPTnBSqjK6KuF1aMZFAGg4S0BzUEGefdNBQHOWuqPG2ZhYNx1sbsHR3JqEJWSPBhOhIF6/MoFKnZYLCQ0SLvivVqFbVJriGTVFGUpcoa6ohpNsHxMrVJxveT8Ml0n2aIzh0WGZ+XUJ70LCBX/VKnSL5hjK1PjYUJofNsdQpsbHhtL8sD2mxbqflw7ZozEG6ekJStU6XhHGmV2X8DYkXB6mWleZOMv2wD/r/lFA58A/6+aHnRIxWPbz0iF7NMagi9b8C5MIBcldERqUnOFh9DbnbsYVumWrtScIuJmc0CtbrTlBoPmx0/TKHmxO2CB7sLFHu2hNhII4PkuW+e7DA0vnEiZCWHo9YnkM2+8XIZ9WmFyLsAYJl4cZC3CoK4przrJfijULZ2kmxdptZ2km5Z2FeJE9GmNoFy0ACPHajOvut35o+ZyZz/+yJUHZfr+I+Jf/xvJ17FyLsA4Jl4fhOQ6XXHKWZvcFueksrewLcstZWtmn5aZ4kT0aY+gkWgAgTIawcE1ApSVZRMFPDkoIBXm8NDUOrm018eh5De8XTy3PnPTXX70cwsszvWNrqgI8OirhtFzHabVu+VqEdUi4PI4bztLqZlY3nKWdzaxOO0s7m4vJHg2ctkcv0dIRJkPnjkUmQ9h7coLDcg3i1cm276I00JhenhnHR6b6J4VcnQrj+z8+wGmVfdblRYSinT7AyQQBuxUYnEwQGKQCg1MJAoNUxCB7NHDKHmZEqxteKA8V4Dm8ZELgCGcg4fIJTjjLQcsGOeEsnSgbNKizdKKME9mjwaD2GES0dLwgXu3LlIR70FftIwZxlk7VuhvEWTpZ686us3Sy9iDZo4FtezggWjpeEC+CDSRcPsOOs3S6QKsdZ+lGgVarztKNgrlkjwZ2xOvRYckR0dKhwrwXAxIuH2LFWbpVVdyKs3SzqrhZZ+lmlXeyRwMr9gC0rECnREunWbweHQ2WnEF4ExIun2LGWbrdCsOMs2TRCqOfs2TRmoTs0cCMPSpn2XcvzYw7Klo6unhVajTlGkVIuHxML2fJqn9TL2fJsn9TN2fJsp8W2aNBX3uo2rjCY+59FxOhYN89WIQ/IeHyOZ2cJeumg52cJeumg8B5Z8m6CSRA9mimlz3GGbVoofqGowltQB4BmjfFHp0122PddLB9U6w+Ltb9m5o3xeobY1k2gQTIHs10s0edcX8zYrSg2xGMRj+uAM8hyDfMGRpC00GtgnhjDOEgz7zpIND62YM8z7wJJED2aMYL9iBGCxIujEY/rnKtbnTKBcC8fxSAs+WxOniOA89xOK3Ume+l0ZejABide1n2j9Ihe2h4xR7EaEHCNQI0x1AmQoGhND9sjqFMhoOYDLPtHwWcT8SYCLHtH6VD9tDoZY8qiRcxABTj8jmdAv8s+0cB3bPVWPWPArpnD7LqH6VD9tDoZ4/njMRTPq0yuU4zrPqFXWRIuHxMr2w1Vs6yV7Yaq+aH/bIHWYkX2UPDjD1qNQ5VaO1DIhOODwEA8PiwxFS4WPYLu+iQcPkUMynWbjtLMynWbjtLsynvbosX2UPDrD3GgoEz4arica2EF6ed3W/1+LCERwclCBNjODrrnOw27f3C+vUKA+z3C7vokHD5ECv7gtxyllb2BbnlLK3u03JLvMgeGnb2zQkTY3j0TCvL5JR46aL10sw4OA54uP/ckfOaQe8XdlqpYe/wBMJEqEOPsGaoJJUdKDnDZ9jZzOpk/yjAXgUGJ/tHAfYrYjjVP0qH7KFh1x7CRAgvzYzj0UEJjw8Hd+LNouX0LM4sp5Ua9p6cYHws0Ee0CLuQcPmIQSowOOUsB6nA4JSzHLQihlPiRfbQGNQeL06POyJeJFoXBxIun+BE2aBBnaUTZYMGdZZOlXEaVLzIHhpO2WNQ8SLRuliQcPkAJ2vd2XWWTta6s+ssna49aFe8yB4aTtvDrni5KVpm+3mRaLGFhMvjuFGg1aqzdKNAq1Vn6VbBXKviRfbQcMseVsXL7ZnWo6NS3++CRIs9JFwepq6orlUVN+ss3awqbtZZul3l3ax4uVnlnezRwKx4sVgerNQUSE9Pun4XJFrDgYTLw9QU1dWq4v2cJYtWGP2cJavWJP3Ei0VrErJHg37ixSqm9fLMOErVekfxItEaHiRcHibIc66XKOrmLFn2b+rmLFn30+omXiz7aZE9GnQTL5aJGKEgj/kXJs+JF4nWcKENyB6G1Y+hfVNsOMjjtFJn2r+pfVPsRCiAck1h2gQSOL9JWf83y35aZI8GujA9OmgIF+vswYlQEPMvTGLvyQmkpyd4aTqM95+dkmgNERIuAkC7s9RaYbBuOqg7y5OyNgaAfRNIoFOFDbZNIAGyRzPt4jWMlHddvN57fAzpSY1Ea8jQUiFBEAThK0i4CADtMZQAFJVt/yigEUNRVBUToQDz/lE6rTEt9v28ALJHM80xLSfLQ1lBj2lNhIMQX5hEpd4725BwFxIuD8PqR9Ee+A/yPPPmh+2B/yDPO1pLzyztiRhO1zY0A9mjQXsihlPloazQnogxNT7WMWGDYAcJl4epKarrzrJbtprThWB70S1bzelCsP3olj3IUrzIHg26ZQ+yFK9KTemYPajHvEi8hsNIClc2m0U+n0cymYQkScMejm2CPOeqs+yXYs3CWfZLsWblLPulvLMQL7JHg34p76zE6ycHpa6JGCRew2PkhEuWZWxvbyMWi2FpaQmpVGrYQ7JNgOdcc5Zm9wW56SzN7gty21ma3aflpniRPRqY3afFQrxCQb5n9iCJ13AYqnAVCgUsLCycOy5JEjY2NpDNZrGxsQFZlk2fUxAEQ6xyuRwSiYRTwx0KbjhLq5tZ3XCWVjezuuUsrW4uJntouGUPq5uL3Ravl6bG+34XJF7sGdo+rmw2C1EUUSgUzj0Xj8exu7sLQBOxW7duIZPJWDp/Pp+HIAgQBMGJ4Q4VJzv32q3A0L4pdpA9RXYrMDjdudduRQyyh4bT9rBbEaN9n5eTe7w4k7f27ZuUaY+XuwxNuFZWVjoeb49JiaKIfD5vPM5msx3jVisrKxBF0Xgci8UQiUSQSCSQy+UcGvXwcMJZDlo2yAlnOWjZIKec5aBlnMgeGk7ZY9AyTm6Kl1lIvNjhucoZ+XwekUik5VgkEkGhUEA0Gu0qeDrpdBqyLGNtbQ2CIPg6OaOdQZylU7XuBnGWTtW6G9RZOlV7kOyhMag9nKo9SOJ1cfBccka3eFaxWDT1/hs3bhiztM3NzZ5LjOVyGYeHhy1/XsdOjMXpAq12YixOF2i1G2NxumAu2UPDrj3k04qjtQeHsc+rHYp5uY/nZlzdMJugIQiCMSuLxWI9X7u+vo4vfvGLgw6NOVbu9N2qKm7lTt+tquJW7/TdqvJO9tCwM/OST6v4qMO1B5tnXpdCbGtM6tDMy108N+MSBOHc7KpYLLqSZHHnzh0cHBwYfw8fPnT8Gm5h5k7f7VYYZu703W6FYfZO3+3WJGQPDbP2qJ59R8LEmCtLevrMSz6tOn5us5iZealKhzcSffGccHWbJS0uLjp+rXA4jOnp6ZY/P9HLWbLq39TLWbLq39TPWbLqp0X20DBjj+rZMWEi5MoYAE28hIkx185vhl7iVVdUPDoaznKm3/GEcDUvAzZnBgJaluHi4uJIpLW7QSdnybLpINDZWbJuOtjNWbJsAgmQPXT62WOM0bLZsIUL6CxedUWF9PQElRpNuewwtBhXPp830tTX19extLRkxKYymQySySSWlpawvb1teQ+X19h+vwj5tGL69dITa3dhzTGWugLUFIVp00GgNcZyUq4BABSVbdPB9hhLkOdRU9g2gQTIHjq97MEFA2C5iPfdhweuvr4fzTGvH/70GABQqSt4eWYch6VjR691EeBUVb3wKS93797F3bt3Ua/X8e677+Lg4MCxZcPt94uIf/lvLL2HH/8Qk3N/jFD5OiZC5u8YTyt11BTtDm5qfIxp00GdmqIYTQcnQgEEefaTehXAUUlzi1pV8+EE6MkeGp3sUcEBTrn38XH1txDBP3ft2vJJBbsPZdvvX7gmQJh0bjnzqFSF9OQEACC+MInn1Tr+4cND3Pu3C1i+/pJj1/Ejh4eHmJmZMeV/fZNV6Ca3b9/G7du3jS/OSfSZ1uuRCUxdMvd1n/Af4AmAYOgUgLnx1BXVcJIAUKnVmXftVXG+5X0gxDN32JWm+FJNUVBXeOYZXWSPBp3swSpIIUyGsHBNwHGlhp8clKDfpgsTY32XEUM876ho1RUVjw7LxuNHh2VMh8kF24G+NUZMXQriI1PmsqeKCOCJhXO3x1AqZ3EEYLByRFZoj6EAcKQckVWaY1qhYMCxckRWIHs06GaPUEgFGOl4aIzHvlyFMBGCeHUST4/LeHRQwswldzIaO6HHtErVOn7mxcsAgL0nJzit1Jhcf9TwRHIGYZ9OgX/WzQ87Bf5Z9o/SaU/EYN0/CiB7NNPLHiVGTTnbm0AGeI75JuVm0Zp/YRIToaAR86LkDHuQcPmYXtlqrJxlr2w1ls6yW/YgS/EiezToaw9OG1e56t530Um0dFiJVyfR0pkIBfHyDPuyVKMACRe05Izr169jaWlp2EMxjZkUa7edpZkUaxbOsl/KOwvxIns0MGOP0Jh2/NFByZXlsl6ipeO2ePUSLZ1QkFywHehbg5ac8c4772B7e9u1azi5JGBlX5BbztLKviA3naXZfVpuihfZo4EVewCa43Y61mNGtHTcEi8zokXYh4SLET9x6M7SzmZWp52lnc2sbjhLq5uL3RAvskcDO5u9X5oex/hYwDHxsiJaOk6LF4mW+5BwMcKJO8tBKjA45SwHqcDgpLO0WxHDSfEiezSwbQ+eg3h10hHxsiNaOk6JF4kWG0i4GPHS1GB3lk6UDRrUWTpRNsgJZzloGScnxIvs0WBQewQcEK9BREtnUPEi0WIHCRfYJGdwPGz/OJ2sdWfXWTpZ624QZ+lU7cFBxIvs0cApewwiXk6Ilo5d8SLRYgsJF9gkZwD2fpxuFGi16izdKNBqx1k6XTDXjniRPRo4bQ87vw8nRUvHqniRaLGHhIsxVn6ciotVxc06Szeriltxlm5VebciXm5WeSd7NMZg9vfhhmjpmBUvEq3hQMI1BMz+OEu1uqtVxfs5SxatMMw4S7dbk5gRLxatScgejTH0+324KVo6/cSLRGt4kHANiV4/Tr2aQIBzvxVGN2fJsn9TL2fJqp9WL/Fi2U+L7NEYQ7ffBwvR0ukmXiRaw4WEa4h0+nGeVmp4dKD9QEJjASbFUNudJeumg0BnZ8m6CWQn8WLdBBIgezSPodPvg5Vo6bSLF4nW8KFvHK39uFij/zilpyd477HWUC40xaMKMHGSOs3ND/VK5iwrqgOtzQ/1/k2sm0C2Nz/Ux8WyojpA9mgeQ/vvYyIcZCZaOnoV+UcHJTw6KIHnORKtIUIzLrDLKuxGgOfw0nTYeDw74VwPICuEmhxSkGffwwrQvovmZochhk5SR6vo3hhDOMi+hxVA9mgeQ/Pv46Xp8FC+i6uXG2O4HA6SaA0REi4PcFqp4f1npxgfC2B8LMCk1UI7+nIUAKPFOosWHO2Ua1rXYN1ZsmrB0Yy2PFYHz3HgOQ6nlTqTlijNkD0atP8+3n92yryPlb48yPMcpi+N4fB5dSi/U0KDbhmaUM/aox4eHjp2zpPjIyjlU1SfB1AJnv+xPa/U8eDpCcJjAbwqTAAA3j2qoP68jmq9Ap6vOjaWbqgAnp/FUCZCQQRUBajVcfq8ivpZA0AWVGrNTQeBgKJqY3hexSVGS3V6TCfAcwif3VE/r9Rw+LxqxL/c5iLYo44q6lwdVfU5Kjju+rpOv4/3n5Xw7sMTzF2dxKWQ+99FXVHx/rNTlKv1s2uqGFPq+PDxM1RL4/jIVLj/SbpQfV6GUj7FyfERDg8nHBy1/9D9ru6He8GpZl51Qfinf/onXLt2bdjDIAiCuLA8fPgQH/vYx3q+hoSrCUVR8OMf/xhTU1PguGFENZzh8PAQ165dw8OHDzE9PT3s4TgKfTZ/Qp/Nn7D8bKqq4ujoCK+88gp4vncUi5YKm+B5vq/S+4np6emR+yHp0GfzJ/TZ/AmrzzYzM2PqdZScQRAEQfgKEi6CIAjCV5BwjSDhcBi///u/j3DYfraTV6HP5k/os/kTr342Ss4gCIIgfAXNuAiCIAhfQcJFEARB+ApKhx8hCoUCbt26hd3d3WEPxXEKhQLy+TwAYHt7G/fu3YMgCMMdlEPon0uWZWxvb+PmzZuIRqNDHpXzJJNJ3LlzZ2TsVigUAADRaBSSJEGW5ZGyWz6fhyRJEEURABCLxYY8oiZUYiTIZDLq7u6uOqomTaVSLf+ORqNDHI2zCIKg7u7uqqqqqpubm6ooikMekfPo/zf39/eHPRTHWF1dVaFV6FJjsdhIfbZcLqeurq6qqqqqe3t7nvs/SUuFI8LKyspI3e01UygUsL6+bjxeWVlBoVCAJElDHJVzZDKZFtuNyoykmeY791FhYWEB+/v72N/fRy6XGym7JRIJpFIpAIAoisjlckMeUSskXITniUajuHfvnvFYlmUAQCQSGdKInKV5CSaTySCRSAxxNM6TzWaxsrIy7GG4giAIIyVYgHaTUSwWIQgCCoUCZFn23E0HxbgIX9Ds+La2thCLxUbKYRQKBWxtbWF5eRmrq6vDHo5jyLI8UnZqRpZlZLNZAFrcNZFIeM7B26FQKCASiSCbzSIWiyGdTkMURU/dfNA+rhGD4zhTbQH8iizLWFhYwO7u7sg5RFmWkUwmsby87CknMQjpdNoQ4vn5+ZGyW7MoFwoFxONx7O3tDXdQDpBOp5FIJLC/vw9BECDLMmZnZz3lV2ipkPAVyWRy5OIJOoIgIB6PIx6PG8uhfiafz+PGjRvDHoZrNMdYRVGEJEkjEXcVRbFlCbRZnL0CCRfhGzY2NpBMJiGKImRZHhnnPjs7azzWl5pGwQECwP3795FOp5FOpyFJEtbX1z3lAO1SKBTwmc985tzxUYi7+mG5k2JcI8goxhWy2Syi0aghWvfv3x+JWFAkEmlJzigUChAEYSQyRNv3/SQSiZGJA4miaGTdAdoNyMrKykj87kRRxOLiouFH9IxQL/2fJOEaEfL5vJGyur6+jqWlpZGJk0iShHg83nJMEISREK5oNIqbN28inU4DAHK53MhtIJdl2fh8qVQKiUTCU07QDoIgYHFxERsbGxAEAXt7e8hkMsMelmNkMhkkk0kjnuy1dHhKziAIgiB8BcW4CIIgCF9BwkUQBEH4ChIugiAIwleQcBEEQRC+goSLIAiC8BUkXARBEISvIOEiCIIgfAUJF0GMGKNQCosgekHCRRAeR5ZlbGxstPxJkmS01GgnmUy2iFehUEAikQDHcUgmk+fOHY/HMTs7i0QiQaJH+AIq+UQQHkaWZdy6daulnJAkSVhYWMBbb73V9X3NNfOi0ajRnDKdTrfU2BMEAffu3UM6ncba2przH4AgXIBmXAThYeLxeIvQAI0iqJ3q/WWzWSwvL587vrOzg1QqBVmWkc/nzz03CnUfiYsDCRdBeJh8Pt+xmnp70WGdra2trsWVBUHAysoKNjc3W45LkjQSVc2JiwMJF0F4GFEUz8WlAHScIcmy3LcfVCKRQDabpVgW4WtIuAjCw2xubiKdToPjOCwvLxvtQTpx//59I5bVjCRJWFxcBKD1yBIEoed5CMLrkHARhIeJxWLY399HJpOBKIpIJBIdZ2CA1surU9wrn8+3HF9dXTWWCwuFwrmGjwThdUi4CMIH6LGpXC7Xcbakd6k1QyKRgCRJKBQK2NnZGYmOxMTFgoSLIDxItzjU4uJixzjW5uZmx2XCTuht2NuTNAjCL5BwEYQH6dYqfWdnp2PWYLcZlyzLHY8nEgncv39/8IESxBAg4SIID7Kzs3NuSVCWZWxubp7b11UoFDru3QK0Dcd6YkYzq6urkGW543ME4XWocgZBeJCbN28iFothY2MDgiAYy4bNFTR0OomZLMtIJpO4f/8+tre3ce/evXN7tVZXVzsmcxCE1+FUVVWHPQiCIOyTSCQoXkVcKGipkCB8TLcSTwQxypBwEYSP6VXiiSBGFRIugvAxtAeLuIhQjIsgCILwFTTjIgiCIHwFCRdBEAThK0i4CIIgCF9BwkUQBEH4ChIugiAIwleQcBEEQRC+goSLIAiC8BUkXARBEISvIOEiCIIgfMX/D8Uqv09DLcXdAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 440x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpl.rcParams['hatch.linewidth'] = 1\n",
    "\n",
    "# Density vs SNR\n",
    "# --------------\n",
    "plt.figure(figsize=(4.4,4))\n",
    "plt.title('Glass-z12')\n",
    "\n",
    "# plt.grid(True, alpha = 0.4, lw=1, ls=':')\n",
    "\n",
    "# plotting the Poisson uncertainty\n",
    "PTD, bin_edges = np.histogram(data_real_P, bins=np.linspace(0, 6, 26), density = True)\n",
    "PFD, bin_edges = np.histogram(data_jack,   bins=np.linspace(0, 6, 26), density = True)\n",
    "bin_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])\n",
    "bin_widths = bin_edges[1] - bin_edges[0]\n",
    "\n",
    "raw_counts, _ = np.histogram(data_real_P, bins=np.linspace(0, 6, 26))\n",
    "uncertainties = np.sqrt(raw_counts)\n",
    "density_uncertainties = uncertainties / (len(data_real_P) * bin_widths)\n",
    "\n",
    "plt.fill_between(bin_edges, np.append(0,PFD), step = 'pre', alpha = 0.7, color='C0', edgecolor='C0', lw=1, label = r'N$_{\\rm pos}$ jacked $\\propto\\mathcal{P}_{\\rm FD}$', hatch='x', rasterized = True)\n",
    "plt.fill_between(bin_edges, np.append(0,PTD), step = 'pre', alpha = 0.7, color='C2', edgecolor='C2', lw=1, label = r'N$_{\\rm pos}$ real $\\propto\\mathcal{P}_{\\rm D}$', hatch='', rasterized = True )\n",
    "plt.step(bin_edges, np.append(0,PFD), where = 'pre', alpha = 1, color='C0',  lw=1,)\n",
    "plt.step(bin_edges, np.append(0,PTD), where = 'pre', alpha = 1, color='C2',  lw=1,)\n",
    "plt.errorbar(bin_centers, PTD, yerr=density_uncertainties, fmt=' ', alpha=0.8, color='C2', label='Poisson uncertainty')\n",
    "\n",
    "\n",
    "plt.xlabel(r'$S/N$', fontsize = 12)\n",
    "plt.ylabel(r'PDF', fontsize = 12)\n",
    "plt.legend(frameon=False, fontsize = 12, loc = 1)\n",
    "plt.semilogy()\n",
    "plt.axis(ymin = 5.5e-4, ymax = 1e0, xmin = 0.1, xmax = 6.5) # added this\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25c370b2-0cd5-4852-9606-10f545b068d7",
   "metadata": {},
   "source": [
    "In the five jackknifed data sets we run, we thus find higher SNR peaks in the data cube. Thus we can't assume that the 4$\\sigma$ peak we find in the data is real. We can be more quantitive. Since we sampled the distribution of false positives with large statistics (preferable use more jack knifed data sets than 5), we can treat the resulting pdf as a sampled posterior distribution. Hence the integral from $\\gamma=3.8$ till the highest bin gives the likelihood of having at least one peak above the detection threshold of $3.8\\sigma$ which is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f081dedf-53a6-43e9-8953-477da28c39f0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0065\n",
      "0.0064\n",
      "0.9839743589743585\n"
     ]
    }
   ],
   "source": [
    "snr = 3.8\n",
    "fd = np.sum(bin_widths*PFD[bin_centers>snr])\n",
    "td = np.sum(bin_widths*PTD[bin_centers>snr])\n",
    "\n",
    "print('{:.4f}'.format(fd))\n",
    "print('{:.4f}'.format(td))\n",
    "print(td/fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "405de4c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Peaks: N= 936\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of Peaks: N=\", len(data_real_P))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74715b82-869c-4467-8b97-5fcbdb24fb53",
   "metadata": {},
   "source": [
    "But with N, amount of peaks in the data we expect to find:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d185abd2-04c3-4c4d-8c7b-e9299661e789",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "~ 6.097719869706839 ± 2.469356165016873\n"
     ]
    }
   ],
   "source": [
    "print('~',     fd * len(data_real_P), \n",
    "'±',   np.sqrt(fd * len(data_real_P)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e6358ce-7fbb-4f5a-a78e-8bbf8b4eb26e",
   "metadata": {},
   "source": [
    "and we find in the real data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2ad510ab-4ff6-4416-a0ab-b79d46c3f193",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "~ 5.999999999999995 ± 2.449489742783177\n"
     ]
    }
   ],
   "source": [
    "print('~',     td * len(data_real_P), \n",
    "'±',   np.sqrt(td * len(data_real_P)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee3721db",
   "metadata": {},
   "source": [
    "Hence, the data is completely consistent with noise."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "interfero",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}