\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/spectrum_models.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",
"[spectrum_models.ipynb](../_static/notebooks/spectrum_models.ipynb) |\n",
"[spectrum_models.py](../_static/notebooks/spectrum_models.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spectral models in Gammapy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"This notebook explains how to use the functions and classes in [gammapy.spectrum.models](https://docs.gammapy.org/dev/spectrum/#module-gammapy.spectrum.models) in order to work with spectral models.\n",
"\n",
"The following clases will be used:\n",
"\n",
"* [gammapy.spectrum.models.PowerLaw](..\/api/gammapy.spectrum.models.PowerLaw.rst)\n",
"* [gammapy.utils.fitting.Parameter](..\/api/gammapy.utils.fitting.Parameter.rst)\n",
"* [gammapy.utils.fitting.Parameters](..\/api/gammapy.utils.fitting.Parameters.rst)\n",
"* [gammapy.spectrum.models.SpectralModel](..\/api/gammapy.spectrum.models.SpectralModel.rst)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"Same procedure as in every script ..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import astropy.units as u\n",
"from gammapy.spectrum import models\n",
"from gammapy.utils.fitting import Parameter, Parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create a model\n",
"\n",
"To create a spectral model, instantiate an object of the spectral model class you're interested in."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PowerLaw\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- --------- ----- -------------- --- --- ------\n",
"\t index 2.000e+00 nan nan nan False\n",
"\tamplitude 1.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n"
]
}
],
"source": [
"pwl = models.PowerLaw()\n",
"print(pwl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This will use default values for the model parameters, which is rarely what you want.\n",
"\n",
"Usually you will want to specify the parameters on object creation.\n",
"One way to do this is to pass `astropy.utils.Quantity` objects like this:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PowerLaw\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- --------- ----- -------------- --- --- ------\n",
"\t index 2.300e+00 nan nan nan False\n",
"\tamplitude 1.000e-12 nan cm-2 s-1 TeV-1 nan nan False\n",
"\treference 1.000e+00 nan TeV nan nan True\n"
]
}
],
"source": [
"pwl = models.PowerLaw(\n",
" index=2.3, amplitude=1e-12 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
")\n",
"print(pwl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you see, some of the parameters have default ``min`` and ``values`` as well as a ``frozen`` flag. This is only relevant in the context of spectral fitting and thus covered in [spectrum_analysis.ipynb](https://github.com/gammapy/gammapy/blob/master/tutorials/spectrum_analysis.ipynb). Also, the parameter errors are not set. This will be covered later in this tutorial."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Get and set model parameters\n",
"\n",
"The model parameters are stored in the ``Parameters`` object on the spectal model. Each model parameter is a ``Parameter`` instance. It has a ``value`` and a ``unit`` attribute, as well as a ``quantity`` property for convenience."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Parameters\n",
"Parameter(name='index', value=2.3, factor=2.3, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n",
"Parameter(name='amplitude', value=1e-12, factor=1e-12, scale=1.0, unit='cm-2 s-1 TeV-1', min=nan, max=nan, frozen=False)\n",
"Parameter(name='reference', value=1.0, factor=1.0, scale=1.0, unit='TeV', min=nan, max=nan, frozen=True)\n",
"\n",
"covariance: \n",
"None\n"
]
}
],
"source": [
"print(pwl.parameters)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Parameter(name='index', value=2.3, factor=2.3, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n",
"Parameter(name='index', value=2.6, factor=2.6, scale=1.0, unit='', min=nan, max=nan, frozen=False)\n"
]
}
],
"source": [
"print(pwl.parameters[\"index\"])\n",
"pwl.parameters[\"index\"].value = 2.6\n",
"print(pwl.parameters[\"index\"])"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Parameter(name='amplitude', value=1e-12, factor=1e-12, scale=1.0, unit='cm-2 s-1 TeV-1', min=nan, max=nan, frozen=False)\n",
"Parameter(name='amplitude', value=2e-12, factor=2e-12, scale=1.0, unit='m-2 s-1 TeV-1', min=nan, max=nan, frozen=False)\n"
]
}
],
"source": [
"print(pwl.parameters[\"amplitude\"])\n",
"pwl.parameters[\"amplitude\"].quantity = 2e-12 * u.Unit(\"m-2 TeV-1 s-1\")\n",
"print(pwl.parameters[\"amplitude\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## List available models\n",
"\n",
"All spectral models in gammapy are subclasses of ``SpectralModel``. The list of available models is shown below."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[gammapy.spectrum.models.ConstantModel,\n",
" gammapy.spectrum.models.CompoundSpectralModel,\n",
" gammapy.spectrum.models.PowerLaw,\n",
" gammapy.spectrum.models.PowerLaw2,\n",
" gammapy.spectrum.models.ExponentialCutoffPowerLaw,\n",
" gammapy.spectrum.models.ExponentialCutoffPowerLaw3FGL,\n",
" gammapy.spectrum.models.PLSuperExpCutoff3FGL,\n",
" gammapy.spectrum.models.LogParabola,\n",
" gammapy.spectrum.models.TableModel,\n",
" gammapy.spectrum.models.ScaleModel,\n",
" gammapy.spectrum.models.AbsorbedSpectralModel,\n",
" gammapy.spectrum.crab.MeyerCrabModel]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"models.SpectralModel.__subclasses__()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting\n",
"\n",
"In order to plot a model you can use the ``plot`` function. It expects an energy range as argument. You can also chose flux and energy units as well as an energy power for the plot"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8FWXa//HPlYQk1EiX3kKHgEgvwUIXBEFX1/azK4og7K79Wd1ldXVLwIIFC1hWRAQUEKSomwDSQTrSkV6kCoR6//7IYZ8sj5BDTpmck+/79crLc+bMmflmd/TKPffMNeacQ0REJLdivA4gIiKRTYVEREQCokIiIiIBUSEREZGAqJCIiEhAVEhERCQgKiQiIhIQFRIREQmIComIiAREhURERAIS53WAcChVqpSrWrWq1zFERCLKokWL9jnnSue0XlQXEjPrAfRITk5m4cKFXscREYkoZrbFn/Wi+tSWc26ic+6BpKQkr6OIiEStqC4kIiISeiokIiISEBUSEREJiAqJiIgERIVEREQCokJyEc45xi7axsnTZ72OIiKSZ0V1ITGzHmY2/NChQ7n6/rxN+/ndmKX0eG0WP2w9GOR0IiLRIaoLSaD3kbSsXpJ372zKoeOn6P3GbF74ahXHT54JckoRkcgW1YUkGDrUK8u0Qanc0rwy78zcROehGXy/YZ/XsURE8gwVEj8USyzAizc0ZNT9LTGDW9+Zx1PjlnE485TX0UREPKdCcgla1SjJ1wNSeSC1OqMXbKVjWjozVu32OpaIiKdUSC5RwfhYnu5Wl/EPt6F4oXju+3Ah/Uct4edfTngdTUTEEyokudSo0mVM6NeWQR1rMWXFTjqkpfPlD9txznkdTUQkrFRIAhAfF0P/a2vyVf92VClZmAGf/sC9Hyxk56HjXkcTEQkbFZIgqFW2KGP7tuZ/utdjzoaf6ZiWwb/mbeHsWY1ORCT6qZAESWyMcW/bakx9LJVGlZJ4ZvwKbn13Lpv3HfU6mohISEV1IQn0zvbcqFyyEB/f24KX+zRk5Y7DdB6awfCMDZw+ozYrIhKdLD9MDjdt2tR58ajd3YczefaLFUxftZuUikm83CeFuuWKhT2HiEhumNki51zTnNaL6hGJ18oWS2T4HVfy+q1XsP3AcXq8Nou06Ws5cVptVkQkeqiQhJiZ0T2lPDMGtef6RuV59Zt1dH91Fot/OuB1NBGRoFAhCZPiheNJu7kxI+5uxtETp+nz5vf8eeIqjp087XU0EZGAqJCE2dW1yzB1YCq3t6jC+7OzmkDOXq8mkCISuVRIPFA0sQCDezVg9AMtiYuJ4bZ35/HE58s4dFxNIEUk8qiQeKhF9ZJMGdCOh9rX4PPF2+iYls7Ulbu8jiUicklUSDyWWCCWJ7vW4YuH21CySAIPfrSIRz5ZzN4jagIpIpFBhSSPaFgxiQn92vD7TrWYvnI3HYekM27xNjWBFJE8T4UkDykQG0O/a2oyeUBbqpcqzKDPlnL3yAVsP6gmkCKSd0V1IfGiRUowJJcpypiHWvN8j3rM37SfTmnpfDRns5pAikiepBYpedzW/cd4evxyZq7bR/OqJXipT0Oqly7idSwRyQfUIiVKVCpRiA/vac7fbkxhza7DdHllJm/+W00gRSTvUCGJAGbGb5pWYsag9lxduzQvf72GXm/MZuWOyDplJyLRSYUkgpQplsjbdzTlzduasOvQCa5/fTZ/n7qGzFNqAiki3lEhiUBdG5ZjxqBUejWuwLDvNnDdqzNZtGW/17FEJJ9SIYlQlxWK55+/acQH9zQn89RZbnxrDs9PWMnRE2oCKSLhpUIS4drXKs3Uganc2bIKH8zZTKchGWSs3et1LBHJR1RIokCRhDj+1LMBYx5sRUKBGO58fz6/H7OUQ8fUBFJEQk+FJIo0rVqCyf3b8cjVNRi/ZDsdhqTz9YqdXscSkSinQhJlEgvE8ofOdZjQrw1liibw0MeL6fvxIvYcyfQ6mohEKRWSKFW/fBJfPNKGx7vU5ps1e+iYlsGYhVvVBFJEgk6FJIoViI3h4auSmTKgHbXKFuEPny/jzvfns3X/Ma+jiUgUiepCEqlNG4OtRukijH6gFX/uWZ/FWw7QeWgGI2dvUhNIEQkKNW3MZ7YdOMbT41eQsXYvTasU56U+KSSXURNIEfm/1LRRflXF4oX44O5m/POmRqzb8wvdXpnJsO/Wc0pNIEUkl1RI8iEzo8+VFZkxqD0d65Xl71N/pOfrs1mxPX+fAhSR3FEhycdKF01g2G1NeOv2K9n7ywl6DpvNy1+rCaSIXJq4i31oZil+bOOUc251kPKIB7o0uJxW1UvywuRVvPnvDUxdsYuX+qTQvFoJr6OJSAS46GS7mR0BlgB2kW1Ucs5VDXKuoNJku/9mrtvLU+OWs+3Ace5sVYXHu9ShSMJF/94QkSjl72R7Tv+FWOKcS81hRxmXlEzytHY1SzP1sVT+OW0tI77fxIxVu3mhd0Ourl3G62gikkdddI4kpyLi7zoSWQonxPHHHvX4/KHWFEqI4+4RCxg0+gcOHD3pdTQRyYP8mmy3LA3NrLOZpZpZyVAHE+9dWaU4X/Vvy6PXJDNh6Q46Dkln0rIdarMiIv/looXEzKqa2RvABmAocDcwCMgws9lmdoeZXWz+RCJcQlwsv+tUmwn92nJ5UiL9PlnCgx8tYs9hNYEUkSw5TbZ/BrwJpDvnzp73WTngNmCfc25kKEMGSpPtwXH6zFnenbWJIdPXEh8Xw/9cV4+bmlZEf0uIRCd/J9vVIkUu2ca9v/DkuOXM37SftsmlePGGhlQuWcjrWCISZEFtkWJmvc2sqO/1k2b2mZk1DjSkRKbqpYvw6f0t+UuvBvyw9SCdh2bw3qxNnFETSJF8yd872593zh0xs9ZAD2A08FboYkleFxNj3N6yCtMGptKiegkGT1rFjW99z7rdR7yOJiJh5m8hOdczozvwhnNuLJAQmkgSScpfVpARdzVj6M2N2bzvKNe9OotXv1nHydNqAimSX/hbSHaa2TDgZmCymcVfwnc9o+eRhIeZ0euKCkwf1J5O9cuSNn0t178+i2XbDnodTUTCwK/JdjMrAnQDljnn1phZeaCRc25KqAMGgybbw2v6qt08+8Vy9h45wf3tqvNYh1oUjI/1OpaIXCJdtZWNCkn4HTp+ir9OXs2nC7ZStWQhXuqTQsvquo9VJJLowVbiqaSCBXipTwqf3NeCsw5uGT6XZ8Yv50jmKa+jiUiQqZBISLVOLsXXj7Xj3rbVGDX/JzoNyeDbNbu9jiUiQaRCIiFXKD6O/+lej7F9W1MkIY57Ri7ksU+XsF9NIEWiQk69tiqY2cdm9p2ZPW5mcdk+Gxv6eBJNrqhcnEn92zLg2ppMWraTDmnpTFiqJpAikS6nEcn7wFzgD0A14DszK+77rHoog0l0SoiLZWDHWkzq35ZKxQvSf9QS7v9wIbsOqQmkSKTKqZCUcc697pxb6JzrC7xLVuffaoD+jJRcq3N5McY93IZnutVl1vp9dExLZ9T8nzQ6EYlAORWSBDP7zx3szrkPyBqdTAcuD2UwiX6xMcb9qdX5ekAq9SsU46lxy7n1nXls+fmo19FE5BLkVEhGAK2yL3DOfQ3cAvwYqlCSv1QtVZhP7mvJX3s3ZMX2Q3QemsG7MzeqCaRIhNANiZKn7DqUybNfLGfG6j00qpjE325sRO3Li3odSyRfCuqd7WZWGegHVAX+c+WWc653ABnDRoUksjjnmLhsJ3+asJLDmad45OpkHr4qmfg4Xa0uEk7+FpK4nFbwmQB8SNbciNq6SkiZGdc3Kk/b5FL8eeJKhs5Yx+TlO/nbjY1oXOkyr+OJyHn8HZHMd841D0OekNCIJLJ9u2Y3z4xfwe7DmdzbthqDOtZWE0iRMAj2qa07gCrAVODEueXOuWWBhAwXFZLIdyTzFC9NWcO/5v1E5RKFeKlPQ1rXKOV1LJGoFuxCMhi4D9jI/57acs651IBShokKSfSYu/Fnnhy7jM0/H+O3zSvzVLc6FEss4HUskagU7ELyI5DinDuR48p5kApJdDl+8gxDZ6zlnZkbKV00gRd6NaRDvbJexxKJOsFuI78M0DWYkicUjI/lqW51+eKRNhQvFM99Hy7k0VFL+PmXiPw7RyTi+XvVVklgjZnN47/nSCLi8l+JTikVL2NCv7a8lb6B175dx6x1e3muR316Ni6PmXkdTyTf8PfU1rW/ttw5903QEwWRmfUAeiQnJ9+/bt06r+NICK3dfYTHP1/GD1sPck2dMrxwQwPKJRX0OpZIRAvFDYl7nHOZvvcFgVLOua0BJw0DzZHkD2fOOkZ+v5l/TP2R2BjjqW51+G2zysTEaHQikhvBniMZx3/fiHgW0PNIJE+JjTHubVuNqY+l0qhSEs+MX8Fv35nLpn1qAikSSv4Wkjjn3H8eZ+e7eivhIuuLeKZyyUJ8fG8L/tYnhVU7D9NlaAZvp2/g9Bk1ZRAJBX8Lyc9m1u3cGzPrDuwPTSSRwJkZv2lWiRmD2tO+Vmn+OmUNvd/8ntU7D3sdTSTq+DtHUhMYRdbVWw7YB9zmnIuIGWzNkeRvzjkmL9/FcxNWcPDYKR6+qgaPXJNMQpzarIhcTFAn27Nt9DIA59zBALKFnQqJABw4epLBX61i3OLtJJcpwst9UriySvGcvyiSTwVlst3MbrFsF+Q75w5mLyJmVtXMWgcWVSQ8iheOJ+03jRlxdzOOnTjNjW99z58nruLYydNeRxOJaDndkFgBWGJm84FFwF4gEUgGrgIOA0+EMqBIsF1duwzTBrXnb1+v4f3Zm5i+ehcv9U6hTbKaQIrkRo6ntswsDugItAHKAceB1cBk59ymkCcMAp3akguZv2k/T4xdxqZ9R7m5aSWevq4uSQXVBFIEQjRHEqlUSORiMk+d4ZVv1jE8YyMlC8czuFcDOte/3OtYIp4L9g2JIlErsUAsT3SpwxcPt6FkkQQe/GgRj/xrMXuPqAmkiD9USER8GlZMYkK/Nvyhc22mr9pNxyHpjFu8jfwwahcJRE5XbTULVxCRvKBAbAyPXJ3M5AFtqV6qMIM+W8pdIxaw/eBxr6OJ5Fk5jUgeNLMlZvaxmd1uZqXDkkrEY8llijLmodY836MeCzbvp1NaOh/O2czZsxqdiJzP3zvbGwBdgc5kXf77LfA1MNc5l+cbGGmyXQKxdf8xnh6/nJnr9tGsanFe6pNCjdJFvI4lEnIhu2rLzAoD15JVWFo455rkLmL4qJBIoJxzjF28ncGTVnH81Bke61CT+9tVp0Csphkleuny32xUSCRY9hzJ5LkvVzJlxS7qly/Gy31SaFAhyetYIiGhy39FQqBM0UTevP1K3rytCbsPn6DnsNn8feoaMk+d8TqaiGdUSERyoWvDcswYlMoNV1Rg2Hcb6PbqTBZu1pMVJH/K6fLfoWbWPFxhRCLJZYXi+cdNjfjwnuacOHWWm96ew3NfruCXE2oCKflLTiOSrcAwM9tgZi/4rt4SkWxSa5Vm2sBU/l+rqnw4dwudh2SQsXav17FEwuaihcQ590/nXDOgE3AMGGVmK8zsaTOrHpaEIhGgcEIcz19fnzEPtiKhQAx3vj+f349ZysFjJ3P+skiEy83lv1cC7wIpzrmIeMScrtqScMo8dYbXvl3HW+kbKV4onr/0qk+XBuW8jiVyyYJ61ZaZxZpZVzP7APgK2AjcHGBGkaiUWCCWP3Suw4R+bShbLIGHPl5M348XsedIptfRREIip8n2q81sOLAd6E/WHe01nXN9nHOfhyOgSKSqXz6JLx5pw+NdavPNmj10TMtgzMKtagIpUSenEcmfgSVAQ+dcV+fcB865I2HIJRIVCsTG8PBVyUwZ0I5aZYvwh8+Xcef789m6/5jX0USCJqfJ9nbOuTedc3vNrKWZ3QlgZiXNrHJ4IopEvhqlizD6gVYM7lmfxVsO0HloBiNnb1ITSIkK/s6RPAs8BzzrW5QIfBKqUCLRKCbGuKNVVaYNak+zqiV4fuIqbnp7Duv3aJAvkc3fO9tvBLoBRwGcc9uBYqEKFSxm1sPMhh86dMjrKCL/UeGygoy8uxlpv2nEhr2/0O2VWbz+7TpOncnzjbRFfpW/heSEy5ohdABmVih0kYLHOTfROfdAUpKa6kneYmb0blKR6QPb07FeWf4xbS3Xvz6bFdv1R49EHn8LyTgzGwYkmdndwDTg/dDFEskfShdNYNhtTXjr9ivZ90tWE8iXpqgJpEQWv29INLOuZN3hbsBU59yUUAYLJt2QKJHg0LFTvDh5NaMXbqVaqcK83CeF5tVKeB1L8rGgPI/EzKY55zoFNZkHVEgkksxat4+nxi9j6/7j3NGyCk90rUORhDivY0k+FKw72/WMdpEwa1uzFFMfS+WeNtX4eN4WOqWl892Pe7yOJXJBOY1INgK/v9DnzrlxoQgVbBqRSKRatOUAT4xdxvo9v9D7igr8T/d6FC8c73UsySf8HZHkNF5OArqTNS9yPgdERCERiVRXVinOV/3bMuzb9bzx7w1krNvL89fX57qG5TD7tX8tRcIvpxHJYudckzDmCQmNSCQarN55mCfGLmPZtkN0qleWwb0aULZYotexJIoFa45Ef/KI5BF1yxVjXN/WPNW1Dulr99IhLZ3RC35SE0jxXE6F5I6wpBARv8TFxvBg+xp8/VgqdcsV44mxy7n9vXn89LOaQIp3ciokL+W0ATObFKQsIuKnaqUK8+n9LflLrwYs3XqIzkMzeG/WJs6oCaR4IKc5koNAxsW+D9R3zuXpx+5qjkSi2Y6Dx3lm/HK++3EvV1S+jL/1SaFm2aJex5IoEKwbEtv7sa+Tzrk5lxIu3FRIJNo55/jyhx38aeJKjp44Q79rknmofQ3i4/ztgiTyfwWlkEQLFRLJL/b9coLnJ6xk0rKd1Lm8KH+7MYWUipd5HUsiVFCf2S4ikaFUkQRev7UJ79zZlAPHTtJr2GxenLya4yfVBFJCR4VEJAp1rFeW6YPac3OzSgzP2EjXVzKYu/Fnr2NJlPL3CYllfmVZ7eDHEZFgKZZYgL/2TuGT+1pw1sEtw+fyzPjlHMk85XU0iTL+jkhmmtlvzr0xs98B40MTSUSCqXVyVhPI+9pWY9T8n+g0JINv1+z2OpZEEX8LyVXAHWY2xswygFpA85ClEpGgKhgfy7Pd6zG2b2uKJsZxz8iFDPh0CfuPnvQ6mkQBvwqJc24n8DXQCqgKfOic+yWEuUQkBK6oXJxJj7ZjwLU1mbx8Jx3S0pmwdIfarEhA/J0jmQ60ABoA3YAhZvaPUAYTkdCIj4thYMdaTHq0HZVKFKL/qCXc/+FCdh3K9DqaRCh/T20Nc87d6Zw76JxbAbQGDoUwl4iEWO3LizKub2ueva4us9bvo2NaOqPmqwmkXDrdkCgibPn5KE+OXc6cjT/TqnpJXurTkColC3sdSzwW1BsSzeyImR32/WSa2Rkz04hEJEpUKVmYT+5vwYs3NGTF9qwmkO9kbFQTSPGLv5PtRZ1zxXw/iUAfYFhoo4lIOJkZt7aozPRB7WmbXIoXJq+m9xuzWbPrsNfRJI/L1Z3tzrkvgGuCnEVE8oDLkxJ5586mvHJLY7YeOE6P12YxZPpaTp4+63U0yaNyemY7AGbWO9vbGKApWc9sF5EoZGb0bFyBdjVL8+eJK3nlm3VMWbGTl/ukcEXl4l7HkzzG3xFJj2w/nYEjQM9QhRKRvKFE4XiG3nIF79/VlCOZp+n95vcMnrSKYydPex1N8hBdtSUifjmSeYqXv17Dx3N/onKJQrzUuyGtk0t5HUtCKFgPtnqNi5zCcs71z1288FIhEQmeuRt/5smxy9j88zFuaVaJp6+rS7HEAl7HkhDwt5DkNEei//qKyH9pWb0kXz+WypAZa3knYyPf/biHv/RqSMd6Zb2OJh7JaUQS55yL+JOhGpGIhMaybQd5/PNlrNl1hO4p5Xj++vqUKpLgdSwJkmDdkDg/2wZfCziViESVlIqXMaFfWwZ1rMXUlbvomJbO+CXb1GYln8mpkFi2121CGUREIlN8XAz9r63JV/3bUbVUYQaOXso9Ixew4+Bxr6NJmORUSPRnhYj4pVbZonz+UGv+2L0eczfup9OQDD6au4WzarMS9XKaIzkGrCdrZFLD9xrfe+ecSwl5wiDQHIlIeG3df4ynxi1n1vp9NK9Wgpf7pFCtlJpARppgXf5b5WJfds5tyUW2sFMhEQk/5xxjFm5j8FerOHn6LAM71uK+ttWIi81VZybxQFAKSbRQIRHxzu7DmTz7xQqmr9pNwwpJvNwnhXrli3kdS/wQ1DbyIiK5VbZYIsPvuJJhtzZh56HjXP/6LP457UdOnD7jdTQJEhUSEQk5M+O6lHJMH9ie6xuV57Vv13Pdq7NYtOWA19EkCHIsJGaW4vtnw9DHEZFoVrxwPGk3N2bE3c04duI0N771PX+auFJNICOcPyOSe8ysJnBvqMOISP5wde0yTBvUnttbVGHE7M10GpLBrHX7vI4luXTRQmJmz/nWmQvEmNkfw5JKRKJekYQ4BvdqwGcPtiI+Nobb35vHH8Ys5dCxU15Hk0t00ULinPsTMAMYDcxwzv05LKlEJN9oXq0Ekwe0o+9VNRi3ZDsdhqQzdeUur2PJJfDn1FYL59zDQLNQhxGR/CmxQCxPdKnDl4+0oXSRBB78aBGP/Gsxe4+c8Dqa+CHP30diZtWBZ4Ak59yNvmUxwGCgGLDQOffBxbah+0hEIsepM2cZnrGRV2aso2B8LH/sXo/eTSpgZjl/WYIqT9xHYmbvm9keM1tx3vIuZvajma03sycvtg3n3Ebn3PkT/T2BCsApYFtwU4uIlwrExvDI1clMHtCO5DJF+N2Ypdw1YgHb1QQyzwr1fSQjgS7ZF5hZLDAM6ArUA35rZvXMrKGZTTrvp8wFtlsbmOOcGwT0DWF+EfFIcpkijHmwFc/3qMeCzfvplJbOh3M2qwlkHhTSQuKcywD2n7e4ObDeN9I4CXwK9HTOLXfOdT/vZ88FNr0NOHcnk26PFYlSMTHGXW2qMfWxVJpUKc4fv1zJzcPnsGHvL15Hk2xyuvw31sweNLPBZtbmvM+ezeU+KwBbs73f5lt2oQwlzewt4Aoze8q3eBzQ2fewrYwLfO8BM1toZgv37t2by6gikhdUKlGID+9pzt9vTOHHXUfo+spM3vj3ek6dOet1NCHn7r/vAoXIelLiHUC673QSZrbYOdckxx2YVQUmOeca+N7fBHR2zt3ne38H0Nw592hgv8qFabJdJHrsOZLJc1+uZMqKXdQvX4yX+6TQoEKS17GiUrAm25s75251zg0FWgBFzGycmSXw309PvBTbgErZ3lcEduRyWyKSz5Qpmsibt1/Jm7c1YffhE/QcNpu/T11D5imd5fZKToUk/twL59xp59wDwA/At0CRXO5zAVDTzKqZWTxwCzAhl9sSkXyqa8NyzBiUyg1XVGDYdxvo9upMFm4+f0pWwiGnQrLQzP7rqivf3e0jgKo5bdzMRgFzgNpmts3M7nXOnQb6AVOB1cBnzrmVuQkvIvnbZYXi+cdNjfjwnuacOHWWm96ew3NfruCXE2oCGU55/obEYNAciUj0O3riNH+f+iMfzNlM+aSC/LV3Q1JrlfY6VkQLyhyJmT2e7fVN5332Yu7jiYgEV+GEOJ6/vj6fP9SKxAIx3Pn+fH4/ZikHj530OlrUy+nU1i3ZXj913mddyOPMrIeZDT906JDXUUQkTK6sUoKv+rej39XJjF+ynQ5pGUxZvtPrWFEtp0JiF3j9a+/zHOfcROfcA0lJujRQJD9JLBDL7zvXZkK/NlyelEDffy3moY8WsedwptfRolJOhcRd4PWvvRcRyVPql0/ii4fb8ESXOnz74x46pKXz2cKt5Ie54XDKqZA0MrPDZnYESPG9Pvdej94VkTwvLjaGvlfVYMqAdtS+vCiPf76MO9+fz9b9x7yOFjVyerBVrHOumHOuqHMuzvf63PsC4QopIhKoGqWLMPqBVgzuWZ/FWw7QeWgGI2dv4oyaQAYs1N1/RUTyjJgY445WVZk2qD3Nqpbg+YmruOmt71m/54jX0SKaComI5DsVLivIyLubkfabRmzcd5Rur8zi9W/XqQlkLqmQiEi+ZGb0blKR6QPb07F+Wf4xbS09XpvF8m26XeBSqZCISL5WumgCw25twtt3XMn+oyfp9cZsXpqiJpCXIqoLiW5IFBF/da5/OdMHtefGJhV5K30DXV+ZybyNP3sdKyKo15aIyHlmr9/Hk+OWsXX/ce5oWYXHu9SmaGL+u1A1WM8jERHJd9okl2LqY6nc06YaH8/bQuchGXz344We/C0qJCIiv6JQfBx/7FGPsX1bUzghjrtHLGDQ6B84cFRNIM+nQiIichFNKhdnUv+29L+2JhOW7qBDWjqTlu1Qm5VsVEhERHKQEBfLoI61mPhoWyoUL0i/T5bw4EeL2K0mkIAKiYiI3+qWK8a4vq15ulsd0tfupUNaOqMX/JTvRycqJCIilyAuNoYHUmsw9bFU6pUrxhNjl3P7e/P46ef82wRShUREJBeqlirMqPtb8sINDVi69RCdh2bw7syN+bIJZFQXEt2QKCKhFBNj3NaiCtMHpdKqRkn+8tVq+rz5PWt3568mkLohUUQkCJxzTFi6gz9NXMWRzFM8ek1NHmpfg/i4yP17XTckioiEkZnRs3EFpg9MpUuDcqRNX8v1r89i6daDXkcLORUSEZEgKlkkgdd+ewXv3NmUA8dOcsMbs3lx8mqOn4zeJpAqJCIiIdCxXlmmD2rPLc0rMzxjI11fyWDOhuhsAqlCIiISIsUSC/DiDQ355P4WOOC378zl6fHLOZJ5yutoQaVCIiISYq1rlOLrAak8kFqdT+f/RKchGXy7ZrfXsYJGhUREJAwKxsfydLe6jHu4DcUSC3DPyIUM+HQJ+6OgCaQKiYhIGDWudBkTH23LwA61mLx8Jx3S0vnyh+0R3WZFhUREJMzi42IY0KEmkx5tR6UShRjw6Q/c98FCdh467nW0XFEhERHxSO3LizKub2ueva4uszfso1NaBp/M+4mzEdZmJaoLiVqkiEheFxtj3NeuOlMfS6UdpruoAAAH5klEQVRBhSSeHr+cW9+dy+Z9R72O5je1SBERySOcc4xesJUXJq/m1Jmz/K5jbe5pW43YGPMkj1qkiIhEGDPjluaVmT6wPW2TS/PC5NX0fmM2a3Yd9jraRamQiIjkMZcnJfLOnVfy+q1XsO3Acbq/Oou06Ws5cTpvtllRIRERyYPMjO4p5Zk+qD09GpXn1W/W0eO1WSz56YDX0f4PFRIRkTysROF4htzcmBF3NeNI5ml6v/k9gyet4tjJ015H+w8VEhGRCHB1nTJMG5jKbS0q896sTXQZOpPv1+/zOhagQiIiEjGKJhbgL70a8ukDLYkxuPXdeTw5dhmHjnvbBFKFREQkwrSsXpKvH0vlwfbV+WzhVjqmpTNt5S7P8qiQiIhEoMQCsTzVtS5fPNKGEoXjeeCjRfT7ZDH7fjkR9iwqJCIiESylYlYTyN91rMW0lbvpkJbO+CXbwtoEUoVERCTCFYiN4dFra/JV/7ZUK1WYgaOXcs/IBew4GJ4mkCokIiJRombZonz+UGue61GPuRv30zEtnY/nbgn5fqO6kKhpo4jkN7Exxt1tqjFtYCpXVC7Out1HQr5PNW0UEYlSzjlOnXHEx+VuzOBv08a4XG1dRETyPDMjPi70nYOj+tSWiIiEngqJiIgERIVEREQCokIiIiIBUSEREZGAqJCIiEhAVEhERCQg+eKGRDPbCxwELnaLe9JFPi8F5I0nyPjvYr9PXt5XINu61O/6u74/6+W0TrQdXxC+Y0zHl3fHVxXnXOkc13LO5YsfYHhuPwcWep0/2L9vXt1XINu61O/6u74/6+W34yvY/7+Haz86vkLzk59ObU0M8PNIE87fJ5j7CmRbl/pdf9f3Z738dnxB+H4nHV95/PjKF6e2AmVmC50f/WZEckPHl4RSOI6v/DQiCcRwrwNIVNPxJaEU8uNLIxIREQmIRiQiIhIQFRIREQmIComIiAREheQSmVkvM3vHzL40s05e55HoYmZ1zewtM/vczPp6nUeik5kVNrNFZtY9GNtTIQHM7H0z22NmK85b3sXMfjSz9Wb2JIBz7gvn3P3AXcDNHsSVCHOJx9dq59xDwG8AXRIsfrmUY8znCeCzYO1fhSTLSKBL9gVmFgsMA7oC9YDfmlm9bKs86/tcJCcjuYTjy8yuB2YB34Q3pkSwkfh5jJlZB2AVsDtYO1chAZxzGcD+8xY3B9Y75zY6504CnwI9LcvLwBTn3OJwZ5XIcynHl2/9Cc651sBt4U0qkeoSj7GrgZbArcD9ZhZwHYgLdANRrAKwNdv7bUAL4FGgA5BkZsnOube8CCcR71ePLzO7CugNJACTPcgl0eNXjzHnXD8AM7sL2OecOxvojlRILsx+ZZlzzr0KvBruMBJ1LnR8/Rv4d3ijSJT61WPsPy+cGxmsHenU1oVtAyple18R2OFRFok+Or4k1MJ2jKmQXNgCoKaZVTOzeOAWYILHmSR66PiSUAvbMaZCApjZKGAOUNvMtpnZvc6500A/YCqwGvjMObfSy5wSmXR8Sah5fYypaaOIiAREIxIREQmIComIiAREhURERAKiQiIiIgFRIRERkYCokIiISEBUSCRfM7MzZvZDtp8nc/5W6JnZZjNbbmZNfe/jzOxFM1uXLeszOWxjpJk9eN6yXmY22cwK+rZx0sxKhfJ3keinXluS3x13zjUO5gbNLM53M1igrnbO7fO9/gtwOdDQOZdpZkWB3+Xw/VHAk8Db2ZbdAoxyzh0HGpvZ5iDklHxOIxKRX+EbEfzJzBb7RgZ1fMsL+x4itMDMlphZT9/yu8xsjJlNBKaZWYyZvWFmK81skm8UcKOZXWtm47Ptp6OZjcshSyHgfuBR51wmgHPuiHPu+Wzr3G5m832jjLd9z6KYAdQxs3LZttMB+CKY/1uJqJBIfnfuFM+5n+xPvdznnGsCvAn83rfsGeBb51wzsp7r8HczK+z7rBXw/5xz15DVCr4q0BC4z/cZwLdAXTMr7Xt/NzAih4zJwE/OuSO/9qGZ1SXraZ1tfKOrM8BtzrkzwDiynrYIcD3w3YW2I5JbKiSS3x13zjXO9jM622fnRgqLyCoKAJ2AJ83sB7LavScClX2fTXfOnXu4UFtgjHPurHNuF/AdZPWJBz4Cbjezy8gqMFMuJbCZ3e0relvNrBJwLXAlsMCX61qgum/1UWSdzsL3z1GXsi8Rf2iOROTCTvj+eYb//XfFgD7OuR+zr2hmLYCj2RddZLsjgIlAJlnFJqf5lPVAZTMr6julNQIY4Xs+d6xvXx845576le/OBsqZWSOgNf9bVESCRiMSkUszFXjUzAzAzK64wHqzgD6+uZKywFXnPnDO7SDruRDPkvWs7Ytyzh0D3gNeN7NE335jgXjfKt8AN5pZGd9nJcysiu+7DvgM+ACYfG6ORSSYVEgkvzt/juSlHNYfDBQAlvlGBIMvsN5Ysh4stIKsq6bmAYeyff4vYKtzbpWfOZ8BdgIrzGwJMJOs4rDDt41nyZrkXwZMB8pl++4ooBFZz+wWCTq1kRcJETMr4pz7xcxKAvPJmgzf5fvsdWCJc+69C3x3M9A02+W/ocoYlv1IdNOIRCR0Jvkmv2cCg7MVkUVACvDxRb67F/jm3A2JwXbuhkSyRldnQ7EPyT80IhERkYBoRCIiIgFRIRERkYCokIiISEBUSEREJCAqJCIiEhAVEhERCcj/B164BYhSJtPWAAAAAElFTkSuQmCC\n",
"text/plain": [
"