{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Optimal Signal to Noise Ratio (SNR) generation and comparison \n", "\n", "This notebook is a full guide on how to use the `gwsnr` package to generate the 'optimal-SNR' ($\\rho_{\\rm opt}$) and 'probability of detection' (Pdet) for a given Gravitational Wave (GW) signal.\n", "\n", "## Contents of this notebook\n", "\n", "$\\rho_{\\rm opt}$ calculation with,\n", " - Noise-Weighted Inner Product \n", " - Partial Scaling Interpolation\n", " - ANN-model and $P_{\\rm det}$ Estimation \n", " - Hybrid SNR Recalculation for $P_{\\rm det}$ Estimation\n", " - JAX assisted Inner product with `ripplegw` as backend\n", "\n", "Note: For more details on SNR calculation methods, please refer to the [gwsnr documentation](https://gwsnr.hemantaph.com).\n", "\n", "\n", "## Requirements\n", "\n", "- `gwsnr` for all SNR calculations\n", "- `jax` and `jaxlib` for Partial Scaling Interpolation, if you want to use the jax version. `\"jax[cuda12]\"` for running on Nvidia GPU.\n", "- `mlx` for Partial Scaling Interpolation, if you want to use Apple Silicon M-series GPU.\n", "- `ripplegw` for JAX assisted Inner product method.\n", "- `scikit-learn` and `tensorflow` for ANN model based SNR calculation. `ml-dtypes` might need to be updated to the latest version for compatibility with TensorFlow.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Below shows GWSNR initialization with all of its (default) arguments\n", "\n", "```python\n", "gwsnr = gwsnr.GWSNR(\n", " #################################\n", " # General settings\n", " npool=4, # Number of processors for parallel processing. Run this to check the number of cores in your machine; import os; os.cpu_count()\n", " snr_method='interpolation_no_spins', # SNR calculation method. Other options: 'interpolation_no_spins_jax', 'interpolation_no_spins_mlx', 'interpolation_aligned_spins', 'interpolation_aligned_spins_jax', 'interpolation_aligned_spins_mlx', 'inner_product', 'inner_product_jax', 'ann' \n", " snr_type='optimal_snr', # Type of SNR to be calculated. 'matched_filter_snr' option will be available in future releases.\n", " gwsnr_verbose=True, # If True, it will print all gwsnr settings\n", " multiprocessing_verbose=True, # If True, it will show progress bar for multiprocessing. \n", " pdet_kwargs=dict(snr_th=10.0, snr_th_net=10.0, pdet_type='boolean', distribution_type='noncentral_chi2'), # Dictionary of Pdet settings\n", " #################################\n", " # Settings for interpolation grid\n", " mtot_min=2*4.98, # Minimum total mass (Mo) for interpolation grid. 4.98 Mo is the minimum component mass of BBH systems in GWTC-3\n", " mtot_max=2*112.5+10.0, # Maximum total mass (Mo) for interpolation grid. 112.5 Mo is the maximum component mass of BBH systems in GWTC-3. 10.0 Mo is added to avoid edge effects.\n", " ratio_min=0.1, # Minimum mass ratio for interpolation grid\n", " ratio_max=1.0, # Maximum mass ratio for interpolation grid\n", " spin_max=0.99, # Maximum spin magnitude for interpolation grid. Note: spin_min= -spin_max, as in aligned spin systems.\n", " mtot_resolution=200, # Number of points in total mass axis for interpolation grid\n", " ratio_resolution=20, # Number of points in mass ratio axis for interpolation grid\n", " spin_resolution=10, # Number of points in spin magnitude axis for interpolation grid.\n", " batch_size_interpolation=1000000, # Number of samples to be processed in each batch for interpolation method\n", " interpolator_dir='./interpolator_pickle', # Directory to save/load the interpolator\n", " create_new_interpolator=False, # If True, it will overwrite the existing interpolator\n", " #################################\n", " # GW signal settings\n", " sampling_frequency=2048.0, # Sampling frequency in Hz\n", " waveform_approximant='IMRPhenomD', # Frequency domain waveform approximant of the GW signal\n", " frequency_domain_source_model='lal_binary_black_hole', # Source model for frequency domain waveform generation. \n", " minimum_frequency=20.0, # Minimum frequency of the waveform in Hz\n", " reference_frequency=20.0, # Reference frequency for spin\n", " duration_max=None, # Maximum duration of the waveform in seconds. Bilby default for IMRPhenomXPHM is 64 seconds.\n", " duration_min=None, # Minimum duration of the waveform in seconds. Bilby default is 4 seconds.\n", " fixed_duration=None, # If a float value is provided, all waveforms will be generated with this fixed duration (in seconds).\n", " mtot_cut=False, # If True, SNR=0 for total mass associated with signal duration < 1.1*chirp_duration\n", " #################################\n", " # Detector settings\n", " psds= {'L1':'aLIGO_O4_high_asd.txt','H1':'aLIGO_O4_high_asd.txt', 'V1':'AdV_asd.txt', 'K1':'KAGRA_design_asd.txt'}, # Power spectral density of the detectors. Other options: psd names from pycbc (e.g. 'aLIGODesign'), psd via gps time (e.eg. 1234567890), or custom PSD as txt file.\n", " ifos=['L1', 'H1', 'V1'], # List of detectors. You can also provide bilby interferometer objects.\n", " #################################\n", " # ANN settings\n", " ann_path_dict=None, # Path to the ANN model for SNR and Pdet calculation\n", " #################################\n", " # Hybrid SNR recalculation settings\n", " snr_recalculation=False, # If True, enables hybrid SNR recalculation for systems near detection threshold. Default is False.\n", " snr_recalculation_range=[6,14], # SNR range [min, max] for triggering recalculation with inner product method.\n", " snr_recalculation_waveform_approximant='IMRPhenomXPHM', # Waveform approximant to use for SNR recalculation. Default is 'IMRPhenomXPHM'.\n", ")\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Noise-Weighted Inner Product " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "Copying interpolator data from the library resource /Users/phurailatpamhemantakumar/anaconda3/envs/gwsnr3/lib/python3.11/site-packages/gwsnr/core/interpolator_pickle to the current working directory.\n", "psds not given. Choosing bilby's default psds\n", "Intel processor has trouble allocating memory when the data is huge. So, by default for IMRPhenomXPHM, duration_max = 64.0. Otherwise, set to some max value like duration_max = 600.0 (10 mins)\n", "\n", "Chosen GWSNR initialization parameters:\n", "\n", "npool: 4\n", "snr type: inner_product\n", "waveform approximant: IMRPhenomXPHM\n", "sampling frequency: 2048.0\n", "minimum frequency (fmin): 20.0\n", "reference frequency (f_ref): 20.0\n", "mtot=mass1+mass2\n", "min(mtot): 9.96\n", "max(mtot) (with the given fmin=20.0): 235.0\n", "detectors: ['L1', 'H1', 'V1']\n", "psds: [[array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,\n", " 4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,\n", " 6.51153524e-46, 6.43165104e-46, 6.55252996e-46],\n", " shape=(2736,)), ], [array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,\n", " 4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,\n", " 6.51153524e-46, 6.43165104e-46, 6.55252996e-46],\n", " shape=(2736,)), ], [array([ 10. , 10.02306 , 10.046173, ...,\n", " 9954.0389 , 9976.993 , 10000. ], shape=(3000,)), array([1.22674387e-42, 1.20400299e-42, 1.18169466e-42, ...,\n", " 1.51304203e-43, 1.52010157e-43, 1.52719372e-43],\n", " shape=(3000,)), ]]\n", "\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|█████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 3.30it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Computed SNRs with inner product:\n", " {'L1': array([46.53022776]), 'H1': array([48.20334347]), 'V1': array([13.20663045]), 'snr_net': array([68.28645184])}\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# loading GWSNR class from the gwsnr package\n", "import gwsnr\n", "import numpy as np\n", "\n", "# initializing the GWSNR class with inner product as the signal-to-noise ratio type\n", "gwsnr = gwsnr.GWSNR(\n", " snr_method='inner_product', \n", " waveform_approximant='IMRPhenomXPHM')\n", "\n", "# signal-to-noise ratio for a BBH with GW150914 like parameters with detectors LIGO-Hanford, LIGO-Livingston, and Virgo with O4 observing run sensitivity\n", "snrs = gwsnr.optimal_snr(\n", " mass_1=np.array([36.0]), # mass of the primary black hole in solar masses\n", " mass_2=np.array([29.0]), # mass of the secondary black hole in solar masses\n", " luminosity_distance=np.array([440.0]), # luminosity distance in Mpc\n", " theta_jn=np.array([1.0]), # inclination angle in radians\n", " ra=np.array([3.435]), # right ascension in radians\n", " dec=np.array([-0.408]), # declination in radians\n", " psi=np.array([0.0]), # polarization angle in radians\n", " geocent_time=np.array([1126259462.4]), # geocentric time in GPS seconds\n", " a_1=np.array([0.3]), # dimensionless spin of the primary black hole\n", " a_2=np.array([0.2]), # dimensionless spin of the secondary black hole\n", " tilt_1=np.array([0.5]), # tilt angle of the primary black hole in radians\n", " tilt_2=np.array([0.8]), # tilt angle of the secondary black hole in radians\n", " phi_12=np.array([0.0]), # Relative angle between the primary and secondary spin of the binary in radians\n", " phi_jl=np.array([0.0]), # Angle between the total angular momentum and the orbital angular momentum in radians\n", ")\n", "\n", "print('\\nComputed SNRs with inner product:\\n', snrs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Partial Scaling Interpolation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "psds not given. Choosing bilby's default psds\n", "Interpolator will be generated for L1 detector at ./interpolator_pickle/L1/partialSNR_dict_4.pickle\n", "Interpolator will be generated for H1 detector at ./interpolator_pickle/H1/partialSNR_dict_4.pickle\n", "Interpolator will be generated for V1 detector at ./interpolator_pickle/V1/partialSNR_dict_4.pickle\n", "Please be patient while the interpolator is generated\n", "Generating interpolator for ['L1', 'H1', 'V1'] detectors\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|█████████████████████████████████████████████████████████| 4000/4000 [00:01<00:00, 2797.16it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Saving Partial-SNR for L1 detector with shape (20, 200)\n", "\n", "Saving Partial-SNR for H1 detector with shape (20, 200)\n", "\n", "Saving Partial-SNR for V1 detector with shape (20, 200)\n", "\n", "\n", "Interpolation results: \n", "{'L1': array([ 7.40707612, 11.5964448 , 31.03520024, 26.59768336]), 'H1': array([ 4.7108495 , 7.37525919, 19.73817404, 16.91594381]), 'V1': array([2.22222811, 3.44022184, 9.33438517, 7.86082025]), 'snr_net': array([ 9.05511885, 14.16711355, 37.94614493, 32.48658816])}\n", "\n", " Inner product results: \n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|█████████████████████████████████████████████████████████████████| 4/4 [00:00<00:00, 14.93it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "{'L1': array([ 7.40686229, 11.5956576 , 31.03469553, 26.59785488]), 'H1': array([ 4.71071351, 7.37475854, 19.73785306, 16.91605289]), 'V1': array([2.22218338, 3.44021659, 9.33421055, 7.8608819 ]), 'snr_net': array([ 9.05486221, 14.16620729, 37.94552222, 32.4868003 ])}\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "import numpy as np\n", "import gwsnr\n", "\n", "# initializing the GWSNR class with default configuration and interpolation method\n", "# for non-spinning IMRPhenomD waveform \n", "gwsnr_no_spins = gwsnr.GWSNR(\n", " npool=4,\n", " snr_method='interpolation_no_spins', # Other options: 'interpolation_no_spins_jax', 'interpolation_no_spins_mlx', 'interpolation_aligned_spins', 'interpolation_aligned_spins_jax', 'interpolation_aligned_spins_mlx'\n", " gwsnr_verbose=False,\n", " waveform_approximant='IMRPhenomD',\n", ")\n", "\n", "# Quick test and comparison between interpolation and inner product methods\n", "mass_1 = np.array([5, 10.,50.,100.])\n", "ratio = np.array([1, 0.8,0.5,0.2])\n", "dl = 1000\n", "print('Interpolation results: ')\n", "print(gwsnr_no_spins.optimal_snr(mass_1=mass_1, mass_2=mass_1*ratio, luminosity_distance=dl))\n", "print('\\n Inner product results: ')\n", "print(gwsnr_no_spins.optimal_snr_with_inner_product(mass_1=mass_1, mass_2=mass_1*ratio, luminosity_distance=dl))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testing speed and accuracy between the interpolation method and inner product method.\n", "\n", "* set up the BBH Binary-Black-Holes parameters\n", "* simulate 50000 events\n", "* Consider IMRPhenomD waveform with aligned-spins" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "Copying interpolator data from the library resource /Users/phurailatpamhemantakumar/anaconda3/envs/gwsnr2/lib/python3.10/site-packages/gwsnr/core/interpolator_pickle to the current working directory.\n", "psds not given. Choosing bilby's default psds\n", "Interpolator will be loaded for L1 detector from ./interpolator_pickle/L1/partialSNR_dict_0.pickle\n", "Interpolator will be loaded for H1 detector from ./interpolator_pickle/H1/partialSNR_dict_0.pickle\n", "Interpolator will be loaded for V1 detector from ./interpolator_pickle/V1/partialSNR_dict_0.pickle\n", "\n", "Chosen GWSNR initialization parameters:\n", "\n", "npool: 4\n", "snr type: interpolation_aligned_spins\n", "waveform approximant: IMRPhenomD\n", "sampling frequency: 2048.0\n", "minimum frequency (fmin): 20.0\n", "reference frequency (f_ref): 20.0\n", "mtot=mass1+mass2\n", "min(mtot): 9.96\n", "max(mtot) (with the given fmin=20.0): 235.0\n", "detectors: ['L1', 'H1', 'V1']\n", "psds: [[array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,\n", " 4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,\n", " 6.51153524e-46, 6.43165104e-46, 6.55252996e-46],\n", " shape=(2736,)), ], [array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,\n", " 4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,\n", " 6.51153524e-46, 6.43165104e-46, 6.55252996e-46],\n", " shape=(2736,)), ], [array([ 10. , 10.02306 , 10.046173, ...,\n", " 9954.0389 , 9976.993 , 10000. ], shape=(3000,)), array([1.22674387e-42, 1.20400299e-42, 1.18169466e-42, ...,\n", " 1.51304203e-43, 1.52010157e-43, 1.52719372e-43],\n", " shape=(3000,)), ]]\n", "\n", "\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from datetime import datetime\n", "import gwsnr\n", "\n", "gwsnr_aligned_spins = gwsnr.GWSNR(\n", " npool=4,\n", " snr_method='interpolation_aligned_spins', \n", " gwsnr_verbose=True,\n", " # waveform_approximant='IMRPhenomD', # default waveform\n", ")" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# gerneral case, random parameters\n", "# add random seed for reproducibility\n", "np.random.seed(42)\n", "nsamples = 10000\n", "mtot = np.random.uniform(2*4.98, 2*112.5,nsamples)\n", "mass_ratio = np.random.uniform(0.2,1,size=nsamples)\n", "param_dict = dict(\n", " # convert to component masses\n", " mass_1 = mtot / (1 + mass_ratio),\n", " mass_2 = mtot * mass_ratio / (1 + mass_ratio),\n", " # Fix luminosity distance\n", " luminosity_distance = 500*np.ones(nsamples),\n", " # Randomly sample everything else:\n", " theta_jn = np.random.uniform(0,2*np.pi, size=nsamples),\n", " ra = np.random.uniform(0,2*np.pi, size=nsamples), \n", " dec = np.random.uniform(-np.pi/2,np.pi/2, size=nsamples), \n", " psi = np.random.uniform(0,2*np.pi, size=nsamples),\n", " phase = np.random.uniform(0,2*np.pi, size=nsamples),\n", " geocent_time = 1246527224.169434*np.ones(nsamples),\n", ")" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duration interpolation: 0:00:00.035357\n" ] } ], "source": [ "# for non-spinning IMRPhenomD waveform \n", "# running it for the first time will take longer, as the jit code generate the compiled code\n", "start_time = datetime.now()\n", "interp_snr_aligned_spins = gwsnr_aligned_spins.optimal_snr(gw_param_dict=param_dict)\n", "end_time = datetime.now()\n", "print('Duration interpolation: {}'.format(end_time - start_time))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|███████████████████████████████████████████████████████| 10000/10000 [00:03<00:00, 2813.75it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Duration inner-product: 0:00:03.628216\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# error might occur for mchirp>95., if f_min=20. \n", "start_time = datetime.now()\n", "bilby_snr = gwsnr_aligned_spins.optimal_snr_with_inner_product(gw_param_dict=param_dict)\n", "end_time = datetime.now()\n", "print('Duration inner-product: {}'.format(end_time - start_time))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* interpolation method is much faster than inner product method (multiprocessing)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Accuracy check plot" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAGACAYAAACtE/1YAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAArrpJREFUeJztnQd4FFUXhk86JCEJPYTee1ekd5SOIAhKkV6EH6RKJ3REqYKgdKTYQAERFBCkSC/Si/TeCQRI3/85FyfsbnY223d25nufZ7K7M7N378xuZu53T/PS6XQ6AgAAAAAAAACF4e3uDgAAAAAAAACAKSBWAAAAAAAAAIoEYgUAAAAAAACgSCBWAAAAAAAAAIoEYgUAAAAAAACgSCBWAAAAAAAAAIoEYgUAAAAAAACgSCBWAAAAAAAAAIoEYgUAAAAAAACgSCBWAAAAAAAAAIoEYgUAAIAm0Ol09P3339Pbb79NGTNmJF9fX8qSJQs1atSIfv31V4N9IyMjycvLi8qUKUMJCQmybWbKlIny5Mlj8r2mFn9/fypUqBCNHTuWYmNj7TqeK1euiDY7duyYvI6f8zreJnHq1CmqWbMmBQcHU968eVNdD16DcwSA+/F1dwcAAAAAZ5OUlETt2rWj1atXU7Zs2ahJkyYUEhJCly5dos2bN9Nvv/1Gw4cPp4kTJxq8759//qEvv/yS+vfvb/VnNmzYkEqWLGmw7tGjR/THH38IQbNv3z7xuSwunEmnTp3o4MGD1KFDB8qcOXOq6wEp6hytWLGC2rdvL54vWLCAunbt6pZ+AG3y008/0V9//UXHjh0T18Nnz55R27Ztxe/SFA8fPqSff/6ZNm7cSCdOnKCbN2+KSRq+FvL/Ey/e3lbaSnQAAACAylm+fLmOb3lNmjTRxcXFGWw7d+6cLmvWrGL7wYMHxboxY8aI17wEBwfrrl+/brLdjBkz6nLnzm2wTnrvkiVLTL4nKipKV7BgQbHPtm3bbD6my5cvizY++uij5HW3bt3SnTlzxuAYfX19dRUqVEjxfrn1QDnn6Nq1a7rQ0FDxG+TvesGCBW7rC9AmpUuXTr4OFilSRDxv27at7P7z5s0T+2TLlk334Ycf6oYOHarr1KmT+B3z+vfee0+XlJRkVR/gBgYAAED1rF+/XjyyRcPPz89gG7tljR8/Xjxfs2aNwTaeUY+Ojqa+ffs6rC9s0fnoo4/Ec561dyRsNSpSpIjBMbIbW9q0aVPsK7ceKOMcsdsiz0Kzy2LPnj2d9jn8++YZc6A9oi347mfMmEHnz5+np0+f0rx581Jtk6+nfL29ceMGrVy5kiZPnkyLFy+ms2fPUs6cOcU1du3atVb1E2IFAACA6mHXBf1HY9gtjAVL5cqVDda3bNmSGjRoINwajONa7IEHoEyaNGlS3ffx48fCDS1XrlwUEBAgBgPjxo2jmJiYFPvqx6xwrIXkYsZuHFJ8i9x6CW6XBxgsevjzsmfPTh9//LFwYdOH2+GYHXaxmzZtmtj//ffft6odqQ0eCPXu3ZvCw8OFywi/h12ejOHvb8SIEZQ/f37RJp+TZs2amRR9lh6HKcydoxcvXgh3wWLFionvL3369CLuadeuXSbbMXeOzDF79mz6888/acmSJRQUFESOFmHsgsjuPFmzZhUDUkvZvXu3OC4W3Rz3lTt3bpowYYI4RqB8Eqz87mvVqkUFCxa02F21du3a4npq7OrF/9uS6N6xY4dVfYZYAQAAoHoqVqwoHjt37kybNm2i+Pj4FDfSkSNHipusPnyDnjNnjhiU9unTh54/f+6Q/uzcuVM8litXzux+7P9dqVIlmjlzppiV7NWrF5UtW1YMwvUFhik+/PBD+vTTT8Vzfi8/50G13HppgM8JCDh+h2M0unTpQvny5RMzqizknjx5kuJz2OrE+/NA/J133rG6ncTERGrcuLGYjX333XdFfAbPynbv3l0kRJCIioqiatWq0aRJk4QA4IFPlSpVRMwRD54PHz6cvK8tx2HJuePvnwdv/Fvh3wb3kWOTtm/fTjVq1BDCwhSmzpE5zpw5Q0OHDqV+/fpR9erVyVHs3btX/I7ZAsfHw+e3atWq4ngtgYUi94fPH5//bt26CaE5atSoFPFeQFnstfO7dwSSxZdFrlU4wJ0NAAAAUDTR0dG66tWrJ8ehhISE6Bo1aqSbOHGi7tChQ7rExESTcScbNmwQrydNmiRef/rppzbHrMTGxuouXrwofLh5e+3atVP13e7evbvYt3///gb7cmxNhgwZUsSs8HNex/EsEvy6Ro0aKdo2tX7EiBFi/ezZsw3Wz58/X6z/3//+l7yO3+vl5aXLkiWL7vTp0za1w23w6zfffFP3+PHj5P3+/vtvsb5hw4bJ6/r16yfW8aP+udi3b5/Ox8dHV758eZuOwxzG52jgwIFiXZcuXXQJCQkGcU/8fQQEBOiuXLli0TmSIz4+XhxLoUKFdC9evDD4TdkSs3L27FndqFGjdPny5Uv+/VesWFGcm7t371rczvTp08V7hwwZYnD+jxw5ItaHhYWl+D8C7uWsg757ie3bt6cas2Lud12iRAnx/s2bN1v1XogVAAAAmoAHUuvWrdM1a9ZMiBXp5s1L/vz5dUuXLk0ehBmLFRYaxYoVEwHXJ06csEismFu6du2qe/Lkidn+xsTE6NKmTavLlCmTeG7MhAkTHCpWePCdPn16XZkyZVKIKH7NgoKDZHnQoS80vvzyS4N9rWlHamPr1q0p+pc5c2Zd8eLFxXPel78zPt8vX75Mse/YsWN1HTt2FIkFrD0Oc+ifI94/Xbp0og+mvruZM2eK/bkvEnLnyBw8uPT29haCTcJascKJFlhcsOiRfnM8UGRxfunSJZ213Lx5UxcUFKSrVKmSSUEiBV7rCzXgHm45+Lt3lFiRhL7+BISlIHUxAAAATcA+1E2bNhUL+22fPHlSuGOtW7dOxAawWxUHkppyZ+E4ivnz5wsXGHZ/4felln7TVOpidp9h94vly5eLWJgWLVrIvv/cuXP08uVL4SLFMRem2md3JEfBAbAcH5MjRw4aNmxYiu3sWsWuWJzumeNmJOrWrWtzOxJvvPFGiv0CAwNFADBz4cIF4W7UqlUrk3E+o0ePNqiNYstxpAZ/Hxwz88EHH1BoaGiK7Rzf9Mknn4iU1MYYnyM59u/fL9zcBg4cKNz/bIVd3ThuieNp2IWNXX1KlSplc3uc8ptd4AYMGGDydy/F1Pj4+JAjYfe+EiVKCFdMV2Lt57qrn6747h0Bx19JMVvffvut1e+HWAEAAKBqOPA3Li5ODKT0faa54CMvHE/AQdQc4zB16lQxIDMFx0twzAtntuHYBI6DMAcPrE3FlUhxJ//73/9EcLjcAI8H5wwXrjQFx9k4Eh7gM1wbgRc5OJ5Ef5DPgeu2tiNhavBvHLtj7lw44jgsbVeu3orUN2k/fYzPkSlYQHP2Of3sdLbCIpkHrNyX33//nTJkyCDOMQfD2wILehYpLJBNwbU0WLBERETY1W+1w9chnvhwJiUd/N3bCws4jr3ihBTbtm0T/bEWBNgDAABQNUeOHBHpZzkTlBwcHM2WDh4w8iy+HCxmOJPXkCFD6P79+zb1p3Tp0sJCc+vWLbp8+bLsfjwzyty7d8/kdrn1tiIJBhZk/7mJm1w4248+6dKlc0g75pBm7iXRYgxbPB48eCCC9Z3x+frHxZ9jCun3EBYWlmKb8TkyBVuR2LLHwfVsPeIAfmkZO3as2IcD2vk1W3DMwckKLl68KLLGsXWOZ9jz5s0rAqrnzp1r1W+H/ycOHDggkhqwtcsYzsR2584dkXhAsrqwpYH/3zipAL+PhdygQYMMMobxd8D/T5zZjf8/eZCtX2iQhT5PInB/pfPAg3BzxMbGinPDWa74HHJiDc5epg/3ja2jPIDm/zFeBg8enNw3az/X3P78WTw5wcfOIpcTQkjrOdjduB22olpyblz13dsLJwbhSRm2OnESClsnWCBWAAAAqBqeqebZTK4cby6bl7RNEgmmYKHyxRdfiPS3PMCxleLFi4vHu3fvyu5TuHBhMcjlmzynyzWG0486Ev48HoyyG5OpNLScSYsHoDw77Ip29GH3Ef4OOT2wcSY3hrMb8YCMB9bO+HypDyya+HdkKgW2VDuiQoUKZAvs6sfWOlMLW+IYHnDya0tcxDj7GWfpYre8Q4cOiUE8u77xIJktIJyVbOnSpcIlzhzsVscigM+XqfPJ7j0MZ0bTh2tssAXz77//FrPrPHDVz+7GLoyLFi0SA+jTp08Ll70ePXqIyufMrFmzxHFyrZnbt2+LhTOzmYMnEfgz2Pp59OhRMcivX7++eK9x3/hYOEPW119/Td98843ony2fm9r+LDJYfPBvl90/LSG1c+Oq794ePvvsM5Fyna3XfA2zxCoqi9VRLgAAAICH0a1bNxHc2bhxY92DBw9SbF+5cqUIai5btqwIwjYOsNeHt0uZxTjTk7UV7KUgat7nt99+M9vvvn37iv169OhhkH2KMzBxoLejs4Hx5/B6zn6mH5zOWYU4wLpWrVopgsdNYWk75trg86p/bjmol/fl86vfJp9DXv/+++/bdBzmMD5H0vfRuXNnEcwvwUkXOBGCn5+fQRCzueOzBnuygenDv6E//vhD/E44WQC3yRnMpkyZIvueRYsWJQdqb9y40WDbnDlzxPr69esbrOfj5qxT+tStW1dkUZOy86VJk0a3c+dOg30401uDBg0M2undu7dFx8Zt8vlftmyZwfFyJizODqffZsGCBQ1+F+PHj9dlz57dps81tz+vL1mypEX783fCGQotPTeu+O5tDbAfN26c2I+D/B8+fKizF8SsAAAAUD08+8szi1zYkQsJcnwKP7LFgoOaObaBZ0NXrVqVavEz3s7B9uzOZWqW3xKkOgPm3MAYdv9hP2+e/WVXHHYfYzckLlLJFiM5tyhb4eBudmlh950NGzYIKwF/3i+//CKO29LigY5qRx92i+HEBnxO+JGDhqXvlC1e06dPd+rnM1z8kN2KeOaefzd16tQRVjauys0uN/y7YAuPUuH4qHr16omFLUx8bnjmn92G5JDq13ANIk4iwAHb7NLEloI9e/aI5Aj8f2OMcVA3z+hLLkhsLeBEB2z10P9/4/+nPHny2HRsfAz8fsnVSjpetnrw5+nD7mH6n8v7sCWC48S42KUjKV++vFX7O+Pc2PrdM/w/wwvD7n4MW6SkeDx282Nrs8SyZctEwgv+PI7z4+B6Y/g4UqsTpQ/ECgAAANXD7lTsvvPdd98JVwweZPHNmge5PLjkgGZ2kTAVb2CKokWLCpcTWwvhcbAp89NPP5mNpeH+8OCY+8f7fvXVV8KdomvXrsIdiAWTI+HgV3bb4eNitybp89jNigcgkvuaq9oxHuxyfASfC/bL58Eynx+uCM/r9IPYnfH50u+IhRIPzjhDFosTjilgEcmxARy34Slwv/nc8cKxPubECrvgsXsVx0Hw/w8LP/6/YfHGmctMZWiTkllI8MBbciOTHvl/kCcNzL3PEVhafd0ZSPFW+nBszyvD3WukiQ9XnJu0Fn73zLFjx4QA0YcnCaRsfhy4ry9WpAkYbldyrTMVI2iNWPFi84rFewMAAAAAAE3AA04WaCzOJQuLral8eXDKIoctYRzzw9YZnuHnWA852ALKQea8X2pwzBmL1AULFoisalL/2QLJ1iApwxr3jeNKOJ5DEjEsalnQclYzaz/X3P5yKY1bt24txIkU58TwfmxxYBFsybnRErCsAAAAAAAAky5J7N5Wrlw5h7bLAogzZPHCc+ZsmeJsaJwUga0OUrA+D97Z/ZGzawUHBwsxIlffiC0YnHmLLVzsmsSWH3b34yQWxtZLzsTHQee8nl1AP//8c4OaRdZ8ri37cyY6/ny2EHJCCHbzvH79umjH0nOjJZANDAAAAAAApECypkjZyBwJWzoiIyOFCxG75XEsBcf+6Mf88ICdXdDYbZKtDdeuXRNZrOTSCXMGKrZasEWCs1AdP36cNm/eTNmyZTPYr23btsLq8tZbb4l00OxSyZmrbP1cU/ubg9NqSwvH2LBAad68uVXnRkvADQwAAAAAAKSAa2SwCxMHVHNQuhIYM2aMiN/6559/khNVWIOt1ebt/VxgO7CsAAAAAACAFHz55ZfCFUkpQkWqL8T1R1wtGNz1uQAxKwAAAAAAwEPgjHBa+lwAy4qmYX9I9r80tYSGhorsFpwy05WwLyh/vn5KO34u5yeqBDggzp7c5644n5YuO3bsIKXAtQs47z3n6uc0i1zduWDBgtS7d2+6evVqiv35O+BjkKopm4Kz0PA+/Nu35P+AfZA5kwzXdeAKzgAAADwbvs9Z6wIG3AssK4AaNmxIJUuWTH7NRZE4r/aWLVvozz//FLUJOHOFpyIJCaWKHWfCopMzo+jDmU/YnM1BfVWrVjXYZpzT3V38+++/wq+Y00hWrlxZBCHGxcWJImycXpILWXGhPC5GZsqvuFWrVlYfi/H/gSSY+PfPgoYzsfB5c2e+fgAAAEBrQKwAMbAzLs7DPqo889C3b18aNmyYGCS6i8mTJ9PQoUMNCn4By0ifPj1NmTLFYB1nNOFBd926dQ0sDEqiR48eQqhwAT/O7GLqd8nZXjjTi7F44Fz7/fr1ExW+7f0/kMQ7iyLOKLN9+3aPFu4AAACApwE3MGASHgCyu03WrFlF7nDO8e0uOOVgkSJFnFLVFigPLhbGFj22cugLFel3ydlp2Npy8uRJOn/+vMF2zm3PFpJffvlF5K93BCEhIfTRRx+J5/BZBgAAAFwLxAqQhYsPSa40T548SXap4llm9t/nAkqc85vzi0uz3t9//z1VqlSJAgMDKWPGjNSsWTPhdmTM48ePRU5zbp9jETguYNy4cRQTE5NiX7mYFfY7feeddygsLEzkKOec6hMmTEiOLZDex/ENvPBz49gSziHPfeQ2uJATZzxZt25dij4kJCTQrFmzRLrDNGnSCAHVtWtXUVjKEvic8fnk4lTGsOWD+8bxFMy9e/dEsSjuK58b/izOv/7333+TM+Hzyf1gywUXAnv33XdFYS12CZTiOkzFtMh9P1u3bqU6deqI74bd0fj5rl27Uu0HW0akx6SkJJP7sKsX56Hn70Ifqf8c48KiRmrLXvi3zBh/HgCeBl8jCxQoIOK/mjZtSg8fPnR3lwAAwCwQK0AWHihKg3x2J5JgUdKyZUuaPXs2VahQITnugWMj2rRpI/z827dvL6qubtq0SQgAqbAUwzdHFjQzZ86knDlzioqzXHCK3b1MueGYgl2ZatWqJXK/8w2X38f95YBsfs59bNSokegTz4zzws/5syR49p37tmfPHmrSpAm1aNGCLl68KAbp8+fPNzjeDz/8UAgIruTLIqVx48ZCXEiVZVPjgw8+EO2YEkK8jsUSJzTgwTUP6lkYcXEprq7LVoSNGzdSjRo16MiRI+Rs+Bzw93P27FnxPWfJksXqNji9IxexunDhAr3//vvUoEEDYaHj70wSZXLw57EIvnTpErVr1060YQyfKxbLuXPnTrGN38tihotycWC8I9i5c6d4dHQVZwBcCU8WrF69WlxH2CrJ/2vTp093d7cAAMA8XBQSaJMxY8ZwQVDdkiVLTG6fNWuW2P7WW28lr8udO7fO29tbV6hQId2tW7eS12/ZskXs27p1a11sbGzy+sOHD+sCAwN1pUqV0iUmJop13bt3F/v2799fl5SUlLzvwYMHdRkyZBDbPvroo+T1/JzXXb58Wby+efOmLm3atLqsWbPqTp06lbwft9+0aVOx74YNGwz6zIs+Dx480IWEhOiKFCmiu3PnTvL6e/fu6UqWLKkLCAgQn8OsWrVKtFmtWjXds2fPDPYtX7682GbcvjHXr1/XeXl56erXr2+wnj+b13fp0kW8/uWXX0R7ffv2NdhvzZo1Yn2fPn109sLfN7fF378+27dvF+v5++3Xr58uISEhxW+F9zHG+Ps5ffq0ztfXV1e9enXd06dPk/e7dOmSLjw8XHxvz58/N9vHnTt36kJDQ0W7vBQoUEDXqVMn3eLFiw1+d/rwd5AxY0bxnH+DxYsX1/n4+OiOHz+evA//LoyPXe7/gNu4ePGibujQoWJ77dq1DX6vADiLb775Rvaawr9B/q3y/wT/n/H/06effqp78eJFqu3++++/uv3794vnMTExuhYtWuimTZvm8P4DAIAjgVjRMNIgrWHDhuJmJy29evXSVa5cWWzjwd62bduS38M3UGMxwLBI4EEuiwBjBg8eLN7DYoRvkCw0MmXKJJ4bM2HChFTFyueffy5ef/HFFyne/88//4j9f/rpJ4M+G9/4p0+fLtr4+eefU7SxceNGsY0/h6lTp454zcLLmN27d1skVpgaNWro/Pz8dE+ePElet2DBAvF+FnvMunXrxOvGjRsbDD5YOPA5N9UHR4uVbNmyGQhOa8UKCy1+ffTo0RT7zp07V2z78ccfU+3no0ePdJMnT9a9+eab4rclCRdemjVrpjtx4oSsWGF27dol9uXfsiSUzYkVc0vXrl0NvjcALPn98iSDHDwJEh0dnfz65cuXuitXroj/keDgYNlryqJFi5Inhn744Qfd8OHDxf9Hhw4dLO7bihUrxGQNTzrxRAoAACgZZAMDIjMUL/qw2xS78bBbVbVq1VK8h12V9OF4Co5N+Pzzz1Psy/EPDLsecN0KdqViNyqOxzCGg6PZvcccUpAzuxYZwzU52EUsNaT4D05LyylpjbM/Sf2VHjNnzmzSBYhdtPhcWQK7kv3111/iXLNbGMNB4OyKwWl6GT7nRYsWFa5S+fLlEy5u7P7FGaj4nLkCdm3j78lW+Nyy6+CqVatENi99rl+/nnxO2cXMHOx6yFngeImKihIZ6TgO5qeffhKuc5xam+sAsQuhKdg9sUuXLrRo0SKxdOvWzerUxRyrxXFYy5cvF783dhUEwBI4xTa7I/74448i5kwf/j2yOym7S7KrJ8O/c3b/NI6T0ocnGDm2j3/b7M7F/2ecxY6ZNGmS2GbKNdKYtm3bCpdddpPk97M7LQAAKBZ3qyWgXDcwU/BsH7tqGcMWmNRmp0eNGpU82927d2+T7bOLT2qWFcnScf/+fYv7bDxLKbVhbmG3L+nY2KVIDp6dtMSywlYndtto1aqVeM2zqmnSpElxLngGf+LEibpy5coZ9IddyNiVytmWlQEDBqR4jzWWlfz586d6btu3by/bv7i4ODHLLFlDjImPjxfnjNtp3ry5rGVFOudsxUufPr3u7t27VrmBSRw7dky46kVERBi4xgFgDv6dssWZXUq3bt2avJ6tLXxNadeunYFbIf//7N27VyzsFmrqmsJuXPxbZRddY4syr1+6dKlu3rx5utKlS6dY2GrL7mX6Vk12dWU3XQAAUDIIsAdWwxYUYzjbE1sC/nMtNLnwrJ8UqM8Zr0wht16foKAg8Wgqi01iYiI9ePBApL81B/eX4SBsuf5KQdXcZ3P9sqTP0kxp/fr1RdIBznrGVgJ+5BlO474NHz5cJCXgttmSwMkCuM4Hz/7zMbr6+zUHWz2M+89ppjmDmty5ZUuFHJwljrN5yc32+vr6ioxG+lY7c+ecK9pz9rnBgweTLZQuXVpYmzjz2+XLl21qA2gP/p2yVY6TeHDSDrYMciFTtrCylW7JkiUGNYI4+x/vy0uOHDlMtnnnzh3xyNm89MmfP7945NpEPXv2FBn8jBe20vL/1WeffSauj/x/+O2334risAAAoGQgVoBDKF++vMgcZiqVL9+g+QZ65swZKly4sBgMc3G9Fy9epNjX2B1NbvDIcC0OY9jNil22OH1tav1lTA2I2VWJ+ytl7nrzzTfp/v37Jmts8PultM6WwO5fnD2MhQq3z4MSdiWTYLcwzjrGgw6Gj+W9996jDRs2CFcxztDljgGzJBDZhU8fFk5Hjx5NcW7j4+MNMsBJ8KCJzy27b8khuWJxtjY5pJTE+lnq5ODMdOxKxwKJf3e2ULx4cfFoKvU0AHJwqmv+n+bU7DzRwKLlrbfeoh9++EGIGWvh/yuGU8PrI72Oi4sz+34WSpxJj10nuU/8f8guaQAAoGQgVoBD6Ny5sxi4sv+1vgjhQeXAgQNpxYoVoqYK36C58jhbPwYMGGBgJeBBL8+CW+JvzTVLuD6Jfm0PtlKw9YZnK3kWUYJf8yy/cRs8+8/xMSyyJHi2ket3fP3118kDYSnWoU+fPiItswQ/50rp1sD94oEFW0s4LoWLHvKxSLDYY7/1efPmpRik8Gwo78uFD5nbt28L8cKPzkbyg+cUyvrnimOU9M8fw98vw+dG3/rF3/WwYcPEueWU1XJwqmMWaV988QUtWLAghSWJv2f+Lhgp9scc/P1zKmr+vjldti1I8QNSPBMAlsIxbfxb5usFT1Twb5AtHI5E30JjDr5+TJw4kf7991+REpwnQcz9LwIAgBJAgD1wCDzo/vnnn8WMIRdOZHcnHmSyyxNbKhYvXpw8O89BnWxt4UEr195gFxsWL/x+nu1LrUgZW2dYlLDQ4GKLHLzKN2EOuGarAw+IpZlwhge+bBXhgFZ2n+BAVh5886Chd+/ewlLDhSHZfYktJYcOHaIOHTokJxbg2VCu98GCiwP4uSYLwzd6HiSYCoSVgwtPsmDhtvj8GLuA8cwnCwAeULAVgI+PhQpbkc6dOydcmSSxwse5bNkyUV3dkqQC9sCzsXycHBDMgxyuAcNB8uzawkkB+NxLcI2WIUOG0NSpU8VvgRMDcMA+HwOLq9GjR5sNAmaRyL8FPs/du3cX7XCb/PkszDgpArt18bnj788SihQpIursSO5j1iLNgsMNDFgL1wviyRH+DbJlki2lXNuJi73a+juULCwSkkWFBTkAAKgNWFaAQ5CyP3355ZdiJpH9sdl6wHEsPBsvzbYzXACR3Q/YssLChLPmcGV0FhM8+LaEESNGiPY5cxZnnOLBOgsBHkwbD0i5+joPjrlP+sUe2QrEcSDstsQuR+wOwTd9Pgae0ZdmK/mR2+f1LBS4HfZF54KYLLr4c62BBQkLFfYzl9zRJPjccTYtHuxzvMo333wjjo8H6gsXLhT+5u6AvzMWJFzUkb877hcLRP7eWMAZw1avlStXCj98zlrE4oyPjX3k+ftIDfajZ2HDgpTdBteuXSt+JyxueRu/5t+bj4+PxcfAcUCSb7+1sDhj+DcHgKWwuGYxzyKC/394YWszi399K62lZM2aVTwaWzNZEDEREREO6jkAACgHL46yd3cnAAAAADXB1j+OleI4Jxb4BQsWTI7b4vgznmjh2DXJ4qwPC3qeINF3c2X4ds1uWxxzwrEw0oQKWx/ZcsjuXbYKcgAAUCqwrAAAAAAOhi2B7ALLbouSUGHKlCkj4tX++ecfiy3JEixOuPYVv58TVbBgmTx5snCtZLdICBUAgBqBZQUAAABwMJzUg+PM9OPn9Dl+/LjIfGcqOF7OssLwLZtdQtnVkt3B2DWVM95xYhDjLGEAAKAGPMqywnEE7ANv6uLNF3aevWLf4PDwcBFErZ9m1ZJ9AAAAAEfAwfByQoXhWC+5LF4sVkwJFYbfwxkKL168KAQRx7ZxFkUIFQCAWlG8WOE0pTx7xCZ1Dsg2BQc8cwA3Bytz0G2XLl1ERiU2k1uzDwAAAAAAAEA5KN4NjAvkcd0JCc7qpD/jxN3PmzevCDrkiuPSTBVni5o0aZLYl+t7pLaPuVSqAAAAAAAAANejeMsKixWufcELW0OM4ZSNbHlp1aqVgUmd634wnFrVkn0AAAAAAAAAykLxRSE5RkWKU+GaGMbcuXNHPHKxP32krCg3b960aB85YmNjxSKRlJQkCnJx7QdLqwYDAICSYIs01xTiGkFcLwfYD98buEI9F0DFvQEA4InoFHpvULxYSQ2pkq9xcKH0mk+6JfvIwWkhueK6Pg0aNKC6deuSlkiTJo2IH9ISWjxmrR63Fo/52bNn1L9/f1GsE9gPC5UZM2bIbud6KlwQEqQE58Y8OD/y4Nw459wo7d7g8WJFDktmtizZZ9iwYQaB/Tyg4axiXG09ICCAtML9+/cpc+bMpCW0eMxaPW61HXOJMb+nWJf12UOau24KFXlwlZ6kDab5n/xPWAGAY5DOJd/kTd0b1PYbcyQ4N+bB+ZEH58ax54Y9iXjSRWn3Bl81pIdkJOuJhGQt4TTFluwjB9909G88/JzbM16vdviHq6Xj1eoxa/W41XTMeYZuJCIfg3XF7l6ixT9FUnj0I7ofFEY9mg2nsvQM7koORDqXcvcGNf3GHA3OjXlwfuTBuXHOuVHavUE5Dmk2kjVrVvHIAfT6cFA9ExERYdE+AACgDqFiSM2Lh+iHVZ8KoXI+Yy5q3n4anQo3jN8DAAAAlIrHixUOms+ePTutW7dOBAZJrF+/XjxWr17don0AAMBTSUzSmRQqTI6ouxQc95J25y5NLdtNpRuhryZvAAAAAE/A493A2FQ1atQoUdyRl0aNGtGpU6dEUHybNm2SM35Zsg8AAHgaG/65Rf9bfVR2+4pyjehx2hD6o1BFiveRd3kFAAAAlIjHixWme/fuIsXalClTaNGiRZQhQwbq3bs3jR8/3qp9AADAk2j/zR7ademJwbqA+FgatOtbmlOpNUWlTSfWbSxazU09BAAAADQkViIjI8ViyrrSrVs3schhyT4AAOAp5Bu6kZKM1mV4EUUL14yjcrfOUeH7V6nD++P44kdK5dq1aw5rKywsTFGpNgEAwNHuvgcuP6J7z2IoS7o0VCFvBvLxVu71XbNixdNy7qupbsPjx48dUiCI61lwsSEAgO2Yik/J9/AGLfkpknI/uUNP0gTTnMqtFS1UGC7466isM2PGjKHRo0c7pC0AAFASm0/eprEbTtPtqNfjymyhaWhMk2JUv0Q2UjsQK04SKj/88AMlJCSQWkhMTCQfH8N0qLbAaZ/ff/99CBYAHChUKlw/Sd+snUBhMdF0NSycOrWMpEsZc8i2cXLsOzRlinyci6tgceEosVKjRg2HtAMAAEoTKr1WHKHX6aFecScqRqyf166c6gULxIoTYIsKC5VatWpR+vTpSQ1wjRpz9Wgstc5s375dnB+IFQCsIy4hiQqN3JRifbNT22nqplkUkJhARyIKU9f3RtOjwFDZdq5MaSQKfykBU269AAAAXrt+jd1wOoVQYXgdT/Xw9nrFwlXtEgax4kRYqGTKlInUgCPECgDANsZtOEWL91xJsZ6D6QfuWiGEysbCVWhAowEU62e6CBjfxi5PaeSC3gIAAHAEHKNyW8/1y5Rg4e28X6X8GUmtQKwAAICCeWv8Zrr7PNHkNhYm7PLV9MxfNLPqh6TzMh1XFuDrRecmNCS1uNkysM4CANQOB9M7cj9PxeOLQgIAgJrjU4yFSkhMNNW4dDj59cVMOWlGtXayQiVHaIAqhMrMmTMpV65cFBoaKpacOXPSjBkzDAr9AgCAmuCsX47cz1OBZQUAADwkkJ6r0S/+cSzleXyLPnp/LO3NXdpsGzPeL03Ny8kH2nsKQ4YMoW+++YYGDx5MlSpVEuv27t1L48aNo9u3b9PUqVPd3UUAAHA4nJ44W2gaEUyvk3HvDQ99lcZYzUCsAACABwiVkrcv0OI1Yynz8yd0OzgjPfmv4KMcFyc1VE3A5cKFC8XSsmXL5HW1a9emwoULU48ePSBWAACqhK/hY5oUE1m/+GquL1ikqztvV8u1Xg64gQEAgMKFSr0L++j71UOFUDmTOQ81bz+NzmTJZzbjl9puXqVKlTK5LinJuDQmAAAoO8PX3osPad2xm+KRX5ujfolsIj0xW1D04ddaSFvMwLKiBE6dkt9WvLgrewIAcBN8w8o//LcU6zsdWkejti0kb9LRjrzlqU+zTyk6INCsUFEbHTp0oLlz59KsWbMM1s+bN4/at2/vtn4BAIArijvWL5FNpCdGBXvgPkqUkN/moODRjz76iDZu3Cj8u/VTEHPQ6tChQ+nWrVuUIYN5n8c7d+5QyZIlqWDBgrRr1y5RJJKDW1u1akU7duyg48ePU0REhEP6C4CW+OnANRq09kSK9TUvHqIx2xaI5yvL1KfR9XpRorfp4qwRIX709/C3SY1wXZhVq1bR77//ThUrVhTr9u/fL65bbdu2pb59+ybvO3v2bDf2FAAAnFPc0cfbS9Xpic0BsaIR2K97+fLltGnTJmratGny+tWrV1Pz5s1TFSpMeHi48Bt/9913ac6cOdSvXz8xgFizZg2tX78eQgUAGyg6YiO9NJ2ZmHbkK09rStSmc5ly0zcVWhDJVHv/Z/TbFBqo3jpIZ8+epXLlyonnV69eTb4e8XLmzJnk/bxkzg8AALgTFHe0D4gVjVC5cmUqUaIErVixIlmsXLp0iQ4cOEATJ04UM5fmUoDyNrbINGvWjLp160bDhw+n8uXLU58+fejjjz+mJk2auPBoAFBvfErm6Mf0LCAtxfilEeJkYMP+siJFrW5fxmzfvt3dXQAAAJtBcUf7QIC9hujevTtt2LCBnj59Kl5/9913lCdPHqpTp47IqpM2bVrZZdKkScntTJ8+XVhRatasKR6/+OILNx4VAOoRKoXuX6Gfvx1AM3+dRt5J/5lbNC5U4uPj6a233qJz5865uysAAGATW0/fcXpxx0QrA/c9CVhWNAQHon766af0008/UefOnYULWKdOnYTrxM8//yysK3JkzZo1+TlXjm7QoAF9+eWXVKxYMSFmAAD2CZUqV47RvJ8nUUjcC4r19acML5/Sg6D0mhYqDFt0L1++DBcvAIDHxqos2nPFqcUdN5sI3K8c4UsdapdSRbYwiBUNERYWRu+//75wBeOZytOnT1PHjh3FtkKFCpl1A9MfKOzbt09k5uGaByx8fvjhB9EuAMA8cQlJVGjkphTrWx3fQpN+n0N+SYm0P0dx6t5iJEWZqaOiFaGinyBkwYIF9Pnnn7u7KwAAYHWsSmrYU9xxs0zg/uPncRYF7nsCECsaDLSvWrWquOnXq1ePcuXKJdYXL148OXDVFCNHjqTx48fTixcvRBpRLsjGIoXjX3r16kXVqlWjbNk8+58BAGcy+ucTtHz/NcOVOh0N2LWC+u79Xrz8pVgNGtLgE4rz9ZO9oV3WmFBhnj9/TitXrqQtW7aIWLmgoCCD7cgABgDwxFgVCZ2NxR0TUwncJ5UE7kOsaIxKlSoJ161ly5YJsSFhqRsYu5HduHGDNm/eLKwtXOeA2+Oge46HgasGACkpOHwjxZuoXThm2zfU6fAG8Xx2pdY0vVo72RiV7KH+tGdYPdIinPFLygbGiUH0wTUHAKBULI1B6Vwlj03WjwMaCdyHWFECJ0+69OM4cxfXW+HMXhJly5ZNNch127ZtImXxtGnTKF++V9Wzc+TIIaw0PXv2pEWLFlHXrl2d3n8APD0+RWJ90RrU8sRWGlenG/1YSr5GyoxWpal5+RykVZANDADgiVgag8KWD2eKoXt2BO4rAYgVJeDiKvVsAWFXLn9/f6vex1nDTMW1sGsZLwCA1IWKb2ICJfi8uvQezV6EqvZcbDY+5eKkhh5tvgcAAK3CMShcoZ4LP+ocHKtijRiyNXBfKSB1sYZgt6/+/fvT+fPnDSo+AwBcI1TK3ThD2xf0oOJ3LyavSy2QHkLlFVzQtlGjRlS0aFG6fv26WMdFatniCwAASoSv3xyLwhhfyaXXtsSqGIshuXfz+mx2iCGlALGiIdhFi4NU2V2L66sAAFwnVBqe3U2rvxtOOaPu0ie7V6bahtYyfpmDr1uccZCzFl65ckW4pTKJiYk0depUd3cPAABk4VgUzsjFFhR9+LW9mbp8nCyGlALcwDSEVAwSAOAcol7EU+lxfxiu1Omox4E1NGzHUvFyS4G3qG+TwWbbgVAxhAUJpy5u06aNsKZIVKxYkUaPHu3WvgEAQGqwIOG4FA505/gRdstia4cjRET9/8SQcZ2V9EH+NK8J6qwAAAD4j8qTttCtp3EG63ySEmnclnnU9thm8XpJ+SY0vnZXSvL2kW0HQiUlFy5cEJkMjeECtZiEAQB4AixMnJWRq74JMZQnKIGyhb8u6O3JQKwAAIAT3L7SxsXQvF8mU83LhymJvGh8na605I3XGfiMCUvjTcciGzi5p55JRESEiLXLnTu3wfqdO3dS/vz53dYvAABQqhi6d+8eqQWIFQAAcEJqYqmw40vfAOrbdDBtKVhRto0jI+tRhmDrsvNpie7du4ukIJILGAfY79q1i4YMGUKRkZHu7h4AAAAnArECAABOqKGS6O1DvZt9Snme3KZTWeVn/+H2lTosSqKioqhevXoUExNDtWrVooCAABo0aBD17t3b3d0DAADgRJANDAAAHCRUal48SCO3LRBB9czzgEAIFQdw7do1Gj9+PD148IAOHDhA+/bto/v379O4cePENgAAAOoFlhUn8vjxY1ILnCrUz++VW4utqOl8AO0Sl5BEhUZuSrG+3dHfaOyW+eSjS6Lj2QrS+mI1zbYDoWI5efPmpdu3b1OWLFnojTfeSF7/8OFDsY1TGAMAgBJJTNI5JQuYloBYcQJp0qQhX19f2r59O6kFHgz4+MhnMLIUPi98fgDwRAavOUI/HrxtsM5Ll0RDdyylHgfWitc/lqhLmwpXMdsOhIp16HQ68vJKeXOPjo7G9QQAoFg2n7ydIqUwF2nk2idqSCnsKiBWnACn0+QCZuxbrRZ4BjNjRvtT7vHAgs8PAJ5G3qEb6ZVz12sC4mNp+sbp1OjcHvH6i2rtaE6l1kQmBtZM1fwZaEW3lCl4gWk4qJ5hoTJs2DAKDAw0mEBhl7AyZcq4sYcAACAvVHqtOJLivnEnKkast7cgpJaAWHESPCBX06A8KSmJMmXK5O5uAKCY+JQML6Jo4ZpxVO7WOYr18aUhDfrRuuK1ZNv48oOy1KR0hJN7qi5OnDiRbFk5c+YM+fu/zpjGz8uVKyeC7AEAQGmuX2xRMRYqDK/j6SzezrVR4BKWOhArAABgQ8avYncvUanbF+hJmmDq0XwE7c9VUraNi5Ma4oZkA5IrbadOnWjWrFkUEhLi7i4BAECqcIyKvuuXKcHC23k/ZxWKVBMQKwAAYENq4t15y9LARv3pRHhBupQxh+x+iE+xnyVLlri7CwAAYDEcTO/I/bQOxAoAAFgoVJqc/ov+yVaIrqV/5Wdszu2LgVBxHN9//z1t27ZNVGVmt1R91q9f77Z+AQCch6dm0uK+OnI/rQOxAgAAekTHJFCJyN8NV+p09L+/v6OBu1fSxQzZ6d0O0+lZQJDZdiBUHMfgwYNp5syZohhkRESEycxgAAB14cmZtFhUcV85mN5U3ApfwcJDX4kvkDoQKwAA8B/1Z2yns3dfGKzzS4ynSZvnUquTW8XrLQXeomj/tLJtcILvixAqDmX58uW0evVqatmypbu7AgBwAZ6eSYutPyyquK8sTPSPQ5pq4e2eYCVSAqhgDwAA/7l9GQuVkJhoWvrjGCFUEr28acTbH9OUWp1J52X60lmtUEYIFSfAbl9IUQyANkgtkxbD23k/JcNiikUVW1D04ddKF1tKA5YVAIDm6bBof4q5m+xR92jJj5FU6OE1YUnp0/RT2pH/dfV0Y+a0KUONy2R3QW+1R/fu3WnFihUUGRnp7q4AABSUSSt/OlI0LEg4PbEnxt0oCYgVAABp3aJSJDTl+jHbvhFC5U5wBurcMpJOZ80n2wZSEzuXJ0+e0KpVq2jLli1UqlQp8vPzM9g+e/Zst/UNAOC+TFr50xleC5QY+M/rLU1P7KkJBZyNr5rcBPiGNX/+fLp8+bIoYPjee+/RpEmTkoszcmGxZcuW0cSJE+nKlSuiInvHjh1pzJgxlDatvA86AEB7qYmHvdNHzOCNqduT7oTIF0RFIL3zOX36dLIb2NmzZw22IdgeAHVhXSatRFJK4P+oRkUpfVCAzULDkxMKOBvViBV2Dxg/fjx98sknVLNmTVHteNy4cfTo0SPhPiDl6u/SpQu1bt1aiJhjx47RlClT6Pbt20LEAAA0nPGLiMrePEtHsxcRzx8GhVGPFiPNtgOh4trikAAA9WNNJq2HD+4rIvCfBcbHq44arLNGaHh6QgFno5oA+y+//JLatGlDM2bMoGbNmtHQoUNp4MCBIoPMixcvhFWFxUvVqlXFulatWgkLC+/HmWauXr3q7kMAALiABjP+SiFUvHRJNGbr1/TzikHU8sSrrF+pAaECAADOy6TFGNsl3JlJy1zgvykkocFCRAsJBZyJaiwr3t7eFBYWZrAuJCREiJT4+HhhPWFBMmDAAAO3AcnKsmPHDvroo4/c0HMAgDvdvtLGxVC/PT9S/cMbxOuwl0/NtuHnRXRhMoSKs2natKmwivN1nJ+bA0UhAVAXUiYtY7eocDe6RaUW+G8MSwsebfIxcJC9nLiyJqFAJQtjX9SGasQKu3+xS9fbb79NNWrUED7Os2bNoq5du1JoaCidPHlS7FegQAGD9+XPn1883rx502S7sbGxYtF/nZCQ4NRjAQC4Rqhkjn5Ei9aMo/DwYIr18aP+jQfSb0Wqyraxb2gdCg9DxWFXwDGF0sQSPwcAaAsWJLWLZKVv916hq49eUO4MgdS+Uh7y9/VWdOC/tULDmoQCWkU1YmXQoEEiU0yLFi2S1xUsWFDEsTBsXWECAwMN3ie9jouLM9nu5MmTaezYsQbrGjZsSB9++CH5+/uTlrLxaA0tHrNaj5tTExtn/OLUxAOOraDM4cF0q0hxGle7G13KlJOKUJLJNpZ3eYso7indu2fe8uIJyF3vlATHGJp67iyOHDlCn376Ke3fv19c2zn2cfr06ZQrVy6nfzYAwLKA84W7L7vNsmJp4L+1QsO6hALaRBVihV29OE7l8OHDQpxUrFiRTp06JWJUqlWrRkePGgY9WZNJZtiwYcJ1TN+ywjfOzJkzU0BAAGmJLFmykNbQ4jGr7bhfWVQMZ+IyPn9CPyzoRyGxz+lS+gia1vAT2pWYiShKG/Ep+tZi8CrbGN8r3nrrLVq8eDE9f/5cTFI1btxYJGJhN2MAgOtQYsB5aoH/tgoNaxIKaBVViJW9e/cKqwqnLe7Ro4dYV7duXSpUqJCwgqxdu5by5s1rYGExnmE0ztsvwYJEX5TwTd7XVxWnDQBVw8GI+Yf/ZnIbZ/pa9EYzqnL1GHVvMZLCOb25RoQKSMmQIUMoT548tGHDBgoKChLr8uXLR+3bt6dz585R0aJF3d1FADRDagHn+nEg7gj8Z7HEfbBEsFgiNMy1686EAkpCFdNFXDOFqVy5ssH6KlWqiMcbN25Q1qxZxXPjrF+XLl0SjxERES7qLQDA2aw7ejOlUNHpKF3s8+SXs6p8QG3bTKQnaUNMthGWxgdCRQNERUXRpk2bqGfPnkKosKWe63axpYXvLRAqALgWawLO3RX4zwIkNawRGnLt8ut5Gk9bzKjCRMAzYszu3bupZMmSyet37tyZHLvCgfXZs2endevWiVorkvuXlEWmevXqbuk7AMCxvPPFNjr3wPBG558QT5//NpNyP7lFbT6YTDF+adgHlOJ9TFtUi2ULot/61XRRj4E9XLt2jXLkyGGzqxbX5GJxwrEpzZs3F8IlMTGRateuTXPmzBH3D1Mg+QoA2qxgz8KBrTr6leYfP4+j8Rvty1xmql1UsFeRWOEYlUaNGom6Krdu3aI33nhDVDnm7GBsbeF4FhYno0aNErNnvPD+HNfCfslcn0XKCgYAUFfGL05F/PXaifTWjVMU7+1Db9w4Q7vzlpVto27RzLTwowpO7ilw5GRV8eLF6auvvhLWEGu5d++eeOzVq5eYtFqzZo24j/D9on79+iKexVR8oqnkK/Xq1aP79++bTL6ixsQVjgLnxjxaOT/s/nXuzjO6eSuKioSaTnSiT9TDB7T3SRyF3npKhcPTuXxQnz8d/SeWEonS+dDaTiVE/6NexlFoWv/kPknXGFvatafo5RMbfjdKTb6iCrHCM2pc6PGzzz4Tj1988YVw6+rYsaO44UjxKN27dxf7sohZtGgRZciQgXr37p2cMQwAoC6hkvvxLVryYyTle3yLnvoHUq/mw2lPnjKybZwZV5/S+vs4uafAkXBA/OXLl0VGSM7kZS3Pnj0Tj2xB4fuHZHVn9y8WP99++61IgW9J8pXZs2ebTb6ipsQVjgbnRtvnJ2XmL/OWUtYlo7fcEKLmbNQ9q6rFO5Ns4a9CDjz1dxOr0OQrqhArTLp06WjChAlikYNvQt26dRMLAEDdQqXcjTO0YO14yvjyKd0IyUydW46h85lfuYyaAvEpnglPSjHGVg5LkYoJSxZ4/ZhHvq8cOnTIpFhB8hUAnJv5yxzGxdzdmSUMOB9VBNgDALSLKaFS6+JBWv3dcCFUjocXoObtp5kVKqKGClB0XAoHvhvD63ibPXC8i9SWqfZZsAAAXJ/5yxqk93Nb3CZQFxArAACPJOpFvEmhwlzImJOeBgTRlgIVqPUHU+h+sHzaSFhUlA+nnudYEGMePXqUnJbeVkqUKCHchn/55RcDwbJ161aKjo6m8uXL29U+AMD2zF/W4M4sYcC5wGYNAPA4qn22ja4/NrrB8UDzPzeeG2Hh9F67z+lGaBZK8jYdg5LWl+jMBAgVT4BFhKkCviwm0qSxr6qzj4+PcB/u3Lkzvfvuu9S2bVsRYM9FhcuUKUMtWrSwq30AgP2Zv9zdJnAvECsAAI+i4LCNFG9k5Q+OfUGz10+lFWUb0p8FXmXyupZe3m/5o4q5aOy7r9OcA2XSt29f8chChQPaAwMDk7dxeuEDBw4IQWEvnTp1EqJl6tSpohBkSEiISGPMSVtMZfYCADiGTMGmk1HYg7lq8cAzgVgBAHgMpty+wp8+oCU/RVLR+1eo1J0LVK3HInrpL3+zOj+hAfn7wgPWEzhx4kSyZYXroegLB35erlw5kQXMEXTo0EEsAADXBdZHrj/lsPYsqRYPPBOIFQCAxwqVYncv0eKfIik8+hHdDwqjLu+NNitUEJ/iWWzfvj3Z8jFr1ixh8QAAKAsOaLe2kKEtGcAcVS0eeB4QKwAAjxQqNS8epDnrp1Jw3Es6nzEXdW41hm6Eyue4h1DxXJYsWWIQp8IEBwe7sUcAANP1USjVmieWZgDjdpqWzkbf7LwsXpvb39pq8cCzgC8EAECx8E3NlFBpd/Q3WrRmvBAqe3KXopbtpkKoqJyZM2dSrly5KDQ0VCw5c+akGTNmmEw5DABwPpJ1xDibl1TzhLfbkwHsi5alaVjDYqJ2CosRfcJDAqh/3YLUq2Z+Wt2tIu3+tDaEioqBZQUAoEh+PnKT+v9wzOS2Qvevko8uiX4sUZeG1+9N8T5+su1AqHg+Q4YMoW+++YYGDx5MlSpVEuv27t0rMnbdvn1bBMYDAFyHOesIr2NHLN5er1h4CrcsS7N1PXj+qpo6ixBux5Sr2b179yhLlowOOSagXCBWAACKo+LEP+jOs3jZ7WPrdqdDOYrS+qI1ktMVG4PUxOph4cKFYmnZsmXyutq1a1PhwoWpR48eECsAuJjUrCNSzZOley5Txyp5DQSLpdm69Pfj91fKD1GiVeAGBgBQFOz2ZSxUMryIoqHbF5NvYoJ4nejtQ+uL1ZQVKp+3LAWhojJKlSplcl1SUpJb+gOAlrHUOjJ+4xmq+tmfBi5hbBXheBS5MHhez9uR1QtIQKwAABSDqfiUfA9v0M/fDqSeB9bS0B2vA63luDipIbV6I6eTegjcAacUnjt3bor18+bNE3VRAACuxZpaJsYxLGwl4WB4xliwIKsXMAXcwAAAihUqFa6fpG/WTqCwmGi6GhZOq8o0MNsG4lPUSWxsLK1atYp+//13qlixoli3f/9+UWmeK85LxSOZ2bNnu7GnAGgDyTrCQiS1FBemYlg4DoUD540ziSGrFzAFxAoAQJFCpdmp7TR10ywKSEygIxGFqet7o+lRYKhsGxAq6uXs2bOiACRz9epV8RgeHi4WLhYpwZXuAQCuoc2buWjG1vMW7SvFsHCsixR7Yi5wHgB9IFYAAG7jZVwiFR292XClTkf/+/s7Grh7pXj5W6HK1L/xQIr1C5D1Zb0EoaKJ4pAAAGXWVrE11gWB88ASIFYAAG7hwwV76O+LT1Ksz/70PvU4sFY8n1+hBX1WsyPpvEyH1+UIS0O7h9Zxel+B+7l7966IWzl9+rSwoBQvXpw+/vhjypIli7u7BoBmsLfyvDWxLu5KyQxLj/KAWAEAKMLtS+JmaBbq3Wwo5Yy6SyvKNpTdr1OV3DSmSQkn9RAoiT179lD9+vUpa9asyXVWVqxYQdOnTxdxLNI6AIDzsLTyvCm8/otHUXKGL1MWI47LQQyN+4FYAQC4Xahkj7pHmZ4/pn8iCovXf+Urb7aN8xMakL8vkhlqhUGDBtEHH3xA8+fPJ2/vV987pyzu2bMnDRw4kP7++293dxEA1WNp5XnywAxfchYjKZMZJwNwhmCBJccyIFYAAG4VKiVvX6DFa8aST1ISNW//BV1NH2G2DQTSa49jx47R0qVLk4UKw88HDBhAZcuWdWvfANAKltZWCUvrR09exntMhi9zFiNTmcwcBSw5lgOxAgBwm1Cpd2EfzV7/OaVNiKUzmfNQrI+/2TYgVLRJaGgoXb58WVSs14fXhYWFua1fAGgJS+NN5n5Yjry9vTzGWpCaxchUJjN7LSiPn8dR71Wut+R4KhArAACn8ig6jspN2JJifadD62jUtoXkTTrakbc89Wn2KUUHBMq2A6GiXdq0aUNdunShqVOnUuXKlZPjWD799FPhHgYAcH9tFSkupWL+jIoWJ7ZajCzdzxILCp8eV1tyPBmIFQCA0yg39nd69DLBYJ13UiKN+nMhdTq8QbxeWaY+ja7XixK9fWTbgVDRNixSdDodde7cmRISXv2e/Pz8qFevXjRlyhR3dw8ATSBVnueZfx5C6zwsLsVei5EtmczkYmGSzGQpcIQlR204JEJ106ZN1LhxYypWrBhdv35drFu4cCFt27bNEc0DADzU7ctYqDCdD61LFiqTanaiEW/3lhUqVfKFQKgA8vf3p1mzZtHjx49F/Aovjx49ohkzZohtAADXIFWeZwuKPvzaU12XJIuRnMTi9dlsyGRmT/Y0eyw5asRuy8rKlStFRpauXbsKcRIf/yqoKjExUcyG1amDGggAaA1zqYm/LduIal88JNIS/1akqux+Z8bVp7T+8tYWoD0CAwOpZMmS7u4GAJpGbZXnnWUxsjV7mjlLTqJGs4fZLVZYkCxYsED4FLM1RaJixYo0evRoe5sHAKgkNfGtkEyiuCNXov+wzUQiL/kLLKwpAACgXNRWeV6yGBnHltiTycxWy4hcTZrNGs4eZrdYuXDhgsmCXMHBwfT06VN7mwcAeBCmhEqVK8do3s+TaFn5JjStevtXKyFUAAAAqNhiZEuMi5wlZ7Ob6sCoJmYlIiKCzp8/n2L9zp07KX/+/PY2DwDwANg0bUqotDq+hZb+OIZC4l5QhRunyC/xde59U0CoAGu5du2aKBAJAACOshg1K5NdPNrjYpVaLAxj3Lyp2J/EVOrAMLyd91MrdltWunfvTn379k12AeMA+127dtGQIUMoMjLSEX0EACiYn4/coP4//GO4UqejAbtWUN+934uXvxSrQUMafELxPn6y7UCoAFvIkycPFS9enL766iuqVq2au7sDAFAwroz5sCQWZs4HZSl9UIDZ/hxwQR0Y1YsVFiVRUVFUr149iomJoVq1alFAQAANGjSIevfu7ZheAgAUSYUJm+ledKLBOv+EeJq6aSa9e/ov8Xp2pdY0vVo7WdevvBn8afuQei7pL1AfixcvFsUh+Z6zf/9+d3cHAKBQ3BHz4YhYmHtOrgOjmTorEydOpBEjRtDp06eFOZ5TGHPMCgBAYxm/dDpa9NNYqnb1GMV7+9Dwd3rTj6Xelm3jZOQ7FJwG5Z6A7XTs2FE8jh071t1dAQAoFHfGfNgbC5PFiXVgPAVfR6aUfOONNxzVHADAE1MTe3nRTyXrUOk7F6jXu8NoT54ysm3A7QtYE5eSM2dO8jKyznGhSHY9zpUrl9v6BgBQNqnFfLiiYrw92dMq/Bf7wsJKZ0X2MDVht1hp2rSp2e3r16+39yMAAAoXKj5JicmFHdcVr0V/5StPT9KGyLYBoQKsIW/evHT79m3KkiWLwXouDMnbuK4XAACYwtNjPnycVAdGU9nAMmbMaLCEhIQI/2HOBpYpUybH9BIAoFih0vDsbvpjUW/KHP04eR2ECnAkbEExtqow0dHRlCaNel0fALDVkrD34kNad+ymeFRzlihLUEPMR/3/Yl/YgpJa9jA1YrdlZcmSJSbXDxw4UAgXAIDn8zIukYqO3my4UqejHgfW0LAdS8XLzofW0Wc1X8UPyAGhAqyBM00yLFSGDRsm3I0l2Jpy4MABKlNG3tUQAK2h9sKBxtm88gTpNBPzUd/BdWA8CadFtvbo0YOqVq1KY8aMcdZHAABcQIdFe2nnhUcp3L7GbZlHbY+9EjBLyjehz6WCjzIm3EsQKsBKTpw4kWxZOXPmDPn7+ydv4+flypUTWcAAAOovHGhKiFWO8KUOtUuZPS41xXz42BH74sk4TaycO3fOWU0DAFxEvqEbybjcXlDsC5q77jOqefkwJZEXja/TlZa80Uy2jeLZ0tHGftWd3legPrZv3y4eO3XqRLNnz6Z06dK5u0sAKBIlBJG7Q4g9fh6XqhBDzIfn4+soM70Ez4BxIOSmTZuoc+fO9jYPAFBQfEqWZw9p6U+RVOzeZXrpG0B9mw6mLQUryraB1MTAEbC7Md9T5s6dS5cuXaLff/9dZAfjYsQcYF+nTh13dxEAt+LpQeTmsLSCuzkh5oh6J8B92D2KOH78uEHgo7e3N2XOnJlmzJgBsQKAylITx/gFkF9iAt0PCqMu742m49kKybaB+BTgKFauXEk9e/akrl270rZt2yg+Pj45bmXq1KkQK0DzqCGI3NlCTMsxH6R1sbJjxw7H9AQAoAg6LNovmyjwaZpg6tQqUjy/EZpVtg0IFeBIWJAsWLCA2rRpI6wpEhUrVqTRo0e7tW8AKAG1BJE7W4hpNeZD06mLeXbrrbfeQnwKACq2qHx4bBN1OLwh+TWLFAgV4EouXLhAlSpVSrE+ODiYnj596pY+AaCklMNSELmcjYDXZ/OQIHJjMgUHOHQ/oDGx4ufnJ2qqmMp/7w6OHDlC9erVEymTucZLy5YtReVj/XiapUuXUsGCBUXfw8PDaejQofTy5Uu39hsAd3Pz0csUQsVLl0RDty+mSb/Ppcit31DxO/+m2g6ECnAGERERdP78+RTruZ5X/vz53dInABwdQF71sz/pgwX7qN93x8Qjv+b1liAFkTPGIzKPDyK3VLNpu5yMqrG7KORHH30kzPPu5vTp01StWjXhw7x48WKaNm2aEC+NGzempKSk5CBNzipTvnx5WrVqFXXp0oU+//xz4QsNgFbJP3QjVZn6p8G6gPhYmrPuM+p5YK14Pb1aWzqVVX5QGOANoQKcR/fu3UUylz179ojX169fp2XLltGQIUOoV69e7u4eUBmuLqooZboyjsuQUg4fuvJQ04UDHzyPdeh+QIMxK8+fPxfBj1u2bBEiICgoyGA7p5t0BXzTypMnD23YsCG5D/ny5aP27dsLN7UiRYrQuHHjRO2X1atXC2tQq1atxH6TJk0S23Lnzu2SvgKgZLevdLEvaPV3Y6jcrXMU6+NLQxr0o3XFa8m2cXB4XcocAvM7cO71PSoqSljOY2JiqFatWhQQECBqrPTu3dvd3QMqwtVFFS1JObxi3zV6540iFllF1BhEruZ4HOAiscKFurgwF8MpJfVxlXsY38Q4reXMmTOFUGF3L17Y0nLlyhWxz8WLF+nq1as0YMAAg361bt1aiBVOFMBWIgC0LFTyPbxBo8/+QCVvnaMnaYKpR/MRtD9XSdk2YE0BrmLixIk0YsQIYUVna3mxYsVEzAoAnlxU0ZJMV4+ex1mVclhtQeRqKuoI3CRWpKJd7oQFE9+8cuXKRc2bNxfChd3BateuTXPmzBExKnfu3BH7FihQwOC9kr/zzZs3TbYdGxsrFv3XCQkJTj0eANyVmrjqlaOUNfoRXQ0Lp04tI+lSxhyybUCoAFcTGBhIb7zxhru7AVSIu4oqqjnlsKNAUUdgt1jhAHYuzmXKisLbWEA4m3v37olH9l2uXr06rVmzhm7dukWjRo2i+vXri5k4KS8/3+z0kV7HxcWZbHvy5Mk0duxYg3UNGzakDz/8kPz9/UkrPHnyhLSGWo+ZUxMXCTW97UDtRrTmSV7aniY3+acJoiIp6te/YnmXt5L/79SAWr9rOeSud0rm+++/FzVW+HcnxSFKrF+/3m39AurAXUUV4eJkGXJFHdMH+dO8JqU8Nh4HuEiscPVgrlifJUsWg/UPHz4U29jC4WyePXsmHtmCIsWjMEWLFhWuYN9++20Ki4qlrmrDhg0TrmP6lhUO1OfCl+wzrSWMv2MtoKZjjo5JoBKRvxvm1dDpqN2xTbSuWA16FvBfvFlYMTob5U0Uqy2Lipq+69TQtxZ7AoMHDxZuvhyrwpnBlJKBEqgHSy0Xd57GuNzFKUOQP1ycZOJx8gQlULZw+VT6QB3YLVY4NsTUjSM6OprSpHHNTEBYWJh4bNasmUFfqlSpQunSpaNDhw6JAHtGsrAYzzByKmNTsCDRFyV8k/f1tfu0AeBS3p7+J52/Z5ii2y8xniZtnkutTm6ld879TR+9P5aSvH3MtqNWoQKUzfLly8VEFKejB8AZWGq5GP/rKUrr5+2wmXxLXJzaVcwFFyeZeBw1WfiBPDaPujmNJMPigK0P+u5VbE05cOAAlSlThlxBjhw5koWTMbyOBUvWrK+UNwfZ6yMlBeDZOgC0Ep8SEhNN836ZRFWuHqcEL2/aVKSKWaFSJEta2jygtpN7CoBp2O3LVfcToE1Ss3BIPHoeb3OwPcfFmMrSJefixEHjLGTKZTE/iQSA2rFZrJw4cSJZDHCAu378Bj/nDGGcVtIVlChRQoiNX375xSDb19atW4WFh1MqsxtY9uzZad26daK+irSP5OvMsS4AaEGoZI+6R0t+jKRCD69RtH9a6tP0U9qRXz5o+WTkOxScBtZE4N46KytWrKDIyEh3dwWoFHMWDlNYG2yfWkpkcymHYT0AWsfX3ixgXGRx1qxZomq8u/Dx8aEJEyZQ586d6d1336W2bduKAHuuncKzcS1atBDihAPuuQAkL40aNaJTp06J4Pk2bdqgCjLQhFApefsCLV4zljI/f0J3gjNQ55aRdDprPtk24PYFlJIAgQv5cj2vUqVKpXDbdVU9L6BuJAvH8J9PinTBjgq2tzQlstpSDgPgKOyeLuVgcyXAoolFy9SpU0UhSBZPnMb4s88+S7b68Oyct7c3TZkyhRYtWkQZMmQQBcXGjx/v7u4D4HSh4p2USDN/nSaEypnMeURq4jshmWTbgFABSoEzOkpuYGfPnjXYhmB74EjYunH8xhP6aodh3Thbg/LdlRIZADWhKt+ODh06iEUOvql169ZNLACokZdxiVR09GaT2zgm5eN3h9L/9nxHQxv0pegAwzTe+kCoACWhhHpeQP2YctWyNyjfXSmRAVATqhIrAGiZjxb+TX/9+ziFNaX43Ut0IltB8fpc5jzU592hZtvhGioAuJvr168LazmSnwBXIOeqRXZWTFdq0Ue5YH8ANCdWuCgkZ+pi1ysAgPPIN3RjivKNaeNiaPaGz6n65cPUrvUEOpizRKoXg3+nNEIwJ3ArEydOpOnTpycX6gwKCqLatWsLN14uyAuAozHnqmWMtRXTlVj0MbVgfwCUhlNVRJ48eah06dK0a9cuZ34MAKT1+BRjoZI5+jF9v3oo1ft3v7i9Znj51GwbB4fXFUIFAHcLFU6E8vjxYypcuLDIKhkcHCyyNjZp0oQaN25MUVFR7u4mUBmpuWrpwxYVa9IWSymR5WQNr+ft1hZ9ZIG19+JDWnfspnjk19ZYkIyPVwr25+0AaMqysnjxYrp8+bJIYbx/Pw+aAADODqQveP8qLfkpknI8vU8P04ZQt/dG0ZHsRWXbQHwKUAoLFiyg9OnT07Zt2wzqqhw8eJCmTZtGP/zwA73zzju0c+dOg3T5ANiDpS5YfWoVoP71ClnlLmVJ0UdzVhoWIadvPaW9t+KT3bW2nL5jk2UEwf5As5YVdvWSK8bIpntODQyhAoBrhEqVK8dozYrBQqhczJCdmrefBqECPIabN2+KVPLGBSDffPNN+u677+jLL78UBYfZTcyZjBkzBlnGNISlLlhVCmSyaRAvpURmq4w1Vhq2clT97E+asukM9fvuGH2wYB+Vn7CFetpoGbEm2B8AVVlW8ubNS7dv36YsWbIYrH/06JHYxtXsAQCOIy4hiQqN3JRifanb52npj2PILymR9ucoTj1ajKAnaeXrH0GoAKURGBhI6dKlk93OqeZZtCxfvpyGDjWfKMJWWAyxOxrQDqlVr7cmoF4Oc0UfUwv4Dw19vf7Ji3iT+1tiGVFqsD8ATressAXF1AwUV45Pk8Z1AWMAaIExv5wwKVSYE+EF6PdClemXYjWofesJECrA4yhbtiz98ccfZvepVq2acC92Bi9evBDp79nVDGgHyVWL8bIzoD61z+H0xM3KZBeP5ly/LA34t8YyosRgfwCcalnp27eveGShMmzYMDEjJsHWFJ6dMjblAwBsp9jI3+hFguHtyz8hnrxIR7G+/qTz8qb+jQdQvLcv/2OabMPfm+j8JAgVoEwmTJhANWrUoOHDh9OkSZNM7nP//n3KlEm+mKk9sLWG71v169en3377zSmfAZSJ5KplHAsS7oYsWdYE/FtjGXGFBQkpkYGixMqJEyeSLStnzpwxCHbk55zFhQPrAQDOiU8Je/mUvl47ke4FZ6C+TQcLsRLv4yfbxmctSlDrCrmd3FMAbKdq1apCMEyePFlkkRwwYAA1aNAg2Uq/adMmWrVqFX3yyScO/+ytW7eKAP6TJ0/Sr7/+mur+sbGxYtF/nZCQ4PB+AddhrauWs7DXDevC3WciQ5hx3+0N9k8NpEQGihMrUkXhTp060axZsygkRN7lBADgWKGS6/FtkfEr/6Ob9NQ/kPI8vk2XM2SXbePipIaY3QIeAceLsOVk9OjR1LJlS1Gni1/HxcWJ2iucwpi3ORJul+9l8+bNs9hqw4KKE8joU69ePWH5MZWpTKobA0jx5yZ/Ol544ieRHj647/LPD6WXVCT0dUL6nEEsK4wT1Muz9cgFsWQI8qd2FXPRG3kyJm8rl8WH5rbIRyv2XaNHz+OS10v78nZbam0duvKQvtz2L4Uaxdh40QuateEQJb0oYNAPtf52lMQTG84NX2dVGWC/ZMkSgzgVhvPiAwCcI1TK3TxDC9aMp4wvn9KNkMzUueUYs0IF8SnA0+jfv7/ICsZWlHXr1onUxZIVY8OGDRQWFkYlSpQQFnyOc+HHChUq2Px5ffr0Ee5nzZs3t/g97P7Mlh8J7t/s2bMpc+bMFBAQYPI9xolowGtwbl6TMVNmGrbpqp67VhKdjbI+xNgrKoH2rr1E89qFGVg2GmbJQu+8UcRhFiR2/Rq/6ATdlukjtzp+2y3a/WkRp0yaWfvb0ZKrWhYrz42+tVh1dVZmzpwpUkly2kkmIiJCXMTZVI/0jwA4Tqg0PLubZvw6jQIS4+l4eAHq8t5ouh8s718MoQI8lWzZstHAgQPFwnGQp06dosOHDycvx48fF48M32dszTy5du1aEdR/7Ngxiol55b4SH/8q4xK/9vHxIT+/lO6VLEj0RQnf5H19nVq6TLPYO7j0tMGpvruWPZjLECYF+zsCa1IiO+ozbcXTXdUSPey37CjsvrIOGTKEvvnmGxo8eDBVqlRJrNu7dy+NGzdOpDSeOnWqI/oJgGaIehFPpcelzIjU4fAGGrf1a/F8S4EK1LfJEHrpL5+1BUIFqAUWDKVKlRILu2sxSUlJdPr0aTp06BAdOWL7oI4LTLLrVvbsKa2TadOmpX79+okJOeCZg0tPHZxKAf9D13B8sO2z3a4QCp6SElk/HbSpGjXmat4ogc0e+ltWhFhZuHChWNi3WIKLQRYuXJh69OgBsQKAFVSetIVuPTXtM3omS16K9fGlVWUa0PjaXSnJ20e2HQgVoHY4loVdwXjp2LGjze2wBwC7nOmzceNGkZmMJ97Cw8Md0FvgjsGl0genqc2Sc9/SpfGj8T/utfuznCkUPCElsrl00JbUqHG3FWSzwn/LzsYhNmue7TK1jme+AAC2u32RTpechvhgzhJUv/Ncs/EpEen86O8RbzuzmwCoijx58ohFn7Nnz4rHihUruqlXwN7BpTsGp86YJa+YL6MIfuf4E2vrrrhKKLgiJbKnuao50gqSqPDfskcUheQCWnPnzk2xnrOqtG/f3t7mAdCsUAl/+oC+XzWUCt+/krzOnFD5Z/TbECoAAFVgzeDSGe93JtIsuXH/pFly3i7Bg0/O0iVXsJKXsEC/FNv098nmZKHgqqKanuKqZs33awkHFPxb9hixwkGFS5cupSJFighTPC9FixalxYsXi5zzXDxSWgAAlgmVYncv0S/fDqC3bpyizzbNemVhMQO7fYUGytdYAQBYDt/HuIYYII8dXCo1jiK1WXKGt/N+Epzyl9182DqhD7/m9VNalHS7UJBibLKGGGbC49e2uCjx8XOtmHXHbopH/fOhZFc1W77f1Lin0N+yR7mBsbmc00YyV69eFY/s48sLF4uUQFYwACwTKjUvHqI56z+j4LiXdCFjTurTbKhsRXoG8SkAALVh7+BSqXEUtrojpVawkgWBsdtRuFuCr+Ukk3sDyV3lquYMd7MsCv0te5RYkYpDAgAsJy4hiQqN3JRifdujv9HYLfPJV5dEe3KXol7vDqenaeTrFkGoAADUiL2DS6XGUdgzS24u3XBqYsbZyAWA331qXQC4swLJLUkH7QgLlDOsIBUU+lv2KDcwAIB1DPv5aAqh4qVLoqHbF9PEP74SQuXHEnWpY6uxECpA03BMJBdfBNpDGlzKDc5SG1wqNY7CmbPkkphpVia7eHTVsTnK9ckZLlT6sMjpXj0vGZ8Wfs3rHWGBcsb366PQ37LHiZW7d+/S6NGjRfriVq1a0ZgxY8Q6AIAhBYZtpNX7b6VY75OURCXuXhTPv6jWjgY37EfxPqZjUPhyBKECtMCKFStEBXugXTh43BiOz7Nkhl2Ko5CL9XBHqldpltydAfGOxlEB4M4OJGerzTc7L5Ox1uHwNF5vbeC7K7/f+gr8LXuUG9iePXuofv36lDVr1uSikCtXrqQZM2bQ77//nrwOAK1jMjXxfyT4+NLH7w6jitdO0B+F5P9nprxbnNpUNEyzCgAAakPOHUgqnGsp7naPMueOxD3QqWCW3FGuT84MJHdV+l9nfr/1FfZb9iixMmjQIPrggw9o/vz5okgXw/VVevbsSQMHDqS///7bEf0EQHVCJd/DG9Tg3B6aW+l9EUDPLl/mhMrFSQ01cVECAGgbcwNLCWsGluZiPdyBNEuujIB4UozrkzNd5FxZZ8WZ36+Pwn7LHiNWjh07JlIXS0KF4ecDBgygsmXL2ts8AKoUKhWun6Rv1k6gsJhoehgYSt+VqW+2Dbh9AQC0gqMHlo6sJO4o1DRL7qgAcGcGkrs6/a+avl9ViJXQ0FC6fPkyFS5c2GA9rwsLC7O3eQBUJ1SandpOUzfNooDEBDqarTBtKWi+SjaECgBASzhyYOmMNLiOQi2z5I5yfXKmC5U70v+q5ftVRYB9mzZtqEuXLiJOhQUKLxwU2bVrV+EeBoAWufMkJqVQ0emoz9/f0axfpwmh8luhytTmg0n0MEhe1EOoAAC0hqMGlo6uJA6cHwDurEByT0ts4OiimKR1y8rUqVNFpd/OnTuLivWMn58f9erVi6ZMmeKIPgLgURQavpHikgzX+SXG06TNc6nVya3i9dcVWtCUmh1J5yU/XwChAgDQIo5wB3JVQLWz0HddC6WXlDFTZkX20xmuT85wofKkxAZKtgZ6rFjx9/enWbNm0eTJk+nixVepV/Pnz0+BgYGO6B8Aqsj4VfbWOXrv5DZK9PKmMfV60oqyDWXbaFshJ01sUcqJvQQAAOXiiIGlKwOqnT1YLRKaRMM2XfWIwaqjXJ+c4ULlCYkNnFUUk7QuViRYnJQsWdJRzQGgqtTEB3KWoNH1etKN0Cy0I/+bsvudn9CA/H1RqxUAoG3sHVi6OqDaUWCw6lyUHPhuaVHMegq1BnqEWAFAy5gSKiVvX6DHgSF0IzSreL2inHm3Lrh9AQCAYwaW7giothdPd13zFJQa+J6aNZDh7XP+vED96hYiLYEpXACcIFTqXdhH368eSkt+jKSQmOhU24BQAQAA+YFlszLZxaOlg3RPC6h2RQV3oGwstfLN2HpBc8khnCpWrl27JgpEAqBGomMSTAqVTofW0ddrJ1JgfCzdCslMSWaC6BkIFQBMkytXLsqePbu7uwE8OO6FMRYsSguo9nTXNeAYrLHyjd1wWlMZwpwqVvLkyUOlS5emXbt2OfNjAHA5b8/YTiUifzdY552USGO2fk1jti0gb9LRqtL1qct7oyk6QD7ZBIQKAPJcuXKFtmzZ4u5uAA/FWWlwnYUnuq4BxyFZAy3htsYsbE6NWVm8eLGouzJo0CDav3+/Mz8KAJdhypqSNi6GZm/4nOr9++p3PrlmR/q6wntEXqZn7YplC6Lf+tV0el8BAEDLKDmg2pUV3IHnWAN7rjhi0f73NGRh83WEq1fOnDnJy2hQxrVXateuLcz4Y8eOtfdjAFB0xq/IrV8LoRLr40cDGg2gjUWrybZxMvIdCk6D3BYAAKDlgGpPrgUCnCeu+9ctRDO2nk913ywasrDZ7QaWN29eun//for1jx49EtsA0EJq4i+qt6cTWfPTB20mmRUq7PYFoQIAAEANrmvA8fSpXYDCQwJkt3spMDmEs7F71MQWFGOrChMdHU1p0mhH9QHtCZXsUffoZmgW8fx+cAZq8tFMWbcvBvEpAAAArHVd4wr21Urlh0VFI/D3HNm0uLCwMTq9bVq1sNksVvr27SseWagMGzbMoGJ9YmIiHThwgMqUKeOYXgLgJjjbRv7hv6VY3+r4Fpr4+1wa3LAfrSte69VKCBUAAAD/3TvsiZPRd127d++epgamwP6iqGrDZrFy4sSJZMvKmTNnyN/fP3kbPy9XrpwIrAfAU9l66hYN//2g4UqdjgbsWkF9934vXla7cuy1WJEBQgUAALQD18AwHmRm0+ggE2gjOYRixcr27dvFY6dOnWjWrFkUEhLiyH4B4Ha3ryKhSQZhXf4J8TR100x69/Rf4vWXlVrT9GptZdsIDfCmf8Y2cEl/AVA7mzZtorlz59KlS5fo999/F4ldFi5cKGIj69Sp4+7uAZAsVNh9xzibF2f44vWIOwFqTA6h+AD7JUuWKFKojBkzxmSGsqVLl1LBggXJz8+PwsPDaejQofTy5Uu39RN4RnxK2Mun9O33I4VQiff2ocEN+tK06u1JJ1PwcVrLUhAqADiIlStX0vvvvy+u3ZwOPz4+PtnleOrUqe7uHgDJrl9sUTGVdlin0WJ+ACimKCTPeDVq1IiKFi1K169fF+t4xmvbtm3kDjheZuLEiSaFFVuCypcvT6tWraIuXbrQ559/Tj179nRLP4FnCJXAuJe0ZsVgeuvGKXrqH0gdW42lH0u9LdvGxUkN6b03cjq5pwBoBxYkCxYsoBkzZpCv72uHgIoVK9KxY8fc2jctwoPtvRcf0rpjN8Wjswffrv48W2F3HX3XL2N0GizmB4AisoHxjBcP9rt27Up//vlnihkvV5vnX7x4QR06dKB33nmHfvvtNwOryrhx46hq1aq0evVqYXVp1aqV2DZp0iSxLXfu3C7tK/CM1MQv/NPS5kKV6d1TO6hTqzF0PnMe2TYQnwKA47lw4QJVqlQpxfrg4GB6+vSpW/qkVVwdj+FJ8R+WFunTUjE/ABRhWVHajBe7dXEWMkmISLCf89WrV8V6ffew1q1bi8cdO3a4vK9A2ULFW8cxK6/rqDTqNAtCBQA3EBERQefPpyyStnPnTsqfP79b+qTleAxj64EUj8HbPfnz7MXSIn1aKuYHgCLEipJmvLZu3Uo//PADzZkzJ8W2O3fuiMcCBQoYrJdudDdv3jTZZmxsrDgO/SUhIcEp/QeuJzomIaVQ0emo+/41NGz7EgqIj321ysubnqSVj82CUAHAeXTv3l2ky9+zZ494ze7Gy5YtoyFDhlCvXr3c3T1N4Op4DE+M/+BMTWz1kcvVpMVifgAowg1MmvEydqFy9YzXkydPRDzKvHnzKFOmTCm2S+5p+vVg9F/HxcWZbHfy5Mk0duxYg3UNGzakDz/80CBds9rh86s2Rvx8nK4/eklFQg2tKR2O/Eq1nxylJ/nzU/voU7Q7j/l6Qcu7vCXy4KsFNX7XqaG1Y5a73ikVFiVRUVFUr149iomJoVq1alFAQIBIj9+7d293d08TWBOPYUv2IuO6JEk6nVM/z1mZm9g9ja0+LExQzA8AhYgVacaLA+qlGa9du3aJm0tkZCS5ij59+lCNGjWoefPmVr3POGOYMVzwcsCAAQaWFg7Uz5w5s7hZaoksWV5Va1cDr60pr42LQbEvaO66z6jm5cOURF60omwDWpi+HFGU6TaC/b3p5Dh1ZvxS03dtKVo6Zr6OeRqcNGXEiBF0+vRpSkpKomLFigkLPvD8eIxDVx7S+EUnDMRJWFo/p32eM0ExP+AJRUc1J1aUMOO1du1a+uOPP0SMDPdB35LCr318fJLjaaT1xjOMnMrYFHws+qKEb/L6sTlAHfEp4U8f0OI1Y6nYvcv00jeA+jUZRNcKviUrVA4Or0uZQ7QlVgFwN2wJf+ONN9zdDY8ayGQKDhBT/A+ex9o1qHFWPAbHnXy57V+6HWXolf7kZbxTPs8VoJgfcCabPSjphKPwVcOMF7uc3b9/n7Jnz55iW9q0aalfv37JwomD7I0D7yV3NqBNoVLk3mVa+uMYCo9+RPcDw6hLy9F0PFshKkKvA+z1QXwKAK6ladOmZrevX7/eZX3xxIGMPrYOaqR4DA5uNxUl4vWf9cCaeAwpLkXPE9dibPk8V4JifsAZbNZo0VGH1FnRn/GqUKGCy03zn3zyCe3du9dgGTlypNjGz3k7B9azmFm3bp1IY2x8k6tevbpL+wyUk5o4xtef/BMT6ELGnNS8wzQhVOSAUAHA9WTMmNFg4ULEXBySJ6pMxShqFbnsWY7IpCXFYzDG9gFb4zFSi4ORA/EfQIskemDSCUVZVr7//ntRAJKDjNmy4uoZrzx58ohFn7NnzyanUJYYNWqUqAnDCxexPHXqlAieb9OmDdJfqhj+x80//HXNHWOuZMhO7VuPp+uhWelpGnmhDaECgHvgOEFTDBw4UAgXYH4gow9v5+E978uuStYM9h0dj2FpvAnHr+i7hSH+w3FoLfbBkzng5CQXqhYrgwcPppkzZ4pYFXalSi1g3Z1wMgBvb2+aMmUKLVq0iDJkyCDcw8aPH+/urgEnseLvSzRy/RmDdV66JBq8cznty1mSduYrL9adymperEKoAKA8evToIQr9jhkzhrSONVYKewY1jozHsDTeZO6H5cjb2wsDagejxdgHT+aehouO2i1Wli9fLirCt2zZkpREx44dxaIPC6lu3bqJBWjT7YvrpkzbOIMan9tN7Y9spBo9FtKjQHmPab4dXoZQAUCRnDt3zt1dUAy2DFBsHdQ4Kh7jdV2SF2bjUirmzwhx4mC0GvvgyWTRcNFRu8UKu31xxXgAlC5UMryIogVrxlP5W2cpztuXRr39sVmhUrdoVtr8bgUn9xQAkBqcHl8fjju8ffs2bdq0iTp37uy2fikJWwYo7h7USHEwszYcQl0SBcU+2OomCJxLBSckudBMgD27Vq1YscIxvQHASUIl38Mb9PO3A4VQeZImWMSo/FK8lmwb5yc0oA6VDeOgAADu4cSJEwYLZ57kFPIzZswQC0i9erpSK6nz7P3/6hQQgyx9+DVm990f+wCUg48TklxoxrLClZ9XrVpFW7ZsoVKlSqWoVzJ79mx7PwIAu4RKhesn6Zu1EygsJpquhWalTq0i6WLGnLJtID4FAOXAtbFevHghXI4LFy7s7u4oFnPV05U+qHkjT0ba/WkRBHq7CC3HPng69TVadNRuscIzXJIbmJSBS0LJwfZAO6mJm57+SwiVIxGFqVuLUfQwKEy2DQgVAJQFT4BxmmLcT2wfyOij1EEN6pK4Di3HPqiB+hosOmqTWLl+/bqoCs/Zv7Zv3+74XgFgBfefxtKbk7bKbo+s24NuhWSmRW80o1g/+arzECoAKJOPPvqIFixYQJ9//rm7u+JxAxlHVbBXK1pM3avl2Ae14KMxce9rbaX66dOnC9cvJigoiGrXri3iVho2bOisPgIgS5HhGynGqNC8b2ICtT22iVaUbUiJ3j6U4ONLX1V6X7YNZPwCQNk8f/6cVq5cKdyNy5cvL+49+sDdWNsDGVvRaupecy6DSnQTBMDXGqHCRRWZIkWKiJvFzZs3RdHHDRs2UIMGDcTNJDRUPrsSAM52+wqJiaavfplMVa/+Q3ke36KxdXuYbaNthVw0sUVJJ/YSAGAvZ86coXLlyonnly5dMtgG9zBgC1pP3avV2AegcrHCJvj06dOLSvX6qYoPHjxI06ZNox9++IHeeecd2rlzJ/n7+zurvwDICpXsUfdoyY+RVOjhNXrul4b+yvtqcGMu45e/r90J8QAATgbuxsCRIHWvdmMfgGdi8UiNrSht2rRJUVPlzTffpO+++46+/PJLOnDggHATA8DVQqXk7Qv0y7cDhFC5E5yBWrWdSjvyv2k2PgVCBQDP4Nq1a6K2itw2AKwBqXtTugw2K5NdPEKoACVi8WgtMDCQ0qVLJ7u9d+/eVKVKFZFeEgBXCpV6F/bR96uHUubnT+hM5jz0bvvpdDprPtk2EEgPgGeRN29eun//for1Dx8+FNsAsAak7gVApWKlbNmy9Mcff5jdp1q1aiLFJACOJi4hyaRQCXv5lKb/Oo0C42NpR97ywqJyJySTbDsQKgB4HmxVMRWbEh0dTWnSIL0qsA6k7gVApTErEyZMoBo1atDw4cNp0qRJJvfhma9MmeQHigDYwphfjtOyfddNbnuSNoQ+aTKIal46LFIUc/YvOSBUAPAs+vbtKx5ZqAwbNkxY+CUSExOF67GxazIAqYHUvQCoVKxUrVqVhg4dSpMnT6Zdu3bRgAEDRAYwaVZr06ZNopL9J5984sz+Ao1RfNRv9Dze8HaSNi6Gsj+9R/9myiVebyvwlljkqJgnlL7rWdXpfQUAOJYTJ04kW1Y4I5h+8hZ+zhnCBg0a5MYeAk8EqXsBUHmdFbacjB49mlq2bEne3t7idVxcnKi90qRJE7ENAEdgyu0rc/RjWrhmHEU8vU/N239BN8LCzbZxZlx9Susvb20BACg/C1inTp1o1qxZFBIS4u4uAZWA1L0AqLiCff/+/UVWMLairFu3TqQujo2NFdu43kpYWBiVKFFCzHhxnAs/VqhQwRl9BxoTKgXvX6UlP0VSjqf36WHaEMr4IsqsWIHbFwDqYMmSJe7uAlAhSN0LgErFCpMtWzYaOHCgWNhv+NSpU3T48OHk5fjx4+JR8jXmfQCwR6hUuXKM5v08iULiXtDFDNmpU8tIupZefuYLQgUAAIClqXs9rU4MBBbQEjaJFX18fHyoVKlSYmFTPZOUlESnT5+mQ4cO0ZEjRxzRT6ABomMSqETk7ynWtzq+hSb9Pof8khJpf47i1KPFCBFYLweECgDAWpYtWyaSx1y5coXy5MlDHTp0EBNyyDYGlMTmk7dTuK5xsgC4rgE145SqeBzLwq5gHTt2pNmzZzvjI4DKaDxrl0mh0vT0X/T5pllCqPxSrAa1bz0BQgUAYFAUkifI7IETx/D9qnbt2vTjjz/Su+++S5GRkdSnTx+H9RMARwgVTgpgXNCSs5rxet4OgBqx27ICgL0UHrGRYmU8BbcUeIv+CS9If+UtR9OrtWO/Qtl2IFQA0B5sBSlevDh99dVXotaXtXCCGBYrnDRm3rx5Yl3Tpk2F1wCvHzNmDOXMmZM8CbgJeS6mvjtpPVtUTKVa5nX87fJ2jsHBdw3UBsQKUFx8SlDsC3run1YIk5f+aej9tp9RrO/rlKXGtHozG33+Xjkn9xQAoEQWL14sihFzCuP9+/db/X52+3r27Bm99957BusrVaokHm/evOlRYgVuQp6L3Hc3qk4EpX/mk8KiYixYeDsLHU+LwQEgNSBWgKKESq7Ht0XGr3XFatLsKh+IdeaEyvkJDcjf1ynejAAABbl6sWAwrmLP9VfYdStXrlw0duxYm9rOmjWrSJFcunRpg/WcKIYpUKCAyfdxFkwpE6b0OiEhgZTgJmQ8+y65CXGqXggWZVpStp6+Q4v2XEmxnb+7L7f9S5VLWObqyBYZANQGxApQjFApd/MMLVgznjK+fErvH99Ci99oRtEBrytWGwO3LwC0Qd68een27duUJUsWg/WPHj0S2+zJOBkaGko1a9Y0WMfZLD/77DNq27atqCVmCnYRMxZI9erVo/v37xsUr5TgWmTOHvQu//M4FQ41PahlmcfbS2fyVpybkLPPjVI5dOUhrdh3jR49jxOvi4Sa3i9HkI6O/3uNish8t/qE0ku6d+8eaQWt/nacdW7YLVaJQKwARQiVhmd304xfp1FAYjwdDy9AXd4bDaECAEi2oBhbVZjo6GiHZuviz/n6669pwIAB9Oabb4o4GDmGDRsm9tO3rHBCmcyZM1NAQIDJ9xiLLUey9+JD+vsWW3bMWJqjEujKc19Fugk589woEbaC9V576T8rWGreAUl0NiqJMgSlocfP40zGrXj9V9CyWqn8ihOjzkZrvx1nnht9a7GSgFgBLoNn/vIP/81wpU5HPQ6soWE7loqXWwpUoL5NhohYFTkgVADQBn379hWPLFRYHAQGvp7AYGvKgQMHqEyZMg75rOvXr4v0+3/99RcNHz6cRo4cSX5+frL7syDRFyV8k/f1dd8t1VL3H7gJuR9zwfLmeLdMBC3Zc0UIE/33StKE45K0JlTUCpJkGAKxAlzCot0XaPyv51OsH7PtG+p0eIN4vqR8ExpfuyslefvItgOhAoB2OHHiRLLF48yZMwbuVfy8XLlyIrDeXrjtOnXqCJcvdgHjumGeBg9oHLkfcB48CDUXLC8HZ/riQatxED5bVJBAQT0gSUZKIFaAW9y+JE5nyUuJXt40oXYXWvJGM9n9quULpW+7V3VSDwEASoQD3xm2eMyaNYtCQuRrLNkK12hp164d5ciRg/78808KDg4mT4QHsTyg4YBsc25CUipc4D6stW7xd8ffrTS7zqIFs+7qBEkyTAOxAlwvVHS65HopP5Z6m45kL0oXM8qnBj0zrj6l9Ze3tgAA1M2SJUsM4lQYR4mKLVu20JEjR0QRyH379qXYXrlyZQP3M6XCg1WeeeUBDdyElI011i1T3x0/KjHuCNgHaunIg5yvwKVCpdjdS/T96mGU8fnrLBXmhAq7fUGoAABmzpwpUhRz9i5eOJXxjBkzhIuYPRw8eFA8sljhbF7Gy61bt8hT4BlXnnllC4o+/FqrM7JKtoJZMtzk7+5/dQrgu9MAqbkH6vRq6aQmejjhxt6LD8Qjv/Z0YFkBLhMqNS8eojnrP6PguJc0fMdiGtjodSYdUyA+BQDADBkyhL755hsaPHhwcrHGvXv30rhx40RK46lTp9rcNgfS86KWIFoe1MJNyLOtYPy6c5U8yTEqDx/cd2Nv1YkSA9gdkSRjs168C6e6Pht1WRXxLhArwKHEJSRRoZGbUqxve/Q3GrtlPvnqkmhP7lI0tk53s+1AqAAAJBYuXCiWli1bJq/jYpCFCxemHj162CVW1BhECzchz7GCIVje9Sg1gN3eJBmbVRzvArECHMaoX07Qt/uuGazz0iXRpzuWUs8Da8XrH0vUpeH1e1O8j3xKUAgVAIAxpjJ08ToOkPc01DyoAJbP7Iem9ach7xQWRSEzBAdQeIgyZvjVjJL/9+xJkpGo8ngXiBXgEAoM20gJRv8lAfGxNH3jdGp0bo94/UW1djSnUuvk4HpjQv296Z9xDVzRXQCAB9GhQweaO3euyAimz7x586h9+/bkSah9UAHsm9l393euRPcorfzv2ZMk44AV8S6eaHWFWAFOS02cNiGWit67TLE+vjSkQT9aV7yWbBtHRtajDMGvaygAAIB+wcVVq1bR77//ThUrVhTr9u/fL4Lf27Ztm1w8kuEq8kpG7YMK4Lkz+0p1j9LS/56t7oH3VF4UFmIFOK2GypO0IdSpVSRljX5EB3KWkN0Pbl8AAHOcPXtWFIBkrl69Kh7Dw8PFwgUdJbjSvdJR+6ACeObMvpJFlNb+92xJkpFF5UVhIVaAQ4VKhesnKffj2/RjqXri9dX0EWKRA0IFAGBpcUg1oPZBBfC8mX0liyit/u9ZmySjgsqLwqLOCnCYUGl2ajt9+/1ImrL5SyFaUgNCBQCgNVKrsaFfrRyoC6XO7DuqvofSUfP/ns9/8S6M8fGpoSgsxAqwiqgX8SmFik5Hff7+jmb9Oo0CEhPo94IV6Z/wgmbbgVABAFjD3bt3afTo0SJ9catWrWjMmDFinaeh9kEF8LyZfaWKKFf+70mibFQjz/3fq6/iorBwAwMWU2niH3T7WbzBOr/EeJq0eS61OrlVvJ5foQV9VrMj6bzkdTCECgDAGvbs2UP169enrFmzJheFXLlypahgz0H30jpPATU2tIlSXXWUKqIscV87fesp7b0Vb3HmMrn/PYnxG0+Tt/er/TyR+vrxLvfuUpYsWVWR0Q1iBVhEh0X76fYzQwESEhNN836ZRFWuHqdEL28aXa8nrSzbULaNNhWy05QWZVzQWwCAmhg0aBB98MEHNH/+fPLmkQSRqK/Ss2dPGjhwIP3999/kaaDSvPawJzWtFkWUJZnLQukFnY3ytipzGW9PStLRx6uOptimhoQCPv/Fu9xLl0hZsqgjo6Cq3MCWLVsmKhoHBASIx4kTJ1JMzGvlrNPpaOnSpVSwYEHy8/MTmWSGDh1KL1++dGu/PTXjV6Ozu4VQifZPS13eG21WqJyf0ABCBQBgE8eOHROiRBIqDD8fMGAAHT2acsDhaYOKZmWyi0cIFfXjKFcdtirsvfiQ1h27KR75tVZcE6XMZcaWEUlo8HZz8Lkav/F1FkF9pLPIQsiecwoci2osK5MnT6bhw4eLmbYGDRoIt4HIyEi6fPkyLVy4UOyzZMkS6tKlC7Vu3ZomTZokboBTpkyh27dvC6EDrEtNvLr0O5T96T36rXBVOp01n+x+cPsCANhDaGiouJbzJJQ+vC4sLMxt/QLAHVY1Z9RD8RTXREdkLlNqVjagcrESFxcnxAoHXnJFY6Zp06bk4+Mj1nMgZo4cOWjcuHFUtWpVWr16tcjHz0GaDAsX3pY7d243H4nyhUrVy0fpaERheh4QKCrRf1G9g9k2IFQAAPbSpk0bMdE0depUqly5sljHE1KffvqpcA8DwNOwNjWtK+qheIJroiOEhlYSCqgJVYiVK1eu0LNnz+i9994zWC8FXd68eVMIGi4mxm4D+oXDJCvLjh076KOPPnJ535VIXEISFRq5KcX6TofW0ahtC+mvfOWo63ujKdHbx2w7ECoAAEfAIoXdeDt37kwJCQliHbvy9urVS1jHtQjPMCt5UAk8sx6KrSLKVThCaHhqQgEtowqxwhliuGhY6dKlDdYfP35cPBYoUIDOnTuX/Fyf/PnzJwsaQDRszXFaffC6wTrvpERqd3QTtdu2QLy+nS5zqu1AqAAAHIW/vz/NmjVLWMovXryYfO0ODAwkLeIMNyCgfBztvuSJgtcRQsMTEwpoHV+1+DPXrFnTYN3hw4fps88+o7Zt21KmTJno5MlXRQqNb27Sa7a8mCI2NlYs+q+lmT0tuH2ljYuh2Rs+p9LBr87B5Jod6esK7wkXMFO8Vy4bTXu/nNP7CgDQHny9LlmyJGkZZ7oBAWXjSPclTxW8jhAaSs3KBlQuVvRhV4Gvv/5auHu9+eab9NVXX5ndX98lzBQ8kzd27FiDdQ0bNqQPP/xQzPapKTVxkVDDdaEx0dT/4ArKFxxL9wsWoi8rtaaDOYtTEfGvnfIysajjm+Tn40337t0jNfDkyRPSIlo8bq0ds9zkDFA2rnADAsrFUe5Lnix4jYUG2Sg0PCWhAFChWLl+/Tp16tSJ/vrrL5EZbOTIkcKvmfH1fXWo8fHxJm/a0n7GDBs2TAgffcsKZxXLnDmzSJGsHouKURZrnY7WrBhH5W+dpYdpQ2he7W70q19uoihtuX1lyZKFtIgWj1tLx6xvLQaeA7IYaRtHWBVcJXid6WKmLzSIXtgsNDwhoQBQmVg5c+YM1alTR7h8sQtYqVKlUsS1MBxkr8+lS5fEY0REhMl2WZDoixK+yUvCR9Wpib28RJHHKZu/pN7NhlJgpqyaEyoAAOVz7do1ke1RvwaLWkEWI23jCPclVwheV7iYSUJj1/GLFEVpbRYaSk8oAF6hiqs7VzJu166duGFxJWNjoSIF1mfPnp3WrVsnXMUk1q9fLx6rV69OWsv4ZUqoRDx97cJ1KrwANfloJl1LL39xgVABALiTPHnyiOQqu3btIrWDLEZAsipkDbGtqKSzBa+9BRutFRrFIkIUV1TVkQU73fkZSkIVJoItW7bQkSNHRBHIffv2pdjOefk5MHPUqFGiaCQvjRo1olOnTol4FM7hL2UF0wKj152g5XuvGa7U6WjArhXU/eDP9EGbSXQ0e5FX683E9ECoAADczeLFi0VxyEGDBtH+/ftJzSCLkbJxbXYtw1+A/iSsuwSv1mOq+Pjn/PkvLdlzmZ68jHeaVWmzhyZHIK2LlYMHD4pHFiumuHDhgrCsdO/eXbgKcF7+RYsWUYYMGah37940fvx40grFR/1Gz+MNLyX+CfE0ddNMevf0X+L1W9dPvhYrJggP9qV9I99xel8BACA1OnbsKB6NE6GoEWQxUi6uGkDKBcfffRprMjjeWECVz53eaYLXVTFV+scUSi8pY6bMbv/N8/cydO0JevLCMC7a0YkLNntwcgTSuljhQHpeUoMzf3Xr1k0sWsSU21fYy6f09dqJ9NaNUxTv7UPD3+lNP5Z6W7aNmS1L0btv5HRyTwEAwDAuJWfOnCmyN/JsMidWyZUrF2kFJWQx8sT6HM7EVQNIay0XcgKqaels9M3OyzYJXnPfvStiqoyPqUhoEg3bdNWtVgXuU88VR2S3O8qqlKhhy5UqxAqwTajkenyblvwUSfkf3aSn/oHUq/lw2pOnjGwbFyc1VN0/AABA+eTNm5du376dImPbo0ePxLbExETSEu7MYqRFFxRzg3RXDiCtsVxEvYyTFVAsVLpXz0vr/7ltleA19d1nCPKnCc1KUMNS2ZweU6VEq4L0/aeGsVXJFsF/QMPZACFWNCxUfv52IGV8+ZRupstMnVqNofOZ88i2sbzLWxAqAAC3wBYUUzWxoqOjKU0abQaTuyOLkRIHi67AnEALTevvsgGkpRaJO1Evaerv58wKKBYqfw2uRYevPrZowCz33T96HkcfrzpCPW7kpSH1izrNxUyJVgXu09I9l81+/8bwubZV8N/TcDZAiBWVI5ea+EZoFjqcoxhle3qfOrccQ/eDM5gNpFdLoUcAgOfQt29f8chChWtecaIUCbamHDhwgMqUkbcGA8ehxMGiEgRapyryk3yOHkBaapFgAWGJgGKhYomAMvfdS3y98zKVzhHmkJgqU1YHS60KLB4ypQtwurXRlOCwhCsPntPMrRdsEvxZNJwNEGJFpbyMS6SiozcbrtTpyEeXRInePpTk7UN9mwwib52OXvinlW0HGb8AAO7ixIkTyZYVrqXl7++fvI2flytXTmQBA87HkS4oxoPRPEE6jxVo647dctkA0tJscBmCAxwqoFL77iVGrjtJB0fUsyumSs7q0LBEuEV9Hb/xjIGL2rtlIoSAdqRwkROwqcHHsfrANZsFfwUNZwOEWFEh7Rf8TbsuPjZY55OUSOO2zKM08bE0sNEAkZI4xs/8xRNCBQDgTrZv3y4eO3XqRLNnz6Z06dK5u0uaxVEuKKYGo5UjfKlD7VKKcyGzRKA9fB5HGYL86PHzeKcPIC3NBseuaY4QUJKo3GRhbZRHz+Np38WHNsdUmbNiLdpzxaI+GPYnjhbvuSIWR8VVWWJlkoMzsf16/LbNgt9Hw9kAVVEUEhi6fRkLlaDYF7Top3HU9thman5qB5W9dc5sG75eECoAAOWwZMkS2r17NzVu3JiKFSsmMoAxCxcupG3btrm7e5rAES4ocgUDeVDp6IKBrhRozctkF4/GQ0RnDCClbHAsgOSKQkoz8HKfyOuzpSKg+Luo+tmf9MGCfbR871WL+9d71avvUYqpsrRgY2pWLIabsPUs8m/OEb8xS61MpjAnVCz93dW34PtXI7CsqDw+JfzpA1q8ZiwVu3eZXvoGUL8mg8zWUNkzpDZlzyDvFgYAAK5m5cqVophv165dhTiJj49PjluZOnUq1alTx91dVD32uqCkNiPN64etPaGomBdLBVrdYuH0Zt4MLksnnZrlwt4ZeFvdnBguhmhLsgVLRIBUpN34mKzB3rgqVwSvp/a7q+/GbIDuAmJFxUKl2N1LtPinSAqPfkT3A8OoS8vRdDxbIdk2YE0BACgRFiQLFiygNm3aCGuKRMWKFWn06NFu7ZsasCSNqr0DYEsGo49fxNOcPy9Qv7ry9ymlCjQ+blcOIFPLBmdrPR573JzsEQWWioDOVfLQppN3bLJu2JuZjc/Ng2ex5EzCQwIschn0MfH9e0osmC1ArKhUqNS4dJjmrptCwXEv6ULGnNSpVSTdCM0q2waECgBAqVy4cIEqVaqUYn1wcDA9ffrULX1SC9akUbWnIKWlg9Ele65Qn9oFFTFLbK1Ac0c6aUfPwO+79NBmNyd7RIGlViw+nhGNiiUfU9TDBzR6yw2nW0dszf5lLR9UyGX2+0mUmVjwpFgwW4BYUWlqYq5GH5AQR3tyl6Je7w6np2mCZduAUAEAKJmIiAg6f/485c6d22D9zp07KX/+/G7rl6djS90UW11QLB2MshuRqUFuXEISfbv3Cl199IJyZwik9pXykL+v88NuUxNofC72XnyoWHccSwUUD4LZqvX1X5cc9tnWiILHz+NS3UeKs9E/ptt3fGnegQey1i9HZGazxy3OWp6+fOXiakqcXHnwXGQTu/M01uCcNC2dTRT61Jk4p2qpfwSxoqbUxHr8nacMtW0zkY5GFKZ4Hz+T+/D19NJkCBUAgLLp3r27qLkiuYBxgP2uXbtoyJAhFBkZ6e7uKQJrK2LbUzfFeADMbaU2YOd1YWn9hBhJDW6H2+RZfm7373/v05HrUQb7TPztDHWrlpeGNSxGzkZOoG05fUcEoRuImJA0YnY8T6ZAt4oXa34PPBgfuvYEPXmR+ndjDZaKAu7r+I2pV4Ef1aioVe6J9mRmk87fnacxNP7XUy4RKsz3h67T8EavrHWWWHNuR8WIGjemkPqshvpHECseyEeL99Ff5x8arAuIj6VxW+bTNxVa0MVMOcW6AzlLyLZRvWBGWt6lotP7CgAA9sKiJCoqiurVq0cxMTFUq1YtCggIEDVWevfuTVrHlorYjqqbYuln80CpU5W8NGPr+VSP58qDF1R+whazg2cOtpYGaa4QLMYCTdYq9TTG4BgdlTLXWb8HW6wGjhQF1mTYSh8UYJX1y7hPlmZmc5XLlymiYxNpzp//UuHwYIdYc3R2xukoBaQu9jAKDd+YQqhkeBFFq78bTq1PbKFvfp4oaqqY48y4+hAqAACPYuLEifTgwQNRtX7fvn10//59Gj9+PGkduXTAd1JJ1eqIuinWfnaf2gUoLNC0pZ/hISRv58G+pbP8C3ZdFi5irsSaIPTUvgdHI/edmErda2swvc6C7aP+sw5Ywp2olw5J6bv709q0ultF6lIlj6h9Y21qXz4fs7ZeoJ4mzp+9sFisUySzRfsu2XOZItc71ppzzwVZzJwJLCseHp+S7+ENWvJTJOV+coeepAmm4e/0FhXq5UB8CgDAUwkMDKQ33njD3d1QDPa4ctlbN8XSz65dJCsdvvo42R1p0rsl6eNVR1K8x9Z0tGxh4ViWLtXykauwptaGJZXJHYW16aGtrRlizXfEbl3e3q9EhDlYPOlXnTdHar9ZyfrFC7tSWeMWyf1ggaAfD+Io+tQqQP3rFRL92Xb2fqr7C1dJy/SbxVgbp6M0IFY8WKhUuH6Svlk7gcJioulaaFaR8etixlcuYKaAUAEAeCrff/+9qLFy7949SkoynElfv349aRF7XLnsrZti6WdXnLxVVDaX4M/sUT0vrf/ndoqA9TZv5qQZWy+QtXDQvSuxdpbaVa441qaHtuU4LMVckgZbXNDSB/pZ7FZmbWY2ZwfQVymQSfSHK9i7Gi8rXfKUCsSKhwqVZqe209RNsyggMYGOZitMXd8bRQ+DwmTbgFABAHgqgwcPppkzZ4pYFc4M5uXluYGijsQeVy5766ZY+tn6QkUaxHLmorkflhUxCNxOKL2kaqXy06/Hb5EtcHYwV2LrLLWzXXGsTQ9tzXFYmhzBUouStS5ozhISjqorY4lQmLfjXyd8ymts+T/2FCBWFAz74RYauSnFei9dEn3wz+9CqPxWqDL1bzyQYv1MB57xz/MyhAoAwINZvnw5rV69mlq2bOnurigKe1257KmbYuuAXRrEsusPxxjwIIqtZfxoS5s8BuM0xq4kNauUM11xzGX5sjY9tDXH0alKHqutXuYsSta6oHEMU2qWKVPnRvosOXcwa/thKcZCgfv29c5LFgvDNH7edPdpbKrfi9d/j91NWCvTB/nTvCaoswKcyIhf/qGV+0wXOtJ5eVOP5iOEYPn6rRbitSmq5w+j5d2qOLmnAADgXNjtq0yZMu7uhuKw15XLnrop0mc7spI4t8mpfzmjlqVw+mJX1FuxNV2uI11xUsvyZW16aEuOgxMeTGlRUvxGvjt43abv25TFxxYrU2rJHozPjZTMQT9Zg3FWNGdZu4ICfOmLVq+FArvevYgzn/xIgv8nioSnE+Iwtd9XuF69n+oFs9DeSw/EL47byBecSNnC5YuBexLIBqZAioz8LYVQCYmJptb//J78OiptOppfsaWsUJnTpgyECgBANXVWVqxY4e5uKA5psMkYSwtrXEAk//5mZbKLR0tcRngfzvhkD8YDRW4zsqllbXIXOfbFFWmLzVmleLBoDltccaS6NeuO3RSP/NqSzGtSemhLkKww0nGEmsjSFhzgKxIi8D76vzVrMWXxscXKJPceuXPDIsU4q5xxdrZMMumQ7cXX20sICIa/P3a9s5RNJ+8IocJiy/h7yRaahvrXLUSz2pQRmc/YOslwvZ+2i/bTnO0Xac72f2nQj//Q0WuPSC3AsqLA1MRxRlkYs0fdoyU/RlKhh9fINymRVpZtaLaNi5Maerx/IgAASDx58oRWrVpFW7ZsoVKlSpGfn+ENfPbs2aRV7HHlshd2M7EHU4NP7u/8duVMFilM6+dNFfJmpOoFM7msgr05jK1SpiqMW/s9mLIQVMrmS2eivCzK+sbpoZf8fVk29bOclcfU/s9jE6j3qiM0z/tVkDwvX31YlvqsPiqysFmCVHXeGGtc0MxZpmyJfZHOF+fpGPdr6sUobUFytWPxz4/WxPtIRP33nfSvW0i2yKhsvZ+oGPpy27/kHRgGNzDg/ED6krcv0OI1Yynz8yd0JzgDHYsobLYNBNIDANTG6dOnk93Azp49a7ANwfa2u3LZiz0uNDxrLOcWJR2PVMGeh5iV8mWiihZafVyJcdYpDly39XuQG3g+fhFHT154W+xWx25bXCvEEiuPNNiXa5eMguQbloqgLtefiPo2lvAyPpG2nL6TYsBsqStdapapc3eeWe2aJp0vUym0nfH/Yev/iSSsvjt4LTm+y5r04WqpXs9ArChYqNS7sI9mbficAuNj6UzmPNSpZSTdCckk2waECgBAjWzfvt3dXVA8qaVqNReYbSv2BIyn9sncN075you7sebcWZMy19FZqaRB8SsLSDkaue4kPXoeZ9bKY0mAub4Q4n7+evy2VdYBuRTGllSe14/JYOHKx5gpOECMxh88j6Wbt6JIqUj/H/b8n9iTpECnkur1DMSKQoVKp0PraNS2heRNOvorbznq3WwoRQfIp2eEUAEAqInr16+Tj4+PSFUMyKmB2a7OiiXV+7B2EOUMweWuc2eMI7JSSYPiV4UWTxsIFa7oPqpR0RR9tnTWn60jfL6X7rlsVT8l68CIn0/Sy7hECg9Na/C96VsFObHCg2exwpLEmyVrGn82x2SY+twioUmKDL/Wd3+z5//EEUkK7nl49XoGYkWBQqXY3Us0ZtsC8XxV6fo0ul5PSvAx/VWl8fWisxPMx7AAAICnMHHiRJo+fbqIU2GCgoKodu3aIsi+YUNc66wd2PNAb7GJ4F5LivalJg6szYplzyDKlGjIEORP75aJEINdc8KFj8MWlzJz8QCpnTtrsWdAqR/TIetK9jyeeq86SvO8vQz6bOms//cHr9Mvx24ZCCBL4b48fB5H/X/4x6TY4+8h6mUcTd181uD7XXPkJjUtnU3U5XFWnRVnYez+ZmvBU4lMwQEpLEsX7kZrono9A7HiJu4/jaU3J201ue101nw0qWYn4nC6ryu8x07ZJverUSgTLev8lpN7CgAArhMqo0aNEs+LFCkihMrNmzdFhfoNGzZQgwYNaOXKlRQaGururioaUwN7a4v2WWpRsMSVx95BlNwAnAfOLMR4CQ8JoA8q5KI8mYIMhBW/1zhYnzMmSSl5zQk1c/EA5s6dLdg6oNSP6WCs6TMfY5JOR0EBPvQ81nxa3edxiWJxBMZiT+775d/T1zsti41RGpL7m6n6J/pkCPSjmIQk2bTG/J2FBfrRwB+OGSRusAQvMwkOPA2IFTdQYvRmijb6YWaOfkzeukS6m+6Vf+43b71nto0z4+pTWn8fp/YTAABcyYIFCyh9+vS0bds2g7oqBw8epGnTptEPP/xA77zzDu3cuZP8/e3LRKVW5AZ+1vi0W2tRMA7wv3T/Oc3aZn4WWW4QZWzNKZ87vUWxHDyQ05+55vZ5Vl5usMvihYPQOfOYKcFiaTzAjC3nRVxNapad1NzXLKlbw/OWIWl8Keplgsk4FJ55tzSGgS0ZtghMR6AvnGoXyeq0CvLuRDqe1MTWI5msbYxkrXxsZh9z71VL9XoGYkUBbl8F71+lJT+NpWcBgdSq7VSzsSkM4lMAAGqErSjs7mVcAPLNN9+k7777jqpVq0b/+9//hJvY0KFDSStYGqtha5D2pv9qTkjiwRaLghRYzn3g+ILUGN6waPIxhdJLypgps3CbSenq5UePnls/WLtj4ax85PpTJq0jlrplcU0LXtglbUKzEtSwVDabLFT8+ebEFaPTkRAq/esWTGFBenXMLy3qM59nrvvhToEgCadv915xi2DyBFiIsjvZExvECr/3f3UiVJG2mIFYcSH5TAiVyleO0fxfJlNI7HO6lD6CQmOiEUgPANAkgYGBlC5dOtntvXv3FqJl+fLlmhEr1gR42xqkvXzvVbFwu+xbb0+GIUv7MGb9qeT4Bw6S/mT9JZODMluEitRPS2CLjKljsdYti4+FU+H2uPG6UKU5CxVbdTpXyUN1imQVCnDbmbsmY4tMsXjPZfqqbflkcckWFRYgPx42LCYtB8eeKMWSsf8yxxEBY9IH+tH7b+RM1UJpTJ9a+alKgczit/HwwX1SC8pLoaBS8g/bSEa1Hqnlia207McxQqgcyFGMWrT/gm6GZpFtA0IFAKBmypYtS3/88YfZfdi6cvmyY/3YdTodLV26lAoWLCgKToaHhwsx9PKlZTPVzsKSquWOzPoj3JosDALe8+99YUUxZsGuixa93zhQ25bZY0fBliWpUrwED/Y4DsZa2DLCrmFxCUnCamOuBgaLE6463nbhfouFCsPWFX5P+QlbxPLBgn3i/c9iXruHmcIr2VJlfZC8s/jj9D13d0GRsOuXtUKF4dgXFt5qcP3SB2LFRa5fifpXLJ2OBuz8lr74bSb5JSXSuqI1qH3rCfQkbYjJ9/t5QagAANTPhAkT6J9//qHhw4fL7nP//n3KlMmxtTeWLFlCnTp1ovLly9OqVauoS5cu9Pnnn1PPnj3JXfDAOXJ96gXf9AfYrsz6w0Hq7O6lL5gm/3aa/jzrebO5bFXiAb/+8fBg78086W1qjweZZcf9YXVAtLWwwLNG5PEvJSgADjVq5pdjt0xOIng6ECtOpuDwlK5f/Xevor57vxfPv6zUmj5pMpBifU0Hi9YomJEuTIZQAQCon6pVqwqLxpQpU4QF5eeff6aYmNfWgk2bNgkx0aFDB4daVcaNGyc+e/Xq1dSqVSuRlYz7we5mV69eJXcw588Lou6EJe5YxkHarppT1bfwsCXB0qrmSkX/eHjAt+uC7S5Kjsqc5WiuP3KvtRA4l0fP4wyuCWoBYsWJNJi5g+KNfb+IaHXpd+hGSBYa3KAvTavennRepr+GLz8oS8u6VHR+RwEAQCGwUODMX8eOHaOWLVuKGJZs2bJRxowZqXHjxlSvXj0aPXq0wz7v0qVLQpCwSPHSSxPfunVr8bhjxw5yNTxYttQdS9/1S6p7wsgJFk557yj0LTzL/r5Mnj6hq388XJflyUv3uaYBYCv3VFAE0hiIFSfAMzK9lh+kM3eeJ68Lin2R/PxOSCaq020+/Vjqbdk2Lk5qSE1Ko3IzAEB79O/fn86fP09Tp06lSpUqiQKRjx8/FlYQrrcSFhYmMoT16NGD5s+fTwcOHLD5s+7cuSMeCxQoYLA+f/78yRnKXImU0ctSjF2/pLonnA1IH7a4cJrenjUMj9McHN+QGpKF5+CVx6REmpQKF3UqLEU6nlcFJAHwPLKooAikMXBedMKM2IAfjtGLuNcmlXI3z9DXayfSmHo96bciVcU6ObcvrpxyEfEpAACNw9aUgQMHiiUxMZFOnTpFhw8fTl6OHz8uHhm2iPA+thAfH5+ciUwf6XVcnOlg5NjYWLHov05IMB/gbAnWZPSSq1ViXPdEP8Uti6GwtH4WWQ1GNS5OF+9Fi9S8qRFkYd2vdGl8KTomwaZsVFLdCWuoWyycZrYpl1zB/vzdZ/TH6bsWvNPDzUQeCv8+UksUAOThFNpqKAJpDMSKg4UKpyPUp8HZ3TRj43RKkxBHnQ6tp98KV5GtSB/k50Wnxjd0UW8BAMAz8PHxoVKlSomFA+GZpKQkOn36NB06dIiOHDG87joCfZcwU0yePJnGjh1rsI5d1DgBgKmClWwdsoR79x6IVL6W0LdOhNn0pPnT8cJWhUSD/Xq9lYl+Ppq6xSi910t6I6uPRf1pXCiIzl5JErVA5OBT+nGNXPTVjn8NpEDOIH5l+jM4ILxL1Tzi+Yp916zOZMU1XPjYC4YQFSybnk7f8qFrNw2zqJmifBYfqhzhS4+fx7ldtpg7P2qje/UclCEogE7fiqL1/9zyyHPDv/NeNQvQdweuufz3807x0OT/dUuvOfrITc64G4gVB8GzVZ+uOf56hU5H3Q+speE7loiXWwpUoH5NBssKlRLZQ+jX/1VzVXcBAMCj8fb2phIlSoilY8eONrfj6+trYGExvmlzKmNTDBs2jAYMGGBgWZk9ezZlzpyZAgJMp7zNkkU+NX3yPs986GxU6oHq/esWooYVCpItdH07M8078EA2kxTfpdiNrFqpV65wwzZdFcHnOjP71ixTkA7cSTBb1LBH9bzUpFJR8ksXRsN/PqFXQyWJzkaZ9kpf2fVNUSGeeeeNIsnWogyB/tT3uyP0+EVCqsegn8aVi09acjzVSxegWL9gEXDPuFewyJ8fJWCLxUuObOHhIvVutVI6+vafJ7Lfk5LPzVcflqWGpSLIPzjM5b+fMcXzU5YsGa265uijby1WEhArDmLOn/+K3OeMT1Iijd0yn9od2yReLynfhMbX7kpJ3qbN5Ccj36HgNPgqAADa4tq1aw5ri+NYQkJMp383R9asWcWjcdYvDrxnIiJMxw6yINEXJXyTl4SPPUgZvcwN0rj+R5/alseeGMOD9yktSposWCgN6zlQXxrk83Pe13hQaryvVAyRs4LpB9tzM92qvS6WyG5qtYtkpYqTt8laSiTRUDHf64EXf4Z+8cbJLUqZHAyaOgb9Niw9Hin+x1RRzqiX8aKmhTVIn8cxNNakHJYz8nHhQG5Pv63gAB+Kjk20W1DI9TFjkD81LRNBOcLSCpej8NC09PBZLPX57qjZ9njfxMQkijLj4qXv1mjue7K0z0EBPiLOTd8t3xFwP0c1KkrjN54xW6xV7vfjLLLJuIWqAYyQHWRVWbLn1WySb2ICLVg7nmpdOkxJ5EXj63SlJW80M/k+P2+iC5MQnwIA0CZ58uRJ1d3KUsaMGWNTljAOrM+ePTutW7dO1FeR+rN+/XrxWL16dXIllgymI5sWt7vom9xAKtxowGXtvixIBr5dhL7de4WuPnpBuTMEUvtKecjf13D2m19Pal4i+TjJQrFh6zHYc+ym4n+4Yryx27c+gf4+4hj1B89S+/rtZQoKoINXHtLSv68axBFxcoPmZbKLmJtcgfF07YUf3Yl6KcRdhuAACg95PTA11Te5Yzt67XEKMWk84GchK/WR02c/io5NFiZS7JMxJ249kbWq8d78XTPWiEtLB/vGfdY/F8y+iw9p76UH4pP4Z/jdwes21cDR7yf37Z0S2UzGhelj/Pu58uA5rT5wzeDzWQCWzhlKx65H2Vyw08uC/xdPxkvHshNYDM+ezZ07l3r37p08q8ZBe1xQSmLEnwup3dFN1K/JIPqjUCWT7RQND6JNn9QkT+HevXtWmxM9HS0es1aPW2vHzNcxrmXCtUTkXJZcQWRkpMPESs2aNalGjRo2vffrr78WBSC7d+9OjRo1EsH8HI/SvHlzUXvFEefU2t8Yx0CamtE3Nwi3dbIttQGXLftac5zL/zxOf99KsPk4be2XvcfDfecq9foDz9C0vtS5Sr5ky5cjzq0t1ydz7XFNHBaTVx4+FzFGHNQurFb5MlFFO6qf/3b8No1cd9JgwG38Xdryu9Y/FhZ3STod7Re1RHQitohd9qzps9Sesfh7+CyGRm84peeeaIij/v/kvhtpPYtNLu6ofx6DA3zIm3T0NDbJ4n7Z8rtRyr3BGIgVB4iVdcduUr/vjiXv46VLoryPbtGljDlMtvFR5Vw0tmlJ8iS0NpjT6jFr9bi1dsxKvSG5C74NLly4UJwTdgfLkCEDtW/fnsaPH58iS5irxIqzxIESuX3nLl157uuRx+mK78iTrk+WnA9HnjNHnxtjYcQmiwfRsS7/XZo6R4wpkSXXLzWJFV+t3ZCWLVsmio5duXJFFBnjwEx2H0ibNq3N7RY6spsW/jSNPn53OMX5+okij3JChf12RzR65bcLAADA/bB1p1u3bmJREsYxGmrFk4/Tk/vurvOh5HOmlL7J9aOSAvrmDpSVQsHJLFmyRKS9LF++PK1atUr4J3/++efC/G8zX39NRbq3pboXD1LnQ+tkd/P6ryI9hAoAAAAAAACW4aslq8q4ceOoatWqwgeZZ9JatWoltk2aNElsy507t+UNJiURffop0dSpQojcbPo+LS7UTDZrxdz/UtkBAAAAAAAALEMzlhVOQ8m+yCxQ9AM6W7duLR537NhhXYNcmGzq1FfPx42j7L98R7M/ektk3DAOfJrfrhyECgAAAAAAAFaiGcvKnTt3ktNU6pM//6uiVzdvpl7N14BffuFqYUSLFxO1a2c2xaGnBAoCAAAAAACgJDQjVqTqxMaZXaTXUrViU5kR9Ct6xsTEUEJCAsVmzky0dClR1aq8k8F7yuUI5kRz4nlCfBzJl0DyHPj8KLWyqbPQ4jFr9bi1dszSsSIZpOOQzqXc70hrvzFrwLkxD86PPDg3jj03Sr03aEasyJFajv/JkyeLfPv6cGzLs44dyXf3biJeVA6Ls127dlG1atUcUqHZE9DiMWv1uLV8zB9//DGlSWPougpsQ5rwmjFjRoptWvyNWQrOjXlwfuTBuXHeueHrmZLuDZqps7J7927xpf3xxx9Ur149gy+Ec0lPmDCBRowYkapl5cmTJ0Ks3L17l0JDQ0kLPH36VOTq5pzdISEhpAW0eMxaPW4tHnNUVBRlzZqVHj9+TGFhYe7ujipISkqi6Oho8vf3TzEJpsXfmKXg3JgH50cenBvHnxuWBDwuDg4OJm9v5YS1a0aK8o2Z4SB748B7JiLCdAA8CxlThXFYcSqpYI4zkY5T7lyoES0es1aPW4vHLM2YKelm5OnwuZQbFGjxN2YpODfmwfmRB+fGOedGSRYVCc3cqTiwPnv27LRu3ToDX7z169eLx+rVq7uxdwAAAAAAAADNWlbYLD9q1ChRAJKXRo0a0alTp0Q8Sps2bZKzggEAAAAAAACUgWbECtO9e3dhqp8yZQotWrSIMmTIQL1796bx48db3Aab08aMGaMpkyOOWTto8bhxzMDZ4HzLg3NjHpwfeXButHNuNBNgDwAAAAAAAPAsNBOzAgAAAAAAAPAsIFYAAAAAAAAAigRiBQAAAAAAAKBIIFYAAAAAAAAAigRixUI4D8HSpUupYMGC5OfnR+Hh4TR06FB6+fIlqYVly5ZR4cKFRfYIfpw4cSLFxMRo6hxw9gzj6tNqPO4jR45QvXr1RBG7TJkyUcuWLenatWuqPWauLj5z5kwqUqSI+H1zzaW+ffuKiuNqO+YFCxZQnjx5Uqy35PjUcg486fqiZVK7DmkVS65XWsSea5uWz48l4zvFw9nAQOosWrSIs6bpWrdurfvhhx90w4cP13l7e+s6dOigUwOTJk0Sx9ezZ0/dunXrdEOGDNH5+vrqunTpoplzsH//fp2Pj484Rn3UdtynTp3SBQYG6mrVqqX78ccfdUuXLtXlzZtXV7JkSV1iYqIqj3nUqFHieD755BPdL7/8ops8ebIubdq0urZt2ybv48nH/PLlS92VK1d0c+fO1QUHB+ty586dYh9Ljs+Tz4GnXl+0iiXXIa1iyfVKKzjq2qbl8zPJgvGd0sFV0wKSkpLED6Bq1ariuQT/Q/APgH8onkxsbKwuXbp0upYtWxqsHzZsmDi+a9euqf4cPH/+XFe4cGFdw4YNDQYTajzuRo0a6YoVK6aLjo5OXrdz505xnKdPn1blMYeFhenatGljsG7kyJHihsbfvacfc79+/UQ/pcX4hmXJ8Xn6OfDE64uWSe06pGVSu15pCUdc27R8fmItGN95AnADs4BLly7R1atXqVWrVgYm/NatW4vHHTt2kCdz5coVevbsGb333nsG6ytVqiQeb968qfpzwCbjMmXKiOPTR23HHRUVRZs2baKePXtSUFCQMJ+zy0G1atXE76Bo0aKqO2aGi8GGhYUZrGPXEz7++Ph4jz/mTz75hPbu3SuWLl26pNhuyfF5+jnwxOuLVrHkOqRlUrteaQlHXNu0fH6uWDC+8wQgVizgzp074rFAgQIG6/Pnz+9RX7YcWbNmpe3bt9M777xjsP748ePJx63mc7B161b64YcfaM6cOSm2qe24z5w5IwYFuXLloubNm1PatGmFDyt/9xcuXFDlMUsX9OXLl9PPP/9Mjx49ot27d9OsWbOoa9euFBoa6vHHzH7KFStWFEuOHDlSbLfk+Dz9HHji9UWrWHId0jKpXa+0hCOubVo+P1ktGN95AhArFiDNZAQGBhqsl17HxcWRJ8MXv5o1a1L69OmT1x0+fJg+++wzatu2rQh8VOs5ePLkCXXq1InmzZsnjtMYtR33vXv3xGOvXr3E4GDNmjX01Vdf0T///EP169en2NhY1R0zM2jQICpfvjy1aNGCMmbMKGZw+XjGjx8vtqvxmPWx5PjUfg6UeH3RKpZch7RMatcr8Bpct+wf33kCECt2oMasLmxmnj9/vrg4li1bVtxA1HwO+vTpQzVq1BCze9bgqcfN5mCGs6asXr2aGjVqRN26daOffvpJmNO//fZb1R0z/6abNWsmLtB8s9+yZYvItPPw4UPxO3/+/LnqjtlSLDk+tZ8DJV5f1I491yG1Y8/1CrwG1y37x3dKwtfdHfAEfH1fnSZjX1FJsXOqPDVw/fp1MQv4119/0fDhw2nkyJHJx6bGc7B27Vr6448/6NixY8kp/KTj49c+Pj6qO27JD5pvhvoX8ypVqlC6dOno0KFDIl2mmo6ZfXn5hs8X6R49eoh1devWpUKFClHDhg3F7yBv3ryqOmZjLPkdq+237gnXF62eU0uuQ+zypEUsuV61b9/e3d1UDLhu2T++8wQgViyAff4YDuLSh2eAmIiICFKDD3GdOnWESZBndEqVKqX6c7Bz5066f/++yGFvDPtQ9+vXj3r37q2q45Z8WnmGxRhexwMFtX3XHGDIVK5c2WA9D4yYGzduCH9fNR2zMZZ8p2r73j3h+sIz5lrEkuuQVrHkegVeg+uW/eM7TwBuYBbAAUh8w1m3bp3BxXX9+vXisXr16uTJcKBju3btxA3k77//NvlDVuM50M+iIS0828Dwc96utuMuUaKEuHj/8ssvBsfDQcBccIz9pNV2zFKRLA5SNR5MSq4oajtmYyw5PrWfAyVeX7SKJdchrWLJ9Qq8Btct+8d3HoG7cyd7CvPnzxc5qbt37y6K6nCRnYCAgBS50D2RzZs3i2OLjIzUbdmyJcUi5XVX8zmQWLJkSYo6CGo77sWLF4vjadq0qe7777/XzZgxQ5c+fXpdmTJlRE52tR0zF5jjmg5cVI1rFXCRtSlTpohaBpUrV9bFxcWp6pjHjBljsjCYJcenlnPgSdcXrWLJdUiLWHq90iL2XNu0en42Wzi+Uzq4aloIFxv65ptvdPny5RNViDNnzqwbMGCAx3zR5hg/frxBUSHj5cKFC6o/B+YGE2o87mXLlumKFy+u8/f312XKlEnXuXNn3f3791V7zE+fPtWNGDFClz9/fl2aNGnEcXF16IcPH6rumOVu6JYcn1rOgVKBWLHuOqRVLLleaRF7rm1aPT/jLRzfKR0v/uNu6w4AAAAAAAAAGIOYFQAAAAAAAIAigVgBAAAAAAAAKBKIFQAAAAAAAIAigVgBAAAAAAAAKBKIFQAAAAAAAIAigVgBAAAAAAAAKBKIFQAAAAAAAIAigVgBAAAAAAAAKBKIFQAAAAAAAIAigVgBAAAAAAAAKBKIFaBadu/eTTVr1qSQkBDy9fWl3Llz04QJEygpKcndXQMAAOBmcI8AwDOAWAGqZMSIEVS9enV68uQJ9ezZk7p160ZPnz6lUaNG0cSJE93dPQAAAG4E9wgAPAcvnU6nc3cnAHAkM2bMoAEDBtCQIUNoypQp5OXlJdYfPXqUypUrR2FhYfTw4UPy9oZWBwAArYF7BACeBcQKUBW3bt2iQoUKUalSpYSJ3/hmU7RoUTp79ixduXJFmPyBc2ncuDFlypSJli5d6u6uAAAA7hEAeCCYNgCqYvXq1fT8+XMxa2ZqViwoKEg8+vj4OOwz2ee5T58+DmtP6Z8LAACeCu4RjtsfAFcBsQJUxbp168QNqGHDhia337x5U9yMIiIiXN43TyQuLs7dXQAAAIeBe4T94L4AXA3EClANCQkJdODAAeF2FBgYmGL7wYMH6c6dO1SrVi1xs+JZpI8//piGDx8u3pMlSxYaNGiQQSYY9pKcOnUq5c+fn9KmTUslS5akFStWJG/v2LEj/fXXXzR37lzh98wLuw+YIzY2lj755BPKmjUrpUmThipWrCjcEfThvnHQZ79+/Sh9+vRiGTx4cHLfbPnc1NqU9unVq5c4D5kzZ6YqVapYdB6YFy9eiH4FBweLY5s0aZLZ/gAAgCvBPUL+c83tL3dfMGWJ4XbY/dfS8wOARXDMCgBq4NixYxx/pQsLC9MlJiam2N66dWuxff369eJ1jRo1dCEhIbpRo0bpzp07p/v+++91Pj4+ulWrViW/Z/jw4bpChQrpNm3apLt06ZJu5cqVusDAQN2vv/4qtj958kRXqVIlXadOnXS3b98WS0JCgtl+9u3bVxceHi7aOH36tK5r1666oKAg3a1bt5L34b4FBwfr+vTpoztz5ozoG/d12rRpNn9uam3q7zNgwACxD/fPkvPA9OrVSxcREaHbvHmz7sSJE7qWLVvq0qVLp/voo49S/e4AAMDZ4B4h/7nm9pe7L/D63r17G7TD1/tGjRpZfH4AsASIFaAaFi1aJG40vGzcuNFg25w5c8T6+vXrJ6/jC23FihUN9qtbt66uS5cu4nl0dLQuTZo0up07dxrs069fP12DBg0M2jG+YMvBbfr5+emWLVuWvI5vCPny5dONGDHCoM2CBQvqkpKSkteNHz9elz17dps+15o2S5YsmaLPqZ2HZ8+e6fz9/XUrVqxI3s7rQkNDIVYAAIoA9wjzyO1v6r5giVix9PwAkBq+ltlfAFA+hw8fFo9NmjShli1b0ocffihM1rt27aI9e/bQG2+8QatWrTJ4D2eE0Yf9lO/duyeenz59mmJiYqh+/frJqS2Z+Ph4ypMnj019vHjxoni/ZEaXAjkrVaokPk8fNv3rfy7vwzUAuBYAFzGzBUvaLF++vMF7LDkPfFzsx8ztSbA7GJv8AQBACeAeYTvG9wVLcMb5AdoEYgWo6kbk7+9P33//PY0bN46WL19ODx48oLx584qqxAMHDhT+v/r4+fkZvOYLquTzKz1u2LCBcuXKZfZ9jkD/Yu5OpGw4Eq4+DwAA4Axwj3DcfYHhuB7j6hcsRCRw7wCOAmIFqILExEQ6fvw4lShRQgTxTZ48WSz2UKxYMQoICKCrV69S7dq1Zffjmx9/viVwkCHvz7N4/Fzq+969e8Usnz779+8XNwLpBrVv3z4xqyfNmFnzuZa2aet54GPhmw+3ly9fPrGO04OePHky+TgBAMBd4B6ROtbuz1ap27dvG6z7559/kq0mlp4fAFIDYgWoAjY3v3z5UlQfdhTp0qUT2U944RtC9erVKTo6WtwQeEape/fuYj++MHOGGc6cwq5PGTJkkK18zLNTnFXl008/FdlleEaPqynfvXtXZJ0xLl7GGWF4/YkTJ+jzzz+nkSNHJm+35nMtbdPW88Cf36VLF3FcfAPjGybPXForpgAAwBngHpH6PcLa/VmA8OevX7+eChcuTF9//TVdv349WaxYen4ASA2kLgaq8kUuW7asQ9sdP348RUZG0hdffEHFixenevXq0Zo1a8QNRIIvxDwjxbNIPFC/du2aqNgulyrys88+o9atW1OnTp2oTJkyYrZv8+bNlC1bNoP92rZtKwb7b731FnXr1k2Igf79+9v8uZa0ac954G2c8rN58+bikWcw+eYEAADuBveI1D/X1P7m6Ny5c/LCMTYsTvj6b+35ASA1vDjKPtW9AFA4//vf/2jOnDnCVM5Bh+5mzJgx9NNPPwmTuK+v9QZMzl/Pg30+Jkd9rq1tAgCAp4N7hGM+FwB3AMsKUAVffvmlMDMr4SbE/Pbbb6K4lqtvBu76XAAAUDK4R7j3cwGwB/xaAXACXAlZS58LAADAcnCPAMBy4AYGAAAAAAAAUCRwAwMAAAAAAAAoEogVAAAAAAAAgCKBWAEAAAAAAAAoEogVAAAAAAAAgCKBWAEAAAAAAAAoEogVAAAAAAAAgCKBWAEAAAAAAAAoEogVAAAAAAAAgCKBWAEAAAAAAAAoEogVAAAAAAAAgCKBWAEAAAAAAAAoEogVAAAAAAAAACmR/wO2y4mxb3bL4wAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# compare the SNRs\n", "plt.figure(figsize=(8,4))\n", "\n", "# Left panel: rho_pred vs rho_true\n", "plt.subplot(1, 2, 1)\n", "plt.plot(bilby_snr['snr_net'], interp_snr_aligned_spins['snr_net'], 'o', color='C0')\n", "# diagonal line\n", "max_ = 100\n", "min_ = 0\n", "plt.plot([min_, max_], [min_, max_], 'r--', label='y=x')\n", "plt.xlabel(r\"$\\rho_{\\rm net, opt, pred}$\")\n", "plt.ylabel(r\"$\\rho_{\\rm net, opt, true}$\")\n", "# plt.xscale('log')\n", "# plt.yscale('log')\n", "plt.xlim(min_, max_)\n", "plt.ylim(min_, max_)\n", "plt.legend()\n", "plt.title(r\"Predicted vs. True SNR\")\n", "plt.grid(alpha=0.4)\n", "\n", "# Right panel: SNR difference\n", "plt.subplot(1, 2, 2)\n", "idx = (bilby_snr['snr_net'] > 4) & (bilby_snr['snr_net'] < 12)\n", "plt.plot(bilby_snr['snr_net'][idx], abs(interp_snr_aligned_spins['snr_net'][idx]-bilby_snr['snr_net'][idx]), 'o', color='C0')\n", "plt.xlabel(r\"$\\rho_{\\rm net, opt, true}$\")\n", "plt.ylabel(r\"$|\\rho_{\\rm net, opt, true} - \\rho_{\\rm net, opt, interp}|$\")\n", "plt.title(r\"SNR difference\")\n", "plt.grid(alpha=0.4)\n", "plt.title(r\"SNR difference for $4 < \\rho_{\\rm net, opt, true} < 12$\")\n", "plt.tight_layout() # Adjust layout to prevent overlapping titles/labels\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Maximum percentage difference between interpolated and bilby SNRs: 0.08%\n", "Maximum absolute difference between interpolated and bilby SNRs: 0.1853\n" ] } ], "source": [ "# percentage difference\n", "percent_diff = 100 * abs(interp_snr_aligned_spins['snr_net'] - bilby_snr['snr_net']) / bilby_snr['snr_net']\n", "print(f\"Maximum percentage difference between interpolated and bilby SNRs: {np.max(percent_diff):.2f}%\")\n", "# absolute difference\n", "print(f\"Maximum absolute difference between interpolated and bilby SNRs: {np.max(abs(interp_snr_aligned_spins['snr_net'] - bilby_snr['snr_net'])):.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note: There can be indeed some outliers with higher SNR difference, which correspond to the edge cases with high total mass and/or low mass ratio. Events that have negative chirp-times, and likely merge outside (below) the frequency band of interest, can also be found in the outliers." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ANN-model and $P_{\\rm det}$ Estimation \n", "\n", "* for more details on $P_{\\rm det}$ calculation, please refer to the [Pdet Description](https://gwsnr.hemantaph.com/probabilityofdetection.html) page." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "psds not given. Choosing bilby's default psds\n", "Intel processor has trouble allocating memory when the data is huge. So, by default for IMRPhenomXPHM, duration_max = 64.0. Otherwise, set to some max value like duration_max = 600.0 (10 mins)\n", "ANN model and scaler path is not given. Using the default path.\n", "ANN model for L1 is loaded from ./ann_data/ann_model_L1.h5.\n", "ANN scaler for L1 is loaded from ./ann_data/scaler_L1.pkl.\n", "ANN error_adjustment for L1 is loaded from ./ann_data/error_adjustment_L1.json.\n", "ANN model for H1 is loaded from ./ann_data/ann_model_H1.h5.\n", "ANN scaler for H1 is loaded from ./ann_data/scaler_H1.pkl.\n", "ANN error_adjustment for H1 is loaded from ./ann_data/error_adjustment_H1.json.\n", "ANN model for V1 is loaded from ./ann_data/ann_model_V1.h5.\n", "ANN scaler for V1 is loaded from ./ann_data/scaler_V1.pkl.\n", "ANN error_adjustment for V1 is loaded from ./ann_data/error_adjustment_V1.json.\n", "Interpolator will be loaded for L1 detector from ./interpolator_pickle/L1/partialSNR_dict_1.pickle\n", "Interpolator will be loaded for H1 detector from ./interpolator_pickle/H1/partialSNR_dict_1.pickle\n", "Interpolator will be loaded for V1 detector from ./interpolator_pickle/V1/partialSNR_dict_1.pickle\n", "\n", "\n", "\n", "Initializing GWSNR class...\n", "\n", "psds not given. Choosing bilby's default psds\n", "Intel processor has trouble allocating memory when the data is huge. So, by default for IMRPhenomXPHM, duration_max = 64.0. Otherwise, set to some max value like duration_max = 600.0 (10 mins)\n", "\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/phurailatpamhemantakumar/anaconda3/envs/gwsnr2/lib/python3.10/site-packages/sklearn/base.py:442: InconsistentVersionWarning: Trying to unpickle estimator StandardScaler from version 1.7.0 when using version 1.7.2. This might lead to breaking code or invalid results. Use at your own risk. For more info please refer to:\n", "https://scikit-learn.org/stable/model_persistence.html#security-maintainability-limitations\n", " warnings.warn(\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import gwsnr\n", "\n", "# it will take a while to gerenerate interpolation data for precessing systems.\n", "# interpolation is used for the partial-SNR computation.\n", "# partial-SNR is one of the inputs to compute the full SNR using ANN models.\n", "gwsnr_ann = gwsnr.GWSNR(\n", " npool=4,\n", " snr_method='ann', \n", " waveform_approximant='IMRPhenomXPHM', \n", " gwsnr_verbose=False,\n", " # ann_path_dict='./ann_data/ann_path_dict.json', # ann_path_dict can be used if you generate your own ANN models\n", ")\n", "\n", "gwsnr_bilby = gwsnr.GWSNR(\n", " npool=4,\n", " snr_method='inner_product', \n", " waveform_approximant='IMRPhenomXPHM', \n", " gwsnr_verbose=False,\n", ")" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████████████████| 10000/10000 [00:14<00:00, 685.96it/s]\n" ] } ], "source": [ "# gerneral case, random parameters\n", "# add random seed for reproducibility\n", "np.random.seed(42)\n", "nsamples = 10000\n", "mtot = np.random.uniform(2*4.98, 2*112.5,nsamples)\n", "mass_ratio = np.random.uniform(0.2,1,size=nsamples)\n", "param_dict = dict(\n", " # convert to component masses\n", " mass_1 = mtot / (1 + mass_ratio),\n", " mass_2 = mtot * mass_ratio / (1 + mass_ratio),\n", " # Fix luminosity distance\n", " luminosity_distance = np.random.uniform(100, 10000, size=nsamples), # Random luminosity distance between 100 and 10000 Mpc\n", " # Randomly sample everything else:\n", " theta_jn = np.random.uniform(0,2*np.pi, size=nsamples),\n", " ra = np.random.uniform(0,2*np.pi, size=nsamples), \n", " dec = np.random.uniform(-np.pi/2,np.pi/2, size=nsamples), \n", " psi = np.random.uniform(0,2*np.pi, size=nsamples),\n", " phase = np.random.uniform(0,2*np.pi, size=nsamples),\n", " geocent_time = 1246527224.169434*np.ones(nsamples),\n", " # spin zero\n", " a_1 = np.random.uniform(0.0,0.8, size=nsamples),\n", " a_2 = np.random.uniform(0.0,0.8, size=nsamples),\n", " tilt_1 = np.random.uniform(0, np.pi, size=nsamples), # tilt angle of the primary black hole in radians\n", " tilt_2 = np.random.uniform(0, np.pi, size=nsamples),\n", " phi_12 = np.random.uniform(0, 2*np.pi, size=nsamples), # Relative angle between the primary and secondary spin of the binary in radians\n", " phi_jl = np.random.uniform(0, 2*np.pi, size=nsamples), # Angle between the total angular momentum and the orbital angular momentum in radians\n", ")\n", "\n", "# ANN model method\n", "ann_pdet = gwsnr_ann.pdet(gw_param_dict=param_dict.copy(), pdet_type='boolean', distribution_type='fixed_snr')\n", "\n", "# inner product method\n", "bilby_pdet = gwsnr_bilby.pdet(gw_param_dict=param_dict.copy(), pdet_type='boolean', distribution_type='fixed_snr')" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Confusion Matrix:\n", "[[6151 209]\n", " [ 126 3514]]\n", "Accuracy: 96.65 %\n" ] } ], "source": [ "# accuracy of the probability of detection\n", "from sklearn.metrics import confusion_matrix, accuracy_score\n", "\n", "pdet_pred = ann_pdet['pdet_net']\n", "pdet_true = bilby_pdet['pdet_net']\n", "\n", "cm = confusion_matrix(pdet_true, pdet_pred)\n", "print(\"Confusion Matrix:\")\n", "print(cm)\n", "\n", "acc = accuracy_score(pdet_true, pdet_pred)\n", "print(f\"Accuracy: {acc*100:.2f} %\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hybrid SNR Recalculation for $P_{\\rm det}$ Estimation" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "psds not given. Choosing bilby's default psds\n", "Intel processor has trouble allocating memory when the data is huge. So, by default for IMRPhenomXPHM, duration_max = 64.0. Otherwise, set to some max value like duration_max = 600.0 (10 mins)\n", "ANN model and scaler path is not given. Using the default path.\n", "ANN model for L1 is loaded from ./ann_data/ann_model_L1.h5.\n", "ANN scaler for L1 is loaded from ./ann_data/scaler_L1.pkl.\n", "ANN error_adjustment for L1 is loaded from ./ann_data/error_adjustment_L1.json.\n", "ANN model for H1 is loaded from ./ann_data/ann_model_H1.h5.\n", "ANN scaler for H1 is loaded from ./ann_data/scaler_H1.pkl.\n", "ANN error_adjustment for H1 is loaded from ./ann_data/error_adjustment_H1.json.\n", "ANN model for V1 is loaded from ./ann_data/ann_model_V1.h5.\n", "ANN scaler for V1 is loaded from ./ann_data/scaler_V1.pkl.\n", "ANN error_adjustment for V1 is loaded from ./ann_data/error_adjustment_V1.json.\n", "Interpolator will be loaded for L1 detector from ./interpolator_pickle/L1/partialSNR_dict_5.pickle\n", "Interpolator will be loaded for H1 detector from ./interpolator_pickle/H1/partialSNR_dict_5.pickle\n", "Interpolator will be loaded for V1 detector from ./interpolator_pickle/V1/partialSNR_dict_5.pickle\n", "\n", "\n", "\n", "Initializing GWSNR class...\n", "\n", "psds not given. Choosing bilby's default psds\n", "Intel processor has trouble allocating memory when the data is huge. So, by default for IMRPhenomXPHM, duration_max = 64.0. Otherwise, set to some max value like duration_max = 600.0 (10 mins)\n", "\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/phurailatpamhemantakumar/anaconda3/envs/gwsnr/lib/python3.10/site-packages/sklearn/base.py:442: InconsistentVersionWarning: Trying to unpickle estimator StandardScaler from version 1.7.0 when using version 1.7.1. This might lead to breaking code or invalid results. Use at your own risk. For more info please refer to:\n", "https://scikit-learn.org/stable/model_persistence.html#security-maintainability-limitations\n", " warnings.warn(\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from datetime import datetime\n", "import gwsnr\n", "\n", "\n", "gwsnr_hybrid = gwsnr.GWSNR(\n", " npool=4,\n", " snr_method='ann',\n", " waveform_approximant='IMRPhenomXPHM', # Use a specific waveform approximant\n", " gwsnr_verbose=False,\n", " snr_recalculation=True, # Enable recalculation of SNRs\n", " snr_recalculation_range=[6,18], # Range of optimal SNR for recalculation\n", " snr_recalculation_waveform_approximant='IMRPhenomXPHM', # Waveform approximant for recalculation. This can be different from the main one\n", ")\n", "\n", "gwsnr_bilby = gwsnr.GWSNR(\n", " npool=4,\n", " snr_method='inner_product', \n", " waveform_approximant='IMRPhenomXPHM', \n", " gwsnr_verbose=False,\n", ")" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Recalculating SNR for 3776 out of 10000 samples in the SNR range of 6 to 18\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████████████████████| 3776/3776 [00:04<00:00, 882.75it/s]\n", "100%|████████████████████████████████████████████████████████| 10000/10000 [00:13<00:00, 718.46it/s]\n" ] } ], "source": [ "# gerneral case, random parameters\n", "# add random seed for reproducibility\n", "np.random.seed(42)\n", "nsamples = 10000\n", "mtot = np.random.uniform(2*4.98, 2*112.5,nsamples)\n", "mass_ratio = np.random.uniform(0.2,1,size=nsamples)\n", "param_dict = dict(\n", " # convert to component masses\n", " mass_1 = mtot / (1 + mass_ratio),\n", " mass_2 = mtot * mass_ratio / (1 + mass_ratio),\n", " # Fix luminosity distance\n", " luminosity_distance = np.random.uniform(100, 10000, size=nsamples), # Random luminosity distance between 100 and 10000 Mpc\n", " # Randomly sample everything else:\n", " theta_jn = np.random.uniform(0,2*np.pi, size=nsamples),\n", " ra = np.random.uniform(0,2*np.pi, size=nsamples), \n", " dec = np.random.uniform(-np.pi/2,np.pi/2, size=nsamples), \n", " psi = np.random.uniform(0,2*np.pi, size=nsamples),\n", " phase = np.random.uniform(0,2*np.pi, size=nsamples),\n", " geocent_time = 1246527224.169434*np.ones(nsamples),\n", " # spin zero\n", " a_1 = np.random.uniform(0.0,0.8, size=nsamples),\n", " a_2 = np.random.uniform(0.0,0.8, size=nsamples),\n", " tilt_1 = np.random.uniform(0, np.pi, size=nsamples), # tilt angle of the primary black hole in radians\n", " tilt_2 = np.random.uniform(0, np.pi, size=nsamples),\n", " phi_12 = np.random.uniform(0, 2*np.pi, size=nsamples), # Relative angle between the primary and secondary spin of the binary in radians\n", " phi_jl = np.random.uniform(0, 2*np.pi, size=nsamples), # Angle between the total angular momentum and the orbital angular momentum in radians\n", ")\n", "\n", "# ANN model method, with SNR recalculation\n", "hybrid_pdet = gwsnr_hybrid.pdet(gw_param_dict=param_dict.copy(), pdet_type='boolean', distribution_type='fixed_snr')\n", "\n", "# inner product method\n", "bilby_pdet = gwsnr_bilby.pdet(gw_param_dict=param_dict.copy(), pdet_type='boolean', distribution_type='fixed_snr')" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Confusion Matrix:\n", "[[6360 0]\n", " [ 1 3639]]\n", "Accuracy: 99.99 %\n" ] } ], "source": [ "# accuracy of the probability of detection\n", "from sklearn.metrics import confusion_matrix, accuracy_score\n", "\n", "pdet_pred = hybrid_pdet['pdet_net']\n", "pdet_true = bilby_pdet['pdet_net']\n", "\n", "cm = confusion_matrix(pdet_true, pdet_pred)\n", "print(\"Confusion Matrix:\")\n", "print(cm)\n", "\n", "acc = accuracy_score(pdet_true, pdet_pred)\n", "print(f\"Accuracy: {acc*100:.2f} %\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## JAX assisted Inner product with `ripplegw` as backend" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "psds not given. Choosing bilby's default psds\n", "\n", "\n" ] } ], "source": [ "import numpy as np\n", "import gwsnr\n", "\n", "# innitialize the class for pdet calculation\n", "gwsnr = gwsnr.GWSNR(\n", " snr_method='inner_product_jax', \n", " waveform_approximant='IMRPhenomXAS', \n", " gwsnr_verbose=False,\n", ")" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/phurailatpamhemantakumar/anaconda3/envs/gwsnr/lib/python3.10/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n", " self.pid = os.fork()\n", "100%|█████████████████████████████████████████████████████████████████| 4/4 [00:00<00:00, 8.52it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "SNR (inner product JAX assisted) : [11.35192063 9.22035654 17.49730142 17.53788394]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "100%|█████████████████████████████████████████████████████████████████| 4/4 [00:00<00:00, 12.69it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "SNR (inner product bilby) : [11.35192174 9.22036535 17.49712218 17.53863404]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# set the GW parameters\n", "mass_1 = np.array([5, 10.,50.,200.])\n", "ratio = np.array([1, 0.8,0.5,0.2])\n", "param_dict = dict(\n", " mass_1 = mass_1,\n", " mass_2 = mass_1*ratio,\n", " luminosity_distance = np.array([1000, 2000, 3000, 4000]),\n", " theta_jn = np.array([0.1, 0.2, 0.3, 0.4]),\n", " ra = np.array([0.1, 0.2, 0.3, 0.4]), \n", " dec = np.array([0.1, 0.2, 0.3, 0.4]), \n", " psi = np.array([0.1, 0.2, 0.3, 0.4]),\n", " a_1 = np.array([0.1, 0.2, 0.3, 0.4]),\n", " a_2 = np.array([0.1, 0.2, 0.3, 0.4]),\n", " geocent_time = np.array([0.0, 0.0, 0.0, 0.0]),\n", " phase = np.array([0.0, 0.0, 0.0, 0.0]),\n", ")\n", "\n", "# jax.jit functions are slow when run for the first time\n", "snr_jax = gwsnr.optimal_snr(gw_param_dict=param_dict)\n", "print(\"SNR (inner product JAX assisted) : \", snr_jax[\"snr_net\"])\n", "\n", "# snr with inner product\n", "snr_bilby = gwsnr.optimal_snr_with_inner_product(gw_param_dict=param_dict)\n", "print(\"SNR (inner product bilby) : \", snr_bilby[\"snr_net\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* as of 27 Oct 2025, jax implementation of inner product is not fully optimized yet.\n", "\n", "* While the waveform generation is fast via `jax.vamp` parallelization, the overall inner product calculation has extra overhead which will be optimized in the future." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Changing the initialization arguments\n", "\n", "To change the initialization arguments, you can modify the `GWSNR` class initialization parameters. Below is an example of how to set up the `GWSNR` class with different parameters:\n", "\n", "What will change in this example? \n", "\n", "* ifos: CE, ET ; interferometers will be changed to Cosmic Explorer (CE) and Einstein Telescope (ET). \n", "\n", "* minimum_frequency: 10 Hz\n", "\n", "* mtot_cut: True ; this will set SNR=0 for total mass > mtot_max, and mtot_max is set according to chirp_time and minimum frequency.\n", "\n", "* Waveform model: TaylorF2\n", "\n", "* multiprocessing_verbose: False ; The progress bar won't be shown but the calculation will be faster.\n", "\n", "* mtot_min: 2*1.0 ; minimum total mass in solar masses Mo. 1.0 Mo is the minimum component mass of BNS systems in GWTC-3\n", "\n", "* snr_method: 'interpolation_no_spins_jax' ; this will use the interpolation method with no spins and JAX acceleration." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Initializing GWSNR class...\n", "\n", "Interpolator will be generated for CE detector at ./interpolator_pickle/CE/partialSNR_dict_1.pickle\n", "Interpolator will be generated for ET1 detector at ./interpolator_pickle/ET1/partialSNR_dict_1.pickle\n", "Interpolator will be generated for ET2 detector at ./interpolator_pickle/ET2/partialSNR_dict_1.pickle\n", "Interpolator will be generated for ET3 detector at ./interpolator_pickle/ET3/partialSNR_dict_1.pickle\n", "Please be patient while the interpolator is generated\n", "Generating interpolator for ['CE', 'ET1', 'ET2', 'ET3'] detectors\n", "\n", "Saving Partial-SNR for CE detector with shape (20, 200)\n", "\n", "Saving Partial-SNR for ET1 detector with shape (20, 200)\n", "\n", "Saving Partial-SNR for ET2 detector with shape (20, 200)\n", "\n", "Saving Partial-SNR for ET3 detector with shape (20, 200)\n", "\n", "Chosen GWSNR initialization parameters:\n", "\n", "npool: 8\n", "snr type: interpolation_no_spins_jax\n", "waveform approximant: TaylorF2\n", "sampling frequency: 2048.0\n", "minimum frequency (fmin): 10.0\n", "reference frequency (f_ref): 10.0\n", "mtot=mass1+mass2\n", "min(mtot): 2\n", "max(mtot) (with the given fmin=10.0): 235.0\n", "detectors: ['CE', 'ET1', 'ET2', 'ET3']\n", "psds: [[array([ 5. , 5.01153 , 5.0230867, ...,\n", " 4977.0194 , 4988.4965 , 5000. ]), array([1.36418639e-44, 1.28941521e-44, 1.21896416e-44, ...,\n", " 3.43993064e-48, 3.46829426e-48, 3.49694151e-48]), ], [array([1.0000000e+00, 1.0030759e+00, 1.0061612e+00, ...,\n", " 9.9387655e+03, 9.9693357e+03, 1.0000000e+04]), array([1.22216783e-33, 1.25637893e-33, 1.28147776e-33, ...,\n", " 2.90395567e-47, 2.92250209e-47, 2.94117631e-47]), ], [array([1.0000000e+00, 1.0030759e+00, 1.0061612e+00, ...,\n", " 9.9387655e+03, 9.9693357e+03, 1.0000000e+04]), array([1.22216783e-33, 1.25637893e-33, 1.28147776e-33, ...,\n", " 2.90395567e-47, 2.92250209e-47, 2.94117631e-47]), ], [array([1.0000000e+00, 1.0030759e+00, 1.0061612e+00, ...,\n", " 9.9387655e+03, 9.9693357e+03, 1.0000000e+04]), array([1.22216783e-33, 1.25637893e-33, 1.28147776e-33, ...,\n", " 2.90395567e-47, 2.92250209e-47, 2.94117631e-47]), ]]\n", "\n", "\n" ] } ], "source": [ "# if snr_method = 'inner_product', interpolator will not be created\n", "import gwsnr\n", "gwsnr = gwsnr.GWSNR(\n", " npool = int(8), \n", " minimum_frequency = 10.,\n", " mtot_min = 2*1, \n", " mtot_cut=True,\n", " waveform_approximant='TaylorF2',\n", " snr_method = 'interpolation_no_spins_jax',\n", " #psds = {'CE':'CE_psd.txt', 'ET':'ET_B_psd.txt'}, # if you want to use your own psd\n", " ifos = ['CE', 'ET'], # this will considet bilby's default psd of CE and ET\n", " multiprocessing_verbose=False,\n", ")" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'CE': array([ 22.20392799, 36.01436234, 112.26751709, 131.33247375]),\n", " 'ET1': array([1.57145751, 2.54927254, 7.9575057 , 9.31898499]),\n", " 'ET2': array([ 4.48666096, 7.27841663, 22.71943855, 26.60659027]),\n", " 'ET3': array([ 4.99195719, 8.09812546, 25.27814484, 29.60307503]),\n", " 'snr_net': array([ 23.2493782 , 37.71058273, 117.56902313, 137.54750061])}" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mass_1 = np.array([5, 10.,50.,100.])\n", "ratio = np.array([1, 0.8,0.5,0.2])\n", "dl = 10000 * np.ones_like(mass_1)\n", "gwsnr.optimal_snr(mass_1=mass_1, mass_2=mass_1*ratio, luminosity_distance=dl)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "gwsnr2", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.19" } }, "nbformat": 4, "nbformat_minor": 2 }