{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", "**This is a fixed-text formatted version of a Jupyter notebook.**\n", "\n", "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.8?urlpath=lab/tree/image_fitting_with_sherpa.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n", "- **Source files:**\n", "[image_fitting_with_sherpa.ipynb](../_static/notebooks/image_fitting_with_sherpa.ipynb) |\n", "[image_fitting_with_sherpa.py](../_static/notebooks/image_fitting_with_sherpa.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting 2D images with Sherpa\n", "\n", "## Introduction\n", "\n", "Sherpa is the X-ray satellite Chandra modeling and fitting application. It enables the user to construct complex models from simple definitions and fit those models to data, using a variety of statistics and optimization methods. \n", "The issues of constraining the source position and morphology are common in X- and Gamma-ray astronomy. \n", "This notebook will show you how to apply Sherpa to CTA data.\n", "\n", "Here we will set up Sherpa to fit the counts map and loading the ancillary images for subsequent use. A relevant test statistic for data with Poisson fluctuations is the one proposed by Cash (1979). The simplex (or Nelder-Mead) fitting algorithm is a good compromise between efficiency and robustness. The source fit is best performed in pixel coordinates." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read sky images\n", "The sky image that are loaded here have been prepared in a separated notebook. Here we start from those fits file and focus on the source fitting aspect.\n", "\n", "The info needed for sherpa are:\n", "- Count map\n", "- Background map\n", "- Exposure map\n", "- PSF map\n", "\n", "For info, the fits file are written in the following way in the Sky map generation notebook:\n", "\n", "```\n", "images['counts'] .write(\"G300-0_test_counts.fits\", clobber=True)\n", "images['exposure'] .write(\"G300-0_test_exposure.fits\", clobber=True)\n", "images['background'].write(\"G300-0_test_background.fits\", clobber=True)\n", "\n", "##As psf is an array of quantities we cannot use the images['psf'].write() function\n", "##all the other arrays do not have quantities. \n", "fits.writeto(\"G300-0_test_psf.fits\",images['psf'].data.value,overwrite=True)\n", "```\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "WARNING: imaging routines will not be available, \n", "failed to import sherpa.image.ds9_backend due to \n", "'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'\n", "WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available\n" ] } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from astropy.coordinates import SkyCoord\n", "import astropy.units as u\n", "from astropy.io import fits\n", "from astropy.wcs import WCS\n", "from gammapy.maps import Map, WcsNDMap, WcsGeom\n", "import os\n", "\n", "# Warnings about XSPEC or DS9 can be ignored here\n", "import sherpa.astro.ui as sh" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Read the fits file to load them in a sherpa model\n", "filecounts = os.environ[\"GAMMAPY_DATA\"] + \"/sherpaCTA/G300-0_test_counts.fits\"\n", "hdr = fits.getheader(filecounts)\n", "wcs = WCS(hdr)\n", "\n", "sh.set_stat(\"cash\")\n", "sh.set_method(\"simplex\")\n", "sh.load_image(filecounts)\n", "sh.set_coord(\"logical\")\n", "\n", "fileexp = os.environ[\"GAMMAPY_DATA\"] + \"/sherpaCTA/G300-0_test_exposure.fits\"\n", "filebkg = os.environ[\"GAMMAPY_DATA\"] + \"/sherpaCTA/G300-0_test_background.fits\"\n", "filepsf = os.environ[\"GAMMAPY_DATA\"] + \"/sherpaCTA/G300-0_test_psf.fits\"\n", "sh.load_table_model(\"expo\", fileexp)\n", "sh.load_table_model(\"bkg\", filebkg)\n", "sh.load_psf(\"psf\", filepsf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In principle one might first want to fit the background amplitude. However the background estimation method already yields the correct normalization, so we freeze the background amplitude to unity instead of adjusting it. The (smoothed) residuals from this background model are then computed and shown." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAEJCAYAAACnhI2ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvWuwLdtVHvbNXr3W3mefc8+9Eg8jJGxwWTEBihgIiMIpFwY7IKAsykAsYwPCKnBSUkggVQZMymAeBdhgx2VUkPCSEHHEMwZjESqG2C6CEcgyZYqYBBXhcS0FGV3de+4++7FWr5750T26vzl6jNlz7fPQlnZ/Vefsufoxe3b37DHH/MZjhhgjFixYsGDB9UH13m7AggULFixIsQjmBQsWLLhmWATzggULFlwzLIJ5wYIFC64ZFsG8YMGCBdcMi2BesGDBgmuGRTAvAACEEN743m7DdcfyjPJYns/DQ1j8mBcAQAjhfozx9nu7HdcZyzPKY3k+Dw+LxrxgwYIF1ww3RmMOIdyMG12wYMG1RowxzB1TP46GXBd8A5Wrfq4QKmBVjduqanqMRtt2/6S878uxBRrZDqDp/wLADsAe6T70f2V7VL9lm0YAsJI2Ypz21PRPfm/68grAUV8+PgY2/Y71BljX4++67p6Jxr4Bmr7R2y1wcdGV758BZw1wvz/uDMAplc/68rb/J/cvbZd2rvvyEbV5Q9vXdLzu1fLcgO7ZyTNr6Jg9umcgdVQYn9Ma47PZ0PXXdP26mj6XSH3Ae8/cnhbjcYD/nlv6q++noX07qkfKUV0DSJ+b3M9KleX+pf/IM/aeu/Us5XmuauCof4ibzdi3Npuuf8kxK3qe+3Z8nttt90+XL7fArn8A8pzlXr3n7H2DW6TPk787Xa9Ani2/H1C5dcqCz3nVq/DNr389SrBQGQsWLFhwzXCjNGbW1ETjrZCOikX1aJVkZnsO1vVFG95jqh1aKDmGwRo/a33yd2603h9wn6xV6Xor9bcUehbRqn2etqOfs/yuMdVSJ9tnnkvpI5E6WtjvOWJ6HdHqKtjvulJ16fPlHNZydVm3j7cxvO3A+Axyx3hYVeNs82HC+zYC0ufJ32GFsvfpHef1k0Oey40UzHqjPMhSAc0CmKdgDKsu/ig9lHQKPY3kc7kj6qmnB03NzB1jHcft4M7O02DpbDLV43tY0XFWx+bL6WfE1BBP8a1pvQj1Gun7kHO4zUyFAEDlCGePOhDha33w3sfLU3KB9R4jlRtVnzfQa0HE5Vptr4zjLAHuoURJYaqQj9eUons+0mfP98OQ57kyyld5nvx+9fent+t+W4qFyliwYMGCa4YbpTG7GjHRGjKdClWnDa+MoYun8pa2XAJLC+Ipaw5aE2WNqOR8YGq8rGY05oaMf9GgQrhtrP1KB2uQaoKW0RL9uZ7mM7QdKeXAmvEOvrEtuR8qB7Xdu75+1pYhT9+j1uhK3g1r2Ayum+uaq5fpCc/4pw2EVv+y2sRtln3cJ65C7zG05uzRYdIGvb2FP0vwnluOGvRmPSWzoUOE7Y0SzNKxPQHNPOIcvWF1vpwgWCEVvNb0NqDsBevppf6ILB7RQ9sLZbmHphk/hqqyvU/0x6an3nLNNVK+ViDbVrRvThgDU0HIApjpC8tDQzA3LecB41KdW1F9pbx+TpDNHQOkA0hNv3lavoZNxQgsD4sVxmewVmWvf10Vlk3CUngsb6BDwDy7ICfgWqTUlrZVaEh/yn2XgN83DuGYFypjwYIFC64ZbpTGbGmie0w1OIC0Z6E5nCEsZwTT2rZHX+hzvNHSOp+pA9FuSrRPDUur0ca+xG/X0Z5ZY+O2cUfTnU7712pqBpg+S4++0H7hnmGO6+Pnyloyt1MMeVbbchqUNhB54D7I03DWkrmtbLBrMD6/nJcAt599x5nWkHqvorElBtmZ/rRvywx8gH0cz0CtfcN1YAs51pblN59nab1syLWuZVGJXM8a5bhRglngWUqznhTKasxldr3zBAHD4wpz7QDsoBLNe2rhkXhIZAYXnlpaH9VeCeMcdxjoL1u9LUs1txsomzprR3/tIpfjmD1+kKEpDLlGjVSASps967xVt7VdW/E1x2wJYAa3yxoMmMrI0RdSl0WHcFstJELOGdAndc0oPVZ9HnJVMJ028OBSd/93hbQ/WQLXGmS979B7/qVYqIwFCxYsuGa4URqzpSnzKMhBCCvneAuVVn/oHDbKCFhL1oab3LW8UViP2jX95imqGFcsDUW0mj2mobK8H7A1GEtLZgOV1K3vQ5/DhsAcWLP0tGc+5iq4RNoH9HTXMgR6hrcc+DmxtqXbz7QGa1/6HA2mMjzKiP2YD9XW5Bvi92sZinmmVdGH53n4eNcC8u/VM+a2GJ+f3O9c+gM+xjLuWf1WxxlcRWO+UYLZE7oCtrpLboVkP71xz4LML6RtU+6PPyq5vuYHLU8CIJ+nQAtp0wpPuQmsnCBmfgz1kEqmlMzD6c6VozKsKeGkPeq3JfS5fuujC7SPB+Qd5qHbxR+lNeiUupoBaZstHlv/tvqNdw0Wujxock4Vqz8JPIHFtNBk0HLsEZpnBrp+uXcE+Fyfs+5Xt9+jABmX8IWh7re5bw/wn6fkIynBQmUsWLBgwTXDjdKYBaxJJZ4YKBupOCMdADe01Avc0NNT1nYsA6Ls47KX/0Br0yvJ5lWNmb1YY5ZscpZfqdZkPLBmyJooUzaeVqnPt/YDvtHpQcB1Woa80muyJpozEHkGWO0/z9oza1872P7gGnq2p/sHMPUdvwp9wWWebbawtWQ2IDfN2B/3SIOWkuvQ9jmqw/Na4tNyGf3mqAw9Y9Jasksh0vGluFGC2QsOsCz18nv4YGi6vyLBNolMYoGtOqVwGR4HtkPeK6GEw+KPP1RdSk9plyWYRSjPWcb3sO+5qoBVmwpdryPmqAyrzMcyhynw+FuuQ87xprcCK+2mDgLRQtsq62tU6FKGAildpAdCGQC5z8g1uR1CueRoopJBPEcZMfRz8YQc1Ha5zr4dI0br2qE1lPLSNFNqA+mhWejvw+vaGyrv1W+BzvHCdfGg6fH3mnsuxY0SzIKIlGu0Htgw4jpCzhNSWWMHqCydEVODzlU1ZskZLB99XaftFMFs5aD2NDlGNOqSWQF/8HOagaUxzx3LHKZ8+DwY5LKrzbWBnydDC6WSurkuEcpFA7ohiJrG51hZe64wCuxsXy5oP+Bn79ODhMfxs5bKmnDT2LYM+W6SKNNmPGdoBz0LK6LTuk/9nGsa9IBUAIqb5AZ2LIAWzFoA10aZ27K4yy1YsGDB+zBulMbM02KP1tCuT562w7TAcG5FARXtVAuQYxtaOsHinq3RMueVwTyy5o+lrL0yLFpDo1XajpzD97KqOk2a84VYYC1qbkYwXB/5hD6elqw5+zl40/qcG1zu2CTAp2CmxTA9F5wQN3azLJlteTMDPs+im3JeLvw7cZeTcx1NOGlfNe1rZo4W+PRJUh+mVJ1goEia7plZvPIlxudgrYIj0Cu95AJ2gIXKmEWJkY9fLuALvLpOjxmu0U6Nf8y1CSzuuYTKyFEs2rAHdMdY7fdc5YCU7wTG9q+o/XKfMkV0E0TRvrkOWuLTy9xf0mbYzzF37RKayPrN010rcnFVd89d3jc/dyB97pZ7WVWlU3mWEnpg4vbneFhLIbEMntazs3h4vq5Xn9yDNxiJkTka30qSCsA+3WyDpVAl6J3CPQEojzogTXQ0t6ybnMPGv0GJOWBFi4XKWLBgwYJrhhulMZdOoRlWUEYJFVBVmHplSDu4IQ2Gt5DLiaw1W6tdopF5mr2lSWv6g6/P00gPtdKqj5rRiFIhTcGpkwJ5KLH86+MszZpfS66jezSR/Obts1oyRi8M0ZatBUhzVEajaK7BXawa+wi7IjKVM/dcLcrHCriaWyiWy6Kle+2R9u8opWxLM0WmwoCUytir/uW5tCXugc5MdwKdIYoQ6BDdz3KuqnKOudjzAVzGjRLM3nTW9a9VwpgftrVdCzc9FY2qwwC9wFSC0PMptlbw1tfXg4ZV1lRGTjBXM4LZ8i09Eos6/A6eEwSaWtBl/XgsX3DrGgxdh0VRWG6J7BLFyf0T6zw957oeqR+mNQD/mQskN7b2gDkUWoBKWWANfpYA9N5fjhpkIV216fdg5f0GUsEcW5vjnqMjrW+IrykYPI3aVBh6fYj7V84TxKJSwkJlLFiwYMH7Lm6UxuytYFI6w2BPBN7mURnsvQBQEiHY56Bwu6UlS/3ae8SiL+aoGJ5uzuUqmORjroBGrsPWdeS1ZIvm0LQE0yLA+N44p7VnZGTfdRj7Lc1c59KtqJ1s7FmDNGll7Nts/JkKI5cbQratqvTeDk3QZBn5gHnt7BCjmz6HtfSEGmumz4KfAdM5uesf4jnDmrmU5ZuoMS4rp9O7ahyiJSfXP0BjvlGCmTtLzjWtBHMC1PrARKA3rX9+7jpeUIv2qrAEM1+/mINTda6pt8j91UYP8qao3lqJOrKLp65SvbUWH/9l1zGO1Jq0jcqWi5uuV9MaXj5jeTZHm04YA6OQZsHMHy+DI/8EIjjM96luzhrM5lBC+cydW4LkW9Pp3TD9VpLc37Td45c1cktU8ffJ/X6P8dmzwJ2cr+qy3qf1zR6KhcpYsGDBgmuGG6UxW9DTOx6ptGGLLcWWEUZriPrcOQ8HjxLR+zxf2JKcFx4sn1rvmqsKiE7PCdV4vs6NINDPSCe6SXx6pV6M70a0V6Y2+B1aXh2TFK7ArPbMx+hc10xlrOpOUwZSY9/RZqoxe1SG9KedMcW/CvQqOQLLE8O6jHfpuVmIdbyV4Cm2vmavv7vcykIPE0l/7/96dCTgG+QfBm6UYL4Kr2YJFi8TVnJ+OyNwMufP0Re5ZDhWO/T53Nk4YU4Omp8buDkVAanzO3j0BV+byzoCbJjiN6MgFk6ZqQ2mP/hedAf3Vsm2eOWcV0ZdTb0vAGC9GWmN2qAyPA+e4ZpVmpNkVaW5oof+RNu8pbV4v763HHJ9wdvn0UIaLRVc4U8DMncZDvYobZeu1wPz93q7Lj+MgXMOC5WxYMGCBdcMN0pjtnxl904ZSOkH9sGsazuVoaUF8dJMZsrDDLRhrlRT9iz7bpCLuqZ3DGvZcv09ps8gR1t47eTntTeec1N1WjPQvQvWIlcYtenk/el7U2Uv9FrTF7KtAsz81nXdacpA6qss3i+eMbYEuXd2FaPdHKyZQwksqoj/WrC+G30+Y2+0y9LOo+YkCzH06Xa6DZi+v9y7LPm+c3isgjmE8AMAPgfAu2KMH9Nv+1AAPwzgeQB/JcZ4GkI4AvBDAD4BwLsB/KUY4+/0x38tgFeje09fEWP8uX7778QYP7y0LfLccssYVUAS/MFTbM+6m1xDCWPezuWSl52zNOtrutvZIm7k7Zg73/M0SYIIMjx7iYtY0wBt3ybNtw5ugY1NbQBTFztGyfJck1y+cs0q5RRrJYCtBFErRfOUQK9Gzn/l/oByDwWNORdBoJyWsM7NbSuF14cOcREsWeDBg1Z6NHVx6OA6vMcDXtrj1phfD+C70AldwVcA+K8B/HEAfxXA96ATvO+JMf6JEMIrAXw7gL8UQvgoAK8E8NEAPhTAPwsh/Ecxxiuvucn8XGWUWRgLtB+xBy3MskLTqUu2xTYVzod2PB5MgGkU4tzgMKs9k7axcXpVbgCT+2EtmX2iTdA7KYn8K9GYLWEsbfH8vzWPnIM3m5F7jlRumm72YIX1A76Q8rRdj1fX55SuU5jb5xlTk/Od2WBkZeVRTAuk6sJZa+73XL2JbeoAwfxYOeYY478E8IzaLAbiFuP7fAWAN/TlHwfw6SGE0G9/U4zxMsb4/wJ4O4BP6o/7D4+y7QsWLFjwuHAdOObvAvBGAM8B+MJ+24sB/D4AxBibEMJzAD6g3/7LdO7T/TbEGD+x9ILeShR6GpxEppEm462Rl1yDuOWhfmOEtgJEeB9r01KXPt5bC21P7WBOVGvPWjM8hB/LaRSlgSya8uEIxSIozTnHb15FS+Yy/7ZycltIco9U9nbWkqUstIaVN6KFHfmYgzdjAPJa8lzdOT7a4+9z/Z63V8Bw00wd88w2h9KZpccra8xRl7p8Va75vS6YY4y/C+DPqM3eSj9zKwDlr0WV84tlYcwLSmpU9FXk8g5bZS8MvGrHemOVToU9+iDndqYbI+fsgTHvs0FdCLVRYmDMUS5czglm9nfme0jyD2PM8wxMDWHSvlghWTJIDjPdnzAeZwkJL6RdC2Yv2pIhvrqecLDyDzdNx60DwG6b/m6cfsc+3qX0AjCfm9r6PSeErfoSwWwMgLwdUO6WfH6hkGvb8ZxQ2UI3J7D1Pm1PmqMkuVxicLfwXhfMDp4G8GEAng4h1ACeREeByHbBSwC8w6skhHD/UTZywYIFCw4Fy6UY423rmOsqmH8awJcA+FcAPh/AL8QYYwjhpwH8oxDC30Nn/HspgF/xKuGbDiGYmjVPCXmqpDUTL5LKO6ZVfzVYS5eXUPWeE0w5WBq0591gIaFCSFtZK1rD0hInWr7SHASeZqw1zIkGTtqL7NL3KxGGWkPWdZdqVVehL3TZgmfUK6VvLCqjaYDLrbMck9+ULErPy9FAAs9jQ1MXlpHvkKg5/layfd04jqvV2jNDz0KlnaVUSElQlZQ9Ycx43O5y/wuATwXwgSGEpwF8fYzx+41Dvx/AG0MIb0enKb8SAGKMvxFC+FEA/xe6iflrHsQjQ2AJVKY4eHsOOurKO4f5Pe7cQU4y/HhZMJRMpYZrWbQDxunxOtMDvE69qqaDhEdlzC7x0x8n7eEPOUeF8O9VNWYG89qf2/YwcKjvNnuisCcGC+mdeGXIuaquHH1RojyY7cvs87hE3R4uay7ZGwCTehy6QNMamjasqJzsMNo3Bxao0rdygp3PsdwbJ+2awWMVzDHGv1x43AWAL3D2fQuAb3mY7VqwYMGC64TrSmU8EsxFM2nVu4K78swAPSLmVhLmdqzoGDlH8jAMlTajtrF2IgyBq1l+B+2i8NwK6fTuEMOgd5zOEQHYmu+hMDUnKI1OafMl0NSEpT2Zmp4xnWXNmGmNXTNN4mQZM72VvS0tes4gPbdaCdfNBvRsgiiiiQQ5w2oJzOc902b34KuggEaxfu/V3xLcKMFcgtI15iwB3GaOAVKPAP4Y5CXsoNDSPhohSgVJDhUJ2ZLjAN9dLyfc5dAkHF4JZZ7W56Ilvd8PEuXFbcpBW+T58KYdn4nlMWOF4utI0mE7ruZVIrC4fDMtAO/HVDjPPVKm4/bw+WagjL/PweTvMVWEHhYO6U4ly3J5VFQOj4pyW7BgwYIFV8SiMSM1+LEW4OaLRToKetqzdR7gT0NlwVIZLXdUEed/jjPDaS5gRSBayMOYxpcEpbAxVagQrmNvTOv3arq/Vxqmm3/Eub609RAwfcTPYHK/Dq2htWSZ+bD22iBPK3hBMVZCJev+rOtX7Ugb6X7uaaJaM5ZLzVGEc14uV6XmcsZ2qZJXvinFHB3pNc/bvlAZVwALQiBvtbaEsRbSgrlpaAn0cRyUkjuezxtcjZTwTIJVnAbp461pvARQWAKeBZn2qtir4yzXMc4TIXkj+HgzutK+lQlKojete7Geoey3ykxTcOa8FiN1xX0IVGa6Qv/mtQR1FKLOO2EF7+yasR9xdj7p2953wEIm15+9R+v1Bw448p4r/9bfnDdO8jObE9KekPeopZLINqueEixUxoIFCxZcM9wojXnOoqynQt5oxxqFNtboUdfSXkHbdT5ohg4dnwtykamuaBW10jy8afzemaKXUBQS9s2+nkN7lFaUoxE4P4RolU0DbLdjWYcts/bpNVN7YpjHFLZRY+7ZaOMk57rYwTcaa7BmKvSFXjVFtH8rjevOcC1q23KqjsGzMt2frUCQklmJ9BnPgGpBvrWcodQCN0eviOLNfDU9aTWv5PoLlTEDdpS3XOSA6cdSIoznpi2WMNYDAIPSJk94aT3FBcap2hDJ2NJxTEWoshzLv4fzDYFdGsiydwSelYRpLgpuu80HXvB7zHkIXCXy7GHBDQzKnJPkwajSZaskOf/RZmy/JZhl3yVtTxJZUQP2mPZvD1qhsFYnr9p8vzsEPLDJdSxK0WqbhlZ0ct+03v44KIkbJZgt7VNrz/qhe1yy1yFyXJe33dM8oLZ7hiHtQpUkiietTSdjGcrONWX/HPQxZrRhRghqob8jjtlK7rNv8gYzz78WmK7bJuU5QWFFp5WCQ8nlfeRWXeFBe4Vx5lPXHbcMdEKZF4AVgSzcOa+8M9TdpgvlsiEwB36+Xt/0VgJiW0DlzcxUfZ5bITB+N/Ldzc06uG9b3zfXy/2JF/rV1yjRfDnNg+AQjXnhmBcsWLDgmuFGaswMjwMW5NzfrPq9oBJrmiXb9W++tkW5aEpDu+FxwIpEE3pTSvktmlRJmkqGZzWXehlzdIacz14ECa0hx6DTaEp4Zb52UFoylx9UY2Z+lLVC/XyGZ6tySMs7rJC6x9UYtWH2xFgrjnmt3OWkD7Anhn62vG4jXzM3XZd2Wv2cZ4GaztDX57bm7B8MrTkL7QJM2+vRGtb3KPvk3ng2xtfI9TlGhZSCtNo0hxslmC1XHy0YPcqByx4HlUPJ1GRuqsNcskBTGQzxiwamU0qGxzfKebmyFaVntt0QinPuZhye3LRT97K5gVbDCwnWbcm138tX7d0LU0hVNQrKVYVhcVkJvdao+8GkNgx+ifGPFnyVbYMBuKb3ru6fsxiuaGBmXlu2zYFTC2hqjI2g3Pc4QRfU8dyvEn/3TBtK7Aya5rIoSY5hECVAyjkqw7IBzSV+8rBQGQsWLFhwzXAjNWZvqpFD6TQk5xEgo6C+pmcI1JFWHh0i2rt+mdoQWGW0W2+bpRHvW0wS7Vht1G1JZibVlM4QJIl+5JpItRietejnac2AQjXVGOdyNej2eUvZe5r3EPlHq37XRM3wVN6ijySij1fdtjRmHWCiITQHpxdd1aPGXtdA07sl5gxkpeBZaOKB1Jalm7UoIKmXcWiAB9NE0jbrWNaeG0wNgd536BlGr4IbJZiZj7IE26HZ5wQOC5DsK4XubFZ2O88Pu0H3Qrk9XGZ+MZdZzc2GJvvhUzkexw5qSwUgtOP0mfdJnbnrNOq4kmcsNIa1ynVJFKAWwp6Qz7kC1sSZs5DWAyAL17oefZfXSkgfmhCIaRFehKFtgbV6uFLNCr73S1K3+i3VsaBn980c9aWfh+fXXKJUaYVAoO+JKQqvfFUK7SpYqIwFCxYsuGa4URqzNpoBU41OUxw5akLOyWluvK20Tgul0yfPn5SPq2nY1wmRmJpgbYGt0xatIPCiE3UCHsv/WoO1k7kILE/D8LwvVrXt02z95m1sZJtbwFXTQKwls/bs5f2Q6/A1B425UEvm/atqDPpho2CtqIOqSWdHHDhyqKFbGwWlAmvZtKsEnRxySk7DbpH2L69cev9MWx6aIwe4wYLZsppawsMLVpBzIvIPvHSJeMbcy2eKQnd85sB40PBQKY44Z5FmrwiL1rDgPdtatdsatLTA1wOGt1bjsIYinS8Ck4WU56HB5+iyl9EtV1fDXK5DZbCQ1m2wEhQd6lWij+MoQh4MqgpoqpFm0jmceaAUlC5T36ofuTzaOpxdwAtMMM3iDsyZ9nhugVoJGdqUqYvbpvtzpfaXYKEyFixYsOCa4cZqzDKirWFP/S1YI96cj+eclpzb9jCs47mgGT7Ooik40Q5bp7X27LU1pzmIoRL9MRIIk4TzztyLFWo7mdkoDZN/W54Mc0ExHMjBHhK5RWfluH1rUxk75aGhr28ZKXWbczkoWEuWduqcGvq5WCHRmnKxNGmgPPSYg5p4gV8PPAO0fK29b7fESMhasoa3nSkLNrpzf+Zv4BDcWMFsRsch717D27mce3G533Pw6ta8tRZGmtrgbcD0I+Ip3Q6pYPaEtJcrRLdTR7EBI40hdaypbn0/JVNm/ij1QKiXNeIkQB79kAsiqY3zteualwuEo+00/dH0ZXZFtNoyR7/w9Sxw3g4WzsPAUPsCWO5BrqETTw2DA8roONDxSXBI5kPJKVGH5KLQCkWOjsnRJIHK0vdY2WOBfQiVcSMFMwtbjo4DppyVRdxfReBeZdT06vZ4aw09yDTGPtFEPWFslZl71tfRbZYOprlrbndEygsfsjioBX4+Wlu0It88dznLcMjuapxqs0Qwc+Sfh6r1BWtOGOtreecLdERgLpLT0p6TpP9tGj6vFwQoQWz9xX4HbZ4qFH655NtzE3Q5vz0hra/FmvFKlaV/10iFdykWjnnBggULrhlulMbM3gpWXL81IjIqZ7uFEjc7r30l8DR4bz1BjwfcI68l7+i4Eo1Zt4WfuXQ20ZCtiEWv/Tl4x2u6gAM0tJbL53jlqrK9Mqy0mwxOXGTxvSWa8CHH5c63fltatuUxoZepYlpjux3r41kBzwAO6d+zq7fLceq3wDtd52AusRF5QVOavliojCvCEiDsamWhor+WIasUJccfUqe3gomGJYzZ3a1Bati7RBmVoQ1zAjaIMC1R2tE0leENgvpe+Zqyj6kLiZRjXpjLQz2O8BJB6tEXnvFPZ1Pj+pjvtQRuiUGvBCyAS90CLTCVsWuATV8WoSyrzfB9Ng2S5Fkl+b1zGNpHhkMgNQaWCmU4v7Vh0aMutdBdqIwFCxYseD/GjdKYD4FoXkxfsPcGaPuhuErkn76ut93TZDmXrKUJb+m4rXEcn681Zq8tNQ4f+XORUvw7Z/gZtBg21lVp3omcV4alNYq2bBn5dBShp3U26qGJJvmgWmQpSiMc57RmIHX3CxVQEZVxuVUniVEQpNU6VInUVwqdvlNDB0Bpo7dn6POidK1AKSlbVAb354dOZYQQ/hiAl8YY/1kI4RaAOsb4/AHX4bo+E8A/QNfO74sxflsI4aMBfB+A/wfAl8YY2xDCNwF4Bbpn+C4Ar4oxviOEEPrzPwvAWb/9bSGEDwfw+hjjpx7cJiprdyvL3zm3/l4pSs7RQr/kxerOtqMy+yrvaP9W/dY0h5S5U3tgKqPB4c+GqQyPb9Z1rpBOF3kpJu3S5rnLWcLAyiY3J8wtocYC2Mtf7S2lxPuuwi97tIpPCy4AAAAgAElEQVT2Y7Z8pa1rSh2revTJrrb2MRM085F6QOqhIb913aIcWG6bejuQngPYgrxE8fFkgqY1rHQDD5XKCCF8GYAfB/A/9pteAuAfH3ANrmsF4HUAXg7gowD85RDCRwH4KgB/AcBbAfzn/eF/N8b4sTHGPwXgZwD8rX77ywG8tP/35QC++yptWbBgwYLrihKN+TUAPgnAWwAgxvhbIYQPvuL1PgnA22OMvw0AIYQ3odOKeUWb0F/nHp13G+Pg9woAPxRjjAB+OYTwVAjhRegGwWdKG2I5hvPoFox9Zg4Go84SuJ4ED1C35O0Q7VdTHKL9srfFtv/X0G9LY2ZtI+dTzMZU7lxaW6lgay18jneMpYEkxhbDQCfl0mi9yTUNjTmXOnUOnM4yt/ioXMs6v9Qo6GncVlDNyrhPsy5q42Z6yQHsycF1aVrHOk8jmyoUaf/yVq73ApZyxj5t9PfoTS8vzDAbOUBAlAjmyxjjNvS1hhDY0+lQvBjA79PvpwG8DB018U8B/BaA75GdIYRvAfDFAJ4D8Gczdbw4xvhWAH8xd/G5qQZbY9f9b45W44ftcVAl36j34nV9lTp2LtKrRUpdbKm+HVIHehbSzCvnBDOMco4H1uew65z+xpj+0MIZyNM6/G5qpHQFr4unBXNJFF3SRkdgzUXhWQJ47whjL2m+1Zbc75J9VoTjHGXjUTGcoY4T8uugHjkmVDZFkYO2bXhZDj1aQ/Z5sL5DLR+sRGcwtguu6t5Yctq/CCH8TQC3Qgh/HsCPAfgnV7ucqfjFGOO/iTG+LMb4V2OMe9rxdTHGDwPwPwN4ba6OK7ZnwYIFC64dSjTmrwHwagC/DuCvA3gzOkPdVfA0gA+j3y8B8I6C8/4ROo366w+pI4Rwn3+zddTShNdI6QrRmuU3n+8R+of63eptenpZYu1n39Kq6XyRge5erEAQNuSJR4bllcH0B59fCt1crcXIfp1vw6IvrOfHSZCGd6O04k0/z173ZfmtNeYS8DsoWfUkt4SXDtbgsoZJ5xRowaWwcnB4z0aHdQP9rKBOtWReKWVv9OF9Rlv2+pmlFXv0hbXdQqmWzNuvrCX3J7JcijHetg6dFcwxxhbA9/b/HhS/CuClIYSPAPDvAbwSwBdaB4YQXhpj/K3+518A8Jt9+acBvLbnp18G4LkY4zudtg83HUKIzAexYJYyP+gaU6dxj08CbdcwtxncHmDzeyXTbV5tuKkwSFadg4LBXhgsgFlo5zwxKtp+xdmaG4k5l6tAKCceUOXdrOtU+LJ7nBbalvA5xCPCozW8abnHJet8zJOglL68p2vNJZbnQaOEGuE2evXyveky/11VZOcopIk8qgywqYireGVosFD2vK44QERzzEDZvU2WG3OEMcMVzCGEX0fmvmKMHzvfpMk5TQjhtQB+Dt19/kCM8Tecw78thPAn0T3n3wXwX/bb34zOVe7t6NzlvvTQdixYsGDBdUZOY/6c/u9r+r9v7P/+FXQC8UqIMb4ZnXCdO+7znO2R2nQQxHrMox5rzBzvLhrzXPy7ph8EOYt9znBU0TTQozK0hsdUxo58SpttqjFbg7sYAllLnvPEyM0MJhrFzHmCOS8Vnc9ZO+6z77I8v80mNQSuHYMfG6L0s9XvxnuvXoDK3KKjgG8I1KhUfZ42r9tWoj1bmvBw3So9RpeB1MhntakEHp2VW/j3KvQFYBvyrGARpj25v5QikvH3EEuYK5hjjL8LACGEPx1j/NO062tCCP8ngG8sv8z1AHOSTEuIwGa6QgS2F/9ucY25CLIcV6fr4qAILltWdM7fu91OP0KZIm8wcs/c8WSaZk0l9YdgwRLElpDWZR1VyeeHmXIN4AhEX2CkLzabNDWnXlXaym8BjGVvReaHDU8A6+08OOhUndYxVwGn2bxqFKL33HLfgBcsAuSj9fgYj2rjpd/0NmDsf3N2p9oYjL2B9xDKaA4lxr/bIYT/LMb4iwAQQvgUdH7F73NgjZmjd6wMUSKYeZ8VUZbjhHMJ1xme4SWX85cxCN8NcHGR7o+n/TFtJ8wA4AJTTfSq33UFuAJ4MrPoy5oj9kLfNccv7d/0dcv71Fnj2OBnLV4KXM1I5qGUY/aQrAaS0ZjbFoP/MB/H1zzEp5qFqTYwct1zMzXPhU5vK30u2hjMmjAfw6jouNzsi/uZnoUNSpkSxiVG4ra1lSXZBwDxYWjMhFcD+IEQwpP972cB/LXySyxYsGDBgkNQ4pXxrwH8JyGEuwBCjPG5R9+sRwOLY9YammhlQmPwKGrl32WKweMw5TfD46U1lcFTcV5vzppOXW47bXFVT/c1p2lqT4788zRbD5aGzPxcSVlTQxa/xzOYIzpm0+/b9I042qRUhpWoSMqe1uMFOzxInopDUKoxJ4EbdI60L7bjGnoaPNWe3C/vI42RaQ6PipDczBZnXgKhKEpXrgG671avUG+9Iou+kGM1pemtgD7UP0P3zM0YDsGsYA4h/C31GwAQY/zGq13yvYc1lS2hoMtMXzB3yeUSIQ2UfdiayrBCiqVsRVPtGuDyootys9rQCq2BNGucZ/DzVkGyBK4lgNnVSBtP5/bp8gZEXaAT1CyM2UXOWjBVQ3Oic1PuQUD3vz3XtbkP0eIh9TkeV9k04/3saf1ALcxXlV+HFYVotSehwwwKiJ+fCGUZNHLUDPtxe/B25ZJ4eQLa8kMGRqqRB53c8mBmO533ptdJ9HjwHEqoDA7SOEbnrfHvDrjGggULFiw4ACVUxnfy7xDCd6AL8nifgxcg4mlodZVqZaKJ8tR5VXe/AVur9UZe1yVKacwyvToiLX1NbmC6zouL7p8+hzXr/Vmaj5mt3dysCqnWbNEdWmNmg6n3bGUKOXecpzFvqlRL3qj3we9gaLvSaLWRzjJg6UAPxsox8nG9c8axueuLVpmscE7HWSuDVNV0lW0zbaZhvLNQ5FXRpiuVeLRGaYRfcv3M9jltGLCjGvl+Soz1WQpoBnzIIXkjSjRmjRMAf/wK573XMccxV3RMrT7+4+OxfHQ8CuPNptsnZYt6APKCQE87LU+M9Qa45VxHOlfTdHTG+dm0zQn98QyAntawsnDJsxH+WfbzR+JRGSVlLZg9KmOjjpHBSAZGy0XuaDPep/atBXzBW8LxJtDZlpxr6HcbDSGlj9fHJPQJ7fOy0wEp5zx3n6X37D1LabNQGbweoKY1ShYHqDD2xzlaQhQXpu28NRjn7D4WtCfMHmPjrEFYylehLxglHDNHAK4AfBCAb7ri9RYsWLBgwQxKNObPoXID4A9ijJ5d6FrDMv5pLwA29rE2fHxcVma6w/OeYA1l1xgaEml8R4bGztdkukO0FqEyTk+Bk5PxfuS4RFM4nQZ4yHPi5EZao9Ya7iEa8yHn5KicwV9ZRfix7zKDjVqcLMib1rO2q/OY1DUGnoe9GPQUeK9Wk/a0RzaKeVP/CnYf2iutTsPTmKVd1lTdo2m8Ork+LjOtkWj5QOKfbPklA1MtWUd4WvSeN6PM+ZvrsncvwGgM1c/IM2haNOEcSgTzN8cYv4g3hBDeqLe9L8Bas6+m3xanzMJQhNzJCXCrL98+GbezwJTOwbzVXi33DnSCea9e/NC22r7+0TGAkz7GZ7MeDgoA1m2L9UUX4/fE2X3cebadtG1Nwrz+Q2D1bDo4Sbz9BVLBLB1MO+azc77n4sYCdw2fymC3QM9dUJ7tnGcMC1IRyt5U3uMRBU3rL7vE14lKePG0fqeE1N4R2Dnu1+SVZ87h9lj3fEiAi/Vb6mKhzsJsR7SGJMjn/mRxr5YNSPcBTe/xt8cUhxeYxffK70mXGduCZbT0oHMVFLAs+Gj+0SfK/4QrXm/BggULFszA1ZhDCF8LQBLkyzJPAZ0S9T89hrY9dHj5lHl05enxreNeO0WnrYrGeudO90+Xb510GjQAhE0/VFf9Vdv9OLxut0OyocvtmHiIp31AqjGfnACrE/px54m+fAvY9GExohJs+6wYZ+d44u5zfTufTegPHapc94ty1RgNoBf9P2Dq06w1Gq0NA2PotJRZky7RjNm3VPuI6ynqSt0PPw6BNf0GptqepzFW1agR68CLXNiyaIy7ra89e3THpA3O9FvDozRMWkFpjjnDpEB7J3j3kOSahr9IKpBqiUkyMWX0BcYZoDVb1VQfGwK5/VozFm14ux3pQNaQmQrTyAXVXGUVD1cwxxi/FcC3hhC+Ncb4tVeo+9rhiMrM41oCQjrBbXrxIoCfuAs8ebcr373b/QOAcOd2JyiBTlhumNUG0PTd8uIc620Xe7e+OEe86N7+djt+xEA6UIST43Q0uPvkWD7u70wOlt52cQncu9e37Qn8keM/6O9lm3hrJBnZngWO+854inTJKe532quFhS4nhRryW6ggEM+Dhd+H5/SfW8l5Lp8Bl0XQ7JWQzHoMkADmBPBe7g2uu2m6gRjohDTTWRYn613buhcul3gb6ET1WqiW0Bx6e/I85Rj4OZM9JIpTldIXx0pRkk+Cywltt0n7k0CEsryD84suOEvK1rOu6+6cuTwrV7lnjZzG/JExxt8E8GMhhI/X+2OMbzvgOgsWLFiwoBA5499XAfhyAN9p7IsAPu2RtOgRQjQ2INV8rMxkMjrzCD1ozHdGLfnJpwA89VT34+6To8YsbhRCZTQ70mQvgLPzvnyE0JePtpc44rlTXY80xcmt8aJ37wJPvWAs83yuCkAbx+uc9g7LpFk/efwfsKqfGy6x0vfdE1fH90Yq4xKpISOrMVNdQgVZQTmWxqxpgatow9Y2TQ+Ix4RsmzO+MXVxSNu0IUzKPDtiiuMQ49uDwrpPae8czTKhf4AklN/SmAF/Ws8LISfb6fvUGvPdu+ks9gmaUEq/u3VMuWPI+rtvUsri/tno/8/xB3z/TJ9pxBamlwk/l4cSYBJj/PK++PIY4wXvCyEcH3CNawPLOssPWwctaMEs8u/u3V4gA51QfuELxx3M/dZr8s/ajfOmi0tg0wvMs3V3HABcnAPbXTp3Ejrkzp2x7rtPjkL6qafGsvDNA5d9OQrmk1tEeRzhTtWxxOv62SS/xKoeO/XxMXDWd9atEh6aCmJu3nPx0wEhc5bzOfcs/s2Wf2+aD8xP0TmCTdelLezcPmt6K+eXWP7ZdU1H7s1RE9Yz856dBxa4XpuTwQxjIqxD1t9j6GARk2vOCOZbJ51ABrrPQOjFO3d6WhHo+vzQuVadrQfAqmlw6+ISt866jBM6mna452Z0gd1ue3tG38AdVB9Sz0Pu/1G5y/0SAE1lWNuuPUQw6A9Ku+BImTvCkXKXs4XkC0hInkzD8i56LfnsjL4g8iKuKqA6V4K5F6Zs6bhzJxXMorGHOwA2Y89Zb4HbZ2N7jm+N9fbXP6oqfGj1TDI4iQvS7ZNRrmvBzE3W4erstsQGGS28LcOeRs5AtaOPpzbcs9hi6bl66bpZMOYMYSUwudcZgafh5VfWnLueZcylAliRIdM61hxM2lQYS5P3/b85wazB6zzWSBdCYI6ZswXyQH/7ZNSSn7wLPPFUfxYrSMdHqXGclZaLy0HxuVPdG1rKnP92m9qgXGUBqQDWgxbwkARzCOFDALwYnVfGx2F8bnfRhWUvWLBgwYJHgJzG/BkAXgXgJQD+Hm1/Hp0b3fscLI1Zc8xHSmOWc24RrXF0Uo1cckIxkIsGj9RAR2VYDWjbkRPWpF9VjVRG4g90lF4/9NfEXXQJAHlhnV5jfuI4Jc8IawAfgmeGZyDa7wlpzBcXtsYslATTF6wlc1CMfrZzaRZZY+X1DMW7YUfuTcISTeppxu25AAkpc84Hza9qmkHXY23XnLVFueTSq3LwSLKPnpmVqMfTmEu26xnE0H6kya9YQ2ReeQc7B4u+N73SvDSBtWymDlbUN9ebqQurOYvV7qSDKnzZ2Xnoxm83z3bt344eGnP9dF/wbKR8SLBJjmN+A4A3hBA+L8b4EwfUeW3hCWaL69SRfzxd7+boYuQjIXl8NPK4J7dTd7lm09EUQO+n0092NpdAvRrrbdvxS9X+YsNxJPTXxxid0m6hW/WLX6uYA+puNwB88JQXkJZ+cPVMMhiJceXiYnT1Ygi/vKZz2BaZRCsOHEfvSijcuuNwvG52w3M62l4O0nd30SZZ9Hg5Ld1Ga5FVC1Z49CRqLeNWl/P3ZcHmTfcFLLA02C1PG0lztNBAhzTjcez2NfdspH0sfPcYp+tSjnRcLqoPSEOtJeH93JR/RbSjGJOlS4XjjR1ocHKLPvyRY8b2aOI/F/oOtT67mLiTAnY35QyA3rPxOPYcZjnmGONPhBA+G10E4DFt/8YDrrNgwYIFCwoxK5hDCN+DjlP+swC+D8DnA/iVR9yuR4Kj4+k2bymnuk4twuy5gHo9aq91PWp+rMnWq3G7QFznEnUnTIdizzBoHZMkLq378i06kPWvXv27fQk81UcHXpx3gS9NN8avAHxAHwa42YxeGUxl8OXFIMP0xdFJv/OEAm6Ob02Db4b7pLUp2j1NN3djFOPFJSAW9LNzrM/uJ5qMjuoCUupAHvmwYCZSeCk02RDoeU+wRZ6nrloj5n28H0in9fymNY3B0YaW9qyXJ0vqqJFQO1JXVLTIXj1Lz8NA2i9Td0uz1vfMtIW+f7ks15ULtuGZQtenZHq7pg5J3kjMZ1WrjkIUDmyzHr7do81FEomqPV14dtWQxszBWGx/ZoqnFCVeGZ8SY/zYEMK/jTH+7RDCdwL4yQOucW0gHKeOjrIEc/LSofilrEPt3il7RGRU82P6DepIyTyawrsRMfUw5vYJzbHHaLO9PXJw4qLH/tP9/T1xfIo7Z930Tkclyke92QCrTTUmVTq5Re4rJJhPbo8fyPGtqWCWZ9W2I80jPCDQCeXTvq7T+8BmjVv1ad/csWHaDU4G09jmXxvD8kiQ+x8STiHlDnMURclUlqf1tdpm+XVzRJx2PeTjkr7epJ4YlleHRmCOge6N71k8D+QteF4JfD6QCh++jHh5SPu9ZbImnK8oMlU1ltmVol6P/aze98oTu1yESb2WV4w1IO/onnlNzasK5pKu2n8ZOAshfGhf/0cccI0FCxYsWHAASjTmnwkhPAXg7wJ4G7pB73sfaaseEbyFOT2NhJEYe/R8tyENr+m1wu1u9LYAeoNfPy1nN4LtdtQQh/ly/5un+NsdaZI0xd9vgZWMxTJucwqY4S6RLH+6pqneyQlwt7doNzSub44QTvqoxGaPo8QdgwyWTFOcnIza853bZJB5ItWYWVsBSE3dUyDO+agxn54OeT9w/Fx3vvhi4x7athmq4URBcold479fz6d54pXRdNoQkDfwaB/eOaPPCjatYXVXr6/qhEq8j93i+T7ld9uO92J5HhRONMa6MX0GvE/qZG8UNgauMD6zprH9yuXdyCzuSHs3JWVrpgl/FkvQOUW0cVjmmTuMUbI72DlmHiqVEWOU1Up+IoTwM+gMgB95wDWuDcRzwEtoHdTHmgQHEF/ZeQj0j/mid1SXcn0+VsBfQtOkQmYgb8/HwJOLy5HzBbqp1TD1Wg0cK85uAadCF5wCTwpdIcHRLf2WT2QL8/MXPk6EKRPJm6NxANAEs9zb8XEqmO88QcmWbhPFcQKshfu+BSSrA4La1mCYpO3PR3+9e8+NXi4klAEAbYtbvauTjqgT6mFd+0EcfRXDX89Dg+kLzSNanOKcUGb6QhDgf5RzQSVS9hQQ2W+dw2Wue1WR0EbKhcu9eV4kAusZ7DFaRvZ93czNyzV5cOSkQ5d9wNNukICXpOzsRpedmvpYQy0RCi8Jv9wP12T3ST0YWAE3W2qzzsoofUMG9RIcNBjGGC9jjM8B+LFDzrsqQgifGUL4v0MIbw8hfE2/7SNCCG8JIfxWCOFHQgibfvs3hBBe9TjatWDBggWPEiVUhoW5QfKBEUJYAXgdgD8P4GkAvxpC+GkA3wDg78cY39R7jLwawHeX1DkQ+VWaD8HCoDFTohkZnePFFkG012NK75mEfB6lXhlb0qzPzjoNEOi1Z9Kkt5fjkFxV6Sg/OFmT9wdf/5a8Thm72Stkj3SyxfWuKJDl1jgN3ByN1AZbkepx1ZTOj1vRF6Ix334CQF8efKyBziC5Rtr95EVcYtAtVs8DT1JSqGTJ7zg+m2Y/zGCOt/cHxYmNYnOBFwzLq6NpO81HFDRO6sTac2K4gq8xX4W+0Ks4z2nPFko0bh2wIss5sYs9m5zlPuUeWLNmaJ9mntqvkBrMpO5LKNaP0qayL/vurMFaZqEcFQKkM1DBdtdTZTJzvRxmh5d0HSunyUCVIdWY5etijZn7xsP2yrBwldzPh+KTALw9xvjbABBCeBOAVwD4NABf2B/zBnSC+rvRpQ8+n1Yzwuqwe4wfop62bLfj+7246DJQAR0lcueYOoFyVAfQCSy9fRDM98cp+imVz+5jv20TubzebMcGsbsdW5MFL2yB2zt6DBukUYAiVs4pjE6oEnH5oxR8m3XKwyUDQ3+8eF6YwvguRte9JzC6wR8jXVAKSJm4ftDDyXj+0Rp4oRzaTrnoi+649cU5NhddXesNsOq/kJyw2qv3zvSF9I0dOiHBAtgS0jrqC/TbIm4E0jzNybqelIpu0Nw5C/Ghbqp4VVGuimrqemcGorT2UKr9giJSmkbK1tp+Uk9D7WwwCrEd7AT2FxedTBW5enoKvGDTKzuJAL4cXeeS73Hf9ZvT+0MF+7PuQpcX9mBw2XvmDIMDxj7Awvg+RtqCBfMhVEYuV8Y/gb8c1wcccI2r4sUAfp9+Pw3gZQCepcVgn+6PQ4zxOx5DmxYsWLDgkSOnMecE3eMQgtZsaGVsK9be2bfTS+2oF18cMltdjFrE0QZY191YeYRnUxVLtOLNujsh8TYYl3zC6fN9+Qz7026sPb/oFFluz7C0VHMfK/aLToIy6Bp3L9OgDms1WDZebi+78wbT/QrDGiTJ+k0hzXQ3hJ73y1wF0ZJvo9OOgU5ztqiMNTo6w9KYLzH6W4uRsMdRf/93L5OAk24G0pc3G9R19zxzy0wBtoFnsuSTGH/7f9I9thg1oB1tZ2OP1n7FyKWhXIWn+w3N+CrIGfzY35u15IkhMWNAlefBp/A9T4J61DHcA6SOBsC238HUxdnZdFUcuZ8n8ew4E7ygXBnaEEheP/H0/jhxPUs1c6ZPdtuxT/Aq8peqzLQG95lS5HJl/IsD6nkUeBrAh9HvlwD4PQBPhRDqXmt+CYB3eBWEEN4I4C/Kb50jSMre9n2bfgjWR34XWxy1nUdAx1v13bOuu44gfC27y5FQOT9tJ52A28G5YO803YFH7bspIIMi5S7Ou3pPOBetij4EkCTtb/rzWehLB67WKcWR0Be98FzfQieImb4QAXwbqWAWWkM4ZkswX2AMitGfa0/RnNzvvFEkXwm3rV4PgtniWy3Kgi3/nMSoadIptdAZQDrdZos8u9FJcAU3QTsw8p3r8qOA9r6QMlMgOa+OoZ6mbGqu783KoaEDVNiV7hK0BuXFSF3oldI1a3jnolN8jk7up7kyhsZ0VNiup73OzoDnhV08TSNeef2/7TZ1hbN45QuMK82z8Ja/IQTh6gDgJ2OMXwSFq3LMjwO/CuClIYSPAPDvAbwSHbf88ejCwt8E4EsA/JRXQX/DXwQAIYTIGbP29PFZgrnpfWCjsY/RtsDJtqvs5OI5rESL44VYAXCi/Es18rNg5vYw18erXNxtt7jVPjM2VPv9WusO6nAwvgGNgb9ep77HUr7zBBBEyN5BKpjvYBTMT2DUfo+RcswrpBMgEWecmZdF4W0M3POqN34OYfFrP9MM3aYWxjujD2h3KI7uY8OeNuqw4Yo1ZoZumXx8XK+GNtJdBYMARtq3WEtmn2Z9LTGaN824vWmAqu92qxaTeFMLLIy14s2aNb/1GqPAW7XAhr4bDpfmtm+33ewTAG4dt9j0J1mLscqnc3Y22pDOz2BqzyKYpW1sc9C+y8w9s68zAMQY5QNx8QATo0eLXiN+LYCfA/DvAPxojPE3AHw1gK8KIbwdHdf9/e+9Vi5YsGDBw8d11pgRY3wzgDerbb+NzmPjYFgrEevUjgzNNzKS6B+y4G42XWV13STec/t2dITQKSt5RN6RJwCvrpK4CjXAk/3N3Nk+gyB88dk5cIem+JwTWvPFrGqwds9eGcmSEUdjvtugtWLmkm/TPvKqwC2Mk9INkKxTASpztog1Rpe/FW1Xvm8ZVVJTF1bwCb9Dzx1qZ/xmjlm0Ih1gwuBWBoxaskE2DceUgu+zxF2Oy5q6qCqgoUfLz6ym/jhMrLbdP++ycp81xmcSkEYHVkgDWXhmwrlD1hI3Uk/vk5fUFKaNeWjreF61nGexFpVxcdHx3Uxb8Xv3KI5H4i4XQvjfAXxBjPHZ/vcLALwpxvgZB1znWkA+Pv2BegthVpVNZbRtGiYqgr0TzF1ZsnyxkdFywWHjgizfZFEZWjDv6Zp3tj2fdnYfODWSBQFpFjzhv4GeEqjTDFxMf3DmvCDc7zFS4asFsFAWIoCBNNKvZOIrYOezjHGARkoOyd4p4Sv7LqnMSwldbkdj05xgHsKG4Rv/5G6BVDBFjILpQX1P99qNrR2vydN36ZNA5/omsDhmbz1C6Xd8jAhJoTaCGpH4Pi0fb90TWozPeQWb1qjP0nP092W5NGs2z4sqnAhjcpcrobB0VKgI7EPsByVfxweKUAaAGON7AHzwAddYsGDBggUHoITKaEMIfzTG+HsAEEL4Y3g8ASYPHTxdtYw9OsXgSk3jvCCEQfPi6Z16sjqFJI/OMp3a0zECXlGF288asxg6bp+0OD5+HsfHnQYdjjepq9CQr5aMgmzUA1KNWX4PDWGNWcpH6jdH9GkDH0/yOaSghe3v4IV07EY3PyB1/2t2yfvg5aeS4ICtU27S6anWhC1NiI09mspgzUffsadBlc4lpA+ugAnVJpqyt5gryCuDjX91k/Zvpta4D3N/lImaG9QAACAASURBVAjLSrTYi5QW5OfBMwZBo7bxc9MGV+lZ5y0QRH1Gqv3yzJWjP/Wz0Avv8qxJ+s19Mv7tm9Q1ktumZ0qtsf2hrmAC4OsA/GIIQdzn/gyALz/gGtcG8oB1BjJvVWSmEriz7lXZ6hCaz9L0A3Nb1rRRwFSG1K15bbmv87N0zb3NZotNHzk45E0GUmF8sgMaSqxf15isQQj0vZsn5SK8JXES8b+JaOGuK4hIGVT2ED7H6Gx0jjGK8T6GiezlRerHLMmf0IXLayu6lL2IrqSM1B/1UpUti7wW2FowW1TGo4LmmLXLYC6rnvzViwBYayDqPlhfpBnuhlfYpgKLfW8038w9ojXOYVojAFg1GLkNpN9hkmaZaBoWznoJML43/j6lLBGeOrLTKlvv+RBttiS73P8WQvh4AJ+M7nl8ZYzxDw+4xoIFCxYsOAC5kOyPjDH+Zi+UgTGQ44/21MbbHn3zHi6GKQkZMbRGoDVmnrpZ0IZAXq2ZvTLaNqUs2KCgfarZ+MdMQhL4QNoBL+vEVmhdPj5u++Pu4+hEFobto/7EyHdM6uPEx7miv1Z5eCr9373axroP72f64hxjrox7/T+gW5y9L9+710VO8uomffn+me3xIkYc1qa5rDVjYN74x1SGpTHLXddUtqb1+il7WnVJKmHpM6y9zq1Grq/h2VVzxjJd93DOWTqVF8JL59do4U/1LePhQGz1lcczW2Ne1TadY7n0Dyl+m/Q+5fuUvNvWKi4PGzmN+avQURbfaeyL6JIJvU9BLLVeZ9OdnaOjSsCRSMzf6euwdwCXmVYRSD1rx1K+24wC/6IPG5dowaOJYO7KJyfA7T6u9E77XOc+N3DOR8D21vRCsc34b7GY2WMUbcwWVup4wLZva8HcB9LgPcDz7+k33xv/AX2oVh9JSa5OOniH+Xy2tp+3UwEMTKO2NJdsCWY9ja3gC2BG6TQ3J5DlL0fv6XUsS7LQ6e+DPYus6f7AMZOwS+iP/h00cH1sujY55aFddI4sB1HxTrITrGkwGuos/JYnXlt0CT14WL3bU1sOCRrJhWQLj/zyGOMF7wshGMuaLliwYMGCh4ES498voQuDntt27SHeC0l6T9gTb6ALAa37nTovwKB5bMdtmgphMP3AVAZPCfU0kq+zJ+syn8N+mud9PIgcd0T72CioF1a9U59SIAmtZrK9xJCPuWmAtWXu2qHzV+YnxMHGDD2Rt4iBM3S0Bfq/vZZ8/gzwTF9+5hng2WdHjfnePTzfF8+UxqwDBS7Ffqj8lTl/ruVtoTVmTV9wmcFGrhx9IdDpM3WSrTlIqLW89/Vm1B45Q62e4uupvZXkizXJS6YLqqnPP3stSeba49afWcj9CoJT1sc3dIy8z7pNZ56HzHqHuqX91GZt9GN/owqpxz6XWzqmFDmO+UPQpdS8FUL4OIzP5y7GBAjvU7gkrkivZKzLFdT0skESDSWCzUuSw/wykLICPA3UfJblsif75Fo7+vAsCzQLZpG3HA3Fnh+rCljXDY42PX1wcqKWupLES+fA2mJiN0hM4xMnMd7Oy0e1SNlc8b44R5daGwCeA+73LvTPPNP9AzqhLP8A7O5dDLkNnj8d8xycGxyz5RbnCeMLpNxzTmizwBFUSF3k9H7Lum/9Ho53duik+UxfcB9Y1eNK8ZrW0PBcQ6Xf1zWwVXSBlYfkmL2OztLlpLxhm6HDkrzjJ2vFex/2gdjR6RIUZFEu2vvGEsxWakwPOY35MwC8Cl0Gt+/EKJjvAfibB1xjwYIFCxYcgBzH/AYAbwghfF6M8SceY5seGUSv09NFreUIdIpGXqRz8KQgbWXfANHRoNvWD6nWGrOlVQGph0ay9A9pLqt6nLpe1qnxj6/DdW42wNGxrORwa6Qyzs5S7Vk8OVa3MOqScsO8nJWlI1ru+KyPksYcZQHWe4NWnGjMz7y7pzKeGw4TLfn+aerXrY1/lr9yznf5gso5jdmiKGR6a+0r0ZJb5OkxxkBd1GmKE89Lx0uZOVzb0ZgHo7WR2W1Pxw19/Xgsr7fAkdiSMU0BqmkgoNMGK6OcGP4MlFBFGtpXSNrm1aXpC3Hc2lKZ54+HJCYqOfYTQgg/r3Jl/Hcxxv/+gOtcC/CEu2RJeUYAUImluU65OC3o5K+XZlLTF5xmki2/E86tr2u1VVMomsau61SAs4eGXJ+DEOQjPur55zsn5yoBvazGTTk4nthgKnw58RDUPmCaXWKL8Y2cYxCH95+nlbGfH3nkZ57pBLKUn3124JVPT9M0jblcuhYZo9M3eh4anMKRk9MwxwxMn0yJV4ZgzjvDsj/ofE4rErrsjbPejFTGeoO0EyisSBqzO6hQcLpNu2YUujsaANZU5j7YYBRew/1IfWqbfItrpLQGC+oSaIGrYQ0M1mDKaxta9MVGnSePzEtWZaFEML88xjhQFzHG94QQPgvA+7Rg9qD5Iysz1q4Z+7I2jgj0toZc4TTfLCsiiOjyNGbplDvdzv6E0HbRUENnqcaPYrexP2oR3vLx3j65jyD5nE9uq0x1lATplnTPFp13Kie3Z7DI4+DaSyASf80L0srqLvfudcIZAJ59T8Ixnz7b4LleMD9/OubStYSxlHm5eW4NB3vPGf+sRTa17zJ/rNrIVSKcPe7ZWziY3Trr3qhnucsdbYD1MR1YG7m6gVSjqGuEvrOuq1R0ySHrTSeUJYyZNfa10t7Xoj036fPfIM0uJwhIhbEuV85xSTv7vyuU8dp6v5srm67JGrPXHw4RzCUDziqEIa0YQgi3MH6FCxYsWLDgIaNEY/5hAD8fQvhBdDOAv4Zuder3ORizsAQyPdLbgFTz0VqycGuVMfxyykR2umdaQS9LxB4jXKWVz1Vr+Mko3gLbXsU7alKNeUWK02YD3KLowSePe1WUc2pwbmd2Obkt8VyiD/AEr4VJEuy2nSserxo+RPGdKY25bwt5YZw/u8Xz9zBQGfdPOz5ZTmeNOYlUQ6rllmrJXC4JKrEi/a4K4ZkPASfu0QEmyeo0Q1kRd20AWnG1SHW3uu3uum07F06g04p3ta2lc5ntH1pjTi5PZU0XaL5Zvldt2bDQYn7FbqsNDP6+OEWXpjIE4pkjbSxFSa6MvxNC+HUAn96365tijD93wDWuDSzBzIJYGxQqpPxgEkJKgrlq0218DId58jkDrYGxg2gqQ0+J5r5P6cTcWYeO2HRhq0BqMBSDkFAZR8fAyUn3pNYnz4+CWWcmHziafXcyW5I47HGQjDuMC8BuO/piEMxnxGuf24L53j2cP7sdNvO6bNbyP7os/qie77FVtgTznO8yf4iWoLiCS+0AL5FWUr9xAbZBpMR0mB40Xm1s8bBEWmOuEyhla5/O7jYI7KqLEwCm9hzdEuZ0Nd+cZPimE7kLDmpCO+WRLcWLkRPg2vjHroB8zCMRzAAQY/xZAD97QL0LFixYsOCKmBXMIYRPBvAPAfzH6LT0FYD7Mca7j7htDx2W1RUYR01e1EhGOhmtObGll+RFBwDo9Imc43Y4h9oVkRqlmAgo0ZzZlUjqTrQRsedcTDVm0T5PaCHKFxzfU/5VNOZLROB21xkFmeawkiaUasyn98cGnI5eGef3mnHzaX4lY/Ec0MtEscMeP3dPS7bojjkXOcD2LtDQhixPky6hMlqnn/F+oE93klio49gAbbnm1K/tnA9Th5WjTTNtlqTgFOMhur5vuasBU4MfY/BOIo09Se2pjmW6MacN5+6YPaKsQBK+F+5bD1tj/i50K1T/GID/FMAXA/gTB1zj2sDzk+S113gausf4gNifsUXq0zxE+7XKJccR4LFNhTE7kYlwln0cDupxzwJNZYi1W2O97YQzMAo1EXKnp+xe1eLOpqcS6vU4pU0yMu2m/DOjSDDfpwakgvn8tB3axTmL7iu3uEuiLzhbn3i86FWuNWXkCWm9nQdKLmsXOQssgD1hnFvnL1ECaNDRvsaJn73yAFpvidCTgZb9J3WF7T7pxLkBoAQVCU/tfeENTh4lpKlHK8ScscdUMM/BEtza6waYCl2mXyzueQ6lVMbbQwirGOMewA+GEH7pgGssWLBgwYIDUCKYz0IIGwC/FkL4OwDeiXElzvcpeEEAUt5hHPnEV5i1qpr21Ybm0LZKc3I05rZN6QrO78rGQK3JceCJpzGz1ViP9nJvW6TpL1ljPjoGjnuFdbMBjjadar2u740VtfuUyjg+QrLQK4M15mEZ48vOK0M8MRw/5tPTzuNCymzs01F9nJyI06jyLIdnGhw8wrQGa8L6+fP5nmFW9ytPS9ZlbYRmWLMzLuv8w5ttSiUw5JxNu0Wom/Eg70C6ENtymY6zcrxY8FJwlhpEc7OJ4NVNDkTJvjb1ay6Z8QDdOTyr9vJfMC3Tqm0lKBHMX9TX+VoAXwngwwB83gHXuDbYG2WdIapV5ZaOT4JNiC+Wl82JhqzIv+EDo2vyNYTWyAkGOcebhq1gZ+mqkK7oyyHhnDT+koS0JN4HgKfq59MPQ1MZnmDmr3d7OZ7DVAYJ5nh6fwgWua+EMUf3aV55cEtsVJnumZ8nr1LtecK0Rrlk9q4FzZww1i5YGlZObw57rus0N/KFiqSyssN1wrzbUVXtrCCXshW9Kh4i1oryh2COm8+5t0XNJ3G9FufcjnRlDrLsFbdPengDnzfmYwQP213ud/viOYC/fUDdCxYsWLDgCsil/fx1ZAapGOPHPpIWPUJYmqS22uoQbP7NxrdhOpMxgmiNmUNqLW1NaA2mObjMRkGPythjHNnZeMl1JZp4k3oynF8AR73Gdf8sDRS40+dJDm1bTmUMGvMu1bKVxrw/6xrAuS48X2VJ5+ktrMr5GHQQCD9rK0Aktz2nCHpace3sm9OSBdp33sqxcm5oyDs6jpch4xSgvMoHvzZvdRMvoZFozx5tx+db5WHbdNMETD1wX6+QUilzq5VU/ZSYKQfvFC0jBDVsjwueedd0/sPSmD/ngHreL9BifFFML8hvi0rgj2BNT1N3TA4IKJ0SM/+sy975Fh1j1cvt53XNdtuRrz0/G++JP9w7uI8gLlTbLXCxUe5yRHqI21XTkDDfjsIZwO6iNZPb67JQHJcX9np+chkWzHKfklfXeod6cLSgvSj0R8xuU1zWLlVSZo+COW8Nj4qwIApAQ+9zLdGfG1swsxtbXaerSXsCjjMkyjPX1Ia0x+LIDwGfwt8n9+/EvlNIaVjuc3I60yqaYrGO26njOO8Hv/NSuIKZKIwFCxYsWPAYcaMCTA4FUwFAOmIyReCFZHvbGHrlhceFgUppp5oYr37M8SWM436FzbXMjwcqYzX1iQU6zXm4yCXitjFXs84tDXVJVAbTLxMqo7+09mRhaoP9wnM5dxms/TSYTmvlry5z0BIbkDyNmcFts9LLAt07PCJa4XI7pvf0lhrjDHA6N7P+bRkGtVfIxOhKPtZWVjymHaxnb2VVBFI6kmc9FVVUQmmsKvUdt1PNWLdDvyN9fe5bfAzTGqV4bAEmIYQA4B8A+Cx0C7u9Ksb4tn7fV6Lz/vj2GOOP9Iu9/kt02XFqAD8eY/z6/tiPAPAmAC8E8DYAXxRj3IYQvgHA78QYX++2gcraG4P/WtB8c8I9Gh1CkKydJtvUMV7Ek25Pzl3IQkm+2qgE865JhZwITC9u5LhpsdlcYFXTgdbXQHl9RZBuidc+JwFsRfRdXIxc6mVf3hmC+RIpl85Cmj+SklWptSBtYX8w+jgOKNggFdTyrgNSr4xsYIncAxPgGN9BXY+0wuW2oyguSQCLYF47wpiFtCw/xft4pfakXYpvZl6b6QsrEIbvqxSWwLSEqQ70yq6U3f/etzCTkFnt5e/YvL5zzMN2l3tYASYvB/DS/t/LAHw3gJeFEO4A+EQAnwTgJwD8CLrv69NijKchhDWAXwwh/GyM8ZcBfDuAvx9jfFMI4XsAvLqvaxaWMOZyLqerhvXgdYdIji/shdoPmd3dQNu9DqFXebCSelscphVRpgUz+4ryahWsla2qFpXRwy1XL9bMLe1Zl1kQ77apu5wkvefcypaQnjPmWcY7OZ4HZN1/WDBvqKw1aDEArVWZBbZA+hb3KYlkbMlXeddMw6H1qjZAukir1pg5Z/KlWvVEztmTtGD/YBHMzDlz2bKzeLmlBZ6yItBrdiZCsk2jALnN3m92gWXwtr1qe9OmodcW33xVlJyfBJj02u1VAkxeAeCHYodfBvBUCOFFGPvioMT0x8iKnNKHY691fxqAH+/3vQHA5/blU4zrEy1YsGDB+yxKA0wqPHiAyYsB/D79fhrAi2OMb+1d894K4O/KzhDCCsC/RkebvC7G+JYQwgcCeDbG2HAdABBj/I65Blg3m0v76U09eGBNcmi0+ZFuTmnWGm2A7WLDLji6jXr5HUsTm7SrtaeeHFHG2kVi9d90K1hY3huak2S6ZLe1KRPP28LabtEXnN9C8lsDY0SltzqINYOq4U9R9XZLY97QP2D6bnKeGALd14K1s53OBoe6q5SKsDRmXgtQgookX0rTpAmvPDc6DvLxPDTmtGSuWg4tdSXU3khJvWo2obfnMMntRPcTaGYApK5z0u+8WdYcDgkwucCDBZhYsiH21/hWAN+qrrsH8KdCCE8B+F9DCB8D4A+8OswLhnDfakAu/ysLRaEG9DEe2uE/Z18PFuaaB+bIPV222mwJCOvjZ1rD4jR52SsroiwXsds0XaJ0OS7J/0v1t+oanmC2tlt5loW+2GH8KFhI6+WjWuRz6wJTYcwflRbMc8a/jfrN9AVvzykH3AZPrlnT/aG/EHe6auj6VSqYWRAfH4/vateMhsS2TXMr89Se/Zr5Xes85B6lx6HOgN2/rwJtvJRBhl0CZZ8FTWXw/VQV0Mh55L7YIKU4BMP3R3IpxmiyD668CSG8IoTwGvr9lhDCb/f/Pt87T9XxmhDCr4UQfg3AO9Bp24KX9Nuy6BeB/ecAPhPAH6KjQKRPZ+uIMd6WfyXtXbBgwYJHjRK5lNOY/wY6bwzBEToj3W0AP4iR58014HUAXgcAIYTPBvDaEMKb0Bn/nosxvtM6L4TwQQB2McZn+zUG/xw6j40YQvg/AHw+Os+MLwHwU3PtEFi0AI9M2vjH2qd3Dntr6GO05lUSl6+1dO8F6fh9OV8vUMnTZete9PRMUxmW5xuQGvK05X/OvUqCEwYqYuuXLTe+7RbYtqlhj41/eqkuKXvaFz8nrTEP7cf0vQtyxj+mMlbquMrYrqH7jBdktFf7GdaKNusW2NCyY5p6EG34lnoYkrxrb9BU7JVheWJ42rJ8K0H9lnIJ+PvS3wwHz2i3wBJqI0nWVI/3ydi3QL/qVtKfDvWkEuQE8ybGyJzwL8YY3w3g3SGEq2igb0bnKvd2dO5yX5o59kUA3tDzzBWAH40x/ky/76sBvCmE8M0A/g2A7y9twFwmKCAVZCWCTUP3PbYaa0syMP2Q+Do6hJz9YT0h49EXGok/qBKanJTJizCTafC+TbnLPXX2oDq6XiWco9MsH1jmoZO8wm2atN6jL/T6iZ47ohcqywNjA19oa8GsvTCYymD6Qvc1gSeUWBhzf9L+2h5/foHxXnlV6n0z5X+TvMnG4Gy5TyZeGSSMvaxzlgIh92bJyBytob9PK2m/tRZhcO4zua7qt4436KziVeKiKcgJ5hcklcb4Wvr5QQdcQ86PAF4ze2B37L8F8HHOvt9G51q3YMGCBe+XyAnmt4QQvizG+L28MYTw1wH8yqNt1qPBHC1RSgtocEIijuXfq2NaZ58H7ZNs5Y7VuIrBklfwZoOGpy2zcaeugbhJtQrLkwOY0iVmfgtnSqyXiWKDn0df5CL6tDask9DIXzlftOUSjVnKYvzbqN/or2e9J8/wJ+WGyhZNI5q0plqAdLWeSR8UT4PeE2drGHMZ2k9YJzE6NAWovoT+dqxj3LqUlmytGC6eRFaEI4Nnk9qXv23HWQKfv2qRpZZKkBPMXwngH4cQvhBdhB0AfAI6rvlz3bOuMawXq4WxQISy5cmR442sDgX401DrY6/oOKuduUxYzJN70INEVIJ5Nsn6Jt3WtukUkTuyGZ3dTj9kKzhhkk1N+E2ky26xt4Uu62drCUN2P9TCeG5xAgELPxHyK6S8LtMcFdLBwOT/+7/eUmMWfSMC18qkqHN187WHBSKa1E2yJk6V32ckISUDted9YbnJhYq8RZDv07nc45biFFQf5jUH133flehGvk8BB8/wup0V3S+Q0ho6QZIeUIF8/9FwBXOM8V0APiWE8GkAPrrf/E9jjL9wQP0LFixYsOBAlPgx/wKA9zthnNOUZb+eYs5NpbQ2krOcW7RGNM6xwLQGX3fOAmxphQ16zcXQinJg7VmnebSmftrvVQecWOcnq4m3o+YozeUZCJc9TJZsou38IXiU0xw1AqQzFiskm8tzlnu5lmXAnKNvLCpDa9GynXOKrJAGiDQNICtQtXWan4MvkoT1q/7gYdA2aSYk7ZW2zX0LFg3EhryV0oo5ECoXGGV5KbVt2qeLvhPk+6SHolwZNwU6R24JLQD4nccSutZ2Pv6Q6U4J+KPmzi5T4gr9Rykco+EKBKQfmLa0106j9dSXz2cqg6kU5rt1oAJz+YcKTOu31/mtqWcpZ6hdLlfwqQwv4IlhLaQg5b2zPScMrMFZ90dvdXdvGTV9XC7Cj709hJqy+GWPV27VdnYFtCL8VtXUEwOYcszadY5tJl7koMZw//DfTSlK+fQFCxYsWPCYcKM0Zh5t9cj7oHVqaA2rxThd1NstzUyOser3DIySS8HSlti4wuWdOq5SGrOXcUvsf5dqn/bYyAWbWFpZSQivpYHwvVkGPks7LfVaEXjBHlw/X0Om5Hxt9rLx2sK0TFS/mYKytrdINU7dbus6rS7Tc9czG4+amKT0LFAPB0qAGzPTZr2NnzVTDJaWzGXRpNkrI1lpG+N2zy/fu1/WkvndPCw/5vc7eNM7ywtCtnsualYSIUuo6o8MxnY+h9Gqcu5Ya7vu3FoIW7jEmHRdw/rYNpgKZ4Z0fH1uLjrMLGM6SFkcrUcLaBvBISleBTkhbUHX7wnp0mt6A/VVXbIETA3Jc7Y4Vva4mbTzQKFckc1h2JZpn2ClyixA14YA1pzwSlEXLHTFe8MLiJGBKRr32rbT7xWYUpiluFGCmSOlLIMLC2L5LbC0M8AXlprrK+WY4fzOXceDvh952XMCWgStJ6A15oTz0J7+wQlHzVqZFSmW+8C9BE/64/YSBOWEpCUgSt+N11arbq9/CVjrkmO87HgangY/B9HwvIHSE1qAvS/7DkUQ9m0bjMYF7ZTjanKFq5wyLzqrIwLZOO1BPwtrNrFv0+/7ECHs3duCBQsWLLhGuFEaM3siCDjH7Zr26YQoczQIYFMPHn1h8aR6lM0lsPHO8dqnNecHpTW85C+NGuot1zntLscBJnqqmKzf1v8N6q/nXqW1U8tjAphqtDltxXv2Vv/IJRTSx1ntLwXXOadplSTV2cOfwci72cP2YrgqQqXaPlNfhSlNwSk9LU+KnEcF9zOdonanKDfLLTC2qStj4W24uJGCWbvBsSGMXZtykXccem1BG248WsOqo+QjL51Ky3G6T9a0XxtDuTzIZSWgPcHsLd9jRRHO+THnPnav42qO2TPE8SDs8c2HCGiLarKab23TxjqBJUStBFe5qDmofR6twQvTissckArjHN+8V8d5747prKFNPd/MA3fSp4x6hBO2EhSxG1zI9E25loZWGrzkW9qt0HqH/D0dgoXKWLBgwYJrhhulMfPIxdFyJVoy1PaSFH+lI6WexnoLUXraFl/T2+dds9C+lxy8qkatAfDdiCQl6HBN0pZ4isxpP7X2PGhhVL9EPnKwBkN7X8g5TFNw5J1Mi/X95AxCcn+CRFtE2uYSg7Ll/aM9f7yVb9bqWrn0tnOamFBu1gyGgy0sQ+AhdIamQfhZu7kxlPbrpfTUASZzyYly+3aG99CcO6dGpf6W4EYJZg47FXieGBHpw2EuWlv3D4XOVKfpBhbQVtt0Xyh1neN7k2exVvtKyk2G1pi0oW+stcwUf/CyhNXe+BAEpVnzrHejswWyMGZXqUPWhbNcxFhgyX1a79qjH7RrWElObq4rx1VbA5Zul5TlFe8oJHtVp/fJYfxJXQfO3TkEGkgHc30clxMB7HhllL4//s1URpIfvCkTzPzdXjVp/kJlLFiwYME1w43SmK1cC7npHfs7Wx4bwGg8m4N3jJ56ehZ6z0CkoQdwb+rMxsuc8c8qs1ZYtVMNGnRcSeSfLDUFTHMwW2BKQiNnyPP8Vj3L/SEJnSztWfyzqzbdLu2Rnzk6qab9mrJgLVnTJlftk2LIHAzVdG/7ZvS68egBXc75PSdtcWYpuTJ7YuiyptcE3G+9Nmvjn4589O7J8hLid7ZQGTNgWsJzSdP8JXtsVEg5aqhjdfmQbfJhANOp9845pyTBjt7n8dFVYbkyOvFQN9EXTBFoK76Xd5nds6xvIBdePccXz330OcHsLZWls5EJ6rp7TiykLIGgbRsa7LLpCeOKjtW0BOi4uaRce3TCZPAsaVNPDHk3vCAC6Fj+q7dreIOeRyfpd6G9MqwczFohGL7nZmpD0Kt+A9O84V6EY1V1CfKBlAbld1OSEG2o74BjFyxYsGDBY8CN0ph5SmFpvDywC42xpnMsrSY3spXkcGCIRsM+xjzdHVaZoOta1EVOy7TQYkqnyPbKKbdKc7Is91VlPzPLui1O/IlvqLoX/QwfpiGv9DiGlfRmDgnV4GmSzm+tJbM30UAzGdfxptIWzSUGR+53K+PdsmeOPOdDvBX4eG+2In/5Oet3awWS6OWfLN/pCl3aUfad5nMsY65s55kSX1NmQ56RdqEyDoSVkEjA/UserLaoWwm+LS8A0D4Lel05fR3mGi1aQ3joB5kGeZy556EB6bgHLE2l3eViO+bm1VNvhg4W4Y9ibbhKlQpYDRYYiYg1GAAAHnNJREFUXI6tz10KSq+pn6HeB4wufXPCGOoYbacoSX7F+3UftAJHclRGUl9GQM/xx4A/0Fp/vWdvtWsQ2A7HbLlsCqyglVVFOcdb+xv2VmO3cKME8yHfqQgoS+NbwR8Fmf+s1W9LSOtr8ofEhh9GA18DL8l6xu2cA2tS1vYKaTgqYAs2YOpGxlreHE+u28ualLeoZqmByttvQml51vmlnOqgEJC2xTyoPJcSYayPsQRDbtDT2nPSB0UT1LOhzLP0DGSm9moIWbYNWMLX+5u7PmvP1jmWYM6txqLbI22u6diraswPolwtWLBgwYJHgBulMedyD2gcMmJZiXFEW+YVk+fSL4qmIi/loKi8Angae+WUPUzSjpKWx5rYJH8CneMlePLapT0wmL7QaR4tTnJPGoy0h7U3z21quKbSzHKadVE0mKM9Mx3mBafI8TlNmhN2yXEN7OdswXo3OrnPcGzle6x40PdfSl+U2AKYctoTZbZS1FQu8o9ndxrStqjaJv2RPVkqjNc5xCvjRglm68HkhJT2ibUEmI7s4TIvXc9+0J6fpZ7i6wjBEijabLJP6i3N0Xvo9a0k4ujb5C0Yqn202ZVxTecw9BSXk9Z4rnFtmwpAL1G/1X5dn/XbQmkUHOcmHupH+mHnBlBNXcj+nTrWMgrmwLQGG2m1YHOVjdKR4AExuY4xqDa0bc5W4CVk4vvmLHZ1nQ5W0h/37fhu567JWKiMBQsWLLhmuFEa89zN5lZ9qGg/a7/eOm6iLR/J+ZVtoGLI6DyMvDRSb6aHA0ipD/ntjbZeYIZOgzpH+cwZDi13N9aKG0ypDAZr9pbh5BBjkW7TbJtnprGMOQ0olwSJ+wC3zUrmw9uioT1zHUJrsNcOp7sdaBKq35qt8Htjak2MlJZ26KHEda4UJVRE6fWv6rUj5+1njImyTbsFluCxCeYQwkcC+EEAHw/g62KM30H7XgngbwD4oRjj/9Bv+wQArwdwC8CbAfw3McYYQnghgB8B8OEAfgfAfxFjfE8I4VUAPjzG+A1eG+Zu1qI1PK8KqatWZeaU1xjXGFsrHnS4RpUKhaQTESmokyrJR8Wuc5Z7FcOjL/S9JQtc0nHWVNoTEBoRo4CQkN+57ypHM+kPIudSxfCs9WyFHyIS+Rirgc5UXso6mTzTLCVCxRTszsedLCjQkgCdP3XWk4f70+ChwddofGVjUpchpLS3RO7ch8Xx6zZYv+ewqpRtYOb8cEAWo8dJZTwD4CsAfIex75UAPhHAJ4cQ7vTbvhvAlwN4af/vM/vtXwPg52OMLwXw8/3vBQsWLHi/wWPTmGOM7wLwrhDCZxu7ZSyJAEII4UUA7sYY/xW6DT8E4HMB/CyAVwD41P74NwD45wC+GsA5gNNcGx4WlcHReVpjluPX6BaK9DwHhmtWSHxYJ6M9GagEbBRjaEOPN+rydpneMn1gHcflSQReZni3Fg/Vi4x6cCmZjEZcqrkBNn2xb8a27mB4oMxgMstq0z40aJxVnupw63fO0RrooEE3qSHQg5f7u3LKe1LFvTSdOeS0ZO0XX1pf7rcFyzistx8Cj54cZnAHaMzXhWP+SQBvBfDDMcbnQwh/EsDTtP9pAC/uy38kxvhOAIgxvjOE8MF9+UfmLuLdrCeAPK8MnVzfFNLVNLGK5TnAF5UIOjmuJbefuk3pC14hOtJ2jztkWPdlUTb6OCt7ljUd9D4KFtKWoJZrlkBTBId+TK6rFFLvEaZfPCHt9Rnh+4fwXCWkk9zLBe0vnXpXNNjLdSeNLgT3pwlIOHvcvtW2oW6HL5YoSxjXLhG4mq6ybAWhSj0mNE2S445Nigwwn+8hXDjjWgjmGOMb0Gm/Am/NygULFix4v8cjFcwhhNcA+LL+52fFGN9ReOrTAF5Cv18CQM79gxDCi3pt+UUA3pW5/n3+bd2sl2gop0nydvZVZopDtDjPc2C4jqEtDO2l5ECxAtaOIVBgBaToFTH4/vhevNWjh/tRx2iNVTDnyfA4oaenuYACXsKKvUfYQUZrj1YiLF69grVlOW6goEh7Fu1NQ2vRybTY0ejkPk2PFBVhIkU9E+DgFc9jJvH+aZCVJId4x0zORUr/eND9rjS0XvJb8DvIfZ9e3Yf0e5ZLMcbb1jGPVDDHGF8H4HVXOO+dIYTnQwifDOAtAL4YwD/sd/80gC8B8G3935/K1DPcdAgh8s1a9IX+wDz6grczr1yT8JX8Dd70zvuw9JSIo6E4exW3i9vS0O9W3edQJ5XlvjwBbPHqXoCMtPMqU+Y5WFOoOfriqsEN7IWg+eY5VoCFmubia6RLeg15N+gYfq6WoGBYfUiUgcF9Tq1m3W7HtrFA1h6aXiZFi6de9xWwkLsKfz4Lg5cuyXdiHTPpO1R3Kb+dcPmZYzQ8Ycx4nO5yH4KOR74LoA0h/LcAPirGeM855b/C6C73s/0/oBPIPxpCeDWA3wPwBY+y3QsWLFjwuPE4vTL+P6T0xNzxbwXwMcb2dwP49Ku0QU8rBTkDnzuVl7KiKLTm4ml1pVbjiuoWjWTVpsY/Bq+GoS+du2eLstH5PVh7H+pRF2Gn+0cdjntIDoLSABMGa8hiCCwxrFoUgYBzbY8XH5+ZlWLSyvUATPua3i5lzpssfvVo7H5SYaox8m85XWvOK3QGaqmD6Yc5zV8jZ4gr6VMlGrPUa3liWM/AwiQXDJWt67cHWMmuhfHvccESPqCyXtLeW+I+l4+hxDsgN/0p4apYkOqVvEug7117mbi8ekYgPyhK3PS8S3Jynau6Ogk8oSTc6yGJl+QDt3hZHaHptZ/7wx7pczeFitEvrci82CJZCklTRVZf0gLac79j7xMedDgx/SE49N1exXVOQy6VazP3O68NvD8ugtkG+/6yZsyrlPDKA8wZe3lh5beFSXRTQRu9kZc186qdumQBNsdstgtT4atnClDb+V5W1fSeH0Q79gywhyz3brXjKq5KbLwz66e/3nH8/Jly93JaJ+c4huBDoAUJv0NWLkTDZV489mXLkByR157ZGMp9KEnCZAhZEX5eFrir4CqaNW9PBjjVFivRlY5B8LTnUjxkvWfBggULFjwobpTGPKcVcoDIqu5Gbmt5dPltITcKl2hDxdwzT5Fo38PkmD365kHBWmQpZld3Jq3G4xRzz3Z4Hw+RF9cJpg51xNdao46I83hYIM1DzZGSns3C0/7142AqRiC0hreKPLve8czAo6ncthj3OYdDNNa5PB6WFjynJedWecnhRgnmkoREnHSoVksWWX7IDI4kGnIrz7zEQ8A5e1kwc0Re0h59PpUT7hg+fWHxaw/DTc2rgq+vbQG563kRXA8Kpoyueu6jQsl0f1Wlg5YcW9e0nYSkuF9yn+LH2NJ2LZy9Z8RCmtfX5LoSQX0AL11iGJxzadOYi/zzvmO+jmkrOmBkXqiMBQsWLLhmuFEaM0NrCEA3MnLSobq203bKsYCiJ9pxhQfB3PQnadMBwyRPvS1DIDDVdqypY6mBLUdjtHOaQgG4+pI28Wwkl1uB4bk3eW3RniI6KCN3rgXvvpL3Ymhn1n2i8o/Jtk8ZcYHumfFsUgyAgL+8GdMSmjpjLw9v4WIrP0qJwbAkCtDTkHN9U69uUqplW3V6/WrxynBg5SNOkvbQVE98PtcOxyxgTpOn1I1BY5ROra2PS29L8sC26XamS63+q7eVUgYCiy4w8xwb19FttNrmJk5yqKRkcERmSn3AO+D6AHuK7x0fVDl37FB2DvI48yT/suqD/Az4WXFUpk4XIK5zIqQtxcV7dNpbA3ROg/F56FBvrShwQn8B902mO4b6r8Ahzx/44IdZ+w6hwxYqY8GCBQuuGW6UxszwNLFEi1DGkpJpohv10/pLFvGKvjqBiq5fQxsCZRtwtWm1B922UgOft/in1iStABPWnkvh+YtrSolnMPsDNOk5rd9KCGUZMxND7AH0g4X/v72zD7akqA7479z3wQKKAkEKBbNSfkRSCQRQVqNxC4SgqYoaJaZMGT6M+SiSIJWYkKqUyH/EIiEhJJAKKlhlqohIEmOqwGUVRSGSFVmQjwCCJJTE1QKUEMPufffkj+6+90y/7pm5bx9v7/OdX9WtO9PT093TM3Pm9Onu0yUzWW5mKKVZ8sOSTA92ZmnfBVxtXiXzhdWKod2BWKJqFhktj1uiVt5pGk4r6fjd2zw3lGDOm0GJ0nCi0pTq0gvSNTzOvvw6Kp9jTSl9WDY7rNLWzH0It1F6kVYDa6vM0y55vhuYeA37P8vrZ9pJBPnknbSto8n115a8KtmXSy9s7s/aMqD8AYLyva89L+ncUv6lZyulXbOx2qFzo1H9mWkzQVkBOu3QuzztmimjzS7dh1K5atfaV6B3vS8rfZ/clOE4jjNjbEiNOW9ujo/nHSWVY23kA8obWlnSgPKTTLMsX8CzT75jjT+mlZsMYLlmkI8trWkeScvquyJ0Q6syBcknvrR1TBbHm6+gOV86bv+HZtmulQwmKdWZ1fjTftFl6Qrys1iTzVIW3jcv2zq0K+XApD5yV7m2ZdGHPF5JQ83rx973XJPO39upJ+107Kc8LX3yKF1Dnr6bMnpgX/g2u24fakPirF3Zvvy1G9SWddtsJCsUrQP2cdgq02bSyVcOLl1s24egz8zDafx05PfE2vnTMes1zpoySi+kzXpUCU+s9gST2gedEZM3Oa5YXRpBVMO+AxrtzSXFxQrJ3N7cR5h1CamaME7k9u7S8b60lXclwrRmS7fl81EZjuM465gNpTHX3H4m2rTlPj4wRiPGqwencLtkUakZaL/yVvtto+ZHNvlCnmasbsnlo6XkXS/P1+7n2ryd9pvKNU97q6FkUip5tLPUzBf5NNlR5X7k7j1LcaBbS+5LzZSTyjmO19LSanTYxudOo7ZsfTCXXIHmNLzOmcTz67X1UelzLo7OSNtd/qwTJU24ZtZIZVvK9ruoacal5cNsnFLatQ7LlbKhBHPN6ftKBXL6HzeJh5mQzgRBqSnTOtmi50uVmBuEAffPhSmjVpYuPwKD0fLwpdhUtgKoZMPuMjH1FcYpPBe0Jd/KeXibIOkS0une2vSKppBR+fxSmRP5KIjxsz0Chs0lzvpgh4xa50uN+0b52aoJaLJwu+RgmxnEkpsvrHK1x+Rty7US80ZtFfTSNdXue9q2k2ps/tO8ixtKMJfWrOu77Hrb+OQ9RktOmspwuHwYVp+OjzZq3sSs9pyvuWe1jd6O9M1Hy858zIcS2riWUkfcSmZBrqTDb9m063Sc/sK4Fm7/l5WlVH6WC6BSh10tgVQW+3Go5Z/KOU/8ONtnoHDf0v6yMg+ara7hsLkySUkoxW9BA7uGYKkFUhLSpZZkLnz7+hFPlIbYUTne9dHI3582p0xt6ffBbcyO4zgzxobSmOfM1SYtIp/R1+ZGsKahLRktOZ9NNtVXclC35eZhDf8QRpMu7U/DgHJrwtZT8lWd55/IzRZQsJu21G0XpXNLoy1K2m+uQTNFONQ1qFrLZ0i/l8xqW7nGVtPmE3PZOQ37q1m9us9IjfRs1dyD5qaNRFo9x15LwpY/xamNhEmkVVFyU0G6hNT6sBq01V5LPjnatNA8/5IJq9YqgnafIG3n19hYgrnSFO/bLE/kzWVrvkhmjaXh8sU7u2ySOX06a/IyJjtzX9JLPH54B82PVnqZc9/UpfrLaZt1B/Up6rUxyfbD0zpEzvzn9T+t+SIXhH2apiUn8Wm/Fr80a66tzDb/vAM3vx3JFDEcUvQp3mfqf8o/HZo3F2+vt4QtvxXMQxNe6zDMzRVDE24FcF4HtTL1Ec65manNjFUyXwxYbtoppduGmzIcx3FmjA2lMZc0vrSf/msz3GrDrobDica8Z2iGy1HuPNobbNmslmnL3zBxZPm2+QMeayWZ+aKkMVuzRlud2XKWtOdax6DVpEtmmZK2XPMrkWu/05ovck25q4MIli+xlI4vUKY2vKr0DNW0tnR+TXOEoOXWOh/b3I7a7dQiHDDpYJwHdttzzLY1U1jN2GrMabv2jliNOV84WQvxyOL1JS93jdroqpL5wqbrpowKVjAn2qZhV8fAZoK55AxnD2UbJTy3zZSGoBx1P2zj3m3TxC2tczg/D4uLMc78crNQn9EZ+ciJUn3mQ+zG9s1Kuim9Rr4t17s39DEr2HyTDbQmkBNtY8m7PhQJa2tNa/AR827Ygk19luqnzTQF5eXN7IehVP7SGPEh8CyTtPoI5vShg3C9Nj1ornBfMjFYSiaOPoKzbcZgzXPeUvbfBzdlOI7jzBgbUmPONeO2xVUtNQc4484/ozHnHTcrxY4p7hM3H63R6BgsnUMz/Xz0RaqzxUWWLVSb8phWY85NGfnEnGGh52TYcl2NPMvBrUzTEdhl8rCk1ogdY5y2F7Pw0siBtjxzchNFerH3FOJC06yRq2f5OPna+Pk02sP6ck5lsOVUE57KYzvGc1NGzRRgtdpccNlHxrYG8ng13xu1lWZqS2vB9GOU3ZRRobY0lGWZjbLFRgqZUKE5ZChv3q6mUxvrUS4vs13eKn/RE/ZhtyMxrPnCCmC7zNbi4vKhc12CuTSSpTSsjuEkTyug84kza01u781H3Nh4mPB8ZEbCmjjskLr8EttszKXJDcmLX0mYWYFl+xVyBaSPvTl/5qzAG2XbVlmxQjq3MdeG241HghSOlaZrj2hXQux5XTPz+k5KybdLabopw3EcZx2zoTTmrqnXuQOe0vH0X3LnWRqLWfpKtnWWdHW+1M5JZbHN/FzLrHaEDJp1U1rRJe8UTB2BbaaMVCagMTnB1t2YpC6ZmQpdE2lWiznqTdZcA24zc5TOtdWSa6yJBZN/rqG1mUys9mnzyBsXNk87cmB8b/IC9SDvZO5DbiYyjaReGjM0vJs2xgv3FWSladt5nnmHee0yc225bfr2tGx4wdxmk7W0+XfI7cqJXCjbpmeJvktLrYT8ZR03dQfltQ7TdmkiyZyxQyebdM2UYYXpyL6JmQAu2i6zD8M0E2empa/PklLTtfZStjnIH1F+4UsCo/QxsNhRHblCkJuySqatAWb0ywqewdIptm76bpeO2bRzO3LNZFKj72STfFalndRSm3zS5ZEu3+7CTRmO4zgzxobSmA879vjW442xtRq2NX7+loawFD+LS0swTNujZodQV/MWmk3KeSZf8vm5oIHOxYMDmYyEkIpKp9rsVFtaMuU05bdjrFMZIPrfnZvkubAA8wuT7WSyWFgIPwidgKlcA5mcC/XWyEgn5dJRqL9U7kZ9Lk0m6dhrGdrrivdiFK9tODT3adRs4qftfJSMZtslD2hWEx5m/zUNrzS+2I6vtduLhThdpoyaf+DGqu9M7q+dtmy3bZ5zg8nzNT8f7mmt9abm/RiPTFoKE0xS3exmMkY5304TUZ4123toTj6pmTLya1lg0oG6aLbnmdStDc/r32Lr1ppZ7DMwoqm158+AvTells3hmzcXQsuI6t5aQ9YHIrIxLtRxnJlGVTstZxtGMDvtiMgzqnrgvi7HLON11I7Xz+rhNmbHcZwZwwWz4zjOjOGC2Ulcv68LsA7wOmrH62eVcMG8jhGRTSJyu4jsFJF7ROSiGP4yEfmqiDwoIteKyGIM3y/uPxSPbzbJ/beI7BCRN5n0z4xpPCgiZ5rwE0Tk7pjOZSKhT19ErhaRrWtx7X0RkaNE5Asicl+so/Ni+LEiclu8jn8RkYNi+KKIfDyG78yu56Oxjj5i0p+qrkVkq4hcvUaX38kK6mdBRK6J4feJyB+b5P5VRO4QkQ+Y9GvPyiEisi3W2zYROTiGnyUiH167GphNXDCvb54FTlbVY4HjgNNFZAvwp8ClqvoK4EngfTH++4AnVfXlwKUxHiLyE/H4zwHnxrBDgAuBk4DXAhemlwe4AvgN4BXxd/pzeZF7yRD4fVV9NbAFOFdEjgGuAi5Q1Z8C/hH4YIz/foAYfirwZyLjwWO/DbwRmDN1NlVdzyDT1s8ZwH4x/ATgN80H/leA1wBbROR5Maz2rFwAbI/1tj3uOxEXzOsYDfxP3E3DOhU4Gbguhl8DvD1uvy3uE4+fEjWY5AfHTiL7eWCbqj6hqk8C2wiC/wjgIFW9TcOQnk+Y9L9P02f6PkdVH1fVO+L208B9wEuAVwFfitG2Ae+M28cQBAWqugt4CjgxHhswGVosse6mrevdhHqaCVZQPwocKCLzwP6E6/lBPCYmjnQ8K7Z+bL39EEjP9IbFBfM6R0TmROROYBfhBfom8JSqprHwjxFeNOL/fwHE498HDlXVe4ADgC8TNJxG3Cydl8TtPBxVPU9Vb13VC1xFomb3M8BXgW8AvxgPnQEcFbd3Am8TkXkReRlBK0zHrgJuBQaqeh9wKNPX9a2qet5zcoF7Sc/6uQ54Bngc+E/gElV9Ih67HtgB7IhCvvqsAIer6uMQPg7Ai+L2tap6yWpf23pjQ838+1FEVZeA40TkhYQm56tL0eJ/1Yuhqv5uFl6L2+UJcSaJTetPAx9Q1R+IyDnAZSLyIeAzTDT9jxHqcAfwKEEQDwFU9UbgRptsIavOup5Fpqif1xImw70YOBi4RURuUtWHVfUaJlowrLM6mCVcMP+IoKpPicjNBDvhC0VkPmpqRwLfjtEeI2g+j8Wm6AuAJ0rpxbhbzf6RwM0x/Mgs/NvMMCKyQBA6n1TV6wFU9X7gtHj8lcAvxPAhcL4591bgwUrS32N16nqfMk39AO8BblDVPcAuEfkKwdTzcCHptmflOyJyhKo+Hk0eu1b5stY1bspYx4jIYVFTRkT2B95MsBF+AXhXjHYm8M9x+zNxn3j881qf+nkjcJqIHBw7/U4DbozNzqdFZEu0mf6aSX/miGX8KHCfqv65CX9R/B8AfwJcGfcPEJED4/apwFBV7y2lHetuNep6nzFt/RDMFydL4ECCInB/Ke2OZ8XWj603B0BV/bdOf8BPA18H7iLYBD8Uw48GbgceAj5F6EUH2BT3H4rHj+5I/5wY9yHgbBN+Yszvm8DlxKn9s/gD3kBoPt8F3Bl/bwXOAx6Iv4vTNQCbgf8gfOBuAn68I/1Vqet1VD/Pi9d1D3Av8MGO9IvPCsE+v53QGtkOHLKv62KWfu4rw3EcZ8ZwU4bjOM6M4YLZcRxnxnDB7DiOM2O4YHYcx5kxXDA7juPMGC6YHcdxZgwXzM6aIyKHi8jfi8jDIvK16F7yHR3nbBaRb6wwv7NE5MVm/6roQa3PuVtF5LMrybcvcXZhusb3rOD8s0Tk8tUvmbOvcMHsrClxBtg/AV9S1aNV9QSCu8gj28/cK84i+HYAQFV/XSuz+fYFqvr6uLmZMOXZ2eC4YHbWmpOB3aqapviiqo+q6l/BWGu8JTpcv0NEXp8n0BZHRP5QJk7uLxaRdxFmn31SRO4Ukf1F5GYROTHGPz2msVNEtve9CBE5RUS+HvP6mIjsF8O/JSIXxTTvTn6b4/T5bTH8b0XkURH5sXgsubm8GHhjLOf5uSYsIp+V6LhfRM4WkQdE5IvAz5o4h4nIp0Xk3+NvfMxZP7hgdtaanwTuaDm+CzhVVY8H3g1c1jeOiLyF4Nf3JA2LB3xEVa8jeIr7VVU9TlV/mBIRkcOAvwPeGeOf0ecCRGQTcDXwbg0O4+cJTvQT34tluwL4gxh2IcFfxvEEL4AvLSR9AXBLLOelLfkfAVxEEMinEnxIJ/6S4Lj/NQQfylf1uSZntnDvcs4+RUT+muCvYXcUJgvA5SJyHMG95CsLp9XivBn4uKr+L4BO/ATX2EIwqTzSM37iVcAjqvpA3L+GsPLLX8T9tPbd14BfittvAN4R87lBRJ7smVeJk4CbVfW7ACJyLc06OCZYjAA4SESer8E/srNOcMHsrDX3MFkNA1U9Nzbpd8Sg84HvAMcSWnT/V0ijFkeYzt/vtPHteW08G/+XmLxjXeeUGNJs1W4y27VyD4DX2ZaBs/5wU4az1nwe2CQitul/gNl+AfC4qo6A9xKWvcqpxfkccI6IHADjdQsBngaeX0jnNuBNElYqsfG7uB/YLCIvj/vvBb7Ycc6XgV+O+ZxGcDKfk5fzW4RFEAYichTBST2EFUa2isihEnwpWxPM54DfSTuxVeGsM1wwO2uKBneGbycIxEdE5HaCKeCPYpS/Ac4UkX8jNM+fKSRTjKOqNxD8/O6QsNxWsu9eDVyZOv9MWb5LWCj0ehHZCVxbKfYpIvJY+hGWXzob+JSI3E1YA/DKyrmJiwj+re8A3kJYmik3L9wFDGNH5PnAV4BHgLuBS4i2eQ1+jj9M+LDcRNNm/3vAiSJyl4jcC/xWR7mcGcTdfjrOGhBHbSyp6lBEXgdcoaquzTpF3MbsOGvDS4F/kLAiyG7g/fu4PM4M4xqz4zjOjOE2ZsdxnBnDBbPjOM6M4YLZcRxnxnDB7DiOM2O4YHYcx5kxXDA7juPMGP8Pc6+CBOT7GdsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sh.set_full_model(bkg)\n", "bkg.ampl = 1\n", "sh.freeze(bkg)\n", "\n", "resid = Map.read(filecounts)\n", "resid.data = sh.get_data_image().y - sh.get_model_image().y\n", "resid_smooth = resid.smooth(width=6)\n", "resid_smooth.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Find and fit the brightest source\n", "We then find the position of the maximum in the (smoothed) residuals map, and fit a (symmetrical) Gaussian source with that initial position:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "yp, xp = np.unravel_index(\n", " np.nanargmax(resid_smooth.data), resid_smooth.data.shape\n", ")\n", "ampl = resid_smooth.get_by_pix((xp, yp))[0]\n", "\n", "sh.set_full_model(\n", " bkg + psf(sh.gauss2d.g0) * expo\n", ") # creates g0 as a gauss2d instance\n", "g0.xpos, g0.ypos = xp, yp\n", "sh.freeze(g0.xpos, g0.ypos) # fix the position in the initial fitting step\n", "\n", "expo.ampl = (\n", " 1e-9\n", ") # fix exposure amplitude so that typical exposure is of order unity\n", "sh.freeze(expo)\n", "sh.thaw(g0.fwhm, g0.ampl) # in case frozen in a previous iteration\n", "\n", "g0.fwhm = 10 # give some reasonable initial values\n", "g0.ampl = ampl" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47566.1\n", "Final fit statistic = 47440.7 at function evaluation 228\n", "Data points = 30000\n", "Degrees of freedom = 29998\n", "Change in statistic = 125.402\n", " g0.fwhm 20.0951 \n", " g0.ampl 0.111272 \n", "CPU times: user 1.58 s, sys: 17 ms, total: 1.6 s\n", "Wall time: 1.62 s\n" ] } ], "source": [ "%%time\n", "sh.fit()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit all parameters of this Gaussian component, fix them and re-compute the residuals map." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47440.7\n", "Final fit statistic = 47429.2 at function evaluation 408\n", "Data points = 30000\n", "Degrees of freedom = 29996\n", "Change in statistic = 11.4937\n", " g0.fwhm 21.5373 \n", " g0.xpos 66.3947 \n", " g0.ypos 69.3172 \n", " g0.ampl 0.107149 \n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAEJCAYAAACnhI2ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztvXvUbUlVH/qb+/vO6YZW5CEo0CgwbEXwIvIOGkNAENAhROBKMEADF0xGExXzEJMMAY0DVJCrwoVceTWot3lGEUHCQ1TkZdvyCEKkw0M6dGiV96u/s/eu+8datfdcc885a9baa39nH776jXHOt3atWlW1alXNmvWbs6oopYSGhoaGhv3B7GwXoKGhoaFhiCaYGxoaGvYMTTA3NDQ07BmaYG5oaGjYMzTB3NDQ0LBnaIK5oaGhYc/QBHMDAICIXnq2y7DvaHXko9XPdKDmx9wAAET0pZTSBWe7HPuMVkc+Wv1Mh6YxNzQ0NOwZTozGTEQn40UbGhr2GiklKsU5PI6C7At+99u7vzMxT+C/5T0Ly6UexsMXSyAtN59ZLu1nrGsJrcza33x9eDgMPzwcxqH+3oFIV74PACzm6/D5HLj2qAsDgKOjLizf49fW+/By8rLRbF0eHkeDVrcLp/4OjDx5PqX8tW/F80zG9/Yg71ttRebj5aG1Fe/d5PfQ0slpZPB2I9tavj443GxfHPzdcnvi78zbUH5/rX9Z1zyt5TLeP7UylvLM4N/p9o+8GA984Ys3IyloVEZDQ0PDnuFEacwZy6WtfeVRL6o582eOCzVlsxB5z6j2yWGlx8O9vMe8m6x/raxpOdTwvDRKbSMariGn7c0eIukdzOxvwtOorc/SzCQCr3/xfCL3d9m3DmbAor9OI/Ph5dTauDdD8HAiBTMwbLhapzxuYash8iE1WmbMoJIfWWDYmDhK0z6vfJ7gi5ZZCpsolaD9XkDHQKgp8cYIrVIdWNRD7uyR9inDrUGPU1YaXaGVbwpFYApEB60cF/DbXWRw1Opfe876ThxUZJbZ8/GoDQ0NDQ3HgROlMS+U6YU1BRljBDxOLTui0XANqQYWFZAR1ZJpBiBQn/JdSmXWDDCWpmxplmmpa80b8fu/eSaxzTf2tC0t/xzfup/b8QLrd/Y0Z69eNWOdfL4GY6mUKDQqp3ZmIhGdLZdmLeZzFRrziRLMGQvRuI+byvCmp5HpmtVxtukEktIAuo5cI5A1aGXSOq0Wz6JUZDpRoSzzznEXKH93jdYoQdZfztcqk7zvURnyuSx0pYDm9wDdE0Pzwii1JSnkF8vy95oCkf4ZoRUk+OBWk08NGpXR0NDQcA7jRGrMHHKkrx1pS+CajDXVkoaGWqOf5YMahZWnp31GNAlu9YaSh/YOEeolqinXzDxCs4GlrYla4HEsy/8YKsP6xppWr3kGWBqzTHsMFWYhLVGtCpZmkHKmwJ/LiNIque5miM1SvDRVg+qBnz/HiRLM/OPxBiepjYxaIS0bEefB5LRWc7MpURkR67nsVBE6QJs+89/yOjqNnLH3L5XdWqwgy6pBCuWIkB18W9TTFLOZLmg9Qebd0ygI2R4i30mrO03oavQFv55KIMtyLir7VI1CYEHra9H0SrRIqV0D67qsUZgaldHQ0NCwZzhRGjNHYposh+a5UQNJV6ws5wX6gl8XrbsOjQHo2k7NlGtKwyf3HFjlKeJY2rJXD9IiX7v4Y0MTyuka+XnfhYd7C1lkm+LvoFEemnY3Zope0phlnKi2bPWhbT2dxiBSZk5ReH1tgXV7KPXHmlnfQaMy4iitCCvB4wd52IBvhf7ho/yyzNdrHFanrEWNC5HJfxodWTbiaPk0KiHqvmgKD8Rpjcg7ewO8dY/nnwJCrsYuUUMpeeB1n/sQV2qscm7jucG/Z3QlKgfN4gKX882lNDMOlP41uG5eGQ0NDQ3nLk68xgzoUzLPLzNKC5hTJaaly+kVoGt6ntEhQgV4U9cxhs3aeDXT7ePGxqylDx/McJz3HzPL8crCv6G2eESWWT5v5VkyFEdgeZZYM09ZTu477mGbhVueVh7VhGX8Ul5W/xqENSpDR42XRXFaUilILOs0n15p8TjGCGN+PeZdIp3iYAbMK63tFjyruRy4NFcpi78v5WMh4kVgtatSGyohC2mPn4+gJDC2FcryvjegaX1A2mHGbiikwaJSOI88RdrSy0Wr8xoa50QJZo5ag8RUAlpDtONpHUgK5Ug5Na5RS1v6XnsaGp8BaKjhzzmsWYRckmvF2+b7aFpg7WzA+ya18GZakfY8Nm9LUEYUHdNnH+t3kQKZuzzWaMslw7EMiwjn0uCuasYsfNB+Gsfc0NDQcO7ixGrMEXi87tRac4bHadUsHrDu5euSZ4QsC9fQrPzzM1OO9p5bYQZfvBPR3Lx4Y6b4UkMq8Y5e/hKRBRqSI9fK58HzTJLasueWmPPhG0RJKmFQLr9YG89FtWevnuXMopZzttK3OPu2H/OOUJqqyHseSq5OUXeiUmePcIqei5q2iiwqcGtXdnnw6BMprCJL3yWs8nn+qFZH9FAaKEvQ2o3cdIenW7tiFUZaVjlK9z1hHIFHZYxxk+Ow6sNr35FBnCs7ljBumxg1NDQ0nMM4URrzNs75/Lmo5hM1ytRijMbsLTzxyqctkPFe5WC2jrBgrmeeEcXK0/tOnpdLRrTKIyv1St/POhnEMxDVIryKMVDXNe1+ypWgkXTHGP+8GaBLp4hZCNC1G48a4s/WeEe1/ZgLGOshwBHlIvMtiwccI7xruO+alUleGIc37Ts8FCu02L2c/1h3qNp6ivqgWijx9yWaSxPMJd66plla59VJH+AI7z6GfvC+R9RfOT9fWqUJbE9j5LxymjLfyPYBQL2r6qhybvd4Q0NDQ8PUOJEaMzAdxVBjubcWS0S0lQj1oGlwGdvum7FBHeR0Cs+pU0JGdWiIatbbfMMSBWWFaRqz5YkhNebINqyWUW8b1PobRzElHSPvTVUez3Bn9T25XS9HlDLR4tQY/06UYLY6sub0vktI+iJqRY80vAjFUkpnzBTXStOz1GvPWDTIVJ8lYmG3wj0r/OHhOlwKZo1/9mANehFhFfWcyOU7DtS6+G0roMciUh9R20yjMhoaGhq+xnCiNOYMOVpvM7ptu21oCbsakSN5RhZllPybLYMnh2cMXRl+sElxRCznWlnGhHM6AtC9L7RrSXNYeWqziHnhvbalPHbhMVSTp8z/OLVjrQzbzFxLz9TiRArmKcA/qLY7XclS7mHshx/LF0fgNVxJOYwRxlrZLIE9dnVhtF4tOsjal8S6zkI5wudr70wzACPpi7NBWYzNp1R+D7VHp3moUUi88CnQqIyGhoaGPcOJ0pi9aVTN8/w5uUeAFd8zcEQwZnSWNIu3cKPGaCRhWbEPZvZ025vKa9cbS3N5Wm7J7Ty1OvWs8JKK8OiL/GwkzzFT+OjRWttirFZ4HEb0bRH1gIremxLHKpiJ6IUAfhjANSml7+rDbgbgtwF8AcCPp5S+SETnAXgJgDsB+AcAP5ZS+lgf/+cAPBZd3/zJlNIb+vCPpZRuuU35lst1J9cc42u4sG0F8Rhogpc7+h8o5a/h1rRn8rU5pQykxw8lsPa6cJ+PZbPxbpFtVOXzGk8sw0v5RpGUQWnKAX1XU/QphLLWBuRWr9H9jT3PoDFlGouU4nGPW2N+MYBnoxO6GT8J4F8DuDWAfwHgeegE72dSSt9GRA8D8MsAfoyIbgvgYQBuB+BmAN5ERN+eUqo6fd7yX/QqPup/uU2cXIZtjTKaoOU7fvEz2WS8Me+ytWsQLw/rfDXJynPyrHLVCGNP27WEcclX2VuGnP8uDEFcI0xq7RGRbzhW4G5zzp9UHKY6M7C2DMee5/gs65FS+lMAnxbBB+hMHEust5J+IIBL++tXArg3EVEffllK6dqU0kcBXAngrn28v9tl2RsaGhqOC/vAMT8bwEsBfA7Aw/uwmwP4BACklOZE9DkAN+rD38mevaoPQ0rpLmMyt6y2XMMcizHTzm3d96wycO+RmvfaxpUpOlW2PDFKFIU2s+BbgHqIruCKXnt5lhY28evFcpPCiNS995228eaZWkvW8irx7mPa4JSugNvao8bgrAvmlNLHAXy/CNYWLyYnfDLIj6AdyR55XoPFj0lDYqlRbs2V8TJh893GpF/DV9Y08MjSb2swkwYyTxhrYRaf7D3vlc2jJfjxSlwgy2e0Zew19NcuKIuazaCi5bCoxtr2PqWiU0qrSFPtMcccxVUAbgHgKiI6BPAN6CiQHJ5xIYBPWokQ0Zd2WciGhoaGWnC5lFK6QIuzr4L5NQAeBeAdAB4C4C0ppURErwHwu0T0a+iMfxcBeLeVCH9pIiqOV3JE5O5mUlOJbEHorczy6ISSYTJqiNTc4uS1pGysmUGNdb9WQ4poNJzW0IyPWpk9Ta5k2KspGzBcZOS1De1eUjTj/FdqzRbOhidF7Wk72m8rrOSa6i1K8SijWoxNy4tnCWOO43aX+/8A3BPANxLRVQCenFJ6gRL1BQBeSkRXotOUHwYAKaUPENHLAfw1gDmAS2o9MqLIjaB22Wu4EbBGtJMXwGbjtnyF5QCkdbixU8KokNvG48Bz/ys9WyqXB889kcNaWq1xyBp9Yb3PPvsK19BE3vOer3GJOojkU0LUAylSlr11l0sp/fNgvK8CeKhx75cA/NKU5WpoaGjYJ+wrlXHWsK2hwUvPQr5FipagTeOi+UeoDDkzkFPxA5FG6V3OdYx5tzHeApHrUrrb+PQC258Iwg3jHFGayAuzZnfb9MkxXhVRKiNi9K8pexPMsHnc6LQl8rFkeqq3BewVemMQpTI4IrvlyXeJ0hzbcoBjuFeOMYsqrDJEnvfaSeQaKHuWyHJE6aIlo6xkHh59V0ufeGWL2B1Kfas0uGiKRQ2iz0xx7BXH17Du09DQ0HBuomnMDqLTFo/CqEXWmqdAVCP0NF5rulpKL+cf9SqxDGHcv7eU3y7g0TeRfGu1ZE/zsrYdtcpXU7aMkrGb3y8ZjD142rNFyY35zmP21qjRfqucA/bV+LdPmJK+iFqHS25oqzL1vznfG1ngolEEHq9cA9n5vHcb05FK9cmpC43KqK2bWtTUmdVWIsJYdnRJHfCFLqUFMFpeFhVTu1udZpc4CLSvKb1ivHrLIOMdpZC2hPGUCkHNApNGZTQ0NDTsGU6sxmzBGwVLlEX0WWvaJjVNPmp69IZlLKxZqsuhncgClL01PJQ0Tk2rlPtGWN4KY6a4YwyBYzxjvN8RmgbwT+CW1zWopeqs2WU2GNfSXp5WbXmpyLKNoV4An65ISt4cY6mzvfVj3kfUehRovyMcoWysUe41wyqittLMmlJ6HUy753loRFe31QpwTxhLwWzV+7Z85zbxoi5dER5Tq3tLMFtHYC2WdtvZRgnxFilpcWX4touNeDspPSMHEAtWejWUmYcaPvpECWarUUS0WSt+1HDjcXDSWGbxwhpkO7M0l5oBqLSSrQbaqjfZOayTXyzBnH9b9e1pbmMNZmPiykFPpsNPfbE0O+vMQC6YrQNfV/nkZ1hYpG158GwW0dWjY2ZzWpratZWPrGetfXv2i9q6GvtujWNuaGho2DOcKI05wxpZ5T3vuaimHCmLpWHVrnTyppTWe0otnd8rlXssahZBWBqzZZH3pqoRbXlXWnSG9j2sPaSztlyiL6IncS8AQGlbM4UOq32XJaNMopTGNhq7VRatbBa8Mzp5WK1Nw+pLjWM2UOJ1o76qUzQkC9EGa8Xh0zW+itDLIzoARI1FJdQY0LhgjrhHlXjEbXx/pxTmFl8r41qCWdIX8r2sAcCiyXJac/HcmLYujzHL8DhqCc2OkN8h51FbNo/O8xSCUpkAe9AZi0ZlNDQ0NOwZTpTGnKFpmVr4lBizck4i7GanGJhkOrV0iWcA2VZj9miU/FvTlkuznZrwbTwzLA39YLZ56rf13SxDmqUxS/rCeh+e38FsqBFrbYBmnbYWaR/ROrNmqtY7W21gl/3TSz+S/9SGzRMlmDUh5R2rNAUiAjn6Ea0plNWJZtCFRpQuiXLqtZ4bnBv1ymBRFx695AlibcpfghXHck/bTGDoJTJTOHNvWm/xyjXg3zo/P4NdF6WzFmsgvZF4fmNtKVOixCtbVErGrnZfbFRGQ0NDw57hRGnM1gYsGWO1510taIjAmxIC407GLnk/eLSGh1XZ+nLmafUYp3/+e0x91mrOsozaYbpWWlxLBmt3PNwrQ8ngp2nZlja+MNq3ZTyUPsnb+rVnlNptqV2VKJfaNhGdOXrw6hMASDtK2sCJEsz8I0Z2ybI6zNipVy0PGvHK4PEsi39GZEMkKYhLVmmPb/Z43AGt4AimiFdM1AruURkRiqB2d7fcTjSOWQq/yB7IB7Pye5ba5moPZgyny4eHm2VU01MUmhK0djeWyrD2kM7peCgJf96+tfJE37l20FfTGPdYQ0NDQ8OucCI1Zg6uPddsxOL5XnpxIunsAt6J3/l+hqY5WFqyZhyJGHv4/ZJPq5ZPFN503/Om8NLK15pmFKUVNt7H+TalpeUcNZqntlRbltNKc1ZoAzWomfHwvDKdYZWDPzcWkdmHlodKZVSU40QJ5gxLgEaOVfJQO5WqSS/Kn3luO9Z0WXbCkmAuTb2t/THkXiHb8oBj6lPrJN4g7PGGEcHMwzzBEaWtIuGRdJZLndaQ9EuOK6+39arQ8qh5n4PZ5gBQ8/wYF7ycZ8m+YLWNxjEb0ITbGIPBcWm50Twt4RFxi5MNVApjbWvKGuPfqrEiZuyrQcRFTrvO8DQeK61i5/O+k59FCFO1PS5MD2YYaOwbh6sGvr3Fy+7aDY7X9zaGyeig6QnlUtuYVQjmxjE3NDQ07BlOpMYc4TQ1jNFWjlO7LuXlOdNbWnKtV4ZVDq69z7D9uYZjZjecV7Y8HErat0ZZSL7aOrZIcqIZrhdEoGw14DOoyPtzFz9AuNsZdID1DiXtOdInefklnZG/wdTnQ1p15mnMWhw0KsNHjbGhdH2cGJOvJ4y135GN6q105T1vAIwI5xouncMzlmn7FnuGGy88yld7Ll4yrVqudSzG0nYarZF9si3edowSpJVBa3ORupXlsH6X8vQEryWYx9J2jcpoaGho2DOcSI25BvuiNU+hLXsGvuhxTlq6Y5Ffab6s0yxqtD2uxXgnTkfzsp6pXf0ZwVTas+Vhod0vPWtpr5H9NaLtxvMK4eWqSW9sWTy6x6IygM39YKy0LIQEMxF9K4CLUkpvIqLrADhMKX0hns0grfsB+HUABwCen1J6OhHdDsDzAfwNgEenlJZE9IsAHghgCeAaABenlD5JRNQ//wAAX+7DryCiWwJ4cUrpnmPKVcKYqe/U+UZh+c1q11ko59/zuU1feP7P23pa5Mc9eiMy3SfHOm7tyFZDa9V2vjHCQ15P2b5qaRoZZrWt2RI4M+/vOflvI0ynSCdCZ0QHqho7A1DnLlf85ET0OACvBPBf+qALAfxePItBWgcAngPg/gBuC+CfE9FtAfwMgB8BcDmA+/bRfzWldPuU0h0AvBbAz/fh9wdwUf/v8QCeO6YsDQ0NDfuKiMZ8CYC7AngXAKSUPkxENxmZ310BXJlS+ggAENFl6LTiAwAJnXZMfT6fZ89d0N9HH/8lKaUE4J1EdH0iuik6ZevTkUJ4I6LcJMYi+rfJo4SS9ilPQeaQ+/9qkJ4XlpZs+THLtGYz/T4/QaVkaFlpYihvMFXjTZFBs00tuZaaGmvU0bxcavx+tzWeleBpyzV5agtWNI8NLTyKDQqF/Y5uDVuDaN1Y2vLYbxYRzNemlI6o18OJ6BBrIVmLmwP4BPt9FYC7oaMm/hDAhwE8L98kol8C8EgAnwPwT500bp5SuhzAj3qZl6yppSOJxjTY2il+ZDmwFb5cdp0issJPcspcYOQpqYxnwep88oihMK/Hrq0NpixYdIV3fl6oTErcUgesmW5bVIbF8daUs5Su9nyU8rDuySOqtDLUYoxQLrXb2vqMCGPtuVpEHv8TIvoPAK5DRPcB8AoAfzAyP41lSSmlv0op3S2l9C9SSgt24z+mlG4B4HcAPMFLY2R5GhoaGvYOEY35SQAeC+D9AH4CwOvQGerG4CoAt2C/LwTwycBzv4tOo35yTRpE9CUt3NKWNG3L0rCm0oRluWp+Z2gWcgCD/WoH95Updf59Zm7su1wqvJInpzhKtIY2xeXPl2YzpHwnLXysxpwx1XJyIDbdlvWxDa1R433Cob2zNivj35CflBL1qoho097pOSVN2ZrZaWXTrr09ub30JLhcSildoMUpCuaU0hLAb/X/tsVfALiIiG4F4H8BeBiAh2sRieiilNKH+58/AuBD/fVrADyh56fvBuBzKaWrjbKvXpqIkuXmojU8rSOv0nKmMFo6JdRyfNb0diVgeypigeE9DVIwD3hQHs8rv4jH+UWL1rDKJDtsRJjI+uLhls3A63AWFTQFj2i5H1reL97+4FxI74p/9mDRFLOZviF/DR3GUdoNUT4/RuDXcMfb1rUljDlMwUxE74dDEaSUbj+iQHMiegKAN6Az+L0wpfQBI/rTieg70PX1jwP4l33469C5yl2Jzl3u0bXlaGhoaNhneBrzD/d/L+n/vrT/++PoBOIopJReh064luI92AhPrExVsLQl7bgeGS/HBeLa0ph7Vv7WKC6XyM7n69MoMAdSgf7I2tqC/85xxDPcx/iAp8XLz35LWmOVjhIW1ZKtONo3kT7NEhYtoe1YNiWF4UFqh5b2LLGt54b3/Jg6yGXVtoC18uawllhHjhrzwiKz0+jMKEotjoEpmFNKHwcAIvrelNL3sltPIqI/B/AL22d/vLAqXuMd5V8pwGWa1u9SuJdGpIFoU+D5fB2PlDJrsKbv1oIPT0jPlGueRy5TrTCWcUsCWG5UxL9vVMhsuz93pJPygRHY/BZkfO9d0BeZLtEEtUcfRBCtj8gp7FG6otTvNIXHO0Isguh7eogY/y4gou9LKb0NAIjoHuj8is85HLK31YSrZhzSdiOLaLUlbliDp5F7aa86kdOAF7PhuW4La7THWuPlwldyWtw1ZoGhcFbTFfxo1P2NwxKy1uBqPRuFdi7klNA6ryaENAOqTGdKIS3tFhxauOaTrWm8JeNchEu2EBGGmgJkccilvpwxhXasISKYHwvghUT0Df3vzwJ4zG6K09DQ0NAQ8cr4SwDfTUTXA0Appc/tvli7QUm7lRqzd1x8SZMtudZEy1jSBDmypjNXysZd5Kzya7CsvzycwDwuIDhmdr16dlmvgWrTy6iLXC0srfW4eOYMj5/n2JZXjkJqpZ6nTy0i9EUJNTMl3rdzvOiMV6OTIvWSKlZbFAUzEf28+N1nkn4hns1+oEYw5w+n3ZOUiLX5+hhuSv7WBI4l5BfLzWdyIznv9DoeX923IeSY8S+KhDW1wWkNyTGPMVCR8Z0sRNOWgjYiCI5LOFv86Qzl/Eu0RkmYWhRGKS2NzojmOQbaO9YcrivbUk2/ramjmngcESqDL9I4H523xgfrs2poaGhoiCBCZTyT/yaiZ6Bb5HHOwdKSeZjUfrN2bGmvpw5tumHsFNOyFFsaOx+R53OmDR/p6fNFJPl69Xu2NiJKjwxt4M9FzbM0SWvwZ2urw6J18m/LWBNdCBDZU0QrCzA0Xlpaquc5IO9F9iTh8SxvmdweasvjhfF7lrbI25Dc33tKeBRhicrQZqTesxZq3mns+4/ZKP+6AG49LruziwgvzDv1oRC6B31tHQiBza8l3RDpIJJfK/Hah4f2tO30EjjqBTIfNGT+0o+Zd6osmA+wFgAWPVYSuBrHXEJJyOZ6tigoD2M7CqcTcj75u3kbNHEhZUF+f6+MUTtD6T0jnGjpORle2kWvJp8ILWEpQiXKspQHh+c37XHuO3eXEysADwDcGMAvbpdtQ0NDQ4OFiMb8w+x6DuBTKaW5FXmfoY2opcUm+ffBYaeBAkNN+vBwrTFLDdvSbqJGE0tLtuiTnFYuz3xeLo/UmJdMw+N+0XNgsHBkkAa7xw2BNbC0Xy080xgRY+CgnI5hLLqfgmaIs3yyx/jylsDz1+iMGsPfVBqe3LM7YgjUFkZZsPz6rRlldDGYVbbBNZsZzbBdPdU8GxHM/zml9AgeQEQvlWHnAg6Ut7WmRlnI5mdOCQGcw887PQyXnDBvDLxRatO+fG8wlVfoEz4YSLcfYL3y7+hILxtH7kQLrVPN10tqCetp0wwxz43SwpMoz5/DpcCOdDrp2jSlhwDPUhOSkme1BKGc+ss4fOCNDEIyXolLjgpJj5KxlA1PCYmcbG15Osn+KenFiILEyxcpvzzXMNqWIu8pEdE1bsd/9Bvl36k+q4aGhoaGCLzd5X4OQN4gPx/zRACOAPy/x1C2yeH5M+ZrPupK+iLfO326+5fD+bU2akssl2utVo7OctMajTI5fdqmT3h6R0drQ6BXFm+fhmU2JAI4k+MgNqJLbVmbauZwSVNo4StaqUJbLmmO2jNRDDRaod3y+5YhzNKkLMphCs+BkqZcsx+GVU7PIBbVHrW+qlGIWj84OFw/b23DoH0brU/OGWk7RvMd8wzgb2L0NABPI6KnpZR+blzy+4WIO43FY3IBzAWzJaQ1wcwbQXZpWzinUuc0819rYOCDAc/n2iPgdC9YrY1+lqJDmh1pvha0C8SPjJmJv6twUR5Ozcjpao5jbXofgUYLlOLWpG2lK+szd3TpUhahMqwwqyxRSM8Rj3v2aJEIPeLRId6+6Pla9ofcDywlSlJjvEzydPicz5k5VnuaHx6uv9nBbHO3PMvjRBXIU6z8I6LbpJQ+BOAVRHTHjTxSuiKeTUNDQ0NDFJ7x72cAPB7AM5V7CcC9dlKiHcLSmLXpspwqHRga6/nnD7VnHl9Ox/IoupgDh/0ozBeEcO1ZK4+lsdPpocp80Cdy6vQcZ47Wz3OKgGvKUmPTtIBrgZUWMcOa1tCQNesZmMY82/Qq0bQa6XHCNelTYmagadMSEU15jEEwalTkcS2jr1XnWjolbbjGQ6VUVi/MigPUnVydMZvp+z57GrPsBxrVJ5/nZZjPh1vk8njXriIp7axQd9EZggePynh8f3n/lNJX+T0iOn+5WHhUAAAgAElEQVRcdmcXEcEsp8v5A5861CmL06c74QwAp05vutTl9OS0iXO/Of/5bDit5FysbITn5S9w+jRw+rw+U8Hqzhc4dXRtn89wGaAUCnwv56LAmHecM6C3UU5fWC5MkhPkQtfzcslxLAFk7SNSoie0+9FN60tpW8KYP28JPy5MaiiKErce9VAp0SyAT4VY+Wtlsd6PtxOpqFjKisY9yzJxwTyfd9SfVuZaagsw9pCechMjAG8HIKkMLWzvYRmMPOGhucidPt0JYXl9HruWLYIAHPStYMYchKWmN1sOO6Jl/Fv/OA84fWp9zXvccrES1gezGa6LbnyVgpg3UMl5m8YLx5Nd05Klu182rgJDV0RPS5aCuSSoShpmRBhbGk/kZA6ZT1QYy2cjhk4OGYfPFiIcuyyzVk5t2fW2wizDOuRUcxnlfSIrSCXjeC4fb/dHTCjLAXTBnj+YrV0jLXtMdJbgweOYvxnAzdF5ZXwP1usGroduWXZDQ0NDww7gacw/COBiABcC+DUW/gV0bnTnHA6Nt/W0Oq6xcpqCX+ctNU+dBuM7TpmqC83mOM1UTj7CcuswLzMvC82wpi9On2LXp4EZW3e3TMBsTW9k7vm6y6OBhsx5bjm9k1y0xEyE87JLjVdqO9psQHLPlvuhxytzRDWWKC86xmPDuo6k5WnIlluhl9Y2mixQ5sjzbw6LSrF+pyVW0y4rnvTK8DylPI35SKEvDg91bw2LEtGgzSIm2Y85pXQpgEuJ6MEppVfFk9xfaB3Z86f1lgRzI9/qgx0edgI5X0shme8tlwMuW/LaGqfIy2JbRA6Go89yuVkGAKfmc5x31GVydNQ14uv008AzR+tGzYX06dO24czbhEkTspoPasn4x+tj7CEEFiJCWROktby1ZSCzYE3pAX9QsoSzNYB6ZeHv7bn4eQOQlW6JU5a/PTpLtidOcfB2Y9Evy+XweW4UjCDk/mhtCaigyDGnlF5FRD+EbgXg+Sz8F+LZNDQ0NDREURTMRPQ8dJzyPwXwfAAPAfDuHZdrJ4is/OMjMvcEION6sAptkACJ4Xa5dmFg8Waz5apcC3bbhVdoy1J0uDYE4vAQ5/UrT649PTSWnDo91JgzxbFcrsM5lktsHAJrzSzkAoAIlWGdFD4FSpqypzF7WuGYU54zohRNjpv/Wk2A52V5anCDVmkTFMso6HlsRKEZJj3KR9J73Htj0IYO1y++mp3Olzg8ZBsUzWHOYr3vob23VTdTe2XcI6V0eyJ6X0rpqUT0TACvjmexP7A45ozotGVrsK+lUQHalHPQuc0vn1tZL4D53Gkw7zuFw8NOMEuuTm7WlIXpgvFuC6fjSH6cUxTSwyWyUsvzSCjxex5ql0Pn97SEzpS8dH7NBexpvUfBbaTHBF6JJ92Gh9bah+fHu3EQAxtAchoRAaVhQHspL02z5Ur5yvGjByzI8gKb7rCqJ0sFlREpwlf6v18mopuhW1twq3gWDQ0NDQ01iAxIryWi6wP4VQBXoFv191s7LdWOUBoJPSZAQtMOhtrrDIN54TINhlFuKEhipLWmR9kgcThfgrI2zC10hwe9GTlgfWHvaNEKGysfFS3A88qw0s3UScT4p6Wby2BNfbXwDE1Ttjbu0SYm3nat2rVVDg7PQ6HmucjzloeENYMYq0GXfHlnzBCXZwa8rmdKWyt5f3Dkewf8xwRuKQvRJuR6gByeKcCxPs0R418+reRVRPRadAbA28Sz2B9YHJsFa6oir/MHOTWb271iuRj4oeXLM+KDzud+J+kfx6ns53N4Cpj3C6SPFIk2ZwJ81aI251Ql75NTh4GNjjCcYkvBfCAEvjYAlPjSnKfsY1rfK/XDEq8sw8d4JWi/c9ms3yuhIj6nxXdqVIZ2ArhHBeW85krdauXcUr5t0CreIMeVIIta4gujFnNgebi+PlA0ityHF+L7avl7nhxaezgzHypbq7QmpjJWSCldm1L6HIBX1Dw3FkR0PyL6H0R0JRE9qQ+7FRG9i4g+TEQvI6LTffhTiOji4yhXQ0NDwy4xllsfc3pQXQZEBwCeA+A+AK4C8BdE9BoATwHwrJTSZb3HyGMBPLc2/VqCn1MJ3hr7UzgaPsgSSv06Xv7MmaN1WtcebU6RuVaXy3ztUefNAQAHs69gAEsVmZ8ZaM9cS5cauvSrztfEtN9BduL5qFeGpjHL/TE0eAeceppc1ENANdws/Y2HrGe8sljvF2mb1mwC2NSq+X4tVj653HmTHstgmJeilwyxVp7aMzMMDZ2WliwNbNxr6JD1qdlMXzgyeI/+2QXrxxq9KN8tibJp9AXX3vlEtcZDZaxgrnD8GI27ArgypfQRACCiywA8EMC9ADy8j3MpOkH9XABfxNpQqSJCZciONWMVzKkEawP6/Pzh4XyDN+NTnbzr29HR8IPyD+l1vtypz5/NQfm1l0vg8AwGmxmtWsgC6Dc0SkfzwcBiruoTwph3kCxI5UpFXiezmX30D810wcxXC/K0+KsAm8f8aKihMUoCV7tnDW6lk5X5O0iM9Qza1qMoP78QVIas/9WgXci/JIS4wF8uMThPb5A/q2e+8GMlmPmJ8P213GBKO1IuU5AapWitfpX0iRwopK0oX692lZxigQkR/QF0AUwAbhTPYjRuDuAT7PdVAO4G4LPsMNir+nhIKT3jGMrU0NDQsHN4GrMn6I5DCGp0iXa259bau6lZMeJ+NhtSDnm0l1OmBZu6DwopRtf8nBydubEg55vTtqbL5/UPnZrPh+onj7xcIh3NV2WW+VvGr9XiF6a9c82Fa9KyzNI3VIZrNInllbGR/pwtyFHKLTFm5y+uCct9q+UU1TKMephtXIwDr6fS4pSSH3NebOL5QwObMxbrW/F40W80W25qqcDmrJX3l69+VU/rlOgOGTndPFs9o/SJnA6f6VozKklrWPvNROHtlfEn8WR2gqsA3IL9vhDA3wK4PhEd9lrzhQA+aSVARC8F8KMy3LMAbyBXMGt4nMPiH+f06aHAknnyD5y5raMj/SPyvPIzebMkLjCusxxOoToKZb56Vpt6Hx2tG3LOnzcqC5og1fqj5hlQcgmzjqjXYPGeJVj0Rf6thUurO48X2dDHg6QPIvFLiyDS0hbOfK8IrRwZVnn4IDxjQm0m7knlwhpAJXXBBbjWbvngnpUb7+R3oFOULE8VKUCPGL0oN/nKz3BqQm6Xq1IhS6ycZnM9ENGXWHFfnVJ6hHyHsRzzceAvAFxERLcC8L8APAwdt3xHdMvCLwPwKAC/byXQv/AjAICIktURZRj/nQUs14w55Co4KbzkznE5Hue2FgHBzP0++YZC3F1PLnXm4GnzhvfVr3YzAF42q/NqmnCE69Ug33PBOmkUpdV7UnjKezy+pQVZ/GKN65ylfeZ4paXmpTqx8hibJl+i7SWbbQGrela4YED/Tht+y7Nh/JlSn2eYxuwdMsz7l9xdjsfh7cMUrIohj3/3hVLO+XytLJzBekqfqyGldMFmqYfYchK1O/Qa8RMAvAHABwG8PKX0AQA/C+BniOhKdFz3C85eKRsaGhqmxz5rzEgpvQ7A60TYR9B5bFSjRmPOo7i2/R/XAA7ZCCpHZ+moz3krTUvWFpdo5/Qtl2sOW2oH/BQWDq7hcY3g2qOhy96ZuV43HNJbhBubram0xpEP3kfE5Xl5ZZH3+LV3yoanMXsUxdK4N1hIUKg/+W5j4kgtk9MCHjQ3NvnM4eFQYy1xv8usLWdpMh967WizK15+WV+c8uD1KvujpBR5P7QWLHFYlIllA9LoLHPW1OeRsEllRBDZXe6NAB6aUvps//sGAC5LKf1gRT57AUswW65O8qgpjnxYo5zCcG5ZUhlySpSv+XRKlnOhCGaeT7eyqbvOhjNresenhFJIcx5NEygk0rXcoyI7ceVOwDvmgnVA7i6nlWW5HL6DRSWUNnPX2kOJorCmrp4wLglNq86nwkowirbtyfCBK1tASHOhCWBFa3Bfds9gq9EZ/B7QpZ/bl+Sx+feQ7pea/YJ/b0D3kQY2V+Zq9JasmyXW77nEWjDXfOIIlfGNWSgDQErpMwBuUpFHQ0NDQ0MFIlTGkoi+JaX0twBARN+K41lgMjlKWvKG1rIEkmLRXS7XG/rwkdaiPvIzWv7y8FOPPuD5DLSVPLrPYhqz1PY4taJRDrwMqzyNId07WVqW5cx88x5HaeWfZkyV1xYVkcNkGb3tGyMGv+jqLmsRx9TgdSjrk0+t8y1Jv1l1Zs0YNGh9gnvSaHmuooq2vjDajJzNaJ4cHuT31FbxZerCqo9V28DQ4FdDYWREBPN/BPA2Isruc98P4PEj8jrr4IJA62D5dwZvMIeH6+dlA89xDg5jy4hlI9D4SQlLMPPpndy/WMKbrlt1YJVFSzdD8+pYoH5J9KIwyGjLYUv74paEsSVwM3ViCX3NTzojfx/uhuYNblPD9QfHujwZvG7kRoVa3cxm5UFW5sGplNJAr/nyZypFUhrA0M6h0RhWmUruj7Kv1CyxrkVkd7k/IqI7Arg7ukUfT0wp/f3uitTQ0NBwsuEtyb5NSulDvVAG1gs5vqWnNq7YffGmxcLQTK2pOx/JB3tCHGJjVeAqDyUsp2WNwpomLyENkVxbWS1UCGjL3nXpHl9UImcWGeZpLKwu+XFVwNCzhdd5iZbhmotFxVhUhPWu1vNy5R/X1qyVfxn5FWgC+oL7e3tpaUZrz3jLy30Iux1YxjILy+WaDuT+zdJX2pulaTPFXP/c+yPSJzwvF35t0VSl2eRU8DTmn0FHWTxTuZcA3GsnJdohLMqgurId1yDLvcv72BatwjGbrRs4EKcvZBmmgjUNlvmp9/owudlM/stXalmuf0BcGAPxBSYyjuSurbQHz7OycqETObOwNFBGzz3UXMTkQFea2mv5r5r9fDhI8LSvdZ7XyhVdpMSfz9elRUZjB0BrcI5gJq4n5ZhTSplHvn9KabASnYjOVx5paGhoaJgAEePf29Etgy6F7T0sbSnDWtyhIRtFFsDAZ9ObKpW05NJ0TtN2rA1ntOct8OelH3JEQ5NWdSv/gRYzX88AFsvh8l6uVVn0TFRLjlA4keclfbLhR83L1v+d9df5FRL/UYEa7U/z27WoIc3nPC/0KVE+0gCtlRfoaKrVjOHQfpeFyDNyUvhyua5OSyuNtnvrGd4eSliVeTncgU3bea0Ej2P+ZnRbal6HiL4H67yuB+C6I/I66+CcIIc1RfGmIKsGVojH40esvvwvhxR+29IXHvfIqQStUUqLOsfBbO0K5VEEUmBoqydl2eR7lbhPmacsr3ZtDZqZxtC8P5aAunhAk8PaAo2xiHxDfi29hmoHdGsPC+2Z1fT/cFj3iUkca7MjDrkLoYdaL4na+rf6Zf7LVziuBjgWt0ZAexrzDwK4GN0Obs/EWjB/HsB/qMijoaGhoaECHsd8KYBLiejBKaVXHWOZdoaIvzCHpQ1L7VVb8ioRWZxQ0kJkGSLxeP4aLSENN4N7iO3vG6lLT4vl+UvjXwYvu7e3ckkjLtFY1h4YWVvm35ArfJFZk1WVlsEtqmHNRD3xZdCz2fAUGYvW8Iy5q+/E6oZrzxm5DvkWoNwP+fBw02BY825TokRbWts0yGc0L6zZbF03h4jt1CcR4ZjvRERvFntl/JuU0n+qyGcv4Ln1SGSBwzksa3FA/nDW1D/HKQljj8+SmwNFGo7sYFojypyktpk4MHT18mB5o1jl8QYgy3IfTU9LWxPM1upPSV8Aa6G8usfKZC2DlSc9LI3B0ULNxurynEaNvtA22cqICOaBnYDROjmtbCdI7BsezNb0BX//FBTMvGweX57zKsHqX2MXi/DyrPavWTLvFZbu1IL5/imlFXWRUvoMET0AwDknmGs5JcvX09ISNQ2zRkv2/Gy7AL2cVkeysNEhmV92SUBb6UV4cVlGLsw1o1DN9/KEsfbbqmupJa+40uVQS15gc59diTzj2kbpG8tDc8GoCY/ZzB5suYDzBgcpZPlMh+epHe4rN8WKvIsWv3Rw71hEBkXr3U6xVcLcJ7zmBOvI6xwQ0Xn5BxFdB8B5TvyGhoaGhi0Q0Zh/G8CbiehF6JSEx6A7nfqcQ4nK8Ebw2eq/TSqDo2b6mdOKLILQflvhEY15UH62YEbyhrUc30ATE1qxhOZlYmnPJUTrxuKnrdnMYIEJhvRFTkJmLYu9NMKjqNGatSm+Nv23Fu9k8G1crWiW9nswG25Xq50beTAbLpiKQtPE5T2JqD0oSgla5eJly7TOmfm6/g4qVObIXhm/QkTvB3BvdNr4L6aU3hDPYn9gVbDGHWvCobQcVhPYY4x8UUNW6b3G+GVLYVlLUVjxtHrS3mkslaHBEsryvkctrcKxSV/k5Dj/fcDC8+8Salwf5cAZSbvWyMa54BruVeOsPTc+69xIbcCW1xmlMxD5PSmgI3aaUno5//lSf+dT3He7wl8udIJJSun1AF4fT7ahoaGhYSwiJ5jcHcBvAvhOAKfRKQFfSildb8dlmxxjRn9tGqRp1vmZqJanGfyytmZp1lraNaN+1uz4SjutPFZ+Hi2haRHeqkSLDvJooghK9a/VkeUZk5Zr7TcfqsnpC81ThC8u8BDVZC1Ysy5tNhWhhraZnUhawqJS+DWnOOQp5xHt16Niao32Hrz40rsp/zw81PtwjYEyojE/G90J1a8AcGcAjwTwbfEs9gfWXrmewC7xVmOpDC+OJZgtIRxtbLloJBrogHJYriOOEZKS+ok8b01d5b1toNEYlmDzllpz+kJzk9NoxPxKM0eYbANJcUg6RlutN8NaIHqukJEpvrfrmhReGiSPbSkEpfSteBZVyekMeQxcLbyTxVeDyZQcMwCklK4kooOU0gLAi4jo7fEsGhoaGhpqEBHMXyai0wDeQ0S/AuBqABfstli7gaaZlTZO4ZrxKq5hoBozVbKm1jJtb8/fakMchotNpGasGTllWpEwT4uwUEOfWPmPma6qVA42vTBKi0p2iYgmLGddfLWePHF6VWeFNqzRaVPNZDSM6Ue16W6TR2kWuAAzSvJvMbHx7xHoaLMnAHgigFsAeHA8i/2EJaQlb6hNgxbKszKOzEc2bu8sPI7andIktPJ74BxxFBEBXlpoUfLKqKFUosK8pj7H9GGtyJbrmreIQqsbayWp5xIqN+r38tQg62jJ2qYXb5ew+PPoQM29T4jJAZ6OSfs51MlqAEXMc0Qi4i738f7yKwCeGk+6oaGhoWEMvG0/3w9nxpZSuv1OSrRn0KY+2skbMo6XVimcT1c1+sIyXNVgF0YowNcKOK3h1dUU080pMDZrXgUH7Hdkj2EJ6zDSfA/opsuWpiy9MvjCj+iBpTI9mb+kTyxsY1zTYG3MlcsE2Jq0V1ZpCLTWNkRmoRENW4OnMf9wPJlzG5GK4x+BUx5RKqOUd005vXBvesWnzvJe6fkSqiiT3hpem9W2QlpO/604ERCGWks2uB+IsDwIAfbublK4RN3aMsd8Zr6pLGiDm6Qv8vUZDMOB8iCirZzkCsU2dEYkvtzUiz8XpTGkssXpRZ52UtKtwaReGYzCaGhoaGg4RpyoBSYWao0fnoHKmobLkXvKaV3Jf9O7pxmAtqE5aqZuJY8N6Xe6jeeFh5BXi/i9gH580IzFrTmxYszijznbdlMe08W/wVxszwlsnrputYdIuTQqY0rvDUv7B3RDfA3FEOn71l7btd1kLxeYEBEB+HUADwDwZQAXp5Su6O89EZ33xy+nlF7WH/b6p+h2sTsE8MqU0pP7uLcCcBmAGwK4AsAjUkpHRPQUAB9LKb24tmy1sDwXpCDR7gPj936NQqMmPE6Rxx/Dg1p0Tk3Dz9EsAVzDN2s8aBQ0w4pcnkHx0jGem7G/xK5ngLo1JP8eUSrJo7J4unyaL09R16gMTWBb8bTySKEslQ7JRXvvUgNJZ0S8JXjcqEJm9WneNrbQZVQc5wKT+wO4qP93NwDPBXA3Ivo6AHcBcFcArwLwMnQnoN8rpfRFIjoF4G1E9PqU0jsB/DKAZ6WULiOi5wF4bJ9WEaUPUSOUcuMbfBxlRB/TAKWbGYCNlVGW8JLplDQk7Xcp3Cqr9jsK/p7WzEQKfG/wq4HmEpXTylqv5+7HwwlD7ZlvVC+FXNTolgXbgTMLs76n5I81gS3LlQ7teBqywN2VAOb5WIO+NARqvDqHF2aV2VMOvPaRn6mRL5GogwUmvXY7ZoHJAwG8JHV4J4DrE9FNsVYwVraUPs4X+5+n+n+p17rvBeCV/b1LATyov/4iOpe+hoaGhnMa0QUmM2y/wOTmAD7Bfl8F4OYppct717zLAfxqvklEBwD+Eh1t8pyU0ruI6BsBfDalNOdpAEBK6RkjyrQxik3lRrYLNy9tlJeag7yOaEiAPfX18o6Ut9ZVSNIaWloatHs1GzypnKVStgzJH2scc9aKLcpA5q2BT9d52fi11KRlmiUqQx4/tVyuTzpJRhviyNqydVqPtWLVQg1tlSq10Ug79GYGsk0PaDNWplo6T6JmgclXsd0CE81ZJPV5PA3A00S+CwB3IKLrA/ivRPRdAD5lpaFmSPSlwe9CBdUe7S7hfYSS204JCwzddjyBp3V+K06ErvCel5BuRmPStugbCxvcnyGQI519Yxm5ENRaEjNxzevf4nKBeN2UBhW5O5uEdqgtL8vhcu3fLI2H/JxAHA5plYwslCP0RY1wjqQVUVZkmvw68g2kIC6VjfdVjuwux+VSSkllH8xiEdEDiegS9vtdRPSR/t9D3DdZP3MJEb2HiN4D4JPotO2MC/swF/0hsG8FcD8Af4+OAskDiptGSumC/C9S3oaGhoZdIyKXPI3536Pzxsg4D52R7gIAL8Ka5/UK8BwAzwEAIvohAE8gosvQGf8+l1K6WnuOiG4M4ExK6bP9GYM/gM5jIxHRHwN4CDrPjEcB+P1SOTIiNMUYjcYyzljxaqAZAgHfOFDSkMYa6Dxo2vnGQpb+b3R652mCprayZR0P0gIG3hbLpWMAZFqyrHN+AGpGZFFJdMGSDPOMgrkNcfqCUxcZeb/u+Zzdm6/Tk9/GcpfztgSV0DTRsbO50iyF03wZmnyQRsWZ8z5aeA0tw+EJ5tMpJc4Jvy2l9A8A/oGIxmigr0PnKnclOne5Rztxbwrg0p5nngF4eUrptf29nwVwGRH9ZwB/BeAFI8qygvXhpRDbVpiNPaKHf1RtwxUPHq8sG+g23Lp32rK2msrLainiz5RwDTVnuGWEeEBW5/L0aO29tTrXhEF0oPSErucn73nq5Pup7/1eW+K0hgfJK9cKI74MOirMavuk5PgtztzKY7m0vS9kG5qzuhgDTzDfgP9IKT2B/bxxbUYppQTgkmLELu77AHyPce8j6FzrGhoaGr4m4QnmdxHR41JKv8UDiegnALx7t8XaDWpGx8hvCW1TlTGanJVnHn2l4UViVxsUcXjGR56/Oj3s//JbWnre9qr8r0Skfr1vybVK0y/doLy0mYlGc4wpl9SQtXoo0RrcqJdXBGpHIfG9nnk8raxenpFvkWdWfNtNa9ZYC2vGcjDbvGfBWlgyw7B9lMqZKjbx9gTzEwH8HhE9HN0KOwC4Ezqu+UHmU+cYolRG5BnApiym2LNW6wjuydcVAlqztnNEeWGPO+UDi6wnK17J46AGFl/rtQH5nUrthV/njh91WbSgtRVLMPP9uzO4u93qnD3GHcud6fh7W2dAanVTK0AtwW7t7jYmbX7NB8bsbVKyKXnvxE/GnsIThcPbxOgaAPcgonsBuF0f/IcppbfUZ9PQ0NDQEEXEj/ktAL6mhHFkGutdc1hGIH4/o2RR34UDfgnZ08A7VivnOXgmYAgrUUGl5eZLnq6jrY1FaTYk3zPqi86vS5sFleDRFyVNOYOHcc8ez5C4mk0ssfpAUa1wuSwbu+U3tAzFJSqktHRcesJYNJOM6yGXm++rwjHFXjihvTK+1jB2NU4NaqiL0u+p4bkjlXhdLTxSl94zC2DVwHmHnbHOKjuyJzCjHisWLFpnFvgulqsi/7tKLziNLnHqfNWdFy9C5eSzBHmeCyVe1JXPKoNGi3CvowxtX2TtnpWP5wkjBXNkD+bavskHuSrbUl02DQ0NDQ27xonUmIGy4acG1kgY9UGtDRsL7Z1XYf1vbRtF6xkZrsGiC+R7cU0qYlCR0LQt7b5Xrsi9Wv/aWmpM5mG1IYsuKHoGZKMejHYQpCumpNNkepbXkeU3np/XDKulOJbWHaFSrANx8+9tcGIFc8YYWiPyTMSVyQor3Yt+dCnkrGkkD5uBrQ4bKYy3gcY9AroF3PLksNKVz8h7NagV0JHfY/Pa9tAFjUvn7XclrEfy7xas7+lxzxylbw0MTyCX97y0tuGJre8xlbvc1xy0xiM5TXmtCTaZlqXhecaVaPki972Ga72bTFMa84DNvaa9cuTb/Lh2K/0oLK3ME9BWXO33GMNPxliNaJtBbde2hxKk/WFXthFrcLbi5r+l7xmte6s/Sa04gtHtZNxjDQ0NDQ27wonSmOUetoA/3fe0aQ6PotiGyvDuR6ZafEGBR2VYfK+Vz7ZUxhgtQu6bEfUMGcP3Rrw6olzscaC0z4TmeVDySuEr/0o0kfZ8DUpl9mZNmsbMuWNvf3Fy2rm2t7SMKykfazHOmPZwogQz31Q76mcZMUpFqAyJyMeaooNrO99FXc9msI2BFjhHzdOO1k2Uu43w/JrwtdyjInzloHwOJynTinDA235razlzLZZLuy7S0qa3tLg5PQ1aG/TiWuElysKjskprEDKkga+W0x9F49U/0tDQ0NCwS5wojTlD7iNg7adsURn5d8YY7TfiEmZpEdH9OAB9n4SIIXMVL+eJ+s2RSrMSvlqthuYZa+wD7Kl81BCo3St5KMiTtjO0/SBKMzPLKJbrcswp5zwvbYZQQ+Ftg1rNUn4zzUXOQq5nPovWFul4lETEGMi/GWlnOBk4UYJZEzt0QoEAABkmSURBVEDyTLUDo1NweN4C0ambLIeVZ0SAl3wtgfIApJXB6ijFxsiuLaHEyya5Oit9a7puldPbCD063Y3AGrS9epJ7alsccYRLl/SDNyBGBlc+UPLyjB2op/TWsH5zwewNTtYgI4Wy1R6jG/9H2oD7/LjHGhoaGhp2hROlMXPwkdPTKktTzBpKI6oda3G45uI53Y9FabrKjXqapmS9j6fNR7RkDXLaGTkiqKQVb+P3WgPNw0HuRczzt2iSbcvmzYZmrH3NjD6Qyw2UV13WliEaX2rMJU8MiagBWqM1orCMzSWcKMGsNWre2KxduDRYjTUK75lSOWUDiS4f1SgbmbYsR5TO2IZzl4I50hE8j4oIzTFVuAZJA5TqkA+0mlvgVMI4CkmNaLTGDMP2xLHtSsScVyk8IvRKdV/yGpI0hrdZ1GqAwuY3BOo45kZlNDQ0NOwZTpTGnBGZHmYNs+QHbGk0U2k3pXympjRknhq0I7Rqpnkli7fUUiRqNOTod5iKGrDoL+8e11A5rTG11uTNciS4Ydb0EmHXllfDGJQoC46SL3rJA6nUBj2azTLoH8z0vcZnzStjHLxppyaMIohMyWobsXSZ4tNgaxo51o2Kg3PxsjwRdy/vd0kYA7ZAntLDwsIYyirqpSPtBxlq51bqVvPk0Op2TF3wVYAyDc1OMxa1nH/EriDh1ZEVr4SSl8yAEm2CWUdtw7RWBnmd1HLh4fAaVXIEmwZL2/Iw5WGtGg8aafARQexB1q0mpOXvWu7bM8by+9GObLnCWVqd2rnzs0rakbKM0WStb1sqv4aSAV3+9jRkoE4g17SFMQa/knypkT+NY25oaGjYM5wojbmEsdNea9rFXXiA2Og+2IuAXdda9zVE94PYBp6GYWnJ3jQ4epZgrSdGiWaxwnhaU0yFtbRkObXvZnkNTfEpJVWmhQ/yHFkXtbRTDR1nafZj4PUpi9rRtPS28q+AqCCyhKrX8fmUmh/EKZ/z3NPyR+VuN/l+6XnpHztGGO9KUFvC2OvwFq8t65qHa9dePhEhPQWiQqr0nT1ovHRtWTj4xkWeYhCt64wpp/0lTP09S44DNdSOmce4ojU0NDQ07AonSmMeM1UqabxSc8vXmrYcmW5LWBpQxEgV0ZYto5mGWs+OkrYcMcBE6ynyHrVT7igtU8q3FhHts5SfRolF6DAPY6iMXWNbd9ESvRjZHhjwqahV29zHBSZEdBsiegcRXUtE/1bcexgRXUFEP83C7kRE7yeiK4noN4g6hoaIbkhEbySiD/d/b9CHX0xET3HLMPP/Ad2H4v8yNG+LSGfkcfO/w0M9XPuXy0FOHFnGUpql8mrXUyAtNzcu0v5NBfltI/DK5f2bz9f/pn6vaFvR2s0BqwOvHXhtI3+3RaEOaup3TLyoS14ub+kf0L2T968GVr2OwXFSGZ8G8JMAnqHcexiAuwC4OxF9XR/2XACPB3BR/+9+ffiTALw5pXQRgDf3vxsaGhq+ZnBsVEZK6RoA1xDRDym3s5KfABAR3RTA9VJK70AX8BIADwLwegAPBHDPPv6lAN4K4GcBfAXAF70yRLwiIhZ9TXvO13IjFY3m8PLxUDLsSOOhhdKoXipnxBCaoWkdUR/SqHYv6zaqIUdojcg+EDI/ubmPhqlnI5F0FxhuL5pRomxKZY22u1K6EUPamIUsFhW07WIs/s7Wdf4NnJteGa8GcDmA304pfYGIvgPAVez+VQBu3l9/U0rpagBIKV1NRDfpr19WyiRKPZSesYQs55XllDEilPMHte5bUyuvIZfgcbNS2NUI5BK2pSymFGx82iyFcmiZscEvzpfDeuPfZhvPi7HPAOvBXe5pbXmmDAQOdDe9Uru1EOFlvXhWWrV5y7SlwJ5i1Wwt9kIwp5QuRaf9ZmhjSzqm4jQ0NDScVexUMBPRJQAe1/98QErpk8FHrwJwIft9IYD87KeI6Ka9tnxTANc4+X+J/64xfGlh1rW2D2zJEKDlVZoSatsK5ue0ND0ruleuEi0wRluL7MscTXtM/t4+IhLWkVdRo9WgPmFrmdYzU6BIMcE+2orDohJK2+JGy2HlU9uma+NkRPKcYhtTDi6XUkoXaHF2KphTSs8B8JwRz11NRF8gorsDeBeARwL4zf72awA8CsDT+7+/76SzemkiSrWCuSSIAd8ljj83RcfjaWh8sxTsXp41Aw4w4hTprzFoXgJhLrz/K08cL/G8YygpCZcaY+WLugYO2qAItzhiC1GOeVtE3fxqaaaxZbOEMcexURlE9M3oeOTrAVj2rnG3TSl93njkXwF4MYDroDP6vb4PfzqAlxPRYwH8LYCH7rLcDQ0NDceN4/TK+N8Y0hOl+JcD+C4l/B8A3HubstR4XtQsy/YoAQ9RrSLibTGF8S/yzrUzAI9K2JWHAkfNdNQ6JJb/5QhrZBgez1XSuGu9U0plkPE0Y7JHpcl3kSeYWHRIxGDntfNtjIrReBGN33o/GWblnyqsZHth/DsulATLbDZOKEUaTmTqqpWnhm+ewl1uW6Ecyb+EyMB5NhDpfNb0WKJmf++IG5mVn+UJAthtKFIWsGf4gRJaOSK2Dk8YHwc9FvGYqRmUo/ctnFjBrIWNXUZdqxlHXX6inZ+feBHVmGWcXWjJ26ImP163fOOjqQ03VufkndrTuLRNmaLfbIyR0NOYaw1p5nv1f6PnZFr5lzTVMahRUqLacSSfrZWT7R5vaGhoaJgaJ0pj5tCm7/LY84jGPGbKNXY01UZxuQdz7RQS2E9tOYN/GwuWJjm1tqzlW7ofoR+mKEfEA8fTUOcjtGdvNrAQbVTL3yvfmHY8pl+FZgTBfErhy8Yx69Aar1ytx+NywezxshE3n6l5Mq0xy6m711m9I3p2OXXkeWzrImblH6WKLKzqximfx93uElGDoAZPwMgDYGsGHVmmJNItHairDRiRtjF1vVt5jhXE26BRGQ0NDQ17hhOlMXNY69+5huxpzMS0A0t75u5WErWajzTiaM/Lw1ijm7RMeThrDSyD2bbG1KkWKmg0RI2WX2pD26LGW8MqX4a3D0utkZCnbe07IfOzNOgpDILbfHvv2dp0m7tcBaxl1LJTAbowjwjCjKi7TY1Xhcc3n43NV8bAet8xZ7xt695kQeNUx3w7jTOfQmBvQ3Hw8njHmZXy1BQHWSZrdz5LUFsCf9f00ZQCeQzOka7b0NDQcHJwIjXmMZqk9HHOiJ5yEF2cwO9bBkdPiyidlL0P8KbH8p1raZaoZd+jmTxY9IsWZ5co+U6XylBjSBvj01tDI2lF1XyiD5R2U3rPKbXbbftVTVlOlGCu2ZQn0xglIW4JwtzxIwKZ34/SGBbfnIUzsB8CmpfF4uUlIjvaeYKJx7F+y+tIp7Hq3SqnpMK8uGPgccxjbBgZ0d0LI2WL5KfWo/KMFNaaoC5hTJ+odbucgkJsVEZDQ0PDnuFEacwadmk15/A05218UjXNMad/tmkNqyxcc+a/JbQFP9o1UDfN9p6Rp5Z4GKONRuLWaoGaUXKskSw/FzmdfRuMopHE7+ie0BYi37mmnJYnyjrDeFonUjBP6R5WI/i2cZ3z3Ia27SxyQ52pFjFIAcHzkK59GbvyVrB+57AxZxNamMLLRIPHy1ttosbDh+dheWlEylMbx3vO5c6D1F8kn7H3M2rsJqH06h9paGhoaNglTpTGPKU2xjUsa/9e+Vt7doz27hmepphq1lr0JUramtSevbRrv1NpVqL91r5lKb0xSMYsoZRfjdE44l3hQaZzMNMPcK31vKiNx+Nr+ZQ8mMYg2nY4Sp5S/F7bK6MCpQ9qbcwi70cwhu8tuWRpnVKjD8ZyzR6vq4VFXbXGCBBP0EfCLWEMxE7DrhU+JeERESoyz+h39NpHNH7+zWmNVXxRRvlNI9demCeMrXJOZSOqoRxrFKS28s9AqcN7H1ZqO1pH9jRmC9qG6bs2FtWkXbrWfkfT3paTiwhQT7BZwjgqPCSva5VlprSbBYar2zx3v5JCwJWGkgtZdNAsxZNCWjNC5+sM61p7P2t3ujHCOBKnZjAutdtJjPv1jzQ0NDQ07BInSmPOKE3PLZQ0LItTBqZxW4tOSbfRoC1teIz2XPIeGcNRlhDVyqLacsSrwwuX4C5V2vl/Eh7/rYVLFzJrUUzJA8h6ZsA9M5ova80laLOAZPSZnCffKzqySEw+781sttWUa2iNGtrzRAnmiHAp8VgSnjBeLrcTyF4ZpzT+WVTKbKaXgYd5xiyerrULn0zbK3/03hjD3rZUBoec0lvfitMaHiL89ypdbNpDShSDZzsI0RpO3cj+keOmZbnO1XKyQWKBIYUTKb82yHh5emWz7lntvIZjblRGQ0NDw56haczG/YySpmSNrpq2XFrpVoNaWsNaBShPb5FaMk8rh2sbOhXrbrZp0KldeRV55ygi2pp2L2Kxry0frzrNwBzR1gcaIgu3tPGoEc2KtxTfsma70Hw/MmNxy8l+D4yRIyibKWG2keaVoSMqjC2UGpJsbBxcCEX8WWswRmBF+TkusL2zEb1y5bJtlDN3Kt4RUXimAhF+ViLKN2/LNQP6dDtKt5TytbhsHr/WW8ASeNF3lsvdJQVYGoA8IT3DpnCWaZQoGx42llLz0KiMhoaGhnMYJ15jjsLTlrURnxs5POziWCdLy/RW22lUhkVfyDgyrZJh0tNIFljXWfQZjsgUegy8b7xrRLVmwPdCkAZBrW636SMWrNmlF1ZKS9V+Ny6UOE4f2AWtMRYnSjBnbONwLu/LeGOmz1PC6pglDwLrnkZfaILZ82yJYMCF92FyCTBP/2x0oqi3B0dkB7ToN/MEtOeFkH+SDB8hhGuppaVRTxrHrNVpifLT3kEOQFFhbHmvRN45MkhX0VzxqA0NDQ0Nx4ETpTHXaAiaVmZpyfn32dCWp9AeLfoCsOkLT2PeBsvlcG+GqetxG41PIuJfPIN9MrRWLi2tGqOjR2VEDIFjNWmrvMDwvTVDedaWtW8tF89IOs56Vw3arK7Wx1mipCkP3ql5ZfjY9mRrL/xs0BdjEOl8chppCe8SHSKRO79ncffSOptcYO2+zRaVYKUh66uGXx6b9i54Za8cEpa3RobmYSLv7ROmkAF7+FoNDQ0NJxsnSmO+/m3v6N7f8CFNht8lC+dO40tlOpbvz2gYTspUakb9vgO0ec/TtuRCBF6mQflZOC+PzHNl8DtYhx8crq/VMlrlMzShjbpN6/KW6jnH0d5Tfg95X+Yn72n58HuDNFmZB2kLTTTXJxGbdVD8+1rvYIF/H1LaFy+bDM+/S+D1tFgM2/5isb6/mAPzXs1dLhjtt+j+Ab3GnMrvR7OuDoGO0uDveXDA6pZdHx6wd6ZYP5Lv5/UvrT1Y+PpvuaUfgYFSjdfzOQwiOhkv2tDQsNdIKRWHvhMjmBt8ENGXUkoXnO1y7DNaHflo9TMdGsfc0NDQsGdogrmhoaFhz9AEc0PGq892Ac4BtDry0epnIjTBfA6DiM4noncT0XuJ6ANE9NQ+/FZE9C4i+jARvYyITvfh5/W/r+zv35Il97+J6HIi+ics/Uf1aXyYiB7Fwu9ERO/v0/kNos5WTkQvJqJ7Hse7R0FEtyCiPyaiD/Z19FN9+HcT0Tv69/gDIrpeH36aiF7Uh79XvM8L+jr6FZZ+VV0T0T2J6MXH9PpFjKifU0R0aR/+QSL6OZbcHxLRFUT00yx9q63ckIje2NfbG4noBn34xUT0lOOrgf1EE8znNq4FcK+U0ncDuAOA+xHR3QH8MoBnpZQuAvAZAI/t4z8WwGdSSt8G4Fl9PBDRbfr73w/gkj7shgCeDOBuAO4K4Mm58wB4LoDHA7io/3e/Xb7klpgD+Dcppe8EcHcAlxDRbQE8H8CTUkr/B4D/CuDf9fEfBwB9+H0APJNo5WT1rwD8YwAHrM6q6noPUVs/DwVwXh9+JwA/wQb4hwG4C4C7E9HX9WFWW3kSgDf39fbm/ndDjyaYz2GkDl/sf57q/yUA9wLwyj78UgAP6q8f2P9Gf//evQZzgM7jOAHIrjw/COCNKaVPp5Q+A+CN6AT/TQFcL6X0jtS59LyEpf85AEfTv+l4pJSuTild0V9/AcAHAdwcwHcA+NM+2hsBPLi/vi06QYGU0jUAPgvgzv29Gbo6WgKgvu5q6/oIXT3tBUbUTwJwAREdArgOuvf5fH+PWBwqtBVeP7zevgIgt+kTiyaYz3EQ0QERvQfANeg60P8E8NmU0ryPchW6job+7ycAoL//OQA3Sil9AMB1AbwNnYYziCvSuXl/LcORUvqplNLbJ33BCdFrdt8D4F0A/juAH+lvPRTALfrr9wJ4IBEdEtGt0GmF+d7zAbwdwCyl9EEAN0J9Xb89pfRTO3nBLRGsn1cC+BKAqwH8LYBnpJQ+3d97NYDLAVzeC3mzrQD4ppTS1UA3OAC4SX/9spTSM6Z+t3MNJ2rl39ciUkoLAHcgouujm3J+pxat/6s5tqc+nX8twq24Zhr7jH5q/SoAP51S+jwRPQbAbxDRzwN4Ddaa/gvR1eHlAD6OThDPASCl9AYAb+DJKlkV63ofUVE/d0W3ZcXNANwAwJ8R0ZtSSh9JKV2KtRYMnGN1sE9ogvlrBCmlzxLRW9HxhNcnosNeU7sQwCf7aFeh03yu6qei3wDg01p6fdx7st8XAnhrH36hCP8k9hhEdAqd0PmdlNKrASCl9CEA9+3vfzuAH+rD5wCeyJ59O4APG0n/Paap67OKmvoB8HAAf5RSOgPgGiL6c3RUz0eUpL228ikiumlK6eqe8rhm4tc6p9GojHMYRHTjXlMGEV0HwA+g4wj/GMBD+miPAvD7/fVr+t/o778l2Us/3wDgvkR0g97od18Ab+innV8gorv3nOkjWfp7h76MLwDwwZTSr7Hwm/R/ZwD+E4Dn9b+vS0QX9Nf3ATBPKf21lnZfd1PU9VlDbf2goy/uRR0uQKcIfEhLu9BWeP3wemsAgJRS+3eO/gNwewB/BeB96DjBn+/Dbw3g3QCuBPAKdFZ0ADi//31lf//WhfQf08e9EsCjWfid+/z+J4Bno1/av4//AHwfuunz+wC8p//3AAA/BeBv+n9Pz+8A4JYA/ge6Ae5NAL61kP4kdX0O1c/X9e/1AQB/DeDfFdJX2wo6fv7N6GYjbwZww7NdF/v0r+2V0dDQ0LBnaFRGQ0NDw56hCeaGhoaGPUMTzA0NDQ17hiaYGxoaGvYMTTA3NDQ07BmaYG5oaGjYMzTB3HDsIKJvIqLfJaKPENFf9ttL/rPCM7ckov8+Mr+Liehm7Pfz+x3UIs/ek4heOybfKPrVhfkdHz7i+YuJ6NnTl6zhbKEJ5oZjRb8C7PcA/GlK6dYppTuh2y7yQv/JrXAxur0dAAAppf8rGav5zgZSSvfoL2+JbslzwwlHE8wNx417AThKKeUlvkgpfTyl9JvASmv8s37D9SuI6B4yAS8OEf17Wm9y/3Qiegi61We/Q0TvIaLrENFbiejOffz79Wm8l4jeHH0JIro3Ef1Vn9cLiei8PvxjRPTUPs33532b++Xzb+zD/wsRfZyIvrG/l7e5fDqAf9yX84lSEyai11K/cT8RPZqI/oaI/gTA97I4NyaiVxHRX/T/Vvcazh00wdxw3LgdgCuc+9cAuE9K6Y4AfgzAb0TjENH90e3re7fUHR7wKymlV6LbKe7HU0p3SCl9JSdCRDcG8FsAHtzHf2jkBYjofAAvBvBjqdsw/hDdJvoZf9+X7bkA/m0f9mR0+2XcEd0ugN+iJP0kAH/Wl/NZTv43BfBUdAL5Puj2kM74dXQb998F3R7Kz4+8U8N+oe0u13BWQUTPQbdfw1EvTE4BeDYR3QHd9pLfrjxmxfkBAC9KKX0ZANJ6n2ALd0dHqXw0GD/jOwB8NKX0N/3vS9Gd/PJ/97/z2Xd/CeBH++vvA/DP+nz+iIg+E8xLw90AvDWl9HcAQEQvw7AObtsxRgCA6xHR16duf+SGcwRNMDccNz6A9WkYSCld0k/pL++DngjgUwC+G92M7qtKGlYcQt1+v7Xx+XMeru3/LrDuY6VnNMwxnNWez66tcs8A/CM+M2g499CojIbjxlsAnE9EfOp/XXb9DQCuTiktATwC3bFXElac/wbgMUR0XWB1biEAfAHA1yvpvAPAP6HupBIev4QPAbglEX1b//sRAP6k8MzbAPyffT73RbfJvIQs58fQHYIwI6JboNukHuhOGLknEd2Iur2UOQXz3wA8If/oZxUN5xiaYG44VqRuO8MHoROIHyWid6OjAn62j/L/AHgUEb0T3fT8S0oyapyU0h+h2+f3cuqO28r87osBPC8b/1hZ/g7dQaGvJqL3AniZUex7E9FV+R+645ceDeAVRPR+dGcAPs94NuOp6Pa3vgLA/dEdzSTphfcBmPeGyCcC+HMAHwXwfgDPQM/Np26f46egG1jehCFn/5MA7kxE7yOivwbwLwvlathDtG0/GxqOAb3XxiKlNCeifwTguSmlps02qGgcc0PD8eBbALycuhNBjgA87iyXp2GP0TTmhoaGhj1D45gbGhoa9gxNMDc0NDTsGZpgbmhoaNgzNMHc0NDQsGdogrmhoaFhz9AEc0NDQ8Oe4f8HKCipN21aI8IAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sh.thaw(g0.xpos, g0.ypos)\n", "sh.fit()\n", "sh.freeze(g0)\n", "\n", "resid.data = sh.get_data_image().y - sh.get_model_image().y\n", "resid_smooth = resid.smooth(width=6)\n", "resid_smooth.plot(vmin=-0.5, vmax=1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Iteratively find and fit additional sources\n", "Instantiate additional Gaussian components, and use them to iteratively fit sources, repeating the steps performed above for component g0. (The residuals map is shown after each additional source included in the model.) This takes some time..." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# initialize components with fixed, zero amplitude\n", "for i in range(1, 6):\n", " model = sh.create_model_component(\"gauss2d\", \"g\" + str(i))\n", " model.ampl = 0\n", " sh.freeze(model)\n", "\n", "gs = [g0, g1, g2, g3, g4, g5]\n", "sh.set_full_model(bkg + psf(g0 + g1 + g2 + g3 + g4 + g5) * expo)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47352.1\n", "Final fit statistic = 47306.6 at function evaluation 218\n", "Data points = 30000\n", "Degrees of freedom = 29998\n", "Change in statistic = 45.4492\n", " g1.fwhm 6.73141 \n", " g1.ampl 0.358694 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47306.6\n", "Final fit statistic = 47295.3 at function evaluation 362\n", "Data points = 30000\n", "Degrees of freedom = 29996\n", "Change in statistic = 11.3247\n", " g1.fwhm 5.88307 \n", " g1.xpos 41.7803 \n", " g1.ypos 81.4247 \n", " g1.ampl 0.45753 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47211.4\n", "Final fit statistic = 47187.6 at function evaluation 215\n", "Data points = 30000\n", "Degrees of freedom = 29998\n", "Change in statistic = 23.794\n", " g2.fwhm 6.77442 \n", " g2.ampl 0.346582 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47187.6\n", "Final fit statistic = 47184.1 at function evaluation 341\n", "Data points = 30000\n", "Degrees of freedom = 29996\n", "Change in statistic = 3.46844\n", " g2.fwhm 6.19417 \n", " g2.xpos 20.8277 \n", " g2.ypos 81.6412 \n", " g2.ampl 0.399908 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47173.4\n", "Final fit statistic = 47116.6 at function evaluation 227\n", "Data points = 30000\n", "Degrees of freedom = 29998\n", "Change in statistic = 56.7946\n", " g3.fwhm 6.43218 \n", " g3.ampl 0.235481 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47116.6\n", "Final fit statistic = 47115.9 at function evaluation 303\n", "Data points = 30000\n", "Degrees of freedom = 29996\n", "Change in statistic = 0.728928\n", " g3.fwhm 6.38175 \n", " g3.xpos 177.458 \n", " g3.ypos 80.2233 \n", " g3.ampl 0.239637 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47121\n", "Final fit statistic = 47090.2 at function evaluation 246\n", "Data points = 30000\n", "Degrees of freedom = 29998\n", "Change in statistic = 30.818\n", " g4.fwhm 7.24022 \n", " g4.ampl 0.112126 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47090.2\n", "Final fit statistic = 47087.3 at function evaluation 330\n", "Data points = 30000\n", "Degrees of freedom = 29996\n", "Change in statistic = 2.93503\n", " g4.fwhm 5.94179 \n", " g4.xpos 135.942 \n", " g4.ypos 59.5875 \n", " g4.ampl 0.156213 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47082.3\n", "Final fit statistic = 47068.5 at function evaluation 249\n", "Data points = 30000\n", "Degrees of freedom = 29998\n", "Change in statistic = 13.7724\n", " g5.fwhm 3.07232 \n", " g5.ampl 0.479624 \n", "Dataset = 1\n", "Method = neldermead\n", "Statistic = cash\n", "Initial fit statistic = 47068.5\n", "Final fit statistic = 47068.5 at function evaluation 320\n", "Data points = 30000\n", "Degrees of freedom = 29996\n", "Change in statistic = 0.0196613\n", " g5.fwhm 3.05445 \n", " g5.xpos 76.0488 \n", " g5.ypos 73.1002 \n", " g5.ampl 0.483818 \n", "CPU times: user 25.7 s, sys: 267 ms, total: 26 s\n", "Wall time: 26.3 s\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAEJCAYAAACnhI2ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztfX+wfElV36fnvV2URVyMgK4VV61lCxJIgFpN1FhJpQxZ4pofwiL+CBpTWxFjafkjKiFZBDWiSUyIVUQlZRmDMUCJJBJF1DKpjRgTA/gziiSy/lpK4rLgLiz7fTOdP6b7zZkz55w+3ffOe/OY86mamnv79q/bt++5pz/ndHfKOSMQCAQCh4PFZVcgEAgEAtsIwRwIBAIHhhDMgUAgcGAIwRwIBAIHhhDMgUAgcGAIwRwIBAIHhhDMAQBASml52XU4dEQb2Yj2mQ8hmAMV0RfaiDayEe0zE6IhA4FA4MCQjmXmX0rpOG40EAgcNHLOqRXn9CIqcij497eu/xdsnEDP+TUNq5UcRsOXKyCvdtOsVnoa7ZhDqrP0X49PT7fDT0+346Ry7YTly+8HAJZnm/CzM+BDj6zDAOCRR9Zh9Ro91u6H1pPWLS029aFxJEhtuzTa70Qpk5bTKl96VrTMrDxvC/y61ld4OVYZUl+x7o0/DymfmkcF7Te8r9Xjk9Pd/kVB7632J3rPtA/V+5feL+2Y5rVa+d9PqY6tMivoc/ozL/hS/I3v/4HdSAKCyggEAoEDw1FpzBWrla591a+eV3OmaS4KPXXT4LlPr/ZJoeVHw62yR+6Nt79U17za1vCsPFp9wxsuoeZtjR48+Z0s9GdC8+htz9bIxAPr/aLleK7v8906WQDVjSQPlkPrKfVxa4Rg4SgFM7DdcaWX8qKFrQTPg5RomZGPSk2yxHZnomgN+6z6WYLPW2cubLxUgnSu+XVtCTUh3ojQarWBRj3Ul93TP3m49tGjlJVEV0j1m0MRmAPej1aNC9j9zvNxlNpfSqc9J4rUZJZJen/UQCAQCFwEjkpjXgrDC20IMmIEvEgt26PRUA2pBxoVUOHVktMCgKM9+b206iwZYDRNWdMs80rWmnfil/86kpjyjC1tSyq/xteu1368xOaeLc3ZalfJWMfT92CUSvFConJ6RyYc3tFya9SipuvQmI9KMFcsWee+aCrDGp56hmvaizPlJeCUBrB+kXsEsgSpTtJLK8XTKBWej1co87Jr3CXaz12iNVrg7VfL1erEr1tUBk9XhS4X0PQaIHtiSF4Yrb7Ehfxy1X5ec8DzfnpoBQ76cesppwdBZQQCgcAVxlFqzBT8S9/7pW2BajLaUIsbGnqNfpoPqhdamZb26dEkqNUbQhnSPXioF6+m3DPycI0GVromqoHG0Sz/I1SG9owlrV7yDNA0Zp73CBWmIa/QrQq2RpB8pEDTVXhpldp2C/hGKVaeokH1xC6f4qgEM314tMNxaqOiV0jzTkR5MD6sldxsWlSGx3rOXyoPHSANn+k5P/YOIxfk/lt11yYr8LpK4ELZI2S3ni36aYrFQha0liCzrkkUBO8PnucktZ0kdCX6gh7PJZB5PZed71SPQqBBete8+bVokVa/BjZt2aMwBZURCAQCB4aj0pgpMtFkKSTPjR5wuuLcct6gL+hx07pr0BiArO30DLnmNHxSz4HzMlkcTVu22oFb5Hsnf+xoQjVfpTzrudBwayIL71P0HiTKQ9LuRoboLY2Zx/Fqy9o7NNXTaQSeOlOKwnrXltj0h9b72DPqOwkqw4/WjLAWLH6Qhm3xrZAfvJdf5uVanUN7KXvR40Kk8p/Ki8w7sbd+EpXgdV9UhQf8tIbnnq0PvHaNlp8dQq7HLtFDKVmgbV/fIarUaPWc4rlBn6d3JipFWvgFLuWbW3lWnAjv19ZxeGUEAoHA1cXRa8yAPCSz/DK9tIA6VCJaOh9eAbKmZxkdPFSANXQdMWz2xusZbl80dkYtJXxrhGPc/8gox6oLfYbS5BFeZ55eK7NlKPZA8yzRRp68ntR33MKUiVuWVu7VhHn8Vlna+7UVFlSGjB4vi+awpFOQaNZpOryS4lGMCGN6PHIvnpfiZAGcdVrbNVhWc/7hklylNP6+VY4GjxeB1q9afaiFKqQtft6DlsCYKpT5deuDJr0D3A4zuqCQBI1KoTzyHHlzLxepzXtonKMSzBS9Bom5BLQE74snvUBcKHvqKXGNUt7c99rS0OgIQEIPf06hjSL4lFwt3pTnI2mBvaMB65n0whppefrzaNmaoPQoOqrPPjb3wgUydXns0ZZbhmMe5hHOrY+7qBmT8K3+ExxzIBAIXF0crcbsgcXrzq01V1icVs/kAe1aPW55RvC6UA1NK7+mmfNrb7kVVtDJOx7NzYo3MsTnGlKLd7TK5/BM0OAcuVQ/C5ZnEteWLbfEWg5dIIpTCVv1squ1k86rPVvtzEcWvZyzlr/G2cd6zHtCa6jCr1louTp53YlaL7uHU7Rc1KRZZF6B2zuzy4JFn3Bh5Zn6zqHVz/JH1V5EC60PZQtSv+GL7tB8e2esQslLq0fruiWMPbCojBE3OQqtPaz+7fmIU2VHE8axiFEgEAhcYRyVxjzFOZ+m82o+XqNML0Y0ZmviiVU/aYKMdSsni02EJXE9s4woWpnWc7K8XCq8Te6Zqdd6ftrOIJaBqBfuWYyOtu7p93POBPXkO2L8s0aAJp3CRiHAut9Y1BBN2+MdFesxNzDqIUDh5SLrJY0HHBHePdx3z8wkK4zCGvadnrIZWuRaLX/UHaq3nbw+qBpa/H2L5pIEc4u37umW2n513AfYw7uP0A/W8/D6K9f0rVmawHQao5ZV8+TlepYPAPpdVYfqOS15IBAIBObGUWrMwHwUQ4/lXpss4dFWPNSDpMFVTF03Y4c6qPk00olDQkJ1SPBq1lOeYYuC0sIkjVnzxOAas2cZVs2oNwW9/sZezEnH8Gtz1ccy3GnvHl+ul8JLmUhxeox/RyWYtRdZcnrfJzh94bWiezqeh2Jp5TMyxNXytCz1UhqNBpnrsXgs7Fq4ZYU/Pd2Ec8Es8c8WtI+eR1h5PSdq/S4CvS5+UwX0KDzt4bXNBJURCAQCH2Y4Ko25gn+tp3zdpi4b2sK+vsieMj2TMlr+zZrBk8Iyhp4bfrBLcXgs51JdRsIpHQHI3hfSMac5tDKlUcRZ476mUh778BjqKZOXf5HasVSHKSPXVppeHKVgngP0gUqr07Us5RZGH/woX+yB1XE55TAijKW6aQJ7dHaht101Okhbl0Q7rkLZw+dL95wWAAbpi8ugLEbLadXfQu/WaRZ6FBIrfA4ElREIBAIHhqPSmK1hVE96mo6vEaDFtwwcHox8nTnNYk3c6DEacWhW7JOFPty2hvLS8c7UXJqXWXO9TKlNLSs8pyIs+qKm9ZQ5MoT3bq01FaNa4UUY0afC6wHlvTYnLlQwp5TOR6E551TCXgjglSXKp+acfzGldDOA38Zmobx/lHP+thL/QQA3lPD/lXO+rYTnmucoVqvNSy45xvdwYVMF8QgkwUsd/U+E+vdwa1KaeqwOKR350U0JtLUuzPS+YnbuzbOMKk8v8cQ8vFWuF1n4KM35Qd/XEH0OoSz1Ab7Uq3d9Y8szaKROo8jZH/eiNeb/CuA9AJ5Hwv4lgO8F8JkA3gTgYwG8FcAq53yaUnoXgG8B8G0ppVdiLZQ/CcCXAHhpSumGnPNDPZXQ/Bethvf6X06JU+sw1SgjCVq64hfdk43HG7mXya5BtD7k5evJlu+Tp9WrRxhb2q4mjFu+ytY05Pq/VARxjzDptUd4nuGowJ2yzx9XHObaM7C3Dhde5niR/cg5/2UAv8yCE4BrADI2GvLHAHhzOX46gJRSSgBeAOChnPO9OeeXleuv22+tA4FA4GJxCBzz3QC+vRx/Bgl/EwDknB9Yy2Q8B8BHALiPxFkBuLXEG6IxNKst1TBHMTLsnOq+p9WBeo/03NcUVybvUFnzxGhRFNLIgi4BasE7g8t7bJXZmthEj5erXQrD0/bWc5rizTO3liyV1eLdR/rgnK6AU+1RI7h0wZxzfjmAlzuiavJkVvaWPwRpS3ZPegkaP8YNia1OOZkro3XC7r2N5N/DV/Z0cM/Ub+1jxg1kljCWwjQ+2Upv1c2iJej2SlQg8zTSNPYe+msflEXPYlDeemhUY29/n1PRaeXVpKk6OOYLpTI6cTsApJRuLOdvAPAwgMeROAsA79QySCnl+ttbLQOBQKADHrl06RqzgvcCeFY5fjuAnHPOKaVXA/h7xWvjS8r1O7VMKL3hEc78i0jdzbim4lmC0JqZZdEJLcOk1xApucXxY07ZaCODHut+r4bk0WgorSEZH6U6W5pcy7DXUzdge5KR1Teka1nQjOs/15o1XIYnRe9uO9K5FtZyTbUmpViUUS9G87LieWjXi3aXW6EY+Iqg/MOc8xOFqJ8G4B1EmL4MAHLOX55S+tsA3lXC397rkeFF7QS9017dnYB0oqlctqcukotbPeYfIOmFGx0SeoXcFI8Dy/2vlbZVLwuWeyKFNrVa4pAl+kK7n0P2Fe6hiaz0lq9xizrwlNOC1wPJU5eDdZfLObuaKef8Tig0S875Bik8EAgEPlxwqFTGpWGqocHKT0O9lAQtQRrGecv3UBl8ZMCH4icsj9a9XHWM3NuIt4DnuJXvFJ9eYPqOINQwTuGliawwbXQ35Z0c8arwUhkeo39P3UMwQ+dxvcMWz8Pi+YneFtBn6I3AS2VQeFbL4/fipTmmcoAj3CvFyKQKrQ6e9FY/8RwDbc8SXg8vXbQilBUvw6LveukTq24eu0Pr3Wp9XCTFogfeNHNse0XxYaz7BAKBwNVEaMwGvMMWi8LoRdWa54BXI7Q0Xm242sqvlu/1KtEMYdS/t1XePmDRN55ye7VkS/PSlh3V6tdTt4qWsZtebxmMLVjas0bJjTznkbU1erTfLueAQzX+HRLmpC+81uGWG9p5nco55Xs9E1wkisDilXvAXz7r3kZepFZ7UupCojJ626YXPW2m9RWPMOYvOqcO6ESX1gQYqSyNiuldrU6yS5w4+tecXjFWu1Uk5R65kNaE8ZwKwYfLBJNAIBA4ShytxqzB+gq2KAtvWm3YxjVN+tW06A3NWNgzVZdC2pEFaHtrWGhpnJJWydeN0LwVRoa4I4bAEc8Y69xD0wD2Dtz8uAe9VJ02uqwG417ay9KqNS8VXrcR6gWw6YoslE0xSp0drB/zIaLXo0A693CEvLN6udcKrYrSTDNtSGm9YNI1y0PDO7utV4BbwpgLZq3dp/KdU+J5Xbo8PKbU9ppg1rbAWq70vjNFCbEmKUlxefjUyUa0n7TS8A+IBi2/HsrMQg8ffVSCWesUHm1Wi+813FgcHDeWabywBN7PNM2l5wPUmsnWA2nWG385tJ1fNMFcz7X2tjS3UYPZSFz+0eP50F1fNM1O2zOQCmZtw9fzcmoaEubpWxYsm4V39ujIaE7KUzrWyuHtLPVvy37R21aj9xYccyAQCBwYjkpjrtC+rPyalc6rKXvqomlYvTOdrCGldp9cS6fXWvUeRc8kCE1j1izy1lDVoy3vS4uukJ6HtoZ01ZZb9IV3J+4lAAh9ayHQYb33siKUiZfSmKKxa3WR6qbB2qOThvXaNLR3KThmBS1e1+urOkdH0uDtsFocOlyjswitMrwfAK+xqIUeAxoVzB73qBaPOMX3d05hrvG1PK4mmDl9we9L+wBoNFnN64ylG+nrfBuzCouj5pDsCPUeahm9dbPoPEshaNUJ0D86owgqIxAIBA4MR6UxV0haphQ+J0ZmznG43ewEAxPPp5cusQwgUzVmi0ap55K23Brt9IRP8czQNPSTxe6u39pz0wxpmsbM6Qvtfmh5J4ttjVjqA2mx1tY8/cPbZtpIVbtnrQ/s8/208veUP7dh86gEsySkrG2V5oBHIHsfojaE0l6iBWSh4aVLvJx6r+cG5UatOmjUhUUvWYJYGvK3oMXR3NN2M9j2ElkInLk1rNd45R7QZ13TL6C3RWuvxR5wbyRa3qgtZU60eGWNSqnY1+qLQWUEAoHAgeGoNGZtAZaKUe15XxMaPLCGhMDYztgt7weL1rBwXrdSzzqsHnH6p+cj7dmrOfM6SpvpanlRLRmk39Fwqw4tg5+kZWva+FLp35rxkPskT/Vrr2j121a/alEuvX3CO3K0YLUnAKTmhlIbHJVgpg/Rs0qW9sKMDr3m7CwWn6V1UM+CSFwQt6zSFt9s8bhbtIIhmDxeMV4ruEVleCiC3tXdaj+ROGYu/DxrIJ8s2vfZ6pvnazBje7h8erpbRzE/QaFpQep3o1SGtoZ0zcdCS/jT/i3Vx3vPvR99MY+xZIFAIBDYF45SY6ag2nPPQiyW76UVx5PPFGjao7Xjd71e09O8WlqyZBzxGHvo9ZZPq1SOF9Zw3/KmsPKqx5Jm5KUVdu7HeDatqeUUPZqnNFWb11PLc9HoAz3oGfFY79hU7xGtbho0OoufUyrKi6MSzBXaw/Vsq2ShdyjlyU/7mHih0TdSnHosDek472wNvbX1MfhaIVOpnZH2lF4S6yNs8YYewUzDLMHhpa084Z58ViuZ1uD0S43Lj6d6VUhl9NxP/bDQfjjyXvTU+2SxS3dqfSA45g5I/OeIwWBuLXdqWZrw8LjF8Q7KhbG0NGWP8e+8s8Jn7OuBx0VOOq6wNB4tr+bLZwl5uwgX5up7VJieLLClse9srup49hovu283ONreUwyT3o+mJZRbfWPRIZiDYw4EAoEDw1FqzB5OU8IUbWU0bY9G7+G9peuWltzrlWHVg2rPU/c1HBndUF5Z83Boad8SZcH5am3bourixWF6QTjq1gM6gvLcP3XxA5i7nUIHaPfg6YtSHaT6VSqG1qc+A+92UL1ucFKbSVqyOlILKsNGL58FtIfF+8Jc/CIFH3bSc89C9fXcyr/WUXvZPMLZukfrGVrGMmndYstwY4V7+WrLxYvnNdI3RzBK20m0RvXJ1njbESXIC0/b8np4YFGdLV6ZXh+l7YLKCAQCgQPDUWrMXmgkPo9z0WiV2RpGStSFdzsnLf8R1Ns46/SG6dH2qBZj7TjtLUtL0zv704O5tGfNw0K63kqrGZM962t4+83cowbrffAar6UwzfgH7K4Ho+WlluuJlFL6ppTSfy/HT04pfbq/iJ283ppSyuX3YAn7nnK+TCldX8IeJvFySumFJTyllK6R8NeU8BenlDqWou6D9BBGOw9N3/M7PV3/tOs8fw7KJ9dfFcpnZ/ZvtQKuna1/y9XmV6/TsPqTytN+ufys9PQ+NCSjbU4G2lx7ZrUs6gaoxfPA80GcA9o9evoaPa/X6++60017eNrT2x6t+699Bdj0IZ5Pb56td1C7lhbbfYJ/rHvc5ZpNlVL6QwDfDuDPlaC/A+At/iK28roBwDMAvAzAJwC4IaX0SgB3AfgCAO8H8N9L9L+Sc0455wTgfQC+u4S/DWtNfwHgPwN43khdAoFA4FDhoTIej7UQ/X0AyDl/Y0rpGwbLe13J4yUAkFJ6CMALyrUzABnFdplzvkep59MA/E7OOQO4o2jNdwH4oLcSni+39IVv0RqjsPLqGuKXf82TQvK8WJ7Z2prkx8zLWCxkSzjdQcXSUhYLYiBCe4GpEUqJanI1fu/zHDXqSF4uPX6/+zSeSfmO9nNpworksSGFe0HTeZeFneJf7W0bzQA8+sxcHHPO+Q9S0cNTSjeOFQUAuBUAbZIHAHwcgNdiLbQzgD9ZL6aUPgTg+nL6BeV/AeAdLN/bc87PAfBdVuFSo/Kh6Y5z/UJP2wNPGm3fNy/4TCgOakVfnm3CVoWSqOetSSVavhX1XvkWQ+6XgRxrC0xp0BaQt/bPc9Wp8Ww06sgDi6qh+Xq4V+u6lq+UvnW9dY1vUSXVwQspzcjEp4sSxlIevfAkv1Y0W6SUfgHAewE8MlieyLLknL+g0BaLnPNDJPxRhcp4CMAPGPlOdYsNBAKBg4FHY34MgN8D8JEAPhXAewA8cbC83wTwKeT8RgAPO9J9K9Y8N7DWuG9l139KSqQZA+egMnrykuDRji2tRfv6052QKbgmTIfU1bgHCB4a1k1IdcP20FvaGNbyDNAoktYGq4loelTr4+GjGnPFXNPJAZ9Wx9tjCq3h1XY5pHuWRmb0GdKdUrTn7aEcOCz6YqrnkKYZz60lU7lUFM8dNAVzzvkRAE/wFdnEFwJ4b0rpbgD/BsANAL5XiphS+o6c8zeW02/EmuYAgF8D8LS05lZ+rNTxVUrdz286pZQ1WsKaDTYXlQH0C2PPkEriKqVrKy5wy3EVypTKoNl5+jet5oqc06G3tHOy54X1CBPeLjScLywvpeHPRaOD5uARNfdDjTu11genQnpf/LMFjaZYLOQF+XvosAprGzNNEHsFtCWI98kfa8J4q27aBeaqtvMbrNADAH4JwEuxNiZ+IOf85Ur0byBl3Qjgq0r4n8X6HV8B+BwArx+pSyAQCBwqLI35ReX/ZeX/35X/L8YETjfn/HRnPI2PzhicGKNpS9p2PR4qw4uRrYlGtXRp2U0KrklTKmOFbS259aBPWPwFOee0xnmeQphXS9biSO2ZFvIzr9BoCWnFsjkpDAtcQ9S0Z46pnhtW+pE2qHVt9UVadoVns1+LumjRQTRccwiokPzULwKqgMs5vxwAUkrfzoTk393nRI59gr6wFd5t4EfRoi+sjjEimK0hoeSqdT6Ro6bHtjCWHjTtDEushfN5+SBUBjkGdl/wXmHM47YEMKem6PP1Cpmp63N7uE7K/wO7gojWYbTNvKjCShLUlsvkHLA2cODlaZRQq27eD5hGb3rafI52cWmeKaWfyjn/lXL8punFXg60aZKaUW8OXtkDi9/SjqUO1jImSYI5r7a15CW2hTHXhsGuJ8jCmcavZXKteco0ZtpOrY9rDe+FtC/knJCel6QdSgZUns+c/ZPm5/FDluwcrUWFerRkS/hONfhR0I92z64x+4BHML8OwJ1MS37tnuoTCAQCRw+PV8bzACCl9FkATnPOP7v3Wu0J1wl3q2moFjc5BzxancdDgw9tvUOtJdWesaEvMnRPDE17rlozsK050zSc0hjRQDWtmNsGJNe5Hmha60XxzBUWP08xlVf2glMzLbpgtPxezwtv3tI71HrP5/R4yR0EcFMwp5Q+wM5LIfnRnfW6dLSE4WW4HLUwUiev3+nW9fJPo2lJFtgWuOfz6LFNa3COeYqBSnJlnCIouaD1bEt0UcJZM14t0C6/RWt4uVcvN06PNSrNK6inuMLNzXdrsKhDj1+6Fx4q44wcnwB4NHzurYFAIBAYgIfKeCw9TyndDOBd+6rQPtGrMR+aBi1peC3jHwW9fj7ZBLsuca2vruR5UTVujdbg2rMXSRhq1uGl5uYo7VIiwbNrOC2TghovNS3V0qj4NTr5xwI1BEqodENvfawwek3TqLmnT8/GqJ7dR3phURNTMJWm8WLEH/gp04q8PPTyuoeGqT61vffGX37N84IKXUpr8HheSEKWh3MXuanuhN50tZi02AgUa4Emj5Digsk7W02rp8dNb5QGsKiQ1ip6PbSGB9ytr5V3i7aUvLa8mJtW8XDMEmX9wPSiA4FAICDBozG/iBw/CODVZWr1lcPI8KbXQGSh9SWWhqBTHf21CRXn56u1JnyGXZxgW2u2KApJM+b+zR7Q+0xMw9E0aQprNOEd4lt1kwxxmk+2NESfqlnR8iU6o8fwN5eWx5fg5FqyVhfvesoUnEqh555Rk+XNMwJPuy0bbSHBI5j/Sc55q9oppRUPuwrwNL71YvcK6VHKQZrdNSJUPB+iKlQp/UD54ypYe+bgU77Zi9akGomuaL1U0uSIqS5eW/mRY6l9OG+sCUI+9OdxeoUPj+edeATYfK9nre96rN2Pd1JNC/xZagqNdE0Lq+AfWp6Xpw9N5c09omPr/SoL5fe8c4FAIBDogKoxl41SbyjHnGceXSj/UjGnD2pLe+4pixswJG3HY9ywcLIAlgIVkBfAgmpv2DXsAf2UhAeSxkuPJd9lYNvw58U+rOlbGq1AOXGtkWuSmlal+jEP9CkrXNLwemgyrZ6WlmytCdLbv3u1ZB4+Ql94PFvm8DKxFjF6DLDe3inn/KjpRV0+Rqytnkb2COFRS682w4/GsSBRAYvFerdfYH1/p9isAnaK7VmAYr1o/sL11nBK21m6HlsuchTeYX2Nuy9Xp5YbGT2u617vbEjgoDK0MK0uXnDPEYt7tmgRDz3i5ZIpvJNeeu7byy1rtKKE1iJM6gslwNKYX5Fz/moAb0gpvWanwJw/319MIBAIBLywjH9fCeCrATxPuX7lBHPP17HCa22X0Kslt3xQW1pZ/d+iKRRNtIZfdwpcO9t0BLp+LvWqaH3spVs9IeEL2BowpzB4uHadg+6UQtNZmuSIQdBrEKJxvb6+Vj6tPjzVn9vSfr1eHb0acs+o7zLhpS2mem8BNpVR38kn5ZzfSa+llG6ZXvTFo/Vier0epixX2YL2wo9yzLRTS4s47ZRLFs2nHPOWsBPSU/rihMQ7F8yLbQd+Tl9Qtz7NE2PERuBtP89LZ235ZEETxjS9Jvx4H/Cih+7i4Lywh6qwBHDLK4TXV4vb81GaSnvQyUPaNQ5tQf/z4zkXMQLwW9ilDaWwg8eI8aRXy5pSplRGi2O2kBbbgpFqybyOlftcLNYaNLA2ClIhfV4vciz5NktasrRTjKTNz7Kn2mpTN2t/QQme3TMA384cFb3CmKftNVhpXLxmWLag1ZNz5Py6RxhbGLnP3vQU2ke/NTpuCWN6PsvqcimlLwNwezn+YWzet5v82QcCgUCgF5bG/K0APr4cP59de8d+qrNfTNVe58prtFyPdZrXq2q0p+RJ0/QnRVuu6ejxaqVrzzUL3gz0nGvCtQ6SxiwtSGTdF78PLa41JOVxKnqH5RY0LZLTBRqsfqYt4mTlNUKHUbQ48nouHXuh0TcjHif7QJOy0DTmjrawOOabACCl9I6c863+LA8XmnCi4B1ihN+T8huNqw09e4d6S2zun3eiLfqC0BqcSrlG5m0vhDwqpP3SJE755NQXrxf0A0KnTXvAhXKv65h1XTOQaeAfEy/fPpXNgVC7AAAgAElEQVR+42lqOsvFb07BzDHSB/a1brbme+25/1kXys8535pS+mUAnwxCJV7FhfIDgUDgKsCzutwZivIB4IMAPhK4mgvlby2K44ivDZ2mGjFGtWjvRAMtjGo+VBPmVEZtp7zStWdN+9N2FV4wKkMyBmp190BrJ+9SqRJ9YVEPllboMSBaNIxX29NGGfSflmW1rdeYKWnJUttMXitCqT8/vkhqQ9KWezTmub0yTnLOKaWUc86PTindBuB/+os4HGwNDckQV3q2lgV7xLuDYg6uT4PHUk354iqkqVdGvUaF9PJMfxEXSttYLnHWztZeId36aNFwS0B76Isarn2Q5uSlaxWX0D0sdtIYbeaxTdDVBjW0OFLap7xpatl84wEKzavEwj53ONc+TvS61BbLjpXAuqqcUnohgPf0pAkEAoFAHzwa8wdTSs/C2nf5lSXsQ/ur0v4gGb8AsjaE8EVuDQNb2lBLk+NhVn6tukjaVaucqj1LmnG9zo+tYb2mMWvUhXau3RMtU6OJrGF8z3ZS0nBd8krg6Xs1Zq9WTMFHHFb6nj5s9cG0gJvE7FkTg44QAHvJTRrO66rd35zas0bhaEbSeo3+e+Ax/lUj361lv79PAfA1/iIOB60Ov8Tui9saBo4Mt/dFY0j10ThHzhdLAphTHHQBHu6VodWh5aHB40/lmK379ORtvXgjXgnSOa+nVG9gd+9CD01lXbP6sMYJz027Wc9J2s/QUmxaH2EJI9uzedah1j7aW/1kX1RGzvnenPPPAvjcnnSjSCm9NaWUy+/BEnY3CVullKpb38Mppfsuol6BQCCwT4xsxnohSCndAOAZAF4G4HsB/H5K6ZUAvhzAvTnnTyoeI78C4E948tS+vlvlVu3EiNODJpXQoZH0UBitNDy+pPGdncn0DzceapM3rKU9PfXU6i7tMOHRnCyDbgtVW24NYzXjYOseKKxnpi2J6knrhUYRLFdy/iPviaX98oWorHp5RqRSnfOqn9bgBmBOXyzPdsM1Q2ALByuYAbwOAHLOLwGAlNJDAF6Atdve00ucNwN4djl+BMDD3sytl16L7xmSUvTwxyPxpDDtJdVcBXl5p6eKcGFCmnLUNWy0bnMIk4qWkNbatiVwpWvUldDroUHr1rtn35S4vZDa62SxqTNVXKR6ULtNC2YbOalFi2/WwkeEMwV97tQ2c8Y8mGgcL9RqEbpg5zdyEwO4FdumhgcAfAQAkM1gf7pezDk/Nuf8yRdUt0AgENgbLI35v15YLWTsbfU6TSvl1lQKz/CoVY5kiOLHVl7ecE1TlsIkT5SK01Nsts8+heiVcW483KP2psHjbWAZ/yyDXQ3jBh1JK5JWWtPypFjsHOwHnhFYpVX4RBN1ZGlc24nXUY8ar6ZJhuZrPd+W9jylv9K+T597HUGdncmGwOUcE0xyzn+ps75z4zex9gCpuBFrquKGlNKNRWv+bCuDlBLtFyL4cJWG03/teusBj3hi9HqCTOlkaeEQzgBwtj2M1eCtM0+jcZfee5M4ey+Xr9EVwK4XRo8VvoUWd9uC1Act6saql8V/qnWrMzlX24LJytfq2+Jzg+5GZz3flpCu6Xu9NDhn3OSbgZ01zRnrsCJr32/q6avOpeALgXMvjJuw3hj21VhvpvH2EudZAO7XMsg5n+ScU855S/vmL5mGEaOGhhFOtfKyPfFHwDslLbf+Tk/X6zhfd7pZzJ5ek9JYhjbNkKLFA9YvUSadX8rbOu79LRvXan2sPKpmzX/83qZAajutHK2eFfzZWr+T8kukH/C+UPuLtYGuVHeK2s5Lof20NJ42lz7AHC1/bLFvnG0E8jWsB51n2AjoKpPKT9zn+GAFc9GIfwnASwH8PoAP5Jy/HMC3ALi5fHUWAJ55ebUMBAKB+XHIXhnIOT9dCHsJgJeM5Kd9Na2tf6SvpKWZjrpj9WJfnK41JKy4jnhuSLy8hzKgefPytCE6daeq6ysAa63NQ1lo4JqXRVGslGt0GG/1HYopfYW2Ryu9pVl66iaF07aoVAYdxdD2qBQE5xWtZ9bqg/X+tTxa1BhvPwvaSIRr0/X+q6YMrIf3K3LshWd1uQzgr+ac31zObwfwE5weuAqQBLNnG3kpH88DHhESFlodSeqs0hCNUxfabMeaZ6vsJYDWVN2p7SDdV2utZW/7c6FS24O7PW29fKSv8HhS/rU+EvIKQ2NXb5tqgllLn0i/4RScJAjP24IYiivoyoWAvYqdp50kvtnKi+dnCfKRTVRpG5zLEWyEcKU0gOYrsgVXd6hCuRy/qSP/QCAQCHTCRWWklP5BzvmfluNv2m+V9gc6vJqyjfxluIaNQKMOvJqB9z65V4Fm3NKGq6PtKSWjkzas2XX1ekVrDYyqLbfieWd3aZRNC9J6EhYsiqBCqrPmpbA1m1N5zgulDegyspzSaHntALtxvCPXFjXmqYOVr3YOrLXkjiUyzuERzG8D8J0ppe8o5wnA/xoo69LRWhmsJZhH/CFH3OV6hmFel7AWPJyiFofOCMNKvmd63DM1VSrf/Giy8zp0pov/S9D6QD2+dra7cQAX2lr+lHJptaeXQmulafUD6xlYW1pJZfBFrc7ONnlk4o1RPTkAnW+W0OqbXkjtN1U52Bea1ck5PxPAnQB+rfzuzDnftu+KBQKBwLFC1ZhTSq/IOX91Suk1JejXy/+dKaU7c86fv//qzYulYjmXtGTrSzoy7LEwZbjv1SStvDzD6pa3hEcL1rTKkaFeNxxGp3pM+0PdgJbP/FuttpdBlWb+VdTmSoP0hTg6gt5uU7Z14kuNase039FFrdw42x1laSPFipENaHvQ+2z2qWVbVMZXAvhqAM9Trl85wSy5NNFj75DwMoY9mgAedQ/zuvxJ4VL5dOgqta22eDiwPWybKqR7PpqtPlDpixquLU6zQ42RMhYkz5YnAa9LFXimN0cBpUt6vQu0CUb1mIZLx5xj5qjtpOWbO/qid43qy3hHz8uc8GE8z0u7QGakPInNVEkAnjS96EAgEAhI8Bj/fgu7CwpJYQcPaS4/1XwseLSdEU3NEy5d99AX2jWJltCMjHzBIx6X5+eBZc2eQ3v2DEm1OmiacNUKuTHw/BrNq/wvynGtBvVXbmmZnjpXcB9cb76e/kq1XC2NdS+r1XpCErC9JgvNp7VeC90FR6v/ZRnvxJEFNsJxtFoWx/xlAG4vxz9MyrhpsKxLBx1WU1hDv7l32OV14R3a+5HQXgQv12t9aBaCME6LtmdDD3j9Je5SWzB9znJ5mOYSV2kMzjkD2JpEsJUvdl9MOqy36iTVT3Pv6vk4SwKMTtapglgSetq7UINPiTSp90nX9z5l0qYqSnUmqdZ3+YYLtJ48fN+g5SZW/nkbroBKN9DZfuKiGAosjflbAXx8OX4+u/aOjjICgUAg0AFr2c+6l947cs63XlyV9gfN+GehfhS1DSJb8JSjxdE007OVrb22ylwstv17Txa72gyNa4VXzY1qEZJW09LkRc+DxbbGTvOwRg29miSfaCJ5XlRtecsYSPLzaPZWtznXMIXwuYy+NI+WJ440ROea63n59dmQjXspfUF3vuGauQeSZtwyTlKM+Eh7Iflo5wXO/bpP4Vsul8PDMd+cUnoWWyvjDTnnj+go5yAgccwWPC5qc6JnQSDN8mu53omuTijO/gJ/TTF1Gx4Oj2DhtEaFJKR5uhasraA4fQFshPL5NZDtk5QyuBGGeq9o4C55fIEmGq+HQ6+Q+oAFbmewvCEyuaZ9qKWJJ546aOVfFp0BrO9Fu+fzZczpe9VRhkcwX8/XykgpPaqjjINBi1/kD3llvPxa/pYw54akesy1Sa/bnnTNE2/HeOfIgxuYaDhvJ6mzSvvFeTRbCfwjwdtPe1aabYHfuzS7rwrl+m1fYnuhGgnnHz3tRjoxt8bcixavK2nD151u2pnyzdb6zJ7yNSOimE4I0/q4J55WNzrqrLaIU5JHj7eEq2lSSjeT4yd35B8IBAKBTng05gcBvCul9IdYKwlPLGFXDlOoDL4044hLVgt8dpmVh0er9lAZ9Suv7QjiHfpq2q+k2VQtcor2ZsHjmaJRRnw0Q7VnTl/ULHhx/HZWSrgH2pKgXipCGxHy53TifM5WXhpNIVEhdGTlhaSx07zVkZKUV1/RW2XSfnte/mq7r1eO/drZpqyTDpW5KZhzzh+VUnobgD+NsoDRVV0rY2ToV7FYtF2dRoaXmnuWds7R4263ZC8QsBY8p2SjVco3804oHUvX6JCWxtmqC+Yb4vdA8vfVaKatPd2wS1/UrOiLf4JtQd3jIqVBeqaa8c5LefHncd4ugw9F6/uUV14yKkJaaY7Xk9bVK5R5e9BrtS3pmtO9oEoI/RhJ5dNNJRYdncG17GfO+Rn+LAOBQCAwBZ4dTN4I4HN4+FXcwaTlwN4y3vUO92pcKcyj2Vgas3QvXg2Au/7RcqgXgEWFcPRoEVI9Kjz30PIOkTQkDo8xNq822u817NIXUv2XmEdL7sGIlqw923qNe4ZY4NqyZtjjVETNl+9q4qHM+GxALY2mPS9X07RmnndaYGuijdS3vGulAD6N+XMAfDuAFwF4DIB3AvgYfxGHA22tXO/W7XS41+vCZAlj/hA1waytztZjOZbK3hoWE06z5wMklUE7q/ai1nIqRmf6aZ2evogWaJtbU60pfSG5yUnaCm0DKkgkNzIvLFrCY5fQ0vBND6x8PV41lsDc8tBwlMPhXTdaQxXOgO+dsqhK7nUk1St1qLKuW8k5/8Py/1DO+eMBXO8vIhAIBAI98G4tdROAnFJ6BMAH91ul/UHSGKSlE1tprOte7dXSfqVhNY3r8drg2NKKQb7ui13NmBo5vf6eEqgWQYd3dYNOT917dtoYhUYRUboiC+f7QKttvQa/UUN3XWqU9sHzyUeGZj7FsD4FPVqy5HOvGQM9tJ2GJcjIiLTl3Ma/FwF4LIC/COAnATwaa2rjSkMT0nwoLXKjBm+nDRF5Z2/ViYNPFfam87q70bxGBJ6ns0oc3D7q4oUmWERBN5C/VHVKX3C+tAW6+FSvAuBFVQak/fx4XnO7OdJ8vfcnfajmLKPlgdRKsyWk5+SYc84vL4e/gbVQDgQCgcAeYS37aY7WrqJXxgjo11Gyro5SGVIZ/Nzasuiyho4SPFoEYBtHeu6nNbz0auPNOP4qbdeFHJ+Qc+/WSFM0Ue+StpqRtVIX0uhMm2TFYY3uPOgZDfSOqDhFwb00gN0242l6RqHed4PD0phf5M/masPTcPQhUMrDS2W0yu6ppxVuDa/o0FlzyB8VCj2dtdIaHqt+heceW+BbQGlxPEjY5pmrlnLCwuh9njBPDMtTRapTy4OIxm3d5wIbNzW6xkOlLuhqd/T9sGi8KXsNWpDeT2lRrV4ByJUtyjdTaBtEeDHilWEt+/ly7VogEAgE9oejmmCioXeyiDZs5sMWy6A0xy4gtNyecHpNmu46xeDmHroRDcWK5vVCmAKXEYidL7Hts3xC4i1YWAtUk+6BtRa3Z6RWy+TrpKxWm7rveGjUiAb9RNPQ/Hn4CHh/kDT4HorB0+4a7dP7yHreK0/UOsEEWE8weTeARzrrhLTGtZRSLr/XkGvvK2G/Xc5vIfFySulhEvduEr4qrnxIKT2cUrqvt14j8FAKmgvTPod8FXzyQv3V7aE4jcEnhYyWJw3T6bX6Oz3F+aSKk8VmpljrR8vzYtnZ3lsL7kj3Wn4n7LcgvwRCY2DT7vx+6HZdHuTV9q+i9qv6q+tI12PpR+MvzzbtVMNqGTzvihqf/ni9at0kzPFxtVxNJWj9tPW77nT9q/2V9lvt2Vr91wNX9JkmmLwNaw19AeA/A3geAKSUbgPwUQAeB+CTStz/A+BTi1b+BACPKpo7AHwzgHvLtRWAX/FWoNVgIy/JknVa3jlGtIRarxPlwfN74emk+60C+cSZlxXO42hpvD9+nz3lc4y88No9ABvBSwUt/52S44RtYX1yKrf76IemggvMqT8qiGn+0lrhNT6FJMS9fd87erQUH62eQLtvU4WhKg2WQkH7rUe5qOiRL66odIJJSul9/uy38DQAv5PXuKPkexeA68r1azViifOL5fSJ5X+VUqrKyNNL2JuxmR7+CIBzzToQCASuKrwTTB6N6RNMFtjdxPX2nPOrUkofwHqN59+pF1JKN2Cz7vMjOee/nlJ6LgDknB8o4T8N4Nkl7LEDddr5iu1jNtlckDQrzpPxY7rYi0Y1AHK8Vtme+vZayulEFM0To95zizrSto/S6srT8Crz5T0hnFfNueZpeWJ4dn+mz5ff41yjNKktz70VWDk0jdammj1Fq5M0685bf+qx4oGnH2pxeL/j9TyfFcqe2ci70zPB5J2Yf4LJspTxGKHch7Cmpp8F4CdTSq8C8DM9mXNf7J5VyUaHl9ZD1eApa4lttx1L4HFha5U5wn9p5fIy82CnlPydrfazhrXS+ssWdhajYYJaymLBjmn7S5QRvQa028iyV/Djet7Dq3ND4IK135bxT/lI8LB67Nk2TYJVjlRvT7ikwHj6JxfErbrRd5WiustRuaQ5UajVSin9ekrpfnJOjXHvNO9kk+b+mgbrPs132/6pVh5lv8EPAvh8AK8p+d5YLn92I22qP099A4FAYN/wyCVLY34KgK9jYZ8H4BYA3+mswPnyoCmlXwbwtMIT/1i5/iopXUrp+QDuzzm/uewx+JEAfjXnXIX827E2FD4LwP1SHhI8NMWIUYbuDGJ94Ue8MWgdttZRNtJIOztI9MWckLRzHkYpCk8d+Dq9FKq2MoPx7zwvEFqilKlVm2rCvM3rBA6NxuiBpCVLHhMt6maJbc3vXCts0A4WvSRRHlZeHJImOjqa00aENNxDY/LZhYuVr21oeSO0pkll5Jy/i9Yx5/yjAJBScglmhj+LtYGvVvP1RtzbAXxJ2kyV+WDO+Wnl+FsA3F0EdAbwzIG6nEN78FyITRVmPQKDd/gKafqoBYtX5h10Creu1YXOpqLxrKJWLP5CCJcwsnGAiwckbU7bSJp5BshtLgkDS3hosIRfRRXK2nWpb1llU1pksdL7plZPLyjP7BVmve8kp5I8/uP8frUt0XgfOms8p2ZdvRFzzq4lQo30Oed8StT45xhxv5Sq+znnR5NrLyHhi5zzvVPqFQgEAocGS9jmlNIf5pyfQANTSv8P+1uOdq/o+Tp6zjmkRVWmbAHFyzy3lFN6Q8hrn0tl8rpIbULLF4eH5Z9ekvKz1mWg/xye9rWepTjE58NshfKSNGGJ5uhBa+hMtVquLXNa4Px+yD30UA8SzWA9D4/xtY6s6CJCHs3cA23EcrLYvaZhx7hc88Z2/2jVM3dITUsw/20Ary6UwRmL/yX+Ig4bXirDkwbQKQtrRtQUDs560XuGetIO2hReXtja7od+WKy1fWk8i2fuBa9/ayjPXzarDTQeUxIG9JqWltexHku8srQKoeY1IX34JKHiST/VnqJ5N2iL1vdAalP6YaSTRaxyrHs8W+n34M1Dg3rLOecfKlbDt2DtT/wggLcUCuEH+4sKBAKBgAceP+bPvIiKXCQ8w1jrmEIzAtHrFV7fzxZGLb0SqqdBz7ZaXHPS2qBFBWmUBtVCzvNVhuRT0BoN9WiIWl+RvDS0sjlalE0r3RxYrbC1tdSIV4YE/gw1Q3GLltJGI1IdLcMsj2uh1jstAMmxfY61cCYZ9K4qvMPyKeihLkbdi0ZhuSO1eF0p3NOWVpolcN7Bt3jQhT7ctgSm12NFg0brLBzPRXNVpP/n+XUKBS78NF7ZC6kNevLpLdPqT6uVvKM5n7BE4Z1AZXnC8I8oL8uiljygz6zLttRXTCAQCAT2jaPUmAGfD6cX3jUDWrSFN2wU0j1vWemxrTV4NNReT5caxu+L0hcegwqHpG1J1616ea71+tf2UmO8jH2Png4FtE9oXkea33hNLxlWW3E0rdtDpUir8VVMfW5HK5grRmgNTxqJb+sVvC03qRa4kJMEHg9bYMP9Wnsb7gsS9wjIFnDNk0PLl6fh13rQK6A951PK8oAKpq11p5318HjS9NaHptW8jrwfW6le0lKynrym8MSaMJ/LXe7DDlLn4ZwmP5YEG89L0/A0oTwqcLXrVsfV7o3nyY15AJsC3ngJ62W6XbuWvxeakc8S0Fpc6XzE8FMxKiynfNS4ViblO4f9hC5kRMvkwmtfthHt46zFrf+t5zn6Aaqw2l/DcD8ZSxYIBAKBfeGoNOal8BW2hvuWNk1hURRTqAzrumeoRWciWlSGxvdq5UzVyEa0CL5uhtczZITv9Xh1WB4ah8oLT31ueQWXKtdLA7T6oDVqkjRmyh1b64sno5/TCTu8fGkPQ+olw+s80h+OSjD3bANvURmW4YpihD/ujdOCtPKd1/VsAd0YqIFy1DRvb9t4uVvPkF0Svpp7lIev3KqfwUnyvDxuUt5nLbUnXwDIQktgWXVbLHR3So7W4vhSH2zRNFJ4i7KwqKzWHIQKbuDr3Uh5iMbrTxIIBAKBfeKoNOaKzLQNbT1ljcqo5xUj2q/HJUzTIrzrcQAbjUCjNXgddiiCWib6F+BpjUpaEyR6NSjpmkZTWO5zPfkDdnsC2yMICmk9iNbITDKKScvA8j7eC6ol0vponjpzoFez5M+Mnrf6ar0vOoqWZulalITHGEifWerYruOoBLMkgOhwhgsvDZa3gHfoxuuhlekR4C1fS6D9AZLqoL0ozc5IjjWhROvGuTotf231Ma2e1kLo3uGuB9pH22onLky1tYi9XLqVN4Xn48q5ZKqcUI8d74d6Tm8N7ZwKZutDpNmJuFDW+uPSIYxp/uGVEQgEAh8mOCqNmULSAiStsjXE7KE0vNqxFIdqUpbT/ShaHifUqOddY4G2rVVmr38oH3Z6tghqacVT/F57IPkE87WIafkaTTK1btZoiBr5Fso7UOsNtGdd9tbBG59rzL2GTa8BWqI1vNCMzS0clWCWOjXtbPSlsIbgNA/v8FJLL6FVT95BvNNHJcqG583r4aUzpnDuXDB7XgTLo8JDc8wVLoHbJVptSD+0klvgXMLYC41L3qoLtvsTRa/XggTPc/AIvVbbt7yGOI1hrRZ5/oHC7jME+jjmoDICgUDgwHBUGnOFZ3hYNcyWH7Cm0cyl3bTKmZvS4GVKkLbQ6hnmtSzeXEvh6NGQvc9hLmpAo7+sa1RDpbTG3FqTNcrhoIZZ1UuEHGteDSNoURYULV/0lgdSqw9aNJtm0D9ZyGuNL8IrYwzWsFMSRh54hmS9nZi7TNFhsDaMHHGZ4qBcPK+Px93LOm8JY0AXyHN6WGgYoay8XjrcflAhvtxC20qeHFLbjrQFXTeD5yHZaUbRy/l77AocVhtp8VrQKEBKa5wjBLOM3o6pzQyyXlLNhYfC6lTZEGwSNG3LwpybtUo8qKfDewSxBd62kpDm573ct2WMpde9L7LmCqdpdeLLXdMKeXvqMqLJas+2VX8JLQM6P7c0ZKBPIPf0hRGDX0u+9Mif4JgDgUDgwHBUGnMLo8NebdhFXXgA39d9a50Jctxr3ZfgXQ9iCiwNQ9OSrWGwdy/BXk+MFs2ihdG85hgKS3nxekrPTfMamuNRcqpMCt8qc7AtemmnHjpO0+xHYL1TGrUjaekx868BryDShKr14tMhNd2Ik6ez3NPqQ6VuN/V6Kz33jx0RxvsS1Jowtl54jdfmbU3DpWOrHI+QngNeIdV6zhYkXrq3LhR5RfI0FANvW1fMOexvYe7n2XIc6KF21DLGqhYIBAKBfeGoNOaRoVJL4+WaWz2WtGXPcJtD04A8RiqPtqwZzST0ena0tGWPAcbbTp776B1ye2mZVrm98GifrfIkSsxDh1kYoTL2januoi160bM8MGBTUed98xCpjJTSKwB8VTn945zzY8m1dwG4GcD7c84fXcJeC+DOEuUMwPU555xSugPAj5Gsn51zflNK6T4Aj8s5f4RaB8fL4xVmIwKjRwhycFqDgnNbU4TEyMfDC22hGKsOUyAtEtSC5jnSk87DcfdganuMfNwppK3GJHjq6XXdk+J510/2PmvPuuw90IT0oVMZ/wPAGwA8Ily7GcBjAHxUSum2EnYngJ/Auo6nAN5awl8P4KGccwLwEID/uM9KBwKBwEXjwjTmnPMPAfihlNLDjajXpZTuKmn+GgCklH4HwJ+p1wH8/XL8NQC+rxw/DFnon8PjFeH1T6bH3OBXy9JoDqscCy3DDjceamhZw1v19BhCKySNxOtD6jUQ8bb1Ui4eWsOzDgQvjy/uI2Hu0Ygn3yXknU5alI1Hu+3VDC2DoWVIG5nIolFBUydjce+Llr/3VfTKuA/Ag1hTHD+fUvp6dv2dAD6xnuScX1X/U0rfV44/uVWI52XwuPBoQpbyyjVOy3OA52u9CNqwy+rILVj0Chd2PQK5hamc5NxWe2lCQWvRms0FuV5nq+12o89miufFaBpg83Hna1prnilbAgfKZJdBCs3Dy1rxtLx6y+Z5c4E9x6zZXhyEYM4538SCLqEpAoFA4DCwV8GcUrofwOPK6VfknP+1M+mPA/ib5PwWEL0kpXRX0ZbvapSf6fmI76THR1ZaB1YyELa05taQUFpWsKaT8vQa1yzDpmeChweedZm9eY+Ub60jwqFteeX1XNhqT+happZmDni8NrStrSg0KqFlCPTWQyunt0/3xqnwlDnHMqYUVC4VW9kO9iqYc84fM5juVSml70spvRHA52JNY/xSuXwNwL8A8Kryf83I5/ymU0q5VzC3BDFgu8TRdHO8eDQPiW/mgt0qs+eDAwzsIv1hBsmLxM2Fl3++43iL5x2hpDhMaozUz+sauNUHFe8DL3/t5Zinwuvm10szjdZNE8ZbdRvLuh8ppS8rX4pHYe19kVNKn2UkeT2Az8FaU14CeEYJ/3wAN5S8bgDweXusdiAQCFw4LtIr4/sBfH9H/Oco4T+KrgX0dtHjeUAHt+oAABAdSURBVNEzLduiBCyMTNXV4s5h/PPcc+8IwKIS9uWhQNEzHNU2iaX/FG6NDNvbc7U07l7vlFYdeDzJmGxRafxe+A4mGh3iMdhZ/XyKUdEbz6Pxa/fHw7Tyc5bDJRyE8e+i0BIsi8WYUPJ0nJH587188xzuclOFsqf8FjwfzsuA5+XThsccPet7e9zItPI0TxBA70OeuoCkkSZqePq6VxhfBD3m8Zjp+Sh7r2s4WsEshY1Oo+7VjL0uP96Xn+544dWYeZx9aMlT0VMebVu68NHchhvt5aQvtaVxSYsyeZ/ZiJHQ0ph7DWnqfZX/3pmBXmE8pd/1KCle7dhTzmTlZFryQCAQCMyNo9KYKaThO9/23KMxjwy5Rr+m0lecr8HcO4QEDlNbrqDPRoOmSc6tLUvltq576Ic56uHxwLE01LMB7dkaDSxZH5XKt+o30o9H3ivXiMBZTit8FRyzDKnz8tl6NC4VzBYv63HzmZsnkzozH7pbL6u1Rc8+h460jKkuYlr5XqpIw3nbGPWzuNt9wmsQlGAJGL4BbM9Hh9cps3xbG+pKHwxP35i73bUyRwXxFASVEQgEAgeGo9KYKbT571RDtjTmRLQDTXum7lYcvZoPN+JI6flmrN5FWubcnLUHmsFsqjF1rokKEg3Ro+W3+tBU9HhraPWrsNZh6TUS0ry1dSd4eZoGPYdBcMqzt9L25hvuch3QplHzlwqQhXnPmr9ed5serwqLb76MxVdGoN3vyB5vU92bNEic6sizkzjzOQT2FIqD1sfazqxVpqQ48Dppq/NpgloT+Pumj+YUyCO4Iq9uIBAIHA+OUmMe0SS5j3OFdwcE7+QEel0zOFpaRGun7EOANTzm99xLs3gt+xbNZEGjX6Q4+0TLd7pVhx5D2ohPbw+NJFVV8ok+EfpN6z7n1G6nvlc9dTkqwdyzKE+lMVpCXBOE9cX3CGR63UtjaHxzFc7AYQhoWheNl+fwrGhnCSYaRzvnx56XRmt3rZ6cCrPijsDimEdsGBXe1Qs9dfOUJ7ajkIYLa0lQtzDyTvS6Xc5BIQaVEQgEAgeGo9KYJezTak5hac5TfFIlzbHmf9m0hlYXqjnTcw5pwo90DPQNs600fNcSCyPaqCdurxYoGSVHjWQ1nWcD1ykYopHYuXdNaA2e59xTT80TZVOgP6+jFMxzuof1CL4prnOW29DUl4UvqDPXJAYuIGgZ3LWvYl/eCtp5DRvZm1DDHF4mEixeXusTPR4+tAzNS8NTn944VjqTO3dSf55yRq9X9NhNXPn1JwkEAoHAPnFUGvOc2hjVsLT1e/m5lHZEe7cMT3MMNXst+hwtbY1rz1bevc+pNSqRzqVn2cpvBFkZJbTK6zEae7wrLPB8ThbyBq69nhe98Wh8qZyWB9MIvH2HouUpRa/FWhkdaD1QbWEWft2DEb635ZIlvZQSfTDKNVu8rhTmddUaESCWoPeEa8IY8O2G3St8WsLDI1R4md7naPUPb/x6TmmN8/isjvyZeo6tMEsYa/Wcy0bUQzn2KEgx809B64W3HizXdqQX2dKYNUgLpu/bWNSTd+tYOvfmPZWT8whQS7BpwtgrPDivq9VlIfSbJbZnt1nufi2FgCoNLRcy70ezFY8LackIXY8rtGPp/rTV6UaEsSdOz8e41W9nMe73JwkEAoHAPnFUGnNFa3iuoaVhaZwyMI/bmndIOkWD1rThEe255T0ywlG24NXKvNqyx6vDCuegLlXS/n8cFv8thXMXMm1STMsDSEuzxT0Tmq9qzS1Io4CsvDO1TLpWtGeSGE9vjWymaso9tEYP7XlUgtkjXFo8FocljFeraQLZquOcxj+NSlks5DrQMMuYRfPVVuHjeVv1914bMexNpTIo+JBee1aU1rDg4b/P88WuPaRFMVi2AxetYbQNfz9q3Lxqt7lYT/KRWGKbwvHUX/rIWGVaddOuaf28h2MOKiMQCAQODKExK9crWpqS9nWVtOXWTLce9NIa2ixAvnsL15JpXjVcWtCp2XaLXYNO78wrzz174dHWpGsei31v/WjTSQZmj7a+pSGScE0b9xrRtHgr9ix7lgut1z0jFrOe5HzLGDlA2cwJtY+EV4YMrzDW0OpIvLNRUCHk8WftwYjA8vJzVGBbeyNa9ap126lnfanoi4hGmg54+FkOL988lWsG5OG2l25platx2TR+r7eAJvC898ynu3MKsPUBsoT0ArvCmefRomxo2CilZiGojEAgELjCOHqN2QtLW5a++NTIYWEf2zppWqY1206iMjT6gsfhebUMk5ZGssSmzbxpKDxD6BFYz3jf8GrNgO2FwA2CUttOeUc0aKNLK6yVl6j97hwIcYx3YB+0xiiOSjBXTHE459d5vJHh85zQXsyWB4F2TaIvJMFsebZ4sMWFlzA+BZjmfxkvkdfbg8KzApr3mVkC2vJCqKeJhw8I4V5qaaW0k8QxS23aovyke+AfIK8w1rxXPPfs+Uh30Vz+qIFAIBC4CByVxtyjIUhamaYl1/PL0Jbn0B41+gLQ6QtLY56C1Wp7bYa523GKxsfh8S9eQN8ZWqqXlFeP0dGiMjyGwFFNWqsvsH3fkqG8asvSs+aTZzgdp92rBGlU1+vjzNHSlLfuKbwybEzd2doKvwz6YgSel48PIzXh3aJDOOrLb1ncrbwukwvsXbdZoxK0PHh79fDLo3nvg1e26sGheWtUSB4m/NohYQ4ZcIC3FQgEAseNo9KYb/xTzzSv7/iQZsXvkoRTp/GVMByr1xdpOzwJQ6lFKusOpN1rlrbFJyLQOm3Vn4TT+vAyzw1+J5vwk9PNsVhHrX6KJrTTtnlT31Y71zjSffLnwa/z8vg1qRx6bStPUuetvJkmWtszJTLqSP7nq92DBvp8ktC/aN14eD1vgbbTcrnd95fLzfXlGXBW1NzVktB+y/UPKBpzbt9fWqzbEFhTGvQ+T05I25Lj0xNyz8n3HvH7s94vqT9o+KhP/CQ7AkHKPV7PVxgppeO40UAgcNDIOTc/fUcjmAM2UkrZ02GOGdFGNqJ95kNwzIFAIHBgCMEcCAQCB4YQzIGKK+Lod6mINrIR7TMTQjBfYaSUbkkpZfJ7uITfTcJWKaWbSvjN5bxeezHJblnC/gvJ/w9I3D8g4a8l4ddSWtvKU0rLlNI9F3P3PqSUvpa10QMl/AdZ+GeV8JvKfdRwej9vKWGPkPy72jqldE9KyTNT+0Iw0D6PZ+3zIMnud0vY+0j+Wl+5g+V/ewm/r/bjo0bOOX5X9AcgAbitHD8e67lFb8Rac3lXCT8D8Efl+I8AnJXjdwFYleNXAHgEwFNJ2B0lvzsA3F6Oby/XMoAfL+VnAG8r4UsA91x2u7A2ugvAa8rxp5f6vrL8/5cS/m4AD5fj+wEsy/EXlXjXl/MVgCcDuAbgFSSsp63vqfkfwm+gfe4l9/LUEu/FpF/cUNrkNhIm9ZVHADxYjh8E8KFyfF8t65h/l16B+M30IDcvyX8q/zeW8B9ff383L0k5vrGcJwDfA+BDRejUl+5eKkCK0L23vMiZhJ/HKwLrTZfdFo12ygB+od57Cfta0kZnAH6XxX9DOV4BuKUIle8mwqanrd8E4Nplt8OE9rm33jOAv1WO7yBpq2D+9EZfyQDuKsd3kfx/G8D7L7sdLvt36RWI38QHuH4Rcvl9CMBz2ctAX6oM4KvItQzgueX4rJz/XDl/P4APkLgfKGE/wvL/GXp+yD8ALy73+Fnl/9fJvdY2ek8RLDcCuLvEe0e59ovl/Fo5H2rrQ/052+fxpX1qn3sPSf8HJez95VztK7zPXJU+dFG/4JivOHLOD+W17+hfBXA9gBd2ZrEs+ZzmnFPO+TNLuOaPeiX7TErpNgDfirWQvQfAvwLwlDLxiM6AfRLWbfJeAC8tYWcAkHO+rbTRdSW8d+bswXDLHB3t8+ry/wQAXwzgY1NKLwWAnPNNpX0eW+Jcyb5yELjsL0P85vthrdX+MTqH10pe3VTGof6w4d/FITKAfwZFYyvpXqlc66YyLrstprYP1h+p32P94v8q6bqpjPiVtrrsCsRvwsMDng/gWeX4yaWz/wp0g9T9EAxSSt7V+Hc7NsY/yiW+kQimt192Wxj3kUp7XGPhX1T+ry/X313OnwrgGeX4F1oCY462vmLt82ARxgnAM8rz/wEjf7GvYNf498hlt8Uh/S69AvGb8PCAH8CG68sonDDWQ/AatgJwcwm/Bdv84Esb+b+bxL2PhP8ICT/DgWqCpa4/xdooA3g7gAfI+YewMXS9mMX9pkb+s7T1FWqf29h9/XEjf7GvYGM4rL87LrstDukXa2UEAoHAgSHI+UAgEDgwhGAOBAKBA0MI5kAgEDgwhGAOBAKBA0MI5kAgEDgwhGAOBAKBA0MI5sCFI6X0ArYkZk4p/WYjzYtH920sS0m+kJwvU0qvdKa9Z9/7Rdb8yz3eO5D+vpTS2fw1C1wWQjAHLhRlPd5/C+D+vF5XIQF4EdZrL+wLHwfgs+tJzvkk5/wVeyyvC3mzT97tAD7xMusSOAyEYA5cNN4CADnnj60BOeeX55wfB2w0Y/J7E8/AipNS+mMS/qGU0jvLpc8rYU8u/z9S4r+V5uW9iZTSz5N0y5TSzSU8p5QeJtdeUcKfz+LnlNJza5qS7V8gebyPa8J04f5yrdb5CSQOLSenlH7Se0+Bw0EI5sBF4ykoq7Up+I8AnlS0yO/EetU8V5yU0tsAPAbA08q15+ScbylpXl809N+omaSUno/1eg8vKfE/13MDKaVbAPx5AP+8pEsAfplE+WAJvx/AV5awVwOoKwH+kpL1fwPWGnTO+aON8u/CehTwBQA+Advv8b8H8NOlnG8A8CzPPQUOC73LFgYCs6JohCfA+ZD+8QB+pm5BpECL81Ss1wf+1ZLfGxvFv3gdLb/MGb/i7hL/68v5L2AtqCv+cfn/OWyE/QnWHxHknJ85kbf+ulLv/wAAZSunx5RrCcBn0/xTSp+ec/75CeUFLhihMQcuGv8bRCHIOZ8CuJNc/wmsNepHAXickocnjgej/b+V7t3l/9pg/hXWyMIS7E+p/H35hVC+YgjBHLhofAYApJT+iIR9Mjk+wXrPt0ewFuIStDi/AuDxKaWnljLuINeeKOTzLeto6W4hvoVvLvH/STn/NADvU2OvscSaWkBK6X8qce5j578B4CSldH1K6WtJ+D8HsEgp/a2U0uMBUNojAzjPP6X0agSuHi57ebv4Hd8P64XR6dKRGcBvl2vfge3lJnMJfzE5FuOUaw/SayXsHSSsrlv9I+Xa22g9hLrew+qZsV4Q/+fJ+RKb5T4zNtt1nW+thM3GrhlrTTrTNOX/8STO+7BZKzmT/3tK3PtYfnXt5+eytj277Ocdv/5fLPsZCFwAitfG+3LOD6SUfhzAs/PGTS4Q2EJQGYHAxeArALy3GOWeDeANl1yfwAEjNOZAIBA4MITGHAgEAgeGEMyBQCBwYAjBHAgEAgeGEMyBQCBwYAjBHAgEAgeGEMyBQCBwYPj/iHsOliODfooAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%time\n", "for i in range(1, len(gs)):\n", " yp, xp = np.unravel_index(\n", " np.nanargmax(resid_smooth.data), resid_smooth.data.shape\n", " )\n", " ampl = resid_smooth.get_by_pix((xp, yp))[0]\n", " gs[i].xpos, gs[i].ypos = xp, yp\n", " gs[i].fwhm = 10\n", " gs[i].ampl = ampl\n", "\n", " sh.thaw(gs[i].fwhm)\n", " sh.thaw(gs[i].ampl)\n", " sh.fit()\n", "\n", " sh.thaw(gs[i].xpos)\n", " sh.thaw(gs[i].ypos)\n", " sh.fit()\n", " sh.freeze(gs[i])\n", "\n", " resid.data = sh.get_data_image().y - sh.get_model_image().y\n", " resid_smooth = resid.smooth(width=6)\n", " resid_smooth.plot(vmin=-0.5, vmax=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generating output table and Test Statistics estimation\n", "When adding a new source, one need to check the significance of this new source. A frequently used method is the Test Statistics (TS). This is done by comparing the change of statistics when the source is included compared to the null hypothesis (no source ; in practice here we fix the amplitude to zero).\n", "\n", "$TS = Cstat(source) - Cstat(no source)$\n", "\n", "The criterion for a significant source detection is typically that it should improve the test statistic by at least 25 or 30. The last excess fitted (g5) thus not a significant source:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=6\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
delstatglonglatsigma
float64float64float64float64
156.52299.66-0.103650.18292
133.86300.150.138470.049966
111.18300.570.142770.052608
68.202297.440.114420.054202
28.633298.27-0.298220.050465
18.816299.47-0.0279940.025942
" ], "text/plain": [ "\n", "delstat glon glat sigma \n", "float64 float64 float64 float64 \n", "------- ------- --------- --------\n", " 156.52 299.66 -0.10365 0.18292\n", " 133.86 300.15 0.13847 0.049966\n", " 111.18 300.57 0.14277 0.052608\n", " 68.202 297.44 0.11442 0.054202\n", " 28.633 298.27 -0.29822 0.050465\n", " 18.816 299.47 -0.027994 0.025942" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from astropy.stats import gaussian_fwhm_to_sigma\n", "from astropy.table import Table\n", "\n", "rows = []\n", "for g in gs:\n", " ampl = g.ampl.val\n", " g.ampl = 0\n", " stati = sh.get_stat_info()[0].statval\n", " g.ampl = ampl\n", " statf = sh.get_stat_info()[0].statval\n", " delstat = stati - statf\n", "\n", " geom = resid.geom\n", " coord = geom.pix_to_coord((g.xpos.val, g.ypos.val))\n", " pix_scale = geom.pixel_scales.mean().deg\n", " sigma = g.fwhm.val * pix_scale * gaussian_fwhm_to_sigma\n", " rows.append(\n", " dict(delstat=delstat, glon=coord[0], glat=coord[1], sigma=sigma)\n", " )\n", "\n", "table = Table(rows=rows, names=rows[0])\n", "for name in table.colnames:\n", " table[name].format = \".5g\"\n", "table" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "* If you look back to the original image: there's one source that looks like a shell-type supernova remnant.\n", " * Try to fit is with a shell morphology model (use ``sh.shell2d('shell')`` to create such a model).\n", " * Try to evaluate the ``TS`` and probability of the shell model compared to a Gaussian model hypothesis\n", " * You could also try a disk model (use ``sh.disk2d('disk')`` to create one)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What next?\n", "\n", "These are good resources to learn more about Sherpa:\n", "\n", "* https://python4astronomers.github.io/fitting/fitting.html\n", "* https://github.com/DougBurke/sherpa-standalone-notebooks\n", "\n", "You could read over the examples there, and try to apply a similar analysis to this dataset here to practice.\n", "\n", "If you want a deeper understanding of how Sherpa works, then these proceedings are good introductions:\n", "\n", "* http://conference.scipy.org/proceedings/scipy2009/paper_8/full_text.pdf\n", "* http://conference.scipy.org/proceedings/scipy2011/pdfs/brefsdal.pdf" ] } ], "metadata": { "anaconda-cloud": {}, "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.6.0" }, "nbsphinx": { "orphan": true }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }