\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.9?urlpath=lab/tree/nddata_demo.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",
"[nddata_demo.ipynb](../_static/notebooks/nddata_demo.ipynb) |\n",
"[nddata_demo.py](../_static/notebooks/nddata_demo.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to use the NDDataArray class"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"This notebook explains how to use the class [gammapy.utils.nddata.NDDataArray](..\/api/gammapy.utils.nddata.NDDataArray.rst)\n",
"\n",
"The NDDataArray is basically an numpy array with associated axes and convenience methods for interpolation. For now \n",
"only the scipy [RegularGridInterpolator](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html)\n",
"can be used, i.e. available interpolation methods are \"nearest neighbour\" and \"linear\". A spline interpolator\n",
"will be added in the future. The interpolation behaviour (\"log\", \"linear\") can be set for each axis individually.\n",
"\n",
"The NDDataArray is currently used in the following classes\n",
"\n",
"* [gammapy.irf.EffectiveAreaTable](..\/api/gammapy.irf.EffectiveAreaTable.rst)\n",
"* [gammapy.irf.EffectiveAreaTable2D](..\/api/gammapy.irf.EffectiveAreaTable2D.rst)\n",
"* [gammapy.irf.EnergyDispersion](..\/api/gammapy.irf.EnergyDispersion.rst)\n",
"* [gammapy.spectrum.CountsSpectrum](..\/api/gammapy.spectrum.CountsSpectrum.rst)\n",
"* Probably some more by now ...\n",
"\n",
"Feedback welcome!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As usual, we'll start with some setup ..."
]
},
{
"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": [
"from gammapy.utils.nddata import NDDataArray, DataAxis, BinnedDataAxis\n",
"from gammapy.utils.energy import Energy, EnergyBounds\n",
"import numpy as np\n",
"import astropy.units as u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1D example\n",
"\n",
"Let's start with a simple example. A one dimensional array storing an exposure in ``cm-2 s-1`` as a function of energy. The energy axis is log spaced and thus also the interpolation shall take place in log."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"NDDataArray summary info\n",
"energy : size = 10, min = 10.000 TeV, max = 100.000 TeV\n",
"Data : size = 10, min = 2.000 1 / (cm2 s), max = 20.000 1 / (cm2 s)\n",
"\n",
"DataAxis\n",
"Name: energy\n",
"Unit: TeV\n",
"Nodes: 10\n",
"Interpolation mode: log\n"
]
}
],
"source": [
"energies = Energy.equal_log_spacing(10, 100, 10, unit=u.TeV)\n",
"x_axis = DataAxis(energies, name=\"energy\", interpolation_mode=\"log\")\n",
"data = np.arange(20, 0, -2) / u.cm ** 2 / u.s\n",
"nddata = NDDataArray(axes=[x_axis], data=data)\n",
"print(nddata)\n",
"print(nddata.axis(\"energy\"))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 10. 12.91549665 16.68100537 21.5443469 27.82559402\n",
" 35.93813664 46.41588834 59.94842503 77.42636827 100. ] TeV\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8FdX9//HXJ4RFEDQCKoICURSR1QQIghjEBSt1p4pLtUotrWvtpm0Vqu33p9WqdamUIrhRioqoRWtRBHGLShQRQQVZJGIlBmQREEI+vz9mEkO4SYaQm8nyfj4e95GZM9vnXi75ZM45c465OyIiIpVJiTsAERGpG5QwREQkEiUMERGJRAlDREQiUcIQEZFIlDBERCQSJQwREYlECUNERCJRwhARkUhS4w6gOrVp08Y7deoUdxgiInVGbm7uV+7eNsq+9SphdOrUiXnz5sUdhohInWFmK6PuqyopERGJRAlDREQiUcIQEZFI6lUbhoh8Z/v27eTl5bF169a4Q5FaoFmzZnTo0IHGjRtX+RxKGCL1VF5eHi1btqRTp06YWdzhSIzcnYKCAvLy8ujcuXOVz5O0KikzO9jMZpvZYjP70MyuCcv3M7MXzWxJ+DOtnOMvDvdZYmYXJytOkfpq69attG7dWslCMDNat269x3ebyWzDKAR+4e5HAlnAFWbWDbgemOXuXYBZ4fpOzGw/YAzQH+gHjCkvsVSH3JXruH/2UnJXrkvWJURioWQhxarju5C0Kil3/wL4IlzeaGaLgfbA6UB2uNvDwBzgN2UOPxl40d3XApjZi8AwYEp1x5m7ch0XTMhhW2ERTVJTmDwqi4yOSctNIiJ1Vo30kjKzTkAf4C3ggDCZFCeV/RMc0h5YVWo9LyxLdO7LzWyemc3Lz8/f7dhylhWwrbCIIofthUXkLCvY7XOISGJ77713pfvcfffdbN68OemxPPTQQ1x55ZUV7jNnzhzeeOONkvVx48bxyCOPJDu0Ss2ZM4fhw4fHHUbyE4aZ7Q1MA6519w1RD0tQ5ol2dPfx7p7p7plt20Z6un0nWemtaZKaQiODxqkpZKW33u1ziEjVVSVh7NixIymxlE0Yo0eP5oc//GFSrlUXJTVhmFljgmQx2d2fCou/NLN24fZ2wJoEh+YBB5da7wCsTkaMGR3TmDwqi+tOOkLVUdLgJas9b86cOWRnZ3POOefQtWtXLrjgAtyde+65h9WrVzNkyBCGDBkCwMyZMxkwYABHH300I0aMYNOmTUAw9M/NN9/MoEGDeOKJJ8jOzubaa6/lmGOOoXv37rz99tsArF27ljPOOIOePXuSlZXFggULdonn3//+N/3796dPnz6ccMIJfPnll6xYsYJx48Zx11130bt3b1599VXGjh3LHXfcAcD8+fPJysqiZ8+enHnmmaxbF3xG2dnZ/OY3v6Ffv34cfvjhvPrqq5HfP8CsWbPo06cPPXr04NJLL+Xbb78F4IUXXqBr164MGjSIp556quRc33zzDZdeeil9+/alT58+PPPMMwB8+OGH9OvXj969e9OzZ0+WLFlSLf92O3H3pLwI7hIeAe4uU347cH24fD3w5wTH7gcsB9LC13Jgv8qumZGR4SISWLRo0W7tP2/FWj/i98975+tn+BG/f97nrVi7xzG0aNHC3d1nz57trVq18lWrVvmOHTs8KyvLX331VXd379ixo+fn57u7e35+vh977LG+adMmd3e/9dZb/Q9/+EPJfrfddlvJuY877jgfNWqUu7u/8sorftRRR7m7+5VXXuljx451d/dZs2Z5r1693N190qRJfsUVV7i7+9q1a72oqMjd3f/xj3/4dddd5+7uY8aM8dtvv73kGqXXe/To4XPmzHF39xtvvNGvueaakjiKj3/uued86NChu3wO5b3/LVu2eIcOHfzjjz92d/eLLrrI77rrrpLyTz75xIuKinzEiBF+6qmnurv7DTfc4I8++qi7u69bt867dOnimzZt8iuvvNIfe+wxd3f/9ttvffPmzbvEkeg7AczziL/Xk/kcxkDgIuADM5sflv0WuBV43MwuAz4DRgCYWSYw2t1HuftaM7sFeCc87mYPG8BFJDkStedV5x13v3796NChAwC9e/dmxYoVDBo0aOcYcnJYtGgRAwcOBGDbtm0MGDCgZPu555670/4jR44EYPDgwWzYsIGvv/6a1157jWnTpgFw/PHHU1BQwPr163c6Li8vj3PPPZcvvviCbdu2Vfpswvr16/n666857rjjALj44osZMWJEyfazzjoLgIyMDFasWBH5/bds2ZLOnTtz+OGHl5z3/vvvJzs7m86dO9OlSxcALrzwQsaPHw8Ed2DPPvtsyZ3P1q1b+eyzzxgwYAB/+tOfyMvL46yzzio5tjols5fUayRuiwAYmmD/ecCoUusTgYnJiU5Eyipuz9teWJSU9rymTZuWLDdq1IjCwsJd9nF3TjzxRKZMSdwhskWLFjutl+0qamYlVT0V7XfVVVdx3XXXcdpppzFnzhzGjh0b9W0kVPzeyntfpfcpvV+iWMuLuZi7M23aNI444oidyo888kj69+/Pc889x8knn8yECRM4/vjjd/etVEhjSYkIEF97XsuWLdm4cSMAWVlZvP766yxduhSAzZs388knn5R77NSpUwF47bXX2Geffdhnn30YPHgwkydPBoK2gzZt2tCqVaudjlu/fj3t2wcdLx9++OGEsZS2zz77kJaWVtI+8eijj5bcbeyJrl27smLFipL3W3zerl27snz5cj799FOAnRLoySefzL333luSbN577z0Ali1bRnp6OldffTWnnXZawrabPaWhQUSkREbHtBrv+HH55Zdzyimn0K5dO2bPns1DDz3EyJEjSxp///jHP5ZU2ZSVlpbGMcccw4YNG5g4MaiQGDt2LD/60Y/o2bMnzZs33ykhFBs7diwjRoygffv2ZGVlsXz5cgC+//3vc8455/DMM89w77337nTMww8/zOjRo9m8eTPp6elMmjRpj997s2bNmDRpEiNGjKCwsJC+ffsyevRomjZtyvjx4zn11FNp06YNgwYNYuHChQDceOONXHvttfTs2RN3p1OnTsyYMYOpU6fy2GOP0bhxYw488EBuuummPY6vLKvolqiuyczMdE2gJBJYvHgxRx55ZNxhJE12djZ33HEHmZmZcYdSZyT6TphZrrtH+hBVJSUiIpGoSkpE6qQ5c+bEHUKDozsMERGJRAlDREQiUcIQEZFIlDBERCQSJQwRSZq6Prx5VJ06deKrr76qamjVfp5kUcIQkVjV5uHNZWdKGCKSdHV1ePP8/HzOPvts+vbtS9++fXn99dcBKCgo4KSTTqJPnz785Cc/STgm1AMPPMCvf/3rkvWHHnqIq666CoAzzjiDjIwMjjrqqJJBBUtbsWIF3bt3L1m/4447Ssa7+vTTTxk2bBgZGRkce+yxfPTRRwA88cQTdO/enV69ejF48ODd/jeKJOqwtnXhpeHNRb6zy1DWE7+36+ut8cG2b79JvP3dYLhs3/TVrtsiqOvDm48cObIkzpUrV3rXrl3d3f2qq64qiWvGjBkOlLyHYmvWrPFDDz20ZH3YsGEl5yooKHB3982bN/tRRx3lX3311U6fxfLly0vej7v77bff7mPGjHF39+OPP94/+eQTd3fPycnxIUOGuLt79+7dPS8vz92DYc8Tqc3Dm4uIlKiLw5u/9NJLLFq0qGR9w4YNbNy4kblz55ZManTqqaeSlrbr+Ftt27YlPT2dnJwcunTpwscff1zyvu655x6mT58OwKpVq1iyZAmtW1c+OvCmTZt44403dhpavXjMrYEDB3LJJZfwgx/8oGS49eqmhLEbcleuI2dZAVnprTUzn9Q9P3qu/G1Nmle8vUXrirdHUBeHNy8qKuLNN99kr732qvSciZx77rk8/vjjdO3alTPPPBMzY86cObz00ku8+eabNG/enOzsbLZu3brTcampqRQVFZWsF28vKipi3333Zf78+ZQ1btw43nrrLZ577jl69+7N/PnzIyWh3aE2jIhyV67jggk5/GXmx1wwIafap7AUaahq8/DmJ510Evfdd1/JevEv6tLX+M9//lMyXWtZZ511Fk8//TRTpkwpuTtav349aWlpNG/enI8++oicnJxdjjvggANYs2YNBQUFfPvtt8yYMQOAVq1a0blzZ5544gkgSLDvv/8+ELRt9O/fn5tvvpk2bdqwatWqcj+3qlLCiCjRbGQisueKhzcfMmQIbdu2LRnevLjRurhRN5Hi4c1Hjx7Ngw8+CARDl8+bN4+ePXty/fXXVzi8+bHHHkubNm1Kyr///e8zffr0kkbve+65p+Rc3bp1Y9y4cQCMGTOGuXPncvTRRzNz5kwOOeSQcuPr1q0bK1eupF+/fgAMGzaMwsJCevbsyY033khWVtYuxzVu3JibbrqJ/v37M3z4cLp27VqybfLkyTz44IP06tWLo446qmRO71/96lf06NGD7t27M3jwYHr16lXZR7/bNLx5RMV3GMWzkdXkBDMiVaHhzaWsPR3ePGltGGY2ERgOrHH37mHZVKB4XsF9ga/dvXeCY1cAG4EdQGHUN5NMxbORqQ1DRBqqZDZ6PwTcBzxSXODuJV0czOwvwPpdDysxxN1r1SOPccxGJiKJaXjzmpe0hOHuc82sU6JtFnQv+AFQvTOUi8hO3D1Sbx6p/6qj+SGuRu9jgS/dfUk52x2YaWa5ZnZ5DcYlUm80a9aMgoKCavlFIXWbu1NQUECzZs326DxxPYcxEkjc0Tow0N1Xm9n+wItm9pG7z020Y5hQLgfK7akg0hB16NCBvLw88vPz4w5FaoFmzZqVPDhZVTWeMMwsFTgLyChvH3dfHf5cY2bTgX5AwoTh7uOB8RD0kqr2gEXqqMaNG5f7BLNIVcRRJXUC8JG75yXaaGYtzKxl8TJwErCwBuMTEZEEkpYwzGwK8CZwhJnlmdll4abzKFMdZWYHmdnz4eoBwGtm9j7wNvCcu7+QrDhFRCSaZPaSGllO+SUJylYD3wuXlwHV/4iiiIjsEQ0NIiIikShhiIhIJEoYIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpEoYSRR7sp13D97KbkrE08QLyJSl8Q1vHm9VzwH+LbCIppoDnARqQd0h5EkOcsK2FZYRJHD9sIicpYVxB2SiMgeKfcOw8z2i3B8kbt/XY3x1BtZ6a1pkprC9sIiGqemkJXeOu6QRET2SEVVUqvDV0UTAjcCNM1dAhkd05g8KoucZQVkpbdWdZSI1HkVJYzF7t6nooPN7L1qjqdeyeiYpkQhIvVGRW0YAyIcH2UfERGpB8q9w3D3rcXLZpYGHARsAVa4e1HZfUREpH6rqNF7H+AKYCTQBMgHmgEHmFkO8Dd3n10jUYqISOwqqpJ6ElgFHOvuR7j7IHfPdPeDgVuB00vN070LM5toZmvMbGGpsrFm9rmZzQ9f3yvn2GFm9rGZLTWz66v43kREpBpVVCV1YgXbcoHcSs79EHAf8EiZ8rvc/Y7yDjKzRsD9wIlAHvCOmT3r7osquZ6IiCRRpQ/umdlAM2sRLl9oZneaWcfKjnP3ucDaKsTUD1jq7svcfRvwL+D0KpxHRESqUZQnvR8ANptZL+DXwEp2vWvYHVea2YKwyipRn9P2BFVhxfLCMhERiVGUhFHo7k7wV/5f3f2vQMsqXu8B4FCgN/AF8JcE+yR6UNDLO6GZXW5m88xsXn5+fhXDEhGRykRJGBvN7AbgQuC5sI2hcVUu5u5fuvuOsFvuPwiqn8rKAw4utd6B4Inz8s45PmyMz2zbtm1VwhIRkQiiJIxzgW+By9z9fwTVQ7dX5WJm1q7U6pnAwgS7vQN0MbPOZtYEOA94tirXExGR6lPp8OZhkriz1PpnRGjDMLMpQDbQxszygDFAtpn1JqhiWgH8JNz3IGCCu3/P3QvN7ErgvwRjVU109w93832JiEg1s6B5on7IzMz0efPmxR1G1X3+LqxdBj3OiTsSEWkgzCzX3TOj7KsJlGqTubfDkhehVXvoqGG6RKR20QRKtckZf4N9D4GpF8K6lXFHIyKyk3IThpm1MrP/Z2aPmtn5Zbb9LfmhNUB7pcH5U2HHdpgyEr7dFHdEIiIlKrrDmETwTMQ04Dwzm2ZmTcNtWUmPrKFq0wVGTIL8xZDzQNzRiIiUqKgN41B3PztcftrMfge8bGan1UBcDdthQ+HiGXBwf3JXrtOsfSJSK1SUMJqaWUqpuS/+FHaPnQvsXSPRNWSdBpK7ch3XTvgPhxUt596Uo5k8KktJQ0RiU1GV1L+B40sXuPvDwC+AbckMSgI5ywr4FY8yLvUuuu34hJxlBXGHJCINWLkJw91/7e4vJSh/wd27JDcsAchKb83/cQlfsi/jGt/JsQcoT4tIfCp9DsPM9gV+CHQqvb+7X528sAQgo2Ma9486idc/2JcfvH8p+8+5DDo/H/SmEhGpYVGew3ieIFl8QDBpUpTJk6SaZHRMY+Twk2k08p9QsBRm/j7ukESkgYrypHczd78u6ZFIxdKPC57RaNc77khEpIGKcofxqJn92Mzamdl+xa+kRya7OvR4aL4fbN8K70yAejQOmIjUflHuMLYRDGf+O76byMiB9GQFJZX48Cl47hfB8CEn3RJ3NCLSQERJGNcBh7n7V8kORiLqNRI+z4U37oG994djroo7IhFpAKIkjA+BzckORHaDGZzyZ/gmP2gEb94Geo+MOyoRqeeiJIwdwHwzm00w8x6gbrWxS2kEZ/0DtqwLksaRw6FpVadaFxGpXJSE8XT4ktomtSmc90/Y8IWShYgkXZSE8SSw1d13AJhZI6BpxYdIjWnaEtqGyeKV2+GgPtDlhHhjEpF6KUq32lnAXqXW9wJ2GTJEYrZ9Cyx+Fv51PizVP4+IVL8oCaOZu5fM5BMuN6/sIDObaGZrzGxhqbLbzewjM1tgZtPDYUcSHbvCzD4ws/lmVocn6a5BjfeCHz4DbQ+HKefD0llxRyQi9UyUhPGNmR1dvGJmGcCWCMc9BAwrU/Yi0N3dewKfADdUcPwQd+8ddXJyIXio74fPQpvDgzuNT1+OOyIRqUeitGFcCzxhZqvD9XbAuZUd5O5zzaxTmbKZpVZzgHOihSmRNd8vuNN49AzYtCbuaESkHqk0Ybj7O2bWFTiCYMrWj9x9ezVc+1JganmXBWaamQN/d/fx1XC9hqNFa/jxbGgU/vNuXQ/N9ok3JhGp88qtkjKzQcXL7r7d3Re6+wfFycLMWplZ96pcNJzutRCYXM4uA939aOAU4AozG1zBuS43s3lmNi8/P78q4dRPxcli2StwV3dYPCPeeESkzquoDeNsM3vDzG4ys1PNrJ+ZDTazS83sUWAGO/eeisTMLgaGAxe4Jx49z91Xhz/XANOBfuWdz93Hu3umu2e2bdt2d8Op/w7sEbRpPH4RvFdefhYRqVy5VVLu/nMzSyNoZxhB0HaxBVhMUE302u5ezMyGAb8BjnP3hMONmFkLIMXdN4bLJwE37+61JFTcpjH1AnjmZ0H11ICfxR2ViNRBFbZhuPs64B/ha7eY2RQgG2hjZnnAGIJeUU2BF80MIMfdR5vZQcAEd/8ecAAwPdyeCvzT3V/Y3etLKU33hvMfh2mj4L83sKSoHTO39SArvTUZHTV7n4hEY+XUCtVJmZmZPm+eHtsoV9EOls+exClzDmJbodMkNYXJo7KUNEQaMDPLjfr4QpTnMKS+SGnE8ynZbCt0OvAlN/gk3l76v7ijEpE6otwqKTNr5+5f1GQwknxZ6a1pkprC4KKFXJz6X9Yv2wqDJmvwQhGpVEV3GBPNLMfMbjWzbDOL8pCf1HIZHdOYPCqLdif8jBUDb2Wf1a/BxFNg/edxhyYitVyFbRhm1oyg4foUYCDwGfAC8IK7f1YTAe4OtWFUwdKX4PFLgobxi2dAm8PijkhEalC1tWG4+1Z3f8HdrwlP+AuCaqz7zOztaohV4nbYCXDpC9A+A1q1izsaEanFdqvR292Xu/vf3P00YFClB0jdcGB3OG8yNGkB326EBU/EHZGI1EJV7iXl7tuqMxCpJd4aB0+Ngv/+Dop2xB2NiNQiasiWnQ26Djblw5v3wboVcOY49aASEUDPYUhZKY3ge3+GYbfBx8/DhBNh7fK4oxKRWqBKCcPM/lPdgUgtkzUaLnwKcGjUJO5oRKQWqOjBvaPL2wT0Tk44UqscOgR++iakpATtGYuegaPOhGCcLxFpYCpqw3gHeIUgQZSVcC5uqYdSwpvQD6fDtMtg8bNw+v1BjyoRaVAqShiLgZ+4+5KyG8xsVfJCklqp+9mwfhW89Af4aknQDTetU9xRiUgNqqgNY2wF26+q/lCkVjODQT+HC54MEsf4bPj05bijEpEaVG7CcPcn3f3jcrY9nbyQpFbrckIwX3haZ2ii7rYiDUlFc3oPr+zgKPtIPdT6UPjxy3Bw32D9nQdhw+p4YxKRpKuoDeN2M/ucxI3exf6PYG5vaWiKe0pt/B+8eBPM/hOcNT4Ym0pE6qWKEsaXwJ2VHL9Lg7g0MC0PDKqonrgYHjsHBv8Ssm8IHgAUkXql3ITh7tl7enIzmwgMB9a4e/ewbD9gKtAJWAH8IJw7vOyxFwO/D1f/6O4P72k8kiRtD4dRs+D5X8Hc2+GLBXD+VD2vIVLPJHtokIeAYWXKrgdmuXsXYFa4vpMwqYwB+gP9gDFmpomna7MmzeGM++GMB6Db6UoWIvVQUhOGu88F1pYpPh0ovlt4GDgjwaEnAy+6+9rw7uNFdk08Uhv1Ph/6XBAsz/8nPHMF7y3N4/7ZS8lducuNpIjUIXEMPnhA8Vzh4c/9E+zTHij9cGBeWCZ1yYbV+HuTSXt0KC+++BwXTMhR0hCpwypNGGbW3MxuNLN/hOtdaqA7baL6jIRzyZrZ5WY2z8zm5efnJzks2S2Df8n0XuNpTCFPNh7L5T6Ntz5dE3dUIlJFUe4wJgHfAgPC9Tzgj3twzS/NrB1A+DPRb5A84OBS6x2AhB393X28u2e6e2bbtm33ICxJho5Hn8gZRbfxfFEW16U+wdDmS+MOSUSqKErCONTd/wxsB3D3LVT8bEZlngUuDpcvBp5JsM9/gZPMLC1s7D4pLJM6JqNjGuNGDWXV8feyePh0jsg6Ndjw+btQVBRvcCKyW6LMuLfNzPYirBIys0MJ7jgqZWZTgGygjZnlEfR8uhV43MwuAz4DRoT7ZgKj3X2Uu681s1sIRswFuNndyzaeSx2R0TGNjI5pwGFBQcGn8OCJ0PEYOO0+SOsYa3wiEo25J2wa+G4HsxMJnofoBswEBgKXuPucpEe3mzIzM33evHlxhyGVcYd3H4b//h5wOOkWyPiRuuKKxMDMct09M8q+FVZJmZkBHwFnAZcAU4DM2pgspA4xg4xL4GdvQIdMmPFzmHyOqqhEarkKq6Tc3c3saXfPAJ6roZikodj3ELjoacidBJvXfjdZk7vuNkRqoSiN3jlm1jfpkUjDZAaZlwZjUAEsfQkmDoM1H8Ubl4jsIkrCGAK8aWafmtkCM/vAzBYkOzBpoLZvga8+hnGDYPb/wfatcUckIqEovaROSXoUIsWO/D4cMgD++1t45TZY+BScdk/Qo0pEYhXlDsPLeYkkR4s2wdwaFz4FO7bBV5/EHZGIEO0O4zmCBGFAM6Az8DFwVBLjEoHDhsLPciC1WbD+/lTwHdDzvO8ayEWkxlSaMNy9R+l1Mzsa+EnSIhIprUnz75YXPglLZkLuw3DqHXBgj/KPE5Fqt9t/prn7u4B6TUnNGzk1eDK8YAn8fTD85zewdX3cUYk0GJXeYZjZdaVWU4AMQMPCSs1LSYGjL4Kup8LLf4S3/g5HnALp2XFHJtIgRLnDaFnq1RSYQTAJkkg8mu8Hw++Eq9/9LlnMmwir58cZlUi9F6UN4w/Fy2aWAuzt7uocL/HbLz34uX0LvHI7bPwimO3v+Jug5QHxxiZSD0WZQOmfZtbKzFoAi4CPzexXyQ9NJKLGe8EVOTDgCnj/X3BvBrx2FxRGGlRZRCKKUiXVzd03EMy9/TxwCHBRUqMS2V3N9oGT/wQ/ews6DQraOL7+LO6oROqVKAmjsZk1JkgYz7j7dvTgntRWbQ6D8/8FV7wNbboEZXNug9XvxRuXSD0QJWH8HVgBtADmmllHYEMygxLZY60PDX5+8xW8PR7GZ8O0UXywcAH3z15K7sp1sYYnUhdVOoFSwoPMUt29MAnx7BFNoCQJbd0Ar/+Vojfuo7CwkId3nMzf7Wz+Pur4cCZAkYar2iZQCk+2j5ndaWbzwtdfCO42ROqGZq1g6I080vcpni4axNmNXmFHYSE5ywrijkykTolSJTUR2Aj8IHxtACYlMyiRZOhxZDduYjRDtt3NltRWZHVOgykj4b3HYEetu2EWqXWizOk93917V1YW+YJmRwBTSxWlAze5+92l9skGngGWh0VPufvNlZ1bVVJSmdyV68hZVkBWemsy9vs2SBir34XWXeD438GRp2tgQ2lQdqdKKspotVvMbJC7vxaefCCwparBufvHQO/wXI2Az4HpCXZ91d2HV/U6IolkdEzbud3ixy/DR8/By7fAE5fAgT3hvH/CvgfHFqNIbRUlYfwUeNjM9iEY4nwtcHE1XX8o8Km7r6ym84nsHjM4cngwJtUHT8D7U6DlgcG2b74K5uYQESDa0CDzgV5m1ipcr84utecBU8rZNsDM3gdWA7909w+r8boiO0tpBL3OC14A276BB46B/btB9g1wSP944xOpBaL0kmptZvcAc4DZZvZXM2u9pxc2sybAacATCTa/C3R0917AvcDTFZzn8uIeXPn5GkRXqok1gmOuhi8XwsST4NEzYdXbcUclEqsorXv/IhjO/GzgnHB5aoVHRHMK8K67f1l2g7tvcPdN4fLzBE+bJ6wbcPfx7p7p7plt27athrBEgMbN4Jgr4Zr34cRb4IsF8OCJ8Pm7cUcmEpsoCWM/d7/F3ZeHrz8C+1bDtUdSTnWUmR1oZhYu9wvjVKd5qXlNWsDAq+HaBXDGODioT1D+7qPw2VvxxiZSw6I0es82s/OAx8P1cwjm+a4mamK0AAANNElEQVQyM2sOnEipqV7NbDSAu48Lr/FTMysk6JF1nlflkXSR6tKkBfQeGSwXboNXboP1q6DTsTD4l9D5uKABXaQei/IcxkaCJ7t3hEWNgG/CZXf3VskLb/foOQypMdu+gdyH4I17g3k42mfC926H9kfHHZnIbqnWoUHcvaW7p7h74/CVEpa1rE3JQqRGNWkRzL9xzfsw/C7Y/BU0bh5s27IOinZUfLxIHRSll9RlZdYbmdmY5IUkUoekNoXMS+Gq92D/rkHZv6+F+/rCvEmwXZNTSv0RpdF7qJk9b2btzKwHkEMwv7eIFCs9nEj3s4MBD2dcC3f3gFf/Alu+ji82kWoS5cG9883sXOADYDMw0t1fT3pkInVVt9PgyO/D8rnw+t0w6+agzWPoTXFHJrJHKk0YZtYFuAaYBhwJXGRm77n75mQHJ1JnmUH6ccHriwXQsl1QvvQl+ODJoP3jwB7xxiiym6JUSf0buNHdfwIcBywB3klqVCL1SbuesHf4UOna5bDoGRg3CB45HZa8COoxLnVElG61rcqOH2VmXdx9SVIjqwJ1q5U6Ycu6oEvuW38PuuQedgJcOC3uqKSBqpZutWb2awiG6TCzEWU2/2gP4hNp2PZKg0E/h2sWwJnjoVfxA4HfwtzbYeP/4o1PpBwVVUmdV2r5hjLbhiUhFpGGJbUJ9DoXepwTrK98A17+E9x1FEwbBXm6W5bapaKEYeUsJ1oXkT116BC4Khf6/hg+fgEmDIV/DIXNa8lduY77Zy8ld+W6uKOUBqyiXlJeznKidRGpDq0PhVNuDaaLnT8FVswldw1c8GAO2wqLaJKawuRRWTvPGihSQyq6w+hlZhvCsaR6hsvF6+oPKJJMTVtC/8vh3MfIWb6WbYVFFDlsLywiZ5kGbpZ4lHuH4e6NajIQEUksK701TVJT2F5YROPUFLLS93j+MpEqiTK8uYjEKKNjGpNHZZGzrICs9NaqjpLYKGGI1AEZHdOUKCR2UZ70FhERUcIQEZFolDBERCSS2BKGma0wsw/MbL6Z7fJIqwXuMbOlZrbAzDT3pYhIjOJu9B7i7l+Vs+0UoEv46g88EP4UEZEY1OYqqdOBRzyQA+xrZu3iDkpEpKGKM2E4MNPMcs3s8gTb2wOrSq3nhWUiIhKDOKukBrr7ajPbH3jRzD5y97mltica4HCXMazCZHM5wCGHHJKcSEVEJL47DHdfHf5cA0wH+pXZJQ84uNR6B2B1gvOMd/dMd89s27ZtssIVEWnwYkkYZtbCzFoWLwMnAQvL7PYs8MOwt1QWsN7dv6jhUEVEJBRXldQBwHQzK47hn+7+gpmNBnD3ccDzwPeApcBmNMufiEisYkkY7r4M6JWgfFypZQeuqMm4RESkfLW5W62IiNQiShgiIhKJEoaIiESihCEiIpEoYYhIrZK7ch33z15K7sp1cYciZcQ9+KCISIncleu4YEIO2wqLaJKawuRRWZppsBbRHYaI1Bo5ywrYVlhEkcP2wiJylhXEHZKUooQhIrVGVnprmqSm0MigcWoKWemt4w5JSlGVlIjUGhkd05g8KoucZQVkpbdWdVQto4QhIrVKRsc0JYpaSlVSIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpEoYYiISCRKGCIiEkmNJwwzO9jMZpvZYjP70MyuSbBPtpmtN7P54eummo5TRER2FsfQIIXAL9z9XTNrCeSa2YvuvqjMfq+6+/AY4hMRkQRq/A7D3b9w93fD5Y3AYqB9TcchIiK7J9Y2DDPrBPQB3kqweYCZvW9m/zGzoyo4x+VmNs/M5uXn5ycpUhERiS1hmNnewDTgWnffUGbzu0BHd+8F3As8Xd553H28u2e6e2bbtm2TF7CISAMXS8Iws8YEyWKyuz9Vdru7b3D3TeHy80BjM2tTw2GKiEgpcfSSMuBBYLG731nOPgeG+2Fm/Qji1FyNIiIxiqOX1EDgIuADM5sflv0WOATA3ccB5wA/NbNCYAtwnrt7DLGKiNRquSvX1dgMhTWeMNz9NcAq2ec+4L6aiUhEpG7KXbmOCybksK2wiCapKUwelZXUpKEnvUVE6qicZQVsKyyiyGF7YRE5y5Jbc6+EISJSR2Wlt6ZJagqNDBqnppCV3jqp14ujDUNERKpBRsc0Jo/Kqr9tGCIiUn0yOqYlPVEUU5WUiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpFYfRqiyczygZVxx7GH2gBfxR1ELaHPYmf6PHamz+M7e/JZdHT3SHND1KuEUR+Y2Tx3z4w7jtpAn8XO9HnsTJ/Hd2rqs1CVlIiIRKKEISIikShh1D7j4w6gFtFnsTN9HjvT5/GdGvks1IYhIiKR6A5DREQiUcKIiZkdbGazzWyxmX1oZteE5fuZ2YtmtiT8WTPDUNYCZtbIzN4zsxnhemczeyv8LKaaWZO4Y6wpZravmT1pZh+F35EBDfy78fPw/8lCM5tiZs0a0vfDzCaa2RozW1iqLOH3wQL3mNlSM1tgZkdXVxxKGPEpBH7h7kcCWcAVZtYNuB6Y5e5dgFnhekNxDbC41PptwF3hZ7EOuCyWqOLxV+AFd+8K9CL4XBrkd8PM2gNXA5nu3h1oBJxHw/p+PAQMK1NW3vfhFKBL+LoceKC6glDCiIm7f+Hu74bLGwl+IbQHTgceDnd7GDgjnghrlpl1AE4FJoTrBhwPPBnu0pA+i1bAYOBBAHff5u5f00C/G6FUYC8zSwWaA1/QgL4f7j4XWFumuLzvw+nAIx7IAfY1s3bVEYcSRi1gZp2APsBbwAHu/gUESQXYP77IatTdwK+BonC9NfC1uxeG63kECbUhSAfygUlhFd0EM2tBA/1uuPvnwB3AZwSJYj2QS8P9fhQr7/vQHlhVar9q+2yUMGJmZnsD04Br3X1D3PHEwcyGA2vcPbd0cYJdG0qXvlTgaOABd+8DfEMDqX5KJKybPx3oDBwEtCCodimroXw/KpO0/ztKGDEys8YEyWKyuz8VFn9ZfPsY/lwTV3w1aCBwmpmtAP5FUNVwN8GtdPE0wh2A1fGEV+PygDx3fytcf5IggTTE7wbACcByd8939+3AU8AxNNzvR7Hyvg95wMGl9qu2z0YJIyZhHf2DwGJ3v7PUpmeBi8Pli4Fnajq2mubuN7h7B3fvRNCY+bK7XwDMBs4Jd2sQnwWAu/8PWGVmR4RFQ4FFNMDvRugzIMvMmof/b4o/jwb5/SilvO/Ds8APw95SWcD64qqrPaUH92JiZoOAV4EP+K7e/rcE7RiPA4cQ/EcZ4e5lG7vqLTPLBn7p7sPNLJ3gjmM/4D3gQnf/Ns74aoqZ9SboANAEWAb8iOAPvAb53TCzPwDnEvQufA8YRVAv3yC+H2Y2BcgmGJX2S2AM8DQJvg9hUr2PoFfVZuBH7j6vWuJQwhARkShUJSUiIpEoYYiISCRKGCIiEokShoiIRKKEISIikShhiMTIzHaY2XwzOygceXW+mX1mZvnh8vxw6JhEx/7RzG4pU5ZpZgvC5VfNbFPYRVdkj6lbrUgVmVlqqbGMqnqOTe6+d5mySwhGZr2ykmOPAqa7++Glyu4ACtz9/4XrrwFXuvv8PYlTBHSHIQ2EmV1oZm+Hf7H/3cwaheWbzOxPZva+meWY2QFheVszm2Zm74SvgWH5WDMbb2YzgUfCp48fD+cdmBreJWSa2WVmdlep6//YzO5MGFy0+E8xszfN7N3wOi3c/UNgq5llhPsYMILgYTaRaqeEIfWemR1J8JTwQHfvDewALgg3twBy3L0XMBf4cVj+V4K5FvoCZxMOux7KAE539/OBnwHr3L0ncEu4DYJf2qeF44VB8KT2pCrGvz/B4IND3f1oYAHB3CEAUwiGU4FgTK7V7r68KtcRqUxq5buI1HlDCX6RvxP8Ec5efDdQ2zZgRricC5wYLp8AdAv3B2hlZi3D5WfdfUu4PIggueDuC4vbD9z9GzN7GRhuZouBxu7+QRXjPwboBrwRxtMEeC3cNgV4xcx+TZA4plTxGiKVUsKQhsCAh939hgTbtvt3DXk7+O7/RAowoFRiCE4U/ML+psy5yzOBYHywj6ji3UWpa7zg7heV3eDuK8xsNXAscCbf3eGIVDtVSUlDMAs4J6zaKZ4LuWMlx8wEShqdK+hp9Brwg3CfbkCP4g3h8OQHA+ezZ3/5vwEcFw7GiJm1MLMupbZPAe4hGPn4f3twHZEKKWFIvefui4DfAzPDKqMXgcqmrLwayAwbsxcBo8vZ729A2/C8vyFoX1hfavvjwOvuvm4P4v+SYL7qqWb2PkECObzULo8D3VFjtySZutWK7IGwt1Vjd99qZocS3M0c7u7bwu0zCBrPZ5Vz/C7daqs5PnWrlWqjOwyRPdMceC38y3868FN332Zm+5rZJ8CW8pJFaEPxg3vVHZiZvUowV8L26j63NEy6wxARkUh0hyEiIpEoYYiISCRKGCIiEokShoiIRKKEISIikShhiIhIJP8fmx3gKXXcYdUAAAAASUVORK5CYII=\n",
"text/plain": [
"