\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.17?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/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 and SkyDiffuseCube](#SkyModel-and-SkyDiffuseCube)\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 :ref:`model-gallery`.\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"
]
},
{
"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.000e+00 nan nan False 0.000e+00\n",
"amplitude 1.000e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 1.000e+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/stable/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\n",
" CompoundSpectralModel : CompoundSpectralModel\n",
" PowerLawSpectralModel : PowerLawSpectralModel\n",
" PowerLaw2SpectralModel : PowerLaw2SpectralModel\n",
" SmoothBrokenPowerLawSpectralModel : SmoothBrokenPowerLawSpectralModel\n",
" ExpCutoffPowerLawSpectralModel : ExpCutoffPowerLawSpectralModel\n",
" ExpCutoffPowerLaw3FGLSpectralModel : ExpCutoffPowerLaw3FGLSpectralModel\n",
" SuperExpCutoffPowerLaw3FGLSpectralModel: SuperExpCutoffPowerLaw3FGLSpectralModel\n",
" SuperExpCutoffPowerLaw4FGLSpectralModel: SuperExpCutoffPowerLaw4FGLSpectralModel\n",
" LogParabolaSpectralModel : LogParabolaSpectralModel\n",
" TemplateSpectralModel : TemplateSpectralModel\n",
" GaussianSpectralModel : GaussianSpectralModel\n",
" AbsorbedSpectralModel : AbsorbedSpectralModel\n",
" NaimaSpectralModel : NaimaSpectralModel\n",
" ScaleSpectralModel : ScaleSpectralModel\n",
"\n"
]
}
],
"source": [
"from gammapy.modeling.models import SPECTRAL_MODELS\n",
"\n",
"print(SPECTRAL_MODELS)"
]
},
{
"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.200e+00 nan nan False 0.000e+00\n",
"amplitude 2.700e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 1.000e+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(emin=1 * u.TeV, emax=10 * u.TeV)\n",
"print(flux)\n",
"\n",
"eflux = pwl.energy_flux(emin=1 * u.TeV, emax=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(emin=energy[:-1], emax=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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xV9f3H8dcnYa+wQVkBUTAQkK1MqUxlWMC9pSJWFKHOqq1WLY4WFUURq+KoOBjKEAEXS0SWsrcoW4aA7JHP74/EX9M0O/fmJDfv5+NxH3q/595z3mmP+eR7xueYuyMiIpIVUUEHEBGR/EfFQ0REskzFQ0REskzFQ0REskzFQ0REskzFQ0REsqxQ0AFyQ8WKFT02NjboGCIi+crixYv3uHul1JYViOIRGxvLokWLgo4hIpKvmNmPaS2L6MNWZtbTzEYfOHAg6CgiIhEloouHu0929wExMTFBRxERiSgRXTxERCQ8VDxERCTLVDxERCTLIrp46IS5iEh4RHTxyOkJc3dnwpKtHD91OsTJRETyt4guHjm1+MdfGPrB9/QYMZclP/0SdBwRkTxDxSMdzWPL88aNLTh0/BR9X/6ax6as4siJU0HHEhEJnIpHBjrWr8yMIe25plVNXpv7A12fm828DXuCjiUiEigVj0woXawwj18az/sDzqdQVBTX/GsB941bxoGjJ4OOJiISCBWPLGhVpwLTBrdjYIezGLdkK52Hz2LGyp1BxxIRyXUqHllUrHA093evz0d/bEOFUkUZ8PZibn93CXsOHQ86mohIrono4hHO+zziq8cwaVAb7u5yDjNX7qLT8FlMXLoVdw/5tkRE8horCL/smjdv7uFsyb7h51+5d9wylvy0nwvrVeKJ38dTrWzxsG1PRCQ3mNlid2+e2rKInnnklrqVS/PhwNb8tWccCzbto8vwWbw9fzMJCZFfmEWkYFLxCJHoKOOmNrWZMaQ9TWqW4+GPV3Ll6G/YtPtQ0NFEREJOxSPEapQvwdv9W/J0v0as2XmQbs/P4eWvNnLqdELQ0UREQkbFIwzMjMub1+CzoR3oWK8ST326hktfmsfK7WrQKCKRQcUjjCqXKcaoa5vx0jVN2XngGL1enMcz09dw7KQaLYpI/qbiEWZmxsXxZzBzSAcuPa8aI7/cyCUj5rD4x31BRxMRyTYVj1xSrmQR/nl5Y968uSXHTibQb9R8Hpm0ksPH1WhRRPIfFY9c1uGcSkwf0p7rz6/Fm/M30+XZ2cxetzvoWCIiWZIvioeZ1TGz18xsXLKxS83sVTP72My6BJkvq0oVLcSjvRvywa0XULRwFNe//i13f/g9+4+cCDqaiEimhL14mNnrZvazma1IMd7NzNaa2QYzuz+9dbj7Jnfvn2LsI3e/BbgRuCLkwXNBi9jyfHJnO27veBYTl26j87Oz+XTFjqBjiYhkKDdmHmOAbskHzCwaGAl0B+KAq8wszszizWxKilflDNb/UNK68qVihaO5p2t9Pr69DZVKFWXgO0u47Z3F/PzrsaCjiYikqVC4N+Dus80sNsVwS2CDu28CMLP3gN7uPgzokZn1mpkBTwLT3H1J6BIHo2G1GD4e1IbRszfx/Ofr+XrjXh7uEUffptVI/FFFRPKOoM55VAO2JHu/NWksVWZWwcxGAU3M7IGk4TuATkA/MxuYyncGmNkiM1u0e3f+OCFdODqK2zvWZdrgdpxTpRR3f/g917/+LVv2HQk6mojIf8mVrrpJM48p7t4w6f1lQFd3/0PS++uAlu5+Rzi2H+6uuuGQkOD8e8GPPDltDQ7c27Ue118QS1SUZiEikjvyYlfdrUCNZO+rA9tDvZFwPs8j3KKijOsuiGX6kPa0iC3PI5NXcfkr89nwsxotikjwgioeC4Gzzay2mRUBrgQmhXoj7j7Z3QfExMSEetW5pnq5Eoy5qQX/vKwxG3Yf4uLn5zDyyw2cVKNFEQlQblyqOxaYD9Qzs61m1t/dTwGDgOnAauADd18Zhm3n25lHcmZG32bVmTmkA53iKvPM9LX0fnEeK7bl759LRPIvPUkwH/p0xU4e/ngF+w6fYED7Ogy+6GyKFY4OOpaIRJi8eM5DcqBbw6p8NqQDfZtW4+WvNnLx83NYuFmNFkUk90R08YiUw1apiSlRmKf7Nead/q04cTqBy0bN5y8fr+CQGi2KSC7QYasIcOTEKZ6ZvpYxX2/mzJjiPPH7hlxYL6Mb80VE0qfDVhGuRJFC/LVnA8YNbE3xItHc+MZChn7wHb8cVqNFEQmPiC4ekXzYKjXNapVj6p1tGdSxLpO+207nZ2cxddkOCsLsUkRylw5bRahV2w9y3/hlLN92gC5xVXj80oZULlMs6Fgiko/osFUBFHdmGSb+sTX3d6/PrHW7uWj4LD5YuEWzEBEJiYguHgXtsFVKhaKjGNjhLKYNbse5Z5Th3vHLuO41NVoUkZzTYasCIiHBeffbnxj2yWoSHO7pWo8bWscSrUaLIpIGHbYSoqKMa8+vxYyhHWhVpzx/m7KKfqO+Zv2uX4OOJiL5kIpHAVOtbHHeuLEFz11xHpv3HOaSEXMZ8fl6TpxSo0URyTwVjwLIzLi0STVmDu1A14ZVGT5zHb1enMuyrfuDjiYi+UREF4+CfsI8IxVLFeWFq5rw6vXN+eXICS4dOY9hn6zm2MnTQUcTkTxOJ8wFgANHT/LktNWM/XYLsRVK8GTfRpxfp0LQsUQkQDphLhmKKV6YYX0a8e4fWpHgcOXob3hw4nJ+PXYy6GgikgepeMh/aV23Ip/e1Y4/tK3N2G9/osuzs/liza6gY4lIHqPiIf+jRJFCPNQjjvG3taZ0sULcPGYRd723lH1qtCgiSSK6eOiEec40qVmOKXe0Y/BFZzNl2Q46DZ/FpO+3q8WJiOiEuWTOmp0HuW/cMr7feoBO5yY2Wqwao0aLIpFMJ8wlx+pXLcP421rz4MXnMnfDbjo/O4v3vv1JsxCRAkrFQzKtUHQUt7Svw6eD29PgzDLcP2E5V7+6gB/3Hg46mojkMhUPybLYiiUZe8v5DOsTz4ptB+j63Gz+NWcTpxM0CxEpKAqlt9DMlmRiHbvdvWuI8kg+YWZc1bImHetV5sGJy3l86momL9vB030bUa9q6aDjiUiYpVs8gKJAr3SWGzAhdHEkv6kaU4x/3dCcyct28MiklfR4YQ63d6zLHy+sS5FCmtiKRKqMisft7r4xvQ+Y2Z0hzCP5kJnRq/GZtK1bkUcnr+S5z9YzbflOnurXiPNqlA06noiEQbp/Grr7VxmtIDOfCYru88hd5UsW4fkrm/DaDc05cPQkfV6axxNTV3H0hBotikSaDO/zMLMWwLVAO+AM4CiwApgKvOvuef5pQrrPI/f9euwkT05bw78X/ETN8iV4sm88rc+qGHQsEcmCbN/nYWZTgEHALOBSoDbQFHgcKAtMNbMeoY0rkaB0scI88ft4xt5yPmZw9asLeGDCcg6q0aJIREh35mFmVdw93a54ZlbZ3X8OebIQ0swjWEdPnObZz9bxrzmbqFS6KI9fGk/nuCpBxxKRDOTkDvMHzKxleh/I64VDgle8SDR/vvhcJv6xDeVKFOGWtxZxx9il7D10POhoIpJNGRWPLcBIM9toZk+YWcPcCCWRqXGNskwa1JYhnc7h0xWJjRY//m6bWpyI5EMZXW31T3dvAXQBjgBjzWyFmf3ZzOrkSkKJKEUKRTG409lMvbMdtSqUZPB739H/zUVs33806GgikgVZ7qprZs2AfwGN3D06LKlCTOc88qbTCc6Yrzfzj+lriY4yHri4Ple1qElUlAUdTUQIQVddM4s2s+5m9iaJl+huAq4IYUYpgKKjjP5tazP9rvY0rhHDgxNXcNWr3/DDHjVaFMnrMrpUt6OZjQa2AXcCXwBnu3tfdx+XGwEl8tWsUIJ3+rfiqb7xrNpxkG7PzeaVWRs5dToh6GgikoaMZh5/A5YC8e7e3d3fzO2bAs2sjpm9Zmbjko2da2ajzGycmd2Wm3kkPMyMK1rU5LOhHWh/TiWGTVtDn5e/ZvWOg0FHE5FUZHTCvJ27v+zuu83sfDO7HsDMKphZzYxWbmavm9nPZrYixXg3M1trZhvM7P4MMmxy9/4pxla7+0DgciDV43GSP1UpU4zR1zXjxaubsH3/UXq+MJfhM9Zy/JRanIjkJZk95/EQ8FfgoaShYsC7mfjqGKBbinVFAyOB7kAccJWZxZlZvJlNSfGqnE6mXsBc4PPM/AySf5gZPRqdycwhHejV+ExGfLGBS0bMZfGPvwQdTUSSZLZndj/gYuAwgLtvA8pk9CV3nw3sSzHcEtiQNKM4AbwH9Hb35e7eI8UrzRsQ3X2Su7cGrsnkzyD5TLmSRRh+xXm8cVMLjhw/Rb9RX/Po5JUcOXEq6GgiBV5mi8dxT7ym1wHMrEQOtlmNxJsPf7M1aSxVSYfIRgFNzOyBpLELzWyEmb0CfJLG9waY2SIzW7R79+4cxJWgdaxXmelD2nNtq1q8MW8zXZ6dzdz1e4KOJVKgZfQ8j99MMLORQIyZ3QT0B17P5jZTu4g/zZtN3H0vMDDF2FfAV+ltxN1HA6Mh8T6PrIaUvKV0scI8dmlDejQ6g/snLOfa1xZwefPqPHhJHDHFCwcdT6TAydTMw92fAqYAk4DGwBPu/lw2t7kVqJHsfXVgezbXlS49zyPytKpTgWmD2zGww1mMX7KNzsNnMX3lzqBjiRQ4GXXVneHuXXK0AbNYYIq7N0x6XwhYB1xE4v0jC4Gr3X1lTraTHt1hHpmWbz3AveOXsXrHQS6JP4NHejWgUumiQccSiRg5ucO8Ug43PBaYD9Qzs61m1t/dT5H4jJDpwGrgg3AVDs08Ilt89RgmDWrDPV3rMXPVLjoNn8X4xVvVaFEkF2Q089gE3J3WcnefEI5QoaaZR+Tb8POv3DtuGUt+2k+Hcyrx9z7xVCtbPOhYIvlaejOPjIrHXuBj0jjJ7e43hyZieKl4FAynE5y352/m6elrMeC+7vW5tlUtNVoUyaacFI8l7t40bMnCzMx6Aj3r1q17y/r164OOI7lky74j/Hnicuas30OL2HI82bcRZ1UqFXQskXwnJ+c88vWfbO4+2d0HxMTEBB1FclGN8iV46+aWPNOvEWt3/kr35+fw0lcb1GhRJIQyKh7X5UoKkRAzMy5rXoPP/tSB39WrzNOfruXSl+axcrsunhAJhYwaI65Ib3lep6utpHLpYoy6rhkvX9OUnQeO0+vFeTwzfQ3HTqrRokhOZPlJgvmRTpgLwP4jJ3h86mrGLd7KWZVK8lTfRjSPLR90LJE8K8dPEkyxsnJm1ijnsURyV9kSRfjHZY156+aWHDuZwGWvzOeRSSs5fFyNFkWyKrMt2b8yszJmVh74HnjDzIaHN1rO6bCVpKb9OZWYMaQ9N1wQy5vzExstzl6n5pkiWZHZmUeMux8E+gBvuHszoFP4YoWGrraStJQsWohHejXgw1svoGjhKK5//Vvu/vB79h85EXQ0kXwhs8WjkJmdQeKT+6aEMY9IrmoeW55P7mzH7R3PYuLSbXQaPptpy3cEHUskz8ts8fgbib2oNrj7QjOrA+iuO4kIxQpHc0/X+kwa1IYqZYpy27+XMPDtxfx88FjQ0UTyLF1tJZLMydMJ/GvODzz72TqKFYri4R5x9GtWHbN8fb+sSLaE9Gqr/EQnzCWrCkdHcduFZzFtcDvqVy3DPeOWcf3r37Jl35Ggo4nkKZp5iKQhIcH594IfeXLaGhy4t2s9rr8gVo0WpcAosDMPkZyIijKuuyCWGUM70CK2PI9MXsVlr8xnw8+/Bh1NJHAZFg8zq29mF5lZqRTj3cIXSyTvqFa2OGNuasHwyxuzcfchLn5+Li9+sZ6TarQoBVi6xcPM7iTxeR53ACvMrHeyxX8PZzCRvMTM6NO0OjOHdKBzgyr8Y8Y6er04jxXbdD5NCqaMZh63AM3c/VLgQuBhMxuctCzPH/jVCXMJtUqlizLy6qa8cl0z9h46Tu+R83hymhotSsGT0cOgVrl7XLL3pYBxwCrgd+5+Xvgj5pxOmEs4HDhykr9/spr3F22hTsWSPNm3ES1rq9GiRI6cnDDfaWb/XyDc/RDQA6gIxIcuokj+E1OiME/1a8Q7/VtxMiGBy1+Zz8MfreCQGi1KAZBR8bge2Jl8wN1Pufv1QPuwpRLJR9qeXZHpd7XnpjaxvLPgR7oMn8WXa38OOpZIWGX0MKit7r4T/tOK3cyamllT4GiuJBTJB0oUKcRfezZg3MDWlChaiJveWMjQ97/jl8NqtCiRqVBmPmRmjwE3AhuB306SOPC78MQSyZ+a1SrH1DvbMvKLDbz01UZmr9/No70acnF8VbU4kYiSqTvMzWwtEO/u+fLPKJ0wlyCs3nGQe8ctY/m2A3SJq8JjlzakSpliQccSybRQ3GG+Aigbukgike/cM8ow8Y+teaB7fWat202n4bP4YOEWCkJLIIl8mZ15NCfxZsEVwPHfxt29V/ii5ZyZ9QR61q1b95b169VBXoLzw57D3Dd+Gd/+sI+2dSsyrE88NcqXCDqWSLrSm3lktnisBF4BlgP/35PB3WeFKmQ46bCV5AUJCc673/7Ek9PWcDrBuadrPW5oHUu0Gi1KHhWK4jHL3TuEPFkuUfGQvGT7/qM8OHE5X67dTZOaZXm6byPOrlI66Fgi/yMU5zwWm9kwM7vgt0t1ky7XFZEsOrNscV6/sQXPX3kem/cc5pIRcxnx+XpOnFKjRck/MnWpLtAk6Z/nJxvTpboi2WRm9D6vGm3qVuTRyasYPnMdnyzfwVN9G9G4hq5NkbxPD4MSyQNmrtrFQx8tZ/evx7mlXR3u6nQOxYtEBx1LCrgcH7Yys7+bWdlk78uZ2eOhCihS0HWOq8LMoR24okUNXpm9ie7Pz+abTXuDjiWSpsye8+ju7vt/e+PuvwAXhyeSSMFUplhhhvVpxLt/aEWCw5Wjv+HBicv59djJoKOJ/I/MFo9oMyv62xszKw4UTefzIpJNretW5NO72vGHtrUZ++1PdHl2Nl+s2RV0LJH/ktni8Q7wuZn1N7ObgZnAm+GLJVKwlShSiId6xDH+ttaUKlqIm8cs4q73lrJPjRYlj8hU8XD3p4HHgXOBBsBjSWNhZ2Z1zOw1MxuXYrykmS02sx65kUMkCE1qlmPqne0YfNHZTF2+g07DZzHp++1qcSKBy+gZ5v9/66u7f+rud7v7n9x9emqfSeX7r5vZz2a2IsV4NzNba2YbzOz+9DK4+yZ375/KovuAD9L7rkgkKFIoiiGdz2HyHW2pUa44d45dyi1vLWLngWNBR5MCLKOZx5dmdoeZ1Uw+aGZFzOx3ZvYmcEM63x8DdEvx3WhgJNAdiAOuMrM4M4s3sykpXpVTW6mZdSLxUbg6ECwFRv2qZZjwxzY8ePG5zN2wh87DZzH22580C5FAZHSTYDfgZmCsmdUG9gPFSSw6M4Bn3f27tL7s7rPNLDbFcEtgg7tvAjCz94De7j6MxEfcZkZHoCSJxeeomX3i7ro9VyJedJRxS/s6dGlQhfvHL+eBCcuZ9N12nuwbT60KJYOOJwVIRk8SPObuL7l7G6AWcBHQxN1rufst6RWOdFQDtiR7vzVpLFVmVsHMRgFNzOyBpFwPuvtdwLvAq6kVDjMbYGaLzGzR7t27sxFTJO+qVaEk797SimF94lmx7QBdn5vNv+Zs4nSCZiGSOzLbngR3PwnsCME2UztHkuYe7+57gYFpLBuTzvdGA6Mh8Q7zrEUUyfvMjKta1qRjvco89NFyHp+6msnLdvB030bUq6pGixJemb1UN5S2AjWSva8ObA/Hhsysp5mNPnDgQDhWL5InVI0pxqvXN2fEVU3Ysu8IPV6Yw7Mz16nRooRVEMVjIXC2mdU2syLAlcCkcGzI3Se7+4CYmJhwrF4kzzAzejU+k8+GduDi+DN4/vP19HxhLt9t2Z/xl0WyIaNLdaeb2RAzq5+dlZvZWGA+UM/MtppZf3c/BQwCpgOrgQ/cfWV21p+J7WvmIQVK+ZJFeP7KJrx+Y3MOHjtJn5fm8fiUVRw5cSroaBJh0u2qa2ZVSbziqhtwDrAA+BT43N0P5UrCEFBXXSmIfj12kmHT1vDugp+oWb4ET/aJp3XdikHHknwkx08STFpJFNCKxPszLgKOAjNy607znFDxkILsm017uX/8MjbvPcKVLWrw50vOpUyxwkHHknwgJMUjlZVWBLq6+79zEi6czKwn0LNu3bq3rF+/Pug4IoE5euI0z322jlfnbKJS6aI8cWk8neKqBB1L8riwFI/8RDMPkUTfb9nPfeOXsWbnr/RsfCZ/7RlHxVJqkC2pC8UzzEUkAjSuUZZJg9oytPM5fLpiB52Hz2Li0q1qcSJZFtHFQ1dbifyvIoWiuPOis5l6ZztiK5ZkyPvfc/OYhWzffzToaJKPZLt4mNlNoQwSDrrPQyRt51QpzbiBrflLjzi+2bSPLs/O5p1vfiRBLU4kE3Iy83g0ZClEJBDRUcbNbWsz/a72NK4Rw0MfreDKV7/hhz2Hg44meVxG93ksS2sRcI675+kzbbraSiTz3J0PFm3h8amrOXEqgaGdz6F/29oUio7oo9uSjmxfbWVmu4CuwC8pFwFfu/uZIUsZRrraSiTzdh08xsMfrWDGql00qh7Dk30aEXdmmaBjSQBycrXVFKCUu/+Y4rUZ+CrEOUUkD6hSphivXNeMkVc3Zfv+o/R6cS7/nLGW46dOBx1N8hDd5yEiafrl8Akem7qKCUu2UbdyKZ7q24hmtcoFHUtyie7zEJFsKVeyCMMvP483bmrBkeOn6Dfqax6dvJLDx9VosaDLqKvukoxWkJnPBEX3eYiERsd6lZkxtAPXnV+LN+Ztputzs5mzXk/oLMgyOmF+FEjvMiUDYty9ZqiDhZIOW4mEzrc/7OP+8cvYtOcwlzevzoMXxxFTQo0WI1F6h60yegxtZp7jobNoIgVIy9rl+WRwO0Z8vp5XZm/iy7W7eax3Q7o1rBp0NMlFOmEuItm2YtsB7h23jFU7DnJJ/Bk80qsBlUrn6du/JAt0wlxEwqJhtRg+HtSGe7rWY+aqXXQaPovxi9VosSCI6OKhE+Yi4Vc4OorbO9blk8FtqVu5FH/68HtufGMh29RoMaJlqniYWVwqYxeGPE2IqTGiSO6pW7k0H956AY/0jGPh5n10GT6Lt+ZvVqPFCJXZmccHZnafJSpuZi8Aw8IZTETyn6go48Y2iY0Wm9Yqx18+XskVo+ezcfehoKNJiGW2eLQCagBfAwuB7UCbcIUSkfytRvkSvHVzS/5xWWPW7TpE9+fnMPLLDZw8nRB0NAmRzBaPk8BRoDhQDPjB3bUXiEiazIx+zaozc2h7LqpfmWemr+XSkfNYsU3nICNBZovHQhKLRwugLXCVmY0LWyoRiRiVSxfj5Wub8fI1Tdl18Di9R87jmelrOHZSt4jlZ5ktHv3d/S/uftLdd7p7b+DjcAYTkcjSPf4MPhvant83qcbILzdy8Yg5LNq8L+hYkk2ZLR4/m1nN5C9gVjiDiUjkKVuiCP+4rDFv3dyS4ycTuOyV+TwySY0W86NM3WFuZssBJ7GXVTGgNrDW3RuEN17O6EmCInnX4eOneGb6Wt6cv5kzY4ozrE887c+pFHQsSSbbTxJMZ4VNgVvd/dachssNak8iknct2ryPe8cvY9Puw/RrVp2HLjmXsiWKBB1LCEN7EndfQuLJcxGRHGkeW55P7mzHoI51mbh0G52Gz2ba8h1Bx5IMZNRVFwAzG5rsbRTQFFAzfxEJiWKFo7m7az26x1fl3nHLuO3fS+jWoCp/692AymWKBR1PUpHZmUfpZK+iwFSgd7hCiUjB1ODMGD6+vQ33davPF2t/ptPwWXy4aIsaLeZBaskuInnSxt2HuH/8MhZu/oV2Z1fk77+Pp0b5EkHHKlCyfcLczCaTeJVVqty9V87jhZ+Kh0j+lJDgvLPgR56atgYH7u1aj+sviCUqyoKOViDkpHh0SG/F7p4v7vVQ8RDJ37btP8qfJyxn1rrdNKtVjqf6xlO3cumgY0W8nBSPmu7+U9iS5RIVD5H8z92ZuHQbf5uyiiPHT3PnRXW5tcNZFI6O6McSBSonl+p+lGwl40OaSkQkC8yMPk2rM3NIBzo3qMI/Zqyj14tqtBiUjIpH8gOLdcIZJM0AZnXM7LXkjRjN7EIzm2Nmo/LDQ6lEJHQqlS7KyKubMvq6Zuw9lNho8clparSY2zIqHp7Gv2eKmb1uZj+b2YoU493MbK2ZbTCz+9MN4L7J3funkusQia1StmY1l4jkf10aVGXm0A70a1qdUbM20v35OSzYtDfoWAVGRuc8TgOHSZyBFAeO/LYIcHcvk+7KzdqT+Ev+LXdvmDQWDawDOpP4i38hcBUQzf8+nfBmd/856Xvj3L1f0r9HuXuCmVUBhrv7Nenl0DkPkcg2b8Me7p+wjC37jnLd+bW4t1s9ShcrHHSsfC+9cx7p3mHu7tE52bC7zzaz2BTDLYEN7r4pKdx7QG93Hwb0yOR6f3sQ1S8k3rQoIgVYm7oVmX5Xe/45Yx2vz/uBz1fv4onfx9OxfuWgo0WsIC5TqAZsSfZ+a9JYqsysgpmNApqY2QNJY33M7BXgbeDFNL43wMwWmdmi3bvVSUUk0pUoUoiHe8Qx/rbWlCxaiJvGLGTI+9+x7/CJoKNFpEz1tgqx1O7uSe9GxL3AwBRjE4AJ6W3E3UcDoyHxsFXWY4pIftS0Zjmm3NmWkV9u5KUvNzB73W4e7d2AS+LPwEw3F4ZKEDOPrUCNZO+rA9vDsSEz62lmow8c0KV8IgVJ0ULRDO18DpPvaEu1csUZ9O5Sbn17MbsOHgs6WsQIongsBM42s9pmVgS4EpgUjg25+2R3HxATExOO1YtIHnfuGWWYcFtr/nxxfWat202n4bN4f+FParQYAmEtHmY2FpgP1DOzrWbW391PAYOA6cBq4AN3Xxmm7WvmIVLAFYqOYkD7s/j0rvbEnVGG+8Yv59rXFvDT3iMZf1nSpK66IlJgJCQ4Yxf+xLBP1nA6wbm7az1ubB1LtBotpirkTxIUEcmPoqKMa1rVYubQ9lxwVriHatgAAAvISURBVAUem7KKvi9/zbpdvwYdLd+J6OKhw1YikpozYorz2g3Nef7K8/hx72EuGTGH5z9bz4lTCRl/WQAdthKRAm7voeM8MnkVk7/fTv2qpXmqbyMa1ygbdKw8QYetRETSUKFUUV64qgmvXt+cX46c4PcvzePvn6zm6Ak1WkxPRBcPHbYSkczqHFeFmUM7cEWLmoyevYnuz89m/kY1WkxLRBcP3echIllRplhhhvWJ591bWuHAVa9+w58nLufgsZNBR8tzIrp4iIhkR+uzKvLp4PYMaF+H9779iS7DZ/P56l1Bx8pTIrp46LCViGRX8SLR/Pnic5nwxzbEFC9M/zcXMfi9pew9dDzoaHmCrrYSEcnAiVMJvPzVRl78cj2lixXmrz3j6NX4zIhvtKirrUREcqBIoSgGdzqbKXe0o0b5Egx+7ztueWsROw8U3EaLKh4iIplUr2ppJtzWmocuOZe5G/bQefgsxn5bMBstRnTx0DkPEQm16CjjD+3qMP2u9jSsFsMDE5Zz9asL2LzncNDRcpXOeYiIZJO7897CLfx96mpOJiTwp871uLlt7YhptKhzHiIiYWBmXNWyJjOHdqBt3Uo88clq+rw0jzU7DwYdLexUPEREcqhqTDFevb4ZL1zVhK2/HKXHiLkMn7mO46cit8WJioeISAiYGT0bn8nMoR3o2fhMRny+np4vzGXpT78EHS0sVDxEREKofMkiPHvFebxxYwt+PXaKPi9/zWNTVnHkxKmgo4VURBcPXW0lIkHpWL8yM4a05+qWNXlt7g90e24OX2/YE3SskIno4qHGiCISpNLFCvPE7+N5b8D5REcZV/9rAfePX8aBo/m/0WJEFw8Rkbzg/DoVmDa4Hbd2qMMHi7bQefgsZqzcGXSsHFHxEBHJBcUKR/NA93P56PY2lC9ZhAFvL2bQu0vYk08bLap4iIjkokbVyzL5jrbc3eUcZqzcRafhs5i4dGu+a3Gi4iEikssKR0cx6HdnM/XOttSuWJIh73/PTWMWsm3/0aCjZZqKh4hIQM6uUppxA1vzlx5xLNi0jy7DZ/H2Nz+SkJD3ZyEqHiIiAYqOMm5uW5vpd7XnvJplefijFVz56jf8kMcbLUZ08dB9HiKSX9SsUIJ3+rfi6b6NWL3jIN2em82oWRs5dToh6GipUlddEZE8ZtfBYzz80QpmrNpFfLUYnu7XiHPPKJPrOdRVV0QkH6lSphivXNeMl65pyo4DR+n5wlz+OWNtnmq0qOIhIpIHmRkXx5/BzCEd6HXembzwxQYuGTGXxT/mjUaLKh4iInlYuZJFGH75eYy5qQVHT5ym36iveXTySg4fD7bRooqHiEg+cGG9ykwf0p7rzq/FG/M20/W52cxdH1yjRRUPEZF8olTRQvytd0M+uPUCikRHce1rC7h33PccOJL7jRZVPERE8pmWtcvzyeB23HbhWYxfso1Oz87i0xW522hRxUNEJB8qVjia+7rV5+Pb21CpVFEGvrOY2/+9hN2/5k6jRRUPEZF8rGG1GD4e1IZ7utZj5qrERovjF4e/0WKeLx5mVsfMXjOzccnGoszsCTN7wcxuCDKfiEjQCkdHcXvHunwyuB11K5fiTx9+zw1vLGTrL0fCts2wFg8ze93MfjazFSnGu5nZWjPbYGb3p7cOd9/k7v1TDPcGqgEnga2hTS0ikj/VrVyKD2+9gEd7NWDR5n10eXY2b83fHJZthXvmMQbolnzAzKKBkUB3IA64yszizCzezKakeFVOY731gPnuPhS4LYz5RUTylago44bWscwY0p7mseVZu/PXsGynUFjWmsTdZ5tZbIrhlsAGd98EYGbvAb3dfRjQI5Or3gqcSPr3vHO/vohIHlG9XAnevKkFJ8LUWDGIcx7VgC3J3m9NGkuVmVUws1FAEzN7IGl4AtDVzF4AZqfxvQFmtsjMFu3evTtE0UVE8g8zo2ih6LCsO6wzjzRYKmNpXhbg7nuBgSnGjgApz4Ok/N5oYDQkdtXNekwREUlLEDOPrUCNZO+rA9vDsSE9z0NEJDyCKB4LgbPNrLaZFQGuBCaFY0PuPtndB8TExIRj9SIiBVa4L9UdC8wH6pnZVjPr7+6ngEHAdGA18IG7rwzT9jXzEBEJAz1JUEREUlVgnySomYeISHhEdPHQOQ8RkfAoEIetzGw3sB9IbwoSk87yikBwT13JnvR+nry8rZysK6vfzeznM/O5jD4TafsX5N4+pv0ruP2rlrtXSnWJuxeIFzA6u8uBRUHnD/XPm1e3lZN1ZfW7mf18Zj5X0PavUP//nlvb0f4VuldEH7ZKYXIOl+c3ufnzhHJbOVlXVr+b2c9n5nMFbf+C3PuZtH/lwf2rQBy2yikzW+RpXHEgklPavyScwrV/FaSZR06MDjqARDTtXxJOYdm/NPMQEZEs08xDRESyTMVDRESyTMVDRESyTMUji8yspJm9aWavmtk1QeeRyGJmdczsNTMbF3QWiUxmdmnS76+PzaxLdtej4gGY2etm9rOZrUgx3s3M1prZBjO7P2m4DzDO3W8BeuV6WMl3srJ/ufsmd0/3QWciKWVxH/so6ffXjcAV2d2mikeiMUC35ANmFg2MBLoDccBVZhZH4sOrfnuMrp6fLpkxhszvXyLZMYas72MPJS3PFhUPwN1nA/tSDLcENiT9JXgCeA/oTeKTEKsnfUb/+0mGsrh/iWRZVvYxS/QUMM3dl2R3m/rll7Zq/GeGAYlFoxowAehrZi8TmS0nJHekun+ZWQUzGwU0MbMHgokmESKt32F3AJ2AfmY2MLsrL5SzbBHNUhlzdz8M3JTbYSTipLV/7QWy/R+0SDJp7WMjgBE5XblmHmnbCtRI9r46sD2gLBJ5tH9JuIV1H1PxSNtC4Gwzq21mRYArgUkBZ5LIof1Lwi2s+5iKB2BmY4H5QD0z22pm/d39FDAImA6sBj5w95VB5pT8SfuXhFsQ+5gaI4qISJZp5iEiIlmm4iEiIlmm4iEiIlmm4iEiIlmm4iEiIlmm4iEiIlmm4iEFnpmdNrPvkr3uz/hb4Wdmm81suZk1N7OJSdk2mNmBZFlbp/HdP5jZ2ynGqiS17S5sZu+b2T4zuzR3fhqJNLrPQwo8Mzvk7qVCvM5CSTdp5WQdm4Hm7r4n2diFwN3u3iOD75YD1gPV3f1Y0tggIN7db016/w6Jz6b5KCc5pWDSzEMkDUl/+T9qZkuSZgD1k8ZLJj18Z6GZLTWz3knjN5rZh2Y2GZhhZlFm9pKZrTSzKWb2iZn1M7OLzGxisu10NrMJOcjZwsxmmdliM5tmZlXc/Rfga+CSZB+9Ehib3e2IJKfiIQLFUxy2Sv50tT3u3hR4Gbg7aexB4At3bwF0BJ4xs5JJyy4AbnD335H41MlYIB74Q9IygC+Ac82sUtL7m4A3shPczIoCzwN93b0Z8A7wWNLisSQWDMysRlKW2dnZjkhKaskuAkfd/bw0lv02I1hMYjEA6AL0MrPfikkxoGbSv890998eytMW+NDdE4CdZvYlJPbETjofca2ZvUFiUbk+m9nPBRoAn5kZQDSJ3VQhsQneCDMrReLjRj9IyiKSYyoeIuk7nvTP0/znvxcj8S/9tck/aGatgMPJh9JZ7xskPkzsGIkFJrvnRwxY5u7tUi5w98Nm9hmJTyi8Ergtm9sQ+R86bCWSddOBOyzpT30za5LG5+aS+NTJKDOrAlz42wJ3307isxUeIvH509m1isQnELZMylLEzBokWz4WuAco6+4Lc7Adkf+i4iHyv+c8nszg848BhYFlZraC/5xjSGk8iYeQVgCvAAuAA8mW/xvY4u6rshvc3Y8D/YDhZvY9sBRolewjn5J4SO297G5DJDW6VFckjMyslLsfMrMKwLdAG3ffmbTsRWCpu7+Wxnc3k+JS3RBn06W6km2aeYiE1xQz+w6YAzyWrHAsBhqReHVUWnYDn5tZ81CHMrP3gTYknnMRyTLNPEREJMs08xARkSxT8RARkSxT8RARkSxT8RARkSxT8RARkSxT8RARkSz7P7WVvRzyjSY2AAAAAElFTkSuQmCC\n",
"text/plain": [
"