{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Removing noise from Kepler, K2, and TESS light curves using Cotrending Basis Vectors (`CBVCorrector`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cotrending Basis Vectors (CBVs) are generated in the PDC component of the Kepler/K2/TESS pipeline and are used to remove systematic trends in light curves. They are built from the most common systematic trends observed in each PDC Unit of Work (Quarter for Kepler, Campaign for K2 and Sector for TESS). Each Kepler and K2 module output and each TESS CCD has its own set of CBVs. You can read an introduction to the CBVs in [Demystifying Kepler Data](https://arxiv.org/pdf/1207.3093.pdf) or to greater detail in the [Kepler Data Processing Handbook](https://archive.stsci.edu/kepler/manuals/KSCI-19081-003-KDPH.pdf). The same basic method to generate CBVs is used for all three missions.\n",
    "\n",
    "This tutorial provides examples of how to utilize the various CBVs to clean lightcurves of common trends experienced by all targets. The technique exploits two goodness metrics that characterize the performance of the fit. [CBVCorrector](https://docs.lightkurve.org/reference/api/lightkurve.correctors.CBVCorrector.html) inherits the [RegressionCorrector](https://docs.lightkurve.org/reference/api/lightkurve.correctors.RegressionCorrector.html?highlight=regressioncorrector) class in LightKurve. It is recommend to first read the tutorial on [obtaining the CBVs](https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-2-how-to-use-cbvs.html)  before reading this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cotrending Basis Vector Types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three basic types of CBVs: \n",
    "\n",
    "- **Single-Scale** contains all systematic trends combined in a single set of basis vectors. \n",
    "\n",
    "- **Multi-Scale** contains systematic trends in specific wavelet-based band passes. There are usually three sets of multi-scale basis vectors in three bands.\n",
    "\n",
    "- **Spike** contains only short impulsive spike systematics.\n",
    "\n",
    "There are two different correction methods in PDC: Single-Scale and Multi-Scale. Single-Scale performs the correction in a single bandpass. Multi-Scale performs the correction in three separate wavelet-based bandpasses. Both corrections are performed in PDC but we can only export a single PDC light curve for each target. So, PDC must choose which of the two to export on a per-target basis. Generally speaking, single-scale performs better at preserving longer period signals. But at periods close to transiting planet durations multi-scale performs better at preserving signals. PDC therefore mostly chooses multi-scale for use within the planet finding pipeline and for the archive. You can find in the light curve FITS header which PDC method was chosen (keyword “PDCMETHD”). Additionally, a seperate correction is alway performed to remove short impulsive systematic spikes.\n",
    "\n",
    "For an individual's research needs, the mission supplied PDC lightcurves might not be ideal and so the CBVs are provided to the user to perform their own correction. All three CBV types are provided at MAST for TESS, however only Single-Scale is provided at MAST for Kepler and K2. Also for Kepler and K2, Cotrending Basis Vectors are supplied for only the 30-minute target cadence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the CBVs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can directly obtain the CBVs with `load_tess_cbvs` and `load_kepler_cbvs`, either from MAST by default or from a local directory `cbv_dir`. However when generating a [CBVCorrector](https://docs.lightkurve.org/reference/api/lightkurve.correctors.CBVCorrector.html?highlight=cbvcorrector) object the appropriate CBVs are automatically downloaded from MAST and aligned to the lightcurve. Let's generate this object for a particularily interesting TESS variable target. We first download the SAP lightcurve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:02.244402Z",
     "iopub.status.busy": "2023-11-03T14:15:02.243753Z",
     "iopub.status.idle": "2023-11-03T14:15:07.470631Z",
     "shell.execute_reply": "2023-11-03T14:15:07.470208Z"
    }
   },
   "outputs": [],
   "source": [
    "from lightkurve import search_lightcurve\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "lc = search_lightcurve('TIC 99180739', author='SPOC', sector=10).download(flux_column='sap_flux')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create a `CBVCorrector` object. This will download the CBVs appropriate for this target and store them in the `CBVCorrector` object. In the case of TESS, this means the CBVs associated with the CCD this target is on and for Sector 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:07.472645Z",
     "iopub.status.busy": "2023-11-03T14:15:07.472525Z",
     "iopub.status.idle": "2023-11-03T14:15:09.219667Z",
     "shell.execute_reply": "2023-11-03T14:15:09.219376Z"
    }
   },
   "outputs": [],
   "source": [
    "from lightkurve.correctors import CBVCorrector\n",
    "cbvCorrector = CBVCorrector(lc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the CBVs downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:09.221582Z",
     "iopub.status.busy": "2023-11-03T14:15:09.221466Z",
     "iopub.status.idle": "2023-11-03T14:15:09.224804Z",
     "shell.execute_reply": "2023-11-03T14:15:09.224509Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[TESS CBVs, Sector.Camera.CCD : 10.1.1, CBVType : SingleScale, nCBVS : 16,\n",
       " TESS CBVs, Sector.Camera.CCD : 10.1.1, CBVType.Band: MultiScale.1, nCBVs : 8,\n",
       " TESS CBVs, Sector.Camera.CCD : 10.1.1, CBVType.Band: MultiScale.2, nCBVs : 8,\n",
       " TESS CBVs, Sector.Camera.CCD : 10.1.1, CBVType.Band: MultiScale.3, nCBVs : 8,\n",
       " TESS CBVs, Sector.Camera.CCD : 10.1.1, CBVType : Spike, nCBVS : 6]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cbvCorrector.cbvs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that there are a total of 5 sets of CBVs, all associated with TESS Sector 10, Camera 1 and CCD 1. The number of CBVs per type is also given. Let's plot the Single-Scale CBVs, which contain all systematics combined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:09.226423Z",
     "iopub.status.busy": "2023-11-03T14:15:09.226327Z",
     "iopub.status.idle": "2023-11-03T14:15:09.834925Z",
     "shell.execute_reply": "2023-11-03T14:15:09.834620Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqUAAAGRCAYAAABc28oUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADESklEQVR4nOydd5wURfqHn+7JM5tzIGfJIIggogJmMJ2CImbFLOrpmbOengFz/KF4BjgwJxAJihKUJDmnZdmc08QOvz9mp9lhAwss7AL1fMTtrq6qfrt7uvpbb73VLem6riMQCAQCgUAgEDQjcnMbIBAIBAKBQCAQCFEqEAgEAoFAIGh2hCgVCAQCgUAgEDQ7QpQKBAKBQCAQCJodIUoFAoFAIBAIBM2OEKUCgUAgEAgEgmZHiFKBQCAQCAQCQbMjRKlAIBAIBAKBoNkRolQgEAgEAoFA0OwIUSoQ1CAQCDBhwgQ2btzY3KYIBILjlLy8PEaPHs2OHTuatN5XX32VZ599tknrPFTWrl3L6NGjqaysbHSZ++67j0WLFh1GqwTNhbm5DRA0HaNHj25w+xVXXMGIESO48cYb69z+0ksv0a1bN1RV5euvv2bevHkUFBRgtVpJS0vjrLPO4uyzzwagrKyMzz//nGXLllFaWkpERATt27fn8ssvp3v37vXa4Ha7+fLLL1m8eDH5+fm4XC7atm3Leeedx+DBg5EkiYceeoh169YZZWJiYujRowfXX389SUlJvP/++6xatYp33323Vv35+fncdNNNPPzwwwwaNKgxpy2MWbNmkZyczAknnGCkrV27lv/973/s2LEDv99PfHw8J5xwAnfccQcWi+WA97Eva9eu5eGHH2batGlEREQccn37ous6s2fPZs6cOezevRuTyURqaiqnn346Z599Nna7vcn32ZwsWrSIH3/8kR07dqBpGikpKQwZMoRRo0YRGRkJBDsf33//Pb/99hvZ2dnYbDZatWrFWWedxemnn47ZbObVV19l/vz5AJhMJiIjI2nXrh3Dhg1jxIgRyPKh9ekzMjL4/PPP2b59O/n5+dx4441ceOGFtfL99NNPfP3115SUlNC+fXtuvvlmunTpcsj11sTv9/P222+zfft2MjMzGThwII8++uhBHVdJSQkzZsxg2bJlFBUVERMTQ/v27bnwwgvp06cPADfccAP5+fkAyLJMTEwMJ554Itdffz0RERE8/fTTqKrKU089Vav+9evX8+CDDzJ48GCWLFnSoC0//PDDQR3D4SY3N5dPP/2UdevWUVFRQVRUFJ06deKaa66hdevWJCQk8MknnxAVFdXcprJkyRK++uorMjMz0XWdhIQE+vXrx0033dRsNo0ZM4bJkyczePDgQ74PBS0LIUqPIT755BNj+Y8//uDzzz/nvffeM9Lsdjvl5eUAPPvss7Rp0yasfOiBPW3aNGbPns3NN99Mp06d8Hg8bN26Nawn+/zzz6MoCvfccw8pKSmUlpayevVqKioq6rWvsrKSBx54ALfbzfjx4+ncuTMmk4l169bx8ccf07t3b0OUnX322Vx55ZXouk5+fj6TJ0/mlVde4T//+Q9nnnkmP/74Ixs3bgwTjwDz5s0jOjqaAQMGHPD503Wdn376iSuvvNJI2717N08++SSjRo1iwoQJWK1WsrOzWbx4MZqmHfA+Die6rqNpGiaTKSx90qRJLF68mLFjx3LzzTcTHR3Nzp07+f7770lKSmLw4MHNZHHdKIqC2XxwTdMnn3zCV199xYUXXsjVV19NXFwc2dnZzJo1i19//ZULLriAQCDAE088wc6dOxk/fjwnnHACTqeTTZs28c0339ChQwc6dOgAQP/+/bn77rvRNI2SkhJWrlzJ//3f/7Fo0SIee+yxWuf6QPD5fKSkpDB06FAmT55cZ54//viDyZMnc/vtt9OlSxe+//57Hn/8cd577z1iYmIOut590TQNm83G6NGjWbx48cEeEnl5efzrX//C5XJx3XXX0a5dOxRF4e+//+bdd98Na4+uvPJKzj77bDRNIysri7feeov333+ff/7zn5x55pm88MILFBYWkpCQELaPuXPn0qlTJ+655x5uvfVWI/3ee+/l7LPPNjrOLRVFUXj88cdJT0/noYceIi4ujsLCQlasWEFVVRUQ7ATFxsY2s6WwevVqXnzxRa666ipOOukkJEli9+7drFq1qlntOvHEE3nzzTdZsWIFAwcObFZbBE2LEKXHEDUbMafTiSRJtRq2kCiNjIyst9FbunQp5513HkOHDjXS2rdvbyxXVlayfv16/v3vf9OrVy8AkpKSGvTeQFAw5Ofn89577xEfH2+kp6enM2zYMKxWq5Fms9kM++Li4jj//PN5++23AejQoQMdO3Zkzpw5YaJU13XmzZvH8OHDMZlMBAIBPvzwQxYvXkxlZSUxMTGce+65XHbZZXXat23bNnJzc8ME7d9//01MTAzXXXedkZaamsqJJ54YVnb9+vV88sknbNu2jaioKE4++WSuueYawwsZCAT4/PPPWbBgAaWlpSQmJnLppZfSp08fHn74YSDoyQYYPnw499xzD4FAgI8++og//vgDt9tNp06duPHGG43zHPKwPvHEE3z22WdkZGTw9NNPG9cEgqLmt99+45FHHuHkk0820pOTkxk0aBButxuALVu28Omnn7J9+3ZUVaV9+/bceOONdOrUySgzevRobrvtNpYuXcqaNWtISkpi4sSJREVF8eabb7J161bat2/PvffeS2pqqlHuzz//ZNq0aWRmZhIXF8eIESMYM2aMIehGjx7NrbfeyooVK1i9ejWXXHIJY8eO5e2332b16tXG+TrvvPO44IIL6rx2oWP44osvuOmmm8LyJScn069fP6NT9f3337N+/XomTZpEx44djXwhIacoipFmsViM32F8fDydOnWia9euPProo8ydO/eQBFCXLl2Ma/nf//63zjzffvstZ599NiNHjgTgtttuY9myZcyZM6fe33Fj6t0Xu93ObbfdBsCGDRsMcXSgvPvuu0iSxKRJk8I88G3btjWOIYTD4Qg7tyNGjGDBggUAnHTSSURFRTFv3jzGjh1rlPF4PCxatIjrrrsOh8OBw+EwtsmybNQ5bdo0Fi5caLQZIe666y5OOukkxo8fz6uvvkpVVRUdO3bkxx9/JBAIcNpppzFhwgRjBETTNL766it+/vlnSktLSUtL4/LLL+eUU045qPMDwY5uTk4Ozz77LElJSUCw/aw5wpSXl8eNN97I66+/TocOHYx7/dlnn+Xjjz9m9+7ddOjQgYkTJ9KqVSuj3PTp0/nhhx/w+/0MHTqUqKgoVq5cyRtvvFGnLfs7vqVLl3LCCSdwySWXGGXS09NrdWSXLl3KtGnTyMjIwG6306NHDx555BEA5s+fzw8//EBWVhY2m43evXtz00031dupgv23pyaTiQEDBvD7778LUXqMIfzeglrExsayevVqysrK6tweehj8+eefBAKBRtWpaRp//PEHp512WpggrVlnfV6niooKFi5cSNeuXY20M888k4ULF+L1eo20tWvXkpeXx5lnngkEh+7++usvHnjgAd577z3++c9/Gg+Buli/fj1paWk4nU4jLTY2lpKSkrBwgn3JycnhySefZMiQIbz55pv861//YsOGDWFeoUmTJrFgwQImTJjAu+++y+23347D4SAhIYGHHnoIgPfee49PPvmECRMmADBlyhQWL17M3XffzWuvvUZqaipPPPFELW/0f//7X6655hreeecd2rVrF7ZtwYIFpKenhwnSEJIk4XK5gODDfvjw4fznP//h5ZdfJi0tjaeeesoQrSGmT5/O8OHDeeONN2jVqhUvvfQSb7/9Npdddhmvvvoquq6HHff69et59dVXueCCC3jnnXe4/fbbmTt3LjNmzAird9q0aQwePJi33nqLkSNHous68fHxPPjgg7z99ttcfvnlfPLJJ/zxxx/1XofffvsNh8PBeeedV+f2kBf+t99+o0+fPmGCNITZbN5vOEOfPn1o3759g0PHoTi5vLy8ButqiEAgwLZt24whbwgKr759+7J58+aDrvdwUVFRwcqVKzn//PPrPIcNhaYUFRWxdOlS4x43mUwMHz6cefPmoeu6kW/RokWoqsqwYcMatOXMM89kz549bNmyxUjbvn07u3btChPHa9asITMzk3//+9/cf//9LFmyhGnTphnbv/jiC+bPn8/tt9/O22+/zYUXXsgrr7zC2rVr69331KlTueGGG+rdHh0djSzLxrEcCJ9++inXX389r776KiaTKUxs/vbbb8yYMYNrr72WV199lcTERGbNmtVgffs7vpiYGHbv3k1GRka9dSxbtoznnnuOAQMG8Prrr/Pcc8+FOShUVeXKK6/kjTfe4JFHHiE/P5/XXnut3voa055CsPO1YcOGBo9PcPQhPKXHKffff3+tWJwvvvgCCMZ7vfDCC1x99dW0adOGbt26MWjQIMODaDKZmDhxIm+99RY///wzHTt2pGfPnpx66qlhHtWalJeXU1lZGdarb4iZM2fyyy+/oOs6Pp+P9PT0sPiy0047jY8++oiFCxcaD5m5c+fSvXt30tPTASgoKCAtLY3u3bsjSVKDgjSUPy4uLiztlFNOYeXKlTz00EPExsbStWtX+vTpw/Dhww3x+sUXX3DaaacZcXtpaWlMmDCBhx9+mNtuu42CggIWLlzIM888Q9++fYGgVy5EKGwiOjraeHB7vV5mzZrFxIkTjfN+5513csMNNzBnzpwwz8WVV15Jv3796jym7OzsRp3zmsIH4I477uDyyy9n3bp1nHTSSUb6yJEjOfXUUwH4xz/+wf3338/ll19O//79Abjgggt4/fXXjfzTpk3j0ksvZcSIEcZxjx8/no8//tjwDAMMGzaslietZhhFSkoKmzZtYuHChcb+6zrW5OTk/Q79Z2dnh3mTD4ZWrVqxa9euerfbbDbS09MPOgwBgveMpmm1RjRiYmLYs2fPQdd7uMjJyUHX9Ubf4x9//DGfffYZmqbh9/vp2rVrmJgbOXIkX3/9NevWrTOu19y5cxkyZIjRmaqPUNzj3LlzDYE0d+5cevbsGXbvmc1m7rrrLux2O23btuXKK69kypQpjB8/HlVV+eKLL3j22Wfp1q0bEPwdbtiwgZ9//rne31BUVFTYPvYlPj6eCRMmMGXKFKZNm0bnzp3p1asXp59+eoPlAK666ipjv5deeilPPfUUfr8fq9XKjz/+yJlnnmncR1dccQV///13WMe9JoFAYL/HN3r0aDZs2MAdd9xBUlISXbt2pV+/fpx++umGN3nGjBkMGzYs7H6t+RwIOQlC9U+YMIF7770Xj8cT5ukOsb/2NDSiFgp70DRNxJUeQwhRepzywAMP1PvwaNOmDW+99Rbbtm1j48aNrF+/nmeeeYYRI0Zw1113AUGxNnDgQNavX8/mzZtZsWIFX331FXfeeWctcQGEeTsaw2mnncaYMWMAKC0t5YsvvuDxxx/n1Vdfxel0EhERweDBg5k7dy4jR47E7XazePFibrnlFqOOESNG8Pjjj3PLLbfQv39/Bg4caIinuvD5fGEhBBAU4HfffTdXXXUVq1evNoaIv/rqK1555RXi4uLYuXMnu3btMoYeQ8eraRp5eXns2rULWZbp2bNno48/JycHRVHChvTMZjNdunQhMzMzLG/nzp3rraex572kpITPPvuMtWvXUlZWhqZp+Hw+CgoKwvLV9MSGht/atm0blub3+3G73TidTnbu3MnGjRvDPKMhEeL1eg2PWl3H8NNPPzFnzhwKCgrw+/0oilJvp+dAjrUp2N++unTpUsuzc6xzoOf/kksuMTorBQUFfPrppzz99NM8//zzmEwmWrduzQknnMCcOXPo1asX2dnZRthQYzj77LN5/fXXufHGG5EkiQULFtSa5Nm+ffswr263bt3weDwUFhbi8Xjw+Xw89thjYWUURTFijuti1KhRjBo1qkHbzj//fM444wzWrVvHpk2bWLRoEV988QWPPvpovR1MCL//Qp2V0tJSkpKS2LNnT61Rgi5durBmzZo668rOzt7v8dntdp544glycnJYs2YNmzdv5sMPP+T777/npZdewm63s2PHjgbDWLZt28bUqVPZuXMnVVVVRix+QUFBrXkNwH7b09atWwNgtVrRNI1AIIDNZqt3/4KjCyFKj1MSEhJIS0urd7ssy0Zs2oUXXsivv/7KpEmTGDNmjNGbt1qt9OvXj379+nH55ZfzxhtvMHXq1DpFaXR0NC6Xq9EeHpfLZdiXlpbGXXfdxdVXX80ff/xhNIBnnnkmjz76KNnZ2axduxZZlsPiYDt16sTkyZNZsWIFq1at4sUXX6RPnz7GcPm+REVF1TtMFR8fz/Dhwxk+fDjjx4/n5ptvZtasWVx55ZV4vV7OOeecOt9+kJiYSE5OTqOO+WBpqEFOT09v1Dl/7bXXKC8vZ8KECSQmJmKxWLj//vvD4iuBsBALSZIAwryBobTQg8fr9TJu3Lg6J1PtG0Nck99//52PPvqI66+/nm7duuFwOPj666/DhmPrOtaNGzfud6JUWlraIXsa9+zZQ3Jy8iHVsT+ioqKQZZmSkpKw9NLS0hYxCWZf0tLSkCSp0ec2Kioq7B632Wzcf//9rF271hhROPPMM3n//fe55ZZbmDt3LqmpqY3u3J100klYLBaWLFmC2WxGVdUDigUNeRgff/zxWiFHTfHWDafTyUknncRJJ53EVVddxeOPP8706dMbFKV13X8H2xk7kONLTU0lNTWVs88+mzFjxnDLLbcYo1QNtT9er5fHH3+c/v37c9999xEVFUVBQQFPPPFErbalZpmG2tMQlZWV2O12IUiPMYTPW9AoQr1Tn89Xb542bdrUO1QkyzLDhg1jwYIFFBUV1dru8XgajK8KDc/4/X4jrXfv3iQnJzN37lzmzp3LsGHDasWyOZ1OTj31VO68807+9a9/sXjx4nrfENCxY0f27Nmz30Y+IiKCuLg441g7duxIZmYmaWlptf5ZLBbatm2Lruv1xqWGBFTN2fypqamYzeawmClFUdi6datxLRrDaaedRlZWFn/++WetbbquGxNaNm7cyOjRoxkwYABt27bFYrEYk+IOhY4dO5KVlVXnuWloyG3jxo1069aN888/n44dO5KWlkZubm6D+zrttNPweDzMnDmzzu2hiU6nnXYaq1evZvv27bXyKIpS7284xOrVq9m1a9chTXZpDBaLhU6dOoV5ujRNY/Xq1WHx1S2FyMhI+vXrx08//VTnOdzfeyjruseHDh2KLMssWLCAX3/9lZEjRxpibH+YTCZGjBhhtA+nnnpqLQGzc+fOsDZt06ZNRqx369atsVgsRhhQzX81xVFTIEkSrVq1arB93R+tWrVi69atYWn7rtfkYI8vOTkZm81mXON27dqxevXqOvPu2bOHiooKrrnmGnr06EHr1q3rnasQYn/taYiMjIwGPdaCoxMhSo9TKioqKCkpCfsXehg8//zzfPvtt2zevJn8/HzWrl3Le++9R3p6Oq1ataK8vJxHHnmEX3/9lZ07d5Kbm8vChQv56quv6pxQE+Kqq64iISGB++67j/nz57N7926ys7OZM2cOEydOxOPxGHl9Pp9h186dO3nnnXcMz2wISZI488wzmTVrFps2bQqLXYLgzOUFCxaQmZlJVlYWCxcuJDY2tt54tF69euH1etm9e7eRNmvWLN555x1WrlxJTk4OGRkZxuzXUKzlP/7xDzZu3Mh7773Hjh07yM7O5s8//zSGb5OTkxk+fDivv/46S5YsITc3l7Vr1xqTdpKSkpAkiWXLllFWVobH48Fut3PeeecxZcoUVqxYwe7du3nzzTfx+XycddZZ9Z7jH3/80Zj1CsGH+qmnnspLL73EjBkz2Lp1K/n5+SxdupRHH33UEDypqan8+uuvZGZmsnnzZl555ZVaoQwHw+WXX878+fONmbmZmZn8/vvvfPrppw2WS01NZdu2baxcuZKsrCw+++yzWg/YJUuWhIVrdO3alX/84x98+OGHTJkyhU2bNpGfn8/q1at54YUXjHeOXnjhhZxwwgk8+uij/PTTT8Zv+I8//uC+++4jOzvbqDMQCFBSUkJRURHbtm1jxowZPPfccwwcOJAzzjijXvu3bNnCLbfcUmcHrGbdO3bsYMeOHSiKQlFRkfH7CXHRRRcxe/Zs5s2bR2ZmJu+88w5erzdsNGLSpElhs+wbU+++vxMIzgrfsWMHlZWVuN1uo44D4dZbb0XTNO69914WLVpEdnY2mZmZfP/999x///1heT0eDyUlJRQXF7NlyxamTJlCdHS0Ed8IwQmQQ4cO5ZNPPqG4uNgY7m8sZ511FmvWrGHlypW12gcIdkLeeOMNdu/ezfLly5k6dSrnn38+sizjdDq5+OKLmTx5MvPmzSMnJ4dt27bxww8/MG/evHr3Wde5rcmOHTt49tlnWbRokdEG/vLLL8ydO/eg3q0cYtSoUfzyyy/MmzeP7Oxspk+fzq5du+oV8Y05vqlTpzJlyhTWrl1Lbm4u27dv5/XXX0dRFMObfcUVV/D777/z+eefk5mZya5du/jyyy+BoGfTbDbz448/kpuby19//cX06dMbPI79tach1q9f36BXWXB0Iobvj1PqejH2/fffz7Bhw+jfvz+///47X375JVVVVcTGxtK7d2/GjRuHyWTC4XDQpUsXvvvuO3Jzc1EUhYSEBM4+++x6X1MDQU/Kyy+/zJdffsn06dPJz88nIiKCdu3acd1114WJxdmzZzN79mwAI88TTzxRKw52xIgRTJ06lTZt2tTyHoWGfLOzs5Flmc6dO/PEE0/U66ELvXrkt99+45prrgH2zvB85513KC4uxm6306ZNGx555BFjwkH79u15/vnn+fTTT3nwwQfRdZ2UlJSwCTm33XYbn3zyCe+99x7l5eUkJiYaMbPx8fGMGzeO//73v7z++uucccYZ3HPPPVxzzTVomsakSZPweDx06tSJp556qsFZzOXl5WEeRUmSuO+++4yX58+YMQOTyURaWhpnnHGGEWN711138dZbb3H33XeTkJDA1VdfzUcffVTvfhpL//79efzxx/nf//7Hl19+idlsNl5S3xDnnnsuO3bs4MUXXwSCE6HOO+88VqxYYeSpqqoiKysrrNy1115Lx44d+emnn5g1a5ZxLU455RSGDx8OBD2QzzzzDN9++y0///wzH330kfHy/NGjR4fFyK5cuZKrr74ak8lkfCBiwoQJDB8+vEFPr8/nIysrq94hSoDi4mImTpxorH/zzTd888039OzZk+effx6AU0891fhQRUlJCR06dOCpp54KG74vKCgIEx6NqXff3wnAU089ZbzQHjDqCL2APvSaopqvgtuXlJQUXnvtNWbMmMFHH31EcXEx0dHRdOrUyXjlVIjPP/+czz//HAiG93Tu3Jmnn3661gvjzzrrLObMmcOAAQPqfHNHQ6SlpXHCCSdQUVFRp3e5d+/epKWl8eCDD6IoCsOGDWPcuHHG9vHjxxMdHc0XX3xBXl4eLpeLjh07NtjO1XVuaxIfH09SUhLTpk0zzndycjLjxo3b70cOGuL0008nNzeXjz76iEAgwNChQxkxYkSDIS/7O76ePXvy008/MWnSJOMjKR06dODpp5822uJevXrxwAMPMH36dL788kucTic9evQAgtf17rvv5pNPPuGHH36gY8eOXH/99TzzzDP12tSY9rSoqIhNmzbxz3/+86DPl6BlIulHcnaAQNDC2blzJ48//jgffPBBnTNDBYLjlTVr1vDvf/+byZMnH5Yvjx0OdF3n5ptv5rzzzuOiiy4K2xZ6T+nBfrnqaOCxxx4jJibmmBNvH3/8MZWVldxxxx3NbYqgiRHD9wJBDdq3b88111xzSO+XFAiORZYvX86YMWOOGkFaVlbGjz/+SElJSZ2TL481vF4v3377rREm8/nnn7Nq1aoDDnk4GoiOjg57BZXg2EEM3wsE+3A8PMAEggPl+uuvb24TDojx48cTFRXF7bffftQI6UNBkiSWL1/OjBkz8Pv9xmdMQ7GfxxIXX3xxc5sgOEyI4XuBQCAQCAQCQbMjhu8FAoFAIBAIBM2OEKUCgUAgEAgEgmZHiFKBQCAQCAQCQbMjRKlAIBAIBAKBoNkRolQgEAgEAoFA0OwIUSoQCAQCgUAgaHaEKBUIBAKBQCAQNDtClAoEAoFAIBAImh0hSgUCgUAgEAgEzU6jPjOqaRrFxcU4HA4kSTrcNgkEAoFAIBAIjhF0Xcfj8RAXF4cs1+8PbZQoLS4u5rrrrmsy4wQCgUAgEAgExxdTpkwhISGh3u2NEqUOh8OozOl0No1lhwlN0ygoKCAxMbFBNS5o2YjreGwgruPRj7iGxwbiOh4bHK3X0e12c9111xl6sj4aJUpDQ/ZOp/OoEKUOhwOn03lUXTBBOOI6HhuI63j0I67hsYG4jscGR/t13F8I6NF3RAKBQCAQCASCY45GeUoFAoHgYAgEAuTl5eH3+8UkyaMUXdcpKirC7/djMpmIjY1t8SNmAoHg6ESIUoFAcFgoLCzkpZdewuPxIMuyEKVHKbquo2la2DU8+eSTueyyy47K4UOBQNByEaJUIBA0OZqm8b///Y+IiAhuuOGG/Qa3C1o2iqJgNptRVZUdO3bwww8/ADB27NhmtkwgEBxLCFEqEAianPLycrZv38748eNp3749ZrNZeEqPUnRdN0SpJEm0a9cOgB9++IHRo0eLoXyBQNBkiLEXgUDQ5FRVVQE0+D46wdFLhw4dACgpKWlmSwQCwbGEEKUCgaDJ0TQNQMQcHqOYTCZg73UWCASCpkA8MQQCgUAgEAgEzY6IKRUIBMcVM2fOZPbs2WRkZDBmzBjGjRvX3CYdUTIyMnj33XfZsWMHCQkJ3HrrrfTq1au5zRIIBALhKRUIBMcXsbGxjBs3jiFDhjS3KUccRVF47rnnGDJkCNOmTWPChAk8//zzlJeXN7dpAoFAIESpQCA4vhg8eDCDBg3C5XI1tylHnKysLCorK7ngggswmUz07duXjh07smTJkuY2TSAQNAGKohh/CwoKKC0tZcCAASxfvryZLWscYvheIBAIjiN0Xa+1vnv37mayRtBYdF2nqqqKyspKkpOTqaqqYuvWrfTu3ZtNmzbRunVroqKigOCX1AKBQNjruva97scSiqKgKArbtm3j9ttvZ8iQIcyZM4cOHTqwY8eOWvm7dOnCtddeS/v27Wnfvj2aplFaWorL5cJut+P1elmyZAnt2rWjrKyMVq1aER0djd1uZ/ny5XTv3h2n00lubi5FRUV06tQJVVVxu90kJCQYH5tQFMW4Tn379iUQCOBwOHC73dhsNkwmE5s2bWLdunXs3LmTCy+8kOLiYhITE5v8HcC33HKLsbxw4ULsdnuT1t9UCFEqEAgExwnp6em4XC6+/fZbRo0axerVq1m3bh0pKSlH3BZd19m6dSubN2/GYrGwY8cOnE4nb731FgBPP/00P/74I3fddRfr1q2jV69eqKpKly5dyM7OJikpCbfbTSAQIDk5GV3Xyc7OJj4+PuyB6/f7sVqtde5/33fnqqqKJElhb41QVdV420CozO7du0lKSsJisbBhwwbKy8uZM2cOAwYMQJZlfD4fo0ePxmw2U1JSwvbt2+nWrRvr1q3jjjvu4IUXXqBt27ZomkZycjI//PADAwcOxOPx4Pf7uf322+nbt68haA4HQ4YMwe12c8kll/DXX3/x008/AfDggw/ywgsv0L9/f8466yzMZjNbt27F7XYbH01oDE899RRbt25lzpw5DB8+nGnTppGUlER+fn5Yvh49erB+/fomPbY5c+YA1ClIAbZs2cLDDz/cpPtsCqZPn35Y6h08eLAxGnLhhRdis9kOy36aAklvRPfJ7XYzduxYpk+f3uJflKxpGvn5+SQlJYnX0RzFiOt4aGiahtfrNe5Xr9dLSUkJKSkpTfoSe0VR8Pv9fP7550iSxKZNmzj55JNJTk5m6tSpjB07lqSkJDp27Ghcx61bt9KxY0cmTZpEfn4+qqoan7IMYbFYkCTJ8DjURVJSEg899FCj7AyJCV3X0XUdWZZ5++23jfjSI82tX64hq8x70OXTo+28e2nvgyq7c+dO3n//fXbv3k2nTp2Ijo4mLS2NK664wsgTeiyE/gYCAcxms5G2Z88eJk2axA033EBycjKaphn/YmNjKSsrw+l0cs0116DrOrm5uQd9rALBwdCvXz+uueYa5s2bx/nnn09paSmJiYn88ssvZGdn8/fff1NZWXlQdZvNZhRF4eKLL0ZRFMrKyli/fj1nnXUW8fHxKIqCruv4/X58Ph+rV6+mffv29O3bl8LCQtavX0/Pnj2xWCxERkbSs2dPfD4fMTExlJWVkZaWRmRkJIqiIEkSJpPJ+Myvruthz8ZQuylJUov+QEljdWSL9pTOnz//gHpmIXw+X4vqCRQVFfHCCy+QlpbW3KYc92iaht/vb3DoIvQgru8Gr7m9uLgYj8dDWloaqqry888/M2rUqIOyLSSYNE3DbDZTWlpKdHQ0mzZtIjIykosuuuig6t2XefPmUVhYSHx8PDExMfh8PkwmE8XFxWzbto0hQ4agqiqlpaW43W5Wr16NzWajb9++PPjgg6xevbrB+n/77TdMJhPt27c30rZv3x6WZ/v27Vx44YWHfCwbNmwAICoqKmyyjtVqxWw243a76y1bUlKCoijs3r2b2NjYWkKs5nJ9fxubtu+2x09JbHT++sjMzNxvnrowm83cfvvtxvpLL71Ejx496q0v9CAMdQ5Cv/vKykq+/PJLIiIijAenz+fj66+/PiB77r77bl577TUuvvhiLrroItxuN6+//jrnnHMOpaWlfPzxxwd1nIeD3r17U1JSQpcuXWjVqhWdO3fGYrGg6zrR0dEkJiZSVlZGXFwcO3bsIBAIkJubS5cuXSguLiY2NhaXy4Wu68TGxpKdnc3ixYu5+OKL2bx5M5IkERkZSWJiIhaLhaSkJLxeLx6Ph8jISMxmMyaTCYvFUqd9VVVVBAIBYmJiyM3NxeVyIcsysixjs9nIz88nPj4eCHqAZVnGarVSUVGB0+mkpKSEqKgoAoEAmqYhSRJWq9U4Rl3XDa9xY6nLI30kGTp0aNh6nz59msmSxtGqVStjOdQRrMm+bcOx5rRp0aJ0+PDhDB8+/IDKtEQP25QpU7jgggtqpXfs2JHPP/8cRVG47bbbyM3N5aOPPkKSJGw2GxaLBYfD0WAPSFVVFixYAMC2bdu47LLLKCsro6ysjA4dOmC1WiktLSUnJ4eePXuybds24zOBiqIYDdi+Q1QhYRRC13VUVTXSQsHUkiRRVVWFrutERkbi9Xp58MEH2bNnD2+88QaFhYW0atWKyMhILBYLBQUFREVFYbPZjF5e6IHn8/mMzxlaLBZUVUXTNObPn88pp5yCxWJBURR27txJamoqmqaxY8cONm3aRN++fdm0aRPJyckoikJMTAz33nsvFRUVxjE8+eSTPPnkk8Z6QkIChYWFB39h62HmzJm0bdsWRVFYsWIFmqYRFRXF0KFD+fPPP1mzZk2T7/NAGDFiRJPXeeqpp7Jw4ULOOeccRo8ezZ133onNZqNNmzbGg7qhGd6h30JkZCTx8fGUlJRQVlbW6P3vW7ff78fv99eZN/S7Cv2TJAm/32/8/kP3Ws37bt+/daU1tO1A66grT1Oxc+dO0tPT0TSNmTNnYrVaOeecc+rNv+9nRiF4DqOjo5kwYQKtW7cOy38ww6Ljx48PW//ss8+M5TvuuOOA62sJ7Hte6qJVq1acdNJJALRt27bOPHa7nZiYmEbts+bkvX1DMkIetZDXraawjYyMBPZ+ga2ucIeD/R22ZO+doOXRokXpscJll12G0+lkzZo1/P3330ZMzfbt2zn55JPD8h6sly3EBx98cEjlm5Km8uw1FTUFKXBYBCnA0qVLWbp0aa30kFfvQImOjjaC7V966SVMJhOxsbGG5ySEqqrGkH1IaMmyHJZn06ZNPPnkk/Tp04evvvqqzv25XC6io6ORJImsrCwSExNJTExk1KhRZGdnExkZyYknnkhcXBxJSUl1ep2//PJLJk2ahMViwW63ExEREeYB2B9Op5P09PQDOEuNZ+rUqUybNs1Ynz9/PhMnTmTkyJGHZX8tjblz5zJv3jw0TaNv37488sgjzW2SQCAQAEKUHhEiIiIYO3Zso2bT6bqO1+tFURQsFgvl5eVs3ryZiooKNE1j+/btREdHG97Mbt264fP58Pl8VFVVoSgKXq+XNm3aGDFjlZWVrFmzhp49e+L3+4mJiWHGjBmccsopdO3aFYvFQm5uLn6/H4fDQc+ePY1YQYfDgcPhMDxKZrMZj8djzDIMDf9IkmQM+VitVkwmU1gcWs3JBiHPaGhYMOSlDR1zaJgoLy+PxMREzGZz2JDFwfS8Q56ehs57qN7mHm46WEwmU5inpC5vR7du3fjf//4H0Oh4zGONcePGHXcvzK/JTTfdxE033dTcZggEAkEthChtYUiShMPhMNbtdjtJSUmHXO++w3Onn3562PqBxNmEhnr2Zd9Yo5qiqL4YzpDQBYzhpJBYrTkL9lBFYkOCdN/6j0ZBKhAIBALB0U7LCLoUCAQCgUAgEBzXCFEqEAgEAoFAIGh2hCgVCAQCgUAgEDQ7QpQKBAKBQCAQCJodIUoFAoFAIBAIBM2OEKUCgUAgEAgEgmZHiFKBQCAQCAQCQbMjRKlAIBAIBAKBoNkRL88XCATHFYFAgHfeeYdVq1ZRVVVFmzZtuPHGG+nWrVtzm3bYmTlzJrNnzyYjI4MxY8aEfdlq7ty5fPbZZ7jdboYMGcLtt98e9n10gUAgONwIT6lAIDiuUFWVpKQk/vOf//C///2PCy64gKeffhqPx9Pcph12YmNjGTduHEOGDAlL37VrF5MnT+bhhx9mypQpFBYWMn369GayUiAQHK8IT6lAIDiusNvtXHHFFcb6sGHD+PDDD8nKyqJTp05Nth9d19EBTdPRdNB0vfofxt9QHl2vLhMsGEwz6glt06n+rzpfKI8eVr6u/KF6krv0BkBd+CfFbj/bCqvQdfhp1lx69j8JOS6dbDcMPecCpk5+l5PPuRhdB13Xqj+/G/wEb36ph/xKH7d9tRYpMsc45iq/yvxthcb68+d144ZBbYh3WpFl8flegUDQMEKUHgG+WZvDJR8vN9ZdVhMnt42lX3o0j5/ZhSq/QoLLitkkHNcCwZFE13UyMrOoqKggOj6RCp8SFIxauHgME5T7bGsIWQJZkoJ/ZanGenBZkiRCUk2SgpJPkqRqARhax1imjvwghecxtknG8t7ywfRIm5lYh4WO8U4AKoty6d27N50SXAAk27rxdlEh6S4Zu92OoiiYzWbDLrvXQVKEjRf+0YvWrVuHnU/5vh+N9YdmbuKhmZvqPDeL7zyF5Egb7WKdbMiroGdqVMMnUyAQHPMIUXoEeH9JRth6lV9l3tZC5m0t5OXfttdZ5vwTkjghOZKXf9tOhM3Ep1f04/ROCcQ4LLj9CjaziUqfgixJVPoVrCaZeJcVn6JS4VWIdliQJQlTtXdC13Wyy72kRzsAcPsV7GYTsiyhV7tWQg+cmjS0rSF0XafcqxBlNxvli91+ou3BGDVZCnpu6vKehPZ5oPs7UBsbQ0DVUDUdu8XU6DI+RcVmNtVpU0DVsFR3PhRVQ9F0bGYZVdMJaDoOiwlV043rpmk65T6FGEfTx/Zpms7a3HLaxjqJtpsp9yo4LCas5v13jjRNP6Ker4CqGaJKRkInKLYCqha0RwerScYTUFE0HQko8wawm03EuSxIgF/VqfAp7C7ZO0wfZYFXX3qZs0ZdREC2oVbfUyERaZHlcGEZEpRycDkkIo82JClc/Ho8HpxOp7E9tOz1erHb7QdQr4T+ymj8ioaOjv2BmYzpk8aM1dm18g55c1HY+qfj+pFX4eOkNjF0TYzAE1A5b/JfbMir5LYh7dhaWMmcLYW16qlJ79Qo1uSUN9re/eG0mkiPstM50UWlT6FDvIsoe1DQn9EpnjKPwuB2sSzLLCXeaSUpwkZepY84p4XlmaVU+BRK3AHSo+10jHeRVe7llHZxzNyYR/s4Jyv2lPHlmhzeuLgnP27IY8aqbB4c3omthVV8vCyTHY+MoKjKT6TNTHQj2gBN0wloGrIkkVfhA+CzFXsY1T0ZkwTuch/ryguJtltxWGQSXFZyK3yUegIkRti4/NMVrMutYEDraE7vmGA8n05qE8MJSRHcOqQdOpAaaSMt2s7uEg8BVcMT0Oic6KLKr3LH12vpmhRBqSdApM3MvK2FvHJBd56ft41TO8QR77SypaCSO4e2Z3uRm4JKH/1aRfP3njIySjzcMbQde8q85Jb76J4cQYTNzJ3frOPCninsKnbz3pIM1udWADD+xHQu75vON2tzGd45nnH9W+33HFX6FCJsZnRdJ6fcxwd/ZtA6xkF6tJ2FO4v5fOUe3rq4FwWVft5ZvIsh7WIZ2SWR0R8uJSnCSn6ln9tPaUduhY9bBrfllPZxbMqrZFtRFcM6xLOn1EOFTyGz1IPdYqJVtJ2lu0tZm1OBO6Ay7e8sIPiM/2ljPiM7JzB3ayERNhO3DG7Hr9sKuf/0jvRIieSpX7ZwVtdEuiVFsHhXCalRNrJKPXy9eg/Pni/hCWhc/PFyhndKYP62Qh4/swvndEvkvMlLKfUEuLJ/Op+vzEJ7eRSSJOENqBRU+mkd6zjIO+LwI+mNUABut5uxY8cyffr0sIbrcLMqq4w/M0oOqIym61RUVBAZGYncQh4WP27I45lzutGvVTQAGcVu2j03D4CHRnTi+XnbmtM8wSHy3yv6cs20VSRH2owHweFi4R2ncEr7uIMuL/3zhyazxWk14far6K+MrrVt645dvPPm69x7772kpKRgNptRNB2zLFHpVzHLEnd/u47sMm+1YRKaHhSTshQUneXVojKgBcX7viS5rDw4ovNB26+qCtPffx27w8Fzjz10RMXlrYu/IstddtDl053RvDvkHwdd/u233zbiSwGeffZZevfuzQUXXABAeXk5V155JTNmzKjTU5qZmcmkSZO49957wzylDaGoGhklHv75/Xq+W5930LYLoGdKJLNuGsT9P2ygfbyTP3YUs3BncXOb1eI4o1M8v24ram4zWhzuF87DcQCOlibZZyN1ZIv2lEbZzbQ9QEWv6Tqlko+YGEeLEaV3DG1Pj5RIY71tnBPff85nR1EV3ZIj+fd5JzD0zYW4rGbuOrU9p7SPY+wnK3juvG4MaB2Drussyyzlh/V5DGoby18ZJfhVjRiHhYfrGRpryXSId7KjyB22rmh6mAdr0gXdWbg1l683hje0KZE2cg+z8DtQrpm2CsAQpKHeNIDNLONTtHrLto11kFHS+Ak2Q99atP9MdXDtwNY8eVaXgypbH26/CsDrv+8gq8zL5f3SiLJb2FFUxSVv/cFZFW4qfAHM7gBFHg/eQPh5uOnktk1mS7s4B7uKPUTazCS4rJhkidxyL5XVNu5Lv7RIXn7lFaLsFh55+F9H3Nt5KILycNC6dWsyMvaO6GRkZJCYmIjD4TiokYu6MJtkOia4+Pb6k4DgiIL9gZlNUvcHl/Xmh/V5/LBhr9jtmRLJutwKOie4OLtrIm8t2sX6+0/HZTUxZVkmT/2yBYAV95xKnNPKqA+Xsjm/khOSI5h782CSn/yFXqmRrM0JeuVynzyLlCd/4bvrBnLhlGWNsmtE5wTmbW3Yu3swrMutoPUzc5u8Xmg6b7PDIuMJaEgSfHPtQJ6du4XlmQffEavJ2L5pTF+VzRsX9eSub9fVm6+pBemEk9vwwZ+7AZh3y2BGvLfE2NYvPYq/s/aet6QIKxaTTFZ1x/u/V/RleWYZr1zQnS0FVXy7LoeHhndm4c5iTntnMQAvj+7OfT9swCxLdXbEm4ojLUgPhBbtKT0YNE0jPz+fpKQkZPn4iNE86/0l/HLzYAC2FVYRbTfjUzRaxThQNd2YTFHmCZAQYTssNmzKq6BrUgSSJDHmk+X8tr2IP24/ha5JEQdVX2Ovo09RKXEHMMsSAU0nNWrvUKMnoJJZ6iGn3MvA1jHYzSZmb87HE9D4dl0ub13SE4tJ5q5v1vHahT3wqRpZZV4+XpbJpAU7AJh788k8OHMjyzPL6JzgYmth1X5tD7x4/n7jgzVN59av1vD+ZX1QVI0Ve8poH+ck+clfwvIdjsbpin7p9EmLomuii4t6pRrp5d4AszbmM7ZfOhD08GeWeqjwKozukczJbyxkYOuYRj1kXf5SzqpYyITbJ5KQknbINqdG2VA1SHBZqPSr7C7xkBZlJ9phptSjkB5de4hZ1XRUXcdqkinzBoiymYOdVo/C1I/eJysri6effhqr1XrI9h0tqKqKqqq8//77xMTEMHbsWEwmE5mZmTz00EM888wzpKam8sILL9C1a1fGjx+PrutN4iltiLwKH3kVPnYWu1mbU85jP2/m55sGoWg6sgQv/7aDT8b1Jb/CT4+USEyyRInbT26Fjx4pkYZdFV6Fcl/ACFM6XOwp9dAqJriPjGI3bWId7Cn1sja3nKQIG71SI7GZgw9+Xdf5a3cpJ7eNxRtQWbGnjMFtY43wl9BwsqrplHkDWE0ydrNMXqUPt19l6e5Sxk/9+6DsvKBHMuNPbMWWgkoenbUZgJ0PD8dlMzN/ayGndYzHYTFhkiWembOFRJeN+87oSGGlD4fFhMsW9Fv5FQ2rWabCq2A1S/gUjYJKP7IksSGvAh0Y1T250XYpqkaVX8VlNRltpaJqYe1mKPypsNLHpyv2cGHPFNrHOdle5DZinyF4n/+2rZB7vl9vdCAaQ4+USP7RK5WHRnRic0Ele0q9nNI+jhiHBV3XWZMTDG0qcQdoF+cwfmNbCyqDYR3RDpbuLuGaaavY+MAZDe4rdP7qI7vMi1/VaBe3V1+FRlTn3nwyc7cW0jbWwS1D2lHm9uMuKyIxMYl1eZVYTDLjPlvJ6vtOQ9V0csq9JEfaMEkSsizx5epsLvtkBVOv7M8V/dMbfX6amsbqSCFKBS2SlnQdVU1nT6mHtjUajGkrs7isTyqZpV5+217I1QNaI0vw7uIMbjulXZPvXwK6/edXthZW0SnBxbZGiOP6qGu4/UDxBlQW7yrh9q/Xsim/stb2xojSBJeVwip/vfs4sVU0G/Mr6ZLowtyEv4H8/HxuuOEGrFZr2G/rySefpEePHk22n5bI1KlTmTZtWljaxIkTGTlyJHPnzuXTTz/F4/EwePBg7rjjDiwWyxERpYL9ExJpX67OpkdKJB/+tZsNeZU8c05Xhr+3hKlX9uf8amGoqBqqrhvCGCCgqOTm5ZGemtLsberhZGNeBVaTTKfn5zOwdQyL7jyFGauycQdUxvVLx2k1UeY9PHH6R4KDeTauzSmnVzNPJBSitAWIGcHBI67j/imq8hPvshoPq+WZpQxoHQMEPbFVfhWvopLgsvLlmhy+XpPDZX3SuKR3asMVHyCegMotX67hrqHt+X1HEVaTjFZRyMaZn3HB1TeTmp7OCUkRWKsfkGuyyzkhOcKY8FWTmpO8BC0DIUqPDUSbemxwtF7HYyKmVCAQ1E+8KzjsHBIKIUEKwZnjkXYzkdW3+GV90risz6EPo9eFw2Liv1f0A+DEahsyM83s+EXmhOQI0hKdYcNyvdPq77ELQSoQCATHL0ePzBYIBEcdLWWyoUAgEAhaPkKUCgQCgUAgEAiaHSFKBQKBQCAQCATNjhClAoFAIBAIBIJmR4hSgUAgEAgEAkGzI0SpQCAQCAQCgaDZEaJUIBAIBAKBQNDsCFEqEAgEAoFAIGh2xMvzBQLBccdbb73F0qVL8Xq9JCUlcfXVV3PSSSc1t1mHnZkzZzJ79mwyMjIYM2YM48aNA6CkpIQ333yTLVu2UFZWxg8//NDMlgoEguMR4Sk9AnT/+kWkKfchTbmPyoCvuc0RCI57LrzwQj788ENmzJjBXXfdxSuvvEJ5eXlzm3XYiY2NZdy4cQwZMiQsXZZlBgwYwD333NNMlgkEAoEQpUeEf7TrbSxHfvYI0pT7uGnRF81okUBwfNO6dWssFgsQ/EyroigUFRU1s1WHn8GDBzNo0CBcLldYenR0NOeddx4dOnRoJsuOLH5VaW4TACj3e6kIeJvbDIOAplLm9zSpTZmVpei6To67HF3XqQh4qarhnFE1DYBCbxWri7PrrGNdSQ6aHsznVxV2VRSzvDCTpQW7mZ+9lYCmouv6QdkXqjdEfcde5K0ylgOaSlXAZ5TVdZ1vM9bVuw9FU42/f+ZnhO37j9wdzM/eiq7reJUAxT4360tyKfRWUeJzA+BVAiiaapwrrcaxLsrbyfbyQj7fvrLWMYXqO5oQw/dHAKfZwqfDrmBbeSFPrZoDwOQtfzF5y1881Hs4z6+Zz7qL7uP3vB3My97KxO6ncnJSW8ySjI6OT1VZUbSHE+NbYZVNmORgX0LTNWTp0PsViqYS0DQcZktYulvxYzeZ8SgBXBZbWH6TJBvfXN8fuq7vN2+xz02s1YEkSQQ0FUnX8agBctzlpLqikJDQ0dF18KoB1pbkkuaMorUrhoqAj4qAj3RXNLquU+b3EmNzUOrzEGW1hZ0jXdepUvyouka01YFHCVCp+JCRsMgmoqz2WufGr6lYZRNLCzIZktyuTvu9SgCbyYwkSZT43PhUhWK/m1RHFA6TBb+m4jBbMFfb4lMVMqpKSLRHEGdz1tqnWTYZ9gKNPtcHgk9VKKxuaF1mKzE2R9g+NV1HliQ0XTd+cwBVAV/Y7+FIUO73YjdZMNewY3dliXHuygNeNB0kKZjXowaMfDFWB/E2J9srioi3uSjyBY/5p0+msWrhnyiBAAMGDKBdu3ZH9JgEB06oLVE1jQJvJVFWO7/n7qA84GVwYjtaR8QAcNIPr7O9ooibugziP2t/BeDlgaO4b9mP9dY9JKkdK4r2MDSpPfNytoZte6DXGWwszeevwt30j0unQvHx08gb2FZeyJ8FGZQFvLSLiGN4aicUTWN7RSH3LP2e01M6clGbnry07jce7j2CHRVFjGrdnc93rGRNcQ7vbV4CgEmSUauFxAO9zuCHzA2U+DwsGXUH83O2sbIoizcGXdRgO1Dic7Mwbyfntz4BWZJZW5yDX1PR0ekZk8LY3z5jecFusr0Vh3IJODu9K7OzNht/AVIckQQ0laJmFkBfnnE1AxJa0coZg1cNf265FT8Ok4X7l/3I9V0GsrIoi1VF2byyfkEzWty0jP99ap3p6rUvkuOuwKMGiLTYSHZEHmHLGo8QpUcAr6rQIyaF8R1P5J4ew0ic9iSB6p7T82vmA9Dz25eN/F9lrG0WOwUHx/TTxzP2t894/sTzeGjFzMOyj1irgxcHjuLStr0N8XgwPPX3Lzy56pcmseniNj2JttqZcurlB1V+X8G9b2cnoKlYZBOKprKlvKDOOgp9VXWm16TU76HU7wEwBCnA+Vdfwbnjx7Jr0xbspZ7DIvyPNzxKAOenDxnr4zr04/SUjtzU9WRUTaPY7+bP/AyeWT2XP0fdiSzJaLqGputGR6wmmq6R76lkxq7VTPzruwO2JyRIgQYFKcDi/F0AtQTpvvXMytoEQPTnj+53/yuLspi0/ncAfsjcUG8+tYa3rua+2n3xb2N5Y2k+83K2Mv+cW+gclUArVwwQvI82luXx8roFTNm6DACLbDKeMU1NSIiG/gLkeg5N6DYVl/76Sdj6xO6n8vqGPxie2on5OduM9GNJiDYG08f/Clt3X/V8LSdUS6FFi9Kvd63l423LDqiMruv4/X6sVmuLechUBnzc1GUQANFWB/5r/oM05T4AbCYzvhYylCQ4OMb+9hnAYROkACV+Dzct+oKbFn1B95hkdlYUs/UfD5Jod7G9IjjsbJVNfL5jJX/m7+br4dcAkOMpJ8JiI9EeAdBkghTgm93B4aqXB44m3h4+HFxWxxBYvqcSt+qnMuDDqyq8sGY++d5KACItdmPYLMbqwKsqeGt4O+siyR7Bg72Hh6WlOqLI8ZTv96EcZbETZ3Oi6BqRvXvz8aQ30WNcdO/bBwkwSRKyJIf9NUkysiTts21vuozU6DYn5+NbCZRkNSpvXVhi00m99t2DLn+42F5RGLY+dcffTN3xNxMWf1krb+hBeUvXwby3eQn5VzzJnqoyOkTGcdbsD1hamHlEbD5aCInl4T+/B4BZklH2GXoOcbgE6YFglU1EW+0UeKto44phd1XpEbfh9Q1/AIQJ0kOhS1Si0UF2ma1UKf5Gl3WZrZzbqhtfnHE1ExZ9wdj2fRme2omt5YV0/fo/ALx20oXcvfTAO18HSksVpNDCRekl7XpxSbteB1RG0zTy8/NJSkpClltuyOyGi+8nz1PB6amd2FFRRMcvnwdgfMf+3NJ1MENnvg2Aft3LYbEykiSh6zp7qsowyRLp0585oP3e3f1UXtvwBwMTWmOVTTzUeziJ9ghOSmzDysI99IpLxSRJFPs8bC0vwCqbOTGhFcU+N5qus6Y4m+FpnVlVlEWJ30P/+HT8mkrStCdp44rBYbawuSzcq5XmjGLPmMcIaCqL83dxemqnOm0r9Xko9FXRKSqh1nWsGQJw6+Kv6BmbwrDkDlz1xzSKfW7Ob3UCM/ds5Jvh1+LXVDpGxZM07ckGz8VFbXry7e7644AOli2XPECczUmV4qdNRCweJYDNZMJb3fl4YPlPTN7yF15VoXjc0/yrer2xbCjNA6DVjPqvvaOGtwrguf7ncku3wWFpvWJT+eWsm0id/nRY+j09hvFqtXenMbT74t9MPuUyZmVt4sYug6hS/Pzjx7c43+2m1O9B9lWRU1bbk7KvoDwUBiS0pjLgI8JiI90VDUCep4LMqjISbE4KfVW4zFairXay3eV0iU40yqY4IrFLJqTSKjpHJaDpuuG9U3Udtcaypmv4NA1ND1SvB7er6Oi6Tl1RbRLUFrNjX8RRQ8hKgIRE8CdevW6k7/UmSzW2+VUluCbVTA+VPjwhH6FjDB13rrsCb1kBPk3BqypsLM03Qj4aS2gIe3/366Ew+ZTLuHHRF+y67GECmsbuyhLG/T6VPE8FGZc9QpuIWHp/+wprS3K4t8cwXjnpAn7N2camsnwmb/mLlUVZ6Ne9zKg5H/LdiOsw/zcoqLtGJ9Zq70J8ccZVnN+qO9vKC+n93Svc2m0w724KHusz/c7hsb9/5qeRN3D+3A8BKBr3NPFTH+fDU8Zww6IZ+z2m+gRpYwmFDPSKTWVtSY6R/umwK7jq92l0iUpkZFpnUh1RPPb3z1SMf47Izx7h9UEX1um19l79AoXeKop9bsoDXiSksHAnXdfxV4vlIT+9yTP9zmFEWmdsJjO7KooZ89unLB09kYqAl7v/+p6Pti7lx5HXM2ruRwB8O/xahqV0IG7q4+y67GHe2riIlwaOZktZgSHqDif/7HEaL580mtNnvcMZKZ14ot9ZhnOppudxR0URyfYIXBYblQEfX+5aw/aKIp7pf45R1wenXGYsd4lOZEhSOwKayp3dT2Fl0R7OTu/K1xlrw0ZNXxxwPv9a/hMAz5wwnMc2BkdZQ0J57UX/pNe3r4TZ/Oagi7jzr2+N9dcHXciHW5Y27YlpYiS9EdHBbrebsWPHMn36dJxO5/6yNytHiyjdl8hPH+F/p1/J+a27H3BZacp9qNe+yMw9m1icv4sz07owZesyJnY/lb5xaZT4PXy0dSmpjiiu6nQiG0pz6R6TchiOIkjvb19hZFpnXhk4ml9ztjE8rfMB13Gg19GjBFB1jYgaMUSVAR+/5W5nSFI77CYzdpOZ/+1YRbHPzQVtetDaFUOZ30u7L5+jZFxQ4O2uKmFtSS6jqxvCO08YypsbFwKQaHdR4N3/cLFyzYthMZj1cebs95lz9s3kuMt5b/MS+salccn8/+633KEyvmN/tpcX0TU6MWzovcTn5pL5/+XXc28F4Orfp9HKFU1FwMel7Xrz1a61xrnYH64KH2etrWDCXbeTkNo0v7UEmwu72cKeGh6XDpHx5HjK6VHH71mvFo0mWWZbeSGdohJQNI3skiJ2rF3PoEGDsFqtLFmyhEmTJvHyyy8flok+eg3xWlPcqrqOVkPMhmKmg0ug69VpVKfVWN8rgOtOr9moSzXWVVVFUzW+/u9nREZHceZFF2Aym5BlmYA/QEVZGc/efT//mfI+kgTm6slguqYjyUGJLAP52Tm8+/qb2M87mYjkBOP+sslmzkrvSrTVTrI9gojPHmHZ6IkUeqs4d85kbus2hHc2LW7S8/vigPO584ShyJKE7ZMHgaDAg6A3MRQn/s7GRdx2wikN1jVrz0bSnNH0iUtrMN/c7C2MTOsCQFZVGbE2BzsritleUUS6M5qFeTuZ2ONUIHhNlhdmMiChNaV+D3FTH6d8/LMsytvFOa26hdUbiteuCvj4KmMtT62awxXt+/LcmnkHfmKq+XTYFXyxcw3fZ65nZGIHvjv7BmwmS1gble0u4/k18xmR2pnzWnXjg81/MjCxNYMS29ZZ5+ayfH7J2sLo1t2xmkzkuCs4MaHVQdt4IPhVBavJHLb+weY/KfV7eezvnw+orhWj76ZffDpVip+AphJbHaOu6Rqmj/+Feu2L9PjmZZaPnoi9+pzVDDvK81Tw6Mqf+b8aIrOpUDQVy38fYPIpl/HJthUsOO82AHKqypAqPMazcX1JLj2/fRn9updxK340XcdptqDqOhbZxJsbFnLXX9+ydNRdDExs0+R2NpbG6kghSo8Bpu34mys69GtuMww0Xav2+hy8p6YlXMcl+bs4KaENhb4q/m/zX9zbcxg57nL+uewHvtu9nhWj72Zh/k4m/vUd7w3+B5e170P81MfRr3t5/5XXg6KpVAb8RtyoVwnwyfYVeJRAkw3rVI5/jvKAD4tsImGfYff98W3GOn7P27FfT2pTiNIoi53ygJdOkQlEWKxhk78qFR8usw35IH5jbrebZ599lh07dqDrOqmpqYwZM6bWa5KORaZOncq0adPC0iZOnMjIkSMZPXp0WHpSUhIffvghuq6jKApms9m4pzMzM5k0aRL33nsvrVu3PihbQl6mAyX38idwma1YZBMShAkUgAeX/8QLA84/qLpbKvOyt/LxtmW4zFbe3/xnvfkKr3gKkyRzyfyPOSEmmbcHX2Jsawlt6uHmhoUzyHKX8c3wa8Nim0O8Oegiesel8tn2lWHeyrrI81S0iAlB+8be1zWKWOirMkK06kKacl+zx5EKUXoM33jHA0fTdfx+93ouaNMDgGUFuw9bb/Sltb/yr+U/8dmwcdz+59eU+WvHbb4ycDTLCjP5385V9dZzKKJ5Xx5f+TPPrJ5bK70xorR9RBw7K4vrrfvE+FZke8pJdUQ2yVsmBAfH4RKlpT4PsVMfq3f7ifGtWFG0h1/PuYW/CnbzQBOGexzNFHqrSLC7+GjLUubnbOOjoWOYuWcTF7Tp3uB9cjS1qU1FSN6U+b1IEkSYbY0axWrJHK3XsbE6skXHlAoERwMhQQoc1uGR+3udwf29zgDgyo79eXvjIpLsEVzWvg8Ao+d+yL09TwOCs06nbF3GsJQOjEjtZMSNfj5sXJPa9HT/c7i12xDSasSlJthcXJbUi9was4g7RMYTbbGDBH8XZdErNhWbyUyczYlfU7GZ6m6K0p3RTWqvoOUQY3OQPfZxFE0j1RmJomksL9pDZlUpj678meUX3E2Z30O01VFvHPrxSGh04/ouJ3F9l+BXyC5q27M5TWqxhDpRh/LGEsGRRYhSgeAo5fZ94uN+GHmDsXxyUltOTtobD6Zf9zJvbljIuI79m9yOVGcU3qtfCBOWmZmZTJq7ihirg+6RiTisNuMB4TRbjfeNSpJUryAVHPukOqOMZbNsYmhyewAjHCnaKsSEQHA8IZ4GAsFxwp3dhx62uhvydFr2ef9k95jkw2aHQCAQCI5ejp6ABIFAIBAIBALBMYsQpQKBQCAQCASCZkeIUoFAIBAIBAJBsyNEqUAgEAgEAoGg2RGiVCAQCAQCgUDQ7AhRKhAIBAKBQCBodsQroQQCwXHLxnXreODhh7n84lFcdvFoJNkEsglJNlf/rbEuHdqnc1sCM2fOZPbs2WRkZDBmzBjGjQt+TGHZsmXMmDGD3bt3Y7fbOfXUU7n22msxm8UjQiAQHDlEiyMQCI4qjC8j6zpUi8SQWNR1HUmS0FUFJBmp+iX9uqYCUjC/roEko3gqmPzRZDq1b4tscyHbI0FT0TUVXfEFy1SvoynoulbLFkky1RCv5hoi1oRk2kfYtoBPpcbGxjJu3DgWLFgQlu52u7niiivo0aMHXq+Xf//733z99deMGTOmmSwVCATHI0KUCgQHiCF8FD+S2Rq+TdPCPGqqpxyTY+9Xa3QlgGS2GPXofk8wv8VOoGg3ktmKe9MCIvuOwrt7Ndakjijl+Xh2LSd68JXIFltYHaF6UAPoqoJkdeDbsw57614AlP7xMa4eZxIoysDR8WQKv38WpTwPR7sB2Nr2xbdnHd6Mv9E8ZehKgIje5+DP307s6RNQyvPx520lULQbR8dBIMl4dy4n99M7SLjoCSL7nE/We1fS8T+bKV86A2enIZhj05FkmfKlXwZFo66j6zq+PWuxJLRDstjRqkoIFGfWPrGyCXNEApLVEdyuqU196cL4ZeFyOqbF4/GrSCYzJnvEAZXXdR10rYZ4VcJEreavArWGqEUPKy8hg8mMZLIgVf/FZKletxiCuikZPHgwAMuXLw9LP+2004xlm83GGWecwdKlS5t8/wJBS0PXNEAPdh5rpqsK6HpYW3vAdVc/KwA0vxfZaq9urxV01Q+SCdlqr7e86ilHrSpBtjrQfG6sie3qzBco2o0W8KH7PegBL7ZWPdAVP7Iz5qgb3RGiVFALXdPQfFXI9ogwD1TIwyRJErqu48/bhikiDnNEPFWbFiDbIij47mliT7uRotmvknLl63h3rybqxIvQVYXy5V/h3bWChFEPITujCRRmYI5NQ1f8eHcsxRyThnf3atybF5Bw8dNU/fACmxZ+AIAtrTu+7A0ARJ18BeV/TiPtxilUbZgXFFa7VxF3zj9xb/4dzVtB1EljsMSm48tahy97E7FnTMCbuYbKNbPQ3KX487YB0PahBVjiW+PLXEvZn1Mp/2v6AZ2r+HPvo2jWy0149usn58Mb9p+pkZTUk162+FMACr5+vMHyhd8+ReG3TwGw8dra4qlAdeGXzsKXvQFVzQfAn7ulYaM0FaU8r+E8TURFlZuZv/3Jc/feyH+/m3dQdUiSBJKp1sOsseiaFhSyagA99C/gRVcDwbQanllJkmuL1up/mMxN/uBZv349bdq0adI6jzR6dadIkmV0TatT5O+bHuoo6pqCyRGFruto3gokk7VO8eDL3Yo1qQOauwzVU44lNg2lsghLTCr+vO3kTr2buLPuJlC4i+jB45CtDnRVCbavVidKWS6V6+cQNfBSqtb9ArpORJ/zkCx23FsWBTs3njI8O1eQdOlzIEkE8ndgSWwHuo5310rKl39F0mX/Bk3Dn7+Nop9fJWXcJHKn3kPC+Q9ijk5G9ZSjVBbj+eMzfKddjT2lC6UL/0vORzcSNfAyLIntiB48HsliQw94yf30TnRdQ/d7SLzkaUp//xBHpyHoviqiBo1BKc9H87nx7lyG6ikjotc5WJM7IVscFP3yGmWLPkG2R6L53SjFe5BsLnRfFQCxI27H1WMkBV89ii9rfTBt5B04OpxE5eqfarXBSZf+m/wvHw4KLNmE7IgiULCz3uvu7HIqKde8g1KWC5pG1fo5ODoNwRKbjikqkeyPbkQtzTWeJ0crxQeQ15rSBX/uFpxdh9H2od9atFCVdGMsrH7cbjdjx45l+vTpOJ3OI2EXAOXLv6b0j48PsJSOz+fHZrMCLePE62qAtBs/whKTesh1aQEfssVW73Zf1gYs8W2QrE58WevZ8WhvUq99H3NcK4pmvYx746+HbIPg2MHZfThKaQ7+7I1NWm++6uIL6SzuuX0C6ckJdeYp/PEF1PL8xlVoMgc7RaEheFXBmtqVlPFvACDX8Y10XdOC3lbZhFKehzkykUDJHsxRKbz7fx/SoUMHzj33XF599VXS0tIYO3bswR7uAZPx7XIC5e4DKKEbIgtdwxJhJX1kh2oBG+6FlWRzuGCtKWBlE5Ik8fbbbxtD+fuyaNEi3n//fd544w1iYmLQdR1FUTCb94rfzMxMJk2axL333kvr1q0P9XTUS/nyr3F0Gkz+9AeIP+dedDVAzie30faB+eh+N9kf3oA1uRMl899FV/wkj3uVvKn3HDZ7BIJjgW7/52nQQ3s4aKyObNGe0qgBlxA14JIDKqNpGvn5+SQlJSEfhuGvgyH/q8dQy/IOWZQWznyJ/On/OuByOR/ffEj7bQymqKTGCQzZFDYsm37LVPKm3YtSlkvyFZPQNZWin14A2Yzj4heJjoqkavWPOLudTtHMF2n30AKU8jx8WRsoX/YlKeNepXzF1yjV+y789imsKV1QKwpxdDqZytUzjX1ZU7ui+apQivcAYG8/EO/OZcHltv3xZqwESQZdw97hJLw7luLofApKSRaBwl2NPhdd3yvDn7sFW2o3JKsDtbyAQGk2jnb9g56YgA/P9j9xdByEbHXgy9mMOTqFnCkTsKZ2o3LVj0T2v5CyxZ+Scs277H7xTOxt+uLdvarRNjRE9//utx8ahq4qwfhIoHzFt/hzNgUFgqbi3rSA3a+c26h6LPFt0HxV6AEvmq+KhFEPYo5OQbZHBj2E6CjFezBFxKGrKraUziiVRZiqPSQHiiTLUN0GhO49a0I7tm/fztatW7nlllsOuM6mou1FAw5LvbquB0MIanpf/W40VQkuawoAqqcM1aLjL8wwQgckk4W1G7fw7jvv8PjjjxMTE9Oktml+DyFHgeYtx71lIc4uQ0HT8O5Zi2fnMlACWBLakj35ulrlQ158gM23RNXaDghBeoiYo1PAZEGpEV5jTe2KP2fzQdXn6Hgynp3LaHX7DAp/fAGlODPoway5z5hUlNIcAOxt+hIoy0EtOzIjJvsjsv+FVKz8zliPGjTW8OSmXvseOR/fQpt/zsKz/S9kZzTlf00n/vwHqFr3CylXvg6ShK6pVKz4Ble30wEwx6SgKwEq/v6eiD7nofmq2PFYX2JPn0DcyDvw7FqBLb0HedPuDfMau3qMpGr9XADsw27F+/u7B31cR1qQHggt2lN6MLREUVrw3bO4ep6Js+OgRpep2vQ7VRvmYUs7gax3rzjofVtTu+HP2WSsnzBFxZu5Gn/2Jnw5G4kZdgPWhLZoAV/w4dSIh3+gNAdzVPJ+Y950JYCu+pFtrgO2uzmvY804oLoIxQZpfi8F3z5J8pgXAFAri8l842LaPbyg3rJNYp/iR1cDSCYrmq8Skyt2v2U0v5et97ah4wsbMUfEHx67NDU4XGoyBz1pr7zCvf/8JykpKWFetubmu+++47PPPsNuDzbMbrcbk8nEkCFDuPvuu5vXuCPE22+/TUxMDFeMvaw6vi3A5s2beO6l17jvthvo0aV9WPiAroNssRniNSu3gFffeIu7b72Rtt1612o3/IUZACil2Wi+KjInnY+u+I/oMQJYEtqRctWbmGPTMUck4C/YiTUxGNssO6LRFZ8xTO/dvQq1shhJkpCsDuxt+6P7PcjO6Fq/XdVdZsQbVm1agOuE4WjuUkp+/5CE8x8wPNIASnkBpQs/Jvb0CaCpqJ4yLLHpwXvYYsez7U8y/jOcEz70AaD53PgLdrDjkV60ffBXMl44o9HHm3rdB+RMmYCj8ymk3/wp/vztRPQYCUDeV4+jdD2H+KQUSma+SNKY/4Cu1Yo71DUNpTQnGN++dRFRJ15kbNN8VUgW+0GHrEAwTlLzVGCJS683z75hFaq7DNkeGQzFUBWQJDRvJSZnNJq3kt2vjiLuzIkEinfjaDeAXc8NBSDhwscp+fU9LPFtsCZ1xN62P0rJHmKG3Ygt7QSUslzMsWmg6yil2Vji2wRjOB3RwbANs7XOURjVU4FsczVZ3LdaWRw8vn3iV0PhI54dSzHZo3B0PCns2VgdXIckm1BKcwkUZ+LoMNCY76D5PSjl+SilOTjaDzCcC81BY3WkEKVHgMKZL+HoMAhXt2G1tum6ji9zLbb07hR89zSF3z3TYF0dn9+ANbUbsHfG8a7/jKDtv+ZS9NN/iDvzLmRb3deopserpdMSr6Og8dQc3m1potTr9eLxeIz1Dz74gOTkZC699FIiIg5sstPRhqqqqKrK+++/T0xMDGPHjsVkMpGZmcmjjz7KnXfeyaBB4Z1nXdcJ+P2YZT0oCNQAmZm7mfTqa1zbHRIpqX67AaCpVK6Z1WT2muNak3z5S3i2/0XSpf8O8/DoqoJaUYg5JqXJ9tfS0XxuKtfMJLLfhdVvdqgRD6upDYpF0aYeGxyt1/GYGL4/VqhaN+eght1DmCIT6fJ6dr2Cst0DwYkaCaMebLCeo0WQCgSHE7vdbnhJAaxWKw6H45gXpADTp09n2rRpxvqMGTOYOHEi69ato6Kigpdf3jtpr3v37jz1VHAymyTLSGYzsiV43kyuSkzOGBIvCY8pVSoK2XJHYq39xpx2I5F9R1O1cT4JFzyKyRVLoCgTc1SiMZKieStRq4qRbRH4cjfj7DTYKB89qHa8r2QyH1eCFEC2OYkaeGmd2w7FeykQtBSESjkCxJx6LVXr5zQqb7vHFgMSzk4nA9XDJWabEJQCwWHinnuOnzjEcePG1Tm5aeTIkU0SumCOTGgwZjmy/wXG8r6vt5HtEcjVr+WqKUgFAsHxg1A6R4DoweOIHjyOQEk2nu1/GpO39he7CBxUPKZAIBAIBALB0YYQpUcQS2walhpvE2gpMXYCgUAgEAgEzc3REyUrEAgEAoFAIDhmEaJUIBAIBAKBQNDsCFEqEAgEAoFAIGh2hCgVCAQCgUAgEDQ7QpQKBAKBQCAQCJodIUoFAoFAIBAIBM2OeCWUQCA47njooYfYvHkzJlPwKzg1v150LDNz5kxmz55NRkYGY8aMMV6kv2nTJt5++20KCgowm82ceOKJ3HLLLTgctb/7LRAIBIcLIUoFAsFxyZ133skZZ5zR3GYcUWJjYxk3bhwLFiwIS09NTeXJJ58kPj4er9fL22+/zbRp07j++usPmy2K24fi9lO+LZeoTsmYnTYUjx97fCT+0iqsMS40RUWSpbBvvAsEgmMXIUoF6Hrws4DH28v8lSofZpetuc0QCI4YgwcHP9+5fPnysPTo6OiwdVmWyc3NPej9FK/aRVTnVDJ/Wkn62X0AyF+8BV0LtjWBMjcl6zIPqu7obmmUbckBrfbnTFNH9ETXNErWZuIrrCBpaFcsLjuVGQW0v3wIu79ZRlTXVCLbJ+EvcyPJMpYIG7LVjGw1o1T5KN+WS1yftshmE76SKsq35JAwoAOqN4DZZUPXdPylVdjiIhq0U9d1AhVedEXFEmFHttZ+3O7b9uqaBpKEUunD7LQimZpGjIfOu67rKB4/Vpe9Ue2+rgePVZJlrNHOBvfhLSjHnhh14LY14suGguMHIUqPMcq35rJ1ym/E9mpNydq9jX5E+0T8pW78JVXE9m5DyZrdpJzRnUC5h6IVOzHZLXS+/nQ2vTMHAEuknfgBHShcup2EgR2xxjgx2SwEKr24s0so/nsXAB2vOpXtn/5Bh3GnYLJZMLtsKG4flgg7WkClfHse5Vty0AIK8f3bI0kSstWM2WUjd8FG4vq2pXT9HvzlHrx5ZVhjnMg2C9YYJ+Wbc8iqtl8yy+iKhiMlGntiFCVrM+l680gyvlmKN78cZ1osHcadgmSSkUwy2z/7g1bn9cOTV4Yl0o6/uIqcX9eTMuwE9sxaZZwXk8NKyrBu5C3egiRLtD6/P47UGLZNWUDiyZ0oXpVBXJ+2SBYTztRYVK+fvEWbaXvxQLz55dgSIrFGOZAtZnRdp2p3ERXb87AnRbFj6iLi+rWjdMMe4vq0o2xzFpIs4y+pqnXdYnq0wldcierxY41xYY1xUrwqg14PXIA7q5iozinIluDtqlT5gg9Ri+mAfx/ewnI2vDGbiLYJ2OIiKFy2HYBW5/fD7LRhi49g65QFOFOiiR/QgarMItJG9sJf5kZXNSwRdqwxLiq25xHZMZnSDXuI7dn6gO1oLKEHlq5pBMo9gIQleu+Qsq5qqL5A8CHusmFyWMO2yWYZQg88HbSAgqao6KrG5MmTmTx5Mu3bt+eGG26gffv2h+04Wiq6roMe/FuQl89d90zE7XZjt9l56P4HCFR60TUdTVXRJD/oQLVQCVR62f3dcnxRO436NEWlYlseACmnd6d4VQbFqzKa1OayTdn1bsuZty5sPX/hZmN51ZNfAlC8ev/2ZHy1NGx993fL68l5ZHCkxhDdJZWKnfm0uWAApRuzMLtslKzOwNk6nortebS/7GRK1mVSlVlE6vAebH5/Xp111X/2Gkd0tzRie7Vmz8xVKFU+INgZqHnuo7ul4c0vx1dcSfLQrsT2asPWjxegevwkD+tG3u+batXbfeK5lG3ORjabsEQ5sMVFYImwU7RqF7ZYF/5SN4kndwJA9QZQqnzYk4OdKU9OKSabGUu0M/gMaKTI1VUNXdeRzbXbUl3TkWSpVpqmqJjq6GAImgZJD3WXGsDtdjN27FimT5+O09lwb6kpqdxVQNnWnAMqo+s67io3TpezxfS+JCSST+2KyW7df+YGUKp8ZM782xCEguObuL5tiWibiOpXKF2fiSXSQen6PQBEdkgi/sQO7PriTwDaXToIW1wEmz+o+0HVFHS8+lQsLjvO9FgyM3Yz6fXX+Oe9/yQ5KQmLzYo3vxzZZkZXNHRNwxYfib+4El3XscVHovkCKB4/stmEpqionsBhs3Xbru20SklHlmV+/m0OPy+Yy7vvvntE27dDoaaYRNPD/uqaDvX9reaD/35ITHQ0Yy66NNhOylLwASxJlFWUM/fXeQwdcgppaakgSaiaitlsCeaTYE9WFpMmTeLee++ldeu9nRJPbikb3vj5gI4lom0ClRmFAMT3b0/Ryp0N5rdEOQiUe3C1TaCqupxAsC9Jg7sgmWVMNjPRJ6Tj3lOMLT4y2PH+eAHevLKw/K7W8QCY7BaQJMq3BLVHRLtEKncVhOXtdM0wZIuJjK+XEdEhiaLlO2g9qj+JgzqhKSq+okpMdgu6rmONdCBZTIYe0RQ16OWXMBwNB4KmaeTn55OUlIR8FIW1NFZHtmi5b0uIJNZ2YCZqmo5eXExMXByy3DJEacGSrfhL3ThSDl6U6prG6ue+aUKrDh5X63iqMosOqEzyqd3I+2MT1hgn/lJ32Dazy2b0uKM6p1C+NThsGDeqJxV/7CBQ5sYS5UDzK6jewydUDjedrzsd2WJi8wfzsCdG4S0oJ/3s3mTNXnNQ9TXkharYkU/FjnxjfdeXf4Vtl61morumhnnTD5Xtn/xhLBd4y/CXuvHklxGQXQQIeoe1gGrk2fm/xcZ1PxjMLhspp3U/qLKd2nU0li88ZxQLli5i8+bN9OvX76DtaQxBManjKX0dTS0MJVb/Df4vzE1Qr88gHqv11jAxWfOvJMnVAlIyhCSSZDwYTXYLZqcNe0JkrZoToxwMPPkkXn3rdSZNmoSu62hKcLRifx39muEwJqeVyHaJuNomkjS4M7LZhOL24c4uIaJNApLZhDu7GFeroBgIeabaXDSAQJkba4xr73ERbAM9OaU40+PCzyd7h6Azf1iJNdZJ8tBuddqnKcHfX8ibVrYlB5PdQkSbBHRVQ1M1JAm8hRUobj9RHZOBoENAsphQKjzouo7qCQSH8BOi2PXln6i+AP7iKqxxLmSTCW9BORD0/pmdNtw5JWT9vJpW5/Yl59f1VO4qIGlIUDSVrs8iom088f07kPnTSlyt4ij6exfpZ/Vmz8xVDZ7vQ6XXAxeQu2AjZpeNnHnrsMa6ao/kVP+OItolUrE9L2yTIy0WpdJbPYoRpGYbvi+2hEh8hRVNfhx1kb9ki7GcPXddAzmD1PdM21eQAmz77+/Gsq+4EoDMH1eS+ePKAzWzFm0uHkhVRiGe/DLce4rpcuNwnOmxIEn4S6pQ/QF0STN++7qmsfGtX+h83WlYIoOjSFpARbaY8JVUIlvMWCLsh2zXkaJFi1JLhP2AT6amaVSa/DiTYlpML8IS5UBTtYMun/v7RrJ+Xm2stzqvL7E9W2OOsFOxIx+l0ouvuBJbXATx/fcOQeq6boi4nPnraH1+/7B6D2csT/n2PKzRDkxWC5ao4I3S6ty+DZbR/IrRowz1Btue3L3B6+grqUIyyVijGp4lHBqmUT0BZKsJJClsCMaTWwqShDXWZaSvePh/9H9mDMgSu778i8j2iWTNXmMIqbSRPets7E789+UN2rLv9thebchdsJGU07vj3lNEZUYhMSeks+XDX+l225lGSEVT0eexSzCHhrmvCN8WqPBUd6CiyfxpFa1H9TOGthS3j6KVuyhdn0llRiGtR/cn84eDb4QbLSglsMVHECj3YIuPRKnwEKj0YYm0o3j86IqGq1UcgUovJpvFCGsI/f5Vjx/JJAdDSlQNXdWM41ervbOSJKH6FRTv3iFqdNCD/wt6Gmul1/BW7ne8ae+xSJKESZqA2VqXmJT2EZOSUeZIoqoqOTkHNkoFYIl0NPj7NzttRHVKMdZDghQwxKdsNmGLry2WJVkOE6RQ+7y0Hh3exu3LvsO00V1S99ZlkjFVx3A6U2PD7a4W26YadoU8a93vPKfBfQJER6Ya+4rqnGJM5AJodU5fI1+orrYXnwRgiOu8hZuI7pJGoMLDlg9/5cR/X447u4SNb82m41WnolT5yPh6Kf2eugxPXilVmUUkDe4CBH/jSBKYJPLz80mMjcdktRjnu80FJwKQNqInAFV7iijdkEXq8B5IslxrCLupcGeXoKuacR5ronr9lG7MJqpjMt6iCiLaJeLJLsES5QjGxEY5UH0KVZlFuLOKUap8FC7fgT0pCm9++WGx90ix+5tlYetbJs+vM9++YRhrnv+u3jq733UO2XPXUrohix73nl9nZ7Sl0KKH7w+Glujazvl1PZEdk4lok9DoMkqVj4Jl2wDI/mWtkb4/wXOs0BKuo65q9U40qCve6HCj+gJsfn8uHa8cStnmHKNX3vP+0RT/vYvsuWtpd+kg4vq2pXjNbrwF5bjS43Cmx1G4bDs589fT6ty+JJ9atxfpYFHcPkrWZobF3RV4y/i+eAV333EX6SlpAJhsZlSfAoAtIQKT1YI7u6Teeh0p0Xhyy3CmxdQ7+zokCA/kWlRWVrJ161Z69uyJrmr88N33fPvDd7z5yuu4XK69IrBa8ISEYdDTyN4Y1X3WW0q4UEOoqoqqqrz//vvExMQwduxYTCYTK1asIC0tjfT0dIqLi3n11VdxuVw89NBDwckxioLZbN7rkczMrHP4XnBkCTkjGkNLaFOPFLqmU7WnCGu0E11RKd2YTcXOfFLP6I4tLoKKHflhsfC6qlG4bDvR3dLw5JUR1TmVsi3ZZM1aTavz+iKZTQTK3cT3a4/mV/AWVeJMjcGTX8b2T/6gx73nUbRyF/H92qH6AuiajslmJlDhZctHv5I2oicmuxWzy4YtLoI1//6WqC6paL4Aqi+AbDVTtae4zsl7h4P+z4xpskl0jeWYGL4/VpBMMvp+PKWBcg+SxYTJbqF41S52fRE+3NrnkYvFTPEjTEM37ZEWpAAmm4Xud50LQNKQSGJ7tja80KnDe5A6vIeRN75vu7CyaSN7kTay12Gxy+y0kTioE4mDOqF6AwTKPRR4Spn96kZs8RFYkiKwWCx1ijZnWmxwooFJrlNgulrF1SpTE0MwHgCqqvLJJ5+QlZWFyWSiQ4cOPPHkk8Qk1/bYHGtMnz6dadOmGeszZsxg4sSJqKrKBx98QGlpKS6XixNPPJFrr722+QwVNIrGCtLjDUmWwpxAyUO7kjy0q7G+7+RMySSTeHJnAMODHdMtnZhu6bXqlq1mnKkxADiSoul53ygAEgZ0AILtYQhbXAS97htdq46GnEuqLzi6abJZ8OSV4ckro3R9JqkjepK/eAuFS7eTestQymdtIqpDMvbkaAqXb6fjuKGUbc5m5/Ql9Z8YMJ4ZLRUhSo8AkkmuFUcJe71tWkBhzQv1u977Pzf2qPDCCI4sLbFxMdktwYkCmWUgSchmU4O/XUmWkNjHI3mYiY6O5tVXXz38O2qBjBs3znhh/r6cffbZR9gagUCwLyabxVh2JEfjSI4mrncbANpeNJDWF5xIfn4+XW4abni8Q9vj+rQlqksqZocVXdNRPX7MLhtlW3KI6pR8VLzvV4jSI0DFjnzKNmax64s/sSdGIZllPDmlDZbpdvtZuNIb9hIJBAKBQCAQhAjFy0uyZIyu1oydbukIUXoEaHvxQNZsDL5xMzQrsz6Sh3bF2SpeCFKBQCAQCATHFUKUHgEsEXb6PzuWlY9Or7Wt532jgp/TCyhhbnuBQCAQCASC4wkhSo8Qkixx4r8vr/c1TEKQCgQCgUAgOJ5p+VGvxxhiwpJAIBAIBAJBbYQoFQgEAoFAIBA0O0KUCgQCgUAgEAiaHSFKBQKBQCAQCATNjhClRxBNLUDX3OhaZXObIhAIBAKBQNCiELPvjwCaWoCm5FNVeJeRFpk8Ddmc1IxWCVoKuq4iSaYmrVNTC5Hk+FoT63TdD9T9yc/wsnFI0t4+q6aVI0lOwISmZmEyt2pSe5sCXQ9+yrem3Q3x1Vdf8eOPP1JVVUVqairPP/98g99kPhaYOXMms2fPJiMjgzFjxtT5dacnnniCVatW8d139X9lTiAQtDw0rRxZjmpuMw4JIUqPABW5VwB6eFreFUSnz2segwR1omtVgIYkR9a5XVMLkE2JDZeX7ATcPyPJ0VgcQ/FX/YQayMAWeTmyae8HEVQlC9W3CrOtPxV543HEPoJsTkWW4/BVfoE9+nYC7pkgOTBZOiLJTlT/BnyV3xGROAld15AkGU0tRvEuQVOLsTiGUZl/PZEpX1ORO5bI5GlI5iTUQCaSZEY2p1KefS4AkhSJydoNe8ydqP7NWBynoSmZ6LqXqoLbcSW+RcA9B03Nx+o8C3fxU9gir8bqPJfKvGuITJ5GRd4VOOMeR/Gtx+o8E8W3AtncBtncGnR3tUjU9p6fGq9DU5UcZFMy6AGQQs2QRrBJCgASuu5HkqyAjKbmoWsVmCwd0dVikMzoWhmyuQ26VoqulqCjAmCytEcN7EKSzEhyTFBMS2Z03YumBD9iMXPWIlYsX8cLzz9GUnJ7du3ahcVy7L+WLTY2lnHjxrFgwYI6ty9ZsgSPx1N7g66BrhjXVNc96LqXgHcx/ipXsLOj+9D1AKAjm9MxWbpisrSuXZdAcIDouo6m7MBk6Vh/Hs0Lkq1Wh1tTS5BNsQe2P82DJO/9jLOmFuItn4Iz9v7gb1wPgGQDdCQpXEYFbd2Nt/z/cMY9hiTZjDoU7wosjmG4S17AGfsQSFZAM+rwln2ILeraGpX5QLKha2UEPL/jLXsTi+tSJP1UFH8+3pJn0bUiHDH3o/j+JuCZi8V5FiZze6wRF4Gu4av4FFvkOKqKHiIi8Y0DOg/NgRClRwCTtSeqf22tdF/l19giLmkGi5qP4ENNMhoOXatE00qRTUGvcVCEgLfsDWS9D7rmBDkKTS3C7/4Fe+QVaGoBIOGr/Aqz7UTMtn4EPPPQdR9WxxkgOVC8f+IpnYTZ1hdd82CydkbxrcMRcye6Voav8ksU72JsUdcjm5LxlDyPydoH1b8as30Yivd3XPEvIZmTkU0plOdcBLrbOA6L40w0JRNV2YXZ2gfZ3Ap/1Vf1HndD2wA8Jc/tk/+bevOWZY2oM91X8TEAFbnB31RF3hX11qHrFSi+ZVTmXV3n/qsK7jCWFe+S6vo/wVfxSVjd7uKnq+0NP77KAgldTUIN7EbSvaiB2jaoWsNfN6sLNbB9n/VtdeTZCYCuB9DVgtrbVY0vv5rHc8/eRnxcAE3ZTfv27Q/YlqORk08+GdBYvuwvdD2AplWArgIKfp+Xzz79mFtuGcejj21GrRbw6CABmmZGQgZJDj6UdR1JsiDJcciSrfrhWYK7+Mla+7VFXo0t4pJ6O3yCI4uu+0BXkGTXfvNqahmyKRpdDyBJdXfcdF1D1yqqO4rp+Cu/xuoaZQg7XVfxlPwHZ9zD6FolkhxhdKyD+ygGZCQ5ApBQfCuw2E+iquhJzLa+mCwdqSq8GwCzbQBW1+jqzq8fT9m7qP7VADjjn8NsO5GqgruwR99KwLsEf+UMLI7hmCyd0dRcNCUbxbecmo4ii+MMJDkaf9W3tY5NkhPQtUIAytw/19putg1C8f1V53kpzz6vVpqn9MXgtpzzjTTZ3AZN2Q2Ar3JqnXWFCFR9iZ0v8RTVrPOlvdvdvxAAvOXv703z/IGmZlGWNQLJlEpUymcN7qM5EaL0CGC29ccefStma1cg2KMrzzkfb9nbWF0XNvnQbUPouoYW2IHJ2mlvmuZFku3hedRsTOZW6JoHb/kH2CKvIuCei7f8fSQpEl2vCGaWXMimRCIS30JVdgM6mpKN2dqNquKnsLpG4a+YgaZmAyao9mY1Bhs/UJkXnuYrnxy27q+cEbbuLX0tbD3g+Q3AaDQq82/Yp76PjOVQw6Z4fwegquj+em0LeOYYy4rvL6inURK0PIqKSvH5/CxZsobvf/gdl8vBJZeM5eyzz25u0w4YXdcJPlwV0NWgt1hXAHXveg1vdRC52tujBD2gkhkJO199M4tTh51BUnI3AEzmdGMfiqJgMpmNzqQku5FkB2bbQCyOvd5QrZ6ORs0OjWxui2xugzP2vmoR0jT4Kr/EYj8V2Zx8SPXougq6P8xTFkJT8pBMcYY403U/mpKHbG61t6Ot6/irvsPiOK1a9JkBHV0tRFV2YbEPCt+f5gbJAbq31j4V/ya8Ze8ckIdL131Ikq3ae60BOt7yD7FF3VZ9DHvwVXyErpXjjHsE2RSPrqvoWjF+9y9Gm2hxnI6uqyjeP4hI/i+VedcA4Ix7Cm/Fp2h1dAhr4i1/r1ZaWda+o4NWwN9gPYr3j/B13/JqUVkbd9EjxnJV4T3GcsAzn4Bnfr37CHh+rXdbSJDWa18TtP0hQXq40NQsY1lXc+r9iE9LoEWL0oDnD/zu2QdURtd1rKofd7G15Zx0XcHq2tsrkmQ7JmsPVP96yrPPwpX4NmZrt8ZVVT2kqesaulqIZEoEFEAHXUPxr8NT+go210WYbf3RUVB8K42GRrZ0RKv2Npls/VF9Kxu1X3/V9zVsqKhhUBWaUkV5zqg6y4WLxMYLUsHhwxn3DKqyEy2wC133ovrXodcjJhwx9+ApfRVJjkfXinAlvBEWG33kkKktrupDAvSgN2afSYWSHEFx8S7cbi85ORW8/+4j5OZV8fgTb9CqVSt69OjR1IYfMMHRBLUBkbnvfSQjYQLJBJiqRaYNZFMwHblWWyjJdiTZgWyKBiAvL49Fi5bw2muvUVJSctC2y3IUUWmzQffi9/xaq5MIoCkZaEoG5TnVYkNyYnEMI1DDC2W2n4rJ3Aqr63zUwDZ8lV+g+tfjiPknFscwNDWPyvwJyKZUNDXHKOct/xjZ3AqLfRCyuS2SHIW76CFAw2TrhxbIxBX/LJUFt2CLuAJkJ77yD8POZeh3ZrKcgCQ70HUN1b8Kq+sf+x3xqIm37M0DOHMNU5Z1NqBgdV1gtMXBEJ951Z5rFTWwucE6/FXf4ACqagweVOSOqTd/qEMPGIIUwF38xMEcQn1WNWFdgsajAC0zXKlFi1KL41QsjlMPqIymaVTm5xMTl4Qst9yXC7gSXqc8eyQAVQW319oelTbLGMrWdRVf5QwkyYm3rHE9Zm/5B3WmazWGPxsrSFs6FudZWByn46uYjsU+uFYP3eI8h4DnNyKT/g/FtwZP6UuYLN1QA5sAkOQYXPHPoyoZgAVv2VvoWgm2qBvxlU/GbD8VXa/CYhuIt/x97NG34i17F8mUWOfwcFTaLyjev5AtbZDlBMMLrfhWga5htvcHCBu+0jVP9d9ydL0SxbcOW8SFQHD4TJKslOeMwh5zNyZzWzQlG0/pS0SmfIGue/FXfouuVyCbUtC1MvxVeyep2KNvQ/GtMzzAFscQLAxp9Pm1usI7HFGp3+IpfRtH7L3oagkVeeOISPoAXSunqvA+bFE3YCnbA6wKK2cyt0NVdhnrvopp6FoptT3oIe9J6K85eA51nZDgBJDkKKwRF2GydKn2bNkaPI7QhDKHywvAFeOuxemKpH17N8OGDWP58uVNIkr3ei9r/KsWmjoa6CofrvdT7K2ZZ1+kvf8kCQmJ4HkyAxJxdhM39mq6ofDJkydz5ZVXYrVaD7kuSTKDFIHNNRpv6WtEpf6IrvuNsJJa6O4wQQpB75gC+CqnhaV7Sl/BU/qKsV5TkAbr8qAFtuILbK21G9X3NwCVBbdAHXVX17g3f2Bj2JYDEaRNj1Jtw17ngLfs7eYy5pjCbD8FxbsIAJOlK6qyB0fMXXhKnsfqugSTpX11vHwrFO8y1MBm45ngiJloCPbIlOl4Sl5C1/04Yu9H18qRTfFIUgR+9yy8ZW8TlTYXdC+amlc9cVTFXfISVufZSKZ4zLbegIWKnItBMlW3j3uRTek4Ez+iJG8KTkcVkmTC4jiDgOdXrM6zg+2rHkDTSsNGAcOQ7NTd5rQMWrQoPZaRJAln3ONGTN6+lGefi8UxAk3NQ/WvO8LWNUxkyv+oyL0cCAo+SbIR8C7BEXM37qKHMdn6YY+6AdmUSkXuP3Alvo2v4n84457AXfw0zrhHqoeqIlB86zBZOxvB4CE0TSM/P5+kpNqdi+AMconQxJhQ+ENoWMxsHwS6H9nSfm9oRGxwKN5iSkXX/Vhd56CrpWFvQDBZuwTzOIYAJiTJhD0yPC7TGvEPJMmELeLSaluCQkfXKoJDe3JkdUNRW/SZbX3D1mvOEg8N2wX/JocF9Ie8WWET42y9sbrOMVYdMbeF1W2PvgNdzUc2pwSzR/wjOCzJoY8eSHIkzrgHg8vm5DC7olJ/QJKd2CIzkc3ZyHIMOu0xW4Iz/o0QFl3DFf/kIduy16iGBSlg/BbS09Mxm0NNnxR8TZvuDYatKAUEG2wNnYZE4373ZvwLCsrqWMxqj+aNPV3GcvC31ryjOmvXrmXTpk289957aJqGpmlcddVVPPvss7Rp0+ag6w39NiQcRKfPQ9d1yrNHVnveywlOahO0JKwRYwhUzULXK7BFXomv4vNaeUIjJ/bou9CUXcimJGyRV1THlpZQkTsGR8y/8JS+SETi+5isnajIuxHZ3Lra0RT03rtLXiEqZSrlORfiiLmvegKdB9mUhKfkBayui7G6zsVX8T/MthNRA5uxR9+Gp/RN7JFXIJmSwsLfVCUbb9n7KN6FRpoz/j+Ybb1QvCurQxX2TqAMD2OrCouxtTpH7l12VceGRl1X61zUbP9cCf+psSXNWLJFXLJ3/ojkwCS3M7ZFJE6qVacj7hFM5lbI5lQgOMlWkhOQJAlN01Dk0dij9z4bzbaewN5nGIA98kq8ZZOxR9+It/wTrM4zjfpaMkKUNiMWx2lEJk9FRzEmnNQk4Dn42fnWiDFh8ZYWx+k4Yh+uFlAekEyo/g2o/s1YXaORZGeNmxUCnoVYHEPrrFvXdSKTpwIqsjl44zkIDunu+0YBV8IrmK3dMFcLEFd89dCPFIwlC91MB0LIxvowWep/iEqShC3iguByPa/kaqj+feN/Q+uSHInJ2nImcUiSjFQtSPemHf7YZUmu+UolCckUA3UIrsa+tulwYLfbOeWUU5gxYwYTJkwgN8fDwoXLePCBfyLJ0dX2hsRk9XJLCQU6RFRVRVVVNE1DVVX8fj8mk4n33nuv2sMLBQUF3H///bzxxhtERTXt62UkSQprIyryb0Gr9mqGxao3AtmUVh2rDvaoW+qMYQxhdf0DXSsm4FkIBPZOarQNQPH9jdV5DqqyB9mcTMD9S511WJznEXDPxOq6MGwkotYxmoKdSsW7uHrfwSF3V+LbaIHteEqDIsSVMImqwntrl5fjq71zy3ElvoUW2G1MjkGKAL3+91xbXRdjcQwLi6eMTp+Hphag+DMoKTURE1WMzTUCXfOiKhnIpkQkOQJJsqIpOcjmVBzRNxvl7VHXG8sB75+4i58jKjU8lt+wXZKRTPFEpf2CJJmwuvbGabsSXq5+I8beOOJox2kARKbMQDbFhx+L80xj2RkXihUNdsSdsbXPG4DJnIYz7nHUwNZaIXEWx+A6yxi2N2LS15HCYh8Ytt7QW18awh59Y/BvVG190VIRorSZCQXlh7wIulZU/QqpuuPnolK/JTicuTcgXtf9+N2zUf1bw27Wmg1LTUJlzba+Yd67mmKsPkEazCchNXIywb7eQYGgJXDLLbfwxhtvMH78eCIjIxk/fjw9e/VrbrMOO9OnT2fatL3D1jNmzGDixImMHLnXK+T3B+P8YmODr9EJidXDQWRSUEiqge3IpmQCnj8wWbugBnYEPWURl2GPuhlf+WT8ngVEpXyGrnnQ1GxkcweCoRHeau/8ZWhqCZIcg66VV4d8+EFyGh0yXVfQtUpkUwy6rtR6nU/wwB8wFhXvCky2vuhacVAYxP6zeouEPfoOyrNHVrfd4W8V2RdHzMTggrUbFscIkCxIkomo1B+RZAcB71+YLN3Q1CzM1u7VtlbPdreegMV5VnCvNd5aUp5zIWDClfAiJmsvasYOO+NfwGzra0zIkk2JmG3x6FI+FkefYF2y3Ri5CLE/T5rFfjLRaT80mCdoZ+0OsGyKqTf/voL0UJAkU6PnaAhaHkKUtiAkSUIyJRCV9gtaYCuqkoHqW4s14hJMlnYNlLNic42GltPREwhaNBERETz88MPNbcYRZ9y4cXW+ML8mycnJR/zF+aFwFavrXGO9pqfMbD/ZGJWRZAcmORTeIoG01zsfeh+lVB3ysm9YhySZg9776uX9YbafWF1fuKfKEXMnUCM84QA8/zXfdBJyEIRCj0KhOsE6LTWW952oFvQ2RqXNoq6JbPt62gSCowUhSlsgkiRhsnYJxofUaJgFAoHgeMRs6wW2Xs1tRotCfHxFcCzScqenCwQCgUAgEAiOG4QoFQgEAoFAIBA0O0KUCgQCgUAgEAiaHSFKjxAL9nhYlus7rDNZBQKBQCAQCI5WxESnI8B/lpWyMj/4mpVnhsTSJbZlft5LIBAIBAKBoLkQntIjQEiQAjy2uIQnlhz8t6UFAoFAIBAIjkWEKD0CjOkS/gLRTcUBft7lbiZrBAKBQCAQCFoeQpQeATQdnh8ay/+NTDDSpqyv/1NxAkFj8Kt745M1XUdrIfHKIm76+GbsT/m8s7qcPLdKuT/8y3RlPo3dFUqj6vkt08OaAn+tdEUL/r7cgbq/etcYynwHX7alElCb7r5TNXEPC5oHEVN6BNB1HQmIsslMOSuB634p3Jt+jHxT+3inyKMS7zi0b8uX+jSirBJyjd/EukI/GmCS4IQ4C6sL/LywrCys3NA0GwuzfcTbZYq8Gk8NjuGJJaVc3tXFiNYOIqwSfhWKvSppEYf/lr98ZgGj470Ue1WyqxRiVcj2KMTbTfhUHatJItIq41d1ApqOX9WJscl4VZ3sShWAJKcJv6qj6DqVfp2OMQcehx26v4KCHcxy8Lxeetll1LzrfD4f1113HRdffHFTHH6LZubMmcyePZuMjAzGjBljfN1p7dq1PPLII9hse7+A9OSTT9KjR48Dqj+7Mig4F+zxsmCPF4AO0WZ0HbrGWfh5lwdZCnbU/zUgmheXl5HokDmjtYPBqTZSXSajTfxpp5vdFSrTzkskt2rvb/fKWQVG2VdPi6PMr9E6wsz6Ij/JThMRVpl4u0xAg00lAZ77q5Sbe0UyrJWdK2cVGLbePyCavolWzLJEvltF0XSKvRo9E6xous6WkgCv/13Otd0jUDQ4Jd3OW6vKuKNvNB5Fw2GWqfBrRFqDvp08t8qH6yp4cGB02D1cF6U+Da+ikeIyk12p/H97dx7eRLm2AfzO0izdd8omlEXKKtCySDkoKCLIdlBWOQoogrIK6gEBleNRQVRAkQMIgiAI+rkhoCiUfd8UWva1WArdlzRtmmTm+yNkIDQtLaRMEu7fdXlJJ7O8mSczuWfmnQmuFFiRbrQit1jAgAb+EETbtnAsoxiBWiXiq+lgEUQoFYBSoUC60YrlJwxIN1rxSvNALD+ej2MZZnzxeDi2pxShdqAaaiVQYBax4YIRU9oE42SWGSGibT/wZ4YJnWrqIIjAwqP56FhTh0ZhN35metCv6ZjSJhh1AtVIL7TC10eJKr4V378Jooi8Ytv2TVQeCrEcpzWMRiP69++PNWvWwNfX93ajy0oQBKSlpSEyMhJKpXtsCKtPGfBwVR1qBdp2qouO5mHzZdsOe81TkXI2zW3dyzqWdnDw7SkD+jXwR2JGMX67WIjX4oKQXWRFisGKd/flYMFjYdCrFbhSYMXkndnoVdcX7apqUStQjdxi22YVrFVCEEX8eNaINlFahOqU+P5MAa4arbAIwOTWwQCALZcLseBoPrrW1uPXi4XwUytQYLn7sxVBWiWq6JU4nWMLCwEaBeZ3CkeRVYShWIBZAML0Shy6ZsIjNfQl1otFAHxUtnUzYH0aPn4kFGE6JQ6nFaNZuAb5ZgF+PkpcLbBi2u4bfaWF7FRodyzBK+NeRZVqNe76fZSmqp8KxYKIzEIBWpVCCr2iCJgFET5K4OYTavZAZJebnYXp41/CwoULERUVVWntdBd79uyBUqnEtm3bUKNGDYdQ+tlnn2HRokUlphFFERaLBWq1WtpOLl++jE8++QQTJkxAzZo1pXGvGCx4dVvWXbezdZQW+6+aSgxvGu6DYxnmu57/3Xi7bTCm782R/m4S5oPEzBttCtIq0Sxcgx0ptn38AwEqJOdbMesfofj5XAFqB/rg65O2K2XDGvvjj+RCXM63OixjcqsgfHDLwafdgyE+OJ3tmnXgq1bAeH0/o1ECxWWcQP76yQi8vScbj9TQoW6wD0J1SujVCmiUClw2WBCsUUKtVMBHpcC3pwzYdcWE9tV1+OW8EV90DodGCejUSuQVCzhw1YTHHtDj73wLqvmr8NmfeXg2xh++agV8lAqIsF0JuphnwdKkfLwXHwKduuT3QJ5JQIBGgfRCAUaLgAcC1A4HBBZBhEoBFFpEXCmwolaAGhlFVlT1u3FwXmwVUWAW4O+jhI9K4fYni9wx45RHeXOkW58pFUURFf1atp8VEUQRcJPLiIII3PwZf6lZoBRK+69Pw7Q2wYgJ9ZG+SH19Kv5Bc4cNySLYznxpVQoYzCICNbZA5uyswfHMYocj85vlFQtQigJyiwF1kRU+KtuOfveVIjSL0MD/+voxmAWoFQrkFwuIuH4UX2wVseKEAc838odaqYBFEJGUWYyHIrQ4nW1GisGCFccNGNEsAJ8czsPKrhF49td0LOsSjrd2ZyO1wIoHQ3yQaxLwt8GK78/e6Pvbf32aQztHbs50+Pvnc0b8fK70vsLfni4oMezWef56sRAAXBJIAdtlypsvVeYXi/jXb+lOx53/V77T4RF6JdILbfOY4ILA4UqpBTe+zE3XL1/e3K3h1iu8t16VPLR7O2rXawC/0EgYigUoFLY+TQqF7ed+pX/DFmiBkr9D7kkefvhhAMDBgwcrZf6u2uM6C6QAZA+kABwCKQCHQArYtjl7IAWA5OuB8/Udtm1n55Ub7+3LUrpxlRZIAbgskAKQAilQdiAFgMHX9xvncsvf9eyX87b94fA/Mkq8tuiY4/5m9xXnNQeA5zdm4B/VdQ7rtTSP1NBBFIEm4T6Y/1c+ukfrse5CodNxB8f4SwcIN/u4QygmbrfVK0KvRKeaeujUCvyVXowIvQptq2qx8VIh2kZpEV9dh+wiKxQKBXJNAmoGqFBkcf49braKsGdrhUKBYquIvalFqBGgRp0gHxzNKEaYTolIvUo6GQDYDvbC9Cpor4fmm+1KKULTCA0CNSWXJ1w/saBRec4+y61D6YFrxdj2t/MPU6lEwGRSQPt3PuAmdVAACLrlA9MsXIOjGdcfE7Uvx+E1jRJY0dX5GVT7JRwAyC4SEKZXQRBFDNyQjjVPRWL1KQPWnjNiVbdInMk2Q6UEogPVEAFkFAqI0CuRV2zbMP5KL8aFXAv+Wc9XurRpEUQsTTLg0DUT5jwahhkHcnAu13aWrVGoD45n2XaI9kvGZZnRPgSTdmYjQq+EAkDa9WBjP7MwoIEfVp8qQLNwHzxd3w9v78m5ZQ5KALazb80jNPjTSf+y0vx+qezPzSeH8wBAupw3ZOONnWZSpvxffO7EHki90cFd2/CPzt2ggC1QCQJghghRBETYDnBt/8Zt++yKsIdXBRSwhVnl9UBrC7sKh4Br/7cCthHsuytpt6VAyWG4cYArveaikJyeno7BgwfD19cXHTt2RL9+/aBSlX7JVhRtZ6ByiqwwC7YD6r/zLYirosGJLDPebB2MUJ0SL99y8EZ0J8oTSAFI3Ua2Xx+/tEAKwGkgBSAFUsC2/1tzywmFP5Jt89x/1YRP/8wrdf4PV9ViT6rte/LFJgFYnOj8wP9WN5/Btqvqp5IOwhsGKnDiUAaWdA7Hp3/m4YlaenSL1mPr5SJYReDJ2noYigWcyjbjl/NGzOsU7mwxbomX72V04KoJHx1yfkT8ykMBaBGpRZrRiv/uy0GPOr5Oz7SR+xjW2L/UMx+VpXddX5zPteDhqlosTsxHi0gNWkbaLn3q1QrsSbX9v9BFZ17Ly9nle51agaJb2lHFV4VrRtuONkSnRHZR5QRgjUrhcAYVAK4kX8TsdyZh+mdL0LR68F0vw74rFa6HWFEEhOvDb/y75DDA8aKOePP/xVuHiSXGvxMrFv8PQUEh6Nl3AAAgNycbxoICVKlaDVevpGDhnFmIf/QxPNG9l+29CCIUSoUU3q+m/I15c2aj/j9fRkiVavBRKuCjBHyUCrSsoi3xLGaDWUBOkQC1Ehi3Vb6z7QsfC8MIhuRyKesMI3mWIK0SdYPU6BbtC41SgQah9/5Z6V5x+d7btYrSlvqa7VLqjaMqBlL3NffRUERd76PUpbZtYxuwPg1zO4Zh7JZMTGsTjLrBahjNIv69Iwv5ZhGzHwmFwSzi4DUTiq0i/lnPT7oZzs5+yea539KxqmsEVMobZ8RmHshB5wf0aFnlxmfo0Zo6qavEYw/Y+ocW7c/B+JaBeH5jyctnpelV19dpN4T6wWqcySn7zmkfJfDxI2FIPGfA0h22I/4oPZBTrEBVPxWUCgXO5ZgRolPihzNGZBXZAqn9LKNw/f9Gsy11aVUK6YYNlRKwXk9xvj4KhOlUeLFpAITrNxIqFLZLWyarCB+lwmF9AbYgWGwVoVUrIYoithzYgVatWrkkkAI3zliWvFJW8kzmBx98gLS0tBLDyysyMhKTJ0++4+l91Ur4a5TS5zbKLwJABACgWoNoPDtwANatW4fn+j99U5/SGzchWfQqBGqVeDbGDzVrBt52ef4+Sqnbjb0f/a0H5c4+d7UD1cguskp9tAHgoQgNxjYPxMubMxCuV6FOkBo7r5jQ+QE9YqtosDgxH7GRWvSq64sCi4gHAhy/5uzLH52QgfRCAe/Fh6BOkFo6YLJ3+/H1UeKPS4VIuFyI5xv5Y0liPpLzrQjQKJBf7HhUMK9jGEZvcQy7jcN8Sr3qMqtDKNafN6JzLT1UCttZvR51fLHpUiEeranHrIM50mX/+Z3C8MnhXJzNsaB2oBoX88r39IKbfflEOD4+lIukTDOWPB4KH7UKmYVWTNyehTfigqBUAO/vz8XyJyOgvekDnG8W0be+HyJ8VfjhTAH2XjXhUp4FtQLVuHSbdox+KBDz/rKdRZzSOhjv7c+pcLvJNXJNAg6nFUvPTF/WJRx6J3103QFDqcwWdw6HwSxAAUhfELf2M6wM9ju171YNfxVCdEqHvl5BWiVq+KvwVLQvfFQKvHdL94S6QWq81TYEP5wtgEalQE1/FT45nIeXmgbg0Ro6iAAOpBYh8WoefH19EV9dB38fJTZeLETfB/3wx6VCRPmpMOtgLgbG+MFkFdGuqg4Tt2fh4w6hqBGgxtgtmdIZuLJMahWEArOIz/7Mw3vxIZiyKxuto7RQKYA9qSYs7hwOjUqBzEIrlArbmb2sItvNPZrrVzad9Zldff2Lb2XXCKlrhF4NfNE5HKk33Ulc1q972cdZ/mREiYD171bBJcZ31o5J12+kWvNUJATRdnfxzAM56F3PD8FaJQotItKNVjQO12BzciF+u1iIQTH++Gc9W0D48Xqf2kCNAk/X98OMA7moHajGC00CYLKK+O9Nte1SS48BDfzg66NEcbAPInxVCNOroFEC1fxvBJpwvRJ6tRIvNg0o9b3fqvh6X6zL+RZE6FUO/bVuft8KhQI6tfPL2QqFAtrrr4miiB3bt2HUqFHlboMr3U2gvBfsAb8ytYrSYs1TkRiwPg3zOoUhXK/CoBh/HLpmwqJj+ehRxxdPReshwnZQ8uIfGWgU6oM3r3+mF3UOh1phu7FmTIsb832mvoDI65+9sDKW3766DjtTilDv+pMdfH1ufG7sffA619Kjcy3bAd6sDjfmllNkxa4rJjxVxxdpRisifFX4plsEUgusMFtFCADqBNnmu+qkAUFaJXJNAs7lmDG8aQCi/NR4+aEbYT76+rj9GvhLy/r0SC52XTEhTK/Ce/Gh2H/VhEahPvj4UK7UjWr2I6F4dVsWnqnvi84P6OHro0RWkRWbk4uwI6UI/40PQZhOCYVCgclxgUhLT4evjxJKpQLV/NX4ptuNbmJPResdAikAvHJTG/vU90PPur5QXb9ZcMOFQiTnWxCkVUp9R99ua+uyEaxVQadWoF6IGt+eLkCzCA2GNPLHrxeNyCsWpSs3q7tFYMUJA7b/XYQQnapcjwvrWEOHvg/6IeFyER6uqnW41F4R9ptKe9X1xYAGfsgrFrHqpAGvPBSI87lmTN6ZjWVdwvH1CQM2JRdhaptg/Jlmks4eP9fIH8uPl+/K2ISWgbhqtGLVyQK3uFnv1jq7E16+d0Nnc8yYsuvGncxT2wQjQq9EsQDpqD+v2Ha38XM33bTyUIQGz8b44+dzBdh1vb9npF6JtEJBCqGrukbg90uF6BrtC6NZgK+PEv3Xp+Hz618KRRYBOrUSBrMAH6Xtkqf9kSf2x6DkmQQEapX4/ZIRT9S68XlYfjwf/R70LxEK7Hdxq5QotQN4jklAgM+Ns1uuqKP9Jqv+69NQM0CFCS2D8L+/8vBSs0AE+CgQqFWWCHIZhVYEaZQOnczvF9lFVixNMmBCbJA07EKuGVX91FJNvziWh/O5FnzQPhSA7fEyL2/OxOjmgWhfTev07uyoqCiHO7fdwZEjR/Dxxx/jq6++KrPfpLexWq2wWq1YuHAhgoOD0b9/f6hUKhw/fhxRUVGIiIjAlStX8MEHH6BDhw7o27dvhe6+d5UCs3D9LPmNz8wflwqRXyygT32/MqYsP4NZQIFZvKNHHd0rZ3PMUmi+2dgtGWhbVYdBMf5Op1t10oCD10z45JEbQbqyvxu/OWlA12jfUh//JIgizIItEFmuPwrO/l0gXr9B+XSOGQ1DNVIotJ/VLraKsAgidqea8PgDeqfzd3Yyx96Ps1NNHRIuFyFYq8TDVbVoHKZBqygtDl8zOVxtsrMKIgb9artPwyrYOs3Yz6DvSTWhfTUtrKLtngT7kxLGtQjE3CN5eD0uCE3CNNCogDPZFgRpFQ4nnNY8FYmsIit0KtuTBvZdNeGzP/PQJMwHVf3UaBWlQcNQDT7YnyMdfAC2AD9gw43ve/tTHUpT1n0Ycjz1p7w5kqHUjdmfPlDWM++GbkyH0SI6fMj+zrdg4vYsfP1kBNRKoMgqlnmq/qvj+Xi+UfnPWt0LrqyjOzyZwFscuGpCVpFV6qYAAEfSTGgQ4uNwsOHuofSTTz6Bn58fRowYIXdT7qlVq1bhm2++cRg2btw45Ofn46effkJBQQGCgoLQsWNHDBw4ECqVSpZQSmVbe/3RUs0inD/BpNAiwGSFQ0D0pu9GZ3ZdKcKnR/IwslkAFhy1dX27+Xsx3WhFqE5Z4qpTaXJNAoIq+HzVC7lm1A6s+L5uWVI+2lTVomGoYz3zigX4qm8cnBVZRFisVmw/l4EnGkTg2d8y8a+G/sgvFjCwlAOUW8N6WTdSVyb2KfUCipvuyC3N0i4R0qUTuxoBasx5NBRqpW0e+lIuZ9q5WyB1NXcKQ57OWT/oFpGl9412VxMmTJC7CbIYNGiQ9GzSW90PPx7gLXrWLftssV6thP4++3aPr6bDp0fy0LGmHq2jtEjMcDxLGFHBM+IVDaTAjW4YFTWksfPv4Fsf86RTKyAolWgZauupPrCBH7rXKftEYe+6vvjpel/tcL0Sn7v5nfj32cfWO/Vw8qG8+eHARERE3q57tO3Svp+PEm2q6mRuTeVSKBToXe/2XVnaVNXip3NGfNMtwl0e3V4mJhciIiLyeP/y8qt+d6JOkM+NbgwecNHQ+zqWEBEREZHHYSglIiIiItkxlBIRERGR7BhKiYiIiDyEyWSSuwmVhqGUiIiIyANkZGTg6aefLjE8NzcXFy9eBGB79FtxcTEuXLiAuLg4mEwm5Ofnl5jGHTGUEt2BdevWybbsq1evwmq1Yvz48bh8+XK5pvnvf/+LZcuWVW7DrktMTERGRsY9WRbR/cxqvf1PKXuCHTt2oLi4GBaLBTk5ObK2JS0tDRbL7X/u1Jlnn30WAJCeno60NMeH1l+5ckX6944dO5y+z+TkZPTq1Qtr167Ftm3bpPmNGjUKhw4dwv/93//BbLb9ylNSUpLDtIsXL8YzzzyDsWPH4vTp01izZg369u0LAIiPj8djjz2GlJQUJCcn39F7u1cYSt2YxWJxu53OtWvXnP4mdkpKCgCgsLAQgiCUOv3UqVORmZkp/W0wGLB9+/ZyLz83N7fM+ZfH2bNnpX/v3bu3xMZ9O3l5eXjnnXeQkZGBPn36ALAFxYKCgrtq163i4+Md/v7jjz8AAN27d8fu3buxc+dOLFq0CIBtR3rs2DGH8Y8ePQoAeOONN/DTTz9h3rx59+Syz7///W9s3rwZgG1HLAgCcnJyYDTaHuDsrH72z5QoiiU+X+62DVD5FBUVoaioyOlrxcXFKC52/hOIni49PR3//e9/bzteZmamtA4sFovT/eqtcnJyHLaHNm3awGq1lghAN++LioqKYLFYpO2vPNatWyfV7v333y/3dHfq1VdfRXJyMtatW4c+ffqgqKgIaWlp+Omnn2AwGGAwGPDrr7/ixIkTDtOdPHmy3PuHy5cvl2vcbt264fXXX0d6ejoEQYDZbMbGjRsBAEOGDEFubi5+/PFHWCwWmEwmPPLII5g5cyaOHz+OU6dOIScnB127dsX06dOlecbFxaFnz57YuXMnBg8ejFdffRUvvfQScnJysHXrVsTFxSEuLg59+vRBSkoKZs+ejYkTJ+LUqVOIi4vDvn37MGLECMyYMQO9evUCAAwdOlSaLi4uTvqltt27dwMA5s6d6/C+BEFAr169pO8sd8WfGXVDY8aMQWpqqnQq/mbr1q1Deno6/P39YTabUa9ePZw9exbR0dHw8XH8NQmDwYAjR47gH//4h8PwzZs349y5c3jppZdgsVhw9uxZhIaGwmKxICQkBIcPH0ZKSgo6deqEJ598EgDQrl07NG/eHPPnzwcAzJ8/H6+88goaNmwIg8GAy5cv46233sKMGTPQs2dPmEwmHDx4EKmpqZg2bRpatmwp/WKMn58fCgoK0KNHD/zyyy8AgA8//BBvvPEGZs2ahddffx3t2rXD+PHjkZaWhqNHj6JLly7SJYvRo0dj3rx5AICXX34ZL7zwAsxmM6xWK7KystCzZ09MmDABcXFxqFGjBnx9fVFYWIikpCQcO3YMn3/+ObZs2YJdu3Zh6tSpAICDBw8CALZu3YrXXnsNgwYNwqpVq7Bv3z4cPHgQsbGxaNu2LYYOHYp9+/bh+PHj0vqcOHEiPv74YwCAr6+vQ8jev38/NBoNEhMTMXjwYAC2L2SDwYDQ0FBpPPuXtL+/P44fPw4fHx8MHDgQAPDEE0/g999/L/XzEhAQ4HBppl69etBoNDh+/DjeeOMNfPjhhw7jL1u2DE2aNHEYJgiC0+0lJycHwcHBTpdrsVjwwgsvwMfHByNHjkR2djZatWqFrl27QhAENGzYEEOGDEFsbKx0Rjc6OhoXLlxAw4YNpS9ie1tv1qhRI1y8eFE6e9KoUSOHn7sUBEH6zfri4mJoNBoYjUZotVppuLP3ZDKZoFKpkJycjPnz5yM5ORmBgYHo27cvunTpUuo69hYbNmzAxo0bcenSJfTr18/h151SUlKwYMECnDx5EjqdDv3790f37t3v+GdG33rrLWzYsAFbt26FTqdDt27dsHbtWuh0OgwfPhwhISE4deoUqlevjuPHj6N58+Z477334OfnJ+37goODodPpoNM5fxB6UlISGjduDFEUcfXqVRgMBtSvX99hnKtXryIgIAAXLlyQxr1y5Qpq1Kghve+oqCjpcwPY9p0nTpxAeHg4oqOjYbFYUFRUBD8/P4iiiIkTJ2L27NkOy9m8eTNmzJiBoqIiFBYWYt++fVCpVNiyZQtWrFiBDh06ICcnByNHjsSMGTOwbt06dO3aFdeuXcPFixelg7ft27dLn1u1Wg21Wo01a9agf//+aN26NWbOnImmTZtK++YxY8bgs88+k/ZhgC0Evfzyy8jLy8PKlSsxbtw4zJ07F0899RQSEhKwY8cOadxff/0V06ZNw86dO1FUVITAwEC0bt0aa9asQd26dREXF4eDBw8iIyMD4eG2XwMSRRHJycmoVasWAOCVV17Bq6++ioCAAFSpUgUKhQJffPEFhg8fjpSUFERGRkrfT3v37kVRURHq1q2L7OxsfPTRR073Vf/5z3/w1ltvSX8HBARgwIABWLt2LdavXw+r1Yo2bdrA19cXf/zxB86ePYtPP/0UXbp0wfvvv4+DBw/CbDbjww8/RIMGDTBjxgwsWbIEe/bsQdu2baXPX+3atWE0GnHixAlpP36zxx9/HJs2bcKiRYvw0ksvOf0cluadd97BO++8U6Fp7oWdO3eWuk1VlvLmSIZSN7Jy5UpotVrMmDHjjucxZ84cjB8/XtqJ2S1fvhyrV6/Ghg0bEBoaiqysLFSpUgXXrl1zRdPpFs7CoDNDhw7F0qVL70GLHG3duhXff/898vPzUbNmTbz77rsAbJeAoqKiYDKZcOzYMbzzzjvYvXs3tmzZgp9++gkzZ87EyZMn8corrzidr6+vL4xGI1QqFaKjozFy5EhERt797yw3atSoRHCtX78+UlJSYDQaERISguzsbABAUFAQCgsLUVxcjJiYGJw8eVI6EFIqlRAEAfPmzUNMTAyGDBmCtLQ0vPnmm/j444+9/nfc9+zZA6VSiW3btqFGjRpSKC0uLsaoUaPw7LPPIj4+HmazGZmZmahZs2aFQqnFYsH27dvxxhtvuLTdzz33HBITE9G5c2fMnDkTnTp1QkJCAgDgySefRFBQENasWQPAtg+cM2cORo8ejUceeQStWrUqdb579+5F27ZtERMTg9q1ayM5ORkfffQRunXrVu62/fzzzzh+/DgmT57sMLxRo0bw9/dHdnY2zpw5U+H3XKdOHZw/fx46nc7hjPMHH3zgsKzq1asjJSUFu3btQnx8PFq3bo39+/ffdv63blNff/21dOAMAI0bN8bVq1cdrmzdasGCBXj33XelK2UVERsbi0OHDlV4OgCoWbNmubsuUUkHDhy45z+/zVDqQaHUYrFgzJgxOHDggNxNIXIJV4dSV5s2bRrGjRuH8PBwREZG4v3330ffvn3x8MMPy900B/buDPb/BEEo979Le10URaxcuRKBgYHo3r07AGD79u04f/48hgwZUmobbt6fXr16FfPmzYO/vz80Go00/OTJk0hPT6/09UJEd+7ms+r3SnlzpFv/zGhCQoJ0ebciTCYTtFptJbTozmRnZ+ODDz5A1apVpWHLli1Dv3798Prrr2Pfvn0yto7o/hMfH4/Dhw/jsccew6FDh5Ceno4HH3yw1AB4N/+VNp/yUigU0n9KpbLMf6tUKumMZmnjKhQK+Pn5ITAwUDrLmZ6ejipVqmDu3Lm4evUqYmJi8PLLLyMsLMzpmVIACA0NLXGmdNOmTZg0aZLrCkVEdyQ4OBgjRozAzJkz5W5Khbh1KO3UqRM6depUoWnc8UzpwoULcfToUSmULliwAIsXL5b6RTozcuRIdOjQAZGRkXjuueewdu1aLF68GAsWLMCHH36IyZMnY8mSJViyZIlD/yBXWLRoEdavX4+ff/5ZGjZp0qRSuxXYuwwAwIABA7B69erbLmPhwoUYMWKEw7CYmBiMHz8eI0eOLHU6+2XxGjVq4O+//8awYcPw5ZdfSpdn7wd6vR6FhYVOX4uOjkZqamqpN5i4k2XLlkmX3O9ESEiI0zN75dGgQQOsWbMGW7ZsAQAMHjwYRqMRRqPRIcgBKBHonP1XnnFuHu+VHxJxJffOa1Q9SIf/PdPsjqe3t8cuMzMTe/fuxX/+8x/Url0bS5cuxSeffIL33nuvQvN1tzPNRO5ApVJh9OjRmDt3Lg4ePIhHH30UBoMBgO2GpLNnz+Kzzz7DtGnT8O6772LTpk0AbP1ZbxYVFYWHH34YP/74IwDgmWeeQUBAAJYuXSrdB9G8eXM0b94cR44cQd++fWGxWFClShW88cYbCAgIwBdffHFv33wFuXUo9RaFhYWYMmUKpkyZcttxe/fujQcffBD9+vWThq1duxYA8OKLL+KBBx5Au3btpLOrs2fPRlxcHCZOnAir1YrHH39cuiTnzAMPPIDk5GQcPHgQFy5cwP/+9z8kJCRgypQpaNu2LbKzs9GoUSO0bNkSEydORH5+PqZNm4auXbuiSZMm+PPPP9G7d2/MmTMHXbp0QfXq1REaGorVq1dj2rRpeO211zB27Fjppp1bZWdnQ6/XQ6fTlXoJ4eDBg4iLi0OvXr0wceJEKJVKaLVa5OXlISgoCElJSXj77bel0DBs2DDodDpkZWVBoVBAq9Xi8uXLMJvNaNKkCeLi4gDYwtyUKVOkm5tu1blzZ+kO95u1adMGAQEB2LRpEzp27IiHHnoIc+bMAWDrm/nNN9/gwIEDGDJkCOrVq4e0tDQMGTIEGzduRJcuXaT5btmyBR07dkT9+vUd+pgFBQVh5MiRmDlzJgYOHIhvvvnGoYP8ggULMHLkSDRr1gwnTpzAjz/+iCeffBKffPIJFi5ciFOnTsHf3x/z589Ho0aNANhudNiwYQO2bNmCOXPmSI8RuRdatGgh3flvp9froVKppB0xgDsOlHfLaDRi6dKleOaZZ9CkSRPk5eVh4cKFaNGiBerVq3dP2rDgLgJlZdBqtWjbti0efPBBAMDAgQMxePBgmEwmh8vzt+Pn54c9e/aUuOny3LlzKCoqwqxZs5CYmOjStnsD9u+XV2hoKLRaLVJTUwHY9sm5ubl44IEHMH78eEyePLnEk0v+9a9/YcWKFSXmtWvXLmRnZ2PTpk1YsWIFXnvtNXTu3FmaBrB9Z+Xn5+O3335DfHw84uPj0b9/f+h0Onz44YcIDg6WHhnVrVs3vP322/jhhx/Qt29fCIKAgIAAHD16FJMmTcK8efNQvXp1TJgwAa1btwZg2//bnzQwcOBA6XFUM2bMuGf7uDvFPqX3wLhx47Br1y6nrw0fPhxbt25F//79Ub169TI75ZfX3r17MXr0aDRq1Ahjx45FrVq1cOLECQQGBqJ58+YlxrffXelOXFnHoqIibNu2Tbq7Oj8/H99//z2GDBkCg8EAPz8/tGrVCvPmzcPo0aMxbdo07N69G1u3bsXcuXPh7++PevXqoX379tJ6ys/Px5EjR9ChQ4cyl22/ez05ORk1a9bEoEGDpEd3HD16FPv27cPw4cNLnb53796YN28eatSoAavVitzcXIe79u1Ku3ve7tVXX8Xx48exbt06bN68GWFhYXj55Zexd+9efPHFFxg5cqT02du7dy/effdd7NixA/3798cXX3yBxx57DMOGDZOem6dSqWC1WhEYGIi8vDwAtvAcERGBiIgIjBw5EoIgYMSIEQgLC0PDhg2hUCggiqL0WJeqVavCYrFIfRDt87LfiAc4v6Hh5tcjIyOlx+FUq1YNV65cQc2aNaHX63H69GlpOfYz6ampqbh8+TJWrFiBFStWSOttxowZaNiwIXr37l1GNb3H559/jpCQEOlGp+XLlyM7Oxvjxo0DYPt8Dx48GN9++y00Gs0d3X1/O/YnXZSlY8eOOH/+PC5dunRHywCAsWPHIioqCm+++SaGDRsGURSxdOlSPP/88/jqq68AAC+99BJ+/vln/PDDD9Kj2G6+CfGVV16RnjwydepUvPfee1IXjCpVquD999/HCy+8UK722K/uDBs2DCNGjIBCoYAgCGjbti3q16+P7t27o2fPnujYsaM0zfvvv48333zzjtfB4sWL8eKLL0p3xJfm7bffxty5c6HT6dC+fXuMHDlSevrG2LFj8eeff8JoNKJXr174+eefsXjxYjRv3lw68Ads26FOp0PVqlWRmpqKNWvWYO/evRgzZozTZXbs2BFbtmxxuDnX3l5nQkNDMWbMGEyfPh0vvPAClixZUuH10a5dO0yYMAHPPPMM1q9fj6CgIKjValy+fBlRUVGwWCzo2LEjpk6dKu0TRFEscXNQRkYGCgoK8N577yE2Nhb9+/d3eFpJUVERNBrNHX2HmUwmfPTRRyVOZgmCgNTUVISGhkKv12PJkiXYunWr04BsZzAYMH78eCxevLjC7XCVcudIsRwKCgrE7t27iwUFBeUZXVZWq1VMTU0VrVar3E2RLF26VIyNjXX6X2V5+umnRbPZXK5x09LSKq0dd+pe19FsNouCIJT5+nPPPXdP2lIZduzYIc6cObPMcYqKisTz58+XGJ6dnS2azWbRaDSKv/zyixgbGysmJyeLP/30k2i1Wp1+jo8cOSKOHj1aPHfunHjy5EmHdWs0GkWj0Sj9nZmZKebm5opms1m0WCyi1WoV09LSxKysLFEUbes+KSlJ/Pvvv8X8/PwSy7py5YqYmZlZYrjBYCgxvtFoFC9evCj269dP3LNnjygIgnjp0iVx8ODB4pEjR8pcP97AYrGIJpNJ/PTTT8Xly5eLJpNJtFgsYnJysjhw4EDx3LlzotlsFhctWiROmTJFFEVRFARBLC4udqhhcnKyOH78eDE5Ofmu2jNp0iQxNjZW7NmzZ4l948qVK6Xl3/rauHHjxFGjRonLly8XBUEQ8/PzxQULFkiv//Of/xQXLFggnjhxwulys7Ozxc8//1zs3r17if3fmTNnpM/0lStXxIyMDPHatWviuHHjxNTUVFEURXHPnj1ibGysmJOT4zBthw4dxGeffVaa3r59xMbGigsWLBC7devmtD2CIIhJSUkOw7766itpPmazWYyNjZVqc+v6OH78eKnfMd98843Ulr///luMjY0VrVareO3aNTE/P1+8ePGi+MQTT4ixsbFl7gPtYmNjxS+++MJhuzcYDGJsbKzDNvTdd9+JU6dOldrfu3dvsbi4WHzuuedEs9ks/vbbb07nbZeTkyNmZ2eLn3zyifRe7PuNAwcOSOPa23Lrf0VFRdK/4+PjK/ydO2/ePLfKEXa3fjcWFxd7RDYrb45kKL1H7Ovu2LFjYmxsrPTlS865ax2pfG4OLbcGmjthNBpd+lk4dOiQOGbMGLFv377i0KFDxe+++85l83ZnK1euFLt37+7w3x9//CGKoihu375dHDZsmNi/f3/xnXfeETMyMkRRrNxQumDBAvGjjz4Sc3NzxdOnT4utWrUSRVEUBw4cKK5bt04ab9u2beLixYvF//u//xOHDx/u9OAkKSlJ/Omnn5y+5sy3334rjh492ulr9vdeGrPZLB47dqzE8O+//148ceKEFF5F0XaAZn9fFfHNN984DVHnz5+XApb9hIe9Dlu3bhVFURRHjRolxsbGir/88os03e32qcOGDStXu+bPny+mpKSIO3bscBi+ePFi6UDyTlUkYxQVFUn//vLLL6V1YrFYxJSUFFEURfGXX34Rv/zyS1EURfHrr7++q7a5C0/9bixvjuTle3JLrKNnu/nyblRUVIk7t8lziHf48Pw7Ye9OYbFYoFKpKvUzY//qc9fPpcViwcaNG/HUU0+VeM1sNuPhhx/G0qVLkZKSgo4dOzo8cWbLli04ffq0w82k3r5PTUlJQa9evdyuK5qreWodveKRUEREdP+4+ZeMKpu7hlE7tVrtNJACgI+PD9auXYuoqCg0bdq0xOsdO3Z06JN6P6hevbrcTSAXYCglIiLyMNWqVZO7CW7H/lQU8lyec+6XiIiIqBTt27eXuwl0lxhKiYiIiEh2DKVE5HL2voEWi0XmllBlsP8Yg0qlkrklRORN2KeUiFwuNDQUarUav//+Ox577DGHO4PJ89jvvhcEAZmZmVi3bh20Wi0iIiLkbhoReRGGUiJyOb1ejxdeeAGLFy9GUlKS9Jvv5HlEUZQe1WSvYd26dTFq1KgSPydKRHQ3GEqJqFLExMRg+vTpOHPmDEJDQxlKPZQoisjMzERYWBhUKhX8/f0REBDgUc9IJCLPwFBKRJVGr9cjKirK4x70TDcIggCNRsMaElGl4x6GiIiIiGTHUEpEREREsnPry/d/nErHj4mpFZpGFEUUFhZCr7/GPmwejHX0Dqyj52MNvQPr6B1cUcfZvRpDq3bPx7m5dShtWysEjaL8KzSNIIjIyMhAeHg4lEpueJ6KdfQOrKPnYw29A+voHVxRRx837hvu1qE0QKdGgK5iTRQEAT4mDSKDdOyU78FYR+/AOno+1tA7sI7ewdvr6H3viIiIiIg8DkMpEREREcmOoZSIiIiIZMdQSkRERESyYyglIiIiItkxlBIRERGR7BhKiYiIiEh2DKVEREREJDuGUiIiIiKSHUMpEREREcmOoZSIiIiIZMdQSkRERESyYyglIiIiItkxlBIRERGR7BhKiYiIiEh2arkbUJa0wnykFuZXaBpBEJCVl4VUtQVKJTO3p2IdvQPr6PlYQ+/AOnoHV9SxSXAUVG76GXDrUHo6LwMbU05VaBpRFFFQUAC/fD8oFIpKahlVNtbRO7COno819A6so3dwRR0bPhQJlZteKHfrUNq+SjTaV4mu0DSCICAtLQ2RkZE8GvRgrKN3YB09H2voHVhH7+DtdfS+d0REREREHoehlIiIiIhkx1BKRERERLJjKCUiIiIi2TGUEhEREZHsGEqJiIiISHYMpUREREQkO4ZSIiIiIpIdQykRERERyY6hlIiIiIhkx1BKRERERLJjKCUiIiIi2TGUEhEREZHsGEqJiIiISHYMpUREREQkO4ZSIiIiIpKdWu4GlCXv4A/I2bGsglOJMJmKUazVAFBUQqvo3mAdvQPr6PlYQ+/AOnqHu69jjVHfQqnRubZZLuLWoTQwrg8C4/pUaBpBEJCWlobIyEgolTwR7KlYR+/AOno+1tA7sI7ewdvr6H3viIiIiIg8DkMpEREREcmOoZSIiIiIZMdQSkRERESyYyglIiIiItkxlBIRERGR7BhKiYiIiEh2DKVEREREJDuGUiIiIiKSHUMpEREREcmOoZSIiIiIZMdQSkRERESyYyglIiIiItkxlBIRERGR7BhKiYiIiEh2DKVEREREJDu13A0oS+HVHBSkZFVoGlEQUZCfh8zLBVAoFZXUMqpsrKN3YB09H2voHVhH7+CKOoY1rw2Fyj3PSbp1KFWolVDpNBWaRhQFKIvVUOl9oFC450qn22MdvQPr6PlYQ+/AOnoHb6+jW4dSXXggdOGBFZpGEAQUp2kQHBkJpdL7Cna/YB29A+vo+VhD78A6egdvr6P3vSMiIiIi8jgMpUREREQkO4ZSIiIiIpIdQykRERERyY6hlIiIiIhkx1BKRERERLJjKCUiIiIi2TGUEhEREZHsGEqJiIiISHYMpUREREQkO4ZSIiIiIpIdQykRERERyY6hlIiIiIhkx1BKRERERLJjKCUiIiIi2TGUEhEREZHs1HI3oCzmwp0oNm6s2ESiCI3VhMIsLaBQVE7DqPKxjt6BdfR8rKF3YB29gwvq6Bs6DQqFxsUNcw23DqU++vbw0bev0DSCICA/LQ3BoZFQKnki2FOxjt6BdfR8rKF3YB29g7fX0fveERERERF5HIZSIiIiIpIdQykRERERyY6hlIiIiIhkx1BKRERERLJjKCUiIiIi2TGUEhEREZHsGEqJiIiISHYMpUREREQkO4ZSIiIiIpIdQykRERERyY6hlIiIiIhkx1BKRERERLJjKCUiIiIi2TGUEhEREZHsGEqJiIiISHZquRtQlmMZxTh4zVShaURRRKFRAX2GAQqFopJaRpWNdfQOrKPnYw29A+voHVxRx8Ex/vBRuednwK1Dad0gNar4qio0jSgIyMg0IjxMD4WSJ4I9FevoHVhHz8caegfW0Tu4oo4qNy6/W4dSXx8lfH0qNo0gKCAagAhfFZTc8DwW6+gdWEfPxxp6B9bRO3h7Hb3vHRERERGRx2EoJSIiIiLZMZQSERERkewYSomIiIhIdgylRERERCQ7hlIiIiIikh1DKRERERHJjqGUiIiIiGTHUEpEREREsmMoJSIiIiLZMZQSERERkewYSomIiIhIdgylRERERCQ7hlIiIiIikh1DKRERERHJTi13A8qSk5ODzMzMCk0jCAKysrKQn58PpZKZ21Oxjt6BdfR8rKF3YB29gyvqGB0d7bafAbcOpefPn8fu3bsrNI0oiigoKICfnx8UCkUltYwqG+voHVhHz8caegfW0Tu4oo4jRoxgKL0TLVu2RMuWLSs0jSAISEtLQ2RkpNuudLo91tE7sI6ejzX0Dqyjd/D2OnrfOyIiIiIij8NQSkRERESyYyglIiIiItkxlBIRERGR7BhKiYiIiEh2DKVEREREJDuGUiIiIiKSHUMpEREREcmOoZSIiIiIZMdQSkRERESyYyglIiIiItkxlBIRERGR7BhKiYiIiEh2DKVEREREJDuGUiIiIiKSHUMpEREREclOLXcDyvLD0VQsO3C5QtOIEFFsKoZGewkKKCqpZVTZWEfvwDp6PtbQO7CO3sEVdfz2uVjofFQubplruHUo7dOsKvo0q1qhaQRBQFpaGiIjI6FU8kSwp2IdvQPr6PlYQ+/AOnoHb6+j970jIiIiIvI4DKVEREREJDuvC6Vmsxlr166F2WyWuyl0F1hH78A6ej7W0Duwjt7B2+volaH0l19+8dqC3S9YR+/AOno+1tA7sI7ewdvr6HWhlIiIiIg8D0MpEREREcmOoZSIiIiIZCdbKF2/fr1ci75jldlmzvve8dT14anzriyeuj48dd6VxVPXh6fOu7J46vqorHl7Yg0B+dvNUFoBnrhhePK8K4unrg9PnXdl8dT14anzriyeuj48dd6VxVPXB0OpI7nbzcv3RERERCQ7hlIiIiIikh1DKRERERHJTl2ekURRBAAYjUaXLVgQBJfOz84+T6PRCKXStZm7strMeZfEOnrHvFlHz593ZdYQ8Lz14anzZh3v7bw9cZ8KVH677XmyNArxdmMAyMjIwNChQ13TMiIiIiK67yxduhTh4eGlvl6uUCoIArKysqDX66FQKFzaQCIiIiLyXqIoorCwEKGhoWWe4S1XKCUiIiIiqky80YmIiIiIZMdQSkRERESyYyglIiIiItmV65FQ90JhYSF++OEHnD59GqdPn4bBYMC4cePw+OOPO4x37NgxvPnmm07nMWvWLMTExGD27NlISEgodVnLli1DWFiY9LfZbMbKlSuxZcsWGAwG1K5dG4MHD0aLFi1c8+buI3LVsTzzo/JxZQ1vduXKFXz99dc4fvw48vPzERERgUceeQT//Oc/odPppPG4PbqGnHXk9ug6lVHHs2fPYsWKFThx4gQAoEGDBhg6dCjq1KnjMB23RdeRq46eti26TSjNy8vD6tWrERERgejoaBw7dqzM8Xv06IH69es7DKtatSoAoGvXrmjevLnDa6IoYv78+YiMjHQIpAAwZ84c7Nq1Cz179kS1atWwefNmTJ8+He+99x4aN25892/uPiJnHW83PyofV9bQLj09HRMmTICfnx+eeuopBAQE4OTJk1i1ahXOnTuHqVOnSuNye3QNuetY3nlS2Vxdx7Nnz+Lf//43wsPDMXDgQIiiiPXr12Py5Mn4+OOPUaNGDWlcbouuI2cdyzM/d+E2oTQ0NBTLly9HSEgIzpw5gwkTJpQ5fuPGjREfH+/0tZiYmBLpPykpCSaTCY8++qjD8NOnT2P79u0YOnQo+vTpAwDo1KkTRo8ejWXLlmHWrFl3/qbuQ3LVsTzzo/JxZQ3ttmzZgoKCAsycORO1atUCADz55JMQRREJCQkwGAzw9/fn9uhCctaxIvOksrm6jitXroRGo8GsWbMQGBgIAHj00UcxcuRILF++XDqrxm3RteSqY3nn5y7cpk+pj48PQkJCKjSN0WiE1Wot17jbtm2DQqHAI4884jB8165dUCqVePLJJ6VhGo0GnTt3xsmTJ5Genl6hNt3v5Krjnc6PSqqMGtp/zSM4ONhheEhICJRKJdRq2/Ext0fXkbOOFZknlc3VdUxKSsJDDz0kBRnAFpgaN26MAwcOoLCwEAC3RVeTq47lnZ+7cJszpRU1d+5cFBYWQqlUonHjxhg6dGiJU9N2FosFO3fuRExMDKpUqeLw2vnz51G9enX4+vo6DH/wwQcBABcuXEBERETlvAlyWR3vZH7kGuVZ502bNsX333+Pzz77DIMGDZIu+/7666/o3r271BeR26N8XFnHisyTXOt269xsNkOr1ZaYTqvVwmKx4NKlS4iJieG2KDNX1bG883MXHhdK1Wo12rVrh7i4OAQGBiI5ORk//vgjJk2ahA8//BB169YtMc3hw4eRn5/v9JJvVlaW06MX+7DMzEyXvwdyfR3vZH50dyqyzmNjYzF48GB8++232LdvnzS8X79++Ne//iX9ze3x3quMOnJ7vPfKu85r1KiBU6dOwWq1QqVSAbAFnNOnTwO4sY1xW5SHq+voaduix4XShg0bomHDhtLfbdq0QXx8PMaMGYPly5dj+vTpJabZtm0b1Go12rdvX+K14uJi+Pj4lBiu0Wik18n1XF3HO5kf3Z2KrvPIyEg0adIE7dq1Q0BAAA4ePIjvvvsOISEh6N69OwBuj3KojDpye7z3yrvOu3Xrhvnz5+PTTz/F008/DVEUsWbNGmRnZwO4sY1xW5SHq+voaduix4VSZ6pVq4a2bdti9+7dDkcNgO0xDPv27UOLFi0c+l7YaTQamM3mEsPtBbVvgFT57qaOFZ0fVY7S1vn27dsxb948LFy4EOHh4QCAdu3aQRAELFu2DB06dEBgYCC3Rzdxt3WsyDyp8jhb5127dkV6ejp+/PFH6ZF79erVQ58+ffDtt99KXTC4LbqPu6ljeefnLrwilAJAeHg4LBYLTCaTQx+YvXv3lnm3dmhoqNPLEPajDWePHaLKc6d1rOj8qPI4W+cbNmxA3bp1pSBj16ZNG2zevBnnz59H8+bNuT26kbupY0XmSZXL2Tp/7rnn0KdPH1y6dAl+fn6oXbs2li9fDgCoXr06AH43ups7rWNF5ucO3Obu+7t19epVaDSaEkcHW7duhV6vR+vWrZ1OFx0djZSUFOmuUrtTp05Jr9O9c6d1rOj8qPI4W+c5OTkQBKHEuBaLBQCkO0K5PbqPu6ljReZJlau0de7v74/GjRujdu3aAIA///wT4eHh0vMtuS26lzutY0XnJzePC6W5ubklhl24cAH79+9HixYtoFQqHcb966+/0LZt21JXfHx8PARBwG+//SYNM5vN2LRpExo0aMC7CyuJq+tYkfmRa1RknVerVg3nzp1DSkqKw/jbt2+HUqmUdqjcHu+9yqgjt8d7727W+Y4dO3DmzBn07NlTGo/bojxcXUdP2xbd6vL9unXrUFBQIF0y2L9/v/Tv7t27w8/PDx9++CE0Gg1iYmIQHByM5ORkbNy4EVqtFs8//7zD/Hbs2AGr1VrmJd8GDRogPj4ey5cvR25uLqpWrYqEhASkpaVh7NixlfZevZkcdazI/Oj2XF3DPn364NChQ5g0aZL0S0AHDhzAoUOH8MQTT0iXArk9upZcdeT26FqurGNiYiJWr16NFi1aICAgAKdOncKmTZvQsmVL9OzZUxqP26LryVFHT9sWFaIoinI3wu6FF15AWlqa09cWL16MKlWqYO3atdi2bRtSU1NhNBoRFBSEZs2aYeDAgahWrZrDNK+99hquXbuGZcuWldmRt7i4GF9//TW2bt3q8Pu+LVu2dOn7u1/IUceKzI9uz9U1BGy/ELNq1SqcP38e+fn5qFKlCjp16oSnn37aoa7cHl1Hrjpye3QtV9YxNTUV//vf/3Du3DkUFhZK9evdu3eJu+25LbqWHHX0tG3RrUIpEREREd2f3KszARERERHdlxhKiYiIiEh2DKVEREREJDuGUiIiIiKSHUMpEREREcmOoZSIiIiIZMdQSkRERESyYyglIiIiItm51c+MEpH36dGjR4XGj4yMxJIlSzB58mQkJiZKv3TiCTZt2oS5c+dKf+t0Onz33XfS39euXcOLL76IJk2a4IMPPnDZcmfPno2EhAS8//77aNq0abmmOXbsGN5880106tQJr776qjT8559/xuLFi6W/7fUgIqpsDKVEVKk6depUYtiJEyeQmpqK6OhoREdHO7wWGBh4r5pWaezvS6PRuGR+9p8n/OWXX1wyv7LUrFlTqllCQkKlL4+IyI6hlIgq1c1n4exmz56N1NRUtG3bFoMGDSp1OpPJhLCwsMpuosuV9b4qw/PPP49nnnkGERERdz2vli1bSr9tzlBKRPcSQykRuaXIyEi5m+AxQkNDERoaKncziIjuCkMpEbml0vqU9ujRA5GRkVi0aBG+++47JCQkIDMzE5GRkXj66afx+OOPAwD++usvrFmzBmfPnoVSqUTr1q3x4osvOu0eYLVasXHjRiQkJCA5ORlWqxXVq1fHY489hu7du0OlUrn8/RmNRnz99dfYs2cPcnNzUaVKFXTp0gU9e/aEUmm7B9Xe7/Pm9253c1/PsvqUXrp0CStWrEBiYiIEQUB0dDT69evnsq4FRESuwlBKRB5p5syZOHr0KJo2bYqoqCgkJiZKNxnp9XrMmjULDRo0QMuWLXHy5Els2bIF165dw4wZM6BQKKT5mEwm/Oc//8HRo0cREBCABg0aQKPR4PTp01i8eLEUDO1B0RXMZjOmTJmC1NRUNGvWDBaLBX/99ReWLFmCCxcuSF0egoOD0alTJ+zevRtFRUUO/XPL0/f2zJkzmDJlCgoLC1GrVi3UqlULV65cwfTp09G1a1eXvR8iIldgKCUij5OWlga9Xo+FCxciKCgIAHD06FFMmTIFK1askEJfq1atANjOSr7++us4fvw4jh07hmbNmknz+vLLL3H06FH84x//wKhRo+Dn5ydNM2vWLOzbtw8bN250aYg7deoUateu7dD+1NRUTJo0CQkJCWjbti0efvhh1KxZE6+++ioSExNRVFTktH9uaURRxOzZs1FYWIgBAwbg2WeflV5bv349FixY4LL3Q0TkCnxOKRF5pOHDh0uBDgCaNWuGOnXqICsrC7GxsVIgBQBfX1906dIFAJCYmCgNz8nJwe+//47w8HCMGzdOCqT2acaOHQu1Wo0NGza4vP3Dhg1zaH/VqlUxYMAAALbQeLeOHTuGy5cvIyoqSpqv3VNPPYUGDRrc9TKIiFyJoZSIPI5arUaTJk1KDI+KigIAtGjRotTXsrKypGHHjh2DxWJBbGwstFptiWlCQkJQrVo1XLp0CSaTyVXNR0BAgNM2dujQAYDtkVmCINzVMpKSkgAA8fHxTvvE2pdFROQuePmeiDxOcHCw06Cl0+kAwOljpPR6PQBbf067tLQ0AMDGjRuxcePGMpdpMBicBtc7Udqjm/z8/ODn54eCggIYDIa7emarPXyXtiw+3YCI3A1DKRF5nNvddFTem5LsZyPr1KmD2rVrlzmuWs3dJRFRZeJelojuW+Hh4QCARo0aYcSIEfdsuenp6U6HG41GFBQUQKPROPRvvRP255aWtiz7WWIiInfBPqVEdN9q1qwZlEol9u/fD4vFcs+Wm5+fj7/++qvE8O3btwMAYmJiHLon2M/SWq3Wci+jUaNGAIDdu3c77Z+6Y8eOCrWZiKiyMZQS0X0rLCwMnTt3RlpaGmbNmoXs7OwS41y5cgW7du1y+bK//PJL5OXlSX9fvXoVq1evBmC7O/5m9rOef//9d7nn37RpU9SoUQOpqalYs2aNw2u//vorTp48eadNJyKqFLx8T0T3teHDh+PatWvYvXs3Dh8+jOjoaERERMBkMiE5ORmpqalo06YN4uPjXbbMBg0awGKx4KWXXkKzZs1gtVrx119/wWQy4dFHH0W7du0cxm/Tpg0SExMxdepUNGvWDFqtFoGBgRgyZEipy1AqlRg/fjymTp2KVatWYffu3XjggQeQmpqKs2fPolu3bpXyqCsiojvFUEpE9zWtVot33nkH27Ztw+bNm3HhwgWcOXMGgYGBiIyMRMeOHV3++CQfHx9Mnz4dy5cvx969e5GXl+fwM6O36tGjBwwGA7Zv347du3fDYrEgMjKyzFAK2MLvrFmzsGLFCiQlJeHq1auoXbs2pk2bBp1Ox1BKRG5FIYqiKHcjiIi8waZNmzB37lwMHDgQgwYNumfLnTlzJnbu3ImPP/4YDz74oMvm26NHD0RGRmLJkiUumycRUWl4ppSIyMX27t2La9euQaPRYNSoUZW6LKvVigsXLkChUEg/EHA3Dh8+jG3btrmgZUREFcNQSkTkYhcuXMCFCxeg0+kqNZTOmzcPJ06cQEpKClq2bHlXD9u3u3z5MhISElzQOiKiiuHleyIiD2X/TfsWLVpgxIgRCA4OlrdBRER3gaGUiIiIiGTH55QSERERkewYSomIiIhIdgylRERERCQ7hlIiIiIikh1DKRERERHJjqGUiIiIiGTHUEpEREREsmMoJSIiIiLZMZQSERERkez+H8QWSRyL1gx+AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 848.5x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cbvCorrector.cbvs[0].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first several CBVs contain most of the systematics. The latter CBVs pose a greater risk of injecting more noise than helping. The default behavior in CBVCorrector is to use the first 8 CBVs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assessing Over- and Under-Fitting with the Goodness Metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two very common issues when fitting a model to a data set it over-fitting and under-fitting. \n",
    "\n",
    "Over-fitting occurs when the model has too many degrees of freedom and fits the data  _at all costs_, instead of just modelling the physical process it is attempting to model. This can exhibit itself in different ways, depending on the system and application. In the case of fitting systematic trend basis vectors to a time series, over-fitting can result in the basis vectors removing intrinsic signals in the times series instead of just the systematics. It can also result in introduced broad-band noise. This can be particularily prominant in an unconstrained least-squares fit. A least-squares fit only cares about minimizing it's loss function, which is the Root Mean Square error. In the course of minimizing the RMS, narrow-band power representing the systematic trends are exchanged for broad-band noise intrinsic in the basis vectors. This results in the overall RMS decreasing but the noise in the time series increasing, resulting in the obscuration of the signals under interest. A very common method to inhibit over-fitting is to introduce a regularization term in the loss function. This constrains the fit and effectively reduces the degrees of freedom.\n",
    "\n",
    "Under-fitting occurs when the model has too few degrees of freedom and fails to adequately model the physical process it is attempting to model. In the case of fitting systematic trend basis vectors to a time series, under-fitting can result in residual systematics. Under-fitting can either be the result of the basis vectors not adequately representing the systematics or, placing too great of a restriction on the model during fitting. The regularization technique used to inhibit over-fitting can therefore result in under-fitting. The ideal fit will balance the counter-acting phenomena of over- and under-fitting. To this end, a method can be developed to measure the degree to which these two phenomena occur.\n",
    "\n",
    "PDC has two **Goodness Metrics** to assess over- and under-fitting:\n",
    "\n",
    "- **Over-fitting metric**: Measures the introduced noise in the light curve after the correction. It does so by measuring the broad-band power spectrum via a Lomb-Scargle Periodogram both before and after the correction. If power has increased after the correction then this is an indication the CBV fit has over-fitted and introduced noise. The metric treats all frequencies equally when measuring power increase; from one frequency separation to the Nyquist frequency. This metric is callibrated such that a metric value of 0.5 means the introduced noise due to over-fitting is at the same power level as the uncertainties in the light curve.\n",
    "\n",
    "- **Under-fitting metric**: Measures the mean residual target to target Pearson correlation between the target under study and a selection of neighboring targets. This metric will find and download a selection of neighboring SPOC SAP targets in RA and Decl. until a minimum number is found. The metric is callibrated such that a value of 0.95 means the residual correlations in the target is equivalent to chance correlations of White Gaussian Noise.\n",
    "\n",
    "_The Goodness Metrics are not perfect!_ They are an estimate of over- and under-fitting and are to be used as a guideline along other other metrics to assess the quality of your light curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Goodness Metrics are part of the `lightkurve.correctors.metrics` module and can be computed directly with calls to `overfit_metric_lombscargle` and `underfit_metric_neighbors`. The 'CBVCorrector' has convenience wrappers for the two metrics and so they do not need to be called directly, as we will show below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Correction with CBVCorrector to Inhibit Over-Fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are four correction methods within `CBVCorrector`:\n",
    "\n",
    "- **correct**: Performs a numerical correction using the LightKurve [RegressionCorrector.correct](https://docs.lightkurve.org/reference/api/lightkurve.correctors.RegressionCorrector.correct.html?highlight=regressioncorrector%20correct#lightkurve.correctors.RegressionCorrector.correct) method while optimizing the L2-Norm regularization penalty term using the goodness metrics as the Loss Function.\n",
    "\n",
    "- **correct_gaussian_prior**: Performs an analytical correction using the LightKurve [RegressionCorrector.correct](https://docs.lightkurve.org/reference/api/lightkurve.correctors.RegressionCorrector.correct.html?highlight=regressioncorrector%20correct#lightkurve.correctors.RegressionCorrector.correct) method while setting the L2-Norm (Ridge Regression) regularization penalty term as the Gaussian prior.\n",
    "\n",
    "- **correct_elasticnet**:  Performs the correction using Scikit-Learn's [ElasticNet](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html), which combines both an L1- and L2-Norm regularization.\n",
    "\n",
    "- **correct_regressioncorrector**: Performs the standard [RegressionCorrector.correct](https://docs.lightkurve.org/reference/api/lightkurve.correctors.RegressionCorrector.correct.html?highlight=regressioncorrector%20correct#lightkurve.correctors.RegressionCorrector.correct) correction. RegressionCorrector is the superclass to CBVCorrector and this is a \"passthrough\" method to access the superclass `correct` method.\n",
    "\n",
    "If you are unfamilar with L1-Norm (LASSO) and L2-Norm (Ridge Regression) regularization then there are [severel](https://www.statlearning.com/), [excellent](https://press.princeton.edu/books/hardcover/9780691198309/statistics-data-mining-and-machine-learning-in-astronomy), [introductions](https://en.wikipedia.org/wiki/Regularization_(mathematics)). You can read how a L2-norm relates to a Gaussian prior in a linear design matrix in [this reference](https://katbailey.github.io/post/from-both-sides-now-the-math-of-linear-regression/).\n",
    "\n",
    "The default method is `correct` and generally speaking, one can use just this method to obtain a good fit. The other methods are for advanced usage.\n",
    "\n",
    "We'll start with `correct_gaussian_prior` in order to introduce the concepts. Doing so will allow us to force a very weak regularization term (alpha=1e-4) as an illustration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:09.836910Z",
     "iopub.status.busy": "2023-11-03T14:15:09.836797Z",
     "iopub.status.idle": "2023-11-03T14:15:09.920700Z",
     "shell.execute_reply": "2023-11-03T14:15:09.919820Z"
    }
   },
   "outputs": [],
   "source": [
    "# Select which CBVs to use in the correction\n",
    "cbv_type = ['SingleScale', 'Spike']\n",
    "# Select which CBV indices to use\n",
    "# Use the first 8 SingleScale and all Spike CBVS\n",
    "cbv_indices = [np.arange(1,9), 'ALL']\n",
    "# Perform the correction\n",
    "cbvCorrector.correct_gaussian_prior(cbv_type=cbv_type, cbv_indices=cbv_indices, alpha=1e-4)\n",
    "cbvCorrector.diagnose();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First note that CBVCorrector always fits a constant term in the model, but the constant is never subtracted in the resultant corrected flux. The median flux value of the light curve is always preserved.\n",
    "\n",
    "At first sight, this looks like a good correction. Both the Single-Scale and Spike basis vectors are being utilized to fit out as much of the signal as possible. The corrected light curve is indeed flatter. But this was essentially an unrestricted least-squares correction and we may have _over-fitted_. The very strong lips right at the beginning of each orbit is probably a chance correlation between the star's inherent stellar variability and the thermal settling systematic that is common in Kepler and TESS lightcurves. Let's look at the CBVCorrector goodness metrics to determine if this is the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:09.927651Z",
     "iopub.status.busy": "2023-11-03T14:15:09.926669Z",
     "iopub.status.idle": "2023-11-03T14:15:38.242232Z",
     "shell.execute_reply": "2023-11-03T14:15:38.241815Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Over fitting Metric: 0.7126149574495612\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Under fitting Metric: 0.9703478826925148\n"
     ]
    }
   ],
   "source": [
    "# Note: this cell will be slow to run\n",
    "print('Over fitting Metric: {}'.format(cbvCorrector.over_fitting_metric()))\n",
    "print('Under fitting Metric: {}'.format(cbvCorrector.under_fitting_metric()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first time you run the under-fitting goodness metric it will download the SPOC SAP light curves of targets in the neighborhood around the target under study in order to estimate the residual systematics (they are stored in the CBVCorrector object for subsequent computations). A goodness metric of 0.8 or above is generally considered good. In this case, it looks like we over-fitted (over-fitting metric = 0.71). Even though the corrected light curve looks better, our metric is telling us we probably injected signals into our lightcurve and we should not trust this really nice looking curve. Perhaps we can do better if we _regularize_ the fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the Goodness Metrics to Optimize the Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start by performing a scan of the over- and under-fit goodness metrics as a function of the L2-Norm regularization term, alpha."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:15:38.244137Z",
     "iopub.status.busy": "2023-11-03T14:15:38.244018Z",
     "iopub.status.idle": "2023-11-03T14:16:03.742560Z",
     "shell.execute_reply": "2023-11-03T14:16:03.739670Z"
    }
   },
   "outputs": [],
   "source": [
    "cbvCorrector.goodness_metric_scan_plot(cbv_type=cbv_type, cbv_indices=cbv_indices);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This scan also plots the last used alpha parameter as a vertical black line (alpha=1e-4 in our case). We are clearly not optimizing this fit for both over- and under-fitting. Let's use the `correct` numerical optimizer to try to optimize the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:16:03.748738Z",
     "iopub.status.busy": "2023-11-03T14:16:03.748591Z",
     "iopub.status.idle": "2023-11-03T14:16:38.922976Z",
     "shell.execute_reply": "2023-11-03T14:16:38.920445Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized Over-fitting metric: 0.9996578702189283\n",
      "Optimized Under-fitting metric: 0.8255781316730026\n",
      "Optimized Alpha: 9.656e+03\n"
     ]
    }
   ],
   "source": [
    "cbvCorrector.correct(cbv_type=cbv_type, cbv_indices=cbv_indices);\n",
    "cbvCorrector.diagnose();\n",
    "cbvCorrector.goodness_metric_scan_plot(cbv_type=cbv_type, cbv_indices=cbv_indices);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better! We see the thermal settling systematic is still being removed, but the stellar variabity is better preserved. Note that the optimizer did not set the alpha parameter at exactly the red and blue curve intersection point. The default target goodness scores is 0.8 or above, which is fulfilled at alpha=1.45e-1. If we want to optimize the fit even more, by perhaps ensuring we are not over-fitting at all, then we can adjust the target over and under scores to emphasize which metric we are more interested in. Below we more greatly emphasize improving the over-fitting metric by setting the target to 0.9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:16:38.930597Z",
     "iopub.status.busy": "2023-11-03T14:16:38.930172Z",
     "iopub.status.idle": "2023-11-03T14:17:13.802163Z",
     "shell.execute_reply": "2023-11-03T14:17:13.801188Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized Over-fitting metric: 0.9996569267677275\n",
      "Optimized Under-fitting metric: 0.8255781330421909\n",
      "Optimized Alpha: 9.662e+03\n"
     ]
    }
   ],
   "source": [
    "cbvCorrector.correct(cbv_type=cbv_type,\n",
    "                     cbv_indices=cbv_indices, \n",
    "                     target_over_score=0.9,\n",
    "                     target_under_score=0.5)\n",
    "cbvCorrector.diagnose();\n",
    "cbvCorrector.goodness_metric_scan_plot(cbv_type=cbv_type, cbv_indices=cbv_indices);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now perhaps biasing too far towards under-fitting, but depending on your research interests, this might be best."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## No Single Best Answer Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at another example, this time where there is no clear single best answer. Again, we will use the Single-Scale and Spike basis vectors for the correction and begin with low regularization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:17:13.809780Z",
     "iopub.status.busy": "2023-11-03T14:17:13.809615Z",
     "iopub.status.idle": "2023-11-03T14:17:37.107344Z",
     "shell.execute_reply": "2023-11-03T14:17:37.106381Z"
    }
   },
   "outputs": [],
   "source": [
    "lc = search_lightcurve('TIC 38574307', author='SPOC', sector=2).download(flux_column='sap_flux')\n",
    "cbvCorrector = CBVCorrector(lc)\n",
    "cbv_type = ['SingleScale', 'Spike']\n",
    "cbv_indices = [np.arange(1,9), 'ALL']\n",
    "cbvCorrector.correct_gaussian_prior(cbv_type=cbv_type, cbv_indices=cbv_indices, alpha=1e-4)\n",
    "cbvCorrector.diagnose();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first sight, this looks good. The long term trends have been removed and the periodic noisy bits have been removed with the spike basis vectors. But did we really do a good job?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:17:37.126914Z",
     "iopub.status.busy": "2023-11-03T14:17:37.126296Z",
     "iopub.status.idle": "2023-11-03T14:20:27.187957Z",
     "shell.execute_reply": "2023-11-03T14:20:27.179389Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Over fitting Metric: 0.2633104475330727\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Under fitting Metric: 0.9995877972880245\n"
     ]
    }
   ],
   "source": [
    "print('Over fitting Metric: {}'.format(cbvCorrector.over_fitting_metric()))\n",
    "print('Under fitting Metric: {}'.format(cbvCorrector.under_fitting_metric()))\n",
    "cbvCorrector.goodness_metric_scan_plot(cbv_type=cbv_type, cbv_indices=cbv_indices);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm... The over-fitting goodness metric says we are severely over-fitting. Not only that, there appears to not be an Alpha parameter that brings both goodness metrics above 0.8. And yet, the fit looks really good. What's going on here? Let's zoom in on the correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:20:27.199372Z",
     "iopub.status.busy": "2023-11-03T14:20:27.198968Z",
     "iopub.status.idle": "2023-11-03T14:20:27.256038Z",
     "shell.execute_reply": "2023-11-03T14:20:27.253013Z"
    }
   },
   "outputs": [],
   "source": [
    "pltAxis = cbvCorrector.diagnose()\n",
    "pltAxis[0].set_xlim(1360.5, 1361.1)\n",
    "pltAxis[0].set_ylim(1.4755e7, 1.477e7);\n",
    "pltAxis[1].set_xlim(1360.5, 1361.1)\n",
    "pltAxis[1].set_ylim(1.475e7, 1.477e7);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see in the top plot that the _SingleScale_ correction has comperable noise to the _original_ light curve. This means the correction is injecting high frequency noise at comperable amplitude to the original signal. We have indeed over-fitted! The goodness metrics perform a _broad-band_ analysis of over- and under-fitting. Even though our eyes did not see the high frequency noise injection, the goodness metrics did. So, what should be done? It depends on what you are trying to investigate. If you are only looking at the low frequency signals in the data then perhaps you don't care about the high frequency noise injection. If you really do care about the high frequency signals then you should increase the Alpha parameter, or set the target goodness scores as we do below (target_over_score=0.8, target_under_score=0.5). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:20:27.258415Z",
     "iopub.status.busy": "2023-11-03T14:20:27.258258Z",
     "iopub.status.idle": "2023-11-03T14:21:01.475436Z",
     "shell.execute_reply": "2023-11-03T14:21:01.473248Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized Over-fitting metric: 0.8001426013816765\n",
      "Optimized Under-fitting metric: 0.4687789542372483\n",
      "Optimized Alpha: 3.340e+00\n"
     ]
    }
   ],
   "source": [
    "# Optimize the fit but overemphasize the importance of not over-fitting.\n",
    "cbvCorrector.correct(cbv_type=cbv_type,\n",
    "                     cbv_indices=cbv_indices, \n",
    "                     target_over_score=0.8,\n",
    "                     target_under_score=0.5)\n",
    "cbvCorrector.diagnose()\n",
    "cbvCorrector.goodness_metric_scan_plot(cbv_type=cbv_type, cbv_indices=cbv_indices);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:01.484240Z",
     "iopub.status.busy": "2023-11-03T14:21:01.483891Z",
     "iopub.status.idle": "2023-11-03T14:21:01.609286Z",
     "shell.execute_reply": "2023-11-03T14:21:01.608977Z"
    }
   },
   "outputs": [],
   "source": [
    "# Again, zoom in to see the detail\n",
    "pltAxis = cbvCorrector.diagnose()\n",
    "pltAxis[0].set_xlim(1360.5, 1361.1)\n",
    "pltAxis[0].set_ylim(1.4755e7, 1.477e7);\n",
    "pltAxis[1].set_xlim(1360.5, 1361.1)\n",
    "pltAxis[1].set_ylim(1.475e7, 1.477e7);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see now that the high frequency noise injection is small compared to the original amplitudes in the lightkurve. We barely removed any of the systematics and are now under-fitting, but that might be the best we can do if we want to ensure low noise injection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ...Or Can We Still Do Better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps we are using the incorrect CBVs for this target. Below is a tuned multi-step fit where we first fit the multi-scale Band 2 CBVs then the Spike CBVs. The multi-scale band 2 CBVs contain intermediate frequency systematic signals. They should not inject high frequency noise. We also utilize the `correct_elasticnet` corrector, which allows us to add in a L1-Norm term (Lasso Regularization). L1-Norm helps snap some basis vector fit coefficients to zero and can result in a more stable, less noisy fit. The result is a much better compromise between over- and under-fitting. The spikes are not well removed but increasing the weight on the Spike CBV removal results in over-fitting. We can also try the multi-scale band 3 CBVs, which contain high frequency systematics, but the over-fitting metric indicates using them results in even greater over-fitting. The resultant is now much better than what we achieved above but more tuning and optimization could possibly get us even closer to an ideal fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:01.611300Z",
     "iopub.status.busy": "2023-11-03T14:21:01.611177Z",
     "iopub.status.idle": "2023-11-03T14:21:02.133939Z",
     "shell.execute_reply": "2023-11-03T14:21:02.123877Z"
    }
   },
   "outputs": [],
   "source": [
    "# Fit to the Multi-Scale Band 2 CBVs with ElasticNet to add in a L1-Norm (Lasso) term\n",
    "cbvCorrector.correct_elasticnet(cbv_type=['MultiScale.2'], cbv_indices=[np.arange(1,9)], alpha=1.0e-7, l1_ratio=0.5)\n",
    "ax = cbvCorrector.diagnose()\n",
    "ax[0].set_title('Result of First Correction to MultiScale.2 CBVs');\n",
    "\n",
    "# Set the corrected LC as the initial LC in a new CBVCorrector object before moving to the next correction.\n",
    "# You could instead just reassign to the first cbvCorrector object, if you do not wish to save the original.\n",
    "cbvCorrectorIter2 = cbvCorrector.copy()\n",
    "cbvCorrectorIter2.lc = cbvCorrectorIter2.corrected_lc.copy()\n",
    "\n",
    "# Fit to the Spike Basis Vectors, using an L1-Norm term.\n",
    "cbvCorrectorIter2.correct_elasticnet(cbv_type=['Spike'], cbv_indices=['ALL'], alpha=2.0e-5, l1_ratio=0.7)\n",
    "ax = cbvCorrectorIter2.diagnose()\n",
    "ax[0].set_title('Result of Second Correction to Spike CBVs');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:02.143721Z",
     "iopub.status.busy": "2023-11-03T14:21:02.143411Z",
     "iopub.status.idle": "2023-11-03T14:21:02.989802Z",
     "shell.execute_reply": "2023-11-03T14:21:02.989314Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Over-fitting Metric: 0.6262886935335121\n",
      "Under-fitting Metric: 0.880152175383505\n"
     ]
    }
   ],
   "source": [
    "# Compute the final goodness metrics compared to the original lightcurve.\n",
    "# This requires us to copy the original light curve into cbvCorrectorIter2.lc so that the goodness metrics compares the corrected_lc to the proper initial light curve.\n",
    "cbvCorrectorIter2.lc = cbvCorrector.lc.copy()\n",
    "print('Over-fitting Metric: {}'.format(cbvCorrectorIter2.over_fitting_metric()))\n",
    "print('Under-fitting Metric: {}'.format(cbvCorrectorIter2.under_fitting_metric()))\n",
    "\n",
    "# Plot the final correction\n",
    "_, ax = plt.subplots(1, figsize=(10, 6))\n",
    "cbvCorrectorIter2.lc.plot(ax=ax, normalize=False, alpha=0.2, label='Original')\n",
    "cbvCorrectorIter2.corrected_lc[~cbvCorrectorIter2.cadence_mask].scatter(\n",
    "                                            normalize=False, c='r', marker='x',\n",
    "                                            s=10, label='Outliers', ax=ax)\n",
    "cbvCorrectorIter2.corrected_lc.plot(normalize=False, label='Corrected', ax=ax, c='k')\n",
    "ax.set_title('Comparison between original and final corrected lightcurve');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, which CBVs are best to use? There is no one single answer, but generally speaking, the Multi-Scale Basis vectors are more versatile. The trade-off is there are also more of them, which means more degrees of freedom in your fit. More degrees of freedom can result in more over-fitting without proper regularization. It is recommened the user tries different combinations of CBVs and use objective metrics to decide which fit is the best for their particular needs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the Goodness Metrics and CBVCorrector with other Design Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Goodness Metrics and CBVCorrector can also be used in conjunction with other external design matrices. Let's work on a famous planet example to show how the CBVCorrector can be utilized to imporove the generated light curve. We will begin by using [search_tesscut](https://docs.lightkurve.org/reference/api/lightkurve.search_tesscut.html?highlight=search_tesscut) to extract an FFI light curve for HAT-P 11 and then create a DesignMatrix using the background pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:02.991660Z",
     "iopub.status.busy": "2023-11-03T14:21:02.991539Z",
     "iopub.status.idle": "2023-11-03T14:21:10.865241Z",
     "shell.execute_reply": "2023-11-03T14:21:10.860541Z"
    }
   },
   "outputs": [],
   "source": [
    "# HAT-P 11b\n",
    "from lightkurve import search_tesscut\n",
    "from lightkurve.correctors import DesignMatrix\n",
    "search_result = search_tesscut('HAT-P-11', sector=14)\n",
    "tpf = search_result.download(cutout_size=20)\n",
    "# Create a simple thresholded aperture mask\n",
    "aper = tpf.create_threshold_mask(threshold=15, reference_pixel='center')\n",
    "# Generate a simple aperture photometry light curve\n",
    "raw_lc = tpf.to_lightcurve(aperture_mask=aper)\n",
    "# Create a design matrix using PCA components from the cutout background\n",
    "dm = DesignMatrix(tpf.flux[:, ~aper], name='pixel regressors').pca(5).append_constant()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [DesignMatrix](https://docs.lightkurve.org/reference/api/lightkurve.correctors.DesignMatrix.html?highlight=designmatrix#lightkurve.correctors.DesignMatrix) `dm` now contains the common trends in the background pixels in the data. We will first try to fit the pixel-based design matrix using an unrestricted least-squares fit (I.e. a very weak regularization by setting alpha to a small number). We tell CBVCorrector to only use the external design matrix with `ext_dm=`. When we generate the CBVCorrector object the CBVs will be downloaded, but the CBVs are for 2-minute cadence and not the 30-minute FFIs. We therefore use the `interpolate_cbvs=True` option to tell the CBVCorrector to interpolate the CBVs to the light curve cadence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:10.873319Z",
     "iopub.status.busy": "2023-11-03T14:21:10.872832Z",
     "iopub.status.idle": "2023-11-03T14:21:12.352954Z",
     "shell.execute_reply": "2023-11-03T14:21:12.352644Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Over-fitting metric: 0.08088984759691169\n",
      "CDPP: 230.75800309519127 ppm\n"
     ]
    }
   ],
   "source": [
    "# Generate the CBVCorrector object and interpolate the downloaded CBVs to the light curve cadence\n",
    "cbvcorrector = CBVCorrector(raw_lc, interpolate_cbvs=True)\n",
    "# Perform an unrestricted least-squares fit using only the pixel-derived design matrix.\n",
    "cbvcorrector.correct_gaussian_prior(cbv_type=None, cbv_indices=None, ext_dm=dm, alpha=1e-4)\n",
    "cbvcorrector.diagnose()\n",
    "print('Over-fitting metric: {}'.format(cbvcorrector.over_fitting_metric()))\n",
    "print('CDPP: {}'.format(cbvcorrector.corrected_lc.estimate_cdpp()))\n",
    "corrected_lc_just_pixel_dm = cbvcorrector.corrected_lc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The least-squares fit did remove the background flux trend and at first sight the resultant might look good, but the over-fitting goodness metric is `0.08`. That's not very good! It looks like we are dramatically over-fitting. We can see this in the bottom plot where the corrected curve has more high-frequency noise than the original. Let's now add in the multi-scale basis vectors and see if we can do better. Note that we are joint fitting the CBVs and the external pixel-derived design matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:12.354648Z",
     "iopub.status.busy": "2023-11-03T14:21:12.354559Z",
     "iopub.status.idle": "2023-11-03T14:21:12.583521Z",
     "shell.execute_reply": "2023-11-03T14:21:12.581823Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Over-fitting metric: 0.3103857337957571\n",
      "CDPP: 55.83742644726604 ppm\n"
     ]
    }
   ],
   "source": [
    "cbv_type = ['MultiScale.1', 'MultiScale.2', 'MultiScale.3','Spike']\n",
    "cbv_indices = [np.arange(1,9), np.arange(1,9), np.arange(1,9), 'ALL']\n",
    "cbvcorrector.correct_gaussian_prior(cbv_type=cbv_type, cbv_indices=cbv_indices, ext_dm=dm, alpha=1e-4)\n",
    "cbvcorrector.diagnose()\n",
    "print('Over-fitting metric: {}'.format(cbvcorrector.over_fitting_metric()))\n",
    "print('CDPP: {}'.format(cbvcorrector.corrected_lc.estimate_cdpp()))\n",
    "corrected_lc_joint_fit = cbvcorrector.corrected_lc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks a lot better! Could we do a bit better by adding in a regularization term? Let's do a goodness metric scan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:21:12.590218Z",
     "iopub.status.busy": "2023-11-03T14:21:12.589121Z",
     "iopub.status.idle": "2023-11-03T14:22:03.612462Z",
     "shell.execute_reply": "2023-11-03T14:22:03.611369Z"
    }
   },
   "outputs": [],
   "source": [
    "cbvcorrector.goodness_metric_scan_plot(cbv_type=cbv_type, cbv_indices=cbv_indices, ext_dm=dm);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a couple observations to make here. \n",
    "\n",
    "First, the under-fitting metric has a very good score throughout the regularization scan. This is because the under-fitting metric compares the corrected light curve to neighboring targets in RA and Decl. that are archived as 2-minute SAP-flux targets. _The SAP flux already has the background removed_ so the neighoring targets do not contain the very large background trends. The under-fitting metric is therefore not very helpful. In the next run we will disable the under-fitting metric in the optimization (by setting target_under_score=-1).\n",
    "\n",
    "We see the over-fitting metric is not a simple function of the regularization factor _alpha_. This can happen due to the interaction of the various basis vectors during fitting when regularization is applied. We see a minima (most over-fitting) at about alpha=1e-1. Once alpha moves above this value we begin to over-constrain the fit which results in gradually less removal of systematics. The under-fitting metric should be an indicator that we are going too far constraining the fit, and indeed, we do see the under-fitting metric degrades slightly beginning at alpha=1e-1.\n",
    "\n",
    "We will now try to optimize the fit and account for these two issues by 1) setting the bounds on the alpha parameter (`alpha_bounds=[1e-6, 1e-2]`) and 2) disregarding the under-fitting metric (`target_under_score=-1`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:22:03.621815Z",
     "iopub.status.busy": "2023-11-03T14:22:03.621319Z",
     "iopub.status.idle": "2023-11-03T14:22:04.499880Z",
     "shell.execute_reply": "2023-11-03T14:22:04.498872Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimized Over-fitting metric: 0.34981996985663844\n",
      "Optimized Alpha: 4.436e-04\n",
      "CDPP: 60.654244871847794 ppm\n"
     ]
    }
   ],
   "source": [
    "# Optimize the fit but ignore the under-fitting metric and set bounds on the alpha parameter.\n",
    "cbvcorrector.correct(cbv_type=cbv_type, cbv_indices=cbv_indices, ext_dm=dm, alpha_bounds=[1e-6, 1e-2], target_over_score=0.8, target_under_score=-1)\n",
    "cbvcorrector.diagnose();\n",
    "print('CDPP: {}'.format(cbvcorrector.corrected_lc.estimate_cdpp()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is looking like a pretty good light curve. However, the CDPP increased a little as we optimized the over-fitting metric. Which correction to use may depend on your application. Since we are interested in the transiting planet, we will choose the corrected light curve with the lowest CDPP. Below we compare the light curve between just using the pixel-derived design matrix to also adding in the CBVs as a joint, regularized  fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:22:04.505934Z",
     "iopub.status.busy": "2023-11-03T14:22:04.505554Z",
     "iopub.status.idle": "2023-11-03T14:22:04.662518Z",
     "shell.execute_reply": "2023-11-03T14:22:04.662211Z"
    }
   },
   "outputs": [],
   "source": [
    "_, ax = plt.subplots(3, figsize=(10, 6))\n",
    "cbvcorrector.lc.plot(ax=ax[0], normalize=False, label='Uncorrected LC', c='k')\n",
    "corrected_lc_just_pixel_dm.plot(normalize=False, label='Pixel-Level Corrected; CDPP={0:.1f}'.format(corrected_lc_just_pixel_dm.estimate_cdpp()), ax=ax[1], c='m')\n",
    "corrected_lc_joint_fit.plot(normalize=False, label='Joint Fit Corrected; CDPP={0:.1f}'.format(corrected_lc_joint_fit.estimate_cdpp()), ax=ax[2], c='b')\n",
    "ax[0].set_title('Comparison Between original and final corrected lightcurve');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The superiority of the bottom curve is blatantly obvious. We can clearly see HAT-P 11b on it's 4.9 day period orbit. Our over-fitting metric settled at about 0.35 indicating we might still be over-fitting and should keep that in mind. However the low CDPP indicates the over-fitting is probably not over transit time-scales. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some final comments on CBVCorrector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Application to Kepler vs K2 vs TESS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CBVCorrector works equally across Kepler, K2 and TESS. However the Multi-Scale and Spike basis vectors are only available for TESS[<sup>1</sup>](#fn1). For K2, the [PLDCorrector](https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-3-k2-pldcorrector.html) and [SFFCorrector](https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-3-k2-sffcorrector.html) classes might work better than `CBVCorrector`.\n",
    "\n",
    "If you want to just get the CBVs but not generate a CBVCorrector object then use the functions _load_kepler_cbvs_ and _load_tess_cbvs_ within the cbvcorrector module as explained [here](https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-2-how-to-use-cbvs.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<sup>1</sup><span id=\"fn1\"> Unfortunately, the Multi-Scale and Spike CBVs are not archived at MAST for Kepler/K2.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applicability of the Over- and Under-fitting Goodness Metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The under-fitting metric computes the correlation between the corrected light curve and a selection of neighboring SPOC SAP light curves. If the light curve you are trying to correct was not generated by the SPOC pipeline (I.e. not a SAP light curve), then the neighboring SAP light curves might not contain the same instrumental systematics and the under-fitting metric might not properly measure when under-fitting is occuring. \n",
    "\n",
    "The over-fitting metric examines the periodogram of the light curve before and after the correction and is therefore indifferent to how the light curve was generated. It simply looks to see if noise was injected into the light curve. The over-fitting metric is therefore much more generally applicable.\n",
    "\n",
    "The Goodness Metrics are part of the `lightkurve.correctors.metrics` module and can be computed directly with calls to `overfit_metric_lombscargle` and `underfit_metric_neighbors`. A savvy expert user can use these and other quality metrics to generate their own Loss Function for optimizing a fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aligning versus Interpolating CBVs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, all loaded CBVS in `CBVCorrector` are \"aligned\" to the light curve cadence numbers (`CBVCorrector.lc.cadenceno`). This means only cadence numbers that exist in both the CBVs and the light curve will have values in the returned CBVs. All cadence numbers that exist in the light curve but not in the CBVs will have NaNs returned for the CBVs on those cadences and the Gap Indicator set to True. Any cadences in the CBVs not in the light curve will be removed from the CBVs.\n",
    "\n",
    "\n",
    "If the light curve cadences do not overlap well with the CBVs then you can set `interpolate_cbvs=True` when generating the `CBVCorrector` object. Doing so will generate interpolated CBV values for all cadences in the light curve. If the light curve has cadences past either end of the cadences in the CBVs then one must extrapolate. A second argument, `extrapolate_cbvs`, can be used to also extrapolate the CBV values to the light curve cadences. If `extrapolate_cbvs=False` then the exterior values are set to NaNs, which will probably result is a very poor fit.\n",
    "\n",
    "**Warning**: *The safest method is to align*. This will not generate any new values for the CBVs. Interpolation can be potentially dangerous. Interpolation uses Piecewise Cubic Hermite Interpolating Polynomial (PCHIP), which can be more stable than a simple spline, but no interpolation method works in all situations. Extrapolation is even more dangerious, which is why an extra parameter must be set if one desires to extrapolate. *Be sure to manually examine the extrapolated CBVs before use!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Joint Fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By including the `ext_dm=` parameter in the `correct_*` methods we allow for joint fitting between the CBVs and other design matrices. Generally speaking, if fitting a collection of different models to a system, joint fitting is ideal. For example, if performing transit analysis one could add in a transit model to the joint fit to get the best transit recovery. The fit coefficient to the transit model is stored in the `CBVCorrector` object after fitting and can be recovered."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hyperparameter Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any model fitting should include a hyperparameter optimization step. The `correct_optimizer` is essentially a 1-dimensional optimizer and is very fast. More advanced hypterparameter optimization can be performed by tuning the `alpha` and `l1_ratio` parameters in `correct_elasticnet` plus the number and type of CBVs, along with an external design matrix. The optimization Loss Function can use a combination of the `under_fitting_metric`, `over_fitting_metric` and `lc.estimate_cdpp` methods. Writing such an optimzer is left as an exercise to the reader and to be tuned to the reader's particular application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More Generalized Design Matrix Priors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main [CBVCorrector.correct*](https://docs.lightkurve.org/reference/api/lightkurve.correctors.CBVCorrector.html?highlight=cbvcorrector%20correc#lightkurve.correctors.CBVCorrector) methods utilize a similar prior for all design matrix vectors as is typically used in L1-Norm and L2-Norm regularization. However you can perform fine tuning to the correction using more sophisticated priors. After performing a fit with one of the `CBVCorrector.correct*` methods, `CBVCorrector.design_matrix_collection` will have the priors set. One can then manually adjust the priors and use `CBVCorrector.correct_regressioncorrector` to perform the standard [RegressionCorrector.correct](https://docs.lightkurve.org/reference/api/lightkurve.correctors.RegressionCorrector.correct.html?highlight=regressioncorrector%20correct#lightkurve.correctors.RegressionCorrector.correct) correction. An illustration is below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:22:04.664859Z",
     "iopub.status.busy": "2023-11-03T14:22:04.664746Z",
     "iopub.status.idle": "2023-11-03T14:22:04.666506Z",
     "shell.execute_reply": "2023-11-03T14:22:04.666244Z"
    }
   },
   "outputs": [],
   "source": [
    "# 1) Perform an initial optimization with a L2-Norm regularization\n",
    "# cbvCorrector.correct(cbv_type=cbv_type, cbv_indices=cbv_indices);\n",
    "\n",
    "# 2) Examine the quality of the resultant lightcurve in cbvcorrector.corrected_lc\n",
    "# Determine how to adjust the priors and make changes to the design matrix\n",
    "# cbvCorrector.design_matrix_collection[i].prior_sigma[j] = # ... adjust the priors\n",
    "\n",
    "# 3) Call the superclass correct method with the adjusted design_matrix_collection\n",
    "# cbvCorrector.correct_regressioncorrector(cbvCorrector.design_matrix_collection, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `cbvCorrector.corrected_lc` will now be the result of the fit using whatever `cbvCorrector.design_matrix_collection` you had just provided."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NaNs, Units and Normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NaNs are removed from the light curve when used to generate the `CBVCorrector` object and is stored in `CBVCorrector.lc`. \n",
    "\n",
    "The CBVCorrector performs its corrections in absolute flux units (typically electrons per second). The returned corrected light curve `corrected_lc` is also in absolute units and the median flux of the light curve is preserved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:22:04.668198Z",
     "iopub.status.busy": "2023-11-03T14:22:04.668081Z",
     "iopub.status.idle": "2023-11-03T14:22:04.670116Z",
     "shell.execute_reply": "2023-11-03T14:22:04.669868Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LC unit: electron / s\n",
      "Corrected LC unit: electron / s\n"
     ]
    }
   ],
   "source": [
    "print('LC unit: {}'.format(cbvCorrector.lc.flux.unit))\n",
    "print('Corrected LC unit: {}'.format(cbvCorrector.corrected_lc.flux.unit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goodness metrics are computed using median normalized units in order to properly calibrate the metric to be between 0.0 and 1.0 for all light curves. The normalization is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-03T14:22:04.671587Z",
     "iopub.status.busy": "2023-11-03T14:22:04.671484Z",
     "iopub.status.idle": "2023-11-03T14:22:04.680876Z",
     "shell.execute_reply": "2023-11-03T14:22:04.680637Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normalized Light curve units:  (i.e astropy.units.dimensionless_unscaled)\n"
     ]
    }
   ],
   "source": [
    "normalized_lc = cbvCorrector.lc.normalize()\n",
    "normalized_lc -= 1.0\n",
    "print('Normalized Light curve units: {} (i.e astropy.units.dimensionless_unscaled)'.format(normalized_lc.flux.unit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}