From 4fa0e04f6e70d7daedb9a4f85245f658e1ced641 Mon Sep 17 00:00:00 2001 From: Richard Stiskalek Date: Sat, 21 Sep 2024 14:48:23 +0200 Subject: [PATCH] Add more MANTICORE (#148) * Add path to fields * Add new Manticore field support * Downsample Manticore to 20 snapshots * Add manticore snapshot * Remove downsampling * Add Manticore * Add names of manti * Add manti * Simplify * Update script * Update CF4 ICs * Updaet LOS settings * Update nb * Update imports * Reorganise funcs * Update script * Add basic void model --- csiborgtools/flow/__init__.py | 8 +- csiborgtools/flow/flow_model.py | 507 +---------------- csiborgtools/flow/io.py | 518 ++++++++++++++++++ csiborgtools/flow/void_model.py | 130 +++++ csiborgtools/read/__init__.py | 10 +- csiborgtools/read/paths.py | 25 +- csiborgtools/read/snapshot.py | 169 +++++- .../flow/reconstruction_comparison.ipynb | 139 ++++- notebooks/flow/reconstruction_comparison.py | 3 +- scripts/field_prop/field_bulk.py | 8 +- scripts/field_prop/field_los.py | 8 +- scripts/field_prop/field_los.sh | 9 +- scripts/flow/flow_validation.sh | 12 +- 13 files changed, 980 insertions(+), 566 deletions(-) create mode 100644 csiborgtools/flow/io.py create mode 100644 csiborgtools/flow/void_model.py diff --git a/csiborgtools/flow/__init__.py b/csiborgtools/flow/__init__.py index 9061201..2fb6057 100644 --- a/csiborgtools/flow/__init__.py +++ b/csiborgtools/flow/__init__.py @@ -12,9 +12,9 @@ # You should have received a copy of the GNU General Public License along # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -from .flow_model import (DataLoader, PV_LogLikelihood, PV_validation_model, # noqa - dist2redshift, get_model, # noqa +from .io import (DataLoader, get_model, read_absolute_calibration, # noqa + radial_velocity_los) # noqa +from .flow_model import (PV_LogLikelihood, PV_validation_model, dist2redshift, # noqa Observed2CosmologicalRedshift, predict_zobs, # noqa - project_Vext, radial_velocity_los, # noqa - stack_pzosmo_over_realizations) # noqa + project_Vext, stack_pzosmo_over_realizations) # noqa from .selection import ToyMagnitudeSelection # noqa diff --git a/csiborgtools/flow/flow_model.py b/csiborgtools/flow/flow_model.py index 5000176..cbb56dd 100644 --- a/csiborgtools/flow/flow_model.py +++ b/csiborgtools/flow/flow_model.py @@ -26,289 +26,21 @@ from abc import ABC, abstractmethod import numpy as np from astropy import units as u from astropy.cosmology import FlatLambdaCDM, z_at_value -from h5py import File from jax import jit from jax import numpy as jnp -from jax.scipy.special import logsumexp, erf -from numpyro import factor, sample, plate -from numpyro.distributions import Normal, Uniform, MultivariateNormal +from jax.scipy.special import erf, logsumexp +from numpyro import factor, plate, sample +from numpyro.distributions import MultivariateNormal, Normal, Uniform from quadax import simpson from tqdm import trange +from ..params import SPEED_OF_LIGHT +from ..utils import fprint from .selection import toy_log_magnitude_selection -from ..params import SPEED_OF_LIGHT, simname2Omega_m -from ..utils import fprint, radec_to_galactic, radec_to_supergalactic H0 = 100 # km / s / Mpc -############################################################################### -# Data loader # -############################################################################### - - -class DataLoader: - """ - Data loader for the line of sight (LOS) interpolated fields and the - corresponding catalogues. - - Parameters - ---------- - simname : str - Simulation name. - ksim : int or list of int - Index of the simulation to read in (not the IC index). - catalogue : str - Name of the catalogue with LOS objects. - catalogue_fpath : str - Path to the LOS catalogue file. - paths : csiborgtools.read.Paths - Paths object. - ksmooth : int, optional - Smoothing index. - store_full_velocity : bool, optional - Whether to store the full 3D velocity field. Otherwise stores only - the radial velocity. - verbose : bool, optional - Verbose flag. - """ - def __init__(self, simname, ksim, catalogue, catalogue_fpath, paths, - ksmooth=None, store_full_velocity=False, verbose=True): - fprint("reading the catalogue,", verbose) - self._cat = self._read_catalogue(catalogue, catalogue_fpath) - self._catname = catalogue - - fprint("reading the interpolated field.", verbose) - self._field_rdist, self._los_density, self._los_velocity = self._read_field( # noqa - simname, ksim, catalogue, ksmooth, paths) - - if len(self._cat) != self._los_density.shape[1]: - raise ValueError("The number of objects in the catalogue does not " - "match the number of objects in the field.") - - fprint("calculating the radial velocity.", verbose) - nobject = self._los_density.shape[1] - dtype = self._los_density.dtype - - if simname in ["Carrick2015", "Lilow2024"]: - # Carrick+2015 and Lilow+2024 are in galactic coordinates - d1, d2 = radec_to_galactic(self._cat["RA"], self._cat["DEC"]) - elif simname in ["CF4", "CLONES"]: - # CF4 is in supergalactic coordinates - d1, d2 = radec_to_supergalactic(self._cat["RA"], self._cat["DEC"]) - else: - d1, d2 = self._cat["RA"], self._cat["DEC"] - - num_sims = len(self._los_density) - if "IndranilVoid" in simname: - self._los_radial_velocity = self._los_velocity - self._los_velocity = None - else: - radvel = np.empty( - (num_sims, nobject, len(self._field_rdist)), dtype) - for k in range(num_sims): - for i in range(nobject): - radvel[k, i, :] = radial_velocity_los( - self._los_velocity[k, :, i, ...], d1[i], d2[i]) - self._los_radial_velocity = radvel - - if not store_full_velocity: - self._los_velocity = None - - self._Omega_m = simname2Omega_m(simname) - - # Normalize the CSiBORG & CLONES density by the mean matter density - if "csiborg" in simname or simname == "CLONES": - cosmo = FlatLambdaCDM(H0=H0, Om0=self._Omega_m) - mean_rho_matter = cosmo.critical_density0.to("Msun/kpc^3").value - mean_rho_matter *= self._Omega_m - self._los_density /= mean_rho_matter - - # Since Carrick+2015 and CF4 provide `rho / - 1` - if simname in ["Carrick2015", "CF4", "CF4gp"]: - self._los_density += 1 - - # But some CF4 delta values are < -1. Check that CF4 really reports - # this. - if simname in ["CF4", "CF4gp"]: - self._los_density = np.clip(self._los_density, 1e-2, None,) - - # Lilow+2024 outside of the range data is NaN. Replace it with some - # finite values. This is OK because the PV tracers are not so far. - if simname == "Lilow2024": - self._los_density[np.isnan(self._los_density)] = 1. - self._los_radial_velocity[np.isnan(self._los_radial_velocity)] = 0. - - self._mask = np.ones(len(self._cat), dtype=bool) - self._catname = catalogue - - @property - def cat(self): - """The distance indicators catalogue (structured array).""" - return self._cat[self._mask] - - @property - def catname(self): - """Catalogue name.""" - return self._catname - - @property - def rdist(self): - """Radial distances at which the field was interpolated.""" - return self._field_rdist - - @property - def los_density(self): - """ - Density field along the line of sight `(n_sims, n_objects, n_steps)` - """ - return self._los_density[:, self._mask, ...] - - @property - def los_velocity(self): - """ - Velocity field along the line of sight `(n_sims, 3, n_objects, - n_steps)`. - """ - if self._los_velocity is None: - raise ValueError("The 3D velocities were not stored.") - - return self._los_velocity[:, :, self._mask, ...] - - @property - def los_radial_velocity(self): - """ - Radial velocity along the line of sight `(n_sims, n_objects, n_steps)`. - """ - return self._los_radial_velocity[:, self._mask, ...] - - def _read_field(self, simname, ksims, catalogue, ksmooth, paths): - nsims = paths.get_ics(simname) - if isinstance(ksims, int): - ksims = [ksims] - - # For no-field read in Carrick+2015 but then zero it. - if simname == "no_field": - simname = "Carrick2015" - to_wipe = True - else: - to_wipe = False - - if not all(0 <= ksim < len(nsims) for ksim in ksims): - raise ValueError(f"Invalid simulation index: `{ksims}`") - - if "Pantheon+" in catalogue: - fpath = paths.field_los(simname, "Pantheon+") - elif "CF4_TFR" in catalogue: - fpath = paths.field_los(simname, "CF4_TFR") - else: - fpath = paths.field_los(simname, catalogue) - - los_density = [None] * len(ksims) - los_velocity = [None] * len(ksims) - - for n, ksim in enumerate(ksims): - nsim = nsims[ksim] - - with File(fpath, 'r') as f: - has_smoothed = True if f[f"density_{nsim}"].ndim > 2 else False - if has_smoothed and (ksmooth is None or not isinstance(ksmooth, int)): # noqa - raise ValueError("The output contains smoothed field but " - "`ksmooth` is None. It must be provided.") - - indx = (..., ksmooth) if has_smoothed else (...) - los_density[n] = f[f"density_{nsim}"][indx] - los_velocity[n] = f[f"velocity_{nsim}"][indx] - rdist = f[f"rdist_{nsim}"][...] - - los_density = np.stack(los_density) - los_velocity = np.stack(los_velocity) - - if to_wipe: - los_density = np.ones_like(los_density) - los_velocity = np.zeros_like(los_velocity) - - return rdist, los_density, los_velocity - - def _read_catalogue(self, catalogue, catalogue_fpath): - if catalogue == "A2": - with File(catalogue_fpath, 'r') as f: - dtype = [(key, np.float32) for key in f.keys()] - arr = np.empty(len(f["RA"]), dtype=dtype) - for key in f.keys(): - arr[key] = f[key][:] - elif catalogue in ["LOSS", "Foundation", "SFI_gals", "2MTF", - "Pantheon+", "SFI_gals_masked", "SFI_groups", - "Pantheon+_groups", "Pantheon+_groups_zSN", - "Pantheon+_zSN"]: - with File(catalogue_fpath, 'r') as f: - if "Pantheon+" in catalogue: - grp = f["Pantheon+"] - else: - grp = f[catalogue] - - dtype = [(key, np.float32) for key in grp.keys()] - arr = np.empty(len(grp["RA"]), dtype=dtype) - for key in grp.keys(): - arr[key] = grp[key][:] - elif "CB2_" in catalogue: - with File(catalogue_fpath, 'r') as f: - dtype = [(key, np.float32) for key in f.keys()] - arr = np.empty(len(f["RA"]), dtype=dtype) - for key in f.keys(): - arr[key] = f[key][:] - elif "UPGLADE" in catalogue: - with File(catalogue_fpath, 'r') as f: - dtype = [(key, np.float32) for key in f.keys()] - arr = np.empty(len(f["RA"]), dtype=dtype) - for key in f.keys(): - if key == "mask": - continue - - arr[key] = f[key][:] - elif catalogue in ["CF4_GroupAll"] or "CF4_TFR" in catalogue: - with File(catalogue_fpath, 'r') as f: - dtype = [(key, np.float32) for key in f.keys()] - dtype += [("DEC", np.float32)] - arr = np.empty(len(f["RA"]), dtype=dtype) - - for key in f.keys(): - arr[key] = f[key][:] - arr["DEC"] = arr["DE"] - - if "CF4_TFR" in catalogue: - arr["RA"] *= 360 / 24 - else: - raise ValueError(f"Unknown catalogue: `{catalogue}`.") - - return arr - - -############################################################################### -# Supplementary flow functions # -############################################################################### - - -def radial_velocity_los(los_velocity, ra, dec): - """ - Calculate the radial velocity along the LOS from the 3D velocity - along the LOS `(3, n_steps)`. - """ - types = (float, np.float32, np.float64) - if not isinstance(ra, types) and not isinstance(dec, types): - raise ValueError("RA and dec must be floats.") - - if los_velocity.ndim != 2 and los_velocity.shape[0] != 3: - raise ValueError("The shape of `los_velocity` must be (3, n_steps).") - - ra_rad, dec_rad = np.deg2rad(ra), np.deg2rad(dec) - - vx, vy, vz = los_velocity - return (vx * np.cos(ra_rad) * np.cos(dec_rad) - + vy * np.sin(ra_rad) * np.cos(dec_rad) - + vz * np.sin(dec_rad)) - - ############################################################################### # JAX Flow model # ############################################################################### @@ -992,235 +724,6 @@ def PV_validation_model(models, distmod_hyperparams_per_model, factor("ll", ll) -############################################################################### -# Shortcut to create a model # -############################################################################### - - -def read_absolute_calibration(kind, data_length, calibration_fpath): - """ - Read the absolute calibration for the CF4 TFR sample from LEDA but - preprocessed by me. - - Parameters - ---------- - kind : str - Calibration kind: `Cepheids`, `TRGB`, `SBF`, ... - data_length : int - Number of samples in CF4 TFR (should be 9,788). - calibration_fpath : str - Path to the preprocessed calibration file. - - Returns - ------- - data : 3-dimensional array of shape (data_length, max_calib, 2) - Absolute calibration data. - with_calibration : 1-dimensional array of shape (data_length) - Whether the sample has a calibration. - length_calibration : 1-dimensional array of shape (data_length) - Number of calibration points per sample. - """ - data = {} - with File(calibration_fpath, 'r') as f: - for key in f[kind].keys(): - x = f[kind][key][:] - - # Get rid of points without uncertainties - x = x[~np.isnan(x[:, 1])] - - data[key] = x - - max_calib = max(len(val) for val in data.values()) - - out = np.full((data_length, max_calib, 2), np.nan) - with_calibration = np.full(data_length, False) - length_calibration = np.full(data_length, 0) - for i in data.keys(): - out[int(i), :len(data[i]), :] = data[i] - with_calibration[int(i)] = True - length_calibration[int(i)] = len(data[i]) - - return out, with_calibration, length_calibration - - -def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, - absolute_calibration=None, calibration_fpath=None): - """ - Get a model and extract the relevant data from the loader. - - Parameters - ---------- - loader : DataLoader - DataLoader instance. - zcmb_min : float, optional - Minimum observed redshift in the CMB frame to include. - zcmb_max : float, optional - Maximum observed redshift in the CMB frame to include. - mag_selection : dict, optional - Magnitude selection parameters. - add_absolute_calibration : bool, optional - Whether to add an absolute calibration for CF4 TFRs. - calibration_fpath : str, optional - - Returns - ------- - model : NumPyro model - """ - zcmb_min = 0.0 if zcmb_min is None else zcmb_min - zcmb_max = np.infty if zcmb_max is None else zcmb_max - - los_overdensity = loader.los_density - los_velocity = loader.los_radial_velocity - kind = loader._catname - - if absolute_calibration is not None and "CF4_TFR_" not in kind: - raise ValueError("Absolute calibration supported only for the CF4 TFR sample.") # noqa - - if kind in ["LOSS", "Foundation"]: - keys = ["RA", "DEC", "z_CMB", "mB", "x1", "c", "e_mB", "e_x1", "e_c"] - RA, dec, zCMB, mag, x1, c, e_mag, e_x1, e_c = ( - loader.cat[k] for k in keys) - e_zCMB = None - - mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) - calibration_params = {"mag": mag[mask], "x1": x1[mask], "c": c[mask], - "e_mag": e_mag[mask], "e_x1": e_x1[mask], - "e_c": e_c[mask]} - - model = PV_LogLikelihood( - los_overdensity[:, mask], los_velocity[:, mask], - RA[mask], dec[mask], zCMB[mask], e_zCMB, calibration_params, - None, mag_selection, loader.rdist, loader._Omega_m, "SN", - name=kind) - elif "Pantheon+" in kind: - keys = ["RA", "DEC", "zCMB", "mB", "x1", "c", "biasCor_m_b", "mBERR", - "x1ERR", "cERR", "biasCorErr_m_b", "zCMB_SN", "zCMB_Group", - "zCMBERR"] - - RA, dec, zCMB, mB, x1, c, bias_corr_mB, e_mB, e_x1, e_c, e_bias_corr_mB, zCMB_SN, zCMB_Group, e_zCMB = (loader.cat[k] for k in keys) # noqa - mB -= bias_corr_mB - e_mB = np.sqrt(e_mB**2 + e_bias_corr_mB**2) - - mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) - - if kind == "Pantheon+_groups": - mask &= np.isfinite(zCMB_Group) - - if kind == "Pantheon+_groups_zSN": - mask &= np.isfinite(zCMB_Group) - zCMB = zCMB_SN - - if kind == "Pantheon+_zSN": - zCMB = zCMB_SN - - calibration_params = {"mag": mB[mask], "x1": x1[mask], "c": c[mask], - "e_mag": e_mB[mask], "e_x1": e_x1[mask], - "e_c": e_c[mask]} - model = PV_LogLikelihood( - los_overdensity[:, mask], los_velocity[:, mask], - RA[mask], dec[mask], zCMB[mask], e_zCMB[mask], calibration_params, - None, mag_selection, loader.rdist, loader._Omega_m, "SN", - name=kind) - elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"]: - keys = ["RA", "DEC", "z_CMB", "mag", "eta", "e_mag", "e_eta"] - RA, dec, zCMB, mag, eta, e_mag, e_eta = (loader.cat[k] for k in keys) - - mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) - calibration_params = {"mag": mag[mask], "eta": eta[mask], - "e_mag": e_mag[mask], "e_eta": e_eta[mask]} - model = PV_LogLikelihood( - los_overdensity[:, mask], los_velocity[:, mask], - RA[mask], dec[mask], zCMB[mask], None, calibration_params, None, - mag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind) - elif "CF4_TFR_" in kind: - # The full name can be e.g. "CF4_TFR_not2MTForSFI_i" or "CF4_TFR_i". - band = kind.split("_")[-1] - if band not in ['g', 'r', 'i', 'z', 'w1', 'w2']: - raise ValueError(f"Band `{band}` not recognized.") - - keys = ["RA", "DEC", "Vcmb", f"{band}", "lgWmxi", "elgWi", - "not_matched_to_2MTF_or_SFI", "Qs", "Qw"] - RA, dec, z_obs, mag, eta, e_eta, not_matched_to_2MTF_or_SFI, Qs, Qw = ( - loader.cat[k] for k in keys) - l, b = radec_to_galactic(RA, dec) - - not_matched_to_2MTF_or_SFI = not_matched_to_2MTF_or_SFI.astype(bool) - # NOTE: fiducial uncertainty until we can get the actual values. - e_mag = 0.05 * np.ones_like(mag) - - z_obs /= SPEED_OF_LIGHT - eta -= 2.5 - - fprint("selecting only galaxies with mag > 5 and eta > -0.3.") - mask = (mag > 5) & (eta > -0.3) - fprint("selecting only galaxies with |b| > 7.5.") - mask &= np.abs(b) > 7.5 - mask &= (z_obs < zcmb_max) & (z_obs > zcmb_min) - - if "not2MTForSFI" in kind: - mask &= not_matched_to_2MTF_or_SFI - elif "2MTForSFI" in kind: - mask &= ~not_matched_to_2MTF_or_SFI - - fprint("employing a quality cut on the galaxies.") - if "w" in band: - mask &= Qw == 5 - else: - mask &= Qs == 5 - - # Read the absolute calibration - if absolute_calibration is not None: - CF4_length = len(RA) - distmod, with_calibration, length_calibration = read_absolute_calibration( # noqa - "Cepheids", CF4_length, calibration_fpath) - - distmod = distmod[mask] - with_calibration = with_calibration[mask] - length_calibration = length_calibration[mask] - fprint(f"found {np.sum(with_calibration)} galaxies with absolute calibration.") # noqa - - distmod = distmod[with_calibration] - length_calibration = length_calibration[with_calibration] - - abs_calibration_params = { - "calibration_distmod": distmod, - "data_with_calibration": with_calibration, - "length_calibration": length_calibration} - else: - abs_calibration_params = None - - calibration_params = {"mag": mag[mask], "eta": eta[mask], - "e_mag": e_mag[mask], "e_eta": e_eta[mask]} - model = PV_LogLikelihood( - los_overdensity[:, mask], los_velocity[:, mask], - RA[mask], dec[mask], z_obs[mask], None, calibration_params, - abs_calibration_params, mag_selection, loader.rdist, - loader._Omega_m, "TFR", name=kind) - elif kind in ["CF4_GroupAll"]: - # Note, this for some reason works terribly. - keys = ["RA", "DE", "Vcmb", "DMzp", "eDM"] - RA, dec, zCMB, mu, e_mu = (loader.cat[k] for k in keys) - - zCMB /= SPEED_OF_LIGHT - mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) & np.isfinite(mu) - - # The distance moduli in CF4 are most likely given assuming h = 0.75 - mu += 5 * np.log10(0.75) - - calibration_params = {"mu": mu[mask], "e_mu": e_mu[mask]} - model = PV_LogLikelihood( - los_overdensity[:, mask], los_velocity[:, mask], - RA[mask], dec[mask], zCMB[mask], None, calibration_params, None, - mag_selection, loader.rdist, loader._Omega_m, "simple", - name=kind) - else: - raise ValueError(f"Catalogue `{kind}` not recognized.") - - fprint(f"selected {np.sum(mask)}/{len(mask)} galaxies in catalogue `{kind}`") # noqa - - return model - - ############################################################################### # Predicting z_cosmo from z_obs # ############################################################################### diff --git a/csiborgtools/flow/io.py b/csiborgtools/flow/io.py new file mode 100644 index 0000000..1a941b7 --- /dev/null +++ b/csiborgtools/flow/io.py @@ -0,0 +1,518 @@ +# Copyright (C) 2024 Richard Stiskalek +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +import numpy as np +from astropy.cosmology import FlatLambdaCDM +from h5py import File + +from ..params import SPEED_OF_LIGHT, simname2Omega_m +from ..utils import fprint, radec_to_galactic, radec_to_supergalactic +from .flow_model import PV_LogLikelihood + +H0 = 100 # km / s / Mpc + + +############################################################################## +# Data loader # +############################################################################### + + +class DataLoader: + """ + Data loader for the line of sight (LOS) interpolated fields and the + corresponding catalogues. + + Parameters + ---------- + simname : str + Simulation name. + ksim : int or list of int + Index of the simulation to read in (not the IC index). + catalogue : str + Name of the catalogue with LOS objects. + catalogue_fpath : str + Path to the LOS catalogue file. + paths : csiborgtools.read.Paths + Paths object. + ksmooth : int, optional + Smoothing index. + store_full_velocity : bool, optional + Whether to store the full 3D velocity field. Otherwise stores only + the radial velocity. + verbose : bool, optional + Verbose flag. + """ + def __init__(self, simname, ksim, catalogue, catalogue_fpath, paths, + ksmooth=None, store_full_velocity=False, verbose=True): + fprint("reading the catalogue,", verbose) + self._cat = self._read_catalogue(catalogue, catalogue_fpath) + self._catname = catalogue + + fprint("reading the interpolated field.", verbose) + self._field_rdist, self._los_density, self._los_velocity = self._read_field( # noqa + simname, ksim, catalogue, ksmooth, paths) + + if len(self._cat) != self._los_density.shape[1]: + raise ValueError("The number of objects in the catalogue does not " + "match the number of objects in the field.") + + fprint("calculating the radial velocity.", verbose) + nobject = self._los_density.shape[1] + dtype = self._los_density.dtype + + if simname in ["Carrick2015", "Lilow2024"]: + # Carrick+2015 and Lilow+2024 are in galactic coordinates + d1, d2 = radec_to_galactic(self._cat["RA"], self._cat["DEC"]) + elif simname in ["CF4", "CLONES"]: + # CF4 is in supergalactic coordinates + d1, d2 = radec_to_supergalactic(self._cat["RA"], self._cat["DEC"]) + else: + d1, d2 = self._cat["RA"], self._cat["DEC"] + + num_sims = len(self._los_density) + if "IndranilVoid" in simname: + self._los_radial_velocity = self._los_velocity + self._los_velocity = None + else: + radvel = np.empty( + (num_sims, nobject, len(self._field_rdist)), dtype) + for k in range(num_sims): + for i in range(nobject): + radvel[k, i, :] = radial_velocity_los( + self._los_velocity[k, :, i, ...], d1[i], d2[i]) + self._los_radial_velocity = radvel + + if not store_full_velocity: + self._los_velocity = None + + self._Omega_m = simname2Omega_m(simname) + + # Normalize the CSiBORG & CLONES density by the mean matter density + if "csiborg" in simname or simname == "CLONES": + cosmo = FlatLambdaCDM(H0=H0, Om0=self._Omega_m) + mean_rho_matter = cosmo.critical_density0.to("Msun/kpc^3").value + mean_rho_matter *= self._Omega_m + self._los_density /= mean_rho_matter + + # Since Carrick+2015 and CF4 provide `rho / - 1` + if simname in ["Carrick2015", "CF4", "CF4gp"]: + self._los_density += 1 + + # But some CF4 delta values are < -1. Check that CF4 really reports + # this. + if simname in ["CF4", "CF4gp"]: + self._los_density = np.clip(self._los_density, 1e-2, None,) + + # Lilow+2024 outside of the range data is NaN. Replace it with some + # finite values. This is OK because the PV tracers are not so far. + if simname == "Lilow2024": + self._los_density[np.isnan(self._los_density)] = 1. + self._los_radial_velocity[np.isnan(self._los_radial_velocity)] = 0. + + self._mask = np.ones(len(self._cat), dtype=bool) + self._catname = catalogue + + @property + def cat(self): + """The distance indicators catalogue (structured array).""" + return self._cat[self._mask] + + @property + def catname(self): + """Catalogue name.""" + return self._catname + + @property + def rdist(self): + """Radial distances at which the field was interpolated.""" + return self._field_rdist + + @property + def los_density(self): + """ + Density field along the line of sight `(n_sims, n_objects, n_steps)` + """ + return self._los_density[:, self._mask, ...] + + @property + def los_velocity(self): + """ + Velocity field along the line of sight `(n_sims, 3, n_objects, + n_steps)`. + """ + if self._los_velocity is None: + raise ValueError("The 3D velocities were not stored.") + + return self._los_velocity[:, :, self._mask, ...] + + @property + def los_radial_velocity(self): + """ + Radial velocity along the line of sight `(n_sims, n_objects, n_steps)`. + """ + return self._los_radial_velocity[:, self._mask, ...] + + def _read_field(self, simname, ksims, catalogue, ksmooth, paths): + nsims = paths.get_ics(simname) + if isinstance(ksims, int): + ksims = [ksims] + + # For no-field read in Carrick+2015 but then zero it. + if simname == "no_field": + simname = "Carrick2015" + to_wipe = simname == "no_field" + + if not all(0 <= ksim < len(nsims) for ksim in ksims): + raise ValueError(f"Invalid simulation index: `{ksims}`") + + if "Pantheon+" in catalogue: + fpath = paths.field_los(simname, "Pantheon+") + elif "CF4_TFR" in catalogue: + fpath = paths.field_los(simname, "CF4_TFR") + else: + fpath = paths.field_los(simname, catalogue) + + los_density = [None] * len(ksims) + los_velocity = [None] * len(ksims) + + for n, ksim in enumerate(ksims): + nsim = nsims[ksim] + + with File(fpath, 'r') as f: + has_smoothed = True if f[f"density_{nsim}"].ndim > 2 else False + if has_smoothed and (ksmooth is None or not isinstance(ksmooth, int)): # noqa + raise ValueError("The output contains smoothed field but " + "`ksmooth` is None. It must be provided.") + + indx = (..., ksmooth) if has_smoothed else (...) + los_density[n] = f[f"density_{nsim}"][indx] + los_velocity[n] = f[f"velocity_{nsim}"][indx] + rdist = f[f"rdist_{nsim}"][...] + + los_density = np.stack(los_density) + los_velocity = np.stack(los_velocity) + + if to_wipe: + los_density = np.ones_like(los_density) + los_velocity = np.zeros_like(los_velocity) + + return rdist, los_density, los_velocity + + def _read_catalogue(self, catalogue, catalogue_fpath): + if catalogue == "A2": + with File(catalogue_fpath, 'r') as f: + dtype = [(key, np.float32) for key in f.keys()] + arr = np.empty(len(f["RA"]), dtype=dtype) + for key in f.keys(): + arr[key] = f[key][:] + elif catalogue in ["LOSS", "Foundation", "SFI_gals", "2MTF", + "Pantheon+", "SFI_gals_masked", "SFI_groups", + "Pantheon+_groups", "Pantheon+_groups_zSN", + "Pantheon+_zSN"]: + with File(catalogue_fpath, 'r') as f: + if "Pantheon+" in catalogue: + grp = f["Pantheon+"] + else: + grp = f[catalogue] + + dtype = [(key, np.float32) for key in grp.keys()] + arr = np.empty(len(grp["RA"]), dtype=dtype) + for key in grp.keys(): + arr[key] = grp[key][:] + elif "CB2_" in catalogue: + with File(catalogue_fpath, 'r') as f: + dtype = [(key, np.float32) for key in f.keys()] + arr = np.empty(len(f["RA"]), dtype=dtype) + for key in f.keys(): + arr[key] = f[key][:] + elif "UPGLADE" in catalogue: + with File(catalogue_fpath, 'r') as f: + dtype = [(key, np.float32) for key in f.keys()] + arr = np.empty(len(f["RA"]), dtype=dtype) + for key in f.keys(): + if key == "mask": + continue + + arr[key] = f[key][:] + elif catalogue in ["CF4_GroupAll"] or "CF4_TFR" in catalogue: + with File(catalogue_fpath, 'r') as f: + dtype = [(key, np.float32) for key in f.keys()] + dtype += [("DEC", np.float32)] + arr = np.empty(len(f["RA"]), dtype=dtype) + + for key in f.keys(): + arr[key] = f[key][:] + arr["DEC"] = arr["DE"] + + if "CF4_TFR" in catalogue: + arr["RA"] *= 360 / 24 + else: + raise ValueError(f"Unknown catalogue: `{catalogue}`.") + + return arr + + +############################################################################### +# Supplementary flow functions # +############################################################################### + + +def radial_velocity_los(los_velocity, ra, dec): + """ + Calculate the radial velocity along the LOS from the 3D velocity + along the LOS `(3, n_steps)`. + """ + types = (float, np.float32, np.float64) + if not isinstance(ra, types) and not isinstance(dec, types): + raise ValueError("RA and dec must be floats.") + + if los_velocity.ndim != 2 and los_velocity.shape[0] != 3: + raise ValueError("The shape of `los_velocity` must be (3, n_steps).") + + ra_rad, dec_rad = np.deg2rad(ra), np.deg2rad(dec) + + vx, vy, vz = los_velocity + return (vx * np.cos(ra_rad) * np.cos(dec_rad) + + vy * np.sin(ra_rad) * np.cos(dec_rad) + + vz * np.sin(dec_rad)) + + +############################################################################## +# Shortcut to create a model # +############################################################################### + + +def read_absolute_calibration(kind, data_length, calibration_fpath): + """ + Read the absolute calibration for the CF4 TFR sample from LEDA but + preprocessed by me. + + Parameters + ---------- + kind : str + Calibration kind: `Cepheids`, `TRGB`, `SBF`, ... + data_length : int + Number of samples in CF4 TFR (should be 9,788). + calibration_fpath : str + Path to the preprocessed calibration file. + + Returns + ------- + data : 3-dimensional array of shape (data_length, max_calib, 2) + Absolute calibration data. + with_calibration : 1-dimensional array of shape (data_length) + Whether the sample has a calibration. + length_calibration : 1-dimensional array of shape (data_length) + Number of calibration points per sample. + """ + data = {} + with File(calibration_fpath, 'r') as f: + for key in f[kind].keys(): + x = f[kind][key][:] + + # Get rid of points without uncertainties + x = x[~np.isnan(x[:, 1])] + + data[key] = x + + max_calib = max(len(val) for val in data.values()) + + out = np.full((data_length, max_calib, 2), np.nan) + with_calibration = np.full(data_length, False) + length_calibration = np.full(data_length, 0) + for i in data.keys(): + out[int(i), :len(data[i]), :] = data[i] + with_calibration[int(i)] = True + length_calibration[int(i)] = len(data[i]) + + return out, with_calibration, length_calibration + + +def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, + absolute_calibration=None, calibration_fpath=None): + """ + Get a model and extract the relevant data from the loader. + + Parameters + ---------- + loader : DataLoader + DataLoader instance. + zcmb_min : float, optional + Minimum observed redshift in the CMB frame to include. + zcmb_max : float, optional + Maximum observed redshift in the CMB frame to include. + mag_selection : dict, optional + Magnitude selection parameters. + add_absolute_calibration : bool, optional + Whether to add an absolute calibration for CF4 TFRs. + calibration_fpath : str, optional + + Returns + ------- + model : NumPyro model + """ + zcmb_min = 0.0 if zcmb_min is None else zcmb_min + zcmb_max = np.infty if zcmb_max is None else zcmb_max + + los_overdensity = loader.los_density + los_velocity = loader.los_radial_velocity + kind = loader._catname + + if absolute_calibration is not None and "CF4_TFR_" not in kind: + raise ValueError("Absolute calibration supported only for the CF4 TFR sample.") # noqa + + if kind in ["LOSS", "Foundation"]: + keys = ["RA", "DEC", "z_CMB", "mB", "x1", "c", "e_mB", "e_x1", "e_c"] + RA, dec, zCMB, mag, x1, c, e_mag, e_x1, e_c = ( + loader.cat[k] for k in keys) + e_zCMB = None + + mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) + calibration_params = {"mag": mag[mask], "x1": x1[mask], "c": c[mask], + "e_mag": e_mag[mask], "e_x1": e_x1[mask], + "e_c": e_c[mask]} + + model = PV_LogLikelihood( + los_overdensity[:, mask], los_velocity[:, mask], + RA[mask], dec[mask], zCMB[mask], e_zCMB, calibration_params, + None, mag_selection, loader.rdist, loader._Omega_m, "SN", + name=kind) + elif "Pantheon+" in kind: + keys = ["RA", "DEC", "zCMB", "mB", "x1", "c", "biasCor_m_b", "mBERR", + "x1ERR", "cERR", "biasCorErr_m_b", "zCMB_SN", "zCMB_Group", + "zCMBERR"] + + RA, dec, zCMB, mB, x1, c, bias_corr_mB, e_mB, e_x1, e_c, e_bias_corr_mB, zCMB_SN, zCMB_Group, e_zCMB = (loader.cat[k] for k in keys) # noqa + mB -= bias_corr_mB + e_mB = np.sqrt(e_mB**2 + e_bias_corr_mB**2) + + mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) + + if kind == "Pantheon+_groups": + mask &= np.isfinite(zCMB_Group) + + if kind == "Pantheon+_groups_zSN": + mask &= np.isfinite(zCMB_Group) + zCMB = zCMB_SN + + if kind == "Pantheon+_zSN": + zCMB = zCMB_SN + + calibration_params = {"mag": mB[mask], "x1": x1[mask], "c": c[mask], + "e_mag": e_mB[mask], "e_x1": e_x1[mask], + "e_c": e_c[mask]} + model = PV_LogLikelihood( + los_overdensity[:, mask], los_velocity[:, mask], + RA[mask], dec[mask], zCMB[mask], e_zCMB[mask], calibration_params, + None, mag_selection, loader.rdist, loader._Omega_m, "SN", + name=kind) + elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"]: + keys = ["RA", "DEC", "z_CMB", "mag", "eta", "e_mag", "e_eta"] + RA, dec, zCMB, mag, eta, e_mag, e_eta = (loader.cat[k] for k in keys) + + mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) + calibration_params = {"mag": mag[mask], "eta": eta[mask], + "e_mag": e_mag[mask], "e_eta": e_eta[mask]} + model = PV_LogLikelihood( + los_overdensity[:, mask], los_velocity[:, mask], + RA[mask], dec[mask], zCMB[mask], None, calibration_params, None, + mag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind) + elif "CF4_TFR_" in kind: + # The full name can be e.g. "CF4_TFR_not2MTForSFI_i" or "CF4_TFR_i". + band = kind.split("_")[-1] + if band not in ['g', 'r', 'i', 'z', 'w1', 'w2']: + raise ValueError(f"Band `{band}` not recognized.") + + keys = ["RA", "DEC", "Vcmb", f"{band}", "lgWmxi", "elgWi", + "not_matched_to_2MTF_or_SFI", "Qs", "Qw"] + RA, dec, z_obs, mag, eta, e_eta, not_matched_to_2MTF_or_SFI, Qs, Qw = ( + loader.cat[k] for k in keys) + l, b = radec_to_galactic(RA, dec) + + not_matched_to_2MTF_or_SFI = not_matched_to_2MTF_or_SFI.astype(bool) + # NOTE: fiducial uncertainty until we can get the actual values. + e_mag = 0.05 * np.ones_like(mag) + + z_obs /= SPEED_OF_LIGHT + eta -= 2.5 + + fprint("selecting only galaxies with mag > 5 and eta > -0.3.") + mask = (mag > 5) & (eta > -0.3) + fprint("selecting only galaxies with |b| > 7.5.") + mask &= np.abs(b) > 7.5 + mask &= (z_obs < zcmb_max) & (z_obs > zcmb_min) + + if "not2MTForSFI" in kind: + mask &= not_matched_to_2MTF_or_SFI + elif "2MTForSFI" in kind: + mask &= ~not_matched_to_2MTF_or_SFI + + fprint("employing a quality cut on the galaxies.") + if "w" in band: + mask &= Qw == 5 + else: + mask &= Qs == 5 + + # Read the absolute calibration + if absolute_calibration is not None: + CF4_length = len(RA) + distmod, with_calibration, length_calibration = read_absolute_calibration( # noqa + "Cepheids", CF4_length, calibration_fpath) + + distmod = distmod[mask] + with_calibration = with_calibration[mask] + length_calibration = length_calibration[mask] + fprint(f"found {np.sum(with_calibration)} galaxies with absolute calibration.") # noqa + + distmod = distmod[with_calibration] + length_calibration = length_calibration[with_calibration] + + abs_calibration_params = { + "calibration_distmod": distmod, + "data_with_calibration": with_calibration, + "length_calibration": length_calibration} + else: + abs_calibration_params = None + + calibration_params = {"mag": mag[mask], "eta": eta[mask], + "e_mag": e_mag[mask], "e_eta": e_eta[mask]} + model = PV_LogLikelihood( + los_overdensity[:, mask], los_velocity[:, mask], + RA[mask], dec[mask], z_obs[mask], None, calibration_params, + abs_calibration_params, mag_selection, loader.rdist, + loader._Omega_m, "TFR", name=kind) + elif kind in ["CF4_GroupAll"]: + # Note, this for some reason works terribly. + keys = ["RA", "DE", "Vcmb", "DMzp", "eDM"] + RA, dec, zCMB, mu, e_mu = (loader.cat[k] for k in keys) + + zCMB /= SPEED_OF_LIGHT + mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) & np.isfinite(mu) + + # The distance moduli in CF4 are most likely given assuming h = 0.75 + mu += 5 * np.log10(0.75) + + calibration_params = {"mu": mu[mask], "e_mu": e_mu[mask]} + model = PV_LogLikelihood( + los_overdensity[:, mask], los_velocity[:, mask], + RA[mask], dec[mask], zCMB[mask], None, calibration_params, None, + mag_selection, loader.rdist, loader._Omega_m, "simple", + name=kind) + else: + raise ValueError(f"Catalogue `{kind}` not recognized.") + + fprint(f"selected {np.sum(mask)}/{len(mask)} galaxies in catalogue `{kind}`") # noqa + + return model diff --git a/csiborgtools/flow/void_model.py b/csiborgtools/flow/void_model.py new file mode 100644 index 0000000..b15c380 --- /dev/null +++ b/csiborgtools/flow/void_model.py @@ -0,0 +1,130 @@ +# Copyright (C) 2024 Richard Stiskalek +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +"""Function to work with the void data from Sergij & Indranil's files.""" + +from glob import glob +from os.path import join +from re import search + +import numpy as np +from jax import numpy as jnp +from jax import vmap +from jax.scipy.ndimage import map_coordinates +from tqdm import tqdm + +############################################################################### +# I/O of the void data # +############################################################################### + + +def load_void_data(kind): + """ + Load the void velocities from Sergij & Indranil's files for a given kind + of void profile per observer. + + Parameters + ---------- + kind : str + The kind of void profile to load. One of "exp", "gauss", "mb". + + Returns + ------- + velocities : 3-dimensional array of shape (nLG, nrad, nphi) + """ + if kind not in ["exp", "gauss", "mb"]: + raise ValueError("kind must be one of 'exp', 'gauss', 'mb'") + + fdir = "/mnt/extraspace/rstiskalek/catalogs/IndranilVoid" + + kind = kind.upper() + fdir = join(fdir, f"{kind}profile") + + files = glob(join(fdir, "*.dat")) + rLG = [int(search(rf'v_pec_{kind}profile_rLG_(\d+)', f).group(1)) + for f in files] + rLG = np.sort(rLG) + + for i, ri in enumerate(tqdm(rLG, desc="Loading observer data")): + f = join(fdir, f"v_pec_{kind}profile_rLG_{ri}.dat") + data_i = np.genfromtxt(f).T + + if i == 0: + data = np.full((len(rLG), *data_i.shape), np.nan, dtype=np.float32) + + data[i] = data_i + + if np.any(np.isnan(data)): + raise ValueError("Found NaNs in loaded data.") + + return rLG, data + +############################################################################### +# Interpolation of void velocities # +############################################################################### + + +def interpolate_void(rLG, r, phi, data, rgrid_min, rgrid_max, rLG_min, rLG_max, + order=1): + """ + Interpolate the void velocities from Sergij & Indranil's files for a given + observer over a set of radial distances and at angles specifying the + galaxies. + + Parameters + ---------- + rLG : float + The observer's distance from the center of the void. + r : 1-dimensional array + The radial distances at which to interpolate the velocities. + phi : 1-dimensional array + The angles at which to interpolate the velocities, in degrees, + defining the galaxy position. + data : 3-dimensional array of shape (nLG, nrad, nphi) + The void velocities for different observers, radial distances, and + angles. + rgrid_min, rgrid_max : float + The minimum and maximum radial distances in the data. + rLG_min, rLG_max : float + The minimum and maximum observer distances in the data. + order : int, optional + The order of the interpolation. Default is 1, can be 0. + + Returns + ------- + vel : 2-dimensional array of shape (len(phi), len(r)) + """ + nLG, nrad, nphi = data.shape + + # Normalize rLG to the grid scale + rLG_normalized = (rLG - rLG_min) / (rLG_max - rLG_min) * (nLG - 1) + rLG_normalized = jnp.repeat(rLG_normalized, r.size) + r_normalized = (r - rgrid_min) / (rgrid_max - rgrid_min) * (nrad - 1) + + # Function to perform interpolation for a single phi + def interpolate_single_phi(phi_val): + # Normalize phi to match the grid + phi_normalized = phi_val / 180 * (nphi - 1) + + # Create the grid for this specific phi + X = jnp.vstack([rLG_normalized, + r_normalized, + jnp.repeat(phi_normalized, r.size)]) + + # Interpolate over the data using map_coordinates. The mode is nearest + # to avoid extrapolation. But values outside of the grid should never + # occur. + return map_coordinates(data, X, order=order, mode='nearest') + + return vmap(interpolate_single_phi)(phi) diff --git a/csiborgtools/read/__init__.py b/csiborgtools/read/__init__.py index ea35f7e..9ae0ed1 100644 --- a/csiborgtools/read/__init__.py +++ b/csiborgtools/read/__init__.py @@ -16,11 +16,11 @@ from .catalogue import (CSiBORG1Catalogue, CSiBORG2Catalogue, CSiBORG2SUBFINDCatalogue, # noqa CSiBORG2MergerTreeReader, QuijoteCatalogue, # noqa MDPL2Catalogue, fiducial_observers) # noqa -from .snapshot import (CSiBORG1Snapshot, CSiBORG2Snapshot, QuijoteSnapshot, # noqa - CSiBORG1Field, CSiBORG2Field, CSiBORG2XField, # noqa - QuijoteField, BORG2Field, BORG1Field, TNG300_1Field, # noqa - Carrick2015Field, Lilow2024Field, CLONESField, # noqa - CF4Field) # noqa +from .snapshot import (CSiBORG1Snapshot, CSiBORG2Snapshot, CSiBORG2XSnapshot, # noqa + QuijoteSnapshot, CSiBORG1Field, CSiBORG2Field, # noqa + CSiBORG2XField, QuijoteField, BORG2Field, BORG1Field, # noqa + TNG300_1Field, Carrick2015Field, Lilow2024Field, # noqa + CLONESField, CF4Field) # noqa from .obs import (SDSS, MCXCClusters, PlanckClusters, TwoMPPGalaxies, # noqa TwoMPPGroups, ObservedCluster, match_array_to_no_masking, # noqa cols_to_structured, read_pantheonplus_data) # noqa diff --git a/csiborgtools/read/paths.py b/csiborgtools/read/paths.py index c6670fb..5bcc902 100644 --- a/csiborgtools/read/paths.py +++ b/csiborgtools/read/paths.py @@ -118,6 +118,8 @@ class Paths: fdir = join(fdir, "2MPP_N128_DES_V1", "resimulations", "R512") files = glob(join(fdir, "mcmc_*")) files = [int(search(r'mcmc_(\d+)', f).group(1)) for f in files] + else: + raise ValueError(f"Unknown MANTICORE simulation `{simname}`.") elif simname == "csiborg2X": # NOTE this too is preliminary fname = "/mnt/extraspace/rstiskalek/MANTICORE/resimulations/fields/2MPP_N128_DES_PROD/R512" # noqa @@ -134,7 +136,7 @@ class Paths: for file in files if search(r'realization(\d+)_delta\.fits', file)] # noqa files = [int(file) for file in files] # Downsample to only 20 realisations - files = files[::5] + # files = files[::5] elif simname in ["Carrick2015", "Lilow2024", "no_field", "CLONES"]: files = [0] elif "IndranilVoid" in simname: @@ -214,9 +216,21 @@ class Paths: f"chain_16417_{str(nsim).zfill(3)}", "output", f"snapshot_{str(nsnap).zfill(3)}.hdf5") elif simname == "csiborg2X": - raise ValueError("Snapshots not available for CSiBORG2X.") + raise ValueError("Snapshots not available for CSiBORG2X based on " + "Stephen's ICs.") elif "manticore" in simname: - raise ValueError("Snapshots not available for MANTICORE.") + fdir = self.manticore_dir + if simname == "manticore_2MPP_N128_DES_V1": + if nsnap != 1: + raise ValueError("Only snapshot 1 is available for " + "MANTICORE 2MPP_N128_DES_V1.") + + fdir = join(fdir, "2MPP_N128_DES_V1", "resimulations", "R512") + nsnap_str = str(nsnap).zfill(4) + fpath = join(fdir, f"mcmc_{nsim}", "swift_monofonic", + f"snap_{nsnap_str}", f"snap_{nsnap_str}.hdf5") + else: + raise ValueError(f"Unknown MANTICORE simulation `{simname}`.") elif simname == "quijote": fpath = join(self.quijote_dir, "fiducial_processed", f"chain_{nsim}", @@ -407,6 +421,11 @@ class Paths: else: raise ValueError(f"Unsupported CSiBORG2X field: `{kind}`.") + if simname == "manticore_2MPP_N128_DES_V1": + basedir = join(self.manticore_dir, "2MPP_N128_DES_V1", + "fields", "R512") + return join(basedir, f"SPH_{nsim}.hdf5") + if simname == "Carrick2015": basedir = "/mnt/extraspace/rstiskalek/catalogs" if kind == "overdensity": diff --git a/csiborgtools/read/snapshot.py b/csiborgtools/read/snapshot.py index aaf70d6..a651dac 100644 --- a/csiborgtools/read/snapshot.py +++ b/csiborgtools/read/snapshot.py @@ -517,6 +517,91 @@ class CSiBORG2Snapshot(BaseSnapshot): } +############################################################################### +# CSiBORG2x snapshot class # +############################################################################### + + +class CSiBORG2XSnapshot(BaseSnapshot): + """ + CSiBORG2X snapshot from the SWIFT N-body simulations provided by Stuart + based on the Manticore ICs. The snapshots are at `z = 0` and correspond to + `manticore_2MPP_N128_DES_V1`. + + Parameters + ---------- + nsim : int + Simulation index. + nsnap : int + Snapshot index. + kind : str + CSiBORG2 run kind. One of `main`, `random`, or `varysmall`. + paths : Paths, optional + Paths object. + keep_snapshot_open : bool, optional + Whether to keep the snapshot file open when reading halo particles. + This is useful for repeated access to the snapshot. + flip_xz : bool, optional + Whether to flip the x- and z-axes to undo the MUSIC bug so that the + coordinates are consistent with observations. + """ + def __init__(self, nsim, paths=None, keep_snapshot_open=False): + nsnap = 1 + flip_xz = False + super().__init__(nsim, nsnap, paths, keep_snapshot_open, flip_xz) + + fpath = self.paths.snapshot(self.nsnap, self.nsim, + "manticore_2MPP_N128_DES_V1") + + self._snapshot_path = fpath + self._simname = "manticore_2MPP_N128_DES_V1" + + def _get_particles(self, kind): + with File(self._snapshot_path, "r") as f: + h = f["Cosmology"].attrs["h"][0] + + if kind == "Masses": + npart = f["Header"].attrs["NumPart_Total"][1] + mpart = f["Header"].attrs["InitialMassTable"][1] * 1e10 * h + x = np.ones(npart, dtype=np.float32) * mpart + else: + x = f[f"DMParticles/{kind}"][...] + + # Convert coordinates to Mpc / h + if kind == "Coordinates": + x *= h + return x + + def coordinates(self): + return self._get_particles("Coordinates") + + def velocities(self): + return self._get_particles("Velocities") + + def masses(self): + return self._get_particles("Masses") + + def particle_ids(self): + raise NotImplementedError("Recovering particle IDs is not implemented " + "for CSiBORG2X.") + + def _get_halo_particles(self, halo_id, kind, is_group): + raise NotImplementedError("Recovering halo particles is not " + "implemented for CSiBORG2X.") + + def halo_coordinates(self, halo_id, is_group=True): + return self._get_halo_particles(halo_id, "Coordinates", is_group) + + def halo_velocities(self, halo_id, is_group=True): + return self._get_halo_particles(halo_id, "Velocities", is_group) + + def halo_masses(self, halo_id, is_group=True): + return self._get_halo_particles(halo_id, "Masses", is_group) * 1e10 + + def _make_hid2offset(self): + raise NotImplementedError("Recovering halo offsets is not implemented" + "for CSiBORG2X.") + ############################################################################### # CSiBORG2 snapshot class # ############################################################################### @@ -771,42 +856,90 @@ class CSiBORG2Field(BaseField): class CSiBORG2XField(BaseField): """ - CSiBORG2X `z = 0` field class. + CSiBORG2X `z = 0` field class based on the Manticore ICs. Parameters ---------- nsim : int Simulation index. + version : str + Manticore version index. paths : Paths, optional Paths object. By default, the paths are set to the `glamdring` paths. """ - def __init__(self, nsim, paths=None): + def __init__(self, nsim, version, paths=None): + self.version = version + if version == 0: + self.nametag = "csiborg2X" + elif version == 1: + self.nametag = "manticore_2MPP_N128_DES_V1" + else: + raise ValueError("Invalid Manticore version.") + super().__init__(nsim, paths, False) def overdensity_field(self, **kwargs): - fpath = self.paths.field( - "overdensity", None, None, self.nsim, "csiborg2X") - with File(fpath, "r") as f: - field = f["delta_cic"][...].astype(np.float32) + if self.version == 0: + fpath = self.paths.field( + "overdensity", None, None, self.nsim, self.nametag) + with File(fpath, "r") as f: + field = f["delta_cic"][...].astype(np.float32) + else: + raise ValueError("Invalid Manticore version to read the " + "overdensity field.") return field def density_field(self, **kwargs): - field = self.overdensity_field() - omega0 = simname2Omega_m("csiborg2X") - rho_mean = omega0 * 277.53662724583074 # Msun / kpc^3 - field += 1 - field *= rho_mean + if self.version == 0: + field = self.overdensity_field() + omega0 = simname2Omega_m(self.nametag) + rho_mean = omega0 * 277.53662724583074 # Msun / kpc^3 + field += 1 + field *= rho_mean + elif self.version == 1: + MAS = kwargs["MAS"] + grid = kwargs["grid"] + fpath = self.paths.field( + "density", MAS, grid, self.nsim, self.nametag) + + if MAS == "SPH": + with File(fpath, "r") as f: + field = f["density"][:] + + field /= (681.1 * 1e3 / grid)**3 # Convert to h^2 Msun / kpc^3 + else: + field = np.load(fpath) + else: + raise ValueError("Invalid Manticore version to read the " + "density field.") + return field def velocity_field(self, **kwargs): - fpath = self.paths.field( - "velocity", None, None, self.nsim, "csiborg2X") - with File(fpath, "r") as f: - v0 = f["v_0"][...] - v1 = f["v_1"][...] - v2 = f["v_2"][...] - field = np.array([v0, v1, v2]) + if self.version == 0: + fpath = self.paths.field( + "velocity", None, None, self.nsim, "csiborg2X") + with File(fpath, "r") as f: + v0 = f["v_0"][...] + v1 = f["v_1"][...] + v2 = f["v_2"][...] + field = np.array([v0, v1, v2]) + elif self.version == 1: + MAS = kwargs["MAS"] + grid = kwargs["grid"] + fpath = self.paths.field( + "velocity", MAS, grid, self.nsim, self.nametag) + + if MAS: + with File(fpath, "r") as f: + density = f["density"][:] + v0 = f["p0"][:] / density + v1 = f["p1"][:] / density + v2 = f["p2"][:] / density + field = np.array([v0, v1, v2]) + else: + field = np.load(fpath) return field diff --git a/notebooks/flow/reconstruction_comparison.ipynb b/notebooks/flow/reconstruction_comparison.ipynb index 5b208f9..f512100 100644 --- a/notebooks/flow/reconstruction_comparison.ipynb +++ b/notebooks/flow/reconstruction_comparison.ipynb @@ -228,13 +228,104 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:18:58\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:18:52\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:20:01\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:20:05\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 14:49:46\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:19:42\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_LOSS_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 21:55:56\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:18:47\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:18:56\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:19:59\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 13:19:50\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 14:50:07\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:19:26\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_Foundation_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 21:55:55\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 15:27:56\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 11:29:14\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 12:03:40\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 12:04:05\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 14:52:59\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:19:10\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 21:59:45\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 15:28:34\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 11:29:14\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 12:05:48\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 30/08/2024 12:06:16\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 14:57:22\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:19:11\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 22:02:16\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:48:49\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:49:07\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 11:13:55\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 11:12:57\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 15:05:21\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:19:57\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 22:15:58\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:48:29\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:49:00\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 11:10:51\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 11:07:55\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_manticore_2MPP_N128_DES_V1_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 15:04:25\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 12/09/2024 10:29:11\n", + "File: /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_w1_mike_zcmb_max_0.05_sample_alpha.hdf5\n", + "Last modified: 20/09/2024 22:07:15\n" + ] + } + ], "source": [ "zcmb_max = 0.05\n", "\n", - "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\",]\n", + "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"manticore_2MPP_N128_DES_V1\", \"CLONES\", \"CF4\",]\n", "catalogues = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n", "\n", "y_BIC = np.full((len(catalogues), len(sims)), np.nan)\n", @@ -255,9 +346,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAGYCAYAAACUM6fLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB6IUlEQVR4nO3dTWwjZ5of8L+IcZAYGatEYZVeKMhERQ8QYLAZu0hNDgGCZFXVkyCwAtgscQUESQ4j0kYGCBCMVa0FgkUuqyHblwCLjEn5sEGQ1Uikexe6rcnuIIccMhKrvYc9BDMsGUjUETRoqtRe+DI7rBy09br4KZJiVbHI/w8w3CwW+ZYk1sP383kXHMdxQEREREREFEGxsC+AiIiIiIhoXGzQEBERERFRZLFBQ0REREREkcUGDRERERERRRYbNEREREREFFls0BARERERUWSxQUNERERERJHFBg0REREREUUWGzRERERERBRZ3wj7AoiIiIjGZds2JEkSj2u1GizLQjweh2VZ2N3dHes4EUXHguM4TtgXQURERDSKSqWC09NT1Go11Ot1cVzTNFSrVQBAoVCALMtIp9MjHyei6OCUMyIiIoqcdDqNXC7XdqxWq7WN1iiKgqOjo5GPE1G0sEFDREREM8E0TcTjcfE4Ho/DNM2RjxNRtMzFGprvfOc7SCQSYV8G0dxoNBr48z//89DKz+VyKBaLAEabH89YQRQcP+LEy5cvJ3J8kL/9t/82/vIv/1I8Xl1dxerq6sjvQ0S9XVxc4OLiQjz+xje+gf/7f//vwNfMRYMmkUjg5OQk7Msgmhubm5uhlV0oFGBZlnicz+fb5sdXKpW+8+MZK4iC40ecWF5ehm3b9z4+iKIojBNEARomVnDKGRHNDMuyurIdcX480fxQFAXNZlM8bjabUBRl5ONEFC1s0BDRzKjValBVVTzm/Hii+aKqatsIrWmayGQyIx8nomiZiylnRDT7arUatra22npbR50ff3Fx0Ta0vb29je3t7YldI9E8Ozw8xOHhoXjsnSM/jlqthmq1CsuyUCqVoKoqZFlGPp8X6ZcBiCmmox4nouhgg4aIZoK7uZ63QTPq/PjV1VXOjSfySWcHwX3X0KiqClVVkc/nex7vd/6wx4koOtigGdLNzQ2++uqrsC9joNdffx2Li4thXwZR4AqFAiRJQqVSgWVZsCwLlUoFiqLg9PRUnOf3/PgoxImwMD4RfS2KsYL3ME0zNmiGcHNzgz/4gz/Ar371q7AvZaDXXnsNP/zhDxlw5lSpVEIqlfKlwm7bNgzDEKmQp403FbNpmqhWq2LaiGEYbc/5NT8+KnEiLIxP02OeY8U0iGqs4D0cPD/v1V6ifP+yQTOEr776Cr/61a/w7rvv4jd+4zfCvpyefvnLX+LJkyf46quvGGzmkGVZaDQayGazAG6DUqlUgqIobVMpDMOAZVmwbRvlcrktAxhwGzy9gcw0TTiOA0mSoGkaCoXCwH1cwmZZForFohihSafTgc2Pj0KcCAvj0/TojBW5XE5M08zn8+I+6Xfcq1c8iUqsCFMUYwXv4eD1ulfdBBb5fF40cnK5HM7OzgAA5XK57V6t1WrQNA3X19dd3/fA7ewGd1sDb53g3/27f4f/9J/+k18/WpuJNaKcOfDOO+/c6/UvXrxwfu/3fs958eLFhK5o8qJwjeSfdDrtNBoN8TibzTqKojjValUcq1arTjabFf9Op9MD37NYLDrFYrHtmCzLQ13Pfe+5sNznunkP9sffzfTwxopGoyFiRLlcFvGh33Gvu+LJMLFiHuOE40TzfojiNUed914tl8vO7u6u4ziOc3197aiqKo5771X3uEtVVUdRFOf6+rrr/ev1urhv3X83Gg1nd3dX3L/ZbNZRVdVRVdWp1+vitdfX104+n2+rYzQaDUeSJHF+Z/3BcW7rFYqiiP/cZki5XHby+Xzf38Uw9xzTNhPNAMuy2nplisUiUqlU2znVahXJZBLA7SLYQemL3R5Xt2fIJUnSyJvQEdH08MYKWZbFCO7R0RE0TRt43OuueMJYQXQ/3ntVURSYpgnbtlGr1cToTDqdFvdq555KhmHAMIy2rQu83NEb97W1Wg2GYSCXy0GSJPyX//JfIEkSqtUqyuVy2/RtwzB67ummqiqq1Sqq1WpX/QEAstks6vU66vU6crmcGJVJp9P3HqFhg4ZoBgxTcdA0DfV6HcBtIPPuvdDJDWqdUqmUGNomoujpjBWWZSGZTMKyrLbpqf2Ou+6KJ4wVRPfjvVdlWYaiKNjY2MD+/j729va6zs/n8+K4aZp9711Xo9Hoauy4jahUKoVYLNazEQX07jR1X28YBgqFwp0/W2en6X07QdigIYo4N13xXdwN5HRd75pn2+n4+Jh7MRDNmF6xQpZl1Ot15PN5bGxs3HncNUo8IaLRdN6rlUoFtm2jXq+jXC533ZOFQgGaponv7f39fXF/np2dQdf1rsaCJEldHRHec37zN3/zzkaUVzweRzweF+f16hR19eo0vW8nCBs095BIJFCpVPo+b1kWNE1DIpFAMplEqVRqe940TSSTSSwtLSGRSHRlY+r3HJHXKL0a7tBxMpnsG2xM0+xbOemc2kbD8d7L7n+6rvtebqlUGit22LaNhYUFThmaMZ2xwjtNLB6Pi8pNv+OdBsUTxorxlEolESO8966maVhaWur5mmQyiUQiAeC2YuuNM51xx/3bhhWTaDid96p3KhnQ3vAwDAOqqrZ1QpbLZdH4kWW5ZxKgTCYjtjUwTRPf/e53xTmWZeF//+//PbAR1euaq9UqJElCOp1GrVbre64fnabMcuYT0zSxsbGBcrkMVVVh2zZ0XUe9XkexWIRt213PHx8fA8DA5ygagk6LKstyW4+OrutiyNmyLDGsq2kaJEkSu2m7EokEGo0GgNvAyQbNZNm2jYODg6kd9er8u0qShHq9PtTIH0WLN1bIsiyymVmWhadPn4pzeh0H2mNFv3gCMFaMw7IsVKtVNBoN2LYtGipu/I7H4yJ7o/c13srt7u6uyC7ndoxeX193lTXtMYna79VsNotcLgdN09BsNlEulwHc1jVKpZJoPMTjcZG1rB/3HlYUBbIsi/v4T/7kT8S6OMuy8Nprr7W97q4OLrczVJIkmKbZt/7Tr9P0vjGDDRqf7Ozs4ODgQMxfdFuuS0tLbekwvc+7Qcsdcuv1HE2/YVMo90vB2Plebi9dJpNBOp3umRZV13UcHx+LMt1g16lfoHMrKED/XbNN0+SX34wyDAN7e3ttn8Gg9j2gYHljhSRJPRfi9jsOtMeKfvGEsWI8tm2L6TqSJCGXy7UtrnYXTnt/t8VisW1xNc2Ozu/1Xn/jbDZ7Z/3QXevm8t7DnR0Rsizjf/yP/4F0Ot23EeVeW2enqSRJ2NnZEQ0f7/nDdJret0HDKWc+sG27b0BXVRVHR0eiwqjretew3KDnaPp1zg3tlQ2kUqn0zR7ipWkayuUyyuVy2+epMyOImznEzylCxWKxK/jR/VUqFSSTya5pOwsLC+LfhUKh7bmlpSUUCgXRg+udElQoFLC0tCQ+O166rovpJW5syeVyqFQq0HVd9M51lt/vGgddB00nxorppShKV0eCd9G2pmldIzKdIzY0O4K4Vzvpuo7/8B/+g7h/i8UiqtUq6vV622ezXC6j0Wi0NbjdqW1uljPvCH9np2nnd9MkOkHYoPHB2dlZ36kasiyLL/3z83MAtx+gpaWltvU4g56j6TZMCuV+KRi9KpWKmPqhadqdaVH97KGzbZsVlHsyDEM0CtzsUaZpYn9/H0+fPhVTvIZZ82Lbtli0rSiKiA/u+52fn6NarXb1duVyOTQaDZTLZTFfvlgsii+Yzp4873v2usZ+10HTjbFi+tm2jWKx2BUP0um0mIJeq9UGZrG6S6+YRNMl6JG3ra0tMdUtSJPoBOGUMx+kUqm+HwZvZVeSJNFKdXtIG42GmIPY7zmabsMEAm8KRgBtc9Rdp6enaDabqFarME0ThmG0TfFwM4J4v9D8WvPAtRT3l8/nu3qgDMNAJpMRv9+9vT2sra0NFdjdv/v6+rro/arVamLoH7gd5n/58mXXaxRFgW3bQ2XIOzo6GniNva6Dph9jxXTrlz0ul8tB13Vks1kUi8U7M08N0ism0fQJ8p6SJAl/+Id/GGiDZlKdIFPfoKlUKmg2m2KRnNtadfPeuxlY3LUE08BdJFkqldrmNtq2jUql0rMXNJ1Oi177zgA26DmaLsOmUPamYLQsCxsbG12fi0QigeXlZQC3FVDu6UBe43zJFQoFkdUmzOsgov50Xe+7rtKtA5imKRZec2SFJi3oRtQkTPWUM+9io3w+j2azKTbryefzyGazondh2qY6uOsi3HnqbpazbDYrdmQtFAqiFew20FRVHfgcTbdhUygPSsHoSqVSYkSm12I5ZhGKvkwmI7IeArcjNltbW+J59/hdWWtcqqqKxnLn6yqVilizdXBw0PY6dz+CXp/Du66RiCZH1/WuBB2d3FGaQft8EM2bqW7Q2Lbdtph6fX0d1WoVtVqtrUWnKErXouug7OzstOVxdxtWiqLg6dOnyOfzYh8aTdPECFMqlcLLly+xtraGpaUlGIaBp0+fQpKkgc/R9Oucf+omd8jn82IvIm8aZXdqgcvdT0BRFGiaBk3TkMvluhbRsUETfYqiIJ/Pt+0j4caIbDaLjY0N6LoupqEO837pdBpra2vide4on7sRYjKZxM7OTttnR9M07OzsYGNjo6tRM+gaiWhyKpUKKpUKNjY2sLS0hKWlpZ57w2SzWciyzOynRB4LjuM4YV/EsHK5HCRJwvLyMhqNhvhSNU1TrDHpZXNzEycnJ2OX+//+3/8TqRF/8zd/c+z38VMUrnFedDZa/GCaJo6Oju4177TVaokpcrHYZPs27nvPheU+1817sD/+bqKLcaLbPNQpOkXxmueVn/dsWIa55yLzk7rZoPb29toWuQ7j4uICm5ub4r/Dw0OfrpIoGmlRT09Psbq6iuXlZayurt57XcXh4WHbPXZxcXGv9yOi8E06ThCRv8K8Z1utFprNJlqtVmBlek19UgDXzs6OyGu9vLw8UmVxdXU1kr1AFF3eNQeTdt+MIK1WC5ubm7i6ugIAXF1dYXNzEy9evGjbe2QU29vb2N7eFo83NzfHvj4iCp8fcYKI/BPmPXt6eorNzU1cXl7iwYMHODk5wfr6uq9ldorECE2hUEA+nxd7uCiK0raoutlscldrmjp+pkW9z3vbto3Ly0vRi9JqtXB5eRl43nkiml6ME0TREtY9268hFfSKlqlv0FQqFaiqKhawuhtJedMUmqaJTCYT1iUSRYokSXjw4IGYWxuLxfDgwQMmnSAigXGCKFrCumenpfNjqhs0lmVB13Ukk0ksLCxgYWFB/ILy+TwKhYLIKsbNoYiGE4vFcHJygpWVFQDAysoKTk5OOI2EiATGCaJoCeuenZbOj6leQyPLct8hK1VVuS8L0ZjW19dxcXGBm5sbSJLESgoRdWGcIIqWMO5ZtyHlrqEJq/Njqhs00+aXv/xl2JfQ1zRfGw0WVorFWCyGpaWlwMqbF7wXu/F3El2ME/6J0n0RpWudd2Hcs9PQ+cEGzRBef/11vPbaa3jy5EnYlzLQa6+9htdffz3sy6ARTENmEJqMqMSJsDA+Ed2KaqzgPUyDhN35wQbNEBYXF/HDH/4QX331VdiXMtDrr7+OxcXFsC+DhsS0qLMlKnEiLIxPFCTTNHtmPzVNE5IkIR6P4+zsDKlUCpIkoVarwbIsxONxWJaF3d1d364tqrGC9zBNMzZohrS4uMgbmSbKzQzi8mYG4RSPaGKcIJoOyWSy7bEkSTg/P8f+/j4qlQokScLe3p5Yi5vP51GtVgFAJBzyM9kQYwXRZE11ljOiWTYtmUGIiGaJaZqo1+twHAeO46Ber6NcLkOSJGiaBsdxcH19LUZharVaW9xVFAVHR0chXT0RjYMjNEQhmZbMILOiVqshHo+j2WyiXC6jWCyK40FNJSGi8HVONTs7O0M2mxWPLcuCbdviPNM0EY/HxfPxeBymaQZzsUQ0EWzQEIVoGjKDzAq35xW4nT5SKpWQzWYDn0pCRNPDMAzk8/m2Y7ZtQ5Zl6LqOg4MDvHz5cqT3vLi4wObmpni8vb2N7e3tiVwvEQGHh4c4PDwUjy8uLu58DRs0RCELOzPIrGg0GuLfzWYTqVSq51SSYrHIBg3RHLBtG5ZltR3zjtRomgbDMJBIJEba1Xx1dRUnJyeTukwi6tDZSeDtQOiHa2iIaCbIsgwAqFQqUFUViqJwKglNvVKp5Ntn0rZt5HI5X947CkqlkogLwO30U03T2s5pNptQFAXNZrPrGBFFBxs0RDQzTNNEs9nE8vIyAIw9lcT9zzvkTTRplmWh0WiIynMul0MymUQymRQjC5ZlYWlpCZqmQdM0lEqlnu9VKBTEOe5og7sIvlAoBPLz3OXw8LDt/hpmGsl9nJ6eIpFIiMeyLLc18Or1OjKZDFRVbRvJMU0TmUzG12sjosnilDMimhmKokBRFOi6jlwux6kkNNW86zsqlQp0XUexWESlUkEulxNrv1RVRblc7vs+pmni9PQU1WoVpmliZ2dHnJ9Op5FIJKYiGcY400juyztCI8syTNMUjcJEIiGmn+bzeRQKBXE+p6USRQsbNEQUeZVKBcViUVQA19fXcXR0BF3XcXp6Ks7jVBKaJpZl9axAd06BsiwLhmFgeXm5Z8PEO5VKURTUarW25yVJgm3bc5cSvlcjsF9DRVVVsScNEUUPp5wRUeRJkgRd18Xj09NTTiWhqddv9DCfz2Nvbw/A7bqveDwuHvdaE9NoNNrWinVKpVI4Ozu7/wUTEU0pjtAQUeSpqopKpYJKpYJmswlZlkVPNqeS0DTqN2LiroVxP6eSJImRx3Q63bWo3T2nM5sXEdE8YYOGiGYCp5JQlLjTwLwMw0Amk2mbFmmaJmRZhiRJME2z55TJTCaD/f19cX4qlWp73ju1jYhoFrFBQ0REFAJZlsVITalUQqlUEutf4vE4qtUqJEnCzs6OaPx414UkEgmRJU2WZWiaBkmSutaOsEFDRLOODRoiIqIQ6LqO4+NjZLNZ8V8nWZb7ZjjzbibrZkvrZJomp1kS0cxjUgAiIqIQZLNZ1Ov1kVKLj6pYLPZt7BARzQo2aIiIiEJSLBZ9e2/bttmYIaK5EIkGTWfvlWmasCwLtm2jVqv52rtFRETkJ7/2h5Ekae72niGi+TTVDZpKpQLDMLCxsdF2fH9/H4lEAmtrazBNkwGbiIiIiGhOTXVSgHQ63XPXY03T+i6SJCIiIiKi+THVIzSDWJYF0zTDvgwiIiIiIgqRbw2aJ0+e+PXWAG7X1ciyDF3XuYaGiIhmXqvVQrPZRKvVCvtSiIimim8Nmp/+9KfIZDJdxw8ODvDRRx/d672z2SwURYEkSdA0DYZhDDz/4uICm5ub4r/Dw8N7lU9E7Q4PD9vusYuLi7AviWimnJ6eYnV1FcvLy1hdXcXp6WnYl0RENDV8a9Ccn59ja2sLmUwGr169EiM2Ozs7WFxcxCeffDLW+9ZqNWia1nas2WwOfM3q6ipOTk7Ef9vb22OVTUS9bW9vt91jq6urYV8S0cxotVrY3NzE1dUVAODq6gqbm5twHCfkKyMimg6+NWgsy8J7772HYrGI3d1dHB0died2dnbGXtQvyzJyuZx4XK/Xe44EERERzQLbtnF5eSmmmrVaLVxeXnK6NRHRX/Ety9nS0hKA2zz4H3/8Md588018+eWX+OY3vwngtsFzl1qthmq1CsuyUCqVoKoqZFmGaZoolUoAgEQigXQ67dePQUREFCpJkvDgwQNcXV2h1WohFothZWWFWxYQEf0V3xo0qqq2PU6n0/j93/997O/vA8BQQ+WqqkJV1a6djtmAISKieRGLxXBycoLNzU1cXl5iZWUFJycnWFhYCPvSiIimgm9TzgzDwAcffCAe7+3tYX9/H3t7e3j16hV7loiIiIa0vr6Oi4sLNJtNvHjxAuvr62FfEhHR1PCtQbO2toZsNovPP/8cALC4uAgAePToEXZ3d7G8vOxX0URERDMnFothaWmJIzNERB183Vjz7bffxltvvdV2bHFxEbu7u6KBQ0RERERENC5fGzT9yLLctS6GiIiIiIhoVKE0aIDbKWlERERERET3EVqDxt1ok4iIiIiIaFyBNGg+/fRTpFIpfPvb38a3v/1tvPnmm9jZ2Qmi6JGUSiWYpunLe9u23bYhKBEREdEk+Fl/6YV1Gpo2gTRoqtUqyuUyzs7OxH/T1qCxLAuNRgOKogC4vVkLhQJqtVrbeYVCAZqmQdO0vrs09zpHkiRomoZCoeDnj0FERERzpLP+AtxuTL6wsNBWT+lXrzEMA7quD6zX5HI5JJNJJJNJWJbFOg1NnUAaNJqmYW1tDYuLi1hcXIQkSfjd3/3dIIoemmEYbb0NhmHg6Oio7RzTNHF6eopqtYp8Pt+zUTbonHQ6jWKx6N8PQURERDBNE5ZlwbZt1Go1UVGv1WoolUqoVCptlfF+x6Ogs/4CAPl8vq2B457XWa9xfzflchmGYfSs11QqFei6jnq9jr29PVEW6zQ0TQJp0FxfX+Ojjz7CkydP8OzZMzx58mQqR2hkWRaPi8UiUqlU2zm1Wg2apgEAFEXp6uUY5hxJkvr2gBAREdH97e/vI5FIYG1tDaZpis288/k8stks0uk0gNvK+qDjUdBZfzEMA4ZhIB6Pt53Xq15TrVaRTCYBAKqq9py2lk6noaoqgNt6TbPZFM+xTkPT4ht+F3Bzc4Mf//jHUFUVv/jFL8TxIOd6DmOYG7LRaIjGyrjnpFIpnJ2dieBARJNRqVTQbDbRaDRg27boOazVarAsC/F4HJZlYXd3N+QrJSK/aZqGcrncdqxWq4mGDXBbOS8Wi5Akqedxt3Ez7bz1F3dkSlXVobbH8P6e3Fg5SD6fx97ennjMOg1NC99HaBYXF1EsFvHxxx+3/Xd8fOx30UOzbbstmPUjSdKdN/sw5xDRZFmWBcuykM1mkc/n0Ww2xbSRKPe8EtH4LMtq6zw1TbNt1CIej8M0zb7HxxXkAn3btvHNb35TTAPb39+HZVnQdR1nZ2fQdX1gh62qquL8crncNtLTyV0fHJWGHs2XQKacbWxsdB1LJBJBFD2UYYdMM5kMTk9PAdwGxs6h22HO6RwaJqL7s227bW74+vo6qtVqzx7ZzjnkRDSbbNuGLMuiUv/y5cue5/U7Po7OBfq5XA66rkPX9bbOzlwuJ5IHeRs/nYvvO5VKJfF8MpnE0tISvvzyS7FAv1wuo16vi8ZJuVy+s8PWTdyUTCb7Zi4zDAOqqnY1ZsKo0zCjG/Xi25SzZ8+eDXy+XC7jJz/5iV/Fj0yW5baRGl3XxdCt2/OrKApkWYamaZAkqW04O5FIiCDW7xyADRoiPyiKgnq9Lh679+Kke16JKBqy2az4t6ZpMAwDiUSiZ+fl8vLySOtALi4usLm5KR5vb29je3sbwG3F353q5Y58qKqKSqWCfD6PYrGISqUCSZJQrVZh2zZ0XUe1WhWL791zcrkcqtVq18/l/mylUgnAbX1KVVUkk8mBU2p71Wvc348kSZBluW2amluvKZVKKJVKYk1wPB4X1xV0ncZtMLrXnsvlxJqefD4PWZZRKpXakhWYpgnHcfqe79XvtW6DkVOWg3F4eIjDw0Px+OLi4s7X+NagyWaz0DRNfIg6eSsf00DXdRwfH4ubpLMh4uo3J7XRaNx5jmmaHKol8pmb1aher2N/fz/syyGigNVqNeTz+bbGQLPZhK7rYgaFe0xRFCiK0vN4P6urqzg5Oen5nLeCL8uy+PfR0REymQyAr9fouLHKLctbP+hcfN/JzUzm/ozHx8ditonbMdtZz+pXr+lsNLnceo23EeUVRp1mmAZjr0bfoPO9+r02nU4jkUiwQRMQbycBgLYOhH58a9AUi8WeU81cz58/96vosWSzWeRyuaHX04yjWCwyxSGRz3Z2dlCtViFJ0kR7XonofsbpdR2HLMttU4Tq9ToymQxUVYVhGOK4aZoDj4+jM964lWgAYuG8LMtQFEXUkZ4+fdr1Pp2L7zt5UzW79Zff+q3fCnSBfhh1mmEajK7ORt9d5w96LYCuBiNNF98aNIMaMwDw9ttv+1X02NweEz/Ytj1UxhEiGl+hUBDTCCzLmmjPKxHdzzi9ruOQZRmmaYoe9kQiIUYS8vk8CoWCqNjedXwUvSq7siyjXq+jVqthY2MD9XodlUoFtm2jXq/Dsixx3DXM4vvj4+O2xkSxWMS/+Tf/ZuRrHldYdZphGoyuXvvzDDr/rtcyo9t08z1tc9T41fJmi57IX5VKBaqqigpJrVZDNpudWM8rEUVHv8aAqqo9K6T9jo+iM8GQaZqiA8VNGw+gayqZ9zWGYSCTyQzseDFNs+e6lYuLi8DWs4RRpxm2wejqbPTddb5Xr9fSdGODhogiz9vr5nJ7DyfR80pENAxvgiF36luz2YRlWWJqmTtFTNM0NJtNsbZl0OJ7d4E+cNsg6tVwmfWkQ8M2GN3nOn8Xg87vPG8ef79RF0qD5smTJ9jf32+bCkJENC5ZlvsmIJlEzysR0TC8CYYkSerby9/reL/F90B74qFeMW1ekg4N02AEejf6Bp3PBmP0BdqgefbsGQzDQKPRwM3NzdCv6xxm5M7fRERENG2CSDDUy7wkHRq2wdir0TfofDYYoy+QjTWfPXuGVCqFdDqNTCaDZrPZtzfVq1KpwDCMrgQD07Dzd6vVQrPZRKvVCrxsIiIimk5BNyzus0A/anWZbDaLer3uWwKnforFIhM7TTlfGzRuQ0bXddGQ+dGPfgQAWFhYuPP16XS6K8vENOz8fXp6itXVVSwvL2N1dZVT54iIiEgIcnRGkqSxyotqXSZKDUYKjm8NmocPH0LXdfzO7/wOXr58iQ8//HAi7xv2zt+tVgubm5u4uroCAFxdXWFzc3OoESciIiKisEW9LhOFBiMFy9eNNd00qpP08uXLkV8zyc3ybNvG5eWleNxqtXB5eQnbtrG0tDTWexJFXVAb5hER0f2xLkOzxrcGzdraGj788EOcn5/j4OAA6+vreOutt+79vqPu/A1MdrM8SZLw4MEDXF1dodVqIRaLYWVlha13mmtBbZhHRET3x7oMzRrfkwKsra1hZ2cHjuPg8ePH+Pzzz+/1foqitG1KddfO35MWi8VwcnKClZUVAMDKygpOTk6GWhNERERE1Cnoxfmsy9CsCSTLGQC8/fbb+PDDD0XDZtx5mqqqdm2eFPTO3+vr67i4uECz2cSLFy+wvr4eaPlEREQ0G8JanD9vdZmoZXSj0QTWoHG5DZt6vX7nubVaDcViEZZloVQqiYaMu/O3m645jNzgsVgMS0tL7M0gIiKisYS9OH9e6jL3bTSWSqVAE1DZtt2V5ZcGC7xB43r77bfvPEdVVeTzeVxfXyObzYodWlVVxe7uLtLp9ExtqjmJG4Y3ARERUTS4i/PdUQPv4nyajPs2Gi3LQqPREMsbcrkckskkkslk24why7Kg6zp0XW/bHzGXy4nj3vNdpVJJvF8ymcTCwgIkSYKmaSgUCvf50edKoA2aTz75JMjiIqXzhul3Y3jZto1CoYBarSaO8SYgIiKKBndxfix2Wx2LxWJ48OABF+dP0H0bjYZhiI7iSqUCXddRr9ext7fX1oGsaRrK5TLK5bKYOeTW5crlMjKZTM/9bNzNQuv1OnK5nNhnJ51OB77nTpT5luWsl0ajEWRxkWIYRtsHXdO0O39fhmHg7OysKylCOp1GIpGYqdErIiKiWeMuzt/c3MTl5SUX5/vgvhndLMsSM4S8Sxy8SaoqlQpkWUYul4NlWcjn81AUBbIsi9ceHR0NXPNt2zbK5TKq1Wrbtdu2zQbuEEKbckbtvDeM98bQNK3vNLRisYhUKtXzOfcmICIiouk1b4vzg3bfjG796lL5fB57e3sAbtfoNJtNFItF5PN5GIYhzrMsS0xPG7Q3o3ckyJVKpXB2djbUdc67QEdoqD/vDePeGNVqFaZpwjCMthb7MNybYNIbmxIREdFkuYvzyR9uo/Hm5gaSJI3UmOk1OlIoFKBpmhixSSQSWF5eBnA7cuNthMiyjHq9jlqtho2Njb5JsY6PjznF7B44QjMFOm+YRCIhhiU7bwwiIiIiGs04Gd16zXYxDAOqqrZNP0ulUqLj2TvjxjvDJh6P90wK4J7nvsbL+140GEdopkDnDZNKpWAYBnZ3d8f+MPMmICIiIrofWZZFx3OpVEKpVBLJmOLxOKrVKhRFgaZp0DQNAFAul8Vrc7kcms0mLMvC06dPxfsmEgmxVrrZbLJBc09s0EwJ7w3T78YA2m8AXddhmiYsy4JlWchms+I83gTdSqUSUqlUVxKFUdi2DcMwOCxMoeJnmYgoGLqu4/j4GNlsVvzXy+7ublcyJkmS+sZYb+InVVW7lgiYphnKPotRFeiUs3g8HmRxkeLeMK7d3V1Uq1VUq9W2hon3BiiXy2g0GqhWq203WBRugqD33PGmxbYsC0tLS6LRWCqVer7GnSOraZoYQWNabApbrz0R3M/poHuqVqthYWGBn2UiohG4aZWDTrTkJhig4QTaoPnwww+DLC5SJnnDTPtN0FkhA7orW516VdpGqZB1Zg9RVVU0GHv1tpimidPTU1SrVeTzeezs7IjnmBuewtS5J4IkSahWqyiXy22ZdTq5aUS9+FkmIrpb0HHStu2prsdNIyYFmCKTuGGicBP0Sk3Yq7LlGlRpG7ZC1jkFz7IsGIbRtzFUq9XElD9FUdo2LwWYFpvC4/0sK4oC0zRh2zZqtVrfe8gwDBiG0XOUnJ9lIqK7BbkXjCRJ3HtmRGzQTJn7foCjcBN0Ni4GVbaAuyttw1TIvM/H43HE43GRP77XtLVGozFwiiRzw1NYvJ9lWZahKAo2Njawv78vPtNe7jq7finc+VkmIqKoY4OGAuetkN1V2QLurrTdVSHrTIvtjvZIkoR0Ot01+uKe0y+9IlFYOj/LlUoFtm2jXq+jXC5jY2Oj6zX7+/uwLAu6ruPs7Ay6rnNEhoiIZgobNBSozgrZMJWtYSptg3SO4LijPe6/e03TyWQyOD09FeekUqm255lFjsLQ+VluNpttz/dqqJTLZXHvyLKMcrncdg/ys0xERFHHtM0UqM4KmTcldTKZ7KpsAXdX2oapkHnTYkuShJ2dHfE+vdJiK4oCWZahaRokSWo7Z9gyifzg/Sxns1mRMKPZbPZN8T4IP8s0iyqVCprNJhqNBmzbFmstTdOEJEmIx+M4OztDKpWCJEmo1WqwLEtsftiZfpdoHK1WS8TrWIxjCH7ibzdiWq0Wms0mWq1W2JcyNrdCdpdEIgEAIguZpmnQdX2sxoU3LbbbS+1mOfM2oLwVwHw+LxIReM+JQlrsedX5uarVaiiVSqhUKr6kJw46/TjQneK9WCyiWq2iXq+3jTb2aszU6/WxP8th/KxE4/DuzZbP59FsNsX9v7+/j0QigbW1NdG4AW7jfTabFfdDpVIJ6/JpRpyenmJ1dRXLy8tYXV0Vsz7IH2zQRMis3BydFTJXZ2XLWyHrV2kbtkI2T2mx51GlUoFhGF3TEf2spHSmH3enTuq63recXC4nznHXaI26H0wYn+Vxf9ZkMolkMjn2zwqwIUWjs20bR0dH4vH6+jqq1SqA244xx3FwfX0tRmFqtVrbd4+iKG2vJxpVq9XC5uYmrq6uAABXV1fY3NyE4zghX9kMc+bAO++8E/Yl3Nuvf/1r58GDB04sFnMAOLFYzHnw4IHTarXCvrSxZLNZ5/r6eiLvM4r7lnl9fT2R6551YdxzjUbDURRFPK5Wq046ne77uJdRrjudTjuNRkM8lmX5zuurVquO4zhOuVzu+uze9fpOQX6WR/1Zy+Vy28+qqmrb88P+rI1Gw9nd3RX/liTJUVXVUVXVKRaLPV+zu7vrpNNpR1XVtp+vXC47+Xx+qHIpGEHFiWw2Kz5HxWLRaTQaTr1eF8/n8/m2+7Ferw/8jM5CnYL89fLlSwdA13/NZjPsS4ukYe45jtBEhG3buLy8FFPNWq0WLi8vI5utKKw9d+YhLTbdMk2zLfV2PB6/d0+/l3eqY6VSgSzLYj1Lr3JkWRbZ/I6OjsQ+R65R94MJ8rM86s+aTqfFz6ooStc6uGF/1lE3xK3VarBtW+xXxQ1xyU33782Oads2ZFkWSWhevnwZ4hXSLJIkCQ8ePBDrZmKxGB48eMD6g4+YFCAi3Jvj6uoKrVYLsVgMKysrkb45JlEhI+rH70qKt0J+enqKZrOJarUK0zRhGIaY4uLlTtUC0JWq3E0/PiiFeVjG+Vld+Xy+b6r1u37WfhviLi8v91y0Xa1WkUwmAdz+fjunmbkNKcaO+bGzs9O2VtLbENY0DYZhIJFIjNSZcHFxgc3NTfF4e3sb29vbk7pkmgGxWAwnJyfY3NzE5eUlVlZWcHJygoWFhbAvLRIODw9xeHgoHl9cXNz5mkg3aPplK5lFvDmIRrO8vDzyCOawFZXOSnEikcDy8jKA2xGJfvsiybKMer2OWq2GjY0N1Ov1ka4vDOP+rABQKBSgadrYSTT6bYhbKpWQy+W6Rlw0TRNJQ9ysVV7T3GicB+NUUu6jUCggn89DlmWRKMBN9uJqNpvQdb1tTWqz2eyZzt+1urqKk5MTX6+dom99fR0XFxe4ubmBJEmsr42g87vX+73cT6SnnPXLVjKr3Juj2WzixYsXWF9fD/uSiKZW51SnuyopwNcVFfe/fr2unVOmUqmUqCT1y7rnnZrlpob1mtb0yeP8rMDtdDFVVXs2Zob5WcfZEFdVVTEK5u67Q9Nje3u77f5aXV31raxKpQJVVcVnoFariamSrnq9jkwmIz43LtM0kclkfLs2mh+xWAxLS0tszAQg0iM03t64eeHeHOQP5oyfHaqqwjAM8XjSlRTvfjCKokDTNLEuptd+MG5lqtlswrIsPH36tO39JtmgmfTneNSftVQqoVQqiUZHPB5v6xUf5mfttSGuLMuQJKnvhrgARDmlUkmkfh+lXIo+79ROlztSY5omSqUSgNvPq9vgzufzKBQK4vPB1PxE0RLpBg1wG7hs276z55WiJYyGxenpqZjS9+DBA5ycnHAULCJqtRqq1Sosy0KpVBI9s35WUtz04+6c/N3d3Z7rOtz045Ik9V2UPsm9jfz4HI/6s2az2Z6L9oHRftZRN8QFIDbDdf/+XmzQzAdZlvumx+332VNVlVMRiaLM/2Rr/ikWi069Xneur6+ddDrdNwUpUyxGy89+9jPnwYMHDgDnwYMHzs9+9jPfy5y1tNhhi+o9N+p1h5V+vB8/P8dh/KzFYrFveuZR1et1kbqXpsO8xAkiup9h7rkFx5mNXX5KpRLq9XrPHtBkMtk2V5cZSaZXq9XC6upqVza3Fy9e+DoHtdlsioXOncc5xe9uvRb7RmHBe6fNzc2RF/veN2uWO+IwiTWAfn+Ow/hZc7kc8vn8vX8/vZIIULjGud+mQVSvmyiqhrnnIjvlrFar9cxW0gszkkSHu9+Oy7vfjp8Ni1lMix2kcTKSzIppSj/u9+c4jJ+1WCzee7+tcfasIiKi6Ijsqud+2Uoo2sLajMpNi72ysgIATItNkTSrn2NuiEtERINEdoRmULYSiq4w99thzniaBfwcExHRvIlsgwZgWsVZFWaFjGmxadLCyNjHzzEREc2TyE45o9nGzahoFpyenmJ1dRXLy8tYXV1t2418FrVaLTSbTbRarbAvhYiI5ggbNEREPmi1Wtjc3MTV1RUA4OrqCpubm333x4i6sBpvbEQREREbNEREPnAz9rkVbW/GvlkTVuNt3kbAiIioNzZoiIh8EFbGvjCE0XibtxEwIiLqjw0aIiIfzGoK5V7CaLzN0wgYERENxgYNEZFP3Ix9zWYTL168wPr6etiX5IswGm/zNAJGRESDsUFDKJVKME3zXu9h23bbRqdEdGteMvYF3XjzsxHFmEg0Pt4/FAY2aOacZVloNBpQFAUAUCgUoGkaNE3rO3Wj1zmSJEHTNBQKhYCunIimTdCNNz8aUYyJROPj/UOhcebAO++8E/YlTK10Ou00Gg3HcRynXq876XS6699ed50jy7LPV0xRENV7LqrXTZPDmBicqN5vUb3uIPD+IT8Mc89xhGbOWZYFWZYBALVaDZqmAQAURUGtVus6/65zJEniolwiiizGRKLx8f6hsLBBM+e8gaLRaCAejw88/65zUqkUzs7OJnV5RESBYkwkGh/vHwoLGzRzzLbttoxAkiTBsqyBrxnmHCKiKGJMJBof7x8KExs0c6xzKDeTyYidtk3TRCqV6nrNXed4h5uJiKKEMZFofLx/KEzfCPsCKFyyLIteFUVRIMsyNE2DJEkol8vivEQiITKX9DsHYPAhomhjTCQaH+8fCgsbNHNO13UcHx8jm80CAPL5fM/zGo2G+He/c0zTRDqdnvxFEhEFhDGRaHy8fygsnHI257LZLOr1+kSyiBSLxb6BiYgoChgTicbH+4fCwgYNoVgs3vs9bNtm4CGimcCYSDQ+3j8UBk45IwBoy0wSxuuJiKYJYyLR+Hj/UNA4QkNERERERJEV6RGaWq0Gy7IQj8dhWRZ2d3fDviQimjKME0Q0DMYKouiK9AhNPp9HNpsVWTAqlUrIV0RE04ZxgoiGwVhBFF2RbdDUarW2OZaKouDo6Ci8C5phrVYLzWYTrVYr7EshGgnjBPmBMXH2MFYEh/cP+SGyDRrTNBGPx8XjeDwO0zRDvKLZdHp6itXVVSwvL2N1dVXs6EsUBYwTNGmMibOJsSIYvH/IL5Ft0Lx8+TLsS5h5rVYLm5ubuLq6AgBcXV1hc3MTjuOEfGVEw2GcoEliTJxdjBX+4/1DfopsUoDl5eWhN266uLjA5uameLy9vY3t7W2frmx22LaNy8tL8bjVauHy8hK2bWNpaSnEK6Npc3h4iMPDQ/H44uIixKv52ihxAmCsoMEYE+9nWuMEwDpFEHj/0LDGiRWRbdAoitI2VNlsNqEoSs9zV1dXcXJyEtSlzQxJkvDgwQNcXV2h1WohFothZWWF+eGpS+cXuvfLPkyjxAmAsYIGY0y8n2mNEwDrFEHg/UPDGidWRHbKmaqqsCxLPDZNE5lMJsQrmj2xWAwnJydYWVkBAKysrODk5AQLCwshXxnRcBgnaJIYE2cXY4X/eP+QnyI7QgPcplgsFAqQZRkARKpFmpz19XVcXFzg5uYGkiQx8FDkME7QJDEmzi7GCv/x/iG/RLpBo6oqVFUN+zJmXiwW4/xWiizGCZo0xsTZxFgRDN4/5IfITjkjIiIiIiKK9AjNsBqNxlQtPiSadY1GI+xLGAtjBVFwGCeIaBjDxIoFhwnAiYiIiIgoojjljIiIiIiIIosNGiIiIiIiiiw2aIiIiIiIKLLYoCEiIiIioshig4aIiIiIiCKLDRoiIiIiIoosNmiIiIiIiCiy2KAhIiIiIqLIYoOGiIiIiIgiiw0aIiIiIiKKLDZoiIiIiIgostigISIiIiKiyGKDhoiIiIiIIusbYV9AEL7zne8gkUiEfRlEc6PRaODP//zPw76MkTFWEAVnUnHCtm1IkiQem6YJSZIQj8dxdnaGVCoFSZJQq9VgWRbi8Tgsy8Lu7i4A9D3eD+MEUbCGiRVz0aBJJBI4OTkJ+zKI5sbm5mbYlzAWxgqi4Nw3TlQqFZyenqJWq6Fer4vj+/v7qFQqkCQJe3t7UFUVAJDP51GtVgEAhUIBlUoF6XS67/F+GCeIgjVMrOCUMyIiIoqcdDqNXC7XdVzTNDiOg+vr67ZRGO8ojqIoODo66nuciKKFDRoiIiKaKZZlwTRN8dg0TcTjcfE4Ho/DNM2+x4koWtigISIiopli2zZkWYau67BtGy9fvux5Xr/jRBQtc7GGhoiIiOZDNpsV/9Y0DYZhIJFIwLbtrnOXl5d7Hh/k4uKibU7/9vY2tre3x71cIupweHiIw8ND8fji4uLO17BBM6Sbmxt89dVXYV/GnV5//XUsLi6GfRlEcykqcSJojEsUlFqt1rbIHwCazSZ0Xcfp6WnbMUVRoChKz+ODrK6u3jspAGPFZDHGzJbOToJhkgKwQTOEm5sb/MEf/AF+9atfhX0pd3rttdfwwx/+kDd2gEqlElKp1J1fgn6zbRuGYaBYLIZ6HfMqSnEiaIxLjBNBkWW5LVFAvV5HJpOBqqowDEMcN01z4HE/MVZM3izFGMaK8bBBM4SvvvoKv/rVr/Duu+/iN37jN8K+nL5++ctf4smTJ/jqq69m4qaOAsuy0Gg0xBQH27ZRKpWgKIpIFQrcpgJ1ewzL5XJbVp1arQZN03B9fS2OG4YBy7Jg23bX+f3eT5IkaJqGQqFw5z4KNHlRiRNBY1zqjhO5XA5nZ2cAbu9fWZb7HrcsC8lkEqlUCgCg63rblKpe8cMts9frZilO1Go1VKtVWJaFUqkEVVUhyzJM00SpVAJwm2LZTcGcz+dRKBTE7/uu435hrJisWYoxw9YpBsUQy7IA3H6u3UbRoDrIzMQKZw68884793r9ixcvnN/7vd9zXrx4MaEr8kdUrnOWpNNpp9FoiMfZbNZRFMWpVqviWL1ed9LpdNe/XaqqOoqiONfX147jOE61WnWy2az4d+f5d72fLMuT+eHu4b73XFjuc928/3rj76U9TpTLZREfyuWyo6rqwOONRqPrHvfqjB+uu17HODG+ealTRMUs/T6HqVMMiiG7u7uO4zjO9fW1OH5XnWFWYgWznBHdg2VZomcEAIrFoujlcLk9qMDtHge1Wk08ZxgGDMNoSxtarVaRTCYBAKqqdqUQHfR+ACBJ0siLXInIP944kU6nRU+roihoNpsDj7uvNwwDhUKh7X17xY/Ocnu9DmCcIJpGw9Qp+sUKRVFgmiZs20atVhOjM3fVGdxyox4r2KAhuodhbvJGo9GzwmGaJizLahtGBm6z8ri7XtdqNTF8fNf7uVKplBiKJqLw9YsT+Xwee3t7A4/H43HE43Hx2F0f0i9+uPq9zsU4QTR9Rm04eGOFLMtQFAUbGxvY398Xx++qM8xKrGCD5p5M00QymcTS0hISiUTb4kLvcfc/XdcB3M7rrVQqYV02TYBt211rW3qRJKmrUQIA+/v7sCwLuq7j7OxM7Jegqqo47p0b632/w8PDqdj8zbbtnjt1U2+2bUPXdREPNE0Tf8dB8cKlaZoYveulUCi0vb7z/YYpK5lMimtzyxsmVlUqFbFugb7WL04UCgVomta1XqPzuCRJqFarkCQJ6XRa9K72ix8u93XHx8f4e3/v7/XslQ0SY8XwZjFO6Lre1fufy+XEyIHLNM22Y7Varevnm1XD1ilcnbGiUqnAtm3U63WUy2VsbGwA6F8HcfWLMVHDpAD3YNs2NjY2UC6XoaoqbNvG8fFx2/MHBwe+LzCkcAw7DJvJZLC/vw/gNli7w8flclmck0wm2xbquYv3SqUSEolE2/v9w3/4D/Enf/InYnj5u9/9LgqFglg06A5Z91sE2O94v8WHLndI2v2Z0ul09BYNhshdeJnP58Xf3jRN0WC9K16Yptk2DamX3d1d8XdwO1uur6+7zrurrHw+L56zbRvJZFL0/vXDONdbrzhhGAYymUzX77PXcfczIkkSTNMUzw2KH+7rYrGY6J1VFKVtIfFf/+t/XdzPve79QckI+i08duVyOfFZzefz4voZK+42q3FC0zSUy+W2v/3Z2VlXRds7PcptpA+qjM+SUaZ29YoVnX9397361UFc/WKMq3Ma3LRig+Ye3C8G9wtAkqS27DPjGNRCH7X1Tv6TZbnt76LrupgKYlkWstksFEWBLMvQNA2SJLVVRPpxz5VlGfl8HsDtqF6j0cB//a//Ff/4H/9jcc7f/bt/F0dHRyIIudnRTk9PUa1WYZomdnZ2UC6XYZpmz+PAbYA8OzvrW2nVNA2NRqPreDqdRiKRYCXlDrlcDtlsti1GjJKWc2dnB3t7e+KLKSiSJEFV1bY52TQab5wolUoolUqiFzQej6NarfY9LkkSdnZ2ROXkrvjhxglJkvCP/tE/wre+9S2Ypol/+S//JVZXV1EsFlGpVPCv/tW/EpWUfve+qqpd5VUqFdGj644kePd8cUeNVFVFpVJBPp8XaV8ZK+42q3Fia2urbYTOrSTLsoxKpSIaRkdHR+Iz535nzcsIDTBcnaJfrMhms2LUq9lsit9jvzqIN1ZsbGzgG9/4Bt54442uez7IBs19UkWzQXMPbkNG13Xkcrm+c5mHZZomDMNo+3IY5jkKj67rOD4+Fl8+/SobbqOkH3fNjKvX39ltTFiW1VWOd169O2TcaxHgoMWBxWKx75SQSqUi9newLKurV9btWWKDu79arTZ2Pv9arSa+/INWq9VwfHzc9Rn1KpVKqNfrkdmvIGjeONFZWXX1Oy7L8p2NGO/fxo0Tsizjb/2tv4X/+T//Z9f5f+2v/TW88cYb4nG/e98dlV1eXhaNEEVRUCwWuxYee6/X/ZweHR117enCWDHYrMYJt4OuVquJxq77XVStVpFOp2HbNmzbjsRogF+GqVP0ixUA+n52etVBvB2UqqqKc2zbFrM+4vG4aGz2SxXda9bHXenmgd7bU9xnJJdraO7p/PwcwO0fa2lpqWsOqWEYSCaT4r95GTqdF9lsFvV6PdAMIIPKKhaLyOfzfRcB3rU4sJ/T01M0m03x/t61YkB0Fg2Gxb3v7/qi7hcvcrncnY3iUQ2KTTs7O2L+fD6fx/n5ed9rr1QqyGazkZ13HYQw4gTQP1YYhoE/+IM/GPjafguF+y089nIrM72SFjBW9DfLcQK4HaFzK+jVahVbW1tQVVVM1XcbO/MsjFhhGEZbh4ZhGDg6OgLwdZ2iUqlA13XU63Xs7e21daK6sz7y+Tx2dnbE+6iqimq1KkaPvGq1mmjIGIbR9rp0Os0RmjB4h+/cP3ij0RA3tXeOaT+VSqVtSsHS0hLi8TgajcbA52g6uL2VQRjUs/kXf/EX4sus3yLAuxYH9pNIJLC8vAzgtoeWFZLRuPHgrqH7XvGiUqmIKQOTTAQxKDa58+ZrtRpyudzA3vR0Ot02x596CzJOAP1jxX/8j/8Rjx49uvN7yZ1WBtz+jd3edO/CY8uysLGx0dUrL8sy6vU6arVaz+ept1mOE8DtWg53+liz2RQ98vF4HJZloVqtztX0sn6CjhW9UkXncrm2OoX3M+BNFT1o1kevEV5X5/YUnSPE44zkcoRmgtLptFioPerrrq+vUa/Xoaoqrq+vRYNl0HM0PYKaPjFo0eDf/Jt/U1xHJpPB6ekpgPZFgP2O3yWVSonKTa8v26gsGgyTqqpj9Tqdnp6KHm93PvWgDEaT5O683mtvAq9B0xXpa0FOs+qXjOCdd97Bv/7X//rO17v7Wbj/dqeW9Vt47H2dy62oejFWDDbLccKtCFcqlbaRGHcKGkdovhZkrBimTuHlTRXdb9bHXamg79qeYpyRXI7Q3EOtVoNpmshms5AkSfxRxr0h3d6KUZ+j+XKfRASDEhT0eh930aCiKCI9J9A9r5eVlLsVi0WR6tQdfnfnDw9a9OudQlKr1WAYRqA93u51u3GuF3fef6lUundiFJqcYZIRAL3v/X7JCPotPHZjhbvWrtlswrIsPH36tO2aGCsGm+U4Adw2Xvb399uuV9d1GIbBz0UI7psqetBskF4jvC53zY6u64jH4xP527NBcw9ur/Xa2hqA2y+Pp0+fjt3wGLT4c5iFoTQf7puIoN/xXu/jHQ30pvr0chMR0GCyLOP8/Bw7OzttUwMPDg5CvrLBZFnG1tbWwMwzbg/r1tZWwFdHgwyTjADofe8P+s7p9TlwY4UkSX0/J4wVd5vlOAF8Pe2sc4TGNM2u1+VyOdFRnEwmsbe3x8/PhN03VXS/lNB3pYIGBm9PMU7HBxs09yBJEvL5fN8K4qCpYZw2RuNye0inJVNQsVhkdqshDUrbPUxMUFV16F5XRVHgOM7IZfV6/7v+vvz7TyfGimia1TgB3PbW9yqv1zF+VoJxn1TR/WZ99BvhdUdygd7bU7jYoCGKoFarJYJJLDbcsragFw32Y9v2xLPqEFFvjBVEdJdR48R9U0X3uq/7jfB6G8n9tiEZdySXSQGIQnR6eorV1VUsLy9jdXVVLNgfxjT0uHJtF1EwGCuI6C7jxImw0sr346aKHhUbNEQhabVa2NzcxNXVFQDg6uoKm5ubfYf/iWg+MVYQ0V3uEyemZXrffUZy2aAhColt27i8vESr1QJwG4wuLy+nppeEiKYDYwUR3eW+cWIaRlDvM5IbyhqazgWKbhYLN2e9m0lp1ONEUSJJEh48eICrqyu0Wi3EYjGsrKxMRVAhounBWEFEd5n3OBHoCE2lUoFhGNjY2Gg7ns/nkc1mxSKgSqUy1nGiKInFYjg5OcHKygoAYGVlBScnJ1hYWAik/FarhWazKXpziGg6hRkrGCeIoiHsOkXYAh2hSafTUBRFpH0DbkdbvK1HRVFQLBa7hp3uOh5EbvJf/vKXvpdxH9N+fdRtfX0dFxcXuLm5gSRJgQWe09NTbG5u4vLyEg8ePMDJyQnW19cDKXvW8T5sx9/HZIQRKxgn/MV7YzL4e/xaWHWKaRB62mbTNBGPx8XjeDwO0zRHPu6n119/Ha+99hqePHniazmT8Nprr+H1118P+zJoBLFYDEtLS4GV12/h4IsXL+Yq+E1alOJE0BiXJiPIWME44R/GisljjPla0HUK1zhp5Scp9AbNy5cvJ3J8kIuLC2xuborH29vb2N7eHvr1i4uL+OEPf4ivvvpq5LKD9vrrr2NxcTHsy6Ap5i4cdHkXDo4bBA8PD3F4eCgeX1xc3Ps6oyZKcSJojEvR40ecoFuMFZPHGBOuaRjNDb1Bs7y83DMDw6jHB1ldXcXJycl4F/hXFhcXebPQTPBj4WBnJ4G3A2GeME7QrJj3BcZ+Y6ygWTEto7mhp21WFAXNZlM8bjabUBRl5ONENJx5XzhIRHdjnCCiYUxLWvnQR2hUVYVhGOKxaZrIZDIjHyei4U37wsFKpYJms4lGowHbtsWmX0zxThScaY8TRBS+aRnNDXSEplaroVgswrIslEolWJYF4DYNc6FQEOmX3Yxlox4nouG5CwenrZJiWRYsy0I2m0U+n0ez2UShUADAFO9EQZvWOEFE02FaRnMDHaFRVRWqqiKfz/c83u/8YY8TUfTZto2joyMxmrK+vo5qtQpFUSKT4p2IiGheTMNobuhTzoiIvBRFQb1eF48bjQYURYlEinciIqJ5FFa6aFF+aCUTEd3Btm3UajXs7e0FkuKdiIiIoocjNEQ0tXZ2dlCtViFJUiAp3u+7ZxUR9cf9qojIL2zQENFUKhQKyOfzkGUZlmVBURScnp6K572p3Ec5Psgk9qwiot64XxUR+YVTzoho6lQqFaiqClmWAdxmSFRVVWRGBNpTuY9ynIiIiGYLR2iIaKpYlgVd19uOuZkR3ZTtbkOnM5X7sMeJiIhodrBBQ0RTRZZlOI7T8zmmeCciIqJOnHJGRERERESRxQYNERERERFFFhs0REREREQUWVxDQ0RERJFl2zYkSRKPa7UaLMtCPB6HZVnY3d0d6zgRRQdHaIiIiChyKpUKDMPAxsZG2/F8Po9sNiuyGlYqlbGOE1F0sEFDREREkZNOp5HL5dqO1Wq1ttEaRVFwdHQ08nEiihY2aIiIiGgmmKaJeDwuHsfjcZimOfJxIooWrqEhIiKimfDy5cuJHB/k4uICm5ub4vH29ja2t7dHfh8i6u3w8BCHh4fi8cXFxZ2vYYOGiIiIZsLy8jJs27738UFWV1dxcnIy3gUS0Z06Owm8HQj9DD3l7MmTJ+NdFREREVEAFEVBs9kUj5vNJhRFGfk4EUXL0A2an/70p8hkMl3HDw4O8NFHH030ooiIiIhGpaoqLMsSj03TRCaTGfk4EUXL0A2a8/NzbG1tIZPJ4NWrV2LEZmdnB4uLi/jkk098u0giIiIir1qthmKxCMuyUCqVRMMkn8+jUCiI9MtuOuZRjxNRdAy9hsayLLz33nvY2NjA7u4urq+v8e677wK4bdR8//vfxw9+8APfLpSIiIjIpaoqVFVFPp/vebzf+cMeJ6LoGLpBs7S0BACQJAkff/wx3nzzTXz55Zf45je/CQBtQ7ZERERERERBGHrKWWfvRTqdxu///u+Lx47jTO6qiIiIiIiIhjB0g8YwDHzwwQfi8d7eHvb397G3t4dXr1617bRLREREREQUhKEbNGtra8hms/j8888BAIuLiwCAR48eYXd3F8vLy75cIBERERERUT9DN2gA4O2338Zbb73VdmxxcRG7u7uigUNERERERBSUkRo0/ciy3JVlhIiIiIiIyG8TadAAt1PSiIiIiIiIgjSxBo270SYREREREVFQht6HxuvTTz/F/v4+bm5uANymbPZutElERERERBSEsRo01WoV5XIZ8XgcwG2D5sc//vFEL4yIiIiIiOguYzVoNE3rWjPzu7/7uxO5ICIiIiIiomGN1aC5vr7GRx99BFmWIUkSbNvG0dERjo6OJn19REREREREfY3coLm5ucGPf/xjqKqKX/ziF+K4aZoTvTAiIiIiIqK7jNygWVxcRLFYxMbGRtvx58+fT+yiiIiIiIiIhjFW2ubOxgwAJBKJe18MERERERHRKIYeoXn27NnA58vlMn7yk5/c+4KIiIiIiIiGNXSDJpvNQtM0OI7T8/l6vT6xi5oHpVIJqVQKiqIEUp5t2zAMA8ViMZDyiIiIiIiCMHSDpte6GS+uoRmeZVloNBrIZrPiWK1Wg6ZpuL6+hiRJsCwLyWQSqVQKAKDrujjfsiwYhgEAyGQySKfTXWXkcjk0m00AQD6fhyzL0DQNhUIBu7u7fv+IRERERESBGLpBM6gxAwBvv/32vS9mXhiGgXw+33Ysn893jdaoqopyudz1ek3T0Gg0+r6/ZVnQdR2qqqJSqSCfz6NYLCKdTiORSLBBQ0REREQzY6ykAHQ/lmVBlmXx2DAMGIaBeDzedZ5hGCgUCuJYpVKBLMvI5XLQNK1numxZlqGqKgDg6OgImqaJ59x9g4iIiIiIZgEbNCHwNihM04RlWaIB4orH44jH49jb2wNwO4UMAE5PT9FsNlEsFpHP58XUs07ulLXO906lUjg7O5vwT0REREREFI6JNGiePHmC9fX1SbzVzLNtG5Ikicf7+/tiitjZ2Rl0XRfnVKtVSJKEdDqNWq0G4DY9diaTAQAoitK3cSLLMur1OvL5/J3TBf1UKpUC3XTVtm3R+CMiIiKi2XevBs2zZ8+wvr6OH/zgB4FWWqOsc8pXuVxGvV5HuVyGLMsol8uQJAmmaYrzTNMU62tSqRSq1SqA7qlrLu/fIh6Pw7Is8bjfa/zgJj/wrg2q1WpYWFho+x3Yto1CoSAaba5cLodkMilGmnoxDAO6rkPTNNEQdJMfEEUdOwSIiIjuNlaD5tmzZ0ilUkin08hkMmg2m33TOVM3WZbvXMciSRJ2dnagaRqKxSIODg4A3I7KaJoGTdOQy+Xakga4m5u6a2x0XcfOzg6ePn0qzgmyQWMYRlflqFfyA8MwcHR01HasUqlA13XU63Xs7e31rGTVajXYto1yuQzDMLCzswMASKfTTE9NkdfZIeCum+tcO9fveKFQEMd7xRvLsrC0tCTOKZVK7BAgIqJIGjrLGXDbkNnd3cX5+TkePXqEDz/8UDy3sLAw8YubVbqu4/j4uC1tM9C+l487WtPL7u5uz0xlbuYzSZJ6VuhN0+yZ4tkv/ZIfdGZ4KxaLXQ0W73UqiiJSUHtVq1Ukk0kAtxnhvO/hjoR5p/cR3UfQe0f9+3//7/E3/sbfAHDbwHenodq2DV3XUa1W+x43TROnp6fi3zs7Oz3jSa9MisyGSEREUTP0CM3Dhw+h6zp+53d+By9fvmxrzEyCuzjetm3R8w7c9sKXSiVUKpW2XsN+x0cVxpSOer2Oer0eeLYxN5FAUIZJfjCMfD4vkiN4aZomGoG1Wq1tWhqTH9AkDTtaAvSeVnnXaEmv1/6f//N/8N5776FQKEBRFDENtVarievod9zd18o9p3M6p/fn6sykCISTDZHT64iIaFxDN2iKxSIePXo0VoV0GPv7+0gkElhbW4NpmqJnPZ/PI5vNih77SqUy8PgoelVSdF2HruuiclwqlcQ6jmQy2TYS1bl+o5fOiow7pcOdHhYU27YDb8wMk/zgLu7vr9fIkqqq4j3dNUhEfvBOn/SOirjTHb06p1V6R0vy+byYGtmL97W2bYvpk7IsQ1EUbGxsYH9/XzTw+x1vNBpdaeA79cukCATfITDK9Lpe6+p6xe5evA1GTq8jIpohzogsy3JKpZLz/PnztuOxWGzUt2pTLBa7jlWrVSedTnc97ne8n3feeafn8XQ67TQaDcdxHKfRaDjVatVxHMcpl8tONpvteY3udVarVXFOv/Lr9bo47v234ziOLMt9r3dW9PsZFUVxrq+v245ls1nx+3ft7u469Xp9qLKKxaKTz+fFY1VVxd+Wgtfvnpt2/a5bURTx70aj4aiq6lxfXzvlctnZ3d0Vz+3u7jrValU87ziOk8/n2+KbJEk9y/C+9osvvhBlKori/OEf/qGIN41GQzznjVXe47u7u233Q78yvT+T937tdT/6yRuLvb/T6+trR1VVcdwbo93jw8Rul6qqXfFnHmLxtJq1OEFE/hjmnhs5KcDa2hp2dnbgOA4eP36Mzz//fGKNK8uy2nrjTNNs62WMx+MwTbPv8XHKc3v1B21GCUAsPnfXvXSu3+hV/qBpH/OwweUwyQ+A2zVFtVoN+XwepVIJwO3IWKlUws7ODpLJZNvfwzu6pWkadF1Ho9Fom/MfZPIDmn3ez3G/UZF+0yqHGS3pfO3i4qIoM5VKdcVZ97nOtWXu8Uwmg9PTU/HeqVSqZ5m9MikCwd8/3vL6TaNLp9Pi9+NdV3dX7Hb128B4HmIxEdGsGykpgNfbb7+Nt99+G8+fP8fjx48nkuXMtm3Isgxd13FwcICXL1/2PK/f8XHK83KnLwHoqpR0ZuzSNE0spu1cv+FqNBp9v1zdKR1+TeGbBsMkPwDQc7FyNpvtep3LTX4AQKSw9go6+QH5o3PaonufuanI3QbsqMfvex2VSkWshbMsCxsbG6jX6z2nVbpp2AdNgwJ6T8n8O3/n74gY9c//+T9HuVyGpmloNpvinslms2J6lve4oiiQZRmapkGSpK5siI1GQ2RSdMvwnhN0g6ZfgxFAW5ZGV+e6ukGxG2hvMHZOvZ2HWExENOvGbtC43IbNfb8MvJVXTdNgGAYSiUTPnrPl5eWRetQuLi6wubkpHm9vb+Of/bN/1pUBy92MslariUqK6/j4uC1zmPvFqOs64vF4zy//YSoys8ytbAWdbaxYLDJtc8AODw9xeHgoHl9cXIz9XpVKBaenp6jVam33YD6fFw3YQqGASqWCdDo98vFRdfbg9xsV8TYIksmkaMxkMhns7+8D6D9a0uu1x8fHOD4+Fo2Lfp/pfsf7rZlzOwT6ZVIMukNg2Aajq9e6ukGxG+jdYHT/PkREFH332ljT6+233x77td6pWa5ms9mVrtc91u94P6urqzg5ORH/bW9vd1VSBm1GaZpmzwaLuyg4mUz2zJYzaNrHuD2grVYLzWYTrVZr5NeGIeiGhZv8gBmTgrW9vd12j62uro79Xul0uut3WavV2iqfiqLg6Oho5OPj8k6fdDtf3OmO/dKre8t2R0v29/d77h3VSzabRb1ex89//vNAR0uCzoY4bIMRuB0pV1W1rTEzKHa7+m1gDHB6KhHRLLhXg+aTTz6ZyEW4G0G66vU6MpmMyGLlMk1z4PFxynW/LAdtRtlsNnt+4fVbv+FWUgZVZMb5Ej09PcXq6iqWl5exuroqGkvTLsheUEmS0Gw278xe5x7vlzGp13GvznOYMclfo66nm9Q6O5c7fdJVLBZRrVZRr9d7dqbU6/W2z707WtQ5KuCdPtnrtblcDv/iX/yLsa97VEFnQ3QN02Dst65uUOweJpskGzQ0C9iJR3PvPlkHHj16dJ+XtymXyyKLmDc7T7VadfL5vFMul4c63ku/7AjerGVBqtfrbZmRhvHrX//aefDggROLxRwATiwWcx48eOC0Wi2frjK6hsle1y9jUr/jXoPOYcakW/fNAuTN2OU4t1m7vNmr6vW6I8vyyMfvoiiK884774j//uiP/kg8l81mu7Lz+W1Qxq5ZEqVYTOP7oz/6o7b7y3uPR8m0ZTlrNBptn+NsNuuk0+m270L3uKqqjqqqbRlEs9msoyiKoyhK3+ygu7u7TjqdbsvgOEz9i2gShrnnpqZB46dBv4ioVFJevnzpAOj6r9ls+nCF0dbvSzKdTjvlcrnreGfl+a7jg87plZJ6Hk26QZPP53s2UEY9ft/rDvJve319PVefpajEYpqcaWsYDGvarnvYTrxR0qF7Ddqmgp14FARf0jbPmrDWeIxKkiQ8ePAAsdjtnywWi+HBgwdc1NpDr+x17tSwXskrOjMm3XV80DlBb0g4L0ZdTzfqOrthBT19cp7u76jEYooGN7Odm/7b/V6o1WoolUqoVCptU4T7HY+CYbagGDUdutegbSqY9pymxdw3aIBoVFJisRhOTk6wsrICAFhZWcHJyQkWFhYmfIXTYdzkB70yqrkZkPL5vEgF6+qVMWnQ8VHPockYdT3dpNbZUbCiEIspGvb395FIJLC2tgbTNMXfOp/PI5vNirhdqVQGHo+CYTrx+u2f5dWvE0/TNJE5sHObCnbi0bRggyZC1tfXcXFxgWaziRcvXmB9fT3sS/LFfZIfjJK9rlfGpEHHhzmHC4zvr1aroVgswrIslEol8TfL5/Mi/TIA8bsf9fg8iFo2RKJJ0zQNjuPg+vq6bW+qILIiBrlA37ZtfPOb32xboN+rE8+bDr1cLg/duQd83aHkJungdxxNo3vvQ0PBisViWFpaCvsyfNNqtbC5uYmrqysAwNXVFTY3N/HixYuhR6PcjEmSJIkMSM1mE5ZliQxIbsakWq0G4LaxU61W+x4Hvt6QcNA5bNDcn6qqPTdAdI/3O3/Y47Pu9PQUm5ubuLy8xIMHD3BycjKznR9Eg7hTztzpVUFkRbQsC41GQ2Trc79/gNtOFvf7YdTjnWUYhgHgdnuIL7/8UmTZVFVV/LzeTry70qFnMpmB03Ld77hSqdSWPZDfeTQt7tWg8QYAokmwbRuXl5ficavVwuXlJWzbHroh56bYzWazkCSp59z8bDbbtpnrXceBr1Ps9jsn6A0JKRparZZoYLtr4Pws674dAkSzwrZtyLIMXddxcHCAly9f9jyv3/FxGIYhOmPcUQ1VVVGpVJDP58Xo8yjHO2ma1pbyvVgsQlVVJJNJsaF1Zyeee1zTNDSbza506IM68dwy3U5Cb2cTGzQ0Le7VoPnwww8ndR2RE2QlZZ64yQ+urq7QarUQi8WwsrIy0lx3N3D3Wk/jp2KxGPjCZppuQY+WTKJDIGoYi6kXb6eTpmkwDAOJRKLnAvbl5eWRFrZfXFxgc3NTPN7e3sb29jaA7gX67r+Pjo7EOr5Rj3tVKhUx88CyLOTzedGJ537f9fseGqVzD2jfJ8tt5HixE4/8cnh4iMPDQ/H44uLiztcw+o8hqhtcRsGkkh8wYxKFrd9oieM4vpU5b9kQGYupl1qtJrJ7udxMh5PIiri6uoqTkxPxn9uYAYbPsjnqcdfp6SmazSaKxSLy+TwMw0A2m0W9Xsdv/dZvBbpA372GoHET0dm3vb3ddo+trq7e+Ro2aEYURiVl3kwq+QEzJlGY3NESd2G+d7TEL/OUDZGxmPpxRzBc9Xo9kKyIo2TZHPW4K5FIiGtTFEU0YOalE89do+Q2OHO5HJLJpGgEutzpdZqmtTV+DMOAruvQNK1vLHYTJLjnSJIk1ijR9GKDZkRhVFLmkZv8YBYrYjQfwhotmZdsiJOKxeztnT3utC13fUgikQgkK+KwWTZHPe6VSqX6JqG5uLgIbD1LWJ14hmGI+6dSqUDXddTrdezt7bUdlyQJ1WoV5XJZJFBw9yNyj+3s7HS9v2maOD09RbVaRT6fF+ek02lOKZ9yzHI2okms8SCi2eeOlrhraIIcLZn1bIjAZGJxZ0Yq4OvpStfX15AkSUwBSqVSAG6Tjrjnd2ab6lUJNgxDZNsql8ttvb1uOmGavH4NEr+zIg6TZXPU48DXC/QVRRGjBwDE4n5gPhboe39G79/YO21QURQUi8WuTUQ7Nwjt1bHgna6oKIpIlgB83WBlfW86sUEzojArKUQULe5oyc3NDSRJYpyYoEnEYm9GKlc+n+9aP6GqalvF0dWZbaqTt0e4VqthZ2cH5XIZ6XQaiUSCDZoZNEyWzVGPA+0L9Hd3d7s+O/OyQL/fCKx3U1DvJqIARMNQ0zRxH3duEOpqNBpd669c7iai87gdQBRwytkY5mVKxzzihoQ0aZw+6Z/7xuLOHm3DMGAYRteWBO5IjHcOvTfbVOc8fVdnj7D3nM7pSTQb3AX6Qf9tw1qgH6R+oyOdm4L220R0mA1C3VFZih42aMbESsrsYcYkoui5TyzuXO/QK7NUPB5HPB4Xvb/uNJVe2aY6aZqGer0OoLtH2O3tpdkTpQX6UerE69UJYBgGVFVtG50atImou64mmUz2nHKWyWTEd79pmmKqKTAfU/qijFPOiMANCYnmTWdv7/7+vui9PTs7E7247uJi4HbOvjsdJZFIYHl5GUB7tikvVVXFPiHxeJyVoTkSdJbNcQS9T9YkeNco9dsUtN8mokD/DUK9a5RkWRbnzdsapShjg4YI87khIdE86+zt9VZcksmkaMyYpglZlsW/3fU1qVQKhmFgd3d3YEXHbQy52bZcrBxRmKLaieddozRoU9B+o2S9NggF2tco9Rrtmpc1SlHGBg0RmL2OaB55e3v7kSQJOzs7ovHjNnwGZZtye3uB/j3CbNBQmKLaieeOvgSdbaxYLDJt85Rjg4YIzF5HNI+8vb1e7roX4LbR0yvDGdA72xTQ3tvbq0eYvb0Utih34rkpmYMS1iaiNBomBSD6K8xeRzRfmJGKpkXQi/PdTryVlRUAiFwnXtBrlKLQ0Jt3bNAQeTB7HdF8iVJGKppNYWXYZCcezRI2aIiIaK6xt5fC0m9xvuM4gZTPTrzhlEqlnntN+cW27Z5ppak/NmiIiIiIQuAuznenmnkX59N0sCxLpHQGbv9mhUJBpIt25XI56LoOXdfb9pwyDAO6rkPTtL5/VzdlvK7rqFQqkCQJmqa1beZLg7FBQ0RERBQCd3F+LHZbHYvFYnjw4AFH8Xww7jolwzDaRksMw8DR0VHbOW6DpFwuI5PJiGmltVoNtm2jXC7DMAzs7Oz0LEPTNJTLZZTLZZEwJJ1OM7PaCNigISIiIgpB1BfnR8V91il1plgvFotIpVJt58iyDFVVAQBHR0cinXu1WkUymQRwu9Fur2lrlUoFsiyLzUC953Tul0X9sUFDREQ0gqAzUtFs4+J8f913ndKwDQrLspBMJmFZlmjcaJom0sDXarW2qWiu09NTNJtNkf3QMAzxXCqVwtnZ2VDlzzs2aIiIiIYUVkYqmm1cnO+f+6xTGmUDT1mWUa/Xkc/nsbGxAeB2VMY7Ha3XZrqJRAKZTAbA7Ya9bMCMhw0aIiKiIYSdkYqIRnefdUrDTvnyThOLx+NtIzHVahXlchnJZLJn5rJUKiU24O2c3tb5mPr7RtgXQEREFAVuT6/L29O7tLQU4pURUT/uOqXNzU1cXl6OvE5JluW2kRpd12GaJizLgmVZyGazYg1Ms9mEZVl4+vSpeL2maZAkCbIst+1BlUgkRPY0TdPEuptyuSzOYYNmeGzQEBERDcHt6b26ukKr1UIsFsPKygozUhFNOXed0s3NDSRJGmlqn67rOD4+RjabBdDe4HBJktQ3I5k7+tKp0WiIf+/u7mJ3d7ftedM0RcYzuhunnBEREQ2BGamIomvcdUrZbBb1ej3wbGNukgAaDkdoiIiIhnSfnl4iiqZisRhog8a2bTZmRsQGDRER0Qjcnl4imh9BTi3lNNbRccoZERERERFFFhs0hFKp1HP32lHYtt0zHSERzQbGCSIaBmMFhYENmjlnWZZIGwgAhUJBpA/sN1+01zmSJEHTNBQKhYCunIiCwjhBRMNgrKCwsEEz5wzDEL0gpmni9PQU1WoV+XweOzs7XecPOiedTvdNW0hE0cU4QUTDYKxo12q10Gw20Wq1wr6UmccGzZzzbtpUq9XExk6KoqBWq3Wdf9c5w+6qS0TRwThBRMNgrPja6ekpVldXsby8jNXVVZyenoZ9STONDZo55w0UjUYD8Xh84Pl3nZNKpXB2djapyyOiKcA4QUTDYKy41Wq1sLm5iaurKwDA1dUVNjc34ThOyFc2u9igmWO2bbelBpQkCZZlDXzNMOcQ0exgnCCiYTBWfM22bVxeXoqpZq1WC5eXl5EdbYoCNmjmWOdQbiaTEUOipmkilUp1veauc7zDzUQUfYwTRDQMxoqvSZKEBw8eIBa7rWbHYjE8ePCA+8v4iBtrzjlZlkWviqIokGUZmqZBkiSUy2VxXiKREJlL+p0DRDf4EFF/jBNENAzGiluxWAwnJyfY3NzE5eUlVlZWcHJygoWFhbAvbWaxQTPndF3H8fExstksACCfz/c8r9FoiH/3O8c0TaTT6clfJBGFinGCiIbBWPG19fV1XFxc4ObmBpIksTHjM045m3PZbBb1en0i8zqLxWLfwERE0cU4QUTDYKxoF4vFsLS0xMZMANigoYnkebdtO/KBh4j6Y5wgomEwVlAYOOWMAODeC9W40I1o9jFOENEwGCsoaByhISIiIiKiyIr0CE2tVoNlWYjH47AsC7u7u2FfEhFNGcYJIhoGYwVRdEV6hCafzyObzYosGJVKJeQrIqJpwzhBRMNgrCCKrsg2aGq1WtscS0VRcHR0FN4FzbBWq4Vmsyl2vCWKCsaJYDFWUFQxVgSHcYL8ENkGjWmaiMfj4nE8HodpmiFe0Ww6PT3F6uoqlpeXsbq6Knb0JYoCxongMFZQlDFWBINxgvwS2QbNy5cvw76EmddqtbC5uYmrqysAwNXVFTY3N+E4TshXRjQcxolgMFZQ1DFW+I9xgvwU2aQAy8vLQ2/cdHFxgc3NTfF4e3sb29vbPl3Z7LBtG5eXl+Jxq9XC5eUlbNvG0tJSiFdG0+bw8BCHh4fi8cXFRYhX87VR4gTAWDEuxgoaxrTGCYB1iiAwTtCwxokVkW3QKIrSNlTZbDahKErPc1dXV3FychLUpc0MSZLw4MEDXF1dodVqIRaLYWVlhfnhqUvnF7r3yz5Mo8QJgLFiXIwVNIxpjRMA6xRBYJygYY0TKyI75UxVVViWJR6bpolMJhPoNXhbj7NYbiwWw8nJCVZWVgAAKysrODk5wcLCQiDlz/rvN+wywyw3KNMQJ4DZ/0zNY6yYp3t21uMEMB2xYtY/U2HHCWB+7p9Z/yz1EtkGDXCbYrFQKIjUim6qxaDMwwdmfX0dFxcX+Kf/9J/ixYsXWF9fD6zsefj9hllmmOUGKew4AczHZ2reYsU83bPzECeA8GPFPHymwowTwPzcP/PwWeoU2SlnwG2PiqqqYV/GzIvFYnjttdcC7UUhmhTGieAwVlCUMVYEg3GC/BDpERoiIiIiIppvC84c5Mv7zne+g0QiMfH3vbi4wOrq6sTfdxrLnZcywyp31n7WRqOBP//zP5/4+/ptlmLFrH2mWGZ45TJOtJulOBFWufxZZ69MP8sdJlbMRYOGiIiIiIhmE6ecERERERFRZLFBEzGvXr0K+xKIKAIYK4joLowTNCvYoJlyr169Ev/d3NxA13V8+eWXDEIR9/z5czx8+BB7e3t49eoVUqkUvv3tb+PP/uzPfCvz/PwcH3zwAf74j/9YlLm8vOxrmRQcxorZE0acABgrZhnjxGxinQKAQ0M7ODgQ/87lcs6bb77pbG1tOTc3N76VubCw4KRSKefhw4eOpmnO0tKSk0wmnVQq5VuZjuM45+fnjqZpTiqVcj799FNx/OHDh76WOy8ePnzo2LbtlMtlJ5VKOZZlOY1Gw9e/69bWllOpVBxd151UKuXUajWnXq/7/lmaN2HECccJJ1YwTvgrjDjhOIwVQWGdgrFiUlincByO0IygXC4DAB4/foxEIoHPPvsMuq5D13Xfyvzss8+wtLQEwzDw2WefQVVVnJ2d4fT01LcyAUDXdeRyORwfH+Ozzz7DkydPAAAvX770tdx54TgOFhcXoWkaHMfB2toaZFnG0tKSr2W+9957OD4+huM42NjYgKIoWFtb863MeRRGnADCiRWME/4KI0645TJW+I91CsaKSWGdglPOxtJsNpHL5bC2toZ0Og3Hx0Rxqqris88+w2effYZMJoPz83PfyvJyP6hra2v4+OOP0Wg08Pnnn/u6EdarV6+wtbWFeDyO5eVlLC8v4/vf/34oQ+GffPKJr+8vyzL29vawtbUFVVXxwQcf4IMPPoAsy76VGY/H8ezZMwBff5G6x2nygowTQDixIow4AcxPrAgjTgCMFUFjncIf0xQngNmLFVMXJ4IeEoqypaUlZ2try0kmk2KY1LZtR9f1QMqv1+uBlZXL5ZynT5+2HTMMw4nH476V+fDhQ8c0zbZj9Xrd2dra8q3Mp0+f9vwviGHwSqXiWJYl/l0qlXwv8/z8vO2xaZq+T4WaN2HHCccJLlaEESccZ75iRRhxwnEYK4IQdqxgncIf8xQrpilOcB+aEXh7MuLxOBYXF/H8+XNIkjSTQ/Gffvop3nvvvbZjBwcH2NnZ8aW8hw8f4rPPPhv6+CQ8fvwYx8fHSCaTbcefPn2Kn//8576UGZYvvvgC2WwW19fX2Nvbw7vvvgsA+P73v48//dM/DfnqZgfjhL9xAmCs8BtjRTAYK2avTgHMT6yYujgRSjNqxngX9k1aWAvpnj9/LlrZ5+fnTi6Xcx4/fuxrmYVCwdna2nIODg6cTz/91Hn8+LGztbXlPHr0yLcybdt2crlc1/Ege9O9/PwsuZ8hy7KcXC4nPk/JZNK3Mulrfv5tHSecWBFGnHAcxgq/P0uMFeFinWIywogTjjM/sWLa4gRHaEbgzhXslM/nfWuNrq+v49GjR1AUBfl8Hg8fPsS7776LVCqFs7MzX8r85JNPcHx8jOvraxwcHGB3dxe5XA7VahWxWAz/+T//Z1/KBW5TD9ZqNTQaDUiShEwmg7ffftu38sISxmep8zPz+PFjaJqGnZ0d3xeEzpMw/rZA8LEizDgBzEesCOuzxFgRDNYpWKeYFNYpgG8EXmKE1ev1nsOIlmX5VqbzVwvpAODjjz/G48ePfV9IVywWxYfx29/+NsrlMt566y289957SKVSvpT55ptvQtM05HI5fPjhh76UMajMt956K5AyXWF8llKpFJ49e4bf/u3fBgB8+OGHePToka9lzqMw/rZA8LEijDgBzFesCOuzxFgRDNYpZqtO0VnurMeKqYsToYwLRVQYw4hhLKTz/jyapnVdj18qlYoYCn/8+HEgC8vCKNNxwhuSrlQqXceCWmQ8L8L62wYdK8KKE44zP7EizKkrjBX+Y51i9uoUYZXLOoXjsEETAUF/YN5//33x785sIN7n/GLbtlMqlUTmlydPnsxkmUSTFmSsCDtOOA5jBdE4WKcI5p5lrAgWGzTUZWFhwYnFYk4sFuv576CYpunouj6TZd7c3Di6rjtLS0tOPB534vG48/Dhw0DTHXoXhBKNalrihOPMbqyYhjjhOIwVdD/TEivCiBNBlTsNsSLsOMGNNe/B3el21spstVr49a9/jV//+tc9/+2nL774Ant7e/j2t7+Nra0tfO9730Oz2Zy5MnVdx97eHprNJl6+fImXL19if3/f11S3nTgfPhhhxIkgyg0zTgDzESumIU4AjBVBYZ1i8sKIE2GUOw2xIuw4wQbNPYTxxwuizC+++KLn8VevXuHzzz/3pcxPPvkEqVQKyWQStm3j+PgYP//5z/GjH/0Ii4uLM1Omy3GcrkwriqLg+vra13IpeGEFeb/LDSNOAPMVKxgn5gvrFJMT1j3LWBEeZjmjLrlcDuVyGW+88Ubb8TfeeAOGYfiSAvCzzz5DPp/HxsbGxN97msp0aZqGTCYDTdMQj8dhWRZOT0+7MpRM2uPHj7GwsADHcVCr1QDcBsKlpSX84Ac/8LVsmi1hxAlgvmJFWHECYKygyZmXOkWY5bJOwQbNyML444VRZmfgcTk+bVt0fHzsy/tOW5muDz/8UOTHPzs7gyRJePToke/58TvTV/7oRz/ytbx5FVaQD7rcoOMEMF+xIqw44ZbtxVjhD9YpZqdOEWa5rFOwQTOyMP54QZf58uVLX9+fbr399tuhbvDl574D8y6sIB9kuYwTwQg7TgCMFX5inYImJexYEXac4BqaewjjjxdEmdlsFplMBl9++aU49urVK2QyGWxtbfle/jz44osv8PDhQ6yvr+OP//iPxfHvf//7gV2Dn73o9LWwgrzf5TJO+G8a4gTAWBEU1iloXNMQK8KOE2zQ3EMYf7wgysxms0gmk/jWt76F9fV1pFIprK2tYX19PbA5kZ9++imePXsGAOL/s1SmruvI5XI4Pj7Gn/7pn4pMM0H2ZGWz2cDKmmdhBXm/y52GOAHMdqyYhjgBMFYEhXUK/4QRJ4IsdxpiRehxwvfE0DPMtu2ZL9M0Tcc0zUDL3Nracmq1mthwK4iNt4IuM5lMtj0uFArO8+fPnVQq5Wu5juM4jx49chzndiO1N9980/noo498L3OehREngi43jDjhOLMfK8KME47DWBE01in8EUacCLpc1im4D83Ibm5uRMt3cXER5+fnWF9f971cN+2h3+mEO4UxJ9O27bYMIY1GY+bKTKVSbb01H374IX76058GkkLTXQT66NEjnJ2d4bPPPvO9zHkTVpwAwokVYc3dnvVYEWacABgrgsA6hf/CiBNBl8s6BZMCjGxxcRE//elP0Ww2Yds2qtVqIFktdnd3cX5+DlVVkclk8NZbb/leZljW1tawt7cHy7Kwt7cHWZZnrsyPP/4Yn376aduxH//4x0gkEr6WCwBLS0v45JNPkEwmsbi4GPq811kUVpwAGCtmqcww4wTAWBEE1in8F0acCLpc1imABYcRaizvv/8+6vU6Tk9PAy33+fPnKBaLKJfLM5055ODgAPV6HclkMrCdbsMoMwzuZ8gwDMTjcZRKpa6sNzQZYcUJgLFi1soMA2NFcFin8FdY9+w8xIppiRNs0Azp4cOHbdlAnL/K365pGiRJwtHRka/lf/755zg6OoJpmlhbW4Ou66FsCBmEjz76KPBc5mGU6fXkyRO8++67oZVPkxF2nAAYK2axTBfjxOwIO1YwTsxuucB8xgpOORtSuVwOtfx0Oo1kMgnDMPDbv/3boV6L3372s5/hv//3/45/8k/+yUyX6RXEPNf19XWcnp62fZE6joOFhQXfdnWfN2HHCYCxYhbLdAW1doaxwn9hxwrGidktF5jPOgUbNEPyLpx79uyZCAA3Nzeo1+u+B4Rf/OIXOD8/x9OnT1EsFgEgkN7esKiqir//9/8+VlZWACCQmyOMMoNUKpUAQHx+vvzyS7zxxhucFz9BYccJgLGCseL+GCv8F3asYJwI5p6d5VgxbXGCDZox5PN5EWwWFxdRLBZ9Dz5ffPEFnj59iv/23/4b/uIv/gKZTMbX8sL0b//tv8X/+l//C9fX1zg/P8fJyclMlvn48WMsLCyIqQbAbe/G0tKSL7n53cwytm1D13UAwPX1daA5+edJGHECYKyYtTKDjhMAY0XQWKfwVxhxIoxy575OEWSO6FmhadrAx374B//gHzi/8Ru/4XzrW99y4vG48/z5c9/LDEsqlXJubm4cx3GcRqMRSB71MMr0KhQKgZWVSqXE3gNh/KzzIow44TiMFbNYpivIOOE4jBVBYZ3CX2Hds/MSK6YlTnCEZgyyLOODDz6ApmmoVquBpAD8y7/8S/z85z/H4uIiLMtCJpMJJXNSEJaWlvDGG28AuP1dOwEMX4ZRppd3cajflpaWxHSHMH7WeRFGnAAYK2axTFeQcQJgrAgK6xT+CuuenZdYMS1xgg2aMXz88cc4ODjAZ599Flgqvng8PhUfmCDIsoy9vT1omoZyuYxUKjWTZXoF+fcM+2edF2HECYCxYhbLdAX9t2SsCAbrFP4K63M8L7FiWuIE0zZHxPvvv4+lpSXxgXEcBx9//HHYl+Wbg4MDVKtVfO973wss7WEYZbpubm4C3bHZ/VnX19e5r8SMYayYzTKB4OMEwFgxqxgnZrvcuaxTBDvDLdrceYGapjkPHz5s+y8IpVLJ0XU98HnUQXv69Kn4t23bbY9nqUxXpVIR5QX9szqO4xwcHPhe5jwJO044DmPFrJXpOMHHiV7lMFZMVtixgnFiNsud1zoFGzQjME3TcRzHsSyr6z+anM5gvrW1NZNluuXUajXn/fffdxzHEf/3u0wvLvSdLMaJ4MxLrAgjTrjlejFWTBZjRTDC+n6fl1gxLXGCa2hG4Kao+/TTTwMZOuzcSbjTLOUz93I6ZkFeX1/PZJnAbbrDjY0NVCoVAECj0fCtrPPzc+i6DsuysL6+DsdxcHNzA1VVfStzHgUdJwDGCtesxoog4wTAWBEU1imCEdb3+6zHimmLE2zQjCGo3V/dzYrmTRgZX8LKSLW2toa9vT1YloW9vT1fy11bW8PZ2RkePXqEH//4x76VQ7eC3CWasWK2Y0WQccItj7EiOKxT+Cus7/dZjxXTFieYFGAMuq7j008/RTKZxNLSEhYWFma2ZyMsBwcHqNfrgWaHCqNMb7mpVMq3jfK8Pvroo8CTHswjxolgzEusCDpOAIwVQWGs8F/Y3++zHCumJU6wQTOGm5ubrmN+ZZNYX1/H6elp21Cx4zgMeDPi2bNnYkfom5sb1Ot133eIzmQyeP/99wMZOZhnQcYJgLFiloURJwDGiqCwTkGTMs91Ck45G8PBwUFgrdFSqQRgPoaKewVal1+BNowyvfL5vAg2i4uLKBaLvgefVquFjY0N9gb6LMg4ATBWuGYxVoQRJwDGiqCwTuGPsO7ZeYsV0xIn2KAZQ5Bz491Fg2tra23Hnzx50nUs6sIItGEH9zAWDX7yySf45JNPfC9n3gUZJwDGilks0xXWombGimCwTuGPsO7ZeYsV0xInYmFfQBQ5joONjQ2sr6/j4cOH+P73vx/4Ncxi74obaM/Pz7G2tib+e/r06UyV6eUuGnzy5Ak++OCDQBYNLi4udv1HkzcNcQJgrIhyma4w4gTAWBGUaYgVjBPRLxeY7zoF19CMIei58b08fPgQn332WaBlBiWTyeDo6Eg8dodvZ61MV9CLBlOpVNsweDwe5zQSH0xDnAAYK2ahTCCcxcWMFcGYhljBODE75c5rnYJTzsawuLiIZ8+ewbZtcezdd98N9BoG5ZKPKjen+fn5eVtO842NjZkq08vNDhJk1hVvL1Gj0cDx8XFgZc+TaYgTAGNFlMt0hREnAMaKoExDrGCciH65wHzXKThCM4ZMJgPHcXB+fo6NjQ08f/7cl9ZoLBbrGWTcjCS//vWvJ17mNAgjp3lYedS3trbwwQcfhJod5IMPPsBPfvKT0MqfVUHFCYCxYtbLnIY4ATBW+IV1Cn+F9f0+r7EirDjBEZox2LaNP/3TP8XBwQG2trbw6NEjX8pptVq+vO+0CyPwhLUplDt3OsjsIN7MKy9fvpzJnrlpEFScABgrZr3MMOIEwFgRFNYp/BXW9/u8xIppiRMcoRnD1tYWjo+P8fz5cxSLRdTr9cDWW8yDJ0+e4OOPP8b19bXoOfL79xtGmcDXc6dt28bS0hIcx/F97vT5+Xnb41nLbDMtGCf8Ny+xIow4ATBWBIWxwl9hfb/PS6yYljjBBs0Ybm5uxAfk4OAAqVRKZLWg+/v2t7+Nzz77DPF4XBzz+4YMo0wA+Pzzz5FOp7GwsIBms4lnz57hu9/9ru/lAsCrV6/wxhtvBFLWPGKc8N+8xIow4wTAWOE3xgp/hfX9Pm+xIuw4wbTNQ/riiy/Ev70fyK2tLQ7DT9jGxgbW1tYCTQEYRpkAsLOzg3q9jp///Oc4PT3FD37wA9/LfP78Od58800oioLl5WX82Z/9me9lzgvGiWDNS6wII04AjBV+YqwITljf7/MSK6YlTrBBM6RcLodXr151HV9cXIRhGCFc0ew6OzvD8vJyoDn5wygTAJaWlkSQk2W5a1MsP2SzWdTrdfziF78ItHI0DxgngjUvsSKMOAEwVviJsSI4YX2/z0usmJY4waQAI+g3lMZZe5NVLpfnokzgNuDs7e1B0zQcHx8jlUr5XmZYlaN5wTgRnHmJFWHECYCxwm+MFcEI6/t9XmLFtMQJNmiG9PLly7AvYW6sra115eT3e5FZGGUCwMcff4yDgwN8/PHH+N73vocf/ehHvpfpDXjlcjmwytE8YJwI1rzEijDiBMBY4SfGiuCE9f0+L7FiWuIEkwIMqVQq4enTp/jkk0/wzW9+E8DtAqidnR1omsah+AkKcv+OMMt0vXr1CkdHR1hfX8dbb70VSJkHBweoVquBVo7mAeNEsOYpVoQRJwDGCr8wVgQnrHt2nmLFVMQJh4aWz+edpaUlJ5VKOclk0onH487jx4/DvqyZ8/DhQ8dxHKdUKjm2bTvvv//+TJX58OFD5/nz547jOI5lWU48Hnfef/99R9M056OPPvKlzPPz857HbdsW10KTwTgRnFmOFWHECcdhrAgSY0UwwogTQZbLOsUtNmjGYJqmY5pm2Jcxs3Rddxzn9vecy+WcVCo1U2V639swDCeXy4nHb775pi9lPnz40Lm5uen7HE0e44T/ZjlWhBEnHIexIgyMFf4KI04EWS7rFLe4hmYMzA/vr4ODAwC3v+dkMolcLjdTZV5fX4t/f/rppygWi+Kxn2kduQA1WIwT/pvlWBFWnAAYK4LGWOGvMOJEkOWyTnGLDRqaOt4bcGdnZ+bKzGazIn3j4uIifvu3fxvAbS537wZck8QFqDSLZjlWhBEnAMYKmj1hxIkgy2Wd4haTAtBU+OijjwY+78ciszDKdD1//hyWZeG9994Tx54+fYp4PO5Lbx0XoNKsmKdYEXScABgraDaEdc/OS6yYxjjBBg1NhVgshmQyiUwmA0mSuoYs/ejdCKPMMBUKBfz4xz9GIpEQmVf29vaYuYgihbHCf4wVFHVh3bPzFCumLU6wQUNTo1KpoFQqYWlpCZlMBu++++5Mlhm258+fA+C8bYouxopgMFZQlIV1z85brJiWOMEGDU2dm5sblEol1Go1yLKMXC7ney71MMokovthrCCiu4R1zzJWBIsNGppaT58+RT6fx9LSEo6Ojma2TCK6H8YKIrpLWPcsY0UwmOWMpsrz589RLBZRr9eRyWRQLpd9T1EaRplEdD+MFUR0l7DuWcaK4HGEhqbC3t4eTNOEoijIZDKBDMuGUSYR3Q9jBRHdJax7lrEiPGzQ0FSIx+Nt+dIXFhYA3G7QtLCwgJ///OczUSYR3Q9jBRHdJax7lrEiPGzQEBERERFRZMXCvgAiIiIiIqJxsUFDRERERESRxQYNERERERFFFhs0REREREQUWWzQEBERERFRZLFBQ0REREREkcUGDRERERERRRYbNEREREREFFls0BARERERUWSxQUNERERERJHFBg0REREREUUWGzRERERERBRZbNAQEREREVFksUFDRERERESRxQYNERERERFFFhs0REREREQUWWzQEBERERFRZP1/MEWTgXAZm2AAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "with plt.style.context('science'):\n", " plt.rcParams.update({'font.size': 9})\n", @@ -289,7 +391,8 @@ " for i in range(3):\n", " axs[1, i].set_xticks(\n", " np.arange(len(sims)),\n", - " [simname_to_pretty(sim) for sim in sims], rotation=35)\n", + " [simname_to_pretty(sim) for sim in sims], rotation=90,\n", + " fontsize=\"small\")\n", " axs[0, i].set_xticks([], [])\n", "\n", " for i in range(2):\n", @@ -722,13 +825,6 @@ " fig.show()" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, { "cell_type": "markdown", "metadata": {}, @@ -738,11 +834,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAD1CAYAAADUHqdoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0T0lEQVR4nO2dd3wcZ53/3zPbV221kixZsmxpJfcu2ymkY5mOgUS2Me1okeACgbsDK4YDjuPAkeF+cHBwSKGEcji2dEkwPVJCAgkJsS33Lq3c5CqtVmX7zszvj9GstVbvkj3v12us3dkpz4xnP/t9nuf7fB5BURQFHR0dHZ0hIU52AXR0dHSmE7po6ujo6AwDXTR1dHR0hoEumjo6OjrDwDjZBRgNixcvpqCgoM/PmpubycnJGfYx9f3Gbr/pUEZ9v+m930jPVV9fz4ULF4a9HwDKNOad73zniD4b6TH1/abuufT9bs39RnquzMzMEe2nKIpy01bPN2/efFPvN1Imspz6PZka+42U6XB9E31PAARFmb55muvXr2f37t2TXYxJ4Va+9v7Q70nf6PelN1lZWVy+fHlE+960kebNzmT8wk519HvSN/p96c1I2kE1dNGcpuhfhN7o96Rv9PvSG100dXR0dCYIXTR1dHR0hoEumjo6OjrDYFqLZnNzM+vXr2fHjh2TXRQdHZ1pwI4dO1i/fj3Nzc0jPsa0Fs2cnBx27949rRu6vV4v5eXlVFVVUVNTQ1VVFfX19VRVVY3J8QsKCnC73f1+XldXR0FBAfX19WNyPq/XG/e+vr4et9uN1+ulrq6u1+dTjaA/RGebb1yXoD80pLK43W7KysoQBKHP/8OamhoEQRizZ6Ungz0340FNTQ2rVq3q9Ty63W5WrVrFhg0bepVpuM/T5s2b2b1796g6gqb1MMrpjtvtZt26dezbtw+HwxFbv2HDBtasWTMm56iursblcvX7eXFxMUVFRaM+T01NDXv27KGuro59+/bF1m/bto2amhocDgdbt26luLh41OcaL4L+EC8/uw9/R3Bcz2NPtnL3u1dhtVsG3M7lclFWVobH46GyspKKioq4zz0eDw6Hg9LS0lGVp66uDqfTGfccDPbcjAclJSUAlJeXx5VFuw89r7O/520imNaR5nRnw4YNVFRUxAkm0OvL0Rd9/cL2tW4sBFGjrq6u34i0pKSEsrKyXuvXrVuHoii0tbWxZcuWMSvLeBAJRfF3BDGZDdiTrOOymMwG/B1BIqHokMtVVlZGTU1N3Lr6+npWr1495GMMFJH19byN9rkZaY2ipKQEj8czaM2nv+dtIpjWojnVq3oD4fV6qa+v7/PhdLlcsV/dqqoq6urqKC8vp66uDlDFKz8/n5qaGsrKymLb3Liuvr6egoKC2H7avtu3b6empobt27f3OndqaiobNmwY03vrdrvHrPo/EZgsJszW8VlMFtOwy+NyuXC5XHH/jx6Pp89IsL/nZdWqVbEmoJ4iXFdXh8fjYefOnbFq/nCeG61Jqef6vp7Fnttq2/dHaWkplZWVcefeuHHjsO/beDGtq+fDcSnxh6KcutQxjqW5zryZydgtA99arW2mvyqQtr6iooLGxkaKi4tJTU2lra2N4uJiiouL8Xg8VFRU4Ha7KSoq6nddz3OWl5fHqjOrVq2Kq/LU1dXxxBNPxAR7rPB6vbhcLjZs2MATTzzRK7LWGZyysjIqKysHbd4Y6Hnxer2UlpZSVFTEhg0bKCkpobi4GJfLxaZNm2I/4EN9btxuNzt37qS6uhogJtYVFRW9nsX6+noqKyupra0F1BpIcXFxn89CWVkZBQUFMeF0u91TqllnWoumz+fD5/ORkJAw6LanLnVwz5f/OAGlgr/++1tYkecccBtNFN1ud5/C6fV6cTgcNDY24vV6+2yUX716NQ6HIy5a7WudRk1NTdzD17MtqK6ujp07d8at83q9bNu2Lfa+L6EfrCmhpyivW7eO8vLyuChCZ2iUlJTw8MMPA+r/lSaCNzLQ86L9vzmdzmHVJPp7bnbu3BnX9r569epYk5P2XnsWy8vLcTgccdGr9sPeVzmLioqoqqpi48aNOJ0Df5cmmmktmoqicPjwYe64445Bt503M5m//vtbJqBU6rkGQ3uY6uvr+xTNuro6SkpKKC8vp6ysjKKionF9eFwuF6tXr6aqqiomdA6HI04U++owGIi6ujoqKipi0QWo1UqdkbFx40aqqqoG7KAZ6HkZLML3er39VvtHi9frZc2aNTHxHSxy1CJrYNQdXWPNtBZNs9nMwYMHhySadotx0Ohvoqmuru6zmqIJV1VVVaxqC8RSd/bu3Tui8xUXF8eiFe14GkVFRZSUlFBQUBCrso0WrddTY9++fWzatGnUx73V0KLCsrIyNmzYEKsO38hgz8tA0aX2Y9aXsPb33GzatCmuJrJ3795+2x7LysooLy+PdQYO1jxVWlpKWVnZlOy3mNYdQSkpKRw6dGiyizFiXC4X+/btY9u2bXF5mtqDp7UL1dXVUVNTw5YtW9i2bRtdXV2xNiLt4dN6tnuuq6+vZ+/evbEvWVFREVu3bo016Gs5lPX19bGIUmt7HG7HTV1dXezcVVVVcc0OWuN/QUHBmLeXjgeRUIRwcHyWSCgy5HLU19ezbds2tm3bFqvKaovb7Wbbtm14vd5YR0t/z4vRaIw1v2hNLh6PJ9YZs3XrVqqrq6mrq8Plcg3pudGq0Js2bYo9u7W1tVRWVvb5LBYVFVFWVtbrGANRWlra7/PS1/M2UUxrP80777yTlJQUdu3aRXLy4FVindEx3Or5dGOq5WnqjB+j8Rid1tVzTSgPHTrE3XffPcmlufmZSj2Y44HVbuHud68aVg7lSDBZjLpgTmOmtWi2tLSQkZFBdXW1Lpo6Y4LVbtEF7SZmx44d7Nix49Yee/6Wt7yFQCAw2UXR0dGZBozF2PNpLZoAy5cv59y5c7S1tU12UXR0dG4BprdoKjLLli0D4ODBg5NcGB0dnVuB6S2aUoi0tDRyc3N10dTR0ZkQpr1oglpF10VTR0dnIpjeoimHAVU0L168yLVr1ya5QDo6Ojc70zrlCEkVzaVLlwJw5MgRHnjggcks0ZDRRkzs2rUr5gpz4wiJgoICamtrY7ZgQ3W6GStqamrweDwxE4i+jDZ6jhHWGXu00TUulyv2fzEUv9WJor9npK6uDrfbjdPpxO129/JS1QxpNLTRaR6PB6/XO7W9V5VpzDvvWRB7/dGPflT5/ve/P4mlGT61tbWKy+Xq9/N9+/bFvS8tLVVqa2vHu1iKoihKY2OjUlFREXtfUlIS915RFKWiokIpLi6ekPLcirS1tSklJSVx6wZ6XiaagZ6Rns9FRUWFUl1drSiKolRXVytbtmxRioqK4o7V8zj79u2LbT9evPOd7xzxvtO7ei4FQZYAWLhwIceOHZvkAo0tEzFcsT83dq/Xy86dO2Pv16xZE+dW5Ha7dV/Mcaav8dRDcSsfa5OL4T4jdXV1cc9GUVFRbLv+HNe1cfHQ93VPJaa3aCoShFR3lkWLFuF2u2+aRPe+3LNv5EbX7JqaGlJTU2Pu3QUFBTEzB83LcqgUFRXFeWs2NjbGibjm6agzfhQVFVFXV8eGDRtiz0HPautwXNMHc14fzE29v/L19YzU19fH2dI5nc5BDWBWr15Nfn5+zHxjKhu7TO82TUWCQAvYMli0aBGyLHPq1CmWL1/ea1O/FOSE7/yEFGtBQi52g3VUx7jRPftG6uvre7lm19bWxtpGtQnTtHbSsrKyET+I2kyS2hdEm37gZvPGDAaDnD8/Mc9Ibm4uVuvgz8i+ffsoLy+PTUFSWVlJaWlpn////bmmD+S8PlQ39cHo+Yz0tIsbKhUVFXg8nlE9pxPFNBdNGbrOgXMhs2fPxm63c+zYsT5F84TvPKtee2RCirXvju9TlDx3xPv35+bek/5cs5944gmqq6tj4qm91hiJG/vDDz9MbW1t7MukNeLfbKJ5/vx5Pv3pT0/Iub73ve8xd+7gz4j2fwjqD+XatWspLi4elmv69u3b+3Ve78tN3eVyjeoZSUtLG1YTgTaNdXV1NW63mw0bNrB9+/Yp2xk0rUWzuTXE+g/9C5vLPGzevHnAds0FCbnsu+P7E1KuBQm5o9p/KF6D/aFNi7Bhwwa2bt3K2rVr4yZwG64b+/bt26moqMDlcuF2u2PT8dbU1MSimJqamikfHQyF3Nxcvve9703YuQZD87fUnoWioqJYlDkWDOSmPppnpKioiD179sQ+93g8A7bP79q1iw0bNgDXPWbXrVs3LqI5FoYd01o0c2Yks/ur98Eb3gWo7ZrPPPMMsiwjivHNtXaDdVTR30RRVVU1pKlZB3LNLi4ujkWYxcXFbNu2rV+374HQ5obRvrR1dXVxD3J9fT21tbU3hWACWK3WIUV/E8mN85273W62bt3aKxocyDW9P+f14bqp90Vfz0hpaWlc+3l9ff2Ajv2a2PZkvDpBN2/ezObNm1m/fv2IjzGtRRPRBIErcPlvkFvMwoUL+cUvfsGFCxeYPXv2ZJduQOrr62PVEa2xft++fVRVVdHY2Bhzzwb1oe/5fvXq1XGu2U6nkz179sRy5MrKyuKmJBioM6k/tGpST2788mrO2TdLpDkV0TrznE4nHo+HTZs2xardff3/93RNLy8vjzmsa87rmrhp/1+am/qN64fCQM9IRUVFn8fV2t615177Ydc6o0CNgKdSLuqNTGvn9vX3zmf3V+4C0Qir/xWfmMaGDRt49NFHectbJmYStenOze7GrjN6bsZn5JZ1bgcgMRc8h+HM70hInUde3hyOHTumi+YQ0dOGdAZDf0bimf6iKQhgz4YLz8GVv7Eov4CDx49Pdql0dHRuUqZ3cruGLR1S5kM0wMKMAOfPn6ezs3OyS6Wjo3MTcnOIJoDRCvZsFiVdBOC4Hm3q6OiMAxNSPdd6xbTGZM21py8XlMHcUQbEkkqW7TypSVaOHTvGbbfdNh6Xo6Ojcwsz7pFmWVkZq1evprS0lD179sTSXyoqKuImg9fGvfa3fkgIAoLVycJsI0ePHh3bC9HR0dFhnEVTG4+qpSpogtifC8pA7ihDxmhjaY6RkydPEA6Hx+AqdHR0dK4zrtVzbRiYFi1q1e3+XFBG4o7SC9HCklkGwuEIp06dYsmSJdc/C7VD1D+qaxoUox0sKeN7Dp1xIxwOE41Gx/UcRqMRs9k85O01v4CCgoJYkvvq1avZu3cvpaWl41jS/ulpkD1aBjIy1q63urp6yphdj6toejwe9u7dG3NR2b59O9u3b6e1tbXP7ftb3x/NVztY/7kdsfeb1y1hc/E8XE4Ju83CkSNHrotmqB32fg2CLSO7mKFiTYfVXxqycA70hVi9enWcu7u2veberUXlQ3FY1xmccDjMgQMHCAaD43oeq9XKihUrhiScbrebdevWsW/fvrha2IYNG+IMO8aC4SSxV1dXj4lgav4FWt9FT7OOdevWoY29qaiooKqqasQ/EtqYc40pO/bc6XTGjaPWxstu2rSpTxeU4bqj5MxIZve3NsevVBQMRFk8dw6HDx/mve99r7o+6lcF02BTo8HxQDtH1D8k0RzsC1FUVERZWVmviGL79u2sXbuWffv2DfjQ6QyPaDRKMBjEaDRiNI7PV0M7RzQaHZJoas5FN9q1VVRUDNv/sic3TjehHXOowxeHOzqoP0HWjIy151UzMt6yZQuNjY2x7bRgYqRoY841RjP2fFzbNPv7JSoqKoqzFdNcUPpbPywEAYAlc2dy7NgxJEmK/9xoB3PS+CzDFOOBvhAD4XK5Yj8ugzms6wwfo9GIyWQal2U4Yuz1euMcqnricrniOktvNCO+0cS6rKyMgoICoG+j4rq6OjweDzt37oxluwz12No4cs3oeDhNagOZXWv6oZmCTJVhnOMqmjeK4J49e9i0aRPFxcVxriaaC0p/6/sjPMCw+aUF6QQCARoaGkZ5FePDUL8QN+J2u+PMZwdzWNeZvgzmOuRyuWJmxFrGSVFREeXl5b1MrHv+EGsmGZpR8erVq2NORZs2bYrVaoZ67MrKytgxRmPconUcb926Nbauvr4ej8dDWlraiI871ox7nuYTTzxBeXl5rP1FC8P7c0Hpb31fvBLt4JwcYrZo6fXZ3CwrFovarjl//vwxvaaxYChfCA3tYdL2c7lcfbpr3+iwrjO90Z6B/kyptVrGQGbEA9HTqLgvhnrssrIyVq1aFfu8tLR0TMyugVgNdMOGDVNm5tNxF03toiFeALVfuxvpb31fSCjc2XmY/7XP5X5TjzZE0YQp2sbChQs5fPgwDz300OguYhwYyhdCe3gcDkfcPamvr2fVqlU0NTXFPWB9PXQ60xdN0PozpR6J5d9Q8Hq9w3LldzqdNDU1UVdXFxO10tLSURkZa/Z2WlPTmjVrhp9+OE5M62GU8wN2CkUrD/iO8hl/E36lu/1SNEOwlSVLlnDkyBFkWZ7cgvZBzy9EXwz0hSgqKsLr9cb8NaH3Q6dzc1BdXU15eXmvDtKqqipKSkrYtGlTnEt6TzNih8MR22+oAuvxeOK8WPs7dk+2bduGw+GgpKSE6urquA6codCXkbHD4Yjz6tSa9qYC09rl6LjHQX14Oc8lnGdr8Bx/jrbzWtJS7EYbBK6xdNE7+OUvuzh79iz5M7onsRrPPM1hHru6urrPyawGS63QHmqtN7E/92ydkTGeeZrDPbY2/cONaWmaeA1mRl1ZWYnD4cDhcMRNV3KjUTHA1q1bqayspKCgINaM1texbzTITktLo6qqKvaD3bNNcjD6MzIuLi6mpqYmlk7ncrmmTEbItDYhdriW43rbp/jr5os0il3c3nmYMnMm3zFlQOAqoRVf4qEP/xOlpaWsf/N90yZPc+PGjXEPtpanqY3Hr62tpaKigqKiItxud6xXVKOiomLKPGDTiamYpzmW9JVmNFFMNSPj0ZgQT2vRfOO9t3PQ9Sm2FF2i/A4/3w5e5J+DZ3gxYRH3dV6AlZ/nnyt2kpaWxhe/+EV9RJDOoEzFEUE6Y88t69ye6JhB6coQ3zmQyUeXneNR20yejrTy0UAjhzCR4LvE0qVLee6551AUBcGSoguazoCYzWZd0HQGZFp3BIHA5x68DYMAj79iwIDCT+2FXJIjfFaMoHQ0sXTpUtra2jh//vxkF1ZHR+cmYFqLZnNzMx955F9ZO7Odn5yaTcOldgoNNrabCvkRIR5vq2fJ4kVYLBZeffXVyS6ujo7OJLNjxw7Wr18/qrHn01o0c3Jy2L17Nz/88j+SlSTy1T3ZoMi8KZxHaTiPL0Su8NSVP3Hbbbfx8ssvT3ZxdXR0JpnNmzeze/ducnJyRnyMaS2aGjazka3vXsSz53M5eTmAK6WLh6VC3hPO4OOnKkldMYvTp09z+fLlyS6qjo7ONOemEE2ATfctZYZd5n8OOTCKkJ8SYYt/DsvFLP7H/hfMZrMebero6Iyam0Y0LSYDpXel8it3Lq1+mTR7BJso8WmKOBI5T9aSOex98Xm4dml8l872IZXX7XZTVlaGIAh9juCpqalBEIQ4x5mxoqCgYMJHDdXU1LBq1SoKCgriRkG53W5WrVrFhg0bYmWqq6vr5ayjozNVmNYpRzfy0bffwTdf+AM/2S/y+btgacIZJGke95gXcnzGOT7zwgGCj7dhtVrHrxCp6VD2JUgaOLXJ5XJRVlaGx+OhsrKyl3mBx+PB4XCMemRPX0nFY2UgOxw03wHNJUdDuw89r7OioiLOuLqmpmZU7jk6OmPJTRNpAmQ409i80k7ViVzCEQnRbKfA0sQnbGtpyvCQQoRLbV5wpI3PYrVBWwsEh55AX1ZW1stMtr6+fliGqwMZN/fldjPaURnDMYruSUlJCR6PZ0C/xTGZJ0pHZxy5qUQT4JF3rOZywMbTJwxgTiFFvshKWxZ3JS4gnAhnrlyDhKTxWazDd4R3uVyxKY01tLG2N6IZvZaXl8e2r6urY9WqVbHqbE8R7stY9kYDWW07LaLrWR3uy4C2LwPbnttq2/dHaWlpnL1XXV1dnAnEmMwTpaMzjtxU1XOABa7ZFOeF+Z/DGbx3sQS+ZhZkwoeu3Icv6TcEL7Xi9/ux28dpyosRoBkrDGaJV1FRQWNjI8XFxaSmptLW1haz0vN6vZSWlsa8B0tKSuKMZbXo8kYDWbfbTXl5ecyDc9WqVZSWlvYyO9bEWjNT0Axs+7Lx6suEpOe1FhQUxITT7XbHlWe480Tp6Ew0vURz69atcS7J2tB0QRDoOUxde79r1644+6ipwMfvzea9P29h/5UWVlpkBM8hHshbzqspTtoud9F07gyLFyya7GLGKCkp4eGHHwZUcdJE8Ea0idP66sTRIlOn0zms6rPmkKShiedgBrQ9DWzLy8txOBxx0avb7e7Xlb6oqIiqqio2btwYF1XC8OeJ0tGZaHqJpsvlin2Bh8JUNLx98733kPP0Ln5y0MT37k2Hs7/Dbm5gSeIcnks+yqETR1g0fyFC93xCU4GNGzfG7LX6o7y8nLKyMoqKinqJzWD/D5qx7Hh0AHm9XtasWRMT38Ei5p4O3Dd2dBUVFcX9CI9onigdnXGkV5vmjYJ54MABDhw4AMCPfvQjtm7dSkdHR7/bTwWMZhsfvnMG1Wdm0SFmQnIhtJ1ghlnAkZOO1BXm0PEjk11M4HqnSllZWcz+rS+qqqrwer1xju89p8EYKDrraSx7I8XFxb0iRLfbPWQDWq3sPSdz047RH6WlpdTX1/dZ5uHOE6WjM9EM2hFUW1tLSkoKTzzxBNXV1Tz22GPs2rVrIso2Kj5UvIKgJLLrmAlMCSBHEUItrJ13O8EUOHLi2KS2n9XX17Nt2za2bdsWq8pqi9vtZtu2bXi93lhHi9aOWFdXR01NDVu2bGHbtm0YjUbq6urYuXNnzJvT4/HEOmO2bt1KdXU1dXV1sYm49u7dGzcx29atW2MdQdrUCj3NbWtqaqitraWyspK6urpYG6YmbtpUwzceYyC0ybr6QpsnSrsGPd1IZyoxqJ/mCy+8wBvf+Ebe/OY3U1ZWxoMPPsjzzz/P2rVrJ6qM/bJq1SpycnJ6zWmssfmrP6apNcKr77+McLkRdr8OuffwUksjLXubyJ9fwJzs2aT1E90Nm6AfggH47DbImDk2x9TR0RkzduzYwY4dO2hubh7xBISD9p673W5SU1PZs2dPLDppamoa0cnGGs2woz8+dk8W73myndcvtXJ7UgbYgK6r3JmczXP2U3iaTpNrtdLWoSaSC4xBG2dq+ohSj3R0dMYfLcBav379iI8xqGiuXbuWyspK9u3bh6IoPPbYY6Snp4/4hBPJG99wD/lP/x8/OmDh9rclw9oFkHY75gUf5sUX/4uLT7zO25a8kSULFnHJbGbhwjHoHLLaBx0NpKOjM33pJZrt7e088cQTlJSUkJeXR35+Po8//njs856vpzqiNZmPvsHJ154z8o0uLxnp+RA8BgYPpW//NJtf+BjR+n0sfMP9+GSZ/Rcus3LlyinVq66jozO16NURlJKSEps/e+vWrTz99NOTUa4x44Nvux9BEPj5AcDiACkIF+qYZ89hxtsW0nHJw7lz57BarciyzOHDhye5xDo6OlOZPnvPNeHctm0b+fn5PPbYY2zdujWWejSdSHOk8NDKVH5yOg/J3wIJs+Dafmhv4LNv+ADnliv85MmfkJGRgcFgIBKJEAqFJrvYOjo6U5RBU45WrlzJ448/zrZt22hsbOSxxx7jW9/6Vlyu5qQRHpq4Pfy22znns/OnBllNP4p0wZndrEtdjvKmmUStAt/97neZN28eACdOnBjPUuvo6ExjhjX2/KGHHuKhhx6ivb2dyspKPB4Pa9as4cEHHxyv8g3MEN2EVhemU5SXwo+alvG23OchKQ+uvI6Y/hc+N38Tjz3wbZRn63n11VdxOp1Eo1Gi0ShG4003NF9HR2eUjMjlqK2tjc9//vNs27aNlStXjnWZhk7AN+RNP168gNrzyTSGZ0OwBYwJcPopPmyfQ/6K+QSWJFBZWUlGRgYAp06dGq9S6+joTGOGJJpPP/00P/rRj2JLeXl57LP8/PxxK9ygdHWAJA1p05I75pBiN/GLK3er1fOEbAi3YzjxM36Q+yCv3d1FUI7w5JNPAhAIBJBleRwLr6OjMx0ZVDQ3btzIc889R0NDQ2xpbGyciLINTiQE1y4OaVOb2UjJHXP41SEFKWEOtJ+GhFzoaGT1qV/xsexlHF4b5ZVXXuHs2bMAU+c6dXR0pgyDNtpt2rSJhx56KG7d/v37x61Aw0KKwnk3ZOUOafMP3lvAj19o4HlxI29KrAbfJRBEkCJ83dNEzVwF2+oZ/PznP+eRRx6J2d/peZs6Ojoag0aaqampvXrKp8owyuauIOs/+Sg7duwY0vZF+U4WzUrhlweicNtXoegxSHZBtJPUpDz+VRL4/R2XSUpJ4NlnnyUajXLu3LlxvgodHZ2JYseOHaxfv57m5uYRH2NQ0WxrayMvL4+5c+cyd+5cCgsLp4wdXI4jmd0l97O55KHBN0Y1Tv7APS5+V3+B1qAJ0pfBkkcgKR86z/Bw0gJSTQLCm6M0NTXx4osv0trayiCeJjo6OtOEzZs3s3v3bnJyckZ8jEFFc8+ePezbt4+9e/fGlqkimh7JxulLMpGzQ49833tXPrKiUP3qGXVFUi4s+UdIdmHraOTzBgc70tp4x32Z/PnPf+bMmTNcunRpfC5AR0dn2jGoaK5bt478/HxSUlJISUnB4XDwhS98YSLKNihh2cDBznRO/PXIkKPBjGQrb12Rwy/+0sMkN8UFq74ArhI+IRtwCCKHl15mYW4S1bt2cubMmfG5AB0dnWnHoKLZ1NTEt771LZ5++mleeOEFnn766SkTaZoMkGiM0HT8Mm1X2oe834fuK+DQuTZeb2i5vtLigHnvI2He+/gXrPzEKPHBNwYJ+LvY/ezTXLt2bewvQEdHZ9oxoGi2t7fz+OOP09DQwHPPPceuXbt47rnnRmzeOR7YLEbCHZ20nr065H3etCybgsxEvv/HG4ZLCgLMeTuPzHuYREQqHVE+da/MgUNH+O2v/2+MS66jozMdGTDlKCUlhcrKyl4u7VMm5QgQrFaM/iDnXj9Klz/KotsLsdjNA+4jigL/+OYFfP4X+zjX4mN2ekKPAwok5b2dL4SvUH7mKf5pvpFVbju/3v1b7r5zDXMXrxrnK9LR0ZnKDFo97zmNq8aUyls0GLGKETqaLtB44Bwn9rj522/2EwlFBtzt/fe4SLGb+GHtyT4/f6TgA2Rb0vhSUipb3hjEahL4r+98GzkSGI+r0NHRmSYMKpo3tl82NTWNuE2zrKws9rquri42adf27dsHXT8QtkQ75q5WhEAXl5qucflMCy0XvQPuk2Ax8uH7C/nZi410BnoLrNVg5quFH6Y6fI1TebfzgbsSaTjfwo4f/SfoKUg6OrcsQ5qNUjMifuKJJ1i3bh2pqanDPtH27dvjpmatqKiIm5FQm3mwv/UDIVgsJOEnJXiNDk8X0UiUK2dbBt2vbN08/OEov/hL38MlP5RdzOKEOTwme3n7HXNYMCuB3XWvc76+/3mJdHR0bm4GFU2Px0NKSgpr1qyhvr4+1ik0HNxuNw6HI/a+rq4u7n1RURE7d+7sd/2QsCdibm0my2kiIdnG5TMtRMPRAXfJcdp58LbZ/OBPJ4lKvc05DIKBb8z9KC94j1A3624+vz4LXyDCL6t/TaipVo84dXRuQfoUzQMHDsQtaWlpbNy4kXXr1nHmzBk++clPDuskdXV1FBcXx97X19fj7DFtrtPppL6+vt/1Q8JqU02JL53FarcQ9Ifo8AxuHffoWxdytsXHr/ec7/Pzd2bcwZ0pi/jCxT+S/YZP8ODtDl4+dIkXn/8TSsMukAcWZh0dnZuLPnvPS0pKWLVqVa+E8T179vDUU08Nq/e8rq6OjRs34vF4YutaW1v73La/9UNDAFsCXDyHcZaLaESi09OFM2vgmSGX5zl5YHEW3/n9MR68fXavTi5BEPjG3I/wwN7P84zUzvs/+A+8cPT77P7bWXJSYIkchrmbQdQNi3V0bgX6/Kb3lWbUk+eff37IJ/B6vTgcjjjRTEtLw+v19tq2v/X9ca2zna3/99PY+7ULllOclY1w5QKCKRPvtU7mDOE4n337Qt61/c/85fgV7luU1evz+53LeVPaKv614UnedecP+dQHGvm3H/yGQ2d9ZJpfICPcDvM/BObkIZddR0dn4tixY0ecsc9oDDt6ieaZM2cGFEwg7vMzZ86Ql5fX53bbt2/H4XBQU1OD2+3G7XZTU1NDUVERe/bsiW3n8XgoKirqd31/ZCSlsO2hj8Sv7GqH5jOYC2dw9YKHgC+ELcEy4PU8sDiLZbNT+c7vjvcpmgBfL/wIa/7+KX5x6QX+4R2fZOFzB3j+4DXmZc8m2f0clnAXLPmEOrJIR0dnSrF582Y2b94ce79+/foRH6tXm2ZlZeWwDjDQ9lu2bIn1hBcXF+NyuWKve/ak19fXs2nTpn7XD4uEJOjqIDHkwXu1g8N/PTnouHRBEPjM2xZSd/gSh8+19bnN6pR5PDTjbr7S+HOCSoQPP/wpmq91ceyanaOB+SjNL8KR/4GQd3jl1dHRmVb0ijSLi4t57LHHhpTAPlSTDLfbTWVlZSzSLCkpoaKigu3bt+NyuQBiKUb9rR8ygggmM4bmJlIWZ3Kp6RptV9pxZjkG3O3B22fz1ZqD/OdvjvLkI3f3uc035n6UJX8r5Ztnqvny8g+wcuVK/nzgAgseup2TgSALLr4EKLDoYbBlDK/cOjo60wJBmcZmkXfNXdy7eg4gy9B2FWV2IdeS8phblMfyexcMeryf/rmBR3/6On/7j7eydHbfuahbT/+Y75x9huN3/YjQuXY+85nPsKHkQVZkSyyN/gEzAUgvgkUfA8e80V6ijo7OOLB+/Xp27x5ZvvWIZqOc8ogiJDoQms9gD3m52HiVUCA86G4fuMdFQWYiX/u/Q/1u88X895FmSuafT1Yyf/587rjjDv768t+QUpdy1PJ2dZbLa3vg0HfBe3osr0pHR2cKcHOKJoDFCqIB28VT+K+10dLcd1tlT0xGkS8+uIw/7G/m76f7toJLNNr4z/mlPHP1Ff7UspcPfehDXLlyhf2HjiKnFRGYvRFsWeA5Bke+D57jY31lOjo6IyAajdLS0sKJEycG33gAbl7RBEhyYPB3wuVzXHZfHtIuD90+h8W5Dr5ac7DfNtuNmfdxf+pyPn3i+2TPmcU999zDyy+/TFSSOR5eAMs+DanzoeUQHPx/cOEFffSQjs4kEI1GaW1t5cSJE+zZs4fDhw9z8eLQZrDtj5tbNAUBUpxY2q/SVn9oSB1Xoijw5ZJl/PX4VV48eqWfwwr8YOGnORO4wjeadvC+972PlpYWDh5UhbYzYQks+ww4F0N7AxyrgsZqUPR51HV0xptIJEJLSwsnT55kz549HDp0iEuXLiHLMunp6SxZsmRUx7/5h7EYTZisJoJnzhA6dgTr4qWD7vLWFTmsLkjj688c4v7FmX1mEixMnM3W/Peyrekp3nvn/dx99928/PLLrFixgtOnT6v5pUsfgWM/gtaD0FgDUgTy14M5aTyuVGeUyLJMJBLB5/Ph8/no6uoiGAzGPlcUJW4ZDYIgxC09MZvNmEymXn9NJhMGgwGDwYAgCLG/tzqKohAKhfB6vVy5cgWfz0coFALAYrGQnp5OZmYmDocDs3lgr92hcPOLJmBKsOFv9eH7026s8+aDaeAbJwgCX3jPUh781ov8+ehl3rhkZp/bbc1/L09dfpHSY9/hZ5s/xSP/+AhHjhxh+fLltLW1kZq6AFZ+Do49AZdfg9O/gtZDMHcTpK9UI2GdCUdRFHw+Hy0tLbS1tSHLMrI8eC1AEzhRFEclVj1FV5blXgIcjQ7sZ9CXYGvlEkW18mi1WrHZbLFFE97Rln2qoP24dXV1cenSJUKhEJFIBEEQsFgsZGVlkZ6eTkpKypgIZU96iWZ7eztVVVVDztMUBIHPfe5zY1qoscYgCEi2JPzu46QdeBXW3DfoPsVLZ6rR5tOHeWBxVp/3w2ow88NFj/LGvVt4MfsUd999Ny+99BJLliyhqalJtdCzZ8Hyf4bk38CFOmjZD/5LMOdtkP9uMFrH4Yp1ehIOh2lpaeHatWtEo9FekaIoihiNRoxGIwkJCSQmJpKQkIDVasVgMEx4eWVZRpIkJEkiEokQiUQIh8NEIhGi0SjRaDT2WpIkQP0uapEyqKLS2dkZO+aN16stBoMBq9UaE1mLxRK7F0ajMSbCk4ksywQCAXw+Hx0dHbS1tcXuh6IosWvIycnB4XCQnJw8rv9v0zpPc17WLDKSUli7cAXFi1bG1kdlCEZFrEYZY/f/+VW/kaVCE/Nnm+GTXwF74qDHrz10kQe/9SLPfv4B1i7tO9oEeO+hr/Oq9zi1OV/h0Uc+zfvf/34WLVrErFmzyMzMvL5h53lVOC/9FYIecC6CRaWQnDfSW6DTB1ovqdaO1TOK1MTC6XSSnp6O3W6/KSKvniiKEhPXUCiE3+8nEAjg9/uJRqOxe9JXlHvjcW4UWC1ytVgsmM3mWHNBz2UoQqv9cGk/Djf+QAQCATo7O4lEIrHPtP9Hg8GA2WwmMTGR1NRUkpOTSUhIGNL/ozYGvbm5ecRznU1r0byjYDHbN8Qnt4clgbaQSKJJoSsikm6NYhChNWAgyx7kzsAeeNeH4YHBx54qisLaf38OQRCo+9K6fv9TjnQ2sfTVMn68+J9p/skeGhsbeeSRRzAYDKxcubL3fu0NcKQS2o5D4iyY+17IvlcdzaQzbBRFob29nXPnzsW+XNo9F0URu91OVlYWycnJUyJymipo0akWuYbDYYLBIMFgkEAgQDgcjkW9msD2JRc3RrFAXFttz9c924S1Y/Z3bK3JwWAwYDKZyMjIiNUCLBbLqH7sRpPcPqI2zQMHDgCq7dtkVs0lJf6mKQp4QyK5SRGWzQjwanMCHWEDqVYJm1HGE7IQSMjA9vIfYdkdkDZjwONrbZvv+daLPHfoIm9entPndkuS8nnPjLv4hvsp/vC+f+XRT32aU6dOsWDBApqbm5k1a1b8DimFsOZLcOp/4eJf1DHr3hNQsAmsw3fFv1Xp7OykqakpVk0D9f/MarWSnZ1NamqqLpIDoHUmGQwGLJaBTW20qFCLBqPRKOFwmFAoFPe3Z/R6ozhqzXnaYjQaY8JoMplISEjAZrPFOr20Zao0E2gMWzQ/8YlPUFBQABDnSDQZKEAwKmA2KASiAl0RA3ajzFxnCLtJweUIU3/ZRkQGq1GhNSDiSZxFzuV6eOWPsP5Dg55j7dKZvGF+Bl+tPsi6pdmIYt+/bv/qeh+rXnuE1+1nuffee/nDH/5AYWEhV69eJScnp/evojkZFn8CUhfBqV+C+9fgOQHzPwAZq/ROon7w+/2cOXOGQCAQ+3KKokh6ejo5OTkYjbdE3+aEo4mc0WjEZrMNe/8bRXM6M+wnbMOGDTFruI6OjjEv0HAwCAr+qEhbUMBqVFiYFiQnKUKqVW0cn50c5lKXkUtdJjLs6rprQRM5Gdmw50VYeTfkugY8hyAI/NuGFbzpP2p5+vWzlNyR1+d2RclzeVv6bXzdvYM/fuDLfKLsE5w6dYpFixbR2NhIYWFhXweHnPtUc49jT4D3FBz8Dsx6IxSU6P6cPQiFQhw9ejTWriWKIg6Hg9zc3DHvHdUZe24GsdQYdsyr2bVt3bqVDRs2jEeZhkyyWeb2bB+LM4LckeNjcXowJpgABhHmOkMYRLWt02pUuOIzEUlOh852eOm3Qxqpc+e8DN6yIpuv1RwiEu0/NeVLrvdz3HeOWuEoxcXF/O53vyMSidDe3o7PN8DUG85FsObfYN4HwGAB9zOw56vqUMxbHEVRaGpq4siRI8iyjNVqZcmSJRQVFVFQUKALps6EM6KGgscff5zS0lJKS0vHujzDwmaSmZkYZUlGkMyEaJ812nSbRJotSmdEwGaU8UVE2kImyJ4Nh/8OJw8O6VxfKVlO07Uuft7PzJUAdzgW8pHsN/PoiR9w+7vvp6uri4aGBgBOnhzE19OSAvM2w6ovQtad4LsI+78JJ34BofYhlfFmo7Ozk/3799Pa2oogCCxcuJAlS5ZgteppWjqTx7BFs6CggPz8fPLz81m1atV4lGlMEQTITY4gySIGAWRF7UknMQWkKDz/rDoh2yAsmZ3KxjvzePzZI/hD/Scff2/BPzLHlsknL1ey7i1v4tlnnyU5ORlZljl+fAjmHSkuWFkOyx5VRw65n4Y9X4FLL0M0MIwrn77IssyRI0diPzQZGRmsXLmShISEyS6ajs7wRfMb3/gGa9asic1OOR1It0WxGmWCUQGTCFf8JrVWPssFDYfV9s0h8MUHl9HaGeIHz53sd5sEo42dy77IaX8zR1YHiUQivPzyy4iiSCAQoK1tcLclRIMaba7+EuQWQ6hNtZp77Qtw7jnwX74pDUAURaG5uZn9+/cTDAYRRZGlS5cyZ86cm6Y9TGf6M+yOoCeeeIKVK9VE8qampjEv0LCIRoa0WYJJxmmNctlnIsEk0xkS8UVEEs0WNeJ8cTcsWDloClL+jEQ+9sZCvv3bY3zk/kLSkvpO01ialM9/LfgkZcf+i8+tu4tnnnmGN7/5zVy6dImmpqahj1iwZcCSfwTvWrWd03sKjv4PWJxq2tKM1arRcWLutM7xlCQJt9tNR0dHrId1zpw5ZGTo7vc6U49hf9M0wQQm/9c/GBhSxCUIMCclgiiAJAuEZYH2ULdoZc6Cqxfhpd8M6Vhb3rUEWVH4z98eHXC7h3PexgdmrqXStQeLzcovfvEL5syZgyzLHDrUv8lxn4VPXQBFj8GaL0PWXYACV/6uRp9//zL87fNw4udqFHrpZbi6T7WlazmoRqVTdG72UCjEoUOH2L9/P+3t7QiCQE5ODkVFRbpg6kxZhhxpfutb3+Jzn/scb3rTm0hNTY31ak5qrqbRBJ5rg0aIANmJEQpTQ5xotSAIcMVnJCcporq8Z+XCnpfUhPfCxQMeJyPZyqNvXci3fnOUT6ybz+z0vtvZBEGgctFnONjpxn23l87f/Zl3v/vdmM1mwuEwly9fJiur75kv+zkgJOfDin+BwBXVFb7tGLQcAF8ztDdejzYFAQSD+iNgSgSjXV3sGZCUD0lz1JFI9qxJiVD9fj8nTpxAkqRYcrPL5cLhcEx4WXR0hsuwh1E2NTWRn5/f6/VksKognxwlyOZ772TzbcsG3b7Fb+Cv5xMxigqiAG+c04nN1H35jcdg7lL4WPmgLkidgQjLP/8b1i2bSWXpnQNue9rXzOpX/5G7fmlggTOPb37zmzHfzaVLl44+ZSbSBZ3nIOSBcIfaWRTuhHAbBFog6gNZAikA0RAoURCNIBjBNgNmv1lNebJljntCfTQa5ejRozFTCYPBwPz587Hb7eN6Xh0djbEYez7sNs39+/eTlpbGrl27WL169YhOOlbkzJvP7rsKIBwcfGPAYZVIMMuEogKBqMgVv4m8lO65g2blq+lH+/4Kdww873uSzcSW9Ysp/996/unti1iQk9LvtnMTcvjxkn/hE3f9B/L/HePVV19l+fLlnDp1iiNHjvQ9Nn04mBJV0RsIWVIFNdyudir5mqH5RdVt6cgPwJyiRp+ziiFjJRiHP+JjILRaicfjibnSLF68eNChezo6Y402//mYzns+GIqisHfvXhRFiZujfFIwmeHut0BH25A6hYwiZNojhCQRo6Bwrt2ErMXZFpvqfPTnX6tV/kH4yAOF5DhtfP3pwdsnS7Lu5d41d9FeYKTqiSrMZnMsDWlCOtNEgzqmPTlPFcW8d8CdFWobaf671YT6loPq1Byvll/vdJKG1tE2ENeuXaO+vp7W1lZEUWTevHkUFRXpgqkzbRm2aDocDiorK9m0adN4lGf43PYAZOZCS99TU9zIrOQIJlHBZFBoDRg55bEQkrojvaxcuHwefr8DBjGCtZgMPPbupTy75zwHz3gGPe9/L/gU7vtFrrZc46mnnqKwsBBRFPF4PLS0tAyp7GOKaFB74Od/AO58HFb8M6SvUKv3p36ljkj62+egYSd0XRj24YPBIPX19Zw5cwaAnJwcVq5cSUpK/1G5js50YNii6XQ62bhxI62trTQ29j86ZsKwJ8KKO6FjcOECcFolMhMihCQBgwCHr1k52dod9RgMau7m/pfh9RcGPdb77s6nMCuJfx9gyl+NbGsa/35HKQ1rZHbu2sW5c+dYtmwZoihy9uxZ/H7/kMo/LpiTYeZdsPpf4favw6KHISFHrc6f/F/Y8+/gfhYiQytjS0tLbJx4SkoKK1euJDs7e/KzLXR0xoABRVOLEgBeeOEFPvnJT9LU1MRDDz1Efn4+n//858e7fENjYRFY7eDvGnRTQYCsxCiyAikWiUSTzNl28/WqekKSujxXA+cH/lEwGkS++OBSnjt4kVdPDV6lf3jW28h50yICKQrf/M5/IooiCxYsAODEiRODTnMwIdhnwOw3qdX3274KrnernUcnfqaOTGo72W9qlqIoNDQ0xJ6bBQsWMG/evCll66WjM1r6fJr379+PwWCgoKCAL3zhC7zwwgsUFxdTWVnJhg0b+PGPfzzR5RyY3AJ1aRnaNL2p1ihmg2riYTcqyArUX7Fzpr27Jztzlhq5/vZ/ITBwdPXgbXNYOtvBl3buH3SyLVEQ2bXqy1x7awoNx0/xzB92Y7fb4/I3p4wntCConUMLP6qOh09fDp1noX4bnOw9Hl6SJPbv309bWxuiKLJ8+XISEwd3x9fRmW70KZrl5eWsXLmStWvX8sMf/pCKigoqKiqorq7m4x//OI8//vhEl3NgDAZYeRcE/TCECbKSzTKJZglfRFRzx60yBkHheIuVEy0WIrIAs+fCyQPw/DMDJr2LosDX31vE30+38Os95wc9d6Yllacf/CYtS4xU/aiKS61XSE9PJy0tDVmWOXhwaAYiE0pKgVp1X/QwGKzqePi/f1Gtsgc9+P1+Dhw4gCzLJCYmsnLlSt3XUuempd960969e3nuuefYtWtXrCr+0EMPUVlZOam5mf0ybzmkOKF98LZNQYDC1BCKAIGo2s6WbJaRFTjSYuVYixXFZIGMbNWs+NjA+VwPLMniTcuz+dLO/YQi0oDbAsxLmEXFo18hisxHvvEZglKI/Px8jEYj0Wg0rllkyiAaIXedKp45D0DEByd/ztWXHuf4356BaIDc3FwWLlyot13q3NT0KZou13Vj3uLi4phTe1+fTxnSM2H+CvBcHdLmuUkRChwhOsMGwpKAIECKRSbJLOP2WrjkM0Jquhq5/u5/B636f/29Kznf6qey7tSQzr929hre8/BmOOrlwR9/hpAcZvny5YiiSGtr6+T0qA+FxFmw7FG8C7ZQb3wXZ7tSoO04C7tqyOx6WU2219G5ielTNG+MFAZ7P2VYepv6NxIedFNBgMXpQWYnh/CGDLEauM2ogAAnW60Eo4Lam37xLPz6ZxDqP4l+QU4KH7m/kO2/PkJL59CS7T/1tn9gyX2riPzmDCXP/ythJcKyZcsQBIGzZ8/GTcE6FVAUhcuXL7Ovvp7TF7uQUxdjyS9mxeo7sZuBpmfhtS/Cmd/q4qlz09LnMMrVq1fH2b7V1dVRXFwc936y5weCPmaUC/jhO1sh4IOsWf3v2IPOsMhfzycgywLJFrU9VJKhJWAkNzlMUVYAc9QPZ0/D2gfhHe/vd7jhtY4gKz7/Gza9IY//9w9rhnT+QCDARz7xMc4JHhyPFPF/RV9BDkY5ceIEiqJgsVhYtGjRpMy/rSHLMufOnaO1tTU202NycjIul+t622XEBxeeV6coDraCNU21t8t7pzpqSUdnCjGa2Sj7FM3BUkQEQYhNUj+Z9Hnhf3gK/rQL5g0+Fl3jVKuFQ9espFklDN2XHpHAEzKQmxRhdkqYmZErCO2tsOmTsOqefo/1vT8c51+fOsDLX3sLS2cPbWbJkydP8tl/+ifO3AbZ71jC7hVfxSQbYlM8aO4/wzL4GAOi0SgNDQ10dXXF5sCeMWMGOTk5/T8jIS9cfAku/lUdD29KhOz71CGaFj2xXWdqMOaiuXHjRioqKvrcQVEUHnvsMXbt2jWiE44lq1atIicnJzaeFIAzp+CH/w7pWWri+xAISwJ/u5BAS8BIoknCZlQQBFU4vSEDogh5yWEKQg0kJRgRPvL5fidkC0cl7vziH5iRYuX3W9cOuSnjqaee4smf/Yzj7zaQu3IuNcu/RJbFSWtrK2fOnEGWZQwGAwsXLhzRbIDDIRwOc/LkSUKhUGyseG5uLunp6UNvmpEi0PwCnPsDBD1gSoIZq2DOO9RcUB2dSWAsDDv6FM3nn38+NuNkXwz2+UTR56+FJKmieb4RZvcxA2Q/dIZFDl214QkYiCoCTosUq4UHowLtYRGrQSHBd4X0rBSy3vd+MubP7lNE6g5d5D3fepGff+pu3nPb7CGdX5ZlvvGNb/D63j0cfp8Z/wyBXy4tZ13aKmRZpqGhITb7Z0JCAnPnzh3ztJ5AIBCzbNPEsqCgYHRDH6MBNfI89ycItnSL523geo8eeepMGmMeaU4X+r3wvz0Hu34IBYvVHM5hcMVnZM8lOyjE2jhBTdUMSQKhiEKoowtzegZz17+RuXfMw2jqfY6N336JI+fa2PP4O0iwDE3cgsEg//Iv/4K3o52Wj2dRGzpMef5G/r3gHzCJRoLBIMeOHYs1jYiiSGJiIjk5OSQkJAy7g06WZTo7O/F4PLS1tcWOazAYmDdv3tgmp0sh1Ry56VnVXcniVI1Dsu8Dgz6jpM7EoovmjbS3wQ/+DbraISdv2Md1e80cuGLDbpSxm/q4PZEwPq8PnyMb15vuYsW6Zb3a+NxXOrn9C7/n42vnsu19RUM+97Vr13j00UeZOXMmGZ9YzZfP/oKipEJ2LPsCLvtMANra2jh//jyRSCQ2D7ggCBgMBrKysrBYLCiKEhtdpP2VJIn29nb8fn8smtTQ9l+4cOH4zvYoheF8LZz7ozq23TZD7SyasQZMuq+mzsQwpqLZ3t5OVVXVkKIWbT6Xz33ucyM6+WgZ8ML/VgvVlTBnLpiHZ0MmK3C8xcpJjwWbQSHB3Mcoo0iYkLeDTudsljxUzNzb5/a6Z9/+3TH+bddBXvjKm1jlShvy+Y8fP86WLVtYvXo1xY88yAePf5OWcAf/b34ZH8l5EwbhemQbDAa5ePEiXq93yJ1zmlu6yWTC6XTidDqxWq0Tm0oWaofGarjyOkT9YHXCzLth1jrVxk5HZxzRI82+CIfgiW9A00lwLRj2sRUFTnksHGuxYjYoJJrk3plGkRBdXj/htFkse++byV+ZHyc8UUnmga/+iYik8JevvhmzcehNBa+//jpf+9rXKCoq4tHH/ol/aXyCn12spSipkP9a8I/cnbqk1z6yLNPe3k40Go2VQxCE2GtRFElISMBkMg37fowb/stw4QW48pra825KhLSlapunfWKzBXRuHUYjmjev/YzZAg+8S50DqKtj2LsLAsxzhliaEUBWoCVguG5YrGGykOiwY2xt5vCu57h4sjnuY6NB5Psfu4MTze18+7fHhnX+2267ja985SvU19fznW/8J1XzPsMrt30bURC5Z88/8476L/GaN34edVEUSU1NJSMjg/T09NiYdi2adDgcU0swQRXGee+DOx6HhR9RRfPqHnXCuL9/CdrdN+V0xTrTl5tXNAEWrFAnS7t4dkRfPEGAQmeYu2b5cNokWgLGPoUzOcWGcu0SJ3b9gVBXIO7jZXNS+ezbF1Lx66McPjeEOc97sHr1ar761a9y8OBBvvKVr7DMOIe/3/5d/nfpY7gDl7jz9c+wdu8Wnr7yMpEpOuPkkDHZ1THtd2yDpZ8GW7raYVS/TRXPlkPqtB06OpPMuFfPa2pq8Hg8NDY24vV6qaysBNRRRW63G6fTidvtZsuWLQOu74shhdgXmqDqP9SpMdIyR3wdnWGRvzcn0BYUsRkVLEYFi+H6rZODAVo6FRa/YR4LP7opbtRQKCJx/7/9CVlReOnf3oLVPLwe/YMHD/K1r32NlJQUvvzlL6tWcorM01de5ptnqnm94yQzzA7+IXsdH8l+MwsTh5bmNKVRFOhww4knIXBVnYbYmqYmyWfeBhbHZJdQZxozZavnbrcbt9tNaWkpFRUVeDwetm/fDkBFRQWlpaWUlJQAqrgOtH7EzMqHO4pVw41RjGJKMsusmunHlRrGalToCou0BFSzDwDRasNuEXH//QTtL8a7vltMBn70iTfQcLmTr9YM3/pt+fLlfPe738VkMvHZz36WV155BVEQKcm6l7/f8T0O3vlD3pt1Pz9u/iOL/vZx7vj7o3z37DM0+JsHP/hURRBUS7rbvwZrvgIzVqtDNU/9L7z6GNRvB89RPfrUmXDGVTS9Xi87d+6MvV+zZg21tbXU1dXFzXFdVFTEzp07+10/au56M2TNVuf/GQWpVolVWQHW5nWyJttPZkKUtpABqTvgTEiyElBMHHvmz0jH48Vxca6Dr2xYzn//8QQvHRuaWXJPsrOz+fa3v82qVav42te+xn//93/j8/kAWJbk4r8W/CMX79tBzfIvkW5K4fOnfsTclz/CvJc/wmdO/IA/tuwhIIVGdf2TRkI2LPkk3PF1tf3TlASeI6pwvrZVnVkz4pvsUurcIoyraBYVFcUNVWpsbKSoqIj6+nqcTmdsvdPppL6+vt/1oybFCQ+sV02KhzAlxmAIAsxKirAqy0+6LUprQHVJEgRIdZi56Ddz7pdPQfOZuP0+9eYF3LNwBh/9n7/R7Bn+nEA2m40vfvGLfPKTn+T555/n4Ycf5uWXX47lW1pEMw9l3sNvi75G6wM1/HrFV3nAuZxnrr7CW+u/iPPPD/HWfV/gv84+wwnfuanjEj9ULA51nvY7H4eVWyB1vjol8bEn1Ojz+E+6O44GN6LW0RkpE5Zy5PV6WbVqFfv27WPbtm1x7Zv19fVs2LCBkpKSPtf3N4GbNvZcI24M+o1Eo/DUD2DPn6FwybBHCvV7XUEDey7Z6QiLpFklRAG8QRFLVwv3LBCwf/Sz4MyIbX+1PcC9X/kTWQ4rf/zCumG3b2pcu3aNH/zgB7z66qusXr2aD37wg8yfP7/PbRVF4bjvHH9o2cOfWvfykucwYSVCjiWdNzpXsNa5grVpK5llzehz/ymN75IaaV75u5osLwhgSYPsuyHrDWCbhtekM+ZoY841xnzs+XiwYcMGKioqcLlcbN++ncbGxl7iWFZW1uf6/kRz2I25rVfhR9ug7dqwxqUPhjdoYO9lO+0hkXSrhAK0+EQKo26Wr8hC+MCjkJgc235/k4c3/Uct77ktl8rSO0eVVP7KK6/w5JNPcv78eYqKinjf+97H4sWLBzymLxrgL22Hed5zgOc9+znY6UZBYZ59Fg84l3Nf6jLuS11GtnXoCfmTjhSGthNw9vfQdV718zRYwDEPCjeo8x3p6HQzmo6gCZnIZfv27THBdLvdFBUVxflxejweioqK+l0/ZqTNgLdsgl99V50WI8U5+D5DwGGVWDHDz2sXE+iKiCSZZZIsCmfFPGYePkjm0z+GjZ8Aq+pOtDLfyX9/7DY+/sNXWTjLwT+9fdGIz33XXXdxxx138Morr/CrX/2Kz33uc8yZM4cHHniA+++/v087uQSjjbdm3MZbM1TT5pZwOy96DvK85wAvth2k8sLvACiwZXN36mLucSzhrtTFzLPPQhSmaJaawQzpy9Ql0qWOcz/7W2g5qIpp6nyY9361fVRHZxSMu2jW1NRQXFwcmyKjrq6O0tJSysvLY9vU19ezadMmiouL+1w/piy7Hdzr4MXdYE+CMUr2TrdLzE0NcfiaOmbdZlLwRQ0cty7Cue8VTDY7vPujsfNtekM+Jy928OWdB8hJtbPxDXkjPrfBYODee+/l7rvvZt++fbzwwgvs2LGDJ598koKCAhYtWsT8+fOZP38+mZmZmM3xBhnp5hRKsu6lJOteAK6E2vhL22Fe9h7hr21H+MXF55GRSTUmcXvKfFYkF7A80cWyJBeF9mzM4hRLmDclQs596lzuLQegoVrN82xvUMe4579Lr7brjJhxrZ673e5e8wtVVFSwZcsW6urqqK+vj0WfPfM0+1rfFyMOsbs64McVcK5BHWI5RmOuIxL8rTmR1oCBdJuEJENr0MiSxFYWBI5D8YPw1veqo5RQ2xrLql6j5rWzPPP5+7lv0dgNGwwGg7z66qvs27ePEydOcOHChdhnSUlJOJ1OkpKSMJvNmEwmTCYTBoMBo9GIwWDAZDJhsViwWCwIZpFrRh/nBA+nucJJ8TJnLG3IJjAKBgrt2cy3z2KuPYe5CTnkWjLItqYx0+wkzZwcN1Z+UpAluPq6Ootm4JpabU+cDXM3QXLBmP3/60wf9LHnI6HxGPz0m2rSe8bMMSvTNb+RV5vtmARIMMt0hUUkBe5IukiG/zy87X1w/ztjwhmJymz4fy+yp7GV3z62lpX5Y9NkcCOdnZ00NDTQ0tKCx+OhtbWVrq4uIpEIkUiEcDiMLMtIkoQkSYTDYUKhEKFQiEAggM/nizkqadiSErA47UhOM52pEpeT/DQlddDplFG66zACAk5TEmmmZFJNiTiMiaQYE0g22kky2kgxJsTWJRqs2A1WErqXRKOVRIMtthjFUYqvFIFre6Hp1+qYd0FUh3EWblLHu0/VpgedMUcXzZHywq/hN79Q7eOG6PI+GIoCJ1otHG2xkWqNYhSgNWggySxzp62JxEALPPBuWPdQrKreGYiwvuIFGq90srv8jazIGx/hHA2yLOPz+Whvb8fj8XD16lWuXr3KlStXuHjxIhcuXKCtTR0majAYSM+ZQUpuOvY5aRjykgimG2hX/LRFumiP+uiM+umQ/LRHfbRHfPjlwXNILaKph4haSTLaSTRYSTEmkGJSxddhTCDVlBT7m2pKJM2UTIY5BYcxUW2TVWRoO64myvsugiJBwiy1zdO5SBfPWwBdNEdKJKKaFb/+gmpYPEZO6FEZXm1O4JrfSLpN6jb8MJKdGOY223mMbZfhzmJ4xwfBpnpIen1h3v3NF3Bf6eI35W9k+RQUzsHw+XycPXsWt9tNY2MjDQ0NuN1uJEnCarWyYMECli1bxrJly5g3b15c22pEjuKTgvilED4piE8K0iUF1CUaoFMKxNZ1RtXXnZKfzmiA9qiPjqgfb7SLtoi6hJVIr/IZBJE0UzLpphTSzclkm9PIQWSm9zSZ4XYyFZHMpDnMyH836Rm3YzRMSD+pziSgi+ZoaG+Dn25X2zcLFo1Z+9alLiOvNSeQYJKxGhUiMrQFjSxIC7LYehXh0hlY8QZ4z8cgRfWP9PrCvGv7CzRd7WLXP9/HHXOnf2dFKBTi9OnTHD9+nCNHjnDkyBF8Ph9ms5mFCxeydOlSli5dyoIFC7BYhud72h+KohCQQ3gjPjyRDlojnbRE2rkWbqc10kFLuJ2rYS8XQx6aQy1cCnnokuKNVgQgw5hIji2THEs6WRYnmWYHWRYnMy1Osi1pZFvSmGlxTr2OMJ1B0UVztJw5BU9+Ux2bPnNszC4UBQ5csdHgtZBmjWIUwR8RCMsid+Z0kSF2wdlTMG85lDwMmWqSvtcXZtN3XqLe7eGn/3gX71g1tKmIpwuSJOF2uzl06BCHDx/myJEjdHV1YTQamTt3LosXL2bx4sUsWLCA1NSJMyP2S0GuBK5w9eJfuNr8AleCV7loMNFsTqLZlsYVJcqVcBtXQt5eUWymOZVZ1nRmWdLJtc5gtjWDObZM8myZuGwzSTMlT6zBs86g6KI5Fvz9BaipUp2QksfmyxqSBF6/aOeKz0S6LYoowLWAgRn2KG/I8WGQwqpJcmaO2qu+/E4QBIJhiYcr/8buvRf45gdX8fDa3q7wNwuyLHPmzBmOHDnC0aNHOXr0KC0tLQBkZWXFUqXmzZtHYWHh+E7FoRHugMuvwtk/gO8CGO1qcvy896OkzMMT7eJSqJXmUCvNoRYuBFu4ELzGhVAL54PXOBu4Sqd0fZhsstGOyzaTfFsW+bYsXLaZFNhnUmDLJs+WiUnUmwEmGl00xwJZhl8/CS/9FvLmD3uKjP7oDIu81qwmvadZJSKyOi3wqqwAeSlh9bzNTWqUe+c6ePMGsCciyTJf+NV+fvDcSTbflc93PrwG+xAnaJvOKIrC1atXOXnyJCdOnODEiRM0NDQQDocRRZHZs2fHBHTevHnk5+ePWbW+F9GAOh3H2d9197YLkJAD8z4Iqf2nqimKgjfaxZnAFdyBS7j9l2gMXKIpcJlG/yXOBq8QVbonsRNE5lgzYyLq6hbWAns2BfaZJBsTxufabnFuWdHsc97z0eD3wS+/A8fqoXBxLC1otFz2GXn9YgJGQSHRLNMWNGA3ydyd24XN2H372z1w5YI6Lv6dH1TnNgKeeqWJR3/6Oq7MJH72yF3Mz771pr2NRqOcPXuWU6dOcerUKU6fPs2ZM2eIRqOIokheXl5MROfOnUt+fn6vBP5RociqDd2pX0HXORCMqm3d/A9Bct7wr0eWuBC6RqP/Eo2BizT41UUT1Z5RaoYphQK7KqYu+0zyrJnMsWVSYJtJrnXG6NOwbjHGbd7z6cKYRpoaV5rhJ9uhvQVmzx2zw55stXDkmpUUi4xRVLgWMDLfGWRpRvB6wBIJw9nTkJQC975djTztiRy74OX93/0r51t9lL9rCZ9528JhzTd0MxIOh2lqaqKhoYHTp0/HhFSSJIxGI3l5ecydO5d58+Yxf/585syZg2G0Ji2KDNf2Q8NONVXJYFWHZ87dDIlj0/asKAqtkQ4a/RdpCFyk0X8Jd+ASjf5LnAlcpjnUioL6lTUJRvK6BVSLTgvt2cy15+CyzcSqT43cL7dspDkuoglwdC/873fBaof0sRmlI8lQf9nGmQ4LGbYooaiAPyqyZqafWck9OhYUBa5dAm+3aD+wHpbejl+Cx589wnf/cJz52clUvH8V9y/WJx7riSakPSPSc+fOIcsyFouFgoKCWPvo/PnzmTlz5sjaiuXo9eGZgSvqCKPUhbDgI2AZ35pAWI5wLnhVjVK7hbUpcBm3/zKNgYv4pCCgDiqYY53BvIRZzLPPYn7CLObZc1iQMJtca8ZN20Y+VHTRHGsUBV54Fn73K8ieM2aJ74GowKvNCbR1D7P0hgxYjTJvmOUj6cZpgqMR1Y9TisLCInjjuyFvHofOefnsk6+zp7GVexdm8qWSZTdFatJ4EQgEaGhoiAnpyZMnuXxZNYFOTExk7ty5saWwsJCsrKyhC4osgecwNOyCYCuIZkjMVYdnJs6e8OGZiqJwOezhtK+Z0/5mTvmbOek7z+nu6r/W659gsDLfPouFibNZlDCHhQm5LE7Mw2WbectU93XRHA8iEdj5P6r/ZsEiMI5NLp4nYOD1S3YCERGnVeJad9L7mpl+TH09r36f2lFks8OcebD8TpS5S/n92TD/8X+HOHLey8o8Jx9+oJANd8whyabnDA5Ge3t7rEqvRaRaj31iYiIFBQUUFhbG/ubk5AxctVcUNfJsrIagR11ndULeesgoUh2YJhlJkTgbuMoJ33mO+87FlmNd5/BGVWNus2BifsIsFnWL6ZLEPJYk5lFgnzn5/gFjjC6a40V7m5q/efb0mCa+X/UZef2SHUUBu0nGEzSQlxJm2Yxg3GRtMRRFNRlpvQKRECQ7Yd4y5HnL+OM1Mz891MFzx1sxGw3cvziTt62cRfHSmcxKs9/y1bCh0tbWFmsfbWxspLGxMRaRWiwWXC4XhYWFsXbS3Nzc3kKqKNB5Vo08fRdUj0+jTa26z36r2u45xf4/tOj0WNc5jnad4bjvPMd8ZznadZbWiDr1tUU0sSAhl8UJc1jcLaRLEvPIs2VOXavAQdBFczw5cwp+UqE6vWeMnRdjk9fMgas2bAYZowhtIQPZiRFWzAiQcGNVvSeyDO2t0HpNHTNtMILFxgWDg/+T8vh9RwqvtZuREci0KKxKgyWpIvOcJualmclzWklNsKidGpKkVjFlWf3Ci2L3YlD/Ct3vDQb1tcGgtvPaEq4vVtuUE4KxorOzE7fbTUNDQ2y5cOECiqJgs9liPfZaO2lmZub1H6lIF5z7E1x+DSIdgAjmJHDMhVnrIGn2lB7jrigKV8NejnSd4UjXGY75znKs6xxHus7EIlO7aGFxYh5Lk1QRXZwwhyWJ+cy0OKf8j7UumuPNC7+G3T9X04AsY5NcrShwtMXKiVYrKRYJg6Coxh4mmUJniDnJ4b6r6zcSiUAooC7BAEQjtEhG/h5ysC+iLsejyVxSbLFdUoQws0U/s0U/cwx+5hh85Il+8g1d5Bn82ISeoq2A9oQoijo+32hS3aFMZvV+pM6AjCxwpEOSAxxp6g9MsuOmE1Sfz9erjfTq1asAJCcnx6LRwsJCtY00cwZCRyOc+yN0nlPFFEWdHM4xD7LvU9OXpkmCu6IoXAy1cririaNdZznSdYbDXU0c6zpHoNt0Jc2UzLLEfJYm5bM0MV8V1MQ5JBntk1z66+iiOd6EQ/Dkt+DEftXYY4yEQJLhwFUbbq8Zp1XCKKjJ8EFJxGmLkm6LkmaTmGGPDE1AB6AjKnI6aOJsyMSZ7r/nQybOhoycC5kIKtejnlxzhEJbmPm2MAttYRbaQyy2hUk2SGp0Gg2rYh0Nq2lSoSCEgyArgAIGEyQkQWoG5M2DrFzInKUu9psvWdvr9XLy5Mm49KfW1lYAEhISKCgoUBeXC1emhdnKYYxdTRD2qpG+JUWd1yjnfjUStWVOux8bSZFoClzmcGcTh7uF9FBnEw3+i8ioP8IFtmyWJ7liy7LEfPJsw+h4G0N00ZwIzrvV+YUgNk58LIhIsPeSnQtdJpLNMjajgiRDR9hARFaNI1KsMnOSw1iMMum2KHbT2P6XyQpcjRhwB024g2YagiZOB82c8JtpDJqRUB/qfEuYpQkhViYEWZkQYmViEKexj6aESBh8neDvVDuyFBlMFnWepKxcdcRV5iw1Ik1MUSNS0+R3lowlHo8n1jaquT1dvHgRAJPJxJw5c3DNSqcwNUBBSgcuRxCbSQHBAOYUSHFB5p2Q7ALrxI3BH2sCUogTvvMc6nRzsMvNgY5GDna58UQ6AUgy2FmWlB8T0hVJBSxJzMNuGN/hsrpoThQv/Rae/SnkFsbm+xkLIhIca7HS6LVgFBVSzHIs0JBk6AwbCMvqCrtRZk5KmAVpQQwT0CQWkgVOBswc8Zs55LNyyG/hoM9Ch6SGvvmWMEWJQYoSQqxODLAiIYT9xs4sRVGjdV+H2qEV6nYUMlvAbFWXpBTVDDprtjqXU1omOGeo6V7TLOrqD5/PR1NTU6ydtLGxkbNnzxKNRhEEgewZDgrToWCGSGGGQkGGiZREC1hS1V74tKWqiBonYPz9OKJV8Q92ujnU5eZgp7qc9F1ARkZEZEFCLiuTC1iZVMjSxDyWJuWTZR67tlJdNCeKSBh+9v/g6B51uOMYfpkVBc52mDnRasEXNqAAJlHBapBjTYpWo4I/IuKPiuQ7Qsx3hkgcqNNonJAVaAya2O+zUt9lZZ/PykGfhYAsYkBhkT3Eqm4hXZUYZKEthOlGgdeENBxSq/ahoCqm0e5Ef7NVFcwUp+o8lZENqelqu2listohZbWP2RxPk0UkEuH8+fMxEdWWQED9YclINlEww0ThDIHCTCOFOQ7S0tIQcu5XRyMl5Eyb9tDBCEghjnad5UBnI/s7G9jf0cjBzsaYQXW6KSUmpFpkOs8+a0SGJ7poTiQXz8IT31C/3GNkI9cTf0SgJWDEHxG54jPRGRYRUYUqKAmkWSWisoA3LJJkklmQrnYaTXYwFpHheMDC3i4r+7qF9Lhf7cW3CDJL7SGWJIRY1v13iT1EUl/pVaAKaigIgS4I+CHoV5P8FcBsVqv6RpO6WG1q+2liiiqmCcnqOqtdfW9P6v48We3tn+wbNQRkWebSpUsxIW04dZKGhtN0dKlj0h12URXQTDNzc5KZ65pFxqK3IKTOB4tzWlzjUJEVWW0r7WriQEcj+7sF9XzwGqDmli5JnBM32d+ypHycpuQBj3vLiuaYG3YMldeeh2d+rH4JZ4xd++aNKN1CaRDAHxE5dNXKFb+JJJOEzajQHhaJygJ5KWFmJUVItUZH3WE0lvgkgUN+C/u7rBzwWTjit3A8YCGqqF9qlyXMYnuIRfYwC20hFtrDFFrDmPtrdlAUNdqPhNUfrV5LVO1YQVCFQ1HUSNRkUdtMrXa1HdWZcT1itSepEa22TNE0KkVRuHbtGg2nT9Nw4hCnTx6joeksbZ1hAFLsIvOyzBRmJzJvTgZzi4pJz1+jzro5xa5lLGiLdHK4s4mD3dX7A52NHOk6Q0hWayqzLOmsSCqIE9MC+0x2PbVLN+yY8EgT1C/jK3+C3/5STbnJnBij4EBE4ESrlfOdZsISOCwykqK2eYqCgt0kk2mPkmKVsBpljIIanEkKGAUwiup/dUgSiEgCCgKioGAzysiKQEQWkBW188kgKhgEMBsUTKKCKCpYDAriKL9/WhvpYb+Fwz4LxwJmjvstXI6oVSwDCgVWtee+5zLXFiaxv8i0P3qKbDikDgzQmgRk6XoqlSheT6EymdW21oRktZ01OVWNYrW8VE1Ye76exCaC1mtXOX3kdU4f38/pE0c51dyJ16fazjkTDczLtjEvN5V5K+9nflExSc7MSSvreBOVJU75L8TaSg90NHKg083lsDpKyyZaWJqYR1HyXJofq7s1I81JE01Qv5Cv1sFvfq5+0bJyJ+y0nqCB0x4LF7tMGAVItqjzEAWiIsGoAAgoKGj6pgCiQOy9HPc/LiAISmy90L0OFARB3c8gqK/NBoV0W5QUi0SyRSLFImMerpD1Q2tE5ETAwomAmRMBMye7l4vh64KUY44w1xpmni1CgTVMgVVNjZpjiWAcjZhL0nVxjUWxWlpVRO39pztyBVUkjZrImlTxTE69vsQS/+3qj6rJogqxxaoORpC7U7ekqPpXFK83NxhN6jEN3fmww3RmUmSZa80NnD78GqcO/51TZ69w+qKfrqDa9j3TaWFBYT4Llq5iwdLVuAoKMU3zduHBuBpq43DXGQ52qiK6v7OB/Mcv66I5KSiK6vi++2fqQz4ObZwDnfpsh5lTHgvtIQN2o6xGhN0RpSp7KtprTSwNwvUam6J0rxfU9T2Pr0WpsiwgAxFJICQLKIqAQVCwGBXSbFGcVolEs0SqVcJqHNvHqVMSONUtoKcDZk4HzZwKmHEHTYS6c0uNgkK+RRXQQmuYudYIc7tfZ5qksa2dKooqpDeKrPZakoj7tRIEVfgMRnURxesjsBT5uhBro7Bu/Gs0Xm9eMFvAYgNLd9aBJsiaeBu7t+m5mCwoso+LDS9z6sQ+Tlzs4OTFII1Xw0QkMBkF5s7OZOGiRSxctoaFi5eRlpY2hjdsanLLtmlOumiC+tDveRGefVKNDLImdk6fQFSgsc3ChU4TUUlAUgSisiqT3VoI3dVwQQARoDuyFLkeScYEtls8RUEV4P6q41EZQpIa2UqKup3FoOCwRmMRqMWg/rUbZWwmtZo/VgImK9AcNtIQNNMQMNEQVEW1oTtxX+6+8kRRptAWZq5VreJrglpgHUF1f0QFla9HlFJUfS+K6o0WRPXGKVwfzjro3x6vJQl6/C/H2nFR1NeisYdgd/8VZIh2EY600BiMcMIvcDxg5ESngatB9ThZSWYW5uWycP58Fi1bTf7CxRgSksbMlHsqMBrtuDlyFSYTQYDbHlCjjad/DB1tYzbH0FCwGRWWZASZ7wwSkVXB9EdFAhFV0BQEJBmCkkhUIpbvCRCRRHUf5fo6RYFgd9umrFz/SqopUGoalFFUMBsUEkwyCd01O0mBYFTkms/Epa7uWwMoCJhEVTBNBoUEs0SiSSbRLJNklrsFVR52W6koQK4lSq4lygM3WFiGZIGmkCkmpqqgmnixw861yPVHfqYpGotOC7ur+gXWMHmWKBZxjARVFFXLuImuAWuiKmuC3aM5wJCMGTsLxSgLDX7eY+4Ai4eWoI/jYRPHwhaOH2/nr4cbiNb8HqsgM98ms9CZwMJZmSwoLCQlKxccMyElXe1cS3GO2RQxUx1dNMeKO4pV1/cXf6NWmyxjl/w+FEwGMBnUKCPZMvTcTUkmTjRlmZj4an+jsoAvIuAJGOiKGAhJAl0REVmJF1KLQcZujO+sVRTteBCKCvjCJi5qNVKBmJimWCRVRE0yiSb1r9WoiqnSHTgJDK0j2CIqLLCFWWALA764z7xRMSaiapRqZm+XlZ0tyfhlNZISUJhljpJvjZBvjZBnUYV0tiXCHEuEDJM06g6xcUczXxmGWqcD98gR7ol0QcBDqL2Z021+jnep0eifLsJTF87Ba+eYKYaZbwqxwBxinl3BlZqENT1TzSbJyIH0XMiYA2mz1UyFm6jdVBfNsUIU1RklW6/AkdfVMeqjnV5hAjCIao/19RVg44YoKxpVE8+FAHI0QiAi0BE24A0ZaQnY6JIthDDgU0Tkbucei0HGYlQwi2A2iJhFA5iNcc4+UrdAR2SBK10mmm8QU6PYs0tLiQmn2aBg7naHMneLrtbLb+xetP2NooLVKMfacR1GmdWJQVYnBuMuUVHgUsSIO2iisXs4aVPQxEGfhV+3JuKVrv9fWgSZWZYoueYIud1immuJMtscYZYlSrZ5DCPViUY0qSOQLKlYHAUsmQNLZAmkAEq4iyudHRz3tHOyXeCkz8YrXYlEOgXEKwqzjVcoNJyh0BCgwBDGZZFIsBnV3NrkBHCmQPosVVAdGd0RarraLisa1YhcNHW/7v6rtf2iqK5eclT9C+qQU83nU5G6t9PWi+r+QrdbF6LqiSB2L6NAF82xxGqD93xUnSTtzClw9T9j4ZREUdRE8q4OdchjONzdNibGeoJFm52EVDsJtgRm2hPBYkWKSARDEoFgFF9Aoa1T4mq7TCCk0BFVj6vIMkbZj1WIYBIkjMgYjAYMBiNWzTnJqD7ksqJ2OkUVATkWBQuxzqmwJCArhlgnVnfftrpN99Zaj79BUDCIqogndEevVqPazmoxKnFtrzNNquDdnRzodWu8UZHzIWPM6ORc2MiFkIljfgt/9CbEVfsBZpii5JjV46l/I2SZJWaaomSZo8wwRXEah98sMSmIBhATEUyJZCVkkZUFD3R/FImEOONto8Hbxen2AKc7QrzkixJRgE7IMkXJM0bIuxIkX7jAHOE02YYQZqOgVo+MBkiygMMOCVYwGdX1ZpO6mAxgMYLFoP7CKzKg1aSE6z/CitJjfY82Y7V+cv29JqijQBfNsSZtBrznI+pwy0vn1OkypipSVDXW6GoHX7dlmcWm5icuXg05+WpbVUp3m1WKU216uOGHwAAkdC/pwBxAikr4vD58rR0EvZ0E2trxNHtob+kkGAghhbrdkUJBCEYwyRHMih8zEQwoWESwGLqFVOt5Ngztge8ppnJ3m60/ItIZMiApXI9cUb+Hxu6sA6tRxmqQsWiRa3dHkdidbuUUFTLMIW63htSONYi170oCXAwbOR8y0Rw2cjFs5ELYyMWwiVc6bVwMJdEmxdc8jIJChlEiwxQlwySRbpJIN0qkmSTSjBLO7iXVKJFqlHEYJRLGsDNtLDCZLMzNyGJuBry1e11UljnfFaSx3UdTh5+mDh+1HX5aQ1FAvZ8zrQZmWyDHKDErFGVWS5hZxggpQhiBKEgRtdPKIHRXhwSwmMBuAbtVzSCwGLtF1qRGkQZR3UeJqmJsMam1G4tBTVSmu5NMGXrzVV/oojkeuBbC298H//cjuHoRZoydefGIURTVb9OnRZEhVYASk9XJ49YshOw8tfd/Rs6ofUMNRgPJ6ckkp8cPZwv5w/g6AgR9IUKBMOFAGL+nA+/FFvyeDtr9QaRgGCESUpsEwmEIRVFkCUFSR7/EeqZipslaBCHE9UxfzwBQE/WNRhkRYvmn2m2RFFAUAV9YpFNRhbVnLqt2Oi2aFboFUzuOKKhNAnajjN2gsMQQYbU1hClBbSIQURNgo4pAqyzSFhVplYxci4q0Ro1cixpoiRg4GzRRH7XiiYq0RQ2xDICemAQFh0EVUodRxmGQSTZIpBjk7uW6wKaZJGaY1YjWYZT6H2k1xhhFkfxkO/nJ8f6ZneEo5zoDnOvyc7YzwPnOAH/1BbnqB6VbihKMBrITreQkWMm2mplpNZBtFsgyglOJIkQj0BZVm4wEUP+JqkIp0B1xCiBEwaSoTbomESxm9VlPTlE9X/nNyK9vxHvqDMya+1Vh+v2v4IJbjdomMkTQ7Nl8HRDwqQ+TxaaOdJm/QrVny8pVl9T0CSubxW7GYu9tA6coCqFAGH9HgFAgQjgQJugPE/IFUQIBLIQxS0EMYT+K34fs64IuL0ZfO2KwCyUaJRJRQIoSiSqEJIEwRsIGGyGjnZBoRcZAFEHN2lHUqrxW5VebAZQ4mYol+gtKLOlVQUHqkYGgZhmISAqgqNsKXE/limUDdS+ayNoEmIPCHCJABMEIglFBtF4/rk8R6JRFfIpAlyziU7pfKyJ+RcQXEWgNGzgnm+hSBHzd24T7EFsAm6AKq8MgkWKUSTVIOLuj2hlmNdJ1mtToNi0W5Y5dE0KS2cjitCQWpyXFrQ9JMhd9QZq7gjT7AjT7glzyBTnc2kFr8PpMrRaDyEy7hawEKzPtiWTaLcy0W8m0W8iyW7D2nNY65vvavQQC0NkO5xqvm8KMkGmdpzlpY8+Hw4FXVTu5rg7VkHc8ct2iEfB3XfewlOXrRsBpMyBvgZp4P0ZR5JQjElHbYrXF74NOL1y5AGdPwbVLKF0dRCNRZAxIqGP2ZQQkgxnZlohkSSBsSiDcnTAfkVSlEMXu3Fa0kVUKJlEdZoqiRo9hSSAYFegKi3RFDES682Vl5XruOlxPo7wu1vGpXXL3PihqNKsJr6g1yfUYrtDzldgj4hUFkFAIKCKdiki7JNIpdy+SSFcP4dVEWHsfUHo/myIKKQaZtG4hTTdJzDBFyTRHmWFSxXaGKUqmSSLTFCVhjHNfg1GJy/4Ql3xBLnb/veIPccmv/o30qBKkmI1k2S1kaotN/Tuj+7XVaGDH3sPs2HuEZlOiPvZ8StNwFGqq4GqzKmDGEQb4MfcfX7cDkE/tgjYar9uozZmrVrNn5Ny0U04MG1+nmg7W1X597Hm0u0310jl1ts+uDvWHB0XtlLL1MPEwdvfi9kxU1zqvbri3iqLmwoYlNVVLuR6kxuI/tfqvCqTU3eaqCm33ehkk1M/CUnzqlxQT4h65tcSLsbZdrFlBAJOg+Qko3U2EapOF1KMMQQnaZQMdkkCHIuKTVeHtVES6FCH2ulMW6VBE2mWxVxOCXZCZYeoWVGO3mJolZpoizDRLzLJEmGWJkGwcfdusrCh4ghEu+4Nc9oe44g9xJRDiavfrq4EwkhIvqpqg/r2lS09un9IULoYP/TPs+iE0HlU7Wpwz1B7p/tB6srUqdjgECGoCsS0B0rMh16VWrzOyVQNfR5oukH2RkKRmMvRHJKKmirVcgmuXVPu/5iY1Wr12URVL6J5krrtDKhpVhRdBjdztiZCQhGC1YzHQ96yiY0BPcexpVK1FrZIsEJQEQlGBUPfgBX9EoCNkIBhV34ej6oAHLUgTUHN8bQaFZLM6jv/GXFs1EiYWQcuKeox2WcArGfDKIh7ZQLss0q6ItEdEGkIW9ski3j6i2ARBJs0gMcOoZhTMNEfJMUXJsUSZZY6QbwuTZpExD9DxJQoC6TYz6TYzS/oY+SkpCp5gmCv+EJf9Ia4GuoXVHxrV/4EumhNF9hz4h3+BvS/BodfUqmM4pDroWG09DHlD18cm27o9IeevgNmFaoeNs9vV/CZyNJ90TN3DX3sOgZVlaGtRxdTXERvHHRvTHfSrn127pLZZX2lWX4cCgHDd4zMhaUzzdbVc1Z4BnmiAHrPf0d8sTD2jYLVJQcQfEekIi3iDqqgGwyKSlsOldGcXdEenRlHBKCjdfW3q+ZKBXKRe54mLfgG/JOKRRFpkA62SgRbJQKss0ioZOBKx8FKXjQ4l/j4lCDJOUSLDIJFplMg2RZhpipBtDJFjCJEtBkgRQhjlKAZFQkTCoEgYkDCgoBgMJIkiyaKBuTYjotWCmARiBN51LD5PdzjoojmRpKbDuofg/vVw7jScPgyH/q5WsxOSISdd9XpMcarCOCNHjSBvtjbI6YAodk+7MaP/bVwLr78Oh1ThvHweLp6BxmPQdk1tklFktRNOM0qepOGGgqBGwNej4OtipyiqZaA/IsbcsnwRka6wiD9qICwJBKKqoEqKoLaqCqidX90+Bz0HJahCq3QPKlBIESRSDBHytWaOmGFJd8iryEQkmVbZwNWokauKiauyhRbFREvExLGwiZex0HmDZNmRSRUkUkUJpyiRapBxGhRSRQkHERxKGAdRzKhTyAhGE4LFArww4vuoi+ZkYDJBwSJ1Wfsetao3BY1vdYaB2aI2l+S6gPtUYfBcgyvn1XZT9/Hr7afRqBp9JiZf9+qc5P97QVCnU7EaJegjcozKEIwKhMIKwWCUSERCiEZQolGkiERUUuiKGumULQRkdXbTqBLf5qkIAqIgqJ1WYnf6pSiqqZhGE0arkSyzmSzN19TYI0/XqNrlhTDQKom0RAVaowItIZmWkExrIMrloMSxUJQ2v6RmM/TAZhBINokkCyJJYQMpo7hXumhONtoDonNzIYqQnqkui1er63ydaiR66Rycb4QzJ6H1qlrVB/WH054EiUmq9dtEC6k2b1MocH1a5nAIIRpV0x1RSDIYu59ZCyRa1OaHJEe3l6gDbAlEBRNB2dDdrtrtY6CIapNAWCEUhVBIJhiSicoKYUVAisrIktzdcSWg9U8LgoAoCgiigCiIiAaRZItAqkFkvkHEYBAxGA0IPfKiZEWhMyzRFpLwBKN4Q1Js6QhLtIelvq5+yOiiqaMzUSQkXa9hgNoBdbVZXa5cgLOn1b+XL6iCBdcnmNMc4g3dxsSDCarS3c0ud7sbRaM9/kauGyxHwtc7ukAdJ27pnmMpYyakzlCblZJS1CakxGT1OrS5l/rIBDECid3LgEWUFcKhCOGgukRCUSKhCJFwlGhYIhqJquvDUaIRiXAgQjgUIdr9Xg5GkaISUrRHMwPdQisIJAmQBOSZBASziJBiRDSIiKLAD08MUrgBmJKiWVdXh9vtxul04na72bJly2QXSUdn7DGZICdPXTT8PrXHvuWy2pt/wa2KaLtHFTmpW/y08dQ9/TRBfa90/xW6zYwNhlj1NmZqnDpDjRBTUq93WGnO8ylOdd1IU+OGiCAKWGxmLLbh1bSkqEQ0LBGJRIkEI4SCEUL+MOGgKqihYJhoWEKRFRQUFFlBispEQlHCoQhS5CaMNCsqKqitrQVg+/bt1NTUUFJSMsml0tGZAOwJaq7tnLnX10XC4G1Vh8GGAterz5FuQ5VY1ClcNxuOTZ9hvO7ybrGpf82WaW0obDAaMBgNWDAzksZJWZb5wdPfHPH5p5xo1tXV4XA4Yu+LioqorKzURVPn1sVkVqvKOmOCOMofjCn3c1NfX4/T6Yy9dzqd1NfXT2KJdHR0dK4z5USztbV1sougo6Oj0y9TrnqelpaG1+sd0rbNzc2sX78+9n5KG3fo6OhMGjt27GDHjh2x983NzSM+1pQTzaKiIvbs2RN77/F4KCoq6nPbnJyc6WHYoaOjM6ncGFD1DLaGy5SrnhcXF+N2u2Pv6+vr2bRp0ySWaGrS81dTR0W/J32j35fejCbSnHKiCWrKkZZqBOg9532gfxF6o9+TvtHvS29uOtEsLi5my5YtlJSUjDixfaQPynTZb6RMZDn1ezI19hsp0+H6JuMHYUqK5lgwHf7DR7PfSJkOD/TNfE8mY7+RMh2ubzJEc1o7ty9evJiCgoI+P2tubiYnJ2fYx9T3G7v9pkMZ9f2m934jPVd9fT0XLlwY9n4wzUVTR0dHZ6K5aavnOjo6OuOBLpo6Ojo6w0AXTR0dHZ1hMOVGBOnEU19fj8PhwOl0snfvXlavXo3D4bhlPUe9Xm+cC1Z/9+FWuj833hP9mYGamho8Hg+NjY14vV4qKyuBMXpeFJ0pTUlJiQIoDodDqaioiK0vLi6Ova6oqFCqq6sno3gTRnV1tbJlyxalqKgobn1/9+FWuD/93ZNb/ZlpbGyMu+6SkpLY+7F4XvTq+RRn3bp1KIpCW1tb3K/ijZ6jO3funKQSTgwlJSWUlZXFrevvPtwq96evewL6M+P1euOubc2aNdTW1o7Z86KL5jTA7XbHeYrqnqMq/d0H/f7c2s9MUVER+/bti71vbGykqKhozJ4XXTSnAV6vF5fLxYYNG/B6vbrnaDf93Qf9/ujPjIbX66Wuro6tW7eO2fOidwRNcUpLS2Ov161bR3l5OQUFBUP2HL2Z6c97dTierDcj+jNznYcffpja2locDseYPS96pDmFqaurY926dXHrNH9Rj8fTa92tRn/34Va+P/ozc53t27dTUVGBy+XC7XaP2fOii+YUxuVyxTX079u3j02bNumeo930dx9u5fujPzMqNTU1FBcX43K5APXHZKyeF33s+RRHyzcDtX2mZ29ofX197Ff0Zs65A/V6a2trqaqqoqKiIvaF6O8+3Ar3p797cqs/M263u5eRT0VFBVu2bBmT50UXTR0dHZ1hoFfPdXR0dIaBLpo6Ojo6w0AXTR0dHZ1hoIumjo6OzjDQRVNHR0dnGOiiqaOjozMMdNHUmdLU1dVRXl4el3w8kQw2vK6qqory8vKJKYzOlEAXTZ0Jw+v1UlZWRmpqKlVVVWzfvp3y8nK2b9/e7z61tbVxQ+HKysoQBKFPEa2pqUEQBKqqqsasvNqxvF5vbBx3T0pLS0lLSxuT8+lMD3TDDp0Jw+FwUFBQwMaNG+NMJQoKCigqKqK4uHjA/bUhgh6Ph8rKSioqKuI+93g8OByOuGOPhl27dlFSUhIru8vl6jWuW+fWQ480dSaU2tpaNmzYEHvv9Xpxu92xMcJDoaysjJqamrh19fX1rF69ute2o3H2aWxsjCuXNn5Z59ZGjzR1JpS6ujqqq6sBVegqKyuprq4elmi6XK7YuHNNxDweTy/RrKurY8OGDZSWlrJmzRr27NlDQUFBLBIdaLxxfX09a9asib3XxFcT6z179vSKdHVuEcZuZg4dnYHZt2+f4nK5lNraWqW6ulopLi5W9u3bN+A+W7Zs6XWMxsZGpbq6WikpKYmtr62tVdra2hSHwxG3fUlJSdx8L0VFRbFj9Jxbp6ioSGlra+v3vFp5NVwuV+x1z/lodG5+9EhTZ8LQIkMtOnQ4HJSXl1NbWzvsY5WUlPDwww/HHbe/qnjPKLa4uJjKykoKCgriqto9p0foi9ra2lgn0I2zP+rcWuhtmjoTRm1tbVxHitvtjjN/HS4bN24cdk95WlraoIJXU1PTy0+xrq4uVv3X2zZvbXTR1JkwbhSbnh0tw5nkS4soy8rKqKioiJsUqy96pift3LkzZjxbV1cXt4223Z49e+Kcu7XzaWKrHWOsUpt0phe6aOqMO263m+3bt1NcXMyuXbti6zdt2hTLhRxKdbe+vp5t27axbdu22PQF2uJ2u9m2bVtcbqXGnj17qKurY/v27WzdujW2z9atW9m+fTs1NTWxDiGv19sr79LtdsdSj0CdEnbv3r199tbr3AJMdqOqjs5A3NghM1xKSkoG7WzqSUVFRVyH0FD30bl10CNNnZue4bSbtra26p08OgOii6bOlGekY8+1PMzq6uohJbm73e643MyhUFVVxc6dO4ddNp3piz5HkI6Ojs4w0CNNHR0dnWGgi6aOjo7OMNBFU0dHR2cY6KKpo6OjMwz+P0w6HVUZNxyKAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\"]\n", + "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"manticore_2MPP_N128_DES_V1\", \"CLONES\", \"CF4\"]\n", "\n", "\n", "with plt.style.context('science'):\n", diff --git a/notebooks/flow/reconstruction_comparison.py b/notebooks/flow/reconstruction_comparison.py index b7a5948..bc57281 100644 --- a/notebooks/flow/reconstruction_comparison.py +++ b/notebooks/flow/reconstruction_comparison.py @@ -153,7 +153,8 @@ def simname_to_pretty(simname): "Lilow2024": "Lilow+24", "csiborg1": "CB1", "csiborg2_main": "CB2", - "csiborg2X": "Manticore", + "csiborg2X": "Manticore V0", + "manticore_2MPP_N128_DES_V1": "Manticore V1", "CF4": "Courtois+23", "CF4gp": "CF4group", "CLONES": "Sorce+2018", diff --git a/scripts/field_prop/field_bulk.py b/scripts/field_prop/field_bulk.py index 8657d56..6c9126f 100644 --- a/scripts/field_prop/field_bulk.py +++ b/scripts/field_prop/field_bulk.py @@ -54,6 +54,8 @@ def get_reader(simname, paths, nsim): kind = simname.split("_")[-1] reader = csiborgtools.read.CSiBORG2Snapshot(nsim, 99, kind, paths, flip_xz=True) + elif simname == "manticore_2MPP_N128_DES_V1": + reader = csiborgtools.read.CSiBORG2XSnapshot(nsim) else: raise ValueError(f"Unknown simname: `{simname}`.") @@ -72,7 +74,7 @@ def get_particles(reader, boxsize, get_velocity=True, verbose=True): if get_velocity: fprint("reading velocities.", verbose) vel = reader.velocities().astype(dtype) - vrad = np.sum(pos, vel, axis=1) / dist + vrad = np.sum(pos * vel, axis=1) / dist del pos collect() @@ -282,14 +284,14 @@ if __name__ == "__main__": parser.add_argument("--simname", type=str, help="Simulation name.", choices=["csiborg1", "csiborg2_main", "csiborg2_varysmall", "csiborg2_random", # noqa "borg1", "borg2", "borg2_all", "csiborg2X", "Carrick2015", # noqa - "Lilow2024", "CLONES", "CF4"]) # noqa + "Lilow2024", "CLONES", "CF4", "manticore_2MPP_N128_DES_V1"]) # noqa args = parser.parse_args() folder = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells" if args.simname in ["csiborg2X", "Carrick2015", "Lilow2024", "CLONES", "CF4"]: main_from_field(args, folder) - elif "csiborg" in args.simname: + elif "csiborg" in args.simname or args.simname == "manticore_2MPP_N128_DES_V1": # noqa main_csiborg(args, folder) elif "borg" in args.simname: main_borg(args, folder) diff --git a/scripts/field_prop/field_los.py b/scripts/field_prop/field_los.py index f164576..1ba948f 100644 --- a/scripts/field_prop/field_los.py +++ b/scripts/field_prop/field_los.py @@ -179,7 +179,9 @@ def get_field(simname, nsim, kind, MAS, grid): simkind = simname.split("_")[-1] field_reader = csiborgtools.read.CSiBORG2Field(nsim, simkind) elif simname == "csiborg2X": - field_reader = csiborgtools.read.CSiBORG2XField(nsim) + field_reader = csiborgtools.read.CSiBORG2XField(nsim, version=0) + elif simname == "manticore_2MPP_N128_DES_V1": + field_reader = csiborgtools.read.CSiBORG2XField(nsim, version=1) elif simname == "CLONES": field_reader = csiborgtools.read.CLONESField(nsim) elif simname == "Carrick2015": @@ -429,10 +431,10 @@ if __name__ == "__main__": Om0 = csiborgtools.simname2Omega_m(args.simname) # r = make_spacing(200, 0.75, 23.25, 34, 0.01, Om0) - r = np.arange(0, 200, 0.75) + r = np.arange(0, 200, 0.5) # smooth_scales = [0, 2, 4, 6, 8] - smooth_scales = [0] + smooth_scales = [0, 8] print(f"Running catalogue {args.catalogue} for simulation {args.simname} " f"with {len(r)} radial points.") diff --git a/scripts/field_prop/field_los.sh b/scripts/field_prop/field_los.sh index dc5efb7..f52d63f 100755 --- a/scripts/field_prop/field_los.sh +++ b/scripts/field_prop/field_los.sh @@ -1,5 +1,5 @@ nthreads=1 -memory=48 +memory=32 on_login=${1} queue="berg" env="/mnt/users/rstiskalek/csiborgtools/venv_csiborg/bin/python" @@ -18,11 +18,8 @@ then fi -# for simname in "csiborg1" "csiborg2_main" "csiborg2X" "Lilow2024" "Carrick2015" "CF4"; do -for simname in "Carrick2015"; do - # for catalogue in "2MTF" "SFI_gals" "CF4_TFR"; do - # for catalogue in "Foundation" "2MTF" "SFI_gals" "CF4_TFR"; do - for catalogue in "CF4_TFR"; do +for simname in "csiborg1" "csiborg2_main" "csiborg2X" "Lilow2024" "Carrick2015" "CF4" "manticore_2MPP_N128_DES_V1"; do + for catalogue in "LOSS" "Foundation" "2MTF" "SFI_gals" "CF4_TFR"; do pythoncm="$env $file --catalogue $catalogue --nsims $nsims --simname $simname --MAS $MAS --grid $grid" if [ $on_login -eq 1 ]; then echo $pythoncm diff --git a/scripts/flow/flow_validation.sh b/scripts/flow/flow_validation.sh index 658d28e..8a897cf 100755 --- a/scripts/flow/flow_validation.sh +++ b/scripts/flow/flow_validation.sh @@ -37,15 +37,17 @@ else fi -for simname in "IndranilVoid_exp" "IndranilVoid_gauss" "IndranilVoid_mb"; do +# for simname in "IndranilVoid_exp" "IndranilVoid_gauss" "IndranilVoid_mb"; do +for simname in "Carrick2015"; do # for simname in "no_field"; do - # for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do + # for catalogue in "LOSS" "Foundation" "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do + for catalogue in "LOSS"; do # for catalogue in "CF4_TFR_i" "CF4_TFR_w1"; do - for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do - # for ksim in "none"; do + # for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do + for ksim in "none"; do # for ksim in 0; do # for ksim in $(seq 0 5 500); do - for ksim in "0_100_5" "100_200_5" "200_300_5" "300_400_5" "400_500_5"; do + # for ksim in "0_100_5" "100_200_5" "200_300_5" "300_400_5" "400_500_5"; do # for ksim in {0..500}; do for ksmooth in 0; do pythoncm="$env $file --catalogue $catalogue --simname $simname --ksim $ksim --ksmooth $ksmooth --ndevice $ndevice --device $device"