{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Restrictions via penalties\n", "===============================" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Parameter restrictions: {'A': (None, 10), 'sig': (4.1, 7.1)}\n", "\n", "Penalty factor: 1e+20\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZwU5bX/8c9hURQ3FBJlEVwwronLiI4anYgLEBWjMgO5mhslkpi4XTWJS+JCov68iSYmaq4aE2JidBBciHKFKAJRB3VwQcHlIkZBUUBwQWQ/vz+eaqdoe2Z6hu6uXr7v12tevdVUne6uOl311HPqMXdHRERKX4ekAxARkdxQQhcRKRNK6CIiZUIJXUSkTCihi4iUCSV0EZEyoYQOmJmb2a4FXqaZ2Z/NbJmZPdPG/+0XxdypmdcvNbM/5ibSgsz3y2Y23cw+MbPrcz3/Vpa93Mx2ztO8rzWz89v4P/9hZpNbeL3GzBZsfHSFY2bfNbMnko6jkMxsvJkNLvRyyyKhm9kjZjY6w/NDzey95hJfwg4DjgZ6u/uAXM7Y3a9x9+9tzDwyJY5czLcZo4AlwFbufmEe5g+AmU01sw3id/ct3H1eHpbVA/gOcGtb/s/d73L3Y2Lz2aidjeg917T3//MhiR2oQjCzeFHPdcAvCx1DWSR04C/AqWZmac+fBtzl7msTiKk1fYF/u/unSQdSBPoCc7y8qty+C0x098+SDqScFenOGu7+DLCVmVUVesEl/wdsBnwEHB57rhuwEvgaMABoAD4EFgI3AZvEpnVg1+j+VOB7sde+CzwRe7w78E9gKfAaUNtCXD2BCdG0c4Ezo+dHRrGtA5YDV2X43w7Az4C3gEXAncDW0Wv9ophHAe9G7+mi2P9eCfwt9vhg4Kno/b8I1MRe2xb4czSfZcADQFfgM2B9FN/y6L18Pl/gf4Gz02J+ETipLZ8TMAZYA6yOlnNU9NwvY9PUAAtij/8NXATMir73eqBL7PWhwAvAx8AbwCDg6ujzXhkt56YM3/3W0ee8OPrcfwZ0iK8HwK+jz+lNYHAL3/0U4NTY42nAydH9Q6PlfjN6PBB4IX19A6ZH030axVyX+iyAC6P1YiFwegtxTE1939H3Ny76vD4BngO+lra+jo/e/5vAuWnr1Njo8/kEmA1UxV6/OPqsPwHmAN/KtA01855eBo6PTd+ZcMS2X4b3k3r/PwXeA/5K2NYfiuJeFt3vnfYZ/JKwDSwH/gFsB9wVrSPPAv3S8sG5wLwojl/F1oNdo+/yo+i1+vj/pcV6O3BFQXNhIReW1zcSPrw/xh5/P7aRHEBIap0IyfAV4Py0L7DVhE5IdPOB06N57Rd9qXs2E9N04BagC7BvtMIdmT7fZv73DMKPwM7AFsB9wF+j1/pFMd8dxbRPNO+jYhtfKvH2Aj4AhhB+JI6OHveIXn+YsIF3izakI+IbTlpM8fl+B3gy9tqehB+MTdvxOY1hwwSe/niDWAgJ/RlCAto2+j5/EL02gLCxHR29317A7pm+2wzf/Z3Ag8CW0Wf8OjAy9n2tAc4EOgJnEX4ErZn3tBg4MPZ4NPD76P6lhOR3Xey1GzOtF/H4Yp/F2uh/Okff6wqgWxbbyJXRezgl+t+LCIm7c/RZzQQuBzYhrHfzgGNj/7syWl5H4FpgRmzew6LvowMhSX8K7JDle/oJGybGocBLzbyH1Pu/jrCubUZIzicDm0ff3b3AA7H/mUrYlnYh/GjPib7bowjr553An9Pie5ywbu0YTfu96LW7gcui99kFOKyFz/sC4L5C5sFEm1zM7E9mtsjMXs5i2gvMbI6ZzTKzx8ysb+y1/yTshZ0eayP9DqEpBnef6e4z3H2tu/+b0K55RDtCPo7QTPLnaF7PE/ZohmWItw9hT+yn7r7S3V8A/hjFlY3/AG5w93nuvhy4BBiedoh5lbt/6u4vEfayR2SYz6mEQ/+J7r7e3f8JNAJDzGwHYDAhGS5z9zXuPi3L+O4H9o19D/9BWHlX0YbPaSP8zt3fdfelhD2ufaPnRwJ/cvd/Ru/3HXd/tbWZmVlHYDhwibt/Eq0n1xOa7VLecvfb3X0dYd3aAfhyM7PchrC3mjKNpnXucEJCTD0+Ino9W2uA0dH3NZGw1/mVLP93pruPc/c1wA2EpHQwcCDhR360u6/2cF7hdsJnkvJEtB6tI+wZfy31grvfG30f6929Hvg/wo9rNv5GWB+3ih6fFs2/OesJe76r3P0zd//A3ce7+wp3/4RwNJa+ff/Z3d9w948IR5dvuPujHppj7yXsdMRd5+5L3f1t4Lc0bVtrCE2EPaPtuqWTvZ8Q1oOCSboNfQwhEWfjecIh3lcJh43/DWBm2wJXEFaut4CrzWw/wsr092ia3czsoegE6cfANUD3dsTbFzjIzD5M/RES2fYZpu0JLI1WsJS3CHuM2egZTR//305smEDmp73es5mYh6XFfBghGfWJYlyWZUyfi97XwzRt8CMIh7CpZWb7ObXXe7H7KwhHMRDe0xvtmF93wp5q+mce/74+X6a7r4jubkFmywh7iykNwG5m9mXCj8+dQB8z605YV6e3IdYPfMPzQvH335rP1xl3X09ovuhJlKTSvrNL2XB9S//Mu6R2MMzsO2b2Qux/9ybLbczd3wWeBE42s20IOxl3tfAvi919ZeqBmW1uZrea2VvR9j0d2Cb6kU55P3b/swyP0z+/5ratnwAGPGNms83sjBbi3JJw1FowiSZ0d59OaGP9nJntEvVamWlm/zKz3aNpH49tRDOA3tH9Y4F/RntqY4BVhBVxkrunvrQ/AK8C/d19q+j19BOoKZ8SDt1S4kloPjDN3beJ/W3h7mdlmM+7wLZmFt+odwTeaWa5mf6/b+zxjoRDzfiK2Cft9XczzGc+oakmHnNXd/9/0WvbRhtRumxOUN4NjDCzasKe3uOxZWb7OWXS0nfQmvmEQ+tMWnpPS2ja+0ppy/eVbhaw2+cLDuvuTOA84GV3X01o072AsLe4pJ3LaavP1xkz60DYjt4lfG5vpn1nW7r7kNZmGB2l3Q6cDWzn7tsQ2sWb28Yy+QvhaHIY0ODuLX3u6d/jhYQjlIOi7fvwVGhtWH66jNuWu7/n7me6e09Cs+4tLfTY2YNwXqlgkt5Dz+Q24Bx3P4DQxndLhmlGEg6bIOxBpX5N7yT8kh5N1NwS2ZJw8mN59APRUmJ5ATgp+tXfNVpWykOEvazTzKxz9Hegme2RPhN3n0/YYK81sy5m9tVoXn9r6c3H3A38l5ntZGZbEI4q6tP2zH4exbkXob26PsN8/gYcb2bHmlnHKJYaM+vt7gsJn+MtZtYtej+pjeF9YDsz27qFGCcSEuDoKLb10fNZf07NeIFwCL6tmW0PtKUv9x2EpreBZtbBzHqldgqi95Sxz3nUjDCWcIS3ZZSkLiD77yvdRL542D+NkPRSzStT0x5n0mzM7XSAmZ0U7VmfT9gBmkE4J/GJmf3UzDaL1pW9zezALObZlZBkFwOY2emEPfTmZHpPDwD7E37w7mzTOwrb92fAh7Ej9o3142ib6BPFVA9gZsPMLLUzuYzwvtc3M48jaMpTBVFUCT1KXIcA95rZC4S27h3SpjkVqCKced5A1O45n3BSZ0LspYuAbxPatG4nc+JL+Q2hx8X7hB+Fzw/9omaGYwjNDO8SDkFTJ2cyGUE4ufYuoc35Cnd/tIVlx/2J0I44nXDiaiVwTto00wgnex4Dfu3uXyhIiX5YhhKOShYTPp8f0/Tdn0bYM32V0Gvi/Oj/XiX8qMyLDqO/0JwTtZffRzi59PfY8239nNL9lbBn829gMi1/X+kxPUP4cfsN4eToNJr2um8ETrFQzPW7DP9+DuHoYB6hR8vfCd9De9xJ+FHaLPbcNELymd7M40yuBP4SfQe17Ywl7kHCSctlhO/+pKgtfh3h3Me+hPVtCeGcT0s/6AC4+xzC+YYGwnazD6EJpTlXkvaePHTvHA/sRFin2uK3hJOjSwg/To+08f8zeZBwRPUCoWnxjuj5A4GnzWw5Icec5xnqGKIfwuXR+lgw5p7NkXUeAzDrBzzk7ntHJ0Vec/cdmpn2KOD3hJ4Yi6LnRhC6ZX0/enwrMNXd7y5E/MXIQpFVb3dvqX1P8szMrgEWuftvk44FwMyuJPQuOTXpWDIxs8uB3ZKOz0KBUH93n7sR8xgP3BGdtC6YotpDd/ePgTfNbBh8Xh7/tej+foQ99hNSyTwyCTgmOjzqRtgznFTg0IuGmRmhC+GbScdS6dz90mJJ5sUuaioZSWhyLXnufnKhkzkknNDN7G7CYdpXzGyBmY0k9IYYaWYvEooXhkaT/4pwJvre6Gz6BIDoZOgvCMUBzxK6cy2lcj1HONF1e9KBiGTDzM4kNAX+b9RRQtop8SYXERHJjaJqchERkfZL7MI23bt39379+iW1eBGRkjRz5swl7t4j02uJJfR+/frR2NiY1OJFREqSmb3V3GtqchERKRNK6CIiZUIJXUSkTCihi4iUCSV0EZEyoYQuIlImSi6hNzTAtdeGWxERaVKUI2Y3p6EBvvENWL0aunSBxx6D6uqkoxIRKQ4ltYc+dWpI5u7hdurUpCMSESkeJZXQa2pg02iIBLPwWEREgpJK6NXVMGUK7Lwz7LADHHxw0hGJiBSPkkroEJL6T38K8+fDiwUdflVEpLiVXEIHOOkk6NgR6rMeaVJEpPyVZELv3h0GDoSxY8MJUhERKdGEDlBXB/PmwcyZSUciIlIcSjahf+tbodnlggtUZCQiAiWc0F99NTS3/OtfoflFSV1EKl3JJvR4UdGqVSoyEhEp2YQeLzLq0EFFRiIiJZvQq6vDtVz22AO22QYOOijpiEREklWyCR1CUv/Zz2DJEnjyyaSjERFJVkkndIDjjw9XXhw7NulIRESSVfIJfcstYcgQGDcO1q1LOhoRkeS0mtDN7E9mtsjMXm7mdTOz35nZXDObZWb75z7MltXVwXvvwfTphV6yiEjxyGYPfQwwqIXXBwP9o79RwB82Pqy2+eY3Q4+Xiy9Wf3QRqVytJnR3nw4sbWGSocCdHswAtjGzHXIVYDZmzYK1a+GZZ1RkJCKVKxdt6L2A+bHHC6LnvsDMRplZo5k1Ll68OAeLDqZObbpIl4qMRKRSFfSkqLvf5u5V7l7Vo0ePnM1XIxmJiOQmob8D9Ik97h09VzCpIqOvfQ023xyqqgq5dBGR4pCLhD4B+E7U2+Vg4CN3X5iD+bZJdTWMHg2ffAKPPlropYuIJK9TaxOY2d1ADdDdzBYAVwCdAdz9f4CJwBBgLrACOD1fwbbm2GNh663DSEaDBycVhYhIMlpN6O4+opXXHfhRziLaCJtuCieeCA88EE6OptrVRUQqQclXiqarq4OPPoLJk5OORESksFrdQy81Rx0F3brBTTfByy+HHi/V1UlHJSKSf2WX0Dt3hkMPhYceCj1fNtkk3Cqpi0i5K7smF4BUF/d162D1ahUaiUhlKMuEfnrUz8Ys7KGr0EhEKkFZJvSvfz30dunUKTS9qLlFRCpBWSZ0gHPPhTVr4IMPko5ERKQwyjahH344bL99KDISEakEZZvQO3aEU06Bhx8OlwMQESl3ZZvQAWprYeVK+Mc/ko5ERCT/yjqhH3oo9OqlAaRFpHg0NMC11+ZnIJ6yKyyK69ABhg0LVaNXXAGDBqnHi4gk54knQjdq93CtqVwXPZb1HjrAHnuE4el++UsNTyciybrzzlDwuH59fooeyz6hL1kSbvP1AYqIZCvVjbpjx/wUPZZ1kwvAN74RCozWrlXVqIgkZ+3a0ORy5JHhIoL5uHBg2e+hV1fD//xPuH/++WpDF5FkTJsGixbBD38Il1ySn1xU9gkd4IwzYOedYebMpCMRkUpVXw9du+Z3NLWKSOhmoU/6Y481tamLiBTKmjUwfjyccEIYyD5fKiKhQ0jo69bBffclHYmIVJopU2Dp0jCiWj5VTELfd1/o31/XdhGRwvv970O/8623zu9yKiahm4Vfx8cfh8suU390ESmM6dPDNaVWrYIhQ/KbeyomoQN85SuhQuvaa1VkJCKFcccdTffzXQtTUQn97bfDrbuKjESkMBYtCrf5KiaKK/vCorh4kVHnzioyEpH8WrkSnnwSjjsODjkkP8VEcRWV0Kurw7UUvv1t+MEPVGQkIvn1yCNhPIazz4Zjj83/8iqqyQVgxAjYZx945pmkIxGRcjd2LGy3XSj3L4SKS+gQers89RTMn590JCJSrlasgAkT4OSTQxNvIWSV0M1skJm9ZmZzzeziDK/vaGaPm9nzZjbLzIbkPtTcSXXuv/feZOMQkfI1cSJ8+mkoaiyUVhO6mXUEbgYGA3sCI8xsz7TJfgaMdff9gOHALbkONJd23RX2319FRiKSP7fcEq7dsskmhVtmNnvoA4C57j7P3VcD9wBD06ZxYKvo/tbAu7kLMT9qa0M7+o9/rP7oIpJbjz0WihhXrAgnQwuVY7JJ6L2AeGvzgui5uCuBU81sATAROCfTjMxslJk1mlnj4sWL2xFu7uy6a7i9/noVGYlIbt1+e7gtdM1Lrk6KjgDGuHtvYAjwVzP7wrzd/TZ3r3L3qh49euRo0e3z+uupmFRkJCK5tXBhuC1EMVFcNv3Q3wH6xB73jp6LGwkMAnD3BjPrAnQHFuUiyHyoqQlnntesUZGRiOTOxx/D00+HAer32y//xURx2eyhPwv0N7OdzGwTwknPCWnTvA0MBDCzPYAuQLJtKq2org59RAG+8x0VGYlIbkyYEC7Edf75+RuZqDmtJnR3XwucDUwCXiH0ZpltZqPN7IRosguBM83sReBu4Lvu7vkKOldOPDGU486YkXQkIlIu6uuhd284+ODCLzur0n93n0g42Rl/7vLY/TnAobkNrTBqa8Mv6auvwu67Jx2NiJSyZctg0iQ45xzokEDZZkVWisYNGxaula4+6SKysR58MJyXy/fIRM2p+ITesyd8/etN7ekiIu3R0ADXXAPbbw8HHphMDBWf0CH8ms6ZE5pe1B9dRNqqoSFcgOv//i8MRJ/UeTkldKBv33D7u9+pyEhE2m7q1NCzBUJtS1J1LUrowKxZ4VZFRiLSHjU14VwcFLaQKJ0SOk1FRhBGNFKRkYi0xS67hB3CmppwHZek6lqU0Akf/oMPhl/Y2loVGYlI24wfHxL6jTcmmz+U0CODB8Mxx8ATT4QvRkQkW2PHhjqWffZJNg4l9Ji6OnjzTZg5M+lIRKRULFwI06aF/JFqR0+KEnrMiSeGtnQVGYlItsaNC0f1hRyZqDlK6DHduoVml7Fj1ewiIq1raAhjKuy8M+yZPo5bApTQ09TVwdtvww9/qP7oItK8VDHRW2+FnFEM+UIJPc0OO4TbW29VkZGINK9YionilNDTPPtsuFWRkYi0pFiKieKU0NPU1DSN0q0iIxFpzg47wPr1YRDoJIuJ4pTQ01RXw8SJIZkfd1xxfEkiUnzuvTfc3nJL8eQJJfQMBg4MXRiffBLWrUs6GhEpRvX1UFUVergUCyX0ZtTWwnvvwfTpSUciIsXmjTdCAWJSA1k0Rwm9Gd/8JnTtqiIjEfmi1IA4w4YlG0c6JfRmbL45HH98uOjO2rVJRyMixaS+PgwCnRpLoVgoobegri6MPjJqlPqji0hwzz3w4otw0EFJR/JFSugt2GabcDtmjIqMRCTkgNNOC/dvvbX4coISegtSX5aKjEQEQg5INcGuWVN8OUEJvQXxIqOOHVVkJFLpdtwx3HboUDzVoXFK6C2orobJk2HTTUOTS7EUD4hIMl5/PZT7//SnxVMdGtcp6QCK3RFHwIgRcP/94UI8m26adEQikgT30LulpgauuSbpaDLTHnoW6urgo4/C3rqIVKZZs+C114pjIIvmZJXQzWyQmb1mZnPN7OJmpqk1szlmNtvM/p7bMJM1cCBsu62KjEQq2dix4VzayScnHUnzWm1yMbOOwM3A0cAC4Fkzm+Duc2LT9AcuAQ5192Vm9qV8BZyEzp3hpJNC/9PPPoPNNks6IhEppFRzy5FHQo8eSUfTvGz20AcAc919nruvBu4BhqZNcyZws7svA3D3RbkNM3l1dbB8OZx5ZvH1PRWR/BozJly/paoq6Uhalk1C7wXMjz1eED0Xtxuwm5k9aWYzzGxQphmZ2SgzazSzxsWLF7cv4oSkTobedZeKjEQqSUNDqBYH+M1vinvbz9VJ0U5Af6AGGAHcbmbbpE/k7re5e5W7V/Uo5uOWDJ54oml0EhUZiVSOxx8v7mKiuGwS+jtAn9jj3tFzcQuACe6+xt3fBF4nJPiyoSIjkcr0peiMYLEWE8Vlk9CfBfqb2U5mtgkwHJiQNs0DhL1zzKw7oQlmXg7jTFx1NTz6KGyxBRxySPEVFIhIfsyZE0Ywu+yy4iwmims1obv7WuBsYBLwCjDW3Web2WgzOyGabBLwgZnNAR4HfuzuH+Qr6KQcdhicfjrMmAGffJJ0NCKSb+vXh6HmBg+G0aOLO5lDlm3o7j7R3Xdz913c/eroucvdfUJ03939Anff0933cfd78hl0kmprYeVK+Mc/ko5ERPKtoQEWLCi+kYmao0rRNjrkEOjVS0VGIpWgvj70cDv++KQjyY4Seht16BD20h95JFwOQETK07p1MG4cDBkCW22VdDTZUUJvh9ra0HXxe98r7j6pItJ+t94KCxfCfvslHUn2lNDbYf360Cd93DgVGYmUo4YGOO+8cP/aa0tnG1dCb4dp05ruq8hIpPxMmdJUTFRK27gSejvEi4w6dCjuQgMRabtUm3kpFBPFKaG3Q3V1+AXv1i20rxV731QRaZuXXoIuXeCKK4q/mChOCb2dDjkEfvADmDkTlixJOhoRyZU1a2D8ePjWt+Dyy0snmYMS+kapqwtdm+67L+lIRCRXHnsMli4tnWKiOCX0jfDVr8Juu6nISKScjB0b2tAHZbwIeHFTQt8IZuFXfOpUeP/9pKMRkY21enUYEP7EE0tzQHgl9I1UVxf6pY8aVTp9VUUks9/9Dj78EPbZJ+lI2kcJfSN9/HHYU58wQUVGIqWsoQEuvjjcv/zy0tyWldA3UrzgoJQKEERkQ48+Gjo5QOluy0roG6mmpqmtTUVGIqWrc+dwW2rFRHFK6BspVWS0/faw++6l1WdVRJq8+GLo3XLVVaVVTBSnhJ4D1dVwzjmhumz+/KSjEZG2WrEiDFozfDj87GelmcxBCT1nUkUIY8cmG4eItN3EifDpp6VZTBSnhJ4ju+wCBxyghC5Siurr4UtfgsMPTzqSjaOEnkO1tfDMM/Dmm0lHIiLZWr4cHn4YTjkFOnVKOpqNo4SeQ7W14fass0qzD6tIJbrhBvjsM9hrr6Qj2XhK6Dm0cGHo8jRpkoqMREpBQ0Po1QJw0UWlv80qoefQ1KngHu6XamGCSCV55JFw6Q4oj21WCT2HVGQkUppKuZgoTgk9h1JFRn37Qp8+pduXVaRSPPcc9OgBv/hF6RYTxSmh51h1NVx4IcybB6+8knQ0ItKcZcvC+a7TToNLLy39ZA5K6Hlx8snhCozqky5SvB58MAw3V+rFRHFZJXQzG2Rmr5nZXDO7uIXpTjYzN7Oq3IVYenr2DAUK9fVNJ0lFpLjU10O/fnDggUlHkjutJnQz6wjcDAwG9gRGmNmeGabbEjgPeDrXQZaiurrQ5DJ7dtKRiEi6Dz4Il8utrQ1H0+Uimz30AcBcd5/n7quBe4ChGab7BXAdsDKH8ZWsVLPLueeWft9WkXLzq1/B2rWwxx5JR5Jb2ST0XkD8GoILouc+Z2b7A33c/eGWZmRmo8ys0cwaFy9e3OZgS8kbb4SE/vjjKjISKSYNDSGhA/zwh+W1bW70SVEz6wDcAFzY2rTufpu7V7l7VY8ePTZ20UVNRUYixemhh8qrmCgum4T+DtAn9rh39FzKlsDewFQz+zdwMDCh0k+MxouMzEq/YEGkXKxeHW47diyPYqK4bBL6s0B/M9vJzDYBhgMTUi+6+0fu3t3d+7l7P2AGcIK7N+Yl4hKRKjLq3z8ULhx8cNIRiQhAY2Mo/hs9ujyKieJaTejuvhY4G5gEvAKMdffZZjbazE7Id4ClrLoaLrkkXLSrsaJ/3kSKw8KFMG0afPe75VNMFJdVG7q7T3T33dx9F3e/OnrucnefkGHamkrfO4878cQw+Gx9fdKRiMi4ceHcVupS1+VGlaJ51q0bHHtsqBpNnYgRkWTU18Pee8OeX6ikKQ9K6AVQWxsGj54xI+lIRCrXggXw5JPlVeqfTgm9AIYODc0uP/5xefV5FSklqb7n/fsnG0c+KaEXwOzZsG4dPPWUioxEktDQADfdFO6ffnr5boNK6AUQL1xYtaq8ChlESsH995dvMVGcEnoBaCQjkWStWBFuy7GYKE4JvQCqq0MBw957w5ZbwoABSUckUllmzAgX4iqXkYmao4ReINXVcMUVYZSUadOSjkakcsydCzNnwsiRodCvXJM5KKEX1JAh0LWrioxECunee8PtsGHJxlEISugFtPnmcPzxMH58uBaziORffX3YK99xx6QjyT8l9AKrqwujpUyZknQkIuXvtdfgxRfLu5goTgm9wAYNCnvql15avn1hRYrFr38dbnfaKdk4CkUJvcCefz70RZ85U0VGIvnU0AB33BHuDx9eGduaEnqBxUcyUpGRSP7U11feqGFK6AWmkYxECuPjj8NtuRcTxSmhF1iqyGj//UNi32+/pCMSKT/u4dpJBxxQ/sVEcUroCaiuhquvDuXIkycnHY1I+Zk1K/RwGTWq/IuJ4pTQEzJwIGy3nYqMRPKhvj40tZx0UtKRFJYSekI6dw4r24QJ8NlnSUcjUj7cwwhhAwdC9+5JR1NYSugJqq2F5cth4sSkIxEpH889B2+8UTnFRHFK6AmqqYEePcLehIjkxg03hMtU9+qVdCSFp4SeoE6d4JRT4MEH4aqrKqPwQSSfnnoK7r47DGbxrW9V3jalhJ6wvfYKBUajR6tyVGRj/e1vlVdMFKeEnrAPPwy369dX5gookktLl4bbSiomiuuUdACV7sgjQ9PL2rWVuQKK5Mr69fDEE3DYYWHsgZqayul/nqI99IRVVzeNRv6Tn1TeCiiSK8zcCckAAAwFSURBVA0N8M47cNZZlVVMFKeEXgTOPDOckX/uuaQjESld9fXQpUsYRKZSZZXQzWyQmb1mZnPN7OIMr19gZnPMbJaZPWZmfXMfavnq0CH0SZ80qalNXUSyt25dGGpuyJAwEHulajWhm1lH4GZgMLAnMMLM9kyb7Hmgyt2/CowD/jvXgZa7urpwUvTBB5OORKT0PPEEvPdeZRYTxWWzhz4AmOvu89x9NXAPMDQ+gbs/7u4rooczgN65DbP8DRgAffvq2i4i7XHjjeFyGpVW6p8um4TeC5gfe7wgeq45I4H/zfSCmY0ys0Yza1y8eHH2UVYAs9DsMnky/Pzn6o8ukq1//Qvuvx/WrIHjjqvsbSenJ0XN7FSgCvhVptfd/TZ3r3L3qh49euRy0WVhjz1CW+DVV6vISCRbY8Y03a/0Wo5sEvo7QJ/Y497Rcxsws6OAy4AT3H1VbsKrLAsXhlt3rZgi2VqyJNxWajFRXDaFRc8C/c1sJ0IiHw58Oz6Bme0H3AoMcvdFOY+yQnzjG01FRp07V/aKKZKNNWvCCdFjjgnbSyUWE8W1uofu7muBs4FJwCvAWHefbWajzeyEaLJfAVsA95rZC2Y2IW8Rl7Hq6qZRys89t7JXTJFsPPZYKPf/0Y8qt5gozjx1JZsCq6qq8sbGxkSWXczcQ1t6z54wZUrS0YgUtzPOgPvug/ffbxp8vdyZ2Ux3r8r0mipFi0yqt8u0aaFfrYhktnp16N1y4omVk8xbo4RehOrqwoWGxo1LOhKR4jV5cqisrvRiojgl9CK0117hTyMZiTTvppvCtVu6dk06kuKhhF6kamtDwcQll6g/uki6qVPDtY9WrYJBg7SNpCihF6nddgu3112nIiORdKneYKrZ2JASepF6881wqxVW5ItSHQZUTLQhjVhUpGpqQnHRmjUqMhKJW7EiHLEOHQoHHaRiojgl9CJVXQ133RXa0keO1AorkjJxInz6KZx3XqiuliZqciliw4bBAQfA008nHYlI8aivhy9/GQ4/POlIio8SepGrq4PGRpg3L+lIRJK3fDk8/DCcckpoP5cNKaEXudracKs+6SLw0EPw2WcqJmqOEnqR69s3nPjRSEZS6Roa4Be/CKMSHXpo0tEUJyX0ElBXBy+8ABdeqP7oUpkaGkI9xpw5odxf55UyU0IvATvvHG5/8xsVGUllmjo1VIVCuM6R6jIyU0IvAXPmhFsVGUmlitdhbLqp6jKao4ReAlJFRqAiI6lMu+8eLi192GFhUAvVZWSmhF4Cqqth/Phw/9vf1sosleeBB8IA6jfcoPW/JUroJeL448Oe+VNPhaYXkUoydmw4l1SVcZweSVFCLyF1dfDqq/Dyy0lHIlI4H3wAjz4aajLMko6muCmhl5CTToIOHdQnXSrLfffB2rVNRXbSPCX0EvKlL8GRR4aErmYXqQQNDWFMgN69Yd99k46m+Cmhl5i6Opg7F849V/3Rpbw1NIQdmDfeCNc/nzEj6YiKnxJ6iendO9zefLOKjKS8xYuJ3FV/kQ0l9BLz/PPhVkVGUu6OOKLpvkYlyo4SeompqQmVchBOkGoll3K1cmXYcRk2TMVE2dKIRSWmuhqmTAnXg+7SBQ4+OOmIRPLjlltgu+3gzjvDui6t0x56CTrkELjmmjCQ9JQpSUcjknvvvhuqQ884Q8m8LbJK6GY2yMxeM7O5ZnZxhtc3NbP66PWnzaxfrgOVDQ0fDj16wI03Jh2JSO7dfnu4quL3v590JKWl1SYXM+sI3AwcDSwAnjWzCe4+JzbZSGCZu+9qZsOB6wCNKZJHXbrAWWfB6NFw0UVw8snh+alTm9rVc3G/ujr0pMn1fCsx1mKPr1hiXbsWrr8eBgyAXXZB2iCbNvQBwFx3nwdgZvcAQ4F4Qh8KXBndHwfcZGbmrvKXfDrwwHB7/fXhz6yp4CgX983CYLzvv5/b+ebjfrHHWuzxFWOszz8fflR0MjR72TS59ALmxx4viJ7LOI27rwU+ArZLn5GZjTKzRjNrXLx4cfsils+99NKG17aI/3zm4r5701++llEpsRZ7fMUY67p16pbbVgU9Kerut7l7lbtX9ejRo5CLLks1NaHppWPH0E93001ze3+zzUKTzmab5W8ZlRJrscdXrLGqW27bZNPk8g7QJ/a4d/RcpmkWmFknYGvgg5xEKM2qrg79c/PdlrrPPsm36ZZDrMUeXzHGquaWtrHWmrmjBP06MJCQuJ8Fvu3us2PT/AjYx91/EJ0UPcndW7w2WlVVlTc2Nm5s/CIiFcXMZrp7xivDt7qH7u5rzexsYBLQEfiTu882s9FAo7tPAO4A/mpmc4GlwPDchS8iItnIqlLU3ScCE9Oeuzx2fyUwLLehiYhIWxT0pKiIiOSPErqISJlQQhcRKRNK6CIiZaLVbot5W7DZYuCtdv57d2BJDsMpBXrPlUHvuTJszHvu6+4ZKzMTS+gbw8wam+uHWa70niuD3nNlyNd7VpOLiEiZUEIXESkTpZrQb0s6gAToPVcGvefKkJf3XJJt6CIi8kWluocuIiJplNBFRMpESSV0M+tiZs+Y2YtmNtvMrko6pkIws45m9ryZPZR0LIViZv82s5fM7AUzK/vrLJvZNmY2zsxeNbNXzKysrwRuZl+JvtvU38dmdn7SceWbmf1XlLteNrO7zaxLTudfSm3oZmZAV3dfbmadgSeA89x9RsKh5ZWZXQBUAVu5+3FJx1MIZvZvoMrdK6LgxMz+AvzL3f9oZpsAm7v7h0nHVQjRQPTvAAe5e3uLDYuemfUi5Kw93f0zMxsLTHT3MblaRkntoXuwPHrYOfornV+kdjCz3sA3gT8mHYvkh5ltDRxOGFcAd19dKck8MhB4o5yTeUwnYLNo4KDNgXdzOfOSSujwefPDC8Ai4J/u/nTSMeXZb4GfAOuTDqTAHJhsZjPNbFTSweTZTsBi4M9R09ofzaxr0kEV0HDg7qSDyDd3fwf4NfA2sBD4yN0n53IZJZfQ3X2du+9LGNt0gJntnXRM+WJmxwGL3H1m0rEk4DB33x8YDPzIzA5POqA86gTsD/zB3fcDPgUuTjakwoial04A7k06lnwzs27AUMIPeE+gq5mdmstllFxCT4kOSR8HBiUdSx4dCpwQtSffAxxpZn9LNqTCiPZmcPdFwP3AgGQjyqsFwILY0eY4QoKvBIOB59z9/aQDKYCjgDfdfbG7rwHuAw7J5QJKKqGbWQ8z2ya6vxlwNPBqslHlj7tf4u693b0f4bB0irvn9Be9GJlZVzPbMnUfOAZ4Odmo8sfd3wPmm9lXoqcGAnMSDKmQRlABzS2Rt4GDzWzzqIPHQOCVXC4gqzFFi8gOwF+is+IdgLHuXjFd+SrIl4H7wzpPJ+Dv7v5IsiHl3TnAXVETxDzg9ITjybvox/po4PtJx1II7v60mY0DngPWAs+T40sAlFS3RRERaV5JNbmIiEjzlNBFRMqEErqISJlQQhcRKRNK6CIiZUIJXcqWmV0WXdluVnRFv4Oisvo9k45NJB/UbVHKUnT52RuAGndfZWbdgU3cPacXQxIpJtpDl3K1A7DE3VcBuPsSd3/XzKaaWRWAmY00s9eja+zfbmY3Rc+PMbM/mNkMM5tnZjVm9qfoOuVjUguIpmmspGvzS3FTQpdyNRnoEyXsW8zsiPiLZtYT+DlwMOGaObun/X83oBr4L2AC8BtgL2AfM9s3muYyd68CvgocYWZfzdu7EcmCErqUpei6+QcAowiXpq03s+/GJhkATHP3pdGFktKv9vcPD+2RLwHvu/tL7r4emA30i6apNbPnCCXcewFqm5dEldq1XESy5u7rgKnAVDN7CfjPNvz7quh2fex+6nEnM9sJuAg40N2XRU0xOR1OTKSttIcuZSkas7J/7Kl9gfiIOM8Smkm6RaPHnNzGRWxFuG75R2b2ZcJlYEUSpT10KVdbAL+PLre8FphLaH4ZB+F662Z2DfAMsJRwGeaPsp25u79oZs9H/zcfeDK34Yu0nbotSsUysy2iAcc7EQbR+JO73590XCLtpSYXqWRXRuPTvgy8CTyQcDwiG0V76CIiZUJ76CIiZUIJXUSkTCihi4iUCSV0EZEyoYQuIlIm/j+aXEmvzpGDGgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully.\n", " Current function value: 371.740238\n", " Iterations: 427\n", " Function evaluations: 721\n", "\n", "Fit result: (array([ 0.75387138, 10.12916274, 7.09539653, 1.00255828]), 371.74023833526115, 427, 721, 0)\n", "\n", "------------------- Parameter summary --------------------\n", " A = 0.753871, free: T, restricted: T, related: F\n", " mu = 10.1292, free: T, restricted: F, related: F\n", " sig = 7.0954, free: T, restricted: T, related: F\n", " off = 1.00256, free: T, restricted: F, related: F\n", " lin = 0, free: F, restricted: F, related: F\n", "----------------------------------------------------------\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU5dXA8d/JAgpRwiaibCq4gLvIYgABUQFxaauorXWjWlzaUmvr0qrYilalRa2+WvelddeK66u8WnADNCiiuAURZFNAFiWyJTnvH88dMklmS3Jn7syd8/185pPMfe7MnJuZnHnus11RVYwxxoRXQdABGGOMSS9L9MYYE3KW6I0xJuQs0RtjTMhZojfGmJCzRG+MMSFnid4kJSKLRGREnLLtReQ5EVkvIk+IyM9E5JVMx5guIjJdRH4Rp2yiiPwr0zH5QUTuEJErmvjYDSKyu98xmfQpCjoA4z8RWQR0AqqBrcDbwHhVXZKGlzvRe632qlrlbft3VCwK9FLVBWl4bVOPiJwJ/EJVByXaT1XHp/h804F/qerdUY8taU6MJvOsRh9ex3r/kJ2Bb4B/pOl1ugOfRyV500wiktYKmIgUpvP5TfaxRB9yqroJeBLoHdkmIseIyPsi8p2ILBGRidGPEZGfi8hiEflWRP4Y77lF5GrgSuBk73R+nIicKSJveuWve7t+4JWfHOd5zhGRT0TkexH5WEQO9rZfKiJfRG3/UdRj6jSbiEgPEdFIkvTiWOg99ksR+VnUvmd7r7dWRF4Wke5RZUeKyKdeU9StgCT5E28nIo95r/OeiBzgPc/vReSpesd5i4jcHOdvsEhELhGReUCliBSJyAAReVtE1onIByIyNGr/BscnIvsAdwADvb/3Om/f+0XkdhF5UUQqgWHetmuinu94EZnrfSa+EJGRIjIJGAzc6j3frd6+KiI9vd/biMiDIrLK+8z8SUQKomJ8U0Qme3/rL0VkVKJjSPK3Nk2lqnYL2Q1YBIzwfm8FPAA8GFU+FNgP90W/P67Gf4JX1hvYAAwBWgJ/B6oizxfjtSbiTu0j988E3oy6r0DPBLGeBCwDDsUl1Z5A96iyXbw4TwYqgc5xXreH91pFQGvgO2Avr6wz0Mf7/XhgAbCPt++fgLe9sg7A97jmqGLgt96x/yLBsW+N2v9i4Evv985evKXevkXASuCQBO/ZXKArsD2wK/AtMNo7/iO9+x2THF+dv7+37X5gPVDmPdd23rZrvPJ+XvmRXvmuwN5e2fT6xx/9ngIPAlOBHbz34HNgXFQsW4FzgELgPGC59z7HPQa7+X+zGn14PePV6CL/wDdGClR1uqp+qKo1qjoPeAQ43Cs+EXheVV9X1c3AFUBNGuP8BXCDqr6rzgJVXezF+YSqLvfifAyowCWlVNQA+4rI9qq6QlXne9vHA9ep6ifqmpuuBQ70avWjgfmq+qSqbgVuAr5O8jpzovb/Oy6JDlDVFcDruC8rgJHAalWdk+C5blHVJaq6ETgNeFFVX/SOfxpQ7sWY6Pjimaqqb3nPtale2TjgXlWd5pUvU9VPkzxfpAnoFOAyVf1eVRcBfwN+HrXbYlW9S1WrcRWOzrg+naYcg2kiS/ThdYKqluISz4XADBHZGUBE+ovIf73T7fW45NfBe9wuwLZOW1WtxNUk06Ur8EWsAhE53WtOWOd9ae0bFWdcXswn445rhYi8ICJ7e8XdgZujnnMNroa5Kw2PXaPvxxG9fw2w1HsecIntNO/304CHUn0uL86TInF6sQ7CndEkOr5Unru+uO9BEh1wZy+Lo7Ytxv0tI7Z9UarqD96vJU08BtNEluhDTlWrVfVp3AicyEiMh4Fnga6q2gbXrhtpi16B+8cHQERaAe3TGOISYI/6G70a9l24L6n23pfWR1FxVuKapSJ2jn68qr6sqkfiapCfes8Veb1fqmpp1G17VX2bhscu0ffjiN6/AOiCa54AeAbYX0T2BcYQNRopjuilZJcAD9WLs7Wq/jXJ8cVbjjbRMrUx34MUHrca1zTTPWpbN1xTXFIJjsH4zBJ9yIlzPNAW+MTbvAOwRlU3iUg/4KdRD3kSGCMig0SkBfBnmvc5+QZINOb6buBiETnEi7Wnl+Rb45LMKu84zsLV6CPmAkNEpJuItAEuizrmTl7nYmtgM67PIdL8dAdwmYj08fZtIyKR5pUXgD4i8mNxnbq/pt4XSAyHRO0/wXu9WVCnI/xh4B1V/SrJc0X7F3CsiBwtIoUisp2IDBWRLkmO7xugi/fepeoe4CwROUJECkRk16jaddz3z2uOeRyYJCI7eO/bRV7sCSU5BuMzS/Th9ZyIbMB1eE0CzohqAz0f+LOIfI8bNfN45EHePhfgktMKYC2uOaKpJgIPeM0PY+sXquoTXnwP4zpCnwHaqerHuPbembhksx/wVtTjpgGPAfOAOcDzUU9bgEs4y3FNM4fjOgJR1f8A1wOPish3uLOEUV7Zalyb+l9xzVW9ol8zjqm4Joi1uLbpH3vt9REPeLEna7apQ92ch+OBy3FfdkuA33vHFvf4gNeA+cDXIrI6xdd6BzgLmILr05lBbS39ZuBEb9TMLTEe/ivc2dVC4E3c+3hvCi+b6BiMz8Q1Qxpj0kFEuuGaJXZW1e+CjsfkJ6vRG5MmXpv9RcCjluRNkGwJBGPSwGt7/gY3CmVkwOGYPGdNN8YYE3LWdGOMMSGXdU03HTp00B49egQdhjHG5JQ5c+asVtWOscqyLtH36NGD8vLyoMMwxpicIiKL45VZ040xxoScJXpjjAk5S/TGGBNyluiNMSbkLNEbY0zIWaI3xpiQs0RvjDEhZ4nemFwxaxbcdx+sXx90JCbHWKI3Jttt3AgXXQSHHQZnnw277ALXXBN0VCaHWKI3Jttdfz1MmQLnnQdvvgk//Sl07hx0VCaHZN0SCMaYei6+GIYMgeHD3f2yMvdTFZ59Fo45BorsX9nEZzV6Y7LVM8/AokVQUlKb5KO98QaccIKr8RuTgCV6Y7LR8uVw2mlwySXx9xkyBE4+Ga6+GubOzVxsJudYojcmG112GWzdCtddl3i/226Ddu3gjDOgujozsZmcY4nemGwzezY8+CA1Ey5i8tO706ED/O1vcfJ4+/bw97/DvHnw8ssZD9XkBkv0xmSbyy+nquPODH7pciZOhG+/hauugkMPhYqKGPufeCKMGAEimY7U5Ajrqjcmm1RXw557cuXs45k1fwdqatzmykr44AM34GblynqPadECpk3LeKgmd1iN3phsUlgIt9/OzEN/vS3JR9TUwL77Jnjs99/DO++kNTyTmyzRG5Mtamrg3XdBlXHj3KjKaCUlbmJsXGedBccd5zpxjYliid6YbPHGG9CvHzz9NMce23AOVFERHHtsgsefdRZ88w1MnZrWME3usURvTLZ44AHYYQcYOZI2bWDtWjf5NXJbuxbatEnw+JEjoXt3uP32jIVscoMlemOyQWUlPPEEnHQStG7dtOcoLIRzz4XXXnMzao3xWKI3Jhs8/TRs2ABnntm85xk71v3873+bHZIJDxteaUw2ePRR2H13GDSoec/TsycsXAi77eZPXCYUrEZvTDZ46il4+WWqa4TJk0k8GzYZS/KmnqSJXkTuFZGVIvJRnHIRkVtEZIGIzBORg+uV7ygiS0XkVr+CNiZ0ttuOCu1J376kNhs2kTVr4JRTbPSN2SaVGv39wMgE5aOAXt7tXKB+l/9fgNebEpwxeWHSJLjhBsrK3JI1lZVuc/Rs2Hiqq2l4BtCmDbz6Kjz+eEbCN9kvaaJX1deBNQl2OR54UJ1ZQKmIdAYQkUOATsArfgRrTOiowq23wpw59OlDSrNhI8m9bVvo1i3GGcDCQhgzBl580SZPGcCfNvpdgSVR95cCu4pIAfA34OJkTyAi54pIuYiUr1q1yoeQjMkRc+fC11/D6NEpzYatqGBb8866dW7Z+phnAMcd53Z4661MHYnJYunsjD0feFFVlybbUVXvVNW+qtq3Y8eOaQzJmCzz0kvu58iRKc2Grd+8U1/kDKB6+JFUFbXk9lHPNr1T14SGH4l+GdA16n4Xb9tA4EIRWQRMBk4Xkb/68HrGhMeLL8Ihh0CnTinNho3VvBOtpARGj4a+Q0u4l3FUbOrS9E5dExp+jKN/FpfQHwX6A+tVdQXws8gOInIm0FdVL/Xh9YwJh5oa2Hln1xaTonHjoLzcza2KpajIXUJ2zRr4Zc1tbmOiJY5NXkhleOUjwExgL2+Y5DgRGS8i471dXgQWAguAu3BNNsaYZAoK4Mkn4dLU6z+xmndKS11zfOQMYN99a2v9RWylIyuTL3FsQi2VUTenqmpnVS1W1S6qeo+q3qGqd3jlqqoXqOoeqrqfqpbHeI77VfXCdByAMTnru++AOEMk64nss8ce8Kc/QVVV/Oad6E7duRzIbVyQfIljE2o2M9aYIKjCXnux5pw/JJ0kFT3SJtFEqsiXwa9/7b4IAMrpy1CmU1xYk3iJYxNqluiNCUD1x5/B119zyd178sEHiSdJpTKRKvrLYO1at5DlQQfB0X8dRkdWs+aN+dtq/amcQZhwsURvTIZVVMB1I2cAMJ3DUa1bXr89PZWJVPG+DI65YajbMH36ttf2ZZkFk1Ms0RuTYWVlsPvS11nBziygZ4Py+u3pqUykivdl0OaAHtCjx7Zli5uyzILJfZbojcmwPr2VIcxgBocD0qC8/iSpVCZSJfwymDIFLnYT1FNdZsGEiyV6YzLsF2fXcF3Lq7mLc7ZtKymBhx6KPYomlYlUCb8MTjgBDjsMSO3swISPaP0GwoD17dtXy8sbjNA0JjTWr3etKevW1W4rLXVX/0t4TdimUoVXXoHWrVm/36DMvrbJGBGZo6oxZ9/ZFaaMybA2c2ewtryLGxSfCSJw4YXQuzdtpg5i7drMvKzJHtZ0Y0ymnXUW/P73mX3NoUNhxgwbS5mnLNEbk0lLlsCXX8Lhh2f2dYcNc21Gc+dm9nVNVrBEb0wmve5dbG3IkMy+buT13n47s69rsoIlemMypLoaPvjHDNZLG/7+f/tnthWlSxfo2hXefTeDL2qyhSV6YzIgMiN1+3de5w0dxJVXF2Z+Ruqbb8J992XwBU22sFE3xmRAWZlbcmCwzqAN6+vMSM3YGvHdumXohUy2sRq9MRkQmZG6kk5UsCcQwIzUdevgggtg2rQMvqjJBpbojcmAcePgrJYP82tu3rYt4zNSW7d2TTcvvJDBFzXZwJpujMmAY4+FHlV30pJKbuE3QMP1atKuuBj69bORN3nIavTGZECb1lUMavkuh/5qYNz1ajLisMPg/ffhhx8y/MImSJbojcmEefNcch04MNg4DjvMXX7K1pPKK5bojcmEmTPdz6AT/cCBbjz9mjXBxmEyytrojcmEb75xwxu7dw82jvbt4auvgo3BZJzV6I3JhD//Gb74wq0kmS2ybIlykz6W6I1JozoX4r65KDsWj3ztNXd28fnnQUdiMsQSvTFpEln24M0//S9PfzuEu69YnB0X4t55Z7eKZqTfwISeJXpj0iRyIe7+m2cwkJks2rhTdlyIe++9YccdYfbsgAMxmWKJ3pg0iSx7MJCZvM9BbGL77LgQd0EBHHqoJfo8YonemDQZN85NlDqUd5mJG1aZNRfiHjCgdmy/CT0bXmmMz6qrYcoUuPZa2GvrR7TmB2YxAAhg2YN4jj7aTc2trIRWrYKOxqSZJXpjfFRRAWPHup+VlVC8XTVv7Diaa58ZyCPDgo4uyuDB7mbygjXdGOOjSAdsZaW7/9amQxi64QX6n9wj0LjqDPP8m3eN8KoqmzyVJyzRG+OjSAdsRGs2BN4BGxnmOXGiu/jJVVe5vtjvfnKW1erzRNJELyL3ishKEfkoTrmIyC0iskBE5onIwd727iLynojMFZH5IjLe7+CNyTbjxrkOV4Ad+I51lPK7Fv8ItAO2/llG5OpWk189yNXov/46uOBMRqRSo78fGJmgfBTQy7udC9zubV8BDFTVA4H+wKUiskvTQzUm+x17rOtwBehLOUVU82XxnoF2wNY/ywB3f+2e/d0dG2YZekkTvaq+DiRa6u544EF1ZgGlItJZVbeo6mZvn5apvJYxua5NGzeYRRVeu9Yl0Ke+OjTz685HiT7LiCgpgcMuPNh9K3mJPmY7vgkFP5LvrsCSqPtLvW2ISFcRmeeVX6+qy2M9gYicKyLlIlK+atUqH0IyJgu88w706gXt2gUaRvRZRkRREYz+yfaw//4we3bcdvzAl2swvkhrLVtVl6jq/kBP4AwR6RRnvztVta+q9u3YsWM6QzImM1RdTbl//6AjqXOWUf/qVtVX/ZnHe17OXnu5dvv67fiBL9dgfOHHOPplQNeo+128bduo6nKvM3cw8KQPr2lMdqupgSuucDX6LFVRAWOvOoaKitgrFgc9Wsj4x48a/bPA6d7omwHAelVdISJdRGR7ABFpCwwCPvPh9YzJfoWFcN55MGJE0JHEVVYGH31QzYDK/+MA5jYoz5rlGkyzpTK88hFgJrCXiCwVkXEiMj5quOSLwEJgAXAXcL63fR9gtoh8AMwAJqvqh74fgTHZaO5cd6GRLNanD1Sr8CQnMp47GpRnzXINptmSNt2o6qlJyhW4IMb2acD+TQ/NmBz2m9/Apk1ZPXRx3DgoLy/gnQ396I+Ls6QEbr8dTj3Vrdezxx5w2WUwYYI7STG5yYY8GuO3qiooL8+KjthEIqNx3qEf+/Eh2/MDRUXQu7eNwAkbW9TMGL99/LFb/rdfv6AjSSgyGofn+sNx1fzw5vtQVsZOO7kEH5lkFT0CZ+XKQEM2TWQ1emP89s477meW1+i3iXwhec1M8WbS2gic3GU1emP8Nnu2myTVs2fQkaSmUyfXedy7NxBpu4cNG2p3sRE4uc1q9Mb4beJEeOYZEAk6ktQdcAAUFwPxZ9LaCJzcZYneGL/tumvuLf/7ySduaM2qVQln0prcZIneGB9EFgQrK53PayfcTPXqtUGH1DirV8PNN2f1cFDTdJbojUlRvNUdoxcEK1v/AsOnTuDIYVW5NRzx4IPdQPlIR7IJFUv0xqQg0eqO0Rf26M9svmB3ZnzcMbcWBGvd2g2rsRp9KFmiNyYF8a7SVFZWdzhif2Yzm/65ORyxXz9Xo4+1wpnJaZbojUlBorHlkQt77MIyurCM2fTPzeGI/fvDdtvZpQVDyBK9MSmId5Wms8+uHY7Ym4+ppoB36JebwxHPPBOWL4fOnYOOxPhMNMtO0/r27avl5eVBh2FMHevXQ48esG5d7bbSUli0qN6ww8pKaNmy4UB0Y9JMROaoat9YZVajNyYFKY8tb906t5P8ddfBT34SdBTGZ5bojfFDdTWMGQPPPRd0JM3z/ffw7LOwcWPQkRgfWaI3xg+ffgovvOAtB5nDBgxwyyy/917QkRgfWaI3pgnqT56qeWumKxgwINjAmisS/8yZwcZhfJXDjYnGBKOiAsaOdT8rK93kqT1azuTY0vYUZvHFwFOy006w++4wa1bQkRgfWaI3ppHKyhpemKNX5SxebTGAo3Jpxcp4Tjkl6AiMzyzRG9NIffrA9Om19wupYildWNBtBEcFFpWPJk0KOgLjM2ujN6aR6k+eqqaIE0teps1VE4ILym81NdvWe4juj7jxRrjhhoYLu5nsZhOmjGmk+pOnhBralBY0nDyVq2pqYJdd4LTTqPjl5Dr9ESLuVlPjpgzsuSc89hjketdEGNiEKWOaKNbSxPUnT9WMGsPao04OR5IHKCiAPfaAWbMaLOamGvui4Sa7WaI3Jo5ESxNvU1MDb78NbdsGFWZ6DBwIc+ZwwD5bGizmFi0nV+nMQ5bojYkj0dLE23zyiWvLGTgwkBjTZsAA2LSJ3x4xr8FibtFycpXOPGSJ3pg4Ei1NvE1kvHkYEz0wbLuZCZfuyclVOvOQJXpj4ki0NPE2M2dCu3bh643s0gWuuYbtR5Q1WMzNLhqee2wcvTH1VFfDlClw7bVu2ZdoDWqwRxzhhp6EYaJUfX/8Y9ARGJ9Yjd6YKNEdsGvXuutlH3QQfP55nBrsqafCH/4QVLjptWkTvP6664mOId7F0k32sURvTJSUOmAjVqyAJUsyGl9GffIJHH44vPxyg6KURiSZrGGJ3pgoKXXARtx2G+y2G2zYkJHYMm6//WCHHeDNNxsUNeoL0QQuaaIXkXtFZKWIfBSnXETkFhFZICLzRORgb/uBIjJTROZ720/2O3hj/JZSB2zEG2+4dp1E4w9zWVGRG0301lsNihr1hWgCl0qN/n5gZILyUUAv73YucLu3/QfgdFXt4z3+JhEpbXqoxqRf5ELf0WIOIdy8GWbPhsGDMxZbIAYNgg8/rHuxXBr5hWgClzTRq+rrwJoEuxwPPKjOLKBURDqr6ueqWuE9x3JgJdDRj6CNSZeUrw1bXu6SfdgTfVmZ+yPUuxBJyl+IJiv4MbxyVyC6R2qpt21FZIOI9ANaAF/EegIRORd3NkC3bt18CMmYNHvjDfdz0KBg40i3gQPdEg8HH1xnc+QL0eSGtI+jF5HOwEPAGaoac9UMVb0TuBPc6pXpjsmYZjv9dNh7b+gY8pPU7bcP36zfPOTHqJtlQNeo+128bYjIjsALwB+9Zh1jwmGXXeCEE4KOIjPmzYPf/x62bAk6EtNEfiT6Z4HTvdE3A4D1qrpCRFoA/8G13z/pw+sYkx0WLoRbb407kSh0KirczKj33w86EtNEqQyvfASYCewlIktFZJyIjBeR8d4uLwILgQXAXcD53vaxwBDgTBGZ690O9P8QjMmwF1+EX/0Kvv8+6EgyIzI4PsZ4+gibJZvd7ApTxjTWySe7DsqvvgrnGjex9OoFvXvD1KkNiioqqHMVKrvyVDDsClPG+EXVjbgZMiR/kjzA0KEwY0bMqrrNks1+luiNaYyFC90aN2EfP1/f8OFQXAyLFzcoslmy2c8SvTE0oo35ww9dTT7fEv1JJ8E338Duuzcoslmy2c8Svcl7jVqJ8YQTYPVq116dT4qK3EXDY0g0S9Y6abODJXqT9xrdxtyuXX61z0c88gjsv3+D8fTxlo1YudKWMs4WluhN3ku5jfmTT+Coo1zzTT5q2dId+7vvprS7ddJmD0v0Ju+l3Mb8yiswbRrsuGPGYssqhx/uzmReey2l3a2TNntYojd5L+WVGF99FXr2hO7dMxZbVmnfHg44AP7735R2t07a7GGJ3uS9lJYm3roVpk+HESOCCjM7DB/uJott3Jh0V1vKOHukffVKY0Lh3XfdkgdHHBF0JME69lh3EZLvv3crWyZgSxlnD0v0xqSiutq1UQ8bFnQkwRo61N1MTrGmG2NSMXiwa7pp3z7oSIKnamMkc4wlepO3Up7Ms2VL/qxUmYrJk2GvvWDVqqAjMSmyRG/yUqNmw772mpskNXt2psPMTsOGuVr9//5v0JGYFFmiN3mpUZN5XnrJDRfZb7+Mxpi1Dj4YOnVy6/KbnGCJ3uSllCfzqMJzz7nRNq1aZSy+rFZQAKNGwcsvQ1VV0NGYFFiiN3kp5ck8H38MX34JY8ZkLLacMHq0GztpzVk5wRK9yUspT+Z57jn30xJ9XUcdBf/5j5spa7KejaM3eSnlyTwnnQQ77QRduqQ9ppzSpo1bstnkBKvRG5PIHnvY4izxLF8Okya5n1FsDfrsY4nemHjeeAMef9w6HONZswb+9Kc6o28aNWzVZIwlemPiuflm+O1vobAw6EiyU58+0LUrvPDCtk22Bn12skRvTCybN7vhg2PG5OfVpFIhAscf7yZObdgA2Br02coSvTGxzJjhkpetqZvY2LGwadO20Um2Bn12skRvTCz/+Y9bhnf48KAjyW5lZbD77rBoEWBr0GcrS/Qm9Jo0CmTePNcsYbNhEysogE8/hcsuAxJfxMVG4wRHVDXoGOro27evlpeXBx2GCYmKCte6UFHhOgZbtYLSUvjhBzdgZMKEOH2tqu4B9dshTHxbtkCLFjGL6r8PrVvDnnvCY49Br14ZjjOkRGSOqvaNWWaJ3oTZTju5YX71OwghQbKpqXE1VZO6k05yiX7q1JjFsd6HggK3vP/KlRmKMeQSJXr7NJtQizUKJCLm0L8NG6BbN3jooYzEFxrdurlVPteti1lso3GCZYnehFqsUSDRGiSbZ5+FZcugR490hxYuY8e6C6jHqdHbaJxgWaI3oRZrFEi0Bsnm4YfdJCCb4dM4/fpB9+6uHSwGG40TrKSJXkTuFZGVIvJRnHIRkVtEZIGIzBORg6PK/ldE1onI834GbUyqokeBrFvnOmKj1Uk2q1e7SVKnnmpt9I0lAj/9qfv7LV3aoDjRaByTfqmsXnk/cCvwYJzyUUAv79YfuN37CXAj0Ar4ZbOiNMYHSVesfPJJt67NT3+asZhC5dxz3Sqf9b9NTeCSVltU9XVgTYJdjgceVGcWUCoinb3HvgrYVZVNbhgyBK68EvbfP+hIclOPHnD++TYkNQv5cX66K7Ak6v5Sb1vKRORcESkXkfJVdmV5E0NGJtv07g1XX21r2zTH5s3wz3+6JSRSZBOp0i8rGiJV9U5V7auqfTt27Bh0OCbLZGTp29tug3ff9fEJ81RBgfuyvOGGlHa3ZY0zw49EvwzoGnW/i7fNGF+kfenb5cvdFNlHHvHpCfNYcbEbS/nSS7B4cdLdbVnjzPAj0T8LnO6NvhkArFfVFT48rzFABibb3Hmnay84/3yfnjDPnXOOa/66++6ku9pEqsxIZXjlI8BMYC8RWSoi40RkvIiM93Z5EVgILADuAs6PeuwbwBPAEd5jj/b9CEzopXWyzZYtrk151Cjo2dOHJzR06+b+nnff7ZYwTsAmUmVGKqNuTlXVzqparKpdVPUeVb1DVe/wylVVL1DVPVR1P1Utj3rsYFXtqKrbe499OZ0HY8Ip1ck2TerUe+op+PpruPBC3+I1wO9+5yaerUh8cm8TqTLDFjUzodDk1RHvuQf+539cR6xNkvJPJK/YCKaMsdUrTeg1a3VEW60yfVavhi++gP79k+9rmsVWrzSh1+hOvaoqd1FrVUvy6XTqqe5Ua8uWoCPJa/YJN6HQ6E69f/3LXfh72rS0x5bXLr4YvvoKHnigzmabJJVZlo+qMjQAABFlSURBVOhNKDSqU2/zZjdDp29fOPLITISXv446yq1sOWnSthE4ySZJ2ZeA/yzRm1Bo1OqId93lJvNce611FqabCFx3Xe3fm8STpGymbHpYojf5ZcUKuOIKGDYMRowIOpr8MHw4/OxnLtmrJuxPifUl8P77sNdedWv3VutvHEv0Jr989ZXLDnfcYbX5TLrvPtdOL5KwPyXepR9Va2v306ZZrb+xbHilyT/V1VBYGHQU+Wn+fCo/W0qXcUfXubxsaSksWgTPPQfnnecu3RtLQYFL+iJ2ofH6bHilCZXo0/Ybb3QLJSY9hV+50rURb9liST5I559P61+extr5y2P2pyS79GNNjdvP1sdpHEv0JivFa4Ot31l3ySVw2WUNT+HrPP7GGvScc93yuV9+Gehx5b3bb4eNG+HEE2OOrY/uVH/oodhNPKecYuvjNJqqZtXtkEMOUZPfPv9c9cADVVu3dvW91q1VDzrIbe/YUbWgILouWPdWUKDatm3dx99QfLkq6MrL/h70oRlV1ccfd2/M+PEJd1u3TrW0tO77W1qqunhx7O3r1mUo/iwFlGucvGpt9CbrJFrOoE8fmD498eOLi12NvqYGzuQ+7uNs7uIc/tjhn6xcZR2wWeHSS+H66+Hpp+FHPwo6mlCwNnqTUxINv4s1YiNaSYlbbbimBkpZy01MYBojOJ/b2Hc/S/JZY9IkuPVWOOaYoCPJC5boTdZJNPwuWWddUZG7WFRJCayjLcN5jZN4gu1Kiq0NN5sUFsIFF0CLFm5uw9SpQUcUapboTdZJtJxBrBmw9UdvnFb9AOdV/wOA9ziE9ZTaGufZ7Ior4Mc/dmPtTVpYojdZp1HLGUSrroarr6bV+Wdyw+Dn0eqaxj3eBOPmm+GII9wp229/m/SqVKbxLNGbcJg/382fnzgRTj/dzbyx5YdzQ+vW7v264AK46Sa3CNqiRRkPI8zLKth/gkmbjP3jfPutu7DFF1/Aww/D/fe7tl+TO1q2dJ2zL7wA7dq5oVeQ0ofGj89Z6BdTizfuMqibjaPPbVVVqjfe6MY177JL7Vj2Vq3c/dJS1cmT3X7N8vXXqnfcUXv/iSdUV65s5pOarFBT435u3Ki6996ql1+uunRpzF0TzbmIFvlctm+vesMNqtdf736PfBZjzc8oKHDbcwUJxtEHntjr3/xO9NFvsC8JJgB+HUO6/xb1/+ni3eL9Mya1dq3qAw+ojhmjWlysKqL62Wf+HoTJHitXqp58svvQiKiOGOHe//Xrt+2SSoKu/7kUqX1M5LPYr1/sz+qwYQEcdxPlbaJP9ds+m/l1DH48T7IvimSzVmP9M8Z9zo0bVefNU122zN2fPl21qMg9uFs31YsuUv3008b9EUxuqqhQvfJK1d13d+//rFlu+wcf6DV7P6QH8L62ZGOdz1dRUeLaev3P4g47qJaU1N1eUqL60EPBHnpjJEr04ZoZW+9Yts2wVAABlAKJs8qdSIPHp1JeXQ1TpsBfrxcuu1SZMCHOmlmJnj+yXG6M8p12gm/XiDeByJU3OIYEj4+U77QTfLtavb8FdZ8nMls0epZS9HMVFlJRAaeeVMWCCuWHH5SSVsqevZSH/l1Arz6uPXzU4A28/WY1BdRQSDVFVLGZlqyjLaDszzxaspnWVNKaSg7a43sWFvZk6rJDKaxcz92F4+lSuJzOVUvorosQVTd78g9/cG/klClujGS/frbEcD5ShdmzXWN6URFcdJH7TADVFLCULiyhK0fwKltoybEtX2FAxy/Y2LKUuV+UsIESvmcH5uAmj3ZkJduxia0U07d/ER99Usjq74rZwA4A7NRmM59/prQpldrPW0FB7djf6M6A6M9jZBBA/Vl/kX1S+H9tikQzYwOvwde/NatGH1lDo96tPzMVVM/inthf6R995B5/002xyxcvduXXXBOzfNdWaxRUJxddErO8U7stOnmyavV55zcsb9myNv6f/7xB+driDtvuPs0JDR+/2261jx8xomH5/vurqurQoaoz6d+g/MMdD6t9fO/eDR9/9NGq6mpFi+naoHxqixO3PXxTq9IG5Xdz9ra7VTSsVt3Cr1wNjC36OT11OkP035yqk4qv0ku7P6wLX13Y9M+DCbctW/S72R/rGa0e14lcqQ/wc32BUds+Xg9yWoPP22rabbv7JD9u+Hnv3r32+WP9P+27b235gAENywcOrC3v06dh+VFH1ZZ3796wfOPGJv85SFCjTzDHMAf17u26zT0ffADPvwBLt3QB4H0OYlKLiYw5Bg44oPZh1e06MmUyTPtLf648ciKq8NZbMHgwDBgABZEB2IMH13n+G26AHzbCtz9sB8BLVSP4ju0BKI76y65aU8BVV8GiTscw8Vc70b59VMzR1f8f/Qj22KPOIS38tBUlz7v1uR/lFOZyIC1aUHsMpaW1O59+OgwaVPdv0qkT4Gab/mvmL3hp86htRS1awPCxXWv3vfBCWLWq9r7Itnj69IHJ0y9mR75DkW234j335jhv95qJf+HyK7dSuamAGgrYSjFLW+/N4o9drD9Z9xRbaEElrSksacUue5bwwns7A1BFMXsSNcRhKxQsgXtOye81xk0CxcXs0G8f7q/cBziJYcPqroP0C+7mD9zA8EO+Y9lnG9ANGyiialv5bVzAi4ymzfZbmXR1Fdu3qK47JXvcOHd1LPVq3qrQsWNt+TnnwOjR28prauCVj7twWge3oupvz7uAgm9X1T4W6v5/T5gA69fXPaZE076bI943QFA3P9vo461+F73KXSodNfHasYcOjVmBj3trSi9+Ksfgx/Mka39/6CH/2jAjr1VS4k5oEv3NcqkzzATLz89oY2VDfyD52hmbilQ6auIl51gfrGS3TCeuVEbapPIhTfRF0ZjRPKmOzMnkP6kJh6ZUivwaiZYNwzPzItE39Q1LpVYeLznH+mBlU+JKtZbRnA9pY2sy8V6rfXtbY9z4K1lO8LMWHi+P1M8d6RziHPpE35w3LFmtPNXknErSz1TtIiLVBJ7qh7Q5r+HHaxmTqlRygp+18FSajdLdvBP6RN+cNyxZgs5UrbL+h6A5M0kjXxiRYefJkmpz2jYbm7iDbEc1+SOVnJDos9vYSlcqzUbpbt4JfaIPQy0xUV9BY775k7WBx0qqzenwbWzi9qtz2ZhEUskJ8T67N96Ynpp3uvNUsxI9cC+wEvgoTrkAtwALgHnAwVFlZwAV3u2MZK+lTUz0maglpqttLVntu7Hf/Mk6lxubVJMdtyVuk41SyQnxPrsdOjS/5h3r/ybdeaq5iX4IcHCCRD8aeMlL+AOA2d72dsBC72db7/e2yV6vKYk+3ckmXW1rjRmBkuo3v5+1hkTHHYY1hEx4NScnNPd/KN7/zZw56c1TzW66AXokSPT/BE6Nuv8Z0Bk4FfhnvP3i3dI5vLKpySldbWuNWRsm1jd/umsN8Y67bdvgxwwbky6x/odatHDbUskb8f6vRdJbKUp3on8eGBR1/1WgL3Ax8Keo7VcAF8d5jnOBcqC8W7duafkjNKdW7nfbWrLmmsGDGz/Rqym1hmRffPGOu7g4+DHDxqRLogEaqeSNREO201kpyvpEH31LV42+ObVyP2vJTeks9ft4YsUR6wMY77j32cffLz5jslWy/7NUz6qT/Z/60RSaKNH7cYWpZUDUgil08bbF2x6IPn0aLiZXUwP77hv/MZEr1/z611BVVbesqRebLiuDefOgsjJ2earPG+941q5N7So79eOorHRrA5WV1e4T7yLdEybUXRIE3P2zz04etzG5JFHeiHdVqt69Ey9ZUz/vZOTqVvG+AaJvJK7RH0Pdzth3vO3tgC9xHbFtvd/bJXutdNXoG9vulq4OWL+agRLVGppzeplKHDbSxuSLRGfzqZxVp9Ia4FcfIM0cdfMIsALYCiwFxgHjgfFeuQC3AV8AHwJ9ox57Nm7Y5QLgrGSvpWlM9I1td0tXB6xfzUDJJnoli9UmLhmTXKz/szZtVP/yl9QmJKZSKfKr8tesRJ/pWyYWNWvurLnmSEdtuCmxWq3cmMbzq49NNfEqrk2pdCVK9OFajz5FffrUXbcaGrabjRsH5eVuHfgIP9qh27Rx7eh+akqs6YjDmLArK/OuWlcTuzzVPraKChg71v2M1V/X1D7AePzojM06kU7UDh1id0yOG5e8MzFeR6Sff3y/5FKsxuSyWJ2zAMOGubr42rWuEpVMrEEZBQXuuiaNeZ5UhS7Rp9KDnUpijNR4o0+o/P7j+yWXYjWmqZJV4DIhlUpiKpoyCrA5QpfoUxk2mMuJMRs+7MZkWkaGIKbAr7Nnv74wUhW6RJ/pb8pMypYPuzGZlkoFLhP8qiRmurk1dIk+09+UmZSpD7udNZhsE7YKXKZbFUKX6MPcMZmJD7udNZhsFOYKXCaELtFna/u7H7XkTHzYs+UU2ZhoYa7AZULoEn028quWnIkPe9hOkU04ZGsFLldYom+CxtbO/aolZ+LDbqfIxoSPJfpGakrtPJdqyXaKbEz4WKJvpKbUznOplmynyMaEjyX6RmpK7dxqycaYIFmib6Sm1M6tlmyMCZIl+kay2rkxJtfk5TLFzWHL+xpjco3V6I0xJuQs0RtjTMhZojfGmJCzRG+MMSFnid4YY0LOEr0xxoScqGrQMdQhIquAxUHH0QQdgNVBBxGAfDzufDxmyM/jzqVj7q6qHWMVZF2iz1UiUq6qfYOOI9Py8bjz8ZghP487LMdsTTfGGBNyluiNMSbkLNH7586gAwhIPh53Ph4z5Odxh+KYrY3eGGNCzmr0xhgTcpbojTEm5CzR+0REficiKiIdvPsiIreIyAIRmSciBwcdo59E5EYR+dQ7tv+ISGlU2WXecX8mIkcHGaffRGSkd1wLROTSoONJBxHpKiL/FZGPRWS+iPzG295ORKaJSIX3s23QsaaDiBSKyPsi8rx3fzcRme2954+JSIugY2wsS/Q+EJGuwFHAV1GbRwG9vNu5wO0BhJZO04B9VXV/4HPgMgAR6Q2cAvQBRgL/IyKFgUXpI+84bsO9t72BU73jDZsq4Heq2hsYAFzgHeelwKuq2gt41bsfRr8BPom6fz0wRVV7AmuBcYFE1QyW6P0xBfgDEN2zfTzwoDqzgFIR6RxIdGmgqq+oapV3dxbQxfv9eOBRVd2sql8CC4B+QcSYBv2ABaq6UFW3AI/ijjdUVHWFqr7n/f49LuntijvWB7zdHgBOCCbC9BGRLsAxwN3efQGGA096u+TkcVuibyYROR5Ypqof1CvaFVgSdX+pty2MzgZe8n4P83GH+dhiEpEewEHAbKCTqq7wir4GOgUUVjrdhKu01Xj32wProio1Ofme26UEUyAi/wfsHKPoj8DluGab0El03Ko61dvnj7hT/X9nMjaTfiJSAjwFTFDV71zl1lFVFZFQjc0WkTHASlWdIyJDg47HT5boU6CqI2JtF5H9gN2AD7x/gi7AeyLSD1gGdI3avYu3LWfEO+4IETkTGAMcobUTMnL+uBMI87HVISLFuCT/b1V92tv8jYh0VtUVXjPkyuAiTIsy4DgRGQ1sB+wI3Ixrdi3yavU5+Z5b000zqOqHqrqTqvZQ1R6407qDVfVr4FngdG/0zQBgfdRpb84TkZG4U9zjVPWHqKJngVNEpKWI7IbrjH4niBjT4F2glzcKowWu0/nZgGPyndcufQ/wiar+ParoWeAM7/czgKmZji2dVPUyVe3i/S+fArymqj8D/guc6O2Wk8dtNfr0eREYjeuM/AE4K9hwfHcr0BKY5p3NzFLV8ao6X0QeBz7GNelcoKrVAcbpG1WtEpELgZeBQuBeVZ0fcFjpUAb8HPhQROZ62y4H/go8LiLjcEuJjw0ovky7BHhURK4B3sd9CeYUWwLBGGNCzppujDEm5CzRG2NMyFmiN8aYkLNEb4wxIWeJ3hhjQs4SvTHGhJwlemOMCbn/BzvAWieNlVjIAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from __future__ import print_function, division\n", "from numpy import arange, sqrt, exp, pi, random, ones_like\n", "import numpy as np\n", "import matplotlib.pylab as plt\n", "from PyAstronomy import funcFit2 as fuf2\n", "import scipy.optimize as sco\n", "\n", "random.seed(1234)\n", "\n", "# Creating a Gaussian with some noise\n", "# Choose some parameters...\n", "gPar = {\"A\":1.0, \"sig\":10.0, \"mu\":10.0, \"off\":1.0, \"lin\":0.0}\n", "# Calculate profile\n", "x = arange(100) - 50.0\n", "y = gPar[\"off\"] + gPar[\"A\"] / sqrt(2*pi*gPar[\"sig\"]**2) \\\n", " * exp(-(x-gPar[\"mu\"])**2/(2*gPar[\"sig\"]**2))\n", "# Add some noise...\n", "y += random.normal(0.0, 0.002, x.size)\n", "# ...and save the error bars\n", "yerr = ones_like(x)*0.002\n", "\n", "# Create a model object\n", "gf = fuf2.GaussFit()\n", "\n", "# Set guess values for the parameters\n", "gf.assignValues({\"A\":3, \"sig\":3.77, \"off\":0.96, \"mu\":9.5})\n", "\n", "# 'Thaw' those (the order is irrelevant)\n", "gf.thaw([\"mu\", \"sig\", \"off\", \"A\"])\n", "\n", "# Restrict parameter range of sigma to within 4.1 to 7.1\n", "# where the upper limit is too small for input values\n", "gf.setRestriction({\"sig\":[4.1,7.1]})\n", "# An upper limit of 10 with no lower limit on the area of the Gaussian\n", "gf.setRestriction({\"A\":[None,10]})\n", "\n", "# Check restrictions \n", "print(\"Parameter restrictions: \", gf.getRestrictions())\n", "print()\n", "print(\"Penalty factor: \", gf.getPenaltyFactor())\n", "\n", "# Inspect the effect of the restriction on the objective function value\n", "sigs = np.linspace(3,8,100)\n", "oval = np.zeros_like(sigs)\n", "\n", "for i, sig in enumerate(sigs):\n", " gf[\"sig\"] = sig\n", " oval[i] = gf.chisqr(gf.freeParamVals(),x,y,yerr)\n", " \n", "plt.title(\"Value of objective function (with 'penalty ramps')\")\n", "plt.xlabel(\"Sigma\")\n", "plt.plot(sigs, oval, 'b.-')\n", "plt.show()\n", "\n", "# Use a minimization algorithm, which itself ignores\n", "# restrictions, and observe effect of penalties\n", "fr = fuf2.fitfmin(gf, gf.chisqr, x, y, yerr=yerr)\n", "print()\n", "print(\"Fit result: \", fr)\n", "print()\n", "\n", "gf.parameterSummary()\n", "plt.title(\"Bad fit caused by restrictions\")\n", "plt.plot(x, y, 'bp')\n", "plt.plot(x, gf.evaluate(x), 'r--')\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.10" } }, "nbformat": 4, "nbformat_minor": 4 }