{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "28f62e7e",
   "metadata": {},
   "source": [
    "# How to Use DeepDISC-Euclid Catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e76fe17b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "from astropy.visualization import ImageNormalize,AsinhStretch\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0591138b",
   "metadata": {},
   "source": [
    "## 1. Basic Information\n",
    "\n",
    "We train three models using the Euclid Q1 Data Release images covering the Euclid Deep Field North (EDF-N). Specifically, we use the background-subtracted mosaic images in 9 bands: 4 Euclid bands ($I_{\\mathrm{E}}$, $Y_{\\mathrm{E}}$, $J_{\\mathrm{E}}$, $H_{\\mathrm{E}}$) and 5 UNIONS bands ($ugriz$) . The models are designed for the following tasks:\n",
    "\n",
    "* Model 1: Source detection and classification\n",
    "* Model 2: Galaxy photometric redshift (photo-$z$) estimation\n",
    "* Model 3: Quasar photometric redshift (photo-$z$) estimation\n",
    "\n",
    "By combining the predictions from these models, we construct the final catalog for EDF-N. \n",
    "\n",
    "Firstly, let's read the catalog. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "69204579",
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = '/data8/yzhe/euclid/Q1/catalog/deepdisc' # replace with your path to the deepdisc catalog\n",
    "tcat = Table.read(f'{filepath}/cat_final.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edea8c9a",
   "metadata": {},
   "source": [
    "You can check the column names here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "43c440fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 pred_id\n",
      "2 tile_id\n",
      "3 ra_bbox\n",
      "4 dec_bbox\n",
      "5 type_prob\n",
      "6 type_pred\n",
      "7 gmm\n",
      "8 z_med\n",
      "9 z_p16\n",
      "10 z_p84\n",
      "11 z_mode\n",
      "12 z_hpd_lo\n",
      "13 z_hpd_hi\n",
      "14 mag_vis\n",
      "15 mag_y\n",
      "16 mag_j\n",
      "17 mag_h\n",
      "18 mag_u\n",
      "19 mag_g\n",
      "20 mag_r\n",
      "21 mag_i\n",
      "22 mag_z\n",
      "23 nsave\n",
      "24 ng\n",
      "25 ni\n",
      "26 nj\n"
     ]
    }
   ],
   "source": [
    "for k,col in enumerate(tcat.colnames):\n",
    "    print(k+1,col)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad17b5cc",
   "metadata": {},
   "source": [
    "* `type_pred` gives the object type with the highest predicted probability from Model 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b4547b5f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of objects: 13355031\n",
      "Number of Stars 2504839\n",
      "Number of Galaxies 10436713\n",
      "Number of Quasars 413479\n"
     ]
    }
   ],
   "source": [
    "print('Number of objects:', len(tcat))\n",
    "print('Number of Stars', len(tcat[tcat['type_pred']==0]))\n",
    "print('Number of Galaxies', len(tcat[tcat['type_pred']==1]))\n",
    "print('Number of Quasars', len(tcat[tcat['type_pred']==2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d4f83c0",
   "metadata": {},
   "source": [
    "* If you prefer not to assign the object type based on the maximum probability and instead apply a custom threshold, you can use `type_prob`, which is a three-element array corresponding to the probabilities of the source being a star, galaxy, or quasar, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "04af62be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "[0.00197251 0.9189384  0.07908909]\n"
     ]
    }
   ],
   "source": [
    "idx = 10000\n",
    "print(tcat['type_pred'][idx])\n",
    "print(tcat['type_prob'][idx])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dd946e1",
   "metadata": {},
   "source": [
    "* The catalog also provides photometric redshift information. Note that stars and objects not matched between Model 1 and Models 2/3 do not have redshift information. The median, 16th, and 84th percentiles of the photo-$z$ PDF can be obtained as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "68715309",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Row index=10000</i>\n",
       "<table id=\"table140045671179056\">\n",
       "<thead><tr><th>z_med</th><th>z_p16</th><th>z_p84</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>1.2597891380441195</td><td>1.007827016582181</td><td>1.4751298966306639</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Row index=10000>\n",
       "      z_med              z_p16             z_p84       \n",
       "     float64            float64           float64      \n",
       "------------------ ----------------- ------------------\n",
       "1.2597891380441195 1.007827016582181 1.4751298966306639"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tcat['z_med','z_p16','z_p84'][idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ffec8c",
   "metadata": {},
   "source": [
    "* The mode and the 68% highest posterior density (HPD) credible interval of the photo-$z$ PDF are also provided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "29f24426",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Row index=10000</i>\n",
       "<table id=\"table140042439353712\">\n",
       "<thead><tr><th>z_mode</th><th>z_hpd_lo</th><th>z_hpd_hi</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>1.2833999999998036</td><td>1.0508999999998292</td><td>1.5100999999997784</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Row index=10000>\n",
       "      z_mode            z_hpd_lo           z_hpd_hi     \n",
       "     float64            float64            float64      \n",
       "------------------ ------------------ ------------------\n",
       "1.2833999999998036 1.0508999999998292 1.5100999999997784"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tcat['z_mode','z_hpd_lo','z_hpd_hi'][idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1627304",
   "metadata": {},
   "source": [
    "* DeepDISC predicts a segmentation magnitude for each object. We sum the flux within the DeepDISC segmentation map and convert it to AB magnitudes across all EDF-N bands (four Euclid bands and five UNIONS bands). Bands with zero or negative flux are masked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "db13b8f2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Row index=10000</i>\n",
       "<table id=\"table140042455442896\">\n",
       "<thead><tr><th>z_hpd_hi</th><th>mag_vis</th><th>mag_y</th><th>mag_j</th><th>mag_h</th><th>mag_u</th><th>mag_g</th><th>mag_r</th><th>mag_i</th></tr></thead>\n",
       "<thead><tr><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>1.5100999999997784</td><td>26.009580612182617</td><td>23.97304344177246</td><td>23.07777214050293</td><td>22.656497955322266</td><td>28.2060546875</td><td>27.010385513305664</td><td>26.820449829101562</td><td>27.93269920349121</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Row index=10000>\n",
       "     z_hpd_hi           mag_vis             mag_y             mag_j             mag_h            mag_u           mag_g              mag_r              mag_i      \n",
       "     float64            float64            float64           float64           float64          float64         float64            float64            float64     \n",
       "------------------ ------------------ ----------------- ----------------- ------------------ ------------- ------------------ ------------------ -----------------\n",
       "1.5100999999997784 26.009580612182617 23.97304344177246 23.07777214050293 22.656497955322266 28.2060546875 27.010385513305664 26.820449829101562 27.93269920349121"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tcat[tcat.colnames[12:21]][idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46497bd3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "05961250",
   "metadata": {},
   "source": [
    "## 2. Retrieve Segmentation Map\n",
    "\n",
    "To store the segmentation maps efficiently, the edge coordinates (remapped onto the Euclid tile coordinate system) are concatenated into 1D arrays, and stored in several separate binary files (`../segs/segEdge_{nsave}_pairs.npy`) to save space. And the corresponding index array (`../segs/segEdge_{nsave}_index.npy`) is provided to locate the edge coordinates for each object.\n",
    "\n",
    "Columns 23–26 provide the indices (`nsave`, `ng`, `ni`, and `nj`) used to retrieve the segmentation map edge coordinates from the binary files for a given source. An example is provided below.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f1d4737",
   "metadata": {},
   "source": [
    "Firstly, let's randomly pick one object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "712035ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "102157301\n"
     ]
    }
   ],
   "source": [
    "idx = 10000\n",
    "pred_id = tcat['pred_id'][idx]\n",
    "nsave = tcat['nsave'][idx]\n",
    "ng,ni,nj = tcat['ng','ni','nj'][idx]\n",
    "tileid = tcat['tile_id'][idx]\n",
    "print(tileid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7f13bee",
   "metadata": {},
   "source": [
    "The edge coordinates can be easily retrieved using `nsave`, `ng`, `ni`, and `nj` with `load_seglist`. The returned `edge_coords` contain the $(x, y)$ positions of the segmentation map edges for the specified object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "818bf94a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from modelTool import load_seglist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0ae774d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "edgeList = load_seglist(filepath=f'{filepath}/segs/', base=f\"segEdge_{nsave}\", mmap=True, n_sublists=10)\n",
    "edgeCoords = np.array(edgeList[ng][ni][nj],dtype=np.int32)[:,::-1] # (y,x) -> (x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3c175af",
   "metadata": {},
   "source": [
    "We can display the segmentation map for this object in the image. We use the $I_{\\mathrm{E}}$ (VIS) band image at `tileid = 102157301`. Since the image is too big, we only display a cutout here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "15422031",
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu = fits.open(f\"{filepath}/demo/VIS_{tileid}.fits\")\n",
    "img = hdu[0].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "83db1434",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a cutout; and convert to local coordinates relative to the cutout\n",
    "srcx = (edgeCoords[:,0].max()+edgeCoords[:,0].min())/2\n",
    "srcy = (edgeCoords[:,1].max()+edgeCoords[:,1].min())/2\n",
    "x0,x1 = int(srcx)-50,int(srcx)+50\n",
    "y0,y1 = int(srcy)-50,int(srcy)+50\n",
    "img_cut = img[y0:y1,x0:x1]\n",
    "edgeCoordsShift = edgeCoords - [x0, y0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "85b5ea49",
   "metadata": {},
   "outputs": [],
   "source": [
    "from modelTool import plot_pred_boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "022e42dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcLZJREFUeJztvXm0FeWd7v9lEMQJFQXFACKizCogCGhsI9GgiUPsDP0z95rc9M1KLiYqt9vWm9a0iYZ0utMxo97YfY29WmMnnTZxiFEbBUdAUBAcAIEoDhiNII5oOPv3h1J53s8+561T5xzcG30+a7msovZQ9dZbu049z3foVqvVamGMMca8y3Rv9A4YY4x5f+IbkDHGmIbgG5AxxpiG4BuQMcaYhuAbkDHGmIbgG5AxxpiG4BuQMcaYhuAbkDHGmIbQs9E7QFpaWuKZZ56JXXfdNbp169bo3THGGFORWq0WL7/8cgwcODC6d2/7OafpbkDPPPNMDBo0qNG7YYwxppOsW7cuPvCBD7S5veluQLvuumun3q83r5133jnZ9rvf/S5Zf+ONN9r9ucOHD0/WV61a1eZr999//2T9ueeeK5aHDBmSbHvsscey3zt48OBi+emnn062bdmyJftehZPgqaeeavd7x4wZk6w//PDDxTL/WHjyySeTdX2KnTRpUrKNY/jiiy+2uQ/77rtv9rW6H6tXr062dabaVM+ef7pEWlpakm1cJ7169SqW33zzzexrd9xxx2K5yrzU+RER0bdv32R92bJl7f4sotfPq6++mn1tjx49Wl2OqJ97vC51H/v06ZNsO/DAA9t8baPYbbfdkvW33nqrWN5pp52SbXvvvXeyrr8Fr7/+erKtynnnOHGf9Hv22GOPZBvPzwsvvFAs77XXXm1u6whlv+dNdwPqrOymj3sc6M58Nj+rvfvA763yOfyszux/7jG4jNw4ln2uvlZ/zKvuU25MuY/c1pkbkH5W1fHPvZ7bOnpuOS5V51eOKvuUG6cq+5g7r81C7tyVHWtXXc9lY9zR13bmd6I9302a7gbUWZ544okOv3fgwIHFMk9E2ZOKsmbNmmRd/wrg09GGDRuSdZ4w/QukyhMPeemll5L1HXbYoVjWv+Ba49FHH03W9S8qPlUSfUpYsWJFsu0Pf/hD9r0Kn/74l75+Nv865F+aOfgX9+OPP97ma/nX3csvv9zmZ23evDnZpn+hRkTsueeexTKfTnPnnU9WS5YsafO1VdG/5l955ZXsa//4xz+2uY3vzT058umoM8ezzz77FMvr16/v8OdMnDgxWR8wYECyrvOY54NPJjpOvB5Iv3792vwePmnlnlT4B9jvf//7Nl+bO4/bAkfBGWOMaQi+ARljjGkI7zkJTo1uNYEjIu6+++5knTKbSjU0D+lf6OMxpRjKNirN3HXXXcm2TZs2RQ79Hu7Dxo0bs+9VKMGpNHDQQQcl2zTIIKJeVqjymP7Rj360WL7pppuSbZ3xZng8SpnkplLrM888k2zLBQswGIOybP/+/dt8L+fELrvskqyrFEvJjbKUSoyUeHg8nUHnNQMJOAeef/75Ynm//fZLtvE6vPfee5N1lcpUiuTnVkWvu6OPPjrZRvmLEp3uswaIRNQH3uj+0/OhJZCT3UaPHp2s54Jyhg0blqwvWLCgzddS5s/B8ec+dEbebg0/ARljjGkIvgEZY4xpCL4BGWOMaQjdap0R4rcBmzZtqgux3X333YvlsjDmZ599ts3PZogzQ2NzUAdWnb53797JtgMOOCBZV6+APgH1foYm50JwGaaZC68k6kEcc8wxybbf/OY37d4H+iLLly9P1lUf5/7Sa1J9uUxbpn+n6/QnRo0alaxrAmxZCLomwObmVlfCY6P3V5bUqtCTyJ1L+kk5f5J+hSaqliVU0v/S9a4MI9d0Ac4n/hZoWkJE6n995CMfSbYxVF99N4acM01h3bp1xTJ9QSYV0+fM+Z58r/5mlqWQ6Djx2uFc5O+tvl73oaWlJZ566ql46aWX6uZV8vnZPTPGGGO2Eb4BGWOMaQi+ARljjGkI20UekOa7sORElaKa9I+YY6AeBYtq5nRt+givvfZasq6aPS03+iI5n4H5FcwPyXlAuaKAN9xwQ7Ke87tImVejn8WiiCRXN4rnjtq6lnfha1euXJmsV8lj6irfhz4h54zuv+YpRdTr8Presv2rUr6JOSA5D4j+nUI/pcxnq1KAWMdx6NChyTZ6SzqmLFBLP4VjrHmAhx56aLKNeUB6Td96663JNo6hlpTiOLFEEXNuch4Qf6/0eOgZMi9L5wiLzlYpkMp9aA9+AjLGGNMQfAMyxhjTEHwDMsYY0xCa1gPq379/oWOq9pnTQcug/0JdXjXMKvkThN6MwrYC9FtyUEunt5GDrQIUatFluq96OcwLYM6KMn/+/Oznqnd25JFHJtuqtMMgzL3S/K+ylhB6rOPHj0+28dgfeOCBZF3nF88zGxxq/hTrudET2lb5SF1VR47zlD4nrw+9tnLzNCKtScc8Gn5u7rcil0cWkZ4verpspaG/GyNHjky2rV27NllXL4reDOcp55P+JvH3iLmHun3cuHHJNnpNzN3Lkcv705y/lpaWduUl+gnIGGNMQ/ANyBhjTENo2lI8e+yxR/FYrDIJH5X5OKmP+5RXKLmxFI+W2GB4d06WYmkOfq6W32AZdYZa8pGdYZHthbIa20vo8Wi5mYh6KSYXZs7xZ6h1LgSd46SyQVnYdQ6eZ0odVaa8yh6UzTgn+D0a5l8WiqwhxGXdR3Wf/uIv/iLZtnDhwmSdMu0hhxxSLC9dujTZxi6bY8eOLZZz5f5JWUv0IUOGJOvacfSXv/xlso2h1Srxcu4x5DknFRPKnHoMn/jEJ5JtlOdvvPHGNveX7Rg0jaRK6aOIiL322qvVz4molx/12mJpsyqdiKvgUjzGGGO2G3wDMsYY0xB8AzLGGNMQmtYDiviTDltlF1XHpu5epsN3FHpADP/U72W5E/oInQmFVX28SrkZhk7TL+K4UX/Oofo+fTW2QVboYXGdobFV0OOlX0TPTcNxuQ8MGS5rZ5xD523ZsammTt+A547lmzQ0tivbd3cGbR3ANvMsg6NpDJyXPJfqbZa1ecj5MQxxVm8sIuL2228vlvVY+DkR6XVJP3jNmjXJes4LZFh/lZI57yb2gIwxxjQlvgEZY4xpCL4BGWOMaQhN6wEdddRRhVavuQwsZZODJSgeeuihrtnJEpgLo54Qcxe2FfR1cp4QX8uciCpl1lnqXfMTynIxNPeHeT/Mczj88MOT9blz57Z7HztK1bwNheNCP0m9p7I8Gv0sbuOcz7WtoD/Ekiy6j1X806plrOjdKMwV6yoGDBiQrHNu6j7zetC8vqpUyU1qRlguqKzcjj0gY4wxTYlvQMYYYxpC00pw7SVXsoWdCxn2u63Csim3sERLRz+Ln8tx0lIxPFZWDs5JgSwnQqlJx5Xb2KWyLPy1ozDcVSU6yrRaIiQirWZcReLheDNkOydzTp48OVlnaRsNV1+0aFG792lbotWNOU45+YjXXWckKy0/E5GGtpfJ2RqqXDVMecSIEcUy0yooKeqcX7FiRbKNkpXOzapltlQSpUxb5TcmJ/F2pgtAa1iCM8YY05T4BmSMMaYh+AZkjDGmIWz3HlBXkmvlwHI7Wmq/SpdWatrUgdk9UrX2Ms04F+KZKwHEbTx2hqyq9s4wTLbL0OMtC+fWLqh333139rVVGD16dLK+atWqYpnaOfdfQ4Rff/317PfQa9p1112L5Ycffjj7Xn1tWVfQKlB/15D7slJBuk/0tzgW6itwvNlBlG0fdL6VtaLQ7+FPF/0L9UzpFzGsn/ukn01P8eijj07WtUMt26lwPunncu7Rey0bi3cDzh+WidJ5oR5VrVaLP/7xj/aAjDHGNCe+ARljjGkIvgEZY4xpCPaAhFyZe3LQQQcVy2y7y5wJ1Zvp47B8Pv0XPT30h9gWQXNjuA8sD5QrxU8vg/lS6h+VlRNRDZw6PMvR6PExj4nlUDrTSkDnF4+NOnyu1IjOgYj69tcnnXRSsXz99ddn90mPnT4hc0tyMD+Kn8V8sPbSlblt5IgjjiiWeS3Rp9J5TX+F80s9oarldHSfHnzwwTb3IaLjbVBIWakn3af58+dnP0uPvTO5PPTKquYu2QMyxhjTlPgGZIwxpiFsdxJcnz59knWGg6pUwJIgq1ev7oI9fJtDDjmkWNZw1Yj6ro4qJ1E2YyVgPsJrKQ+GZVYJ16XcpWGnrIbNsOtHHnkkWc+Fq+fgOLFMiZ4fdoukjMBj132ibMbz0VFGjRqVrHNcyAc+8IFimfOU8kquLBTHbciQIcUyJSpKk+ycqRJjlc62ZXKjQrmXMi2PR6Uzzv/OVMPuKmmM1yj3aerUqcUyQ7YpKeag3EU5mOHsypgxY5J1lfbvuOOOdu9DGQxX1/mlc7pWq8WWLVsswRljjGlOfAMyxhjTEHwDMsYY0xC2Ow+IUJPUUhH8nLKSOerrLF26tN37QOihaEhk1fBh3Sf6HmvWrGnzfQylZjjusmXLimWG7tJvYQiujjn1XXbg7Cr/pSvJhdvTl1KPhT5Cla6z++67b7Ku5VsiOt46oKzNQ2fo6D4xtYBe0wEHHJCsq/c3adKkZBvneBXfanuH16weO/0h+kdaKqysbYX66rkSSxH15Y/KsAdkjDGmKfENyBhjTEPwDcgYY0xD6Fn+kuYmp2+WeT7jx49P1svK0yvawoDvY7z+oYceWizTA9LWvxERjz32WLKuXpTmf0TUezfqt0yYMCHZxvyEnJZLL4M5Kvo91OinTZuWrD/wwAPFMvVlenTqnfHc5XIgIlK9nMfGXKVcmaVc6R3CcaJmr/kizCXJeZfMXyNaRqbM86Efqd9TlkeW832YF6TXA8vEVMnHY6kgehAHH3xwscwSRbnSPF1Zyok+YZU5o14gWxvwOtOW6BHpfGM+JEtXVfFq1OPlNbqtQwT8BGSMMaYh+AZkjDGmIfgGZIwxpiFsFx6Q6rfUbtmGV8l5JBGpP1EGWyGor1BWqj63j2Vl4ZUqdaWuu+66ZF3rkkWknsTEiROTbXfeeWeyzhLs+++/f7FMH0HbXUekPg/Hhfq+fg89n7JaZKqf83wQ9WqYV8L8qbI24go/S3V6+k7U//W1++23X7LtoYceavc+EI4j67C1F+Zy0HvNtQeoMsfLfNhc7UGdlxFpLgw9tzJyrb8744to/hd/U+i/cD7pvC1rq6Gfzbw+nrucV172G9pZ/ARkjDGmIfgGZIwxpiFs96V4iJZH4SM55YhcGCplEIY9drSzZGdLW+RQ6YkhqVVKqVCm4Tjp2LDNA8ftnnvuaff3agkXlU8i6qUxHp/uY1mLCA2jZXgu5a5cmH8ulJowhJ6ylH4WP6eslIrCc8dyQfpZ3FaltQbbZWho9YEHHphsq3KtMFx96NChybqGG/NaooSlc57nmaHTHIvc9cJzqdL4uHHjkm0M1dd5S+mrTKrMzZHO0FVtKzREfsuWLfH444+7FI8xxpjmxDcgY4xpMs4444yo1WpRq9Xia1/7WqN3Z5vhG5AxZrtg2rRpsWHDhvjkJz/Z4c+YNWtWPP300/HWW28V/73++uvxxBNPxG9+85v49Kc/3YV7bMrYLsKw1Wdg2XEtcxMR8eijjxbL1J5z4dCEpWvodVRBy9ozFJkeEMMetZwLNW6in616bER92HLOe6IeTlSH5/nIhXSyJAvDQ3PlURgOTa9DtfWyFgXDhw8vltlemZ+rWjtLpVTR4RlCT39S21iwJQHnou4jSzex9E5HvUrC+c9yOrpP9HGqwPOhc2Krj3PffffFxo0bk/I/EXnfZuvntOb79uzZM/r37x/9+/ePY445Jg4++OA4//zz2/wsnku9vvX3J6J+zuRaXOgc1uVu3bpF9+7dk1bnTDXgdZe7vlniR9MfOD78LaAXrq01OtLCxk9AxpjtgtGjR8fmzZvrfuQ7yv333x9/9Vd/Feeff37ceuutxb9/5jOf6ZLPN+X4BmSM2S4YM2ZMPPbYY52K1FI2btwYDz/8cCxZsiS++93vFv/OhnpHHXVU/PznP4+VK1fGhg0bYvPmzfH000/Hv//7v8fYsWOT1/7zP/9zvPnmm/Hmm2/G9OnT46KLLop169bF66+/HnfffXfd6yMijjnmmFi4cGG8+uqrsXLlyvjSl77UJce3PbBdSHDGmPc3O+ywQwwfPjyuvfbaLv/snj17xjHHHFOsU66eOnVqfOITn0j+beDAgfHJT34yTjzxxJg2bVpdZe6IiB/84AdJuPq0adPi5z//eYwZM6aQn6dMmRI333xzIVsOGzYsfvjDH3aqI/P2RNPegPbee+9CB1XdlK1o6RvoaxmTzxwPxqe/8sorxTI1VGq5OXIauPoPEfU+AnMBVHOt0gK6tQuivTDHhrk96knkWhtEpGX6mc+S0+zpYfF46HXouSvLp8jlJtGD02NnmRiOf84T4piydbmyZMmSZJ3epY4/c5Ho+dCjU4+FXsBTTz3V5j5xH7iu3gE9qzL0R5o+59br+5BDDolevXol13CV9imt8eEPfzg+/OEPJ//2hz/8Ib75zW8m/7Zw4cI488wz44knnoiXX345evbsGePHj49vf/vbsfPOO8ff/u3fxje+8Y2ISOf4oEGD4txzz41Vq1bF9773vRg8eHDsv//+cfTRR8dvfvObiIj4zne+U/xe3HbbbfG9730vDj300Pi7v/u74nNqtVq0tLQk5bToazI3SZ/iFi1alGzjnMi1OefvLT06fa/m17W0tLSrfXrT3oCMMWYrW/+Y3NZPBm+88UbdD/T8+fPjqKOOii984QsxbNiwuu2jR49u9bN+/OMfxz/8wz9ExNu1DP/+7/8+Iv6UqLv33nvHlClTiu/91Kc+FRs2bIibbropRowY8b7woip5QFu2bIkLLrgghg4dGn369Ilhw4bFN77xjeRpoVarxYUXXhj77rtv9OnTJ6ZPn15XoNIYY6pwyCGHRMSfVIy1a9fGq6++Gi+//HLxX1Xv5NZbb42PfOQj8bGPfSy+973vRUtLS+y3337xgx/8IHnC/tnPfhYXXXRRjBs3ru7mE9F2gdd58+YVy1plYutTtkaQrV69OjZs2FCsL1y4sNKxbK9UegL6+7//+7jsssviqquuitGjR8eiRYvic5/7XPTt2ze+8pWvRETEt7/97fj+978fV111VQwdOjQuuOCCOP744+ORRx5J5LEyXn311aLUhoaA8vH+pJNOStZVCmC30SOOOCJZZzilSjOcaLmuoJSW+JiqlIXu8ns1tJqyIF+rj92UeBiyqvIjX1tWkkUlRh4rJSwtl1JFmil7bS68OFf6owxW++2q6r9VyumUVQ1XyuYTQ991zlMWZAksDe3NbYtIpW+VQ9uDzj/tGKqMGzcu1q9fn5TQOe6447JyqoYqc0wj3pb3brnlloh4+wd/zJgxceyxx0afPn3ipJNOiiuuuCIGDRoUJ598ckS8/dty7rnnFnPvtttuKz6rtfOgNxSVa7eeTw2r7tGjR7K/KnF27949evToke1+mwsNp3TPEPoxY8YUy8uXL0+2cf7nqsXrHGlvibFKN6B77703Tj755DjxxBMj4u1chp/97GfF3bpWq8Wll14af/u3f1uctH/913+NAQMGxK9+9SsneRljOsS4ceO2ufymN/qtf7CpB3jLLbfE5ZdfHr17967L1eoI+uM9ZMiQ2H333Ysf/MMPP7zTn789UEmCmzp1asyZM6eIElm6dGncfffdMWPGjIh4+7F4/fr1MX369OI9ffv2jcmTJ8d9993X6mdu3rw5Nm3alPxnjDFbGTBgQAwYMKDLb0D9+/ePadOmxdFHHx1nnXVWHH300cW2rb9x+mTxoQ99KD796U/HaaedFldeeWWnv//5558vAgT69OkT//Zv/xYzZsyIc889t1PVHrYnKj0BnXfeebFp06YYMWJE9OjRI7Zs2RKXXHJJnH766RHxpwx5RigNGDAgyZ5XZs+eHRdddFFH9t0Y8z6A/s9Wbr755qSaxsknn1zXTDHHCSecECeccELdvy9dujRuuOGGiHi7gdyNN94YH/3oR2PPPfeMn/3sZxHxthrEiuAd4aKLLopf/vKX0atXryQqb9WqVXURs+9FKt2Afv7zn8fVV18d11xzTYwePTqWLFkSZ599dgwcODDOOOOMDu3A+eefH7NmzSrWN23aFIMGDYq99tqr0KjVV+CjKW92H/zgB4tllg+hb/DII48k66q7qn4bUW806mur6PsMTWQEDUO4c14HO5UyXFqhh6Uh0PzOsjYDej7Y4iLnmVCLnjBhQrKuGndZ2SGi3SKpPx933HHJuuZ5MJSabSByfh7pTKsN1c+pAlT5HM1niYi44447knUtAaTdOSOiUDK2otcH/SLOmSqdY4neRFoLBd8aAccb0IwZM9odUs/2C+S1116L3/3ud3HzzTfHD37wg2RefPazn43vfOc7ccIJJ0TPnj1j7ty5cfHFF8fixYsj4u05tNVX0e/cf//9k0AE8vrrr8fcuXPjtNNOi69//esxatSoeO655+Jf/uVf4oUXXogf/ehHEfGnUjw6TvQ5OWdyaRgsp6NpC/RweZ55fee6C7eHSjegv/7rv47zzjuv8HLGjh0bTzzxRMyePTvOOOOMIq/gueeeS8zE5557rq5m21Z69+7dqdpRxpj3Nv/4j/8Y//iP/9gln3XRRRfFRRddVBeglLtBbdy4MT7/+c8X61t7FG19QlFT/6tf/Wp89atfrfuMq666Kq666qpWP3/u3LnxwQ9+sO4PtCuuuKLkaLZ/KnlAr732Wt1fQj169CieAIYOHRr77LNPzJkzp9i+adOmWLBgQRHvbowxxkRUvAF97GMfi0suuSRuuumm+N3vfhfXXXdd/NM//VOceuqpEfH2o+LZZ58dF198cVx//fWxbNmy+O///b/HwIED45RTTtkW+2+MeZ9y6623JnlAF198caN3yVSkUkvul19+OS644IK47rrr4ve//30MHDgw/uIv/iIuvPDCIoZ9awOln/zkJ7Fx48Y48sgj48c//nEcdNBB7fqOrS2599hjj+Jpa+TIkcV2eibUJDW5i14G9eVly5Yl65p7wrh6FihUr4AlTdhmgL5CFfR7+fRJ/0WhrEkvQ/eZASKq60Z0TNttDY5TW4EpEfUtlLkP1LE1L4X5UWx9oN9LKYa+IL3AbUVZzkp7YR4NfR6VeehZUf/XecsxpS+o55bfSZhTpOuaGxNR7y1pjhF9Wc4ZzS+iB8rfBl5bmgtXNm91HJknw9+RKnBu6k81rweWRspJiszHVD9YPfTWPpc+HH1o3dfNmzeXtuSu5AHtuuuucemll8all17a5mu6desWX//61+PrX/96lY82xhjzPsPtGIwxxjSEShLcu8FWCW7PPfcsHou1CixDjdkZUMMnKdOweyRDnLVcTa4yMOE+dEZCoQShYebc/1yVaj6+8705mN/A0h1VUJmktW6U7YVjnJNbqkA5QsuSRKSSHGVMylK5JOqtBSi3wkrgKiVXLWWjsNzU/Pnzk/Wc9EqpRKVkSi28tvS9OWk4on5uqtxVJldr1WemWfC6U9mtTBYkeh0yVD+X0tDR74iol7ty4dJl6Ht57eQ6+3JMq3QBaI0yCc5PQMYYYxqCb0DGGGMagm9AxhhjGkLTekAHHnhgoYmq1r527drk9SzZoho3PQdqrtTAWQJdYXirlt9hvyOGYStlXVqJ7jO7HrJzprK1ftZWqhRynDp1arLODG0tmUNfhOHfqiFXKVnEsH22Sc7BLo6c4m2FjraGhqRzHA477LBk/eabb07W1f/iPOW51HDjzujuZRq+emfawTKiPi1BPcayUjua/rBmzZpkW66Ef0R6bTFEu4p3mUuVaBT0dfS3gWHknfFIc/A850K0OX8456uWyLIHZIwxpinxDcgYY0xD8A3IGGNMQ2haD2jgwIGFD1MlJ0dhHD3Lb7D1dK4cCvNONHaeOjs9BtW12T6CuUj0SbSMDHNdciU1ylBtmp4Vxy2n+1Ij5nTK+WH0BrSkCcvEcIyZa6L5OzlvjHAMqYHnWryzHAqPVbV3+iI8z/palpdivpF+L+eTlq1qbR8fffTRYpmeKF9bJYdL95F+KT0fXlvq27LsDa+PXI6U5gtG5HOr2KqcOTb6WZzjvA713PG1vCZffPHFVva8OiwvRTqaF9cZNI+yVqvFhg0b7AEZY4xpTnwDMsYY0xB8AzLGGNMQmtYD0lpwWmOL+jHzdXLtC7Rkf0S9J1Sl7YBq3szpoN5fxZthHTbNq2FeQ5U8h7L2DDkmT56crKt3Ro0+V1tt0aJF2e9R74CaPT0h5nyoX/HMM88k28raMSvMQ9HLg95YWc0zHQv6E/SANCeNrQ5yuSSs/cbcNl4DOq4cF75Wc7HoAWndQcLux6xZeNNNN7X5XnqtPJd63ZXlc6mfx/HmeWbeGWvotReeK45pZ2urbYW/OfSW9HzRl+I+9evXr1guq5nH41PfXK/BlpaWeOGFF+wBGWOMaU58AzLGGNMQKjWkezfZe++9i8c9fexjaXGGT5bJIgqljiqoVEa5iI/zuXI7ZZKDSmWUEdg6YPny5W1+DyUsDZksCw2l/KLyhX5Oa5+VC5ulXKdhwZQCWFqI5XZUXqXkwFIkKr0yxJ/yhM6RXLuFiHpZUOctpUrKv7m5mGv78OCDD7b5nRERxx13XLKuZZQIQ7gffvjhNl+bg2HwvEZzZa5y0l5EtTJKObmL85TycC5NISdnszwT51euC3AVeH3wt0DnCOcef3N0vWrbCr229Dvb6+z4CcgYY0xD8A3IGGNMQ/ANyBhjTENo2jDso446qtDy77jjjg59FrXzKvoxoe6rOjZDYefOnZusb6shpg+ix1f2nbkSP2X+hGr21IwZHqq+VVlLYQ3fpY9AH4ch0eq7sT2GhplGRCxevLjNfehMe3V6G0OHDi2WO9PWvAoMeeW40Q9oBOPHj0/W1Z/MnZttCcsHVfVCtsJrhy0XNCWDaSFlpbV0zrM8Fq8HnfMdLWXWGpxf6rO1VrLLYdjGGGOaEt+AjDHGNATfgIwxxjSEpvWA2oJ6PmP0WcpeYQlztmPIlZzvqva59G3++Mc/JutVyvYQzQXI5QRFpDk3GzZsSLbRy2BuiZY4ofbMMa1SOkVbCxx44IHJNvpHudYBzAXjPuo5YDmjKl5AWUn/jsJjZ+6IXh9leTM5WCKH7RnK5lBbVPVTNMeLJaJ4Pahfyf3vzFiwBJZ+T66lCOH1ncuD23333bOv5W+DUqVlCmGenHpwzJ3ibwFzEcuwB2SMMaYp8Q3IGGNMQ2jaUjyDBg0qHv9UIqIUViW0mmVj+Lifk9VYZkUfYydMmJBsmzNnTrKuKmdZZVp2uNTjYzg0x6JKNW/Kbgofs7muj+zch9xYUHKgXKr7xNIolGU1xDkilUkowVGeyMkVlIsOOOCAYplVzsskH+3uyXPHuafSR1modE7WKUPlU849Sm56nqnUU87WY2XoOqUkyp5Lly4tllk1PBdCTBmK6DXLUk7/9V//laxXCZPnPFaJmvIjyxlpmDbnRE5yi8jPRaLde3lslNl0n7iNZcUee+yxNr9TQ/63VsMuw09AxhhjGoJvQMYYYxqCb0DGGGMaQtN6QD169Gi1I+qpp56avO6uu+5K1lV/Zvl8eg4Mr8y1KKDer7ov9eQcuTDxiPqwUz0e+kMMedYyMnvttVeyjXqshkbS3+I+8Nh1n1jGh/6LjjHDMemH6ZjSJ2D3zipdTon6SfQQWYqHpVWUspBU9ZrKWjkoHJcyb6AKGrpclkqgIfT0nXjeNfyY4egsOZODvkiu9UGZV6Zz6IEHHmj3PhDOW17D6ocxDJuh+k8++WSxzGMr+171fcpKRuk1W9aFVX0rhvzzdySH2zEYY4zZbvANyBhjTEPwDcgYY0xDaFoPaMuWLYWOqDr8vHnzktcxb4P5FUpZ3obGsdMDouaa08+1rUBEfWuBHLny/2XtxtUrKGs3rnotdWrqt7n8nbKW4grzsNgmXP0Y+ivMA6riKxAtOcPPKcuNqYJ+NssZ5fK/6AER9XE4p7nOc0nvRmE7AJ1P3MaSWAcffHCxzDFjrhiPT/0lXnf8Xh23KmVhON6DBw9O1ul7qsdIH4R+nuYq8dhGjx6drOucp4/J6yHnG/J4iP5mMl+KOXSaJ6TnMSLi7rvvTtbpS+k+dqSMmJ+AjDHGNATfgIwxxjSEppXgdthhh+IxU+Ukhh7nJDfCEEl+lso8DCWlDKUlTSgbVJHcuhJ9/NfK2BH15ThU7nriiSeSbexqmgsDzklupCy8O/c9ZZJbLoSecgvPZXv3gXIKy6zkpEqG/HOMVYZilecFCxYk6xq+S0lNS+JE1IfV5qgioTAUWbsWa8mYiHq5iPukY8HQ6tw+UdKdNGlSsv7b3/62zfcy/J5SpUp/ZWVvdM5wnt58881tvo/ljCjf8fdJ51BZaL6G0NM+4D5qugDnNKlSdbs9+AnIGGNMQ/ANyBhjTEPwDcgYY0xD2O46opJcV0SGf9IvyoU8d2XYr6LeUUR9WCN1X9X7+Vq2VFBdmOHR7HZZpvV2lFyXUJb84fnRMWboKH0cenS5bqTU2nOhyNTh+b1KmTegrQVybQUi0rEYNWpUso0ekMJx0lSCiPqSRmXh+W3BMHKGDOvcYygvw37vv//+ZF09IHplq1atStb1+JiWQF9N5wRDtumhTJ8+PVmfO3duscyfyCpeGeePluLhuaJ3pj5ORFpS584772z3PtArY1ml3O8g30vvTOdBa+PijqjGGGOaEt+AjDHGNATfgIwxxjSEps0DUjS/gvkHXFcNf9iwYck25hjktM8qng+1W+r9qpPSh6LnoK10I9I8G/pS9FQeffTRYpmtjXPtC+hLcZ/o1eQ+i16M5oCwnQRLnKxcubLV90XUl67PeT6EZWNUk37ooYeSbfQUdR8XLlyY/R62bshp30RLO+U8H8LcI3o8HfV8CL2ZXNkelpDRHKGIfDty5pmwlYB6EH/+53+ebLv66quT9Y9+9KPFMv2UG264IVnP5etwntIPoy+i0A/T65ItIvhbxjw5nfP04HJWfhXPh5S1ktffIJfiMcYYs93gG5AxxpiG4BuQMcaYhrBd5AFpLD3zWZgnMHz48GKZ+rh6JBH1uTKq2Wu8fkS950BNXMm1w6YfQR8n590wp4N5QDpu/Fwej+r0HG96Gcw10RYY1OxZg07PAfX8ZcuWJeuqTdOHor5MHV7JnRt+NvNB6Hfl2m5w/nCfNKeL45S77Dj+1OH1fJV5lZwHnEPtRds2R9S3QVHoo3FMH3vssTbfS++Pc0/PXVkrjSrtGnL7UfY5uTqEud+CRpGb4/QtOcfpd+nvouZhtbS0xLp165wHZIwxpjnxDcgYY0xD2C7CsLV8O6Uklt/QMiuUDVjihG0T9PGeYZt8vD/ppJOK5Tlz5rT5OYSP8xp6HFH/iK4hlPxchpnr4zFlNZb9UEmLr+UYn3LKKcm6Sg6LFi1KtlFe0fBphntTHsqVeqcEl5PZWHqHIdzafoLnNSe5cU5UaUVRhkrLlFYZcquyDmUoyh1sIaGpCGXlgXS+lclQuk8cl7IOr8rYsWOTdYb862fxXHVGciP621FWikdlN8qnDHlWSZrj0hl5jnNcpX7KriyN9MgjjxTLPFZK1AzDVjlPv6e9zo6fgIwxxjQE34CMMcY0BN+AjDHGNITtIgy7qzjwwAOTdZYaV31z6dKlyTa2+1VvRjXUiHptV3VUlqNhSfnOoN4HQ9CpEauvQN+AXg01ZA2N5Rg++OCDybr6AdTDcyHb1NJvvfXWZD0Xqlx27HruGDrNMjLaxp3zkqVSqsDW2ar/l7VbrhIizGOv0sJe4RzhGNM3zLGtWp0QPfay46a/p/OLng9bjuj1z7BrzqcqZXByHHroock654HOTU1Niaj3zvQapf/Iccl5Qvrb1tLSEs8//7zDsI0xxjQnvgEZY4xpCL4BGWOMaQjbRR6Q+gwsS0ItVHN7qN2yHcPEiROT9Vyb6lwpfnom9Hk0B4f6d5kHlNP7mQOlORP0J1gGR+E4UPdlvoiWRqLfRY1b11lmiO3Hly9f3uY+EObGKNS4WeZex42lXqjZa55DVc9HPQharZrbVga9y7Vr17b7vTnvI3ftRKRzlfOH51lbkmieVUS9B8HP2lYeUO7Y2SZhxYoVybrm1HGu8XqYNm1au1+bawfPVizMqTvyyCOLZc4fzluFOYAsF0TfR2HpnVx7+4742ZWfgJ5++un4zGc+E/369Ys+ffrE2LFjk2TEWq0WF154Yey7777Rp0+fmD59el29NmOMMabSDWjDhg0xbdq02GGHHeLmm2+ORx55JL7zne8kWdzf/va34/vf/35cfvnlsWDBgth5553j+OOP71CzImOMMe9dKoVhn3feeXHPPffEXXfd1er2Wq0WAwcOjP/9v/93/NVf/VVEvP14PWDAgPjpT38an/70p0u/oywMO/cISMpCUCdMmJCsL168uHT/WoNhhpQYtGwMq8syjJmP3So5UtqjNHPPPfe0c49TKOXx8Z6hyirBsfIx5ZY1a9YUy+wce/nll7e5Txynrix7UwWVDctKFnF+dVXIbVeikiNDhlnxvVevXsUyj4Xv1XnKEHNeo5SS9ZrOXc9dCStCM7xYw7ApFZNx48YVy6tXr85+rob1cxw4f2g36PXOc8XfupwkR1TmZ3g91/V6jkivj9ZKO3VpGPb1118fEydOjE984hPRv3//OOyww+KKK64otq9duzbWr18f06dPL/6tb9++MXny5LjvvvuqfJUxxpj3OJVuQGvWrInLLrsshg8fHrfcckt86Utfiq985Stx1VVXRcSf/nLmX+oDBgxo03TdvHlzbNq0KfnPGGPMe59KUXAtLS0xceLE+OY3vxkREYcddlgsX748Lr/88jjjjDM6tAOzZ8+Oiy66qEPvNcYYs/1S6Qa077771rU0GDlyZPzyl7+MiD9pv88991xSzvy5556rC/ncyvnnnx+zZs0q1jdt2lRX6kJhWRiW09f3qt4aUa8vM/RSwyBzXRsjUt33oYceSrbRM9F9ZEgwn/hYKl3DILn/hxxySLKuujytPYYba+gl9/dDH/pQss4AEtV06QXkQqAZbnvEEUe0+Vr1HyLqzyXDQ6uESGupeur7DIXVEjnUuDlPeXwanFMlRJUeCt+roe+58PrW0JYXZd6G+iD0gHKtA+gTlKEeKceUXkdHoSpD/4VeTRX/Tq//HXfcMdnGc8kQdYWeLj9L0xQIW8/o8TDN4s4770zW9VqjD1V2LstaepRRSYKbNm1a3Y/2ypUrCxN76NChsc8++yT9cTZt2hQLFiyIKVOmtPqZvXv3jt122y35zxhjzHufSk9A55xzTkydOjW++c1vxic/+clYuHBh/OQnP4mf/OQnEfF2EtbZZ58dF198cQwfPjyGDh0aF1xwQQwcOLCuqZkxxpj3N5VuQIcffnhcd911cf7558fXv/71GDp0aFx66aVx+umnF68599xz49VXX40vfOELsXHjxjjyyCPjt7/9bd3jpDHGmPc3TduOYffddy/KWmj+CFv0UiNWzZV+C70MloZRXbgs70RzF7hPuXJBLHfCMj7U9D/ykY8Uy9pGIKLe/3r00UeLZfU5Iur1fvWaci0IIurbJKuPwD8smDegpUmoJ7Odt/pq1ORZTSNXYp6e3OTJk5N1nfJz585NttHTqtImmfNJfTZ6VlW+p+xcVkHboNNPzVGWU6etv+mfco7k2mOobxaR9xirtIBgiRx6JrwO1dvgezk3dZ/pNbHNvPqE/D3iOLHtvB475wTbk+jv4kEHHZRs42+b2h5lv3vcx7ICA27HYIwxpinxDcgYY0xDaNpq2CqfqfzCci6U4HJVhvnYzTIy+sjLx2x+r8puZY+h+lqWAOE+jBw5MlnXEFxWsaUMpY/lPFaGYmrY8uGHH55s4z5SqlGpqUzGURmKEhWlPpVPGRpOiYEylM4DhvIyNJwVvLsKhoJTalIouek8oLzVlRKcni9KuJRq9BpgJfncPGbIP6VWypEKQ8455/Wz2VGX51UjdimFEYbQU/7KoRLiSSedlGxjeSz9XWNZMZ5XzluVt3k8uTmxcuXKNrdF1Ev7OXK/dXp9t7S0tKt8lp+AjDHGNATfgIwxxjQE34CMMcY0hKb1gIYPH174AKol5nT1MvjeBQsWJOsaislSMCzFr6HWkyZNSrbNnz8/WVevgz4IobauoabUaqlb6z6yBAhbLug+cRv3gceuWjTHlCVONAST4awaEhyR+j4s8cHvoSekobDPPvtssi3XiZXb+L167NwHhu7m5iZD3Zn9oOts0ZFr6Fi2/zl4ruifql9BT47+xK233los0wPqDExLUL/y3nvv7fDn0p/IddjledZ0gYi0k/LWwsxbYVi5fg9DtunJlXX2zX2PXtNM/chBX5a/Obzu1K/UOZIbT8VPQMYYYxqCb0DGGGMagm9AxhhjGkLTekCqe+e8AZZ5yDW0K8vXUb155513TrZRR1X9nHka9AJU26WGynwjatPq5dA34LFqyZxjjjkm2cbvVQ8r10Ygor49turCzL1gGXtdZxkPeku6H9zGvAeW/FHfhC0tqK1rSfwyz0Q/S3OaIspbaeSo0rKAfp6OBb0ZXg8sZaPHzlyXnMfA19K/0PwdeksshUTvVSnLeeqM76Ow9QG/Rz1g5iYtXrw4Wc/5nLwe9BrOjUNV+L2c8+2lLCeIZcaUXKuJtvATkDHGmIbgG5AxxpiG4BuQMcaYhtC0HtDOO+9c6JqqS7LeFmu0aUtl+jgPPvhgsk6dtH///m1+Lj0gzUeYN29eso35LcOGDSuWqTVT86Yfo9outXT6IKovMxeD3oD6VLkW4hH13oD6Cqwfxu/VmmfcxhwPPV/8TnoQ3Gd9PfMneDw6bsy5Ye0x/R56YzyXzNfR46Hnw3067bTTimW2k+Cx0xNS6DkcfPDBybrmh9x///3Jtpz+z/lD71JhPcMyn019xLI6d1pzrqwNu44xfSm+l3l/WkuxrJ36scceWyxzTDnH1Ycua/PAfcq1CR88eHCyrnOR112uLh79a/qabLXRWfwEZIwxpiH4BmSMMaYhNG1HVEXLrlPm4GNrlTLqlEFUWmLoMcu7jB8/vlhmSwjukz7us1w7wxr5Xi21QrmOjBo1qlhm+G2ulA0fycumhJZsYWg1pQ6VV8q6i6pUwH3g+LOk/+9+97timWOYC31nSC3HTY+V43/HHXck65TzVCpbu3Ztsu1Tn/pUsq7SDEs5Uf5SiY6lU3QORNS3QtB95PnIybYs+UNZSqUmXlccFy1dE5FK1KtXr44c+tvAY89JcpTUOU91/nQGtohYvnx5u99b1nVWUashov76OOSQQ4plytWct1WYMWNGsq62htoHtVot3njjjdKOqE3rARmzvfG5N96IL8qNpId4KvQCdv3P/0zW9Yb7Bv4g2PqD/kqPHnHukCGROkTGbL/4BmRMF7FHrRYH6l/WGdM4MgnTOf559eqYFhG/69C7jWkufAMyZhvwSkRsESlqCyQfRjgplBC7desWvVtaYsdaLQa+9VbMjYg/C9+EzPbPdnED0vBQlgDRdgWE2iNDhh944IFkPddm+8Mf/nCyruXnGfY7YcKEZF3L4lCPpddEv0I9FP5o8XhUa6c3Q71cv4fjxONh+KfuM/eX2rr6S9wnhjwrDEfnec+FDDPslyHduh9lbYP1+CijnTJjRgyQ4xu9dm3EO+VI/nHChLhXzge9JXoqGrL9q1/9Ktk2ZcqU2PW11+Ls66+PfTdsiCERsaBPn/jq1Klx5dy5yWvpBXzwgx9M1ul5KfQCdZ1lh+iHqSfEdvBMSzjyyCOT9Vy7Cc5b/Sz6gGwjotda7obfGnoM+rsQUV8mSq9RXs/0wxhir9DzYZi2yrT8I4Weqc7rMu9VvSfuA0P+6Xfrb5B6cC0tLXW+Z2tsFzcgY5qR/V59NX54553vyne9vNNOcelJJxU3of6vvx6X3Htv/LpWixfb2XvFmGbDYdjGbCdsvQk9/c7TUv/XX4+2U1KNaX78BGRMF/D0zjvHI1Jh4QXIjV3FyzvtFI/06xf7QX4yZntku7gBdTRGn7o1WwdQM1YtlDo2/QqNs6eXRK9AtWpq0SzvwrbI6hWwdTbfq7o8j41lezgWCn0qegPqi1CHpzat6/xOtjpWzZs6+/PPP5+sU1tXT4jHSu65555imRo387J0HPfee++YvmZN7PTO+X1VclZub2mJn6BsyZviY7GVw2GHHZas6/FoaZeI+vInb8r86tOnT+ws54saPb2yXG4JUb+F458rr8P5w7wf5uTwelG4/+qV8Vzxc3Tecgy5D5wHOjeZ20OfTX1Pllx6/PHHo6PwetHj4xjzt07LitF/pCeamxP8PeJ6Z9kubkDGNAunPvZY9M/0RDHGtB97QMYYYxpC0z4BTZgwoXjM1BDcG264ocOfyUd0hofqYy3DKVkeRcNbp06dmmzLdbBkWRKWP9FqvxGpVKCP1RH1j+gqPTHkmWGauY6JlLBYVVzDshmGTblO95Gh1dwHHQseG9cpaemYM4Sex65lcFjlnLKgfs/8+fMLueKlbt3iEjk3z/TqVVdKZeHChcUyS5hQKtPjYzmaOxFp90k5vmHDhkV3OQccU4bCctxyMKQ+h8rDZZ0xKRdVkQWXLVtWLLMCNM+7XkuU65hqwH1QqYmSW04yLAt5rgK/R681lh1iyLZu5+9cDnaKHT58eLK+//77J+uXXXZZuz+7NZr2BmRMM/NGt25xBzy4vm281hjTOpbgjOkAe7e0xIyS3jXGmDy+ARlTgVvekQm7R8TfrV3rm5AxnaBpJbjVq1cXumYuDJvhuqpbU+elpkrdNxe6yM6S6rEsXbo02abdUiPSkFWWtaEvwlBr9aK4jaHWuk/0TOiD0K/IfS59BdXwuf9sh6G6NUO0+bnqwfFYOW4MB9WSTPwehqyuWLGiWGYoOzV8fe2zzz4bX6zVYnO3bjGzVituQn/csiV+s8ceMRdlcaZMmVIss2wSw4JVs+e8HDtiRPSRed1Nyum8+OKL8XsJVWYYcO48E855jlsOHUf6EQzVZ0i3zlWWhWL6g17fPM8TJ05M1nWMc6WbyuA48HdEfRH6W7we9HrhONGrYQmgQYMGFcv0gOgla0pGrlxZRDpHGDbO72GZsbb2r6WlpbRjbUQT34CMaUq6dYuz3lncehO6+J0LvP2dX6ox9aWX4h9WrtxGn25M47AEZ0xVunWLs7p1i2vf+Stz603o/6sQNWaM8ROQMR2jW7f41jthvp/+wx+ie0T8tFaLaGmJazIh7p3lyR13jIcl3PhVFyI12zFNewPabbfdCo8g5wHRK1BfgaXcCb0OLSl/1113Jdu0bXBEmhvA0vTMhVFtnSVMWPJHddSItAwI8x6I5jpQe6aOrSWBcjlBraE6PfVw5o6o90QdnuPf1v5F1OeWcIz1HDBvg8eee22uDQTn4VNPPx1n1mrxSu/e8ZebNxc3oZYtW+LqSD0Kauv0J/XccX93lfG+Y8iQ+JrcgMpybnI5NvSH6Bvmyq7wutO8pj8gMIPzi7kkOq65Fh2Ec4RjWsX3Yd6c7nNZi3r1DavMcebXleXJqR9J+FvHPDOF46/7zzw4tv5mOwn1nvS3bMuWLfaAjNnmdOsW573zQ7L1JvSvEdEtIjpeBcyY9wf2gIzpLO/chH74zmr3iLgqIj6CaC9jTErTPgGVhQ5uRUtzRKSPhMcff3yyjWGNlI+0BA0r+E6ePDlZ13BLhoLffvvtybpWzqbswdBLhqFyPUeueydlBIaAdhTKIFzX88FtDM/VMaXcRTmI0qVKUZQYKNdpqDIlB5ZgUmmDY8jw+y+/8/8z4+2b0AWrV8fjq1fHv7XyuRwL3UfKmOvltX948cWITDdbdlqlrJML91YJOiKd1xxvdpJl+LHCY+W5zZXxyVWwZ8h5FTj/uY96PAzV53xS2Y3zieHRKhNWDQ3X3ydWtGapoRwcU00/YboDx5jnWcdRS0+1e18qv8MY0ya8CV0VEX8eET3f8TZ+uMsu8VAr7aG//Oqrccg7fzTwRveBCvk4xmxP+AZkTBfDm9DJERHvBA78AsVbt3L4m2/Gcfjr05j3OvaAjNkGfDmi8ISMMa3TtE9AO+20U6Ev0itQGAKt3sd1112XbKNezhBJ1TdZcoKlU1TrZYuFXKkRti9g+CehNt1eqEVvK3g81NY1pJjl8nNdWamd0ztjOLuOU1npF91Ob2PRokXJum6nZ1LWRuCqCRPiN+/Mx6XvhK+++OKL8UbUd+RUPjN2bGyS86fh3K/26BFPPPhgm/vPY+V6DvpUGhLN+V8FnneSCyVnCxK9htlagtezzkV6F5xP/B49t/wehslrygPbqbTXy24Pud+KKq0zmBKgodZs31EWgl62vYymvQEZs93TrVs8/86N8ZnWttdq8eg7F7z+mfGHXr1io96AtmFiqzGNxDcgYxpIx55vjXlv4BuQMQ1k+TtPSFtEQtni8jrmfUK3WmdFvC5m06ZN0bdv3zjqqKMKj0DzLajdUl/OtRxmSXzqprp+2GGHJduY86HtGcrKV+h2vpaeyejRo5P1o48+ulim17StYJn7nFfD6UNfRP0W5qQQ9RlYWqSszH2uhAtLgugcov/AfVRfgbkvzJmg36K5ZBxD7n+VPI5mh/4d82bY8l1zWuhzsL29lk5ifk6udJDm4kXU5/lxzudg7l4Vn03L67DUTpV96Azcfy0zxhI+vHYmTZqUrJfl/rz00kvZXEaLy8YYYxqCb0DGGGMagm9AxhhjGkLTBiGsXLmy0EhVP2cJc6J5ESy1z5wJxsOrH7N69epk2/Dhw5N11Z+p7zOWXlszl1luzKtRb2NbeUD0I6rYgvSw6NVoaX5+LmvzadsK1rqizk5vQL2DXFvwiNSjYN4Ga+ipL8VS9Gw7wPmln818ENZsUz+JOU5VKKstqPOaNQxzHt348eOTdfqnS5YsKZbp+ZBcPgvbgtPTVQ913LhxybacB0QPl/tQ5RrIeT4HHnhg9nuWL/9T31z6nKwNR7+4LO+svXD/dY5wHAg9U72W9DpraWnJ+vFb8ROQMcaYhuAbkDHGmIbQtBKcloZXmYfhqixpwlBZhXKEdhuNSMNDKXdRLlLJhHJKrvQIQ1T5OMzH8Fz59rKOr+2Fj/ZlZT20FAlDSVk6SMeRodLabTQilUxZVoUll9g6QMeRMhrR42MoL2UolRFGjhyZbKPUpDJURMSRRx7Z5j4wFFmPj+eVkq5KHZwTnE+UjlWCKwuL1+suJ7mVQdk8J5VRfuR7TzrppGL5+uuvz36vjhN/N3iNUqbVOUQ5juHSei4p61MSVShXcx8602JByyqVlQPSOcFyTCxRxOtOKZNeW8NPQMYYYxqCb0DGGGMagm9AxhhjGkLTluJpqx0DS6FT41YNnDo7QzyJluegN0BvQ7XdadOmJdteeOGFZF11VHpW1KJZikf9AG1d3No+aYkZeig5yrRmhlqr/0WNnr6CvpbHyjB59RlY5obtANjWWcNdGfpK7081fPXYIiLuu+++ZF33OdeKOaLeu5k/f36xPGPGjDb3ISIN8abOzhbjeuz00XgueXw56LPpHOfPxIoVK9r9udxHXof6WVXK0fA857zXXAuFiHovM+eDcv8HDx5cLNPXoael+1i2T7n2GDxW7r/OxVy79G2NS/EYY4xpSnwDMsYY0xB8AzLGGNMQmjYPqC1fgjopyz1oORSWGrnnnnvafG1EqulTY6Xuq2VAuK8sI6MeEHVffi59BfV1mK+T08uZF8D8I9WFy/IcuK6aLj0G+iQ6FsxzoOegHhz9LR4PP0vPF70/eirqL9G3oa+g84vnRsuqRNTng6k/9sADDyTb+Fm5tvP0u3TccnkZVWHpFz32Kl4S4T7SJ9yhjfbjraF+GPeJ16x6prm23xHVWlrTp9Xrkt/D61s9ah4r5zznhG6nrzJq1Khkfe7cua3teqfhdaf7kcvvaotOPQF961vfim7dusXZZ59d/Nsbb7wRM2fOjH79+sUuu+wSp512WpdeJMYYY94bdPgGdP/998f//b//t64g4DnnnBM33HBD/OIXv4h58+bFM888Ex//+Mc7vaPGGGPeW3RIgnvllVfi9NNPjyuuuCIuvvji4t9feuml+Jd/+Ze45ppr4kMf+lBERFx55ZUxcuTImD9/fhxxxBHt/o7ddtutCP/VqsKUOW6//fZkXSWVXPhfRH3JFu2ayJI+7MT6m9/8pliePn16so0yG8NQFZa6YPi0So7sVshQZQ2bpeTG71HZrawDKkM8Verga1nCSKUNPr4zPFSPnRIVQ8EpV+j5oqRIaUClMcq0gwYNStb12Cn5UIJjuLRWy6bEk5PcymAVbqWsnIuOMfeJ81a7YX70ox9Ntj344INtvrYspJlSmV6nPM/cf5VXeZ4pfev38jxTyud8y8HSSBqGzerqlE9zcPw5R3Q7uy5TctP94LHlSjBRMuRvJN+r46oV6ltaWurKl7VGh56AZs6cGSeeeGLdD+/ixYvjrbfeSv59xIgRMXjw4Lr8CmOMMe9vKj8BXXvttfHAAw/E/fffX7dt/fr10atXrzoTeMCAAW0aVJs3b05MvEYmTRljjHn3qPQEtG7dujjrrLPi6quvzjaVqsLs2bOjb9++xX+UQIwxxrw3qVSK51e/+lWceuqpiRa5ZcuW6NatW3Tv3j1uueWWmD59emzYsCF5ChoyZEicffbZcc4559R9ZmtPQLwJTZw4sVimDk9PSEtfsFwIPaiFCxcm61omh74INWQtt8Nwbvo4WkKe3gyfFqkD5zwIvla/lyHO1KZVa8+FWUfUa9Hq+/BY+YeJhp1S32d7BvU2WMKEr2Vkpe5TWUl/fcpm+CrHWPefev7KlSuTdZ2nERH/8R//EdsCLRPFNATOJ/oiSplXo+NGbyAX1lzW6ZPnQ/2YXDkdwnPH79E5U8XjKSNXQoftFxjarsfHOb2toF9K71ihj0kfp0qLiIjyUjyVJLhjjz02li1blvzb5z73uRgxYkT8zd/8TQwaNCh22GGHmDNnTpx22mkR8fZN4Mknn0zaUiu9e/euM/iNMca896l0A9p1113rmrjtvPPO0a9fv+LfP//5z8esWbNizz33jN122y2+/OUvx5QpUypFwBljjHnv0+WVEL773e9G9+7d47TTTovNmzfH8ccfHz/+8Y+7+muMMcZs53T6BsT48x133DF+9KMfxY9+9KPOfnTBmjVrimXGpTNfJxccoeXxI+pbLqj3xDa2ug8REZ/5zGeKZUb40TPRtgPUj3k81Gf1s6gv87PUv+A4MHdEdXrmjjBHgr6Vfg/3get6PNSPc22S6UuV+RU6xvS7eD40Z4V+BNFjpa9Gcvk5nYHnUr20Qw89NNmmxxZR711qHofmbUTUX0s6bmV5S+oR0Ydim+oqOVC8RtXzosdLO1vnAceQLVOIeoqci7w+FLZEZ1v2nAfE6zDnb9NPzY0pLY6cB8TrjMfO+ab+kvqyLS0tdXmLrX5f6SuMMcaYbYBvQMYYYxpC01bDVvRRlZVotYJyRFpypqwrKGU1reRcFi6tlR34OM9HZ5VxKLmxlA3RME5KDJSENJydj9mUDTS8kmGzlKVYhVvHiVJlrussQzoZSq0SBKWARx99tM19iEirk+fCuyPSckHsysr3alfKJUuWJNsYlk3pb/LkycXyggULkm0spaISEUv8MDRZzyXHn/vPcVQouVXhwAMPTNZ1HMvkrRycP7n9L5OKcxJQWah4lc6sCveXsr9+L2VBStI8Ph0bzlv+1umc4fwZP358sq6V2rkP/B3h72Bnq277CcgYY0xD8A3IGGNMQ/ANyBhjTEOoVIrn3WDTpk3Rt2/f2HnnnYtQQ9VnGdZInVdLhFBnf+yxx5J11fcj0jBnasAskqr+y+GHH55sYzkX9TrYmoF6OfdJNVe+l2HZqgNTDyfqLdHjoY7NUh56fGVemYaLclzoh6lWnfPnIuq9Mz32hx56KHLouWQ4K0OTVWtnGX56SyxPo15TWdkVLVfCuVbmVygcF/qRVdoD5Bg2bFiyrueH5X/oK+RaRrC/GN+b83U4n9ST4/jT52xP64CtsLyOXi/0TDj++rvC88xj5U+zXlu5EkukrEXHtqSsFI+fgIwxxjQE34CMMcY0BN+AjDHGNISmzQPS0hKq7zN3IVcWviwfgdq66sTUjFnqQveJ2jPzOLRUDz0SasQsPaL6Mttqc/9Vh6dnksslYbkW6svMV9B95LjQT1KtnZ+7evXqNr+XfhfzXXKtA6idszyKnjseG1sOa54ZS5owF4Nad5Vy+7lGjPR89Hs4Jzi/mGulnh3nSBXYBkJbmZTlF+VaLpT5dzk4/jm/aNKkSck6c250frG8FH1P9Y849/je9pSn2Qr9PPWXmAOYKwPFa2VboZ5bS0tL3fXRGn4CMsYY0xB8AzLGGNMQfAMyxhjTEJrWA1I0f6SKbk0PiDkebKl86qmnFstLly5NttEb0BLzjLNnvTr1iKjHslbUjBkzknX1eXjsudwYejOssaW+AvNzWOb+qKOOStbV0+I+0INQXZ46++DBg5N13Wf6HtTSOY567Mzt4WtVh6eWzjnB3DGFpemZG3PnnXcWy2wzr/W3CL0M+kN6vjTXKCJi0aJFbX5uROp5lbW4UHjtcM6oF8u8Ma4PGTIkWVdPiB4K88FycJx0HvBz6U/kPBT6kZzHev0zD4h5V+rZleUeMRdRcwT5m0NfSv1H5jx98IMfTNbVd2NrGULPUc+tjmF700v9BGSMMaYh+AZkjDGmITRtKZ5ddtmleITm477CEGJ9DKe8csQRRyTrixcvTtZz5XUefPDBZF1lKj52U8qYNm1asXzIIYck23hs3GeVlhjGzDBslWP4Wsog+njP8j+U1RgqrpIi5S7KOioNMPxWOzxGpI//uTDriPpQXw19p7TH9+p+aPhwRMQtt9ySrOt8otzFSycnkTL8ludDpcsqZWHKoGSi56dMblEpmV1NOaa5FAC2CuC45UoL5ch1S41IryXKWSyn05kWEjk4v/bdd99imXIpw65Zakvl05w0HJGeO153TH/QOVEWss1xGzVqVLGsvzlbO6K6FI8xxpimxDcgY4wxDcE3IGOMMQ2haT2gwYMHF9qkeiydaSPMsiuqX0akZco1JDsi4uabb07WtewEtdzDDjssWVd9mb4N9VFqrPo9DP9kuLS+d8qUKck2hp2qLsywcZac4bip/k/Ph+HfLFOk8L3qt7AFN1/L49GwZtXZI+pDuNXLoY+QKxOjbdgj6svcVLmU6Ccp3P9HHnmkze+hv8JQcIYB58r4M9Ra5x7HkGHAudbM9AlzpW1yfm8ZubYVbEPNcllMU6jyO7P//vsXyyxFxTbtuk/01ei/5PyYMv+rq6DXR1+qrL2HPSBjjDFNiW9AxhhjGoJvQMYYYxpC05bieeWVVwrdX/VylpxgmQzVbpmfQ32fPoPmL7BUysiRI5N11eGHDh3a5udEpG3Er7rqqmTbn/3ZnyXrTz/9dLKuOjA1eur56pPQH2KpEfUZuL/MmaAnpK+nN8O21epL8XuoYys8zxwXfq+usxU79f5ly5YVyzw27r/mcNFjYHmmXJmoiRMnJuv0bu69995iWVtAtIbOmXvuuSfZRs2ecyRXnoZ5QernsXwO91/zjdiGgvPphBNOSNaZe5VD5xP9ReYTaa5VrvRRRL13o9Cv4/HlfBD6X/o9nHvMAaQfpj4P30v0+m5PW4S24O9elXYS7cFPQMYYYxqCb0DGGGMaQtOGYe+5556FrKLh0dxdhipqqDIf/cvQ6sbsokkJSEMvy0IRP/e5zxXLDOVl2CylJX38Z6glJSCVI/mInqtSreMbUS/XDR8+PFlXOWP69OnJNsqeej74+M7wbn0vw7lZnobnVkNauQ8sR6PSDWVZylA6FiyVwmrYuUrUPNaDDz44WVepjF1+KXflQpXLOmVeeOGFxTKlL0qVPAcKJdKc/Mi5yNBk3UeWKOJ1qNIS5wTTFFSyLis7ROlSyY1DRFrpnJ/DclO5skMcF/7mKJT2+Lt44IEHFstlJbBy8PeIx5e7liIchm2MMaZJ8Q3IGGNMQ/ANyBhjTENoWg+oo6hGzA6J1N2phVYpg6++AsNKf/rTnybr2s4g5/FE1Gusqp/SC6CvoDo3S/FwLPRz6RMwJJXtGrSUCnV3fo/C8ia5EFW2gKBfwVBrbZdBX5ClbbScCEssER1jeklPPvlk9r3qu1H7p4dVJWxWQ7o5Z+kLMoxW94PhxCxN9V//9V/Z/VDUE2W6ALv+MnVC5wzDvdl2QP0ktjZZuHBhu/e3DL3WysoD6TXNa6WK31KGzid6brxm9XzwuqNHqnAOsAwU0XOt+1Sr1WLTpk32gIwxxjQnvgEZY4xpCL4BGWOMaQhN6wHtvPPOhUegJVvKSkFofgJzbAg9FNV62baWMfqaq0GPgXk12p6B+ndZ2wHVb2fMmJFsYwl59YC4v7mSMywxQ3+CeRD62Yz9p16upVPoNdG/0JJGbEuh5XMi6ksN6THQ62OelurlLFWTy2chZeddL62cNxaR+jq5fKLOoj4C/Ub+FHS09QnnE32RlStXJusdLRvDnCH6gnot8XM593js48aNK5Yfeuihdu9TGXr9s51HGVV8KR0blrzKfW9ZfhG9WX29zqctW7bEqlWr7AEZY4xpTnwDMsYY0xCaVoJTVPLJlaeIyD+m8pGdVaxV1uFjI6WmXCjjX/7lXybrKrMdddRRyTaW12ElZO2mSomBkk+uajilJpXv9Dsi6sONKVeoJMd9oHynkkNOMoyIGDx4cLFMuY4lWVgJ+e67725z/xmue9ttt0VbsNOnzidKF5xPDEnXcPXcfIlIQ8N5rvg9udIwGn4bUS+h6HwrkwWroOeZoeAcFy1dExGxZs2aDn0nu3NyjuRgV1Ze77kOo5Tu9feJ0j059thji+VVq1Yl2xjWz7QLHUdedzn4e8rz0xn0s3U+bZ13luCMMcY0Jb4BGWOMaQi+ARljjGkI24UH1FXQg6Cnol4HQ0WJhhxSi2YpEh1itjagv0IdWP0MhhNPmDAhWW+tHPpWqI+rZk/PpKxzqXoUDItnuLT6LwyP5rGrJ8RpmfOWItJj57kr8w0Vzgn9XI4/w40POOCAZD0XvsvWDnp8uW0RqV9BH+eggw5K1ullqHdGb6YstLe98Dqjx8hxU09Lw58jIpYuXdru782FENPzYZg/rx29JnjtqFcZkZZ+4ufS39I5wpJRLNuTayfDMcy1eegMPFZ6mXqudX8dhm2MMaap8Q3IGGNMQ/ANyBhjTEPoWf6SxnDAAQcUpTS0bDzLh3Bd/ZclS5Yk26hFswUxtWuFWqjqmvwc6rH6vWy3TM1YW+lGpDkH3H96NaqBs6QPvSfVa+nj0F9hboO2fWZJf5bTpyehsBWF5mXde++9yTZtAx5Rf941F4P7RNTnYbl/ruvx0Efgec7lZjB3hH6Fng/mItHnyeXv0MehB6TjxvHnnFEvtqz1hMJ5SnJ+RVkejeYQlfktCv0urhMdtzJPuqw8mKI+D/O7cp5PROr7sCQO90F9NV4rzC/SvCD+prCMEnP3dJ8ff/zxYrm9oQV+AjLGGNMQfAMyxhjTEHwDMsYY0xCaNg+oe/fuhU6umj01SNbNKtN2uwr1gKjlcp/GjBlTLFOjp+bNWl7qS9EHYd6Jvpa5C8zl0Zwc+lLU6Bn7r9/LMvfUsdU7o6/DlhHqPWm7iIj6vAf6e7k24aRKWXv1HKilM4+D80C9QXpj9IT0fDGvjF6H1gvk5cv5Qz8v12Ih52XQc2CukuZIsfYefSh6WOo98bW5moys5Uj/rjPkatvl8u0Ix0LHmF4l8wfpKeo54G8Mfxc5N9sLfXC2+ub1UnbdOQ/IGGNMU+IbkDHGmIbQtGHYw4YNKx7v9DGQZfpzaPfHiPrHVg0n5mdTYhg9enSyrqX5GTbLx1+VsNiZlKVfKCOoVEMpRsMeuX7ooYdm9+mmm24qljlODDvldoVyHaU/3Wc+znOfVM7jY3tZOwOVDso6TebC1Sl/aVhqWUkZjpueO0qglBB1HlDKq9KRk+WCcl1/y8ZUYfsLlo1Rpk6dmqxzjn/nO99J1nMh3rwO9XiqtHGg3FtWnqlqt9KtcE5wzuRSBLiPLCn1yCOPtPlehlbnpOXca7kPua7FEenx6LHWarW66701/ARkjDGmIfgGZIwxpiH4BmSMMaYhNG0YtqK6alnZcWqYSi70det3txf1KNjOgPqyeg48NrYomDt3brKuYbXUl1kCSPeJXlOuXfExxxyTrLOlgoaRR6RhqPRQ2C5a9WaeO/ogOViiiOdKj5d+BcuJaKg4w+AZLq2hsAzDZvt0jpOGqDJcndp6FdS7LPNEOW7qt9AD4vWh53nGjBnJNpZs0eOhV8FQcPorWhLosMMOS7Zxfmk7dfp1I0eOTNY1fJqh1PTKqrDvvvsm6zr3OvO57xY5b5wtXhYvXpys04fmtUYchm2MMaYp8Q3IGGNMQ/ANyBhjTENo2jwgpUq72T59+hTLzJNhS4Wcfsm8DeYjqKdS5h3pe+kFcJ2lVHKaMrXcD3/4w8Xyo48+mmzLeUD0AugfrVq1KlnX4+E+MPZfcxn03ETUlx7RciK5ltvch4h8GRmWNNF1lg5iDo7mZdFP4bHTk1PYfoHo3My1UIhIx4n+FseYHsqkSZOKZeakDBs2rM19Zi4PvSf1+ugBleXUaA4X/dNcSRnOac69Kp4u56JeA/RtmZfFc9Be2N6jShkxnmd6cjmYD6n7UeZhcS6qJ5T7jWmLSk9As2fPjsMPPzx23XXX6N+/f5xyyimxYsWK5DVvvPFGzJw5M/r16xe77LJLnHbaaR1O6jLGGPPepdINaN68eTFz5syYP39+3HbbbfHWW2/Fcccdl0Q+nXPOOXHDDTfEL37xi5g3b14888wz8fGPf7zLd9wYY8z2TafCsJ9//vno379/zJs3Lz74wQ/GSy+9FHvvvXdcc8018ed//ucR8baEMXLkyLjvvvviiCOOKP3MrWHYe+yxRyEBqGxw6623Jq+vUvW1irxFKAnpY2xZyQkNbaQ0xtBqltDIPVozhFtDbLmN33PfffcVy5QfykJ7VaqhbMCqvNoRleHd06ZNS9a1tArLrJSVLFKpiaVgGEqqY0qplfKXQtmPcy83FymF5V7LisSUW6p04CQaws1w4pNPPjlZ15QGyuCszK7yy5w5c5JtHGN+loZha/XxiPpK7fo9nZGhSFnZrhx6bssqpGuoexWJsDOwVBjlYP1t4BjmKpdHpJXxaXlEbOMw7K2x9Vsn0OLFi+Ott96K6dOnF68ZMWJEDB48OPnBM8YYYzochNDS0hJnn312TJs2rUjAW79+ffTq1avur4kBAwbU/WW8lc2bNydPEO/WXwXGGGMaS4efgGbOnBnLly+Pa6+9tlM7MHv27Ojbt2/xHx/BjTHGvDfp0BPQmWeeGTfeeGPceeedSZe+ffbZJ958883YuHFj8hT03HPP1fkMWzn//PNj1qxZxfqmTZti0KBB0a9fvyI88Le//W2b+8LQavVBqHF3pkwGdWv1dRjKu3DhwmRdw1JZCp0h0NRYc/oyy7mofsuQ4JxfwXOjvk1ExN13352sswRNjqVLlxbLLFVDvyvni+TaVESkWjtL/OS8gZznQxiCyjBafo+u069gySKFYeMMj9aQdHoM9FvY6kCvCbYj4WfpXOW85bWlYzNu3LhkG31PlvTX42F5nUWLFkV7yZUSIhwn7pN6XCw/w3U9X0xhoJ+n13OZ2lOlhQTPj56PttSnrWiYedk+8RrNtZdoD5WegGq1Wpx55plx3XXXxe2331530iZMmBA77LBDYkKuWLEinnzyyZgyZUqrn9m7d+/Ybbfdkv+MMca896n0BDRz5sy45ppr4te//nXsuuuuxZ21b9++0adPn+jbt298/vOfj1mzZsWee+4Zu+22W3z5y1+OKVOmtCsCzhhjzPuHSjegyy67LCIi/uzP/iz59yuvvDI++9nPRkTEd7/73ejevXucdtppsXnz5jj++OPjxz/+cZfsrDHGmPcOTduOYffddy/i1VWnHz9+fPJ65uBo2Rjq7izJXgXq4+p9kbVr1ybrqk2X5Sow10T9JeruLCWk+jM1bvpf+j0sJUKtmeM2ceLEYpn5IIcffniyvmDBgni3qVoyvquokuvTGXQucr5Q6qZ3Nnbs2GKZ84mekHoDHFN6olpSinlA9KE4Lup1UNan/6itT+iB0nvSMlBHHnlkso2+VFmL7q5Cr8sy/5H+tubz8LzSa1Iv59RTT022MR9PfejWcnlyaG6lVryp1WrxxhtvuB2DMcaY5sQ3IGOMMQ2haatht1UKg5V2WXJG5SENAa4Ky1VQxtF8pXvuuSfZxtIXKmUwVJSvZRikHg/lxhdeeCFZV8kh15U1IpUyGQ5dJlVqaCxlTko1CsvpULLScNbOhMwzrJzHp3OL5YDmz5+frKtcxFBX0pkkapUpOEcoNalkwjD+sorQWsrmqaeeSrZxLFRmY3kmXg8q3TDslx13KdHpXC2TwnKdZFml+vzzzy+WZ8+enf3czqCyIKV6jrF23KX7wd88Spesxq7wt0Bh11/+hqrUV1VG1utUr4/2Ojt+AjLGGNMQfAMyxhjTEHwDMsYY0xCaNgy7LVjolBq36pBVw67VJxk9enSyjaVtVKen70HPoa3viKgPm7333nuTdS1Vz2Nny4KOwtI76iVF5NtYqKYdUa9j6zmgB/TMM88k63re+Tlsh0GvQ700NknsDHpuWRKHfgTLkuh7y0LBNTyXx8qQW203wbJII0eOTNbpaSlMaWBqgXoQTB/gPFYPiKWpGCLM60M/i69lOZ2c18G5qK0c6J+wE7GW1oqIuPPOO9v8HqKfzVqW999/f7KuHhF/uzhuXQV9KfqGOv4sJZQrZ0R0Dtdqtdi4caPDsI0xxjQnvgEZY4xpCL4BGWOMaQjbnQfE1tIs/XLHHXe0+V6Wz+ehd7R0SpXYeZY7KWt/3VE0NyGi3jNRryDnE0RE3flQX4f6LjVwXS/LN9LSJPRTmC9FvVy/h/o+tfaOluZhHhDPJfdZx4aeCdtqq6eovl9EvWaf80HK0GuA7SU4TtpumeeqM/O2SvtrXu865mzXzePRlhedKc/EVuW8XrQETRlaOom/R7xGc+TatkTk/SR6Quo58reLn1u1Hbw9IGOMMU2Jb0DGGGMagm9AxhhjGkLT1oLr06dPoZtrjTCWgc/pr8w7oeZKHVs/i14G61upLs+S+Kzlpdp6lRbQZXAfVZ9ljgr1cR1H5n8w54PejXpC9Abog2ieCn0ceii5OmDUtKnhaw4C2wwwt0HHhttY206/t0x35zxgHTCFnsShhx5aLNMjYWuNm266qVimb3PIIYck66yHqC2vc/M0Ij1fXZVzxs+NyHtAbE+u85btF5hXpp4K21tXqTX461//Ollnawedb/TnOAf0XDJvqaz9uNY4LGuznYOeoo4x6yjyN7OqB1SGn4CMMcY0BN+AjDHGNITtIgxbw/goB7EdA8uWKJSaKN+NGjWqWGZIJKWlXFn4zsCQxVyJf3Yy1fdS1sg9srNc/oMPPpis5yQSou0jIlLphnIdpTJtFUCZYPLkyck6S5xouC67ajKUV7tsUgahjKblaKqGb2vpHkozHAsNveZcY3kalUwYbk9pifKdSlEsuZQr0cJxoTSjcl6V+cJ94jylHFkl/FvPO0OPOcbbO5TjNWWAv1VVwvjHjBmTrLP0k14/On9aWlpi/fr1DsM2xhjTnPgGZIwxpiH4BmSMMaYhNK0H1K1bt0KDrlIiR8ONtRRHRH3IbZXSF+39zoh6vVlDr6nJl4VlaysEhq/ye1RrXbZsWbu/J1fCJKLeZ6jif2mpei3Z39r3agjx4sWLk20MW+Z71T+ih9iVoe85uE/qTzKMmS0wFi5c2Obn5nxBejMMV6fvqZ4KW4xUgV5fWQhxe+G4MIyZYdlKZ0r8ML2jSsizemX8OaVPq3Nc52xE+Rjqb0FZGLl+L0tT8VrS3xWGYXNOs/xRmZ9kD8gYY0xT4huQMcaYhuAbkDHGmIbQtB5Qe8m1i+4MWrIkol6vrdJuOQc1VubGaLkX5m10Bi39Qi+AOjDXH3rooWKZLQq6qkw/ywHR72LewyOPPFIss0QOS7YsWbKkWKYnR59Q81s4tzhPeS51jlA7z5VR4rFyfml+Eb0l+qXM49D8o7KyKqrd5/LRyuD54DVbpcyPXpe8JnPwXDGPif6Rzvl777233d9DeKyaf8RtPO/0TPVa6+qSOG0xevToZJ3fq55cay0g7AEZY4xpSnwDMsYY0xC2OwmOYacqR0SkYZtl8hCrAWuZnzJpT7+XUgAf5/V7WJX3vvvuS9ZZFicXmslw6VyIKtHw3DfeeKPNbRH1obAqLTHEM0eVrpSUI1gmhuGt+vjPMHGeSw3PZaiuho0TylvcJ4b25spCMaRepSWGtrL68ooVK4plzjV+J8+twmuJ54cy6LZCpUtK0BwnHXMea1lIukKJnZKoVtbuTOVpovOac5y/Zaxk3tbnRNTPRT0eymZl17fCjgLaPTUiYt68eW2+N8ISnDHGmCbFNyBjjDENwTcgY4wxDWG78IC0VD2pUladGj39GPWI6K9QD9cwToaZMnSRXpNCj4cehIZiUtOuUloo97llMAxbfRNq9PTdNPyYIc/UhrUcCv06lmhZu3Ztsq7lUOjV5LwAQh9Bu7ZSOz/uuOOSdc7ba6+9tlhmt07OCR03lrwnGkKvIeUR9eHFHOMqZYn0ezTMPaLec9BOshxDtlDpKjin6edVaTvAc8d2GTn03JV5SwrTKhjGz31SP4++M0sJDR48uFjOeZER6e9XlbJnEWmYf2ulnewBGWOMaUp8AzLGGNMQfAMyxhjTEJrWA9ppp50KXV/1ZUJ9v0o74AkTJiTrbAHQCBjfr1oufZuyckEKW1rruG3ZsiXZRh+B7RhWrVpVLNPLoE+i+VNlZXq0fYG2zY6o5mEx74d5HLncmBwsa3Pqqacm6xy3BQsWFMscY+ryOo6c77l8EHpjq1evTtZPPPHEZF3b0Ov+lUE/gsejHgQ9UfoKnOPqMWoL9Ih6T1G/p+ynS3OkRo0alWxjeR2eOx5fDm2ZzmOnR6p+GMelzJ/TOUK/l2Oh/jBfy3wv9eyYm1fFP1WvvlarxYYNG+wBGWOMaU58AzLGGNMQmlaC2xZQ2mAYcxX5LgdDL7VcDR/1eawsI6MyAvePUlOuejHlI91HPnbzEf3RRx9N1rV0R5Vw7rJSPLkK42USiYa+l4Ux5z6XY6qyB8PI+V6WKdGwVO4vLzuVU5kuwO62ep4pwd18883JOkOiNfWA2yjf6fGw5A9fmyvbU9ZxV4+hrOL7pEmTimV2kWW6hl7fuTDliPycYTV1no9t9RNK2ZCh8Dn0t+7dKqk0ZMiQYrmlpSXWrVtnCc4YY0xz4huQMcaYhuAbkDHGmIbQs/wljWGfffYpQhpVH6f2TA/l6aefbvMz6c1Qb86FLvK9qo9T+3/ggQeSdQ1lZFgjPR+S86Vy3V/LvCXVlzWMNKLeg2D4dJXOrFpen14A/aNp06a1ua2s1UEV30ehN8PSTjqO3Cd6igxB1/JAPFf0j3Jjyvmlob78XIbbs9S+zqeyUjU6V1lmhceuHhH9O76WHpB6BHytlkKKyPsgPJe5ViEM787BNic5z6czvg3Pc5X2KkS9Pl47PJe5UGt6xwwd1/nFrr/twU9AxhhjGoJvQMYYYxqCb0DGGGMawnaRB6Q6MHVRxrirX5QrTdMaqgt3poR8Lo+prI1zlX1mC1/NWaEOzzyOGTNmFMvMJSFsGz537tximdOH6+qlMffijjvuaPM76UuxPH6unXeV1sxljB8/vlimP8QxZWuEdwOOEzV6zmM9B/RTSa5MP+etwpwblmvqaCmkMnKtKOgtcf/ZtrpK7swBBxxQLK9ZsybZxjmS83SrlNYqQ8sdcVxyx0avmx5i1duF84CMMcY0Jb4BGWOMaQi+ARljjGkITZsHpGgeiuqtEfWaa2d0U9XLNYcjot6/UF+BuUesw6a6KT+XbcGpl+Z0etZW07bPQ4cOTbaxBL56UfRXqBEzR0W3c/8GDhzY5j6VtXnQsaF2zjL3/CytQzV//vzIoePGfBD6PHremdPBNg/aTiIizRsqK++vviHPB3Nh9HqgFzZlypRknf4d50wO9X04/9k2QSmr/ZaDviavF/WP6ONwnNS/oA9IX4pzvoofTM9L6UrPR/e5zEej19xe6IWxvh7HSee1XvstLS3ZObIVPwEZY4xpCL4BGWOMaQhNG4a93377FbKLPhZq98SIeqlGw1JZnoWP6HyczJXx0ZIyEfkyH9sKdldkaGzusZz7ryVDKIURPpazW2lH4eP91KlTi+Ubb7wx+16OhR4fzzPPq5agYcmlKrA9Q1lHyxw6b3lJUj5txNz72Mc+lqzfddddybpKTZSRKUlXIfdZ3MY5oa/ltTF8+PBkfe3atcl67mcxlwJAKF3myigxHaJKyasqdOW8VbRtxZYtW+Lhhx92GLYxxpjmxDcgY4wxDcE3IGOMMQ2hacOwVbfXkuDUX1mWPwdDL6mla8gnvSb6Rao3M7SVYbS5VtllaAkRegG5EE+WhWcpEtWiGeLM0GQeu3os3AeGG2toJsv/M+T5nnvuiY6i38v2EWxFkdPwy8L8lbJy+drem62+tZxRRDqmK1asSLZpu/GIvAdU5r+oJ1EWJqtzhm3Z6bfonKenWGXOcy7yvVpmidty30NfkGPM3wYdR84flqvRNi5s2cFx0s9l+HZZqH6OXKg+qVKair8b/A068MADi2X1U+lPt4WfgIwxxjQE34CMMcY0BN+AjDHGNISm9YAUtpDtKGX5Exr7X/ZajdnP6a0Rqe5Lv4hlS6hFay5AWYn4ESNGFMvMa+B7Ne+EJU1YOoXlglSnL9OttYwJ835Y4kTzEehhMReJ2rp6LCyZQ89HzwG35Twfttym/8WcIi0Fw/NK6EkoLPmTgx4Kyfk+LIOj51Z9joj6c6ftpNkOowz1Cbn/nDP6PYsWLUq20f/S3J8yT4J5Qlrqqay9vZ53jlPOo+Y4lfktufI6ud8g+tlV8n7oP/I61BYkHfGwttkT0I9+9KPYf//9Y8cdd4zJkyfHwoULt9VXGWOM2Q7ZJjegf//3f49Zs2bF1772tXjggQfikEMOieOPP77uLwdjjDHvX7ZJKZ7JkyfH4YcfHj/84Q8j4u3H30GDBsWXv/zlOO+887Lvba0jaq4zI9EyE2Xhhtyek7hyXTYpB40dOzZZX7BgQXY/FIZAd6Yza3spC91laLJu5+M95VKVJFgCpIq0xDHWTqURqXSmFbirQqlM1ymxlUmiKi1xnw499NBkXctG8ZLkPFWptUoaAmF5JsotWnLpiCOOSLatWrWqzX187rnnOrxPJDc32Q2WYf46Z1iyq2yOaFmZhx56KPtavWY5R6pUpWb4epU/2KvIdRw3nW9ldgLRKvT629zS0hJPP/30u1+K580334zFixfH9OnT//Ql3bvH9OnT60rDR7x9wjZt2pT8Z4wx5r1Pl9+AXnjhhdiyZUtiFka8bR629hfv7Nmzo2/fvsV/TPIyxhjz3qThUXDnn39+zJo1q1h/6aWX6irIVlEJ9bVl7+vo55Ztq5Jp3Jl96irKvpOyp66XRb7oZ7c3O7rscyLqx7gzn537nirzieT2ieNW5Xs6kzGvcP84profZePdVeNPcmNR9p2dmXtVxrgzc0Tpyuujyvd01T63tlz22V3uAb355pux0047xX/8x3/EKaecUvz7GWecERs3boxf//rX2fc/9dRTfgoyxpj3AOvWratLX1C6/AmoV69eMWHChJgzZ05xA2ppaYk5c+bEmWeeWfr+gQMHxrp166JWq8XgwYNj3bp1WRPr/c6mTZti0KBBHqcSPE7tw+PUPjxOeWq1Wrz88stJIE5rbBMJbtasWXHGGWfExIkTY9KkSXHppZfGq6++Gp/73OdK39u9e/f4wAc+UAQj7Lbbbj7B7cDj1D48Tu3D49Q+PE5tw2jm1tgmN6BPfepT8fzzz8eFF14Y69evj0MPPTR++9vf1gUmGGOMef+yzYIQzjzzzHZJbsYYY96fNG0x0t69e8fXvva1utpKJsXj1D48Tu3D49Q+PE5dwzaphGCMMcaU0bRPQMYYY97b+AZkjDGmIfgGZIwxpiH4BmSMMaYhNO0NyA3t/sTs2bPj8MMPj1133TX69+8fp5xySl0HzTfeeCNmzpwZ/fr1i1122SVOO+20Li2Lvz3yrW99K7p16xZnn3128W8ep7d5+umn4zOf+Uz069cv+vTpE2PHjk06jNZqtbjwwgtj3333jT59+sT06dPrWjC819myZUtccMEFMXTo0OjTp08MGzYsvvGNb9TVfnu/j1OnqDUh1157ba1Xr161//f//l/t4Ycfrv3P//k/a7vvvnvtueeea/SuNYTjjz++duWVV9aWL19eW7JkSe2EE06oDR48uPbKK68Ur/niF79YGzRoUG3OnDm1RYsW1Y444oja1KlTG7jXjWXhwoW1/fffvzZu3LjaWWedVfy7x6lWe/HFF2tDhgypffazn60tWLCgtmbNmtott9xSe/zxx4vXfOtb36r17du39qtf/aq2dOnS2kknnVQbOnRo7fXXX2/gnr+7XHLJJbV+/frVbrzxxtratWtrv/jFL2q77LJL7Xvf+17xGo9T52jKG9CkSZNqM2fOLNa3bNlSGzhwYG327NkN3Kvm4fe//30tImrz5s2r1Wq12saNG2s77LBD7Re/+EXxmkcffbQWEbX77ruvUbvZMF5++eXa8OHDa7fddlvt6KOPLm5AHqe3+Zu/+ZvakUce2eb2lpaW2j777FP7h3/4h+LfNm7cWOvdu3ftZz/72buxi03BiSeeWPsf/+N/JP/28Y9/vHb66afXajWPU1fQdBJc1YZ270deeumliIjYc889IyJi8eLF8dZbbyVjNmLEiBg8ePD7csxmzpwZJ554YjIeER6nrVx//fUxceLE+MQnPhH9+/ePww47LK644opi+9q1a2P9+vXJOPXt2zcmT578vhqnqVOnxpw5c2LlypUREbF06dK4++67Y8aMGRHhceoKGt4PiOQa2nWm/fB7hZaWljj77LNj2rRpMWbMmIh4u7V1r1696toOt9UE8L3MtddeGw888EDcf//9dds8Tm+zZs2auOyyy2LWrFnxf/7P/4n7778/vvKVr0SvXr3ijDPOKMaivU0l36ucd955sWnTphgxYkT06NEjtmzZEpdcckmcfvrpEREepy6g6W5AJs/MmTNj+fLlcffddzd6V5qOdevWxVlnnRW33XZb7Ljjjo3enaalpaUlJk6cGN/85jcjIuKwww6L5cuXx+WXXx5nnHFGg/euefj5z38eV199dVxzzTUxevToWLJkSZx99tkxcOBAj1MX0XQS3F577RU9evSoi0x67rnnYp999mnQXjUHZ555Ztx4441xxx13JE2e9tlnn3jzzTdj48aNyevfb2O2ePHi+P3vfx/jx4+Pnj17Rs+ePWPevHnx/e9/P3r27BkDBgzwOEXEvvvuG6NGjUr+beTIkfHkk09GRBRj8X6/Bv/6r/86zjvvvPj0pz8dY8eOjf/23/5bnHPOOTF79uyI8Dh1BU13A9KGdlvZ2tBuypQpDdyzxlGr1eLMM8+M6667Lm6//fYYOnRosn3ChAmxww47JGO2YsWKePLJJ99XY3bsscfGsmXLYsmSJcV/EydOjNNPP71Y9jhFTJs2rS6Mf+XKlTFkyJCIiBg6dGjss88+yTht2rQpFixY8L4ap9deey26d09/Inv06FG0m/Y4dQGNjoJojWuvvbbWu3fv2k9/+tPaI488UvvCF75Q23333Wvr169v9K41hC996Uu1vn371ubOnVt79tlni/9ee+214jVf/OIXa4MHD67dfvvttUWLFtWmTJlSmzJlSgP3ujnQKLhazeNUq70dot6zZ8/aJZdcUlu1alXt6quvru200061f/u3fyte861vfau2++67137961/XHnroodrJJ5/8vgsvPuOMM2r77bdfEYb9n//5n7W99tqrdu655xav8Th1jqa8AdVqtdoPfvCD2uDBg2u9evWqTZo0qTZ//vxG71LDiIhW/7vyyiuL17z++uu1//W//ldtjz32qO200061U089tfbss882bqebBN6APE5vc8MNN9TGjBlT6927d23EiBG1n/zkJ8n2lpaW2gUXXFAbMGBArXfv3rVjjz22tmLFigbtbWPYtGlT7ayzzqoNHjy4tuOOO9YOOOCA2le/+tXa5s2bi9d4nDqH2zEYY4xpCE3nARljjHl/4BuQMcaYhuAbkDHGmIbgG5AxxpiG4BuQMcaYhuAbkDHGmIbgG5AxxpiG4BuQMcaYhuAbkDHGmIbgG5AxxpiG4BuQMcaYhuAbkDHGmIbw/wNekZ4y6QNaaQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "norm = ImageNormalize(img_cut, stretch=AsinhStretch(),\n",
    "       vmin=0, vmax=np.percentile(img_cut, 99.5))\n",
    "ax.imshow(img_cut, origin='lower', cmap='gray', norm=norm)\n",
    "plot_pred_boundaries(edgeCoordsShift, ax, img_cut.shape,\n",
    "                         color='red',linewidth=2,alpha=1.0)\n",
    "ax.text(0.7,0.8,r'$I_{\\mathrm{E}}$ Band', transform=ax.transAxes, \n",
    "            color='white', fontsize=12,fontweight='bold')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4acb4085",
   "metadata": {},
   "source": [
    "## 3. Convert `gmm` to PDF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dd8dac0",
   "metadata": {},
   "source": [
    "The `gmm` column in the catalog provides the parameters required to reconstruct the full posterior photo-$z$ PDF. Given a photo-$z$ grid and the corresponding `gmm` entry, the PDF can be reconstructed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2174192c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from modelTool import gmm2pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "f3ea315a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmm = tcat['gmm'][idx].tolist()\n",
    "zgrids = np.arange(-2, 6, 0.1) # the grid is sufficient to cover redshfit range\n",
    "pdf = gmm2pdf(gmm, zgrids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "89d0e0f5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAG4CAYAAABmaLwhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAARHBJREFUeJzt3Xl4lPW9/vF7ZpJMFkggBLJA2EShKpsgEZcKNRr4WSrtKSLaglSx5UirRmulp4KtPUatC1apuLHYVhGrxbYqLlHgoCwCpsW1BMGwJGELmWSALDPP749kJpkSIJPMzDPL+3Vdc0lmnnnmM+cUuPl8N4thGIYAAADCmNXsAgAAAE6HwAIAAMIegQUAAIQ9AgsAAAh7BBYAABD2CCwAACDsEVgAAEDYI7AAAICwF2d2AYHgdru1b98+de3aVRaLxexyAABAOxiGoZqaGuXk5MhqPXUPJSoCy759+5Sbm2t2GQAAoAN2796tPn36nPKaqAgsXbt2ldT0hVNTU02uBgAAtIfD4VBubq737/FTiYrA4hkGSk1NJbAAABBh2jOdg0m3AAAg7BFYAABA2COwAACAsEdgAQAAYY/AAgAAwh6BBQAAhD0CCwAACHsEFgAAEPYILAAAIOwRWAAAQNjzO7CsXbtWkyZNUk5OjiwWi1auXHnK66+//npZLJYTHuecc473mnvuueeE14cMGeL3lwEAANHJ78DidDo1fPhwLVy4sF3XP/bYYyovL/c+du/erfT0dE2ZMsXnunPOOcfnunXr1vlbGgAAiFJ+H344ceJETZw4sd3Xp6WlKS0tzfvzypUrVVVVpZkzZ/oWEhenrKwsf8sBEEA7Dzq18uO9mj3uDCXG28wuBwC8Qj6H5bnnnlN+fr769evn8/z27duVk5OjgQMH6rrrrlNZWdlJ71FXVyeHw+HzANB58177RI8Vb9fLm3ebXQoA+AhpYNm3b5/efPNN3XjjjT7P5+XlaenSpVq1apWefPJJ7dy5U5dccolqamravE9RUZG3c5OWlqbc3NxQlA9EtbpGlzbtPCxJ2nHAaXI1AOArpIFl2bJl6tatmyZPnuzz/MSJEzVlyhQNGzZMBQUFeuONN3TkyBGtWLGizfvMnTtX1dXV3sfu3fxrEOisrV8fUV2jW5JUdvioydUAgC+/57B0lGEYWrx4sX74wx8qISHhlNd269ZNZ511lkpLS9t83W63y263B6NMIGZ9uOOg99cEFgDhJmQdljVr1qi0tFQ33HDDaa+tra3Vjh07lJ2dHYLKAEjSB6UtgWX34aNyuw0TqwEAX34HltraWpWUlKikpESStHPnTpWUlHgnyc6dO1fTp08/4X3PPfec8vLydO65557w2h133KE1a9Zo165d+vDDD/Xd735XNptN06ZN87c8AB1Qc7xB/9xT7f25rtGt/TV1JlYEAL78HhLavHmzxo8f7/25sLBQkjRjxgwtXbpU5eXlJ6zwqa6u1iuvvKLHHnuszXvu2bNH06ZN06FDh9SzZ09dfPHF2rBhg3r27OlveQA64KNdh+VyG+qbnixDhnYfPqayw0eVlZZodmkAIKkDgWXcuHEyjJO3ipcuXXrCc2lpaTp69ORj4suXL/e3DAAB9EHpIUnSRYN6qOzwUW9gGTMg3eTKAKAJZwkB8M5fufCMDPVNT5EklR1iaTOA8BGyVUIAwtPB2jp9UdG059GFZ/TQnqpjklgpBCC80GEBYtz6HU3DQUOyuqpHF7v6pidLIrAACC8EFiDGefZfufCMDElSvx4EFgDhh8ACxLgPd7RMuJWk3OYOy8HaejnrGk2rCwBaI7AAMWxP1VF9feiobFaLd0VQWlK80pLiJUm7q+iyAAgPBBYghn3YvJx5eJ80dU2M9z7vGRb6+hCBBUB4ILAAMeyD5vkrFw3K8HneMyy0m3ksAMIEgQWIUYZheOeveCbcerBSCEC4IbAAMap0f60O1NTJHmfVef26+bzWL50hIQDhhcACxCjP7rbn90+XPc7m81pfhoQAhBkCCxCjPvAMBzUvZ26tb/Ok2z1Vx+Ryn/zsMAAIFQILEIMaXW5t+Kp5/5X/mL8iSdlpSYqzWlTvcqvCcTzU5QHACQgsQAz6ZJ9DNccblZoYp3N7p53wus1qUZ/uSZKkMuaxAAgDBBYgBn1e7pAkndevu2xWS5vX9O3RdGoz81gAhAMCCxCDDjvrJUk9u9hPek3f9OYOC4EFQBggsAAx6MjRpsDSPSXhpNd4Vgp9TWABEAYILEAMqjraIEneM4Pa0je9aUiIDguAcEBgAWKQt8OSfPoOC3NYAIQDAgsQgzwdlu7Jp+iwNO/FcthZr5rjDSGpCwBOhsACxKCq5g5Lt1N0WLrY49SjeY4Lw0IAzEZgAWLQEU+HJeXkHRaJU5sBhA8CCxBj3G6jXXNYpFYrhdg8DoDJCCxAjKmpa5TneKBup5jDIkn9muexMCQEwGwEFiDGeLoryQm2E05p/k+eISECCwCzEViAGNOyQujUw0FSy5AQgQWA2QgsQIxpWSF06uEgqWVIaG/VMTW63EGtCwBOhcACxJgjfgSWzK6JSrBZ1eg2VF59PNilAcBJEViAGFPlbBoSOtUeLB5Wq0V9OAQRQBggsAAxpmVJ8+k7LJLUj3ksAMIAgQWIMf5MupWYeAsgPBBYgBjTnm35W/MubWbzOAAmIrAAMab62OkPPmytX48USXRYAJiLwALEmKp2bsvvwZAQgHBAYAFiTMsqofZ1WHKbVwlVH2tQdfP8FwAINQILEGPae/ChR3JCnFIT4yRJB511QasLAE6FwALEkPpGt5z1LkntDyySlJrU1I1xHKPDAsAcfgeWtWvXatKkScrJyZHFYtHKlStPef3q1atlsVhOeFRUVPhct3DhQvXv31+JiYnKy8vTpk2b/C0NwGl4uitWi9S1uWvSHl0TmwPL8cag1AUAp+N3YHE6nRo+fLgWLlzo1/u+/PJLlZeXex+9evXyvvbSSy+psLBQ8+fP19atWzV8+HAVFBRo//79/pYH4BQ8e7CkJcXLarW0+32eIaGa43RYAJij/f/EajZx4kRNnDjR7w/q1auXunXr1uZrjzzyiGbNmqWZM2dKkhYtWqTXX39dixcv1l133eX3ZwFom78rhDxahoTosAAwR8jmsIwYMULZ2dm6/PLL9cEHH3ifr6+v15YtW5Sfn99SlNWq/Px8rV+/vs171dXVyeFw+DwAnJ4/Bx+25hk+ctBhAWCSoAeW7OxsLVq0SK+88opeeeUV5ebmaty4cdq6dask6eDBg3K5XMrMzPR5X2Zm5gnzXDyKioqUlpbmfeTm5gb7awBRwd9t+T1SE5l0C8Bcfg8J+Wvw4MEaPHiw9+cLL7xQO3bs0KOPPqo//vGPHbrn3LlzVVhY6P3Z4XAQWoB2OHK0/Sc1t+YZEqph0i0AkwQ9sLRlzJgxWrdunSQpIyNDNptNlZWVPtdUVlYqKyurzffb7XbZ7fag1wlEG39PavZIZUgIgMlM2YelpKRE2dnZkqSEhASNGjVKxcXF3tfdbreKi4s1duxYM8oDopZ30m0KQ0IAIovfHZba2lqVlpZ6f965c6dKSkqUnp6uvn37au7cudq7d6+ef/55SdKCBQs0YMAAnXPOOTp+/LieffZZvffee3r77be99ygsLNSMGTM0evRojRkzRgsWLJDT6fSuGgIQGFVH/duW3yM1ybOsmSEhAObwO7Bs3rxZ48eP9/7smUsyY8YMLV26VOXl5SorK/O+Xl9fr9tvv1179+5VcnKyhg0bpnfffdfnHlOnTtWBAwc0b948VVRUaMSIEVq1atUJE3EBdI6/2/J7tGwcR4cFgDkshmEYZhfRWQ6HQ2lpaaqurlZqaqrZ5QBhK/+RNSrdX6sXbszThYMy2v2+bXuqNemJdcpKTdSGX14WxAoBxBJ//v7mLCEghrTsw+LvKiF2ugVgLgILECMMw/Aua+6e4u/GcU3XO+tdanS5A14bAJwOgQWIETV1jWp0N40A+z+HpWW6GxNvAZiBwALEiCPOpu5KYrxVifE2v94bb7MqOaHpPQQWAGYgsAAx4sixjq0Q8uA8IQBmIrAAMaKqg9vye7B5HAAzEViAGNHRbfk9POcJORgSAmACAgsQI6qcDAkBiFwEFiBGeIaE0jraYWFICICJCCxAjOj8kBDnCQEwD4EFiBGeDkvHh4Q4TwiAeQgsQIyo6uC2/B4tQ0J0WACEHoEFiBHebfk7OSREhwWAGQgsQIzobIfFMyTEAYgAzEBgAWJEdWc7LJ5lzQwJATABgQWIAQ0ut2rqmoJGRyfdtmwcR4cFQOgRWIAY4Jm/YrG0BA9/pXqHhOiwAAg9AgsQAzx7sKQlxctmtXToHp4hoZrjDXK7jYDVBgDtQWABYoD34MMOdlekls6M25Cc9XRZAIQWgQWIAZ1dISRJ9jirEmxNf2QwLAQg1AgsQAzo7Lb8kmSxWDgAEYBpCCxADOjstvwe3pVCLG0GEGIEFiAGBGJISPKdeAsAoURgAWJAZzeN8+AARABmIbAAMcDbYUnp7JAQu90CMAeBBYgBVQHqsKRynhAAkxBYgBjQskqocx2WllVCdFgAhBaBBYgB3o3jAtRhcRyjwwIgtAgsQJQzDMPbYen0KqEkzhMCYA4CCxDlnPUuNbiazv7p/CohNo4DYA4CCxDlqpxN3ZWEOKuS4m2duhdDQgDMQmABotyRViuELJaOndTs4d3pliEhACFGYAGiXFWAVghJLUNCLGsGEGoEFiDKHTkWmBVCku9ZQoZhdPp+ANBeBBYgygVqDxap5SyhepdbdY3uTt8PANqLwAJEuSqnp8PS+cCSkhAnzzQYVgoBCCUCCxDlWuawdH5IyGq1qKud84QAhJ7fgWXt2rWaNGmScnJyZLFYtHLlylNe/+qrr+ryyy9Xz549lZqaqrFjx+qtt97yueaee+6RxWLxeQwZMsTf0gC0oWXTuM4HFqn1SiE6LABCx+/A4nQ6NXz4cC1cuLBd169du1aXX3653njjDW3ZskXjx4/XpEmT9PHHH/tcd84556i8vNz7WLdunb+lAWhDy7b8nR8SkqSuiex2CyD04vx9w8SJEzVx4sR2X79gwQKfn++77z699tpr+vvf/66RI0e2FBIXp6ysLH/LAXAagZx0K7VMvGXzOAChFPI5LG63WzU1NUpPT/d5fvv27crJydHAgQN13XXXqays7KT3qKurk8Ph8HkAaFtVq43jAoEhIQBmCHlgeeihh1RbW6urr77a+1xeXp6WLl2qVatW6cknn9TOnTt1ySWXqKamps17FBUVKS0tzfvIzc0NVflAxKkK0MGHHi2bxzEkBCB0QhpYXnjhBf3617/WihUr1KtXL+/zEydO1JQpUzRs2DAVFBTojTfe0JEjR7RixYo27zN37lxVV1d7H7t37w7VVwAiSqPL7Q0WAeuwcJ4QABP4PYelo5YvX64bb7xRL7/8svLz8095bbdu3XTWWWeptLS0zdftdrvsdnswygSiSnWrUJGWxJAQgMgVkg7Liy++qJkzZ+rFF1/UlVdeedrra2trtWPHDmVnZ4egOiB6eQ4p7GKPU5wtML/dUxkSAmACvzsstbW1Pp2PnTt3qqSkROnp6erbt6/mzp2rvXv36vnnn5fUNAw0Y8YMPfbYY8rLy1NFRYUkKSkpSWlpaZKkO+64Q5MmTVK/fv20b98+zZ8/XzabTdOmTQvEdwRilueQQk/ICASGhACYwe9/cm3evFkjR470LkkuLCzUyJEjNW/ePElSeXm5zwqfp59+Wo2Njbr55puVnZ3tfdxyyy3ea/bs2aNp06Zp8ODBuvrqq9WjRw9t2LBBPXv27Oz3A2KaZzdaz94pgZCa1LysmQ4LgBDy+59d48aNO+UprUuXLvX5efXq1ae95/Lly/0tA0A7eDosXQPYYWnZOI4OC4DQ4SwhIIp5JsamBmjCrdR6SIgOC4DQIbAAUcwzMTaQHZaWISE6LABCh8ACRDHPPJPUQM5hab7X0XqXGlzugN0XAE6FwAJEMc9KnkB2WLq0ulctE28BhAiBBYhiniGhQM5hibdZlZxgk8SwEIDQIbAAUcwRhFVCEhNvAYQegQWIYi0bxwWuwyK1PgCRDguA0CCwAFGsZeO4AHdYOE8IQIgRWIAoVlPnGRIKbIfFs9U/Q0IAQoXAAkQxT6BISwpsh8UTgOiwAAgVAgsQpQzDaLU1f4A7LJwnBCDECCxAlHLWu+RuPvYr0JNuObEZQKgRWIAo5emuxFktSowP7G/1lgMQ6bAACA0CCxClWm8aZ7FYAnpvzhMCEGoEFiBKBWNbfg+GhACEGoEFiFI1QTj40KNl4ziGhACEBoEFiFLB2pZfYuM4AKFHYAGilGfJMUNCAKIBgQWIUp4wEYwhIc9OtzV1jXJ71k4DQBARWIAoVePtsAQhsDQPCRmG5KxnHguA4COwAFHKM78kNcDb8kuSPc6qBJu1+XMILACCj8ACRKlgdlgsFot3bgzzWACEAoEFiFKenW5TgzDpVmoZFmJpM4BQILAAUapl47jAd1ikliBEhwVAKBBYgCjVsjV/cDosniDEXiwAQoHAAkQp76TbYHVYktjtFkDoEFiAKFUTxI3jJDaPAxBaBBYgCjW43Dpa75IUvA6Ld5UQQ0IAQoDAAkSh2lbDNF2C3GFhSAhAKBBYgCjk6XokJ9gUbwvOb3MOQAQQSgQWIAoFe/5K63s7jtFhARB8BBYgCgV7hVDre9fQYQEQAgQWIAp5uh4h6bAwhwVACBBYgCjk3ZY/KXgdlq5MugUQQgQWIAo5gnjwoYenw8KQEIBQILAAUcgTIkIxJFTX6FZ9oztonwMAEoEFiEqeOSzBnHTbxd4ShmrrGBYCEFx+B5a1a9dq0qRJysnJkcVi0cqVK0/7ntWrV+u8886T3W7XoEGDtHTp0hOuWbhwofr376/ExETl5eVp06ZN/pYGoFkoOixxNquSE2w+nwcAweJ3YHE6nRo+fLgWLlzYrut37typK6+8UuPHj1dJSYluvfVW3XjjjXrrrbe817z00ksqLCzU/PnztXXrVg0fPlwFBQXav3+/v+UBUKtlzUGcdCu1dFmYeAsg2Pz+59fEiRM1ceLEdl+/aNEiDRgwQA8//LAk6Rvf+IbWrVunRx99VAUFBZKkRx55RLNmzdLMmTO973n99de1ePFi3XXXXf6WCMQ8T4BIDWKHRWrq4OyvqSOwAAi6oM9hWb9+vfLz832eKygo0Pr16yVJ9fX12rJli881VqtV+fn53mv+U11dnRwOh88DQIuWwBLcDktXNo8DECJBDywVFRXKzMz0eS4zM1MOh0PHjh3TwYMH5XK52rymoqKizXsWFRUpLS3N+8jNzQ1a/UAkcoRgDkvr+9NhARBsEblKaO7cuaqurvY+du/ebXZJQFjxdliCPIfFE1hYJQQg2IL7zy9JWVlZqqys9HmusrJSqampSkpKks1mk81ma/OarKysNu9pt9tlt9uDVjMQyQzDkONYaDosLZNuGRICEFxB77CMHTtWxcXFPs+98847Gjt2rCQpISFBo0aN8rnG7XaruLjYew2A9jve4Faj25AU3J1uW9+fISEAweZ3YKmtrVVJSYlKSkokNS1bLikpUVlZmaSm4Zrp06d7r//JT36ir776Snfeeae++OIL/eEPf9CKFSt02223ea8pLCzUM888o2XLlunzzz/X7Nmz5XQ6vauGALSfZ/6K1SKlNO+TEizeOSwMCQEIMr/7xZs3b9b48eO9PxcWFkqSZsyYoaVLl6q8vNwbXiRpwIABev3113XbbbfpscceU58+ffTss896lzRL0tSpU3XgwAHNmzdPFRUVGjFihFatWnXCRFwAp9eyaVy8LBZLUD+LfVgAhIrfgWXcuHEyDOOkr7e1i+24ceP08ccfn/K+c+bM0Zw5c/wtB8B/qPZsy58U9Clq3mXTtcxhARBkEblKCMDJeTss9uDOX5FY1gwgdAgsQJRpWdIc/A5LFwILgBAhsABRxtFqDkuweT6DfVgABBuBBYgyodqWX2oZEnIwhwVAkBFYgCgTqk3jJKmrvWWnW7f75JPxAaCzCCxAlAnVtvxSy5CQYUhHG1xB/zwAsYvAAkQZz/BMagg6LInxVsVZm/Z6YXt+AMFEYAGijKfDEoohIYvFwkohACFBYAGijGcOSygm3UrsxQIgNAgsQJRp6bCEKLDYPQcgMiQEIHgILECU8c5hCcHGcVLL5nHsxQIgmAgsQJQJdYcllSEhACFAYAGiiMtteDsdoVglJLU+sZkhIQDBQ2ABokhtqy5HyOaweE9spsMCIHgILEAU8cxfSYy3KiEuNL+9W7bnJ7AACB4CCxBFQnnwoQf7sAAIBQILEEVCuWmcR8uJzcxhARA8BBYgioR607imz6LDAiD4CCxAFDGjw9KySojAAiB4CCxAFGnZNC50HZaWISECC4DgIbAAUcTT5QjVHixS67OEmMMCIHgILEAU8YSGUM5h8QwJsawZQDARWIAo4jgW+jksnnBU3+hWXaMrZJ8LILYQWIAoUlMX+jksXVqFI3a7BRAsBBYgipjRYbFZLUpOsEli4i2A4CGwAFHEM4elqz10HRap9cRbAguA4CCwAFHEM/E1lENCUuuJt6wUAhAcBBYging7LCEcEmr6PE5sBhBcBBYginjmsIS6w8KQEIBgI7AAUeJ4g0v1LrckMzosbB4HILgILECU8HQ3LBapS0KIA4ud7fkBBBeBBYgSngmvXexxslotIf1shoQABBuBBYgSLecIhXb+itSyeRzb8wMIFgILECUcx8xZIdT0mQwJAQguAgsQJczssDDpFkCwEViAKOGZw5KaZEKHxc4cFgDBRWABokTLpnFmdFjYOA5AcHUosCxcuFD9+/dXYmKi8vLytGnTppNeO27cOFkslhMeV155pfea66+//oTXJ0yY0JHSgJjl3TTOlDksDAkBCC6//2R76aWXVFhYqEWLFikvL08LFixQQUGBvvzyS/Xq1euE61999VXV19d7fz506JCGDx+uKVOm+Fw3YcIELVmyxPuz3W73tzQgppnZYfGsEqph0i2AIPG7w/LII49o1qxZmjlzps4++2wtWrRIycnJWrx4cZvXp6enKysry/t45513lJycfEJgsdvtPtd17969Y98IiFHeSbdmzGFpDiy1dY1yu42Qfz6A6OdXYKmvr9eWLVuUn5/fcgOrVfn5+Vq/fn277vHcc8/pmmuuUUpKis/zq1evVq9evTR48GDNnj1bhw4dOuk96urq5HA4fB5ArHOY2GHxrEwyDMlZT5cFQOD5FVgOHjwol8ulzMxMn+czMzNVUVFx2vdv2rRJn3zyiW688Uaf5ydMmKDnn39excXFeuCBB7RmzRpNnDhRLperzfsUFRUpLS3N+8jNzfXnawBRyWHismZ7nFVxzbvrshcLgGAIae/4ueee09ChQzVmzBif56+55hrvr4cOHaphw4bpjDPO0OrVq3XZZZedcJ+5c+eqsLDQ+7PD4SC0IOZVHzVvWbPFYlHXxDhVHW1QzfFGZaeFvAQAUc6vDktGRoZsNpsqKyt9nq+srFRWVtYp3+t0OrV8+XLdcMMNp/2cgQMHKiMjQ6WlpW2+brfblZqa6vMAYt3ho02T29NTEkz5/C6sFAIQRH4FloSEBI0aNUrFxcXe59xut4qLizV27NhTvvfll19WXV2dfvCDH5z2c/bs2aNDhw4pOzvbn/KAmGUYhqqc5gYWz4nNbB4HIBj8XiVUWFioZ555RsuWLdPnn3+u2bNny+l0aubMmZKk6dOna+7cuSe877nnntPkyZPVo0cPn+dra2v185//XBs2bNCuXbtUXFysq666SoMGDVJBQUEHvxYQWxzHG9XYvDqne7JJgYUTmwEEkd+D3VOnTtWBAwc0b948VVRUaMSIEVq1apV3Im5ZWZmsVt8c9OWXX2rdunV6++23T7ifzWbTv/71Ly1btkxHjhxRTk6OrrjiCt17773sxQK00+Hm7kpKgk2J8TZTaiCwAAimDs3OmzNnjubMmdPma6tXrz7hucGDB8sw2t6bISkpSW+99VZHygDQzBNY0ruY012RWp/YzBwWAIHHWUJAFPAGFpOGgyQ6LACCi8ACRAGzJ9xKUhdObAYQRAQWIAocag4s3U0MLJ4hIQILgGAgsABRoKp5D5YepgYW9mEBEDwEFiAKHKoNhw5LywGIABBoBBYgCoRXh4XAAiDwCCxAFPDOYTF1lZBnDgtDQgACj8ACRAHPKqEeJu7D4lklxJAQgGAgsABR4HBYdFiaAouDISEAQUBgASJcXaPL29XokWLecRaeIaH6RrfqGl2m1QEgOhFYgAhX5WyaM2KzWrxdDjN4hoQkqZYuC4AAI7AAEa71cJDVajGtDpvVopSEpoMXWSkEINAILECE854jlBJvciVSF5Y2AwgSAgsQ4Q4fNf8cIQ/v0mZObAYQYAQWIMIdrq2TFC6BhQ4LgOAgsAAR7vDRpm5GOAQW714sBBYAAUZgASLcYWdzh8XEPVg8UtntFkCQEFiACOdZ1hwOHRaGhAAEC4EFiHCHmjssZp7U7MH2/ACChcACRDhPh8XMXW49PKuE2J4fQKARWIAI5z2pOQz2YWkZEmIOC4DAIrAAEcwwDFU178MSDh0Wz8ZxDAkBCDQCCxDBHMca5XIbksKjw5LKpFsAQUJgASKYZ5fbLvY42eNsJlfTaqdbhoQABBiBBYhgh70rhMzvrkhsHAcgeAgsQAQ77N2Dxfz5KxL7sAAIHgILEMFadrkNjw6LZ0iotr5R7ua5NQAQCAQWIIKFa4fFMCRnPV0WAIFDYAEimLfDEiZzWOxxVsXbLJIYFgIQWAQWIIKFW4fFYrGwPT+AoCCwABEs3DosEkubAQQHgQWIYIePhleHRWqZx8J5QgACicACRLBw7LCwFwuAYCCwABGsKszmsEith4QILAACh8ACRKi6Rpd3Ymt6coLJ1bRI5cRmAEFAYAEi1GFn0zlCNqtFqUlxJlfTogu73QIIgg4FloULF6p///5KTExUXl6eNm3adNJrly5dKovF4vNITEz0ucYwDM2bN0/Z2dlKSkpSfn6+tm/f3pHSgJjhCSzdkxNksVhMrqZFt+ZuT1XzwYwAEAh+B5aXXnpJhYWFmj9/vrZu3arhw4eroKBA+/fvP+l7UlNTVV5e7n18/fXXPq8/+OCD+v3vf69FixZp48aNSklJUUFBgY4fP+7/NwJihCew9EgJn+EgSerZtWk+zYGaOpMrARBN/A4sjzzyiGbNmqWZM2fq7LPP1qJFi5ScnKzFixef9D0Wi0VZWVneR2Zmpvc1wzC0YMEC/epXv9JVV12lYcOG6fnnn9e+ffu0cuXKDn0pIBZ4OyxhtEJIknp2aQpQB2sJLAACx6/AUl9fry1btig/P7/lBlar8vPztX79+pO+r7a2Vv369VNubq6uuuoqffrpp97Xdu7cqYqKCp97pqWlKS8v75T3BGJdS4clfFYISa06LAQWAAHkV2A5ePCgXC6XT4dEkjIzM1VRUdHmewYPHqzFixfrtdde05/+9Ce53W5deOGF2rNnjyR53+fPPevq6uRwOHweQKypag4s6WE2JJTRpWVIyDA4sRlAYAR9ldDYsWM1ffp0jRgxQpdeeqleffVV9ezZU0899VSH71lUVKS0tDTvIzc3N4AVA5HhkHdIKDwDy/EGt5z1LpOrARAt/AosGRkZstlsqqys9Hm+srJSWVlZ7bpHfHy8Ro4cqdLSUknyvs+fe86dO1fV1dXex+7du/35GkBU8KzCCbdJtyn2OKUk2CQx8RZA4PgVWBISEjRq1CgVFxd7n3O73SouLtbYsWPbdQ+Xy6Vt27YpOztbkjRgwABlZWX53NPhcGjjxo0nvafdbldqaqrPA4g1h2rDs8MiSRnN81iYeAsgUPzebaqwsFAzZszQ6NGjNWbMGC1YsEBOp1MzZ86UJE2fPl29e/dWUVGRJOk3v/mNLrjgAg0aNEhHjhzR7373O3399de68cYbJTWtILr11lv129/+VmeeeaYGDBigu+++Wzk5OZo8eXLgvikQZcK1wyJJPbvY9fWho3RYAASM34Fl6tSpOnDggObNm6eKigqNGDFCq1at8k6aLSsrk9Xa0ripqqrSrFmzVFFRoe7du2vUqFH68MMPdfbZZ3uvufPOO+V0OnXTTTfpyJEjuvjii7Vq1aoTNpgD0KL1xnHhxjOPhQ4LgECxGFEwjd/hcCgtLU3V1dUMDyEmuN2GzvzVm3K5DW385WXKTA2vcH/3yk/0xw1f66ffGqTbrxhsdjkAwpQ/f39zlhAQgRzHG+RyN/1bo1tyeG0cJ/kubQaAQCCwABHIMxzU1R4ne5zN5GpO1JNJtwACjMACRKDDYboHi0dG8/b8dFgABAqBBYhAh8N0l1uPlg4LJzYDCAwCCxCBIiWwsD0/gEAhsAAR6PDR8A4snkm39S63HMcbTa4GQDQgsAAR6HBteAeWxHibuiY2bfPEPBYAgUBgASJQuHdYpKbdbiVWCgEIDAILEIG8c1jCcJdbj4yu7MUCIHAILEAEqgrzSbdSS4eFwAIgEAgsQAQ6FOb7sEhsHgcgsAgsQATydFjC8aRmDzaPAxBIBBYgwhxvcMlZ75JEhwVA7CCwABGmqnmFUJzVotTmpcPhyLt5HIEFQAAQWIAIc6i2Zf6KxWIxuZqT82wed7CG7fkBdB6BBYgwng5LOM9fkXyHhNxutucH0DkEFiDCeE9qDuM9WCSpR0pTYGl0GzpyrMHkagBEOgILEGG8m8Z1Ce/AkhBnVbfkeElMvAXQeQQWIMJUOpr+8s8I8yEhqWUeC0ubAXQWgQWIMF8dqJUkDchIMbmS0+M8IQCBQmABIkxpc2A5o1cXkys5Pc4TAhAoBBYggtQ3uvX1oaOSpEEREFi85wnRYQHQSQQWIIJ8fcgpl9tQSoJNWamJZpdzWj3psAAIEAILEEF2tBoOCudN4zw4TwhAoBBYgAhSur8psAzqGf7DQVLrzePY7RZA5xBYgAjiCSyRMOFWYlkzgMAhsAARZMcBp6TImHArSb2aOyyHnXVysT0/gE4gsAARwu02WuawRMiQUHpKgiwWyW207NALAB1BYAEiRLnjuI7WuxRntahfj2Szy2mXOJtV6c1nHrF5HIDOILAAEcIzf6V/RoribZHzW5elzQACIXL+1ANi3A7PhNue4b8lf2sZbM8PIAAILECE8GzJHykTbj3osAAIBAILECG8e7BEWGBh8zgAgUBgASLEVxG2QsijZfM4AguAjiOwABHgyNF6726xkRZYMjgAEUAAEFiACOAZDspJS1SKPc7kavzj7bDUsA8LgI4jsAARINK25G+NDguAQOhQYFm4cKH69++vxMRE5eXladOmTSe99plnntEll1yi7t27q3v37srPzz/h+uuvv14Wi8XnMWHChI6UBkSlSNvhtjVPh6XqaL0aXG6TqwEQqfwOLC+99JIKCws1f/58bd26VcOHD1dBQYH279/f5vWrV6/WtGnT9P7772v9+vXKzc3VFVdcob179/pcN2HCBJWXl3sfL774Yse+ERCFInWFkCR1T06QzWqRwfb8ADrB78DyyCOPaNasWZo5c6bOPvtsLVq0SMnJyVq8eHGb1//5z3/Wf//3f2vEiBEaMmSInn32WbndbhUXF/tcZ7fblZWV5X107969Y98IiEKRugeLJNmsFqWnsLQZQOf4FVjq6+u1ZcsW5efnt9zAalV+fr7Wr1/frnscPXpUDQ0NSk9P93l+9erV6tWrlwYPHqzZs2fr0KFDJ71HXV2dHA6HzwOIVscbXNpTdUxSZA4JSVJP5rEA6CS/AsvBgwflcrmUmZnp83xmZqYqKiradY9f/OIXysnJ8Qk9EyZM0PPPP6/i4mI98MADWrNmjSZOnCiXy9XmPYqKipSWluZ95Obm+vM1gIjy1QGnDENKS4r3bsIWaTLY7RZAJ4V0feT999+v5cuXa/Xq1UpMTPQ+f80113h/PXToUA0bNkxnnHGGVq9ercsuu+yE+8ydO1eFhYXenx0OB6EFUav1cJDFYjG5mo7pyXlCADrJrw5LRkaGbDabKisrfZ6vrKxUVlbWKd/70EMP6f7779fbb7+tYcOGnfLagQMHKiMjQ6WlpW2+brfblZqa6vMAolVphB562FpGV+awAOgcvwJLQkKCRo0a5TNh1jOBduzYsSd934MPPqh7771Xq1at0ujRo0/7OXv27NGhQ4eUnZ3tT3lAVNoRwRNuPVo6LKwSAtAxfq8SKiws1DPPPKNly5bp888/1+zZs+V0OjVz5kxJ0vTp0zV37lzv9Q888IDuvvtuLV68WP3791dFRYUqKipUW9v0h3Btba1+/vOfa8OGDdq1a5eKi4t11VVXadCgQSooKAjQ1wQi144IXtLs0XJi83GTKwEQqfyewzJ16lQdOHBA8+bNU0VFhUaMGKFVq1Z5J+KWlZXJam3JQU8++aTq6+v1/e9/3+c+8+fP1z333CObzaZ//etfWrZsmY4cOaKcnBxdccUVuvfee2W32zv59YDI5nIb+uqgU5I0qGdXk6vpODosADqrQ5Nu58yZozlz5rT52urVq31+3rVr1ynvlZSUpLfeeqsjZQBRb0/VUdU3upUQZ1Xv7klml9NhPVklBKCTOEsICGOeCbcDM1Jks0bmCiGp5Tyh6mMNqmtse7sCADgVAgsQxiJ5S/7W0pLiFW9rClyHGBYC0AEEFiCMtSxpjuzAYrVa1COlqcuyn2EhAB1AYAHCWDQsafYYkNG0j8zn5RylAcB/BBYgTBmGETVDQpI0ql/TgaZbvq4yuRIAkYjAAoSpA7V1chxvlMXS0p2IZJ7AspXAAqADCCxAmFr95QFJTWElMd5mcjWdN7JvN0nSVwedOsSZQgD8RGABwpBhGFrywS5J0pRR0XGwZ7fkBO/Q1tayI+YWAyDiEFiAMLRx52F9Xu5QYrxV08ZER2CRpNHMYwHQQQQWIAwt+WCnJOm7I/uoW3KCydUEznnewHLY5EoARBoCCxBmdh8+qnc+q5Qkzbyov7nFBJhn4u0/91SrvtFtcjUAIgmBBQgzz6/fJbchXXJmhs7KjNwDD9syMCNF3ZPjVd/o1qf7qs0uB0AEIbAAYcRZ16jlH+2WFH3dFUmyWCzsxwKgQwgsQBh5dese1Rxv1ICMFI07q5fZ5QTFeQQWAB1AYAHChNttaMmHuyRJM8b2kzWCT2c+lVF9mwLL5q+rZBiGydUAiBQEFiBMrN1+QF8dcKqrPU7fHx09S5n/0/DcboqzWnSgpk57qo6ZXQ6ACEFgAcKEd6O40bnqYo8zt5ggSoy36ZzeaZIYFgLQfgQWIAyU7q/Vmn8fkMUiXX9hf7PLCTrPsBCBBUB7EViAMLCsee7KZUMy1bdHsrnFhMDo/gQWAP6J3r4zcBoV1cd13bMbZLFY9N2RvTV5ZG/17pYU8jre/axSf974tSTpR1G4lLktnqXNX1Q4VFvXGNVDYAACgw4LYpLLbejWlz7WjgNOle6v1e/e+lIXP/Cern1mg/6yZY9q6xpDUseWr6s058WtchvSlFF9NPaMHiH5XLNlpiaqd7ckuQ2phIMQAbQDgQUx6Yn3SrXhq8NKTrBp3rfP1tiBPWQY0oc7DumOl/+p83/7rl4r2RvUGkr31+qGZR/peINb4wf31H3fGyqLJTqXMreFDeQA+IPAgpiz4atDeqz435Kk//3uufrRxQP04k0XaN0vxuvnBYM1MCNFxxpcuuuVbdp10BmUGiodxzVj8SYdOdqg4bndtPC68xRvi63fjp55LJs5CBFAO8TWn5CIeYed9bp1eYnchvRf5/XRd0f28b7Wp3uybh4/SO8WXqqxA3voWINLd7z8T7ncgd3czHG8QTMWb9LeI8c0ICNFi2eMVnJC7M3hOK95pVBJ2ZGA/98YQPQhsCBmGIahn7/8T1U4jmtgzxT95qpz2rzOarXowe8PUxd7nDZ/XaVn/++rgNVQ1+jSTc9v1hcVNcroYtfzPxqjHl3sAbt/JBmS1VXJCTbV1DVq+/4as8sBEOYILIgZiz/YpeIv9ishzqrHp41UyilWpuSmJ+vub39DkvTw2//Wvys7/xdqg8utn734sTZ8dVhd7HFaOvN85aZH/xLmk4mzWTWybzdJ0uZdzGMBcGoEFsSEbXuqdf+bn0uSfnXlN3ROTtpp33P16Fx9a0gv1bvcKlxRogaXu8Of3+hy69blJXrr00ol2Kxa9INROrf36WuIdp4N5LYy8RbAaRBYEPWON7h0y/KP1eAyVHBOpn54Qb92vc9isej+7w1VWlK8Ptnr0BPvlXbo811uQ7e//E+9vq1c8TaLFv3wPF18ZkaH7hVtPCc3b/jqkOobOx4IAUQ/Agui3qPv/FtfHXQqM9WuB/5rmF9Lh3ulJureyedKkp54v1Tb9lT79dkut6Gf/+Wfeq1kn+KsFv3hulH61pBMv+4RzUb3T1e35Hjtqz6u3xdvN7scAGGMwIKo9nFZlZ5pnjR733eHqltygt/3mDQsW1cOzZbLbahwRUm7N5Vzuw3NffVfenXrXtmsFj0+baQuP5uw0loXe5zu++5QSdIfVpdq8y6WOANoG4EFUauu0aU7//IvuQ1p8ogcXfaNjoUFi8Wieyefq4wudm3fX6srf/9/+rjs1HMuGlxu/c/KT7Ri8x5ZLdKCqSM0cWh2hz4/2v2/odn63nm95Tak21aUqOZ4g9klAQhDBBZErSfeK9X2/bXK6JKg+ZPaXsLcXukpCXpm+ijlpCXq60NH9f1F6/XYu9vV+B8TcRtdbq3YvFvfeni1XtxUJotFeuTqEZo0PKdTnx/t7vnOOerdLUm7Dx/Tb/7+mdnlAAhDBBZEpU/2VusPq3dIku696lx1T/F/KOg/jezbXW/e+k19Z3iOXG5Dj777b019eoPKDh1Vo8utv2zZo8seWaM7//Iv7T58TBldEvT7a0Zq8sjenf7saJeaGK9Hp46QxSK9vGWPVn1SbnZJAMKMxTCMiN9i0uFwKC0tTdXV1UpNTTW7HJisweXWVU98oM/KHfp/Q7P0h+tGBfwzVn68V3ev/EQ1zScNZ3RJ0K5DRyVJPVIS9ONLB+oHF/SLyR1sO+P+N7/QojU71D05Xm/d+k31Sk00uyQAQeTP3990WBB1nlqzQ5+VO9QtOV6//s65QfmMySN7641bLtH5/burtq5Ruw4dVffkeN01cYjW3jleN33zDMJKBxRefpbOzk5V1dEG/fwv/1IU/HsKQIDQYUFU+WRvtb73hw9V73JrwdQRQR+OcbkNvbDxa9W7DF1zfu4pd89F+2yvrNG3H1+nuka3po3J1S2XnaWsNDotQDTy5+9vAguigtttaMmHu/Tgqi9U1+hW/jd66Znpo/3acwXh4/n1uzTvtU8lSQk2q6aM7qOfXHpGTB9lAESjoA8JLVy4UP3791diYqLy8vK0adOmU17/8ssva8iQIUpMTNTQoUP1xhtv+LxuGIbmzZun7OxsJSUlKT8/X9u3s4kU2mdP1VFd++wG3fuPz1TX6NalZ/XUg98fTliJYNPH9tfzPxqjMf3TVe9y688byzTuodW6fcU/Vbq/1uzyAJjA7w7LSy+9pOnTp2vRokXKy8vTggUL9PLLL+vLL79Ur169Trj+ww8/1De/+U0VFRXp29/+tl544QU98MAD2rp1q849t2l+wQMPPKCioiItW7ZMAwYM0N13361t27bps88+U2Li6VvBdFhik2EYennLHv3m75+ptq5RSfE2/erb39C1Y/oSVqLIxq8O6Yn3S/V/2w96n8tOS9SwPmka1qebhvfppqF90pSWFG9ilQA6IqhDQnl5eTr//PP1xBNPSJLcbrdyc3P105/+VHfdddcJ10+dOlVOp1P/+Mc/vM9dcMEFGjFihBYtWiTDMJSTk6Pbb79dd9xxhySpurpamZmZWrp0qa655pqAfmFErkaXW3uqjql0f612HKjVutKD3r/ERvXrroenDFf/jBSTq0Sw/HP3ET3xfqmKP6+Uu40/tXLSEpXdLUnZaYnKTktUVlqSctIS1T0lQamJ8eqaGKfUpHh1tcfJaiXQAuHAn7+//ZohWF9fry1btmju3Lne56xWq/Lz87V+/fo237N+/XoVFhb6PFdQUKCVK1dKknbu3KmKigrl5+d7X09LS1NeXp7Wr1/fZmCpq6tTXV2d92eHw+HP12i3Rpdbd726LSj3jnWtY7Ihw+f5BpdbjS5DjW63Gpr/e6CmTrsOHlX9f2zUFm+zqPDywbrpmwNl4y+hqDY8t5uemT5azrpGfbK3Wv/aU61/7jmif+45ot2Hj2lf9XHtqz5+2vtYLFKXhDjZ422yx1llj7fKHtf064Q4q2wWi+JsFtmsFtksTf+1WiyyWiWLLLJYJKul6b8W7z0tTb9ufsKiE/+3SNMPkS7OatH9/zXMvM/35+KDBw/K5XIpM9N3i/PMzEx98cUXbb6noqKizesrKiq8r3ueO9k1/6moqEi//vWv/Sm9Q9yG9Jcte4L+OWg/e5xVA3t20Rk9U3RGzy66cli2zsrsanZZCKEUe5zyBvZQ3sAe3ucOO+u165BT5UeOq7z6mMqrj6uiuunXR441yHGsUTXHG1TX6JZhSDV1japp55lQAJokxFkjJ7CEi7lz5/p0bRwOh3JzcwP+OTarRb+YMCTg90WT1v/itLR6Ls5qVbzNojibVXFWi+JtVqUlx2tQzy7q3S2Jdj5OkJ6SoPSUBKnvqa873uBSzfFG1dY1qq7RpboGt+oa3T6/dhmGXO6mLp/bMNTgMmQYTX1At7v5v4a8e8R4uoWeTmFbg+wnG3eP/DWaiCU2k3du8yuwZGRkyGazqbKy0uf5yspKZWVltfmerKysU17v+W9lZaWys7N9rhkxYkSb97Tb7bLb7f6U3iE2q0Wzx50R9M8BEBqJ8TYlxtvUs2vw//wAEFh+5aWEhASNGjVKxcXF3ufcbreKi4s1duzYNt8zduxYn+sl6Z133vFeP2DAAGVlZflc43A4tHHjxpPeEwAAxBa/h4QKCws1Y8YMjR49WmPGjNGCBQvkdDo1c+ZMSdL06dPVu3dvFRUVSZJuueUWXXrppXr44Yd15ZVXavny5dq8ebOefvppSU2T1W699Vb99re/1Zlnnuld1pyTk6PJkycH7psCAICI5XdgmTp1qg4cOKB58+apoqJCI0aM0KpVq7yTZsvKymS1tjRuLrzwQr3wwgv61a9+pV/+8pc688wztXLlSu8eLJJ05513yul06qabbtKRI0d08cUXa9WqVe3agwUAAEQ/tuYHAACm4LRmAAAQVQgsAAAg7BFYAABA2COwAACAsEdgAQAAYY/AAgAAwh6BBQAAhD0CCwAACHsEFgAAEPb83po/HHk263U4HCZXAgAA2svz93Z7Nt2PisBSU1MjScrNzTW5EgAA4K+amhqlpaWd8pqoOEvI7XZr37596tq1qywWS0Dv7XA4lJubq927d0ftOUV8x+jAd4wOfMfoEQvfs7Pf0TAM1dTUKCcnx+fg5LZERYfFarWqT58+Qf2M1NTUqP0fnAffMTrwHaMD3zF6xML37Mx3PF1nxYNJtwAAIOwRWAAAQNgjsJyG3W7X/PnzZbfbzS4laPiO0YHvGB34jtEjFr5nKL9jVEy6BQAA0Y0OCwAACHsEFgAAEPYILAAAIOwRWAAAQNgjsLTTrl27dMMNN2jAgAFKSkrSGWecofnz56u+vt7s0gLqf//3f3XhhRcqOTlZ3bp1M7ucgFm4cKH69++vxMRE5eXladOmTWaXFDBr167VpEmTlJOTI4vFopUrV5pdUsAVFRXp/PPPV9euXdWrVy9NnjxZX375pdllBdSTTz6pYcOGeTfgGjt2rN58802zywqq+++/XxaLRbfeeqvZpQTMPffcI4vF4vMYMmSI2WUF3N69e/WDH/xAPXr0UFJSkoYOHarNmzcH9TMJLO30xRdfyO1266mnntKnn36qRx99VIsWLdIvf/lLs0sLqPr6ek2ZMkWzZ882u5SAeemll1RYWKj58+dr69atGj58uAoKCrR//36zSwsIp9Op4cOHa+HChWaXEjRr1qzRzTffrA0bNuidd95RQ0ODrrjiCjmdTrNLC5g+ffro/vvv15YtW7R582Z961vf0lVXXaVPP/3U7NKC4qOPPtJTTz2lYcOGmV1KwJ1zzjkqLy/3PtatW2d2SQFVVVWliy66SPHx8XrzzTf12Wef6eGHH1b37t2D+8EGOuzBBx80BgwYYHYZQbFkyRIjLS3N7DICYsyYMcbNN9/s/dnlchk5OTlGUVGRiVUFhyTjr3/9q9llBN3+/fsNScaaNWvMLiWounfvbjz77LNmlxFwNTU1xplnnmm88847xqWXXmrccsstZpcUMPPnzzeGDx9udhlB9Ytf/MK4+OKLQ/65dFg6obq6Wunp6WaXgVOor6/Xli1blJ+f733OarUqPz9f69evN7EydEZ1dbUkRe3vP5fLpeXLl8vpdGrs2LFmlxNwN998s6688kqf35fRZPv27crJydHAgQN13XXXqayszOySAupvf/ubRo8erSlTpqhXr14aOXKknnnmmaB/LoGlg0pLS/X444/rxz/+sdml4BQOHjwol8ulzMxMn+czMzNVUVFhUlXoDLfbrVtvvVUXXXSRzj33XLPLCaht27apS5custvt+slPfqK//vWvOvvss80uK6CWL1+urVu3qqioyOxSgiIvL09Lly7VqlWr9OSTT2rnzp265JJLVFNTY3ZpAfPVV1/pySef1Jlnnqm33npLs2fP1s9+9jMtW7YsqJ8b84HlrrvuOmGC1H8+vvjiC5/37N27VxMmTNCUKVM0a9Yskypvv458RyBc3Xzzzfrkk0+0fPlys0sJuMGDB6ukpEQbN27U7NmzNWPGDH322WdmlxUwu3fv1i233KI///nPSkxMNLucoJg4caKmTJmiYcOGqaCgQG+88YaOHDmiFStWmF1awLjdbp133nm67777NHLkSN10002aNWuWFi1aFNTPjQvq3SPA7bffruuvv/6U1wwcOND763379mn8+PG68MIL9fTTTwe5usDw9ztGk4yMDNlsNlVWVvo8X1lZqaysLJOqQkfNmTNH//jHP7R27Vr16dPH7HICLiEhQYMGDZIkjRo1Sh999JEee+wxPfXUUyZXFhhbtmzR/v37dd5553mfc7lcWrt2rZ544gnV1dXJZrOZWGHgdevWTWeddZZKS0vNLiVgsrOzT+j8feMb39Arr7wS1M+N+cDSs2dP9ezZs13X7t27V+PHj9eoUaO0ZMkSWa2R0aDy5ztGm4SEBI0aNUrFxcWaPHmypKZ/HRQXF2vOnDnmFod2MwxDP/3pT/XXv/5Vq1ev1oABA8wuKSTcbrfq6urMLiNgLrvsMm3bts3nuZkzZ2rIkCH6xS9+EXVhRZJqa2u1Y8cO/fCHPzS7lIC56KKLTthW4N///rf69esX1M+N+cDSXnv37tW4cePUr18/PfTQQzpw4ID3tWj6l3pZWZkOHz6ssrIyuVwulZSUSJIGDRqkLl26mFtcBxUWFmrGjBkaPXq0xowZowULFsjpdGrmzJlmlxYQtbW1Pv9627lzp0pKSpSenq6+ffuaWFng3HzzzXrhhRf02muvqWvXrt75R2lpaUpKSjK5usCYO3euJk6cqL59+6qmpkYvvPCCVq9erbfeesvs0gKma9euJ8w7SklJUY8ePaJmPtIdd9yhSZMmqV+/ftq3b5/mz58vm82madOmmV1awNx222268MILdd999+nqq6/Wpk2b9PTTTwd/1CHk65Ii1JIlSwxJbT6iyYwZM9r8ju+//77ZpXXK448/bvTt29dISEgwxowZY2zYsMHskgLm/fffb/P/ZzNmzDC7tIA52e+9JUuWmF1awPzoRz8y+vXrZyQkJBg9e/Y0LrvsMuPtt982u6ygi7ZlzVOnTjWys7ONhIQEo3fv3sbUqVON0tJSs8sKuL///e/Gueeea9jtdmPIkCHG008/HfTPtBiGYQQ3EgEAAHROZEzCAAAAMY3AAgAAwh6BBQAAhD0CCwAACHsEFgAAEPYILAAAIOwRWAAAQNgjsAAAgLBHYAEAAGGPwAIgaM4991xZLJaTPn7961+bXSKACMHhhwCCZtq0aWpoaPB5rr6+XgsWLNCxY8f0zW9+06TKAEQazhICEDJ1dXX63ve+pzfffFN/+MMf9JOf/MTskgBECIaEAITEsWPH9J3vfEerVq3SM8880+6wMm7cON1xxx1Brg5AuCOwAAi6o0ePatKkSXr33Xe1ZMkS3XDDDSH9fEIPEPmYwwIgqJxOp6688kqtW7dOf/zjH3XttdeaXRKACESHBUDQ1NTUaMKECfrggw/04osvthlWxo0bp5/97GeaNWuWunbtquzsbD3++OM+17hcLt12223q1q2bevfurQULFvi8fuzYMc2ZM0c9e/ZUYmKixo8fr23btkmSrr/+eq1Zs0YPP/ywd3USgMhDYAEQFNXV1briiiu0ceNGvfzyy5oyZcpJr12yZIl69eqlzZs36+6771ZhYaHeffdd7+uLFy9WZmamPvroI/3P//yPbrvtNm3dutX7+p133qnXXntNf/7zn7V582b16tVLEyZM0NGjR/XYY49p7Nixmj17tsrLy1VeXh7U7w0gOFglBCDgqqqqdMUVV2jbtm36y1/+om9/+9snvXbcuHGqrq7Wxx9/7H3u2muvVW1trf72t79p3Lhxslqteu+997yvn3XWWbrpppt0xx13qLa2Vunp6frTn/6kq6++WlLTnJm+ffuqqKhIs2bN0rhx4zR69Gg99NBDwfvSAIKKOSwAAu4HP/iBNm/erEsvvVSbN2/W5s2bfV7v1auX/vu//9v78/nnn+/z+gUXXOAzLDR06FCf17Ozs7V//35J0o4dO9TQ0KCLLrrI+3pycrJGjhypzz//PGDfCYC5CCwAAsrtdmvt2rWSpDVr1mjNmjUnXHPVVVf5BJbTzSuJj4/3+dliscjtdp/wHIDoxRwWAAFltVpVU1MjwzBO+li5cqXPez766COfnzdu3KghQ4a06/POOOMMxcfH68MPP/Q+d+zYMX388cc6++yzJTUFHpfL1bkvBsBUdFgAmG779u26++679cMf/lDvvfeeVqxYoddff71d7+3SpYt+/OMf6/bbb1daWppycnL0m9/8Rna73bsqqX///tq0aZPKysqUnJysjIyMYH4dAEFAYAFgupkzZ2rv3r0677zzlJKSot/97ne64oor2v3+Bx98UIZh6Nprr1VNTY0uuOACvfnmm0pOTpYk3X777Zo+fboGDx6s48ePi7UGQORhlRAAU7GCB0B7MIcFAACEPQILAAAIewwJAQCAsEeHBQAAhD0CCwAACHsEFgAAEPYILAAAIOwRWAAAQNgjsAAAgLBHYAEAAGGPwAIAAMIegQUAAIQ9AgsAAAh7BBYAABD2CCwAACDs/X9QR+uW5lynCAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(zgrids, pdf)\n",
    "plt.xlabel(r'$z_{\\mathrm{phot}}$',fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "668d76cc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "euc",
   "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.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}