{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Gaussian fit and parameter handling\n", "===================================" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAZsklEQVR4nO3dfawc11nH8d8PlzR/lN6mtSkljn0NGFTzIre6MpUisEUhdUNVQ0EogUIKkaKIBvFSBAlF9XWgAgTiTYSiFKK2CAhRC9RCQSFtc8MfNOCb5gWSktaEtLEbqFGpQQK1JDz8sbP1eDqzO7s7szNz9vuRru7u7O69Z3Znn3nmOWfOOCIEAEjXl3TdAABAuwj0AJA4Aj0AJI5ADwCJI9ADQOKe13UDinbu3Bnr6+tdNwMABuXBBx/894jYVfZY7wL9+vq6tre3u24GAAyK7U9UPUbpBgASR6AHgMQR6AEgcQR6AEgcgR4AEkegB4DEEegBIHEEegBI3EoE+iNHRj8AsIpWItCnoKmdFTs9YPUQ6AEgcQT6BpAlA2lI9btMoG9YcUNJdcMBMBwrE+ifemry4/mATHCuxnsDDM/KBPpPVE7gORzTdlYAutXXRGhlAn0KmtpZVe0wGNkDpCmpQF8MMJubki3df//ovj362dxs/n/3Odsuvi/LPLoh6KOPJm2XbX+Xu/hOJBXoizY3pQjp8OHR/YjRTxuBPh88ixvKLBtOlzsrAGmUeYt6dynBLuUD8iJ79eKGMu+GMw74EaPb998/uj3La7e2Rr+femq0gxgb3z5+fPadxrzvTbFNy349sKqSzujz9u6d/px8QK4TnMuy7fHtvllfb+7oJsWMB+mbVDJp8si5j2XclQn06+uL/42ysko+eOaNN5Q6G07dml2dndWiFsnW8+vQx40dqNJkmbePiVBygX7WADMpK591r162oTTZP7DIzir/vkzaYVRtpLNmPJP6LIA+GG+XVYnWop2mXXb4FiVXo5+2Ny3WeTc3Rz/5Gvis9XCp2Wx71s5baXrdOv++zLrDWKSvoPi/gb6o2i7n+S5vbkonTly4P+4Dq/pby/5OJJfRdyUfPIsfbtmHPWlvv+xsuKn65LjDd1KfxaLrw9EB2jbPkXNZ6efw4Qt/q+thxklk9FV70/GIkrqjNPIBeZEMvbihLFJyWWTP//DD5SNt8us2T7Ze9d6sr48C8fjv5OXbsQiODjA2zyisOt+Jpv5X0bSRb22OKpua0du+w/anbf9jxeO2/du2T9t+1PYrc49dZ/vj2c91TTY8r6mOlHxArtoTt5VRjjfAqmx4a2v6BlBs6/nzk7OMedV9b8r+N9ClNr4Txe/BpERoWef1FNUp3bxL0tEJj79W0v7s5wZJ75Ak2y+WdFzSN0s6JOm47csWaWzT5jmcaiujLG6AeU2eJDVth7HIkUz+vSn+nUVKQ5w0hrZVJXDzJHb5nUZfSo1TSzcR8Te21yc85Zik90RESHrA9otsv0zSEUn3RsRnJMn2vRrtMP5k0UZPsowhiEVtHGrNUkopnuhVdni4tlbv/1ZlNrNusMU+i3xJZ56Tvood5kBdVaXd/HdinKQUv8uLJnZNdvguooka/eWSns7dP5Mtq1r+RWzfoNHRgPbs2bNQY+ocgs0atGY5q7S4oUzaCVTVC+sG5bHiiJpiUG2iE6hqg616b/Ibctlnkg/inPGKJtUdWbds+bjTxHk9s+hFZ2xE3C7pdkna2NhoPV8rC1r5IFMMOGXBswnjcs2kDbBvJ0nVeW8mfYkWXZ8ujtjQX4uURqqSlMOHL+4fmzSIoSw5qfq7ZcqmKmlDE8Mrz0q6Ind/d7asavlSNDmcqcs6W9Wev+pEr4cfvvCcukFxkZOkFinpzGPZmRD6bbzt1j3pqXikWdY5urU1+0CC8fdga2t0e96BCG31ATYR6E9K+qFs9M2rJJ2PiGck3SPpKtuXZZ2wV2XLOlEnMFaZ1MnYVDtmLddUjds9f/7Cc2YNisWO2jqjmcremzojhIAutDXUedJjZd/3ZSePdYZX/omkD0v6OttnbF9v+0bbN2ZPuVvSk5JOS3qnpB+VpKwT9hckncp+bh13zLZlUoCpExjrWDSjrGrHwYP1/8asRytl70sbI1nK3huCPppWlSzNWoLMb5dVCVwTpcayZKlsAEWbo8rqjLq5dsrjIenNFY/dIemO+ZrWrVlPaW7KPEFxbW32k0DKOqgmyf+tpt6bOtNC92V4GvqjTudqkyXFSd+taYM06iSGbfUB5q3kFAjjwDhpDzrtlOYmTAuMdbPhgwfbb2v+bzX13pRNC93l1bCQjqa2m6o6/nhbn+UEqKrv+zKSmV6MumlL1Rt48OCFYX1djstus2NxliODqrl4Zv070yz6tygBoSi/7U4LmE1vP7PuTKq+7032AVZJOtA3tVfPdzIuoip4NrkBzrOhLLLDmfW9qSr75G9XlaLmuRoW0pbfdj/xiWbOS6ky6bvVVIBuK/lbydJNXp0PaEhD+pbd1ln/X1XZJ3+7rBS1rDlBMGzzDLioKpEWl6+vz3ZVuWmDQ5Y5Gie5QD/raJJ5A2PX047mx+1S0kDKpl3Ao27gbUJT/VPL6APMS650M200SV+C4qLtWGYnZVvvWdW00HWvhgVUnZldNb1IV+W/ruNOcoEes5m0AbY9GiCfvRRrrWXLgbxJ2+4yBlw0lYQ01Qc4SdKBvokPouk3f5HgOe0CK02bdtTQZZbCRGgoavPor2w7ayoJWUYyk3Sgn/YGdhEkFim5DKUsNY86s2ACk1R93/t4kuOyJdcZu0wpnrXZ1UU+yk5MabNzCsOyyHeNbSjxjL5tdbPzNkoubWUps06NMK9J00JPW47VM8uRMNvNFyPQL0EbwTPlLIVyDfquqZ3JsnZKyZZu2hpfvkrXL+0q4NbdiaVYOsPIkSOj7aDud43zSSYjo5/Rotn5kLLVvh81MOFZOspGUTUxqyPBfyTZjL6v+h48AaSHQL+ArrLzVT5MXaXSGS4Y0pFwH1G6WQDZefOm7cCWNSoI/cJ3bTFk9Bi8rieYQzvobG8OGT1KDaE0tKwLpmA5ioGdzvbmEOiXiODTLA7n00Jgbw+Bfk4EbaB5k84ix/wI9AA6Uwzs49t7944yfDrbm0FnLJJAx11/TeosL7vSUgRluaaR0WOQiqWzfH2XoJ8Gxs43h4weyaFTb5iKgZ2svjkEegxW2VmybV4YGu0isLeHQI/BKtZ385gaAbiAGj2SEcHUCH1Fv0m3CPRIAh13/TZrvwnnqTSLQI8k5Ou7TI3QX1WfA59Luwj0SA6dev1QdZYrR1/LV6sz1vZR20/YPm375pLH99r+oO1HbW/Z3p177DnbD2c/J5tsPID+KjsZ6vBhdsRdmJrR294h6TZJ3yHpjKRTtk9GxOO5p/2apPdExLttf5ukX5L0g9lj/xMRBxtuNwCgpjqlm0OSTkfEk5Jk+05JxyTlA/0BST+V3b5P0l802UhgEuq7/ddUvwl9LfOpU7q5XNLTuftnsmV5j0h6Q3b7uyV9me2XZPcvtb1t+wHb37VQa4EpJl1mkSF+3cmXa/gclq+pE6Z+WtJh2w9JOizprKTnssf2RsSGpO+X9Ju2v7r4Yts3ZDuD7XPnzjXUJOBiTI3QD3wOy1cn0J+VdEXu/u5s2RdExKci4g0R8QpJb82WfTb7fTb7/aSkLUmvKP6DiLg9IjYiYmPXrl3zrAcAoEKdQH9K0n7b+2xfIukaSReNnrG90/b4b90i6Y5s+WW2nz9+jqQrdXFtH2hV2Xw4TI2wfE89VT4vEWWc5Zga6CPiWUk3SbpH0kcl3RURj9m+1fbrs6cdkfSE7Y9Jeqmkt2fLXy5p2/YjGnXS/nJhtA7Qqqr5zgn0y7O1NQroTQ21ZOcwu1onTEXE3ZLuLix7W+72eyW9t+R1fyvpGxdsI4BEzRO0qfHPjtkrsTI4I7Mf8p8DQXs5CPRIVvESdpyR2Q/zfA70tSyGQA9gqco6ZqcFbfpaFsOkZgCWan19FOy5dsDykNEDGBT6WmZHRo+VwNwo/VD8HOYJ2vS1zI6MHkkrG75X7KRFdwjay0GgR9IYvgcQ6LHiyO6xCgj0SA5jrtM1aRpqVKMzFsnZ3Bz9MHyv3wjYy0NGj5XExFjNowzWXwR6JK1q+B6dtFglBHokrc7wPbJ7pI4aPVbG5qZ04sSF+3ZnTUkeF/HuFzJ6rIyyibHGtyVqzEgXGT1WVj6jH99mHpXFUAbrJzJ6JGvSmOu9e5u7tB0uoJO7nwj0WEkEdKwSSjdYeWXlGjoT6zlyZFSuyWfy4zLY2loXLUIZMnqsvHx2T415duvr5WWw8+c7bRZyCPRADjVmpIjSDVYO5Zj2rK2Vj2Y6fpxJ5bpEoMfKG1+seoyhlvM7eHC0I2VCuX6hdIOVtrU1CvSTTqTCZPRr9B+BHqhAvb6eqveJI6L+INADOQSn5nCuQn8Q6IECrk41XdlVvO6/n7mC+opAD+SUjQmPINAXlU0Qx/vUXwR6AEgcwyuBCtTr6yl7nzhXoV/I6IEKdCaWK87bz/vUf7UCve2jtp+wfdr2zSWP77X9QduP2t6yvTv32HW2P579XNdk44EmTZrWGBiyqYHe9g5Jt0l6raQDkq61faDwtF+T9J6I+CZJt0r6pey1L5Z0XNI3Szok6bjty5prPtA+rjyFoauT0R+SdDoinoyIz0u6U9KxwnMOSPpQdvu+3OOvkXRvRHwmIv5D0r2Sji7ebABAXXU6Yy+X9HTu/hmNMvS8RyS9QdJvSfpuSV9m+yUVr728+A9s3yDpBknas2dP3bYDraGEUw/v0zA01Rn705IO235I0mFJZyU9V/fFEXF7RGxExMauXbsaahKAtjC/zbDUyejPSroid393tuwLIuJTGmX0sv0CSd8TEZ+1fVbSkcJrtxZoL9Aprjw1wjxAw1Inoz8lab/tfbYvkXSNpJP5J9jeaXv8t26RdEd2+x5JV9m+LOuEvSpbBgwKGSyGbGqgj4hnJd2kUYD+qKS7IuIx27fafn32tCOSnrD9MUkvlfT27LWfkfQLGu0sTkm6NVsGDAoZbPn8NswDNAyOnl0ZYGNjI7a3t7tuBvAF+YtoULrhoiJ9ZfvBiNgoe4wpEIAKm5vSiRMX7nPlKQwVUyAAFcpmaOTKUyPs7IaFQA/MiHo989sMDYEeqIEMFkNGoAdqYsQJhorOWKCG9fXRWHpGnGCICPQAZrLKQ0uHitINMKNVqdczPXM6yOiBKYoZLCNOMDRk9ACmIrsfNgI9ACSOQA8AiaNGD8xg1UacMD1zGsjoAVRiuoc0EOgBIHGUbgBcpGp65rW1TpqDBpDRA7hI1fTM58932iwsgEAPAIkj0AOotLbGrJ0poEYPoNLBg6MhpczaOWxk9ACQOAI9gC+YNKfNqszamSICPYBamLVzuKjRAyi1atM9pIyMHlhBTDu8Wgj0wJwIlhgKSjcALlI2YyVlnGEjowdwEWasTA+BHgASR6AHVlS+RLO5yVQHKSPQAwsY8hWY8iWashkrIwj0qSDQAwugno0hqBXobR+1/YTt07ZvLnl8j+37bD9k+1HbV2fL123/j+2Hs5/fa3oFgD7q69DLOiUapjpIz9RAb3uHpNskvVbSAUnX2j5QeNrPS7orIl4h6RpJv5t77J8j4mD2c2ND7QY6M+R6dp0SDVMdpKdORn9I0umIeDIiPi/pTknHCs8JSS/Mbq9J+lRzTQT6pSpYjqfzBfqmTqC/XNLTuftnsmV5m5LeaPuMpLsl/VjusX1ZSed+299S9g9s32B72/b2uXPn6rcewNwo0ayOpjpjr5X0rojYLelqSX9o+0skPSNpT1bS+SlJf2z7hcUXR8TtEbERERu7du1qqElA+4YcLCnRrI46UyCclXRF7v7ubFne9ZKOSlJEfNj2pZJ2RsSnJX0uW/6g7X+W9LWSthdtONAHqQVLpjpIU52M/pSk/bb32b5Eo87Wk4XnfFLSqyXJ9sslXSrpnO1dWWeubH+VpP2Snmyq8UCfDXmMPdIyNaOPiGdt3yTpHkk7JN0REY/ZvlXSdkSclPQWSe+0/ZMadcy+KSLC9rdKutX2/0r6P0k3RsRnWlsboGP54M4Ye/RFrdkrI+JujTpZ88velrv9uKQrS173PknvW7CNwGAMJbhTolktTFMMzKkqWNpffPv48f6Osx8PCSX4p4tADyxoc1M6caL8sYilNgUoxVw3wILKTqAa3wb6gEAPtGTIY+yRFgI90KB8cB+Pse/rBGdYHQR6oEF9OoFqlh0MY/7TRqAHMJhhoZgPgR4AEkegB1qwtVU+Lr1P9fohz6uP2TCOHmjIpBOO+lgD39wc/Rw5Mgr2jPlPFxk9sARlNfBlZPd93MFg+Qj0QMLqdrIy5j9tlG6AlhSnRhjPezNLUF3WPDR9GhaK5pHRAy3pamoEOllRRKAHlqztMetVFy8n0K8uAj2wBNTA0SVq9MCSlM1Tv7bW7v+ss4NhHvr0kdEDS7C+Xl6vP3++/f8LEOiBnimOr2csPBZF6QZYsrW18jLO3r3lQZ0Jx7AoMnqgZcV5bw4eLC/jzFtmyR8B9GkuHfQHGT3QQ089tfhFxulkxRiBHuhQ1aiY9fVRsGfCMTSB0g3QoapyDR2waBIZPbBEdacyznfAcrIVFkVGD/RE1eiaOp20+Z0ERwMoIqMHeqRq2GVRcVbL/E6iaoexrJkw0T8EeqAjW1sXZposE8FQSTSD0g3QoUWmMh4PwcxPR8zUxChDRg/00LhcM6nMUjYEk+GYKENGD/REvhbPZGRoUq1Ab/uo7Sdsn7Z9c8nje2zfZ/sh24/avjr32C3Z656w/ZomGw+kZNbgnh9dk99JTBqOyYic1TS1dGN7h6TbJH2HpDOSTtk+GRGP557285Luioh32D4g6W5J69ntayR9vaSvlPQB218bEc81vSLAqsmPrsnvJCbtMJggbTXVyegPSTodEU9GxOcl3SnpWOE5IemF2e01SZ/Kbh+TdGdEfC4i/kXS6ezvAQCWpE6gv1zS07n7Z7JleZuS3mj7jEbZ/I/N8FrZvsH2tu3tc+fO1Ww6kKbibJdjR46MsvWykTaTSjJcLBxNdcZeK+ldEbFb0tWS/tB27b8dEbdHxEZEbOzatauhJgHDURXci6quVDWpXMPFwlFneOVZSVfk7u/OluVdL+moJEXEh21fKmlnzdcCWABnumKaOln3KUn7be+zfYlGnasnC8/5pKRXS5Ltl0u6VNK57HnX2H6+7X2S9kv6+6YaD6yystE1044MmCBtNU3N6CPiWds3SbpH0g5Jd0TEY7ZvlbQdESclvUXSO23/pEYds2+KiJD0mO27JD0u6VlJb2bEDdCMecbaMz5/NdU6MzYi7taokzW/7G25249LurLitW+X9PYF2gggwzh4zIMzY4EBYRw85kGgB4DEMakZ0HObm9KJExfu5y8UPgtG56wuMnqg5xgHj0UR6AEgcQR6YEAYB495EOiBAWEcPOZBoAeAxBHoASBxBHoASBzj6IGBYBw85kVGDwCJI9ADQOII9ACQOAI9ACSOQA8AiSPQA0DiCPQAkDgCPQAkjkAPAIlzRHTdhovYPidpiFfG3Cnp37tuRAdWcb1XcZ2l1VzvIa3z3ojYVfZA7wL9UNnejoiNrtuxbKu43qu4ztJqrncq60zpBgASR6AHgMQR6Jtze9cN6MgqrvcqrrO0muudxDpToweAxJHRA0DiCPQAkDgCfUNsv8V22N6Z3bft37Z92vajtl/ZdRubZPtXbf9Ttm5/bvtFucduydb7Cduv6bKdTbN9NFuv07Zv7ro9bbB9he37bD9u+zHbP54tf7Hte21/PPt9WddtbYPtHbYfsv2X2f19tv8u+8z/1PYlXbdxVgT6Bti+QtJVkj6ZW/xaSfuznxskvaODprXpXknfEBHfJOljkm6RJNsHJF0j6eslHZX0u7Z3dNbKBmXrcZtGn+0BSddm65uaZyW9JSIOSHqVpDdn63mzpA9GxH5JH8zup+jHJX00d/9XJP1GRHyNpP+QdH0nrVoAgb4ZvyHpZyTle7aPSXpPjDwg6UW2X9ZJ61oQEX8dEc9mdx+QtDu7fUzSnRHxuYj4F0mnJR3qoo0tOCTpdEQ8GRGfl3SnRuublIh4JiI+kt3+L42C3uUareu7s6e9W9J3ddPC9tjeLek7Jf1+dt+Svk3Se7OnDHK9CfQLsn1M0tmIeKTw0OWSns7dP5MtS9GPSPqr7HbK653yupWyvS7pFZL+TtJLI+KZ7KF/lfTSjprVpt/UKGn7v+z+SyR9NpfUDPIzf17XDRgC2x+Q9BUlD71V0s9pVLZJzqT1joj3Z895q0aH+n+0zLahfbZfIOl9kn4iIv5zlNyORETYTmpstu3XSfp0RDxo+0jX7WkSgb6GiPj2suW2v1HSPkmPZF+C3ZI+YvuQpLOSrsg9fXe2bDCq1nvM9pskvU7Sq+PCCRmDX+8JUl63i9j+Uo2C/B9FxJ9li//N9ssi4pmsDPnp7lrYiislvd721ZIulfRCSb+lUdn1eVlWP8jPnNLNAiLiHyLiyyNiPSLWNTqse2VE/Kukk5J+KBt98ypJ53OHvYNn+6hGh7ivj4j/zj10UtI1tp9ve59GndF/30UbW3BK0v5sFMYlGnU6n+y4TY3L6tJ/IOmjEfHruYdOSrouu32dpPcvu21tiohbImJ39l2+RtKHIuIHJN0n6Xuzpw1yvcno23O3pKs16oz8b0k/3G1zGvc7kp4v6d7saOaBiLgxIh6zfZekxzUq6bw5Ip7rsJ2NiYhnbd8k6R5JOyTdERGPddysNlwp6Qcl/YPth7NlPyfplyXdZft6jaYS/76O2rdsPyvpTtu/KOkhjXaCg8IUCACQOEo3AJA4Aj0AJI5ADwCJI9ADQOII9ACQOAI9ACSOQA8Aift/TA3il/LwTMYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from __future__ import print_function, division\n", "# Import numpy and matplotlib\n", "from numpy import arange, sqrt, exp, pi, random, ones_like\n", "import matplotlib.pylab as plt\n", "import scipy.optimize as sco\n", "# ... and now the funcFit2 package\n", "from PyAstronomy import funcFit2 as fuf2\n", "\n", "# For reproducability\n", "random.seed(1234)\n", "\n", "# Create some mock data\n", "# Choose parameters...\n", "gPar = {\"A\":-5.0, \"sig\":10.0, \"mu\":10.0, \"off\":1.0, \"lin\":0.0}\n", "# ...and 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.01, x.size)\n", "# ...and save the error bars\n", "yerr = ones_like(x)*0.01\n", "\n", "plt.errorbar(x, y, yerr=yerr, fmt='b+')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fitting and parameter management\n", "---------------------------------------------------------------" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "List of available parameters: ['A', 'mu', 'sig', 'off', 'lin']\n", "\n", "\n", "Parameters and guess values: \n", " A : -10.0\n", " sig : 15.77\n", " off : 0.96\n", " mu : 7.5\n", "\n", "------------------- Parameter summary --------------------\n", " A = -10, free: F, restricted: F, related: F\n", " mu = 7.5, free: F, restricted: F, related: F\n", " sig = 15.77, free: F, restricted: F, related: F\n", " off = 0.96, free: F, restricted: F, related: F\n", " lin = 0, free: F, restricted: F, related: F\n", "----------------------------------------------------------\n", "\n", "Parameter name/value dictionary: OrderedDict([('A', -10.0), ('mu', 7.5), ('sig', 15.77), ('off', 0.96), ('lin', 0.0)])\n", "Names and values of FREE parameters: {'A': -10.0, 'mu': 7.5, 'sig': 15.77, 'off': 0.96}\n", "Names and values of FROZEN parameters: {'lin': 0.0}\n", "\n", "Optimization terminated successfully.\n", " Current function value: 98.642163\n", " Iterations: 154\n", " Function evaluations: 261\n", "\n", "Fit results: (array([-5.02400025, 9.88993313, 10.0142805 , 1.000591 ]), 98.64216303710529, 154, 261, 0)\n", "Optimization terminated successfully.\n", " Current function value: 98.642163\n", " Iterations: 89\n", " Function evaluations: 147\n", "------------------- Parameter summary --------------------\n", " A = -5.024, free: T, restricted: F, related: F\n", " mu = 9.88993, free: T, restricted: F, related: F\n", " sig = 10.0143, free: T, restricted: F, related: F\n", " off = 1.00059, free: T, restricted: F, related: F\n", " lin = 0, free: F, restricted: F, related: F\n", "----------------------------------------------------------\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3hU1daH353eIAkpKpAQQGroBCwUKSKoVAsQ5QoKIvApWKgXFUS9oqBYLjaKIHApoiIqCEgRRVSIhF4CSAk1hCSQnpnZ3x8nCZOQMkkmc2Ym+32eeTJnn/abk5l19ll7rbWFlBKFQqFQOC8uegtQKBQKReWiDL1CoVA4OcrQKxQKhZOjDL1CoVA4OcrQKxQKhZPjpreAwgQHB8uIiAi9ZSgUCoVDERMTc0VKGVLUOrsz9BEREezevVtvGQqFQuFQCCFOF7dOuW4UCoXCyVGGXqFQKJwcZegVCoXCyVGGXqFQKJwcZegVCoXCyXF6Q280wuzZEBwM776rLSsUCkVVwqkNfVwcREXB9OmQmAjTpkG7dlq7I2Gtm5W66SkUVROnNvQdOsC+fZCWpi2npcHevVq7tahs42mtm5Wz3PQUisrCqTtCUkq7erVt21Zaiy5dpISbX127Wuf4x45J2aqVlL6+2nF9fKSsWVPKgAApZ8+WMitLylmzpAwK0pYNhrKfIyRESheXgvpdXLR2PY6jUDgjhX/Lvr5Stm6ttTsKwG5ZjF116h798OHg51ewzc8PnnqqYJv5nXzWLHjnHcvu6oWfGNLT4fx5SE6GV16BgAB49dWK9aAjI8FkKthmMkGzZrY/jlP3eBRVGls8/etKcXcAvV7W7NEnJ2u9a/NebECA1p5H4Tu5EDd6vqXd1Yt7YijuVZ4e9JIlUvr5FTyOn5/WbsvjWNrjMRgq/hSjUNgaazz96/3dp4Qeve6GvfCrvIa+vBe5KJeGpca5KONZ2qusbiNLblYlkXddAgOl9PIq/jilXT9LXD+lubKU0VfYA0V9152hI+T0hr4i/jVLeuXFGeeijHBJr9K+ONb+Elh6XSzZzpIeT0k3TUf0eSqcj+K+6zExFetQlacjZO3fhNMb+ooMNJbWK7f0rm6J0S/pi2PN3nDeDUMI7VXadbHk+lnS47HkpilE+W5iej8WK5yDygpKKG9HyJoBERUy9MBC4DJwoJj1AvgQOA7sA9qYrRsKxOW+hpZ2LllOQ18R/1ppBjrPOFeWoSnJKJenN1z4hmHJdbHk+pXkQsr7DH5+Unp6lm7sy9qTcYaICIW+5H1H3dzKbytKoiIdIWtFAVbU0HcG2pRg6B8A1uca/DuBP3PbawAnc/8G5r4PLO185TH01hqwLI7KMjSWGOWy3vlLG3Mo6rpU5PqV5TOUtyejQkMVFaG076iHh/Z9N+/AlbVjV1RHyMtLyho1rDcOUBolGXqhrS8ZIUQE8IOU8qZgPCHEZ8A2KeXy3OWjQJe8l5TymaK2K46oqChZ1olHUlIgIgIaJ+/kRd4DwN0dHnxA+8ubb0LDhhg3beH4+E84cgQaN4aGDUAI4L33ICwMfvwRFi266fhNtn7MsaQQ+pu+ZjArtM8NeHhA377AF19ocZvLlsGaNTcLXL4c3NxgwQLYsCG/+bu1kJblzuMsA2AMc+nM9vz1EsE1qvMMnwMw9/Y5jGm/WxMtBCYEMWdCuf/gbKZMgReYw4/vHiH+ggtGXDHiSjy1mc0EAJ5kIfW8LvDSvz3w9vcET0+oVYuUTr2JiID6ydp1v0Z1RPXq/BUXiF+QJ3PmwMyZMGUKPP88uLre+GihoVr4qHnoposLBAZq4ZfJycX/37p2hS1bil9vvt22beXfX1G1Keo7WhS+vtCwIbz9NkycqIVCp6XdaF+5Eho0MNtBSrh+XfuSX7umvc/OJq7mPQwcCLWPbCIsM45qHlnUCsrikX45fLSoGjMznwdgFJ/QyuMww54w4ulq1AR+/LFmK8qBECJGShlV1DprzDBVCzhrthyf21Zce1ECRwIjAcLDw8sswN8fkpKADdfgxUM3VuTFrGdkEBcHn45MZsTpQ9wuwbQfTh6D2rXBMytL2+7qVTh0qPDhadrQwJGdEMwVmnJjva8bcIgbAeWXLt28f97NG7Qg+wMHkMDVRGiQBdl45G8azhlasE+7JkgEkiQCAe0+0r1OHPz5JwA5WSYuX5ZkG8JINGlx+s3d/6Cb4RfShQkXqZn6QzTlU78JfPIJDPnwU9i1C14109exI/69e2vXr8m/4MgRrf0acAv8Ur03043fk5YGgZOfYc1bcM+/wgiOqgt163JXw0as3RFU4CObTNCixQ0jvHQpjB4Nqak3tikqn6E4hg+H3bvLv7/C+TAaKbEDYk5kZNEdBXd37Th5N4C82PmePbW+1K2mc3RmL3XSTlNnzxn2NT9Dg25JsG6dtsOTT8LixQUPGhhIB7erJCbCVNNnPMLXkA1cAD6FKQHhvJH6vKZ1wEZN2Bo3rXfk6gofflhuQ18ixXX1zV9ABMW7bn4AOpotbwaigPHAy2btrwDjSzuXNePozbH2gG15H7lKe4z085Pys89KjwAo7vMEBZWwr8mkpetevy7llSvScOacnPvKhRuPp7/+LuXatVIuXSrl3LnyTZ835BPiy/zj/Mbd8iKhBQ6+zO1fuX55k3ybCTKaZbK5z3G55EtTvtaKhohWdH+Fc1FWV2pxv98mTXJdLKTLDvwqxzFHLuFxeat/ugQpZzIxf4cs3OUJ6spzde6ShrRM7cDff6/5d+bPl3LVKinXr5fy99/zffE1uCJDuSj9SZJepEtXcqSvj6nSwi6paNRNKYb+MyDabPkocBsQDXxW3HbFvSrL0Bc3EOLmVvyFtTQGvSyU5kO39Ljl+TzmVCSksnO7dDmg8SE5wPNHGcVf+V/qVHzyNzKG3iLlkCFS/v572S+SQlECZe20FddRWDPhV7nNtavM4EYEQbyoJV8dfFT6+Ul5O8fkXeyQt3FOumAo8ndiaUx+SVqtNQZY2Yb+QQoOxv6V214D+AdtIDYw932N0s5VWYa+pItf1IWtrAFYa428l/XzFKYiIZXVqhW9763BOVLGxmqPJI89ph1s1SrtYCdPSvn551JeuWLR51PhlIriKPdvKCZGyhdflPLXX6WUUl7f+LuMdW0tZ/GS7M1aeSvnZUCAlKdPlxyJl/c7KUtMfklarRVsUCFDDyxH8zDloPnZhwOjgFG56wUwFzgB7AeizPZ9Ci3s8jjwZGnnkpVo6EsLoyx8YSsr0sNabqCyfp7CVCSksmNHC39oRqOUOTna+/ff1zby9JTGJ4bKL8ftlkFBUr7zjpRvv13QoKtwSkVJlPQbuqmDcC1NygULpGzXTtvQ3V3KDz+0+Fwl/U4qmn9i7ZDPCvfobfmyhqEvrTdoiZGrrJjXyvA3l0drRW445drXZJIyJkYmPTZGprloFnwr90gPsm6qLVSjhgqnVBRPcb+hmJiCHQQ/H6M86dlIW2jaVMqPPpIyKalM5yrqu+7rK+XgwZYZaEu1Fn6Vp/NXpQy9Jb1BSwxVZce8WpPyaK3IDaci+4aESBkokuQ45siZTMzfvy4n8g26u3vl3GQVzol50qEbOXIoX0iBUYKUA1khO/OLnD3LVC73X1HfdRcXLXO9IgbaWmN15lQpQ2/J45QlhspeIz2KelqxV61FUdTTR2tipAEXOZfRsjrJskkTx7nJKvTFvGN3FztkLC2kBNmNn2/qhVvD/VfZwRQV6cxUKUNf2WnGeuIMvuuinj6qkSLf5QVpwEWeEzXlT6PXOMyNS6EvISFS+olU+QnPSAnyDLVlf76RYLrJBtiqpo0lVIbHoEoZekdyuZQVZygFUNIgcjv+lPtEcylBxnQcqyJtFKXSpYuUG+ghjQg5mxelL9eL7W1bo8NXmcEUFe3MlGTonW6GqT59bk4sc3PT2vXEGrMzWWu2qdKozJmk8rKYC/8Ejx2DnFbt6eQdwztMYOVfEWpOW0XJSMnw4TDb+1V6sInxvEsafvj5wahRls0uV1asZV+K+h0kJWntlUJxdwC9XpUVXqkn1nK52OJpRS/3UHFPK9HVf5By+/bKPbnCsTAapXzpJSmnTCm2Z1xULLyzu/+oSj16W1DWHq+15qO0xdOKXnNnFv20IpnKG3DffbB2beUKUNglN/3W0jIhOlpbuH4d/+qyyJ5xeLiNe8x2jjL0ZSQuDqKiYPp0yyf9tpbLxRaPe7ZyDxWm6IncBYdmfg/Nm8NDD2lVQhVVhsK/tVmvXif21p6wahXMmqUVABNCb5kOgTL0ZaQ8Pd6ijZh9Vl/US2txTyv3PRaslcHs1k0T8emnlTqGoLAfCv7WJMvS+9MydQfPVPsfjB+vjHwZUIa+jJSnx2uvA8RFoZfWEp9W/Pzghx+gd29SNu8iqq0s0xOVwjEp+FsTfMoohrCUuKhoPWU5JJVQ+Ni5KU9t9Px6+Q6A3Wr18IDVq2lY250rVwWYjIBrgSeqy5f1FqmwJsOHQ+yuHBqn7eYP7mI1j+LnB5/Y4ZOwvaN69GXEkXrnToenJ02buRBmOkUsrbiL3wHbjCEobE+fPjAr53m205m6nATUb628KENfRmwe/6oowPDhIHx98SKTb3iImpyz2/EORfnIG4N5o9bHjMj+GNeXXuCkrKd+axVAGXqFQ9GnDyS7h9CfNfiSxjc8hK9rpurlOTDmg+uTJ0PbtrDtlc28lTaWda69ab/5LTUGU0EsmhzclpRncnBFFeXbb7Wwy6eegvnzVRSGAxIXBwMH3piIG6AW8eyjBeepyd38TppLdYKC1BhMaZQ0Objq0SsclwED4JVXtAnZ86yEwqEoHK4McJ6avMlU+rKW61RXYzBWQBl6hUNjfGU6c/ptIzjCT8XUOyCFw5W9SUfiwnu8xD/UA+w358SRUIZe4bDExUFUexdeecMTkZiAmDyJemE5BAaqRCpHwTxB7x628Q91aUNMgW1UpE3FUYZe4bCYP/Z3YAcvGt5hxIUZJCerRCpHIS9cOZCrLGUIyQRwyb8Ryckqqs2aOJWhj78Wz9eHvuZa1jUALqVeYvvp7aTnpAOQnJlMXGIcOcYcAAwmA0aT6vY5KuaP/d/Rn4U8yb/5D3ezw2bF2BQVIy9c+erQF6ntdolGMcuJT/ZTht3KOJWh//3s7zzy1SPEX4sHYNPJTdyz6B7OXz8PwDeHv6HhfxtyIfUCAPP/no/b625cuK4tL4pdRJO5TUjOTAbgh2M/MGzNMDJyMgCIvRjL//b/D4PJAEBGTkb+e4XtKVyXZxwfcIZw5jMCD7LUIJ4dYFFdop9/hsWLYcIEaNPG5hqrAk5l6HvW78neUXupF6gN4txb7142P7GZWtVqAXBPnXtYMmAJwT7BALSr2Y7XurxGgFcAACE+IbS4pQXebt4AnE05y9ZTW/Fw9QBg9aHVPPHtE7gI7bLN+GUG1d6qRl6I6rJ9y5iwcUK+nivpV/JvEgrrUzhLOZVqjOFjmnCEF5ijBvF0xuJKrz/+CA0aaBFUikpBxdGXgWtZ17iUeokGQQ0A2PLPFnaf383EDhMBmLhpIj8d/4l9o/cBMOSbIfx65ldOP38a0G4UAsHDTR/W5wM4MSkpEBEByckwhrms5hGyA27h1Cnl39WL0FDNwJtH1bi4UHRM/NWrUKOGTfU5GyXF0Vtk6IUQvYAPAFdgvpRyZqH1dYCFQAhwFRgipYzPXWcE9uduekZK2bekc9mzoS8rm05s4lLaJYa0GAJA5y864yJc2DZsGwDjN44nrHoY4+4cp6NKJ8Zk0pKoVCKVLnTtCtu2Fd2+ZQtw+LBm+Rs1srU0p6RCCVNCCFdgLnA/0BSIFkI0LbTZbOBLKWULYAbwltm6DCllq9xXiUbe2ehRv0e+kQfYMnQLqx5dlb98KOEQJ5NO5i/3Wd6HhXsW2lSj03L5MnTsCCtW6K2kylLi3AZSwujRcO+9YFDjXJWNJT769sBxKeVJKWU2sALoV2ibpsCW3Pdbi1ivANxc3Aj1Dc1fXvf4Oj64/wMA0nPSyTRkkm3MBrSB3odXPcyOMzt00erwBAVBdjZMnAjp6XqrqZKUWOn122/hl19g6tSbN6pCLNyzkA///LDSz2OJoa8FnDVbjs9tM2cv8FDu+wFANSFEUO6ylxBitxDiDyFE/wqpdWJ83H3Y9K9NjIoaBcDJpJPsPr+bDIM2mBt/LZ4VB1aQZcjSU6bj4OoK770H8fGY3p2jZqTSgcKVXg0Gza43qZdF8tPjkc2awYgResu0KanZqayPW5+/vOHEBr498m3ln7i4WcPzXsAjaH75vOV/Af8ttE1N4BtgD5ovPx4IyF1XK/dvPeAUUL+Ic4wEdgO7w8PDK22WdEfDZDJJo8kopZTy/Z3vS6Yjjycel1JKaTAa9JTmMFy/t79MdfGTdb0vSJDS11fK1q2lPHZMb2VVi2PHpGzVSrv+E3hbSpDP1N9U5f4PkzdNlq6vucqEtAQppZRp2WlWOzawWxZnx4tbIW8Y4buADWbLU4ApJWzvB8QXs24R8EhJ52vbtq3VPrgzYTQZ5Z/xf+Yvj1w7UvZf0V+aTCYdVdk/dwQeldm4ybeZkD+DgIuLlCEheiurWoSEaNcdpJzBy/IrHq4S/4cjCUfk3Qvuln+f/1tKKeWZ5DNyx5kdlfK7LcnQW+K62QU0EELUFUJ4AIOBteYbCCGChRB5x5qCFoGDECJQCOGZtw3QAThkwTkVhXARLrSv1T5/uVFwIyJDIhG5ESV/xP+BSZqK273K4t2yIfeznmm8lt+mEqlsj3kW86u8zkBWYTJprh1L3WmOMil8jjEnPwnzFr9bSMtOIyE9AYAw/zDuDrs7/3drM4q7A8iCPfEHgGPACWBqbtsMoK+84d6Jy91mPuCZ2343Wmjl3ty/w0s7l+rRl53DCYelmC7k7B2z9ZZidyxZIqWfX25PHoMEbXnJEr2VVS2WLJGysc9p2ZltsuD8bJa508xdP/bsgjOZTPLuBXfLbou72fzclNCjVwlTToDRZGTFgRX0vL0nwT7BHE44TFpOGlE1iwyprVLkJVKFJ+/lWwYwmBXEBbRXiVQ2JiUF1oYO59HsZYRxliuEFFhfbCJVLmVKvtKBv879Rbua7RBCsPLASnw9fHmwwYM27bmriUecHFcXVx5v8Xh+aYdp26Zx/7L7VfkFbkR+7L1Wj3pB1/ir1zRVDVEH/C8d41/GxXg9P5pmXUJuWl+aO61w3XpL9rEVPx77kTvm38EPx34AYFCzQfRu2Nv27pkSUIbeCZnfdz7fDf4Ob3dvpJRsPLERe3tysznVqmlFs376CXbu1FtN1WP6dPD0hMmTS06kKoby7FOZGEwGTiWfAqDn7T355MFP6FG/hz5iLEAZeiekumd17g67G9AqePZc2pOVB1fqrMoO+L//g5AQrbqWwnbs369lKI8bB7fcUnIiVTGUZ5/KJPrraLp/2Z1MQyZuLm6MihqFl5uXPmIsQBl6J+feevfyv4f+xyNNHwEgKSNJZ0W2ocgIDT8/mDQJNm2CP//UW2LV4eRJqF8fxo8Hbk6ksmRykfLsY22klPmRbS/c+QJvdX8LT1dP2wmoAGowtgqRmp1Kq09b8XCTh3m7x9t6y6k04uJg4EDtb1oa+PpCw4awciU0qJWuuW/699dG8xS2wWRy6Oudlp3Go189SteIrkzoMKH0HXRADcYqAPBy82JIiyH0bthbbymVivkUg0DB2aZ8fOChhxza6DgUf/2lPU45+PX2cfch1DeUap7V9JZSLhz76ivKhJuLG9O7TKdTnU4AfLb7M7478p3OqqyPRREas2fD00/bVFeVIz5eqyA6fbreSsrNl3u/5FLqJYQQLOq/KL8WlaOhDH0VxWgysnjvYhbGLnS6iJySIjTyfPfvvZqEXLAA46Gj+oisCrz7ruZMd9DCZaeSTzHqh1F88OcHekupMMpHX4XJNGRiNBnx9fDletZ1vNy8cHd111tWhTGfbSqPgADYvFm7CcTFgW/aJU4RwfqgITTfOY8GDXST65wkJGj/hEce0eaDdSBM0pQ/Xejei3tpFtoMVxdXnVWVjvLRK4rEy80LXw9fTNLEgJUDeHjVw07Ruy8uQqNXrxu++8vcwhc8yYOJX9Kl4Xm7rp3ikHz0EWRkwOTJeispE5fTLnP3grvzk59a3trSIYx8aShDr8BFuBDdLJpHmj5iV9l81qaw734243HDwPPMKX7iakXZkVJ7fOrTB5o00VtNmfB288bXw1dvGVZHuW4UN/HXub+ICIgoMBuWM7B0qTZ7XWrqjbaxfMBftOcP7rKr2imOhtEIc+bAzJkwZQo8P9aE6/Vkh5nwO/5aPLf43oK7q7tWBMwBOzzKdaOwmIycDPqt6MfI70fqLcXqFJVd+SHj+IO7APupneJoxMVBVJQWXJOYKJn5ajrt7nAhLtExjHxyZjJ3zL+DcT+NA3BII18aytArCuDt7s1Xj37Ffx/4r95SrI65737JkhuROfU5zhtMpbqvUbfaKY6Med5CTzZwJD0MYmO1vAUHIMArgCkdpzCm3Ri9pVQaytArbqJjeEdqV6+NlJKP/vzIKcsmmPfu2/A3U/kPPbJ/ZNw4+57Uwh4xH/t4kffIwpMDsmmZJhXRg7MpZ4lL1AZlnm3/LM1CnfdxThl6RbEcvnKY8ZvGs3DPQr2lWB3z3v2bhx7ignsYzxnf5+pV1MBsGcnLW4jkAPexiY94jhw8MBjs91pKKYn+Opo+y/tgNNnpnciKqMFYRYnsv7SfZqHNnNJvmUdoKDx55R3elpNoSSz7aKkGZstAXt7CrOQRRLOcMM6SxA3/vL1ey8MJh0nOTOausLv0lmIV1GCsotw0v6U5QgguXL/AG9vfcIo4+8JERsI8OYI0fBiHlgVZ1vlMqyJ5Wcb168PrLyUz3GsZW24bUsDIg30NcucYc/jx2I8ANAlp4jRGvjSUoVdYxIoDK5j520yOJR7TW4rVGT4ccvxq8DFjSOVG7QR7dj3oTcFIG5j8lj9Phm/h4hMT7WqCkMJ8svsTei/vTcz5GL2l2BTlulFYhJSSU8mnqBtYV28pVqeokgnm2KvrQU+Km8M1MFDr6RcuP2Evc/TmGHNYf3w9fRv11VuK1VGuG0WFEULkG/lvD3/LzrPOMx2f+cBsl3sk7fgLF274a+zJ9WAvmEfa9GAjn/IM1U1JtGih/wQhRbHz7E5Ss1Nxd3V3SiNfGsrQK8pEliGLiT9P5O0dzjlxyStR6/mLO+jJhvw2e3I92AvmFUKf53368D34+tnldUrJTOH+Zffz7Lpn9ZaiG8p1oygzp5JPcavfrXY9R2Z5SUnIJuOWOsTINvRGG7SzJ9eDvZDn7gpKPs5xGjCdaXwQMN1ur9Pmk5tpHNyYWtVr6S2l0lCuG4VViQiIwMvNi0xDJkv2LtFbjlXxD/Hg1ldG8qBYjzxx0m5cD/ZGnrvr+IufgJsb08+NtLvrZJImDiccBqB7ve5ObeRLwyJDL4ToJYQ4KoQ4LoS4qe6oEKKOEGKzEGKfEGKbEKK22bqhQoi43NdQa4pX6Mv8v+fzxJonnC+C4emntZHFTz/VW4l9k5EBX3wBAwZAzZp6q7mJD/74gFafteJQwiG9peiOW2kbCCFcgblADyAe2CWEWCulNL96s4EvpZSLhRDdgLeAfwkhagDTgChAAjG5+zpfTn0VZHTUaFre0pK2NdvqLcW61K4N/frBt9/C22+DEyeLVYj0dIiOhkGD9FZSJE+0fAKAJsGOVSq5MijVRy+EuAuYLqXsmbs8BUBK+ZbZNgeBXlLKs0JLoUyRUlYXQkQDXaSUz+Ru9xmwTUq5vLjzKR+9Y3I44TDh/uHOU8v7/HktVtDbW28lijJyPes6fh5+Tp3NXRQV9dHXAs6aLcfntpmzF3go9/0AoJoQIsjCfRFCjBRC7BZC7E5ISLBAksKeSEhLoP389kz+2bFmEyqRmjU1I29nwQp6k5cN2zbwJMuf3YHRYF/XxyRNPLTqIQZ/Pdgps7jLi7UGY8cD9wgh9gD3AOcAixPHpZSfSymjpJRRISEhVpKksBUhviF8+uCn/LvTv/WWYl127UI2acLCCYcJDlblEMyzYYclz2HA3O50b5NkV1nDAkHfhn3pWb9nlevRl4Qlhv4cEGa2XDu3LR8p5Xkp5UNSytbA1Ny2ZEv2VTgHj7d4nNuq3YaUkvScdL3lWIWTxjoYjp0k/f3PSUxU5RDy6s7LtDSe4EtW8wi/HqxhV3XnhRA8d8dzPNXaDgP6dcQSQ78LaCCEqCuE8AAGA2vNNxBCBAsh8o41Bcira7sBuE8IESiECATuy21TOCkj1o6g34p+mKSp9I3tnDv7hvKNfIjHDIvxIoO0NNi7F7sybLYkLxt2ECvx5xqfMspusoYzDZncv+x+tp3aprcUu6RUQy+lNADPohnow8AqKeVBIcQMIUReLnEX4KgQ4hhwC/Bm7r5XgdfRbha7gBm5bQonpUN4B7pFdHMK/2hkJHzGSGqQxMN8DVTtcgh52bDP8BkHacoOOthN1vCF6xc4lXyKjJwMvaXYJSozVqEohqVLYfQoSUxaIy5xC535FQ8P8PDQ/NTPPw+urnqrtB0pKRAVfpmd15ryOq/wIePsKmvYYDLg5lJqxLjTojJjFTZnx5kdDF0z1KFn7+nTB9zcBZOZyQeMAyTZ2ZCaWjX99f7+EJcSSnDmOT5IH2kXWcPpOel88McH5BhzqrSRLw1l6BWVwomkE2w/vZ1z1x137D0vzf8b+RDbQx7BxeVGFEeV9NcbDJrvytPTbvILvjn8Dc9veJ6/zv2ltxS7Rhl6RaXwrxb/4uCYg4T7h+stxSp0vP0iE01v4UlmfluV89cvWwa33w7n7OfmPaTFEGKfiaVDeFW645YdZegVlYIQAh93H4wmI8v3L3f4wdkxnQ/wFv9mAN/mt9nLQKTN+PxzcHe3i7o2WYYszl8/D0DLW1vqrMb+UYZeUamsObKGx755jPXH1+stpUK0m9SNUy51eZp5+W1ubpofv0pw8CD8/rtW8M0OEpFm/T6LJnObcO6a/Txd2DNq9EJRqQxoMoD1j6+nZ/2eekupEP6BLi/GdAIAACAASURBVPi/PoKIqVORx+KgQQO9JdmW+fO13vxQ+yhAO7jZYNxd3Kt06eGyoHr0ikrFRbjQ6/ZeCCFIy07TW07FGDZMi6ecP19vJbYlMxO+/BL69wc7KVFye43bmdRxkt4yHAZl6BU2IeZ8DHU/qOvYmYs1a2rli69WsZw/NzdYuBAm6W9Yl+1bxsjvRzpNmQ1boQy9wiY0CWlC93rdCfGxjx5hufnqK5g3r/TtHJy8KpXBwfDuB24Ye/eDtvrPO3A65TSHrxx2ymksKxOVGatQlIeLF+HWW/VWUSnExcHAgdrfW9OOM8JtMcuCxhKfFcLLL+ufEVzVM2CLQ2XGKuyGa1nXmPzzZMeOlli0CGrV0nL/nZC8KpVpafA08xhveIsrlwwkJ+uXEXw44TB7L+4FUEa+HChDr7ApCWkJfPjnh/x0/Ce9pZSfbt20vwsW6KujksirUulONsNYxA/05iK3AfplBE/ZPIWeS3uSZciy7YmdBGXoFTalfo36nBx3kuFthustpfyEh0OvXtoApcGgtxqrk1elsg/fcwuX+ZyRBdbrkRG8sN9CVg9cjaebp21P7CQoQ6+wObf6ab7tk0knHbfo2dNPa/PKrluntxKr06ePFmjzNPM4QxgbKJgDYcuMYIPJgJSSGt416Bje0TYndUKUoVfowoHLB2gytwkL9jio++PBB+G225zSfePvD0lXjPQa6E/oy89QPaDgyKstM4Jn/DKDnkt7km3Mts0JnRQ1qqHQhciQSKbdM43eDXvrLaV8uLvDqlXQqJHeSioHV1dYuRIvIOl1/WTUqlaLBjUa4OHqoZ8IJ0CFVyoUioLk5MA//0DDhnorUZQBFV6psFvir8UT/XU0p5JP6S2lfGzZAgMGONeg7HffaU8qO3boJmHPhT38cOwHh696ai8oQ6/Qnc0nN7Pnwh69ZZSP69dhzRr44Qe9lVSYvGzYbY99xrWAcIzt7tRNy5w/5jBi7QhV6sBKKNeNQncycjLwdrePGYvKjMEAERHQvDmsd9xSzHnZsMajx9mX0YDX3WfwbbNXWLlSn0KdOcYcjiUeIzI00vYnd1CU60Zh1+QZ+diLsY73qO7mpoVabtgAJ0/qrabc5GXDPp4xDwOufJIzXJfEKIPJgMFkwN3VXRl5K6IMvcIuWBe3jtaftXaYCUrMi359ZhiOFMJhip0VKFj2rrYcGQnSZGIQK/mePlygpi6JUfNi5tHs42ZcTrts2xM7OSq8UmEX9KjXg/d7vs89de7RW0qpmBf9SkuDl+bUxjv4Re4LisTey5wV1j5tmjYV7GOPwe7dLrRJ/ZsAkgF9pkqsF1iPjuEdHb/KqZ1hkY9eCNEL+ABwBeZLKWcWWh8OLAYCcreZLKVcJ4SIAA4DR3M3/UNKOaqkcykfvcLeCQ2FxEStFEAeLi4QFASX7bwjWpz2wECtZ5+cfKM9IECr2+bvb3OZinJQIR+9EMIVmAvcDzQFooUQTQtt9jKwSkrZGhgMfGy27oSUslXuq0Qjr1DsOreL/iv623W0RV7RL3NMJohqnKolUdkxxWl/uG4MSQ3aIw8dRkqQEpKSbGfkL6VeYu5fc1UGbCVhiY++PXBcSnlSSpkNrAD6FdpGAtVz3/sD560nUVGVyDBkEHMhhuNXj+stpVjyin6Z4+cHr9VZCIMGQUyMPsIsoDjtE3zmwqFD2ixaOrD8wHLG/TSO08mndTm/s1Oq60YI8QjQS0o5Inf5X8AdUspnzba5DdgIBAK+wL1Syphc181B4BhwDXhZSvlrEecYCVqJvPDw8LanT6t/dlUm25ht1ynvKSlaROVNbo69Kfg3rQWPPKLVrLdDitJet3oiJ7JrI4YNg08+0UWXlJLDVw7TNKSws0BhKbYIr4wGFkkpawMPAEuEEC7ABSA816XzIvA/IUT1wjtLKT+XUkZJKaNC7GTyYYV+eLh6YJImDl4+qLeUIvH319waeS6OfDdHuD888QSsWAEJCXrLLJKitJ98eSEiMxP+7/900ZRtzEYIoYx8JWKJoT8HhJkt185tM2c4sApASrkT8AKCpZRZUsrE3PYY4ASgCmgoSmXSpkncueBOEtLs02AWy7PPQlYWvz4xr0D4ot1iNMLHH0PnzraPpUSbOSpsTphjTxrvAFhi6HcBDYQQdYUQHmiDrWsLbXMG6A4ghGiCZugThBAhuYO5CCHqAQ0Ax80qUdiMEW1G8OmDn1LDu4beUiwiLzY9sENTfvXozuVNe0lM1G/qPYsxmWDqVO2lA0IIOoR1oFmo7W8yVQlLwysfAN5HC51cKKV8UwgxA9gtpVybG4UzD/BDG5idKKXcKIR4GJgB5AAmYJqU8vuSzqXCKxWORuHYdB/SSMc3f72jhF4qHJuSfPSq1o3Crlmydwknk04yrcs0vaUUS1Gx6QDBJHCFYEDQtatW6NKuiInRKlQ+/TR4277W0MI9C+nXqB9BPkE2P7czomrdKByWP+L/YMOJDRhM9lsGuKjY9M78Qjy16cx2XTJMLWLmTM23pMMgwvGrxxmxdoTjzjDmYKgevcKuycjJwMvNCyGE3lKKZelSGD0aUlNvtHmRwWnqsIt2DAn40f4yTE+e1MpSTpwIb72li4SDlw9SJ6AOfh5+pW+sKBXVo1c4LN7u3gghSM9Jt9sInLzJtM3xCvCm+r+f40HWkfTrAfsy8gDvvw+urhjHPHdTgbPKxiS1x5/I0Ehl5G2EMvQKu8dgMtDikxa8uPFFvaUUSXFx9V4vjgEfHy0cx564ehUWLOBa78eI6luT6dOxWYSQlJL7l93Pa9teq7yTKG5CGXqF3ePm4saUjlMY2Wak3lLKRlCQVnNgxQrN8tsLCQnQujW9t41n3z4tUgi0v5Vdgz7bmE2EfwShvqGVdxLFTSgfvUJRiRjjL/DFf9OYPP92pkyB558HV1e9VWl07QrbthXdbncRQopSUT56hVNwPes6M36ZweGEw3pLsYi4OIjqcxvP//d2EhPhtVeNNkueKmpyEUCbCevSJaD4AmeVFSG09+Jejl45WvqGCqujDL3CYcgyZjHr91n8dPwnvaVYRN70fGlpkhUM4p30/7PJ9HxxcRAVxU2+9xO7rmqZXS9qYx1FDSK7uWntlcGETRPoubQnRpM914RwTtQMUwqHIdgnmBNjTziMfzcyMs81IrhCMCP5nHdME4loVq9Sz9uhQ8EErjzf+5ou7/NS+jWYPBm4MYhsK5Y+tJS4xDhcXezEd1WFUD16hUORZ+SvZ13XWUnpmLtG3mQqBtyY5vI6u3ZVbihjUQlc/qarPJP1gVZCuXnzyjlxKYT6htIh3MazjSsAZegVDsi6uHXc9u5tdlvGOA9z18gFavIJoxli+pKaqUcrNZSxKN/76+4z8DGlan4cG/PzyZ/pt6IfF1Mv2vzcCg1l6BUOxx217iC6WbTdJ9uYx9eHhMA7YjIZeDOVN0lLgz17oFEj6/fub/a9S2qQRM6TI3UpRXwx9SKnkk8R6BVo83MrNFR4pUJhA/JCGXuwkd1EkcSN8su+vtCwIaxcqVUlyMNohDlztJI0VgnNNBp1i+2UUtp1GQtnQIVXKpySf5L+YeGehXrLsIg8d8om7iOJGrhgxB1tIuyiEpWKi5wps6tnxw7Yv197r4ORjzmvzZ+rjLy+KEOvcFg+3vUxz61/jqsZV/WWUirm7hQf0viL9kxmZv56k6mgV+VGaKa2XK6s1awsGDYMHntM8x/ZmO2ntxM1L4qVB1ba/NyKgihDr3BYJneczLFnjznELFTm/vrPlvjyj1sDpvAWdXMnXCucqFRU5Ezhm0GpvPEGHD8Os2aBDj3qdjXb8fEDH9O3UV+bn1tREGXoFQ5LkE8QtarXAm5URHQE+vSBV33exYAbCxiOC8abEpUqnLX622/wn/9oPfpevawlvUx4u3szut1ovN1tP6mJoiDK0CscGiklQ9cMZeT3jlPwzN8fDqXUotqi/9KVbeS8PIOpU6F+fa3z/c47MHYsGArNtWJx1mpKCgwZAhER8OGHlfERSuWFn15gw/ENupxbcTPK0CscGiEEYdXDqF29tt5Sys7QoVwbMJRj767lP9OySEyESZO0CJukJG3stHVrOHbsRunjwnXti6xp4+2t+eWXLoVq1Wz+sZIykvj+2Pfsv7zf5udWFI0Kr1QodKROSDpXEgXpsmj3RkkTi980KbkP1KyeypVMP15+Wd9KmQaTAaPJiKebpz4CqiAqvFJRJfgj/g/OXTunt4wyUa+ZD+nSG2/S+YCxBFIwgqikAdjCkTl905ez/WIDQpMrN/O2JK6kX8FoMuLm4qaMvB2hDL3CKbicdpnOX3TmvZ3v6S2lTOQNurZmD8/wGZvpTjA3pkwsaQDWPDKnD2tZwr84RkPOEmaTSUSKYvja4XT8oiP25imo6ihDr3AKQn1D+T76e17r6lhT1OXF1/9OB/qylsYcYRtduAWtLkxJA7DDh4Ofr2QQK/iKR4mhLb35gQx8gHKEY1qBYS2HMbLNSJUgZWdYZOiFEL2EEEeFEMeFEJOLWB8uhNgqhNgjhNgnhHjAbN2U3P2OCiF6WlO8QmFOz9t72n39m8KYx9dvkD3x3rKOSJ9TXAxtifx5c4kDsGPHwqDsL1lBNHtozf2sJ5Ubg6/mTwPFTkRiZQY0GcCTrZ+snIMryo+UssQX4AqcAOoBHsBeoGmhbT4HRue+bwqcMnu/F/AE6uYex7Wk87Vt21YqFOVlz4U9suuirvLC9Qt6Syk/sbFS3nuvlKdOacv790t5+LCU2dkybn+GHFF/i+zvuU6ClEE+6fI/tefKmD9zZECA+fTkUgYESJmcLOWxY1K2aiWlr6/W7usrZevWWru1OHDpgPzvn/+VGTkZ1juookwAu2UxdtWSHn174LiU8qSUMhtYAfQrfL8Aque+9wfO577vB6yQUmZJKf8BjuceT6GoFHzdfTmVfIoTV0/oLaX8tGwJmzZBnTra8uTJ0KQJ+PoS1tyfeSe68WrWvwFITPfm5fNj6NXbLf/JIO+V9zRglXIKpbDy4EqmbJ5Cek669Q6qsBqWzDBVCzhrthwP3FFom+nARiHEc4AvcK/Zvn8U2rdW4RMIIUYCIwHCw8Mt0a1QFEmDoAYcH3scF+FEw09vvQWDBsGhQ3y7zMDSs/fwK53yV5fmi78x09UNrO2/f63LazzV+imHKEdRFbHWVILRwCIp5btCiLuAJUIIi79GUsrP0dw/REVFqeF6RYVwES6YpInDCYeJDI3UW07Fad48f1YoQyT8MhpSU2+sLq00wvDhsHt32fYpC9nGbDxcPYgIiLDOARVWx5JuzzkgzGy5dm6bOcOBVQBSyp2AFxBs4b4KhdWZtGkSd8y/g8T0RL2lWJXyTOhdmZOAH71ylLA5YWz5Z0vFD6aoNCwx9LuABkKIukIID2AwsLbQNmeA7gBCiCZohj4hd7vBQghPIURdoAHwl7XEKxTFMbzNcOb1mUeAV4DeUqxCXtRM/frw8staHZzCvvjiMI/ssXQfS5FIOoR1oFmo7WeuUliORSUQcsMl30eLwFkopXxTCDEDbZR3rRCiKTAP8EMbmJ0opdyYu+9U4CnAADwvpVxf0rlUCQSFoiCFSx0UNyOVompTUgkEVetG4dR8secL0nLSeLb9s3pLKTehodosU+b16UuqgWMrvjr4FV3rdiXYJ1g/EYp8VK0bRZVl3fF1rDmyxqFT8q0yCYmVSUhL4PFvHuft397WT4TCYqwVdaNQ2CUL+y7Ez8PPoVPyKztqpjyE+Iawd9RegnyC9BOhsBjVo1c4NdU8qyGEICMng9Ts1NJ3sEMqM2qmPOQ9HTUJaUKob6g+IhRlQhl6hdOTkplC/Q/rM/v32XpLKReVGTVTHsb8OIaXNrykz8kV5UK5bhROj7+XP6OjRtOtbje9pTg8UkrcXNycK/O4CqCibhQKRZmRUjr0uIczoqJuFAo0F86MX2ZwJf2K3lIcktPJpzl65SiAMvIOhjL0iirDuevneO2X11gXt05vKQ7J9F+m025eO4cd1K7KKB+9osrQNKQpJ8eepE5AHb2lOCRv3/s2jzZ91OEmd1GoHr2iipFn5FXd9LIT6hvKAw0eKH1Dhd2hDL2iyrHq4CrC5oRx7poqpGoJsRdj6bO8D2dSzugtRVFOlKFXVDmiakbRp2EfNaBoISeunuDg5YNU86hW+sYKu0SFVyoUilIxmoy4urjqLUNRAiq8UqEoglPJp5j/93y9ZdgtUkp2ndsFoIy8g6MMvaLK8smuTxj30zgVV18MW09tpf389qw+tFpvKYoKolw3iipLUkYSaTlp1K5eW28pdkmmIZPFsYsZ1moYnm6eestRlEJJrhsVR6+osgR6BxLoHQhAliFLGbNCeLl58UzUM3rLUFgB5bpRVHnGrR/HvUvudejJSaxJtjGbfiv6sf30dr2lKKyE6tErqjxtbmuDv5c/BpMBd1d3veXozunk0xxKOERadpreUhRWQvnoFQrFTeQYc3BzcVO5Bg6ECq9UKCxg59mdfH/0e71l6ErM+Zj8Jxtl5J0HZegVCrSY8Yk/T2T6L9OrrK/+UuolOn3RiUmbJuktRWFllI9eoUCrr/5l/y8J8Q2psj3ZUN9Qlj60lDa3tdFbisLKWGTohRC9gA8AV2C+lHJmofVzgK65iz5AqJQyIHedEdifu+6MlLKvNYQrFNambmBdQOvdp2anUs2zatV2EULwUJOH9JahqARKdd0IIVyBucD9QFMgWgjR1HwbKeULUspWUspWwEfAN2arM/LWKSOvsHeklNy39D6GfTdMbyk2wyRN9Fneh//t/5/eUhSVhCU9+vbAcSnlSQAhxAqgH3ComO2jgWnWkadQ2BYhBP0b9cfH3afKzIualJHEtaxrGEwGvaUoKglLDH0t4KzZcjxwR1EbCiHqAHWBLWbNXkKI3YABmCmlXFPEfiOBkQDh4eE3HTcnJ4f4+HgyMzMtkKtQVIxu1boBcOTIEQC8vLyoXbs27u7OGWMf5BPEtqHb9JahqESsPRg7GFgtpTSatdWRUp4TQtQDtggh9kspT5jvJKX8HPgctDj6wgeNj4+nWrVqREREVIkelkJ/pJQkZSZhNBkRGYL4+Hjq1q2rtyyrs+bIGrpGdMXfy19vKYpKxJLwynNAmNly7dy2ohgMLDdvkFKey/17EtgGtC6ryMzMTIKCgpSRV9iUK+lXuJpxlRo1ajjl0+SF6xcY+NVAXt/+ut5SFJWMJT36XUADIURdNAM/GHis8EZCiMZAILDTrC0QSJdSZgkhgoEOwDvlEaqMvMKWCCGoF1gPV+HqtN+926rdxo6ndlC/Rn29pSgqmVINvZTSIIR4FtiAFl65UEp5UAgxA9gtpVybu+lgYIUsmG3SBPhMCGFCe3qYKaUsbhBXobAr3Fy0n4fJZCLHmKOzGuuSnpOOj7sP7Wq101uKwgZYlBkrpVwnpWwopawvpXwzt+1VMyOPlHK6lHJyof1+l1I2l1K2zP27wLrybYcQgiFDhuQvGwwGQkJC6N27d5mOExERwZUrJU90Udw2X331FU2aNKFr165F7FV+3n//fdLT08u836JFizh//nz+8ogRIzh0SJ/7uJ+fX6Ud+3jScS6nXXaaqJTLaZep/2F9vtjzhd5SFDZClUCwEF9fXw4cOEBGRgYAmzZtolatWjbVsGDBAubNm8fWrVst2t5gsMwwlcfQG43Gmwz9/Pnzadq0aQl7OSY1/WoS5BOU38N3dNxd3OlZvyd31C4yeE7hhDikoe+yqAuLYhcBWpW9Lou6sHTfUkB7JO2yqAsrD6wEICUzhS6LuvDNYS2H60r6Fbos6pJfvOpi6kWLz/vAAw/w448/ArB8+XKio6Pz1129epX+/fvTokUL7rzzTvbt2wdAYmIi9913H5GRkYwYMaJAHZWlS5fSvn17WrVqxTPPPIPRaKQ4ZsyYwW+//cbw4cOZMGECmZmZPPnkkzRv3pzWrVvnG/9FixbRt29funXrRvfu3QscIy0tjQcffJCWLVvSrFkzVq5cyYcffsj58+fp2rVr/pPC6NGjiYqKIjIykmnTbqREREREMGnSJNq0acPy5cvZvXs3jz/+OK1atSIjI4MuXbqQV3nUz8+PqVOn0rJlS+68804uXboEwIkTJ7jzzjtp3rw5L7/8cpE98cmTJzN37tz85enTpzN79mxSU1Pp3r07bdq0oXnz5nz33Xc37btt27YCT1nPPvssixYtAiAmJoZ77rmHtm3b0rNnTy5cuADAhx9+SNOmTWnRogWDBw++6Zh+nn54uXkBOEWvPtA7kEX9F9E0xPluyoqicUhDrxeDBw9mxYoVZGZmsm/fPu6440aPaNq0abRu3Zp9+/bxn//8hyeeeAKA1157jY4dO3Lw4EEGDBjAmTNnADh8+DArV65kx44dxMbG4urqyrJly4o996uvvkpUVBTLli1j1qxZzJ07FyEE+/fvZ/ny5QwdOjQ/MuTvv/9m9erV/PLLLwWO8dNPP1GzZk327t3LgQMH6NWrF2PHjqVmzZps3bo1/2bx5ptvsnv3bvbt28cvv/ySf9MCCAoK4u+//2bIkCH5emJjY/H29i5wrrS0NO6880727t1L586dmTdvHgDjxo1j3Lhx7N+/n9q1i57Cb9CgQaxatSp/edWqVQwaNAgvLy++/fZb/v77b7Zu3cpLL71kcQGynJwcnnvuOVavXk1MTAxPPfUUU6dOBWDmzJns2bOHffv28emnnxZ7jC/2fEGbz9qQnlN2N5c9cCn1EtFfRxN/LV5vKQob45DPotuGbct/7+7qXmDZx92nwLK/l3+B5WCf4ALLt/rdavF5W7RowalTp1i+fDkPPPBAgXW//fYbX3/9NQDdunUjMTGRa9eusX37dr75RnuaePDBBwkM1Kau27x5MzExMbRrpw2GZWRkEBoaarGW3377jeeeew6Axo0bU6dOHY4dOwZAjx49qFGjxk37NG/enJdeeolJkybRu3dvOnXqVOSxV61axeeff47BYODChQscOnSIFi1aAJoRtgQPD4/8nnXbtm3ZtGkTADt37mTNGi1n7rHHHmP8+PE37du6dWsuX77M+fPnSUhIIDAwkLCwMHJycvj3v//N9u3bcXFx4dy5c1y6dIlbby39f3j06FEOHDhAjx49AM31dNtttwHa//Xxxx+nf//+9O/fv9hj1A2sS/0a9cnIycDH3cei62BP7Lm4h40nNvJK51f0lqKwMQ5p6PWkb9++jB8/nm3btpGYmFju40gpGTp0KG+99ZYV1Wn4+voW2d6wYUP+/vtv1q1bx8svv0z37t159dVXC2zzzz//MHv2bHbt2kVgYCDDhg0rEENe3LEL4+5+o565q6urxeMFeTz66KOsXr2aixcv5t9cli1bRkJCAjExMbi7uxMREXFTfLubmxsmkyl/OW+9lJLIyEh27txJYX788Ue2b9/O999/z5tvvsn+/ftxc7v5p9ElogtdIrqU6XPYE71u78WZ58/g62HZ/1DhPCjXTRl56qmnmDZtGs2bNy/Q3qlTp3zXy7Zt2wgODqZ69ep07tyZ//1PKxa1fv16kpKSAOjevTurV6/m8uXLgObjP336tMU6zM937Ngxzpw5Q6NGjUrc5/z58/j4+DBkyBAmTJjA33//DUC1atW4fv06ANeuXcPX1xd/f38uXbrE+vXriz2e+X6Wcuedd+Y/+axYsaLY7QYNGsSKFStYvXo1jz76KAApKSmEhobi7u7O1q1bi7xederU4dChQ2RlZZGcnMzmzZsBaNSoEQkJCfmGPicnh4MHD2IymTh79ixdu3bl7bffJiUlhdTU1BI/Q2J6IsPWDONS6qUyfXa9OHLlCOvi1gEoI19FUT36MlK7dm3Gjh17U/v06dN56qmnaNGiBT4+PixevBjQfPfR0dFERkZy991359fyadq0KW+88Qb33XcfJpMJd3d35s6dS506dSzSMWbMGEaPHk3z5s1xc3Nj0aJFeHp6lrjP/v37mTBhAi4uLri7u/PJJ58AMHLkSHr16pXvq2/dujWNGzcmLCyMDh06FHu8YcOGMWrUKLy9vYvsKRfF+++/z5AhQ3jzzTfp1asX/v5Fp95HRkZy/fp1atWqle9iefzxx+nTpw/NmzcnKiqKxo0b37RfWFgYAwcOpFmzZtStW5fWrbVEbA8PD1avXs3YsWNJSUnBYDDw/PPP07BhQ4YMGUJKSgpSSsaOHUtAQECJn+FC6gXWHFnDgMYD6Ne4n0WfW0/e2P4GG09s5OS4k/h5VF4YqsJ+cYg5Yw8fPkyTJk10UqSwJunp6Xh7eyOEYMWKFSxfvrzI6Bl7o/B3MCUzxWHqw2QaMjl65Sgtb22ptxRFJVLSnLGqR6+wKTExMTz77LNIKQkICGDhwoV6SyoXeUZ+88nNuAgXuta1bhKbNfjtzG+0q9kOLzcvZeSrOMrQK2xKp06d2Lt3r94yrIJJmhi/aTzVPKrRJaKLXdXEOX/9PD2W9GBM1Bje7fmu3nIUOqMMvUJRTlyEC2sHr8Xfy9+ujDxAzWo1WfHwCjqGd9RbisIOUFE3CkUFCPMPo7pndQwmAx/9+RHZxmxd9VzLusbBywcB6Ne4H0E+QbrqUdgHytArFFZg6z9bGfvTWL47ou/A8v+t+z86fdGJ5MxkXXUo7AvlulEorECP+j3Y9fQuomoWGfRgM/7T7T8MbDqQAK+SQ0QVVQvVo7cQV1dXWrVqRcuWLWnTpg2///57uY5TWqXI6OhoWrRowZw5c3j11Vf5+eefLdrPWgwbNozVq1dXeJuqSJ6RP3H1BOPWj7NZAbS07DQW/L0AKSVh/mH0adTHJudVOA6qR28h3t7exMbGArBhwwamTJlyU9EwS8hLGPLxublWysWLF9m1axfHjx8v034K+2LTyU0s3b+Ul+5+iXD/mye7tzZfxH7BhP+5PAAACdNJREFUuJ/G0bZmW1rd2qrSz6dwPBzT0HfpcnPbwIEwZgykp0OhgmMADBumva5cgUceKbhu27Yynf7atWv5xckAZs2axapVq8jKymLAgAG89tprpKWlMXDgQOLj4zEajbzyyitcunQpvyRwcHDwTXXl77vvPs6dO0erVq346KOPWLBgAb179+b8+fMl7hcREUF0dDTr16/Hzc2Nzz//nClTpnD8+HEmTJjAqFGjkFIyceJE1q9fjxCCl19+mUGDBiGl5LnnnmPTpk2EhYXh4eGRf9yYmBhefPFFUlNTCQ4OZtGiRflZqoriGRU1iv6N++cXzEtISyDEN8Tq5zGYDLi5uDGm3Rja3qaMvKJ4HNPQ60BGRgatWrUiMzOTCxcusGXLFgA2btxIXFwcf/31F1JK+vbty/bt20lISKBmzZr59etTUlLw9/fnvffeY+vWrQQHB990jrVr19K7d+/8J4cFC7QJucaOHVvifgDh4eHExsbywgsvMGzYMHbs2EFmZibNmjVj1KhRfPPNN8TGxrJ3716uXLlCu3bt6Ny5Mzt37uTo0aMcOnSIS5cu0bRpU5566qn8sr7fffcdISEhrFy5kqlTpzpsgpOtyTPyi2MXM+6ncex6ehcNghpY7firDq7i9e2v88uwX6jhXYO7wu6y2rEVzodjGvqSeuA+PiWvDw4ucw8eCrpudu7cyRNPPMGBAwfYuHEjGzduzK+pkpqaSlxcHJ06dbKoJLC16Nu3L6CVIk5NTaVatWpUq1YNT09PkpOT+e2334iOjsbV1ZVbbrmFe+65h127drF9+/b89po1a9KtWzeg5LK+CsvpEtGFQZGDqBdYD7jRC68otarVola1WhbX41dUbRzT0OvMXXfdxZUrV0hISEBKyZQpU3jmmWdu2q60ksDffvstr732GqBNw1dcb90S8gqaubi4FChu5uLiUuYSwVByWV+F5dQJqMNnfT4DtEHT9vPbM6nDJJ5o+USZjiOl5Ln1z1G7em0md5xMh/AO/DTkp8qQrHBCVNRNOThy5AhGo5GgoCB69uzJwoUL80vbnjt3Ln/SjNJKAg8YMIDY2FhiY2OJiio5LK88JYHN6dSpEytXrsRoNJKQkMD27dtp3749nTt3zm+/cOFCvv+/uLK+ivKTnpNOi1tacJuf9mSUmJ7Il3u/LDbm/cTVE/lTYAohuJpxtUxTXyoUeagevYXk+ehB610tXrwYV1dX7rvvPg4fPsxdd2k+Uj8/P5YuXZo/EFpaSWBLKe9+eQwYMICdO3fSsmVLhBC888473HrrrQwYMIAtW7bQtGlTwsPD8z9HcWV9IyMjy3xuhUaIbwjLH16ev7x472Je2vgSJ8aeIMArgMWxi3nn93f4Y/gfVPOsxoI9C5j1+yyuTb6Gt7s3yx5aZnelFhSOgSpTrFBYQGV8B40mI0euHKFpSFOEEPxw7AcW713Mon6L8PXw5WzKWTIMGTSo0UAZeEWplFSm2CLXjRCilxDiqBDiuBBichHr5wghYnNfx4QQyWbrhgoh4nJfQ8v/MRQK58LVxZXI0Mh8I967YW++evSr/FmgwvzDaBjUUBl5RYUp1XUjhHAF5gI9gHhglxBirZTyUN42UsoXzLZ/Dmid+74GMA2IAiQQk7tvklU/hUKhUCiKxZIefXvguJTypJQyG1gBlDR/WjSQ54jsCWySUl7NNe6bgF7lEWpvLiZF1UF99xSOjiWGvhZw1mw5PrftJoQQdYC6wJay7lsSXl5eJCYmqh+cwuZIKUlMTMTLy0tvKQpFubF21M1gYLWU0liWnYQQI4GRQP7k2ebUrl2b+Ph4EhISrCJSoSgLXl5e1K5dW28ZCkW5scTQnwPCzJZr57YVxWD4//buL0TKKozj+PeHbXYR4ZZk0kq70EJsfy68CMEbqbA00S4ijCgrIYINDIRKvbMuiiAromCpwEAw6Q8uUdRmdqllmolaEURBrFnYPxCCpV8X5wjTsLvutjP77nvm+cDCnPe8A8/DsM+875mZ8zDY9NwVTc/9tPlJtoeAIUjfumme7+rqoq+vbwqhhhBCaDaVpZvPgX5JfZIuJBXz4eaTJF0DdAONP6X8EFgpqVtSN7AyHwshhDBLzntFb3tM0iOkAj0PeN32cUnbgUO2zxX99cBuNyyk2z4j6UnSmwXAdttnWptCCCGEydTiB1MhhBAmN9kPpuZcoZf0C/BD1XH8DwuBX6sOogKdmHcn5gydmXedcr7K9riND+Zcoa8rSYcmejctWSfm3Yk5Q2fmXUrOsXtlCCEULgp9CCEULgp96wxVHUBFOjHvTswZOjPvInKONfoQQihcXNGHEELhotCHEELhotC3iKTNkixpYR5L0ou5WctXkpZWHWMrSXpW0tc5t3clLWiY25Lz/kbSrVXG2Wrna8JTAklLJO2XdELScUmb8vFLJY3kJkIjeVuT4kiaJ+mIpPfyuE/Swfyav5m3gqmVKPQtIGkJaR+fHxsOrwL6899DwCsVhNZOI8B1tm8AvgW2AEgaIG2HcS2p98DLuXlN7TU04VkFDAB353xLMwZstj0ALAMGc55PAPts9wP78rhEm4CTDeNngB22rwZ+AzZWEtUMRKFvjR3AY6QuWuesA95wcgBYIGlxJdG1ge2PbI/l4QHSzqSQ8t5t+2/b3wPfkZrXlGC6TXhqyfao7cP58V+konclKded+bSdwB3VRNg+knqA24FX81jATcBb+ZRa5h2FfoYkrQN+sn20aaolTVdq4kHgg/y45LxLzm1cknpJrUEPAotsj+apU8CiisJqp+dJF23/5PFlwO8NFzW1fM1b3XikSJI+Bq4YZ2obsJW0bFOcyfK2vTefs410q79rNmML7SfpYuBt4FHbfzY2KbdtSUV9N1vSGuC07S8krag6nlaKQj8Ftm8Z77ik60mtE4/mf4Ie4LCkG5lew5Y5aaK8z5F0P7AGuLlhe+ra5z2JknP7D0ldpCK/y/Y7+fDPkhbbHs3LkKeri7AtlgNrJa0GLgIuAV4gLbtekK/qa/max9LNDNg+Zvty2722e0m3dUttnyI1Z7kvf/tmGfBHw21v7Um6jXSLu9b22YapYWC9pPmS+kgfRn9WRYxtMKUmPHWX16VfA07afq5hahjYkB9vAPbOdmztZHuL7Z78v7we+MT2PcB+4M58Wi3zjiv69nkfWE36MPIs8EC14bTcS8B8YCTfzRyw/XBuSrMHOEFa0hmcbg/huWqiJjwVh9UOy4F7gWOSvszHtgJPA3skbSRtJX5XRfHNtseB3ZKeAo6Q3gRrJbZACCGEwsXSTQghFC4KfQghFC4KfQghFC4KfQghFC4KfQghFC4KfQghFC4KfQghFO5f0UfoSkLhp9IAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a fitting object representing a Gaussian and set guess parameters.\n", "gf = fuf2.GaussFit()\n", "# See what parameters are available\n", "print(\"List of available parameters: \", gf.availableParameters())\n", "print()\n", "\n", "# Set guess values for the parameters\n", "gf[\"A\"] = -10.0\n", "gf[\"sig\"] = 15.77\n", "gf[\"off\"] = 0.96\n", "gf[\"mu\"] = 7.5\n", "print()\n", "\n", "# Let us see whether the assignment worked\n", "print(\"Parameters and guess values: \")\n", "print(\" A : \", gf[\"A\"])\n", "print(\" sig : \", gf[\"sig\"])\n", "print(\" off : \", gf[\"off\"])\n", "print(\" mu : \", gf[\"mu\"])\n", "print()\n", "\n", "# More convenient overview of parameter status\n", "gf.parameterSummary()\n", "print()\n", "\n", "# Exporting and assigning parameter values via dictionary\n", "ps = gf.parameters()\n", "print(\"Parameter name/value dictionary: \", ps)\n", "# Assigning values from dictionary (identical values here)\n", "gf.assignValues(ps)\n", "\n", "\n", "# Evaluate the model for the starting values\n", "startmodel = gf.evaluate(x)\n", "\n", "# Which parameters shall be varied during the fit?\n", "# 'Thaw' those (the order is irrelevant)\n", "gf.thaw([\"A\", \"sig\", \"off\", \"mu\"])\n", "\n", "# Show values and names of free/frozen parameters (i.e., parameters not varied in a fit)\n", "print(\"Names and values of FREE parameters: \", gf.freeParameters())\n", "print(\"Names and values of FROZEN parameters: \", gf.frozenParameters())\n", "print()\n", "\n", "# Use scipy's fmin to minimize chi square \n", "fr = sco.fmin(gf.chisqr, gf.freeParamVals(), args=(x,y,yerr), full_output=True)\n", "print()\n", "print(\"Fit results: \", fr)\n", "# Set the parameter values to best-fit values\n", "gf.setFreeParamVals(fr[0])\n", "\n", "\n", "# Use a convenience function to carry out optimization\n", "# using scipy's fmin algorithm (parameter values are set to\n", "# best-fit result)\n", "fr = fuf2.fitfmin(gf, gf.chisqr, x, y, yerr=yerr)\n", "\n", "# Get a summary of current parameters from the model\n", "gf.parameterSummary()\n", "\n", "# Let us see what we have done...\n", "bestfitmodel = gf.evaluate(x)\n", "\n", "plt.plot(x, y, 'bp')\n", "plt.plot(x, startmodel, 'g:', label=\"Model for starting values\")\n", "plt.plot(x, bestfitmodel, 'r--', label=\"Best-fit model\")\n", "plt.legend()\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 }