csiborgtools/notebooks/flow_mapping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using a calibrated flow model to predict $z_{\\rm cosmo}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright (C) 2024 Richard Stiskalek\n",
    "# This program is free software; you can redistribute it and/or modify it\n",
    "# under the terms of the GNU General Public License as published by the\n",
    "# Free Software Foundation; either version 3 of the License, or (at your\n",
    "# option) any later version.\n",
    "#\n",
    "# This program is distributed in the hope that it will be useful, but\n",
    "# WITHOUT ANY WARRANTY; without even the implied warranty of\n",
    "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General\n",
    "# Public License for more details.\n",
    "#\n",
    "# You should have received a copy of the GNU General Public License along\n",
    "# with this program; if not, write to the Free Software Foundation, Inc.,\n",
    "# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from h5py import File\n",
    "\n",
    "import csiborgtools\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_calibration(catalogue, simname, nsim, ksmooth):\n",
    "    fname = f\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/flow_samples_{catalogue}_{simname}_smooth_{ksmooth}.hdf5\"  # noqa\n",
    "    keys = [\"Vext_x\", \"Vext_y\", \"Vext_z\", \"alpha\", \"beta\", \"sigma_v\"]\n",
    "\n",
    "    # SN_keys = ['mag_cal', 'alpha_cal', 'beta_cal']\n",
    "    SN_keys = []\n",
    "    calibration_samples = {}\n",
    "    with File(fname, 'r') as f:\n",
    "        for key in keys:\n",
    "            calibration_samples[key] = f[f\"sim_{nsim}/{key}\"][:]\n",
    "\n",
    "        for key in SN_keys:\n",
    "            calibration_samples[key] = f[f\"sim_{nsim}/{key}\"][:]\n",
    "\n",
    "    return calibration_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test running a model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23:16:25: reading the catalogue.\n",
      "23:16:25: reading the interpolated field.\n",
      "23:16:25: calculating the radial velocity.\n",
      "Selected 100/100 galaxies.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mnt/users/rstiskalek/csiborgtools/csiborgtools/flow/flow_model.py:113: UserWarning: The number of radial steps is even. Skipping the first step at 0.0 because Simpson's rule requires an odd number of steps.\n",
      "  warn(f\"The number of radial steps is even. Skipping the first \"\n"
     ]
    }
   ],
   "source": [
    "# fpath_data = \"/mnt/extraspace/rstiskalek/catalogs/PV_compilation_Supranta2019.hdf5\"\n",
    "fpath_data = \"/mnt/extraspace/rstiskalek/catalogs/PV_mock_CB2_17417_small.hdf5\"\n",
    "\n",
    "simname = \"csiborg2_main\"\n",
    "catalogue = \"CB2_small\"\n",
    "# nsim = 0\n",
    "loader = csiborgtools.flow.DataLoader(simname, 19 - 5, catalogue, fpath_data, paths, ksmooth=1)\n",
    "\n",
    "calibration_samples = load_calibration(catalogue, simname, 17417 - 500, 1)\n",
    "flow_model = csiborgtools.flow.get_model(loader, zcmb_max=0.07)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = csiborgtools.flow.Observed2CosmologicalRedshift(calibration_samples, loader.rdist, loader._Omega_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "# zcosmo_mean, zcosmo_std = flow_model.predict_zcosmo_from_calibration(\n",
    "    # calibration_samples[\"mag_cal\"], calibration_samples[\"alpha_cal\"], calibration_samples[\"beta_cal\"])\n",
    "zcosmo_mean = loader.cat[\"zcosmo\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Marginalizing:   0%|          | 0/5000 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Marginalizing: 100%|██████████| 5000/5000 [00:02<00:00, 2262.93it/s]\n"
     ]
    }
   ],
   "source": [
    "# n = 2 is a very good test\n",
    "# CONVERT DEGREES TO RADIANS!!\n",
    "n = 40\n",
    "zcos, post = model.posterior_zcosmo(\n",
    "    loader.cat[\"zobs\"][n], np.deg2rad(loader.cat[\"RA\"][n]), np.deg2rad(loader.cat[\"DEC\"][n]),\n",
    "    loader.los_density[n], loader.los_radial_velocity[n])\n",
    "# zcos, post_bad = model.posterior_zcosmo(\n",
    "#     loader.cat[\"zobs\"][n], loader.cat[\"RA\"][n], loader.cat[\"DEC\"][n],\n",
    "#     loader.los_density[n], loader.los_radial_velocity[n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfLElEQVR4nO3deVxU5eIG8GdmWEZANlEWxVBERRFRQFwyKylod8vletPIrDS9FqVFN6VSL2pmlJqW95prLm1mZmiXxFzAFVRwSU1FlAGFYARkcTi/P/jNXElQloH3zMzz/XzORxzOnHmOp/DxvOc9RyFJkgQiIiIiMnlK0QGIiIiIyDhY7IiIiIjMBIsdERERkZlgsSMiIiIyEyx2RERERGaCxY6IiIjITLDYEREREZkJFjsiIiIiM2ElOoApqKysxNWrV9GyZUsoFArRcYiIiMiCSJKEGzduwMvLC0rl3c/JsdjVwdWrV+Ht7S06BhEREVmwy5cvo127dnddh8WuDlq2bAmg6g/U0dFRcBoiInGKy4vh9ZEXAODqG1dhb2MvOBGhuBjwqjomuHoVsOcxMTdarRbe3t6GPnI3LHZ1oB9+dXR0ZLEjIoumKlcB6qqvHR0dWezkQKX639eOjix2Zqwul4Nx8gQRERGRmWCxIyIiIjITLHZEREREZoLX2BEREZkpnU6HiooK0THoHqytraG6/VrJRmCxIyIiMjOSJEGj0aCgoEB0FKojZ2dneHh4NPp+uSx2REREZkZf6tq0aQM7OzveXF/GJElCSUkJcnNzAQCenp6N2h6LHRERkRnR6XSGUteqVSvRcagOWrRoAQDIzc1FmzZtGjUsy8kTREREZkR/TZ2dnZ3gJFQf+uPV2GsiWeyIiIjMEIdfTYuxjheLHREREZGZYLEjIiIiMhMsdkRERGRS8vLy0KZNG1y8eFF0lDoZPXo0Pvroo2b5LFkWu6VLl8LHxwdqtRphYWE4ePBgret+9913CAkJgbOzM+zt7REUFIS1a9dWW+f555+HQqGotkRGRjb1bhAREVETmDt3Lp555hn4+PiIjlIn7777LubOnYvCwsIm/yzZFbtNmzYhOjoasbGxOHr0KHr27ImIiAjD/V3+ytXVFf/85z+RnJyM48ePIyoqClFRUdixY0e19SIjI5GdnW1YNmzY0By7Q0REREZUUlKC//znP5gwYYLoKHUWEBAAX19frFu3rsk/S3bFbtGiRZg4cSKioqLQrVs3LF++HHZ2dli5cmWN6z/44IMYOnQo/P394evri2nTpiEwMBB79+6ttp6trS08PDwMi4uLS3PsDhEREdWDRqOBQqHAJ598gl69ekGtVqN79+6Gv9e3b98OW1tb9O3bt9r7/vWvf90xOqdQKBAfHy9gL+701FNPYePGjU3+ObIqduXl5Thy5AjCw8MNrymVSoSHhyM5Ofme75ckCYmJiThz5gweeOCBat9LSkpCmzZt0KVLF0yaNAl5eXlGz09E5uHmzZt466238OOPP4qOQmRx0tLSAAArV65EfHw80tLS0L59e4wdOxaVlZXYs2cPgoOD73jf1KlTq43MTZw4Effddx9GjBjRzHtQsz59+uDgwYMoKytr0s+R1ZMnrl+/Dp1OB3d392qvu7u74/Tp07W+r7CwEG3btkVZWRlUKhU+++wzPPLII4bvR0ZGYtiwYejQoQPOnz+Pd955B4899hiSk5NrvLtzWVlZtT94rVZrhL0jIlPxj3/8A//+97+xYsUKaDQa2NjYiI5E1Cj6x1aJUN9Hmh07dgzW1tb44YcfDNfQzZkzByEhIbhy5QouXboELy+vO97XsmVLtGzZEgAwc+ZM7Ny5E0lJSWjXrp1R9qOxvLy8UF5eDo1Gg/vuu6/JPkdWxa6hWrZsibS0NBQVFSExMRHR0dHo2LEjHnzwQQBVs1H0evTogcDAQPj6+iIpKQmDBw++Y3txcXF4//33mys+EcnI+vXr8e9//xsA8Oeff2LHjh146qmnBKciapySkhI4ODgI+eyioiLY29vXef20tDQMGzas2sQIR0dHw9c3b96EWq2u9f2zZs3C2rVrkZSUJKvJFfrHhjV1wZbVUKybmxtUKhVycnKqvZ6TkwMPD49a36dUKtGpUycEBQXhjTfewIgRIxAXF1fr+h07doSbmxvOnTtX4/djYmJQWFhoWC5fvtywHSIik3LmzBm8/PLLAIC2bdsCAL766iuRkYgsTlpaGoKCgqq9lpycDDc3N7Rt2xZubm74888/a3xvbGws1qxZc0epW7lyJQIDA9GzZ0+8+eabAID58+cjICAAPXr0wPr16wFUldDIyEj06NEDPXr0wI4dO3Dx4kX07NkTY8eOhZ+fHyZNmoQtW7YgLCwMAQEBOHv2rOFzatqmXn5+PgCgdevWjf0juitZnbGzsbFBcHAwEhMTMWTIEABAZWUlEhMTMWXKlDpvp7Ky8q5j2FlZWcjLy4Onp2eN37e1tYWtrW29shORabt58yZGjhyJ4uJiPPTQQ4iLi0Pfvn3xww8/4MaNG4YhHiJTZGdnh6KiImGfXVc3b97E2bNnodPpDK9VVlYiPj4e48ePh1KpRK9evWqcXRobG4vVq1ffUepOnDiBjz/+GHv27IGzszPy8/Nx6NAhbN68GYcPH0ZJSQlCQ0Px0EMPITk5Ga1atUJCQgIkScKNGzeQn5+PU6dOYfPmzejUqRMCAgLg4OCAAwcO4PPPP8eSJUvwySef1LpN/bBxeno62rVrBzc3t4b/YdaBrM7YAUB0dDRWrFiB1atX49SpU5g0aRKKi4sRFRUFABg3bhxiYmIM68fFxeGXX37BH3/8gVOnTuGjjz7C2rVr8fe//x1AVfuePn06UlJScPHiRSQmJuKZZ55Bp06dEBERIWQfiUh+XnvtNRw/fhxt2rTB+vXr0adPH/j5+eHmzZv44YcfRMcjahSFQgF7e3shS32urztx4gQUCgXWrVuH5ORknDp1CqNGjUJBQQHeffddAEBERAQyMjKqnbWbM2cOli1bho0bN0KtVkOj0UCj0aCsrAy7du3CqFGj4OzsDKDqNmn79u3D8OHDoVar4erqisGDB+PQoUPo0aMHfvvtN8yYMQMpKSmGIeAuXbqgS5cuUKlU8Pf3N0zy7NGjh+EmybVtU2/Pnj149NFHG3MY60R2xW7UqFFYuHAhZs2ahaCgIKSlpSEhIcEwoSIzMxPZ2dmG9YuLizF58mR0794dAwYMwLfffot169bhxRdfBACoVCocP34cTz/9NDp37owJEyYgODgYe/bs4Vk5IgIAbNiwAV988QUUCgXWr18PT09PKBQK/O1vfwPA4Vii5pKWloauXbvinXfewfDhwxESEgKdTofdu3cbilmPHj3Qu3dvbN68GUDVxJAPP/wQ165dQ79+/eDp6WlYjh8/Xq/P79y5M9LS0tC9e3dER0djyZIlAFCtLyiVSsPvlUpltbOLtSktLcWWLVswceLEeuVpEInuqbCwUAIgFRYWio5CREZ25swZycHBQQIgzZw5847vAZBUKpWUm5srKKG8FJUVSXgPEt6DVFRWJDoOSZIkFRVJElC1FBVJN2/elE6ePCndvHlTdLJ6mzx5sjRmzJh7rrdt2zbJ399f0ul091z3xIkTUkBAgFRQUCBJkiTl5eVJhw4dkoKDg6XS0lIpPz9f8vX1la5evSpduXLF8Oe2ceNGacKECdKFCxek4OBgw/aGDx8u7dq1S5IkSUpOTpaeeOIJSZKkWrcpSZL02WefSY888shdc97tuNWnh8jqGjsiouZUWlqKkSNHoqioCIMGDUJsbGy173fu3BkhISE4fPgwvv76a0yePFlQUiLLkJaWVqdZ6E888QTOnj2LK1euwNvb+67rBgQEYNq0aRgwYACsrKzw6KOPYsGCBXj22WcRHBwMhUKB999/H56entixYwfefPNNqFQqtGjRAv/5z3/qnD0kJKTGbQKAtbU1Fi9eXOdtNYZCkiSpWT7JhGm1Wjg5OaGwsLDalGsiMm2TJk3C8uXL0bp1a6SlpdV4b6yPP/4Y0dHR6N+/P/bt2ycgpbwUlxfDIa7qthlFMUWwt6n7bSyoiRQXA/pbmRQVoVSlwoULF9ChQ4e73hZEbiRJgpOTEzZu3IjHH39cdJxmV1paWutxq08Pkd01dkREzWHTpk1Yvny54ULtmkodUHXdr0KhwP79+3HhwoVmTklkORQKBbRarUWWOmNisSMii3P27FnDRczvvPPOXWeqeXl54aGHHgKAZnnOIxFRY7DYEZFF0V9Xd+PGDQwcOBDvvffePd8zduxYAJwdS0Tyx2JHRBbljTfeQFpaGtzc3LBhwwZYWd17DtmwYcNgY2OD9PR0nDhxohlSEhE1DIsdEVmMzZs347PPPgMArF271vDYsHtxdnbGE088AQB3PCaIiEhOWOyIyCKcO3fOcOPymJgYREZG1uv9+psVb9iwAZWVlUbPR0RkDCx2RGT2br+u7v7778cHH3xQ72088cQTaNmyJTIzM7F///4mSElE1HgsdkRk9t58802kpqaiVatWdb6u7q9atGiB4cOHA+AkCiKSLxY7IjJr33zzDZYuXQqg6rq6du3aNXhb+uHYzZs3o6Kiwij5iIiMicWOiMzW+fPnMWHCBADAW2+9hccee6xR23vooYfg7u6OvLw87Ny50xgRiYiMisWOiMxSWVkZRo0aBa1WiwEDBmD27NmN3qaVlRVGjRoFgMOxRCRPLHZEZJamT5+OI0eOwNXVFRs2bIC1tbVRtqsfjt2yZQuKi4uNsk0iImNhsSMis/Pdd99h8eLFAIA1a9bA29vbaNvu06cPfH19UVJSgq1btxptu0RU5V//+hcUCsUdS3x8vOhoJoHFjojMyh9//IEXXngBQNVZO/2NhY1FoVAYztpxOJbI+KZOnYrs7GzDMnHiRNx3330YMWKE6Ggmof5z/omIZKq8vByjR49GYWEh+vXrh7lz5zbJ54wZMwazZ89GQkICrl+/Djc3tyb5HCJjkSQJJRUlQj7bztoOCoWizuu3bNkSLVu2BADMnDkTO3fuRFJSUr1mtF+8eBEjRozA4cOH653X1LHYEZHZmDFjBg4dOgQXFxds3LjRaNfV/ZW/vz969eqF1NRUfPPNN3jllVea5HOIjKWkogQOcQ5CPrsopgj2Nvb1ft+sWbOwdu1aJCUlwcfHx/jBzBSHYonILKSkpOCTTz4BAKxevRrt27dv0s/jcCxR04mNjcWaNWvuKHUrV65EYGAgevbsiTfffBMA8OSTTyI4OBgBAQF8ljN4xo6IzERiYiIAYNiwYXjqqaea/PNGjx6NGTNmYM+ePcjMzGzyIknUGHbWdiiKKRL22fURGxuL1atX31HqTpw4gY8//hh79uyBs7Mz8vPzAVRNkHJ1dUVxcTFCQ0Mt/lo8FjsiMgupqakAgL59+zbL57Vr1w6DBg1CUlISNm7ciBkzZjTL5xI1hEKhaNBwaHObM2cOli1bhq1bt0KtVkOj0QAAXFxcsGvXLowaNQrOzs4AAFdXVwDAxx9/bJihnpmZiczMzCa7DMMUcCiWiMyCvtj16tWr2T5TPxzL4R+ixpMkCR9++CGuXbuGfv36wdPT07AcP368xvfs2rUL+/btw4EDB3Ds2DF07doVZWVlzZxcXljsiMjkFRYW4o8//gDQvMVu+PDhsLa2xvHjx5Gent5sn0tkjhQKBQoLCyFJ0h1LaGgoHn74YWzatAmFhYUAgPz8fGi1WrRq1QpqtRppaWk4duyY4L0Qj8WOiExeWloaAMDb2xutWrVqts91dXU1PH92w4YNzfa5RJYoICAA06ZNw4ABAxAUFIR58+YhMjISN27cQLdu3TB37lwEBweLjikcr7EjIpMnYhhW729/+xu2bt2Kr776CnPmzKnX/bqIqH5efPFFvPjii9VeS0hIqHFdS7yHHcAzdkRkBkQWu6eeegoODg64ePEiUlJSmv3ziYhux2JHRCZPZLGzs7PD0KFDAXASBRGJx2JHRCattLQUJ0+eBCCm2AH/mx27efNmVFRUCMlARASw2BGRiUtPT4dOp4Orqyu8vb2FZBg8eDBat26Na9euGW6UTEQkAosdEZm024dhRU1csLa2xsiRIwHwEWNEJBaLHRGZNJHX191u7NixAIDvv/8eJSUlQrMQkeVisSMikyaXYte3b1/4+PigqKgIP/74o9AsRGS5WOyIyGTpdDrDo4ZEFzuFQmGYRMHhWCIShcWOiEzW77//jpKSEtjZ2aFz586i4xiK3c8//4z8/HzBaYjIErHYEZHJ0g/DBgYGQqVSCU4DdO/eHYGBgaioqMC3334rOg4RWSAWOyIyWXK5vu52+kkUHI4lIhFY7IjIZMmx2I0ePRoAsHv3bmRlZQlOQ2R6/vWvf0GhUNyxxMfHi45mEljsiMgkSZIky2LXvn17DBw4EJIkYePGjaLjEJmcqVOnIjs727BMnDgR9913H0aMGCE6mklgsSMik3T58mXk5+dDpVIhICBAdJxqODuWZEeSgOJiMYsk1Stqy5Yt4eHhAQ8PDyxduhQ7d+5EUlIS2rVrh6ysLAwbNgy+vr4ICQnBs88+i5ycHABVM9MnT55s2E52djZUKhXee+89AICVlRWCgoIQFBSE0NBQpKWlGetPV1asRAcgImoI/dm6bt26Qa1WC05T3YgRIzB16lSkpqbi1KlT8Pf3Fx2JLF1JCeDgIOazi4oAe/t6v23WrFlYu3YtkpKS4OPjA0mS8Mwzz2Dy5Mn47rvvAAB79uzBtWvX4O7uDldXV6SkpECn00GlUuGbb75B9+7dDdtzdnY2lLlvv/0WH3zwgWE75oRn7IjIJMlxGFbPzc0NkZGRAIANGzYITkNkemJjY7FmzRpDqQOAxMREODg4YMKECYb1Bg4caDhjr1AoMHDgQOzevRtA1VNghg0bVuP2tVotnJ2dm3QfROEZOyIySXIudkDVcOy2bduwfv16vP/++8KeY0sEALCzqzpzJuqz6yE2NharV6+uVuoA4OTJk+jdu/dd3zty5EisXbsWXbt2hY2NDdzc3HD9+nUAQEFBAYKCglBSUoK8vDzs37+/3rtiCljsiMgkyb3YPf3007Czs8Mff/yBgwcPIiwsTHQksmQKRYOGQ5vbnDlzsGzZMmzduhVqtRoajQYA4OLiUqf39+/fH1OnTsXGjRsxYsQIlJaWGr53+1DsN998g1dffRX//e9/jb4PonEolohMTl5eHi5fvgwACAoKEhumFvb29hgyZAgATqIgqgtJkvDhhx/i2rVr6NevHzw9PQ3L8ePH4e/vb/gHXW0UCgUeeOABzJs3D0OHDq11vSeffNJsz9ix2BGRydH/cPf19YWTk5PgNLXTz47dtGkTbt26JTgNkbwpFAoUFhZCkqQ7ltDQUISHh0Or1WLVqlWG9+zduxfp6enVtvPqq69i/vz5aNWqVa2ftX//fnTs2LGpdkUoWRa7pUuXwsfHB2q1GmFhYTh48GCt63733XcICQmBs7Mz7O3tERQUhLVr11ZbR5IkzJo1C56enmjRogXCw8Nx9uzZpt4NImoich+G1Xv00UfRqlUr5OTkYNeuXaLjEJk0hUKBLVu2YMuWLfD19UX37t2xePFitG7dutp6fn5+iIqKuuP9+mvsevbsienTp+OLL75orujNSnbX2G3atAnR0dFYvnw5wsLCEB8fj4iICJw5cwZt2rS5Y31XV1f885//NFwouW3bNkRFRaFNmzaIiIgAACxYsACffvopVq9ejQ4dOmDmzJmIiIjAyZMnZXebBCK6N1MpdtbW1hg5ciSWLVuG9evX45FHHhEdiciktW/fHlu2bKnxe/pJErebMmWK4WtLOWsuuzN2ixYtwsSJExEVFYVu3bph+fLlsLOzw8qVK2tc/8EHH8TQoUPh7+8PX19fTJs2DYGBgdi7dy+AqrN18fHxePfdd/HMM88gMDAQa9aswdWrV2v9j4OI5M1Uih0AjBkzBgDwww8/oLKyUnAaIjJ3sip25eXlOHLkCMLDww2vKZVKhIeHIzk5+Z7vlyQJiYmJOHPmDB544AEAwIULF6DRaKpt08nJCWFhYbVus6ysDFqtttpCRPJQXFyMM2fOADCNYte3b1+o1WoUFBTg3LlzouMQkZmTVbG7fv06dDod3N3dq73u7u5umPJck8LCQjg4OMDGxgZPPPEEFi9ebBjy0L+vPtuMi4uDk5OTYfH29m7MbhGRER0/fhySJBkeOSR31tbWhpm7hw4dEhuGiMyerIpdQ7Vs2RJpaWk4dOgQ5s6di+joaCQlJTV4ezExMSgsLDQs+tsqEJF4pjQMq9enTx8AuOtEMCIiY5DV5Ak3NzeoVCrDA331cnJy7vovc6VSiU6dOgGouqfVqVOnEBcXhwcffNDwvpycHHh6elbbZm33v7K1tYWtrW0j94aImoIpFrvQ0FAAPGNHRE1PVmfsbGxsEBwcjMTERMNrlZWVSExMRL9+/eq8ncrKSpSVlQEAOnToAA8Pj2rb1Gq1OHDgQL22SUTyYMrFLjU1FRUVFYLTkKWQJEl0BKoHYx0vWZ2xA4Do6GiMHz8eISEh6NOnD+Lj41FcXGy4J824cePQtm1bxMXFAai6Hi4kJAS+vr4oKyvD9u3bsXbtWixbtgxA1X1vXnvtNcyZMwd+fn6G2514eXkZ7gpPRKahoqICJ06cAGBaxc7Pzw9OTk4oLCxEenq6SWUn02NtbQ0AKCkpQYsWLQSnoboqKSkB8L/j11CyK3ajRo3CtWvXMGvWLGg0GgQFBSEhIcEw+SEzMxNK5f9ONBYXF2Py5MnIyspCixYt0LVrV6xbtw6jRo0yrDNjxgwUFxfjpZdeQkFBAe6//34kJCTwHnZEJubUqVMoLy+Ho6MjOnToIDpOnSmVSoSEhCAxMRGHDh1isaMmpVKp4OzsjNzcXACAnZ0dFAqF4FRUG0mSUFJSgtzcXDg7O0OlUjVqewqJ52rvSavVGv617ejoKDoOkcVavXo1nn/+eTzwwAPYvXu36Dj1EhMTg3nz5uHFF1/EihUrRMdpsOLyYjjEOQAAimKKYG8j/wfLm73iYsCh6pigqAiwt4ckSdBoNCgoKBAajerO2dkZHh4eNZbw+vQQ2Z2xIyKqjSleX6fHmbHUnBQKBTw9PdGmTRte12kCrK2tG32mTo/FjohMhikXO/0EioyMDJSUlMDOzk5wIrIEKpXKaIWBTIOsZsUSEdWmsrISaWlpAEyz2LVt2xYeHh7Q6XSGgkpEZGwsdkRkEi5cuACtVgtbW1v4+/uLjlNvCoWCw7FE1ORY7IjIJOjPcgUEBDT6dgCi8EbFRNTUWOyIyCSY8vV1eix2RNTUWOyIyCSYU7E7d+4c8vPzBachInPEYkdEJsEcip2rqyt8fX0BAIcPHxachojMEYsdEcmeRqOBRqOBQqFAYGCg6DiNwuFYImpKLHZEJHv6s3VdunSBvb1pP+mAM2OJqCmx2BGR7JnDMKwez9gRUVNisSMi2TOnYterVy8olUpkZ2fjypUrouMQkZlhsSMi2TOnYmdvb4+AgAAAHI4lIuNjsSMiWSssLMT58+cBmEexAzgcS0RNh8WOiGTt2LFjAABvb2+0atVKcBrjYLEjoqbCYkdEsmZOw7B6+pmxhw4dQmVlpeA0RGROWOyISNbMsdgFBARArVajsLAQ586dEx2HiMwIix0RyZo5Fjtra2sEBQUB4HAsERkXix0RyVZZWRlOnjwJwLyKHcAbFRNR02CxIyLZSk9Px61bt+Dq6gpvb2/RcYyKEyiIqCmw2BGRbN0+DKtQKASnMS59sUtNTUVFRYXgNERkLljsiEi2zPH6Oj0/Pz84OTmhtLQU6enpouMQkZlgsSMi2TLnYqdUKhESEgKAw7FEZDwsdkQkSzqdznBzYnMsdgCvsyMi42OxIyJZOnv2LEpKSmBnZ4fOnTuLjtMkODOWiIyNxY6IZEk/DBsYGAiVSiU4TdPQn7HLyMhASUmJ4DREZA5Y7IhIlsz5+jq9tm3bwsPDAzqdzrC/RESNwWJHRLJkCcVOoVBwOJaIjIrFjohkR5Ikiyh2ACdQEJFxsdgRkexkZWUhLy8PKpUKAQEBouM0KRY7IjImFjsikh392bpu3bpBrVYLTtO09MXu3LlzyM/PF5yGiEwdix0RyY6lDMMCgKurK3x9fQEAhw8fFpyGiEwdix0RyY4lFTuAw7FEZDwsdkQkO5ZW7DgzloiMhcWOiGQlLy8PmZmZAICgoCCxYZoJz9gRkbGw2BGRrKSlpQEAfH194eTkJDZMM+nVqxeUSiWys7Nx5coV0XGIyISx2BGRrFjaMCwA2NvbG27rwuFYImoMFjsikhVLLHYAh2OJyDhY7IhIVljsWOyIqOFY7IhINkpKSnDmzBkAllfs9DNjDx06hMrKSsFpiMhUsdgRkWwcP34clZWV8PDwgIeHh+g4zSogIABqtRqFhYU4d+6c6DhEZKJY7IhINix1GBYArK2tDbd34XAsETUUix0RyYYlFzuANyomosZjsSMi2bD0YscJFETUWLIsdkuXLoWPjw/UajXCwsLu+q/XFStWYODAgXBxcYGLiwvCw8PvWP/555+HQqGotkRGRjb1bhBRPVRUVODEiRMAWOxSU1NRUVEhOA0RmSLZFbtNmzYhOjoasbGxOHr0KHr27ImIiAjk5ubWuH5SUhLGjBmDXbt2ITk5Gd7e3nj00UfvuHt7ZGQksrOzDcuGDRuaY3eIqI5Onz6NsrIyODo6okOHDqLjCOHn5wcnJyeUlpYiPT1ddBwiMkGyK3aLFi3CxIkTERUVhW7dumH58uWws7PDypUra1x//fr1mDx5MoKCgtC1a1f8+9//RmVlJRITE6utZ2tra5hp5+HhARcXl+bYHSKqI/0wbFBQEJRK2f1oahZKpRIhISEAOBxLRA0jq5+e5eXlOHLkCMLDww2vKZVKhIeHIzk5uU7bKCkpQUVFBVxdXau9npSUhDZt2qBLly6YNGkS8vLyat1GWVkZtFpttYWImpalX1+nx+vsiKgxZFXsrl+/Dp1OB3d392qvu7u7Q6PR1Gkbb731Fry8vKqVw8jISKxZswaJiYmYP38+du/ejcceeww6na7GbcTFxcHJycmweHt7N3yniKhOWOyqcGYsETWGlegAxjRv3jxs3LgRSUlJUKvVhtdHjx5t+LpHjx4IDAyEr68vkpKSMHjw4Du2ExMTg+joaMPvtVotyx1RE5IkCWlpaQBY7PRn7DIyMlBSUgI7OzvBiYjIlMjqjJ2bmxtUKhVycnKqvZ6Tk3PPu9AvXLgQ8+bNw86dOxEYGHjXdTt27Ag3N7da7+5ua2sLR0fHagsRNZ0LFy6gsLAQtra28Pf3Fx1HqLZt28LDwwM6nc5wFpOIqK5kVexsbGwQHBxcbeKDfiJEv379an3fggULMHv2bCQkJBguPL6brKws5OXlwdPT0yi5iahx9AUmICAA1tbWgtOIpVAoOBxLRA0mq2IHANHR0VixYgVWr16NU6dOYdKkSSguLkZUVBQAYNy4cYiJiTGsP3/+fMycORMrV66Ej48PNBoNNBoNioqKAABFRUWYPn06UlJScPHiRSQmJuKZZ55Bp06dEBERIWQfiag6Xl9XHSdQEFFDye4au1GjRuHatWuYNWsWNBoNgoKCkJCQYJhQkZmZWe1WCMuWLUN5eTlGjBhRbTuxsbF47733oFKpcPz4caxevRoFBQXw8vLCo48+itmzZ8PW1rZZ942IasZiVx2LHRE1lEKSJEl0CLnTarVwcnJCYWEhr7cjagJeXl7Izs7G/v3773rZhaXIz89Hq1atAAB5eXl33L5JpOLyYjjEOQAAimKKYG9jLzgRobgYcKg6JigqAux5TMxNfXqI7IZiiciy5OTkIDs7GwqF4p4TnyyFq6srfH19AQCHDx8WnIaITAmLHREJpR+G7dKlC+x5psGAw7FE1BAsdkQkFK+vqxlnxhJRQ7DYEZFQLHY14xk7ImoIFjsiEorFrma9evWCUqlEdnY2rly5IjoOEZkIFjsiEkar1RqeAMNiV529vT0CAgIAcDiWiOqOxY6IhDl27BgAwNvb23B7D/ofDscSUX2x2BGRMGlpaQB4tq42LHZEVF8sdkQkzIkTJwCA96+rhX5m7KFDh1BZWSk4DRGZAhY7IhJGX+z015JRdQEBAVCr1SgsLDRci0hEdDcsdkQkhCRJSE9PBwD06NFDcBp5sra2RlBQEAAOxxJR3bDYEZEQly5dQlFREWxsbODn5yc6jmzxRsVEVB8sdkQkhH4YtmvXrrC2thacRr44gYKI6oPFjoiE4DBs3eiLXWpqKioqKgSnISK5Y7EjIiE4caJu/Pz84OTkhNLSUkMZJiKqDYsdEQnBM3Z1o1QqERISAoDDsUR0byx2RNTsKioqcPr0aQAsdnXB6+yIqK5Y7Iio2Z05cwYVFRVwdHSEt7e36Diyx5mxRFRXLHZE1Oz0w7ABAQFQKBSC08if/oxdRkYGSkpKBKchIjljsSOiZseJE/XTtm1beHh4QKfTITU1VXQcIpIxFjsianb6Ysfr6+pGoVBwOJaI6oTFjoiaHWfE1h8nUBBRXbDYEVGzunHjBi5cuACAQ7H1wWJHRHXBYkdEzerkyZMAAE9PT7Rq1UpwGtOhL3bnzp1Dfn6+4DREJFcsdkTUrDhxomFcXV3h6+sLADh8+LDgNEQkVyx2RNSsOHGi4TgcS0T3wmJHRM2KEycajjNjieheWOyIqFlxKLbheMaOiO6FxY6Imk1ubi6uXbsGhUKBbt26iY5jcnr37g2VSoXs7GxkZWWJjkNEMsRiR0TNRn+2ztfXF3Z2doLTmB47OzsEBgYCAFJSUgSnISI5YrEjombDiRON17dvXwAsdkRUMxY7Imo2nDjReGFhYQBY7IioZix2RNRsOHGi8fRn7I4cOYKKigrBaYhIbljsiKhZVFZWIiMjAwDP2DWGn58fXFxcUFpaiuPHj4uOQ0Qyw2JHRM3i4sWLKC4uhq2tLTp16iQ6jslSKpUcjiWiWrHYEVGz0A/D+vv7w8rKSnAa08ZiR0S1YbEjombBiRPGw5mxRFQbFjsiahacOGE8+keLnTt3Dnl5eYLTEJGcsNgRUbPgPeyMx9XVFV26dAEAHDhwQHAaIpITFjsianJlZWX4/fffAfCMnbFwOJaIasJiR0RN7syZM7h16xacnJzQrl070XHMAidQEFFNWOyIqMndPnFCoVAITmMe9GfsDhw4gMrKSsFpiEguWOyIqMlx4oTx9ejRAy1atIBWq8WZM2dExyEimWCxI6Imx4kTxmdlZYXQ0FAAHI4lov+RZbFbunQpfHx8oFarERYWhoMHD9a67ooVKzBw4EC4uLjAxcUF4eHhd6wvSRJmzZoFT09PtGjRAuHh4Th79mxT7wYR/T/9UCzP2BkXJ1AQ0V/Jrtht2rQJ0dHRiI2NxdGjR9GzZ09EREQgNze3xvWTkpIwZswY7Nq1C8nJyfD29sajjz6KK1euGNZZsGABPv30UyxfvhwHDhyAvb09IiIiUFpa2ly7RWSxtFotLl26BIBn7IyNEyiI6K8UkiRJokPcLiwsDKGhoViyZAmAqgeHe3t7Y+rUqXj77bfv+X6dTgcXFxcsWbIE48aNgyRJ8PLywhtvvIE333wTAFBYWAh3d3esWrUKo0ePvuc2tVotnJycUFhYCEdHx8btIJGFSU5ORv/+/dG2bVtkZWWJjmNWrl69irZt20KpVKKwsBAODg5N/pnF5cVwiKv6nKKYItjb2Df5Z9I9FBcD+mNfVATY85iYm/r0EFmdsSsvL8eRI0cQHh5ueE2pVCI8PBzJycl12kZJSQkqKirg6uoKALhw4QI0Gk21bTo5OSEsLKzO2ySihuPEiabj5eUFb29vVFZW4vDhw6LjEJEMyKrYXb9+HTqdDu7u7tVed3d3h0ajqdM23nrrLXh5eRmKnP599dlmWVkZtFpttYWIGoYTJ5oWr7MjotvJqtg11rx587Bx40Z8//33UKvVDd5OXFwcnJycDIu3t7cRUxJZFk6caFq8zo6IbierYufm5gaVSoWcnJxqr+fk5MDDw+Ou7124cCHmzZuHnTt3IjAw0PC6/n312WZMTAwKCwsNy+XLlxuyO0QWT5IknrFrYrefsZPZJdNEJICsip2NjQ2Cg4ORmJhoeK2yshKJiYno169fre9bsGABZs+ejYSEBISEhFT7XocOHeDh4VFtm1qtFgcOHKh1m7a2tnB0dKy2EFH95eTkIC8vD0qlEv7+/qLjmKXevXvDysoKOTk5yMzMFB2HiASTVbEDgOjoaKxYsQKrV6/GqVOnMGnSJBQXFyMqKgoAMG7cOMTExBjWnz9/PmbOnImVK1fCx8cHGo0GGo0GRUVFAACFQoHXXnsNc+bMwdatW3HixAmMGzcOXl5eGDJkiIhdJLIY+rN1nTp1QosWLQSnMU8tWrRAUFAQAA7HEhFgJTrAX40aNQrXrl3DrFmzoNFoEBQUhISEBMPkh8zMTCiV/+ujy5YtQ3l5OUaMGFFtO7GxsXjvvfcAADNmzEBxcTFeeuklFBQU4P7770dCQkKjrsMjonvjMGzz6Nu3Lw4fPoyUlBSMGjVKdBwiEqhRxa6iogIajQYlJSVo3bq14RYjjTVlyhRMmTKlxu8lJSVV+/3FixfvuT2FQoEPPvgAH3zwgRHSEVFdceJE8wgLC8OSJUt4xo6I6j8Ue+PGDSxbtgyDBg2Co6MjfHx84O/vj9atW+O+++7DxIkTcejQoabISkQmhmfsmod+AkVqairKysoEpyEikepV7BYtWgQfHx98+eWXCA8Px5YtW5CWlobff/8dycnJiI2Nxa1bt/Doo48iMjKSz2MlsmA6nQ4ZGRkAWOyamq+vL1q1aoWysjIcO3ZMdBwiEqheQ7GHDh3Cb7/9hu7du9f4/T59+uCFF17A8uXL8eWXX2LPnj3w8/MzSlAiMi0XLlzAzZs3oVar4evrKzqOWVMoFOjbty9++uknpKSkoE+fPqIjEZEg9Tpjt2HDBkOp+/TTT3H16tUa17O1tcUrr7yCF154ofEJicgk6Ydhu3XrBpVKJTiN+eMTKIgIaMTtTl577TUMHDjwjpv36p/3SkSWjRMnmhefQEFEQCPvYxceHo5BgwZVK3d//vknhwGIiBMnmlmfPn2gUChw4cIF5Obmio5DRII0uNgpFArMnj0bY8eOvaPc8bE2RMRi17ycnJwMT/c4cOCA4DREJEqjnzwxe/Zs/P3vf69W7hQKRaODEZHpKi0tNcyK51Bs8+F1dkTU4GJ3+1m5Dz74AM899xwGDRqES5cuGSUYEZmu06dPQ6fTwcXFBV5eXqLjWAxeZ0dEDX7yxNy5c2Fvb2/4/fvvvw8AeOqppxqfiohM2u0TJ3gGv/noz9gdPHgQOp2Os5GJLFCDz9jFxMRUK3ZAVbmbOnUqWrZs2ehgRGS6eH2dGN27d4e9vT2Kiopw6tQp0XGISIBGX2P3V++++y4KCgqMvVkiMiEsdmKoVCrDXQk4HEtkmepV7DIzM+u18StXrtRrfSIyD7yHnTicQEFk2epV7EJDQ/Hyyy/j0KFDta5TWFiIFStWICAgAN9++22jAxKRaSkoKDDMkGexa36cQEFk2eo1eeLkyZOYO3cuHnnkEajVagQHB8PLywtqtRp//vknTp48iYyMDPTu3RsLFizA448/3lS5iUimMjIyAADt2rWDs7Oz2DAWSF/sTp48Ca1WC0dHR8GJiKg51euMXatWrbBo0SJkZ2djyZIl8PPzw/Xr1w33qxo7diyOHDmC5ORkljoiC8Xr68Ty8PCAj48PJEm66+gKEZmnBt3upEWLFhgxYgRGjBhhuI6ubdu2Rg1GRKaJxU68vn374uLFi0hJScHgwYNFxyGiZtTgWbH79u1Dhw4d0L59e7Rv3x7u7u546623oNVqjZmPiEwMJ06IxwkURJarwcXu5Zdfhr+/Pw4dOoQzZ87gww8/xH//+1/07t2bs2GJLJQkSTxjJwO3T6Dgs7uJLItCauD/9S1atMCxY8fQuXNnw2uSJGHkyJEAgK+//to4CWVAq9XCyckJhYWFvBCZ6C6uXLmCdu3aQaVSoaioCGq1WnQki1RWVgZHR0eUl5fj/Pnz6Nixo9G2XVxeDIc4BwBAUUwR7G3s7/EOanLFxYBD1TFBURFgz2NiburTQxp8xs7f3x+5ubnVXlMoFPjggw+QkJDQ0M0SkQnTD8P6+fmx1Alka2uLXr16AeBwLJGlaXCxe/755zF16lTD/ar0eFaLyHJxGFY+eJ0dkWVq0KxYAHjttdcAVP3LfNiwYQgKCoJOp8O6deuwYMECY+UjIhPCiRPywRsVE1mmBhe77OxspKWl4dixY0hLS8OqVatw9uxZKBQKLFiwAD///DMCAwMRGBiIyMhIY2YmIpniGTv50J+xS0tLQ2lpKYfGiSxEgydP1KS0tBQnTpyoVvjS09NRUFBgrI8QgpMniO5Np9PBwcEBpaWl+P333+Hn5yc6kkWTJAkeHh7Izc3F/v370a9fP6Nsl5MnZIiTJ8xefXpIg8/Y1UStViM0NBShoaHG3CwRmYDz58+jtLQULVq0MOosTGoYhUKBvn37YuvWrUhJSTFasSMieWvw5Akiotvph2G7desGlUolOA0BnEBBZIlY7IjIKPQTJ3h9nXxwAgWR5WGxIyKj4MQJ+QkNDYVCoUBmZiays7NFxyGiZsBiR0RGoS92vNWJfLRs2dJwPA4cOCA4DRE1BxY7Imq0mzdv4ty5cwB4xk5ueJ0dkWVhsSOiRjt16hQqKyvh6uoKDw8P0XHoNix2RJaFxY6IGu32iRMKhUJwGrqdfgLFoUOHcOvWLcFpiKipsdgRUaNx4oR8+fv7w9HRESUlJcjIyBAdh4iaGIsdETUaJ07Il1KpRJ8+fQBwOJbIErDYEVGj8R528sbr7IgsB4sdETXKn3/+iStXrgAAunfvLjgN1YQ3KiayHCx2RNQo+mHY9u3bw8nJSXAaqom+2J0+fRoFBQViwxBRk2KxI6JG4TCs/LVu3Rq+vr4AgIMHDwpOQ0RNicWOiBqFEydMA6+zI7IMLHZE1Cg8Y2caWOyILAOLHRE1mCRJPGNnIvTX2R04cACSJAlOQ0RNhcWOiBosKysLhYWFUKlU6Nq1q+g4dBc9e/aEra0t8vPzDc/1JSLzw2JHRA2mH4bt0qULbG1tBaehu7GxsUFwcDAADscSmTPZFbulS5fCx8cHarUaYWFhd53BlZGRgeHDh8PHxwcKhQLx8fF3rPPee+9BoVBUW3hmgcg4OAxrWnidHZH5k1Wx27RpE6KjoxEbG4ujR4+iZ8+eiIiIQG5ubo3rl5SUoGPHjpg3bx48PDxq3W737t2RnZ1tWPbu3dtUu0BkUThxwrSw2BGZP1kVu0WLFmHixImIiopCt27dsHz5ctjZ2WHlypU1rh8aGooPP/wQo0ePvuswkJWVFTw8PAyLm5tbU+0CkUXhGTvTop9AcezYMZSUlAhOQ0RNQTbFrry8HEeOHEF4eLjhNaVSifDwcCQnJzdq22fPnoWXlxc6duyIsWPHIjMzs7FxiSzerVu3cOrUKQA8Y2cqvL294enpCZ1Oh6NHj4qOQ0RNQDbF7vr169DpdHB3d6/2uru7OzQaTYO3GxYWhlWrViEhIQHLli3DhQsXMHDgQNy4caPW95SVlUGr1VZbiKi6c+fOoaysDHZ2dujQoYPoOFQHCoWCw7FEZk42xa6pPPbYY3j22WcRGBiIiIgIbN++HQUFBdi8eXOt74mLi4OTk5Nh8fb2bsbERKZBPwzbvXt3KJVm/6PEbLDYEZk32fw0dnNzg0qlQk5OTrXXc3Jy7joxor6cnZ3RuXPnu97HKSYmBoWFhYbl8uXLRvt8InPBiROmSX+dHYsdkXmSTbHT32MpMTHR8FplZSUSExPRr18/o31OUVERzp8/D09Pz1rXsbW1haOjY7WFiKrjxAnTFBISAqVSiStXriArK0t0HCIyMtkUOwCIjo7GihUrsHr1apw6dQqTJk1CcXExoqKiAADjxo1DTEyMYf3y8nKkpaUhLS0N5eXluHLlCtLS0qqdjXvzzTexe/duXLx4Efv378fQoUOhUqkwZsyYZt8/InOiL3Y8Y2da7O3tERgYCKDq8WJEZF6sRAe43ahRo3Dt2jXMmjULGo0GQUFBSEhIMEyoyMzMrHYtz9WrV9GrVy/D7xcuXIiFCxdi0KBBSEpKAlD1yKMxY8YgLy8PrVu3xv3334+UlBS0bt26WfeNyJyUlJTg/PnzAFjsTFHfvn2RlpaGlJQUDB8+XHQcIjIiWRU7AJgyZQqmTJlS4/f0ZU3Px8fnng+z3rhxo7GiEdH/O3nyJCRJgpubG9q0aSM6DtVT3759sXz5cl5nR2SGZDUUS0Sm4fZhWIVCITgN1Zd+AsXhw4dRUVEhOA0RGROLHRHVm35GLCdOmKbOnTvD2dkZpaWlhpJOROaBxY6I6o0TJ0ybUqnkbU+IzBSLHRHVG+9hZ/p4o2Ii88RiR0T1kpeXh+zsbABVT50g08RiR2SeWOyIqF70w7A+Pj5o2bKl4DTUUH369AEAnD17Frm5uYLTEJGxsNgRUb3wiRPmwdXVFUFBQQCAhIQEsWGIyGhY7IioXvRDd7179xachBrrySefBAD8+OOPgpMQkbGw2BFRnUmShN27dwMAHnjgAcFpqLGeeuopAMCOHTtQXl4uOA0RGQOLHRHV2R9//IErV67A2toa/fr1Ex2HGikkJAQeHh64ceOGobATkWljsSOiOtP/5R8aGgo7OzvBaaixlEolnnjiCQAcjiUyFyx2RFRn+mI3aNAgwUnIWPTDsdu2bbvns7eJSP5Y7IiozljszE94eDhsbW1x4cIFnDx5UnQcImokFjsiqpNLly7h0qVLUKlU6N+/v+g4ZCT29vZ4+OGHAXA4lsgcsNgRUZ3oz9YFBwfzxsRmRj8cy2JHZPpY7IioTn777TcAHIY1R/r72SUnJ+PatWuC0xBRY7DYEVGd8Po68+Xt7Y2goCBIkoTt27eLjkNEjcBiR0T3dPXqVZw7dw4KhQIDBgwQHYeaAIdjicwDix0R3ZP+bF1QUBCcnZ3FhqEmoR+O5VMoiEwbix0R3ROHYc2f/ikURUVFfAoFkQljsSOie2KxM398CgWReWCxI6K7ysnJwenTpwEAAwcOFJyGmtLt19nxKRREponFjojuSn+bkx49eqBVq1aC01BT0j+F4uLFi8jIyBAdh4gagMWOiO6Kw7CWw97eHoMHDwbA4VgiU8ViR0R3xWJnWfSzY7dt2yY4CRE1BIsdEdUqLy8P6enpAIAHHnhAcBpqDnwKBZFpY7Ejolrt2bMHAODv7482bdoITkPNgU+hIDJtLHZEVCsOw1omPoWCyHSx2BFRrVjsLJO+2O3YsQNlZWWC0xBRfbDYEVGNCgoKkJaWBoDX11ma4OBgPoWCyESx2BFRjfbu3QtJktCpUyd4eXmJjkPN6PanUHB2LJFpYbEjohpxGNay8SkURKaJxY6IasRiZ9n4FAoi08RiR0R3uHHjBo4ePQqAxc5S8SkURKaJxY6I7rBv3z7odDr4+Pigffv2ouOQILztCZHpYbEjojtwGJaA/z2FIiUlBbm5uYLTEFFdsNgR0R1+++03ACx2lq5du3aGp1D8/PPPouMQUR2w2BFRNSUlJTh06BAAFjvicCyRqWGxI6JqkpOTUVFRgXbt2qFDhw6i45BgfAoFkWlhsSOiam6/vk6hUAhOQ6LxKRREpoXFjoiq4cQJup1SqTRMouBwLJH8sdgRkUFpaSkOHDgAgM+Hpf/hUyiITAeLHREZHDhwAGVlZXB3d0fnzp1FxyGZGDx4MGxtbXHp0iVknORTKIjkjMWOiAx4fR3V5PanUCT8nCA4DRHdDYsdERnw+jqqjX44dvvP2wUnIaK7kV2xW7p0KXx8fKBWqxEWFoaDBw/Wum5GRgaGDx8OHx8fKBQKxMfHN3qbRJaqvLwcycnJAFjs6E76CRT8+Ukkb7Iqdps2bUJ0dDRiY2Nx9OhR9OzZExEREbU+yqakpAQdO3bEvHnz4OHhYZRtElmqw4cP4+bNm3Bzc0O3bt1ExyGZadeuHXr16gVw7gSRrMmq2C1atAgTJ05EVFQUunXrhuXLl8POzg4rV66scf3Q0FB8+OGHGD16NGxtbY2yTSJLpR+GfeCBB3h9HdVIPxxLRPIlm2JXXl6OI0eOIDw83PCaUqlEeHi4YXioubZZVlYGrVZbbSEyd7y+ju5FPxxLRPIlm2J3/fp16HQ6uLu7V3vd3d0dGo2mWbcZFxcHJycnw+Lt7d2gzycyFbdu3cK+ffsAsNhR7YKDg9HGvY3oGER0F7IpdnISExODwsJCw3L58mXRkYia1NGjR1FUVAQXFxf06NFDdBySKaVSiccee0x0DCK6C9kUOzc3N6hUKuTk5FR7PScnp9aJEU21TVtbWzg6OlZbiMyZfhh24MCBUCpl82OBZOjxxx43fM2nUBDJj2x+gtvY2CA4OBiJiYmG1yorK5GYmIh+/frJZptE5ojX11FdPfTQQ4av+RQKIvmxEh3gdtHR0Rg/fjxCQkLQp08fxMfHo7i4GFFRUQCAcePGoW3btoiLiwNQNTni5MmThq+vXLmCtLQ0ODg4oFOnTnXaJpGl0+l02LNnDwA+H5buzc7OzvD1z9t/RlivMIFpiOivZFXsRo0ahWvXrmHWrFnQaDQICgpCQkKCYfJDZmZmtWGiq1evVt1X6f8tXLgQCxcuxKBBg5CUlFSnbRJZumPHjkGr1aJly5YICgoSHYdMyPaft+O9f74nOgYR3UYh8SKJe9JqtXByckJhYSGvtyOz8/HHHyM6OhqPPfYYtm/n46Lo7orLi+EQ51D1m38BOZdz0KYNZ8oKVVwMOPz/MSkqAuztxeYho6tPD5HNNXZEJAavr6MGk8B/DBDJDIsdkQWrrKw0XF/HYkcN8eOPP4qOQES3YbEjsmAZGRnIz8+Hvb09goODRcchE7Rz506UlZWJjkFE/4/FjsiC6Ydh+/fvD2tra8FpyNS0cW+DoqIiw2Q1IhKPxY7IgvH6OmoM/VMoOBxLJB8sdkQWSpIk/PbbbwBY7Khh9E+h+PHHH/kUCiKZYLEjslCnT59Gbm4u1Go1QkNDRcchE/TQQw9BrVYjMzMT6enpouMQEVjsiCyWfhi2X79+sLW1FZyGTJGdnR0GDx4MgMOxRHLBYkdkofTFjo8Ro8Z46qmnALDYEckFix2RBZIkiRMnyCieeOIJAMCBAweQm5srOA0RsdgRWaBz584hOzsbNjY26Nu3r+g4ZMLatWuHXr16QZIk/PTTT6LjEFk8FjsiC6SfDdunTx+0aNFCcBoydRyOJZIPFjsiC8RhWDImfbHjUyiIxGOxI7JALHZkTL1794aXlxeKi4vx9ddfi45DZNFY7IgszMWLF5GZmQkrKyv0799fdBwyA0qlEpMnTwYAzJ8/nzcrJhKIxY7IwujP1oWEhMDe3l5wGjIXkydPhoODA9LT07F9+3bRcYgsFosdkYXhMCw1BRcXF7zyyisAgHnz5glOQ2S5WOyILAyLHTWV1157DdbW1ti7dy/27dsnOg6RRWKxI7IgWVlZ+OOPP6BUKjFgwADRccjMtG3bFuPGjQNQda0dETU/FjsiC6I/W9e7d284OjoKTkPmaPr06VAoFPjxxx+Rnp4uOg6RxWGxI7IgfD4sNbUuXbpg2LBhAIAFCxYITkNkeVjsiCwIr6+j5vDWW28BAL766itcunRJcBoiy8JiR2QhNBoNfv/9dygUCgwcOFB0HDJjoaGhePjhh6HT6bBo0SLRcYgsCosdkYXQPx82MDAQLi4ugtOQuXv77bcBACtWrMD169cFpyGyHCx2RBaCw7DUnMLDw9G7d2/cvHkTixcvFh2HyGKw2BFZCBY7ak4KhcJw1m7x4sUoKioSnIjIMrDYEVmA69evIyMjAwBnxFLzGTZsGDp16oQ///wT//73v0XHIbIILHZEFkB/fV337t3h5uYmOA1ZCpVKhenTpwMAPvroI5SXlwtORGT+WOyILMDPP/8MgMOw1PzGjRsHDw8PZGVl4auvvhIdh8jssdgRmblr165h3bp1AICRI0cKTkOWRq1W4/XXXwdQ9ZixyspKwYmIzBuLHZGZW7p0KUpLSxEaGsrr60iIV155BU5OTjh9+jS2bt0qOg6RWWOxIzJjJSUlWLJkCYD/PcOTqLk5Ojpi8uTJAIB58+ZBkiTBiYjMF4sdkRlbvXo18vLy0KFDBwwdOlR0HLJg06ZNg62tLQ4cOGCYzENExsdiR2SmdDodPvroIwBAdHQ0rKysBCciS+bu7o6oqCgAVWftiKhpsNgRmaktW7bg/PnzcHV1NfyFSiTSm2++CaVSiYSEBKSlpYmOQ2SWWOyIzJAkSfjwww8BAJMnT4a9vb3gRESAr6+vYWb2/PnzBachMk8sdkRmaN++fThw4ABsbW0xZcoU0XGIDN566y0AwObNm/HHH38ITkNkfljsiMyQ/mzd+PHj4e7uLjgN0f8EBQUhIiIClZWVWLhwoeg4RGaHxY7IzJw5cwZbt26FQqFAdHS06DhEd3j77bcBACtXrkROTo7gNETmhcWOyMzoZ8I+/fTT6NKli+A0RHcaNGgQwsLCUFZWhk8++UR0HCKzwmJHZEZycnKwZs0aAFUzEInkSKFQGM7affbZZ9BqtYITEZkPFjsiM7JkyRKUlZWhb9++GDBggOg4RLV6+umn0bVrVxQWFuLzzz8XHYfIbLDYEZmJ4uJifPbZZwD4+DCSP6VSiRkzZgAAFi1ahNLSUsGJiMwDix2RmVi5ciXy8/PRqVMnPPPMM6LjEN3T2LFj0bZtW2g0Gqxdu1Z0HCKzIMtit3TpUvj4+ECtViMsLAwHDx686/pff/01unbtCrVajR49emD79u3Vvv/8889DoVBUWyIjI5tyF4ia1a1bt/Dxxx8DqHp8mEqlEpyI6N5sbGzwxhtvAAAWLFgAnU4nOBGR6ZNdsdu0aROio6MRGxuLo0ePomfPnoiIiEBubm6N6+/fvx9jxozBhAkTkJqaiiFDhmDIkCFIT0+vtl5kZCSys7MNy4YNG5pjd4iaxXfffYcLFy7Azc0N48ePFx2HqM4mTpwIFxcXnDt3Dt99953oOEQmT3bFbtGiRZg4cSKioqLQrVs3LF++HHZ2dli5cmWN63/yySeIjIzE9OnT4e/vj9mzZ6N3795YsmRJtfVsbW3h4eFhWFxcXJpjd4ia3O2PD3v11VdhZ2cnOBFR3Tk4OGDq1KkAqh4zJkmS4EREpk1Wxa68vBxHjhxBeHi44TWlUonw8HAkJyfX+J7k5ORq6wNARETEHesnJSWhTZs26NKlCyZNmoS8vLxac5SVlUGr1VZbiOTqt99+w+HDh6FWq/Hqq6+KjkNUb1OnTkWLFi1w5MgRJCYmio5DZNJkVeyuX78OnU53xyOQ3N3dodFoanyPRqO55/qRkZFYs2YNEhMTMX/+fOzevRuPPfZYrddzxMXFwcnJybB4e3s3cs+Imo7+bF1UVBRat24tOA1R/bm5ueHFF18EAMybN09wGiLTJqti11RGjx6Np59+Gj169MCQIUOwbds2HDp0CElJSTWuHxMTg8LCQsNy+fLl5g1MVEcnT57ETz/9BIVCgddff110HKIGe+ONN6BSqZCYmIjDhw+LjkNksmRV7Nzc3KBSqe54dmBOTg48PDxqfI+Hh0e91geAjh07ws3NDefOnavx+7a2tnB0dKy2EMmR/vFhQ4cOhZ+fn+A0RA1333334W9/+xuAqmvtiKhhZFXsbGxsEBwcXO0ai8rKSiQmJqJfv341vqdfv353XJPxyy+/1Lo+AGRlZSEvLw+enp7GCU4kQHZ2NtatWweg6obERKZOf8Pib7/9Fr///rvgNESmSVbFDqi6B9eKFSuwevVqnDp1CpMmTUJxcTGioqIAAOPGjUNMTIxh/WnTpiEhIQEfffQRTp8+jffeew+HDx/GlClTAABFRUWYPn06UlJScPHiRSQmJuKZZ55Bp06dEBERIWQfiYzh008/RXl5OQYMGIC+ffuKjkPUaAEBAXjyySerzfQmovqRXbEbNWoUFi5ciFmzZiEoKAhpaWlISEgwTJDIzMxEdna2Yf3+/fvjq6++whdffIGePXvim2++wZYtWxAQEAAAUKlUOH78OJ5++ml07twZEyZMQHBwMPbs2QNbW1sh+0jUWDdu3MDy5csB8GwdmZe3334bALB69WpcvHhRbBgiE6SQeNOge9JqtXByckJhYSGvtyNZiI+Px+uvv47OnTvj1KlTUCpl9280MlPF5cVwiHMAABTFFMHext7on/Hggw9i9+7dCAkJwZ49e6BWq43+GWaluBhwqDomKCoC7I1/TEis+vQQ/m1AZGIqKioMjw974403WOrI7KxatQqurq44fPgwJk2axJsWE9UD/0YgMjHffPMNMjMz0aZNG4wbN050HCKj8/HxwaZNm6BUKrFq1SosW7ZMdCQik8FiR2RCbr+ofMqUKRyiIrMVHh5uuFnxtGnTsHfvXsGJiEwDix2RCfn111+RmpoKOzs7TJ48WXQcoib15ptvYuTIkbh16xZGjBiBK1euiI5EJHssdkQmZOHChQCAF154Aa1atRKchqhpKRQKrFy5Ej169EBOTg5GjBiBsrIy0bGIZI3FjshEnDhxAgkJCVAqlXx8GFkMe3t7fP/993B2dkZKSgr+8Y9/iI5EJGssdkQmQn+2bvjw4ejYsaPgNETNx9fXFxs2bIBCocAXX3yBL774QnQkItlisSMyAVlZWfjqq68A8IbEZJkiIyMxZ84cAFUTh5KTkwUnIpInFjsiE/Dpp5/i1q1beOCBBxAaGio6DpEQMTExGDZsGCoqKjBixAhoNBrRkYhkh8WOSOa0Wi0+//xzADxbR5ZNoVBg1apV8Pf3x9WrV/Hss8+ivLxcdCwiWWGxI5K5FStWQKvVwt/fH48//rjoOERCtWzZElu2bIGjoyP27t2L6Oho0ZGIZIXFjkjGKioqEB8fD4CPDyPS69y5M9avXw8AWLp0Kb788kvBiYjkg39LEMnYxo0bkZWVBXd3d/z9738XHYdINp588km8//77AIBJkybh0KFDghMRyQOLHZFMSZJkuMXJP/7xD9ja2gpORCQv7777Lp5++mmUlZVh2LBhyM3NFR2JSDgWOyKZ+uWXX3D8+HHY29tj0qRJouMQyY5SqcSaNWvQuXNnZGVlYeTIkaioqBAdi0goFjsiGcrMzMTLL78MAHjxxRfh4uIiOBGRPDk5OWHLli1wcHDA7t27MWPGDNGRiIRisSOSmatXr2Lw4MG4ePEiOnXqhHfeeUd0JCJZ8/f3x5o1awAA8fHxWLduneBEROKw2BHJSE5ODgYPHoxz587Bx8cHv/76K9q0aSM6FpHsDR06FO+++y4AYOLEiUhNTRWciEgMFjsimbh+/TrCw8Nx+vRptGvXDr/++iu8vb1FxyIyGe+99x4ef/xxlJaWYujQobh+/broSETNjsWOSAb+/PNPPProo0hPT4enpyd+/fVXdOjQQXQsIpOiUqmwbt06+Pr64tKlSxg9ejRu3bolOhZRs2KxIxJMq9UiMjISqampaN26NRITE+Hn5yc6FpFJcnFxwZYtW2Bvb4/ExEReo0oWh8WOSKCioiI8/vjjOHjwIFxdXfHf//4X/v7+omMRmbSAgADD0yg+/PBDbNq0SXAioubDYkckSElJCZ5++mns27cPTk5O+OWXXxAYGCg6FpFZePbZZ/HWW28BAF544QX89NNPghMRNQ8WOyIB9Bd379q1Cy1btsSOHTvQu3dv0bGIzMrcuXMRERGBkpISPPnkk3jppZdw48YN0bGImhSLHVEzKy8vx7PPPoudO3fCzs4O27dvR1hYmOhYRGZHpVLh+++/R3R0NBQKBVasWIGePXvit99+Ex2NqMmw2BE1o1u3bmHMmDHYtm0b1Go1fvzxR9x///2iYxGZrRYtWuCjjz7Cr7/+ivvuuw8XLlzAgw8+iOnTp6O0tFR0PCKjY7EjaiY6nQ7PPfccvvvuO9jY2GDLli14+OGHRccisggPPvggjh8/jhdeeAGSJGHhwoUICQnB0aNHRUcjMioWO6JmUFlZiQkTJmDjxo2wsrLCN998g4iICNGxiCyKo6Mj/vOf/2Dr1q1o06YNMjIyEBYWhjlz5vB+d2Q2WOyImpgkSZg0aRJWr14NlUqFjRs34qmnnhIdi8hiPfXUU0hPT8ewYcNw69YtzJw5E/fffz/OnDkjOhpRo7HYETUhSZIwbdo0fPHFF1AqlVi7di2GDx8uOhaRxWvdujW++eYbrF27Fk5OTjhw4AB69eqFJUuWoLKyUnQ8ogZjsSNqIpIkYcaMGVi8eDEAYOXKlRgzZozgVESkp1Ao8Pe//x0nTpxAeHg4bt68ialTp+LRRx/F5cuXRccjahAWO6ImMmvWLCxcuBAA8Pnnn2P8+PGCExFRTby9vbFjxw4sWbIELVq0QGJiInr06IG1a9dCkiTR8YjqhcWOqAnMmTMHc+bMAQB8+umneOmllwQnIqK7USqVePXVV5GWloawsDAUFhZi3LhxGD58OK5duyY6HlGdsdgRGdnChQsxc+ZMAFXPqZw6dargRERUV507d8bevXsxZ84cWFlZ4fvvv0dAQAB++OEH0dGI6oTFjshIsrKy8Prrr2P69OkAgNmzZ+PNN98UnIqI6svKygr//Oc/cfDgQXTv3h25ubkYMmQIoqKiUFhYKDoe0V2x2BE10sGDBzFmzBj4+PggPj4eAPDPf/4T7777rthgRNQovXr1wuHDhzF9+nQoFAqsWrUKgYGB+OSTT5CTkyM6HlGNWOyIGuDWrVv4+uuvMWDAAISFhWHjxo3Q6XR48MEH8eOPP2L27NmiIxKREajVaixYsAC7d+9Ghw4dkJmZiddeew1eXl6IjIzE2rVrcePGDdExiQxY7IjqoaCgAAsXLoSvry9GjhyJ/fv3w9raGuPHj8fRo0exa9cuPPnkk1AoFKKjEpERDRw4EMePH8enn36KsLAwVFZWYseOHRg3bhzc3d0xZswY/PTTT6ioqBAdlSycQuJc7nvSarVwcnJCYWEhHB0dRcchAc6dO4dPP/0UK1euRHFxMQDAzc0NkyZNwuTJk+Hh4SE4IVHzKC4vhkOcAwCgKKYI9jb2ghOJcfbsWXz11VdYv349zp49a3jdzc0NI0eOxNixY9GvX7/m+UdecTHgUHVMUFQE2FvmMTFn9ekhLHZ1wGJnmSRJwu7du/Hxxx/jxx9/NNzPKiAgAK+99hrGjh0LtVotOCVR82Kxq06SJBw+fBjr16/Hhg0bkJuba/hex44d8be//Q1jx45F165dmy4Ei53ZY7EzMhY7y1JWVoaNGzciPj4eaWlphtcff/xxvP766xg8eDCHWslisdjV7tatW0hMTMT69evx3XffGc7uA0BwcDDGjh2L0aNHw9PT07gfzGJn9ljsjIzFzjLk5uZi+fLl+Oyzzwwz3lq0aIHnn38e06ZNQ5cuXQQnJBKPxa5uiouLsXXrVqxfvx47duzArVu3AFTdCPnhhx/G3//+dwwdOtQ4f6ew2Jk9FjsjY7EzT5IkISsrCxkZGfjmm2+wbt06lJWVAQDatm2LqVOnYuLEiXB1dRWclEg+WOzq79q1a/j666+xbt06JCcnG163tbVFYGAgunfvblgCAgLQrl27+o0KsNiZPRY7I2OxM22SJEGj0SAjIwMZGRlIT083fK3VaqutGxoaitdffx0jRoyAtbW1oMRE8sVi1zh//PGHYdLF6dOna1zH0dER3bp1u6PweXp61lz4WOzMnskXu6VLl+LDDz+ERqNBz549sXjxYvTp06fW9b/++mvMnDkTFy9ehJ+fH+bPn4/HH3/c8H1JkhAbG4sVK1agoKAAAwYMwLJly+Dn51enPCx2puPatWs1Frj8/Pwa17eyskLnzp0RHByMV155pflmsRGZKBY745AkCWfPnsWJEyeq/az6/fffDcO2f+Xs7Fyt7OkLXxt7eyhatqxaicXOLJl0sdu0aRPGjRuH5cuXIywsDPHx8fj6669x5swZtGnT5o719+/fjwceeABxcXF48skn8dVXX2H+/Pk4evQoAgICAADz589HXFwcVq9ejQ4dOmDmzJk4ceIETp48WadZjSx28lFZWYmioiIUFBQgMzPzjgJ3+4y02ymVSnTq1KnaD8Pu3bujc+fOsLGxaea9IDJdLHZNq7y8HGfPnq32cy0jIwPnzp2DTqer8T3erq7I/P9/vEa/9BJaenigVatWdyyurq5wcnLiP15NkEkXu7CwMISGhmLJkiUAqv4i9/b2xtSpU/H222/fsf6oUaNQXFyMbdu2GV7r27cvgoKCsHz5ckiSBC8vL7zxxhuG53YWFhbC3d0dq1atwujRo++ZicXOeCRJMhSz25c///zzjtdqWgoLC1FZWXnXz+jYseMdBa5r1668NQmREbDYiVFWVoYzZ87cUfjOnz+PFpIE/fxbewAld9mOSqWCq6trjcXv9sXBwQFqtfqui7W1NUtiM6lPD7Fqpkx1Ul5ejiNHjiAmJsbwmlKpRHh4eLULTm+XnJyM6Ojoaq9FRERgy5YtAIALFy5Ao9EgPDzc8H0nJyeEhYUhOTm5TsWuqeXn52Pnzp2QJMlwrzT913/9fV2+1i+VlZWGX++23GudW7duoaKiAuXl5SgvL2/U1zdv3qz1X531YWNjAw8PjzsKnL+/P+w5DEFEZkY/0SIwMLDa6zdv3sTvqanAgAEAgBnTpyOnqAh5eXl3LCUlJdDpdLh27RquXbvW6EwKheKe5U9fAFUqlWGxsrKq9vv6fE+pVEKhUBiWv/6+Pq/pF/2+1PTr3b73118HDx4MNze3Rv+5Npasit3169eh0+ng7u5e7XV3d/daLzLVaDQ1rq/RaAzf179W2zp/VVZWZpgdCVSd4QNwx4X2xnLixAmMGTOmSbYtV1ZWVnB2doaTk5Ph19uXv77219/XdvZNp9M12XEioqozdiit+lqr1UJn0/h/qFHjdOjQAfqfeq+//nqt19jdvHkTf/75J/Lz8++53Lx5E6WlpSgrK6v2a3l5uWF7kiTh5s2buHnzZjPspfz98ssvd50P0Bj6v9fqMsgqq2InF3FxcXj//ffveN3b21tAGvN069YtXL9+HdevXxcdhYgayGuel+gI9FdePCaiPPLII03+GTdu3ICTk9Nd15FVsXNzc4NKpTLcHFYvJyen1mdxenh43HV9/a85OTnV7vadk5ODoKCgGrcZExNTbXi3srIS+fn5aNWqFa8nuAutVgtvb29cvnyZ1yKaKB5D08bjZ9p4/ExbUx4/SZJw48YNeNWhuMuq2NnY2CA4OBiJiYkYMmQIgKpSlZiYiClTptT4nn79+iExMRGvvfaa4bVffvkF/fr1A1B1itrDwwOJiYmGIqfVanHgwAFMmjSpxm3a2trC1ta22mvOzs6N2jdL4ujoyB9KJo7H0LTx+Jk2Hj/T1lTH715n6vRkVewAIDo6GuPHj0dISAj69OmD+Ph4FBcXIyoqCgAwbtw4tG3bFnFxcQCAadOmYdCgQfjoo4/wxBNPYOPGjTh8+DC++OILAFUXNb722muYM2cO/Pz8DLc78fLyMpRHIiIiInMgu2I3atQoXLt2DbNmzYJGo0FQUBASEhIMkx8yMzOhVCoN6/fv3x9fffUV3n33Xbzzzjvw8/PDli1bDPewA4AZM2aguLgYL730EgoKCnD//fcjISGBt78gIiIisyK7+9iR6SorK0NcXBxiYmLuGMom08BjaNp4/Ewbj59pk8vxY7EjIiIiMhPKe69CRERERKaAxY6IiIjITLDYEREREZkJFjuqZunSpfDx8YFarUZYWBgOHjx41/W//vprdO3aFWq1Gj169MD27dtrXfeVV16BQqFAfHx8tdfz8/MxduxYODo6wtnZGRMmTEBRUZExdsfiiDh+c+fORf/+/WFnZ8f7PTZScx+/ixcvYsKECejQoQNatGgBX19fxMbGVntsFNWdiP//nn76abRv3x5qtRqenp547rnncPXqVWPsjsURcfz0ysrKEBQUBIVCgbS0tEbsBYsd3WbTpk2Ijo5GbGwsjh49ip49eyIiIgK5ubk1rr9//36MGTMGEyZMQGpqKoYMGYIhQ4YgPT39jnW///57pKSk1HjX7LFjxyIjIwO//PILtm3bht9++w0vvfSS0ffP3Ik6fuXl5Xj22WdrveE31Y2I43f69GlUVlbi888/R0ZGBj7++GMsX74c77zzTpPsozkT9f/fQw89hM2bN+PMmTP49ttvcf78eYwYMcLo+2fuRB0/vRkzZtTpqRJ1IhH9vz59+kivvvqq4fc6nU7y8vKS4uLialx/5MiR0hNPPFHttbCwMOnll1+u9lpWVpbUtm1bKT09Xbrvvvukjz/+2PC9kydPSgCkQ4cOGV77+eefJYVCIV25csUIe2U5RBy/23355ZeSk5NTo/bBkok+fnoLFiyQOnTo0LCdsGByOX4//PCDpFAopPLy8obtiIUSefy2b98ude3aVcrIyJAASKmpqY3aF56xIwBVZ12OHDmC8PBww2tKpRLh4eFITk6u8T3JycnV1geAiIiIautXVlbiueeew/Tp09G9e/cat+Hs7IyQkBDDa+Hh4VAqlThw4EBjd8tiiDp+ZBxyOn6FhYVwdXVtwF5YLrkcv/z8fKxfvx79+/eHtbV1A/fG8og8fjk5OZg4cSLWrl0LOzs7I+wNh2Lp/12/fh06nc7whA89d3d3aDSaGt+j0Wjuuf78+fNhZWWFf/zjH7Vuo02bNtVes7Kygqura62fS3cSdfzIOORy/M6dO4fFixfj5ZdfruceWDbRx++tt96Cvb09WrVqhczMTPzwww8N3BPLJOr4SZKE559/Hq+88kq1kxuNJbtHipH5OHLkCD755BMcPXoUCoVCdByqJx4/01bf43flyhVERkbi2WefxcSJE5shId1NfY7f9OnTMWHCBFy6dAnvv/8+xo0bh23btvH/W4HqcvwWL16MGzduICYmxqifzTN2BABwc3ODSqVCTk5OtddzcnLg4eFR43s8PDzuuv6ePXuQm5uL9u3bw8rKClZWVrh06RLeeOMN+Pj4GLbx14tTb926hfz8/Fo/l+4k6viRcYg+flevXsVDDz2E/v3744svvjDejlkI0cfPzc0NnTt3xiOPPIKNGzdi+/btSElJMd4OmjlRx+/XX39FcnIybG1tYWVlhU6dOgEAQkJCMH78+AbvD4sdAQBsbGwQHByMxMREw2uVlZVITExEv379anxPv379qq0PAL/88oth/eeeew7Hjx9HWlqaYfHy8sL06dOxY8cOwzYKCgpw5MgRwzZ+/fVXVFZWIiwszNi7abZEHT8yDpHH78qVK3jwwQcRHByML7/8Ekol/1qoLzn9/1dZWQmg6vYZVDeijt+nn36KY8eOGb6vv13Kpk2bMHfu3IbvUKOmXpBZ2bhxo2RrayutWrVKOnnypPTSSy9Jzs7OkkajkSRJkp577jnp7bffNqy/b98+ycrKSlq4cKF06tQpKTY2VrK2tpZOnDhR62fUNCsoMjJS6tWrl3TgwAFp7969kp+fnzRmzJgm2UdzJur4Xbp0SUpNTZXef/99ycHBQUpNTZVSU1OlGzduNMl+misRxy8rK0vq1KmTNHjwYCkrK0vKzs42LFQ/Io5fSkqKtHjxYik1NVW6ePGilJiYKPXv31/y9fWVSktLm2xfzZGon5+3u3DhglFmxbLYUTWLFy+W2rdvL9nY2Eh9+vSRUlJSDN8bNGiQNH78+Grrb968WercubNkY2Mjde/eXfrpp5/uuv2a/sPOy8uTxowZIzk4OEiOjo5SVFQUS0EDiTh+48ePlwDcsezatctIe2U5mvv4ffnllzUeO/6bv2Ga+/gdP35ceuihhyRXV1fJ1tZW8vHxkV555RUpKyvLmLtlMUT8/LydsYqdQpIkqeHn+4iIiIhILngxBREREZGZYLEjIiIiMhMsdkRERERmgsWOiIiIyEyw2BERERGZCRY7IiIiIjPBYkdERERkJljsiIiIiMwEix0RERGRmWCxIyIiIjITLHZEREREZoLFjoioifzrX/+CQqG4Y4mPjxcdjYjMlEKSJEl0CCIic3Tjxg0UFxcbfj9r1izs3LkTe/fuRbt27QQmIyJzZSU6ABGRuWrZsiVatmwJAJg5cyZ27tyJpKQkljoiajIciiUiamKzZs3C2rVrkZSUBB8fH9FxiMiMsdgRETWh2NhYrFmzhqWOiJoFix0RUROJjY3F6tWrWeqIqNnwGjsioiYwZ84cLFu2DFu3boVarYZGowEAuLi4wNbWVnA6IjJXnBVLRGRkkiTB2dkZWq32ju8dPHgQoaGhAlIRkSVgsSMiIiIyE7zGjoiIiMhMsNgRERERmQkWOyIiIiIzwWJHREREZCZY7IiIiIjMBIsdERERkZlgsSMiIiIyEyx2RERERGaCxY6IiIjITLDYEREREZkJFjsiIiIiM8FiR0RERGQm/g89P7deBI6rFAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "mask = post > 1e-5\n",
    "plt.plot(zcos[mask], post[mask], color=\"black\", label=r\"$p(z_{\\rm cosmo})$\")\n",
    "# plt.plot(zcos[mask], post_bad[mask], color=\"black\", ls=\"dashed\")\n",
    "\n",
    "\n",
    "plt.ylim(0)\n",
    "plt.axvline(zcosmo_mean[n], color=\"green\", label=r\"$z_{\\rm cal}$\")\n",
    "# plt.fill_betweenx([0, plt.ylim()[1]], zcosmo_mean[n] - zcosmo_std[n],\n",
    "                #   zcosmo_mean[n] + zcosmo_std[n], color=\"green\", alpha=0.2)\n",
    "plt.axvline(loader.cat[\"zobs\"][n], color=\"red\", label=r\"$z_{\\rm CMB}$\")\n",
    "\n",
    "plt.xlabel(r\"$z$\")\n",
    "plt.ylabel(r\"$p(z)$\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"../plots/zcosmo_posterior.png\", dpi=450)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_csiborg",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}