{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![HELP-Logo](https://avatars1.githubusercontent.com/u/7880370?s=200&v=4)\n",
    "\n",
    "\n",
    "# DMU 24 - SSDF: Photo-z Selection Function\n",
    "\n",
    "The goal is to create a selection function for the photometric redshifts that varies spatially across the field. We will use the depth maps for the optical masterlist to find regions of the field that have similar photometric coverage and then calculate the fraction of sources meeting a given photo-z selection within those pixels.\n",
    "\n",
    "1. For optical depth maps: do clustering analysis to find HEALpix with similar photometric properties.\n",
    "2. Calculate selection function within those groups of similar regions as a function of magnitude in a given band.\n",
    "3. Paramatrise the selection function in such a way that it can be easily applied for a given sample of sources or region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This notebook was executed on: \n",
      "2018-06-28 14:38:39.108611\n"
     ]
    }
   ],
   "source": [
    "import datetime\n",
    "print(\"This notebook was executed on: \\n{}\".format(datetime.datetime.now()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rc('figure', figsize=(10, 6))\n",
    "\n",
    "import os\n",
    "import time\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.table import Column, Table, join\n",
    "import numpy as np\n",
    "import healpy as hp\n",
    "\n",
    "import warnings\n",
    "#We ignore warnings - this is a little dangerous but a huge number of warnings are generated by empty cells later\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from astropy.io.votable import parse_single_table\n",
    "from astropy.io import fits\n",
    "from astropy.stats import binom_conf_interval\n",
    "from astropy.utils.console import ProgressBar\n",
    "from astropy.modeling.fitting import LevMarLSQFitter\n",
    "\n",
    "from sklearn.cluster import MiniBatchKMeans, MeanShift\n",
    "from collections import Counter\n",
    "\n",
    "from astropy.modeling import Fittable1DModel, Parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def coords_to_hpidx(ra, dec, order):\n",
    "    \"\"\"Convert coordinates to HEALPix indexes\n",
    "    Given to list of right ascension and declination, this function computes\n",
    "    the HEALPix index (in nested scheme) at each position, at the given order.\n",
    "    Parameters\n",
    "    ----------\n",
    "    ra: array or list of floats\n",
    "        The right ascensions of the sources.\n",
    "    dec: array or list of floats\n",
    "        The declinations of the sources.\n",
    "    order: int\n",
    "        HEALPix order.\n",
    "    Returns\n",
    "    -------\n",
    "    array of int\n",
    "        The HEALPix index at each position.\n",
    "    \"\"\"\n",
    "    ra, dec = np.array(ra), np.array(dec)\n",
    "\n",
    "    theta = 0.5 * np.pi - np.radians(dec)\n",
    "    phi = np.radians(ra)\n",
    "    healpix_idx = hp.ang2pix(2**order, theta, phi, nest=True)\n",
    "\n",
    "    return healpix_idx\n",
    "\n",
    "class GLF1D(Fittable1DModel):\n",
    "    \"\"\"\n",
    "    Generalised Logistic Function \n",
    "    \"\"\"\n",
    "    inputs = ('x',)\n",
    "    outputs = ('y',)\n",
    "\n",
    "    A = Parameter()\n",
    "    B = Parameter()\n",
    "    K = Parameter()\n",
    "    Q = Parameter()\n",
    "    nu = Parameter()\n",
    "    M = Parameter()\n",
    "    \n",
    "    @staticmethod\n",
    "    def evaluate(x, A, B, K, Q, nu, M):\n",
    "        top = K - A\n",
    "        bottom = (1 + Q*np.exp(-B*(x-M)))**(1/nu)\n",
    "        return A + (top/bottom)\n",
    "\n",
    "    @staticmethod\n",
    "    def fit_deriv(x, A, B, K, Q, nu, M):\n",
    "        d_A = 1 - (1 + (Q*np.exp(-B*(x-M)))**(-1/nu))\n",
    "        \n",
    "        d_B = ((K - A) * (x-M) * (Q*np.exp(-B*(x-M)))) / (nu * ((1 + Q*np.exp(-B*(x-M)))**((1/nu) + 1)))\n",
    "\n",
    "        d_K = 1 + (Q*np.exp(-B*(x-M)))**(-1/nu)\n",
    "        \n",
    "        d_Q = -((K - A) * (Q*np.exp(-B*(x-M)))) / (nu * ((1 + Q*np.exp(-B*(x-M)))**((1/nu) + 1)))\n",
    "        \n",
    "        d_nu = ((K-A) * np.log(1 + (Q*np.exp(-B*(x-M))))) / ((nu**2) * ((1 + Q*np.exp(-B*(x-M)))**((1/nu))))\n",
    "        \n",
    "        d_M = -((K - A) * (Q*B*np.exp(-B*(x-M)))) / (nu * ((1 + Q*np.exp(-B*(x-M)))**((1/nu) + 1)))\n",
    "\n",
    "        return [d_A, d_B, d_K, d_Q, d_nu, d_M]\n",
    "    \n",
    "class InverseGLF1D(Fittable1DModel):\n",
    "    \"\"\"\n",
    "    Generalised Logistic Function \n",
    "    \"\"\"\n",
    "    inputs = ('x',)\n",
    "    outputs = ('y',)\n",
    "\n",
    "    A = Parameter()\n",
    "    B = Parameter()\n",
    "    K = Parameter()\n",
    "    Q = Parameter()\n",
    "    nu = Parameter()\n",
    "    M = Parameter()\n",
    "    \n",
    "    @staticmethod\n",
    "    def evaluate(x, A, B, K, Q, nu, M):\n",
    "        return M - (1/B)*(np.log((((K - A)/(x -A))**nu - 1)/Q))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0 - Set relevant initial parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "FIELD = 'SSDF'\n",
    "ORDER = 10\n",
    "\n",
    "OUT_DIR = 'data'\n",
    "SUFFIX = 'depths_20180221_photoz_20180612'\n",
    "\n",
    "DEPTH_MAP = '../../dmu1/dmu1_ml_SSDF/data/depths_ssdf_20180221.fits'\n",
    "MASTERLIST = '../../dmu1/dmu1_ml_SSDF/data/master_catalogue_ssdf_20180221.fits'\n",
    "PHOTOZS = 'data/master_catalogue_ssdf_20180221_photoz_20180612_r_optimised.fits'\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I - Find clustering of healpix in the depth maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_map = Table.read(DEPTH_MAP)\n",
    "\n",
    "# Get Healpix IDs\n",
    "hp_idx = depth_map['hp_idx_O_{0}'.format(ORDER)]\n",
    "\n",
    "# Calculate RA, Dec of depth map Healpix pixels for later plotting etc.\n",
    "dm_hp_ra, dm_hp_dec = hp.pix2ang(2**ORDER, hp_idx, nest=True, lonlat=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The depth map provides two measures of depth:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_values = Table(depth_map.columns[2::2]) # Mean 1-sigma error within a cell\n",
    "p90_values = Table(depth_map.columns[3::2]) # 90th percentile of observed fluxes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the photo-z selection functions we will make use of the mean 1-sigma error as this can be used to accurately predict the completeness as a function of magnitude.\n",
    "\n",
    "We convert the mean 1-sigma uncertainty to a 3-sigma magnitude upper limit and convert to a useable array.\n",
    "When a given flux has no measurement in a healpix (and *ferr_mean* is therefore a *NaN*) we set the depth to some semi-arbitrary bright limit separate from the observed depths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_clustering = 23.9 - 2.5*np.log10(3*mean_values.to_pandas().as_matrix())\n",
    "dm_clustering[np.isnan(dm_clustering)] = 14\n",
    "dm_clustering[np.isinf(dm_clustering)] = 14\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To encourage the clustering to group nearby Healpix together, we also add the RA and Dec of the healpix to the inpux dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm_clustering = np.hstack([dm_clustering, np.array(dm_hp_ra.data, ndmin=2).T, np.array(dm_hp_dec.data, ndmin=2).T])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we find clusters within the depth maps using a simple k-means clustering. For the number of clusters we assume an initial guess on the order of the number of different input magnitdues (/depths) in the dataset. This produces good initial results but may need further tuning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "NCLUSTERS = dm_clustering.shape[1]*2\n",
    "km = MiniBatchKMeans(n_clusters=NCLUSTERS)\n",
    "\n",
    "km.fit(dm_clustering)\n",
    "\n",
    "counts = Counter(km.labels_) # Quickly calculate sizes of the clusters for reference\n",
    "\n",
    "clusters = dict(zip(hp_idx.data, km.labels_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we illustrate the different clusters with each colour corresponding to a different group of similar sources (although the colours may not be unique to a single cluster of sources). You can see structures and patterns within the field, e.g. the outline of the HSC coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAHwCAYAAABOjq0vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XecXHW9//HXZ2ZLdtMbJSGNJNRQpAiBEJYmVcGKgghiQa/e37VcL1evCnZRkeu1IagUQZQmcEFQAUMooYSWQiC997K9zcz5/P44Z9q2LFd2Ttx9P3mMM6dNPrPzePiZ7/d8v5+vuTsiIiLSPyXiDkBERET6jhK9iIhIP6ZELyIi0o8p0YuIiPRjSvQiIiL9mBK9iIhIP6ZELyIi0o8p0YsIAGY2y8yeMbM6M9tpZk+b2bFmVmFm15rZejNrNLPVZvbfBdetNrMWM2sws9roPT5lZomCc242s/bo+uzjwng+qcjAUhZ3ACISPzMbBjwIfBq4E6gATgLagC8DxwBvBzYBk4DZHd7ine7+qJkNB04GfgIcB3y04JwfuPtX+/JziEhnSvQiAnAAgLvfEW23AH8FMLPvAn9y943RsdXRoxN3rwMeMLPNwLNmdq27L+rDuEVkN9R1LyIAS4GMmd1iZmeb2ciCY88CXzCzfzGzw8zMdvdm7v48sJ6wV0BEYqRELyK4ez0wC3DgRmCbmT1gZnsD3wOuAS4G5gMbzOzSXrztRmBUwfa/R/fwa81s+1v7CUSkO0r0IgKAuy9x98vcfT9gBjAO+G93z7j7z939RGAE8B3gt2Z28G7ecjyws2D7R+4+InqM6ZMPISKdKNGLSCfu/jpwM2HCL9zf4u4/B3YBh3R3vZkdS5jon+rDMEWkF5ToRQQzO8jMvmhm+0XbE4APEQ6o+5yZ1ZhZlZmVRd32Q4GXu3ifYWZ2HvAH4DZ3X1jKzyEinWnUvYgANBBOh/uCmY0Aagmn230J+CBwLTCN8B7+UuC97r6y4Pr/NbM0EACvAT8Gri9d+CLSHXP3uGMQERGRPqKuexERkX5MiV5ERKQfU6IXERHpx5ToRURE+jElehERkX6sX0yvGzNmjE+ePDnuMERERErmxRdf3O7uY3d3Xr9I9JMnT2b+/PlxhyEiIlIyZramN+ep615ERKQfU6IXERHpx5ToRURE+jElehERkX5MiV5ERKQfU6IXERHpx5ToRURE+jElehERkX5MiV5ERKQfU6IXERHpx5ToRURE+jElehERkX5MiV5ERKQfU6IXERHpx5ToRURE+jElehERkX6sLO4ARERE9mi/+mX+9cwT4PAj4ovl/0AtehERke4sW1q8Pe+ZeOL4ByjRi4iIdGf9urgj+Icp0YuIiHRn6VK8YNO7PXHPpUQvIiLSC9kk/+i6zbHG8WYp0YuIiLwJm5atiTuEN0WJXkREpJcceNeyl+IO401RohcREenGzug5221vwF8B6upiief/QoleRESkG4OjZ4ueHWg89Hja5zweU0RvngrmiIiI9JIDxy1+9p+qlfzPFKuIiEhJVRS8NiAJHAykAPd/jsl2SvQiIiI9sC72lQHt23aVOpT/EyV6ERGRbmQIu+sLH7XARqDiT3+MMbLeU6IXERHpRmP0bORb9vcddzr7YTgQBEE8gb0JGownIiLSjaFd7Buzawe3nnwuZDK8a2cTo8d0ddaeQy16ERGRbtx6yruAfLd9HbB9/ARIJKC8nAe27uzp8j2CEr2IiEh3NmwA8l33y/bdP9xvUUd+MhlLWG+GEr2IiEg3Prr0RSDfoj9m08rwQDS1LtHQAA0N8QTXS7pHLyIi8iZc9vcHuP/QYzjk9SVMyzTBi8AVn447rG4p0YuIiHRlTbhKnZGvdd9MWETn/MXzAVgItIydydva26moqOjiTeKnrnsREZGuFFS+y06tq6a4hXwYkEmU8/KDS0oY2JujRC8iItKVluaiTevwunA7kyiDTKYUUb1pSvQiIiJd2bYNyHfbQz7Bd6xyP75+Bcz5e4kCe3N0j15ERKQr5RWkKE6UHRO8ATO3zAXgjaaJHFii0N4MtehFRES6Ul/HbbPPJaDrBN/x9aq9x7N95fLSxPYmKNGLiIh0Zecu2LYtV+8+m+y7WujGgQM3rOaReXtepTwlehERka5UV3HZkue5f5/prCVcgx7gdzXvpCHabov2BcD4rZsgmWTT356KI9pu6R69iIhIV4YOgc3wkc3Lcq12gMTrr3NPVAOfIOCyJx4MX447jvdveo6yOnA/EbOuVrIvPSV6ERGRLgQHHYotW1Z0f34zYeJn87Kic1NAfXMzKWA4MOjXN1L+iU+WLtgeqOteRESkC+llK9hF2JIPgFtsGI9QkRt8VziXfsr+xxIEAYtxBrOLZJCice2aOMLuRC16ERGRLlS8vpD7a96ZX6muvZ19n34Ed8/ty7X2177KzGUjWTUFmvZuoHrLcJY89iTHfnRSLLEXUqIXERHpxkmJO3lyThLIcFJN2C2PVXU67x5auZgMHhjz1sJp5ZAO6kocbddi6bo3s6vNbIOZvRI9zik49mUzW25mb5jZmXHEJyIi0gY0NcNJNWGSZyvUHQdtjCqaW98G2JijeK56O0EA/7OxirPWJDhg6Pq4Qi8SZ4v+Onf/UeEOMzsE+CBwKDAOeNTMDnD3PbOAsIiI9FteOZ6m1fvTlLsT74y8/j/4+75ncUpNTcGJztCJj3HcvGm8VL6BO45v4ePPruXhbfDhOALvYE8bjHc+8Ad3b3P3VcBy4O0xxyQiIgPR3tVMeGYneApIM+GZnWx7YRVNqU24kbtPHyTgjK01gDHymLXMePYEnuU4rql4MMbg8+Js0X/WzD4CzAe+6O67gPHAswXnrI/2iYiIlFSL701Ly0YmzAu3M5kMwfVnc2b9+wkISJAAgwYC2LaNfVhPSwtAI9eygLYNZ8cZfk6ftejN7FEzW9TF43zgl8BU4EhgE3Dt/+H9P2lm881s/rZohSEREZG3SjlQWV2Zm0eXLEvy4m+WcsPtX2dqpTOlMsOUygxHVMJjV3+OzYznkGdOAIbyRY7je1MejvkThPos0bv76e4+o4vH/e6+xd0z7h4AN5Lvnt8ATCh4m/2ifV29/w3ufoy7HzN27Ni++hgiIjJAlRk8+tFpBAQ4ThNN8N4r8ieY5brva8/+LDcl55FdxHblqfN494r3xhJ3R3GNut+3YPPdwKLo9QPAB82s0symANOB50sdn4iISGtbM97czGOXT+fRy6fxzOWHwx2/yp+QrYsbGEHNbE46CSANpGlvT7BvsJOV8/8WS+yF4hqM9wMzW2hmC4BTgM8DuPti4E7gNeAR4DMacS8iInEwM07/4xaC+XNpampi19eu4NwpVxac4dAAv2xy3l/RxLQ5k/kDo/ntwTfjL9/Jp8puYcqD72PVL+IthRvLYDx3v6SHY98BvlPCcIrUPbyQ8qij5pFVOxm/z37MfM8BcYUjIiIx8cBYc8ww3jF/P2zBFoIpV7Lm2OGwCshWxxsGn3bAE8B4th9ofJmVLMhsjzrxYd8djzL3Nz9nxfiD+OhZp5X8c+xp0+ti1bq9mQocw0iS4Mwpo9iya1PcYYmISAxWrlxJbW0t644fxdrjR7H++FE0NzeHB4MAsKhqjkGmAjD+8kYdDy35GQHluYI6O6sO4YH9xzBtw+useOO1kn8OJfpC28NyhRYVR0hiDN1XVYJFRAaiZ57/Hcv/9kOampoIgoCWlhYWP/jt8GAiAXh+VZugjbOBS4AzWMz3+CTrK49k4ZDzuKHxGKaMOhoDdqxaXvLPoURf4JZBt+HRfwCLK1fy+iNtBEEQc2QiIlJql7/3KxAELLz/Gzz/xy/z6n1X5w8W5QWHIOBh4NWxAT897VrSVPCbtlO4t3E6ABvXPEPzoCEsnnpMKT8CoEVtiq2GWw/6a24zkwbajyaR0O8hEZGBJmEWtdy7EASQzKZQg1QVtx8/mRufXc13HvsJlw2aX3T6flNmsb6yMtcBUErKYAVmrKgOkzvRj7VrpwGggjwiIgNP+c5Xej4hu7JNBj5FmlQqFa5uB7RS3EW/5Z5bMYsjzSvRd3LE3CpmPF7F4XOqmHZs+EW1tbXFHJWIiJRakw3m/G+d2uWx0QsXQSKaSJ90rm9roby8nHLCLv2zJ8zk9Aljcw8SQ8OR+jFQou9gsbeQIsVWWqjKVMcdjoiIxGXvQ6n87bFdHrrSxrD/S42QSjF4bSP3VI7i8ZdWcXd0fPms82kE2oGlJ13APZNuLlHQnSnRd3CIV1PuFezl1Sx8KdxXVVUVb1AiIlJyN/7+Kto77oxa5duev4Frpo/gntbB3DpyFI208Zvcycbql19m40kXsPqkCwiCgLNWnVXCyItpMF4HheMujjwWXnkBRo8eHV9AIiISj6Gjcl3xOWYcf9EPAHjvsMLbus6g6Bkg/fozrHk9PyBvDWs4pE+D7Z5a9AWOHje7aPvIvU8EoL290286ERHp5w4c/jYa7/tUp/35IXU7yY7Ge3K/J4qOnrnPfqRJ4zgttPCZCScTzx16teiL3HLwI3zEz8xtb968DTDq6urQCnkiIgPLyNNPp/2hG6C9HcrLw51mOE5AwNzxr2YXryOVyl7VBlTxTDmcPWFc7r2WUVnK0Iso0RdaD7ce/Beqr4F6oOzfgeePjjsqERGJwdIpI+CBR+HC98MBoxg+6kDqlszl+7//Ji/QDld8nXc/dzcPb4PWihMYBPyULTzLeprGnc2yqVPhyfvgpAsAsJYWqC79IG8l+gIz3qhm0T7NNF8JZQ4HP1rNR2jhPel03KGJiEiJnbV8B49MH8OoP95FIp3m/De28toZl8Efv0wN8MP6ev503PvCk3ftghda+FcmA5PhiXqumgrMfnd43D3fK1BiSvQdzHii4NdWEo5Bo+5FRAai8pYWPrhiZ2577nNPM/SIWRx/0TUAVMxbTYLGDldlKBr+lp077w7JZF+G2y0NxivQcTDe28bNooZB0BBTQCIiEpuKigrq6+vJZDJs2LCBE06YzT4v1JEhw8/bNpIIsqXxsqIV7bJW7IRMJryBb1Z4I7+k1KIvcMvBj3AJ78itXreVrRgJqvdW4RwRkYGmtaKCoWVlmBnjx4/nj3+6ganlIxlz5d38F3Duu35QcHY24Uer2pGGqaPIjdbLrl8fAyX6Qm90XtTGOYqKiooYgxIRkThUQVF9+vec91FebEvw8I6Xwx1BAImChW2KlqyJ9heWvY3pHr267gvMWF9Npj38XjIZ4Nq9AS1qIyIyEN2w9lXa29txdzKZDK/dtZSX7vla/oSyzm3lh4jq37M93GGWe2ge/R7iiKfyA+92HbuF1S+M6+FsERHpr0adeCK//uZn8zui6XEfmHIlAD+lnjCp51vyl9g3+f3xNwH7s8Q9bOfH1GWfpRZ9gf3GHcrR+86mnhbW0c4pe51OQn8iEZEB6QOr6zh7wUqor6dswUrOWrYJCMdnP9rpbAcC3uHnQLqei+atAoo78+NK92rRF1i/cTFj9j6RU/Y9E3fnxY1PAQdTWTks7tBERKTErv/dVzkDOHt12A2/qTJcbf65y/cP7/He9Cod0/cfOJY/vLAD8E7NRO90dmko0Xfw8uanc68zmfB56FANxhMRGYj+dsTUrg+YcfIkeGJNbgfZNL5k9kQAbkmlwsHcMa1Dn6V+6S5kb6ck9NcRERnQhp2eb4PPvOiUMGkHAaRSnHzakV1eEwQBTU1NJKMCOWYW6316teg7KPzhtfCV8PmFFbt4+7RR8QQkIiKxmJgYzIxZR1FWE6bKdDoNacMTCUgkaGtr6+KqbVRXPkV1JYxafD51h0/AY27RK9F3YIn8F3Lo0WkWveAox4uIDDyNV1/LrtdqiwrdnDf+CLYvqev2msmsZspj4RLnTxzUyF7Ld+SOpVIp7jlk374LuBvqnO4gRViiMCAgkQmLG+zcEMQZkoiIxCAA2HZP2NXrDvX1PLj+x51K4xS6+ugLclc/aHeFvQCE3fl3HbR33wfdBbXoCyxKNnNoJpxHnyBBQ7KZNGkmHzwi5shERKTUUqkUI2ZVUl7+l/zOJyrwqPSNd1EC5+YX17LqtHnhcR/FzZPzyd2DeBqNatEXqoFFtNBOOytoYfUp8N9DtEStiMhAdMnqOmqXncv27dtpa2vjga8+y8x3nUs9KVpI8b6hnfPDHKC1Nayts//jJ3LdvLNoa2ujsbGRnzx3Tsk/A6hF39mpsMwyYU9NAMc2QkLD70VEBpzrf/dVZi4bjF/1deq3wsyL3knVlR/gPb+YGZ7QZQPdOPjpE6PXzQD84qXzSxFut5ToCz0LHFewfHAaaqhi165djB49OtbQRESk9AY1L2DQle/Lba8A7mnIL04TvH1iUWPw4meXser053LbN098JNfBH9cUOzVVC8xoqsYao0VtUsB11QQEPL2mMe7QREQkBn8/YirthFXtmoE3jphKmjSOs5MUFz+/lkwmQxAEPL9kKZAku0p99jqIt969WvQdHPpcNWZRTYSj4KmX2jh5qgbjiYgMVI91qI534dAALByJXw5c8sK6aBR+BZAhO0u73MEK1qF391gSvhJ9D8zghhHwmbgDERGRkhszbArb61aGG1EL0Dsk6oCOXeNJbp5UsOJdzMVyQF33vVJVVbX7k0REpF95/wUfh3TByHqz3i1Kk03uqVTRbt2j3wNlpzyWl5f3fKKIiPQ76XSakd//bVFlvFw52+g508V1bgZmeHl5fnW0GKnrvgcxjp0QEZGY3fj7qzrtS2QTQ3YwVw8se17MyUQt+g6ylQ4Bvr0ifN62bVt8AYmIyJ4lCPKJ4h3jOxx03MNH4A7RCna4E8RUGU8t+g4Kf3h9bRpcsiIsgygiIgPLF049N/fagesefwiA/3ljB+ue+Ant+53JT86bHd23z7fuv3XT84BTzXa+9MFToboagHVz5zKppqZU4eeoRd+Fwp4Z0D16EZGB6OePP4RjBJbglcb8oOydyf8BwDb8teDswqVuKoBKFtSs5pLt+Xv8cSR5UIu+R9memeHDh8cbiIiIlFwbcN3jD3baP7J6OE1Wx38fVNPFVdlk77z6Kux/NLFPsVOLvheS2XssIiIysLlzyHdfCl/XnNLjqYcdFlbHA2IdkKdE34Psj7CyMnV8iIhI6K/ZSnlBT8k74IC5xzMaNOp+T5RN8O3RT7H6+vr4ghERkT3Ceyd8sThp54rpZKvbQ77KfSNQzy6ATZtKF2QXlOi7kMlAUxO8sSjcbm5ujjcgERGJVSV7c8+6a8ONTZsIgNt/dWM3S9WWA+UsOv41RgCUlcV6n16JvgtlZTB4MBx2VLgd19xHERHZM7SxJXwRBAz+xd0kggA+8KnoaGEqzbb4X2LGs8cD8N66RL7UagyU6LuQ/eGV7aEZNGhQfMGIiMieI5EIe3kTCWprawta9IX34Q04GiijDWhraIDGRnbs2FHqaAEl+i51HDexrU2j7kVEBppjLri6aPvtF34PgHOH/wWAzxwwqocsGtZfuWvER3hw4lCoquLON77UR5H2TMPJe5Bt2U/fe2i8gYiISMnNv+9qjv/QNbnt7du3A7CpuevFbIoFwIvsX/saK9aMo/XY4/n0Cb/tq1B7pBZ9D/aAZYRFRCRGE+ftoLa2lhV3/Joj3wiTwpgP3Zg/oYc88cbpLazc+73YSTWxLneuRN9Bdq0Cd7h0xWQAbpy7It6gRESk5EYAmUyGw5dkqJl0AXev/REEAY8UnmTQ2sW1i2c+z20TPgOVlbmBeB5T61GJvoNE9Bcxg1unrQbg/IPUdS8iMtBcfuq5pK54mdaPz6ftEy9y7jdOxRMJOGJqfgxeBv4VCJv22Ucb350H37ppPjQ1cf6qWnDHYiqco0Tfhex3kU36iYT+TCIiA447T742iCCAVAqeSVyYG1ufywpJqJk5GdhAeF++FSjndo5j5anPQUUF908bHWt1PA3G60G2l2XIkCHxBiIiIiWXyWTIHHA+T1sSyqA9KpfadcrelzD9DwLCaqq7dgFTKkoSa0+U6LuQTfDZ54qK+L8oEREprWRZkvqr/6XLY5lOLwqnYQ9nMq/xeuO44lHd6rrfc3T8LrSojYjIwLOCdLfHkhAmiy7KrJyIs5rpbNy+MTwn+4iJEn0Xsj/AGhvD561bt8YXjIiIxGKfcSd1eyyAMFkEm/AOc+yexoAFfGnXLMLTvOi51JToOzCgpQX+tByGa7C9iMiA1d1A7MJ0vYy5/Iw0tx8/ObfvPQAczdcvP7IPo+s99Ul3oboa3jONXA1jLWojIjLwvL7htW6PGUCQZrpdyD2Us3PnzlyyT6VS3PvihuJzIbbpdbEkejO7GvgEsC3a9RV3/7OZnQF8H6ggXND3S+7+eClj89z/5NcqGDpUTXsRkYFm2KxN8FjxvmxrPgWQKIOWcOe/Lt0J1BWcuT26oGB09wAcjHedux8ZPf4c7dsOvNPdDwMuBX4XR2DZ7yL7XF5eHkcYIiISpy4yZDZVl2c3ymFN3UpWnPIskCZsIjbx5KFvRBfEPxhvj+q6d/eXCzYXA1VmVunubfHEEz5rep2IyADUw8o1YXowKHO+sN94DlsNq057LjzmsP/js7jqOKCtLSyDG6M4E/1nzewjwHzgi+6+q8Px9wIvdZfkzeyTwCcBJk6c+JYGlk3w2Vvzzc3NVFdXv6X/hoiI7OHauz+UhoLueOOdy46BZZV0KqfT2opHjcV+VwLXzB41s0VdPM4HfglMBY4ENgHXdrj2UOAa4Iru3t/db3D3Y9z9mLFjx77l8QcBXLZ6AgA7dux4y99fRET2bFOfr+L8b5xatM8JG/o/u+JbYJkw2Wec12fNj46ELcWrLj8SMhkurr01tgSf1Wctenc/vTfnmdmNwIMF2/sBfwI+4u4lXzbOov9JJuD2qeu4eMVkkskuKiKIiEi/N/WJKr5w6rkArDixhfu/9Tj/8fEZfLr21yz/yuOc87OvQfJgLhsEq057NnfdrleeonootJfHN38+K5bBeGa2b8Hmu4FF0f4RwEPAf7r703HEVvh1ePQjTN32IiIDU2Fb3CvDDLFixQpu/MNqEsBqFoIZc9iLOXOgvR0eXwIjR0JFMp9k42zTxzXq/gdmttDMFgCnAJ+P9n8WmAZ83cxeiR57xRGgo0VtREQkn6Qrlh/B2VfN4v6/t4T7r3sHk5MXhgcTNcyeDRUVcOrB0NQUrniXJL5781mxDMZz90u62f9t4NslDqeIAXj4bNFgPHXdi4gMTF7wvIU5ZOqiofhm7Pr8X+EXl0UnOFP/fmLBlQHPHTYP9qZ4YZsYqARuBw65n2+Z3Hz6eH+NiYhIib34IuD5qnbAyGkjeea/F+ZO8W/+vqjwyv4sIZxHn2HVafMYPhyGR8fipETfBfcw4a9fH25v37491nhERKTEtm4hQyZasMZJkyKdTlMxK39K5ddvyLfW02n+dspOVp02Lzcor7IyTPvuHZe9Ka09qmDOnqDwd9ekCcAKdd2LiAw4mzaz6tQU0I6Z4e6UUcZxxx3Jk0+9AoAzh6icC5SVkcj36ofHHVqJeoVj7L5Xi74r2WwffS+pVCq2UEREJAaJMBEU3rptbW2lqqoqt/3eH/+2U3nbol56g6I5W22xFHlVou+ocOBFEH1hI0aMiCscERGJQ0UF5vmGuJkxaFAV9fX1YX4AGDSIT2f+GN3vdYIgU5RDAB6G/JvEVApXib4XVOteRGSAaWwkoPMS5RUVFbQDT0w/A4DPNRNOnn9sU6eB241Ei9+kUrB+bl9H3C0l+h7EPCNCRETi4gHunVPkyMRILplyJT9PH8UD/3IRBw75QDR5fp/w/IK8UQacBuHxCSeXKvJONBivCwWDKAFYv7OR/UapaI6IyECSyUAiUdzoa//rXbTfdyMA9x4wA+ZuZFBreCw1K+ydz54/CFhhFPb/lyz2QmrRd8E97Gm5fPk4AF5aUxdzRCIiUmqJRJpUKncLnlSqndTs2TjQBNz8if+Cw6sLzqfofHcYB/kBezF1EyvRdyGRgPJyuPWgjQCcMm1YzBGJiEiptbWVUV7umIWPh7//FI3b96b2mrtpv+Zujgg8rHULgLN2a5g7Eol8Xi9a0j6mFr267nuQ/U7Ky8vjDUREREpqF+FsuOpqyyXtc758Kt/YZMyxgmb7svw1L2+BqePD2/RmkEhCMxSsWx8Pteg7sA7PAIMGDYojFBERickmYNSq8D69OwQB7FoKQ4N8v/wMfzKsZR/50q4TwIrzx1AorqATAyX6Qiuj5+i7SKRji0RERGK04djzmVo7iGQyLIGbSDhjD3G+MOcuNgV3ssnv4iE2w2LID7VPgxcvd94A+dY/MOquq0v4KUJK9AUmZQ6H5UAaaID91hwad0giIhKDoD1f+jw7P97d+QZ7kQbqgUmJD4RD67NJI5w1X9wj3OF9T2//Y5/F3B3do+9gAjOw1eHX9DV2AuX87ysbeOeR4+MNTERESubM4/cm82pAJpPJrXeS2baNOadezYSiYvYZwlQ6FLpYuqZopZSWFl4Aavo08s6U6Duwgt9iX2coc2hl+tiqHq4QEZH+ZtHWrRxKkrJkPk2W77U3rKPoXvuYg76aex320J8XbkSp5I3CN62qYlpx9fuSUNd9B8004zhp0iSj32KbG3SzXkRkIFlclwY8XKbWACtorxcsZLP99dMIgnCwXtMbZ+aXp49OnQ7kTlj/E/528Yul/BiAWvSdVEe/tsooYxM7gApqDtor3qBERKSkKnYkyZAp7rrPZCDIhPPmANyp4Ayalp4RJXbH/emimXQpyE+s3+/fYN3XYPpvSvpZ1KLvwq0sZA4N7MOouEMREZFYBLRPnE1ZogzDSM95IuzGtyhtBnOK5saHrf3Oc+Ubcydk+wPm9FXA3VKLvguXcjiOd7lykYiIDABtGTawElsJ7A/lNTXQXnA8eUrR6dkU37E2zsjCnek0THiyT8PuihJ9Fzz6bdZKuFJBe3u7lqoVERlA6kdBZu9VVN58JGVzwlRpAKenwcpzLfR2IFxlPtxOpaByUP5+/j3MyGf+sjJYdzNM/0rJPgeo676TtSwkTZpdNHBmVKU4CNSyFxHQWOOxAAAgAElEQVQZSIKyBM/WXkj9h16hmWbSpGnd/wUgCUETBOnwcfpetAEp0sDdJBKAhzPr5zIaxs2C7dsL3nlFyT+LWvQdTOQwDGMkZTxBwMk08fL6RmZOUxlcEZEBw8BHGS/ZB+FyoAlOSr4CWx0Sg8Nzsvfdz9yLjJeRWfUibW3nUV4els6Z7TuYveF66g3u9c+ELfsJ3yr5R1GLvgvZrvtM1KI/aB+tXiciMqBkZ1Vn++BbgI1AdEsXiJaoC3L7xpbBI0ujyzz//Apj80XzY6BE30ELLZ3m0Y8covvzIiIDyeBGwmTvQAuMcGAarOAhCNqjlW4ybOQeyAwF4HPHP8yXdqUAz8+nN1g58cJwKTwtU7tnqCKsgldGGZvZDAyONyARESm5ZmBEXcHC8tEsrE3AJu7DosVrHKAJGAq7du0iLHD7VPQuBS34aC4+606C6aW9T68WfYE1UVf9WhYCsA/7ABqMJyIy0DgUT4tPJvL7OzzfFfX+/uy+NeFOs9y8ejeDtSugthaAiX0ZdDeU6AtNX0yGDBM5DMhPr2tubo4zKhERiUM2mxd2uU8aU3T7fjwjeD8tRZe1Rptu8DJDeQ8Pw9Ch0NoKc7dTauq6L7QS1k9fnNtsoQWWHUJ1dekXIRARkRiNAHaRH0Bn4TK1kE+cBqynliVXXsoZP7uN7C+Da14sJ98d0MoHZgG126FyMMweU6pPkKMWfYFJmcNZtayFNGnqaeGAZccCMHfu3JgjExGRkrJsgs8uSZtflx7CcXqthMvQnvGL3zNm/WUMOegrDDnoK/znce8jrHLfzpXHvY/RAOMnwaBB4f38ElOLvoOTeTu2LPwyX+M54JB4AxIRkdLb1cUI+Sj3pwmTZxn5yXZm4fFMAE1NzfzncR8ke3Yu0VZVsfYglcCN1Ro2MTEagAcwnaOAVmpqamKLSURESq+KNFBRNCE+O8Bugr2/+L695dewTSTg9brnoS7fYb5x0uTwRXbkfYmp677Q9G1F8+jX7/86AM8991zMgYmISClVphJQuzWaLx8wqG5b7q77HX5XuL68OwSbw53RwcZ6WHriu4EMYfbfRePE82KbQw9q0ReZtOxw1kxbkBtouX1r+FxWpj+TiMiAEmQ4gwzJum3hJgGbgjCfnwBs4m6I5tKP48LwGgM2fheCgKUnva9gDn6Bdb+F6V8t4QdRi76TicsPw5bDr5bDUfW6Py8iMhAFlRVsnPp31vMoG3iUjdMehwT435xJyfexL3vzfWBc4sKi635Y+QxtHWuvZFv/AOHQvJJSou8gO4/+O8xgA68B0NQUwzBJERGJ3zRgetTzbmBnRl3wyRp+krwwbMVn83o7LM5UkslkOr9P1HU/kR+VIOhiSvQdlFGGYSRIMInDAVXGExEZaByn8P/6i3vg24tL40Uj7selhnL3tHv44cM/zF9kFo7Qcw8fq0oSfhEl+g7WRK14CO/JiIjIwGNYD+PnkoUngkFZHWwc3EDjZnj9iE9Hx6yL56P6JuAeKNF3MJGDc68tGkY5eLAWthERGUjq8aJWfHHSbymugw+kR8AH+SSvH3E78Mdu3zdJzVsYZe9oOHkXnOJRkscee2xMkYiISCyGQCoZtt0LquACcCyv8ILPLjr9+uF/IJ1OU1l2B7lb9B1H3LuzpS9j7oZa9B1k59EXJvvfPv9yjBGJiEipVTdCIhW+Nito0afgukUpaGnJ3XevXruWTCbDZ5s+TCYD0+Yc2vnCIICn/kjzlJEl/yxK9B1UUUWCRK7bHqCxsTHGiEREpNTGDU+yoi6/HTbOHcodTm7l+SczvPC3Jl74WxMPvz6KZDLJTwbdQjIJq05bHF3TYT36kz4IfK2UHwNQ132Psq36t5XHHIiIiJRUfQMcUpvAx4SDsnNNv9XbmRLAszNHsXTeTgCmzxzJp2o/iDks4LzwPPeiRXCy3fjDSxR/IbXoe+HFNo2+FxEZSNrTAeM5NZ/goxf7O2z6O7TSzgEzRzF95kgeJRUejs5Jpbp507a2WBK9WvQF1rCAiRyWa8lnCEdUnDhc69GLiAwohUVvcnPmHTPYpwaOKXew9tz+MbX5Vn9ZGRAEeDKZ/6Gw/kEm8APWcmQpoi+iFn2h6bCGhWTI0Ewz6/ZfAsDL3f06ExGRfikoT7KeR8ON7OJ0ZrS2wqqpYxmz/tKw6Z5KMXb9pWFCd0inYe7c8NyiGXhlx7AOYNT3S/o5QC36TmwqrLewaM63l4f7musqYoxIRERKzTAyU6LSOAVj6jYcNBZrhRZgzOaPkT1sDmTC1vzJJ8OqguVtAdhnH2AurH8CmF6qjwEo0RfLAMl898tV0+HiZfC5s4+IMyoRESkxx3mt9qLcdh1tnDjqHp7mQqiEpvEn05TId4qPXfeRooza2NrKkMJia0EQTrWbUFOC6Iup677ApJWHF23vtyycC7lhZ0sc4YiISEzqBzsBAY5TNziAUVFJXHf4w+8Yu+Gy/Kp0W7fmrnOH5mYYkt1ROPLeDFpbS/gpQkr0BcacPIKJyw7LPbazHYDnXo+jlpGIiMQlkYCGUVA/yrEk0P4q7nDg63UceEQ4hW7MhssYu/5SxrT9e9F1gwfDZdtvCndku+43rQ6f73lbyT5DVrdd92a2oBfXb3P3097CeGI3ZvYIXpr7JEyHI/Y6AZ7eyFHj4o5KRERKaUiDkXIHgwpPkMkcBXu9AdzNzMGtUFt8fjvwzOAbgfuLK9/OuRsYA7NngzsTS/YJ8nq6R58EzunhuAEPvLXhxC+TyXDU7JNwd3Y+WQ/ALS/t5KrJk+MNTERESiYRQIWFVVIdZ/FTc3nbe+CEwWHX+7Z9fsOozR8jQbgcfd24m6jZ+FEWjghb++6Eo/Jr3hdurNwFU0fB7DEl/yw9Jfor3H1NTxeb2b+8xfHE6hYe4SPJMwEwM6qPrYAX4JwJmd1cKSIi/Ukq4ZRFSR6DQ2pOBn4PhKVtqzfPY+eEW/MX1NZGx6LrMxCUlUWFdAymjsKCeIqvdXuP3t2f2t3FvTnnn8rB8NKWl0in09TW1nJX3VwA/rIx5rhERKSkmkZAMwFpApoJaOpw/FB+DetWhRPnt29nbMP/45fD7iCdhrZ0OGp/3R2/Jp1Ok06nWfvSS1BX1+W/1dd2O73OzBZCh3VboQ6YD3zb3Xf0RWBxmPFYNYtO2cYieyzcsQFYCh8+4cBY4xIRkRKrhdQoSAE4DNuVHz1vZjww7jwO33gVbMpf8kr988waCuXkW/Yb77wpf8L00s6fz+rNPPqHCWeY/z7a/iBQDWwGbgbe2SeRxWTG4wXlbqOKePsOq4wnGBERiUUlMGhnvtM7Qwbci1ekIyqWE73+FdfxGT833O/OpA99vHh6HYTN5BLrzfS60939y+6+MHr8F3Cyu18DTO7b8GLQTPjNpYGoIF5lpRK9iMhAUkWSVDSPPkNAGYmwrG3RinTQIY3j7gRBwJGb/pzfmb2moQHmbu/z2DvqTaJPmtnbsxtmdixRVUDCdPimmdnVZrbBzF6JHud0OD7RzBrN7N+7e4++MGnZ4TCY8K9SsDRtW1tbKcMQEZE9QAUJEiTCJF90B9vCh3W+r21mJJNJFow/Dzq0/kmn97hR91kfB35rZtlCPw3Ax81sMPC9f+Dfvs7df9TNsR8T3jIoqTVkmLjsMNaykIkcFq1e11zqMEREZI9RkKzdCxey6/4Kjzr0O57UtBV2vdXx7d5uW/Tu/oK7HwYcCRzp7oe7+/Pu3uTud77VAZnZBcAqYPFb/d67NX0xGTK5pWproiT/1XsXljwUERGJW7bT2jl01B/yXfBGLnt27LqHsFWfLZfrhM9kMozjE6Sq7+vroDvZbaI3s73N7DfAH9y9zswOMbOPvQX/9mfNbIGZ/dbMRkb/1hDgSuAbb8H7v2mLks2sn76YNdMXsHb6Qm6evgyAA6d09VWKiEi/5mVhg96NxVs+kGuhG9FqdT1dCrBxY34efTLJxvFPsK1PA+5ab+7R3wz8BcgWgl0KfG53F5nZo2a2qIvH+cAvgamEvQSbgGujy64m7NJv7MX7f9LM5pvZ/G3b3qI/XQ0spoU0aZpo4Y2aFA9jzFu023BERKSfGVS7Nbyv3tTEWY07w527SfBZBlDRYYlzM9LTRr+VIfZKb+7Rj3H3O83sywDunjaz3ZaKc/fTexOAmd0IPBhtHge8z8x+AIwAAjNrdfefdfH+NwA3ABxzzDG9/NP37Nzrj+ehTz7LGxbNq3sczsZ59jzNoxcRGWhmkaYsSvCN9KLB1yETXdZ2Fzfb/ws3Ok6zK6HeJPomMxtN9BHM7Hj+wZmAZravu2fLDLwbWATg7icVnHM10NhVku8ra9jEOTccl9tOkeJ7tFJdXd3DVSIi0h9tnjanaNtSYRu3aIxdF83Morn2ha9jSva9SfRfIFy8ZqqZPQ2MBd73D/67PzCzIwn/RKuBK/7B93trTN9G87JhVFFFhgwbpy+BZVMYNmxY3JGJiEhMLDuAPpks3s6d0PH8cIe7x9qSz9ptonf3l8zsZOBAwo/zhrun/pF/1N0v6cU5V/8j/8b/1bbpy/MxRD/V7pi3nA/NnBZHOCIiEpPu8rR7+DiK9/KS3dOL9/HiQjsl1tN69O/p5tABZoa739tHMcUuuyxhQLjS0MThQ2OOSEREYhetPpfN2TXjz+TkiZ8gyolkMhlY/4v8PPqu5tLHoKdR9++MHh8DfgNcHD1+DVze96GV3qRlh5MmnUvyq5e1A5CsXxVzZCIiEiczIFGcMi9JP0A6nc6Vvf2f58+NzrV8T0DHsrkx6LZF7+4fBTCzvwKHZAfPmdm+hFPu+qWpy47Ktegn41zVm5GWIiLS72QyUBZlSXewTKbofnyiDH76wnmdrsu36Cm4OL5k35t59BMKRsgDbAEm9lE8IiIie4ROuTnRm5TZzbUx6k3Uj5nZX8zsMjO7DHgIeLRvw4rHQ596Fo/+A1jJSoju04uIyMASdPy//2RyNyd0mFq3h+jNqPvPmtm7gdnRrhvc/U99G1Z8/nzF87nX7e3ATXvFF4yIiMSmoQFGjcqPp/OoOz6byoPCFn72nGiEfS7fRy/iHHnfq34Id/+Tu38+evTbJA9gjeH3kkrBj2+aTjjNX0REBhR3NvrFtLaGOaG9PXzeSpjTs7frL+ZCMpAbZH/kpofD0fdAc3P4Pm3Zi2PSbaI3swe7O/ZmzvlnMuOxahgCCYOKcvjyh5cR1gcSEZEBJWp8L2+6mMU7LmLplneTAOonjWHZhFEky8K17a4c38TOib8ruvQ7z4WrumcycNm6n7P57luwdS+XNv4CPXXdzzKzB3o4bsAhb3E8sUoRJfusDMBgdsYUj4iIxCkaPW9AdSsZh6RBMgjIpKG9DMJla9rzlxjU1NRA+yNkq6efNGEsGZ5hA0eX+gMAPSf683txffvuT/nn8cZpzcx4pDr8qwSwuq4ZWE1mW2XcoYmISCk5HYbOV5HJhAPvg0SCbz+b5OsnpPBEgo6DtrNd989t2IeZkzZjk9azYeJnSxZ6Rz3No3+ilIHsCQZfD4uuaM7vaAJum8wW34PmSYiISIkUzoevzQ26D0hiZLBEostVa+fOncspJwJsxjNwXBLWxDiXvveTAgeAGo6HXxH2U2yAs247FoAZew2ONS4REYlbGW1t4SvzNF+dmT/yxeu/nd9wqJkwAoATJjpmcGeMrXlQou/kLI7lnJuO4+wH3856FgPwwft2xByViIiUSm1tbfSqcDpdgkQiP0/+W/Py52/pcP2/lv0marwbz7bsHc633xNH3Q9USZIYRoIEkzgcgE+NiDkoEREpKTeKCqjV8ghbKqC2tpG2bNOe8EfAbZ8KW/RP73UzGPx8ys9zeX1m9RYuXf+LN1VV762223/ZzE40s7+Z2VIzW2lmq8xsZSmCi4MX3HHJvq6pGR1XOCIiUmK1+TyOZe/Rjz6G0c3w5z/Dn/6U4epZ2eOQHYy3cFgrAJ9Z+ZmoRe+A8wTA+q+XIvQu9eYnxm+AHwOzgGOBY6Lnfq3wl5y67kVEBo6xw4eSisrbOk6GgBE7FjCkHcaNg3AydoHGcBD3uvXrOr1X2tKcBMCcPoy4Z71J9HXu/rC7b3X3HdlHn0cWo2yCzz7P0Ow6EZEBo7wcKqK2vAFlJIB2fKTztrcBBPmyuBn40m3fB4fvNH+nqBQuQBnl/G7iZ5m+6gul/yCR3iT6v5vZD81sppkdlX30eWQxMoqnQHz1bHXdi4gMFBXJZEF2LM4HVVUVXHhhOHArAMaXvR+OmJUrgZt9hHXt84P5lpUk8q71JtEfR9hd/13g2ujxo74MKk7ZtegBmmgCWvjRnH7dgSEiIgV2NrSRT49Rz27U3qusrGTBgpb8kUQ7x5xYExXAJ/fs3mGg/ex3lyDyru020bv7KV08Ti1FcHFpo40naeBCALbwSu1uLhARkX6jPZUp2GoFfs+iHfsAEAQBS5fCkiXh0V8CFwZR7blsCbq10NIC7hl+udfHsd/XlCbwbvRm1P1wM/uxmc2PHtea2fBSBBeXQQxiNsN4kMHAvtx2gbruRUQGinlLGwq2qoCLIH0y7s4ddzSRyQRs2Rl2zH+aynCevMEQwsVsfn7qz6mshD+tr+CKe49i39lRidyYptL3puv+t0AD8IHoUQ/c1JdBxckL/gvv1WsknojIQDIEYNeOaHach88NDfz+943RGSlqTgxXrwu76sN78RVUgMGuXbuYt3483+B+/j57DLXZ82LSm0Q/1d2vcveV0eMbwP59HVicLDeAwtF69CIiA0sGOIEMNDWFrfWmJs4knFx/0UVDcuclsy+iYjhjGcs0pvHVLf/JVrbyysSzaQIqBz1V0vg76k2ibzGzWdkNMzsRaOm7kOJVOH8+LIE7mVQq1fNFIiLSb7QDz4zcizNTjZxVv52zUk3UUcfFFw+loaGx0/n73nsvELbol3s4vv7dk9p5iqcIqp6G5rrwxLZOl5ZET8vUZn0auCW6L2/ATuCyvgwqToVT68ZxMNDCtmYlehGRgSI7FG/DtEeL9t9+e33n88zY9P73M3btA8xgBgtZCBVlPG53gcMFTT8hAO7f+W84Tn2Tw/RSfIq83oy6f8XdjwAOBw5z97e5+6t9H1o8UqRwnICAv37oRX5zxiC+8Pd+24EhIiJdcae1tRV3x6PX4Y32qOFnBV33kTu4g9wdX4dT7f0EwAPJK8IkjzOi01V9r9sWvZl92N1vM7MvdNgPgLv/uI9ji0UFFbnX595xPLOo46OHaD16EZGBprKyMpfzysrKuPjiQQDcfvt2cAg6pIaTW85k9ObBbGUrp3o4Oe2BIf8GQUCCBMNjGpHXU4s+uwj70C4eQ7q76J/aycWbI04MP+aZB2h6nYjIQGNm1Nc3sHZtA2Vl3bWLHWjlwHkXcPAr5/ACL/DhMz6cP9y0nQtafhq+H/E0Grtt0bv7r6KXj7r704XHogF5/c8TMHLWUBKJBO7O9idVKUdEZKBqbW1l2LChDBsG6XSasrKyoiVqC1euO4VTOHrWUt62/1f51tJv8fB+R+LuXMDtNMUSfV5vBuP9FOhY276rff3Czqfygy0CAsBobtE9ehGRgeaee1J0WqkO2C96bgfCjvE0/7LPu0kth6PTEyAICIe2wX2D/600wfag2677aBGbLwJjzewLBY+r6TwGoV/Jdq9kn8u77bIREZF+yYyLLhqW2xw/Pn+oJppwXg5AK/dVGK3A/Hc8DKkMd/7Ac/f2sfjHePWUwSoI78WXEd6Xz6oH3teXQcXJCwZL/I4mYD3NzZNji0dEROLR3AwXXxwm+yCARMIxd6auyXff/2kQOAleOfMRGurrufPacL+7k8kU1Mx3jy3p93SP/gngCTO72d3XlDCmWBUOlvgIQ/g1+/H6hjRnxhiTiIiU3qrWi8I1bQBI8eojvy46/vlZMGbMPIIALnj2JsqAP0bHXvjX+8MW8k1fDHfEmOh7UxmvOVqP/s9m9nj20eeRxaiwBO5jVHLd0pgDEhGRGOTXlM9kkrnXUMZVJ4bL3QQBnLzhBBLAjfvfmDtjCPAY7BFd971J9LcDrwNTgG8QlgJ4oQ9j2qN8mjfiDkFERGIUri3v5FemyfDNp8NheidvPBWopBmgQ3XcadFwPSDWhN+bRD/a3X8DpNz9CXe/HOiX69GvIVO07ThqzIuISJjo26NHBaOJSuB6mPyrgPTQdOEVlH8nvxV4fMvX9SbRZ+cWbDKzc83sbcCoPowpNos+9QIZMrmFbf7yifnAZG47b3jcoYmISIzCUfSVQCWDCEelV4QHgHAq2qeXXRGeG10z7r/CSqsxrlAL9G4e/bejBW2+SDh/fhjw+T6NKkZ/uWJ+p3076upiiERERGITDZ7zqCXe0NDAVSeGo+0zGUgmYV32PKDNobKg1Z4b2O0e3unfk7vu3f1Bd69z90Xufoq7H+3uD5QiuNhkokdY8IhRw1QCV0RkILsw00omA62t8O2bFgEwIg2k07DuI9QP/QlYFyVmMpnwEaPdtujNbCzwCWBy4fnRvfp+5dzrj+ehTz6bLwd0d/h0yUM7WH1SbGGJiEhMzAzcuXf+fWxZshAoqBhXBiTLYcKtUFvbuem8F2HTH6C1lUNaG3ht5F6lCbxAb7ru7weeBB4F4v1ZUgLn3HBc7rXjfI96ejeUQURE+pVsV7w77znmAn5YdLe9necLz2m4Nzy18JytcFbt1tzmI3vaojYFqt39yj6PZA9iWMGXFf8cSBERiUHUmseMH86/D6L1T6ZGY9SPCgDzsD9//GWw7vGiomuOExBgGC20wMgpcXyKXiX6B83sHHf/c59HE7M10xcwcdlhuSTfSivgnDZSyV5EZCBxomZe9H//H/jAIOrrorXM3aHpf1maIPwxUFUFW7dG1xW2+o2/jtgn9qI5vUn0/wZ8xczaCKfaGeDuPqzny/45rZm+AFtm+LRo/MSqCTy2K+6oRESkpIzciHt35847WznrxHAcevYG74INTwCvM5blYXr37M+D3M+Eglfx6c2o+6HunnD3KncfFm33yyQPYB4meXe4dNVELmVF3CGJiEiJNUfTrtydlStXAvB2zyf5dBqmMZaxLM9fZBYVzc0n+T1Bty16MzvI3V83sy7XnXf3l/ourJisBKZE1Y0Nbp+8lkXlaS7YEHdgIiJSSqmRUGdhsh8zegrThr6Hbz5zD/kWu3P9lJWMoaDF7l7Ugm8vfMM9cfU6wgI5nwCu7eKY0w/L4E7KHM6a5QvCiYStsN/mQ2mZPp9DTum3HRgiItKF4TsT1A0JE/2IxgS7Kiq6OKsGY06HfflUX4kXd9vHlOx7Wqb2E9HzKaULJ15rpi9gwrIZ2Orwi1jLQh7Eue+hes2jFxEZQMyc4Y0JzMJZWIlEV3e653S8KP8SCu7bR1rqobr0JdV76rp/T08Xuvu9b3048Vs3bVHR9n3LDogpEhERiUu2Ez47IG/1vL/mjoQCRtPhPnwmKLg2NLR1Vv4w0MwzfRRx93oajPfOHh7n9X1oMcovQcyVn9T6dSIiA0+YHrM16kfPmEl4192BNF+flcl1y+fa8cnomug/gIbcFdBcXfokDz133X+0lIHsUTz/3NAQayQiIhIjd8fd2bHof/n6CWCJsFhOPV1MnSvops+ugkr1Mx2XqS+53tS6/y7wA3evjbZHAl9096/2dXCxMximcXgiIgNWtkV/yGkf55uP3VhwxPEpnRN9tiWfa9UXDMBz91hWsetNEfezs0kewN13Aef0XUh7gOz3EITfkdajFxEZYNqKJsdFg/Gy+5yTSXUuhFPWRds55qp40LvKeEkzq3T3NgAzqwIq+zasmGV7X5ohlYINO7QevYjIgNJcBxVjc4Pxamtr+fossCjZB0B6A5QXXuPFBXC9YN8evR49cDvwmJl9zMw+BvwNuKVvw4pZO7ACqIbycrhyXtwBiYhIqVXt2kp7eztLlizhQ97OuMQHaAIagfHJCztfEISj7tuAJrK18oun3MVhty16d7/GzF4FTo92fcvd/9K3YcUoACqA/SGzAp58Em5/54i4oxIRkVIaOZbZ2zeTaKqFfUZzzfP3wknvYlpBgk/RoUUf7auMHq2li7ZHvem6B1gCpN39UTOrNrOh7t7/xqPXA0Oj7haD9ER4Zsl0Lh/fEnNgIiJSauHStAW8uHp9M4dSzeJ8d33Gqf/SZ3Mt92zXvXVxbSnttuvezD4B3A38Kto1Hriv+yv+eU3acjjsDHtf2trg8rXjuJW1zJ/fFndoIiJSckNzr2ac8cmoG74gYU/8D3YRdgS3ASSN0T/8GRnC4jjbGRkmlEwmPH8PrHWf9Rng7cBzAO6+zMz2+kf+UTO7mrCO/rZo11ey692b2eGEPyqGEf79jnX3kvSA5Erg7gy/jNtZyKLpbcw6UKPuRUQGmiPO+jCOExCQac+ESTuRjI6GCT8z8XfsBAgCBq2/FCfN2B/+DAPGAm9MeRK+9Hmoim+udm8SfZu7t2dHDJpZGW/N6nvXufuPCndE730bcIm7v2pmowlveZRMUQlcB46BTG2mlCGIiMgeIDuG3jAqWlshyEAySS4FBg5JBzLgxevPZ5PklFVzWFV1VXR+AF3WzO9bvfkXnzCzrwBVZnYGcBfwv30UzzuABe7+KoC773D30mfZwt6V+XDpw3HXNRIRkVIqLGPrOJnBgwsOWvEz6f/f3p3HyVWV+R//PFW9Zw+JBAgJAZIAabYQHEACIWFTnAEHHEVxZUAFF0BxdHQcxuWnM86A46goKIgbKgqCK0vCoixiCBAaSEhCQkIWsnen96q6z++Pe6vqdnV3FpKu26l837xquurcpZ/TN85T59xzzwm753Px40NpSPT+POxcov8MYRf7c8CHgD8Ae2JWvI+a2UIzuyWabQ9gCuBmdq+ZLTCzT/d3sJldbmbzzWz+hg0b+ttt9zgEZ8A/Dx2Y078sY3gAACAASURBVIuIyOCUf34eh5Sn2LhxIzzxx55JvtA6t/BzB2ykcBgAqwFao8ZiVzLjvXaY6N09IBx8d4W7X+TuN7vv+OuJmT1gZk19vM4HbgQOA44D1lJc874KOBV4d/TzbWY2p5+4bnL3Ge4+Y+zYsTtT113m6bCap5+ueXBFRPZl656+kyP6ynwOh+LQ3c2Cqst49du/Z/HE09kGLB4xkY5bH2S/VeHKd15XV9aY87a3TK0B/w58lOgLgZnlgP9z9y/u6MTufuaO9onOeTPwu+jjq8Aj7r4x2vYHYDowd2fOtUdFX9pS62HevBZmaj16EZF9RnwmuyAIoGYqs7oXs6iPfdfn4JXVC5ledTOXffZGhtwwlzX51n5XF2sZTo1ZYl3422vRXw28iXDU+2h3Hw38HfAmM7t6d36pmR0Q+/g2ID8C7l7g6OhZ/SrgdOCF3fldr4uHLwuMr/5mMukxZY9AREQSVJjM1iCbhmmnndbvMPSblj3NwwBZ2NTiTL10DlM/cAZTP3AGkz58LjXTzi5X2H3aXqJ/D3Cxuy/PF7j7y8AlwHt38/f+l5k9Z2YLgTMIv1TkF8y5Hvgb8AywwN1/v5u/a9f1WGQ4DVvKHoGIiAwSafIt/L4z/XeB944/lgWZy3g7sHjseHJAB7D81gfZibvdA2p7j9dV57vQ49x9g5mVzvq3S9z9PdvZ9hPCR+ySEb8enQCvceqphycUjIiIJMkw0qQICPrd5wTgEWC63cw/HWA0fvk2lua77rNZrK9V7cpoe7+9+3Vu27vFH61Lw2cv30Lzc22JhSMiIsmIr0WXy+UA4ysHjulx//5fX13PDZNPom7lSo6squeXF27jwkvnxM4BL936IJDcCnbb67o/1sxa+nhtA44uV4Bl10l4ZXKQqwkfihwzRovaiIjsa+LP0qfT4Yx4X1q9AXfH3Zm3Onq0u6uL57uW8FMu4fZvweIJU/JDvXjp/dftkRnmdke/LXp3T/e3rVI1zWln2tz68EMKNtMFnfDi8q286U3JxiYiIuWXb9W7Oxh0uvO51RvC5+bdwZ11K5/GGcIPuh+mEeC67/KSWWEfc09snnvYuQlz9h0r4PnZHTxPB89XdfDa6dD4owa+3LTDI0VEpELkk/tWcjQT4IQt+Hxi78GMXwMHHP533Bcs732yaJ8kKdHHNC5rCKcwnA3MhMbvN9A9uZufX7Bf0qGJiEiZtEcD70aSZkSUJlPxOepLpsC9cvJJAPw0/W56Ta8WS/JJjb5PdijgINT4UEPxw4GwlGwFD0gQEZFSwTBgW8/BeNlsNkzanZ3Q0FDsus8zY01mDQBdmQy1NTVhuRe7/nO5HFUJjMBXi76EL3XIApshNyS3nQcqRESkEg3bFg5Rs9hjWKlUinPOOSdava5vn04/BEBtvKs+et/R0cHzD9y854PdCUr0JXyy49WO7+e8SBecmHREIiJSbrno3jxAC051FWQyGdhYMr1MbKDdj+3isKyzs9f5ljxyGzBsIEPul7ruS6Ri330aaQjn6Dus90UTEZHKlY7lgvA+fRvV1dVMXTafxQcdVOy2j42+b86vUle/DmdYj679Y8/5cBmj70kt+hLeUrwwIziQ5mmT4NZbE4xIRESSEESt+i4CCOpYtmwZZ+Q35h+fAwgCXl76V4YPGULLwQdz+GWX09nairsTBAG2dUPP/ctMiT6maU47PtzJ/7eapawa8jwceWTSoYmISBllCUiRwjCqgK6ss2Tp0p6T3+Rb7Ok0h04+iSAIGL5qFZ3ABdkO3tK8kfNaNpH04qfqui/xwuxYN30Wzv3Biaw9b5+bO0hEZJ+Wjgbi5ee6b+t13z2W8qOE3xntMwRYsOaRnnuOumjggt0BtehLZaKfARyytJ4/jfgbVO/WGj4iIrKXaaHYu5sjwOvqgPhyKCVd8W1t1NfXx44vrHjORqK1zvUc/SAwl3CyHIAUrDisg2OHnZxkRCIikoQh0FKbT9Xgy9rDn9ArYa9d8gQAF9fV8e9RWdeR/0hXgvfl49Sij5lKQ5jsO4F10PjnBg78XQ6y2aRDExGRMhrRloLNhPlgc1To9Dm3ShdwMfDQ/kf1mGSnVFKL26hFH7N4TjuN9zXAo4RXZBn8hSeZ9s+3wKIXkw5PRETKxDBGksbbw/TcWlWFO/y4q2RHdw6JpsAdu/hdOM7at4f3gJ1iV39+BH56OxPuDBQl+riHoems9uLnmdB4ZwMdyUUkIiJJiprhLy9YQCoF7691vpvLQXWUPvNd8w6BTcRYwRvuSNFyVharru4xv30ul0sk0avrPqYx2wDrCftmOuHIR2thVDdDkw5MREQSYVEif8OUKcXCwgI3RqHN7s6mcZ8lQ4Yqqgr7mFnhHDX5+e/LTC36mKY57Rw1t64wv3EXXaTXvkbrpEkJRyYiIknasmVL+MYMcjlI55N9NGCvvR2CgBqiZJ7J9JgX3xIclKdEH7e493P0jY9MhF/+NrmYRESk7HoMqjPIrHy4+LmjA2qqKQ6vC3h1zXMcvq2meNzWrTBuXM+pchOirvuYxlcbwuGTDuTgyEdqSS9phFt+kHRoIiJSRhb7D0qeqNu2rbBXKMX4ySfROWlS4Zipn31ntEvyj9ipRV+i8S+x9egdxpBi46ZNyQUkIiJlM3/+az0+lz4u5w6XvTSfm8ePj0ry9+mdodu2xPZM/vn5PLXoS7g5rANf6QSpgLkfeoIR13466bBERKSM8rPi5d+Xuva7n4eWjbD21fA9cOAry/s4UVJPzxepRd+XA8IfjsPqZEMREZHyaWgAWih22eO0lCT6fE/8tT/7BgCZSY0ArJlYHLhtpTsnSIm+D/lvb110MXVxNe072F9ERCpD6RNw5vn16IuCAL7+4S/So3ve4bdf+3z846AYiAfquu/F3MiRo4UOGnINHLrkeJ794bVJhyUiImWwsXTmuz6StBncVddNuApalk8QHTT92O0f6w6drXsq1J2mRF/KocqrGOEN8DK8xmtk0GA8EZF9wdptfZeX5uzOTrirLuCuuhwzq4Fcjq8OH9lzp/j9eXdwZ3hw9h6Nd2eo676Ep4oXJjM5Q9OsVyCXYEAiIlJ2pc/R48Wn5rNZuJhq6Iy1ldPw2Ycejh1rYR9/Ol1I+MM6T4XlmgI3UScceBrWajhOliwvnZ7l2JvHF1cuEhGRfYbF7sG3dYb5OpOB2zMzuLM2Q9gKjGbGC7pKjvPYPPgOHVvoAOr4RBlrEFKLPuapNY8wffLMwufxD67hNZYlGJGIiAwGdzADui1M2mlobt7CXSNHAWHR1ueeLjnCoLsb6urCee8bRpP1RznhwA1lj12JvsSCNX8ufgiAgwAtRy8iIrF77nOWrqC6ZkW+7U4AvFS6f2dnmOhjx25o3sDBZQg1Tl33MScceFoxqQcwcfkxjB82FaoTDUtERBIUEIRvYtPZls6DUzo23wC++b+xgnCPlW0vDkiM26NEH3Mbf2L6xJlMP2gm0w+eyfDTG8gtMo5tOS7p0EREJCE5KIyaxx2CgNtXUcj2BlhJ5ndgyvKHog/FYzu5vnyBR5To40bCj7iXbrrZyEZunzKPKqpYOuy1HR8rIiIVIz4Qr9pT4X32trbwlU7z0Os9cQKrnivRx5x390kwBX5+xIP8YepTDLkZ1pCiLbM26dBERCQBxYQfwNCh4QuYl56wwxnvDMJn8fI619HS8NiAxLk9GowX8worOO/mk3qUZYFNLUclE5CIiCTH40/TB8TvxOcKy9f2XrTG4z/Tsefm6w/Y8zHuBLXo4ya3sIHiow8vH/o0qSDg1dFzEwxKREQSYYYVWu09f6YKu8T36b2n57clON+9WvQl2g9byysWdtWnl6bZlkox5BtPwckJByYiIgnynm8Lc+F4r0QfZ/lH7PJT4Ha8Cco8P4sSfUyLwfD49ToEtlYvpO3StyYVkoiIJKDXFLiFdBmWBxTzfb77vvA5tucQn0Oqo3iqFjQFbqJWzm6HbsKrk4Oms9rDq9miRW1ERPYlFvuvdAsWS55GoUXfVyd/W91fCl8KMgANsUnZykQt+hJNb46tPr8ZyMDR774bFn0tsZhERCQZjhMQcAh/YwUz8oWkLYtFTXh37zEoL98b4Dh0rKSt4dHiCdtPodxd92rRx70S/VwP/Ah4Bo49+BSe+80HEwxKREQS45DyFLPq4P0184H10UQ5ThA11c2MxQ+HK9dxbCMALVG7fvjqd/U43fDyRV6gRB/TuLQhnALpDcB7oHFJAxsebYWu5qRDExGRBJiFXfXTzrkcgPfVrgTAqe4xDe4Rs2Yx9RMfZcqRD4Qt+f/5byBFZvNl+RMBJPIcvRJ93KvQ+FADjXMbaJzXAAdF09y/49dJRyYiIklxIOcQhEvMk0qBO80lbUADnhjxeXKk4JOfAmDYd2+DXK44DW57+W8DK9HHTOw4BpYStupbYOLqRjpPb+WkHR0oIiIVLZ1O9xzVZkZVFT1Ws3GAujqqcIZGZbnw4MI+w7lnoEPtRYPxYsacPhIePgaWF8vmjm7i2B/emlxQIiKSuEwmA8CdmWLZ1M0TWDsiGtxl4Wv8K6vYSoY0acBIkwpb9Ok0mKnrPmm/YyH7nTaix+vcu06CTGbHB4uISMVwvMdo+urqcL3yf4wvWz77DYCRX8POSNF+5920UkOOHEaK9SefDd2dxWPanylTDYqU6GM2HbmGtra2wsVdvXY1a8/sIrjs8qRDExGRMjKsML1t/jn5xsZGfpw5EQjXQal7YB2bo/nwHdjqDqxhCN1UUxN+AXj8PuBz4Uk7OhjOFWWvi7ruS/x6xl+KH0bA27edxh0//CCnJBeSiIgkKN+qb2pqAk5kWg00dUMulWJkEGCpsM3c0N1NMwdSS3fh1v1kYEndDeGHmppwUrYyU4s+bkXJ53Hw/JonYUwSwYiISOLc87ffOfvsswE4tT7aBKSie++YUVNbywigLt+ap8cM+QzrPo3WsgVepEQf07isIeyPAQigcV4DL9AJ7S2JxiUiIgmJkjjAfXfeCUBHNHe9AatffbXQdb+65FAHlt76IHSeCoTpJdB69AlbAo3eUPzcDo3DGmh67y3wzJeSi0tERBIRn9p22KpFMOkMftQdDsFLk+MXLz0LLz1b2OctfZ2kNpzfvqP+UWj/V+AHAxt0CbXoYyYSPUffBSyF8eumMfrr34LvaQpcEZF9Xe2oq8M3seR/9ezzCu+vOaPPNF98jr7jcwznoQGKrn9q0ZeYyDH4ynDUfds9H2ZzVS14EndVREQkaWaGu9PS0kI33bEk7yw/ZA3B2ju4ZvbbwZ1cLtfXCcJjcjmG8BB/fuQ9zJld1iqoRR933xELwzmKDSxlVF/wTXjrZXD7H5MOTUREEpDvuh82bBgbRv01XwikmLRiAt08BWs/g637LOkNn2fqO/crHpvf1wzSadrqHwWeKncVlOjjMrPaWcc6AgI66WTuZc/QdGk7x77jG0mHJiIiCXF3Ft77PWqrg7AgepwOh9Tokrnr9784f1T4iJ31XM9+5mkvDGSofVLXfUzj3AYWXP4KWDil4ZBF0PYwbGt8GU7Rk/QiIvsqw8NZUvOJ24EO2FoLbxj31WIXfaFr3/o7VdklkujN7DrgMmBDVPSv7v4HM6sGvg9Mj2L7kbt/tZyxNc6LjbrPwniOYe7hS8oZgoiIJMzxwkPwQRCEiby2FvKz2ZpDAxwLrKXYxe/u4bH5qXHja9kmJMmu+xvc/bjo9Yeo7O1ArbsfDZwAfMjMDilrVJ1ESxIC1bD22f9HZuimsoYgIiLJyk+BC5BK9ZUqnXO4hTVRb35+XzPDotZ8oevekm3dD7Z79A4MMbMqoJ5wssCyzVbTNKc9/K0pCn0ddctXcMjqY8oVgoiIDDKlC9sAdJLi7rM+T1sbPRa/yZQughbvzs9kSD+ysSwxxyWZ6D9qZgvN7BYzGxWV/QpoI+wJWQn8t7tvLltEK6BpdjtNtENqCE2ntzP6gpuovuqqsoUgIiKDi7uzrRPa2opJvI6VfOqmLzD5wEMwMxYtWMCib/wfNTU1xePiJ+l4E1RXk0SzccASvZk9YGZNfbzOB24EDgOOI0zq/xMd9kbCTvMDgUnAJ83s0H7Of7mZzTez+Rs2bOhrl13WuKwh/O2zoen0DTTOa2ANsOTXmjBHRGRfZWY0VMPPghMAGAoMX76BDFC3cCutra0cMX06R1z1MVavLk6EW+iwd4f6R6F9Iw+f+1K5wx+4wXjufubO7GdmNwO/iz6+C/iTu2eA9Wb2KDADeLmP898E3AQwY8aMPTbaofHh2GC8NOwPZbx5ICIig00QBPw4M4Nw2tS6cGGaSdP5HtO5kVW8+v1b+zyukJhyOYZlTscNWtv/AbihHGEXJNJ1b2YHxD6+DWiK3q8EZkf7DAFOAhaVNbj4V4YctP/mcnhDUNYQRERk8AgH2m2FoA6AqQBkOYlVvHlSuI/Tu71pAB2tUFVFBmjt+g3DuacsMccl9Rz9f5nZcYRpdQXwoaj828CtZvY84d/oVndfWK6gmua003hfQ/Hrz8uwDaCzcztHiYhIJTMz8FHQAre/bzTZbJaZLzRt/xiidmPDMHCns/5RqA1oSe8jq9e5+3v6KW8lfMQuMU1ntRc/zITzuIlnq59ILiARESm7QgvdoDvt/Pulx4X32jevoqqqqo/2e2mr3gigMOLeAVIprP3PwEUDGXovg+3xuuS1AAHQDWffGg68IK2uexGRfYlF/+FQDWRzOQ5fv5wgCFi3bhMBYaooPSYuBRAE5ILinsP5lwGOvDcl+pjGuQ0wjPCvUgP3XfIUnYe1wnr9mURE9kUGVLnTuHU1Zsa7f7uVq5+ANOFrR1PhPHv/TTTdfxMW7dvSsI903Q9WGXpPgbsiu5j9mJRYTCIiUn75bvgwkYfPz6fT6WJXfHxq27pa6OwsHOM4FnXdH3vOh2Mn9URmyVNTNWbxnPZwLj4HAmha0c6ymRk2jViedGgiIlJG+a77cL766sKa9AAEQTjVbZS0p37osh7H5LvwU9BzCtyEpsJVi75E05tjg/FaoPGxKpqOzyYXkIiIJKqjK8u77t1SLIimvLWSxJ1vyRc/x+RXt0uAWvRxr0Q/txE+vT8EDn5lOoxNMCYREUlUEJQMu6vqeXe+70VvYhJewU6JPqZxaTQF7jBgKpz9/RN4eujT0N2RdGgiIpKQdDpd/NCr+z383BX9dJyAAPDiV4GEV69T132JaQ/XF94vm7yAdWc43P4kfC7BoEREpGxKu+Cbm5t7bC28i94+9TSEw7jDY9KEXwxy+Z20TO3gssw7yJGjjQ7qvZ73LD6b6T9elnRYIiJSJqXPw9fV1fXYihsdHZDLwbp18ATj+znP4KAWfYnDPfxeNpQq6IRmmsl9+tMJRyUiIkkws9736M25dfGU6INzzvF/JPNQH8eGJxjYAHeCEn3M2jSMi3XLZOozLFnzDKO+fiN88AMJRiYiIuVSmMo2+vHAw0t77XP2J58tvO/IHkTVDT0OAWJd9wlT133MplnttNOB42TJ8tKsLC96F+kdHyoiIhWi8DR89Kz8yTMOoNgRH85x19UV5vBMFu67bl7s2OKrkGDjz9InQC36mMa5DTTNbgei1epWQ2NXQ5+LF4iISOVzdx6fvy76lE/WOe790rz+DunvRHsyrF2iRB+TAabNK466d5wJq46muf9DRESkwuRH3TsOBqWz08/KdHHNGecB+XlwnEXP/l+v8/TqDe7qgh4D+8pDXfcxi+e0kyWcBc9xXqATw9C8eCIi+478qPv4CnZxx7UND7cXZrY1jrz6471PVFVf8jmZtrVa9HFzYfHsDPkFDGiF3JKM7tGLiOyj3B2yFLOlQzVV5PO/RfuUTocLMOYL38qfJPy5ehFMbCxD1D2pRR/TSAPMBbqAZdD4RAMcESjRi4jsw65tq4dmoAOubakPp0XvOZF9r2MceO0Ll/Lsvf8Zdtk/8j2GjL28PAGXUIu+RCMN8CjhVXJ4Ivc0U5MOSkREEpFvqV9Lfbi6KdC5oQsfD5Zv0hfa9rHjop+nnHYPqeAe/ERo3fAhmFimwGOU6PuS/6aWgfHU0jlmTKLhiIjIwMtPdVsYjFcYKd+GR5PcAjyPE2QhPgW+l4yqj1Y7p2H5NYWyWmDTwIS+Xeq6L7GW9vDqdAKrgOtv5NmNGxOOSkREyqUwGC9qzZ/+0SW000lAwGY6+fRVi0sei7deiR7CBNvJHQSEI782TTpjoEPvk1r0JQ7whvDqRE9ANP3qg5x60S2JxiQiIgOvO1dDYR6VPIOXWk/hgmseKyT3v/8BPDce4qvT5peqLS6IE3bnZyetAq4HoJbr6er1sN7AU6IvFf+Wdgi8+SdvpAUlehGRSjc0ehqudArca36zkdsvOLmw3/ffmqPrr89RVRUNqDeiGVXDZ+fzx+fopqPuL4nPd6+u+xLeSXhxcxCsCPjjW59MOiQRESmD+vr6Qos8PwUuAO60t7fj7uRyOS75XTO5XLjJosF4YSs+nHUl3/Wf78z34o6JUKKPaZrTDvWOpxyvcrond3PO72bo8ToRkX3Atu7i9Ggee37uY698j0vv6+Bd92zhkt+FA/aO3Tyhx7HuTlXUSZ4/dtMpgJUuelt+SvRxK+D52R08TwfPpzpYenrAvZfPJ9jhgSIisrdr6wyfn8snasejhvhhUYs83kaP9ikZhGcYa6NJ18Y+VgPZ5OdWVaKPaVzWEPa8zAZmASmoX9DXVAgiIlJpVq7roitq2hUH1cEpV18c7WE9fnZHz9V79PmIq8JpcA+ghoCA9d/pZlj36RAE4c38fH9/mWkwXonGh4vPStIO+6+ZQldy4YiISJms3txBbdT+NYxOctSSYhET6P0EvPH0Gzt4/Df11GwdzRXHbeb6B1t5S5T206Q46IosHdSy7eubYcyYMNEnMN+9WvQlfKmHDzyuBx/ivDb5paRDEhGRMqipguaGgEw6oLkhCGe4AYYV9oh30zvV1bAS57Jx4dS2Q0edEKX5DLAFz39taIhGelVXl6EWvSnRlwgmB3iN4/s7bXQwZslhWo9eRGQf0J1Nk0pD+wiwaqjrChN0f+O0bnwUJgDbaq7FHd49+aloSw0wGoDlt8wDrxno0LdLib5E2tOYG+bGUG9g/ujn8mvZiYhIhRu+Lc2IzSmGN6cKg/LaeuxRHIh35Uy46hD44bNPYOZ878mHC1sLDUR3yLSH7xMamKdEX8Lx4niLAA5rnKnBeCIi+4j4ILx8oh/hUap0I0wQ4fa3/SJc8uwt180E4Jozzisk1ULeMINUeA/g/HnlX6IWlOh7aJrTjrkVv4oZTOkwOk85ebvHiYjI3m/rfjk8+g+gNZb0ewvLu6yL6oXVgNHdXezmzy9qw7ZtjHr1fgDaBzD27dGo+xJNZ8YuRRaCR56kevLk5AISEZGyGA60jI51vGfAW2KjtEpy/s/+YRS1z9XCyk5yw2HTfkcRQGGStRRwwCfOZ+33HwDg/nMWD2wF+qEWfaluClPgkoIX6IR16xIOSkREBtqFkw/skQNGbktjFk6Fe/35J1A6Mvtdd28B4N5lR5EC7vpzd2Fb/jtBPTC08zRNgTtYNM5tCAdLhg9BQjOctWQ6XHBBwpGJiMhA27p1KyNb04zckmZkc7ow6527EwQBhUxfHGmHBwYsx3MA6wvd5PldXj7v3aRay1SBfijRl5g2rx5WA3PhqKfquP/yBZDRuHsRkX1FfJ77fLIPl6G1aLk6L2TyYkM9xxXjrigctw5YPGkWh//+pwTDyxF1/5To+9C4uIFpVo+nnLpFMOSrX0s6JBERGWAjR44EiqvP5bvtAbq7u6MB9xa9wj2zWbjieEibYZzLEVd9jBwwDpiy/CGW3jKPVk5MojoFGowX05TuZlounLnIMNKkuTB3Jq+s1Xr0IiL7ivwjdu5eSPY1NTXhMPqSW+1V+cnuUil83OkALL31wbDMoh6AuhuinoBkqEUfNyvL83SQiabIOfLAGTy75DG23apELyKyr+k1fK7nmjYc2ZwjC+S2bOl1bHydu5x7OFlOQreBlehjGuc2wCxYPDtD0+x27hj2CE37tUNb2w6PFRGRytKrDV5SMJQ06dWf4WcrngEodPN77D3A5s2bw8VsEhp5r677uBw0Phhbva6bcNnaw+uTikhERBJQWJM+3uUeZCFdVUj4swiwA7/KJQd4j/1K16hf89RDjD3nIuhshepRAx57KbXoYya+fAy0EF7EDBy8spH91hwDKf2ZRET2BX9/4lgs/5/1HJBHlUE2CO/VN8PXR24jCAJueOiPPfYrTIMbfT628WTo6GA455a/QijR9/DK5IUcvL6RCcuOZsIrR9NOO1vPryK9YEHSoYmISBlUpcO06L077jmI+VBtkHYYGfBA+iN84+E/9XmeHrfz99+f4T4Hlve+l18O6rovserwpuKHLIz6wLfC91demUxAIiJSVhlyVJGKHpcvJvx7s//LkCHF/To74ZrZ5xX3y+8aBJBOF0faZ/6d2hWfxNMkQi36vkSrEYxdcTgONJ91ZtIRiYhImXTuBwFBYVBdvgt+7syfkc2G+TuXgyPtZ/nH6aOu/qgd744HQbhjdzcBbTiwkWRmzlGLPqZpTjuNDzQUViHcMGYp6Vtv0TK1IiL7kH84eDR3+WZGbOpZ/onHnTecX+yq/15LC6yAaJx9cRDea2vgoAlghtfU4MH1GHcxmhY2l6kOcWrRx60Ikz1bYcWWdpouTGpRQRERScrIujpqN+fnu4mPqN/IU7e1s3XrVr5422ZeuivfF289R9qvXdmjgZhOpdgAbD7qwvJUoIQSfUzjsgbIQdOF7bReBGwm0RWHREQkGTWkcO/ZdX/t1glkyLDs7hreQi2/prOwf3y/I779b2HfvnuY8LduhmkXJVCLkLruS0x7uL4wzzEAR3miUxeKiEj5Daky2ruD06D+zQAAEZtJREFUHmUdW67nW5M+TDjBSk+lz86fu63YSZ8D7h+IIHeSWvQlzIuDKcgBK1Yw/IOXJhqTiIiU15tn7F9ooed/5nbh+D8Rfh3oAO4f9YY9Hd4uUYu+lBf+D1z/LYbTxzSIIiJS8VpbWxk6dGi4uI1nOe2qz3Hj3fERehncwS0o9AR3d3eHOWPkWB6I3/r1qHc4gdvBatGXyj8rgcE1H+XPd31QiV5EZB8UjKuheb+A5v0Cto5J8cCvfks4LV5eNZNWjODJaV2FtFFTW8Ooqy8rJnaP3f5NaN0UJfoSjhemNMqRg+qel1VERPYN7512MKmOTNhnv2Qtiw44m9Il7CbRDPZ3UUlY9gLfLp4kv349EMRn2ykjJfpS+Vv05pAOP+uPJCKy70mljCcfXsXEFzfjm435m9spNv0ccL6ZApoWRiXhf0ye0Ocg7lRHR7lC70H36GMyBlWxKQ8DAlive/QiIvuq369p5vdrmmMlWcJWIIDx98FIXg4205ZK0XTMcQAM2TiBluogXBAtlvDrs2fRwWNliz1PjdWYxbPbyZABwiS/eFYnx35xnhK9iMg+avF/nNmjtffC/uspjr/PsvyQrdzw0B+pye+QzTL82k8XE3ys6/6Jx8sUdAm16GNO+8ef88iv3wFRsqelleHLV9CSaFQiIpKU2traHon+qHWwfNKawud8Pr/7MZh2wwfDMoCvXgdf+EpxB3f+7jSo/9MUOCneQzDw1KKPS6WY+eZbYds2+MavmPmOX4flEyYkG5eIiCTm2pZ6aAZaYUbzBL79NLS2QksLfOeZKdSNOo/1wXhaga3AQ7O+yfAZD0L7KdDRCu1bGd7xJjrdmZVA/Im16M3sY8CVhH0gv3f3T0flnwUujco/7u73li2o1lYATvunX0XDLKKxef9xXdlCEBGRweXsjx0G/7cUctFobZvKj5ZGbw1gCgzJ8rdZ3wKcEbwTjgsYzkjgnDCROAzv7Ozz/AMtkRa9mZ0BnA8c6+7TgP+Oyo8C3glMA84FvmNmCa3gi1atExERDj9qAtBJsQ/f6T1MO99uNpr5BYwcGe1afI5+9t/2refoPwJ8zd27ANx9fVR+PvBzd+9y9+XAUuCNZYuqurrwtkeS/8AHyxaCiIgMLuEUuDubLqPssXITqT9vgmwWstnwfUKSSvRTgJlm9lcze9jMTozKDwJWxfZ7NSrrxcwuN7P5ZjZ/w4YNeyaqTKbP4mSmOBARkcHigo+eQKEPfiec8QqcAsx+vJnZjzdzykAGtwMDdo/ezB4AxvWx6XPR7x0NnAScCPzSzA7dlfO7+03ATQAzZszYo0/AlZ4smSkORERksDioMb8wTb6/t7T7/mng+ELZgzP3gyB6lj4yO6FW/YC16N39THdv7ON1N2FL/U4PPUk41dAYYDVwcOw046OysopfRlCLXkREgIbSgny2CPjIuC8RjiF34Clobw9H6kXr0vPiJtrLGWtMUl33vwHOADCzKUANsBG4B3inmdWa2SRgMvBkQjHihN9AUgmsNiQiIoPL5V87namnjKb3evTdbM3CFePezhXjLuQj474CT3XAps2QTocJ/4jRiTUak3q87hbgFjNrArqB97m7A8+b2S+BFwj/kle6+64sAbx7hg+HlpZCSz5FOHXOtOcWli0EEREZnKpr0qx77EnCjuc0xX7fFobkKGRUA2YBqRcBkhuEl5dIonf3buCSfrZ9BfhKeSOKtIRz4OWHWzjRujZVmkBQRETgkusvovOLk/jB1jsIM0SWK/b/MPn1auLD9fLvk+4TVgaLa2gI76tEHE0dKCIiMQ2j+fHWi7li/7fnp7DHPUwf8eF5W4H9EgqxlPJYXEp/DhER2b7LvvzlHp/NojlxHNqAeYSPleUlvTCaMltca2thzKQDrTfflHBAIiIy6IwZQ+kY7fwidUOBORS7643ku+6V6EukKV6YIZddnnA0IiIyKA0ZX2gUJt1i3xEl+rjYFLgQJn0REZFern2eVWMv6fWg3WCkRB8XmwLXGBxdLiIiMjiNG1LNGg5POowdUqLvg5e8REREStW8/xtMZOmgbxAq0fdBLXkREdkpVy9KOoId0nP0JUpb8H2vZyciIgKMOCDpCHZILfq4sWN7tOQdrVwnIiI7cF1z0hFslxJ9XLSufY/V637w/aSiERER2W1K9HEN4RqE+e77VoDNm5OKRkREZLcp0cdFU+C2Ai3HHM0wYPi1n040JBERkd2hwXhxra0ADANY+JwerRMRkb2eWvT9UJIXEZFKoEQfVzIFroiIyN5OiT4u0/Op+QBNnCMiIns3Jfo+5Lvtc6gLX0REdkZd0gH0S4m+xDbC5J5Bk+WIiMhO+lRT0hH0S6PuSwwj7K6vAXTHXkREdsrQsUlH0C+16OOOmAr0XLVOXfciIrJTBulUuEr0ce1hZ318CtxtJ5+UWDgiIrKX+dTy7W9PN5Qnjhgl+pgj77sXKLbiswCNjUmFIyIie5uho6OWfT/p9dKnyhoOKNH3ctSiF2kBWqZOpWPSIYy78MKkQxIRkb3NdVvgrd/pWXb5i3DggWUPRYPx+jB94bOsWrWKVCrFwQcfnHQ4IiKyN5rx7vCVMCX6PtTU1HDYYYclHYaIiMhuU9e9iIhIBVOiFxERqWBK9CIiIhVMiV5ERKSCKdGLiIhUMCV6ERGRCqZELyIiUsGU6EVERCqYEr2IiEgFU6IXERGpYEr0IiIiFUyJXkREpIIp0YuIiFQwJXoREZEKpkQvIiJSwZToRUREKpgSvYiISAUzd086ht1mZhuAV/bwaccAG/fwOQcT1W/vV+l1VP32fpVex6TrN9Hdx+5op4pI9APBzOa7+4yk4xgoqt/er9LrqPrt/Sq9jntL/dR1LyIiUsGU6EVERCqYEn3/bko6gAGm+u39Kr2Oqt/er9LruFfUT/foRUREKpha9CIiIhVMib6EmZ1rZovNbKmZfSbpePYEM1thZs+Z2TNmNj8qG21m95vZkujnqKTj3BVmdouZrTezplhZn3Wy0Deja7rQzKYnF/nO6ad+15nZ6ug6PmNmb4lt+2xUv8Vmdk4yUe88MzvYzB40sxfM7Hkz+0RUXknXsL86VsR1NLM6M3vSzJ6N6vcfUfkkM/trVI9fmFlNVF4bfV4abT8kyfh3xnbq+EMzWx67hsdF5YPz36m76xW9gDSwDDgUqAGeBY5KOq49UK8VwJiSsv8CPhO9/wzwn0nHuYt1Og2YDjTtqE7AW4A/AgacBPw16fhfZ/2uAz7Vx75HRf9Wa4FJ0b/hdNJ12EH9DgCmR++HAS9F9aika9hfHSviOkbXYmj0vhr4a3Rtfgm8Myr/LvCR6P0VwHej9+8EfpF0HXajjj8ELupj/0H571Qt+p7eCCx195fdvRv4OXB+wjENlPOB26L3twEXJBjLLnP3R4DNJcX91el84EceegIYaWYHlCfS16ef+vXnfODn7t7l7suBpYT/lgctd1/r7gui99uAF4GDqKxr2F8d+7NXXcfoWrRGH6ujlwOzgV9F5aXXMH9tfwXMMTMrU7ivy3bq2J9B+e9Uib6ng4BVsc+vsv3/Ye4tHLjPzJ4ys8ujsv3dfW30fh2wfzKh7VH91amSrutHoy7BW2K3W/bq+kVduMcTtpYq8hqW1BEq5DqaWdrMngHWA/cT9kJsdfdstEu8DoX6Rdubgf3KG/GuK62ju+ev4Veia3iDmdVGZYPyGirR7xtOdffpwJuBK83stPhGD/ucKurxi0qsE3AjcBhwHLAW+J9kw9l9ZjYU+DVwlbu3xLdVyjXso44Vcx3dPefuxwHjCXsfjkg4pD2utI5m1gh8lrCuJwKjgX9JMMQdUqLvaTVwcOzz+Khsr+buq6Of64G7CP8H+Vq+Syn6uT65CPeY/upUEdfV3V+L/p9OANxMsVt3r6yfmVUTJsCfuvudUXFFXcO+6lhp1xHA3bcCDwInE3ZXV0Wb4nUo1C/aPgLYVOZQX7dYHc+Nbsu4u3cBtzLIr6ESfU9/AyZHo0ZrCAeM3JNwTLvFzIaY2bD8e+BsoImwXu+LdnsfcHcyEe5R/dXpHuC90YjYk4DmWPfwXqPkXt/bCK8jhPV7ZzSqeRIwGXiy3PHtiuje7A+AF939+timirmG/dWxUq6jmY01s5HR+3rgLMJxCA8CF0W7lV7D/LW9CJgX9doMWv3UcVHsy6gRjkGIX8PB9+806dGAg+1FOGryJcJ7TZ9LOp49UJ9DCUfyPgs8n68T4b2xucAS4AFgdNKx7mK9bifs9swQ3ge7tL86EY6A/XZ0TZ8DZiQd/+us34+j+BcS/j+UA2L7fy6q32LgzUnHvxP1O5WwW34h8Ez0ekuFXcP+6lgR1xE4Bng6qkcT8IWo/FDCLyhLgTuA2qi8Lvq8NNp+aNJ12I06zouuYRPwE4oj8wflv1PNjCciIlLB1HUvIiJSwZToRUREKpgSvYiISAVTohcREalgSvQiIiIVTIlepEzMLBetdNVkZr+NPZ97oJn9aieOb+2n/AIzO2oHxz5jZj9/fZHvGTtbz508V35Fxhl9bJtlZr/bjXM/aGatfZ1bZG+kRC9SPh3ufpy7NxIuWHMlgLuvcfeLtn/odl1AuPJZn8zsSMKVGWdGkyYlYg/Us9QZ7j5/D54PAHc/A9jj5xVJihK9SDIeJ1rswswOsWjdeTNrMLNfWriG+V3Rut2FlqWZfSVaG/sJM9vfzE4B/gH4etRqP6yP33Ux4SQt9xFbjdHMPh79noX51r6ZDTWzW6PW8kIzuzAqP9vMHjezBWZ2RzR/e75l/R9R+XNmdkRUfroV1+p+2syGldSzLvZ7njazM6Ly95vZnWb2JwvXpP+vnfljmtm5ZrbIzBYA/xgrH2LhwjFPRr/n/J35O4tUEiV6kTIzszQwh76nV74C2OLuRwH/BpwQ2zYEeMLdjwUeAS5z98ei81wb9RYs6+Oc7yBccvl2wqSf9xngeHc/BvhwVPZvhNN2Hh2VzzOzMcDngTM9XBxpPnBN7Dwbo/IbgU9FZZ8CrvRwMZCZQEdJTFcSrltzdBTTbWZWF207Lor5aOAdZnYw2xEddzPw99Hfa1xs8+cIp1p9I3AG4ReiIWz/7yxSUZToRcqn3sLlLvPLr97fxz6nEiZl3L2JcOrNvG4gf+/5KeCQHf3CqJW60d1XEk4te7yZjY42LwR+amaXAPllRc8knMKTKIYtwEmEtwYejeJ/HzAx9mvyC9LEY3oUuN7MPg6M9OKypfF6/iT6HYuAV4Ap0ba57t7s7p3ACyW/qy9HAMvdfYmHU33+JLbtbOAzUdwPEU7DOoHt/51FKooSvUj5dEQt3ImEc2JfuYvHZ7w4Z3UOqNrezpGLgSPMbAXh/NvDgQujbecRJvXpwN+suOJYKSNch/u46HWUu18a295VGpO7fw34Z6Ce8AvCrixf2hV7v7P17I8BF8Zin+DuL+7G+UT2Okr0ImXm7u3Ax4FP9pFcHwX+CSAaSX/0TpxyGzCstNDMUtG5jnb3Q9z9EMJ79BdH2w529wcJ19IeAQwl7GW4MnaOUcATwJvM7PCobIiZTWE7zOwwd3/O3f+TcFXI0kT/Z+Dd0b5TCFvZi3eirn1ZBBwSG58Qvz1xL/CxaJUxzOz4qPz1/J1F9kpK9CIJcPf8ilgXl2z6DjDWzF4Avky44mDzDk73c+DaaLBZfDDeTGC1u6+JlT1C2A1/EPATM3uOcHWub3q43vaXgVHRI4DPEo5s3wC8H7jdzBYSDiTcUQv9qugcCwlX4PtjH/VMRb//F8D7PVzbe5dFXfyXA7+PBuOtj23+ElANLDSz56PP+d+/q39nkb2SVq8TGUSigXrV7t4ZJe0HgKnu3p1waINKdCtihrtvfJ3Hb/fvbGYPAZ8aiMf3RMptd+59icie1wA8aGbVhPeXr1CS79MGYK6ZXfo6k3G/f2cze5BwTfXMHotWJEFq0YuIiFQw3aMXERGpYEr0IiIiFUyJXkREpIIp0YuIiFQwJXoREZEKpkQvIiJSwf4/ndDjMBVOMDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4446338048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Fig, Ax = plt.subplots(1,1,figsize=(8,8))\n",
    "Ax.scatter(dm_hp_ra, dm_hp_dec, c=km.labels_, cmap=plt.cm.tab20, s=6)\n",
    "\n",
    "Ax.set_xlabel('Right Ascension [deg]')\n",
    "Ax.set_ylabel('Declination [deg]')\n",
    "Ax.set_title('{0}'.format(FIELD))\n",
    "Fig.savefig('plots/dmu24_{0}_sf_hpclusters_illustration.png'.format(FIELD), \n",
    "            format='png', bbox_inches='tight', dpi=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II - Map photo-z and masterlist objects to their corresponding depth cluster\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now load the photometric redshift catalog and keep only the key columns for this selection function.\n",
    "Note: if using a different photo-$z$ measure than the HELP standard `z1_median`, the relevant columns should be retained instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "photoz_catalogue = Table.read(PHOTOZS)\n",
    "\n",
    "photoz_catalogue.keep_columns(['help_id', 'RA', 'DEC', 'id', 'z1_median', 'z1_min', 'z1_max', 'z1_area'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we load the relevant sections of the masterlist catalog (including the magnitude columns) and map the Healpix values to their corresponding cluster. For each of the masterlist/photo-$z$ sources and their corresponding healpix we find the respective cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "masterlist_hdu = fits.open(MASTERLIST, memmap=True)\n",
    "masterlist = masterlist_hdu[1]\n",
    "\n",
    "masterlist_catalogue = Table()\n",
    "masterlist_catalogue['help_id'] = masterlist.data['help_id']\n",
    "masterlist_catalogue['RA'] = masterlist.data['ra']\n",
    "masterlist_catalogue['DEC'] = masterlist.data['dec']\n",
    "\n",
    "for column in masterlist.columns.names:\n",
    "    if (column.startswith('m_') or column.startswith('merr_')):\n",
    "        masterlist_catalogue[column] = masterlist.data[column]\n",
    "\n",
    "masterlist_hpx = coords_to_hpidx(masterlist_catalogue['RA'], masterlist_catalogue['DEC'], ORDER)\n",
    "\n",
    "keys = np.array([*clusters])\n",
    "incl = np.array([hpx in clusters.keys() for hpx in masterlist_hpx])\n",
    "if incl.sum() != len(incl):\n",
    "    notincl = np.where(np.invert(incl))[0]\n",
    "    for i in notincl:\n",
    "        masterlist_hpx[i] = keys[np.argmin(np.abs(masterlist_hpx[i] -keys))]\n",
    "    \n",
    "masterlist_catalogue[\"hp_idx_O_{:d}\".format(ORDER)] = masterlist_hpx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "masterlist_cl_no = np.array([clusters[hpx] for hpx in masterlist_hpx])\n",
    "masterlist_catalogue['hp_depth_cluster'] = masterlist_cl_no"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = join(masterlist_catalogue, photoz_catalogue, join_type='left', keys=['help_id', 'RA', 'DEC'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constructing the output selection function table:\n",
    "\n",
    "The photo-$z$ selection function will be saved in a table that mirrors the format of the input optical depth maps, with matching length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "pz_depth_map = Table()\n",
    "pz_depth_map.add_column(depth_map['hp_idx_O_13'])\n",
    "pz_depth_map.add_column(depth_map['hp_idx_O_10'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III - Creating the binary photo-z selection function\n",
    "\n",
    "With the sources now easily grouped into regions of similar photometric properties, we can calculate the photo-$z$ selection function within each cluster of pixels. To begin with we want to create the most basic set of photo-$z$ selection functions - a map of the fraction of sources in the masterlist in a given region that have a photo-$z$ estimate. We will then create more informative selection function maps that make use of the added information from clustering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NCLUSTERS # Fixed during the clustering stage above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 In cluster: 596592 With photo-z: 425630    Fraction: 0.713 \n",
      "1 In cluster: 211678 With photo-z: 189034    Fraction: 0.893 \n",
      "2 In cluster: 667672 With photo-z: 487879    Fraction: 0.731 \n",
      "3 In cluster: 890995 With photo-z: 642956    Fraction: 0.722 \n",
      "4 In cluster: 562339 With photo-z: 522927    Fraction: 0.930 \n",
      "5 In cluster: 197157 With photo-z: 188706    Fraction: 0.957 \n",
      "6 In cluster: 64899  With photo-z: 40065     Fraction: 0.617 \n",
      "7 In cluster: 811444 With photo-z: 527891    Fraction: 0.651 \n",
      "8 In cluster: 426919 With photo-z: 296813    Fraction: 0.695 \n",
      "9 In cluster: 665740 With photo-z: 441831    Fraction: 0.664 \n",
      "10 In cluster: 83455  With photo-z: 59243     Fraction: 0.710 \n",
      "11 In cluster: 68772  With photo-z: 46715     Fraction: 0.679 \n",
      "12 In cluster: 76238  With photo-z: 58039     Fraction: 0.761 \n",
      "13 In cluster: 839498 With photo-z: 626420    Fraction: 0.746 \n",
      "14 In cluster: 77293  With photo-z: 57192     Fraction: 0.740 \n",
      "15 In cluster: 107165 With photo-z: 94739     Fraction: 0.884 \n",
      "16 In cluster: 577869 With photo-z: 409561    Fraction: 0.709 \n",
      "17 In cluster: 241038 With photo-z: 194667    Fraction: 0.808 \n",
      "18 In cluster: 57413  With photo-z: 49282     Fraction: 0.858 \n",
      "19 In cluster: 138182 With photo-z: 94066     Fraction: 0.681 \n",
      "20 In cluster: 154254 With photo-z: 105926    Fraction: 0.687 \n",
      "21 In cluster: 96444  With photo-z: 69522     Fraction: 0.721 \n",
      "22 In cluster: 83242  With photo-z: 59300     Fraction: 0.712 \n",
      "23 In cluster: 546857 With photo-z: 401100    Fraction: 0.733 \n",
      "24 In cluster: 113801 With photo-z: 71302     Fraction: 0.627 \n",
      "25 In cluster: 215231 With photo-z: 132240    Fraction: 0.614 \n",
      "26 In cluster: 99679  With photo-z: 73718     Fraction: 0.740 \n",
      "27 In cluster: 169650 With photo-z: 122592    Fraction: 0.723 \n",
      "28 In cluster: 84213  With photo-z: 58959     Fraction: 0.700 \n",
      "29 In cluster: 68793  With photo-z: 53956     Fraction: 0.784 \n",
      "30 In cluster: 43780  With photo-z: 30373     Fraction: 0.694 \n",
      "31 In cluster: 513089 With photo-z: 392497    Fraction: 0.765 \n",
      "32 In cluster: 85659  With photo-z: 58855     Fraction: 0.687 \n",
      "33 In cluster: 159564 With photo-z: 109789    Fraction: 0.688 \n",
      "34 In cluster: 86970  With photo-z: 81054     Fraction: 0.932 \n",
      "35 In cluster: 129861 With photo-z: 94664     Fraction: 0.729 \n",
      "36 In cluster: 52994  With photo-z: 38451     Fraction: 0.726 \n",
      "37 In cluster: 70968  With photo-z: 50992     Fraction: 0.719 \n",
      "38 In cluster: 462542 With photo-z: 334900    Fraction: 0.724 \n",
      "39 In cluster: 474188 With photo-z: 314212    Fraction: 0.663 \n",
      "40 In cluster: 571123 With photo-z: 414266    Fraction: 0.725 \n",
      "41 In cluster: 51055  With photo-z: 45392     Fraction: 0.889 \n",
      "42 In cluster: 866730 With photo-z: 615039    Fraction: 0.710 \n",
      "43 In cluster: 98858  With photo-z: 67972     Fraction: 0.688 \n"
     ]
    }
   ],
   "source": [
    "cluster_photoz_fraction = np.ones(NCLUSTERS)\n",
    "\n",
    "pz_frac_cat = np.zeros(len(merged))\n",
    "pz_frac_map = np.zeros(len(dm_hp_ra))\n",
    "\n",
    "for ic, cluster in enumerate(np.arange(NCLUSTERS)):\n",
    "    ml_sources = (merged['hp_depth_cluster'] == cluster)\n",
    "    has_photoz = (merged['z1_median'] > -90.)\n",
    "    \n",
    "    in_ml = np.float(ml_sources.sum())\n",
    "    withz = np.float((ml_sources*has_photoz).sum())\n",
    "    \n",
    "    if in_ml > 0:\n",
    "        frac = withz / in_ml \n",
    "    else:\n",
    "        frac = 0.\n",
    "    \n",
    "    cluster_photoz_fraction[ic] = frac\n",
    "    print(\"\"\"{0} In cluster: {1:<6.0f} With photo-z: {2:<6.0f}\\\n",
    "    Fraction: {3:<6.3f}\"\"\".format(cluster, in_ml, withz, frac))\n",
    "    \n",
    "    # Map fraction to catalog positions for reference\n",
    "    where_cat = (merged['hp_depth_cluster'] == cluster)\n",
    "    pz_frac_cat[where_cat] = frac\n",
    "    \n",
    "    # Map fraction back to depth map healpix \n",
    "    where_map = (km.labels_ == cluster)\n",
    "    pz_frac_map[where_map] = frac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The binary photo-$z$ selection function of the field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAHwCAYAAABe2J4CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XmcZHdd7//Xp7bel5nuWZLJJBOykLAIYhDBoAFxARdEQMAVN9xAELy/q957NSjee+UCgv68ahBk8V4wqPwMiMID2RRlSYBAQlhCkpnMTDLTM9N7134+vz++51Sd7umerpnp7upT/X7mUY+uc+qc6m9NV+p86vP9fj9fc3dEREREsi7X7QaIiIiIbAQFNSIiItITFNSIiIhIT1BQIyIiIj1BQY2IiIj0BAU1IiIi0hMU1IiIiEhPUFAj0sPM7EYz+3czmzWzM2b2STN7opmVzOz1ZnbUzBbM7AEze2PqvAfMrGxm82Y2Ez/HL5lZLnXM28ysFp+f3F7QnVcqIgKFbjdARDaHmY0C7wd+GbgVKAFPBarAbwE3AN8KPARcAXzHiqf4QXf/sJmNAd8JvAl4EvAzqWNe6+7/dTNfh4hIpxTUiPSuawHc/V3xdhn4EICZ/Xfgve5+PH7sgfh2FnefBW4zs4eBT5nZ6939rk1st4jIBVH3k0jv+hrQNLO3m9kzzWxX6rFPAa80s18xs8eama33ZO7+GeAoIdsjIrLtKKgR6VHuPgfcCDjwZmDKzG4zs33A/wD+EPhx4HbgmJn9dAdPexzYndr+jXjMzYyZndrYVyAicn4U1Ij0MHe/x91f7O6XAY8BLgXe6O5Nd/9Td/92YBz4A+CtZnb9Ok95ADiT2n6du4/Ht8lNeREiIh1SUCOyQ7j7V4C3EYKb9P6yu/8pMA08aq3zzeyJhKDm3zaxmSIiF0xBjUiPMrPrzOxVZnZZvH0QeBFhsO8rzOwmMxsws0Lc9TQCfH6V5xk1sx8A3g38tbt/aStfh4hIpzT7SaR3zROmYL/SzMaBGcIU7/8EvBB4PXA1YczN14Dnuvt9qfPfZ2YNIAK+DLwB+POta76IyPkxd+92G0REREQumrqfREREpCcoqBEREZELZmZvNbOTZrZqUU4L/tjM7jWzL5rZEzarLQpqRERE5GK8Dfi+czz+TOCa+PYS4M82qyEKakREROSCufsnWF6/aqVnA+/w4FPAuJldshltUVAjIiIim+kA8GBq+2i8b8P1xJTuyclJP3ToULebISIisswdd9xxyt33bMXv+t6nDfnpM80Nf947vli9G6ikdt3i7rds+C/aAD0R1Bw6dIjbb7+9280QERFZxswOb9XvOn2myWc+ePmGP2/+kq9X3P2Gi3iKY8DB1PZl8b4Np+4nERGRHuBAtAn/bYDbgJ+KZ0F9GzDr7g9txBOv1BOZGhEREekOM3sXcBMwaWZHgd8FigDu/ufAB4BnAfcCS8DPbFZbFNSIiIj0BKfpG5JZOb/f6v6idR534Fe3oi3qfhIREZGeoEyNiIhIDwhjanb2eo4KakRERHrEBg3szSx1P4mIiEhPUKZGRESkBzhO03d295MyNSIiItITlKkRERHpERooLCIiIpnnQHOHBzXqfhIREZGeoEyNiIhIj9jp3U/K1IiIiEhPUKZGRESkBzjs+CndCmpERER6xM6uJ6ygRkREeshb7v1NHqzc09q++TH/0MXWyFZTUCMiIj0jHdDsNI5rSne3GyAiIrJZ/vq+13S7CbKFlKkREZEe5dy79B/dbsTWcWju7ESNMjUiItJrdviVfQdTpkZERHqY0fSIvPX+d3hHs596/68sIiI7jKXuP8DhhXu71pKtZTQ34ZYlCmpERKQnlBsLRBGEnEXSBXWITz/8ou41SraUup9ERKQnTNdOYgbtTI1zw8ADRBnLNlwoB6IdPpxIQY2IiPSEW+779fiekwQ21Src19zTtTbJ1lL3k4iI9IhpzNoBDQ5fbF7FvI/wlnt+u6st2yo7fUyNMjUiItIj5oDx9qYZFmdtPj/3+W41ass4ZC4I2WjK1IiISE+4Pt+I73nqpxNFUCotdKlVspUU1IiISE+4p3klZxfecyKOMFoc4fDCV7vRrC0VuW34LUsU1IiISI+4H87qfpnm24bqfEv/Ef7qgZ/rRqNkC2lMjYiI9ITHFeDORhPIk2Rsnjg423r8iYMz3WnYFtGYGgU1IiLSI+5sXAU04632xd2SyVA9XsPFMZo7vANmZ796ERHpMfnU/TBI2D3cogjeeM93d61lsvmUqRERkZ51R+URXMpxKjQ5w0FGc0vdbtKmytrA3o2mTI2IiPQYJ73+03Eu5QwHgTKnKkMcmTnSzcbJJlJQIyIiPSa5tLWzFu6G+yCLS/Cyu/6wO83aZMlAYVUUFhER6SlGumaNWRhXMzgAlGG+Os1I366utW5zGE3f2bmKnf3qRUSkx6SnOBkr69ZEEQzny/zJ139pS1slW0OZGhER6TGrBzZmMDDgXMMp9hTPdKVlm8mBaIfnKnb2qxcRkR7iwHx8P52lWRnkwERhgQ/e966ta5psCQU1IiLSE6LIuGFgiuWzn5KMzPIZUSfqo/zH3B90o5mbaqcPFFZQIyIiPeOzS1fx2Nz9QJX5eeNxhWngCFCPbzXAOV7dBQW4+a5f7WZzZYNpTI2IiPSEXC7CPceXoqtwh5GROnc2HsENA/dh9iDuIegxi6g5DABRdLTbzd4w7pr9tLNfvYiI9JAIs9DFFH4u4W58dukqPrN4FZ9dugqARiPHQg3w/SzU4ffuen43G72hImzDb1mioEZERHrCt/QfJixo6UCNKznVWszSrL2w5ZfmIZcb4DG5hylG0IhqnFia61KrZSMpqBERkZ5w5gw8cfAwTxy8nycOHuVvp69adWXucvkyHq7BXRHcFvc+ff9Hbt7Stm6GUFE4t+G3LNGYGhER6Ql/W/12OFZkko9xipuARdzPtKoJJx5cnGBsuIZzFX2lA8AxGh51qdWykboSgpnZzWZ2zMy+EN+elXrst8zsXjP7qpl9bzfaJyIi2TNb6weiOKCJmKydYbYBtWb7mEYTYIYrcieJIgDnc1NQbQx0ocUbLQwU3uhblnQzU/NH7v669A4zexTwQuDRwKXAh83sWndvrvYEIiIibUeYrV3V2pplH1dX5sgNOpXUVeTqyTqHo/1cWX0YGOJo5VJgactbu9FUUXj7jal5NvBud6+6+/3AvcC3drlNIiKSAQeG+oCThMHCc1w66Mwzhru1igq7QyEXRgx/aOqxXD2+yKFh2NPfrVbLRupmpualZvZTwO3Aq9x9GjgAfCp1zNF4n4iIyDkN1uDSwX7MFnEPdVvAqdehWAyBTaMB0CBkZsZYXISH6lCtdbXpG6bp2ZqCvdE2Lagxsw8D+1d56L8Afwb8PiF2/n3g9cDPnufzvwR4CcDll19+UW0VEZHsm7YxSnhrULA7YBHH67ug3r7Y9+Ua1JsjXDU2jxkcasB0uTttlo21aUGNuz+jk+PM7M3A++PNY8DB1MOXxftWe/5bgFsAbrjhhlUm7YmIyE7SV2ywtJSnvz8ENCcrfezpr0EEbmAWMjfHl2CuAlePh/MeWNjNUDH7UY1jmZuCvdG6NfvpktTmc4C74vu3AS80sz4zuxK4BvjMVrdPRESyaJFZH+ZEeYSTlREgB1aNr3SGe1j0EkaAXTSbzbg7qkql0eSOY/d0se0bI/Lcht+ypFtjal5rZo8ndD89APwigLvfbWa3Al8mdHr+qmY+iYhIJ/b1gfsCD8cTmfYMQFjhCcyc+YpRw5noXySaW+Irp6GvVKJcHQHgxz59KzfuvYS33PiS7rwAuWhdCWrc/SfP8dgfAF1ZD/4HPvIqFomAiN84+Pc8etetHJp8cjeaIiIi52mxDLU67B8M23NLMJSUn3EYiWc4RZGxfy/kKgPUm1CuAoSuqc+cPsL17/1v/MKhb+aV3/wjW/8iLkJSUXgn29mvPuWOU19jESe8sfO87sEf5q7pn+92s0REpEOD/XCqMcLxpVGOL42ywEh70HDc/eRueKtLpcyBgWQsTTiwlK8Dc/zlA7dvdfNlA2iZhNjHHvoC7tYqp21eYCyf/YFjIiI7xfHKGPsGFjlRHiB8QV0Kg4PjwCa9XMLJCuzrh8MzMD4wx0J1kL58naGBBkt1cPLdehkXzDFN6e52A7aLv3/4Doq59hu+5jnua+7iqd1tloiIdCiKnFzO2DdQJhkYvJZ9AxG1CA6O55mqQF8p+RIb1oAaz2vV7ixSUJNSj9KReQModaspIiJyntxzccG9kJUJ2Zr04+2ftagARPER6eNyXLE7u4tbapkEAWBqbpgQyABEPG/0U+c6XEREtpk6uVWzM3mg0QgPRBHMNsKEkMTCQj2+5+CNOCDKHnd2/IKW2WrtJpuaG2dqbpipuVH+7OhTaXXEiohIBiwSLfvYblcXLhbDtO583pnoa5K+/JWj5LN/hKn5sXN2W8n2pu6nZRYJXU41Qm2DjIbrIiI7UAnDSQKWdjeUO0Ru5KjTdDArAhHF+CP+it3THD4zGN8vLxtQnC3W6lDbqRTULDNI+B+hyPX8M+n0pIiIbG9NmhRXTFpyh2Yc50QUW11LxdS1f6lhXLE7DBSOIshnb+KTxBTULNN+l9/D04F/6V5TRETkvBTzA6mtkGqp1WAmGmJ3X7tEh8clydyTAcV5lhrt7qjBXJ0scsjcGJiNtrNffcqBXTNAmfC2aHBg10KXWyQiIuer3lpYJ9QdK8Rf3ZuedEO1c/BJ1mZqbgRoAk6zXOfwmS1t8oZqktvwW5YoUxOrNQsc2FUFqvGeh7rZHBEROW81ivmkFIfHFYSTxyKiVkG9ENbU602KxTwwz9TcWOtZrtg9ndExNaKgJjEH1ZEiffavDDlMNZ8KHOl2q0REpEMv3n87b3s43L8UuLPyFMYLeSb6F8Gh4TUACrkQ3MzPwpkCXLarQt6aHD4DV+wOgdCR6WKXXsWFc4xIFYUFYIphdjUXqOafSsVhemkAmO12s0REpEMzlRIv3h8Cl0/PXkpUHqFQWArdTiTBTFi4EmBiD0wAM7WwfWiindk5NJHNcTU7nYKalhmml8ZS20vAeLcaIyIi5+kzlWuxeARBCE5qrYHCE8VF8nmPx9GEyGWhAUOpmU5JQBMB9SqZlLUxMBtNQU3swC44Nn0EOAAscclYAxXfExHJDjNoNEJwcrIBe0tVmp6HaJ5CwXGPIxhPPt3DdiEHjVQFj7zB0YXRrW7+RXMg2uGznxTUxGrNApeMjWEWZj3Vo539xhARyZpaDUrxOOF9Bl8sw76BJpcWw+yncNk3PC7fMZxv4uQpReDmYUVvkoyNFrTMIgU1iTo0ivlUckb9qSIiWVJMje0tFGBPESAXBzRJKqbdBRVFIbtTLEK1sXxV7yt2ZzFTbzR3eEVhpSNiU+VhqNaIyxfx9Manu90kERE5D8croSKwO5wpQy5PnHaJqFZJrcht4AbxsQ/PhPPN0redHRxklTI1LTNMVXe3ytT8Dd/Gf+Hfu9skERHpnI/yULX9XX26ErFrwMkb9BWTzicAB4OpOrgPhyX/qLUqDGe1Ro3G1ChT03JgF4Qp3PPADPtG4mheREQyIXxiR0BE3iJqQKUCx5fgZDXfOqJdlG+I+WqV/YPtCvJZDWgkUKYmdmx6nP2jM+RyEe7GQ7OjaEFLEZHsaERQzIdsSxTBeC4fupKaQzgNoEJYPiGMn9nTt8Te/lwrkEn3OGU1uNnpY2oU1KQ8PJeuU9NAiSwRkeyo1CDX357wMROFnqUw3ynf7noiCXycuUYfALkeKOHhbup+6nYDtp8kytU/jYhIluRLUI4HC89VYKwEgzkIWff5eCp3uIXxM3lwp1puBzTJQGHJJmVqlmm/sSf4QhfbISIi52NhYYHZWpOxUok60FcCaBBFRviSOkQUVci1vq8aS1WYqowA0BfNk8+nqgpndPRBU5kaaWuH56f5pi62Q0REzsfw8DC7uZ/ZWv+yWxmIMCLyqYAmWASuHj/N1eOniTyi6da6VTIa1HSDmX2fmX3VzO41s99c5fHLzeyjZvZ5M/uimT1rs9qioGaZRvwzQv80IiLZUS4f4zkHZmgXTo0YK9UZNEh/YV3evdQfjowgn8vRztY7Q9lbpDtM6cY2/HYuZpYH/hR4JvAo4EVm9qgVh/1X4FZ3/2bghcD/3vhXH6j7qWUBGIrv5wglslVVWEQkC/L5ASaKNcZKTaAZ762z6HmSYCWKaGVrQjfTEvfO7AFgNDfP4CCYeWrKd9ZYN7qfvhW4193vAzCzdwPPBr6cOsaBZDGtMeD4ZjVGQU1szyhMzc0CA0CZPaM5JljqdrNERKQDNWp8eO4xzNaMsVyVegRLjNDPEg1qQBOzMIC4SKhI1peHUHG1wdBQDfdSN19CVh0AHkxtHwWetOKYm4EPmdnLCNmDZ2xWYxTUpOwZzROyM3kg4jRj65whIiLbQb1xinodIM9slCQF5hgo5Rm0Ju4wVYW9oceJUhNO13NcPR4K783MQGE4PJbV2U+hovCmNH7SzG5Pbd/i7recx/kvAt7m7q83sycD7zSzx7j7ho9cUlCTmCNOjiU5xwb0QN0CEZGdoNoY5uHaLi4dbGfY3eGhMoQMPECDk5X2ZW+iuMi9MxPx1jz7YdVCfMIpd79hjceOAQdT25fF+9J+Dvg+AHf/DzPrByaBkxvdUI2GjU0xzMKCEQKZOlNzd6KgRkQkGw5Xv4xZqCocebg1PXQ1JZ/rk32hqylsL8bja6J4Ow9kv05Nk9yG39bxWeAaM7vSzEqEgcC3rTjmCPBdAGZ2PWGE9tQGv3RAmZplytEQ5bkksPlO4NYut0hERDpx/di1uEeY5VtfR82c0X4IyyM4Duztr7aCllCLJrlo9xMmeWeXY5vV/bT273RvmNlLgQ8SIsO3uvvdZvZ7wO3ufhvwKuDNZvbrhAvsi903Zyi2gpo1GfsLytSIiGRBrQk5A2y1z+2wz9vLdLe29w/PLz9S3U/nzd0/AHxgxb7fSd3/MvDtW9EWBTXnpKBGRCQL3nfXe5iq5NgzEJFeCiEIEYqtCHjMiJdLgGoV+vu3ssWbI9rho0p29qs/p+b6h4iIyLawt34tUKDZzOGerO3ULqbn7q0gJ/0zl3PMoK/PiSIIvSL6QptVytSsKQfrVFIUEZHt4ascY89AlSgKn922opIwtLuUlv1MZXPCIOHsfu67Q3OLx9RsN8rULNOO0J+jQcIiIplxevBfgVAxOGRonGZzedbFV1TbXWpVSXHaJVOcTRrDKltAmZpl2hHue3kBT+fdXWyLiIh0amlpF48aOdHaXqzkuL28n90lyOcjomZEMxdP8I5jlqE81Ot1oggeLAOUuGQkCYqymfHY6tlP242CmrMk+cid/cYQEcmSo+Q5WDWGS06tAQ/U9jPZT2sdp5l6npG+cGzSw+QOxWThyvIYoYZNrQut3xhhSvfO7oBRULOmsN6piIhsf48aPMLdi5djqZjkZAX29AHYsqkfyYwnSC9yGaFFjLNPQc2aHA05EhHJhkoUFqNcPhzGyeUi3POMrLFWZfv49JfY7H72N3d4L0N2/3JbQtO6RUSyoJbrA3zZMgd7+yO8NcakfbFPT3Ayg2YTDg0vEhZV2NlBQdYpU3OWJGyvkqwFIiIi29tCq9spfIabhQxNWxgrGXlcsCPV/bS0BDMRQIkHTtc5NEEmbeIq3ZmhoOYsTaAMDHe7ISIi0qG+HIQsS8iwt5c7CNmagVxEzXNn5WEKhXCr1SCs/dRaOWoLWr3RNFB4Z7/6VRUIAU2ycquIiGx3k5QJNWbyhCx7HjOnUQtdUg9XwnHpT/UkW3O8BoeGAXJctgvCLCjJIgU1Z0lH6VmM1EVEdp4TVVrjZ9rLIBjTzSTrPrxsTmsy1Rvgsv6wjjdEHJ2G/gxfGSNsw29ZkuE/3WZJ/wGVqRERyYJTBlPVdrCyGGdmIsoATPSXl31VTQ8oLhahzyG5JFZUzSOzNKZmTRGa/SQikg3XDsKJch9T1eSy5gyxxHC+AF5vZXEizp4CEkVweDE8emjCyOqlUWs/ZfUvt2X0zyMikgUN4JKBJR4KiRkuHYRGlGO0r0ojymO2soZNW9g/wCUjzfi+kdV1LXf6QGFdtZeJCOlH502Pf0+3GyMiIh0qAA+VQzAD0IigXjesVXSvHaW0qwiHgObI0gCXjNRwz9N0KOzsuCDT9KdbJvnnMF7+hR9FyySIiGRDAxjJFak0w63hRaabtGrVuHs8SjLXCmgAGg2AMqWoSbM1wHhLm75hwtpPG3/LEgU1Z0n+gDn0zyMikg0DBViMilSrIRNzpgq51Ge4WfvT/Rsz4ZhqNdSoOTQMhxd3AyFLk9WuJ1H30zk4mv0kIpINS0sA85QZp1yHHM5k/xJJKJPOvlw5Grqf+vpgoQKNCsAZCrk1FojKkKxNwd5oCmrO4qmfO/vNISKSFYODsKdRJVQFDtKBjNEeUJDufhrqg28sQDXaA8yw/AzJGvWvnCX9RtaUbhGRLJiNiwAns5eSz/KkK8mx1qf7yUq7ns1UFQ6NA0ylzstmQJOs/bSTx9QoU3OWJLSfA5a62RAREenQGHCUBsu7mwz3MH4GC4OC8wWYLKVX8oZ7Z+C6iaSqcDg3q+NqdvqU7p396ldV5hCfIPwvMtLtxoiISAcawIGBs/dPV9vdTc24/+lUbYjpCpQrIWtzcADKXmydY6bxlFmlTM1ZBnmA70ALWoqIZEehAI0yrBwPubu/3dVUbF3xFqkzRB2AOvlCjcayj/uMpmky2F200bqSqTGzm83smJl9Ib49K97/3WZ2h5l9Kf759K1vXfqdrTo1IiJZcay2m6TLKZEENGYeVxU2Joqwt38xvtU4M5988jua+Zpt3czU/JG7v27FvlPAD7r7cTN7DPBB4MDWNstoR/rqnRMRyZrW4GBPbccZjMVFON3cxRVD0+TzUC6HEZRD4citb+wGCquQZ/s1XKxt1f3k7p9Pbd4NDJhZn7tXu9AaFK2LiGRPO1OT7ooK90eGndOz0xxenGjtv3zkDDW3nggH1P3UPS81sy+a2VvNbNcqjz8X+NxaAY2ZvcTMbjez26empjawWcn/DRFh6JmIiGTD6l9EfcW9ywfh6vHTXLPrNFePn6EWH/DQ/KY3UDbZpgU1ZvZhM7trlduzgT8DrgIeDzwEvH7FuY8G/hD4xbWe391vcfcb3P2GPXv2bHDrI974uL/d4OcUEZFNEwEcp0mSZ49oEuEQ1nSK0zdRBEeWRmg2w65mE6BEvd5+Ks/o4k+qU7OJ3U/u/oxOjjOzNwPvT21fBrwX+Cl3/8YmNe9cLQJyvOLOF/AnT3jX1v96ERE5fw7fNApmx3EPNWm+uHApC/UBIGK4uMRCvch4H8AC989PpE4+A/SRzHo1i4D81r8GuWhdGVNjZpe4+0Px5nOAu+L948A/Ar/p7p/sRtvashWdiojsZE2Wdz20i+eVmegLM58W6ukzThNqkc0DE4SAJlliIbuf/1nLrGy0bg0Ufq2ZPZ6QLXuAdjfTS4Grgd8xs9+J932Pu5/c+iZmM/0oIiJBIQeT/eH+VLkd6Fw5DvfP7CKEQRPAAtQM4vUss1pR2Mled9FG60pQ4+4/ucb+1wCv2eLmrGxF/DNCaz+JiGRDzaGf9synWlxmLFkuYaSUo5Brf1m9eny6dT+K4BuzkxhzuGd37SdRIZZVWOqn3tgiIlmQjIBJMiylHITxMeHnbM2WrQf1jZkQzDQacN/cBMkX2ixmaNIibMNvWbKt6tRsD+GNXeR+NFBMRCQb7lkMA4WTTM3di5AEMCH7kovXdHLMHGcX9821v9enC/ZltftJFNSsqc4jgE93uxkiItKhO2cPxMFIXBk+ciz+bjpWqhMueUm0kgQ0SRV5SDI77efIGNdAYXU/nSVdfVJERLJj+bABt1yoRROFQcPRWcv5rRxuoGEHWadMzTkpsBERyYLHjsAX59JLJBjuERjkc6EAX7Np5HLtrExyXDpTk2VJ8b2dTEHNmnb2G0NEJEtanUmpsTG5nNGMl7spGjTiQcOrFwyep700TnEzm7qpdnpQo+6nNWU7YhcR2UnO6lmK9eXCrQ4sNME9h/tal74iWQ5oRJmaVSTBTAPFfCIi2ZBkX9JZmCjy1sd4wWEo16TY+lhvEAKY5IQR9g/PxzOf0qt7Z4eK7+mqvQoH6rzu0X/f7YaIiEiHcjkwi4iiENhEUdSqTdNoQLlRZDHK0x4MnCMENh7fIh5eyPDMJwEU1KwiBxT5jbt/FHVBiYhkhzvk804u5+RzoRDdQ+VRTtZGmWvmlx0HVUJnRY7lg4Wzzd02/JYl6n5ak6b2iYhkhXu0YqxM+Pze3Vdu7Wk2091TFWAwdXy62Go2u5+AzFUA3mjK1JxTb0TuIiK9L8dsLWoFLYtV2N3X/gyv1KBQiFLF9dKf7ysDAWP1GVKy3SlT03IC2IsGCouIZE8UNRgteqsOzXA/QJOZ6lB8hFNs1MgXwjypK0fg/vmwf6VkqYSscVUU1lU7cf0+CIFNA5jj+n2nutsgERHpWMEMs3SFYAgBSxWoMVaqYQb1OpTLkG/1Nq0+1CCLQY0oU7PM9fscOA3ANXwcgHq9TrGougUiIttZw+MsjRtYmM3UaMBYKUzRdof5Bkz0QbG4dtCS7D88vWVN31BZG9i70RTUpKQHmd1VezrfzUcU0IiIZEAtLseB9cV7chQK0VldSYfnxped165LA5eMtA98xORa5fy2M9WpUfdTytISrZoGfX1NABqNxjpniYhIt+WJaM9gWt4NlfRKjRQg1B52YCb1WHJ8UremeY6qw7Kd6a+WMjQU0pfFogMPAVAoKJklIrLdVZglrN8EyeDfZjP8TDI1842IK0bnODQ2yxWjLDs2KBACnHwrEMqaXqhTY2b7zOwtZvZP8fajzOznOjlXQc0K95z4eHzvkq62Q0REOncJDxIWqzSgES7I8QiL2VrIyIz3rTyrN1bn7kFvAz4IXBpvfw14RScnKqhZ4fp9N3W7CSIicp6+dxf056HpDdyLITtjYfjA2cHMSunup5hnb+iBE6Z0b/TB5PqoAAAgAElEQVStCybd/VbidUrdvQE0OzlRfStnSSL2MuBEUUQup9hPRGQ7azSg0tzNUKFBWJMb6lGNZhLgJCw9cDgZh5N0U7UnhqTvy5ZbNLMJ4j+MmX0bMNvJibpap9xzIqQtYZrnlz7T7eaIiEiHTuTGeCK3UY2gEdUxmyeij9laCHgaDYiiHI14MsjhOVieoZknmReSLIqZOd4uHLiRty54JXAbcJWZfRJ4B/BrnZyoTE3KdXsBirjv5t3lp/My/gXP5DtbRGTnefyl4zyefwecfdEsv//QsxgphBW8IYytGS3msDxMluChSkSyiDEUOFVx9g8vxNO8O+rt2HZ6ZO2nu4HvBB5JiDy/SodJGGVqVkhimBCxG/l8/lyHi4jItuGt28ICmBVIjx4YK7UHGIQOqnLqPIiHcFBvQj6fxTo1PeM/3L3h7ne7+13uXgf+o5MTlalJKZehvz+s5BpmcitLIyKSDbPASHx/nrExYCHJzoS9c3UYju8vHzFjhM/7EAEV84Bn7/LoZLuisJntBw4AA2b2zbT7B0dZvqT6mrL3V9tEAwPhZ6EAx04Cl3e1OSIi0qEDOTCbx+KBwN84kzzSx1y9fVzUqJIvhmz81eM17p05+7mOzGS1onDmfS/wYuAy4A2p/fPAb3fyBApqEg0gD185CdftdQ7szW60KyKyUyVDCK4cB184+/FGFLI0tSji2FzrrGU/czZGsg5gtmR7mQR3fzvwdjN7rrv/3YU8h4Ka2D2n93HNxIl4sLCxtARhSneDXE7/TCIiWZBkasySJRCqQFKopkohFwYHL9FHGFlTJ4Q5DiwCcGjiFAMZHU7ZC3Nb3P3vzOz7gUcD/an9v7feubpat5zg66f3praXAFNAIyKSAVEE6XkdUav3KD16pkSptEgU5QjBDoz0Ra37zWYJqFFcdr5sNTP7c8IYmqcBfwk8D+iozopmP8Wu3wdwklC0cIZrJ1fJW4qIyLa0cq2m9rYRgpY6UMO9SLkMMA6MY9YIhfuqMFgKhWrqXqCawYHC0BtrPwFPcfefAqbd/dXAk4FrOzkxm3+1TfLIPWAW+lG/chK4PLt9kyIiO0kUsWz6djvTUgEGkr1APdVF5VgThvpCMFOtQmGAeM5N9pZJ6CHJXPslM7uUMMCpowUZFdQkToFNtDev2tW9poiIyPmpVsPM1WRMydJSyFoM5418vto6brYGC40kk2OcXByFVMG6oYHp+F72Lo+hAnBPfBl/v5mNA/8L+BxhwNNfdnKiup9i9zT3UamEN0WjAd+Ynux2k0REpEO3Vr+Fej18htdqcGv1RgAWmkWazbC/WgvH7mqNvVkEZggZHAeSgCa7gUEvLGjp7r/v7jPxDKgrgOvc/b91cm72QtFNcv2+E9xzIj1Q+FjX2iIiIudp8Q7e6Te2upaOVkL3Eiwx3xjGLCmwV6XaWgGhzsHxErncbOu8pFuqJ6YRZZSZ5YHvBw4Rxylmhru/4VzngYKaZa7be5KvnITv4GM8NHFTt5sjIiIduhL46NJu4AywO97bYLyvgHuVmbgHaqwEx6vhsXo0RBhz05421Wwm08GzOf2pR2Kx9xH+MF8iWbuiQwpqUprNsKjlCb+Jr08Rkl4iIrLt3XQQPvogQDI40mnXn4FdcbUT97CYwkwzRzFfJx3QQBhgnIy3ka65zN2/6UJO1JialKTGgRlcPtrdtoiIyIXwVe+nMxhzDuV6HuhjsVJr7Tdz+vrOPjdLemRK9z+Z2fdcyInK1KR85SRcMwGLizA2ms3Uo4jIzrb6RThdx6ZQheH+CJilr1Ra49zsZWqcrgUhG+1TwHvNLEcoMGSAu/u66QZlalKu2xumBI6NgeufRkQkMxpnlZUx1gxM+iEU5OvnTAVWzco0z94lW+YNhIJ7g+4+6u4jnQQ0oEzNWXpkkJWIyI5y3Mfjeys/xB2wZZ/txSLkowKL9T7CAtDhuHaWw6g3s7n4U49cwh4E7nI//yuygpqUSoVWf2ou1yNvDRGRHSCXm0ltJdO3I/CII1NNDk6GNaDK5ToRRWpNgCZ7czUgdEEl07obDTg+r2tAF90HfMzM/olkYS7QlO7zNTAQIvqeiXVFRHaI/Q5JVibhnsc9z9V76+RydQBKg1CuQCGCOnlORhPsZ37Zyt7FonNoYtVfs731TkXh++NbiSTi7JCCmjUpsBERyYp3zj35rH3zNYjoo+kDwFy8d5RKvYx78azj15o5JVsrXsTygiioSXG3OP1ocSVKvalFRLLh7AxFCai0ei+ScabJWonJ53vI4LRr0ySf/RkdKZzhy5aZvdHdX2Fm72OVV+LuP7TecyioiX3jBDxib3ugcK2m2U8iIllRreYIgUi4rLlD2QH6GStUyOfLrWOP1UZprZpAofW5H77MwkC+mdlJI93ofjKz7wPeRKhk+Jfu/j9XOeZHgZsJ/+p3uvuPrfJU74x/vu5C26KgJlZjH185eYJrJsJqr0cWJkPdbRER2fYeaoyxMluTB+pUeHhphIniPIUCnChDOyPTAGYxS4ZthP2NJtR1eexIvE7TnwLfDRwFPmtmt7n7l1PHXAP8FvDt7j5tZntXey53vyO++3h3f9OK3/Ny4OPrtUfpiGX28PXTezmysJfv41bCG15ERLazucpneBr/SPuSFgYMj/Q5E/0RMM/p+ignyqOEbqh4ZhQFYKKVlfF4oG3Ni5znkkPbRrIo50be1vGtwL3ufp+714B3A89eccwvAH/q7tOhjX5ynef86VX2vXjdlqBMzQrtKP+feS7fz991sS0iItIJqx/kpoNwE7cBUG3CHxx7NjNVAwbY01ehVJprHX9sHsLlLwwWjqIG+Xz7cri0BEND2QxquuAAoa5M4ijwpBXHXAtgZp8kJNBudvd/XvlEZvYi4MeAK83sttRDo4SVSteloCZ2/b4T3HNiFOgDmlw7Od3tJomISAcW7WucaZbYlatRjXKcYij1aI2p6iD7ckvk82F4QTtb40AZs3ApTOrUADSbq82O2t6cTRtTM2lmt6e2b3H3W87j/AJwDXATcBnwCTN7rLvPrDju34GHgEng9an988AXO/1FAszV+rhubyqSPwlc3r32iIhIZ/KFx7LEAEvRAAC7mQUWiBgmfFFdirue0nXIqsAwMESjMUep1A5oTpThiv4tfhEbwYHNCWpOufsNazx2DDiY2r4s3pd2FPi0u9eB+83sa4Qg57Ppg9z9MHDYzJ4BlN09MrNrgeuAL3XS0DWDGjPrJCqacvfv6uQXZcF8vcSx6VBd8sCujE7nExHZYfpyJS7Lz1JrQmnF6gZ5KxMCmzRnuAjT1XA/inua7j8NV07A5budfE7XgA59FrjGzK4kBDMvJHQhpf1/wIuAvzKzSUJ31H3neM5PAE81s13Ah+Lf8QLgx9drzLkyNXngWed43IDbzvF4BjU5sKtEpmsUiIjsMIdPhIk2ffEVzR1+3v6F19uzgX6G+5yFapNwWXNgnvlae33ESgNKHgKayCFv0GxmszLvVk9Fd/eGmb0U+CDhH/it7n63mf0ecLu73xY/9j1m9mXCxfU/ufvpczytufuSmf0c8L/d/bVm9oVO2nOuoOYX41TQ2r/V7Fc6+SVZ0GxCPp+E+Nl8M4uI7EzhSp6+oL+ltahz2HlgZLH1WL0OJyvtEyrAePyxn4sLsObzGijcKXf/APCBFft+J3XfgVfGt06YmT2ZkJn5uXhfRyuMrjml293/bb2TOzkmK5oUqDZqhGncNSr1gW43SUREOnIHR5vQiG9Hm/DY1nfTcOfYfFiocnERTlZGCZ/1DWARyLHUCI83GlBuFujPasET34Tb1nsFoa7Ne+OszyOAj3Zy4roDhc3sS5z9smaB24HXrJNCyoypuWEmhqEal6Y5vdDd9oiISKfeCQzyMMmMpUo8qjQZ7WvAKCfKSaSTXNKKtC+DS9Ra3/NrVDKZqLGeWNDS3T8OfNzMBuPt+4Bf6+TcTmY//ROhD+z/xtsvBAaBh4G3AT94nu3dtk4vLJ8GGGXyTS0istMc5EDuwXj9psT3cvb38fR2IbXtDBaamKULrnbU2yGbIO56egthetrlZvY4wpCYdYe8dBLUPMPdn5Da/pKZfc7dn2BmP3FhTd6uyoTIvr1+iIiIbHd3YAZRFGrNzMdfSJNFipcHMxY/lux3QhdU+zFwHjiV0ckiGV2zaoU3EqLS2wDc/U4z+45OTuyk1zBvZt+abJjZE2mHsBe0joCZ3Wxmx8zsC/HtWSsev9zMFszsNy7k+S/E9ftOAAOEN3UBOEEu1xvvDhGR3jaEO+Ry4TbaSrIkA4jT9WnSjHA5G42Po3V8xOTmNlnOyd0fXLGroyizk3TEzwNvNbPheHse+HkzGwL+R+dNPMsfuftaK3G+gdDttXUaSVVh4/p9ybogIiKy/bUr5Xm7R2mVLA3tB5fNcvU4KPLW+Y/ed2pTWrqpvDurdG+CB83sKYCbWRF4OXBPJyeuG9S4+2eBx5rZWLw9m3r41gto7DmZ2Q8D99POB26Je07v45qJE3FA4/zi/n/dyl8vIiIX7DgQAhqzMGUbqxGK7tmKsTawvLJwsh2ewMyABku1TW6ynMsvAW8irCt1jFCA71c7OXHd7icz22dmbwHe7e6zZvaouCDOxXqpmX3RzN4aVw0kzgb9Z+DVG/D852mBr5/eyz0n9nHPif284s7ncu+xnoh4RUR6XrVKK3gpFmGSfyYEKxZnX9K9F45ZtGw7ioifIAxBsHxGP/97YEq3u59y9x93933uvtfdf6LTmdadjKl5G6Ea4KXx9tcIc8jPycw+bGZ3rXJ7NvBnwFXA4wmLVyULV91M6JZad0K1mb3EzG43s9unpqY6eBnntmcUQnIo1C3YPVTjQ2FhURER2ebeMfskZmZCnZm/PQ353I2triRrBStrsXgxyyRjA9mdLGKbcMuOTv5qk+5+q5n9FrRKIq87YMfdn9FJA8zszcD7480nAc8zs9cC40BkZhV3/39Xef5bgFsAbrjhhouOJScHFomiIcyqAJyab/LqR371Yp9WRES2wLFakb/L34jNhm6oExEsvyBHpKdprxx7EkXQLiqfrQu5tHUS1Cya2QRx7Gpm30YovnfBzOwSd38o3nwOcBeAuz81dczNwMJqAc1mqJyAPXvbw3hGVaJARCRDGhwp7z1rb3vpBE8FMukMRDigXIFCAcw8dUwG7fBJu510P72SMFf8KjP7JPAO4GUX+Xtfa2ZfilcCfxrw6xf5fBftfvaxtBT+B6jX4b6ZvWgGlIhIVtTjn8u7TJYnXRwz4psvO95bx2avy6XXmFnTzP6nWfuvZ2af6+TcTmY/fc7MvhN4JOEv/VV3r69z2nrP+ZMdHHPzxfyOC3FkYS+0RvNEKKgREcmGVx+8g9998MBZ+5NMzUAJyjVopzKWd0ctMcIY85vdzM3XG5mauwlJlw+Z2Qvc/QwdRpprBjVm9iNrPHStmeHuf3/+7cyCZKpfhHuRKIrI5bK6spmIyE6z+lU9BwwW29/H3aFeK5Lulkqv8p3JYTUO9Eadmoa7/z9m9gLgX83sp+gwXDtXpiZZ02kv8BTgI/H204B/B3oqqAmF9yYIkXvEntEwe6wRPUQpd3b0LyIi21F6rExSdyYELLU6lIrhfrnRPj7MfMpoINObQtlE978xs7sJa09e3smJawY17v4zAGb2IeBRycBeM7uEMM2751y/7xSQvLsd94hG5JS63TAREelQ+wt9uqKwGTQo0ajD8uJ7KzM1nhpbkz3eG91PP5/ccfe7zOypwLM7ObGTfpWDqZlKACfoMGLKmtYbGyOK4nVEqHa5VSIictHOebHvjUigV7j7Hcl9M7vF3Wfd/R2dnNvJlO5/MbMPAu+Kt18AfPj8m7m9nSoPsbuvPaX7TDUHOP2lq7rXKBEROU/Lp2q3VuU+6/H0DCnD3ZcX4Muq3ovPbjifgzuZ/fRSM3sOkCz7fYu7v/dCWrbdnakMpt7nNS5wEXIREemKlVf0ZJVuSHc5tQOXsLBle3Xu9s/Mjq/JelB2tpPnc3BHdaDjIKYnA5k0ZwljAGjwm5d/iM5650REZDt4Ou/jI/wAAO4zjJfGman1k9SeGSpVWayVMMsDTgFouAN5avVQVTis1p3hoKaHxOtBPu98zlnzqm1m71/rsfM5Jium5oaBIdxzuJf47bt+EMgxs/BP3W6aiIh04KbL4dUH38+rD74fs91MVyPG+2q0MzhN0hl4M6eYK7e2Ty6NcHQWTiySWeYbf9vy12D2WDP7PKFezd1mdoeZPaaTc8+VqbnRzG471+8FHnUe7dz2Ts0PpbYagBHVxrvVHBERuUChi6kP94gQzKTXvklfqUOR1Zyl9zfI7oKWPeEvgFe6+0cBzOwmwlqPT1nvxHP91TqZPlXrpHVZsGd0gam5PsI/ScTuoTMAHG98lN08uattExGRziQ1ZwBGSk3C9+/Q3QRF2sFNBOQo5YvUmhA54E5faYC9gwtnP3EWOL0yUHgoCWgA3P1jZjZ0rhMS56pT8/GNaFl2LDI5AsRTuOfjrM2IH+pai0RE5PykF7Ccrxm7+h2YBibj/clgmRDUzC5Cvpg85EAZN42n7LL7zOy/Ae+Mt38CuK+TE/WXi00OwKn5Rep1ODVf5YrJRcAZGND6TyIiWRNCl3CJOzSQfiT5TC/gDk0GAKcvrrK6fxgKmb0yWpj9tNG3rfezwB7CygV/R4hIf6aTE9VpmHLtJORyi+wfga9OAZfnqNa/1u1miYhIh5LuJ8fopwo0OFzZ1do/Vqgy2ygQRk/0Mdxfph6FwAbmqdfzlEpNMnt57I3up2e4+6+ld5jZ84H3rHdiZuPRzZCsW2kWAhyIKNQe180miYhIh9zTt4gKJUbiGU9Jt9RsI6RtRkoRZk4pHiDsDlG0i9PVQbK6REIP+a0O951l3VDUzL4duBm4Ij7eAHf3R5xHAzPC48qS4T5AtfBx4LndbJSIiHQgPUjYzNg3UGHeB4G+1GMhYGlEYcBwNSq29pmFmjV9uYhqlF/59NmQ4UyNmT0TeBZwwMz+OPXQKB1Ww+0kv/YW4NeBOwjz4npWUmzJzMnnQ6XJKOqZCV4iIj0rvZDjVBQ+z3NWB6/jPsDyAcL5eJVuC4GMJwFNBTCqdVs+A1y2ynHgduCHCDFHYp4Qh6yrk6Bm1t13RgU6J1VpyGk2jVzt0d1skYiIdCBdAXhPPJSgEg0yWVpiptoOevYMVphaGiJM725/3ofP/372D8/zwJk9HJqY3rrGb6QMZ2rc/U7gTjP7v4Qo9Nr4oa+6e72T5+gkqPmomf0vwijk1pLV7v6582zv9tdabr69EFo9/xHgZV1rkoiIdObspQ2MU7XwE2Cx1u6CStZ3cg9Z+eWlczMc0PTG2k9PAd4BPED44x00s59290+sd2InQc2T4p/plTIdePp5NnLb81YKEmCBxUWn0renq20SEZH1pbufmk1oVxFOZjYZUCCKLB4UDPm8x91P7ap17nBoIjleuuQNwPe4+1cBzOxa4F3At6x3YierdD/topuXGQZUyHMfT208zOgoFOnrdqNERGQdZrC0BNUqvHHhO2gHNWXCGNMxoEmlWQCcfD4U34My7mGRy9FcmShqz4TNom6s1bQJiklAA+DuXzOz4rlOSKz7pzOzMTN7g5ndHt9eb2ZjF9Pa7Sqs0NpPPXoUt0VPByKcu7rdLBER6cDQEOzeDa8++Al2FxaAOgPAwVIyHCNcF0M3VSjC15czoqifaq3MifIoZmHdP+mq283sL83spvj2ZsIA4nV1Eo++lTDy+Efj2xzwVxfc1G0sXadmKA+Qo0FHC4OKiEgXpWvUANQbMFmqUQaOVCeWdU/FZwBQjQo0mnkm++D6vWeoVpNxlRlNefgm3LbeLwNfBn4tvn053reuTsbUXOXu6UItrzazL5x3EzMiXecAnAL3dLM5IiLSgbeevpGfHv83cjl4+zT0leBMDSZLMF1ZbRAxQJ2QlZlm3gfJL84xMpDULtESOd3i7lXCuJo3nO+5nQQ1ZTO70d3/DVrF+Mrn+4uyIB3Jf+0UcHlGI3URkR1mqQzvsBuB8Fl+qtYHDAJLK76sVgiDhyF0RzXoK40AZxgdzHFkDi4Z2eLGyzKrFP0FoJOiv50ENb8MvD0eR2PAGeDFF9LQ7S4dyV+9GyCiwRXdao6IiHToFGOcWkpXzHMmS0utYKYd1PSvOLPAI8bPUK3CkblxoEbRlmhk9DttjwwUvuCiv53MfvoC8DgzG4235y6khVlQr0OhEN78sw0YoMzJxte73SwREVmXEbqMkm+nVdzhdD08trz7Kam470CNZhOOlXcDTfpZIvLMjqjpFRdc9HfNoMbMfsLd/9rMXrliPwDuft59XdvZXA1G4gljZjBsUGaEPno2hhMR6SFlYCi13eR0fYgQwDjeKkrnhEvfPAOlfoiaHJ7bDQaXj87TaEDOIJ/VCVAZLr5nZk+I715w0d9zZWqSd8dqvYs9GMT2MV+vkgOaDpYLacyJgorviYhsfwVgMbU9eI5j5+kr9oHXgWIruXNkDmCMy0dnN6uRcm6vX7F93kV/1wxq3P0v4rsfdvdPph+LB/H0lNFSlblagYhcXFW4Svg3VPeTiMh2t6dvmqnqPtrdT4vAIBPFKjPVkbj7yYEmfcV+FueajIwUWwtbhsfGgZk1ZkplQIZnosPGFPvtZKDwnwBP6GBf5o0W66maS01C/+y6g61FRKTLHsNn+FLpu1rb1Rr0lZbiJROS7icDckzNDQFNCo064XM+RAJXjKXXfMpoZJPhoCZhZhPA7wI3El7RvwG/5+6n1zv3XGNqnkxYVGrPinE1o/TqouzLFrTM0Z7sJSIi29lNe+GjD8LuAlQbsMggVltiJhpOHVVjsGjMlevsH62wVIdH7prhrtP748eTz/+MBjS9493AJ4CkRt6PA38DPGO9E8+VqSkBw/Ex6XE1c8DzLqiZ21wUGblcHLHzccA43sjx2O42S0REOjLImUY7IJn1YebqqUeL4bH9oxWqcUBTbz3eaB+Y4cCmR6Z0X+Luv5/afo2ZvaCTE881pubjwMfN7G3ufvhiW5gFoUBTeCPf709ncfEjDA1NdLlVIiLSiVcffF/r/u8++F3s6YM9fVCpLHC0vpuhQghcTiwOkwQxxSJcUXqYwzUgWcA4uzFNr/iQmb0QuDXefh7wwU5O7GRMzVI8terRpKoWufu6o5CzKKk86Q5vn3kivzL0+W43SUREOpAM8HUP3QvJ5/nR+u7W/pNLydHN1uMjI3BdA8oYHlfpy2xM0xuZml8AXgG8M97OA4tm9ouAu/voWid2EtT8H0Jf1g8AvwT8NDB1Uc3NiF858Gk6W/NTRES2k3mgr7XApRNFeU5Wk5EUZ08TqkbEH/eG+1prRWVADwQ17n7BC1V0csWecPe3AHV3/7i7/ywdzBXPnMbyzfA/Qlbf1SIiO1OyHMLuQpOH50lN5V5pxXyX1CHN8yrML9tJJ0FNMozqITP7fjP7ZmD3JrapK07Vh2g2211Ppyuh9uAuVRQWEdn2ksxK8jNnxv4RqFTgquGZdc/rS8U4yYSRrDHfnFuWdNL99Jp4MctXEerTjBIWmuo507V2ie1cLgwaHih2sUEiItKRpMsoydRU6nWOVi8F4FDhdGpBy0ShdR7AiVkYHXPMIJ/XSOGs6mRBy/fHd2eBi672t52dmoddg+FNnktyWI1zniIiIttMFEGp0AfVJjDDyfpqY2SW77hkFyxG8UaWY5oMr/2UZmZ5YB+pOMXdj6x33rpBjZntIYxEPrTiyX/2Qhq6XU0OLBJFQ603/qn5BQCOM6o6NSIimVNlsBABI0zm51lY46gku/PlabhiLOxbahqDnfRjbEcZ6y5ajZm9jFBR+ASh5DOEV/ZN653byZ/tH4B/BT5MWDugZ+0ZXGgFNXsGPdsj4EVEdpikK6leh0YDLh0IYyIfngvZm+VXfF/WJfXo3XCymh5voDR9F70ceGQnyyKs1ElQM+ju//n825RFIYJxB2u92zVQWEQkC5KsS6EAcwwz4XVqNZhnIgwpaCbLIADYKl9ak8fqZLWcR9YG9q7hQcKQl/PWSVDzfjN7lrt/4EJ+QVbcc2If1+090YrcF8u5VQaWiYjIdhc+u+vcOz+Z7Ek9mkQy0Vmf8QP5ZhzoZDOgybrUOpP3AR8zs38Eqsnj7v6G9Z6jk6Dm5cBvm1mVEL4a61T0y6qvnNwLnAT2As34Dd9zL1NEpKcsNWAgNSU7imDMYBpvZWPa3U/pbE0PyvZLS4ruHYlvpfgGHb6yTmY/XXBlv2zaCzhvevx7ePkXfog/ecJfd7tBIiJyDqfJc2nUJJcLwcstMzeRy6XHRDpmRjtLE36mu59UcK/73P3VAGb2fHd/T/oxM3t+J8+xZo7NzK6Lfz5htdvFNHx7OpG6b7z8C88FHmZorcNFRGSbaHLcxzjaHOO4j/HDE58kn4erR07HtzPU6tD+sh+WSUi6n9yhXE4/nlG9U3zvtzrcd5ZzZWpeRZjK/fpVHnN6bKmE6/fBPSdOEooll7l2sszXTu3uqH9ORES657I8HG3OAkPAIpfFXVHpBS7Xm8s0NATzrYMyHNxkuOlm9kzgWcABM/vj1EOjdDgdbc1rtrv/QvyzpwvuJe45sY9H7jmB2RkAvnISnnPFPzOrMTUiItveZXkwWwQ470ke6UrEmY4Ksu84cDvwQ8Adqf3zdLiSwZpBjZn9yLlOdPe/7+QXZMlXp/Yu275xh40mEhHJsnR30kwNxkurHrXsZ/ocW+WozMlsw8Hd7wTuNLP/4+4XVCjoXL0rP3iu3w30XFATJG9r5x/mb+D5E7d3tTUiItKZJONiBkOWzHiCbyxA+3KXzH6ys85ZHg+o8upWM7Nb3f1Hgc+bnT2ax90vvKKwu//MRbYvo5J/x4jlg4dFRGS7iqL2mj+HS7cAACAASURBVH2NRpjN9EB5svVYGJKxauoGSHc/ZTuYyXjxvZfHP3/gQp9g3QpDZvbfzWw8tb3LzF5zob8wO3KEtbRERGS7y6WuZoUC9PXBNaOnuGb0FI8cPxU/0lnA4n7+43Lk4rn7Q/HdZwAldz+cvnXyHJ2UTXymu8+kfuk0YXRyj2pXmwSY0DIJIiLbWjLLyax9v1artgKTSuXc56X2rLFfttjlwF+Y2X1m9h4ze5mZPb6TEzuZsZw3sz53r/L/t3fnYXKd1Z3Hv6eqel+0dUvYkpAtYnlhd4QhiTFLgLAMGEgyMRkImSGQCVuSCcyQkCEJDBMgCUyYhyVOAmGGDGYHAyZmMxAIGGzj3RhvYEte1C21VL3Xcs/88d5bdat6Vau7q2/17/M81113rfdWlauOzrsBZtYDdJ1CYTe4JDyfAKY56dm0RERkXbmHKqdCob4+OdnFvbNDqaOaf7asIRvjHv7jGa9+ynJD4YS7/xnU4o1XAm8E/heQX+w8WF5Q88/A183sw/H6fwQ+srKiZkEJOA7sALpbXBYREVmOjxy/kBfnv0NPD/xD8ULMYLIUfuTm70YTfv1nZ6FUgYE+8Gp9n7I1rWNmfwr8EtAP/Ah4A/Cvyzl3OdMkvNPMbiDUcQG8zd2vXGFZN7iI0JBsGBhhB9eyo8UlEhGRxZnBdASftQthMuRk8nmg2kllTkft+nqlEtredHXBxARYd3oahQymPFo3AvBqezEhFv0S8C3ge0lt0VKWO2DubUDF3b9mZr1mNuDu4ysr60Y1AbVJEQzYwoUD0KUhhUVENrSkGmmqWu/ddHeptrf56NrfG0b2NGw/q2s0+xmaNghq3P18MxskZGueCVxqZkfc/cKlzl1O76dXAp8C/i7etBv43CmUd0M6d9ckMEbI1szy14/8LH9087LmzxIRkRZ73a7vALOE7/BjPKy2Jz2RZSl1hgGH4uOrHBh+iGrDrJZZj26yy8weBfwH4OXAbwCHgW8s59zl5CFeA1wAXA3g7neY2c7FT1mcmf05ofHPSLzpT9z9injfYwgB1CDh0/YEd1+g7frqqU+TELr+veGWpwOjzK5oTEMREVkvSXbldbu+izscqsJn/UIeGGnsyVScLjDYU18//2Fgdn9tAL7ZKpTI427Zzdi0QaYGeAfwbeC9wA/dvbzcE5cT1My6e8ks6epmBVbnZXuPu/91ekN87Y8CL3P3G8xsB7DsmzlVjdMkOMODdzG1Xk8uIiIrlu7JNFx7ZE2D6lWADprbzCTHlPKpx/oHbcu4+4oH31tOUPMtM/sToMfMngm8GvjCSp9wCc8Cboznf8DdW9Cjuv5hP3cbTK9/AURE5CSlZ+QOVQBdhJH2rbZ9W5dRrTUWtoZz3aHLYc2rBdaQ0TYNhVdsOYPvvYnwGbkJ+F3gCuBPV+G5X2tmN5rZh8xsW7ztAOBmdqWZXWdm/3Whk83sVWZ2jZldMzIystBhp8ABZ9uSx4mISKs8+OCDc7btzqVn3g7f5WZOteEXb4bpuD9NbZA+4I7REOwoUZNNy+nSHZnZ54DPufuyowcz+xqk2mrVvRn4APA2wqftbcDfAP8pLs+FwBOAKcL4ONe6+9fnKdelwKUABw8eXIPYNHz6o9W/sIiIrIHmqQ3Sk1WmjgIitnU7M56nWkk3Du4CBoBxepYc5m2D2uSZmgWDGguNaP4MeC3xL7yZVYH/7e5vXerC7v6MpY6Jr/n3wBfj1UPAt919NN53BXA+MCeoWVuhCupcruLE+j6xiIichNnyR2uPkyDm06l/flstonHCOGSTAIzNwJl9VWbppD49wiynD46nsjwZ0ybj1JjZAcIowvtIxSnu/vSlzl2s+ukPCX3En+Du2919O/BE4JfM7A9PscCnpVZfBNwcP74SeHQ8Fk4BeApw66k818nz2nJW/+D6PrWIiJyUoe3PqU1AmQQiT9oxED+OqH+nQ2Mao4N7JmGgUGKgMMtAYRYc+jozGtC0l08C1xGaurwxtSxpseqnlwHPTLImAO5+t5m9FPgK8J4VFxfeFU9O5cBPCW11cPcxM3s38MN43xXu/qVTeJ4VSBoKG4UCDGW1W5+IyCYwy71AY0PhPfnxeN0IHWhzhACng1DFNAuUOH0AxivQazDlgHWBO2aW3cCmBeU2s2cDf0uYm+kf3P0dCxz3q4Rx757g7tcscsmKu39gJWVZLKjpSAc0CXcfMbOOlTxZ6hovW2TfRwndulsg/WmY4vpj8PPDCx4sIiItNlu6vaHNTGPVUdKFG8LvbRLg9ADTFMdhcKCTKU8PzlcflViWZmZ54H2EkX8PAT80s8vd/dam4waA3yce824JXzCzVwOfJUSgALj7saVOXKz6qbTCfRmWTst0cIf/fMtKIiIiS6tOdxJFzKmCcodOK9PbUV+CPCHYMSbY0lD91Jer/8M204PvrfayuAuAO939bncvAZcBF89z3NuAd7K8XvMvJ1Q3/RtwbbwsltmpWSxT81gzK86z3Wjb6atnCKnJKpDHLM/QEmeIiEjrTPEvtQAkXQVlBiXvwqoRuRzx6PAdQJVt3WXMYKjrBMUyDBTCOZNR9n/aWtBQeDdwX2r9EKH9bY2ZnQ/sdfcvmdmSbWPc/cyVFmbBoMbds9qhbUWGBycYKSYTWhaABxkfh1FVP4mIbFiz1aO1ICadpQmPp5mN+iFK94BycjkYm4LeKtDZxXjFaudHUTzDt6QNmVk6U3JpPKzKkswsB7wb+O3lPlncxOX3gIviTd8E/m450yVoDuqUENhUgaMMD+7ghuM7Wl0kERFZRBRtB34GhDmf9uShVBtgbAbojx8nKYyIo1O9wBRjhEkG03LLGZJ2I1ubTM2oux9cYN9hYG9qfU+8LTEAPAr4Zty9/mHA5Wb2gkUaC3+AkFZ7f7z+snjb7yxVUAU1sZFiP8ODxxkezAPDjBSngS5O1yskIrJhzeJ0xY/3xtVIXbVMS/Os3E5oYuAM9VXZ0QknSl4byyaz49O01g+Bs8zsTEIwcwnwm8lOdz8B9ZYcZvZN4A1L9H56grs/NrX+DTO7YTmFyXpMuqpGilsYKfYzUhwA+oC7W10kERFZRBd9VJqGfp+MBwnuq3VkSkcqYfC9yckQ+XgqikmmVXCH4izZsxaNhJcI8ty9Qhik90rgNuAT7n6Lmb3VzF6wwjupmtkjkhUz209o7Lok5SEaHAG2A0VgKxfuzX6jMRGRdpbL3Ys1BTV54K4iUJu9b25XpumGn79kkkvHDMayGNDEWjGisLtfQZgXMr3tLQsc+9RlXPKNwFVmdjfhjdkH/MfllEVBTYNhwuu3HZhkZAQ4bfEzRESkle4L4wanGgmPAhE7gKOEMWnS+gBnW3clPt5rIU8YdC9UUW3ryvJ83dnm7l83s7OAs+NNt7v7skJNBTUNGsepuX3m4Uwu2dZaRERaqTPVW8ksNBYO3+db463N1U+9VCqhvc1stRB6O9WCIst2QJPhNkFm9nR3/4aZvbhp18/FAednlrqGgpoG4cMeBmaa5t27Ps5Yi0skIiJLm65Cdy49K3eFA1uO86OR3viIpKEwwBTjlTyHJst05MuUSh10dIRzZx16c1kdeS/zngJ8A3j+PPscUFCzXI3j1HQAef5y+lm8ia+0slgiIrIM3bkQzJjBkSoc2HK81vB3bpua0KV7Tx9c+0CVPdt6qMbtcsqTFdiyrkVfVVmepdvd/yx++FZ3vye9L+5dtSQFNSnDgxOptQpDPZ1o7D0RkY0vPa3BzvxiXbOTzE1oU/Pzp8G1D2xN7Y/YwpJTDMna+jRwftO2TwFLzl2koCbFvUqYmysiiow93dfRd0pTd4qIyFozg+kKdMfBzEi01JgzJbZ1W8N6mMjSOTCcBDQZTXlktNgAZnYO8EhgS1O7mkGWOT2TgprYSBGGBgrx/wR5ymV4xfb7W10sERFZQujxNAjVJFCppgKaiNDJOzHFjl44OpVnZ0+YFuGsoSKhO7fF1/NsTmi5vAkoN7KzgX9HaOGdblczDrxyORdQUFPTz+j4BAOdMF4C6OYt972Qv9z/uVYXTERElrC3UOS+SgGoxL2fFnZ0CkJmZjLe0hwJWMOgfLI+3P3zwOfN7Bfc/XsruYZGFI6F9jQ9jJf6CHWuR/nJKPxUXbpFRDa0JKuyt9AY0JTLEDI1SQojIny/9xACF5j2pIGxNVyvVMpeqsbWaGmB/2xmtYZOZrbNzD60nBOVqYmVizA8OJXakmOkGC14vIiIbCzp5IoZdHQAjACnx1uTAWkqQIQZ9BhMeoR7Lh58LxyZ/jWQdfcYdz+erLj7mJk9fjknKlMTO04/I8UpQiQ/w0hxiDC56PbWFkxERJaluR1M43pzzqFAuUoYeA+rHZv83dZFNq3z3E9rJGdmyRwXmNl2lpmEUaYmFsapGWakmHzwZ3gMXwJ17RMRyazHDcH1o8la4y90ZyG9vbH6KauyPE5Nyt8A3zOzTxLemF8D3r6cExXUJMowPDjZsOHG4ouAf2xViUREZBma2/S6z9ed21N/yyRtaoKIUHHh8YjEGY5q2oC7/x8zuxZ4Wrzpxe5+63LOVVATG5nuZzh3FPIdQMRIcQvn7voEZ2icGhGRDW+hqqd64JJMk2CEn778PBkZy3SWBsh6l+4ad7/FzEaIx6cxs4e7+71Lnac2NSkjk9sZKQ4wUtxC+OBfoN5PIiIZkWRoKtX6evyIxl97A9KTPlfXo3iyTGb2AjO7A7gH+BbwU+DLyzlXQc0cR+IFdnSCGgqLiGSHe+jblK5+Gu6dIoxJM8lw7+QCZ7aJ9mgo/DbgScBP3P1M4JeB7y/nRAU1c+wCdgLOj+5rdVlERGQpzVVGnVaf3BLg+BQM94alOAVh0uLO2vED7dIQw0ND4dVeWqDs7kcJvaBy7n4VcHA5J7bLW7mKkndwhkn28MCx62B/SwskIiLLkA5u0pmaMj2MTKX/DV8GTsxzheYqKmmR42bWD3wb+GczO0J9+OdFKVMzR5Xw2nVz7q5j/IC9rS6QiIicgp5ChRDIVOO/0DgfVMJSfzPa4rY9qp8uJox/+IfAvwB30TgX1IKUqWnghA96H3AEpqegu6fFZRIRkZOV7v00XTF6OwCieL0A9APH5zuTzAY0bcDM8sAX3f1phL72HzmZ8xXUNEinG7cySifDjLesNCIisrTljFMzVW4en+NE0zHeMAZfViufsj74nrtXzSwysy3uPl8d4aIU1MR2bzvO4bEOwkRnFYYHpxgplvjgOVe3umgiIrKEdGYmcojiqfvun4ZOK1NyqP/kVRjqy2NWqZ0DhmNYnKWJMh4cZNwEcJOZfZVUWxp3f/1SJyqoiR0e28rubWNAKWwolVpaHhERWZ7m3k8OPDCzpdYDqqMjouAQOnvDVLlMUhVVPzdJ09g8oxFnSJbLXveZeDlpCmpSDo9tTa1VgHGmW1UYERFZliQISQKUfFOQ01MwnPpIqt0FCM01Gq4y7zWzJsvVT8mowe5+Uu1o0tT7KbZ723GSKB4idm85xmO2fZ2OymJniYjIRpDO1tzvjdst/bg2fs180yTUT3RvDnpkHXwueWBmn17JBZSpiZWqcPrWejf4UtRFufMJHJ64kse3sFwiIrK4uVkVm7PPUpNVmjlQbap+KmPWFY41J5/Ff/K3rgv2akmHmSsaIS6Lb9saKVCOqC1BF5+bvaiVhRIRkSXkmn7Jdufm/rKnOjYtyU/mYFlNvsDjZVOmJjbUM8nodDfJ9PMwzVO4nF/Q1E8iIhte86zc6ezNbAU6CtYUqMwdj8bd4+opz2ybmoxnah5rZkXCm9MTPyZed3cfXOoCCmpiYw/Bjp0ztfXZWfgWz+Os0pdaWCoREVmOdBASRZBLZWu6CkbVAbypt1P9vL4cTFRz9e7dGQwOjGw3FHb3+YZ5Pimqfoo9yC6OHQv/M8zMwD0nhoA8ZVOqRkRko0tPYFkPXIwkPRO21dfT59X/1oOe3DxVWLLxKVOTcqQyzJHR5AN/BHiIK8bhxa0slIiInLLGRsHzszjNkQQ4mZTRYq8WBTUN0p/4bdz20A6G9+klEhHZyJqriqII3JMMTWj1G4KaeRoQp9rg5LMczAigoKZmeHCCkWIH0AFEbO8rcWzyBzxv8GiriyYiIotIZ2Dc672h6lVOtSOToxrOTbI41ajeJiOro9RYFhsDrSIFNSnDgyXcw6iTo+MAT+QnxStaWiYREVlcc9XSMQezKkmIEkVGPud4qsfTjs76uakrEdXCmgwGB9kfp+aUqaFwE7NJRscngQq7t01xB7tbXSQREVlEc1uZrQ6P7P08UASKzMRpF6NeFXV0Fu4aDdsPTSbBTZlNHxVknIKa2EixH6gCvQwP9gCzDDLN8wbvanHJRERkucygUIAX9cFgZ/iJS5rKeO0Yw+nlWGlXw3m9bTDgnvnqL1mi6qeaI4wUh1Pr+RDjF6utKpCIiJykpDrpw5OAR+ADxEPP0Dj+jAETVKv186Yy9gMucylTEzt3lxO6cVeBImcPTwBFtm9XUCMikjW/uw36UkO5pSe1DJycQT4Pp/ck2zuZe0bG+BosGaJMTaxY6uLs4VnMQm+n4zNw7rYfc4am6RYRyYykN9O/jsHd49BdyNFbgHIZ8h0NR3LGjmmq1ebB+pKRhjP2ax7LWnXRalOmJlGF8XIXxVJYcrkuHpw5v9WlEhGRJbg3LgAXboWIbbX9J8r1YxN3j0bMzECl1nsqio8x9POYTXrXYiUKTE7W/6coR/D7vV8kyupgBSIim0QyRUJ6qoRCAfb3jeEO45WOpN9T7Tt+puLADrq6oDMftvVamfq0ChlNeWzy6icFNSmdPXkqHhaAtx1/Nt/V2HsiIpnjDndPJm1kIoZ7w+9zzpxczukuwGlbjuG5ejDktUkvM/hrLoDa1DRIhtUGqFQKjBT3cN3gk1tYIhERWYmQZS/h3kPBUtVK9WYzgNOZakyc7haSy2I74Qx2wV5tytTERor9VOI2wVEEx6d7uCj/HS7efkNrCyYiIictl4NOtoTxZwrhl362mq5+AggRTbKtlOr9FG3y4CCrFNTUPMTx6T5Gx/s4NtkHTPDt6hMJI1KKiMhGNV9HpSgK487Uq5VC1UTenLx5HOhUm64x3zxRGbPJ29So+il27i647aEjwBbgBGftqPCqPV/RJAkiIpnVQZj6IMjljKmKkyP0c+rpaD4+NCeGbIY0hqqflKlJOTAE5+46zrm7Is4u/BufvgEO09vqYomIyCKaZ+l2h9+/5xKSgCaKQhDzszHoyhtdHdDTEbG9o/ka9Z9EVT9lk4Ka2Ox0qIN1N9xz3Dx9EQ8M/iIw1eqiiYjIIpLqpeSvGbxn32XACQCmovqP3U+Pbefu0R3cPTrMNQ/01K4RGhanMjVZTXk0D9qzGkuGKKiJjdPHkSPhgz09DZN0M0UXe7KYgxQR2YSSjM3MDFx6/Jfpi3/h8vGObV0A44S8TZH9Q/Uv+Jx+DduC2tTERor9DA9EjJVCl79qqcquLSO87yj8w5mtLp2IiCzEmv7xeTju3zHrW+IGBE4Uwdiss3+oRKiWylYGYrmymmBaLS2JTc3sz83ssJldHy/Pjbd3mNlHzOwmM7vNzP54Pcs1Mj7ASLGfkWI/x2Z6CY2GH7meRRARkZPUXFNyRTyAaiVe784lmZikjaTX/ibnVGsdoTIcFaxFz6eMvRytzNS8x93/umnbrwNd7v5oM+sFbjWzj7n7T9enSLOEWVojIM/vPOwL3HV8fZ5ZRERWprmh8GkGDzgQ93Myg6kpGOwpJWeQtJ9Junzn88Q9vBtG55OM2Wi1iA70mVkB6AFKrNNAMcODE0AX4QOdB4pMTm9n33o8uYiIrJoH/KlUq1CwkH6ZrcIkfUCFxvRDpXZOvT1stgMai1Z/yZJWBjWvNbMbzexDZrYt3vYpYBJ4ALgX+Gt3P7ZeBRoamASOMDRwhG29HXzsxC/xJZ6yXk8vIiIr0Nz76Y5inHlhBgi596HeSbZ2QPiH61S8FOZcI5N1LlKzZkGNmX3NzG6eZ7kY+ADwCOBxhADmb+LTLiAkAE8HzgT+yMz2L3D9V5nZNWZ2zcjIyCmXd6TYT7UKQwN9uPcxNtXDaXyBR/GtU762iIisnaT6Kfl71iBMTsGZvSUqFegtgEdwvAwhmOmNl/qQHW3T+0ltataGuz9jOceZ2d8DX4xXfxP4F3cvA0fM7LvAQeDuea5/KXApwMGDB1flZR+b6kutTfIAz+Wl3VesxqVFRGSNeNN0CPdPD4LB9GSR7g4oVaCzANu6YXsHmE3N6TEVpapZDKOasWqXhHo/tYCZnZZafRFwc/z4XuDp8TF9wJOAH69fyU4Q6lhPAD38zun/xj/NPHX9nl5ERE5ac6YGyrgXKbGFqQpUyDE1E6ZJWPwajrsxU2mjzM0m06reT+8ys8cREls/BX433v4+4MNmdguh4vPD7n7jehRoeHCCkeJg/LRbgCMwdZwQZ4mISFa4dxJFnRSsg/AzU6Gnx+J9c8e1qQu9oboKmRtIN3AyWvDV05Kgxt1ftsD2CUK37pYYHpxMrRX44OgLOa3/1NvriIjI2mn+HZ+qpGeqNMBqx5TL0NU195xa9ZODZ3I6SwGNKNzgRLHKlsEcUGZLocJUGcK4NSIislE1j1OTZGfC8BxVevIQeQhscoX6OenAJpcDr4SJLHO57M79lNFirxrVGsZGiv0M9Bdwz+Hexb3H+jhc7AKuanXRRETkpBjQAVTpLQDxiML5PBQX+dE/Xu6hWOmIA5pNHh1klDI1KUcn0r2fKsTDS4qIyAY3f1OSMKJwtQz5gs97fNLGJopgW9f0Whdz7W3yWExBTSw0FO4kvCTO8OAxtvNjHtXqgomIyLwOHz5ce5xUQdWrn5y+jihUOeVSB8Tj0zQ3Fs7lgCjZ6Jlsb2uo+klBTcrw4Cxh/ifwahe3Tz6Rx25bxx7lIiKybOXojnm39xbKWDpqqTW08FomI90LqlSCQurX0D00LpbsUVAzRzV82PNw7q4xHtQ0CSIimdKQhbGm7MU8sUoUtUlP6PRU5ZuUGgrHwjQJYXZusxyD+VmgSA83tbpoIiIyjyk+uMAem3fV3cLgejPJeghoenoaA6GFx7GRjU6ZmpRjk4O1xyMUGB7s5Be7W1ggERFZUJWjDeuNk1Ja7WE6d1EqQ1fcJ8Qs9IhKgpt2SHKoTY2kTALdhJld++jMz3De9nWbJFxERE5K4y9489gzyTFeNSJgNoJogZ+9XA4sNd9TZgOcrJZ7laj6qUEvYbCmPsLsrQ+1tjgiIrKgrnzvksckVUnTUScRnRBVFjym0SaPDjJKmZoG6U92B4fHDsAZ17asNCIisrDZahjxff6sitf2FQrOVkJDmiiCajX0dkqfV6nSFnU3bXALp0SZmtjw4AQhOxOG197eV6I+KqWIiGw81wMh05IskAQrVlsqlbCtWoWpqjXM6p0shXzjOerSnU3K1MRGiv0MD47jPgXA6DjUgxwREdl4ZhrW6pkXa5zXKW+4haE6unMwW4XefON59UH7DPD6BJdZ4oTJqzYxBTUpI8X+1FrEubuOAFn8ZIuIbG5Jo+FqFXJxD6ckRTNNaEGZFkXEyZkQFBRnyKbNHdOo+ikRqp+SBmQOTFAuA4y3rEwiIrJ88zX4tXrP7gWPAThWrVdbucNDkzvWoISy1pSpiY0WYWigPk1CcTzPtJeBnpaWS0RElidd5ZQ8nqqGgTpyOattN+sDJudpYGyp9jhHm3dmghoKCwBOP6PjMDsLo+MVSvRyeGwQ6Gp10UREZJnMoFiG2XL4G9G5rPN2FODBEWdmBm4fcR6xbY0LKmtCmZoGfYyXQiOx0JbmKHC8tUUSEZFlqVctdTJb673kqb8273QIyUzepw2DmXN2F9xxdPu6lHnVZXbUwNWhoGaO5ANRAvaxq5VFERGRZUv3fkptJZ+naXZurx2bBDaVCvzk6HDqPLWnzCJVPzWYJGRopoEir9/+WR5ia4vLJCIiJ2cGqAIVtnSWwqbmOS6b1vN5CJn5CJjkwPD0WhdyTZiv/rLkc5o928xuN7M7zexN8+z/L2Z2q5ndaGZfN7N9a3HvoKCmSS/hJQmNgz9WuZBN3z9ORCRj+gy2dpbZ2lmhku7UShLM2Ly1NOcMlzh35yjn7pyvEXEG+BotizCzPPA+4DnAecBLzOy8psN+BBx098cAnwLedUr3uQgFNQ3Sofs2hnqcEO2LiMhGVR88L/ydpYOJcicT5U5mvTNV9eS1Y9PHA5SBCKMSGVXXiMIn4QLgTne/291LwGXAxekD3P0qT0a2he8De9aqMGpT02CG0NspAsZwn4ZltpwXEZH1F9rI1B8HFcLkxAZUyFkYaNfmGaSmNkifJ9fKbjBjgK1/imk3cF9q/RDwxEWOfwXw5bUqjIKa2PDgBCPFvngtDwzW05YiIrLhjFVh6zz1Db0FgOo8A+153AnKmK5AX2djZyEzUlFSFuuf1syQmV2TWr/U3S892YuY2UuBg8BTVq1kTRTUpITAJnyYhwernKho4D0RkY0qnwpa0sFJtZo0/J27L9mQnDpZhb58mLp4NjX3U2bzNWszs8+oux9cYN9hYG9qfU+8rYGZPQN4M/AUd59d/SIGCmpiI8V+hvuOMzyY/J+QAybQ3E8iIhtTP1CpJjNs16uSgiRASa0CkTtmxhTQHYWAJoqgRCfgRG7kLFw3i1pQ/fRD4CwzO5MQzFwC/GZDmcweD/wd8Gx3P7KWhVFD4ZSRyS2MFPvjJUd3NY+CGhGRjWk6slpAA2GsGTMoFEbmVD3NVENcMzubA4yJci+5+Bcwl4NCpYR7mTuPhqBA3/zL4+4V4LXAlcBtwCfc/RYze6uZvSA+7K8IMegnzex6M7t8rcqjTE2DI8A2YAzYyX2THagFvIjIxjRJD1PVKr2UOUYHp5HUajT/tDm9BahitUAmsVBrLgAAF7JJREFU+W4vl+HECRgaglIZDgyFvZ1Z/HVcRhfsNXla9yuAK5q2vSX1+BnrVRZlahoMEWpWdwKTnDF4G4rXRUQ2ph02xSxdHGOAPnIU4kDkL8+4quE4MyjFHT8Oj4df/Z5C+NvZCcPDyXGdjJX0D9ksy2IsuobSMV4fk5wGXN2qwoiIyCLyedhDcc6YM3X1DR0FcIP9O8K2wpxu4IBHdOULgDM5sValXks+34uwqShTM0cSpUd05jtaWhIREVlYtdo4Tg00D6qXHkQv/M3FD9P9m2rnm3FobAfucN/UjjUt+1ppxTQJG4mCmtjw4ARzKyTH2adclojIhnR/U4bmULW5B1Syvz5jdxSvVyqNAVC1CmFc4SlCAHRsjUsva0E/2SlDA1O1x5UKHB7bs9kzeSIiG1g/h6qp7k+cAOYbYRimKtDTYVQjxwzK1VxDliefByrwiO0hqDl7OKNf/pv8R0uZmpRSqR6x53JVhgeP8EBVVVAiIhvR3kIylpgDs+wtNFZFmTXNxu1QquYwg2ouX9vXfIyZ+r1mlYKa2Eixn87O8GHO5+HYZJ6hnhI785orQURko9pbGGdvocjewgxQT1RYPHklhO/13kIIfSq1Y0pMxF/v6eTGXWPxtnUo+6pzsGj1lyxR9VMsB4yO9wGT8ZZeRqcjqorXRUQyxT3kb5q/vSvlMGUxOCPF7Zw+FJocjMfBTUcFfm5H2C/ZpExNrB6M9gN9JDH9kepAq4okIiKLOH3gQwvuy83z65bPwyzJTNxjjNNNFMFAAfpyUM53xeele01ljPvqLxmiTE3NBPVgBmAGi6oaek9EZIPatuWp3D9e79bd2L27TC0vk1RD5aA/F47ZPxQ2TlS76lUsFkYYLhQy91tel9VyrxJlamLDgxCqnirACXb0lzkyEbGvcKK1BRMRkXnlrBNg/ga/TT9vUQTVCkxPJ1ssdUzIzITEhFEuOzOR/kmbRcrUxEaK/ezon8BsFigwOg5QYErthEVEMqhEmPamzoAHp5z9vda0tf74eKnAcF+FQpTNlEcLZuneUJSpieUoc3Sij9HxvrjBcAdQpju/1JkiItJK8zUBOZNravsAylXIF2B/rSFw8uPvDedNTHdg5qmMjmSJgprY2buOAUXCB73E2cPHOXfXz8hltK2YiMjmsLVhPJrk7yv3HWW8HAKWUhm6Oz2VlJk7fUJy3hk7pghzKK1P6VedGgoLwG0P7eLs4YcwC2MdjI/D3z7qao5XulpcMhERWch5p1/NrfefPe9v70ABLOd0dUJpFjq65nb1do/I5UKQY+YYTg64bzqDcz85bPbeLQpqUm4f2Zlaq/CzyjayG66LiLS/fL6vIaBJz+d017Gh1JHOmZ3HCCHNfN/rHp9n/Hh0mFz8D1zJFgU1DcqEl8TZ23cMcAY53uIyiYjIYg5VYXeuPi1CvWv3cWBLfNSx1GSX9QkuobFLeKVqQJVtjK/zXZw6wzd9Q2EFNbHhwQlGin3xmnHfZAf7CgpoREQ2vj6iaDJMSply2hagFpzkgWpqb/jxz+WquNd7SUWRA3nG6Fm74sqaUUPhlKGBSUJkf4Shgc5WF0dERJbhnP7nzqlQcoeK56l4OV7S/4avHz1biRq2d3QYMMKB4Yx2f9rkDYUV1MRGiv1UqzA00MHQQB+j480DOYmIyEa0f/g9c3qqhu/vEmFU4S7CBAm1vSRVUFXCwHzJwHszERwYzvCX/yYPalT9lDI21Ue9EVk1i++niMimlM83fl+7w2AHwGxoLwPUx9NLf7E7s15I1UyV176wsmaUqVmQMVnOcW91W6sLIiIiK9A808F82feH9UAIZKpAhd4ChCkTMpitSbp0r/aSIQpqGtRHmXwmn2GEraAu3SIiG96BnbfWqpEAPnLXWbWgpvbNXv9PrIw79Bagt+C1gAacSkXf/VmkoKZBPTL/Ki/m8rFH0NhaXkRENqKurkHu9y0cqoblOfvvp5BqYGHJf6qOu8VLJ3ePhF6v7knoE4KZckbn/TP3VV+yREHNHOmUYy+hG6CIiGx0H79vPzkiOinxt/c9tra99q3uUPHws5dURW3Zko5e6rN1d3ZksPpJ1FB4cTkyV6EoIrJJ7eYi3n9oH1Dv8jFfnsHMa3+ny9DZ5aQH5DODiWq2MhQ1GcusrDZlamK9tUfJByJiJ99BmRoRkWz4nxe9karnqEQ5oiQjAw3RTUfBqVSM6XIH0+VO5huQr/FRlqxBd+6MBUkKamIH991FvStfxEX77uFJ28rAVAtLJSIiK5FUL0WAW7zE2Zhc7Zevwp6+Sur4evXTVo2/mkmqfopt4Tou2peOSIvAduBYi0okIiIrVZ/YsrFtzHj8b9cz+k+kjksG6atXS2WSk7nMympTpib21O0hsAnzhFzHxdvvxszZp7BPRCQz3n/+n3KiBBNlKJZhstLJeLw+USbu9VRgYgKOj8M949vJN9RROTmgkFN7yixq2U+2mb0OeA2hQvNL7v5f4+1/DLwi3v56d79yvcr01O0Ad8ZrEYr5RESyZV//TqKoG3L1KiizroZjoghGcttxh6MTTl/UycBAfabuML7Nj9a97Ktik8diLQlqzOxpwMXAY9191sx2xtvPAy4BHgmcDnzNzA64ewsGiwkBzSbP5ImIZM6xyW6GBmYavr/Townn817bNjwIA/GE3Onj3/vw+9ahpKsva+PKrLZWpSJ+D3iHu88CuPuRePvFwGXuPuvu9xDSJhesf/Hqn/57q4Pr//QiIrJiX37W7wGLT0oc2tqEA257KCzVKpRKcNvITk1onFGtCmoOAE82s6vN7Ftm9oR4+24gHR4firfNYWavMrNrzOyakZGRNSxqaQ2vLSIiq21//2m1qiQAS3XQrjcCDn/djYidRAzz49Fd3DG2C3hofQu8mjZ5l+41q34ys68BD5tn15vj590OPAl4AvAJM9t/Mtd390uBSwEOHjy4iq9686VmVu/SIiKybhp7MaUflzDrJBlo75G7jixyrGTJmmVq3P0Z7v6oeZbPEzIwn/HgB4SmTUPAYWBv6jJ74m3rqD6qZLB9fZ9eRERO2RVPeQuhv0nyXV6vT7rs7E9SqST7JrGoGO8J01L/rz2Xr2dRV48Dka/+kiGtqn76HPA0ADM7AHQCo8DlwCVm1mVmZwJnAT9oTREdqHA2b2jN04uIyIoNdvTTW4iYm213ogg+/ajL+Pg5/4/Lzv4cz8t/gzD4ag7I8QeHnr/u5V0dGlG4VV26PwR8yMxuJjRaebmHKVJvMbNPALcCFeA169/zKXkDc0CJA/t/a32fXkREVsXzdzyejz90c7xWnwIniiAfz4BjBk8chF/ZekUriiirrCVBjbuXgJcusO/twNvXt0Rp6WnQupY4VkRENqrXPfrlfHPkP/FQNECoEIj4+DmXEUWhp1M+D5UKDA7Wx6jJfK+njGVWVpvGy12Qo8H3RESy7dz+fbz39L8DGgfXK8S/foUCjI7C8HD9nLYIbjYp/WovaHNHuyIi7eB/XPBnc7YlyYzJSfjUvTA01Lg/0wHNJm9To6CmQURS9bSb4lIHi4hIBpilp0yAzs7wt78ffn1f4z7JNgU1DXIkU88fRiMJi4i0g61b3z8n8ZCxBMTyqEu3gpqF5VtdABERWQW9fS/kqhNnUqm0uiRrzcGj1V8yREHNvCy1iIhI1r3ojL/kT+5/Tm1dVU3tSb2f5tBQ2SIi7aav7yL+at8lrS7G2mvLerXlU6ZmjnSGRqG8iEi7GN75s1YXQdaYMjUNNJmliEi76ujooD7AahtKGgpvYsrUNEhnZiJgulUFERGRNXD67nWeI3m9aZwaaVQPbPao9klERCQzFNTMkUSlRQ5lK0AVEZHNbpNnatSmZo4JwhT0W1FDYRERkexQUDPHQPw3W9GpiIhsdtnLrKw2BTXz2twfChERySAnTEG+ialNjYiIiLQFZWoWFKH5n0REJFM2efWTMjVzJB+ICqqGEhERyQ4FNQ0mCIFMGehqcVlERGQt7Nx1bauLsHbUpVvq+gnduDtRlkZEpD0VCqe1ugiyRpSpqdkS/3XqAY0CGxGRdnT67vtbXYQ14GHup9VeMkRBTczmPHLgeEvKIiIia2/7jm+0ugiry8E9WvUlSxTUxP7d/n+NHyVRablVRRERkXXQ3X1Om2ZsNi8FNSnP338TUATGgEmeNPTmFpdIRETWWghsfmuebRm0yauf1FC4ya88/Ds8eOJddOZPZ8fAK1pdHBERWQen734H8I5WF0NOkYKaJp2FvTx8x/9udTFEREROXsa6YK82BTUiIiLtwF1zP7W6ACIiIiKrQZkaERGRdrHJq5+UqREREZG2oEyNiIhIm/BN3qZGQY2IiEhbyN4ElKtN1U8iIiLSFpSpERERaQdO5kYAXm3K1IiIiEhbUKZGRESkXWRsVu3VpkyNiIiItAVlakRERNqAA77J29QoqBEREWkH7qp+anUBREREJLvM7NlmdruZ3Wlmb5pnf5eZfTzef7WZnbFWZVFQIyIi0iY88lVfFmNmeeB9wHOA84CXmNl5TYe9Ahhz958D3gO8cw1uHVBQIyIiIit3AXCnu9/t7iXgMuDipmMuBj4SP/4U8MtmZmtRGLWpERERaRfr36ZmN3Bfav0Q8MSFjnH3ipmdAHYAo6tdmLYIaq699tpRM/vZKl5yiDV4sVusHe8J2vO+2vGeoD3vqx3vCdrzvlp1T/vW64nGGbvya/6poTW4dLeZXZNav9TdL12D5zllbRHUuPvwal7PzK5x94Orec1Wa8d7gva8r3a8J2jP+2rHe4L2vK92vKdm7v7sFjztYWBvan1PvG2+Yw6ZWQHYAhxdi8KoTY2IiIis1A+Bs8zsTDPrBC4BLm865nLg5fHjXwO+4b4204m3RaZGRERE1l/cRua1wJVAHviQu99iZm8FrnH3y4F/BP6vmd0JHCMEPmtCQc38NmRd4Slqx3uC9ryvdrwnaM/7asd7gva8r3a8pw3B3a8Armja9pbU4xng19ejLLZGGSARERGRdaU2NSIiItIWFNSkLDXUc5aY2U/N7CYzuz7pimdm283sq2Z2R/x3W6vLuRgz+5CZHTGzm1Pb5r0HC94bv3c3mtn5rSv54ha4rz83s8Px+3W9mT03te+P4/u63cx+pTWlXpyZ7TWzq8zsVjO7xcx+P96e6fdrkfvK7PtlZt1m9gMzuyG+p7+It58ZD2F/ZzykfWe8fd2GuD8Vi9zXP5nZPan36nHx9kx8BuUkubuWUAWXB+4C9gOdwA3Aea0u1yncz0+BoaZt7wLeFD9+E/DOVpdziXu4CDgfuHmpewCeC3wZMOBJwNWtLv9J3tefA2+Y59jz4s9iF3Bm/BnNt/oe5innacD58eMB4Cdx2TP9fi1yX5l9v+LXvD9+3AFcHb8HnwAuibd/EPi9+PGrgQ/Gjy8BPt7qezjJ+/on4NfmOT4Tn0EtJ7coU1O3nKGesy49VPVHgBe2sCxLcvdvE1rKpy10DxcD/8eD7wNbzey09SnpyVngvhZyMXCZu8+6+z3AnYTP6obi7g+4+3Xx43HgNsIoopl+vxa5r4Vs+Pcrfs0n4tWOeHHg6YQh7GHue7UuQ9yfikXuayGZ+AzKyVFQUzffUM+LfXltdA58xcyuNbNXxdt2ufsD8eMHgV2tKdopWege2uH9e22cBv9Qqmowc/cVV088nvAv5bZ5v5ruCzL8fplZ3syuB44AXyVklI67eyU+JF3uhiHugWSI+w2n+b7cPXmv3h6/V+8xs654WybeKzk5Cmra14Xufj5h5tTXmNlF6Z3u7iz+r5gNrx3uIeUDwCOAxwEPAH/T2uKsjJn1A58G/sDdi+l9WX6/5rmvTL9f7l5198cRRn+9ADinxUVaFc33ZWaPAv6YcH9PALYD/62FRZQ1pqCmbjlDPWeGux+O/x4BPkv44nooSa/Gf4+0roQrttA9ZPr9c/eH4i/kCPh76lUWmbkvM+sg/PD/s7t/Jt6c+fdrvvtqh/cLwN2PA1cBv0CofknGLkuXu3ZPtsZD3K+W1H09O65CdHefBT5MRt8rWR4FNXXLGeo5E8ysz8wGksfAs4CbaRyq+uXA51tTwlOy0D1cDvxW3KPhScCJVLXHhtdUl/8iwvsF4b4uiXugnAmcBfxgvcu3lLiNxT8Ct7n7u1O7Mv1+LXRfWX6/zGzYzLbGj3uAZxLaCl1FGMIe5r5X6zLE/alY4L5+nAqqjdBOKP1ebfjPoJykVrdU3kgLoTX8Twj1y29udXlO4T72E3pg3ADcktwLoR7868AdwNeA7a0u6xL38TFCar9MqO9+xUL3QOjB8L74vbsJONjq8p/kff3fuNw3Er5sT0sd/+b4vm4HntPq8i9wTxcSqpZuBK6Pl+dm/f1a5L4y+34BjwF+FJf9ZuAt8fb9hADsTuCTQFe8vTtevzPev7/V93CS9/WN+L26Gfgo9R5SmfgMajm5RSMKi4iISFtQ9ZOIiIi0BQU1IiIi0hYU1IiIiEhbUFAjIiIibUFBjYiIiLQFBTUiK2Rm1XjW35vN7AupMTJON7NPLeP8iQW2v9DMzlvi3OvN7LKVlXx1LPc+l3mtZFb5g/Pse6qZffEUrn2VmU3Md20RaS8KakRWbtrdH+fujyJMUPkaAHe/391/bfFTF/VCwmzP8zKzcwmzyj85HlyxJVbhPps9zd2vWcXrAeDuTwNW/boisvEoqBFZHd8jngzPzM4ws5vjx71m9gkzu9XMPmtmV6czBmb2djO7wcy+b2a7zOwXgRcAfxVnYx4xz3O9hDD421dIzSRvZq+Pn+fGJItjZv1m9uE4C3Kjmf1qvP1ZZvY9M7vOzD4Zz22UZEz+It5+k5mdE29/Slye683sR2Y20HSf3ann+ZGZPS3e/ttm9hkz+xczu8PM3rWcF9PMnm1mPzaz64AXp7b3WZhA8gfx81y8nNdZRDYHBTUip8jM8sAvM/+0Gq8Gxtz9POC/Az+f2tcHfN/dHwt8G3ilu/9bfJ03xlmgu+a55m8AlxFGJn5JavubgMe7+2OA/xxv+++E4d8fHW//hpkNAX8KPMPDpKfXAP8ldZ3RePsHgDfE294AvMbDZIFPBqabyvQawpyVj47L9BEz6473PS4u86OB3zCzvSwiPu/vgefHr9fDUrvfTBim/wLgaYTgr4/FX2cR2SQU1IisXI+ZXQ88COwCvjrPMRcSAhDc/WbCEO6JEpC0FbkWOGOpJ4yzD6Pufi9h+oHHm9n2ePeNwD+b2UuBSrztGYSh4InLMAY8iVC99d24/C8H9qWeJpmMMl2m7wLvNrPXA1vdvUKjCwlD0OPuPwZ+BhyI933d3U+4+wxwa9Nzzecc4B53v8PDkOcfTe17FvCmuNzfJAzh/3AWf51FZJNQUCOyctNx5mIfYR6Z15zk+WWvz1NSBQqLHRx7CXCOmf2UMGfNIPCr8b7nEQKY84EfWn3G5WYGfDXOBD3O3c9z91ek9s82l8nd3wH8DtBDCIbOWc4NNl2v4ZorZMCvpsr+cHe/7RSuJyJtREGNyCly9yng9cAfzRNIfBf49wBxj6ZHL+OS48BA80Yzy8XXerS7n+HuZxDa1Lwk3rfX3a8C/huwBegnZI9ek7rGNuD7wC+Z2c/F2/rM7ACLMLNHuPtN7v5Owoz2zUHNvwL/IT72ACF7cvsy7nU+PwbOSLUnSlexXQm8Lp5xGTN7fLx9Ja+ziLQZBTUiq8Ddk9mBX9K06/3AsJndCvwPwqzpJ5a43GXAG+OGsOmGwk8GDrv7/alt3yZUJe0GPmpmNxFmKn6vux+Pn3Nb3O38BkIPoxHgt4GPmdmNhEbOS2Ve/iC+xo2E2cW/PM995uLn/zjw2+4+23yR5YirqV4FfCluKHwktfttQAdwo5ndEq8nz3+yr7OItBnN0i2yhuJGxB3uPhMHKF8Dznb3UouLtqHE1WkH3X10hecv+jqb2TeBN6xFl3ER2ThOpW5bRJbWC1xlZh2E9iCvVkAzrxHg62b2ihUGHgu+zmZ2FbCfkGESkTamTI2IiIi0BbWpERERkbagoEZERETagoIaERERaQsKakRERKQtKKgRERGRtqCgRkRERNrC/wcIorXzIkTzBAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f44441c8a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Fig, Ax = plt.subplots(1,1,figsize=(9.5,8))\n",
    "Sc = Ax.scatter(dm_hp_ra, dm_hp_dec, c=pz_frac_map, cmap=plt.cm.viridis, s=6, vmin=0, vmax=1)\n",
    "\n",
    "Ax.set_xlabel('Right Ascension [deg]')\n",
    "Ax.set_ylabel('Declination [deg]')\n",
    "Ax.set_title('{0}'.format(FIELD))\n",
    "CB = Fig.colorbar(Sc)\n",
    "CB.set_label('Fraction with photo-z estimate')\n",
    "Fig.savefig('plots/dmu24_{0}_sf_binary_map.png'.format(FIELD), \n",
    "            format='png', bbox_inches='tight', dpi=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add the binary photo-$z$ selection function to output catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "pz_depth_map.add_column(Column(name='pz_fraction', data=pz_frac_map))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## IV - Magnitude dependent photo-z selection functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The binary selection function gives a broad illustration of where photo-$z$s are available in the field (given the availability of optical datasets etc.). However, the fraction of sources that have an estimate available will depend on the brightness of a given source in the bands used for photo-$z$s.\n",
    "Furthermore, the quality of those photo-$z$ is also highly dependent on the depth, wavelength coverage and sampling of the optical data in that region.\n",
    "\n",
    "To calculate the likelihood of a given source having a photo-$z$ that passes the defined quality selection or be able to select samples of homogeneous photo-$z$ quality, we therefore need to estimate the magnitude (and spatially) dependent selection function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the photo-$z$ quality criteria\n",
    "\n",
    "A key stage in the photo-$z$ estimation methodology is the explicit calibration of the redshift posteriors as a function of magnitude. The benefit of this approach is that by making a cut based on the width of redshift posterior, $P(z)$, we can select sources with a desired estimated redshift precision.\n",
    "Making this cut based on the full $P(z)$ is impractical. However the main photo-$z$ catalog contains information about the width of the primary and secondary peaks above the 80% highest probability density (HPD) credible interval, we can use this information to determine our redshift quality criteria."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parse columns to select the available magnitudes within the masterlist:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20 magnitude columns present in the masterlist.\n"
     ]
    }
   ],
   "source": [
    "filters = [col for col in merged.colnames if col.startswith('m_')]\n",
    "print('{0} magnitude columns present in the masterlist.'.format(len(filters)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaled_photoz_error = (0.5*(merged['z1_max']- merged['z1_min'])) / (1 + merged['z1_median'])\n",
    "\n",
    "photoz_quality_cut = (scaled_photoz_error < 0.2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the magnitude dependent selection function in a given masterlist filter, for each of the Healpix clusters we do the following:\n",
    "1. Find the number of masterlist sources within that cluster that have a measurement in the corresponding filter. (If this is zero - see stage 3B)\n",
    "2. Calculate the fraction of sources at that magnitude that haThis relation typically follows a form of a sigmoid function the declines towards fainter magnitudes - however depending on the selection being applied it may not start at 1. Similarly, the rate of decline and the turnover point depends on the depth of the optical selection properties of that cluster.\n",
    "3. \n",
    "    1. Fit the magnitude dependence using the generalised logistic function (GLF, or Richards' function). Provided with conservative boundary conditions and plausible starting conditions based on easily estimated properties (i.e. the typical magnitude in the cluster and the maximum point), this function is able to describe well almost the full range of measured selection functions.\n",
    "    2. If no masterlist sources in the cluster have an observation in the filter - all parameters set to zero (with the GLF then returning zero for all magnitudes).\n",
    "4. Map the parameters estimated for a given healpix cluster back to the healpix belonging to that cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m_ap_vista_j\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "75ea9b28bd1b49568af284c86ae863f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_vista_j\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1ffb16840fae47648f6d032c5fd97a72",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_vista_h\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "deda9b0876eb49919d0a8ce7b97c1fac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_vista_h\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "160bdbd9351545b3a57738b412eea70c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_vista_ks\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6dd013d30ede4740a05edc98cdb1c3e8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_vista_ks\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eb3b5c05736d4f27a300393264f0e1c2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_irac_i1\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4c1a9be939994f589ea3b25989ac0259",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_irac_i1\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f6fce6ba4aae4f778afd97de77282801",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_irac_i2\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3d763cdd40994bb6ab9aab0e46cf6c5a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_irac_i2\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be5ea4dff2d544d3bdca63033c0354ee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_decam_g\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d13c4818bc594e2dae6299fa17051b6f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_decam_g\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "335157ca044248f29d7aece09b31e56a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_decam_r\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b21a8b705a614cf4ba10992ae66f4238",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_decam_r\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "82621199d9b542b9a395a1675789f2bd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_decam_i\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "739c55b1b66b4551841e00b6e952c4f0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_decam_i\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0489fbc01d9a4915b315451343bfdffe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_decam_z\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "392e244ac7104af294f4ab82fd04869f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_decam_z\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d34559d0bfaf4c7881af4f1abecf6689",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_ap_decam_y\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b7f82cad79284e4e89856ee2f4455f02",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "m_decam_y\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "34a492d5957944f3ab5a762b43714e48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FloatProgress(value=0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "for photometry_band in filters:\n",
    "    print(photometry_band)\n",
    "    \n",
    "    pz_frac_cat = np.zeros(len(merged))\n",
    "    pz_M_map = np.zeros((len(dm_hp_ra),6))\n",
    "\n",
    "    if np.isnan(merged[photometry_band]).sum() < len(merged):\n",
    "        m001, m999 = np.nanpercentile(merged[photometry_band], [0.1, 99.9])\n",
    "\n",
    "\n",
    "        counts, binedges = np.histogram(merged[photometry_band], \n",
    "                                        range=(np.minimum(m001, 17.), np.minimum(m999, 29.)),\n",
    "                                       bins=10)\n",
    "\n",
    "        binmids = 0.5*(binedges[:-1] + binedges[1:])\n",
    "\n",
    "\n",
    "        with ProgressBar(NCLUSTERS, ipython_widget=True) as bar:\n",
    "            for ic, cluster in enumerate(np.arange(NCLUSTERS)[:]):\n",
    "                ml_sources = (merged['hp_depth_cluster'] == cluster)\n",
    "                has_photoz = (merged['z1_median'] > -90.) * photoz_quality_cut\n",
    "                has_mag = (merged[photometry_band] > -90.)\n",
    "\n",
    "                in_ml = np.float(ml_sources.sum())\n",
    "                withz = (has_photoz)\n",
    "\n",
    "                frac = []\n",
    "                frac_upper = []\n",
    "                frac_lower = []\n",
    "\n",
    "                iqr25_mag = (np.nanpercentile(merged[photometry_band][ml_sources*has_photoz], 25))\n",
    "\n",
    "                if (ml_sources*has_photoz*has_mag).sum() > 1:\n",
    "                    for i in np.arange(len(binedges[:-1])):\n",
    "                        mag_cut = np.logical_and(merged[photometry_band] >= binedges[i],\n",
    "                                                 merged[photometry_band] < binedges[i+1])\n",
    "\n",
    "                        if (ml_sources * mag_cut).sum() > 0:\n",
    "                            pass_cut = np.sum(ml_sources * withz * mag_cut)\n",
    "                            total_cut = np.sum(ml_sources * mag_cut)\n",
    "                            frac.append(np.float(pass_cut) / total_cut)\n",
    "\n",
    "                            lower, upper = binom_conf_interval(pass_cut, total_cut)\n",
    "                            frac_lower.append(lower)\n",
    "                            frac_upper.append(upper)\n",
    "\n",
    "                        else:\n",
    "                            frac.append(0.)\n",
    "                            frac_lower.append(0.)\n",
    "                            frac_upper.append(1.)\n",
    "\n",
    "                    frac = np.array(frac)\n",
    "                    frac_upper = np.array(frac_upper)\n",
    "                    frac_lower = np.array(frac_lower)\n",
    "\n",
    "\n",
    "                    model = GLF1D(A=np.median(frac[:5]), K=0., B=0.9, Q=1., nu=0.4, M=iqr25_mag, \n",
    "                                  bounds={'A': (0,1), 'K': (0,1), 'B': (0., 5.),\n",
    "                                          'M': (np.minimum(m001, 17.), np.minimum(m999, 29.)),\n",
    "                                          'Q': (0., 10.),\n",
    "                                          'nu': (0, None)})\n",
    "\n",
    "                    fit = LevMarLSQFitter()\n",
    "                    m = fit(model, x=binmids, y=frac, maxiter=1000,\n",
    "                            weights=1/(0.5*((frac_upper-frac) + (frac-frac_lower))), \n",
    "                            estimate_jacobian=False)\n",
    "                    parameters = np.copy(m.parameters)\n",
    "\n",
    "                else:\n",
    "                    frac = np.zeros(len(binmids))\n",
    "                    frac_upper = np.zeros(len(binmids))\n",
    "                    frac_lower = np.zeros(len(binmids))\n",
    "\n",
    "                    parameters = np.zeros(6)\n",
    "\n",
    "                # Map parameters to cluster\n",
    "\n",
    "                # Map parameters back to depth map healpix \n",
    "                where_map = (km.labels_ == cluster)\n",
    "                pz_M_map[where_map] = parameters\n",
    "\n",
    "                bar.update()\n",
    "        \n",
    "        \n",
    "    else:\n",
    "        pz_M_map = [np.zeros(6)] * len(dm_hp_ra)\n",
    "        \n",
    "    c = Column(data=pz_M_map, name='pz_glf_{0}'.format(photometry_band), shape=(1,6))\n",
    "\n",
    "    try:\n",
    "        pz_depth_map.add_column(c)\n",
    "    except:\n",
    "        pz_depth_map.replace_column('pz_glf_{0}'.format(photometry_band), c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The selection function catalog consists of a set of parameters for the generalised logistic function (GLF, or Richards' function) that can be used to calculate the fraction of masterlist sources that have a photo-$z$ estimate satisfying the quality cut as a function of a given magnitude. e.g.\n",
    "\n",
    "$S = \\rm{GLF}(M_{f}, \\textbf{P}_{\\rm{Healpix}})$,\n",
    "\n",
    "where $S$ is the success fraction for a given magnitude $M_{f}$ in a given filter, $f$, and $\\textbf{P}_{\\rm{Healpix}}$ corresponds to the set of 6 parameters fit for that healpix.\n",
    "\n",
    "\n",
    "In practical terms, using the GLF function defined in this notebook this would be `S = GLF1D(*P)(M)`. Similarly, to estimate the magnitude corresponding to a desired photo-$z$ completeness one can use the same parameters and the corresponding inverse function: `M = InverseGLF1D(*P)(S)`.\n",
    "\n",
    "### Save the photo-$z$ selection function catalog:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "pz_depth_map.write('{0}/photo-z_selection_{1}_{2}.fits'.format(OUT_DIR, FIELD, SUFFIX).lower(), format='fits', overwrite=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ![HELP LOGO](https://avatars1.githubusercontent.com/u/7880370?s=75&v=4)\n",
    " \n",
    "**Author**: [Kenneth Duncan](http://dunkenj.github.io)\n",
    "\n",
    "The Herschel Extragalactic Legacy Project, ([HELP](http://herschel.sussex.ac.uk/)), is a [European Commission Research Executive Agency](https://ec.europa.eu/info/departments/research-executive-agency_en)\n",
    "funded project under the SP1-Cooperation, Collaborative project, Small or medium-scale focused research project, FP7-SPACE-2013-1 scheme, Grant Agreement\n",
    "Number 607254.\n",
    "\n",
    "[Acknowledgements](http://herschel.sussex.ac.uk/acknowledgements)\n",
    "\n"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}