\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/v0.18.2?urlpath=lab/tree/models.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",
"[models.ipynb](../_static/notebooks/models.ipynb) |\n",
"[models.py](../_static/notebooks/models.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Gammapy Models\n",
"\n",
"\n",
"This is an introduction and overview on how to work with models in Gammapy. \n",
"\n",
"The sub-package `~gammapy.modeling` contains all the functionality related to modeling and fitting\n",
"data. This includes spectral, spatial and temporal model classes, as well as the fit\n",
"and parameter API. We will cover the follwing topics in order:\n",
"\n",
"1. [Spectral Models](#Spectral-Models)\n",
"1. [Spatial Models](#Spatial-Models)\n",
"1. [SkyModel](#SkyModel)\n",
"1. [Model Lists and Serialisation](#Model-Lists-and-Serialisation)\n",
"1. [Implementing as Custom Model](#Implementing-a-Custom-Model)\n",
"\n",
"The models follow a naming scheme which contains the category as a suffix to the class name. An overview of all the available models can be found in the [model gallery](https://docs.gammapy.org/dev/modeling/gallery/index.html#spectral-models).\n",
"\n",
"Note that there is a separate tutorial [modeling](modeling.ipynb) that explains about `~gammapy.modeling`,\n",
"the Gammapy modeling and fitting framework. You have to read that to learn how to work with models in order to analyse data.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from astropy import units as u\n",
"from gammapy.maps import Map, WcsGeom, MapAxis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spectral Models\n",
"\n",
"All models are imported from the `~gammapy.modeling.models` namespace. Let's start with a `PowerLawSpectralModel`:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from gammapy.modeling.models import PowerLawSpectralModel"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PowerLawSpectralModel\n",
"\n",
" name value unit min max frozen error \n",
"--------- ---------- -------------- --- --- ------ ---------\n",
" index 2.0000e+00 nan nan False 0.000e+00\n",
"amplitude 1.0000e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00\n"
]
}
],
"source": [
"pwl = PowerLawSpectralModel()\n",
"print(pwl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To get a list of all available spectral models you can import and print the spectral model registry or take a look at the [model gallery](https://docs.gammapy.org/dev/modeling/gallery/index.html#spectral-models):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Registry\n",
"--------\n",
"\n",
"ConstantSpectralModel : ['ConstantSpectralModel', 'const'] \n",
"CompoundSpectralModel : ['CompoundSpectralModel', 'compound'] \n",
"PowerLawSpectralModel : ['PowerLawSpectralModel', 'pl'] \n",
"PowerLaw2SpectralModel : ['PowerLaw2SpectralModel', 'pl-2'] \n",
"BrokenPowerLawSpectralModel : ['BrokenPowerLawSpectralModel', 'bpl'] \n",
"SmoothBrokenPowerLawSpectralModel : ['SmoothBrokenPowerLawSpectralModel', 'sbpl'] \n",
"PiecewiseNormSpectralModel : ['PiecewiseNormSpectralModel', 'piecewise-norm'] \n",
"ExpCutoffPowerLawSpectralModel : ['ExpCutoffPowerLawSpectralModel', 'ecpl'] \n",
"ExpCutoffPowerLaw3FGLSpectralModel : ['ExpCutoffPowerLaw3FGLSpectralModel', 'ecpl-3fgl'] \n",
"SuperExpCutoffPowerLaw3FGLSpectralModel: ['SuperExpCutoffPowerLaw3FGLSpectralModel', 'secpl-3fgl'] \n",
"SuperExpCutoffPowerLaw4FGLSpectralModel: ['SuperExpCutoffPowerLaw4FGLSpectralModel', 'secpl-4fgl'] \n",
"LogParabolaSpectralModel : ['LogParabolaSpectralModel', 'lp'] \n",
"TemplateSpectralModel : ['TemplateSpectralModel', 'template'] \n",
"GaussianSpectralModel : ['GaussianSpectralModel', 'gauss'] \n",
"EBLAbsorptionNormSpectralModel : ['EBLAbsorptionNormSpectralModel', 'ebl-norm'] \n",
"NaimaSpectralModel : ['NaimaSpectralModel', 'naima'] \n",
"ScaleSpectralModel : ['ScaleSpectralModel', 'scale'] \n",
"PowerLawNormSpectralModel : ['PowerLawNormSpectralModel', 'pl-norm'] \n",
"LogParabolaNormSpectralModel : ['LogParabolaNormSpectralModel', 'lp-norm'] \n",
"ExpCutoffPowerLawNormSpectralModel : ['ExpCutoffPowerLawNormSpectralModel', 'ecpl-norm'] \n",
"\n"
]
}
],
"source": [
"from gammapy.modeling.models import SPECTRAL_MODEL_REGISTRY\n",
"\n",
"print(SPECTRAL_MODEL_REGISTRY)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Spectral models all come with default parameters. Different parameter\n",
"values can be passed on creation of the model, either as a string defining\n",
"the value and unit or as an `astropy.units.Quantity` object directly:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"amplitude = 1e-12 * u.Unit(\"TeV-1 cm-2 s-1\")\n",
"pwl = PowerLawSpectralModel(amplitude=amplitude, index=2.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For convenience a `str` specifying the value and unit can be passed as well:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PowerLawSpectralModel\n",
"\n",
" name value unit min max frozen error \n",
"--------- ---------- -------------- --- --- ------ ---------\n",
" index 2.2000e+00 nan nan False 0.000e+00\n",
"amplitude 2.7000e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00\n"
]
}
],
"source": [
"pwl = PowerLawSpectralModel(amplitude=\"2.7e-12 TeV-1 cm-2 s-1\", index=2.2)\n",
"print(pwl)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The model can be evaluated at given energies by calling the model instance:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[2.70000000e-12 2.40822469e-13 1.70358483e-14 1.51948705e-15] 1 / (cm2 s TeV)\n"
]
}
],
"source": [
"energy = [1, 3, 10, 30] * u.TeV\n",
"dnde = pwl(energy)\n",
"print(dnde)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The returned quantity is a differential photon flux. \n",
"\n",
"For spectral models you can computed in addition the integrated and energy flux\n",
"in a given energy range:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.108034597491956e-12 1 / (cm2 s)\n",
"4.982075849517389e-12 TeV / (cm2 s)\n"
]
}
],
"source": [
"flux = pwl.integral(energy_min=1 * u.TeV, energy_max=10 * u.TeV)\n",
"print(flux)\n",
"\n",
"eflux = pwl.energy_flux(energy_min=1 * u.TeV, energy_max=10 * u.TeV)\n",
"print(eflux)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This also works for a list or an array of integration boundaries:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1.64794383e-12 4.60090769e-13 1.03978226e-13] 1 / (cm2 s)\n"
]
}
],
"source": [
"energy = [1, 3, 10, 30] * u.TeV\n",
"flux = pwl.integral(energy_min=energy[:-1], energy_max=energy[1:])\n",
"print(flux)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In some cases it can be useful to find use the inverse of a spectral model, to find the energy at which a given flux is reached:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.0 TeV\n"
]
}
],
"source": [
"dnde = 2.7e-12 * u.Unit(\"TeV-1 cm-2 s-1\")\n",
"energy = pwl.inverse(dnde)\n",
"print(energy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a convenience you can also plot any spectral model in a given energy range:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3zV5fn/8deVhD0iU9lhCQIB2bK1ZSpDAQdupVKsKMK3KtbaatXiaHEiiFWxWrEKqIAi4CKCiAxlr4CMiAwR2TvX74+E/tI0O+fk5Jy8n49HHuTc53PuzxX9kIv7vj+f6zZ3R0REJDeiQh2AiIiEHyUPERHJNSUPERHJNSUPERHJNSUPERHJNSUPERHJtZhQB1AQKleu7HFxcaEOQ0QkrCxbtuwnd6+S0XtFInnExcWxdOnSUIchIhJWzGxbZu9F9LSVmfUzs0kHDhwIdSgiIhElopOHu89092GxsbGhDkVEJKJEdPIQEZHgUPIQEZFcU/IQEZFci+jkoQVzEZHgiOjkkd8Fc3dn+vIkTpw+E+DIRETCW0Qnj/xatm0/o99ZQd/nFrB8+/5QhyMiUmgoeWShTVxFXru5LYdPnGbQhK94ZNZajp48HeqwRERCTskjG5c0rsrcUV25rn1tXlnwPb2eSWBh4k+hDktEJKSUPHKgXMliPHp5PP8edhExUVFc94/F3Dd1JQeOnQp1aCIiIaHkkQvt61Vi9sguDO9Wn6nLk+gxbj5z1+wKdVgiIgVOySOXShaLZkyfxrz/u05UKluCYW8s4463lvPT4ROhDk1EpMBEdPII5nMe8TVjmTGiE7/veT7z1uym+7j5vPdtEu4e8HOJiBQ2VhR+2bVp08aDWZI9cc8h7p26kuXbf+HiRlV47Ip4apxTKmjnExEpCGa2zN3bZPReRI88CkqDquV4d3hH/tyvCYu3/EzPcfN5Y9FWkpMjPzGLSNGk5BEg0VHGLZ3qMndUV1rWrsCDH6zhmklfs2Xv4VCHJiIScEoeAVarYmneGNqOJwc3Z/2ug/R+9ksmfLGZ02eSQx2aiEjAKHkEgZlxVZtafDK6G5c0qsITH6/n8hcXsmanCjSKSGRQ8giiquVLMvH61rx4XSt2HThO/xcW8tSc9Rw/pUKLIhLelDyCzMy4NL4a80Z14/ILazD+881c9tyXLNv2c6hDExHJMyWPAlKhTHH+flULXr+1HcdPJTN44iIemrGGIydUaFFEwo+SRwHrdn4V5ozqyo0X1eH1RVvp+XQCCRv3hjosEZFcCYvkYWb1zOwVM5uapu1yM3vZzD4ws56hjC+3ypaI4eEBzXjntx0oUSyKG1/9ht+/u4Jfjp4MdWgiIjkS9ORhZq+a2R4zW52uvbeZbTCzRDMbk1Uf7r7F3Yema3vf3W8DbgauDnjgBaBtXEU+uqsLd1xSn/e+/YEeTyfw8eofQx2WiEi2CmLkMRnonbbBzKKB8UAfoAkwxMyamFm8mc1K91U1m/7/mNpXWCpZLJp7ejXmgzs6UaVsCYa/uZzb31zGnkPHQx2aiEimYoJ9AndPMLO4dM3tgER33wJgZm8DA9x9LNA3J/2amQGPA7PdfXngIg6NZjVi+WBEJyYlbOHZTzfx1eZ9PNi3CYNa1SDlRxURKTxCteZRA9iR5nVSaluGzKySmU0EWprZ/anNdwLdgcFmNjyDzwwzs6VmtnTv3vBYkC4WHcUdlzRg9sgunH9uWX7/7gpufPUbdvx8NNShiYj8lwKpqps68pjl7s1SX18J9HL336S+vgFo5+53BuP8wa6qGwzJyc6/Fm/j8dnrceDeXo24sUMcUVEahYhIwSiMVXWTgFppXtcEdgb6JMHczyPYoqKMGzrEMWdUV9rGVeShmWu56qVFJO5RoUURCb1QJY8lQEMzq2tmxYFrgBmBPom7z3T3YbGxsYHuusDUrFCaybe05e9XtiBx72EuffZLxn+eyCkVWhSRECqIW3WnAIuARmaWZGZD3f00MAKYA6wD3nH3NUE4d9iOPNIyMwa1rsm8Ud3o3qQqT83ZwIAXFrL6h/D+uUQkfGknwTD08epdPPjBan4+cpJhXesx8tcNKVksOtRhiUiEKYxrHpIPvZudxyejujGoVQ0mfLGZS5/9kiVbVWhRRApORCePSJm2ykhs6WI8ObgFbw5tz8kzyVw5cRF/+mA1h1VoUUQKgKatIsDRk6d5as4GJn+1leqxpXjsimZc3Ci7B/NFRLKmaasIV7p4DH/u15SpwztSqng0N7+2hNHvfMf+Iyq0KCLBEdHJI5KnrTLSuk4FPryrMyMuacCM73bS4+n5fLjyR4rC6FJECpamrSLU2p0HuW/aSlb9cICeTc7l0cubUbV8yVCHJSJhRNNWRVCT6uV573cdGdOnMfM37uXX4+bzzpIdGoWISEBEdPIoatNW6cVERzG8W31mj+zCBdXKc++0ldzwigotikj+adqqiEhOdt76ZjtjP1pHssM9vRpxU8c4olVoUUQyoWkrISrKuP6iOswd3Y329Sryl1lrGTzxKzbtPhTq0EQkDCl5FDE1zinFaze35ZmrL2TrT0e47LkFPPfpJk6eVqFFEck5JY8iyMy4vGUN5o3uRq9m5zFu3kb6v7CAlUm/hDo0EQkTEZ08ivqCeXYqly3B80Na8vKNbdh/9CSXj1/I2I/WcfzUmVCHJiKFnBbMBYADx07x+Ox1TPlmB3GVSvP4oOZcVK9SqMMSkRDSgrlkK7ZUMcYObM5bv2lPssM1k77mgfdWcej4qVCHJiKFkJKH/JeODSrz8d1d+E3nukz5Zjs9n07gs/W7Qx2WiBQySh7yP0oXj+GPfZsw7faOlCsZw62Tl3L329/yswotikiqiE4eWjDPn5a1KzDrzi6M/HVDZq38ke7j5jNjxU6VOBERLZhLzqzfdZD7pq5kRdIBul+QUmjxvFgVWhSJZFowl3xrfF55pt3ekQcuvYAFiXvp8fR83v5mu0YhIkWUkofkWEx0FLd1rcfHI7vStHp5xkxfxbUvL2bbviOhDk1ECpiSh+RaXOUyTLntIsYOjGf1Dwfo9UwC//hyC2eSNQoRKSpisnrTzCrmoI9kd1ddiyLGzBjSrjaXNKrKA++t4tEP1zFz5Y88Oag5jc4rF+rwRCTIslwwN7PjwE4gq7rd0e5eO9CBBZIWzIPL3Zm58kcemrGGQ8dPccclDfjdxQ0oHqOBrUg4y2rBPMuRB7DO3Vtm0/m3eY5MIoKZ0b9FdTo3qMzDM9fwzCebmL1qF08Mbs6Ftc4JdXgiEgTZ/dOwQw76yMkxIaHnPApWxTLFefaalrxyUxsOHDvFwBcX8tiHazl2UoUWRSJNjp/zMLMKQHXgGLDV3cNmAwhNWxW8Q8dP8fjs9fxr8XZqVyzN44Pi6Vi/cqjDEpFcyPNzHmYWa2Z/MLNVwNfAS8A7wDYze9fMLgl8uBIJypUsxmNXxDPltoswg2tfXsz901dxUIUWRSJCdtNWU4EdQBd3b+Tund29jbvXAh4HBpjZ0KBHKWGrQ/1KfDyyK8O61uPfS7bTY9x85q1VoUWRcKfyJFJgVuz4hfumrWT9rkP0a1Gdh/o1oVLZEqEOS0QykZ9pq7Vm9oCZ1Q9OaFKUtKh1DjNGdGZU9/P5eHVKocUPvvtBJU5EwlB201ZDgLLAXDNbbGZ3m1n1AohLIlTxmChGdm/Ih3d1oU6lMox8+zuGvr6Unb8cC3VoIpILWSYPd1/h7ve7e31gJFAH+NrMPjOz2wokQolI559bjmm3d+TBvk1YtHkfPZ9O4F+Lt5GsEiciYSHHjwC7+9fuPgq4EagAvBC0qKRIiI4yhnauy5y7u9KiViwPvLeaIS9/zfc/qdCiSGGXo+RhZm3NbJyZbQMeBiYBNYIamRQZtSuV5s2h7XliUDxrfzxI72cSeGn+Zk6fCZtHiUSKnOwWzP9qZpuBCaTUuOrk7t3cfYK7/1QQAZpZPTN7xcympmm7wMwmmtlUM7u9IOKQ4DIzrm5bm09Gd6Pr+VUYO3s9Ayd8xbofD4Y6NBHJQHYjjxNAn9RnO/7m7km56dzMXjWzPWa2Ol17bzPbYGaJZjYmqz7cfYu7D03Xts7dhwNXARneRibh6dzyJZl0Q2teuLYlO385Rr/nFzBu7gZOnFaJE5HCJLsF84fdfaOZlTazB83sZQAza2hmfXPQ/2Sgd9oGM4sGxgN9gCbAEDNrYmbxZjYr3VfVzDo2s/7AAuDTHMQhYcTM6Nu8OvNGdaN/i+o891kilz23gGXb9oc6NBFJldMF89dIGYWcLYKYBDya3YfcPQH4OV1zOyAxdURxEngbGODuq9y9b7qvPVn0PcPdOwLX5fBnkDBToUxxxl19Ia/d0pajJ04zeOJXPDxzDUdPng51aCJFXk6TR313fxI4BeDux8h6j4+s1CCl5MlZSWSx+G5mlcxsItDSzO5PbbvYzJ4zs5eAjzL53DAzW2pmS/fu3ZvHUKUwuKRRVeaM6sr17evw2sKt9Hw6gQWbCmTJTUQykd1+HmedNLNSgAOkPnF+Io/nzCjpZHpzv7vvA4ana/sC+CKrk7j7JFLuCqNNmzZ6eCDMlStZjEcub0bf5tUYM30V17+ymKva1OSBy5oQW6pYqMMTKXKyu9tqbuq3DwEfA7XM7F+krDPcm8dzJgG10ryuScqdXAGn/TwiT/t6lZg9sgvDu9Vn2vIf6DFuPnPW7Ap1WCJFTnbb0H57didBM6sEXETKyOHrnN6qa2ZxwCx3b5b6OgbYCPwa+AFYAlzr7mvy/mNkTYURI9OqpAPcO20l6348yGXx1Xiof1OqlFOhRZFAyXNhRCDWzAaa2UCgG1ACKA50TW3L7sRTgEVAIzNLMrOh7n4aGAHMAdYB7wQrcWjkEdnia8YyY0Qn7unViHlrd9N93HymLUtSoUWRApDdyGMf8AGZrFO4+63BCiyQNPKIfIl7DnHv1JUs3/4L3c6vwl8HxlPjnFKhDkskrGU18sgueSx391ZBi6yAKHkUDWeSnTcWbeXJORsw4L4+jbm+fR2iovJ6Y6BI0Zafaauw/lunaauiJTrKuLlTSqHFVnUq8KcP1nD1pEVs3ns41KGJRJzsRh7N3H11pgeECY08ih53Z+qyJB6ZtZbjp5O5u3tDhnWpR0x0jgtJixR5eR55RELikKLJzLiyTS0++b9u/KpRVZ78eAOXv7iQNTs1ChUJhIj+Z5imraRquZJMvKE1E65rxa4DJ+j/wkKemrOe46dUaFEkP7KctsrwA2YVgFruvjI4IQWepq0E4JejJ3n0w3VMXZZE/SpleGJQc9rEVQx1WCKFVn4WzM928IWZlTezisAK4DUzGxfIIEWC7ZzSxfnblS34563tOH4qmStfWsRDM9Zw5IQKLYrkVk6nrWLd/SAwEHjN3VsD3YMXVmBo2koy0vX8Kswd1ZWbOsTx+qKUQosJG1U8UyQ3cpo8YsysGimbL80KYjwB5e4z3X1YbGxsqEORQqZMiRge6t+Ud3/bgRLForjx1W/4/bsr+OXoyVCHJhIWcpo8/kJKOZFEd19iZvWATcELS6RgtImryEd3deGOS+rz3rc/0H1cArNX/RjqsEQKvVwvmIcjLZhLTqzZeYB7p65kzc6D9G56Hn8Z0JSq5UuGOiyRkMn3grlIUdC0eizv39GJ+3o35rMNe+g+bj7vLt2hQosiGYjo5KEFc8mtYtFR3H5xfWaP7ELj88pzz9SV3PjqN+z4+WioQxMpVDRtJZKJ5GTnX4u38fjs9Thwb69G3NghToUWpcjI17SVmTU2s1+bWdl07b0DFaBIYRQVZdzQIY65o7vRNq4iD81cy5UvLSJxz6FQhyYSctltQ3sXKft53AmsNrMBad7+azADEyksapxTism3tGXcVS3YvPcwlz67gBc+28SpM8mhDk0kZLIbedwGtHb3y4GLgQfNbGTqexq7S5FhZgxsVZN5o7rRo+m5/G3uRvq/sJDVP2g9TYqm7JJHtLsfBnD3raQkkD6ppUkKffLQgrkEWpVyJRh/bSteuqE1+w6fYMD4hTw+W4UWpejJLnnsMrMLz75ITSR9gcpAfDADCwQ9YS7B0qvpecwb1Y3BrWoycf5mLn32S775/udQhyVSYLJLHjcCu9I2uPtpd78R6Bq0qETCQGzpYjwxuDlvDm3PqeRkrnppEQ++v5rDKrQoRUB2m0ElufsuSCnFbmbNzayVmbUCjhVIhCKFXOeGlZlzd1du6RTHm4u30XPcfD7fsCfUYYkEVY6e8zCzR4Cbgc3A2Q+4u/8qeKEFjp7zkIKybNt+7pu2ksQ9hxnYsgYP9m1ChTLFQx2WSJ5k9ZxHTpPHBiDe3cOy5KiShxSkE6fPMP6zRF78YjPnlC7Gw/2bcWn8eZgV+ntMRP5LIGpbrQbOCVxIIpGrREw0o3s2YuadnakWW4o73lrOb99Yxu6Dx0MdmkjA5HTk0YaUhwVXAyfOtrt7/+CFFjgaeUionD6TzCsLvmfcvI0Uj4niwcuacGWbmhqFSFgIxLTVGuAlYBXwn8dq3X1+oIIMBjPrB/Rr0KDBbZs2afsRCZ3vfzrCfdNW8s33P9O5QWXGDoynVsXSoQ5LJEuBSB7z3b1bwCMrIBp5SGGQnOy89c12Hp+9njPJzj29GnFTxziiVWhRCqlArHksM7OxZtbh7K26qbfrikgORUUZ119Uh7mjunJRvYr8ZdZaBk/8ik27VWhRwk9ORx6fZ9CsW3VF8sjdmbFiJw/NWMORE2cY8asGDO9Wn+IxEb3FjoSZfE9bhTslDymsfjp8godnrmXmip00Pq8cTwxqTotaurFRCod8T1uZ2V/N7Jw0ryuY2aOBClCkqKpctgTPD2nJyze2Yf/Rk1zx4kLGfrSOYydVaFEKt5yOkfu4+y9nX7j7fuDS4IQkUvT0aHIu80Z34+q2tXgpYQt9nk3g6y37Qh2WSKZymjyizazE2RdmVgookcXxIpJL5UsWY+zA5rz1m/YkO1wz6WseeG8Vh46fCnVoIv8jp8njTeBTMxtqZrcC84DXgxeWSNHVsUFlPr67C7/pXJcp32yn59MJfLZ+d6jDEvkvOV4wT92zvDspm0DNdfc5wQwskLRgLuHq2+37uXfqSjbtOczlF1bnT/2aUlGFFqWA5PluKzMzzya75OSY/DCzesADQKy7D07TXgZIAP7s7rOy6kPJQ8LZydPJjP88kRe/SKRcyWI81L8p/ZpXU4kTCbr83G31uZndaWa103VY3Mx+ZWavAzdlceJXzWyPma1O197bzDaYWaKZjckqAHff4u5DM3jrPuCdbOIXCXvFY6IY1eN8Zt7ZmVoVSnHXlG+57Z9L2XVAhRYldLJLHr2BM8AUM9tpZmvN7HtgEzAEeNrdJ2fx+cmpffyHmUUD44E+QBNgiJk1MbN4M5uV7qtqRp2aWXdgLaCJYCkyGp9Xnum/68QDl17AgsSf6DFuPlO+2U5ReFZLCp+YrN509+PAi8CLZlaMlL3Lj6W9bTebzyeYWVy65nZAortvATCzt4EB7j6WlP3Rc+ISoAwpyeeYmX3k7snZfEYk7EVHGbd1rUfPpucyZtoq7p++ihnf7eTxQfHUqVQm1OFJEZLjWgjufsrdf8xp4shCDWBHmtdJqW0ZMrNKZjYRaGlm96fG8oC73w28BbycUeIws2FmttTMlu7duzefIYsULnUqleGt29ozdmA8q384QK9nEvjHl1s4k6xRiBSMUBTSyWiVL9Mr3t33uftwd6+fOjpJ+97kzBbL3X2Su7dx9zZVqlTJZ8gihY+ZMaRdbeaN7kbnBpV59MN1DJzwFRt2qdCiBF8okkcSUCvN65rAzmCcyMz6mdmkAwcOBKN7kULhvNiSvHxjG54b0pIdPx+l7/Nf8vS8jZw8rZlcCZ4sk4eZzTGzUWbWOIDnXAI0NLO6ZlYcuAaYEcD+/8PdZ7r7sNjY2GB0L1JomBn9W1Tnk9HduDS+Gs9+uol+zy/gux35nWUWyVh2I4+bgP3AQ2a23MwmmNkAMyubk87NbAqwCGhkZklmNtTdTwMjgDnAOuAdd1+Tj58hq/Nr5CFFSsUyxXn2mpa8enMbDh4/xcAXF/LorLUcPXk61KFJhMnNE+ZRQHtSbrH9NXCMlCfNnwxeeIGhhwSlKDp0/BRjZ6/nrcXbqV2xNI8PjKdjg8qhDkvCSCB2EsTdk919kbv/yd07kTLd9EOgghSRwCpXshh/vSKet4ddRJTBtf9YzJhpKzmoQosSABG9GZSZ9QP6NWjQ4LZNmzaFOhyRkDl28gzPfLKRl7/cQpVyJXjs8ni6Nzk31GFJIaedBDVtJQLAih2/cN+0lazfdYh+Larz535NqFxWuytIxgIybSUi4a9FrXOYMaIzo3ucz8erf6THuPm8922SSpxIruU5eZjZLYEMJBh0t5XI/yoeE8Vdv27Ih3d1Ia5yGUb9ewW3Tl7Czl+OhTo0CSN5nrYys+3uXjv7I0NP01YiGTuT7Lz+1VaemrOB6ChjTJ/GXNuuNlFRKvcu+dvPY2VmbwHnu3tYTJYqeYhkbfu+o9z/3koWJu6jXd2KPDGoOXUrq9BiUZef5LEb6EXKg4L/9RbwlbtXD1iUQaC7rURyzt15Z+kOHv1wHSdPJzO6x/kM7VyXmGgtjRZV+VkwnwWUdfdt6b62Al8EOM6AU3kSkZwzM65uW5tPRnej2/lVGDt7PQMnfMXanQdDHZoUQrpVV0T+h7vz0apd/HnGan45eorbL67PiF81oERMdKhDkwKkW3VFJFfMjMuaV2PeqG70v7A6z3+WyGXPLWDZtvQz2FJUZVdVd3l2HeTkGBEJTxXKFGfcVRfy2i1tOXriNIMnfsXDM9dw5IQKLRZ12S2YHyNlv/JMDwFiC+stu1owFwmcwydO8+TH6/nnom3UrFCKsQPj6dJQG61FsvzcbVUnB/2fcfekvAZXELTmIRI433z/M2OmrWTLT0e4qk1NHri0CbGli4U6LAkC1bZS8hAJqOOnzvDcp5t4KWELFcsU55EBzejd7LxQhyUBpgVzEQmoksWiubd3Yz64oxNVypZg+JvLuONfy9l76ESoQ5MCouQhInnWrEYsH4zoxD29GjFv7W66j5vPtGUqtFgU5Ch5mFmTDNouDng0AabCiCLBVyw6ijsuacBHIzvToGpZ/u/dFdz82hJ+UKHFiJajNQ8zWw28ATwJlEz9s427dwhueIGhNQ+RgpGc7Pxz0VaenLMBA+7r05jr29dRocUwFYg1j/ZALeArYAmwE+gUmPBEJFJERRk3d6rLnLu70qpOBf70wRqunrSIzXsPhzo0CbCcJo9TwDGgFCkjj+/dPTloUYlIWKtVsTT/vLUdf7uyBRt3H6bPs18y/vNETp3Rr41IkdPksYSU5NEW6AwMMbOpQYtKRMKemTG4dU3mje7KrxtX5ak5G7h8/EJW/6A1yEiQ0+Qx1N3/5O6n3H2Xuw8APghmYCISGaqWK8mE61sz4bpW7D54ggHjF/LUnPUcP3Um1KFJPsTk8Lg9Zpa+BMn8QAcjIpGrT3w1OtSvxKMfrmP855uZvXoXTw5qTpu4iqEOTfIgp3dbrQKclFpWJYG6wAZ3bxrc8AJDd1uJFC4JG/dy//RV7DxwjJs6xHFPr0aUKZHTf8tKQcn33VbuHu/uzVP/bAi0AxYEMshg0HMeIoVT1/OrMHdUV27qEMfri7bS8+kEEjbuDXVYkgt5rm1lZsvdvVWA4wkKjTxECq+lW3/m3mkr2bL3CINb1+SPl13AOaWLhzosIeuRR47GiWY2Os3LKKAVoH8miEi+tYmryEd3deGFzxKZMH8zX2zYyyMDmtInvlqoQ5Ms5PRuq3JpvkoAHwIDghWUiBQtJYtF8/tejZgxohPnli/B7f9azvA3lrHn4PFQhyaZUEl2ESlUTp9J5uUvv+fpTzZSMiaKB/s2YXDrmpipxElBy89mUDNJucsqQ+7eP//hBZ+Sh0j42bz3MGOmrWTJ1v10aViZv14RT62KpUMdVpGSn+TRLauO3T0snvVQ8hAJT8nJzpuLt/HE7PU4cG+vRtzYIU6FFgtIfpJHbXffHrTICoiSh0h4++GXY/xh+irmb9xL6zoVeGJQPA2qlgt1WBEvP895vJ+mk2kBjUpEJIdqnFOKybe0ZdxVLdi89zCXPruAFz7bpEKLIZRd8kg7NqwXzEBERLJiZgxsVZN5o7rRo+m5/G3uRvq/oEKLoZJd8vBMvi8wZlbPzF5JW8XXzC42sy/NbGI47GgoIoFTpVwJxl/bikk3tGbf4ZRCi4/PVqHFgpZd8mhhZgfN7BDQPPX7g2Z2yMwOZte5mb1qZntSdyJM297bzDaYWaKZjcmqD3ff4u5D0zcDh0mps5WUXRwiEnl6Nj2PeaO7MbhVTSbO30yfZ79k8ZZ9oQ6ryAjqcx5m1pWUX/L/dPdmqW3RwEagBym/+JcAQ4BoYGy6Lm519z2pn5vq7oNTv49y92QzOxcY5+7XZRWHFsxFItvCxJ8YM30lO34+xg0X1eHe3o0oV7JYqMMKe/kuT5JX7p5gZnHpmtsBie6+JTW4t4EB7j4W6JvDfs+uku0n5Yl3ESnCOjWozJy7u/L3uRt5deH3fLpuN49dEc8ljauGOrSIldPyJIFUA9iR5nVSaluGzKySmU0EWprZ/altA83sJeAN4IVMPjfMzJaa2dK9e1WGSyTSlS4ew4N9mzDt9o6UKRHDLZOXMOrf3/HzkZOhDi0ihaKAfkZP92T1FPs+YHi6tunA9KxO4u6TgEmQMm2V+zBFJBy1ql2BWXd1Zvznm3nx80QSNu7l4QFNuSy+mkqcBFAoRh5JQK00r2sCO4NxIu3nIVI0lYiJZnSP85l5Z2dqVCjFiLe+5bdvLGO3Ci0GTCiSxxKgoZnVNbPiwDXAjGCcyN1nuvuw2NjYYHQvIoXcBdXKM/32jvzh0sbM37iX7uPm8+8l2ykKBWGDLajJw8ymAIuARmaWZGZD3f00MAKYA6wD3nH3NUE6v0YeIkVcTHQUw7rW5+O7u9KkWnnum7aK619ZzPZ9R0MdWlhTSXYRKTKSk50pS7Yz9qP1nEl2ft+rETd3jCNahRYzlO89zEVEIsjGN3sAAAu/SURBVEFUlHFd+zrMG92VDvUr8cistQya8BUbdx8KdWhhJ6KTh6atRCQj1WJL8cpNbXj2mgvZtu8Ilz33Jc9+somTp1VoMac0bSUiRdq+wyd4aOZaZq7YSePzyvHEoOa0qHVOqMMqFDRtJSKSiUplS/D8kJa8fGMb9h89yRUvLuSvH63j2EkVWsxKRCcPTVuJSE71aHIu80Z34+q2tZmUsIU+zyawaLMKLWYmopOHnvMQkdwoX7IYYwfG89Zt7XFgyMtf84f3VnHw+KlQh1boRHTyEBHJi471K/PxyK4M61qPt7/ZTs9xCXy6bneowypUIjp5aNpKRPKqVPFo/nDpBUz/XSdiSxVj6OtLGfn2t+w7fCLUoRUKuttKRCQbJ08nM+GLzbzw+SbKlSzGn/s1oX+L6hFfaFF3W4mI5EPxmChGdm/IrDu7UKtiaUa+/R23/XMpuw4U3UKLSh4iIjnU6LxyTL+9I3+87AIWJP5Ej3HzmfJN0Sy0GNHJQ2seIhJo0VHGb7rUY87dXWlWI5b7p6/i2pcXs/WnI6EOrUBpzUNEJI/cnbeX7OCvH67jVHIy/9ejEbd2rhsxhRa15iEiEgRmxpB2tZk3uhudG1ThsY/WMfDFhazfdTDUoQWdkoeISD6dF1uSl29szfNDWpK0/xh9n1vAuHkbOXE6ckucKHmIiASAmdGvRXXmje5GvxbVee7TTfR7fgHfbt8f6tCCQslDRCSAKpYpztNXX8hrN7fl0PHTDJzwFY/MWsvRk6dDHVpARXTy0N1WIhIqlzSuytxRXbm2XW1eWfA9vZ/5kq8Sfwp1WAET0clDhRFFJJTKlSzGY1fE8/awi4iOMq79x2LGTFvJgWPhX2gxopOHiEhhcFG9Sswe2YXfdqvHO0t30GPcfOau2RXqsPJFyUNEpACULBbN/X0u4P07OlGxTHGGvbGMEW8t56cwLbSo5CEiUoCa1zyHmXd25vc9z2fumt10Hzef975NCrsSJ0oeIiIFrFh0FCN+1ZAP7+pM3cplGPXvFdwyeQk//HIs1KHlmJKHiEiINDy3HFOHd+RPfZuweMvP9Bw3nze+3kZycuEfhSh5iIiEUHSUcWvnusy5uysX1j6HB99fzTUvf833hbzQYkQnDz3nISLhonal0rw5tD1PDmrOuh8P0vuZBCbO38zpM8mhDi1DqqorIlLI7D54nAffX83ctbuJrxHLk4Obc0G18gUeh6rqioiEkXPLl+SlG1rz4nWt+PHAMfo9v4C/z91QqAotKnmIiBRCZsal8dWYN6ob/S+szvOfJXLZcwtYtq1wFFpU8hARKcQqlCnOuKsuZPItbTl28gyDJ37FwzPXcOREaAstKnmIiISBixtVZc6ortxwUR1eW7iVXs8ksGBT6AotKnmIiISJsiVi+MuAZrzz2w4Uj47i+lcWc+/UFRw4WvCFFpU8RETCTLu6FfloZBduv7g+05b/QPen5/Px6oIttKjkISIShkoWi+a+3o354I5OVClbguFvLuOOfy1n76GCKbSo5CEiEsaa1YjlgxGduKdXI+atTSm0OG1Z8AstFvrkYWb1zOwVM5uapi3KzB4zs+fN7KZQxiciEmrFoqO445IGfDSyCw2qluX/3l3BTa8tIWn/0aCdM6jJw8xeNbM9ZrY6XXtvM9tgZolmNiarPtx9i7sPTdc8AKgBnAKSAhu1iEh4alC1LO/+tgMP92/K0q0/0/PpBP65aGtQzhXskcdkoHfaBjOLBsYDfYAmwBAza2Jm8WY2K91X1Uz6bQQscvfRwO1BjF9EJKxERRk3dYxj7qiutImryIZdh4Jynpig9JrK3RPMLC5dczsg0d23AJjZ28AAdx8L9M1h10nAydTvC8/z+iIihUTNCqV5/Za2nAxSYcVQrHnUAHakeZ2U2pYhM6tkZhOBlmZ2f2rzdKCXmT0PJGTyuWFmttTMlu7duzdAoYuIhA8zo0RMdFD6DurIIxOWQVumtwW4+z5geLq2o0D6dZD0n5sETIKUqrq5D1NERDITipFHElArzeuawM5gnEj7eYiIBEcokscSoKGZ1TWz4sA1wIxgnMjdZ7r7sNjY2GB0LyJSZAX7Vt0pwCKgkZklmdlQdz8NjADmAOuAd9x9TZDOr5GHiEgQaCdBERHJUJHdSVAjDxGR4Ijo5KE1DxGR4CgS01Zmthf4BchqCBKbxfuVgdDtupI3Wf08hflc+ekrt5/N6fE5OS67YyLt+oKCu8Z0fYXu+qrj7lUyfMfdi8QXMCmv7wNLQx1/oH/ewnqu/PSV28/m9PicHFfUrq9A/38vqPPo+grcV0RPW6UzM5/vh5uC/HkCea789JXbz+b0+JwcV9SuLyi4n0nXVyG8vorEtFV+mdlSz+SOA5H80vUlwRSs66sojTzyY1KoA5CIputLgiko15dGHiIikmsaeYiISK4peYiISK4peYiISK4peeSSmZUxs9fN7GUzuy7U8UhkMbN6ZvaKmU0NdSwSmczs8tTfXx+YWc+89qPkAZjZq2a2x8xWp2vvbWYbzCzRzMakNg8Eprr7bUD/Ag9Wwk5uri933+LuWW50JpJeLq+x91N/f90MXJ3Xcyp5pJgM9E7bYGbRwHigD9AEGGJmTUjZvOrsNrraP11yYjI5v75E8mIyub/G/pj6fp4oeQDungD8nK65HZCY+i/Bk8DbwABSdkKsmXqM/vtJtnJ5fYnkWm6uMUvxBDDb3Zfn9Zz65Ze5Gvz/EQakJI0awHRgkJlNIDJLTkjByPD6MrNKZjYRaGlm94cmNIkQmf0OuxPoDgw2s+F57Twmf7FFNMugzd39CHBLQQcjESez62sfkOe/0CJpZHaNPQc8l9/ONfLIXBJQK83rmsDOEMUikUfXlwRbUK8xJY/MLQEamlldMysOXAPMCHFMEjl0fUmwBfUaU/IAzGwKsAhoZGZJZjbU3U8DI4A5wDrgHXdfE8o4JTzp+pJgC8U1psKIIiKSaxp5iIhIril5iIhIril5iIhIril5iIhIril5iIhIril5iIhIril5iABmdsbMvkvzNSb7TwVfmriqm9ni1O+3m9neNLHGpfvMxWa2KF1bjJntNrNqZvaUme0ys98X5M8ikUW1rURSHHP3CwPZoZnFpD6olR9p42qf2u/NQBt3H5HJZxKAmmYW5+5bU9u6A6vd/UfgHjM7ks+4pIjTyEMkC2a21cweNrPlZrbKzBqntpdJ3YBniZl9a2YDUttvNrN3zWwmMNfMSpvZO2a20sz+nTp6aGNmQ83s6TTnuc3MxuUhvvpm9rGZLTOzL82ssbsnA+/y3xv9XANMydd/DJE0lDxEUpRKN22V9hfvT+7eCpgAnJ3qeQD4zN3bApcAT5lZmdT3OgA3ufuvgN8B+929OfAI0Dr1mLeB/mZWLPX1LcBreYh7EnCnu7dOje3F1PYppCQMzKwEcCkwLQ/9i2RI01YiKbKatpqe+ucyUrYhBuhJyi//s8mkJFA79ft57n52Y57OwLMA7r7azFamfn/EzD4D+prZOqCYu6/KTcBmVhboCLxr9p/q2yVS+19iZmXNrBFwAfC1u+/PTf8iWVHyEMneidQ/z/D//84YMMjdN6Q90MzaA2nXEzLaU+GsfwB/ANaTt1FHFPBLFknvbVJGHxegKSsJME1bieTNHOBOS/0nv5m1zOS4BcBVqcc0AeLPvuHui0nZb+Fa8vDL3d0PAt+b2ZWp/ZuZtUhzyBTgeuBXqNy7BJiSh0iK9Gsej2dz/CNAMWClma1OfZ2RF4EqqdNV9wErgQNp3n8HWJiPKaXrgKFmtgJYQ5p90N19LXCUlLUZ3V0lAaWS7CJBZGbRpKxnHDez+sCnwPnufjL1/VnA0+7+aSafP+zuZYMQ10PAYXf/W6D7lqJBIw+R4CoNLEgdGbwH3O7uJ83sHDPbSMpCfYaJI9XBsw8JBiogM3uKlOksjUYkzzTyEBGRXNPIQ0REck3JQ0REck3JQ0REck3JQ0REck3JQ0REck3JQ0REcu3/ASBNICJ7CsMVAAAAAElFTkSuQmCC\n",
"text/plain": [
"