diff --git a/maf/science/sigma8tomography_demo.ipynb b/maf/science/sigma8tomography_demo.ipynb new file mode 100644 index 0000000..59cd3fd --- /dev/null +++ b/maf/science/sigma8tomography_demo.ipynb @@ -0,0 +1,763 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.4.1.dev0+ge6735eb76.d20251021\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/metrics/base_metric.py:32: UserWarning: Redefining metric TotalPowerMetric! (there are >1 metrics with the same name)\n", + " warnings.warn(\"Redefining metric %s! (there are >1 metrics with the same name)\" % (metricname))\n", + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/metrics/base_metric.py:32: UserWarning: Redefining metric StaticProbesFoMEmulatorMetricSimple! (there are >1 metrics with the same name)\n", + " warnings.warn(\"Redefining metric %s! (there are >1 metrics with the same name)\" % (metricname))\n", + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/metrics/base_metric.py:32: UserWarning: Redefining metric NestedRIZExptimeExgalM5Metric! (there are >1 metrics with the same name)\n", + " warnings.warn(\"Redefining metric %s! (there are >1 metrics with the same name)\" % (metricname))\n" + ] + } + ], + "source": [ + "import os\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import healpy as hp\n", + "import pandas as pd\n", + "import rubin_sim\n", + "import rubin_sim.maf as maf\n", + "#from rubin_sim.data import get_baseline\n", + "print(rubin_sim.__version__)\n", + "from os.path import splitext, basename\n", + "from rubin_scheduler.scheduler.utils import SkyAreaGenerator\n", + "\n", + "from rubin_sim.data import get_baseline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# fixing a weird issue\n", + "#import os\n", + "#os.environ['RUBIN_SIM_DATA_DIR'] = '~/rubin_sim_data'\n", + "#os.environ['RUBIN_SIM_DATA_DIR'] " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demo of Tomographic Sigma8 bias Metric" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "from rubin_sim.maf.metrics.uniformity_metrics import NestedLinearMultibandModelMetric\n", + "from rubin_sim.maf.metrics.cosmology_summary_metrics import TomographicClusteringSigma8biasMetric" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from rubin_sim.maf.metrics.tomography_models import DENSITY_TOMOGRAPHY_MODEL\n", + "# this contains the current model.\n", + "# the first set of keys are the years (year1, ..., year10) since this would change typical depth and galaxy catalog cuts.\n", + "# in what follows we have 5 tomographic bins.\n", + "# the second nested dictionary has the following:\n", + "# sigma8square_model is the fiducial sigma8^2 value used in CCL for the theory predictions\n", + "# poly1d_coefs_loglog is a polynomial (5th degree) describing the angular power spectra (in log log space) in the 5 tomographic bins considered, thus has shape (5, 6)\n", + "# lmax contains the lmax limits to sum the Cells over when calculating sigma8 for each tomographic bin. thus is it of shape (5, )\n", + "# dlogN_dm5 contains the derivatives of logN wrt m5 calculated in Qianjun & Jeff's simulations. It is an array of 5 dictionaries (5 = the tomographic bins)\n", + "# each dictionary must have keys that are the lsst bands. If some are missing they are ignored in the linear model.\n", + "# they are the ones which will be fed to LinearMultibandModelMetric. Everything else above is going into the modeling \n", + "# The notebook I used to make this dictionary is https://github.com/ixkael/ObsStrat/blob/meanz_uniformity_maf/code/meanz_uniformity/romanrubinmock_for_sigma8tomography.ipynb" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# a simple wrapper around the metrics, to store the results, but not critically needed\n", + "def extract_sigma8_tomography_metric(\n", + " opsim_fname, run_fname,\n", + " years, \n", + " percentage_uncorrected,\n", + " density_tomography_model, \n", + " lmin = 10,\n", + " mag_range_tolerated=1.0,\n", + " n_filters = 6, \n", + " extinction_cut = 0.2, # sky cuts\n", + " nside=32,\n", + " convert_to_sigma8=True\n", + "):\n", + " surveyAreas = SkyAreaGenerator(nside=nside)\n", + " map_footprints, map_labels = surveyAreas.return_maps()\n", + " slicer = maf.HealpixSubsetSlicer(\n", + " nside=nside, \n", + " hpid=np.where(map_labels == \"lowdust\")[0],#hpid=np.arange(hp.nside2npix(nside)), \n", + " use_cache=False)\n", + " \n", + " # prepare empty arrays to fill in the results\n", + " n_bins = 5 # set to 5\n", + " results_spuriousdensitypower = np.zeros((len(years), n_bins))\n", + " results_sigma8_squared_bias = np.zeros((len(years), ))\n", + " # loop over years\n", + " all_depth_map_bundles = []\n", + " for iy, year in enumerate(years):\n", + " print('year', year)\n", + " \n", + " # constraints\n", + " days = year*365.25\n", + " constraint_str = 'scheduler_note not like \"DD%\" and night <= XX and scheduler_note not like \"twilight_near_sun\" '\n", + " constraint_str = constraint_str.replace('XX','%d'%days)\n", + " \n", + " # leave empty if not specified\n", + " mean_depth = {}\n", + " min_depth_cut = {} \n", + " max_depth_cut = {}\n", + " \n", + " ##############################\n", + " # now converts depth fluctuations to density fluctuations\n", + " ##############################\n", + " metric = NestedLinearMultibandModelMetric(\n", + " density_tomography_model['year'+str(year)]['dlogN_dm5'], \n", + " extinction_cut=extinction_cut, n_filters=n_filters, # cuts going into ExgalM5WithCuts\n", + " mean_depth=mean_depth, min_depth_cut=min_depth_cut, max_depth_cut=max_depth_cut,\n", + " )\n", + " # summary metric measures total power via angular power spectra of healpix map (thus needs nside)\n", + " # _but_ has a bin-dependent lmax to consider same scales to consider the same scales as a fct of redshift\n", + " summary_metrics = [\n", + " TomographicClusteringSigma8biasMetric(\n", + " density_tomography_model['year'+str(year)], convert_to_sigma8=convert_to_sigma8,\n", + " power_multiplier=percentage_uncorrected, lmin=lmin),\n", + " ]\n", + " # then standard way of packing MetricBundles into a MetricBundleGroup\n", + " depth_map_bundles = [maf.MetricBundle(\n", + " metric=metric,\n", + " slicer=slicer,\n", + " constraint=constraint_str,\n", + " run_name=run_name,\n", + " summary_metrics=summary_metrics\n", + " )]\n", + " bd = maf.metricBundles.make_bundles_dict_from_list(depth_map_bundles)\n", + " bgroup = maf.MetricBundleGroup(bd, opsim_fname)\n", + " bgroup.run_all()\n", + "\n", + " # compute bias\n", + " # should probably also return fsky\n", + " results_sigma8_squared_bias[iy] = depth_map_bundles[0].summary_values['TomographicClusteringSigma8bias']\n", + " all_depth_map_bundles.append(depth_map_bundles[0])\n", + "\n", + " return results_sigma8_squared_bias, all_depth_map_bundles" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'/Users/yoachim/rubin_sim_data/sim_baseline/baseline_v5.0.0_10yrs.db'" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "baseline_file = get_baseline()\n", + "baseline_file" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "run_name: baseline_v5.0.0_10yrs\n", + "Healpix slicer using NSIDE=64, approximate resolution 54.967783 arcminutes\n", + "year 1\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/maps/dust_map.py:46: UserWarning: Slicer value of nside 64 different from map value 128, using slicer value\n", + " warnings.warn(\n", + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/metrics/cosmology_summary_metrics.py:73: HealpyDeprecationWarning: \"verbose\" was deprecated in version 1.15.0 and will be removed in a future version. \n", + " data = hp.remove_monopole(data, verbose=False, bad=self.mask_val)\n", + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/metrics/cosmology_summary_metrics.py:75: HealpyDeprecationWarning: \"verbose\" was deprecated in version 1.15.0 and will be removed in a future version. \n", + " data = hp.remove_dipole(data, verbose=False, bad=self.mask_val)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "spuriousdensitypowers: [0.00568099 0.00622485 0.00399589 0.00128107 0.016763 ]\n", + "fskys: [0.40405994 0.40405994 0.40405994 0.40405994 0.40405994]\n", + "0.6400000000000001 0.6623500832062111 0.005453857495016782 4.098032122517393\n", + "0.8 0.8138489314401114 0.003350657157813148 4.1331985899596\n", + "year 2\n", + "spuriousdensitypowers: [0.02120177 0.01728666 0.00693721 0.0046805 0.00184558]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6500411741131265 0.005413417437544279 1.8548678776343113\n", + "0.8 0.8062513095264568 0.0033571526480520027 1.8620867687038734\n", + "year 3\n", + "spuriousdensitypowers: [0.00661991 0.00819832 0.00337406 0.00193849 0.00055764]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6439435189551981 0.005392792880939481 0.7312572617309481\n", + "0.8 0.8024609142850498 0.003360159220804027 0.7323802603797281\n", + "year 4\n", + "spuriousdensitypowers: [0.00316623 0.0044001 0.00232425 0.00151433 0.00107296]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6432237706359281 0.005373997667820726 0.5998831475554484\n", + "0.8 0.8020123257381573 0.003350321120610334 0.6006366750273329\n", + "year 5\n", + "spuriousdensitypowers: [0.00426516 0.00698654 0.00229943 0.00166337 0.00146881]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6439688491092934 0.005365226082378959 0.7397356697284772\n", + "0.8 0.8024766969260188 0.0033429170609757812 0.7408789631459689\n", + "year 6\n", + "spuriousdensitypowers: [0.00312654 0.00468064 0.00204403 0.00143537 0.00144758]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6433850445700525 0.005378806347953882 0.6293300689919952\n", + "0.8 0.8021128627381888 0.0033528986990634603 0.630160028031551\n", + "year 7\n", + "spuriousdensitypowers: [0.00268651 0.00407459 0.00239563 0.00092605 0.00152737]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.643165557951361 0.005358097025898052 0.5907989228377768\n", + "0.8 0.8019760332774047 0.0033405593207077038 0.5915276717750584\n", + "year 8\n", + "spuriousdensitypowers: [0.00290831 0.00571121 0.00246014 0.00084173 0.00107826]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6428582456260991 0.0053571332017877645 0.5335401451554515\n", + "0.8 0.8017844134342467 0.0033407566373370866 0.5341345174035117\n", + "year 9\n", + "spuriousdensitypowers: [0.00252918 0.00447619 0.00254608 0.00087256 0.00129307]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6429274281752003 0.0053559232541604464 0.5465776928237704\n", + "0.8 0.801827555135891 0.0033398223968822967 0.5472012935768552\n", + "year 10\n", + "spuriousdensitypowers: [0.00290333 0.00386475 0.00230076 0.00070763 0.00093397]\n", + "fskys: [0.40415986 0.40415986 0.40415986 0.40415986 0.40415986]\n", + "0.6400000000000001 0.6423523213141789 0.00534966111623347 0.4397140796527548\n", + "0.8 0.8014688523668146 0.003337410493517875 0.4401173813252805\n" + ] + } + ], + "source": [ + "#baseline_file = get_baseline()\n", + "sim_list = [baseline_file\n", + "]\n", + "name_list = [splitext(basename(sim))[0] for sim in sim_list]\n", + "\n", + "years = range(1, 11) # [1, 2, 4, 7, 10]#\n", + "percentage_uncorrected = 0.1\n", + "\n", + "# large mag_range_tolerated and no min depth in order to make comparison fair between strategies etc\n", + "results_fsky = {}\n", + "results_sigma8_squared_bias = {}\n", + "for opsim_fname, run_name in zip(sim_list, name_list):\n", + " \n", + " print('run_name:', run_name)\n", + " results_sigma8_squared_bias[run_name], _ = extract_sigma8_tomography_metric(\n", + " opsim_fname, run_name,\n", + " years, \n", + " percentage_uncorrected,\n", + " DENSITY_TOMOGRAPHY_MODEL,\n", + " nside=64,\n", + " lmin=10, n_filters=6, extinction_cut=0.2, mag_range_tolerated=2.0,\n", + " convert_to_sigma8=True\n", + " )\n", + " \n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'Top panel minus baseline_v5.0.0_10yrs')" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnoAAAJaCAYAAAC4Buo6AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAc05JREFUeJzt3Xl4TGf/BvD7ZE9kkSAbEVFLRMSSUBFBX8RWS6soam0ttYb2LaotRSltUfWzVqutavUt2rRFpZbYYgshiNQSEiRiy0JIYub8/hgZRraZZGbOnMn9ua65MnPmmTPfkYjbc55FEEVRBBERERGZHQupCyAiIiIiw2DQIyIiIjJTDHpEREREZopBj4iIiMhMMegRERERmSkGPSIiIiIzxaBHREREZKYY9IiIiIjMlJXUBZgDpVKJGzduwMnJCYIgSF0OERERmTlRFJGTkwNvb29YWJTcb8egpwc3btyAj4+P1GUQERFRJZOamopatWqV+DyDnh44OTkBUP1hOzs7S1wNERERmbvs7Gz4+PioM0hJGPT0oPByrbOzM4MeERERGU1ZQ8Y4GYOIiIjITDHoEREREZkpBj0iIiIiM8UxenKgVAC39gMP0wB7L6BGOGBhKXVVRETlIooiHj9+DIVCIXUpRCbL0tISVlZWFV62jUHP1KVuAeImA7nXnh5zqAUEfwn4vCpdXURE5ZCfn4+0tDTk5uZKXQqRyXNwcICXlxdsbGzKfQ4GPVOWugXY/xoAUfN47nXV8fBfGfaISDaUSiWSk5NhaWkJb29v2NjYcJF5omKIooj8/HzcunULycnJqF+/fqmLIpeGQc9UKRWqnrznQx7w5JgAxEUCNXvzMi4RyUJ+fj6USiV8fHzg4OAgdTlEJs3e3h7W1ta4evUq8vPzYWdnV67zyHYyRmpqKq5de3o58+jRo4iMjMSaNWskrEqPbu3XvFxbhAjkpqraERHJSHl7JogqG338XZHt37ZBgwZhz549AID09HR07twZR48exfvvv485c+ZIXJ0ePEzTbzsiIiKqdGQb9M6cOYNWrVoBAH755RcEBgbi0KFD2LhxI9avXy9tcfpg76XfdkRERFTpyDboFRQUwNbWFgDwzz//oFevXgAAf39/pKWZQS9XjXDV7FqUNFBZABx8VO2IiCoTpQK4uRe48pPqq9Lwy7R06NABkZGRBn+fkgwfPhx9+vQxmXpIPmQb9Bo3boxVq1Zh//79iI6ORteuXQEAN27cQLVq1SSuTg8sLFVLqAAoMewFL+VEDCKqXFK3AFF1gF0vAYcGqb5G1VEdr0S2bNmCuXPnSl2G2pUrVyAIQpHbjh07Sn3dvXv3MGTIELi4uMDFxQVDhgxBZmZmqa8RRRGzZ8+Gt7c37O3t0aFDB5w9e1brWtesWYMOHTrA2dkZgiAU+37lqctUyTboLVy4EKtXr0aHDh0wcOBANG3aFAAQFRWlvqRbHgsWLIAgCKbxPyWfV1VLqDjU1Dxuac+lVYio8ilccur5iWqFS05VorDn5uYGJycnqcso4p9//kFaWpr69p///KfU9oMGDUJ8fDx27NiBHTt2ID4+HkOGDCn1NYsWLcLixYuxfPlyHDt2DJ6enujcuTNycnK0qjE3Nxddu3bF+++/r9e6ylJQUFCh15ebKGOPHz8W7969q3EsOTlZvHnzZrnOd/ToUbFOnTpiUFCQOHnyZK1fl5WVJQIQs7KyyvW+ZVI8FsX0PaJ4dpEo/ghR/NFCFHNvGOa9iIgM5OHDh+K5c+fEhw8fPj2oVIpiwf2yb3lZoril5pPfgcXdBFHcUkvVTpvzKZU61d6+fXtx/Pjx4vjx40UXFxfRzc1NnDlzpqh8cp4ffvhBDA4OFh0dHUUPDw9x4MCBGv8W3b17Vxw0aJBYvXp10c7OTqxXr574zTffqJ+/du2a2L9/f7Fq1aqim5ub2KtXLzE5OVn9/LBhw8TevXtr1PPsv1O+vr7iJ598Io4YMUJ0dHQUfXx8xNWrV2t8hrLeoyQ7duwQbW1txXv37mkcnzhxotiuXTtRFFX/9gIQT548Web5Cp07d04EIB4+fFh9LDY2VgQgnj9/vtjXKJVK0dPTU/z000/Vxx49eiS6uLiIq1at0vq9RVEU9+zZIwIo8rnKqkupVIovvPCC+Nlnn2m8LiEhQRQEQbx48aIoiqIIQFy5cqXYq1cv0cHBQfzoo4/K/Dl4XrF/Z57QNnvItkcPUG0P4urqqnGsTp06cHd31/lc9+/fx+DBg7F27doi55SchSXg0QEI+C9QIwyAErj0jdRVERFVnCIX+MWx7NuvLsDD66WcSAQeXlO10+Z8Ct135vjuu+9gZWWFI0eOYNmyZViyZAm+/vprAKo1AufOnYtTp07ht99+Q3JyMoYPH65+7Ycffohz585h+/btSExMxMqVK1G9enUAqh6ml156CY6Ojti3bx8OHDgAR0dHdO3aFfn5+VrX98UXXyAkJAQnT57EuHHj8Pbbb+P8+fMVfo9OnTqhatWq2Lx5s/qYQqHAL7/8gsGDB2u07dWrF9zd3REWFoZff/211PPGxsbCxcUFL774ovpY69at4eLigkOHDhX7muTkZKSnpyMiIkJ9zNbWFu3bty/xNboqqy5BEDBy5Eh8++23Gq/75ptvEB4ejhdeeEF9bNasWejduzcSEhIwcuTIUn8ODEXWCyb/+uuv+OWXX5CSklLkB/XEiRM6nWv8+PHo0aMHOnXqhHnz5pXaNi8vD3l5eerH2dnZOr1XhbwwGrh1ELi0Fmg8AxBkndWJiGTDx8cHS5YsgSAIaNiwIRISErBkyRKMGjUKI0eOVLerW7culi1bhlatWuH+/ftwdHRESkoKmjdvjpCQEACqTolCP//8MywsLPD111+rdwr59ttvUbVqVezdu1cj1JSme/fuGDduHABg2rRpWLJkCfbu3Qt/f/8KvYelpSUGDBiAjRs34s033wQA7Nq1C/fu3UO/fv0AAI6Ojli8eDHCwsJgYWGBqKgoDBgwAN999x3eeOONYs+bnp5ebMeMu7s70tPTS3wNAHh4eGgc9/DwwNWrV0v749GaNnWNGDECH330EY4ePYpWrVqhoKAAGzZswGeffabxmkGDBmn8bJT2c2Aosg16y5Ytw8yZMzFs2DD8/vvvGDFiBC5duoRjx45h/PjxOp3r559/xokTJ3Ds2DGt2i9YsAAff/xxecquuNr9VDtmPLgKpO0EvLtKUwcRkT5YOgD975fdLmMfsLd72e06bAPc22n3vjpq3bq1xpZtoaGh+OKLL6BQKHD69GnMnj0b8fHxuHv3LpRKJQDVP+wBAQF4++230bdvX5w4cQIRERHo06cP2rRpAwCIi4vDxYsXi4y5e/ToES5duqR1fUFBQer7giDA09MTGRkZenmPwYMHIzQ0FDdu3IC3tzd+/PFHdO/eXX0FrHr16pgyZYq6fUhICO7du4dFixaVGPQK63yeKIplbo33/PPavEYXZdXl5eWFHj164JtvvkGrVq3w559/4tGjR+rgW6gw0BUq7efAUGTbHbRixQqsWbMGy5cvh42NDd577z1ER0dj0qRJyMrK0vo8qampmDx5MjZs2KD19iIzZsxAVlaW+paamlrej6E7K3vAb6jq/kUz2QWEiCovQQCsqpR984zQbskpzwjtzqfHUPDo0SNERETA0dERGzZswLFjx7B161YAUF9t6tatG65evYrIyEjcuHEDHTt2xLvvvgtAtQdwcHAw4uPjNW7//vsvBg0apHUd1tbWGo8FQVAHzoq+R6tWrfDCCy/g559/xsOHD7F169ZSAxygCsYXLlwo8XlPT0/cvHmzyPFbt24V6bF79jUAivT4ZWRklPgaXWlb11tvvaX+8/j2228xYMCAIlv7ValSReNxaT8HhiLboJeSkqJOwfb29urZNkOGDMFPP/2k9Xni4uKQkZGB4OBgWFlZwcrKCjExMVi2bBmsrKygUBRdn8nW1hbOzs4aN6OqN1r19XoUd8Ygosqh1CWnnjw28JJThw8fLvK4fv36OH/+PG7fvo1PP/0U4eHh8Pf3V/ekPatGjRoYPnw4NmzYgKVLl6q37GzRogUuXLgAd3d31KtXT+Pm4uKil9r18R6DBg3Cjz/+iD/++AMWFhbo0aNHqe1PnjwJL6+SF/UPDQ1FVlYWjh49qj525MgRZGVlldjL5efnB09PT0RHR6uP5efnIyYmRm89Y9rW1b17d1SpUgUrV67E9u3bNS7RlqaknwNDkW3Q8/T0xJ07dwAAvr6+6r+AycnJUE120U7Hjh2RkJCg8T+ckJAQDB48GPHx8bC0NMF16qo2Vk3KEBWclEFElUdJS0451DLKklOpqamYOnUqkpKS8NNPP+Grr77C5MmTUbt2bdjY2OCrr77C5cuXERUVVWSNu48++gi///47Ll68iLNnz+LPP/9Eo0aNAKgui1avXh29e/fG/v37kZycjJiYGEyePFljT/eK0Md7DB48GCdOnMAnn3yC1157TeMq2HfffYeNGzciMTERSUlJ+Pzzz7Fs2TJMnDhR3ebo0aPw9/fH9euqSTWNGjVC165dMWrUKBw+fBiHDx/GqFGj8PLLL6Nhw4bq1/n7+6t7SAuXP5s/fz62bt2KM2fOYPjw4XBwcNC69zM9PR3x8fG4ePEiAKgzwN27d3Wqy9LSEsOHD8eMGTNQr149hIaGlvnepf0cGEypc3JN2JtvvinOnj1bFEVRXLlypWhvby926tRJrFq1qjhy5MgKnfv5aetlMfjyKsW5/L1qSYHffEVRqTDe+xIRlVNpS0XopHDJqeSNqq+Kx/oor1Tt27cXx40bJ44dO1Z0dnYWXV1dxenTp6uXV9m4caNYp04d0dbWVgwNDRWjoqI0lhuZO3eu2KhRI9He3l50c3MTe/fuLV6+fFl9/rS0NHHo0KFi9erVRVtbW7Fu3briqFGj1P+uaLO8ypIlSzRqbtq0qThr1iyt30MbLVu2FAGIu3fv1ji+fv16sVGjRqKDg4Po5OQkBgcHiz/88INGm8LlTJ5d0uXOnTvi4MGDRScnJ9HJyUkcPHhwkeVOAIjffvut+rFSqRRnzZolenp6ira2tmK7du3EhIQErT/DrFmzRABFbs++hzZ1iaIoXrp0SQQgLlq0qMhzAMStW7dqHCvr5+B5+lheRXhSjOwolUoolUpYWanmk/zyyy84cOAA6tWrh7Fjx8LGxqbc5+7QoQOaNWuGpUuXatU+OzsbLi4uyMrKMt5l3McPga3eQEEm0GE7J2UQkcl79OgRkpOT4efnp/WYaCJTdvDgQXTo0AHXrl3T2xjBZ5X2d0bb7CHbWbcWFhawsHh65bl///7o37+/Xs69d+9evZzHoAonZfy7TDUpg0GPiIjIKPLy8pCamooPP/wQ/fv3N0jI0xfZBj1AlXRPnz6NjIwM9cyiQr169ZKoKiOqN1oV9AonZdiXPOiViIioJI6OjiU+t337doSHhxuxmvL78ccfMWbMmGKf8/X11WlP3NL89NNPePPNN9GsWTP88MMPejmnocg26O3YsQNDhw7F7du3izwnCEKxs2XNTuGkjFsHVZMyAmdKXREREclQfHx8ic/VrFmzxOdMTa9evTR2tHjW88vPVMTw4cM1dj4xZbINehMmTEC/fv3w0UcfmXSXqcHVG8OdMoiIqELq1asndQl64eTkVGRR6MpOtqkgIyMDU6dOrdwhDwB8XgNsXJ/ulEFEZOJkOgeQyOj08XdFtkHvtddek8ekCUPjThlEJBOFl85yc3MlroRIHgr/rlTksrNsl1fJzc1Fv379UKNGDTRp0qTIH8KkSZOMVosky6s8K/MssC0QECyB3imAg7fxayAi0kJaWhoyMzPh7u4OBwcHve5PSmQuRFFEbm4uMjIyULVq1WJ3GDH75VU2btyIv//+G/b29ti7d6/GLwtBEIwa9CT37KSMy99yUgYRmazCvUqL2yKMiDRVrVpV/XemvGTbo+fp6YlJkyZh+vTpGuvpSUHyHj0ASP4BiB0KVPEFel3mpAwiMmkKhQIFBQVSl0FksqytrUvdhtXse/Ty8/MxYMAAyUOeyfB5DYib/HRSBhdQJiITZmlpaZp7iROZGdmmpGHDhmHTpk1Sl2E6OCmDiIiIniPbHj2FQoFFixbh77//RlBQUJHJGIsXL5aoMgnVGw0kfanaKSP3BidlEBERVXKyDXoJCQlo3rw5AODMmTMaz1XaWVwuAUCNtsCtA5yUQURERPINenv27JG6BNNUb7Qq6F1aCwRMByw4BoaIiKiyku0YPSrBsztlpEdLXQ0RERFJSLY9eq+88kqxl2gFQYCdnR3q1auHQYMGoWHDhhJUJ6HCSRlJX6omZXD2LRERUaUl2x49FxcX7N69GydOnFAHvpMnT2L37t14/PgxNm3ahKZNm+LgwYMSVyqBeqNVXwsnZRAREVGlJNug5+npiUGDBuHy5cvYvHkztmzZgkuXLuGNN97ACy+8gMTERAwbNgzTpk2TulTjK5yUISpUkzKIiIioUpLtzhg1atTAwYMH0aBBA43j//77L9q0aYPbt28jISEB4eHhyMzMNGgtJrEzxvOe3Smj5yVOyiAiIjIj2mYP2fboPX78GOfPny9y/Pz581AoFAAAOzu7yrvUCidlEBERVXqynYwxZMgQvPnmm3j//ffRsmVLCIKAo0ePYv78+Rg6VLVDRExMDBo3bixxpRLhpAwiIqJKT7aXbhUKBT799FMsX74cN2/eBAB4eHhg4sSJmDZtGiwtLZGSkgILCwvUqlXLoLWY5KVbAMg6B/zVGBAsgd4p3CmDiIjITGibPWQb9J6VnZ0NAJKFLJMNegAQHa5aQDloHnfKICIiMhNmP0bvWc7OzqYXsExF4VIrl9YCSoW0tRAREZFRyWqMXosWLbBr1y64urqiefPmpU60OHHihBErM2E+rwFxk59OyuBYPSIiokpDVkGvd+/esLW1BQD06dNH2mLkQmNSxmoGPSIiokrELMboSc2kx+gBnJRBRERkZsx+jF5qaiquXbumfnz06FFERkZizZo1ElZlorhTBhERUaUk26A3aNAg7NmzBwCQnp6OTp064ejRo3j//fcxZ84ciaszQZyUQUREVOnINuidOXMGrVq1AgD88ssvaNKkCQ4dOoSNGzdi/fr10hZnirhTBhERUaUj26BXUFCgnpjxzz//oFevXgAAf39/pKWlSVmaaSqclAGoJmUQERGR2ZNt0GvcuDFWrVqF/fv3Izo6Gl27qmaT3rhxA9WqVZO4OhNVePn2+h9A7g1payEiIiKDk23QW7hwIVavXo0OHTpg4MCBaNq0KQAgKipKfUmXnqMxKeMbqashIiIiA5P18ioKhQLZ2dlwdXVVH7ty5QocHBzg7u5utDpMfnmVZyX/AMQOBar4Aj0vARaWUldEREREOjL75VUAwNLSUiPkAUCdOnWMGvJkh5MyiIiIKg1ZBz0qB07KICIiqjQY9CojTsogIiKqFBj0KiNOyiAiIqoUZBX03NzccPv2bQDAyJEjkZOTI3FFMlZvjOrrpa+5UwYREZGZklXQy8/PR3Z2NgDgu+++w6NHjySuSMZ8+j4zKWOn1NUQERGRAVhJXYAuQkND0adPHwQHB0MURUyaNAn29vbFtv3mG16SLFXhpIykL4GLawDvblJXRERERHomqx69DRs2oHv37rh//z4EQUBWVhbu3btX7I20wEkZREREZk22Cyb7+fnh+PHjJrHdmawWTH5edDhw6wAQNBcI/EDqaoiIiEgLZr9gcnJyskmEPNnjpAwiIiKzJdugBwAxMTHo2bMn6tWrh/r166NXr17Yv3+/1GXJCydlEBERmS3ZBr0NGzagU6dOcHBwwKRJkzBhwgTY29ujY8eO2Lhxo9TlyYfGThlrpK2FiIiI9Eq2Y/QaNWqE0aNHY8qUKRrHFy9ejLVr1yIxMdFotch6jB4AZJ0D/moMCJZA7xTAwVvqioiIiKgUZj9G7/Lly+jZs2eR47169UJycrIEFckYd8ogIiIyS7INej4+Pti1a1eR47t27YKPj48EFckcJ2UQERGZHVktmPysd955B5MmTUJ8fDzatGkDQRBw4MABrF+/Hl9++aXU5cmPT18gbtLTSRlcQJmIiEj2ZNuj9/bbb+Pnn39GQkICIiMjMXnyZJw5cwabNm3CmDFjdDrXypUrERQUBGdnZzg7OyM0NBTbt283UOUmysoe8Bumus9JGURERGZBtpMx9OmPP/6ApaUl6tWrB0C1j+5nn32GkydPonHjxmW+XvaTMQpxUgYREZEsmP1kDH3q2bMnunfvjgYNGqBBgwb45JNP4OjoiMOHD0tdmnFxUgYREZFZYdB7jkKhwM8//4wHDx4gNDS02DZ5eXnIzs7WuJkNTsogIiIyGwx6TyQkJMDR0RG2trYYO3Ystm7dioCAgGLbLliwAC4uLuqbWc3y5U4ZREREZoNB74mGDRsiPj4ehw8fxttvv41hw4bh3LlzxbadMWMGsrKy1LfU1FQjV2tAnJRBRERkNsxmMoZCoUBCQgJ8fX3h6upa4fN16tQJL7zwAlavXl1mW7OZjFEoKxH4K4CTMoiIiEyU2U/GiIyMxLp16wCoQl779u3RokUL+Pj4YO/evRU+vyiKyMvLq/B5ZMmlEVAjnJMyiIiIZE62Qe/XX39F06ZNAaiWR0lOTsb58+cRGRmJmTNn6nSu999/H/v378eVK1eQkJCAmTNnYu/evRg8eLAhSpeHeqNVXy+u5aQMIiIimZJt0Lt9+zY8PT0BANu2bUO/fv3QoEEDvPnmm0hISNDpXDdv3sSQIUPQsGFDdOzYEUeOHMGOHTvQuXNnQ5QuD4WTMnJTOCmDiIhIpmS7BZqHhwfOnTsHLy8v7NixAytWrAAA5ObmwtLSUqdzFV4CpmcUTspIWqqalMEt0YiIiGRHtj16I0aMQP/+/REYGAhBENS9b0eOHIG/v7/E1ZmJwsu31/8Acm9IWwsRERHpTLY9erNnz0ZgYCBSU1PRr18/2NraAgAsLS0xffp0iaszE4WTMm7tV03KCPxA6oqIiIhIB7Lt0fv+++/Rs2dPTJkyBbVq1VIfHzhwILKysiSszMxwUgYREZFsyTbojRgxothAl5OTgxEjRkhQkZmq/RonZRAREcmUbIOeKIoQBKHI8WvXrsHFxUWCisyUpR13yiAiIpIp2Y3Ra968OQRBgCAI6NixI6ysnn4EhUKB5ORkdO3aVcIKzVC90arZt4WTMrhTBhERkSzILuj16dMHABAfH48uXbrA0dFR/ZyNjQ3q1KmDvn37SlSdmeKkDCIiIlmSXdCbNWsWAKBOnToYMGAA7OzsJK6okqg3WhX0Lq4FAmYAFrqtVUhERETGJ9sxesOGDWPIMyZOyiAiIpIdWQU9Nzc33L59GwDg6uoKNze3Em+kZxqTMlZLWwsRERFpRVaXbpcsWQInJycAwNKlS6UtpjJST8r4E8i9DjjUlLoiIiIiKoUgiqIodRFyl52dDRcXF2RlZcHZ2Vnqcgwrup1qrF7QXE7KICIikoi22UNWPXrPUyqVuHjxIjIyMqBUKjWea9eunURVmTlOyiAiIpIN2Qa9w4cPY9CgQbh69Sqe75QUBAEKBbfrMojarwFxk55OyvDuJnVFREREVAJZTcZ41tixYxESEoIzZ87g7t27uHfvnvp29+5dqcszX5yUQUREJBuy7dG7cOECfv31V9SrV0/qUiofTsogIiKSBdn26L344ou4ePGi1GVUToU7ZYgK4NI3UldDREREJZBtj97EiRPxzjvvID09HU2aNIG1tbXG80FBQRJVVkkUTsq49DXQ+H1OyiAiIjJBsl1excKiaGekIAgQRdHokzEq1fIqhRSPgK3eQP49oMM2TsogIiIyIrNfXiU5OVnqEiq3wkkZSUtVkzIY9IiIiEyObIOer6+v1CUQJ2UQERGZNNkGve+//77U54cOHWqkSiqxwkkZt/arJmU0+VDqioiIiOgZsh2j5+rqqvG4oKAAubm5sLGxgYODg1HX0quUY/QKJf8IxL4BONQGel3mpAwiIiIj0DZ7yHZ5lWcXSL537x7u37+PpKQktG3bFj/99JPU5VUetfsCNq5Pd8ogIiIikyHboFec+vXr49NPP8XkyZOlLqXy4E4ZREREJsusgh4AWFpa4saNG1KXUbnUG636Wjgpg4iIiEyCbCdjREVFaTwWRRFpaWlYvnw5wsLCJKqqkuKkDCIiIpMk26DXp08fjceCIKBGjRr4z3/+gy+++EKaoiqzemO4UwYREZGJkW3QUyqVUpdAz6rdF4ibpJqUkfY3ULO71BURERFVemY3Ro8k8uykjEtrpK2FiIiIADDokT7VG6X6ykkZREREJoFBj/SncFKGqFBNyiAiIiJJMeiRftUbo/p66WtAqZC2FiIiokqOQY/0q3ZfwMbt6aQMIiIikoxsZ90CQGZmJo4ePYqMjIwis3CHDh0qUVWVXOGkjKQlqkkZnH1LREQkGUEURVHqIsrjjz/+wODBg/HgwQM4OTlBEAT1c4Ig4O7du0arRduNhSuNrETgrwBAsAR6XwUcakpdERERkVnRNnvI9tLtO++8g5EjRyInJweZmZm4d++e+mbMkEfF4KQMIiIikyDboHf9+nVMmjQJDg4OUpdCxeGkDCIiIsnJNuh16dIFx48fl7oMKgknZRAREUlOtpMxevTogf/+9784d+4cmjRpAmtra43ne/XqJVFlBICTMoiIiEyAbCdjWFiU3BkpCAIUCuNdLuRkjBJwUgYREZFBmP1kDKVSWeLNmCGPSuHSCHBvx0kZREREEpFt0COZeGG06uultZyUQUREZGSyGqO3bNkyjB49GnZ2dli2bFmpbSdNmmSkqqhUtfsCcZOA3FTVpAyO1SMiIjIaWY3R8/Pzw/Hjx1GtWjX4+fmV2E4QBFy+fNlodXGMXhnipqomZdTqDbT7TepqiIiIZE/b7CGrHr3k5ORi75OJqzdKFfSu/wnkXuekDCIiIiPhGD0yPE7KICIikgSDHhkHJ2UQEREZHYMeGYd6p4xU7pRBRERkJJU+6C1YsAAtW7aEk5MT3N3d0adPHyQlJUldlvkp3CkDUO2UQURERAZX6YNeTEwMxo8fj8OHDyM6OhqPHz9GREQEHjx4IHVp5qfeKNXXwkkZREREZFCyWl6lUEZGBs6ePYvg4GA4Ozvj5s2b+O6776BUKtGjRw80adKk3Oe+desW3N3dERMTg3bt2mn1Gi6vooN/2gMZ+4Amc4AmH0pdDRERkSyZ7RZoe/fuRd26ddGxY0f4+/vj9OnTCAkJwddff43169ejZcuW2LlzZ7nPn5WVBQBwc3PTV8n0LE7KICIiMhrZBb0PPvgAw4cPR3Z2NqZOnYoePXqgd+/e+Pfff3H+/HlMnDgRH3/8cbnOLYoipk6dirZt2yIwMLDEdnl5ecjOzta4kZY4KYOIiMhoZBf0EhISMGXKFDg6OiIyMhLp6el466231M+PHj0aZ8+eLde5J0yYgNOnT+Onn34qtd2CBQvg4uKivvn4+JTr/SqlZydlXFwtbS1ERERmTnZBz8bGBo8ePQIA5OfnQ6lUqh8DwMOHD2Ftba3zeSdOnIioqCjs2bMHtWrVKrXtjBkzkJWVpb6lpqbq/H6VWr0nl29vcFIGERGRIcku6IWFhWH69Ok4ePAgpkyZghYtWmDevHl48OABcnNzMXfuXISEhGh9PlEUMWHCBGzZsgW7d+8udQ/dQra2tnB2dta4kQ5c/J/slKHkThlEREQGJKu9bgHgs88+Q/fu3REeHo6AgADs3LkTb7/9NqpWrQoAcHV1xY4dO7Q+3/jx47Fx40b8/vvvcHJyQnp6OgDAxcUF9vb2hvgIBKgmZWTsAy6uAaq3AfIyAHsvoEY4YGEpdXVERERmQZbLqwDAnTt3UK1aNfXjXbt24eHDhwgNDdU4XhZBEIo9/u2332L48OFanYPLq5SD4hGwuQbw+L7mcYdaQPCXgM+r0tRFREQkA9pmD9n16BV6Psx17NixXOeRac6VvxvbioY8QDVmb/9rQPivDHtEREQVJNugV+jcuXNISUlBfn6+xvFevXpJVBGVSakA4iaX8KQIQADiIoGavXkZl4iIqAJkG/QuX76MV155BQkJCRAEQd0zV3gpVqHgYrwm69Z+IPdaKQ1E1Tp7t/YDHh2MVRUREZHZkd2s20KTJ0+Gn58fbt68CQcHB5w9exb79u1DSEgI9u7dK3V5VJqHafptR0RERMWSbY9ebGwsdu/ejRo1asDCwgIWFhZo27YtFixYgEmTJuHkyZNSl0glsffSbzsiIiIqlmx79BQKBRwdHQEA1atXx40bNwAAvr6+SEpKkrI0KkuNcNXsWhQ/4xkQAAcfVTsiIiIqN9kGvcDAQJw+fRoA8OKLL2LRokU4ePAg5syZg7p160pcHZXKwlK1hAqA4sOeCDT/nBMxiIiIKki2Qe+DDz6AUqkEAMybNw9Xr15FeHg4tm3bhmXLlklcHZXJ51XVEioONZ974knwu7nb6CURERGZG9kumFycu3fvwtXVtcRFkA2FCyZXgFKhml37ME01Jq8gB9jXG4AItFwF1B8jdYVEREQmx+wXTC6Om5ub1CWQriwsiy6h0vQT4NT7QNxEoGogUCNMktKIiIjkTrZB79GjR/jqq6+wZ88eZGRkqC/jFjpx4oRElVGFBUwH7p4AUn9V7ZLR9Xgxl3iJiIioLLINeiNHjkR0dDRee+01tGrVyuiXa8mABAFo/S2QkwRkJgD7XwU6xQCWdlJXRkREJCuyHaPn4uKCbdu2ISxM+st6HKNnIPcvAztCgPx7QN2RwItfq0IgERFRJadt9pDtrNuaNWvCyclJ6jLIkBzrAmGbAMECuPwNcGGF1BURERHJimyD3hdffIFp06bh6tWrUpdChuTVGWi2UHU/LhK4GSNpOURERHIi26AXEhKCR48eoW7dunBycoKbm5vGjcyI/zuA70BAfAwc6Ac8SJG6IiIiIlmQ7WSMgQMH4vr165g/fz48PDw4GcOcCYJqfF52InAvHtj3CtD5AGBlL3VlREREJk22Qe/QoUOIjY1F06ZNpS6FjMHKAWj3m2pyxr0TwNHRQOj3nJxBRERUCtleuvX398fDhw+lLoOMqYov0PYXQLAErmwAkpZKXREREZFJk23Q+/TTT/HOO+9g7969uHPnDrKzszVuZKY8XgJaLFbdP/kukL5L2nqIiIhMmGzX0bOwUGXU58fmiaIIQRCgUCiMVgvX0TMyUQQOjwCSvwNsqwFdjgGOflJXRUREZDRmv9ftnj17pC6BpCIIQKtVQNY54O4x1eSMiIOAVRWpKyMiIjIpsu3RMyXs0ZNI7jVgRzDwKAOoPQAI+4mTM4iIqFIw+x6906dPF3tcEATY2dmhdu3asLW1NXJVZFQOtYC2m4FdLwEpmwC3FkDAe1JXRUREZDJkG/SaNWtW6tp51tbWGDBgAFavXg07OzsjVkZG5d4WCPkKOPY2ED8dqBoEeHeVuioiIiKTINtZt1u3bkX9+vWxZs0axMfH4+TJk1izZg0aNmyIjRs3Yt26ddi9ezc++OADqUslQ6s3BnjhLQAicHAgkHNR6oqIiIhMgmzH6LVq1Qpz585Fly5dNI7//fff+PDDD3H06FH89ttveOedd3Dp0iWD1sIxeiZAkQf80wG4cxhwaQxExALWTlJXRUREZBDaZg/Z9uglJCTA19e3yHFfX18kJCQAUF3eTUtLM3ZpJAVLWyB8M2DvBWSdBWKHAaJS6qqIiIgkJdug5+/vj08//RT5+fnqYwUFBfj000/h7+8PALh+/To8PDykKpGMzcEbCN8CWNgA17YCZ+dLXREREZGkZDsZ4//+7//Qq1cv1KpVC0FBQRAEAadPn4ZCocCff/4JALh8+TLGjRsncaVkVNVbAyH/BxwdBZz+CHBtBtR8WeqqiIiIJCHbMXoAcP/+fWzYsAH//vsvRFGEv78/Bg0aBCcn447N4hg9E3RsPHBhBWDtDHQ5Cjg3lLoiIiIivdE2e8g66JkKBj0TpMgHdncCbu1XhbyII4CNi9RVERER6YVZLpgcFRWFbt26wdraGlFRUaW27dWrl5GqIpNkaQO0/R/wdwiQnQTEDgHa/QYIsh2WSkREpDNZ9ehZWFggPT0d7u7usLAo+R9sQRCgUCiMVhd79EzYneNAdFtAmQcEfgQEfSx1RURERBVmlsurKJVKuLu7q++XdDNmyCMTVy0EaLVGdf/MHCB1q7T1EBERGZGsgl5ZMjMzpS6BTFHdoUDDyar7sUOBzLPS1kNERGQksg16CxcuxKZNm9SP+/XrBzc3N9SsWROnTp2SsDIySc0/AzxeAh7fB/b1AfLvSV0RERGRwck26K1evRo+Pj4AgOjoaPzzzz/YsWMHunXrhv/+978SV0cmx8IaCNsEVPEF7l8EDg4ClLzET0RE5k22QS8tLU0d9P7880/0798fEREReO+993Ds2DGJqyOTZFcDCN8KWNoDaTuA0x9KXREREZFByTboubq6IjU1FQCwY8cOdOrUCQAgiiInY1DJ3JoDL65T3T+3ALj6i7T1EBERGZBsg96rr76KQYMGoXPnzrhz5w66desGAIiPj0e9evUkro5MWp2BQKMnl/cPjwDunZa2HiIiIgORbdBbsmQJJkyYgICAAERHR8PR0RGA6pIu97elMjVdAHhGAIpc1eSMvDtSV0RERKR3slow2VRxwWSZyrsL/N0SuH8Z8OwEdNgOWMhqsxgiIqqkzHLBZCK9snVTbYtmVQVI/weIny51RURERHrFoEeVW9UmQOvvVPfPfwEk/yhtPURERHrEoEdUuy/QeKbq/tG3gLsnpK2HiIhITxj0iACgyceAd3dA8QjY9wrw6JbUFREREVWY7Eee5+fnIyMjA0qlUuN47dq1JaqIZMnCEmjzI/B3KyDnAnCgP/CfnaodNYiIiGRKtj16Fy5cQHh4OOzt7eHr6ws/Pz/4+fmhTp068PPzk7o8kiObqkC73wErJyBjL3DiXakrIiIiqhDZ9ugNHz4cVlZW+PPPP+Hl5QVBEKQuicyBSyOgzQ+qtfX+XabaSaPucKmrIiIiKhfZBr34+HjExcXB399f6lLI3NTqDTSZDSTMBo6OBZwDgOqtpK6KiIhIZ7K9dBsQEIDbt2/r5Vz79u1Dz5494e3tDUEQ8Ntvv+nlvCRjgR+qAp8yD9j/KvAwXeqKiIiIdCbboLdw4UK899572Lt3L+7cuYPs7GyNmy4ePHiApk2bYvny5QaqlmRHsABCvwecGwEPrwMHXgMU+VJXRUREpBPZboFmYaHKqM+PzRNFEYIgQKFQlOu8giBg69at6NOnj9av4RZoZiz7X9VM3IIsoN5YoNVKqSsiIiLSOnvIdozenj17JHvvvLw85OXlqR/r2oNIMuLcAGizEYh5Gbi4SjU5o95oqasiIiLSimyDXvv27SV77wULFuDjjz+W7P3JyGp2B5p+Apx6Hzg+AXAJBGq0kboqIiKiMsnq0u3p06cRGBgICwsLnD59utS2QUFB5XoPbS7dFtej5+Pjw0u35kwUVYsop/4K2HkCXeMAB2+pqyIiokrKLC/dNmvWDOnp6XB3d0ezZs0gCAKKy6kVGaOnDVtbW9ja2hrs/GSCBAFo/S2QkwRkJqhm4naKASz5c0BERKZLVkEvOTkZNWrUUN8nMiprRyB8K/B3S+DOEeDYOODFr1UhkIiIyATJKuj5+voWe7+i7t+/j4sXL6ofJycnIz4+Hm5ubtwzlzQ5vQCE/Qzs7QZc/gZwCwYajJO6KiIiomLJdh09fTp+/DiaN2+O5s2bAwCmTp2K5s2b46OPPpK4MjJJXhFAs4Wq+3GTgYx90tZDRERUAllNxjBVXEevEhJF4NBg4OpPgG0N1eSMKj5SV0VERJWEttmDPXpE5SEIqvF5rs2AvFvA/leAxw+lroqIiEgDgx5ReVk5AO1+A2yrA3fjgKOjVT19REREJkK2QS81NRXXrl1TPz569CgiIyOxZs0aCauiSqeKL9D2F0CwBK5sAJK+lLoiIiIiNdkGvUGDBqm3QUtPT0fnzp1x9OhRvP/++5gzZ47E1VGl4vES0GKx6v7Jd4H03dLWQ+WnVAA39wJXflJ9VRpuPU4iImOQbdA7c+YMWrVqBQD45ZdfEBgYiEOHDmHjxo1Yv369tMVR5dNgIuA3DBAVwMH+QM4lBga5Sd0CRNUBdr0EHBqk+hpVR3WciEimZLWO3rMKCgrUu1P8888/6NWrFwDA398faWlpUpZGlZEgAK1WAVnngLvHgD8bAWLB0+cdagHBXwI+r0pXI5UsdQuw/zUAz42xzL2uOh7+K793RCRLsu3Ra9y4MVatWoX9+/cjOjoaXbt2BQDcuHED1apVk7g6qpQs7YB6o1T3nw15wNPAwN4h06NUqNZDfD7kAU+PxUWyV5aIZEm2QW/hwoVYvXo1OnTogIEDB6Jp06YAgKioKPUlXSKjUiqAMyWND2VgMFkZ+4Dca6U0EIHcVODWfqOVRESkL7JeMFmhUCA7Oxuurq7qY1euXIGDgwPc3d2NVgcXTCYAqrF4u14qu13gh0CtPoBTfcDaydBVUSFFPnD/EpB9HshOevo18zSgyC379VaOgFMD1UzrKnWefH3mZuPKfY+JyGi0zR6yDnqmgkGPAKgmXhwapNtr7DxUgU/j1gBwqqdap4909+i2KsTlJGmGuvuXVZNlDMXKqWj4ezYQ2nkwCBKR3mibPWQ7GQMAfv31V/zyyy9ISUlBfn6+xnMnTpyQqCqqtOy9tGvn3BjIy1DtqPHopup260Ax56tZTAisDzi9oBoPWJkpC1TB7fneuezzQP7dkl9n5Qg4NwSc/QGnhoCLP+BYD4jpCTy8geLH6QmAvTfQ/g/VJdwHVzVvuVeBRxnA4xwg64zqVhwLW6BK7eJ7A6v4qr7fFgb8laxUqC4/P0xT/azWCAcsLA33fkRkEmQb9JYtW4aZM2di2LBh+P333zFixAhcunQJx44dw/jx46UujyqjGuGq2bW511FiYHCoBXQ/pfoHNj8LyLnw3O1f1df8e8DD66pbxt6i56lSu/gQWMUPsLQx/Gc1lrw7RYNcTpJq+Rrxccmvq+KrCnLO/k+DnXNDVWArrlctZNmTWbcCNL93wtPn3ZqrbsV5nAs8SHka/B5cBe5feXo/9zqgzHv6fS6OYKn6+ajiCzg8CX+OdZ55XBuwtC3zj6xYqVtUE06eHYvImeBElYJsL936+/tj1qxZGDhwIJycnHDq1CnUrVsXH330Ee7evYvly5cbrRZeuiU19TIdQLGBQdtlOvLuFBMCn9wKskt+nWD5JOQ0KCYE+uqnx0jfPUPKx09655KKXnLNu13y6ywdigY5Z3/VZy3PZe9iw5APELy04mFIWaA674MrRXsEH1xR9RQqC8o6C2DnWUKP4JNj1o7Ff67ilo7R9WfSVLGnkiopsx+j5+DggMTERPj6+sLd3R3R0dFo2rQpLly4gNatW+POnTtGq4VBjzQYMjCIouoyYUkhsLRJBRbWqh6/wuDn/EwYdPABBC0m4VekZyj/HpBV3Ni5S6WHHAef4gOdfU39j3mTKjQoFcCj9KIB8NnH2kwYsXHTDH8OPsDZ+UB+Sb8Pn/Qy90qWZzhiTyVVYmY/Rs/T0xN37tyBr68vfH19cfjwYTRt2hTJycmQaXYlc+HzKlCzt2ECgyAA9h6qm3tbzedEUfV+z14CVt8uPrl0+K/q9jwLW9XYv+ImhhRe7tRmUeGavVQBpbixc3m3Sv5clg6q4Fk4dk4d6BoAVlUq+qemPQtLwKOD8d7v2fd1qKm61WhT9HlRVPXyPh/+cp9cIn5wFSjIVI1PzL8L3Dup5Rs/WTpm10uqYGhVpfibpUMJxwu/2hl/ogkXuSbSimx79N566y34+Phg1qxZWLVqFaZOnYqwsDAcP34cr776KtatW2e0WtijRyZPVKp6PQqDX/aTIHj/guqyaWk9apYOgOMLwP2LgOJhye2EJ/9vLG3snEOt4sfOOdTSrkeRSlaQXbRHMGM/cOeIEd5cUF0ufzb86RoWrRyKOfbk9vyQA6VCtT1diesfyrynsjLgJfcKM/tLt0qlEkqlElZWql8Av/zyCw4cOIB69eph7NixsLEx3oB0Bj2SNeVjIDcFyC5mUsiDK7ovSWJp9yTMPT92rkHxY8jIcLRd27HBZFVv4uMHgOKB6uvjB6pJJkWOPXNT5hn8IwAALGw0w6KoVP3HoyzNFgHu7QGbqqqbdVV5TVYy1zDES+56YfZBz5Qw6JHZUhYA95OBS18DiZ+V3T74S6DBBPbOmQp1z1cZM8HL2/OlVKjGDj4fAIsEw9ySw6L6NblFzyEqK/gHUAxLe9Xi1oXBT/3VVTMQFt63cX2mnYthl8B5lrmGIXOfHGREZjlG7/Tp0wgMDISFhQVOnz5datugoCAjVUVkxiysVePkvLtrF/SqBjHkmRILS1UwKG3pmOCl5e8lsrAELJwMs8OLKKp6DIsLg7digVPTyz5HlRdU+04XZD6dra54CDx8+GTdxHKwciwmAFZ9euz5oPhsiLR21n7SkzmOPyxzX2lBtU1kzd7m0XNpImTVo2dhYYH09HS4u7vDwsICgiAUO/FCEAQoFMbbT5Q9emT2DN0zRIZlyJngUijPz6NSATzOBvIzn9zuPZnA8uRW8MzxZx8Xfn18Xw+FC6qw93xQ1OhddFHtmZ1/r+Rz2LkD7X5XPRQfq4ZfiI9Vwyyeva9+rqT72rZ7cl98rPpzVN/Xsl1hXQXZqrVBy9L0U6Dmy6q1I7lNZInM8tLt1atXUbt2bQiCgKtXr5ba1tfX10hVMehRJaGvNQJJGuY23svYP4/KAtUi588HwLJCY+H90iYyUclsXAGH2qrQp/H1ySLidp7y/jmuALMMeqaKQY8qDXPrGSJ5k9PPoyIPKMgqvdcwPxO4Fw/cOVz2+WyqqXoHLaxUM94Fy5LvC1ZPHpdxv6zzlPecFk/u30sATkwq+7NVqataJqggs+y2gtWTHWWeC4DPhkIznQRm9kHvzp07qFatGgAgNTUVa9euxcOHD9GrVy+Eh4cbtRYGPapUzK1niOTN3H4etZ0p3XGPNGs+VoSul9wLsoEHqapVAdRbDD65n5uiCvjarApg4/pkK8ESegbtPfU7tthIP5NmG/QSEhLQs2dPpKamon79+vj555/RtWtXPHjwABYWFnjw4AF+/fVX9OnTx2g1MegREZFemPt4WH1eclcqgEdpT9aNTHkmEBbev6rqRS2LhTVgX0szBFap/XSPaQcf7XsFjThb2myDXrdu3WBlZYVp06Zhw4YN+PPPPxEREYGvv/4aADBx4kTExcXh8GEtur71hEGPiIj0xtzHwxrzknt+lmr3F40g+EzP4MPrWvYKupV+edjeE7j2m1GXjjHboFe9enXs3r0bQUFBuH//PpydnXH06FGEhIQAAM6fP4/WrVsjMzPTaDUx6BERkV7JafxheZjKJXflY1UNJV0efpCiXa+gYAVALCU06r8n1izX0QOAu3fvwtPTEwDg6OiIKlWqwM3NTf28q6srcnJypCqPiIio4gy5Z7YpkGpf6edZWAFVfFS3GmHFt8nPKhr+nr08/PB66Vs/AlDvK31rv9E/t+yCHqBaJ6+0x0RERLJnKmGosrNxAWyaAFWbFP+88jFwYSUQp8WM4odp+q1NC7IMesOHD4etrS0A4NGjRxg7diyqVKkCAMjLM9Lei0REREQWViWHwOfZexm2lmLILugNGzZM4/Ebb7xRpM3QoUONVQ4RERFVdjXCVWPwypotXcO4y78BMgx63377rdQlEBERET1l6H2lK1Ka0d+RiIiIyNz4vKpaQsWhpuZxh1qSLokjux49IiIiIpNkgrOlGfSIiIiI9MXEZkvz0i0RERGRmWLQIyIiIjJTDHpEREREZopj9PSgcLvg7OxsiSshIiKiyqAwcxRmkJIw6OlB4d66Pj4+EldCRERElUlOTg5cXFxKfF4Qy4qCVCalUokbN27AycmJ++6WQ3Z2Nnx8fJCamgpnZ2epyyEt8HsmT/y+yQ+/Z/JkjO+bKIrIycmBt7c3LCxKHonHHj09sLCwQK1ataQuQ/acnZ35i0xm+D2TJ37f5IffM3ky9PettJ68QpyMQURERGSmGPSIiIiIzBSDHknO1tYWs2bNgq2trdSlkJb4PZMnft/kh98zeTKl7xsnYxARERGZKfboEREREZkpBj0iIiIiM8WgR0RERGSmGPSIiIiIzBSDHhEREZGZYtAjIiIiMlMMekRERERmikGPiIiIyEwx6BERERGZKQY9IiIiIjPFoEdERERkphj0iIiIiMwUgx4RERGRmWLQIyIiIjJTDHpEREREZopBj4iIiMhMMegRERERmSkGPSIiIiIzxaBHREREZKYY9IiIiIjMFIMeERERkZli0CMiIiIyUwx6RERERGaKQY+IiIjITDHoEREREZkpBj0iIiIiM8WgR0RERGSmGPSIiIiIzBSDHhEREZGZqnDQy87Oxm+//YbExER91ENEREREeqJz0Ovfvz+WL18OAHj48CFCQkLQv39/BAUFYfPmzXovkIiIiIjKR+egt2/fPoSHhwMAtm7dClEUkZmZiWXLlmHevHl6L5CIiIiIykfnoJeVlQU3NzcAwI4dO9C3b184ODigR48euHDhgt4LJCIiIqLy0Tno+fj4IDY2Fg8ePMCOHTsQEREBALh37x7s7Oz0XiARERERlY+Vri+IjIzE4MGD4ejoCF9fX3To0AGA6pJukyZN9F0fEREREZWTIIqiqOuL4uLikJKSgs6dO8PR0REA8Ndff6Fq1aoICwvTe5FEREREpDudgl5BQQEaNmyIP//8EwEBAYasi4iIiIgqSKcxetbW1sjLy4MgCIaqh4iIiIj0ROfJGBMnTsTChQvx+PFjQ9RDRERERHqi8xi9V155Bbt27YKjoyOaNGmCKlWqaDy/ZcsWvRZIREREROWj86zbqlWrom/fvoaohYiIiIj0qFyzbkmTUqnEjRs34OTkxPGLREREZHCiKCInJwfe3t6wsCh5JJ7OPXqzZ8/GiBEj4OvrW6ECzcmNGzfg4+MjdRlERERUyaSmpqJWrVolPq9zj15wcDBOnTqF9u3b480338Srr75a6XfEyMrKQtWqVZGamgpnZ2epyyEiIiIzl52dDR8fH2RmZsLFxaXEduW6dHv69Gl8++232LhxI/Lz8/H6669j5MiRaNmyZYWKlqvs7Gy4uLggKyuLQY+IiIgMTtvsofPyKgAQFBSEJUuW4Pr16/jmm29w/fp1hIWFoUmTJvjyyy+RlZVV7sKJiIiISD/KFfQKKZVK5OfnIy8vD6Iows3NDStXroSPjw82bdqkrxqJiIiIqBzKFfTi4uIwYcIEeHl5YcqUKWjevDkSExMRExOD8+fPY9asWZg0aZK+ayUiIiIiHeg8Ri8oKAiJiYmIiIjAqFGj0LNnT1haWmq0uXXrFjw8PKBUKvVarKniGD0iIiIyJm2zh87Lq/Tr1w8jR45EzZo1S2xTo0aNShPyiIiIiEwVF0zWA/boERERkTEZpEfv2rVrWLlyJQ4dOoT09HQIggAPDw+0adMGY8eO5aLBRERERCZE6x69AwcOoFu3bvDx8UFERAQ8PDwgiiIyMjIQHR2N1NRUbN++HWFhYYau2eSwR4+IiIiMSdvsoXXQa9myJdq2bYslS5YU+/yUKVNw4MABHDt2rHwVyxiDHhERERmT3hdMPnPmDMaOHVvi82PGjMGZM2d0q5KIiIiIDEbroOfl5YVDhw6V+HxsbCy8vLz0UhQRERERVZzWkzHeffddjB07FnFxcejcuTM8PDwgCALS09MRHR2Nr7/+GkuXLjVgqURERESkC62D3rhx41CtWjUsWbIEq1evhkKhAABYWloiODgY33//Pfr372+wQomIiIhIN+VaR6+goAC3b98GAFSvXh3W1tZ6L0xOOBmDiIiIjMlgO2MAgLW1NcfjEREREZk4rSdjlOXSpUv4z3/+o6/TEREREVEF6S3o3b9/HzExMfo6HRERERFVkNaXbpctW1bq89evX69wMURERESkP1oHvcjISHh5ecHGxqbY5/Pz8/VWFBERERFVnNZBz9fXFwsXLixxCZX4+HgEBwfrrTAiIiIiqhitx+gFBwcjLi6uxOcFQUA5VmohIiIiIgPRukdvzpw5yM3NLfH5gIAAJCcn66UoIiIiIqo4rXv0AgICEBISUuLz1tbW8PX1VT8+ePAg8vLyKlYdEREREZWb3pZXeV63bt04E5eIiIhIQgYLehyvR0RERCQtgwU9IiIiIpKW7ILeihUr4OfnBzs7OwQHB2P//v2lto+JiUFwcDDs7OxQt25drFq1qsS2P//8MwRBQJ8+ffRcNREREZHxySrobdq0CZGRkZg5cyZOnjyJ8PBwdOvWDSkpKcW2T05ORvfu3REeHo6TJ0/i/fffx6RJk7B58+Yiba9evYp3330X4eHhhv4YREREREYhiAYaTOfs7Iz4+HjUrVtXb+d88cUX0aJFC6xcuVJ9rFGjRujTpw8WLFhQpP20adMQFRWFxMRE9bGxY8fi1KlTiI2NVR9TKBRo3749RowYgf379yMzMxO//fab1nVlZ2fDxcUFWVlZcHZ2Lt+HIyIiItKSttlDNpMx8vPzERcXh4iICI3jEREROHToULGviY2NLdK+S5cuOH78OAoKCtTH5syZgxo1auDNN9/Uqpa8vDxkZ2dr3IiIiIhMjdYLJj9LoVDg9u3bEAQB1apVg6WlZZE2OTk5FS7uWbdv34ZCoYCHh4fGcQ8PD6Snpxf7mvT09GLbP378GLdv34aXlxcOHjyIdevWIT4+XutaFixYgI8//ljnz0BERERkTDr16G3duhVhYWFwcHCAt7c3vLy84ODggLCwMJ0udVaEIAgaj0VRLHKsrPaFx3NycvDGG29g7dq1qF69utY1zJgxA1lZWepbamqqDp+AiIiIyDi07tFbvXo1Jk2ahJEjR+K///0vPDw8IIoiMjIy8Pfff+P111/HV199hVGjRhmk0OrVq8PS0rJI711GRkaRXrtCnp6exba3srJCtWrVcPbsWVy5cgU9e/ZUP69UKgEAVlZWSEpKwgsvvFDkvLa2trC1ta3oRyIiIiIyKK2D3meffYYVK1YUO46tT58+aNmyJT755BODBT0bGxsEBwcjOjoar7zyivp4dHQ0evfuXexrQkND8ccff2gc27lzJ0JCQmBtbQ1/f38kJCRoPP/BBx8gJycHX375JXx8fPT/QYiIiIiMROugd/36dbRt27bE59u0aYMbN27opaiSTJ06FUOGDEFISAhCQ0OxZs0apKSkYOzYsQBUl1SvX7+O77//HoBqhu3y5csxdepUjBo1CrGxsVi3bh1++uknAICdnR0CAwM13qNq1aoAUOQ4ERERkdxoHfQaN26MNWvW4Isvvij2+bVr16Jx48Z6K6w4AwYMwJ07dzBnzhykpaUhMDAQ27Ztg6+vLwAgLS1NY009Pz8/bNu2DVOmTMH//d//wdvbG8uWLUPfvn0NWicRERGRKdB6Hb2YmBj06NEDvr6+iIiIgIeHBwRBQHp6OqKjo3H16lVs27atUi44zHX0iIiIyJi0zR5a9+i1b98eZ86cwcqVK3H48GH1JAdPT0+8/PLLGDt2LOrUqVPhwomIiIhIPwy2M0Zlwh49IiIiMibJd8YgIiIiImnpLegNGzYM//nPf/R1OiIiIiKqoHJtgVacmjVrwsKCHYREREREpoJj9PSAY/SIiIjImDhGj4iIiKiS0+nS7bVr17By5UocOnQI6enpEAQBHh4eaNOmDcaOHcstw4iIiIhMiNaXbg8cOIBu3brBx8dHvWCyKIrIyMhAdHQ0UlNTsX37doSFhRm6ZpPDS7dERERkTNpmD62DXsuWLdG2bVssWbKk2OenTJmCAwcO4NixY+WrWMYY9IiIiMiY9D5G78yZMxg7dmyJz48ZMwZnzpzRrUoiIiIiMhitg56XlxcOHTpU4vOxsbHw8vLSS1FEREREVHFaT8Z49913MXbsWMTFxaFz587w8PCAIAhIT09HdHQ0vv76ayxdutSApRIRERGRLrQOeuPGjUO1atWwZMkSrF69GgqFAgBgaWmJ4OBgfP/99+jfv7/BCiUiIiIi3ZRrweSCggLcvn0bAFC9enVYW1vrvTA54WQMIiIiMiZts0e5tkCztrbmeDwiIiIiE6e3nTFWrFiBOXPm6Ot0RERERFRBegt6mzdvxvr16/V1OiIiIiKqoHJdui3Orl279HUqIiIiItIDvfXoEREREZFpqVCPXnx8PC5cuAAvLy+EhYVBEAR91UVEREREFaR1j96gQYOQk5MDALh//z66dOmCFi1a4I033kC7du3QqlUrZGZmGqpOIiIiItKR1kFv06ZNePjwIQDg448/xoULF3D8+HHk5eXh9OnTePDgAWfdEhEREZkQrYPes+sqb9++HZ9++ilatGgBAAgMDMTnn3+OP//8U/8VEhEREVG56DQZo3AM3s2bNxEYGKjxXOPGjZGamqq/yoiIiIioQnSajPHhhx/CwcEBFhYWSE9PR0BAgPq527dvw9HRUe8FEhEREVH5aB302rVrh6SkJABAQEAAkpOTNZ7ftm0bGjdurN/qiIiIiKjcBPHZwXcVcPnyZdjY2KBWrVr6OJ2saLuxMBEREZE+aJs9dF4w+cqVK8Uer1u3bqUMeURERESmSuegV7duXbRt2xarV6/G3bt3DVETEREREemBzkHv+PHjCA0Nxbx58+Dt7Y3evXvjf//7H/Ly8gxRHxERERGVk85Br0WLFvjss8+QkpKC7du3w93dHWPGjIG7uztGjhxpiBqJiIiIqBz0MhnjxIkTePPNN3H69GkoFAp91CUrnIxBRERExmSwyRiFUlNTsWjRIjRr1gwtW7ZElSpVsHz58vKeTmsrVqyAn58f7OzsEBwcjP3795faPiYmBsHBwbCzs0PdunWxatUqjefXrl2L8PBwuLq6wtXVFZ06dcLRo0cN+RGIiIiIjELnoLdmzRq0b98efn5++O6779C/f39cunQJBw4cwNtvv22IGtU2bdqEyMhIzJw5EydPnkR4eDi6deuGlJSUYtsnJyeje/fuCA8Px8mTJ/H+++9j0qRJ2Lx5s7rN3r17MXDgQOzZswexsbGoXbs2IiIicP36dYN+FiIiIiJD0/nSrY+PD15//XUMHjwYzZo1M1BZxXvxxRfRokULrFy5Un2sUaNG6NOnDxYsWFCk/bRp0xAVFYXExET1sbFjx+LUqVOIjY0t9j0UCgVcXV2xfPlyDB06VKu6eOmWiIiIjMlgl25TUlLw2WeflRnyxo0bh9u3b+t6+hLl5+cjLi4OERERGscjIiJw6NChYl8TGxtbpH2XLl1w/PhxFBQUFPua3NxcFBQUwM3NTT+FExEREUlE56AnCIJW7TZs2IDs7GydCyrJ7du3oVAo4OHhoXHcw8MD6enpxb4mPT292PaPHz8uMYROnz4dNWvWRKdOnUqsJS8vD9nZ2Ro3IiIiIlNT7skYZdHTzmpFPB80RVEsNXwW17644wCwaNEi/PTTT9iyZQvs7OxKPOeCBQvg4uKivvn4+OjyEYiIiIiMwmBBT9+qV68OS0vLIr13GRkZRXrtCnl6ehbb3srKCtWqVdM4/vnnn2P+/PnYuXMngoKCSq1lxowZyMrKUt9SU1PL8YmIiIiIDEs2Qc/GxgbBwcGIjo7WOB4dHY02bdoU+5rQ0NAi7Xfu3ImQkBBYW1urj3322WeYO3cuduzYgZCQkDJrsbW1hbOzs8aNiIiIyNTIJugBwNSpU/H111/jm2++QWJiIqZMmYKUlBSMHTsWgKqn7dmZsmPHjsXVq1cxdepUJCYm4ptvvsG6devw7rvvqtssWrQIH3zwAb755hvUqVMH6enpSE9Px/37943++YiIiIj0yUrqAnQxYMAA3LlzB3PmzEFaWhoCAwOxbds2+Pr6AgDS0tI01tTz8/PDtm3bMGXKFPzf//0fvL29sWzZMvTt21fdZsWKFcjPz8drr72m8V6zZs3C7NmzjfK5iIiIiAxBL1ugFeftt9/G3LlzUb16dUOc3qRwHT0iIiIyJoNugbZ//3688cYbCA0NVe8g8cMPP+DAgQPqNitXrqwUIY+IiIjIVOkc9DZv3owuXbrA3t4eJ0+eRF5eHgAgJycH8+fP13uBRERERFQ+Oge9efPmYdWqVVi7dq3GzNU2bdrgxIkTei2OiIiIiMpP56CXlJSEdu3aFTnu7OyMzMxMfdRERERERHqgc9Dz8vLCxYsXixw/cOAA6tatq5eiiIiIiKjidA56Y8aMweTJk3HkyBEIgoAbN27gxx9/xLvvvotx48YZokYiIiIiKged19F77733kJWVhZdeegmPHj1Cu3btYGtri3fffRcTJkwwRI1EREREVA7lXkcvNzcX586dg1KpREBAABwdHfVdm2xwHT0iIiIyJm2zR7l3xnBwcNBqX1giIiIikobOQe/Bgwf49NNPsWvXLmRkZECpVGo8f/nyZb0VR0RERETlp3PQe+uttxATE4MhQ4bAy8sLgiAYoi4iIiIiqiCdg9727dvx119/ISwszBD1EBEREZGe6Ly8iqurK9zc3AxRCxERERHpkc5Bb+7cufjoo4+Qm5triHqIiIiISE90vnT7xRdf4NKlS/Dw8ECdOnU09rsFwP1uiYiIiEyEzkGvT58+BiiDiIiIiPSt3Asm01NcMJmIiIiMSdvsofMYPSIiIiKSB60u3bq5ueHff/9F9erV4erqWuraeXfv3tVbcURERERUfloFvSVLlsDJyQkAsHTpUkPWQ0RERER6wjF6esAxekRERGRM2mYPrXr0srOztX5jBh0iIiIi06BV0KtatWqZe9qKoghBEKBQKPRSGBERERFVjFZBb8+ePYaug4iIiIj0TKug1759e0PXQURERER6Vq519Pbv34833ngDbdq0wfXr1wEAP/zwAw4cOKDX4oiIiIio/HQOeps3b0aXLl1gb2+PEydOIC8vDwCQk5OD+fPn671AIiIiIiofnYPevHnzsGrVKqxduxbW1tbq423atMGJEyf0WhwRERERlZ/OQS8pKQnt2rUrctzZ2RmZmZn6qImIiIiI9EDnoOfl5YWLFy8WOX7gwAHUrVtXL0URERERUcXpHPTGjBmDyZMn48iRIxAEATdu3MCPP/6Id999F+PGjTNEjURERERUDlotr/Ks9957D1lZWXjppZfw6NEjtGvXDra2tnj33XcxYcIEQ9RIREREROVQ7r1uc3Nzce7cOSiVSgQEBMDR0VHftckG97olIiIiY9I2e5RrHT0AcHBwQEhICPz9/fHPP/8gMTGxvKciIiIiIgPQOej1798fy5cvBwA8fPgQLVu2RP/+/REUFITNmzfrvcDnrVixAn5+frCzs0NwcDD2799favuYmBgEBwfDzs4OdevWxapVq4q02bx5MwICAmBra4uAgABs3brVUOUTERERGY3OQW/fvn0IDw8HAGzduhVKpRKZmZlYtmwZ5s2bp/cCn7Vp0yZERkZi5syZOHnyJMLDw9GtWzekpKQU2z45ORndu3dHeHg4Tp48iffffx+TJk3SCKSxsbEYMGAAhgwZglOnTmHIkCHo378/jhw5YtDPQkRERGRoOo/Rs7e3x7///gsfHx8MHToU3t7e+PTTT5GSkoKAgADcv3/fULXixRdfRIsWLbBy5Ur1sUaNGqFPnz5YsGBBkfbTpk1DVFSUxmXlsWPH4tSpU4iNjQUADBgwANnZ2di+fbu6TdeuXeHq6oqffvpJq7qMNkbv8YOSnxMsAUs77drCArCyL2fbXAAl/cgIgJVDOds+BKAsuQyrKuVrq3gEiAr9tLV0AAThSds8QHysp7b2gPDk/1yKfEAs0E9bCzvAwlL3tsoCQJlfSltbwMKqHG0fA8q8UtraABbW5WirAJSPSm4rWAOWNrq3FZWA4qGe2loBlrZP2oqAIldPbXX4e8/fEcW35e8I3dvyd4Tqfll/75/92TEAbbOHzrNufXx8EBsbCzc3N+zYsQM///wzAODevXuws7Mr49Xll5+fj7i4OEyfPl3jeEREBA4dOlTsa2JjYxEREaFxrEuXLli3bh0KCgpgbW2N2NhYTJkypUibpUuXllhLXl6eeus3QPWHbRS/lDLhxbs70OGvp483u5f8D4R7e6DT3qePf68D5N0uvq1bCND12NPHfwUAD64W39YlAOhx9unjv1sCWeeKb1vFF+h95enjf9oBd48X39a2OtD31tPHe7sBGTHFt7V0AAY884/S/r7AjW3FtwWAQc/8I3NoCJD6a8lt+99/+hf36Bgg+buS276aAdjVUN0/MRW4sKLktr2SAcc6qvunZwKJn5fctvsZoGpj1f2z84EzH5fctstRoFpL1f2kL4H490pu23EP4NFBdf/iGuB4KTPo2/8J1Oyhun/lR+DwiJLbtv0FqN1Pdf/aVuBA/5Lbtv4WqDtcdT/tbyDm5ZLbhiwHGoxX3b+1H9j1Usltmy0CAv6run/vBPB3q5LbBs4Cgmar7mclAtsCS27b6F2g+Weq+w9SgCi/ktvWHwe0/D/V/bzbwBb3ktv6DQNC16vuK3JL/3vv8xoQ/r+nj/k7QoW/I1T3+TtCdV+q3xGDyjXXVe90vnQbGRmJwYMHo1atWvD29kaHDh0AqC7pNmnSRN/1qd2+fRsKhQIeHh4axz08PJCenl7sa9LT04tt//jxY9y+fbvUNiWdEwAWLFgAFxcX9c3Hx6c8H4mIiIjIoMq1vEpcXBxSUlLQuXNn9bIqf/31F6pWrYqwsDC9FwkAN27cQM2aNXHo0CGEhoaqj3/yySf44YcfcP78+SKvadCgAUaMGIEZM2aojx08eBBt27ZFWloaPD09YWNjg++++w4DBw5Ut/nxxx/x5ptv4tGj4rtvi+vR8/Hx4aVbXpYpZ1tellG1ldllGV66LaEtf0eo2vJ3hO5tzex3hFwv3QJAcHAwgoODNY716NGjPKfSWvXq1WFpaVmkpy0jI6NIj1whT0/PYttbWVmhWrVqpbYp6ZwAYGtrC1tb2/J8jIrR5YfGYG0dym5Trrb2ZbcpT1tLHYYT6NTWFoCWPwM6tbUBYCNtWwvrp78g9drW6ukvdL22tQQstPwZ1qWtYKH93w2d2gqGaQuYSFv+jlC15e8I3dua8e8ICZUr6F27dg1RUVFISUlBfr5mUl+8eLFeCnuejY0NgoODER0djVdeeUV9PDo6Gr179y72NaGhofjjjz80ju3cuRMhISGwtrZWt4mOjtYYp7dz5060adPGAJ+CiIiIyHh0Dnq7du1Cr1694Ofnh6SkJAQGBuLKlSsQRREtWrQwRI1qU6dOxZAhQxASEoLQ0FCsWbMGKSkpGDt2LABgxowZuH79Or7//nsAqhm2y5cvx9SpUzFq1CjExsZi3bp1GrNpJ0+ejHbt2mHhwoXo3bs3fv/9d/zzzz84cOCAQT8LERERkaHpPBljxowZeOedd3DmzBnY2dlh8+bNSE1NRfv27dGvXz9D1Kg2YMAALF26FHPmzEGzZs2wb98+bNu2Db6+vgCAtLQ0jTX1/Pz8sG3bNuzduxfNmjXD3LlzsWzZMvTt21fdpk2bNvj555/x7bffIigoCOvXr8emTZvw4osvGvSzEBERERmazpMxnJycEB8fjxdeeAGurq44cOAAGjdujFOnTqF37964cuWKgUo1XdzrloiIiIzJYHvdVqlSRT3j1NvbG5cuXVI/V7hkCRERERFJT+cxeq1bt8bBgwcREBCAHj164J133kFCQgK2bNmC1q1bG6JGIiIiIioHnYPe4sWL1duczZ49G/fv38emTZtQr149LFmyRO8FEhEREVH5lGvBZNLEMXpERERkTAZdMBkAjh8/jsTERAiCgEaNGhVZQJmIiIiIpKVz0Lt27RoGDhyIgwcPomrVqgCAzMxMtGnTBj/99BP3fSUiIiIyETrPuh05ciQKCgqQmJiIu3fv4u7du0hMTIQoinjzzTcNUSMRERERlYPOY/Ts7e1x6NAhNG/eXOP4iRMnEBYWhocPS9ng10xxjB4REREZk8HW0atduzYKCgqKHH/8+DFq1qyp6+mIiIiIyEB0DnqLFi3CxIkTcfz4cRR2Bh4/fhyTJ0/G559/rvcCiYiIiKh8tLp06+rqCkEQ1I8fPHiAx48fw8pKNZej8H6VKlVw9+5dw1VronjploiIiIxJr8urLF26VF91EREREZGRaBX0hg0bZug6iIiIiEjPdB6jR0RERETywKBHREREZKYY9IiIiIjMFIMeERERkZmqcNDLzs7Gb7/9hsTERH3UQ0RERER6onPQ69+/P5YvXw4AePjwIUJCQtC/f38EBQVh8+bNei+QiIiIiMpH56C3b98+hIeHAwC2bt0KURSRmZmJZcuWYd68eXovkIiIiIjKR+egl5WVBTc3NwDAjh070LdvXzg4OKBHjx64cOGC3gskIiIiovLROej5+PggNjYWDx48wI4dOxAREQEAuHfvHuzs7PReIBERERGVj1Y7YzwrMjISgwcPhqOjI3x9fdGhQwcAqku6TZo00Xd9RERERFROOge9cePGoVWrVkhNTUXnzp1hYaHqFKxbty7H6BERERGZEEEURVHqIuQuOzsbLi4uyMrKgrOzs9TlEBERkZnTNnvo3KM3cuTIUp//5ptvdD0lERERERmAzkHv3r17Go8LCgpw5swZZGZm4j//+Y/eCiMiIiKiitE56G3durXIMaVSiXHjxqFu3bp6KYqIiIiIKk4ve91aWFhgypQpWLJkiT5OR0RERER6oJegBwCXLl3C48eP9XU6IiIiIqognS/dTp06VeOxKIpIS0vDX3/9hWHDhumtMCIiIiKqGJ2D3smTJzUeW1hYoEaNGvjiiy/KnJFLRERERMajc9Dbs2ePIeogIiIiIj3T2xg9IiIiIjItOge9mzdvYsiQIfD29oaVlRUsLS01boZy7949DBkyBC4uLnBxccGQIUOQmZlZ6mtEUcTs2bPh7e0Ne3t7dOjQAWfPnlU/f/fuXUycOBENGzaEg4MDateujUmTJiErK8tgn4OIiIjIWHS+dDt8+HCkpKTgww8/hJeXFwRBMERdRQwaNAjXrl3Djh07AACjR4/GkCFD8Mcff5T4mkWLFmHx4sVYv349GjRogHnz5qFz585ISkqCk5MTbty4gRs3buDzzz9HQEAArl69irFjx+LGjRv49ddfjfK5iIiIiAxF571unZycsH//fjRr1sxAJRWVmJiIgIAAHD58GC+++CIA4PDhwwgNDcX58+fRsGHDIq8RRRHe3t6IjIzEtGnTAAB5eXnw8PDAwoULMWbMmGLf63//+x/eeOMNPHjwAFZW2uVg7nVLRERExqRt9tD50q2Pjw90zIYVFhsbCxcXF3XIA4DWrVvDxcUFhw4dKvY1ycnJSE9PR0REhPqYra0t2rdvX+JrAKj/wEoLeXl5ecjOzta4EREREZkanYPe0qVLMX36dFy5csUA5RQvPT0d7u7uRY67u7sjPT29xNcAgIeHh8ZxDw+PEl9z584dzJ07t8TevkILFixQjxV0cXGBj4+PNh+DiIiIyKh0DnoDBgzA3r178cILL8DJyQlubm4aN13Mnj0bgiCUejt+/DgAFDsWUBTFMscIPv98Sa/Jzs5Gjx49EBAQgFmzZpV6zhkzZiArK0t9S01NLeujEhERERmdzpMxli5dqrc3nzBhAl5//fVS29SpUwenT5/GzZs3izx369atIj12hTw9PQGoeva8vLzUxzMyMoq8JicnB127doWjoyO2bt0Ka2vrUmuytbWFra1tqW2IiIiIpKZz0NPnNmfVq1dH9erVy2wXGhqKrKwsHD16FK1atQIAHDlyBFlZWWjTpk2xr/Hz84Onpyeio6PRvHlzAEB+fj5iYmKwcOFCdbvs7Gx06dIFtra2iIqKgp2dnR4+GREREZH0tAp62dnZ6hkdZU08MMSs00aNGqFr164YNWoUVq9eDUC1vMrLL7+sMePW398fCxYswCuvvAJBEBAZGYn58+ejfv36qF+/PubPnw8HBwcMGjQIgKonLyIiArm5udiwYYPGxIoaNWoYdF1AIiIiIkPTKui5uroiLS0N7u7uqFq1aqnj5RQKhd6LBIAff/wRkyZNUs+i7dWrF5YvX67RJikpSWOx4/feew8PHz7EuHHjcO/ePbz44ovYuXMnnJycAABxcXE4cuQIAKBevXoa50pOTkadOnUM8lmIiIiIjEGrdfRiYmIQFhYGKysrxMTElNq2ffv2eitOLriOHhERERmTttlD5wWTqSgGPSIiIjImbbOHzpMxAODRo0c4ffo0MjIyoFQqNZ7r1atXeU5JRERERHqmc9DbsWMHhg4ditu3bxd5zpBj9IiIiIhINzovmDxhwgT069cPaWlpUCqVGjeGPCIiIiLToXPQy8jIwNSpU0tcqJiIiIiITIPOQe+1117D3r17DVAKEREREemTzrNuc3Nz0a9fP9SoUQNNmjQpsl3YpEmT9FqgHHDWLRERERmTwWbdbty4EX///Tfs7e2xd+9ejcWTBUGolEGPiIiIyBTpHPQ++OADzJkzB9OnT4eFhc5XfomIiIjISHROavn5+RgwYABDHhEREZGJ0zmtDRs2DJs2bTJELURERESkRzpfulUoFFi0aBH+/vtvBAUFFZmMsXjxYr0VR0RERETlp3PQS0hIQPPmzQEAZ86c0Xju2YkZRERERCQtnYPenj17DFEHEREREekZZ1QQERERmSkGPSIiIiIzxaBHREREZKYY9IiIiIjMFIMeERERkZnSatZtVFSU1ifs1atXuYshIiIiIv3RKuj16dNHq5MJggCFQlGReoiIiIhIT7QKekql0tB1EBEREZGeVWiM3qNHj/RVBxERERHpmc5BT6FQYO7cuahZsyYcHR1x+fJlAMCHH36IdevW6b1AIiIiIiofnYPeJ598gvXr12PRokWwsbFRH2/SpAm+/vprvRZHREREROWnc9D7/vvvsWbNGgwePBiWlpbq40FBQTh//rxeiyMiIiKi8tM56F2/fh316tUrclypVKKgoEAvRRERERFRxekc9Bo3boz9+/cXOf6///0PzZs310tRRERERFRxWi2v8qxZs2ZhyJAhuH79OpRKJbZs2YKkpCR8//33+PPPPw1RIxERERGVg849ej179sSmTZuwbds2CIKAjz76CImJifjjjz/QuXNnQ9RIREREROUgiKIoSl2E3GVnZ8PFxQVZWVlwdnaWuhwiIiIyc9pmD50v3RbKz89HRkZGkV0zateuXd5TEhEREZEe6Rz0Lly4gJEjR+LQoUMax0VR5F63RERERCZE5zF6w4cPh4WFBf7880/ExcXhxIkTOHHiBE6ePIkTJ04YokYAwL179zBkyBC4uLjAxcUFQ4YMQWZmZqmvEUURs2fPhre3N+zt7dGhQwecPXu2xLbdunWDIAj47bff9P8BiIiIiIxM5x69+Ph4xMXFwd/f3xD1lGjQoEG4du0aduzYAQAYPXo0hgwZgj/++KPE1yxatAiLFy/G+vXr0aBBA8ybNw+dO3dGUlISnJycNNouXboUgiAY9DMQERERGZPOQS8gIAC3b982RC0lSkxMxI4dO3D48GG8+OKLAIC1a9ciNDQUSUlJaNiwYZHXiKKIpUuXYubMmXj11VcBAN999x08PDywceNGjBkzRt321KlTWLx4MY4dOwYvLy/jfCgiIiIiA9P50u3ChQvx3nvvYe/evbhz5w6ys7M1boYQGxsLFxcXdcgDgNatW8PFxaXIWMFCycnJSE9PR0REhPqYra0t2rdvr/Ga3NxcDBw4EMuXL4enp6dW9eTl5RnlcxMRERFVhM49ep06dQIAdOzYUeO4ISdjpKenw93dvchxd3d3pKenl/gaAPDw8NA47uHhgatXr6ofT5kyBW3atEHv3r21rmfBggX4+OOPtW5PREREJAWdg96ePXv09uazZ88uMzAdO3YMAIodP1cYLkvz/PPPviYqKgq7d+/GyZMndSkbM2bMwNSpU9WPs7Oz4ePjo9M5iIiIiAxN56DXvn17vb35hAkT8Prrr5fapk6dOjh9+jRu3rxZ5Llbt24V6bErVHgZNj09XWPcXUZGhvo1u3fvxqVLl1C1alWN1/bt2xfh4eHYu3dvsee2tbWFra1tqXUTERERSa3cCybn5uYiJSUF+fn5GseDgoK0Pkf16tVRvXr1MtuFhoYiKysLR48eRatWrQAAR44cQVZWFtq0aVPsa/z8/ODp6Yno6Gg0b94cgGqR55iYGCxcuBAAMH36dLz11lsar2vSpAmWLFmCnj17av05iIiIiEyRzkHv1q1bGDFiBLZv317s84YYo9eoUSN07doVo0aNwurVqwGolld5+eWXNWbc+vv7Y8GCBXjllVcgCAIiIyMxf/581K9fH/Xr18f8+fPh4OCAQYMGAVD1+hU3AaN27drw8/PT++cgIiIiMiadZ91GRkbi3r17OHz4MOzt7bFjxw589913qF+/PqKiogxRIwDgxx9/RJMmTRAREYGIiAgEBQXhhx9+0GiTlJSErKws9eP33nsPkZGRGDduHEJCQnD9+nXs3LmzyBp6REREROZIEEVR1OUFXl5e+P3339GqVSs4Ozvj+PHjaNCgAaKiorBo0SIcOHDAULWaLG03FiYiIiLSB22zh849eg8ePFAvdeLm5oZbt24BUI1tM+QWaERERESkG52DXsOGDZGUlAQAaNasGVavXo3r169j1apV3FWCiIiIyIToPBkjMjISaWlpAIBZs2ahS5cu+PHHH2FjY4P169fruz4iIiIiKiedx+g9Lzc3F+fPn0ft2rW1WirFHHGMHhERERmTttmj3OvoAapdJuzt7dGiRYuKnIaIiIiIDEDnMXoAsG7dOgQGBsLOzg52dnYIDAzE119/re/aiIiIiKgCdO7R+/DDD7FkyRJMnDgRoaGhAIDY2FhMmTIFV65cwbx58/ReJBERERHpTucxetWrV8dXX32FgQMHahz/6aefMHHiRNy+fVuvBcoBx+gRERGRMRlsHT2FQoGQkJAix4ODg/H48WNdT0dEREREBqJz0HvjjTewcuXKIsfXrFmDwYMH66UoIiIiIqq4cs26XbduHXbu3InWrVsDAA4fPozU1FQMHToUU6dOVbdbvHixfqokIiIiIp3pHPTOnDmjXk7l0qVLAIAaNWqgRo0aOHPmjLqdIAh6KpGIiIiIykPnoLdnzx5D1EFEREREelaudfSIiIiIyPQx6BERERGZKQY9IiIiIjPFoEdERERkphj0iIiIiMxUudbRS0pKwldffYXExEQIggB/f39MnDgRDRs21Hd9RERERFROOvfo/frrrwgMDERcXByaNm2KoKAgnDhxAoGBgfjf//5niBqJiIiIqBwEURRFXV5Qt25dvPHGG5gzZ47G8VmzZuGHH37A5cuX9VqgHGi7sTARERGRPmibPXTu0UtPT8fQoUOLHH/jjTeQnp6u6+mIiIiIyEB0DnodOnTA/v37ixw/cOAAwsPD9VIUEREREVWczpMxevXqhWnTpiEuLg6tW7cGABw+fBj/+9//8PHHHyMqKkqjLRERERFJQ+cxehYW2nUCCoIAhUJRrqLkhmP0iIiIyJi0zR469+gplcoKFUZERERExlGudfRIU2GnaHZ2tsSVEBERUWVQmDnKujBbrqAXExODzz//XL1gcqNGjfDf//630k7GyMnJAQD4+PhIXAkRERFVJjk5OXBxcSnxeZ3H6G3YsAEjRozAq6++irCwMIiiiEOHDmHr1q1Yv349Bg0aVOGi5UapVOLGjRtwcnKCIAhSlyM72dnZ8PHxQWpqKsc4ygS/Z/LE75v88HsmT8b4vomiiJycHHh7e5c6f0LnoNeoUSOMHj0aU6ZM0Ti+ePFirF27FomJieWrmCotTmaRH37P5InfN/nh90yeTOn7pvM6epcvX0bPnj2LHO/VqxeSk5P1UhQRERERVZzOQc/Hxwe7du0qcnzXrl0co0ZERERkQrSejDFy5Eh8+eWXeOeddzBp0iTEx8ejTZs2EAQBBw4cwPr16/Hll18aslYyU7a2tpg1axZsbW2lLoW0xO+ZPPH7Jj/8nsmTKX3ftB6jZ2lpibS0NLi7u2Pr1q344osv1OPxCmfd9u7d26DFEhEREZH2tA56FhYWSE9Ph7u7u6FrIiIiIiI90GmMHpcOISIiIpIPnXr0XFxcygx7d+/e1UthRERERFQxOu2M8fHHH5e6+jIRERERmRBRS4IgiDdv3tS2OVGp5s+fL4aEhIiOjo5ijRo1xN69e4vnz5+XuizSwfz580UA4uTJk6Uuhcpw7do1cfDgwaKbm5tob28vNm3aVDx+/LjUZVEpCgoKxJkzZ4p16tQR7ezsRD8/P/Hjjz8WFQqF1KXRM2JiYsSXX35Z9PLyEgGIW7du1XheqVSKs2bNEr28vEQ7Ozuxffv24pkzZ4xao9Zj9Dg+j/QpJiYG48ePx+HDhxEdHY3Hjx8jIiICDx48kLo00sKxY8ewZs0aBAUFSV0KleHevXsICwuDtbU1tm/fjnPnzuGLL75A1apVpS6NSrFw4UKsWrUKy5cvR2JiIhYtWoTPPvsMX331ldSl0TMePHiApk2bYvny5cU+v2jRIixevBjLly/HsWPH4Onpic6dOyMnJ8doNXLWLZmEW7duwd3dHTExMWjXrp3U5VAp7t+/jxYtWmDFihWYN28emjVrhqVLl0pdFpVg+vTpOHjwIPbv3y91KaSDl19+GR4eHli3bp36WN++feHg4IAffvhBwsqoJIIgYOvWrejTpw8A1V603t7eiIyMxLRp0wAAeXl58PDwwMKFCzFmzBij1KV1j55SqWTII4PJysoCALi5uUlcCZVl/Pjx6NGjBzp16iR1KaSFqKgohISEoF+/fnB3d0fz5s2xdu1aqcuiMrRt2xa7du3Cv//+CwA4deoUDhw4gO7du0tcGWkrOTkZ6enpiIiIUB+ztbVF+/btcejQIaPVodNkDCJDEEURU6dORdu2bREYGCh1OVSKn3/+GSdOnMCxY8ekLoW0dPnyZaxcuRJTp07F+++/j6NHj2LSpEmwtbXF0KFDpS6PSjBt2jRkZWXB398flpaWUCgU+OSTTzBw4ECpSyMtpaenAwA8PDw0jnt4eODq1atGq4NBjyQ3YcIEnD59GgcOHJC6FCpFamoqJk+ejJ07d8LOzk7qckhLSqUSISEhmD9/PgCgefPmOHv2LFauXMmgZ8I2bdqEDRs2YOPGjWjcuDHi4+MRGRkJb29vDBs2TOrySAfPz3EQRdGo8x4Y9EhSEydORFRUFPbt24datWpJXQ6VIi4uDhkZGQgODlYfUygU2LdvH5YvX468vDxYWlpKWCEVx8vLCwEBARrHGjVqhM2bN0tUEWnjv//9L6ZPn47XX38dANCkSRNcvXoVCxYsYNCTCU9PTwCqnj0vLy/18YyMjCK9fIak084YRPoiiiImTJiALVu2YPfu3fDz85O6JCpDx44dkZCQgPj4ePUtJCQEgwcPRnx8PEOeiQoLC0NSUpLGsX///Re+vr4SVUTayM3NhYWF5j/RlpaWUCqVElVEuvLz84Onpyeio6PVx/Lz8xETE4M2bdoYrQ726JEkxo8fj40bN+L333+Hk5OTeiyDi4sL7O3tJa6OiuPk5FRkDGWVKlVQrVo1jq00YVOmTEGbNm0wf/589O/fH0ePHsWaNWuwZs0aqUujUvTs2ROffPIJateujcaNG+PkyZNYvHgxRo4cKXVp9Iz79+/j4sWL6sfJycmIj4+Hm5sbateujcjISMyfPx/169dH/fr1MX/+fDg4OGDQoEHGK9Koq/YRPQGg2Nu3334rdWmkg/bt23PBZBn4448/xMDAQNHW1lb09/cX16xZI3VJVIbs7Gxx8uTJYu3atUU7Ozuxbt264syZM8W8vDypS6Nn7Nmzp9h/y4YNGyaK4tMFkz09PUVbW1uxXbt2YkJCglFr1HodPSIiIiKSF47RIyIiIjJTDHpEREREZopBj4iIiMhMMegRERERmSkGPSIiIiIzxaBHREREZKYY9IiIiIjMFIMeERERkZli0CMiKgdRFNGpUyd06dKlyHMrVqyAi4sLUlJSJKiMiOgpBj0ionIQBAHffvstjhw5gtWrV6uPJycnY9q0afjyyy9Ru3Ztvb5nQUGBXs9HROaPQY+IqJx8fHzw5Zdf4t1330VycjJEUcSbb76Jjh07olWrVujevTscHR3h4eGBIUOG4Pbt2+rX7tixA23btkXVqlVRrVo1vPzyy7h06ZL6+StXrkAQBPzyyy/o0KED7OzssGHDBly9ehU9e/aEq6srqlSpgsaNG2Pbtm1SfHwikgHudUtEVEF9+vRBZmYm+vbti7lz5+LYsWMICQnBqFGjMHToUDx8+BDTpk3D48ePsXv3bgDA5s2bIQgCmjRpggcPHuCjjz7ClStXEB8fDwsLC1y5cgV+fn6oU6cOvvjiCzRv3hy2trYYPXo08vPz8cUXX6BKlSo4d+4cnJ2d0a5dO4n/FIjIFDHoERFVUEZGBgIDA3Hnzh38+uuvOHnyJI4cOYK///5b3ebatWvw8fFBUlISGjRoUOQct27dgru7OxISEhAYGKgOekuXLsXkyZPV7YKCgtC3b1/MmjXLKJ+NiOSNl26JiCrI3d0do0ePRqNGjfDKK68gLi4Oe/bsgaOjo/rm7+8PAOrLs5cuXcKgQYNQt25dODs7w8/PDwCKTOAICQnReDxp0iTMmzcPYWFhmDVrFk6fPm2ET0hEcsWgR0SkB1ZWVrCysgIAKJVK9OzZE/Hx8Rq3CxcuqC+x9uzZE3fu3MHatWtx5MgRHDlyBACQn5+vcd4qVapoPH7rrbdw+fJlDBkyBAkJCQgJCcFXX31lhE9IRHLEoEdEpGctWrTA2bNnUadOHdSrV0/jVqVKFdy5cweJiYn44IMP0LFjRzRq1Aj37t3T+vw+Pj4YO3YstmzZgnfeeQdr16414KchIjlj0CMi0rPx48fj7t27GDhwII4ePYrLly9j586dGDlyJBQKBVxdXVGtWjWsWbMGFy9exO7duzF16lStzh0ZGYm///4bycnJOHHiBHbv3o1GjRoZ+BMRkVwx6BER6Zm3tzcOHjwIhUKBLl26IDAwEJMnT4aLiwssLCxgYWGBn3/+GXFxcQgMDMSUKVPw2WefaXVuhUKB8ePHo1GjRujatSsaNmyIFStWGPgTEZFccdYtERERkZlijx4RERGRmWLQIyIiIjJTDHpEREREZopBj4iIiMhMMegRERERmSkGPSIiIiIzxaBHREREZKYY9IiIiIjMFIMeERERkZli0CMiIiIyUwx6RERERGaKQY+IiIjITP0/v300VTzZu8IAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(2, 1, figsize=(7, 7), sharex=True)\n", + "\n", + "colors = ['orange', 'blue', 'black', 'red']\n", + "for i, run_name in enumerate(results_sigma8_squared_bias.keys()):\n", + " axs[0].plot(years, results_sigma8_squared_bias[run_name], label=run_name, marker='o', color=colors[i])\n", + "axs[0].legend()\n", + "\n", + "results_sigma8_squared_bias.keys()\n", + "axs[1].set_xlabel('Years')\n", + "axs[0].set_ylabel('Bias in sigma8 in units of sigmas')\n", + "\n", + "run_name_ = list(results_sigma8_squared_bias.keys())[0]\n", + "for i, run_name in enumerate(results_sigma8_squared_bias.keys()):\n", + " axs[1].plot(years, np.array(years)*0, ls='--', c='orange')\n", + " if run_name != run_name_:\n", + " axs[1].plot(years, results_sigma8_squared_bias[run_name]-results_sigma8_squared_bias[run_name_], label=run_name, marker='o', color=colors[i])\n", + "axs[1].set_ylabel('Top panel minus '+run_name_)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demo of AreaAtRisk metric" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Healpix slicer using NSIDE=64, approximate resolution 54.967783 arcminutes\n", + "year 1\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/maps/dust_map.py:46: UserWarning: Slicer value of nside 64 different from map value 128, using slicer value\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 2\n", + "year 3\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n", + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n", + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 4\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/maps/dust_map.py:46: UserWarning: Slicer value of nside 64 different from map value 128, using slicer value\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 5\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n", + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n", + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 6\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/maps/dust_map.py:46: UserWarning: Slicer value of nside 64 different from map value 128, using slicer value\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 7\n", + "year 8\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n", + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n", + "/Users/yoachim/anaconda3/envs/rubin12/lib/python3.12/site-packages/sklearn/utils/validation.py:2739: UserWarning: X does not have valid feature names, but StandardScaler was fitted with feature names\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 9\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/maps/dust_map.py:46: UserWarning: Slicer value of nside 64 different from map value 128, using slicer value\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 10\n" + ] + } + ], + "source": [ + "from rubin_sim.maf.metrics.cosmology_summary_metrics import UniformAreaFoMFractionMetric\n", + "from rubin_sim.maf.metrics.uniformity_metrics import NestedRIZExptimeExgalM5Metric\n", + "from rubin_sim.maf.metrics.weak_lensing_systematics_metric import RIZDetectionCoaddExposureTime, ExgalM5WithCuts\n", + "from rubin_sim.maf.metrics.exgal_m5 import ExgalM5\n", + "\n", + "nside = 64\n", + "\n", + "sim_list = [baseline_file]\n", + "name_list = [splitext(basename(sim))[0] for sim in sim_list]\n", + "\n", + "years = range(1, 11) \n", + "surveyAreas = SkyAreaGenerator(nside=nside)\n", + "map_footprints, map_labels = surveyAreas.return_maps()\n", + "slicer = maf.HealpixSubsetSlicer(\n", + " nside=nside, \n", + " hpid=np.where(map_labels == \"lowdust\")[0],#hpid=np.arange(hp.nside2npix(nside)), \n", + " use_cache=False)\n", + "\n", + "results_allruns= {}\n", + "for opsim_fname, run_name in zip(sim_list, name_list):\n", + "\n", + " results_allyears = np.zeros((len(years), ))\n", + " # loop over years \n", + " for iy, year in enumerate(years):\n", + " print('year', year)\n", + "\n", + " days = year*365.25\n", + " constraint_str = 'scheduler_note not like \"DD%\" and night <= XX and scheduler_note not like \"twilight_near_sun\" '\n", + " constraint_str = constraint_str.replace('XX','%d'%days)\n", + "\n", + " \n", + " metric = NestedRIZExptimeExgalM5Metric(\n", + " depth_cut=25.0 # what depth cuts to apply year after year?\n", + " ) \n", + " summary_metrics = [UniformAreaFoMFractionMetric(\n", + " year,\n", + " nside=nside,\n", + " verbose=False, \n", + " metric_name='FoMRatio'\n", + " )]\n", + " depth_map_bundles = [maf.MetricBundle(\n", + " metric=metric,\n", + " slicer=slicer,\n", + " constraint=constraint_str,\n", + " run_name=run_name,\n", + " summary_metrics=summary_metrics\n", + " )]\n", + " bd = maf.metricBundles.make_bundles_dict_from_list(depth_map_bundles)\n", + " bgroup = maf.MetricBundleGroup(bd, opsim_fname)\n", + " bgroup.run_all()\n", + " results_allyears[iy] = bd[list(bd.keys())[0]].summary_values['FoMRatio']\n", + " \n", + " results_allruns[run_name] = results_allyears\n", + "\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'FOM ratio')" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAFzCAYAAACHCIXLAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZdVJREFUeJzt3Xl8VNX5P/DPzGTfJhskIAHCIgQiiwlLglitGtwQ65ZWxaKgUqwS+bZfm687Wim2QtygRQXcwSq0tD+0YFsFZI+EXQQJJEBiyB7IPrm/P07uJJN1JpmZc+/M5/16zSuTyc2dZ5gh88w55zmPQVEUBURERESkeUbZARARERGRfZi4EREREekEEzciIiIinWDiRkRERKQTTNyIiIiIdIKJGxEREZFOMHEjIiIi0gkmbkREREQ64SM7AC1qamrCuXPnEBoaCoPBIDscIiIi8nCKoqCqqgr9+/eH0dj5uBoTtw6cO3cOcXFxssMgIiIiL5Ofn48BAwZ0+nMmbh0IDQ0FIP7xwsLCJEdDREREnq6yshJxcXHWHKQzTNw6oE6PhoWFMXEjIiIit+luiRaLE4iIiIh0gokbERERkU4wcSMiIiLSCSZuRERERDrBxI2IiIhIJ5i4EREREemE1MRty5YtmD59Ovr37w+DwYC//e1v3f7O119/jaSkJAQEBGDIkCH485//3O6Yzz77DKNGjYK/vz9GjRqF9evXuyD6XmiyAD9+BZz6WHxtssiOyHk8+bF5Kj5n+uTJz5snPzZP5cnPmcYem9R93C5evIixY8fi/vvvx+23397t8bm5ubjxxhvx4IMP4oMPPsA333yDefPmoU+fPtbf37FjB9LT0/HCCy/gZz/7GdavX4+77roL27Ztw6RJk1z9kLqXvw7Ing9Un2m5LWgAkPQqEHebvLicwZMfm6fic6ZPnvy8efJj81Se/Jxp8LEZFEVRpNxzGwaDAevXr8ett97a6TFPPPEENmzYgKNHj1pvmzt3Lvbv348dO3YAANLT01FZWYnPP//cesz111+PiIgIfPzxx3bFUllZCbPZjIqKCuduwJu/Dth6B4C2/+TNm+1N/VS/L3JPfmyeis+ZPnny8+bJj81TefJz5ubHZm/uoavOCTt27EBaWprNbdOmTcM777yDhoYG+Pr6YseOHXj88cfbHZOVleXGSDvQZBFZe7sXAFpu2/0QoFgAg8mdkfWeYgF2/wqdPzYDkJ0BXDIDMOrssXmqbl+PfM40iX9H+JrUEr4epbwedZW4FRYWIiYmxua2mJgYNDY2ori4GP369ev0mMLCwk7PW1dXh7q6Ouv3lZWVzg0cAM5vtR1q7TCQEmDbXc6/b+kUoDpf/BvEXCU7GALseD3yOdMk/h3ha1JL+HqU8nrUVeIGtO/hpc70tr69o2O66v21aNEiPP/8806MsgM1BfYdF3opENDHtbE4W+15oOr77o+z99+AXM/e54LPmbbw7whfk1rC16OU16OuErfY2Nh2I2dFRUXw8fFBVFRUl8e0HYVrLTMzEwsWLLB+X1lZibi4OCdGDiCwn33HTfyL/j5N/vgV8O+ruz/O3n8Dcj17nws+Z9rCvyN8TWoJX49SXo+62sctJSUFmzdvtrlt06ZNSE5Ohq+vb5fHpKamdnpef39/hIWF2Vycrs9UUYmCzkb+DEBQnDhObzz5sXkqPmf65MnPmyc/Nk/lyc+Zhh+b1MTtwoULyMnJQU5ODgCx3UdOTg7y8vIAiJGw++67z3r83Llzcfr0aSxYsABHjx7FypUr8c477+A3v/mN9Zj58+dj06ZNWLx4Mb777jssXrwYX375JTIyMtz50NozmkT5MID2L4Tm75Oy9Lno1pMfm6eyec7a4nOmWZ78vPHviP7w9SjlsUlN3Pbu3Yvx48dj/PjxAIAFCxZg/PjxeOaZZwAABQUF1iQOAOLj47Fx40Z89dVXGDduHF544QW89tprNnvApaamYs2aNVi1ahXGjBmD1atXY+3atdrYwy3uNlE+HHSJ7e1BA/RdMg149mPzVOpzZvC1vd0vks+ZlsXdBkx8u/3tnvB/jX9H9Ed9zkwBtrd7wnOm0dejZvZx0xKX7eOmarKISpSaAjE/3meqPj+RdKTJAhR8AXx9s/j+1rNAUH+5MVHnGi4Afw0DoACRSUBpNjByAXD5K7Ijo66c3Qh8fRMQNBAY9wfP/Dvyw1vAnl8BvmbgtmLApKsl2d7nb4OA6jxg9FNA7DWe93p0w3u2R+7j5jGMJv0t1LSX0QRcchNgTgQqDgElu4GgW2VHRZ0p3QNAEWs1hs8Dds0GyvbJjoq6U5otvva9Ehj8C7mxuILRBAyZBez9NdBQAdQVNq83Ik2qKxFJGwAk/A/gFy41HKfT2Hu2rooTSEeiU8TX4u1y46CunW9+fqJTxIgbAJR+C3AgXtvKmhM39TnzRKYAwDxaXC/9Vm4s1DX1w17IUM9L2jSIiRu5hjVx2yE3Duqa+vxEpwDmUYDRX4xwXPhBblzUtVIvSNwAIPJy8ZWJm7apz4/6fJFLMXEj1+jTvP1K6V7AUi83FuqYogAlO8X16BTA6AtEjBXfq4kBaU9tUfNu9QYgYrzsaFwrojkRKGPipmlq4hbBxM0dmLiRa4ReKqoTLbVA+X7Z0VBHqo6LtSlG/5YEwDpdysRNs9TnJmwE4BsiNxZXs4648fWoaWUccXMnJm7kGgYDED1ZXOd0qTapz0tUMmDyE9eZuGmft0yTAkD4WAAGoOYcUNN5v2mSqL5CfAgEPH8EWCOYuJHrqOvczrNAQZOKWxUmqFigoH3elLj5hoiRRYDVzlpVliO+BsXprx+pTjFxI9dhgYK2tS5MUJlHNxcolAMXTkoJi7rhTYkb0LJuigUK2sRpUrdj4kauEzURMBjF/j7V52RHQ601VALlh8T11omb0RcIHyOul+51f1zUtdoioDofXlGYoFITVBYoaJO1MMFLPkhoABM3ch3fUMB8mbjOUTdtKdkNQAGCB4mdwFvjOjftshYmXCr+f3kDbgmibRxxczsmbuRanC7VpvPqNGlq+58xcdMu9TnxptGNiHHi68VTQF2pzEiorcaLQOV34joTN7dh4kauxcRNmzpa36ZigYJ2edv6NkDsxB8yVFxngYK2lB0AlCYgILb9yD25DBM3ci01MSjdC1jq5MZCgtLUdeJmHg0Y/VigoEXemLgBnC7VKk6TSsHEjVwrdBjgHw001fPTslZUHhNJmSmwpVNCaya/VgUKnC7VjNrzzYUJACK9pDBBxQ4K2sSOCVIwcSPXMhg4Xao16vMQmSyqSDvCdW7aoz4XoZcCvmFyY3E3dlDQJo64ScHEjVyPiZu2dDVNqmLipj3eOk0KtGx9UnVcbGVD8llqW7YUYuLmVkzcyPWYuGmL+jz06aCiVNV67ywWKGiDNyduAX3EzvxAy079JFf5IUBpFD2pgwbKjsarMHEj14uaABhMQPUZ4GK+7Gi8W305UHFYXO9qxM2cKAoU6suAi7luCY264c2JG8ACBa1pPU1qMMiNxcswcSPX8wluWezOUTe5ineJryFDgIC+nR9n8gPCmzdP5nSpfLXFogMJ4D0dE9qKaLVNDcmnPg/e+kFCIiZu5B7qRq9M3OSyZ32bKjJZfGXiJp+1MGE44GeWG4sskaws1RRWlErDxI3cg+vctMGhxI0FCppR5uXTpEBL4lZ5FGislhuLt2tqAMoPiOssTHA7Jm7kHn2aE4Wyb0U1Ermf0gSU7BTXO2p11VbrxI0FCnJ5+/o2QOzMHxArXsdq0kByVBwFmurEtjQhQ2RH43WYuJF7BMeLNVVNDRzBkaXiiNhKwSe4Zf1aV2wKFE65PDzqAhM3gQUK2qBOV0eMBwxMI9yN/+LkHtyIVz7rxrsTAKNP98ezQEEb6kqAi6fFdW9fT8QOCtrA9W1SMXEj92GBglyOrG9TcZ2bfCxMaMEOCtrAjglSMXEj92k94sY1U+7HxE2fOE3aQk0Uyg8Bljq5sXirJgtQ2tx3mombFEzcyH0ikwGDD1BT0LInFblHXSlQ+Z243tPEjcm2HEzcWgQNFDv1K41AxSHZ0Xinqu8BSzVgCgRCR8iOxisxcSP38QkEIsaJ6+c5XepWJc0b74YOBwKi7f89c6JoRF9f2rLOityLiVsLg4EFCrJZ17eNA4wmqaF4KyZu5F7W6dLtcuPwNueb/70dGW0DAJM/YGaBgjR1JS0VvVwILkSyg4JUZeyYIBsTN3IvVpbK0ZP1bSquc5NHTU5ChrEwQcXKUrlYUSodEzdyrz7NlaVlOUBjjdRQvEaTpWWqlImbvnCatD1r66v9Yl9Ich+liRWlGsDEjdwraKDYAV1pBEr3yo7GO1QcBhovAD4hYs2ao9SkoYwFCm7HxK29kCFix/6mupaCG3KPC7liE2+jH2AeJTsar8XEjdyLG/G6n/rvHDWpZ4uJwy8TBQp1JawGdjcmbu0ZjGLHfoDr3NxNHW0LHyP+JpAU0hO3ZcuWIT4+HgEBAUhKSsLWrVu7PP7NN99EQkICAgMDMWLECLz33ns2P1+9ejUMBkO7S20t+2NqBgsU3Ku4h4UJKpN/y0gdp0vdp64UuJgrrnNaylYEN+KVopTTpFpgR98b11m7di0yMjKwbNkyTJkyBX/5y19www034MiRIxg4cGC745cvX47MzEy89dZbmDBhAnbv3o0HH3wQERERmD59uvW4sLAwHDt2zOZ3AwICXP54yE5tN+I1GOTG4+l6U5igikwCyvaJN8q425wTF3VNHd0IGQr4hUsNRXMiWaAgBQsTNEHqiNuSJUswe/ZszJkzBwkJCcjKykJcXByWL1/e4fHvv/8+Hn74YaSnp2PIkCH4+c9/jtmzZ2Px4sU2xxkMBsTGxtpcSEMik8Qwe21Ry4gCuUZtMVB1XFyPntzz87BAwf04Tdo5a+KWI4pvyPUURaxzBTjiJpm0xK2+vh7Z2dlIS0uzuT0tLQ3bt3c8hVZXV9du5CwwMBC7d+9GQ0NLddGFCxcwaNAgDBgwADfffDP27dvXZSx1dXWorKy0uZALmQJaPrFxI17XKtkpvoaNAPwje34edlBwPyZunQsdIXbub7zY8sGEXKs6X6xzNZjEuleSRlriVlxcDIvFgpiYGJvbY2JiUFhY2OHvTJs2DW+//Tays7OhKAr27t2LlStXoqGhAcXFxQCAkSNHYvXq1diwYQM+/vhjBAQEYMqUKTh+vPP/3IsWLYLZbLZe4uLinPdAqWMsUHAP6zRpau/OE36ZaFdWV8wCBXcpaa66ZuLWntHU0oWF06XuoU6TmkeLD98kjfTiBEOb9U2KorS7TfX000/jhhtuwOTJk+Hr64sZM2Zg1qxZAACTSVTLTZ48Gffeey/Gjh2LqVOn4pNPPsGll16K119/vdMYMjMzUVFRYb3k5+c758FR55i4uYcz1rcB4g91OAsU3IaFCd1jBwX3YscEzZCWuEVHR8NkMrUbXSsqKmo3CqcKDAzEypUrUV1djVOnTiEvLw+DBw9GaGgooqM77r9oNBoxYcKELkfc/P39ERYWZnMhF1MTifL9YrqDnK+pESjuxca7bXGdm/tYCxOGAH4RcmPRKnZQcC8WJmiGtMTNz88PSUlJ2Lx5s83tmzdvRmpq19M6vr6+GDBgAEwmE9asWYObb74ZRmPHD0VRFOTk5KBfv35Oi52cIDgOCBoAKBagZI/saDxT+UHAUi02K3XGZplM3NyH69u617rZPNdduh47JmiG1O1AFixYgJkzZyI5ORkpKSlYsWIF8vLyMHfuXABiCvPs2bPWvdq+//577N69G5MmTUJZWRmWLFmCQ4cO4d1337We8/nnn8fkyZMxfPhwVFZW4rXXXkNOTg7efPNNKY+RuhCdAuT9VUznxVwlOxrP03rjXYMTPqNFtClQ4DYursPErXvmUWIH/4YKMa0cMkR2RJ6rpkBcYAAixsqOxutJTdzS09NRUlKChQsXoqCgAImJidi4cSMGDRoEACgoKEBeXstCaIvFgldeeQXHjh2Dr68vrr76amzfvh2DBw+2HlNeXo6HHnoIhYWFMJvNGD9+PLZs2YKJEye6++FRd1onbuR8zlrfpooY06pAIR8Ibr/XIjkJE7fuGX3FDv6le8WoGxM31ylt3pkhbCTgEyw3FoJBUTjG3FZlZSXMZjMqKiq43s2VincCm1IA/2jgtiKO4DjbhmHAhR+Aq74A+k9zzjk3jhPrEqeuA+J+5pxzkq36MuDT5q1bbi/p3TYunm73w8CJFcCo3wHjFsmOxnMdehE48DQw+B4g9QPZ0Xgse3MP6VWl5MUixoupjrpioOqE7Gg8S22RSNoAIHqS887LdW6upy4CD45n0tadiFbr3Mh1WJigKUzcSB6Tf0siwOlS51L/Pc2jnNsuiYmb63Ga1H6tW19x8sh1StkxQUuYuJFc6sawTNycy9nr21TsoOB6TNzsF36Z2Mm/rhioPiM7Gs9U22rTbXXTY5KKiRvJxY14XcNViVv4mOY3yvN8o3QVJm72MwWInfwB7ufmKmXNhQkhQ507ek89xsSN5FITi4qDQEOV3Fg8RVNDy954vW111ZZPYMsbJadLna++vGVtIqel7MMOCq7Fjgmaw8SN5ArqDwQNBJQmoGS37Gg8Q9l+wFID+IaL5vLOxnVurmMtTBgM+EdJDUU3WKDgWqXceFdrmLiRfJwudS7rNOlk52y82xYTN9fhNKnjItn6yqVYUao5TNxIvj4sUHAqV61vU0Umi69lLFBwOmviliw3Dj2JGAvAANScA2oKuz2cHFBfAVxo3qopYrzcWMiKiRvJZx1x28lEwBlcnbipBQq1RUDNWdfch7fiiJvjfILFjv5Ay0J6co6yHPE1aCAQEC01FGrBxI3kCx8rqsPqS4Gq72VHo281hcDFUwAMzt14tzUWKLhG69ENridyjLXhPF+PTsXG8prExI3kM/m1TA1xurR31H+/8ETA14Xt2rjOzfnKWJjQYyxQcA2ub9MkJm6kDeq03vntcuPQu+Lmfz9XTZOqmLg5H6dJe44FCq7BjgmaxMSNtIGVpc7h6vVtKnZQcD4mbj2nLpy/eBqoK5Ebi6dovAhUfieuM3HTFCZupA3WjXgPi7U+5DhLPVCyV1x3deIWPra5QOFHUc1HvcfEref8zEDIMHGdBQrOUbYfgAIE9hMX0gwmbqQNgbFAcDwAhRvx9lRZDtBUB/hFAqGXuva+fAJFA3uA06XOUF8BVB0X15m49Uwk17k5Fde3aRYTN9IOTpf2TutpUoPB9ffHdW7Oo44SBQ9iYUJPMXFzLlaUahYTN9IOa+LGAoUecVdhgiqCiZvTcJq09yJYoOBUHHHTLCZupB19Wm/E2yQ3Fj1yV2GCiiNuzsPErffUAoWq40BDpdxY9M5SK9YbAxxx0yAmbqQd4WMAUxDQUNFSzUT2qT4LVOeL3qRRE91znxFjxf3VFgLVLFDoFTVxi2Di1mMB0WKHf6Blx3/qmfJDgNIopu2D4mRHQ20wcSPtMPoCURPEda5zc4z672W+DPANcc99+gQBYSxQ6LWGypaOIRxx6x12UHCOslbTpO5YL0sOYeJG2sIChZ5R/736pLr3fjld2nulzYUJ7AfZe+yg4BylLEzQMiZupC1M3HrG3evbVEzceo/r25yHHRScwzp1z8RNi5i4kbZETxZfK44A9WVyY9ELS13LH1ppidte996vJ1H/7Zi49Z6auFV+J3b+J8c1NQDlB8R1jrhpEhM30paAvi07oBfvkhuLXpR+CzTVA/7RQMhQ9953xDgWKPQWR9ycR93lX2kCyg7IjkafKo6Ivye+Zvf/PSG7MHEj7eF0qWPcvfFuaz5BQFiCuM7pUsexMMH5uJ9b71j3bxvPwgSNYuJG2tOHiZtDZK1vU3GdW89ZCxPigIA+cmPxFOyg0DvsmKB5TNxIe9QEpGQX0GSRG4vWKUqrjglurihVMXHrOU6TOh9H3HqHHRM0j4kbaY85EfAJFtNIlUdkR6Nt1flAzTnAYAKikuXEoCYdZUzcHMbEzfnUkaLyQ6Jwh+zXZGnZvJgjbprFxI20x+jTsvs/p0u7pv77hI8Vya4MaoFCTYG4kP3KmLg5XVCc2PFfaQTKD8qORl+qvgcs1aKDTeilsqOhTjBxI21Sp/2YuHVN9vo2QCSMYSPFdU6X2q+hCqhkYYLTGQycLu0p6zTpOMBokhoKdY6JG2kTK0vto4XEDWjpscnEzX5l+wAoQNAAsQ0OOQ8LFHqGhQm6wMSNtEndiLfyGFBXKjcWrbLUNr/5w/2trtpigYLjuL7Nddj6qmfYMUEXpCduy5YtQ3x8PAICApCUlIStW7d2efybb76JhIQEBAYGYsSIEXjvvffaHfPZZ59h1KhR8Pf3x6hRo7B+/XpXhU+u4h/VssaieKfcWLSqNFvsch4QAwQPlhsLEzfHWd8kmbg5nbVA4YD4P0LdU5paPghyxE3TpCZua9euRUZGBp588kns27cPU6dOxQ033IC8vLwOj1++fDkyMzPx3HPP4fDhw3j++efxyCOP4B//+If1mB07diA9PR0zZ87E/v37MXPmTNx1113YtYu78OuOdbp0u9w4tOq8ug2IhI1324oYB8AgKlxrCuXGohcccXOdkKFi5/+mOqDiqOxo9OHCSVHJb/QHzKNkR0NdkJq4LVmyBLNnz8acOXOQkJCArKwsxMXFYfny5R0e//777+Phhx9Geno6hgwZgp///OeYPXs2Fi9ebD0mKysL1113HTIzMzFy5EhkZmbimmuuQVZWlpseFTlNHxYodEkr69sAwDeEBQqOaKgSywAAJm6uYDCInf8BFijYS51WDh8DGH3lxkJdkpa41dfXIzs7G2lpaTa3p6WlYfv2jkdY6urqEBAQYHNbYGAgdu/ejYYGMRy+Y8eOduecNm1ap+dUz1tZWWlzIQ2wbsS7mxvxtqUo2krcAE6XOqIsB9bChMAY2dF4JhYoOIaFCbohLXErLi6GxWJBTIztH62YmBgUFnY81TJt2jS8/fbbyM7OhqIo2Lt3L1auXImGhgYUFxcDAAoLCx06JwAsWrQIZrPZeomLi+vloyOnCBsF+IQCjReAikOyo9GWi6dFY3eDDxApaePdtpi42Y/TpK7HLUEcU8rETS+kFycY2qzNURSl3W2qp59+GjfccAMmT54MX19fzJgxA7NmzQIAmEwte844ck4AyMzMREVFhfWSn5/fw0dDTmU0AdGTxHVOl9pS/z0ixgM+gXJjUTFxsx8LE1xPTUDKcjhi3x1FaUlwWVGqedISt+joaJhMpnYjYUVFRe1GzFSBgYFYuXIlqqurcerUKeTl5WHw4MEIDQ1FdHQ0ACA2NtahcwKAv78/wsLCbC6kEeo04HkWKNgoblWYoBUR4yEKFM4CNT/KjkbbOOLmeqGXig4AjRdFRwDqXHU+UFciRvDDE2VHQ92Qlrj5+fkhKSkJmzdvtrl98+bNSE3tek8qX19fDBgwACaTCWvWrMHNN98Mo1E8lJSUlHbn3LRpU7fnJI1iB4WOaW19G8ACBXs1XAAqvxPXmbi5jtHUXO0MrnPrjvrvYx4NmAK6Ppak85F55wsWLMDMmTORnJyMlJQUrFixAnl5eZg7dy4AMYV59uxZ615t33//PXbv3o1JkyahrKwMS5YswaFDh/Duu+9azzl//nxceeWVWLx4MWbMmIG///3v+PLLL7Ft2zYpj5F6SZ0qvXACqD0PBPSRG48WNFYDZfvF9T4aStwAkYhUHhWJ2yU3yo5Gm9TChMBLWJjgapGXi9Hpsm+B+HtkR6NdLEzQFamJW3p6OkpKSrBw4UIUFBQgMTERGzduxKBBgwAABQUFNnu6WSwWvPLKKzh27Bh8fX1x9dVXY/v27Rg8eLD1mNTUVKxZswZPPfUUnn76aQwdOhRr167FpEmT3P3wyBn8IoCwBJEMFO8EBkyXHZF8pXtFA+3AfkDQQNnR2IpMAk590NI8ndrjNKn7sIOCfdgxQVekJm4AMG/ePMybN6/Dn61evdrm+4SEBOzbt6/bc95xxx244447nBEeaUF0SnPitoOJG9BqmjRV/sa7bbFAoXtM3NwnslVlqdIEGKTX42kTK0p1ha9i0j52ULB1XoOFCSq1QKH6DFBbJDsabSpj4uY25lGiE0BDJXAhV3Y02lRT0Ly1kBGIGCs7GrIDEzfSPutGvHuApka5scimxY13W/MNAcJGiOscdWuv8SILE9zJ6Cs6AQDcz60z6mhb2EjAJ1huLGQXJm6kfeYE0XfQUi2aRnuzCyeBuvPiDUmr0xqcLu1cWY6YsgvsDwTGyo7GO7CDQtdKuX+b3jBxI+0zGIHoyeK6t28LYt1493Ltlu0zcesc17e5HxO3rrGiVHeYuJE+WNe5MXEDoM1pUhUTt84xcXO/1q2vFEVuLFrEETfdYeJG+sDETVAffx8NbyhtLVDIF3vvUQsmbu4Xnig6AtQVi9cktagtBqqbt9xSNysmzWPiRvoQNQmAQazx8tZ2Sg0XgPLmjXe1POLmGwqEXSquc9StReNFsa0NwMTNnUwBoiMAwOnStsqat9cKGQb4meXGQnZj4kb64Gdu+ePrraNupXvEwvagAeKiZRGcLm2nbH9zYUI/cSH3ab2fG7Xg+jZdYuJG+uHt06V6WN+msq5z2ys3Di1R/y0iONrmduyg0DHr1D0TNz1h4kb64e2J23k9Jm4ccbPi+jZ5OOLWMRYm6BITN9IPNWEp3QM0NciNxd0UBShp1epK6yLHi68sUGjBxE2eiLFiW6GaAnEhoL4cuPCDuM4RN11h4kb6EXapaDpvqRUbmXqTquNAXYlo3xMxXnY03fMNA0JZoGDFwgS5fIJFZwAAKO2+37VXUP+GBg8C/KOkhkKOYeJG+mEweu90qfp4I5MAk5/cWOzF6dIWamFCQCwQ1F92NN4pgtOlNjhNqltM3EhfvD1x08P6NhUTtxacJpWPHRRssaJUt5i4kb4wcZMbhyOYuLVg4iaftbKUr0cAHHHTMZ+e/FJ5eTneeecdHD16FAaDAQkJCZg9ezbMZm7gRy4WNVFMmV48LRYZe8N+WA1VQMUhcV1PiZu6Fq86T+zQHhAtNx6ZmLjJp3YG4Ouxec3ld+I6R9x0x+ERt71792Lo0KFYunQpSktLUVxcjKVLl2Lo0KH49lsOQZOL+YYC5kRx3VtG3Up2i/VRwYP0tT7KzwyEDhfXvXmUo7EaqDwirjNxk8fPLDoEAC0dA7xV2X4ASvNm0LGyoyEHOZy4Pf7447jllltw6tQprFu3DuvXr0dubi5uvvlmZGRkuCBEojbUUafz2+XG4S7q49TTaJtKTVTKvDhxsxYmxACBOkq8PRH3cxM4TaprPRpxe+KJJ+Dj0zLL6uPjg//93//F3r3cJZ3cQN3HzFtG3PS4vk3FdW6206QGg9xYvB0LFIQydkzQM4cTt7CwMOTl5bW7PT8/H6GhoU4JiqhL1o14swFLvdxYXE1pAkp2iutM3PSpjOvbNIOtrwSOuOmaw4lbeno6Zs+ejbVr1yI/Px9nzpzBmjVrMGfOHPziF79wRYxEtkKHAf7RQFOd569VqfweqC8DTAFA+FjZ0ThOfWO4eFpsIOyNrCNuyXLjoJYRpgsngPoKubHIYqkFKg6L6/wwoUsOV5X+6U9/gsFgwH333YfGxkYAgK+vL371q1/hD3/4g9MDJGrHYACiJgPn/immEaMnyY7Idawb707Qz8a7rakLwi+cEAlMvzTZEblXYw1QwcIEzfCPEkU+F0+LzgExP5EdkfuVHwQUi/jwGzRAdjTUAw6PuPn5+eHVV19FWVkZcnJysG/fPpSWlmLp0qXw9/d3RYxE7fVR93Pz8AKFYh0XJqi8ebq0fL94k2RhgnZ4eweF1tOkXHOpSz3egDcoKAiXXXYZxowZg6CgIGfGRNQ9bylQ0HNhgsqbEzcWJmiPtxcosGOC7tk1VXrbbbdh9erVCAsLw2233dblsevWrXNKYERdipoAGExA9Rlx8cQh//qKlmk2Jm76xI13tcfbOyiUMnHTO7sSN7PZDEPzp8WwsDDrdSJpfIKB8DGiOKF4BzDwTtkROV/JLgAKEBwPBMbIjqbn1DeIi6dEgYJ/lNRw3IqJm/aor8fK70QHAZ9gufG4U1MDUH5AXGdFqW7ZlbitWrXKen316tWuioXIMdEpInE776GJmzpN2idVbhy95RcOhAwFLvwgPu33u052RO7RWMPqPS0KjBUdA2oKxObIev//5YiKI0BTPeBrBkKGyI6GesjhNW4//elPUV5e3u72yspK/PSnP3VGTET2ifbwAgU9d0xoyxunS8sPNBcm9AUCL5EdDbXmrfu5WQsTxnPNpY45nLh99dVXqK9vv+lpbW0ttm7d6pSgiOyiflIu+1bsTeRJlKbmqVJ4SOLWvIeZNyVu6mONYGGC5nhr66tSdkzwBHbv43bgwAHr9SNHjqCwsND6vcViwRdffIFLLuGnSnKj4HgxmlFbJD5JetKUR8VRoKECMAWJtXx6540jblzfpl3eOuJWxo4JnsDuxG3cuHEwGAwwGAwdTokGBgbi9ddfd2pwRF0yGMRo1Jm/i/VgnpS4qevboiYARof3ydYea4FCLlBXCvhHyo3HHZi4aZf6nFQcFqP1pgC58bhDk0VsOgzwNalzdr8j5ObmQlEUDBkyBLt370afPn2sP/Pz80Pfvn1hMplcEiRRp1onbp7EE/Zva611gULZt0DstbIjci22FdK2oAGic0BdMVB+CIjygnZkVccAS42oog0dLjsa6gW717gNGjQIgwcPRlNTE5KTkzFo0CDrpV+/fj1O2pYtW4b4+HgEBAQgKSmp23VyH374IcaOHYugoCD069cP999/P0pKWnogrl692joy2PpSW+tha6BIsBYo7AAURW4szmRN3DxoFNGbpkvLDgBKI+DfxzP3GNQ7g8H7OihYCxPGAUYOsuhZj+dgjhw5gry8vHaFCrfccovd51i7di0yMjKwbNkyTJkyBX/5y19www034MiRIxg4cGC747dt24b77rsPS5cuxfTp03H27FnMnTsXc+bMwfr1663HhYWF4dixYza/GxDgBUPh3igyGTD4ADXngOo80YdQ7+pKgcqj4nr0ZLmxOFNkEpD3iZckbuyYoHmRlwOFm7xnnVsp17d5CocTt5MnT+JnP/sZDh48CIPBAKV5lEPdlNdisdh9riVLlmD27NmYM2cOACArKwv/+te/sHz5cixatKjd8Tt37sTgwYPx2GOPAQDi4+Px8MMP4+WXX7Y5zmAwIDY21tGHRnrkEyQ+QZbuFfu5eULiplaThgwDAvp0fayeeNOIG9e3aV+kl3VQYKsrj+HwdiDz589HfHw8fvzxRwQFBeHw4cPYsmULkpOT8dVXX9l9nvr6emRnZyMtLc3m9rS0NGzf3vG+XKmpqThz5gw2btwIRVHw448/4tNPP8VNN91kc9yFCxcwaNAgDBgwADfffDP27dvXZSx1dXWorKy0uZCOtJ4u9QSetr5Npb5hXDgJ1JfJjcXVmLhpnzryVH5AdBTwZEqT2Kwc4IibB3A4cduxYwcWLlyIPn36wGg0wmg04oorrsCiRYusI2H2KC4uhsViQUyMbSufmJgYm61GWktNTcWHH36I9PR0+Pn5ITY2FuHh4TbVrCNHjsTq1auxYcMGfPzxxwgICMCUKVNw/PjxTmNZtGgRzGaz9RIXF2f34yAN8NTErY+HJW5+ES27tXvy9JSlVix4B5i4aVnIENFBoKm+pSewp7pwEmioBIz+gDlBdjTUSw4nbhaLBSEhIQCA6OhonDt3DoAoXmi7rswebfueKorSaS/UI0eO4LHHHsMzzzyD7OxsfPHFF8jNzcXcuXOtx0yePBn33nsvxo4di6lTp+KTTz7BpZde2uVWJZmZmaioqLBe8vPzHX4cJJGauJXtE22G9KzJAhSrG+96UGGCyjpdulduHK5kLUyIBoL4IVCzDAbRQQDw7A8SQMvjCx8DGH3lxkK95vAat8TERBw4cABDhgzBpEmT8PLLL8PPzw8rVqzAkCH29z6Ljo6GyWRqN7pWVFTUbhROtWjRIkyZMgW//e1vAQBjxoxBcHAwpk6dihdffBH9+vVr9ztGoxETJkzocsTN398f/v7+dsdOGhM8CAiIBWoLRULQd6rsiHqu4jDQWAX4hADmRNnROF9kEpD3V89eV8TCBP2IvBwo+qp5/df9sqNxHXZM8CgOj7g99dRTaGpqAgC8+OKLOH36NKZOnYqNGzfitddes/s8fn5+SEpKwubNm21u37x5M1JTOx5pqK6uhtFoG7K6DYnSyVYQiqIgJyenw6SOPITB0LL5rt6nS60b7070zJJ9byhQ4Po2/fCWDgosTPAoDo+4TZs2zXp9yJAhOHLkCEpLSxEREdHpFGdnFixYgJkzZyI5ORkpKSlYsWIF8vLyrFOfmZmZOHv2LN577z0AwPTp0/Hggw9i+fLlmDZtGgoKCpCRkYGJEyeif//+AIDnn38ekydPxvDhw1FZWYnXXnsNOTk5ePPNNx19qKQn0SlA/jrPSdw8rTBBFdGmQMEvQm48rsDETT/U56gsRyxT8MQPS4rSkpjyNekRHErcGhsbERAQgJycHCQmtkzjREb2rH1Neno6SkpKsHDhQhQUFCAxMREbN27EoEFiS4eCggLk5eVZj581axaqqqrwxhtv4H/+538QHh6On/70p1i8eLH1mPLycjz00EMoLCyE2WzG+PHjsWXLFkycOLFHMZJOtN2IV69TVJ6euPlHih6zF3PFm0nsNbIjci4WJuhL6HDRSaDxIlD1vWcu3K/OA+pLxX6Xnrj8wgsZlM7mGDsxdOhQrFu3DmPHjnVVTNJVVlbCbDajoqICYWFhssMhe1hqgb+GibL+W04CIfGyI3JcXQnwWbS4fnsx4B8lNx5X2XonkP8pMG4xMOp/ZUfjXCV7gH9NFM/dbef1+wHCm2y+Ajj/DZDyARB/j+xonC9/PbD1NrHf5Q1db41Fctmbe/RojVtmZiZKS0t7FSCRU5kCWirE9DpdWrxTfA0b4blJG+DZ69zUxxTBwgTdiPDwjXjZMcHjOLzG7bXXXsOJEyfQv39/DBo0CMHBwTY///ZbD1/kSdoVnQqU7AbObwcG3y07GscVN2887anTpCpvSNw4Taof6oJ9T+1ZysIEj+Nw4nbrrbe6IAwiJ4hOAY5l6XjEzcPXt6msHRR+AOrLAb9wmdE4FxM3/bE2m98nOgwYHJ6I0jaOuHkchxO3Z5991hVxEPWemvCU7xeLjX2Cuz5eS5oaxWgh4PmJm38UEDwYuHiquUDhp7Ijcg5LHVDBwgTdMSeIjgINlaLaOXSY7Iicp6ZA7G9pMAIRY2RHQ07iYR8tyKsFxwGBlwCKBSjR2c78FYeak81QIGyU7GhczxOnS8sPiuIYv0ixKTTpg9FXdBQAPG8/N/XxhI3U1wdZ6hITN/Iseu1bap0mneyZe0m15YmJm3WaNJmFCXrjqevcOE3qkZi4kWexJm7b5cbhqPNeUpig8ujEjdOkuhPpoR0UytjqyhMxcSPP0rr1lWNbFMrlLYUJKjW5uXACqK+QG4uzMHHTL2sHhW/19XejO+yY4JGYuJFniRgPGP2AumJRtagHtUUtsUZPkhuLu/hHtawD84TpKUsdUHFQXOebpP6YE0VngboSoDpfdjTOUXu+5bFEjJMaCjmX3VWlCxcutOu4Z555psfBEPWayV+8cRbvEBc9VIhZN95N8MzenZ2JTAIunhYjVTFXy46mdyoOsTBBz0z+QHii6Fla+i0QPFB2RL1X1twlIXQ44MsOQJ7E7sTtueeeQ//+/dG3b1901iXLYDAwcSP5olNaErf4mbKj6Z63TZOqIpOA/HWesc6t9TQpCxP0KeLy5sQtG4i7VXY0vcfCBI9ld+J2/fXX47///S+Sk5PxwAMP4KabboLJ5AXVb6Q/agJ0XicFCmohhbo+z1tEeFCBAte36V/k5cDJlZ4xdQ+wY4IHs3uN28aNG3Hy5ElMmjQJv/3tbzFgwAA88cQTOHbsmCvjI3JcdHMCVHEQaKiSG0t3mhpEY3LAO0fcAKDquP4LFJi46V+Eh1WWljJx81QOFSf069cPmZmZOHbsGNauXYuioiJMmDABU6ZMQU1NjatiJHJMUH8gaKBoX6MmRVpVfgCw1AC+4WKTTG8SEC2eJ6BlPY4eWerF5rsAEzc9ixgjOgzUFoqOA3pWX95S8BQxXmoo5Hw9riqdMGECrr76aiQkJGDfvn1oaGhwZlxEvaOXjXjPq+vbJnlej0R7eMJ+bhWHgKZ6UVgSPFh2NNRTPsEtH570PupWliO+Bg8SFdzkURx+p9ixYwcefPBBxMbG4vXXX8cvf/lLnDt3DmFhrFohDdFL4uathQkqT0jcWJjgOTxlupSFCR7N7uKEl19+GatWrUJJSQnuuecebNu2DZdddpkrYyPqudaJm6Jo9w3Vmrh5WWGCytMSN9K3yMuBUx/ov0ChlB0TPJndidvvfvc7DBw4EHfddRcMBgNWrVrV4XFLlixxWnBEPRYxDjAFAPWlQNX3QNgI2RG1V1MIXMwFYPCejXfbshYofA80VOpzvykmbp7D+kFC54lbGTsmeDK7E7crr7wSBoMBhw8f7vQYg1ZHNcj7mPxEs+/z28SolhYTN3W0zTxanwmLMwT0EQUK1XlA6T4g5ieyI3KMpV4UmAB8k/QEaoeB6jygtlgU0OhNwwWgsnm3B06VeiS7E7evvvrKhWEQuUB0SkviNmSW7Gja8/b1barIpObELVt/iVvF4VaFCfGyo6He8g0TnQaqjotK537XyY7IceX7AShAYH8gMEZ2NOQCPS5jKy4uRklJiTNjIXIurRcoMHET9LzOTY054nLtrqMkx1gLFHT4egRYmOAFHErcysvL8cgjjyA6OhoxMTHo27cvoqOj8etf/xrl5eUuCpGoh9SEqPyQWD+lJZZ6oHSvuO5tHRPaUhO3Mh2+UXJ9m+dRF/TrtUCBHRM8nt1TpaWlpUhJScHZs2dxzz33ICEhAYqi4OjRo1i9ejX+/e9/Y/v27YiI8KIm2aRtgbFiX62Lp4DiXdqa9ijLASy1oil56KWyo5FLTXoqdVigwMTN80TqfEsQdkzweHYnbgsXLoSfnx9++OEHxMTEtPtZWloaFi5ciKVLlzo9SKIei05tTtx2aCtxs06TTuYUW0AfICgOqM7XV4GCpb55PRGYuHkStdPAhR9EBwK/cJnROMZSK9ZdApwq9WB2T5X+7W9/w5/+9Kd2SRsAxMbG4uWXX8b69eudGhxRr2l1nRvXt9nS4zo3tTDBNxwIGSI7GnIW/yjRcQBo6UCgF+UHAcUC+EcDQQNkR0MuYnfiVlBQgNGjR3f688TERBQWFjolKCKn6aMmbjtF71KtYOJmS4+JW+tNTr191NTT6LWDQuvCBL4mPZbdiVt0dDROnTrV6c9zc3MRFcWeaKQx4WMAUyDQUN6yt5Fs1efE9hcGIxA1UXY02qDHAgWub/Ncei1QYMcEr2B34nb99dfjySefRH19fbuf1dXV4emnn8b111/v1OCIes3oC0RNENe1Ml1q3Xj3MsA3VG4sWmFToFAlNxZ7MXHzXHrtoMCOCV7B7uKE559/HsnJyRg+fDgeeeQRjBw5EgBw5MgRLFu2DHV1dXj//fddFihRj0WnAkVbgOLtwNAHZEcj4gA4TdpaQF+xJqf6jNj4tO+VsiPqWlMDOyZ4MnWqtPI7oPEi4BMsNx57WOrFGjeAI24ezu7EbcCAAdixYwfmzZuHzMxMKIoCQLS5uu666/DGG28gLi7OZYES9ZjWChS4vq1jkUkicSvN1n7iVnEYaKoDfM1AyFDZ0ZCzBcaIzgM154Cy/frYa7HySHOxjJldPDyc3YkbAMTHx+Pzzz9HWVkZjh8/DgAYNmwYIiMjXRIckVNETxZfK47IL++31LVMsTFxsxWRBJz5uz4KFFiY4PkiLheJW2m2PhK31vu38TXp0exe43by5EnrKFtERAQmTpyIiRMnMmkj7Qvo2zIqUrxLbixl+8SnYv9oIHSY3Fi0Rk+VpVzf5vn0VqDAVldew+7Ebfjw4Th//rz1+/T0dPz444+9DmDZsmWIj49HQEAAkpKSsHXr1i6P//DDDzF27FgEBQWhX79+uP/++9v1TP3ss88watQo+Pv7Y9SoUdxfjrQzXdp6mpSfim1ZCxSOab9AwdqjlImbx9JbBwW2uvIadidu6mibauPGjbh48WKv7nzt2rXIyMjAk08+iX379mHq1Km44YYbkJeX1+Hx27Ztw3333YfZs2fj8OHD+Otf/4o9e/Zgzpw51mN27NiB9PR0zJw5E/v378fMmTNx1113YdcuySMtJJc61aEWBshynoUJnQqMAQIvAaBoe+PTpgax7gkAopLlxkKuo45cVRwWHQm0rMnS8n+GI24ez6Em8862ZMkSzJ49G3PmzEFCQgKysrIQFxeH5cuXd3j8zp07MXjwYDz22GOIj4/HFVdcgYcffhh79+61HpOVlYXrrrsOmZmZGDlyJDIzM3HNNdcgKyvLTY+KNElNlEp2yd2Il4UJXdPDdGnFERYmeIOgAWJJg2JpqdbUqqpjgKVGVL+GDpcdDbmY3YmbwWCAoc3UTtvvHVFfX4/s7GykpaXZ3J6Wlobt2zseFUlNTcWZM2ewceNGKIqCH3/8EZ9++iluuukm6zE7duxod85p06Z1ek5A7ENXWVlpcyEPY04Uf9QaKsUbrwwX84Gas4DB1LK3HNnSQ+LGwgTvYDDop4OCdX3bOMBokhoKuZ7dVaWKomDWrFnw9/cHANTW1mLu3LkIDrbd32bdunV2na+4uBgWi6Vd79OYmJhOW2elpqbiww8/RHp6Ompra9HY2IhbbrkFr7/+uvWYwsJCh84JAIsWLcLzzz9vV9ykU0Yf0aXgx/+KUa/wRPfHoI62hY/Rx75QMugqceP6No8XmQQUbtJ+gQLXXHoVu0fcfvnLX6Jv374wm80wm82499570b9/f+v36sVRbUftFEXpdCTvyJEjeOyxx/DMM88gOzsbX3zxBXJzczF37twenxMAMjMzUVFRYb3k5+c7/DhIB2QXKHCatHvWAoXvgIYLcmPpDN8kvYdeChRYmOBV7B5xW7VqlVPvODo6GiaTqd1IWFFRUbsRM9WiRYswZcoU/Pa3vwUAjBkzBsHBwZg6dSpefPFF9OvXD7GxsQ6dEwD8/f2tI4nkwayJm6QCBWvHBB3sCSVLYGyrjU9zgL5XyI7IVlMjUN5cmMARN8+nJkLlB0RRitFXbjwdUZqA0n3iOhM3ryCtOMHPzw9JSUnYvHmzze2bN29GamrHb2zV1dUwGm1DNpnEfL5a9ZqSktLunJs2ber0nORFopo34q08BtSVuve+LbViDzcA6MMRty5pebq04oh4Ln3DgFAWJni84HhRhNJUL6pLtajqB6CxCjAFAGEJsqMhN5BaVbpgwQK8/fbbWLlyJY4ePYrHH38ceXl51qnPzMxM3Hfffdbjp0+fjnXr1mH58uU4efIkvvnmGzz22GOYOHEi+vfvDwCYP38+Nm3ahMWLF+O7777D4sWL8eWXXyIjI0PGQyQtCYgGQi8V14t3uve+S7PFJ/aAvmxH0x0tJ27WadLLAYPUP5/kDgaD9qdL1WnS8DFiLS95PKl/edLT05GVlYWFCxdi3Lhx2LJlCzZu3IhBgwYBAAoKCmz2dJs1axaWLFmCN954A4mJibjzzjsxYsQIm4KI1NRUrFmzBqtWrcKYMWOwevVqrF27FpMmTXL74yMNkrXOjRvv2k9N3Mo0nLhxmtR7aL2ylB0TvI709HzevHmYN29ehz9bvXp1u9seffRRPProo12e84477sAdd9zhjPDI00SnALnvyk3cqGutCxQaL2qrApeJm/fReusrFiZ4HY71k3dpvRFvk8U996korTomcK1ltwL7iYvSpK0OCixM8E7qSFZZjvv+ZthLUWyby5NXYOJG3sU8GvAJBRovABWH3HOfF08DtYWAwQeIZIsku0RocJ1b5dHm3elDgdBhsqMhdwkdLkZ9LTWiQ4GWVOcB9aXib4tZwt6UJAUTN/IuRhMQ3bze0V3Tper9RIwDfALdc596p8UCBZuOCfzT6TWMJvF/F9DeOjc1nvBEwMQtrbwF//qQ93F3gQLXtzlO04kbp0m9jnUEWGuJG1+T3oiJG3kfJm7aZy1QOCoKFLSAb5LeS6sFCqwo9UpM3Mj7RDdvxFt1HKgtdu19NVa3LLBn4ma/oP5AQGxzgcJ+2dGIwgT1eWTi5n2sids+8ZrUAkVp2TKHhQlehYkbeR+/CCBspLju6lG30r2A0iiqJIMHufa+PI2Wpksrv2tVmDBcdjTkbmEJojNBQ6XoVKAFNQVAbZFYbxk+RnY05EZM3Mg7qdtyuDpx48a7PadW4GohcbNOk45nYYI3Mvq0JEdamS5V4whLAHyC5MZCbsW/QOSd3LXOjevbes464rZXbhytY4jgNKnX0loHBa5v81pM3Mg7WTfi3S3WL7mCojBx6w0tFSiwMIG0VqDAjglei4kbeSdzAuBrBizVQPlB19zHxVyxBsXoyzf8ntBKgQILEwiwbTavKHJjUeMAmLh5ISZu5J0MRiDKxRvxnlc33r1cLGwmx2mhQMFamBAChF0qLw6Sy5woOhTUl4qOBTLVngeq88V1dXNg8hpM3Mh79VELFLa75vzqeTlN2nNaSNzU+45gYYJXM/mLDgWA/HVuZfvE19DhgG+Y3FjI7fhXiLyXqwsUuL6t97SUuHGalCI10kGBr0mvxsSNvFfUJAAG4MJJsRbNmRovAuUHxHUmbj1nLVA4IjYzloFvkqSK0EiBAitKvRoTN/JefmbAPEpcd/aoW8keQLEAgZcAwXHOPbc3CewPBMTIK1BosrAwgVpEamRLEFaUejUmbuTdXDVdqp5PXUdHPWMwyJ0urfxOVB77BAOhLEzweuFjxDrH2kKg+pycGOrLxCwBINZdktdh4kbeTe2gcN7JBQrnWZjgNGriViYhcWtdmGA0uf/+SVt8gkSnAkDedKk6Ahw8GPCPlBMDScXEjbybmliV7gWaGpxzTkUBSnbanp96TuaIm3V9W7L775u0SXYHBe7f5vWYuJF3C7tUNJ231DhvDVXVCaCuGDD6cSrDGdTEreII0Fjj3vsuY2ECtSG7gwILE7weEzfybgYjEDVZXHfWOjf1PJFJYu8n6p3AS4CAvqLYo9yNBQpNFqC0eb8sJm6kkl2gwMIEr8fEjcjZBQrcv825DIaW5u7unC6tOsbCBGpP7VRQnS86GLhTwwWg8lhzHEzcvBUTN6I+auLmpAIFa8cEVpQ6jYx1bixMoI74homOBUBLBwN3Kd8PQBHb5ATGuPe+STOYuBFFTRRTphdPAzUFvTtXQxVQcUhc54ib88hM3DhNSm3J6qDA1ySBiRuR+ARtbu5B2Nvp0pLdYrPYoIFAUP/ex0aCtUDhsPsKFPgmSZ2R1UGBhQkEJm5EgrPWuXF9m2sEDQD8+zQXKBxw/f01WVqmwZi4UVvWAgU3b1HDwgQCEzcigYmbtrm7g0LV96LfrE8wEDrC9fdH+qJu83PhpOhk4A6NNWJLHICJm5dj4kYEtCRaJXsBS33PzqE0AcXNG++y1ZXzuTNxsxYmjGNhArXnHyk6FwAtnQxcrfygGHH27yO2yCGvxcSNCBBVYv5RQFNdzyvFKr8H6ksBUwAQPta58ZGkxI3TpNQJd+/n1nqa1GBwz32SJjFxIwLEH8KoXk6XWjfeTQZMfs6Ji1q0LlCw1Lr2vliYQN1xd+srFiZQMyZuRKo+TkrcuL7NNYLiAP9oQGkEylxYoKA0sTCBuufu1lcsTKBmTNyIVL0tUGDi5lqtCxTKXDhdWvk90HgBMAUBYSNddz+kb+rIV+Ux0dHAlSz1Yo0bwMSN5Cduy5YtQ3x8PAICApCUlIStW7d2euysWbNgMBjaXUaPHm09ZvXq1R0eU1vr4qkV0r/ICWIj3up8oPqMY79bXyGm8AAmbq7kjnVuLEwgewTGiA4GUFzfQ7fyCNBUD/iGA8Hxrr0v0jypidvatWuRkZGBJ598Evv27cPUqVNxww03IC8vr8PjX331VRQUFFgv+fn5iIyMxJ133mlzXFhYmM1xBQUFCAgIcMdDIj3zDWkpKnB01K1kFwBF/FENjHV6aNTMnYkbp0mpO+7qoGB9TbIwgSQnbkuWLMHs2bMxZ84cJCQkICsrC3FxcVi+fHmHx5vNZsTGxlove/fuRVlZGe6//36b4wwGg81xsbF8IyU7qaNl5x1M3DhN6h7qG2X5IdcVKJQxcSM7uauDQinXt1ELaYlbfX09srOzkZaWZnN7Wloatm+3r9n3O++8g2uvvRaDBg2yuf3ChQsYNGgQBgwYgJtvvhn79nW9vUNdXR0qKyttLuSlerrOjYmbewQNFNu2KI0ta36cSWkCSlmYQHZyVwcFVpRSK9ISt+LiYlgsFsTExNjcHhMTg8LCwm5/v6CgAJ9//jnmzJljc/vIkSOxevVqbNiwAR9//DECAgIwZcoUHD9+vNNzLVq0CGaz2XqJi4vr2YMi/VMTr7JvAUudfb9js/EuEzeXMhha9lZzxZtl1XGgsQowBbIwgbqnJm4VR1zXQ7epsWUNHUfcCBooTjC0ma9XFKXdbR1ZvXo1wsPDceutt9rcPnnyZNx7770YO3Yspk6dik8++QSXXnopXn/99U7PlZmZiYqKCuslPz+/R4+FPEDIELEzeVO9/etWKr8DGirEm334GNfGR63WFe11/rlLms8ZMQ4w+jj//ORZAi9p1UPXBSPAgKhatdQAPiFio3DyetISt+joaJhMpnaja0VFRe1G4dpSFAUrV67EzJkz4efX9UanRqMREyZM6HLEzd/fH2FhYTYX8lIGg+PTpepxURMBo69r4qIWrixQYGECOcJgcP1+bup5I8aJqnfyetJeBX5+fkhKSsLmzZttbt+8eTNSU7vu8/j111/jxIkTmD17drf3oygKcnJy0K9fv17FS15E7TNabN9aS5xvPo7r29wjKll8dUWBAgsTyFGu7qDA9W3UhtS5gAULFmDmzJlITk5GSkoKVqxYgby8PMydOxeAmMI8e/Ys3nvvPZvfe+eddzBp0iQkJia2O+fzzz+PyZMnY/jw4aisrMRrr72GnJwcvPnmm255TOQBWo+4KUr35fcsTHAvtUChrkRMT0VNcM55WZhAPeGuETeub6NmUhO39PR0lJSUYOHChSgoKEBiYiI2btxorRItKChot6dbRUUFPvvsM7z66qsdnrO8vBwPPfQQCgsLYTabMX78eGzZsgUTJ050+eMhDxGZDBh8gJpzYjPe4IGdH1tfBlQeFdejJ7snPm+nFigUbhJTm85K3GwKExKcc07yfGpCVX5QdDhwZp9imw8TTNxIkL76dt68eZg3b16HP1u9enW728xmM6qrqzs939KlS7F06VJnhUfeyCcIiBgrkoLiHV0nbsW7xNeQoUBAX/fER2JETE3cnEU9V/hYFiaQ/YLjRUeDhnLR4SBinPPOXfVD84eJAH6YICuudCTqiL0FCpwmlcMVBQosTKCeaF2g4Ox1bvwwQR1g4kbUkejmAoXz3RQoqAUMfbouqCEnU5OrikP277fXHSZu1FOuSty4vo06wMSNqCPWjXj3db6xZpOlZaqUI27uFTwI8IsEmhqcs3+W0tSqrRATN3JQhIs6KLCilDrAxI2oI8GDgIBY0Vqpsz/GlUfE+hOfYMDcvsKZXMhgcO50adWJlrVE5lG9Px95F2uBwn7R6cAZFIUjbtQhJm5EHbFnI16bjXe5/sTtnJm4WdcSjeNzSY4LHS46G1hqRKcDZ7h4WlStG30B82jnnJM8AhM3os7Ym7hxmlQOVyRunCalnjAYW6pJnbWfm3oecyJg8nfOOckjMHEj6kzrDgqK0v7n1o4JLEyQwlqgcLD3BQpM3Ki3nN1BoZTTpNQxJm5EnYlMEtMUtT8CF0/Z/qyuBKj6XlznxrtyBA8G/CJEgULFoZ6fR2lqtZaIiRv1kLM7KLAwgTrBxI2oM6YAIGK8uN52urR4p/gaeqlov0Tu56wChaofgIZKFiZQ71i3BNknPgz0hqK06pvLxI1sMXEj6kpn69y4vk0bnJG4cZNTcoawBJH8N1aJDwO9UVMA1BaJtXPhY5wTH3kMJm5EXWHipm3OSNzKuL6NnMDoI5J/oPfTperrOWyUaMFH1AoTN6KuWDfizQEaL4rrTY1AyW5xnR0T5FKTLbXBd0+wMIGcJdJJG/Fy/zbqAhM3oq4ExQGBlwCKBSjZK26rOAQ0XgB8QsUnYpInOL65QKG+ZwUKisKOCeQ8zqosZWECdYGJG1FXOtqI1zpNOgkwmuTERYLB0Lt2Qxd+ABoqAKM/CxOo91pXlna0hZC9OOJGXWDiRtSdtonbea5v05TerHNTfydirNj6hag3zKPF66i+THQ+6InaIqD6jLiubupL1AoTN6LutE7cFIWFCVrjjMSN06TkDCb/lr7FPS1QKN0nvoZeCviGOicu8ihM3Ii6E3k5YPQD6s4DJbuACyfE7dx4VxusBQoHHC9QYOJGzhbZy3VunCalbjBxI+qOyb/ljf27JeJrWIJYFE/yhQwBfMObCxQO2/97LEwgV+htgQILE6gbTNyI7KFOi+Z9Kr4GDwaaLNLCoVYMhp5tw3DhJNBQ3lyYMNoloZEXshYoZPesQIEjbtQNJm5EdlH/qzT/IS74HNgwGMhfJysgaq0n69ysHRPGsDCBnCd8jOh4UFskOiA4or5MfKAAWtrtEbXBxI2oO/nrgO9eaX979Vlg6x1M3rSgN4kbp0nJmXyCWvZ3dLRAQS1MCI4H/COdGxd5DCZuRF1psgDZ82EdabPRfFt2BqdNZWtdoNDUYN/vMHEjV+lpBwVOk5IdmLgRdeX81pY9lTqkANX54jiSJ2Qo4GsGmursK1BQFCZu5Do9LVAoZeJG3WPiRtQVe9eoOLqWhZyrdYGC2pqsK9bCBD8WJpDzte6g4IgyVpRS95i4EXUlsJ9zjyPXcWSdW+vCBJOf62Ii76R2PKg+I4oU7NFQBVR+3/z7LEygzjFxI+pKn6lA0AAAhk4OMIhG9H2mujMq6khEDxI3TpOSK/iGis4HQEvBQXfK9gNQgMBLgMAYl4VG+sfEjagrRhOQ9GrzN22Tt+bvk7LYbF4LHClQYOJGrubodCkLE8hOTNyIuhN3GzD1UyDoEtvbgwaI2+NukxMX2QodCviGdV+goCit3iSZuJGLOFqgwI4JZCcf2QEQ6ULcbcAlM0T1aE2BWNPWZypH2rTEYBSJ2I//FSNq6jqjti7mio1OjX4tDcGJnI0jbuQiHHEjspfRBMRcBQz+hfjKpE177ClQsBYmXMbCBHIdtcDgwknxQaErjTVAxRFxnYkbdYOJGxF5DnsKFLi+jdzBP1J0QACAspyujy0/ACgWIKCvKE4g6gITNyLyHGoyVra/8wIFa+KW7J6YyHvZ20Gh9f5ths4q2IkE6YnbsmXLEB8fj4CAACQlJWHr1s53oJ81axYMBkO7y+jRthtofvbZZxg1ahT8/f0xatQorF+/3tUPg4i0wKZA4Uj7n7NjArlTpJ0FCuyYQA6QmritXbsWGRkZePLJJ7Fv3z5MnToVN9xwA/Ly8jo8/tVXX0VBQYH1kp+fj8jISNx5553WY3bs2IH09HTMnDkT+/fvx8yZM3HXXXdh165d7npYRCSLwdiqmq+DUY6Lp1iYQO4TYWeBAitKyQFSE7clS5Zg9uzZmDNnDhISEpCVlYW4uDgsX768w+PNZjNiY2Otl71796KsrAz333+/9ZisrCxcd911yMzMxMiRI5GZmYlrrrkGWVlZbnpURCRVVwUKLEwgd1ILFCq/F50ROmKpByoOiusccSM7SEvc6uvrkZ2djbS0NJvb09LSsH37drvO8c477+Daa6/FoEGDrLft2LGj3TmnTZvW5Tnr6upQWVlpcyEinbInceM0KblDYExzsYHS3BmhAxWHxXpM33AgeLAbgyO9kpa4FRcXw2KxICbGtrVHTEwMCgsLu/39goICfP7555gzZ47N7YWFhQ6fc9GiRTCbzdZLXFycA4+EiDTF2kFhP9DUaPszJm7kbt3t59Z6/zYWJpAdpBcnGNq8UBVFaXdbR1avXo3w8HDceuutvT5nZmYmKioqrJf8/Hz7gici7QkdBviEApZa2wIFFiaQDN11UGBhAjlIWuIWHR0Nk8nUbiSsqKio3YhZW4qiYOXKlZg5cyb8/GzXqcTGxjp8Tn9/f4SFhdlciEinDMaOt2G4eBqoLwWMvixMIPfpbsSNhQnkIGktr/z8/JCUlITNmzfjZz/7mfX2zZs3Y8aMGV3+7tdff40TJ05g9uzZ7X6WkpKCzZs34/HHH7fetmnTJqSmpjoveIjksbGxERaLxannJfIkJpMJPj4+do2iO1VkElD0tUjchjYXL6lJnPkywOTv3njIe6mJW8UR0SHBJ7DlZ02NYkq/9XFE3ZDaq3TBggWYOXMmkpOTkZKSghUrViAvLw9z584FIKYwz549i/fee8/m99555x1MmjQJiYntPzXPnz8fV155JRYvXowZM2bg73//O7788kts27bNaXHX19ejoKAA1dXVTjsnkacKCgpCv3792o2Ou1RHBQqcJiUZAi8RHRFqi4Dyg0D0xJafVX4HWGoAnxAgdLi8GElXpCZu6enpKCkpwcKFC1FQUIDExERs3LjRWiVaUFDQbk+3iooKfPbZZ3j11Vc7PGdqairWrFmDp556Ck8//TSGDh2KtWvXYtKkSU6JuampCbm5uTCZTOjfvz/8/PzcP5pApAOKoqC+vh7nz59Hbm4uhg8fDqPRTasz2hYoGH2YuJEcBoOYBi34AijLtk3crNOk48UUP5EdpCZuADBv3jzMmzevw5+tXr263W1ms7nbka477rgDd9xxhzPCa6e+vh5NTU2Ii4tDUFCQS+6DyFMEBgbC19cXp0+fRn19PQICAtxzx6HDRYFCYxVQeVSsaStj4kaSRDYnbm0LFMpYmECOY4rfQ24bOSDSOSn/VwxGILJ589PSbKA6D6grEYUJ4Ze5Px7ybp1VlrIwgXqA2QcReaaIVuvcrIUJiSxMIPezFigcFJ0SAEBpAsr22f6cyA5M3IjIM0V2kLhxmpRkCB4sOiM0NYhOCQBQdQJovACYAoCwkTKjI51h4iZTkwX48Svg1Mfia5Nrtxa56qqrkJGR4dL76MqsWbNsNkyWHQ95ODVJK8sBSnbb3kbkTgZD+/3c1GnS8LGieIbITkzcZMlfB2wYDPz7amD73eLrhsHidi+xbt06vPDCC7LDsDp16hQMBkO7yxdffNHl75WVlWHmzJnWlmkzZ85EeXl5l7+jKAqee+459O/fH4GBgbjqqqtw+PBhu2NdsWIFrrrqKoSFhcFgMHR4fz2Jy6OEXSq2WbDUAD/+V9zGxI1kiWyzzo2FCdRDTNxkyF8HbL0DqD5je3v1WXG7lyRvkZGRCA0NlR1GO19++SUKCgqsl5/+9KddHn/33XcjJycHX3zxBb744gvk5ORg5syZXf7Oyy+/jCVLluCNN97Anj17EBsbi+uuuw5VVVV2xVhdXY3rr78e//d//+fUuLrT0NDQq993K4MRCB8nrisWAEYgbJTMiMibtS1QYGEC9RATN2dQFKDxon2X+kpg72MAlI5OJL7snS+Os+d8Skfn6VxjYyN+/etfIzw8HFFRUXjqqaegNJ/jgw8+QHJyMkJDQxEbG4u7774bRUVF1t8tKyvDPffcgz59+iAwMBDDhw/HqlWrrD8/e/Ys0tPTERERgaioKMyYMQOnTp3qNJa2U6WDBw/GSy+9hAceeAChoaEYOHAgVqxYYfM7jt6H6l//+hcCAgLajTg99thj+MlPfmJzW1RUFGJjY62XrjaOPXr0KL744gu8/fbbSElJQUpKCt566y3885//xLFjxzr8HUVRkJWVhSeffBK33XYbEhMT8e6776K6uhofffRRt48FADIyMvC73/0OkydP7lFciqJg2LBh+NOf/mTze4cOHYLRaMQPP/wAQPT9/fOf/4wZM2YgODgYL774YrevA83IXweU57S6oQn4fyO95oMRaUzbvQWtI24cBSbHMHFzBks18EmIfZdPzUDN2S5OpgA1Z8Rx9pzP4lj3hnfffRc+Pj7YtWsXXnvtNSxduhRvv/02ALFH3QsvvID9+/fjb3/7G3JzczFr1izr7z799NM4cuQIPv/8cxw9ehTLly9HdHQ0ADECdPXVVyMkJARbtmzBtm3bEBISguuvvx719fV2x/fKK68gOTkZ+/btw7x58/CrX/0K3333Xa/v49prr0V4eDg+++wz620WiwWffPIJ7rnnHptjb7nlFvTt2xdTpkzBp59+2uV5d+zYAbPZbLPB8+TJk2E2m7F9+/YOfyc3NxeFhYVIS0uz3ubv74+f/OQnnf6Oo7qLy2Aw4IEHHmiXcK1cuRJTp07F0KFDrbc9++yzmDFjBg4ePIgHHnigy9eBZqij2o0XbG/3slFt0pDQYS1T9wVfAPVlzX1zR8uOjHSGKyK9TFxcHJYuXQqDwYARI0bg4MGDWLp0KR588EE88MAD1uOGDBmC1157DRMnTsSFCxcQEhKCvLw8jB8/HsnJyQDECJlqzZo1MBqNePvtt62dJFatWoXw8HB89dVXNklKV2688UbrhsxPPPEEli5diq+++gojR47s1X2YTCakp6fjo48+sva4/fe//42ysjLceeedAICQkBAsWbIEU6ZMgdFoxIYNG5Ceno53330X9957b4fnLSwsRN++fdvd3rdvXxQWFnb6OwAQExNjc3tMTAxOnz7d1T+P3eyJ6/7778czzzyD3bt3Y+LEiWhoaMAHH3yAP/7xjza/c/fdd9u8Nrp6HWhCkwXIno/OR7UNQHYGcMkMwGhyb2zkvQxG0SHh/FbgxFviNvNlgMmNreDIIzBxcwZTEHDXhe6PA4CiLcBXN3Z/3FUbgb5X2nffDpg8ebJNi66UlBS88sorsFgsOHDgAJ577jnk5OSgtLQUTU1NAMQb9ahRo/CrX/0Kt99+O7799lukpaXh1ltvRWpqKgAgOzsbJ06caLdmrba21jrtZo8xY8ZYrxsMBsTGxlqna3t7H/fccw9SUlJw7tw59O/fHx9++CFuvPFGREREAACio6Px+OOPW49PTk5GWVkZXn755U4TNzXOthRF6bYVWtuf2/M7jugurn79+uGmm27CypUrMXHiRPzzn/9EbW2tNZFVqQmaqqvXgSac39p+/agNBajOF8fFXOWuqIhEIcL5rcC5/9fyPZGDOFXqDAYD4BNs3yU2DQgaAKCzN2gDEBQnjrPnfE56o6+trUVaWhpCQkLwwQcfYM+ePVi/fj0AWKchb7jhBpw+fRoZGRk4d+4crrnmGvzmN78BIHq4JiUlIScnx+by/fff4+6777Y7Dl9fX9t/DYPBmkD29j4mTpyIoUOHYs2aNaipqcH69eu7TMgAkegeP36805/Hxsbixx9/bHf7+fPn242otf4dAO1G5IqKijr9HUfZG9ecOXOs/x6rVq1Cenp6u1ZuwcHBNt939TrQhJoC5x5H5CxqIYLSvPWTT4jLt4Eiz8PEzd2MJiDp1eZv2iZdzd8nZblsCmfnzp3tvh8+fDi+++47FBcX4w9/+AOmTp2KkSNH2hQmqPr06YNZs2bhgw8+QFZWlrV44PLLL8fx48fRt29fDBs2zOZiNpudErsz7uPuu+/Ghx9+iH/84x8wGo246aabujx+37596NevX6c/T0lJQUVFBXbv3m29bdeuXaioqOh0FCo+Ph6xsbHYvHmz9bb6+np8/fXXThu5sjeuG2+8EcHBwVi+fDk+//xzmynRrnT2OtCEwM6frx4dR+Qs9aW23x/L8rptoKj3mLjJEHcbMPVTIOgS29uDBojb425z2V3n5+djwYIFOHbsGD7++GO8/vrrmD9/PgYOHAg/Pz+8/vrrOHnyJDZs2NBuj7VnnnkGf//733HixAkcPnwY//znP5GQkABATENGR0djxowZ2Lp1K3Jzc/H1119j/vz5OHOmq2kr+znjPu655x58++23+P3vf4877rjDpun5u+++i48++ghHjx7FsWPH8Kc//QmvvfYaHn30Uesxu3fvxsiRI3H2rCgwSUhIwPXXX48HH3wQO3fuxM6dO/Hggw/i5ptvxogRI6y/N3LkSOsIpsFgQEZGBl566SWsX78ehw4dwqxZsxAUFGT36GRhYSFycnJw4sQJAMDBgwetU9yOxGUymTBr1ixkZmZi2LBhSElJ6fa+u3odaEKfqfaNaveZ6s6oyNvlrwO+XdD+dhbMkIOYuMkSdxtwyyngmv8CqR+Jr7fkujRpA4D77rsPNTU1mDhxIh555BE8+uijeOihh9CnTx+sXr0af/3rXzFq1Cj84Q9/aLdVhJ+fHzIzMzFmzBhceeWVMJlMWLNmDQAgKCgIW7ZswcCBA3HbbbchISEBDzzwAGpqahAWFuaU2J1xH8OHD8eECRNw4MCBdtWkAPDiiy8iOTkZEyZMwJo1a7By5UqbdW/V1dU4duyYzX5mH374IS677DKkpaUhLS0NY8aMwfvvv29z3mPHjqGiosL6/f/+7/8iIyMD8+bNQ3JyMs6ePYtNmzbZva/dn//8Z4wfPx4PPvggAODKK6/E+PHjsWHDBofiAoDZs2ejvr7e7tG2rl4HmiB5VJuonW4LZiAKZjhtSnYwKIqDG4F5gcrKSpjNZlRUVLRLCGpra5Gbm4v4+Hib0Roivfrmm29w1VVX4cyZM05bY9eatP8z+evEm2XrQoWgOJG0ufgDEpGNH78S3XG6c81/WTDjxbrKPVpjVSmRl6qrq0N+fj6efvpp3HXXXS5J2qSKu01s+XF+qyhECOwnpkc50kbuxoIZciJOlZLHCAkJ6fSydetW2eHZ7cMPP+z0cYwe7bzNOj/++GOMGDECFRUVePnll512Xk0xmsQIxuBfiK9M2kgGFsyQE3HEjTxGTk5Opz+75JJLOv2Z1txyyy02HQ9aa7tdSm/MmjXLpjMGEbmIWjBTfRYdr3MziJ+zYIbswMSNPMawYcNkh+AUoaGhdhcpEJEOqAUzW++AKJBpnbyxYIYcw6nSHmJNB5F9+H+FCFK3gSLPwhE3B6lTVdXV1QgMDJQcDZH2VVdXA3DuNC+RLrFghpyAiZuDTCYTwsPDrV0FgoKCnNpfkshTKIqC6upqFBUVITw8HCYT35yIrAUzRD3ExK0H1F6THbWEIiJb4eHh1v8zRETUO0zcesBgMKBfv37o27evzQ76RGTL19eXI21ERE7ExK0XTCYT35SIiIjIbVhVSkRERKQTTNyIiIiIdIKJGxEREZFOcI1bB9QNQysrKyVHQkRERN5AzTm627SciVsHqqqqAABxcXGSIyEiIiJvUlVVBbPZ3OnPDQr70bTT1NSEc+fOITQ0lJvr9kBlZSXi4uKQn5+PsLAw2eGQHfic6ROfN33i86Y/7njOFEVBVVUV+vfvD6Ox85VsHHHrgNFoxIABA2SHoXthYWH8o6QzfM70ic+bPvF50x9XP2ddjbSpWJxAREREpBNM3IiIiIh0gokbOZ2/vz+effZZ+Pv7yw6F7MTnTJ/4vOkTnzf90dJzxuIEIiIiIp3giBsRERGRTjBxIyIiItIJJm5EREREOsHEjYiIiEgnmLiRUyxatAgTJkxAaGgo+vbti1tvvRXHjh2THRY5aNGiRTAYDMjIyJAdCnXh7NmzuPfeexEVFYWgoCCMGzcO2dnZssOiLjQ2NuKpp55CfHw8AgMDMWTIECxcuBBNTU2yQ6NWtmzZgunTp6N///4wGAz429/+ZvNzRVHw3HPPoX///ggMDMRVV12Fw4cPuzVGJm7kFF9//TUeeeQR7Ny5E5s3b0ZjYyPS0tJw8eJF2aGRnfbs2YMVK1ZgzJgxskOhLpSVlWHKlCnw9fXF559/jiNHjuCVV15BeHi47NCoC4sXL8af//xnvPHGGzh69Chefvll/PGPf8Trr78uOzRq5eLFixg7dizeeOONDn/+8ssvY8mSJXjjjTewZ88exMbG4rrrrrP2OHcHbgdCLnH+/Hn07dsXX3/9Na688krZ4VA3Lly4gMsvvxzLli3Diy++iHHjxiErK0t2WNSB3/3ud/jmm2+wdetW2aGQA26++WbExMTgnXfesd52++23IygoCO+//77EyKgzBoMB69evx6233gpAjLb1798fGRkZeOKJJwAAdXV1iImJweLFi/Hwww+7JS6OuJFLVFRUAAAiIyMlR0L2eOSRR3DTTTfh2muvlR0KdWPDhg1ITk7GnXfeib59+2L8+PF46623ZIdF3bjiiivw73//G99//z0AYP/+/di2bRtuvPFGyZGRvXJzc1FYWIi0tDTrbf7+/vjJT36C7du3uy0ONpknp1MUBQsWLMAVV1yBxMRE2eFQN9asWYNvv/0We/bskR0K2eHkyZNYvnw5FixYgP/7v//D7t278dhjj8Hf3x/33Xef7PCoE0888QQqKiowcuRImEwmWCwW/P73v8cvfvEL2aGRnQoLCwEAMTExNrfHxMTg9OnTbouDiRs53a9//WscOHAA27Ztkx0KdSM/Px/z58/Hpk2bEBAQIDscskNTUxOSk5Px0ksvAQDGjx+Pw4cPY/ny5UzcNGzt2rX44IMP8NFHH2H06NHIyclBRkYG+vfvj1/+8peywyMHGAwGm+8VRWl3mysxcSOnevTRR7FhwwZs2bIFAwYMkB0OdSM7OxtFRUVISkqy3maxWLBlyxa88cYbqKurg8lkkhghtdWvXz+MGjXK5raEhAR89tlnkiIie/z2t7/F7373O/z85z8HAFx22WU4ffo0Fi1axMRNJ2JjYwGIkbd+/fpZby8qKmo3CudKXONGTqEoCn79619j3bp1+M9//oP4+HjZIZEdrrnmGhw8eBA5OTnWS3JyMu655x7k5OQwadOgKVOmtNtq5/vvv8egQYMkRUT2qK6uhtFo+5ZrMpm4HYiOxMfHIzY2Fps3b7beVl9fj6+//hqpqalui4MjbuQUjzzyCD766CP8/e9/R2hoqHUtgNlsRmBgoOToqDOhoaHt1iEGBwcjKiqK6xM16vHHH0dqaipeeukl3HXXXdi9ezdWrFiBFStWyA6NujB9+nT8/ve/x8CBAzF69Gjs27cPS5YswQMPPCA7NGrlwoULOHHihPX73Nxc5OTkIDIyEgMHDkRGRgZeeuklDB8+HMOHD8dLL72EoKAg3H333e4LUiFyAgAdXlatWiU7NHLQT37yE2X+/Pmyw6Au/OMf/1ASExMVf39/ZeTIkcqKFStkh0TdqKysVObPn68MHDhQCQgIUIYMGaI8+eSTSl1dnezQqJX//ve/Hb6X/fKXv1QURVGampqUZ599VomNjVX8/f2VK6+8Ujl48KBbY+Q+bkREREQ6wTVuRERERDrBxI2IiIhIJ5i4EREREekEEzciIiIinWDiRkRERKQTTNyIiIiIdIKJGxEREZFOMHEjIiIi0gkmbkREEP12r732WkybNq3dz5YtWwaz2Yy8vDwJkRERtWDiRkQEwGAwYNWqVdi1axf+8pe/WG/Pzc3FE088gVdffRUDBw506n02NDQ49XxE5PmYuBERNYuLi8Orr76K3/zmN8jNzYWiKJg9ezauueYaTJw4ETfeeCNCQkIQExODmTNnori42Pq7X3zxBa644gqEh4cjKioKN998M3744Qfrz0+dOgWDwYBPPvkEV111FQICAvDBBx/g9OnTmD59OiIiIhAcHIzRo0dj48aNMh4+EekAe5USEbVx6623ory8HLfffjteeOEF7NmzB8nJyXjwwQdx3333oaamBk888QQaGxvxn//8BwDw2WefwWAw4LLLLsPFixfxzDPP4NSpU8jJyYHRaMSpU6cQHx+PwYMH45VXXsH48ePh7++Phx56CPX19XjllVcQHByMI0eOICwsDFdeeaXkfwUi0iImbkREbRQVFSExMRElJSX49NNPsW/fPuzatQv/+te/rMecOXMGcXFxOHbsGC699NJ25zh//jz69u2LgwcPIjEx0Zq4ZWVlYf78+dbjxowZg9tvvx3PPvusWx4bEekbp0qJiNro27cvHnroISQkJOBnP/sZsrOz8d///hchISHWy8iRIwHAOh36ww8/4O6778aQIUMQFhaG+Ph4AGhX0JCcnGzz/WOPPYYXX3wRU6ZMwbPPPosDBw644RESkV4xcSMi6oCPjw98fHwAAE1NTZg+fTpycnJsLsePH7dOaU6fPh0lJSV46623sGvXLuzatQsAUF9fb3Pe4OBgm+/nzJmDkydPYubMmTh48CCSk5Px+uuvu+EREpEeMXEjIurG5ZdfjsOHD2Pw4MEYNmyYzSU4OBglJSU4evQonnrqKVxzzTVISEhAWVmZ3eePi4vD3LlzsW7dOvzP//wP3nrrLRc+GiLSMyZuRETdeOSRR1BaWopf/OIX2L17N06ePIlNmzbhgQcegMViQUREBKKiorBixQqcOHEC//nPf7BgwQK7zp2RkYF//etfyM3Nxbfffov//Oc/SEhIcPEjIiK9YuJGRNSN/v3745tvvoHFYsG0adOQmJiI+fPnw2w2w2g0wmg0Ys2aNcjOzkZiYiIef/xx/PGPf7Tr3BaLBY888ggSEhJw/fXXY8SIEVi2bJmLHxER6RWrSomIiIh0giNuRERERDrBxI2IiIhIJ5i4EREREekEEzciIiIinWDiRkRERKQTTNyIiIiIdIKJGxEREZFOMHEjIiIi0gkmbkREREQ6wcSNiIiISCeYuBERERHpBBM3IiIiIp34/+1rdIrrL9Q+AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(1, 1, figsize=(7, 4), sharex=True)\n", + "axs = [axs]\n", + "colors = ['orange', 'blue', 'black', 'red']\n", + "for i, run_name in enumerate(results_allruns.keys()):\n", + " axs[0].plot(years, results_allruns[run_name], label=run_name, marker='o', color=colors[i])\n", + "axs[0].legend()\n", + "\n", + "axs[0].set_xlabel('Years')\n", + "axs[0].set_ylabel('FOM ratio')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demo of meanz metric" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Healpix slicer using NSIDE=32, approximate resolution 109.935565 arcminutes\n", + "year 1\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/yoachim/git_repos/rubin_sim/rubin_sim/maf/maps/dust_map.py:46: UserWarning: Slicer value of nside 32 different from map value 128, using slicer value\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "year 2\n", + "year 3\n", + "year 4\n", + "year 5\n", + "year 6\n", + "year 7\n", + "year 8\n", + "year 9\n" + ] + } + ], + "source": [ + "from rubin_sim.maf.metrics.cosmology_summary_metrics import MultibandMeanzBiasMetric\n", + "from rubin_sim.maf.metrics.uniformity_metrics import MultibandExgalM5\n", + "from rubin_sim.maf.metrics.weak_lensing_systematics_metric import RIZDetectionCoaddExposureTime, ExgalM5WithCuts\n", + "from rubin_sim.maf.metrics.exgal_m5 import ExgalM5\n", + "from rubin_sim.maf.metrics.tomography_models import DENSITY_TOMOGRAPHY_MODEL, MEANZ_TOMOGRAPHY_MODEL\n", + "\n", + "nside = 32\n", + "\n", + "sim_list = [baseline_file]\n", + "name_list = [splitext(basename(sim))[0] for sim in sim_list]\n", + "\n", + "years = range(1, 10) \n", + "surveyAreas = SkyAreaGenerator(nside=nside)\n", + "map_footprints, map_labels = surveyAreas.return_maps()\n", + "slicer = maf.HealpixSubsetSlicer(\n", + " nside=nside, \n", + " hpid=np.where(map_labels == \"lowdust\")[0],#hpid=np.arange(hp.nside2npix(nside)), \n", + " use_cache=False)\n", + "\n", + "results_allruns= {}\n", + "for opsim_fname, run_name in zip(sim_list, name_list):\n", + "\n", + " results_allyears = np.zeros((len(years), ))\n", + " # loop over years \n", + " for iy, year in enumerate(years):\n", + " print('year', year)\n", + "\n", + " days = year*365.25\n", + " constraint_str = 'scheduler_note not like \"DD%\" and night <= XX and scheduler_note not like \"twilight_near_sun\" '\n", + " constraint_str = constraint_str.replace('XX','%d'%days)\n", + "\n", + " \n", + " metric = MultibandExgalM5() \n", + " summary_metrics = [MultibandMeanzBiasMetric(MEANZ_TOMOGRAPHY_MODEL, year=year, \n", + " metric_name='MultibandMeanzBias')]\n", + " depth_map_bundles = [maf.MetricBundle(\n", + " metric=metric,\n", + " slicer=slicer,\n", + " constraint=constraint_str,\n", + " run_name=run_name,\n", + " summary_metrics=summary_metrics\n", + " )]\n", + " bd = maf.metricBundles.make_bundles_dict_from_list(depth_map_bundles)\n", + " bgroup = maf.MetricBundleGroup(bd, opsim_fname)\n", + " bgroup.run_all()\n", + " # take the lowest z-bin I think\n", + " results_allyears[iy] = bd[list(bd.keys())[0]].summary_values['MultibandMeanzBias'][0][\"y1ratio\"]\n", + " \n", + " results_allruns[run_name] = results_allyears\n", + "\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'baseline_v5.0.0_10yrs': array([-0.05167519, -0.05229187, -0.03561607, -0.03225271, -0.02639071,\n", + " -0.02410579, -0.021529 , -0.02272353, -0.0194467 ])}" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "results_allruns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}