{ "metadata": { "name": "", "signature": "sha256:1a43a8da38a744c132dcd876607636b7276532690229ee93b95ea26435386404" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting Data to Theory\n", "\n", "From an experimentalist's point of view, comparing experimental data to theory is a fundamental skill. In this tutorial we will introduce fitting methods in python, simulating some 'experimental data' using random numbers.\n", "\n", "We will be using the *curve_fit* method, which is part of *scipy* and based on a Marquardt-Levenburg method of least-squares minimisation (see Ch. 7 of *Measurements and their uncertainties*, Hughes and Hase (2010), OUP). \n", "\n", "Let's generate the data - imagine we have some data that we think fits a power-law and we want to know exactly what the power is." ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from IPython.display import display\n", "\n", "# x-data is evenly spaced - this might represent the time axis \n", "# on an oscilloscope trace, for example\n", "x = np.linspace(2,80,100)\n", "\n", "# y-data roughly follows a quadratic dependence \n", "# but with some extra noise added\n", "a = 1\n", "b = -0.5\n", "c = 2000\n", "y = a*x**2 + b*x + c + np.random.randn(len(x))*450\n", "\n", "#y-error bars also have some random element\n", "y_er = 450 + np.random.uniform(-50,50,len(x))\n", " \n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "ax.errorbar(x,y,yerr=y_er,\n", " # options for colors of errorbars and point markers...\n", " color='k',marker='s',mfc='w',mec='k',mew=1.25,linestyle='None')\n", "plt.draw()\n" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+QHHd14D/PrIkt7GXlmGhly6A9LBLEwQYcEEmseEkp\nlMMltst3JeOquDDypYrzOdoNl5xlrnwrF3UOdiVk11wFX3InYnOxg0IIZYJP2E5Yal0pew0Y2WAU\n5LtdwiqS+CUxJiZ3O+bdH909+m5v//h2T89Oz8z7VE2p5zvf7n49o33f933f931PVBXDMAxjsDir\n2wIYhmEY648pf8MwjAHElL9hGMYAYsrfMAxjADHlbxiGMYCY8jcMwxhAcpW/iEyKyLMi8lURmQzb\nLhCRR0XkGyLyiIiMOP1vE5GjInJERN7ptF8WXueoiMx25nEMwzAMHzKVv4j8S+DfAm8FxoFfE5HX\nAvuAR1X1dcDfhO8Rke3AdcB24Ergj0REwst9FLhJVbcB20Tkyg48j2EYhuFBnuX/M8CTqvrPqvoS\n8AXgXwNXAfeFfe4DrgmPrwYeVNUVVV0Cngd2iMhm4HxVXQj73e+cYxiGYawzecr/q8DO0M2zAXgX\nsAXYpKonwz4ngU3h8UXAsnP+MnBxQvuxsN0wDMPoAkNZH6rqERG5C3gE+CfgK8BLsT4qIpYjwjAM\no4fIVP4AqnoAOAAgIv+FwII/KSKjqnoidOl8O+x+DLjEOX1L2P9YeOy2H0u6nw0khmEYxVFVye91\nBp9on58K/301cC3wAPAQ8J6wy3uAT4fHDwHvFpGXi8gYsA1YUNUTQENEdoQLwDc45yQ9RK1f09PT\nXZfB5DQ5Tc7+kXF6ejpRF05PT3Pq1CkWFxfXvE6dOtU6vwy5lj/wSRH5SWAFuFlVfyAiHwIOishN\nwBKwO1Taz4nIQeA5oBn2jyS7GfhT4FzgYVU9VEpiwzCMPmNqaoobb7wRgLGxMRYXFwEYGRlhZmaG\nO+64Y80509PT7N+/v/Q9fdw+v5TQ9n1gV0r/O4E7E9q/BLyxhIyGYRh9zcjICCMjre1SbN26tXWc\nNTC0g4/lb8SYmJjotghemJzVYnJWSy/IWQcZswaGdpCy/qJOISJaN5kMwzDWCxFJ9eOnfRa2V7vg\naxiGYfQfpvwNwzAGEFP+hmEYA4gpf8MwjAHEon0MwzC6SLPZZHl5eVXb0tJS63jLli0MDVWvqk35\nG4ZhdJHl5WXGxsZWtbnv5+fn2bLlTHacqgYGC/U0DMPoIktLS4yNja1R8svLy+zcuTPz3MXFRbZu\n3Voq1NMsf8MwjBqwZcuW1A1cZQaGPEz5G4Zh1Jy0gSG+VlAEU/6GYRg9SjvWvyl/wzCMGhC34n2s\n+sgdFF8w9sGUv2EYRg3IsuLTBoasdYI8LNrHMAxjnTh9+jSnT59e1dZsNnnxxRcZHh4GVqdtzlvY\nbSfax5S/YRjGOrF///7cwixu5s74BjB3YIAzcf6m/A3DMGqMa/nHC7NEOfvXK6WzKX/DMIwukKPI\ny5xj+fwNwzCMbHKVv4jcJiJfE5FnReQBEfkJEblARB4VkW+IyCMiMhLrf1REjojIO532y8JrHBWR\n2U49kGEYRr/TbDZZWlpqvcqQ6fYRka3A3wKvV9X/KyKfAB4G3gB8V1XvFpFbgY2quk9EtgMPAG8F\nLgYeA7apqorIAnCLqi6IyMPAPap6KOGe5vYxDKNnSYrogbW1eNtx+0T5gFyqdvs0gBVgg4gMARuA\nfwSuAu4L+9wHXBMeXw08qKorqroEPA/sEJHNwPmquhD2u985xzAMo2+YmZlhbGxszWtmZqbye83P\nz6+K/ilC5iYvVf2+iPwB8A/Aj4DPqeqjIrJJVU+G3U4Cm8Lji4AnnEssE8wAVsLjiGNhu2EYRl8x\nNTXFjTfeCKyN6EkjPluIXDnROfHPqtjklan8ReS1wBSwFfgB8Bci8htun9ClY34awzAM1rp3Ir77\n3e/y1a9+FaC1oStS8jMzM8zOnlkKjVw609PTAKv2BpRJ5ZBEXnqHnwP+TlW/ByAinwJ+HjghIqOq\neiJ06Xw77H8MuMQ5fwuBxX8sPHbbj6XdNNrsADAxMcHExITPsxiGYdSOLGWdVbQFzlj+0Uwi4jOf\n+Qx79+5lZmYmc0aRRd6C7zjwZwQLuP8M/CmwALwG+J6q3iUi+4CR2ILv2ziz4HtpODt4Etgbnv9Z\nbMHXMIw+R0Rabp9I0afl5o9SNfgQLwAzNjZWbTEXVT0sIvcDXwR+DHwZ+GPgfOCgiNwELAG7w/7P\nichB4DmgCdzsaPKbCQaPc4GHkxS/YRhGL+Ab0QOsUejt+OnjtJPS2Xb4GoZhFMQnRw+sDtsUCQzz\nuIUfWfFFLP+knD+W3sEwDKMi0ix8V6mn5eiJ96tS+afIYzV8DcMwqmBmZsbLwo+UdrTz1iX+vkzR\nlk5glr9hGEYKRbNwJu289cXH8k+biVS+4GsYhjHIxN04aco5UsqRFZ8W0eO2J+XmzyNtJlIGy+pp\nGIbRJlFKhyj6ZsuWLV7KfOvWra3X0FC+LT41NcXi4mJr0HCPi2JuH8MwDA+yErGdOnWqZflHMfuQ\nvxu3HV0XX0w2t49hGEYCRWLzi5J1jTQXEKxdDIb0QSbpHmXTOYMpf8MwBgTfyJ084jH2cEYJR+3u\n51mbupJmBldccQVf+MIXvORsJ8+PuX0MwxgIfCJ3ssjKpZ9FWlx/9FlcniJ7CNzzze1jGIaRgG/k\nTkSahZ8X0QOBUnffp+HKkCZPlpztpIkw5W8YhpHA8vLyGgvffZ/lznHb67KpK44pf8MwjAziFv5T\nTz3F7t27vc/Psv7bWbBtF1P+hmEYGcQtfF/LPRowXL/85ORkYtEWH7IWmstgyt8wDKMEkSJuNpuc\nOHGC48ePtz5zlfKWLVtaG7j279/P1NTUmmv5DAJ5bqiimPI3DKPvybKaG40GGzZsWLPDttFoZF4z\ny53jKuXl5eXWzKGKPQVZC81FMOVvGEbfU8ZqnpycbJ3r0mw2AVZV0XLfu/fcuXNn5Qu8VRWDMeVv\nGMbAkJdwzY2rbzQazM7OplrVcSWcppSTzk+biUQupNHR0dZMJL6BrCpM+RuGMTCkKWi33c3N7yZN\nK5OFE9bOEKDYTKQdv34WpvwNwzASGBoaWjNQlHG3ZLlp0sJIs/z6Ve0byFX+IvLTwJ87Tf8CuB34\nn8AngNcQFnFX1dPhObcBe4CXgL2q+kjYfhlBEfdzCIq4T5aS2jAMo0MkhU/6LNKmKWW33V1khvQw\n0qwBo52i7S65yl9V/x54M4CInAUcA/4K2Ac8qqp3i8it4ft9IrIduA7YDlwMPCYi28KEPR8FblLV\nBRF5WESuVNVDlTyJYRiGg5vLJ0kRJ72HZDfL9PQ0kB1Xn6aU3fbo2tFiclHSisGUcQ0VdfvsAp5X\n1W+JyFXAFWH7fcAcwQBwNfCgqq4ASyLyPLBDRL4JnK+qC+E59wPXAKb8DcNYQ7spmJOyePpYzUnJ\n1kZGRrjjjjtSlWza2kBWGKm72cuX+IxgPXP7vBt4MDzepKonw+OTwKbw+CLgCeecZYIZwEp4HHEs\nbDcMw1hDuymYp6amuPHGG4FAIc/PzwMwPDzM+Ph46uJtlnJ1B4boeu6AEg1KeUo5LYIn2ii2HvmA\nvJW/iLwc+HXg1vhnqqoiUlkeZveHnZiYYGJioqpLG4bRI8SVt2uF+xCfIVx++eWrPvexmtNmH5Ds\nzoncQ774uIqSmJubAyhUhyCOdz5/Ebka+HeqemX4/ggwoaonRGQz8HlV/RkR2Qegqh8K+x0CpoFv\nhn1eH7ZfD1yhqu+L3cfy+RvGAJLl5tm4cWNlJQ/d9/F4+7gf/fbbb+eDH/xg4jWTauf6ypp233ic\nf1J4aeRCareMYxHl/+fA/1LV+8L3dwPfU9W7QoU/oqrRgu8DwNsIF3yBS8PZwZPAXmAB+CxwT3zB\n15S/YQwm+/fvT3Xz3HHHHR1R/nmFWQ4fPszw8PCa9rB4ite9ysjnc611Uf4i8goCy31MVV8I2y4A\nDgKvZm2o5wcIQj2bwKSqfi5sj0I9zyUI9dybcC9T/oYxgGRV2vK1/NNmD3FlHVf+aXH18Spc8fOT\n6BXl7+XzV9V/Ai6MtX2fIPonqf+dwJ0J7V8C3lhEQMMwBoOilbaSSFskzqOqfDmdJj64tZPS+awK\n5DEMw6gFU1NTLC4utmYN7nE/MDMzw9jYWMtV5R4XxdI7GIbRN2TNHpKsZnfR1V2ETdoU5i62dgs3\nAsqlzADgveC7XpjP3zAGh6I++iLEfeLRwnESSTH7cVz/f1yerPUKn9DUMj7/hH7V+/wNwzA6QVkf\nfR5Rzn3XJ37ttdeya1ewTBmFUkYLu67S983L7177wIEDq0JC3bj/tFj8NP99FQVffDDL3zCMrpFm\nMcfj27Pi3ZPIC+GMrPi4qycpwifvWpAeEpqlyNNCW2+99Vbe975g+5PvLMIsf8MweoosH31Wvvu0\nEMw4aVZ8o9EoFSmTdr3h4eHC0UJp/vt777131bP6zCLKYMrfMIza0m692rQQzgMHDlSSWK0d0iz5\nffv2tSz/eP8qMeVvGEbXSfLRQ+fi7/fs2cPU1BRwxrWS5tvvRFK1LNbL52/K3zCMrhMpWN+Qxai/\nqyiTQjXTiLtp3OOqiqXUHVP+hmHUhsjNk+feiT5z/eBJdXGLWPGReylp0bkus4IqMeVvGEZtSCtr\nGH8fDRJJ7hE3Zr+IFR9FD2Xl8++nWYEpf8Mwakuass1aC9iyZQtbtmxZY7W7i8fuIrIP7vUgOfS0\n1zDlbxhG7YiUrev+cd0xecp2aGhozeCQNWC4C81JIaBJ1+uFRHBZmPI3DKMr5OXSieMq26pz7CTF\n1fc7pvwNw+gKSQu0VfjUyyzKJmX+7PdBwJS/YRhdZX5+ntHRUU6cOAEERcx37969ykdfRBGXGUB6\n3YVTBlP+hmF0lcgXf+mllwJnfO6+G7yi/EDNZrMV6QO08vTA6nWCeP1c957Rfbudunk96P8nNAyj\nr8nKDJq0TpDkbiqTN6jXMeVvGEZP4yZIi2cGzcIndfN6p1leT7yUv4iMAP8deAOgwHuBo8AngNew\ntoD7bQQF3F8C9qrqI2F7VMD9HIIC7pMVPothGD1I3gJtUt6fuJsmLTNo0jmNRqN1XpKF764ZdCqj\nZh3wtfxnCZT1vxGRIeAVwH8CHlXVu0XkVmAfsE9EtgPXAduBi4HHRGRbmKT/o8BNqrogIg+LyJWq\neqjypzIMo2dIW6CN0i4n5f3xddMknTM5mW1zxmcEEFj+3S6+UjW5yl9EXgnsVNX3AKhqE/iBiFwF\nXBF2uw+YIxgArgYeVNUVYElEngd2iMg3gfNVdSE8537gGsCUv2EMIFm7ZsfGxvjkJz+5qjpWmfTO\nSSGcjUYjM51z2owgXnyl12cFPpb/GPAdEfkYMA58CZgCNqnqybDPSWBTeHwR8IRz/jLBDGAlPI44\nFrYbhjGA5O2aff/738+ePXtaSr5Meuek/pHFXnQ/QFrxlV60+sFP+Q8BbwFuUdWnRGSGwMJvoaoq\nIpXVXnRH0YmJCSYmJqq6tGH0LGnFznvV7ZBH2edK8vMnXavofoA6fc9zc3PMzc21dxFVzXwBo8Ci\n8/5y4LPA14HRsG0zcCQ83gfsc/ofAnaE1/m60349cG/C/dQwjLVMT08rQcDFqtf09HS3RasEQBcX\nF3VxcXHV8eHDh1vvXdx+7jXyvqOVlZXWteP3Wlxc1JWVlfV54AoJ9WauPndfuZa/qp4QkW+JyOtU\n9RvALuBr4es9wF3hv58OT3kIeEBEPkzg1tkGLKiqikhDRHYAC8ANwD2+g5RhDDppIY11sUarIGuB\n1tdNk+Tnd7+jfkzSVgYJBo2cTiLjBKGeLwf+N0Go58uAg8CrWRvq+QGCUM8mMKmqnwvbo1DPcwmi\nh/Ym3Et9ZDKMQUZE6KW/Ex+XlYikLtCOj4+nXtuN9inzvfTad5lE+AxS6Jy6PbQpf8PIp9cUVjxS\nJsKNlEl7png6hrGxsVYah2gPwNDQEMPDw4yPj6/Js5+XqqHXvsskTPkbxoCQpbDquDDsyhR3WbmW\nv6cnwvu+PqkaBlX5W3oHw+gz0nLddDMePT7wRAr59OnT3tE5LvEUDmX2AAw6ZvkbRg/ia/mnWdll\nqWJW4cru4w7KOj+aBcQt/KWlpdazm+WfjFn+htFnpFnZVVD1rKLfNk71Eqb8DaMHqEsO+qrDTeu0\ncWrQMOVvGD1AXXLQd3JW0Q5FUzX0W5K2MpjyN4weohcXNtdj1lL0O4i7r3o9SVsZTPkbRg9RJrlZ\nt+n0rCUtMyiwJjVzhK01mPI3DGOdaGfWkuamgbWuJ5+BZJDcO2mY8jeMPqIuC8NJpM1aInmzFHKa\nmwZWP5/hjyl/w+gh8hY267IwXITI+s/yt0dumqTZQl6tXiMZU/6G0UP4ukl6aWE4kjXLDROfFfTS\n89UVU/6G0QNklTyMPo/3r9rCb9ellDZrKSNrLy581w1T/obRA6xHDvq81A3tupTMMq8XZ3VbAMMw\n6sHMzAxjY2NrXjMzM6v6zc/Ps7i42HpF6ZXTiGYt0QtYdZwWjml0FrP8DaMPKbrjFfxTNxR1ueTN\nWspEH2U93yDu1i2DKX/D6EPKuFjqmrohiaznG8TdumUw5W8YfUDkr282m6vcMDt37kxdGM7y8XeK\nqPJWUg7/SL6smUDSwvfhw4dpNBoAjI6Ots43qz8bU/6G0QekpVqGdAs+Kz1zFmVcSvG+SYXaIX/R\nOMmF9KY3vcn7/sYZfAu4LwEN4CVgRVXfJiIXAJ8AXsPaAu63ERRwfwnYq6qPhO1RAfdzCAq4Tybc\ny4q5GEYO8QIkaQVcxsbGUvs1Go2WxezOEEZGRti4ceOaAidRgZQ0fDaQRddIi9MvugmtHwqxVEEn\ni7koMKGq33fa9gGPqurdInJr+H6fiGwHrgO2AxcDj4nItlCjfxS4SVUXRORhEblSVQ8VEdgwBpW8\nNMQ+/voyM4SIInsN4rJGA80LL7zQ6lvnNYVBoIjbJz6qXAVcER7fB8wRDABXAw+q6gqwJCLPAztE\n5JvA+aq6EJ5zP3ANYMrfGDjKlEOsIg1xWkRPkkWflDOnikHGqAdFLP/HROQl4L+p6p8Am1T1ZPj5\nSWBTeHwR8IRz7jLBDGAlPI44FrYbxsBRphxiFWmIs5R3XNknDQhZ8rkD2rXXXsuuXbuAwKUULcq+\n8MILvOtd7/KW1+gcvsr/F1X1uIi8CnhURI64H6qqikhljjf3P9fExAQTExNVXdroM6ooKN4NypRD\n7PQzZfnzffLvZFn70aJsNMC0s2hswNzcHHNzc21dw2vBd9UJItPAD4HfJFgHOCEim4HPq+rPiMg+\nAFX9UNj/EDANfDPs8/qw/XrgClV9X+z6tuBreLN///5KC4p3g6oXLUVklTvHHVhc5e3eVyTw6hZZ\niC2z6FzFonGWDINKmQXfXOUvIhuAl6nqCyLyCuAR4A5gF/A9Vb0rVPgjqhot+D4AvI1wwRe4NJwd\nPAnsBRaAzwL3xBd8TfkbRUhTOHW3/F06ofyTiA+ISco/rnwjZe2j/NM+c4/jyeGSFo2L7Pg15R/Q\nqWifTcBfhf85hoA/U9VHROSLwEERuYkw1BNAVZ8TkYPAc0ATuNnR5jcThHqeSxDqaYu9Rlv00q7U\nTpGUbdPd6BVtfOrEYJjmdkujigR1Vny9GnKVv6ouAj+b0P59Aus/6Zw7gTsT2r8EvLG4mIZhpJGU\nbdNNf5BktSfttM27B6zeGwCB2212draE1OWx4uvVUNjn32nM7WOUpVddAO3KXWbjVJrvvWyRlMnJ\nSWZnZxP9/L7uIV96dZG/k3Ryk5dhGDUnbeNUo9FgaWmJZrPJiRMnADh+/DhwRtlHSj5N0Uf9okEG\nYHh4mOHhYSBQvLOzs2zdurXt/D15DLKSrxJT/oaxzqx3kfUDBw6kumaiASPaveta+pGS37lzJ6Oj\no61z3v72t3tV7Cqbv8dYH0z5G8Y6s95F1vfs2cPU1FRLsbtKPSJpIdb9fNu2bavk95HP6uzWG1P+\nhtElqlaOaRunhoeHVylr38pZReRzXT15tXkjN1SERet0B1P+Rs+y3u6Tqqk6uVmZgSNrp20R+ZJc\nPWnE3VAWrdMd6vuXYRg5rLf7pK4UybYZp52Zhjv4Rv8muZTiRG6oOGb1ry+m/I2eZ9B9y0U3Ti0v\nLzM6OtpS1MePH2f37t1rIn/ySBp8479DvD+sdUMZ3cGUv9Hz9GpueJ/kZp2IaU9T7PHvcXl5OTE8\n1NfaH6QBuBcx5W8MHHXZJOSjHMukfk4jssqTkr6NjY2tcQ+lyRdvj5/nuqHcCKOon++Cs9FZTPkb\nA0eVCrUMeUXIm81ma+H6ve99b+HUzy5Ji+JxWaJF8R/+8IetIvAQWPRxl1BEmmso7V69OjvrZ0z5\nGz1P0dzwZXLpV0mSj/5Tn/qU14BUVIHmLYq7MfvxQdFV7r7K21w9vYPl9jF6lipyw9chH5CIcOrU\nqVWpqZNSKMSLsWddL54/Py/vT1yGtH4R7nc/Pz/P6OjoqrUBd7YwNjbGyspKrcNuex3L7WMMFO2E\nONaN+HpD1Ra0j+Veds3DZwHZFH/9sF/E6FmqyA1fV9IWZdNIy3Hvpl9uhzTXmjujyFtANuqFKX/D\nqCHuIOYzoKXluJ+cnMw8z1XqWWkWfMND3WOz9uuN+fyNvqGM/74uPn9XhrQSiFmypoWvNhoNxsfH\nC22Em56eZmpqqhX5s23btlWx/GllF+Oy+tQRNqqhIzV81xtT/vWjLnHxeayn8q/yO6lC+aeRtyge\nHxQgeIa0cFjAu25vEpa/pzOY8jc6wv79+7saF+/Leir/Mt+JO2C45RBda3pkZISNGzdWpvyTCqa7\nVvzRo0db9X3dQSttcMuKOEqz/F3qZjD0Cx1T/iLyMuCLwLKq/rqIXAB8AngNYfF2VT0d9r0N2AO8\nBOxV1UfC9ssIirefQ1C8PdEZacq/friKoM5TeB/lmKQM3Q1WUbFzyH4+3+/E7TczM5Nb73Z6epo7\n7rgDVaXZbHL22Wcnuk+gXNbSdi1y35KMdXCnDRKdVP7vBy4DzlfVq0TkbuC7qnq3iNwKbFTVfSKy\nHXgAeCtwMfAYsE1VVUQWgFtUdUFEHgbuUdVDCfcy5V9j3D/qurmDfBROnhvEpQqFmDZDgDM5cTZs\n2MCLL74IBHH94+Pja6pqJVEma2m7Fnme8jc/f3coo/xR1cwXsIVAib8D+EzYdgTYFB6PAkfC49uA\nW51zDwFvBzYDX3fa3w3cm3I/NeqL+/tMT08rsOY1PT29bvKcOnVKFxcXdXFxUYHW8alTpxL7R/3m\n5+dbfRcXF3V+fr4lf9414mT9n02Tzz0nep/2SpN1cXGx0HeVJ2uZ8+PP1+3/D4NK+Lvk6nP35TNn\n/EPgd4Fhp22Tqp4Mj08Cm8Lji4AnnH7LBDOAlfA44ljYbvQw3U6TAOkhjnlWe9ampyr3CsSt3qxr\nx3fKHj58mFtuucW7IlbaPTtJ2sLw5ORkK2e/Wf31JFP5i8ivAd9W1adFZCKpj6qqiFTqp3H/aCcm\nJpiYSLy10WWKKLZO4Q5ALr2ocKKoG9/dvWmF2ateiE/bQDYyMpL5/ffib9ArzM3NMTc319Y18iz/\nXwCuEpF3ESzUDovIx4GTIjKqqidEZDPw7bD/MeAS5/wtBBb/sfDYbT+WdtM6RZAY9cZ3UTYiL+lb\n3vl596yKohWxOjnzyptdmZJff+JGcdq6UhaZyl9VPwB8AEBErgB+R1VvCBd83wPcFf776fCUh4AH\nROTDBG6dbcBCODtoiMgOYAG4AbinsLSGUYCsWHWfTKDdTP1cpiJW2swry3L3Udz9NLsyzuAd5x8q\n//+gQbTPBcBB4NWsDfX8AEGoZxOYVNXPhe1RqOe5BKGee1Puo74yGZ0nKTQyKdywjqF9SeGYPiUK\no+fwDef0ffa0UEg38yasjvv3zVpaJuKobvs0jPLYJi+jcnzTJtdhp2wWkXx5g1naJibf+HYfGeLH\nSd/x4uJiq54upFfESkqtEKduIblG9VhKZ6NjdKpI+nq7VpIygcZxI2iSNlIlVcdys2hu2LCBoaGh\nVfVvo5z87rWbzSZDQ0OJZQ+rrIhlSt5IwpS/4UWnyvDVIVw0q9JV0kaqvOpYvvd64oknEq14sIpY\nRucx5W+UIqrzmpQSGPxTD9QhXBTKzWzSznErWEX9YK1Cd99Hg0zRIuuGURZT/kYpIneGq8DyLOYi\nrLefuszMJu2ceLursPMGmWjAzMuRn+V6iu5p+fSNLOx/h+FFXNEcP34c6J+1gDhZM5siewXiVOU+\ny3M9tTv4Gv2PKX/DC99KTlURXwtwC5r7xKm3axnnzWzqQqcGX6P/MeVvZJJWJL3TSsanoHnWLKCo\nZVxmZuOzUSyrvQo6Nfga/Y8pfyOTuhRJT0sVnEea8o4r5DIzG9/Bz+3nO2AYRqcx5W+0RTvKrIhr\npmgag+i6aco7SXFPTk4yOzvrNbM5ePAgmzdvXpV/35U7Ih6/X2a25H4n7n4Cw2gHU/5GW7Tj+ql6\n0TIrl0+ctNq1s7OzXvfcvXs3ELieIBiQosEsaQAcHR0FkmcwAOedd16iko/6xo8nJ4NCeDaTMMpi\n6R3WiX7ZYu+bJsEn1NDNaZPkmklKHeGbxiBSkvEBJLpn2sASVaNKSqfgyha1uzV3faqEpT1HVsWv\npMpbjUaD8fHx1PtYtM9gYekdakwnQxe7MbBUuRZQ1aJl0vOWsYxdBV5mLaBMBI4b3RQv7h7hPl+z\n2UxciHflM4xMipb+6vSLPi3jWLTcYBE6XU7RR/Yyv5t7vax299q+90n6PtxXWgnE6LOjR4+2Sii6\nZRMj2VZWVtbIlPc8vs9R5vfs178bww86VMbRqIBOpjHodH6csqUSy5KXOiJKiJZH3DKO9gpE11ha\nWkqcLbi/zeWXX57Y3snds5Y/31gPTPn3AZ3Oj+OrjMrWk427Yp5++mkgfYPV8vKy1zPG+6TtFZia\nmkosdtI5LBPNAAAVWElEQVQJfK7da+tARm9iyt/IxVcZJS12uoW8066V5g+vevdq2l6BtJlNWbLW\nGeq4S9gYTCzapwt0supVtypqRVEyECi4KGY+jusqSooYcila0CQuj/s9ZEUMpS2Yu4Vd4uefOnWq\ncJUwV+60ojFlqWMlNWP9sEpe60S70TX9qvzzlGOWfFHIpKtEi5QyLFKhq2jlrWazydlnn+01uLkz\nnazw16p/J1P+g42Feq4T3c442QukrUNk1ZN1+8XJS9OQZHVX4WJZWlpq3cO9nqv4Xdmy1l/aLaQe\np+rrGQNGVigQcA7wJPAV4Dng98L2C4BHgW8AjwAjzjm3AUeBI8A7nfbLgGfDz2Yz7llZ+FOnaDds\ns8pnXFlZad0/Ls/i4uKqkMSqSfse4s/nvs/67oiFTMbDLF3c8En3NT8/v+r5o/Oz5EvD99pZYaMu\nVYfkdjrE1+gdKBHqmev2EZENqvqiiAwBjwO/A1wFfFdV7xaRW4GNqrpPRLYDDwBvBS4GHgO2qaqK\nyAJwi6ouiMjDwD2qeijhfponU5Wshwunit2wafgWWO8EWbtS3e8k7TtK88unFTRP2qkbfRbfFZzU\nL02+NOK7fcvsFnbvU/VmvH7ZNW60T0fcPqr6Ynj4cuBlwCkC5X9F2H4fMAfsA64GHlTVFWBJRJ4H\ndojIN4HzVXUhPOd+4BpgjfJfb9bDhZOXwybLbeBLN/K6p4WAtutucUsZ+jyHz+AWJWKDILooKYFc\nUqK5KqlaKZuSN9ohV/mLyFnAl4HXAh9V1a+JyCZVPRl2OQlsCo8vAp5wTl8mmAGshMcRx8L2rrOe\nBcTTFLSr3MoOOt3I694p5ZOUybPdBGZRIraIpARyZQuzNxqNju4NMIxO4GP5/xj4WRF5JfA5EXlH\n7HMVkUr9NK7ym5iYYGJiwus832lwWj/obK76LAVddNCJniFPCbqfd0pZ56VmziJt0TKO7ywma5Dw\nnR3F+z311FPs3r079doHDhxIjPwBW4Q1OsPc3Bxzc3PtXaTIAgFwO4HP/wgwGrZtBo6Ex/uAfU7/\nQ8AOYBT4utN+PXBvyj1KL3r4LoCl9Stzb59ziuR88SX+DD4Lop1aCExbfM37XiOZ0s5xc+wcPXo0\ncUE7uraPDGnfUdSe9ju5i8ZJr8OHD6+Sx+f/oGFUCVXn9hGRC4Gmqp4WkXOBXwHuAB4C3gPcFf77\n6fCUh4AHROTDBG6dbcCCqqqINERkB7AA3ADck3XvMvi6cNL69dLuy+gZ0ipTxS1e8K9/W5bIYm42\nm5w4cYLjx4+33C2PP/44zWaT73znO7zqVa9quXauvfZadu3atapvhPt7uFlEk2ZPaeUmof31jygX\nf/R8eQv27mcRZvUbdSPP7bMZuC/0+58FfFxV/0ZEngYOishNwBKwG0BVnxORgwRhoU3g5nBUArgZ\n+FPgXOBhTYj0aRffHDedzoWTho/f2td15VPjNuszdzOSr1ssb8CI3FpLS0tr7pkmXzwnfdmF606W\nm4wUu+u2y7q25dE3eoFM5a+qzwJvSWj/PrAr5Zw7gTsT2r8EvLGcmP2BjxIrGn3kRsZA8s7WpJnN\n7Oxsy08dv3ZVEVBpvvMsBd/JheuyBdetOpbRj9gO34IUqTvrtqUp5Oh9RNHoI58at+77pEHC1y1W\n1HURV+R5dXV9iS8K533/EWUKrhtGv2LKvyBl6s4muSTi+BQur4KsQSKirFvMTbVQNVkJ4PK+/7zB\n100S59PPMPqBgVb+UdGQtMLZWVZku5uqqixcXheqsJjTXC7xaxf5/n3XAzq5bmAYdWOglX9Swq40\nJRyPq6/CN92NXbmdIMlihrWK/Pjx44nt7vu83P7RtbuxqS2JMm5Aw6gDA/e/0v1jjf71UcJpi6Bx\n0nZ7JkXK1EWBtUuaWytNkfsUb3FdLo1Ggw0bNnjJklT6sYqw1rTNaI1GY03EUj/M4oz+p2+Uv68F\nluSz91HC8bj6NNJ2e1aRKyhrZ3KS66qbseVJvvNms8m2bds4evQoQ0NDufHyeWmgk0gq/VjFd59W\n7WtychLon1mcMTj0jfIvuhAbbXzy/QONW49pros9e/YkFvPwUcR5hcsPHDjABz/4wcRzk1xXkWIq\nU/w8Dd89AFm+80svvTS1XxJu9FHepra4Eo5ka5e0JHaNRoPZ2dm+mcUZg0PfKP8IXwus3ciNtEFj\neHjYK6Imks0lr3D54cOH2bNnj7dVGc1A3Gs88cQTrWcvMxAU3QOQt6juI0ORur+dUsJpriNL6Gb0\nKj2t/F0rNG8hNisMscimnrxNVWmDSpLFXLRweXxgyRro0mY27nt3IMhKa+wquFtuuaXQHoC8RfXl\n5eVCyrpsOKZVvTKMGEWTAXX6RYEkZ0lJwXwrPuUl4kq6VoQrY5a87mdpCcwmJycrS0DmPo/7fG5F\nrLxnLfKduM+XVqHr8OHDLRnchGxFq2D5fMdZlKl65XNt9zsu8nyGUSVUndit7sR3omYRDxVcXl5m\ndHS0ZSFHicVca7rKTT1pPuORkRFmZ2e9C5eUwX0On1q4Wf3SyFsQ7bZPPOv7dyk7Q7DFXaPX6Gnl\nn/QHmebCiSufuviMI9ZLMaY9X/z7KPo95C2IphH9PtF30ynXjO810gaxrNxKtivY6EV6Wvkn4WOB\nZSU+q8MfaxUJyMqkWJifn09MrexD2QXR6Peanp4GKKR4O4HvDCHCdgUbvUpfKX9fCyzrjzUp8iRv\nETQrYsV3w5dL2gAWd834pnH2HQg66brIC82Mvo8iircT2AKwMSj0lfIvElZZhLw9BFHEStIgkbQW\nUSQ98+TkZMttkpTfJmqPNk5FoZXRhrZ4jeCs2YK77tBuWuN4OgxfN5spXsNYH0RbtVbqgYhoXCaf\njUUignueiKTeI94veh+/RsTS0hJjY2Opi6DRBrKoXxqulesquTQZRIRTp061lGiWZZ6URqDZbHL2\n2Wd7LyarqvczQPLsKpoBZe3KzdrVW5S036xb1E0eYzAI/9+lK70EesLyL1tcJCkWv51SjXl7CBqN\nBpAeKeOe71u4PD5QFM1mCcGMyKfoS/SMaf2S7pW2izrNdz42Nta2T9xi9g2jfXpC+cdDOufn52k2\nm/zoRz/i8ccfb/WLK9DIFROn6l2ZkUIsEtaYFlWSR9loJJ+iL779fAagTiriohE5ncYGI6MXyVX+\nInIJcD/wUwQbY/5YVe8RkQuATwCvIazjq6qnw3NuA/YALwF7VfWRsP0ygjq+5xDU8Z30EdK3Xm2S\nAs3bYQrpSdF8/3AjZZgX1ugSDWhRsfOInTt3VpqLpxP0Ssz+elG3wcgwfPDRKivAb6vqV0TkPOBL\nIvIo8F7gUVW9W0RuBfYB+0RkO3AdsB24GHhMRLaFjvyPAjep6oKIPCwiV2qJQu6uD9vXDZLVL2lg\nSPrDzdtDUGRGEQ0uScXO20l/MAjUzaKu22BkGD7kKn9VPQGcCI9/KCJfJ1DqVwFXhN3uA+YIBoCr\ngQdVdQVYEpHngR0i8k3gfFVdCM+5H7gGKKz8XWXoa4Vm9YuiZoaHhxkeHgaS/3B9QyHLRMr4DmJV\nFhfPqoVbxxlHXanbYGQYPhT66xaRrcCbgSeBTap6MvzoJLApPL4IeMI5bZlgsFgJjyOOhe0dY2lp\nyUs5uhuNiu7ibDQaNJvNVfcqEy/vu5hcZSy+T/hqElUOQIZhdAdv5R+6fP4SmFTVF9xQSlVVEald\nfJvPIqprcefleXHJKzTyW7/1W3zkIx9ZFYsfDRLgH+IYKfvbb78dWB2Fk3TtrDQJSczPzzM6Otpa\nd4h2+LoKPT4LsDw2htH7eCl/ETmbQPF/XFU/HTafFJFRVT0hIpuBb4ftx4BLnNO3EFj8x8Jjt/1Y\n0v1c63tiYoKJiYlU2fJ2juZtWspyB/mUboxHIrkzgo985COrFOW2bdtax77l/Q4ePMjmzZsZHh5e\nU8gl6dpZaRJg7UAQDXxZqZ+jWYClUzaMejA3N8fc3Fxb18jd5CWBiX8f8D1V/W2n/e6w7S4R2QeM\nqGq04PsA8DbCBV/g0nB28CSwF1gAPgvcE1/wjW/yiu+ajRSO76anvE1LSUXaIxqNBo1Gg+HhYcbH\nx1dZ3Unfm7vBp+jGsDJlALOqVkXPMTMzkxmBFJ9JZMma9bxZpM2OLBrGMKqhU5u8fhH4DeAZEXk6\nbLsN+BBwUERuIgz1BFDV50TkIPAc0ARudrT5zQShnucShHrmLvbmpVZwUzVHiqzZbLbcIFGcf1Iq\nBFhtrWZtJoNyCbt8F6SzFL3PprE40SCwf//+VlnJvIEgb92hrKVu0TCGUT9qn97BtYyTfNOu8l9Z\nWWFoaKiQte/iWv7xKlUbN27MTQORZPnH7xVvT5vZQLrlnXbtPOKVz6JQWSDzPhFxS91SGRhGPejb\n9A5QzDcdUdSVkleEvBOsZ0rgIrVwXeKZNw3D6H16Rvm7FCnS3u0NUnUNi3QXb9tZFDcMozfpSeXf\nLWUUTwPRbDZ58cUXWxvDAJ555hkajQbHjx8HqgmL7MQAkjTjsBBOwxgcekb518FizsoPFDE+Pp56\nftnyfnnFXdoNmczL+OnKamGbhtEf9Izy97FK47th27WY4zHx0eIo0Fp8zsovFM99X3S2khRX75K3\nM9kX34yfYEnMDKNfqH20jxsN45vMza1+lURaUZWYHJlyxq3kvGgc38iYrH4ismowiKjC6vYpagN+\nhXUMw1hf+jLap4hvOlLq5513Xiu2Hfws5qQSjO41I3w2X3WSbi+8mpI3jP6g9srfxcc3HbkwLrzw\nwlXnplnMEUmbyaJrdlrhmh/dMIz1pqeUfxHfdBzffvGcQD60u7bQLT+6bylJwzD6j55S/utBGUs/\nbZBoNBpeVcK6lf6gbClJwzB6H1P+OWRZ9XlZLg8cOLAqE2eaRd8t94476ERJ7GB1KUlzPRlGf1L7\naJ+Ez72iUtLOSSMeoVMmP1D8PlVHxnQyl45l3jSM3qUvo31gfX3TkWXfbDaZn59vJZCD4pu0qrCa\n12sx2DJvGsZg0ROWf5pVClRu+WfhUfugcsvcLHLDMPIoY/n3hPJ3rd+4b9pNu+xaqVnpmZOs2azU\nytH7bih/21RlGEYeZZQ/qlqrVyBSOtPT0wqseU1PT5fql0Zcjjy5fPsYhmFUTah7CunanrD8XXwt\n4XYt5rgV72PVW3ETwzC6Qd+6fbokh5fyL+peMgzDqJoyyv+sTgkzKMzMzDA2NtZaLI6OZ2ZmuiyZ\nYRhGOrmWv4gcAP4V8G1VfWPYdgHwCeA1hMXbVfV0+NltwB7gJWCvqj4Stl9GULz9HILi7ZMp9+tZ\ny9/FLH/DMNaLTln+HwOujLXtAx5V1dcBfxO+R0S2A9cB28Nz/kjO5Eb+KHCTqm4DtolI/Jpd5/Tp\n0ywtLbVi6Z955hkef/xxnnnmGYDWZ3/913/dOmdkZIStW7euedVB8c/NzXVbBC9MzmoxOaujF2Qs\nS67yV9V54FSs+SrgvvD4PuCa8Phq4EFVXVHVJeB5YIeIbAbOV9WFsN/9zjm1Ie7CGR8fZ+fOna3q\nXNFnv//7v99NMb3plf+4Jme1mJzV0QsylqXsDt9NqnoyPD4JbAqPLwKecPotAxcDK+FxxLGwvVbE\nd7lGewqGh4dX1em99957uyCdYRhGdbSd3kFVVUS676SvAF8//TnnnLMO0hiGYXQQn80AwFbgWef9\nEWA0PN4MHAmP9wH7nH6HgB3AKPB1p/164N6UeyVuzrKXvexlL3ulv4pu8ipr+T8EvAe4K/z30077\nAyLyYQK3zjZgIZwdNERkB7AA3ADck3ThoivWhmEYRnFylb+IPAhcAVwoIt8C/jPwIeCgiNxEGOoJ\noKrPichB4DmgCdzsxG3eTBDqeS5BqOehah/FMAzD8KV2O3wNwzCMzlObHb4icqWIHBGRoyJya7fl\niRCRAyJyUkSeddouEJFHReQbIvKIiHQ9qF9ELhGRz4vI10TkqyKyt26yisg5IvKkiHxFRJ4Tkd+r\nm4wuIvIyEXlaRD4Tvq+dnCKyJCLPhHIu1FjOERH5pIh8Pfztd9RNThH56fB7jF4/EJG9dZMzlPW2\n8G/9WRF5QER+oqictVD+IvIy4L8SbAzbDlwvIq/vrlQtPobnJrcuswL8tqq+AXg78O/D77A2sqrq\nPwPvUNWfBd4EvENELq+TjDEmCVyY0fS4jnIqMKGqb1bVt4VtdZRzlsDd+3qC3/4INZNTVf8+/B7f\nDFwGvAj8FTWTU0S2Ar8JvCXMuvAy4N0UlbPoCnEnXsDPA4ec96uihrr9IjnaaVN4PEoY7VSnF8Ei\n/K66ygpsAJ4C3lBHGYEtwGPAO4DP1PV3BxaBn4y11UpO4JXA/0lor5WcMdneCczXUU7gAuDvgY0E\n67afAX6lqJy1sPwJIoO+5byPNofVlbRNbrUgtAzeDDxJzWQVkbNE5CuhLJ9X1a9RMxlD/hD4XeDH\nTlsd5VTgMRH5ooj8ZthWNznHgO+IyMdE5Msi8ici8grqJ6fLu4EHw+Nayamq3wf+APgH4B+B06r6\nKAXlrIvy79lVZw2G2drILyLnAX8JTKrqC+5ndZBVVX+sgdtnC/BLIvKO2Oddl1FEfo0gkeHTQGLo\ncR3kDPlFDdwUv0rg6tvpflgTOYeAtwB/pKpvAf6JmEuiJnICICIvB34d+Iv4Z3WQU0ReC0wReCQu\nAs4Tkd9w+/jIWRflfwy4xHl/CavTQdSNkyIyChDmLfp2l+UBQETOJlD8H1fVaO9FLWVV1R8AnyXw\nrdZNxl8ArhKRRQLr75dF5OPUT05U9Xj473cI/NNvo35yLgPLqvpU+P6TBIPBiZrJGfGrwJfC7xTq\n933+HPB3qvo9VW0CnyJwnRf6Puui/L9IkOlzazjqXkewYayuRJvcYPUmt64hIgL8D+A5VXWLCdRG\nVhG5MIpAEJFzCfyUT1MjGQFU9QOqeomqjhFM//9WVW+gZnKKyAYROT88fgWBn/pZaianqp4AviUi\nrwubdgFfI/BV10ZOh+s54/KBmn2fBL79t4vIueHf/S6CwIRi32e3F1acRYxfJVjEeB64rdvyOHI9\nSOBX+38E6xLvJVhweQz4BvAIMFIDOS8n8E9/hUChPk0QpVQbWYE3Al8OZXwG+N2wvTYyJsh8BfBQ\nHeUk8KV/JXx9Nfq7qZucoUzjBAv8hwks1VfWVM5XAN8lyEIctdVRzv9IMIA+S5BZ+eyictomL8Mw\njAGkLm4fwzAMYx0x5W8YhjGAmPI3DMMYQEz5G4ZhDCCm/A3DMAYQU/6GYRgDiCl/wzCMAcSUv2EY\nxgDy/wGGEzElI9a70AAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 388 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's try and fit that data. We first need an analytic function to fit to, so let's define a quadratic function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def quaddy(x,a,b,c):\n", " return a*x**2 + b*x + c" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 389 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we use curve_fit to fit the function. In most cases, there are few options that need to be passed, but there are a lot of optional keyword arguments that may help in certain cases - check out the official documentation here for more information.\n", "\n", "For now though, we just need to pass the x-data, y-data and y-error bars, along with the function we want it to fit to. For more involved functions it's a good idea to pass an initial guess at the correct parameters, but in this simple case we should be fine without.\n", "\n", "*curve_fit* returns a list of the optimised parameters, and the covariance matrix which can be used to find the error on the fit parameters (which we will do) and information about how correlated the parameters are (which we won't go into this time)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.optimize import curve_fit\n", "\n", "p_opt, p_cov = curve_fit(quaddy, x, y, \n", " sigma=y_er, #error bars array\n", " absolute_sigma=True)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 390 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Important Note!!** The *absolute_sigma* option is relatively new (May 2014), and requires scipy version 0.14.0 or later. You may have to update your version to be able to take advantage of this. You can check what version of scipy you are using by running an interactive python session (type `python` and press enter on a command line / terminal window) then `import scipy` and `scipy.__version__` (note the two underscores on either side).\n", "\n", "Setting the *absolute_sigma* option to True means that the error bars represent one standard deviation of the data. Otherwise, curve_fit uses the *sigma* option only as a relative weighting, i.e. if all errorbars are the same size then the result is the same as not having error bars at all. Needless to say, **you should use *absolute_sigma=True* in most cases!!**\n", "\n", "The errors on the function arguments are given by the square-root of the diagonal elements of the covariance matrix (Hughes and Hase, Ch 7). So we quote results as `p_opt[i] +/- np.sqrt(p_cov[i][i])`, or we can use the `.diagonal()` property of numpy 2D arrays to make things a bit easier: " ] }, { "cell_type": "code", "collapsed": false, "input": [ "p_errs = np.sqrt(p_cov.diagonal())\n", "\n", "for paramlabel, p, err in zip(['a','b','c'],p_opt,p_errs):\n", " print paramlabel+':', p, '+/-', err" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "a: 1.0057988965 +/- 0.0960793019782\n", "b: -3.16885611971 +/- 8.10014558537\n", "c: 2062.33133492 +/- 143.880942065\n" ] } ], "prompt_number": 391 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's add the best fit line to the plot. We make use of the shorthand when calling functions that arguments can be passed as lists:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x_fit = np.arange(0,x.max(),0.1)\n", "y_fit = quaddy(x_fit, *p_opt) # shorthand for passing arguments using *\n", "\n", "ax.plot(x_fit,y_fit,'b--',lw=3)\n", "display(fig)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt8HVW1+L+rTaEtJaQ8bFoCNEoBqxIFeYhWwhUEX4Co\nBe5VgfZ65ZZHcr162+LF1BcvFVNQAR+FolKpeK8CQqEI4Re8Qnm1RUqlhRMltS2lpJwij+a06/fH\nzJzMmTMzZ87JSXJOsr6fz3wys2fPzJpJsvbea6+9lqgqhmEYxshi1FALYBiGYQw+pvwNwzBGIKb8\nDcMwRiCm/A3DMEYgpvwNwzBGIKb8DcMwRiAFlb+ItIjIUyLyZxFpccv2FpHlIvKsiNwrInW++vNF\nZJ2IrBWRD/vKj3Tvs05EFg7M6xiGYRhJiFX+IvJO4F+Bo4Am4OMi8jZgHrBcVQ8B/uAeIyLTgTOB\n6cApwI9ERNzbXQfMVtVpwDQROWUA3scwDMNIQKGe/2HAI6r6hqruBB4EPgWcCix26ywGTnf3TwOW\nqGqvqnYB64FjRGQysKeqrnDr3ey7xjAMwxhkCin/PwMzXDPPeOCjQAMwSVU3u3U2A5Pc/SlAt+/6\nbmD/kPINbrlhGIYxBNTEnVTVtSJyJXAv8A9gJbAzUEdFxGJEGIZhVBGxyh9AVRcBiwBE5Ns4PfjN\nIlKvqptck86LbvUNwAG+yxvc+hvcfX/5hrDnWUNiGIZRPKoqhWv1kcTb5y3uzwOBM4BbgNuBc9wq\n5wC/dfdvB84Skd1EpBGYBqxQ1U1AWkSOcSeAP+e7JuwlKnpra2sbchlMTpPT5Bw+Mra1tYXqwra2\nNnp6ekilUnlbT09P9vpSKNjzB24TkX2AXmCOqr4iIlcAS0VkNtAFzHSV9hoRWQqsATJufU+yOcBN\nwDjgLlVdVpLEhmEYw4zW1lbOPfdcABobG0mlUgDU1dXR3t7O17/+9bxr2traWLBgQcnPTGL2+WBI\n2cvAiRH1LwMuCyl/HHhXCTIahmEMa+rq6qiryy6XYurUqdn9uIahPyTp+RsBmpubh1qERJic5cXk\nLC/VIGclyBjXMPQHKdVeNFCIiFaaTIZhGIOFiETa8aPOueXlnfA1DMMwhh+m/A3DMEYgpvwNwzBG\nIKb8DcMwRiDm7WMYhjGEZDIZuru7c8q6urqy+w0NDdTUlF9Vm/I3DMMYQrq7u2lsbMwp8x93dnbS\n0NAXHadcDYO5ehqGYQwhXV1dNDY25in57u5uZsyYEXttKpVi6tSpJbl6Ws/fMAyjAmhoaIhcwFVK\nw1AIU/6GYRgVTlTDEJwrKAZT/oZhGFVKf3r/pvwNwzAqgGAvPkmv3jMHBSeMk2DK3zAMowKI68VH\nNQxx8wSFMG8fwzCMQWLbtm1s27YtpyyTyfDaa69RW1sL5IZtLjSx2x9vH1P+hmEYg8SCBQsKJmbx\nR+4MLgDzNwzQ5+dvyt8wDKOC8ff8g4lZvJj9SUM679oFX/wifOQj8KlPmfI3DMOoCgrE5o+9RhX+\n/d/hhhtg9GjYudOUv2EYRlXQH+X/8MNw3HHQV20AkrmIyHwReVpEnhKRW0RkdxHZW0SWi8izInKv\niNQF6q8TkbUi8mFf+ZHuPdaJyMJihDQMwzD6eO97M1x99RZGj1ZOO+3Vku4R2/MXkanA/cDbVfVN\nEbkVuAt4B/CSql4lInOBiao6T0SmA7cARwH7A/cB01RVRWQFcKGqrhCRu4BrVHVZyDOt528YRtUS\n5tED+bl4+9Pz9+IBwXHAI8DOsvf800AvMF5EaoDxwN+BU4HFbp3FwOnu/mnAElXtVdUuYD1wjIhM\nBvZU1RVuvZt91xiGYQwb2tvbaWxszNva29vL/qzOzitJpdaXdG3sIi9VfVlEvgf8DXgduEdVl4vI\nJFXd7FbbDExy96cAD/tu0Y0zAuh19z02uOWGYRjDitbWVs4991wg36MniuBowQvb7F3z3HNp9tln\nV/ZcORZ5xSp/EXkb0ApMBV4Bfi0in/XXcU06ZqcxDMMg37zj8dJLL/HnP/8ZILugy1Py7e3tLFzY\nNxXqhWtoa2vj8ceP5M47Pwh8NOdcfykU3uG9wP+p6lYAEfkf4H3AJhGpV9VNrknnRbf+BuAA3/UN\nOD3+De6+v3xD1EO9xQ4Azc3NNDc3J3kXwzCMiiNOWcclbQG48863cOed4wEYP76TX/5yE+9+9w7u\nuOMOLr74Ytrb22NHFHEUmvBtAn6JM4H7BnATsAI4CNiqqleKyDygLjDhezR9E74Hu6ODR4CL3et/\nj034GoYxzBGRrNnHU/RRsfm9UA0eixfDeef1uXO+972wfDnU1eUngGlsbCxvMhdVXSUiNwOPAbuA\nJ4AfA3sCS0VkNtAFzHTrrxGRpcAaIAPM8WnyOTiNxzjgrjDFbxiGUQ0k9egB8mzySez0v/hFruJ/\nz3vg3nsdxe9nQEM6q+pVwFWB4peBEyPqXwZcFlL+OPCuEmQ0DMOoKNrb2wvG6OkPu+8Oo0bBzp3Q\n1OT0+CdO7Dvf0NCQE+OnlHkAW+FrGIYRQVQP3++LHxWjJ1hPxLHKBM07ngknWP6b38Dll8OyZbDv\nvvFyWg5fwzCMMpK0h+8p7Uwmk/Xg8QgeJ03a8qlPwemnO7F7BgLr+RuGYURQbBTOvpW3xRPs+ReS\nx0/ZJ3wNwzBGMkEzTpRy9pSy14uP8ujxl3uNyf33j2P8+F00NOxXUJ6okUgpFAzsZhiGYcTjhXTw\nvG8aGhryfPbDWLNmKv/+75OYPXsynZ2F++Ktra2kUqnsCMS/Xyxm9jEMw0hAXCC2np6ebM/f89mH\nQl44H2fMmDvo7XWO3vlOWLkyuY0/OJlsZh/DMIwQivHNL5a4e4SbgL4HLM0q/gMO6OW66zbxwgs7\nIxuZsGcEJ5OLwZS/YRgjgnL55gfz6kKfEvbK/eeDi7q2bh0F/BwY45as54UXTmDGDOea448/ngcf\nfDCRnP2J82NmH8MwRgRJPHfiyI+ln4xwv/6LgN8wdapwyy2bOO64A7PyFLOGwG9eMrOPYRhGCEk9\ndzyieviFPHrAUfj+43zuBE7lj39cxpQpB8TKEydnqeGcwZS/YRhGKN3d3Xk9fP9xXIwef3n0oq57\nmDKlHJKWhil/wzCMGII9/EcffZSZM2cmvn7GjM8Am0LP9WfCtr+Yn79hGEYMXg/f2yZPnpz4OvgC\nY8du4KabNgLQ0tKSU8dL8ZgEL3SEtwE5+8ViPX/DMIwS8Mw3mUyGTZs2sXHjxuy5rq4ufv7zPYEf\n88YbcP759cD7WbBgAa2trXn3StIAFDJDFYspf8Mwhj1x7pnpdJrx48dTU5OrDtPpdOw942LpNzZe\nC3wve/zOd8Jjjz1dljUFcRPNxWDK3zCMYU8pvWbPRBNsNDKZDEBOFi3/8Y9+tBff+U5f8P3p09P8\n7GfbaGrKX2BWCv1J2u7HlL9hGCOGQgHX/H716XSahQsXRvaqg0rYOz75ZLjmGnjzTYAHWbPmEzQ1\nbc+5Nmok4pmQ6uvrsyOR4AKycmHK3zCMEUNUr9lf7o/NH8yW5T+OCtz2oQ/BbbfB97//Gvff/xE6\nO+/NGSFAcSOR/tj14zDlbxiGEUJNTU1eQ5HU3PLxj8M73vEib33r67Fmmig30ji7ftJkMIUoqPxF\n5FDgV76itwKXAr8AbgUOwk3irqrb3GvmA7OAncDFqnqvW34kThL3sThJ3HP9ngzDMIaYMNfJqEla\nfySaoBLesCE/zo9/khnyRyJe3bgGoz9J2/0kSeD+F+A9ACIyCtgA/C8wD1iuqleJyFz3eJ6ITAfO\nBKYD+wP3icg0N2DPdcBsVV0hIneJyCmquqwsb2IYhuHDH8snLOBa2DGEm1na2tqA3IbhjTeEiy7a\nDzgd+G2kUvaXe/cO+vsnJSwZTJTMhSjW7HMisF5VXxCRU4Hj3fLFQAdOA3AasERVe4EuEVkPHCMi\nfwX2VNUV7jU343w1U/6GYeTR3xDMYVE8k/Sa/QrVH1Tt61//uk/JTgB+h2P4uJWbburh+ONfz7su\nzo104cKFBWUJEhwRDGZsn7OAJe7+JFXd7O5vBia5+1OAh33XdOOMAHrdfY8NbrlhGEYe/Q3B3Nra\nyrnnngs4CrmzsxOA2tpampqaIidv45RrKpUinR5FU9MG4H1u6W6sWrUbxx//erZRKqSUozx4vIVi\n5bLrx5FY+YvIbsAngLnBc6qqIlK2OMz+X2xzczPNzc3lurVhGFVCUHn7e+FJCI4QPvCBD+ScT9Jr\nDo4+tmwZzec/Pwk40Ffrv/j+97/D97/fZx5KShJTURgdHR0AReUhCJI4nr+InAb8u6qe4h6vBZpV\ndZOITAYeUNXDRGQegKpe4dZbBrQBf3XrvN0tPxs4XlXPDzzH4vkbxggkzswzceLE0OxWSQlmx/KO\ng/72QTv6pZdeyje/+U3fnY7EsXBPAODrX9/K5z/f58OfVNao5wb9/MPcSz0TUn/TOKKqiTYcj59z\nfMdXAXPd/XnAFe7+dGAlsBvQCDxHXyPzCHAMIMBdwCkhz1HDMEYebW1tCuRtXnl/CF7vHadSqdBn\netuqVas0lUrlbD//+UaFV3Xx4mTPKkW+JPfyn3P3E+tzVU3W8xeRPXB67o2qut0t2xtYijP+6SLX\n1fMSHFfPDNCiqve45Z6r5zgcV8+LQ56lSWQyDGN4EZdpK2nPP2r04Ga6yh4Hs3JF+dUHs3D1Xf8W\nVF8MlSEqB28cccnho+41KAncVfUfwL6BspdxvH/C6l8GXBZS/jjwrmIENAxjZFBspq0woiaJC1F8\nvJwtRT+jHAQbt/7kA7B4/oZhDBtaW1tJpVLZUYN/v3jO5rbbJpRPuDLQ3t6ekwOgmHwAQSyBu2EY\nFUfUBG2p9xARenp68sxKfvNOQ0NDdhK2vX0nCxe+jVGjlOuu28KHP/xa5GRrIdlLed+o8gJmraLM\nPqb8DcMYMoq10RdDUPm3tbVFmoS8NQCOi2UbsMB39jHgaFKp57OmoTilHJyvSOKaWorNP6Re+W3+\nhmEYA0GpNvpCeDH3/TbxM844gxNPdKYpPVdKr+fvKP1RwLXAhdlr3vOeN/jWtzJ87GMamQwGYNGi\nRTkuoZ4pJm5BWpT9vhwJX5JgPX/DMIaMqB5z0L89zt89DM+LJwrPi8fvb//441v49KffCuwDwMkn\nw29+A1u2xN8LYNWqVdTW1uaVxynyBQsWhDZ8c+fO5fzzneVPSUcRpfT8TfkbhlERBM00cUS5YHoU\ncuEMU9bOuf9i990f4vTTR3HzzbDbbrkNSbEuoXFEmbyuv/56rrzyyrzyuFGEmX0MwxhW9DdfbZQL\n56JFiyIDq/3P/2zk5JP3Z/To5Pcrhaie/Lx587I9/2D9cmLK3zCMISfMRg/lVbZ+Zs2aRWtrK5Dv\n+VNbm+KFF3qzdQciqFocg2XzN+VvGMaQ4ynYpD7rXn2/ovTb7wsp7DVr9uOllyZzoTu3629gypUs\npdIx5W8YRsXgmXkKmXe8c347eFhe3PDQyJ/hnHPq2bED9tmn75xnXgqbdPbkGYxQy4OFKX/DMCqG\nqLSGwWOvkQgzj/h99vMbkIuBX7FjhzM3On8+wO4AWe+huHj+w2lUYMrfMIyKJUrZxs0FNDQ00NDQ\nEOi1f5Czzvobv/pV3+TxoYfCsmXQ2PhmQTn894Nw19Nqw5S/YRgVh6ds/eYfvzmmkLKtqakJNA77\n8Yc/1GePjjsObr+9z+zjn2gOC5aWf7/+pVCsBEz5G4YxJIRN0MbZ0P3KNm6BVzgvctNNL3LWWVP4\np3+CJUtg3Li+s/65glIDpVUbpvwNwxgSwiZoy2FTj5onmD59B3/6Exx2GHk+/GGRP4d7I2DK3zCM\nIaWzs5P6+no2bdoEOEnMZ86cmbPAqxhFHNeAvOMd4eXVbsIpBVP+hmEMKd7k7cEHHwz02dyTLvDy\nwiRkMhk6Ozu577596eoaz+LFB4XOEwTz5/qf6T23eLNS9TH839AwjGFNbmTQrwLfcvdnh84ThJmb\n/MelxOmpRiyTl2EYVU1rayt/+UuKT31qO32KH+BC3KgRoXR2dmYzfaVSqez6gO7u7uxIoKuri66u\nrtAAbNVOop6/iNQBPwXegZPV/jxgHXArcBD5Cdzn4yRw3wlcrKr3uuVeAvexOAncW8r4LoZhVCGF\nVs2Gxf3x7++xRwPnn1/HAw/0XfOhD8Ef/tBMd/fKvGvS6TQQbVbyzxkkictfrSQ1+yzEUdafFpEa\nYA+c8dVyVb1KROYC84B5IjIdOBOYDuwP3Cci09w4zdcBs1V1hYjcJSKnqOqysr+VYRhVQ9QEbTqd\npqurKzTuj3//0Uf/SlfXgdnjWbPg+utht91eCb2mpSW+zxmMJApODKGhTr5SdlQ1dgP2Ap4PKV8L\nTHL364G17v58YK6v3jLgWGAy8Iyv/Czg+pD7qmEYw5/e3l5NpVLZDcjZv/TSSxXH0qCAdnZ25tTv\n7OzMXvPMM6p7760K83TXLuf+/vv5t1WrVmXP+fHLEEZbW1uOPN7W1tY2oN8pCa7eLKjP/VuSnn8j\nsEVEbgSagMeBVlfxb3brbAYmuftTgId913fjjAB63X2PDW65YRgjkEKrZr/0pS8xa9as7CrfOO+f\nww6DZ5+Fffe9ApHLQ+/n4fXYiw3S1trayrnnnptXXpW9fpKZfWqAI4ALVfVREWnHMfFkUVUVkbKl\n3/Lb1pqbm2lubi7XrQ2jaonK/FS1ZocChL1XJgOqMGZMfv24UA1h9yp2QVklfeeOjg46Ojr6d5NC\nQwMck07Kd/wB4PfAM0C9WzaZPrPPPGCer/4y4Bj3Pn6zz9mY2ccwElPJZodyQMDsEzTTrFz5Vz3p\nJNV/+zfVXbvCzTRh3yf4jeLMTalUSnt7ewf/5fsJA2H2UdVNIvKCiByiqs8CJwJPu9s5wJXuz9+6\nl9wO3CIiV+OYdaYBK1RVRSQtIscAK4DPAdckbaQMY6TjNzsEE3sPF6InaKfxiU/sywsvOOcmT97K\niSeGm2nCQjX4v9FwDNJWCokSuItIE46r527AcziunqOBpcCB5Lt6XoLj6pkBWlT1Hrfcc/Uch+M9\ndHHIszSJTIYxkvEnO68GkpisRCRUcd977y6++MWJwERfaRvwDSB3UVYp36XavmUYpSRwT6T8BxNT\n/oZRmGpTWAsWLPCtwu3D7z8f9k533gmnn67s3Onotd1338Wbb55NZ+dFQN8agJqaGmpra2lqasqL\ns18oVEO1fcswTPkbxgghTmFV4sSwX6agycrf8w++UzrtxN5/+mmYMgV+9zs46qjkOi5JqIaRqvwt\nto9hDDNyY930MZSrVIMNj6eQt23bFuudU1vrJF258EL46U+dBgByE7tA/sKsQjmADev5G0ZVkrTn\nH9XLLpVyjCr8sicxB8VdL+J0doM9/K6uruy7W88/HOv5G8YwI6qXXQ7KParwezDdddd4TjjhdcaN\n02HlwVSpmPI3jCqgUmLQl9vdtK6ujvHj62htheuug3/+Z/jFL0CK6sMapWDK3zCqgEqJQV/uUcWm\nTfCZz8BDDznHt9ziROScNau4+xQbqmHYBWkrAVP+hlFFVOPEZtSo5cknd2POnLewaVOfGpo5E848\ns/hnFPsNguar4Ry6OQpT/oZRRSRNbVhJRI9afoazFhRGjYLLL4evfKV4k4/fr99vigLyQjN7DLcg\nbaVgyt8wjEEhOGpZv/7vnHTSE+y1VxNLl47mwx+OvjbKTAP5pqckjeNIMu9EYcrfMIYRlTIxHEb4\nqOV4fvzj2zjkkP3Yti1aIUeZaSD3/YzkmPI3jCqi0MRmpUwM5zMuovxvnHnm0UC8vd0z04TNcQTf\n10iGKX/DqCKSTmxWysRwJgPf+U4d8DivvhpuzPdkjTPDBM00lfJ+1Ywpf8OoAhoaGoqa2ByIieFi\nTUqbN8PZZ8MDD9QBdVx88RZ+9rOu7ISud69SZK3Gie9Kw5S/YVQBgxGDvlDohmJMSg895LhtbtzY\nV/eBBx7nrW89A3i9rHIbpTFqqAUwDKMyaG9vp7GxMW9rb2/PqdfZ2UkqlcpunZ2dOefXroUTTuhT\n/CJKa2sP69cfRiq1BiB7LUS7YxoDi/X8DWMYUuyKV0geuqGQyeWww5wVuj/+sZNX95e/FE4+OTcZ\ni//6UryP4t5vJK7WLQVT/oYxDCll8rOcoRsWLoRdu+DSS+HAA0u+TSRx7zcSV+uWgil/wxgGePb6\nTCaTY4aZMWNG5MRwnI2/v4wdCz/5SX65l3krLIa/J1/cSCBs4nvVqlWk02kA6uvrs9dbrz8eU/6G\nMQyICrUM0T34uPDMcXgmlq1bRzF//r6ceWZPYjm9a8MStUPhdQhhE9+HH3544ucbfSRN4N4FpIGd\nQK+qHi0iewO3AgeRn8B9Pk7Qjp3Axap6r1vuJXAfi5PAvSXkWZbMxTAKEExAEpXApbGxMbJeOp3O\n9pj9I4S6ujomTpyYl+DES5DicDxwCzAFWAccQSr1VEFTkXePKD/9YhehDYdELOVgIJO5KNCsqi/7\nyuYBy1X1KhGZ6x7PE5HpwJnAdGB/4D4RmeZq9OuA2aq6QkTuEpFTVHVZMQIbxkilUBjiJPb6UkYI\nHg0NDaxfn+Laa+u49tq92LXL0zXTuPba52ho2DtSVq+h2b59e/Ze5qc/tBRj9gm2KqfiNP8Ai4EO\nnAbgNGCJqvYCXSKyHjhGRP4K7KmqK9xrbgZOB0z5GyOOUtIhliMMcZRHT1iIhLCYORddtD933z0m\ne7zffnDzzXDKKW+JldWoPIrp+d8nIjuBG1T1J8AkVd3snt8MTHL3pwAP+67txhkB9Lr7HhvccsMY\ncZSSDrEcYYjjRghBZR/WIHzqU78GPu0ePcDKlSdkk6r7G7QzzjiDE088EXBMSt6k7Pbt2/noRz+a\nWF5j4Eiq/N+vqhtFZD9guYis9Z9UVRWRshne/H/8zc3NNDc3l+vWxjCjHAnFh4JS0iEO9DvFBUjz\nx9/Zbz+YPBna2k5kypSd2TpxvX1vUtZrYEpZh2D00dHRQUdHR7/ukWjCN+cCkTbgVeALOPMAm0Rk\nMvCAqh4mIvMAVPUKt/4yoA34q1vn7W752cDxqnp+4P424WskZsGCBWVNKD4UlHvSUkRyzDn+hsXf\nePifK27AnSQTsapOwpVSJp1zJ43zsQnf0ihlwreg8heR8cBoVd0uInsA9wJfB04Etqrqla7Cr1NV\nb8L3FuBo3Alf4GB3dPAIcDGwAvg9cE1wwteUv1EMUQqn0nv+fgZC+YcRbBDDlH8qleKAA6byne/A\nHnvAJz7Rlf2uQaUcJ3fw3t5+MDhcWIC6Ylb8mvJ3GChvn0nA/7p/HDXAL1X1XhF5DFgqIrNxXT0B\nVHWNiCwF1gAZYI5Pm8/BcfUch+PqaZO9Rr8od0LxaiQs2qZ/oZe38ClJY/j3v4/m3HPhwQdht93g\n4IPHxNaPMrtFUY4AdZZ8vUyoakVtjkiGUTzV+rfTX7lTqZTiOGWEbqlUKu+a3t7e7Lm+6z+jtbUZ\ndQw7zvbpT6cV0M7OTk2lUrpq1Srt7OzUzs5OBbSlpSXyuUner5R3b2trC31eW1tb0fcaLrjfsShd\nW7TNf6Axs49RKtVqAuiv3KUsnMq3vf8n8N3s0ahRykUXvcLpp/+ZE04oHCeopaWFhQsXhtr5k5qH\nklKtk/wDyYDY/AcbU/5GqYx05R9U8l75qlWrqK2tJZPJsGnTJgA2btzIzJkzsw3GihVbOPPMw4A9\ngRTwL8Cfsvfy6nmNDEBtbS21tbVA7qrgTCbDmDFjQiedIdeuX62/s0pjIFf4GoZRJgY7yfqiRYtY\nuHBh6DlvpW1DQwOXX76Njo7XueeeJmB7VsnPmDGD+vr67DXHHntsrHz9jd9jDA6m/A1jkBnsJOuz\nZs2itbU1awbq6Ohk9OjcsMg1NTXMnbsvZ53VRWOjE4LBf37atGk58ieRz/LsVjam/A1jiCi3coxa\nOFVbW8vUqVN5800BruQb3ziCn/70xbzrgx6ixcjnD9VcKDdvOp0ODek8km32Q4Epf6NqGWzzSbkp\nd3CzuIbjiSfg7LMnA//F/ffDr341AYhfaVuMfGGmniiCZihLvjI0VO5/hmEUYLDNJ5VKWIIT73jH\nDrjxxgO47DLIZHbL1nnwwXFAaRm/PPyNr/fTP08QhWeGCmK9/sHFlL9R9Yx023Lcwqnrr4dvfMN/\n5h986UsbOf982LzZUdRBz5+k3y+s8Q3+HoL1oc8MZQwtpvyNqqdaY8MnCW7WX5/22bPhhhtg5Up4\n73vf4LHHmrj66ue4+ur8usHv2N3dnece6pezUG9/JDXA1Ygpf2PEUSmLhJIox1JCP/sZMwYWLYKO\nDpgzp4axY58L9b9vbGzM6bXHyRcsD17nN0N5owj/6CxY3xgaTPkbI47+KtT+UigJeSaTyU5cn3fe\neYlCP2/fDuvWwRFH5D4rk8kwcWI3n/wkuB33PFm8SfFXX301mwQenB590CTkEWUaigrNXK2js2FN\nsfEgBnqjSuOzGIOPF5PGizvjbV7cmbCYNqqqPT09OTFtvP2enp7BfQEX3Lg0FIhXE/W/cc89qgce\nqDppkurWrbn1ksb9iZMh7Fv675sbHyh5fCGjfGCxfYyRRDliw1dCeAERoaenJyc0dVgIhWAy9q1b\n4ctfhptu6rvXOefA4sV975Q07k9Qhqh6Hv5v39nZSX19fWToiMbGRnp7eyva7bbasfAOxogizsXR\nO18tBOcbCs0H3HEHzJoFL73UV7bPPnDSSbB4cX79JGaXUuc8omT1P9MUf+VhvxGjailHbPhKJWpS\n1mOffXIV/yc/uYNLLtnEvvvuAvoWu3nzCP0lyjPJP6IoNIFsVBam/A2jAvE3YmEN2nHHwb/9G9x9\nN/zgB/DEE5dx1FF9k9heQ9HS0hL7HL9SjwuzkKR3H5TVevuVjdn8jWFDKfb7SrH5+2UIpkDcuVMZ\nNSq/XjoNo0bBhAnR7qvpdJqmpqaiFsK1tbXR2tqa9fyZNm1aji9/kvDMSfMIG+XB4vkbA0Kl+MUX\nYjCVfzmRGtQoAAAgAElEQVS/SZTy37IF3vKWm7nggs/zgx+UJmuhSfFgo+C9Q5Q7LJAoMUvSPMJG\neTDlbwwICxYsGFK/+KQMpvIv5Zv4G4x0Op21x/t7005SlH24/vqdXHIJvPyyE23zT3+CY48tXtaw\nhOn+Xvy6deuy+X39jVZU4xb0OPIT1fP3U2kdhuHCgCl/ERkNPAZ0q+onRGRv4FbgINzk7aq6za07\nH5gF7AQuVtV73fIjcZK3j8VJ3h5qjDTlX3n4FUElD+GTKPIwZehfYOUlO4f490v6Tfz12tvbI5Oq\neMya9RMWLXoPcGRO+QUXbOOHP5wYaXJJSn975ElTMlaCOW0kMZDK/0s4f417quqpInIV8JKqXiUi\nc4GJqjpPRKYDtwBHAfsD9wHTVFVFZAVwoaquEJG7gGtUdVnIs0z5VzD+f+pKMwclUTiFzCB+yqEQ\no0YI0BcTZ/z48bz22msALFz4dm67bR9frRQwB8j7Vykpaml/e+SFlL/Z+YeGUpR/khW3DThK/ATg\nDrdsLTDJ3a8H1rr784G5vmuXAccCk4FnfOVnAddHPK/gajZj6PD/fpKsSh1oil2tW2hVcJJ7BIn7\nm42Sz39N7urYvRVeUnhN4WsKY4tewVyqrKVcH3y/of57GKlQwgrfJGPG7wNfAWp9ZZNUdbO7vxmY\n5O5PAR721evGGQH0uvseG9xyo4ppbW1NFHdmIAlOTCZNDBK36KmcawWCvd64e3srZX/3u40ccMDr\nbNnyFi688I3EGbGinjmQRE0Mt7S0ZGP2W6+/MolV/iLyceBFVX1SRJrD6qiqikhZ7TT+f9rm5maa\nm0MfbQwxxSi2gcLfAPmpZIXzwgsAv+Dhh+HYY/vKPa+bL3/5XYnuE5WYvdwT8UHznn89QNz3r+Tf\nQbXT0dFBR0dHv+5RqOd/HHCqiHwUZ6K2VkR+DmwWkXpV3SQikwEvIegG4ADf9Q04Pf4N7r6/fEPU\nQyvJg8SobJJOynpERZ1Men2hZ8azB5deCt/9LsC/cOGF8Mgj4TWLzYg1kCOvQqMrU/KDT7BTHDWv\nFEes8lfVS4BLAETkeODLqvo5d8L3HOBK9+dv3UtuB24RkatxzDrTgBXu6CAtIscAK4DPAdcULa1h\nFEGcr3qSRCrlCv28cyfcfDPAs3zrW33ljz8O998P06blX1NKRqyokVdczz2J4q7G0ZWRgKSTA8Dx\nwO3u/t44k8DPAvcCdb56lwDrcSaFT/aVHwk85Z67JuY5AzQlYpRCb29vzmQjvknLVCqlvb29qtr/\nicSBIGyy1T+xG7XFXR82EVzo3V95RXW//VTBvz2mDz7onPdPQvsnc4sJkxwnQyVMzBsDCxbS2Sg3\nScMmV8JK2Tg8+cL8/IORQcPeI6l/exQ33ADnnw9TpsBll8G5545CtS8IW/Abp1KpbD5diM6IFRZa\nIUilueQa5cdCOhsDxkAlSR/srFphkUCD+D1owhZSBRsQ/zWvvJJmjz3GU1NTk5P/9qijaoGfsnz5\nfzJ+vNP5zmQy1NTUhKY9LGdGLFPyRhjW8zdi8XqlUck8ytnzH8iFQVHyRa149Qh7v/DR0HigBTgN\neD/OAvd4wnrxpSaosRW1Ixvr+RuDhpfnNSwkMCQPPVAJ7qJQ2sims7OTSZMa+PWvJ7BwYR0vvui8\n71e+8gxz5ozJyXQF+Z47/mNPqXsyRMXztxj5Rrkw5W+UhGfO8Cswf4+1lNADfgbbTl2KOeW55w5m\n9ux6nn02t/yhhxq46qpxOff2KNTIeA1moRj5caYn75kWT9+Iw/46jEQEFc3GjRuB4TMXECRuZOOV\nbd8+KkfxT5qUYfPm87n55v8Gpobet5RGJozu7u4881A5G19j+GPK30hE0kxO5SIYOsKf0DyJn3p/\ne8aFRjYAH/3oa9x0E6xfD/Pnw6mnbmD69J9RU/PfSV6xLAxU42sMf0z5G7FEJUkfaCWTJKF53Cig\n2J5x3Mimp+etNDb2MmGC5rz33//ezZVXjmGffXZSV7cr0kOnmFXFxTJQja8x/DHlb8RSKUnSo0IF\nFyKqZxxUyOEN2eF84xtHsHz5eL79bbjkEhJck4+/XpKVxYYxGJjyN/pFf5RZMaaZYsMYePeN6hmH\nKe6WlhYWLlzIPfd0c/nlu9PRsS/LlzvnvvtduOCCvrpLly5l8uTJ1NbW0tTUlJdkxSPov1/KaMn/\nTbx9L/mMYZSKKX+jX/TH9FPuScu4WD5BonLXLly4jFNO2Z+gy/yHPgT/+Eff8cyZMwHH9AROg+Q1\nZmENYH19PRA+ggGYMGFCqJL36gb3W1qcRHg2kjBKxRZ5DRLDZYl90jAJSVwNvQVNUaaZsAVWScMY\neEqy0OK0sPd7//tf449/dFw1Z8zYyrx5bzB9em+ObJ7MTs7diahqoixhUe8Rl/ErLPNWOp2mqakp\n8jnm7TOysEVeFcxAui4ORcNSzrmAck1ahr1vXM+4txfGjMm/zx//+E/AfwNfo7PzCVxHo8Qyl+KB\n4/duCiZ39/C/XyaTCZ2I98tnGLEUGwluoDcqMDpkOSg23WAxDHTUxiSyl/J7898vrtx/76TPCfse\nfduH9NhjX9MLLgi/LpVK6bp162IjbXrRTP0yFXqfpO9Ryu9zuP7fGMnAonpWB+WOwzLQ8XGiTBL+\nUUsp7xRlglm/fj3Tpk3L9qCDIQ56e3sLmpSCicobG9/K5Zf/mZtvPoBnntkTgN12U1avTnPooXvl\nXOe9R5J9/3GhOEgQbfbxU8pIzmL7jGzM7DNCGej4OEmTeZSaTzZomnnyySeB6AVW3d3did7Rq9Pb\nC/AE8+dPzzm/Y8dOvvrVe/npT08KTXYyECS5d7XNAxnViSl/oyBJlVHYZKc/kXfUvaLs4eVaverY\n9VcB73ZL3uCzn93BF7+Y5p3vPCkyTWGpxM0z9PfehlEuzOwzBAzkEH2ohv9+M0tjY2PWZz6I31QU\n5jHkp9iEJuCkTBw9Osw0czi1tav5whfge9+bjOrG7LkoM4s/sUvQ7NPT05NnaivUOPnljkoaUypm\n9hnZlGL2MeVfAv31rhmuyr+QcoyTz3OZ9CtRf2MSRSqV4qCDpvLHP2a47LI32bULrr9+S2iGrn/8\nQxk/Pvk38ru1jhkzJlHj5h/pxLm/lvv3ZMp/ZGM2/0FiqCNOVgNR8xBxk8f+ekHCTUAf4qc/3cEd\nd7zJ6tW74/05NzYe4f7MbTTGjy/+Pbq6urIjD//9/IrfL1vc/Et/E6kHKff9jBFGnCsQMBZ4BFgJ\nrAEud8v3BpYTnsB9PrAOJ4H7h33lXgL3dcDCmGf2z+dpEOiv22Y53zFpgvWBIOo7BN/Pfxz37Qi4\nTAbdLP08/3xK4WnNTYrubF/72kvZ+/oTtkfJFwUBV8vOzs6cbxslW9h7q5bfJdcSsxseDISrp4iM\nV9XXRKQGeAj4MnAq8JKqXiUic4GJqjpPRKYDtwBHAfsD9wHTVFVFZAVwoaquEJG7gGtUdVnI87SQ\nTOVkMEw45VgNG0Wpaf/KQdyqVP83iUuhGFYvKqF5vvvkrcBcAHbbbRc7dtzE3XefzCmn7B+ol/99\nkvyNefMYwRXHwXvHrRb2P6fci/GGy6pxo/8MiNlHVV9zd3cDRgM9OMr/eLd8MdABzMNJYLpEVXuB\nLhFZDxwjIn8F9lTVFe41NwOnA3nKf7AZDBNOoRg2cWaDpAxFXPcoF9D+erT4Uxk+++xGTj55dkTN\nHwEzufzyRj784W6OPHI2hx2WCq3pBWIDx7soLIBcWKC5clJupWxK3ugPBZW/iIwCngDeBlynqk+L\nyCRV3exW2QxMcvenAA/7Lu/GGQH0uvseG9zyISeYNKSYcMHFEqWg/Uq61EZnKOK6D5TyGT26BjiK\nb31rKkuWHAgs4YUXwtwn/wa8jXnzdtHVtSv2nl4gNo+wAHKFGuko0un0gK4NMIyBIEnPfxfwbhHZ\nC7hHRE4InFcRKaudxq/8mpubaW5uTnRd0mFwVD0Y2Fj1cQq62EbHe4dCPVX/+YFS1oVCM8eR+7uo\n4Rvf2MqSJXsCK/jZzwBGAU188INHA4+G3CH3Ty/Oxz7p6ChY79FHH2XmzJmR9160aFGo5w/YJKwx\nMHR0dNDR0dG/mxQzQQBcimPzXwvUu2WTgbXu/jxgnq/+MuAYoB54xld+NnB9xDNKnvRIOgEWVa+U\nZye5ppiYL0kJvkPcvaO+Q7kIe1aS7+rJlFt3bcgk7lN6003doRPa3r2TyFAohlDU78k/aRy2rVq1\nKkeeJH+DhlFOKGHCN7bnLyL7AhlV3SYi44CTgK8DtwPnAFe6P3/rXnI7cIuIXI1j1pkGrFBVFZG0\niBwDrAA+B1wT9+xSSGrCiapXTasvvXeIykwV7PFC8vy3peL1mDOZDJs2bWLjxo1Zc8tDDz1EJpNh\ny5Yt7Lfffq5pB8444wxOPPFEX92fAN8FXgNudY//xPHH902qho2eotJNQv/nP7xY/GGxhrxn+yfs\n/ec8rNdvVBqFzD6TgcWu3X8U8HNV/YOIPAksFZHZQBcwE0BV14jIUhy30Awwx22VAOYANwHjgLs0\nxNOnvySNcTPQsXCiSJJ4I6npKkmO27hz/sVISc1ihRoMz6zV1dWV90zneA/gDJz+wv8DyItJf/vt\nn2blyq2cdto/qK09ge7uaYkU90Cmm/QUu99sF3dvi6NvVAXFDhUGeqNMPvBJ7+OvV8qzk1xTjFmk\nWN/t3t7eHFOIfz/M/z7JvYuVIcp80tnZqevXp/QXv9ion/zkdt19996sKWfKlF597rlwX/y4ewe/\nd9T395f75Ynz009aL+53Xq6/X8MoBspt9jHyKSbvrL8syiThHXsU632UJMet/9hv4irWLFas6cIx\nAU3ls5/NP7dpUw3p9FTe/e78c4UITiYX+v4epSRcN4zhiin/Iikl72yYSSJIksTl5SCukfAo1Szm\nuWP6G8eDD4b3vQ/+9Cev5Gnmzp3MRRftzf5FOPvGBYAr9P0LNb7+IHFJ6hnGcGBEK/9MJgPkKt6k\nvcj+LqoqZ+LyoeSZZ+CGG+qANXzwg/8K/F9enQsugKOOgsMPf5J//dcjOP/8FPvvv3fo/aLmRYLf\ntpjvn3Q+YCDnDQyj0hjRyj8sYFeUEg761ZdjUdVQrMotB2vXwi9+Ab/9LTz9NEAdUMc559zNggUv\nZ7+h963e/35ne/TR9TnlHv7jQrH9vXsPxaK2MEoxAxpGJTDi/ir9/6zezyRKOCoMRJCo1Z5hnjKV\nosCK5YEH4Nvfzi//4x9rOeig2uxxlCJPkrzFb3JJp9OMTxiS06+Iy+nWGhVBM51O53ksVesozhhZ\nDBvln7QHFmazT6KEg371UUSt9ixHrKC4lclhpqv+KLzt22HduvBzp50Gc+Y4+2PHwsc+Bmee6fwU\nN7RUmO08k8kwbdo01q1bR01NTUF/+UJhoMMIS/1Yjm8fle2rpaUFqN5RnDFyGTbKv9iJWG/hU9J/\n0GDvMcp0MWvWrNBkHkkUsTcHEdZ7Badh+eY3vxl6bZjpylNM/ntkMplQM4SqY7+/6y64+27o7ITa\nWnBi+fWxbds2duzYRktLHYcdtoMPfvB1xo9X6urqGD8+epLYf3zwwQdH1gvD731UaFFbUAlDeRZY\nRQWxS6fTLFy4sGpHccbIZdgof4+kPbD+em5ENRq1tbWJPGo82fwUSly+atUqZs2albhX6Y1A/Pd4\n+OGHs+/uNQQ7d8L06fDss7nXb90K8P6csmKjoBaaVI9qjPwUk/d3oJRwlOnIAroZ1UpVK3+/GaTQ\nRGzQzh92LurYjz/cMCR3Bwwz2RSbuDzYsMQ1dOEjmxpmzPgQsAPIbQgOOuhAnn12VI4c06e/yZo1\ne+QouAsvvLCoNQCFJtW7u7uLUtalumNa1ivDCFDsqrCB3ihihWTYStQkAc68eoVW3gbvFVhNl7cf\nVy9q1WxLS0vZApD53yeVSulzz6UUDtcLLnhOjz12q0Ja4axQOb72tZe0tlb1Ix95VeE8hckFv4n/\n/aIydK1atSrRqtm4bxdH0nqlZL1Kcu+kq4INYyBhpK3wDa5EjSPoKtjd3U19fX22h+wFFvP3psu5\nqCfKZlxXV8fChQtzMkaVg1//egJXXDERWMUPf9hX/rGPfZ8f/ODy7HFfToG/8NnPTmHz5m7uvvvG\noicwC02IDrVNPO77+yl1hGCTu0a1UdXKP+wfMsqEE1Q+lWIz9ijlma++Kjz/PLz1rfnn9thjFy+/\nPDqvfMuWOqZOHZtXftJJubb9Yr9DoQnRKLzfj/dtBso0k/QeUY1Y1JyGrQo2qpWqVv5hJOmBRf2z\nNjY2VsQ/a3gDNpo1a8Zw331w//37AI/T1HQgRx/9OkuWbM67bsqUFPAW4EU+/vE9eN/7XuerXz2G\nX/3qD8DU0Od2dnbmhGEuhlInRL3fV1tbG0BRincgSDpC8LBVwUa1MqyUf9IeWNw/a5jnSaE1BHEe\nK0kXfPkJb8AO5WMf84Lh7Akcwa5d8PDDSmPjITiZMvv45CePBQ4D/sLcuf/PLX0+64ef/LnloZBr\npvc9ilG8A4FNABsjhWGl/ItxqyyGQmsIPI+VsEYibC7C35N1rBwf4Ic/hJUrDwTu57nnGhk1yrm2\npaXFNZs8C7wC7JW9j4hy8MG7WLduMuvW/SGbhBzILmgL5giO82zyzzsU4wEVRjAcRlIzmylewxgc\nxJkorhxERIMyJUkuIiL4r5OYLm6wnnccvIdHV1cXjY2NkZOg3gIyr14U/l7uXnvV8Y53OAurgqxf\nD297myNPT09PVonOmNGFkwvnUXd7DEgD4WEEMpkMY8aMSTyZrKqJ3wHCR1feCChuVW7cqt5iifqd\nDRWVJo8xMnD/7mLG9flURc+/2IVFHmG++P1J1VhoDcHGja8C7+Kyy35FOl3P88+P4fnnx/ClL/2Z\nf/7no3Ku37ZtG7vvPg7YPe9+q1c7yh+CZogZrvI9Hzg/++y4aJbgjIgKrU/wB0yLqhf2rKhV1FG2\n88bGxn7bxM1n3zD6T1Uo/6BLZ2dnJ5lMhtdff52HHnooWy9oX/dMMUFKW5UpvPTSKF55BQ48ECZO\n7DvjKcQDDlgLrOaSS3Kv3L79gLy7tbe3s3LlFJyUhmuA1dnthBOWR0pRqjdSkqQvSeslcQMdSEVc\nrEfOQGONkVGNFFT+InIAcDOO64gCP1bVa0Rkb5wM2wfh5vFV1W3uNfOBWcBO4GJVvdctPxInj+9Y\nnDy+LUmETJqvNqxXH7/CdCwwhdWr/0Zt7S4g9x/3Jz+p49e/hu7uBuBNjjpqDABLlsBZZ/U9w1OG\nX/pSLS+8kC/X88+PyStrbW3ljDPSjBnzAlu3/gN4G/A2Zsy4mm3buvB0SZLwB4NNtfjsDxaV1hgZ\nRhKSaJVe4D9UdaWITAAeF5HlwHnAclW9SkTmAvOAeSIyHTgTmA7sD9wnItNcQ/51wGxVXSEid4nI\nKRqSyP2222DUqL7tkEPgsMP6zvfZsOfw1a/exp577k0mI/T2wksvpfnhD88GHsm5Z2dnJ8uXv51b\nb92TV14ZRU+P8Oabjh98U9NFwA+A3H/cnp4FPPpo/md68skejj32lbw1BIce+jLwFz7wgQM5/PBx\nHHqoI/fee7/KlVfmvqPXoIUlO+9P+IORQKX1qCutMTKMJBRU/qq6Cdjk7r8qIs/gKPVTgePdaouB\nDpwG4DRgiar2Al0ish44RkT+Cuypqivca24GTgfylP9nPpN7/NWvwre+1Xfcpww/xbe/fUjg6n2A\nEwkq/4aGBnbt2oe//CX/Hc877yvMmnUmtbW11DqhLKmrq+P664M1e4C/c9VVP+Gqq/IXLn3xi2mu\nuOIwLr883yxSiKQravvrheMnLhduJY44KpVKa4wMIwlF/XeLyFTgPTiadZKqbnZPbQYmuftTgId9\nl3XjNBa97r7HBre8IKNGRZ3pjSifADjKzK8cJ0wI1tsBbOXGG2/gxhsvyxumn302zJgB++2XYdeu\nDYwdq8AeNDYuzEk0kslkcp5Vir98ocnkdDpd8r2jSOK+GkY5GyDDMIaGxMrfNfn8BmhR1e1+V0pV\nVREpm3/bGWfArl1926GHRtV8itNPf5WJEydQUwNjxsAbb7zCokXOJHBQuX3+83Dyyc5k7auv/pV3\nvnOq2+P+AvCFvN7bXnttw53GyKNQopGLLrqIa6+9Nie6ptdIQHIXR0/ZX3rppUCuF07YvePCJITR\n2dlJfX09mzZtAvpiHPkVenAUYHFsDKP6SaT8RWQMjuL/uar+1i3eLCL1qrpJRCYDL7rlGwC/e0sD\nTo9/g7vvL98Q9rx3vWtBdr+5uZnm5uYIyW7nP/9zZZ65ZNGi32fNKP5FSw0NUFfnJC7p6SmcizdJ\n6sagJ5J/RHDttdfmKMpp06Zl95Om91u6dCmTJ0+mtrY2L5FL2L3jwiRAfkPgfbugQg8uDPPcRS2c\nsmEMPR0dHXR0dPTrHgUXeYnTxV8MbFXV//CVX+WWXSki84A6VfUmfG8Bjsad8AUOdkcHjwAXAyuA\n3wPXBCd8g4u8gqtmPYVTaNFS0oVXYUnaPdLpNOl0mtraWpqamnJ63WHfzb/Ap9iFYaWkAYzLWuW9\nR3t7e2xgteBIIk7WuPeNI2p0ZN4whlEeBmqR1/uBzwKrReRJt2w+cAWwVERm47p6AqjqGhFZiuO8\nngHm+LT5HBxXz3E4rp55k71BCoVW8Idq9hRZJpPJmkE8P3+/icTLIQu5vdW4xWRQWsCupG6RcYo+\nSinH3dtrBBYsWJBNK1moISg071BqT928YQyj8qj48A7+nnGYbdqv/Ht7e6mpqSmqt+/H3/MPZqma\nOHFiwTAQYT3/4LOC5VEjG4jueUfduxDBzGczZszIaRSjnuMR7KlbKAPDqAyGbXgHKM427VGsKSXY\nsx0M//rBDAlcTC5cP8HIm4ZhVD9Vo/z9FJOkfagXSFWqW6R/8rZQJM9K+I6GYZSXqlT+Q6WM/J4y\n3nzCa6+9ll0YBrB69WrS6TQbN24EyuMWORANSNiIw1w4DWPkUDXKvxJ6zOHxgXJpamqKvL7U9H5R\nSrm/E7FBWaIifvplNbdNwxgeVI3yT9IrDa6G7W+POegT702OAtnJZ2/SNMwMFYx9X+xoJcyv3o8/\nBWJ/XCaTRvwEC2JmGMOFivf28XvDhClVf7lHX/arcPwTmFG91bhkMJDfSy7kjZPUMyaunojkNAYe\n5eh1J0lqA8kS6xiGMbgMS2+fYmzTnlKfMGFC1rcdkvWYw1Iw+u/pkWTx1UAy1BOvpuQNY3hQ8crf\nTxLbtGfC2HfffXOujeoxe4QtJvPuOdAK1+zohmEMNlWl/IuxTQdJWi8YEygJ/Z1bGCo7elSjYxjG\n8KeqlP9gUEpPP6qRSKfTee6hkN+jH6rwB1GNjmEYwx9T/gWI69UXinK5aNGinEicUT36oTLv+Bsd\nL4gdOI2ZmZ4MY3hT8d4+IecTeaVEXRNF0EOnlPhAweeU2zNmIGPpWORNw6hehqW3Dwyubdrr2Wcy\nGTo7O7MB5KD4RVrl6DUP1mSwRd40jJFFVfT8o3qlQNl7/nEkyH1Q9p659cgNwyhEKT3/qlD+/t5v\n0DbtD7vs76XGhWcO683GhVb2jodC+duiKsMwClGK8kdVK2pzRIqmra1Ngbytra2tpHpRBOUoJFfS\nOoZhGOXG1T1F6dqq6Pn7SdoT7m+POdiLT9Krt+QmhmEMBcPW7DNEciRS/sWalwzDMMpNKcp/1EAJ\nM1Job2+nsbExO1ns7be3tw+xZIZhGNEU7PmLyCLgY8CLqvout2xv4FbgINzk7aq6zT03H5gF7AQu\nVtV73fIjcZK3j8VJ3t4S8byq7fn7sZ6/YRiDxUD1/G8ETgmUzQOWq+ohwB/cY0RkOnAmMN295kfS\nFxv5OmC2qk4DpolI8J5DzrZt2+jq6sr60q9evZqHHnqI1atXA2TP3Xnnndlr6urqmDp1at5WCYq/\no6NjqEVIhMlZXkzO8lENMpZKQeWvqp1AT6D4VGCxu78YON3dPw1Yoqq9qtoFrAeOEZHJwJ6qusKt\nd7PvmoohaMJpampixowZ2exc3rnvfve7QylmYqrlD9fkLC8mZ/moBhlLpdQVvpNUdbO7vxmY5O5P\nAR721esG9gd63X2PDW55RRFc5eqtKaitrc3J03v99dcPgXSGYRjlo9/hHVRVRWTojfRlIKmdfuzY\nsYMgjWEYxgCSZDEAMBV4yne8Fqh39ycDa939ecA8X71lwDFAPfCMr/xs4PqIZ4UuzrLNNttssy16\nK3aRV6k9/9uBc4Ar3Z+/9ZXfIiJX45h1pgEr3NFBWkSOAVYAnwOuCbtxsTPWhmEYRvEUVP4isgQ4\nHthXRF4AvgZcASwVkdm4rp4AqrpGRJYCa4AMMMfntzkHx9VzHI6r57LyvophGIaRlIpb4WsYhmEM\nPBWzwldEThGRtSKyTkTmDrU8HiKySEQ2i8hTvrK9RWS5iDwrIveKyJA79YvIASLygIg8LSJ/FpGL\nK01WERkrIo+IyEoRWSMil1eajH5EZLSIPCkid7jHFSeniHSJyGpXzhUVLGediNwmIs+4v/tjKk1O\nETnU/Y7e9oqIXFxpcrqyznf/158SkVtEZPdi5awI5S8io4Ef4CwMmw6cLSJvH1qpstxIwkVuQ0wv\n8B+q+g7gWOAC9xtWjKyq+gZwgqq+GzgcOEFEPlBJMgZowTFhesPjSpRTgWZVfY+qHu2WVaKcC3HM\nvW/H+d2vpcLkVNW/uN/xPcCRwGvA/1JhcorIVOALwBFu1IXRwFkUK2exM8QDsQHvA5b5jnO8hoZ6\nI9zbaZK7X4/r7VRJG84k/ImVKiswHngUeEclygg0APcBJwB3VOrvHUgB+wTKKkpOYC/g+ZDyipIz\nINuHgc5KlBPYG/gLMBFn3vYO4KRi5ayInj+OZ9ALvmNvcVilErXIrSJwewbvAR6hwmQVkVEistKV\n5anJFnsAAAKESURBVAFVfZoKk9Hl+8BXgF2+skqUU4H7ROQxEfmCW1ZpcjYCW0TkRhF5QkR+IiJ7\nUHly+jkLWOLuV5Scqvoy8D3gb8DfgW2qupwi5awU5V+1s87qNLMVI7+ITAB+A7So6nb/uUqQVVV3\nqWP2aQA+KCInBM4PuYwi8nGcQIZPAqGux5Ugp8v71TFTfATH1DfDf7JC5KwBjgB+pKpHAP8gYJKo\nEDkBEJHdgE8Avw6eqwQ5ReRtQCuORWIKMEFEPuuvk0TOSlH+G4ADfMcHkBsOotLYLCL1AG7coheH\nWB4ARGQMjuL/uap6ay8qUlZVfQX4PY5ttdJkPA44VURSOL2/fxKRn1N5cqKqG92fW3Ds00dTeXJ2\nA92q+qh7fBtOY7CpwuT0+AjwuPtNofK+53uB/1PVraqaAf4Hx3Re1PesFOX/GE6kz6luq3smzoKx\nSsVb5Aa5i9yGDBER4GfAGlX1JxOoGFlFZF/PA0FExuHYKZ+kgmQEUNVLVPUAVW3EGf7fr6qfo8Lk\nFJHxIrKnu78Hjp36KSpMTlXdBLwgIoe4RScCT+PYqitGTh9n02fygQr7nji2/WNFZJz7f38ijmNC\ncd9zqCdWfJMYH8GZxFgPzB9qeXxyLcGxq+3AmZc4D2fC5T7gWeBeoK4C5PwAjn16JY5CfRLHS6li\nZAXeBTzhyrga+IpbXjEyhsh8PHB7JcqJY0tf6W5/9v5vKk1OV6YmnAn+VTg91b0qVM49gJdwohB7\nZZUo53/hNKBP4URWHlOsnLbIyzAMYwRSKWYfwzAMYxAx5W8YhjECMeVvGIYxAjHlbxiGMQIx5W8Y\nhjECMeVvGIYxAjHlbxiGMQIx5W8YhjEC+f9paqWV54VkFwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 392 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Not too bad, but we should also look at the residuals, i.e. the difference between the data and the fit line. Structure in this can hint that the fit is bad, even when the reduced chi-squared value is small.\n", "\n", "To highlight this, let's look at what happens when we try to fit a curve of the form $\\sqrt{a^2+bx^2}$ to the data:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def lin(x,a,b):\n", " return np.sqrt(a**2 + b*x**2)\n", " #return a* x + b\n", "p_opt_lin, p_cov_lin = curve_fit(lin,x,y,sigma=y_er, absolute_sigma=True)#, p0=[15,1])\n", "y_fit_lin = lin(x_fit,*p_opt_lin)\n", "ax.plot(x_fit,y_fit_lin,'r--',lw=2.5)\n", "display(fig)\n", "\n", "# calculate reduced chi-squared for both quadratic and linear fits:\n", "chisq_quad = (y-quaddy(x,*p_opt))**2 / y_er**2 #as array\n", "chisq_quad_reduced = chisq_quad.sum()/(len(chisq_quad)-3)\n", "print 'Reduced chi-squared (quadratic fit):', chisq_quad_reduced\n", "\n", "chisq_lin = (y-lin(x,*p_opt_lin))**2 / y_er**2 #as array\n", "chisq_lin_reduced = chisq_lin.sum()/(len(chisq_lin)-2)\n", "print 'Reduced chi-squared (linear fit):', chisq_lin_reduced\n", "\n", "\n" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXt4FNX5+D9vCIhcYgCRAAGJiiJWqVrFS9FYrVK1XlvU\n/qwXrNV6C1otFy9BbUHxQrBWra0oWG9o+7VqFcVqbGxFvBEviIJulCAgYkIARbLk/P6Ymc3s7Mzu\n7GZDNsn7eZ55mDlzzsw7s+Q9Z97znvcVYwyKoihK5yKvrQVQFEVRtj2q/BVFUTohqvwVRVE6Iar8\nFUVROiGq/BVFUTohqvwVRVE6ISmVv4iUich7IvK+iJTZZX1FZIGIfCwiL4hIoav+ZBFZJiJLReRo\nV/n+9nWWicis1nkcRVEUJQxJlb+IfA/4FXAAMAo4XkR2BSYBC4wxuwP/to8RkZHAacBIYCxwl4iI\nfbm7gfOMMcOB4SIythWeR1EURQlBqpH/COB1Y8xmY8xW4BXgVOAEYI5dZw5wkr1/IvCIMabRGFMD\nLAdGi8hAoLcxZpFdb66rjaIoirKNSaX83wfG2GaeHsCxQDEwwBizxq6zBhhg7w8Cal3ta4HBPuUr\n7XJFURSlDchPdtIYs1REbgZeADYBi4GtnjpGRDRGhKIoSjsiqfIHMMbMBmYDiMgfsEbwa0SkyBiz\n2jbpfGlXXwkMcTUvtuuvtPfd5Sv97qcdiaIoSvoYYyR1rWbCePvsZP87FDgFeBh4CjjbrnI28KS9\n/xRwuoh0E5ESYDiwyBizGmgQkdH2BPAvXW38HiKnt/Ly8jaXQeVUOVXOjiNjeXm5ry4sLy+nrq6O\nSCSSsNXV1cXaZ0LKkT/whIj0AxqBi4wx60XkJmCeiJwH1ADjbKW9RETmAUuAqF3fkewi4AFge+BZ\nY8z8jCRWFEXpYEyYMIFzzjkHgJKSEiKRCACFhYVUVFRw/fXXJ7QpLy9n6tSpGd8zjNnnMJ+yr4Gj\nAupPA6b5lL8F7J2BjIqiKB2awsJCCgtjy6UYNmxYbD9Zx9ASwoz8FQ+lpaVtLUIoVM7sonJml/Yg\nZy7ImKxjaAmSqb2otRARk2syKYqibCtEJNCOH3TOLs/uhK+iKIrS8VDlryiK0glR5a8oitIJUeWv\nKIrSCVFvH0VRlDYkGo1SW1sbV1ZTUxPbLy4uJj8/+6palb+iKEobUltbS0lJSVyZ+7iqqori4ubo\nONnqGNTVU1EUpQ2pqamhpKQkQcnX1tYyZsyYpG0jkQjDhg3LyNVTR/6Koig5QHFxceACrkw6hlSo\n8lcURclxgjoG71xBOqjyVxRFaae0ZPSvyl9RFCUH8I7iw4zqHXOQd8I4DKr8FUVRcoBko/igjiHZ\nPEEq1NtHURRlG1FfX099fX1cWTQa5ZtvvqGgoACID9ucamK3Jd4+qvwVRVG2EVOnTk2ZmMUdudO7\nAMzdMUCzn78qf0VRlBzGPfL3JmZxYvaHDenc1AQXXAA/+Qmceqoqf0VRlHZBitj8SdsYA7/5Dfz5\nz9ClC2zdqspfURSlXdAS5b9wIRxyCDRXa4VkLiIyWUQ+EJH3RORhEdlORPqKyAIR+VhEXhCRQk/9\nZSKyVESOdpXvb19jmYjMSkdIRVEUpZkf/CDK7bevpUsXw4knbszoGklH/iIyDHgJ2NMY852IPAY8\nC+wFfGWMmSEiE4E+xphJIjISeBg4ABgMvAgMN8YYEVkEXGKMWSQizwJ3GGPm+9xTR/6KorRb/Dx6\nIDEXb0tG/k48IDgEeB3YmvWRfwPQCPQQkXygB/AFcAIwx64zBzjJ3j8ReMQY02iMqQGWA6NFZCDQ\n2xizyK4319VGURSlw1BRUUFJSUnCVlFRkfV7VVXdTCSyPKO2SRd5GWO+FpHbgM+Bb4HnjTELRGSA\nMWaNXW0NMMDeHwQsdF2iFusLoNHed1hplyuKonQoJkyYwDnnnAMkevQE4f1acMI2O20++aSBfv2a\nYueyscgrqfIXkV2BCcAwYD3wuIic6a5jm3TUTqMoikKiecfhq6++4v333weILehylHxFRQWzZjVP\nhTrhGsrLy3nrrf155pnDgGPjzrWUVOEdfgD8zxizDkBE/gEcDKwWkSJjzGrbpPOlXX8lMMTVvhhr\nxL/S3neXrwy6qbPYAaC0tJTS0tIwz6IoipJzJFPWyZK2ADzzzE4880wPAHr0qOKhh1bz/e9v4emn\nn+ayyy6joqIi6RdFMlJN+I4CHsKawN0MPAAsAnYG1hljbhaRSUChZ8L3QJonfHezvw5eBy6z2/8L\nnfBVFKWDIyIxs4+j6INi8zuhGhzmzIFzz2125/zBD2DBAigsTEwAU1JSkt1kLsaYahGZC7wJNAFv\nA/cCvYF5InIeUAOMs+svEZF5wBIgClzk0uQXYXUe2wPP+il+RVGU9kBYjx4gwSYfxk7/t7/FK/59\n94UXXrAUv5tWDelsjJkBzPAUfw0cFVB/GjDNp/wtYO8MZFQURckpKioqUsboaQnbbQd5ebB1K4wa\nZY34+/RpPl9cXBwX4yeTeQBd4asoihJA0Ajf7YsfFKPHW0/Essp4zTuOCcdb/ve/w/TpMH8+7Lhj\ncjk1h6+iKEoWCTvCd5R2NBqNefA4eI/DJm059VQ46SQrdk9roCN/RVGUANKNwtm88jZ9vCP/VPK4\nyfqEr6IoSmfGa8YJUs6OUnZG8UEePe5ypzN56aXt6dGjieLi/inlCfoSyYSUgd0URVGU5DghHRzv\nm+Li4gSffT+WLBnGb34zgPPOG0hVVeqx+IQJE4hEIrEvEPd+uqjZR1EUJQTJArHV1dXFRv6Ozz6k\n8sI5nq5dn6ax0Tr63vdg8eLwNn7vZLKafRRFUXxIxzc/XZJdw98EdBswL6b4hwxp5O67V7NixdbA\nTsbvHt7J5HRQ5a8oSqcgW7753ry60KyEnXL3ee+irnXr8oAHga52yXJWrDiCMWOsNocffjivvPJK\nKDlbEudHzT6KonQKwnjuJCMxln44/P36LwX+zrBhwsMPr+aQQ4bG5ElnDYHbvKRmH0VRFB/Ceu44\nBI3wU3n0gKXw3ceJPAOcwH//O59Bg4YklSeZnJmGcwZV/oqiKL7U1tYmjPDdx8li9LjLgxd1Pc+g\nQdmQNDNU+SuKoiTBO8J/4403GDduXOj2Y8b8HFjte64lE7YxtmzJqJn6+SuKoiTBGeE728CBA0O3\ng/Pp3n0lDzywCoCysrK4Ok6KxzA4oSOcDWD9eeex8fTTQz+LGx35K4qiZIBjvolGo6xevZpVq1bF\nztXU1PDgg72Be9m8GS68sAg4lKlTpzJhwoSEa4XpAPzMUGfMns2zGcqvyl9RlA5PMvfMhoYGevTo\nQX5+vDpsaGhIes1ksfRLSv4I3BY7/t734M03P8jKmgKvGerT226DO+9M+zqq/BVF6fCkmrz1wzHR\neDuNaDQKEJdFy3181107cMstzcH3R45s4L776hk1KnGBWWjefhvZbjsgcaK55re/VeWvKIqSjFQB\n19x+9Q0NDcyaNStwhO9Vws7xMcfAHXfAd98BvMKSJT9l1KgNcW2DvkQcE1JRURH5+fn0Bzacfjq9\n5s1DfvnLbLyCGKr8FUXpNAS5Z7rL3bH5vdmy3MdBgduOPBKeeAJmzvyGl176CVVVL8R9IUDqL5F8\n4BLgY6D3Y48BUDh3Lr3DP2pKVPkriqL4kJ+fn9BRhF1UdfzxsNdeX7LLLt8mXQ/g50b6q3HjWLXz\nzvT47LNY+dcHH8xBr73GBsIng0lFSuUvInsAj7qKdgGuBf4GPAbsjJ3E3RhTb7eZDIwHtgKXGWNe\nsMv3x0ri3h0riXu835OiKEob4+d7HzRJ645E41XCK1cmxvlxTzJD4pdIbW0tDYAZMQI++wyGD4eK\nChpGjmSZ/WXQkqTtbsIkcP8I2BdARPKAlcD/AZOABcaYGSIy0T6eJCIjgdOAkcBg4EURGW4H7Lkb\nOM8Ys0hEnhWRscaY+Vl5EkVRFBfuWD5+Adf8jsF/Iri8vByI7xg2bxYuvbQ/cBLwZKBSdpc71/b6\n+3upu/pqeh55JJSVQbduYN/XLxlMkMypSNfscxSw3BizQkROAA63y+cAlVgdwInAI8aYRqBGRJYD\no0XkM6C3MWaR3WYu1ltT5a8oSgItDcHsF8UzzKjZrVDdQdWuv/56l5LtBfwTy/DxGA88UMfhh3+b\n0C6ZG+msWbPo9sEH4GMSig4ZAj6yer8UtmVsn9OBR+z9AcaYNfb+GmCAvT8IWOhqU4v1BdBo7zus\ntMsVRVESaGkI5gkTJnDOOecAlkKuqqoCoKCggFGjRgVO3iZTrpFIhIaGPEaNWgkcbJd2o7q6G4cf\n/m2sU0qllL946in+CxSddBK1zz9PdJddAGILxbJl109GaOUvIt2AnwITveeMMUZEshaH2f3DlpaW\nUlpamq1LK4rSTvAqb/coPAzeL4Qf/vCHcefDjJq9Xx9r13bhrLMGAENdtX7HzJm3MHNms3kokC+/\nhClTGDh7NoMAolFePvJIzvJUS/WFUllZCZBWHgIvoeP5i8iJwG+MMWPt46VAqTFmtYgMBF42xowQ\nkUkAxpib7HrzgXLgM7vOnnb5GcDhxpgLPffReP6K0glJZubp06ePb3arsHizYznHXn97rx392muv\n5cYbb3RdaX8sC3cvAK6/fh1nndXsw59U1ueeg9NPB3uy13TtSsO551J/ySUM22cfIpFIgp+/n3up\nY0JqaRpHjDGhNiyPn7NdxzOAifb+JOAme38ksBjoBpQAn9DcybwOjAYEeBYY63MfoyhK56O8vNwA\nCZtT3hK87Z3jSCTie09nq66uNpFIJG578MFVBjaaOXPC3SvG558bs/32xoAxxx5rzEcfpWyT7Lnd\n5+z90PrcGBNu5C8iPbFG7iXGmA12WV9gHtb3Tw3xrp5TsFw9o0CZMeZ5u9xx9dwey9XzMp97mTAy\nKYrSsUiWaSvsyD/o68HOdBU79mblClr5683C1dx+J4z50leGoBy8ANx7LxQXw7HHhmqT7FrbJIG7\nMWYTsKOn7Gss7x+/+tOAaT7lbwF7pyOgoiidg3QzbfkRNEmcimQLsfxZG3imF8CKFTBkSOLJX/86\nTcni8XZuLckHoPH8FUXpMEyYMIFIJBL7anDvp88ZPPFEr/DVm5pgzhw+BjjnnPgVYFmioqIiLgdA\nOvkAvGgCd0VRco6gCdpMryEi1NXVJZiV3Oad4uLi2ORvRcVWZs3albw8w913r+Xoo78JnGwF4PXX\n4bLLYNGi5rIFC+AoX+NIyucNKk9h1krL7KPKX1GUNiNdG306eJV/eXl5oEnIWQNguViWA1NdZ98E\nDiQS+TRmGoqT5+KL4a67YrVXAt1uv51NJ55IYd++oVxTM7H5+9TLvs1fURSlNcjURp8KJ+a+2yZ+\nyimncJQ9EndcKZ2Rv6X084A/YsXTtNh33838/vdRjjvOBCaDKejdm77AZqz0LdOBTVdcAVdckXRB\nWpD9PhsJX8KgI39FUdqMIA8fr397Mn93PxwvniAcLx63n/9bb63lZz/bBegHwDHHwN//DmvXJr9W\nN2Dlaaex+aKLiA4dGncumSKfOnWqb8c3ceJELrzQWv7k9XoKulYmI39V/oqi5AReM00yglwwHVK5\ncFZXV1NQUBDXxjr3O7bb7lVOOimPuXOdmGrWtXYD7v/Pfyh2efGkcglNRpDJ65577uHmm29OKE/2\nFaFmH0VROhRByjssQS6cs2fPZtasWb5t/vGPVRxzzGC6dLGO8+rruQP4DfB1dTU7ZSmkctBIftKk\nSbGRv7d+NlHlryhKm+Nno4dM/O/DMX78eCZMmAAkev4UFERYsaIRolF6P/wwA2+7jUvtdn1mzIDf\n/IZYz9AKbCubvyp/RVHaHMfuHtZn3anvVpRu+32qKJhLlvTnq68Gcok9t+vuYMaMGcNQ4BlgmKvN\nv4ERf/kLg1tR8W9LVPkripIzOGaeVOYd55zbDu6XF9c/NPLPOfvsIrZsgX79ms855qVIJAKNjXy7\n++4ANA4ZwrILL+Soq6+mqmdPGl1fJ60RanlbocpfUZScwS+toRvn2Okk/Mwjbp/9xA7kMuBRtmyx\n5kYnTwbYDiDmPeTc/yjgxWnT6Hr55fRYvRquvjprKRRzAVX+iqLkLEHKNtlcQHFxMcXFxTEXSesr\n4jBOP/1zHn20efJ4jz1g/rNNHLHrd77X+Tc4vUPc9cDf9bS9ocpfUZScw1G2bvOPew1AKmWbn5/v\n6Rz68+9/F8WODjkEnp26iB3+Xxn/Az57/31MLyuOj1+wtMTrtSyFYi6gyl9RlDbBb4I2mQ3drWyT\nLfDy50seeOBLTj99EKcesoq/9J9C/tEPADAQKN97b26wa2YaKK29ocpfUZQ2wW+CNhs29aB5gpEj\nt7Bk8oMMnnYRsnGjdbJbN27asoXx777Lub17x7Xr6J2AKn9FUdqUqqoqioqKWL16NWAlMR83blzc\nAq90FHGyDqT44CHgKP4TToDbbmPy8OFM2rvzpRlR5a8oSpviTN7utttuQLPNPewCLydMQjQapaqq\nihdf3JGamh7MmbNz4jzBsGE0XXklX+69N5sPOyx2DbedP1XcoI5Cx39CRVE6NPGRQa8Gfk8hdfTm\nl77zBJ9ffHHCl4T7OJM4Pe0RzeSlKEq7ZsKECXz0UYRTT91AHtdzAfewjOFMpQk7aoQvVVVVsUxf\nkUgktj6gtrY29iVQU1NDTU2NbwC29k6okb+IFAJ/BfbCymp/LrAMeAzYmcQE7pOxErhvBS4zxrxg\nlzsJ3LtjJXAvy+KzKIrSDgmaoHXwi/vj3u/Zs5gLLyyk6eVK3qaMUbwLwKU8xOrXLiBqR+F02jQ0\nNADBZiX3nIHzRZAsomZ7JazZZxaWsv6ZiOQDPbG+rxYYY2aIyERgEjBJREYCpwEjgcHAiyIy3I7T\nfDdwnjFmkYg8KyJjjTHzs/5UiqK0G4ImaBsaGqipqfGN++Pef+O1T/jtokkcx+OxMrPzzpz+2Wf8\nw2XXd9qUlSUfc3ojiYIVQ6itk69km5RmHxHZARhjjJkNYIyJGmPWAycAc+xqc4CT7P0TgUeMMY3G\nmBpgOTBaRAYCvY0xTpLLua42iqJ0MpyFXEEJ15944glKSkpinUOQmWbHojx++KNuAGyiK+aGG5EP\nP+Qfruu5t/Hjx6eUa9iwYXFbYWFhYPL0ioqK1ng9rU6YkX8JsFZE7gdGAW8BE4ABxpg1dp01wAB7\nfxCw0NW+FusLoNHed1hplyuK0glJtWr2iiuuYPz48bFVvsm8f3a46yY298xnj0fnUHvtNb7Xc3BG\n7KnMTV4mTJjAOeeck1DeHkf9EE755wP7AZcYY94QkQosE08MY4wRkayl33Lb1kpLSyktLc3WpRWl\n3RKU+am9mh1S4fdc0cgKTPEQunb1VC4upvsjD7Dy0Tm+cwN+10p3QVkuvefKykoqKytbdhFjTNIN\nKAIiruMfAv8CPgSK7LKBwFJ7fxIwyVV/PjDavs6HrvIzgHt87mcURUmkvLzcYDlcxG3l5eVtLVpW\nAEwkEjGRSCRuv7q62gwEs+74U02j5JvrT602TU0mrp77GqneUWNjY+za3ntFIhHT2Ni47R++hdh6\nM6U+d2+hcviKyH+AXxljPhaRqUAP+9Q6Y8zNIjIJKDTGOBO+DwMHYk/4ArsZY4yIvI4VU3WR3YHc\nYTwTvprDV1H8CUp2nksj0pbgl7e3O/DEwQdz+GsL6YWlF+ZzDAvLH+Kooz5MyJ8rInHRNh1CJD/P\n2nO0Ba2WwF1ERmG5enYDPsFy9ewCzAOGkujqOQXL1TMKlBljnrfLHVfP7bG8hy7zuZcqf0VJQXtT\nWGFMVl7F3W3xYna65BLyV66Mlf2Dk7mKW/iUuWCHYvMq/3TfS3t7l360mvLflqjyV5TUtDeFNXXq\nVNcq3Gbc/vMJz7RmDY0lw+n67Qaq2YfLmcn/tivlu+/OoKrKyqrrrAHIz8+noKCAUaNGJcTZTxWq\nob29Sz9U+StKJyGZwsrFieEwJiu/Z/p25j3MuEW4YdWvKBrUhX/+Ew44ILyOCxOqQZV/jqDKX1FS\nk0xhhRlltyVu2evr61m/Zg1569Yx9NBDfTuFTz+FSy6Bv/4VBg2KNw85PvfehVmOe6gq/yRtcu2h\nVfkrSmrCjvyzPTGcja+KmOzGMO+00zjg8cf5AsuN0CFZR+V+dmeS2Kvka2pqYs+uyt8fjeqpKB0M\nryLOZoTK+AiazaT9VfHOO3D55Yx75RXAWkn67K/vZ9iEI9h+e9MhvJdyHVX+itIOcKc8dGiLGPTu\nVa7er4qw3Ayw//5gj7ZNv348NvIGzrz3TE7bmM/f/gY+Xp9KllHlryjtAL+Uh20Rgz4bXxUrwFL8\n+flsOvcSfvbedcyv6gPAww/DkUdCivA7CaQbqqGjBWnLBFX+itKOCJrYzGW8Xy33AH8480xeP/B8\nzpp2CKtXN6uhcePgtNPSv0e678BrvurIoZuD0AlfRWkHBE1gpprYbM3JzJTX/vBD2H13alasCMjB\nex/WWlDIy4Pp0+Gqq1KbfLwTvm6/frcpCoLNYbnoDtsSdMJXUZS2Z+1auPZa+MtfrO1HPwISv1qW\nL/+CH//4bXbYYRTz5nXh6KODLxlkpoFE01MYU1R7VfLZRJW/onQg2nRieMsWuPNOuOEGWL/eKisv\nB9sk4x+S+XDuvfcJdt+9P/X1wQo5yEwD8c+nhEeVv6K0I1JNbLbZxPDnn8NRR8GyZc1lxxwDt9+O\nFX95+6CGnHbagUBye7vjZeQ3x+FvUlJSocpfUdoRYSc2t/nE8ODB0MMO9rvHHpbS/8lPiG4Vbimr\nB95i40Z/k7QjazIzjNdM0x4nvnMNVf6K0g5wUh46+E1seutne4Sf0qRUUUH+4sVw8cXQtStr1sAZ\nZ8DLLxcChVx22Vruu68mNqHrXCsTWVvj+TobqvwVpR2QKuVhNkjlAVNbW8vwkhL2BN6zzyWYlOys\ne6++arltrlrVfJ2XX36LXXY5Bfg2q3IrmZEygbuiKJ0Dd4Jy9+YkKN/+5Zd5D3hrhx34bPHihETq\nDkuXwhFHNCt+EcOECXUsXz6CSGQJEJ+o3fvVomwbdOSvKB2QdFe8QnDohr5ffAFjxzLg+ecZALB+\nPUOffNLy5PFhxAhrhe6990K/fvDQQ8Ixx/QB+sTquL9aMvE+SvZ8nXG1biao8leUDkgmk5++oRtm\nzYIrroCmJgDWA1unTKHv5MlJrzVrltXk2mth6NC0RUlJsufrjKt1M0GVv6J0ABx7fTQajTPDODHt\nHdwmlmQ2/hgHHWRp8bw8Gs44g+EPPcSi88+nb7duSeXp3t1a3+XFybzlnihOZx2C38R3dXU1DQ0N\nABQVFcXa66g/Oar8FaUDEBRqGYInhpOFZ44xejTceiuMHcvXPXuy9qGHYiaWdevymDx5R047rS60\nnE5b90RxOusQ/Ca+99lnn9D3V5oJm8C9BmgAtgKNxpgDRaQv8BiwM4kJ3CdjBe3YClxmjHnBLncS\nuHfHSuBe5nMvje2jKCnwxtUJSuBSUlISWK+hoSE2Yp4wZgz/eOMNmnbckcLCQvr06ZMQt8eJI2Rx\nOPAwMAhYBuxHJPJeSg8k5xotybyV7D10Vlozto8BSo0xX7vKJgELjDEzRGSifTxJREYCpwEjgcHA\niyIy3NbodwPnGWMWicizIjLWGDM/HYEVpbOSKgxxmFDL3tH+YGAa8CZYdpo//znw/sXFxSxfHuGP\nfyzkj3/cgaYmR9cM549//ITi4r6BsjodzYYNG2LXUj/9NsYYk3IDIkA/T9lSYIC9XwQstfcnAxNd\n9eYDBwEDgQ9d5acD9/jcyyhKR6eurs5EIpGEra6uLrBNeXm5wRqIxW3l5eVx9dx/Q96/J+e+NR98\nYKaC2dq9u51QEWO6dzfmyy9j7fzk+8lPtsSqgzH9+xvz3HPhZXW2SCQSVz8SifiWp0L1hYX9HkLp\nc2dLZ+T/oohsBf5sjPmLrfjX2OfXgOUFhvUduNDVthZrgNFo7zustMsVpdORSTpEtyumm3QmNgsL\nCyns2hX23JNygM2bAfgHcOBzzxHdtAk2bQL8Y+aceurjwM/so5dZvPgIBg2yjtyj/VNOOYWjjjoK\nsCadnUnZDRs2cOyxx4aWV2k9wir/Q40xq0SkP7BARJa6TxpjjIhkzfDm/s9fWlpKqb1qUFG8tNe4\n7JmkQ8zaM/XsCWPHWmaeffeF22/n1COOsFZmBeCOv9O/PwwcCOXlRzFo0NZYnWSTzs6krGOqymQd\ngtJMZWUllZWVLbtIup8KQDnwWyyzT5FdNpBms88kYJKr/nxgNJZpyG32OQM1+ygtJKwpJJfJ9v95\nXCYb936cSWn1anM2GLN1a6wNYKqqquLMPFVVVQnmmKYmf7ndpiyvDA7OcdCmZp/MIAOzT0pvHxHp\nAXQxxmwQkZ7AC8D1wFHAOmPMzSIyCSg0xjgTvg8DB2JP+AK7GWOMiLwOXAYsAv4F3GE8E77q7aOk\nQ5CXS66P/N1k22NF7MhpBcCxwKN2udek5M2IBZar5ZAhw7jlFusD4ac/Dc4Ulkxu77WdfW9wuLCZ\nt5I9q+qL1vP2GQD8n/2fIx94yBjzgoi8CcwTkfOwXT0BjDFLRGQesASIAhe5tPlFWK6e22O5eqqn\nj9IispFQvL3jVahdgKVXXMHO991H1/XrmfHEE2zdf/9QneEXX3ThnHPglVegWzfYbbeuSesHmd2C\nyEaAOk2+niXS/VRo7Q39jFMypL3+32mp3G5TytFg3nO744BZf9ZZCW0aGxsTTDPwc1NQEI3z5vnZ\nzxrizEHV1dWmqqoqZg4qKysLNOGEeb5Mnr0jmPqyDa1h9tnWqNlHyZT2agJoqdzOwqllV17Jbrfe\nGivfPHAgZ61axYxPP2WYx3MnfsEWWNN4zW3z8gyXXrqek056nyOOSB0nqKysjFmzZvkuLgtrHgpL\ne53kb00yMfuo8lc6DJ1d+X/2zjsM/dGPIBqFq6/ms5NOYtiIEVRXV1NQUEA0GmX16tUArFq1inHj\nxsW8eBaWrRDaAAAgAElEQVQtWstpp40AemMt6/l/wGuxezj1nNW5AAUFBRQUFADErQqORqN07do1\nriMIsuu3198s12jNFb6KomSJ1kqy3lRYCE88AXvtBQMGYOxrzp49m1mzZvm2cVbaFhcXM316PZWV\n3/L886OADTElP2bMGIqKimJtDjrooKTytTR+j7Jt0JG/0mFoL6PIRJNLPEmVozHw979bjvaHHhp3\nPW87p9wZ+Tvxcyorq+jShYRYOsbAZ58lly2ZfM77Tyd+T3v5zXIdHfkrSjsi7STkCxfCb38L//uf\ntTjrzTchLy+uvRvnuKCggGHDhvHddwLczA037Mdf//plwuXFozrSkc8dqjlVbt6GhgbfkM6d2Wbf\nFqjyV9otrWU+2VaEDm4WicCkSTBvXnPZ55/DsmWwxx6xomQdx9tvwxlnDAR+x0svwaOP9gKSr7RN\nJ/ian6knCK8ZSpOvtA25+5ehKCmora1NUDYdzra8dSv86EfgdGrdukFZGUyZAvYo2S/BiXO8ZQvc\nf/8Qpk2DaLQ5Acsrr2wPZJbxy8Hd+Tr/uucJghg/fjwTJkxIKNdR/7ZFlb/S7knbfNKe6NIFrrkG\nfvUrOP10mD4dPB1asoVT99wDN9zgPrOJK65YxYUXwpo1lqL2ev6EfX9+na/3d/DWh2YzlNK2qPJX\n2j3tNTZ8mOBm9fX11JeW0u3pp9nyve9ZhTU1oe3j551nhehfvBh+8IPNvPnmKG6//RNuvz2xrvc9\n1tbWJriHuuVMNdrvMB1wB0WVv9LpyJVFQm7leABwMYl/kJmEfnbTtSvMng2VlXDRRfl07/6Jr/99\nSUlJ3KjdK1+ycm87txnK+Ypwf5156yttgyp/pdPRUoXaUtzKMb+2lqoxYzjDPvfDyy8nGo3GJq7P\nPffcUKGfN2yw5n/32y/+XtFolD59ajn5ZLAH7gmyOJPiGzdujCWBB2tE7zUJOQSZhoJCM7fXr7MO\nTbrxIFp7o53GZ1G2PU5MmjBhiN0EhR5OlkUr69TVGXPVVcZ06xYLpLMZzJQk8WqC/jaef96YoUON\nGTDAmHXr4uuFDaHs3CtVPb/rxscHyk6oZiU9aMVMXoqSs6RrW86JSKAvvAC33BI73HLyyay9/HLO\nHzKEaZ4QCm73VTfr1sGVV8IDDzSXXXGF/+3CjNz9MoWFmfytra2lqKgoJrP3a8HPpKS0Par8lXZL\nMhdH53zO8vOfw+23Q5cuHPS//7HwH//ALW0qhfv00zB+PHz1VXNZv37w4x/DnDmJ9cOYXTKd8wiS\n1X3PXF5v0VnRX0Rpt2QjNvw2wZjE5bMi8K9/Qd++vO5apesQNCnr0K9fvOI/+eQtTJmymh13bAKa\nF7s1NDRk5RGCPJPcXxSpJpCV3EKVv6K0Fp99BldfDfvsA7/7XeL5fv0Cm7o7Mb8O7ZBD4Ne/huee\ngzvvhLffnsYBBzRPYjsdRVlZWVIR3Uo9WZiFMKN7r6w62s9tNLCb0mHIJEhYqwQWW78ebroJZs6E\n776DggJYvhz69w8lgzcF4tathry8xHoNDVZon169gt1XGxoaGDVqVFoL4crLy5kwYULM82f48OFx\nvvxhwjOLiO/Xi8bvaR00nr/SKuSKX3wqtqXy930nTU3s9Pjj9LjlFms21uHnP4eKChg0KJQMzvHa\ntbDTTnO5+OKzuPPOzGRNFUHU2ymA9bsGucMCoRKziNfMZaPxe1oHVf5KqzB16tQ29YsPy7ZU/kHv\nZPkuu7Drp59aBwcdBLfdZtloiO8wGhoaYvZ492jaSorSj3vu2cqUKfD119b0wGuvwUEHpS+rX8J0\n9yh+2bJl5OfnJ3TkQR2+O0OXl6CRv5tcGzB0FFpN+YtIF+BNoNYY81MR6Qs8BuyMnbzdGFNv150M\njAe2ApcZY16wy/fHSt7eHSt5u68xUpV/7uFWBLn8CR9Gkfspw+rq6pgiLioqipkxkj1f0Dvp9/nn\n9D7/fPjDH+DUU6lfvz5Wr6KiIjCpisP48X9h9ux9gf3jyi++uJ4//alPoMklLC0dkYdNyahx+rct\nran8r8D639jbGHOCiMwAvjLGzBCRiUAfY8wkERkJPIy1Wn0w8CIw3BhjRGQRcIkxZpGIPAvcYYyZ\n73MvVf45jPuPOtfMQWEUTioziJukCnHdutiEbcJ9m5picfaDvhCgOSZOjx49+OabbwCYNWtPnnjC\nPREcAS4CEv5UMopa2tIReSrlr3b+tiET5R9mxW0xlhI/AnjaLlsKDLD3i4Cl9v5kYKKr7XzgIGAg\n8KGr/HTgnoD7hV3UprQB7t8naEWoe1Vqa5Puat1Uq4JTXqO21pjx443p2dOYzz83xgSvvE0mn7tN\n/OrYvga+MvCNgesMdE97BXMyWvr35W3vfb62/v/QWaGVVvjOBK4CClxlA4wxa+z9NcAAe38QsNBV\nrxbrC6DR3ndYaZcr7Rj3itBkcWdaE+/EZNjEIMkWPfmWr18PM2ZYHjzffmuVXXcd3H9/UvnSWU1c\nVVVFUVER//znKoYM+Za1a3fikks2h86IFXTP1iRoYrisrCwWs19H/blJUuUvIscDXxpj3hGRUr86\nxhgjIlm107j/aEtLSykt9b210sbkQpgEv5AEkGWF89xzcNZZ8auqfvpTuOqqjC63YgXA31i40JoT\ndnC8bq68cu9Q1wlKzJ7tiXivec+9HiDZ+1el33pUVlZSWVnZsosk+ywApgErsAyPq4BNwINYZp8i\nu85Ams0+k4BJrvbzgdFYpiG32ecM1OzTLgn6fXLxd3ObJFKZTNxmC2/7Ff/+t2nKz7eCr40aZRqe\neSaubdhnh57mmmuM6d7diuW2//7GRKPN93abT6qqqlLKWl1dnXaAukx+p1ww7ynJIdtmH2PMFGAK\ngIgcDlxpjPmlPeF7NnCz/e+TdpOngIdF5HYss85wYJExxohIg4iMBhYBvwTuCN1DKUoGJPNVD5NI\nxd1+ErAceKK6mvI33mDqcceFlmPrVpg7F+Bjfv/75vK33oKXXoLhwxPbZJIRK+jLK9nIPczofJt8\nXSnbnrC9BHA48JS93xdrEvhj4AWg0FVvCtbfyVLgGFf5/sB79rk7ktynNTtIJU0aGxvjRs64RpiR\nSMQ0NjYaY3J/5O/I7Z7Y9W57ginxGfmnGlmnevb1643p398a7Tdvb5pXXrHOB4320wmTnEwGHbl3\nfMhg5K+LvJSkpHKNdNwNs7pSluzbjB35/Pz8P1+4kMKZM+n1+OP8X1MTp/g8R1j/9iD+/Ge48EJr\nke+0aXDOOXkY0xyEzfuOI5FILJ8uBGfE8gut4CXXXHKV7JOJq6dGXlJC0VpJ0rd1Vq24SKD19fwB\nGFxaSt7mzQCcANT+5z9Ehw4F/BdSeTsQaDalrF/fQM+ePcjPz4/Lf3vAAQXAX1mw4Lf06GENvqPR\nKPn5+b5pD7OZEUuVvOKHjvyVpDijUu+CIm95Nkb+rbkwKEG+LVtg113BpWSfwrLtf+hq5/d8/l9D\nPYAy4ETgUKwF7snxG8WH/dJK+XxKp0JH/so2w8nz6hcSGMKHHmgzd9Fu3eCXv4Tp03kN6PnHP7LP\n8cfzrH06zJdNVVUVAwYU8/jjvZg1q5Avv7Se96qrPuSii7rGlLg7lo4b97Gj1J3OwG+lrMbIV7KJ\nKn8lIxxzhluBuUesmYQecLNN7NS/+x0nTp/OU0Dk+OPTlveTT3bjvPOK+Pjj+PJXXy1mxoztY8du\nhZ3KfOZ0mKli5CczPTn31Hj6SjL0f4cSCq+iWbVqFdAO5gJeegmefNL/XGEhTwU0S/Zl45Rt2JAX\np/gHDIiyZs2FzJ17DTDM97qZ2Oz9qK2tTTAPZbPzVTo+qvyVUITN5JQtvKEj/BKaJ/sKiL7+Oo1X\nXsn2r74KwDGkNzJO9WUDcOyx3/DAA1aelsmT4YQTVjJy5H3k51+TxpO2jNbqfJWOjyp/JSlBSdJb\nW8l4FbvfvXy/ApYtg2uuIX/evNh/7k3AUJKPjJN92dTV7UJJSSO9epm45/7ii1puvrkr/fptpbCw\nKdBDJ6g8G7RW56t0fFT5K0nJlSTpQaGCE/jnP2HePMCKJrj2lFOITprElP79mUKwK6V/R7YPN9yw\nHwsW9OAPf4ApUwjRJhF3vTArixVlW6DKX2kRLVFm6Uxahg1jsP744xk8cyZf77UXBy1YwEu33eZb\n309xl5WVMWvWLJ5/vpbp07ejsnJHFiywzt16K1x8cXPdefPmMXDgQAoKChg1alRCkhUHr/9+Jl9L\n7nfi7DvJZxQlU1T5Ky2iJaafFk1afvMNbLcddOkSK3ImiXsBG7/4Ium9g3LXzpo1n7FjB+N1mT/y\nSNi0qfl43LhxgGV6AqtDcjozvw6wqKgo9kzOc7o7jF69evkqeaeud7+szEqEp18SSqboIq9tREdZ\nYp8sTIJ39JvK1dBZ0BQ0aem3wGo7Eb678074/e+t4fj/+3+xdt4FY5DYgQQtWnM/36GHfsN//2u5\nao4Zs45JkzYzcmRjnGyOzFbO3T4YY0JlCXOew7soK1nGr4hP5q2GhgZGjRoVeB/19ulc6CKvHKY1\nwxi0RceSzbmAUJOW0SjMnctHAJdcYpWVl8O4cdC1K+D/vMlGxo2NsaZx/Pe/PwKuAa6jquptbEej\n0DJn4oHj9m7yJnd3cD9fNBr1nYh3y6coSUk3Elxrb+RgdMhskG66wXRo7aiNYWTP5HdzXy9Z+c5g\nzO67x4fF3GUXY+bOtQLi++D3Ppq3I81BB31jLr7Yv10kEjHLli1LGmnTiWbqfvZUz+N+R8neVya/\nZ0f9u1HCQQZRPdtc2ScI1An+E2f7GVuzYzEmnDLKpvJ3K95IJGLywHy3667GgFkBJvqnPxmzZUvS\na7vfg3UfMdOnf2D23LMh1n9069Zkli6tT2iXzr77OFvK3y8JTarfszP83SjBZKL81ezTAWjt+Dhh\nk3lkmk/Wa5p55513gHiTxymffMJw4G5g6bHHMszPXuPBeQ+NjQBvM3nyyLjzW7Zs5eqrX+Cvf/2x\nb7KT1iDMtdvbPJDSPlHlr6QkrDLym+x0J/IOupaj5A8B+gP/tMv9bOcVGXgXWf1ENfB9u2QzZ565\nhQsuaOB73/txYBL4TEk2z9DSaytKtlBvnzagNcPvtlVoXxGJc2N0fOa9uCe4HY+hbu+9R+Htt9Oj\nspI1wC7AN6Sf0ASslIlduiS+B5F9KCh4l/PPh9tuG4gxq2LngibMS0pKfL1zRIS6urqEUNSpJnbd\ncruvnQ00pHPnJhNvH1X+GdBS75qOqvxTKUdvva4ff0zhzJn0nD8/dh3TrRvv33wz+1x+eVxnEkQk\nEmHnnYfx3/9GmTbtO5qa4J571iZ4v5SUlLBpk6FHj/DvyO3W2rVr11Cdm/tLJ5n7a7Z/J1X+nZtM\nlH+bT/B6N9rBxFVLvWta8xnb6v0RYkLUmPh394rLeyeal2fuAmNqaxMmrnFN/jqb5YHTzVx99Udm\nn302xzkCwe6+v0+QfMmeqflewd5DbtmSeUBle2K+tSf6lfYD2fb2AboDrwOLgSXAdLu8L7AA/wTu\nk4FlWAncj3aVOwnclwGzktxzG7yqltHSP7psPmPYBOutQdB7SKb83W1+BKYpL89sOPVUs37x4gSv\nGa+bpZtPP40Y+MCj9K3tuuu+8nQSJJUviGRK3n1tr2x+z21M9l1yNTG74pCJ8k9p9hGRHsaYb0Qk\nH3gVuBIr1elXxpgZIjIR6GOMmSQiI4GHgQOAwcCLwHBjjBGRRcAlxphFIvIscIcxZr7P/UwqmbLJ\ntjDhZGM1bBCZpv3LBslWpbrfSaEI9T7vSEQwkQjY8jnvMiiheeJK3ceAiQB069bEli0P8NxzxzB2\n7GBPvcT3E+b/mDOP4V1x7L120Dv2/t/I9mK8jrJqXGk5rbLC1xjzjb3bDegC1GEp/8Pt8jlAJVb6\n0xOBR4wxjUCNiCwHRovIZ0BvY8wiu81c4CQgQflva7ZFAvFUMWzcE4GZ/uG2RVz3IBfQ2LN99BHc\neCOfAaxdC/37J17ER2m6Uxl+/PEqjjnmvAAJ7gLGMX16CUcfXcv++5/HiBER35pOIDawvIv8Asj5\nBZrLJtlWyqrklZaQUvmLSB7wNrArcLcx5gMRGWCMWWNXWQMMsPcHAQtdzWuxvgAa7X2HlXZ5m+NN\nGhJJFi64hQQpaLeSzrTTaYu47kHKZzhY+XEffhiamtgB4Lbb4KabQl23S5d84AB+//thPPLIUOAR\nVqzwc5/8HNiVSZOaqKlpSnpNJxCbg18AuVSddBANDQ2tujZAUVqDMCP/JuD7IrID8LyIHOE5b0Qk\nq3Yat/IrLS2ltLQ0VLuwn8FB9aB1Y9UnU9DpdjrOM6QaqbrPt9ZI0T1i7j1nDh8C/O1vABgRHjWG\nM846y7dt/G+Rzw03rOORR3oDi7jvPoA8YBSHHXYg8IbPFeL/6yXzsQ/7deSt98YbbzBu3LjAa8+e\nPdvX8weaF3XpKF3JJpWVlVRWVrbsIulMEADXYtn8lwJFdtlAYKm9PwmY5Ko/HxgNFAEfusrPAO4J\nuEfGkx5hJ8CC6mVy7zBt0ln2HxbvMyS7dtB7yBbue+1jz7puBfMomJFJns+RKV7OpT6TuO+ZBx5I\n9AJyT976Pa93C3pHTnnQ75TK26e6ujpOnjD/BxUlm5DBhG/Skb+I7AhEjTH1IrI98GPgeuAp4Gzg\nZvtfJ0P2U8DDInI7lllnOLDIGGNEpEFERgOLgF8CdyS7dyaENeEE1WtPqy+dZwjKTOUd8UL4/LeZ\n4oyYv7r3Xj4fMYK87bdniW1uefXVV4lGo6xdu5b+/fvbph045ZRTOOqoo1i1apVtmvkLcCvWMq/H\n7OPXOPzw5klVv6+noHST0PL5DycWv/N8qSbsIz4hmHXUr+Qaqcw+A4E5tt0/D3jQGPNvEXkHmCci\n5wE1wDgAY8wSEZmH5RYaBS6yeyWAi4AHgO2BZ42Pp09LCRvjprVj4QQRJvFGWNNVmBy3yc65FyOF\nNYvF6r31FkybBjffDLvtFjvvmLVqfv1r9vd0pJYMPYFTsMYL/wFIiEn/1FM/Y/HidZx44iYKCo6g\ntnZ4KMXdmukmHcXuNtslu7bG0VfaBel+KrT2RpZ84MNex10vk3uHaRPGJOGQru92Y2NjnCnEve/n\nfx/m2kEy/PXcc40ZO7bZHnPeeXHP5zWfVFVVmeXLI+Zvf1tlTj55g9luu8ZY00GDGs0nn/j74vu9\nO6fc+76D3r+73C1PMj/9sPWS/ebZ+v+rKOlAts0+SiLp5J11lwWZJJxjh3S9j8LkuHUfu01cYc1i\ntS+8QN/rrmP7++9339iKmBb7sEukuLiYaHQYZ56ZeG716nwaGobx/e8nnkuF17Mm1ft3yCThuqJ0\nVFT5p0kmeWf9TBJewiQuzwZhEqF7zUCDR4ywTD1g5c391a/gqqtg553j2jnumO7Ocbfd4OCD4bXX\nnJIPmDhxIJde2pfBaTj7JouMmer9p+p83UHiwtRTlI5Ap1b+0WgUiFe8YUeRLV1UlXHi8rZgyBC4\n6CJrtP/b34K9WArgww/hz38uBJZw2GG/Av6X0Pzii+GAA2Cffd7hV7/ajwsvjDB4cF/fWwXNi3jf\nbTrvP+x8QGvOGyhKrtGplb+jWNyKN0gJe/3qs7Goqi1W5Qby7bcwdy7suy8ceGDi+YqK2O7SpZYb\n/5NPwgcfABQChZx99nNMnfp17B067+rQQ63tjTeWx5U7uI+Dnt/taQNts6jNj0zMgIqSC3S6/5Xu\nP1bn3zBKOCgMhJeg1Z5+rpU5ocDq6uCuu+COO+DLL+H44+Hpp5M2efll+MMfEsv/+98Cdt65IHYc\npMhTKXiIN7k0NDTQo0ePME8Tp4iz6dbq9YByrt3Q0JDgsZTTX3GK4pDuDHFrb2ToLRE2uqWf10sq\nDxNjmqNRBkVydNqUlZWl9KgJe1/vuwjK7ep+Vvd+XV1doDdTHzDm8suN6dmz2XsHjBk2zJiGBtPQ\nYMxbb/l7r6xc2Vy9e3djTj3VmHnzjNm0qfk+fvI5uXmXLVuW9DfyyurnfRT07lK9e7/3GkQqGdy/\nd7oRPxUlm9CZvX3SnYh1Fj6FNbN4R49Bpovx48f7JvMIM/J05iD8Rq9ghRG48cYbfdv6ma7KysoS\nrhGNRq0gZtYFYdMmAMyoUaz8xe94zIzj2ZPyqaqCggKwYvk1U19fz5Yt9ZSVFTJixBYOO+xbevQw\nFBYW0qNH8NoJ9/FurrUBYUbFbu+jVIvavF9xkJ0FVkFB7BoaGpg1a1ZufMUpShp0GOXvENaO3lLP\njaBOo6CgIJRHjSObG7/E5W6lXl1dzfjx40PPDTjxZtzXWLhwIcXFxWwAmi64gLw33mDrlRMZOeFo\nPp4YHxF23TqAQ+PK0o2CmmpS3emMkpEs76+X1lLCQaYjDeimtFfatfJ322FTTcR67fx+54KO3bjD\nDUN4d0C/VbNhbN+OPGPGjEnoWJJ1dFVVVciWLfz5yCNZg5VxB/IZM+ZIYAsA/zv2WIp/8xsAdh5m\n+HhZvPIfOfI7lizpGafgLrnkkrTWIaSaVK+trU1LWWfqjhlks9eAa0qnJV07UWtvpGHzb4kt2G2P\nDtqC7LVuGZPJSxo24zDP4F2JGlRvRzBfX365aezf3xgwkZ32MQeN/spAg4HTfeW47rqvTEGBMT/5\nyUYD5xoYmPKduJ8vKKtXdXV1RjbxsP8PwtbLJOtVmGs7z6s2f6UtobPZ/L0rUZPhdRWsra2lqKgo\nZvt3Aou5R9PZXNQTZDMuLCxk1qxZcRmjMiWvoYF7saLmdZ85M1be/cs11Hy5BejHccfN5M47p8fO\nNecU+IgzzxzEmjW1PPfc/Wm7oXrNQd55h7a2iSd7/24y/ULQVcFKe6NdK3+/P8ggE45X+eSKzdgh\nk3tu3Ch8+inssot13NSjBz/GSrwM8A7fZyaX8xinsYXtAFi7tpBhw7onXOvHP4637af7HlJNiAbh\n/D7Ou2kt00zYawR1YkFzGroqWGmvtGvl70eYEVjQH2tJSUlO/LH6d2BdWLKkKy++CC+91A94i1Gj\nhnLggd/yyCNWUrXa1av5O1Z+zR1v/BNjrv0NsJbjj49y8MEbuPrq0Tz66L+BYb73raqqcoVWTo9M\nJ0Sd36u8vBwgLcXbGoT9QnDQVcFKe6VDKf+wI7Bkf6x+niepVnEm81gJu+DLjX8HtgfHHTeY/XiL\ni7iLfnyfO5v2Y+FCQ0nJ7liZMi0qAK69GPgj8BETJ/7HPvMpkiTFc2uaLlK5ZjrvIx3F2xroBLDS\nWehQyj8dt8p0SLWGwPFY8esk/OYi3CNZy8rxQ/70J1i8eCjwEp98UkJentW2rKyMe2bNYhwfcTE/\nYDRWgLVPKeEuLmLX4U0sWzaQZcv+HUtCDlYn5pcjOJlnk3veIR0PKD+84TDCmtlU8SrKtkGsieLc\nQUSMV6YwCU5EBHc7STLE9dZzjr3XcKipqaGkpCRwEtRZQObUC8I9yt1hh0L22ssKjOZl+XLYdVdL\nnvolS+h16KF0qauLnY8i/JMdOZ/N1LEB8A8jEI1G6dq1a+jJZGNM6GcA/68r5wto6tSpgeEwgtpk\nQtBv1lbkmjxK58D+f5fkuz6RdjHyT3dhkYOfL35LUjWmWkOwatVGYG+mTXuUhoYiPv20K59+2pUr\nrnifX/zigLj29fX1bLfd9mBPxLp5911L+QPsMGKENaP71lusBracfTZywQXsP3Agb5M6miVYX0Sp\n1ie4A6YF1fO7V9Aq6iDbeUlJSYtt4uqzrygtp10of69LZ1VVFdFolG+//ZZXX301Vs9rX3dMMV4y\nW5UpfPVVHuvXw9Ch0KdP8xlHIQ4ZshR4lylT4ltu2DAk4WoVFRUsXjwIK6XhB/yAf7KGL1hBDUcc\nscB1W4Ebb+TLTz5h6KWX8vHUqRkpzzBJX8LWC+MG2pqKOF2PnNZGOyOlXZJqIQAwBHgZ+AB4H7jM\nLu8LLAA+xlpAWuhqMxlYBiwFjnaV7w+8Z5+bFXC/lIsZwm7GpEqh2N3AIFNd/ZlvULQZM4w54ABj\nBg5sNLDFOMHMHnnExF3bWeBz8snrYnXc28SJ6xIW/NTV1Zn3X37brJp8jdm4yy7GgFlx+ukJC9C8\nAelSLQALWPyRdD9svUxlCLpPpgQFuKurq2vxtTMhkwVkipJNaKVFXo3A5caYxSLSC3hLRBYA5wIL\njDEzRGQiMAmYJCIjgdOAkcBg4EURGW4LeDdwnjFmkYg8KyJjjU8i9yeegLy85m333WHEiObzzTbs\ni7j66ifo3bsv0ajQ2AhffdXAn/50BvB63DWrqqpYsGBPHnusN+vX51FXJ3z3nRW0bNSoS4E7gfhR\nZF3dVN54A7wfSO+8U8dBB61PWEOwxx5fAx/xwx8OZZ99tmePPSy5+/bdyM03uy6wZAmFv/sdhc8/\nD64vky6PPkoe/pPJSjO5NqJO1z1UUXKBlMrfGLMaWG3vbxSRD7GU+glYLuUAc4BKrA7gROARY0wj\nUCMiy4HRIvIZ0NsYs8huMxc4CUhQ/j//efzx1VfD73/ffNysDE/lD3/Y3dO6H3AUXuVfXFxMU1M/\nPvoo8RnPPfcqxo8/jYKCAgqsUJYUFhZyzz3emnXAF8yY8RdmzEhcuHTBBQ3cdNMIpk9PNIvE0bs3\nPPus9WEAfAJ0Of98uowfzydFRbE2frb8lnrhuEmWCzdMwDXFItc6I0UJQ1p/3SIyDNgXS7MOMMas\nsU+tAQbY+4OAha5mtVidRaO977DSLk9JXl7QmcaA8l6ApczcyrFXL2+9LcA67r//z9x//7QEm/EZ\nZ8CYMdC/f5SmppV0726AnpSUzIpLNBKNRuPuNWbMGAQ4COtFNXlvO2QInHwy7Lgjq448kuGnncan\nU6YwJMlkckNDQ+za2SKM+6of2eyAFEVpG0Irf9vk83egzBizwe1KaYwxIpI1/7ZTToGmpuZtjz2C\napivdJUAAAxlSURBVL7HSSdtpE+fXuTnQ9eusHnzembPtiaBvcrtrLPgmGOsydqNGz/je98bZk9e\nng+cnzB622GHeoxJdDGF5q8Pt0tjHnAY8DPgrJ492WHTJt674w7W77svY8aMiXUSAMWPPUZ+fj7f\n1dSQ7MU5yv7aa68F4r1w3DkJnGsnC5PgR1VVFUVFRaxevRpojnHkVujerwCNY6Mo7Z9Qyl9EumIp\n/geNMU/axWtEpMgYs1pEBgJf2uUrsSaJHYqxRvwr7X13+Uq/++2999TYfmlpKaWlpQGSPcVvf7s4\nwcQye/a/Yh4p7kVLxcVQWGhZW+rqUufiDZO60bH39r7/fqI33BD7/HGSpLxy2WVcahcNHz481s7P\nL9+PefPmMXDgQAoKChISubiVsHPtZGESILEjcN6dV6F7F4Y57qIaTllR2p7KykoqKytbdpFUM8KA\nYNnnZ3rKZwAT7f1JwE32/khgMdANKMEyaTuLyV4HRtvXfBYY63O/uFnsoPSMTsjcoM0b/jhVPWMS\nvUiqq6tNVVVVLCyxW4YE7rqr2b2nSxfzzaGHmgvALHzyyaShflOFBE62eds4Hi/u5whKK+l+/lQy\nhPEkSoZ6wyhK60IG3j5hlP8PsczWi4F37G0slqvni/i7ek4BlmO5eh7jKndcPZcDdwTcL+6hUilv\nR2G5FdmyZctieWKdvLFVVVUxZeaci3jyxgYpqeuvucYcAsZce60xU6f6K71Vq8y/wJj77jNm7drQ\nbpGpni8TpewmVUfg7VT9ZHXu73WlDKv8c801U1E6Gq2i/Lf1FqT8q6qqzLJly2JKfN68eQnKP1mS\n9qDRvhu3kuoNZt3kyWZTaalpcic579fP5AUoPbfsYZV/ssTzqZLFh1H+Qc/nXNvdKQbdJ2ikHlb5\nK4rSumSi/NuNL186tmmHdBOSuG3QUaDvrbdCo8ejaOhQiqzktllhW4YETicXrhtv5E1FUdo/7Ub5\nu0knSbtXkXYDui1ZAq+8Am+8Aa+9Bs89BzvtFFfvW4CDD4bPPoMjj+QXs2fz8BdfwMCBfJEsLrKH\nXHWLdE/epork2dZZuBRFyT7tUvlnqoz6X3opG4Guxx0Xf2LhQjjhhMQGzzxjLcgCHpk9m2nffQe2\np0pNTQ3RaJRvvvkmtjAM4N1336WhoYFVq1YB2XGLbI0OxO+LQ104FaXz0G6Uv6/CW7ECPv6YnosX\nMwnoe911sHYt3X/xC99rNHXvTld3Qffu8IMfQNBKVlvxO7hdJoOig44aNSrwGSIZpvcLUsruFIgt\nMcmkivjpllXdNhWlY5Cbyv/ccy1be2Mj/evrOZkABThjBtx5J/2B6QAPPgiAsVeFeTuMb/bdl2ee\neIJzZ86k/1FHWavHunYlCK9PvLOoCogtjBozZkygGcob+z7drxU/v3o37hSILYlmGTbiJ+ReRE1F\nUTIjN5O5eMrqL72U+iuuSFCqO9x9N31mzGiuB9QAy444gnEvvxx4D/cEZtBoNVkyGEgcJbuVpZMQ\nxV0eNslHsnoiEtcZOGRj1B0mqQ2ES6yjKMq2peMkcxkyxBqRd+sGXbtSuPPOFLqUqzPiHQbsDHyB\ntVT4eVupH9GrF5GNG2P1w4yY/VIwQvoeQ61NW0+8qpJXlI5Bbir/zz/3LfazTb8ckBJwxx13jGsb\nNGJ28MvT61yztRWu2tEVRdnW5KbyDyAd27SXsPW8MYHC0FJvnLayowd1OoqidHzalfLfFmQy0g/q\nJBoaGuIUatCIvq2SgQR1OoqidHxU+acg2ag+VZTL2bNnx0XiDBrRt5V5x93pNDQ0xOUMUNOTonRs\nctPbJ4lMYb1SgtoE4fXQcY6D8Hr4+N0n254xYZ83E9x5CdyoC6ei5D4dx9vHw7a0TTsj+2g0SlVV\nVSy5CaS/SCsbo+ZtNRmseWgVpXPRLkb+QaNSIOsj/2SkukZrjMx1RK4oSioyGfm3C+XvHv16bdPO\naNw7Ena3cdvig0bMXj9/r/2+pKSkTZS/LqpSFCUVmSj/No/f791IESM+bFaolmaP8sqRSq6wdRRF\nUbINGcTzbxcjfzdhR8ItHTF7R/FhRvWtOSGrKIoSRIc1+7SRHKGUf7rmJUVRlGyTifLPay1hOgsV\nFRWUlJTEJoud/YqKijaWTFEUJZiUI38RmQ0cB3xpjNnbLusLPIYVV60GGGeMqbfPTQbGA1uBy4wx\nL9jl+wMPAN2BZ40xZQH3a7cjfzc68lcUZVvRWiP/+4GxnrJJwAJjzO7Av+1jRGQkcBow0m5zlzTH\nRr4bOM8YMxwYLiLea7Y59fX11NTUxHzp3333XV599VXeffddgNi5Z555JtamsLCQYcOGJWy5oPgr\nKyvbWoRQqJzZReXMHu1BxkxJqfyNMVVAnaf4BGCOvT8HOMnePxF4xBjTaIypAZYDo0VkINDbGLPI\nrjfX1SZn8JpwRo0axZgxY2LZuZxzt956a1uKGZr28h9X5cwuKmf2aA8yZkqmK3wHGGPW2PtrgAH2\n/iBgoateLTAYaLT3HVba5TmFd5Wrs6agoKAgLk/vPffc0wbSKYqiZI8Wh3cwxhgRaXsjfRYIa6fv\n3r37NpBGURSlFQmzGAAradZ7ruOlQJG9PxBYau9PAia56s0HRgNFwIeu8jOAewLu5bs4SzfddNNN\nt+At3UVemY78nwLOBm62/33SVf6wiNyOZdYZDiyyvw4aRGQ0sAj4JXCH34XTnbFWFEVR0iel8heR\nR4DDgR1FZAVwHXATME9EzsN29QQwxiwRkXnAEiAKXOTy27wIy9VzeyxXz/nZfRRFURQlLDm3wldR\nFEVpfXJmha+IjBWRpSKyTEQmtrU8DiIyW0TWiMh7rrK+IrJARD4WkRdEpM2d+kVkiIi8LCIfiMj7\nInJZrskqIt1F5HURWSwiS0Rkeq7J6EZEuojIOyLytH2cc3KKSI2IvGvLuSiH5SwUkSdE5EP7tx+d\na3KKyB72e3S29SJyWa7Jacs62f5bf09EHhaR7dKVMyeUv4h0Ae7EWhg2EjhDRPZsW6li3E/IRW5t\nTCNwuTFmL+Ag4GL7HeaMrMaYzcARxpjvA/sAR4jID3NJRg9lWCZM5/M4F+U0QKkxZl9jzIF2WS7K\nOQvL3Lsn1m+/lByT0xjzkf0e9wX2B74B/o8ck1NEhgHnA/vZURe6AKeTrpzpzhC3xgYcDMx3Hcd5\nDbX1hr+30wB7vwjb2ymXNqxJ+KNyVVagB/AGsFcuyggUAy8CRwBP5+rvDkSAfp6ynJIT2AH41Kc8\np+T0yHY0UJWLcgJ9gY+APljztk8DP05XzpwY+WN5Bq1wHTuLw3KVoEVuOYE9MtgXeJ0ck1VE8kRk\nsS3Ly8aYD8gxGW1mAlcBTa6yXJTTAC+KyJsicr5dlmtylgBrReR+EXlbRP4iIj3JPTndnA48Yu/n\nlJzGmK+B24DPgS+AemPMAtKUM1eUf7uddTZWN5sz8otIL+DvQJkxZoP7XC7IaoxpMpbZpxg4TESO\n8JxvcxlF5HisQIbvAL6ux7kgp82hxjJT/ATL1DfGfTJH5MwH9gPuMsbsB2zCY5LIETkBEJFuwE+B\nx73nckFOEdkVmIBlkRgE9BKRM911wsiZK8p/JTDEdTyE+HAQucYaESkCsOMWfdnG8gAgIl2xFP+D\nxhhn7UVOymqMWQ/8C8u2mmsyHgKcICIRrNHfj0TkQXJPTowxq+x/12LZpw8k9+SsBWr/f3v3y9Jg\nFMVx/HuCYQ7RYBW02MRiMAiWlb0CMRh8EQr6XswWMRhFjBb/bDo1iM2yoG9AhGM498EJlhV34P4+\nMNhu+nHhOc927r3P3P26fD4hbgbDZDkbXeC2zCnkm8814MrdP9z9CzglWudjzWeW4n9DPOlzsdx1\nt4gDY1k1h9zg9yG3iTEzA46AZ3cf/TOBNFnNbL7ZgWBmLaJP2SNRRgB3P3T3BXdfIn7+X7r7Dsly\nmtm0mc2U922iTz0gWU53HwJvZrZchjrAE9GrTpNzxDY/LR9INp9Eb3/dzFrluu8QGxPGm89JL6yM\nLGJ0iUWMV+Bg0nlGch0TfbVPYl1il1hwuQBegHNgLkHODaI/3ScKao/YpZQmK7AC3JWMD8BeGU+T\n8Y/Mm8BZxpxEL71fXo/NdZMtZ8m0Sizw3xPfVGeT5mwD78RTiJuxjDn3iRvogHiy8tS4OXXIS0Sk\nQlnaPiIi8o9U/EVEKqTiLyJSIRV/EZEKqfiLiFRIxV9EpEIq/iIiFVLxFxGp0DdMVntT6ZXrsgAA\nAABJRU5ErkJggg==\n", "text": [ "" ] }, { "output_type": "stream", "stream": "stdout", "text": [ "Reduced chi-squared (quadratic fit): 1.25260026575\n", "Reduced chi-squared (linear fit): 1.85413915105\n" ] } ], "prompt_number": 393 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both reduced chi-squared values are reasonable in this case (the quadratic is slightly better), so we turn to the (normalised) residuals to try and identify which function is more representative of the data:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "fig = plt.figure()\n", "axQ = fig.add_subplot(211)\n", "axL = fig.add_subplot(212)\n", "\n", "# horizontal lines marking 0, +/- 2\n", "for pos in [-2,0,2]:\n", " axQ.axhline(pos,linestyle='dashed',color='k',lw=1.25)\n", " axL.axhline(pos,linestyle='dashed',color='k',lw=1.25)\n", "\n", "# normalised residuals\n", "axQ.plot(x,(y-quaddy(x,*p_opt))/y_er,'ro')\n", "axL.plot(x,(y-lin(x,*p_opt_lin))/y_er,'bo')\n", "\n", "axL.set_ylabel('Norm. Res. (linear)')\n", "axQ.set_ylabel('Norm. Res. (quadratic)')\n", "\n", "plt.draw()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEBCAYAAAB13qL/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXuYFdWV6H9LECECdkedaFTQ9M3kJqiJioKjQ3dI6NZo\nnJvkRmSMouMd8zDgJPGBNKTb1x0fkQlyJ+ZtTO5EzcPJRMg1EKFBTURiUMBXkICKiZgYWoPGB7Lu\nH1WnOY8659Tz1D7d6/d99XWfql27Vu2q2muvtdfeW1QVwzAMY+iyR94CGIZhGPliisAwDGOIY4rA\nMAxjiGOKwDAMY4hjisAwDGOIY4rAMAxjiGOKwDAMY4gzPEwiEXk3cCiwC3hKVR9PemERGQmsBPYC\nRgD/paqXJc3XMAzDiIZUG1AmIocBnwM+BDwL/B4Q4EDgYGAx8G+quiX2xUXeoqqviMhw4F7gIlW9\nN25+hmEYRnRqWQTXAt8AvqCqbxQfEJE9gfcD1wGnx724qr7i/zsCGAb8OW5ehmEYRjyqWgQNubjI\nHsBvgDbgJlW9JDdhDMMwhih1O4tF5LMi0lr0u1VEPpPGxVV1l6q+D8/VNEVEOtLI1zAMwwhPXYtA\nRB5W1feW7XvIr8DTE0RkPvBXVf1S0T6bEc8wDCMGqiph04YJH93Dd+EAICLDgD3jCFaMiOwnIi3+\n/6OAacDa8nSq6vzW09OTuwyNlHPl4sV0d3bS095Od2cnKxcvdkbOnvZ2FCq2nvZ2Z8szyXOY29ZW\ncp9z29pKnkfeMjZLWQ42OaMSJnz058BtIvI1vKihTwJ3Rb5SJQcCt/hKZg/ge6p6dwr5GhmyaskS\nfn7hhVy9adPAvm7//ymnnJKXWAPs3GuvwP1vjhzZYEmyZ+mNN5Y8B4CrN21i/qJFTjwLo3kIYxFc\nCqwAPg18CvgFkLhTV1XXq+rRqvo+VT1SVa9PmqeRPdUqn2WLFuUkUSmds2fT3dZWsm9uWxvTZs3K\nSaLsGP7aa4H7h736aoMlMZqduhaBqr4J3ORvRgAdHR15ixCKNORsROWTRM5CS3j+okUMe/VV3hw5\nkpNmzcqkhZz3cw9j/eQtY1hMznypNaDsh6r6cRHZgOd+LEZV9cjMhRPROP4uIzvmdXVx1dKlFfvn\nd3Vx5V1peAyNsAS56ea2tXHSwoXmGhriiAgaobO4lkVwof/3FLy+gWKsdh6idM6eTfemTZWVzyB0\nvbhOI60fY3ATJnz0WlW9tN6+LDCLwE1WLVnCsqLKZ5pVPobhFFEtgjCKYK2qHlW2b72qHhFTxtCY\nIjAMw4hOaq4hEfk08BmgTUTWFx0aA9wXX0TDMAzDJWp1Fu8DtALX4IWQFrTLX1T1hYYIZxaBYRhG\nZFJ3DRVl/DfAQFyaqj4dXbxomCIwDMOITlRFEGbSudNEZCOwGW8hmS3A/4stYWneh4jIChF5REQ2\niMjsNPI1DMMwwhNmZPFVwPHAb1X1MOADwOqUrv8G8DlVnQBMBi7wV0MzDMMwGkSYuYbeUNU/icge\nIjJMVVeIyMI0Lq6qzwHP+f/vEJHHgLcDj6WRv2EYHquWLGHpjTcy/LXX2LnXXnTOnm0hv8YAYRTB\ndhEZA9wD/IeIPA/sSFsQETkUOIr0rA0jB6zCcQ/XJwo08ieMIvgH4FW89YvPBMYCl6cphIiMBn4E\nXKiqJUpmy5YtFelbWlpoaWmp2N/f309/f7+lzyn9z26/nb5LL+W6p54a2HfJE0+wY8cOPjR9uvPy\nD9b0SxYs4NqAiQLnLFjAuAkTnJc/q/Rrli/n/u98h1G7dsHee1c0WlyXv176SNSZ03o4sCLjebP3\nxJvq+l8CjgVNLa89PT0aRE9Pj6XPMf30d7xDFSq2M9ramkL+wZp+5vjxgc9l5vjxTSF/FulHg55e\nVh5z29p05eLFTSF/efoVK1ZoT0/PwAaoRqiHw4wsvhv4mKpWqqKEiIgAtwAvqOrnAo7r5s2bK85r\nFo081NLPO/FErrqvcqzhvBNO4Kp773Ve/sGa/qbzzuPa5csr9s+ZOpVPfetbzsufRfpFZ5/NDffc\nU5GmePJEl+Wvlz6LKSZ+iue7Xwq84u9WVU0c6ikiJwKrgHUwMJHdZap6l39c68nXzAw2f7rNTOom\nNktpJb0dHfSuXFm5v72d3r6+xguUMmnOPlrgDn8rJpXaWVXvJVwI66BjMHbg2cykbmKzlFbSDCvZ\nNbKhGHpkcR4MZotgsLaebWZSIyvSrBhdt5ICG4ptbXSFlC/NSefWVzsGjVmYZjCT9Upfebmdppxy\nihMfkjE4KLzHf3z2WeR3v+Orf/3rwLEkFrTrVlKj16Ou5Rr6sP/3M/7f7+FNPHdm6lIMQbI0TQej\n28kYehS/x/PwpjgoJmnF6HKjpdHrUVdVBKq6BUBEOlX1fUWH1onIWrwZSY2YZOlPb3RrYjDgQse9\nCzK4RPF7XK2iyqpizJtG92GE6SwWETnR79hFRE6gculKIyJZmqaNbk00Oy5YUC7I4BrF7/HOKmka\n3bnbKGXd8MCLegMNgGPwwjuf8reHgaOjDFaIu3niGVHp7uwMHEA0r6srb9GcxIXyckEGV1i5eLF2\nd3bq9NbWgXJYCTq3rGwuKxsA1gi55ra11RyElvb15nV1aU97u87r6op0HSIOKKtrEajqg8CRItLi\nZ/5iNirJSAsL44yGCxaUCzK4QLFltAroBq4GpvjHp48axYFtbYw56KCGd+422uXayD6MMK4hRORU\n4D3ASG8wMKjqFUkvLiLfBk4BntcGrIGcJS75d12PiKhFHuXoQky5CzK4QHFlW6j85wNPt7Yy7rjj\nuCDH93gwK+u6ikBEvgaMAqYC3wA+TnozhN4MLAK+m1J+qRC1MnLRv+tyREQ18ipHFywoF2RwgfLK\ndoq/9R55JL05j68Z1Mq6nu8IWO//Xef/HQ3cG8X/VCf/QwvXCDgW2ieWFnH8gObfTYc8yzGJP3Yw\nyZA3Ln9LQXVDo/spCnJ0d3ZqT3u7dnd2Bl6ftPsIgMIIjldE5CDgBeCANJVR3hRbAI9t2MDtL7xQ\ncryeHzANk9El11Je5Gl6u2BBuSBD3rhsGSVxuab1fWdlNYdRBItFpBW4HnjQ3/eN2FdsIGEKv7xg\ne6vkVasySmoyuuhayoNBbXoboXC9fyuOsk7z+86qwzpM1FChU/jHIrIEGKkZTEldjd7e3oH/Ozo6\n6OjoCHVe2MIvL9g48cpJWzE2AMzD5dag0Tia1TKq1vBM8/uuZjU/89xzJXVlVMJ0Fs8s+qn+PlS1\nIR28cW8ubOGXF2wnu0PWCtSrjJK2YgZzNEIUXG8NGkY1ajU80/y+q1nNhxxwQEldefnl0RaRDOMa\nOpbd004Xood+QwqRPiJyK9AO7CsizwBfVNWbk+YL4Qu/vGALIWtn7Lsv//3ww0NXRklaMeYS2U2z\ntgZdxPqdGkethqfW+L6jPqOsrOYwrqHPFv/2B5bdnuiqu/OekUY+QYStXIMK9q62Nj7TwOlozSXS\nGMJ+dIOhArV+p8ZSq+E59eKLA7/vgydPjvyMMrOao4QYeVFJjAB+G/W8OBsJwkejhHq5ELbnggyD\nmbBhwY2eRiArXA7DHIzUK++g7zvLZ0TE8NEwlfGdRdsSYDNwbZSLxN2SKAJVq1yN3YT96AZLBdrT\n3h54Hz3t7XmLlithYvDj5ht1jEGWzyiqIgjTR3BD0f87gadU9ZlkdkhjMH/z4CKJyyZsn9Fg6bi3\nfqfdxFncJuq7Fsdl49QziqI1Gr1hs4+qanatmGYiqctmqFkEjRwF6/L7WVwO3QHPNejZNso9mOUz\nIgPX0F9qbC9FuVjUzRTB4PFZJyVpBR32o3NlGoE0aIRr1PX3s/i96amiCMpdMY1sDGT1jKIqgjCu\noYXA74H/6/8+E3i7qs5Pzy4xqtHMg83SjL5J6rIJa7oPprEMjXCNuv5+xlncZrCuJ16LMIrgNC1d\nqP4mEVmHNzuskTFRXkqXXrC0wxfT8KeGrRitbyk8rvepFL83YQeLDsn1xOuZDMCvgE8Aw/ztTOCX\nUcyOuBvmGgptprpmoqdtXpff30rQ00eN0gsPP9w5v/RQIu8ZY+v1TdR6b6q5YrJ0DzaqvMjANfSP\neO6hL/u/7/P3JUZETvLzHQZ8U1WvTSPfRtGIFnjYwWa1TPTC8UZaCmm3FItdNs9v3Yr87nfc/te/\nwoYNsGGDG62qQUaY9zuvwZBhW9ZBrr56i9sMyfXEo2iNNDe8yv9JvPUI9gQeAt5dliZVLZkmjWyB\nh+lQqhaTfP6ECblYClm2fFyM7HE5ciYOUd7vOB2eScvLxXegmGr356pFUKui7gXeVuP4gcDlUS5W\ndv7xwF1Fv+cAc8rSRC6ARn2QUVw2ecpz+r775vLBZOnKyXOwVNDzdM0tlwZZVlhplJfLA+Zq3V+j\notKiKoJarqFfA7eJyAi8Seb+AAjeojRHA68BX0pgjBwEFA9M2wpMSpBfpI6YpG6dMCZeIzuGqpno\nB44aBWUL7awCNj7wAL0dHZm5irJ05eQ1EKfa89w+dixfcThyph7F38LWl15iBPDK734XmDYNF0Ya\nkUZODcYqo9b9Xekvt+laVFpVRaCqi/EWpTkEOAEY5x+6F2+Kia0Jr61hEm3ZsqViX0tLCy0tLRX7\nlyxYwLUBD2DOggUlBR30QV/yxBNs27aNY6dOrZp/f38//f3eUgx/2bUrUN7iFzGsPEH5h7nf4vTj\nJkzgfXPnctEttzDyzTeR0aM5adYslt54o1f5Fu4d+Dlw2/btsHIlAHM2bgQqK+Uk8hRkOu+rX+Wm\n887j2kceCSyHcRMmRM7/vdOnc8kTT3DdU08NHAvySyeVvzx9tQ/8rLFjK86BykozbXnSSL/sjjtY\nc8UVXPfUUwPvxtXAvMA7gh2qFd9kVHl4+eXAvHXHjtDfe+fs2czZuJFrNm8e2HfxuHEcd/rp9Pf3\n51ae/f397HzxxcD7K7wPhai04vyL7zsteSIRxXxIcwMmU+oaugy4tCyNBm09PT2B5tDM8eMDzcWZ\n48eXpKtm9k6sk39PT8+ADKNBTy87v9zECytPUP5h7jdM+nJTtHx05Up/39mtrRUum7TkqVUOhbKc\nCNru/z1/xoxQ5V8454wqpnXa5VnNHTEtYF+QGyWL55s0/cQieYvfjZWgc8vu55OtrTo6BXmmv+Md\ngeV1RltbJPnPnzGj5L0ZHVOeeumjvp8TA+4t6/dhxYoV2tPTM7B5VXv4+ljUq3AbjogMB54APoA3\nYO0BYIaqPlaURjcXafwC1TTgpR/4ANcuX16xf87UqVxz990Dv3s7Ouj1W8Ml6SZN4lO33RZaI69Z\nvpzVRS3waWUmXlh5quVf737Dpl+1ZAnLFi1Cd+zgqfXr+d5LL3n72d0CLNDd1kaXPwV3WvLcdN55\nVcvhmDPPHGiRDuw/7DA+5LsJit0Wfx02jIlnnsmxU6cmkidu+nldXVy1dGnF8bMnTOCAHTsqLZSy\nqcyjyLNqyRKWLFiAvPwyb4wYweRzzhm47zTv95qTTuKa1asBr1Owt1gGYBnwzD77cMjkyRx/7rm8\nZ1Kl9zaqPI+uXs193d0VbswTr746lfyzspgKFL+fQfmvWb684pyk70PU9P7iYVKRuBpRtEbaG3Ay\nnjJ4Ergs4HigZqxG2I6YRvXcuzhdQfG9h517JSm1yqHWs0ijUzHNzvpa95HmVAGN7HzO431QbZ6Z\ngePWFXnfHxEtglwVQV3hIioC1XAPoNETcrn0whffe0+VDz+LyItq5VAr+iOL+YWSVqiNeJ6Nnuum\nUEZB7qC8Gy5543J0Ui2iKoIwA8oqEJEPq+qdcc7NmjDTAzRyPpmspytIMl3uxgcegO3bK9IERV4k\njbKqVg61oj+GV4lQCRu5ksU8OEmfZ5hyrBeRluZAxvJvYdtLL3GBCPuPGeNMREueuBydlCpRtEZh\nI8H4gYjXSVlPDi6StniTzMiZlqsirtsoDK615sKWY9buMiM8Lrp3w8BQdw0NFsL4ttNwIYRxdWTt\nqqgmQ9KP0LXRp0nmjTr3gAP000cdpdNbW526p6GAa+7dWhTqjaiKoK5rSEROxwvzfElE5uMNJrtS\nVX+TmZkSgjjmsUuzc9Yi7EC0NOYtCePqyHp+lGoyJHXh5TUPTjWCyrHW4L7CfW996SVa/vAHFqxd\nWxLVU0yUZ9Es34ErNMtstMX1xtX1k5dST1MA6/2/JwJ9wKnAA1G0TdyNKhZBHPO4mUxq11bTcq1l\nHQWXWnPl5RjUOVvPVZQ0sqeZvgMjGsXvCWm7hoCH/L/XAGf6/6+NcpG4WzVFEKdicqUyC+PyCevb\nbpT/sln9pK5Rb3BftXey+H1IGtmTJBxyME2qNxgpfk+iKoIwUUPPisjXgWnANSIyEtgjquWRJuUm\n9ipgKbD1/vuZ19UVKxKjEYR1+YSNVGhU9NNgWrUrT8rL8Zl16wKjtsrfyeL3YYr/dz7wdGsr4447\nLtKziPMdOLuYShUa6fpyyc1Wrd4IRT1NAewNfBR4p//7QKAziraJuxHCIohjXudlESTpLBxsLXBr\nYebzPjSzNR2GJK6vqO+ka262YnlI2zWkXoX898C5/v/7A++IcpGA/D4OPAK8CRxdI13dGw5rXrtQ\nuUYJZ3TJt502rn1AeRHlnUzrfYjzHbgWhluLJK6vqO+kiwqy8J5EVQRhooZ6gWOAdwE3AyOA7+HN\nSBqX9cBHgK/FObnYxN56//0QMNtfuanrgnsjyuCUapEKLpmicUl7oFfcMsm7LKO8k2lFrsT5Dppp\nUFVcF3Ccd9IFd3M5hffkKgk/zRCEW7z+I8BRwIMAqvqsiIyJLuJuVPVx8CZGikvhhud1dUHARGBR\nKtdGkTScsdl8tdVI8wOKWyaNLMtaCiePdzLqNV0Lw61FXKUV551sJgVZl3omA36oKH6kEF6fwboo\nZkeNvFcQwzVUbgrl7fKJQhIT30VTNA5p3kfcvPKceLAZ3WDN4qqMWx/EeR9crntI2zUE/FBEvga0\niMj5wD8B36x3kogsw1vNrJy5muI8RS64fKKQpAXooikahzRbmHHLJMp5SVxIWcx3lAeuzZlVjbj1\nQZx3Mou6Jy93ZV1FoKrXi0gn8Bfgb4H5qrosxHnTUpCP3t7egf87Ojro6OioSJO3y6dRuGiKxnlx\n0/yA4pZJ2POSupAGi/LOkrTddHHqg7jvZJp1T5Jy6Ovro6+vL/7F65kMeFFChQVsBPhHYEMUs6NG\n3iuAY2ocT99mSkDeIY9Zm6LNGD4Xt0zCnldvArhGzAc12LEy8kizHEjLNSQiHwW+DrwBvCkin8Fb\nwOhp4Oz4qgdE5CPAjcB+wBIRWauqJyfJM2uiaOuszLss3WBxWiMuuD2StOTCnFetRf/81q2hyquZ\nOlrzwqwmj1zLoZqGwAvx/G/+/8cAO4EPR9EySTccsgiSDP5phs7BOK2RZoovj0u1cjl9331Dl1ez\ndLRGIU3rOK1ZdJt9gGKeFkGtSnht2e9U3EGRhHNIEYSt9JrVzA26v5Wg01tbq35czXqvUajmQrrw\n8MMHvRKsRtqNnaQuz2ZtfJWTpus3qiKo1Vm8v4h8Hq9fALyoocJvVdUFaVsnLlJw8zy5bl3g8fLO\nxWY1c8s7TwsL29+2fTusXAlUuj6Ggtujmgtp6Y03woYNFembMoY8Imm7BJO6PF1wUaZBnhGQtRTB\nN4ExNX4Peor95quAbiiZ5zuo0nMxsicM5ZX6UqiY07z842q20N24VIsMGexKsBpZhN5aWLVHXhGQ\nVRWBqvY2UA4nKW5phJ31sVlbyXFnxhwqobvlDBUlGESjQm/TlqeZyXx8QRQ/UqM3cu4jiNsZOhg6\nB4eC/9+IR9zQ25V4k0Se3dqaaoeuyyN80yBOHwgp9hEMeeK2NAZDK7lZLRsje2pZQ8Ut12cefnjg\nnEKf09XgWZpLl6ZmHQx266whfSBRtEajN3K2CAZ7S6Meg8GyMRpHrRXYki6xOZSJ45mgERaBiByj\nqg+mo4rcZbC3NOoxGCyboUTe02qXt1w72R1gUa2iacYO3UbTiD6QuK6hTwH/nOTCInI9cCrwOrAJ\nb+GbyoUFcsYqQ6MZcGGK8vLonUKAxYzWVhQCgw8GU4duVjTCTRtLEahqIiXgsxS4VFV3icg1wGXA\nnBTyNYwhhwux9EEt1ynAsuOOY9qsWXSXKSrrcwpHPc9EkCUYlVCKQETeCxxalF5V9Y7IVytCS2cw\nXQ18LEl+hjGUcSGWvlbLdai7WZNSa8XCapZgFMIsVXkzcATeGsO7ig4lUgRl/BNwa4r5GcaQwoVY\n+nqVvblZ06eaJfi/I+YTxiKYBEzwe6IjEWZxGhHpBl5X1e8H5bFly5aKfS0tLbS0tFTs7+/vp7+/\n39Jb+iGXvnP2bOZs3Mg1mzcP7Lt43DiOO/10+vv7GyZPcWVfSF/+DTdDeTZL+p0B67XHol5YEXAL\nniLIIjz0HOA+YGSV4xq09fT0BIdZ9fRYeks/ZNOfP2OGTgRtB50IOrrJ5Lf00dNPxAslXQHaU7RB\ntPDRwoIzVRGRDuCnwHNAwRGpqnpkzRPrICInATcA7ar6pyppdHNRC6eASxrZ0lt6S2/p80q/Zvly\nHrzqqhJLcG5bG/+6aROqKuV5VCOMItgEfA7YQFEfgapuCXuRKvluBEYAf/Z3/UpVP1OWRuvJZxiG\nMZRZtWQJy4r6ZabNmkX7qaemrgh+parHJ5Y2BqYIDMMwoiMiqSuCm4B9gDvxBn8BycNHw2CKwDAM\nIzpRFUGYqKGReH0DnWX7M1cEhmEYRvbUVAQiMgz4s6p+oUHyGIZhGA1mj1oHVfVN4AQRCW1iGIZh\nGM1FGNfQQ8B/icgPgVf8fQ3pIzAMwzCyJ2wfwZ+BqWX7TREYRhXynhLaMKJQVxGo6jkNkMMwBg0u\nTAltGFGo2UcAICKHiMh/isgf/e3HInJwI4QzjGak2kRgyxYtykkiw6hNXUUA3Iw3xcTb/e1Of19s\nRORKEXlYRB4SkbtF5JAk+RmGS7gwJbRhRCGMIthfVW9W1Tf87TvA3yS87nWq+l5VfR/wE6AnYX65\n0tfXl7cIoTA506WanC5MCV2g2cvSNZpFzqiEUQQviMhZIjJMRIaLyCeAwEniwqKqfyn6OTppfnnT\nLC+HyZku1eTsnD2b7ra2kn1z29qYlsNqXM1elq7RLHJGJUzU0D8Bi4AF/u9fAucmvbCIXA2chReS\nOjlpfobhCrYal9FshIka2gJ8OGrG9RalUdVuoFtE5gD/RgrKxTBcwVbjMpqJqpPOiUg1v70CqOoV\nqQggMg74maoeHnDMZpwzDMOIQVqTzr2MX+kXsTdwHrAfEFsRiMg7VXWj//MfgLVB6aLciGEYhhGP\nutNQA4jIWGA2nhL4AXCDqj4f+6IiPwLeBbwJbAI+nSQ/wzAMIz41FYGI7Iu3OtmZwHeBL6vq9gbJ\nZhiGYTSAquGjIvIl4AHgL8CRqtrTSCUgIieJyOMislFELm3UdeshIt8WkW0isr5o31tFZJmI/FZE\nlopI5YKjjZXxEBFZISKPiMgGEZntqJwjRWS1P7DwURH5VxflLOCHUK8VkTv9387JKSJbRGSdL+cD\nDsvZIiI/EpHH/Gc/yTU5ReRdfjkWthdFZLaDcl7mf+vrReT7IrJXVBlrdRbvwluR7I2Aw6qqY5Pf\nQhWhvHUQngA+CDwLrAFmqOpjWV0zLCLy98AO4LuqeoS/7zrgT6p6na+0WlV1To4yHgAcoKoPicho\n4EHgf+BFZjkjpy/rW1T1FREZDtwLXASc5pqcACLyeeAYYIyqnubac/dl3Awco6p/Ltrnopy3ACtV\n9dv+s98b6HZNzgIisgdeXXQcMAtH5BSRQ4HlwLtV9TURuR34GTAhkoyq6twGHA/cVfR7DjAnb7mK\n5DkUWF/0+3Hgbf7/BwCP5y1jmbw/wVOqzsoJvAVP4U9wUU7gYOAXwPuBO1197sBmYN+yfU7Jibf0\n7e8C9jslZ5lsncA9rskJvBWv0dyKF/xzJzAtqoxhRhbnwUHAM0W/t/r7XOVtqrrN/38b8LY8hSnG\nbzEcBazGQTlFZA8ReciXZ4WqPoKDcuKNdbkY2FW0z0U5FfiFiPxaRP7Z3+eanIcBfxSRm0XkNyLy\nDRHZG/fkLOYM4Fb/f2fkVM/yuwF4Gvg90K+qy4goo6uKoGnHD6ingp2Q33cL/Ri4UEun9XBGTlXd\npd6cUwcDU0Tk/WXHc5dTRE4FnlfVtUBgSLMLcvqcoKpHAScDF/iuzAEckXM4cDTwFVU9Gi9UvcRt\n4YicAIjICLxBtT8sP5a3nCLSBvwLnpfi7cBo8aYBGiCMjK4qgmeB4hlJD8GzClxlm++XR0QOBHIP\nhRWRPfGUwPdU9Sf+bufkLKCqLwJL8Hzwrsn5d8Bpvv/9VmCqiHwP9+REVf/g//0j8J94Pm3X5NwK\nbFXVNf7vH+Ephucck7PAycCDfpmCW+U5Efilqr6gqjvxFgw7nohl6aoi+DXwThE51NfG0/GmwnaV\nnwIz/f9n4vnkc0NEBPgW8KiqfrnokGty7leIZhCRUXi+zbU4JqeqzlXVQ1T1MDwXwXJVPQvH5BSR\nt4jIGP//vfH82utxTE5VfQ54RkT+1t/1QeARPP+2M3IWMYPdbiFwqzwfByaLyCj/u/8g8ChRyzLv\nTpganSAn43WCPAlclrc8RXLdiueLex2vH+NcvA6bXwC/BZYCLTnLeCKeL/shvIp1LXCSg3IeAfzG\nl3MdcLG/3yk5y2RuB37qopx4vveH/G1D4btxTU5fpvfiBQc8jNeK3cdROffGmx15TNE+p+QELsFT\npOuBW4A9o8oYamSxYRiGMXhx1TVkGIZhNAhTBIZhGEOc3BVB+bB9wzAMo7HkrgiAC/F6ua2zwjAM\nIwdyVQQicjDwIeCbVBmoYxiGYWRL3hZB0LB9wzAMo4HkpgjCDNs3DMMwsie3cQQi8r+Bs4CdwEhg\nLPBjVT27KI31GxiGYcRAIyz1m5tFoMHD9s8OSOf81tPTk7sMJqfJ2awympzpb1HJu4+gGGv9G4Zh\n5MDwvAUIWQQ0AAAel0lEQVQAUNWVwMq85TAMwxiKuGQRNC0dHR15ixAKkzNdmkHOZpARTM68Cd1Z\nLCIj8dY4eC1bkUquqXH8XYZhGEMZEUEjdBZXdQ35izX/D7y5uP8Oz3oQEXkT+BXwH8BPrKY2DMNo\nbqpaBCKyCrgHbxGGhwqWgIjshbcG7mnAiao6JTPhzCIwDMOITFSLoJYi2KueGyhMmhrnjsTrIN4L\nGAH8l6peVpbGFIFhGEZEoiqCqp3FqvqaiAwXkcdrpYkqYNG5rwLvV2/h8iOB94vIiXHzMwzDMOJR\nM3xUVXeKyBMiMl5Vn0r74qr6iv/vCGAY8Oe0r2EYhtGMLFmyihtvXMprrw1nr712Mnt2J6ecko0n\nPsw4grcCj4jIA8DL/j5V1dOSXtzvkP4N0AbcpKqPJs3TMAyj2VmyZBUXXvhzNm26emDfpk3dAJko\ng7rhoyLSEbRfVftSE0JkH+DnwJzifK2PwDCMoUhX1zyWLr0qYP987rrryrrnpxY+WiDNCr/GNV4U\nkSXARKDker29vQP/d3R0DNoBHYZhGAVeey24an711WGB+/v6+ujr64t9vTAWwfHAjcC78SJ8hgE7\nVHVs7Kt6+e4H7FTVfhEZhWcRXK6qdxelMYvAMIwhRy2LYNasaXX7DlK3CID/gzc76A/wWuxnA+8K\ne4EaHAjc4vcT7AF8r1gJGIZhDFVmz+5k06bukj6Ctra5TJ58cCZ9B2EsggdV9RgRWaeqR/r7HvLD\nPjPFLALDMIYqS5asYtGiZbz66jBGjnxzwBII03eQhUXwsj+a+GERuQ54DltRzDAMI1NOOWVKRSv/\n+uuXB6at1ncQljCK4Gw8181ngc8BBwMfS3RVo6lpZHyzYRi72WuvnYH7R458M1G+YaKGtojIW4AD\nVLU30dWMpqfR8c2GYeymWt/BrFknJco3TB/BacD1wF6qeqiIHIUX3ZN4QFld4ayPwDmSxjcbhpGM\noL6DRkQN9QKTgBUAqrpWRN4RQe6qiMghwHeBv8FbqvLrqnpjGnkb2RA1vjkPzHVlhKGR70ma1wrq\nO0hKGEXwhh/rX7xvV0rXfwP4nKo+JCKjgQdFZJmqPpZS/kbKZOWjTAtzXRlhiPuexKnQm+KdLKx6\nX20Dvg2cCawH3gksAr5a77w4G/AT4ANFv3Uws3jxSu3s7Nb29h7t7OzWxYtX5i1SXRYvXqltbXMV\ndGBra7vMGdk7O7tLZCtsXV3z8hbNcIg470nwuz+37rtf61pZ1QF+3Rm67g1jEcwCuoHXgFvxRgCn\n7gwWkUPxFrxZnXbeLtIUrYQACrItWjS/yEd5kjMyN4PrKk3MDRaPOO/JjTcuLfleATZtuppFi+bX\nLPNq19q69Xln6oAwUUMvA3P9LRN8t9CPgAtVdUfxsS1btlSkb2lpoaWlpWJ/f38//f39TZG+2kt1\n3XUXMWHCOKflP+GEIwNfVBfKv5rrSnXHwLvkWnnGTR/UmHjiiUvYtm0bU6ce67z8eabftesvFceh\n1MVZnv+LLwa/WwXlUU2eYcP+Gnje73+/ne3bv1ayb9Omq1mwYE4q31ck6pkMeNNJfANYhtdhvAJY\nHsXsqJP/nnhWxr8EHNOgraenJ9Ac6unpKUs7WmGijh8/M9Ds2p3eSwftChN1ypT/GWiuVeYfVZ7d\n6dvbewLNRU+G5PkP1fRB5jt83H/G7ssfJX01l4P3Lrsvf77pRyucXlJu5S7OyvwnBpZ3wZ1UTZ4Z\nM84PdKfuv/8ZgfmNHz8z8v2uWLFCe3p6BjaI5hoKEz66DrgJb92AgrpUVX2w5okhEK8H+hbgBVX9\nXMBx3bx5c8V5YTT+8uVruOKKNTz11HUDx9vaulm4sGtA2/b393PHHcvK0q1i+PDvs3PnVyvOO+GE\nI1NrsUyf/qXAMMwpUy7ills+mzj/oZy+EF63Y4cybNirzJw5aVC2kDs6elm5srfi+KRJc7jttk85\nL3/e6ZcvX8Mtt6zmzTdHMnq0VIRhBqWvrFPmsnCh5xq9/faf8e//fjevv74nI0a8wTnnTGbq1GNp\naWnhvvvWVYR8LliwhOXLr62Qc+rUOdx99zWJ7jdq+GiYFvuDUTRLlA04ES8C6SFgrb+dVHQ8UDOG\nIWxnUGW6eJ2NUTt9XO90NdzHOsYbz+LFK7Wra562t/cMdPYW9kftSM6yDiCiRRCms/hOEbkAuAOv\nw7igQBIvK6mq91Jj3eQkhO0MqkwXvRMpTsdv1p2ueXUiWudl48hqlKlRnWox/FE6kou/kbFjt3H0\n0RcwZsz+uQZehFEE5+D5oi4q239Y6tKkSNh498p00ePk40YTZDEwBPKLSGrWSKhmxfUIrrRxuZER\ntuEZ9I20tXVzxRVTU7mXQhlFJor50OiNBK6hsGZXZbqVOnz4JyOZa9U6ftvbe2LLn4S8XAbmqjCy\nIm4Mf6OI74pO7xspLaOUXEMi8gFVvVtEPoZnEZQrkDuiq53GEba1FJRu8uQjuf/+8K2sNEbbptna\nySuWfqjF8BvZU/gu1qx5ku3bbys5FsbqbhRh3XRZfiNBnomw1HINTQHuBj5MgCLA6zNIhIh8GzgF\neF5Vj0iaXzlhXS9JXTRJfbVpu1TymgbC9eknhjIuu1WqUfpd9AamcaWREbbhmeU3Uk3JhCKK+ZD2\nBvw93mji9VWOJzaXGkW1aIIwpG0u5hWRZJFQbuK6W6Uapd/F4HA7ZvmNlJZXeq6hLwTpDbzVyVRV\nF8RXPwNK6B5/aommJ4lVkba5mFcn4lDrvGwW4gYz5EXBelm9emvR3k68mW6aO0Iqy28kyDMRllq2\nxBiCXUJO04wmcBbmYlYRSa5e16hOM/XdlLqD5hUdKbxT82ltfZrjjhuXaSMjaT1S6/ysvpFiJfPz\nn0c8OYr5kMUGHEpKrqFmNYHNpWJkSTNFc5XKulLBFRdn+HrEhXqIFF1DvcBNqrqtyvEDgU+pak9E\n3ROJ3t7egf87Ojro6OiomjbuoI68LQdzqRhZ0kwDz0qtl91WwD77PMPkyYck/i7CfPdJXWl5uOL6\n+vro6+uLfX4t19CvgdtEZATePEN/wOsfOAA4Gm+U8ZdiXzkkxYqgHkkGdcSJ0nF91SHDgOZqaFS6\nSacAU5g8OflSqGG/+6SutDxcceWN5MsvvzzS+VUVgaouBhb7y0meABTmRr4XuFZVt1Y7NywicivQ\nDuwrIs8AX1TVm+PmF9bXnobGbuQoWpesF6M5aZaGRhbWS9SxCEn77JoyjDqKH6nRG6n0EVT6FNMY\nCdwov6sL/sahQDOuFjdYSRKKHZTX7u8n3HeftM/OhT4/0uojaEYaOagjihsqSWu+2UL/kpKH9ePK\nHElm+Xmkab2Ufj/hvvukrrRmcsUNEEVrNHojowFlaWjsMBZBGq151+YxypK8rB8XomrM8suG0u8n\nnyikPGCoWARJWk9paOwwvsw0WvNN6W+MSV7Wjwtx9kPN8msUpd9PvLEIQ8FSi6UIROTDqnpn2sKE\nJQ1TPqn5GUaZpFHBROk8a/YXNq8K2QVl64IyGoxUfj9TaGu7i4ULzwv1bTTabZjXNxzXIpgIJFYE\nInIS8GVgGPBNVa1cty2AytbTKjZtEs4661sce+zSqoWXdiHXUyZpVDBhrRdX/NxJyKtCdiHO3gVl\n1GyE+Z6TWv+NtNRy/Yaj+JHS3PAq/yfxRhbvibdc5bvL0gT6v+r7/Sp9q3n4YBsZPeCCnzspeUZb\npBmpEvf6eUeaNBON+p4b2UeX5jdM2n0EInI6cJeqviQi8/EGk12pqr9JqIOOA55U1S3+dW4D/gF4\nrN6Jpa2npRRPRAXBGjsPH2wjowdquRaaxWWUZ7RF3nH2TRlpkiON+p4baanl6R4M4xqar6o/EJET\ngQ/gjSb+Kl5FnoSDgGeKfm8FJoU5sdSUj7s2cXC6tGlUBVPthX3ppa1N5TLKu0LOk7TvvVkaAHFo\n1Pdcz22YZhnXUjq1rhN0LCphFEFB9Z0KfENVF4tIsrHeHhom0ZYtWyr2nXDCkSxc6LWeHnhgI9u3\nB2SuOwbObWlpqVrIxekKtLS00NLSUpG2v7+f/v7+iv0upJ8+/b088cQlPPXUdQPH29rmAiMCW07X\nXXcREyaMc0Z+S59u+uXL13DFFWtK3ofyBkBa8qxe/Sjf/vYvKyqpLO+31vfc39+fWnlOmDCOuXPf\nxy23XDRwf5/5jLe+cJBP/4knLmHbtm1MnXps5Put9g1Pnnxw1esAXHXVg2zefM3AscJzjkQ93xGw\nBPg6sBloAUYCD0fxP1XJdzKey6nw+zLg0rI0GrT19PQM+MKCfIXwcYXRJemD0rW2frIkXVD+xfT0\n9NSVJ9/0oxUm6vjxMwf83ME+zpUK0xTaFSbqjBnnZy5P0GjdqPnPmHG+wsQBuQvPLi953E4/sa6/\nOR15Rmtr66dKrlHw1Wd5v4sXr6y4buG7b1T5V/Ppe2Ufvzzb2s4o6auqfZ3Cc16h3sjpwhatj0DU\nq3CrIiJ7A114U0Vv9GcdPUJVl9Y8sQ4iMhx4As/d9HvgAWCGqj5WlEY3b95ccW65hl2yZBWLFi1j\nxw5l2LBXmTlzUqBGLqQr+GDPPfd4Jk16T938C6TRwlm+fA3/8R+/5s03R1WYeOUtuu98535ef31P\n9t5b+fznT6kwOfv7+7njjmUD6UaMeINzzpnMRz86bUCerq55LF16VdFZq4CfU9yvcthhc1i06EOB\n+afZIm1r62bhwq5YLdIlS1Yxa9bPSlo+48dfwhe/eGzJ/aYhT3HZjxjxBhdc8AGmT/9Qzfyjlk/W\n6c8446usXn1NRZr29l76+npTk+fssxdxzz03VKTt6prPbbd9oWr+9923rsKdccIJR0aS5/bbf8ZX\nvrJ8II/Cd9+o8u/o6GXlyt6K45MmzeG22z6VOP8Cta4DBD5nEFRVAg4EE0Zb4C0pea7///7AO6Jo\nmxr5noynDJ4ELgs4HqhJm5WwkQ5ppqtMU966WKnQra2tZ6c6x07aUUxJ8wt7/mAZ4duoKLI4UTVW\nxuldp7q1EM0iCFNZ9+KNGfit//sg4L4oF4m7DTZFEPbFSTtdcWhka+vZZUogmw+ysoLwFM4++wS7\nZbLOL2yFNRjCcFUbF44ap7xcLeOoEw82qoxrXafasaiKIExn8UfwFph/0K+ZnxWRMaFNDmOAsJEO\naacrjkbxXEWFI+FCb+NQ2pm32x314ouwdGn0yKWk+YUNAxwsI3zTCEcNExETZzCei2UcZzBXo0J+\nw1yn/Nipp/5rtIvU0xTAA/7ftf7fvYF1UbRN3A3HLIKkUxXnZRGU30PUaXnjUHqd5C3ApPmFbb25\n2lptNFHcN1EH47lYxi7KlAQycA1dDHwNL2rofOB+YHaUi8TdXFAEhcp/woTzddSoTyZyo4StjNJO\nFyRHV9c8bW2dHvrlj6MEC9fZZ5+ZqSicpPmFqbCG+gjfwnOO8m7EuUaWZRznXR1ss/ymrgi8POnE\nG0j2JWBalAsk2fJWBHFbobVexLCtp7TT1b+/qAonvBJ0reO4HnHKtJkXtglu7GRbMWY1pUfcd9Us\ngvqV8f4wEGYqwD8CG6JcJCDPjwOP4A1WO7pGuqzKKRSlL0e4DyPriIi0K5wwH2TSj6SewnG1ky4s\nzRwFU72x05wVY9x31bV3KilRFUHVzmIR+SjeQLI3gDdF5DN4EURPA2dH64moYD1eJ/TXEuaTKaWd\nWo1bD7kaWcxOGGZag6Sde7U6u+rdU60OS1fm5WnmtQRKZS9+zp1AN8XBBI2ekTUOcd9V196pRlMr\nauhyYLKqPikixwCrgY9oCusQqOrjACLhxzvkQWmkSbgPI8uIiLwqnLSm0w6SsdY9ATWVRB4faZBi\nyuKZN2qeoOqNnXiLuORNknd1KM9zVUsR7FTVJwFU9UEReTwNJdBMlIbGeS/IqFHTaWs7kIMOGpPZ\nesjVyCvsLsv5+mvdUxqKL80KtZr1MnZswGRXxH/mjZyXvnZjJ9oiLi7gwtoSzUgtRbC/iHwer18A\noKXot6rqgloZi8gy4ICAQ3ObRaEEm4sX1PwosnwR81q8JEuzudY9vfpqMsWXdoVaTTEdddT/oq0t\nvWdeTwGmqdziNHZcZqi7eOJSSxF8ExhT43dNVHVaXKGK6e3tHfi/o6ODjo6ONLINTVRzMcsXMc/W\nTlZmc617uvHG4Omswiq+tF1p1ayXsWMP5sorp6b2zOutLZGmcovT2HGdoeji6evro6+vL/b5VRWB\nqvbGzjUaNTsKihVBs5DViziYWjvFrdqxY7dx9NEXMGbM/hX3lETxpe1Kq2W9pPnMg6+zig0bHuOs\nsx5n+/bbSo6k4S67667aM8u7vraB6/JlTXkj+fLLL490ftw1ixMhIh8BbgT2A5aIyFpVPTkPWZqN\nLFs7jfqYglq1bW3dXHHF1JLrJVV8abvS0rDIysv4+OPfzq9+9fuSMq+8ziqGD/8+L7xwO17gXiVZ\nussa0WeR5N0bDOt1506UWNNGb+Q8jmAo0chY+EYN3skiNjzdAXwrdfjw4NHqxdfZd9/TU4vvd3GS\nONcGLA4GSGscgTG0aGRoaqOin7JwpSWxyCrLeCk7d361JE2hzO+668qB63jz0RdSJIvvj1P2WT+v\npO+ei5PYNRuxFIGIHKOqD6YtjJEfcT+mOCZ92i6bWjK41HFYWcbhyry0vJLF98cp+6yj1ZJW5HlF\n0w0qopgPhQ1v7WJzDQ0igs3rlbrvvqdXnfohrkmfpsummaZ3qCzjJIvlpFletfPKevqFrKcwGYoQ\n0TWUeWVe9cJwPfAY8DBwB7BPQJosysgIIIr/ukCSDzitSceayT8croyDK7A0J2mLO6leFpPEFfJO\nWpFnKV8zElUR1F2zGEBE3gscym5bVlX1jiSWiIhMA+5W1V0ico2f6ZyyNBpGPiMditd03rDhMT9K\npZSurvkDoYbV1lItXhc3a5LIENatlfbo5OJ1sydPPpD77/9DUR/GNGdcWY2kvFyGajmkhUi0NYvr\n9hGIyM3AEXizhe4qOpRIEajqsqKfq4GPJcnPSE6xP720g3I3xX5bF3yzcWUIG3KYRmjiYIxxT1uJ\nutSXMySpZzIAj+JPQ53Vhrcm8j8G7E/PVjIiEcbl4oJvNq4MWa4CV18+N/swwhL2ngbjvTcLZBA+\nugZ4D55FEIkw8w2JSDfwuqp+PyiPLVu2VOxraWmhpaWlYn9/fz/9/f2WPoX0s2d3snHjHDZvvmZg\n37hxF3P66cfR399PS0tLRXjmsGGvcuaZxzBhwriS55al/BMmjOPqq0/g5psrQ0Rr5V8tUqW//40S\n2V98MdjiePXVYaHkrxYaed11FzFhwrjI9+tC+muv/S82bbqh4p7Kwz0XLFjCpk3XFqVaxaZNwpln\nfpNJk5ZWWAeu3m+zpo9EPU0BdAAvAb/FW0dgPSmtWQycA9wHjKxyXIO2np6eQC3Y09Nj6VNMP2PG\n+QoTFdr9v6ObSv5a6au19L37LE4/sapFEKZ8qi2B6J3jbvnUTt8eeE/lizSNHz+z6PhKhdrWgbv3\n6376FStWaE9Pz8AG0SyCup3FIrIJ+BywgaI+AlXdUvPEOojIScANQLuq/qlKGt28eXPF/mbRyJbe\n3fT33beuwvd/6KGXMn/+RKZOPXZg3/Lla7jqqgdLLKO2trl84hMH893vPl2yf/z4S/jiF4/lox+d\nNiBPV9c8li69qkKGKVMu4pZbPtuw+00z/dlnL+Kee26oON7VNZ9Zs6YN9AmsW/co27f/wD86D6gs\nh+LgA1fvtxnTR+0sDtNq/1UUzRJ2AzYCTwFr/e0rAWkCNaNhpEGSdaHD9h240I+SNtXuqafn32uE\nxwZbRs26OLzrENEiCNNH8JCIfB+vQ/f13fojWfioqr4zyfmGkZSwkSpB6a6/fnlg2vLRsINpxtgC\n1e6psj9kCjt3wr77nsGuXbA9YP0eG/3rBmEUwUjgNbxJTopJpAgMo5mJErY6GEMjwyvHKRx++HIu\nvngqF15oK4e5Sk1FICLDgD+r6hcaJI9hNAW2JGIl9dZrgMFlGQ0mwnQW3w8cr/USZoCNLDZcxkbD\nlhK8zsRcFi60Cr/RRO0sDqMIvgq8Hfgh8Iq/O3EfQRhMERhGeFwYwWzK0Q2yUATf8f8tSaiq50aW\nbneeVwKn+Xm+AJyjqs8EpDNFYBghqLbq28KFXVYRD0FSVwRZICJjVPUv/v+zgPeq6v8KSGeKwDBC\nUG28QnGcvjF0iKoI9giR4SEi8p8i8kd/+7GIHJxEyIIS8BkNBA4oMwwjHLZKl5GEuooAuBn4KV4/\nwdvxxhPcnPTCInK1iDwNzASuqZfeMIzquDATrNG8hFEE+6vqzar6hr99B/ibeieJyDIRWR+wfRhA\nVbtVdRzwHeDfktyEYQx1Zs/upK2tu2SfF846LSeJjGYizICyF0TkLOD7gABnEMKVo6ph38DvAz+r\ndrC3t3fg/46ODjo6OkJmaxhDB4vTH9r09fXR19cX+/wwUUOHAouAyf6uXwKzVPXp2BcVeaeqbvT/\nnwUcp6pnBaSzzmLDMIyINEvU0I+AdwFvApuAT6vq8wHpTBEYhmFEJDVFICI9Vc4pTAt6RXTxomGK\nwDAMIzpprln8MmWDyIC9gfOA/YDMFYFhGIaRPaFcQyIyFpiNpwR+ANwQ5MpJG7MIDMMwopOmRYCI\n7Iu3OtmZwHeBo1U1YFZxwzAMo1mpqghE5EvAR4CvA0eWjQY2DKMGLkwAZxhhqdVZvAtvRbI3Ag6r\nqo5NfHGRLwDXA/up6p8DjptryGg6bAI4I29Sm2tIVfdQ1ZGqOiZgS0MJHAJMw1u3uKlJMpCjkZic\n6VJNzsolG2HTpqtZtGhZA6QqpdnL0jWaRc6ohJliIisWAJfkeP3UaJaXw+RMl2pyujQBXLOXpWs0\ni5xRyUURiMg/AFtVdV0e1zeMLLEJ4IxmI8xcQ7EQkWXAAQGHuoHLgM7i5FnJYRiNxtYzNpqNhk8x\nISKHA3eze9nLg4Fn8eYber4srfUUG4ZhxMD5uYZKBBDZDBwTFDVkGIZhZE+encUFrNVvGIaRI7lb\nBIZhGEa+uGARBCIiJ4nI4yKyUUQuzVueAiLybRHZJiLri/a91V+R7bcislREWnKW8RARWSEij4jI\nBhGZ7aicI0VktYg8JCKPisi/uihnAREZJiJrReRO/7dzcorIFhFZ58v5gMNytojIj0TkMf/ZT3JN\nThF5l1+Ohe1FEZntoJyX+d/6ehH5vojsFVVGJxWBiAwD/g9wEvAeYIaIvDtfqQa4GU+uYuYAy1T1\nb/E6wuc0XKpS3gA+p6oT8BYUusAvP6fkVNVXgfer6vuAI4H3i8iJOCZnERcCj7LbneminAp0qOpR\nqnqcv89FORcCP1PVd+M9+8dxTE5VfcIvx6OAY/ACXP4Th+T0Fw77Z7x54I4AhuGtIhlNRlV1bgOO\nB+4q+j0HmJO3XEXyHAqsL/r9OPA2//8DgMfzlrFM3p8AH3RZTuAtwBpggoty4kW3/QJ4P3Cnq88d\n2AzsW7bPKTmBfYDfBex3Ss4y2TqBe1yTE3gr8ATQijcc4E68GRsiyeikRQAcBDxT9Hurv89V3qaq\n2/z/twFvy1OYYvwWw1HAahyUU0T2EJGHfHlWqOojOCgn8G/AxcCuon0uyqnAL0Tk1yLyz/4+1+Q8\nDPijiNwsIr8RkW+IyN64J2cxZwC3+v87I6d60ZY3AE8Dvwf6VXUZEWV0VRE0bQ+2eirYCflFZDTw\nY+BCLZs91hU5VXWXeq6hg4EpIvL+suO5yykipwLPq+paqgx+dEFOnxPUc2WcjOcS/Pvig47IORw4\nGviKqh6NtwhWievCETkBEJERwIeBH5Yfy1tOEWkD/gXPS/F2YLSIfKI4TRgZXVUEzwKHFP0+BM8q\ncJVtInIAgIgcCGS+aE89RGRPPCXwPVX9ib/bOTkLqOqLwBI8X6xrcv4dcJo/5uVWYKqIfA/35ERV\n/+D//SOeP/s43JNzK94UM2v83z/CUwzPOSZngZOBB/0yBbfKcyLwS1V9QVV3AnfgudYjlaWriuDX\nwDtF5FBfG08HfpqzTLX4KTDT/38mnk8+N0REgG8Bj6rql4sOuSbnfoVoBhEZhefbXItjcqrqXFU9\nRFUPw3MRLFfVs3BMThF5i4iM8f/fG8+vvR7H5FTV54BnRORv/V0fBB7B8287I2cRM9jtFgK3yvNx\nYLKIjPK/+w/iBTREK8u8O2FqdIKcjNcJ8iRwWd7yFMl1K54v7nW8foxz8TpsfgH8FlgKtOQs44l4\nvuyH8CrWtXiRTq7JeQTwG1/OdcDF/n6n5CyTuR34qYty4vneH/K3DYXvxjU5fZneixcc8DBeK3Yf\nR+XcG/gTMKZon1Ny4s3i/Aie0r8F2DOqjDagzDAMY4jjqmvIMAzDaBCmCAzDMIY4pggMwzCGOKYI\nDMMwhjimCAzDMIY4pggMwzCGOKYIDMMwhjimCAzDMIY4/x/GQg8WsOslJwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 394 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we clearly\\* see some structure in the residuals from the second fit, whereas the quadratic fit has no structure and most of the points are within +/- 2 error bars.\n", "\n", "\\* Usually - though the simulation uses random numbers so occasionally it's not so clear...\n", "\n", "We can go further by plotting histograms of the normalised residuals. If the function fits well and the central limit theorem applies then the normalised residuals should follow a gaussian (normal) distribution, which we can also fit to using an appropriate function and *curve_fit*. However, this is left as an excercise for the reader in this particular case..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## More information and further reading\n", "\n", "Another example in the Lineshape Comparison and Analysis tutorial.\n", "\n", "Wikipedia entry on the covariance matrix\n", "\n", "Wolfram Mathworld entry on Covariance\n", "\n", "Chapter 6 of *Measurements and their uncertainties*, Hughes and Hase, OUP 2010\n", "\n", "Official scipy curve_fit documentation\n" ] } ], "metadata": {} } ] }