{ "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://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/master?urlpath=lab/tree/spectrum_analysis.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/docs/tutorials).\n", "- **Source files:**\n", "[spectrum_analysis.ipynb](../_static/notebooks/spectrum_analysis.ipynb) |\n", "[spectrum_analysis.py](../_static/notebooks/spectrum_analysis.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Modeling and fitting\n", "\n", "\n", "## Prerequisites\n", "\n", "- Knowledge of spectral analysis to produce 1D On-Off datasets, [see the following tutorial](spectrum_analysis.ipynb)\n", "- Reading of pre-computed datasets [see the MWL tutorial](analysis_mwl.ipynb)\n", "- General knowledge on statistics and optimization methods\n", "\n", "## Proposed approach\n", "\n", "This is a hands-on tutorial to ~gammapy.modeling, showing how the model, dataset and fit classes work together. As an example we are going to work with HESS data of the Crab Nebula and show in particular how to :\n", "- perform a spectral analysis\n", "- use different fitting backends\n", "- acces covariance matrix informations and parameter errors\n", "- compute likelihood profile\n", "- compute confidence contours\n", "\n", "See also: [Models gallery tutorial](models.ipynb) and docs/modeling/index.rst.\n", "\n", "\n", "## The setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:31.501039Z", "iopub.status.busy": "2020-10-16T10:10:31.500372Z", "iopub.status.idle": "2020-10-16T10:10:32.850119Z", "shell.execute_reply": "2020-10-16T10:10:32.850573Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from astropy import units as u\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as st\n", "from gammapy.modeling import Fit\n", "from gammapy.datasets import Datasets, SpectrumDatasetOnOff\n", "from gammapy.modeling.models import LogParabolaSpectralModel, SkyModel\n", "from gammapy.visualization.utils import plot_contour_line\n", "from itertools import combinations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model and dataset\n", "\n", "First we define the source model, here we need only a spectral model for which we choose a log-parabola" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:32.860666Z", "iopub.status.busy": "2020-10-16T10:10:32.860119Z", "iopub.status.idle": "2020-10-16T10:10:32.863807Z", "shell.execute_reply": "2020-10-16T10:10:32.863330Z" } }, "outputs": [], "source": [ "crab_spectrum = LogParabolaSpectralModel(\n", " amplitude=1e-11 / u.cm ** 2 / u.s / u.TeV,\n", " reference=1 * u.TeV,\n", " alpha=2.3,\n", " beta=0.2,\n", ")\n", "\n", "crab_spectrum.alpha.max = 3\n", "crab_spectrum.alpha.min = 1\n", "crab_model = SkyModel(spectral_model=crab_spectrum, name=\"crab\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data and background are read from pre-computed ON/OFF datasets of HESS observations, for simplicity we stack them together.\n", "Then we set the model and fit range to the resulting dataset." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:32.870777Z", "iopub.status.busy": "2020-10-16T10:10:32.868616Z", "iopub.status.idle": "2020-10-16T10:10:33.155628Z", "shell.execute_reply": "2020-10-16T10:10:33.155036Z" } }, "outputs": [], "source": [ "datasets = []\n", "for obs_id in [23523, 23526]:\n", " dataset = SpectrumDatasetOnOff.from_ogip_files(\n", " f\"\$GAMMAPY_DATA/joint-crab/spectra/hess/pha_obs{obs_id}.fits\"\n", " )\n", " datasets.append(dataset)\n", "\n", "dataset_hess = Datasets(datasets).stack_reduce(name=\"HESS\")\n", "\n", "# Set model and fit range\n", "dataset_hess.models = crab_model\n", "e_min = 0.66 * u.TeV\n", "e_max = 30 * u.TeV\n", "dataset_hess.mask_fit = dataset_hess.counts.geom.energy_mask(e_min, e_max)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting options\n", "\n", "\n", "\n", "First let's create a Fit instance:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:33.161277Z", "iopub.status.busy": "2020-10-16T10:10:33.159863Z", "iopub.status.idle": "2020-10-16T10:10:33.161914Z", "shell.execute_reply": "2020-10-16T10:10:33.162491Z" } }, "outputs": [], "source": [ "fit = Fit([dataset_hess])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By default the fit is performed using MINUIT, you can select alternative optimizers and set their option using the optimize_opts argument of the Fit.run() method.\n", "\n", "Note that, for now, covaraince matrix and errors are computed only for the fitting with MINUIT. However depending on the problem other optimizers can better perform, so somethimes it can be usefull to run a pre-fit with alternative optimization methods.\n", "\n", "For the \"scipy\" backend the available options are desribed in detail here: \n", "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:33.201681Z", "iopub.status.busy": "2020-10-16T10:10:33.185259Z", "iopub.status.idle": "2020-10-16T10:10:33.733055Z", "shell.execute_reply": "2020-10-16T10:10:33.733522Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "No covariance estimate - not supported by this backend.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 566 ms, sys: 4.11 ms, total: 570 ms\n", "Wall time: 562 ms\n" ] } ], "source": [ "%%time\n", "scipy_opts = {\"method\": \"L-BFGS-B\", \"options\": {\"ftol\": 1e-4, \"gtol\": 1e-05}}\n", "result_scipy = fit.run(backend=\"scipy\", optimize_opts=scipy_opts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the \"sherpa\" backend you can choose the optimization algorithm between method = {\"simplex\", \"levmar\", \"moncar\", \"gridsearch\"}. \n", "Those methods are described and compared in detail on http://cxc.cfa.harvard.edu/sherpa/methods/index.html. \n", "The available options of the optimization methods are described on the following page https://cxc.cfa.harvard.edu/sherpa/methods/opt_methods.html" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:33.739600Z", "iopub.status.busy": "2020-10-16T10:10:33.739039Z", "iopub.status.idle": "2020-10-16T10:10:34.388033Z", "shell.execute_reply": "2020-10-16T10:10:34.386982Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "No covariance estimate - not supported by this backend.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 647 ms, sys: 0 ns, total: 647 ms\n", "Wall time: 645 ms\n" ] } ], "source": [ "%%time\n", "sherpa_opts = {\"method\": \"simplex\", \"ftol\": 1e-3, \"maxfev\": int(1e4)}\n", "results_simplex = fit.run(backend=\"sherpa\", optimize_opts=sherpa_opts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the \"minuit\" backend see https://iminuit.readthedocs.io/en/latest/reference.html for a detailed description of the available options. If there is an entry ‘migrad_opts’, those options will be passed to [iminuit.Minuit.migrad](https://iminuit.readthedocs.io/en/latest/reference.html#iminuit.Minuit.migrad). Additionnaly you can set the fit tolerance using the [tol](https://iminuit.readthedocs.io/en/latest/reference.html#iminuit.Minuit.tol\n", ") option. The minimization will stop when the estimated distance to the minimum is less than 0.001*tol (by default tol=0.1). The [strategy](https://iminuit.readthedocs.io/en/latest/reference.html#iminuit.Minuit.strategy) option change the speed and accuracy of the optimizer: 0 fast, 1 default, 2 slow but accurate. If you want more reliable error estimates, you should run the final fit with strategy 2.\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:34.395587Z", "iopub.status.busy": "2020-10-16T10:10:34.393540Z", "iopub.status.idle": "2020-10-16T10:10:34.597443Z", "shell.execute_reply": "2020-10-16T10:10:34.596853Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 201 ms, sys: 72 µs, total: 201 ms\n", "Wall time: 200 ms\n" ] } ], "source": [ "%%time\n", "minuit_opts = {\"tol\": 0.001, \"strategy\": 1}\n", "result_minuit = fit.run(backend=\"minuit\", optimize_opts=minuit_opts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit quality assessment\n", "\n", "There are various ways to check the convergence and quality of a fit. Among them:\n", "\n", "- Refer to the automatically-generated results dictionary" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:34.601668Z", "iopub.status.busy": "2020-10-16T10:10:34.601092Z", "iopub.status.idle": "2020-10-16T10:10:34.605266Z", "shell.execute_reply": "2020-10-16T10:10:34.605869Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "OptimizeResult\n", "\n", "\tbackend : scipy\n", "\tmethod : L-BFGS-B\n", "\tsuccess : True\n", "\tmessage : b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n", "\tnfev : 148\n", "\ttotal stat : 75.60\n", "\n" ] } ], "source": [ "print(result_scipy)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:34.609677Z", "iopub.status.busy": "2020-10-16T10:10:34.609120Z", "iopub.status.idle": "2020-10-16T10:10:34.616951Z", "shell.execute_reply": "2020-10-16T10:10:34.616443Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "OptimizeResult\n", "\n", "\tbackend : sherpa\n", "\tmethod : simplex\n", "\tsuccess : True\n", "\tmessage : Optimization terminated successfully\n", "\tnfev : 170\n", "\ttotal stat : 30.35\n", "\n" ] } ], "source": [ "print(results_simplex)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:34.620819Z", "iopub.status.busy": "2020-10-16T10:10:34.620237Z", "iopub.status.idle": "2020-10-16T10:10:34.624965Z", "shell.execute_reply": "2020-10-16T10:10:34.624445Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "OptimizeResult\n", "\n", "\tbackend : minuit\n", "\tmethod : minuit\n", "\tsuccess : True\n", "\tmessage : Optimization terminated successfully.\n", "\tnfev : 52\n", "\ttotal stat : 30.35\n", "\n" ] } ], "source": [ "print(result_minuit)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Check that the fitted values and errors for all parameters are reasonable, and no fitted parameter value is \"too close\" - or even outside - its allowed min-max range" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:34.635818Z", "iopub.status.busy": "2020-10-16T10:10:34.635176Z", "iopub.status.idle": "2020-10-16T10:10:34.638655Z", "shell.execute_reply": "2020-10-16T10:10:34.639233Z" } }, "outputs": [ { "data": { "text/html": [ "Table length=4\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namevalueunitminmaxfrozenerror
str9float64str14float64float64boolfloat64
amplitude3.814e-11cm-2 s-1 TeV-1nannanFalse3.140e-12
reference1.000e+00TeVnannanTrue0.000e+00
alpha2.196e+001.000e+003.000e+00False1.122e-01
beta2.265e-01nannanFalse6.962e-02
" ], "text/plain": [ "\n", " name value unit min max frozen error \n", " str9 float64 str14 float64 float64 bool float64 \n", "--------- --------- -------------- --------- --------- ------ ---------\n", "amplitude 3.814e-11 cm-2 s-1 TeV-1 nan nan False 3.140e-12\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\n", " alpha 2.196e+00 1.000e+00 3.000e+00 False 1.122e-01\n", " beta 2.265e-01 nan nan False 6.962e-02" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result_minuit.parameters.to_table()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Plot fit statistic profiles for all fitted prameters, using ~gammapy.modeling.Fit.stat_profile(). For a good fit and error estimate each profile should be parabolic" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2020-10-16T10:10:34.652271Z", "iopub.status.busy": "2020-10-16T10:10:34.651636Z", "iopub.status.idle": "2020-10-16T10:10:35.236577Z", "shell.execute_reply": "2020-10-16T10:10:35.237225Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEWCAYAAABfdFHAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA74UlEQVR4nO3deXwU9fnA8c+TGwhHAuFKSLjvmwgo4q2AR73QSsHbUttqrffxs9VaW6u23lVrPbDe4lW1gieHigLhvu8A4UxIIIQjkOT5/TETXdYcm5Dd2WSf9+uVV3Z3Zmee2Z2ZZ77f+e73K6qKMcYYEyxRXgdgjDGmYbNEY4wxJqgs0RhjjAkqSzTGGGOCyhKNMcaYoLJEY4wxJqjCKtGIiIpIV/fxsyLyhzpc9iQRub+ulmeMMSYwYZVofKnqtar6ZwAROUlEcryOqZyITBORXBEpFJFFInJuFfPGu0lzh4jki8hHIpLqM/3PIrJEREpE5F6/957sTtstIrtE5H2/96aKyH/d5eaIyLWVxHC5m8Sv8XttnrsNOSLykIjE+L3vEhFZISL7RGSdiIz0mdZYRJ4WkTwR2SMiM32m3SoiS0Vkr4hsEJFb/ZZb6Tb7zfeS78WH+9rfRWSNu+yVInKZz7Tu7ueR634mn4pIj0C3WUReFZFt7vTVvp9XTYlIX3f9eSLykx+rich1IpIlIsUiMqm266lhTK1E5Ft3X9otIt+JyIgq5r9YRGaJyH4RmV6TZbn7/aMislVECtx9JdZvGVXtX6e63+9+93jL8HvvYBGZKSJF7rF1g8+0bBE54E4rEpHPfKbd5fN6kTtfmYi0CmCbR/q9t8jdPy90p4uI3C8iW9xjYrqI9An2NotIeiVx3VzZd1sVcc4709xtyPab1lpE3nC/1z3uPjCs2oWqatj8AQp0reD1k4Cco1z2JOD+OoqzPxDjPh4G7AXaVTLvbcAioA2QALwCvOcz/XJgDPBf4F6/97YB2ruP44GHgA99pk8DHgNigQFAPnCy3zKSgJXAUuAan9d/DYwE4oBUYB5wh8/004GNwHCcC5JUINVn+qvAm0AKEA0M8dvmwUAM0MNdziWBbLPPPMcDM/33CeBPQE83pmFAAXCcO20ocDWQ7H4mfwZW1mCb+wDx7uOewHbf7arhPtLDjeVc5zD7yfQLgPOAZ4BJNT1OahlTghtXFCDu+vPL9+UK5j8NuBj4IzC9JssC7gG+dr+LFOB74E+B7F9AK2APcJG7noeB733e2wrYCYzHOS6aAr18pmcDpwX4mdwLfBXINlfw3pNwjv0m7vOLga1AZ5xj4gFgfii22S+uTkAp0LGW+8lQ4FJgIpDtN60zcBPQzt3GiUAekFjlMgNY6R3AOvcDXQ6c7zPtCuBb4FFgN7AeOM59fbP7wVzuM/8k4Fngc3d5M4AM3wMI96Tizns/0AQ4AJQBRe5fe/wSB37JCBgEzHfX8xbOSdF3/rOBhW7cs4D+R/GlHASGVjL9GeAhn+dnAasqmO9VKjnputPj3R13ufs80f28UnzmeQ54xe99zwK/Aabjk2gqWP5NwEc+z2cBV1cybw+gEGgW4Gf0BPBkoNuMk6AW4CT0Ci8+fOb9ELi5kmnJ7vtbBrLNFWzjNuDio9lngK5UkRjcfXxSDfe5SpdXg2VEAee4n0/raua9hipOuhUtC8gCLvKZ5xfA5gD3r4nALJ/n5eeAnu7zv/rv537vzyaARIOTINfhc44KdJvdeV4CXvJ5fjvwts/zPsDBUGyz37LuAab5vXYVsALnwuxTfM67VSznNPwSTSXzFVLNBVkgVWfrcK4Cm+NcTb4qIu18pg8DFgMtgddxTujH4BxgE4CnRCTRZ/7xOFearXAO2teqWrmq7sO5+t2qqonu39aq3iMiccAHOKWHZGAycKHP9MHAi8Cv3Lj/BXwoIvHu9KdF5Olq1vGxiBwEZuOcxLMqmfUFYISItBeRxjjbP6WqZfutJ11EduPsdLfglGrAOUh8/5c/7uvz3qFAJk6yqc4JwDL3fdHu+1JEZK041UxPiUgjd95hOFdmfxKnamhJefVBBfELzv6zLIAYyt0IzFTVxVXN5MZzTBXLPgHYrqq7qph+xHvd734/TilwG/CJ+3qV+0x9IiKLcS6OPgSeV9WdQViW8NN9M01Emgewf/XBqQUAfjgHrHNfB6dEkO9Wce0Upzo63S+018SpQv1MRAZUEv5InFqDd2ux3Y2BscDLPi+/CXQVpwo3FqfkPtWdPxTbXO4y37hE5DzgLpxSdApOSfONmm5zRURkIE4NwdoqZ6zFldBC4Fz38RXAGp9p/XCuatr4vLYLGOg+ngS86TMtEaeI18F9/pMSjfv4JPyqzqiiRINzAtkKiN/VRPnyngH+7Le8VcCJNfwsYnGS4I1VzNMM50tVoATnSj25gvmqK9Ek41wxDfd57RvgSZyi9mCcqotV7rRonOR3rPt8OpWUaIArgRyglfu8vRtvFk4RuRVOyfUv7vS73On34uxkJ+KUNH9SlMe5OFmEWyVV3TYDHXB22ub++0QF738Z50CWCqalAVuAcYFss9+0aJyqu7uB2KPZZwjfEk0CMI4KruYrmLe6Es1PluVu17c4J7a2OBdk6u5P1e1fLwB/81vHt8AV7uPVOKXKY9x1PwF86zPvCKAR0Bi4E6cKtEUFcb9Q2WcfwDZfCmzgyHNMHPA4Px7rG4BOAR5TR7XNPu8ZiXMsJvq8NgWfkhROCXQ/1ZRqqKZEg3NuWwLcWd0+VG2JRkQuE5GF4tzw241zxdzKZ5YdPo8PAKiq/2u+JZrN5Q9UtQjn5Ni+ujhqqD2wRd1Pw7XR53EGcHP5Nrnb1aGmcajqYVWdAowSkZ9VMtszODtGS5zi8HvUoETjs658nBPrf+XHG9jjcepjN7vreQ3n5AlOddliVf2uquW6Vzt/A8aoap778gH3/5Oqus19/RHgTJ/ph3ES9yFVnYFzv+gMv2Vfh3N1dZaqFge4qY8B96nqnmrifhhnX7zY73tGRFKAz4CnVfUnV26VbPMPVLVUVb/BSVa/dl+udJ8RkfE+N2Fr/N1WR0SO91svvs9F5Hj3tWU+cYysapmqetD9bO6o4oo/IJUs6y84F1ULcS7yPsDZZ3ZS/f5VhHMS89UMpxoc9/3vq+pcVT2IczFznIg0d+P5VlUPqOp+VX0A5wR9xOfhliQu4sgSSU1cDvzHb9+7BycRdMA55v8EfOWWfoK6zX5xveueW8tlAI/77D/5OCXMVDmycUQgNR/AD5/fRzj3kR6obv4qE404rR7+DVyHU8/dAuemslT1vmp08Fl+Is6VepVVYThXAv724VyxlGvr83gbzofoG6dvMXMzzpVEC5+/xhWdlAIUA3SpZNoAnKumfPdk+yQwVNxWLrVYT2vcHVJVN6rq2aqaoqrDcJLZHHfeU4HzRWS7iGzHuXf2DxF5qnxhIjIa5/s9R1WXlL+uqgU4Cauizx2cqtIqichVOPf3TlXVmrQYPBV42CdugO9E5Bc+y/4TTknyDFUt9FtvEk6S+VBV/1JBXBVucyV8v9dK9xlVfU1/rNYdU4NtDYiqfuO7Xvc13zi+cV/r4xPH1wEuPhbnBm9d+GFZ7on+OlVNVdXOODUb89wkXt3+tQznuAFARJrgfA/l1ZyL/d5b/riy85JWMO0CnBPu9AC26wgi0gGnBuU/fpMGAG+pao6qlqjqJJzGOL1Dsc1VJM/NwK/89plGqjpLVf/qs89U2Gq1gu2Px7lw2IJTlVy9aopOvXHqX3vgVCdciVMkvMadfgXwjc/8P6kmwPlwj3cfT8K5cXQ8TjHzUY68AVZZ1VlPnIze3GfeX+LUoyfjJJnv+bHqLA7YBNyAc7K4gB+vwMGpK92Mc69BcEoaZwFNq/o8fGIZg1M0j8W5D3UIGFzJ/C/h1AE3d+e/C6e0VT49Fufq53Wc6oYEINqddgE/tuxJAd7myFYsvXBan8S5ceThNg4AWrifS/nfLJyb383d6afgHPwnVBL3fcBcnMSWhFOv+2efmNcCf3A/3xE4V17lNy7H41RXVNYqpqptbu0Xt+LUTzdyp98JrKGCVn44CXgO8FQl6610m931XoJT+o4GRuFczJxbm33GnScB5xhS93G8z/QY97UHcO4lJlBJ668Klq2BzFfB+4bz47HXCKcqdi9uy8YK5o9247oWpwVgAj9WJVa5LJwWVe3dz2G4+9mdEeD+lYLTAutCd50PcmQLrFNwbmoPdPelR4Gv3Wnp7v4Y5773ViAXvwYhOBcj99Vkm33muQvnHqL/e+/Bqc5ug3PMXuruQy2Cuc0+8/wCp+ZG/F4/H6eA0Md93hyfhhoVbEeUG8MYd3kJQJzPsfsRTqIJaH9V1YBanf0FJ/OXF/VmcHSJprzVWZH7RXbyPYCoING4z1/EOUnsxtmBE3BakxXiZPsbObLVWSZO0b281dlbfssb7X7pu3FKQJNxTxpujM9W8nn0wqlv3uu+dy5HtsQbCRT5PG+JU6W1053/G3xaqLnbqX5/V7jTrsep592Hc+J+kyNb6f0e5yDa5y43s4rvcTpHNm+ehnPRUOTzN8VneizwtBvzdpw64QSf6X2A79x1+7dG3ICT2H2X/Wwg21xB3D/sEz7Pi/2WfZc77XJ3+j6/6enVbTPOgT7D3d5CnLrnX/rFUuk+U0HcHSvYxmyf6fdWMP3egA7a2ieaE3Hul+3FOaZn4JN0cS4Qlvk8v6KCGCcFuKwTcFp/7ce5lzXeL5bq9q/TcC4kD+Dsux393v9rnCvqApwTX/l93j4454N9OOeLL/E7LnCSYAkV/5Si0m32mWclFbQewzkn/dPdNwpxWr2ODvY2+0z/FL/7iD7TLsXZpwtxkv6LVewnJ1XwGUz3+d7V/V59j6ORVe174r45JMT5YVqOqt4dspUaY4zxVNj2DGCMMaZhsERjjDEmqEJadWaMMSbyWInGGGNMUMVUP4s3WrVqpR07dvQ6DGOMqVfmzZuXp6opXsfhK2wTTceOHcnKqqz7MGOMMRURkY3VzxVaVnVmjDEmqCzRGGOMCSpLNMYYY4LKEo0xxpigskRjjDEmqCzRGGOMCSpLNMYYY4KqwSWa+ZsK+Punq7wOwxhjjKvhJZqNBTw1bS3LtxZWP7Mxxpiga3CJZuyQNOJjonh1dtj9ONYYYyJSg0s0LRrH8bMB7flgwRb2HjzsdTjGGBPxGlyiAZgwPIP9h0p5f8EWr0MxxpiI1yATzYAOLeiX2pxXv9+IjbdjjDHeapCJBmDC8HRW7yhibnaB16EYY0xEa7CJ5mcDUmmaEMOr31ujAGOM8VKDTTSN4qIZOySNKUu3kbu32OtwjDEmYjXYRAMwflgGh0uVt7M2ex2KMcZErAadaLq2TuTYzi15ffYmSsusUYAxxnihQScagEuPzWDL7gNMX7XT61CMMSYiNfhEc3rvNrRuGm+NAowxxiMNPtHERkdxyTEdmL46l835+70OxxhjIk7IEo2ItBCRd0RkpYisEJFjQ7XuccPSiRLh9TmbQrVKY4wxrlCWaB4HpqpqT2AAsCJUK27XvBGn9mzNW3M3U1xSGqrVGmOMIUSJRkSaAScALwCo6iFV3R2KdZebMDyD/H2HmLp0eyhXa4wxES9UJZrOQC7wkogsEJHnRaSJ/0wiMlFEskQkKzc3t04DOL5rKzJaNrZGAcYYE2KhSjQxwGDgGVUdBOwD7vCfSVWfU9VMVc1MSUmp0wCiooQJwzKYm13Ayu02KJoxxoRKqBJNDpCjqrPd5+/gJJ6QGjskjbiYKCvVGGNMCIUk0ajqdmCziPRwXzoVWB6KdftKahLH2f3b8f78LRQVl4R69cYYE5FC2erseuA1EVkMDAT+GsJ1/2DC8Az2HSrlAxsUzRhjQiJkiUZVF7r3X/qr6nmq6slAMYM6tKBP+2Y2KJoxxoRIg+8ZwJ+IMGF4Biu372XeRhsUzRhjgi3iEg3AuQPb0zTeBkUzxphQiMhE0zguhguHpPHJku3sKrJB0YwxJpgiMtEAjB+WzqHSMt7OyvE6FGOMadAiNtF0a9OUYZ2SeX3ORspsUDRjjAmaiE004DR13px/gBlr6ra7G2OMMT+K6EQzqk9bWiXG85o1CjDGmKCJ6EQTF+MMivblyp3kFNigaMYYEwwRnWjAGRRNgDdsUDRjjAmKiE80qS0acYo7KNqhkjKvwzHGmAYn4hMNOI0C8ooOMXWZDYpmjDF1zRINcEK3FNKTbVA0Y4wJBks0OIOi/WJYOnM25LN6x16vwzHGmAbFEo3r4swOxMVEWVNnY4ynVm3f2+B6lrdE40puEsdZ/drx7vwt7LNB0YwxHlizYy9jHp/Jy7OyvQ6lTlmi8TFheDpFxSX8d+FWr0MxxkSgB6eupElcDD8bmOp1KHXKEo2PwelJ9Gzb1AZFM8aE3JwN+XyxYifXntSF5CZxXodTpyzR+BARLj02g+XbClmwebfX4RhjIoSq8sCUFbRpFs9VIzp5HU6ds0Tj57yBqSTGx/Dqd9YowBgTGp8u286CTbu58bTuNIqL9jqcOmeJxk+T+BjOH5TKx0u2UbDvkNfhGGMauJLSMh6auoqurRMZOyTN63CCwhJNBSYMz+BQSRmT5232OhRjTAP3VtZm1uft47ZRPYiJbpin5JBtlYhki8gSEVkoIlmhWm9t9GjblKEdk3lt9iYbFM0YEzT7D5Xw2BdryMxI4vTebbwOJ2hCnT5PVtWBqpoZ4vXW2Pjh6WzctZ+v1+Z5HYoxpoF6/usN5O4t5s4zeyIiXocTNA2znFYHRvdtS8smcdb/mTEmKHYVFfOvGesY1acNQzKSvQ4nqEKZaBT4TETmicjEEK63VuJjovn5MR34csUOtu4+4HU4xpgG5smv1nKwpIxbR/X0OpSgC2WiGaGqg4ExwG9F5AT/GURkoohkiUhWbm5uCEOr2Lih6Sg2KJoxpm5t3LWP12Zv5OLMDnRtneh1OEEXskSjqlvd/zuB94GhFczznKpmqmpmSkpKqEKrVIfkxpzcozVvzt3M4VIbFM0YUzf+/tlqYqKiuPG0bl6HEhIhSTQi0kREmpY/Bs4AloZi3UdrwvB0cvcW89myHV6HYoxpABbn7OajRVu5+vhOtG6W4HU4IRGqEk0b4BsRWQTMAf6nqlNDtO6jcmL31qQlNeKV77O9DsUYU8+pKn+bspLkJnH86sTOXocTMiFJNKq6XlUHuH99VPUvoVhvXYh2B0X7fn0+a3faoGjGmNqbsTqXWet2cf0pXWmaEOt1OCFjzZsDcHFmB2KjhVe/t0YBxpjaKStzSjMdkhsxfliG1+GElCWaALRKjGdM33a8Oz+H/YdsUDRjTM19sHALK7fv5ZYzehAXE1mn3sja2qNw6bEZ7D1Ywoc2KJoxpoYOHi7lH5+tpl9qc87p397rcELOEk2AMjOS6NGmKa/OtkHRjDE188p3G9my+wB3jOlJVFTD7WqmMpZoAiQiTBieztIthSzK2eN1OMaYemLPgcM8NW0tJ3RPYUTXVl6H4wlLNDVw3qBUmsRFW/9nxpiAPTN9HYUHD3P76B5eh+IZSzQ10DQhlvMGpfLRoq3s3m+DohljqrZ19wFe+nYD5w1MpU/75l6H4xlLNDU0YXgGxSVlvDMvx+tQjDFh7rEvVqMKN53e3etQPGWJpoZ6tWvGkIwkGxTNGFOl1Tv28s68HC47NoMOyY29DsdTlmhq4dLhGWzI28e362xQNGNMxR6cspIm8TH89uSuXofiOUs0tTCmX1uSbVA0Y0wlZq/fxZcrd/Lrk7qQ1CTO63A8Z4mmFuJjorkoM40vVuxk+56DXodjjAkjqsoDU1bStlkCV43o5HU4YcESTS2NH5pBmaoNimaMOcLUpdtZuHk3N57ejYTYaK/DCQuWaGopvWVjTuyewhtzNtmgaMYYAA6XlvHQp6vo1jqRCweneR1O2LBEcxQmDMtg595ivlhug6IZY+DNuZvZkLeP20f3JCbaTq/l7JM4Cif3bE1qi0a8OtsaBRgT6fYVl/D4F2sY2jGZU3u19jqcsGKJ5iiUD4r27dpdrMst8jocY4yHnv96A3lFxdw+picikddxZlUs0Ryl8kHRXrNB0YyJWHlFxTw3cx2j+7RlSEaS1+GEHUs0RymlaTyj+rTlnXmbOXCo1OtwjDEeePLLNRwsKePWCO44syqWaOrAhOEZFB4s4aPFNiiaMZEmO28fr83exCXHdKBLSqLX4YQlSzR1YFinZLq1TrSeAoyJQA9/torY6ChuOLWb16GELUs0dcAZFC2DxTl7WJyz2+twjDEhsmjzbv63eBu/HNmJ1s0SvA4nbIU00YhItIgsEJGPQ7neUDh/cCpN42N48qu1XodijAkBp6uZFbRsEsfEE7t4HU5YC3WJ5gZgRYjXGRLNEmK59qQufL58B3M25HsdjjEmyKavzuX79fn87tRuJMbHeB1OWAtZohGRNOAs4PlQrTPUrhrRibbNEvjrJytQtbFqjGmoSsuUB6esJKNlY8YNTfc6nLAXyhLNY8BtQKUdg4nIRBHJEpGs3NzckAVWVxrFRXPTGd1ZuHk3nyzZ7nU4xpgg+WDBFlZu38stZ/QgLsZudVcnJJ+QiJwN7FTVeVXNp6rPqWqmqmampKSEIrQ6d+HgNHq0acpDn67kUIl1tmlMQ3PwcCmPfL6a/mnNOatfO6/DqRdClYpHAD8TkWzgTeAUEXk1ROsOqego4Y4ze7Jx135etz7QjGlw/vNdNlt2H+COMT2JirKuZgIRcKIRkQwRae7z/GQReVxEbhKRKoeQU9U7VTVNVTsClwBfqeqEWkcd5k7qnsJxXVryxFdrKTx42OtwjDF1ZM/+w/xz2jpO7J7CcV1aeR1OvVGTEs3bQBMAERkITAY2AQOAp+s8snpMRLhzTC/y9x3iXzPWeR2OMaaOPD3DuXi8Y0xPr0OpV2qSaBqpankfKxOAF1X1H8CVwNBAF6Kq01X17Bqst17ql9accwe25/mvN7BtzwGvwzHGHKWtuw/w0rfZnD8olV7tmnkdTr1Sk0TjWxl5CvAlgKraHe9K3HJGD1Th0c9Xex2KMeYoPfL5alC46fTuXodS79Qk0UwTkbdF5HEgCfgKQETaAYeCEVx91yG5MZcdm8E783JYub3Q63CMMbW0cnsh787P4fLjMkhLaux1OPVOTRLNDcB7QDZwvKqW3+VuC/xfHcfVYFx3SlcS42N4cMpKr0MxxtTSQ1NX0TQ+ht+e3NXrUOqlmiSaT1X1TVV9VFW3lL+oqgtU9dMgxNYgtGgcx29P7sq0VbnMWpvndTjGmBr6fv0uvlq5k9+c3JUWjatsYGsqUZNEUz9/QRkGLj+uI6ktGvHAlJWUlVnXNMbUF07HmStp1zyBK47r6HU49VZNeoJrLiIXVDZRVd+rg3gapITYaG4Z1Z0b31rER4u3cu7AVK9DMsYEYMrS7SzavJuHxvYnITba63DqrRolGuBsjmx9Vk5x7t+YSpw7IJV/z9zAw5+uYnTftsTH2E5rTDg7XFrGw5+uokebplw4OM3rcOq1miSajap6VdAiaeCiooS7zuzFhBdm88p3G7lmZGevQzLGVOHNOZvYkLePF6/IJNq6mjkqtf0djamF47u14oTuKTz51Vr27LeuaYwJV0XFJTz+5RqGdkrm5B6tvQ6n3qtJork0aFFEkDtG96Tw4GGenm4jcRoTrp6Zvpa8okPcOaYnInaNfbQCTjSqujSYgUSK3u2bccGgNF6alU1OwX6vwzHG+Fm0eTfPzljPBYNTGZSe5HU4DYKN2OOBm89wurB45DPrmsaYcHLwcCk3T15ESmI895zTx+twGgxLNB5o36IRV43oxPsLt7B0yx6vwzHGuB75fDVrdxbx0Nj+NG8U63U4DUaNE42IdBORd0RkuYisL/8LRnAN2a9P6kLzRrE8ONW6pjEmHMzNzuffX6/nF8PSOaG7/T69LtWmRPMS8AxQApwM/Ad4pS6DigTNG8Vy/Snd+HpNHjNX53odjjERbf+hEm6ZvIi0pEbcdWYvr8NpcGqTaBqp6peAqOpGVb0XZ9gAU0MThqfTIdnpmqbUuqYxxjN/m7KSTfn7eXjsABLja/LzQhOI2iSagyISBawRketE5HzAGprXQnxMNLeN6smKbYV8sGBL9W8wxtS5b9fm8Z/vNnLlcZ0Y3rml1+E0SLVJNL8HGgO/A4bgjLZ5WR3GFFHO6teOAWnN+cdnqzh4uNTrcIyJKIUHD3PbO4vpnNKE20b38DqcBqs2iaajqhapao6qXqmqFwLpdR1YpIiKEu4Y04utew4yaVa21+EYE1Hu/3g52/Yc4B8XDbBOM4OoNonmzgBfMwE6tktLTu3Zmn9OW0vBPhus1JhQ+GrlDt7OyuHaE7vYDzODLOBEIyJjRORJIFVEnvD5m4TTAs0chdvH9GRfcQlPTbOuaYwJtoJ9h7j93SX0bNuUG07r5nU4DV5NSjRbgXnAQfd/+d+HwKiq3igiCSIyR0QWicgyEflTbQNuqLq3acrFmR34z3fZbNplXdMYE0z3fLiMgn2H+MfFA2zIjhCoSV9ni1R1EtBFVV/2+XtPVQuqeXsxcIqqDgAGAqNFZHito26gbjy9O9FRwsOfrfI6FGMarE+WbOPDRVv53and6NO+udfhRISAG4yLyBKcAc4q7M1UVftX9l5VVaDIfRrr/tkPR/y0aZbAL0d25smv1nLN8Z0Y0KGF1yEZ06DkFRVz9wdL6Z/WnF+f1MXrcCJGTX6ZdPbRrEhEonGq2roC/1TV2RXMMxGYCJCeHpkN2Sae0JnXZ2/ir5+s4M2Jw62LcmPqiKpy13tLKCou4R8XDSA22rp6DJWaVJ1tLP9zX+rmPt4J5Afw/lJVHQikAUNFpG8F8zynqpmqmpmSEpl9DTVNiOWG07oxe0M+01bt9DocYxqMDxZu4bPlO7jljO50a9PU63AiSm061fwl8A7wL/elNOCDQN+vqruB6cDomq47Uowbmk6nVk144JOVlJSWeR2OMfXe9j0H+eN/l5GZkcTVx9sw6qFWm7Ljb4ERQCGAqq6hmi5oRCRFRFq4jxsBpwHWbXElYqOjuG1UD9bsLOLd+Tleh2NMvaaq3P7uYkpKlb9fNIDoKKuODrXaJJpiVf3hV4UiEkP1N/bbAdNEZDEwF/hcVT+uxbojxui+bRmc3oJHPl/N/kP2MyVjauvNuZuZsTqXO8/sScdWTbwOJyLVJtHMEJG7gEYicjowGfioqjeo6mJVHaSq/VW1r6reV5tgI4mIcNeZvdhRWMyL32zwOhxj6qXN+fu5/+PljOjakgnDMrwOJ2LVJtHcAeQCS4BfAZ8Ad9dlUMaR2TGZUX3a8OyM9eQVFXsdjjH1SlmZcus7ixARHho7gCirMvNMjRONqpbh3Pz/jaqOVdV/u7+TMUFw2+ieHDhcypNfrvE6FGPqlZe/y+b79fn88ezepLZo5HU4Ea0mfZ2JiNwrInk4N/JXiUiuiPwxeOGZLimJjBvagddmb2JD3j6vwzGmXlifW8SDU1dySs/WXJSZ5nU4Ea8mJZrf47Q2O0ZVW6pqMjAMGCEiNwYjOOO44dTuxMVE8fCn1lDPmOqUlik3T15EfEw0f7ugn/3oOQzUJNFcBoxT1R/uTKvqemzgs6BLaRrPr07owidLtjN/U3XdyhkT2Z6buZ4Fm3Zz37l9aN0swetwDDVLNLGqmuf/oqrm4vRdZoLompGdSGkaz1//twK7JWZMxVZt38ujn69mTN+2/GxAe6/DMa6aJJqqRuSy0bqCrEl8DDee1p2sjQV8tnyH1+EYE3YOl5Zx09sLaZoQw/3n9bUqszBSk0QzQEQKK/jbC/QLVoDmRxdnptElpQkPTlnJYeuaxpgj/HPaWpZtLeSvF/SjZWK81+EYHzXpVDNaVZtV8NdUVa3qLARioqO4Y0wv1uft4625m70Ox5iwsSRnD099tZbzB6Uyqk9br8Mxfqyf7HrmtF6tGdoxmce+WE1RsXVNY0xxSSk3T15Iy8Q47j2nj9fhmApYoqlnRIQ7z+xJXtEh/j1zvdfhGOO5Rz9fw+odRTx4YX+aN7bKlXBkiaYeGpSexFn92/Hvr9ezs/Cg1+EY45l5G/N5buY6xg3twEk9quxE3njIEk09dduoHhwuLeMx65rGRKgDh0q5ZfJi2rdoxP+d1dvrcEwVLNHUUxktmzB+WAZvzd3M2p17vQ7HmJB7cOpKNuTt4+GxA0iMr8mo9CbULNHUY9ef0pXGsdE8OHWV16EYE1Kz1uUxaVY2VxzXkWO7tPQ6HFMNSzT1WMvEeK49qQufL9/BnA35XodjTEgUFZdw6+TFdGrVhNtH9/Q6HBMASzT13FUjOtG2WQJ//ng5h0rsR5ym4fvL/5azbc8B/n7RABrFRXsdjgmAJZp6rlFcNPf+rDdLtuzhzx8v9zocY4Jq2qqdvDFnMxNP6MKQjCSvwzEBskTTAIzu245fndCZV77fyOQs6zHANEx79h/mjncX071NIjee3s3rcEwNWKJpIG4d1YMRXVvyfx8sZUnOHq/DMabO3fvRMnYVHeKRiwcSH2NVZvWJJZoGIiY6iifHDSYlMZ5rX53HrqJir0Myps5MXbqd9xds4bpTutI3tbnX4ZgaCkmiEZEOIjJNRFaIyDIRuSEU6400yU3ieHbCEHKLirn+jQWUWA/PpgHYVVTM/72/hL6pzfjtyV29DsfUQqhKNCXAzaraCxgO/FZE7Ke8QdAvrTl/Pb8fs9bt4qFP7fc1pn5TVe7+YCl7D5bwyMUDiY22Spj6KCTfmqpuU9X57uO9wAogNRTrjkRjh6Rx6fAMnpu5no8Xb/U6HGNq7cNFW5mydDs3ndGd7m2aeh2OqaWQXx6ISEdgEDC7gmkTRSRLRLJyc3NDHVqD8oeze5OZkcRt7yxm1XbrosbUP6u27+UPHyxlcHoLfjmys9fhmKMQ0kQjIonAu8DvVbXQf7qqPqeqmaqamZKSEsrQGpy4mCieHj+YJvEx/OqVLPYcOOx1SMYEbHP+fi59YTaN4qJ5/JJBREfZsMz1WcgSjYjE4iSZ11T1vVCtN5K1bpbAM+MHk1NwgBvfWkhZmXodkjHVyt1bzIQXZlNcUsYrVw+jQ3Jjr0MyRylUrc4EeAFYoaqPhGKdxpHZMZl7zunNVyt38rgNKWDC3J4Dh7nsxTnk7i3mpSuPsfsyDUSoSjQjgEuBU0Rkoft3ZojWHfEmDM9g7JA0Hv9yDV8s3+F1OMZU6MChUq55eS5rd+7l2QlDGJxuXcw0FKFqdfaNqoqq9lfVge7fJ6FYt3GGf77/vL70TW3GjW8tZEPePq9DMuYIh0vL+O3r88naWMCjPx/ICd3tHm1DYo3SI0RCbDTPThhCbEwUE/+Txb7iEq9DMgaAsjLltncW89XKndx/Xl/O7t/e65BMHbNEE0HSkhrz5LhBrMst4tZ3FqFqjQOMt1SV+z5ezvsLtnDrqB6MH5bhdUgmCCzRRJgRXVtx++iefLJkO/+aud7rcEyEe/KrtUyalc3Vx3fiNyd18TocEySWaCLQxBM6c1b/djw0dSXfrMnzOhwToV75LptHPl/NhYPT+L8ze+E0TjUNkSWaCCQiPHRhf7q2TuT6N+azOX+/1yGZCPPfhVv444fLOK1XGx68sB9R9oPMBs0STYRqEh/Dvy7NpKRMufbVeRw8XOp1SCZCTF+1k5vfXsTQjsk89YtBxFhHmQ2efcMRrFOrJjz284Es21rIXe8vscYBJujmbczn2lfn0aNtU/59eSYJsTaAWSSwRBPhTu3Vht+f1o335m/hle83eh2OacBWbi/kypfm0q55I16+aijNEmK9DsmEiCUaw+9O6capPVtz30fLmZud73U4pgHatGs/l74wh8ZxMbxy9VBaJcZ7HZIJIUs0hqgo4ZGfD6RDcmN+89p8dhQe9Dok04Ds3HuQCS/M5nBpGa9cPZS0JOskM9JYojEANG8Uy7MThrCvuITfvDafQyU2DLQ5ensOHOayF+aQV1TMS1ccQzfrJDMiWaIxP+jRtikPje3PvI0F3PfxMq/DMfXcgUOlXD1pLutyi/jXpUMYZJ1kRqwYrwMw4eXs/u1ZkrOHf81cT/+0Flyc2cHrkEw9dLi0jN+8No95mwr45y8GM7KbdZIZyaxEY37i1lE9GNG1JXd/sJTFObu9DsfUM2Vlyi2TFzFtVS5/Pb8fZ/Zr53VIxmOWaMxPxERH8eS4waQkxnPtK/PYVVTsdUimnlBV/vTRMv67cCu3je7BuKHpXodkwoAlGlOh5CZxPDthCHn7DnH9GwsoKbXGAaZ6j3+5hpe/28gvR3bi1ydaJ5nGYYnGVKpfWnP+en4/Zq3bxUOfrvI6HBPmJn27gce+WMNFQ9K4yzrJND6sMYCp0tghaSzavJvnZq6nX2pzzhlgg1KZn/pgwRbu/Wg5Z/RuwwMX9LMkY45gJRpTrT+c3ZvMjCRue2cxK7cXeh2OCTPTVu7klsmLGN45mSfGWSeZ5qdsjzDViouJ4unxg0lMiOFXr8xjz4HDXodkwsTcbKeTzF7tmvHvy6yTTFMxSzQmIK2bJfDM+MFsKTjA799cQFmZ9fQc6ZZvLeSqSXNJbdGISVceQ1PrJNNUImSJRkReFJGdIrI0VOs0dSuzYzL3nNObaatyeezLNV6HYzy0cdc+LntxDonxMbxyzTBaWieZpgqhLNFMAkaHcH0mCCYMz2DskDSe+HINXyzf4XU4xgM7Cp1OMkvLnE4yU1s08jokE+ZClmhUdSZgfdDXcyLC/ef1pW9qM258ayHrc4u8DsmE0J79TieZu4oOMenKoXRtbZ1kmurZPRpTYwmx0Tw7YQixMVHWOCCC7D9UwpWT5rAhbx//viyTAR1aeB2SqSfCKtGIyEQRyRKRrNzcXK/DMVVIS2rMk+MGsSFvH2c98TULNhV4HZIJokMlZfz61fks3LybJ8YNZETXVl6HZOqRsEo0qvqcqmaqamZKivX2Gu5GdG3FW786FlW46Nnv+NeMddYarQEqLVNunryIGaudTjJH97VOMk3NhFWiMfXPkIwkPvndSE7v3YYHpqzkiklzyd1rnXA2FFt2H+CKl+bw0aKt3DGmJ5dYJ5mmFkLZvPkN4Dugh4jkiMjVoVq3Ca7mjWN5evxg7j+vL9+v38WZT3zNN2vyvA7LHAVV5Y05mxj16EzmbSzgL+f35VrrJNPUkqiGZ1VHZmamZmVleR2GqaGV2wu57vUFrMst4jcndeH3p3Un1rokqVe27j7AHe8tYebqXI7r0pIHL+xPh+TGXodlAiQi81Q10+s4fFmnmqZO9WzbjA+vG8GfPlzOP6et47t1u3hi3CDSkuxEFe5UlbezNnP/xysoVeXP5/Zh/LAMoqKsg0xzdOxS09S5xnExPDi2P0+MG8TqHUWc+fjXTF26zeuwTBW27j7A5S/N5fZ3l9AntRmf/v4ELj22oyUZUyesRGOC5mcD2jMgrTm/e2MB1746nwnD07n7rN7W8WIY8S/F3HduHyZYKcbUMUs0JqgyWjZh8rXH8ffPVvHczPVkZRfw1C8G2S/Kw8C2PQe4490lzFidy/DOyTx04QDSW1oVp6l7VnVmgi4uJoq7zuzFS1cew869xZzz5Le8PXcz4doQpaErL8Wc8ehM5mzI508/68Pr1wy3JGOCxhKNCZmTe7Rmyg0jGZTegtveXcwNby5k70HrviaUtu05wJWT5nLbO4vp1a4ZU38/ksuPs3sxJris6syEVJtmCbxy9TCemb6WR79Yw6Kc3Tw5bhD901p4HVqDpqpMnpfDnz9eTkmpcu85vbnMbvabELESjQm56CjhulO68ebE4RwuKePCZ2bx/NfrrfuaINm+5yBX+ZVirhjRyZKMCRlLNMYzx3RM5pMbRnJyj9bc/78VXP3yXHYVWfc1dUVVeWdeDqc/OoPv1u/innN68+Yvh5PRsonXoZkIY4nGeKpF4zj+dekQ7ju3D9+u28WYx79m1jrrvuZobd9zkKtfzuKWyYvo1bYZU284gSutFGM8YonGeE5EuOzYjnzwmxEkJsQw/vnZPPLZKkpKy7wOrd4pL8Wc8egMZq3L449n9+bNicPp2MpKMcY71hjAhI3e7Zvx8fXH88f/LuOJr9by3fpdPH7JINrbUMEB2VF4kLveW8KXK3dyTMckHho7gE6WYEwYsBKNCSuN42L4+0UDeOznA1m+tZAxj3/NZ8u2ex1WWFNV3pufw+mPzODbdXn84ezevDnxWEsyJmxYicaEpfMGpTKgQwuuf2M+E1+ZxxXHdeSOMT2t+xo/vqWYzIwkHr7ISjEm/FiiMWGrU6smvPvr43ho6ipe+GYDczbk8+QvBtElJdHr0Dynqry/YAv3friM4pIy7j6rF1eO6ES03ew3YciqzkxYi4+J5g9n9+aFyzPZtucA5zz5De/Oy/E6LE/tLDzIL/8zj5veXkS3Nk2ZcsNIrhnZ2ZKMCVuWaEy9cGqvNky54QT6pTbn5smLuPGthRQVl3gdVkipKh8s2MLpj87k6zW53H1WL97+1bF0thKeCXNWdWbqjbbNE3j9l8N56qu1PP7laqYs3cbADi0Y2jGZYzolMzg9iSbxDWuXLiouYcGmAuZtLODbtXnMzS5gcHoLHr5ogFUhmnrDhnI29dKCTQV8tGgbc7PzWbZ1D2XqdG3Tp30zMjOSGdopicyOybRKjPc61ICpKjkFB5i30UksWRsLWLW9kDIFEejRpiljh6TZvRhTpXAcytkSjan3iopLmL+xgLnZ+czZkM/CzbspLnF+7Nk5pYlT4umYzNBOyaQlNUIkPE7Sh0rKWL6tkKzsfOZvKiAru4Cde50ueJrERTMoPYnBGUlkZiQxML0FzRJiPY7Y1AeWaGrAEo2preKSUpZu2cOcDU7yycrOp/Cgcz+nbbMEjumUzNCOSRzTKZnurZuGrFuWgn2HnITillgW+STE1BaNyOyYxJAM569Hm6bERNstVFNzlmhqwBKNqStlZcrqnXuZuyGfOdkFzN2Qz/bCgwA0S4gh84cSTxL9UlsQF3P0J3hVZV3uPuZvLCBrYz7zNhawLncfADFuFd+QjOQfEkvb5glHvU5jIMITjYiMBh4HooHnVfVvVc1vicYES/m9kDkb8p3qtux81rtJID4mymlg0MlJPoMzkkgMoIHBwcOlLNq8m3mbCpiXXcD8TQUU7HcGdWveKPaHhDIkI4kBaS1oFGc/PDXBEbGJRkSigdXA6UAOMBcYp6rLK3uPJRoTSnlFxWRl5zM3u8BtYFBIaZkSJU4fbMd0TP6hdVurxHh2Fh78oQosa2MBy7bsocQdT6dzShOGpCf9UBXWuVWi9ZpsQiaSE82xwL2qOsp9fieAqj5Q2Xss0RgvlTcrdqrb8lmw6cf7KclN4sjfdwiAuJgoBqQ1Z0hGMpkZzs375CZxXoZuIlw4JppQ/eggFdjs8zwHGOY/k4hMBCYCpKenhyYyYyqQGB/DyG4pjOyWAjgtxJZs2cPc7HzW7iyiZ9umDM5Iom/75nVyT8eYhixUiaaieoOfFKVU9TngOXBKNMEOyphAxcVE/XCPxRhTM6G6FMsBOvg8TwO2hmjdxhhjPBSqRDMX6CYinUQkDrgE+DBE6zbGGOOhkFSdqWqJiFwHfIrTvPlFVV0WinUbY4zxVsh6IFTVT4BPQrU+Y4wx4cGayxhjjAkqSzTGGGOCyhKNMcaYoLJEY4wxJqjCtvdmEckFNtbR4loBeXW0LC/U9/ih/m9DfY8f6v82WPyByVDVlBCsJ2Bhm2jqkohkhVvfPzVR3+OH+r8N9T1+qP/bYPHXX1Z1ZowxJqgs0RhjjAmqSEk0z3kdwFGq7/FD/d+G+h4/1P9tsPjrqYi4R2OMMcY7kVKiMcYY4xFLNMYYY4KqQSQaEUkQkTkiskhElonInyqYZ7yILHb/ZonIAC9irUwg2+Az7zEiUioiY0MZY1UCjV9EThKRhe48M0IdZ1UC3I+ai8hHPvNc6UWsVRGRaBFZICIfVzBNROQJEVnrHguDvYixKtXEH9bHcbmqtsFnnrA7joMlZL03B1kxcIqqFolILPCNiExR1e995tkAnKiqBSIyBufG3E+Gk/ZQINuAiEQDD+IMuRBOqo1fRFoATwOjVXWTiLT2KNbKBPId/BZYrqrniEgKsEpEXlPVQ55EXLEbgBVAswqmjQG6uX/DgGcIr+MAqo4/3I/jclVtQzgfx0HRIEo06ihyn8a6f+o3zyxVLXCffo8zymfYCGQbXNcD7wI7QxVbIAKM/xfAe6q6yX1PfdwGBZqKiACJQD5QErooqyYiacBZwPOVzHIu8B93W78HWohIu5AFWI3q4g/34xgC+g4gTI/jYGkQiQZ+KKouxPniPlfV2VXMfjUwJSSB1UB12yAiqcD5wLMehFetAL6D7kCSiEwXkXkiclnIg6xGANvwFNALZyjyJcANqloW2iir9BhwG1BZTKnAZp/nOe5r4eIxqo7fV1gex1SzDeF+HAdDg0k0qlqqqgNxrnCGikjfiuYTkZNxdtDbQxheQALYhseA21W1NNSxBSKA+GOAIThXe6OAP4hI99BGWbUAtmEUsBBoDwwEnhKRCqtHQk1EzgZ2quq8qmar4LWw+I1DgPGXzxuWx3GA2/AYYXwcB0ODSTTlVHU3MB0Y7T9NRPrjFGfPVdVdoY0scFVsQybwpohkA2OBp0XkvFDGFogq4s8BpqrqPlXNA2YCYXkzt4ptuBKn+k9VdS3OPYOeoY2uUiOAn7n7x5vAKSLyqt88OUAHn+dpOKWzcBBI/OF+HAeyDfXiOK5Tqlrv/4AUoIX7uBHwNXC23zzpwFrgOK/jre02+M0/CRjrddw1/A56AV/ilGwaA0uBvl7HXsNteAa4133cBtgCtPI69gq25STg4wpePwunukmA4cAcr2OtYfxhfRwHsg1+84TVcRysv4bS6qwd8LLbkiMKeFtVPxaRawFU9Vngj0BLnKsHgBINr55UA9mGcFZt/Kq6QkSmAotx6q+fV9Wl3oX8E4F8B38GJonIEpyT9e3qlM7Cll/8nwBn4pys9+OU0MJaPTuOK1SPjuOgsC5ojDHGBFWDu0djjDEmvFiiMcYYE1SWaIwxxgSVJRpjjDFBZYnGGGPqiIi8KCI7RaROWlOKyFQR2e3fOaeIXOd2jKoi0qou1hVMlmhM2Av04BWRdiLyWQWvNxKRGW6z5aOJo4OITBORFW7PzTfUYhn9xOm9eqGI5IvIBvfxFxXMe6+IPOD32kARWeE+/kJEkmq/RSYIJlHBj8WPwsPApRW8/i1wGrCxDtcVNJZoTH0wicAO3tFU3BvuVTi/5j/aLj9KgJtVtRfOjx1/KyK9a7IAVV2iqgPV6ebmQ+BW9/lpFcz+BvBzv9cuAV53H78C/KYm6zfBpaozcTpa/YGIdHFLJvNE5GsRCbgnCVX9EthbwesLVDX7qAMOEUs0JuxVdPBWYjQVd7I4Hvhv+RMRuU1Elogzpszf3Nemi8ijIjLTLbEcIyLvicgaEbnfjWObqs53H+/F6Qb+Jx1SishFIrLUXf7MQLZRRM4Qke9EZL6ITBaRRFVdBewWEd9u8C/G6doEnEQ1LpDlG089B1yvqkOAW3CGyogoDaVnABPh3GqxHqq63O/1OKBz+dWfOGOYnAcMU9X9IpLsM/shVT3BrRL7L04HoPnAOhF5VH361RKRjsAgoKJewv8IjFLVLeKMwVNd7K2Au4HTVHWfiNwO3ATch1OquQSYLSLDgV2qugZAnTFZ4kWkpYZfn18GEJFE4DhgstuTAUC8O+0CnO/Y3xZVHRWaCEPDEo1pKIZR8Um/FbDb5/lpwEuquh9AVX1LSh+6/5cAy1R1G4CIrMfpiHKX+zwRZyyR36tqYQXr/Banm5q3gfcCiH040Bv41j0ZxQHfudPeBGaJyM04CecNv/fuxOlJ2hJNeIoCdrtVpUdQ1fcIbP+o9yzRmIZiDDC1gtcPAAk+z4XKu8Uvdv+X+Twufx4DIM7Im+8Cr7knip9Q1Wvd6q6zgIUiMrCaEofgjH3zk2owVd3s9vJ7InAhcKzfLAk422jCkKoWug0+LlLVyeJcSfRX1UVexxZKdo/GNBSn4vQMfQR1RmOMFpHyZPMZcJWINAbwqzqrknuSeAFYoaqPVDFfF1Wdrap/BPI4slv+inwPjBCRru77G8uR4/S8ATwKrFPVHL942gLZgW6DCS4ReQOnNNpDRHJE5Gqce4RXi8giYBnOKKeBLu9rYDJwqru8Ue7rvxORHJxhHhaLSFWjeXrOSjQm7LkH70lAK/fgukdVX/CZngIcrKQaC5zkcjzwhapOFZGBQJaIHMLpzfiuAEMZgdPUdIk4o3AC3KWqn/jN97CIdMMpqXwJVHn1qqq5InIF8IaIxLsv3w2sdh9PBh7HGf7X1xDge1UNm6GkI11FpVJXrZo8q+rISl5/AniiNsv0gvXebOo9EZkApKnq3yqZPgi4SVUr+j1CvSUijwMfuk1gjQlbVqIx9Z6q/mQURr/pC9wfWkbXwW9pwslSSzKmPrASjTHGmKCyxgDGGGOCyhKNMcaYoLJEY4wxJqgs0RhjjAkqSzTGGGOC6v8BduQoqevidJAAAAAASUVORK5CYII=\n", "text/plain": [ "