{
 "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",
    "\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": [
    "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\",]\n",
    "catalogues = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "y_BIC = np.full((len(catalogues), 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, 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_BIC[i, j] = get_gof(\"BIC\", fname)z\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": [
    "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.5 * figwidth))\n",
    "    fig.subplots_adjust(hspace=0)\n",
    "\n",
    "    x = np.arange(len(sims))\n",
    "    y = y_lnZ\n",
    "    for n in range(len(catalogues)):\n",
    "        i, j = n // 3, n % 3\n",
    "        ax = axs[i, j]\n",
    "        ax.text(0.1, 0.875, catalogue_to_pretty(catalogues[n]),\n",
    "            transform=ax.transAxes, #fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='left',\n",
    "            bbox=dict(facecolor='white', alpha=0.5),\n",
    "            )\n",
    "        ax.scatter(x, y[n], c=\"k\", s=7.5)\n",
    "\n",
    "        y_min, y_max = ax.get_ylim()\n",
    "        y_offset = (y_max - y_min) * 0.075  # Adjust the fraction (0.05) as needed\n",
    "\n",
    "        for k, txt in enumerate(y[n]):\n",
    "                ax.text(x[k], y[n, k] + y_offset, f\"({y[n, k]:.1f})\",\n",
    "                        ha='center', fontsize=\"small\")\n",
    "\n",
    "        ax.set_ylim(y_min, y_max + 2 * y_offset)\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=35)\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",
    "        axs[i, 0].set_ylabel(r\"$-\\Delta \\ln \\mathcal{Z}$\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/lnZ_comparison.pdf\", dpi=500, bbox_inches='tight')\n",
    "    fig.show()"
   ]
  },
  {
   "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": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_LOSS_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:24:46\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_Foundation_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:24:56\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:26:54\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:28:43\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:20:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:19:39\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_LOSS_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:25:09\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_Foundation_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:25:47\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:27:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:28:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:19:59\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:26:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_LOSS_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:29:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_Foundation_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:42:40\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:54:29\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:10:19\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:27:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:58:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_LOSS_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:33:23\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_Foundation_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:39:06\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:59:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:16:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:39:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:01:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_LOSS_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:46:28\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_Foundation_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:49:33\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:12:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 14:07:26\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:23:47\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:19:35\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_LOSS_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 18/09/2024 15:44:36\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_Foundation_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 18/09/2024 15:48:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:49:16\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:50:20\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:52:12\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:51:03\n"
     ]
    }
   ],
   "source": [
    "sims = [\"Carrick2015\", \"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)\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": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAIBCAYAAABp6FWdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAADtKUlEQVR4nOz9d3hc53mgD99nesUU9EoUFrFJIkBJVnNskXAcJ5YTmxAj7ya7m0SknP3l+3a/jUTRSb7dFEsiN3E28SY2Scdx7NiyBNqSJUu2RFC9kwAp9oZBJ/pgKqbP+f1xMAMM6gCYQSHPfV24BnPmnPd9z8x7nvOc532KIIqiiIyMjIzMqkGx3AOQkZGRkZkfsuCWkZGRWWXIgltGRkZmlSELbhkZGZlVhiy4ZWRkZFYZsuCWkZGRWWXIgltGRkZmlSELbhkZGZlVhiy4ZWRkZFYZsuCWkZGRWWXIgltGRkZmlSELbhkZGZlVhiy4ZWRkZFYZqsUc3DXkZ9gXytRYpiXXpKU8z5j2/i6Xi0ceeYSWlhYAqqurOXDgALW1tQvq3+VyYbPZGBkZwWq1zrifIAgsNtFiU1MThw4dorGxMaXd6urq5PtDhw6xc+fORfWTNt5OCA5ltw9dHpgrstvHTY57yMuoN5jVPgxmHZY8c9r7Z/o6XSjpXt/pcPjwYQ4cOADArl27kv8fPXqUp556CpfLxc6dOzl06NBih71wwd015KfuiV8QCMcWPYjZ0GuUND/9W2kJb4fDQV1dHQcOHEgKv5aWlhTBN1+sVivNzc2L/lHnoqGhAZfLhcPhmPJZa2trVvueFm8nPLMRoqPZ7UdlgIcvpi28Z7vg6+rqcLlcVFdX43A4sFqtHDhwIOVGN9PFdaPiHvLyncefJRqOZrUflUbFowd3pyW8s3GdTtdHOu1l6vp2OBwcO3aM1tZWXC4XdXV11NTU8NBDD/Hss8/S3NwMQE1NDUePHmXXrl2L6m/BgnvYFyIQjvHYFzdTnmdY1CBmomtolP/90nmGfaG0BPfevXvZs2cPe/bsSW7LxB08U1pAU1MTdrt92vYaGxtpaWmhoaEhI30tmuCQJLRrvw7mNdnpw9sBLU9KfaUhuNO54A8cOJC8KBL7Nzc3J4X5dBfXxPlyozHqDRINR7n3wW3k5Jmy0odnyMd7L55i1BtMS3Bn6zqdyL59+9i/f39a7c6n75muYZfLxf79+wHpZrB3716OHTvGnj17Up6gM3VzWpSpBKA8z8DaopxMjGXRJEwNM9HQ0JDU1CaaHGpqajh06BD79u2jsbGR+vr6lPc1NTVJM0hTUxN79+4FpB/o+PHjU+7We/fupaGhIWMmjYQmuSwaonkNWNcvbZ8zMN8Lvrq6mu3btyeF+2wX141OTp6J3CLrcg8DmPs6TZgWALZv357cd6I58uDBg7S2tiY/m3gNr1+/nqNHj9LS0pLUqNNp1+FwJJ98QVKm0hXq0+1nt9unnNfJkydTBPlCuWEWJxMmhtnuaHv37qW1tZXGxsYUzdbhcHDgwAGOHz+e1Mwmvp+4X0NDA83Nzcl2JgvthoYGrFZrxoS21WqlsbGR1tZWmpqaOHr0aEbaXY1MvGnOhcvl4uDBgzidzqQGXltbO+UCm3xxyWSXua7TlpYWnnrqKY4fP540Yezbty+tdhPX7DPPPMPOnTtpbGxMCu102j169Cg7d+6ktbWV1tbWRa2LJW4iCerr63nkkUfYv39/Rsyui9a4VwqJiTCbbSshTGtra3G5XLhcruSXeOjQoZQvdPJ7kH7Yhx56KLl9cj979+7F4XAk76gulyt5h0+MbfJxc2nQIyMjKeN/9tlnF20fW42kc2MG6RE5sRCUeCKajsTFdezYsYyPVWZm5rpOn332WXbv3p28xvbv309VVVVaT5rTXbPzabe2tpaGhgYcDgd79+5l586dC7qGGxoaaGxsTNknMc/q6+uxWq2Lfsq7YQQ3kFyxnelHPnjwICdOnJj2s8mTaDG2qKamJnbu3JlcHJu4fSYbt8zspHNjhqk27h07drB///4pN7vpLi6ZpWGu63ShLPa33LlzJ83NzRw9epSGhgaOHDkyxTw51zXc0NAwq3dMwgy7WMF9w5hKQLrjHj58mMOHDye3ORwOWlpaOHr0KMeOHaOxsZEjR44sqP1du3bx3HPPJW1gDocj+b/Vak268u3duze5fTEcPXqUgwcPprxP11RwIzJfV6rq6mr2798/5Zi5Li6Z7DLbdbp7924OHTqUvH727dvHQw89lNwvsX2uJyWr1Zpyfc7VbmIM1dXVPP744+zfv39GJW8mGhoapiyIHj16NMW8eezYsYyYUW8owV1dXU1bWxvHjh2jpqaGmpqapC17586dSS+DRx55ZEF35+rqao4cOZL0Rti7dy9Op3PKPvv27eORRx6ZV9uJBc3EGBMuQ62trdTU1FBXV8e+ffuWzod7BTLbBT8dCXNIfX19ctt0F5fM0jLbdVpbW8uBAweS1xiQvPHu2bOHHTt20NDQQHV19ay24oRNeceOHbhcrlnbTdDS0pIcz7Fjx+alJCUE9I4dO7DZbNhsNhoaGti1axcnTpxI9ptwX10sgrjAqJHT7U7u////akncAd/5q89ze+XqX0RaVaaSwRZorFsad8CGZshP7zuZ7MdttVo5cuRI0o8bpAVHp9OJy+Vi7969PP744wDJR+CJF3xiEetGpbdtkO/9xU+XxB3wD/76KxRX5Welj5XCSrmGFyy4V2IAjkwGWaEBODLzYyUG4MgsngULbliZIe8yGUQOeb8hWIkh7zKLY1GCW0ZGRkZm6bmhFidlZGRkbgZkwS0jIyOzypAFt4yMjMwqQxbcMjIyMquMRYW8dwYGGIq4MzWWaclTW6jQF8y5X01NTUq4c4LJidIzUfBgPsyV/7m+vh6n05lMhrOi6O2EkSx7ldjyoFj2KskmAwMDuN3ZvU4tFgsFBXNfpyAFujzyyCM4HA7sdvuKyos+0/U6MeXydGHwkwugHD16NJlkKlPFEyayYMHdGRhg43t/yGg8u+6ABoWWi/f+S1rCezqWqhDCdMyV/7mlpWVK5OWKobcTHtwIwSz7cesM8OLFOYV3fX09J0+eTEm6lSCR9jaTBScyWRllORkYGOCRRx4hFMrudarVajly5MicwtvlcrFjxw4aGxuTSZyee+65rI4tXWa7Xnfs2JHM615fX5/MRzRdAZREkFhbW1syffDhw4czmj54wYJ7KOJmNB7i61UPs0a3MKE6Fx3BAZ5se4ahiHvBghsyn6R9MgtJrg4k0zxOzD62YhgZkoT2I1+HkixFTl7vgCNPSn2loXXb7fYp1UMm5qNYDJOTVy3nDT+TuN1uQqEQv/u7v5u2RjxfBgYG+MlPfoLb7Z6zj5MnTwLjmTozkSlvvsz3ek3kdU/Mj7179ybz+U9XAMXhcLB9+/bk3Kmvr894FspFZwdcoytgvbEsE2PJGjOZR6ZLrN7Q0JBM6djU1ERDQ0NSy9u7dy91dXVpT7TZkqs3NTVRXV298rPTlayBNSujkMKuXbs4dOhQiuA+dOhQ8kJaDNNVTFnusOZMUlBQQGlp6XIPIymwJ15nE5mp2MHkYieJ1KuQWtBkooliOtPpbMx0vU6uspRIYDVbO06nM3mNZyMT4k27ODlTYvX6+vqkraqxsTFZQQUkYTs5o1i6TE6uvnfv3hVj11st1NfXT9Gwp6vf19DQkEwW1NTUlNxus9k4ePBg8hE4cfHt3bs3mcckke8EpBt+gqampmSbCdNMov+6ujrq6uqSgiTRb11dHQ6Hg6NHjyaPvZkLYSRoa2sDpN/JZrMlv5PZih1MLJSQOHZyQZOWlhb27duX3L6Ym/nk63Xik1eimtJsHDlyJDmfpivgsVhuqHzc82GmxOptbW1JgXry5En2799PU1MTVqs1+bfY5OpHjx6ltraW6urqGTPbyUxPIrXunj17knbGySQ0uZaWFnbs2JF8YkoUEm5ubqahoYGjR4/y+OOPc+jQoaRgmO4CS1Q+StgsE4WIJwqahJDZt29fipBxOBxJYWK1Wqmvr78pC2FMJFHVCcYTf7W2ts5Z7CBRKOHw4cPTFjR59tlnAZJmC4fDQVNTE9u3b190MYSJgnpiAZbpmDxf9u7dy8GDB5PJzjLBTSu4ZyIhnBOPOYkFiInlyOZbIGFy/ucTJ04k07cmFjYSRW1lZieR/nbPnj0cOnQoaZOcyGyVjhKf3XHHHWkvZs5U+WgmQQPjQmYmYXIzp+edyK5du6itrU1LgZnLrJibmzujh8pirtfa2toUG/VcxTyampqora1NzouGhgb27duXUcF905pKZkusvnPnTg4cOJAsM+R0OmlsbGT37t3z7me6/M8HDhygubmZ5uZmDh06RG1trSy00yRxwbS0tNDS0jLtxXfw4EEaGhpSFowSLNViY2KcCWFy7NixpMfCzSy0m5qaOHjwYPK6a2pqwuFwsHPnzrSKHcDMBU127drF0aNHpy10ki7TXa87d+7k5MmTSU39qaeemlUWbN++naampmTfEwuTZ4obSnA/8sgjSVviXPbE2RKr7969O8WevXv3bk6ePDlvO9VMydVlFkdC654u0f1CKx1NrpgykZkERTqCJhPC5EZi+/btDA8PU1VVhc1mY9++fUlTUzrFDmDmgiaJIgU7duxIWXNIl9mu18bGRurr66mpqWHnzp1Jc9d0BVBqa2tTxme32zO+nrXg7IAtnqvUffhfl8QdsPlT/0Rtzrqs9JEJVkpy9YxyoQV21y2NO+CzzbBp9u+uvr4+pQJQYhE5IXDr6+tT/G+tVmtyDSFhX57oXXTw4MGUBazDhw+zb98+qqurk4Jk4v4TvRUSngKJ9YrJgRaTvZgmekrY7fbkuJeCq1ev8id/8idL4g74rW99i3XrVu51muBGuF4XLLhXSwCOzAJZYQE4MgtjpQXgyGSGReXjXkkh7zJZQA55vyFYaSHvMotHLqQgIyMjs8q4oRYnZWRkZG4GZMEtIyMjs8qQBbeMjIzMKkMW3DIyMjKrDFlwy8jIyKwyZMEtIyMjs8qQBbeMjIzMKkMW3DIyMjKrDFlwy8jIyKwyVlU+7rKysowkhunp6Vl0GadMtLHS2llJY8lUOz09PQSDQc6fP7+odhY69xZ6DvJxN85xra2ti55/UxBXERaLRfziF78o/vjHP15UO1/84hcXPZZMtLHS2llJY8lEOz/+8Y/FwsJCsba2dtFjKSwsXNBxCz0H+bgb57hMXQ8TWVUat06n48UXX1x0Ow8//PCKaGMltpMJVso5JY5/5plnMjGcRY1hpR+3UOTzy+xx6bKqkkwVFRXR19e33MO4YXnwwQczcmNcaWTivG70uXej/vYJlvP8stH3qlqczITdVGZmVpLWvtLQ6/U8+OCDy6q9Z5Mb/bdfjvN75plnePDBB+np6cl426tK487KXbPzVxANQPXvZLZdmRVDJubNqtRIe94Edyts/AMQhOUezU1LNubOqrJxZxznBfjFb0j//34XmMqWdzwyMpli5BL8/LPS//o8qPrS8o5HJqOsKlNJxml+ElQG6f/upuUdi4xMJjl1EPQFkFMDl36w3KORyTA3r+B2t8K1Z2DjH4F5DQy2LPeIZGQyQ9iDePUnUPlFKLkfrr8Oq8ciKpMGN6/gPv8dUJthzRfAVAEjF5d7RDIyGaH9/LcRYwF+qbOBbROEXOC+utzDkskgN6/gdjwPJZ8GpRaMZeCSJ7bMjcHgxSOc0Zj4zvDHxCxrxzbKT5Q3Ejen4B7tB08r5G2T3hsKwd8D8djyjktGZpE4R3vZ6m6j3boWd9RPa9Qv2bqHzyz30GQyyKryKunp6eHBBx9Mvn/44YcX5p85dEp6tW6QXg2FIEbBfx3M5RkY6c2HIAjM5Fm6b98+Wlokja+2tpYDBw4kPzt48CCHDh3CarVit9s5duzYrNvT4Zlnnknxt86GH+1K5ZPL3+ezYpz80s+h6H+dq6M9rDdXLpvgFgSBnTt3Jt9XV1dz6NChjPfjcrnYsWMHzc3NGW97JbKqBHdpaWlm/CFHLoNCIwlskDQSAH+3LLgzzN69e7FarSkCee/evRw6dIiWlhaOHTtGa2srQFK4z7Q9XSbf0Cfe7G90fB2vMKTSYbSspWCkhWuj1yGnCvo/WrYxzeemmy0cDgfV1dWANJ8eeeSRVS3kb05Tifuq5LMtKKX3ujzp1X/zaGZLgcvloqmpKUXDfvzxx2lqasLhcOBwOFL2T2Tfm2n7cpJ42lvJkZPBWJhC5zkGTGUgCBRorHQGB8BcCd52iPiXe4jLRkNDQ/L/2traJRHa2YycvDkFt7dzXMsGybtEqQWfLLgzSVNTU1LLmUhtbS0tLS3s2rULp9NJTU0N+/btSwrsmbYvJ4mnvZUcGv5873FuC7oRbFsAJgjuKmkH54VlGVdDQ0PyL/H0tHfvXurr66mrq0tuczgc1NfXA9JNv66uLrm9rq6OvXv3UlNTw9GjR5NtJ9rYt29f8riGhgbq6+uTwjoxhxoaGpJKQ6Kf2cYyU5/p8vDDD/Piiy9mJVXHqjKVZAxfh6SFJBAESesevb5sQ7pZaW5upqmpicbGRurq6mhra8Nqtc64XWZ6/NEAjRf+iYcRUdo2ApLgHgy7CRiK0COA8xwU3rHkY2tsbEx5f/jw4aT5LCGgE2axmXA4HDQ3N+NwONi7dy+7du3i4MGDyTWThGC2Wq3J/hKC+sCBA8m5lGgrnbFM1+dK4ebUuH1doM9P3aa1g793ecZzg1JbW8vJkyenbG9qakoxf+zcuZNDhw7x0EMP8dxzz825XWYqey/8A6WeDuKCgoBJWqfJ01gA6I2NgrEUnGeXc4hJjh07ltR4Ezdjl8s16zGJJze73Y7T6QSkAgUTNecELS0tHD58GIfDMWe7s41luj5XClkR3Lt37055397eno1uFkYsBKER0OambtfZZY07w1RXV/PQQw8lH2NBemzduXMn1dXVNDU1pSw8OhwOtm/fPuN2mel5Z+QsP+p7nf+oziVkLEFUqgHIHxPc3aEhyKmGodPLOMpx7rjjjhTt1+l0JoVmQkBONo/Z7fYp7dTU1CQXPhPC9vDhw8mb/URvlpkE+Gxjma7PlUJWBPdkL4CJCwPLzuiA9KqbJLhljXvR1NfXJ/8SNsGE61diO4w/Otvtdvbt20ddXR01NTXU19dTW1s743aZ6flO1y+o0BVwS2CYgHE8UZpJqUev1NAdHARLjSS4V0Do++OPPw5AXV0dDQ0NHD9+HBjXcPfu3cuzzz47p+B8/PHHaWlpob6+nkOHDmG329m+fTsnT55McUEF6emtrq6Ow4cPpzWWlc6S2LhXVObY0bFk+Fpb6nZdLozKgnuhzPYbT/QqmUhtbe20rmIzbZeZSjQe4+Whj/ntvLsxXvoZAxVfSH4mAAUaG53BQbCsk540fZ1Sbp4lYqZ5MZMv93TeHtXV1cn5kFj/SDDdPJmujcn9TTxuurHM1udKYEls3MJKygUcGNO4pxPcoRHJlCIjs0po8V7FHfXzgMaGKhogYEr1YCjQWOgI9IF1vbRhYGUJIJmFkRWNe3h4mOeffz55t3W5XCnvAb785S/Pu92EL+2CIyYBgoPSq9aSul039lg22rekGolM9khEUN7IkZPvjpxDp1CzJRIAIGQoTvm8SGOnyddCXGtDocuDwZNQM/9rT2ZlkRXBXVVVxWOPPZaybeJ7QRAWJLgzEjkZGAS1CRTq1O3JIJxeWXDfICRu8JmInMxYuoUM84HrIhsM5ZxtPslmUcWgaGKiSlKqzSUQC3M9NEyZdb0kuGWyylKkXMiK4F5p9qAUAoOgsU7dnlis9MueJTJTyVi6hQzzsecya+JriPR9QqfeyovN1/m9+8eDnsp0ktvrRX+nJLjbfi4tUK4k8+UNxlKkXMjq4qTH40lGKtXU1LBz507MZnM2u5yb4BBoLFO3q3MkLVx2CZRZJQyF3XQGBzAOrOc2wxBhfRFnO0fodwcptOgAMCl1FGlsnPG2UW9ZvywLlDKZJ2uLk0eOHKGyspInn3ySjz/+mG984xtYrVb+7u/+LltdpkfICZppbh6J6Ek57F1mldDikXLId7UqWKcegJwSTFo1b1/sT9lvraGEk54riHJu7huGrGjcP/3pTzl06BDNzc1UVVUltzscDh566CGqq6v5nd9ZpqrqgSEpN8l06PLkRFMyq4ZT3lZ0gha9X8QqeriiK2BrhZWPrg5y59o81uQZAdhorOBd13m6iFOhtcHQJ1C9TNefTEbIisb99NNP8/rrr6cIbSAZLffkk09mo9v0CDlBkzP9Z7pcKRxeRmYV8InXQU7ERq3JDYBXU8DmcisFFj2Hmq7wUks3cRHWGkpRC0o+9lyWigcn8tHLrFqyIrhFUSQnZ3rhaLValzcgJzg8vakEpPwlsuCWWSWc9rYScRm4O88DgE+Tj1JQ8JvbSqmttPPupX7euzyAVqGiSl9Ei/vaWATlJ8s8cpnFkhXBnZubu6jPs4YoSoVT1TNo3PpCSXCL8SUdlozMfAnGwlzxd+Md0HGbcYSAykJUKS1IqpRKtlXlsqHYwlsX+4mJItWGEs74HMTNVVJ2zJB7mc9AZjFkxcbd2to6JeBmIsuWXzkagHhY8uOeDkOR9PloPxiLp99HRmYFcMHfQYw4Gq+Z9QX9eCmYss+GEguXet10DvmpNBTwq1iQAZ2dIpBSvBbfu+TjlskMWRHcFouFb3zjG7N+viyEXdLrTII7Iaw9Dllwy6xoznjbQIQKbSEFoW7cupIp+xRadOjVSq72ebhvkyTYzxGnSFDC8FlZcK9iVlUAzqJD3kMj0utMNm7D2OR3X5Mn9Q1AJkPeM5JuIYN87LyKYtTItlI7Vk83PTm3TdlHEASKrHpa+718bmsJeZocroT62WmqWDG5uW9ksplyIWt+3O3t7Xz3u9+lo6MjY20uunxUyCW9qmbQuFU6qaSZ68rC2pdZUWSydNRKK112rOc8gt/EZ0vCqMQIHm3RtPsVWfR0DwWIxUVKNXlc8fdIwTfD55Z4xDcf2SxdlhXBffz4cWpra3nuueeora3ljTfeyEY38ychuGfSuAG0pfDPP4A6Hfz5f1kR+YtlZCYSicZwxLooVRZSHu8GmFFwF1j0ROIx+lwBynR5XBvtQTRXSTZueW6vWrIiuJ944gmam5t57bXXOHHiRDJZ+bKTsHGrjDPv84YLzvTA/V+An38fXk2/ZFYkHqU/YY6RkckS32/5hLg6xJ35VeSOthFSGgiqpveUyjNrUAA9zlFKdXn4YkGc+lwpniGR4lhm1ZE1P+5E8E11dfXKKaQQcoFCJVV0n47uLvioDTaKsHsPbLgNfvrdtJoOxELUffhfKXprNwfans3cmGVkJvGDCx8BsM1WSd5oKx5tyYxJo1RKJVajlh7nKBVjCaeuKsYu+5FLSzJemcyzJH7cK6aQQtgthbvPNJ4XnweDAdYA3g6o+zQ0vw2B0TmbPtT9Mhf9HXzWdhv7r36Pj1wXMzt2GRkgEo3TMnoZY9SMSWEk33cNt252D6i8HC3dzlFMSj0FGisfR30gKGFEnqOrlSXx43a5XHz3u6ma6x/90R9lo+vZCblmNpOEw/D6Mai7A5RvSS6Bt2yDSBhOvw9375z+uDG+1/Mq91m38GfVX6XjwgB/fu37HNs+fckuGZkUIhF44V+hvAY+tWPWXZvbhgnk9bBZVY4m6sMa6uFq7q/NekyuScvJ/mFiosg6Qwnve67wJ8YSBNflTJ6FzBKyJH7cFouF73znO8n3giAsj+AOu0E9g+D+8D3w+eDOe6DzvCS4N/866E1w7uNZBfdVfw9nfW38Vc3voxQU/G7RZ3iy7Rku+Tu5xViRpZORWTVEIqBWz/z5Nx+Hf/8/0v+HXoN76mfc9dUr1xCtbm43fZoCv+T9NKIvn7X7PLOWSDzOoCfIreZq3nNdwKO1YnFfne+ZyKwQVpUf96IJu2fWuI+9ChVrIL8AhorA2wYKBVSuh3MnZm3254PvoxXUbM/ZAMCv2W7l/3b+nB9cb+LJdX+Q6bOQWS24R+BPvig9sf3Gw/BX/wJaXeo+7Vfgx9+CXXvg/En4y0fg+XNgmN5l9dXBZiiBW7TVFA0/T1ihx6uZGjU5kVyT1Od1Z4BtVWWs0Rfyvu88nxu5hDIjJyqz1CxJseAVQ9gDKsPU7cPD0HwC6u6U3hsKwTPmf75mHVycPZvaCwPvU5ezDr1SA4BGoeI+2xae63tr5SzMyiyKRADOxJJUc/LdJ+HSafjtP4Cmn8K+r0I0mrrPt/8SLHao/wr8/v8PhgfgkXr41/8NLe9Ocdm7EL+COWzDrDBS4j3LsH4Nze85+Pdvvc4HTRennW9atZIcvYbrI6MIwH8oeoBrCgG87RCPzfu7kEmPZ555hgcffHB1BeCsSEKu6U0lr78GSiXctk16ry+QKuVEA5LdsbcDPK5pm+wODvK+6wL327akbL/PuoXWQC+XR+VsgzcC8w7AEUV4+cdw7+fhi/8RHv2f8OZLsLsOvncQ+rrh3V/BKz+GL/4eqDVQWAr//QD4PZJA/0/3w5/9Z4hLSc9G/GE8lh7KKEUQY5R4zvKrnkpa3r1Gjs3AuZMdXPqke9rh5Jo0XB+RCgoXaCxYc9ahFGPE5WyYWWPVBeBkiwVpPROZTuMWRclMsmkL6PXSNkOh9OrtlAQ3wNXpQ4R/1Ps6GkHFfdZUwV2bsxatoOblwY8XNlaZRZNNjWdOulph8DpsuUN6f/vdsO//gMEM//Q/ob4cvvYbsPUu+PRvjh+3fiv82T/B/30J/vAJ+MUP4Yd/D8CHHT2IZi9r9WXk+a8RD4d56YSCkjW5bL9/HaWVuTS/c5VodKoWnWvS0TMySkIfz7NKZr22/vez+CXIZItVJbgXHXYcdk9NMNXRDp0dsK1ufJt+zGbo74aiclCp4MqZKc3FxTiHu1/h07atmFT6lM+0CjW3mat5dViuqr1cZFPjmZOrYyHlFWvHt63dDP+fv4G/Pwp7/kz6+5O/kdZSJiMIcO+vwwO/Dd/5a/C6Od5zAYCNpjJKvGf5ZVc5gWCMDbdK57duSwkBf5ir56bWTc3L0RIIRxnxhwGwW9YSB/oGV+h6lMysrCrBvWgi02jcH7wHOh2s2zC+TamX/L19XaBSQ0kVXJ6afP7YcAuOQC/vHdPw7PvtUz6vy1nHOyPnCMRCGT4RmRWP4wIYcyT79WT0RvjUTulPNYd/wBcehoAfXvwBp0baQRQoVOdS5D7HC21VFJfbMZikgDKjWUdRuY0zH7URj6fauvPN0j49w1JMglZlxKnU4pfjDVYlWfEqeeKJJ2YNusnNzeVP//RPs9H1zIgihL1TbdzNJ2DtuqkXkD5/vHBwWRVcPj2lySfPvYDgNbFGVcL337qGRqWgIteAOxDh7g0F3GHZwLe7f8FbI2f4fN4d2TkvmZVJTxsUzBzRmDa2fLj9HvjpEa799m+i05tQCUoGW7vo963j7k+lepTUbCzivdcu4rjUy9pN46leDVo1Jq2aHqefrRVWAEY0JtSyjXtVkhWNu6amhurqaqqrq2lubmZkZAS73Y7dbqe5uZnW1tZsdDs7UT8gproDhkJw6QJUr526vz4PfJ3S/xVr4crZFI8AT8TPO8Fm8r2V7N2xgXvX53P4+BX+/LnT/O+XzvNnz7RQrsmnWGPnF4MfZffcZFYefV1gzctMW/f+Olw9S+HAJXJEC4awk19etmK3KLHlpZr+rLkmCkqsnHjzCtFIqq07P0dLp3M8CtivsWAdlfOVrEayonE/8sgjyf+PHj2aEnzz2GOPsXv37mx0OzthqS4f6gmmEkerJIzXVE7dX18wnkFtzToIB8FxUVo8Av7u9HFERYxfL9iGIAg8uL2Cu9bmIyjAF4xyqOkKb18c4B7rJo72v8P/2fA1VArZa/amobcLqjbMvV86bLmTuMnK7166wHc3fYHo5RM0D+Zzx10F0z7Zbqot5+1fnuPjN69wT/3G5Pb8HD1nOpzERVAIENTaWePrQRTFlZOWQiYtsm7jPnHiBF6vN/ne7XbT1NSU7W6nkhDcEzXuq5clN8CiqdVD0OVDZFQqvrBmg7SAdHZcc/7utdfRjFqoKypLbiu06inI0VNdYGZtkZlffXKd+txa+sMj/Gp49iAemRuMwR7JzJEJVCqGt9zP751to6ItyvGff0K52U9+5fSLrkazjltuK+d8cwdtl/uT2wstWoLRGAMeyS0wZiikOBZiODiYmXHKLBlZF9xPPPEElZWV7N+/n/3797N9+/YUjXzJSAruCRp3exsUFE2/QJTwLPF2gk4PZdVw6j0ATrQOct3YykZ19Yyaym0VNs53jlAkFHKLoZxvdf48k2cjs5KJhMHrBrM1Y01+VHUn39TfCi9dxayJ8dtb3QiKmbXkyvUFFJXbeOvlM7icfgDyc3QoBYH2Qem9qJdyeA/J1XBWHVkX3I8//jjHjh1L2rife+45nn766Wx3O5XpBHdHGxQWTr+/Pg/J7jFm515/K3z8BogiT739NuhC3Je7cfpjkQq1xoHT7U4eLLib14abaRvtnXZfr9fLoUOH+Ju/+RtOnpTdB1ci84ohGBmSXjMouF/q6OITnZ2vhjr4L5s7UZhts+4vCAK33lWFVqfm2M9aiISjqJVKck062gekJ2ClSXrS9DrlajjZYNVHTtbW1lJXV8djjz1GTU0NHo9nKbpNJWnjnmAq6eyEghkEt6CSzCXeNun9pu3Q20HwynmODZ1CGVdRrSmb/ljAYtBQaNFzpmOEz9huRSuoaex/Z8p+gUCAxx57jF/96le0tbXx53/+5/z4xz9e6FnKZIl5xRCMjJkezJkpiu0Phujs6kKZH2SLfxjzmQFG1bMLbgC1WkntfWvxugJ80CS5/RXZ9LT2+xABjaGECAJhuVRfVljVkZM/+9nP+NznPsfevXsBGB4epqGhIdvdTiUySeP2eqS/vFnskMYicDuk/zduA42OzsYfE7D3UqYsRiXMvrZbmW/ik84R9EottTlr+eXQVDv39773Pa5fv86jjz7Kn/zJn/C5z32OH/zgB/z857JpZdXiHBPcOXML13R488wFxHgcf3mUns1bEBwiYpqFlswWPZu2lXP5TA/XO4YpselxB8IMeYMoFSoGVDoUnraMjFNm6ci64N63bx+vvfZasiJOVVUVTqdzQW0tKuQ97AGlTqqAA3B9LLosbxaXLUMJeFolzxKtDrbeieHN5xELBtigmz2VJkBVgYke5yguf5janHV84L5AMBZOft7V1cUvfvEL6uvrKSoqQhAEduzYwX333ceRI0fo7p4+74RMeixbyLtrzFRiyozG/WrLGUS9Eo1WjbeqCAxQ9GH61WvKa/Kx5hr56I3LFFkkO/e1PslcMqQ2opd9uVcdWRfcoiji8XiSi3hut3vBGfMWFfI+OU/J9bGL2T6L4DaVQ8QP/jEhv/3XKOu/RIXPRaV6ZjNJgqoCycf2XJeLzcY1hOIRzvrGtZuf/OQn5OTkcM8996Qc9/nPfx6z2cy///u/p3duMtOybCHvnhFQKKVF7UXSN+LiSk8v8RzQKTToYx7ia5SY+wbIuT79mslkBEFg3ZZShvo8OPs8FFh0XOuTnkCdmhxyR/sWPU6ZpSXrgvvpp5+mtrYWh8OR9Cr5+te/nu1upxLxpgru/l4wGMcTS02HeawIgvO81MSWOwgLCr54YZBy9fRVtSdiNWjINWk51zVCjaEYlaDgpEeyJw4PD/Pmm29y3333oZrk1aJWq7n//vt5++23GR4ent95yiw/nhEw5Sw+ahJ46+xFlAoloi2GXqFFH3ERyjMRMpkoOp9+uHp+cQ4Gk5YrZ3oosRm41uclJoq4dXYKg045vesqI+uCe9euXTz33HN85StfwW6389prr/HlL385291OJeSeJLj7wT5NHomJqIySnXvoNAAdXpEz9ny+eHGEHMX0ie4nU11g4pOOETQKNZW6Is6MLXa+/PLLqFQq7rhj+lD4bdu2oVAoePPNN9PqR2YF4RkBo3nRzYiiyOufnMdisqE0xNALWgxRN2G1CVd5KXmOVoSxlK9zIQgCJWvstF/pp8SiJxCJ0eMMMKrLRyPGpdzcMquGrAvu06dPU1tby9NPP011dTWHDx+mvb09291OJeJJ9Sjp7wNrGotHlnUwcAJEkfZBH6eKLdzb5kQxOSH+DKwtyqF90MeQJ0iFvoCzvjbC4TAvv/wytbW16GfQ+A0GAxs2bOCdd6Z6osiscNzOGSvYzAdH3wAdA0MIihxEdRidQosu4iakNOIuKkYVCpPTm565BKCozEY4FCXuDaBWKGgb8BI0Sk+OovPCoscrs3RkXXAnPEiOHz/Onj17sNvtSQ+TJWWyjbuvNz3BbdskFVXwtNIx6OeTSj2GcJT8jvRcqNYX56AUBD66NkSVvohzvnbeeust3G73FNv2ZDZv3sylS5dkc8lqwzMi1SpdJK+2nMGs1zESVBJVhDGgRhvzEVIaCdisRLVabJ3pL2Dn2AxodWqutw9TYNXjGPCi0BfiFxQEhmev8iSzsliytK6HDh1i//79PPbYYzgcjqXqdpyQa1xwiyIMDoItHY17reSN0vserYMuLhUriKhUlFxJL2jBoFVRXWjmrQt9rNEV4I74afzZUdavX09Bwey1Am+55RYEQeDjjz+m51o/J147R1/7UFr9yiwjnhEwLk5whyIRjp8+x/qycuLKKAiQG49InymNIIA/Lw/LPDyPBEEgt9BMT/swRRYd7YN+TCoD7So94aGpaYtlVi5ZF9wWi4X9+/fT1NTErl27AJanDmPYM15Ewe2SwpLT0bgFFdg2QP9HOPwDRFUKBoqKKbp2Pu2ub6+0ca7LhTFqwXodOts6uO++++Y8zmg0UlpSxk9/+HO+/7+e57UfvMu//PlRWl6XH2uXmnm5onpcizaVvHP+Mv5giDxrAWgkgW2PSnndw2P5dnx5eeT0D6RttgPILchhuN9Drl5FIBwlFlTRrtahdKY/n2XSY1VHTjY2NmK322lsbKSyspK2tjYOHDiQ7W6nEp6wODkwlsoyHY0bwLIW0XUZP9Jx7tJqiq+mP9G3lNtQKxVcuhphzSkFeruJdevWzXpMPBbnxGtniQ4puD7SxT0P3spX9/0m62vX8Or335E17yVmXq6oPveiTSUvf3yKdSVFuEbBaJIUHUssQBwFEYVUtd2Xn4sQj5PTm747n73AjCiC6JVuAm63iEOlx+BxQHzuG4AoinIB7DTJpjtqVtK6ejwecnJyACng5rHHHkt+VlVVlQzGWVLC7vHFyYGxiZ6Oxg1grkIQ49QYe7iEFndZDaaP3sI0PIAvd3ZzB4BOrWRDqYWPzvRQ2AraT+WhmK5c1QTefaGF7qt93Hb7bbz6UQeu+BAK5Rrq6rcw0O3kl//6Nv/5f/2OnI5zJeJ1SVVuFkhb3wAXu3r4D5+9l9cv+DDmiowA5qifsNKYdDMM5uQQ1WmxdvfgKp87rgDAaNai1akZ7nNjNWjoHQ4xqrOg9HSD6yrYp+bfcQ95uXyyjcsn2+hpHUChENiwvYoHfvdTmG0LP0+ZhZMVjdtqtXL69GmpA4UCpVKZ/Eu8XwiJx9XEX9oRlKIo+XEnBHd/P2g0kh93OhgKiaOgWtePVWFmuLQSgOKr6Sfn2VJmxdPqQIiDa83swvZqSwedl66z9b51bK7bQEFeASdPSUWHlUoF23du4XrrABc+XIaCFKuAxCNq4m9JIydFUarSvghTycsnTpFj0LOupJh+VwCtMYZW0GCIeiTBnUCQzCW2zs602xYEAVu+ib6uEfLMOnqcowzrx9I+DKfaudvOdfP9//U8//e//Yjjz3xILBrjtk9vYNNdNbR+0sl3v95IT2v/NL3IZJusaNwjIyNYLFK4bzxNP9N0SDyuzpuID8Q4qMYupoF+sNnTD5AQVPgUdiqVI1xUbiRkNOO15VF87TxXP/VAWk2sK87hnHuAWL6eHmFmL5FwMMLpNy9SUlNAQXkuAOvX3sKHJz9gNDCKQW+gqDKPsnWFvPHsh2yoq0SlycrPuGp5+OGHU0waDz744NJ1HvBDLLZgjdsbCHLs1Fnu33wLPSNBqeiBJopOkFwB/ZrclP09hYVUtJxCPRogYkgvUtOWZ+LK2R5q71ZzqstFqcaCW23EMnwG1v0uAB+98glNP/6A/DI7932plpK1BWi06mQb67at4a2fnuRHT75Ew3//PFVb0tP4ZTJDVjTuhNBeMYTd0mvSVNIPVuu8mhiO5VCKF6tSCqwYLqui5PLUyu8zoRXi2PxOQrkW+sOulJwlE7l8so1IOMra2yqS2zZt2EwsFuPDkx8kt2377EY8Tj8nXpNTcq4ovGNzLd2nuUm88P4J4nGRT21cS8egD61aSUQZQi+o0UU9hJWpxa49xUWIgkD+latp92HLMxGLxlGFI4SjMbSiji61CYalvNwXPrxG048/YPPda/nc791D5ebSFKENoDNq2fHwXeSX2Xn2b1/h0oll8BS7icmK4J5oHlEoFMm/xZpKFsxkwd17HWy5M+8/DT0BA+WxABbFuODO63KgDgbSOt7p6EQQRZwqyfbfFZpadSQajnHxYwelawvQGjTJ7SajiXXV6zj+9rHkE4wlz8z6bWt494Vm/O70xiCzBHhd0usCTCV9Iy6OvvcR92xah1mvp33AR6FFx6g4Sn4cFMQJK1PbjWk1eIqLKf3kDEIsvbB1i82AUqkgPCLVnxTCGlpVGhg+i3vIy8vffYs1m0q5/TO3zLqGolKr+LVdd1C2vpif/eMxTr0hV4xfKrIiuOPxOLFYjFgsRjweT/4l3sfSnGAZIzQmuFVGyQbZ1wf29AW3LxThekxNXjyCTZAeR4fLqlDE4xS2pueaN3StHYVez7BfOr4tMNUT4NrpDiKhCGs2TV2FvmPbnQwNDyZt3QC33i/VNHz3581pn4tMlvGNzTXd/DTucCTKkz95AYNWywO3bSEeF+kY8lNk0eONByiMSR4fIdXUG0Lfxg3oPB7WfJReeTyFUoE1z8hgjwuTVk00oOaiQgm+Dl791zdRaZTc9fmtaS18K5UK7vvSNtZtW8Mr//IWH/8q/adQmYWzZAE4Ho8n5W9JCbukV7UJ3G4IBiA3fcHd4xxlSKlGAZSNmTi8uQWEdXqKWtPTMoZbOzAXF4CowiKYaR29nvJ5LBbn/IfXKKrMQ2/UTjm+sKCIqjXVvHzspaTWrTVo2PSpGlqOX8Az7Ev7fGSyiG9sbs/TVPIvr72Bo3+Q//jAveg0anpdAUKRGIU2LQExSH4sSBwFoUmmEoCgJYe+zZuoONmMtSu9gBx7vpn+bhc2o4ZRn4J2lY6rzi1c/aSH7Ts3o9Gp525kDEEQuOPXt7D57rUc+/f3ef9FOQoz22RdcB8/fhy73Y7VasVms2G1Wqmrq8t2t6mEXNKr2gTXxyb2PDTuziE/IxppIlsTBRkEBSNF5RRdm1vjDvsDePsGsZcVYdSq0EWNUwR36+lOAr4gVZtn9vm8q+5u+gb6OHP+dHLb+rpKVGolJ16V6wZmk7Q9mhIa9zwWJ7uHnLz4UQu/XruVsjxpXrb1+1ApFBjGfLjtsQAhpXnGBfWB9evw5eVS/e57afVpzzcTCkYwxOO4RwQ6BQPH2hsoKlZQvqE47bEnEASB2z9zC1vvW88bz33EO8/fvE+BS+HVlHV3hEcffZTjx4+zbds2HnroIY4cObL0xYLDLqmAglIHiRDh2QooTKJj0EfIIhAOCRjD44EvI8UVlF2cW7twtkuJ6i0lBeQH3Lj8Bi5rxn/MWDTG2feuUFSZj9EyVaNKUFJUQklRKa+/+zq3b60FQKNVU3NrOZ+8fZnPPHQnStUSrx/cJKTt0eR1S8JVm34u7l+ePI1Bq+GeTeuT29oGvBRYtIwKkh3aHvYQVM+ScVCA/ltuoebd98jp7cVTPLvwteUZpWLDngA+UaDYcSuuYB6/uWlowbEBgiBw26c3oFAKvP3TEwgKuO9LS6ykrQCWwqtpSQopbNu2DZB+WIvFQlvbEpdKCo6Aekxb6e6SAm80U80R0xGLx2kf9BHXh3ArDZhDA8nPXEWlmFzD6D2z15FytnWiMRnRmk3kW/T4hjT4YgGuhyS3wGunOgn6QtRsndulasvGrVy5dokR13ifNbdVEPAFuXqqI61zkskiPrekbc8RYJVAFEXeOX+JW6sqUCUW7UVoG/BRZNXjjY2CKJITcREaW9ieseuCfMJ6fVoeJkqVEluuCf+QF01EQ83FWsoLrmIVFu8dsvXe9dx6/wbeajzBJ2+nX6lHJn2yLrhra2t5/fXXAdi+fTuf+9zncLlc2e42lbBrPE9JVwfkzx3tmKBreJRwLE5I6cerMmEOTxDcBZJZI7/j2qxtDLd2klMs9Zlv0RL1SNrYJX8X8Viccx9cpagqD0PO3Frauur1KBQKTp0ZfxS15puxFeZw8WPZJWvZ8brn5VFy3TnCgMvDhrJxDXnAHcAfilJsM+CN+ykUQSWGCajmcLMVwF1STJ6jPa2+7QUmhq+7qRkuQRREciuugC/9YJ7Z2HrfOmpuq+CX33ub/g45PUOmybrgfu6553jgASlI5bHHHmPv3r2cPHky292mEhrTuAHaHFCUvg3P0e9Fo40TIcqoxo45NB4p5rflElVryO2e+QkiPBrA2zdATokkuAvNeoSoihzBzHlfOx0XrhPwBqmcxpNkOrRaLWUl5Zy5kBrlVr6hmKst7UQjciWTZSWhcafJufYuBKC6aFyZcPT7UCgECi06vHE/NWMeJQG1de7uiwrQer0YnHNXE7bnmwkFIhg9SpxruhjQxqQyfeLig+YEQeDOX99CTq6Jn3/7dWJReV5mkiX3Kqmvr1+wDW3BxYJDI5IPdzAgFVAonLvsWIJrfR4sedKkC6nz0Ue9aKJjHhyCAk9+EbldM2u6ww7JfJFTUgiARq3AZtSiC+Vw1tfGlZZ27EUWTNaZbduTqaqo5krrZUKhUHJb+foiIqEonZeuz3LkzcWyFAv2uuaVYOpS93UKbVa06nEvDsegj3yzDrVSgTvuozISIqw0ElHObd7z5eUjKhRpeZfkWEyAgFPRh8vsok0UIR6G0cyEsStVSu7+rdsZuj7Ch6/IaWMzSdYF9xNPPIFCoaCqqorKykoqKysXnGRqwcWCg05J4+5ol/y4i0vSOiwaj9Ex5MdgiSIgENVKx1mD44LAk1tEbk/7jG0MX+tAZ8lBaxrXwoqsOiIjBgYHXQx2Oylbm/6NBGBNeSWxWIxrbeO2TGu+GaNFL9u5J7AsxYJ97nm5Al7u7qUsb0KyM1HyKCmySmYzb9xHecTHaBraNkBcpcRvt2OdI0+3GAf3GS0GpRGvapCQX0V3Qp/yZ67qu73Qwi13VPHu8824h7wZa/dmJ+uC+8iRI4yMjDA8PIzT6cTpdC59RZfQmOB2tEoLlIWFaR3WPRwgEhcR9EFMCj1BtY2oQo01OD6xvXmF2K93SDeESYiiyOAVB5bS1P4KrXp8gzpqruej0CjIL08zS+EYdpsdk9HMxSvjroiCIFBaU8C1Ux1y2s3lxD2StqkkFovTOThMiX389x/2hvAEwhRb9cQRGY15yY948avmqI86AV9+HtbunmnnJIAYg6FTKoIjCmz5JnxeD7GQmhGlmrigAl/6xRnSYet961FrVRx/5oO5d5ZJi6wL7h07dix/7pLQCGjMkn07vwDUmrmPQXIDVCsE/EofOYIRBAVeTT72wPgCjievCHUoiHl46uPl6PAIgRE3lrJUm3qRVY8Q0FLTW0gkX0ShnN/PIAgC5aXlXL6aGvxTsrYQ16CX4V7XvNqTySC+9Bcnu4aGicZiFE8Q3O1Dkhmu0KrHHw9QFA2gJM6oJv2buy8/H1UohGlw6qJgNAj9H6oJDiiwrIthLzQTCoSI+wEEApoc8GVO4wbJZXXbZzdy8SMH7eeX0Gx1A5N1wb13715yc3PZvXs3X/va15J/S0picbKtdV727fZBHwUWHc64C/NYjgivphB7YNwc4cmTtGn79amr8YOXHQhKRXJhMoFZp6ZUqcEc0NJpW9jTR3lpBV3Xu/D7xyMmi9bkolQrudoim0uWjXnk4u4YkARrkX1csWkf8GE3adCpFXjiXsqjUsRkUGVNewijdhsxlWpKuteoX6D/PQ3RUQHrphham4jFKs1rYTSEDj0utR68mRXcAFVbyigot/PL778jL1RmgKwL7kcffZQdO3awfft2qqurk39LRiwspXVVm6GtLW2PEhGRziE/NquSgBgkR5AmuE9bgD7iQTcWQTlqsRFVa7BdnyosBy63klNUgFKdGj4sCLAhriGijHNZ20MwPn2mwNmoKFuDKIpcujbuJ6tSqyipypcztWWBtBfGva60Ne6OgUFyDHoM2vFFx84hP4VjbqGuuJeKaIigykw8Tb9wAFGpwJ+Xh71jXADHw9D/kQoEEduWKOqxiEytVoNOr0UIB1BFdfQqVBlzCZyIFBa/lZF+Nx/84nTG21+JZHNxPOuRkxaLheeeey7b3cxMyCm9BpXg86YtuN2jEbyhKNocSaiaFQmNW9KercEu+tSbJfOJvUCyc08gFongbOukrG7rlLZFEQr9In3qGHGFSGugl83GNfM6rRxzDnabnYuXL1B32/bk9vINRbz/0mncQ14sebNE2snMi7QiJ6NR8HvnIbiHKLCMB9WEwjH63QE2lEgauCvmpS4aJqhJP+4ggaeoiNIzZ1CGw8Q0GobPqhAjAratUZSTLIUWiwnP6CgEbLRq42wPDkI0AKr0oz/TwVaQw6a7anj3hWY2faoGe5E1o+2vNBIRlKsycnL37t08//zzS59YKkFwTHD3j61op2kq6RryS//oAygQMCmkSRxUmQkrdNgCExYocwuw96QKbmdbF/FoDEv5VA8WcSSMKhinXwEGQccV38LuyGvKqzh78ZOUxcjy9cWoNErOvHN5QW3KLIJESlfj7BGOCToGhiiwjZtJup2jiCIU5EgauDvmpjAWJDBHxOR0eIoKEOJxbJ2dBAYVBPoUmCpjTOdRaLGYCHr9BNwquhVjriXezGvdIAXm6E06Xv23d+VF9EWQdcH91FNP8ZWvfAWr1bo8+biDYzbkfhcolWknl+pxjmLSKPEJHkwKIwphbMyCgE+TP6dnyeCVNjRGA3rr1Isu1j0KagU+rYAxlsPl0YXZFKvXVONyu+i+Pu4FoNaqWLOxlJbjF4iG06/+LZMBEqkPjHM/6USiMa4Pj1BonSC4h/2oVQqsY9khleFBlIhzR0xO177RQMBiwe7owHVBidoios2dXlBarCZEUSTmjDOQUMe92VknUalV1O3cjONsN9dOZ+fmcDOQdcHd3Ny8vPm4E4K7zwl5+ZLwToNup5+8HB3DMRdmRepik0+Tn6Jxe/KK0Pm9GNzj0WpDV9uwlBZOCTYSRYh2jaLI02I3a4h7DfSHR/BER+d9amUl5Wi1Wk6dTc3EtulTNfjdAZqPp1+JXiYDuMee7tIwlfQMO4mLIoUTbuw9zgB5Zi0KAeLEMYal9ubKUTIT3sJCfANqIj4BU3l8xkp9JrMBhUJBzBdCo7TgUxnA176gPtOhbF0hhWvyeOPZDxHjsta9ELIuuB966CG83mV0vA8OAQL09Kado0REpMs5Sq5ZizPuxqJIvRC92gL0US+6iAuQBDeAvUcKfQ/5/PgGhsgpnWqWEUfCiL4oygItdpMW35C0cHlldP6+s0qlkprKtZw8fSLlsdOSa2LdtgreOnoCZ59r3u3KLJCExm2aW9AmPEoKU0wlfvJMOgDcMR+FsQBhQUN0slE6TdxFhVww1qE3hZOLkdOhUCiwWEyIgVFMEQvdCiVxd/YWuAVB4PZf28Bg9wiXTy5xwrkbhKwL7qeffppdu3bR0bH4R68FhbwHh0CTA11daQvuEX+YQDhGjlkgKIanaNwe7ZgL4JhboN+WS0ylSuYscbZJ2ngisdREYl2joFEgWNXk5+iJR5XkCGYu+xcW9HDLuo30D/TR3pl6Adz+2Y3oTToav/krQoH5e63cCCx5yLtr7OkuDVNJW/8glgkeJYFwjGFviLwx+7Yz5qYkFiaoXpi2DXBdWYFHY6OEuU0SFpuZqN+PwpNDt1JNyJ3dNZL8MjtFlXm8/4tTsq17ASyJH/exY8eoqqpKqT25EBYU8h4cBswwPJS24O4ZlswWCqNUyzFHSBXcYZWJgMpE7mg7AKJCiSeviLyxnCXDjk50FjMaY2r+EVGEaKcfRZ4WQSFg0qnQa5Towiau+LtZyPStKFtDjjmHN997I2W7Rqvm176yHfewj1e+9/YCWl79LHnIu2sINLq0cnG39w9QaLMm3/eN1X/MN0uCeyg+QkksTGgB9m2Q5tpAvw1bfJjCwbmVJpvNTDwaw9kbJKzLRx9yEQ9nt6rSLXdW0+sYpOdqZnKj3ExkXXBfu3Ytpe7kktu4AwPgG1tKT1dwjy1Mjir9KBAwKqYmgPJoS8jzj6dzdeeXJNO7Otu6MBfmTzlGHA4jjsZQFkjjEQQosOjxOzV4YqMMJkqszQOFQsHtW2s5ceojhoZTCxBb8szc+fmtXPjgmpzDZCkYGQJzeoL22vV+SuzW5Pvu4VFUCgUWo2QW8USHsS9C43aNmAiF1OQbhzC4XCgikVn3t1jNCIKA3+UhxygVdDjf+8asxyyW0poCzDajvBazAJYkO+D+/fvJzc0lNzeXP/7jP16KLscZHQD3mLt63lRhOh09Tj+5OTqcYwuTgjD1a3LpS7EFu5OZAl1FpeT2tBHz+fD1D2IumtpXtMOHoFUiWMYDcoqtekJuDQoEro4u7JH+ts23odPq+Okvjk75rHJTKYUVubzx3EfyI+kiSMtMNzIIZuucbbn8owx5vJTkjucf6RlbU1GMrSAqgj0IQGAeEZMJ4jGB/l47RlMA7BoEMY5pcHDWY5RKhWQu8fnwRgqJIdDWn93cIoIgsHZbBRc/amXUG8xqX8tBNk11WRfcX/va1xBFEYfDQWtrK/F4fGlD3gP94BLBbAbD3KlTRUQ6naPkmbU4Y25MiunDl516KWCmxHsOkMqYKaNR4i0tAJgKU0ujiTGIdQVQFGikklFj5Jq16FUa9HETraNTK7+ng1qt4dP3fIaWT05y+mxLymeCIHDr/esZ7HLSdi6zyYNuJtIy07mG0vLhvtwtpd6dmBWwa3iU/DH7dkAMkhdxjYW6zz+IanDQSjSqwGb3EdHriej1mOcQ3AB5eVZifj+dw6N4NDkoPQ5C8dk19cVSs7UcUYRz713Jaj/LQTZNdVkX3CdPnuTpp5/GYrFgtVr5zne+Q1NTU7a7HSfQD8NhyEvPTDLsCxEIxyjI0eGMuzFPYyYBCCuNjOjKqHB9DICrsJSYUkXw4iVUWg06S+oFF+8LIIZiKIpS7Z+CAMU2PWGPhrZA7wJOUOKWdRupqVrLvzf+MCV/CUBBRS62whyam+RH0qwy1Ac51jl3u9DZQ45Bj90seSsFQjGGvEHycySPkt7oEJWRAH61BXGap73ZCAXVDPZbsFj9qDWSSXLUYsHcn4bgzreCKNLWMUDUUEBV2McZb3bTJ+iMWsrWFXL6zYvyE+E8WJKakxPdAT0ez9L9QGJcSgo/4Evbvt09JC0Smc0CATGUzFEyHX2mTeSPtpET7CWuUuMsWcPw9QGMBblT/Lcj7T4EkwqFaWqWgXK7kZhfhyvqwxXxz+MExxEEgR2friccCfHCKz+b8lnNrRVcO9WB3xNYUPsyaTDYC9a5A7zOtHWypiAvOUe6xjICFlikm/r1aD9rIwEC6vSCxRKIInR35aFSxbDaxm/eAasVrc+LKjD7b2806tHptXiHnPiVhVRHg5z2zl6WLxOsvb2Cwe4RrrcOzL2zDLBEhRQeeOAB/vZv/5a//du/pa6ujoMHD2a7W4nAAMRjMOCCgvRycHcM+rDqNQSU0sSf7Ao4kSFjNWGFjgrXCQD616yjMyhizkvNnRwPxohfD6As0k3bjtmgxhyXNPS2wMLMJQAmo4m7t9/DOx++TW9/aiWcys0lIAhc+CD7F+JNiShKgtsyu7D1BYJc7r7Ougk+/u2DPvQaJRaDtPYRDDgwi1F82vnlKBnstzLq15FX4Gaioh4Yi840D8ytdRcU2Yl5vVzxGdCJMfqGz85rDAuhqDIfo0XPqTcuzr2zDLAEgnvXrl0cOnSIoaEhhoaGeO655/jyl7+c7W4l/NfBC0SiUJDeRdA+5KPQIkVMSjlKZraLi4KSIUM1pZ4zAFwoW49PoaYslhoFGWuTtGjFDIJbEKDMkoMQ0SxKcAPcuuV2zEYzrxz7Rcp2nUFLSXU+596fuwK4zALweSAUmFPj/ujyNeKiyIay8Rw2bQM+Cq16BMAvBikJ9hJFiX8eGrfPq6e/z4bV5kOnT7VLxzQaQiYTpoG5NdrCwlzEWIz32iTff6WnlbCY3dQJCoXA2tsqOP/BNYL+0NwHyCyNV0ltbS1PP/00Tz/9NNu2bVuKLiV8XTAWhZxOcqlAOErvSIAiqyS4U3KUzMCIfg2myDDG8DBtccmV67b2cVuyGIeow4eiQIegnvnrLrYZIKDjkmdxNSNVShV1t9/BiVMfM+xMTaS/ZlMp11sHGOl3L6oPmWno6yKGwMuDfnqdrhl3e/PMBdYU5GEd8/EPR+J0DPkptUnvOyPXuS3kxaPNJ65IL94hMKqho70AvSGE1T697/Wo1UpO/8CchYDNZgNavZbhQQ9+pYmqsI9r/uzXMa25vYJ4LC4nR0uTrAju7du3c8cdd8z4d+edd2aj26l42sGtBJ0OLNY5d28f8hFHEqJD0+QomQ6XrgQRyPW3MjDsxqIWqG49h7VP8uCI9wakEPfS2YMyzHoV+oiJwZiTUHxxGs7mjVvQarRTgnLK1xWi0ig5975sLsk419tpNFXyrQ/P8di//IjINMUCBt0eTl5ro27teM3V1n4vsbhIWa4kuJ2Ba5THQnh15Wl16/PqcbQWo1ZHKSgcmTEfyajdhiocQu+a/aYtCALFxblEvV66ozY2Rka54M9+DIDBpKPilmI+/tUZ4rHFV5m/0cmK4D5y5AiHDx+e8vf0009js9k4efLkgtqdd8i7tx1cWsm+nUZl+dY+LyaNkhy9iqGYE6tiblesqFKHX51LXqCN/kEXplwbfoud2391FOIikcsehBw1ihz1nG2VaHNBELm6wDSvCTRqDZtv2cK7H75NZELghUqjonx9MWffuXxTrOAvZci72N3GK6ZyKvJzGfJ4+eDiVPe2lz5qQatWcXvNeO71C10urEYNFoOGOCK5o9eIoUimVZixPxGGBi20tRah1UQoKnYym4IeNJuJqdTk9M5tiisqykOMxXinx8CGsJ+LvqUJ3tp0Vw3uIR8XPpQVi7nIiuDetm1byl9NTQ3PPvssDQ0N1NTUMDIyMncj0zDvkHf3NRgGCtMrntDa56XEbsAd9xEhhlWZXtSaR1eI0dPOsNOLzWbm6p2fwX69g6IPTxIfDKGqmNt/HKDKkgsRFR8MLH7i3rr5dkYDozR/knqTrLmtnJEBD52XFu56uFpYypD39ksXGFDqqa/dSnmenbfOpS60jYZCvHziFNvXVaMdq4gUi4mc63KxJt+EADhjLjYF3YxocokrZq5xEo8p6GwvoLfHjsXqo7BkZFahDSAqlYzarViuz232MJr0mHOMXB5UoBdj+J3n5jwmE9iLLJSuLeCd55tlrXsOsmrj9ng8PProo1RVVSEIAm1tbXz7299euuLBw5dgKAjFcwtuXyhCjytAqd1IX0xafU9H4wbwaIvpdYURRRGb1YSnoJi+mk2Eroyi0CsRctPL7mbSq9AHLVwNdxJfUOaScWxWGxVla3jng7dSthdW5JKTa5J9uufJXE97Z1rbUCFSVZjP5jXlnLzqIDThaecXH50iGI5w/5Zbktsu93oYDUdZWyTNs8FQB2tiAfy6shnHEQqquXa1GJ/XQEHRCPY8XzoPkwD4cvPQezxo0yhqUlpWgMcToi+oocLXgzOyNBk+b/30Bpx9bk69ufo9TFZd5OREgW2z2XA4HDz11FNLW+09FoIOB8TiUDK3xtU6ViGn1K6nLzqEWWFEI6QncD3aYs55bWhUAmaTZMu+suFT9OgryVcOpERKzkWFppiIMsSJoda0j5mJLRu3cq3tKn3949q1IAhsqKvk0gkH7qFlTLe7ypjrae/ssIdynRK1SsXmNWWEIlGar0oZG/3BEI3vfsj2ddXJRUmAj68MkpejI88kRUzq/JcQAZ92ekXD4zZw7WoJ8biC4rIhjKb5eWBIRYTV2DvnzhZYWJSLUqnkhc4C7gm6OJ/F/NwTyS2yUn1rOW81fsyod3XHHKy6yEmbzUZbWxuNjY3s3r2btrY2Tp8+nfKXdVxXYHDscSsNwX2114PdqMGoVdEV7SdXYU27q4AqhxZfPqWmWDKoYmjYjIYIG9o/RJhHVet11kIUAQMvD35MfJF26LXV69DrDbz9YarWXX1rORqtmg9f+WRR7cuM4R7hYlxDxVhu7QJrDqW5NppOSyaGZ9/+gGA4wgO3b04eMuAOcKnXzeYyKwAiUBzoZkRlmpJ/W4wL9PbY6WgrRK8PU1I6hEYz/0RtolKJtyAfe3snQnx2U4RKpaSkvIBP+jVU+UbpGF66ubLtMxuJReO89sP3lqzP1UZWigXffvvtDA8P8/jjj0/7uSAInDhxIhtdjzP0CQwA+fmgnd5/OoGIyOXrHiryjPjFAK64h/Xq9Iv3hmMil305PFAk+R5GAwKjQwqs9ija66Pkd15joHpDWm2plQJl0TV06i7yysAJfqtw4R44KqWKTRs288HH7/Hbv/FlNBpJIKg1KjZsr+LUGxe554vbMNvm9p6RmZnBj95mWKmjomzcxHHnhrW88P4Jfvj6Oxx99yN2bNuSom2/frYPg1bFujEziSfqZHPYS5ehOqXtUZ+O7u5cwiE19jwPOZbRtE0j0+EuKcZ6/Tq2zk6clZWz7ltZWUx3Zz8vd9kpzzsOa//jwjueB3qTlu31m3n/pdPcsr2aW+6snvugm4ysCO7m5ua5d8o2w5/AoAoqqubctXckgDsQoSLXSEdYMivkq+xzHDVOx0CAqChwq7aHa7EAQz1mFGpQF2jw2goovno+bcENsDm3iJ7rQxwXWqgxFbHRWJH2sZO5bfNttJw+yYlTH3HvXfcnt99yRxWXTjp47+ctfP4/3z9LCzJzcem9NwGoqBqfa9vXVXGqtZ0fvfEemypK+eytm5Kf9ThHOd3h5L5bClCOmdHi3hY0iAT0FYgiBEa1DA5Y8LiNaHVhSsqG0GgXnw45YjDgt9vJu9qKs3INMPNdQKNRk1+czzs9cf6ytBOf+zImS/rzeDFUbSmj+2o/L3/3TYqr87HkzT/Z1o3MkgTgLAudH8JgFCrnFtwXul1olALFdj0d0R7sCkva9m2Aa71+cnQK8jVBtCODjA4pMOTGEBTgLF2DZagXg8s5d0Nj6DRK1qrWoPAb+Unvm0QW4ddttdioqqzm2JuvEZ/weKzRqdn8qbWcev2CHJCzSC5dvIhViJNjGn9yUSmV7PmNB3hs12/x+zvuR6mQLjUxDi983InNqOGWkvE1nwLvBfqURgY9ZVy7Ukrr1RICo1L4enGpMyNCO4GrrBS914O5b+5IyvXrShEVCn7Wnov/3KGMjWEuBEHgrt+4FaVKyQv/fFz2MpnEjSm4xTh80iwZDmvWzr4rIp90jIy5ZIl0Rq5TpMqb9ZiJxOMil6/7KM0z4Fda6e+xoNSI6K2SfdpdUEJEo6P08pl5ncK6YgsGZwmeaIAPXZfmdexk7th2J73916ekfN2wvRKdSUfTM9nNu3xDE4lwfthDhXlqgJVSoSAvx5yScOzj1iE6h/zcu6EwmXtbE+xB69Py1sBX6enKQ6GIU1jspLRiEHNOYFGmkekIWCwEzWbyr8ydSlWrUWMpyONUv4Fox3kp/88SodVruPdLtfRc6+eto1k2ra4yVpfgDofgT74E7782+34jl8ERAJsF7LPne+ga9tPvCbK20Mz12CBhohQq0xfcHYMBAuE4Ffk6rsW2MRzIJyc/kkzyIyoUDJVXU+i4hCaQfiV3tVLgjooiFD4Tv+w7tSjnwNLiMtaUV/Kzl39KNDquvavUKrZ9diNXTrZz7bRcIWchhJvf5ZrKRGVpyZz7uv1hXmnpZkOJhRKbJOiFWIChq5c45mwApZLS8kEKi10YjOGMC+wkgoCrtBTz0BB619wxFZWVRQgaJc+32Qm1/TxLg5qegnI7t316A++/dIrWT+b2hrlZWFWCu+fCWR78txd55v/7h7Pv2PsetAObtswZMfnhlSHMWhXluUbawt3oBW3a/tsA57u8mPUqLHotl4a2kKe+TqGqPWWf4bJq4goF5efmZ/u3GrVUa8oJKP283b24oJxP3/MZhoaHeO3NV1O2V24qoaQ6n18ceROfK/0by2pgKSInr/zqeaKCgsq1sz/ZIcLRjzpRKgQ+tU6qjhQJx2n9pINz7lpEUwclJZk1icyGLy+XiF5HwaW5k47lmnUY8uycHTZw+ZM5lKYssPnutZTUFPDzb7+OZzi7dTBXC6tKcJeKIV5cAw8Hu2FolgKj7x6FUeDW2lnb8wTCnO5wsqlMsjU6ot0UKfPSCo8HCEXiXOr2UVmgp7cnj1hcQYWlnQL/VYQJdumYWs3gmvWUXD2H3uNKq+0Em23FKCNaXus7w2K8A/Nz89l++3Z+8erP6bk+XglHEATu/q3bicdFnjnwC7wjC8sHvhLJpB/tTAE4Z5ub0QkixbmzL2a/e2mAK9fdfHpTIWoFXO8JcrrZyXDAisl+FlXuMIo0k0plBIWCkbIyLNd70Ltcc+wsUFNejNYAz5/XEhpomWP/zCIIAvd88XYEhcDPvnWM2Dzca5eTVReAkzV8bvj1h6T/T741/T6iCK+/DRY9rJl9YfL9y4MoBNhYZmUwPown7qNEnV7eboCL3V4iMZE8XR6uERP2PC8+UxHqeJASb2pk4mBFNWGdnnUfvjGnD+1EFAqBKnUZo8ZB3m3rSvu46bj7jnuxWe185/v/jH90XEDrTTp27L4Ln3uU7369kXPvX70pcpnMh2kDcK538Ik3RKXFhEIx86XUPuDjlVPdbC2zgU9B8wkPHe0B8tU9bM07SZ/RiV2xhMFpY3gLCggbDJSdOj3nnMy36NGVFNIX0PDvjd9b8vmhM2j59O/U0ds2yKv/9u6qmJ+rLgAna8RisOUOKCiFT96ffp/Lb8ClANTeOqvmHIzE+ODqIBtKLGhVCs6GrqEXtOQrbDMeMxFRFGludVNkMeDsL8ZoCmAyB4ioDAzrq8gLOLAExu+0olJF16Zacob72Pzmy6hC6Ue9bTKuQSEqeXn4QwKRhWsbKpWKL/76l/D4PPyf7/wdbs+4N4m1IIff+M+fJq/Uxs//+TjPHHhZfiydg/CrR7mgtVFTPbOC0OMc5QdvOqg05qAYUXO9J4jFomZb0RU2GU5wTgtGQYcxjUyUmUZUKhlcvw6d203V+x+gCs5WsFdgc1kpleVRPm4N8INnDhOYx5pNJsgrtXHX57dy6o2LfPDS6SXte6WxugQ3QFk11GyC5nem//wfHpO80z/9hVmb+ejqIOFojNvW2OiNDnEp3EqNumLaiu7T0dY/yoA7jDG2BrUmQn6+O3mfcOuK8WnyqHCfxDbaQcLG4bfl0X7bp8gZ6mPr6y8hRNNz81MKCrZqNhDKGebvP3mNUJrHTYfNamPXgw04R4b5m7/7S9q72pOf6U1aPv3l7Xz2oTvp7xziyNcbudLcPmNbNzMvNXfx2g//nbCgZG15am6RPleAN8728cMmB0ePdVGpsGCMarBaVWxYJ7BR+w65oU9o1dpxCTFKlUVpm+cyTTAnh95Nm9C7RtjQ9DrGoaEZ99Vr1KjK1/DZMjcftXzME3/9OL86/krKgne2qbmtgq33ruON5z7iw5dPL1m/K43VJbhVarDY4Zbb4dIpcA2nfv7GC/BuC9yxBgwzZ/YLhKO8eaGPDSUW2mnned8xbMocatRz50CORhUMDpj51UkfBqUZm8FAYfEIKfUWBIEBw3r8mjwqPC0Ue8ezq3lzC3Fsuweje5ialvRDeqv1xVTFqxg0dfA/z/yUK/2utI+dTEFeIV/d9XsYDUa++U8HaW1LXfgsXVvIb/7Br5FXaqPx73/FL//1bUKB8IL7u9Fo7ffyl3/bSH9oFJUIV3qD9DsD/OrdXr7/fBu/auqn53IYwaUiV60nN1fDuvVGSgvjFA2/jCY8wCV9MedUSipUJWiV2mU9n4DdRldtLWG9nup338M0S0X4UmMp3lI9f3BLHwY1vPDK8xz4x28wMo84hcVy66c3sOWetRx/5kNe+Ze3iIaX7saxUlhdgls/FjK85U5JQzn20/HPXvsR/OkuKAJ2NMzYhIjIS81dRGJxbBU+Xg98RIWqmHt1dVOq3YgiBAIaBgcsOK4Vcaq5mg/f28jrLRHcIS+3FOdRXOpCqZxqb4srlPSbNjCkr6Zg9Br5vnGf2UCOjetrN1Ny9Rz57emXErs9p4ZtylsJ6918p+tlXj3bteAFS5PRxK4HG8jPK+AfDv89l6+m+oprDRp+7SvbuePXt3Lm7ct8+09/QvPx80RuwotkMv/n5Qv8kest3tcXYsnJpfmci9deH2SgT7q5mS0qSsq01Kw1sGmTiZISLRq1iM35BohR3tXnckkRZ42qGLNy5mLUS0lMo6F3yyaCFgtV73+AuX/6xX8FCkLGjVw0mvjP1a3cXhWjf2iYp//hyZRF72wiCAK3f2Yjn/rCbZx55zL/8udH6bk2i7PCDYggrgYr/xh1RXaaXxqrYP5P/xPaLsFX/gjeeREunYMiBfzuw1C8fdrjY2KcVz/p5c0Lfdy9ycIH+jcpUNm5Q7s1+agaj4NrxMTQoIURp5lIRMoKoNFE0BtDjES6eb/9ElUFFraUpefvnTvahi3YTZ9xAwOmWxAFBYgi5eebsQ70cOVTOxmoWpf29zAcdfHu6GniYRU1vi38x223YTXMXahhOsLhMC+9+gKdXZ2IcfjM/Z8lLzeP6jU1VFZUoVQq8blHOf3mJdov9KDTa9h633q2PbCJ/LL00wIsJw8++CAvvvjiotqoq6ujubkZlz9M7X/9Ea9ceIz/Zt/ORvsGVNE8dDkKigo1WFRB9P5LqCIu4ko9IV0pcaURk+cUmlAv7+gsBFUmSpSF6JWzV0VaSk6ev8r2zetQxGIUXriEccSJq6QEd3kZ7uJixEmLr864C+3oNe70d3MhksP7naVEQ6P8Rv1v8em7fw2TcWluSK4BDx+88gnD111suXcdn/7ydmyFUxd6n3nmmfTz+GeYTMy/yawqwV2UY6Lv+FgRXNcwHHkSuq+AaRTWWqH+/wHt1MXFSCzG2U4X71zs57orwF3rcnHYTzEcc/KA/m7iYSM+nw7XiInhQQuRiAqtNow5ZxS9YZSIwsvIqB/HgBvHoJsSq5HaqvHItzkRRWzBLuyBDsJKI72mLbj1JQjxOGUXWrD3dTFUVk3vus24CkoRVXO7hXnjfj70ncOn8KIaKGIrm9lQZGN9SQ42Y/rh+gDxeJzzl87xre/8I/fdex8u9wjRaBSD3sDmW7Zy66Zb2bRhM2JEwdVTHTjOdhP0h6jcXMpdX7iNmq3l80pdu9Rk4sIpKiqir6+Pf371ErZ/fITr/gBv6AvZmnMr64q8aKyF6P1XMPnOIgpqImoLytgoqpi0wBsV1HysMRLR5FOqKoA011KWiu80vsKjDWPrQqJITl8/lt7raH1+IjodfZs24lxTMWXc8WA7t7g/wR1T82rfFkac/QiCwLatdfz6A5+noiz9ZG0LJR4XuXa6g7PvXiHgD1FzawUb76qmaksZOXbpBpIN4Zku2eg7K0mmsoZKDbEw+DrBdRm2e2GdDwrugOovg1KHKEIwDO5Rgfb+KOc7w/QOKxDFPBQUYVGK/OSN99m8th6raOJERIMYlyajRhMhx+pHY3DS5R7gbI+bQW8gmV41R6/htop8KnLNNF+4xvbNaWrJgsCIvoJRlQ17oINK90cMhavoN93Ci6KWz27eTlHbJfLecBBVqXGWVeHJLSSuUhFVa/HmFREypnodmBVGdprv5Gqokwv513j75GnO+e9F0WJha2kev7GtFLsxPdupQqFg66ZbKS4o5vd3/2fi8Th9A320dTpo63Bw4tRHUi3CwhKq1lSz5o5KDGIBg1c8PPu/X8GSZ2LD9iqKKvMx24288tpL7PrKQ2j1aowWAyr1wvyTM6ElpV3mLg0C4SgnfvIT/pM/zg91BRTpitiWfwVr6Cr0g4iCgL4Sj76SKAIaQcVH585y6/p8zsb70Ap6qtIU2gkNeL4s9LgUBAFPcRGe4iI0fj+2jk7KW05RePESnqJiAjYLIZOJsMHAx5e8BLfeSZX7FL9X2szbBaW0hyq52HqRk6c/ZsO6jdy66TZKikooyC/EbrWjUCj4+csv8KXf/O15D2264xQKgfW1lVRvLaftXDeOs9384vCbgOQttWZjCd4RPwFfEL1p9kyhk1noHMy2hr+qBHc0FODyvz1OQDTijtnpje6gJ/JfCHbkMPrBKKG4FxCJE096cgiCAkHpBaWXiBgiHIzx4ZW3yclVIwBKRRy1Io5CESHgDdHfHyE0ltAmR6OkwqzFrFZi1qjQKASIRfAOOPnw1HnW583fVDBIMaa4AYtvkLzBQS59dInf+MwtuAqAsIn4qAKu92HuHi8xlQuEtVrcNiujejNRlCAoEAXYCFQLBp758Brb6i3E9RBzKfjlr9QoYypEFEQUMeLKOApBQCtoMam1qJQa4igJKjV49RaUgoIe1zDPnHgbpUKBRqFGrdViXX8LxjWV+J3DuEdG+ODkB7z3UapHjzAk8PIlAQEBpaDkZEsz773+MXqlAbWgRSEoMduNWHNzyLGZMZj0aLUaBIUCQZBslgqFEkRQKdVUllWhUis59E9HuKXodqKRGCCiUitRqVWoNEqUSiUKpYBCqUAQBARBoHRdIWpN6pTOlOAO+Ub5x4e/hjlu5qDOSpwYa61h2oJaFPG1IAbwCAKj4QC4x6u3vHLqFJrcbehRYhH0eEivbN9C59dCj4uGwnj6h6f9bCi3ALXJgnnEh7Z7EE3XAIlnukvnLvFpVx1DpiJcCh8bRTcblacZNqq4ptLQ3n6By1enVrNpaWmh6b1fojfoMBhzUFsKqCkspdxWiFo1VtotHkMURQRBQKmQfvMfPvMDNq3bQnHB9CkGiqvyKa7KJzgaor9jmIGuYc69d5W+tkG++ej30Ro0FFXmkWM3odaqiMdFoqEokXCUeCyOUqVEa9BgMOvQm3R851uH2VJeh1KlQDFhzinGnjATc6+4Oh+tfvxJN9uCe1WZSgxaI8oJ0WVmvRWLwUpMjOGKzLwSPplgMIhON787bzbaWGntrJSxGJRm9Eoj7lEXFoM17eNMNiPFVfn09PQko9WCwSAmk4nu7sUtnOk0ekRhPEhFq9WmdY4L/S7k42Y/zqLOQyWkr3fOdy7N9zi9WU/ZOil4r6enh9bW1uT5qVSqRc+/yawqwS0jIyMjs9rcAWVkZGRkZMEtIyMjs9qQBbeMjIzMKmNVeJU0NTXhcDiw2+04HI4ZixDLzE5LSwtWqxW73c7JkyfZvn07Vqt1xu9X/t5vnO/g6NGjOJ1OWltbcblcHDoklSG7EX/7vXv33tDnB4C4Cti5c2fy/wMHDoiNjY3LOJrVy65du0RAtFqt4oEDB5LbZ/p+5e/9xvgOWltbU37vXbt2Jd/faL/9gQMHUsZ+o51fghVvKmlqasJqtSbf19bW8uyzzy7fgFYx9fX1iKLIyMhIiuYx3fcrf+83ztxzuVwp477jjjs4duzYDffbOxyOlHHfaOc3kRUvuFtaWrDbxwMK7HY7LS1LW4HjRsLhcKR8fzN9v/L3fuPMvdraWpqbx8vmtba2Ultbe8P99k1NTezcuTP5/kY7v4mseME9PDx9NJfMwnC5XFRXV9PQ0IDL5Zrx+5W/9xvzO3C5XDQ1NbF///4b6rdvamrioYceStl2I53fZFa84M7Nnb1Ku0z67Nmzh9raWqxWK/X19ezbt2/G71f+3m/M7+CRRx7h2LFjWK3WG+q3d7lcKeYPmPk8VuP5TWbFC+7a2lqczvEk7U6nk9ra2YsAy0ylqamJ+vr6lG2J73K671f+3m+8uXfw4EEOHDhAdXU1DofjhvntDx48iNPp5OjRoxw9ehSHw8HRo0dvmPObjhUvuHfu3InD4Ui+b2lpYffu3cs4otVJdXU1e/fuTb5vbm5m9+7dM36/8vd+Y829o0ePsnPnTqqrq4Fxe/CN8Ns//vjj7Nmzh127diXPMfH/jXB+07EqcpU0NTXR0tKS1BRWnc/lCiHhywvSo+VEz5Lpvl/5e78xvgOHw0FNTU3KtgMHDvD444/fUL+9w+HgwIEDNDU1ceDAAXbt2nVDnd9EVoXglpGRkZEZZ8WbSmRkZGRkUpEFt4yMjMwqQxbcMjIyMqsMWXDLyMjIrDJkwS0jIyOzypAFt4yMjMwqQxbcMjIyMqsMWXDLyMjIrDJkwS0jIyOzypAFt4yMjMwqQxbcMjIyMqsMWXDLyMjIrDJkwS0jIyOzypAFt4yMjMwqQ7XcA5gPmzdvnpJXeL709PRQWlq67G3IY8leG5PbaW1t5fz584tqr6ysbFFVUhZzXov9Tm7GvlfSuDMx/6YgriK++MUv3jBtZKqdG20s2TifTLRZWFi4qOMXM4bFjv9m7HsljTtTc3oiq0rj7unp4cEHH+Thhx/m4YcfXlAbCz0u021kikyN5Ub8Xp555hmeeeYZenp6Ft1eMBjkwQcfTGl/Pue7mO9msd/rcv4uy3Xey3nOa9asSZkrmZh/U8j4rSCLZOPOtZzcaOeTCbLxnawEjXu1crPO0Uyedza+Q3lxchlZSRrqSmGlfieZsLuvRlbq75FtVvp5y4J7MktYgnOlT47lYKV+JytacIvxrDW9Un+PbLPSz1sW3BPxdsL38qD9peUeiYxMevS+C/9ig4GTyz0SmSVEFtwT6X0HQk5o+/lyj0RGJj36PoCwBzpeXu6RyCwhsuCeiP+69Br2LO84ZGTSJTgkvfq6l3ccMkuKLLgnEnRKryHn8o5DRiZdQi7pNTi8rMOQWVpkwT2R0MjYq2tZhyEjkzaJp8PE3JW5KViVATgJFhOIMy1h99irN3Ntyiw5icCbBFkJgFgpRHxjr/KcvZlYVYK7tLSUF198MXsdJLQX+SJY1Uy+oU+82S+UTETtZoWkxu1e3nHITCGTkbuTWVWCO+tEEoLbt7zjWOXYbDb27NnDgQMHUrYLgsCBAwd4/PHHOXjwIMeOHQPg5MmTVFdXY7fbATh27BiCILBz587ksdXV1Rw6dGjpTmISWVcaFkpCyYjcvAvqBw8e5NChQ1itVux2e3JezTSHBEFAXIJ4jcRNPhOKw2RkwT2RsAcEJUT9UiCOICz3iFY0e/fu5eRJyX+4sbGR6upqAOx2Oy0tLSn7Hj16NCW73uOPP87jjz8OQF1dHUeOHJmSfS9xAcrMQtgNCs1Na95raWnh2LFjtLa2Jt9P5EadQ/Li5ETCHtDlSpFoseByj2ZFc/ToURoaGmhubmb//v3s3bs35fPa2tqUi+jQoUPs3r17qYd54xNygaEIYgGI3nxz1uFwpLxfSOpdl8s1r+0rAVlwTyTsBl3e2P83pwaTLrt27Uo+htbW1uJ0prpQ7t69O2nacLlcWK1WrFbrvPpoaGhI/k3WpGSQFIywG0xl0vuET/dNxK5du3A6ndTU1LBv374pgjydOfTII49MK6Rn2r4SkE0lCURREtb2rdL7iA8oWNYhrRYOHDjA/v37U7bV1tYmzSiHDx9m9+7dU4T7XDQ2NmZsjDckQSeIMciphr73YbR/XIjfRDQ3N9PU1ERjYyN1dXW0tbUllYTZ5tDRo0d56qmncDgc7NixI2kfn2n7SiIrGvfkR+L29vZsdJNZogEQo6Af07hlz5K0OHjwIPX19ezatWvKZ7t376apqYljx45N+7nMIklE+lrWSa+jfcs3lmVm586dHDp0iIceeojnnnsurWN27dpFc3MzO3fu5Pjx40nhPNP2lURWBPfkR5KGhoZsdJNZEj7cCVOJ7FkyJ/v27WPnzp0zCuWEZ0li0VImw4yOCe6cse/3JoyebGpqSpE3DoeD7du3L+OIloYlsXFnyvUm4Us7MbgiYyQEtz5fepU17lk5fPgwhw8f5pFHHqGuro76+vop+yTs2pMXLpeKZ555hgcffPDGDcBJCGpdHii1EHYt63CWA7vdzr59+6irq6Ompob6+vp5L1AeOXJk2vWXmbavBJbExi1kyK0uq760iUAG3ZjglhcnZ2XPnj3s2bNn2s8SrlmQamOcaf/m5uYp2zJxs8+mH+2KIOQChUoS2irjTTlna2trZzRlzDSHJm+fSTivVKENWRLcw8PDPP/888kvyOVypbwH+PKXv5yNrhdOQsNOepW4lm0oMiuPFRk5GfFKAlsQQKWXzXsrjFUXOVlVVcVjjz2Wsm3ie0EQVqDgHpv0aiOoTXKiKZkUVmTkZHRUEtgASp0UOCazYlh1kZPTPfqueBKPmSo9aCw35UKPzCojMgoKrfS/UgcRWXDfLGTVxu3xeGhqasLhcFBTU8POnTsxm83Z7HLhRP0gKKTwYY0FAoPLPSIZmdmJBUCpkf5XaqT3MjcFWfMqOXLkCJWVlTz55JN8/PHHfOMb38BqtfJ3f/d32epycUT8oNRL9kJZcMusBqIBSdEAaYEyOrq845FZMrKicf/0pz/l0KFDNDc3U1VVldzucDh46KGHqK6u5nd+53ey0fXCifpBpZP+1+TclOHDMquMWBAUaul/WXDfVGRF43766ad5/fXXU4Q2SGkVm5qaePLJJ7PR7eKIjEp2QpAWJ8NyfmOZFU4sNG4qUWikOSxzU5AVwS2KIjk5OdN+ZrVaF+yjm9UAnOiopLUAqAxyweAbgBs+ACc6SeOWM1reNGTFVJKbm7uoz2ciqy5ZEwW37Fp1Q3DDB+DEghMWJ7WSzVvmpiArgru1tXVKwM1EJqdeXBFEAxM0btm1SiaVFRmAEwtJcQcgmUpkwb2iWHUBOBaLhW984xuzfr7iiAVSHzvjYYjHQKFc3nHJrAhWZABOLAhaq/S/bCpZccgBOEtBNDghmGHsNRYEhXH5xiQjMxux0AR3QFnjvpnImh93e3s73/3ud+no6MhWF5klOpq6Qg/yhSCzsolOsHEr1JIgl7kpyIrgPn78OLW1tTz33HPU1tbyxhtvZKObzDJ5oQfkC0FmZROfpHHHQ1IlJ5kbnqyYSp544olk8I3D4WD37t2cOHEiG11ljtgEU0nC1i0LbpmVTDSQ+pQoxiEeBaV6ecclk3Wy5sedCL6prq7OWCGFrBKdlPcB5NwPMiubiTbuxGtcVjZuBrIiuCf7aWeqkEJ2A3AmBDMkNe5w5vuRWTJu+ACc2DRzNip7ltwMLIkft8vl4rvf/W7KPn/0R38073az6pIVC47bthPai2wqWdXc0AE48RjEIxOCxmTz3s3EkvhxWywWvvOd7yTfC4KwIMGdVWITMq3JNm6ZlU7CZ1s5ac7G5afEmwHZjztBSsIeWXDLpJIw0yVY9gjKhOCesqAuC+7lJhExmWDVRE6uOkRxTHBPugjkhR6ZMVZc5GTClj15zsrKxrIz+aaeDVNd1gJwVhWJyS6bSmRWC0mNe9LipGwquSmQBTeMR0hOMZXIF4HMCiU5Z8c0bmHs4TkeWZ7xyCwpsuCGCdqLrHHLrBJm0rhlZeOmQBbcMFV7USilwsGy4JZZqUwR3GNZLOU5e1OQtZD32YJucnNz+dM//dN5t5u1lf3YJFMJSNq3fBGsSpZiVX/ZSczNKQvqsqnkZiArgrumpib5f2NjI9XV1cltx44dS/l8PmRtZT+pcevGtynUslfJKmUpVvWXnehkjTth45ZNJTcDWRHcjzzySPL/o0ePpgTfPPbYY+zevTsb3S6chOBWTNC4lbLGLbOCSWrcY3M2sTgp27hvCrJu4z5x4gRerzf53u1209TUlO1u58dkGzeMlYKS8z7IrFBkd8CbmqwH4DzxxBNUVlayZ88eQNLAJ2rkK4LYDIJb1rhlxlhxNSenxB6MLU7KNu4Vw6qrOTmRxx9/nJ07d3L8+HEAnnvuObZt25btbudHdFR6nSi4lRo5ratMkhUXORkLgqAcF9iCUvqTTSUrhlVXc3IytbW1uFwuHnjgATweDx6Ph5ycnKXoek7euzzAT7//Nt8sZ6pXiWwqkVmpxELj5pEEcvmym4as27h/9rOf8bnPfY69e/cCMDw8TENDQ7a7TZtDr12GaIBIXIkvFB//QKmRq2bLrFwmFlFIoFDLNu6bhKwL7n379vHaa68lK+JUVVXhdDqz3W1aiKLIu5cH2FKsIySqON0+Mv6hQq6aLbOCmZjNMoFCLZtKbhKyLrhFUcTj8SQDctxu94JLmWW6Ak738CiDnhDrC1RERDXnuyYIbqUsuFc7N3QFnGlNJSo59uAmIes27qeffpra2loEQWD//v0cPXqUAwcOLKitTC8Qne6QNP98A4guDVd6PeMfKrUQ9WesL5ml54augBMLTmMqkddlbhayLrh37dpFdXU1zz33HHa7PcVsstyc6RjBatRgVEURFRo6h0fHP1RqIeRatrHJyMxKLDS1mrscNHbTkHVTyenTp6mtreXpp5+murqaw4cP097enu1u0+Js5whV+SYUomQv9AUj+ENR6UOldtxNMEE0uvSDlJGZjmkXJ2XBfbOQdcGd8CA5fvw4e/bswW63Jz1MlpvzXW4qC0wo4iGEsYtg2DuhJNREwf3D/wOfKYQux9IPVEZmMjO6A8rrMjcDS5bW9dChQ+zfv5/HHnsMh2P5hV8gHKVjyEdFnhFFPIwwtkLv9I2tyiu1qe6AL/wruJ3w6rMztjkaC/K/rv2A3zt7gFHZlfCGItML44smFhpPLJVAXlBfUWRzcTzrgttisbB//36amprYtWsXwIK9SjJJa58XUYSyXAOKeBBBJWkv7tGE4J5wEUSj0H5Z+v/912Zs8y+u/Rt/5fgR/957nMa+t7M5fJklJrEwviLC3UHyHpmscSu1suBeQTz88MO8+OKLlJaWZrztrAvuxsZG7HY7jY2NVFZW0tbWtmCvkkxyrU9KfFVmNyDEwwgKDSqFMEFwT9C4e9ogHIK7HoCTb8E//vmU9mJijH+7foyHCj/NWn0J77rOL9WpyNyMRIOZE9zx+Nz7yKwosuJVMjGkvaqqisceeyz5WVVV1YrwKrnW78WkU5GjV6OIBxEVaoxaFd7AhMXJWAjEOLRdkrZ95RGw5MKRb8CnfxNuvzvZ3kn3FYYjHu61bsYV9XHa27oMZyVz0xCbQXDP1xPqH/8cnv8XeP4cWHMzNjyZ7JIVjdtqtXL69GmpA4UCpVKZ/Eu8XwiZtDM6+r2U2AwIgoAiHkYU1Bi0qnGNO7FiHw1KZhKtHnIL4aFHwZoHTT9Nae/tkbPoFVo2Giuo1hdzwddBTIwtepwyC+fGD8CZ7FUyT4171C8pIUN98M4vMzs+maySFcE9MjLC7bffDkA8HicWiyX/Eu8XQibtjG0DPoqsUsUbRTxAXKHBqFXhCY6lxUxkCowFJI27uAIEARQK2LhNMplM4KTnKusMpagUSmoMJYzGQ1wdvQEFxioimzbGZWc6P26VbqoL62y8+qw0p1VqOH8is+OTySpZEdwWiyUbzWaUjkE/hRY9AIpYCFFQo5+ocScEdzQArRegqGz84Mr1cPVsil/3OV8bVfpCADYYyhAQ+MB1cUnOReYmZFp3wHl6lRw9AlvugNvvgcufZHZ8MlklK4J7onlEoVAk/xZrKskUsXic3pFR8nMSGncQUVBh0ijx+Cdr3GOmkqKK8QZKq6TFym7JrTEuxmkN9FKmywfApNJTYyjmnZFzS3ZOMjcZsfAMi5Nppmno74EzH8KndkJZFVw5AyvA20smPbIiuCeaR+LxePJvsaaSTDHkCRGNi+TlSMJZEQ8RV2gw6FS4A2OCO2E/dPZJ/tvFEwR3Ubn02nkVgJ7QEKF4hFJtXnKXTcYKPnZfyvq5yNykxKexcScW1NPhrV9IRRi23gUV68AzAtc7Mj9OmaywZAE4iQIKib/lpNclPU7aTRMFtxaTVo0nEJH8zBMa95hwpmCCndSaB2oNdEmeI9dGrwNQphsX3FX6Yi6PdhONywuUMllgOlNJIvYgHc359edh/a1gyoHqTdK2lncyP06ZrJB1wX38+HHsdjtWqxWbzYbVaqWuri7b3c7KoEfyz7YZJY1FEZds3EadikgsLuUrSQjuRIh7fvF4AwoF5BVBdxsAl/xdqAQFRRpbcpdSbS5RMUZncCD7JySTddL1aBJFkdhS+EVPa+PWAuLcxRT8Xvjoddh2r/Q+xwpr1s0aXCYzf1Z15OSjjz7K8ePHicfjfOUrX2FkZGTZa04mBLfFoAZRRBEPEldoMOulC8HpC40L7t4uMOaAwZTaSG6RFJgDnPG2UaErQD0hBLlEK/nEtgf6snw2MktBuh5ND/zla/zet97N/oBi4WlC3hML6nOkW2h5F6IR2Hrn+LaNtfDxG5kd403Oqo6cFEUxKagFQcBisdDW1pbtbmdl2BtCr1GiUSkRxAgCcUSFhpyk4A5PENw9qdp2gryi5OLk+67zrDeUpXxcoLEiINAmC+6bBpc/zEnHMC81dxOOZtlENq2Ne+z9XImmzn4EZisUTvSU2gADPeAczOgwZbJD1gV3bW0tr7/+OgDbt2/nc5/7HC6Xa0FtZSoAZ9gXSgppRVzSTuKCBsvYtkFPcFxw9/VCbsHURgpKoLsNV9jLWV87t5qrUz5WK1Tkayw4ZMG9bCx1AE7HoC/5v6PfN8ueiyQelSJ6pzWVMLdL4KXTULFW8uFOUFopvbbJC+qrgawXUnjuueeS/z/22GNUV1ezc+fOBbWVqQo4Tl+YHMOY4B7LRxJXaFGrFJh0aklwJ7SZ/n7YtnlqI/klEPBxvvMEIiKbjBVTdinX5XPZ37Xo8cosjKWugNMzMh780j7o45bSLMUzJDxHpiukAHObSq6dk4LIJlJQKgnyzmtQd39mximTNZbcq6S+vj5Zf3K5cPpCmHUJjVvSTsQxQW0zaBjwBKVJLGhhyDm9qaRQslsNXjuJSlCkuAImqNIX0ey5mqWzkFlpDHrGXfG6h+cRwThfEoJ7Oq8SmN1UEgpKazMllanb1Rqw5Sc9pWRWNlkX3E888QQKhYKqqioqKyuprKxc9iRTw95QciFSmTCVjAluq0lD38jYxA9opcxp+SVTGxlzD4y0X6RAY0OlmBpUdJupmvZgP61j7oIyNzZDniBmvZpCi45uZxbrlSYF9zS5SmB2U0nnNWlOF6U+If5L9y+5ZIjyi1M/5r9d+vaKSL0sMzNZF9xHjhxhZGSE4eFhnE4nTqeT4eHhbHc7K07vBI17TDuJC9Kktxs19LvHHjV9Y8J4Oo1bo4XcQjRdbeRrpn8krs1ZiwKB152nMzp+mZWJazRMjl5NrlnLdWcW82In0g1PF4ADs+crSeSVTwSRAV2BQX7Ye5xhq4n1g6P8Q+fzNDlbMjhgmUyTdcG9Y8eOFZe7xOkPjdu4J5lKLHoNTl9I0ji8Y19PbuH0DRWUYrneS646Z9qPDUoda/SFfOJd/oo/MtnHPRrBqFVhN2npdS2BqWSKjTsNjbvjChjNYB6/Jt9znUOjUFGxppbqQR9VukL+redYhgctk0myvji5d+9ecnNz2blzJ3a7Pbn929/+dra7nhZRFBnxhZNeJcqExq2Q8pZYjGrCsTieQASLFzDqJPvfdOQXU+j4AJvKNP3nSIE4V0a7M3oOMktPwqMpQWLhcyJuf1gS3EYNV3qzGB2c1Li1qduV0hyeVePuvDa+EDlGi+ca1fpiIkWlqEb9/JZYzo+GTxAX4yiEJVsGu2F45plnUjzfsuHVlHXB/eijj7Jjxw62b9+e7a7SwheMEorGpeAbxjXuhI07YUJxekNYfHEwzSC0AfKLKTnhxaY2z7hLkcbOGZ+sca920vFoco+G0WuV2EzaZFqFrJDwGpmS1jUNjburNcX0F4vHOOtr4wH7NoJ6ad3ms4MC38rxct7XwVbz8hc9WW1Mvqlnw6sp64LbYrGkuAQuN0Ne6TEzRy8JZGVM0k7iY9pLMgjHH6bKGwODeppWJOL2AiyBCIXhmbWSfI2F66HltenLLA2eQASrUUOuSYt7NEIwHEOnyUImzJk0bkE5ltp1Fo27p11K4zrGpdEuAvEw6wwlhLW5xLQ6NvV6UFkUvOs6JwvuFUrWn4N2797N888/n5HEUonH1cTfQgJx+t2SNmIb06QVsQAiAqIgCWiTXrqXjfhD4I2CfuavyGeVbNsVIzNrOHkaC75YEE+66TZlFk0i8Cbxt1QBON5ABINGhd0sza0+d5a07oTgVk7zNKjUQWSG4J9YTIqOnBBQ1uy+il6hoVyXDwoFwaIScjo7WWco4105LfGKJesa91NPPZWMlBQEAVEUEQRhQaldMxGA05/IDGiUtBVlbFSyb4/Z/NRKBXq1Epc/DJ4wVMxsvx7JMZIDlIzMrOHYVZIZpS80Qo7KuKixy6THUjyqTocnEMGgVZI7lnWyzxWgMn/m+bNgEqaQ6QS3SjdzTu7hfohFJX/tMT5yX2adoRTlmC07lF+AobuDzabtvD1yNtMjl8kQWde4m5ubV1Q+7u7hUbRqBeYxzVoRHyWu1KfsY9Kp8Y54IRQF7cxjHdariAmQ7/TOuI9dIwnu3pAzA6OXWcn4Q1F0GlUyXXDvLE9iiyIhuCebSgCUeojMoEgMjD152CWN2xXxccHfwUbjmuQu4dx89D3d3Gaqpjs0RNtobyZHLpMhsi64H3roIbzemQXbUtM1LJUsS0RvKmP+pH07gUmnIj7QL73RRic3kcQlBhg0asgdHplxn6TGHZYF942MKIr4g1EMGiUmnQqtWsF1Z5ZcAhM27MmLkzC7xj0wFghmk6J833ddAGCzaYLgtueiGRmmVluGAoEm56mMDVsmc2RdcD/99NPs2rWLjo6VUV2jfdBHgUWXfJ80lUzApFMhDI1lSdPMnNt4JOJj0KzBNDA04z5GpQ6dQi1r3Dc4oUicaFxEp1YiCAJ5Zl1K7pKMEhvLpSNMs/Cp1M68ODnQAyoVmCQf7jedZ6jRF5OjGn/iDNskl91cp5tNxjW8OnQy48OXWTxZF9x79+7l2LFjVFVVpdSeXC7a+n0Upgju6TVutXNMGGtmLgU1EvXhzDGi75158UsQBHLVObLgvsHxhaSSd3qtZILLz9HSk618JdHAeLDNZJRaiMykcfdI1ZsUCrzRUZq9V7h1ktdI2Cblkdf297Hdso5jzha5itMKJOuC+9q1ayl1J5fTxi2KIh1Dfoqs4xqGJLgna9xqDO5hMOlBiEtJ66fBGfHgtpjQ9c2eiyRXnSO7BN7g+IOSSU2nli6pPLOOjqEseRLNJrgVs2jc/WOCG3hv5DxxUeQ2c03KLhGbVMVJ19fL9pz1eKKjNHuuZGzoMplhScKi9u/fT25uLrm5ufzxH//xUnQ5LU5fGH8omqpxR33EJglus05NjncEMWfMIyA2fZrM4bAXr82CZngQITyzSSVPY6ErKCeoX83MlQveH5IEt0EjadwFFh2dQ1nKyR0dnVlwq2YT3N1glTTqD9wXWaMrwKIypOwiqtSELVZ0fdfZYChHp9DI3iULZFWXLvva176GKIo4HA5aW1uJx+N87Wtfy3a309IztliUnzPJVDLJq8SsV1MQdBMzjQnuGSLRhiJugnY7giiiG5h59T1fbaFLrj25qpmrdJkvoXGPBdwUW/UMekJJgZ5RooGpCaYSKHUzm0p6O8GejyiKnPJeY72xbNrdwvY8dH3XUSmUbDSWJxcxZebHqi5ddvLkSZ5++mksFgtWq5XvfOc7NDU1LaitxVbASQjuPPNEwe2bYuM261UUhNyETWOh7DNoMIMRN7GxBFS66zPfVYu0drqCQ8RE2Va4lCxlBZyEqUSvHhPcNkkZcPTP7FEliiK+YGT+ncWCs9u4p1M0RFHSuG35dAUH8ERHqdYXTdtEODcPfU8nABuNa/jQfXFZ07wGAgF+9atf8eKLLxIMzlEk4iYh6wE4oiji9XoxmyUh6PF4FjwJFhuA0zsyikIAq3FcW1HERokrJmncOhW5IS9evRkDTFtRJBgL44mOIuQXIArCrHbuYq2diBilOzjEGv0MmQaXAVEUGR0dxWAwLHtxi2ywlBVwEouTCY27zC4FW1257mFrhS25XygS4z/84zucaB3CqFXR7RzliS9t4etfvjX9zmZdnNRNr2i4hiEUgNxCLo5VZVqjm34uhvLyyf3gHQA2Gsv5cd/rdAUHqdBPU8Ivi3i9XhobG3nllVfw+/0IgkBTUxMHDx5Ep9PN3cANTNYF9xNPPMEDDzzA7t27ATh06BAHDx7MdrfT0usKYDdpUSrGhdR0i5PWaACNGMOrNZELEJt6IQyEJd9tm85KxGpDNzBzbcmyseo4V0a7V4zg/uSTT/jHf/xHenp6sFqtiKJIKBTirrvu4pFHHiEvb2pFH5mZGV+clAR3jkGN3aThfLeLrzDuJ/3t1y7TdLaX376jnHA0zp1r83jqhXPYTFq+9rkN6XUWC0wfNQljGvc0ppLrY+64uYVcHr1IgdqKfoY2QgVFaNwuVF4Pm8Z8vN93nV9SwX39+nWeeOIJ3G43d955J/feey9+v5/Dhw/z1FNP8Rd/8ReoVFkXXyuWrJ/5rl27qK6uTiaaeu6555JV35eaPlcgmaMkgTI2SkyZKrh1bsl1z6kxUwnTPnr2jQluuzqHsNWGtn9mwV2ktaEWVFzyd1GfW7eoc1gsHR0dPPPMM7z11ltUVVWxe/duhoaGEAQBlUrF+++/zx//8R/zl3/5l2zcuHFZx7qa8IeiKBUCGtW49bGm0EyLY9ybKBiO8Y+/vET9rcX84QPrUo5//N+bsRjUfPW+1KLT0zKXjXs6jft6u/SaW8Tl7ibK9DPfmIOFUvZAQ4cD+5bbKdfm847rHL9b/Nm5x5YBXC4XX//61xFFkf/+3/87tjFPF5vNxle/+lV+8IMf8Nd//dd8/etfR6ud4cnjBmdJblm1tbXU1tYuRVezMuAOYjVMmPBiHGU8gDhJ41aPSBdbrzi2OBmZKrh7Q06UggKLykDEYkM72D9jv0pBSbkun4u+zsWfxCI4ffo0f/EXf4HJZOJLX/oSd911FwpF6jLHnXfeyQ9+8AO+/vWv881vfnPZy8ytFnzBKHqNMsXktKnMytEPOwhFYsTiIt/8xQWGvUG+fOealGP/4LNrcY2GeeyHzfzGtjJsxhmEcoJZvUr00uKkKKZWce9pA52BmMlEa6CHz82iQIQKihAFBca2VjxbbmeLuZK3nGfm/A4yQTgc5i//8i/x+/388R//cVJoJ7jlllv4/d//fX70ox/xZ3/2Z/zN3/zNTWk2ycri5Pbt27njjjtm/Lvzzjuz0e2c9LuDKfbtRErXye6AGucQMUFBV0wLgmLa4qvXg8PYVTkoBQURixXt4OxeIxW6As77ly961O/3c+DAASorK/kf/+N/cPfdd08R2gBGo5E/+IM/wGaz8Vd/9VfyYlCa+IMR9JNSuN6zIR9/KMpDf/8WFV87yoGfn2P3PZWU5aa64AmCwB98di3BSIzvvZ5GcemIf7xowmRURkCcmiGwqxXyS+gKDRGMRyjT5U97OICoVhPKL8DYLhUOvt1cw3l/B0Nh99xjWyT/+q//ytWrV/n93//9lMIrE7nlllv4wz/8Q65cucI3v/nNm7I+ZlY07iNHjky73el0cuDAgQV7lSyWAXeQ6sLxbG3KmDS5Jy9OaoYHGDWYpdzdOdPnN+4ODZI3lkAqYrWhGZ5LcOfzq+HFhw+Losi1051cPdWB3qhlfV0lJTUFcy4uNjY24vP5ePTRR1GrZ84xDqDVavnqV7/KP/zDP/CDH/yAPXv2LHrcNzq+UDRp305QnmvkK3dV8Pq5Pr5y1xru2ZBPdeH0RTdsRi2f2VTIoWNX+JPfuAWNapbo4tk0bvXY/A65QDOhr86rkF/MZb9UjalcO7PgBslcYmiXCoDcPhak89bIGb5SeP+sxy2G1tZWXnjhBb7whS9QUVEx676VlZU0NDTwox/9iPvvv5/778/euFYiWRHck23YHo+HJ598ksOHD7N7924aGxuz0e2cDHmDWCaYSpRRyVUrNsmPWzPYR8hsxROMElfoUEwjuLuCg1SOuVNFciyoAgGUfh8x4/RpPMt1BQyEXbgiPqzqhaX6jEVjvPK9tznz9mUseWZCgTDvv3QKQ46eoso8ytYVcdunN5CTm9r+4OAgP/vZz7jvvvuwWq1p9ZWfn8+OHTt44YUX+PznPz/nhXSzI5lKpl5Of/jAuin27Jn4nTsrOHa2l++/2cqenetn3jHik0wi06EZq38aHALzeEFg2q/Atvu45O+iSGObcWEyQaiwCOvpZgAKNFbKtfk0DZ/KquD+/ve/T15eHvfee29a+9966620tLRw5MgR7rrrLjSaOUxMNxBZ9eP2eDw8+uijVFVVIQgCbW1tfPvb316W4sH+UJRAOIZtguBWjQnuuDL10VXb200sV1q8iaKeonFH4zGuh5zJ6u4Rq2SHm83OvWZsRf7ymCvWfPGO+Pn3J1/i3HtXueeLt/PFPZ/hy39Sz86v3k31llJCo2E++MUp/vl//JjjP/6AgH88x8qhQ4fQ6XR85jOfmVef999/P1ardcYnqJuJuWIIfIHIoqvdrMk3sXNrMX/z0zNSPviZiPik9K3TobVKr4EJkbqjfujrguIKzvvapaIJcxAsKEI7NIAiIJkJt1vW8/LQR1kzS1y9epUTJ06wY8eOeeUy+sIXvsDQ0NCyKYOzseoiJycKbJvNhsPh4Kmnnlq0wF5MAM6AW7LVpti4o5LNLqYc11CFcAjtQC8Ul6AEQqJ6ilfJ9dAwUTFGoUYS2OOCe2ZzSbk2HwGBc772eY+9t22Qf/nzowz3uqj/D3dTvVXSpBQKgaLKPG7/zEY++9CdfPn/qWfz3es42XSOf/7vP6K56TwtLS28++67fOELX5j3Io5KpeLzn/88J06coLm5ed7jXm4yeeHMGTkZiiaDbxbDf/q1GkbDUQ4duzz9DqIIYQ+oZyjKMaZMEJyQsbJdaitYWEzraG9aLqmhfGkfQ7e0LnOvdRNdwUFavGnY4BfAj3/8Y/Ly8rjtttvmdVxBQQGf+cxn+PGPf8yZM0uzgJouqy5y0maz0dbWRmNjI7t376atrY3Tp0+n/C2EuS6e2RgYKyM1UXCro5JLX3SC4NZ3dyCIIuHCEixGDYGYckpi+s6gpFknBbfFCoB2YGaNW6eUykOd8bXNb9xdTn705EvoTTq+8F/uJ79s+gUbALVWxa33r+dLjz5A2foifvmvb/N3T/89lZWV3H777fPqN8Gtt95KdXU13/rWtwgEslgANwtk88KZjDcDGjeA/f9t783D2yrvvO/v0b6vXuMltpw9IYvjhC0ECA6FltJp68DQPrP0bUlon+edd95OCWSmc80802kT80LnKV0ggbYDpaSJQ9gDxEoCISFksQnZN8t2vC/abe065/3j6GixdGRJlmzLPp/ryqXo6Og+t+T7/PS7f6tCjA1LS/Dy4RsIBMn4EwJugPRFbNlj4UsAQgD4bJFjpssAgKsqAkGQqGbJmIzGW0jvECW9tE18pbIGGoECe/uPpvV5UuHatWs4ceIE7r333oQO8/Gor6+HwWDAz372MwwNzY6aQDkR3CtXroTZbMbWrVvxgx/8IO7f448/novLJmUgpHHrouK4hX4LSEIIKirlXdLdAYrgwa8vhE4ugtPHj0toaHcPQM6TQBUysVACIXwqTdJ6JQAwTzYHZ+ypV1pzOT3Y+9wByFQSbPjrWyFVpKYxSxUS3PbVFShdrYDZMYQiYQVIMrMtLkEQ+Pa3vw2LxYJf/epXs9KDnwoOtx9ycXZcRl9ZWYZ+mxtHLibIDWCKlYlZdq8EQWvj3qgIENMlQFeELwMDkPLEKBWx//gzBOUKBGRyyLpojZtP8HGHZgneHDye7sdJfp1gEL/97W9RUlKScX4Hn8/HY489Bh6Ph9/85jdZnd90JSfOyem4re63uSHgE1BKIxEVQp8ZAYE6Jt5V0nMTfq0elEAIvVIMu4MPyj+K6JiNDvcASsTamDBZn04PSZK63ACwQFaGV3uNCFJB8BMVwY8iGAjizV83w+Py4cHv3QWROHkkyFgoisLpy8dRrC+Bdxg49uYZrP/WGhC89FPbCwoK0NDQgNdffx1LliyZtB6O+YTT7YcsS4J7fokSlQVy7D7ejo3L58S+6GI6MyURvkIFEB2613YJKK3ElyPtqJaWgBdauJ5RCv2dFIZ6STjMgDcU+SkUAWIpcL7qf0J4XgTv+zehqFDg9oJFODB8Cu2uPlTLSrPyWV999VVcu3YNTzzxxITq9Mvlcnz961/Ha6+9hlOnTk1ZyPFkMSllXacDvVYXCpSS8KIFAKF/mBbcUYj7uuHT0Y7JAqUEroAAfm9soSCTuw/F4tjEAJ9eD2lPcsfjfFk5XKQXN1zJ63eTJIl3dx3Bzat9WP/N1VCoZUnPT8SFy+fQ0dWBdXfcheXrF6L7+gBOvHc27XEYVqxYgTvvvBMvvfQSOjo6Mh5npuJ0x8dxZwpBELjvllK8c6YLlpExjTxGQ2snmeAWyGg7OMONCwiWVuCisx3Vviq0XyRx4kAAh/YGcekUCc8ooCsB5i4kMHchgaJyAjIF4BfJYBsVo//YAK68fBXEMyP4m0O34/3PPgOV4Q4ums8//xx79uzBAw88gKqqqgmPt2zZMtTU1ODFF1+EL0mZ5ZnArBHc3WYX9MrY2FehbxgBgSrmmKSvC3497XUvVEngJoWgvJFkhiAZRLdnCKVjbhxvYTFkN5Mn2DC2xYsj7Of5vX7s/3UzLn3ehju/vgrFc9OvGRIMBrDvnSaUz6lAZflcFJZpsfSO+TCd70Lr4cxLdD744IPQarX45S9/OaUNn6cjTrcfiixp3ACw8ZZSkCTwhyM3Yl9w9dEty5jokUQIZBGN2+3CSP8w3gnegltPPgLfJ/Nx+QyJYBAwLCOwegOBJWt5qFzAQ8lcAiVzCZTPI1C9lIflRT24v/svWPR3NVjwd/NQuq4EJaMa2P97CC/9SxOutrRnbDrzeDx4/vnnsXjxYtx9990ZjTEWgiDw8MMPo7+/H2+++WZWxpyuzBrB3TE0giLVmAxJ3yACQk34OX/UCaHDBp+edsyIhTxQfGk4UQcAerzD8EdFlDB4SssgdNggtLB3utEI5FDypbjCEhJoG3Tglf94C21nb2L9t+owd8mchOeNxwfGAxgYGsC96zaEE3NKqwqwsK4al07cwIXPMosMEAqFaGhowPXr12f8jZEOHl8Q3gAJuSR7glsjF+H+5aV4/sBlWKNDA0d7AYmezuhlQygHfHZQJIUzfzmEFxb8G67YCuFQD2JhHYE19xFYXMdDUTkBgZDddOYpmQO+zwdpVydEKhH0K3SwfSWIw7ddAY9HYN9/fYQ//tt+mM6nH+L6zjvvwG634+GHH85qZcqSkhLccccd2L17N2w2W9bGnW7MGsFtGnCGayQziLwDCAg04eeS7g4AkVAoAOCLFBDADwTpLWuHmw75KxkjuN3ldIKK6spF1jkQBIEKSRGuubrjXuu83Ivf/+sbcDk8+Mrf3omKBeN7/hNxo/063m9+F2tX34rCgthqbpWLSlGzvAJnj1zG+WOZtaOaO3cu1q1bh//+7//GjRs3xn/DLMDuogVrtpyTDI+tq4bHH8T2N6PC3Fx9tOBOhlCBvn4Cr/zHW/joUA8qXTfQt7gVIoMV2gICPH5qgtJbXIKgRAr1xS/Dx5YoK9GuGkTRwyW477Hb4PcGsLvxfby+4z0M91pTG9frxRtvvIG6ujrWtPaJsGHDBgDAnj17sj72dGFWCG7LiBdDDi8qomtEUBREvgH4hZGFI2u/AVIogl8bMU+IZaEkGxetSXe4+6HgS6Ack7Tj0xXAr9JAfa416VzKxPpw2jHD9S86sbvxPWgKlXjg79dBW5xZvLvL7cLLr+5EafEc3F53R8Jzqm8pR82KSnz5yRWc+vA8yEQhZ+PwwAMPoLi4GP/5n/+JkZEctefKIxiNOJsaN0CHBv71HVXYZbyOG/0hm7WrHxBpWN9jt/rx9qnb8IcjX8WI3YWNiquoFVxFF4ZgkKbpUOTxMFI9D5ovTtPx46BreGsFChgtX6C0uhBf+ds7cfe36zDcY8XL/9yEY2+1IBhIbkY7ePAgHA5H1kwkY5HL5Vi3bh3ee++9GRsemFeCO9MEnHOdtCZQXRSp3SD0D4NH+eEXRrQXedsVeErLgahYUpWKFuKdobCoDk8/SsQ6xO3uCAIj8xZA03oy6VwqJIW4OtoVtg22X+jGG7/6CHMMRbj30Vshlmaetrvv7T0YdY/iwfqvssbDEgQBwy3lWHxrDa5/0YGPXj0Gy0B6xYMEAgG++93vwm6347/+67+mbYhgNhNwkq09xoGokqYX+ZMK31hTAa1chMa3LtAHXAPAGMc4AASDFI4fHsYLz5rQ1luItWXH8LXvr0dx2zH0FKggIPiozKAW/GjNfEiGBsNli3kEgTr1QhwyfwF30EfvIheW4qHH78aitQYc3X8Gv//pPnRfS1zmOBgMYt++fVi+fDn0+nF2DhPgrrvugkgkwquvvpqza4xH3mVO5opME3A+uzoIpVQYU5VN5KVjrv2i0OIhSchuXIFnTmxNDrmGdlT299KJMzdcvShhiYN1LlgM5bUrEDjYBWG1tATWwAh6vWb0tQ+h6b8+RPHcAqz75mrw+Zn/OW6YruP4qWO467a7oVKOr7GXzy/GmvuXwTPqwwe/P4qWQ5cQTEP71uv1aGhowPHjx3HgwIGM551LspmAk2ztMRq3MgeCWyTg49u3zkXT551oHxwB3INxjsn+Hg/+8OsOfHJwGAuWKPBwvQULNKfB67oKOO04pQ6iWloMUTK7OAuu8gqQfAFUVy6Ej92uXgwP6UOzJRL2yxfwseqexXjwe3eBJCm88h9v4cNXPoVvTGu2w4cPY2BgIO3yC+kikUhw//33o7m5GV9++eX4b8gBeZc5Od04fLEft1RoYkIBxYzgDmnckr4uCFwjcJdXxb6ZL8EoJYHHehOuoAddnmGUiRNHejgXLQFBUdCdPsE6l0VyOl398PVW/OWZ96HWK7B+gkI7EAjgtaZXUVo8B8uXpp4yrC5Q4vavrcC8lZW4etqE5j99Bk+yGhljWLZsGW677Ta8+OKLuHYtM5v5TGDI4QFB5EbjBoCvrJwDlVSIX753kU5lDwnuYIDCJweH8Idfd8DvI/HAN4ux+jYNhKqQMtJ6CKRQiNMqPxbIKtgvkAyBEO6yciiifDdaoQKL5ZU4MHQ67nRdsRpf+dt1qKtfirMfX8HL/9KEnjY69tzr9eJPf/oTli1bhjlzMnO8p8PatWtRU1ODX/ziF7h5c2pr4WebGS+47S4fzrSZUWuI1ZLF3l5QBD8cxy2/ch4knw9vaXzn61FCB5m3CxednaBAoZply+nX6OCqmIuCo4dY51MgUmOdfRGu//oahGIB7tm0FoIEVeXS4QPj++gf7Md9d9en7aEneASqlpahbuNSOC0j+OCPn8Dca0v5/Q899BBKS0vx05/+FFevstTXmOEMOeiqk/wM0rVTQSLk46/WVuLN45fogmdiLQb6PPjjbztw/IgZy1ap8MBfFUNfEDKzyebQIYOtx9BdoodUrEClJPO2Y66KKiivXwXPE4kpX62aj6uuLvR6LXHn83gEFq014Kv/13rw+Dy88u9v4aNXjuH3L/0BZrMZDzzwQMZzSQcej4fvfve7kEql+PGPfzyj1ueMF9yfXBpAkKRQWx1rTxO7u+ATFobDqpQXvoCnbC6oBLWq/eJizBEOonngPNQCeVwoYDTW2rUoOPYxhLb4BT3SNYorf7iKB95eAqtoFEUPl0Iin1jrpQuXz+P95ndxW91tKCrIvJ+lukCJtQ/cAr6Ajw9f+RQtxotx29xECIVCfO9734NWq8VPfvITHDhwYNravHNFn9U9fteaCfL12nIYFCOwegrwwdFi/P5XHfC6STzwV8VYUacGPzpShC8CRPOBngEcLeRhhcKAFANJEjJqmA9ewA/l5fPhY0sVcyEkBPjEyl7YSa1X4Ct/eydqbp2Dve/sxjvvvY1lhpWQizMra5wJcrkcW7ZsgV6vx7/8y7+gqyuz6pzTjRkvuJvP9aJMJ0OJJjYUUOo2wSemBR1/1AnF1fNwGRI3a+UrSlAmsuCo/TTqVAviHZNRWNbeAYpHoOyNv4SP+Rw+XP3va/jy/zsHZ7sTc+4twcW7+vGy5QOQyEzIkSSJoyc+wQt//A2q5xpw6+rbMxonGolcjLVfuQU1yytwrbUDb/3uEC6dbBs3SkAmk+Hxxx/H6tWr8fzzz2P79u1wOp1J3zOT6Da7UKjKXe9DiqLg7jbjIUqA37X+DOcvC7FyjRoPfitKyx4D2S8AxQM65+jCDX8zJaBWw1NYDM2XEZu2mCfEEsVcNJtbWVdwIBBA09t78MqBFzDk78Yt1avAs8jx5m+N+OSN0+i5MZBRVFO6yGQyfO9734NcLse//uu/zohIkxktuIccHuw/eRN3LIyvPywfuQSPpAKgKBQefBsAMLJwWcJxAtJyCAkSS0dd2KBdmfSaQbkC5jvvQXnTa+CbLRg4MYgvfnEW1ss2zNlQivn/Yx50y3T4etFtOHvoMxyxpO846ezqQOPzv8Cfm17F4gVL8NBXHs6oqloiCB6B6mXluPPhVSgq1+GLw5fw9u8O4+qZdvi9Adb3CYVCfPOb38R3v/tdnD59Gt///vexf//+tFufZVKydzJIFhnQNuBEsZqlPvYEoCgKlrYBnH3lE5zffRziQBBr5zSjT2/HvGXyWC07iiAZhOOyCYQG2KBfmJFTkuHwiVMAgJGa+VBdPAde1N/zVvUitLv7cNJ+Je59Pp8Xv3n5V/jksyO489a7sPnvnsDGB+tx17fqsHB1FWwDDhzZcxL7f92M1sOX4LTkNqxUKpXi7//+7+H1evGjH/0I7777Lvx+9h3ldF2HDDO6v/0v9p8HBeCba2IjRYS+IUg9NzFU8Feo+s0voLrQCsvt97J2r+kSSCHgC3G/1YmPzw1j4/LSmG7e0bhGKJyu+RooUwmu/+eXcAUlUC9Uo/SuYgikka/bICsFcXoIf7rXiHt1K8DD+HtZ54gT7370No5+9jEKdAXY9I2/RkVZhk6ncRBLRVh8qwFzF5ei7Xw3zjRfQIvxIoor9SitLkTZvGKoC+PbcC1fvhzV1dX46KOP8PLLL+P111/Hhg0bcNttt2Hx4sWQSpMLuN27d2dUtjfXsAlujy8I04ATG5dnp+iSb8QDZ68V9i4zhq/2wmMdhaxAieoNS7HUtRN8jxXmdj9+99E1PLhqDhaWqmN2gFQggO5Xf4W5Qy5gCTA3SCJBz/eUOXLiNDbcvhbORUuhO30C5U2vof/+h+ArLsEiWTkWyirwM9Of8YOyB/C1glsh4gkwZB7C7/+0C929Xfjm176NyvKIxi8Q8lGxsBTlC0owYnWhp20Q11s7cOnEDehLNSitLoRKr4BUKYFMIYFCJ5+Q4z4anU6H//W//hcOHDiA3/3ud3jjjTfwne98B+vXr4+rVT9d1yHDjBXcZ9qG8erRNjx2Z3VMDW7je2/gh/OOgwIP8uaLUF48i76HH4Nr3qK4MSiKwinveZz1XMZDYh2+EuxEe1czfnxYhX/b8lXopGKMOgCHmYJ1iIK5j8KoAwBEUBSsRMXwZShXlYCsLYBPEl+AqFCsRqdnAP/R9hru0a7AOu1SCKKqBno8bgxbhmGxWnDp6kUcP/UpCILA+tvvxspbasHn83Hk6GHcu37DhL8vtnFkKiluuXM+5q+sxMBNMyz9dnx59CpaD1+CXC1DSXUBdMVqKLRyfN5yHN/4q4ehVCrR0NCA++67DydOnMAnn3yCd955BzweDzU1NVi4cCGqq6tRUVGBwsJCaLXarHfqnqwb78S1IQRICovKkodgGt9/Exse+AY81lG4raPwjbgR8PgR9AUR8PrhtY9ipN8Or4OueS6UiaAo1WDO6mq0XD6OW9QyqAdOY7jwIXx9dQWOXRnEH4+0YW6hHHcvKUFZoRA3ey5A3mzE3PZhXL1tCQziTqgs5/H+hT7cu27dhD5nUKHE4L33o/DoYWhbT2FgwwOwr6zDdwrW4+2OT/HyqX14w/8uFjl06L1wA3KZHA3feBRXrl6OEdwMBEFAqZNjka4a81fNxVC3BYNdFtz4sjMc2XT55jksrV4JbZEKBeVaaIvUUGhkkMrFEElFEIkF4LM0r/jggw/w4IMPxh1XKBR45JFHsH79enz00Uf45S9/ieeffx5VVVUoKyuDwWDA0qVLQZKZm3AmY+0RVB55kh5++GG88847cccPX+jDN545kvA9dywsxA8u/BGbLtBbn693Evibu7fgujq9ThvRvPbx7/A/7vlRwtd4VABKDEBA+aAaYW88EOTz8O9nD2Dh0kUgA6mH4BUoiyEURH6IDh09iPvW35/65FlIZxyKouBz++MW96enj+CuNfcmfI8n6IKHHE34WjStra2ora1NaR5jqaurw09/+lNIJJKYtcK2btKhpKQE/f2RpJJRbwDf/dVRHLpAH7t3aQmKzTZU9SbugpRszURDEoBbLEJAIAjbjt/58Fn829fuAQA4KSVIgoSfT68ZoYeEuj/WB0ERQEDMgwAkCADPHXsP/7TuIYwIxCDTtI7+7pM3cEfdWnQkiB5JhlCmBMHj4dRnx7D2jtgfDZVABiHBrjNSFAUySOLwsWasq7snresyJFuL0fhIL1xBR9zxdNfhypUr8a//+q+Qy+Vx6y0b628seSW4y8vLEQhE7KxlZWUoKytDr9WFsx3xdRI0MiEIgsB863UUuOmU9W4/gaBmHnz8zDU8u8sGtUwTfs6nvCBCN8lYCIoCj+Ur7vWPguKN/8tO8HgQ8AQgElzB5XZBJk2/7Gu2x6EoCi6XC1JJqmNQIKkgSMR+fo/Hk7H2LZVKsXTpUgwMDKCtrS08jkAgQHd3fH2YdNDr9RBGRRyVlM5Bx4gQI54AxAIeZGIBZD4fZL7EdtPoNRPk8UCmEbbpGh1GsUIECgRI8GiBTtDfGxEE+P6o9UUgxllIALC6R6CVKkASvLRd4Ta3E6SQBx/F7t8AAD4IEAQPFI8HKuoqHrcbkjHmMT7BS7iWx5JoTVIUBSoFZdjtSWctRs1NwEOA9GN0dDStdSiRSLB06VIIhUKcP38+pltUNtbfWPJKcHNwcHBwzPCoEg4ODo6ZCCe4OTg4OPKMGRtVMtXs27cPFosFbW1tsNls2LlzJwDAaDTCZDJBp9PBZDJh69atSY/PZLZs2TLtv5epvn4u4dZofqzBhFAcWaetrY1qbGwMP29oaAg/r6+vDx9vbGykmpqakh6fqTQ2NsZ85un6vUz19XMFt0bzZw0mgjOV5ACbzRbTfWPNmjVobm6G0WiERqMJH6+trcWePXtYj89UTCZTzOedrt/LVF8/l8z2NZova5ANTnDngNraWrS0ROo6tLW1oba2Fq2trTGtmnQ6HVpbW1mPz1SMRiPq6+vDz6fr9zLV188ls32N5ssaZIMT3DnGZrPBaDRi27ZtMJsTNxJmOz4TMRqNeOSRR2KOTdfvZaqvP1nMtjWaT2uQDU5w55jHH38czc3N0Gg0rK2actnCabphs9litp4A++ef6u9lqq8/Wcy2NZpPa5ANTnDnkGeeeQaNjY0wGAwwmUyora2FxRJJHbZYLKitrWU9PtN45plnYLFYsG/fPuzbtw8mkwn79u2btt/LVF9/MphtazTf1iArU+0dnak0NTVRLS0t4ec7d+6kKIqiamtrw8eiPdRsx2cqLS0tMd766fq9TPX1c8lsX6P5sgYTwaW85wCTyYSampqYY42Njdi6dSuMRiNaW1vDGk50rGii4zMRk8mExsZGGI1GNDY2oqGhYdp+L1N9/Vwx29doPq3BRHCCm4ODgyPP4GzcHBwcHHkGJ7g5ODg48gxOcHNwcHDkGZzg5uDg4MgzOMHNwcHBkWdwgpuDg4Mjz+AENwcHB0eewQluDg4OjjyDE9wcHBwceQYnuDk4ODjyDE5wc3BwcOQZnODm4ODgyDM4wc3BwcGRZ3CCm4ODgyPP4AQ3BwcHR54hmOoJpMPSpUvjir8DQE9PD8rKyiY8frbGyeZY022cbI41WXNqa2vDxYsXJzQ+29rLhGx+7lyPO5vHzNa42Vh/cUxtA570+PrXv57W8WyNP5VjTbdxsjnWZM0pG9eZjt/fZIw7m8fM1ri5mFteadxsPPbYY1M9hTiyNafpNk42yeaccv35enp68PDDD8dcb7p9p7mYT76MmSsymevu3buxe/fu8POenp5sTokm6z8FOSRXv6qTNf5MId++p+mmNeXb95cPTOfvNBdz45yTUeSTJjCVcN/TxOC+v+wz277TvBLczHY1ehsyLkE/cP43gN817qmz7Y+fKfnyPe3evRsPP/xwbraqmTJ8Do+t8k31LGYc+bIms0Ve2bjLysrwzjvvpPemroPAp/83IFQAi/4+J/PimJ4wduho2/SU81EDYL8O1GwChLKpng1HnpJXGndGuAfpR59jaufBwQHQQhsARqfRLoAj75j5gtt8jn4kA1M7Dw6Oke7I/139UzcPjrxn5gtuRxv9GBid2nlwcAx8Hvm/Z3jq5sGR98x8we210o8pOCc5OHKKowPgSwAQgMc81bPhyGPyyjmZEf6Qph10T+08ODg8Q4CkAPA7AffQVM+GI4/JK8HNhAOmlbUWCGnagdktuHft2oWdO3dCp9MBAHbu3AmDwQCCIFBfXx8+z2AwYOfOneOO98wzz6C5uRkAcObMGRgMhvDYzHG2sQmCwObNm2GxWGCz2bBx40Zs3bo1a5+Vgclgy0Y4YEZrbyweCyBSAjxhXtu483Etbdy4EU899VR4jH379mH79u1oaWkJn7N69eqY5xMlm+svjqyn9OSQjDKQXplLUb8FRTX/Tdbnky80NjZSDQ0N4edWqzX8/2RLoLm5mQIQc34iamtrqZaWlrjjbGOPPd7Q0EA1Nzezjr9z506qra0t6RySMW0yJz/cRFF7V1NU01qKem/6ZvolI1/XUmNjI7V169bw882bN1O1tbXh+VitVqq+vj78/8bGxqRrMh24WiWZQHrpx1lsKtm+fTva29vDzzUaTUrva2xsRG1tbY5mFcFmsyV9ffPmzTmfw6TgcwICGSCfA/Qdn+rZZES+rqWGhgZs2rQp/NxisaC+vh5GoxENDQ0wGo3YuHEjAOCpp57CmTNnJmW+mTLznZMBT+zjLKO1tRUGgyHpDbZp06bwv9bWVgD04n3qqafCW9ZMSTQ2c3zjxo2oqanBxo0bY7bB0ezatQtbtmyZ0BymDX4HIJACmkXAaDfgzq/IknxeSwaDASaTCTabDa2trVizZg02btwYNsU0NzejoaEBAG36qaurm9Bcc80s0Lj99CMXDshKU1NTzPPW1laYTCbU19ejsbExq2NHH7fZbFi9ejWrfdtoNOKRRx5JyU6aF/hHAEkhIA/Vd3a0A9KCqZ1TlpmuawkAHnnkERiNxvB8amtrw0oBY1vPF2a+xh0MmUpmqXOytrY2rGmkyvbt22EymbBp0yacOXMGmzZtSuv9qaLRaNDQ0IBnnnkm4evMVvbRRx/N+rWnBP8orXFL9PRzV+/UzidN8nktAQhr2KdPnw6bQQwGA4xG47TXsMeSE8E99kbr6OjIyriMZ5/5N26xKYoEqADtxZ+lghsAtm3bhscffzzmWLKbp6mpCS0tLWhqaoLBYEBTU1PKtsxM5rZ9+3bW13fu3InNmzfDaDSmPCZTXIr5N22KTPlHAIEEEKno53lmKgHyey01NDTgzJkzMSabjRs3orGxMWzfzhdyIrij7U8AYpwCE4EpMsX8GzcsizGTCBWzWnBv3boVGzdujPlnMpmmeloAaE1p8+bNeOqppxK+bjAYsHfv3rQ0osceeyxmneSipVVGBFwAXwrwBIBAEUkOyyPyeS0xrF69Ovx/xjEZbRfftGkTjEYjGhsbsWvXrpzNdyJMio2boqjJuEw8wVD5TIEiEs89S9m8eXPC6Izx/japxLWyncM29tjjyWyfM8a+TVEhwS2hn4uUeSm4gfxdS4nGNxgMcWOw2dKnE5Ni4yYIYjIuEw9j3xbKZ73g5pgYGdWCj4b0AVQQ4Ivp50I54LNnb4Ic045c1oPPicZtNpvx5ptvhn/JbDZbzHMA+Na3vpWLS8dChjRuoQJwduT+ehwzloxqwUfDlF4QSEOPnOCe6eSyHnxOBHd1dTWefPLJmGPRzwmCmBzBzZhKREog6KGdlcTMD6ThmIb4R+hHQchUIpBxgpsjY3IiuLOZ7z8hwqYSBf0YcNNbVA6OyYbRuPkhjVso45p7cGRMTtVPh8OB/fv349lnn8Wbb74Jp9M5ofHStjOGTSVK+pGzc88qplXPSX9o7XOmEo4skDPB/dJLL6Gqqgq/+MUvcOrUKfz85z+HRqPBc889l/GYjJ0x5epsYY07JLj9XPbkbIIJC5wW4YBhwR3qMymQ0rVLODgyICemkjfeeAM7d+5ES0sLqqurw8dNJhMeeeQRGAwGfPOb38zFpWOJtnEDXNo7x9ThGyu4ZRFhzsGRJjnRuHfs2IHDhw/HCG0gkl76i1/8IheXjYfkNG6OaYLPRj8yPhaBLOKwnGL6rC7cHObujXwiJ4KboiioVKqEr2k0mslLyBlrKuFs3BxThddGJ9/whPRzgZRej2RwSqcFAF9vPIwVT74Dt49rqJ0v5ERw6/X6Cb2eNRjBLeI0bo4pxmuN1CgBIk7KbJrvKAo49BZgS72fpW3Uh6u9DgSCFI5czN+uPLONnNi429ra4hJuopm02gZjNe5Z3EyBY2JMuHWZ1xoJSwUiYYH+kViBPhG+OA784zeBb30f+N8vp/SWG/2RkMTTN8z46qry7MyFI6ety3IiuNVqNX7+858nfT0T0r55mOYJYRs3ZyqZTWTzxplw5qTHElmHQKRmSTbNd6eO0I+Xv0j5Le2DtJ19QakK1/q4uPJskneZk7lKwEn75iG9AIjQtpTgNO5ZRi5vnLTx2WKTv5iaJdk035ku049dN1J+S9uAE2qZEDXFCrQPclEu+ULO4rg7Ojrw8ssvo7OzM1eXGJ+Ah75BCIJONeackxxThc9OJ90whG3cWVyT3SZAKAJGHIDDNu7pQZLEofN9qC5SokAlQY+Fuz/yhZwI7kOHDqG2thZ79+5FbW0tjhw5kovLjE/QE/Hi80Sztu8kxzTA54jEcAMAX0Q/ZrNOfP9NYN4y+v99N8c9/f3WHnx+fRgP15VDrxDDMuKDL5BalIsr6MFZR9tEZssxAXIiuJ9++mm0tLTg4MGDOH36dNI+cDkl4I5sSflizlTCMXX4nLGmEl5oXWZLcAcCgHkAmLeUfj7QNe5bTrcNo1gtwW3zC6FV0D8kg/bUlJvvX/glVn3+Q9xwTYNyArOQnMVxM8k3iQqVTxpBD61pA5zGzTG1BEYjDkkgolBkS5kY7gNIEqheDPAFKWncnUOjKNXSJhutPHXBHaSCeGvoMwCA0Zy6I5Qje0xKHHe2GimkXWQq4I5sSfkiWpBzzBqmV5Gp0UhJVyAiuLOlcfd00I+FJYCuEOgd37fUZ3VBK6fnoWEEt2P8e+SM/Ro8oQJuXzqnR9uy2cakxHHbbDa8/HJsXOkPfvCDtMdNO6ok6IncIDzRrO47ORuZNlEl4bZl0sgxxveSLWXixgVa0y6cAxSUAD3t475lyOHFirm03V0tFYWOjT+fd4Y+h1ogxwqlAZdGpzD4YBYzKXHcarUaL774Yvg5QRAZCe60iTaV8EWRhBwOjskk6AVARZQIgI50ypb5jqKAPS8Ai1YCIjEtvDuvjfs2y4gXSin9AyIU8CAXC2B2jn+PHLF8iVXKeZgrKcL7w6cmOnuODMirOO60iTaV8IScqYQjYxgzHUNaGZRMWnu0jRvInvmuqw24fh74h1DxttJKOhmHJAFeYmsoSVKwu3xQSCIiQCMXYXgcwU1SJL50tuHv5mxEkUiDYb8dQz4bCkWaiX+OGQKT+MWQN5mT04aYcEAhp3FzZMyEMieZjN1ojRugI0uyYb67do5+rF5IP86pAjwu2s5dXp3wLSPeAEgKYY0bANQy4biCu8M9ABfpRZW0BCUiLQDgvLMdG/SrJvwxZgpjf9RzYaqb2Q0YA9FRJZzGzTFFMJEjYwV3tsx3XW2ARAaoaEGK0kr6sf0K61vso7RzUS6O6G4qqRDD49i4r47SYYaVkiKUSfQQEgJc5Ozck87M17ijo0q4cECOqSDAIrh5WRLc/V2Avpi2mwOArggQCOhMShYcbj+AMYJbJhrXOXnN1QMxIUSRSA0ewUOFpABXR7sn/hk40iKvNO70wwGnh6nEHyCn5LqznWkTDhhgMZVky8Y91AuodZHnPB6gKwb62DVhRnDLogS3WibEkCP5PdLm7sUciR48ghYdpWI92t1cOdjJJq8Ed9o9J0lvbMo7OfmC+0qPHXP/5xvY+1nHpF97tjNtek6yatzC7Ni4h/tjBTdAPx/qY32LMyy4+eFjGpkI5pHkPyTtrn4URzkiC4Vq3PQMpj9njgmRE1PJ008/nTTpRq/X4yc/+UkuLh1L0DvGxj35gvvTywNwuv3Yf+omHrmjatKvzzENyLWpZHgAWFIbe0yloVPgWQgLblFUVIlMBJc3iFFvIMaEEk2nZwDV0tLwc71IhaO285nPnSMjcqJx19TUwGAwwGAwoKWlBVarFTqdDjqdDi0tLWhrm6TiNMFojVsQaR6cTQIBwGlnfbljiA4F6xyaHv0FOaYAVuekMDsp77YhWlBHo9QAliHWtzCmEokoonGr5fS9ksxB2e0ZRqEoUk9fL1TB4nfCR/rTnzdHxuRE43788cfD/9+3b19M8s2TTz6JRx99NBeXjSfopW8OIHvazVj+9+PAx+8BBzsBqSzu5QE7fWN2s5TM7POa8cuON/CzeX8PCeNI5ZhZMOGAvDF/32wk4Pj9tOKgGNOcRKEG7GdZ3zbi8UMq4oMXtTPWyCLZk3MLFXHv8QR9sAZGUCCMdOzRCujzhnx2lEkKJvBBONIh5zbu06dPw+mMFGi32+0wGo25vixN0Jd75+Rb/w3YhoHzJxO+zAhu26gPHl98ycyXuj/As537cNA8TZKWOLJP0A0QvMhaZOCLJl6P2x7qLxknuFWA3cr6thFPIMYxCQDqkOBmi+Xu89LXKhBGrqUNtWMbZLrYc0wKORfcTz/9NKqqqrBt2zZs27YNdXV1MRp5TiHHmErILJtKRqM6hrDEzA7aPCjX0Zp4ny1+Wzzgo2+uFsf17M6NI6ukHdEUDVNeeKzfh5+FBBymMbBiTN9KuRJwjwD+xGve4fbH2bFVspCphE1w+ywAaLs2gzqscdvSnfmMJ5dRTTmP4966dSvq6+tx6NAhAMDevXuxalVmWVZppx0H/WM07iwLbqYiG8BajW3A4cGKuVp0W1zoHBpBdVHsFrTfSwvuk3b2ZAmO9MhFyvGEMicD7kj97Wj4konbuK3D9ONYjVsW6m9ptwIFxXFvc7ppU0k0Qn7yeiV9Xlpw66J6Z6pDXX2G/Ox+ntlK3vWcHEttbS1sNhs2bNgAh8MBh8MBlSr9ztYZ9ZzkhT4iT5j9cMD+UM3jghK6HvIY/AESlhEvlpSr8fn1IZy4NoR7lpbEnMNo3F84U+8TyJGcyUg5TougO94xCdDHJtrA2kELU8gTaNzM6wkEt8Ptj4koYVBJhUkFt5AQQMWP+HKkfBEkPCGGfVyj4ckk56aS/fv34/7778eWLVsAAGazGZs2bcr1ZQEyCFAkQERp3BRJH88WvZ10Kc1yA2COj2UdsLtBUUChSoLVBj0OXYgX7n1eMwqEagz6bHAwxYg4ZhZ+V3yBKSCkcU/Uxs0IbmXscUbjZuk9aXf5YmK4GZRSIayj7IJbL1TGhfpqBArOVDLJ5FxwP/XUUzh48GC4I051dTUsFkuuLxtxRPKjBHf08WzQ006nGqu0gDU+9IqxaesUYswtUOBGf2wXbZIi0eM1Y7mS/m5uurlEhhlJMKpKZTR8ycS7vNstgFQB8McIYVnIJOe0JX6bK97GDdCC2zKS2KTY57VAK1TGHVcL5Bj2cxr3ZJJzwU1RFBwOR/hX2m63T04rM0ZA83IsuAtK6JvEEe/BH7DRoV46hQgFKjHMTm9MM9YvnSZ4ST/WqBbQw4W89hwzjACLqUQgzY7gVsQL0/EEt8Plg1wijDuukAhgGUl8j/R7LTH2bQaVQIYhH2fjnkxyLrh37NiB2tpamEymcFTJP//zP+f6shF7NiOwGY0nm3buLhMtuOXKhIJ7yOEBQQAqqQjacIxs5PqfWM9BTAhxt3Y5gIjzh2OGEXCx2LglABWYmDLhsEbMItGIJQCPD4wkFqhsGrdKKoQlSVSJnk3j5gT3pJJz52RDQwMMBgP27t0LnU4XYzbJKWGNOyrlPfp4Nui/CSxeRWs3ow66E0mU/W/Y6YVKKgSfR0RiZB0elIXCA2+4elEmKYCUL4aKL0O/jxPcM5KAOz75BgAEISefzwlIEwj2VHDaAJk8/jhBAFI5a1bv2CYKDEqJEDZXYlPJgNeKFQpD3HG1QI7Lo+M3J+bIHjnXuM+ePYva2lrs2LEDBoMBu3btQkdHR0ZjpRVLy2SkhTvgiGKPTxSfl46h1RTQN0gwCLhjHU2WEVpwA5GC9dbRyE3R6zWHNRidUBkODeTIDtOnOuBoYo1byAjuCdiHnfaIWWQsMnlCU4nbF4A3QCYU3AqpIKGNm6IoDPnt4YSbaDhTyeSTc8HNRJAcOnQImzdvhk6nC0eYpEta1QGZcpljTSXZ0riHQ6UstXpacANx21LrqC8ssBkBHm0/HPBaoQndCDqhkjOVZJlsVgecUAKOn8VUwmjcfmf8a6nisEbW31ikCmAk/kfBGhLMigQ2bqVECI8/CLcvEHPcFhhBgApCI4gX3GqBHNbAyOT4rvKIXCoOk1bWdefOndi2bRuefPJJmEzsBd6zxtiKbIzGna0uOIzgVunoGwSIF9wj3rBWIxXzwSPo1PfwEH57+EbQCpXhlGKO6UfaJYWjCbCEA2ZDcI/YI+tvLBJZQhs3s+tTSeMFd0TBiNW6mThttTD+R0IlkCNABeGYaPr+DCOXZYVzLrjVajW2bdsGo9GIhoYGAJikqJKQ4GYy1rKtcVtCoXtqXcTG6IzXuOVi+kbgEQTkEiFsrkgVtWG/A5pQ5plOqES/jzOVzEj8o4AgieD2TUBwO23sphKpLLYsQwgmTjuhjTvBzhCglQwgkikZjSr0OSxcSOCkkXPB3dTUBJ1Oh6amJlRVVaG9vR2NjY25vmxE42ZumGxr3EzctlzFGnplHfFCKY3cHAqxAPaQ4ydIBWH1j0AVJbgHOME9MwmymUqk9ONEEq9GHElMJbKENm5Gm1Ym0LjVMhbB7WMEd3wFTOaYeSI7B460yElUSXRKe3V1NZ588snwa9XV1ZMTVTK2XRSjcWej4whA14iQKejefuH04ljBa3f5Y+yIMrEAjpDGbfE7QYEKa9wagQKOgAte0gdxoggEjvzFP8qeOQkA/gxrtQcCdDd3VsEtB6zxNXQYoZzQxs1iKmGEspIfL7iZY5zGPXnkROPWaDQ4e/YsfQEeD3w+P/yPeZ4JaTmImBoQzM2R7XBA2zBdrB4ARBJAKIoU/Alhd/ugiIqVlYv5cLjpG4KxGTLbTGYLytV8yB7TJqqETXDzBPS6zFRwj4bWSoI68AAAiTyhc9I2SocC8nnxXarkYvr4WI3b4ndAzpdAyEsQ+81p3JNOTjRuq9UKtZquVkaS2WuUm1aRqYALIPhRRaaybCqxmSOaNkHQae/mSNNUf4CEyxuEPMqOKBUJ4PTQ3vqxNkN12E7o5ArSZ4lcVmdLmaCPTrJhzCJj4Uszr8nN+FTYnJMyeULnpM3lgzKBtg0ABEGE0t7HCm5nTHGpaCQ8EQQEH1ZOcE8aOdG4GaE9pYzNVsu6xm2OrcimLQQGusNP7SHNOjo7TRZl4zaHtpVMVElEa+E07hkFo02zCW6BJPMKgYxQTpSAA9CmkpFQYlgU9lF/jEIxlkQVAs1+J5QJ7NsALezVAhksAU5wTxY5EdzR5hEejxf+N1FTSVoEPbGCmyBC7cuyqHFHF68vKAa6Ir00GVt2dJcRmYgfbtJqDplEFCHnKWMn5LSWGQbz92QReuBLMndOMo5H1gQcBRAMAJ5Yv47d5WNtBgwASkl8Eo7F74SSz/LjA3r9mnNk5nM5Pfi46RT62tl7aAL0Tv+5554L1/6fyeREcJMkiWAwiGAwCJIkw/+Y58FgFkursjFWcAOhVlFZEtx2S+wNU1QO3IzU1GaasY7VuJnjZr8TCr4UfIL+EWMEuI0r7TqzYEL9WE0lE+iCEzaVsGjczPoc4zS3uXxxbcuiUUjiTSVmv51V4wboHWOubNzNrx3H8bdb8eft78JhSewPoCgKzz77LJqbm/Hcc89NvV8jx0xaAg7TQIH5l3MCngTNWbPYd9Jhie06UlgCmAfCae9OFsHNHLf4nTExsXyCDxlPzGnc05SMMyeZVPAE8c8AQoI7Q1MJI5ATFZkC2AX36DgadwIb97DPETbnJXyPQJYTM5/P48flk21Yesc88HgEml/7LOF5X3zxBVpaWvCd73wHMpks825FWSSvMycPHToEnU4HjUYDrVYLjUaD1atX5/qytMY9tjlrNrvgOKyxGrcu1GUkZOdmNGtpVLF6mZgfJbgdUI7RwpQCGadxT1MyzpxkzAcJanwACHXByfBv7rDSVQAFLEKY8cHYY0spON3+hE0UGFRSYUxNHYD2vSSK4WZQC+Q5aRjccakHwQCJmuUVWLVhCa6cMsF0vivuvN27d6OiogLLly9HbW0tjhw5gkAgkGDEySOvMyefeOIJHDp0CCRJ4tvf/jasVuuEe06mpPUEvQkEtyg7GrfHTReZirZxa0ORIEO9AFg0bpEA/iAFjy+Y0Gao4EtgzTQ0jCOOaREOmIrGnWnfSbslvtdkNMxr9thSCokaBUfDGlXC9hkAaATynHTB6bzUA7laCqVWjuqlZSiu1KP5T8dBkRGHa3t7O86fP4/169eDIAisWLECDocDFy5cyPp8pguT0kiBEdQEQUCtVqO9vT2jsdIuMpVI485Gw+DwFjVKi1Lr6MdQDROH2w8Bn4CQH/mKmRRju8uX0Euv4Eth5zTurJFLjSdlvDaA4CW3cWcaVWK3xPeajEauoJ3ytljBTWvcyQS3ALZRH8iQcHQHvXCR3qQat1aoxKDPlvVyFu0Xe1BcqQdBELRQvnsRhnttaL8YieD64IMPoFQqsXTpUgC0nFCr1Th16lRW5zKdyLngrq2txeHDhwEAdXV1uP/++2Gz2XJ9WRaNO0s2bkZwR980EhmdhBOqYcIUqo/uz8fcLDaXj3b2jNG45XwJZ+OeafjsgEARU6c9Br4486gS23B8r8loeHx6VxjVVo+iKIx4AwkbBTMoJUKQFKIc6UyyGLvGrRMq4SZ9cE60h2YU9mEnhrosKDUUhY8Vlmuh1Mpx9Qyt/AUCAXz88cdYuXJlOFqNIAjMnz8fZ86cydpcphs5b6Swd+/e8P+ffPJJGAwG1NfX5/qykyS4o24agqAzKUPZkw53vAOI0bgdbj/MfidWKufFvi6QwsIJ7pmFzw4kqKgXZiJ9J61DgHKcnAmlBrBEBPeoNwCKAqQidht3dO14jVwUFtyJCkwxMC3N+ryWpAI+HY6/3QqhWIDy+ZEu9QRBoLS6EO0XaI375MmTcDgccX4zRnAPDw+joGDmJbRNelTJxo0b47pE54REgpufJcHNOHvGxs/KVeFtqd0Vn+TA1IawjfpCNsPYraeSL4WVS2KYWfjs7I5JIBTHnaGWahlKbioBAIUmUskSEd9LUlOJJLbQ1NjyDIkoENLz6M1SaeKrLe344shlrLp3MYRjdgdFFTpYBxwYtbvx5ptvYu7cuSgtLY05Z/78+SAIAidPnszKfKYbORfcTz/9NHg8Hqqrq1FVVYWqqqpJal3miXR4Z8hWVEkijRugm7aGHEGOBEkOzPOh0RF4SF+cBqMUyDjn5EzD52BPvgFCmZMZ/s2tQ5F6OWwoVbGCO1RyISWNOyS4w6YSlpR3ANBnUXBTFIXDfzmJOTVFmL9qbtzrBWVaAMDRQ8dx4cIF3H333XHnyOVyzJs3D0eOHJnwfKYjORfcL730EqxWK8xmMywWCywWC8zmSWgYwGYqyUYCjt0SCsMaM75MGdbG7S5/nB1RKuKDzyPQ7aLPGevsoTVuTnDPKHwOdsckEKpVkoGphKLo3V2appIRRuNOYuMe20zB7HdAQPAgT1QoK4SUL4aMJ85K+73etkFY+mxYcmtNwt25XC2FWCbCQeOHKCoqwpIlSxKOU1dXhwsXLmQcDDGdybngvu+++7JWu2Ti4YDC7KS828yJw7Ciur3bXb44UwlBEJCLBehz2wAA6jFtoJQCGbykH65speXPcqZNOGAyG7dARtu4qTSLsTntdDq7Spv8PLmKThZj3uYZ31QiFvIg5BNRphI71AL5uCZOvUiVlfZ7bee6IJIIUVSpT/g6QRBQFypwreMKVq1axTqvW265BWq1Gm+//faE5zTdyLlzcsuWLdDr9aivr4dOpwsff+GFF9IeK63qgEFvgszJLNUqsY+pU8IgUwDtV+hTXP5wN/doFBIBBrw2QIRwLW4GJsrE6h+BLIl2w5Ea2awOyCgNY8ceF58DUFSwvy6QAaBo4S1KEiEyFluohHCyOG4g5HeJCFNH2MbNbiohCAIqmSgsuM3jxHAzaAWKrDQD6b7Wj8JyHXgJys4y+EQjCJIBLFq0iPUcPp+PNWvW4OOPP8YTTzwBiWRy7qndu3fHKJe5UBxyLrifeOIJ3Hfffairq8v1pWIJJIjj5osiCRETYWxlQIYojTtRVAkAyCUCDAUGANCxr9FEVwjkSrtOL9JSGqIZz1TCaOM+e3qCm6n9Pp6pRK4E3CN00wWBIGHxs0SopMKwqWTYb08aUcKgESombCqhKAp9pkEsqEvuB7P5hsEnBFCPs+NYtWoVjEYjTp8+jbvuumtCc0uVsT/quSgrnHPBrVarY0ICJw2SLXMyG6aSYRaNW0lXbCPJcBz3WORiIcwBG9QCOURjitIzCTmcg3IG4XPQcdxsMD/eXiugKE99XCYbcryokui2etoC2F2+kCkkuZVUKRXCHNK4B322pJUBGbQCBdrd/eOelwyHeQQelw+6kuQ/SEP2Psj4Slh77ZAvYJ9bQUEB5syZg2PHjk2a4J4Mcm7jfvTRR/Hmm29OTmGpaILeSLsyBp4wO63LrMPsGjdJwm+3w+0LJqx5rJAIYIUtHD4VTdhUwoUEzhz8zuQaN6Nle9K0DTMhqYkUiGjG9EOl090TN1GIRiUVYthJKzlDPju0yUIaQ2iFSgxMMO29v5PeSWiL2D8XRVHo6r0JpUSDgZvjBzosW7YMJ0+ehNebpTpF04CcC+7t27fj29/+NjQazeTX4x5r4+aLsxPHHd22LJrQttXR1wcAkCfoMqKQCODg2cLhUzFvF0Rs3BwzgKCPXofJnJOi0DrwJK81HYfdEmmZlwymyUKoBKxt1AdlkiYKDCqpEMOOiMY91pGeCK1AgWGfHWS6jtYoBm9aIJaKIFOy26PNlmGMukZRUlSCnraBccdcvnw5PB4PWlpaMp7XdCPnppKWlpasxW2n5SAK+hKbSiaqcVMUrXEnsi2GhPloPy24E90gCokQLoEThaL474RP8CHnS7guOBNkMpxDKTFegSmANpUQfMA9zH5OIuyW8bVtINLWLKRxW0fjo50SoZbRXXBIisSQzx7OjEyGVqhAECTMfgcKRZrx55aA/s4haItVSSNYOrs6AACGedW4+nknbEMOaArZv4vCwkIUFRXh1KlTuOOOOzKa13Qj54L7kUceweHDh6FUpuF4YSFlBxFFsTdSmKiN22kHAv7EGnfIUeLu7QYgDMfDRqOUCOATjUKXQOMGmIL0nOCeCJPhHEoJRnAnMzMQPECsAdzja44xOKzJ65QwMBp3qM2ZZcQbTrBJhkYmgnnEC4vPCT8VSMlUwgj3Aa81c8HdPoyyeUVJz+ns7oRSoUS5YQ5MrT3ouNSLlXcn/xGrqanBuXPnMprTdCTnppIdO3agoaEBnZ2dub5UBMYckiiqZKKmEnPoBlMn8GYr1ABfAG8freGpZfHbWIVUAFLsCfeaHItaIM9ZCyiOScbLCO5xIjLEWsCdpqlkbD14NsIaNz2XIYcX6hQEt1ougi9A4oaTzrpMZNobC6OM9Pkyi+UedbjhMI9AX6pJel5HVweKCovB4xHQl2nRc2P8H73Kykr09vZiZGRmmCFzLri3bNmC5uZmVFdXx/SezCmMVj1W4+aJANIHkBNoncYIbpUu/jUeD9DoQfZ3g88jEmo2AkkQ4JGQUokdVmq+HEP+LIQsckw9jKNuPG1VqAI8GZhKUhHcfD4globDVAfsbmgV4nHeRGvcAHDFSkeJpCa4I4WmMqHnOn1vMSntiSBJEp1dHSgpLKGvWayGbcABvy950wSmlklHR0dGc5tu5Fxw37hxI6bv5KT0nGTs2HGmktDziZhLhmn7dbj+9ljUemCoD1q5CPwECQSElNb4ecHETiWNUJ6V7LOppL+/H8899xyOHj061VOZWrw2+nE8+7BImb6N22Zmb1k2FjkdpkqSFAbtHujk4zg0AWjktNLRFta4x7+WmCeEki/NeP12Xe2DTCWFXMUehTM0PAiPx43iIlpwqwsUoCgK5j5b0rELCwvB4/Fw8+bNjOY23ZiU6oDbtm2DXq+HXq/Hj370o9xfcDzBPREH5VAf7c1na9CqVINvHYReyaLViGjBzQ8kvnn0QhX6p7ngdlpHYem3JXytv78fTz/9NA4dOoQdO3bMmDoRGfWcZJJRhEmKTAF0ZEm6f/NUnZMAHbpqt2DQ4YEvQKJQPX4GIWPm6xwdhpIvhWis2ZGFApEaPd40f4RCdFzuQVGFLqljsv2mCQBQUkxr0HKVFDwBH9b+5LtUgUAAvV6Prq74tme5Iq97Tv7whz8ERVEwmUxoa2sDSZL44Q9/mNuLMmUyEzkngYkJ7sEeuk0Z2+JSqiEesUDLotWQIcFNeRObi3RCFfp91qx3EskWXxy5hF//P6/hhZ/8BXuf+wC+UO0Li8WCJ598Et/73vfg9XrxT//0T9DpdHjllVemeMbZIaOek14bbSYhxjENCpXpx3E7LKlr3KGqlTeH6WJWxSkIbqVUCD6PQK/HnJKZhEEnUKLXk34ROfeoFwMdZpTMTZ4xbOo0QafVQSKmPwPBI6DUSGEbHN8vpNfrJzXCKK97Tp45cwY7duyAWq2GRqPBiy++CKPRmNFYKWs9TLW1OBs3YyqZgOAe6AE0SRaXQg25ywYdix2RFNJpxEFv4oCeApEKXtI/LRsq9HcM44M/fop5Kypx5zdq0XGpB03/9SHIIInnn38enZ2daGhowD/+4z+ioKAA9913Hz7//HNcuXJlSuY75UWmfLbx7dsAfU46iSvBIG2zVqYoUJUawDyArpDgLlSNL7h5BAGtXIQBnzWlUECGQpEaXenGpIPuLUlRFEqqxxHcHW0oKYqtva3QyGEZmH6CO5dMSs9JpzMihBwOR8baZMpaD9NRZGyabjZMJYM9gIbFvg0ACjWUXgerxu0l3ABJwOtNrLEXCun48Ey3m7nk8F8+h0qnwJr7l6F6aRnu/vYadF7uxd6X3sHnn3+Or371q6irqwsX81m5ciVKSkrw8ssvT8kOYsp7Tnpt49u3gZDgdqReITBUVmHcAlMMKi0w1I8u8yhkYn64ocd4aBUimIPWuJo6ySgQqdGdwdo1neuCSq+AQs1uVvJ4Pejp68ackti/p1Irg33YCTKY/PsrKChAf39/7n1sk8CkNFLYsGEDnn32WTz77LNYvXo1nnnmmdxelNG4x6YahwX3BPrimQcSR5SEIGVKqIIuaKWJt8dOygVeUBiu0jaWiXrmc8XATTPaL3Rj2Z3zwQvVuSipKsDS2+fhzff3Qa/TY8WKFTHv4fF4ePDBB3HhwgWcPn16KqY9tXhtyZsoMAgVACjAl+Iui2mMMF4TBQZtITDQjT6LC4VJMhLHoleI4SAcKTkmGQqFagx4rQikEblFURTaznWhtLow6XkdN9tBkiTmlI4V3AqQQRL24eTfX0FBAYLBIAYHB5Oelw/kXHA3NDRg586dGB4exvDwMPbu3Ytvfetbub0okzIeJ7hDi3YiGvc4ff7cYtppWUwkjlxxkqPgB4WwjybuNh9OYshCecxscqb5AmRKCeYuit2mdnuuwe4fxhypARQZr1UvXLgQ1dXVeO2116at3T5n+Gzjx3ADEXNKquYSc0jwjFeLm6GwFHA5MdLfl1IoIINeKYZbMJKWjbtQpEYQJPrTiOW2DtjhMI+MK7jb2m9AIpZAr42t063QykAQBMx9yR2UTO/J7u7upOflA5MSVVJbW4sdO3Zgx44dWLVqVe4v6B8BQCSJKsmwOWswSG9Tk2xRR0MaVhGV+BpOchRCUgwLi+AW8YRQ8CUYyEInkWzhcnpw4dg1zF81N6xtA8DA0ADeb34Xq5auBs8lwbG3WuK2qwRB4J577sG1a9emzNY9ZXhtydPdGRjB7U0xft8cqsDHFpI6lmK66qCwp43VhJcItYYAyQ+gQJR6IxQmY7I7jbj09os9IHgEilkaJzDcaL+BkuLSuKgTgZAPuUYGc2/ye0aj0UAoFE5qZEmuyEnKe11dXdKQHoIgcOrUqVxcmsY/QvfyI8b8LjEatz9DU8mInbYtJtG4HTwZigDogokztJzkKMQQw+xgz+DUCJQYnGCVtWxypvkCKCCu/9+hTw5CKpHirjvXwzrgxLmjV3Hqw3O47WsrY85bsGABNBoNDh06hMWLF0/exKcarw2QJk/fBhBb2jUVxgtJHUsRbVpQDnZAMy/1uvhiDb1rlFMpOFhDFIaEfLdnCEBqf+uuK33Ql2ogTFIjnCRJtHe2YdXy2oSvq/RyDPfakl6Hx+OhuLh4crO4c0ROBPdLL72U8LjFYkFjY2PGUSUp4x9JbFsUiCOvZ0K4uzu7zc9K0OYZnc+GRKLZTjohgSRc6zgRWqF82ghuj8uLUx+ew7wVlZDIIzsYv9+Pk60nsXzJCggEAhSWabGwrhqXT7ahelk5iqPCung8HpYtW4bjx4/jRz/6EXi8SdnoTT3jdXhnCAvuFM0L44WkjkUsAbSF0Fu7YJel5pgEAEIxCpAA3KmbV1R8GSQ8Ibq8qUeW3Lzah/L5xUnPGRwagNvjRmnxnISvq/VK9JmGEPQHwReyh18WFxfDZDKlPLfpSk7uoFWrVsX8q6mpwZ49e7Bp0ybU1NTAas3MDJByOKDPmVhwE3xa656o4E6S+GAlxCBBQBXqKxk3BDkCBV8K26gPQTKxF1wtmD6Cu8V4EQFfAEtvnxdz/NLVi/B43Fi8INKotWxeEVR6BS5+dj1unGXLlsFqteLSpUs5nzPD1IcDpii4RQoABJBq/PNgL52hmwakvhglowNQp2Eq8UocAEkgMJr6ewiCQLFIm3JIoMM8AqdlFIVlyc0+ps42AERcKCCDukABiqRgGUhubiorK0N7ezv8/sTBAflCTlUfh8OBJ554AtXV1SAIAu3t7XjhhRcybh6cejjgSMQsMhaBDMi0+h4juJN0HXF4Sdj5MkiciW9CW9ABFV+OIEUljSzpnwbOyWAgiNMfnUf1svK4+shfnG+FXquHXhcRIARBoHx+CXrbh+ByxjqA586dC51OhwMHDkzK3IHshgOmnTlJBukmCqk4Jwk+nT3pTjHaYbAH0KQnuP3qQlT6hsI1SFLBwhsG4ZbB7Ezsj2GjUKTBzRQ/C1MgqqA8uaO1/WY79FodxOLE2r9CLQOPz4N5HHNJeXk5AoHApNQsybvMyWiBrdVqYTKZsH379qx1ex8XNlMJQGtAqTqBxmJj2kWxm0qcbj8cIiVEtviFS1EUHORIuEmwdSTxDaEXqtDrTT/7LNtcb+3EqN2NhWP6/5EkiXMXv4Shal7ce5iU5a6rsS2seDwe7rrrLhw5ciQvnZRpZ06GC0ylGEon0QOuvtTOHQiZStJgVKlHpW8Y6jRMJV2Bfoi8Cgw706uoWSzSwuRO7bN0X++HQiODTJE8TLGt/QZKWMwkQCiDUiuDZZzU9zlz5oDP5+Py5cspzW8i5F3mpFarRXt7O5qamvDoo4+ivb0dZ8+ejfmXU/xOdo1bqEy/LgSDdYi2F4rZF5nd5cMIi+B2kqMIgoQ+1K7KymLn1gvVGPLZ4SOndjt39uPLKCjTxrWRar9pwqhrBDVVNXHvEYoF0BWr4wQ3ANx2220oLS3F73//+5zNedrApLCnLLgLgJEUwtQoChjqTVvjdkp1KPVboRWnaBcHcNPfC3lACbMzvaJsc8Q6tLn7Ugr/7Lzch8JxtG2X24Xe/h6UlSYXgCq9AsM9yXeqQqEQZWVlkyK4c0lOnJMrV66E2WzG1q1bE75OEERuEzL8TjqqJBFidfq1jxksg4Ay+SJzuP0YkahQYYkXXBaS1gYKxSoAw6whgYUiNShQ6PWaUSUtyWyuE8RhGYHpfDfWfGVZ3GvnL52DVCINF/oZS1GFDlfOtMPj8kEStTXn8Xi455578Prrr6Onp2fqMhonAyYcTqxJ7XxpEeDsGP88px3wuNLWuM0SLeaCQrHHDArjx2X7qQD6g8OoJOfCnKbGXSYugCPgwrDfnrShgsvpxsDNYdz+1RWs5wDAddM1UBSF8jnJmymr9Qp0Xe2H1+2DWMpuEqqsrMTFixeTjjXdyYngnvLebj4nvfVMhEiTfrcRhqG+cWNnbaM+uGVaiM0X4t8epLUwvVANmcgGC4uppDi02Ls8Q1MmuL/85Ar4Ah6qlsQL1y8vfIGqymrW6JCiSj2unGlH58UeLFwTa2ZZvHgxRCIRPvvsM2zatCknc58WMMpBqjHQshKg/9j45/WFypLqkkdhjGVAQM9DZe2BvWz8VoK9gUGQoKAjtBhOEgGViHIJ/aNybbQnqeC+3kqH5c2pSR4yefX6FaiUKqhV7GMBgKaY/owDnWZULkqsVABAVVUVjh07hqGhIRQWJk/6ma7MzLisZDZusQZwZZjyOtA1rqZjc/ngUeghsg0Bwdji7v2BYRAgoOYpoZIKYWG5IZhY2EyK9WQDn8ePMwcvoHpZOURj6loMDA2gt78X86rns75fJBGiqEKPq2fagTG7ZZFIhHnz5uU2jn864B4AQKQuuKVFdNz3eBFPvR30oz49wd0F2mQjGUot+aQrQO8YS4Q6WJzetLJeyyW0MLzmSm76OX/sGoorCyAdx759+dolVJRVJs0NAQCpXAyFRoabV3qTnsf0wM3nVmZ5JbhT9uz7nfHp7gwideY27t5OQJdcO7CP+uDTFIAggxBbY80lnf4+6HkaCAg+1DIRhlhshzK+BHK+BD3pdkXJEsffboXH5YsLAQSAz898BpFIhKq5ybW2ysWlcFhGcPNqvJNqwYIFuHTpEkZHM8xgTZEpDQcc6QEkOoCXYrcnaWHkfcnoaqOTb9K0cfePBmAVqyEZSK2RQE9gACIIUSzVwO0PYtSbvMNMNGKeEMUiLW642AWodcCOzsu9qFme3PzhHHGit78HFWWVKV17Tk0Rbl7uxYiVfW0pFAqUlZXhzJkzKY05HckrwZ2yZ9/nYNe4RSpasAfTC3ECRQH9XUk1HZ8/CJcviKCWPkcyGHuT3PB3Yo6AvkHVciGG7OxOnwKhGj1TEFnSfb0fJ947i1vWzYdCE/sder1eHP3sYyxesARCQfLoBE2BEroSNc59ejVO6160aBGCwWDOb5wprQ440pVa1iQDI7hHx3FQmi4DpRWpJ9+E6LO6YZcXQNKfWmOLHv8ACvha6GR0+N1AkrWaCNpByS64vzx6FUKxAJWL2CNFAODaDToCKVXBXT6/GCKpCK1HkjsflyxZghMnTsDlmkDBuSkkrwR3SlBUclMJs3VNt9i7eRDweYECdsFtddE/BpSetktLBmJTa9v93SgV0DezWirEUJK094IpCAkkSRIH/nAUulJ1Qm37yLFDcLndqFu5JqXxalZUwjboQPuFWGGk1WpRUVGBjz/+OBvTnp7Yb9B261SRMBr3OKaMa+eAOVVpTSVIkui1uuFSF0PWeyOl93QF+lHA14Q7OfVZ0xNwJWIdq8ZNURQunriBykWlECTJcgSAK9cvQ6vRQalILTqHL+DDcEsFbl7uTVotsK6uDj6fDx988EFK4043Zp7gDowCVJA98YER3OlGlvSENBU9+83IeN/lShl86gJI+jvCr7lIN6ykA0V82rmplotgd/vgCyTOntSLVJNu475yyoShLgvqNi6Lczz6fF4cPPIhli1eNq6TiEFTqERxpR4thy/BP2arvXr1apw8eTJvSmymnYBjuwIoKlK/AF8EiPWAM0kdjWAQuH4eqIj/UU1Gj9kFX4BEoKQS0t62lJpldwf6UcjXQyYSQCYSoNuSnuAuExfghqs3oW18oNMM26ADcxcn17YpisL5S+cwt6IqrWuXGgohkYtw6QT7j5RGo8Hq1auxZ8+enGndeZeAM6V4mD5/LL/QklBUSKrJDgztoaSRIvZt97CD3k6qZSJ4daUx2k1PgBZQBXxt6Bza1MAWalUq0qWcxJANKIrC5we+RElVAQoTdNk+ceYEXG4X1qxam9a4C1ZXwe/x48ujsUk3tbW1kEgk2L9//4TmPVmklYDjGqCzIJXjR2/EICsBHEnqaLRdokMBqxakNeylHjt4BCCsXgC+zw1ZT3xJgmjcpAcDQTOKBXoQBIESjQRt/el1ZKqQFMIeGE1YuuHyyTaIpaJx25TdaL8Oq92K+Yb0Pi+PR6By0Ry0X+iGw8Lu7K2vr4fL5cI777yT1vipkncJOFMKo6WyhSFJCuiqgck0m0S0HKW3qBL2DtR9NjekIj6kIj68BeWQdV0Nv9YTCKX2hgQ3k3o8zOKgrJQWYdBnw9Ak1SzpuTGAPtMQFq0xxL1GkiQOHzWipmpeyto2g0QuhmF5Oa6ebo/JahOLxbj99ttx4MCBjGvXTFuGz9KPavbIm4TIS2kTCxtfHAP4AqB6YVrDnr9pRZlOjuDcBaB4fKgvnUh6/g0/7ZuZw6fNehV6OS732NOKLDGEwljPOttijlMkhYsnrqNiYUlMieBEHDpqhFajHTd+OxHl80sglolx6oPzrOcwWvf+/fvh9aYX8jjV5JXgTmm76gw5BKUs8Zk8ASCbA1jTyJwKBIBDbwKr7kh6WufQCIpUEhAE4CmqhLS/HTwfXbPD5O+GgpBBwaNt74zg7rclbuqwXEFrax8OT47n+8S7Z6HSK1A2L96hdqr1JPoH+1K2bY9l7qI5UGhlOPHe2Zh63evWrQOfz89ZJuWURZUMtdJ1uOXsscQJkZcnF9wfvwPMWwaI2ZWHsVAUhbMdVhiKFCBFUoxWLoK29VDS95z2XICUkGBOyB9TU6yEecSLjqHUo4DKxAXQCZQ4ZPki5njnlV7Yh0dgWJZcGPcP9uPsuVasXrlm3DDARPAFPCxea0B/xxBuXmZ3kt59991wOp04dCj5dzLdyCvBndJ21XIeEKoAcZIMR5UhohWlwvXzdAOFW25jPYUkKZzrtKGigLatu0sNIMgg5J30D8RF7w1UCiM3skjAh0YmYrUdFojUWCavwt7+T1KfZ4Z0XurBtdYOLLtjXtxNMjQ8iN37/4wFNQvjWkalCsEjsOTWGtgGHbgYZXeUyWR46KGHYDQaU7cdp8GURZUMfwGo58XXgx8PRTltYvHa4l+zW4ETRqBufVpD9lpdsIx4UV1MVyl0LKiD9vxR8DzsQvio+wyWimrAD82/plgJmUiAwxdSN90RBIE7NEuwt/+TGE29xXgRKr0ChRXJE9mMHx+ETCbHkoVLU77mWPRzNCgo16Hl0CUEWfpR6vV6LFmyBG+++SZIlmqd05G8EtwpYT4PqGuSh0upDbSAT5Xzp+h43LnsW9+L3TaYR7xYWq4BAHiK54IUCKG8dgZ+KoDzvmuoEcaGNM3RSnG1h70oznrtLThoboHdn7t4Z/uwE2/97hCKKvSoHqMFBYNBvPSnnZCIJdh4z1cmdB2VXoGqpXNw/tOrMfUk6urqsHHjRrzyyit49913J3SNacNQS/pmEgBQhhpVWBKkYx//kE7oql2X1pBn263g8whUF4UE96K14AV80J09nPD8oYAFJn8XlooiDlABn4dVVVocutDHWoo4ERt0K9HpGcQpO+3fsPTbcfV0OxbWVSfVoh1OB06c+Qwrlq2EgD+x5O4Fqyrhcrpx+WQb6znr169HV1cXPvkk90pStpiBgvtc5AZgQ15OhwOOqRLIqvV9eQKoMCQtLnX4fB90ChGqQjcIJRDCVb4Q2i8/xkXvDXgoLxYIq2LeU12sxOUeO2tkyd265fBRAbw3/Hnyz5MhLqcHr+94DwBw1zdr426mDw8dwM3uTjxY/9WYcppHjia+6cejZnklVHoFPtl3GiP2yE6jvr4ed9xxB3bu3Jn/Re49FtrBqGF3qO1+i2VbrpxLF0frS5D6fvR9oHIe3fg3DU5cH0JVoQKSUNidT1sCd3EVdGc+Snj+Sc85ECCwQBTrWF0+VwfrqA+XkygaY1muNKBAqMIrvc2gSAoH/3QcUqUYNcuTR9u8++Fb4PEIrFi2MuVrsa1JuVqGyoWlOPfpNdiHEzsqq6qqcMstt+DXv/41rl9P7ridLswswe0xA7ZrgGYc540sZLJwxiYjJBTcfj/w6fvA4sQtkwBg1BvA0csDWFmlAy9K9lmX3w3d2cOwt7wOFU+BckFsKOG8YiV8QRIXuhI754pEGiySVeDtwc+Sf54MsA068Oft78Ll9GDDX98ak3bscrtw5NhhvPvR21hbe1tc15GPjx7J6JoEj8Dy9QtBEAQOvnI8Js72a1/7GgoKCrB9+/acZ1TmlP7Q30q3hPWU3W+z/PDxhEDpOuDc/6GTyBgcNuDw22lr25YRL852WHBLpSbmuGNBHfRnDoLwxTrGSYrEe6MfY56wMuyLYagskEMhFuLk9dSzefkEDw8V3orXOg/h1V++g7ZzN7H2geWssdskSeKA8X0cPfEJ1t22HtIkgQBjSbYma1ZUQKYQw/jnz9DfnjjEtqGhAXq9Htu2bcsL4T0jBHdY4HYdBEABReM40RQhm6f1avLzAODI24B1GLjjftZT3j59E/4ghdvmR8KbPjj2Kawr74Vt7iI83PQX1AmXgDdGoy3VSFGqkeK3H12B3ZW4hOucE368O3QS7emGL45h9+7doCgKln4bDr1+Ai8+tQejdjfue+w2qHT0LoEkSRg/OYin/v2f8Jf9f8ayxbfg9jXJHbLpIpaKULdxKXgCHp77t9/gwvHr6DUNASSB73znOxgeHsY//MM/4N1334XTmXoIWi5s5GlDUUDbPjqsT5Y8RpmVJY/TCsjFF+nngQDwnz+iY6/v/nrSt+7+KKLJB4IkfvX+ZYgEPKyYG2tPtq68B3zPKObv/KewrZuiKPzOthuXfG24X3Zn+NyWD+kxeQSBWyo0MJ7vC4e9psJD5Bp878N16LjQg1u/tRLl84rx9vtvxZ137cZV/OzZf8PbB/bj1rrbsWLpypSvMR58AR+19y2BRCGG8fUTaGm+ENfUWiKR4Pvf/z60Wi1+/OMfY9euXejs7Jwe6yoBOakOONnsfv3PeOw2GdDyc4AnAoJeYDR5oRkQfODQ39Dhg0VrgeJQfDJFAS2fAl03gGvngdf+D510I5LQLaNCUKBwo9+JtgEnmj/rwK0lSgjMAxg1A37Kjz3N70Bb7oaxhsSPDzvw61+/hZtVl9BdXo0b8xfDFWrGcGetAntPdOKHfz6MB1aVoVwnw4JSFXgh1X344FV4Vqhxy4kt+Ls5G7FAVo5b1Ytwq3oRCIJAwBdA17V+kEESJBlxAgUDJLwuL0ZsLtiHndix4zmY3qdDukYCNugqFdAa5Pj8i2MYtg6js7sDXb10iOTcsiqsWFoLhUyBwf54DcXj8aK/N8MKiyEqbynEnz86h1OHIzW9lTo5Vs1bi4sdX+K3v/0tfvvb30Kn0aNAVxiuDqfT6SCTyiBViSGWiiCRSCCRSPDCCy9g3bp1qKhII+klG1Ak0HOENo+0vw10vg8YGpLnCQQ8idcnRQFXrwN9c4BXngL27wHODwNtN4H6e4Hhy4CXPTb81bfewwJNEeyjfnxyuR83h0fRsLwU1GAvovcwowDab/82DJ/uRfGnb+B6ZQX2Ly+BRjKC/1e3DOQSIW7w6b/vp+99APXtdGnfuVU8HO+x4bu/P4hSnwAbDEXQS0SgCECvlqBQIYbf4oenzQ1vlwdBO510JQuQODr3DN7/4lMsu1KJT197F1a+Bw6HHS6rHYNdPRh1jUAqkaL+rgdQpC/CQF96iVmprMnyJQXgS0h8+fklfPn5JRRV6qEtUkEoFiDgD0Ig5GPN0jtwXngW+/fvx/79+9Ha2orTp09j0aJFUKvV0Ov1kEql4PP5IAgCJElCIBBAr9fDYIgPpc0VBJVOcOYUw7QdYigrK0NZWRl6rp5CmXhiggSl69Az7EIZnwSunJ3YWAB6/EAZSzkPu4SPozXJ63qHGfYCBfHtmtaqF6FYpEV/xzCcSZIMwtd02aCWaeAjPXAGbKldmwWPxwOJJHlFt8kcJ3qs2tpaDA0NxYQACgQCdHen0KQgCWxrD452wBJfwjcZPTagTJPghREAE6grlmzNpUObXopLJaHM4wTrT+YRo3IoeVQIg5/yw+GPLd2Qzb97LsdMd9xly5ZBoVCgp6cn6+tvLHkluDk4ODg4ZoiNm4ODg2M2wQluDg4OjjyDE9wcHBwcecaMiCpJlX379sFisaCtrQ02mw07d+4EABiNRphMJuh0OphMpnCTY7bjs4UtW7Zw31EGcOsst3DrEgA1S2hra6MaGxvDzxsaGsLP6+vrw8cbGxuppqampMdnA42NjTGfn/uOUoNbZ7mFW5c0s8ZUYrPZsGfPnvDzNWvWoLm5GUajERqNJny8trYWe/bsYT0+GzCZTDGfnfuOUodbZ7mDW5cRZo3grq2tRUtLS/h5W1sbamtr0draCp0uEpOq0+nQ2trKenw2YDQaUV9fH37OfUepw62z3MGtywizRnBHY7PZYDQasW3bNpjNifs6sh2f6RiNRjzyyCMxx7jvKDO4dZY9uHUZy6wU3I8//jiam5uh0Wig1+sTnsN2fKZjs9litpkA+3cxW7+jVOHWWfbg1mUss05wP/PMM2hsbITBYIDJZEJtbS0sFkv4dYvFgtraWtbjM5lnnnkGFosF+/btw759+2AymbBv3z7uO8oAbp1lD25dJmCqvaOTSVNTE9XS0hJ+vnPnToqiKKq2tjZ8LNoDzXZ8NtDS0hLjmee+o9Th1lnu4NYlzaypVWIymVBTUxNzrLGxEVu3boXRaERra2tYO4qOBU10fKZjMpnQ2NgIo9GIxsZGNDQ0cN9RinDrLHdw6zLCrBHcHBwcHDOFWWfj5uDg4Mh3OMHNwcHBkWdwgpuDg4Mjz+AENwcHB0eewQluDg4OjjyDE9wcHBwceQYnuPMEo9GIp556CiaTKSfj22y2hMd37dqFp556KifX5MgfuPU3veAEd57Q3NwcTqE2Go2oqanJWsUzm82GXbt2JXxt8+bNs6L2A0dyuPU3veAEdx5SX1+f1doLe/fuRUNDQ9bG45jZcOtv6uEE9wyEbdvJRltbGwwGQ24mwzHr4NZf7plVPSenii1btuDMmTPh2gkajQaHDh2KK1OZKVqtFvX19XjppZdw5swZbNmyBU899RQMBgOam5uxZs0aaDSacHeWpqam8HtbW1uxZs0aAHQtCKPRCIPBgNbW1qxrVhxTA7f+Zh6c4M4xmzZtCjc3feaZZwAgqwVvjEYjXnrppfBWs76+HvX19bDZbKivr4fBYMDq1avR3t4OjUaD5uZm7Nu3L3z+nj170NjYCADYuXMnHn30UdTW1sJgMKStOXFMP7j1NzPhTCU5hHHeMO2Wamtrs9qdg/H0J7IPMltPnU4HnU4X1q40Gk1MreJotmzZgvvuuw8bN26E0WjktJ08h1t/MxdOcOcQo9GIjRs3hp8z28ZsYTAYUFdXl9Ajn8o2eN++fXj00UfDz3U6Hdrb27FlyxY0NTWxevo58gNu/c1cOMGdQzQaTUzT0tbW1rB2woRAMYtzy5YtaY9fW1uLnTt3orGxMS6+NpVt5unTp2O0mu3bt0Oj0aChoQFNTU1oa2tLe04c0wdu/c1cOBt3Dtm8eXM4ecBkMsU4ZSwWS7jIOwCsXr065XGZIvGNjY3YuXMnDAYDNm3ahJdeein8usViQX19PbZv3x5u+8TE4JpMJqxduzYuPlav12PXrl3heW3btm2iXwHHFMKtvxnM1Dbgmd1s3ryZslqtFEVRMa2uErF169asXruxsTF87VTO5Zh5cOsvf+E07ilk06ZNMBqN0Gg0YQfSZGE2m7MWDsaRn3DrL3/hbNxTSH19PRoaGlK+abJVK8JkMqXspNq1axf27Nkz4WtyTD+49Ze/cD0nOTg4OPIMTuPm4ODgyDM4wc3BwcGRZ3CCm4ODgyPP4AQ3BwcHR57BCW4ODg6OPIMT3BwcHBx5Bie4OTg4OPIMTnBzcHBw5Bmc4Obg4ODIM/5/7Qq0ff8cqEEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 350x477 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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\", \"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": "markdown",
   "metadata": {},
   "source": [
    "## Void testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evidence comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"no_field\", \"IndranilVoid_exp\"]\n",
    "cats = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "neglnZ = {}\n",
    "kfound = []\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        sample_alpha = sim not in [\"IndranilVoid_exp\", \"no_field\"]\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"mike\",\n",
    "            sample_alpha=sample_alpha, zcmb_max=zcmb_max)\n",
    "        \n",
    "\n",
    "        neglnZ[f\"{sim}_{cat}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simA = sims[0]\n",
    "simB = sims[1]\n",
    "\n",
    "print(f\"lnZ_({simA}) - lnZ_({simB})\\n\")\n",
    "for cat in cats:\n",
    "    lnZ_A = - neglnZ[f\"{simA}_{cat}\"]\n",
    "    lnZ_B = - neglnZ[f\"{simB}_{cat}\"]\n",
    "    print(f\"{cat:15s} {lnZ_A - lnZ_B:.1f}\")\n",
    "\n",
    "\n",
    "print(f\"\\n(Positive -> preference for {simA})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Goodness-of-fit comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "no_Vext = True\n",
    "\n",
    "sims = [\"IndranilVoid_exp\", \"IndranilVoid_gauss\", \"IndranilVoid_mb\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "neglnZ = {}\n",
    "kfound = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        kfound[f\"{sim}_{cat}\"] = []\n",
    "        for ksim in range(500):\n",
    "            sample_alpha = False\n",
    "            fname = paths.flow_validation(\n",
    "                fdir, sim, cat, inference_method=\"mike\", nsim=ksim,\n",
    "                sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "                sample_beta=True,\n",
    "                no_Vext=no_Vext, verbose_print=False)\n",
    "\n",
    "            if not exists(fname):\n",
    "                continue\n",
    "\n",
    "            kfound[f\"{sim}_{cat}\"].append(ksim)\n",
    "            neglnZ[f\"{sim}_{cat}_{ksim}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "\n",
    "neglnZ_no_field = {}\n",
    "neglnZ_dipole = {}\n",
    "sim = \"no_field\"\n",
    "for cat in cats:\n",
    "    sample_alpha = False\n",
    "    fname = paths.flow_validation(\n",
    "        fdir, sim, cat, inference_method=\"mike\",\n",
    "        sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "        no_Vext=True, verbose_print=False)\n",
    "\n",
    "    if not exists(fname):\n",
    "        continue\n",
    "\n",
    "    neglnZ_no_field[f\"{cat}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "    fname = paths.flow_validation(\n",
    "        fdir, sim, cat, inference_method=\"mike\",\n",
    "        sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "        no_Vext=None, verbose_print=False)\n",
    "\n",
    "    if not exists(fname):\n",
    "        continue\n",
    "\n",
    "    neglnZ_dipole[f\"{cat}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    figwidth = 8.3 \n",
    "    fig, axs = plt.subplots(2, 2, figsize=(figwidth, 0.65 * figwidth))\n",
    "\n",
    "    for n, cat in enumerate(cats):\n",
    "        i, j = n // 2, n % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            x = kfound[f\"{sim}_{cat}\"]\n",
    "            y = [neglnZ[f\"{sim}_{cat}_{ksim}\"] / np.log(10) for ksim in x]\n",
    "            x = np.array(x) * 0.674\n",
    "            ax.plot(x, y, label=simname_to_pretty(sim))\n",
    "        \n",
    "        # if no_Vext is None:\n",
    "        #     y_no_field = neglnZ_no_field[cat] / np.log(10)\n",
    "        #     if cat != \"CF4_TFR_w1\":\n",
    "        #         ax.axhline(y_no_field, color=\"black\", ls=\"--\", label=\"No peculiar velocity\")\n",
    "        y_no_field = neglnZ_no_field[cat] / np.log(10)\n",
    "        ax.axhline(y_no_field, color=\"black\", ls=\"--\", label=\"No peculiar velocity\")\n",
    "\n",
    "        y_dipole = neglnZ_dipole[cat] / np.log(10)\n",
    "        ax.axhline(y_dipole, color=\"black\", ls=\":\", label=\"Constant dipole\")\n",
    "\n",
    "        ax.text(0.5, 0.9, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, #fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='center',\n",
    "            bbox=dict(facecolor='white', alpha=0.5),\n",
    "            )\n",
    "\n",
    "        if n == 0:\n",
    "            ax.legend(fontsize=\"small\", loc=\"upper left\")\n",
    "\n",
    "        ax.set_ylabel(r\"$-\\Delta \\log \\mathcal{Z}$\")\n",
    "        ax.set_xlabel(r\"$R_{\\rm offset} ~ [\\mathrm{Mpc} / h]$\")\n",
    "        ax.set_xlim(0)\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fname = f\"../../plots/void_goodness_of_fit_observer.png\"\n",
    "    if no_Vext:\n",
    "        fname = fname.replace(\".png\", \"_no_Vext.png\")\n",
    "    print(f\"Saving to `{fname}`.\")\n",
    "    fig.savefig(fname, dpi=450)\n",
    "    fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Single parameter radial dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "key = \"beta\"\n",
    "# key_label = r\"$\\sigma_v ~ [\\mathrm{km} / \\mathrm{s}]$\"\n",
    "# key_label = r\"$|\\mathbf{V}_{\\rm ext}| ~ [\\mathrm{km} / \\mathrm{s}]$\"\n",
    "key_label = r\"$\\beta$\"\n",
    "no_Vext = True\n",
    "\n",
    "sims = [\"IndranilVoid_exp\", \"IndranilVoid_gauss\", \"IndranilVoid_mb\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "data_mean = {}\n",
    "data_std = {}\n",
    "kfound = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        kfound[f\"{sim}_{cat}\"] = []\n",
    "        for ksim in range(500):\n",
    "            sample_alpha = False\n",
    "            fname = paths.flow_validation(\n",
    "                fdir, sim, cat, inference_method=\"mike\", nsim=ksim,\n",
    "                sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "                sample_beta=True,\n",
    "                no_Vext=no_Vext, verbose_print=False)\n",
    "\n",
    "            if not exists(fname):\n",
    "                continue\n",
    "\n",
    "            kfound[f\"{sim}_{cat}\"].append(ksim)\n",
    "            with File(fname, 'r') as f:\n",
    "                x = f[f\"samples/{key}\"][...]\n",
    "            if key == \"Vext\":\n",
    "                x = np.linalg.norm(x, axis=-1)\n",
    "\n",
    "            data_mean[f\"{sim}_{cat}_{ksim}\"] = x.mean()\n",
    "            data_std[f\"{sim}_{cat}_{ksim}\"] = x.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    figwidth = 8.3\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(figwidth, 0.65 * figwidth))\n",
    "\n",
    "    for n, cat in enumerate(cats):\n",
    "        i, j = n // 2, n % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            x = kfound[f\"{sim}_{cat}\"]\n",
    "            y = [data_mean[f\"{sim}_{cat}_{ksim}\"] for ksim in x]\n",
    "            yerr = [data_std[f\"{sim}_{cat}_{ksim}\"] for ksim in x]\n",
    "            x = np.array(x) * 0.674\n",
    "\n",
    "            ax.plot(x, y, label=simname_to_pretty(sim))\n",
    "            ax.fill_between(x, np.array(y) - np.array(yerr), np.array(y) + np.array(yerr), alpha=0.5)\n",
    "\n",
    "        ax.text(0.5, 0.9, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, #fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='center',\n",
    "            bbox=dict(facecolor='white', alpha=0.5),\n",
    "            )\n",
    "\n",
    "        if n == 0:\n",
    "            ax.legend(fontsize=\"small\", loc='upper right')\n",
    "\n",
    "        ax.set_ylabel(key_label)\n",
    "        ax.set_xlabel(r\"$R_{\\rm offset} ~ [\\mathrm{Mpc} / h]$\")\n",
    "        ax.set_xlim(0),\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fname = f\"../../plots/void_{key}_per_observer.png\"\n",
    "    if no_Vext:\n",
    "        fname = fname.replace(\".png\", \"_no_Vext.png\")\n",
    "    print(f\"Saving to `{fname}`.\")\n",
    "    fig.savefig(fname, dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}