{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Correlated analyses with MultiTrialRunner\n", "\n", ".. contents:: :local:\n", "\n", "In this tutorial, we consider the use case of multiple, potentially partially-correlated analyses. Here, we want to perform multiple likelihood tests on the same \"trials\" — the same realizations of the datasets." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import histlite as hl\n", "import csky as cy\n", "\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/mnt/lfs7/user/ssclafani/software/external/csky/csky/plotting.py:92: MatplotlibDeprecationWarning: Support for setting the 'text.latex.preamble' or 'pgf.preamble' rcParam to a list of strings is deprecated since 3.3 and will be removed two minor releases later; set it to a single string instead.\n", " r'\\SetSymbolFont{operators} {sans}{OT1}{cmss} {m}{n}'\n" ] } ], "source": [ "cy.plotting.mrichman_mpl()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "ana_dir = '/data/user/mrichman/csky_cache/ana'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll test with one year of data, PS IC86-2011." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Setting up Analysis for:\n", "IC86_2011\n", "Setting up IC86_2011...\n", "Reading /data/ana/analyses/ps_tracks/version-003-p02/IC86_2011_MC.npy ...\n", "Reading /data/ana/analyses/ps_tracks/version-003-p02/IC86_2011_exp.npy ...\n", "Reading /data/ana/analyses/ps_tracks/version-003-p02/GRL/IC86_2011_exp.npy ...\n", "<- /data/user/mrichman/csky_cache/ana/IC86_2011.subanalysis.version-003-p02.npy \n", "Done.\n" ] } ], "source": [ "ana = cy.get_analysis(cy.selections.repo, 'version-003-p02' , cy.selections.PSDataSpecs.ps_2011, dir=ana_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Everything that follows will use this same `Analysis`, so we put it in `cy.CONF`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "cy.CONF['ana'] = ana\n", "cy.CONF['mp_cpus'] = 15" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Individual Trial Runners" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We begin by setting up an arbitrary catalog of two sources:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
decdec_degsindecrara_degweightextensionextension_deg
0-0.785398-45.0-0.7071070.0000000.00.50.00.0
10.78539845.00.7071073.141593180.00.50.00.0
\n", "
" ], "text/plain": [ " dec dec_deg sindec ra ra_deg weight extension \\\n", "0 -0.785398 -45.0 -0.707107 0.000000 0.0 0.5 0.0 \n", "1 0.785398 45.0 0.707107 3.141593 180.0 0.5 0.0 \n", "\n", " extension_deg \n", "0 0.0 \n", "1 0.0 " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ra = np.r_[0, 180]\n", "dec = np.r_[-45, 45]\n", "src = cy.sources(ra, dec, deg=True)\n", "src.as_dataframe" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We get a trial runner for each:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[,\n", " ]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trs = [cy.get_trial_runner(src=s) for s in src]\n", "trs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now for the tricky part: we need to make sure our background trials have events for every source, but we don't want to scramble the whole sky. So we set up our \"injection\" trial runner to know about the whole catalog.\n", "\n", "We also need to decide on how the signal injection should work; here we'll just inject signals for the 0th source." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "tr_inj = cy.get_trial_runner(src=src, inj_conf=dict(src=src[0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To confirm that signal events wind up where they should, let's look at an example trial:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Trial(evss=[[Events(13188 items | columns: dec, idx, inj, log10energy, ra, sigma, sindec), Events(5 items | columns: dec, idx, inj, log10energy, ra, sigma, sindec)]], n_exs=[123056])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trial = tr_inj.get_one_trial(n_sig=5, poisson=False, seed=2)\n", "trial" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
declog10energyrasigmasindecidxinj
0-0.7891575.5772820.0007940.003491-0.704913773766<csky.inj.PointSourceInjector object at 0x7f91...
1-0.8315545.1182370.0082260.011637-0.736841311727<csky.inj.PointSourceInjector object at 0x7f91...
2-0.7870045.2416690.0038330.003491-0.720210698751<csky.inj.PointSourceInjector object at 0x7f91...
3-0.7890605.3979880.0016290.005621-0.702275499808<csky.inj.PointSourceInjector object at 0x7f91...
4-0.7733555.2758706.2791550.018069-0.697965746488<csky.inj.PointSourceInjector object at 0x7f91...
\n", "
" ], "text/plain": [ " dec log10energy ra sigma sindec idx \\\n", "0 -0.789157 5.577282 0.000794 0.003491 -0.704913 773766 \n", "1 -0.831554 5.118237 0.008226 0.011637 -0.736841 311727 \n", "2 -0.787004 5.241669 0.003833 0.003491 -0.720210 698751 \n", "3 -0.789060 5.397988 0.001629 0.005621 -0.702275 499808 \n", "4 -0.773355 5.275870 6.279155 0.018069 -0.697965 746488 \n", "\n", " inj \n", "0 \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
gamma_0000gamma_0001mlog10p_0000mlog10p_0001ns_0000ns_0001ts_0000ts_0001
01.9726293.25000010.4003460.26760614.9347830.00000040.8149300.000000
12.0284844.0000007.7199880.26760613.1356720.00000029.3091750.000000
22.0759802.8425994.3509880.3355038.7859780.91964215.0359110.046768
31.9239784.0000009.2931340.26760613.9247850.00000036.0522410.000000
42.0259772.0000005.0807750.2676067.9435360.00000018.0980950.000000
\n", "" ], "text/plain": [ " gamma_0000 gamma_0001 mlog10p_0000 mlog10p_0001 ns_0000 ns_0001 \\\n", "0 1.972629 3.250000 10.400346 0.267606 14.934783 0.000000 \n", "1 2.028484 4.000000 7.719988 0.267606 13.135672 0.000000 \n", "2 2.075980 2.842599 4.350988 0.335503 8.785978 0.919642 \n", "3 1.923978 4.000000 9.293134 0.267606 13.924785 0.000000 \n", "4 2.025977 2.000000 5.080775 0.267606 7.943536 0.000000 \n", "\n", " ts_0000 ts_0001 \n", "0 40.814930 0.000000 \n", "1 29.309175 0.000000 \n", "2 15.035911 0.046768 \n", "3 36.052241 0.000000 \n", "4 18.098095 0.000000 " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trials.as_dataframe.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once again, since we've injected plenty of signal (10 events per trial, but now with Poisson smearing), we see that the results are much more significant for source `0000` than for source `0001` (which only has background nearby)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the resulting distributions like so:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7MAAALDCAYAAADHS1caAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAASdAAAEnQB3mYfeAAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdf4wbdZ7/+ZdJvt0HyQS3uycjRpCQMkH7PXIs407E8P3qq1lNuyesOHYyg50Idpj7YwcbCCM02tku2K/u1Oi+WnDPH6vRBth2RrfSwsEm9gArbr+npd0rcV9pGUS6N/Pd8L0RpCskMNwXknbXDAls9yZT90ePjd2//KvsqrKfD8mCKtvlD1Sct19Vnx8hx3EcAQAAAAAQIFd53QAAAAAAABpFmAUAAAAABA5hFgAAAAAQOIRZAAAAAEDgEGYBAAAAAIFDmAUAAAAABA5hFgAAAAAQOIRZAAAAAEDgEGYBAAAAAIFDmAUAAGuamJjQ6OiohoeHZZqm180BAKDKZq8bAAAA/Me2bUnS1NSUJGlgYECHDh1SLBbzslkAAJRxZxYAgC43OzuraDS65nPZbFYDAwMaGBhQOp0u7w+HwxobG5MkWZYl27ZlGEZH2gsAQD0IswAAdKmJiQmFQiElk0lZlrXq+Xw+L9M0NT09rZmZGZ04cULJZLLqNclkUtFoVLlcTuFwuFNNBwCgppDjOI7XjfAL27b1+uuv64YbblB/f7/XzQEA+NTi4qLef/99fe1rXwtEwMvn80omk1pZ8qPRqEzTVCqVkrR8BzYajWphYaH832XbtizLUjKZVC6Xa7ibMbUVAFCvhuurg7JXXnnFkcSDBw8ePHjU9XjllVe8Ll11yeVyjlRd8hcWFhxJztzcXNX+cDjs5HK5VcdIJBJOKpVq+LOprTx48ODBo9FHvfWVCaAq3HDDDZKkV155RTfddJPHrQEA+NXp06d14MCBct0IolK345XjYA3DkGVZyufzMgyjfCd2dna25ozGH3/8sc6fP1+1b2lpSRK1FQBQW6P1lTBbodT96aabbtItt9zicWsAAH4X5G6zxWJxzf2RSETz8/MaGxuTaZoyTVPFYlGJRKLcHXk9zzzzjJ544ok1n6O2AgDqVW99JcwCAIA1ZTKZhl7/8MMPr5pAqnSVHQAAtxFmJY2Pj697JRkAgG4UiUTW3F8sFjU4ONjUMbdv367t27e30iwAAOrG0jxaDrOO4+jUqVNeNwUAgI4ojZVduWSPZVmBmKEZAADCLAAAPSgcDsswDOXz+fI+y7Jk27YOHjzY8vHHx8cVCoW0Z8+elo8FAMBaCLMAAPSo0gRPs7Oz5bVkE4mEK3dm6fUEAGg3xswCANClstms0ul0eTsUCkmSHMeRJKVSKdm2rZGREdm2rUQioVwu50lbAQBoFHdmRVcoAEB3SqVSchxn1aPS2NiYFhYW5DgOQRYAECiEWdEVCgAAAACChjALAABcR68nAEC7EWYBAIDr2tHr6cLFRb3z0Sd656NPdOHiomvHBQAEExNAtcGFi4uafH1On/3rFUnSd++4UTd/6QsetwoAgGCbfH1OR//LGUnS7bsievoPYxra2u9xqwAAXiHMtkHx0lK52EqSHOk/fet/8q5BAAB0gc+WrpT//c0zRRUvLRFmAaCH0c0YAAAAABA4hFm5P0lFZEufbt8VKW9f3bfJleMCABAU7ZgA6rv/7kZ956s79J2v7tAD/2GXIlv6XDs2ACB4Qs7KBed62Ntvv609e/bo1KlTuuWWW1o61oWLiypeWpK0HG7pBgUA3cPNetHt+H8FAKhXozWDMbNtMrS1nwALAAAAAG1CN2MAAAAAQOAQZgEAAAAAgUOYBQAAAAAEDmEWAAC4rh2zGQMAUIkwKwouAABuGx8fl+M4OnXqlNdNAQB0KcKsKLgAAAAAEDSEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAA4DpWCgAAtBthFgAAuI6VAgAA7UaYBQAAAAAEDmFWdIUCAAAAgKAhzIquUAAAAAAQNIRZAAAAAEDgEGYBAAAAAIFDmAUAAAAABA5hFgAAAAAQOIRZAADgOlYKAAC0G2EWAAC4jpUCAADtRpgFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmxSQVAAAAABA0hFkxSQUAAAAABA1hFgAAAAAQOIRZAAAAAEDgEGYBAAAAAIFDmAUAAAAABA5hFgAAuI6VAgAA7UaYBQAArmOlAABAuxFmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmtbyweygU0p49e7xuCgAAAACgDoRZLYdZx3F06tQpr5sCAAAAAKgDYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAALiOlQIAAO1GmAUAAK5jpQAAQLsRZgEAAAAAgUOYBQAAAAAEDmEWAAAAABA4hFkAAAAAQOAQZgEAAAAAgUOYBQAAAAAEDmEWAAAAABA4hFkAAAAAQOAQZgEAAAAAgUOYBQAAAAAEDmEWAAAAABA4hFkAAAAAQOB0RZg1TVPDw8OKRqOamJjwujkAAAAAgDbb7HUDWpXP52XbtmZmZiRJ0WhUsVhM8Xjc45YBAAAAANrFV3dmZ2dnFY1G13wum81qYGBAAwMDSqfT5f3xeFyZTKZqe3Z2tu1tBQAAAAB4xxdhdmJiQqFQSMlkUpZlrXo+n8/LNE1NT09rZmZGJ06cUDKZlCSFw2GFw+HyawuFAndlUdPl+Xmd/4sj5cfl+XmvmwQAAACgAb7oZjw2NqaxsTHl8/lySK1kmqYymYxisZgkKZfLKRqNyrbtqiBrmqYSiUT5dcB6rhSLuvD00+XtbXfu1+bBwbZ81uX5eS288GJ5e+C+e9v2WQAAAECv8EWY3Yht27Isq+puq2EYCofDKhQKSiQSkpaDbDQaVSqV8qqp8KnL8/O6UixW7Vs6e27DbUnaFIm4Ejo7GZwBwE2maapQKMi2baXTaY2NjXndJAAAynwfZkvdjg3DqNpvGEb5uXQ6rWQySfdirGnhhRerwuRaPnjkkVX7hg4f1he/v3r/RrwOzgDgFiZYBAD4ne/DbHFFMCiJRCKan59XPp9XNptVoVAoP1fP1eOPP/5Y58+fr9p3+vTp1huMnuZmcKZ7MgC3zM7OKplMam5ubtVz2WxWpmlKkg4ePKjJyUlJyxMqVgbX0gSLhFkAgF/4PszWkkgk5DhOw+975pln9MQTT7ShRdUuXFzUc2+cLW/ff8dODW3tb/vnIvjongygVRMTEzJNs6o3U6XKCRbD4bCSyaSSyaRyuVzVnBTS8gSLlasJoLdxwRWAH/g+zEYikTX3F4tFDbbwl+bDDz+8arKp06dP68CBA00fcy3FS0v68fS75e27br2OMNthA/fdq2137q/at3T2XNUd0uuPHFHfzh1Vr9m0zp+9dvjNp59q8d13q/bRPRlAq5hgEW5YbwhN5QXX/+Hf/ts16yg1CkA7+T7MlsbKWpZVNW7WsqxVV40bsX37dm3fvr3l9sH/Ng8O1iymfTt3qH/37pY/q9ngbL/8sqy7/2DDY7s1rhcAJCZYRP06OfcEADTC92E2HA7LMAzl8/nyOFjLsmTbtg4ePOhx6xBUmyIRDR0+XLXthmaD81VXX+PK5wNAvdoxwSLzUaCWTndPpjs00N18H2al5avC6XRa8Xi8PKYnkUi0dGe20vj4eEfGz8I/Ng8OcrUYQE9rxwSLnZqPAsHVzvkg6A4N9B5fhNlsNls1qUQoFJKk8sROqVRKtm1rZGREtm0rkUgol8u59vnj4+MaHx/X22+/rT179rh2XHSfZq7w1nMXOAjjegH0lmYmWOzUfBTorGZrlBMKdXQ+CLpDA73HF2E2lUrVHItTmsQC8FIzV5TruQvcyXG9ACC1Z4JF5qPoTs3WqPN/cYRwCaCtfBFmAT9ar7vSRttudlVqdlwv44MA1KNdEywiWJqtGe2aewIAGkGYRaC1M7g1013JzavJzY7rZX1aAPUI8gSLrOHunmZrhh/nnmDIDtB7CLNq7wRQkS19enRkudvNp0uXlTvxvq7pW/7fTvFtXa8Ht2buHktMdgFgWTsnWGxnbWUN9+Z0uma4GS7ruXjNkB2g9xBm1d4JoIa29usHozdLkt756BN948//n/JzFN/GENxWY7ILABvxcoLFTk6uuPDpUluP3y06XTPcDJe9fvEawNoIswiMThfhZq4oB7mrEmNtge7TKxMsvnryQ92+i7+vuoWbF68Z2wt0N9fC7MmTJ9fcf9ttt7n1EYHHleNg6bXuSlz1BvyH2opuU0+4dPPitR/H9gJwT8th9ujRo3rwwQcladV6dKFQSFeuXGn1I7rGqyc/9LoJ6DJBWfsPQGO6oba2c8wsmuOHCZIIlwDc1HKYNU1TDzzwgEzTXHfNOr+j4AaDH4qw37orsfYf0J26pbZ2aszs3bd9ua3H7xa91uMIQPdzpZvxY489phtvvNGNQ3mikwUXzfNDEeaKsjsYnwvUFvTa2kkD1/R53YTA8ttFWskfF6+pU0AwtBxmU6mUfvrTn+qP//iP3WhPV7v7ti/r+TdXd90E3OTHHyYrMT4X2Bi1FZ3ix4u0frh4TZ0CgsGVO7N/8id/oj/7sz/T3r17V61Nd+zYMTc+oitw5Xh9zV4BDUJw67R6fpj44ao3gI1RW9dXuYZ7aRvdrZ31nqX/gOBqOcxalqVEIlHeXjlRBVCPZq+A+vGKchB08qo3PxKAxnVDbW3nfBSVa7ijN7Sz3rs5ezLdk4HOajnMHj9+3I129ASuJKPb1FO0O70+MNANuqG2Mh8FehHdk4HOcm2d2SDr1GzGXElexp06f2qmCxdFGwDQq37z6acscwd4rOEwu3//fiWTSX3ve9+TJD300EMbvv7ZZ59trmUdxNXjzuJOnT/RZRvwTjfWVnQW3Vub1+w8EvbLL8u6+w82PDbdk4H2ajjMzs/Pb7gNoHutvKtezxXoL/z+nUw2BdRAbUWr6CnTvGbnkbjq6mtcawPnD2hOw2H2xIkTVdvdMK4HQH1q3VVv9o56u5dYAPyO2gqgWdzVRS9jzCw6jmVhAKD7dWo+il7E3BPtV888Es3+nnFCoabG2pbasfIcclcXvcyVMPvSSy9pcnKyfGV57969mpiY0O/+7u+6cfiuc+Hiop5742x5+/47dmpoa7+HLeosPyyGDgB+F/TaynwU7cPcE+1XzzwSzf6eOf8XR5o6fxLnEFip5TD7ox/9SKZpKpFI6KmnnpK03F0qFospn8/rW9/6VsuN7DbFS0v68fTnV+TuuvW6ngqzCK6VV6GbvaPezMzJQC+htgJYCzMoA9VaDrNPPvmkJiYm9MMf/rC874EHHtDo6KjGxsYouKgL4SYYal2FrveOOjMnAxujtgJYy2enTqn4V3+14Wu4K49e4ko340Qisea+xx57zI3Dt53X43oWPl3y7LP9gnADANWCXlvRnHom82HuCf9q51hbaXk5oM/eesvFFgPB1nKYffzxx1UoFMpr45W89957Gh4ebvXwHeH1uJ5XT36o23fR9QPBwx11oD26obZ2SrfNQ1HPZD7MPeFf7RxrK7m7HBDQDRoOsysXcp+fn9dPf/pTTU1NVe0vFAoyDKO11gHwNe6oA+6gtjaPeSjQS7grD1RrOMyutZD7PffcI8dxqvaNjIw036oe8/MPbP351DuSgn9FWWK9M7QPf7bQrbqxtno9hCcI3Fxih54yvYG78kC1hsMsC7m7759/+Wv98y9/Lak7riiz3hnahT9b6FbdWFu9HsITBG4usUNPmWDjYgTQHFcmgAIAAPADv0yqSE8SNKKVixEEYfQywqzP/PUb72lwy/Kd2SB0OW6mixRrnQEA2uU//V//TV//nS9J8raO0pMEncJdefQywqzPPP+zz4Of112O67mq3EwXKdY6Qz3cHEsGoHcEaegOk/kAQGsIs+r8JBWRLX16dGR5YP78xUU9/+bqH+R+wFVleMnNsWQA4EdM5gMArSHMqvOTVAxt7dcPRm+WJL3z0Se+DbMAAKA2epIAgDcIs2hJM12k6B4FAOgm9CRBEDFJGbpBy2H2Jz/5iQ4ePKht27a50R4EDF2k0C6MJUMvo7b2LmamRacwnAzdoOUw+9RTTykUCumP/uiP3GgPPEIXKfgNF0rQy7qhtnZ6Popuwcy0AFA/V8JsOp1WPB7Xzp073WgTPEAXKQDwj26orZ2aj8IPkyrSkwQAvNFymD1x4oR27dolwzAUj8dlGEbV888++2yrHwEAQE+httbPD5Mq0pMEfkcPPHSrlsOsZVkyDKNcaOfn51tuFIKN8T4A0BpqKwA30QMP3arlMHv8+HE32tGz/NA9SnK3ixTjfdAuXChBr6C2AgBQmytL87z00kv6m7/5G508eVLvvPOOJOmhhx5SMpnU17/+dTc+omv5oXuURBcpBAMXStBLqK0AAGys5TD7ox/9SNlsVn/5l3+pb3zjG+X9IyMjymQyFFwAABpEbQ0+epLAT5ikDN3KldmMp6enddttt1XtHx4e1qFDh1o9PAAAPYfa2pzKoTulba/QkwR+Qg88dKuWw6zjOBoYGFi137Is7dq1q9XDd4Rf1sLzUxEGAHinG2qrFyqH7gAAul/LYfaBBx5QKpVSLpcr7ztz5owefPBBPfjgg60eviM6tRZeLX4rwnSRQre4PD+vhRdeLG8P3HcvSw3A17qhtgIA0G4th9lMJqN0Oq1wOCxJ2r17tyzLUiqV0g9/+MOWGwjv0EUK3eJKsVi1JMG2O/cTZuFr1FYAAGpzZTbjyclJmaapf/qnf5IkxWIxukEBANCCoNdWvwzhAbA2euChG7gSZiVVLe4OAABaF+Ta6pchPADWRg88dIOr3DjISy+9pP3792twcFCDg4Pav3+/fv7zn7tx6J534eKi/nzqnfLjwsVFr5sEAOgAaisAABtrOcw+9thjSiaTchxHTz31lJ566in95je/USwW08svv+xGG3ta8dKSfjz9bvlRvLTkdZMAAG1GbQUAoLaWuxlns1mNjY3pySefLO974IEHZJqmxsbG9K1vfavVjwCAul2en9eVYrFq39LZcxtuS8tjhZgUCn5BbQUAoLaWw2wkElE6nV61/0//9E919OjRVg8PAA1ZeOHFqpmL1/LBI6vHCA0dPszYIfgGtRUAgNpa7macSCQ0PT29av/09LQOHTrU6uEBAOg51FYAAGpr+c6sZVn60Y9+pNdee628z7ZtFQoFxePxqqJ77NixVj+u5y18yphZAOh21FYAAGpzZWmee+65R47jlLevvfZa3XPPPZJUtR+te/Xkh7p9F+P6gPUM3Hevtt25v2rf0tlzVV2Lrz9yRH07d1S9hvX14DfUVndcuLio5944W96+/46dGtraX/f7L8/Pa+GFF8vbA/fdy/h6APCJlsPs8ePH3WgH2ohCjF6yeXCw5p/vvp071L97d4daBDSO2uqe0qoAJXfdel1DYfZKsVg1Dn/bnfupoQDgE67cmYW/UYgBAADgBm6SwE8IswFz921f9roJAAAA6FHcJIGfEGYDZuCaPq+bAACA75XGys5fXKza/9dvvKfBLcvdjBsdPwsA8BfCrKTx8XE98cQTXjejynpFGEBnNNuNqtPvA7C2lWNlS57/2bnyvzc6fhYA4C+EWS2H2fHxcb399tvas2eP182RtH4RBtAZzXaj6vT7ALjn8vy8rhSLVfuWzp7bcFtang2d7yuwMS7aoh1aDrMvvfSSJOnb3/62JOmhhx7S8ePHZRiGcrmcbrzxxlY/AhVqdY+iEAOrbYpENHT4cNU24GfUVm8svPBi1UWltVQu81UydPiwvvj91fsBfI6LtmiHlsOsaZqanJyUJB09elTZbFbHjx/Xa6+9pmQyqbfeeqvlRuJztbpHUYiB1TYPDvLnG4HSDbXVj0N4ADSGmyTwu5bD7NzcnPbu3StJyuVySqVSuueee/SVr3xFu1nHEQCAhnVDbfXjEB4AjeEmCfyu5TBrGIYWFha0bds2FQoFPfjgg5KkX/3qVwqHwy03EACAXkNtbd3Cp0sNv2fgvnu17c79VfuWzp6r+rF+/ZEj6tu5o+o1DF0A2oeJFbGRlsNs6Urx4OCgDMMoj+8pFAoaGRlpuYFY9p2v7pAc6fk3V3flqEQhBoDgo7a27tWTHzb8ns2DgzV/7Pbt3KH+gNwdB7oBEytiIy2H2Uwmo0OHDunMmTOKx+Pl/YZhVG2jMSuvKH/3jhslVYfZtSaDohADzVk5LqieMUFOKKSQ42z4Ojffxxik3kFtBeAHzd4kcUIhLb777qr3bbQtUefQOFeW5onFYorFYlX77rnnHjcO3bPquaLMWnmAe2qNC1prTNDV+/bpsxoT8bj5PsYg9RZqKwCvNXuT5PxfHGGsLTqi4TD70EMPNfT6Z599ttGPAACgp1Bb3Xf3bV+uOTQHABBsDYfZ+fn5drQDAICeRW1138A1fV43AQDQZg2H2ePHj7ejHViBK8pAZ60cF1TvmKC1xr62631M1Na9qK0AugkTkqJTXBkzC/dxRRnorFrjgpqdOK3T7wMAwGvNjrW9PD/f1MRRTKzYu1wJsydPntSxY8dkWdaq544dO+bGR0DSZdv2ugkAgA6htrrr3t+5Vr/+z/9Zf2f8+4betykS0dDhw1XbANqj1mSMEhMrotpVrR7g6NGjisVimpqakuM4yufzchxHU1NTsglfrvrNr37ldRMAAB1AbXXfH/6PYd195h8bft/mwUF98fuPlB/cxQEA/2j5zuxjjz2mqamp8iLukUhEx48fV6FQ0NGjR1tuIFrHVWUACBZqKwC/4ncl/KTlMLuwsKB9+/aVtw3D0CeffKJ9+/bp0KFDrR4eLihdVQYABAO1FYBftfN3ZbMTRzGxYu9qOczGYjEVCgV9+9vfliTt3btXx48f18DAAF2hAABoArUVQC9qduKoejCxYndqOcxmMhlNTU2VC65pmopGowqFQkqlUi03EACAXkNtbc3l+Xlt/ehjPfKVz++6bL3w33Xt4iX94S9ekyRt+/3f19aPPtDirz//KcTMpgAQLC2H2ZGRkfKYHknatWuXFhYWVCwWtWvXrlYPDwBAz6G2tmbhhRf1q6ef1l0V+z6TFJb0nd+GWdv6Rz39f/+78vP/s/WPuumB/4VhOQAQIG1ZZ/baa6/Vtdde245Dr8uyLGUyGU1OTnb0c9slsqVPj47srtr+7x62B+g1zU5w0en3oXd4UVu72a/6t+j//J1vlLf/wy9/7mFrgO5HnUM7NBxm9+/fr2Qyqe9973uSpIceemjD1z/77LPNtawBo6OjKhaLinTRl2Joa79+MHpz1T7CLNA5zU5w0en3oTv4sbYCgJuoc2iHhsPs/Pz8htutmJ2dVTKZ1Nzc3KrnstmsTNOUJB08eLDqDuzU1JTy+XzX3JWVlsf7XCkWq/b964f/37qvXzp7Vou/7me8DwAEUDtra6uC2POpnhlRv/Qf/6P0s6Xy9vVPH9FA9LqOtREA0LqGw+yJEyeqto8fP95yIyYmJmSapgzDkGVZq57P5/MyTVPT09MKh8NKJpNKJpPK5XItf7ZfLbzwoi48/XTVvn/p26o/NJbH99h9W/R3xr8vP/fB4Ue06ZOPNHT4MFe9ACBg2lFb3RDUnk/1zIj6b758naSz5e2+nTu1efALbW4ZAMBNV7V6gJ/85Cf69a9/3dIxxsbG5DiOMpnMms+bpqlMJqNYLCbDMJTL5ZTP53tueYLw0kV95xev6Tu/eE13n/lHr5sDAGgTN2prpdnZWUWj0TWfy2azGhgY0MDAgNLpdNVzU1NTevzxx11rBwAAbmo5zD711FNtvUNq27Ysy1I8Hi/vMwxD4XBYhUKhbZ8LAIBX3KqtExMTCoVCSiaTNXs+zczM6MSJE0omky1/LgAAndDybMZPPfWU0um04vG4du7c6UabqpSKr2EYVfvX65LcLeoZ71Pp+qePyBjoZ2Y4AOgCbtXWsbExjY2NKZ/PrxlSK3s+SVIul1M0GpVt2wqHw01/LgC4hVUCsJGWw+yJEye0a9cuGYaheDy+KnS2OuNiccUkSCWRSKQ8QUY6ndaJEydkWZZGR0erCvN6Pv74Y50/f75q3+nTp1tqq5vqGe9TqW/nTvV/ibE+ANAN2l1bpdo9nxKJRMufESR//cZ7GtzSr/vv2Kmhrf1eNwfAb7FKADbScpi1LEuGYZQLrRczMDYzw+IzzzyjJ554og2tAQCgNZ2orb3a86nE/pcrVdvP/+ycJOmuW68jzAJAQLQcZts94+J6MygWi0UNtrAEzcMPP7yqy9Xp06d14MCBpo/pJxcuLuq5Nz6fpZErzQAQHJ2YzbgdPZ/83uup0t9Z7k2wBQDwRsth9h/+4R/09a9/fdX+kydPSpJuu+22lo5fumJcukpdYllWS+N5tm/fru3bt7fUNi9du3hJj3wlok2R5UAf2dJX9Xzx0pJ+PP1ueZsrzQAQHO2urfVqtOcTvZ4AAJ3UcpgdHR3VlStXVu2fn5/XxMSE/v7v/76l44fDYRmGoXw+r7GxMUnLQda2bR08eLClYwdZeOmivh8bUv/u3V43BQDgsnbXVqk9PZ+6qdcTPZwAwP9aDrOO46z73MpF4JtlmmZ5VsdwOKxkMqlEIuHaTIvj4+NcSQYA+EYnams7ej75udfTyplN/+C26/XiL/7buq+nhxMA+F/TYTYSiSgUCikUCq26gmvbthzH0fDwcF3HymazVQu1h0IhSZ8X81QqJdu2NTIyItu2lUgkXF3bdnx8XOPj43r77be1Z88e144LAEAj3KyttfRaz6eVM5sufPSJh60BALih6TA7PT0tx3G0d+/eNSeqMAxDu3btqutYqVRKqVRqw9eU1spDtXq7QZWWHNjoNQCC5/L8vBZeeLG8PXDfvTWX9WrmPegMN2trPdrZ84leTwCAdms6zH7lK1+RJMXjcY2MjLjWIDSm3m5QpSUHNnoNgOC5UizqwtNPl7e33bm/ZjBt5j3oDLdrq2FxVfwAACAASURBVJc9n+j1BABot6taPcBrr73mRjsAAMBvuVVbU6mUHMdZ9ag0NjamhYUFOY7j6hAeAADareUw2w3Gx8cVCoW4cgwAQI9b+HTJ6yYAAOpEmNVymHUcR6dOnfK6KRsqzcRYemxaZ1kFAADQnFdPfuh1EwAAdWp5aR50zsqZGCVJzMYIAPAhJoACALQbd2YBAIDrgtLrCQAQXIRZAACA37r7ti973QQAQJ3oZiy6QgEAgGUD1/RJWj0RVGm9dtZqB7pbs2uxd/p9WEaYFWvhAQCAaisngiqt185a7UB3a3Yt9k6/D8voZgwAAFzHsncAgHYjzAIAANf5fQKoyJY+PTqyW9+5fYfXTQEANIluxl2Gxd6B7nV5fl5XisWqfUtnz224LUlOKKSQ4zT0nk2RCN2c0NWGtvbrB6M3652PPtHzb67+Dtx925fX3A8A8A/CbJd59eSHun0XP0CBbrTwwotV42rW8sEjj6zad/W+ffrsrbcaes/Q4cOr17UGekhpIigAgH8RZtVdsxn//ANbfz71juYvLnrdFABdhhkXAQCAnxBm1V2zGf/zL3+tf/7lr71uBoAuxIyLAADATwizABAQA/fdq2137q/at3T2XFU34euPHFHfzuoJbdYaM1vrPZsiETebjh4U1F5PpfVk6+3hdOHiop5742x5m3VoAaBzCLMAEBCbBwdr3gnt27lD/bt3N3TcZt4D1BLUXk+l9WTrVby0pB9Pv1veZh1aAOgcluYBAAAAAAQOYRYAAAAAEDiEWQAAAABA4DBmVsGdpAIAAHTWwqdLXjcBgEsuz8/rSrFYtW/p7LkNt6XVEyu2+32bIhFWD1gHYVbBnaQCAAC0JrKlT4+OLE+ANn9xUc+/ufEEUK+e/FC37+JHJdANFl54sWrJubVUzv5fcvW+ffrsrbc69r6hw4f1xe+v3g/CbM+qXErg06XLkqRr+pb/OLCsAACgVUHp9TS0tV8/GL1ZkvTOR5/UDLMAAP8gzPaolUsJVGJZAQBAq+j1BABoN8IsAABAnf7l8hX9+dQ7kpZ7MgEIroH77tW2O/dX7Vs6e66qq+/1R46ob+eOqtesN/a1Xe/bFIk08F/VWwizAAAAdcrP/LL873fdep2HLQHQqs2DgzUnVurbuUP9u3c3fOxOv69XEWa72He+ukODW/rrmtACAAAAAIKEMNvF7v7dL+v2XYNMaAGgYc0uV8DyAegm3/nqDskRNRQAfIow28VYPgBAs5pdroDlA9BNvnvHjZIaC7OsQwsAnXOV1w3wg/HxcYVCIWZbBAAALXn15IdeNwEAegZhVsth1nEcnTp1yuumAADQFbhQDABoN7oZd7G7b/uy100AEFDNLlfA8gEoYZ1ZAEC7EWa72MA1fV43AUBAtXO5AqCbcSEZADqHbsYAAAAu4UIyAHQOd2YDLrKlT4+OLN8ZWW892crXlLaLl5htEegGmyIRDR0+XLXdjvcAAAD4DWE24Ia29usHozdL0rrryZZec+Hiop5746yee+Os5i8udrqpANpg8+Bgw0vhNPMeL1yen9fCCy+Wtwfuu5c1bAEAQBlhtocULy3px9Pvet0MAKjLlWKxaq3bbXfuJ8wCAIAywiwAAIBL/vqN9zS4pV+SdP8dOzW0td/bBgFAhW7r9USYBQAAcMnzP/t8uM9dt15HmAXgK93W64nZjAEAAAAAgUOYBQAAAAAEDmFW0vj4uEKhkPbs2eN1UwAAAAAAdSDMajnMOo6jU6dOed0UAAC6AheKAQDtRpgFAACu40IxAKDdCLMAAAAAgMBhaR4AAIB1RLb06dGR3ZKk+YuLev7NczXeASDINkUiGjp8uGrbj+/DMsIsAADAOoa29usHozdLkt756BPCLNDlNg8O6ovff8T378MyuhkDAAAAAAKHMAsAAAAACBzCLAAAAAAgcAizAAAAAIDAIcwCAAAAAAKHMAsAAAAACBzCLAAAAAAgcAizAAAAAIDA2ex1AwAAuDw/ryvFYtW+pbPnNtyWpE2RiDYPDra1bQAABFEv1FbCLADAcwsvvKgLTz+94Ws+eOSRVfuGDh/WF7+/ej8AAL2uF2orYVbS+Pi4nnjiCa+b0bLIlj49OrK7arsZf/3Ge7r632ySJF3Tt/xH5P47dmpoa78uXFzUc2+cLb+2tL+k8vlPly6veQwAwbQpEtHQ4cNV2/W4PD+vhRdeLG8P3HdvW6/4dvrzsLYg1tZG62ipXt6+K6Jbr79Wny1d0fNvnlv1/H/94Fe69fprdU3f5rrraaXSaz9durzmsQCg3fxaWwmzWi644+Pjevvtt7Vnzx6vm9O0oa39+sHozS0f5/mfre5ucNet12loa7+Kl5b04+l3V+0vWfn8WscAEEybBwebulJ7pVisujK87c79bS2Anf48rC2ItbXROlpZL//3A8v/jZVhtvL5N88sd/Wrt55WWvnalccCgHbza20lzAIAPDdw373aduf+qn1LZ89VdX+6/sgR9e3cUfWaeu8OAwDQa3qhthJmAQCe2zw4WPMKb9/OHerfvXvD1wAAgGW9UFtZmgcAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4Gz2ugHwzne+ukODW/r16dJlSdJnS1f0/JvnPG4VAAD+FNnSp0dHdkuS5i8uUjMBwGOE2R723Ttu1M1f+kJ5+52PPqEwAwCwjqGt/frB6M2SqJkA4Add0c04n89reHhYw8PDmpiY8Lo5AAB0BeorAMDPAn9n1rZtmaapubk5SdLw8LDi8bhisZjHLQMAILiorwAAv/PVndnZ2VlFo9E1n8tmsxoYGNDAwIDS6XR5f6FQqCqshw4d0rFjx9reVgAAgoL6CgDoRr4IsxMTEwqFQkomk7Isa9Xz+XxepmlqenpaMzMzOnHihJLJpCTJsiwZhlF+bTgcXvMYAAD0GuorAKCb+SLMjo2NyXEcZTKZNZ83TVOZTEaxWEyGYSiXyymfz8u2bUkq/1OSIpFIR9oMAIDfUV8BAN3M92NmbduWZVmKx+PlfYZhKBwOq1AoKBwOq1gslp9beSUZ7lj4dKnqn7VeB6A3XZ6f15WKv5MlaensuQ23JWlTJKLNg4Md+bxmP6vbUF9b02y9q7ee1vNZ1FwA7dDpWt4K34fZUpemlQXUMAxZlqVUKiXTNGXbtsLhsI4dO6ajR4/WPO7HH3+s8+fPV+07ffq0ew3vMq+e/FC37xrUqyc/rPk6AL1r4YUXdeHppzd8zQePPLJq39Dhw/ri91fvb8fnNftZ3aYd9bWXamuz9a7eelrPZ5WOBQBu6nQtb4Xvw2xxxVWBkkgkovn5eYXDYeVyOY2MjEhanqCinpkWn3nmGT3xxBOuthUAgKBoR32ltgIAOsn3YbYe8XhcMzMzDb3n4YcfLk9yUXL69GkdOHDAzaYBABBYjdZXaisAoJN8H2bXm3CiWCxqsIU+2du3b9f27dubfn+vufu2L5f/+fybq/vIV75uo+cBdLeB++7Vtjv3V+1bOnuuqjvS9UeOqG/njqrXbGpycqFmPq/Zz+o27aivvVRbm6139dbTej6rdCwAcFOna3krfB9mS2N5Vk48YVmWwuGwV83qOQPX9FX9s9brAPSmzYODNSd/6Nu5Q/27dwfy87oJ9bU1zda7eutpPZ9FzQXQDkGqrb5Ymmcj4XBYhmEon8+X91mWJdu2dfDgQQ9bBgBAcLW7vo6PjysUCmnPnj0tHwsAgLX4PsxKy+vgmaap2dlZWZalZDKpRCLh2pVjCi4AoBe1s76Oj4/LcRydOnXKhZYCALCaL8JsNptVKBQqTxoRCoUUCoXKz6dSKWUyGY2MjCgajZYXdncLBRcA0I28rq8AALSTL8JsKpWS4zirHpXGxsa0sLAgx3EotAAA1IH6CgDoZr4IswAAAAAANIIwK8bMAgDgNmorAKDdCLNizCwAAG6jtgIA2o0wCwAAAAAInM1eN8BPFhcXJUmnT5/2uCXt8d6FS1o6f7a8/e4v/l/964Ut6z5fqfRaN44BoLcsnT2r93/796sk/cu776pvacl372tEqU4sVnwO1tattXWteihp3RpY+bp66ulGn1XPewBgLX6urVLj9TXkrJzWsIf97d/+rQ4cOOB1MwAAAfHKK6/om9/8ptfN8DVqKwCgUfXWV8JsBdu29frrr+uGG25Qf39/zdefPn1aBw4c0CuvvKKbbrqpAy2EGzhvwcM5C6ZuPm+Li4t6//339bWvfU3hcNjr5vhao7V1Pd3856mbcJ6CgfMUDL14nhqtr3QzrhAOh5u6wn7TTTfplltuaUOL0E6ct+DhnAVTt563WCzmdRMCodnaup5u/fPUbThPwcB5CoZeO0+N1FcmgAIAAAAABA5hFgAAAAAQOIRZAAAAAEDgbBofHx/3uhFBtmXLFv3e7/2etmxhavwg4bwFD+csmDhvcBN/noKB8xQMnKdg4DxtjNmMAQAAAACBQzdjAAAAAEDgEGYBAAAAAIFDmAUAAAAABA5hFgAAAAAQOITZJmWzWQ0MDGhgYEDpdNrr5qAOyWRSoVCo6hGNRr1uFlaYnZ1d97zwvfOv9c4b3zu4ge++/9i2rXQ6veF54bz5i2VZymazsm27aj/nyT9M01Q0GlUoFNLo6GjVc5yntRFmm5DP52WapqanpzUzM6MTJ04omUx63SzUIZVKaW5urvyYmpryukn4rYmJCYVCISWTSVmWtep5vnf+VOu8SXzv0Bq++/5UOgczMzM6evSojh8/XvXjm/PmP8lkUul0WsVisbyP8+Qfo6OjKhQKmpqa0sLCgkzTLD/HedqAg4YZhuFMTk6Wt+fm5hxJzsLCgoetQi2JRMIZGxvzuhmoIZfLOWv91cT3zt/WO29879Aqvvv+Mzc354TD4ap9MzMzjiRnZmbGcRzOm99MTk46iUTCkeTMzc2V93Oe/GFycnLVd6oS52l93JltkG3bsixL8Xi8vM8wDIXDYRUKBQ9bBnQvvndAb+K770+GYejo0aNV+2KxmCTpxIkTnDcfMk1Tjz/+eNU+zpN/ZDKZVeenhPO0McJsg0rd6AzDqNpvGMa6XezgH6VxfQMDA3TPCBC+d8HG9w7N4rvvX4lEYs39kUiE8+Yzpmlq79695QsOJZwnfyiFVcMwlEwmNTAwoGg0qtnZWUmcp1oIsw2qHGdQKRKJaH5+vsOtQaMKhYImJyc1MzMjy7JWDa6HP/G9Cza+d2gW3/3gyGazCofDSiQSnDcfsW1bExMTymQyq57jPPlDKZCapql0Oq2ZmRnFYjGNjIxI4jzVQphFzzh06JCmpqYUj8dlGIYymYwKhQJXtYA24nsHdD/LsmSapnK5nNdNwQoPPPCAEonEqruy8I9SWM1kMuVamcvlZNu28vm8x63zP8JsgyKRyJr7i8WiBgcHO9waNCKRSFSNN9i7d68kMd4gAPjeBRffO7SC777/2bat0dFRHT16tPxd57z5w+zsrPL5/KrxzSWcJ38odR9eecHBMAy99dZbnKcaCLMNKv2BW3lXwbIshcNhL5qEJpWuhK33lwT8g+9d9+B7h0bw3fc327Y1PDysycnJqjG0nDd/OHbsmCRp165d5fVJJSkajWp4eJjz5BOl87By/V9JGhwc5DzVQJhtUDgclmEYVbf9LcuSbds6ePCghy1DLSu7apQG1tP1xv/43gUX3zu0gu++f9m2rZGREU1OTlb1vpA4b36RyWS0sLCgM2fO6MyZM5qenpYkTU1NaXp6mvPkI7FYbFWPJcuyFIvFOE+1eL02UBBNTk6W11Kbm5tzYrGYk0gkvG4WNlBaE29yctJZWFhwZmZmnHA47KRSKa+bhhXWW6+U752/rXXe+N7BDXz3/SkWizljY2POzMxM1aO07iXnzX9Ka5NWrjPLefKHqakpR5IzNTXlzM3NOYlEwonH4+XnOU/rI8w2KZPJOOFw2JHEH6aAmJubc+LxuCPJMQzDyWQyXjcJFUp/Ua98VOJ75z+1zhvfO7iB776/lELRWo/K7zjnzV/WCrOOw3nyi1wu5xiG4YTD4TXPA+dpbSHHcZw23/wFAAAAAMBVjJkFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFekAymVQ6nfa6GQAAAIBrCLMAmmLbtvL5vKLR6JrPT0xMKBqNamBggCANAAAA1xFmATTMNE3t2rVLTz75pCzLWvV8NpvV5OSkcrmczpw5I8uyNDo66kFLAQAA0K0IswAalslktLCwoMcff3zd5zOZjGKxmMLhsHK5nAqFgmzb7nBLAQDwh9nZWQ0MDHjdDKCrEGaBgEmn05qYmJBpmuVuvNlstnz3MxQKaXh4uGZwrOwGnEwmV70+mUxqYGBA0WhU6XS6/NpaLMuSZVmKx+PlfeFwWIZh6Pjx45KW7+yWHqXjTkxMNPF/AwAAAL2KMAsETLFYLIfAmZkZpVIppdNpJZNJTU5Oam5uTrZtyzTNdY+RTqd17NgxTU1N6cyZM4pEIhoeHq56XpLOnDmjTCajbDarmZkZLSws1GxfqdtxOByu2m8Yhubm5sqvmZiY0L59+zQ1NaVUKiXTNJXP5xv+/wEAgN+Njo6WLzSHQiGFQiEVCgVJKtf0UCikaDRKLQQaQJgFAsgwDKVSKYXD4XLwPHTokAzDkGEYSiQSOnHixJrvtSxL2WxW09PTMgxD4XBYk5OTkpbHukrS8ePHlU6nFQ6HlUgkFA6H1xwb24pEIqFEIiHDMJTJZBSPx8vtAACgm0xNTWlqakrhcFiO48hxHMXjcc3OziqbzWpqakoLCwvUQaBBhFkggGKxWPnfDcNYtW9wcHDd9xYKhXKIrRSPxzU1NeVySz9XLBY3fD6ZTK4bwAEA6EaWZSkSiZTrcjweVyKR8LpZQGAQZoEAikQide1bSz2TMB08eFCZTEa2bWtiYkKRSKQqLG+kFK5Xfo5t2+su41NS738DAADdoDS/RCgU0ujoaLmHFID6EGaBHhOPx8uTNFUqFArat2+fpOUrxeFwWLt27SqPra1X6epyabInaTnIrpwUaqVcLld3YAYAoBuEw2HNzc0pl8vJMAyl02nWZgcaQJgFekwsFlMikVAymdTs7Kwsy1I6nVaxWNTY2Jik5S7Bhw4d0szMjHK5XMN3TB9//HGZpqnZ2VnZtq1kMlkeH1uSz+fLy/WYpqlCobDuUj8AAHSD9XpHJRKJ8vrslReDAWxss9cNANB5uVxOpmkqmUyqWCwqHo/rzJkz5ecNw1Aymax6TzweVy6XUzgcVjabrbpyHAqFJEmO40hSORSXjn/w4MFVk1qUjpdMJhWJRLgzCwDoaqULuqVZjCORSNVyduFwWFNTU1UXfgFsLOSUfn0CgJaLbDKZ1MzMTLmgzs7O6oEHHtDevXtdmWmxFJRzuVzLxwIAICiGh4dlWZb27t2rTCYjSXryySfLPZVisVi5yzGA2rgzC6DK7Oys9u7dW1VIY7GYDh061NbZjgEA6HYzMzOr9nFhF2geY2YBVInH4yoUCspms+WxPYVCQZOTk6u6HgMAAABeIcwCqBKLxcoTP+3atUuhUEimaco0TaVSKa+bBwAAAEhizCwAAAAAIIC4MwsAAAAACBzCLAAAAAAgcAizAAAAAIDAIcwCAAAAAAKHMAsAAAAACBzCLAAAAAAgcAizAAAAAIDAIcwCAAAAAAKHMAsAAAAACBzCLAAAAAAgcAizAAAAAIDAIcwCAAAAAAKHMAsAAAAACBzCLAAAAAAgcAizAAAAAIDA2ex1A/zEtm29/vrruuGGG9Tf3+91cwAAPrW4uKj3339fX/va1xQOh71ujq9RWwEA9Wq0vhJmK7z++us6cOCA180AAATEK6+8om9+85teN8PXqK0AgEbVW18JsxVuuOEGScv/82666SaPWwMA8KvTp0/rwIED5bqB9VFbAQD1arS+EmYrlLo/3XTTTbrllls8bg0AwO/oNlsbtRUA0Kh66ysTQAEAANeNj48rFAppz549XjcFANClCLMAAMB14+PjchxHp06d8ropAIAuRZgFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBQ5gFAAAAAAQOYRYAAAAAEDiEWQAAAABA4BBmAQAAAACBs9nrBgDwl8vz81p44cXy9sB992rz4KDr7wEAAEAwXLi4qOfeOFvevv+OnRra2u9hi5Z1RZg1TVOFQkG2bSudTmtsbMzrJgGBdaVY1IWnny5vb7tzf81g2sx7AAAAEAzFS0v68fS75e27br2OMOuGfD4v27Y1MzMjSYpGo4rFYorH4x63DAAAAI2itw/grQsXF1W8tFS1770LlzbclqTIlr6OB1xfhdnZ2Vklk0nNzc2tei6bzco0TUnSwYMHNTk5KUmKx+NVwTUej2t2dpYwC9Th8vy8rhSLVfuWzp7bcNsJhRRynIbeI0mbIhF+jAAAanKztw/BGGjcc2+crboLu5bUczOr9j06sls/GL25Xc1aky/C7MTEhEzTlGEYsixr1fP5fF6maWp6elrhcFjJZFLJZFK5XE7hcLjqtYVCQel0ulNNB9qunYV44YUXq34wrOWDRx6p2r563z599tZbDb1HkoYOH9YXv796PwAAlRdXN7pAWnlhtJ76yDAYoLv5IsyOjY1pbGxM+XxeyWRy1fOmaSqTySgWi0mScrmcotGobNuuCrOmaSqRSJRfB3QDCjEAoNttdHG18gJp5YXR9epjM8EYQDD5IsxuxLZtWZZV1W3YMAyFw2EVCgUlEglJy0E2Go0qlUp51VQAAACsw42eRr/59FMtvrvc/XG9oGq//LKK/8dfrfn+9YIxgM/df8dO3XXrdVX73rtwqaprcfb+Yd04tKXqNZEtfR1pXyXfh9lSt2PDMKr2V3ZJTqfTSiaTjJMFGjRw373aduf+qn1LZ89VFfvrjxxR384d5e31xsxu9B5p+Qo4AKB3udHT6LNTp1T8q42D6tX79jXfSAAa2tpfcyKnG4e26OYvfaFDLVqf78NsccXkNCWRSETz8/PK5/PKZrMqFArl5+pZnufjjz/W+fPnq/adPn269QYDG6jnqvTKSZmanVypns/aPDhY84dE384d6t+9e8PXuPEeAED3qbfL7xd+/87yxdWNLpDaL79cc94GAO6LbOnToyO7q7b9wPdhtpZEIiFnxV2iejzzzDN64okn2tAiYH31XJWuNSlTvZMrMdYWAOC1ZsbCrlR5gfSqq6+p+ZlX79mj6/63/1XSxsGYHkNA/Ya29nd8puJ6+D7MRtb5i6ZYLGqwhR/mDz/88KrJpk6fPq0DBw40fUygV22KRDR0+HDVNgAAbqscHrNeUN1oYid6DgHdxfdhtjRW1rKsqnGzlmWtWpanEdu3b9f27dtbbh/gB5UTYpQ02z25mWC6eXCQSTQAAK7YqA5tNDxmraDKxVagu/k+zIbDYRmGoXw+Xx4Ha1mWbNvWwYMHPW6duy5cXNRzb5wtb99/x86ag6/hXyvHvkr1BczKcUOl19SaXMl++WVZd//Bhu2pt3sywRQA4JZ67qRKqwOrW3WImgZ0N9+HWWl52Z10Oq14PK5wOKxkMqlEItHSnVk/Kl5a0o+nP7+7dtet1xFmA6zW2Fep/oBZaa0rz/WMIQIAoNMavZPaCO66AvBFmM1ms0qn0+XtUCgkSeWJnVKplGzb1sjIiGzbViKRUC6Xc+3zx8fHmQwKAAAgQLjrCsAXYTaVSimVSm34mrGxsZrL7TRrfHxc4+Pjevvtt7Vnz562fAbQrHquPDezXux6xwIAoB24kwrAbb4Is73swsVFFS8tSZLeu3Cp6rnK7ciWProcB4xbAbOeK8/tWi+23epZCxcA0B24kwrAbYRZjz33xtmqcbKVUs/NlP/90ZHdvlzbCesLasDsJNbCBQAAQLOu8roBAADAn0zT1PDwsKLRqCYmJrxuDgAAVbgzC3QpxiYBaEU+n5dt25qZWe4lFI1GFYvFFI/HPW4ZAADLCLMeu/+Onbrr1uskLY+RrexanL1/WDcObZG0PGa2hPVog6uTAdNvY5OaXXd3UyRC12OgRbOzs0omk5qbm1v1XDablWmakqSDBw9qcnJSkhSPx6uCazwe1+zsLGEWqzD/AQCvEGbl7dI8Q1v71w2iNw5t0c1f+sKq/axHG1x+C5id1K51dwGsb2JiQqZpyjAMWZa16vl8Pi/TNDU9PV1exz2ZTCqXy61ay71QKFQtoweUMP8BAK8wZlbLYdZxHJ06dcrrpgAA4JqxsTE5jqNMJrPm86ZpKpPJKBaLyTAM5XK5cvfila9LJBKKxWKdaDYAAHXhziwAAD3Itm1ZllXVbdgwDIXDYRUKBSUSCUnLQTYajdZcDx69pXLoyEZDRhgqAqCdCLM+EtnSp0dHdldtl7AeLYLOrXV3GZsFuKPU7dgwjKr9lV2S0+m0kslk3eNkP/74Y50/f75q3+nTp11oLfxmo6EjlX+vM1QEQDsRZn1kaGv/umvJsh4tgs6tdXcZmwW4o7hiQraSSCSi+fl55fN5ZbNZFQqF8nPpdFpjY2PrHvOZZ57xbA4KAEDvIcx2KWY8bi/uDgLodolEQo7jNPSehx9+WMlksmrf6dOndeDAATebBgCAJMJs12LG4/bi7iCAoIusszRYsVjUYJN/n23fvl3bt29vpVkIiMqhIxsNGWGNcwDtRJiVt0vz1KuZ9WgBv6tn3d2V69OyNi3gjtJYWcuyqsbNWpa1alkeYKWNho7UM2QEANxAmNVymB0fH9fbb7+tPXv2eN2cNTWzHi3gyBbwnQAAIABJREFUd/Wsu1trfVrWpgWaEw6HZRiG8vl8eRysZVmybVsHDx70uHUAANRGmO0Axq8G28o7gxJ3BwF0B9M0lU6nFY/HFQ6HlUwmlUgkuDMLAAgEwmwHdHL8amkJH5bvcU+tO4MSdwcB+FM2m1U6nS5vh0IhSSpP7JRKpWTbtkZGRmTbthKJhHK5nCufHYQhPHBHPUNGAKAdCLNdZr0lfFi+B0G1cn3aZtamBXpVKpVSKpXa8DVjY2MbLrfTrCAM4YE76hky4nesUgAEE2E2gCJb+vToyO6qbaBb1VqflolGAACtYpUCIJgIs21S6u4ryfUuv0Nb+7mz2kEr7wxK3B0EAAAAvEaYbZP1uvtK7e3yW1rCh+V73FPrzqDE3UEAAIKmcoLHjSZ2ZEJHwL8Is+quSSrWW8KH5XsAAAA+t9EEj5W9r5jQEfCvq7xugB+Mj4/LcRydOnXK66YAqKE0a2bpQXduwJ/Gx8cVCoWY/AkA0DbcmW2TUndfSXT5BVzUDbNmAr2A2YwBAO1GmG2T9br7Sv7q8nvh4qKee+Nsefv+O3ayBm0dWFMveFh2AQBQqXKCx40mdqTGA/5FmO1S9S7fU7y0VDVR1V23XkeYrQN3B4OHZRcAAJU2muCRiR2BYCDMdim3l+/hDi4AAAAAPyHMoi7cwUWvonsyAACAPxFmAWADdE8GmtNNy96h+zEXBhBMhNkOqHf8aqdcuLio4qUlScszLVeq3I5s6ePuKwLp8vy8rhSLVfuWzp7bcFta/vFCUAXcwWzGCBLmwgCCiTDbAW6PX23Vc2+creoyXKlyCaFHR3b7qt210B0UJQsvvFh1N3UtlbNWlgwdPsyPGQAAgIAgzIquUOv5dOmy3vnoE0nBuINLd1AAANBOXDgH/IUwK7pCree/fvArHf0vZ9Z8Lsh3cIH1/ObTT7X4bnWvhVrdk+maDAC9gwvngL8QZnvQ/Xfs1F23Xidp+Q5rZTDN3j+sG4e2SJJyJ97Xm2eKax4D8LOB++7Vtjv3V+1bOnuuqmvx9UeOqG/njqrX2C+/LOvuP9jw2Cu7J9M1GQAAwBuE2R40tLV/3W7BNw5t0c1f+oIk6Zq+7v/jQXeh7rR5cLDmeezbuUP9u3dX7bvq6mva2SwAAAC4qPvTCppW7x1cr2dnbgXdhQAAAIBgIsxiXfXewfUCS6+gHZrpnsxahMDamFwR3aLyN8dGvzX4jQF0HmEWgcTSK2iHZrsnA1iNyRXRLTb6zVH5W4PfGEDnXeV1AwAAAAD8/+3df2xb553n+w9tjbQd2Q5JKfE4N2PZR5W3iwncllIw6wssZjCi6gBZYNKGtGFf+J9tQyZVAqMXHTGe/Yfei52EygUWxTrNkO5/3k3GIlt0gP2rpAYYXOB6B5E02YsYu0isIyv1zk4SiTpp7GakWuH9QyVDShRFkYc/Dvl+AYJ9Dn89h+fHw+95nuf7ANgvWma7nLe/V5fHR0qWgU500OvV4ORkyXKjkFgMQCfgWgag3RHMdrnBQ32OnCO21qlXci5XyTyijLPtHj0DA03r/kViMQCdgGvZluLfHJV+a5BDAWg+glk4Uq1jGz/5j9cqjrVlnC0AAChW6TcHeRSA1iKYRVXs7o68cn9dN24tF5YvnRnaNXMyAAAAAGxHMCumD6iG3d2Rsw829KPZL7v7PnP6GMEs2lIzx9oCQKsxDQ0AJyGYFdMHdJPtY22rGWdL8NLdmjnWFgBajWloADgJwSy6yl5jbRn7gv0qbsXII7EYQK8ndCZ66wDthWAWAOpQqRUjj8Ri6Eb0ekInorcO0F4IZtExuFsKAEB9mIYGgJMQzKJpVu6vK/tgQ5J0d+VByWPFy97+3pqSQXG3FACA+jANDQAnIZhF09y4tVySwbhY6MZ84f+Xx0dszZwMNNL2pGISicUAAACagWAWAOqwV1IxidYMAACARiCYRVdjnC0AAADgTASzaJpLZ4b0zOljkrbGyBZ3LU5cGtWJwX5JW2Nmm4VxtgAAAIAzEcyiaQYP9e2a2OnEYL9OHT3c5BIBAIDd0HsJQLsjmAUAAMAO9F6q3cPVVa299XZh2XPxwp75FQDsH8EsAAAAYKPNbFYrb7xRWD7y9FmCWaABDrS6AAAAoPNEo1G5XC49+eSTrS4KAKBDEcwCgM3y48zyf4wzQzeKRqPK5XJ67733Wl0UAECHopsxANiMcWYA4Bwr99d149ZyYfnSmaFdE1YCaC8Es9q6e3z16tVWFwMAAABNln2woR/NflBYfub0MYJZwCEIZrUVzEajUd2+fZuxPU3i7e/V5fGRkmUAAACneri6qs1sVpK0sfxhyWPFywe9XpJBATYhmEVLDB7q0w8mTrW6GAAAALZYe+vtkgzGxe699OXQk8HJSYaiADYhARQAAAAAwHEIZgEAAAAAjkM3YwAAAHSFlfvryj7YKFl3d+VBxWVpK7fHXkmhPBcv6MjTZyVtjZEt7lr8xLVr6h06LklM1wbYiGAWAAAAXeHGreWSzMXlhG7M71h3eXxkz1wfPQMDuyZ26h06rr6RkbKPAagd3YwBAAAAAI5DMAsAAAAAcBy6GQMAAKArXDozpGdOHytZd3flQUnX4sSlUZ0Y7C95jre/tynlA7A/BLMAAABd5OHqqtbeeruw7Ll4Ydexnp1m8FDfnomcTgz269TRw00qEey0cn9dN24tF5YvnRnac3/D2WwLZt99992y67/xjW/Y9RHoQlyUAHQz6lY0wmY2q5U33igsH3n6bNcEs+hs2QcbJQm+njl9jN+NHa7uYPb69et64YUXJEm5XK7kMZfLpc3NzXo/Al2MixKAbtQJdWs0GtXVq1dbXQw4TKfcxD7o9WpwcrJkGYD96g5mI5GInn/+eUUiEXk5UQEAqFsn1K3RaFTRaFS3b9/Wk08+2eriwCE65SZ2z8CAHn35pb2fCKAutnQzfuWVV3TixAk73goAAIi6FdgPp7bodvP4ZcAOdU/NEwqF9NOf/tSOsgAAAFG3AvuVb9HN/2UfbLS6SFXJj1/O/21ms60uEuAotrTM/tmf/Zn+4i/+QmNjY3K73SWP3bx5046PAACgq1C3AgDKcWpPhEaoO5g1TVOBQKCwvD1RBQAA2B/qVtjt4epqodVvY/nDkseKlw96vV3XzdXb36vL4yMly2h/K/fXd7TA3115UHFZ2tq/Tg/8OmVsuR3qDmZnZmbsKAcAAPgt6lbYbe2tt0um4yl276UvExUNTk52XeKiwUN9+sHEqVYXA/t049ZySUBXTujG/I51l8dH2N8dxLZ5ZgHJnm4PH334j1pK/bUk6ZcPeyT9XuGx4jtsnXBnDQCAbtXNLWsA7LHvYPbs2bMKBoP63ve+J0l68cUXKz7/zTffrK1kcCQ7uj3cuLWkax8/Ufax4jts3FkD0CmoW9GNaFlrLsZZNgbfa2vtO5hdXV2tuOxETOwOAGilTqxb0V48Fy/oyNNnJW2NkS3uWvzEtWvqHTouaWvMbLv79cZDvf/RZyXrqmnRldqjVbdV45c7bZzlpTNDeub0sZJ1d1celNwASVwa1YnB/pLn2D0mutO+V6fZdzA7NzdXstwJ43qY2B0A0EqdWLeivfQMDOwaGPUOHVffyEjZx9rR/3fvU13/f5YqPqdci65kb6turS1y9YxfphXwS4OH+vbc9hOD/Tp19HCTSoRWYMws2s7Ff+HWk//3FUnSP/QP6N/9y39TeKz4DhvZBgEAcK5aW9aSc7/U3y21fj7WVrTI0QrYfRhbXpktwezPfvYzxePxwp3lsbExTU9P6+tf/7odb48uM/CVHg199lHZx7jDBqBbULei09Xasva7vbTFoHswtryyuq8Gr7/+uiKRiAKBgF577TVJW92lfD6fUqmUvv3tb9ddSAAAugl1K7C7Wlt0pfbo1dVJ45eBVqs7mH311Vc1PT2tH/7wh4V1zz//vCYmJjQ1NUWFCwDAPlG3Artz+ljJThq/DLSaLf00AoFA2XWvvPKKHW8PAEDXoW6Fk5GoCKis2nOkXbI2t6u6g9krV64ok8kU5sbLu3v3rkZHR+t9ezhA8cD0SgPSyw1EL05Pn7c9TX3pY8ta/1Wf7enqgXbwcHVVa2+9XVj2XLzAcd6lqFvhdCQqap1yCYMkkgbZwc5kTNWeI07vidBo+w5mt0/kvrq6qp/+9KdKp9Ml6zOZjAzDqK90cIRKA9OL7xqVG4heKT19OfcmX9LBzz4qm64ecLrNbLbkfDjy9FmC2S5B3QrALtUkDJK6O2lQrUjG1H72HcyWm8j9ueeeUy6XK1k3Pj5ee6kAAOgi1K1opoNerwYnJ0uWUZnd06Pku5huZle19rVv6V+b/6/cG/dt+cxfbzysvDEdzNvfq8vjIyXL6Gz7DmaZyB0AAHtRt6KZegYG6N20T3a3yJV0Mf3at/Sv/ud/2xHM1vqZf3iye29ODB7qowW0yzBRF+pWPDC90oD0cnfHitPT521PU1/siTeuyfD07biLTKIJAGgv0WhUV69ebXUx4DC0rNXv9BOP6P969skd60kaVD+nJ2PqxN/LBLOoW6WB6XsNSK+Unr6c3qEh9ZV5PxJNAEB7iUajikajun37tp58cucPa6AcWtbq97u9PVUlA+rmpEG1cnoypk78vUwwCwAAAFTQihY5p7cCAs1AMIu29sj6A730Ta8Oerdab7lAAwCA7RrdPbnRLXLuCxc0+Du/KRlG5fRWwHp0YndYNAbBLNqae+O+XvYNqm9kZO8nAwCApmqXOU2d3j3Z+39c1KMdGJTWqtO6w9aamVpi7t+9EMwCAACgJsxpCuyt1szUUvnzhERpXzpQ7xv85Cc/0a9+9Ss7ygIAAETdCgDYXb4nQv6vm1tu626Zfe211+RyufTd737XjvIAQFd4uLqqzWy2ZN3G8ocVlyXpoNe7rwzgcCbqVqAztEs3bKBT2RLMhsNh+f1+DQ0N2VEmoGGqTShA4gE02tpbb2vljTcqPqfcfMuDk5N69OXS9Q9XV7X21tuFZc/FCwS8DkfdCqcol3FXIutuHt2wIdWemVrqjvOkHnUHs3Nzczp58qQMw5Df75dhGCWPv/nmm/V+BGCbahMKdFriAXS2zWy2JDA+8vRZglmHo26FU1STcVfq3Ky7QDW6OTN1o9UdzJqmKcMwChXt6upq3YUCAKCbUbcCsFutSYNq7a3WLb3c2jUZU7dkUK47mJ2ZmbGjHMC+FZ+klU5OJ52Q6B6eixd05OmzJes2lj8s6Vr8xLVr6h06XvKc4jkI0bmoW4HO0E7dsGudvqjW3mrd0sutXaeFsjuDcruyZWqen/3sZ/qrv/orvfvuu3r//fclSS+++KKCwaD+5E/+xI6PAHaodJIWn5xOOiHRPXoGBvbsCtw7dLxhcywzzrb9UbcCzkc3bKCx6g5mX3/9dSUSCf3lX/6lvvWtbxXWj4+PKxaLUeF2GTu6Whz0ejU4OVmyDMBejLNtb9StQPtrRffSdu3SWo9au8PS8w6STdmMZ2dn9Y1vfKNk/ejoqM6fP1/v28Nh7Ohq0TMwsCNbKwB0E+pWoP21ontpu3ZprUet3WHpeVdZt2RQrjuYzeVy8ng8O9abpqmTJ0/W+/bAropP0konp5NOSACQqFsBAPXplgzKB+p9g+eff16hUEi/+tWvCuuWlpb0wgsv6IUXXqj37YFdDR7q06mjh3Xq6OEdd5XyJ+epo4fpggLAcahbAQDYW90ts7FYTOFwWG63W5I0MjIi0zQVCoX0wx/+sO4CAgDQbahbUQsSu3WuavetE4+BWrvD0vMOkk3ZjOPxuCKRiP7+7/9ekuTz+egGBQBAHahbsV8kdutc1e5bJx4D3dIdFo1hSzArqWRy91YwTVOxWEzxeLxlZQAAwE6trluBWnVi1l3ATpwj9rBtntl4PK65uTlJ0tjYmKanp/X1r3/djrff08TEhLLZrLxM4YIyilO+V5Pqfa/nkQoerfRwdVWb2WzJuo3lDysuH/R62/7OPHZqdd0K1KMTs+52g1qnySmn2td16+8qzhF71B3MvvLKK3r99dc1Pj6u1157TZI0MzMjn8+nVCqlb3/721W/18LCgoLBoBYXF3c8lkgkFIlEJEnnzp0raYFNp9NKpVK0yqKsSinfi8dj/OFJr/5uKbvn80gFj1Zae+vtki5k5dx7qXRqq8HJSaa7chg761agViv313Xj1nJh+dKZoa4MOrpJrdPkVPoNVel1Uvv8ruJ4d6a6g9lEIqGpqSm9+uqrhXXPP/+8IpGIpqamqqpwp6enFYlEZBiGTNPc8XgqlVIkEtHs7KzcbreCwaCCwaCSyWS9xQcAoO3YUbcC9co+2CgJbJ45fYwf9+hYHO/OVPfUPF6vV+FweMf6P//zP9fq6mpV7zE1NaVcLqdYLFb28UgkolgsJp/PJ8MwlEwmlUqlZFlWXWUH6rFyf13/If1+4W/l/nqriwSgQ9hRtwJoX/nxkvk/xksCtam7ZTYQCGh2dlbf/e53S9bPzs7q/Pnz9b69LMuSaZry+/2FdYZhyO12K5PJKBAI1P0Z6GzFKd+rSfW+1/PyFQ538NAKnosXdOTpsyXrNpY/LOla/MS1a+odOl5YPkg+AcdpdN0KoLXadbxkrdPklFPt6wjkUY+6g1nTNPX666/rF7/4RWGdZVnKZDLy+/0lle7Nmzdren9JO7I57tYlGdiuUsr3alO9kxIe7aJnYGDPZE69Q8fVN/JlhsSHq6ta/6B0DNReSaMkEke1UqPrVgAop5HT5PBbCo1gSzbj5557TrlcrrD8yCOP6LnnnpOkkvW1yGbLDyb3er2FrlbhcFhzc3MyTVMTExOFLsmVfPzxx/rkk09K1t25c6eusgJAPQ56vRqcnCxZtkMtSaMkEke1WiPrVgBA9+nE6YDqDmZnZmbsKEddasli/OMf/1hXr15tQGkAoDY9AwMEj5DUHnUrAKCztGv39nrY0jLbSLvNHZvNZjVQR/e373//+woGgyXr7ty5o2effbbm9wQAoNOYpqlYLMb0d22qeO7pSsMHGDbgPNXu25zLJddve2vs9jz2PzpV2wez+bGypmmWjJs1TVNut7vm933sscf02GOP1V0+AGh3tSSNkkgcBWliYkLZbHbXG8tovUrDCIrPcYYNOE+1+/YrTz2lz995p+Lz2P/oVG0fzLrdbhmGoVQqpampKUlbgaxlWTp37lyLS4d20YljAAC71JI0Cp1lYWFBwWBQi4uLOx5LJBKKRCKSpHPnzpW0wKbTaaVSKVplATQVv+tQrbYPZqWteWbD4bD8fr/cbreCwaACgUBdLbPoLJ04BgAA6jU9Pa1IJLLrDACpVEqRSESzs7OF+jUYDCqZTLagtACwhd91qFZbBLOJRKJkcniXyyXpy2yNoVBIlmVpfHxclmUpEAhQ0QIAsIepqSlNTU0plUrtyBMhbd0sLp4BIJlManh4WJZlccPYIYqHEVQaPsCwAeepdt9uHzNb7nnsf3SquoPZn/3sZ5Kk73znO5KkF198UTMzMzIMQ8lkUidOnNjzPUKhkEKhUMXn5CvkRohGo2Q2BgC0DTvq1r1YliXTNOX3+wvrDMOQ2+1WJpNRIBCo+zPQeJWGEexn+MDK/XVlH2yUrLu78qDisrTV/XOveUlRGzv2LUNIdip3rEsc705VdzAbiUQKY2muX7+uRCKhmZkZ/eIXv1AwGNQ7uwxIbyfRaFTRaFS3b9/Wk08+2eriAAC6XDPq1ny34+Lkivnlcl2Sq8Ec7s5149ayfjT7QcXnhG7M71h3eXyE7qBwlGqOdYnj3SnqDmYXFxc1NjYmaat7UigU0nPPPadvfvObGuFOEAAA+9aMujX72yk/tvN6vVpdXZUkhcNhzc3NyTRNTUxMlHRJLoc53AEAzVR3MGsYhtbW1nTkyBFlMhm98MILkqRPP/2U8TbQw9VVrb31dmHZc/EC85wBDsY53RztUrfuN4sxc7gDAJqp7mA2f6d4YGBAhmEUxvdkMhmNj4/XXUA422Y2WzJH2pGnzzr+h29+rEWlsRWMqUCn6sRzuh01o27dbe7YbDargRr3KXO4N06jbyRdOjOkZ04fK1l3d+VBSVfLxKVRnRjsL3kOU6bAacod6xLHu1PVHczGYjGdP39eS0tLO5JIFC8DnWK3sRbFF0DGVACoRzPq1vxYWdM0S8bNmqZJz6o21OgbSYOH+va8CXtisF+njh627TPhDLXO+dquc8VWc6xLHO9OYcvUPD6fb8cYmueee86Ot24KshkDANpNo+tWt9stwzCUSqUKswWYpinLsnTu3DnbPgeAs9U65ytzxaIZ9h3Mvvjii/t6/ptvvrnfj2g6shkDAFqpVXVrJBJROByW3++X2+1WMBhUIBCgZRYA4Aj7DmbzGQ4BJ6q2y0ul5+XHWlQaW9EuXWnQHQ56vRqcnCxZhrM0qm5NJBIKh8OFZZfLJUnK5XKStuZ5tyxL4+PjsixLgUBAyWTSls+m1xMAoNH2HczOzMw0ohxAU1Tb5aXS83Yba8HYCrRKz8CAHn35pVYXA3VoVN0aCoUUCoUqPmdqaqrQzdhO9HoCADTagVYXAAAAAACA/bIlAdS7776rmzdvyjTNHY/dvHnTjo8AHGvl/rpu3FouLF86M8S0PQD2RN0KoBGYLxydpO6W2evXr8vn8ymdTiuXyymVSimXyymdTsuyLDvKCDha9sGGfjT7QeEv+2Cj1UUC0OaoWwE0Sn6ap/zfZjbb6iIBNas7mH3llVeUTqc1NzenmZkZPfLII5qZmdHMzIxjsiFGo1G5XC7G9AAA2kIn1K0AADRa3cHs2tqannrqqcKyYRj67LPP9NRTTymTydT79k0RjUaVy+X03nvvtbooANAU+QzI+T8yILeXTqlbuVEMAGikuoNZn89XUrGOjY1pZmZGs7OzdIUCgDaVz4Cc/2O8VHvphLqVG8UAgEarOwFULBZTOp3Wd77zHUlbE7APDw/L5XLtOR0AAADYiboVtWDO6c5V7b7lGEC3qTuYHR8f1/j4eGH55MmTWltbUzab1cmTJ+t9ewAAug51K2rBnNOdq9p9yzGAbmPL1DzbPfLII3rkkUca8dYAAHQl6tbu83B1tZBpdmP5w5LH8ssHvV6GCQDoWvsOZs+ePatgMKjvfe97kqQXX3yx4vPffPPN2koGAECXoG7tXpXm/Fx7622tvPFG2dfde2mr9W1wcpKWOABda9/B7OrqasVldK/iO8h5u91JzuOOMtC+ajmnJc7rWnRi3RqNRnX16tVWF6Pt5ef8zDvy9FnOHwCo0r6D2bm5uZLlmZkZ2woDZ6t0Bzkvfyc5jzvKQPuq5ZyWOK9r0Yl1azQaVTQa1e3bt5mex6G8/b26PD5Ssgx0Ko53Z6p7zOxPfvITnTt3TkeOHLGjPC3B3eO9VeoGheZaub+uG7eWC8uXzgxp8FBfC0sEwG6dULeifp6LF3Tk6bOStnpBFN88euLaNfUOHW9ottrBQ336wcSphr0/0E443p2p7mD2tddek8vl0ne/+107ytMS3D3eG92g2kf2wYZ+NPtBYfmZ08cIZoEO0wl1K+rXMzCwa13bO3RcfSMjZR8DgG5hSzAbDofl9/s1NDRkR5ngUMV3kPN2u5Ocx/xnQPuq5ZyWOK/tQN0KAMDe6g5m5+bmdPLkSRmGIb/fL8MwSh4n42L3qHQHOY87yYBzcE63DnUrAAB7qzuYNU1ThmEUKtpOyMAIAEArUbcCsFM1cxZLZKOH89QdzHZCxkUAANoJdSsAO1UzZ7FENno4z4F63+Bv/uZvyq5/99139e6779b79gAAdJ1OqFuj0ahcLheJFW2y+vlD/aevfavwt/r5w6pfu3J/Xf8h/X7hb+X+egNLCgDNU3cwOzExUXb96uqqIpFIvW8PONLK/XW9/9Fnev+jz3R35UHJY3dXHhQe4wcFgHI6oW6NRqPK5XJ67733Wl2UjpD9p0395699q/CX/afN6l/72yz4+b/sg40GlhQAmqfubsa5XG7Xx7ZPAg90kkqTa9+4tVwyfU6x0I35wv8vj48wpxmAHahbAdipmjmLJbLRw3lqDma9Xq9cLpdcLpcGtg0UtyxLuVxOo6OjdRewGaLRqK5evdrqYsBhmFwbgN06qW4F0D6YsxidquZgdnZ2VrlcTmNjY2UTVRiGoZMnT9ZVuGaJRqOKRqO6ffs2Y3sAAC3TSXUrAACNVnMw+81vflOS5Pf7NT4+bluBgE5w6cyQnjl9TNLWGNnirsWJS6M6MdgvqbRrMgBQtwIAUL26x8z+4he/sKMcQEcZPNSnwUN9ZR87MdivU0cP71i/cn9dN24tF5YvnRna9T0AdDbqVgAA9lZ3MAvAHvlsk3nPnD5GMAsAHejh6qo2s1lJW8l4ihUvH/R6dx3nCAAgmAUAAGiqtbfe1sobb5R9rDjL7ODkpB59+aWyzwMA2DDPLAAAwHbRaFQul4vEigCAhiGYBQAAtotGo8rlcnrvvfdaXRQAQIeimzEAAEATeS5e0JGnz0raGiNb3LX4iWvX1Dt0XNLWmFkAwO4IZgEAgOM8XF3V2ltvS5Kymwf0X4z/XQe+8ruS2j8bfM/AwK6JnXqHjqtvZKTJJQJ2Kj7HpK2bMCQkQ7shmMUOxVkW8yplW8wj62LjrNxfV/bBhqSteWuLFS97+3vb+gccANhlM5stJFFaPnxU/3F8qPAY2eCB+hWfY5J05Omz/M5D2yGYxQ6VsizmFXeJyiPrYuPcuLVcMm1PsdCN+cL/L4+P6AcTp5pVLAAAAKBlSAAlMi4CAAAAgNMQzIqMiwAAoL0deOSRissA0I3oZowdirMs5lXKtphH1sXGuXRmSM+cPiZpa4xscdfixKVRnRjsl7Q1ZhYA0Hl63O6KywDQjQhmsUOlLIt5ZFtsrsFDfbsmMzkx2K9TRw83uUQAgEYpTvqXVyn5317hGyY5AAAeO0lEQVSqeS0JBAE4EcEsAABAG6mU9C+vuIdOsT886dXfLWXLPlbptSQQBOBEjJkFAAC2I7kiAKDRCGYBAIDtSK4IAGg0uhmjoQ56vRqcnCxZBuBcnNNA4xUn/curlPxvL9W8lgSCAJyIYBYN1TMwoEdffmnvJwJwBM5poPEqJf3Lqyf5H4kDAXQKglkHe7i6qrW33i4sey5e2DMLcS2vQWdbub+uG7eWC8uXzgyR0RK2qvW6w/UKAABUQjDrYJvZrFbeeKOwfOTps3v+0KvlNehs2QcbJVkznzl9jGAWtqr1usP1CgAAVEIwCwAAAHQJch+gkxDMAgAAAF2C3AfoJASzAAAANtttzPf2PAXnj262oniAbchvgFYimAVaaOX+urIPNiRtTZ1QrHjZ29/LOFYAcJDdxnxvz1Mw8Z0hHWxFAQGbkN8ArUQwC7TQjVvLJT9qihXPCXh5fEQ/mDjVrGIBAAAAbe9AqwvQDqLRqFwul5588slWFwUAAAAAUAWCWW0Fs7lcTu+9916riwIAAAAAqALdjIEWunRmSM+cPiZpa4xscdfixKVRnRjsl7Q1ZhYAut1HH/6jPvnHFUnSb/7hf+mjw0clSf/QXzo+b3HpH3e8ltwDwN4erq5qM5uVJG0sf1jyWPHyQa+XcbFoCwSzQAsNHurb9cfVicF+nTp6uMklAoD29ZMbGV3/zPPlivE/K/u8F3/+vqT3S9aRewDY29pbb5ckcyp276Uvp/MZnJxkeh+0BYJZoMG8/b26PD5SsgwAnS4ajerq1autLgYAoIMxZhZosMFDffrBxKnCH93cAHQD8lEAABqNllkAAOAI37vk158Wj5n99/9e0taY2X/3L/9N4XlvPntKwyd/r+S19IoB9ua5eEFHnj4raWuMbHHX4ieuXVPv0HFJW2Nmm2Fzc1P37t3T559/rs3NzaZ8Jux34MAB9fb2yuPxyOPxyOVy2fbeBLMAAMARjh7/PR09vhWkrn/Qp3/22Udlnzd88vfIOQDUoGdgYNfETr1Dx9U3MlL2sUbY3NzUnTt3dP/+ffX09Kinh7DFqX7zm9/o17/+tSzL0ieffKLh4WHb9idHBQAAQIsceOQRDU5OFpab1eIFtLt79+7p/v37euyxx/TEE0/Y2pqH5nv48KHu3bun1dVVffzxx3r88cdteV+CWQAAgBbpcbvJCguU8fnnn6unp4dAtkP09PRoaGhIn376qT799FPbglkSQAEAAABoK5ubm+rp6SGQ7SAul0s9PT364osvbHtPglkAAAAAgOMQzAIAAAAAHIcxswAAAAC6wsPVVa299XZh2XPxwq4ZnNH+CGYBAAAAdIXNbFYrb7xRWD7y9FmCWQejmzEAAAAAwHEIZgEAAAAAjkMwCwAAAADYl+npaQ0PD8vj8SgcDrekDASzAAAAAICqJRIJxeNxJZNJLS0tyTRNTUxMNL0cJIACAAAAAFQtFospFovJ5/NJkpLJpDwejyzLktvtblo5CGYd4uHqqjaz2ZJ1G8sfVlyWpJzLJVcut6/XHPR6yerW4Vburyv7YEOSdHflQcljxcve/l4NHupratngfNuvV7Vcq6p9HdcrAEA55X47S9XVLXmNqmMikYhSqZRM05RhGIpEIgqFQpK2uu7G43Fls1n5/X5dv369EBwGg0F5vV7F43FJkmVZ8ng8mp+fLwSV4XBYw8PD8vl8CofDZd87/7mxWEyBQKBQpkQiIUm6cuWKpqamdi2/aZoyTVN+v7+wzu12yzAMzczMFD6vGQhmHWLtrbdL0oiXc++ll3as+8pTT+nzd97Z12sGJyf16Ms716Nz3Li1rB/NflD2sdCN+cL/L4+P6AcTp5pVLHSIva5XtVyrdnsd1ysAQDnV/HaWytcteY2oYxYWFpRIJLS0tCS3262FhQVlfxt0h8Nhzc3NKZ1Oy+v1KhKJaHR0VIuLi1W/fzab1dzcnF599VXFYjGdO3dO0pcBdDwe19jYmObm5rSwsCBpK0i2LEtLS0vKZrOamJiQz+crCVaLmaYpSTtaYA3D2FdZ7UAwKykajerq1autLkZbO+j1anBysmQZreHt79Xl8ZGS5WZYub+uG7eWC8uXzgw5rtW22m3ohG3tZg9XV7X21tuFZc/FC7TeoimKW4I2lj+U1XtIn/b1S5L+6fZd/c7aupZ/tVHymu29YyR6xaD9bT/WixUv03tmJ9M05S36HZ1vUTVNU4lEQmtra4UgMR6Pa3h4WIlEYl+tnQsLC1pcXJRhGJK2WnCnp6dLWnD9fr/8fr9M01QqldLi4qLcbrfcbrfi8bji8fiuwWw7IZjVVjAbjUZ1+/ZtPfnkk60uTlvqGRig9aNNDB7qa0lrafbBRklr7jOnjznux1a129AJ29rNNrPZkrvxR54+y48pNMX2lqD/8rVv6T9/7VtbC/91Q9LyjtcU94bJo1cM2l2lVs/ilk56z+yUDxA9Ho/8fr+CwaBCoZAymYwMw9jR2un3+5VOp/cVzPr9/kIgK0mZTEbSl4FzsXzr7OjoaMl6bw0NV9ky3bobjWDWITwXL+jI02dL1m0sf1hywXji2jX1Dh0veU65MbN7vYZW18536cyQnjl9TNJWq0Dxj6nEpVGdGNxqSWhWqy86y/brVS3Xqmpfx/WqfTmh15OTemHU0yunVT16gGLN7jVT7rezVF3dkteIOsbtdmtxcVGJRELpdFrhcFjz8/MaHh627TOKA9lqn7+f7sHFLb7FwbdlWbZuRzUIZh2iZ2BgzxO+d+i4+kZGKj7HjtfA+QYP9e36g+3EYL9OHT3c5BKhk+x1var1usP1ylmc0OvJSb0w6umV06oePUCxZveaqea3s9S6uiUUCikUCimVSun555/X7OysIpFIITlTXiaTKZnDtbj1s9qW0HyL7MLCwo7W2XxX4+2fW0m+Bbk42ZNlWTuSQjUDwSwAAIANiluCNpY/1L/+P1/Rv/qf/02SdPTf/lv9zuPHtPyrDU1m/lfhNcW9YfJoOUW7236s79bSSe+ZnfJZjAOBgLxer9LptAzDkM/nUyAQUDAYLGQwjsViymazhczCXq9Xc3NzhfeKRCJVfaZhGAqFQgoGg0omkzIMQ3Nzc4V5YqempkoeS6VSSqfThazJ5Vy5ckWRSERjY2MyDEPBYFCBQGDfrcL1IpgFAACwwfaWIPfGfbk37kuSjD84ob6REfV+9JlUFMzSGwZOVKnVk140lRmGoXg8rldffVWWZcnn8ymZTEramqs1EokoGAwWpuZZWloqvDYcDmtmZkYej0djY2OKRCJaWFioal7XeDyu6elpBYNBmaYpn8+nK1euSNqaM3Z6eloTExPKZrOF964kH2Dny3ru3LmKwW+jEMwCAAAAQBP4fD6l0+ldH4/FYorFYru+dm1trWTd9rGu+cC4nKmpqV3nj630WC3v1ywHWvrpAAAAAADUgJZZAAAAAF3hoNerwcnJkmU4F8EsAAAAgK7QMzDA3LcdhG7GAAAAAADHIZgFAAAAADgOwSwAAAAAwHEIZgEAAAAAjkMwCwAAAABwHIJZAAAAAIDjEMwCAAAAAByHYBYAAAAA4Dg9rS4AAAAAADTDyv113bi1XFi+dGZIg4f6Wlgi1INgFgAAAEBXyD7Y0I9mPygsP3P6GMGsg9HNGAAAAADgOASzAAAAAADHIZgFAAAAAOyLZVlKpVIaHh5uWRkYMwsAAAAAqFokElEikZBhGDJNs2XlIJgF2oS3v1eXx0dKlgEAALB/K/fXlX2wsWP93ZUHFZeLeft7SQ61i1gsplgsplQqpWAw2LJyEMwCbWLwUJ9+MHGq1cUAAEc46PVqcHJSknRg84BeNv43HfjK70riZiBgh+JzLL/sJDduLZdkLd5N6Mb8ro9dHh9pyG+zSCSiVCol0zRlGIYikYhCoZAkaXp6WvF4XNlsVn6/X9evX5fb7ZYkBYNBeb1exeNxSVvdfD0ej+bn5+Xz+SRJ4XBYw8PD8vl8CofDZd87/7mxWEyBQKBQpkQiIUm6cuWKpqambN/uRiCYBQAAjtMzMKBHX35JkvSopH/e2uIAHaf4HIN9FhYWlEgktLS0JLfbrYWFBWWzWUlbgejc3JzS6bS8Xq8ikYhGR0e1uLhY9ftns1nNzc3p1VdfVSwW07lz5yR9GUDH43GNjY1pbm5OCwsLkraCZMuytLS0pGw2q4mJCfl8Pvn9fvu/AJsRzAIAAABAE5imKW9RK3e+RdU0TSUSCa2trRVaYuPxuIaHh5VIJAqtq9VYWFjQ4uKiDMOQtNWCOz09XdKC6/f75ff7ZZqmUqmUFhcX5Xa75Xa7FY/HFY/HCWYBAAAAoNkunRnSM6eP7Vh/d+VBSdfixKVRnRjsL/sejRiykA8QPR6P/H6/gsGgQqGQMpmMDMMoBLLFz0+n0/sKZv1+fyGQlaRMJiPpy8C5WL51dnR0tGS91yHdyglmAQAAAHSUwUN9VSVvOjHYr1NHDzehRFvcbrcWFxeVSCSUTqcVDoc1Pz9v6/Q2xYFstc/fT1fmdtIR88ymUimNjo5qdHRU09PTrS4OAAAdgfoVABojFAopmUwqmUxqZmam0OV3+zQ3mUxGTz31VGE5P752+/8rybfI5lthi+32uU7h+GDWsixFIhHNz89rfn5eN2/eLLujAABA9ahfAcB+qVRK09PTMk1TlmUpnU7LMAz5fD4FAgEFg0EtLCzINE2Fw2Fls9lCZmGv11sSdEYikao+0zAMhUKhwntblqVMJqNgMCi3262pqSkFg8HCe6dSKYXDYfs3vgHaKphdWFjYtYk9kUjI4/HI4/GUfLmZTKak//f58+d18+bNhpcVAACnoH4FgPZgGIbS6bRGR0fl8Xg0NzenZDIpSUomk4VxtKOjo8pms1paWiq8NhwOyzRNeTweTUxMKBwOlx1nW048Hlc4HFYwGJTH41EkEtH58+clbc0Ze/78eU1MTMjj8Sgej+85d2wikZDL5So8z+VyyeVy1fq11KwtxsxOT08rEonIMIyyTdypVEqRSESzs7Nyu90KBoMKBoNKJpOFeZLy3G633nnnnWYWHwCAtkT9CgDtxefzKZ1O7/p4LBZTLBbb9bVra2sl67aPdc0HxuVMTU3tOn9spcfKCYVC+0pK1Sht0TI7NTWlXC63646LRCKKxWLy+XwyDEPJZFKpVEqWZUlS4V/JOZm3AABoNOpXACjl7e/V5fGRwl8jMhajedqiZbYSy7JkmmbJPEf55vRMJiO3210y+Hn7nWQAtVu5v67sgw1JW6nsixUve/t7q8oY2ArVbkOl9U7Z1m70cHVVm9sSYGwsf1hx+aDXq56BgYaXrd11W/1afC3Iq+aawDkP1K74Gl3p2tzM6/LgoT79YOJUUz4Ljdf2wWy+W9T2CjTfZSoUCikSiciyLLndbt28eVPXr1/f830//vhjffLJJyXr7ty5Y1/BgQ5w49ayfjT7QdnHiudouzw+0rYVQ7Xb8Icnvfq7pfJZAZ2yrd1o7a23tfLGGxWfc++ll0qWBycn9ejLL+3y7O7RiPq1nevWSteCvOJzPY9zHqhdpWt08bWZ6zJq1fbB7G4pp71er1ZXV+V2u5VMJjU+Pi5pK0FFuQmBt/vxj3+sq1ev2lpWAACcohH1K3UrAKCZ2j6YrYbf79f8/M67qZV8//vf35Gl686dO3r22WftLBoAAI613/qVuhWAXQ4cOKDf/OY3rS4GbPbFF1+op8e+ELTtg9ndEk5ks1kN1NG3/rHHHtNjjz1W8+uBbnDpzJCeOX1M0tZYsuIueIlLozox2C9JbZ08odptKObUbe1GnosXdOTpsyXrNpY/LOm+9sS1a+odOl5YPkgiI0mNqV/buW4tvhbkVXNN4JwHald8ja50bS53Xe7t7dWvf/1rPXz40NbgB62zvr6ujY0NHT582Lb3bPsjIz+WZ3viCdM0q5pTCUDtBg/17Zr45MRgv04dte9i1Ch2bINTtrUb9QwM7Jk0pHfouPpGRppUIufotvq10rUgj3MdsFela/Re12aPxyPLsnTv3j0NDQ21ZA5T2Gd9fV13796VtLVv7dL2wazb7ZZhGEqlUoW5j0zTlGVZOnfuXItLBwCAM1G/AmhnHo9Hn3zyiVZXV/Xpp5/SOutgX3zxhTY2trLJP/bYYzpy5Iht7+2IoyISiSgcDsvv9xcmdQ8EAh155xgAgGahfgXQrlwul4aHh/Xxxx/r008/1RdffNHqIqFGPT09Onz4sDwej44cOWJrK3tbBLOJRELhcLiwnN/AXC4nSQqFQrIsS+Pj47IsS4FAQMlk0rbPj0ajZF8EAHScVtav1K0A6tXT06PHH39cjz/+eKuLgjZ1oNUFkLYq01wut+Ov2NTUlNbW1pTL5WwNZKWtCjeXy+m9996z9X0BAGilVtav1K0AgEZri2AWAAAAAID9IJgFAAAAADgOwSwAAAAAwHEIZgEAAAAAjtMW2YxbbXvGxTt37rSwNNXbWF7WL9fXC8v/9MEH6v3tHE52vgad7e7KA218slxY/uB//Hf9ZqW/5ue1s27a1nZS63Wn2a/bj3w9sV70OSjlhLq1kef6bsch1xd0mmqvufwGRTX2W7+6ctvTGnaxv/7rv9azzz7b6mIAABzi5z//uf70T/+01cVoa9StAID9qrZ+JZgtYlmW/vZv/1a///u/r76+vprf586dO3r22Wf185//XF/96ldtLGFn4XuqDt9TdfieqsP3VJ29vqf19XX98pe/1B/90R/J7Xa3oITOQd1au27bZra3s7G9nc2u7d1v/Uo34yJut9vWO+xf/epX9Qd/8Ae2vV+n4nuqDt9TdfieqsP3VJ1K35PP52tyaZyJurV+3bbNbG9nY3s7mx3bu5/6lQRQAAAAAADHIZgFAAAAADgOwSwAAAAAwHEORqPRaKsL0Yn6+/v1x3/8x+rvJ91+JXxP1eF7qg7fU3X4nqrD99R+unGfdNs2s72dje3tbK3YXrIZAwAAAAAch27GAAAAAADHIZgFAAAAADgOwSwAAAAAwHEIZgEAAAAAjkMwa7NEIiGPxyOPx6NwONzq4rSlYDAol8tV8jc8PNzqYrWFhYWFXb8Ljq0v7fY9cWx9ybIshcPhiscMx9Te3xPHVHNVugaW4/RjeD/b6/RjsZprUjlO3ce1bK/T93Hx9gaDQVmWtedrnLp/pf1vr9P3b55pmkokEm2zfwlmbZRKpRSJRDQ7O6v5+XnNzc0pGAy2ulhtKRQKaXFxsfCXTqdbXaSWmp6elsvlUjAYlGmaOx7n2Nqy1/ckcWzl5Y+P+fl5Xb9+XTMzM5qYmCg8zjG1Za/vSeKYaoZqzu3tnHwM17K9krOPxWrOte2cvI9r2V7Juft4dHRUkgr7yrIsjY+PV3yNk/dvLdsrOXf/FgsGgwqHw8pmsxWf17T9m4NtDMPIxePxwvLi4mJOUm5tba2FpWo/gUAgNzU11epitKVkMpkrd1pybJXa7Xvi2NqyuLiYc7vdJevm5+dzknLz8/O5XI5jKper7nvimGqu3c7tcjrhGN7P9jr5WKzmXCvHqfu41u118j4u3k+5XC6XTqf33FdO3b+5XG3b6+T9mxePx3OBQCAnKbe4uFjxuc3av7TM2sSyLJmmKb/fX1hnGIbcbrcymUwLSwan49jCfhmGoevXr5es8/l8kqS5uTmOqd/a63tC++IYdpZazjUn7+NuvLaEQqGS5XQ6LbfbLbfbXfb5Tt6/0v63t1NEIhFduXJlz+c1c/8SzNok303IMIyS9YZh7KsLUbfIjxPKjzPA7ji29odja0sgECi73uv1ckwVqfQ95XFMtZ9uPYadfCxWc64Vc/o+3u/25jl5H+dlMhlNT09XDHqcvn+LVbO9eU7ev5FIRGNjY4UbM5U0c/8SzNpkt37jXq9Xq6urTS5N+8tkMorH45qfn5dpmlWNI+lWHFv7w7FVXiKRkNvtViAQ4JiqoPh7yuOYaj/degx30rFY7lwr1mn7eK/tzXPyPl5YWJDL5dLExIT8fr+mpqZ2fW4n7N/9bG+eU/evZVmanp5WLBar6vnN3L8Es2i68+fPK51Oy+/3yzAMxWIxZTIZx92JQ/vh2CrPNE1FIhElk8lWF6WtlfueOKbQLjrpWOy2a1K12+v0fezz+bS2tqb5+XlJckygVqv9bq+T9+/zzz+vQCBQVatssxHM2mS3biPZbFYDAwNNLk17CwQCJX3ox8bGJMkRYyRagWOrehxbO1mWpYmJCV2/fr3w3XBM7VTue5I4ptpVNx7DnXIs7naubdcp+7ja7ZU6Yx+73W75fD4lk0llMhmlUqmyz+uU/Vvt9krO3b8LCwtKpVI7xoFX0sz9SzBrk3yf8O13V0zT7PjB4PXKd0XYaxxJt+LYql23H1uWZWl0dFTxeLykaxvHVKndvqdyuv2Yahccw848FvdzrnXCPt7P9pbjxH2cl0+G9M4775R9vBP2b7G9trccp+zfmzdvSpJOnjxZmDNWkoaHhwtTFG3XzP1LMGsTt9stwzBK7siYpinLsnTu3LkWlqz9bL9rtbCwIElt2XWhHXBsVY9j60v5Oe/i8fiO1gCOqS9V+p4kjql21Y3HsNOPxb3Ote2cvo/3u72Sc/exaZo7ghbLsmRZloaHh8u+xsn7t5btlZy7f2OxmNbW1rS0tKSlpSXNzs5K2srgnP//dk3dv7ZO9NPl4vF4YQ6xxcXFnM/nywUCgVYXq63k516Lx+O5tbW13Pz8fM7tdudCoVCri9YWdptzkGOrVLnviWOrlM/ny01NTeXm5+dL/vLzu3FMban0PXFMNd9u18DFxcUd8zp2wjFc7fZ2wrG41zWp0/bxfrfXyft4bW0tZxhGLhaLFfaV3+/PGYZReE4n7d9attfJ+3e7/HyxxfPMtnL/EszaLBaL5dxud06SI07IVsif9JIKF4Nulz/ht/8V49ja+3vi2NqSr2jK/RV/J91+TFXzPXFMNcde53b+8e2cegzXsr1OPharOdc6aR/Xur1O3sdra2u5QCCQc7vdObfbnQsEAoXAPZfrrP2by9W2vU7ev8XKBbOt3L+uXC6Xq6lJFwAAAACAFmHMLAAAAADAcQhmAQAAAACOQzALAAAAAHAcglkAAAAAgOMQzAIAAAAAHIdgFgAAAADgOASzAAAAAADHIZgFAAAAADgOwSwAAAAAwHEIZgEAAAAAjkMwCwAAAABwHIJZAAAAAIDjEMwCAAAAAByHYBaAJCkYDCoSiSgSicjj8cjj8SiRSJQ8JxKJaHh4WC6XS8PDw0qlUi0qLQAAjWVZliYmJgp1Xjgc1vDwsDwejyQpk8lodHRULperbJ25W71a/L6jo6MyTbOu1+xVDqCTEcwCKJientbAwIDm5+d17tw5hcPhwmMLCwtKJBJKp9NaW1tTPB5vYUkBAGisYDAot9ut+fl5RSKRQh24tLQkSUqn07py5YrW1taUTCYVDoe1sLBQ8h7l6tXx8XHFYjHNz8/LsixFIpG6XlNNOYBO5crlcrlWFwJA6wWDQWUyGa2trUmSTNPU8PCw5ufn5fP5lEqlFIlEtLi42OKSAgDQeB6PR7Ozs/L5fIXlZDIpv99f9vmjo6Py+/2KxWKSdq9XY7GYpqamJG31eMpkMpqfn6/5NXuVA+hktMwCKBgbGyv83+v1ljyWr7xdLpcmJiboxgQA6Gh+v1/xeFyWZSmTyciyrJJ6cjvDMEq6/0rl69V8cCxJAwMDsiyr7tfsVQ6gUxHMAihwu90VH1tcXFQymZRhGAqHwyXdkAEA6DSmaerkyZMKh8NKp9Ml9aRpmgqHwxodHZXH4ymbR6Jcvbr9ZnG9r6mmHECnIpgFsC+BQEDxeFzJZFIzMzOtLg4AAA2xsLBQyBOxuLhY0r3YsiwNDw9rdHRUs7OzWltbUyAQaHoZ26UcQKsQzAKoSiqV0vT0tBYWFmSaptLptAzDaHWxAABoCMMwND09rUwmU6j78rLZrKQvW0wzmYwymUzTy9gu5QBapafVBQDgDIZh6ObNm3r11VdlWZZ8Pp+SyWSriwUAQEO43e4dmYYNw1AymZTP51MoFCpkPPb7/RXH0zaKYRhtUQ6gVchmDAAAABRJJBKFqXDyY1hN01QkEpFlWUqn0y0uIQCJbsYAAABACa/Xq2w2q7m5OVmWVfLHEBugfdDNGAAAACiST6IUiURkmmYhiA2Hw4X5XgG0Ht2MAQAAAACOQzdjAAAAAIDjEMwCAAAAAByHYBYAAAAA4DgEswAAAAAAxyGYBQAAAAA4DsEsAAAAAMBxCGYBAAAAAI5DMAsAAAAAcByCWQAAAACA4xDMAgAAAAAc5/8HcM+TnHcY6BQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(2, 2, figsize=(8,6))\n", "axs = np.ravel(axs)\n", "names = 'mlog10p ts ns gamma'.split()\n", "for (i, name) in enumerate(names):\n", " ax = axs[i]\n", " mask0 = trials.ns_0000 > 0 if name in 'ns gamma' else trials.ns_0000 >= 0\n", " mask1 = trials.ns_0001 > 0 if name in 'ns gamma' else trials.ns_0001 >= 0\n", " hl.plot1d(ax, hl.hist(trials[name+'_0000'][mask0], bins=25), crosses=True, color='C1', label='source 0')\n", " hl.plot1d(ax, hl.hist(trials[name+'_0001'][mask1], bins=25), crosses=True, color='C0', label='source 1')\n", " ax.set_xlabel(name)\n", " ax.set_ylabel(r'trials per bin')\n", " ax.semilogy()\n", "ax.legend()\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stacking" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simplest way to implement a stacking analysis is to simply add the TS values. Let's do some background trials for that setup:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Performing 10000 background trials using 15 cores:\n", " 10000/10000 trials complete. \n" ] } ], "source": [ "bg_trials = multr.get_many_fits(10e3, seed=1)\n", "bg_trials['sum_ts'] = np.sum([bg_trials[name] for name in bg_trials.keys() if 'ts' in name], axis=0)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "sum_bg = cy.dists.Chi2TSD(bg_trials.sum_ts)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAHUCAYAAADftyX8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAASdAAAEnQB3mYfeAAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde1xUdcLH8e+AggIqIOqa5mVQsiRTsMw2swK0stxaQbc2d2sz0MpbGeOW1aClorVm62phrbtrl0clu21WCpl2se3x0gVb8zJpaaYCYgIKCvP84eMkzYjADJy5fN6vF6/Hc86PM9/HF1vfDr/z+5nsdrtdAAAAQAAKMjoAAAAAYBTKMAAAAAIWZRgAAAABizIMAACAgEUZBgAAQMCiDAMAACBgUYYBAAAQsCjDAAAACFiUYQAAAAQsyjAAAAACVjOjA3iTkpISrVu3Tueff75CQ0ONjgMAAIB6qqio0Pfff6/BgwcrMjLynOMpw2dYt26dbr75ZqNjAAAAwE2vv/66fvOb35xzHGX4DOeff76kU395PXr0MDgNAAAA6mvnzp26+eabHb3uXCjDZzg9NaJHjx7q3bu3wWkAAADQUHWd8soLdJKsVqtMJpPi4+ONjgIAAIAmRBnWqTJst9tVUFBgdBQAAAA0IcowAAAAAhZlGAAAAAGLMizmDAMAAAQqyrCYMwwAABCoKMMAAAAIWJRhAAAABCzKMAAAAAIWZVi8QAcAABCoKMPiBToAAIBARRkGAABAwKIMAwAAIGA1MzpAoCssrVBxWWWj3Ds6PEQxEaGNcm8AAAB/QBnWqTnDWVlZhnz20g17ND9/R6Pce2JST01OiWuUewMAAPgDpkmIF+gAAAACFWUYAAAAAYtpEgYbPbCrhvXp2Cj3jg4PaZT7AgAA+AvKsMFiIkJ5yQ0AAMAgTJMAAABAwKIMAwAAIGBRhnVqNQmTyaT4+HijowAAAKAJUYbF0moAAACBijIMAACAgEUZBgAAQMCiDAMAACBgUYYBAAAQsCjDAAAACFiUYQAAAAQsyrBYZxgAACBQUYbFOsMAAACBijIMAACAgEUZBgAAQMCiDAMAACBgUYYBAAAQsCjDAAAACFiUYQAAAAQsyjAAAAACFmUYAAAAAYsyDAAAgIBFGRbbMQMAAAQqyrDYjhkAACBQUYYBAAAQsCjDAAAACFiUYQAAAAQsyjAAAAACFmXYR1VUVGjixIn6xz/+YXQUAAAAn9XM6ACov507d2rUqFHavHmzwsLCdNlll+miiy4yOhYAAIDP4cmwD/rmm2+0efNmSVJ5ebn+9Kc/yW63G5wKAADA91CGfdCwYcP0wAMPSJL69Omjf/7znzKZTAanAgAA8D1Mk/BRM2fOVLt27TRhwgS1bNnS6DgAAAA+iTLso0JCQmSxWIyOAQAA4NOYJuGHmD8MAABQN5RhP5Obm6urr75a5eXlRkcBAADwepRhP2G32/XEE08oLS1N69ev1x//+EdVV1cbHQsAAMCrUYYlWa1WmUwmxcfHGx2lwSoqKrRy5UrH8auvvqr//Oc/BiYCAADwfpRhnSrDdrtdBQUFRkdpsBYtWujNN99Ux44d1axZM7388ssaOHCg0bEAAAC8GqtJ+JFOnTrprbfe0sGDB3X99dcbHQcAAMDrUYb9TGJiotERAAAAfAbTJALEwYMH9f777xsdAwAAwKtQhgPAwYMHdc011+j666/XqlWrjI4DAADgNSjDfu7YsWO69tpr9fXXX6uyslK33HKLPvzwQ6NjAQAAeAXKsJ9r2bKl/vCHPziOBwwYoH79+hmYCAAAwHvwAl0AyMzMlCT9+9//1qpVqxQREWFwIgAAAO/Ak+EAkZmZqfz8fIowAADAGXgy7McKSytUXFb5i7PHncZ99slHunTgr2Uymep03+jwEMVEhHogIQAAgLEow35s6YY9mp+/o9YxP332mg6vfUGtL/utIq++s06FeGJST01OifNUTAAAAMMwTSKAlX6Vp8NrX5Ak/fTZSh3Oe9bgRAAAAE2LMhzAmrVuL1NIy1MHQc3UsscAYwMBAAA0MaZJ+LHRA7tqWJ+OtYy4Sl/dcYUybk/TIzOf1PXDb6nTfaPDQzwTEAAAwGCUYT8WExF6zhfd4oYO1hDbLrVq1aqJUgEAAHgPpkngrEW4oKBAq1evbuI0AAAATccvyvCcOXOUkpKixMREWSwWo+P4hf3792vYsGG64YYb9NxzzxkdBwAAoFH4/DSJkpISSdKaNWskSVFRURo1apQSEhKMjOXTKisrddNNN+m7776TJI0dO1adO3fWsGHDDE4GAADgWV71ZHjz5s2KjY11eS0nJ0dRUVGKiopSRkaG43xkZKRju2GbzaaSkhKZzeYmyeuvQkJCdMcddygo6NSPx/Dhw3X99dcbnAoAAMDzvKIMz5kzRyaTSWlpabLZbE7Xc3NzZbFYlJ+fr02bNmnjxo1KS0urMSYtLU2xsbFasWKFIiMjmyq637rvvvv05ptvauDAgVq6dKmjGAMAAPgTk91utxsd4rTc3FylpaXpl5FiY2NlsViUnp4u6dQT4NjYWB0+fNhRfEtKSmSz2ZSWlqYVK1Y0aJrE1q1bFR8fr4KCAvXu3dv9/4f8gN1ud7kr3dnOAwAAGKm+fc7rH/edLrnJycmOc2azWZGRkcrLy3Oci4yMVEJCghISEnjhy4POVoR///vfy2q1qrq62oBUAAAAnuH1Zfj0tIlfzgM2m82y2WzKzc3V5s2bHec3b96sxMTEJs0YaGbPnq1XXnlFWVlZGjZsmIqKioyOBAAA0CBev5pEcXGxy/PR0dEqKipSZmamLBaLLBaLiouLlZqa6phOUZuDBw/q0KFDNc7t3LnTI5n92f79+zVjxgzH8ZdffqmqqioDEwEAADSc15fhusjOzq739yxcuFBZWVmNkMa/dezYUevXr1dqaqr27dunFStWqH379kbHAgAAaBCvL8PR0dEuzxcXF6tt27YNvu8999zjtCLFzp07dfPNNzf4noGif//+2rRpkz766CNdccUVRscBAABoMK8vw6fnCttsthrzhm02m1tLqLVv354nmg1UWFqh4pMhuvDya7X9wFGn66vfflMXXNRbXbu7XjP6bKLDQxQTEeqpmAAAAOfk9WU4MjJSZrNZubm5TptrjBw50uB0gWnphj2an7/D5bWKH3fqxxcflCm4udpeN17hFw6q830nJvXU5JQ4T8UEAAA4J69fTUKS4wW5zZs3O9YSTk1N9djmGlarVSaTSfHx8R65X6CqOl6qQ6/PkqpOyF5ZrsK35upE8T6jYwEAAJyVV5ThnJwcxw500qm1bc9c3zY9PV3Z2dlKSkpSbGyszGazVqxY4bHPt1qtstvtKigo8Ng9A1FQs1C1NP+8rF2bK0apeXQnAxMBAADUzqt2oDMaO9DVTWFphYrLKs96fdXrr+rNlcu14O8vqVmzus/EYc4wAABwV337nNfPGYb3iYkIrbW0xmXcoUkZd7i8tmvXLrVr106tW7dupHQAAAB15xXTJIzGnOGmUVZWphtvvFGJiYnauHGj0XEAAAAowxJzhpvK+PHjtW3bNu3cuVMDBw7Uxx9/bHQkAAAQ4CjDaBLHjh2TzWZzHA8YMEADBgwwMBEAAABlGE2kZcuWys/PV1ZWlqKjo7V06dJ6vVwHAADQGCjDaDLBwcF69NFHZbPZ1L17d6frpaWlOnLkiAHJAABAoKIMixfomlqbNm1cnp80aZIuueQSffjhh02cCAAABCrKsHiBzhu89tpreuGFF7Rnzx5dffXVWrJkidGRAABAAKAMwyssW7bM8edWrVrp2muvNTANAAAIFJRheIWXX35Zf/nLXxQSEqKFCxeqa9euRkcCAAABgDIMrxAUFKTJkydr+/btuu2225yu2+127d+/34BkAADAn7G2lU7NGc7KyjI6RsArLK1QRYtobT9w1Ona26/n6uHJ9+nBR6brupF/0Ntf/ei4dtMl5ykqLKTWe0eHh9S6hTQAAAhMJrvdbjc6hLfYunWr4uPjVVBQoN69exsdJ+DMW7Nd8/N3OJ0/WVqs/S/co+rjpZKksAuvUrvhmfW698SknpqcEueRnAAAwHvVt88xTQJe79iu/3UUYUkKi7vCwDQAAMCfME0CXq/VJUPVrE0HFa2ar9DOFym815VGRwIAAH6CMgyvMXpgVw3r0/EsV6/S0Rl/VHV1tdpERjld3fvdHnU6v4tMJpPL744Or31OMQAACEyUYXiNmIjQ2l9y69DK5ek9e/bo5qQrlJycrL/97W/q2PFshRoAAKAm5gzDp1VXV+tPf/qTjh49qtdee029e/fWoUOHjI4FAAB8BGVYp5ZWM5lMio+PNzoK6slms+nLL790HKempqpdu3YGJgIAAL6EMqxTZdhut6ugoMDoKKinHj166Ouvv9att96qTp06ae7cuUZHAgAAPoQ5w/B57dq108svv6xDhw6pTZs2Ttd37NihoKAgxcbGGpAOAAB4M54Mw2+4mh5RVVWl0aNHq0+fPpo/f76qqqoMSAYAALwVT4bh1+bPn6///Oc/kqRJkyapsLRSo8eMdTn2cHml3vriB8cx2zwDAOD/KMPwa/v27XP8uVnUefpncayWzltfp+998dPvzjmGbZ4BAPBtTJOAX3vqqae0du1axZzXVW2vn6Cg5jzFBQAAP6MMw+9dffXVsrywSi3Od146z36yUiUfv6LqinIDkgEAAKMxTUKnllbLysoyOgYa0R1Xxmp4v/Odzj+dPUOLPnpJYd+u0wPWbB2Oudhxra5zhgEAgO8y2e12u9EhvMXWrVsVHx+vgoIC9e7d2+g4aGQFBQXq16+fTp48KUlKTEzUZ599pqAgfmECAICvqm+f49/6CFgtW7bUlVdeKUkKDg7W888/TxEGACDAME0CASs2Nlbvv/++lixZov3796tv375GRwIAAE2MMoyAZjKZ9Kc//ems1++//37FxMTowQcfVPPmzZswGQAAaAr8Thg4i48//ljz5s3Tww8/rMTERH311VdGRwIAAB5GGQZcqK6u1tixP+9Ut3PnTrVq1crARAAAoDFQhgEXgoKC9Pzzz+vii08ttfbYY4+pW7duxoYCAAAeRxkGzmLAgAHatGmTFi5cqPvvv9/lmDO3ewYAAL6HMgzUonnz5ho3bpzLl+feeOMNxcbGKisrSxUVFQakAwAA7qIMAw1QWlqq8ePHq6KiQlarVddee63YvwYAAN9DGdap7ZhNJpPi4+ONjgIfkZeXp7179zqOx4wZI5PJZGAiAADQEGzHfAa2Y0Z9bNy4URkZGYqIiNAHH3xw1jJcWFqhpRv2OI5HD+yqmIjQpooJAEBAqW+fY9MNoIH69++vVe9/qG/3HtCOg6VO199f/Y46nd9Foe26aX7+Dsf53ue1VreY8FrvHR0eQmEGAKAJUIYBN7z82d4aRfe0qrIS/fD8WFVXHlPr/r9Rm1/fqqCQlpKk9KWbznnfiUk9NTklzuN5AQBATcwZBhrB4bUvqPp4qVRdpZ8+W6mKfduMjgQAAFygDAMeZrfbFdwqRgoKliSFxV2hlt37GZwKAAC4wjQJwA2jB3bVsD4dnS/cP1jb//uA5kyfpikzntJnh35+ue6mS85TVFiIqqurFRTk+r9Ho8NDGisyAAA4A2UYcENMROhZX3SL6zBAN169RpI0/BfXduzYoVuGD9ecOXN00003NXJKAABwNkyTAJqY3W7Xfffdp23btmn48OG65ZZbVF1dbXQsAAACEmUYaGJr1qzR6tWrHcedO3c+63QJAADQuPg3MNDEkpOTtXjxYrVt21YdOnTQjBkzjI4EAEDAogwDTSwoKEhjxozRN998o5UrVyoyMtJpzNatW/X9998bkA4AgMBCGQYM0rZtW11xxRVO56uqqjR69GhdeOGFmjt3rk6cOGFAOgAAAgNlGPAyixYt0pYtW1RWVqbMzEzl5OQYHQkAAL9FGQa8zI4dP2/vHBcXpzFjxhiYBgAA/0YZlmS1WmUymRQfH290FEDz58/X+vXrdfHFF2vBggUKDXW9jjEAAHAfZVinyrDdbldBQYHRUQBJ0qBBg7RlyxalpKQ4XSsrK9O1116rt99+24BkAAD4F3agA7xUcHCwy/OzZs3S2rVrtXbtWg257gZdnTFDLcIjJJ3aHvpsO+IBAABnlGHACxWWVqi4rNLp/Pd7dmvO3LmO430HC/Xshh9kMpkkSb3Pa61uMeG13js6PITCDADA/6MMA15o6YY9mp+/w+m8vbpKra65WyXr/qXqijIV97ldIf9fhCUpfemmc957YlJPTU6J82heAAB8FWUY8CGmoGC16nu9wi64Usd3b1FIB7PTmKqyElWVH1FIu64GJAQAwLfwAh3gg4JbtlL4hVe5vHZ43T+0f8l4FecvVnVFWRMnAwDAt/BkGPBCowd21bA+Hc85bndhWY2pEZP7SJOy8yRJRze+oaFxbfTEXxbU+J7o8BDPhgUAwIdRhgEvFBMRWqeX3KLDQzQxqafjuHLre2rWrJlOnjypli1b6qlZM9SlQ6vGjAoAgE+jDAM+LCYitObLcClxuum6FE2YMEGDBw9Wly5djAsHAIAPoAwDfqZXr1567733VF1d7XTNbrfr1ltv1aBBg5SRkaFmzfhHAAAgsPECHeCHTCaTy007Xn/9dS1btkz33XefEhIStG3bNpffX1haoXlrtju+CksrGjsyAACG8Nhjoc8//9zl+b59+3rqIwC4obKyUvfff7/j+If9P6osOELbDxx1Gru7sKzGOsds5gEA8Fdul+HFixdr7Nixkk79CvZMJpNJVVVV7n4EAA8ICQnRc889p9F3jdXBvd/K3v9WjXjhizp9L5t5AAD8ldvTJCwWi+6++27t3LlThw8frvFVXFzsiYwAPGTIkCF6MOcttb1hsiL6JDtdt1dXqXz7J07/YQsAgL/yyDSJqVOnqlu3bp64FYBGFtysuSIuTnJ5rfTzd1S85lmFdrpQUUnpCu3Y0+U4AAD8hdtlOD09Xa+++qoeeOABT+QB0MjOtqHH4eIiDX3udklSxb7/qmjVPHX80wKZTEHKGZ1YpznDAAD4Go88GX7wwQc1c+ZM9e/fX5GRkTWuLVu2zBMfUSuLxaK8vDyVlJQoIyNDmZmZjf6ZgK8624Ye//vdNrVu1UpHSkokSXfe/6h69b9AkpTQNYqX4wAAfsntMmyz2ZSamuo4buq5hrm5uSopKdGmTade8ImNjVVCQoKSk53nQwI4u0svvVTbtm3TnDlztG3bNj3357ucxpw8eVLBwcEymUwGJAQAwPPcLsPLly/3RA5J0ubNm5WWlqZdu3Y5XcvJyZHFYpEkjRw5Us8995wkKTk5uUbxTU5O1ubNmynDQAOEhYXJarWe9T9qp02bpi+++ELz5s1Tr169mjgdAACe5xWbbsyZM0cmk0lpaWmy2WxO13Nzc2WxWJSfn69NmzZp48aNSktLkyRFRkbWmJqRl5dHEQbc5OrJ786dOzVv3jy9++67uvjii5WTk2NAMgAAPKveT4aHDh2qtLQ0jRkzRpI0bty4WscvWrTonPfMzMxUZmamcnNzHSX3TBaLRdnZ2UpISJAkrVixQrGxsSopKalRhC0Wi1JTUx3jAHjO9OnTVVlZKUmqrq7W5ZdfbnAiAADcV+8yXFRUVOuxp5WUlMhms9V42ms2mxUZGam8vDzHfGWLxaLY2Filp6c3ah4gUM2fP1/R0dFasGCB0tPT1adPH6MjAQDgtnqX4Y0bN9Y49uScYVdOT5swm801zpvNZse1jIwMpaWlMT0CaERRUVF6+umnlZ6ervbt27scs2jRIg0bNkxdunRp4nQAADSMR5ZWa0xn28UuOjpaRUVFys3NVU5OjvLy8hzX6rK82sGDB3Xo0KEa53bu3Ol+YMDPXXTRRS7Pb9iwQffcc49atGihBx98UBaLReHhta9NDACA0TxShleuXKnnnnvO8dS4f//+mjNnji655BJP3L5WqampDVrObeHChcrKymqEREDgqa6u1oQJEyRJx48f11/+8heNGzeOMgwA8HpuryYxd+5cpaamqk2bNpo9e7Zmz56tbt26KSEhQa+99prbAaOjo12eLy4uVtu2bRt833vuuUcFBQU1vl5//fUG3w8IZOXl5YqPj3ccP/TQQ+rY0XmXOwAAvI3bT4ZnzZqlOXPmaMqUKY5zd999t1JSUpSZmalbbrnFrfufnitss9lqzBu22WxOu93VR/v27c867xFA/URERGjJkiW677779OSTT2ry5MlOY+x2u3744Qd16tTJgIQAALjmkXWGz9yB7sxzntiNLjIyUmazWbm5uY5zNptNJSUlGjlypNv3B+A5iYmJeuWVV9SyZUuna7m5uYqNjdXDDz+so0ePGpAOAABnbpfhP//5zzVeXjtt9+7dSkxMdPf2kk4tm2axWLR582bZbDalpaUpNTXVrSfDZ7JarTKZTDV+zQvAc44fP64HH3xQFRUVmjlzpq655pom37odAABX6j1N4pebbBQVFenVV1/VmjVrapzPy8tzWg7tbHJycpSRkeE4Pr371el/Waanp6ukpERJSUkqKSlRamqqVqxYUd/oZ2W1WmW1WrV161YKMdAI3n//fX333XeO47hrR+rpvB0aPbCrYiJCDUwGAAh0Jns9H8/Ud2pCY69D7Emny3BBQYF69+5tdBzAr3z++efKuGe8Pt99QL8a/ZRMpiCtnnyV4jq0MjoaAMCP1LfP1fvJsC+VWwBNp7C0QsVllWe9HtYxVtNzlmvM4nUymU7N0NpdWOa4/uILz+nH/T9o7MQHFNGqteN8dHgIT48BAI3G6zfdaApWq5U1hwE3Ld2wR/Pzd5xzXHDLn58Epy/dJEmqKj+ifTlW2SvK9MLfl6jt0HsVFjdQkjQxqacmp8Q1TmgAQMDzyGoSvs5qtcput6ugoMDoKEBAKvnoJdkrTj0lri4vUVAom3UAAJoGZRiA4Vr1vU4tuvaRJLWMG+j4MwAAjY1pEgA8YvTArhrWp/Zd53YXljmmRkhSzuhEdYsJl3SV7PY/6v3V76jnBb3UpdvPK9FEh4dIkt58800NHjxYbdq0aZT8AIDARBkG4BExEaH1ftGtW0x4jdUkLvjDKJfjvvnmG40YMUJRUVGaMWOGxowZo+DgYLfyAgAgMU1CEptuAE0lOjxEE5N6Or5OP/U9lylTpujkyZM6dOiQ7r33Xu3atauRkwIAAoXbZfj555/XTz/95IkshuEFOqBpxESEanJKnOOrLk+Sy8vLVVFR4Ti+5557FBfH6hIAAM9we5rE7NmzZTKZdNddd3kiDwDUEBYWpvfee0+rVq3S448/rscee8xpTGFphRbnbVWL8AhJYmc7AECduf1kePbs2crMzNSePXs8kQcAnJhMJg0bNkwbNmxQ27Ztna6/szpfD48apOkzZ+np97bWuvkHAABncvvJ8MaNG9W9e3eZzWYlJyfLbDbXuL5o0SJ3PwJAADvXznbV1dWaPs0ie2W5Sj74h8q+Xi/bHz+s073Z3Q4A4HYZttlsMpvNjhJcVFTkdqimxg50gPc61852x3Z/roP//Xm+f/hFgzX2pS11uje72wEA3C7Dy5cv90QOQ1mtVlmtVm3dupUVJQAf07JbX3W4dZYOr31B1ceOqnXiTUZHAgD4EI8srbZy5UqNHDmyxhve48aN0/vvv++J2wNArVp0uVi/+sNf1OHWWTI1c16urezrdTr8wRJV//+WzwAAnOb2k+G5c+cqJydHzz77rIYMGeI4n5SUpOzsbF177bXufgSAAFafne2atWkv6cyd7aRj5eW67soMFf6wT8E7PtC0J+boxlvSJKnO6xwDAPyXR5ZWy8/PV9++fWucT0xM1KhRrneTAoC6cndnu8cfn68ff9gnSTpcXKTO0RE1dr0DAAQ2t6dJ2O12RUVFOZ232Wzq3r27u7cHALcMHz5cycnJkqSBAwcqLS3N4EQAAG/i9pPhu+++W+np6VqxYoXj3LfffquxY8dq7Nix7t4eAM7p9DbPZx6f1qdPH61evVrvvPOOOnToIJPJ5PT9q1atUr9+/dSxY+3TMQAA/sftJ8PZ2dnq1q2bIiMjZbfb1bNnT/Xo0UPJycmaMmWKJzI2OqvVKpPJxEoSgI861zbPJpNJN9xwgxITE52+98cff9SoUaPUs2dPzZgxQ+Xl5U0VGwDgBUx2u93uiRvZbDZt2XJqbc+EhASfnCJxemm1goIC9e7d2+g4AJpAenq6Fi9e7Dj+5JNPNHDgQAMTAQDcUd8+5/Y0idPO3HgDAHxBdXV1jSfBaWlpFGEACDAeW2d46NChatu2rdq2bauhQ4fqiy++8MStAaDRBAUF6cUXX9Qnn3yiwYMHa/bs2S7H+eLOmgCAunG7DE+dOlVpaWmy2+2aPXu2Zs+ererqaiUkJOi1117zREYAaFQDBw7UBx984PK3W1u3blXnzp01efJkFRcXG5AOANCY3J4zHB0drYyMDM2aNavGeYvFopUrV2rHjh1uBWxKzBkG8EtJQ4bq/TWrJUmtottp67btOr9dpMGpAABnU98+5/aT4dNl+JceeughfrUIwKdt27ZN697PdxybYq/UsepgAxMBADzN7TKcmpqq/Px8p/P5+fnsQAfAp/Xq1Utvrf1ULWMvVVCLCLW+gn+mAYC/cXs1CZvNprlz52r16tWOcyUlJcrLy1NycnKNQrxs2TJ3P65RWK1WZWVlGR0DQBMrLK1QcVllrWOCozurfepjOnm0SMEtIrS7sKzG9Q/X5mn5i//QAw9b1c3cQ9KpTT/qu4U0AMAYbs8ZHjlyZJ3HLl++3J2PanTMGQYCy7w12zU/v+HvNdirq7R/yQSdKNwjBQWrzeUjFTno95qY1FOTU+I8mBQAUFdNvs6wtxdcAGgsZQX5p4qwJFVXKSistbGBAAD15pF1hgEgELWMvUytEoZJpiA1izpPrfpeZ3QkAEA9eWwHOgDwNaMHdtWwPh1rHbO7sEzpSzc5jnNGJ6pbTPgZI4bLtnOHSg4XK+HSAZJOzRk+7Y033lDr1q11zTXXeDQ7AMAzKMMAAlZMRGi9X3TrFrAvmBUAACAASURBVBOuuA6tapyL65DgcmxpaakyMjJ04MAB3XjjjZo7d6569erV4LwAAM9jmgQANJInn3xSBw4ckCT9+9//1vbt2w1OBAD4JcowADSSmJgYtWnTRpJ01VVX6aabbjI4EQDgl9wuwytXrtTKlSsdx+PGjVPbtm116aWXavfu3e7eHgB81n333addu3Zp0qRJevLJJ2UymZzGbNmyRVVVVQakAwBIHijDFotFkZGRkqTFixcrJydHOTk5SkhIUFpamtsBAcBI0eEhmpjU0/F15stxddG2bVvNmzdPl156qdO1vXv36oorrlBCQoLWrFnjqcgAgHpw+wW6Xbt2qX///pKkFStWKD09XSNGjFC/fv3Us2dPtwMCgJFiIkIbbQONRx55RMePH9eXX36pIUOGaNu2bbrgggsa5bMAAK65/WTYbDbr8OHDkqS8vDylpKRIko4cOeJ4YuztrFarTCaT4uPjjY4CIECUl5frk08+cRzffvvtFGEAMIDbZfjMp8Bms1m//e1vJZ0qxklJSW4HbApWq1V2u10FBQVGRwEQIMLCwvTVV19p/vz56tSpkx5//HGX4yorK5s4GQAEFrfLcHZ2tvLy8jR79mxt2vTzwvRms1l//vOf3b09APitkJAQTZgwQd9++626du3qdH3z5s3q2rWrnnvuOZ08edKAhADg/zyytFpCQoJGjBjhWEJI+vmJMQCgds2bN3c6Z7fbNWXKFP34448aO3asLr/8cladAIBGUO8X6MaNG1ev8YsWLarvRwCA3yssrdDSDXscx6MHdq2xG97HH3+stWvXOo6HDh2q4OBgt+8LAKip3mW4qKioMXIAQEApLqvU/PwdjuNhfTrWKK2//vWv9cYbbygzM1PFxcWyWCweuS8AoKZ6l+Hly5c3Rg4AwBlMJpOGDx+u66+/Xtu3b1fr1q2dxvz73//WBx98oIceekjR0dEGpAQA38d2zADgxZo3b67evXs7nT958qSmTJmip556SrGxsVq8eLEB6QDA97m96YYkff7551q2bJlsNpvTtWXLlnniIwAAZ3jhhRf0zTffSJJKSkpkt9sNTgQAvsntMrx48WJlZGQoISFBZrNZr776qkaMGKG8vDyX248CgL8rLK1QcVnt6wPvLiyr9fhcOvfqp8FJQ7Quf7Vie16gK4elafuBow2+b3R4CHOLAQQkt8vw1KlTtWbNGscGG9HR0Vq+fLny8vL4tR2AgLR0w54aL7HVRfrSTeccM6B7tP7zbfHPJ/pPUId2V6vUZNINf/3EaXz5jv/ojidtCulgPue9Jyb1bLRtpwHAm7ldhg8fPlzjCbDZbNbRo0d16aWXatSoUe7eHgBQixZd+7g8X11RpqJ35qv62FGF975akVeNVrPW7Zs4HQB4P7dfoEtISFBeXp7juH///lq+fLny8/NVUlLi7u0BAA1w5NNcVR/7SZJdZVvX6kTh90ZHAgCv5PaT4ezsbK1Zs0a//e1vJUkWi0WxsbEymUxKT093OyAA+JrRA7tqWJ+OtY7ZXVhWY2pEzuhEdYsJd/uzT9/XFNxcpmYhsp+s1BVXXaO/PzNJJpPprN8XHR7i9mcDgC9yuwwnJSU55gtLUvfu3XX48GEVFxere/fu7t4eAHxOTERovV9G6xYTrrgOrTyWIfLK2xTRJ0UDivM0zfKALviV8zrFW7Zs0SWXXKKgIFbZBBC4GuWfgG3atPGpImy1WmUymRQfH290FADwmGat22nW/EXq16+f07Vvv/1Wl19+ufr376/8/HwD0gGAd6j3k+GhQ4cqLS1NY8aMkSSNGzeu1vGLFi1qWLImZLVaZbVatXXrVgoxgIAwbdo0VVZWasuWLRoyZIh27dqlbt26GR0LAJpcvctwUVFRrccAAO929OhRbdiwwXF85513UoQBBKx6l+GNGzfWOF6+fLnHwgAAGl+rVq309ddfa+HChfrLX/6irKwsl+PKy8sVFhbWxOkAoGm5PWf4+eef108//eSJLACAJtKiRQvdf//9+vbbb9WpUyen6xs2bFCXLl309NNPq6KiwoCEANA03F5NYvbs2TKZTLrrrrs8kQcAAkJ0eIgmJvWscWzEfZs3b+50zm63KzMzU0VFRZo8ebKWLFmiLVu2sOoEAL/kkTKckZGh5ORkde3a1ROZAMDvxUSENsr2x56470cffaSPPvrIcTxixAiKMAC/5XYZ3rhxo7p37y6z2azk5GSZzeYa131hNQkAwM8GDRqk9957T5mZmTpw4IDuv/9+oyMBQKNxuwzbbDaZzWZHCWZ1CQDwfUOGDFFycrJsNpsiIiKcrufm5urdd99VVlaWyznHAOAr3C7DrCYBAP4pKChIPXr0cDpfWVkpi8Uim82ml19+WTNnztSkSZMa/DmFpRVaumGP43j0wK713sGvKe8LwL+4XYbff/99XXvttU7nP//8c0lS37593f0IAIAXefbZZ2Wz2SRJx44dU1RUlFv3Ky6r1Pz8HY7jYX06eqS0NtZ9AfgXt9+ISElJcXm+qKhIFovF3dsDALzMTTfdpNtuu02S1KdPH91+++0GJwKAhnP7ybDdbj/rtV9u0AEA8H3du3fXSy+9pPvvv18nTpxQcHCw05iXXnpJ5513nq655hoDEgJA3TW4DEdHR8tkMslkMqlt27Y1rpWUlMhutysxMdHtgAAA73S2f8YXFxfr3nvv1ZEjR3Tdddfpqaee0kUXXdTE6QCgbhpchvPz82W329W/f3+XL9GZzWZ1797drXAAAN8za9YsHTlyRJL07rvvaurUqQYnAoCza3AZ7tevnyQpOTlZSUlJHgsEAPBtXbp0UVRUlA4fPqwbbrhBgwcPNjoSAJyV23OGV69e7YkcAAA/MX78eN1+++3Kzs7W9TenavuBo05jPlybp/4DrlDLsDDtLiyrce2Xxw3V0PtGh4ew6gQQQNwuw97CZrMpOztbzz33nNFRACDgRUVFafbs2Zq3Zrvmz1tf49qJ4n364YV7FBzWRm1+fZsi+qTIFPTzS3jpSzed8/4DukfrP98W1ytTXe4rSROTejbKVtkAvJNflOGUlBQVFxcrOjra6CgA4Lc8tYlFyfqlUnWVqkqLVbx6oVp07aPmUed5MioA1Jnb6wx70ubNmxUbG+vyWk5OjqKiohQVFaWMjIwa19asWaM///nPTRERAALW6U0sTn8Vl1XW+x72k5WqOlroOI7oex1FGIChvOLJ8Jw5c2SxWGQ2mx27Gp0pNzdXFotF+fn5ioyMVFpamtLS0rRixQoD0gIA6mr0wK4a1qdjjXP2KUnKe+ffWvT0k3run/NVFhRRYwpDzuhEhRwv1nmdz3frs3cXljndt1tM+Dm/Lzo8xK3PBeBbvKIMZ2ZmKjMzU7m5uUpLS3O6brFYlJ2drYSEBEnSihUrFBsbq5KSEkVGRjZ1XABAHcVEhLqcSnHBnbfpnjtulclkcnrB7sA3m3VH2k268847lZWVpU6dOnkkS7eYcMV1aOWRewHwH141TcKVkpIS2Ww2JScnO86ZzWZFRkYqLy/PwGQAAHeYTCanc3a7XU8+/qiqq6v1wgsv6NJLL1VlZf2nYwBAXXl9GT49bcJsNtc4f7YpFQAA33Xi0G79t+BLx/HEiRMVEsK0BQCNxyumSdSmuNj10jnR0dEqKiqSJGVkZGjjxo2y2WxKSUmpMaXibA4ePKhDhw7VOLdz507PhAYANEhI++5656NNWvJMttavX68JEyY4jamurpbdbldwcLCLOwBA/Xh9Ga6LhqwtvHDhQmVlZTVCGgCAO87v2k0vvviijh49qpYtWzpdf/HFF/Xkk08qOztb1113ncvpFgBQV14/TeJsawcXFxerbdu2Db7vPffco4KCghpfr7/+eoPvBwDwrFatnF92O378uB555BF99dVXuuGGG3TXXXcZkAyAP/H6J8On5wrbbLYa84ZtNptbK0m0b99e7du3dzsfAPiDwtKKc64b3JDtjT29tfGSJUv03XffOY6HDBnisXsDCExeX4YjIyNlNpuVm5urzMxMSaeKcElJiUaOHGlwOgDwD0s37NH8/B31+p66bG/s6a2N77rrLlVVVWn69Onq0qUL/x4A4DavnyYhnVpn2GKxaPPmzbLZbEpLS1NqaqrH1hi2Wq0ymUyKj4/3yP0AAI0jJCRE9913n3bt2qVXXnlFQUHO/xp76qmn9Nhjj6m09KiLOwBATV5RhnNycmQymRwbbphMphovRKSnpys7O1tJSUmKjY2V2Wz26O5zVqtVdrtdBQUFHrsnAKDxtGrVSj179nQ6f+jQIVmtVk2fPl0pAy5R+fZPDEgHwJd4xTSJ9PR0paen1zrm9C51AADPc7Vt8i81ZHvjumxtHB0eoolJPWscN9SMGTNUWloqSSouKtT4qy5W99493b4vAP/lFWUYAGCss22bXBtPbW8cExHqsXnFt99+u7788kutW7dOt9xyi56Z9DuP3BeA//KKaRJGY84wAPiHyy67TGvXrtWqVas0e/Zsl2P++te/asuWLU2cDIC3ogyLOcMA4E9MJpOuv/56xcU5P23++uuvNWnSJCUkJOj222/X3r17DUgIwJtQhgEAAeOhhx5SdXW1JGn58uU6ceKEwYkAGI0yDAAICNXV1erVq5datGgh6dROpN27dzc4FQCjUYYBAAEhKChIs2fP1o4dOzR27Fg9/PDDTmPsdrteffVVVVbWvhsfAP9BGRYv0AFAIOncubMWLVqkdu3aOV1bs2aNUlNTdeGFF+rll192TKkA4L8ow+IFOgDAqWkUU6dOlSTZbDZNmDDBsWYxAP9FGQYAQNKBAwd08uRJx/G0adPUunVrAxMBaAqUYQAAJHXs2FFbtmzRv/71Lw0ePFjjxo1zGnPy5Ent2LHDgHQAGgtlGACA/xccHKzRo0frgw8+UGio8458S5Ys0YUXXqj09HTt27fPgIQAPI3tmHVqznBWVpbRMQDAq0WHh2hiUs8ax4GkvLxcVqtVVVVVWrx4sTZu3KhNmzbJZDIZHQ2AG3gyLF6gA4C6iIkI1eSUOMdXTITzk1N/9vHHH+vQoUOO42nTplGEAT9AGQYAoA5SUlK0bds23Xbbbbr88st1yy23OI05duxYjZfwAHg/yjAAAHVkNpv10ksv6YMPPnD5VPiJJ57QxRdfrJUrV8putxuQEEB9UYYBAKgnVy/X7d+/X/PmzdO2bds0YsQITZ482YBkAOqLMgwAgAcsWrRI5eXljuPRo0cbmAZAXVGGAQDwgEcffVQ5OTk677zzNHLkSCUmJjqNYeoE4H0owzq1moTJZFJ8fLzRUQAAPqpZs2a6++67tWPHDv31r391OWbMmDEaP368Dhw40MTpAJwNZVgsrQYA8JywsDC1b9/e6fwXX3yhJUuWaMGCBYqNjdWbb75pQDoAv8SmGwAANIFp06Y5pklUVVUpISHBaUxhaYWWbtjjOB49sGvArecMNDXKMAAATWDu3LkKCQnRypUrNX78eHXu3NlpTHFZpebn73AcD+vTkTIMNDLKMAAATaBXr1569dVX9dlnn6lHjx5O1+12ux59cKLKj3VRS3N/drcDmghlGACAJnTZZZe5PL9q1SotW7pEkhTaubdibry/KWMBAYsX6AAAMFhVVZWmTp3qOD5xeJ+CWrY2MBEQOCjDAAAYzG63a+zYsWob006SFHnF7xQU0tLgVEBgoAyLdYYBAMZq1qyZ7r33Xq35zxeKuuYuRVwy1GlMZWWlHn/8cRUVFRmQEPBflGGxzjAAwDuEh0eo9WW3yBTc3Ona4sWL9cgjjyg2NlYzZ85UZWWlAQkB/0MZBgDAy5WWlmr69OmSpCNHjuiVV15RcHCwwakA/0AZBgDAyx05cqTGKhSzZs2iDAMeQhkGAMDLderUSW+99ZbWr1+v8ePHa9iwYU5jioqKtGbNGgPSAb6NMgwAgI8YNGiQnnnmGZcbcmRnZ2vIkCFKSkrSZ599ZkA6wDex6QYAAD5u7969+utf/ypJev/995WVlaW33367ST67sLRCSzfscRyPHtiVLaThUyjDAAD4uP/+978KCwvT8ePHJUkzZ85sss8uLqvU/PwdjuNhfTpShuFTKMMAAPi4lJQU2Ww2PfXUUzpw4IAuueQSpzH79u1TaGioYmJiDEgIeC/KMAAATaCwtELFZbWvDby7sKzWY1eiw0MUExGqNm3aOJZfc2Xy5Ml69913NWXKFE2ePFmtWrWqW3DAz1GGAQBoAks37KkxnaAu0pduOueYiUk9NTklrtYxGzdu1IoVKyRJjz32mH766Sc9+eST9coC+CtWkxDbMQMA/NuqVascfw4LC9MDDzxgYBrAu1CGxXbMAAD/9uijj+rjjz/WoEGDNGnSJHXs2NFpTFFRkex2uwHpAGMxTQIAgCYwemBXDevjXELPtLuwrMbUiJzRieoWE17r90SHh9Tp86+44gqtW7dOJ0+edLpmt9t1ww03KCgoSLNnz9bgwYPrdE/AH1CGAQBoAjERofVecqxbTLjiOnjuRTeTyaTmzZs7nX/99dcdG3VcffXV+p//+R+NGjXKY58LeDOmSQAAEOAWLFjg+PN5552n4cOHG5gGaFqUYQAAAtybb76pWbNmqU2bNrJarWrZsqXTGFfTKwB/QBkGACDAhYeHa+rUqbLZbLrjjjucrh8/flzx8fGaOnWqDh8+3PQBgUZEGQYAAJKk6Ohol3OKFy1apG+++UbZ2dkym8365ptvDEgHNA7KMAAAOKuKigrNmjXLcWw2m9WzZ08DEwGeRRkGAABnFRoaqvfee0/XXXedJGnWrFkKCqI+wH/w0wwAAGrVr18/vfPOO9q0aZNSUlKcrleVlejHlzJ1bNdGNu6Az6EMAwCAOklISJDJZHI6f2TDclXs/VoHc636w2+H6fjx4wakAxqGMgwAABps3/ff6ejnqxzHbaKi1KJFCwMTAfXDDnQAAMClwtIKFZdV1jqmLDhCUdfcpSOf/I+qjx1VWsYD2n7gaK3fEx0eUu/d+IDGQhkGAAAuLd2wR/Pzd5xzXOvEmxQRn6Tj332pJz4p1ROfrK9xvfLQHpV+8a7aDByl4PBITUzqqckpcY0VG6gXyrAkq9WqrKwso2MAAHxEYWmFlm7Y4zgePbCrVz/pbIq8QaFhCut5uctrJR8u1bEdn6r0yzVqc3malDTNo59dH435d8HPhW+iDOtUGbZardq6davi4+ONjgMA8HLFZZU1npgO69PRq0uEkXkr9m3TsR2fSpLsJ46r6thPTfK5Z9OYfxf8XPgmyjAAAHBp9MCuGtanY61jdheWKX3pJsdxzuhEdYsJdxzv+baTnjr8G7339hsKC4/Q6n/MU4+u5zVaZqC+KMMAAMClmIjQej8p7BYTrrgOrRzHcR36KuXfr+uzzz7Tjh07dHnv7k7f8+GHH+rgwYP67W9/63LpNqAxsbQaAABodJdddpl+//vfO52vrq7WpEmTlJqaqgEDBmj9+vUuvhtoPJRhAABgmBUrVmjz5s2SpP/93/91/BloKpRhAABgmGbNmqlTp06SpC5dumjs2LEGJ0KgYc4wAAAwzIgRI3TDDTfob3/7m84//3yXu9etW7dOZrNZ559/vgEJ4e94MgwAAAzVsmVLTZkyRaNGjXK6Vl5erttuu009e/bUlClTVFhYaEBC+DPKMAAA8FoLFizQDz/8oIqKCj311FP66KOPjI4EP0MZBgAAXmvv3r2O5dYGDBig3/zmNwYngr+hDAMAAK/1zDPP6PPPP9eNN96omTNnulyHeNOmTTp58qQB6eAPKMMAAMCr9enTR2+99ZauvfZap2v79+/XoEGDdPHFFys3N1d2u92AhPBllGEAAOCzHn/8cR07dkzbtm1TWlqatm7danQk+BjKMAAA8EknTpyosWPdyJEjFR8fb2Ai+CLKMAAA8EnNmzfXli1b9Pzzz6tbt26aPn26y3Hff/99EyeDL2HTDQAAvER0eIgmJvWscYzaNWvWTHfddZfuuOMOBQcHO13//PPPlZiYqN/97nenynJEewNSwpv5xZPh3NxcJSYmKjExUXPmzDE6DgAADRITEarJKXGOr5iIUKMj+QxXRViSHn74YVVXV+vll19WQkKCyspKmzgZvJ3PPxkuKSmRxWLRrl27JEmJiYlKTk5WQkKCwckAAICRvvvuO+Xn5zuOMzIyFB4eYWAieCOvejK8efNmxcbGuryWk5OjqKgoRUVFKSMjw3E+Ly+vRvEdNWqUli1b1uhZAQCAd+vSpYt27typu+++W23btpXFYnEaY7fbVVZWZkA6eAuvKMNz5syRyWRSWlqabDab0/Xc3FxZLBbl5+dr06ZN2rhxo9LS0iRJNptNZrPZMTYyMtLlPQAAQODp3LmzcnJy9O2336pt27ZO19eueVdms1nPPPOMKioqDEgIo3lFGc7MzJTdbld2drbL6xaLRdnZ2UpISJDZbNaKFSuUm5urkpISSXL8X0mKjo5ukswAAMB3tGrVyumc3V6tebOm6+DBg5o4caKuvvpqNu0IQF4/Z7ikpEQ2m03JycmOc2azWZGRkcrLy1NkZKSKi4sd1375pBgAAMCV47ZN+u6/P2/Scfvtt7vc7hn+zevL8OkpD78suGazWTabTenp6bJYLCopKVFkZKSWLVumxYsXGxEVAOAHCksrVFxWWeuY3YVltR67Eh0e0iirQzRWXqnxMjeW+v5dtDD314yFL+qlv2Wr9OhRDR7+O20/cLTG+OrqagUF1e0X6YHwc+FrPxN14fVl+MynvmeKjo5WUVGRIiMjtWLFCiUlJUk69QJdXVaSOHjwoA4dOlTj3M6dO90PDADwaUs37NH8/B31+p70pZvOOWZiUk9NTolraKyzaqy8UuNlbiz1/bswmUx6fk+k7MOekP2nQt34t0+dxvy06S01/36jgi67TaHnXVCvPP74c+FrPxN14fVluC6Sk5O1aVPd/od92sKFC5WVldVIiQAAgK8wmYLUrI3zZhzVlcd05JNlqi4vkb7ZpNYDRijq6jsNSIjG5PVl+GwvxBUXF7t8K7Su7rnnHseKFKft3LlTN998c4PvCQAA/MfRLW+fKsL/L/T8eAPToLF4fRk+PVf4ly/G2Ww2RUZGNvi+7du3V/v2bMkIAKhp9MCuGtanY61jdheW1fiVcs7oRHWLCa/1expra+XGyiv53nbQnv67KD3aV39f1EFLnvubevWO18t/vb/WF+wC4efC134m6sLry3BkZKTMZrNyc3OVmZkp6VQRLikp0ciRIw1OBwDwNzERofV+QahbTLjiOjgv3dUUfC1vY/L430WHVkp4arYem/qAiouLdcGvWjsNefTRR3XixIlTHeUXRZKfC9/gFesMn4vFYpHFYtHmzZtls9mUlpam1NRUt54Mn8lqtcpkMik+nl9/AACAmtq1a6cLLnB+ee7777/XnDlzNHv2bJnNZq1++00D0sFdXlGGc3JyHDvQSafe7jzz1xDp6enKzs5WUlKSYmNjHRtveIrVapXdbldBQYHH7gkAAPzbjBkzHLvWHTlyRF3NsQYnQkN4RRlOT0+X3W53+jpTZmamDh8+LLvd7tEiDAAA0BAZGRkaOnSoJOm2227TBRf2NjgRGsLr5wwDAAB4o8TERL377rtat26dunTpohMuxkydOlWXX365fvOb37C7nZfyiifDRmPOMAAAaKjBgwere/fuTue/+nyzsrOzdcstt2jgwIH6+uuvDUiHc6EMiznDAADA8+bN+nlzr6+++sqt/RHQeCjDAAAAHma325Vyw3B17Hhqrd9JkyapQ4cOBqeCK5RhAAAADzOZTLr1j3dp586devLJJ/Xggw86jbHb7bJardq9e3fTB4QDZRgAAKCRhIWF6YEHHnC5N8Jrr72mrKwsxcXFacKECTpy5IgBCUEZFi/QAQCAplVVVaVp06ZJkk6cOKGVK1cqJMT/tjr2BZRh8QIdAABoWmVlZRowYICCgk5Vsccee0wtW7Y0OFVgogwDAAA0sdatW2vJkiX66quvdO+99+qOO+5wGnP8+HG98MILOnHC1QrG8BTKMAAAgEEuuugiLViwQM2bN3e6tnDhQo0ZM0a9evXSSy+9pOrqagMS+j/KMAAAgJf56aefNHPmTEmSzWbTvHnz2MGukVCGxQt0AADAu+zdu1e/+tWvHMczZ86kDDcSyrB4gQ4AgIaKDg/RxKSejq/ocFZE8ISLLrpIX3zxhf71r3/pj3/8o1JSUpzGHDx4UBs3bjQgnX+hDAMAgAaLiQjV5JQ4x1dMRKjRkfxGcHCwRo8erX/84x8unwo/8cQTuvTSSzVixAh9/fXXBiT0D82MDgAAAID62bNnj5599llJ0sqVK2UymZSbm2twKt/Ek2EAAAAfs2XLFscKFEFBQZo+fbrBiXwXT4YBAAB8zM0336xdu3Zp1qxZOnbsmC666CKnMd9//73CwsLUtm1bAxL6Dp4MAwAA+KAOHTro6aefdkyX+KXx48fLbDbr8ccf19GjR5s4ne+gDIul1QAAgO9y9XLdp59+qjfeeEM//fSTHnnkEc2aNcuAZL6BMiyWVgMAAP5l1apVjj9HRERo8uTJBqbxbpRhAAAAPzN9+nR9+OGHGjRokB544AG1a9fOaczBAz/KbmeLZ16gAwAA8ENXXnml1q1bp6qqKqdr1dXVSv99qvYfPKrIq/6glub+BiT0DjwZBgAA8FMmk0nNmjk/+8zNzdV/C77UiYPf6lBulsr/u96AdN6BMgwAABBgnn/+ecefg1u3V1jcFQamMRZlGAAAIMC89dZbmvb4HAWFRSryyttkatbcaUxZWZkByZoeZRgAACDAhIaGavSYseqU8bzCe1/jdL28vFwXXHCB7rjjDu3evbvpAzYhyrBYZxgAAASmoJAWMgUFO51fsGCB9u3bp3/+85+Ki4vz6+VnKcNiOwKHagAACXlJREFUnWEAAIDTTpw4oXnz5jmOExMT1bt3bwMTNS7KMAAAAByaN2+uTz/9VHfeeaeCgoI0c+ZMl7vcuVqyzRdRhgEAAFBD165d9fe//102m03XXOM8p/iHH35Qjx499Mwzz6iiosKAhJ5DGQYAAIBLXbt2dXn+iSee0O7duzVx4kRdcMEFOnr0aBMn8xzKMAAAAOps7969ysnJcRwPHDhQrVq1MjCReyjDAAAAqLNOnTrpzTffVN++fRUcHKysrCyjI7mFMgwAAIA6M5lMuv7667Vp0yZ9+umniouLcxqzZcsWXXXVVfrwww8NSFg/lGEAAADUW1BQkPr37+/y2sMPP6wPP/xQV111lUaMGCG73d7E6eqOMgwAAACP+eijj/TOO+84jnv27OlyaTZvQRkGAACAx/Tt21dPPPGE2rRpo9atWyszM9PoSLWiDIvtmAEAADwlIiJCDz30kGw2m1577f/au2PdNLo0jOMP2k9K4WYgdZqxu6Qy7hMlUMcFJFcQfAeM3NGO72DsSzDSOjU06Y3bXWk9o12l3C/2NFt4tzhbWBCImQE7wHh4/z8JKSYz8jmvz7x6DOPDX1Wr1YoeUi7CsPg4ZgAAgFWr1Wp6//590cNYiDAMAAAAswjDAAAAMIswDAAAALMIwwAAADCLMAwAAACzCMMAAAAwizAMAAAAswjDAAAAMIswDAAAALMIwwAAADCLMAwAAACzCMMAAAAwizAMAAAAswjDknq9niqVit68eVP0UAAAALBBfxQ9gOeg1+up1+vp6upK9Xpd19fXRQ8JAPCM/fPP/+i///7X5Ot//P1v+t+fOwWOKF/ZxrtO66xF2epctvEua5zj7u7uljq+4pxz6xxQmXz9+lWHh4dFDwMAAAC/6eLiQh8/flx4HGF4Spqm+vbtm169eqUXL16s/ftdX1/r8PBQFxcX2tvbW/v3Kxvqk43a5KM+2ahNNmqTj/pkozbZiqjN3d2dvn//rrdv38rzvIXHc5vEFM/zlvoNYtX29vb0+vXrjX/fsqA+2ahNPuqTjdpkozb5qE82apNt07XZ399f+lj+gA4AAABmEYYBAABgFmEYAAAAZv2l1+v1ih6EZTs7O3r37p12dsq/lck6UJ9s1CYf9clGbbJRm3zUJxu1yfbca8NuEgAAADCL2yQAAABgFmEYAAAAZhGGAQAAYBZheM1OT09VrVZVrVZ1dHS0tnPKJk1THR0dPWqe7XZblUpl5rG7u7uB0W7eU+e67WtnOBw+qMv4kTffbV47V1dXmXOh/2TXhx6UXRv6z/za0H/uLXPtlK73OKzN+fm58zzPjUYjF8ex29/fd61Wa+XnlFGj0XCdTsfFcTyZc6PRyD2n1WpNzpl+bKOnzNXC2rm9vXWj0ejBw/d9NxgMMs/bxrUThqGT5Hzfd/NaufX+s6g+lnvQotpY7j95taH/3Ft07ZSx9xCG18j3fRdF0eTrOI6dJHd7e7vSc8omjmPned7Mc6PRyElyo9Eo87xWq+W63e66h/csPGWuFtbOPGEYuk6nk3vMNq+d8/PzuYGG/nNvXn3oQfey1g79J7s2v7LWf5a5dsrYe7hNYk3SNFWSJGo0GpPnfN+X53kaDocrO6eMfN/X2dnZzHPjzxC/vLwsYkilZ2Xt/CpJEgVBoDAMix7Ks0L/yUcPWi1La2eaxf6z6Nopa+8hDK9JkiSS7n+g03zfn/zfKs4pq1arNff5Wq2We974Pq5qtap2u72OoT0bj5mrpbUzLQgCdTodeZ638FhLa4f+sxg9KB/9ZzGr/Sfv2ilr7yEMr8nNzc3c52u1mn78+LGyc7bF6empPM/LvMjGhsOhoijSaDRSkiRqNpsbGuHmPWauFtdOkiTq9/sKgmCp4y2tHfrP49GDZtF/8tF/fpq+dsraewjDKNz4rabz8/Pc4z5//qzBYKBGoyHf9xWGoYbD4Va+8mBprk8VRZF833/wasI81BN56EGzrMzzd9B/7i177Tx3hOE1yXqr7ebmRi9fvlzZOWWXpqmazabOzs5m7heap9VqzRxzcHAgSVt5T9pj52px7fT7/YWv4o1ZWjsS/ecx6EEP0X8Wo//Mv3bK2nsIw2sy/m3x19/8kiTJvL/oKeeUWZqmqtfriqJo6aYybfzWyqJ7/LbBorlaXDu/81bjtq8d+s9y6EHLof/Mov9kXztl7T2E4TXxPE++76vf70+eS5JEaZrq06dPKzunrNI01YcPHxRF0cJXY8am6yLd/0GC9PMvWbfJY+dqae1IP//if5m3KCVba0ei/yyDHpSN/pPPev/Ju3ZK23s2soGbUVEUTfbem7eJdBzHM/vqLXPOttjf33fdbvfB5uXjPQV/rc14b8MoiiYbn3uet3B/xzJaZq6W145zP/cAnbcHpbW1k7UfKv3nXlZ96EH5ezBb7z95+wxb7z+Lrp0y9h7C8JqFYeg8z3OSHvxgxz/8x5yzDcabac97hGHonJtfmziOXaPRmHw60PjYbbRorlbXztj4U6LmsbJ2xvP89THNcv/Jq4/1HrRo7VjuP8teV1b7zzLXjnPl6z0V55x70kvKAAAAQMlxzzAAAADMIgwDAADALMIwAAAAzCIMAwAAwCzCMAAAAMwiDAMAAMAswjAAAADMIgwDAADALMIwAAAAzCIMA8CWOD09VaVSWfhoNpuTc/r9vur1uiqViqrVqtrttq6urgqcBQBsFh/HDABbJE3Tyb+TJFG9XlcYhup0OjPHeZ6nk5MTBUGgKIp0cHCgJEk0GAwkSVEUbXTcAFAUwjAAbKk0TVWtVhVF0YMwLEnVanVuUAYAS7hNAgCMStNUtVqt6GEAQKEIwwBgVKvV0pcvXzQcDoseCgAUhjAMAEadnZ3p4OBAzWZTlUpF7XabYAzAHMIwABjleZ4Gg4HiOFYYhkqSRM1mUycnJ0UPDQA2hjAMAMb5vq9ut6vRaKRut6sgCIoeEgBsDGEYADBxfHwsSew1DMAMwjAAGJUkyYPnLi8vJd2/WgwAFvxR9AAAAJuXJIl2d3fV6XTUbDbl+74uLy8VBIE6nY48zyt6iACwEbwyDAAG+b6vOI51c3OjIAgmn1R3fHzMp88BMIVPoAMAAIBZvDIMAAAAswjDAAAAMIswDAAAALMIwwAAADCLMAwAAACzCMMAAAAwizAMAAAAswjDAAAAMIswDAAAALMIwwAAADCLMAwAAACzCMMAAAAwizAMAAAAswjDAAAAMOv/8t/k3KY5kgcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "h = sum_bg.get_hist(bins=25)\n", "hl.plot1d(ax, h, crosses=True)\n", "norm = h.integrate().values\n", "x = np.linspace(.5, 20, 300)\n", "ax.plot(x, norm * sum_bg.pdf(x), 'k:')\n", "ax.set_xlabel(r'TS')\n", "ax.set_ylabel(r'trials per bin')\n", "ax.semilogy()\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Hottest source analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We might be interested in the hottest source in each trial. Then the \"post-trials\" test statistic is just the maximum `mlog10p` from each trial. This can be obtained like so:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "mlog10p = np.array([bg_trials[name] for name in bg_trials.keys() if name.startswith('mlog10p_')])\n", "ts = np.array([bg_trials[name] for name in bg_trials.keys() if name.startswith('ts_')])\n", "ns = np.array([bg_trials[name] for name in bg_trials.keys() if name.startswith('ns_')])\n", "gamma = np.array([bg_trials[name] for name in bg_trials.keys() if name.startswith('gamma_')])" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "imax = np.argmax(mlog10p, axis=0)\n", "r = np.arange(len(imax))\n", "bg_trials_hottest = cy.utils.Arrays(dict(\n", " mlog10p=mlog10p[imax,r], ts=ts[imax,r], ns=ns[imax,r], gamma=gamma[imax,r]))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7MAAALDCAYAAADHS1caAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAASdAAAEnQB3mYfeAAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdQWwb953//Q/7FDZQuc2IFly4QOBoZOdiwUkpB4UPD7qoqO3BKOo2HPnkU2syycUI0OUkexpdNqF6KHJJVnRv7mFtskEWPe2KWiDP/+AtVuKmf0j4P0iksd38UaCJTE0bq4D8/It5DioZUaIkihxqZsj3CxASDjnDLylTX35nfr/vL+H7vi8AAAAAAGLkK2EHAAAAAADAUVHMAgAAAABih2IWAAAAABA7FLMAAAAAgNihmAUAAAAAxA7FLAAAAAAgdihmAQAAAACxQzELAAAAAIgdilkAAAAAQOx8NewAosTzPH344Yd69tlndfLkybDDAQBE2NbWlj799FN997vflWEYYYcTWeRWAEC7jppbKWZ3+PDDD3Xt2rWwwwAAxMgHH3ygH/7wh2GHEVnkVgDAUbWbWylmd3j22Wclbb9558+fDzkaAECUra6u6tq1a43cgdbIrQCAdh01t1LMSnIcRzMzM43b58+f18WLF0OMCAAQFwydPVj9/SG3AgDa1W5upQGUtotZ3/e1vLwcdigAAPQFx3GUSCQ0Pj4edigAgD5FMQsAAALHiWIAQK9RzAIAAAAAYodiFgAAAAAQOxSzAAAgcMyZBQD0GsWsSLgAAASNObMAgF6jmBUJFwAAAADihmIWAAAAABA7Xw07gH61/mRLtc2nR94vOXRCI6faWyQYAIBBQm4FAOxEMdsjd+4/0jsLnxx5v1uTF/T61PM9iAgAgOPjOI5mZmYCPSa5FQCwE8OMAQBA4OhHAQDoNYpZAAAAAEDsMMxYvRkKdePKOV29dPbI+yWHTgQaBwAA/YLcCgDYiWJW28Ws4zhaWVkJbK3ZkVMnaTYBAECAyK0AgJ0YZgwAAALnOI4SiURgJ4kBANiNYhYAAASOBlAAgF6jmAUAAAAAxA7FLAAAAAAgdihmAQAAAACxQzELAAAAAIgdilkAAAAAQOxQzIrlAwAACBq5FQDQaxSzYvkAAACCRm4FAPQaxSwAAAAAIHYoZgEAAAAAsUMxCwAAAACIHYpZAAAAAEDsUMwCAAAAAGKHYhYAAAAAEDsUswAAAACA2KGYBQAAgXMcR4lEQuPj42GHAgDoUxSzAAAgcI7jyPd9LS8vhx0KAKBPUcwCAAAAAGKHYlYMhQIAAACAuKGYFUOhAAAAACBuKGYBAAAAALFDMQsAAAAAiB2KWQAAAABA7FDMAgAAAABih2IWAAAAABA7FLMAAAAAgNihmAUAAAAAxM5Xww4AzdafbKm2+bSjfZNDJzRy6mTAEQEAcHSO42hmZibsMAAAfYxiNmLu3H+kdxY+6WjfW5MX9PrU8wFHBADA0TmOI8dxtLKyovHx8bDDAQD0IYYZAwAAAABih2IWAAAAABA7DDOOmBtXzunqpbMd7ZscOhFwNAAAAAAQTX1RzM7Ozmp+fl61Wk3pdFqFQiHskDo2cuokTZwAAAAA4BCxL2Y9z5Mkzc/PS5KGh4d1/fp1pVKpMMMCAAAAAPRQpIrZarUqy7K0tra2575isSjbtiVJ09PTmpubkyQZhqF8Pi9Jcl1XnufJNM3jCzpCOl3WhyV9AAAAAMRNJIrZ2dlZ2bYt0zTluu6e+8vlsmzb1sLCggzDkGVZsixLpVKp8RjLslQul1UqlWQYxnGGHxmdLuvDkj4AgH7GGu4A0J8iUczm83nl83mVy2VZlrXnftu2VSgUGkOHS6WSxsbG5Hleo3C9ffu23nzzTVmWJdM0GWYMAAAksYY7APSryC/N43meXNdVOp1ubDNNU4ZhqFKpNLYZhqFUKqVUKtUYggwAAAAA6E+RuDJ7kPqw493zYOtDksvlctOV2Gq12phbe5DPPvtMn3/+edO21dXVgKIOR6fL+rCkDwAAAIC4iXwxW6vVWm5PJpN6/Pix8vm8bNuWbduq1WrKZDLKZrOHHvfdd9/VzMxM0OGGimV9AADYizXcAaA/Rb6YbUcn68q+9tpre+bnrq6u6tq1a0GFBQBAbLGGOwAg6iJfzCaTyZbba7WaTp8+3fFxz5w5ozNnznS8fz+hyyMAYCfWcAcAxEHki9n6XFnXdZvmzbquO7BL8ASNLo8A0L9Ywx0A0K8i383YMAyZpqlyudzYVk+s09PTgTyH4zhKJBIaHx8P5HgAAIRtdnZWiURClmUduob70tKSFhcX90y/sSxLY2NjA72GOwAguiJfzEpqNHiqVqtyXVeWZSmTyQSWWB3Hke/7Wl5eDuR4AACELZ/Py/f9fee67lzD3TRNlUollcvlxhBjaXsN96WlpUYOBgAgSiIxzLhYLCqXyzVuJxIJSZLv+5KkbDYrz/M0OTkpz/OUyWRUKpVCibUf0eURAAbLYWu4ZzIZSXvXcGcddwBAlESimM1ms4cup5PP5xvzd4LmOE7fLdNzFHR5BIDBwhruAIB+EIliNmyO48hxHK2srDBvFgDQ91jDvX2ddvyn2z8A9B7FLAAA2IM13Ld12vGfbv8A0HsUs+gKZ6wBIH5Ywx0A0A8oZsWc2W5wxhoA4uc41nAntwIAeo1iVsyZBQAMlp1ruNebK/ZiDfd+yK2ddvyn2z8A9B7FLAAAA8i2beVyOaXTaRmGEfga7v2Cjv8AEF0Us+gKZ6wBIJpYwx0A0O8oZtEVzlgDQDSxhjsAoN99JewAosBxHCUSiVjP6QEAIEocx5Hv+1peXg47FABAn6KYFQkXAAAAAOKGYhYAAAAAEDsUswAAIHBM4QEA9BrFLAAACBxTeAAAvUYxCwAAAACIHZbmEcsHAACAYK0/2VJt82lH+yaHTrDsHQC0gWJW28Ws4zhaWVlhbs8xaZXkN/7yVL/53R8at3/wwrc0/LUTe/YlyQNA9A36ieI79x/pnYVPOtr31uQFvT71fMARAUD/oZhFKNpJ8r/6z9+33E6SB4Do40QxAKDXmDMLAAAAAIgdrswCAAAE7MaVc7p66WxH+yaH9k6xAQDsRTGLUHST5CXp4z9+ceR9mGsLADguI6dOknMAoMcoZhGKbpL8L+Y/7qipBnNtAQAAgP7BnFltN6lIJBI0qAAAICDkVgBAr1HMajvh+r6v5eXlsEMBAKAvkFs7t/5kSx//8Ysj/6w/2Qo7dAA4VgwzRuy0mm/bzhq1NNQAAMRBp2vUMp0GwKChmEXs7Dff9jujp0OIBgAAAEAYGGYMAAAAAIgdrswCAABECNNpAKA9FLMAACBwjuNoZmYm7DBiiek0ANAehhkDAIDA0c0YANBrFLMAAAAAgNhhmLEYCjUo1p9sqbb5tGlbu3OQWg33AgAAABAeilltF7OO42hlZUXj4+Nhh4MeaWfdvl/95+/3bGPdPgAAACB6GGYMAAAAAIgdilkAAAAAQOwENsz4o48+arn9xRdfDOopgK6wbh+AuCG3AgCwv66L2du3b+uVV16RJPm+33RfIpHQX//6126fAggE6/YBiAtyKzrRqtFhu2h2CCCOui5mbdvWzZs3Zdu2kslkEDEBADDQ+iG3slLA8Wun0eF+aHYIII4CmTP7xhtvaHR0VM8888yeHwAAcHRxz62O48j3fS0vL4cdCgCgT3VdzGazWf36178OIhYAACByKwAA7QikAdQ//MM/6J/+6Z90+fJlGYbRdN/du3eDeAoAAAYKuRVH1WmjQ4lmhwDiqeti1nVdZTKZxu3djSoAAMDRkFvRCRodAhg0XRez9+7dCyIOAADwN+RWAAAOF0gDqLhzHEeJRELj4+NhhwIAAAAAaMORr8x+//vfl2VZ+ulPfypJevXVVw98/HvvvddZZMfIcRw5jqOVlRUKWgDAsevH3AoAQK8duZh9/PjxgbcBAMDRkFsBADi6Ixezi4uLTbeZ14NBtf5kS3fuP2rcvnHlXMvGGwBwGHIrAABHF8jSPEA/W3+ypdrm0z3bH65v6p2FTxq3L37rG3puZKjpMcmhExS4AAAAQA8EUsy+//77mpuba5xZvnz5smZnZ/XCCy8EcXggVHfuP2oqWveTvbO0Z9utyQt6fer5XoQFoM+RWwEAOFjXxezPf/5z2batTCajt99+W9L2cKlUKqVyuawf/ehHXQcJAMAg6Yfc6jiOZmZmwg4DbWo1CmnjL0/1m9/9oXH7By98S8NfO9H0GEYgAQhT18XsW2+9pdnZWf3sZz9rbLt586ampqaUz+djkXABAIiSfsitrBQQL+2MQvrVf/5+zzZGIAEIUyDDjDOZTMttb7zxRhCHB0J148o5Xb10ds/2ds9YA0AnyK0AABys62L2zTffVKVSaayNV/fw4UNNTEx0e3ggdCOnTu47hOo7o6ePORoAg4DcCgDA4Y5czO5eyP3x48f69a9/rfn5+abtlUpFpml2Fx0AAAOA3IqwtRqF9HB9s6m5YfHGRMuu/QAQliMXs60Wcn/55Zfl+37TtsnJyc6jAgBggJBbEbZWo5CSQyd0a/JC43bq3DDNngBEypGLWRZyB9pHd0gA7SC3IopGTp2kuROASAukARSA1ugOCQAAAPTGV8IOIAi2bWtiYkJjY2OanZ0NOxwAAAAAQI/F/spsuVyW53laWtpuUDA2NqZUKqV0Oh1yZAAAAACAXolUMVutVmVZltbW1vbcVywWZdu2JGl6elpzc3OSpHQ63VS4ptNpVatVillEAt0hAQAAgN6IRDE7Ozsr27ZlmqZc191zf7lclm3bWlhYkGEYsixLlmWpVCrJMIymx1YqFeVyueMKHTgQ3SEBAACA3ohEMZvP55XP51Uul2VZ1p77bdtWoVBQKpWSJJVKJY2NjcnzvKZi1rZtZTKZxuOAKKI7JAAAANC9rhtA/fKXv9Sf//znIGJpyfM8ua7bNGzYNE0ZhqFKpdLYZtu2xsbGVCgUehYLAADHode5tR00VwQARF3XV2bffvttJRIJ/eQnPwkinj3qw45N02zavnNIci6Xk2VZR5on+9lnn+nzzz9v2ra6utpltAAAdK/XufUwNFdEN9afbOnO/UeN2zeunGM6DYCeCKSYzeVySqfTOnfuXBAxNanVai23J5NJPX78WOVyWcVisekqbS6XUz6fP/C47777rmZmZgKNFeg1viAAgyHI3EpzRfTS+pMt1TafNm17uL7ZtMb6xW99o2WjQ/IXgG51XcwuLi5qdHRUpmkqnU7vuYL63nvvdfsUB8pkMvJ9/8j7vfbaa3vm566ururatWtBhQZ0rNWXA4kvCMCgCCK30lwRx+HO/UdNeamVnR38625NXqB/BICudV3Muq4r0zQbifbx48ddB7VTMplsub1Wq+n06dMdH/fMmTM6c+ZMx/sDvdTOlwOJLwhAvwoit9JcEQDQ77ouZu/duxdEHPuqJ/J6Yq9zXXfPmWMAAPpBr3PrYc0VM5mMpC+bK2az2Z7GAwBAJwJZmuf999/Xv/zLv+ijjz7Sxx9/LEl69dVXZVmWvve973V1bMMwZJqmyuVyYx6s67ryPE/T09Ndxy5JjuMwfxYAECm9zK29aq6IwXPjyjldvXS2advD9c2mkUPFGxMtp8QAQLe6LmZ//vOfq1gs6p//+Z/193//943tk5OTKhQKXSdcafvMcL0RRn1eTyaTCezKrOM4chxHKysrGh8fD+SYQDdafTmQ+IIADIpe59ZeNVdkpYDBM3Lq5J4+DcmhE7o1eaFxO3VumF4OAHoikG7GCwsLevHFF5u2T0xM6Pr1620do1gsNjWWSCQSktRo7JTNZuV5niYnJ+V5njKZjEqlUrehA5HV6suBxBcEYFAEkVu70WlzRVYKgLSdw+jdAOA4dF3M+r6v4eHhPdtd19Xo6Ghbx8hms4fOx6k3sugFhhkjLviCAAyGIHLrQXrVXJGVAtAtlqADcBRdF7M3b95UNpttulL64MEDvfLKK3rllVe6PfyxYJgx+h1fDoB46XVu7VVzRVYKQLtYgg5AELouZguFgnK5XCP5XbhwQa7rKpvN6mc/+1nXAQI4GhawB+Kv17mV5ooIG0vQAQhCIN2M5+bmZNu2/vu//1uSlEqlAhkGBeDoWMAe6A+9zq00VwQAxF0gxaykpsXd44azxwCAKOomt9JcEYOE6TTAYApsndm5uTktLi5Kki5fvqzZ2Vm98MILQRy+5zh7DACImm5zK80VEWVBL0FX23zaNCrp6qWzFLPAAOi6mH3jjTf085//XJOTk3r77bclSffu3VMqlVK5XNaPfvSjroME0D4WsAfirx9yKyeKcZBulqBbf7Klj//4RdO2h+ubB96uH5sCF+gvXRezxWJR+Xxeb731VmPbzZs3Zdu28vl8LBIu0E/2+4Kw03MjQ3r+m18/9FgM2wLCQW7FoGpnCTp6QwCo67qYTSaTTXNy6v7xH/9Rt2/f7vbwAAKw+0x3u1dhGbYFhIPcCgDA4b7S7QEymYwWFhb2bF9YWND169e7PfyxcBxHiUSCYVDoW/Uz3fWfVgVpfdjWzp9Ww7Z2P+bjP36h9Sdbx/VSgIFAbgUA4HAJv97WsEPT09P69a9/rUwm09jmeZ4qlUqj3X/d3bt3u3mqnqvP61leXtbFixfDDgc4Vr+Y/7itNf9aYegWBlEvcwa5Fdjffuupt9MbgtFFQLQdNWcE0s345Zdf1s6a+JlnntHLL78sSeqyVgYAYCCRW4HWWvWGaKdxFID+03Uxe+/evSDiAAAAf0NuBY6mncZRAPpPIFdmAcRfp0v6SCzrAwAAgONHMQtAUrBL+gCA4ziamZkJOwwAQB/ruptxP6DjIgAAwXIcR77va3l5OexQAAB9imJWJFxgP/WGGvUfhhMDAAAgKhhmDGBfNNQAAABAVHV9Zfb999/X+++/37j96quv6vTp03rppZf08OHDbg8PAMDAIbcCx2f9yZZ+Mf9x42f9yVbYIQFoU9fFrG3bjcXbb9++rWKxqGKxqFQqJcuyug4QAIBB0w+5lX4UiIva5lO9s/BJ46e2+TTskAC0qethxmtra7p8+bIkqVQqKZvN6uWXX9a3v/1tXbhw4ZC9AQDAbv2QWx3HkeM4WllZoaBFJKw/2WpZqD5c3zzwtrTdQ+Kwjv8Ajl/XxaxpmtrY2NA3vvENVSoVvfLKK5KkP/3pT42zygAAoH3kViB4d+4/0jsLnxz6uJ3rq9fdmrxADwkggroeZrzzTLFpmvrxj38sSapUKpqcnOw6wOPAUCgAQJT0Q24FAKDXur4yWygUdP36dT148EDpdLqx3TTNpttRxlAoAECU9ENuBQCg1wJZmieVSimVSjVte/nll4M4NAAAA4ncCgTrxpVzunrp7J7tD9c3m4YWF29M6LmRoabHsM46EE1HLmZfffXVIz3+vffeO+pTAAAwUMitQO+NnDrZVhOn50aG9Pw3v34MEQHo1pGL2cePH/ciDgAABlY/5lbHcTQzMxN2GACAPnbkYvbevXu9iAMAgIHVj7mVfhSIi+TQCd2avNB0G0A8BDJnFgAAAIijkVMnO1p2Z/3Jlu7cf9S4fePKOdaiBY5ZIMXsRx99pLt378p13T333b17N4inABAzJHmgO+RWINpqm0+b1q29eukseQ44Zl2vM3v79m2lUinNz8/L932Vy2X5vq/5+Xl5nhdEjABiqJ7k6z+1zadhhwTEBrkViJb1J1v6+I9fNP08XN9seszD9c09j1l/shVSxMBg6PrK7BtvvKH5+fnGIu7JZFL37t1TpVLR7du3uw7wONCkAgAQJf2QW4F+cuf+o6arsK3sXN6n7tbkhY6GMANoT9dXZjc2NvTSSy81bpumqS+++EIvvfSSKpVKt4c/Fo7jyPd9LS8vhx0KAAB9kVsBAOi1rq/MplIpVSoV/fjHP5YkXb58Wffu3dPw8DBDoYABsf5ka88w4lbDr3ZLDp3YM7+om7m2zNNFvyC3AgBwuK6L2UKhoPn5+UbCtW1bY2NjSiQSymazXQcIIPqCHH7VTUMNmnGgX/RDbmUKD/rJjSvndPXS2aZtD9c3m3Jb8caEnhsZanpMq2V+OPEKBKfrYnZycrIxp0eSRkdHtbGxoVqtptHR0W4PDwDAwOmH3Mo6s+gnI6dOHlpwPjcypOe/+fVDj8WJVyA4PVln9plnntEzzzzTi0MDADCQyK0AADQ7cjH7/e9/X5Zl6ac//akk6dVXXz3w8e+9915nkQGIjU6HX0nSx3/8Ys9+B90+SKfzdIGwkVsBADi6Ixezjx8/PvA2gMHT6fCrX8x/3NFcW0n6zmhSv31QO/K+LJOAKCK3AvGTHDqhW5MXmm7v1mmDxPrxOPkKHOzIxezi4mLT7Xv37gUWDAAAg4jcCsTPyKmTh54c7bRBosTJV6AdXa8z+8tf/lJ//vOfg4gFAACI3AoAQDu6bgD19ttvK5FI6Cc/+UkQ8QAYIN3MtW2l02USgKghtwIAcLhAitlcLqd0Oq1z584FEROAPtDOXKIglzoIel8gTORWoD90c9KWk6/A4bouZhcXFzU6OirTNJVOp2WaZtP9cei4yMLuQPDamUsEoLV+yK0Aen/SFhh0XRezruvKNM1Goo1jB0YWdgcAREk/5FYAAHqt62KWjosAAASL3AogSOtPtnTn/qPG7RtXzrHsD/pC18Xsf/zHf+h73/venu0fffSRJOnFF1/s9ikADJB25tr2Yl8gSvohtzKFB2gtjFxV23zatETQ1UtnKWbRF7ouZqempvTXv/51z/bHjx9rdnZW//Zv/9btUwAYIN3MtWWeLvpFP+RWpvAArZGrgOB0Xcz6vr/vfbsXgQcAAIcjtwJoheHCQLOOi9lkMqlEIqFEIqHTp0833ed5nnzf18TERNcBAkCvdfPlgC8WCBK5FcBB2hkuvP5kS7XNp03bHq5vHni7Ljl0ghyGWOm4mF1YWJDv+7p8+XLLRhWmaWp0dLSr4ADgOHQzl4h5SAgSuRVAt+7cf9SUl1rZuc7tTrcmLzAEGrHScTH77W9/W5KUTqc1OTkZWEAAAAwqcisAAO37SrcH+Pd///cg4gAAAH9DbgUA4HBdN4ACABwdc20BAAfpdO7r1UtndfXS2T2P2zm0uHhjQs+NDO3ZlyXtEDcUswAGSjeNMVppZ99WDTWYawsAOEinc1/bmff63MiQnv/m17uKD4gCilkAA6XTLwffGU3qtw9qhx6/0y8WAAAAOJqu58wCAAAAAHDc+ubKrOu6KhQKmpubCzsUAAAAoCs3rpzraO4r814xSPqimJ2amlKtVlMymQw7FAAR1+mXg/3wxQIA0Asjp04e2kuBua8YdJEqZqvVqizL0tra2p77isWibNuWJE1PTzddgZ2fn1e5XOaqLIBD9frLQat9159s6eM/ftG0rdPGUQAAHFVy6IRuTV5out2uMLrv0/Ef7YpEMTs7OyvbtmWaplzX3XN/uVyWbdtaWFiQYRiyLEuWZalUKoUQLQAcTS87UgK9xBQeoD+MnDrZcT4Jo/s+Hf/Rrkg0gMrn8/J9X4VCoeX9tm2rUCgolUrJNE2VSiWVy2V5nnfMkQIAMBimpqZkWVbLk8wAAERBJIrZg3ieJ9d1lU6nG9tM05RhGKpUKiFGBgBA9FWrVY2NjbW8r1gsanh4WMPDw8rlck33zc/P68033zyOEAEA6EgkhhkfpH5G2DTNpu37DUlu12effabPP/+8advq6mrHxwOA/dCREmFgCg8AoN9Fvpit1WottyeTST1+/FiSlMvltLi4KNd1NTU11RiSfJB3331XMzMzgccLALvRkRJhyOfzyufzKpfLsixrz/07p/BIUqlU0tjYmDzPk2EYxx0ugDZ008gJ6EeRL2bb0Uljitdee21Pcl9dXdW1a9eCCgsAgEg6bApPJpMJMToA++mmkVM71p9sqbb5tGlbL7vvt3q+Xj8n+kvki9n91o6t1Wo6ffp0x8c9c+aMzpw50/H+APpHN2e6OUuOOGIKD4BWjrv7fjvPF/Rzor9EvpitJ1rXdZuSruu6DIMCEIhuznT3+iw50AtM4QEA9IPIF7OGYcg0TZXLZeXzeUnbhazneZqeng7kORzHIfkCALADU3gAAFEX+WJW2m5SkcvllE6nGx0XM5lMYFdmHceR4zhaWVnR+Ph4IMcEACCqmMIDoJXj7r7f6vl6/ZzoL5EoZovFYtP6dolEQpLk+74kKZvNyvM8TU5OyvM8ZTIZlg4AAKBDxzGFh1FPQPwcd/f9dp4v6OdEf/lK2AFI28Wq7/t7fnbK5/Pa2NiQ7/uBF7KO4yiRSHBVFgAwEHZO4anrxRQe3/e1vLwcyPEAANgtEsVs2Ei4AIBBY9u2bNtWtVqV67qBT+EBAKDXIjHMGAAGDUv6oNeYwgMA6HdcmQWAENSX9Kn/sPA7gsYUHgBAv6OYFQkXAICgMYUHANBrFLMi4QIAAABA3FDMAgAAAABihwZQAAAgcKwzC/SHMBoW0iQR7aKYBQAAgXMcR47jaGVlhZ4UQIzVGxb2+3MinhhmLBpAAQAAAEDcUMyKBlAAAAAAEDcMMwYAAIFjziww2NafbOnO/UeN2zeunGNNdQSOK7MAACBwjHoCBltt86neWfik8VPbfBp2SOhDFLMAAAAAgNihmAUAAAAAxA7FrOhmDAAAAABxQwMosRYegHjptKlGnJpxxClWAAAQDopZAIiZelONuquXzrZV6HW6XxjiFCtao5sxAKDXGGYMAAACRzdjAECvUcwCAAAAAGKHYhYAAAAAEDvMmQUAAADQsfUnW6ptPm3a9nB988DbkpQcOkE/BHSFYlY0qQAAAAA6def+o6amfa1k7yzt2XZr8oJen3q+V2FhADDMWDSpAAAgaKzhDgDoNYpZAAAQOE4UAwB6jWHGAAAAADp248o5Xb10tmnbw/XNpqHFxRsTem5kqOkxyaETxxIf+hfFLAAAAICOjZw6eWgjp+dGhvT8N79+TBFhUDDMGAAAAAAQO1yZBQAAABB760+2dOf+o8btG1fO9eXSP4PyOttBMQsAEdVq3T6pvbX7WonqmikAUNgAACAASURBVH+drk8osUYhAOBLtc2nTUsEXb10ti9zxKC8znZQzAJARLWzbp/Ueu2+74wm9dsHtSPvF8aaf52uTyixRiEAAIOMObNiLTwAAAAAiBuKWbEWHgAAQeNEMQCg1xhmDAAR1WrdPqm9tftaieqaf52uTyixRmGUOY4jx3G0srJCQQsA6AmKWQCIqHbW7ZM6X7svKmv+sT4hAADoBMOMAQAAAACxQzELAAAAAIgdilkAAAAAQOxQzAIAAAAAYodiFgAAAAAQOxSzAAAAAIDYoZgFAAAAAMQO68wCAAAAiJX1J1uqbT5t2vZwffPA25KUHDrR1hruUTEor7NTFLOSHMfRzMxM2GEAAAKw/mRLd+4/aty+ceXcQCR0ABgkd+4/0jsLnxz4mOydpT3bbk1e0OtTz/cqrMANyuvsFMWstotZx3G0srKi8fHxsMMBAHShtvm0KfFfvXSWYjYEnCgGBlty6IRuTV5ouh1lnZ4IjdMJ1G5ijerrpJgFAACB40QxMNhGTp2M1ZXBTk+ExukEajexRvV1UswCAAAAiJUbV87p6qWzTdserm82Dbkt3pjQcyNDTY+J+hXi3QbldXaKYhYAAABArIycOnnolcHnRob0/De/fkwR9cagvM5OsTQPAAAAACB2KGYBAAAAALFDMQsAAAAAiB2KWQAAAABA7FDMAgAAAABih2IWAAAAABA7FLMAAAAAgNihmAUAAAAAxA7FLAAAAAAgdr4adgBBKJfLeuuttyRJ169fVz6fDzkiAOid5NAJ3Zq80HS7l/uFIU6x9ityKwAg6mJfzHqeJ9u2tba2JkmamJhQOp1WKpUKOTIA6I2RUyf1+tTzx7ZfGOIUaz8itwIA4iBSw4yr1arGxsZa3lcsFjU8PKzh4WHlcrnG9kql0pRcr1+/rrt37/Y8VgAA4oDcCgDoV5EoZmdnZ5VIJGRZllzX3XN/uVyWbdtaWFjQ0tKSFhcXZVmWJMl1XZmm2XisYRgtjwEAwCAhtwIA+l0kitl8Pi/f91UoFFreb9u2CoWCUqmUTNNUqVRSuVyW53mS1PivJCWTyWOJGQCAKCO3AgD6XeTnzHqeJ9d1lU6nG9tM05RhGKpUKjIMQ7VarXHf7rPJAID+tf5kS7XNp03bHq5vHnhb2m4oNXLqZE9jizJyK4BB1SpvSO3ljlaimnM6zY/7ierrjHwxWx/WtDuJmqYp13WVzWZl27Y8z5NhGLp7965u37596HE/++wzff75503bVldXgwscANBzd+4/0jsLnxz4mOydpT3bbk1eGOgGU+RWAIOqnbwhtc4d3xlN6rcPai0effB+YeScTvOjFK/XGflidueZ4Z2SyaQeP34swzBUKpU0OTkpabtJRTvdFt99913NzMwEGisAAHFAbgUA9IPIF7PtSKfTWlpqfWZhP6+99lqj0UXd6uqqrl27FmRoAADEErkVABB1kS9m92s6UavVdPr06Y6Pe+bMGZ05c6bj/QEA4btx5ZyuXjrbtO3h+mbT8KfijQk9NzLU9Jjk0IljiS+qyK0ABlWrvCG1lztaiWrO6TQ/7ieqrzPyxWx9Ps/u5hOu68owjLDCAgBEwMipk4c2m3huZEjPf/PrxxRRPJBbAQyqdvKG1HnuiErO6XV+jMrrjMTSPAcxDEOmaapcLje2ua4rz/M0PT0dyHM4jqNEIqHx8fFAjgcAQJSRWwEA/SDyxay0vRaebduqVqtyXVeWZSmTyQR29thxHPm+r+Xl5UCOBwBA1JFbAQBxF4lhxsViUblcrnE7kUhIknzflyRls1l5nqfJyUl5nqdMJqNSqRRKrAAAxAG5FQDQ7yJxZTabzcr3/T0/O+XzeW1sbMj3/cCTLUOhAAD9htwKAOh3kShmw8ZQKAAAgkVuBQD0GsUsAAAAACB2KGYBAAAAALFDMSvm9QAAEDRyKwCg1yLRzThsjuPIcRxVq1VNTExodXU17JAAAB16uL6pp58/atz+5P/9X/r/1ocCf556rtja2gr82P2A3ArguHXz97/TfY8r5wTxnGG8P0d11Nya8He3Nhxg//qv/6pr166FHQYAIEY++OAD/fCHPww7jMgitwIAjqrd3Eoxu4Pnefrwww/17LPP6uTJk2GHE2urq6u6du2aPvjgA50/fz7scAYC7/nx4v0+flF7z7e2tvTpp5/qu9/9rgzDCDucyAoqt0bt94+D8fuKH35n8dKvv6+j5laGGe9gGAZn1wN2/vx5Xbx4MewwBgrv+fHi/T5+UXrPU6lU2CFEXtC5NUq/fxyO31f88DuLl378fR0lt9IACgAAAAAQOxSzAAAAAIDYoZgFAAAAAMTO/+U4jhN2EOhPQ0ND+ru/+zsNDfW2PTm+xHt+vHi/jx/v+WDj9x8v/L7ih99ZvPD7opsxAAAAACCGGGYMAAAAAIgdilkAAAAAQOxQzAIAAAAAYodiFgAAAAAQOxSzAAAAAIDYoZhF4KrVqsbGxsIOYyB4nqdcLqfh4WENDw8rl8uFHdJA2PmeW5Ylz/PCDmlguK6rYrHIez5gisUif+diwrIsJRKJph++E0TLQd/T+KxF036/Mz5vFLMI0OzsrBKJhCzLkuu6YYczECzLkiQtLS3p9u3bunfvnqampkKOqr9NTExIkhYWFrS0tCTP8zQ5ORlyVIPDsizlcjnVarWwQ8ExKZfLsm278ZlbXFxs/O1DNGWzWa2trTV+5ufnww4JOvx7Gp+16Gnnu/Wgf95YZxaBK5fLsixL/NPqLdd1NTExoY2Njca2arWqiYkJLS0tKZVKhRhd/yoWi8pms43blUpFU1NT2tjYkGEYIUbW/4rFoubn51Uul7W2tibTNMMOCcdgbGxMtm03Pneu62psbIzPXERZliXTNFUoFMIOBfvY73san7Xo2u93xueNK7NAbJmmqdu3bzdtqxewi4uLYYQ0EHYWspI0Pz8vwzBI9MfAtm29+eabYYeBY+R5nlzXVTqdbmwzTVOGYahSqYQYGdBf+KwhrihmgRjLZDIttyeTyWOOZDBVKhXNzs5SYB0D27Z1+fJlRhwMmPqwut1X4U3TZDpLhNXn99X7CiD6+KzF16B/3ihmgT5SLBZlGMa+RS6CUa1WlUgkNDU1pXQ6rXw+H3ZIfc3zPM3Ozg70MKpBtd/c6GQyqcePHx9zNGhXpVLR3NyclpaW5LouvRxigM9afA36541iFugTruvKtm2VSqWwQ+l7qVRKGxsbWlpakqSBSxzH7ebNm8pkMlyVBWLg+vXrmp+fVzqdbszlq1QqXN0DeoDPG8Us0Bc8z9PU1JRu377dNN8FvWMYhlKplEqlkiqVisrlctgh9aVqtapyubxnfjgGw35TJmq1mk6fPn3M0aAdmUymKQ9dvnxZkph3GXF81uKJzxvFLBB7nudpYmJCc3NzDC8OQb3503/913+FHUpfunv3riRpdHS0sfahtN11s75MEvpXff7e7qsMruvSdC0m6sNX6eUQbXzW+sMgft4oZoEYq69xOjc3xxXZY+C67p5E73mePM8buEXKj0uhUNDGxoYePHigBw8eaGFhQdJ2F+n6/6N/GYYh0zSbRj64rivP8zQ9PR1iZNjP7lEq1WpVkpgmEHF81uKJz5v01bADANC5yclJpdNpJZPJxh8w6ct2+ghWMpnUxMSEcrmc0um0DMNQLpeTaZp7luxBcHb+W67/P//GB4dt202fOcuylMlk+P1HkOu6unnzpmq1mqanpxu3s9ks60LHAJ+1eOHz9jc+EJC5uTlf0p4f9Mba2lrL91uSXygUwg6vb21sbPiZTMY3DMM3DMPPZDL+xsZG2GENjPq/+7W1tbBDwTEqFAq+YRi+JD+TyYQdDg6wtrbmp9NpX5Jvmib5KELa+Z7GZy1aDvud8Xnz/YTv+/5xFc4AAAAAAASBObMAAAAAgNihmAUAAAAAxA7FLAAAAAAgdihmAQAAAACxQzELAAAAAIgdilkAAAAAQOxQzAIAAAAAYodiFgAAAAAQOxSzQB+qVqvyPC/sMLBDpVIJOwQAGAjkwGgh/6GXKGYx0CzLUi6XCzuMQFUqFdm2LcMwJPXna4wjz/M0NTUVdhgA0NfIgdFD/kMvUcwCAfM8T+VyWWNjYy3vn52d1djYmIaHhwNPsJ7nybIszc3NBXpcdC+TycgwDM3OzoYdCgD0JXJgNJH/0EsUs0CAbNvW6Oio3nrrLbmuu+f+YrGoubk5lUolPXjwQK7rBnq28q233lI6nZZpmoEdU9oesjU8PNz29iCfo58UCgXZth12GAAQGUH+7ScHRhf5D71CMQsEqFAoaGNjQ2+++ea+9xcKBaVSKRmGoVKppEqlEtjcntnZ2X2fG+EzTVOmaapYLIYdCgD0HXJgdJH/0CsUs+gbuVxOs7Ozsm27MYy3WCw2rn4mEglNTEwcWjjuHAZsWdaex1uWpeHhYY2NjSmXyzUeexjXdeW6rtLpdGObYRgyTVP37t2TtH1lt/5TP267w3LqDRZSqVTL+23bbsTdKpns97oty2q8b4lEQolEQtVqdd/tu59v92uov7ZEItGI5bBj7afVsQ57PfXXtHOId/15dz5n/d9TpVLZ99j15y2Xy4e+7rpMJqNSqXToawOAfnfQ3/6D/r63Mmg5kPwH/I0P9IlMJuNL8ufm5vyNjQ0/n8/7kvxUKuWvra35a2trvmmafjabbdpn5+1sNtt4/MbGhp/NZn3TNJvuz2Qy/sbGhl8qlXxJ/sbGxp5Y6vftND8/v2eb7/t+Op328/l802solUr+2tpa4zWUSqVDX3+hUPBTqdSh70s2m90Tx2Gve35+3jcMY8+x99ueyWT8dDrtb2xsNN73+fl5f2lpyTcMo/GeLS0t+fPz8wceaz8HHeuw17P7976xseFL8peWlpoek0qlfMMwGu+d7/t+Pp9vvJ6NjQ1/fn7eLxQKB77unebm5ppiAYBB1upv/0F/3/czSDmQ/Ad8iWIWfSOTyTT9kVxbW/MlNf7Q+v72H+KdyW7nH/X643cXp6Zp+nNzc77v+75hGE1/nA3DaEoAdd0Us5lMZs/96XT64Bfvf1lo77bf+7K2ttZ0+6DXfZREvvv49cdlMhm/VCr5pmm2PAFw1GJ2v2O183raTea7X0erx7Xzune/Ts4jAsC2Vn/7D8oV+xmkHEj+A77EMGP0lZ3Di+oNIHZuO3369L77VioVmabZaOdfl06nNT8/H3CkX6rVagfeb1mWFhcX2zpOMplsed/O96D+mPqwo6Bfd32o0sTERGO4kWVZqlarjSHWw8PDmpqa6mruzH7HCvL17G4kctAwtoNe9067338AQLNOcsUg5UDyH/Aliln0lVaJbL/ktls7f1ynp6dVKBTkeZ5mZ2eVTCb3nZ+zWz0p7H4ez/P2Xcanrt3XsF9hfND+vUgqpmlqY2Oj6WdtbU2GYWhtbU1zc3MyDEO5XK7j5Yn2O1aQr+eoHTH3e90AgPZ1misGJQeS/4AvUcwCf5NOpxtNmnaqVCp66aWXJG03cTIMQ6Ojo7p79+6RznTWz5bWmz1J20l0d1Oo3UqlUlsFs2maLZcDOkw7rzuI4+2UzWZVKpVUKpWa3o9O7D5Wu69n55eew66O19V/D60ac7Tzunc+1+4z5wCAZkfJFYOYA8l/AMUs0JBKpZTJZBpDY1zXVS6XU61WUz6fl7T9h/j69etaWlpSqVRq+4pp3ZtvvinbtlWtVhuLu2cymaYzoOVyubFcj23bqlQqbS01MDY21tFZ2XZet2ma8jxP1Wq18Zj9thuGoXw+L8uyGo8rl8vK5XIql8uanZ2V67ryPE/z8/ON177fc+xnv2O183qSyWTT8dtd+840TWWz2caxPc9TpVKRZVkHvu6dXNcNfA1EAIirVn/7D8oV+xmkHEj+A3YIe9IuEJTdTQ1839/TrGB3t8NW+9S79RmG0ehcvPPxkpp+6t37fH+7U9/u+3d/zAqFQuP4u5+73g0wm836hmH4pmm21cnY97e7Gbb6SLfT7OGw1+37fqOzYSqVatp3v+07X2c6nW50ckyn075hGE2dpg871n6v96BjHfR66p0gd8ZmmmbT/q3+bex+bfXn3fk7avW6dzrouAAwiHb/7T/s73srg5QDyX/AlxK+7/vHWj0DMVU/A7m0tNQ4s1itVnXz5k1dvnxZc3NzXT+HZVmS1PE6bGNjYyoUCspkMl3Hgt4YHh5WqVQ6cGg5AODoyIHRRv5DLzDMGGhTtVrV5cuXm4bIpFIpXb9+vaN5Or2Qy+UCKarRG+VyWaZpksgBoAfIgdFF/kOvUMwCbUqn06pUKioWi00t/efm5hpXVMOWz+dVq9UaLfQRLbZt6/bt22GHAQB9iRwYXeQ/9ArFLNCmVCrVaPw0OjqqRCIh27Zl27ay2WzY4TWUSiXZth2Zq8XYZlmWbNtueyknAMDRkQOjh/yHXmLOLAAAAAAgdrgyCwAAAACIHYpZAAAAAEDsUMwCAAAAAGKHYhYAAAAAEDsUswAAAACA2KGYBQAAAADEDsUsAAAAACB2KGYBAAAAALFDMQsAAAAAiB2KWQAAAABA7FDMAgAAAABih2IWAAAAABA7FLMAAAAAgNihmAUAAAAAxA7FLAAAAAAgdr4adgBR4nmePvzwQz377LM6efJk2OEAACJsa2tLn376qb773e/KMIyww4kscisAoF1Hza0Uszt8+OGHunbtWthhAABi5IMPPtAPf/jDsMOILHIrAOCo2s2tFLM7PPvss5K237zz58+HHA0AIMpWV1d17dq1Ru5Aa+RWAEC7jppbKWZ3qA9/On/+vC5evBhyNACAOGDo7MHIrQCAo2o3t9IACgAAAAAQOxSzAAAAAIDYoZgFAAAAAMQOxSwAAAAAIHYoZgEAAAAAsUMxCwAAAuc4jhKJhMbHx8MOBQDQpyhmAQBA4BzHke/7Wl5eDjsUAECfopgFAAAAAMQOxSwAAAAAIHa+GnYAaLb+ZEu1zacd7ZscOqGRUycDjggAAADAoOi0HgmjFqGYjZg79x/pnYVPOtr35v89Kuvys0fejyIYAAAAgNR5PXJr8oJen3q+BxHtj2K2j/zP//0n3f4fD468Xxj/8AAAAACgG8yZBQAAAADETl9cmZ2dndX8/LxqtZrS6bQKhULYIXXsxpVzunrpbEf7lhY/1W8f1AKOCAAAAMCg6LQeSQ6d6EE0B4t9Met5niRpfn5ekjQ8PKzr168rlUqFGVbHRk6d7Hj+au67Y3vmzG785al+87s/NG7/4IVvafhrzf/QwviHBwAAACB6uqlHjlukitlqtSrLsrS2trbnvmKxKNu2JUnT09Oam5uTJBmGoXw+L0lyXVee58k0zeMLOkL2+4f3ndHTIUQDAAAAAL0TiTmzs7OzSiQSsixLruvuub9cLsu2bS0sLGhpaUmLi4uyLKvpMZZlaWxsTKVSSYZhHFfoAAAAAIAQRKKYzefz8n1/37mutm2rUCgolUrJNE2VSiWVy+XGEGNJun37tpaWlmTbtqrV6nGFDgAAAAAIQSSK2YN4nifXdZVOpxvbTNOUYRiqVCqNbYZhKJVKKZVKNYYgAwAAAAD6U6TmzLZSH3a8ex6saZpyXVflclmmaTYaPlWr1cbc2oN89tln+vzzz5u2ra6uBhR1vKw/2VJt82lH+yaHTsRmgjgAAACA/hH5YrZWa73UTDKZ1OPHj5XP52XbtmzbVq1WUyaTUTabPfS47777rmZmZoION5bu3H+kdxY+6WjfW5MX9PrU8wFHBAAAAAAHi3wx245O1pV97bXX9jSRWl1d1bVr14IKCwAAAADQI5EvZpPJZMvttVpNp093vuTMmTNndObMmY73BwAAAACEJ/LFbH2urOu6TfNmXddlCZ6A3LhyTlcvne1o3+TQiYCjAQAAAIDDRb6YNQxDpmmqXC4rn89L2i5kPc/T9PR0IM/hOE7g82c7baoURkOlkVMnaeIEAAAAIFYiX8xK2+vM5nI5pdNpGYYhy7KUyWQCuzLrOI4cx9HKyorGx8cDOWanTZVoqAQAAAAAh4tEMVssFpXL5Rq3E4mEJMn3fUlSNpuV53manJyU53nKZDIqlUqhxAoAAAAACF8kitlsNnvocjr5fL4xzBjREafh1AAAAAD6RySK2bD1Ys5sp02V4tZQieHUAAAAAMJAMavezJmlqRIAAAAA9M5Xwg4AAABEz+zsrKampjQxMSHbtsMOBwCAPbgyi64MynBqABgknudJkubn5yVJw8PDun79ulKpVJhhAQDQhGIWXWE4NQBEW7ValWVZWltb23NfsVhsXHWdnp7W3NycpO013nev7W6a5vEFDQBAGxhmrO05s4lEIrD5sgAAhG12dlaJREKWZcl13T33l8tl2bathYUFLS0taXFxUZZlNT3GsiyNjY2pVCoFtrY7AABBSfj1xVzRaAC1vLysixcvhh1OX+t0SR+JZX0ARENccka5XJZlWdqd7sfGxmTbdmNpPNd1NTY2po2NjUbh6nmeXNeVZVkqlUodDTOOy/sEAAjfUXMGw4wRik6X9JFY1gcAulUvUtPpdGObaZoyDEOVSkWZTEbS9nDjVCqlVCqlubm5xjBkAACigGHGAAAMmPqw493zYE3TlOu6KpfLqlarje3ValUTExPHGiMAAIfhyiwAAAOmVqu13J5MJvX48WPl83nZti3btlWr1ZTJZBrDkQ/y2Wef6fPPP2/atrq6GkjMAADsRjGLUHS6pI/Esj4AcBwKhcKR93n33Xc1MzPTg2gAANiLYlbb3YxJvseLJX0AIDzJZLLl9lqtptOnT3d83Ndee21PR+TV1VVdu3at42MCALAf5sxqu5j1fV/Ly8thhwIAQM/V58ruXrLHdd2uluA5c+aMLl682PRz/vz5rmIFAGA/FLMAAAwYwzBkmqbK5XJjm+u68jxP09PTIUYGAED7GGYMAMAAsm1buVxO6XRahmHIsixlMpmurswCAHCcuDILAEAfKhaLSiQSjTmsiURCiUSicX82m1WhUNDk5KTGxsZkmqZKpVJgz+84jhKJhMbHxwM7JgAAO1HMAgDQh7LZrHzf3/OzUz6f18bGhnzfD7SQlehHAQDoPYYZi27GcbP+ZEu1zadN2zb+8lS/+d0fGrd/8MK3NPy15iV8kkMn6KAMAAAA9AmKWW0Xs47jaGVlheFQMXDn/iO9s/DJgY/51X/+fs+2W5MX9PrU870KCwAAAMAxYpgxAAAAACB2KGYBAEDgaAAFAOg1hhkjdm5cOaerl84eeb/k0InDHwQACARTeAAAvUYxi9gZOXWSRk4AAADAgGOYMQAAAAAgdihmAQAAAACxwzBjsc7soGB9WgAAAKB/UMyKJhWDgvVpAeD4cKIYANBrDDMGAACBcxxHvu9reXk57FAAAH2KYhYAAAAAEDsMM8bAaLU+bbtzZgEAAABEC8UsBsZ+69N+Z/R0CNEAAAAA6AbDjAEAAAAAsUMxCwAAAACIHYpZAAAAAEDsUMwCAIDAOY6jRCLB+u0AgJ6hmBUJFwCAoLHOLACg1yhmRcIFAAAAgLihmAUAAAAAxA7FLAAAAAAgdr4adgBA1K0/2VJt8+me7Rt/earf/O4Pjds/eOFbGv7aiabHJIdOaOTUyZ7HCAAAAAwailngEHfuP9I7C58c+rhf/efv92y7NXlBr08934uwAAAAgIHGMGMAAAAAQOxQzAIAgMCx7B0AoNcYZgwc4saVc7p66eye7Q/XN5W9s9S4XbwxoedGhpoekxw6sXs3ABgIjuPIcRytrKxQ0AIAeoJiFjjEyKmTLZs4JYdO6Nbkhcbt1Llhmj0BAAAAx4RiFujQyKmTNHcCAGBArT/Z0p37jxq3b1w5x0lt4JhRzAIAAKDvdbPUXisP1zebVju4+K1vtJxuRIEL9A7FLAAAAPpeN0vtfWc0qd8+qB24384+GnUs0Qf0Ft2MAQAAAACxQzELAAAAAIgdhhlre/mAmZmZsMMAAABAj+y31F6nc2bb2a/VEn00jgKCQzEr1sIDAACIi1aNnNotLJ//5tdbHvM7o6c7iuWw/dafbOnjP37RtK2dxlESzaOAdlDMAgAAIDbaaeTUqolTGM2Y2om1VeMoqXW8XNUFmlHMAgCAwDGFB+jOX57+n46u6nJFF4OEYhYAAASOKTxAd/7n//6Tbv+PBwc+huWAMOgoZgEAABAbrRo5ddqMqdc6jVWSSoufHrq2LTDoKGaBHuqmSQVDhAAA2Gvk1MmWObLTJk691E2sXzvB13TgMHxKgB6KU5MKAAAQHXG6Ag2EhWIWAAAAiJg4XYEGwvKVsAMAAAAAAOCouDIL9FCrIUIP1zebug8Wb0y0bKsPAAAAYH8Us0APtRoilBw6oVuTFxq3U+eGafYEABgorRokSjRJDEKnzScl3lvED8UscMxGTp2kuRMAYKC10yBRokliJzptPinx3iJ+mDMLAAAAAIgdrswCAAAAaGn9yZbu3H/UuH3jyjmGIiMyKGYBAABwrFo1SJRokhiETptP1n38xy/27Ltz2PLFb32j5e+EAhdhoJgFYoSzowDiwnEczczMhB0GImq/NVRpkti9bppP/mL+40Pn2+4siuuYa4uw9EUxa9u2KpWKPM9TLpdTPp8POySgK/t1eeTsKIC4cBxHjuNoZWVF4+PjYYeDmKBJYm+E8b5yAh7HIbBi9qOPPmq5/cUXXwzqKVoql8vyPE9LS9tnicbGxpRKpZROp3v6vEAvtdvlkbOjQH8LK7cCwFH85en/YXgyQtF1MXv79m298sorkiTf95vuSyQS+utf/9r2sarVqizL0tra2p77isWibNuWJE1PT2tubk6SlE6nmwrXdDqtarVKMQsAiK0gcysAtKvT+balxU/197/4fw48dr+egOcKdLi6LmZt29bNmzdl27aSyWRHx5idnZVt2zJNU67r7rm/XC7Ltm0tLCzIpLWmagAAIABJREFUMAxZliXLslQqlWQYRtNjK5WKcrlcR3EAABAFQeRWADiqTufbfu1EX8xcPFSraWBxugLdj4V3IP/y3njjDT333HMd75/P55XP51Uul2VZ1p77bdtWoVBQKpWSJJVKJY2Njcnz/v/27ie2jfTM8/iPswMbGLnTJUpww4Md2Cq2+xLDmVBuDPqwOwOIQh+MxTgxKZ186piMfWkEk2G150RjgWlTfQj64obo3HyyyfRksLeQCjDYQ2MxEqdnYWMXaalsI4sFxpaoSmINIO9iag9aMqJEiv+KrCry+wGIpIrF4vN2mfXoqXrrfZ2mYtayLCWTycZ2QFgxyiOAQXMrAHiB55h/r5vHwIJyB7rfwlsKTvHdjYGL2XQ6rZ/97Gf6q7/6Ky/iOcZxHNm23dRt2DRNGYahSqWiZDIp6aCQjcViSqfTQ4kDGKV2ozwedWF2Su+989YIIgIwSsPOrQDgpVYX4Xf/9Y3+yz//78byf/rOH2v6j5ovuLe7AD+OdxBHrd/CWwpX929P7sz+9V//tf72b/9WV65cOdbt99GjRwPtu97t2DTNpvWHuyRnMhmlUimek8XYO9rVh7uwwPgaZm4FAC+1uwj/Z3MzJ35u+/X+sYGjpHB13fUDxf7vDVzM2rbduDsqHR+oYlC1Wq3l+mg0qp2dHZVKJRUKBVUqlcZ73UzP8/LlS7169app3ebm5uABA0NEVx9gMgw7twJAEHg9e8Owi7x+B8hqdfOh21j77S78r2/+bxctCr+Bi9nHjx97EUffkslkX0n+/v37TOYOAAgkv3MrAASdH9MB9TtAVqs70N0+v1pc/7Ue/NdnJ8bVqti/+R/m9Isf/cdj39mp8K63KSwCP/RYu1Eca7WaZmZO7rpwktu3bx8bbGpzc1PXrl3re58AAACTji6QGIX//r9+01eR5/XzoN30mhvk+dU/m+tvRPs/OvWHx8ZV6abwDpuei9kPP/xQqVRKP/jBDyRJt27dOnH7L774or/I/r/6s7K2bTc9N2vb9rFniHpx9uxZnT17dqDYAADwwqhzK+CVsE9VAn8NMntDcf3X+m/PWj+OiNbG8XG1novZnZ2dE5e9ZhiGTNNUqVRqPAdr27Ycx9HS0pIn35HL5ehyjLHGVXIg2EadWwGvhGmqEgRPu4Gj/JjbNqh/K13+92/rP1+71LSOqRp/r+d/Bevr603Lo3iux7IsZTIZJRIJGYahVCqlZDI50J3Zw3K5nHK5nJ4+fapLly51/gAQMrW9N01/bFy9fC4QJ2gAB/zIrQAQVN3cQex3MCZJIx9BeZBYW33nOHYX7lcgnpktFArKZDKN5UgkIun3ozem02k5jqOFhQU5jqNkMqlisehLrEDQtevyddKyRJcvAAAQHu3u6h52YXbq2HOjPyn/ytMRlLvR78BRJ+2Png0HPClmv/zyS62urjauLF+5ckUrKyv6zne+09Xn0+m00un0idtks9mO0+0AoMsXMC4Gza3AKHg5VQkwqKMFYpD/nVGQemPgYvazzz6TZVlKJpO6d++epIPuUvF4XKVSSd/73vcGDhIAgElCbkVYeH3HCRgEBeLkGbiY/fTTT7WysqIf//jHjXU3b97U4uKistlsKBIuA0ABAIJkHHIrJhcFBYJskBGUg3ynd1J50s04mUy2XPfJJ594sfuhYwAojBO6fAHjIey5FQCCqJtnbaXWz9sieAYuZu/cuaNKpdKYG6/u+fPnmp+fH3T3AHrU74AIAIKD3AoAoxWm523xez0Xs0cnct/Z2dHPfvYzlcvlpvWVSkWmaQ4WHQBPcIIGgm0ccyuP8AAIE7rHh1PPxWyridyvX7/emEanbmFhof+oRoyEi3HHCRoItnHNrTzCAwAYpp6L2XGcyJ2ECwDw0zjmVgAAhs2TAaAAjKft1/t6+NWLxvKND84zvQIATAhyAICgo5gF0FZt740+X/umsXz18jn+kAGAMbT9el+1vTdN655v7zXlgG//8bdajoRPXgDgF4pZAACACffwqxdNhWsrh6d4q/t44SJjMgDwDcUsAEntr8qftFzHlXkAAACMGsWsGM0YkPq/Ki9xZR4AAACjRzErRjMGAACT7cYH53X18rmmdc+395ouYhZuzLd8ZhYA/DJwMfvTn/5US0tL+ta3vuVFPAAATDxyK0Zt9szpjo+LXJid0nvvvDWiiACgs4GL2Xv37ikSieijjz7yIh4APun3qrzElXnAa+RWBEF06pQ+XrjYtAwAQeJJMZvJZJRIJHT+/HkvYgLgA67KA8FBbkUQzJ45zXgIAAJt4GJ2fX1dc3NzMk1TiURCpmk2vf/FF18M+hUAAEwUcisAAJ0NXMzati3TNBuJdmdnZ+CgRo3RjIHW6GIG+GMccisAAMM2cDH7+PFjL+LwFaMZA63RxQzwxzjkVgAAhu0PvNjJl19+qaWlJb333u//6L1165Z++ctferF7AAAmDrkVAICTDVzMfvbZZ7IsS5lMRltbW431CwsLyufzg+4eAICJQ24FAKAzT0YzXltb05/+6Z82rZ+fn9fy8vKguwcAYOKQWwEA6GzgYtZ1XU1PTx9bb9u25ubmBt09gJDafr2vh1+9aCzf+OB8x6l/ABwgtwIA0NnA3Yxv3rypdDqt3/72t411z5490w9/+EP98Ic/HHT3AEKqtvdGn69903jV9t74HRIQGuRWAAA6G7iYzefzunDhggzDkOu6unjxot59910lEgn9+Mc/9iJGAAAmCrkVAIDOBu5mLEmrq6uyLEv/9E//JEmKx+Oh6gbFPLPAYLZf7x+78/p8e+/EZelg3lq6HgOthT23AgAwbJ4Us5KaJncPG+aZBQbz8KsX+nztmxO3ST/cOLbu44WLzGMLnCDMuRUAgGHzbJ7ZDz/8UDMzM5qZmdGHH36of/7nf/Zi1wAmzPbrff2k/KvGa/v1vt8hAb4gtwIAcLKBi9lPPvlEqVRKruvq3r17unfvnv7t3/5N8Xhcf/d3f+dFjAAmCANHAeRWeIcLhADG2cDdjAuFgrLZrD799NPGups3b8qyLGWzWX3ve98b9CsABNyND87r6uVzTeueb+81dS0u3JjXhdmppm2iU6dGEh8QNkHIrZZlqVKpyHEcZTIZZbPZoX8nvFe/QFh39fI5xioAMDYGLmaj0agymcyx9X/zN3+jBw8eDLp7ACEwe+Z0xz+OLsxO6b133mpat/16X7/6l981retm4CiJwaMw3vzOraVSSY7jaGPj4IJULBZTPB5XIpEY+nejf/0OxidxTgUQTgMXs8lkUmtra/roo4+a1q+trWl5eXnQ3QMYY/0OHCUxeBTGm5e5tVqtKpVKaWtr69h7hUJBlmVJkpaWlrS6uipJSiQSTYVrIpFQtVqlmA04zqkAJs3Axaxt2/rss8/0i1/8orHOcRxVKhUlEommpPvo0aNBvw4AgLHnRW5dWVmRZVkyTVO2bR97v1QqybIsra2tyTAMpVIppVIpFYtFGYbRtG2lUml5pxgAAD95Mprx9evX5bpu4/X222/r+vXrevvtt5vWAwCA7gyaW7PZrFzXVT6fb/m+ZVnK5/OKx+MyTVPFYrHRvfjodslkUvF43NP2AQAwqIHvzD5+/NiLOACMmejUKX28cLFp+ah+B45qtz9gXAw7tzqOI9u2m7oNm6YpwzBUqVSUTCYlHRSysVhM6XR6qPHAG5xTAUyagYtZAGhl9szpjs9f9TtwFIDB1Lsdm6bZtP5wl+RMJqNUKtXTc7IvX77Uq1evmtZtbm4OGC26xTkVwKShmJWUy+V09+5dv8MAAGAkarVay/XRaFQ7OzsqlUoqFAqqVCqN97qZnuf+/fvkUwDAyFDM6qCYzeVyevr0qS5duuR3OAAA+CqZTPY11sXt27eVSqWa1m1uburatWtehYYedfPIBwCEFcUsAAATJhqNtlxfq9U0MzPT937Pnj2rs2fP9v15eK+bRz4AIKwoZgEECncRgOGrPytr23bTc7O2bR+blgcAgKAaeGqeL7/8Ul9++WVj+datW5qZmdH777+v58+fD7p7ABOmfheh/uo0mAkwjoadWw3DkGmaKpVKjXW2bctxHC0tLQ28fwAARmHgYtayrMZV3AcPHqhQKKhQKCgejx97bgYAAHQ2itxqWZYsy1K1WpVt20qlUkomk57dmc3lcopEIoxFAQAYmoG7GW9tbenKlSuSpGKxqHQ6revXr+u73/2uLl682OHTAADgKC9ya6FQUCaTaSxHIhFJagzslE6n5TiOFhYW5DiOksmkisWiZ21gcEUAwLANfGfWNE3t7u5KkiqVihYXFyVJv/nNb3juBgCAPniRW9PptFzXPfY6LJvNand3V67relrIAgAwCgPfma1fKZ6ZmZFpmvr+978v6SD5LiwsDBwgAACThtwKAEBnAxez+Xxey8vLevbsmRKJRGO9aZpNywAAoDvkVgAAOvNkap54PK54PN607vr1617sGgCAiURuBQDgZD0Xs7du3epp+y+++KLXrwAAYKKMY27N5XK6e/eu32EAAMZYz8Xszs7OMOIAAGBijWNuZTRjb22/3tfDr140lm98cJ55uAFMvJ6L2cePHw8jDgAAJha5FZ3U9t7o87VvGstXL5+jmAUw8QaemmccMLE7AAAAAISLJwNAff3113r06JFs2z723qNHj7z4iqGiKxQAIGjCnlsBABi2ge/MPnjwQPF4XOVyWa7rqlQqyXVdlctlOY7jRYwA0JXt1/v6SflXjdf2632/QwL6Qm6dbNuv9/Wrf/ld0+v59l7TNs+3945twzkPwKQZ+M7sJ598onK53JjEPRqN6vHjx6pUKnrw4MHAAQJAt3imDONiHHIroxn37+FXL5rOZa2kH24cW/fxwkX9aPG9YYUFAIEz8J3Z3d1dvf/++41l0zT1u9/9Tu+//74qlcqguwcAYOKMQ27N5XJyXVdPnjzxOxQAwJgauJiNx+NNifXKlSt6/Pix1tbW6AoFAEAfyK0AAHQ2cDfjfD6vcrms73//+5Iky7IUi8UUiUSUTqcHDhAAgElDbp1sNz44r6uXzzWte76919S1uHBjXhdmp5q2iU6dGkl8ABAUAxezCwsLjWd6JGlubk67u7uq1Wqam5sbdPcAAEwccutkmz1zuuPz/hdmp/TeO2+NKCIACCZPpuY56u2339bbb789jF0DADCRyK0AADTruZj98MMPlUql9IMf/ECSdOvWrRO3/+KLL/qLDABOsP16X7W9N03rWk1dcVR06hQjHCNwyK0AAPSu52J2Z2fnxGUAGAWmrsA4GcfcytQ83opOndLHCxeblgFg0vVczK6vrzctP3782LNgAACYROOYW3O5nHK5nJ4+fapLly75HU7ozZ45zYU4ADhi4Kl5fvrTn+q3v/2tF7EAAACRWwEA6MbAA0Ddu3dPkUhEH330kRfxAEBXmLoC44zcCgBAZ54Us5lMRolEQufPn/ciJgDoiKkrMM7IrQAAdDZwMbu+vq65uTmZpqlEIiHTNJveZ8RFAAB6Q24FAKCzgYtZ27ZlmmYj0Y7DCIwAAPiJ3AoAQGcDF7PjMOIiAABBQm4FAKCzgUcz/uUvf9ly/ddff62vv/560N13zbZtZTKZkX0fAADDEpTcCgBAkA1czC4uLrZcv7OzI8uyBt191zGkUinZtj2S7wOAuu3X+/pJ+VeN1/brfb9DwhgIQm4dVC6XUyQSYY5ZAMDQDFzMuq7b9r2jk8B3Uq1WFYvFWr5XKBQ0PT2t6enpY3dgy+Wy7ty509N3AYAXantv9PnaN41Xbe+N3yFhDHiZW/2Sy+Xkuq6ePHnidygAgDHV9zOz0WhUkUhEkUhEMzMzTe85jiPXdTU/P9/VvlZWVmRZlkzTbHl3tVQqybIsra2tyTAMpVIppVIpFYvFfsMHMIaiU6f08cLFpmUgTLzMrQAAjLu+i9m1tTW5rqsrV660HKjCNE3Nzc11ta9sNqtsNqtSqaRUKnXsfcuylM/nFY/HJUnFYlGxWEyO48gwjH6bAGDMzJ45rR8tvud3GEDfvMytAACMu76L2e9+97uSpEQioYWFBc8COspxHNm2rUQi0VhnmqYMw1ClUlEymRzadwMAMEqjyq0AAIyDgafm+cUvfuFFHG3Vux0fnTC+XZfkbr18+VKvXr1qWre5udn3/gAA8MqwcysAAONg4GJ22Gq1Wsv10Wi0MYl8JpPR+vq6bNvW4uJiU5fkdu7fv6+7d+96Hi8AAAAAYPgCX8x2Y3V1tefP3L59+9jzuZubm7p27ZpXYQEAAAAAhiTwxWw0Gm25vlarHRvpsRdnz57V2bNn+/48gMmy/Xq/5bQ7z7f3TlyWDkZVnj1zemixAQAATKLAF7P1Z2Vt2256bta2bUYyBjAyD796oc/Xvum4XfrhxrF1Hy9cZJRlAAAAj/2B3wF0YhiGTNNUqVRqrLNtW47jaGlpyZPvyOVyikQiunTpkif7AwAAAAAMV+CLWelgnlnLslStVmXbtlKplJLJpGd3ZnO5nFzX1ZMnTzzZHwAAk44LxQCAYQtEN+NCoaBMJtNYjkQikiTXdSVJ6XRajuNoYWFBjuMomUyqWCz6EiuAyXTjg/O6evncsfXPt/eauhYXbszrwuxU0zbRqVNDjw8Imlwup1wup6dPn1LQAgCGIhDFbDqdVjqdPnGbbDarbDY7oogAoNnsmdNdDeJ0YXZK773z1ggiAgAAmGyh6GY8bHSFAgAAAIBwoZgVz8wCAAAAQNhQzAIAAAAAQodiFgAAAAAQOoEYAAoA/LT9el8Pv3rRWL7xwfmuBnsCAACAfyhmdfDM7N27d/0OA4BPantv9PnaN43lq5fPUcwCAAAEHN2MxQBQAAAAABA2FLMAAAAAgNChmAUAAAAAhA7FLAAAAAAgdChmdfDMbCQS0aVLl/wOBUDIRKdO6eOFi41XdOqU3yEBgUBuBQAMG8WsGAAKQP9mz5zWjxbfa7wYBRk4QG4FAAwbxSwAAAAAIHQoZgEAAAAAoUMxCwAAAAAInT/0OwAAGKXt1/uq7b1pWvd8e+/EZelgoCeehwUAAAgOilkAE+XhVy/0+do3J26TfrhxbN3HCxf1o8X3hhUWgD5tv97Xw69eNJZvfHCeC08AMCEoZnUw4uLdu3f9DgMAAPSotvem6QLV1cvnKGYBYELwzKyYPgAAAAAAwoY7swAmyo0Pzuvq5XNN655v7zV1LS7cmNeF2ammbaJTp0YSXzfC1K0yTLECAIBwoZgFMFFmz5zuWExdmJ3Se++8NaKIehembpVhihXoFhdpACAYKGYBAAB6wEUaAAgGnpkFAAAAAIQOd2YBAEAoME80AOAwilkAABAKzBMNADiMbsY6mJonEono0qVLfocCAAAAAOgCxayYZxYAAAAAwoZuxgAAIBT8mCea53QBILgoZgEAQCj4MU80z+kCQHDRzRgAAHiO8SgAAMNGMQsAPth+va+flH/VeG2/3vc7JMBTjEcBABg2uhkDgA9qe2+aui5evXyO5+uAAPLjOV0AQHcoZgEAANrw4zldAEB3KGYBIKBajaIqBXckVUZ9BQAAo0QxCwAB1c0oqlJwRlJl1FcAADBKFLM6GKTi7t27focBAADG2PbrfT386kVj+cYH5+mVAAADoJjVQTGby+X09OlTphAAACBEolOn9PHCxabloGLgNwDwFsUsAARUq1FUpeCOpMqor/DD7JnTdFMHgAlFMQtg4gX1zk43o6hKwRlJlVFfAQDAKFHMAph43NkBAAAInz/wOwAAAAAAAHrFnVkAAIAeBPXRBACYNBSzAAAAPeDRBAAIBroZAwAAAABChzuzAAAAHtt+va/a3pumdc+3905crotOnWL+WQDoAsUsAACYSNuv9/XwqxeN5RsfnPesiHz41Qt9vvbNidscnoP5sI8XLtKNGQC6QDELAAAmUm3vTVPBefXyOe6IAkCI8MwsAAAAACB0uDMLAEPW77NzPDcHhNeND87r6uVzTeueb+81dS0u3JjXhdmpY59lqh8A6A7FLAAMWb/PzvHcHBBes2dOd7wYdWF2Su+989aIIgKA8UMxKymXy+nu3bt+hwEAQ9XvYDfDHCRnGMIWLwAA6A/FrA6K2Vwup6dPn+rSpUt+hwMAQ9HvYDdhGyQnbPECAID+UMwCwJD1++wcz80BAAC0RzELAEPGs3MAAADeY2oeAAAAAEDocGcWAAC0ZNu28vm8VldX/Q5lIK2mx5KYIgsAwo5iFgAAHLO4uKharaZoNOp3KAPrZnosiSmyACBs6GYMAMAYq1arisViLd8rFAqanp7W9PS0MplM03vlcll37twZRYgAAPSFYhYAgDG0srKiSCSiVCol27aPvV8qlWRZltbW1rSxsaH19XWlUikfIgUAoD90MwYAYAxls1lls1mVSqWWRaplWcrn84rH45KkYrGoWCwmx3FkGMaowx2qVtNjSUyRBQBhRzELAMCEcRxHtm0rkUg01pmmKcMwVKlUlEwmfYzOe91MjyUxRRYAhA3FLAAAE6be7dg0zab1pmm27JLcrZcvX+rVq1dN6zY3N/ve37iJTp3SxwsXm5YBAP2jmAUAYMLUarWW66PRqHZ2diRJmUxG6+vrsm1bi4uLTV2S27l//77u3r3rebzjYvbMaUZGBgAPUcwCAIBj+plb9vbt28eez93c3NS1a9e8CgsAgAaKWQAAJky7uWNrtZpmZmb63u/Zs2d19uzZvj8PAEAvmJoHAIAJU39W9ujzsbZtj91IxgCA8cWdWQAImTANIhOmWCeJYRgyTVOlUknZbFbSQSHrOI6WlpZ8jg4AgO6MRTFbKpX06aefSpKWl5cbiRkAxlGYBpEJU6yTxrIsZTIZJRIJGYahVCqlZDLp2Z3ZXC7HYFAAgKEKfTdjx3FkWZY2Nja0sbGhR48eqVqt+h0WAAC+KhQKikQijQGZIpGIIpFI4/10Oq18Pq+FhQXFYjGZpqlisejZ9+dyObmuqydPnni2TwAADgtUMVutVhWLxVq+VygUND09renpaWUymcb6SqXSNFXA8vKyHj16NPRYAQAIsnQ6Ldd1j70Oy2az2t3dleu6nhayAACMQiCK2ZWVlcbV41aTtZdKJVmWpbW1NW1sbGh9fb1xpdm27aZJ3w3DGGjCdwAAAABA8AWimM1ms3JdV/l8vuX7lmU1Jmuvd4MqlUpyHEeSGv8rtZ9uAAAAAAAwPgI/AJTjOLJtW4lEorHONE0ZhqFKpSLDMFSr1RrvHb1TCwCTaPv1vmp7b5rWPd/eO3G5nW4/F506pdkzp3uIcnD9ttOPWAEAgLcCX8zWuwwfLVBN05Rt20qn07IsS47jyDAMPXr0SA8ePOi435cvX+rVq1dN6zY3N70LHAB89PCrF/p87ZsTt0k/3Di27s/movpvz2ottj75c5L08cLFkY9c3G87/Yh10jCaMQBg2AJfzB6+63pYNBrVzs6ODMNQsVjUwsKCpIMBoA4PCNXO/fv3SbIAAAxJLpdTLpfT06dPdenSJb/DaYl5kAEg3AJfzHYjkUhoY6P1nYJ2bt++3RhEqm5zc1PXrl3zMjQAABBQzIMMAOEW+GK23YBOtVpNMzMzfe/37NmzOnv2bN+fB4Agu/HBeV29fK5p3fPtvaYut4Ub87owO9VxX91+zo+7Wv22kztwAACEX+CL2fqzskcHdrJtW4Zh+BUWAATa7JnTHQc4ujA7pffeeavnfff7uWEYZjsBAECwBWJqnpMYhiHTNFUqlRrrbNuW4zhaWlry5DtyuZwikUhgn+kBAAAAADQLfDErHcwza1mWqtWqbNtWKpVSMpn07M5sLpeT67p68uSJJ/sDAGDScaEYADBsgShmC4WCIpFIY0CmSCSiSCTSeD+dTiufz2thYUGxWEymaapYLPoVLgAA6IALxQCAYQtEMZtOp+W67rHXYdlsVru7u3Jdl0IWAAAAACZcIIpZv9EVCgAAAADChWJWdIUCAAAAgLChmAUAAAAAhA7FLAAAAAAgdP7Q7wCCZH9/X5K0ubnpcyQAxt3z7T29efWisfzN//wf+j/bU4H8Tj9iHcSo4q3ninruQLNcLqe7d+82lsmtAIBOes2tEffosMET6GjCBQCgWz//+c/1l3/5l36HEVh///d/r2vXrvkdBgAgRLrNrRSzhziOo3/4h3/Qn/zJn+j06dN972dzc1PXrl3Tz3/+c7377rseRohecByCgeMQDBwH7+3v7+vXv/61/vzP/1yGYfgdTmCRW/tDe8fbpLVXmrw2097+9Jpb6WZ8iGEYnl5df/fdd/Xtb3/bs/2hPxyHYOA4BAPHwVvxeNzvEAKP3DoY2jveJq290uS1mfb2rpfcygBQAAAAAIDQoZgFAAAAAIQOxSwAAAAAIHT+XS6Xy/kdxDiamprSX/zFX2hqKrjTV0wCjkMwcByCgeOAsJu0f8O0d7xNWnulyWsz7R0+RjMGAAAAAIQO3YwBAAAAAKFDMQsAAAAACB2KWQAAAABA6FDMAgAAAABCh2IWAAAAABA6FLMeKxQKmp6e1vT0tDKZjN/hTKRUKqVIJNL0isVifoc19qrVatv/zvwuRqfdceB3gaA76RzSStjPK720dxx+v47jKJPJ9HzMwnqc+2lv2I/z4famUik5jtPxM2E9vlLv7Q378a2zbVuFQiEwx5di1kOlUkmWZWltbU0bGxtaX19XKpXyO6yJlE6ntbW11XiVy2W/QxpbKysrikQiSqVSsm372Pv8Lkaj03GQ+F0gmLr5t3tUmM8r/bRXCv/vt358NjY29ODBAz1+/FiLi4snfibMx7mf9krhPc7z8/OS1DhWjuNoYWHhxM+E+fj2014pvMf3sFQqpUwmo1qtduJ2Izu+Ljxjmqa7urraWN7a2nIlubu7uz5GNXmSyaSbzWb9DmPiFItFt9Uphd/FaLU7DvwuEHTt/u22Mg6lpaWTAAALT0lEQVTnlV7aG/bf79bWlmsYRtO6jY0NV5K7sbHR9nNhPc79tjfMx/nwcXJd1y2Xyx2PVViPr+v2194wH9+61dVVN5lMupLcra2tE7cd1fHlzqxHHMeRbdtKJBKNdaZpyjAMVSoVHyMD/MPvAoDXOK+Ej2maevDgQdO6eDwuSVpfX2/5mTAf537aG3bpdLppuVwuyzAMGYbRcvswH1+p9/aOC8uydOfOnY7bjfL4Usx6pN5VyDTNpvWmafbUjQjeqD+LVH+OAf7gdxEs/C4wDib1vBL2328ymWy5PhqNtlwf9uPca3vrwn6cJalSqWhlZeXEoifsx/ewbtpbF+bja1mWrly50rgwc5JRHl+KWY+06zcejUa1s7Mz4mhQqVS0urqqjY0N2bbd1XMq8B6/i2Dhd4FxMKnnlXH7/RYKBRmG0bboG7fj3Km9dWE+ztVqVZFIRIuLi0okEspms223HYfj20t768J6fB3H0crKivL5fFfbj/L4Usxi7CwvL6tcLiuRSMg0TeXzeVUqldBd6QO8xO8CCK9x+/3ati3LslQsFv0OZSS6bW/Yj3M8Htfu7q42NjYkKTSFWr96bW+Yj+/NmzeVTCa7uis7ahSzHmnXbaRWq2lmZmbE0Uy2ZDLZ1Ef/ypUrkhSKZzDGDb+L4OB3gXExieeVcfr9Oo6jxcVFPXjwoKlNR43Lce62vdJ4HGfDMBSPx1UsFlWpVFQqlVpuNy7Ht9v2SuE9vtVqVaVS6dhz4CcZ5fGlmPVIvU/40asrtm2P/cPgQVfv6tDpORV4j99FcPG7QFhxXgnv79dxHM3Pz2t1dbVjd9txOM69tLeVsB5nSY3BkP7xH/+x5fvjcHwP69TeVsJyfB89eiRJmpuba8wZK0mxWKwxRdFRozy+FLMeMQxDpmk2XZGxbVuO42hpacnHyCbP0ati1WpVkgLZNWLc8bsIDn4XGBeTeF4Zh99vfR7O1dXVjncopfAf517bK4X3ONu2faxocRxHjuMoFou1/EyYj28/7ZXCe3zz+bx2d3f17NkzPXv2TGtra5IORnCu//+jRnp8PZ3oZ8Ktrq425hDb2tpy4/G4m0wm/Q5rotTndltdXXV3d3fdjY0N1zAMN51O+x3a2Gs3ZyK/i9FqdRz4XSAM2p1Dtra2js3pOA7nlW7bOy6/33g87mazWXdjY6PpVZ9zctyOc6/tDfNx3t3ddU3TdPP5fONYJRIJ1zTNxjbjdHz7aW+Yj+9R9fliD88z6+fxpZj1WD6fdw3DcCWF4gc5juonFUmNkw2Gp36yOvo6jN/F8HU6DvwuEFSd/u3W3z8qrOeVftob9t9v/Y/fVq96W8bpOPfb3jAf593dXTeZTLqGYbiGYbjJZLJRuLvueB1f1+2vvWE+voe1Kmb9PL4R13VdL+7wAgAAAAAwKjwzCwAAAAAIHYpZAAAAAEDoUMwCAAAAAEKHYhYAAAAAEDoUswAAAACA0KGYBQAAAACEDsUsAAAAACB0KGYBAAAAAKFDMQv4qFqtynEcv8PAIZVKxe8QAAABQZ4OHvI0DqOYBXxSqVRkWZYMw5AkpVIpZTIZn6OC4zhaXFz0OwwAgM/I08FEnsZhFLOADxzHUSqV0urqqt+h4IhkMinDMLSysuJ3KAAAn5Cng4s8jcMoZgEffPrpp0okEjJN09P9VqtVTU9Pd73ey+8YJ/l8XpZl+R0GAMAn5OlgI0+jjmIW8MHKyoru3LnjdxhowzRNmaapQqHgdygAAB+Qp4ONPI06ilmgg0wmI8uyZFmWpqenFYvFjp08LctSLBZTJBJp+f5h9YEL4vF4y/dP+h7pIMHGYjFNT08rlUo1BqZIpVKan5+X4ziKRCKKRCKqVqtt1x/9vunp6aYuO63a1Glf7Zz036dde+ptOvx8Uv17D39nJpPRysqKKpVK233Xv7dUKnVsd10ymVSxWOzYNgCYRKlUqpGrMplM4zwuSaVSSfPz8y3PvdLvz9v13DA9Pa1CoSDbtrW4uKhIJNLINYN8plMc7ZCnydMIERfAiZLJpCvJXV1ddXd3d910Ou0e/ulsbGy4hmG4u7u7jeVyudx2f/l83o3H4z1/j+u6bjqdduPxuLu1tdXYxjTNxvvlctk1DOPYvtutTyaTbiKRcHd3d92trS3XNE23XC6f2KZ2+2rnpH11ak8ymXTT6XRjeXd315XkbmxsNG0Tj8ddwzAa/+1c13Wz2WyjPbu7u265XHbz+fyJ7T5sdXW1KRYAwIF0Ou0mk0l3d3fXLRaLrqTGudd1D86/xWLRdd2DnNHqvH0432WzWVdSIx/Uz8uHz//9fKZTHO2Qp8nTCA+KWaCDZDLZdLLc2tpyJblbW1uu67pusVh0TdNsSuQnqf8R0Ov31JePfo9pmu7q6qrrur0lyaP7r2+XTCZPbFOvSbLdvrppT7dJ8mg7Wm3XTbuPtpPrfQBwnGEYTYWFYRgnFonxeLxRpLhu+3x3eJtsNttUUPbzmU5xtEOebv5+8jSCjG7GQBcOdzWKRqOS1Ohmk0gkJEnT09NaXFzs+PxGrVZr7KOX76lUKjJNszFFQF0ikVC5XO6lOZLU6AI0Pz/f6MaTSqVUrVZ7btNJ2u3Ly/YcHaTjpC5iJ7X7sKP//QEA/YlGo9ra2mpad/j8XD9/H143MzNzbD/9fKZTHK2Qp8nTCA+KWaAL7ZKaJBmGoa2tLa2ursowDGUymY7z0NVqtZ6/Zxgna9M0tbu72/Ta2trqq03ttNuXl+3pdbTJdu0GAHS2tLSkfD4vx3G0srKiaDTaVJQ4jiPLsjQ/P69YLNYoXA5rle9OyoH9fKabONohT5OnEQ4Us4BH0um0isWiisWiHj9+3HY70zRl23bP+08kErJt+9hnK5WK3n//fc/2d1i3berG0X11257Df1C0++PiqPofVa0Gveim3Ye/6+gVaQCYdLZtyzAMzc3N6dGjR0136RzH0dzcnGZmZrS2tqatra3Gnb9RGiQO8jR5GuFBMQsMqFQqaWVlRbZty3EclcvlE69AxmKxvq52xuNxJZPJRlcb27aVyWRUq9WUzWYlHSRgx3FUrVYb27RbbxiGstmsUqlUY7tSqaRMJnNim9p9R6//fbppTzQabdp/t3PKmaapdDrd2LfjOKpUKkqlUie2+zDbtj2fXxAAxkGtVtPy8rI2NjZULBab7lbWajU5jtPonlqtVrW+vu5LjP3GQZ4mTyM8KGaBAZmmqXK53Hi2Y319/cSh4q9cudLXFV9JKhaLSiQSjaH3a7Wanj171hRLPB7XwsKCbt682UjG7dbn83ktLy9rcXFR09PTWl1dVSqVOrFN7fbVz3+fTu3JZDKybbvxHE8mk2n5/E4rq6urymQyjekjLMvS8vLyie0+rFwu+3I3AQCCzjRNpVIpxWKxxpQti4uLjeKxXqTUz6+JRGLkd88GiYM8TZ5GeERc13X9DgKYNLFYTPl8Xslk0u9Q0Mb09HQjiQMADtTvnm1sbDTuilWrVd28eVNXrlzR6uqqzxF6gzwdfORpSNyZBXyRyWTGJuGPo1KpJNM0SZAAcES1WtWVK1eaunfG43EtLy/3fTcziMjTwUaeRh3FLOCDbDarWq3W08iKGB3LsvTgwQO/wwCAwEkkEqpUKioUCk1T0rTqBhpm5OlgI0+jjm7GgE9s21YqlVKxWGQAgwBJpVJaXFxUOp32OxQACKRqtSrLsrS+vi7HcRSPx5XJZMbuvEmeDibyNA6jmAUAAAAAhA7djAEAAAAAoUMxCwAAAAAIHYpZAAAAAEDoUMwCAAAAAEKHYhYAAAAAEDoUswAAAACA0KGYBQAAAACEDsUsAAAAACB0KGYBAAAAAKFDMQsAAAAACJ3/ByvjE4DcP6hgAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(2, 2, figsize=(8,6))\n", "axs = np.ravel(axs)\n", "names = 'mlog10p ts ns gamma'.split()\n", "for (i, name) in enumerate(names):\n", " ax = axs[i]\n", " mask = bg_trials_hottest.ns > 0 if name in 'ns gamma' else bg_trials_hottest.ns >= 0\n", " hl.plot1d(ax, hl.hist(bg_trials_hottest[name][mask], bins=25), crosses=True)\n", " ax.set_xlabel(name + ' (hottest source)')\n", " ax.set_ylabel(r'trials per bin')\n", " ax.semilogy()\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## p-values from single batch" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Above we described the most intuitive way to set all this up. Technically, the initial trial batches may be redundant if equal statistics are required for the per-hypothesis background distributions and the final distributions used for post-trials calculations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can do all the work with one batch; let's see how. First we run some trials:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Performing 10000 background trials using 15 cores:\n", " 10000/10000 trials complete. \n", "CPU times: user 52.9 ms, sys: 149 ms, total: 202 ms\n", "Wall time: 12.2 s\n" ] } ], "source": [ "%time t = multr_nobg.get_many_fits(10000, seed=1, mp_cpus=15)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
gamma_0000gamma_0001mlog10p_0000mlog10p_0001ns_0000ns_0001ts_0000ts_0001
03.2500003.2500000.0000000.0000000.0000000.0000000.0000000.000000
14.0000004.0000004.9381290.0000005.1841950.0000004.9381290.000000
23.2362522.8424870.0306010.0469030.3320280.9210090.0306010.046903
32.0434734.0000001.2553030.0000002.6146850.0000001.2553030.000000
42.5000002.0000000.0000000.0000000.0000000.0000000.0000000.000000
\n", "
" ], "text/plain": [ " gamma_0000 gamma_0001 mlog10p_0000 mlog10p_0001 ns_0000 ns_0001 \\\n", "0 3.250000 3.250000 0.000000 0.000000 0.000000 0.000000 \n", "1 4.000000 4.000000 4.938129 0.000000 5.184195 0.000000 \n", "2 3.236252 2.842487 0.030601 0.046903 0.332028 0.921009 \n", "3 2.043473 4.000000 1.255303 0.000000 2.614685 0.000000 \n", "4 2.500000 2.000000 0.000000 0.000000 0.000000 0.000000 \n", "\n", " ts_0000 ts_0001 \n", "0 0.000000 0.000000 \n", "1 4.938129 0.000000 \n", "2 0.030601 0.046903 \n", "3 1.255303 0.000000 \n", "4 0.000000 0.000000 " ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t.as_dataframe.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We fit each TS distribution:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[Chi2TSD(10000 trials, eta=0.524, ndof=1.119, median=0.005 (from fit 0.007)),\n", " Chi2TSD(10000 trials, eta=0.500, ndof=1.016, median=0.000 (from fit 0.000))]" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new_bgs = [cy.dists.Chi2TSD(t['ts_{:04d}'.format(i)]) for i in range(len(src))]\n", "new_bgs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally, we patch the `mlog10p` values:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "for i in range(len(src)):\n", " suffix = '_{:04d}'.format(i)\n", " t['mlog10p'+suffix] = -np.log10(new_bgs[i].sf(t['ts'+suffix]))" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
gamma_0000gamma_0001mlog10p_0000mlog10p_0001ns_0000ns_0001ts_0000ts_0001
03.2500003.2500000.2807520.3005960.0000000.0000000.0000000.000000
14.0000004.0000001.7413440.3005965.1841950.0000004.9381290.000000
23.2362522.8424870.3293300.3715630.3320280.9210090.0306010.046903
32.0434734.0000000.7956450.3005962.6146850.0000001.2553030.000000
42.5000002.0000000.2807520.3005960.0000000.0000000.0000000.000000
\n", "
" ], "text/plain": [ " gamma_0000 gamma_0001 mlog10p_0000 mlog10p_0001 ns_0000 ns_0001 \\\n", "0 3.250000 3.250000 0.280752 0.300596 0.000000 0.000000 \n", "1 4.000000 4.000000 1.741344 0.300596 5.184195 0.000000 \n", "2 3.236252 2.842487 0.329330 0.371563 0.332028 0.921009 \n", "3 2.043473 4.000000 0.795645 0.300596 2.614685 0.000000 \n", "4 2.500000 2.000000 0.280752 0.300596 0.000000 0.000000 \n", "\n", " ts_0000 ts_0001 \n", "0 0.000000 0.000000 \n", "1 4.938129 0.000000 \n", "2 0.030601 0.046903 \n", "3 1.255303 0.000000 \n", "4 0.000000 0.000000 " ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t.as_dataframe.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Remarks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we've demonstrated how to perform an analysis that includes multiple hypothesis tests. In general, one can either:\n", "\n", "* Characterize the background for each hypothesis test, independently:\n", "* Run more trials using a `MultiTrialRunner` so that correlations between hypothesis tests are properly accounted for.\n", "\n", "Or:\n", "\n", "* Run trials with a `MultiTrialRunner` and, if needed, compute per-test p-values from distributions fitted to those trials.\n", "\n", "Either way, `MultiTrialRunner` provides a mechanism for efficiently performing trials that account for possible partial correlations between multiple hypothesis tests." ] } ], "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.7.5" } }, "nbformat": 4, "nbformat_minor": 4 }