{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Producing and loading survey masks\n",
    "\n",
    "It is never the case in optical and NIR imaging that data exists without any abnormalities; there are always regions where the output flux is not necessarily reliable. The primary culprits of these are image edges and stellar diffraction spikes, however there may well be other image artefacts present. Over the first few years of JWST/NIRCam operations, which observational astrophysics have been using extensively, these artefacts consist of wisps, claws, and snowballs which are not always properly accounted for in data reduction pipelines. For more information on these types of artefacts, please visit the official JWST user documentation following the hyperlinks for [wisps, claws](https://jwst-docs.stsci.edu/depreciated-jdox-articles/nircam-claws-and-wisps#gsc.tab=0), and [snowballs](https://jwst-docs.stsci.edu/depreciated-jdox-articles/data-artifacts-and-features/snowballs-and-shower-artifacts#gsc.tab=0). In the segmentation process, these artefacts may be picked up they can often be bright, however we choose to also mask these specific areas just to make sure.\n",
    "\n",
    "Galfind utilizes two different masking methods which we will cover in the following examples. As usual, our first code block will be to instantiate our Data object with the JOF data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading GALFIND config file from: /nvme/scratch/work/austind/GALFIND/galfind/../configs/galfind_config.ini\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Aperture corrections for ACS_WFC loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/ACS_WFC_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for WFC3_IR loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/WFC3_IR_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for NIRCam loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/NIRCam_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for MIRI loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/MIRI_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for ACS_WFC loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/ACS_WFC_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for WFC3_IR loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/WFC3_IR_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for NIRCam loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/NIRCam_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for MIRI loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/MIRI_aper_corr.txt\n",
      "INFO:galfind:Aperture corrections for NIRCam loaded from /nvme/scratch/work/austind/GALFIND/galfind/Aperture_corrections/NIRCam_aper_corr.txt\n",
      "INFO:galfind:Loaded aper_diams=<Quantity [0.32, 0.5 , 1.  , 1.5 , 2.  ] arcsec> for F277W+F356W+F444W\n"
     ]
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "from copy import deepcopy\n",
    "import shutil\n",
    "import astropy.units as u\n",
    "from galfind import Stacked_Band_Data, Data, config\n",
    "from galfind.Data import morgan_version_to_dir\n",
    "from galfind import useful_funcs_austind as funcs\n",
    "\n",
    "survey = \"JOF\"\n",
    "version = \"v11\"\n",
    "instrument_names = [\"NIRCam\"]\n",
    "aper_diams = [0.32, 0.5, 1.0, 1.5, 2.0] * u.arcsec\n",
    "forced_phot_band = [\"F277W\", \"F356W\", \"F444W\"]\n",
    "\n",
    "JOF_data = Data.from_survey_version(\n",
    "    survey = survey,\n",
    "    version = version,\n",
    "    instrument_names = instrument_names, \n",
    "    version_to_dir_dict = morgan_version_to_dir,\n",
    "    aper_diams = aper_diams,\n",
    "    forced_phot_band = forced_phot_band,\n",
    ")\n",
    "JOF_data_2 = deepcopy(JOF_data)\n",
    "JOF_data_3 = deepcopy(JOF_data)\n",
    "JOF_data_4 = deepcopy(JOF_data)\n",
    "JOF_data_5 = deepcopy(JOF_data)\n",
    "JOF_data_6 = deepcopy(JOF_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: Manually masking the data\n",
    "\n",
    "In this example, we will convert our pre-created .reg mask paths made in ds9 (or other image analysis package which produces an identical .reg output) to pixel masks using the `Data.mask()` method. We explicitly passing in `manual` as our preferred method. This requires masks to the relevant filters to be included in the following sub-directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/reg\n"
     ]
    }
   ],
   "source": [
    "print(f\"{config['Masking']['MASK_DIR']}/{JOF_data.survey}/reg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The internal code uses glob.glob to extract the required paths which contain the band name labelled in 1 of 4 ways. For example, the F444W filter could be written as either:\n",
    "1. F444W (obviously)\n",
    "2. f444w (lower case)\n",
    "3. F444w (capitalized)\n",
    "4. f444W (upper case with lowered F)\n",
    "It is worth noting that including more than 1 mask path per filter in this directory will cause the code to crash.\n",
    "\n",
    "The first thing the code will do post path extraction will be to 'clean' these mask paths for regions with zero size. This may sound a little odd but for certain versions of ds9 (don't ask me which as I don't know), it is possible to accidentally create these regions, which then cause the entire region to be masked. Obviously this is not the expected behaviour, and so to get rid of these explicit cleaning of these .reg files is required. Following cleaning, these .reg masks will be converted to a pixel mask on the same scale as the SCI/ERR/WHT imaging and saved as a .fits mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fcef92b02b0> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_manual.fits\n"
     ]
    }
   ],
   "source": [
    "JOF_data.mask(method = \"manual\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look to see if we can find the newly created fits maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'F090W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F090W_v11_manual.fits', 'F115W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F115W_v11_manual.fits', 'F150W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F150W_v11_manual.fits', 'F162M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F162M_v11_manual.fits', 'F182M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F182M_v11_manual.fits', 'F200W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F200W_v11_manual.fits', 'F210M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F210M_v11_manual.fits', 'F250M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F250M_v11_manual.fits', 'F277W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F277W_v11_manual.fits', 'F300M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F300M_v11_manual.fits', 'F335M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F335M_v11_manual.fits', 'F356W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F356W_v11_manual.fits', 'F410M': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F410M_v11_manual.fits', 'F444W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/manual/F444W_v11_manual.fits', 'F277W+F356W+F444W': '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_manual.fits'}\n",
      "F090W mask exists\n",
      "F115W mask exists\n",
      "F150W mask exists\n",
      "F162M mask exists\n",
      "F182M mask exists\n",
      "F200W mask exists\n",
      "F210M mask exists\n",
      "F250M mask exists\n",
      "F277W mask exists\n",
      "F300M mask exists\n",
      "F335M mask exists\n",
      "F356W mask exists\n",
      "F410M mask exists\n",
      "F444W mask exists\n",
      "F277W+F356W+F444W mask exists\n"
     ]
    }
   ],
   "source": [
    "print(JOF_data.mask_paths)\n",
    "\n",
    "for filt_name, path in JOF_data.mask_paths.items():\n",
    "    if Path(path).is_file():\n",
    "        print(f\"{filt_name} mask exists\")\n",
    "    else:\n",
    "        print(f\"{filt_name} mask does not exist\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And also see how this impacts the print statement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****************************************\n",
      "DATA OBJECT:\n",
      "----------\n",
      "SURVEY: JOF\n",
      "VERSION: v11\n",
      "****************************************\n",
      "MULTIPLE_FILTER\n",
      "----------\n",
      "FACILITY: JWST\n",
      "INSTRUMENT: NIRCam\n",
      "FILTERS: ['F090W', 'F115W', 'F150W', 'F162M', 'F182M', 'F200W', 'F210M', 'F250M', 'F277W', 'F300M', 'F335M', 'F356W', 'F410M', 'F444W']\n",
      "****************************************\n",
      "****************************************\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(JOF_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with other examples we've seen, for example with the [segmentation maps](source_extraction.ipynb), attempting to load masks into an object which already contains them will not work. The original mask paths will remain. This is to prevent confusion with which masks have been used to create any products that require these mask paths, for example the band [depths](running_depths.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fcef92b02b0> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_manual.fits\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****************************************\n",
      "DATA OBJECT:\n",
      "----------\n",
      "SURVEY: JOF\n",
      "VERSION: v11\n",
      "****************************************\n",
      "MULTIPLE_FILTER\n",
      "----------\n",
      "FACILITY: JWST\n",
      "INSTRUMENT: NIRCam\n",
      "FILTERS: ['F090W', 'F115W', 'F150W', 'F162M', 'F182M', 'F200W', 'F210M', 'F250M', 'F277W', 'F300M', 'F335M', 'F356W', 'F410M', 'F444W']\n",
      "****************************************\n",
      "****************************************\n",
      "\n"
     ]
    }
   ],
   "source": [
    "JOF_data.mask(method = \"auto\")\n",
    "print(JOF_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Plotting the mask for a Stacked_Band_Data object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To be clear, the `Data.mask()` function we have used in example 1 loops through the stored `Band_Data` objects and individually masks each one of them. But what if we want to make a mask for a `Stacked_Band_Data` object? To have a look at this we'll instantiate a couple of fresh new JOF `Data` objects and produce a `Stacked_Band_Data` object for the F277W, F356W, and F444W filters commonly used for selection when performing forced photometry (as we will do in the notebook on [cataloguing the data](cataloguing_the_data.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Loaded aper_diams=<Quantity [0.32, 0.5 , 1.  , 1.5 , 2.  ] arcsec> for F277W+F356W+F444W\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<galfind.Data.Stacked_Band_Data object at 0x7fcfc02dbac0>\n"
     ]
    }
   ],
   "source": [
    "LW_nircam_stack = Stacked_Band_Data.from_band_data_arr(\n",
    "    [band_data for band_data in JOF_data_2[\"F277W+F356W+F444W\"]])\n",
    "print(LW_nircam_stack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fcfc02dbac0> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_manual.fits\n"
     ]
    }
   ],
   "source": [
    "LW_nircam_stack.mask(method = \"manual\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To conclude this example, we will plot the mask for the stacked data. We explicitly write the input arguments used by default here for clarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LW_nircam_stack.plot(ext = \"mask\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we have implemented the default mask plotting which is used for both `Band_Data` and `Stacked_Band_Data` objects. Let's have a peek at some of the plotting options which we can implement here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Automated masking of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fcef92d8400> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_auto.fits\n"
     ]
    }
   ],
   "source": [
    "JOF_data_3.mask(method = \"auto\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, let us have a look at the masks paths that have been produced here and plot the mask for the F444W filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****************************************\n",
      "DATA OBJECT:\n",
      "----------\n",
      "SURVEY: JOF\n",
      "VERSION: v11\n",
      "****************************************\n",
      "MULTIPLE_FILTER\n",
      "----------\n",
      "FACILITY: JWST\n",
      "INSTRUMENT: NIRCam\n",
      "FILTERS: ['F090W', 'F115W', 'F150W', 'F162M', 'F182M', 'F200W', 'F210M', 'F250M', 'F277W', 'F300M', 'F335M', 'F356W', 'F410M', 'F444W']\n",
      "****************************************\n",
      "****************************************\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAEtCAYAAAD+0e+cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQ80lEQVR4nO3deXwM9/8H8Nfskc25iSOHI4hq3bcipVQbglS1oo4q2mrRhjr6Q5VS2qL00lJ6ol/U0aoiRVOKqrhScR9thSg5XMnm3Ozx+f2RWlYSdlc2s5u8no9HHrUzn5l9z0x39r2f+RySEEKAiIiIiO5KIXcARERERO6CiRMRERGRjZg4EREREdmIiRMRERGRjZg4EREREdmIiRMRERGRjZg4EREREdmIiRMRERGRjZg4EREREdmIiRMRERGRjZg4EZHbWLp0KSRJgiRJ2L17d5H1QgiEhoZCkiQ8/vjjRdZnZGTA09MTkiTh5MmTJb7Pxo0b0blzZwQFBcHb2xt169ZFv379sGXLFkuZc+fOQZIkvP/++0ViGDFiBCRJwltvveX4wRKRS2LiRERux9PTEytXriyyfOfOnfj333+h0WiK3W7t2rWQJAkhISFYsWJFsWXef/99PPHEE5AkCZMnT8ZHH32E6Oho/PXXX1i1atUd4xJC4JVXXsEXX3yBN998k4kTUTmkkjsAIiJ79ezZE2vXrsUnn3wClermbWzlypVo3bo1rly5Uux2y5cvR8+ePVG7dm2sXLkS77zzjtV6o9GIt99+G127dsUvv/xSZPv09PQ7xjV69GgsXrwYU6ZMwcyZMx04MiJydaxxIiK3M3DgQFy9ehVxcXGWZQUFBfj+++/xzDPPFLtNcnIyfv/9dwwYMAADBgxAUlIS9uzZY1XmypUr0Ol06NChQ7H7CAoKKjGmMWPGYOHChZg8eXKRhIyIyg8mTkTkdurUqYPw8HB89913lmWbN29GZmYmBgwYUOw23333HXx8fPD444+jbdu2uO+++4o8rgsKCoKXlxc2btyIa9eu2RzPuHHj8Mknn2DSpEmYNWuWYwdFRG6BiRMRuaVnnnkG69evR15eHgBgxYoV6Ny5M6pXr15s+RUrVqB3797w8vICAPTv3x9r1qyB0Wi0lFEoFJgwYQISEhJQq1Yt9OzZE7NmzcKff/5ZYhwLFizAxx9/jAkTJmDOnDmleIRE5IqYOBGRW+rXrx/y8vKwadMmZGVlYdOmTSU+pjty5AiOHj2KgQMHWpYNHDgQV65cwdatW63KzpgxAytXrkTLli2xdetWTJkyBa1bt0arVq2K7YmXlpYGAHjggQdK8eiIyFUxcSIitxQYGIiIiAisXLkS69atg8lkQt++fYstu3z5cvj4+KBu3br4+++/8ffff8PT0xN16tQptnfdwIED8fvvv+P69ev45Zdf8Mwzz+DQoUPo1asX8vPzrcpOmjQJDz74IEaMGIHvv//eKcdKRK6DveqIyG0988wzeOmll5CamooePXogICCgSBkhBL777jvk5OSgUaNGRdanp6cjOzsbvr6+RdZptVp07doVXbt2hVqtxrJly7Bv3z507tzZUsbX1xebN29Gp06dMGjQIGi1WnTr1q1Uj5OIXAdrnIjIbT311FNQKBTYu3dviY/pboztNHPmTKxdu9bq74svvkBubi7Wr19/1/dq06YNACAlJaXIuipVquCXX35BtWrV0KdPH8THx9/TcRGR62KNExG5LV9fXyxatAjnzp1Dr169ii1z4zHdhAkT4OnpWWT9vHnzsGLFCjz77LPIzc3F4cOHER4eXqTc5s2bAQD169cv9n1q1KiBuLg4dOzYEVFRUdi5cyeaNm16D0dHRK6IiRMRubWhQ4eWuE6v1+OHH35A165di02aAOCJJ57A/PnzkZ6eDoVCgYceegjt27dH9+7dERoaioyMDKxfvx6///47nnzySbRs2bLE97v//vuxdetWPPLII4iMjMTu3btRt27dez5GInIdfFRHROVWbGwsMjIySqyNAoBevXrBaDRi1apVCAgIwJdffomQkBAsWbIEr7zyCt58801kZ2dj3rx5WL169V3fs0WLFti0aRMyMjIQERGBS5culeYhEZHMJCGEkDsIIiIiInfAGiciIiIiGzFxIiIiIrIREyciIiIiGzFxIiIiIrIREyciIiIiGzFxIiIiIrIRB8C0gdlsxqVLl+Dn5wdJkuQOh4iIiGwghEBWVhaqV68OhaJ06oqYONng0qVLCA0NlTsMIiIicsCFCxdQs2bNUtkXEycb+Pn5AQCO/RMKPz8+3SSisqMwqHCyfzSEYG23q1IozAgJu4RL/xT9Yvb2yUVujrcMUVVcVUKuInBRHAAgK8uMJvddsHyPlwYmTja48XjOz08BrZaJExGVnWsjHkWfAwFyh0F3s79SCSv8yzQMAvrom2H+im1Wy0qzmQ2zACIiF/bc8rZyh0BEt2CN0z3Q7vSG4e9AAICkMkHhZbi5UhJQ+BZYb+BlAFRm62XeRggJyKsPKHMFPC5IgEZAKItOIWguobbXrJEgSkiBhVKySo+FBBi0BigMakgmQFEACDv+LzCrbS9b/PaF711srIrbjlkyQyisz5fCoIZQmopdfq8UBuvXRm8DlPlqSP+9lUkjoNRLxZa9ocDfAI9MNSRD8VNA6qsa4ZWihGQsfH3jfN7YnyK38L9CU3itFPmAZJKA/95XXxfQnANgBqScWy6ctxEwKIBrXlbvZ85RA+ab/wMIswRzjsaqjCHTehu97ub/aGaDEoZ8D8vrvCxvQNzcV272zbK5Wd6Wx0n6fA8YjYXxmUxK5OZ44tE/F0AojcWeFyrZgsdPoMOWQLnDIKL/MHFy1KymaDurB1IVOQAABSSobstePKC0fi0UUN1Wyef13yXobNLishA4ocqCp1BCgaLVin7m4i+Xl1CWeCE9oIDHLa/VkPBk+BkcPVEHV7PVuGaS4GtHDaa/h8n2wreRJKBLp6PYs7cx8vVFMz1PjRkatQkeahNatT6JZpM3QtemAJ6pSmBDKHZ80wOHT9bA2I++QN4LyYXbpCqhD5Sw5+EXcfhkDYdjM0NCVoF1TOGNL2H/iWow/ZcMhFbKx4XrnhAAdCXMjf1u/huYEDwDORAQsC4jQcJ7Yzdi9vze+AuFSbWXKPx/JE8ywSwBOf9lVN7/Lc+VTDDCjIL/srcWRj8cVmYVJtu4mYR4QQUjzMiUrJN1I8ww3xKHGQJGyTrpLID1a4NkfY2NuNd5wE0AcpCaI8GgvcddVUAnjjwA4LrcYRDRfyQhSvgGIAudTgd/f3+cT69taePUxm8M0m5UD1CpqGb2wVjfAAyZsApeXf+GIb4mdv+vK5YcqI2dHimWpCDzl29hCjQit47AJ9p3MGXzOzAH6fF683fwueasbPFrhBIXT36Jqo1euGOZ6mYfJCl1ZRiZa0hLXooC/xKq6qhErfxexWVFntxhELmNPvramJ8/FQCg05lRO+g8MjMzodWWzi831jg5KEPSyx1CuTLPtzr6jfwJwqTAr99G4rspz+APdVphsqQ5b1U2u5UZOW+0hYePHqFB2Xi4xxhsfXsd3v73TdRpMhWT8y7KcgzVzT5Inf3wHcvoJVOFTJoA4LaKLbJRnsTHm0T2yFA49zPDxMlBIcIH56UK+gVYyjRCiVNXPdHrvR44rbgOvZQLeJRcm2dWG/Dlol74HlfxboMsHLmaj0bTu2LT3oZ44eS7qN55FIaezy7DIyhU16zFxK8fBTwulPl7uwNVthkFJXU8IiIqJbmS401KbMFedQ46r2DSVFr0kglfa87iiPIK9Hf5Hz74vxbyy6V0JCl1GHam8JfFdUmPDlsCsbzpZDyx9hPEPlgAVTHtxJzpmqTHz0yaSsYaJ4eoeZsmcin8RJJb6WwIAgCk/Nco//bHGK9lpWBQi+lo9+x2/NqlbBvUHlJdLtP3o4phlCJE7hCI6BZMnMit9G+dDL8Td37C/LPHBTQe3QdmswIr7/cso8iInOO+0Gtyh0BEt2DiRG6lXa8/8Pfcrnctl6bIxaM7tfjhGH+tuwplBm83jhh1Pl/uEIjoFryTkVvxiziDmas62Fz+B8055wVDdpEMnGvNEfnsVUdkF/M9jz13Z0ycyK3kPSDwsyZZ7jCIiMhF3T4QcGlj4kRuRZNiLoWRrImIqLwyOrkLLxMnB/mKe58bjez369DhcodAjsrmZ8YR/KFA5FqYODnIxJuZLIackDsCcpTIZA9HR8zzrS53CER0C5dJnObMmQNJkjB27FjLsvz8fMTExKBKlSrw9fVFdHQ00tLSrLZLTk5GVFQUvL29ERQUhAkTJsBotG5MuWPHDrRq1QoajQb16tXD0qVL7zleToMgD553qmhe2jRL7hCI6BYukTgdOHAAn3/+OZo1a2a1fNy4cdi4cSPWrl2LnTt34tKlS+jTp49lvclkQlRUFAoKCrBnzx4sW7YMS5cuxbRp0yxlkpKSEBUVhS5duiAxMRFjx47Fiy++iK1bt5bZ8REROcpQg7XbRK5E9sQpOzsbgwYNwpdffolKlW5OZJWZmYmvv/4aH374IR599FG0bt0aS5YswZ49e7B3714AwC+//IITJ05g+fLlaNGiBXr06IG3334bCxcuREFBYav6xYsXIywsDB988AEaNmyIUaNGoW/fvvjoo49kOV6iikoUKOUOwS3tG/is3CEQuZX88j5XXUxMDKKiohAREWG1PCEhAQaDwWp5gwYNUKtWLcTHxwMA4uPj0bRpUwQHB1vKREZGQqfT4fjx45Yyt+87MjLSso/i6PV66HQ6qz8iujcF6Vq5Q3BL3/9RX+4QiNxKrmRw6v5lTZxWrVqFP//8E7Nnzy6yLjU1FR4eHggICLBaHhwcjNTUVEuZW5OmG+tvrLtTGZ1Oh7y8vGLjmj17Nvz9/S1/oaGhDh0fEd3CzAEwHTFqwC65QyCiW8iWOF24cAFjxozBihUr4OnpWr1tJk+ejMzMTMvfhQuc8Z6I5FGr3Wm5QyCiW8iWOCUkJCA9PR2tWrWCSqWCSqXCzp078cknn0ClUiE4OBgFBQXIyMiw2i4tLQ0hIYXzj4WEhBTpZXfj9d3KaLVaeHl5FRubRqOBVqu1+iOieyMEa5wcsf/7TnKHQES3kC1xeuyxx3D06FEkJiZa/tq0aYNBgwZZ/q1Wq7Ft2zbLNqdPn0ZycjLCw8MBAOHh4Th69CjS09MtZeLi4qDVatGoUSNLmVv3caPMjX0QUdnISQ2QOwS3VOt+1ngTuRKVXG/s5+eHJk2aWC3z8fFBlSpVLMuHDRuG8ePHo3LlytBqtRg9ejTCw8PRvn17AEC3bt3QqFEjDB48GHPnzkVqaiqmTp2KmJgYaDQaAMDIkSOxYMECTJw4ES+88AK2b9+ONWvWIDY2tmwPmKiCMxnYq84RC5dGAJqzcodB5DacPcmvbImTLT766CMoFApER0dDr9cjMjISn332mWW9UqnEpk2b8PLLLyM8PBw+Pj4YOnQoZs6caSkTFhaG2NhYjBs3DvPnz0fNmjXx1VdfITIyUo5DIiKyy6+qq3KHQORW8uHc4QgkIQRHV7sLnU4Hf39/nE+vDa228OlmK79XcVlRfK88IirqzPBj8JpX8jAgVLww/+Gcr47IDhqhxN+6xQAAnc6M2kHnkZmZWWrtlWUfx8ldjdYEyh0CkVvJuuIvdwhuqZ6p0t0LEVGZYeLkoAIDewgR2UOfp5E7BLcUaOZ5I3IlTJwcNAfJcodARBXA7+oUuUMgolswcXIQ2xwQERG5HpNkdur+mTgRUZnI0fnIHQIRVQDOrthg4uQAyczTRmQvXaaf3CEQEd0zZgBERC6KP9KIXA8/lURErkrwFk3kavipJKIykZ/PbvVE5P6YOBFRmUjnAJhEVA4wcXKAUDi3qyNReST4sSGicoCJkwMkk0vPjUxE5YRSzxkKiFwNEycH8GZGZD8h+LkhIvfHxMkBZrXcERC5n0oBOXKHQER0z5g4OUAoON0Kkb2S/q0kdwhERPeMiZMD+KiOyH66AqXcIbgdhUHuCIjodkyciKhMmFlRazdFHrsiErkaJk5EVCbSOB6B3Yy+vEUTuRp+Kh2gzOVPZyJ7FTh5xvLySMV7DZHLYeLkAIk/nIns5sPbjd0kExMnIlfDOxkRlQk2DbefIlfuCIjodkyciKhM/KsokDsEt6PIZQ9eIlfDxImIykS2ZJQ7BCKie8bEiYiIiMhGTJwcoMxmg00iIqKy1shUWe4QmDg5QpHPdgdE9sqXTHKH4Haky55yh0DkUn4c+Tu+CdXCV9x50ljJ7Lz0hokTEZWJLLBxuN2MvEUT3ar74nBEz/sap6bsQidDdVli4KeSiIiI3EKSUocOT4+EX8tkbDj4Lj7UhkCFsn0KxMSJiIiI3MYx5VU83CcG0kUvDPh7Dg70PYdQs1+ZvT8TJyIqE0aJnSrsJQy8RRMV54jyCh5+ahS8t2hR9etfcOSb5XhCX9uyXjI5b8hdfiodIOlUcodA5HayJbZxspcxTSt3CEQu64jyCh4c9AK8v62B7KfTseT0dCyoGggVJKdOjcbEyQH5Ddk7iIjKgGAPXqI7+UuZgWYx/aD8oCFya5vw1Kl52B993qnvycTJAerrnOWXiIjIFVxUZKPpWxHQT24Hs9qAwG+2wqQxOO39mDg5QPmPl9whELkdM9jGiYic46oiH00XtkLacz2cOoYTwMSJiMpIPviI217GbI3cIRC5jTzJiPY/huJ4z2ec+j5MnBxh4mkjIufLv1Z2XayJygMjBCLjveCReeeRxe8FMwAHGFN5MyMiInJVktF5TQOYODlAOPn5KVF5ZHJm/2AiojIiawawaNEiNGvWDFqtFlqtFuHh4di8ebNlfX5+PmJiYlClShX4+voiOjoaaWlpVvtITk5GVFQUvL29ERQUhAkTJsBoNFqV2bFjB1q1agWNRoN69eph6dKlZXF4RHQLIxuHE1EZUWaX0xqnmjVrYs6cOUhISMDBgwfx6KOPonfv3jh+/DgAYNy4cdi4cSPWrl2LnTt34tKlS+jTp49le5PJhKioKBQUFGDPnj1YtmwZli5dimnTplnKJCUlISoqCl26dEFiYiLGjh2LF198EVu3bi3z4yUiskf2VTYLIHKEwonj7UpCCJf6GVi5cmXMmzcPffv2RWBgIFauXIm+ffsCAE6dOoWGDRsiPj4e7du3x+bNm/H444/j0qVLCA4OBgAsXrwYkyZNwuXLl+Hh4YFJkyYhNjYWx44ds7zHgAEDkJGRgS1bthQbg16vh16vt7zW6XQIDQ3F+fTa0GoV8FoaiqpjujvxLBCVTxcyv5Q7BLdy+YVItPqhltxhELmd6we/Qfb9Juh0ZtQOOo/MzExotaUzEr/LNNYxmUxYtWoVcnJyEB4ejoSEBBgMBkRERFjKNGjQALVq1UJ8fDwAID4+Hk2bNrUkTQAQGRkJnU5nqbWKj4+32seNMjf2UZzZs2fD39/f8hcaGmq1vuCq7z0fLxHR3dw3cI/cIRDRbWRPnI4ePQpfX19oNBqMHDkSP/74Ixo1aoTU1FR4eHggICDAqnxwcDBSU1MBAKmpqVZJ0431N9bdqYxOp0NeXl6xMU2ePBmZmZmWvwsXLlitN+U7r5sjEdENpjQ+qiNyhGRw3nRFss9WW79+fSQmJiIzMxPff/89hg4dip07d8oak0ajgUZT8sBzPgOPA3PblWFERFQRJW1rJncIRG5JuqYGYLxrOUfInjh5eHigXr16AIDWrVvjwIEDmD9/Pvr374+CggJkZGRY1TqlpaUhJCQEABASEoL9+/db7e9Gr7tby9zeEy8tLQ1arRZeXo5NnZJbx6WahRFROaV24nxbROQY2R/V3c5sNkOv16N169ZQq9XYtm2bZd3p06eRnJyM8PBwAEB4eDiOHj2K9PR0S5m4uDhotVo0atTIUubWfdwoc2MfjvC45nKnjYjKoX//qSl3CER0G1lrnCZPnowePXqgVq1ayMrKwsqVK7Fjxw5s3boV/v7+GDZsGMaPH4/KlStDq9Vi9OjRCA8PR/v27QEA3bp1Q6NGjTB48GDMnTsXqampmDp1KmJiYiyP2kaOHIkFCxZg4sSJeOGFF7B9+3asWbMGsbGxDsed/XGL0jh8IqI7Mpud106DqFwzldM2Tunp6RgyZAhSUlLg7++PZs2aYevWrejatSsA4KOPPoJCoUB0dDT0ej0iIyPx2WefWbZXKpXYtGkTXn75ZYSHh8PHxwdDhw7FzJkzLWXCwsIQGxuLcePGYf78+ahZsya++uorREZGOhx3Zkplxw+aiIiInMqc7gsg1yn7drlxnFyRTqeDv7+/ZRynK8O6oeX3teUOi8jtcBwn+xzoOBx9jvIWTWSva59vRM6A1PI9jhMREVm7fDlA7hCI6DZMnIiIXJQQbONE5BAnfnaYODkgL8exYQyIiIjI+fSppfNYrjhMnBygu+68C0JERET3xmxwXt83Jk4OYPU5EZWF6zrWbhO5GiZOREQuKjOH82ISOcKZFRxMnIiozEhm3nKIyPmMeR5O2zfvYg7I1vnIHQIRERGVwCtQ57R9M3FyQKbOV+4QiIiIqARXTtVw2r6ZOBFR2RG85dgj1cDzReQIt2njlJvrnHlhiIgqonxwuhUiR2SmV3Lavu1OnB577DFcvHixyPL9+/ejRYsWpRGTy+NwBESOkcz87BCR8xXondcj1e7EydPTE82aNcPq1asBAGazGW+99RY6duyInj17lnqArijtKts4ETlCMskdARHRvbF7aM3Y2FgsXLgQL7zwAn766SecO3cO58+fx6ZNm9CtWzdnxOhy8vVsd0DkCMksdwTuxcBHdUQOMTtx6BOHxiSPiYnBv//+i/feew8qlQo7duzAQw89VNqxuawaQdlAutxRELkfZZ6A0VvuKNzHZUWB3CEQuaUraVXwgJP2bXdKdv36dURHR2PRokX4/PPP0a9fP3Tr1g2fffaZM+JzSdm5GrlDIHJPrECxi+AJI3KIwah02r7trnFq0qQJwsLCcOjQIYSFheGll17C6tWr8corryA2NhaxsbHOiNOlpGdqAOddEyIiInKARijRxVAdHXv95LT3sLvGaeTIkdi1axfCwsIsy/r374/Dhw+joKBiVCunGOWOgMg9sY2TfQLNzps2gqg8qW3W4h1NKM5M2IsVZ2ZAOSPBae9ld43Tm2++WezymjVrIi4u7p4Dqui8hAp5EjMzKp9UmQIIkjsK9zG5XzzW/VRd7jCIXJJGKPGYoTrGdTmFtmNXI7tLHsxqA5w9oqRNidORI0fQpEkTKBQKHDly5I5lmzVrViqBVVRbHruMztudN3AXkaxY42SXA3uaAbgidxhELqWG2RcvKKtg2Ngf4fdiInJDTbgxM50yXw2fPzTIPVIDxnGnnfL+NiVOLVq0QGpqKoKCgtCiRQtIkgQhbjZavPFakiSYTOV/oBZnNtds/vx2YHu0E9+BiNzFlymSg32ficqfcGMIJrRMw6MTvkVudx1MmsLaJVWuGt67vXD6f52wet1DWCddxTBvfwwdN80pcdj0kUxKSkJgYKDl3xXdVSc9SnvYUA3ZJ1gtT+UYO4nZJZeP7amC8xIqRBeE4v8G/4Y645Ygq5ERWQA8MtXw3loZJ7/riFU/hWO98grOK3SAqjBHydVXdlpMNiVOtWvXLvbfFVW+k1q4vhN5Aid38VEnlV+Ka2oA5b9WmojuTaDZC6+oAjF88iqoh52GvooRhlQlvL+tgUNrOuG7X5thizodFxXZgPpske3TnThQtUOVwKdPn8ann36KkydPAgAaNmyI0aNHo379+qUanLsLNntDJxXY3Ni7+bSf8O7DMwAla/WonDJw1H17zG19Fb0Tec6o4gg0e2FaoC8GzVwK/ZOX4XVOIPurBtjzQ2d8d6Qaflen4aoiH9AUTZbKit2J0w8//IABAwagTZs2CA8PBwDs3bsXTZo0wapVqxAdzfY5N3hAif2DTuOVZR3xuzrlruWzm5uwFVllEBkRuYP2fXcBiY/IHQZRmehZEIpFb6yDT/XrOPNza/wyti025OXjqPIq9FIBoDkvd4gAHEicJk6ciMmTJ2PmzJlWy6dPn46JEydWiMQp38aZSi8osjB9yaPY+MP7+OKlMXgjJwXGEhp51DD7QijMOK24XpqhErkWIckdgVvp9Xp/QJUmdxhEZSJdocfTcx7HX8oMXJf0AJId7hxxzYlNAuyuA05JScGQIUOKLH/22WeRknL3WpXyIEOyfaDPdZrzGBY9AcNjZ2HPE5cQaPYqttyLqsKGbHpOH0/lmPkKJ6qzR5JCd/dCROXEQVU69qvS/kua7k22E0fbtTtxeuSRR/D7778XWb579248/PDDpRJUefOD5hz6PfQm6sf8ihMfbkAHQzXLug0tCxOlF8aug0em2qb9eQkVGpmc12OAyFmEnn3r7WFmN0Qil2P3XeyJJ57ApEmTkJCQgPbt2wMobOO0du1azJgxAxs2bLAqWx450lTzF48L6PrY/+GXlYsQGz8bb4XPwMfqJHT8ehk6NZ0K3xeOQHmHgS9rmH3R1RCE6PAzeKD1cSTuaon+zhnbi+ieVDF7opbZD4dUl+UOxe0FC29cRb7cYRC5HWdOkC2JW0eytIFCYVvaUJ4Gw9TpdPD398f59NrQahUY7fkO1jvYSO0+kz/+mLUB0tAk7Og6HG0iDiCgwUXkDrmIf558Go/8FmApW9usxVOmKujfey8aDvgDkpcBhxd3w4wtTbFN/W8pHR1R6dIIJY6/koiN33XFhIw0q3Z91xbGIufZSzJG516e8HmDCSiRA1oaA7EhZxZ0OjNqB51HZmYmtFptqezb7hons5lzJsx87jes/66uQ9v+o8xE4ymR2P3PYXSI/xyey6ojd8hFAMDSLa1QQ52Op4xBGNznDzQYuhqGpjkQm0LxwxvP46PzKhxTXgWYNJEL00smdFn4EA7NX4Zuh8LwxJIH8Y8yEwBQcM1X5ujcyyklO4sQOcKZj7nZ4MABu39ti3uZP+qqIh8tljTCxjOhaLxhjWX52MHb8fHjfyK/iw5eZyScXfQIFjzRBes8LuK6lAkoSyF4ojJwQZGFlmP64s8P1iFh1RGM7DsBqzRJMGRr5A7Nrdxn9i/8sURELoMjqzlgXqrtvepKYpLMeGdnPfhuu9nLqMqXccj9Owg/tn4NnTq/hpbf18bXmrOl0sOAqKxdUGShyWu9kX28BhYkvYnPQypBUrKxsz0izKXzaIGoonHmdEWscXJAtmRweNvaZi2GKCth6MsbUeWFr5EfZJ27fvPeQEzVX+CVoXIhTZGLpu8+gvi0ADx+7CN4nzMjW+6g3MhXKj6WJ3KEMxMn1jg5wNYpVG5QQUL3glD88lAujqz5GqOSZiCw8ymcmNoLU4JnQjIXXgbJrMCwc29j1QMaqMCBAql8uC7p0fKrpkh6phdy7mONkz3u5UcaETmHzYnT9u3by00vuXtlsrHRWRWzJyaYwnA6JhHfHZqBduM24u817fFezRlo9vSL6LAlEF9pzkIozNBu8YcyV4mcma3x8J7PsbN7OrwEq52ofMiTjHjkl6pQfvyA3KG4lftM/nKHQOSWnNk43ObE6cUXX0RgYCCeeeYZrF69GjodR7QtSRWzJz4KCMGxGXGYvHgBUv+ugbdavYNm0S+h9bpQfKhKwvn/RgQONfsBAJ5/cgJMXgItPmuFxB7PoO7XsTgyMrHEkcaJ3I0RAhnng+QOw61ksH0jkUNynVhba3PidPbsWezYsQONGjXCBx98gODgYHTt2hWffvopkpOTHXrz2bNn48EHH4Sfnx+CgoLw5JNP4vRp61Ed8/PzERMTgypVqsDX1xfR0dFIS7Oeuyk5ORlRUVHw9vZGUFAQJkyYAKPR+nHajh070KpVK2g0GtSrVw9Lly51KGZbvKwIhrenHjFvPIuGLz+NDlurYr76rCVZulVPQyAAYIMmGSnPR0INBaIOeOCzmtNRddh+HJu1GWEmNhCl8qFKowtyh+BWrio4+CWRI2x9MuQIu9o4NWvWDFOnTsX+/fvxzz//IDo6Gps3b0b9+vXRokULTJs2DQcPHrR5fzt37kRMTAz27t2LuLg4GAwGdOvWDTk5OZYy48aNw8aNG7F27Vrs3LkTly5dQp8+fSzrTSYToqKiUFBQgD179mDZsmVYunQppk2bZimTlJSEqKgodOnSBYmJiRg7dixefPFFbN261Z7DtzDcZT65d3AeI1KvY53mPNIUuXcsO6hnApT5ahgh0GV9bWT9Nw/e2+I8+racBnWlHPy55H9oYqriUKxEriTrfKDcIRAR3RO7Rw4vTk5ODrZs2YKffvoJP//8M8aPH4833njD7v1cvnwZQUFB2LlzJzp16oTMzEwEBgZi5cqV6Nu3LwDg1KlTaNiwIeLj49G+fXts3rwZjz/+OC5duoTg4GAAwOLFizFp0iRcvnwZHh4emDRpEmJjY3Hs2DHLew0YMAAZGRnYsmXLXeO6feTwUP+X7D62klz7fCMUGiMCnnuq2PXVzD7YGH0a9frsQ8TAl3FQlV5q701U1rZ1ysIDG1fJHYbbKM17DVFFohFK/K1b7JSRw0ulV52Pjw+io6Px7bffIi0tDS+95NiHPTOzcHThypULJ7BNSEiAwWBARESEpUyDBg1Qq1YtxMfHAwDi4+PRtGlTS9IEAJGRkdDpdDh+/LilzK37uFHmxj5up9frodPprP6cxfToFcS993SJ61MUOWj/Yyg+H/Eqfn7/OzTg5L7kxoRgb1Eicj6T5LxZTkp9OAKlUonAQPur481mM8aOHYsOHTqgSZMmAIDU1FR4eHggICDAqmxwcDBSU1MtZW5Nmm6sv7HuTmV0Oh3y8vKKxDJ79mz4+/tb/kJDQ+0+HlvlBxkx62jJk/sChY1qJ+VeRM//G4iYIA4fTu6ntlmLt5S10OaFOLlDISK6Jy7T3z0mJgbHjh3D7t275Q4FkydPxvjx4y2vdTqdU5MnWyfxPKhKx0HOvkBuQgUJXQpq4LVOf6P9q2uR+1gO9OwkRkRuzubE6dKlS6hevbpTghg1ahQ2bdqEXbt2oWbNmpblISEhKCgoQEZGhlWtU1paGkJCQixl9u/fb7W/G73ubi1ze0+8tLQ0aLVaeHkV7e6v0Wig0Th/Ti1/4QGfs6xBovKlktBgqLE6Rr26AZWHfw0AyF1XH3GTIxFS/TIabVkuc4REVN6phfO+W21+VNe4cWOsXLmyVN9cCIFRo0bhxx9/xPbt2xEWFma1vnXr1lCr1di2bZtl2enTp5GcnIzw8HAAQHh4OI4ePYr09JuNpuPi4qDVatGoUSNLmVv3caPMjX3IpbZZi1NvPC5rDESlJdjsjQ+1ITj93s+Y8cvbUHoYsO7x8ejT5E3cN7Mzhp7PRuol9qoj9/VNKIeGcRf2zvBhD5sTp3fffRcjRozA008/jWvXrpXKm8fExGD58uVYuXIl/Pz8kJqaitTUVEu7I39/fwwbNgzjx4/Hb7/9hoSEBDz//PMIDw9H+/btAQDdunVDo0aNMHjwYBw+fBhbt27F1KlTERMTY6k1GjlyJM6ePYuJEyfi1KlT+Oyzz7BmzRqMGzfOobif09ctleMPNXnj1Y1NSmVfRHLyFWqsibiIB9uewCevv4BOkWNx//zWePFiBnaoLzr1JkZUVqKXLJA7BHIBNidOr7zyCo4cOYKrV6+iUaNG2Lhx4z2/+aJFi5CZmYlHHnkE1apVs/ytXr3aUuajjz7C448/jujoaHTq1AkhISFYt26dZb1SqcSmTZugVCoRHh6OZ599FkOGDMHMmTMtZcLCwhAbG4u4uDg0b94cH3zwAb766itERkY6FPelUhqR9Jgqg8MLULmQLRkQtS0ID/9aGe/gPA6rrsBYzAB0uXnOfwRO5AwtjYFAhofcYZALcGgcpwULFmDcuHFo2LAhVCrrZlJ//vlnqQXnKpw5jhNRRTJTXQvPX3lT7jDcBu81rmN7Zx0e6HAcNWfJ28SDbHch80unjONkd6+68+fPY926dahUqRJ69+5dJHGqCJ7W18FazTm5wyByO2bnDa1C5FStZv2AhMl95Q6DXIBdWc+XX36J1157DRERETh+/LhD4zWVBy39jVjLKaSIiCqMrCZGbNvVFMA5uUMhmdmcOHXv3h379+/HggULMGTIEGfG5PKUCudNHkhUnumNpT7mLpHT3Zix4VdjnguNfkhysfl/AZPJhCNHjliNs1RR/XFdBbCNK5HdUguYONlDI5TQ32VScXK+yXULJ18/oSydHuXk3my+i8XFxTFp+s8GzXm5QyCiCmCYobbcIbilXgW1SnV/T0wtHMMwu5R6VJN7488/IiIXtU+RLXcIbmna03uxpX0+As1FZ4ZwRN4T16DUq0tlX+T+mDgRUZkpKGZsJypZAsd5c0iX9bVRv90JnJy/Hj0L7m2eUY1QwuRpgPeJUgqO3B4TJyIqM+l81EFlIFsyoMXH7ZF1LhDLj8/AZ4FVoYLk0L46GwrnaL3+Y7PSDJHcGBMnIiozdQVHXrZHfVMluUNwW5lSAZp+3A5p8zug98n3caDvOVQz+9i9n0kRhVVNv6x+rLRDJCeTTM7pAsnEyU7OuhBEFYGOj+psJpkVuMg2TvckWzKgxZJGuDC0OwK//BXHF63BYwb7Ojm1GbcJALDiQum0lyL3x8TJTpLZsepeIgJSFAVyh+A2hMKM+00Bcofh9vSSCR1jq+Fwt8HI638Z3x+aiemKknvdNfpvzKYbsh/WAwAOsL0Z/YeJk53MarbRIHLUNYVe7hDcSobE81UajBCIOuCBLS3GIL+GAiMvzsSW9vnwFUV7yu1Z/bmlPZRGKC33/DzJeM9xqCAhzKRFVCkPl0DFU+qdU9HBxMlOCgO7pBI5ysxHdTZTZauRpNTJHUa5MuzfTCwOnQaVzozGW/+HY+P2FanVy+6hwweVggAAHYwhAACP647f9ysJDboXhOKrGgE4OeIo9ry9CX3qZjq8P5IfEyc7KfPkjoCIKgIFBwx3ihnmZLxW5234HVVBOSMBB1Z/aVUD5JegxqD42QCAV9v8CwBQb61i13sEm73xrD4Mm9vp8ddHG7Hy5Aw8/eEXOBbfBN2n9sGwf5k4lQWhdM5+mTjZSbBtOJHDWONEruBbTRKebDcVfrEByIrKwP+Ov4WZ6sLkKX5qP+RVM+F1cx10HrURAHBw2d171FUz++BFfV3seuw6Tn3zHT5ImY5Wq7+DKV+NZZ3fQIu+L+KpI8Bh1RWnHhvdpMx3zn6ZBtjJ6G2Ar1Bz6H0iB2RKbBxOrmGH+iLa9R+OXe+tQ07M33gpeSZa9ByMKbtrYT2AydtnIrehGQDw1W+NAM25IvsINHshylANz3X/Ey2G/Ij8iEwYtAZIR1VImdQF8799FOs8/kWmlAo4qfaDSqbMNQP2j0BxV0yc7OSRyaSJyFGscbKDWe4Ayr/Tyuto9HoP7D5zGL7zd6HV9qVYPr4DAED34M2G+dvUqZZ/VxIaRBRUw7DOp9Bu2AYURF5Dgb8B+To11D8G49f5ffDxaT8cVKUDmrNlfkzkfEyc7KTK4N2MiJxPpeO9pixcVeSjxZJG2HC6FppsXA3fj38vUuY5czUcNRoxvOUlPDJsI8y9L0Jf1YgcswJ+f6pwaVE3LFrRBT+oL+GqIg9QsTGsK1Bdl1jjRETujTVOtpN4qsqMXjIhMt4Ln9afgL7H5sHobf1UYULKdKgzzcgNNSEPgGeqEqoPG2D9oiexKEXCEeUV1i65Iid1sGDiZKdSGMqDqMLK4WNucmGjr17GkBMSdG2slxu0Bhh9VPDbUBn7F/XAZ781xK8el5AtXWXbJVfmpAGrmTjZSXmdnxIiovLEX3igZ0ENvBKVgOzmxVf1KefXQ+PpPQqnwdGcL+MIyRHSVQ0QVvqDyDJxIiJyQZKTulJTIRUktDEGI+Z+HSJf3gD0vQB9FWOJbfKVPnoUOOvZDzmH0TkjLnEcJ3sJzlVH5CgT2zjZTJHJ37XOUM3sgwmmMBx97hS2HJiNRxLnQz8iCRCA55d1sLf9y9AsqltkO/2IJBx/fyOamOwbDJNk5KRHdUyc7HXVU+4IiNyWQeIvdip7Kkh4zFATG1qacPLb5Rh3ZTp85+9Cfk0JvmuCcaL7s4ipOQP1XovC0yeNGPDaEKvtNVdV0FxVQT8iCX9sWICeBaEyHQnZw5jm55T98ieNnYSep4yIyB0Em73xPILw4qgNqDzia+TUMSEvWw3vX31x+OvH8O3m1tikTkWaItdqgMunaxQ+J/X5R4mc+0wwaCXM9JuJ9+KnQfdoNpYfm4EZTd/BfDV70rkyYXROm2RmAURErojPAxzWxhiEiU2uo+uElch7PAOKAkCxS4vkiY9g+Ya22KBO/6+Rd/GJz+PP/QwAuPxZewQ++A9yBqTid1UGenR8HZsWfYOcIRfxxsVpqP/g/+GVy5xCxVUJEx/VuQRzrofcIRC5LSPbONlMX5ePNe2hgoQn9bVxoPcl/Lr3PXTZ+AWEUYkzA57E1KCZaDzwOXTYWhWLNGcLk6Y78BlwHADw9edRGPx8DABgRhMd4lWpaBLzNHSjO8Pgb0Lvk+/j57YFUIFtX11RblqAU/bLxMlOhiu+codQavhhJ3JdmvP8fNrqaX0dHH3uFJb9PAc1W/yDbSOH4ZVaM1D/pX545LcAfGVDsnSr3LDCBH+FdBmxHsnQ7vRG1zdWAwBSFDlotawh9nYcDmUe0Gzr/3CwXxL8BX9UuxpnPapj4lSBLQwJkDsEIiqBOO8vdwhuwUuo0LfZJaxY2g1de41F3be6IPqYGWs153BZYf/UJypIEIrCQQluJFsvRE5B/iNZljJ6yYS+x02YW20GfI8rUOXLOByduBvVzE6Y34McZuZwBK5BHZh190Juov+aj+QOgYjonuRJRvQ/XYB3cB77VWnIu8fpHRqYKgMAPNNvNgH+QXMO+3oNKVL2Q1USnmg3BV5LQyFNPYyj89dZtif5ZaZWcsp+mTjZSfnAVblDKBUqSDBV5iSiRK7KnK+WO4QKKVpdWGtk+PYBq+W9E4v/utylvoTGo/vgyrBuyB96Cfu+X4z2xhCnx0l3ZzLxUZ1LuPJ9S7lDKBXdCmpCedA52TgR3buCy84Zg4buLLr/DgDA8g+ftnmbNEUuHvy+DmIbjYehaTa27pjLsZ5cgMnAxMklGAtcewQHjVDipxZm1DDfuRH7tKcO4Nzm8pEEEpVL7IAoi9r99gEAPtfl2LWdEQIjUq/jlftmQqT4YEXiDPTT13FChGSrq+nOGeWdiZOddJcD5A7hjvSSCVv21cOh99ff8UPbYNIWbN3YoewCIyK7eD2QLncIFVJu+8JJYf9RZjq0/RrNObR85nn8+3FnfH50GpMnGZk55YprMOhdv93BIs1ZTB8zHF/+PhPfhGqthh248e/spib8qOM4MUQuy8sgdwQVktHbAO8L9/aIJ0mpQ8uV9bDk0Tfwxffvozsf28kiP885U6QxcbKTUukeycYizVnEPDwdfVbNx6Fn/rZ0k/3j8UuoYvaEUJhxWFnyiLf3mfzxnL4ufy0RyeT6bw3kDqHC8RKFTTGuf9bmnvdlhMB4XSq69Y/Bqx3/uef9kf2ys7ydsl8mTnZqMGyH3CHYbLkmCf0enIZaz+zFsYVr8bChGmov2ooxmiAAsOq2qxFKtDEGYZ5vdZwYega7396I3g+exWGVY9XVRHRv9NmcULysNTAVdpj5dP6TpbbP/ao09NzPwTHlUPf+ZKfs17VbOrsgXXxduUOwyy8eF9Cp+xj89vkybNw3C/gXGLXvLZjSVVBBQnNjIIYESug9dDMq9TwG4+mq2PxxNOZ+8xCOKK8AyutyHwIRUZmIVBTWzH+ruihzJFQazpwKQ7gT9itrjdOuXbvQq1cvVK9eHZIkYf369VbrhRCYNm0aqlWrBi8vL0REROCvv/6yKnPt2jUMGjQIWq0WAQEBGDZsGLKzrYfWP3LkCB5++GF4enoiNDQUc+fOdTjmA5vbO7ytXA6rrqDpy/2R8X1zZDUyIuc+E1TXBE4MP4pte+Zi6KY5yMv0xoxHp6NBTDSe/SevMGkiIqpAekQcgsKgxnVJL3coVAqM5XHk8JycHDRv3hwLFy4sdv3cuXPxySefYPHixdi3bx98fHwQGRmJ/Px8S5lBgwbh+PHjiIuLw6ZNm7Br1y4MHz7csl6n06Fbt26oXbs2EhISMG/ePLz11lv44osvHIo5pGaaQ9vJLR9GzJ3bF2pdYeP27AYmVGmajE3Pv4we7Saj6ZKG+FCV5NAUBURU+swmtqQoaw27J8BnE6e6oTuT9VFdjx490KNHj2LXCSHw8ccfY+rUqejduzcA4Ntvv0VwcDDWr1+PAQMG4OTJk9iyZQsOHDiANm0KG/N9+umn6NmzJ95//31Ur14dK1asQEFBAb755ht4eHigcePGSExMxIcffmiVYN1Kr9dDr7/5i0On01n+/eeBxgAul9IZcC4VJLQzhODVJlcRMWodRNdUqPf4wtC9sN2SPt0fU/5S47wqVeZIieh211M5dUdZ83joAlb3HQ+ATRSoZC77kyYpKQmpqamIiIiwLPP390e7du0QHx8PAIiPj0dAQIAlaQKAiIgIKBQK7Nu3z1KmU6dO8PC42TgvMjISp0+fxvXrxX84Zs+eDX9/f8tfaOjNrqQFBc4ZibQ0VRIavFJQF0efO4XYP2aj++xlOPdbU7x93zt4qs84SznvnmeQ+PVytDUGyxgtERVHmF329lxu5YVJeOfSvc11R+Wfy34yU1MLa0GCg62/1IODgy3rUlNTERQUZLVepVKhcuXKVmWK28et73G7yZMnIzMz0/J34cKFez+gMlDN7IOPAkJw4q1f8e63HyH7uh/efXgGWvUeidbrQjFffRZdtYBkUsH3LyXebP02hEmBuO3z8KS+ttzhExHJSmECLijKz0Tu5BzsVVcMjUYDjUYjdxh2UUHC4lZZMJtz8MYbQ/Gr+nLhDUCZZFWuz6A4+B1SYM6j0/FAtSw0HT4Acc//iW+OT0f1Ju/gM4+zMh0BVRSSSQWh5K96cj0XxkTcvRBVeC5b4xQSUji7dFqadWPstLQ0y7qQkBCkp1tPS2A0GnHt2jWrMsXt49b3KA+MEIhOVOKpI8ASzdkSfzWF9DuI0/O7YbZ0Dt+kARcV2Wi5tCE294vBu4lT8abEmiciV5DjpMH7qGRd190ndwjkBlw2cQoLC0NISAi2bdtmWabT6bBv3z6EhxeOzBAeHo6MjAwkJCRYymzfvh1msxnt2rWzlNm1axcMhpvTF8TFxaF+/fqoVKmS3XG1aHnG0UNyOqMNs4LmNDdj3trC83dIVdjIXS+ZMOjvPIxv/g7GfzUfc32qOzVOIrq7Wg3Pyx1ChcNhCMgWsiZO2dnZSExMRGJiIoDCBuGJiYlITk6GJEkYO3Ys3nnnHWzYsAFHjx7FkCFDUL16dTz55JMAgIYNG6J79+546aWXsH//fvzxxx8YNWoUBgwYgOrVC7/8n3nmGXh4eGDYsGE4fvw4Vq9ejfnz52P8+PEOxdxs0felceiyMXkasMnj32LXLdGcRdvnnkOXiIN4Tu9eA30SlTfGArakIHJFsn4yDx48iC5dulhe30hmhg4diqVLl2LixInIycnB8OHDkZGRgY4dO2LLli3w9Lw5FcGKFSswatQoPPbYY1AoFIiOjsYnn3xiWe/v749ffvkFMTExaN26NapWrYpp06aVOBTB3RyL6ePg0cpPI5RQ5qutplq53V/KDIT/VAOdFPzlRc6h1Esw8inUXZ082BAA24IRuRpZE6dHHnkEQpT8eEmSJMycORMzZ84ssUzlypWxcuXKO75Ps2bN8Pvvvzsc563+TGgAwD3HPapj1kL9vxp3LWeEwHY1pxwgJzHLHYB7EEKSOwQiKobLtnGi0tfE6I8FE4fJHQYR2SDtcoDcIRBRMZg4VSDVJAU+NrlnbRmVHxJrnGxSwDZORC6JiVMF8qMqHZlSgdxhUAWnzGXmRETui4lTBXJRkS13CESQTHJHQETkOCZOREQuKFfPR3VEroiJExGRC0rLVssdAhEVg4kTEZUp1XV2syci98XEiYjKlomJExG5LyZORERERDZi4kREZYq96myTLjhsA5ErYuJERGXripfcEbiFPI4USuSSmDgRUdky8bZDRO6LdzAiIiIiGzFxIqIyJdirziYGCLlDIKJiMHEiojJlTNPKHYJbuKLQyx0CERWDiRMRlSmzkbcdInJfvIMRERER2YiJExGVLfaqIyI3xjsYEZWp3LQAuUNwC9mSUe4QiKgYTJyIqEyxjZNt8pg4Ebkk3sGIiIiIbMTEiYjKlNmklDsEIiKHMXEiojKVkVpJ7hDcgpkDYBK5JCZORFSmWONkm2ypQO4QiKgYTJyIiIiIbMTEiYjKlIm96mzSyFhZ7hCIqBi8gxFRmbqSWkXuENxCZaF2aLv7TQGlGwgRWWHiRERlSghJ7hDcwoejf3Zou18nxZVyJER0KyZORFSmlCqT3CG4hfUrujm0nbbf0VKOhIhuxcSJiMpUrQeS5Q7BLUzNSrN7m0amypBOBZR+MERkwcSJiMrUn7tbyB2CWzBJZsu/a5h9bdpmaoNcXNrWxFkhERGYOBFRGcvK8ZI7BLczJdgTszxrQoU7tw/r+dZy7Pj5oTKKiqhiYuJEROSCjLeMHD4u/Toef2oXEgacRW2ztsRtciKz8cMlz7IIj6jCYuJERGWqoIAjh9tLL5nQbuUDMBsVOLxsGZ7U1y5SxkuoYFYbcFCVLkOERKUv0Fy0drqa2Qcacfd7SCWhQYs2J5wRFlRO2SsRUQlSrvNRnS18hRrZksHyOk8yov3auoi7rsU3p6bj0S6v49Vr6ahi9sRr3lURryv8HZx5h6laAs1eaGsMRGOlEnOVSQ7HphFK6CXn9I68/bjlikENBa5L+nvaz60xNjdWha9QI18y4bTy+l2P0V94oJrZF3VNPtimvgQ1FMiWDFBBgp/wQNYt11kpFAgRPuhuqIoN6jRck/LRxhiE/sEmfJsOnFJet5RVF1Nf4lFMIuJ5W3qgEIAGKnhAgWCTJy4oc+EllFDc9vjY+5Z9tYYn/hEmmG6pPVVCgpcojGFw+N/Q6zXY8qf1DwF/pUDzeunoNeU7fD99KFIu+xXGrjLh5UULkL7nfmxe+xgCtDnI0PnA3y/Xsq0kCfj65aLzsA0oeCYF150wcxETJyIiF3Ri4h8Y9m6fIsvf3lwNs/xyEX1mLlo/G4W3N7TCS4dn4NHxPaEwqNG1oGaRbRpJHuj3+H407BMPc+fLyPikDQ593POuMQQIFcwAlAAm9P8D+3a3QFKqL8ZPX4G4/0UiOaUSOnU8gk+3tIC4w6TEviU83FBAQuMq+Th6VQMAaF41H4PfWIElbw/GvxmaEvenlATGTV6DxF9bY/uBevBTWydQje5Pw8m/gyFuCcnf986JikIS8PXJBwD0Gr4RHpVy8MN7Ayzr1WoTvL1KTqR8tdlWr2vdfwGVw9KQ+Esb+PrloOWUb2EOMELKl5Ad2wB/720IL588SMWcGg9PPao1OwfPBy/AUNuA/GUNofbLR0ZSENQaA7T3pyLjdHVLeY1fLvwapCCvWyZm/1QVOeerwq/XSWQ3NmPoXg/kH7r5/4TKL6/I+yn9ih6XpM2/bQEAHwOgNsMYaIbqkgrC0wwora+72fvm6/zqCqgzzJBuKSIkQKgLky19FSMAoEO29WCvQgWY1QLZSiN6RH8IyVxYXigEspRGeD11BdHv7YNQmCGZFRAKM25nOUonJE6SEIJTcN+FTqeDv78/zqfXxtbG0zFelyp3SERua4IpDK9mvyF3GG7B+3wJjyQUQG6oyVImt7YJqlw1jN6GYrcpCFTA6G2dOJS471uYfG7WJuirGqEwqKHQA0ZfAySTCgqDBJOnAZorKkhFv7sszJqSG7UbfAXU2YXrTZ6ASWOAUq/GnSqThAIweRogmRVQFChxe4WJWW2AwnDbl7HChq+6/w7ixhexZL6Z1RT35UyuT6czo3bQeWRmZkKrLbl9oD0qVI3TwoULMW/ePKSmpqJ58+b49NNP0bZtW7v2oVKZoYIEI4Tlv45SQUILYyBSFLm4qMi++wYlqGH2RZoixyoWZ1ZjE92LPBNHDrdVbu27f4ZvlLmRGBW/TdFltuz7dma1Aeb/8hGhNML0X8Kir2q0e1+3KvC3fm3S3P1RHVCYzJg8i09ozGrb9nG3/RPdrsIkTqtXr8b48eOxePFitGvXDh9//DEiIyNx+vRpBAUF2byfIR9+jicO14I+2xNqzwIY9WoIISFP52NVTp+jgclk/TMo55YykkIgqGY6ao36BuJEFexc0AsAkHVLmQcfTUB+thc03vn4K/F+ZGb4WdZVqpKJwJAr+Od0GHrOWIILsS2RuKeZZf0jT/+G9L+rQ+VhxNWUKriYHFLkWCSFQNvHDmD/tgchzM77MpMUAiqlGQZD4fkIqX4ZqZcC7d6PUmWGEBLM/OJ1aw8+Git3CEREDqswj+ratWuHBx98EAsWLAAAmM1mhIaGYvTo0Xj99dfvuO2tj+q0WnZEJCIicgd8VOeggoICJCQkYPLkyZZlCoUCERERiI+PL1Jer9dDr7/ZWC4zMxMAkJXFalsiIiJ3ceN7uzTriCpE4nTlyhWYTCYEBwdbLQ8ODsapU6eKlJ89ezZmzJhRZHmT+y44LUYiIiJyjqysLPj7+9+9oA0qROJkr8mTJ2P8+PGW1xkZGahduzaSk5NL7cSTfXQ6HUJDQ3HhwoVSq24l+/AayI/XQH68BvKz5xoIIZCVlYXq1avfsZw9KkTiVLVqVSiVSqSlWc82npaWhpCQoo2mNRoNNJqiY4j4+/vzgyIzrVbLayAzXgP58RrIj9dAfrZeg9Ku8KgQLZ09PDzQunVrbNu2zbLMbDZj27ZtCA8PlzEyIiIicicVosYJAMaPH4+hQ4eiTZs2aNu2LT7++GPk5OTg+eeflzs0IiIichMVJnHq378/Ll++jGnTpiE1NRUtWrTAli1bijQYL45Go8H06dOLfXxHZYPXQH68BvLjNZAfr4H85L4GFWYcJyIiIqJ7VSHaOBERERGVBiZORERERDZi4kRERERkIyZORERERDZi4mSDhQsXok6dOvD09ES7du2wf/9+uUNyS7Nnz8aDDz4IPz8/BAUF4cknn8Tp06etyuTn5yMmJgZVqlSBr68voqOjiwxcmpycjKioKHh7eyMoKAgTJkyA0Wi0KrNjxw60atUKGo0G9erVw9KlS519eG5nzpw5kCQJY8eOtSzj+Xe+ixcv4tlnn0WVKlXg5eWFpk2b4uDBg5b1QghMmzYN1apVg5eXFyIiIvDXX39Z7ePatWsYNGgQtFotAgICMGzYMGRnZ1uVOXLkCB5++GF4enoiNDQUc+fOLZPjc3UmkwlvvvkmwsLC4OXlhfvuuw9vv/221VxmvAala9euXejVqxeqV68OSZKwfv16q/Vleb7Xrl2LBg0awNPTE02bNsXPP/9s/wEJuqNVq1YJDw8P8c0334jjx4+Ll156SQQEBIi0tDS5Q3M7kZGRYsmSJeLYsWMiMTFR9OzZU9SqVUtkZ2dbyowcOVKEhoaKbdu2iYMHD4r27duLhx56yLLeaDSKJk2aiIiICHHo0CHx888/i6pVq4rJkydbypw9e1Z4e3uL8ePHixMnTohPP/1UKJVKsWXLljI9Xle2f/9+UadOHdGsWTMxZswYy3Kef+e6du2aqF27tnjuuefEvn37xNmzZ8XWrVvF33//bSkzZ84c4e/vL9avXy8OHz4snnjiCREWFiby8vIsZbp37y6aN28u9u7dK37//XdRr149MXDgQMv6zMxMERwcLAYNGiSOHTsmvvvuO+Hl5SU+//zzMj1eV/Tuu++KKlWqiE2bNomkpCSxdu1a4evrK+bPn28pw2tQun7++WcxZcoUsW7dOgFA/Pjjj1bry+p8//HHH0KpVIq5c+eKEydOiKlTpwq1Wi2OHj1q1/EwcbqLtm3bipiYGMtrk8kkqlevLmbPni1jVOVDenq6ACB27twphBAiIyNDqNVqsXbtWkuZkydPCgAiPj5eCFH4AVQoFCI1NdVSZtGiRUKr1Qq9Xi+EEGLixImicePGVu/Vv39/ERkZ6exDcgtZWVni/vvvF3FxcaJz586WxInn3/kmTZokOnbsWOJ6s9ksQkJCxLx58yzLMjIyhEajEd99950QQogTJ04IAOLAgQOWMps3bxaSJImLFy8KIYT47LPPRKVKlSzX5MZ7169fv7QPye1ERUWJF154wWpZnz59xKBBg4QQvAbOdnviVJbnu1+/fiIqKsoqnnbt2okRI0bYdQx8VHcHBQUFSEhIQEREhGWZQqFAREQE4uPjZYysfMjMzAQAVK5cGQCQkJAAg8Fgdb4bNGiAWrVqWc53fHw8mjZtajVwaWRkJHQ6HY4fP24pc+s+bpThNSsUExODqKioIueI59/5NmzYgDZt2uDpp59GUFAQWrZsiS+//NKyPikpCampqVbnz9/fH+3atbO6BgEBAWjTpo2lTEREBBQKBfbt22cp06lTJ3h4eFjKREZG4vTp07h+/bqzD9OlPfTQQ9i2bRvOnDkDADh8+DB2796NHj16AOA1KGtleb5L697ExOkOrly5ApPJVGR08eDgYKSmpsoUVflgNpsxduxYdOjQAU2aNAEApKamwsPDAwEBAVZlbz3fqampxV6PG+vuVEan0yEvL88Zh+M2Vq1ahT///BOzZ88uso7n3/nOnj2LRYsW4f7778fWrVvx8ssv49VXX8WyZcsA3DyHd7rnpKamIigoyGq9SqVC5cqV7bpOFdXrr7+OAQMGoEGDBlCr1WjZsiXGjh2LQYMGAeA1KGtleb5LKmPv9agwU66Qa4mJicGxY8ewe/duuUOpMC5cuIAxY8YgLi4Onp6ecodTIZnNZrRp0wazZs0CALRs2RLHjh3D4sWLMXToUJmjqxjWrFmDFStWYOXKlWjcuDESExMxduxYVK9endeAbMIapzuoWrUqlEplkV5FaWlpCAkJkSkq9zdq1Chs2rQJv/32G2rWrGlZHhISgoKCAmRkZFiVv/V8h4SEFHs9bqy7UxmtVgsvL6/SPhy3kZCQgPT0dLRq1QoqlQoqlQo7d+7EJ598ApVKheDgYJ5/J6tWrRoaNWpktaxhw4ZITk4GcPMc3umeExISgvT0dKv1RqMR165ds+s6VVQTJkyw1Do1bdoUgwcPxrhx4yy1sLwGZassz3dJZey9Hkyc7sDDwwOtW7fGtm3bLMvMZjO2bduG8PBwGSNzT0IIjBo1Cj/++CO2b9+OsLAwq/WtW7eGWq22Ot+nT59GcnKy5XyHh4fj6NGjVh+iuLg4aLVayxdSeHi41T5ulKno1+yxxx7D0aNHkZiYaPlr06YNBg0aZPk3z79zdejQocgQHGfOnEHt2rUBAGFhYQgJCbE6fzqdDvv27bO6BhkZGUhISLCU2b59O8xmM9q1a2cps2vXLhgMBkuZuLg41K9fH5UqVXLa8bmD3NxcKBTWX31KpRJmsxkAr0FZK8vzXWr3JruakldAq1atEhqNRixdulScOHFCDB8+XAQEBFj1KiLbvPzyy8Lf31/s2LFDpKSkWP5yc3MtZUaOHClq1aoltm/fLg4ePCjCw8NFeHi4Zf2N7vDdunUTiYmJYsuWLSIwMLDY7vATJkwQJ0+eFAsXLmR3+BLc2qtOCJ5/Z9u/f79QqVTi3XffFX/99ZdYsWKF8Pb2FsuXL7eUmTNnjggICBA//fSTOHLkiOjdu3exXbNbtmwp9u3bJ3bv3i3uv/9+q67ZGRkZIjg4WAwePFgcO3ZMrFq1Snh7e1fIrvC3Gzp0qKhRo4ZlOIJ169aJqlWriokTJ1rK8BqUrqysLHHo0CFx6NAhAUB8+OGH4tChQ+L8+fNCiLI733/88YdQqVTi/fffFydPnhTTp0/ncATO8umnn4patWoJDw8P0bZtW7F37165Q3JLAIr9W7JkiaVMXl6eeOWVV0SlSpWEt7e3eOqpp0RKSorVfs6dOyd69OghvLy8RNWqVcVrr70mDAaDVZnffvtNtGjRQnh4eIi6detavQfddHvixPPvfBs3bhRNmjQRGo1GNGjQQHzxxRdW681ms3jzzTdFcHCw0Gg04rHHHhOnT5+2KnP16lUxcOBA4evrK7RarXj++edFVlaWVZnDhw+Ljh07Co1GI2rUqCHmzJnj9GNzBzqdTowZM0bUqlVLeHp6irp164opU6ZYdWPnNShdv/32W7H3/qFDhwohyvZ8r1mzRjzwwAPCw8NDNG7cWMTGxtp9PJIQtwyXSkREREQlYhsnIiIiIhsxcSIiIiKyERMnIiIiIhsxcSIiIiKyERMnIiIiIhsxcSIiIiKyERMnIiIiIhsxcSIiIiKyERMnIiI7SJKE9evXyx0GEcmEiRMRuQ2TyYSHHnoIffr0sVqemZmJ0NBQTJky5a776NKlC7766iuHY0hJSUGPHj0c3p6I3BunXCEit3LmzBm0aNECX375JQYNGgQAGDJkCA4fPowDBw7Aw8OjxG2vXbuGkJAQXLhwAcHBwWUVMhGVI6xxIiK38sADD2DOnDkYPXo0UlJS8NNPP2HVqlX49ttv75g0AUBsbCxatWpVYtJUp04dvP322xg4cCB8fHxQo0YNLFy40KrMrY/qvv32W/j6+uKvv/6yrH/llVfQoEED5Obm3tuBEpFLYuJERG5n9OjRaN68OQYPHozhw4dj2rRpaN68+V2327BhA3r37n3HMvPmzUPz5s1x6NAhvP766xgzZgzi4uKKLTtkyBD07NkTgwYNgtFoRGxsLL766iusWLEC3t7eDh0bEbk2PqojIrd06tQpNGzYEE2bNsWff/4JlUp1x/J6vR5Vq1bF3r170bhx42LL1KlTBw0bNsTmzZstywYMGACdToeff/4ZQGGN048//ognn3wSAHD9+nU0a9YMvXr1wrp16/Dqq6/ijTfeKJ2DJCKXwxonInJL33zzDby9vZGUlIR///33ruW3b9+OoKCgEpOmG8LDw4u8PnnyZInlK1WqhK+//hqLFi3Cfffdh9dff922AyAit8TEiYjczp49e/DRRx9h06ZNaNu2LYYNG4a7VZ5v2LABTzzxhFPi2bVrF5RKJVJSUpCTk+OU9yAi18DEiYjcSm5uLp577jm8/PLL6NKlC77++mvs378fixcvLnEbIQQ2btx41/ZNALB3794irxs2bFhi+T179uC9997Dxo0b4evri1GjRtl+METkdpg4EZFbmTx5MoQQmDNnDoDCdknvv/8+Jk6ciHPnzhW7TUJCAnJzc9GxY8e77v+PP/7A3LlzcebMGSxcuBBr167FmDFjii2blZWFwYMH49VXX0WPHj2wYsUKrF69Gt9//73Dx0dEro2JExG5jZ07d2LhwoVYsmSJVa+1ESNG4KGHHirxkd1PP/2Enj173rUBOQC89tprOHjwIFq2bIl33nkHH374ISIjI4stO2bMGPj4+GDWrFkAgKZNm2LWrFkYMWIELl686OBREpErY686Iir3mjVrhqlTp6Jfv353LFenTh2MHTsWY8eOLZvAiMjtsMaJiMq1goICREdHc5oUIioVd6+3JiJyYx4eHpg+fbrcYRBROcFHdUREREQ24qM6IiIiIhsxcSIiIiKyERMnIiIiIhsxcSIiIiKyERMnIiIiIhsxcSIiIiKyERMnIiIiIhsxcSIiIiKy0f8DuMsikzFOth0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(JOF_data_3)\n",
    "JOF_data_3[\"F444W\"].plot(ext = \"mask\", show = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are plenty of arguments that can be passed into the `Data.mask` method which impact the automated masking process. These can either be passed in as:\n",
    "1. A single value, which will then be used for all `Band_Data` objects within `Data`\n",
    "2. A list of values with length equal to the length of the `Data` object, which will be passed to `Band_Data.mask` elementwise\n",
    "3. A dict of {filt_name: value} containing all filter names in the `Data` object, which will be passed to `Band_Data.mask` explicitly by filter name\n",
    "\n",
    "Below we explicitly produce dictionaries of the default values for each filter that are used in the masking process. These can be changed on a band-by-band basis for additional mask personalisation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fce5d2c8f10> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_auto.fits\n"
     ]
    }
   ],
   "source": [
    "default_star_mask_params = \\\n",
    "{\n",
    "    \"central\": {\"a\": 300.0, \"b\": 4.25},\n",
    "    \"spikes\": {\"a\": 400.0, \"b\": 4.5},\n",
    "} \n",
    "star_mask_params_dict = {band_data.filt_name: default_star_mask_params for band_data in JOF_data_4}\n",
    "edge_mask_distance_dict = {band_data.filt_name: 50 for band_data in JOF_data_4}\n",
    "scale_extra_dict = {band_data.filt_name: 0.2 for band_data in JOF_data_4}\n",
    "exclude_gaia_galaxies_dict = {band_data.filt_name: True for band_data in JOF_data_4}\n",
    "angle_dict = {band_data.filt_name: -70.0 for band_data in JOF_data_4}\n",
    "edge_value_dict = {band_data.filt_name: 0.0 for band_data in JOF_data_4}\n",
    "element_dict = {band_data.filt_name: \"ELLIPSE\" for band_data in JOF_data_4}\n",
    "gaia_row_lim_dict = {band_data.filt_name: 500 for band_data in JOF_data_4}\n",
    "overwrite_dict = {band_data.filt_name: False for band_data in JOF_data_4}\n",
    "\n",
    "JOF_data_4.mask(\n",
    "    \"auto\", \n",
    "    star_mask_params = star_mask_params_dict, \n",
    "    edge_mask_distance = edge_mask_distance_dict, \n",
    "    scale_extra = scale_extra_dict, \n",
    "    exclude_gaia_galaxies = exclude_gaia_galaxies_dict, \n",
    "    angle = angle_dict, \n",
    "    edge_value = edge_value_dict, \n",
    "    element = element_dict, \n",
    "    gaia_row_lim = gaia_row_lim_dict, \n",
    "    overwrite = overwrite_dict\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the second time these masks are made automatically is much faster since they are simply loaded from the masks produced from the `JOF_data_3` object. Let us now check that these two implementations are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data objects are the same\n"
     ]
    }
   ],
   "source": [
    "if JOF_data_3 == JOF_data_4:\n",
    "    print(\"Data objects are the same\")\n",
    "else:\n",
    "    print(\"Data objects are different\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this metadata does not change the path that the mask for each band is saved in, should we attempt to load in the mask again using a different set of input parameters, the previously created mask will be loaded in instead of the one that would be implemented using the input arguments. Should you wish to overwrite the mask with the new implementation, you must pass `overwrite=True` into the `Data.mask` method. In either case, the saved mask arguments will match those that were used to create the mask the object points at if the ones the user inputs are different.\n",
    "\n",
    "As an example of this, we will change the `edge_mask_distance` for the F444W band and reload the mask for a previously unmasked JOF `Data` object. When masking the JOF_data_4 object, the code populated the dictionary value for the F277W+F356W+F444W stacked band, although this becomes an issue when attempting to update one or more of the automasking parameters in one of the bands used for stacking. We therefore explicitly include F277W+F356W+F444W into our dictionary below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fcef92bfd00> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_auto.fits\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50 != 100\n"
     ]
    }
   ],
   "source": [
    "edge_mask_distance_dict = {**{band_data.filt_name: 50 if band_data.filt_name != \"F444W\" else 100 for band_data in JOF_data_5}, **{\"F277W+F356W+F444W\": 50}}\n",
    "JOF_data_5.mask(\n",
    "    \"auto\", \n",
    "    star_mask_params = star_mask_params_dict, \n",
    "    edge_mask_distance = edge_mask_distance_dict, \n",
    "    scale_extra = scale_extra_dict, \n",
    "    exclude_gaia_galaxies = exclude_gaia_galaxies_dict, \n",
    "    angle = angle_dict, \n",
    "    edge_value = edge_value_dict, \n",
    "    element = element_dict, \n",
    "    gaia_row_lim = gaia_row_lim_dict, \n",
    "    overwrite = False\n",
    ")\n",
    "if JOF_data_5[\"F444W\"].mask_args[\"edge_mask_distance\"] == edge_mask_distance_dict[\"F444W\"]:\n",
    "    print(f\"{JOF_data_5['F444W'].mask_args['edge_mask_distance']} = {edge_mask_distance_dict['F444W']}\")\n",
    "else:\n",
    "    print(f\"{JOF_data_5['F444W'].mask_args['edge_mask_distance']} != {edge_mask_distance_dict['F444W']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the saved mask arguments are not the same as those input by the user since the mask has not been made using the user inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 4: Loading fits masks directly\n",
    "\n",
    "In this final masking example, instead of producing the masks by using the `method = \"auto\"` or `method = \"manual\"` arguments, we will explicitly load them using the `fits_mask_path` argument. This is useful if you have masks saved in other locations outside of the `GALFIND_WORK` directory set in the config file. We will first copy our manually created masks from the default paths into a directory one step up from `GALFIND_WORK` and rename them for the sake of it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_auto.fits\n",
      "dict_values(['/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F090W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F115W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F150W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F162M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F182M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F200W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F210M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F250M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F277W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F300M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F335M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F356W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F410M_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/auto/F444W_auto.fits', '/raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_auto.fits']) dict_values(['/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_1.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_2.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_3.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_4.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_5.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_6.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_7.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_8.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_9.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_10.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_11.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_12.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_13.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_14.fits', '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_15.fits'])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_1.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_2.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_3.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_4.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_5.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_6.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_7.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_8.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_9.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_10.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_11.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_12.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_13.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_14.fits',\n",
       " '/raid/scratch/work/austind/GALFIND_WORK/new_masks/confusing_name_15.fits']"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create directory for new masks\n",
    "print(JOF_data_5.forced_phot_band.mask_path)\n",
    "new_dir = f\"{config['DEFAULT']['GALFIND_WORK']}/new_masks\"\n",
    "new_mask_paths = {**{band_data.filt_name: \n",
    "    f\"{new_dir}/confusing_name_{str(i + 1)}.fits\"\n",
    "    for i, band_data in enumerate(JOF_data_5)}, \n",
    "    JOF_data_5.forced_phot_band.filt_name: \n",
    "    f\"{new_dir}/confusing_name_{str(len(JOF_data_5) + 1)}.fits\"}\n",
    "funcs.make_dirs(list(new_mask_paths.values())[0])\n",
    "# recursively copy masks across to new directory\n",
    "print(JOF_data_5.mask_paths.values(), new_mask_paths.values())\n",
    "[shutil.copyfile(src, dst) for src, dst in zip(JOF_data_5.mask_paths.values(), new_mask_paths.values())]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we shall load them into our `Data` object using `Data.mask()` as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:galfind:Combined mask for <galfind.Data.Stacked_Band_Data object at 0x7fcef92d3bb0> already exists at /raid/scratch/work/austind/GALFIND_WORK/Masks/JOF/combined/JOF_F277W+F356W+F444W_auto.fits\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****************************************\n",
      "DATA OBJECT:\n",
      "----------\n",
      "SURVEY: JOF\n",
      "VERSION: v11\n",
      "****************************************\n",
      "MULTIPLE_FILTER\n",
      "----------\n",
      "FACILITY: JWST\n",
      "INSTRUMENT: NIRCam\n",
      "FILTERS: ['F090W', 'F115W', 'F150W', 'F162M', 'F182M', 'F200W', 'F210M', 'F250M', 'F277W', 'F300M', 'F335M', 'F356W', 'F410M', 'F444W']\n",
      "****************************************\n",
      "****************************************\n",
      "\n"
     ]
    }
   ],
   "source": [
    "JOF_data_6.mask(fits_mask_path = new_mask_paths)\n",
    "print(JOF_data_6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can see that the paths to our masks match the confusingly named paths set previously. Let's now delete the copied directory to clean up what we have done in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "shutil.rmtree(new_dir)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "galfind_test",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}