{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

## Minimization

" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import chi2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Minimization, not using least squares

\n", "\n", "Let's try a general minimization function, for Poisson uncertainties on a Gaussian profile. First write out obligatory routine to return a Gaussian. Here, we'll add another constant background/continuum value. We will also add an option to return the vector of derivatives of the Gaussian with respect to each parameter.\n", "
\n", "if \n", "$$f = A \\exp{-0.5(x-x_0)/\\sigma^2} + c$$\n", "then the derivatives are:\n", "$${df\\over dA} = \\exp{-0.5(x-x_0)^2/\\sigma^2}$$\n", "$${df\\over dx_0} = f {(x-x_0)\\over \\sigma^2}$$\n", "$${df\\over d\\sigma} = f {(x-x_0)^2\\over \\sigma^3}$$\n", "$${df\\over dc } = 1$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1fn48c/JSkJISCBhT4JCtVY0SkDcU7GKfN3w64IdEQSNa39Ctf0iWGtrU7V1Qa2KQRC0I4o7WBURG2y1ARGBIGJBSULY1yQQyDI5vz9mgkmYLbPcuXPzvF+vvMjcuTP34c69T+6c+5xzlNYaIYQQ0S8m0gEIIYQIDUnoQghhEZLQhRDCIiShCyGERUhCF0IIi4gzcmM9e/bUubm5Rm5SCCGi3ldffbVHa53paz1DE3pubi4rV640cpNCCBH1lFIV/qwnTS5CCGERktCFEMIiJKELIYRFSEIXQgiLkIQuhBAWIQldCCEsQhK6EEJYhCR0IYSwCEnoIjoVFDh/hBBHSUIXQgiLkIQuhBAWIQldCCEsQhK6EEJYhCR0IYSwCEnoQghhEZLQhRDCIiShCyGERUhCF0IIi5CELkR70gtVRCmfCV0pNUAp9U+l1LdKqW+UUne7lmcopZYopTa6/k0Pf7hCCCE88ecKvQm4R2v9U2AEcKdS6iRgKrBUaz0YWOp6LIQQIkJ8JnSt9Xat9SrX77XAt0A/4Apgnmu1ecCV4QpSCCGEbx1qQ1dK5QKnAcuBXlrr7eBM+kCWh9cUKqVWKqVW7t69O7hohRBCeOR3QldKpQBvAZO11jX+vk5rXay1ztda52dmZgYSoxBCCD/4ldCVUvE4k7lda/22a/FOpVQf1/N9gF3hCVEIIYQ//KlyUcBs4Fut9ROtnloIjHf9Ph54L/ThCSGE8FecH+ucDYwDypRSq13LpgGPAAuUUpOASuCa8IQohBDCHz4Tutb634Dy8PTI0IYjhBAiUNJTVFhSwdwCCuYWRDoMIQwlCV0IISxCEroQQliEJHQRfex2KC2FZcsgN9f5WAghCV1EGbsdCguhvt75uKLC+ViSuhCS0EWUmT4d6uraLqurcy4XopOThC6iS2Vlx5YL0YlIQhfRJTu7Y8uF6EQkoYvoUlQEycltlyUnO5cL0clJQhfRxWaD4mJITHQ+zslxPrbZIhuXECbgz1guQpiLzQazZjl/LymJaChCmIlcoQshhEVIQhdCCIuQhC6EEBYhCV0IISxCEroQrck4MSKKSUIXooWMEyOinCR0IVrIODEiyklCF6KFjBMjopwkdCFayDgxIspJQheihYwTI6KcdP0X1rR6dcdf0zIezKRJzhujOTnOZC7jxIgoIQldiNZknBgRxaTJRQghLEISujCnggLnTwDsZXZKU2tY1r2a3Bm52Mukjlx0DtLkIizFXmancFEh9bEagIrqCgoXFQJgGxKatvCCuQUAlEwoCcn7CREqcoUuLGX60unUNbbtHFTXWMf0pdI5SFifJHRhKZXV7jsBeVoeFkE0FwkRDEnowlKy09x3AvK0XAgrkYQuLKVoZBHJ8W07ByXHJ1M0UjoHCeuThC4sxTbERvFlxSQ6FGjIScuh+LLikN0QFcLMpMpFWI5tiI1Zs+8EoOTB8sgGI4SBJKGL6CS9OIU4hs8mF6XUHKXULqXUulbLHlRKbVVKrXb9jA5vmEKETsHcgqO15EJYiT9t6HOBUW6WP6m1znP9fBDasIQQQnSUz4Sutf4M2GdALEIIIYIQTJXLXUqpta4mmXRPKymlCpVSK5VSK3fv3h3E5oQQQngTaEJ/HjgeyAO2A497WlFrXay1ztda52dmZga4OSEESPu/8C6ghK613qm1dmitm4FZwPDQhiWEEKKjAkroSqk+rR6OAdZ5WlcIIYQxfNahK6XmAwVAT6VUFfB7oEAplQdooBy4NYwxCiGE8IPPhK61vt7N4tlhiEWI6BDIfKVCGEDGchFCCIuQhC6EwaRSRYSLJHQhhLAISejCeux2XptTw6dPVUNuLtg7OEl0SYkM/iWikiR0YT52O5SWwrJlHU/IdjsUFtK7VjsP7ooKKCzseFIXIgpJQhfm4krI1Nc7H3c0IU+fDnVtJ4mmrs65XAiLk4QuzCXYhFzpYTJoT8uFsBBJ6MJcgk3I2R4mg/a0vIPsZXZKU2tY1r2a3Bm52MuiqylHKmysTRK6MJdgE3JRESS3nSSa5GTn8iDZy+wULiqkPlaDgorqCgoXFbZN6sG0/wsRJEnowlyCTcg2GxQXs6ObohkgJweKi53LgzR96XTqGts2B9U11jF9qas5KNj2fyGCJAldmIsrIZOY6HwcSEK22Rg7MZUL7k6D8vKQJHOAymr3zT5Hl4fghqw0iYhgyCTRwnxsNpg1y/m7ierBs9OyqaiucLsckBuyIuLkCl1ERLivREtW51GyOi+k71k0sojk+LbNQcnxyRSNdDUHhfmGrBC+SEIXwk+2ITaKLysm0aFAQ05aDsWXFWMb4mrSCeMNWSH8IU0uQnSAbYiNWbPvBKDkwfJ2T7oS+6RJzhujOTnOZB6iNnwhfJGELkQombT9X3QO0uQiOhV7mZ3SqlKWVSyLuo5B0Ry7MIYkdNFpHO0Y5HDWibvtGGRS0Ry7MI4kdBEeBQXOHxPx2THIxKI5dmEcSeii0/DZMcjEojl2YRxJ6KLTONoByM/lZhLNsQvjSEIXnYbPjkEmFs2xC+NIQhedxtGOQbHOcWKO6RhkYtEcuzCO1KGLTsU2xMasr5x14iUTSjr+Bq75SrNqNbyba2jHoaBjF5YnV+hC+CsE85VGspZc6titTxK6EP4KcnjcSNaSSx175yBNLqLTa3I0s/XAYX7Yc4jyPYfYvOcQNYcbyc5IZmBmVwb2TGFgz66kBTk8rrda8nC3hUdy28I4ktBFp1N/JItDNSdy87yVbN5zkMp9dTQ69NHnuyXGkZoUz8I122j+cTFfpGbSt3rXsW/o5/C4kawllzr2zkESuugUmhzNLP5mJ3O/2Mz28vGgmtgSX8fgrG5c9LPeDOzZ9ehPj64JKKWob3KwZd9hNu85xOY9B/nsyGSufP4PdGmoP/q+DYldqJ32e3r4EYPPCTLCKJLbFsaRhC4sbd+hBuavqOTvpRVsrz5CdkYy6Y2vkOL4jMVTtnl9bWJcLIOyUhiUlQL0gvPug2HZ7Lh9HFm1mv09e/PQWTfw3uYsLnx5JTedlcuZx/dAeXi/opFFFC4qbNP0YVQteSS3LYwjCV1Y0rfJmbzU+3TefXgpDU3NnDOoJw9dcTI/PzGLkb8eE/gb22yM/dI1HvqM7fxf9WH6l1by6opKlqzfyU96pTAh6xSu2r2eLu1f6mqrnvTGOOpjNDndcygaWWRIG/bRbb83iXpHPTlpxm1bGEcSurCUww0O/rJ4Ay+dMoEkRwPXDO3PhLNyGdyrW1i21yctiXsvPoG7LhjEojXbmPtFOdOOu5gX+g7n8fJ95OdmtFnf6wQZYSZ17Nbns2xRKTVHKbVLKbWu1bIMpdQSpdRG17/p4Q1TWIlf9dAlJR2eIGJV5X7+5+l/8dLn5YzfsYrSVTMpGjMkbMm8tS7xsVyTP4D3f3UOf1+/gGYU17zwH4r+sZ4jjQ6/3kPqxEWw/KlDnwuMardsKrBUaz0YWOp6LIRP4aiHrm9y8OhHG7j6+S+ob2rm1ZvP4A9zf0fa0sUBvV8wE1grpbj/uPdQ/JpfDs9m1r82c+kz/2bNlgNeXyd14iIUfCZ0rfVnwL52i68A5rl+nwdcGeK4hEWFelzvdVurufyZz3m+5HuuzR/AR5PP5axBPUMRalBiqKdozBBenjicQ/VNXPX8Fzz+8Xc0NDW7Xd+w8c5Xr3b+CEsKtA29l9Z6O4DWertSKsvTikqpQqAQINvPel1hXaGqh250NPPsPzfxt083kdE1gZcmDOPnJ3o8DCPmvJ9k8tHk8/jjovU88+kmPvl2Fw1qAAl6S5v1pE5chELYu/5rrYu11vla6/zMzMxwb06YXCjG9d5Ve4RrZv6HGZ9s5NJT+vDxlPNMmcxbpCXF8/i1pzLrxnx219azLeFP1MZe0GYdGe9chEKgCX2nUqoPgOtfN93nRKdlt0NpKSxbBrm5bQavCnZc7027ahnz7Bd8t6OW52ynM2PsaXRPTghl9GHzi5N6sWTKeSQ1l7E3fhKPfrSBZldXVBnvXIRCoE0uC4HxwCOuf98LWUQiurlGJKTe1ZuyZURCAJstqHro0h/2UvjyShLjY1lw65kM6Z8Wrv+FVyWr87yvkOf5+fSuCWQ1PsE+fSPPl0DV/sM8ds0pIasTb7mZK2WJnZPPhK6Umg8UAD2VUlXA73Em8gVKqUlAJXBNOIMUUcTbiISuccMDqYd+b/VWfvPGWrJ7JDP3pmH0T0/2/SKTUjST0TSXOy/7FY98uIGdNUcoHjdU6sRF0HwmdK319R6eGhniWIQVBDkiYXtaa54r+Z6/Lv6OEcdl8MIN+aQlxwcRoDko4Lbzj6df9yTuWbCGq57/gnk3DY90WCLKyXjoIrQ8VTIFUOHU5Ghm2jtl/HXxd1yR15d5E4dbIpm3dtmpffn7zWew92ADY577nPrDvSMdkohiktBFaBUVQXK75pDkZOfyDjhY38SkeSuZv2ILd/18EDOuyyMxLjaEgZrH8IEZvHX7WSQlxLKj8jrqao/3uG7JhBLvzTFSZ96pSUIXoWWzQXExJDonMyYnx/m4A/NuVh9uxDarlH9v2sPDVw3h3otPQClPYxhaw6CsFN6+/WziE/eya+uVvL2qyv2KBQXOHyHckMG5ROjZbDDLeXOvo+Ox1Bxp5MbZy1m/vYaZNwzlFyf1Cn18JpXZLZHe2a+zq2oM97wRg1Iw5rT+kQ5LRBG5QhemUXOkkXGzV7B+ew3P2zpXMm8RE9NIVv+3GTGwB/csWMN7q7dGOiQRReQKXZiC88p8Beu3VfPsL0/nwjAmc5915BEWE9PE7An5TJz7JVNed7aHX5HXL8JRiWggV+giLAryVlOQ59/NudojjYyfs4J1W53J/KKfdeJKD9dNzeSEOOZMGMbwgRlMeX21XKkLv0hCFwErmNydgsndg3qP2iON3DhnBWVV1Txr6wTJ3G7ntTk1fPpU9THDIrTXktSH5TqT+sI13qfM80fJ6jzTf0MRgZOELiKm2ZHAeFcy/9svT+fiTpDMKSykd612nngtwyL4SOov3TSM/NwMJr/2NYt6nGBYuCL6SEIXEdHsiGdn1dWsrarmb788jVEnWzyZg/dhEbxITojjpQnDyM/JYPKgS3k/Q5K6cE8SujDckUYHO6uuov5wH565/jRGndwn0iEZI4hhEbomOq/Uh9Zu5e7Bl/LJ+p0hDk5YgSR0YahGRzN32FdRf3gAPfv+g0uGdJJkDkEPi9A1MY45G97i5EM7uOPVVXzx/Z6Obd/LsMbCGiShC8M4mjX3LFjDpxt2kdFrCSmpG4wPIpJJLQTDIqQ0NzJ3w1vk9kjmlnkrfc5VepSnYY0lqVuKJHRhCK01D7y3joVrtvHbUSeQmr7G+CBCldRKSjrcAxYIybAIAOlNR3hl0hlkpCQw/qUV/Hdnre8XBdh+L6KLJPQoFszs9Eb76+LvsC+v5Nbzj+OOgkGRCcIMSc1mgxEj4Pzzoby8w8m8Ra/ULtgnjSAhNoZxs5ezZV+d9xeEeFhjYU6S0EXYzVz2Pc+VfM/1w7OZOurEyAVisaSW3SOZVyadwZHGZmwvLmdXzREvK4duWGNhXpLQRVi9urySRz7cwKWn9OFPV54c2VET/Uhq9jI7pVWlLKtYRu6MXOxlHWuOCfb1HXVC727MmzicvQfrGfl8EaWpB1nWvfrYbYdoWGNhbpLQhUfBNukcihnB9HfLKDghkyeuzSM2JsJD4PpIavYyO4WLCql3ONvYK6orKFxU6HdSDvb1gcob0J1rzi3nm8N/pT7WAcrNtkPUfi/MTRK6CIu6mFPZHX87w3IyeN42lIS4toeaz4kawsFHUpu+dDp1jW3bousa65i+1L829mBfH4y53zyCVvXetx2i9nthXjLaogi5FZv3sTt+Mgl6Cy9OGE1SgolmGvIyVntltfu2dE/L/V3P39dD4CNBhmLbIvrJFboIqXVbq5k090ti9R56NTxKapfomQM0O819G7un5aF+fTA8bWNA2oCwb1uYhyR0ETKbdh3kxjkrSE2Kp/fghcTmeZ4b04yKRhaRHN+2jT05Ppmikf7dOAz29cEoShxNcmPbZUkNcOX2M9Fah337whwkoYuQqNpfx7jZy4lR8PebzyAu3o/OLiZjG2Kj+LJiEmOdbew5aTkUX1aMbYh/bc3Bvt4nL71cbY9+QPFCyDkASjv/nbUIfjPzU578ZGNoti9MT9rQTa6lysTwG4gdsLu2nnGzV3CwvonXC89kYM+ukQ4pYLYhNmZ95WxjD2SfB/t6jzz1cgXnfYHKSmwabGVtX6bVHp5eupG0pHgmnTPQ93ZaJqAOpCesiDi5QheeuWbPcavVRA3NOTkM/fwD5t40jJP6phobY2fhq5erpxr7AQO45OTePPT+ehas3BLeGEXESUKPUkZ3YGm7cTv2J29ixCRN3O/hjPE7ObViBkP//WHkY7MqX71cPdTYqz//mRlj8zh3cE+mvrWWDzN+Et44RURJQo9CfndgKSj48St0KLf/4t0UXtxIRXfQCiq6w22XNGJ/8e6Ida6xPF+9XF019ju6KZqhTY19YlwsL4wbymnZ6fy/QZfyWVquQUELo0lCj0KR7MACMC1vL3UJbZfVJcD0vL2GxRZNA5OFhD9d9202xk5M5YK7047pOJScEMec8cMYdHgvt5xwJV9s6uBY6iIqSEKPQpHsRNLoaGZLmvvnKtOkg0tQvI3VHoKu+2nJ8fz92wXkHDnApHkrKf1hb8fiC9M3PhE6ktCjUKQ6sDQ5mrn7ta9JaejmfvvxPSLauSaq+TNWewi67vdoOoz92wX0S09i4twv+bJ8X2jiF6YgCT0KRaIDS5OjmSkL1vBB2Q6uPe0PJKu2bS7JKoGiy5+KaOeaqGbgWO2ZjXW8essZ9E7rwoQ5K/iqYn/ItyEiQxJ6FDragcWhQIehA0s7jmbNvW+sYdGabdx3yYm8eM0UisfMabv9MXOwDbGFv3ONVRk8VntWty7Mv2UEmd0SmTBnBav9ncpOmJok9ChlG2JjRE0q5x9Io3xyedgSpkbx2zfX8u7qbfzm4hO49fzjfW7fNsTGiP4jOD/nfM+xSXtsWxGYgKJXahfmF44gvWsC42Yvp6yqOmzbEsYIKqErpcqVUmVKqdVKqZWhCko4RbKe215mpzS1hs+6H+Cpb37Beadu4M6fR2jquM4gQhNQ9ElLYn7hCNKS4rlh9nK+Sc4K6/ZEeIXiCv3nWus8rXV+CN5LuESynvvotmM1KHDE7OatzQ9ILXk4GTUBhZsJrvt1T2L+LSPomhDLDT+9xpJJvbOUuUqTi0lFstZ82ifTIlrn3mlFcAKKARnJzC8cQVJzE2NPuo4Vm6X6JRoFm9A18LFS6iulVKG7FZRShUqplUqplbt37w5yc51HpOq56xqaqKxxP+aHkbXklr6icnOVbAY5Pbryxjevktl4iHGzl/Pphp2RDkl0ULAJ/Wyt9enAJcCdSqnz2q+gtS7WWudrrfMzMzOD3FznEYl67gN1Ddzw4nJim3savm1hDv0aannjm/n8pFc3bnn5K975usrv11r6j3CUCCqha623uf7dBbwDDA9FUML4WvOdNUe47oVS1m2t4ddn/F5qyTuxHk2HmV84guG5GUx5fQ0vfb7ZmA1L5VPQAk7oSqmuSqluLb8DFwHrQhVYZxd0rbm3buTtlO85xNUzv6Bqfx0v3TSMv4y+k+L08fSvdk2WcDCW4vTxUktuda2OmZQTBjGvy0YuOqkXf1i0nieW/BeZ9ygIBv2xCmaCi17AO0qplvd5VWv9UUiiEoBrsoTZdwJQ8mC5/y/0NRlCK+u31XDjnBU4mpt59ZYRnDqgO9jt2O6dh+3ofVEHJM+D1LNlpnircnPMJNx+G8/PfIH78ofw9NKN7M+9kD+UfxK2SoqCPOfY+yUeVyhw/mvC+w9mEXBC11r/AJwawlhEqHjrRt4qIa+4zMakE64ipUd3Xis8i0FZKR16vbAQD5957O/u59HNm0lPTuCFz+BAXBcea3KQGBcbmTgD0NKfo95RT+6MXIpGFln226ZMQRfFSlbnuX/Cj27k81dU8sBPr2NA/QFeuX00/bondej1wpw8HhO+ePnMlVLcN/qnpL/yEo/knM/W4lJm3jCUrNQugQdqEE/9OYC2Sd0iV/9Shx6tvLWRe+lG3tDUzLR3yrjv7TLOrKnknXX2tsncx+stwaRlgxHl6zO327ntnafZ/OilPDPtKp6Z9AdWVbYb1MvblIUREum5A4wmCd3MWs3b2SZp+xpq1UM38urfPcgvZ5Xy6vJKbjv/eF7a8BZpriuXNiLUDV1EkLfPvNXxpoB+1buY9t6T/P1XD/P6lx341haBKpbONj6/JHSzcp1EvWu180NqnbR9DbXqpht5+cMzuHjnANZtq+aZ609j6iUnEuupbsHLdGbC5AL99uHtM3dzvCU11jPt36/wf2+V8bt319HQ1ByC4EPPFOPzd6DiLFjShm5W3pK2P23cNhvMmgXAgsde4f5315HVTfH27WdzUt9U39u32Rj7pavCZkZ5AP+BzqtkQkmkQwiMp8/cw/HWY/9OCs87juLPfuC7HbU4SCWWGgMC9V/RyCIKFxW2aXYxtE9FByrOQiE6rtB9fVWL5g4JnmL3lrT9bONuVDE8mHMBv31zLfk56Sy86xz/knknVzKhJHqTcjh4ON5UdjbTRv+Up8bmsXbrAbYl/ol6NdDj2xTkrT5amngMT82LQYr4+PwGTlwC0ZLQOyNvSduPNu5Nu2q5+me/ZG6fodx8zkBenjicjK7tZnYOUsmMA5TMkIkRLM/H8XZFXj/evO0sQLM94fc8X/I9TY4ONMF4a15svY6vZgsPF0c+x+cPZ5OIwRVjktDNyttJ5GWoVUez5oVl3zP66X9TmZjGs/99j/svPYm4WPmoLSESFTp+DO17cr80+tZPJ7n5ax79aANXz/wP3+8+6N/7+7qK9We+1UCF873B8IoxOct9iVRzjq+TyM1Qqz/sPsg1M7/g4Q83UPCTTD5e+xL/s++/xsdOeJstIjnxR6flY2hfe5mdL1O38W3XP1KXcSur9ixi9FP/4sV//YCj2cegAZWV2IdA7mSI+b3zX/sQfryKDWezRaje21OeMLhiTG6KBqlldLmwJK9WNza9XZU1N2vmflHOXxZvICE2hhnX5XFFXl/UnDqPr4lWfncUEYZpMyEKsPvwVpLinmF4z6786R/NLP5mB3+92nOncvv5GRSetZc6V4tgRXcovAzokYENwttsEe4mkZY/fJMmOb8F5OT8+C07DKL/Cj0E7V/ehv30eiPHBCoT0xg7q5Q/vr+eM4/rwZJfn8+Vp/XDNcaOOQXxmXW2jiJm4ulccPeZHG6qY3Pjizx+zals2FHLqKc+oyb2IjTHHpfTL+RoMm9Rl+BcDoS32cKIJhEDJy4xf0K329mx6jOa3Z384W7/MrG6hib+1ncEo06ZwLfbavjL1acwZ8Iwepm9O3aQn1ln6ygSDTzt+y3VW/jfof35eMp5nDGwB/vix7Mj4X6+btfDtLLJ/exIR5eHs9nCYp3ozJ3Qfd39NrgkyAyaHM28tqKSgr+W8Fj2uZxdXcHiKedxbf6AY67KTfntIsjPzBQdRUQbvj6TPmlJzL1pGD0ai2lUfRjz3BfcYf+KzXsO+fX6UHR083hPx6i5XA1i7oTu6+TvRINIaa1Zsn4no576F1PfLqN/ehJvrnuVWf99l77tx2IJkZLVeYEP9uRJkJ+Z0RN/CN+KRhaRrNq2mSSrhDafiVKKbo5l9K//NXePHEzJd7v5xRPLeOC9dUw9+w++P1ObjbETU7ng7rTAmi28FTdEcC7XUDN3Qvd18oeq/cvToEJh6uzQ+v39aUteldKHa1/4D7e8vJLmZs3MG4by1u1nkX9wa9i3HXJBfmYR7ygSblE4cJhtLRQv1OQccE2IcsD52La21Uquc6nkqZ1MueUi/pOzg+uGDcC+vJJnFvXiqtw/kuiICWwyF9f7R+R4jvS22zF3lUt2trOZxd1ycLZzFRa2vYoPVftXS3NPnavkKtRddn10CdZa8/WWAxQPvoKPevyEnnvq+NOVJ3PdsAHE+1tTnufh6trg7shthOAzsw2xMesrZ/WP9Og0genTsVU0Yvuq9cLGH8fPd3Mupf2/OykqLmbilCv4y0cbWLz2RPrpeaQ2vc+yqa+Q2iXe/+1H8nj2d9sG/ZE29xV6URH2ofFt61OHxv948oeg/cteZqc0tYZl3avb1jRPn479+Lq22z6+bVtvUPXQHpqTmqdN442VW7j8b59z1XNf8HlaDpO3fM6y3xRww4gc/5N5ANs25N6DH5+Z1JlHhsdzwRdf36S9nEvHZ6bwwrh83rr9TOL0DvbHj2PEn5dy/7tlbNxZ619skTyeTXYfz9RX6PZToPByxdE/7N2djzkFjp7+ftZqu33/dvWzrWuaSa2g8DKOrY1dVIGt9Wt91EN7rFP3eBJs4TdvrmVwVgoPXXkyV027ma7NjZAYwo8q0vcevHxmUmceGd7OBZ/73cc3abuPcwlgaE4GfRr+RL3K5eyh81iwsoq/l1Zy1vE9GNB/FX/7+jeeY4vk8Rzpc6kdU1+hT186nTrd0GZZnW4IWc2xt5rm6RfHuq+NvTjW52v94qHNeH+PXrx68xl8POU8xo3IoeunSzz/oQq0vdXEE1hInXlkBLXffZT++TqXWkvU5Tx2zan8Z+oF/HbUCZTvOcTjpQ96j82A49ljXxWTnUumTujhrjn29v6VKQ73z7mW+x1buxuuRxodfLJ+J69efiuH4xPbrNqclESPpx7jrEE9w9sxyN/a2wjcoJM688gIar/7KCv0dS650yMlkTsKBvHZb39Oc8we77GFopY80GPdz21767wYSqZO6OGuOfb2/tlpOR6ey+lwbP1MX14AAA6WSURBVM10YdGabdz56iqGPrSEm19eycMZp/PmrQ9wOCnFOc1ETg4xs2YZUzJl4gkspM48MoLe717KCn2dS97ExcZ4jCGmuSejZnzGk1nDqHr0Kc/HczirUEx2Lpk6oYe75rgocTTJjW2XJTc6l/vatrfnD9U38fmmPcz45L/siJ9KZeJMfjX/a5b/sJfL8/oxb+Jwvrr/F4x7ZhpJw4eiIlH/Gmxdb5hInXlk+LXf8/I8V04F+95e3t/d65Pikrjp5Gmkdonn6U83ck5VH865bTan3HMvC978F5svvhKttTG9yU10Lpn6pmjLzZhJ702i3lFPTloORSOLjr1JE2CzgO3RDyAVpo+EyjTIroaipWCr+QBuf8657TfGUR+jyenedttHY3M9n5ncj3OzfsVrJTk8MP9jHM0apSBOpZLqWMLsux7i9Ox0YmNMPMaKCfj9mYuQCud+b3+utD+Xgo1td209n3y7kz++MZe62Hx++6azAL5nSgKLZ9xLD09VKCa5iAklUyd0CHPNcWUlNg22snbLVeWP257tmpLrwXL2Hqyn9Ie9bNx1kE07a9m463j61b+MQ6XDEfimNoa8AbHcUXA8Q3PSOT0nncunZgIwLPfZ0MZuYVJnHhnh3O/tz6WAXu8htsxuiVw/PJsXXn0CjeLF+yr5snw/Kyv2kb53p9v305WVfFi2ncFZKeT06EpCnKkbK/xm+oQeVh7KrQ736cdbpRVsrz7M7vhbaVJZnP7QEvYd+rHipmtCLBMr/8Njb8+gd20NTf36E/Pww8SNu+ToOi21s/UxmtwZuW6vSlrGWikJ8X+tpZa73lHvcdtC+MuI4ymYPyKtz7VfzB9C0cginrjWBr91f45v7daTO+yrAIiLUeT27MrAnl3pm9aFPt2T6JPWhT5pzn8/rXyb0h8+83ge+zrPjTwXLZfQGx3N1NU7ONTQRF1DEwfrHdTVN3HgcCP76xo4UNfI/kMN7K9r5ISLJ3HjS0V0aaw/+vq6uESmnn4dC99d52weiTmJOL2Xi07qxaCsFAb36sbgrBT6/ONt1OOPHe1UkLC1Cm67FWIU2GzB1fUGSWq5RSiZ/Xjyeq556Jmc+fTjvH/hOWzcVcvGnQfZuOsgFXsPUfrDXmqPNB1d9WDsP9kX/zd0q/ee8M7NvPHlFs7tP4YN1R/yynfT2mz7loW3UHu4iRvzbLyz4XUK3zdu30VFQq89MITDB4/jppdW0OBopr6xmfqmZhqamqlvctDQ1MzhRgeH6h00+DGXYVJ8LOnJ8aw/qYCGG5sZ88aT9Kmp5lCvvmy6exoTbL/kvrQuZKYkcuE9GQA88r+/bvsm3nqI2Wxe63rDfRJEctvCevw9nkI+kJufvMY3udy5oN0EE4k2GyfjnDoPXB0A06FsSgm1RxrZUX2EbdVHuPqdW9FH6tu8d5M+wgdbnmTtxiFUJvwZR8zhNs8fbjrMXf/4DY+83ZOqxCk4Yow7F6MioTuaUmhsSGfPwQYS42JIiIuhW5c4EuJiSIyL5dPyxaiEJiYOu5quCbEkJ8aRkhhLckIcXRNjmfbpPcTE1PP29a/QPTmeLvGtOzScS0HKk0AaJTO2clr7jXu6q++jh1gk66mllluEktmPJ5/xdbA3ebcu8XTrEs/gXt3Ye2Sb23Wa2M2motHEPeS+Rt4Rs4epl5zI7SU+auhDLCoS+up7H/b6fMHcKQBMveR+t893Ka0CoHea+8kfArqy8NHdOTsug4qmvcc+HZfR8W11UHZaNhXVx8YmtdwiEGY/nsIZn7f3jolRHp/PScvmtvOP55Gvjd131ri1G0YeB8b30UOs6BNIbmj3dINzeRsB1vZ6I7XcIpTMfjyFM75g+qOEOzZ3JKEHyseogbZl+yheRNsxohc5l4c9NKuPGS4MZfbjKZzx2YbYKE4fT/9q13l8MJbi9PFt+qN427bR+y4qmlyC5m7yilDw1jaXnY2trOLYGvccY76mSi23CKWQHE9BfhP1tt2wHe92O7Z752E7el/TAcnzIPXsHy/efGzbyHNRrtAhPINQWWzyWSF8isLZlnwy2XjnvgSV0JVSo5RS3ymlNimlpoYqqI7wNRlCwIP2B8uPQXsiOZFDxPaL6LQ83o8yZOMlHv/YeD0PTTbeuS8BN7kopWKBZ4FfAFXAl0qphVrr9aEKzhdfHR4i2bnHuREbY790dXeeUd6h2MMp4vtFCJPweR76mgbTZIK5Qh8ObNJa/6C1bgBeA64ITVj+8TUov5knS4hkbGbeL0IYyee5EGVNp8HcFO0HbGn1uAo4o/1KSqlCoBAgO8R/1Xx1KDBzhwjpeCRE5PnVKQmO6WnafqRGX01JRjU1BZPQ3Y0Dq49ZoHUxUAyQn59/zPPB8NWhwJAOEQHeBIpkZw1TdBQJ8uaZVO6IUPDrXAhi3mKjBdPkUgUMaPW4P+C+n2yY+Czq9zKBhVFKVue57Ykayc4aZu8oIiInojcugxRI7FY7F4JJ6F8Cg5VSA5VSCcBYYGFowvKPz6L+Rz+geGG7zj0LXRNbRFgkO2uYvaOIEEax2rkQcJOL1rpJKXUXsBiIBeZorb8JWWR+8lq072MCi0iLZOcf6XgkhJOVzoWgeopqrT8AIn+564kZSo5M3uYmhLAOa/cUjbKSIyGESUVJL1hrj+XiKi3acfs4smo1MR5Kjqwq2r8+CiE6xhIJ3Wvi8tJbU5iX/DESouMskdB9idTUWEJYiZX/yFrl/2btNnQhhOhEOsUVuplZ5cpACBF5coUuhBAW0Tmu0KOg3EgIIYIlV+hCCGERktCFEMIiJKELIYRFdI42dOGWVNgIYS1yhS6EEBYhCV0IISxCEroQQliEJHQhhLAISehCCGERktCFEMIiJKELIYRFSEIXQgiLkIQuhBAWobTWxm1Mqd1ARYAv7wnsCWE4oSSxBUZiC4zEFphoji1Ha53p600MTejBUEqt1FrnRzoOdyS2wEhsgZHYAtMZYpMmFyGEsAhJ6EIIYRHRlNCLIx2AFxJbYCS2wEhsgbF8bFHThi6EEMK7aLpCF0II4YUkdCGEsAhTJXSl1DVKqW+UUs1Kqfx2z92nlNqklPpOKXWxh9cPVEotV0ptVEq9rpRKCFOcryulVrt+ypVSqz2sV66UKnOttzIcsbjZ5oNKqa2t4hvtYb1Rrn25SSk11aDY/qqU2qCUWquUekcp1d3DeobtN1/7QSmV6Pq8N7mOrdxwxtNquwOUUv9USn3rOifudrNOgVKqutVn/YARsbm27fUzUk5Pu/bbWqXU6QbFdUKr/bFaKVWjlJrcbh3D9ptSao5SapdSal2rZRlKqSWuPLVEKZXu4bXjXetsVEqN92uDWmvT/AA/BU4ASoD8VstPAtYAicBA4Hsg1s3rFwBjXb/PBG43IObHgQc8PFcO9DR4Hz4I3OtjnVjXPjwOSHDt25MMiO0iIM71+6PAo5Hcb/7sB+AOYKbr97HA6wZ9jn2A012/dwP+6ya2AuB9I48vfz8jYDTwIaCAEcDyCMQYC+zA2SknIvsNOA84HVjXatlfgKmu36e6Ow+ADOAH17/prt/TfW3PVFfoWutvtdbfuXnqCuA1rXW91nozsAkY3noFpZQCLgDedC2aB1wZznhd27wWmB/O7YTBcGCT1voHrXUD8BrOfRxWWuuPtdZNroelQP9wb9MHf/bDFTiPJXAeWyNdn3tYaa23a61XuX6vBb4F+oV7uyF0BfCydioFuiul+hgcw0jge611oL3Tg6a1/gzY125x62PKU566GFiitd6ntd4PLAFG+dqeqRK6F/2ALa0eV3Hswd0DONAqYbhbJ9TOBXZqrTd6eF4DHyulvlJKFYY5ltbucn3NnePh65w/+zPcJuK8gnPHqP3mz344uo7r2KrGeawZxtXMcxqw3M3TZyql1iilPlRK/czAsHx9RmY4xsbi+WIrUvsNoJfWejs4/3ADWW7WCWj/xYUkvA5QSn0C9Hbz1HSt9XueXuZmWft6S3/W8ZufcV6P96vzs7XW25RSWcASpdQG11/soHiLDXgeeAjn//0hnE1CE9u/hZvXhqR+1Z/9ppSaDjQBdg9vE5b95i5cN8vCelx1lFIqBXgLmKy1rmn39CqczQkHXfdK3gUGGxSar88o0vstAbgcuM/N05Hcb/4KaP8ZntC11hcG8LIqYECrx/2Bbe3W2YPza12c60rK3Tp+8xWnUioOuAoY6uU9trn+3aWUegfnV/ygE5O/+1ApNQt4381T/uzPgPix38YDlwIjtaux0M17hGW/ueHPfmhZp8r1madx7FfosFBKxeNM5nat9dvtn2+d4LXWHyilnlNK9dRah30AKj8+o7AdY366BFiltd7Z/olI7jeXnUqpPlrr7a5mqF1u1qnC2dbfoj/Oe4teRUuTy0JgrKviYCDOv6YrWq/gSg7/BK52LRoPeLriD4ULgQ1a6yp3TyqluiqlurX8jvOG4Dp364ZSu3bKMR62+SUwWDmrghJwfjVdaEBso4D/Ay7XWtd5WMfI/ebPfliI81gC57H1qac/RKHkaqefDXyrtX7Cwzq9W9rzlVLDcZ7Pew2IzZ/PaCFwo6vaZQRQ3dLMYBCP354jtd9aaX1MecpTi4GLlFLprmbTi1zLvDPiTm8H7giPwfmXqR7YCSxu9dx0nBUJ3wGXtFr+AdDX9ftxOBP9JuANIDGMsc4Fbmu3rC/wQatY1rh+vsHZ5GDEPnwFKAPWug6cPu1jcz0ejbNy4nsDY9uEs11wtetnZvvYjN5v7vYD8Eecf3QAuriOpU2uY+s4g/bVOTi/Yq9ttb9GA7e1HHfAXa59tAbnTeazDIrN7WfULjYFPOvar2W0qlozIL5knAk6rdWyiOw3nH9UtgONrtw2Cec9mKXARte/Ga5184EXW712ouu42wTc5M/2pOu/EEJYRLQ0uQghhPBBEroQQliEJHQhhLAISehCCGERktCFEMIiJKELIYRFSEIXQgiL+P9Bhj0WXKLMQQAAAABJRU5ErkJggg==\n", "text/plain": [ "