\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.11?urlpath=lab/tree/mcmc_sampling.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",
"[mcmc_sampling.ipynb](../_static/notebooks/mcmc_sampling.ipynb) |\n",
"[mcmc_sampling.py](../_static/notebooks/mcmc_sampling.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting and error estimation with MCMC\n",
"\n",
"## Introduction\n",
"\n",
"The goal of Markov Chain Monte Carlo (MCMC) algorithms is to approximate the posterior distribution of your model parameters by random sampling in a probabilistic space. For most readers this sentence was probably not very helpful so here we'll start straight with and example but you should read the more detailed mathematical approaches of the method [here](https://www.pas.rochester.edu/~sybenzvi/courses/phy403/2015s/p403_17_mcmc.pdf) and [here](https://github.com/jakevdp/BayesianAstronomy/blob/master/03-Bayesian-Modeling-With-MCMC.ipynb).\n",
"\n",
"### How does it work ?\n",
"\n",
"The idea is that we use a number of walkers that will sample the posterior distribution (i.e. sample the Likelihood profile).\n",
"\n",
"The goal is to produce a \"chain\", i.e. a list of $\\theta$ values, where each $\\theta$ is a vector of parameters for your model. \n",
"If you start far away from the truth value, the chain will take some time to converge until it reaches a stationary state. Once it has reached this stage, each successive elements of the chain are samples of the target posterior distribution. \n",
"This means that, once we have obtained the chain of samples, we have everything we need. We can compute the distribution of each parameter by simply approximating it with the histogram of the samples projected into the parameter space. This will provide the errors and correlations between parameters.\n",
"\n",
"\n",
"Now let's try to put a picture on the ideas described above. With this notebook, we have simulated and carried out a MCMC analysis for a source with the following parameters: \n",
"$Index=2.0$, $Norm=5\\times10^{-12}$ cm$^{-2}$ s$^{-1}$ TeV$^{-1}$, $Lambda =(1/Ecut) = 0.02$ TeV$^{-1}$ (50 TeV) for 20 hours.\n",
"\n",
"The results that you can get from a MCMC analysis will look like this :\n",
"\n",
"\n",
"\n",
"On the first two top panels, we show the pseudo-random walk of one walker from an offset starting value to see it evolve to a better solution.\n",
"In the bottom right panel, we show the trace of each 16 walkers for 500 runs (the chain described previsouly). For the first 100 runs, the parameter evolve towards a solution (can be viewed as a fitting step). Then they explore the local minimum for 400 runs which will be used to estimate the parameters correlations and errors.\n",
"The choice of the Nburn value (when walkers have reached a stationary stage) can be done by eye but you can also look at the autocorrelation time.\n",
"\n",
"### Why should I use it ?\n",
"\n",
"When it comes to evaluate errors and investigate parameter correlation, one typically estimate the Likelihood in a gridded search (2D Likelihood profiles). Each point of the grid implies a new model fitting. If we use 10 steps for each parameters, we will need to carry out 100 fitting procedures. \n",
"\n",
"Now let's say that I have a model with $N$ parameters, we need to carry out that gridded analysis $N*(N-1)$ times. \n",
"So for 5 free parameters you need 20 gridded search, resulting in 2000 individual fit. \n",
"Clearly this strategy doesn't scale well to high-dimensional models.\n",
"\n",
"Just for fun: if each fit procedure takes 10s, we're talking about 5h of computing time to estimate the correlation plots. \n",
"\n",
"There are many MCMC packages in the python ecosystem but here we will focus on [emcee](http://url), a lightweight Python package. A description is provided here : [Foreman-Mackey, Hogg, Lang & Goodman (2012)](https://arxiv.org/abs/1202.3665)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"from gammapy.irf import load_cta_irfs\n",
"from gammapy.maps import WcsGeom, MapAxis\n",
"from gammapy.spectrum.models import ExponentialCutoffPowerLaw\n",
"from gammapy.image.models import SkyGaussian\n",
"from gammapy.cube.models import SkyModel\n",
"from gammapy.cube.simulate import simulate_dataset\n",
"from gammapy.utils.fitting import Fit\n",
"\n",
"import emcee\n",
"import corner"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulate an observation\n",
"\n",
"Here we will start by simulating an observation using the `simulate_dataset` method."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"irfs = load_cta_irfs(\n",
" \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SkyModel\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- --------- ----- -------------- ---------- --------- ------\n",
"\t lon_0 0.000e+00 nan deg -1.800e+02 1.800e+02 False\n",
"\t lat_0 0.000e+00 nan deg -9.000e+01 9.000e+01 False\n",
"\t sigma 2.000e-01 nan deg 0.000e+00 nan False\n",
"\t index 2.000e+00 nan nan nan False\n",
"\tamplitude 3.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n",
"\t lambda_ 5.000e-02 nan TeV-1 nan nan False\n"
]
}
],
"source": [
"# Define sky model to simulate the data\n",
"spatial_model = SkyGaussian(lon_0=\"0 deg\", lat_0=\"0 deg\", sigma=\"0.2 deg\")\n",
"\n",
"spectral_model = ExponentialCutoffPowerLaw(\n",
" index=2,\n",
" amplitude=\"3e-12 cm-2 s-1 TeV-1\",\n",
" reference=\"1 TeV\",\n",
" lambda_=\"0.05 TeV-1\",\n",
")\n",
"\n",
"sky_model_simu = SkyModel(\n",
" spatial_model=spatial_model, spectral_model=spectral_model\n",
")\n",
"print(sky_model_simu)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Define map geometry\n",
"axis = MapAxis.from_edges(\n",
" np.logspace(-1, 2, 30), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
")\n",
"geom = WcsGeom.create(\n",
" skydir=(0, 0), binsz=0.05, width=(2, 2), coordsys=\"GAL\", axes=[axis]\n",
")\n",
"\n",
"# Define some observation parameters\n",
"pointing = SkyCoord(0 * u.deg, 0 * u.deg, frame=\"galactic\")\n",
"\n",
"\n",
"dataset = simulate_dataset(\n",
" sky_model_simu, geom, pointing, irfs, livetime=20 * u.h, random_state=42\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVEAAAEMCAYAAAB9ZoVrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXu8VWW1978/uYiiiIIXBGxjoh4kUzJE8xgqJfIa9JqZeUKxzGNqN+t47WIn63jMzMyULA21TCnJyBCvUa8pKKKCF1RSUryGhCIoCoz3jzmXLNYcz2btvde+rL3H9/NZH9Yaa8y5nnnh2c98xvMbQ2ZGEARB0Dw2ae8GBEEQ1DPRiQZBELSA6ESDIAhaQHSiQRAELSA60SAIghYQnWgQBEELaNNOVNJVkl6R9EiZbUdJd0n6g6Qtctumkm6QtEjSHEkNZf5n5fYnJB1aZl/chocSBEEAtP1IdAowtsL2JeCLwC+Az+S2zwH/MrNdgB8B/wsgaRhwNLBHvp/LJHVrSgMkXdvcxndk4rjqj856bJ31uFK0aSdqZn8FllWYuwHr8pdy2wTg6vz974BDJCm3X29mq83sGWARMDL3+2eVzTiimc3v6MRx1R+d9dg663G5dG/vBgCXAtcCrwHH5LaBwHMAZrZG0mtAv9w+u2zbJbkNM/tgWzU4CIKgRLt3omb2D+DACrM810bsjZI/XpT+Om6+Y79+G2yzNrWdY9u6h++7iTOmX716Yy1bz5sJ+6aOrVfPoq3v5j15/3v6Fc7Fm6uc7Xv5v7XS8U2dG28OJXUhvMedNQnfdyo+b9azJ0P6FY8r9QjVyzlhbyaug3ce163zfbs5B/xOZWNzulf5v6p3z540OMfmnVvvuADWOu3dxLtxgZVvF21bJu6FTZxGePcSwDsVRyAJSSvLTNPMbKK/tc/YsWNt6dKlVfk+8MADt5pZ5TRhm9HunWiCJcBgYImk7sBWZNMAJXuJQcALG9tZfgEnAvSWbOCyDWcUKucXSnj317d293379Cna7rvP9+3mnPVbE72KN7w+aETRtvvB4B3JzJlF3wOc7QGmTS/aXvVd2cqxJfoftnSO94XE8S5wbCOWFY/LOd0AHLh/0XbHPb6vdx6cnwJgu+2KtkWLfN9Bjq/Xsd7ykv+D3rk9cKRjBF512rtV4uRMnVe0nTDa991mm6Jtxgzfd9FbG36+f/hwHlqwoLfvXR1Lly5l7ty5VflK6t+S32opHXWJ03TguPz9kcBdlmVKmQ4cnUfvhwBDgURXFQRBfbOuylf70qYjUUm/AUYD/SUtAb5tZlc6rlcC10paRDa8OhrAzB6VNBV4jOyJ8BQzSz1xBkFQtxjpSZ+ORZt2omb26Sr93gI+mfjue8D3atmuIAg6Iu0/yqyGjjon2mr07w4nVMygpAIBNy4p2hYv9n379i3aUpO1eznzTZ9JTKz84qWibb+3irY5iUmNVU4wIDVXu53TrqmJOcJrnGn8dxIDhxedE/G1w33fqVOLNu8YvLlAgIceKtqc0wXA3c4cYYMznwmwww5F2wPO9gCPzC/anMvICYm5aW+udb6zT4D+zmzg9MRc7XBnkn9e4hg237xoW5g4kc6pqQFGdKJBEAQtIjrRIAiCZhIj0SAIghYSnWgQBEEzieh8EARBC4mRaIeke/diJP2tRNTxY06k9o03fN/ly4u21xNt6OnIDec7KwEAjnQiqjMeK9rOneRvP3t20TZnoe/rKWVOSqhfLnSUUE5wH4ARTvg2dc7fu0vR9oTT3oYGf/u1zuDl446KCeCmm4q2ma/4vrs6qwFWJAZKuziR7d7OCoPLE5Hxrx1QtF1wt+/7Xue+OW6M73vlHUXbUXv6vt59/stT/e7i8ks2PBEJ1WkTiTnRIAiCFhKdaBAEQTOJkWgQBEELicBSEARBM4mRaIfl1bfguopARWJu3U2rdmkirdoYJ6pyxM6+771OYKghEcDZ3Um997Qj25w1y9/eS2l2p+/KJ527IRXA6eXoGBckgjIvOb5HX+L7evnThjq2lJTTiekwa7Hv653yVHDMk0c+lfB9xQkijRtWtN3p3Afgyz7/O5Et8xdOgM8LIIE/rpuakAB7Us7TvuGPDC+rkK8qkWe16dRHJ9pRU+EFQdClKY1EW54KT9JgSX+W9LikRyV9ueL7r0uyUl5SZVySF8ScLymR5SCjy41EgyCoF2o2El0DfM3M5knaEnhA0u1m9pikwcBHgGfL/A8jewAaCuwLXJ7/6xIj0SAIOii1GYma2YtmNi9/vwJ4nLw2G1k14dPZsLrNBOAay5gN9JU0ILX/GIkGQdABaZLss7+k8loiV5jZFZ6jpAZgb2COpPHA82b2cFZM+F3eLZSZUyqI+aK3z+hEgyDogDQpOr/UzPbZmJOkLYAbga+Q9dDnAB/1XBMNculynegOW8J/Vcxu/CkRzXzBSSbsReEBHnFkgY8nEgcf5ERqn37a973IiZ6e6Gy/JCEbXbi4aLsiEen9rBPpHZNIBryjI0fdt8H3XeS04ahBvm+qUFwlhyeSOl/jJHU+NZUA+uaibUSD7+sl4/7yaN/3hllF25+dSPwof3N3NcPFzrUBX7a5446+r3c/p1ZfXOYULTwtsYLkvytWLqwZ7vs1ndpF5yX1IOtAf21m0yS9DxgClEahg4B5kkbSxIKYMScaBEEHpWbReZHVbXvczC4CMLMFZradmTWYWQNZxznCzF4iK4h5bB6lHwW8Zmbuozx0wZFoEAT1QE0X23+IrGT6AkmlAjJnm1miCDQzgHHAImAVcHxjO49ONAiCDkjt8oma2d1sJLlUPhotvTfglGr3H51oEAQdlPpQLHW5TvTNN+HRRza0NeUkDHOCOgBLnHyPwxp83+HOxHsqUDJwWtE2wAkcTEtICPdyAkC3JwJpZzm5PLdKBBO8AFCfhO8sx/aZxCAjVXm1klTF0hHO9bnSCSABDHXa6wXBAI45qmhLVRz96jFFm5e7dJ4jDwUY7gQvT0pokx9zrvsrCfntE07A6vuJwKFTvJbErc/3K+TRZ6xIODaZ6ESDIAiaSSQgCYIgaCHRiQZBEDSTGIkGQRC0kEjKHARB0ExiJNph6dcPjqmInk6Z4vtu51T7vDtRdXFXR8b48GLf9zHHfkAiY+EAJzvuH53o+mgviy7Qy4nOr0pEhV90hG3znATBAD0cW0oXd6hzly1xIsXg35DjxxdtUxxZIvjJrT+SiGz/3IlMfyIhR/WSXr+eKOe6zMkY7QXMR3gZpIE9nfZOc34f/IqyeyZWSXgJr09t8H09meuMVPnaitUutlPCr8lEJxoEQdBMYiQaBEHQQqITDYIgaAHRiQZBEDST2mnnW5su14m+9DJccNGGtuSlcoIqryVct3u7aEvlzjp2ZNHWs6fvO3du0fY+R575KUeWCHDs94u21Yl2He8Et3ZJnJz7nKBM6jxu279oSwW3nksFLyr3mbhzJ00q2m67zfd1CqkmeciJDI1NBKw86aonxVy61N9+8qyiLdVWL/9pSrrqMSfh28+xjU/k0n2mQv6azF7cJGJONAiCoIVEJxoEQdACohMNgiBoJvE4HwRB0AIisBQEQdBCYiTaIenRHQZVRIsPPtj3nelUWPSkcwB/dKKvZ4/xfR+cV7RtnpAADnRkiJ4k77rr/O09aeG9icj4HU673vFdeY8jJ90r4euxMtGGlY7NqzzpKCsBuOSSoi1VQLSbY0tJOT/urKjwZMEAbzqNm+qc24ZEuzzV5MsJ3xFOxNxZvAHAcZcUl18cP95pWIK3Eif9+CM2/PyXhFS4acTjfBAEQQuJTjQIgqAF1EcnGnXngyDogJQe52tSd36wpD9LelzSo5K+nNt/IGmhpPmSfi+pb9k2Z0laJOkJSYc2tv/oRIMg6ICUovPVvDbKGuBrZvZvwCjgFEnDgNuB4Wa2J/AkcBZA/t3RwB7AWOAySd4UOtAFH+ffWQMvVeSyvDARlPmnYzsiMWu/oxOQuDZRVXNMqmyiw+6O3m8PZ3sv3yXAp53Kk/N/4ft6sa1jE3LSKVOLthEJbeI9C4u23RNBmZVOgO6EzxRtmzqBLYD7nSqgqeqXezhVV72gHcCWTo7OOxLXd9y4os07t9slZJQjnSBWqkLrGqcPSck+51xWDCKNTFyHUaOKtiVLfN8ZMzb8vC4V2WoytXmcN7MXyVXYZrZC0uPAQDMrFwTPBo7M308Arjez1cAzkhYBI4F7vf3HSDQIgg5KbR7ny5HUAOwNzKn46rPALfn7gcBzZd8tyW0ubd6JShqbzzMsknRmbttD0r2Srpa0SW77bj5X8ZCk2yTtmNsl6ZJ8+/mSRuT2Bkmz2vp4giBoDZo0J9pf0tyy14neHiVtAdwIfMXMXi+zn0P2yP/rkinRIJc27UTzeYWfAocBw4BP5/MPpwHjgbnAR3P3H5jZnma2F3Az8K3cfhgwNH+dCFzedkcQBEHbUXUnutTM9il7XVG5J0k9yDrQX5vZtDL7ccDhwH+YWamjXAIMLtt8EOnqN20+Eh0JLDKzp83sbeB6svmHbqz/0yOA8r8UZGvcSwc4AbjGMmYDfSUNANaSXlcdBEFdUdPovIArgcfN7KIy+1jgDGC8mZXLP6YDR0vaVNIQsgGbM9ue0daBJW+uYV/gx8CfgKeAyaUvJX0POJYsjedBjexjoJnNBSq0E0EQ1C81085/CJgILJD0UG47G7gE2BS4PetnmW1mJ5nZo5KmAo/ljTjFzNamdt7Wnag712BmD5J1ppVfnAOcI+ks4FTg26l9NPqj0rXkHWwviomVnbzHAJzoVNB8PhGhvMeRxCUCn9z7WNG2fSLa7CXu3cFpV0pG+ftpRduFJ/m+l04u2m50tgfYxmnvEYk/YWOclQuPPFK0gR9tvnd20bYiIc9saCjatklEwf/kRLxTjzJ3OlLGcYlVFq85bfNWZKRko2+8kWiEg7eawMkPDvjnISU3vtaR2nqSXIAdKz7/47nnkFTuPs3MJiY2T1A72aeZ3Y3fb8xwbKVtvgd8r5r9t/XjfJPmGsq4DvhEc/dhZhPNrLeZ9U5Ukw2CoEbsNHgwpf9v+auJHWjOunXVvdqZtu5E7weGShoiqSfZgla3grikoWUfxwOl1YbTgWPzKP0o4LV8HVgQBJ0Fs7rpRNv0cd7M1kg6FbiVLJh0lZk9mnA/X9JuZGP6fwClh9AZwDiyCkirgONbt9VBELQLHaCDrIY2VyyZ2QwamYso8/tEwm7AKbVuVxAEHQizTF5YB3Q52edmPWCP7Te07ZHwffGloq1fIkgx0wksJeI3eGG+EU6lTfCDAV5QZlnifuvnXOGUtNFLA3m6I2EEmOekoUzl4rznnqLt8MN93z/MKtrGji3a7rrL395rwzxHdgowbnTRtmNllCRnmhNg2yUhb/Tyn57gHMN9iUUzw5wg1KOJ63uEk2+WRPDzgoS82cOb5xueuPfHVOTNnZ44303GYiQaBEHQPEpzonVAdKJBEHRMohMNgiBoJkZ0okEQBM0nHuc7LCvfgdkVE+9HOPkbwc9Z+dNEEa7THdskJw8m+MGiKVN830FO4MALfqxKKJa84MdDDxVtAEftWbRNdlfxwtmTirabb/Z9H1tctHnHBbCbo8bazFHVvJ2Q5YwfX7R96jQ/i9nf//p8wXb33f5+H3IChwsT52Z0g2+vpGdP3/4npw0JQRt/dtRvX0goxxY6wbFUbtt5zn49NRkUCzqu8yrtNZWIzgdBELSQGIkGQRA0E6NuljhVJfuU9B5JY/L3m0nasnWbFQRB16Z+ZJ8b7UQlfR74HfCz3DQIuKk1GxUEQdBpOlEyieWHgNcBzOwp0lnegiAIWk4nS0Cy2szezpOWIqk7G8nf2ZHpQTaULucBR8II8AcnOJgoaOlWbqysgvjuPpyd7L+/7+tJPD3Z5upEtPqvjuTy4oQe1asoOdaJ2IO/GiC1QuAbpxVtqUjvXnsVbX9zjuGLX/K3/9WvirYXXyhG4QH2c8556hiOdVZwfCsh29x6cdHWzfmfllphcL9ju8RZdQDFyrUAv0zkgPVypT79tO/r0SuxROCFimqqNekcjPRN0sGophP9i6Szgc0kfQQ4Gfhj6zYrCIKujcHa9h9lVkM1j/NnkpVgXwD8J1kGpm+0ZqOCIOjilBRLneFx3szWAT/PX0EQBG2A1c0Sp2QnKmkBjUxvmFlitiwIgqAGdIBRZjU0NhItZXwsJUC+Nv/3P8gyytcl6yg2/tnE/PVJjjRxZiJX4zgn7+aZTpADYE9H7pfK8ekVqnvakSCOcOSSAHstL9q8ABLAWuc8zJrv++7kSPtSBeG8vJkjE1JbnXRiwTb0lUIZ8WSUY8AOxZMzOCEx9Yra/S0h+/RiHJ9NrFH5wslF2zfPLdpGH+Bv/z6nDY86MkyAvzgy5GGJQmIDnGP4L0+vDEw8t2g7drTvW3l95ZWEaypm9R9YMrN/AEj6kJl9qOyrMyX9Dfjv1m5cEARdmDoZiVYTWOot6d2/mZL2B3q3XpOCIOjyWB6dr+a1ESQNlvRnSY9LelTSl3P7NpJul/RU/u/WuV2SLpG0SNJ8SYm6ExnVLHH6HHCVpK3yz8uBz1axXRAEQfOp3Uh0DfA1M5uXS9YfkHQ7MAm408zOl3Qm2UqkM4DDgKH5a1/g8vxfl2qi8w8A75fUB5CZvdbCAwqCIGicGpYHyUuqv5i/XyHpcWAgMAEYnbtdDcwi60QnANfkRTFnS+oraUCqNPtGO1FJ36r4XGpYzIkGQdB6VL/Eqb+kuWWfrzAzJxoJkhqAvYE5wPaljtHMXpRUChUOBJ4r22xJbmteJwqsLHvfiyxq/3gV23VItu4DR1TI/bzKlVBM3gzQN7Hf/3Ii8V8e7ft6Ur0XXvB9r3Ii8Rc6MzSpaLcnG01VtPSSBO+ZWMi291c+XLB9YFEiY/WwfyvavGUHAL/7XcG0/alHFv32+Ii7+SG7XFo0LvMEj/DA9KIcdPRov1n7jisuPbj8PH+/P7igaNvGiZinqqP2c2xbJSLuPXyzi7eg4T/P9X2dPNikLm9lIuyf/KUJjUrRtOj8UjPbZ2NOkrYAbgS+YmavK72MwPsiudyzmsf5H1Y05EIgkdM7CIKgRtQwOi+pB1kH+mszK2UXeLn0mC5pAFDKArAEGFy2+SAgMcypMp9oBZsDOzdjuyAIguqoYRYnZUPOK4HHzeyisq+mA8fl748D/lBmPzaP0o8CXkvNh0J1c6LlyqVuwLbAdzfa8iAIgpZQu5Hoh4CJwAJJpQpjZwPnA1MlfQ54Fvhk/t0MYBywiEybc3xjO69mTvTwsvdrgJfNrD6kBEEQ1CdGzbI4mdnd+POcAIc4/sZ6peZGqaYTPc/MJpYbJF1baasXVq4sytRSuTzvcqpXfitRSdHLDXl+ovqlE1cioZjkMkcaeJ4jCxydCI5961tF2/OJ2Z3eTjThgoR09dLhxehBUwIlfP4u3/kpZ9HHdGcK/o03/O2327Zg+vtNC1zXDxxVnJWyRX6CzQduKwaRZr7iOJI9rlXyaaeq5iuJ7f/zhKLtm7/wfb2zsH3iOhxzTNF283W+b4NjSwUkT79kw89rhvt+TcN8HXIHpJpOdI/yD3lS5g+0TnOCIAio6TrR1iYZWJJ0lqQVwJ6SXs9fK4CXWT8BGwRB0DrUST7RZCdqZv9jZlsCPzCzPvlrSzPrZ2ZntWEbgyDoitRJJ9pYPtHdzWwh8FtPgG9miVm4IAiCFlJHj/ONzYmeBpwI/ND5zoCDW6VFQRAEUDc1lhrLJ1rKjnuYmW0gPpSUqPvX8Xl1LUypCLRu/pDvO9ixXZqopHji4UXbwYmkvW85Us5liYjqjU4kfnxCAujhqCh54knfd4JTUfJn1yR0FY5Otd9JTvgXYNm/HGMiND3U0a8Oe6po22ILf/u7ilH/9x78Ht93zdqCSbvu4rq+MrOoeTwhkSDNu76eynWTxGTaP51T4yhyAX+1yNTEPXqVE4lP3HZurss7Esm8K2Potan2WT9JmatRLDkFa11bEARB7bB11b3amcbmRHcgy1yymaS9Wb9YtQ9+foIgCILa0EnmRA8lS1o6CCjXm64gk0wFQRC0HvXeiZrZ1cDVkj5hZje2YZuCIOjqdJKRKABmdqOk/0OmXOpVZq/LpMxbkWUWKOfuRAVPL4X/+xP7nTu3aLsvETvxTvoBDb6vVwX0wAOLtt8lJKbjneDH8IQsb5VTw/Wdhb4M0guebLno7/6OD/DKWnriSIAvFU2HjS7aLBE92WFA0TbzFt93r72KNicwBTDMkW2m8msOcqqLHumkRH3aP7Wc70htU/Nnpzmn4YTdfd/+/Yu2WYnqpkeMKdpuTQSWjqqIBT78pu/XZDpLJyppMtk1PAj4BXAk4BTBDYIgqBFm8E7nic7vb2bHAv8ys+8A++Gv/gmCIKgNRv0rlsooDc5XSdoReBUY0npNCoIgsA6xfKkaqulEb5bUF/gBMI/sb8TPW7VVQRAEHWCUWQ3VBJZKWexvlHQzWXApMXUdBEFQA0qP83VANSPRdzGz1cBqSb8FdmqdJrUu3TaBbSoUg30S89cjnMj2CYlo5kedTMsHOFFagMXOaoDHFvu+8x3b5k6odoQTPQbo65QnTUWV9z25eMDvzPbzzGw5wpFH3vM3f8cHXVi0vXOT79tj06LNnKTMSmTS7vdA0fYf5/i+33GqPuw40HX1KmWmknl7asVfO5LLFQnN5emOevZHieTJR+5QtE1e6Ps6wXmuSSh1v+/83ta+K9tUZBRXsqRbU+hcsk+PZK3RWiJprKQnJC2SdGZuGyJpjqSnJN0gqWduP1fSpLZoVxAErYxZloCkmlc709xOtCY5BhpDUjfgp8BhwDDg05KGAf8L/MjMhgL/Aj7X2m0JgqAdqPfovKQ/4neWIlE2p8aMBBaZ2dN5e64HJpCl4Cs9hFwNnAtcTlZuplbLfIMgaE86yZyoM5FV1Xe1YiDwXNnnJcC+wPKyaqNLcj/MrC3aFARBm1C7JU6SriKrWvyKmQ3PbXsBk8kC5WuAk83svrxG/Y/JhI2rgEkbS0DfmHa+WM6xbfHmXT2t4EanFiRdCxwB2fzF5IoJ/fMSQZmpThDpgoTvfk6Q4YnEBL9Hag69jxOwGj61aLswkWPUUzE+9pjvu++kou6zx6BEUtTKaALAxz/u+7K8aFr6T9/1rUuKtiGHOo6pmaj9CpaXv3CQ67n9pLFF43y/Muj2g5z/Kt39/z4/+n5RE+tJQXskSrw+5OS3TT1ieYFDt5wscJRz76ZymjqqXgYn7rGZMzf8/MxWzyFpZZlpWrOqA9duJDoFuBS4psx2AfAdM7tF0rj882iy6cOh+WtfsqfcfRvbeZOi823MEjZURg0CngX6Suqej0YHARuNBeYXcCJAL6nV53ODoCuz0+DBLFu+3MvrXD0Ga9fU5r+qmf1VUkPxFyj9WdiK9f3IBOCavPb8bEl9JQ0wsxdT++/Inej9wFBJQ4DngaPJ5kJHkOn3rweOIyqPBkGnwwzWFQsPpOgvqTwF0BVmdsVGtvkKcKukC8kea0rPkt404kCg/jpRM1sj6VTgVrLH+KvM7FFJZwDXSzoPeBC4sj3bGQRB67Cu+oHoUjPbp4m7/wLw1TxL3VFk/cgY/GnERltSTRan24FPmtny/PPWwPVm5k1U1RQzmwHMqLA9TRa5D4Kgk9IGwfnjgC/n739LlqEO/GnERqcMqxmJ9i91oABm9i9JiWhDx2croDKcsFNCe/WcE4B5KTFpP2VK0fb3RLDoJCeecenMog3gtD2LtmOcNpxwgr/9qd8v2ryiegAvzypGwhKxE/p5gaV5iSDmvztZWK+5pmgDOOOPTsOcnODbbe9vv2plwbT95X5JsFe/WIwG9jviw67vP+55vmB76SX/Anv5Wh95pGjzTmG236KtCbUJGZ2we53SjBlFG8CIJhQAOruizsUf/1T9tkma9jjfHF4APgzMIls2WaqGOB04NV9SuS/wWmPzoVBdJ7pO0k5m9iyApPfQBovtgyDo2tRqJCrpN2R/W/pLWgJ8G/g88GNJ3ckWI5SqG88gW960iGyJk6MN3pBqOtFzgLsllZY8HVj2g0EQBDXHzK1o3cx92acTX33A8TXglKbsv5osTjMljQBGkU26ftXMnCraQRAEtaGOBEuNyj53N7OFeQcK6ydXd8of7xtdxR8EQdBsWn9OtGY0NhI9jeyx/YfOd0Y2GRsEQVBzjCYtcWpXGpN9luY9DzOzDVRgkpzsivXBGmBZhc2rqAl+lpVnKjfO+fxRRdv3HXkmQM+eRdtXEhFzL6q70KkierIThQe/imifRKjXi8Q/npCuDlhe/OK9pzo/BnDLT4q2M27wfR3ZJtuvLtpeTpQ39Z4BN/Fzl/Yb41T7nP+w63utU4Fzv1F+E7xFCi84OspJk/ztvZUewzx5J3Cdc33OHe/7/txJy/q+xDqbJ53Kr/s4uUsBrrpqw89rvcSlTaV+KiZXlQrPWx/irxkJgiCoAWaw5p3qXu1NY3OiO5DJnTaTtDfrV/L3IV0GOwiCoCbU/eM8cCgwiWzF/g9Z34m+Dpyd2CYIgqDFdIrovJldDVwt6RNm5khGgiAIWok6mhOtZrH9ByTdWaGd/5qZfaN1m9Y6bNEdDqiY+L4rETwZ6cjytk5I9T7rBJEcdScAk50J/s8mnL2gl5dUtUfit97r1JNLqTP3dCSmPRJ3iCsHTVXA67VZ0fbzk31frzrgKy8XbUsTS5UnFtdJ/+vLX3NdlzgFA983xo+0nOI090Qn9SnAFo6twbHdcYe/vXcKUqfWI+X7Eef6ft2rhMj6lEblLElIniuFtlaLwBL1s8SpmsDSYZXaeTJZVBAEQatQWuJUzau9qWYk2k3Spnm5ZCRtBjh1bYMgCGpDKTpfD1TTif4KuFPSL8n+QHyWrEBcEARB69CZ5kTN7AJJC4BDyCL03zWzW1u9ZUEQdGk6wqN6NVSV2d7MbgFuaeW2BEEQAHl5kM4yEpU0CvgJ8G9AT7Lg8Eoza0qe2A5D796wT0UhgZ0cGSX4FRpTkc/TnEqKyxIS0becKOdViaTMXrWvvzu20YmrcacTAe4b800EAAAcbUlEQVSVEO1u7kgoUkmZ9x1d3MnfZ/on570ff1/RuKdjA3jLqWvpZTlO5Elb8bViJD61ouJup5rrIE9Tm/A9vMHf71hnpcVFk4u2118v2gBWv120TXVWEgDs6NhWOZJN8E/jfQf6vhMvLdqOSdSTqLxHZq3w/ZpKp+lEyUqNHk2WQn8f4FjAWTgTBEFQG5pYqK5dqfZxfpGkbma2FvilpNDOB0HQqtQqKXNrU00nukpST+AhSReQlQ5tWU3pIAiCRqinkWg1i+0nks2DnkomThgMfKI1GxUEQbBuXXWv9qaaJU7/yN++CXyndZvT+jz/GnyzIhXlQYmcVGucYo7bJvIvPurk/dzRm/VP2N9y8k2CH3x4v3PVkkEsx5ZKBjvSCRzcd5/v+5eZxT1v4ekdgXfmLSi2656iLfV7+zsaxFTwxDsPl13m+w5zgoEpKaYXlLk4kS/2ZieIdKJzDC8kCvF6eUp7J+7R3Xcv2n6SaJcn673JT7Xq3iPPPuv7Lq4IlK5zzlVT6RRJmfO1ocnDMDPnkgRBENSATrLEKZFrPQiCoHWpJ9lnck7UzP7R2KstGxkEQdejVglIJF0l6RVJj1TYvyjpCUmP5kHzkv0sSYvy7w7d2P43GliSNErS/ZLekPS2pLWSEsuEgyAIWk4pKXONAktTqMhMKekgYAKwp5ntAVyY24eRrYvfI9/mMkle9sl3qSY6fynwaeApYDPgBDIFUxAEQeuQL3Gq5rXRXZn9lWJ9yi8A55ey05lZSao2AbjezFab2TPAIiCh1crocovtu1Os4pmSNt7iJC/+bWK/X3b28VIiia0nJ52VkPUd6mjD7nTUlUclLnMfJwL8ZOK3bnci0wMSFR7vmFW0HeNUPAW4+OKibZuEFPNzJxVP5F/uKC6T+PAxA93t5015vmDrn0gS/LCTkHi3XX3fLR1Z7c/O9X0nOvY33nB+f7G//XaJFSAe050E3/+WOLeeZDm1ysG77LskdIqVKwRuTawUaSpNCCz1lzS37PMVZnbFRrbZFfh3Sd8jW8TydTO7n6yu3OwyvyW5LUkstg+CoMNhTUu4vNTM9tm42wZ0B7YGRgEfBKZK2pn1teQ2aE5jO6p2sf0mxGL7IAjaCCNbp13Nq5ksAaZZxn3AOqB/bh9c5jcISKzozWjKYvu36ASL7YMgqANaX/Z5E3AwMEvSrmQZ6pYC04HrJF1EliRrKJCQnGQ0tth+AjDIzH6af54DbJt/fbqZ/a6lRxEEQZCiVoolSb8BRpPNnS4Bvg1cBVyVL3t6GzjOzAx4VNJU4DFgDXBKHgtK7z/bzv3hvwFHm9lz+eeHyLLb9wZ+aWaH1OD42pxtJft4hS01p+HFhRJz9pzk5JD0ZJQAZ/530bZzIoAzw2mEo/RL4qkFvWqS4MsrL3bySgKcelLRlqoiuqsTrHnsMd/Xy3Xqbf+e/f25/pfnFQNLKUnt207eztQxeMGplNR2nnNsrzl+70sEkG53Upqm1tgc6NyQf0606yPO770voTu82gkypp6cx1fs47s2nHnzF3hzi1Wzcy/Z9xqq8z3mCR5oxpxozWjscb5nqQPNudvMXgVelRSBpSAIWpV6yeLUWCe6dfkHMzu17OO2BEEQtBJm9ZNPtLHo/BxJn680SvpPNjLRGgRB0BJqrFhqVRobiX4VuEnSMUBppugDZDXnK6cVgyAIakdnyOKUy6D2l3QwmY4U4E9mdlebtCwIgi5NZ5gTBSDvNDtNx7ldX/jSwRvaUlFWj/Nm+faHHqp+H1s6Z/0PCYno/o7vWidMmpJRerK+hQur9900cYd41SsPSET9+/Yt2lIJiXfaqWh78smibdGiYhQ+tX2qquY9jng5JYNsaCjaUisMHnZsPRzbQQk56mlOxPyiRLJor0JrQ+J+PsqR5V6YWH1x/qlFWypBd+W5UeL+agqdqmRyEARBe1D3me2DIAjai3pKyhydaBAEHZJ4nA+CIGgmnaJQXWfl2eVw8rQNbYk0ifR0bAl1JgcfXLQ94lQABXjYCQx9N1HRyssB6VVdTOVEfcf5rX6JIJQXcEqlr/m4I2mdlQg8eNVNU8Gens5J94I9qQDQ3XcXbV7+VoBnnADMiMTN4MlB+zg5RsEP7Gzl+KXOwZ+dgNUeRRMAi53csNsnyrke7gSRxiT2+03Hd0XC95KK6qR6IuHYFOoosFRNKryaoIxL8tol8yWNKPvuq5LmSfpU/rmXpPskPZzXP/lOme8QSXMkPSXphjzXKZLOlTSprY4nCILWw6hdZvvWps06UeAwsrRSQ4ETgcsBJG1BlhR1JHBM7rsaONjM3g/sBYyVVPp797/Aj8xsKPAv4HNtdgRBELQNVj+KpbbsRCcA1+RJUGcDfSUNYH0m6XdnQHKfUkGFHvnLJIksB2ApDd/VrFdPvQG82crHEARBG2DW6kmZa0ZbdqIDgfKsUEuAgWa2AlgAzAVuKH0pqVuefu8V4HYzm0NWHmm5ma0p3weAmV1oZjcQBEGnIEaiRZK1S8zsf8xsbzO77t0vzNaa2V5k6flHShre2D4a/WHpWkkrJa1MpJYMgqBGPPPsc5T+v+Wva5u6j86SgKTFSDoFKGWCup8m1i4BMLPlkmaR1YD+Idk0QPd8NFrtPiaS1YpigGQHVHy/SyJ660kpvUS+AGdcV7Ttm4qSNhRtqYqU1cpJFzqJfAE+6YRfb01ICIc72Z63S0gIvRUCH608sY2wLBFdP9eJro92/LZJRMYPcNqQWrmwm3O8nqQW4DknCv4nZ+UE+AmYPWYnqq56p/FwJ+k3wHUzi7ZtvGUlwHBnBJFKx+YlDL7yS75vZcXRnQYN5tVly1ucc7gD9I9V0aojUTP7qZntlY8obwKOzaP0o4DXzOxFbztJ20rqm7/fjGwlxsI8ff+fgSNz1+OAP7TmMQRB0D6sq/LV3rTlOtEZwDhgEbAKOL4R3wHA1ZK6kXX0U83s5vy7M4DrJZ0HPAhc2XpNDoKgPTDS5Ug6Gm3WieajyFOq9J0P7J347mmy5VBBEHRSjI4xyqyGLqdYCoKgPohOtIPSZwsYW1EXMJVf8/b5RdvQREDjm+OLtsnTizaAlxcXbcft4l+KRVOKDzVjRxf97p7l/9Z1ThDpw00oF+pVngQ4x5Gp3uIEOQC2dYJmiZgbk5wgn5fvdXFCMvmcE5g6IqFt9IJIF83yfU9wcqUO8F051nlOesCRjQ5NSExvcu5HTw4LMMjJJzorcW7GOZrl5Yk8tu93bsfKAFKJhYs3/GzDfb+mUqtOVNJVwOHAK2Ybtk7S14EfANua2dJ8LfqPyaYeVwGTzCxRAzajLZc4BUEQVEXpcb5GgaUpZKt7NkDSYOAjQPlaE1dZ2RjRiQZB0CGpVSdqZn8FvMV6PwJOZ8O15illZZIu9zgfBEHHp4nR+f6S5pZ9vsLMrmhsA0njgefN7OHsCf5dXGUl4C7HhOhEgyDooDRhTnSpme2zcbcMSZsD5wAf9b52bI2qIrtcJ/rqGzBl1oa2rRO+n3OCJ3clSvZ915l075fY70FOkOLqX/h/d70bySuW977Eb23q2C5NFapzbCcn1FwX31y0nX6E7/uwE6C7xndlZ0fFMy4RzPPwgki/TSi0PugEdhKpVhk9umhLBSQ91ZSXW/aVRNDuGOf+WJJQN3myxwbflblOEGlxwne4czu+NxEIuye1kxbQykuc3gsMAUqj0EHAPEkjyUaeTVJWxpxoEAQdktZSLJnZAjPbzswazKyBrOMcYWYvAdOpUllZIjrRIAg6HLWMzkv6DXAvsJukJZIay0E8A3iaTFn5c+Dkje2/yz3OB0FQH9Tqcd7MPr2R7xvK3letrCwRnWgQBB2O0M4HQRC0kJB9dlDeAiqDqmc70VCA+5xki0sSeTB/79iSay4cEdnzKV+HUY4tEShmqWNLVSxd7m3v7QC4+NSiLSX7/KcThT4q0YYxTmTbq3i6beI6rHZyZo5NpKvxpJjDHBklwHkXFW0JV2YmzkMl0xK5Wnd17HOLJgActXGyXaP3LNo+nDiPzzurAbxzC+kVDS0hEpAEQRC0kOhEgyAIWkB0okEQBM0kAktBEAQtIOZEOzDbbQqnDN7QdmkiW+DpTpAjVVBu1WNFW1JC6AQ6Nk9EA7o5V2jx4qJtx0QOybOdoM7kRD7R3znRKU/uCPD5S4u2wUUTALs4J2KkE+QAON3JB/oZR/b5p8Qw5QNOIcGLE9XY9nBspxzt+947u2jzrgP4wTgvIDk2kVR1rhPAmZKQ1J4xrWj7riNXBnjBES/+IVFsbwvHNse5NpCWN7eU6ESDIAiaSYxEgyAIWkh0okEQBC0gOtEgCIJmEtH5IAiCFhBzoh2YHj1gQIXucUQiQjnfSSa8LBEFP2r/om2TRKLBWU6U0wkqA7DCsT3l2M4rlOHKuNiJuKeiyl6Fx6cS5+ZTToT/noT29GJHxnhgQvJ4oqNJneckE77M0zviJ8feyXdlK8c2L7FS4+/OeRg2zPft7pxHb1XHXGdFB8Bax+ZVPAWYcX4xxD/xTF+f2eDYJiQSLXtS3dS9X6mI3dl3azLRiQZBELSA6ESDIAiaSTzOB0EQtJDoRIMgCJpJROc7MGbFConbJCSXrzgT6QcmclP+7p6ibUQicefwhOzSw5ODejkg30rkevQqSi5L+H7jtKLtS04eTfArpA5KnMf3O+3dzHd1q1futV3R5gWQwJdybpq4y19w/pdul8if6l2H2xKBodOPKdouv65o6+FvzoSGou3cWb7vuPuKF3N/53wB3OTcC1skqoiOcpLWevl1AS6qkOVO9t2aTIxEgyAImknMiQZBELSQeulEo2RyEAQdjhqXTL5K0iuSHimz/UDSQknzJf1eUt+y786StEjSE5IO3dj+oxMNgqBDsqbKVxVMASrlKLcDw81sT+BJ4CwAScOAo8mm18cCl0nq1tjOoxMNgqDDUcuRqJn9FVhWYbvNzEp98GxgUP5+AnC9ma02s2eARUAinJzR5eZEV60qSvsaGnzfBudPzBtv+L4Xn160jb/A9z3TkdqlpJgeMx155ciELHCEU8n0vERy3U85kfiffMb3Xe6UBp2bKEnZ24nOfyxxW85yIsAHOMfQz4k0A3jKxNEJeeZQR2u7T6JEqyflnDHD9/Wqkw52EksvTMgoPcYlVj70chI7e5F1gBecFQ3LEys1Fjr32DOJyqBU2C2RuLyptOGc6GeBG/L3A8k61RJLcluSLteJBkHQ8WlidL6/pPI/4VeY2RXVbCjpHLJZgV+XTInmJIlONAiCDkkTOtGlZpZ4hkgj6TjgcOAQMyt1lEvYsNLNIMAprLKemBMNgqBDUqs5UQ9JY4EzgPFmVj4hMR04WtKmkoYAQ4GEzCAjRqJBEHQ4ain7lPQbYDTZY/8S4Ntk0fhNgdslAcw2s5PM7FFJU4HH8iacYmZedsL1+18/iu0aNPSUfXP7DW2rUhPmDtMTAZzxTkXLRQnfd6r/Ob59UtE22dHVOapTAHo7ttR6jbGDirY7ErJAjwMTckMvr+rDTo5QgI84VUBvd/K67p/Ig+kF6JYl/jcucGxjEhU4vYDTjYkA3T8dm5e3c2UiSOkF7fr2LdoA3naCYz17+r5jnZyzl03xfb1g3sPOdQD4fxXn96Xhw3l4wQJvbrFqekhVx6deggea8zhfK2IkGgRBh6ReFEvRiQZB0CGJTjQIgqCZRAKSIAiCFhKdaBAEQTOJpMwdGLNiRNOT9AFc40gLfzDG9/3jHUWbV5UT4P2OLRHI56abijavqGZKQXiEkxh6G2clAcB+TsXSVdN8312caPMfE6vphjiSxdQKhWudCPDgoonnE6sGeju/5a1QANjN8Z2SWDXwhhOJ944LYHtntUdlhVmA2xLRfU+JOSqRyNuTmG6xhe/7nuFbFmzz3Hqy0OCc336J+2Zkxf+TP/puTaZeRqJttthe0u6S7pW0WtLXK747WtI8SV8ps31A0oI8JdUlyhdzSdpG0u2Snsr/3Tq3T5J0blsdTxAErUctE5C0Nm2pWFoGfAm40PnuaOCDwChJpb+jlwMnkikGhrI+ldWZwJ1mNhS4M/8cBEEnIzrRCszsFTO7H/9JrrQw1wBJGgD0MbN7c03rNcDHc58JwNX5+6vL7G8CieXLQRDUE/U0Eu0oc6LTgLnAr8xshaTdyBIBlChPR7W9mb0IYGYvStouf38DVdB9s83pv2d1leJ2cyYqew/1fQc6vm8m9rujY0tMY9G/X9E2xNnBpk3Yvs9Wvm+PnYq2bb0JXKCPo24anIgEbO+ogFKF9bymbe/Ytk0ccI8m3NFeGrndvIuDf816J9RN3rFtuWvRNjChlFvt2LZI3Hf9nIp/m6eqAG5XnB0eOmKl69o/Mf/pUnHv652maPLS1Etgqc1ln/m85Rtm5j3Wl3w+CPyPmY3JP/87cLqZfUzScjMrT+X/LzPzik+W7+9a4AiATTbZZPO99967BkfSsXjmmWcYMmRIezej5nTW44LOe2wPPvgg69atK/8TMc3MJjZlH5JmAtUqP5eamSNqbRtadSQq6RTg8/nHcWbWaEqpMpawPtM0bJiO6mVJA/JR6AAgkZ53PfkFnJi3aeXcuXNTAdu6RdLKV199NY6rjuisxyZppZm16Ljas1NsKq06J2pmPzWzvfJXtR0o+eP6Ckmj8qj8scAf8q+nA8fl748rswdBELQ5bfY4L2kHsnnPPmTzwW8Aw8zMXeIoaR+yAlObAbcAXzQzk9QPmArsBDwLfNLMUsssvf22+K9kRySOq/7orMfWWY8rRZsFlszsJTZ8RN+Y/1xguGN/FTikBU1JLB+ve+K46o/Oemyd9bhculw+0SAIgloS5UGCIAhaQHSiQRAELaCuO1FJV0l6RdIjZbYdJd0l6Q8lCWledOqGXIc/R1JDmf9Zuf0JSYeW2Re34aEUkDQ2b9MiSWfmtj3y/ANXS9okt31X0nxJD0m6TdKOuV15zoFF+fcjcnuDpFntdmAJEsc7JL9eT+XXr2duP1fSpHZqZy9J90l6WNKjkr6T26dIeia/Dg9J2iu3by3p9/k1uE/S8LJ9eTkjZpbte7Kkbrm9w+WMqJdr1uqYWd2+gAOBEcAjZbbzgT2AjwEn5baTgcn5+6OBG/L3w4CHyQQ/Q4C/A93y7xa343F1y9uyM9Azb+Mw4EpgW+CLwNjct0/Zdl8qO85xZKsaBIwC5uT2BmBWe1+7Ko93KnB07jMZ+EL+/lxgUju1VcAW+fsewJz8/E4BjnT8fwB8O3+/O1neh9J3N+XHfn3ZPvuU/c6NZcd/AXBm/v5M4H/z95OAc+Oatd+rrkeiZvZXilnkurFeVlvS5Jfr7X8HHJKvP50AXG9mq83sGWARMDL38+qNtRUjgUVm9rSZvU32n2wC2bGVZMUCsA2XiPXOvyf3v8YyZgN9c3HCWtKZ99qL1PEeTHa9YMM8CW+QVtW2Kvn5LOVo6JG/GovODiNLlIOZLQQaJJWUrBvkjMh9StezO1nnVH49O1LOiLq5Zq1NXXeiCS4FfgacBPwqtw0EngMwszXAa0C/cnvOuxp9M/tgG7XXI9WuHwN/AvYDbit9Kel7kp4D/gP4VmP7MLPnzOyIVmx7c0gd7/L8epXbMLMLrcpcCa2BpG6SHiJTy91uZnPyr76XP7b/SFJJ3f8wueRY0kjgPaxf6lfKGTHXzFaU7f/WfN8rWN8hbZAzAng3Z4Q1IqFuRerqmrUmna4TNbN/mNmBZvaxshvTK99qjdjbG7ddZvagme1rZp+xslrYZnaOmQ0Gfg2c2tg+WqGttcBrq1fZuUO038zWmtleZJ3hyHye8yyyx/UPAtsAZ+Tu5wNb553uF4EHyXNrmNnVZra3mf2wYv+HAgPIppkOboNDag51dc1ak07XiSZYQp4gXVJ3smRBy8rtOeUa/fakue26DvhEC/fRHnhtfZZsCqJ7ma1Dtd/MlgOzyOanX8wf9VcDvySfFjKz183s+LzTPZZsTvuZKvb9FpnEeUJuejmfjkFV5oxoZerymrUGXaUTLdfbHwncZdls93Tg6Dx6P4Qs+XOiyEWbcj8wNI909iQLhk33HCWVJ0kbz/rqIdOBY/Mo/SjgtdLjYAckdbx/Jrte0EHyJEjaVlLf/P1mwBhgYVkHJ7J5wEfyz31LEWrgBOCvlpY6b1G2n+5kwcHy69mRckbUzTVrddo7stWSF/Ab4EWyRM9LgM8l/HoBvyULHN0H7Fz23TlkUcYngMPa+5jK2jUOeDJv2zmN+N1I9h92Pll5m4G5XcBP8+0XAPu09zE19XjJIr/35dftt8CmHaCde5I9ks/Pz/u3cvtd+Xl+hGwuvhRt34+s3NZCsjnQrRvZ9/ZkndN84FHgJ0D3/Lt+ZAGqp/J/t+kA56Iurllrv0L2GQRB0AK6yuN8EARBqxCdaBAEQQuITjQIgqAFRCcaBEHQAqITDYIgaAHRiQZBELSA6ES7EJK2l3SdpKclPZCn1fu/G9mmQWWpBpv4e5NKqfnyz7+QNKzKbUdLurk5v1stku7J/22QdEwztp8k6dLatyyoJ6IT7SLkSpqbyBQzO5vZB8hUJlXXvWoGk4B3O1EzO8HMHmvF32sSZrZ//rYBaHInGgQQnWhX4mDgbTObXDJYlqzlJ/DuaOz/5UmC50nav3IHjflIOl3Sgjyh8PmSjgT2AX6dJyneTNIsZVVcSwl95+X+d1Z7EJIOkfRg/ltXlbIlSVos6Tv5PhdI2j23b5snMZ4n6WeS/iGpf/5dKYXc+cC/5+38auUIU9LNkkbn74+X9KSkvwAfKvPZVtKNku7PX+9+F3Ry2lsyFa+2eZElbP5RI99vDvTK3w8lS88G2SjtkY34HAbcA2yef94m/3cWZXLT0meyJBzPAUPK/SvaMxq4ucLWK99u1/zzNcBX8veLycpqQ5aE+xf5+0uBs/L3Y8myCvXPP7/h/RbZCPrSss835z4DyJJsbEuW6/NvJT+y5C8H5O93Ah5v72ser7Z5tVnJ5KBjIemnwAFko9MPkiUXvlRZWYu1wK7OZimfMcAvzWwVgJltLOnzKLJphWeq9C+xG/CMmT2Zf74aOAW4OP9cKtX7AHkOz/wY/2/+OzMl/avK3/LYl6wqwD8BJN3AhudgWDZrAkAfSVtaWZ7QoHMSnWjX4VHWp8nDzE7JH2vn5qavAi8D7yeb5nnL2UfKRzQtb2RT/cu3a4zV+b9rWX9vb2wbjzVsONXVq+x9qt2bAPuZWafM3h6kiTnRrsNdQC9JXyizbV72fivgRTNbB0zET7Cb8rkN+KykzSErqpbbVwBbOvu5F/hwnn6w3H9jlMpr7JJ/ngj8ZSPb3A0clf/OR4GtHZ/Kdi4G9pK0iaTBrC8ZMwcYLamfpB7AJ8u2uY31CbHJR+tBFyA60S6CmRlZnssPK6tKeR/Z43ApA/tlwHGSZpM9oq50duP6mNlMslySc5VlcP967j8FmFwKLJW15Z/AicA0SQ8DqbIRh0haUnoBewPHA7+VtICs1tTkxLYlvgN8VNI8srnbF8k6zXLmA2vyINdXyeY6nyFLbXchMC9v94tkBdfuBe4o2XO+BOyjrDzIY2TlaYIuQKTCCzo1efR+rZmtkbQfcLllWeaDoCbEnGjQ2dkJmCppE+Bt4PPt3J6gkxEj0SAIghYQc6JBEAQtIDrRIAiCFhCdaBAEQQuITjQIgqAFRCcaBEHQAv4/bdDVuWglOpQAAAAASUVORK5CYII=\n",
"text/plain": [
"