diff --git a/gen_samples_py3/.ipynb_checkpoints/99-whiten_wavepic_draw-checkpoint.ipynb b/gen_samples_py3/.ipynb_checkpoints/99-whiten_wavepic_draw-checkpoint.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..975ed1e1ecda9e5770312ef1bd8bc5c1c739eb12 --- /dev/null +++ b/gen_samples_py3/.ipynb_checkpoints/99-whiten_wavepic_draw-checkpoint.ipynb @@ -0,0 +1,347 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "import time\n", + "import numpy as np\n", + "import pandas as pd\n", + "import utils.samplefiles\n", + "\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "# Start the stopwatch\n", + "script_start_time = time.time()" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "hdf_file_path = './output/train.hdf'\n", + "plot_path = './plots/signalnoise_id0.png'\n", + "# which sample is used to shown\n", + "sample_id = 0" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "key: config_file values: default.json\n" + ] + } + ], + "source": [ + "data = utils.samplefiles.SampleFile()\n", + "data.read_hdf(hdf_file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "# when split_injections_noise is set to be True, GWs contain signals and pure waves are splited\n", + "df, noise = data.as_dataframe(injection_parameters=True, \n", + " static_arguments=True, \n", + " command_line_arguments=False, \n", + " split_injections_noise=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Index(['approximant', 'bandpass_lower', 'bandpass_upper', 'coa_phase', 'dec',\n", + " 'delta_f', 'delta_t', 'distance', 'domain', 'event_time', 'f_lower',\n", + " 'fd_length', 'h1_output_signal', 'h1_signal', 'h1_snr', 'h1_strain',\n", + " 'inclination', 'injection_snr', 'l1_output_signal', 'l1_signal',\n", + " 'l1_snr', 'l1_strain', 'mass1', 'mass2', 'noise_interval_width',\n", + " 'original_sampling_rate', 'polarization', 'ra', 'sample_length',\n", + " 'scale_factor', 'seconds_after_event', 'seconds_before_event', 'spin1z',\n", + " 'spin2z', 'target_sampling_rate', 'td_length', 'tukey_alpha',\n", + " 'waveform_length', 'whitening_max_filter_duration',\n", + " 'whitening_segment_duration'],\n", + " dtype='object')\n" + ] + } + ], + "source": [ + "print(df.columns)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "7.340279606636548\n" + ] + } + ], + "source": [ + "print(df.injection_snr[sample_id])" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "36.217808319315374\n", + "76.55000144869413\n", + "0.7305299539277823\n", + "0.5974611672286425\n", + "0.017164934231118905\n" + ] + } + ], + "source": [ + "print(df.mass1[sample_id])\n", + "print(df.mass2[sample_id])\n", + "print(df.spin1z[sample_id])\n", + "print(df.spin2z[sample_id])\n", + "print(df.scale_factor[sample_id])" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.9801424781769557\n", + "0.48680334688549004\n", + "3.776917009710014\n", + "0.821770111935331\n", + "4.448951217224888\n" + ] + } + ], + "source": [ + "print(df.coa_phase[sample_id])\n", + "print(df.inclination[sample_id])\n", + "print(df.ra[sample_id])\n", + "print(df.dec[sample_id])\n", + "print(df.polarization[sample_id])" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "sample = df.loc[sample_id]" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.5\n", + "2048.0\n", + "8.0\n" + ] + } + ], + "source": [ + "# Read out and construct some necessary values for plotting\n", + "seconds_before_event = float(sample['seconds_before_event'])\n", + "seconds_after_event = float(sample['seconds_after_event'])\n", + "target_sampling_rate = float(sample['target_sampling_rate'])\n", + "sample_length = float(sample['sample_length'])\n", + "print(seconds_before_event)\n", + "print(target_sampling_rate)\n", + "print(sample_length)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a grid on which the sample can be plotted so that the\n", + "# event_time is at position 0\n", + "grid = np.linspace(0 - seconds_before_event, 0 + seconds_after_event, int(target_sampling_rate * sample_length))\n", + "\n", + "# for time from -0.15s to 0.05s\n", + "#grid = np.linspace(0 - seconds_before_event, 0 + seconds_after_event, int(target_sampling_rate * sample_length)+1)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "det_name = 'H1'\n", + "det_string = 'h1_strain'" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "145.16235\n" + ] + } + ], + "source": [ + "maximum = np.max(sample[det_string])\n", + "print(maximum)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.173865814288047e-21\n" + ] + } + ], + "source": [ + "maximum = max(np.max(sample['h1_signal']), np.max(sample['l1_signal']))\n", + "print(maximum)" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "23.74179629063074\n" + ] + } + ], + "source": [ + "maximum = max(np.max(sample['h1_output_signal']), np.max(sample['l1_output_signal']))\n", + "print(maximum)" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAF9CAYAAACUKdl5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAB7QklEQVR4nO3dd5gT5fbA8e/ZwtJ771XpTYqIIogoil2vvWDv9V6Va+8Xe9d70fsTr4rYGxaKiiigAoIiIL23pfft7++PJEs2mzJJZjLJ5nyeJ89uksnkpM2ceed9zyvGGJRSSimlVPrIcDsApZRSSimVWJoAKqWUUkqlGU0AlVJKKaXSjCaASimllFJpRhNApZRSSqk0owmgUkoppVSayXI7gGRWv35907p1a7fDUEoppZSKyZw5c7YaYxoE3q4JYBitW7dm9uzZboehlFJKKRUTEVkd7HY9BayUUkoplWY0AVRKKaWUSjOaACqllFJKpZmkTgBFZJCIfC4i60XEiMjIgPvHem/3v/wcsEyOiLwoIltFZJ93fc0T+kKUUkoppZJIUieAQHXgT+Bm4ECIZaYATfwuJwbc/xxwJnAecBRQE5ggIpkOxKuUUkoplfSSehSwMeYr4CvwtPaFWCzfGLMp2B0iUgu4HLjUGDPZe9tFwGrgWGCi3TErpZRS6WD37t3k5uZSWFjodihpKTs7m4YNG1KzZs2YHp/UCaBFR4pILrAT+AG42xiT673vMCAbmORb2BizVkQWAUegCaBSSikVtd27d7N582aaNWtGlSpVEBG3Q0orxhgOHDjA+vXrAWJKApP9FHAk3wAXA0OBvwP9gO9EJMd7f2OgGNga8LjN3vvKEZGrRGS2iMzesmWLM1ErpZRSKSw3N5dmzZpRtWpVTf5cICJUrVqVZs2akZubG/kBQaR0C6AxZrzf1fkiMgfP6d0RwMcxrnMMMAagT58+Ju4glVJKqQqmsLCQKlWquB1G2qtSpUrMp+BTvQWwDGPMBmAd0MF70yYgE6gfsGgj731KKaWUioG2/Lkvns+gQiWAIlIfaAZs9N40BygEhvkt0xzoBMxIeIBKKaWUUkkgqU8Bi0h1oL33agbQUkR6Atu9lweAj/AkfK2BfwG5wCcAxphdIvJf4AnvQJFtwDPAH3jKxyilVFoxxrBy6z7aNqjudihKKRclewtgH2Cu91IFeND7/0N4Bnd0Az4DlgBvAouBAcaYPX7ruAVPQvgeMB3YC5xsjClOzEtQSqnkMXbGKo55+gd+W7PD7VCUSmpjxozh008/tXWdY8eORUTYu3evreuNRVK3ABpjpgLhTnAfb2Ed+cCN3otSSqW1eWt3ArBm2356t6zjbjBKJbExY8bQtWtXTjvtNNvWOWLECGbOnEnVqlVtW2eskjoBVEoppZRKZgcOHLA8IrpBgwY0aNDA4YisSfZTwEoppZRSjliwYAHDhw+nbt26VKtWjU6dOvHyyy8zePBg5syZw5tvvomIICKMHTsWgNatW/P3v/+dhx9+mObNm5cWYZ45cyannHIKTZo0oVq1avTs2ZN33nmnzPMFngJetWoVIsL777/P1VdfTa1atWjevDn3338/JSUljr52bQFUSimlVFo6+eST6dSpE2+//TY5OTksXryY3bt388orr3DmmWfStm1b7r33XgDatWtX+rhx48bRpUsXXnnlFYqKigBYvXo1AwcO5JprrqFy5cpMnz6dSy+9lIyMDM4777ywcdxxxx2ceeaZfPjhh3z77bc89NBDdOnShbPPPtux164JoFJKKaXi9uAXC1i4Ybcrz925aU3uP7lLVI/ZunUrK1eu5LPPPqNbt24ADB06tPT+atWq0aBBAw4//PCgj58wYQKVK1cuvX7uueeW/m+MYdCgQaxbt47XXnstYgI4aNAgnn76aQCGDRvGN998w8cff+xoAqingJVSSimVdurWrUuLFi245ppreO+996KaUm3o0KFlkj+AHTt2cNNNN9GqVSuys7PJzs5mzJgxLFmyJOL6jjvuuDLXO3fuzLp16yzHEwttAVRKKaVU3KJtgXNbRkYGkyZN4u677+ayyy7jwIEDDBw4kBdeeIFevXqFfWyjRo3K3TZy5Eh+/vln7r33Xjp37kzNmjV59dVX+eyzzyLGUrt27TLXK1WqRF5eXlSvJ1qaACqllFIqLXXs2JGPPvqIwsJCfvzxR+68805GjBgRsfUtcAq2vLw8JkyYwMsvv8w111xTervTAznioaeAlVJKKZXWsrOzOeaYY7jtttvYuHEjO3fujKoVLj8/n5KSEnJyckpv27NnD59//rlTIcdNWwCVUkoplXb++OMP/vGPf3DOOefQtm1bduzYweOPP06PHj2oW7cuHTt2ZOLEiUycOJF69erRpk0b6tWrF3RdtWrVom/fvjz00EPUrFmTjIwMRo8eTa1atdi9252BMZFoC6BSSiml0k7jxo1p1KgRjz76KCeccALXXXcdnTp1Km21u+eee+jUqRNnn302ffv25Ysvvgi7vnHjxtG2bVsuvvhibr75Zs4880wuvvjiRLyUmIgxxu0YklafPn3M7Nmz3Q5DKaVsc/P4uXw2bwPPndOT03o1czsclaIWLVpEp06d3A5DEfmzEJE5xpg+gbdrC6BSSimlVJrRBFAppZRSKs1oAqiUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaUYTQKWUUkqpNKMJoFJKKaVUmtEEUCml0ohO/qSUNYMHD+ass85yNQYR4aWXXnJk3VmOrFUppVRSE3E7AqWS2yuvvEJ2drbbYThGE0CllFJKqQCdO3d2OwRH6SlgpZRSSqWlBQsWMHz4cOrWrUu1atXo1KkTL7/8MhD8FPAHH3xAhw4dqFKlCkOGDGHu3LmICGPHji1dpnXr1vzjH//g2WefpXnz5tSpU4dzzz2XnTt3li6zb98+brjhBg499FCqVq1KmzZtuP7669m9e3ciXjagLYBKKaWUSlMnn3wynTp14u233yYnJ4fFixeHTMJmz57Nueeey1lnncWLL77IokWLOOecc4Iu+/7779O9e3fGjBnDunXruO2227jrrrt45ZVXANi/fz/FxcU8+uijNGjQgLVr1/Loo4/yt7/9jYkTJzr2ev1pAqiUUilk5/4CRn/9Fw+c0oXK2Zkxr0cHgyjbfT0KNs1357kbd4MTRkf1kK1bt7Jy5Uo+++wzunXrBsDQoUNDLv/444/TqVMnxo8fj4gwfPhwCgsLufPOO8stm52dzaeffkpWlifNWrhwIePHjy9NABs0aMCrr75aunxRURFt2rThyCOPZM2aNbRs2TKq1xILPQWslFIp5KlJixk/ay0fzFkX0+N18IdSHnXr1qVFixZcc801vPfee+Tm5oZdftasWZx88smI34/olFNOCbrskCFDSpM/8PQnzM3NpbCwsPS2t956i169elG9enWys7M58sgjAViyZEk8L8sybQFUSqkUoi13KmlF2QLntoyMDCZNmsTdd9/NZZddxoEDBxg4cCAvvPACvXr1Krf8pk2baNCgQZnbAq/71K5du8z1SpUqYYwhPz+f7OxsPvnkEy6++GKuvfZaHnvsMerWrcvGjRs5/fTTycvLs+01hhNVC6CIVBKRZiLSTkTqOBWUUkoppZTTOnbsyEcffcTOnTuZMmUKeXl5jBgxgpKSknLLNm7cmC1btpS5LfC6VR988AH9+/fnlVde4YQTTqB///7UqZPYtCpiAigiXUTkcRGZA+wF1gBLgK0ikisin4rIhSJSxelglVJKKaXslp2dzTHHHMNtt93Gxo0by4zY9enbty9ffPEFxq8Z/vPPP4/p+Q4cOEBOTk6Z2955552Y1hWrkKeARWQg8AgwCJgF/AC8AGwF8oHaQGugD/As8KKIPAM8a4zZ62jUSqmUYIzh1R+Wc1bv5jSsWdntcJRSqtQff/zBP/7xD8455xzatm3Ljh07ePzxx+nRowd169Ytt/ydd95J//79Offcc7n00ktZtGgRr732GuA5nRyNYcOGcf311/Poo4/Sv39/vvrqK7799ltbXpdV4foAfown4bvIGBO2t7GIZALHArd4b3rYluiUUilt0cY9PPHNYqYu3sL7Vw9wOxyllCrVuHFjGjVqxKOPPsqGDRuoXbs2Q4YM4fHHHw+6fJ8+fXj33Xe56667+Oyzz+jTpw+vvvoqw4YNo2bNmlE999VXX82KFSt4/vnnycvLY9iwYYwbN47DDz/cjpdmSbgEsJUxxlJPRGNMMTARmCgiepivlAKguMRzqmR/QZHLkSilVFkNGzbkrbfeCnn/1KlTy9129tlnc/bZZ5def/vttwHo0aNH6W2rVq0q97iRI0cycuTI0uuZmZk89dRTPPXUU2WWMwGjvAKv2ylkAmg1+bPrcUoppZRSyezaa69l2LBh1KlTh99++41HHnmEESNG0KZNG7dDi1rcZWBEpB7QxRgzzYZ4lFJKKaWS0rZt27juuuvYtm0b9erV45xzzuGJJ55wO6yY2FEHcDDwPhB7SXqllFJKqST3/vvvux2CbXQmEKWUYwxatVgppZJRuDIwKyyuo6pNsShgzbb91KteiWo5OkmLqjgEnX9MKaWSSbgsoxkwF/gxwjoOAU6yLaI0N+jJ7zmsVR0+uvYIx57DGMOTExdzYrcmdG1Wy7HnUUopVXEZY8rMi6sSL55RwuESwHnAJmPM7eFWICJnogmgreas3uHo+otKDK9MXc5rP65g6aMnOvpcKnkt37KXoU//wJTbjqZ9w+puh6OUSiHZ2dkcOHCAqlX1JKCbDhw4QHZ2dkyPDdcH8Begv8X16CFACtJJ5dPbF79vAOBz71+llLKqYcOGrF+/nv379ztaq04FZ4xh//79rF+/noYNG8a0jnAtgA8B/7UQxEfoYBKllFIqbfhmvtiwYQOFhYUuR5OesrOzadSoUdSzkPiEKwS9Fc+8v0opFRNtGLBfRXhL84uKeWriYm4a2oEalWM7faXcV7NmzZiTD+U+bblTSjlO+4knj2RIyt+btZbXflzJC98udTsUpdJWuDIw0VQ7NMaYc2yIRymVaMmQESjL7Mql3UzKi4o937nCYv3uxeLm8XPp0LA6NxzTwe1QVAoL1wewQcB1AY7CUxpmj2MRKaWSzupt+1ixdR9DDo2ts7FSyj6fzfMM3NIEUMUj5ClgY8wQ/wtwLJ4k8KrA+7z3O0JEBonI5yKyXkSMiIwMuF9E5AER2SAiB0Rkqoh0CVimjoi8JSK7vJe3RKS2UzGrimfxpj1JPdJt6uJc/ly/K7YHW2gKOvrJqVz6xqzY1u+gtdv3U1ySvJ+LUkolq2j6ALq1la0O/AncDBwIcv8dwN+BG4G+QC4wWURq+C0zDugNDPdeegNvORhzUnjsq0UMf26a22G4btvefMb/uibmx/+8YhvHPzeNt35ebWNU9hr5xixOevEnt8MIadOuPF6ZuszWJHr9zgMc9cT3PDHxL9vWaafnpyyl9agvKSoucTuUkEZ99AetR33pdhhKKRck/SAQY8xXxpi7jDEfAmW2pOIpQX4LMNoY85Ex5k/gEqAGcL53mU54kr6rjDEzjTEzgauBk0Tk0AS+lIQbM20Ff20qf7Y+iRuyHHHDuLmM+ng+K7fui+nxq7d5HhdzC1sa833Vcvfk88Q3i1m+Za9t6966Jx+Amcu32bZOO/37h+UAFCRxAjh+1lpXn3/xpj22fidUeM9PWcpbM1e5HUYZBwqKufJ/s1m3Y7/boaSdpE8AI2gDNAYm+W4wxhwApgG+udQGAHuBGX6Pmw7s81smLYU683f9O79VqFaBrXs9iUJhEu+IXZXAIwI9W6v8zVyxjaFP/+B2GGnj2SlLuPezBW6HUcakhZuYvHAzj3+zOOZ17Msv4oHPF3CgoNjGyOyTuyePl75bmnTdiGJJAJPpFTT2/t0ccPtmv/saA1uM3zvv/T/Xb5lSInKViMwWkdlbtmxxIOTk9+X8jW6HYKtk+sLabXdeIXvzi9wOQ6moaFkgZad//7CcsTNW8WaStW763Pbe7zw1aQm/r0uus0jhysDMIvi+8y0RKddWa4zpZ2dgbjHGjAHGAPTp06ci5w5pw5f7V8R9TvcHJpGZEecr072xI0yFPvSIT5I1hKgUV+Q9tZCsA8L2FXgO0pMtvnAtgAuCXN4EZoW4zw2bvH8bBdzeyO++TUADb39BoLTvYEO/ZVQaqKh5TrJtVFRZYuHQ469Nu5mxLL6Jl/KLirl5/NxyfalWb9tHXmH5U2PpnoQdKCjmvVlrEnZabs22/XxVwc6uqNQWbiq4kQmMI1Yr8SRxw/AkpohIZTz1Cm/3LjMTz0jiARzsBzgAqEbZfoHKAmMMizfvoWPj1Jn+J833cyoFDH/uRwBWjR4RdrmCopKQ/Zx+XLKVz+ZtYG9eEf8d2Rfw/F6PfnIqRx/SgDcv85ykqagHQtH619eL+N/M1TSsWTkh9S1PeH4a+wqKI37G6SrZ+selg6QfBCIi1UWkp4j0xBNvS+/1lt6+fM8Bd4rIGSLSFRiLZ9DHOABjzCLgG+A/IjJARAYA/wEmGGNi73Wapt6Yvorhz/3Iryu3ux1KDHTPF5RueFPGuWNm8vHc9VE/7oclFbc/c0mJ4ZWpy9h1oNDyY777azMbdnqqiu3PT8zAgX1JOkDBbaJHJK4JmQCKyL0iUiualYnIMSJycvxhldEHz+wjc4EqwIPe/x/y3v8E8CzwMjAbaAIcZ4zxr39yPvA7MNF7+R24yOY408KfGzydWNdsT6Eh+5rfWHLyiz9xykv21RJ8bdoKFmxIrk7PieREXv3bmp32rzTFTVu6hSe+Wcz9n/1pafm5a3Zw2djZTFmU63BkqWdPXiH5RZqopotwU8H1BdaKyGfAB8BMY0yZw0gRyQa6AScA5+CZPu4SOwM0xkwlTNONtxXwAe8l1DI7gAvtjEslTmFxCYc/9i0PntqFk7o3jXk9eqAZgveNmW9zncNHv1pk6/pSlX7vyrPzPSko8pR32muxJW9nFC2F6abbA5Po2LgG39wyyO1QKqjkao0INxXcKXj61gnwLrBJRDaLyAIR+U1EVuCZE3gWnuTv/4B2xphJodap7LM7r5BvFwVWv7Em1UYn7thXwLZ9BTz4xcKYHp9ar1ZZpZ+rSmXFJYbXf1yRdLXrgk0eEOijOeuiOuVuRUX+PSfrMWDYPoDGmF+MMRfiGVV7Ip7+dlPwJH3vAdcAHY0x3Y0xzxljUui8YGq7Zfw8Ln9zNut3BpsdzxoroxMrEqde7czl29i8O8+htSfejOVbeX7K0rDL/Lh0C61Hfcmy3Mg7C6el17c4NXw0Zx25Feg34YQJf2zgkS8X8czk1OqKvnDDbv7+we/c8eHvtqzP6u/3s3nr+W3NDsBTWHl/wcH6p7oNiE24U8CljDF7Odh/TiWBVd5pzYKVd1BlOT267LzXfqZetUrMuXeYo8/jmID35/zXfgHg5mM7hHzIl394ylnMWrWD9g1rhFwuXjv3F7A0dy99W9d17DmUvbbtzefvH/xO5yY1+ermo9wOp7wkaWra723525OXWoXcD3j3ObneqRgT5ebx8wDPSPl+j36rp6ptYCkBVBVTqp0KjpeTo8227StwbN3p7ML//sKf63ez4rETyQhR8DpZv8VOxxXtt/meT+eTkYAOib66lFv2hk4QnDkmi22l6bYdrCisnKpW4WkCmIbS7dSvbt4jSOJRCn+u3w0EDzF5o06MaL/Xb/+8BoBTe8Y+kCoZaRkRd9iexOuGOuGSvg6gSjLeH+l/fljubhxR8G2odDeh7LJzfwE792urb6w0Z0td/p/dlj35XPx/v8b1W7DjuxApdywsLmH8r2socXnWpGQruaoJYIpz6wu1NHevO08cB93pKLv0fGgyPR+a7HYYZbi5b9mTVxiyP3JJiSkt1aIqltd+XMG0JVt4b9Zat0MJa8y0FYz6eD4fzEnuOBNNE8BUZXMyY4zh6UmLWb7F/sSupMTwr68WlVbeT7RU7+Nz47tzmbrYwaK1yXZY6rC12/c7+376uPS2unGg0+2BSQx/blrQ+x6asJBD7vla56xOIa1Hfcmu/dGXeVmyeU/Ug+7s2Pz8sDj8TDc7vH20dx9wd8BNsjVCaAKoAGjzz6948btlXPT6L7asb9XWffR5ZDLrdx5g7tqd/GfaCm55b54t645VIvs+Tl2cy9c2Tfz+xe8bGPnGLFvW5c+NvqDJkGse8/RUR97PUNzc6Ad7v+P9DEpKDON/XUNhcdlWvVXbglcBe/vn1Z7HJcOH74INOw/QetSXUT3mnV9W03rUl67OyrF4c/hBFoGf5g9LtnDcs9P4YM4654IKYsPOAyzcuDuhz1lRhJsKbouI5Fq9JDLoiqS4xNhSymXumh0c9cR37MmL/qjt+GcPHrkX2nSU/u6sNWzdW8DA0d+VHvkXlxg+nLOOi//vV1uewyo39jsj35jFte/8lvgnjkI8LaOpvC8vLE5s8Ke8OD2hz2dVrInpR7+tY9TH8+n36BRe/DZ8vchEsvqdTNQZAWMM89buZMbybUHv35fvaY0aP2strUd9SZFfQv3s5CWA+y1WwYT62iz3dgtauCG6ZCzeA6QDUew/U/1skN3CjQJ+GR2X47jr3pnDxAWbWTV6RFzreXrSEtZuP8C8tTs5qkODsMsG/ggiHenFa+XWg6eV//GBPcVDY5Fsze9OmLhgE2//vJq3Lu9v/UFxvDGJekuNKR+mkxsn3867V8s6ca/L6d9XovlmgNixv5CnJy/hxqGh60UmgpXv4IadB2hau4rjsfh7f/Za7vxoPsd1bhT0/ke+LDtVYkFxCVmZyXNSzun6qYmwJ6+QM16ZQcOaOW6HkpRCJoDGmAcSGEfamrggtuncDor9R2r1FODO/QXUqJxNZog6bFafZ87qHaW3GWMSVr4h0nZswYZdjHjhJz6/YSDdm9dOSEyxaD3qS5Y+egLZYXYSV781J+w6pi/biggc0a6+LTHZuYv4aelW3vp5Ff++8LCovhvBlvxp6Vaa1alCm/rVYorlfzNXc//nCxh7aV8GH9owpnW41dpQAfbbtvp07npueW8e7111OP3b1nP8+YwxPDdlKau2eYr1r9leMSfIWrJpDz2b1wKSt6Xo5xXbWZq7N6ZBi3+s28kpL03nxzuG0KJuVQeic1/yHG6oqDiVOgWud3deIT0fmsy/vloU1RP7J5fBdoT//mFFjBHa79tFnh4MkxfGm4zb7+EJZec/3hvHrAG79hdyweu/cP5rv5Q53RRLxuBE7j7yjV+ZuGAzRTZ0Q7jwv78w5KmpMT9+ibfVbu2O8AOX3GwluffTP8tcT4MG7pj4Djx9LbFO9339bc1Onv92KZ/N2+Do8zgt1Dfbd3B2oLCYN2eujus5/tq0m6Wb93qfL7nSyPHekc0/LAk/wCSUa96aww3jynYDSraDM8sJoIgMEJHXRWSaiPwaeHEyyHTw4rdLk3KU3G7v6Z6v/9zkucGmED+bt96eFVVw//1ppW3r8u8rY0i/guAqNfy8YltU88xGu0lyOtGIdTuejr/H4c/9yPMO9iO1eoBWVFxi+8HcNws2McE7ZWayFiu3lACKyDBgGtAcOBLYAuwFegD1gD9DP1pZ8fTkJUz44+AR49TFubQe9WXpnL+pxs3v+9fzN/L6j+VbGJP0N+iaZDvitoseXKS2c8f8zPuz1/HXpoODCT6du55bI1QRqIi/71T6jfp/Xm6JJYnbtjef9nd/zf9NX2X5MRt2HmDMtNSZDCEUqy2ADwHPA76RCvcaY44BDgEKgan2h5Z+8v2KpfpOH/y2ZkfQZZ3aLCTjRjTa13rtO7+V62Addv1RPsHm3Xmc9vJ0tiR4MnTHJOOHHgc7W03jkWynexIhr7CYKYvsKQrh393hlvfm8cnc8Il9xXq/7f9NRkqgA0Xzfm7alcf7s62Vfykqdq8oeLDXtGFnHgCfzI0c//Z9BWzfV8AVb87msa/+Yq2F/p3jf10TdZyJYjUB7Ax8DZTg2R9XAzDGrAYeAO52IjgV2WVjZ/PolwsjLxgjOzaqdm2Yg20SpyzczFKbRlla3eS+OWMV89bu5P3ZsVeV/+6vzbQe9SUrQhTezi8qZtwv9m44kqE1wWquaYwJefATjdajvuT9JJ+lwC7uf7rw2FeLuOuT+Ql9TicOX1Zt3cdF//2F/QWJK8PS6b5vGP31X1E95tMISbG/SAm0VcHeb6tTwc1cvo32d3/N7FXby93nRAJv9yp7PzyZ3g9PZq+3hI+V2pajPk7s7yEaVhPAPCDDeNpXNwLt/O7bjefUsIrTHR/+Ufp/pKZs349wzfb9vPZjdC0e/s8Tbw3Cmcu38cb02Fpc9uYXRfX8uUFa3K7432yGPRt8BgKfZCxn8MXvnr4h89buDHr/81OWJm5HGsP74/Rb+vnvGzjjlRnWYolw/2tBugOE4tuw2yEw2d25v4A1IYol28nNBt3VMb6+1dv2MWd1+aTACYHf3WDf5ce+WsSPS7cybcnWGNYf+4/j31HOsX7Le/PYvDuv9HpxieHP9btifn5/oQ4Y4/np/7TMM6Di5xXBayOC5/2Lph9lEm7ey/GvgJFMrCaAvwOHev//FviniAwTkaPxnB5O3hQ3xTmxMfcfnfbUxMVxreu8137mwS9ia4Hsev9Ejnz8u7ie3wrf9iFZO+IGs31f7JOrWzFj+baInc5Pe3k63/xZdjaTTbvy+P6vg6f4nHhHjYGVAX1f2971VbmBAaGeO9aY5q/bRdf7J/LlH36v2ca9y7HP/MCgJ7+3bX2JNm/tTlZu3RfzW+L/uMDP6Ognp3LmqzNjji0eS3P3Mmf1jqCJkzGG9S5NYWnVWzNXs9s7AcDzU5Zw0os/2ZYEWmXntvXm8fNod9dXrsdhRbS/hWTLVa0mgM9xMPa7gH3AROB7oCFwve2RqTICd4h22exQPzarP8Ote8MnOpt25YW9PxqhYrJ6arTIO4OEb2kr/T+enbyE3D32vQbwlCi58d25cfWleXbykoive97andwwbm6Z205/ZTqXjnVmGrVI2+7APkZ2b0z/3ODZaf64dEsUp6qt3VdYXBLxuw5woMCeqb82785jd15hTO9R7u48+jwymWW5ZbtWnPby9JhK6yTbcVdgPC98u5QzX53BSS/+xIxlZVv8XvxuGQNHf8fqbaG3v3+u30VJmBYru1uoAstVvfT9Mv7pPc0435v4zVi+1bFBhE5/nJ//bl/5HFu6MAW5Ldbv9F+bkqsovKUE0BjzlTHmZe//64HD8LQI9gTaG2PCV59VUfP/0n01fyNDnprKFLvr1CXZhjnQp3PXc/i/vmX2qoPN54c9PLl0JgKrLG8EIvyqPw7oQzPeQt+y579dyj8+OHjK/dhnfgg7eiyvsJi7P5kf9jXePH4eX/y+wb4ZJqLYmm0MSMiT4Yi2XPRxZhxOnFJ61OKgpIKi2JL6XQcKy/RH6v/Ytwx9+oewj1m1NfgBzDcLNrF1bwFvzoivxluyiKbv69odZd8T3xyzgd97nzmrd3DSiz/xaphTt3H1vQ3y0Cv/N7vcbdv2lj2Q//g3z7bqyxDzkXe9fyKjv/YMYvBv0Y/03In04Zx1UZ9OfzLOM1r+otmKjJ2+0lKDQGDtTrfFVAjaeCw1xvxhjHH2XFWaE4QF3paJRAyzj3Ra0KZcypJZ3o7Ci/1e97Z9BfyxbmdM64tmAIJ/HxSr26Cte/N54PMFFAa0zPn3c1yWu5fHvgrd0fuj39bxzi9rDtZdtFGsiU2wvlmxfr6LbJi0fX9BUVQ7Bsvf2SC3bYtwKn5WkM7swcwN0dfTDtv3FdDjwUk8P6VsPbVIo9SnL4++f5tdnMgtIq3T9/nameBv8J4eXhjme71kc/SzUDjdaro3v4h//7CcoU//4FiLfjC+9/6pSUtC3ufzjw9+5xsHtoN2MngmS3jgi4WcO+Znt8OJWjSFoJuKyFUi8pCIPBFwedzJINPR9GWeTrKhjh4D+zr8tMy9jXmyi3QEHrjh+b/pq2h311fljqojeeDzBYydsSquGUUi9X12ZMcQYY9oZ9+sm8fPi+vx63bsp/N9E3nr59CtU/G+Rf7fl+emhC9Se04SbPTv8p7+C5eEBOVAFhbv9zOaUa2Bz5lfVBz3oLYnJy6h9agv3W78SpiCMN1IQs8E4kwswUT9nS7H+if55/rdnPrydEvL+r8FxvsW7smL7sxUMrBaCPp0YAXwMnA58LcgF2WjrVEmH7FyaoRsNFXt/UexWWEMTFywiVemLrM1Jt+9H83x9DXznfbx3+DNXbMjZL404Y/gp1vsFMvHtXl3Hh8ElKwRsfYZJduO0DfK1ImWAd/nbNdPIpbVxHK68JsF0b0Xi+Psh+RkOaH5foMXok00pi/bRvcHJpVeLykx5VpBv/lzE5e/Wf4Uqo9vuxv4HUiG3jIzl2/j49+C16oThFVb98U8CtsJizbu5s0Zq0qv780vivqb8+J30W3jA1n5LRcUHzxo+N2vpX6bhT67ZZ4rqqWTQ5bF5R4DJgEjjTGJGauvgMRPD2TX0Z3/TiLSD+OeT//ktYv7WF73Nws2RVUjz+oO3cpiX/6xkaxMd6fQjjZpv+T/fuWvTXuYcOORpbc58a3adaCQG8b9xlN/60GjmpVjXo/BRJWEbbdYgyySVJyK64HPF5S7LVIrsm19R+Ng5Z2OJRH3tWjl7s7j/dlreWrSEgYf2qD0/icnRldnrzQWK8sYw6xVOxxLBM57LXxr8+AYBujYJdg26YTnfwTgkiNak7s7j36PfUtGmA8+1n1PvPusUIXjrRxYJWOJsWhY3ZO1AF7Q5M9dvu9aQVGJpQKUkVcY38NfsGkOx0i/38Awt8Y4cjnUhiLZRilaFSlhySsspqTElLZqBNbW8iXpL8RwlB3s6/fRnHX8uHQrr061VsvsuSlLePn7ZUxZuJkZcfRHW7vd3jIdqbRJH+vXwuLz/WJ7ZuJIVZt2eZINXz+zjd6ZHjbsymP5FntHxfp/V8bPWsvZ/5nJ1yEGXkTD13IZqk6oVcmSoPjOpoQ7OAkVauDgisBW3bAj8SNcn7VqO1/Nj/5sQrByM6m4G7HaAjgDz6jfKQ7GokIIPEI55J6vbVt3tHWT/Jd+ZnL5jrwHl4v/5xBNaF/P30jNKtkMbF+/3H2RNoGR+nmlogMFxXS67xuuPrqto88Tz6cc+L5nZ3rWVlgc305rVZiSHWGJ/7/xfX9HffRHVMVs7RTPPn/igk3c91n5VkWr/N813446mnjK1AoUz4CfhydEV2c01OAXqwcm3kjKXIv0bRAondVn3Q77Dkjen72O47o0tm19dgj22wj3Ec9etZ09edYKrPu63/j7YM46bht2SOn168f9Zmldwfznh+VcPagtIkJJieFv/46tf3OwxDo5Uu3oWG0BvA24SkQu8Q4GqRp4cTLIdFZQXEJeoXNzJwZ+kePZ7eUVFrt2xHntO79xweu/BL3Pt0OwshF3wq79heXKe/i/T5MXbo5y5+QxPczAH9+MFsE2qGB9YxXp8/zAu/65a3aU1tGLV9f7J5aZFztaO/fH1xnbjq+wlRJBwfR79FvL02rF69cgI5ivfyf2nWsksbS0/2/mat799eB7Gaz2pTMt+NGt1HDwe+PGGQUnnjOa34H/5AKBzvr3TMsjjf/+we8Rl9kVx+97x/7C0sLell9emDdCPJ2pU5bVBPAPoBvwBrAW2BPkohzgP21bou3cX2C55l7u7jw63vtNaWul/wYp2O8n2lp+sYqldM66HeE7UkfqBB+4PVi8eQ/Xvl22VKZvAnIRT12vx7+Jvm9SsB14ODMDpl+yo5XWN8XR6a/MKK09ZocDEeZg3WfjlG0+ybAdLygu4ecV22k96kvuTvCculZFmyAXxdESGvhcxbFMWxhD28yURdZG8gf7zvy6Mnl6Sll9u16dujzo4KCVW/dFXX82kY0Ay7fsZdqSLWFisXZbrPxf65686MpTJQOrCeBlwKXey2UhLipJFRWXlKtNZ8VRj3/PnR95EtBIR02++7+wOBLWv7Dq0ty9fPPnRo579ocyp83e/tkz0COekWBlTgdZ2MN/8+dGdls8XRGNbwOKrfqSt3Lbiyg2IMu3WKsv5ltltBPNO8WO+XZ/X7cr7o3tu7+uCXqaNljCUFRcwjOTFgc9cHn5+2URZ2WJvj0J3olioFO8tu7Nj2kbEQ0rH9f+CIl/LOuNpQ5fICtdZXxhxLv98J83dsqizTFNRRdt157Hv/kraBHlIU9N5YqAwtNOtXDGkqgPffoHHoqyi0C0gkUV6v21s3tWIljqA2iMGetwHMqCWHd37e/2fClXjR5R9o4IP+Q9+UWlc0xGcvorM6JZdRnb9xXw9/d/Z19BMQcKi6meY7VranSstHhd83bZU2C78wrL7PTjSTnimbrNZ29+Uel0Qivi7NDuZDmPcDbsPMAhjWok7PlCJYr//Hg++YXFjBzYBvDbqAdZfOKCzbzw3TI2787n8bO6l7nvyYmLqVO1Euf3bxkyhlCd+V//cQWdm9SM/CIc1ueRKZzWs6mlZZ08xXkgzjp+Tlm1bR/92tQtd7v/V8Wuxp8zXy27LR33y2puP75j2Mc48ZlMXhhbqSW75+MVytbYs2O7dcHrP9OnVfnPM1oFxSVlpkyMt/9yojmzp1VJy44kJJzf1+4sV4w1VUbZBovTv64YePqXfRSiFlckkU4NTV2cW6YOWjBHPfF9TM8dKxEJumdz4jP1JOie5wq1GS07SCD+IHb6Jffh1uZrHcsrCp6gxFqA+JEgU8Q5dYoqksDRkJ/MXc/Dp3UFoPWoLy2tw+6df7K448M/GDNtBTv3FzL7nmPL3f/lHxsZ0b2JC5GF9/TkJfy2ZgdvXNov6scu2pg8PbtK/HZb0bTohkoWpy/bVjrZQqDhz03jq5uOKr3+4BcLqVUlmyODDDD858fzmbtmp+V4ko3VQtArRWRFiMsyEflNRN4QkcOcDlhZE6pQbrT7llj6iS3YULbfXaR9UNn+gp6FF26If8qwES/8GPc6AsWa/AERm0VHvjGL92fHsX4HhEogwn2mY2esspww2M2OAupOH8Pvziuk7T+/tFSu5aZ35/Ly9/EVw42VHafqK5JluXvDfr++jKIQ/Kqt1lvvww2wsOL7xVvYHmFKw6ACtldFxSVlTk+Hkmr94AL9tWlPmQO94hLDzePnlTkV7ntrUjn5A+t9AD/C01pYA/gFmOD9WxPIBmYDhwM/i8jxDsSpgB1RjA68JmDQgVWBR/CJOqAPfN4TbUjeFmzYnZLFfRNBiG8QyHt+M4ucHWMphXCiHmhQXMIb01fS55HYK1VZ+a5HimvTrryIO8Alm/ZQYuClEH1b/R/9+e8beHLiYnKjnC3HSWHrriV4579o4+64Roy7xWpZFIitrEzgV/nyNyOPwo1Ui/OpSUs489UZ/BnhLEWskj1v3BBDX8xQ/li3kx+Xhh68kihWTwHnAkuAk4wxpVsiEakCfAGsAboCnwMPAhNtjlMBb0xfxYC29eJaR2J+ZLE/yYGCYtfqKdn13qTCWTBD2SmQ4hHtaGQnvDtrLQ9+YV9n8MDP0MpnOnfNDk5/ZQZPBPQRtEO/x77llmM72L7eQIGvc+32/SzLjX8QRazCve879xdwwvM/UrVSZuICChDrqX+7+99GOpizkkTuLwjouhNw/yLvvLyRWtkLSwz5IbpKxCKe7amd+7vNu+2bnvWUlzxzDk8fdQzNalexbb3RstoCeBPwjH/yB2CMOQA8C1xvjCkGXsNTLkY5ZLaFJng7+f/2Ji7YxCcWJmvPLyqJ+lfrW3rwU1PL9bsrt2yYVYdrgVi9bR+tR30Zci7U522a2SQVzF2zk5e/j1x70KVaxiGF+uz32jxyO9odhwgs9fZNmuVQGZD/zVztyHp9CoL0Dz7qie8t13Czau32/TxgQ7LuS1gCE5dEcrNMl7+fAmqC2nEQGus6xv2yhkPv+Sb+AGwQ6+Yr4tkRmw7ybx0/z54VxchqAlgbaBTivkZAde//u4DkHMZVQVgp1RCuxl60P+pVfpOLh5ozMdCHAcWHoznatbIxD7dz/u6v0P2qvvb2i/x4rn397HbE0rcmgbZanNA8VMFoJ9z9yXzHByMF822IkhqtR33JFW/ODvIIj10HCvnNe+BlZXBKJB9HOIhK9lNh4VgZBBLtzB6B/vbvGXS+L3SCkQot8KnMjYE+blUscML5fnM6x1LX0k5WTwFPAJ4QkV3ABGNMgYhUAk4BnvDeD57Wv+inNFCWvTF9VcRl+jwy2fL6nPj6lZSYMv1EChPYR2dbQEIW9DSnjS+618PW3+tk7o9opQK/XWat2sHcYGVRHH57Ln9zNrWqZAe9b8qizRzfJfgx7sg3fo2qs3ekr9c4b32/bVEOWImpI78Dwr2+cAdg4KkSMCnKwsKBZq0KfxbkpBd/imv9TnPjoDGWn1ait1ehvlfxxJGMA1JmLA8++tgNVlsAr8EzH/CHwAER2QkcAN4HpgPXepfbANxlc4wqCoKErEVkjClTZytwejI7+e8IIp26ivY3uiiK2T0uG1t+5JbPlj35nPby9OieXNki6Cbd+P8b+UuR7/ddttpCEK51/PYQp/Pmr4vc6V08o2qi4t+67s/O1o7HvkqO4t/geV2n6u+ttLh+sgv8HiZfKuWcdGlFtloIeidwqoh0AfoAjYFNwGxjzAK/5T50IkhlXbidxytTl5er9u7E9zwwglg7SoeydnuYDs0WtlK+Rcb/uiZkgd74WH9XndrQfBxPuZokEOqgwP/2x74qX0Mv6GNieH7/eZb9pzIL1aLw4BcL6djYU9y6zOwzLvu/6da6bURjaYwza+zPj347EPh2+/ffdWsn/fqPK7jiqLYxP95qcX1/4eb9ThSnp+9Mppwrls8oFm63UFptAQTAGLPAGPOmMeZx798FkR+lksXncdaTsvoDjbZuVSytJ8ks1I4pWHLu1O//X37TviXjkXuknfeEEDXVLvzvL6X/r9kefs7meKwIUact3MGMb4aW39ftjOu53dwnWDndFlgL81MLA8Mglu9h+Vh8oyfdFKx4dzQOFBSzc390CcYFr/8SeSE/dnyHAr8LvztysHxQvF0DfJ6auDjuxOrROD5jO2rYJkrIFkAR6QwsN8bke/8Pyxjj7IR8SlkQrgU019sy49s4PD15Sdh1/fPj+THFcPVbsdVgVAdZaW3Y4bcTtTNpCpcCRerjlo5ueW+eQ2s2FfJU3O/rdnHx//3q6HMs3Fg2CcmNoVXa/723UrTcSdH0A3zp+2XUqVaJy49sE/PzHYgwGDHcNJxTFm2mc1P3p3e0Itwp4D/xFHf+1ft/6D6anvvcK8akLFm8uXz5k8APNdzE40tjrAdmZdRYIrbzn/8eXctkpGnZohVsI5aEfZQTwu7ivckyQMIOFfUrYeU37v97GDNtBRMXhG4VStffTqL4b7YvfcPeUkDRirZf7MMTFoZMAK3UtbSrNTLZhUsAhwAL/f5XaS7WnWy4pBKI+nRIOHd+FLnVLpl2HJuSaIaHRDr/tSCntOI4CghXoijc0Xowy7aE3kFY+erY+X1OtGC1AN0SLvkDkmImhXSU6FbZZ6cs4d8/xFZcJNi23u66loE++m0dmRnW3iS3d0UhE0BjzA/B/ldKxWbsjPJJSuCgnHTm5Kj0aPwRYdTvN39an/M13bW/66vS/+1OHAJLPil7hTrt6kY5qwM2DyR00upt+1Nmu261DmAZInIy0BHPSOBPjTHBp1ZQSW9/Gk747tZRV6QWDZUarnn7N0fXn8i6mU7zH0H9hIWdotstIuqgwNlFUpFz/VMrhpCjgEXkThH5MeC2bBGZBnwKPA68CfwpIk0djTIMEXlAREzAZZPf/eJdZoOIHBCRqd5yNhVStKc3v1+cnqdRkn0GD5UYyThi7+5PYxt8lOx8BbDtksyF1VVkYx0oUaSiE64MzOl4ijz7uwk4EngEqImnJmAxcLcj0Vm3GGjid/Gfj/gO4O/AjUBfIBeYLCI1Eh1kItw0fq7bISQ9Y6KbwUNVXMnY2TuvsOK0ADqpIo4QTid2zAcdztn/meno+iuCcAlgO+DngNvOBVYaY+43xuw1xvwGjAaGORWgRUXGmE1+ly3gaf0DbgFGG2M+Msb8CVwC1ADOdy9c56RyB3Slkp3mHM7avi95imir1Pbryu1uh5D0wiWAVYGdvisiUh3oBUwJWO4voJntkUWnrfcU70oRGS8ivjLtbfDMWjLJt6Ax5gAwDTjChThVEnBidgSVHkq0k5qjpi+zPk+qm6P5Z1SA/nGxyi9KnQEZKrxwCeAKoJ/f9WF4DoADE8BagJsdaX4BRgLDgSvxJHwzRKSe93+AwPM8m/3uK0NErhKR2SIye8uW9Owfp1RFc/9nf7odgrKZm6eAz49yZo6K5LkpS90OQdkk3CjgN4AHRKQIT8L0ILAF+CpguSF4+uC5whjztf91EfkZT/J6CeVPYVtZ3xhgDECfPn30eF+lPLfnm0wGb85c7XYIqoKZtUpPMarUFq4F8AXgXeBfwFtANnCeMaa0qqqI1MKTaAUmha4xxuwFFgAd8JSpAWgUsFgjv/uUUkqlmP0Rputy2qZd6VnEXdnH7WPzkAmgMabIGHM1UBtoaIxpZYz5PmCxfcAhwLPOhRgdEamMp0bhRmAlnkRvWMD9RwEzXAlQqQTboQODVAX0wrd6KlKltjyXC1xHLATtHTQRdC4vY0wRYL3XrgNE5CngC2AN0BC4F6gGvGmMMSLyHHCXiPwFLAHuAfYC49yJWCmllFLp7q9N7s6hEdNMIEmmOZ5T1fXx9FH8GTjcGOPr9PMEUAV4GaiDZ9DIcTp7iVJKKaXSVcongMaYcyPcb4AHvBellFIqbt8vznU7BKXiEm4QiFJKKaWC+Pi39W6HoFRcNAFUSimllEozmgAqpZRSSqWZkH0AReS7aFZkjDkm/nCUUkoppZTTwrUAbgu4HIKnfl5VPGVUqgJH4im4nL4TIyqllFJKpZiQLYDGmL/5/heRy4FDgSOMMWv8bm8JTAAmOxmkUkoppZSyj9U+gHcD9/knfwDe6w8Ad9kcl1JKKaWUcojVBLAxkBPivkp4ZuBQSimllFIpwGoCOBV4XET6+N8oIn2Bx4EfbI5LKaWUUko5xGoCeBWwHfhFRDaIyDwR2YBn2rXt3vuVUkoppVQKsDQVnDFmHdBbRE4E+uI5JbwJmGWM+crB+JRSSimllM2imgvYm+xpwqeUUkoplcKiSgBFJAdoBlQOvM8Ys9CuoJRSSimllHMsJYAi0hQYA5wQ7G7AAJk2xqWUUkoppRxitQXwdaA3cBuwEChwLCKllFJKKeUoqwngQOBKY8z7TgajlFJKKaWcZ7UMTC5wwMlAlFJKKaVUYlhNAO8D7hSRmk4Go5RSSimlnGf1FPAZQEtgtYjMAnYG3G+MMefYGZhSSimllHKG1QSwPrDc+3820MCZcJRSSimllNOszgQyxOlAlFJKKaVUYljtA6iUUkoppSoIyzOBiEgN4FTgEILPBHKHjXEppZRSjukvi6gru/m6pL/boSjlCqszgbQDZgBVgGrAFqCu9/E7gF2AJoBKKaVSwns5DwPQOm+cy5Eo5Q6rp4CfBWYBjfBM/XYinmTwQmAvoCOAlVJKKaVShNVTwP2AK4B87/VKxphiYJyI1AeeB45wID6llFLKVrXZ43YISrnOagtgZWC3MaYE2A409bvvT6CH3YEppZRSTjgv83u3Q1DKdVYTwCVAK+//c4FrRKSyiGQDlwMbnAhOKaWUslsx4nYISrnO6ing8UBP4C3gXmAisBsoATKBkQ7EppRSSimlHGC1EPQzfv//LCJdgeF4BoJ8Z4z506H4lFJKKVtp+59SUdQB9GeMWQu8ZnMsSimllOOyKHY7BKVcpzOBKKWUShv12cXt2e+XXq9CnovRKOUeTQCVUkqljcayrcz1q7MmuBSJUu7SBFAppVTaKAnY7WVT5FIkSrlLE0CllFJpIYcCLs/6usxtGRiXolHKXZoAKqWUSgvXZX3OmZk/lrlNRwSrdGV5FLCINAVOAprjmRnEnzHG3GlnYEoppZSdqpTOZqqUspQAisjpwLt4ij7nAgUBixhAE0CllFJJqQnbaBIwAMRDTwGr9GS1BfAxYBIw0hiz3cF4lFJKKdvNrHyj2yEolVSsJoAtgBs1+VNKKVWR1Gav2yEo5Qqrg0BmAIc6GYhSSillt5rs47as90Pef27W1MQFo1QSsdoCeBvwjojsBSYDOwMXMMbstzEupZRSKmZN2cpLlV5ghWnKWZnT3A5HqaRjNQH8w/v3DUL3mM2MPxyllFIqem9kP86QzN/pkPc/Csmib8Zf9M5YRm+WuR2aUknJagJ4GTpUSimlVJIakvk7AG1lA4tNSzIpcTkipZKbpQTQGDPW4TiUUkqpmOT4VSarI3vBQFXRmn9KhWO5ELRSSikVqLnkUmCyyaWOK89fn100koMFKmqw3/v3gCvxKJUqQiaAIvIrnrp/C0VkFhFOARtj+tkdnJ1E5DrgdqAJsAC4xRjzY/hHKaWUAmgn69lnKrOJeqW35VDATzm3ANA173X2UjXu58mkmCyKyadSmdsrUUh9drGB+n63Gr7JuZP6srv0lnqymxMzfqaG6LhEpcIJ1wK4AEoPoRaQwn0AReQc4HngOuAn79+vRaSzMWaNq8EppVTSM3ybczsA/fJeLm3t652xtHSJ+7Le4o6iq2Nau1BCC9nCGtOIl7Jf4ITMWbTOG8ftWePZamrxRvEJLKl8CQDnF9zFcRmz+b6kFwtKWpdJ/gBGZb1LbdnH8pImUT2/QTgj40cmlfSxJZFVKtmFTACNMZf6/T8yIdE45zZgrDHmNe/1G0VkOHAt8E/3wlJKqcTxJDq+8q+GHApLW9pqso+uGSv5paQTguGCzG8xwP+Kj2dQxh+l6/gh51Y65Y8lgxLerfRo6e1nZ/3AnUVXYsigMvmckzmV5aYpLSWXx7L/y1UFt/JzSSeezX6VaSXd+aj4KI7M+JPvS3qyuPJIAO4svJITMmcBcF3mZ1yf9TkA52ROLX2eazK/YFDmfEYyiZPyHyn3GmvLPgDaZWy0/L5ckzmB27I+IFuKAfh30UkcnfE7txReD8DEnFEAbDU1+bOkDRtMPTaZutyW/SEAE4oPp2fGMsYVDeXN4uO4OusLJhb3ZbupyduVHqNdxkbeLzqax4vOpU/GEjaZOuynMq1kM6dkzuCnkq60lU08V3QGed7P4/CMRfxV0oIShGIyKSKTARkLmFbSg8NkCQeoRI+MFUwuPoz7s9/kwcJL2E4NqpLPPipTTAYglt8DlX7EmJRt2LNERCoB+4HzjDEf+N3+MtDVGHN0qMf26dPHzJ4929H4/nX3NdSWyJXojcUfspVP0/q6Ii9n9dtj7TltfI0mOV+jnZ+j9fcr0a/RznVV7PdLMDSSHTSXLaw39alGHnuoSglCXfaQgWEjdekkaxAMy01TarGP3VTj8IwF9MxYwbiiIXTPWElT2Upd77ZknanP3JL2nJz5c+lzzStpR8+M5UHjeL/oaM7O+qHc/wDvFA3lgqxvyyzfLu8t2sv60sTojaLjuTRrIg8XXsiMki58nWP9uHpRSUs6ZcR2IubKgtt4rdIzpdf3mxwd/GGjpSXNWGRacljGEpp551FeVtKUtrKRfLIZW3w8lSiiCvnspQqNZTs7TXX2k0NnWc06U5+9VKUShQAUkA1AIZlkYsihgAJvO1QelciimExKEO8vp5BMDEIWxYT7/aZqFrPR1OWhR19w/HlEZI4xpk/g7ZYHgYhIa+BC4BCgcuD9xpiz4wnQQfXx1CjcHHD7ZuDYwIVF5CrgKoCWLVs6HtzZmVNpLlsjLGV19xZ5OavHg9bWZd9uN0NS9SesVHxKjFCC51LJ2wJVYoR8ssmhkAwxfonNHLab6tRmX+lv5vys78uts7lspXlm2e1KqOQPKJPw+f8PlEn+RhVewejs12ktm+iesQKAY/KfYq+pwqVZE7k3++3oXjzEnPwBtJKym/U9VKEqmgDapUPGejqwvsxt7TM2AFCFAq7N+oK9pjKVKESAbCmmyGSQgSn9fvon5QUmkyKyyKEAg7CfHLIpRjBUoogiMigms3Tf4kn8oChMmeHw+xcTcQk3zTdtXH1+SwmgiBwGTAPW4EkA/wBqAa2BdVBxKm0aY8YAY8DTAuj08w0teNrpp0hDFT9htr4u+77Cqfwa7fwc7WxvFAz7qVzaCgJCpnenV0wm1ThAY9nOctOURuwgn2x2UoNMimkmW1lrGlCNPPZSlbrspnvGCmaWdKaDrKMSRfxmOtBJ1tAv4y++Ke7LAXLIoYDuGSs4KfNnXik6lXyyuTHzEz4sPpqVpjHTcm7hl5JOXFV4G6dn/sTo7NcBGJj3PLW8p1c7ylp6yHJ2myqsNI2Dvrbni07n5qxPACgyGWSJpy7fl8X9aCC7+L2kHZOLD+P9nIfLPXZ+SWu6Zawq9///ioZxcdbk0uX6ZCxht6lCTTnA0pJmADSSnRbeeWtmFnemc8YqMjBsMPU4NGMdJUbKHbAuLGnF7yVtOS/re34tOZR+GYtDrnN1SUP+MG1Zb+rTUdbSQHbSJWM1ADOKO/NU0dl8nPMALxadxkWZk9lpqvNy8amckfETS0wz9lEFgE6ymh4Zy6kre9lqavJF8QA2mzrUln1sNrVZbFpQiSL2mcrkk802U5OWGbnMK2lHNfLJoogacoAmso1lJc0oIJsShKqST1O2stI0Zgc1qEQRxWTQQnJZb+qTRQk12M82alJEJsVkehNAQyFZlCAIhhocYC9VKCHDr+uBf0KW3MlZoqxy8bktnQIWke/wJH+XA4VAH2PMbyJyBPAucLUx5htHI41Rsp8Cbj3qS0fXr5RS0SjbTxBuzfqQ0zJ+4uiCZ6lEEQtzLuW14hGcmjmdRSWtuLzQMzikm6ygnuyiV8YyXio6nUKyAENz2cp6Uw9DBtXZzz4ql1l/dfazl6oMzZjDfys9zYn5j7HItGRl5QsBOCv/Pj7MeQiAQ/Le5ILMKdSSfdyS9TEAf5W04KyC+wH4X6XR9M6Ivj1icUlzfizpRlXyeK34JJrLFtaaBqwyBweSVCaf+7L+xxvFJ7DUNKcBOylB2EatkOutxV4qUUjfjMUsNi1YbppFHVskvtOo/u+pSh2rRo9w/DniPQXcE3gcSkurVwYwxswQkQeB0UBSJoDGmAIRmQMMAz7wu2sY8JE7USmlVHIKTCSeLTqLZzkTEArIZrVpxLVZXwDwaPGFpcvNN23BwNSSXn6PFtaZBqXXgo2u9d32bclhtM4bV3r70Pwn2W5qsIOaXF1wKz+VdKWAbN4oPoHjM2aVLrfD1Chdxz5TrneSJRcV/LNMHcOVpvwI4jxyuKvoytLrW6gdcb27qA7AVyWHxxSXFYHlcpSyyuohgwEKjKe5MBdo5XffWqCD3YHZ7BlgpIhcISKdROR5oCnwb5fjUkqpFHDwVJ3vNDDA9JIujj3jctOMHdQEYGJJ39JTnwCTSg4r/X+HN8kCyk3/trqkoaXncquItVJuspoALgTaef+fCdwqIh1EpBVwBxC6d3ESMMa8B9wC3APMA44ETjTGrHYxLKWUSjmfFg8EPCNEd1LDlRj8Wymb+M0CUlk8U8LNL2kNwDZvAqmUKs/qKeAxeAZ8ANwFTAL+8l7fB5xlb1j2M8a8ArzidhxKKZXK/lV0Pv8uOoWtLidXE4r7c1LmL+R7S4sANPQOALmn8DI6ZKynAbti6hOoVDqw1AJojHnLGPOw9/9FQCdgOHA60N4YM8m5EJVSSiWLYjLZSi3cHsH5fvFgTzzm4G7sp+KuACw1zfmw+Gi+KDmc6cXOnaZWKpVFTABFpLKITBKRwb7bjDF7jTGTjTGfG2NynQxQKaWUCrTQe5p3TPFJpbfdW3QZR+S9wH5vqdp1piEXFN7tRnhKJb2Ip4CNMXki0hfCVGJUSimlEmgrtcqMGgYoJIsN1HcpIqVSi9VBIJ8DpzkYh1JKKeWIrUYHgygVyOogkInAkyLSBPgKzzRqZSpIG2O+sjk2pZRSKm7nFNxLT1nO05W08pdSPlYTQN8Ej2d4L4EMeopYKaVUElpumrHcNKN70XIu8ZtKDuCuwstdikopd1lNAN2dsVgppZSK0/1Fl5ZLAH8q6epSNEq5y2oCaICNxpjCwDtEJAvPrBpKKaVUSjGRF1GqQrI6CGQl0CvEfT289yullFJKqRRgNQEMV/GzMpBvQyxKKaWUUioBQp4CFpHuQE+/m04UkY4Bi1UGzgaW2B+aUkopZa8LC/7J25X+VXrdWG4HUapiCdcH8HTgfu//BrgvxHIrgavtDEoppZRywk8l3cpcLzHuTmmnlFvCHfo8BtQAauI5BXyM97r/JccY084YM8XpQJVSSim7lbg8p7FSbgnZAugd8esb9att5EoppSocowmgSlMhEzsRqSoi5cq7iEhTEXlaRL4Ukf8Tkf7OhqiUUko5Q1sAVboK17L3NJ4p4EqJSCPgN+AmoBFwIvCDiBzmWIRKKaWUQ/LIcTsEpVwRLgE8Cngz4LY7gAbAicaYPkBr4BfgXkeiU0oppWw2JP/p0v/3UNXFSJRyT7gEsAXwR8BtpwFzjDGTAYwxecCLQG9HolNKKaVsttI0cTsEpVwXLgEswa8AtIg0wTMn8NSA5TbhaRVUSimllFIpIFwCuAA4xe/6GXjqAX4dsFwLINfmuJRSSilHzSo5xO0QlHJNuELQjwOfiUhLPK185wO/U74F8GQ8A0OUUkqplNAjbwx5VHI7DKVcE7IF0BjzBXABUAcYAHwEnGKMMb5lRKQB0BEY73CcSimllG12UZ18TQBVGgvXAogx5l3g3TD3b0EHgCillFJKpRSd4UMppZRSKs1oAqiUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaSZkGRgRKcEz84clxphMWyJSSimllFKOClcH8CYOJoDZwN+BvcBneKZ+awScClQDnnYwRqWUUkopZaOQCaAx5iXf/yLyDPAL8LeAmUBGAR8AbZwMUimllFJK2cdqH8CLgdf8kz8A7/XXgAvtDkwppZRSSjnDagKYCXQKcV+XKNajlFJKKaVcFnYuYD/vAI+JSBbwOZ4+gA3x9AF8CPivM+EppZRSSim7WU0AbwMK8SR7j/vdng/8B7jD5riUUkoppZRDLCWAxpgC4FYReRjojmcE8CZgvjFmu4PxKaWUUkopm1ltAQTAm+xNdSYUpZRSSimVCJYHb4hIdxF5T0SWi0i+iPT23v6oiJzgXIhKKaWUUspOlhJAb4I3B2gM/A9PYWiffOBG+0NTSimllFJOsNoC+C9grDHmaODRgPvmAT1tjCmtnN2nudshKKWUUirNWE0AOwLvef8PnB94N1DXtojSTMu6Vd0OQSmllFJpxmoCmAu0DXFfF2CNPeEopZRSSimnWU0AxwMPiciRfrcZETkEuBNPoWillKrwalfNjryQUkolOatlYO4FOgM/4Kn/B/AZnkEhk4DH7A9NKaWUUko5wWoh6HzgJBEZCgwF6gPbgW+NMZMdjE8ppWJ22cA2/N/0lW6HoZRSSSfaQtDfAt86FEvURGQqcHTAze8ZY871W6YO8AJwivemz4EbjTE7ExFjJCZwSI1SyjZtGlRzOwSlVAXRql5VVm/b73YYtrFcCBpARHJEpK2IdA68OBWgBW8ATfwuVwfcPw7oDQz3XnoDbyUyQKWUOwa0red2CEqpCuKMXmXLttWsHFUbWtKxWgi6qYhMAPYDS4H5fpc/vX/dst8Ys8nvsst3h4h0wpP0XWWMmWmMmYknQTxJRA51K+B41a+e43YISqUEEfvX+fiZ3e1fqUp5/dtoNbR007FJTbdDiIvV9PV1PC1ntwELgQLHIoreuSJyLrAZ+Bp40Bizx3vfAGAvMMNv+enAPuAIYHEiA7VLszpV2Lo33+0wlLKkTtVsduwvDHpfm/rVWLl1n2PP7UD+R9v6elo5nYkE77qjB+YVnxMHlG6yegp4IHCTMeZ5Y8xkY8wPgRcngwxjHHABMAR4GDgT+Mjv/sbAFmMO/ly9/+d67ytHRK4SkdkiMnvLli2OBR6Pk7s3cTsEpSyrlhP6OHNox4YJjESp+N19Yqegt5tycyQoldyiKQR9wMlAfETkERExES6DAYwxY4wxE40x840x44FzgGEi0jvW5/eus48xpk+DBg3seVE2u/zINm6HoBKgbQUZwBDuqNnpXaZUtEN2lbR0QF/66dmittshxMVqAngfcKeIJOKE93NApwiXX0M8djZQDHTwXt8ENBC/vYD3/4YcrGeYcnSnlrzaWUjarj461KQ6ZQ0+pOK3jnV2uA9Ndqb+VpS9dPurAL686UjuOD5lhxIA1vsAngG0BFaLyCxgZ8D9xhhzjh0BGWO2AltjfHg3IBPY6L0+E6iOpy+grx/gAKAaZfsFKmULKzuHRjUqJyCS5BGuZaRXy9qOPneDGtovS6lEe+n8Xtwwbq7bYTiqS9NabocQN6stgPWB5cA8IBtoEHBJeFOFiLQTkftEpI+ItBaRE/FMWTcXz0APjDGLgG+A/4jIABEZAPwHmGCMSckBICr1+edDWRnampDsmtWuUvp/q3pVtaeXCqp5nSqRF0pCfVvXsX2dJ3VvyrGdKv4ZjFRnKQE0xgyJdHE60CAK8MxKMhHPaN4X8ExLd6wxpthvufOB373LTfT+f1FiQ1VOOqVHU7dDiFm4ZELPNCWHZn479kEdyvYLPqN3M/56eHiiQwrriHb1eOpvPdwOo8IyIZq0K2VFVVY3aRzZPjn7uivnRf2NFY+mIuJqBURjzFpjzNHGmHrGmBxjTHtjzM3GmO0By+0wxlxojKnpvVyYLLOAgPOd4Cui6wa3K3PdyUTp6Sh3pFZCqVstO7ZgVNIRhMrZmW6HUcYNx7TnzN7NXHv+bs1S/9RYOJkhWu1TdRCIjl5OX5YTQBE5UUR+AfKAtUB37+2viciFDsWnVDmJ7HvhxKbx2E6NHFirSoQ+Dpwuq2huHtoh8kJxcLtlPEOEn+5046RXZKtGj3A7hLRwzdHtIi+UAqzOBHIxnjl0/wKuomxDxxLgcvtDU3Y5vZd7rQHpJtTO6de7h/otk17ndlP+5fodBZzas1nAXZ47q1VKrlbAimzlv9xPcprXqep2CLZp4X0tPZpHd2Cd8r/rOFSU6gJWWwDvBp40xlwCvB1w3wLAzbmA08rDp3aJ+jHPntPT/kBc5t8xP5mEGnXa0G/kr38folD9iWJVPUzRZbeEe4mO1wF0YC6QirHpV7Gy+zfrtkMb12Da7UO4dKDWl43kcO/c4hVljnGrCWArYHKI+/KA1J4QL4VUT5LJp0/sFnQiFUdMuW1Qudu+uumo0v/t2CFbKRfynIVE+p8ndOKx07vZEFFsquVoS5RVfzuseeSFgHYNI9d2TLdW3XCceiuOPqQBP96RnKdendShYXXHn6NlvapR9wW0+jFH27KYzPq1qcviR4ZzRPv6bodiC6sJ4FqgV4j7+gDL7Akn/VyWorN6PHlW4kYZtm9Yo9xttap6BlIcadMPsWmt4LX5hneNLtGtnJ3B+f1b2hFShZHqudHlR5Yt3F2x2n9Sx5uX9aNFXedPvU67PXyS2S4BCVk83rysX0yPc6ph84ZjnO0TGsiJsjb+crLKHmRflsItp1YTwP8C93sHe/jOvYmIDAXuAF5zIrh0kIyn7KwINb/rGXGOPrzpmPaWi/f+cPtgXru4j6OtL258Pr5BhimeN1kS6jW+fnEfR55vzj3HHnxui29wKiSw3StQK4vb6lWvFPK+KbcN4qgOwcumNEmSbilHHxJbWZeSCAlgm/rVqJkkZ6DCyUjwDzaVi81bTQAfB94C3gR8ZVZm4Kmr954x5gUHYlMuGN4lvlO7/zwh+ETpoRzVoWwLXnZmBoe1tHYE16peNapUykyyRCm6aLo1r11+DamQcdigcnbozc+xne0ZKe3/Vr55WT/qVT+4sXaif6BKfaEObiH42QjwtDpd0M9ay/9L54c6mRbcxUe0jmr5WJVEaAIUYM69wyx1hUlmE2480tZWwlQuo2O1ELQxxlwPHALcANwD3Ax09t6uEiTcTuuIdvF1TG1dryovnBfdxilegc3pEHlD5AgLSdeAON/fYP4X4+maaCT6M01WTk664nYqeduwQw5ecXl/lCbHL2UMPrQhGRa/YCd1j65w/UWHt4olpOhZ+N5kZ2YwonuT6FZbwQbNBIrl5QU2fLjFahmYQSJS3Riz3BgzxhjzmDHm38aYJSJSTUTK99JXCRdvU3SdapUsVbN3esqjSKci3FK3WuhTQ9Hwf3m1qpQvCh3PBvOl83uX/t+6XvT9pYLF4xS39wtWExUrcbZ1uV/Y4EMbOj4y0c73K15uTKF414kdbVtXrKdSnaw/F+7AOytDeO2Sst0yIp+pSMMjAYveury/2yEA1k8Bf0/oUi8dvfcrl7StH3mUohWB01wlwgVBBky40QLo21T1blnb0vKhBo2EO60ZrcDt67JHT4j4mL6t65b+/8l1A5lw45GWnqtjY8+pLbv6kt0zohP1vadbrxvc3pZ1RuvR07uSnRn687CzpWrsyL5ceVTydAZ3ohuB72cZaZBERXXR4a1tW9e0GEczO5n4htvq3n9yZ9o1SPxBTtsGB/dt5/RpYekxSdp+kJSs7q3CfeuqA/ttiEVZEKy/gW9HWyWOKak+v2FguQr+oVr6rhrUNujtsRjSseyE4Qa4dGDr6FZi4zaxl8X+h52aBK985GSB2KwwyUwwdapVomuzWpZaFI/r0phVo0fYVl/xiqPa0rNFbQDqBbScPnDywWPJYJG9ZsMAkDN6N+OC/pFOndn3xalTrRJHH9Iw8oJhPHFWd/51hnslhKxq5vAZACt8+a0LDYG2qF01trMJTp5ed/LA+87hHRl3RXStXn1b1ylTUqtRiIPuQG6dck7k2RO7hNyjeE/73ici93lvusJ33e/yGPACMD8h0argvBuFrs1qcUavZuV2uFY0rlm5XB+Wn+48ptxyq0aP4OIBrUOux2C4fkg77j4xusEg/o7q0CCm05f+nNhQpui+xhKnSyf4y/EeqITaTA+LMADkPAud7QcfWj4ZC2zBsPIdefT0ruWWs7J/iaU0RP82dS29tmSXiD6Avs/gxfN6h1/QRuFe14hunn5xdao6mwTYefAdKNq8KZqP+drB7WKqnRfLvswKKwPARkY5+CYVy3+Fa1LoD9zovRjgb37XfZdLgK14BoYom3x+w8CYH/vMOT2ZcJO1035Ouf34jlwZ54Yq3mO4O4dH7q/zyXVHxPksqSHS6cB3rugfsrSFE+LND+4cfmhUy9f27pQDW8itxFGjcnbYHaNvxw/QqYnnNPrDp3XlvpOTZ3KkBQ8eb2k5pwv2Ln30BKaPKn9QGY+hnRpy5VFtElIsOZzW3m443/19sKPPU6NyNt/cclTkBWNgteUs1GJOlG7q0KgGg6Isa2NXA+ADp3QJ2dWn7PN5njAVGwhCJoDGmCeNMQ2MMQ2ANcAQ33W/SzNjzFBjzG+JC7ni6x6kNIgVvv18k1pVuGdE7C1wyc7KaNxKmRkc0ij8TsH/dK9/jmRXcelUkYhah3b2jbQiERvjQxpVZ2ing62V9arnsGr0iMSN2rQoXFkTu6waHXl+3uzMDNuncMzMEO4e0Tls7T67WGnZrOPXYhVpO/LBNQNiiqNjY2cm3op28F3g+2FX6aZAVvtl+8SS//Xz6zsdq1Qc/W61DEwbY8w8h2NRfk7rGbxUQLCm60pB+oZdcVTZFrhI/QMzk6QzTXsLR/KBo3FDNefHUufNGBh7ad/g6/P7hcdzkFmzcvDTREe0q0eGwJm9rU1RZhf/DVeojVi4o/ubhkau9O/fP6b0OaJ4E/2T1Mpx9HX1Z2WDHa5VJNaWhosHBE8QG9W01scplETVI4s89tP+bYmTB2WNI7zvsQ686Bshqejbui4tEzCziVVhv+uxrTHWUOISSx/Ad67sz6fXlz/zVtFrsobrA3iiiNT0+z/sJXEhp4fnzu3FOxY6zT50ahdL0yNFmtLMv0Cuv6n/GBxx3XY6sVt0NaYCHdrIcxqud6vo+rT5/8yjHWxhl5Z1q7LiXyPoEKHlMlp2dIoOd3Qfzc65X5u6DO8a/Wd8eNt6rBo9glWjR5BjoVSRFbEmKvHsExrXrMzAIO/XJQNa2ZbY2p1/3XrsIZEXcshhreow995h/HekMzPDQORBR75fjxuFwwe2d7a0jz+rLYBRzxkc5Q9m8KGeU76RWq5/u3dY0INPS68jIKTszAx6tqjN8sfKpjJWtp3+i6TambdwW9IJeEq8+P7/wvs32OULB2NMGw1jqON3bl9rHU9j3XS1jrbEjItj8J88qzsTbx3E4keG07NF7ag2VL6NVCpXdU+UwBbYJhb6ybx0fm8GH9qAd688PO7TwSJiqSRO5PVYXdK+78Qn1x/hev3DaN18bAdWjR5hvQ6gje/Xh9cMoE61SuUKxrvxFkab+Ft5HyI1LtrRN/ftIDXngtU0jfY9DZYQX9C/JaPjGMl+Xr8WnO8dCJWV4dlOhPq91K1WqWwBdK9Ir6N2mIE6gWfCojktLki5M2/JLtyWuA0wz+//tt6/wS6p9aqT1EOndrXUlyYaTnaOtnLqLxFuGtqeni1qc5x3Grtgs4skm3tGdErZQSi/3TuszHUrLdB9W9dl7KX9yMwQsjM8fcFGnxl5RxGq/48drbS+TX21Ss58X4K1VDapVbYP3KRbB9GgRg4XhRlZHzWHsqNkTFwT0SbnZFmRZxIwrdqRFmed6N/GWj+4cG/Ho6d349w4RrIfHqaYeajPOtqZSaIp1WUliQ9cwr90TbILNwhktTGmwO//sJfEhZx+QtWci+SH2wfzsY1JxhsBfeOCHX2FEthPsV0De4pXg2dO4E+vH1iuDlOsfQBDsXNnc8VRbS3XHEws53epGRnC9FHHcEbv5mXe7+o5WZzS42Df11WjR4Sd+SDcpPfhjvJ9fP1No5083mo+EKmExXGdG3FIoxrMuvtYS31fI7HjFKWVucAjvV2pMMfyfy46rPT/SDv5g6eArfF1ifD/npwUIklxelalUIIltV2bRTkKPEk+5pfPDygFZOEHWtdiHUZfC2DVMAeJvqfz/S7O798yZEPON7ccxY8xFgF3QlSH0SKSIyJtRaRz4MWpAFXstaVa1atGDe+Ag2A/iRO7Rd7Y+xsc5XB8f12blU1i376iP7/ePTTk8snY2uDPyVYBu1ddtVL4vjTJssP+7d5hUc1b/GaQeZRP79WMf1/Y29Kps1N6NgOgYc1IXS9CDDJKjrfNVvEUs33zsn62lAI5o3ez0v+d6oR/vIVEN5DVWNp4u834H4Sc2rNZqMVTQp0YC1eDtW1luDNV0WwOOzauUa4EU7C+9I+f2Z2HT+0ScX2+0L+66Shm33Ns2GWtfDs6Nq5p6YxJolidC7ipiEzAM+PHUjyFn32XP9FC0AkTaTBHNC4/MnFn7gM3nlWzs2hYI75RjxVFqD500e783hgZfPTysZ0a0qKu+7M3JEKlzIxyg0yO8RaFzgnoe1irSjbPntOD/0Wcl9OZMZBJfowTtaMPaWCpFEikOYufObtn1M8dba24aER7QHb3iE6MvbSvpWkV7Tr4stIP1ydUa6Tdj4mmifBQ71SUYddmYXVvXtaPNvUjt6bXqpptsduF58OvlpNVOuNWRWK1BfB1oA9wGzAcOMbvMsT7V8XIN/OBlS948JGCqbUreepvPajlcMX8aFnZVJUtl2Jfy0Q0fSkPa1WHGiEmkg91GlFEXKlNVy3H811N1Hy8odY3+szuTB91TNCW0NN7NY9p8FU0vv370WWuJ7Ll0L81LV59ohxZH8qZhwUvc9SsdhUm3jLI0joCW5WuDdNVwC5WP7bK2ZkMPrRhQs5i+LoHvXrhYfx+33GWHhNLya9g2zsnvsZln8f9/VpJwOndYGId9PTGyL5Muc3a990pViuEDgSuNMa872QwKj7RH00m7gfmv8E+K8QOIF1Fs6P46FpPn87Wo7505DnCbejqVqvE3rwiy8/50Cldad+geti+enYLFn+lLPsLEEcjnjm64zH/geNsfe7/u7Qva7btj3zwE2NmEMssIb5YAmPq2LgGf23aE/X6wg0ISuZT/lkZYvmgul+bukxbsiXsMoc2qkHzOlX49q9cS+u0M5Gx+7R/PGuLapaPKOMe0jG+ucPtYLUFMBc44GQgSrkpGTbuTsYQLv9r7Hf6KFyi+MtdQ1nwkLVpxcBzmuXGoR3KzTGd7C4e0CrqkYXBuN2PtUbl7KhHS4f7DtasnB39QIE05n77lefgx9/Mfx5TWmYlnMqVMvlviC4lOVkZnNC1cZkuJ+0bRj6Fa1Vg626b+p4+c63rWRs4GPj4eM42+dYUbqCY27/zeFjdOtwH3OkrDK2c4eQXKZ6Op74p1VL5ix5Kl6Zlv9JODu5wU7iXVd/iNFrZmRlhT+fGy476cX/r0yLudVTOziyzwY/3K+H0wYXV06aJEunl9m2djKPfQ2vkHSSUirNCLHmkbL3MJrWqWHsdYb70IsKrFx7GETHMzvJiFIO8fN+k03o246NrB3Bqz7IVAqwYd2V/ujStxRi/Ud/RKPGeAw5/Ctg/2tQS8hSwiASe7m0JrBaRWcDOgPuMMeYcm2NLG6G+OFZ2PCLWlrvpmPb0bFGLy8bOjio2gHevPJxFG/fE3ZIzonsTx2quResob22sVvWqsmDDbkunz8tMBWcxK/j1rqHkFZbEFqSN0qHA9ZVHtaG3hdI6x3dpxKqt+2N+HqvvZbzvec8WtZm3dmfE5ax0oE+kSDPDeKoTZLEniu4E0bjo8Fa0b1id+z9fYGn5SD/lj68byG+rd5S7/eubj+KE53+Meb2Q/L/Lj64dwI9Lt9q2vpN7NGXMtBXMX7+rzO3h3gUR4bBW1ufq9V+Xb8DRcTGM+vZfl7X9Q0xP4apwh/MNgQZ+l+V4CkNnB9zewLusclm4L2BWZgbHdLQ2WXdgyZZ61XMsFxMNRUR4+fzePHFWj5DL3DCkfcj73rvq8Liev2ws8PolkctVnBpiPuZoNKxZmZb1Ym99jWbWjLBHqQ7tZ8KV8rHCzristtD856I+TLw1dKuZ3XMxJ0uZnUAzgvS5Oy+OIr4+iZhKcbB3ZHewd/bh07oG/d1EU7fUX7PaVTi5R/ltQTT1WZ06sxDrN8tqNIe1qsstNk8F+P7VA0LeF+8vxf/xdatVir/V1jcIJMxXuqP3AOyQRmUPxM7v35JjkqCfXzghWwCNMYMTGIcKwq0jii9uODLmx8azmfv7cQc3NIFHxlVsbDmsW7X81FJ/69Oc5nWqcLN3Y/fXw8MdPd3pL1wrwJx7htHl/olRrzPSad0R3Zrw5fyN5W6P9juXTKV84v25ZGUIyx4LPa25rx6n1YSuU5OabN69JeTcxU4l5eFWO/ueY+nzyBQAmgYZGNOjRW2WPnoCP6/YxkX//dWZAIn/s3rlgt7k7s6PeFZiYPt6TF+2rfR612Y16ROkNemdK/pzweu/xBlVWYc2Sq6WWSustjrHw/q2vPw3+Yh29cImVXb/pEosDAI5sVsTJt4yqFxLvG9GkGgH7CVSyD2ciBwuIslVq0PZKtQOSESiPnKyMvNCPLo53Pm8aqUsnjmnZ+kcmZWzM2MqlxAf33zEB1XLyeLmoR2Cto6GKu3y18PDy42oDGyBOMymkh6JE3nTXifCzBuh+D7lSC069aKsA/bieb0Yf9Xh1KlWiaqVMktnBon3WzX20r78cPvgcrdb+claqWWWnZlhyxy0Vg05tAE1cqwWpPConJ1pqWW9We0qZX47E248igdOKVsA2OBJfO1Wq2o2x3aKrwUo0kGC3QlP4HvjhnD7nnFXHm55vl07tt5WBoFA8nXDsCrcr24GkCcic4Dp3ssMY8y2MI9RaSpwqrdA8Z4CcaoDtm/mg2hOtSaK7xXfGuTUla8T9Fs/l5+FMVityNoBlfyTu+fRQYMPbcg5fVoEfQ/8/euMbjGXF8rKzOC9qw63fSNeo3J26dymf9xvrUZbOB9dO4CGNSo7N5OAS6cc3ri0/Iwu8apXzZPoNq1dhaLiKEp5RMlqX8ZQ269ougeEOyBNxf5nTom3Zf2NkX1ZsMHTR7G0BbCCvr/hEsDjgQHAEcDVwB2AEZGleJNBYLox5i/Ho1Rh3XBMe1Zs2ceJ3aIvXfHgKV0cL4SbzO4Z0ZkODWtE3VfDyQQq1g1YuCT5vH4tqZSZweRFm5m8cHOMkSVepawMHj+re8Tl4u271j/C7BTxiqVfXODHGU1HeKfdMKQ9v6zcxqxV5QdH2OGz6wcyaeGmmB8/tFND/nPRYQzt2JDnpiyNeT3h+qsBfHPLIJZujlxvMJoD4FpVstl1oLDMbc+e04PuzWtbXoeKPWkb0rFhaY2+0nl+k7Qfb7xCbpWMMZONMQ8ZY4YDdYHuwHXAz8BRwGvAAhHZKiKfJyTaCurGYzrQsEYOh7e1voHP8jsabF6nKu9fM8DSPJ4XHl52R3nJEa05IYbEMVrRtuANbOcZdHKCjVPfBVMtJ4vLjmwTNr7nz+1ZrlyMz2k9mzLltqOD3hc3G7c5mRnC2X1blH5vQu2QKuZmLj6p0loai8m3DuJfZ3SL+nH/OP5QPrjmCAci8ujRoja3H98x5seLCMd3aWw58a6anclRHerz6oVly4X0axN+m9ysdpXSASkhIilz7bx+kcsUzRh1DOf2Lbvc6b2a065B5CnO7GLHd963r7GUuNr4I7NzZHXpKOAKumG01PHCePYWf3ov/wEQkcHAP4ATAGtFeVRQ3ZrX4te7w080bRf/I5loRrElmq8/VzIUnj21Z7OQk7mf2rNZyCnYkpFTG7IMOThtUkXle+uSvXRHNDo0qsHefO8pzIpaAzPC52WMISNDeCtgTuiWNp9qX/mv0AOMfE7p0ZRqOVmWDuaT3eBDG1qu1+dj7/Yp/pV1bVqT39bsjNgHMFVZSgBFpBrQH8/p4COAw4EawEI8LYEznQpQOeOhU7tQLcqO11bEsw9JVKHVZPstn96rGTOXb4u5TIVPNC/L7s/+m1sG8cvK7bzz8+qYpuBKRfGfFgr9Y4l2UEQ83Chw7M5zWl928q2DaGBz1xgrrzmwH2qkh8Tat7pbs1ph6/sl2SbSMWf2bs5Hv60Lef8bl/ZjWe6ecjOqVBThCkGfz8GErxuwC8/p3xnA08Cvxpj02NInsWQ5aA+2oZo+6hiW5+7l4v+zp5xEjcpZdG/ufoug3apWyuKl83vHvZ5qlawnDTUrH2xhsGNnfEijGhzSqAan92rGjn0FUTwySb7AAS4Z0Iovft/A4W3rleuPFa9I7/cPtw+mRuVsrnhzlq3Pmw56t6xtedkeLWrze4iSJx1SrIRLtAcjr154GF3DlJdKhl9lzPu2EI975uwebNmTX+a2J8/qzuNnhu4CUatKdlL1vbVbuD3G28A+4H/ARcYYa2XVVUIkWytWMM1qV2HjTvumkJ7/gPV5aMOpbkPrSrIk3j7Tbh8S15yXPvF+r6rnZMX8/ibTe9qndd3S01e+0z9HH5KY0iitLM55qsqadfex1Khs/bt3fr8W/L52p+2nelNB9ZwsujStyYINu8vcngy7FbtiCNyWnRGkwHtGhpCRFK/aHeF+LU/iGQV8KTDSWw5mpvcywxiTm4D4VIprU9+zM7NjhgE73HtSZ47rbG1GFEuSZNsRz2wj7jr4BgYrXxPJm5f144PZa+0MqJwGNXKY+c9jaBBlHcB08PiZ3Vi9LfZp9ewU7Snbc/q25Jy+ybFdCiWZDopShe9AON4uNekg3CjgO40xg4CawGDgQ6AV8AKwSUSWi8jbInK9iMR//kpVSPWq57Bq9IiYa7TFo32j8oMzLj+yjXN11FLABf09xaN7t6rtbiAB2jaIrdXr6EMa2HL6PJImtaokZIqzeFxyRGsg+qK08dToPKdvS+4YXn607hPe0j01o2iRc8rFA1rTo3ktzu4beQRu0rB4YOkbpex0IX6n+aYajXekswA5WZmsGj0iaRodklnELZoxpsgYM8sY84Ix5lxjTEugBTAKqA08Dzg3Z5AKKtZttq/ivd0jV6t6+58l06npJ8/qztsBI/vS3cD29Vk1egRNapWfBkyltuO7NGbV6BGWZvsIysYf79l9WrBq9IiQJUAiTVVop0Y1K/PZDUe6Om2h1c31gHbR1aO896TOTL51UNBp/VLJuX1bMOeeY5O6MkVFZPnwTERygL4cHBgyAPB1illjf2jKimi32Wf2bka/1nVtP2X45qX9mDB/A41qRr+RPb5LIyYusL84cdVKWaVHlkrZpbH3O55SLUouqRci0Xvs9G6cM+ZnOqboFFpWRdo+++6vV60Sc+4dFvX6szMz4h6w4sZp5qM61OfHpVv532X9+O6vXEQk6qkWVfzCjQJuysFk7wigJ1AJKALmAePwjAieYYxZ73Sg6cg3UrNFnapA8Bn4ov3xiogj/cVa1qvKdYPLz1drxQvn9WLn/rIjLX2tGHVjnN810UTs25BWpDpzFVHtqpWirm8WTDr073rktK5s31dQbv7Wni1r06dVHf55YieXIktt467ozwdzQpcvscLNszVvXtqPEmPIysxgUMDgquZ1Pa2Zh7VOtfnKU0+4FsB1eFqud+AZ+PEgningZhlj7BvaqULq1rwWYy46jEGHNOC9gI7uyXSqNV45WZk0qll2AMDII1pTt1o2p/YIXoDZbef1a8kPS7bQ2XvKYtrtQ1izPXxn+GsHt+PLPzZafo5knn7olQt68/WfsU/VpSJzo1ae3WpUzi5XYBk8v/kPr3VuJpFk4VSSf0T7+hzRPnXPboQbfduxcU2+/8dgWlWwvtrf/d2hGaPiEC4BvAJP657O9eui47o4OxVassrMEE7vlfiBI1YN79q4TCtQi7pVIw4uuXN4R+4M0mE+HpNvHRTT6Nl4nditSUxzT5eXBs1gIVSA/E5VUG63TvuqR1QkbRM4lZ9VIRNAY8z/JTIQpVT04un/07BGDrkBhVHdormQqohiTfJzvCPOszMT+8vQg5L04v4YfaUSZG4Mnawrsil/P5p9vnlgvZL5tLNymNvNPqrUtYPbU1hiuGhAq4Q+r34F0osmgClOBwxYV8emASWJLGHhpJqVs8tMCadSy70ndaZf6+inqZox6hjyi0pKr4fra9i/TcWdBiuRQiVWoW6vUinT9u4iSgVK6gRQRK4CzgN6AbWANsaYVQHL1MFTnPoU702fAzcaY3b6LdMNeAnoB2wH/gM8bOKpgOoyT0tNyoafsibeMqjCJIAqtV1+ZJuYHme1Ztzyx07U9uA4WX3/kuXUa7LEoRIjqRNAoCowCfgMeDbEMuOAlsBw7/XXgbeAkwFEpCYwGZiGp45hR+ANPPMcP+1U4Cp5/Oeiw6hayZ6BEtHOsqBUKMly+BnqODgzQ7OBeFn9iJPlu5DqcrISPyAulSV1AmiMeQ5ARPoEu19EOuFJ/I40xsz03nY18KOIHGqMWQxcgCeRvMRbvuZPEekI3CYiz6RyK6Cy5vg0HUkdC20BcJ7lViHv3/tP7kxP7ww+jtIP3zGh3lp9y+1VK8WnxEu05J7cMrIBwF48Bal9puNp3TvCb5kfA2oXTgSaAq0TEKNSSsWsa7Na9GqpRXFTWTI2M/gOKoLNI6x9y9NDUrcAWtAY2OLfimeMMSKS673Pt0xgyfTNfvet9L/D2+/wKoCWLXUyaRWb03o2ZdaqHW6HEbOmtSqzYVee22EoldKSuYHvvpO6cG7flhHrl6aCl87vxSbdXkUt4QmgiDwC3B1hsSHGmKkJCKccY8wYYAxAnz59kvYwqGaVbLbuTY4abqq8587tFfNjff1YKme710B/1aC2jBzYhtajvnT0edrWr86lA1tz0eGJLXeRSpKx9UilvkpZGXRtVivofalWDuqk7k3dDiEludEC+BzwdoRl1lhc1yaggYiIrxVQPDUNGnrv8y3TKOBxjfzuS0kfXTuAaUu3aqfXCujCw1uya38B1xzdLuHPnejNfkaGcP/JXRL8rKlB+4cpt+gp4PSQ8ATQGLMV2GrT6mYC1fH08/P1AxwAVPO7PhN4XEQqG2N8bcTDgA3AKpviSLhW9apxUb2KN12O8rQA3nbcoW6HoZSK0/FdGjNp4WYOaZR804AFk2otfyo+Sd0HUEQa4+mnd4j3ps4iUhtYY4zZboxZJCLfAP/x9t0DT42/Cd4RwOApE3M/MNZ7+vkQYBTwoI4Ajt/NQzuwv6Ao8oJKqTIibXzqVPXUm6yUlepj9dLXmYc156QeTfRMjUpKSZ0AAtfgSd58fB2SLgXGev8/H3gRz8he8BSCvsH3AGPMLhEZBrwMzAZ24Kn/94xjUaeRW4cdEnkhpVQpq6d2nzyrB0f9sYEezYP301KpwY3k74L+LRnaqWGZ29o3TI1WSJU4SZ0AGmMeAB6IsMwO4MIIy8wHBtkWmFIV1MiBbZiyKJcTuzUB4NPrB7J9nw42ckOtqtm2DY751xndWLl1X9D7qniLpDeskWPLc6no2X0q6tHTu5W5PueeY0s/53Ba1feMCL56UOL7H6vES+oEUCmVWG3qV2P6qGNKryekALFy3Hn9Qpe06ti4Jk//rQfHdg4cK6eclqged/WqW0vua1bOZtXoEQ5Ho5KFJoBKKZXmzjysudshKKUSTBNApZRSSikHnNm7OfWqV3I7jKA0AVRKKRdoEQKlKr6nz+7hdgghaX0BpZRKIC3wrJRKBpoAKqWUUkqlGU0AlVJKKaXSjCaASimllFJpRhNApZRSygVZmZ5dcKt6VV2ORKUjHQWslFJKuaBWlWz+e0kferWs43YoMXnlgt7szdO54FOVJoBKKZVAbep75mQd0rFhhCVVOhjaKXVnYPFNGalSkyaAKWLCjUfyy8rtboehlIpTm/rV+P2+46hZRTe/Sin36BYoRXRtVouuzWq5HYZSyga1qma7HYJSKs3pIBCllFJKqTSjCaBSSimlVJrRBFAppZRSKs1oAqiUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaUYTQKWUUkqpNKMJoFJKKaVUmtEEUCmllFIqzWgCqJRSSimVZjQBVEoppZRKM2KMcTuGpCUiW4DVbseRAPWBrW4H4SJ9/fr60/X1p/NrB339+vrT4/W3MsY0CLxRE0CFiMw2xvRxOw636OvX15+urz+dXzvo69fXn96vX08BK6WUUkqlGU0AlVJKKaXSjCaACmCM2wG4TF9/ekvn15/Orx309evrT2PaB1AppZRSKs1oC6BSSimlVJrRBFAppZRSKs1oAqhKichUETEBl/Fux5Vo4vG19/Wf5XY8iSIir4nIchE5ICJbROQzEenkdlyJICJ1ReRFEfnL+/rXisirIlLP7dgSRUSuEpHvRWSn97vf2u2YnCQi14nIShHJE5E5InKU2zElgogMEpHPRWS993Me6XZMiSQi/xSRWSKy27ud+0JEurodlxs0AVSB3gCa+F2udjccV/wdKHE7CBfMBkYCnYDjAQGmiEi2m0ElSFOgGXAH0A24EBgEvOtmUAlWFZgEPOByHI4TkXOA54HHgF7ADOBrEWnpamCJUR34E7gZOOByLG4YDLwCHAEcAxTh2c7VdTMoN+ggEFVKRKYCfxpjbnA7FreISF/gY+AwYDPwN2PMh+5G5Q4R6Q78DnQ0xix2O55EE5ETgQlAbWPMbrfjSRQR6QPMAtoYY1a5HI4jROQX4A9jzJV+ty0FPjTG/NO9yBJLRPYCNxhjxrodi1tEpDqwCzjNGPOF2/EkkrYAqkDnishWEVkgIk+JSA23A0oU72sdB1xljMl1Ox43iUg14FJgDbDK3WhcUxPIB/a7HYiyj4hUwnOANyngrkl4WoVUeqmBJxfa4XYgiaYJoPI3DrgAGAI8DJwJfORqRIn1b+AbY8zXbgfiFm+/qL3AXuAEYKgxJt/lsBJORGrj+Q28ZowpcjkcZa/6QCaeFn5/m4HGiQ9Huex5YB4w0+U4Ek4TwApORB4JMrAj8DIYwBgzxhgz0Rgz3xgzHjgHGCYivd18DfGw+vpF5CKgB3C72zHbKZrP3+sdPH2ijgaWAB+ISFUXQrdFDK/fd0roC2A9nj6BKSuW169UuhCRZ4AjgTONMcVux5No2gewghOR+niOeMNZY4wpd5pLRDKAAuACY8x7TsTnNKuvH0+n4IspO/gj03t9pjHmSGcidFacn38lPKdFrjHGvOVEfE6L9vV7k7+v8AyAOcEYs9fhEB0Vy+df0fsAer/X+4HzjDEf+N3+MtDVGHO0a8ElWDr3ARSRZ4FzgSHGmL/cjscNWW4HoJxljNkKbI3x4d3wJEEb7Ysosay+fhG5G3gq4Ob5wD+AzxwILSHi/PzFe8mxL6LEiub1e/uAfo3nNQ9P9eQP4v78KyRjTIGIzAGGAR/43TWM9OrykrZE5Hk8Z7jSNvkDTQCVl4i0w9P/7ys8O4zOwNPAXGC6i6ElhDFmPZ5TfqVEBGCtMWaFK0ElkIi0x9PncwqwBWgOjMIzCGKCi6ElhDf5m4Rn4MdpQDXvQBiA7caYArdiSxQRaYynD9wh3ps6e/tCrjHGbHctMGc8A7wlIr/i2b5dg6cU0L9djSoBvK3c7b1XM4CWItITz/d8jWuBJYi3pfciPL/zHd7vPcDeinDQFw09BawAEJEWwNtAVzx1otYCXwIPVsCNvyUiYkiTMjDez38MntGRtfF0iJ8GPJwOR8jefnDfh7h7iDFmasKCcYmIPADcH+SuSyviKUIRuQ5PH88meOri3WqMmeZuVM4L811/0xgzMqHBuMC7XQ/mQWPMA4mMxW2aACqllFJKpRkdBayUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaUYTQKWUUkqpNKMJoFJKKaVUmtEEUCmV9kTEWLgMFpFVIhI4ZaAb8R4mIjtEpKbF5c8SkcUikul0bEqp1KCFoJVSaU9EDve7WgX4DngEz2w4PguBdsA2t6fMEpGvgXnGmH9aXD4D+At4rCLO6qGUip4mgEop5cc7V+oeknQKNBHpACwBDjHGLI3icfcApxtjDnMsOKVUytBTwEopZVHgKWARGSsis0VkhIgsFJH9IvKliNQVkfYi8r2I7PMu0z1gXRkiMkpElolIvogsEZFLLIRxCfCHf/InItki8pSIrPGua4OIfCIilfwe9xHQW0S6xPs+KKVSnyaASikVn5bAQ8A9wFXAEcAYYLz3chaQBYwXEfF73Ivex4wBRgCfAP8nIidFeL6hwIyA2/4JXADcCwwDbgF2AaV9/owxi4AdwLHRvkClVMWT5XYASimV4uoCA4wxywG8LX23A5cYY/7nvU3w9CfsCCwSkfbAtXhOM7/pXc8UEWkC3A9MCPZE3vX0At4OuKsfMM5vXQDvB1nFH95llVJpTlsAlVIqPqt8yZ/XMu/f74Lc1sz7dyhQAnwiIlm+C/At0DPMaN06QA6wNeD2ecBIEblDRLoHtDT62wo0jviKlFIVnrYAKqVUfHYGXC8Icrvvtsrev/XxnJ7dFWKdTYB1QW73PT4/4PZH8CSU1wGPA+tF5EljzPMBy+X7rUMplcY0AVRKqcTbDhQBA/EkboFywzwOoLb/jcaYPOA+4D7vKOFrgOdEZLEx5hu/RWv7rUMplcb0FLBSSiXed3haAGsZY2YHuRQEe5A30VsDtAm1Yu/o4H/gae3rHHB3azwlZJRSaU5bAJVSKsGMMYtF5N94RgY/AczGc2q2C576fleEefh0oEwtPxH5BJgDzAUOcHDk8TS/ZarhGYRyr40vRSmVojQBVEopd1yPpzXuSjxlZHbjmW3kvxEe9zHwhohUMcYc8N42AzgHz+jjDO96zjTGzPZ73HHAfmCiba9AKZWydCYQpZRKId7izuuA640xH0TxuHeBfRFaF5VSaUITQKWUSjEicjtwqjHmSIvLtwAWA92NMcsiLa+Uqvj0FLBSSqWel4CqIlLLGBOqlIy/5sA1mvwppXy0BVAppZRSKs1oGRillFJKqTSjCaBSSimlVJrRBFAppZRSKs1oAqiUUkoplWY0AVRKKaWUSjP/D2eMXjWISRJiAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<Figure size 648x432 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(grid, sample[det_string], color='C0',label = 'strain')\n", + "plt.plot(grid, sample['h1_output_signal'], color='C1',label = 'signal')\n", + "#plt.xlim(-1.5, 0.5)\n", + "#plt.ylim(-150, 150)\n", + "plt.xticks(fontsize=14)\n", + "plt.yticks(fontsize=14)\n", + "plt.ylabel('Whitened Strain and Signal ({})'\n", + " .format(det_name), fontsize=15)\n", + "plt.xlabel('Time (s)', fontsize=15)\n", + "plt.legend(fontsize=15)\n", + "\n", + "# Adjust the size and spacing of the subplots\n", + "plt.gcf().set_size_inches(9, 6, forward=True)\n", + "plt.tight_layout(rect=[0, 0, 1, 0.9])\n", + "plt.subplots_adjust(wspace=0, hspace=0)\n", + "\n", + "# plt.show()\n", + "plt.savefig(plot_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.9.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/python3_samples/generate_sample.py b/gen_samples_py3/1-generate_sample.py similarity index 73% rename from python3_samples/generate_sample.py rename to gen_samples_py3/1-generate_sample.py index e343b0abdc12e79b37caf5c5c2584f17f005ba25..ef0e8d5ec9da380e81e5b6afe53b227f2cf67a66 100644 --- a/python3_samples/generate_sample.py +++ b/gen_samples_py3/1-generate_sample.py @@ -24,24 +24,8 @@ from utils.waveforms import WaveformParameterGenerator if __name__ == '__main__': - # ------------------------------------------------------------------------- - # Preliminaries - # ------------------------------------------------------------------------- - - # Disable output buffering ('flush' option is not available for Python 2) - #sys.stdout = os.fdopen(sys.stdout.fileno(), 'w', 0) - - # Start the stopwatch script_start = time.time() - - print('') - print('GENERATE A GW DATA SAMPLE FILE') - print('') - # ------------------------------------------------------------------------- - # Parse the command line arguments - # ------------------------------------------------------------------------- - # Set up the parser and add arguments parser = argparse.ArgumentParser(description='Generate a GW data sample.') parser.add_argument('--config-file', @@ -50,9 +34,7 @@ if __name__ == '__main__': default='default.json') # Parse the arguments that were passed when calling this script - print('Parsing command line arguments...', end=' ') command_line_arguments = vars(parser.parse_args()) - print('Done!') # ------------------------------------------------------------------------- # Read in JSON config file specifying the sample generation process @@ -89,8 +71,6 @@ if __name__ == '__main__': # Define some useful shortcuts random_seed = config['random_seed'] - max_runtime = config['max_runtime'] - bkg_data_dir = config['background_data_directory'] # ------------------------------------------------------------------------- # Construct a generator for sampling waveform parameters @@ -107,46 +87,14 @@ if __name__ == '__main__': waveform_parameters = \ (waveform_parameter_generator.draw() for _ in iter(int, 1)) - # ------------------------------------------------------------------------- - # Construct a generator for sampling valid noise times - # ------------------------------------------------------------------------- + print('Using synthetic noise! (background_data_directory = None)\n') - # If the 'background_data_directory' is None, we will use synthetic noise - if config['background_data_directory'] is None: - - print('Using synthetic noise! (background_data_directory = None)\n') - - # Create a iterator that returns a fake "event time", which we will - # use as a seed for the RNG to ensure the reproducibility of the - # generated synthetic noise. - # For the HDF file path that contains that time, we always yield - # None, so that we know that we need to generate synthetic noise. - noise_times = ((1000000000 + _, None) for _ in count()) - - # Otherwise, we set up a timeline object for the background noise, that - # is, we read in all HDF files in the raw_data_directory and figure out - # which parts of it are useable (i.e., have the right data quality and - # injection bits set as specified in the config file). - else: - - print('Using real noise from LIGO recordings! ' - '(background_data_directory = {})'.format(bkg_data_dir)) - print('Reading in raw data. This may take several minutes...', end=' ') - - # Create a timeline object by running over all HDF files once - noise_timeline = NoiseTimeline(background_data_directory=bkg_data_dir, - random_seed=random_seed) - - # Create a noise time generator so that can sample valid noise times - # simply by calling next(noise_time_generator) - delta_t = int(static_arguments['noise_interval_width'] / 2) - noise_times = (noise_timeline.sample(delta_t=delta_t, - dq_bits=config['dq_bits'], - inj_bits=config['inj_bits'], - return_paths=True) - for _ in iter(int, 1)) - - print('Done!\n') + # Create a iterator that returns a fake "event time", which we will + # use as a seed for the RNG to ensure the reproducibility of the + # generated synthetic noise. + # For the HDF file path that contains that time, we always yield + # None, so that we know that we need to generate synthetic noise. + noise_times = ((1000000000 + _, None) for _ in count()) # ------------------------------------------------------------------------- # Define a convenience function to generate arguments for the simulation @@ -170,23 +118,43 @@ if __name__ == '__main__': samples = dict(injection_samples=[], noise_samples=[]) injection_parameters = dict(injection_samples=[], noise_samples=[]) + for sample_type in ['injection_samples', 'noise_samples']: - print('Generating samples containing an injection...') - n_samples = config['n_injection_samples'] - arguments_generator = \ + # --------------------------------------------------------------------- + # Define some sample_type-specific shortcuts + # --------------------------------------------------------------------- + + if sample_type == 'injection_samples': + print('Generating samples containing an injection...') + n_samples = config['n_injection_samples'] + arguments_generator = \ (generate_arguments(injection=True) for _ in iter(int, 1)) - print('Number of samples:',n_samples) + + else: + print('Generating samples *not* containing an injection...') + n_samples = config['n_noise_samples'] + arguments_generator = \ + (generate_arguments(injection=False) for _ in iter(int, 1)) - sample_type = 'injection_samples' - for i in range(n_samples): - print(i) + # --------------------------------------------------------------------- + # If we do not need to generate any samples, skip ahead: + # --------------------------------------------------------------------- + + if n_samples == 0: + print('Done! (n_samples=0)\n') + continue results_list = [] - arguments = next(arguments_generator) - print(arguments) - result = generate_sample(**arguments) - results_list.append(result) - + # --------------------------------------------------------------------- + # Loop over injections or noise: + # --------------------------------------------------------------------- + for i in range(n_samples): + print('Generating the',i,'th samples\n') + arguments = next(arguments_generator) + print(arguments) + result = generate_sample(**arguments) + results_list.append(result) + # --------------------------------------------------------------------- # Process results in the results_list # --------------------------------------------------------------------- @@ -202,7 +170,10 @@ if __name__ == '__main__': injection_parameters[sample_type] = \ list([injection_parameters[sample_type][i] for i in idx]) - print('Sample generation completed!\n') + if sample_type == 'injection_samples': + print('Signal+noise generation completed!\n') + else: + print('Noise generation completed!\n') # ------------------------------------------------------------------------- # Compute the normalization parameters for this file @@ -248,8 +219,8 @@ if __name__ == '__main__': static_arguments=static_arguments) # Collect and add samples (with and without injection) - for sample_type in ('injection_samples', 'noise_samples'): - for key in ('event_time', 'h1_strain', 'l1_strain'): + for sample_type in ['injection_samples', 'noise_samples']: + for key in ['event_time', 'h1_strain', 'l1_strain']: if samples[sample_type]: value = np.array([_[key] for _ in list(samples[sample_type])]) else: @@ -258,10 +229,10 @@ if __name__ == '__main__': # Collect and add injection_parameters (ignore noise samples here, because # for those, the injection_parameters are always None) - other_keys = ['h1_signal', 'h1_snr', 'l1_signal', 'l1_snr', 'scale_factor'] - print(injection_parameters['injection_samples']) + #other_keys = ['h1_signal', 'h1_snr', 'l1_signal', 'l1_snr', 'scale_factor'] + other_keys = ['h1_signal', 'h1_output_signal','h1_snr', 'l1_signal','l1_output_signal', 'l1_snr', 'scale_factor'] + for key in list(variable_arguments + other_keys): - #print(key) if injection_parameters['injection_samples']: value = np.array([_[key] for _ in injection_parameters['injection_samples']]) diff --git a/gen_samples_py3/99-loopplot.py b/gen_samples_py3/99-loopplot.py new file mode 100644 index 0000000000000000000000000000000000000000..5381fdda79da6a7bccac0da8292326f638be9d77 --- /dev/null +++ b/gen_samples_py3/99-loopplot.py @@ -0,0 +1,10 @@ +import os + +for i in range(10): + print(i) + os.system('python plot_sample.py --hdf-file-path output/train.hdf --sample-id '+str(i)+' --plot-path "plots/'+str(i)+'.png"') + + +for i in range(100,110): + print(i) + os.system('python plot_sample.py --hdf-file-path output/train.hdf --sample-id '+str(i)+' --plot-path "plots/'+str(i)+'.png"') diff --git a/gen_samples_py3/99-whiten_wavepic_draw.ipynb b/gen_samples_py3/99-whiten_wavepic_draw.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..975ed1e1ecda9e5770312ef1bd8bc5c1c739eb12 --- /dev/null +++ b/gen_samples_py3/99-whiten_wavepic_draw.ipynb @@ -0,0 +1,347 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "import time\n", + "import numpy as np\n", + "import pandas as pd\n", + "import utils.samplefiles\n", + "\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "# Start the stopwatch\n", + "script_start_time = time.time()" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "hdf_file_path = './output/train.hdf'\n", + "plot_path = './plots/signalnoise_id0.png'\n", + "# which sample is used to shown\n", + "sample_id = 0" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "key: config_file values: default.json\n" + ] + } + ], + "source": [ + "data = utils.samplefiles.SampleFile()\n", + "data.read_hdf(hdf_file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "# when split_injections_noise is set to be True, GWs contain signals and pure waves are splited\n", + "df, noise = data.as_dataframe(injection_parameters=True, \n", + " static_arguments=True, \n", + " command_line_arguments=False, \n", + " split_injections_noise=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Index(['approximant', 'bandpass_lower', 'bandpass_upper', 'coa_phase', 'dec',\n", + " 'delta_f', 'delta_t', 'distance', 'domain', 'event_time', 'f_lower',\n", + " 'fd_length', 'h1_output_signal', 'h1_signal', 'h1_snr', 'h1_strain',\n", + " 'inclination', 'injection_snr', 'l1_output_signal', 'l1_signal',\n", + " 'l1_snr', 'l1_strain', 'mass1', 'mass2', 'noise_interval_width',\n", + " 'original_sampling_rate', 'polarization', 'ra', 'sample_length',\n", + " 'scale_factor', 'seconds_after_event', 'seconds_before_event', 'spin1z',\n", + " 'spin2z', 'target_sampling_rate', 'td_length', 'tukey_alpha',\n", + " 'waveform_length', 'whitening_max_filter_duration',\n", + " 'whitening_segment_duration'],\n", + " dtype='object')\n" + ] + } + ], + "source": [ + "print(df.columns)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "7.340279606636548\n" + ] + } + ], + "source": [ + "print(df.injection_snr[sample_id])" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "36.217808319315374\n", + "76.55000144869413\n", + "0.7305299539277823\n", + "0.5974611672286425\n", + "0.017164934231118905\n" + ] + } + ], + "source": [ + "print(df.mass1[sample_id])\n", + "print(df.mass2[sample_id])\n", + "print(df.spin1z[sample_id])\n", + "print(df.spin2z[sample_id])\n", + "print(df.scale_factor[sample_id])" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.9801424781769557\n", + "0.48680334688549004\n", + "3.776917009710014\n", + "0.821770111935331\n", + "4.448951217224888\n" + ] + } + ], + "source": [ + "print(df.coa_phase[sample_id])\n", + "print(df.inclination[sample_id])\n", + "print(df.ra[sample_id])\n", + "print(df.dec[sample_id])\n", + "print(df.polarization[sample_id])" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "sample = df.loc[sample_id]" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.5\n", + "2048.0\n", + "8.0\n" + ] + } + ], + "source": [ + "# Read out and construct some necessary values for plotting\n", + "seconds_before_event = float(sample['seconds_before_event'])\n", + "seconds_after_event = float(sample['seconds_after_event'])\n", + "target_sampling_rate = float(sample['target_sampling_rate'])\n", + "sample_length = float(sample['sample_length'])\n", + "print(seconds_before_event)\n", + "print(target_sampling_rate)\n", + "print(sample_length)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a grid on which the sample can be plotted so that the\n", + "# event_time is at position 0\n", + "grid = np.linspace(0 - seconds_before_event, 0 + seconds_after_event, int(target_sampling_rate * sample_length))\n", + "\n", + "# for time from -0.15s to 0.05s\n", + "#grid = np.linspace(0 - seconds_before_event, 0 + seconds_after_event, int(target_sampling_rate * sample_length)+1)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "det_name = 'H1'\n", + "det_string = 'h1_strain'" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "145.16235\n" + ] + } + ], + "source": [ + "maximum = np.max(sample[det_string])\n", + "print(maximum)" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "5.173865814288047e-21\n" + ] + } + ], + "source": [ + "maximum = max(np.max(sample['h1_signal']), np.max(sample['l1_signal']))\n", + "print(maximum)" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "23.74179629063074\n" + ] + } + ], + "source": [ + "maximum = max(np.max(sample['h1_output_signal']), np.max(sample['l1_output_signal']))\n", + "print(maximum)" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAF9CAYAAACUKdl5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAB7QklEQVR4nO3dd5gT5fbA8e/ZwtJ771XpTYqIIogoil2vvWDv9V6Va+8Xe9d70fsTr4rYGxaKiiigAoIiIL23pfft7++PJEs2mzJJZjLJ5nyeJ89uksnkpM2ceed9zyvGGJRSSimlVPrIcDsApZRSSimVWJoAKqWUUkqlGU0AlVJKKaXSjCaASimllFJpRhNApZRSSqk0owmgUkoppVSayXI7gGRWv35907p1a7fDUEoppZSKyZw5c7YaYxoE3q4JYBitW7dm9uzZboehlFJKKRUTEVkd7HY9BayUUkoplWY0AVRKKaWUSjOaACqllFJKpZmkTgBFZJCIfC4i60XEiMjIgPvHem/3v/wcsEyOiLwoIltFZJ93fc0T+kKUUkoppZJIUieAQHXgT+Bm4ECIZaYATfwuJwbc/xxwJnAecBRQE5ggIpkOxKuUUkoplfSSehSwMeYr4CvwtPaFWCzfGLMp2B0iUgu4HLjUGDPZe9tFwGrgWGCi3TErpZRS6WD37t3k5uZSWFjodihpKTs7m4YNG1KzZs2YHp/UCaBFR4pILrAT+AG42xiT673vMCAbmORb2BizVkQWAUegCaBSSikVtd27d7N582aaNWtGlSpVEBG3Q0orxhgOHDjA+vXrAWJKApP9FHAk3wAXA0OBvwP9gO9EJMd7f2OgGNga8LjN3vvKEZGrRGS2iMzesmWLM1ErpZRSKSw3N5dmzZpRtWpVTf5cICJUrVqVZs2akZubG/kBQaR0C6AxZrzf1fkiMgfP6d0RwMcxrnMMMAagT58+Ju4glVJKqQqmsLCQKlWquB1G2qtSpUrMp+BTvQWwDGPMBmAd0MF70yYgE6gfsGgj731KKaWUioG2/Lkvns+gQiWAIlIfaAZs9N40BygEhvkt0xzoBMxIeIBKKaWUUkkgqU8Bi0h1oL33agbQUkR6Atu9lweAj/AkfK2BfwG5wCcAxphdIvJf4AnvQJFtwDPAH3jKxyilVFoxxrBy6z7aNqjudihKKRclewtgH2Cu91IFeND7/0N4Bnd0Az4DlgBvAouBAcaYPX7ruAVPQvgeMB3YC5xsjClOzEtQSqnkMXbGKo55+gd+W7PD7VCUSmpjxozh008/tXWdY8eORUTYu3evreuNRVK3ABpjpgLhTnAfb2Ed+cCN3otSSqW1eWt3ArBm2356t6zjbjBKJbExY8bQtWtXTjvtNNvWOWLECGbOnEnVqlVtW2eskjoBVEoppZRKZgcOHLA8IrpBgwY0aNDA4YisSfZTwEoppZRSjliwYAHDhw+nbt26VKtWjU6dOvHyyy8zePBg5syZw5tvvomIICKMHTsWgNatW/P3v/+dhx9+mObNm5cWYZ45cyannHIKTZo0oVq1avTs2ZN33nmnzPMFngJetWoVIsL777/P1VdfTa1atWjevDn3338/JSUljr52bQFUSimlVFo6+eST6dSpE2+//TY5OTksXryY3bt388orr3DmmWfStm1b7r33XgDatWtX+rhx48bRpUsXXnnlFYqKigBYvXo1AwcO5JprrqFy5cpMnz6dSy+9lIyMDM4777ywcdxxxx2ceeaZfPjhh3z77bc89NBDdOnShbPPPtux164JoFJKKaXi9uAXC1i4Ybcrz925aU3uP7lLVI/ZunUrK1eu5LPPPqNbt24ADB06tPT+atWq0aBBAw4//PCgj58wYQKVK1cuvX7uueeW/m+MYdCgQaxbt47XXnstYgI4aNAgnn76aQCGDRvGN998w8cff+xoAqingJVSSimVdurWrUuLFi245ppreO+996KaUm3o0KFlkj+AHTt2cNNNN9GqVSuys7PJzs5mzJgxLFmyJOL6jjvuuDLXO3fuzLp16yzHEwttAVRKKaVU3KJtgXNbRkYGkyZN4u677+ayyy7jwIEDDBw4kBdeeIFevXqFfWyjRo3K3TZy5Eh+/vln7r33Xjp37kzNmjV59dVX+eyzzyLGUrt27TLXK1WqRF5eXlSvJ1qaACqllFIqLXXs2JGPPvqIwsJCfvzxR+68805GjBgRsfUtcAq2vLw8JkyYwMsvv8w111xTervTAznioaeAlVJKKZXWsrOzOeaYY7jtttvYuHEjO3fujKoVLj8/n5KSEnJyckpv27NnD59//rlTIcdNWwCVUkoplXb++OMP/vGPf3DOOefQtm1bduzYweOPP06PHj2oW7cuHTt2ZOLEiUycOJF69erRpk0b6tWrF3RdtWrVom/fvjz00EPUrFmTjIwMRo8eTa1atdi9252BMZFoC6BSSiml0k7jxo1p1KgRjz76KCeccALXXXcdnTp1Km21u+eee+jUqRNnn302ffv25Ysvvgi7vnHjxtG2bVsuvvhibr75Zs4880wuvvjiRLyUmIgxxu0YklafPn3M7Nmz3Q5DKaVsc/P4uXw2bwPPndOT03o1czsclaIWLVpEp06d3A5DEfmzEJE5xpg+gbdrC6BSSimlVJrRBFAppZRSKs1oAqiUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaUYTQKWUUkqpNKMJoFJKKaVUmtEEUCml0ohO/qSUNYMHD+ass85yNQYR4aWXXnJk3VmOrFUppVRSE3E7AqWS2yuvvEJ2drbbYThGE0CllFJKqQCdO3d2OwRH6SlgpZRSSqWlBQsWMHz4cOrWrUu1atXo1KkTL7/8MhD8FPAHH3xAhw4dqFKlCkOGDGHu3LmICGPHji1dpnXr1vzjH//g2WefpXnz5tSpU4dzzz2XnTt3li6zb98+brjhBg499FCqVq1KmzZtuP7669m9e3ciXjagLYBKKaWUSlMnn3wynTp14u233yYnJ4fFixeHTMJmz57Nueeey1lnncWLL77IokWLOOecc4Iu+/7779O9e3fGjBnDunXruO2227jrrrt45ZVXANi/fz/FxcU8+uijNGjQgLVr1/Loo4/yt7/9jYkTJzr2ev1pAqiUUilk5/4CRn/9Fw+c0oXK2Zkxr0cHgyjbfT0KNs1357kbd4MTRkf1kK1bt7Jy5Uo+++wzunXrBsDQoUNDLv/444/TqVMnxo8fj4gwfPhwCgsLufPOO8stm52dzaeffkpWlifNWrhwIePHjy9NABs0aMCrr75aunxRURFt2rThyCOPZM2aNbRs2TKq1xILPQWslFIp5KlJixk/ay0fzFkX0+N18IdSHnXr1qVFixZcc801vPfee+Tm5oZdftasWZx88smI34/olFNOCbrskCFDSpM/8PQnzM3NpbCwsPS2t956i169elG9enWys7M58sgjAViyZEk8L8sybQFUSqkUoi13KmlF2QLntoyMDCZNmsTdd9/NZZddxoEDBxg4cCAvvPACvXr1Krf8pk2baNCgQZnbAq/71K5du8z1SpUqYYwhPz+f7OxsPvnkEy6++GKuvfZaHnvsMerWrcvGjRs5/fTTycvLs+01hhNVC6CIVBKRZiLSTkTqOBWUUkoppZTTOnbsyEcffcTOnTuZMmUKeXl5jBgxgpKSknLLNm7cmC1btpS5LfC6VR988AH9+/fnlVde4YQTTqB///7UqZPYtCpiAigiXUTkcRGZA+wF1gBLgK0ikisin4rIhSJSxelglVJKKaXslp2dzTHHHMNtt93Gxo0by4zY9enbty9ffPEFxq8Z/vPPP4/p+Q4cOEBOTk6Z2955552Y1hWrkKeARWQg8AgwCJgF/AC8AGwF8oHaQGugD/As8KKIPAM8a4zZ62jUSqmUYIzh1R+Wc1bv5jSsWdntcJRSqtQff/zBP/7xD8455xzatm3Ljh07ePzxx+nRowd169Ytt/ydd95J//79Offcc7n00ktZtGgRr732GuA5nRyNYcOGcf311/Poo4/Sv39/vvrqK7799ltbXpdV4foAfown4bvIGBO2t7GIZALHArd4b3rYluiUUilt0cY9PPHNYqYu3sL7Vw9wOxyllCrVuHFjGjVqxKOPPsqGDRuoXbs2Q4YM4fHHHw+6fJ8+fXj33Xe56667+Oyzz+jTpw+vvvoqw4YNo2bNmlE999VXX82KFSt4/vnnycvLY9iwYYwbN47DDz/cjpdmSbgEsJUxxlJPRGNMMTARmCgiepivlAKguMRzqmR/QZHLkSilVFkNGzbkrbfeCnn/1KlTy9129tlnc/bZZ5def/vttwHo0aNH6W2rVq0q97iRI0cycuTI0uuZmZk89dRTPPXUU2WWMwGjvAKv2ylkAmg1+bPrcUoppZRSyezaa69l2LBh1KlTh99++41HHnmEESNG0KZNG7dDi1rcZWBEpB7QxRgzzYZ4lFJKKaWS0rZt27juuuvYtm0b9erV45xzzuGJJ55wO6yY2FEHcDDwPhB7SXqllFJKqST3/vvvux2CbXQmEKWUYwxatVgppZJRuDIwKyyuo6pNsShgzbb91KteiWo5OkmLqjgEnX9MKaWSSbgsoxkwF/gxwjoOAU6yLaI0N+jJ7zmsVR0+uvYIx57DGMOTExdzYrcmdG1Wy7HnUUopVXEZY8rMi6sSL55RwuESwHnAJmPM7eFWICJnogmgreas3uHo+otKDK9MXc5rP65g6aMnOvpcKnkt37KXoU//wJTbjqZ9w+puh6OUSiHZ2dkcOHCAqlX1JKCbDhw4QHZ2dkyPDdcH8Begv8X16CFACtJJ5dPbF79vAOBz71+llLKqYcOGrF+/nv379ztaq04FZ4xh//79rF+/noYNG8a0jnAtgA8B/7UQxEfoYBKllFIqbfhmvtiwYQOFhYUuR5OesrOzadSoUdSzkPiEKwS9Fc+8v0opFRNtGLBfRXhL84uKeWriYm4a2oEalWM7faXcV7NmzZiTD+U+bblTSjlO+4knj2RIyt+btZbXflzJC98udTsUpdJWuDIw0VQ7NMaYc2yIRymVaMmQESjL7Mql3UzKi4o937nCYv3uxeLm8XPp0LA6NxzTwe1QVAoL1wewQcB1AY7CUxpmj2MRKaWSzupt+1ixdR9DDo2ts7FSyj6fzfMM3NIEUMUj5ClgY8wQ/wtwLJ4k8KrA+7z3O0JEBonI5yKyXkSMiIwMuF9E5AER2SAiB0Rkqoh0CVimjoi8JSK7vJe3RKS2UzGrimfxpj1JPdJt6uJc/ly/K7YHW2gKOvrJqVz6xqzY1u+gtdv3U1ySvJ+LUkolq2j6ALq1la0O/AncDBwIcv8dwN+BG4G+QC4wWURq+C0zDugNDPdeegNvORhzUnjsq0UMf26a22G4btvefMb/uibmx/+8YhvHPzeNt35ebWNU9hr5xixOevEnt8MIadOuPF6ZuszWJHr9zgMc9cT3PDHxL9vWaafnpyyl9agvKSoucTuUkEZ99AetR33pdhhKKRck/SAQY8xXxpi7jDEfAmW2pOIpQX4LMNoY85Ex5k/gEqAGcL53mU54kr6rjDEzjTEzgauBk0Tk0AS+lIQbM20Ff20qf7Y+iRuyHHHDuLmM+ng+K7fui+nxq7d5HhdzC1sa833Vcvfk88Q3i1m+Za9t6966Jx+Amcu32bZOO/37h+UAFCRxAjh+1lpXn3/xpj22fidUeM9PWcpbM1e5HUYZBwqKufJ/s1m3Y7/boaSdpE8AI2gDNAYm+W4wxhwApgG+udQGAHuBGX6Pmw7s81smLYU683f9O79VqFaBrXs9iUJhEu+IXZXAIwI9W6v8zVyxjaFP/+B2GGnj2SlLuPezBW6HUcakhZuYvHAzj3+zOOZ17Msv4oHPF3CgoNjGyOyTuyePl75bmnTdiGJJAJPpFTT2/t0ccPtmv/saA1uM3zvv/T/Xb5lSInKViMwWkdlbtmxxIOTk9+X8jW6HYKtk+sLabXdeIXvzi9wOQ6moaFkgZad//7CcsTNW8WaStW763Pbe7zw1aQm/r0uus0jhysDMIvi+8y0RKddWa4zpZ2dgbjHGjAHGAPTp06ci5w5pw5f7V8R9TvcHJpGZEecr072xI0yFPvSIT5I1hKgUV+Q9tZCsA8L2FXgO0pMtvnAtgAuCXN4EZoW4zw2bvH8bBdzeyO++TUADb39BoLTvYEO/ZVQaqKh5TrJtVFRZYuHQ469Nu5mxLL6Jl/KLirl5/NxyfalWb9tHXmH5U2PpnoQdKCjmvVlrEnZabs22/XxVwc6uqNQWbiq4kQmMI1Yr8SRxw/AkpohIZTz1Cm/3LjMTz0jiARzsBzgAqEbZfoHKAmMMizfvoWPj1Jn+J833cyoFDH/uRwBWjR4RdrmCopKQ/Zx+XLKVz+ZtYG9eEf8d2Rfw/F6PfnIqRx/SgDcv85ykqagHQtH619eL+N/M1TSsWTkh9S1PeH4a+wqKI37G6SrZ+selg6QfBCIi1UWkp4j0xBNvS+/1lt6+fM8Bd4rIGSLSFRiLZ9DHOABjzCLgG+A/IjJARAYA/wEmGGNi73Wapt6Yvorhz/3Iryu3ux1KDHTPF5RueFPGuWNm8vHc9VE/7oclFbc/c0mJ4ZWpy9h1oNDyY777azMbdnqqiu3PT8zAgX1JOkDBbaJHJK4JmQCKyL0iUiualYnIMSJycvxhldEHz+wjc4EqwIPe/x/y3v8E8CzwMjAbaAIcZ4zxr39yPvA7MNF7+R24yOY408KfGzydWNdsT6Eh+5rfWHLyiz9xykv21RJ8bdoKFmxIrk7PieREXv3bmp32rzTFTVu6hSe+Wcz9n/1pafm5a3Zw2djZTFmU63BkqWdPXiH5RZqopotwU8H1BdaKyGfAB8BMY0yZw0gRyQa6AScA5+CZPu4SOwM0xkwlTNONtxXwAe8l1DI7gAvtjEslTmFxCYc/9i0PntqFk7o3jXk9eqAZgveNmW9zncNHv1pk6/pSlX7vyrPzPSko8pR32muxJW9nFC2F6abbA5Po2LgG39wyyO1QKqjkao0INxXcKXj61gnwLrBJRDaLyAIR+U1EVuCZE3gWnuTv/4B2xphJodap7LM7r5BvFwVWv7Em1UYn7thXwLZ9BTz4xcKYHp9ar1ZZpZ+rSmXFJYbXf1yRdLXrgk0eEOijOeuiOuVuRUX+PSfrMWDYPoDGmF+MMRfiGVV7Ip7+dlPwJH3vAdcAHY0x3Y0xzxljUui8YGq7Zfw8Ln9zNut3BpsdzxoroxMrEqde7czl29i8O8+htSfejOVbeX7K0rDL/Lh0C61Hfcmy3Mg7C6el17c4NXw0Zx25Feg34YQJf2zgkS8X8czk1OqKvnDDbv7+we/c8eHvtqzP6u/3s3nr+W3NDsBTWHl/wcH6p7oNiE24U8CljDF7Odh/TiWBVd5pzYKVd1BlOT267LzXfqZetUrMuXeYo8/jmID35/zXfgHg5mM7hHzIl394ylnMWrWD9g1rhFwuXjv3F7A0dy99W9d17DmUvbbtzefvH/xO5yY1+ermo9wOp7wkaWra723525OXWoXcD3j3ObneqRgT5ebx8wDPSPl+j36rp6ptYCkBVBVTqp0KjpeTo8227StwbN3p7ML//sKf63ez4rETyQhR8DpZv8VOxxXtt/meT+eTkYAOib66lFv2hk4QnDkmi22l6bYdrCisnKpW4WkCmIbS7dSvbt4jSOJRCn+u3w0EDzF5o06MaL/Xb/+8BoBTe8Y+kCoZaRkRd9iexOuGOuGSvg6gSjLeH+l/fljubhxR8G2odDeh7LJzfwE792urb6w0Z0td/p/dlj35XPx/v8b1W7DjuxApdywsLmH8r2socXnWpGQruaoJYIpz6wu1NHevO08cB93pKLv0fGgyPR+a7HYYZbi5b9mTVxiyP3JJiSkt1aIqltd+XMG0JVt4b9Zat0MJa8y0FYz6eD4fzEnuOBNNE8BUZXMyY4zh6UmLWb7F/sSupMTwr68WlVbeT7RU7+Nz47tzmbrYwaK1yXZY6rC12/c7+376uPS2unGg0+2BSQx/blrQ+x6asJBD7vla56xOIa1Hfcmu/dGXeVmyeU/Ug+7s2Pz8sDj8TDc7vH20dx9wd8BNsjVCaAKoAGjzz6948btlXPT6L7asb9XWffR5ZDLrdx5g7tqd/GfaCm55b54t645VIvs+Tl2cy9c2Tfz+xe8bGPnGLFvW5c+NvqDJkGse8/RUR97PUNzc6Ad7v+P9DEpKDON/XUNhcdlWvVXbglcBe/vn1Z7HJcOH74INOw/QetSXUT3mnV9W03rUl67OyrF4c/hBFoGf5g9LtnDcs9P4YM4654IKYsPOAyzcuDuhz1lRhJsKbouI5Fq9JDLoiqS4xNhSymXumh0c9cR37MmL/qjt+GcPHrkX2nSU/u6sNWzdW8DA0d+VHvkXlxg+nLOOi//vV1uewyo39jsj35jFte/8lvgnjkI8LaOpvC8vLE5s8Ke8OD2hz2dVrInpR7+tY9TH8+n36BRe/DZ8vchEsvqdTNQZAWMM89buZMbybUHv35fvaY0aP2strUd9SZFfQv3s5CWA+y1WwYT62iz3dgtauCG6ZCzeA6QDUew/U/1skN3CjQJ+GR2X47jr3pnDxAWbWTV6RFzreXrSEtZuP8C8tTs5qkODsMsG/ggiHenFa+XWg6eV//GBPcVDY5Fsze9OmLhgE2//vJq3Lu9v/UFxvDGJekuNKR+mkxsn3867V8s6ca/L6d9XovlmgNixv5CnJy/hxqGh60UmgpXv4IadB2hau4rjsfh7f/Za7vxoPsd1bhT0/ke+LDtVYkFxCVmZyXNSzun6qYmwJ6+QM16ZQcOaOW6HkpRCJoDGmAcSGEfamrggtuncDor9R2r1FODO/QXUqJxNZog6bFafZ87qHaW3GWMSVr4h0nZswYZdjHjhJz6/YSDdm9dOSEyxaD3qS5Y+egLZYXYSV781J+w6pi/biggc0a6+LTHZuYv4aelW3vp5Ff++8LCovhvBlvxp6Vaa1alCm/rVYorlfzNXc//nCxh7aV8GH9owpnW41dpQAfbbtvp07npueW8e7111OP3b1nP8+YwxPDdlKau2eYr1r9leMSfIWrJpDz2b1wKSt6Xo5xXbWZq7N6ZBi3+s28kpL03nxzuG0KJuVQeic1/yHG6oqDiVOgWud3deIT0fmsy/vloU1RP7J5fBdoT//mFFjBHa79tFnh4MkxfGm4zb7+EJZec/3hvHrAG79hdyweu/cP5rv5Q53RRLxuBE7j7yjV+ZuGAzRTZ0Q7jwv78w5KmpMT9+ibfVbu2O8AOX3GwluffTP8tcT4MG7pj4Djx9LbFO9339bc1Onv92KZ/N2+Do8zgt1Dfbd3B2oLCYN2eujus5/tq0m6Wb93qfL7nSyPHekc0/LAk/wCSUa96aww3jynYDSraDM8sJoIgMEJHXRWSaiPwaeHEyyHTw4rdLk3KU3G7v6Z6v/9zkucGmED+bt96eFVVw//1ppW3r8u8rY0i/guAqNfy8YltU88xGu0lyOtGIdTuejr/H4c/9yPMO9iO1eoBWVFxi+8HcNws2McE7ZWayFiu3lACKyDBgGtAcOBLYAuwFegD1gD9DP1pZ8fTkJUz44+AR49TFubQe9WXpnL+pxs3v+9fzN/L6j+VbGJP0N+iaZDvitoseXKS2c8f8zPuz1/HXpoODCT6du55bI1QRqIi/71T6jfp/Xm6JJYnbtjef9nd/zf9NX2X5MRt2HmDMtNSZDCEUqy2ADwHPA76RCvcaY44BDgEKgan2h5Z+8v2KpfpOH/y2ZkfQZZ3aLCTjRjTa13rtO7+V62Addv1RPsHm3Xmc9vJ0tiR4MnTHJOOHHgc7W03jkWynexIhr7CYKYvsKQrh393hlvfm8cnc8Il9xXq/7f9NRkqgA0Xzfm7alcf7s62Vfykqdq8oeLDXtGFnHgCfzI0c//Z9BWzfV8AVb87msa/+Yq2F/p3jf10TdZyJYjUB7Ax8DZTg2R9XAzDGrAYeAO52IjgV2WVjZ/PolwsjLxgjOzaqdm2Yg20SpyzczFKbRlla3eS+OWMV89bu5P3ZsVeV/+6vzbQe9SUrQhTezi8qZtwv9m44kqE1wWquaYwJefATjdajvuT9JJ+lwC7uf7rw2FeLuOuT+Ql9TicOX1Zt3cdF//2F/QWJK8PS6b5vGP31X1E95tMISbG/SAm0VcHeb6tTwc1cvo32d3/N7FXby93nRAJv9yp7PzyZ3g9PZq+3hI+V2pajPk7s7yEaVhPAPCDDeNpXNwLt/O7bjefUsIrTHR/+Ufp/pKZs349wzfb9vPZjdC0e/s8Tbw3Cmcu38cb02Fpc9uYXRfX8uUFa3K7432yGPRt8BgKfZCxn8MXvnr4h89buDHr/81OWJm5HGsP74/Rb+vnvGzjjlRnWYolw/2tBugOE4tuw2yEw2d25v4A1IYol28nNBt3VMb6+1dv2MWd1+aTACYHf3WDf5ce+WsSPS7cybcnWGNYf+4/j31HOsX7Le/PYvDuv9HpxieHP9btifn5/oQ4Y4/np/7TMM6Di5xXBayOC5/2Lph9lEm7ey/GvgJFMrCaAvwOHev//FviniAwTkaPxnB5O3hQ3xTmxMfcfnfbUxMVxreu8137mwS9ia4Hsev9Ejnz8u7ie3wrf9iFZO+IGs31f7JOrWzFj+baInc5Pe3k63/xZdjaTTbvy+P6vg6f4nHhHjYGVAX1f2971VbmBAaGeO9aY5q/bRdf7J/LlH36v2ca9y7HP/MCgJ7+3bX2JNm/tTlZu3RfzW+L/uMDP6Ognp3LmqzNjji0eS3P3Mmf1jqCJkzGG9S5NYWnVWzNXs9s7AcDzU5Zw0os/2ZYEWmXntvXm8fNod9dXrsdhRbS/hWTLVa0mgM9xMPa7gH3AROB7oCFwve2RqTICd4h22exQPzarP8Ote8MnOpt25YW9PxqhYrJ6arTIO4OEb2kr/T+enbyE3D32vQbwlCi58d25cfWleXbykoive97andwwbm6Z205/ZTqXjnVmGrVI2+7APkZ2b0z/3ODZaf64dEsUp6qt3VdYXBLxuw5woMCeqb82785jd15hTO9R7u48+jwymWW5ZbtWnPby9JhK6yTbcVdgPC98u5QzX53BSS/+xIxlZVv8XvxuGQNHf8fqbaG3v3+u30VJmBYru1uoAstVvfT9Mv7pPc0435v4zVi+1bFBhE5/nJ//bl/5HFu6MAW5Ldbv9F+bkqsovKUE0BjzlTHmZe//64HD8LQI9gTaG2PCV59VUfP/0n01fyNDnprKFLvr1CXZhjnQp3PXc/i/vmX2qoPN54c9PLl0JgKrLG8EIvyqPw7oQzPeQt+y579dyj8+OHjK/dhnfgg7eiyvsJi7P5kf9jXePH4eX/y+wb4ZJqLYmm0MSMiT4Yi2XPRxZhxOnFJ61OKgpIKi2JL6XQcKy/RH6v/Ytwx9+oewj1m1NfgBzDcLNrF1bwFvzoivxluyiKbv69odZd8T3xyzgd97nzmrd3DSiz/xaphTt3H1vQ3y0Cv/N7vcbdv2lj2Q//g3z7bqyxDzkXe9fyKjv/YMYvBv0Y/03In04Zx1UZ9OfzLOM1r+otmKjJ2+0lKDQGDtTrfFVAjaeCw1xvxhjHH2XFWaE4QF3paJRAyzj3Ra0KZcypJZ3o7Ci/1e97Z9BfyxbmdM64tmAIJ/HxSr26Cte/N54PMFFAa0zPn3c1yWu5fHvgrd0fuj39bxzi9rDtZdtFGsiU2wvlmxfr6LbJi0fX9BUVQ7Bsvf2SC3bYtwKn5WkM7swcwN0dfTDtv3FdDjwUk8P6VsPbVIo9SnL4++f5tdnMgtIq3T9/nameBv8J4eXhjme71kc/SzUDjdaro3v4h//7CcoU//4FiLfjC+9/6pSUtC3ufzjw9+5xsHtoN2MngmS3jgi4WcO+Znt8OJWjSFoJuKyFUi8pCIPBFwedzJINPR9GWeTrKhjh4D+zr8tMy9jXmyi3QEHrjh+b/pq2h311fljqojeeDzBYydsSquGUUi9X12ZMcQYY9oZ9+sm8fPi+vx63bsp/N9E3nr59CtU/G+Rf7fl+emhC9Se04SbPTv8p7+C5eEBOVAFhbv9zOaUa2Bz5lfVBz3oLYnJy6h9agv3W78SpiCMN1IQs8E4kwswUT9nS7H+if55/rdnPrydEvL+r8FxvsW7smL7sxUMrBaCPp0YAXwMnA58LcgF2WjrVEmH7FyaoRsNFXt/UexWWEMTFywiVemLrM1Jt+9H83x9DXznfbx3+DNXbMjZL404Y/gp1vsFMvHtXl3Hh8ElKwRsfYZJduO0DfK1ImWAd/nbNdPIpbVxHK68JsF0b0Xi+Psh+RkOaH5foMXok00pi/bRvcHJpVeLykx5VpBv/lzE5e/Wf4Uqo9vuxv4HUiG3jIzl2/j49+C16oThFVb98U8CtsJizbu5s0Zq0qv780vivqb8+J30W3jA1n5LRcUHzxo+N2vpX6bhT67ZZ4rqqWTQ5bF5R4DJgEjjTGJGauvgMRPD2TX0Z3/TiLSD+OeT//ktYv7WF73Nws2RVUjz+oO3cpiX/6xkaxMd6fQjjZpv+T/fuWvTXuYcOORpbc58a3adaCQG8b9xlN/60GjmpVjXo/BRJWEbbdYgyySVJyK64HPF5S7LVIrsm19R+Ng5Z2OJRH3tWjl7s7j/dlreWrSEgYf2qD0/icnRldnrzQWK8sYw6xVOxxLBM57LXxr8+AYBujYJdg26YTnfwTgkiNak7s7j36PfUtGmA8+1n1PvPusUIXjrRxYJWOJsWhY3ZO1AF7Q5M9dvu9aQVGJpQKUkVcY38NfsGkOx0i/38Awt8Y4cjnUhiLZRilaFSlhySsspqTElLZqBNbW8iXpL8RwlB3s6/fRnHX8uHQrr061VsvsuSlLePn7ZUxZuJkZcfRHW7vd3jIdqbRJH+vXwuLz/WJ7ZuJIVZt2eZINXz+zjd6ZHjbsymP5FntHxfp/V8bPWsvZ/5nJ1yEGXkTD13IZqk6oVcmSoPjOpoQ7OAkVauDgisBW3bAj8SNcn7VqO1/Nj/5sQrByM6m4G7HaAjgDz6jfKQ7GokIIPEI55J6vbVt3tHWT/Jd+ZnL5jrwHl4v/5xBNaF/P30jNKtkMbF+/3H2RNoGR+nmlogMFxXS67xuuPrqto88Tz6cc+L5nZ3rWVlgc305rVZiSHWGJ/7/xfX9HffRHVMVs7RTPPn/igk3c91n5VkWr/N813446mnjK1AoUz4CfhydEV2c01OAXqwcm3kjKXIv0bRAondVn3Q77Dkjen72O47o0tm19dgj22wj3Ec9etZ09edYKrPu63/j7YM46bht2SOn168f9Zmldwfznh+VcPagtIkJJieFv/46tf3OwxDo5Uu3oWG0BvA24SkQu8Q4GqRp4cTLIdFZQXEJeoXNzJwZ+kePZ7eUVFrt2xHntO79xweu/BL3Pt0OwshF3wq79heXKe/i/T5MXbo5y5+QxPczAH9+MFsE2qGB9YxXp8/zAu/65a3aU1tGLV9f7J5aZFztaO/fH1xnbjq+wlRJBwfR79FvL02rF69cgI5ivfyf2nWsksbS0/2/mat799eB7Gaz2pTMt+NGt1HDwe+PGGQUnnjOa34H/5AKBzvr3TMsjjf/+we8Rl9kVx+97x/7C0sLell9emDdCPJ2pU5bVBPAPoBvwBrAW2BPkohzgP21bou3cX2C55l7u7jw63vtNaWul/wYp2O8n2lp+sYqldM66HeE7UkfqBB+4PVi8eQ/Xvl22VKZvAnIRT12vx7+Jvm9SsB14ODMDpl+yo5XWN8XR6a/MKK09ZocDEeZg3WfjlG0+ybAdLygu4ecV22k96kvuTvCculZFmyAXxdESGvhcxbFMWxhD28yURdZG8gf7zvy6Mnl6Sll9u16dujzo4KCVW/dFXX82kY0Ay7fsZdqSLWFisXZbrPxf65686MpTJQOrCeBlwKXey2UhLipJFRWXlKtNZ8VRj3/PnR95EtBIR02++7+wOBLWv7Dq0ty9fPPnRo579ocyp83e/tkz0COekWBlTgdZ2MN/8+dGdls8XRGNbwOKrfqSt3Lbiyg2IMu3WKsv5ltltBPNO8WO+XZ/X7cr7o3tu7+uCXqaNljCUFRcwjOTFgc9cHn5+2URZ2WJvj0J3olioFO8tu7Nj2kbEQ0rH9f+CIl/LOuNpQ5fICtdZXxhxLv98J83dsqizTFNRRdt157Hv/kraBHlIU9N5YqAwtNOtXDGkqgPffoHHoqyi0C0gkUV6v21s3tWIljqA2iMGetwHMqCWHd37e/2fClXjR5R9o4IP+Q9+UWlc0xGcvorM6JZdRnb9xXw9/d/Z19BMQcKi6meY7VranSstHhd83bZU2C78wrL7PTjSTnimbrNZ29+Uel0Qivi7NDuZDmPcDbsPMAhjWok7PlCJYr//Hg++YXFjBzYBvDbqAdZfOKCzbzw3TI2787n8bO6l7nvyYmLqVO1Euf3bxkyhlCd+V//cQWdm9SM/CIc1ueRKZzWs6mlZZ08xXkgzjp+Tlm1bR/92tQtd7v/V8Wuxp8zXy27LR33y2puP75j2Mc48ZlMXhhbqSW75+MVytbYs2O7dcHrP9OnVfnPM1oFxSVlpkyMt/9yojmzp1VJy44kJJzf1+4sV4w1VUbZBovTv64YePqXfRSiFlckkU4NTV2cW6YOWjBHPfF9TM8dKxEJumdz4jP1JOie5wq1GS07SCD+IHb6Jffh1uZrHcsrCp6gxFqA+JEgU8Q5dYoqksDRkJ/MXc/Dp3UFoPWoLy2tw+6df7K448M/GDNtBTv3FzL7nmPL3f/lHxsZ0b2JC5GF9/TkJfy2ZgdvXNov6scu2pg8PbtK/HZb0bTohkoWpy/bVjrZQqDhz03jq5uOKr3+4BcLqVUlmyODDDD858fzmbtmp+V4ko3VQtArRWRFiMsyEflNRN4QkcOcDlhZE6pQbrT7llj6iS3YULbfXaR9UNn+gp6FF26If8qwES/8GPc6AsWa/AERm0VHvjGL92fHsX4HhEogwn2mY2esspww2M2OAupOH8Pvziuk7T+/tFSu5aZ35/Ly9/EVw42VHafqK5JluXvDfr++jKIQ/Kqt1lvvww2wsOL7xVvYHmFKw6ACtldFxSVlTk+Hkmr94AL9tWlPmQO94hLDzePnlTkV7ntrUjn5A+t9AD/C01pYA/gFmOD9WxPIBmYDhwM/i8jxDsSpgB1RjA68JmDQgVWBR/CJOqAPfN4TbUjeFmzYnZLFfRNBiG8QyHt+M4ucHWMphXCiHmhQXMIb01fS55HYK1VZ+a5HimvTrryIO8Alm/ZQYuClEH1b/R/9+e8beHLiYnKjnC3HSWHrriV4579o4+64Roy7xWpZFIitrEzgV/nyNyOPwo1Ui/OpSUs489UZ/BnhLEWskj1v3BBDX8xQ/li3kx+Xhh68kihWTwHnAkuAk4wxpVsiEakCfAGsAboCnwMPAhNtjlMBb0xfxYC29eJaR2J+ZLE/yYGCYtfqKdn13qTCWTBD2SmQ4hHtaGQnvDtrLQ9+YV9n8MDP0MpnOnfNDk5/ZQZPBPQRtEO/x77llmM72L7eQIGvc+32/SzLjX8QRazCve879xdwwvM/UrVSZuICChDrqX+7+99GOpizkkTuLwjouhNw/yLvvLyRWtkLSwz5IbpKxCKe7amd+7vNu+2bnvWUlzxzDk8fdQzNalexbb3RstoCeBPwjH/yB2CMOQA8C1xvjCkGXsNTLkY5ZLaFJng7+f/2Ji7YxCcWJmvPLyqJ+lfrW3rwU1PL9bsrt2yYVYdrgVi9bR+tR30Zci7U522a2SQVzF2zk5e/j1x70KVaxiGF+uz32jxyO9odhwgs9fZNmuVQGZD/zVztyHp9CoL0Dz7qie8t13Czau32/TxgQ7LuS1gCE5dEcrNMl7+fAmqC2nEQGus6xv2yhkPv+Sb+AGwQ6+Yr4tkRmw7ybx0/z54VxchqAlgbaBTivkZAde//u4DkHMZVQVgp1RCuxl60P+pVfpOLh5ozMdCHAcWHoznatbIxD7dz/u6v0P2qvvb2i/x4rn397HbE0rcmgbZanNA8VMFoJ9z9yXzHByMF822IkhqtR33JFW/ODvIIj10HCvnNe+BlZXBKJB9HOIhK9lNh4VgZBBLtzB6B/vbvGXS+L3SCkQot8KnMjYE+blUscML5fnM6x1LX0k5WTwFPAJ4QkV3ABGNMgYhUAk4BnvDeD57Wv+inNFCWvTF9VcRl+jwy2fL6nPj6lZSYMv1EChPYR2dbQEIW9DSnjS+618PW3+tk7o9opQK/XWat2sHcYGVRHH57Ln9zNrWqZAe9b8qizRzfJfgx7sg3fo2qs3ekr9c4b32/bVEOWImpI78Dwr2+cAdg4KkSMCnKwsKBZq0KfxbkpBd/imv9TnPjoDGWn1ait1ehvlfxxJGMA1JmLA8++tgNVlsAr8EzH/CHwAER2QkcAN4HpgPXepfbANxlc4wqCoKErEVkjClTZytwejI7+e8IIp26ivY3uiiK2T0uG1t+5JbPlj35nPby9OieXNki6Cbd+P8b+UuR7/ddttpCEK51/PYQp/Pmr4vc6V08o2qi4t+67s/O1o7HvkqO4t/geV2n6u+ttLh+sgv8HiZfKuWcdGlFtloIeidwqoh0AfoAjYFNwGxjzAK/5T50IkhlXbidxytTl5er9u7E9zwwglg7SoeydnuYDs0WtlK+Rcb/uiZkgd74WH9XndrQfBxPuZokEOqgwP/2x74qX0Mv6GNieH7/eZb9pzIL1aLw4BcL6djYU9y6zOwzLvu/6da6bURjaYwza+zPj347EPh2+/ffdWsn/fqPK7jiqLYxP95qcX1/4eb9ThSnp+9Mppwrls8oFm63UFptAQTAGLPAGPOmMeZx798FkR+lksXncdaTsvoDjbZuVSytJ8ks1I4pWHLu1O//X37TviXjkXuknfeEEDXVLvzvL6X/r9kefs7meKwIUact3MGMb4aW39ftjOu53dwnWDndFlgL81MLA8Mglu9h+Vh8oyfdFKx4dzQOFBSzc390CcYFr/8SeSE/dnyHAr8LvztysHxQvF0DfJ6auDjuxOrROD5jO2rYJkrIFkAR6QwsN8bke/8Pyxjj7IR8SlkQrgU019sy49s4PD15Sdh1/fPj+THFcPVbsdVgVAdZaW3Y4bcTtTNpCpcCRerjlo5ueW+eQ2s2FfJU3O/rdnHx//3q6HMs3Fg2CcmNoVXa/723UrTcSdH0A3zp+2XUqVaJy49sE/PzHYgwGDHcNJxTFm2mc1P3p3e0Itwp4D/xFHf+1ft/6D6anvvcK8akLFm8uXz5k8APNdzE40tjrAdmZdRYIrbzn/8eXctkpGnZohVsI5aEfZQTwu7ivckyQMIOFfUrYeU37v97GDNtBRMXhG4VStffTqL4b7YvfcPeUkDRirZf7MMTFoZMAK3UtbSrNTLZhUsAhwAL/f5XaS7WnWy4pBKI+nRIOHd+FLnVLpl2HJuSaIaHRDr/tSCntOI4CghXoijc0Xowy7aE3kFY+erY+X1OtGC1AN0SLvkDkmImhXSU6FbZZ6cs4d8/xFZcJNi23u66loE++m0dmRnW3iS3d0UhE0BjzA/B/ldKxWbsjPJJSuCgnHTm5Kj0aPwRYdTvN39an/M13bW/66vS/+1OHAJLPil7hTrt6kY5qwM2DyR00upt+1Nmu261DmAZInIy0BHPSOBPjTHBp1ZQSW9/Gk747tZRV6QWDZUarnn7N0fXn8i6mU7zH0H9hIWdotstIuqgwNlFUpFz/VMrhpCjgEXkThH5MeC2bBGZBnwKPA68CfwpIk0djTIMEXlAREzAZZPf/eJdZoOIHBCRqd5yNhVStKc3v1+cnqdRkn0GD5UYyThi7+5PYxt8lOx8BbDtksyF1VVkYx0oUaSiE64MzOl4ijz7uwk4EngEqImnJmAxcLcj0Vm3GGjid/Gfj/gO4O/AjUBfIBeYLCI1Eh1kItw0fq7bISQ9Y6KbwUNVXMnY2TuvsOK0ADqpIo4QTid2zAcdztn/meno+iuCcAlgO+DngNvOBVYaY+43xuw1xvwGjAaGORWgRUXGmE1+ly3gaf0DbgFGG2M+Msb8CVwC1ADOdy9c56RyB3Slkp3mHM7avi95imir1Pbryu1uh5D0wiWAVYGdvisiUh3oBUwJWO4voJntkUWnrfcU70oRGS8ivjLtbfDMWjLJt6Ax5gAwDTjChThVEnBidgSVHkq0k5qjpi+zPk+qm6P5Z1SA/nGxyi9KnQEZKrxwCeAKoJ/f9WF4DoADE8BagJsdaX4BRgLDgSvxJHwzRKSe93+AwPM8m/3uK0NErhKR2SIye8uW9Owfp1RFc/9nf7odgrKZm6eAz49yZo6K5LkpS90OQdkk3CjgN4AHRKQIT8L0ILAF+CpguSF4+uC5whjztf91EfkZT/J6CeVPYVtZ3xhgDECfPn30eF+lPLfnm0wGb85c7XYIqoKZtUpPMarUFq4F8AXgXeBfwFtANnCeMaa0qqqI1MKTaAUmha4xxuwFFgAd8JSpAWgUsFgjv/uUUkqlmP0Rputy2qZd6VnEXdnH7WPzkAmgMabIGHM1UBtoaIxpZYz5PmCxfcAhwLPOhRgdEamMp0bhRmAlnkRvWMD9RwEzXAlQqQTboQODVAX0wrd6KlKltjyXC1xHLATtHTQRdC4vY0wRYL3XrgNE5CngC2AN0BC4F6gGvGmMMSLyHHCXiPwFLAHuAfYC49yJWCmllFLp7q9N7s6hEdNMIEmmOZ5T1fXx9FH8GTjcGOPr9PMEUAV4GaiDZ9DIcTp7iVJKKaXSVcongMaYcyPcb4AHvBellFIqbt8vznU7BKXiEm4QiFJKKaWC+Pi39W6HoFRcNAFUSimllEozmgAqpZRSSqWZkH0AReS7aFZkjDkm/nCUUkoppZTTwrUAbgu4HIKnfl5VPGVUqgJH4im4nL4TIyqllFJKpZiQLYDGmL/5/heRy4FDgSOMMWv8bm8JTAAmOxmkUkoppZSyj9U+gHcD9/knfwDe6w8Ad9kcl1JKKaWUcojVBLAxkBPivkp4ZuBQSimllFIpwGoCOBV4XET6+N8oIn2Bx4EfbI5LKaWUUko5xGoCeBWwHfhFRDaIyDwR2YBn2rXt3vuVUkoppVQKsDQVnDFmHdBbRE4E+uI5JbwJmGWM+crB+JRSSimllM2imgvYm+xpwqeUUkoplcKiSgBFJAdoBlQOvM8Ys9CuoJRSSimllHMsJYAi0hQYA5wQ7G7AAJk2xqWUUkoppRxitQXwdaA3cBuwEChwLCKllFJKKeUoqwngQOBKY8z7TgajlFJKKaWcZ7UMTC5wwMlAlFJKKaVUYlhNAO8D7hSRmk4Go5RSSimlnGf1FPAZQEtgtYjMAnYG3G+MMefYGZhSSimllHKG1QSwPrDc+3820MCZcJRSSimllNOszgQyxOlAlFJKKaVUYljtA6iUUkoppSoIyzOBiEgN4FTgEILPBHKHjXEppZRSjukvi6gru/m6pL/boSjlCqszgbQDZgBVgGrAFqCu9/E7gF2AJoBKKaVSwns5DwPQOm+cy5Eo5Q6rp4CfBWYBjfBM/XYinmTwQmAvoCOAlVJKKaVShNVTwP2AK4B87/VKxphiYJyI1AeeB45wID6llFLKVrXZ43YISrnOagtgZWC3MaYE2A409bvvT6CH3YEppZRSTjgv83u3Q1DKdVYTwCVAK+//c4FrRKSyiGQDlwMbnAhOKaWUslsx4nYISrnO6ing8UBP4C3gXmAisBsoATKBkQ7EppRSSimlHGC1EPQzfv//LCJdgeF4BoJ8Z4z506H4lFJKKVtp+59SUdQB9GeMWQu8ZnMsSimllOOyKHY7BKVcpzOBKKWUShv12cXt2e+XXq9CnovRKOUeTQCVUkqljcayrcz1q7MmuBSJUu7SBFAppVTaKAnY7WVT5FIkSrlLE0CllFJpIYcCLs/6usxtGRiXolHKXZoAKqWUSgvXZX3OmZk/lrlNRwSrdGV5FLCINAVOAprjmRnEnzHG3GlnYEoppZSdqpTOZqqUspQAisjpwLt4ij7nAgUBixhAE0CllFJJqQnbaBIwAMRDTwGr9GS1BfAxYBIw0hiz3cF4lFJKKdvNrHyj2yEolVSsJoAtgBs1+VNKKVWR1Gav2yEo5Qqrg0BmAIc6GYhSSillt5rs47as90Pef27W1MQFo1QSsdoCeBvwjojsBSYDOwMXMMbstzEupZRSKmZN2cpLlV5ghWnKWZnT3A5HqaRjNQH8w/v3DUL3mM2MPxyllFIqem9kP86QzN/pkPc/Csmib8Zf9M5YRm+WuR2aUknJagJ4GTpUSimlVJIakvk7AG1lA4tNSzIpcTkipZKbpQTQGDPW4TiUUkqpmOT4VSarI3vBQFXRmn9KhWO5ELRSSikVqLnkUmCyyaWOK89fn100koMFKmqw3/v3gCvxKJUqQiaAIvIrnrp/C0VkFhFOARtj+tkdnJ1E5DrgdqAJsAC4xRjzY/hHKaWUAmgn69lnKrOJeqW35VDATzm3ANA173X2UjXu58mkmCyKyadSmdsrUUh9drGB+n63Gr7JuZP6srv0lnqymxMzfqaG6LhEpcIJ1wK4AEoPoRaQwn0AReQc4HngOuAn79+vRaSzMWaNq8EppVTSM3ybczsA/fJeLm3t652xtHSJ+7Le4o6iq2Nau1BCC9nCGtOIl7Jf4ITMWbTOG8ftWePZamrxRvEJLKl8CQDnF9zFcRmz+b6kFwtKWpdJ/gBGZb1LbdnH8pImUT2/QTgj40cmlfSxJZFVKtmFTACNMZf6/T8yIdE45zZgrDHmNe/1G0VkOHAt8E/3wlJKqcTxJDq+8q+GHApLW9pqso+uGSv5paQTguGCzG8xwP+Kj2dQxh+l6/gh51Y65Y8lgxLerfRo6e1nZ/3AnUVXYsigMvmckzmV5aYpLSWXx7L/y1UFt/JzSSeezX6VaSXd+aj4KI7M+JPvS3qyuPJIAO4svJITMmcBcF3mZ1yf9TkA52ROLX2eazK/YFDmfEYyiZPyHyn3GmvLPgDaZWy0/L5ckzmB27I+IFuKAfh30UkcnfE7txReD8DEnFEAbDU1+bOkDRtMPTaZutyW/SEAE4oPp2fGMsYVDeXN4uO4OusLJhb3ZbupyduVHqNdxkbeLzqax4vOpU/GEjaZOuynMq1kM6dkzuCnkq60lU08V3QGed7P4/CMRfxV0oIShGIyKSKTARkLmFbSg8NkCQeoRI+MFUwuPoz7s9/kwcJL2E4NqpLPPipTTAYglt8DlX7EmJRt2LNERCoB+4HzjDEf+N3+MtDVGHN0qMf26dPHzJ4929H4/nX3NdSWyJXojcUfspVP0/q6Ii9n9dtj7TltfI0mOV+jnZ+j9fcr0a/RznVV7PdLMDSSHTSXLaw39alGHnuoSglCXfaQgWEjdekkaxAMy01TarGP3VTj8IwF9MxYwbiiIXTPWElT2Upd77ZknanP3JL2nJz5c+lzzStpR8+M5UHjeL/oaM7O+qHc/wDvFA3lgqxvyyzfLu8t2sv60sTojaLjuTRrIg8XXsiMki58nWP9uHpRSUs6ZcR2IubKgtt4rdIzpdf3mxwd/GGjpSXNWGRacljGEpp551FeVtKUtrKRfLIZW3w8lSiiCvnspQqNZTs7TXX2k0NnWc06U5+9VKUShQAUkA1AIZlkYsihgAJvO1QelciimExKEO8vp5BMDEIWxYT7/aZqFrPR1OWhR19w/HlEZI4xpk/g7ZYHgYhIa+BC4BCgcuD9xpiz4wnQQfXx1CjcHHD7ZuDYwIVF5CrgKoCWLVs6HtzZmVNpLlsjLGV19xZ5OavHg9bWZd9uN0NS9SesVHxKjFCC51LJ2wJVYoR8ssmhkAwxfonNHLab6tRmX+lv5vys78uts7lspXlm2e1KqOQPKJPw+f8PlEn+RhVewejs12ktm+iesQKAY/KfYq+pwqVZE7k3++3oXjzEnPwBtJKym/U9VKEqmgDapUPGejqwvsxt7TM2AFCFAq7N+oK9pjKVKESAbCmmyGSQgSn9fvon5QUmkyKyyKEAg7CfHLIpRjBUoogiMigms3Tf4kn8oChMmeHw+xcTcQk3zTdtXH1+SwmgiBwGTAPW4EkA/wBqAa2BdVBxKm0aY8YAY8DTAuj08w0teNrpp0hDFT9htr4u+77Cqfwa7fwc7WxvFAz7qVzaCgJCpnenV0wm1ThAY9nOctOURuwgn2x2UoNMimkmW1lrGlCNPPZSlbrspnvGCmaWdKaDrKMSRfxmOtBJ1tAv4y++Ke7LAXLIoYDuGSs4KfNnXik6lXyyuTHzEz4sPpqVpjHTcm7hl5JOXFV4G6dn/sTo7NcBGJj3PLW8p1c7ylp6yHJ2myqsNI2Dvrbni07n5qxPACgyGWSJpy7fl8X9aCC7+L2kHZOLD+P9nIfLPXZ+SWu6Zawq9///ioZxcdbk0uX6ZCxht6lCTTnA0pJmADSSnRbeeWtmFnemc8YqMjBsMPU4NGMdJUbKHbAuLGnF7yVtOS/re34tOZR+GYtDrnN1SUP+MG1Zb+rTUdbSQHbSJWM1ADOKO/NU0dl8nPMALxadxkWZk9lpqvNy8amckfETS0wz9lEFgE6ymh4Zy6kre9lqavJF8QA2mzrUln1sNrVZbFpQiSL2mcrkk802U5OWGbnMK2lHNfLJoogacoAmso1lJc0oIJsShKqST1O2stI0Zgc1qEQRxWTQQnJZb+qTRQk12M82alJEJsVkehNAQyFZlCAIhhocYC9VKCHDr+uBf0KW3MlZoqxy8bktnQIWke/wJH+XA4VAH2PMbyJyBPAucLUx5htHI41Rsp8Cbj3qS0fXr5RS0SjbTxBuzfqQ0zJ+4uiCZ6lEEQtzLuW14hGcmjmdRSWtuLzQMzikm6ygnuyiV8YyXio6nUKyAENz2cp6Uw9DBtXZzz4ql1l/dfazl6oMzZjDfys9zYn5j7HItGRl5QsBOCv/Pj7MeQiAQ/Le5ILMKdSSfdyS9TEAf5W04KyC+wH4X6XR9M6Ivj1icUlzfizpRlXyeK34JJrLFtaaBqwyBweSVCaf+7L+xxvFJ7DUNKcBOylB2EatkOutxV4qUUjfjMUsNi1YbppFHVskvtOo/u+pSh2rRo9w/DniPQXcE3gcSkurVwYwxswQkQeB0UBSJoDGmAIRmQMMAz7wu2sY8JE7USmlVHIKTCSeLTqLZzkTEArIZrVpxLVZXwDwaPGFpcvNN23BwNSSXn6PFtaZBqXXgo2u9d32bclhtM4bV3r70Pwn2W5qsIOaXF1wKz+VdKWAbN4oPoHjM2aVLrfD1Chdxz5TrneSJRcV/LNMHcOVpvwI4jxyuKvoytLrW6gdcb27qA7AVyWHxxSXFYHlcpSyyuohgwEKjKe5MBdo5XffWqCD3YHZ7BlgpIhcISKdROR5oCnwb5fjUkqpFHDwVJ3vNDDA9JIujj3jctOMHdQEYGJJ39JTnwCTSg4r/X+HN8kCyk3/trqkoaXncquItVJuspoALgTaef+fCdwqIh1EpBVwBxC6d3ESMMa8B9wC3APMA44ETjTGrHYxLKWUSjmfFg8EPCNEd1LDlRj8Wymb+M0CUlk8U8LNL2kNwDZvAqmUKs/qKeAxeAZ8ANwFTAL+8l7fB5xlb1j2M8a8ArzidhxKKZXK/lV0Pv8uOoWtLidXE4r7c1LmL+R7S4sANPQOALmn8DI6ZKynAbti6hOoVDqw1AJojHnLGPOw9/9FQCdgOHA60N4YM8m5EJVSSiWLYjLZSi3cHsH5fvFgTzzm4G7sp+KuACw1zfmw+Gi+KDmc6cXOnaZWKpVFTABFpLKITBKRwb7bjDF7jTGTjTGfG2NynQxQKaWUCrTQe5p3TPFJpbfdW3QZR+S9wH5vqdp1piEXFN7tRnhKJb2Ip4CNMXki0hfCVGJUSimlEmgrtcqMGgYoJIsN1HcpIqVSi9VBIJ8DpzkYh1JKKeWIrUYHgygVyOogkInAkyLSBPgKzzRqZSpIG2O+sjk2pZRSKm7nFNxLT1nO05W08pdSPlYTQN8Ej2d4L4EMeopYKaVUElpumrHcNKN70XIu8ZtKDuCuwstdikopd1lNAN2dsVgppZSK0/1Fl5ZLAH8q6epSNEq5y2oCaICNxpjCwDtEJAvPrBpKKaVUSjGRF1GqQrI6CGQl0CvEfT289yullFJKqRRgNQEMV/GzMpBvQyxKKaWUUioBQp4CFpHuQE+/m04UkY4Bi1UGzgaW2B+aUkopZa8LC/7J25X+VXrdWG4HUapiCdcH8HTgfu//BrgvxHIrgavtDEoppZRywk8l3cpcLzHuTmmnlFvCHfo8BtQAauI5BXyM97r/JccY084YM8XpQJVSSim7lbg8p7FSbgnZAugd8esb9att5EoppSocowmgSlMhEzsRqSoi5cq7iEhTEXlaRL4Ukf8Tkf7OhqiUUko5Q1sAVboK17L3NJ4p4EqJSCPgN+AmoBFwIvCDiBzmWIRKKaWUQ/LIcTsEpVwRLgE8Cngz4LY7gAbAicaYPkBr4BfgXkeiU0oppWw2JP/p0v/3UNXFSJRyT7gEsAXwR8BtpwFzjDGTAYwxecCLQG9HolNKKaVsttI0cTsEpVwXLgEswa8AtIg0wTMn8NSA5TbhaRVUSimllFIpIFwCuAA4xe/6GXjqAX4dsFwLINfmuJRSSilHzSo5xO0QlHJNuELQjwOfiUhLPK185wO/U74F8GQ8A0OUUkqplNAjbwx5VHI7DKVcE7IF0BjzBXABUAcYAHwEnGKMMb5lRKQB0BEY73CcSimllG12UZ18TQBVGgvXAogx5l3g3TD3b0EHgCillFJKpRSd4UMppZRSKs1oAqiUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaSZkGRgRKcEz84clxphMWyJSSimllFKOClcH8CYOJoDZwN+BvcBneKZ+awScClQDnnYwRqWUUkopZaOQCaAx5iXf/yLyDPAL8LeAmUBGAR8AbZwMUimllFJK2cdqH8CLgdf8kz8A7/XXgAvtDkwppZRSSjnDagKYCXQKcV+XKNajlFJKKaVcFnYuYD/vAI+JSBbwOZ4+gA3x9AF8CPivM+EppZRSSim7WU0AbwMK8SR7j/vdng/8B7jD5riUUkoppZRDLCWAxpgC4FYReRjojmcE8CZgvjFmu4PxKaWUUkopm1ltAQTAm+xNdSYUpZRSSimVCJYHb4hIdxF5T0SWi0i+iPT23v6oiJzgXIhKKaWUUspOlhJAb4I3B2gM/A9PYWiffOBG+0NTSimllFJOsNoC+C9grDHmaODRgPvmAT1tjCmtnN2nudshKKWUUirNWE0AOwLvef8PnB94N1DXtojSTMu6Vd0OQSmllFJpxmoCmAu0DXFfF2CNPeEopZRSSimnWU0AxwMPiciRfrcZETkEuBNPoWillKrwalfNjryQUkolOatlYO4FOgM/4Kn/B/AZnkEhk4DH7A9NKaWUUko5wWoh6HzgJBEZCgwF6gPbgW+NMZMdjE8ppWJ22cA2/N/0lW6HoZRSSSfaQtDfAt86FEvURGQqcHTAze8ZY871W6YO8AJwivemz4EbjTE7ExFjJCZwSI1SyjZtGlRzOwSlVAXRql5VVm/b73YYtrFcCBpARHJEpK2IdA68OBWgBW8ATfwuVwfcPw7oDQz3XnoDbyUyQKWUOwa0red2CEqpCuKMXmXLttWsHFUbWtKxWgi6qYhMAPYDS4H5fpc/vX/dst8Ys8nvsst3h4h0wpP0XWWMmWmMmYknQTxJRA51K+B41a+e43YISqUEEfvX+fiZ3e1fqUp5/dtoNbR007FJTbdDiIvV9PV1PC1ntwELgQLHIoreuSJyLrAZ+Bp40Bizx3vfAGAvMMNv+enAPuAIYHEiA7VLszpV2Lo33+0wlLKkTtVsduwvDHpfm/rVWLl1n2PP7UD+R9v6elo5nYkE77qjB+YVnxMHlG6yegp4IHCTMeZ5Y8xkY8wPgRcngwxjHHABMAR4GDgT+Mjv/sbAFmMO/ly9/+d67ytHRK4SkdkiMnvLli2OBR6Pk7s3cTsEpSyrlhP6OHNox4YJjESp+N19Yqegt5tycyQoldyiKQR9wMlAfETkERExES6DAYwxY4wxE40x840x44FzgGEi0jvW5/eus48xpk+DBg3seVE2u/zINm6HoBKgbQUZwBDuqNnpXaZUtEN2lbR0QF/66dmittshxMVqAngfcKeIJOKE93NApwiXX0M8djZQDHTwXt8ENBC/vYD3/4YcrGeYcnSnlrzaWUjarj461KQ6ZQ0+pOK3jnV2uA9Ndqb+VpS9dPurAL686UjuOD5lhxIA1vsAngG0BFaLyCxgZ8D9xhhzjh0BGWO2AltjfHg3IBPY6L0+E6iOpy+grx/gAKAaZfsFKmULKzuHRjUqJyCS5BGuZaRXy9qOPneDGtovS6lEe+n8Xtwwbq7bYTiqS9NabocQN6stgPWB5cA8IBtoEHBJeFOFiLQTkftEpI+ItBaRE/FMWTcXz0APjDGLgG+A/4jIABEZAPwHmGCMSckBICr1+edDWRnampDsmtWuUvp/q3pVtaeXCqp5nSqRF0pCfVvXsX2dJ3VvyrGdKv4ZjFRnKQE0xgyJdHE60CAK8MxKMhHPaN4X8ExLd6wxpthvufOB373LTfT+f1FiQ1VOOqVHU7dDiFm4ZELPNCWHZn479kEdyvYLPqN3M/56eHiiQwrriHb1eOpvPdwOo8IyIZq0K2VFVVY3aRzZPjn7uivnRf2NFY+mIuJqBURjzFpjzNHGmHrGmBxjTHtjzM3GmO0By+0wxlxojKnpvVyYLLOAgPOd4Cui6wa3K3PdyUTp6Sh3pFZCqVstO7ZgVNIRhMrZmW6HUcYNx7TnzN7NXHv+bs1S/9RYOJkhWu1TdRCIjl5OX5YTQBE5UUR+AfKAtUB37+2viciFDsWnVDmJ7HvhxKbx2E6NHFirSoQ+Dpwuq2huHtoh8kJxcLtlPEOEn+5046RXZKtGj3A7hLRwzdHtIi+UAqzOBHIxnjl0/wKuomxDxxLgcvtDU3Y5vZd7rQHpJtTO6de7h/otk17ndlP+5fodBZzas1nAXZ47q1VKrlbAimzlv9xPcprXqep2CLZp4X0tPZpHd2Cd8r/rOFSU6gJWWwDvBp40xlwCvB1w3wLAzbmA08rDp3aJ+jHPntPT/kBc5t8xP5mEGnXa0G/kr38folD9iWJVPUzRZbeEe4mO1wF0YC6QirHpV7Gy+zfrtkMb12Da7UO4dKDWl43kcO/c4hVljnGrCWArYHKI+/KA1J4QL4VUT5LJp0/sFnQiFUdMuW1Qudu+uumo0v/t2CFbKRfynIVE+p8ndOKx07vZEFFsquVoS5RVfzuseeSFgHYNI9d2TLdW3XCceiuOPqQBP96RnKdendShYXXHn6NlvapR9wW0+jFH27KYzPq1qcviR4ZzRPv6bodiC6sJ4FqgV4j7+gDL7Akn/VyWorN6PHlW4kYZtm9Yo9xttap6BlIcadMPsWmt4LX5hneNLtGtnJ3B+f1b2hFShZHqudHlR5Yt3F2x2n9Sx5uX9aNFXedPvU67PXyS2S4BCVk83rysX0yPc6ph84ZjnO0TGsiJsjb+crLKHmRflsItp1YTwP8C93sHe/jOvYmIDAXuAF5zIrh0kIyn7KwINb/rGXGOPrzpmPaWi/f+cPtgXru4j6OtL258Pr5BhimeN1kS6jW+fnEfR55vzj3HHnxui29wKiSw3StQK4vb6lWvFPK+KbcN4qgOwcumNEmSbilHHxJbWZeSCAlgm/rVqJkkZ6DCyUjwDzaVi81bTQAfB94C3gR8ZVZm4Kmr954x5gUHYlMuGN4lvlO7/zwh+ETpoRzVoWwLXnZmBoe1tHYE16peNapUykyyRCm6aLo1r11+DamQcdigcnbozc+xne0ZKe3/Vr55WT/qVT+4sXaif6BKfaEObiH42QjwtDpd0M9ay/9L54c6mRbcxUe0jmr5WJVEaAIUYM69wyx1hUlmE2480tZWwlQuo2O1ELQxxlwPHALcANwD3Ax09t6uEiTcTuuIdvF1TG1dryovnBfdxilegc3pEHlD5AgLSdeAON/fYP4X4+maaCT6M01WTk664nYqeduwQw5ecXl/lCbHL2UMPrQhGRa/YCd1j65w/UWHt4olpOhZ+N5kZ2YwonuT6FZbwQbNBIrl5QU2fLjFahmYQSJS3Riz3BgzxhjzmDHm38aYJSJSTUTK99JXCRdvU3SdapUsVbN3esqjSKci3FK3WuhTQ9Hwf3m1qpQvCh3PBvOl83uX/t+6XvT9pYLF4xS39wtWExUrcbZ1uV/Y4EMbOj4y0c73K15uTKF414kdbVtXrKdSnaw/F+7AOytDeO2Sst0yIp+pSMMjAYveury/2yEA1k8Bf0/oUi8dvfcrl7StH3mUohWB01wlwgVBBky40QLo21T1blnb0vKhBo2EO60ZrcDt67JHT4j4mL6t65b+/8l1A5lw45GWnqtjY8+pLbv6kt0zohP1vadbrxvc3pZ1RuvR07uSnRn687CzpWrsyL5ceVTydAZ3ohuB72cZaZBERXXR4a1tW9e0GEczO5n4htvq3n9yZ9o1SPxBTtsGB/dt5/RpYekxSdp+kJSs7q3CfeuqA/ttiEVZEKy/gW9HWyWOKak+v2FguQr+oVr6rhrUNujtsRjSseyE4Qa4dGDr6FZi4zaxl8X+h52aBK985GSB2KwwyUwwdapVomuzWpZaFI/r0phVo0fYVl/xiqPa0rNFbQDqBbScPnDywWPJYJG9ZsMAkDN6N+OC/pFOndn3xalTrRJHH9Iw8oJhPHFWd/51hnslhKxq5vAZACt8+a0LDYG2qF01trMJTp5ed/LA+87hHRl3RXStXn1b1ylTUqtRiIPuQG6dck7k2RO7hNyjeE/73ici93lvusJ33e/yGPACMD8h0argvBuFrs1qcUavZuV2uFY0rlm5XB+Wn+48ptxyq0aP4OIBrUOux2C4fkg77j4xusEg/o7q0CCm05f+nNhQpui+xhKnSyf4y/EeqITaTA+LMADkPAud7QcfWj4ZC2zBsPIdefT0ruWWs7J/iaU0RP82dS29tmSXiD6Avs/gxfN6h1/QRuFe14hunn5xdao6mwTYefAdKNq8KZqP+drB7WKqnRfLvswKKwPARkY5+CYVy3+Fa1LoD9zovRjgb37XfZdLgK14BoYom3x+w8CYH/vMOT2ZcJO1035Ouf34jlwZ54Yq3mO4O4dH7q/zyXVHxPksqSHS6cB3rugfsrSFE+LND+4cfmhUy9f27pQDW8itxFGjcnbYHaNvxw/QqYnnNPrDp3XlvpOTZ3KkBQ8eb2k5pwv2Ln30BKaPKn9QGY+hnRpy5VFtElIsOZzW3m443/19sKPPU6NyNt/cclTkBWNgteUs1GJOlG7q0KgGg6Isa2NXA+ADp3QJ2dWn7PN5njAVGwhCJoDGmCeNMQ2MMQ2ANcAQ33W/SzNjzFBjzG+JC7ni6x6kNIgVvv18k1pVuGdE7C1wyc7KaNxKmRkc0ij8TsH/dK9/jmRXcelUkYhah3b2jbQiERvjQxpVZ2ing62V9arnsGr0iMSN2rQoXFkTu6waHXl+3uzMDNuncMzMEO4e0Tls7T67WGnZrOPXYhVpO/LBNQNiiqNjY2cm3op28F3g+2FX6aZAVvtl+8SS//Xz6zsdq1Qc/W61DEwbY8w8h2NRfk7rGbxUQLCm60pB+oZdcVTZFrhI/QMzk6QzTXsLR/KBo3FDNefHUufNGBh7ad/g6/P7hcdzkFmzcvDTREe0q0eGwJm9rU1RZhf/DVeojVi4o/ubhkau9O/fP6b0OaJ4E/2T1Mpx9HX1Z2WDHa5VJNaWhosHBE8QG9W01scplETVI4s89tP+bYmTB2WNI7zvsQ686Bshqejbui4tEzCziVVhv+uxrTHWUOISSx/Ad67sz6fXlz/zVtFrsobrA3iiiNT0+z/sJXEhp4fnzu3FOxY6zT50ahdL0yNFmtLMv0Cuv6n/GBxx3XY6sVt0NaYCHdrIcxqud6vo+rT5/8yjHWxhl5Z1q7LiXyPoEKHlMlp2dIoOd3Qfzc65X5u6DO8a/Wd8eNt6rBo9glWjR5BjoVSRFbEmKvHsExrXrMzAIO/XJQNa2ZbY2p1/3XrsIZEXcshhreow995h/HekMzPDQORBR75fjxuFwwe2d7a0jz+rLYBRzxkc5Q9m8KGeU76RWq5/u3dY0INPS68jIKTszAx6tqjN8sfKpjJWtp3+i6TambdwW9IJeEq8+P7/wvs32OULB2NMGw1jqON3bl9rHU9j3XS1jrbEjItj8J88qzsTbx3E4keG07NF7ag2VL6NVCpXdU+UwBbYJhb6ybx0fm8GH9qAd688PO7TwSJiqSRO5PVYXdK+78Qn1x/hev3DaN18bAdWjR5hvQ6gje/Xh9cMoE61SuUKxrvxFkab+Ft5HyI1LtrRN/ftIDXngtU0jfY9DZYQX9C/JaPjGMl+Xr8WnO8dCJWV4dlOhPq91K1WqWwBdK9Ir6N2mIE6gWfCojktLki5M2/JLtyWuA0wz+//tt6/wS6p9aqT1EOndrXUlyYaTnaOtnLqLxFuGtqeni1qc5x3Grtgs4skm3tGdErZQSi/3TuszHUrLdB9W9dl7KX9yMwQsjM8fcFGnxl5RxGq/48drbS+TX21Ss58X4K1VDapVbYP3KRbB9GgRg4XhRlZHzWHsqNkTFwT0SbnZFmRZxIwrdqRFmed6N/GWj+4cG/Ho6d349w4RrIfHqaYeajPOtqZSaIp1WUliQ9cwr90TbILNwhktTGmwO//sJfEhZx+QtWci+SH2wfzsY1JxhsBfeOCHX2FEthPsV0De4pXg2dO4E+vH1iuDlOsfQBDsXNnc8VRbS3XHEws53epGRnC9FHHcEbv5mXe7+o5WZzS42Df11WjR4Sd+SDcpPfhjvJ9fP1No5083mo+EKmExXGdG3FIoxrMuvtYS31fI7HjFKWVucAjvV2pMMfyfy46rPT/SDv5g6eArfF1ifD/npwUIklxelalUIIltV2bRTkKPEk+5pfPDygFZOEHWtdiHUZfC2DVMAeJvqfz/S7O798yZEPON7ccxY8xFgF3QlSH0SKSIyJtRaRz4MWpAFXstaVa1atGDe+Ag2A/iRO7Rd7Y+xsc5XB8f12blU1i376iP7/ePTTk8snY2uDPyVYBu1ddtVL4vjTJssP+7d5hUc1b/GaQeZRP79WMf1/Y29Kps1N6NgOgYc1IXS9CDDJKjrfNVvEUs33zsn62lAI5o3ez0v+d6oR/vIVEN5DVWNp4u834H4Sc2rNZqMVTQp0YC1eDtW1luDNV0WwOOzauUa4EU7C+9I+f2Z2HT+0ScX2+0L+66Shm33Ns2GWtfDs6Nq5p6YxJolidC7ipiEzAM+PHUjyFn32XP9FC0AkTaTBHNC4/MnFn7gM3nlWzs2hYI75RjxVFqD500e783hgZfPTysZ0a0qKu+7M3JEKlzIxyg0yO8RaFzgnoe1irSjbPntOD/0Wcl9OZMZBJfowTtaMPaWCpFEikOYufObtn1M8dba24aER7QHb3iE6MvbSvpWkV7Tr4stIP1ydUa6Tdj4mmifBQ71SUYddmYXVvXtaPNvUjt6bXqpptsduF58OvlpNVOuNWRWK1BfB1oA9wGzAcOMbvMsT7V8XIN/OBlS948JGCqbUreepvPajlcMX8aFnZVJUtl2Jfy0Q0fSkPa1WHGiEmkg91GlFEXKlNVy3H811N1Hy8odY3+szuTB91TNCW0NN7NY9p8FU0vv370WWuJ7Ll0L81LV59ohxZH8qZhwUvc9SsdhUm3jLI0joCW5WuDdNVwC5WP7bK2ZkMPrRhQs5i+LoHvXrhYfx+33GWHhNLya9g2zsnvsZln8f9/VpJwOndYGId9PTGyL5Muc3a990pViuEDgSuNMa872QwKj7RH00m7gfmv8E+K8QOIF1Fs6P46FpPn87Wo7505DnCbejqVqvE3rwiy8/50Cldad+geti+enYLFn+lLPsLEEcjnjm64zH/geNsfe7/u7Qva7btj3zwE2NmEMssIb5YAmPq2LgGf23aE/X6wg0ISuZT/lkZYvmgul+bukxbsiXsMoc2qkHzOlX49q9cS+u0M5Gx+7R/PGuLapaPKOMe0jG+ucPtYLUFMBc44GQgSrkpGTbuTsYQLv9r7Hf6KFyi+MtdQ1nwkLVpxcBzmuXGoR3KzTGd7C4e0CrqkYXBuN2PtUbl7KhHS4f7DtasnB39QIE05n77lefgx9/Mfx5TWmYlnMqVMvlviC4lOVkZnNC1cZkuJ+0bRj6Fa1Vg626b+p4+c63rWRs4GPj4eM42+dYUbqCY27/zeFjdOtwH3OkrDK2c4eQXKZ6Op74p1VL5ix5Kl6Zlv9JODu5wU7iXVd/iNFrZmRlhT+fGy476cX/r0yLudVTOziyzwY/3K+H0wYXV06aJEunl9m2djKPfQ2vkHSSUirNCLHmkbL3MJrWqWHsdYb70IsKrFx7GETHMzvJiFIO8fN+k03o246NrB3Bqz7IVAqwYd2V/ujStxRi/Ud/RKPGeAw5/Ctg/2tQS8hSwiASe7m0JrBaRWcDOgPuMMeYcm2NLG6G+OFZ2PCLWlrvpmPb0bFGLy8bOjio2gHevPJxFG/fE3ZIzonsTx2quResob22sVvWqsmDDbkunz8tMBWcxK/j1rqHkFZbEFqSN0qHA9ZVHtaG3hdI6x3dpxKqt+2N+HqvvZbzvec8WtZm3dmfE5ax0oE+kSDPDeKoTZLEniu4E0bjo8Fa0b1id+z9fYGn5SD/lj68byG+rd5S7/eubj+KE53+Meb2Q/L/Lj64dwI9Lt9q2vpN7NGXMtBXMX7+rzO3h3gUR4bBW1ufq9V+Xb8DRcTGM+vZfl7X9Q0xP4apwh/MNgQZ+l+V4CkNnB9zewLusclm4L2BWZgbHdLQ2WXdgyZZ61XMsFxMNRUR4+fzePHFWj5DL3DCkfcj73rvq8Liev2ws8PolkctVnBpiPuZoNKxZmZb1Ym99jWbWjLBHqQ7tZ8KV8rHCzristtD856I+TLw1dKuZ3XMxJ0uZnUAzgvS5Oy+OIr4+iZhKcbB3ZHewd/bh07oG/d1EU7fUX7PaVTi5R/ltQTT1WZ06sxDrN8tqNIe1qsstNk8F+P7VA0LeF+8vxf/xdatVir/V1jcIJMxXuqP3AOyQRmUPxM7v35JjkqCfXzghWwCNMYMTGIcKwq0jii9uODLmx8azmfv7cQc3NIFHxlVsbDmsW7X81FJ/69Oc5nWqcLN3Y/fXw8MdPd3pL1wrwJx7htHl/olRrzPSad0R3Zrw5fyN5W6P9juXTKV84v25ZGUIyx4LPa25rx6n1YSuU5OabN69JeTcxU4l5eFWO/ueY+nzyBQAmgYZGNOjRW2WPnoCP6/YxkX//dWZAIn/s3rlgt7k7s6PeFZiYPt6TF+2rfR612Y16ROkNemdK/pzweu/xBlVWYc2Sq6WWSustjrHw/q2vPw3+Yh29cImVXb/pEosDAI5sVsTJt4yqFxLvG9GkGgH7CVSyD2ciBwuIslVq0PZKtQOSESiPnKyMvNCPLo53Pm8aqUsnjmnZ+kcmZWzM2MqlxAf33zEB1XLyeLmoR2Cto6GKu3y18PDy42oDGyBOMymkh6JE3nTXifCzBuh+D7lSC069aKsA/bieb0Yf9Xh1KlWiaqVMktnBon3WzX20r78cPvgcrdb+claqWWWnZlhyxy0Vg05tAE1cqwWpPConJ1pqWW9We0qZX47E248igdOKVsA2OBJfO1Wq2o2x3aKrwUo0kGC3QlP4HvjhnD7nnFXHm55vl07tt5WBoFA8nXDsCrcr24GkCcic4Dp3ssMY8y2MI9RaSpwqrdA8Z4CcaoDtm/mg2hOtSaK7xXfGuTUla8T9Fs/l5+FMVityNoBlfyTu+fRQYMPbcg5fVoEfQ/8/euMbjGXF8rKzOC9qw63fSNeo3J26dymf9xvrUZbOB9dO4CGNSo7N5OAS6cc3ri0/Iwu8apXzZPoNq1dhaLiKEp5RMlqX8ZQ269ougeEOyBNxf5nTom3Zf2NkX1ZsMHTR7G0BbCCvr/hEsDjgQHAEcDVwB2AEZGleJNBYLox5i/Ho1Rh3XBMe1Zs2ceJ3aIvXfHgKV0cL4SbzO4Z0ZkODWtE3VfDyQQq1g1YuCT5vH4tqZSZweRFm5m8cHOMkSVepawMHj+re8Tl4u271j/C7BTxiqVfXODHGU1HeKfdMKQ9v6zcxqxV5QdH2OGz6wcyaeGmmB8/tFND/nPRYQzt2JDnpiyNeT3h+qsBfHPLIJZujlxvMJoD4FpVstl1oLDMbc+e04PuzWtbXoeKPWkb0rFhaY2+0nl+k7Qfb7xCbpWMMZONMQ8ZY4YDdYHuwHXAz8BRwGvAAhHZKiKfJyTaCurGYzrQsEYOh7e1voHP8jsabF6nKu9fM8DSPJ4XHl52R3nJEa05IYbEMVrRtuANbOcZdHKCjVPfBVMtJ4vLjmwTNr7nz+1ZrlyMz2k9mzLltqOD3hc3G7c5mRnC2X1blH5vQu2QKuZmLj6p0loai8m3DuJfZ3SL+nH/OP5QPrjmCAci8ujRoja3H98x5seLCMd3aWw58a6anclRHerz6oVly4X0axN+m9ysdpXSASkhIilz7bx+kcsUzRh1DOf2Lbvc6b2a065B5CnO7GLHd963r7GUuNr4I7NzZHXpKOAKumG01PHCePYWf3ov/wEQkcHAP4ATAGtFeVRQ3ZrX4te7w080bRf/I5loRrElmq8/VzIUnj21Z7OQk7mf2rNZyCnYkpFTG7IMOThtUkXle+uSvXRHNDo0qsHefO8pzIpaAzPC52WMISNDeCtgTuiWNp9qX/mv0AOMfE7p0ZRqOVmWDuaT3eBDG1qu1+dj7/Yp/pV1bVqT39bsjNgHMFVZSgBFpBrQH8/p4COAw4EawEI8LYEznQpQOeOhU7tQLcqO11bEsw9JVKHVZPstn96rGTOXb4u5TIVPNC/L7s/+m1sG8cvK7bzz8+qYpuBKRfGfFgr9Y4l2UEQ83Chw7M5zWl928q2DaGBz1xgrrzmwH2qkh8Tat7pbs1ph6/sl2SbSMWf2bs5Hv60Lef8bl/ZjWe6ecjOqVBThCkGfz8GErxuwC8/p3xnA08Cvxpj02NInsWQ5aA+2oZo+6hiW5+7l4v+zp5xEjcpZdG/ufoug3apWyuKl83vHvZ5qlawnDTUrH2xhsGNnfEijGhzSqAan92rGjn0FUTwySb7AAS4Z0Iovft/A4W3rleuPFa9I7/cPtw+mRuVsrnhzlq3Pmw56t6xtedkeLWrze4iSJx1SrIRLtAcjr154GF3DlJdKhl9lzPu2EI975uwebNmTX+a2J8/qzuNnhu4CUatKdlL1vbVbuD3G28A+4H/ARcYYa2XVVUIkWytWMM1qV2HjTvumkJ7/gPV5aMOpbkPrSrIk3j7Tbh8S15yXPvF+r6rnZMX8/ibTe9qndd3S01e+0z9HH5KY0iitLM55qsqadfex1Khs/bt3fr8W/L52p+2nelNB9ZwsujStyYINu8vcngy7FbtiCNyWnRGkwHtGhpCRFK/aHeF+LU/iGQV8KTDSWw5mpvcywxiTm4D4VIprU9+zM7NjhgE73HtSZ47rbG1GFEuSZNsRz2wj7jr4BgYrXxPJm5f144PZa+0MqJwGNXKY+c9jaBBlHcB08PiZ3Vi9LfZp9ewU7Snbc/q25Jy+ybFdCiWZDopShe9AON4uNekg3CjgO40xg4CawGDgQ6AV8AKwSUSWi8jbInK9iMR//kpVSPWq57Bq9IiYa7TFo32j8oMzLj+yjXN11FLABf09xaN7t6rtbiAB2jaIrdXr6EMa2HL6PJImtaokZIqzeFxyRGsg+qK08dToPKdvS+4YXn607hPe0j01o2iRc8rFA1rTo3ktzu4beQRu0rB4YOkbpex0IX6n+aYajXekswA5WZmsGj0iaRodklnELZoxpsgYM8sY84Ix5lxjTEugBTAKqA08Dzg3Z5AKKtZttq/ivd0jV6t6+58l06npJ8/qztsBI/vS3cD29Vk1egRNapWfBkyltuO7NGbV6BGWZvsIysYf79l9WrBq9IiQJUAiTVVop0Y1K/PZDUe6Om2h1c31gHbR1aO896TOTL51UNBp/VLJuX1bMOeeY5O6MkVFZPnwTERygL4cHBgyAPB1illjf2jKimi32Wf2bka/1nVtP2X45qX9mDB/A41qRr+RPb5LIyYusL84cdVKWaVHlkrZpbH3O55SLUouqRci0Xvs9G6cM+ZnOqboFFpWRdo+++6vV60Sc+4dFvX6szMz4h6w4sZp5qM61OfHpVv532X9+O6vXEQk6qkWVfzCjQJuysFk7wigJ1AJKALmAePwjAieYYxZ73Sg6cg3UrNFnapA8Bn4ov3xiogj/cVa1qvKdYPLz1drxQvn9WLn/rIjLX2tGHVjnN810UTs25BWpDpzFVHtqpWirm8WTDr073rktK5s31dQbv7Wni1r06dVHf55YieXIktt467ozwdzQpcvscLNszVvXtqPEmPIysxgUMDgquZ1Pa2Zh7VOtfnKU0+4FsB1eFqud+AZ+PEgningZhlj7BvaqULq1rwWYy46jEGHNOC9gI7uyXSqNV45WZk0qll2AMDII1pTt1o2p/YIXoDZbef1a8kPS7bQ2XvKYtrtQ1izPXxn+GsHt+PLPzZafo5knn7olQt68/WfsU/VpSJzo1ae3WpUzi5XYBk8v/kPr3VuJpFk4VSSf0T7+hzRPnXPboQbfduxcU2+/8dgWlWwvtrf/d2hGaPiEC4BvAJP657O9eui47o4OxVassrMEE7vlfiBI1YN79q4TCtQi7pVIw4uuXN4R+4M0mE+HpNvHRTT6Nl4nditSUxzT5eXBs1gIVSA/E5VUG63TvuqR1QkbRM4lZ9VIRNAY8z/JTIQpVT04un/07BGDrkBhVHdormQqohiTfJzvCPOszMT+8vQg5L04v4YfaUSZG4Mnawrsil/P5p9vnlgvZL5tLNymNvNPqrUtYPbU1hiuGhAq4Q+r34F0osmgClOBwxYV8emASWJLGHhpJqVs8tMCadSy70ndaZf6+inqZox6hjyi0pKr4fra9i/TcWdBiuRQiVWoW6vUinT9u4iSgVK6gRQRK4CzgN6AbWANsaYVQHL1MFTnPoU702fAzcaY3b6LdMNeAnoB2wH/gM8bOKpgOoyT0tNyoafsibeMqjCJIAqtV1+ZJuYHme1Ztzyx07U9uA4WX3/kuXUa7LEoRIjqRNAoCowCfgMeDbEMuOAlsBw7/XXgbeAkwFEpCYwGZiGp45hR+ANPPMcP+1U4Cp5/Oeiw6hayZ6BEtHOsqBUKMly+BnqODgzQ7OBeFn9iJPlu5DqcrISPyAulSV1AmiMeQ5ARPoEu19EOuFJ/I40xsz03nY18KOIHGqMWQxcgCeRvMRbvuZPEekI3CYiz6RyK6Cy5vg0HUkdC20BcJ7lViHv3/tP7kxP7ww+jtIP3zGh3lp9y+1VK8WnxEu05J7cMrIBwF48Bal9puNp3TvCb5kfA2oXTgSaAq0TEKNSSsWsa7Na9GqpRXFTWTI2M/gOKoLNI6x9y9NDUrcAWtAY2OLfimeMMSKS673Pt0xgyfTNfvet9L/D2+/wKoCWLXUyaRWb03o2ZdaqHW6HEbOmtSqzYVee22EoldKSuYHvvpO6cG7flhHrl6aCl87vxSbdXkUt4QmgiDwC3B1hsSHGmKkJCKccY8wYYAxAnz59kvYwqGaVbLbuTY4abqq8587tFfNjff1YKme710B/1aC2jBzYhtajvnT0edrWr86lA1tz0eGJLXeRSpKx9UilvkpZGXRtVivofalWDuqk7k3dDiEludEC+BzwdoRl1lhc1yaggYiIrxVQPDUNGnrv8y3TKOBxjfzuS0kfXTuAaUu3aqfXCujCw1uya38B1xzdLuHPnejNfkaGcP/JXRL8rKlB+4cpt+gp4PSQ8ATQGLMV2GrT6mYC1fH08/P1AxwAVPO7PhN4XEQqG2N8bcTDgA3AKpviSLhW9apxUb2KN12O8rQA3nbcoW6HoZSK0/FdGjNp4WYOaZR804AFk2otfyo+Sd0HUEQa4+mnd4j3ps4iUhtYY4zZboxZJCLfAP/x9t0DT42/Cd4RwOApE3M/MNZ7+vkQYBTwoI4Ajt/NQzuwv6Ao8oJKqTIibXzqVPXUm6yUlepj9dLXmYc156QeTfRMjUpKSZ0AAtfgSd58fB2SLgXGev8/H3gRz8he8BSCvsH3AGPMLhEZBrwMzAZ24Kn/94xjUaeRW4cdEnkhpVQpq6d2nzyrB0f9sYEezYP301KpwY3k74L+LRnaqWGZ29o3TI1WSJU4SZ0AGmMeAB6IsMwO4MIIy8wHBtkWmFIV1MiBbZiyKJcTuzUB4NPrB7J9nw42ckOtqtm2DY751xndWLl1X9D7qniLpDeskWPLc6no2X0q6tHTu5W5PueeY0s/53Ba1feMCL56UOL7H6vES+oEUCmVWG3qV2P6qGNKryekALFy3Hn9Qpe06ti4Jk//rQfHdg4cK6eclqged/WqW0vua1bOZtXoEQ5Ho5KFJoBKKZXmzjysudshKKUSTBNApZRSSikHnNm7OfWqV3I7jKA0AVRKKRdoEQKlKr6nz+7hdgghaX0BpZRKIC3wrJRKBpoAKqWUUkqlGU0AlVJKKaXSjCaASimllFJpRhNApZRSygVZmZ5dcKt6VV2ORKUjHQWslFJKuaBWlWz+e0kferWs43YoMXnlgt7szdO54FOVJoBKKZVAbep75mQd0rFhhCVVOhjaKXVnYPFNGalSkyaAKWLCjUfyy8rtboehlIpTm/rV+P2+46hZRTe/Sin36BYoRXRtVouuzWq5HYZSyga1qma7HYJSKs3pIBCllFJKqTSjCaBSSimlVJrRBFAppZRSKs1oAqiUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaUYTQKWUUkqpNKMJoFJKKaVUmtEEUCmllFIqzWgCqJRSSimVZjQBVEoppZRKM2KMcTuGpCUiW4DVbseRAPWBrW4H4SJ9/fr60/X1p/NrB339+vrT4/W3MsY0CLxRE0CFiMw2xvRxOw636OvX15+urz+dXzvo69fXn96vX08BK6WUUkqlGU0AlVJKKaXSjCaACmCM2wG4TF9/ekvn15/Orx309evrT2PaB1AppZRSKs1oC6BSSimlVJrRBFAppZRSKs1oAqhKichUETEBl/Fux5Vo4vG19/Wf5XY8iSIir4nIchE5ICJbROQzEenkdlyJICJ1ReRFEfnL+/rXisirIlLP7dgSRUSuEpHvRWSn97vf2u2YnCQi14nIShHJE5E5InKU2zElgogMEpHPRWS993Me6XZMiSQi/xSRWSKy27ud+0JEurodlxs0AVSB3gCa+F2udjccV/wdKHE7CBfMBkYCnYDjAQGmiEi2m0ElSFOgGXAH0A24EBgEvOtmUAlWFZgEPOByHI4TkXOA54HHgF7ADOBrEWnpamCJUR34E7gZOOByLG4YDLwCHAEcAxTh2c7VdTMoN+ggEFVKRKYCfxpjbnA7FreISF/gY+AwYDPwN2PMh+5G5Q4R6Q78DnQ0xix2O55EE5ETgQlAbWPMbrfjSRQR6QPMAtoYY1a5HI4jROQX4A9jzJV+ty0FPjTG/NO9yBJLRPYCNxhjxrodi1tEpDqwCzjNGPOF2/EkkrYAqkDnishWEVkgIk+JSA23A0oU72sdB1xljMl1Ox43iUg14FJgDbDK3WhcUxPIB/a7HYiyj4hUwnOANyngrkl4WoVUeqmBJxfa4XYgiaYJoPI3DrgAGAI8DJwJfORqRIn1b+AbY8zXbgfiFm+/qL3AXuAEYKgxJt/lsBJORGrj+Q28ZowpcjkcZa/6QCaeFn5/m4HGiQ9Huex5YB4w0+U4Ek4TwApORB4JMrAj8DIYwBgzxhgz0Rgz3xgzHjgHGCYivd18DfGw+vpF5CKgB3C72zHbKZrP3+sdPH2ijgaWAB+ISFUXQrdFDK/fd0roC2A9nj6BKSuW169UuhCRZ4AjgTONMcVux5No2gewghOR+niOeMNZY4wpd5pLRDKAAuACY8x7TsTnNKuvH0+n4IspO/gj03t9pjHmSGcidFacn38lPKdFrjHGvOVEfE6L9vV7k7+v8AyAOcEYs9fhEB0Vy+df0fsAer/X+4HzjDEf+N3+MtDVGHO0a8ElWDr3ARSRZ4FzgSHGmL/cjscNWW4HoJxljNkKbI3x4d3wJEEb7Ysosay+fhG5G3gq4Ob5wD+AzxwILSHi/PzFe8mxL6LEiub1e/uAfo3nNQ9P9eQP4v78KyRjTIGIzAGGAR/43TWM9OrykrZE5Hk8Z7jSNvkDTQCVl4i0w9P/7ys8O4zOwNPAXGC6i6ElhDFmPZ5TfqVEBGCtMWaFK0ElkIi0x9PncwqwBWgOjMIzCGKCi6ElhDf5m4Rn4MdpQDXvQBiA7caYArdiSxQRaYynD9wh3ps6e/tCrjHGbHctMGc8A7wlIr/i2b5dg6cU0L9djSoBvK3c7b1XM4CWItITz/d8jWuBJYi3pfciPL/zHd7vPcDeinDQFw09BawAEJEWwNtAVzx1otYCXwIPVsCNvyUiYkiTMjDez38MntGRtfF0iJ8GPJwOR8jefnDfh7h7iDFmasKCcYmIPADcH+SuSyviKUIRuQ5PH88meOri3WqMmeZuVM4L811/0xgzMqHBuMC7XQ/mQWPMA4mMxW2aACqllFJKpRkdBayUUkoplWY0AVRKKaWUSjOaACqllFJKpRlNAJVSSiml0owmgEoppZRSaUYTQKWUUkqpNKMJoFJKKaVUmtEEUCmV9kTEWLgMFpFVIhI4ZaAb8R4mIjtEpKbF5c8SkcUikul0bEqp1KCFoJVSaU9EDve7WgX4DngEz2w4PguBdsA2t6fMEpGvgXnGmH9aXD4D+At4rCLO6qGUip4mgEop5cc7V+oeknQKNBHpACwBDjHGLI3icfcApxtjDnMsOKVUytBTwEopZVHgKWARGSsis0VkhIgsFJH9IvKliNQVkfYi8r2I7PMu0z1gXRkiMkpElolIvogsEZFLLIRxCfCHf/InItki8pSIrPGua4OIfCIilfwe9xHQW0S6xPs+KKVSnyaASikVn5bAQ8A9wFXAEcAYYLz3chaQBYwXEfF73Ivex4wBRgCfAP8nIidFeL6hwIyA2/4JXADcCwwDbgF2AaV9/owxi4AdwLHRvkClVMWT5XYASimV4uoCA4wxywG8LX23A5cYY/7nvU3w9CfsCCwSkfbAtXhOM7/pXc8UEWkC3A9MCPZE3vX0At4OuKsfMM5vXQDvB1nFH95llVJpTlsAlVIqPqt8yZ/XMu/f74Lc1sz7dyhQAnwiIlm+C/At0DPMaN06QA6wNeD2ecBIEblDRLoHtDT62wo0jviKlFIVnrYAKqVUfHYGXC8Icrvvtsrev/XxnJ7dFWKdTYB1QW73PT4/4PZH8CSU1wGPA+tF5EljzPMBy+X7rUMplcY0AVRKqcTbDhQBA/EkboFywzwOoLb/jcaYPOA+4D7vKOFrgOdEZLEx5hu/RWv7rUMplcb0FLBSSiXed3haAGsZY2YHuRQEe5A30VsDtAm1Yu/o4H/gae3rHHB3azwlZJRSaU5bAJVSKsGMMYtF5N94RgY/AczGc2q2C576fleEefh0oEwtPxH5BJgDzAUOcHDk8TS/ZarhGYRyr40vRSmVojQBVEopd1yPpzXuSjxlZHbjmW3kvxEe9zHwhohUMcYc8N42AzgHz+jjDO96zjTGzPZ73HHAfmCiba9AKZWydCYQpZRKId7izuuA640xH0TxuHeBfRFaF5VSaUITQKWUSjEicjtwqjHmSIvLtwAWA92NMcsiLa+Uqvj0FLBSSqWel4CqIlLLGBOqlIy/5sA1mvwppXy0BVAppZRSKs1oGRillFJKqTSjCaBSSimlVJrRBFAppZRSKs1oAqiUUkoplWY0AVRKKaWUSjP/D2eMXjWISRJiAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<Figure size 648x432 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(grid, sample[det_string], color='C0',label = 'strain')\n", + "plt.plot(grid, sample['h1_output_signal'], color='C1',label = 'signal')\n", + "#plt.xlim(-1.5, 0.5)\n", + "#plt.ylim(-150, 150)\n", + "plt.xticks(fontsize=14)\n", + "plt.yticks(fontsize=14)\n", + "plt.ylabel('Whitened Strain and Signal ({})'\n", + " .format(det_name), fontsize=15)\n", + "plt.xlabel('Time (s)', fontsize=15)\n", + "plt.legend(fontsize=15)\n", + "\n", + "# Adjust the size and spacing of the subplots\n", + "plt.gcf().set_size_inches(9, 6, forward=True)\n", + "plt.tight_layout(rect=[0, 0, 1, 0.9])\n", + "plt.subplots_adjust(wspace=0, hspace=0)\n", + "\n", + "# plt.show()\n", + "plt.savefig(plot_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.9.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/python3_samples/config_files/default.json b/gen_samples_py3/config_files/default.json similarity index 70% rename from python3_samples/config_files/default.json rename to gen_samples_py3/config_files/default.json index 271ec00b3cb89b33fbe77ea088d338a611c5eca9..b830806fe5f8385aff498f36f092f09480cf703e 100755 --- a/python3_samples/config_files/default.json +++ b/gen_samples_py3/config_files/default.json @@ -5,8 +5,8 @@ "inj_bits": [0, 1, 2, 4], "waveform_params_file_name": "waveform_params.ini", "max_runtime": 60, - "n_injection_samples": 4, - "n_noise_samples": 16, + "n_injection_samples": 100, + "n_noise_samples": 100, "n_processes": 4, - "output_file_name": "default.hdf" + "output_file_name": "train.hdf" } diff --git a/python3_samples/config_files/waveform_params.ini b/gen_samples_py3/config_files/waveform_params.ini similarity index 75% rename from python3_samples/config_files/waveform_params.ini rename to gen_samples_py3/config_files/waveform_params.ini index c6a3477608477cce042a8af37110bc5fbcd8e477..3dadf30be35681ef5a3c998df9de306039d9a75f 100755 --- a/python3_samples/config_files/waveform_params.ini +++ b/gen_samples_py3/config_files/waveform_params.ini @@ -1,9 +1,4 @@ -; ----------------------------------------------------------------------------- -; DECLARE ARGUMENTS -; ----------------------------------------------------------------------------- - [variable_args] -; Waveform parameters that will vary in MCMC mass1 = mass2 = spin1z = @@ -15,9 +10,7 @@ inclination = polarization = injection_snr = - [static_args] -; Waveform parameters that will not change in MCMC approximant = SEOBNRv4 domain = time f_lower = 18 @@ -57,60 +50,39 @@ seconds_after_event = 2.5 tukey_alpha = 0.25 -; ----------------------------------------------------------------------------- -; DEFINE DISTRIBUTIONS FOR PARAMETERS -; ----------------------------------------------------------------------------- - [prior-mass1] -; Prior for mass1 name = uniform min-mass1 = 10. max-mass1 = 80. - [prior-mass2] -; Prior for mass2 name = uniform min-mass2 = 10. max-mass2 = 80. - [prior-spin1z] -; Prior for spin1z name = uniform min-spin1z = 0 max-spin1z = 0.998 - [prior-spin2z] -; Prior for spin2z name = uniform min-spin2z = 0 max-spin2z = 0.998 - [prior-injection_snr] -; Prior for the injection SNR name = uniform min-injection_snr = 5 max-injection_snr = 20 - [prior-coa_phase] -; Coalescence phase prior name = uniform_angle - [prior-inclination] -; Inclination prior name = sin_angle - [prior-ra+dec] -; Sky position prior name = uniform_sky - [prior-polarization] -; Polarization prior name = uniform_angle diff --git a/python3_samples/output/default.hdf b/gen_samples_py3/output/default.hdf similarity index 100% rename from python3_samples/output/default.hdf rename to gen_samples_py3/output/default.hdf diff --git a/gen_samples_py3/plot_sample.py b/gen_samples_py3/plot_sample.py new file mode 100755 index 0000000000000000000000000000000000000000..12fc8e883f73d1035c5874a01d11950d14cc2773 --- /dev/null +++ b/gen_samples_py3/plot_sample.py @@ -0,0 +1,197 @@ +""" +Plot the results produced by the generate_sample.py script. +""" + +import argparse +import numpy as np +import os +import sys +import time + +from utils.samplefiles import SampleFile + +# We need to load a different backend for matplotlib before import plt to +# avoid problems on environments where the $DISPLAY variable is not set. +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt # noqa + + +# ----------------------------------------------------------------------------- +# MAIN CODE +# ----------------------------------------------------------------------------- + +if __name__ == '__main__': + + # ------------------------------------------------------------------------- + # Preliminaries + # ------------------------------------------------------------------------- + + # Start the stopwatch + script_start_time = time.time() + print('') + print('PLOT A GENERATED SAMPLE (WITH / WITHOUT AN INJECTION)') + print('') + + # ------------------------------------------------------------------------- + # Parse the command line arguments + # ------------------------------------------------------------------------- + + # Set up the parser + parser = argparse.ArgumentParser(description='Plot a generated sample.') + + # Add arguments (and set default values where applicable) + parser.add_argument('--hdf-file-path', + help='Path to the HDF sample file (generated with ' + 'generate_sample.py) to be used. ' + 'Default: ./output/default.hdf.', + default='./output/default.hdf') + parser.add_argument('--sample-id', + help='ID of the sample to be viewed (an integer ' + 'between 0 and n_injection_samples + ' + 'n_noise_samples). Default: 0.', + default=0) + parser.add_argument('--seconds-before', + help='Seconds to plot before the event_time. ' + 'Default: 0.15.', + default=0.15) + parser.add_argument('--seconds-after', + help='Seconds to plot after the event_time. ' + 'Default: 0.05.', + default=0.05) + parser.add_argument('--plot-path', + help='Where to save the plot of the sample. ' + 'Default: ./sample.pdf.', + default='sample.pdf') + + # Parse the arguments that were passed when calling this script + print('Parsing command line arguments...', end=' ') + arguments = vars(parser.parse_args()) + print('Done!') + + # Set up shortcuts for the command line arguments + hdf_file_path = str(arguments['hdf_file_path']) + sample_id = int(arguments['sample_id']) + seconds_before = float(arguments['seconds_before']) + seconds_after = float(arguments['seconds_after']) + plot_path = str(arguments['plot_path']) + + # ------------------------------------------------------------------------- + # Read in the sample file + # ------------------------------------------------------------------------- + + print('Reading in HDF file...', end=' ') + data = SampleFile() + data.read_hdf(hdf_file_path) + df = data.as_dataframe(injection_parameters=True, + static_arguments=True) + print('Done!') + + # ------------------------------------------------------------------------- + # Plot the desired sample + # ------------------------------------------------------------------------- + + print('Plotting sample...', end=' ') + + # Select the sample (i.e., the row from the data frame of samples) + try: + sample = df.loc[sample_id] + except KeyError: + raise KeyError('Given sample_id is too big! Maximum value = {}'. + format(len(df) - 1)) + + # Check if the sample we have received contains an injection or not + if 'h1_signal' in sample.keys(): + has_injection = isinstance(sample['h1_signal'], np.ndarray) + else: + has_injection = False + + # Read out and construct some necessary values for plotting + seconds_before_event = float(sample['seconds_before_event']) + seconds_after_event = float(sample['seconds_after_event']) + target_sampling_rate = float(sample['target_sampling_rate']) + sample_length = float(sample['sample_length']) + + # Create a grid on which the sample can be plotted so that the + # event_time is at position 0 + grid = np.linspace(0 - seconds_before_event, 0 + seconds_after_event, + int(target_sampling_rate * sample_length)) + + # Create subplots for H1 and L1 + fig, axes1 = plt.subplots(nrows=2) + + # If the sample has an injection, we need a second y-axis to plot the + # pure (i.e., unwhitened) detector signals + if has_injection: + axes2 = [ax.twinx() for ax in axes1] + else: + axes2 = None + + # Plot the strains for H1 and L1 + for i, (det_name, det_string) in enumerate([('H1', 'h1_strain'), + ('L1', 'l1_strain')]): + + axes1[i].plot(grid, sample[det_string], color='C0') + axes1[i].set_xlim(-seconds_before, seconds_after) + axes1[i].set_ylim(-150, 150) + axes1[i].tick_params('y', colors='C0', labelsize=8) + axes1[i].set_ylabel('Amplitude of Whitened Strain ({})' + .format(det_name), color='C0', fontsize=8) + + # If applicable, also plot the detector signals for H1 and L1 + if has_injection: + + # Get the maximum value of the detector signal (for norming them) + maximum = max(np.max(sample['h1_signal']), np.max(sample['l1_signal'])) + + for i, (det_name, det_string) in enumerate([('H1', 'h1_signal'), + ('L1', 'l1_signal')]): + axes2[i].plot(grid, sample[det_string] / maximum, color='C1') + axes2[i].set_xlim(-seconds_before, seconds_after) + axes2[i].set_ylim(-1.2, 1.2) + axes2[i].tick_params('y', colors='C1', labelsize=8) + axes2[i].set_ylabel('Rescaled Amplitude of Simulated\n' + 'Detector Signal ({})'.format(det_name), + color='C1', fontsize=8) + + # Also add the injection parameters + if has_injection: + keys = ('mass1', 'mass2', 'spin1z', 'spin2z', 'ra', 'dec', + 'coa_phase', 'inclination', 'polarization', 'injection_snr') + string = ', '.join(['{} = {:.2f}'.format(_, float(sample[_])) + for _ in keys]) + else: + string = '(sample does not contain an injection)' + plt.figtext(0.5, 0.9, 'Injection Parameters:\n' + string, + fontsize=8, ha='center') + + # Add a vertical line at the position of the event (x=0) + axes1[0].axvline(x=0, color='black', ls='--', lw=1) + axes1[1].axvline(x=0, color='black', ls='--', lw=1) + + # Set x-labels + axes1[0].set_xticklabels([]) + axes1[1].set_xlabel('Time from event time (in seconds)') + + # Adjust the size and spacing of the subplots + plt.gcf().set_size_inches(12, 6, forward=True) + plt.tight_layout(rect=[0, 0, 1, 0.9]) + plt.subplots_adjust(wspace=0, hspace=0) + + # Add a title + plt.suptitle('Sample #{}'.format(sample_id), y=0.975) + + # Save the plot at the given location + plt.savefig(plot_path, bbox_inches='tight', pad_inches=0) + + print('Done!') + + # ------------------------------------------------------------------------- + # Postliminaries + # ------------------------------------------------------------------------- + + # Print the total run time + print('') + print('Total runtime: {:.1f} seconds!' + .format(time.time() - script_start_time)) + print('') diff --git a/gen_samples_py3/plots/0.png b/gen_samples_py3/plots/0.png new file mode 100644 index 0000000000000000000000000000000000000000..9aca8236ebbbb3416adeffee4726631e266c9d7b Binary files /dev/null and b/gen_samples_py3/plots/0.png differ diff --git a/gen_samples_py3/plots/1.png b/gen_samples_py3/plots/1.png new file mode 100644 index 0000000000000000000000000000000000000000..3669dc758f090e10f1212833a1215ca2100bf915 Binary files /dev/null and b/gen_samples_py3/plots/1.png differ diff --git a/gen_samples_py3/plots/100.png b/gen_samples_py3/plots/100.png new file mode 100644 index 0000000000000000000000000000000000000000..c3364890b2a350722c6665e369035ad5c3539392 Binary files /dev/null and b/gen_samples_py3/plots/100.png differ diff --git a/gen_samples_py3/plots/101.png b/gen_samples_py3/plots/101.png new file mode 100644 index 0000000000000000000000000000000000000000..4a3c60da55cbbc7a469fb207952c4333c9f519d0 Binary files /dev/null and b/gen_samples_py3/plots/101.png differ diff --git a/gen_samples_py3/plots/102.png b/gen_samples_py3/plots/102.png new file mode 100644 index 0000000000000000000000000000000000000000..952f21e5aaf0f6bb774b8dc05ed7ffc906276dea Binary files /dev/null and b/gen_samples_py3/plots/102.png differ diff --git a/gen_samples_py3/plots/103.png b/gen_samples_py3/plots/103.png new file mode 100644 index 0000000000000000000000000000000000000000..f99caf5446b03427b4a512d2d7379022186b619d Binary files /dev/null and b/gen_samples_py3/plots/103.png differ diff --git a/gen_samples_py3/plots/104.png b/gen_samples_py3/plots/104.png new file mode 100644 index 0000000000000000000000000000000000000000..07ac49c4928e0d3bbd4402927b9ae7f7e2dc27d4 Binary files /dev/null and b/gen_samples_py3/plots/104.png differ diff --git a/gen_samples_py3/plots/105.png b/gen_samples_py3/plots/105.png new file mode 100644 index 0000000000000000000000000000000000000000..10b110eba71fe605295612edec6aaf592eb7d4d2 Binary files /dev/null and b/gen_samples_py3/plots/105.png differ diff --git a/gen_samples_py3/plots/106.png b/gen_samples_py3/plots/106.png new file mode 100644 index 0000000000000000000000000000000000000000..3476bf4e4a9fd66a87fdeb1ce115ce06824084a6 Binary files /dev/null and b/gen_samples_py3/plots/106.png differ diff --git a/gen_samples_py3/plots/107.png b/gen_samples_py3/plots/107.png new file mode 100644 index 0000000000000000000000000000000000000000..a7c762ef949d6812a87bfcc26d327eeab41d781d Binary files /dev/null and b/gen_samples_py3/plots/107.png differ diff --git a/gen_samples_py3/plots/108.png b/gen_samples_py3/plots/108.png new file mode 100644 index 0000000000000000000000000000000000000000..e27678b829fb93a627e0846646acdfb0bf2dc67b Binary files /dev/null and b/gen_samples_py3/plots/108.png differ diff --git a/gen_samples_py3/plots/109.png b/gen_samples_py3/plots/109.png new file mode 100644 index 0000000000000000000000000000000000000000..68fea19a4812d761548584c6ae054abcae19516a Binary files /dev/null and b/gen_samples_py3/plots/109.png differ diff --git a/gen_samples_py3/plots/2.png b/gen_samples_py3/plots/2.png new file mode 100644 index 0000000000000000000000000000000000000000..bb012b4ba99e299e884947aebc039af6673f85d7 Binary files /dev/null and b/gen_samples_py3/plots/2.png differ diff --git a/gen_samples_py3/plots/3.png b/gen_samples_py3/plots/3.png new file mode 100644 index 0000000000000000000000000000000000000000..52a42c84886f57b8872e92df7c79729668d43620 Binary files /dev/null and b/gen_samples_py3/plots/3.png differ diff --git a/gen_samples_py3/plots/4.png b/gen_samples_py3/plots/4.png new file mode 100644 index 0000000000000000000000000000000000000000..c29128f97a4697c4f9589039e6189408687629a3 Binary files /dev/null and b/gen_samples_py3/plots/4.png differ diff --git a/gen_samples_py3/plots/5.png b/gen_samples_py3/plots/5.png new file mode 100644 index 0000000000000000000000000000000000000000..c150fe76cf172f5c7e3f8faa49616260adfd7e28 Binary files /dev/null and b/gen_samples_py3/plots/5.png differ diff --git a/gen_samples_py3/plots/6.png b/gen_samples_py3/plots/6.png new file mode 100644 index 0000000000000000000000000000000000000000..5d1b18f69429f8b6cf3504ecdebb3de43eafb0ad Binary files /dev/null and b/gen_samples_py3/plots/6.png differ diff --git a/gen_samples_py3/plots/7.png b/gen_samples_py3/plots/7.png new file mode 100644 index 0000000000000000000000000000000000000000..568f30da484172be95614ad8fca486af0817eeaf Binary files /dev/null and b/gen_samples_py3/plots/7.png differ diff --git a/gen_samples_py3/plots/8.png b/gen_samples_py3/plots/8.png new file mode 100644 index 0000000000000000000000000000000000000000..fcfcc7a9139045709b76a38a94f11999e7951069 Binary files /dev/null and b/gen_samples_py3/plots/8.png differ diff --git a/gen_samples_py3/plots/9.png b/gen_samples_py3/plots/9.png new file mode 100644 index 0000000000000000000000000000000000000000..8c62fea54283e8d918fce97db47a6353f93b5643 Binary files /dev/null and b/gen_samples_py3/plots/9.png differ diff --git a/gen_samples_py3/plots/signalnoise_id0.png b/gen_samples_py3/plots/signalnoise_id0.png new file mode 100644 index 0000000000000000000000000000000000000000..dfa8aa83b488aad66c3226052803336f797c180b Binary files /dev/null and b/gen_samples_py3/plots/signalnoise_id0.png differ diff --git a/python3_samples/utils/__init__.py b/gen_samples_py3/utils/__init__.py similarity index 100% rename from python3_samples/utils/__init__.py rename to gen_samples_py3/utils/__init__.py diff --git a/python3_samples/utils/configfiles.py b/gen_samples_py3/utils/configfiles.py similarity index 100% rename from python3_samples/utils/configfiles.py rename to gen_samples_py3/utils/configfiles.py diff --git a/python3_samples/utils/hdffiles.py b/gen_samples_py3/utils/hdffiles.py similarity index 100% rename from python3_samples/utils/hdffiles.py rename to gen_samples_py3/utils/hdffiles.py diff --git a/python3_samples/utils/progressbar.py b/gen_samples_py3/utils/progressbar.py similarity index 100% rename from python3_samples/utils/progressbar.py rename to gen_samples_py3/utils/progressbar.py diff --git a/python3_samples/utils/samplefiles.py b/gen_samples_py3/utils/samplefiles.py similarity index 99% rename from python3_samples/utils/samplefiles.py rename to gen_samples_py3/utils/samplefiles.py index 4a6de14b4e1f5b9ce9976a607c947355df18b12b..526802117d3c87617461b277ac98338a7b475fb4 100755 --- a/python3_samples/utils/samplefiles.py +++ b/gen_samples_py3/utils/samplefiles.py @@ -139,15 +139,16 @@ class SampleFile: # Read in dict with command_line_arguments self.data['command_line_arguments'] = \ dict(hdf_file['command_line_arguments'].attrs) + self.data['command_line_arguments'] = \ - {key: value.decode('ascii') for key, value in + {key: value for key, value in iteritems(self.data['command_line_arguments'])} # Read in dict with static_arguments self.data['static_arguments'] = \ dict(hdf_file['static_arguments'].attrs) self.data['static_arguments'] = \ - {key: value.decode('ascii') for key, value in + {key: value for key, value in iteritems(self.data['static_arguments'])} # Read in group containing injection samples diff --git a/python3_samples/utils/samplegeneration.py b/gen_samples_py3/utils/samplegeneration.py similarity index 85% rename from python3_samples/utils/samplegeneration.py rename to gen_samples_py3/utils/samplegeneration.py index 76be58a1245be67e48e3105c457e8f1dc663de98..31751875dd89ba6257a82242b073af1541351045 100755 --- a/python3_samples/utils/samplegeneration.py +++ b/gen_samples_py3/utils/samplegeneration.py @@ -12,7 +12,7 @@ from __future__ import print_function import numpy as np from lal import LIGOTimeGPS -from pycbc.psd import interpolate +from pycbc.psd import interpolate, inverse_spectrum_truncation from pycbc.psd.analytical import aLIGOZeroDetHighPower from pycbc.noise import noise_from_psd from pycbc.filter import sigma @@ -25,6 +25,22 @@ from .waveforms import get_detector_signals, get_waveform # FUNCTION DEFINITIONS # ----------------------------------------------------------------------------- +def signal_whiten(psd, signal, segment_duration, max_filter_duration, + trunc_method='hann', low_frequency_cutoff=None): + # Estimate the noise spectrum, no need for segment_duration, because we already get in line 193. + psd = interpolate(psd, signal.delta_f) + max_filter_len = int(max_filter_duration * signal.sample_rate) + + # Interpolate and smooth to the desired corruption length + psd = inverse_spectrum_truncation(psd, + max_filter_len=max_filter_len, + low_frequency_cutoff=low_frequency_cutoff, + trunc_method=trunc_method) + + # Whiten the data by the asd + white = (signal.to_frequencyseries() / psd**0.5).to_timeseries() + return white + def generate_sample(static_arguments, event_tuple, waveform_params=None): @@ -133,6 +149,7 @@ def generate_sample(static_arguments, if waveform_params is None: detector_signals = None injection_parameters = None + output_signals = None strain = noise # Otherwise, we need to simulate a waveform for the given waveform_params @@ -154,6 +171,10 @@ def generate_sample(static_arguments, event_time=event_time, waveform=waveform) + # Store the output_signal + output_signals = {} + output_signals = detector_signals.copy() + # --------------------------------------------------------------------- # Add the waveform into the noise as is to calculate the NOMF-SNR # --------------------------------------------------------------------- @@ -199,7 +220,7 @@ def generate_sample(static_arguments, # injection SNR strain[det] = noise[det].add_into(scale_factor * detector_signals[det]) - + output_signals[det] = scale_factor * output_signals[det] # --------------------------------------------------------------------- # Store some information about the injection we just made # --------------------------------------------------------------------- @@ -210,7 +231,7 @@ def generate_sample(static_arguments, 'l1_snr': snrs['L1']} # Also add the waveform parameters we have sampled - for key, value in waveform_params.items(): + for key, value in iter(waveform_params.items()): injection_parameters[key] = value # ------------------------------------------------------------------------- @@ -231,6 +252,13 @@ def generate_sample(static_arguments, max_filter_duration=max_filter_duration, remove_corrupted=False) + if waveform_params is not None: + output_signals[det] = \ + signal_whiten(psd = psds[det], + signal = output_signals[det], + segment_duration = segment_duration, + max_filter_duration = max_filter_duration) + # Get the limits for the bandpass bandpass_lower = static_arguments['bandpass_lower'] bandpass_upper = static_arguments['bandpass_upper'] @@ -241,14 +269,20 @@ def generate_sample(static_arguments, strain[det] = strain[det].highpass_fir(frequency=bandpass_lower, remove_corrupted=False, order=512) - + if waveform_params is not None: + output_signals[det] = output_signals[det].highpass_fir(frequency=bandpass_lower, + remove_corrupted=False, + order=512) # Apply a low-pass filter to remove everything above `bandpass_upper`. # If bandpass_upper = sampling rate, do not apply any low-pass filter. if bandpass_upper != target_sampling_rate: strain[det] = strain[det].lowpass_fir(frequency=bandpass_upper, remove_corrupted=False, order=512) - + if waveform_params is not None: + output_signals[det] = output_signals[det].lowpass_fir(frequency=bandpass_upper, + remove_corrupted=False, + order=512) # ------------------------------------------------------------------------- # Cut strain (and signal) time series to the pre-specified length # ------------------------------------------------------------------------- @@ -267,6 +301,7 @@ def generate_sample(static_arguments, # Cut the detector signals to the specified length detector_signals[det] = detector_signals[det].time_slice(a, b) + output_signals[det] = output_signals[det].time_slice(a, b) # Also add the detector signals to the injection parameters injection_parameters['h1_signal'] = \ @@ -274,6 +309,11 @@ def generate_sample(static_arguments, injection_parameters['l1_signal'] = \ np.array(detector_signals['L1']) + injection_parameters['h1_output_signal'] = \ + np.array(output_signals['H1']) + injection_parameters['l1_output_signal'] = \ + np.array(output_signals['L1']) + # ------------------------------------------------------------------------- # Collect all available information about this sample and return results # ------------------------------------------------------------------------- diff --git a/python3_samples/utils/staticargs.py b/gen_samples_py3/utils/staticargs.py similarity index 100% rename from python3_samples/utils/staticargs.py rename to gen_samples_py3/utils/staticargs.py diff --git a/python3_samples/utils/waveforms.py b/gen_samples_py3/utils/waveforms.py similarity index 100% rename from python3_samples/utils/waveforms.py rename to gen_samples_py3/utils/waveforms.py diff --git a/python3_samples/.DS_Store b/python3_samples/.DS_Store deleted file mode 100644 index 749e222af6d2280dd99e24f4c7942a3c7586b65e..0000000000000000000000000000000000000000 Binary files a/python3_samples/.DS_Store and /dev/null differ