{
 "cells": [
  {
   "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",
    "from os.path import exists\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from corner import corner\n",
    "from getdist import plots\n",
    "from astropy.coordinates import angular_separation\n",
    "import scienceplots\n",
    "from os.path import exists\n",
    "import seaborn as sns\n",
    "import matplotlib as mpl\n",
    "\n",
    "\n",
    "from reconstruction_comparison import *\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
    "fdir = \"/mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = paths.flow_validation(\n",
    "    fdir, \"Carrick2015\", \"2MTF\", inference_method=\"mike\",\n",
    "    sample_alpha=True, sample_beta=None,\n",
    "    zcmb_max=0.05)\n",
    "\n",
    "\n",
    "get_gof(\"neg_lnZ_harmonic\", fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogue = \"CF4_TFR_i\"\n",
    "simname = \"Carrick2015\"\n",
    "zcmb_max=0.05\n",
    "sample_beta = None\n",
    "sample_alpha = True\n",
    "\n",
    "fname_bayes = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"bayes\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "fname_mike = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"mike\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "\n",
    "X = []\n",
    "labels = [\"Full posterior\", \"Delta posterior\"]\n",
    "for i, fname in enumerate([fname_bayes, fname_mike]):\n",
    "    samples = get_samples(fname)\n",
    "    if i == 1:\n",
    "        print(samples.keys())\n",
    "\n",
    "    X.append(samples_to_getdist(samples, labels[i]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = [f\"a_{catalogue}\", f\"b_{catalogue}\", f\"c_{catalogue}\", f\"e_mu_{catalogue}\",\n",
    "          \"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"alpha_{catalogue}\"]\n",
    "# params = [\"beta\", f\"a_{catalogue}\", f\"b_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "# params = [\"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"mag_cal_{catalogue}\", f\"alpha_cal_{catalogue}\", f\"beta_cal_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "\n",
    "\n",
    "g = plots.get_subplot_plotter()\n",
    "g.settings.figure_legend_frame = False\n",
    "g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "plt.gcf().suptitle(catalogue_to_pretty(catalogue), y=1.025)\n",
    "plt.gcf().tight_layout()\n",
    "# plt.gcf().savefig(f\"../../plots/method_comparison_{simname}_{catalogue}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# catalogue = [\"LOSS\", \"Foundation\"]\n",
    "catalogue = \"CF4_TFR_i\"\n",
    "simname = \"IndranilVoid_exp\"\n",
    "zcmb_max = 0.05\n",
    "sample_alpha = False\n",
    "\n",
    "fname = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"mike\",\n",
    "    sample_mag_dipole=True,\n",
    "    sample_beta=False,\n",
    "    sample_alpha=sample_alpha, zcmb_max=zcmb_max)\n",
    "\n",
    "\n",
    "samples = get_samples(fname, convert_Vext_to_galactic=True)\n",
    "\n",
    "samples, labels, keys = samples_for_corner(samples)\n",
    "fig = corner(samples, labels=labels, show_titles=True,\n",
    "             title_kwargs={\"fontsize\": 12}, smooth=1)\n",
    "# fig.savefig(\"../../plots/test.png\", dpi=250)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paper plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0. LOS velocity example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpath = \"/mnt/extraspace/rstiskalek/catalogs/PV/CF4/CF4_TF-distances.hdf5\"\n",
    "\n",
    "loader_carrick = csiborgtools.flow.DataLoader(\"Carrick2015\", [0], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_lilow = csiborgtools.flow.DataLoader(\"Lilow2024\", [0], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_cb2 = csiborgtools.flow.DataLoader(\"csiborg2_main\", [i for i in range(20)], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_cb2X = csiborgtools.flow.DataLoader(\"csiborg2X\", [i for i in range(20)], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_CF4 = csiborgtools.flow.DataLoader(\"CF4\", [i for i in range(20)], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_CLONES = csiborgtools.flow.DataLoader(\"CLONES\", [0], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angdist = angular_separation(\n",
    "    np.deg2rad(loader_carrick.cat[\"RA\"]), np.deg2rad(loader_carrick.cat[\"DEC\"]),\n",
    "    np.deg2rad(csiborgtools.clusters[\"Virgo\"].spherical_pos[1]),\n",
    "    np.deg2rad(csiborgtools.clusters[\"Virgo\"].spherical_pos[2]))\n",
    "k = np.argmin(angdist)\n",
    "print([loader_carrick.cat[\"RA\"][k], loader_carrick.cat[\"DEC\"][k]])\n",
    "print(csiborgtools.clusters[\"Virgo\"].spherical_pos[1:])\n",
    "print(csiborgtools.clusters[\"Virgo\"].spherical_pos[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loaders = [loader_carrick, loader_lilow, loader_CF4, loader_cb2, loader_cb2X, loader_CLONES]\n",
    "simnames = [\"Carrick2015\", \"Lilow2024\", \"CF4\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\"]\n",
    "\n",
    "\n",
    "with plt.style.context(\"science\"):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    plt.figure()\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "    for i, (simname, loader) in enumerate(zip(simnames, loaders)):\n",
    "        r = loader.rdist\n",
    "        vrad = loader.los_radial_velocity[:, k, :]\n",
    "\n",
    "        if simname == \"Carrick2015\":\n",
    "            vrad *= 0.43\n",
    "\n",
    "        if len(vrad) > 1:\n",
    "            ylow, yhigh = np.percentile(vrad, [16, 84], axis=0)\n",
    "            plt.fill_between(r, ylow, yhigh, alpha=0.66, color=cols[i],\n",
    "                             label=simname_to_pretty(simname))\n",
    "        else:\n",
    "            plt.plot(r, vrad[0], label=simname_to_pretty(simname), c=cols[i])\n",
    "\n",
    "    plt.xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "    plt.ylabel(r\"$V_{\\rm rad} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "    plt.xlim(0, 90)\n",
    "    plt.ylim(-1000, 1000)\n",
    "    plt.legend(ncols=2, fontsize=\"small\")\n",
    "    plt.axvline(12.045, zorder=0, c=\"k\", ls=\"--\", alpha=0.75)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/LOS_example.pdf\", dpi=450, bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Evidence comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:02:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:17:56\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:24:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:41:00\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:14:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:52:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:33:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:03:00\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:18:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:25:22\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:43:31\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:17:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:55:20\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_LOSS_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:36:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 02/10/2024 17:47:56\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:18:31\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:25:04\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:46:23\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:20:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:58:23\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:38:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:04:09\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:18:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:25:46\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:56:16\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:23:54\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:01:40\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_Foundation_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:40:28\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:04:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:18:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:33:12\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:58:05\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:37:40\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:12:23\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:42:56\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:05:03\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:20:08\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:34:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:05:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:40:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:22:45\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:45:29\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:05:44\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:20:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:35:56\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:07:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:44:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:27:38\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:48:03\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:06:41\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:21:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:38:16\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:11:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:50:18\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:31:36\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:50:51\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:07:54\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:22:46\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:50:31\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:26:06\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:20:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:59:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:52:35\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:09:23\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:23:08\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:09:25\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:41:38\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:20:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:58:36\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:55:16\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:07:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:22:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:54:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:27:36\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:12:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:48:05\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:53:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:08:41\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:23:09\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 11:58:58\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 12:30:30\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:17:33\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:50:41\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_w1_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 26/09/2024 13:55:53\n"
     ]
    }
   ],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"manticore_2MPP_N128_DES_V1\", \"CF4\", \"CLONES\"]\n",
    "catalogues = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "y_BIC = np.full((len(catalogues), 2, len(sims)), np.nan)\n",
    "y_lnZ = np.full_like(y_BIC, np.nan)\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    for j in [0, 1]:\n",
    "        for k, simname in enumerate(sims):\n",
    "            fname = paths.flow_validation(\n",
    "                fdir, simname, catalogue, inference_method=\"mike\",\n",
    "                sample_alpha=simname != \"IndranilVoid_exp\",\n",
    "                zcmb_max=zcmb_max, smooth=j)\n",
    "\n",
    "            # y_BIC[i, j] = get_gof(\"BIC\", fname)z\n",
    "            y_lnZ[i, j, k] = - get_gof(\"neg_lnZ_harmonic\", fname) / np.log(10)\n",
    "\n",
    "        y_lnZ[i, j] -= y_lnZ[i, j].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzcAAAHfCAYAAABtWSYcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC//ElEQVR4nOzdeVxUVf8H8M+wb8IwuKS4MbibpgPY4pLK4FqWykCYZouAW9svBWl5nsoUhxazXBqw3SJkNLPFjNFKsyxgNMsljcElcmUYQRYRuL8/eJgktkFmgeHzfr18Pc6dM3O+jzcO93vPud8jEgRBABERERERURvnYOsAiIiIiIiIzIHJDRERERER2QUmN0REREREZBeY3BARERERkV1gckNERERERHaByQ0REREREdkFJjdERERERGQXmNwQEREREZFdYHJDRERERER2gckNERERERHZBSY3RERERERkF5jcEBERERGRXXCydQBE9k4QBOTn5+PKlSvw8vKCn58fRCKRrcOi6/ActS08X20DzxMR2QJnbogsxGAwYM2aNRg4oC86deqEgIAAdOrUCQMH9MWaNWtgMBhsHWK7V3OOBgwYUOscDRgwoNWcI4VCYesQWg3j+erfr/b56t/PpueL56i2mvPUv3/tsa9/f459RGR5IkEQBFsHQWRvdu7ciYiImSgpKcHMMGDmBAG+3kBBIbDlGxG2ZAAeHh7YvHkLJk6caPF4DAYDEhMToVQqLd5XW1F9jiJQUlKCMWPGYMyYMfDy8sKVK1ewZ88e7Nmz53/naLNVzlFDFAoF0tPTbdZ/a7Fz505EKMJRUlKC8QPFGD/AB97ujigsrcTuY5ex+6ih+nylq61+vniO/rFz504oFOEoLilB4FBPSId4wNXDAVdLqqD7rQQ5h4rh6eGBdBucJyJqH5jcEJnZzp07cdddUzFxpICNy6twU6e6bc5dBOY954Cd+0T44osvzfJLXqPRQKVSGS+ykpOTkZOTA4PBgMDAQKSlpSE2NhYxMTE39P0ikQj2MlxUn6O7EBISgiVLlsDPz69Om/z8fLzyyivIzMzEF198YZZzlJycbEwww8PDoVQqjedNIpFALBZDqVRCp9NBoVBALpdDo9EgOzu72X3Z3/maitukHfDc3T3Q0cu5TptLV65h+ednsF9XZJafKbVajcTERBgMBsjlcqhUqlrnKj4+vsXnCLDP89S9nzvGRvjB07vuyvfiwgp8tzkffx0vbfF5CgoKgsFggFQqBQDo9XokJCQgPDz8hr+zKQaDAb6+vigoKIBYLG5xOyIyPz5zQ2RGBoMBEREzMXGkgG1rq+DUwE/YTZ2AbWurcO9iB0REzMSpU3+16BegQqGAwWCATqczHsvJyQEAxMfHQyKRID8//4YTG3tSfY4iEBISguXLl8OpgZPk5+eH5cuX47nnnkNERAROnTrVonOk0+mQkZFhTDiDgoIQGBgIqVSKkJAQxMXFISwsDACMSapUKjUea68MBgMiFOG4TdoBr0YGwMmh/mc2Ono549XIADyVlosIRThOnT5zw+fLYDAgLS3NmLAEBgZCrVYbv0+lUiE+Pp7n6DoGgwEKRTi693PH5Ic6w8Gx/vPk6e2EyQ91xo53L0ChCMfpFpwnAFAqlcZkpubnSiqVQiaT3fB3NkYsFiM7O7vJmE1tR0Tmx2duiMzo/fffR0lJCTYubzixqeHkBKS8WIWSkhJ88MEHLeo3PT29zpIzpVKJyMhIY+JD1WrO0ZIlSxpMbGo4OTlhyZIlZjlHBoMBCQkJAKovfGJjY5GRkWF8TfWrOV/P3d2jwcSmhpODCM/e1aPF50ssFtdaZlYzMwCAiUwD3n//fRSXlGBshF+DiU0NB0cRxkb4odgMP1fXE4vFxpk0SzI1cbJUgkVEjWNyQ2QmgiBgw/o3MTMM9S5Fq0/XzsAMObB+3RtmX5qSlJQElUoFqVSK3r17Q6vVIj4+vtHPBAYGIj4+HoGBgfD19YVWq623nUKhQGBgIAIDA2tdSGg0GuPxmiUjQPUSn5rjarXabP8fm0sQBKxfvx6jR4+udylaffz8/DB69GisW7euRedIJpPVudiRSCT1to2NjUVsbCySk5Mb/c52cb7WrcX4AeJ6l6LVp1MHZ4wb4IN1a980y8+UWq1GVlYW5HJ5reOmniOgfZyndeveROBQz3qXotXH09sJ0iGeWLvWfGOfRqPB5s2bjTM5Df37+Pr6IikpyTh7qtFooFAo4OvrW2uMbOh8XF/x7d/fdf3suantkpOTERQUhKCgICgUilrnmIhugEBEZnHx4kUBgLB5NQThqOl/0l6DAEC4dOlSi/rPzs4WpFJpi74DgJCRkSEIgiAolcpa33f9cFHTJjs7WxCLxYIgCEJOTo4gFouFgoIC4+vr46o5LpfLWxRjS9Sco+eff1747rvvTP7z3//+1yznqEZBQYEglUqN/0Y3qr2cr1XhAUL2f2Um/0kMDzDL+ZLL5YJYLBaUSmWLvqe9nKeJD3QRFr0WaPKfCQ90adF5kslkglgsNv6vXC43/ns09u8DQEhPTxcEQRDCw8ON/9YFBQVNno+azzf0Xdf/t2JKu4KCAuN3Z2dn2/Q8EtkLztwQmcmVK1cAAL7ezftcTfuioiIzR3Rjau5Qx8XFQafT1XsHsaaNTCaDwWCAwWCAWq1GRESEcYlVzVKetLQ0ANV3QcPCwqDT6Sy+bKQhNefIy8urWZ/r0KEDAPOdo5rqWtcvd7pR7eF8ebs7Nutz3m7V7Vt6vjIyMlBQUICMjAyTZmga0x7Ok6tH8y4p3Nyr27fkPKWkpCA7Oxvp6enQ6XTGf6em/n1q/q1DQkKMfxeLxRCLxcZzU9/5qM/131XzrGNz2l0/g6vX6039v05EDWBBASIzqblgLihs3udq2tdcQLc29T0PkpSUhMzMTJM+7+fnZ6wMZms156jmYsxUNRdf5jhHCoUCSqXSog88/1tbP1+FpZXN+lxhWXV7c/1MhYWFQaVSmbUghz2ep6slVc36XFlpdXtznCe5XA6pVIqkpCTExcU1+e9z/b9/Q8tDTT0fpj4zV1+7moQqKCgIQHWyRkQtw5kbIjPx8/ND/36B2PJN83bg3vKNCP37BTb4C9baap4HSEpKqvcCXK1WIyMjA+np6bV+EYeHh2Pz5s3Gu5s1d6fDw8OhVqvrHLcFPz8/9OvXD3v27GnW5/bs2YN+/fq1+BwpFAokJCSYNbGx+/PVtw92H7vcrM/tPmpAv759bvh8qdXqWs9oZGRk1Hnmprns/Tz17RsI3W8lzfqc7lAJ+vY139inUqmM5btb+u/T0Pkwt5pS1tnZ2cjOzmYRAiIzYHJDZCYikQgLFj6KLRnV+9iY4uwFYKsGWLjosVoPnzZXbGwsFAoFdDodgoKCWvRwsUqlQlBQENLS0urdmFAulxv7iY6ONi6TkUqlSElJMT4wGxsbC71eD6lUCqVSidDQUAQFBSE2NvaGY2spkUiEhQsXYs+ePcjPzzfpM/n5+di7dy8WLVrUonNUc8EcGhoKX19f+Pr6mmVne7s/X4sWY/dRAy5duWbSZy4WXcO3xy5j0eJHb/h8hYeHIzMzs1a57pbOkNj7eVq06FHkHCpGcWGFSZ8pLqyA7rdiLF7csrHvelKpFBEREYiPj2/xv09D58PcxGIx1Go1RCIRfH19ERQU1GDBCSIyDTfxJDIjg8GAXr26Y/Tw0kb3uQGAigrg3sUO2HvAvcX73JiLPW0o2JDqc9QLgwcPbnSfGwCoqKjAc889h8OHD7d4nxtLaDfnq2cPDO3q1Og+NwBQUSXgqbRcHDpb0aJ9bsytvZynnj17wK8HGt3nBgCqKgV8/e4FXDqDFu9z09bFx8fDz88PcXFxAKorvsXGxjb67A4RNY4zN0RmJBaLsXnzFuzcJ8K9ix1w9kL97c5eqE5sdu4TIT19a7v+5W5t1edoMzIzM/Hcc881OIOTn5+P5557DpmZmUhPT+c5shGxWIzN6Wrs1xXhqbRcXCyqfwbnYtE1PJWWi/26IqSrt/B8WVn13kBq/HW8FDvevdDgDE5xYQV2vHsBZ46XQs3zVIelZoiI2hPO3BBZwM6dOxERMRMlJSWYIQdmThDg611dPGDLNyJs1QAeHh5IT9+KCRMm2Dpco/Zwh7lG9TmKQElJCUaPHo0xY8agQ4cOKCoqwp49e7B3797/naP0VnWOrtfuzpciHCUlJRg3wAfjB4rh7eaIwrJK7D5qwLfHLlefL/WWVne+2tt5UijCUVxSAukQT0iHesDN3QFlpVXQHSqB7rdieHp4QN0Kz5OtxMbGQqPRGBO9lJQUPntD1AJMbogsxGAw4IMPPsD6dW/gj+P/LDHo3y8QCxc9hrlz58LHx8eGEVLNOVq3bh2OHz9uPN6vXz8sWrSI56iVMZ6vtW/i+Ik/jcf79e2DRYsf5flqJWrO09q1b+DEiX/Gvr59A7F4Mcc+IrIsJjdEFiYIAvR6PYqKitChQwdIJBKzPUBL5sFz1LbwfLUNPE9EZAtMboiIiIiIyC6woAAREREREdkFJjdERERERGQXmNwQEREREZFdYHJDRERERER2gckNERERERHZBSY3RERERERkF5jcEBERERGRXWByQ0REREREdoHJDRERERER2QUmN0REREREZBeY3BARERERkV1oF8lNYGAg1Gp1g+/rdDqEhYUhMDAQQUFBSE5OrvW+VqtFUFAQfH19ERgYiPj4eJPeI6K24/qf45o/CoXC4v0mJyff0LhhMBggEolgMBjMHxQRNSo5Odk4Tlz/8xsWFgZfX996PxMUFITAwEAAQFJSUq2x5t9jj1arNX7GFuMSUVvmZOsAbE2r1SI0NBTp6emQy+UwGAxQKBTIzs6GSqWCwWCo8/7mzZsBoNH3iKhtMRgMSElJQXh4uK1DqZdOp4NUKjW+FovFyM7Ohlgstl1QRO2QTqdDRkYGcnJyYDAYjElLTEwMAEAikUCtVtcaS3Q6Xa0bEXFxcYiLiwPwz03SgoKCOn219nGJqDVqFzM3jYmOjkZKSgrkcjmA6guGjIwMbN68GVqtFllZWQBQ6/2aAayx94iIzCk+Pt54N7eGTCazUTRE7ZfBYEBCQgKA6t/7sbGxyMjIML4fHh4OlUpV6zMqlQqxsbFWjZOovWrXyY3BYIBWq633johcLkdaWpoxcVEoFNBoNHXaNPQeEdkPtVqNoKAgBAUF1bpAEYlExr8nJSXVes/X1xdJSUnGu7o6na5WW19fX4SFhSE9Pb1WXwqFwrj8pGZciY2NhVqthkKhQFBQUL39NxRjY3EQUfPJZLI6NxYkEonx72FhYXVmav49k0NEltOuk5usrKwGl3RIpVLjRUBubi6A6osOX1/fWs/vNPYeEbUt8fHxxgQhKCgIOp0OWq0WiYmJ2LVrl3EZmCnPyBgMBkilUmRnZ0MmkxnHhprvy83NRUZGRq2lZkB1IpOTk4P09HTj2nqVSgW5XI709HRkZ2fX6auxGBuKg4hazmAwQKVS1RkTwsPDjcvUNRqN8WbojahvXCKihrXr5CY4OLjBh3GvX98uFouRnp6OgoICpKSkQKFQGAeXxt4jorZFqVQiOzvb+EcqlSItLQ2RkZHGGyEJCQl1io40pOaCJiQkBDk5OQCqL3RiYmKM31fzgPG/PyOTyWAwGEwqGNBUjPXFQUQtp1AokJ6eXu9NipqlaS1dklbfuEREDWvXyY1YLIZUKq1zoWIwGKBWqxEZGVnnM+Hh4ZDJZHXWvjf1HhG1PzfysH9SUhIUCoVZKyKx6ACR+SkUCiiVynqffatJQLRaLbRaLZ+PI7Kidp3cAEB6ejri4+ONa9trqqXFxMRAJpNBo9EgKSnJePdUo9FAp9NBLpc3+h4R2YfIyEhj5USgeolIRESE8f2a49c/UNwYuVwOtVpd7+fUajUyMjKQnp6OlJSUWp8Ti8V11vGbGiMRmZdCoUBCQkKjSUtsbCwUCgULCRBZWbtJbqKjo2vVia9Zdy6TybBr1y4olUrjPjdhYWHG6eTg4GDk5+cjICAAvr6+iI+Px65duyAWixt9j4jsg0wmg1KprLVHRc34EBMTg9DQUCgUCkilUpN+9mUyGcLDwxEQEGD8nJ+fH4DqxEen0yEoKAjR0dG1lp+EhYUhOjoaoaGhdRKcxmIkIvNSq9VQq9UIDQ2Fr68vfH19651pjYmJgVQqZRVVIisTCYIg2DoIIiIiIiKilmo3MzdERERERGTfmNwQEREREZFdYHJDRERERER2gckNERERERHZBSY3RERERERkF5jcEBERERGRXXCydQDmNHjwYOMeD9T2FBQUIC8vD2VlZXBwcICPjw969OgBFxcXW4dGNygnJweHDx+2dRh1cKygG3Hp0iX8/fffuHr1KhwdHeHr64sePXrAycmufpXaRGscKzhO0I0oLS3FmTNnUFRUBEEQ4OHhgR49eqBDhw62Dq3NM3WcsKsROTAwENu3b7d1GNRMhw4dwsMPP4Ds7F/Ro6sTbg6oxNVrDtj/61nk5JxA9LwYvPraa0xy2qBp06bZOoR6cayg5vjhhx8QE/Mwjh49gcCeTugfWIkrJQ74+dDfOHkyB08+uQTPP/88HBy4GOJGtcaxguMENUdRUREWLIjFV199CTdXEUYMEeDoCBz6wwHf5eRg3LgxePfdD9CrVy9bh9pmmTpO2FVyQ21PdnY2xo+/EwHdyvDFBmDS6Ao4OgJAJS7qgeTNVXhxwwbk5JzAZ9u/gLOzs61DJqJ2RKPR4K67piB4cCV2vweMHVEBkQgAKvHXOWDtR6VYseIlnDp1Eu+99z5E1W8SUTty5coVjB8/Bsf/+A2rl1Vh7r2At1f1e+XlVfhUAySs/hG33x6Cfft+RkBAgE3jtXe8zUQ2U15ejunT78aA3mXYu6kSU8fif4lNtU4S4Jn5wBfrq6DZpUFiYqLNYiWi9qegoADh4dMxbkQldr9bhXG3AtfnLt1vAlY9BXyUJOCDDz5ESkqK7YIlIpt54vHH8cex3/D9+5V4dPY/iQ0AuLgAkVOAn1Ir4OlSgIiIGRAEwXbBtgNMbshmPv30U5w5cxbvvFSJDp4NtwsbCcRGVGH9+jdQXl5uvQCJqF17//33UVJSjHdXVKGxVbH3TQVmThBhzZpXedFC1M5cunQJmz76EM/EVGLYwIbbdekIrH22AllZB/HTTz9ZL8B2iMkN2UxKylu4c4QDBvdtuu2C+4Dz5/Px+eefWz4wIiIAKckbMEMO3NSp6bYLowQcOXIcP/74o+UDI6JWY9OmTRCESjwS3nTbsDuAwJ5OSE5Otnxg7RiTG7KZ438cw5igKpPaDuoDSMROOHHihIWjIiICBEHA8RM5GBNs2kzM6KDq/+UYRdS+HD9+HIMCHdHRt+m2Dg7AKFkFThw/avnA2jEWFCCbae7yDZGo+Z8x1dWrV3HixAmUlpaiU6dO6N27t0X6ISL7xmVprYdWq4VOp4Ner4fBYEBcXByA6iIROp0OEokEOp3OeJzsh8FgwPvvv49t27ahoKAAXl5ekMvliI6Ohr+/v63D4zhhYZy5IZvp06cffjxo2n+Cf+QC+QUVZt9z4MyZM1i2bBn8/f0xZMgQjBgxAgEBARgxYgTef/99XLt2zaz9EVHbIBKJECjtjR8Pmtb+p/+1474orYdGo0F4eDhiYmIgl8uhVqsBAEqlEjExMQgPr15HVHPcGqqqqlBQUID8/HxUVlZard/2ZM2aNejWrRuWLFmC8vJy9O7dGy4uLkhKSkKvXr3wxBNPmPXfPjAwEMd0ldAbmm4rCMD+Q04I7NPfbP1TXUxuyGYemReDXT9V4Y/cptu+9QnQsaPYrHsh7Nu3D0OHDsW6deswbtw4vPHGG0hJScHzzz8PkUiEBx98EFOnTkVxcbHZ+iSitmNe9Hyk73TARX3TbdenitCvrxSjRo2yfGBkkrS0NBgMBgCATqcDUJ3wiMViYxuZTIa0tDSLx3LmzBk888wz6NqlMyQSCTp27IiOfhI8+eSTOH78uMX7by9WrlyJJ554ApMmTUJaWhpWrlyJJ598Ev/973+hVqvxyCOPYO3atXjwwQfNNnsye/ZsVFY54L1tTbfdvR/4Q1eBefPmmaVvqh+TG7IZhUKBrl07Yd5zjigpbbjdnkxgQ5oDYmMXwc3NzSx9//HHH5g6dSp69+6N1NRULF68GEOHDkXfvn0xduxYJCUl4bXXXsO+fftw3333cQqZqB166KGH4OLihuj/iFBR0XC7TzOAzV8LePSxJ7mRZysSHByMgIAAJCcnQ6fTITw8HFqtFhKJxNhGIpFAq9VaNI5PP/0U/fr2xZrXXsa4AAcoFQF4OSIA9wx2x3spGzBo0ECWETeDAwcO4JlnnsHcuXPx+OOPw8/Pr9b7np6emDVrFp555hls2rQJqampZum3S5cuiIqKwvINjjjyZ8Pt8guAxS854pahgzFmzBiz9E314yhMNuPq6ootWz6D9qgzxs51xO791VO2NQyFwOr3gEkxDhg5cjSee+45s/X90ksvwcPDAytWrECHDh3qbSOTyZCQkIAvvvgCe/bsMVvfRNQ2+Pn54ZNPNuPL7x0wKcYB+w/WHqMu6oEX1wGKJ0VQKMKxcOFCm8VKdSmVSsjlcsTGxiIzMxMAkJ+fb/Ln8/LyMG3aNOOfG7kYzsjIgEIRjtF9PLDjiUGIm9ID8kG+GD/QF4+F+eOrJwZixnA/xMTE4OOPP27299M/1q1bhy5dumDOnDmNths/fjyCg4Px5ptvmq3vN954Ez179cfoOY54W41aN2wrK4Htu4GRs51wqdAbm9O3crNfE6Wmptb6GczLyzPpcywoQDZ1++234/vv92Lu3PsR+tBx9AtwwkBpBcquirA3W4TyayLMfWAu1q5bB1dXV7P0efHiRWzevBmPPPIIPD0b2WAHwKhRo9CzZ0+sX78ed955p1n6J6K2Y+rUqfj6652InvcQbo86g5v7OaFPjwpcKRVhT5YIIpEjnnj8UaxSKjlr04oYDAbEx8cjPT0dOp0OCoUCSUlJ8PPzMy5Va4q/vz+2b99+wzEIgoDFixZC1tMLL83oDSeHuhe0rk4OiJ/SHYVllXj8sUcxc+ZMs/2ua0+uXr2Kjz/+GFFRUXByavrSdtq0afjPf/6D48ePo1+/fi3u38fHB999/wMeefghRP9nO5a87IA7hlXCyQk4cNQJZ85W4LZbh+GLnano06dPi/trL6KiohAVFWV8beqjCRyJyeaCg4Px++/H8O2332L0+Lm45joZ7h3vwbKE53H69BlsfPttsy1HA4Bvv/0W5eXlCAsLa7KtSCRCWFgYduzYYbb+iahtCQ0NxYk/c/Hll19CdtssXHWZBO+bpmPFCiXy8s7ilVdfNemCiqxn8+bNUCgUAACpVIrs7GxkZGRAJpNBr//nISq9Xg+ZTGaRGL799lscP/Enosd0qTexqSESiRB95024lK/Hli1bLBKLvbt06RJKS0vRt68JG+cBxnanT582Wwy+vr7Y+uk2nDhxAvMXLoWTeBoq3afg7ukxyMzMxE/7M5nYWAlHY2oVRCIRxo4di7Fjx1q8r8LCQgCo9VBpY3x9fVFUVARBEDiVTNROOTo6YsqUKZgyZYqtQyETSKVSYxGBGjKZDHK5HPHx8cZjWq0WkZGRFonh008/RXc/D8h6eTXZNqCjG27p2QFbt27FrFmzLBKPPau5uWBqFbSadpa4KREYGIjExESzfy+ZjskNtTve3t4AgIKCgjoPHNZHr9fD29ubiQ0RURshl8uRnJxs3AneYDBAqVQCqH4WJykpCVKpFACMJaHNTa/Xo4u3s8m/O27ydkJ+/iWLxGLvOnbsiM6dOyM7O9ukioVZWVlwdHRE//4syWyPmNxQuxMaGgpXV1d88803tdZy1qeqqgoZGRmYOnWqlaIjoqbo9XocO3YM5eXl8Pf3N3kpCrUvMTEx9R6Xy+WQy+UW79/T0xNFZabvp1JYWgk/z6ZneaguR0dHxMTEYPXq1Zg3b16jz9NWVVVh+/btmDZtGrp27WrFKMla2sQzNxqNBlqtFhqNBrGxsbYOh9o4Pz8/3Hfffdi2bRuKiooabbtnzx6cOXOGVZCIWoGDBw9i7ty56NatG0aOHIlx48ahX79+uPXWW/HBBx9wU0RqVcaNG4fjZ68g92Ijex38T/6Va8g6eQXjxo+3QmT2KTY2FoIgYOXKlQ1uwC0IAlQqFXJzc/Hkk09aOUKyljaR3ISFhRnXyup0OuM0M5mfIAj4+eef8frrr2PVqlV45513mlU6s6149tlncfXqVSQkJDRYOefnn3/GqlWrMH36dIwcOdK6ARJRLampqRgxYgQyMjIwd+5cvPPOO/jwww/x4osvQhAEzJ07FwqFAuXl5bYOlQgAMGPGDHTq6IcPf7rQZNu0Xy7CydkZDz30kBUis0/du3dHeno6MjMz8cQTT2Dfvn3GGx6CIODQoUN49tlnkZaWhtdffx2jR4+2ccRkKW1iWVpOTo7x73q9HsHBwTaMxn5t3boVy198AQd/PQQ3Fye4OjuiqKQcCxcswH1RUVi5ciW6detm6zDNok+fPti5cyemTJmCWbNmQS6XY8yYMXB3d8fff/+NL7/8EgcPHsTUqVOxadMmPm9DZEO7d+/GnDlzEBoairi4uFoPAffo0QNjxozBvn378PzzzyM2NhbvvvuuDaMlqubq6or/Pv8CFi9ejN5+bphzR+d6f5d8fjAf7/xwHs888wx8fX1tEKn9mDJlCnbt2oUnn3zS+O/p5+eHoqIinD9/HoGBgUhLS0NERIStQyULahPJTc1Df2q1GnK53GJlG9uzl19+GXFxcbgt0AdvzArE7X284SASoaD4Gj47kI/UT9OwS5OB777fg8DAQFuHaxYjRozAb7/9BpVKBZVKVWs/g9GjR+OTTz5BeHg4HB0dbRglET3zzDMYOHBgncTmeiNHjsTixYuxevVqxMfHY8CAAVaOkqiuhQsX4uzZs1ixYgW+/aMQimAJZL06wEEE/J5XjPQsPX7RXca8efPwwgsv2DpcuzBq1ChkZmYiMzMTn332GQoKCuDl5YXQ0FDI5XLuR9UOtInkBqgu16jX6xutblWzm3CNf2/+Q/X7/PPPERcXh4dHdcHC8d1q3Vny9XTGg6NuwpShEsz/UIcpkyfh98NH4OzsbMOIzadr1654/vnn8cwzz+DkyZMoLS1Fp06d+JDhDUpNTa21i7epuwkTNeTgwYPYv38/li9f3mTZ1smTJ+Pdd9/Fhg0bsGbNGitFSNQwkUiEl156CbfffjtWv/Yanvt0d633b7/tVqSueAKRkZFcIWBmISEhCAkJsXUYZANtJrmRyWSQyWRQKBSIjY2FSqWq06aluwm3V4krVyCot3edxOZ6nb1dkDizJ2apjmHbtm3GzdHshbOzMysumcGN7iZM1JCMjAx4eHjg9ttvb7Kti4sLxowZg4yMDCtERmS6qVOnYurUqdDpdDh+/DiqqqoQEBCAgQMH2jo0IrvT6ufm1Gp1rZ3kQ0JCkJWVZcOI7MuhQ4fw0/6fcd8IvybvGvW/yQPDe3tj3dq1VoqOiNq7oqIieHl5mbzZno+PT5NVEIlsRSqVYtKkSZgyZQoTGyILafUzN2KxuNYsQWZmpsV2E26PahLF0f18TGo/pk8HvLM/25IhEbUpJSUlSE1NhVqdjksXL8DDwxNj7hyLmJgY9OjRw9bhtXm+vr4oLCzE1atX4erq2mT7S5cuQSwWWz4wIiJqlVr9zI1cLodEIoFarUZycjKkUini4uJsHZbduHr1KhwdRHByMG2tr6uzCOXl9dePJ2pvNm3aBP9uXREdPQ/64z+he9VpuOT/jtWvKNG7d2/ExsayNHEL3XXXXSgrK8P333/fZNuSkhLs2bMH9957r+UDIyKiVqnVz9wAQHh4uK1DsFvdunVDZZWA0/qr6OXn1mR73cUydOnSyQqREbVuGzduRHR0dHWxjbGD4e/7z6xC8dVKfHYgH2+8vRHnz5/Dli1bWXXvBvXt2xdhYWFIS0vDmDFj4ObW8Di1ZcsWlJaWNrgzPRER2b9WP3NDljVx4kRIxD7YknWpybYl5ZXY8ftlzHngQcsHRtSK5ebmYv78+QgP7oQX7+1VK7EBAE9XR8y6rTNejgjA559/Xm8BFDLdqlWrcPbsWTz77LP1Pk8jCAK2bduGd955B/Hx8VwOSETUjjG5aefc3NzwSHQMth0sQM6F0kbbqr47i7JrlbwrSu2eSqWCh4sDnpjg32ghjtH9fDBugBhvvrEGgiBYMUL7IpPJ8MUXX+D48eOIjIzE6tWr8dNPPyErKwvp6emYO3cuXn/9dTz66KN46aWXbB0uERHZkFWXpb3yyiu1XguCAI1Gg507d1ozDPqXZ599Fl/v+AoLNv2J56f1wO2BHWpdsBWWViD5+7NI/fkiVq9ejZ49e9owWiLbe/ftjZg6RAx356bvD4UHd8SCD47jl19+wa233mqF6OzTuHHjcOTIEeOmu5999hkAwMnJCdOnT8f777+PsWPH2jZIIiKyOasmN4IgQC6XIysrC8HBwQCA/Px8a4ZA9fD29sau3d9i+r334NGPfkJAZ0+M7ecFd2cHnNaXQ3P0MioFYM2aNXjsscdsHS61Q1qtFjKZzKR2YrEYEonEOM6Yu3JWWVkZLlzKx6BRvUxqP7CrBwDg5MmTTG5ayN/fHy+++CL+85//4Pz58ygvL0enTp3g5eVl69CIiKiVsGpys3TpUgBAQUEBhg8fbvw72V6nTp2w94d9+P7777F+3Tpk/LgPZVfL0LlzZzzz3BOYN28ebrrpJluHSe1UUFBQrddisRi5ubl1EpfExESo1WqIxWIkJCRALpebPZaawgCVVaa1r/rfcjRT92mhpjk5OcHf39/WYRARUStkk9+22dnZMBgMAKr3rRk/frwtwqB/EYlEGDt2LJd2UKui1WqRnZ1tnLnRarXQ6/X1zsiEhYUhPT3dovE4OzujT6AUv+QaMG24X5Ptf9ZVPwA/aNAgi8ZFRERENioosHTpUuTk5ECn0yExMdEWIRBRGyGTyWotScvKymp0Rkan00Gr1Vo0ptj5C6A5ehn64sb3fBIEAersfIwaeQd3IyciIrICm62TqFmiRkRkqvj4eCiVykbbGAwGSKVSKBQKpKSk1DvDk5eXh2nTphlfR0VFISoqyuQ4HnroIax4aTme2Xoaq+8LgFsDhQU+/PECsnML8elrHO+IrpeamorU1FTj67y8PBtGQ0T2xCbJTUpKCjQaDdLS0rBx40bMmzfPFmEQURtiMBig0+kabXN9mfKwsDDEx8fXu8eMv78/tm/ffsOx+Pn54dNtn2HK5EmIfv9PPDyqM0b384GTQ3WVwT/OleCjny7gy0N6PPvss7j33ntvuC8ie/TvGwrX32wgImoJmyQ3UqkUy5YtAwAEBATYIgQiamOSk5MhlUobfF+j0UCpVCIjI8N4TK/XWyyesWPH4vs9e/Ho4kVYkpYJSQdXdPF2QfHVKpy+VIxuXW/CW2+9hdjYWIvFQERERLXZJLnRarUQiUQQiUQ4cOAAQkNDbREGEbUhmZmZCAsLq3VMp9MZEx6pVForkcjOzkZkZKRFYwoJCcH+n3/BgQMHoFarcenSJXh4eGD06NGYNm0aK6QRERFZmU1+8y5duhQvv/wy0tLSjDM4RERN+ffMTWxsLOLj4yGXyyGVSqHVapGcnAwACAwMRHh4uFXiGj58uLG8PREREdkOCwoQUZtQX4nn65egAbBaMkNEREStk02SmwMHDmDz5s0Aqpeo7dy50xZhEBERkR2rmcmVSCSQyWSQSqXQaDTQ6XSQSCTQ6XSIi4uzcZREZE422ecmKysLMTExiImJ4Z1WIiIiMrvY2FgEBwcjJiYGmZmZ0Gg0AAClUlnr+kOtVtsyTCIyM5vM3Pj5+RmrpP37AeHWprCwEEeOHEFZWRluuukm9O/fHyKRyNZhERERUQMMBgM0Go2xFHzN/lgajabW3lcymQwqlYo3WonsiFWTmwkTJsDX1xcFBQVITExEQEAADhw4gBMnTlgzDJMcO3YMr7/+Oj788AOUlJQajw8dOgSLFz+KBx98EM7OzjaMkIiIiOqj0WgglUqNszI1y8+0Wi0kEomxnUQigVartVWYRGQBVk1uVCpVnX1tcnNzrRmCSb788kuEh8+Esxtw82gPBNzsBycXEQwXr+HI/lzExsYgbXMaPtv2GTw9PW0dLhEREV1Hr9cjKyvLWHQkKSkJSUlJyM/PN/k78vLyam0u+u+NR4nIslJTU5Gammp8nZeXZ9LnrJrc/DuxOXjwYKOb8tlCVlYWZs6cAf++Lgib3QlOLv88liTu5ILegzzx1wlv7Hjne0TNisJn2z7jMjUiIqJWRCKRIDg42PhaKpUiMTERkZGRMBgMJn2Hv78/tm/fbqEIiagp/76hcP3NhsbYpKDA1q1bAQDDhg0zPuDXWvz3v/9BBz9HTHigc63E5nrd+7pjbKQfPt/+OX7++WcrR0hERESNaejGqUwmg16vN77W6/WQyWTWCouIrMDqyc38+fOxcuVKREZGIjIyEjqdztohNCg3Nxc7dnyNIaO84OjU+GxMn6GeEHd0w9p1a60UHREREZni30lMZmYmIiMjIZfLa113aLVaREZG2iJEIrIQq1dLUyqV0Ol0rXI37927d0MQBPQd7tVkW5GDCIHD3PDNzq+tEBkRERE1R0pKCuLj4xESEgIAxv1slEolkpKSjLM7rJRGZF+sntz4+PhAKpVi2bJlEIlESEhIgLe3t7XDqNeVK1fg7OIIZ1fTJrTcPB1x5UqxhaMiIiKi5pLJZMYlZ9cnMHK5HHK53FZhtRslJSVITU3FW2+tx9FjxyBUCejZqyei58XgwQcfrFW1jsicbPLMTXp6OmJjYxEREWHcPbg1EIvFuFZeibKSSpPaF1+ugI/Yx8JREREREbUdv/zyC3oH9EJ09DzklxzHLePcMVzugUq3vxEXvxQ9enTHZ599ZuswyU7ZZBNPqVRqrJx2/ZpYW5s4cSKcnJ3wR2YRbrlT3GjbykoBfx4owwP332+d4IiIiIhauUOHDiE0dDw6dBJwf3RP+PjV3hOwpKgCe7bkY+bMGfj88y8wefJkG0VK9somMzfZ2dk4ePAgdu/ejQMHDtgihHrddNNNmDF9Bn7fV4yrpY3P3hzdX4grl69iwYIFVoqOiIiIqHWLjY2Bu08V7orpUiexAQCPDk6YMKczuvdzx8MPP4Rr167ZIEqyZzZJbpYuXYqMjAxkZGQgOjraFiE06MUXX0TVNWd8ufECSq/UTXAEQcAfWUX4YZse8+bNw80332yDKImIiIhal4MHD2L//p8RFOYNl0aeX3ZwFOG2Kb44d+489xIis7PJsjSgOsFpjfr374+MbzSYMnUyPnzpDPoM80DAEE84OYtguHgNx34uxoW/SjHngTlYv369rcMlIiIiahU++eQTePm4ImCwZ5NtO/q7omtvD3yc+jFmzpxpheiovbBJcnP58mUkJia2umppNUaMGIEjh49i48aNWL9+Hb765W8AgIODA6ZOnYJFGxdjwoQJEIka3wuHiIiIqL04f/48vCVOcHA07frIu6MDzp09a+GoqL2xSXJTUy3NYDAgOTkZS5YssUUYjercuTOefvppxMfH4/z58ygrK0PHjh1bXSJGRERE1Bq4ubmh4ppgcvtrVwW4+7hbMCJqj1gtrQmOjo7o1q2brcMgIiIiatVuu+02qFRv4fKla/DpWLeYwPXKr1Yh78+riJp2u5Wio/bCJslNdnY2JBIJ9Ho9Dhw4gNDQ0Ebbq9Vq6PV65OTkwGAwQKVSWSlSIiIiIjJFREQEHn/iMRzaexmjp3dstO2xX4pw7WoVYmJirBQdtRc2r5bW1H/UOp0OOp0OMTExUCqV0Ov1SEpKslKkRERERGQKd3d3xC2Nx28/XMYfWUUNtsv7sxT7vyzAnDmz0bNnTytGSO2BTZIbAJDL5YiIiEBycnKj7QwGA9LS0oyvQ0JCkJGRYenwiIiIiKiZli1bhoceehiajy9g5wcXkPdnKYQqAYIg4OJfV/Ht5ov4PPkcxowZi7fe4kocMj+bLEubP38++vTpA0EQkJmZ2WhbmUyG7Oxs4+ucnBzIZDJLh0hEREREzeTg4ICNGzdixIgRePXVl7FtfQ5EDiKIREBVpYCuXW/CC8/HYenSpXBxcbF1uGSHbJLcKBQK43M2hYWFJn/OYDBAo9HUSnaul5eXh2nTphlfR0VFISoqqmXBEpHJUlNTkZqaanydl5dnw2iIiMgWRCIRYmNjERMTgz179uDYsWOorKxEQEAAwsLC4ORks20WqR2wyX9dWq0WycnJkEql0Gq12Llzp0mfi46ORkZGBsRicb3v+/v7c6dbIhv69w2F6282EBFR+yISiXDnnXfizjvvtHUo1I7YLHVetWoVgOpExxRJSUlQKpWQSqXQ6XSQSqWWDI+IiIiIiNoYmxQUCAwMREBAAAICAhAUFNRke7VaDblcbkxoNBqNpUMkIiIiIqI2xibJzcqVKxESEoKIiAiEhYU12lan00GhUCAoKAgikQgikQgGg8E6gRIRERERUZthk2VpKSkpGD58OAAgNze30bZSqRSCIFgjLCJqpbRaLcRiMSQSCbKyshAcHFzvs3cajQY6nQ4SiQQ6nQ5xcXHWD5aIiIhsxiYzN7m5uSgsLMTGjRtx+fJlW4RARG1IYmKicTlrTaJTH6VSiZiYGISHhwOoXtJKRERE7YdNkhtBEJCVlQVBEKDT6WwRAhG1IWFhYRAEAQUFBQ3Oxmg0mlpJj0wmq7UBMBEREdk/myQ3YrEYKpUKkZGRtuieiNognU7XaHVFrVYLiURifC2RSEyuxkhERET2wWLP3GzduhUzZsyo9z2JRIKIiAjk5+cjJyfHUiEQkR0xGAyQSqVQKBRISUmpszQtPz/f5O/ihr9EtsUNf4nIUiyW3HzyySdIS0ursywkJSUFly9fxpIlSwAAS5cutVQIRGQnYmJijH8PCwtDfHw8VCpVrTZ+fn4mV1K0pw1/BUHA5cuXce3aNfj6+nLnb2oTuOEvEVmKxZal5ebmIiIiApGRkSgsLMTWrVsBANHR0fDx8cHGjRst1TUR2RGNRlOnZLxer6/TTiaT1Tqu1+shk8ksHp+tXLx4EatWrUKvXj3g6+uLzp07w9u7A+bNm4cDBw7YOjyiViU2Ntb4d41Gg+TkZKjVaiQlJdkwKiKyBIslNzqdDjNnzoRKpUJcXFytGZzo6Gikp6dbqmsisiNSqbTWhUl2drbxeb3rC5LI5fJar7Vard0+17d371707dcH//nvs/C8qRAT5nTB5Ae7YMgYd6Rv3QSZTIbnn3+eZfSJACQlJdUaG1hVkci+WWz9gq+vL4Dq4gFvvfUW+vTpg6KiInTo0AEAWCWNiEwilUqh1WqRnJwMAAgMDDRelMTGxiI+Ph5yuRxA9UVLUlISpFIpABjb2ZODBw9i8uRJ8O0mwsw5PeDu5Wh8TzoUCJL7QrvbgBdeeAHu7u6Ij4+3YbREtqXT6Wo9n1dfVUWVSmWXYwVRe2Wx5KbmYqNGeHg4Vq5cicTERADgHUUiMllDFx4ZGRm1Xsvl8jpjj71ZsnQJ3H0ETHm4C5xd606+OziKEBzmi/KyKjz3n+fw0EMPoXPnzjaIlMj2NBoN5HK5cbUIqyoS2T+LLUuLj4/HggULjK8TEhKQmJiIhIQEFBYWNrgJHxER1e/48ePYpdmFYWM71JvYXE82XgygCu+8845VYiNqbTQaDSIiImodu5GqijV/rq/uRkSWl5qaWutn0NSqihabuQkICEBMTAwOHjyIYcOGwcfHBwCwbNkyxMXFwc/Pz1JdExHZpW3btsHF1QmBt3g22dbN0xG9B7tDrd6MZcuWWSE6otbFYDBALBbXKjTSXqsqErVFN1pV0aKbeA4fPhzDhg2rdczHxwdxcXHGZIeIiEyj1+vh6e0MJ2fThm4vsRPy9abfqSayF0lJSdDr9VCr1VCr1dDpdFCr1e2uqiJRe2STDRGkUimUSqUtuiYiarO8vLxwtbQSQpUAkYOoyfZXSyvh6ellhciIWpe4uDjj37VaLTIyMozP7l1fZMOeqyoStVcWnblpTEBAgK26JiJqk8aNG4eSK+U4c6K0ybYV16pw6vBVyEPDmmxLZK90Oh1UKpVx5gb4p6pizWtWSiOyL9zKmoiojbjjjjtw882D8ev3J9GjnztEooZnb/7IuoLiovJahV2I2hupVAqVSlXrWHuoqkhkbWfPnsXRo0dRUVGBXr16oX///jaLxWYzN0RE1DwikQgrVqzE6WPF2PtpPoSq+kvqnzpagh+26TFnzmyb/oIhIiL7tnfvXkyffi969OiO0NBQTJw4EQMGDMCtt47ARx99hKqqKqvHZLOZmy1btiArKwsAEBISghkzZtgqFCKiNmPatGnYsGEDFi5ciHO6cgy63RM9B3rA0VGES3+X48j+IuQeLsbkyZORkrLR1uESEZGdeuONN/DEE0+gY1c3jLpXgu793OHgKMLFv67iyE9HMHv2bHz11Vd4//334eRkvZTDJslNSkoKJBKJcVrYYDDglVdewZIlS2wRDhFRmzJ//nzcfPPNePXVV7D9089RteWS8b3Bgwdhw/rH8Mgjj1j1lwnZp4qKCuzduxcXLlyAq6srQkJC4O/vb+uwiMjGNm/ejMcffxzDxvrgjrv8ahW58ZY4I3CoF04cvIJPNqWiY8eOWLNmjdVis1m1tNDQ0FrHdu3aZYtQiIjapFGjRmHUqFHIy8vDkSNHUF5eDn9/f9xyyy2NPotDZIorV67gtddew1tvrcfZs+eNxx0dHTBt2jTExy/DrbfeasMIichWqqqqkJAQj4CbPXHH3X4N/s7pO8wLRfprWLduLeLi4qx2Y8QmyY1Op0Nubi6kUikMBgP0ej0MBkOdhIeIiBrn7+/PO+lkVpcuXYJcHoojRw+jr8wTo+7rDnEnZ1wrr0Lub8XY++NOfD7qc3zwwYe1NtgjovZBo9FApzuJ6Yu7NXkz7eY7fJCdUYjk5GS88MILVonPJgUFoqOjERAQgOzsbOTk5CAgIIBL0oiIiGysqqoK06bdjT91xzDjsa4YF9EJnXu4wsXNAZ7eTrh5pA/C/68rAod5YM4Dc/DDDz/YOmQisrLdu3fDW+KGrgFuTbZ1cXNAz4FuyMj4xgqRVbPZguzQ0FDO1BAREbUiGRkZ+Omn/Zg2vys6dnOtt42jowjjIzuh4PxZvLj8RXyz03oXLURke1euXIGbh6PJS6DdPBxQVFRo4aj+0WpKQW/dutXWIRAREbVr69atRefu7uje173Rdg6OIgwZ5YWMbzLw559/Wik6ImpMRUUFfv75Z+zYsQN79+5FcXGxRfrx9fVFceE1VDWwHcG/FRdWws+vo0ViqY9NZm6Cg4Ph5+cHQRAgEokgCAJyc3NZDpqIiMiGvt/zPQbe0fgGsTUCb/HCrtSL+OGHH9CnTx8rREdE9bl8+TLWrl2LDRvWIS/vrPG4VwdPPPTgw/i///s/9O7d22z93XPPPXjppZdw6kgJAm72bLRtcWEFTh0txROvzTRb/02xSXKTkJCAmTNr/59ktTQiIiLbKi0tg4tb/cvR/s3JWQQHRxFKSkosHBURNSQvLw9yeShydH+i73BP3DrDH15iJ1wtrcSfB6/gnfdU+HDTB9jx1de47bbbzNJncHAwgoJk0O4+hp4DPODo1PDNkGyNAa4uLnjggQfM0rcprLIsrbCw9jq7fyc2APj8DRERkY1JfMUo1F8zqW3x5UpUVQrw8/OzcFREVJ+ysjJMmjQRf184CcX/dcO4yE7oGuCGDr5O6NjNFbdN8UPUsm7w9L2GyVMmITc312x9r1nzBvL/voadH1xAeVlVnferKgX88rUev/1wGYmJqyAWi83Wd1OsktwkJibWOZabm4vg4GBMnDgRBw8etEYYRERE1IiIiPvwp7YUlRVNr6U/8nMhPDzcMWnSJCtERkT/lpaWht9/P4xJD3WCb2eXetu4ujti8iOdUVFVildffdVsfY8cORKfbfsM53WV+GD5GXy/5SJ0vxXj5OFiZH6jx0cr85D5TQFWrFiBxx9/3Gz9msIqyY1EIkFwcDBCQkKMiUxSUhJSUlKwc+dOZGVlWSMMIiIiasSCBQtQXFSOQ3svN9qu+HIFjvxUjNmz58DHx8dK0RHR9daufRO9Bnihk3/jS0ld3R0xYIQn3nv/PRQVFZmt/8mTJ+OPP44jbkkCzh93wY53z+HLt8/h0HclUMyYDa1Wi6efftps/ZnKatXS0tPTkZaWhrS0NADVG3kOHz4cQHXyQ0RERLY1cOBALFmyBD99kY9f9xhQVVl3BqfgQjk+V12Aj5cf/vvf/9ogSiK6cuUKsrKyETis8cqGNfoO90LxlWKzTyh0794dL774Is6ePYcLFy7g77//xuXLhdi4caPxOt/arFJQQCQSISAgAAAQEhICANDr9cb3dTqdNcIgIiKiJiiVSlRWVmL16tU4tOcK+od4QNzRGdfKBZz8vRSnjl1B79698fXXO9GtWzdbh0vULtWUeXbzcDSpvauHQ63PmZujoyM6depkke9uLqskN4IgIDg4GCKRCFKpFDk5OfD19cXWrVshk8lw6dKlJr/DYDBY9WEkIntXUFAArVaL0tJSdOzYESEhIXB0NG2QJCL75eDggNdeew33338/1q9fj08+SUVJSQEAQCYbjv+8/SgiIyPh4eFh40iJ2i8fHx+IRCJcuVxhUvsrhup27eFa2irJzdKlSxEeHg6RSGSss7106VJs2bIFKpUKq1atavCzarUamZmZ0Gg0yM7Otka4RHbt6NGjePnll5Ga+jHKyq4aj/fo4Y+FCxdj8eLF8PLysmGERNQaBAUF4e2338bGjRtRUlICV1dXODnZZAcJIvoXNzc3TJw4AQcz92LIyKafezv2SxE6d+6EW2+91QrR2ZbVnrlRqVSQyWTw8/PDggULAFSXhK6vktr1wsPDERsba40Qiezezp07ERQkw5bPPsYt4z0xa1kPPPh8L8x4tBs8byrEf/77LEaNGomLFy/aOlQiaiVEIhE8PT2Z2BC1MosWLca50yU4eaTxpWaXL13DCW0J5s9fAGdnZytFZztWSW5eeeUVhIWFITc3FzqdDuHh4UhISLBG10T0P4cOHcL06feiS4Aj7ovrhmC5L3w7u8DT2wldA9wRGtUZMx7ripyTx3DX3VNRUWHaVDcRERFZ35QpUzBlymRkfHipwQRHf64cXyRfQHf/nlYvyWwrVrkNExAQUGuTztDQUBgMBmt0TUT/s2LFCrh5ARMf6Awnl/rva3Ts5oqwOR3x6dpMfPnll7jnnnusHCURERGZwsHBAZs3pyMiQoEvN+7ATb080D/EA15iJ1wtrULOwRKcPHIF/fr1w86d37Sb6sRWSW4KCgpMOtZSeXl5mDZtmvF1VFQUoqKizN4PUVtz7tw5bN26BbdNFTeY2NToJnVH114eeHPtm81OblJTU5Gammp8nZeXd0Px/ptarYZer0dOTg4MBgNUKlW97bRaLcRiMSQSCbKyshAcHNwuHp4kIqL2ydPTE9u3f46vvvoKa9etxTfqb4zvDR06BMnJj2HWrFntqgCIVZIbX19fREZGGstAZ2ZmIjIy0uz9+Pv7Y/v27Wb/XiJLEgQBV65cgYuLC1xdG9+I60bt27cPFRWV6DPctEIB0mHu+P6r7yEIAkQikcn9/PuGwvU3G26UTqeDTqdDXFwcAEChUCApKcn4+nqJiYlQq9UQi8VISEiAXC5vcf9EREStmaOjI+6++27cfffdKC4uhsFggKenZ7u9uWeVZ25mzpyJZcuW4dKlS7h06RKWLVuGGTNmWKNrolbr8OHDWLRoEbx9vOHt7Q03NzcMGNAPb7zxBi5fbnx38OZqdj18dwdUXKvAtWvXzBrHjTAYDMbNf4HqvbIyMjLqbRsWFgZBEFBQUFBv8kNERGTPPD094e/v324TG8BKMzcAMHz48Fo7lRYWFsLb27vJz2k0GmRkZECn0yE5ORlyuRxSqdSSoRJZlCAIWLFiBZ577jl4+bii/63u8OvaGZUVAk4dPYcn/+9JvLRiOb76cgeCg4PN0qevry8AoKigAuJOTVdKuWKogIeHO1xcXMzSf0vIZLJaZeBzcnIgk8kabK/T6WAwGBptQ0T2r6HlrBqNBjqdDhKJpNasMBHZB4slN7t37270/fT0dGzYsKHJ75HL5ZDL5VAqleYKjcimXnnlFTz33HMImeCLILkvHJ3+WfY1IMQbVwwV+ObDiwiVj8fP+3/BgAEDWtznuHHj4NXBE0d/KcTtU/0abStUCTieVYrpM2a2uF9zMxgMTe55ZTAYIJVKoVAokJKSUu/dKz6fR2Rblno+r0Zjy1mVSqVx9jcpKQlqtRrh4eFm7Z+IbMdiyU1MTIxxiUh9uCEntUfnz5/H0888jeHjxBgxqf6qJV5iJ0yd1xlb15zDkiVP4Ysvvmxxv15eXnjowYfx9rsqDBnpAy9xwz/6x7KKYLhUhkULF7W4X3OLjo5GRkZGg9PtMTExxr+HhYUhPj6+3uIDfD6PyLYs8Xze9WqWs9YkNzXLWWUyWa3xQyaTQaVSMbkhsiMWS25UKlWt8s//duDAAUt1TdRqbdy4ESKRAFmouNF2ru6OGHpnB3yl3oGTJ0+id+/eLe776aefxtatW/BF8gVMeqhTneVpgiDghPYK9qjzMWfOHNx2220t7tOckpKSoFQqIZVKodPp6ixP1Wg0te7IAoBer7d2mETUCjS0nFWr1dYqhyuRSKDVam0RIhFZiMWSm8YSGwC1nr8hai82p6ch4GZ3kx7s7yvzwr7P9Ni2bRueeOKJFvd900034dtvv8OECWH4WHkaAYM9IR3iAWc3BxTmX8Mfv5Tg0tlSzJo1639JmOlV0ixNrVbXet5Oo9EgJiamVpIjlUoRGxtr/Ex2drZFqjISUdty/XLWxMREkz/H5atEtnWjy1etVlCAiID8/Eu4aYBpP3Yurg7w8HJGfn6+2frv27cvfvvtd3z00UdYu/ZNaD4+DABwdHTAPffcg4ULF2H8+PGtKrHR6XRQKBS1jtU8gxcbG4v4+Hhj4qPVapGcnAwACAwM5FITIqq1nNXPz8/kTcS5fJXItm50+SqTGyIr8vDwQHlZsUlthSoBV8sq4enpadYYvLy8EBsbi9jYWBQVFaGkpARisdhie+y0lFQqbfDZvX+XhGYyQ0TX+/dyVplMhszMTOP7er2elRWJ7IxV9rn5t61btxo39CRqT8aPk+Pk71dRWVn/xfr1Tv9RirKSaxgzZozF4unQoQO6dOnSahMbIqIbVd9yVrlcDp1OZ2yj1Wq5fJXIzlg1udm9ezdCQkIwb948PsBH7dKCBQtQZLiKPw9cabSdIAg4tKcQQ4cOwe23326l6IiI7EPNctagoCCIRCKIRCLjcjSlUmksAQ1wxpfI3lhlWdru3bsRFxcHnU6Hp59+GkuWLIGDg00mjYhs6pZbbsGMGdPxxZbP0cHXCd0C3eu0EaoE7Nuej9N/FOONbctb1fMvRERtQWPLWWv2zyMi+2TR5KYmqcnNzcWyZcuwdOlS43u8YKP26oMPPsRdd9+F7W/tQeAtnhh8RwdIurigslLAqaMlOPxjMS7+VYq1a9finnvusXW4RERERG2GxZKbCRMmIDs7GwkJCViyZImluiFqczw9PbHz651Yt24d1q59A5+uPWl8TyQSYfLkSVj6QRzGjh1rsxiJiIiI2iKLbuJZ8zAfEdXm4uKCJ598Eo8//jh+/vln5OXlwdXVFUOGDDHLhp1EliYIAv744w9cvHgR7u7uGDRoEDw8PGwdFhERtXMWS24CAgKwdOlS5ObmIiUlBSEhIRg2bJiluiNqkxwcHFgwgNqUq1ev4p133sHatW/iyJGjxuPePh3wyMPz8NhjjzFBJyIim7H4U/0BAQGIjo6GIAh4+eWXcfDgQUt3SUREFnD58mWMDx2HRYsX4arTGUyddxNmxfdA+OP+CJQ5QpWyDkNvGYK9e/faOlQiImqnrLaJ5/DhwzF8+HAcOHAAL7/8coNVTIiIqPWpqqrCvdPvwYEDWZixuBtu6u1W6/0uvdwQJBdj53sXMWXqZGT+koUBAwbYKFoiImqvrF6Pefjw4Vi6dCmys7Ot3TUREd2gjIwMfPft95DP7lgnsanh6u6ISQ91hqNLBVYmrrRyhERERDZIbmoMHz7cVl0TEVEzrVu3Fp27e6BH/7p7M13Pxc0Bg273RNonn+DSpUtWio6IiKgad9IkIqJGCYKAnd98gz7D3U3ao6x/UAeUl1/Dnj17rBAdERHRP5jcEBFRoyoqKlB+tRxunqb9ynDzdAQAFBUVWTIsIiKiOpjcEBFRo5ydneHu7oaSwkqT2hcXVgAAxGKxBaMiIiKqi8kNERE16e67p+FEdqlJlS6PZRbBw8MdY8eOtXxgRERE12FyQ0RETVq0aBHyz5fiz1+LG21XeqUSR/cXY86cB+Dj42Ol6IiIiKoxuSEioiaNHj0aM2ZMx7efXELu4foTnOLLFfgi5TzcXDogISHByhESERFZcRNPIiJqu0QiETZt+gj33ReJ7W9/jm4BHug/whM+fs4ov1qFnF+LkfNrCTr6dUSG5hv06tXL1iETEVE7xOSGiIhM4u7ujq1bP8Vnn32GdevXYXfabuN7PXt2x4qXnsUjjzwCPz8/G0ZJRETtGZMbIiIymaOjI2bMmIEZM2ZAr9cjPz8f7u7u6NatGxwcuNKZiIhsi8kNERHdEIlEAolEYuswiIiIjHibjYiIiIiI7AKTGyIiIiIisgtcltZKHT16FB9//DHOnTsHV1dXBAcHIyIiAh4eHrYOjYiIrKikpASffPIJsrOzUV5ejq5du+L+++9H//79bR0atQGCIGDPnj34/PPPYTAY4OXlhbCwMEyaNAmOjo62Do/I7JjctDLHjx9H7PwYfPft9/DwcoGPnzMqKgSsX78OTzz5OP7vyafw7LPP8sFdIiI7V1lZiRdffBGrX38NV4quoGM3Dzg6iXD5UjmWL1+O8aHjoXpLhT59+tg6VGqlMjIy8Pjjj+Ho0WPwlrjBs4MjykqrsGbNGvTs2QOrVikRFRVl6zCJzIrJTSvy22+/Ycydo+HgchVh93dG4C1ecHQSAQAuX7qG3/ZdxvPP/xfHj/+BDz74kAkOEZGdqqqqwv33z8Lm9HTcMsYbQ0b2hLefMwCg4lr1vkLZmh9x620jsHfPDxg0aJCNI6bWRq1W4777ItFV6oZp87uie193iEQiCIKAC6evQvutHrNmzcKFCxfw+OOP2zpcIrNpE8mNRqOBTqeDRCKBTqdDXFycrUMyu/Lyctx111S4eJbj7vk3wc2j9lSxT0dnjLqnI7r0dMPHmz5GSMgIDkbUbpg6BrSHsYLah9deew2bN2/GxAe6IPAWr1rvOTk7oH9wB/Qa6IHtb53HXXdNxfHjJ+Dk1CZ+pdtcexgncnJycP/99yPwFk+ERnWCg6PI+J5IJEKXXm6YNLcLfvw8H08++SSCg4MxcuRIG0ZMZD5t4ta/UqlETEwMwsPDAVTfjbA327Ztw+nTZzA+yq9OYnO9vsO90C+oA1avfhWVlZVWjJDIdkwdA9rDWEH2r6KiAqtXv4aBI7zrJDbXc/N0xLhIP+TmnsTnn39uxQjbtvYwTmzYsAFOLgLGRXSsldhcTyQS4Y67/CC5yQ2vv/66dQMksqBWn9xoNBqIxWLja5lMhrS0NNsFZCEbNqxH9z6e8Ovq2mTbm0d649SpM9BoNFaIjMi2TB0D2stYQfZvx44d+Pvvs7h5pHeTbTt1d0W3AE9seGu9FSJr+9rDOFFeXo6Nb6egX7AHnFwav8wTOYgw6DZPbNv2Kc6fP2+lCIksq9UnN1qtttYmcRKJBFqt1oYRWcbhI4fRNdDFpLZderrCxdUJR44csXBURLZn6hjQXsYKsn9HjhyBu6czOnVv+mYXAHQNdMbhw79bOCr70B7GifPnz+OyoRDd+7qb1L57X3dUVFTizz//tHBkRNbR6hfo5ufnm9w2Ly8P06ZNM76OiopqM1VAqqqq4OBQ/9Txv4lEIjg4iLgsjVqd1NRUpKamGl/n5eW1+DtNHQPay1hB9q+qqgoikWm/D4Dqu++VlVUWjMj8LDFWmKI9jBM11wYiU68pHGp/jqi1uNFxotUnN35+fjAYDCa19ff3x/bt2y0bkIX07tULF8+cMKmt4UI5ykqvoVevXhaOiqh5/v3L//oLgxtl6hjQXsYKsn+9evVCyZVyXM6/Bp//VUhrzKW/ytG7V18rRGY+lhgrTNEexokuXbrA1dUFF85cRc/+Te+Nd+HMVQDgNQW1Ojc6TrT6ZWkymQx6vd74Wq/XQyaT2TAiy3jkkWjkHinGFUNFk21//6kQvhIx7r77bitERmRbpo4B7WWsIPt37733ooO3Fw7/WNhk28L8azh1tBjz5kVbIbK2rz2ME+7u7oiKmoVjP5egqkposv3Rn4oxdtydTG7IbrT65EYul0On0xlfa7VaREZG2jAiy5g9ezY6eHlhz5Z8VFU2PBidO1WGIz9dQWzMfLi5uVkxQiLbaGwMuP54exkryP55eHggJjoWh3+8YryrXp/KSgF7tubDR+zTJpZLtQbtZZxYvHgxLueXIfObgkbbHfm5EHm6Yjz2KLeWIPvR6pMboLpsY1JSkrFcY035RnvSoUMHfPJJGs78UYov3z6PS3m1f6FVlFfhyP5CfJF8HkGyYPznP/+xUaRE1tfQGBAbG1uramB7GCuofXjxxRdxy9Bh+Pytczj6SyEqrtV+pubiX1fxZcp55P15Femb1fD09LRRpG1PexgngoKCsHLlSmR9U4Dv1RfrrAopvVKJn3fo8d3mi4iNjcW9995rm0CJLKDVP3MDVN9pkcvltg7D4iZPnowdO77GnAdmI+3Vv9A1wAM+HR1RWSEg78RVlFy5BoUiHO+88y7c3U2rgkJkDxoaAzIyMkxqR9TWeHh4YNeu3Zj74Fx8+smn2P/FZfj3dYGjkwiGC5U4d6oE3bp1xTc7P8K4ceNsHW6b0l7GiYSEBHTo0AHLlsXjyP7T6N7PAx7eDrhaUoUzf5RBJHLA008/gxdffLFZBSyIWrs2kdy0J3K5HKdPncH27duxadMm/H02D27ebpgWOwKxsbHo27dtPTRKREQ3pkOHDti6ZSuOHz+Ot956C1nZmbh69SoGyrpjzuo5uPvuu+HkxF/j1LDFixdj7ty52LRpE7Zt+xT6Aj38u/ggdvYEPPzww+jUqZOtQyQyO5EgCE0/bdZGDB48GIGBgbYOg4j+JycnB4cPH7Z1GHVwrCBqXVrjWMFxgqh1MXWcsKvkhoiIiIiI2q82UVCAiIiIiIioKUxuiIiIiIjILjC5ISIiIiIiu8DkhoiIiIiI7AKTGyIiIiIisgtMboiIiIiIyC4wuSEiIiIiIrvA5IaIiIiIiOwCkxsiIiIiIrILTG6IiIiIiMguMLkhIiIiIiK7wOSGiIiIiIjsApMbIiIiIiKyC0xuiIiIiIjILrSb5Ear1SIoKAi+vr4IDAxEfHy88b3rj9f8USgUAIDAwECo1WpbhU1EVmYwGKBQKIxjQVhYGLRaLYDGx4oaYWFhCAoKavD7k5KSan3+399nSl9BQUHG2Gr6M2WcUqvVSE5OvtF/GiK6jj2OFQqFAklJSbWOxcbGIiwsrNYxrVZb65hGo6nz/4/IZoR2oKCgQBCLxUJGRobxtUqlMr4vlUqF9PT0ej/b2HtEZF9ycnIEsVhca3zIzs4WCgoKBEFoejzIzs4WZDKZIJPJTOovOztbaGgYbqwvmUxW672CggJBKpUK2dnZJvVLRC1jr2OFSqUS5HJ5ne8Qi8W1jimVSkGpVAqCIAjh4eGCXC4XpFKpSf9fiCytXczcZGVlAQDkcjkAQCwWIyYmpkXfaTAYbug9Imq9YmNjERMTU2t8kMlkEIvFJn0+OjoaCQkJFoquYWKxGHK5HBqNxup9E7VH9jpWRERE1HpPp9NBKpVCLpfXmvFJS0tDeHg4ACA9PR1KpdKygRM1Q7tIbmqSGoVCYZZf/lqttsHp18beI6LWTaPRIDY29oY/K5VKIZVKzRyVaX1v3rzZeLFRn+Tk5Bv+/0ZEtdnrWCEWiyGVSo3XSmq12rikLSMjA0D1DVyDwWCT+IlM0S6SGwDIzc0FUJ3g+Pr61llzGh8fj6CgIOMfnU5nizCJyEZqfuab+oXd0FgRGxtr9ruXjY1L0dHRxrX2SqUSubm5DcauVqsRExPDmR0iM7DnsQIAwsPDkZ6eDgDIyMhAREQE5HI5Nm/eDKA6Qaq5aUzUGjnZOgBrEYvFxh9WtVoNhUKBnJwc4w+4Uqls9K5nzeeio6ONy858fX0hkUiQk5PT6HtE1PrVjAU1yzAaUt9YoVarIZPJIJVKjQ/5mkNj41JKSgrCw8ONd5AbWw4THh4OrVbLO61EZmDPYwUAREZGGleg6PV6iMViiMViSCQS6HQ6ZGRkcIUKtWrtZubmeuHh4ZDJZM0eWMLDw1FQUIDs7GzI5XIUFBQYk5fG3iOitkEul0OlUjX7c5mZmdDpdAgKCoJCoTBWZ7QGuVwOqVRap8LRv6lUKi5LIzITex4rZDIZ9Ho91Gp1rRmamuduOHNDrV27SG40Gg2SkpKMsyoajQY6ne6Gfzhr7mI09z0iat1UKhWSk5NrlUvW6XRN3ghRKpXIzs5GdnY2VCoVZDIZsrOzLR2ukUqlQmJiYqPFTDQaDcLDw1kKmsgM7HmsAKoTmcTExFrlnhUKBdLS0jgDTK1eu0hugoODkZ+fj4CAAPj6+iI+Ph67du264SREKpUal7g15z0iat2kUilyc3ORkZHR4N4UrZFUKkVERESt/bv+reaua0REhBUjI7JP9jxWANVL07RabZ2Zm/qKJsXGxkKhUBhnpLg3INmaSBAEwdZBEBERERERtVS7mLkhIiIiIiL7x+SGiIiIiIjsApMbIiIiIiKyC0xuiIiIiIjILjC5ISIiIiIiu8DkhoiIiIiI7IKTrQMwp8GDByMwMNDWYVAbU1paipMnT6JAr0dlVQVcXNzg7++Pbt26wcGB+X9L5OTk4PDhw7YOow6OFdRcVVVV+Pvvv/H333m4erUMjg5O8JVI0Lt3b7i7u9s6vDavNY4VHCfoRhgMBpw6dQpFRYUAAHd3D/Ts2RMdO3aESCSycXRtm6njhF0lN4GBgdi+fbutw6A2oqysDAsXzMd773+ADl4OmDSyEh08gT9PX8H3mRdw5sxJvPVWCmbOnGnrUNusadOm2TqEenGsoObYsmUL5s+PxqVLBbgzxAF9elahqBj4et9l5OTkYO7cOdiwQQU3Nzdbh9pmtcaxguMENcfZs2dxX6QCe/fuQ4+uTrhrdAWcnICsw0XYv/9v9OsrRdrmLRg2bJitQ22zTB0n7Cq5ITLVtWvXMP3eafjuu1148xkBD06vhKdHzbtVOJ4LPP26AQqFAps2bcKsWbNsGS4R2Uhqairuv/9+zAgDVj4B9AuoMr5XXFKJ9z4Flry8CWf//huff/EVnJ2dbRcsEdnEhQsXMHr07SgrzsOnbwJ3ja1ObABAECrx4wHg8ZWnMGbMSOzZs48JjoVxzQ21S2vWrIFmlwZfbKjCovtxXWJTrV8AsHm1gAfuEfDwww/i7NmztgmUAACxsbHGv2s0GiQnJ0OtViMpKanJ40Q36ty5c3joobmYfXf1eNAvoPb7nh7AovuBLzZUYdfuXXj99ddtEicR2daiRQtw5XIefthUgXvlMCY2ACASASNlwHfvV6Jfz6uIjJiJqqqqhr+MWozJDbU7lZWVWL/uDdx/l4DQ2xtu5+AAvJ4AODpU4u2337ZegFRLUlISdDqd8bVSqURMTAzCw8MBAGq1utHjRDdq48aNcBBV4o1nBDT2+F3o7cD9dwlYv+4NVFZWWi9AIrK5v/76C59+ug3/WVCB3v4Nt/PyBNY8XYnjJ3TQaDTWC7AdYnJD7c7333+P3JNnsOC+ptuKvYGoKVV45+1kywdGdeh0OojFYuNrjUZT67VMJkNaWlqDx4la4t13UjBrahXE3k23XXAfcPLUX/juu+8sHhcRtR6bNm2Cm6sIs014HOSO4cCQ/o545x3eMLUkJjfU7pw6dQoAIBtkWvvgm4HTZ/IgCIIFo6L6aDQayOVy42utVguJRGJ8LZFIoNVqGzxOdKMEQcCp038haLBp7Wva1YwvRGQ7V65cgUqlwuTJk3HrrbdCLpdj1apVuHjxotn7OnXqFPr1doC3V9NtRSIgaFAlTp3UNd2YbhiTG2p3aso7V5mYq1RWgiWhbUCj0SAiIqLWsfz8/HrbNnS8IXl5eZg2bZrxT2pq6g3HSfbLwcEBlSYuja9ZQu/o6Gi5gOxIampqrZ/BvLw8W4dEduKtt95Ct27dsHDhQuTn56Njx464evUqnn/+eXTv3h3Lli0z6zMvzRkngOprCo4TlsVqadTuDB5cfYt1935g8pim2+/aL8Kggf1Zn97KDAYDxGIx9Hq98Zifnx8MBkOdtg0db4i/vz9LvFKjRCIRBg8agN37D2Px/U3fCdm1v/p/Bw0ycUq4nYuKikJUVJTxdWssBU1tz6uvvoolS5bgrrvuwpw5c9ClSxfje4WFhfj000/x8ssv48KFC3j77bfN8nt90KBBSE6uwNkLQNfOjbetrAS+y3TCpLuGtrhfahhvR1O7ExQUBNnwoVj7cdP/+Z85C3y2G4idv8gKkVGNpKQk6PV6qNVqqNVq6HQ6qNVqyGSyWsmOXq+HTCZr8DhRS8TELsRnu4HTfzfddt3HDhg+bAiCg4MtHxi1OYIgQKPRIDIiArcMGYwhgwfhnnumYfv27SxCYSaHDh3CkiVLMGvWLCxZsqRWYgMA3t7emDt3LpYtW4Z3330XmzdvNku/999/P1xcXJGc3nTbL74DzpytqFUBlMyPyQ21GmVlZfjss8+wYcMGvP322zh48KBF+hGJRPi/p+Lw1fdVWP9xw+2uFAOzljrCz88Xs2fPtkgsVL+4uDhj5TO5XA6pVGr8+/WV07RaLSIjIxs8TtQSs2fPhp+fL+6Pc8SV4obbbUgFvvy+Cv/3VBxneKmOP//8E0OH3IywsDBk7/kK/d0uYJDnJfyp/Q733HMP+vYJxIEDB2wdZpu3bt06dOrUCQ8//HCj7SZMmACZTIY333zTLP2KxWLMmxeDVSkO2JvVcLuTecCil5wwZvRIBAUFmaVvqh+TG7K5oqIixMfHw9/fH/feey8effRRzJs3D8OHD8ett96Kbdu2mb3PWbNm4YnHH8ei5UD0c8DhE/+8d+0aoN4JjLzfEQf/cMFnn32BDh06mD0GappOp4NKpTLO3ADVJZ+TkpKMr2tKPzd03NJOnjyJhIQEBEp7Q+zdAd26dsF9kZHYs2cPi1CYWXl5OTZv3oxx48bBx8cHHh4e6NOnD5YvX26Rvag6dOiA7du/xME/XDDyfkeod1aPDzUOnwBi/gMsfBF4/LHHcP/995s9BmrbcnNzMfKO23Hlwikkz+2LtNh+WDa1J+Kn9MCHj/TFB/P6w/2aHneOGY1ff/3V1uG2WeXl5fjoo48wdepUODk1/cTFtGnTsG/fPpw4caLJtqZQKpW4/Y5RmDDPAS+sA85e+Oe9wivA+o+B2+5zhKtHN3ySZsIUD7WITZ65qVlLX0Or1UIsFkMikSArKwvBwcEQi8XQaDTQ6XSQSCTQ6XSIi4uzRbhkQfn5+QgNDcWJEycwdepU3H333ejRowcqKyvx008/YevWrZg+fTpWrlyJhIQEs/UrEonw2urV6B0QgMTE5diozkd/qTO8PIDTZwVczK/AqJEjsPeT9dxJ2IakUilUKlWtY3K5vFYFtaaOW4ogCHj++eexfPlyeLk5YcIgH/gHeuNKWSV2ffsl0jZvxrixY7Fl61b4+vpaLS57deLECUyZMgV//vknbrnlFtx3331wdnZGbm4uVq5cieXLl2PdunWIjo42a7+33XYbfvjhJyxetACKJ35CJz8n9OomwpUS4FjONXTp4ofXX38Ojz32GGdtqI55jzwCl6oSpDwUCImnc533B/t7QvWAFNHv52BW1H34/fAR/nd0Ay5duoTi4mL079/fpPY17U6dOoW+ffu2uH83Nzd89dVOxMfHI2ljMl566yoG9XGCkyPwR24Vyq4KmD79HqxbtwGdOzfxYA61mFWTG7VajczMTGg0GmRnZxuPJyYmQq1WQywWIyEhwXiBolQqkZGRAQDGO7LWuhtLlicIAmbOnIlTp05h7dq1kEqlxvecnJwwevRojBo1Cu+99x6efvppBAYG1qme1RIikQiPP/44FixYgG3btmHfvn0oLS3FpM6dERERgaFD+cAfNey5557DihUrEHtnV8y5ozPcXf6pfrNwvIAfThTiv9t/xMQJYfju+z3w8PCwYbRt25kzZ3DnnXfC2dkZGzduRJ8+fWq9v3DhQiQnJyMmJgYODg545JFHzNr/Lbfcgr0//IjffvsNaWlpuHDhAtzd3fH8HXdg+vTpcHFxMWt/ZB+OHj2K3d9+i5dm9K43sanh4eKIx+VdseCDY9izZw/uvPNOK0ZpH4xVUE2sglbznJM5q5a5ublhzZo1ePHFF/HRRx/hyJEjqKiogGJ2b8yZMwf+/o3s8ElmZdXkJjw8HDKZrM7OrGFhYUhPrz1NV9+mfCqVismNHdm3bx++//57JCYm1kpsricSifDggw/i+PHjWL58ORQKhdnvarm4uCAiIsKsiRPZt0OHDmHFihVYHNoND426qc77IpEIo/v5YN39znjk3YNYvXo1nnnmGRtEah/i4uJQUVGBtWvXws/Pr877Xl5eePLJJ1FVVYVHH30U06dPr7XvkbkMGTIEQ4YMMfv3kn368MMP4evlitCB4ibbhvT2Qq9OHnjvvfeY3NyAjh07omPHjtBqtRg5cmST7bVaLRwcHNCvXz+zx+Lj44OFCxea/XvJdK3mmRudTldr0z1uymf/NmzYgB49euDWW29ttJ1IJMLMmTPx+++/Y9++fVaKjqhhGzZsQCdvN8y+vUuj7QZ29cCkm8V4a8N6VFRUWCk6+3Lu3Dls2bIFERER9SY2NUQiER5++GFcu3YN77//vhUjJKrfmTNnENDRFS5OTV9qiUQi9O3kgjOnuQnsjXBycsK8efOwc+dOlJSUNNpWEARs374dU6dO5WyKnWo1yY3BYIBUKoVCoYDBYGj2pnyA+Tfmq6ysxFdffYX//ve/iI+Px+rVq/HXX3+16DvpHz/++CNGjhxp0gaZMpkMHh4e2L9/vxUioxvVHjbmEwQBH3+0CdNuEcPZselZxJlBHfFX3t/44YcfrBCd/dm6dSsEQcCkSZOabCuRSDBq1Ch8/HEjZRCJrMTJyQkVzdncsUqAszOXON6o+fPno6KiAkqlssGbSYIg4O2338aff/6JJ5980soRkrW0ik08Y2JijH8PCwtDfHw8AgMDm7UpH2C+jfkEQUBKSgpWrngJp06fQUdvN7i7OOLC5TIsXboE0+6+G6+tfh29e/ducV/tWVlZGdzc3Exq6+DgADc3N5SVlVk4KmqJ9rAxX1lZGQqLrqB3x4ZnEa7Xu2P1f+Pnzp2zZFh26/z585BIJCZXLOzevTu+//57C0dF1LQhQ4bgo00foqD4GnwbeeYGAMquVUF7phQxd/NZzxvVq1cvpKamQqFQ4KmnnsLs2bMRFBRkvIF67NgxfPLJJ/juu++QlJSEcePG2ThishSbz9xoNBqEhYXVOlazAZ+tNuWLj49HbGwsBvoU4/15/fH1EwOxbVF/ZDx1M5ZO9MfP33+DW0eE4I8//rBKPPaqS5cuJs+EXb58GZcvX2aVEbK5mofHy66ZVua57Fr1rVtXV1eLxWTPam5qmFpWuzk3TYgsae7cuXBwcMSn2qZXonxzuACXi8tr3eyl5rv33nuxc+dOVFZWYunSpbjvvvuwcOFCzJ49G/Pnz4dOp8MHH3yApUuX2jpUsiCbJzdSqbTWTq3Z2dk23ZTvww8/xMsvv4ynJnbHihm9cbO/p/EBdk9XRyhCOuH9R/rAy6EUUyZPQnl5ucVjspaaZXhT75oKPz9feHp5ondALzz99NM4dcr864CjoqKwd+9eXL58ucm2O3bsgKOjI6ZPn272OIiaw9HREcOH3YI9xwtNav/9H5chEomsdnPG3tx2220oLCw0aQ+Qqqoq/Pjjj7jtttusEBlR4/z8/PDQww9j497zOHDqSoPtTpwvxesZZ3HvPffUqQRIzTd+/HgcOnQIP/zwAx544AHcfvvtmDlzJrZv3w6dToc5c+bYOkSyMKsmNxqNxrghX3JyMnQ6nbFKVnJyMpKTkxEYGGizTfkEQUCSchXu7C/GrNsaniGQeDpj5fSe0OWetMgGk7Zw9uxZjBgRgqlTpyL70HfoM0KE4XJ3dOhmwOrXX4ZUKsWqVavMuinhww8/DEEQ8O677zbaLj8/H1u2bIFCoUCnTp3M1j/RjVq4aDF+OGHAGX3jyyQrqwSos/WYOnUKevXqZaXo7MvYsWPRr18/bNmypcm2P/30E/7++29WKqJWY/Xq1bj9jpFY+FEO1u7Kw9nL/9wQzb9yDW/vOYt57/+J3n3649333rNdoHZGJBJh5MiRePnll7Fx40a8/vrruPvuu81a+plaL6s+c1OzyZ5Sqax1vKGkxdqb8v3000/4/fARrJvd9J2Tvl3cERTgg/Xr1rb5EsIFBQUYO/ZOnLt4GtMXd0PXALda5ZZvv7sK2l0FSEhIgCAIZttMs1OnTnj99dexcOFCODo64pFHHqmzF8iff/6JF198Ec7Ozli1apVZ+iVqqaioKCx/8QUsST+FDbOl9e5hUSUIUH51BifOlyAlLt4GUdoHkUiEZ599Fg888AA++ugjzJo1q95y8DqdDi+//DLGjRvXZAVGImtxc3PD1zu/wTPPPIOUZBXe//EC/CUecBABf+tL4ejkhNlz5uLVV1+Ft7e3rcMlsgutoqBAa/Hrr7/C0UGEEVLTHly9XeqJj7IPWTgqy1u+fDlO/3USMx+7CeLOdSu1uLg64LYp1Q9PP/PMM7jvvvsQEBBglr4XLFiAyspKPPHEE9ixYwdCQ0MREBCA8vJy/PTTT/j111/Rp08f7Nq1C927dzdLn0Qt5enpiR1f78T4cWMxO+UEIkMkmDa8I3w9nHCtUsB3xwxI/eUSfvurGCkpKRg9erStQ27T5syZg5ycHLzwwgv49ddfMX36dIwYMQIODg7466+/sH37dnz55Zfo27cv1Go1d3inVsXV1RWvvPIKnn/+eaSlpeH48eOoqqpCQEAAoqKi4Ovra+sQiewKk5vrVFRUwNHRAQ4m/mJ0chC1+b0rSkpK8PY7GzHwVs96E5vrBcl9ceTHYqhUKrPOoixevBj33HMPUlJS8MEHH+Cbb76Bi4sLZDIZ0tLScO+993IHcGp1Bg0ahJ9/ycRzzz4L1ebNeEPzNzzdnHH1WiUqKqswauQd+Hrj83UKptCNef755zF48GAkJiYiISEBIpEIjo6OqKiogEQiweLFi/Hss8+aXFWNyNq8vLzwyCOP2DoMIrvH5OY63bt3R/m1SpzKL0Mvv6ar7Zy4UAp//25WiMxydu7cicLLRRh8e88m2zq7OKCPzB0fbvrA7EvEevTogRdffBEvvviiWb+XyJJ69eqFDz78EK+tXo0vvvgC+fn58PDwwKhRo7iTvQUoFAqEh4cjKysLBw8eRHl5Obp164bJkyezQhoREQFgclPL5MmT4ecrxpasS/i/iY0vgTKUVEBz5DJeWB5npegs49y5c3BwEMHbz7T/FMSdXHDsl0sWjoqobenYsSMefPBBW4fRLohEIoSEhCAkJMTWoRARUStk81LQrYmbmxvmxcTi0wN6/HGupMF2giBgTUYeIHLEww8/bMUIzc/NzQ1VVQIqK0yrglZxrQqurlwiRkREREStD5Obf3n22WcxcPAQLPhQB82RAlRU1b7oP3e5HM99egrbD+YjZePGNl+aeMSIEQCA3MMNJ3PXO3W4DCNGsBIREREREbU+XJb2L15eXtDs2o37IiMQn56BLmJ3jO7jCTdnB5zKL8e+E5fh6elhLEna1g0ePBijRo/E7z8cQJ+hnhA5NFxM4fzpMvydW4K1ry62YoRERERERKa54eTmlVdeqfVaEARoNBrs3LmzxUHZmlgsxtc7v4FWq8X69evxy/6fUFZQhi433YR1T87B/fffb1cVeZ55+llMnjwZ+z7Px8hpfvWWUS3UX4NmUz4GDRqIu+++2wZREhERERE17oaTG0EQIJfLkZWVheDgYADVO8nbE5lMho0bN9o6DIubNGkS3nzzTTz66KO4lFeBISO90PtmTzg6inDFUIHD+wtx5MdidO7YFTt2fA0nJ074EREREVHrc8NXqUuXLgVQvbv98OHDjX+ntmnx4sWQSqV46aXl+Pr9/RA5iODk5IBr5ZXw9PTA3AcewfPPP4/OnTvbOlQiIiIionq1+BZ8dnY2DAYDACAzMxPjx49v6VeSjUyZMgVTpkzBoUOHsH//fpSWlqJLly6YMmUKvL29bR0eEREREVGjWpzcLF26FC+//DJEIhESExPNERPZ2NChQzF06FBbh0FERERE1CxmeXiiZokaERERERGRrbR4n5uUlBRERkYCQLt4+J6IiIiIiFqnFs/cSKVSLFu2DAAQEBDQ4oCIiIiITKFWq6HX65GTkwODwQCVSgUA0Gg00Ol0kEgk0Ol0iIuLa/Q4EdmPFs/caLVa7Nq1CwcPHsSBAwfMERMRETQaDbRaLTQaDWJjY2sdT05OhlqtRlJSUpPHicg+6XQ66HQ6xMTEQKlUQq/XG3/2lUolYmJiEB4eDqA6CWrsOBHZjxYnN0uXLoUgCEhLS0N0dLQ5YiIiQlhYGGQyGeRyOXQ6HZKTkwHwooWIqhkMBqSlpRlfh4SEICMjAxqNBmKx2HhcJpMhLS2tweNEZF9YUICIWqWcnBzj3/V6PYKDg+u9OFGpVBCLxfUer0l0iMj+yGQyZGdnG1/n5ORAJpNBq9VCIpEYj0skEmi12gaPE5F9aXFyc+DAAWzevBlA9RK1nTt3tjgoIiKpVAqgegZGLpdDJpMhKSmJFy1EVIfBYIBGo0F2dnaD21Lk5+c36zvz8vIwbdo04+uoqChERUW1KE4iMl1qaipSU1ONr/Py8kz6XIuTm6ysLMTExACoXvNORGQuWq0Wer0efn5+ABq+OOFFC1HbcqMXLQ2Jjo5GRkYGxGIx/Pz8jJuLX6+h4w3x9/fH9u3bWxQXEd24f/9uvv73dmNanNz4+fkZq6SFhYW19OuIiIxkMhlkMhkUCgViY2MRGBjIixYiO3CjFy31SUpKglKphFQqhU6ng0wmQ2ZmpvF9vV5vHEvqO05E9uWGCwpMmDABkZGReOuttxASEoKIiAgmN0RkFmq1utZ4EhISgqysLMhkMuj1euPx6y9a6jtORPatZtlqzTJWjUZjLEJSQ6vVIjIyssHjRGRfbnjmRqVS1dnXJjc3t8UBERGJxWIoFArj68zMTOPFSXx8vPH49Rct9R0nIvul0+lqjRNAddXEmv9NSkoyJj01xUUaOk5E9uOGk5t/JzYHDx40DhZERC0hl8uhVquNG/RJpVLjZnu8aCEioLroiCAI9b4nl8shl8tNPk5E9qPFz9xs3boVM2bMwLBhw4x/JyJqqYaSE1602F5ZWRnUajU2b96M8xfOwcPDE2NGj0F0dDS6d+9u6/CIiKgda1FyM3/+fGRlZRk3wQoJCTFLUERE1Dqp1WrEzo+BPr8A3ft4wMvXEfn5VVAm/YAVK17CI4/Mw5tvvgkXFxdbh0pERO1Qi5IbpVIJnU6H4cOHmyseIiJqpTZt2oQ5c+agzy1emBTdA76d/0lgysuqcOTnQrzzzkb8ffZvfLr1Uzg5mWWfaCIiIpPdcLU0APDx8YFUKsWyZcuQkJCAwsJCc8VFREStSF5eHh5+5GEMGNEBE+Z0rpXYAICLmwOG3SnGpIc646svv8T69ettFCkREbVnLUpuACA9PR2xsbGIiIhAcnKyOWIiIqJWJiUlBQ4OAkbf2xEiB1GD7XoN9ETgLV54443XUVVVZcUIiYiIzJDcSKVSBAQEYPjw4VyeRkRkpzZuTEYfmQdc3Jr+tXHzHR2Qk5OLH3/80QqRERER/aPFyU12djYOHjyI3bt348CBA+aIiYiIWpHKykrk5Z1F5+6uJrXv3KO63cmTJy0YFRERUV0tTm6WLl2KjIwMZGRkIDo62hwxERFRKyISiSASidDAliJ1VP2vnaOjo+WCIiIiqodZStksXbq0We0NBgPEYrHxtUajgU6ng0QigU6nM27W19BxIiKyHgcHB/Tv3xd5f57F4Nu9m2yfd6IUADBw4EBLh0ZERFRLi2duLl++bHK1NLVajfj4eISGhtY6rlQqERMTY9y0T61WN3qciIisa/78hcg5VIziwoom2/7+YxGCg4MwbNgwywdGRER0HatWSwsPD0dsbGytYxqNptYsjkwmQ1paWoPHiYjI+ubOnYsOHbyw6+NLqChvuAraob2XcfpYMZ56aokVoyMiIqrW4mVpNdXSAECv1zf781qtFhKJxPhaIpFAq9U2eJyIiKxPLBZj+2efY9Kkidi2/hxkch/0HugBB8fqstD6c+X4dY8BR/YX4f/+7/8QGRlp44iJqDUoKCiATqdDVVUVevTogZtuusnWIZGda3Fyk52dDYlEAr1ejwMHDtRZctaU/Pz8Zh1vTF5eHqZNm2Z8HRUVhaioqGZ/DxHdmNTUVKSmphpf5+Xl2TAaMrcxY8Zgz569WLBgPna8kw0vH1d08HXCtasCLp0tQadOHbFmzUt49NFHIRI1vBcOEdm/X375BWveWIP0zem4du0agOriJJOnTMZjjz6GiRMn2jhCslctTm6WLl2Kl19+GXq9HgkJCc3+vJ+fHwwGg8nHG+Pv74/t27c3OwYiMo9/31C4/mYD2Yfg4GBkZmYhOzsb6enpuHTpEtzd3TFy5EjMmDEDLi4utg6RiGxs/fr1WLx4McQdXREy2Rvd+7gDIuDC6avI3v8dJk36Ck899RRefvll3gghszNLtTS5XA4ASE5OxpIlzVtnLZPJkJmZaXyt1+shk8kaPE5ERLYXFBSEoKAgW4dBRK3M5s2bsWjRIgwd7YNR9/hB5PBP8tLJ3xWDbuuAQ3sv49VXX4Wfn98N3RgnakyLCwrMnz8fu3btgkajwS+//NLsz8vlcuh0OuNrrVaLyMjIBo8TERERUetTWVmJuLglkA7xxKh7ayc2NUQiEW4ZI8awsT5YvvzFZq/SIWpKi2duFAqF8TmbpkpBazQaZGRkQKfTITk5GXK5HFKpFEqlEklJSZBKpQBgLP3c0HEiIiIial127tyJU6fOIPxx/yaXmw0bK8ahvWfwwQcf4LHHHrNShNQetDi50Wq1SE5OhlQqhVarxc6dOxtsK5fLIZfLoVQq6z3eUHsiIiIiat2+/vprSLq4o3NP1ybbeno7oXtfN+zYsYPJDZmVWZ65WbVqFQCwVDMRERFRO1VYWAg3TweTiwS4ezng8mWDZYOidqfFyU1gYKBxnxtWvCAiIiJqn3x8fFBaVAlBEEy6JiwpqoJU6muFyKg9aXFBgZUrVyIkJAQREREICwszR0xERERE1MbcddddKLhYhrO5ZU22vWKowF8nSnH33XdbITJqT1o8c5OSkoLhw4cDAHJzc1scEBERERG1PaGhoQgMDIB21wV07e1Wb7W0GtrdBfBwd8fs2bOtGCG1By2eucnNzUVhYSE2btyIy5cvmyMmIiIiImpjHBwc8Nprr+P0sRJ8l34JlZVCnTZClYCsjAL89kMhVqxYiQ4dOtggUrJnLZ65EQQBWVlZEAQBOp0Ow4YNM0NYRERERNTWTJs2De+99x4eeeRhnPnjKgaM8IB/H3eIHIALp6/i6P4S6C+U4oUXXmCVNLKIFs/ciMViqFQqbrBJRERERHjggQdw8OCvmB31MI7su4pt6//Gp2v/xs87LmNi6L3Yt28f/vOf/7AQFVmESTM3W7duxYwZM+p9TyKRICIiAvn5+cjJyTFrcETUfqnVauj1euTk5MBgMEClUgGo3gxYp9NBIpFAp9MhLi6u0eNERGR9gwcPxvr16/HKK68gLy8PlZWV6NatG7y9vW0dGtk5k2ZuPvnkk3pnZlJSUrBr1y7MnDkTAQEBWLp0qdkDJKL2R6fTQafTISYmBkqlEnq9HklJSQAApVKJmJgYhIeHA6hOgho7TkREtuPh4YG+fftiwIABTGzIKkxKbnJzcxEREYHIyEgUFhZi69atAIDo6Gj4+Phg48aNFg2SiNoXg8GAtLQ04+uQkBBkZGRAo9FALBYbj8tkMqSlpTV4nIiIiNoXk5IbnU6HmTNnQqVSIS4urtZFQ3R0NNLT0y0WIBG1PzKZDNnZ2cbXOTk5kMlk0Gq1kEgkxuMSiQRarbbB40RERNS+mPTMja9v9e6xYrEYb731Fvr06YOioiJj+T6dTme5CImoXTMYDNBoNMjOzkZiYmK9bfLz85v1nXl5eZg2bZrxdVRUFKKioloUJxGZLjU1FampqcbXeXl5NoyGiOyJScmNXC6v9To8PBwrV640XmgIQt065kRE5hAdHY2MjAyIxWL4+fnBYDDUadPQ8Yb4+/tj+/bt5guSiJrl3zcUrr/ZQETUEiYtS4uPj8eCBQuMrxMSEpCYmIiEhAQUFhbWWutORGQuSUlJUCqVkEql0Ol0kMlk0Ov1xvf1ej1kMlmDx4mIiMiy8vPzsXr1ajzyyCOYO3cunnvuOZw4ccJm8ZiU3AQEBCAmJgYHDx4EAPj4+AAAli1bhri4OPj5+VksQCJqn9RqNeRyOaRSKYDqUs9yubzWMlitVovIyMgGjxMREZFllJaWYsGCBejWrSvi4pZgx+407PphC15drUS/fv0wcdJEnDlzxupxmbQsDQCGDx9e55iPjw/i4uKwbNkyswZFRO2bTqeDQqGodUypVBr/NykpyZj01JR+bug4ERERmVdZWRkmTAzDzz/vR5DcB4Nu84a7lyMAoKK8Cn/+Woyfd36PW28dgR9//Am9e/e2WmwmJzcNkUqlxosOIiJzkEqlDT7LJ5fL6zwH2NhxIiIiMq9ly5bh55/34+7YLuga4F7rPScXBwwI6YAe/d3x2brzmDlzBrKysiESiawSm0nL0poSEBBgjq8hIiIiIqIbUFlZCa1WC41Gg59++gklJSUW6aewsBApKcm4Zax3ncTmep7eThg1wxda7QH8+OOPFomlPmZJboiIiIiIyPqKior+V3ynN4KCghAWFoY77rgD3fy74qmnnsLp06fN2l9qairKysow+DbvJtv26OsOSWc3vKV6y6wxNKbFy9KIiIiIiMj6zp49C7k8FMdPHEefYR6Yfnc3eImdcLW0CicOXMFbyWvx7nvv4OsdOzFixAiz9Hns2DFIurjDS9x0GiFyEOGmAGccOXLYLH2bwizJzZYtW5CVlQUACAkJwYwZM8zxtURERGQjubm5UKlU+OjjTbh08RJcXF1wx+13YOHCRZgyZQocHR1tHSJRu3b16lVMmjQRf53VQfF/3SDp4lLr/U7dXSEbL8ZX71zApMkToc0+YJYH+6uqqtCsx2dEIlRVVbW4X1O1eFlaSkoKgH8e5hUEAa+88kqLAyMiIqJqVVVVuHDhAv766y+UlpZatC9BEJCYmIjAwECsefM1+PYsQvAkLwwc6YyDR/Zi2rRpCAqSIS8vz6JxEFHjNm/ejEOHfsOkhzrVSWxquHk6YsrDnVF+rQSvvfaaWfoNCAhAwcWrKCuubLKtIAi49FcFAqV9zNK3KVqc3EilUsycOROhoaEIDQ3FzJkz6y0bTURERM1z4cIFrFy5Er169UCXLl3Qo0cPePt447777sMPP/xgkT5XrVqFp59+GrJQHzzwXHeMmdkJt4wRI1jui5lPdMWMR7sh98wfGDv2TuTn51skBiJq2tq1b6Jnf0906u7aaDs3T0cMuNUD7773Dq5cudLifmfNmgURHHD0l6Im254/dRUX/irBvHnzWtyvqVq8LE2n0yE3NxdSqRQGgwF6vR4GgwGhoaHmiI+IiKhd+vHHH3HX3VNx5UoR+gzzwGR5Fzg6i6A/fw0Z336GtLQ0PP7443jttdfg4GCe+kA6nQ7PPPMMguRi3Dal/g26uwa44+7YLtj6xmksX74cr7/+uln6vlEGgwFisdj4WqPRQKfTQSKRQKfTIS4urtHjRG1RcXExfvklE+Pv62RS+76yDsjWnEFmZibGjRvXor47d+6M+6Lug3rLJ+g1yKPBWaPysir88GkBAgOlmDBhQov6bI4Wj4bR0dEICAhAdnY2cnJyEBAQgCVLlpgjNiIionbp8OHDmDhpAjwk1zDn2R4Yf19nSId6oddATwwfK8Z9cd0wZkZHvPHGGjz99NNm61elUsHN3QlBct9G24k7OWPgrZ545923UVxcbLb+m0OtViM+Pr7OzVSlUomYmBjjRr5qtbrR40RtUc0MjKuHac++uXk41PpcS615fQ169wzEZ+vP44+sIlRW/LM3nSAI+OtECT5bfw4lBgds2bLVbDdgTGGWggI1S9KIiIio5ZYti4eLRyWmPNwVLm51LwpEDiIMGeWDq2VVSEpKwvz5883yoPCmjz5EH5k7nF2avhAZfLs3tLtP4+uvv8bMmTNb3HdzhYeHQyaTQaPRGI9pNJpaszgymQwqlQpisbje4zWJDlFb4+1dXYa55HKFSe2L/9fO17fxGxem8vX1xQ8/7MOcObPx1cc78NPnBtwkdYaDgwj6vyuRf74UgwYNxCdfpWHIkCFm6dNUFkmjtm7daomvJSIisnsnT57El19+haFjOtSb2FzvltE+cHV3gkqlMkvfly7mQ9yp/iUm/9ZB4gQHBxHOnz9vlr7NQavVQiKRGF9LJBJotdoGjxO1Ve7u7gibEIbj2aZt1HnslyJ06tTRbOWggeqfoy+//ApHjx5FzLxFCOh0G7r5yHDXpAh89913+P33w1ZPbAAzzNwEBwfDz88PgiBAJBJBEATk5uayHDQREdEN+PLLL+HgKEI/mVeTbZ1dHSAd6o4tW9VITExscd8urs6oKDetZGtVJVBVJcDNza3F/ZpLQwUObqTwQV5eHqZNm2Z8HRUVhaioqBuOjcjcFi9ajHvuycCpo8XoNdCzwXaF+ddwPLsES5c8BhcX025eNMeAAQPMVonteqmpqUhNTTW+NrVCY4uTm4SEhDrT0bt27Wrp1xIREbVLBQUFcPdwhrOraYsrvHwcocspMEvft912Ow4f/hEyE1aa634r/t9nbjNL3+bg5+cHg8Fg8vHG+Pv7Y/v27eYJjMgCpk6diokTJyBj025MeECEnv096rQpuFCOHW9fhH+37njiiSesH2QL/PuGwvU3GxrT7OSmsLDQuM4PQL3rbPn8DVHTKisr8fXXX+PDDz/EmTOn4eLqiiBZEGJiYtCvXz9bh0dENuLl5YWrZRWorBTg6Nj0TnllJVXw8upglr4XLVyE6dM1OHeqDDf1anhGRqgS8Pu+IowaPRKDBg0yS9/mIJPJkJmZaXyt1+shk8kaPE7Uljk6OkKt3oKZM2fgc1UGukk90T/YA16+TrhaUoU/DxYj93AxAqVSfPNNBvz86q+AaG+a/cxNfdPeubm5CA4OxsSJE3Hw4EFzxEVk1/bt24fAQCnuuusufLtvO/LLf8NfBZnYoHoT/fv3x/QZ01FYWGjrMInIBsaPH49r5ZXI/b3pKmSVlQJyfyvDhLBJZun7rrvuws03D4Zm0yVczr9WbxuhSsC+7fk4m1uKZ595ziz9motcLodOpzO+1mq1iIyMbPA4UVvn5eWFr77aga1bt6J/7xB8u/kiPledxTcfnod7VU+sX7ceBw4cREBAgK1DtZpmz9xIJBIEBwdDJBIhJSUFw4YNQ1JSElJSUjB8+HBs3LgRw4YNs0CoRPZh7969kMvl6NTdCTMf90eXnq4QiarvzlZcq8KJA1fw9fYvMW7cWOzZsxeeng2voyUi+zN06FDcccftOPT9r5De7AmHRmZvTmiLUGS4ioULF5qlbycnJ3z11Q7ceecYbHk9D4Nu98Tg273RwdcZlRUCcn8vxu8/FCFPV4K1a9di4sSJZun3Rmg0GmRkZECn0yE5ORlyuRxSqRRKpRJJSUmQSqUAYKyI1tBxorbO0dER06dPx/Tp1TdGCwoK4OnpCT8/P+P1RXtyQ8/cpKenQxAEY3Kj0+kwfPhwAKhVjYSIaisrK8PM8Bno3MsZd0V3gaNT7UHHydkBA0d4o2M3V3y2/jc8/fTTWLNmjY2iJSJbeemlFQgLk2N32kWMi+xU7/K0M3+UYM8WPSIiFLjlllvM1nePHj3wyy+ZeOGFF/Due+8iW3MaTk4OqKioLjQwctQd2LjuOUyaZJ7Zohsll8shl8uhVCrrPd5QeyJ75u3tXevxkfao2cmNSCQyTm2FhIQAqF67WuP6ad/m0mq1EIvFkEgkyMrKQnBwMMRiMXcVJruhVqtx8cIlzHq4R53E5nqdurvi5lFeePudt7FixQp4eTVdNYmI7Me4cePw4Yeb8MDcB3Dh1N8YeJsneg/ygKOzCPpz5Tiy/wpOHr6CiRMn4r333jd7/x07dsSbb76JxMREfPXVV7hw4QLc3Nxw66232qS0KxGRqZqd3AiCYFyWJpVKkZOTA19fX2zduhUymQyXLl264WASExOhVqshFouRkJBgvMOiVCqRkZEBAEhKSoJareZ0MrVJG99OQY9+nvDt3HQpxuoN8s5gy5YtmDt3rhWiI6LWJCoqqrrE6urXkJaWhh8//6ec8dChQ6BSPYqHHnoITk5m2Y+7Xl5eXoiIiLDY9xMRmVuzR8SlS5ciPDwcIpHIuBvy0qVLsWXLFqhUKqxateqGgwkLC0N6enqtYw3tNszkhtqi3FwdOgU6m9S2g68zvLxdkJuba+GoiKi1Gj58OD784EOsfm01Dh8+jPLycnTt2hWDBw9ul2vpiYiackO3e1QqFZKTkyESiRAREYENGzZg5syZ9ZaFbi6dTgeDwWAs0chdhcmeODg4QBBMby8I1Q8KWkJlZSXy8vJQWlqKjh07tpsSkURtUceOHXHnnXfaOgwiolav2aWgX3nlFYSFhSE3Nxc6nQ7h4eFISEgwW0AGgwFSqRQKhQIGg6FZuwrX7CZc8+f6XU2JWoNBg27GWV25SW3158px5fJVDBgwwKwx5OfnY9WqVQgI6IVevXphwIAB6NixI8aNH4stW7agsrLyhr87NTW11s+gqbsJExEREZlDs2duAgICam3SGRoa2uxdfxsSExNj/HtYWBji4+MRGBho8vdzN2Fq7ebHzse0aV/h/KkydGlkgzwA+H3fZXTs5Id77rnHbP3/+uuvmDRpIvL1lxB4iwemTrgJLm4OuJx/Dcd+yUR4eDim3jUV6ZvT4e7u3uzvv9HdhImIiIjModkzNwUFBSYday6NRoOwsLBax2p2EL6+Ght3Faa2bMqUKejTJxDfbdajrLjhGZKTR4px+KciPPbo43Bxabr4gClOnz4NeVgoqpwLcX9Cd4RGdUbvwZ7oFuiOgSO8MX1xV0yddxMyMnbi/vtnQWjO+jkiIiKiVqDZyY2vry8iIyPxyiuv4JVXXkFkZKRZ9raRSqWIjY01vs7OzuauwmR3HB0d8fnnX0C45o6tb57DH1lFqKz4J4koKqjA/q/y8fW7FzD1rrvMuuRz5cqVKCsvwl3RXeDpU/+kbe9BnhgX6YdPP92GPXv2mK1vIiIiImto9rK0mTNnQiqVIi0tDQCwbNky4waeLSGVSqHVapGcnAwACAwM5K7CZJcGDBiAn/f/guiYaGg+/hY/bjfAp6MzqiqBC3kl8HB3x1NPLcGKFSvMVuL18uXL+PDDDzB4tCfcvRovUBB4iyf8urhj3bp1Nn+A2WAw1KqW2NCeV9wLi4iIiIAbrJY2fPjw/2/vvuOaOts+gP8AQVFEBBHBKrVW3FoQR1G0YlFRtK7WR1HrAq1V68JR7eur1oei1oWjIu6Bo46CuF7kwVWGonUgLoYKCjJkhpXkev/wSQqOKnKSE5Lr+5ck+eQ6gZyf133Ofe5TbkCTm5sryN1Q3zZo4bsKM23TtGlThJ0Lw71797B3716kpKTAyMgI9vb2GDlyJGrXri1ovfDwcEgkhWjR0fKdr9XT00OzDsYIDhbv+rXff/8dV65cQWhoKGJiYpSPv+2eV3wvLMYYY4wB7zm4CQsL+8fnDx8+jM2bNwuyQYzpkubNm2PZsmUqr5OTkwMAqGX6fsczapoaoKioGCUlJYJd81MRw4YNg4ODA0JDQ5WPve2eV2ZmZnwvLMYYY4wBeM/BjZeXF1xdXd96gXHZI6uMMc2jOBNUmC+Didm7d/uiAjkMjQxhaPh+NxxVh7fd84rvhcUYY4wxhfca3GzZsqXc8s+vun79umAbxBgTnrOzM4yMDHH/Wh4cXOr+42uJCA+vF6JP794adQf0t93zqiL3wlJQ3BNL4dUlrBljqhUYGFjuXnR8TyzGmFDea3DzTwMbAIIsKMAYU5169eph+L/+haATh9DGqQ6Marx9ocQn9wrxPFmC77dOVeMWvpuFhcUb73n1tsf/Cd8TizFx8T2xGGOqUuGloBljVdPCHxdCLjXEqe1pKC588z12UpOKELo3Az1deqJ3795q3sJ/9rZ7XvG9sBhjjDGmwIMbxnRE8+bNcTLkFHLTDbDfJwWXgzKQ+qgIWaklSIotwKkdaTi64SkcHDri2NFj0NfXrHh42z2v+F5YjDHGGFOodPdy9OhRdOzYUYhtYYypWLdu3XDzxi18N2k6Ev8iHFmXgsAVTxCyLRXG8sbYtHETzoWGoU6dOqJuZ2hoKLZs2YKEhAT4+/srBy+Ke179/vvvAPDavbBefZwxxhhjuuWD7xAYFhaGefPmIT4+XrnMLGNM89na2mLVqlVYtmwZ4uLiIJFIYGlpCTs7O41ZQEBxbytfX983Pv621zPGGGNMt1X4zE1YWBgcHR0xbNgwDB8+HFlZWW9dIpoxprmMjY3h4OCAbt26oXnz5hozsGGMMcYY+1DvPbhRDGq+/vpr5aBmzpw5AMBNEWOMMcYYY0x07zUtrXfv3oiJicGCBQuUAxrGGGOMMcYY0yTvfRPP33//nee0M8YYY4wxxjTWe01La9KkCby9vVGnTh1s3boVf/31l4o3izHGGGOMMcYqpkILCjRp0gSenp4gIqxcuZIHOYwxxhhjjDGN8UFLQdvb28Pe3h7Xr1/HypUrebU0xhhjjDHGmOg++D43wN+DHL4WhzHGGGOMMSa2Ct/n5k3s7e2FeBvGGGOMMcYY+2CCDG4YY4wxxhhjTGw8uGGMMcYYY4xpBR7cMMYYY4wxxrQCD24YY4wxxhhjWoEHN4wxxhhjjDGtwIMbxhhjjDHGmFbgwQ1jjDHGGGNMK/DghjHGGGOMMaYVeHDDGGOMMcYY0wo8uGGMMcYYY4xpBR7cMMYYY4wxxrQCD24YY4wxxhhjWqGa2BvAXieXyxEWFoZ9+/bh2bNnqFGjBhwdHTFhwgRYW1uLvXmMMQ3x+PFjBAQE4MaNGygpKcFHH32EMWPGoFu3btDT0xN78xhjGqC4uBhHjhxBUFAQXrx4gdq1a8PV1RUeHh4wMTERe/MYExyfudEwERERaNHCDq6urgg+fRAPnl7C9Xv/h2U//y8aN24ET09PFBUVib2ZjDER5efnY+TIEWjSpAlWrPLBjQehuJd8EUeC9qJ79+5o27YN/vrrL7E3kzEmst27d6PhRzbw8PDAhahgPHx2GZE3TuG7Kd+hgXUDrFy5EkQk9mYyJig+c6NBLly4gN69XWFhUw2Dp9rAukkN5dHX4kIZ4qLzsGv3DiQmJeJkyEkYGRmJvMWMMXWTSCT40rUX/vrrGpwHm8POsTaMqr88TkVyQvLDQkSGJMLZuRvCw8+jQ4cOIm8xY0wMGzduxNSpU2HnYII+ExqhrtXfPUPei1L8FZ6DuXPnIiMjA76+viJuKWPCqhKDm9DQUCQkJMDc3BwJCQmYO3eu2JskOIlEgsFDBqG+rSH6TaiPaoblT6pVNzbAZz3MUM+mOkK2/gc+Pj5YvHixSFvLmGbShaxYuHAhrl+PwcDvGsCqcY1yz+np66GRXU1Y2dZA8G9pGDx4EBISElGtWpWIesbUQhdy4tatW5g2bRrad6+Drl9ZvDZNtXZdQzgPrgdT82pYsWIFvvjiC7i5uYm0tYwJq0pMS/P19YWXlxeGDRsGAPj9999F3iLhHThwAC+yXuCLry1eG9iU9VEzY7TobIJNmzaipKREjVvImObT9qzIz8/H1oCtaOtc+7WBTVlG1fXRfZg5njxJRnBwsBq3kDHNp+05Abw8a2NSxwhOA14f2JTVvocZGjSuifXr16lx6xhTLY0f3ISGhsLMzEz5s4ODAw4ePCjeBqmIv/8WNG5RC6YWhu98bRsnUzx/no6TJ0+qYcsYqxp0ISsOHz4MiUSC1p+bvvO1lg2rw6ZJTfhv9VfDljFWNehCThQVFWHPnt1o0akm9A3evbBIK6daOHPmLJKTk9WwdYypnsYPbq5duwZzc3Plz+bm5rh27ZqIW6QaD+Mfwsr2/a6hsbCujurG1RAfH6/irWKs6tCFrIiPj4epWXXUrvvugyAAYNnYEA8f3lfxVjFWdehCTqSlpUEiKUSDj99+dresBrY1QERITExU8ZYxph4aPxE7MzPzvV+bkpKCgQMHKn8eMWIERowYoYrNEtzL08YVWLGEwEu9Mo0TGBiIwMBA5c8pKSlqq60LWaGnpwdUcLfnnGCaSKys0JmcqABF58FZwTTNh+aExg9uLCwskJ2d/V6vbdiwIYKCglS7QSrS3K45niTeeK/XpqcUo7hIimbNmql4qxirmFf/8y/bGKiaLmSFnZ0dcrOKkJNZijrvMYU1LbEUHdq1UMOWMVYxYmWFLuSElZUVapnUwtOEQjRuUfOdr3+WUAR9fX188sknatg6xt7fh+aExk9Lc3BwQFZWlvLnrKwsODg4iLhFqjFp0mQ8vleA7OfvXiTg9uUcWFs34JVNGCtDF7Ji6NChqG1aG7F/5r7ztWmPi5D6WILJkyarYcsYqxp0ISeqV6+OcWPH4W60BDLpP88IISLcichH//79YGNjo6YtZEy1NH5w8+WXXyIhIUH587Vr1zB8+HC1b0fZ02Kq8PXXX8OqQX3852AmSorl5Z67fy1P+e+kOwW4G52P6dN/UPnyrqr+zJpYmz9z1aULWVGzZk18N3kKbl3Kw9P4wnLPlc2JIokM5w9noWnTT9RyEEQXv7v8masmXcgJAPj+++9RVCDFhaMZIHn5AU7ZrIgJzcbz5EL88MMMlW4PoJvfXf7M4tD4wQ3wctnGFStWKJdrVCzfqE6q/mPVqFEDQX8EI+c58MemVCTdKYD8v4H04Ho+JHlSRJ/Jwqkdz9Hf3R1z5sxR6fYAvFPqSm1NCCKh6EJWLFmyBN26dcOJrWm4FvYCRQUyAC9zQi4jJNzMx/ENaSiVGOGPP4JgYGCg0u0BdPO7y5+56tKFnGjRogUCArYhLioPp3Y+x/MnRcrnHlzPx4u0Epw78BxRp7KwdOlS9OrVS6XbA+jmd5c/szg0/pob4OWRli+//FLszVC5Tp064eLFS5gwYTxCAv5CHYsaqGNhgOfJxdi97AkMqxnih+k/wNfXl2/Kx9gb6EJWVK9eHadOnsaMGTOwfcd2XDmTgwa2NZDxrBh7l6cgL7sYXbp0xs6du9C8eXOxN5cxjaMLOQEAY8eOhampKWbMmI7Da1JgaWMM49r6SH1UjP2+T2BRzxybN2/G5Mk8dZVplypx5kaX2NvbIybmGqKiovCthxc6t+8Hi7oN8Ouq1Xj69BlWr14NQ8P3WwaWMaadatSogd9++w0pySn493If9Og8CHVN62PiuCm4fv06IiIieWDDGMOQIUOQmPgIwcHB+KqfBzq26QsrSxvs27cPKclPeWDDtJIeEVVg/WHN1rp1azRt2lQl752SkoKGDRuq5L01tTZ/Zt2orcq68fHxiI2NVcl7V4Y2ZgV/d3WjtrZ+Zk3MCm3MCTFr82fWjdqakBNaNbhhjDHGGGOM6S6elsYYY4wxxhjTCjy4YYwxxhhjjGkFHtwwxhhjjDHGtAIPbqoguVz+7hcxxnQeZwVj7F04J5i24QUFKiEhIQGpqalwcnJSS72CggLUqlVL+bNcLoe+vvDj0/z8fAQHByMxMRGOjo5wcXHh++roKIlEAiMjI/77V4K6cwLgrGDqx1lRedxTMG2nrpzgwU0ldO3aFatXr0bnzp1BRCguLoZUKoWJiYlK6rm7u4OI4OHhgX/9618qCSEA8PLywv3792FtbY2aNWuiTZs2mDlzpkpqfahHjx4hPT0deXl56NChA0xNTdVaX1X/CbzJhQsXEB4ejubNm2Po0KFq+U9BJpNh1qxZKCkpQVFREZYuXYpGjRqpvK42UndOAJwVCmLnBKC+rBAjJwDOCiFxTyEesbOCewph8bS0D5ScnAwzMzN07twZmZmZmDZtGvr27Yvhw4dj2rRpyMnJEbRecXExSktLkZKSgmPHjsHJyQk///wzZDIZli5dColEIsip5YyMDFy9ehXh4eEICAjAxIkTsWnTJty+fVv5mtzc3ErXqYzc3FyMHz8eS5YsQUhICFxdXTF+/Hg8evRIpXWLiopw/PhxAFCGkFwuhyqPD+Tl5WHatGn466+/4O/vj/nz5wN4+XcqLCxUWd0jR46gtLQU48ePR5cuXbBmzRqkpqbi+PHjSExMVFldbaPunAA4K8rWFiMnAPVnhVg5AXBWCIV7CvFwT6GFPQWxDzJ06FBauHAhERH5+PjQhAkT6OLFi7R9+3YaMmQIzZkzR/CaEomEJk+eTMeOHaPTp0/Tt99+S3Z2dlS9enWKj48XpEZMTAxNnjy53GMnTpygxYsXExGRVColR0dHkkqlgtT7EAsWLKBdu3YREVFqairFxsaSv78/9e7dm7y9vSk/P18ldW/dukVNmjShNm3a0OLFiykzM7Pc83K5XPCamzdvpv379xMR0YULF6hVq1bk6+tLvXr1olatWtGRI0cEr0lE1KNHD5JIJEREVFhYSL169aJJkybR999/T926daMTJ06opK62ESMniDgriMTLCSL1Z4VYOUHEWSEU7im4p+CeQjh85uYDtW7dGuvWrUPz5s2xceNGTJ48Gd26dcO3336L5cuXIy4uDvfv3xesXmlpKYyNjeHo6Iht27ahTZs22LlzJ5o0aYJOnTqhT58+2Lx5c6XrtGzZEtbW1vjzzz9BRCgpKUH//v2Rnp4OAPj111/h5uYGAwODStf6UCYmJsqjHFZWVmjVqhU8PT2xefNmJCcn47ffflNJXQsLC9jZ2WHu3Lm4ceMGOnXqBC8vL9y8eRMAcOvWLeW/hWJoaIi0tDQAgLOzMyQSCerWrYvFixdj+PDh2L17NwoLCwU/0uPo6IihQ4ciOjoaFy5cQFhYGNauXYsVK1Zg+PDhOHToECQSiUqPMGkDdecEwFmhIFZOAOrPCrFyAuCsEAr3FNxTcE8hYE4IOlTSQdu3b6dPPvmE+vbtW+5xBwcHun//vmB1ZDKZ8t9r166lr776isLDw6ljx46UlpZGly5dorS0NEFq5eTkUEJCQrm6y5Yto9jYWGrZsiVlZ2cLUudDhYeHk5WVFc2bN49u3LhR7rlbt25Rr169qLi4WPC6sbGxtGHDBiJ6+Ts6f/48TZgwgT799FMaOnQoNWnShLKysgStmZKSQsHBwURElJeXRzNnziz3fMeOHen69euC1iQiKigooAULFlCXLl1o/vz55OXlpXzu+fPn1KFDB8FrajN15QQRZ4WCWDlBpP6sECsniDgrhMY9hfpxT/E3bekpeEEBgaSmpqJBgwa4ceMG9u/fj9u3byMkJAREBD09vQ9+X7lcjkOHDuHatWuQSqVwd3eHi4sLFi9ejJMnT+Lrr7/G3LlzIZPJKn3kQ3FBW3FxMapXr17uuaSkJDg5OeHzzz/HkSNHKlWnMnJzc2Fqaoo7d+5g8+bNKCkpQaNGjeDo6Ii+ffvi0qVLWLVqlXIeq9B1X1VUVISsrCxMnToVRIRjx44JVlMikSAhIQFSqRSfffbZa89fvXoVnp6euH79umA1X61fs2ZNAMA333wDAwMDfP7557h69Srq1KkDPz8/tV4EqQ1UlRMAZ0VZYuVE2dqvUlVWiJ0Tim3grBAW9xTqwT3F37SqpxB0qKQDrl69ShMnTqQ///yz3OOK+aJ79uyhnTt3UlxcHBFVfs7k0aNHydXVlYYNG0bz58+nsWPHUklJCT169IiWLl2qnKNZ9ijMh0hKSqLJkydTx44dlUcSXn3PFStW0LNnzypVpzJu375NNjY2FBYWRkRE165do4CAAJoxYwZ98803VLt2bfrmm28oMTFR8LoNGzZU1i0pKXntNaNHj6aHDx8KVvP+/fvUo0cPatu2LfXr14+WLVtW7vlTp06Rl5cXHThwgIgq//dXOH36NN25c+e1x588eUILFy6kHj160N69e6mgoECQetpK3TlBxFmhIFZOKGqrMyvEygkizgqhcE/BPQX3FMLnBA9uKujy5cvUokULsrOzI2dnZ9q3b1+55x8/fixovWbNmin/nZKSQl9//TX5+PgoH3v27BmlpqZWuo6DgwMFBwfTnj17qGXLlnTv3j26evUq/fzzz7R161YqLS2lvLy8SteprNWrV9O0adPKBfzTp08pNzeX0tLSVHbh3+rVq2n69OlvvOhRLpfT2bNnBa3XpUsXOnv2LEVFRdGZM2eoR48eFB4eTkQvg/DMmTOCnyaXyWT07bffkrW1NfXt25dCQkLKPS+RSFQ2jUfbqDsniDgryhIrJxS11ZUVYuQEEWeFkLinEBf3FNrZU/A54grKyspC165dERERga+++grr169Hu3bt4OPjAwDYvn07IiIiBKlVUFCAzz//HHFxcSAi2NjY4ODBg4iIiEBWVhYAYNy4cbhw4UKl6jx8+BB2dnZwd3fHqFGjMHr0aPj7+yuXhTxy5AgeP36s0vtyvK/Ro0cjPz8fdnZ2+M9//gMAsLa2Rq1atVC/fv1yNyQTum5eXh5atmyJ8+fPA3i5bjsA6OnpwdXVVbBaubm5qFu3LlxdXdGpUyf07t0bM2bMUH6vDA0NERMTI/ipY319fVhaWmLSpEno2bMnFixYgM6dO2PLli0AAGNjY2zfvh0FBQWC1tVG6swJgLPiVWLlhKK2OrJCrJwAOCuExD2FuLin0NKeQiVDJi0ll8spIyNDeYpSKpXS8+fPac+ePTRw4ECqX78+WVpaClrzjz/+oJ07d1JJSYny1KWvry/duXOH4uLiyNzcvNI1Jk6cSJcuXVL+PG/ePBo1ahQRvRzRT5o0iXbs2FHpOkLauXMnzZw5U+0XIr6prtDLNUokEpoyZQrt2bOHiouLSS6XU3FxMQ0dOpSIiHJzc8nc3JyeP38uaF25XE4bNmxQXsCYlJREmzZtom7dupG9vT25uLjQgAEDBK2pjcTICSLOijcRKyfeVlvIrBArJ4g4K4TCPYXm4J5Cu3oKHtx8gFe/eIq5iX369KFFixYREQm2ZrtMJqN79+6VeywiIoIOHDhAY8eOpdmzZ1eqXlZWFhkbG1NgYCDduXOHiouLycnJie7evat8zaxZsygoKOjDP4RAys4BTU9PJ09PT7Kzs6PTp08TkWrWhBer7vPnz+nPP/8s996+vr506dIlWrZsGQ0cOPC1bROCYh16BalUSunp6XT27FnS09NT+e9am6gzJxTvz1khXk6IUVusnCDirBAS9xTi4J5Ce3sKHtwIpKSkhIYMGSL4sn1v069fP9LT06OMjIxKvU9ubi7NmzePunTpQm5ubuTh4UEuLi7lXuPs7EylpaWVqlMZDx48oJSUFCKi15bCPHnyJM2dO5eKioq0pq6C4gI7xU4vk8lo5syZVKtWLbp27Vq55ypr6dKlNGrUKBo9ejQNHjz4tXnf0dHRr30vWMWpOyeIdCcrxNxfxaytzpwg4qxQF+4pVId7Cu3vKXgp6Ep48OABmjVrplwyMS0tDVZWVoIs6/ommZmZsLCwAABcu3YNkZGRmDJlygcvnSeXyyGRSGBiYoKCggLs2LEDQUFBSEtLw6BBgzBgwABERESgVq1aGD9+vNAf571IpVI4OTkhJiYGbdq0gbGxMWJjY/H555+jsLAQRUVFiImJgY+PD+bNm1fl677L9OnTcf/+fZw+fVqw79nly5fh5+eHcePGwcTEBOfPn8epU6dARPDy8sLIkSPx+PFjZGZmomPHjrykawWpOycA3csKMfdXTcwKVeQEwFmhatxTqB73FOVpa0/Bg5sPFB8fjw0bNmDNmjVqqXf//n1s27YNvr6+kMvl0NPTg0wmQ7Vq1T74PQ8fPoynT5/C09NTufY4ABw9ehS7d+9GSkoKnjx5gvj4eJVegPtPJBIJIiMj0bFjR5w/fx5mZmaQy+WIioqCjY0N0tLSkJmZiR9//FHQbRSrroIiZLKzsxESEgIPDw8AQHp6OkpLS2FjYyNYEE2ePBljxoyBk5MTiAhyuRxFRUUIDg5GcHAwvLy80KNHj0rX0UXqzglAN7NCzP1VzNrqzAmAs0KVuKdQD+4pdKSnUNk5IS308OFD5Rrw27ZtUy7Vp4r5zEREycnJytOIAQEBynpCzL2VyWTUokULevDggfKxV9/3ypUrtH///krXYhVT9nT9gQMHKCAggIhUN//36NGj1KRJEzp37ly5x/Py8ujf//43ubu7V3qqgi5Rd04QcVboInXnBBFnhdC4p2DqoIs9xYcP0XVMeno6bt68CWtra0ilUly5cgWjRo0CAJWcdk9NTcXVq1fRuHFjyGQy3Lx5U3kat7J3DQaAsLAw9OrVC59++qnyFLjifYkIcXFxaN++PRwdHStdqzKozJEExb/LnsIU4i7KmlQ3Ozsbt2/fxscff4w6derg1q1bmDt37mvbJKTBgwcjKysLe/fuxaVLl9CrVy98/vnnMDExwYIFC9CuXbvX7i7N3kzdOQFwVii2Q4z9VazaYuQEwFkhJO4pxME9hW70FDy4eU+Wlpbo1KkT8vLysH37dpiYmCA5ORlGRkaoX78+jIyMBK/XokULSCQSHDx4EDVr1kRqaioMDQ1Rr169Sr9/dHQ0OnToAADlvtyKne3+/fs4deoUZs+eXelalVF22xT/Lhv8qmpYxKqbn5+PvLw8PHnyBIGBgSgtLYWpqelr9YUgk8lw8eJFxMXFoaCgAK1atUJ6ejq2bNmCLVu2oGbNmpBKpWjWrBlMTEx4/vx7UHdOKGrqelaItb+KVVudOQFwVqgC9xTi4J5CN3oKvubmA8yaNQuDBg0CEaGkpATGxsZo27Yt6tSpI3gtmUyGH374Ae7u7jAyMoKenh5MTU3RsmXLcnNaK+rEiROYMWMGrl69CjMzM+VRg5KSEhgZGeHrr7/GvHnzRD/KUpYqj0hqUl0iQmpqKhYsWIBu3bqhXbt2MDY2hpWVFerXry9YnS1btmDHjh0wMTHBJ598go8//hheXl64fPkyiouLERoaCnd3d3Tv3h1mZmai/f6rKnXmBMBZoSDm91SdtdWVEwBnhapxTyEO7im0uKdQ6aQ3LaKYO3rmzBnaunUrEb1cVi8uLo7CwsIEvV8F0d9zJM+ePUvbtm0jopfrlF+5coXOnDlTqTm5MpmMCgsLqV+/frRy5crXnn/48CG1a9fug9+ffTjF9+j27dv0yy+/UF5eHkVGRtLp06fpzJkzVFxcLEid0tJSsrGxISKinJwcioqKotatW9Ply5dfey3fp+L9qTsniDgrdJG6coKIs0JVuKdg6qCrPQUPbipo0aJFr90AS5XmzJlT7gI9IhLsy7h3714yNTWlVq1a0caNG+nIkSM0cuRIGjFihPJCQ01w4cIF2rNnD508eVJ5MaQ21yUiWrJkCV25ckX5c1ZWFsXHxwv2/o8ePaIxY8aUeyw6OprmzJmj/E/uyy+/VK59zypG3TlBxFkh5v4qVm1V5wQRZ4WqcU+hftxTaH9PwdPS3oJeDvygr6+PGzduoFGjRjA3N1euOw8IfwGYYv5hREQEbG1tYWNjg9jYWLRu3VqwGq8qKSmBj48PAgMDYWtrCwsLC6xZs0b5GcWi+N3+/PPPePz4MQoKCpCamopz586VW5tfG+rSKxcX5ufnY+/evZg8ebLK5qSWlpZi+fLl6NChA/r16weZTAYjIyNMmTIFmzZtwunTpzFlyhQkJCQIXlubiJETAGeFglg5IUZtMXIC4KwQCvcU3FNwT6HGnFDLEKqKkkqlJJVKadasWcql81S5nKvC999/T97e3oJOLXgfOTk5aq33Prp160bp6em0ceNG5an7X375hUJCQrSu7tq1ayk7O5uI1HPaNicnh2JjY8s9tnLlSoqJiaG+ffvSunXriEg93/mqTKycIOKsUBArJ8Sore6cIOKsEAr3FOLjnkJ1NCkneHDzBseOHaO7d+8qf7506RL179+f8vLy6NatW3TkyBGaMGECZWVlCVpXMTcyLi6ORowYQc7OzhQUFERyuVwlc/Wrgl9//ZXc3d2pR48eysd69uxJL1680Iq6ir/r8ePH6bPPPiOil6d3d+7cSdu2baOSkhJB673q1cCTSCTk6OhINWrUUGldbSBWThBxVrxKrJxQV22xc4KIs6IyuKfQHNxT6EZPwUtBv8HBgwdx+PBhdOrUCT/99BPc3NxgbW2Ntm3b4osvvkBOTg7GjBmDunXrCrrag4GBAXJyctCiRQsEBAQgPDwc9+7dQ1JSEpo0aSJIjark4MGDcHNzw4kTJ1C/fn0sXboUJSUlaNasGczMzKp8XfrvFAXg5Sojc+bMwY0bN+Dt7Y2srCx0794dDx8+RMuWLQWr+apXv7vGxsbw8PBQ3qValfcHqerEygmAs6IssXJCXbU1IScAzorK4J5CM3BPoTs9BV9z8xYymQzLli1TzhWtV68e2rRpA39/fwCqWcovMjISTk5OWLBgAczMzHDu3DkkJSUhISEBt2/fhp2dnaD1NJHi95qdnY0xY8YgKCgI4eHhiIyMRGlpKaysrDBq1KhKLVmpCXVzc3OVa83L5XIsX74cqamp2Lt3L5YvX44pU6Zg4MCBGDduHIYOHSpIzbeJiIhATEwMpk6dCuDluvjGxsbcqLwHMXIC4KwQKyfUXVuTcgLgrKgM7inEwT2FjvYUaj1PVAXIZLLXTtvt37+fOnToQCYmJrRkyRKSSCQqmb8okUho3rx51Lt3b4qKiqKbN2/S//7v/9LGjRsFr6WpFL9Xf39/cnd3p4MHDxKR6udoqrtu27Zt6cCBA8qfpVIpbd26lXbt2kVERPHx8dSwYUOV1JbL5eVOha9atYqioqKUz7F3EzMniDgrxMoJddcWMyeIOCuEwD2FuLin0M2eggc3r1B88e7fv0979+4t91xkZCR17tyZAgMDBaun+MPL5XIqKSmh0tJS2rhxIw0cOFD55VDModSl/0x8fX2pRo0a1LlzZ7p48aLW1c3LyyOZTEZZWVlkZWVF/v7+yucU3wE/Pz8iEj4MHz58SCdPnqT4+HiKjIykhQsXCvr+ukDdOUHEWfEmYuWEumqLmRNEnBVC4J5CM3BPoVs9BQ9uXqHY6T08POiHH34gopdfBFVffOfv70+LFi0iLy8vOn78OO3cuZOGDBmi9ntliEnxO758+TL98ccf1LdvX9LT0yM9PT1q3LgxbdiwQSvqKv5DUdxUbffu3dSlSxcyMjKiqVOnUm5ubrnXCS0jI4MiIiLo8uXLNGnSJFq9ejXl5eWRRCJRST1tJFZOEHFWiJUT6q4tdk4QcVYIgXsK8XBPobs9hWoWxq/CDAwMkJGRgdDQUKxduxbAy7mTAODn54fExESV1DUzM4ONjQ0GDRqEAwcOICgoCHFxcQgMDFRJPU1kYGAAqVQKb29vxMbGYvny5Vi3bh1cXFwwfvx4NG/eXCvqKuZVKy6wc3JywsmTJ5GYmIgXL16gTp06+O2331RyrQYAWFhYoEuXLvjoo49QWlqKli1bIioqCjdv3sSDBw8gk8lUUlebiJUTAGeFWDmh7tpi5wTAWSEE7inEwz2F7vYUvFpaGVTmgj5nZ2fk5OSgTp06MDAwQGZmJjZv3oyxY8cKVk+xakR0dDT09fVx4sQJNGrUCH5+fqhduzYAwMjISLB6mkzxu9+5cyeaNWuGBQsWAAAcHBwQFxeHTp064csvv9SKuoq/e3BwMAIDA1FaWoqzZ8+iZ8+eOHDgAFauXKn8HpIKLjJV1L9y5Qr69euHvn374t69e0hOToZEIkGzZs0Eradt1J0TAGeFglg5IUZtsXOi7DZwVnwY7inEwz2FbvcUfOYGQEFBAR49eqT8g9etWxe1a9fG8OHDcfjwYdy7dw9r165F+/btUbt2bcjlckHqGhgYgIgwe/Zs3L17F/Pnz4e+vj569+6NtLQ0VK9eXaVH5TSJ4nN+8sknMDAwgEQiUT7XsmVLxMfHa01dxYohPj4+GD16NAIDA5GTkwMLCwusXbsW1tbWqFevXrntE7q+TCZDWFgYevbsCQBo3rw5XFxc4OjoKHg9bSFWTgCcFQpi5YQYtcXOCcU2cFZUHPcU4uOeQsd7CnXPg9NEW7ZsodGjRxMRUWZmpvLx//mf/yEPDw+ysrKi+fPn0+PHj4lImHmLivfYvn07jRo1qtxzkydPpm3btglWqyrJzs6mUaNG0YoVK+jEiRMUHh5OrVu3Lvd30Ya69+/fJzc3t3KPxcTEkJubG6WlpamkZll5eXkUERFBRHxX8fclRk6UfR/Oir+JlRPqri12ThBxVnwI7ik0B/cUutlT8LQ0AF5eXpgwYQIAoEePHigoKMBvv/2GJUuWKF9TWloKQ0NDAMKMfOVyOfT09NC4cWMYGxujpKREebq4TZs2iIqKwvjx43XmKAsAxMXFoWXLlpg+fTq2bduGu3fv4tmzZ1izZg3Mzc21qm6DBg1gbGwMX19feHp6wtzcHAYGBnj+/Dnq168vaC3672loqVSKBw8eQF9fH/r6+ujcuTMAKG/6xf6ZGDkBcFa8SqycEKO2OnMC4KwQCvcUmoF7Ch3uKUQeXGkExRr0CQkJ9OTJE1qxYgVZWFjQRx99RGvXrhW01tWrV8uN3DMyMmjYsGG0YsUKOnPmDEVGRlKbNm3o5MmTRKT9R1kUI/zIyEhq27YtZWRkENHL9flVeWRFjLqv/i2jo6PJy8uLli1bRm5ubtS/f3/l8o1CHvlQrNwyf/58mjVrFtnZ2dH3339PRKRcRYW9mzpzgoizoiyxckKM2mLlBBFnhVC4pxAP9xTcUxDxUtDleHt706FDh5Q/Hz16lGxtbcnBwUGwGgMHDiQ9PT0aOHAg3bhxg4iIkpKSyNPTk0aOHEm9evWiffv2CVZP0yl2uP79+9Pq1auJ6O/lDG/evEnJyclaVZeIaM2aNZSUlESlpaV0/Phx2rRpEwUEBND9+/dVVpOIyNramoiIpk+fTsePHyciIh8fH/rzzz9VWlfbqCMniDgryhJzfxWrtlg5QcRZIRTuKdSPewruKYh4cKMc+SYnJ5Ovry95enq+9pqsrCwiEm7km56eTpMmTSJDQ0Oyt7en8PDwf9w2baX4fRYXF9PAgQMpLy+PiP7+3P/617/o1q1bWlP38ePHFBISQra2tmo/ulFQUEBTpkyhZcuW0YABA5SPOzo6qqVZqurEyAkizgoi8fZXsWqLmRNEnBWVxT2FeLinUB9NzwkNmBgnLsX807i4OKxevRoBAQHw9/cv95q6desCqPw8wsGDB+PRo0eoV68efvvtN5SUlGDIkCEYNGgQmjRpglWrVr1x27SV4vdpZGQEe3t7dO3aFREREcq5nNevX0fr1q21oq5UKsVff/2FsWPHombNmso5uADw9OlT+Pn5Ke99oAo1a9ZE9+7dERAQgIYNG+LUqVP49ddfYWNjg2bNmqm0tjZQZ04AnBVliZUTYtQWOycAzorK4p5CPNxTcE+hJObISlOkp6cr7xi7b98+atq0KZmbm9OUKVNIKpUKcrSjtLSUFi1aRMXFxZSWlkZbt24t9/yBAweobt26tGLFikrX0nTp6enUqFEjKiwsLPf4nDlzaOzYsdSsWTNyd3dXnuas6nXL+uabb2jRokU0YMAAGjlyJO3bt4/69u1LQ4cOJSJhj6wp5sUmJSVRTk4OyWQy8vPzo7lz59KAAQPol19+oYSEBMHrait15AQRZ4WCmPur2Fmhzpwg4qwQGvcU6sU9BfcUr9LpwY3iVGJiYiJNnTq13HPh4eHk7OysklN90dHRZGRkRBYWFvTjjz8q52USif+FUJfo6GgiItq4cSN98803lJmZSXK5nCIiIig8PFxl00zEqktEdOjQIeUSncXFxbR3716aOHEifffdd5SYmEhEwv79Fe/l4uJCGzduVD6umBLB3o9YOUHEWSHm/ipWbXXnRNn346yoHO4pxMM9BfcUZen04Ibo5R/L3d2dLC0tKTIy8o2vUdWa3efOnaMuXbqQsbEx9enThzIyMsqFki64desW9ezZkwwNDWnYsGF0+/Ztra179epVmjt3rnIlHYWCggLBaylCKDo6mlq0aEFEfx91yc3NpYsXL+rcd60yxMwJIs4KsXJCjNrqzAkizgqhcU8hLu4puKcg0vFrbui/a3UPHDgQ9evXh5ubG3x9fV97ndBrdkulUgCAi4sLIiIicPPmTZiYmCA3NxfVqmn/rYdkMhkAoLi4GG3atEFYWBhSU1NhaWmJjh07omHDhsrfkTbUBYDY2Fg4Oztj5cqVOHv2bLnnatasKXg9xdzqR48ewc3NDcDfdzG+ceMGfvrpJ534rglBrJwAdDsrxNxfxaqt7pwAOCuExD2FOLin4J7iNSIPrjTKH3/8QR06dCA9PT1at26dyusVFRUREb02X1NX9OnTR3lKt6wTJ05oVV3FEY/s7GxasGABGRsbU7t27VS2PGfZ09AZGRlkZWVFvr6+lJSURERE48ePp8WLFxORZtxJuKpRd04Q6XZWiJUT6q6t7pwoW5OIs0IVuKdQL+4puKdQ0MnBjeJUWmpqKu3fv5+8vb3LrUF+/fp15SnFys5XjIqKouPHj9OzZ89e+6Mr3nvQoEEasXSeOshkMjp37hx999131LZtW+Xjit+FqnYMMeoq3jszM5Pu3btX7rm1a9eSnp5euXsgCF3Xx8eHzp8/TxcuXKAJEyaQu7s7tW7dmtzd3SknJ0fwutpGnTlBxFlRllg5IUZtsXKibG3OisrhnkI83FO8xD1FeTo5uFH8obp27Uqenp7UvXt3MjIyouHDhyvnyAp1Edb8+fOpcePGNHz4cAoMDKT4+HgqLi5WPh8SEkLt27cXpFZV8eTJE7K1taVPP/2Ufv75Zzp//jwVFhZScHAwTZo0Sevq+vv704ABA8jHx4eioqJUVudVhw4doiNHjhAR0dmzZykiIoIuXrxI+fn5RKQ5R1g0lTpzgoiz4lVi7a9i1RYrJ4g4KyqLewpxcU+hHlUpJ/SIxF6MWhw3b96Eh4cHbt26BQAYOHAgiAghISGIiIhA586dBatVUFCAXbt2YfXq1WjSpAkGDRoEZ2dntGvXDm5ubhg5ciRGjx4NuVyuknn7mmjVqlXo168fDhw4gIiICDRt2hRRUVHw8/NDt27dtKrupUuXcP36dcTFxSE/Px+2trZwdnaGq6sriEglf/NDhw5h5syZMDc3V37HWcWpMycAzopXiZUTYtQWIycAzgqhcE8hLu4puKcoR7xxlfplZ2cr/33u3Dny9vYmIqKgoCD65ZdfiIhoxIgRKqsvlUrpxIkT5ODgQE5OTjRu3Dj6+OOPVVavqsjKyqI9e/YojwhoQ903HcG4efMmbdq0iezt7cnLy6vc0TYhyeVykkgktHr1arKxsaGPP/6Ydu/erZJa2kjsnCDirHgTsXJClbXFzAkizorKEjsrOCfejHsKYVXFnNCZMzcpKSnYt28fpkyZAhMTEwAvVxipVq0aTp48idDQUCQmJsLJyQne3t7KVU9U5caNG5gxYwb69u2LefPmaf0RFsXvUyqVws/PD/r6+jAxMYGrqysaN26sdXUVtR0dHTFjxgyMHj1a+fjSpUvRtGlTeHh4CPo9U3yHFLt0XFwcWrVqhVOnTmHWrFm4d+8ekpKSVP65qzJNywlAt7JC7P1VrIxSZ04AnBVC0LSs0KWcALin4J7in+nM4GbMmDFwcXHBmDFjoK+vX27Hz8/Px+LFi3Hnzh0cPnxYGVTqoI7mSBMoPufKlSuxfv16uLm5QSaTQSaToXXr1ujduzfat2+vNXUV/8mtWLEC27dvR2pqKmbPno0ZM2Zg8ODBCAgIwMcffyxoTcVnXbhwIR4/foyUlBT8+eefWLVqFaZOnYqLFy/C2dlZ0JraRlNzAtCNrBBrfxWrthg5AXBWCEFTs0IXcgLgnoJ7indQ9akhTZCVlUUdO3Z86/NPnjxRLmlHpFkXRWmbX3/9lRISEoiI6Pz587R27VoaMWIEBQUFaWXd3NxcKikpoVu3blHXrl3JysqKZs6cSUSquXPw1atXycnJiR49ekRERAkJCdSjRw/lSj2KVX3Y6zgnNIdY+6tYtdWVE2Xfj7Piw3FWaA7uKbineBOdGNwcP36cFixYQET02p1ciYgePHhAy5YtU+ncZkYUHBxMVlZWtGbNGuVjUqmUoqKi3vh3qYp1FTt6ZGQkeXp60uDBg8nU1JTGjRtHUqmUnj17Ri9evCAi4YJILpcr32v9+vX073//m4j+vtfB4sWLld9/9nacE5pBrJxQZ20xckLxXpwVlcdZoRm4p+Ce4m20d0JmGfb29jh69CiePXsGQ0PDcneVBYCQkBBER0fDyMhIzM3Ues7Ozvjqq6+wcOFCfPbZZzhx4gQMDAzQqVMnGBoaakVdxbSEJUuWwNHREQEBAcjJyUFxcTF8fHzQoEED1KlTBwAEmzqgp6enfK+uXbsiKCgI4eHhyueTk5NRu3ZtAFDOnWWv45zQDGLlhDpri5ETivfirKg8zgrNwD0F9xRvJfLgSuXkcjkVFhZSv379aNWqVW98jaOjI4WFhSlfz1Tv559/pvr165Oenh7l5eVpVd3k5GTq3r17ucdiYmKof//+lJKSImit48eP0/nz58s9tm7dOpo1axbNnj2bhg4dSk5OThp3gy1NwzmhmcTKCXXUVmdOEHFWCIWzQjNxTyEMbckJnVlQ4ODBg5gwYQI6d+6MRYsWwdraGv/5z3/w8OFD3L17FyEhIWJvolZSXGQZGxuL/fv346OPPsJ3332nfP7s2bPo3bt3la+bk5OjPHpSWFiIiRMnwtHRERMnTkTt2rWRkZGBL774Ardv3xasZlZWFvr164elS5eid+/eyov/Xrx4gW3btqFevXowMDBA3759YWlpqfWr5wiBc0IcYuWEumuLkRMAZ4UqcFaIg3sK7inei6hDKzXLzs6mMWPGkImJCbVv357atWtH+/fvp4yMDCLii/6EpjhiFR8fT7a2tjRnzhyytbUlGxsbmjdvHqWnp2tF3eTkZPL19S131Oby5cvk6elJa9asoTFjxtDgwYPpp59+IiLhvmdz586lTZs2KX++efMmzZgxgzp27PjWI4rs3Tgn1EusnFB3bbFygoizQlU4K9SLewruKd6Xzpy5edWdO3fQqlUrsTdDa8nlchARDAwM4O3tDZlMhtWrV+Pu3buYPn06nj59irp16+LixYtVvu6rS4ICQGlpKY4ePYqEhAQYGBigZ8+ecHBwgIGBgSA1ZTIZOnfujP/7v/9D3bp18ezZM3Tp0gUuLi5o0aIFdu3aBTc3N/z666+C1NNVnBOqJVZOiFFbjJwAOCvUhbNCtbin4J6iIqqJvQFi4RBSrbI7pJmZGezs7AAAGzZswL59+xATE6OSiy3VXffFixe4e/cudu/eXe5xQ0NDDB8+HI8fP4aZmRlMTU0FqwkABgYGsLe3x4wZMzBmzBj4+fnB0tISO3bsAAAMGTIEP/74I/Lz89V+PxZtwjmhWmLlhLpri5UTAGeFunBWqBb3FNxTVISGTpZjVZm/vz/2798P4OUOuXDhQri5uaG0tBTm5uYICAiAt7c3unfvXuXrXrhwAV9++SWAl+GnIJfLAQBFRUVYv349SkpKBKupMHPmTDx48ACurq6ws7ODv7+/8rnY2FjcvXu3SoQQ001i5YQYtcXMCYCzglVt3FNwT1Fh4s2IY9qopKSE9PT0SF9fn9q0aUPe3t6UlpamfD4sLIyGDh1KO3bs0Iq6jx49oubNm9PTp0+J6O816YuKiojo5SojAwYMELRmWRKJ5LX5ts+ePSN7e3vauHEjEfG8b6Z5xNpfxaotdk4QcVawqol7Cu4pPoTOXnPDVOfYsWM4deoUnJ2dcf36dQQHB6Nr166YOXMm2rdvrzV1iQjFxcUYOnQoXFxcMHv27Nde07FjR6xYsQI9e/ZUrjyiCopVS54+fYrNmzcjOzsbfn5+KqnFmBDEygl119aknAA4K1jVwz3F37ineD88uGEqsWXLFhw9ehQLFiyApaUltmzZgtDQUFhaWmL9+vUqCyQx6mrSkqBFRUV48eIFTE1NUatWLZU3SoxVhlg5IUZtTcoJgLOCVS3cU3BPUSGinC9iOiEwMJCmTp1Kjx8/JiKip0+f0qJFi8qd2tWWurwkKGMfRqycEKM25wRjH457Cs6K98VnbphK7dq1Czt27MD8+fPRt29fra8L8JKgjFWUmPurWLU5JxirOO4p2PvQ2aWgmbAyMzORnp6OR48ewcTEBJ06dUJSUhI8PDxgYWGBzZs3w97eHlZWVlpR959wCDH2ZmLur5qWFZwTjL0d9xR/46yoOD5zwwTh6uqKS5cuYdCgQbhz5w5sbGzQrFkznDlzBj169EBoaChcXFwQEBCgFXUZYxUn5v7KWcFY1cE9BasMHtywSiMiHD16FJs2bYKFhQXGjx8PIsJnn32G6tWr48aNG2jbti2MjIwEvfGUWHUZYxUn5v7KWcFY1cE9BassHtwwwcjlcly4cAHr1q2DiYkJunTpgj59+sDW1haGhoZaV5cxVnFi7q+cFYxVHdxTsA/FgxumEnfv3sW6deuQkpKCNm3aYNq0abC2ttbauoyxihNzf+WsYKzq4J6CVQQPbphKPX/+HH5+fli4cCFq1Kih9XUZYxUn5v7KWcFY1cE9BXsfPLhhjDHGGGOMaQV9sTeAMcYYY4wxxoTAgxvGGGOMMcaYVuDBDWOMMcYYY0wr8OCGMcYYY4wxphV4cMMYY4wxxhjTCjy4YYwxxhhjjGkFHtwwxhhjjDHGtAIPbhhjjDHGGGNagQc3jDHGGGOMMa3AgxvGGGOMMcaYVuDBDWOMMcYYY0wr8OCGMcYYY4wxphV4cMMYY4wxxhjTCv8P8mDqsqCWpRUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 830x456.5 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.lines as mlines\n",
    "\n",
    "plot_smoothed = True\n",
    "\n",
    "\n",
    "def make_colours(y, n, k):\n",
    "    sorted_indices = np.argsort(y[n, k])[::-1]  # Sort in descending order\n",
    "    # Initialize colors to gray\n",
    "    colors = ['olivedrab'] * len(y[n, k])\n",
    "\n",
    "    # Assign specific colors to top 3 values\n",
    "    if len(sorted_indices) >= 1:\n",
    "        colors[sorted_indices[0]] = '#FFD700'  # Highest value gets gold\n",
    "    if len(sorted_indices) >= 2:\n",
    "        colors[sorted_indices[1]] = '#C0C0C0'  # Second highest gets silver\n",
    "    if len(sorted_indices) >= 3:\n",
    "        colors[sorted_indices[2]] = \"#CD7F32\"  # Third highest gets bronze\n",
    "    return colors\n",
    "\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    figwidth = 8.3\n",
    "    fig, axs = plt.subplots(2, 3, figsize=(figwidth, 0.55 * figwidth))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    x = np.arange(len(sims))\n",
    "    y = y_lnZ\n",
    "\n",
    "    for n in range(len(catalogues)):\n",
    "        i, j = n // 3, n % 3\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        ax.text(0.025, 1.075, catalogue_to_pretty(catalogues[n]),\n",
    "                transform=ax.transAxes,\n",
    "                verticalalignment='center', horizontalalignment='left',\n",
    "                bbox=dict(facecolor='white', alpha=0.85, edgecolor='none', pad=3),\n",
    "                zorder=5)\n",
    "\n",
    "        k = 0 if not plot_smoothed else 1\n",
    "        colors = make_colours(y, n, k)\n",
    "        ax.scatter(x, y[n, k], c=colors, s=75, zorder=-1, marker=\"o\", edgecolors=\"k\")\n",
    "\n",
    "    for i in range(3):\n",
    "        axs[1, i].set_xticks(\n",
    "            np.arange(len(sims)),\n",
    "            [simname_to_pretty(sim) for sim in sims], rotation=66,)\n",
    "        axs[0, i].set_xticks([], [])\n",
    "\n",
    "    for i in range(2):\n",
    "        for j in range(3):\n",
    "            axs[i, j].set_xlim(-0.75, len(sims) - 0.25)\n",
    "\n",
    "            axs[i, j].tick_params(axis='x', which='major', top=False)\n",
    "            axs[i, j].tick_params(axis='x', which='minor', top=False, length=0)\n",
    "            axs[i, j].tick_params(axis='y', which='minor', length=0)\n",
    "\n",
    "        if plot_smoothed:\n",
    "            axs[i, 0].set_ylabel(r\"$\\Delta \\log_{10} \\mathcal{Z}_{\\rm smoothed}$\")\n",
    "        else:\n",
    "            axs[i, 0].set_ylabel(r\"$\\log_{10} \\mathcal{Z}$\")\n",
    "\n",
    "    kwargs = {\"marker\": \"o\", \"linestyle\": \"None\", \"markersize\": 10, \"markeredgecolor\": \"k\"}\n",
    "    gold_circle = mlines.Line2D([], [], color='#FFD700', **kwargs, label=r'1\\textsuperscript{st} place')\n",
    "    silver_circle = mlines.Line2D([], [], color='#C0C0C0', **kwargs, label=r'2\\textsuperscript{nd} place')\n",
    "    bronze_circle = mlines.Line2D([], [], color='#CD7F32', **kwargs, label=r'3\\textsuperscript{rd} place')\n",
    "    olive_circle = mlines.Line2D([], [], color='olivedrab', **kwargs, label=r'Remaining')\n",
    "\n",
    "    fig.legend(handles=[gold_circle, silver_circle, bronze_circle, olive_circle], \n",
    "               loc='upper center', bbox_to_anchor=(0.5, 1.05), ncol=4, frameon=False)\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fname = \"../../plots/lnZ_comparison.pdf\"\n",
    "    if plot_smoothed:\n",
    "        fname = fname.replace(\".pdf\", \"_smoothed.pdf\")\n",
    "    fig.savefig(fname, dpi=500, bbox_inches='tight')\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evidence plot for Harry's talk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"Carrick2015\", \"Lilow2024\",  \"manticore_2MPP_N128_DES_V1\", \"CF4\", \"CLONES\"]\n",
    "catalogues = [\"CF4_TFR_w1\"]\n",
    "\n",
    "y_BIC = np.full((len(catalogues), 2, len(sims)), np.nan)\n",
    "y_lnZ = np.full_like(y_BIC, np.nan)\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    for j in [0, 1]:\n",
    "        for k, simname in enumerate(sims):\n",
    "            fname = paths.flow_validation(\n",
    "                fdir, simname, catalogue, inference_method=\"mike\",\n",
    "                sample_alpha=simname != \"IndranilVoid_exp\",\n",
    "                zcmb_max=zcmb_max, smooth=j)\n",
    "\n",
    "            # y_BIC[i, j] = get_gof(\"BIC\", fname)z\n",
    "            y_lnZ[i, j, k] = - get_gof(\"neg_lnZ_harmonic\", fname) / np.log(10)\n",
    "\n",
    "        y_lnZ[i, j] -= y_lnZ[i, j].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "    x = np.arange(len(sims))\n",
    "    y = y_lnZ\n",
    "\n",
    "\n",
    "    ax.text(0.025, 1.075, catalogue_to_pretty(catalogues[0]),\n",
    "            transform=ax.transAxes,\n",
    "            verticalalignment='center', horizontalalignment='left',\n",
    "            bbox=dict(facecolor='white', alpha=0.85, edgecolor='none', pad=3),\n",
    "            zorder=5)\n",
    "\n",
    "    k = 0\n",
    "    colors = make_colours(y, 0, k)\n",
    "    ax.scatter(x, y[0, k], c=colors, s=75, zorder=-1, marker=\"o\", edgecolors=\"k\")\n",
    "\n",
    "    labels = [simname_to_pretty(sim) for sim in sims]\n",
    "    j = sims.index(\"manticore_2MPP_N128_DES_V1\")\n",
    "    labels[j] = r\"CSiBORG2\"\n",
    "    ax.set_xticks(\n",
    "        np.arange(len(sims)),\n",
    "        labels, rotation=30,)\n",
    "    # ax.set_xticks([], [])\n",
    "\n",
    "    ax.set_xlim(-0.75, len(sims) - 0.25)\n",
    "\n",
    "    ax.tick_params(axis='x', which='major', top=False)\n",
    "    ax.tick_params(axis='x', which='minor', top=False, length=0)\n",
    "    ax.tick_params(axis='y', which='minor', length=0)\n",
    "\n",
    "    ax.set_ylabel(r\"$\\Delta \\log_{10} \\mathcal{Z}$\")\n",
    "\n",
    "    kwargs = {\"marker\": \"o\", \"linestyle\": \"None\", \"markersize\": 10, \"markeredgecolor\": \"k\"}\n",
    "    gold_circle = mlines.Line2D([], [], color='#FFD700', **kwargs, label=r'1\\textsuperscript{st} place')\n",
    "    silver_circle = mlines.Line2D([], [], color='#C0C0C0', **kwargs, label=r'2\\textsuperscript{nd} place')\n",
    "    bronze_circle = mlines.Line2D([], [], color='#CD7F32', **kwargs, label=r'3\\textsuperscript{rd} place')\n",
    "    olive_circle = mlines.Line2D([], [], color='olivedrab', **kwargs, label=r'Remaining')\n",
    "\n",
    "    fig.legend(handles=[gold_circle, silver_circle, bronze_circle, olive_circle], \n",
    "               loc='upper center', bbox_to_anchor=(0.7, 1.075), ncol=2, frameon=False)\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fname = \"../../plots/logZ_harry.pdf\"\n",
    "    # if plot_smoothed:\n",
    "    #     fname = fname.replace(\".pdf\", \"_smoothed.pdf\")\n",
    "    fig.savefig(fname, dpi=500, bbox_inches='tight')\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Dependence of the evidence on smoothing scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "ksmooth = [0, 1, 2, 3, 4]\n",
    "scales = [0, 2, 4, 6, 8]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\"]\n",
    "catalogues = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\"]\n",
    "\n",
    "y = np.full((len(sims), len(catalogues), len(ksmooth)), np.nan)\n",
    "for i, simname in enumerate(sims):\n",
    "    for j, catalogue in enumerate(catalogues):\n",
    "        for n, k in enumerate(ksmooth):\n",
    "            fname = paths.flow_validation(\n",
    "                        fdir, simname, catalogue, inference_method=\"mike\",\n",
    "                        sample_alpha=True, smooth=k,\n",
    "                        zcmb_max=zcmb_max)\n",
    "            if not exists(fname):\n",
    "                raise FileNotFoundError(fname)\n",
    "\n",
    "            y[i, j, n] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "        y[i, j, :] -= y[i, j, :].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, simname in enumerate(sims):\n",
    "    for j, catalogue in enumerate(catalogues):\n",
    "        print(simname, catalogue, y[i, j, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "    plt.figure()\n",
    "\n",
    "    ls = [\"-\", \"--\", \"-.\", \":\"]\n",
    "    for i, simname in enumerate(sims):\n",
    "        for j, catalogue in enumerate(catalogues):\n",
    "            plt.plot(scales, y[i, j], marker='o', ms=2.5, ls=ls[i],\n",
    "                     label=catalogue_to_pretty(catalogue) if i == 0 else None, c=cols[j],)\n",
    "\n",
    "    plt.xlabel(r\"$R_{\\rm smooth} ~ [\\mathrm{Mpc} / h]$\")\n",
    "    plt.ylabel(r\"$-\\Delta \\ln \\mathcal{Z}$\")\n",
    "    plt.xlim(0)\n",
    "    plt.ylim(0)\n",
    "    plt.legend()\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/smoothing_comparison.pdf\", dpi=450)\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. External flow consistency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CF4\", \"CLONES\"]\n",
    "# sims = [\"Carrick2015\", \"Lilow2024\", \"CF4\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "# cats = [[\"LOSS\", \"Foundation\"], \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "# cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_not2MTForSFI_i\"]\n",
    "\n",
    "X = {}\n",
    "\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"bayes\",\n",
    "            sample_alpha=True, zcmb_max=0.05)\n",
    "\n",
    "        if not exists(fname):\n",
    "            raise FileNotFoundError(fname)\n",
    "\n",
    "        with File(fname, 'r') as f:\n",
    "            X[f\"{sim}_{cat}\"] = np.linalg.norm(f[f\"samples/Vext\"][...], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = paths.flow_validation(\n",
    "    fdir, \"CF4\", \"CF4_TFR_i\", inference_method=\"bayes\",\n",
    "    sample_alpha=True, zcmb_max=0.05)\n",
    "\n",
    "with File(fname, 'r') as f:\n",
    "    x = f[\"samples/Vext\"][...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(3.5, 2.65 * 1.25))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    for k, cat in enumerate(cats):\n",
    "        i, j = k // 2, k % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            sns.kdeplot(X[f\"{sim}_{cat}\"], fill=True, bw_adjust=0.75, ax=ax,\n",
    "                        label=simname_to_pretty(sim) if i == 0 else None)\n",
    "\n",
    "        ax.text(0.725, 0.85, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='center',\n",
    "            bbox=dict(facecolor='white', alpha=0.5, edgecolor='none'))\n",
    "\n",
    "        ax.set_ylabel(None)\n",
    "        ax.set_yticklabels([])\n",
    "        ax.set_xlim(0)\n",
    "\n",
    "    handles, labels = axs[0, 0].get_legend_handles_labels()\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 1.1),\n",
    "               ncol=3)\n",
    "\n",
    "    for i in range(2):\n",
    "        axs[-1, i].set_xlabel(r\"$|\\mathbf{V}_{\\rm ext}| ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "        axs[i, 0].set_ylabel(\"Normalised PDF\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/Vext_comparison.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. What $\\beta$ is preferred by the data? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CF4\", \"CLONES\"]\n",
    "cats = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "# cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_not2MTForSFI_i\"]\n",
    "\n",
    "X = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"bayes\",\n",
    "            sample_alpha=True, zcmb_max=0.05, sample_beta=True)\n",
    "\n",
    "        if not exists(fname):\n",
    "            raise FileNotFoundError(fname)\n",
    "\n",
    "        with File(fname, 'r') as f:\n",
    "            X[f\"{sim}_{cat}\"] = f[f\"samples/beta\"][...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "\n",
    "    fig, axs = plt.subplots(3, 2, figsize=(3.5, 2.65 * 1.8))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    for k, cat in enumerate(cats):\n",
    "        i, j = k // 2, k % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            sns.kdeplot(X[f\"{sim}_{cat}\"], fill=True, bw_adjust=0.75, ax=ax,\n",
    "                        label=simname_to_pretty(sim) if i == 0 else None)\n",
    "\n",
    "        ax.text(0.1, 0.85, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='left',\n",
    "            bbox=dict(facecolor='white', alpha=0.5, edgecolor='k')\n",
    "            )\n",
    "\n",
    "        ax.axvline(1, c=\"k\", ls=\"--\", alpha=0.75)\n",
    "        ax.set_ylabel(None)\n",
    "        ax.set_yticklabels([])\n",
    "\n",
    "    handles, labels = axs[0, 0].get_legend_handles_labels()\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 1.075),\n",
    "               ncol=3)\n",
    "\n",
    "    # for i in range(3):\n",
    "    for j in range(2):\n",
    "        axs[-1, j].set_xlabel(r\"$\\beta$\")\n",
    "\n",
    "    for i in range(3):\n",
    "        axs[i, 0].set_ylabel(\"Normalised PDF\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/beta_comparison.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What $\\sigma_v$ is preferred by the data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"manticore_2MPP_N128_DES_V1\", \"CF4\", \"CLONES\"]\n",
    "cats = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "# cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_not2MTForSFI_i\"]\n",
    "\n",
    "X = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"mike\",\n",
    "            sample_alpha=True, zcmb_max=0.05, smooth=1)\n",
    "\n",
    "        if not exists(fname):\n",
    "            raise FileNotFoundError(fname)\n",
    "\n",
    "        with File(fname, 'r') as f:\n",
    "            X[f\"{sim}_{cat}\"] = f[f\"samples/sigma_v\"][...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "\n",
    "    fig, axs = plt.subplots(3, 2, figsize=(3.5, 2.65 * 1.8))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    for k, cat in enumerate(cats):\n",
    "        i, j = k // 2, k % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            sns.kdeplot(X[f\"{sim}_{cat}\"], fill=True, bw_adjust=0.75, ax=ax,\n",
    "                        label=simname_to_pretty(sim) if i == 0 else None)\n",
    "\n",
    "        ax.text(0.9, 0.85, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='right',\n",
    "            # bbox=dict(facecolor='white', alpha=0.5, edgecolor='k')\n",
    "            )\n",
    "\n",
    "        ax.set_ylabel(None)\n",
    "        ax.set_yticklabels([])\n",
    "\n",
    "        xmin = ax.get_xlim()[0]\n",
    "        if xmin < 0:\n",
    "            ax.set_xlim(0)\n",
    "\n",
    "    handles, labels = axs[0, 0].get_legend_handles_labels()\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 1.075),\n",
    "               ncol=3)\n",
    "\n",
    "    # for i in range(3):\n",
    "    for j in range(2):\n",
    "        axs[-1, j].set_xlabel(r\"$\\sigma_v ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "    for i in range(3):\n",
    "        axs[i, 0].set_ylabel(\"Normalised PDF\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/sigmav_comparison.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Bulk flow in the simulation rest frame "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"manticore_2MPP_N128_DES_V1\", \"CLONES\", \"CF4\"]\n",
    "\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "    plt.figure()\n",
    "    for i, sim in enumerate(sims):\n",
    "        r, B = get_bulkflow_simulation(sim, convert_to_galactic=True)\n",
    "        B = B[..., 0]\n",
    "\n",
    "        if sim == \"Carrick2015\":\n",
    "            B *= 0.43\n",
    "\n",
    "        if sim in [\"Carrick2015\", \"Lilow2024\", \"CLONES\"]:\n",
    "            plt.plot(r, B[0], label=simname_to_pretty(sim), color=cols[i])\n",
    "        else:\n",
    "            ylow, yhigh = np.percentile(B, [16, 84], axis=0)\n",
    "            plt.fill_between(r, ylow, yhigh, alpha=0.5,\n",
    "                             label=simname_to_pretty(sim), color=cols[i])\n",
    "\n",
    "    plt.xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "    plt.ylabel(r\"$|\\mathbf{B}_{\\rm sim}| ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    plt.xlim(5, 200)\n",
    "    plt.legend(ncols=2)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/bulkflow_simulations_restframe.pdf\", dpi=450)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Bulk flow in the CMB frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\"]\n",
    "# cats = [[\"LOSS\", \"Foundation\"], \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "\n",
    "data = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"bayes\",\n",
    "            sample_alpha=True, zcmb_max=0.05)\n",
    "        data[f\"{sim}_{cat}\"] = get_bulkflow(fname, sim)\n",
    "\n",
    "def get_ax_centre(ax):\n",
    "    # Get the bounding box of the specific axis in relative figure coordinates\n",
    "    bbox = ax.get_position()\n",
    "\n",
    "    # Extract the position and size of the axis\n",
    "    x0, y0, width, height = bbox.x0, bbox.y0, bbox.width, bbox.height\n",
    "\n",
    "    # Calculate the center of the axis\n",
    "    center_x = x0 + width / 2\n",
    "    center_y = y0 + height / 2\n",
    "    return center_x, center_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    nrows = len(sims)\n",
    "    ncols = 3\n",
    "\n",
    "    figwidth = 8.3\n",
    "    fig, axs = plt.subplots(nrows, ncols, figsize=(figwidth, 1.15 * figwidth), sharex=True, )\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "    # fig.suptitle(f\"Calibrated against {catalogue}\")\n",
    "\n",
    "    for i, sim in enumerate(sims):\n",
    "        for j, catalogue in enumerate(cats):\n",
    "            r, B = data[f\"{sim}_{catalogue}\"]\n",
    "            c = cols[j]\n",
    "            for n in range(3):\n",
    "                ylow, ymed, yhigh = np.percentile(B[..., n], [16, 50, 84], axis=-1)\n",
    "                axs[i, n].fill_between(\n",
    "                    r, ylow, yhigh, alpha=0.5, color=c, edgecolor=c,\n",
    "                    label=catalogue_to_pretty(catalogue) if i == 1 else None)\n",
    "\n",
    "\n",
    "    # CMB-LG velocity\n",
    "    kwargs = {\"color\": \"mediumblue\", \"alpha\": 0.5, \"zorder\": 10}\n",
    "    for n in range(len(sims)):\n",
    "        axs[n, 0].fill_between([r.min(), 15.], [627 - 22, 627 - 22], [627 + 22, 627 + 22], label=\"CMB-LG\" if n == 0 else None, **kwargs)\n",
    "        axs[n, 1].fill_between([r.min(), 15.], [276 - 3, 276 - 3], [276 + 3, 276 + 3], **kwargs)\n",
    "        axs[n, 2].fill_between([r.min(), 15.], [30 - 3, 30 - 3], [30 + 3, 30 + 3], **kwargs)\n",
    "\n",
    "    # LCDM expectation\n",
    "    Rs,mean,std,mode,p05,p16,p84,p95 = np.load(\"/mnt/users/rstiskalek/csiborgtools/data/BulkFlowPlot.npy\")\n",
    "    m = Rs < 175\n",
    "    kwargs = {\"color\": \"black\", \"zorder\": 0, \"alpha\": 0.25}\n",
    "    for n in range(len(sims)):\n",
    "        axs[n, 0].fill_between(\n",
    "            Rs[m], p16[m], p84[m],\n",
    "            label=r\"$\\Lambda\\mathrm{CDM}$\" if n == 0 else None, **kwargs)\n",
    "\n",
    "    for n in range(3):\n",
    "        axs[-1, n].set_xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "\n",
    "    for n in range(len(sims)):\n",
    "        axs[n, 0].set_ylabel(r\"$|\\mathbf{B}| ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "        axs[n, 1].set_ylabel(r\"$\\ell ~ [\\mathrm{deg}]$\")\n",
    "        axs[n, 2].set_ylabel(r\"$b ~ [\\mathrm{deg}]$\")\n",
    "\n",
    "    for i, sim in enumerate(sims):\n",
    "        ax = axs[i, -1].twinx()\n",
    "        ax.set_ylabel(simname_to_pretty(sim), rotation=270, labelpad=7.5)\n",
    "        ax.set_yticklabels([])\n",
    "\n",
    "    # Watkins numbers\n",
    "    # for n in range(len(sims)):\n",
    "    #     rx = 150\n",
    "\n",
    "    axs[0, 0].set_xlim(r.min(), r.max())\n",
    "\n",
    "    axs[0, 0].legend()\n",
    "    handles, labels = axs[1, 0].get_legend_handles_labels()  # get the labels from the first axis\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 0.975), ncol=len(cats) + 2)\n",
    "\n",
    "    fig.tight_layout(rect=[0, 0, 0.95, 0.95], h_pad=0.01)\n",
    "    fig.savefig(f\"../../plots/bulkflow_CMB.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Full vs Delta comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogue = \"CF4_TFR_i\"\n",
    "simname = \"csiborg2X\"\n",
    "zcmb_max=0.05\n",
    "sample_beta = True\n",
    "sample_alpha = True\n",
    "\n",
    "fname_bayes = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"bayes\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "fname_mike = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"mike\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "\n",
    "X = []\n",
    "labels = [\"Full posterior\", \"Delta posterior\"]\n",
    "for i, fname in enumerate([fname_bayes, fname_mike]):\n",
    "    samples = get_samples(fname)\n",
    "    if i == 1:\n",
    "        print(samples.keys())\n",
    "\n",
    "    X.append(samples_to_getdist(samples, labels[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = [f\"a_{catalogue}\", f\"b_{catalogue}\", f\"c_{catalogue}\", f\"e_mu_{catalogue}\",\n",
    "          \"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"alpha_{catalogue}\"]\n",
    "# params = [\"beta\", f\"a_{catalogue}\", f\"b_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "# params = [\"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"mag_cal_{catalogue}\", f\"alpha_cal_{catalogue}\", f\"beta_cal_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 11})\n",
    "    g = plots.get_subplot_plotter()\n",
    "    g.settings.figure_legend_frame = False\n",
    "    g.settings.alpha_filled_add = 0.75\n",
    "    g.settings.fontsize = 12\n",
    "\n",
    "    g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "    # plt.gcf().suptitle(catalogue_to_pretty(catalogue), y=1.025)\n",
    "    plt.gcf().tight_layout()\n",
    "    plt.gcf().savefig(f\"../../plots/method_comparison_{simname}_{catalogue}.pdf\", dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guilhem plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manticore vs linear comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"Carrick2015\", \"csiborg2X\"]\n",
    "catalogues = [\"LOSS\", \"Foundation\", \"2MTF\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "y_lnZ = np.full((len(catalogues), len(sims)), np.nan)\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    for j, simname in enumerate(sims):\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, simname, catalogue, inference_method=\"mike\",\n",
    "            sample_alpha=simname != \"IndranilVoid_exp\",\n",
    "            zcmb_max=zcmb_max)\n",
    "\n",
    "        y_lnZ[i, j] = - get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "    # y_lnZ[i] -= y_lnZ[i].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayes_factor = y_lnZ[:, 1] - y_lnZ[:, 0]\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    plt.figure()\n",
    "\n",
    "    sns.barplot(x=np.arange(len(catalogues)), y=bayes_factor / np.log(10), color=\"#21456D\")\n",
    "    plt.xticks(\n",
    "        np.arange(len(catalogues)),\n",
    "        [catalogue_to_pretty(cat) for cat in catalogues],\n",
    "        rotation=35, fontsize=\"small\", minor=False)\n",
    "    plt.ylabel(r\"$\\log \\left(\\mathcal{Z}_{\\rm Manticore} / \\mathcal{Z}_{\\rm linear}\\right)$\")\n",
    "    plt.tick_params(axis='x', which='both', bottom=False, top=False)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/manticore_vs_carrick.png\", dpi=450)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All possible things"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dipole magnitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "sim = \"IndranilVoid_gauss\"\n",
    "\n",
    "X = []\n",
    "for cat in cats:\n",
    "    fname = paths.flow_validation(\n",
    "        fdir, sim, cat, inference_method=\"mike\",\n",
    "        sample_mag_dipole=False,\n",
    "        sample_alpha=False, zcmb_max=0.05)\n",
    "    \n",
    "    if not exists(fname):\n",
    "        raise FileNotFoundError(fname)\n",
    "\n",
    "    samples = get_samples(fname, convert_Vext_to_galactic=False)\n",
    "\n",
    "    # keys = list(samples.keys())\n",
    "    # for key in keys:\n",
    "    #     if cat in key:\n",
    "    #         value = samples.pop(key)\n",
    "    #         samples[key.replace(f\"_{cat}\",'')] = value\n",
    "    \n",
    "    samples = samples_to_getdist(samples, catalogue_to_pretty(cat))\n",
    "    X.append(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# params = [\"Vmag\", \"l\", \"b\", \"a_dipole_mag\", \"a_dipole_l\", \"a_dipole_b\"]\n",
    "params = [\"Vx\", \"Vy\", \"Vz\"]\n",
    "# params = [\"Vmag\", \"l\", \"b\"]\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    g = plots.get_subplot_plotter()\n",
    "    g.settings.figure_legend_frame = False\n",
    "    g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "    g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "    # plt.gcf().suptitle(catalogue_to_pretty(cat), y=1.025)\n",
    "    plt.gcf().tight_layout()\n",
    "    plt.gcf().savefig(f\"../../plots/vext_{sim}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow | catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogues = [\"LOSS\", \"Foundation\", \"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "params = [\"Vmag\", \"beta\", \"sigma_v\"]\n",
    "\n",
    "for catalogue in catalogues:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim)\n",
    "          for sim in sims]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{catalogue}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "     plt.gcf().savefig(f\"../../plots/calibration_{catalogue}.png\", dpi=500, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow | simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogues = [\"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "params = [\"Vmag\", \"l\", \"b\", \"beta\", \"sigma_v\"]\n",
    "\n",
    "for sim in sims:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim, catalogue)\n",
    "          for catalogue in catalogues]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{sim}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "     plt.gcf().savefig(f\"../../plots/calibration_{sim}.png\", dpi=500, bbox_inches='tight')\n",
    "     plt.gcf().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stacking vs marginalising CB boxes\n",
    "\n",
    "#### $V_{\\rm ext}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"csiborg2X\"\n",
    "catalogue = \"2MTF\"\n",
    "key = \"Vext\"\n",
    "\n",
    "X = [get_samples(sim, catalogue, nsim=nsim, convert_Vext_to_galactic=False)[key] for nsim in range(20)]\n",
    "Xmarg = get_samples(sim, catalogue, convert_Vext_to_galactic=False)[key]\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharey=True)\n",
    "fig.suptitle(f\"{simname_to_pretty(sim)}, {catalogue}\")\n",
    "fig.subplots_adjust(wspace=0.0, hspace=0)\n",
    "\n",
    "for i in range(3):\n",
    "    for n in range(20):\n",
    "        axs[i].hist(X[n][:, i], bins=\"auto\", alpha=0.25, histtype='step',\n",
    "                    color='black', linewidth=0.5, density=1, zorder=0,\n",
    "                    label=\"Individual box\" if (n == 0 and i == 0) else None)\n",
    "\n",
    "axs[i].hist(np.hstack([X[n][:, i] for n in range(20)]), bins=\"auto\",\n",
    "            histtype='step', color='blue', density=1,\n",
    "            label=\"Stacked individual boxes\" if i == 0 else None)\n",
    "axs[i].hist(Xmarg[:, i], bins=\"auto\", histtype='step', color='red',\n",
    "            density=1, label=\"Marginalised boxes\" if i == 0 else None)\n",
    "    \n",
    "axs[0].legend(fontsize=\"small\", loc='upper left', frameon=False)\n",
    "\n",
    "axs[0].set_xlabel(r\"$V_{\\mathrm{ext}, x} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1].set_xlabel(r\"$V_{\\mathrm{ext}, y} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[2].set_xlabel(r\"$V_{\\mathrm{ext}, z} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[0].set_ylabel(\"Normalized PDF\")\n",
    "fig.tight_layout()\n",
    "fig.savefig(f\"../../plots/consistency_{sim}_{catalogue}_{key}.png\", dpi=450)\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\beta$ and others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"csiborg2_main\"\n",
    "catalogue = \"Pantheon+\"\n",
    "key = \"alpha\"\n",
    "\n",
    "X = [get_samples(sim, catalogue, nsim=nsim, convert_Vext_to_galactic=False)[key] for nsim in range(20)]\n",
    "Xmarg = get_samples(sim, catalogue, convert_Vext_to_galactic=False)[key]\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f\"{simname_to_pretty(sim)}, {catalogue}\")\n",
    "for n in range(20):\n",
    "    plt.hist(X[n], bins=\"auto\", alpha=0.25, histtype='step',\n",
    "             color='black', linewidth=0.5, density=1, zorder=0,\n",
    "             label=\"Individual box\" if n == 0 else None)\n",
    "\n",
    "plt.hist(np.hstack([X[n] for n in range(20)]), bins=\"auto\",\n",
    "         histtype='step', color='blue', density=1,\n",
    "         label=\"Stacked individual boxes\")\n",
    "plt.hist(Xmarg, bins=\"auto\", histtype='step', color='red',\n",
    "         density=1, label=\"Marginalised boxes\")\n",
    "\n",
    "plt.legend(fontsize=\"small\", frameon=False, loc='upper left', ncols=3)\n",
    "plt.xlabel(names_to_latex([key], True)[0])\n",
    "plt.ylabel(\"Normalized PDF\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(f\"../../plots/consistency_{sim}_{catalogue}_{key}.png\", dpi=450)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SN/TFR Calibration consistency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# catalogues = [\"LOSS\", \"Foundation\", \"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "catalogues = [\"Pantheon+\"]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "\n",
    "for catalogue in catalogues:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim)\n",
    "          for sim in sims]\n",
    "\n",
    "     if \"Pantheon+\" in catalogue or catalogue in [\"Foundation\", \"LOSS\"]:\n",
    "          params = [\"alpha_cal\", \"beta_cal\", \"mag_cal\", \"e_mu\"]\n",
    "     else:\n",
    "          params = [\"aTF\", \"bTF\", \"e_mu\"]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{catalogue}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "     # plt.gcf().savefig(f\"../../plots/calibration_{catalogue}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $V_{\\rm ext}$ comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogues = [\"LOSS\"]\n",
    "# sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "sims = [\"Carrick2015\"]\n",
    "params = [\"Vmag\", \"l\", \"b\"]\n",
    "\n",
    "for sim in sims:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim, catalogue)\n",
    "          for catalogue in catalogues]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{simname_to_pretty(sim)}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "    #  plt.gcf().savefig(f\"../../plots/calibration_{sim}.png\", dpi=500, bbox_inches='tight')\n",
    "     plt.gcf().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Bulk flow in the simulation rest frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Carrick2015\", \"csiborg1\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "convert_to_galactic = False\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
    "cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "for i, sim in enumerate(sims):\n",
    "    r, B = get_bulkflow_simulation(sim, convert_to_galactic=convert_to_galactic)\n",
    "    if sim == \"Carrick2015\":\n",
    "        if convert_to_galactic:\n",
    "            B[..., 0] *= 0.43\n",
    "        else:\n",
    "            B *= 0.43\n",
    "\n",
    "    for n in range(3):\n",
    "        ylow, ymed, yhigh = np.percentile(B[..., n], [16, 50, 84], axis=0)\n",
    "        axs[n].fill_between(r, ylow, yhigh, color=cols[i], alpha=0.5, label=simname_to_pretty(sim) if n == 0 else None)\n",
    "\n",
    "axs[0].legend()\n",
    "if convert_to_galactic:\n",
    "    axs[0].set_ylabel(r\"$B ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    axs[1].set_ylabel(r\"$\\ell_B ~ [\\degree]$\")\n",
    "    axs[2].set_ylabel(r\"$b_B ~ [\\degree]$\")\n",
    "else:\n",
    "    axs[0].set_ylabel(r\"$B_{x} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    axs[1].set_ylabel(r\"$B_{y} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    axs[2].set_ylabel(r\"$B_{z} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "for n in range(3):\n",
    "    axs[n].set_xlabel(r\"$R ~ [\\mathrm{Mpc}]$\")\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../../plots/bulkflow_simulations_restframe.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bulk flow in the CMB rest frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"csiborg2_main\"\n",
    "catalogues = [\"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharex=True)\n",
    "cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "# fig.suptitle(f\"Calibrated against {catalogue}\")\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    r, B = get_bulkflow(sim, catalogue, sample_beta=True, convert_to_galactic=True,\n",
    "                        weight_simulations=True, downsample=3)\n",
    "    c = cols[i]\n",
    "    for n in range(3):\n",
    "        ylow, ymed, yhigh = np.percentile(B[..., n], [16, 50, 84], axis=-1)\n",
    "        axs[n].plot(r, ymed, color=c)\n",
    "        axs[n].fill_between(r, ylow, yhigh, alpha=0.5, color=c, label=catalogue)\n",
    "\n",
    "\n",
    "# CMB-LG velocity\n",
    "axs[0].fill_between([r.min(), 10.], [627 - 22, 627 - 22], [627 + 22, 627 + 22], color='black', alpha=0.5, zorder=0.5, label=\"CMB-LG\", hatch=\"x\")\n",
    "axs[1].fill_between([r.min(), 10.], [276 - 3, 276 - 3], [276 + 3, 276 + 3], color='black', alpha=0.5, zorder=0.5, hatch=\"x\")\n",
    "axs[2].fill_between([r.min(), 10.], [30 - 3, 30 - 3], [30 + 3, 30 + 3], color='black', alpha=0.5, zorder=0.5, hatch=\"x\")\n",
    "\n",
    "# LCDM expectation\n",
    "Rs,mean,std,mode,p05,p16,p84,p95 = np.load(\"/mnt/users/rstiskalek/csiborgtools/data/BulkFlowPlot.npy\")\n",
    "m = Rs < 175\n",
    "axs[0].plot(Rs[m], mode[m], color=\"violet\", zorder=0)\n",
    "axs[0].fill_between(Rs[m], p16[m], p84[m], alpha=0.25, color=\"violet\",\n",
    "                    zorder=0, hatch='//', label=r\"$\\Lambda\\mathrm{CDM}$\")\n",
    "\n",
    "for n in range(3):\n",
    "    axs[n].set_xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "\n",
    "axs[0].legend()\n",
    "axs[0].set_ylabel(r\"$B ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1].set_ylabel(r\"$\\ell_B ~ [\\mathrm{deg}]$\")\n",
    "axs[2].set_ylabel(r\"$b_B ~ [\\mathrm{deg}]$\")\n",
    "\n",
    "axs[0].set_xlim(r.min(), r.max())\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(f\"../../plots/bulkflow_{sim}_{catalogue}.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smoothing scale dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simname = \"Carrick2015\"\n",
    "catalogue = \"Pantheon+\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Goodness-of-fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scales = [0, 4, 8, 16, 32]\n",
    "\n",
    "y = np.asarray([get_gof(\"BIC\", simname, catalogue, ksmooth=i)\n",
    "                for i in range(len(scales))])\n",
    "ymin = y.min()\n",
    "\n",
    "y -= ymin\n",
    "y_CF4 = get_gof(\"BIC\", \"CF4\", catalogue) - ymin\n",
    "y_CF4gp = get_gof(\"BIC\", \"CF4gp\", catalogue) - ymin\n",
    "\n",
    "plt.figure()\n",
    "plt.axhline(y[0], color='blue', label=\"Carrick+2015, no smoothing\")\n",
    "plt.plot(scales[1:], y[1:], marker=\"o\", label=\"Carrick+2015, smoothed\")\n",
    "\n",
    "plt.axhline(y_CF4, color='red', label=\"CF4, no smoothing\")\n",
    "\n",
    "plt.xlabel(r\"$R_{\\rm smooth} ~ [\\mathrm{Mpc}]$\")\n",
    "plt.ylabel(r\"$\\Delta \\mathrm{BIC}$\")\n",
    "plt.legend(ncols=1)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../../plots/test_smooth.png\", dpi=450)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"Carrick2015\"\n",
    "catalogue = \"Pantheon+\"\n",
    "\n",
    "\n",
    "X = [samples_to_getdist(get_samples(sim, catalogue, ksmooth=ksmooth), ksmooth)\n",
    "     for ksmooth in [0, 1, 2]]\n",
    "\n",
    "params = [\"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\"]\n",
    "# if \"Pantheon+\" in catalogue or catalogue in [\"Foundation\", \"LOSS\"]:\n",
    "#      params += [\"alpha_cal\", \"beta_cal\", \"mag_cal\", \"e_mu\"]\n",
    "# else:\n",
    "#      params += [\"aTF\", \"bTF\", \"e_mu\"]\n",
    "\n",
    "\n",
    "\n",
    "g = plots.get_subplot_plotter()\n",
    "g.settings.figure_legend_frame = False\n",
    "g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "plt.gcf().suptitle(f'{catalogue}', y=1.025)\n",
    "plt.gcf().tight_layout()\n",
    "plt.gcf().savefig(f\"../../plots/calibration_{catalogue}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}