{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using Minuit in Python\n", "\n", "This tutorial focuses on how to use Minuit, a standard program for finding the minima of functions, from within Python. \n", "\n", "It assumes a general familiarity with Python and with the concept of ${\\chi}^2$ in statistics. The last example assumes some knowledge of the analysis package [ROOT](https://root.cern.ch/).\n", "\n", "You can copy the original of this notebook:\n", "\n", "On the Nevis particle-physics cluster:\n", "bash \n", "cp ~seligman/root-class/minuit-class.ipynb \n", "\n", "Elsewhere:\n", "\n", "wget https://www.nevis.columbia.edu/~seligman/root-class/files/minuit-class.ipynb or \n", "curl -O https://www.nevis.columbia.edu/~seligman/root-class/files/minuit-class.ipynb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Useful links:\n", "\n", "[iminuit tutorial](https://iminuit.readthedocs.io/en/latest/tutorials.html) \n", "... from which I have stolen liberally for this notebook.\n", "\n", "[iminuit complete docs](https://iminuit.readthedocs.io/en/stable/reference.html) \n", "This has everything, but it's a little tricky to read if you're not familiar with Python [structures](https://docs.python.org/3/tutorial/datastructures.html) and [function returns](https://www.geeksforgeeks.org/g-fact-41-multiple-return-values-in-python/). After going through this tutorial you'll be able to understand those docs and access all of Minuit's functionality.\n", "\n", "[Minuit's Wikipedia page](http://en.wikipedia.org/wiki/MINUIT) \n", "The links on the page will take you to the original Minuit documentation (which fully discusses the minimization techniques and formulae) and the original Minuit paper if you want to include it in your bibliography. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Installing iminuit\n", "\n", "In the Nevis particle-physics cluster, you don't have to do this. It's already installed when you do:\n", "\n", "module load python (or module load root/06.08.06 or earlier)\n", "\n", "module load python/3 (or module load root (06.10.xx or later))\n", "\n", "or when you're on the [notebook server](https://notebook.nevis.columbia.edu) and you select one of the Python kernels; here's more on the [notebook server](https://twiki.nevis.columbia.edu/twiki/bin/view/Main/IPython). \n", "\n", "*Exception:* The above doesn't apply if you're using a python distribution from your collaboration's software suite (e.g., LArSoft for MicroBooNE; Athena for ATLAS).\n", "\n", "Everyone else:\n", "\n", "pip install iminuit\n", "\n", "If you don't have admin access to your python installation:\n", "\n", "pip install --user iminuit\n", "\n", "(pip may be pip3, pip3.6, pip3.7, pip-3.6, pip-3.7, or something like that. It depends on the particulars of your python installation.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What is Minuit?\n", "\n", "Minuit finds the minima of functions. I'll leave the details of how it's done to the complete Minuit documentation linked above. Basically, Minuit repeatedly evaluates a function, steps in intervals, numerically computes gradients, and follows the slope of the function until it reaches a minimum. \n", "\n", "Super-fast summary (more below): \n", " * Migrad - find the function's minimum \n", " * Hesse - more accurate matrix of correlation coefficients (if minimum is parabolic/symmetric) \n", " * Minos - more accurate errors (if minimum is _not_ parabolic/symmetric) " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Basic setup for this notebook\n", "import math\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pprint\n", "from iminuit import Minuit " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Let's define the coefficients of a polynomial (easy to change and play with later)\n", "coeff = np.array([1, 2, 5]) # x**2 + 2x + 5" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Simple polynomial evaluation with numpy\n", "def polynomial(x):\n", " return np.polyval(coeff,x)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkfklEQVR4nO3dd3yV5f3/8dcnO4QMSAKEMAKELTsCigtcVHFVbRW1bqzrq61WW9v61dYO60RtVRwVt7iKioPhQIZogLBXmCEECIQEAhkkuX5/JPrlZ4EkkJP7jPfz8ciDkzvn5H4f0Heu3OO6zDmHiIgEnjCvA4iIyJFRgYuIBCgVuIhIgFKBi4gEKBW4iEiAimjOnaWkpLiMjIzm3KWISMCbP3/+Dudc6o+3N2uBZ2RkkJ2d3Zy7FBEJeGa28WDbdQhFRCRAqcBFRAKUClxEJECpwEVEApQKXEQkQKnARUQClApcRCRAqcBFRHyorLKa+z5YRl7Rvib/3ipwEREfent+Hi/N2UBBSXmTf28VuIiIj1RV1/Dc1+sY3CmJYzNaNfn3V4GLiPjIx0u3kldUxg0nd8PMmvz7q8BFRHzAOcczX66lW2ocp/du65N9qMBFRHzg6zU7WF6wmxtO6kZYWNOPvkEFLiLiE8/OXEvbhGjOG9TeZ/tQgYuINLElm0uYnbuTa0Z0IToi3Gf7UYGLiDSxZ75aS3x0BGOHdfLpflTgIiJNaMOOvXyytIDLhncmPibSp/tSgYuINKHnvl5HRFgY14zI8Pm+VOAiIk2kcE8Fb8/fzE8Hp9MmIcbn+1OBi4g0kZfmrGd/dQ3jTuraLPtTgYuINIHSiipembuRM/u0o2tqy2bZZ70FbmYxZvatmS0ys2Vmdn/d9pfMbL2Z5dR9DPR5WhERP/XGvE3sLq/ihpObZ/QNENGA51QAo5xzpWYWCcwys0/qvvYb59w7vosnIuL/yvdXM+HrdRzXNZlBnZp+0qpDqXcE7mqV1n0aWffhfJpKRCSAvD1/M4V7KrhlVGaz7rdBx8DNLNzMcoDtwDTn3Ly6L/3FzBab2WNmFn2I144zs2wzyy4sLGya1CIifmJ/dQ3PfLmWgR2TOL5bcrPuu0EF7pyrds4NBDoAQ83sGOB3QC/gWKA1cPchXjvBOZflnMtKTU1tmtQiIn5ics4W8ovLuHVUpk+mjD2cRl2F4pwrBr4ARjvnCuoOr1QA/waG+iCfiIjfqq5x/OvLXHqnJTCqV5tm339DrkJJNbOkusexwOnASjNLq9tmwPnAUt/FFBHxP58u3cq6wr3cPNI3CzbUpyFXoaQBE80snNrCn+Sc+8jMPjezVMCAHOCXvospIuJfnHM89UUuXVPj+MkxaZ5kqLfAnXOLgUEH2T7KJ4lERALAF6u2s6JgNw9d1J9wHy3YUB/diSki0kjOOZ76PJf0pFjOH5TuWQ4VuIhII81dt5MFm4r55cldiQz3rkZV4CIijfTPL3JJjY/m4qyOnuZQgYuINMLCTbuYnbuT60/sQkyk75ZLawgVuIhIIzz1eS5JLSK5bFhnr6OowEVEGmrx5mJmrNzOdSd0IS66IVdh+5YKXESkgZ6YsYbE2EiuPD7D6yiAClxEpEGWbC5h+ora0bevFytuKBW4iEgDjP9+9N0MixU3lApcRKQeS/NLmL5iG9ee0IUEPxl9gwpcRKRe42esISEmgqv8aPQNKnARkcNaml/CtOXbuO7Ern41+gYVuIjIYfnr6BtU4CIih/T96PvaE/xv9A0qcBGRQ3pixhri/XT0DSpwEZGDWralhKnLa688SYz1v9E3qMBFRA7q+9H31SO6eB3lkFTgIiI/smRzCZ8t8+/RN6jARUT+yyPTVpHUIpJrT/Df0TeowEVE/j/ZG4r4clUhvzy5m9/MeXIo9Ra4mcWY2bdmtsjMlpnZ/XXbu5jZPDPLNbO3zCzK93FFRHzHOcdDn60iNT6aK4/L8DpOvRoyAq8ARjnnBgADgdFmNhx4EHjMOZcJ7AKu9VlKEZFmMDt3J/PWF3HLyExio7xdbach6i1wV6u07tPIug8HjALeqds+ETjfFwFFRJqDc46Hpq4iPSmWS4Z6u9ZlQzXoGLiZhZtZDrAdmAasBYqdc1V1T9kMpB/itePMLNvMsgsLC5sgsohI05u+YjuL8or5n1MziY7w/9E3NLDAnXPVzrmBQAdgKNCroTtwzk1wzmU557JSU1OPLKWIiA/V1DgembqKjOQWXDi4g9dxGqxRV6E454qBL4DjgCQz+35RuA5AftNGExFpHlOWFLBy6x5+dXoPIsID5+K8hlyFkmpmSXWPY4HTgRXUFvlFdU+7Epjso4wiIj5TVV3DY9NW07NtPOf0b+91nEZpyLLKacBEMwuntvAnOec+MrPlwJtm9gCwEHjBhzlFRHzi/YX5rNuxl2evGEJYmHkdp1HqLXDn3GJg0EG2r6P2eLiISECqqKpm/Iw19O+QyBl92nodp9EC52CPiEgTe33eJjbvKuPOM3piFlijb1CBi0iI2lO+nyc/z2VEZjIndk/xOs4RUYGLSEh6buY6ivZWcvfoXgE5+gYVuIiEoO17ynnu6/WM6Z9G/w5JXsc5YipwEQk546evYX91DXee0dPrKEdFBS4iIWVdYSlvfpfH2GGdyEiJ8zrOUVGBi0hIeXjqKqIjwrh1VHevoxw1FbiIhIyFm3bx8ZKtXH9iV1Ljo72Oc9RU4CISEpxz/P2TlaS0jOL6k7p6HadJqMBFJCR8ubqQeeuL+J9Tu9MyuiGziPg/FbiIBL3qGseDn6ykc3ILLjm2k9dxmowKXESC3vsL81m5dQ93nNGTqIjgqb3geSciIgexr7KKhz5byYCOSYzpl+Z1nCalAheRoDZh5jq27a7gj2f3DrjpYuujAheRoLW1pJxnv1rH2f3SyMpo7XWcJqcCF5Gg9fDUVVTXOO4e3eBlfAOKClxEgtLS/BLeXbCZq0dk0Cm5hddxfEIFLiJBxznHA1OW06pFFDeNzPQ6js+owEUk6Exbvo1v1hXxq9O6kxgb6XUcn1GBi0hQqayq4W+frCSzTUsuHRo8N+0cTL0FbmYdzewLM1tuZsvM7La67feZWb6Z5dR9nOX7uCIih/favI2s37GXe87qRUR4cI9RGzIhQBVwh3NugZnFA/PNbFrd1x5zzj3su3giIg1XvK+Sx6ev4YTMFEb2bON1HJ+rt8CdcwVAQd3jPWa2Akj3dTARkcYaP2MNu8v38/uzewfsOpeN0ajfL8wsAxgEzKvbdIuZLTazF82sVVOHExFpqFVb9/Dy3I2MHdqJ3mkJXsdpFg0ucDNrCbwL3O6c2w08DXQDBlI7Qn/kEK8bZ2bZZpZdWFh49IlFRH7EOcf9Hy6jZXQEdwT4OpeN0aACN7NIasv7NefcewDOuW3OuWrnXA3wHDD0YK91zk1wzmU557JSU1ObKreIyA8+XbqVOWt3cscZPWgdF+V1nGbTkKtQDHgBWOGce/SA7QdO63UBsLTp44mIHF5ZZTUPTFlBr3bxjA3yywZ/rCFXoYwArgCWmFlO3bZ7gEvNbCDggA3ADT7IJyJyWM98tZb84jLeHDc86C8b/LGGXIUyCzjY6dyPmz6OiEjD5RXt45mv1jKmfxrDuyZ7HafZhdaPKxEJKn/9eAVhZtxzVm+vo3hCBS4iAWl27g4+WbqVm0d2o31SrNdxPKECF5GAs7+6hvs+WEan1i247sSuXsfxjApcRALOy3M3smZ7KX8c04eYyHCv43hGBS4iAWXb7nIem7aak3ukclrv4J/v5HBU4CISUP780XIqq2u4/9y+ITHfyeGowEUkYMxcXchHiwu4ZWQmGSlxXsfxnApcRAJC+f5q7p28lK4pcdxwcuieuDxQQ+7EFBHx3NNfrmXDzn28dt0woiNC98TlgTQCFxG/t66wlKe/XMt5A9szIjPF6zh+QwUuIn7NOce9k5cRHRnG788OzTsuD0UFLiJ+7YNFW5iVu4O7zuxJm/gYr+P4FRW4iPitkrL9PDBlBf07JDJ2WGev4/gdncQUEb/1yNRV7Cyt4MUrjyU8LLSv+T4YjcBFxC/N37iLV77ZyC+Oy6Bfh0Sv4/glFbiI+J2Kqmp+++5i0hJiuPPM0FnjsrF0CEVE/M6/vljLmu2l/PuqY2kZrZo6FI3ARcSvrN62h399mct5A9szsldoT1ZVHxW4iPiN6hrHXe8spmV0BPeO6eN1HL+nAhcRvzFxzgZy8or533P6ktwy2us4fk8FLiJ+Ia9oHw99topTeqZy3sD2XscJCPUWuJl1NLMvzGy5mS0zs9vqtrc2s2lmtqbuz1a+jysiwcg5xz3vLyHM4C8X9Av5eb4bqiEj8CrgDudcH2A4cLOZ9QF+C8xwznUHZtR97jOVVTW+/PYi4qH3FuTz9Zod3DW6F+khukDxkai3wJ1zBc65BXWP9wArgHTgPGBi3dMmAuf7KCOPT1/Nz56dS1W1Slwk2BTuqeDPU5YzpHMrrhiu2+Ubo1HHwM0sAxgEzAPaOucK6r60FWh7iNeMM7NsM8suLCw8opCZbVqSk1fMM1+tPaLXi4h/+v7Qyb7Kah68sB9hul2+URpc4GbWEngXuN05t/vArznnHOAO9jrn3ATnXJZzLis1NfWIQo7p354x/dMYP2MNy7fsrv8FIhIQ3l+Yz7Tl2/jNGT3JbBPvdZyA06ACN7NIasv7Nefce3Wbt5lZWt3X04DtvolY68/nHUNibBS/npSj4+EiQaCgpIz//WAZx2a04poTungdJyA15CoUA14AVjjnHj3gSx8AV9Y9vhKY3PTx/k+ruCj+9tN+rNy6hydmrPHlrkTEx5xz3P3uEqqqHQ9fPEAzDR6hhozARwBXAKPMLKfu4yzg78DpZrYGOK3uc586vU9bLhrSgae/WktOXrGvdyciPvLGt3nMXF3IPWf1onOyVpc/UvXOEuOcmwUc6sfjqU0bp373ntOH2bk7uGNSDlP+50RiIrW4qUggySvaxwNTljMiM5nLtEjDUQm4OzETYiJ58ML+rC3cyyNTV3kdR0QaoabGcefbiwgz4x8XDdBVJ0cp4Aoc4KQeqVw2rBPPz1rPt+uLvI4jIg300pwNzFtfxL3n9NENO00gIAsc4J6zetOhVSx3vr2IvRVVXscRkXqsLSzlH5+t5NRebbh4SAev4wSFgC3wuOgIHrl4IHm79nH/h8u8jiMih1FZVcPtb+YQExnOX3+quU6aSsAWOMDQLq256ZRuTMrezMdLCup/gYh44pFpq1iSX8KDF/anbUKM13GCRkAXOMDtp/VgQMckfvvuYrYUl3kdR0R+ZE7uDibMXMelQztxZt92XscJKgFf4JHhYYz/+UCqaxy/npRDdc1B7+gXEQ/s2lvJryctoktKHH8c09vrOEEn4AscICMljvvO7cs364p4dqYmvBLxB845fvfeEnbureCJSwbRIkqLEze1oChwgIuGdODsfmk8OnU1i3SXpojn3vouj0+XbeXOM3pyTHqi13GCUtAUuJnx1wv60SY+mtvfytGlhSIeWltYyv0f1t5tef2JXb2OE7SCpsABEltE8ujPB7Jh515dWijike8vGYyODOORiwfqbksfCqoCBxjeNfmHSws/XLTF6zgiIefhqf93yWC7RF0y6EtBV+BQe2nh4E5J/O69JazfsdfrOCIhY/rybUyYuY7Lh+uSweYQlAUeGR7GU2MHExFu3PTaAsr3V3sdSSTobd61jzveXkTf9gn84ew+XscJCUFZ4ADtk2J57GcDWVGwm/s/XO51HJGgVllVw82vL6SmxvGvywZrmudmErQFDjCyVxtuPKUbb3y7if8szPc6jkjQ+vsnK1mUV8w/LuqvBRqaUVAXOMAdp/dgaEZr7nl/CbnbS72OIxJ0Pl1awIuz13PV8Rn8pF+a13FCStAXeER4GE9cOojYyHBufm0BZZU6Hi7SVDbt3Mdv3lnMgA6J3HOWbpVvbkFf4ADtEmN47OcDWb19D/dOXup1HJGgUL6/mpten48BT40dTFRESNSJXwmZv/GTeqRy68hM3p6/mUnZeV7HEQl4D0xZztL83Tzys4F0bN3C6zghqd4CN7MXzWy7mS09YNt9Zpb/o1Xq/d5tp/VgRGYyf/jPUs2XInIUJn2Xx6vfbGLcSV05vU9br+OErIaMwF8CRh9k+2POuYF1Hx83bSzfCA8znrx0MKkto7nhlfkU7qnwOpJIwFm4aRd/+M9SRmQmc9eZPb2OE9LqLXDn3EwgaFYObh0XxbNXDKG4rJKbX1tAZVWN15FEAsb2PeX88tX5tEmI5qlLBxMRHjJHYf3S0fzt32Jmi+sOsbQ61JPMbJyZZZtZdmFh4VHsrukck57Igxf259sNRTwwRTf5iDREZVUNN726gJKy/Uy4IotWcVFeRwp5R1rgTwPdgIFAAfDIoZ7onJvgnMtyzmWlpqYe4e6a3nkD07n+xC68PHejTmqKNMCfPlpG9sZd/OOiAfRpn+B1HOEIC9w5t805V+2cqwGeA4Y2bazmcffoXrUnNd9fSo5Oaooc0lvfbeLVbzZxw0ldOXdAe6/jSJ0jKnAzO/B2qwuAgLy4OiI8jKcuHUybhGh++cp8tu8p9zqSiN9ZsGkXf/zPMk7snsJdo3t5HUcO0JDLCN8A5gI9zWyzmV0L/MPMlpjZYmAk8Csf5/SZVnFRTLgii+KySm58VTMXihyooKSMG1+dT9vEaJ68dBDhWpzBrzTkKpRLnXNpzrlI51wH59wLzrkrnHP9nHP9nXPnOucKmiOsr/Rpn8AjFw9k/sZd3P3uYpzTyvYieyuquPalbPZWVPPcL7JIaqGTlv5G1wDVObt/Gr85syeTc7YwfsYar+OIeKq6xnHbmwtZuXU3T44dRK92OmnpjyK8DuBPbjqlG+sK9/L49DVkJMdx/qB0ryOJeOKvH69g+ort/Om8vozs2cbrOHIIGoEfwMz420/7MaxLa+56ZzHZG4Lm/iWRBnvlm428MGs9V4/I4BfHZXgdRw5DBf4jURFhPHvFENJbxTLulfls3Kk1NSV0fLW6kPs+WMaoXm20LFoAUIEfRFKLKF686lhqnOPql76jZN9+ryOJ+NyqrXu4+bUF9GgbzxO64iQgqMAPoUtKHM9ePoS8on388tX5VFTp8kIJXtt3l3PNS9/RIiqcF67MomW0To8FAhX4YQzrmsw/LurP3HU7+fWkRVTX6PJCCT4lZfv5xYvfsmtfJS9ceSztk2K9jiQNpB+z9bhgUAcK91Tw149XkhwXxf3n9sVMv1pKcCjfX831L2eztrCUF686ln4dEr2OJI2gAm+AcSd1Y0dpJRNmriM5LprbTuvudSSRo1ZVXcOtbyzkuw1FjL9kECd295/J5qRhVOAN9NvRvdhRWsFj01eT3DKKy4d39jqSyBFzzvH795cybfk27junjyaoClAq8AYKCzMevLA/xfv288fJS2kdF8VZ/dLqf6GIH3p46ireys7jlpGZXDWii9dx5AjpJGYjRIaH8c+xgxncqRW3v5nDnNwdXkcSabQXZ63nn1+s5dKhHbnjjB5ex5GjoAJvpNi6y6wyUlow7pX5mkdcAsq78zfzp4+Wc2bftjxwfj+dkA9wKvAjkNQiionXDKVVXCS/eGEeS/NLvI4kUq/JOfn85p1FHN8tmfGX6EadYKACP0JpibG8ft1w4mMiufyFeSzfstvrSCKHNGVxAb+etIisjNY8f2UWMZHhXkeSJqACPwodW7fgjeuHExMRzuUvzGP1tj1eRxL5L58t28ptby5kUMck/n3VsbSI0rULwUIFfpQ6JbfgjXHDiQgzxj43j9ztpV5HEvnBjBXbuOX1BRyTnsi/rz6WON0iH1RU4E2gS0ocr18/HHCMfe4b1u/QDIbiva9WF3LjqwvonZbAxGuGEh8T6XUkaWIq8CaS2aYlr18/nKqa2hLXNLTipVlrdjDu5Wwy27Tk5WuGkhir8g5GKvAm1KNtPK9eO4yy/dVc/MxcHRMXT0xfvo1rJn5Hl5Q4Xr1umNayDGINWZX+RTPbbmZLD9jW2symmdmauj9b+TZm4OjTPoG3xh2HA37+7FyWbNYlhtJ8Jufkc8Or8+ndLp43xw2ndZzKO5g1ZAT+EjD6R9t+C8xwznUHZtR9LnV6tovn7RuOo0VUBGOf+4Zv12tpNvG9N77dxO1v5ZDVuRWvXT9cI+8QUG+BO+dmAj9uoPOAiXWPJwLnN22swJeREsc7Nx5HakI0v3hxHjNXF3odSYLY81+v43fvLeHkHqlMvGaoFmQIEUd6DLytc66g7vFWoG0T5QkqaYmxTLrhOLqmtOS6idl8unSr15EkyDjneHz6ah6YsoKz+6Ux4QrdpBNKjvokpnPOAYdcqsbMxplZtpllFxaG3ig0pWU0b4wbzjHpCdz8+gLezs7zOpIEiZoax58/WsHj09dw8ZAOPHHpIKIidF1CKDnSf+1tZpYGUPfn9kM90Tk3wTmX5ZzLSk0NzQnjE2MjeeXaYRzfLZnfvLOYx6evpvbnnsiRKd9fzS1vLODF2eu5ekQGD17YX3ObhKAjLfAPgCvrHl8JTG6aOMErLjqCF648louGdODx6Wu48+3FVFbVeB1LAtDO0grGPvcNnyzdyh/O7s29Y/oQpvIOSfWe6TCzN4BTgBQz2wz8L/B3YJKZXQtsBH7my5DBIioijIcu6k+n1i14dNpqCkrKePryIbrJQhps/Y69XPXvb9laUs7Tlw1m9DFaVCSUWXP+Kp+VleWys7ObbX/+7L0Fm7n73cV0SYnj31cPJV0rgUs9sjcUcf3L2ZgZz1+ZxeBOuv0iVJjZfOdc1o+364yHR346uAMTrx5KQUk55/9ztuYUl8OasriAsc/PI6lFFO/fdLzKWwAVuKeOz0zh3RuPJyo8jIuemcPknHyvI4mfqa5xPDJ1FTe/voABHRJ578bj6Zwc53Us8RMqcI/1aBvPf24eQf8OSdz2Zg4PfLScqmqd3BQo2befayd+x5Of5/LzrI68cu0wWunWeDmAbtfyA6nx0bx23TD+MmUFz89az7Itu3lq7CCSW0Z7HU08snLrbm54ZT5bisv4ywXHMHZoJ61fKf9FI3A/ERkexn3n9uXRnw1gwaZdnPPkLBZvLvY6lnjgo8VbuOCfcyirrObNccdx2bDOKm85KBW4n/np4A68e+PxmBkXPTOXSbpzM2RUVdfwt49XcMvrC+nTPoGPbj2BIZ11slIOTQXuh45JT+TDW0/g2IxW3PXOYn71Vg57yvd7HUt8KK9oHz97di7PzlzHFcM788b1w2mTEON1LPFzOgbup1rHRTHx6qH868u1jJ+xhuyNRYy/ZJAuHwtCk3Py+cP7tdPtj79kIOcNTPc4kQQKjcD9WER4GP9zancm3TAc5+DiZ+by1OdrqK7RPCrBoLSiil9PyuG2N3Po0S6ej287UeUtjaICDwBDOrfm49tO5Kx+aTw8dTWXPvcNW4rLvI4lRyEnr5izxn/Nfxbmc9up3Xlr3HA6tm7hdSwJMCrwAJEQE8kTlwzkkYsHsCy/hNGPz+Sd+Zs1q2GAqaiq5tFpq7no6TlU1zjeuuE4fnV6DyLC9b+iNJ6OgQcQM+PCIR3IymjFHZMWcefbi5ick89fL+in0VsAmL+xiLvfXULu9lIuGJTOfef21URmclQ0mVWAqqlxvDZvI3//ZCU1Du44owdXj+iiOaH9UGlFFQ99upKXv9lI+8RYHrjgGEb2bON1LAkgh5rMSgUe4LYUl/GH/yzl85XbGdAxiQcv7Eevdglex5I6X6zczu/fX0LB7nKuPC6DO8/sqfUqpdFU4EHMOccHi7Zw/4fL2V22n2tO6MItozJJiNGv517JK9rH3z9dyZTFBWS2acmDF/bXTTlyxA5V4BoKBAEz47yB6ZzYPZW/fbyC575ex7vzN/Or03twybEddYKsGZVWVPGvL3J5ftZ6wgx+dVoPfnlKV6IjtNCwND2NwIPQ0vwS/vTRcr5dX0SPti35w9l9OKlHaK5H2lyqaxxvZ+fx8NTV7Cit4IJB6fzmzJ6010Id0gR0CCXEOOf4dOlW/vrJCvKKyhjZM5XfndWbHm3jvY4WVJxzzFm7kwemrGBFwW4Gd0ri3nP6MrBjktfRJIiowENURVU1L83ewJOf57K3soqz+qVx66hMneg8Ss45ZuXu4IkZa/huwy7Sk2L57U96MaZ/mmYOlCanAg9xRXsreWHWOibO2UhpRRWj+7bj1lMz6ds+0etoAcU5x5erChk/Yw05ecWkJcZw4ynd+FlWR2IidZxbfEMFLgAU76vkxdkb+Pfs9ewpr+K03m24eWQmgzRJ1mFV1zhmrNjGk5/nsiS/hPSkWG4a2Y2LhnTQCUrxOZ8UuJltAPYA1UDVwXZwIBW4/ygp28/EORt4YdZ6Ssr2M6BDIpcP78w5A9prJHmAor2VTMrO47V5G8krKqNT6xbcMjKTCwanE6mre6SZ+LLAs5xzOxryfBW4/ymtqOK9BZt5ee5GcreXktQikp9ndeSyYZ3plBy6t+fn5BXz8twNfLS4gMqqGoZ1ac0Vx3XmzL7tVNzS7FTgcljOOeau28krczcydfk2apzj5B6pnDewPaf1bkt8CNwUtH13OVOWFPD+wnwWby4hLiqcCwanc8XwDHq209U74h1fFfh6YBfggGedcxMO8pxxwDiATp06Ddm4ceMR70+ax9aScl7/dhPvZOexpaScqIgwRvZM5ZwB7RnVqw0tooLn/q+dpRV8snQrHy7awrcbinAOerWLZ+ywTlwwKD0kfnCJ//NVgac75/LNrA0wDbjVOTfzUM/XCDyw1NQ4Fubt4sNFBUxZUkDhngpiI8MZ1bsNI3u2YURmMmmJgXWjinOOtYWlzFqzgxkrtzNn7U6qaxzdUuMY07895wxII7ONRtviX3x+FYqZ3QeUOucePtRzVOCBq7rG8e36Ij5avIXPlm1lR2klAN1S4zghM4XjM1MY3jXZL6dH3VpSzuzcHcxeu4PZuTvYtrsCgIzkFpzdP40x/dvTq128rt8Wv9XkBW5mcUCYc25P3eNpwJ+cc58e6jUq8OBQU+NYuXUPc9buYFbuDuatK6JsfzVhBt3bxNOnfQJ92yfQJy2BPu0TSGoR1Sy5nHNs213Bsi0lLN+ym2VbdrO8YDebivYBteuMHt8tmRMyUxiRmaI51CVg+KLAuwLv130aAbzunPvL4V6jAg9OlVU15OQVMzt3B0vya8tz6+7yH76enhRL19Q42ifGkpYU88OfaYmxpMZHExsZTmS4HXYE7JyjoqqGsspqtu0pZ0txGVuKyykoKaOguJz84jLWbC+laG/lD6/pkhJHn7QEBnZM4vjMZHq3SyBM86VLAGry2Qidc+uAAUeVSoJCVEQYQ7u0ZmiX1j9s21lawfKC2lHwsi272bRzLysK9rCjtOKg3yM8zIiNDCcmMoyYyHCiwsNqC3t/NeX7qynbX83BxhrhYUbb+GjSkmI5vXfbH0b/vdISNO+2BD39Fy4+kdwymhO7p3Ji9/9/FsSKqmq2lVSwpaSMgpIydpZW/lDQ5fvrCruymorqGmIiwomNCqsr9tqP2Mhw2iREk5YYS/ukGNrEx2gVIglZKnBpVtER4XRKbhHSNwmJNBXdUiYiEqBU4CIiAUoFLiISoFTgIiIBSgUuIhKgVOAiIgFKBS4iEqBU4CIiAapZ18Q0s0IgECcETwEatGhFkNH7Dj2h+t79/X13ds6l/nhjsxZ4oDKz7PrW+wxGet+hJ1Tfe6C+bx1CEREJUCpwEZEApQJvmP9a6zNE6H2HnlB97wH5vnUMXEQkQGkELiISoFTgIiIBSgXeSGZ2h5k5M0vxOktzMLOHzGylmS02s/fNLMnrTL5kZqPNbJWZ5ZrZb73O0xzMrKOZfWFmy81smZnd5nWm5mRm4Wa20Mw+8jpLY6nAG8HMOgJnAJu8ztKMpgHHOOf6A6uB33mcx2fMLBz4J/AToA9wqZn18TZVs6gC7nDO9QGGAzeHyPv+3m3ACq9DHAkVeOM8BtwFhMyZX+fcVOdcVd2n3wAdvMzjY0OBXOfcOudcJfAmcJ7HmXzOOVfgnFtQ93gPtWWW7m2q5mFmHYCzgee9znIkVOANZGbnAfnOuUVeZ/HQNcAnXofwoXQg74DPNxMiRfY9M8sABgHzPI7SXB6ndlBW43GOI6JFjQ9gZtOBdgf50u+Be6g9fBJ0Dve+nXOT657ze2p/1X6tObNJ8zGzlsC7wO3Oud1e5/E1MxsDbHfOzTezUzyOc0RU4Adwzp12sO1m1g/oAiwyM6g9jLDAzIY657Y2Y0SfONT7/p6ZXQWMAU51wX3jQD7Q8YDPO9RtC3pmFklteb/mnHvP6zzNZARwrpmdBcQACWb2qnPuco9zNZhu5DkCZrYByHLO+fPsZU3CzEYDjwInO+cKvc7jS2YWQe2J2lOpLe7vgLHOuWWeBvMxqx2VTASKnHO3exzHE3Uj8Dudc2M8jtIoOgYu9XkKiAemmVmOmT3jdSBfqTtZewvwGbUn8iYFe3nXGQFcAYyq+zfOqRuVip/TCFxEJEBpBC4iEqBU4CIiAUoFLiISoFTgIiIBSgUuIhKgVOAiIgFKBS4iEqD+H/viU4D88z4nAAAAAElFTkSuQmCC\n", "text/plain": [ "