{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting problem 3\n", "\n", "This is a linear model. (It is not a linear function of $x$, but the model\n", "is linear in the fit parameters $a_1$ and $a_2$.)\n", "\n", "The data for this problem is available in the file `fit_3.dat` which can be downloaded from http://www.eg.bucknell.edu/~phys310/hw/hw4.html. One option is to download the data file to a local directory, and then import it using `np.loadtxt()`. Another option is to look download it directly into the Jupyter notebook. I will use the second option here." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import optimize\n", "\n", "import urllib # for importing data from URL\n", "\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Following is an Ipython magic command that puts figures in notebook.\n", "#%matplotlib notebook\n", "#for making pdf-file of this notebook, use instead following line\n", "%matplotlib inline\n", " \n", "# M.L. modification of matplotlib defaults\n", "# Changes can also be put in matplotlibrc file, \n", "# or effected using mpl.rcParams[]\n", "plt.style.use('classic')\n", "plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures\n", "plt.rc('axes', labelsize=16, titlesize=14)\n", "plt.rc('figure', autolayout = True) # Adjusts supblot params for new size" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Options for getting data into notebook\n", "+ Start by downloading the data files into the working directory. The `np.loadtxt` \n", "function imports the content of the data file into a `numpy` array.\n", "The `unpack = 'True'` option transposes the array so that each column is in \n", "a separate array.\n", "+ Download directly from a URL" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def f(x, a1, a2):\n", " '''Linear model function with adjustable parameters a1 and a2'''\n", " return a1*np.sin(2.*np.pi*x)+ a2*np.sin(4.*np.pi*x)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "g = urllib.request.urlopen('https://www.eg.bucknell.edu/~phys310/hw/assignments/fitting_3/fit_3.dat')\n", "data = np.loadtxt(g)\n", "\n", "x, y, u = data.T" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "#Uncomment to look at the data.\n", "#data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### (a) Is the suspected relationship relationship between $x$ and $y$ a linear model, or a nonlinear model?\n", "\n", "\\begin{equation}\n", " y=f(x)=a_1 \\, \\sin(2 \\pi x) + a_2 \\, \\sin(4 \\pi x)\n", "\\end{equation}\n", "This is a linear model. (It is not a linear function of 𝑥 , but the model is linear in the fit parameters $a_1$ and $a_2$.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### (b) Perform a fit to these data using the generic optimize.curve_fit routine from scipy that we demonstrated in class, using the assumed functional form shown above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Perform fit" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.96501629 0.87923822]\n", "[[6.39999997e-03 1.10687288e-11]\n", " [1.10687288e-11 6.39999999e-03]]\n", "0.07999999982280558\n", "0.07999999993543226\n" ] } ], "source": [ "popt, pcov = optimize.curve_fit(f, x, y, sigma=u, absolute_sigma=True)\n", "print(popt) # best fit parameters\n", "print(pcov) # covariance matrix\n", "for i in range(len(pcov)):\n", " print(np.sqrt(pcov[i,i])) # uncertanties in fit parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Results\n", "$a_1 = 1.96 \\pm 0.08\\quad\\quad\\mbox{and}\\quad\\quad a_2 = 0.88 \\pm 0.08$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFbCAYAAACOHWQYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAxOAAAMTgF/d4wjAABG+0lEQVR4nO3dd3iN5xsH8O85GRKyJEKEBBVixogiqB2htNQqpWaL1Cgxa9RoqZWYkRq19yaCtDRoNT+rRltKY0VJrJCThOQkOef3R5oQTnL2ec+b8/1clwvvvPNeSe7zPO/z3I8kJiZGCSIiItKZVOgAiIiIxI7JlIiISE9MpkRERHpiMiUiItITkykREZGemEyJiIj0xGRKRESkJyZTIiIiPYk6mU6dOhWtWrXChQsXhA6FiIgsmGiT6ZEjR5CRkSF0GEREROJMpomJiVi/fj0mTJggdChERETiS6YKhQJz587FgAED4O7uLnQ4RERE4kumu3fvhr29PTp06CB0KERERAAAa6ED0Mbdu3exc+dOfP/99xodr1Ao8PTpU9jb20MikRg5OiIiEhOlUomXL1/Czc0NUql+bUuJmJZgO3r0KObPn58vMSoUCkilUrRq1QpTp07Nd/zjx4/Rs2dPU4dJREQisnPnTr1fG4qqZdqsWTP4+vrm2zZo0CCEhISgYcOGbx1vb28PALh69SrKlStnkhiLismTJ2POnDlChyFKfHa64XPTHZ+dbmQyGby8vPJyhT5ElUwdHBzg4ODw1nYPDw+VnypyW7COjo5wcnIyenxFia2tLZ+ZjvjsdMPnpjs+O/0Y4jWg6AYgERERmRtRtUxViYmJETqEIikoKEjoEESLz043fG6647MTHlumpBJ/OHXHZ6cbPjfd8dkJj8mUiIhIT0ymREREemIyJSIi0hOTKRERkZ6YTImIiPTEZEpERKQnJlMiIiI9MZkSERHpicmUiIhIT0ymREREehJ9bV5LEx2d8yc7Gzh7FmjYELCyAoKCcv4QEZHpMZmKTG7SlMkAZ+ecxMqVl4iIhMVuXiIiIj0xmRIREemJyZSIiEhPTKZERER6YjIlIiLSE5MpERGRnphMiYiI9MRkSkREpCcmUyIiIj0xmRIREemJyVSk5HI5gIT//iYiIiGxNq/IKBQKjBsXisjIeABeCAi4hw8+8MbChWMhlfKzERGREJhMRWbcuFBERLyL9PTxAIC4OCAi4gSAUISFjRc0NiIiS8WmjIjI5XJERcUjPb1lvu3p6S0RFXWXXb5ERAJhMhWRp0+fQibzUrkvJcULSUlJJo6IiIgAJlNRcXNzg5PTPZX7HB3/haurq4kjIiIiQITvTLdu3YqjR4/i0aNHKFasGGrVqoVhw4bBy0t1i60osbW1RceO3oiPP5Gvq9fO7gQ6dvSGra2tcMEREVkw0SVTT09PfPnll/D09ERaWho2bNiAr776Cps3bxY6NJNYuHAsgFBERu5GXJwXfHz+zRvNS0REwhBdMm3ZsmW+/w8cOBCDBw9GUlKSRXRzSqVShIWNx+TJcri7JyE21hWlSrFFSkQkJNEl09dlZGTg6NGj8PLygouLi9DhmFROl64H2LNLRCQ8USbT2NhYzJo1CxkZGShfvjzmzZvHggVERCQYUWagunXrYs2aNViyZAkqVKiAb775BllZWUKHRUREFkqULVN7e3uUK1cO5cqVQ7Vq1fDhhx/izJkzaNq0qcrjZ86cCUdHRwBAUFAQgoKCTBmuVqKjc/5kZwNnzwINGwJWVkBQUM4fIiLSXXR0NKKjowHAoIVuRJlM36RUKmFlZVXg/unTp6N8+fImjEh3uUlTJgOcnXMSq5OT0FERERUNrzeoZDIZwsPDDXJd0SXTlStXolmzZnBzc8OzZ8+wbds2ODs7o1atWkKHRkREFkp0yfTRo0eYMWMGkpOT4ezsDD8/P4SGhsLBwUHo0IiIyEKJLplOmzZN6BCIiIjyEeVoXiIiInMiupYpGY+6kcQcaUxEpBqTKeVRN5KYI42JiFRjMhWZ11uHjRsD06Zp3jpky5KIyDiYTEVGn8THliURkXFwABIREZGemEyJiIj0xGRqpnJqRiYYtHYkEREZB9+ZmhmFQoFx40IRGRkPwAsBAffwwQfeWLhwLJeZIyIyU0ymZmbcuFBERLyL9PTxAIC4OCAi4gSAUISFjRc0NiIiUo1NHTMil8sRFRWP9PSW+banp7dEVNRddvkSEZkpJlMz8vTpU8hkXir3paR4ISkpycQRERGRJphMzYibmxucnO6p3Ofo+C9cXV1NHBEREWmCydSM2NraomNHb9jZnci33c7uBDp29Iatra1J4lA3kpgjjYmI8uMAJDOzcOFYAKGIjNyNuDgv+Pj8mzea19jUjSTmSGMiItWYTM2MVCpFWNh4TJ4sh7t7EmJjXVGqlGFbpDktyqeQy90AvLq2upHEHGlMRKQak6mZyunS9YC2PbtKpRKXEi8h5k4M4pPjkZCagAcpD/Ds5TO42rkica8Ejy6WAWyqod67CejepQpCQ8chKyvrv5HE+ZNizkji3Zg1K7XQ/XPnyk3WDU1EZG6YTE3MGCu3KJQKHPnnCA5cP4DD/xxGckYyWldqjcolK6NB2QYoW7UsStqVxLJvduH2T32QJQ8EAPx7B1i8/AiO/NMaI8Z0h0xWXuX1U1K8cPPmTbUjjT08PHT7AoiIRI7J1MQMuXJLtiIbO//aiW9/+RbJ6cnoUaMH1ndZj+YVmsPWKn8rUS6XI+Ts0bxEmierAx5f3oBll5bhcVZlABPfuo+j47+oXLkynJzuITHx7Tg40piILB2TqQgplUrs/GsnZpycgReZLzC52WQMqDsAxayLFXhOYXNYi2XXQ8wnYRj+57c4uCUaisxXTeTckcQODg7o2NEb8fEn8hWVMPVIYyIic8RkKjKP0h5hSOQQnLl/BrNazkL/uv3faoWqkjuHtaCWZSm3Utjzw3KMKxmK3ft3416CK6QOV9AsyB3zF6wFoH6kMRcfJyJLxWQqIpHXI/FZ5Gd4z/s9/Bn8J9yKu2l8bu4cVnUty1cjiZ9i/cnTmPbbODRb1wzfd/oedT3qFjrSmIuPE5GlYjIVgczsTIw6Mgpb/9yK5R2Wo69fX0gkEq2vo+kc1pzEWhadq3VHt3rvY/ap2Wi6tinmtpmLEQ1H6DzSmC1XIiqqmEzNXKo8FT0PdkdiaiKuDLuCCi4VdL6WLnNYi9sUx+w2s9GhSgf02t0LJ+6ewKKWPwBw0fr+bLkSUVHFsjXmrMQjdNrdCpmKTJwccFKvRPq6Vy1LzZuWzbyb4dKwS0jPSsd7W+oBnucNEgsRUVHAZGqmbj2/CQxugndcfHD4k8NwtnMWOiSUKl4Kkb0jMchvGDCgJaJuHhA6JCIis8BkaobuPL+D93c2B65/iDUdthQ65cXUpBIpxrw7Edi3EZ8f6YPws+FCh0REJDi+MzUzOw49wrCz7eCd2RXeyaGY/rVE40E6Jh3gc60r9nX1QO+DHyI+OR7ftf0OUgk/mxGRZWIyNSOyDBnmP+iADnUaYHPXJZBqOWLX1AN8Gnk2wW+Df0OHLR2QmJaItR+uhZXUyng3JCIyU2xKCOTNNUHTs9LRZXsXuBd3x/ou60XTyqvqVhWnB53GufvnMODAAGQrsoUOiYjI5ETXMt28eTNOnTqFe/fuoXjx4mjYsCGGDh0KFxcXoUPTiKo1QTt18sKDgAt4kfkCkb0jNapoZE48HDwQ0z8GrTe2Rr/9/bCs9QaI8FuLiEhnovuN9+eff6JHjx7w9fVFWloali5dilmzZiEsLEzo0DSiak3Q5St+gt3va3Dr8GmUsC0hcIS6KeNQBj/3+xmtN7bG0KP9AOlGiPDbi4hIJ+LoS3zN3LlzERgYCG9vb1SvXh0jRozAxYsXkZqaKnRoasnl8v/WBG2Zb3uWPBCu/74HZxvhp7/oo4xDGcT0j8GfiZeAoF5Iz0gXOiQiIpMQfdMhOTkZtra2sLe3FzoUtQpbuSXzZVWzXxP09dHCjRsD06blHy2sUCgwd9oGpB9oASSWQNW6/TGwpz9CQ8dBKhXd5zYiIo2JOpnK5XJs3LgRQUFBsLIy/1Gk6lZuEXpNUHXJUt0Umze7sJNf5HRhSyShCAsbn3dczqCrp5DL3QCI6/0wEZEqok2m2dnZmDNnDgAgODi40GNnzpwJR0dHAEBQUBCCBKqqrunKLcaib7IszKsu7PH5tmfJA7F1z0bMnSuHtbX1W4Ovcgvts+VKRKYQHR2N6OhoAMibTWEIkpiYGKXBrmYiCoUCc+fORVxcHBYvXgynAiZTpqWloVOnTrh37x7Kly9v4ihVUygU+PjzL7HnwA0ok5rDp3KiwRNK7jzT5GTTFZJPSEhA/fqbkJg44a19EofR+OFwHfyx78l/LdeWefvs7E4gOPhcvpYrEZEpyGQyODs749ChQyhRQr/Bn6JrDiiVSixYsABXr17FwoULC0yk5io1MxXnakfiux3vA8rBiI3N6QIVe8sstwtblbKlnmPEsRHYc+DGW4Ov0tNbIirqrkE/IRIRmZrofoOHhYUhNjYWU6ZMAQAkJSUhKSkJ2dniKBYw5ugYVHGrgqENRkLblVvMWW4Xtp3diXzb7exO4OOPamKs/1j8+1j1B5+UFC8kJSWZIEoiIuMQ3TvTQ4cOAQC++OKLfNu3bdtm1iNhASDyeiT2XNuDP4L/MGqFI6EG+BS2+HhmZiaWOvZFcsrb55nD4CsiIn2ILpnGxMQIHYJOnrx4gs8jP8eS9kvg5ewFmczw91BVXcmUA3wKW3y8WLFi6N/DH8tXREOR+WqU0+uDr0xaqJ+IyIBEl0zNXUEJ4VLl0WhUvhH61elntHurqq4UEXECQKhJB/i8Wnw8//ZFYROQpZiD1duXIjO5Dny8U/KSPWD6Qv1ERIbCZGpgqhLC+ac/Y832g7j26TVItFwJRlMFTU3JGeCzG3PnygV/PyuVShG+dCoGjLqGhosD8GXPUIxoPljQmIiIDEF0A5DERp4txxdRX2Bmy5ko51TOaPcprLqSuQ3w8S1dHTiyGxN/GYUz/54ROhwiIr2xZfoGQ7+3W3phIYpZF8PIRiMNH+xrzL260ltutcW0prPx0Y6PcH7IeXg6egodERGRzphM32DQ93YutxF6ZjZ+6vcTrKXGfdRCV1fSRXC9L/H380vouqMrTg44iWLWxYQOiYhIJ0ymRqJUKoEOo9C92ido4tXEJPcsbGqKOZJIJPi+0/dosb4FRh4ZiVUfrBIsFo4kJiJ9MJkayeFbBwGvWMxsdt1k9yxsaoq5srO2w+4eu1F/VX00KtcIg+sLMyCJI4mJSB8cgGQE8mw5pp4aBxyfA1d7N5Pf/9XUFPNOpLm8nL2wo/sOjDo6CucfnBc6HCIirbFlagQR5yJQzKoYcHHQW/vUrdxiqVpXao3pLaaj285uONH7AoBSQodERKQxJlMDe/byGWadmoXV7begm+Ltx2vpSbMw45uMx5n7ZzDocG9AchSA+a9RS0QEMJlqTd1Aldm/zIZ/WX+0qcCMqS2JRIJ1ndehwcqGQMsZAL4ROiQiIo0wmWqpsIEqt57dQvi5cJz57IzRKh2ZM0N0YTsVc8LadtvxXlxTRN1ogN4NOhs3aCIiA2AyNaBJxyahT+0+8CvjhydPhFm5RUj6dmHnK9R/ZwL6bv0ex3qew+rls0S/3isRFW1MpgYSey8Wh/85jL+H/42QkAWCrdwiZm8W6lc8BdatOQoH23lYsvgrgaMjIioYk6kBKJVKfHX8K4QEhCBs5jazWLnFHBXWDdyqlepC/crM9tiw6wcsmC98oX4iooIwmRZAmwW2f779M648vIKdXXfivWHfmPXKLUIqrBs4IaHgQv2ylPJY/etqDG893IjRERHpjv2Ob1AoFAgJWYCAgLEANiEgYCxCQhZAoVCoPF6pVGJqzFSMbzIe2WnZolm5xdzkFupXxdPtOSaenoirj6+aOCoiIs0wmb4h971dXNwyABMQF7cMERHvYty4UJXH/3j7MG4m3cTIRiMLTQhmuXKLGckt1G9ndyLfdju7E+j5UQ182eRLdN/ZHanyVGECJCIqBLt5X6P1AtsSBb6NnYavmn0FB1sHABB05RaxV1cqrFC/AgrE/huLYYeGYdNHm4y6yLqljcImIv0xmb5GkwW2PTw8Xm2stg+P0x5iWINheZuEXLlFLEmzIIUV6pdCiq3dtqLeynpYdWEVhjYYatB755uWw1HYRKQlJtPXaLPAdrYiG2j1NcY1mgp7G/u87WJcucXcvCrUn3+7h4MHdnTfgQ5bOqCBZwP4e/rn7dN3CbU3p+VwFDYRaYPJ9DXaLLC998YOwDYN/WqpXjKsoIRA+mleoTmmt5iO7ru648KQC3C1z/mAo88Salp37xMRvYH9V29YuHAsgoPPwcdnBIB58PEZieDgc/m6aRVKBRac/hY4PhzIFi5WSzW+yXj4lfFDv339oFCqHmWtDU2694mICsNk+obcbtrY2DAA/REbm9PNl/veTKFQoPOAYYibUxH4Q6F26gwZnkQiwYYuG3D18VXM/XWu3tfjKGwi0heTaQEKWmB77NiFiNrWFdlPDgOYqHbqDBmHi50L9vTcg9m/zMbxW8f1ulZh03JMMQqbiMSPyVQLcrkcO/dfgzKzfb7tOe/W7v43rYJMpV7ZeljSfgl67+mNf2X/6nUtTbr3iYgKwgFIWnjy5AkeP3NXuU/l1BkyusH1BuO3e7+hx64eODngJHSdG8pR2ESkDyZTLVyRXUG2zR8q9/HdmvGpnv4iQZfAcFzOboaQ6BDMabZci3PfnjrDUdhEpAsmUy3MOzMPjVq54mKkMBWOLF3B01/s4fd8D/xX+cPPtRGAT7U4l4hIf6JLpqdOncL+/ftx48YNpKWl4dixY7CysjL6fU/Hn8aFBxdwZ/0dfDv1B0EqHFHBKrpUxJauW9B9Z3egjB+AOkKHREQWRHTJNCMjA/Xr14e/vz/WrFljsvvO/20+vnj3C7gWd+W7NTPV3qc9RjeYiNkfd0VS+jk4ObHbnYhMQ3TJNDAwEABw6dIlk93z2uNriI6LRkTHiLxtfLdmeIYo1D+u0RTMXncBg6J64cf+h2EtFd23OBGJEH/TaGDhbwvRp3YfeDp6qj1W7Cu3CMkQzygrMwvYOx/3G3bG5OOTMT9wvmGCIyIqBJOpGgmpD7Dljy24OPSiRsczaQoj36ovci+kL2uOpeV2wq+0H/rW6St0eERUxDGZqvH9xaUI8glCdffqQodChXhz1Zf424DNg+MYOHwgauyqgfpl6wscIREVZRaRTGfOnAlHR0cAQFBQEIIKaTq+3k3boKkMK85GoMuLw4guyRanuSpo1ZfMjDYoFd8UnTd3xrngc/BwYEENIksXHR2N6OhoADBo1TqLSKbTp09H+fLlNTr29W7ahb+tgu3ftbBjUFMjRkf6KmzVF5usOmjoKsdHOz5CTP8YAHamDY6IzMrrDSqZTIbw8HCDXFd0tXllMhni4uJw//59AEBcXBzi4uLw8uVLg95Hni3Hov8twoQmEwx6XTK8wld9uY91vddBAgkGHRgEpVJp4uiIyBKIrmX622+/Yd68eXn/HzZsGABg0aJFqFu3rsHus/3P7XCwdcAHvh8Y7JpkHOoWdXcq7oR9H+9DwzUNseDMtwCmvXUNfUdha1qukIiKJtEl0/bt26N9+/bqDyyAJr/0lEolFv1vEUY3Gg2pRHSN9yIv5z3HU8jlbsgtbJ9TgSq0wMpUZRzK4FDvQ2iyqglQtTTk8oF4vSi+vkmP5QqJLJvokqm+NPmld/LuScQnx6NfnX7CBEkq5Zv+Ai8EBNzLS5jqVn1RKBRYN+8oHPd2R+qTv1G3QTB6dq2Wdy4RkT4sLplqIiw2DEP9h6KEbQmhQ6HXvDn9JS4OiIg4ASAUYWE52wqqTPXmuffTgBUrfs53LhGRrviR/A3/PP0H0TejMfzd4Sr3R0cDISE579Ry362FhORsJ+N5Nf2lZb7tmizMXtC5GRmtcSDyJhd1JyK9sWX6hiVnlqBHjR4o51RO5X4OKBFGYdNf1C3MXti58Y8cEJ8QD58KPgaLlYgsD1umr3n28hnWXVqHMY3HCB0KvaHw6S+FL8xe2Lm2xeMw9OehyMjKMEicRGSZmExfs+rCKjTwbAB/T3+hQ6E35E5/sbM7kW+7JguzF3buwJ7vIjU7Fb329EKWIssIkRORJWA3738yszOx/NxyLOuwTOhQqADqpr/odu5XeJ4RjBbrW2DggYHY0GUDp0MRkdaYTP+z99pe2Eht8EFVFmkwV+qmv+h6rqu9K3769Cc0X9ccw6OGY0XHFZBIJMb8UoioiGEy/c/Ss0sxsuFIWEmthA6F1NBnYfaCzvVw8MCxfsfQbG0zjP9pPBYELhBNQmX1JSLhMZkCOP/gPC4nXkbUJ1FCh0IC8nb2xvF+x9F6Y2tkK7IRFhQmioTK6ktEwuPLIQDLzi5D/zr94WLnInQoJLAqblVwcsBJ7P17L0YdGaV1YfycOasJnLtKZGG0apnOmTMHxYsXh5+fH2rXrg13d3djxWUyj188wo4/d+Di0ItCh0Jm4p2S7+DkgJNovaE1MhWZWNFxhdpBSepKHRJR0aZVMvXz80NMTAyio6Mhl8vh4eGB2rVro3bt2vDz84OXl+qJ8eYot1j6D7+vRPMKzVHdvbrQIZEe9F315U0VXSrmJNSNrTHowCCs/mA1bKxsCjxek1KHRFR0SWJiYrRe4DE7Oxt///03/vjjD1y5cgUXL16EXC6Hu7s7unXrhm7dupnFp/G0tDR06tQJ9+7dy1sc/PUWRFxceViVOolOHStg79pws4iZ1Mt9N5icrP27QW3PTUhJQIctHeDp6ImdPXbCwdbhrWPkcjlq1x6LGzfenlZVteoI/PFHWKHzYA1Fn+dCZIlkMhmcnZ1x6NAhlCihXy12nbKHlZUVatasiV69emHOnDlYsWIFgoKC0LZtW+zfvx+TJ09Gdna2XoEZS24LIi5uGYCJyH5yGNE7emLcuFChQyMzVNaxLE4NPIWM7Ay02tAKj9IevXWMJqUOiaho0yqZymQy/Prrr3jy5Em+7ZUqVYKXlxc+++wzbNy4EX5+fti8ebNBAzWEgoult1JbLJ0sl1MxJxz+5DB8XH3QdG1TxCXF5duvSanD3AUSvvwSCAjI+fv1BRLU7Sci86bVO9NvvvkG9+/fx6NHj1CvXj20aNEC1apVQ2ZmJuLicn7BWFlZ4ZNPPkFERIRRAtaHPsXSybIVsy6GLV23YMJPE9BwdUPs6L4DgZUDAbwqVxgffyLfB7XXSx2qm75SlKe36DMPlnNoSSy0Sqa1a9fGggULcOvWLURHR2PTpk14/PgxrK2tMXZsTkm3//3vf0hOTkbJkiWNErA+clsQiYlv71NXLJ1IKpFiYbuFqF26Nrrs6IJvWn2DMY3HQCKR6FXqsKjT54NCUf6QQUWLVsm0WrVq2LRpEwIDAxEcHIzg4GCkpKTA2toa9vb2AICrV69i+/btGDp0qFEC1kduC+JufAwy0lvlbdekWDpRrv51+6O6e3V03dEVFxMvYlWnVbC3sde51CERiZ/Wo3llMhnOnz+P1q1bF3hMcnIynJ2d9Q5OXwWN5m3RuxfOn8xA+sMm+VoQHM1r3sytuzAxNRHdd3ZHckYytnXbhlqla6kdUavv/sKY+2jeovy1kTgZcjSv1uUEnZycCk2kAMwikRZECSXuBZxFxJh5GBjQgi0IEdEn8RnjHZuHgwdi+sfgm1PfoNGaRpjXdh4+9R0OQJgShLlzp+VyNwCvvqf53pHI+CyuNm/UP1HIVmbjw2rdAFjrVCydKJeNlQ1mtZqFtu+0RZ+9fRD1dzRQYg2AMiaLQV31Jb53JDI+i+vXXH52OYb5D4O11OI+R5ARNa/QHFeGXUFxmxLAiOrY+OcPWtf11VX+udMTEBe3DBER73LuNJEJWVQy/fvJ3zh19xQ+9/9c6FBIZDSZB1rSviTWddwO7NuAef+biVYbWuH6k+tGjavgudMtOXeayIQsqnm24twK9KjZA6VLlIZMJnQ0JCZadZXe+ABn+rXCgvPTUG9lPYxsOBJfvfeVUVYlMsTcab5TJdKfxSTTlIwUrL+0Hj9++qPQoZAFcLB1wKL2i9CvTj+M/2k83lnyDqa8NwWfVhsOwE6raxWW7Fq10n/uNN+pEunPYpLp5iubUdWtKhqVayR0KGRB6pWth58+/Qk/3vwRE49NxNIzy4AGk/Ayqz+cYK/RNQpPduqrLxGR8VnEO1OlUonwc+EY/u5wSCTCTFsgyyWRSBDkE4Tfh/6OyQ1nATWWoNbKCph9ajaevXyW71hdFhdfuHAsgoPPwcdnBIB58PEZieDgc6y+RGRCFtEyPXv/LB6kPECvWr2EDoUsVL7pK3EDYFPhOr4/eQBzms1Bb7/eGFR3EHYt+hWHDt2DtouLS6VSo1ZfMod3qgXNoSUyF6JNplu3bsXevXuRmpoKf39/jB07tsD3Q+svr8egeoNgb6NZtxqRob25eHjCXcDu4Ql08z4E23ov0OqTTyA/8z2QNQGAbouL53Tpehh87rSh3qnqkhDVzaElMhei/G48cuQINm3ahFGjRmH58uVIS0vDzJkzCzz+aNxRBDcINmGEZGlyp85MmwY0bpzzd+7UmcKmr5w7mY75LebD+1FHIKv9W/v3HriBjIwME34lhqdQKBASsgABAWMBbEJAwFiEhCyAQqFQey7n0JJYiLJlum/fPnTr1g3NmzcHAEyYMAF9+vRBXFwcfHx83jq+mVczVHatbOowyYIU1uWZkFD49JWbN28iNaWCyv13H5VAxe8qIrBOIAI8WgPOraBUekOokoWqqOsGfrNVrmmr+9WHkPzH5Myh3Y25c+UcYEVmQ3TJVC6X4+bNm/lWpfH09ISHhweuXr2qMpla/zUAISGvftnl/uDntiA4p46MSd3Sf5UrVy5wf5WyWVjWez1+uf8L1v2xEvjyM/iuckcjr4Zo6NkQDTwboHaZ2ijrUBbqEqyx3jsW1g2sT0Lk+sMkJqJLpjKZDAqF4q31Ul1cXPD8+XOV5ywa1RK1ar76P5MmmZK6xcMdHBwK3N+pUwUE+QYhyDcIExp8C+dSadhw+SL+en4WZ++fxbpL63Dr2S2UtC+J6q61gI41EfG7L/zKVUVVt6qo4FIBUkgFe++oT0Lk+sNkLLm9KS9eGO6aokumutQ7ffniBWQseUQGkPttpO2309dff46MjOU4enQHbt3yxjvv3EP79p74+usRkMlkavfn3TMTqOnoh4Byfvis5mcAgDR5Gq4/vY4L967i9MtrOH4jGmvPr8CtZ7eghBIlTvkg9cwSKDJfdbOGrziO1JczsXDeq+kz6r42Xfbb2NjAweGuyuNLlLgLa2vrQn8227Ytg/j4GKTnW384Bm3blkF6ejrS09MLPJeoIAEBQOPGSsw8PBdYaZhrar2eqdDkcjk6dOiA+fPnw9/fP29779690bt3b3z44Yd523LXMyWybIEAVFX+CgRwzAT39wCwFkCH17YdATAIgIpmp8rzawGoAeAagD80PI9IM4KsZyo0W1tbVK5cGZcuXcpLpgkJCUhMTESNGjVUnnP27FX4+pYzZZhUxCgUCkyZshxHjiTg9m0vVKp0Dx06lMXs2SO06iqVyQAvL+DevYIXBy9ovy7nJiYmonnzvXj48O3j3Uu3RPiuIUiySkJcUhyuPYxD1Jm/YVX6H9hZ26FaqWqo61EX/mX94e/pjzLWVVCxgpXWsec+u6NHD+Zrdc+efU3jZ/fkiRyVKz/DzZsluf4wGcSQyCGwSrfH1s/WG+R6okumANClSxcsX74cVatWRdmyZbFixQr4+fmpHHwEAA4OjnBisVHSQ0jIAqxd2yzvnebt28DatSdQrNhqjeeBvs7JqfC5mqr25w4gsrNzg5NTwQnl9XPt7Ozg7HxfZTIt6ZKIzo3H5w0AevJEDvfPnyI+0RFPcRt/PPoDvyf8jm3/bMO4U+NgJbEC+jbBuhtt8H611qjrURdWUiuNYg8Pn5pzffcknDmja1GJUihVinWDSf9CIo/THmP/7f048uEv2Ir1BolJlMn0/fffx7Nnz7B48eK8og3jxo0TOiwqooSeoqFP4QJ1g59sbW3fun6LZq+u/0ntTwAAWYosnLn9F5r1OYn/3T+OhWe/hZXECh2qdEDXal3R3qc9gMK7yYxVVIIsj76FRNZeXIum3k1RxdXXYDGJMpkCQJ8+fdCnTx+hwyALIPQUDV3naebKqdEbisjI3YiL84KPz795yVLT61tLrVHbvQ5wpg62/TgKxR2ycOHBBRy4fgBTfp6Cvvv6ok2FIKBGX8izPwRL/pGQCmu5tg3MRsT5CIQFhRn0nqJNpkSmIuQUDUO0igur3avN9V+fp+oktUWj8o3QqHwjzGkzB9ceX8Pmi7sQFTgB1Vd/gX51PsXg+oNRw131OAYiYyqs5Rp5/TCyFFn40PdDJD4w3NwYUZYTJDKl3K5SO7sT+babYpkzTVrFhZUyfN2rbtZX8WpyfU3KAVZ3r46Jjb8GlsZh7fvbkZiWCP9V/gjaHITjt47rNKWNyBjCz4VjqP9QWEsN25Zky5RIA+q6So1Fk1axPtW7NLm+Vt3MSilaeLfGB7Va43HaY4SfC0fP3T1RyaUSRtSbAEi6g5/hSSj/PP0HMXdisL7LeoNfm9/VRBrI7SqNjQ0D0B+xsTmJxNgVhIzdKlZ3fQAFFumPirqbb93VN9didS/hjhktZyB+dDz61+mPab+MA4b44+e7P+kVM5GuIs5HoGv1rvBwMPwYB4tomcrlmUKHQEWEECNSjd0qLuz6Dx8+VNsNXLp06UJHG5ewLYGRjUbi4yqfo0yncAyM+hgrLjXAvLbzUK9sPYN8DUTqpMnTsO7SOkT2jjTK9S0imXbuPAddu1blGogkSsZe/Luw6xuyG9jO2g6IHYtL+wch/PJ3aLq2KQbXG4zZbWbDqRgnj5Jxbf1jKyo4V0BTr6Z52zIz5YWcoR2LyCx374ZyDUQSPVUDiIx9fUN2A+cqaVcS8wPn4/Kwy/jr8V+oEV4D+//eb+CvhugVpVKJ8HPhGNFwBCQSSd6gurZtpxjsHhbRMgW4BiIJ6/V5b6qW/lO3X0j6dgMXNAe3ilsVHO93HBsvb8RnBz/DxssbsbLTSriXcDfiV0OW6H8PTuNu8t28IiSvelOGAlhlkHtYTDIFuAYiCUddUjSHpFkQfbuBCyORSNC/bn+8X+V9BEcFw+97P2zsshGBlQMBmPeHDBKPVZeWY1DdQShuU/yNudWGW03MopIp10Ak0p2qwVealCvUJCG6l3DHrh67sPbiWnTd2RXD/IdhdpvZCAqyZdKkAmm04L3jAxy6uQ9/Bf0FoPC51fqwmGRqign2RJZI3WhjTVuREokEg+sPRjPvZui9pzdifojB3o/3wtvZ27hfAImOVvWq/VehpXdb+LjmLIRSWG+KPiwimVasGIKPPvI1+gR7Iktk6NHGvqV8ETs4FqOPjkaDVQ2wu+duNK/Q3IARk9hpOoI89WUq4LsCA2usyduWvzelvsFisohkun//FNSpU0noMIgEYar3jrrMwS24IHkxRHSKQL2y9dBhSwcsCFyA4AbBkEgkhgu40Pvznay50qSetLW1NcaNC8XWPX8CT/ph6vEfcfaD63kt19zelAMHtuLWLcPEZRHJ1NbWRugQiARjzolB3VJaQ/yHoIZ7DXTf2R2XEy8jvGO4QWuq6ruUF5meJvWk58/flK/levONlmtub8rIkU/wzjurDRKXRcwzJSLxaubdDOeHnEfsv7Hosr0L0uRpQodEAsp956mKo+O/cHBw0Hjus42NAYufGOxKREWYpiuzkHGUdyqPUwNPIS0zDa035hTRJ8ukrpBISkqK2parMVhENy+Rvsy5q9RSuNi54Gifo+i3vx+arm2K6L7RqFSSYyEsUWEjyLOyslDc4bbK84w5PZLJlMjMGXsAkZgKIxSzLoZt3bZhzNExaLq2KX7u/zOqlaomdFhkYoWNILe1tYVHvce4czcaisxX38DGnh7JZEpk5oyd1MwxaRZGKpFicfvFcLB1QJNVLdEp6ThKZtbkaFwLpGoEeWZ2Jm43/A2dipfC1dORJlt/mMmUiERHIpHg29bfwtbKFuHnWmF/12NY6uPH0biEvdf2wt7WHnvXhuNZUrZRVlpShcmUiAplrt3AEokE01tOh7XUGp12twI8jgHg+qiWbunZpRjx7ghYSa1ga2sFU60/zGRKRIUSOmmqM6X5FGRlWmNGv0Bce3ISjZxqCh0SCeT8g/O4nHgZUZ9EmfzenBpDRKI35t2JwNkR6LI3EDeTbgodDglk2dll6F+nP1zsXEx+byZTIioaTkxHN99eaLupLf6V/St0NGRiCSkJ2PHnDoxsNFKQ+7Obl4j0Yoh3qhotpaWWBLObhyJDmYLATYE4OeAkSpcordf9WbtXPFacW4E277QRbKqURSTTBQsAFxf+ABAZgz4/V1otpaUBiUSC7zt9j777+qLDlg44OeAkHGwddL4/a/eKw8usl/j+wvfY1m2bYDFYRDKdNQsoX17oKIjoTZoupaUNK6kVNnTZgPe3vI8eu3rgYK+DsLFSvdiFMe5Pprfj2maUdSiLNpXaaHR8bo/DixeGi4HvTIlIEK+W0mqZb7uqguTasrWyxZ6ee/Ag5QGGHhoKpVJp0vuTKSkR8ftijG48WuMl+oKCgLAwYP58w0XBZEpEgtBkKS19ONs54/Anh3Hs1jHMPDnT5PcnE6n8I56mP8EntT8RNAwmUyIShLqltAxRkLycUzkc6XMES84swQ+//2Dy+5MJBIThM78vYGdtJ2gYokqmly9fxqRJk9ClSxe0atUK9+/fFzokItKRuqW0tClIntMlm6Cya7Zm6ZrY23MvRh0dhRN3Xt3LkPcn08pdEvGLGX9BUukUHh8JFnxJRFENQEpPT0fVqlXRrFkzhIaGCh0OEempsKW0NKFuNO6rqS2t4PFgCTpkdEPvtDP4ONAHQUH635+EkTvK+vODizFQ+QlWd84/BUqIEpiSmJiYt9/Mm7nExET07t0bmzdvRrly5Qo8Li0tDZ06dcK9e/dQnsN5iczWkyc5S2k9fqxdQfKQkAX/jcZtmbfNzu4EgoPP5RuNmzu1ZcTBsfjpbhRiB8eipH1Jje+fe35yMqfGmIvE1ERUWlIJ5z8/j5qldSshKZPJ4OzsjEOHDqFEiRJ6xSOqbl4iKppeLaWlXdeutqNxZ703H1XcqqDn7p7IzM7U6/4krKVnlqLtO211TqSGxmRKRKKky2hcK6kVtnbdioepDzEmeoyxQyQjSclIwYpzKzChyQShQ8ljFu9Mw8LCEBkZWeD+OnXqYPHixTpff+bMmXB0dAQABAUFIYhlkIhEL3c0bmLi2/sKG43rWMwRB3sfRINVDVC/bH0MqjfIyJGSoa3+fTWqlaqGZt7NtD43Ojoa0f+NVDLkXGKzSKZDhgzBp59+WuB+GxvV1Us0NX36dL4zJSpickfjxsefeOudqbrRuBVdKmJH9x34cPuHqOFeAzWcGhd6L8PUDiZDyMzOxKL/LcKS9ks0LtLwutcbVDKZDOHh4QaJyyySqYODAxwcCq6fSUSkij6jcdu80wazW89G1x1dcaL3BQBl3zrG0LWDSX/b/9wOe2t7dPbtLHQo+ZhFMtXUy5cvcf/+fTx58gQAcPfuXbx8+RKlS5eGE4fYEVkcqVSKsLDxmDw5ZzRubKx2o4G/bPQlfk/4HZ8e6gZYxQAolm8/a/cKo6DVetq1U2L+7fkY12QcrKRWQoeZj6iS6fXr1zFmzKtBA1OmTAEATJw4Ee3btxcqLCIS2KvRuNqdJ5FIsLLTSjT9oTnQYRSAlXn7Xo0Wzp80c0YL78bcuXKO/jWSglbrOfLPUTz68xH61ekndIhvEVUyrVu3LmJiYoQOg4iKEHsbe2zqtAe1btfH5r8a44uAgQA0Gy3s4eFhylAt3tzTczGq4SjBSweqIqpkSkRkDF5O3sDu7Rjv1BmNK9ZB/bL1dR4tTOrpsuj6qbuncDnxMg72OmjaYDXEZEpEBAC32mJco6notrMbLgy5AFd7V51HC1PhdFl0ffYvszGy4Ug42zmbJkgtMZkSUZGn6dSWMe9OxOUnZ9Bnbx8c6n2ItXvNxIXEszgdfxpbum4ROpQCMZkSUZGleSH8nILo07+WopzVBhwo8y6+OfUNZrScoddoYTKMhWdmI7hBMEoVLyV0KAViMiUiwRh7dQ91U1tU38cZwx7uQcAPAWhcvjHa+7TXebSwkHR5L2mWylxGTPxP+OGjleqPFRCTKREJxpi/2PWZ2lK7TG2s6LgCfff2xe9Df4eLxNs4QRqRLu8lzdJ7c/BprcHwcDDvkdMs4UFERZIuhfBf169OP3St3hU9d/WEPNtwNVxJczeS/gaqHcAof/MvkMFkSkRFUu7UFlU0ndqytMNSZCoyMeUUBxwJYcGZb4ErfXOmLpk5JlMiKpJyC+Hb2Z3It12bqS121nbY3WM3dlzbDNTabqRISZW/Hv2FA//sBk5NFToUjfCdKREVWYaY2lKpZCWsDNqIXrI+uJFUFw2cqhkvYMoz/cR0fOI7AOueF/tvapN5j/5iMiUi0VI3GljfQvi5OlT+ADj3Bfp5dse5IWdQwraEEb4aynX+/nnsX3IGXomdAGwSxWo9TKZEJFqajgY2yNSWn7+F68exCI4KxoYuG3RaS5M00/PzUcC5H3Ansx0AcazWY54pnojI3Cissfb97fjx5o9Y8/saoaMpEnK6bxP++zvHqZuncOesE7L/S6S5cqY03c13rDlhMiUi0pCHQ1ls67YNo6NH42LCRaHDES2FQoGQkAUICBiLnG7csQgJWQCFQoEpUVNQXNlA5XmaTGkSCrt5iYi00KpSK0x5bwq67+qOC0MuwMXOReiQRKegylT3kr/ElUpX4OFaGzdV5ExzXq2HyZSIqAAFDXAKbDcJ1Uv9hg6r+qPxnX1QZEtNXrJPrOUCC6tMFRkVislbx+J5kj0iIsS1Wg+TKRFRAQpOTFI0erkR/qv8UbbbQgyrPcHkJfs0LReo6Yo5plJYZarMF1XQt0pfVFxYEWJbrYfJlIhIB672rtjVYxdarG+BWiUbAWiRb7/QLUd1K+YIpbBF10u7PkH5MuUNNqXJlJhMiYh01MCzAcLahWFg1MeA4+8APPP2CV1oXt2KOULJrUz15qLrEpuj6PVRrXzduGJarYejeYmI9DDEfwjaVAgCephPQfxX7yVb5ttuLtNLFi4ci+Dgc/DxGQHgO0hLtUe3Tw8jNHScoHHpg8mUiEgPEokEYW0iANtUTD1lHslA3xVzjC23Gzc2NgxocwPt5ttg1w9Lzba6kSbEGzkRkZkoblMc2LEX269twuYrm4UOxyAr5pjCndRbQOPt+KbVAqFD0RvfmRJRkaWudq9BPXsHazpswYBDPeFXxg9+ZfwMfAPNFfRe0pymlyiVSoyPGQn8/jmquop/8QAmUyIqskw957Jdpfcxvsl4dN3RFWc/PwtXe+FagIZYMceYdv61E1ef/AH8vFvoUAyCyZSIyICmtZiGCwkX0Gt3LxzucxhC/Zo15+klsgwZxkSPwZwWYfgsw1nocAyC70yJiAxIKpFic9fN+Ff2LyYdmyR0OK9NLzGPRAoA036ehuru1dHdt7fQoRgMW6ZERAbmVMwJB3odQMM1DeHrXBdAX6FDMhsXEy5i9e+r8fvQ34vUMnZMpkREelJVsq+KWxVs77Yd3XZ2AzyrAVC9Eoox7i20gqo/BbZTYOa/wQgJCEG1UtUgkwkdqeEwmRIR6Uhdyb4gnyBMajwD0552QULqWTg5eaq9pqHuLaSCqj+Fn43Aw+sPMfm9yYLGZwyiSqaHDh3C0aNHcefOHVhbW8PPzw/Dhg2Dp6fhvkGJiDSlScm+obVHYtr8c+ixpyNih/6KErYlTHZvcxKXFIeJxybiYO+DOfNyC2HSKU0GIqpkevnyZbRr1w41a9YEAKxevRqTJk3C2rVrYW0tqi+FiESusKXEoqJ2Y86cdEyevCyn5RhXH/9cyoDfkSD8fTAGNtY2Rr333LlysxpwlK3IRv/9/TGo3iC0rtRa7fHmnDQLIqoMNGXKlHz/Hz9+PLp37467d++icuXKAkVFRJZIXcm+L7+ci40bW+YlvPSHwK2fjqBxz49wYe8ho947KSkJHh4eet3DkJZdCMWTF08wt+1coUMxGlEl0zclJycDAJxMuRQDEREKX0rMweEuYmKkbxWaR1YHXPplCZb9tgwjm4w0yr01LReozxJxWp1b+g/M/d8M/Nz/57zuXTF246oj2mSqVCrxww8/4N1334W7u7vQ4RCRhSmsZF/Llm6IjHRUeV5JaRNMjJwIL1cvdKnWxeD31rRcoD5LxGm8MHm2HPioH4LrjUbj8o3fOr8oMYtkGhYWhsjIyAL316lTB4sXL863bcWKFbh9+zaWLVtm5OiIiFQrqGTfnDlTcPLkeJUtRzeXx4jovQF99/bFoU8OoWXFlga9t7mUCwSA6b9MBCRKTGo8XehQjM4skumQIUPw6aefFrjfxib/y/rVq1fjxIkTWLp0Kdzc3NRef+bMmXB0zPmUGBQUhKCi9pGIiARRWMm+wlqOPfx6ICU7BZ23d8aJ/idQr2w9g97bHOz4cwe2XF0H7DyPYvOKCR1OnujoaERHRwOAQdd1lcTExCgNdjUT2LBhA/bv34/FixejQoUKhR6blpaGTp064d69eyhfvryJIiQiS5Pb3Zmc/Kq789U80LtvtRxz54HOPz0fobGhOD3oNHxcfQx2b0Dz95oFna/Pvf969Bca/9AYa9pvRa/6H+h0bVOQyWRwdnbGoUOHUKKEflOWzKJlqqmtW7di27ZtmDVrFhwdHfMWuHV0dHyr9UpEJCRNWo4Tmk7A47THaLuxLU4MOIGKLhUNdn+h3kvKMmTourMrRjcajQ6VPzB9AAIRVTI9ePAgMjIyMHHixHzbFy1ahLp16woTFBFRIV4Vmle9f37gfKRnpaPl+pYGT6implQqMfDAQFR0qYgZLWcgLVXoiExHVMl0+/btQodARGRQEokESzsshVQiRYv1LXCi/wlUKllJ6LB0mjoz5ecpuPDgAs4POQ8rqZVZ1g02FlElUyKioiJ/spLg3YaL4VZCgkbft8D/hp3AP2ff0XkeqCFoO3VmUewirP59NX4d+Ctc7VwRErLALOsGGwuTKRGRAN5OVhI4Oi7C2B+leG/de4j6JAphQXV1mgdqatuvbsL0E9Pxc/+f4VvKFyEhC0RVN9gQiuZHBCIiEZJIJAhtF4oR745Ai/Ut8NPNn4QOSb0qURhzfBj2frwXDTwbvFY3uGW+w3LqBt816HQUc8JkSkRkRiQSCb567yuEvx+OLju6YMtf64UOqUDH7hwFuvfCiqD1aPtOWwCa1Q0uitjNS0Rkhvr69YWnoye67ugKtLqJbMUMAFYGv4+ug4TWX1qPLw58AewNRceQznnbDVE3WIzYMiUiMlOtK7VGdM/TQM2d6LK3HRJTVWQoHSkUCoSELEBAwFgAmxAQMBYhIQugUCgKPU+pVOLbk99iyIjpKLmuN3A9Jd+5uXWD7exO5DtPm7rBYsSWKRGRGateqiaw6jzKbByGut/XxdZuWzVaE1QdXRYXz8zOxOijo7F+XjQk59fiQUYbleeKoW6woTGZEhHpyGRLickdsbr9Zuy6+QM+3PYhRjQcgQD51zh5rLhO99ZlcfEbSX9j2Pa+SM9IR5kHbXH7v0Ra0LnmXDfYGJhMiYh0ZMqSfRKJBJ/V/wyNyjXCkENDsC2lOpYMWYLOvp0hkUi0upY2i4srlUqgYThabp2E4AbB+KL6F2gyZ5dG56qr/lSUMJkSERmBsVqttcvUxulBp7H+0np8Hvk5Vl1YhdB2oajuXl3ja2g6SOhS4iWEHJkANLmOHZ0PoWPNlpDL5RY5wEgdDkAiIjKCoCAgLAxYsgSIjc35OyzMMC1ZqUSKQfUG4fqI66jkUgn1VtZDp62d8PPtn3NakmqoGyT0x5M/0Hl7ZzT5oQmqudUAIq7gPa+WGp1bVAcYqcOWKRGRSLnauyK8YzimNp+K8HPh6LmrJ8o5lcOAOgPQrnI71HCvUWAX8JuDhCq+cwc+AWn4q+5DNF8/A8P8h2Hllythne6KiIynkMvtkTt1xhIHGKkjuvVMtcH1TInI3KlbT/TJEznc3Z/i8WM3tYN4XmS+wJYrW7Dv7304efckXOxcEPhOIGqVroUyJcrAw8EDpUuURlpmGh6kPEBCSgKuP7iD8L0/QVrpLzSu2BidqnTCoHqD4F7c/b/1WOP/S5hv19fNiS0Jjx+rHmCkz1qppmCx65kSERU1BRVNeLW4uObF4ovbFMfn/p/jc//PkZGVgd/u/Yafbv2Es/fP4mHaQySmJuJR2iOUsCmBso5lUdahLEoV8wQuT8TNhe1RsYxb3rU0qa9rSQOM1GEyJSISgLpkqcs80NcVsy6GVpVaoVWlVoUeJ5MB67oBrvavtukydcbScQASEZEAcpNlXNwyABMQF7cMERHvYty4UMGLxVtqfV19MJkSEZmYumSZkJAgaDLLnTqjiiVPfykMkykRkYmpa/lJJBJBkxmnv2iP70yJiExMXdEEDw8PdOzojfj4E/lar6ZMZpz+oh0mUyIiE8tt+RWWLIVOZlKpVOf6uiarWWxGmEyJiASgLlnqk8wMSZfpL0U5aRaEyZSISACaJkvO5RQHJlMiIgExWRYNHM1LRESkJ7ZMiYgskD6DhCxxgJE6TKZERBZIn8RnyUmzIOzmJSIi0hOTKRERkZ6YTImIiPQkqnemR48exa5du5CQkACpVIqqVatiyJAhqFatmtChERGRBRNVMnV1dcVnn30Gb29vZGVlYc+ePZgwYQK2bdum9yrpREREuhJVN2/Dhg0REBCAcuXKoUKFChg2bBhSUlJw9+5doUMrcqKjo4UOQbT47HTD56Y7Yzy76GggJCRn2kvu9JeQkJzt9DZRtUxfl5WVhUOHDsHJyQne3t5Ch1PkREdHI4hj33XCZ6cbPrf8tJnLaYxnx+kv2hFdMr116xaGDx8OuVyOkiVLYv78+XBwcBA6LCIirahLlkxm4mIWyTQsLAyRkZEF7q9Tpw4WL14MAPDy8sKaNWsgk8kQFRWFWbNmISIiAk5OTm+dp1QqAQApKSmQyWRGib2oksvlfGY64rPTjaU9t4CAnD+qaPsYLO3ZGUruM8vNFfqQxMTE6H8VPaWmpuLly5cF7rexsYGLi4vKfZ9++im6dOmCbt26vbXv8ePH6Nmzp6HCJCKiImjnzp1wd3fX6xpm0TJ1cHDQuatWoVDAyspK5T43Nzfs3LkT9vb2kEgk+oRIRERFjFKpxMuXL+Hm5qb3tcwimWpq48aN8PPzg4eHB1JTU3Hw4EEkJyejcePGKo+XSqV6f9ogIqKiy1BjbkSVTFNSUjBv3jw8ffoUDg4O8PX1RWhoKDw8PIQOjYiILJhZvDMlIiISM1EVbSAiIjJHourmVWXr1q3Yu3cvUlNT4e/vj7Fjx8LV1VXlsS9fvsTSpUtx6tQpWFtbo127dhg2bFiBA5iKMk2fm0wmw9q1a3Hu3Dk8efIEpUqVQlBQEPr06WORzw3Q7nsuV1paGgYPHoyHDx/i2LFjfHYaPrvjx49j69atuHfvHpycnNC9e3f06tXLhBGbB22e2+3bt7FixQpcu3YNVlZWqFOnDoYPH44yZcqYOGphnTp1Cvv378eNGzeQlpam9udO3/wg6pbpkSNHsGnTJowaNQrLly9HWloaZs6cWeDxixcvxtWrV7FgwQJMnz4dMTEx2LBhgwkjNg/aPLenT5/i+fPnGDlyJNauXYvhw4dj37592Lx5s4mjNg/afs/lWrp0qcVX6tL22f34449YunQpevTogXXr1mHOnDnw9fU1YcTmQdvnNnXqVDg4OGDFihUIDQ1Famoqvv32WxNGbB4yMjJQv3599O7dW6Pj9c0Pok6m+/btQ7du3dC8eXP4+PhgwoQJuHLlCuLi4t46NiUlBceOHcPIkSNRo0YN1K9fH4MGDcKBAweQnZ0tQPTC0ea5VapUCTNmzEDjxo1Rrlw5NGnSBN27d8fp06cFiFx42jy7XL/88gvi4+Px8ccfmzBS86PNs8vKysL333+P4OBgtG/fHuXKlUPVqlVRr149ASIXljbP7fnz53jw4AH69OkDb29v+Pj4oHv37rhx44YAkQsrMDAQffv2Rc2aNdUea4j8INpkKpfLcfPmzXw/XJ6envDw8MDVq1ffOj73m6lu3bp52+rXrw+ZTIb79+8bPV5zoe1zUyU5ORmOjo7GCtFs6fLskpKSsHz5ckyaNMliu3YB3X5enz17huzsbAwcOBA9e/bEd999h+TkZFOGLThtn5uTkxPKly+PH3/8EXK5HC9fvsTx48fRoEEDU4YtOobID6JNpjKZDAqFAiVLlsy33cXFBc+fP3/r+GfPnsHBwQHW1tb5jgWg8viiStvn9qYHDx7g8OHD6Nixo5EiNF+6PLvQ0FB07doVFSpUMEGE5kvbZ5eYmAgg513h0KFD8fXXXyM+Pt7iuiu1fW5SqRQLFizA+fPn0aFDB3Ts2BEPHjzAV199ZaKIxckQ+UG0yVTbWoqqjrfEqkj61KB89uwZJk2ahNatW6N169YGjEoctH12R44cQXJyMnr06GGkiMRD22enUCgA5JQLbdy4MWrVqoWxY8fi/PnzePTokTFCNEu6PLfFixejQoUKCA8Px5IlS1C8eHGL+xCiLUPkB9GO5nV2doZUKsWzZ8/ybX/+/LnKOr6urq5ITU1FVlZW3qeP3HMLqvtbFGn73HIlJydj3Lhx8PX1xejRo40bpJnS9tldvnwZ165dQ2BgYL7t7dq1w+jRo/HBBx8YM1yzou2zy22JvT5oK/ffjx49QunSpY0XrBnR9rldvHgRFy9eRGRkJGxtbQEAX331FXr06IFbt27hnXfeMUXYomOI/CDalqmtrS0qV66MS5cu5W1LSEhAYmIiatSo8dbxVapUAZDzCy7XxYsX4eTkhHLlyhk9XnOh7XMDcl7Ojx8/HmXLlsWkSZMglYr220Yv2j67wYMHY82aNXl/xo0bBwBYuXIlWrZsaaKozYO2z87X1xfW1tb53lfl/tuSpnho+9zS09MhkUjy/Yzm/ju3tU9vM0R+EPVvxS5dumDPnj345ZdfEBcXhwULFsDPzw8+Pj54/Pgx+vXrh2vXrgHIeTHfpk0bLFu2DNeuXcPFixexdu1adO7c2eIGhmjz3NLS0jBhwgRYWVlh5MiRSE5ORlJSkkW9Z36dNs/O3d0dlSpVyvtTtmxZADkjpC1xAJc2z87BwQFBQUFYt24drly5gps3b2Lx4sVo1KiRxdXb1ua51axZEzY2NggNDUV8fDxu3ryJhQsXwtPT0+Le28tkMsTFxeV9CIuLi0NcXBxevnxplPwg2m5eAHj//ffx7NkzLF68OG8yc+6n/+zsbNy7dw8ZGRl5x48ZMwZLlizBuHHjYGVlhXbt2qF///5ChS8YbZ7bP//8g7///hsA8k2WL1OmDLZv32764AWm7fccvaLtsxs5ciTCw8MxZcoUWFlZoWHDhhgxYoRQ4QtGm+fm4uKCuXPnYtWqVfjiiy9gZWWFGjVq4LvvvoONjY2QX4bJ/fbbb5g3b17e/4cNGwYAWLRoETw8PAyeH1ibl4iISE+i7uYlIiIyB0ymREREemIyJSIi0hOTKRERkZ6YTImIiPTEZEpERKQnJlMiIiI9MZkSERHpicmUiIhIT0ymREREemIyJSIi0hOTKRERkZ6YTImIiPQk6iXYiEi9gwcPQiaT4c6dO2jXrh0eP36MpKQkxMXFYcSIERa3PiiRMbBlSlSEHTlyBJUrV0bfvn3RuXNnzJo1CzY2NvD19cWpU6dw69YtoUMkKhLYMiUqwmQyGWrWrAkAePr0KSQSCdq0aYPMzEwsWrQIdevWFTZAoiKCyZSoCPv444/z/n3lyhX4+fnBysoKVlZWTKREBsRuXiILcfHiRSZQIiNhMiUqorKzs3H+/HlkZ2fj2bNnuHPnDvz8/PL2b926VcDoiIoWJlOiIioqKgrjx4/Hw4cPERMTg2LFisHZ2RkA8Ouvv6JixYrCBkhUhEhiYmKUQgdBRIZ3+/ZtbN68Gd7e3qhcuTJevHiBCxcuwNPTEx4eHggKChI6RKIig8mUiIhIT+zmJSIi0hOTKRERkZ6YTImIiPTEZEpERKQnJlMiIiI9MZkSERHpicmUiIhIT0ymREREemIyJSIi0hOTKRERkZ6YTImIiPT0f+FYINjJADzZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.errorbar(x,y,u,fmt='o')\n", "xc = np.linspace(0,1,101)\n", "yc = f(xc,*popt)\n", "plt.plot(xc,yc)\n", "plt.xlabel('$x$')\n", "plt.ylabel('$y$')\n", "plt.axhline(0, color='k')\n", "plt.axvline(0, color='k');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### (c) Plot your residuals. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Residuals:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "r = y - f(x,*popt)\n", "r_norm = r/u" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFbCAYAAACOHWQYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAxOAAAMTgF/d4wjAAA6wklEQVR4nO3de1hU9boH8O8wgIDjgBiGopZKeL/ig1rGtjwgXkgTtUjTcm9NU1ORVPTs9q721sx7Ql7qsaO27WZaeSESw2Oiu5M3rCNqg5akkqjAAHJRhvMHh0lkBmbWWjOz1prv53l8ZGbWzLy/mVnrXb/r0mRkZFSDiIiIBPNwdQBERERKx2RKREQkEpMpERGRSEymREREIjGZEhERicRkSkREJBKTKRERkUhMpkRERCJ5ujoAKRw+fBhffPEFLly4gNLSUqSnp0Or1Vrdfu7cucjKyqpz38yZMzF27FhHh0pERCqkimRaUVGBvn37Ijw8HO+//75Nzxk7dizi4+PNt/38/BwVHhERqZwqkmlUVBQA4PTp0zY/x8fHB4GBgQ6KiIiI3IkqkqkQe/fuxVdffYWgoCBER0cjLi6uwaZhIiIia9wymUZFRaFVq1YICAjA2bNnsXnzZpSUlGDKlCmuDo2IiBTILZPpiBEjzH936NABHh4eSE5OxosvvgiNRlNnW5PJhJs3b8LX17feY0REpH7V1dUoKytDixYt4OFheRKMWybT+4WFhaGsrAxFRUUICAio89jNmzcxfvx41wRGRESy8emnnyIoKMjiY0ymAHJycuDj4wN/f/96j/n6+gIAcnNzodfrnR0aAGDx4sVYunSpS95bLtz9M2D53bv8AD8DV5bfaDSibdu25nxgiSqSqdFoxPXr13HlyhUAgMFggFarRUhICEpKSjB//nwkJSWhS5cuuHLlCjIyMhAREYFmzZohOzsbGzduxOjRoy0249bep9frXZZMvb29XfbecuHunwHL797lB/gZyKH8DXX1qSKZHj16FMuXLzffnj59OgBgzZo1CA4ORm5uLioqKgAAXl5e+OGHH/DJJ5+gsrISwcHBGD9+PMaNG+eS2ImISPlUkUxjYmIQExNj9fGMjAzz3y1btsS6deucEZZkhg4d6uoQXM7dPwOW373LD/AzkHv5NRkZGdWuDkLOSktLMXLkSBQVFbm8iYGIiJzPaDTC398fe/fuRdOmTS1uw4XuiYiIRGIyJSIiEonJlIiISCQmUyIiIpGYTImIiERiMiUiIhKJyZSIiEgkJlMiIiKRmEyJiIhEYjIlIiISicmUiIhIJCZTIiIikZhMiYiIRGIyJSIiEonJlIiISCQmUyIiIpGYTImIiERiMiUiIhKJyZSIiEgkT1cHQOqRlpaGtLQ0VFVV4X/+538QEREBrVaLoUOHYujQoa4Oj4jIYZhMSTK1SdNoNMLf3x9paWnQ6/WuDouIyOGYTImIFIatQPLDZEpEpDBsBZIfDkAiIiISicmUiIhIJCZTIiIikZhMiYiIRGIyJSIiEonJlIiISCQmUyIiIpGYTImIiERiMiUiIhKJyZSIiEgkJlMiIiKRuDYvkRNwYXIidVNFMj18+DC++OILXLhwAaWlpUhPT4dWq7W6fVlZGd555x0cPnwYnp6eiI6OxvTp0xt8DpEYXJicSN1U0cxbUVGBvn37Ij4+3qbt165di7Nnz2LFihX429/+hoyMDGzdutXBURIRkVqpomYaFRUFADh9+nSj2xYXFyM9PR3Lly9H165dAQBTpkzBpk2bMHnyZNZOiYjIbqqomdrjwoULAIDevXub7+vbty+MRiOuXLnioqiIiEjJ3C6ZFhQUQKfTwdPzj0p5QEAAAKCwsNA1QRERkaK5XTKtrq6ud59Go3FBJEREpBaq6DO1R2BgIEpKSnD37l1z7bSgoADAHzVUSxYvXgxvb28A4HQGIiKVq53OBgCVlZWNbu92yfSRRx4BAGRlZSE8PBwAcOrUKej1eoSEhFh93tKlSzmVgYjITdxbaTIajUhJSWlwe1U08xqNRhgMBvMAIoPBAIPBgLKyMuTn52PSpEnIzs4GAOj1egwZMgTr169HdnY2Tp06hS1btmDUqFEcyUtERIKoomZ69OhRLF++3Hx7+vTpAIA1a9YgODgYubm5qKioMD8+b948rFu3DomJidBqtYiOjsbkyZOdHjcREamDKpJpTEwMYmJirD6ekZFR57avry8WLVqERYsWOTo0IiJyA6pIpkTkXrjWMckNkykRKQ7XOia5UcUAJCIiIldiMiUiIhKJyZSIiEgkwX2m165dw6VLl1BQUACNRoOAgAC0b98erVq1kjI+IiIi2bMrmf7222/Yu3cvDh06hPz8fAB/rHVbu77tAw88gMGDB2PEiBFo166dxOFSQzjCkYjINWxKplevXsV7772Hw4cPw9vbGz179sSwYcPQunVr+Pv7o7q6GkajEVevXsXZs2fx1VdfYefOnYiMjMS0adNYW3USuYxwrF3H0pb1LInIOp4gK4dNyfSFF15Ahw4dsGjRIgwaNAi+vr4Nbl9WVobvvvsOu3fvxuTJk/HNN99IEizJm8lkQmJiIvbs2QMAGDhwIGJjY7Fy5Up4eLB7nshejZ0gK+nEVe0nBjYl09dffx0DBw60+UV9fX0RHR2N6OhoHDt2THBwpCyJiYnYsGEDysvLAdSskbxhwwYAwOrVq10ZGpGqKPHEVS4tZ45i06duTyKV8rmkHJWVldi3b585kdYqLy/Hvn37FHHmTKQUtSeuBoMBwB8nromJiS6OzH2JPoW5e/cucnNzcfbsWRgMBty4cQNVVVVSxEYKcvPmTRiNRouPFRcX49atW06OiEideOIqT4KmxuTn52P//v3IzMzEpUuXoNFooNfrUVVVhZKSEmg0GoSFhWHQoEGIjY1F06ZNpY6bZKZFixbQ6/XIy8ur91izZs0QGBjogqiI1MeWE9fg4GAnR0V2JdPKykps3rwZP/30Ex577DHMnDkTYWFh9QYkGY1GZGdn48yZM5g7dy6ioqIQFxfH64WqmLe3N0aMGIHLly/XOWP28fHBiBEj4O3t7cLoiNSDJ67yZHMyvXHjBtasWYORI0di1qxZDW6r1+vRv39/9O/fH1OmTMGBAwfw5ptvYsGCBfDz8xMdNMnTypUrAQB79uyBwWBAaGioeVAEEUmDJ67yZHMy3blzJ5KSkqDT6ex6A61Wi5iYGPTp0we7du3CxIkT7Q6SlMHDwwOrV6/G4sWLERQUhGPHjuGBBx5wdVhEqsMTV/mxOZlOnz5d1Bs9+OCDTKQq0dh8sdozY54hEzkGT1zlh9czJbupfb4YkVLwxFU+5Dm7l4iISEEkSaaff/55nblNtYvf3759G99++y1KSkqkeBsiIiJZkqSZt6qqCsXFxWjRogXefPNN5OTkoFOnThgyZAgiIyOxb98+jBo1Soq3IiIikh1Jaqbjx49HixYtAADdu3fH5s2bERsbi8zMTEyZMgXZ2dlSvA2R4ilpYXIisp3kA5CGDRuGo0ePon///ujevbvUL0+kSEpcmJyIbCcqmd64cQMXL15ERESE+T4fHx8MHjxYbFxEqsIr6hCpm6hT4uTkZCQlJeHMmTPm+z766CNedo3oHlyYnEj9RNVMQ0NDMWLECHTq1Ml8X3x8PI4cOYIDBw4gKipKdIBESseFyeVN7RetJucQlUybN2+OwsLCehOGBw0ahG3btokKjEgtuDB5w1ydzLgICUlBVDINCgrCG2+8geTkZPTs2RO9evVCr1694O/vj2vXrkkVI5GicWHyhjGZkRqISqbp6elYtmwZbt26haysLOzfvx/vvvsumjRpgoSEBKliJFI8LkxOpG6ikmm7du3Qo0cPAMCf/vQnAEBBQQEOHTrES60R3UPOC5O7upmVSA1EJVN/f3+cPXsWXbt2Nd8XEBCAgQMH4sCBA3jsscdEB0ikJnJcmJzNrJbxJIPsISqZxsbG4uuvv8ZPP/2E8ePHAwBOnjyJV199FUOGDJEkQCIiV+BJBtlD9ApIMTExdW736dMHc+fORc+ePcW+NJFF7lJjcJdyEqmB5MsJenh44KmnnpL6ZYnM3KXG4C7lJFIDm1ZA2rlzJ4qKigS/SUlJCd577z3BzycisoQXDiC5sCmZjh49Gv/1X/+F9PR0mEwmu94gMzMTb7/9NsaMGSMoQCKi+5lMJiQkJGDgwIEAai4ckJCQYPfxiUgqNjXzenp6Yvbs2fjyyy/x8ssvY+DAgejevTs6d+6Mpk2b1tm2vLwcFy5cwJkzZ3D06FH069cPixcvho+Pj0MKUGvHjh3YtWsXSkpKEB4ejvnz51tdWWbu3LnIysqqc9/MmTMxduxYh8ZIRNLghQNci/359dncZ+rh4YGnn34aQ4cOxTfffIPPPvsMp0+fRnV1NXQ6HTw8PGA0GmEymdC5c2dERETgjTfecMpcutTUVGzfvh1JSUlo3bo1kpOT8frrr2PdunVWnzN27FjEx8ebb3NeLJEyNHbhgLfeektWU4/UiP359dk9AMnPzw+jR4/G6NGjUVVVhYKCAty6dQsmkwkBAQEIDAx0+g959+7diIuLQ2RkJABgwYIFmDBhgnmlGUt8fHzcfk1UIiXihQNIjuxKplu3bsWVK1ewePFiAIBWq8UDDzzg0pVcKisrkZOTg5deesl8X+vWrREcHIyzZ89aTaZ79+7FV199haCgIERHRyMuLg5ardZZYRORQLxwQMPYBOsadiXT69evo6CgwHx73759GDFihORB2aO2abl58+Z17g8ICEBhYaHF50RFRaFVq1YICAjA2bNnsXnzZpSUlGDKlClOiNjxOMKR1MwdLxxgzz7NJljXsOvi4CUlJYiNjUVxcTEA4KeffnJIUPaorq62+zkjRoxA37590aFDB4wcORLTp0/Hzp07Bb2WnHCEI7mLlStXYsaMGeaWp9DQUMyYMUN1Fw7gPq0cdtVMp06diiVLluC3337DQw89BADYtWsXOnTogA4dOrjk7Mff3x8eHh51aswAUFhYiICAAJteIywsDGVlZSgqKrL6nMWLF5vPeOXaXMIRjuQu5HzhACmpcZ9WSstZbXM5YFusdiXTNm3aYOvWrTh//jx+/PFH7Ny5E19//TV+/fVX3L17F4GBgejYsSPat2+PsLAw9O3bF/7+/sJKYiNvb2907NgRp0+fRnh4OADg2rVryMvLq7MAf0NycnLg4+PTYKxLly6VdVMJRziSO5LjhQOkYss+rSQmkwmJiYnYs2cPgJpadu1lCD087GokdYp7K01GoxEpKSkNbi9oOcFOnTqhU6dOyMnJwcKFC2EymfDbb7/h4sWLuHjxInJycpCRkYEbN26gc+fOGDNmDJ588kkhb2WT0aNHIzk5GWFhYWjVqhXeffdd9OzZE6GhocjPz8f8+fORlJSELl264MqVK8jIyEBERASaNWuG7OxsbNy4EaNHj4ZGo3FYjI4m9QhHtQ1iUFt5SP1s2aeVNKVPjbXse4lam/fZZ58FUNPk0q5dO7Rr1w6DBw82P15SUoILFy7ghx9+gMFgwLRp00QFa83w4cNRUFCAtWvXmhdtSExMBABUVVUhNzcXFRUVAAAvLy/88MMP+OSTT1BZWYng4GCMHz8e48aNc0hsziL1CEc5DGKQMgHKoTxE9rBln76/1ipX7tByJiqZ1vabWqPT6XD9+nU0adJEzNvYZMKECZgwYUK9+4ODg5GRkWG+3bJlywYXc1AqNY5wZAIkd2bLPq2UZOoOc4Mlv2rM/e7cuYPPPvsMcXFxjn4rt1c7knHPnj3mBStq+ySISHnUsk+7w9xghyfT2NhYxMbGOvptCK4f4Xhvs+yAAQPw17/+lf2SRCK4ep+Wihpbzu7n8GRKzuesEY73D3Fn0iRyDDWMWlZLLdsaJlMFcvXIVKUNcSdSGyW2Aqmllm0Nk6kCuXpgjtqHuBPJnZyTZmPUUMu2hNUIsktjQ9zlvqoJEZEjSFYz/fXXXxudKkPK5w5D3KmGq7sTiJREsmSanJyMFStWSPVyJFPuMMRdybjQBZFrSJZMlX7FFbKN1EPcWfuRFhMgKZEajgOSJVMlr2tL9pFyiDsP/kSkhuMAR/OS3dQ+xJ2IyF5MpiSYWoe4i6GG5ioish+TKamGHBKZGpqr3JVSLlpN8sRkSqrBRKZsrkpmXNGLpMBkSqQwcqiBS8nVyUzKFb2U+N0oMWY5YjJ1E9xhLFPi5yKHGriUn5srl6e05aLV9pDDd2MvJcYsR5xn6ia4w1jmrM9FCQuT29PMKtXnZksyc+QAN1tW9PL09DTHSmSNZMl0zJgxUr0UkawJqZXJKWnez5XNrK5enrKhFb10Oh2WLVuG/fv3A1B+X6qzWmGU2NojBcmS6aOPPirVSxHJmtpq+a7sM3T18pQNreil0+mwefNm1VwdyVm/W7XtH7ZinykplrtMZXBkOV3dZyj18pRCWFrRa/jw4UhNTXVZ8zMpj/LaKsjtmUwmJCQkYODAgQBqmt8SEhJgMplcHJm0nFFOW5pZHW3lypWYMWMGQkNDAQChoaGYMWOGoOUphahd0evYsWMAgGPHjmHRokUoLi62uL2zPhdSFtZMSXHc5eLkziinLc2s99fOpCaX5SnvXdFLr9fz6khkF9ZMSVHc5eLkzipnbTOrj49PnfulbmZNS0tDQkIC5syZg4EDB2LOnDlISEhAWlpanVju/d+VnPW5kHpIUjP9/PPPERsba/6BVVdXQ6PR4Pbt2/j3v/+NiIgI6HQ6Kd6K3JwtzZJ+fn5Ojkp6ziynlFcBskaJg1Kc8bkonbuMW7CFJMm0qqoKxcXFaNGiBd58803k5OSgU6dOGDJkCCIjI7Fv3z6MGjVKirciNyeHZklncGY55dLMKjf8XKyPzo6KisKBAwe4BOM9JEmm48ePN//dvXt3LFy4EBcuXMCBAweQnJyMrl27MpmSJGwZ/amGZCp1OW2ZsiKnZlY5UernIkWt0VqLQkJCgluMW7CH5AOQhg0bhqNHj6J///7o3r271C9PElPiBGt3aX7jRdhJCEcvwuHqVavkSlQyXbp0Kfz8/NCzZ0/06NEDQUFB8PHxweDBgyUKjxxNiQdZd2l+c5dykrQcPQrc1atWyZWoZNqzZ09kZGQgLS0NlZWVCA4ORo8ePczJtW3btlLFSVSPUpvf7OUu5STxpF6EwxJXr1olV6KS6ciRIzFy5EhUVVXh3Llz+PHHH3HmzBmsX78elZWVCAoKQlxcHOLi4ty2U5qIpKeECwe4gjNGgcth1So5kqTPVKvVolu3bujWrRueffZZXLp0CZ999hkCAwPxxRdf4MSJE/jnP/8JrVYrxdsRkZtTU9KUctyCs0aBu8u4BXuISqZGoxFnzpxB586d6/TltG/fHm3btkV8fDxefPFFfPLJJ/jwww8xefJk0QGTc6ltHpnaykPKJ+W4BWeNdmd/fn2ikumbb76JK1eu4Pr16+jTpw/+9Kc/oXPnzrhz5w4MBgOAmlrrc889Z+4AJ2Vw5WW57idFApRTeYgcyZm1Rvbn/0FUMu3RowdWrFiBixcvIi0tDdu3b0d+fj48PT0xf/58AMC///1vFBUVoXnz5pIETM4hh/VvpUyAcigPqYecp5Sx1ugaopJp586dsX37dkRFRWHGjBmYMWMGiouL4enpCV9fXwDA2bNn8fHHH+Oll16SJGCyTqpBGXKZRyZVAnTGCEdyL0qYUsZao3OJSqYRERHo3Lkzjh8/bp5X1KxZszrbTJkyBXFxcfD39xfzVmQDqc6K5TCPTMoE6C7r+ZK8a4xCcNSycogezavX6/Hkk082uI0zEumOHTuwa9culJSUIDw8HPPnz7c636msrAzvvPMODh8+DE9PT0RHR2P69OmyGm1sy0HBUYNp5DCPTMoE6C7r+ZIyaoz2YNJUDpuSaXZ2Nrp06SLoDcQ811apqanYvn07kpKS0Lp1ayQnJ+P111/HunXrLG6/du1anDt3DitWrEB5eTmWLl0KX19fTJkyxaFx2qOhg4KjB9PIYR6ZlAnQXdbzJZIztdeybUqms2bNQkREBMaMGYPw8PBGD9gmkwnHjx/H559/jhMnTiA9PV2SYK3ZvXs34uLiEBkZCQBYsGABJkyYYB7Jdq/i4mKkp6dj+fLl6Nq1K4CapuhNmzZh8uTJsqqdWuOMwTSunkcmdQJ0dXmI3J2QpKmkZnubkunGjRuxYcMGLFy4EM2bN0d4eDi6dOmCkJAQNGvWDNXV1SgpKcGVK1eQnZ2NEydOoLCwEL1798bGjRsdWoDKykrk5OTUGeDUunVrBAcH4+zZs/WS6YULFwAAvXv3Nt/Xt29fGI1GXLlyBe3atXNovGI5a3CQHEYESpkA5VAeIrKPkprtbUqmjzzyCFavXo0ff/wRe/bsQWZmJtLT06HRaOpsV11dDV9fXzz22GMYNWoUunXr5pCg72U0GmEymepNvQkICEBhYWG97QsKCqDT6eDp6VlnWwAoLCyUfTJ19uAgV44IdEQCdJcRjlL2p3OhC6LG2TUAqUePHujRoweqqqrw888/45dffkFhYSE0Gg38/f3x8MMP45FHHnFqU2l1dbXo7e8/KbDESv5yuJr3bWZ+fy+vFtDp9ADq9yU2bdoMnp6BDcZ6/+sJvb+xx+zV0GuVl3sDaIbycu86j0sZs5CyuPqzsXS/yWTCkiWJSE3dA0CDiIiBGDYsFv/8Z01/uj0xN/Za9sbW0P2NPWYvKb9POfyehHDGdyP193njRiWApv//v7jXkpIt76vJyMiwLxvJTGVlJYYNG4a3334b4eHh5vvj4+MRHx+Pp556qs72J06cwIIFC5CWlmaunebl5SE+Ph5bt26tVzMtLS3FyJEjAcwEUFubGfr//1wlAcAGAPc29foAmAFAqQsQVAK4CaAF/vicyX5S/jbU+DsjeTIBSASwD4ARgB7ACAArAbhqhbK0//8H1ByfUrB37140bdrU4taSXxzc2by9vdGxY0ecPn3anEyvXbuGvLw88wCjez3yyCMAgKysLPP2p06dgl6vR0hIiNX3yc1d6pK2eqPRiLZt2yA397d7RvOuxJIlwNdf78HFizno0KEjYmJqawz2v56Q+xt7zFb31n4uXcpB+/YdLdR+HB+zkLI4+rOxN+bKykoMHLgPBsP9A7PKERq6D8eOvYXy8nKbYvbx8Wn0tby9vWX3exLzWnL+PQnh6v3GnvIkJSViy5YN94wFyYOPz2VMmQIsW7Za8s/GNn9UmmreP6XBrRWfTAFg9OjRSE5ORlhYGFq1aoV3330XPXv2RGhoKPLz8zF//nwkJSWhS5cu0Ov1GDJkCNavX4+FCxeivLwcW7ZswahRoxpsntbra/65RvF97++BlJTVuHGjpi/x++/t7Uu8//Xsu7/uEPduWLFC+BD3hIS6O9GlSwZs2bIBTZrcPzJZXMzin2ONkPcRovGYr127iZISy+1RpaXFuHv3FvR6P5tiLi1t/LUeeCC4znMc+90IIeX3KYffkxCu3m8aL09lZSXS0y0PqkxP34c1a2oXaJH6s5GWKpLp8OHDUVBQgLVr15oXbUhMTAQAVFVVITc3FxUVFebt582bh3Xr1iExMRFarRbR0dGKvKKNqwbTSDUs3VEjk+UwYMYVMUg5N1cOC3dITUnTLNyJWlYoU0UyBYAJEyZgwoQJ9e4PDg5GRkZGnft8fX2xaNEiLFq0yFnhkQVSj0yWw5VhXBmDlHNz5bBwh9SUNM3CnahlhTLVJFNSHqlrP3K4MowzY7BU+5Vybi4XuhBODq0jSqGWFcp4IUc3I6edvHYn8vHxqXO/kNpPY03Gziivs2IwmUxISEjAwIEDAdTUfhMSEmAymcxzc48dOwYAOHbsGFavXi2oVmzLa8np9yQHDX03cufK73LlypWYMWOGeZGd0NBQzJgxQ1EnbpLWTE+fPo3i4mKEh4croo3bncihCdQSqWo/YpuMpTiQOGtBDVtqv1L2p1t6Lbn+nlytoe9m6NChsuyzlcN3qYYVygQl023btiErKwurVq0y3/faa68hMzMTQE3zXUpKCoKCgqSJkkSTQxOoJVLtREKbjKU8kDhj0I7arjWrJrZ8N3Los73/xFFO36WSVygTdNpx6NAhPPTQQ+bbJ06cwJEjRzBy5EjMmTMHt2/fxr/+9S/JgiRx5NAE2hixO5HQJuPaA4nBYADwx4GkdjR4Q+4/KEnZbG2NLbVfR1PC78kV5PDdNMRSE/ScOXOwd+9eh3yX7tYFIKhmmp+fj7Zt25pvZ2ZmomXLlpg3bx4A4Pfff8ehQ4ckCZDEU8vQ88bY22QstJbXUG3W0YN25DBlxZbfU+3qYmIOpEqbyiKH76YhlmqgmzZtsjq/XmjXhByajV1BUDK9c+dOnYPMyZMn6yzl16ZNG9y8eVN8dCQJtQw9b4y9TcZC+zgbaxZzZN+PHKasNPR70ul0WLZsGfbv3w9A3IFUaVNZ5PDdWGPtxLGiogJeXl4WnyP0BEBOzcbOJCiZtmzZEhcvXgRQs67t5cuXMXHiRPPjhYWFaNKkiTQRkmhqGXpuK1ubjIXUJGytzTqy78fVU1Ya+j3pdDps3rzZrgOpqy8aLeX7u/q7saahE0cvLy94eHjUWdhG6AmAXPr0XUFQMh08eDA+/vhjmEwmZGdnw8/PDwMGDDA/bjAY0Lp1a8mCJPHkupO7kpCahLMvgWeJHEY+Wvo9DR8+HKmpqQ0eSC31o7m62VbK95fDd2NJQyeOISEhGDZsGPbv3y/62CCH/cNVBDVgT5w4EUOGDMG3336L4uJiLFmyBDqdDgBQUlKCzMzMOs2+5HpSzj9UE3vnt9UelCxxdr+YHK41e+/vadGiRSguLra4vdFoxOzZsxU5/1IIuY1KbWhw3MiRI7Fu3TpJjg1y2j+cTVDN1NvbGwsXLrT4mJ+fH3bu3FnvSyN5kNtO7moN1SSsDYDp1q2bpP1irh5oI6aZ897fk16vt1r7KS8vx7Zt29yuH01OGmudkuLYIOd+Y0eTfDlBDw8Pcy2VSCksHUisDYC5d7SiFE3mrh5oI1XStnYgrR0/obR+NFf35UpNTBO0PdNc3LVLyaZkmpWVJejFe/XqJeh55HzuNidMDLn2i9nC0d+zpQPpE088YZ4mcT8596MpNWk2xp4aqJBpLkreP8SwKZnOmzcPGo3G5hetrq6GRqPBwYMHBQdGzqHEOWFyqTE4uslcynI663u2dCDV6/X47//+b9nOvyTrxExzcbcuJZuS6YIFCxwdB7mIEueEqbXGcD8py+ns7/neA6k796MpmS3TXOgPNiXTmJgYR8dBLsCdxT3IYe6fu/ajKZm7rJwmFXm245FTyH0tUZKGHL5nTs1SHnee5iKEqNG8t27dwvnz51FSUmJxvpg7NMUpmTOXGXT19A93Jqc1Y92tH80echsE6G4rp4klKJlWVVVhzZo1+Prrr1FdXW11Ox4k5c2ZO4urp3+4M/ZZypucBwFaa56PiopCQkKCywcByomgZPrRRx9h//79iI6ORnh4OJYtW4Zp06bBz88Pu3btgq+vL6ZOnSp1rOQA7MtyD/yeHUOK2qScBwE2NM1l2LBhLo1NbgSd9hw4cAADBgzAokWLEBERAQAICwvDU089hU2bNqG0tBTnzp2TNFByDPZluQd+z9KydG1QIcsjOuLasGlpaUhISMBf//pXc60xISEBaWlpdr9WLTbPN05QzTQvLw+jR48GAPPOePfuXQA1q50MHToUqampeO6556SJkhyOO4t74PcsDalqk45YGN6dm1pdSdBpqY+PjzmJ+vn5QaPR1BkRqNfrcf36dWkidHNyG5RApASO3G+krE1yxKx6CEqmrVu3xm+//QYA0Gq1aNeuHQ4fPgygZvWj7777zi2Wj3IkqZqRiNyJM/YbKacaNXQ1Fw4OUxZBzbz9+vVDWloaZsyYAQ8PD8TGxiI5ORkTJkwAUNMM/OKLL0oaqLuR86AEIrlyxn4j9ZQyDg5rnBJa6ATVTJ977jmsXLnSPC1mzJgxeOmll9C0aVM0a9YMU6ZMMSdWsp8jBiUQqZ2z9hupa5McHGadklroBNVMfX190a5duzr3PfPMM3jmmWckCcrdyeFq9XJZTF4qaisP1efM/cYRtUkODqtPSS10kl/PlMSTw4o1ckgyUiZAOZSHHMuZ+427XmbMmeSwprQ9BCXT5cuXN7qNRqPh1WYE4oo1NZgAyR6u2G9Ym3QcObTQ2UNQMrVl8i+TqTgclEBkP+436iGHFjp7CEqm3377bb37TCYTrl27hs8++wwXLlzA22+/LTo4d8ZmJCL7cb9RD6W10EnWZ+rh4YGQkBDMnTsXf/3rX7Fx40YkJiZK9fJuy1IzEgfTEDWsoeZXJUyzoBpKamlwyACkfv364YMPPmAydRAmTSL7yfnqLGSZkloaHJJM8/PzUVFR4YiXJiISREnTLKguJQz0EpRMf//9d4v3FxcX49SpU/j888/RrVs3UYERkWupqTtBadMsSHkEJdP4+HhoNBqLj1VXV6NNmzZ45ZVXRAVGRPbj3FzL5DTNgn226iQomU6aNMliMtXr9QgJCUG/fv2c2geRmpqK7du34+bNm+jcuTMSExPRtm1bq9u/9dZb9ab3xMXFYdasWY4Olcih1JQApSSHaRbss7VODa0ggpLpCy+8IHEYwp08eRKrVq3C3Llz0a1bN2zbtg1JSUn44IMP4OXlZfV5kZGRmDNnjvn2/etsknOoYSci+ZPDNAv22Vqnhv1d8csJfvHFFxg8eDBGjhwJAFiwYAGefvppfP/99xg0aJDV53l7e8tu0q87UsNORMrgymkWtvTZkrLZlEwTEhIEvbgzzrays7Pr1JR9fX3RpUsXZGdnN5hMv//+e4wePRoBAQEYNGgQnn/+eTRp0sTh8ZK6saYtPan6GF05zcKWPls/Pz+nxEKOYVMyvXr1ar0+0vLychQVFQEAdDodqqurUVpaCgDw9/eHr6+vxKFaVlhYiObNm9e5z9/fHwUFBVafExERgSeeeAItW7bEpUuXsGnTJly/fh2LFy92dLjkQHJIZEya0nFUH6MrpllIfQ1Ukh+bkunHH39c5/bVq1eRkJCAuLg4xMfHm5tLb926hR07duDIkSNYtWqVqMBWr15t3oks6dWrF9auXSvotZ988knz3+3bt0dAQADmz5+Pl19+GQEBAYJek1yPiUxd1NTHaEufLZOpsgnqM01JSUG3bt0wc+bMOvcHBgZi1qxZuHXrFlJSUvCPf/xDcGDTpk3D888/b/Xx2sFFAQEB9WqhRUVFCAkJsfm9wsLCAAB5eXlWk+nixYvNZ7I8aJMryaEG7mhKnRfa0Hcj16Xx3OH3JETt5wLY1sUgKJlmZWVh2rRpVh/v3bs3Nm/eLOSlzXQ6HXQ6XaPbdenSBadPn8aIESMA1Oxs2dnZiIuLs/m9cnJyAKDBeWZLly6FXq+3+TXlhjuMcHKbFyiH78zRvyc5zQu1R2Pll+PSeHL4PcnRvZ+L0WhESkpKg9sLHs17+fJlq4/98ssvQl/WbqNGjcLChQvRu3dvdO3aFdu3b0eLFi3Qv39/8zaTJk3C1KlT8fjjj6OsrAzbtm1DZGQkmjdvjkuXLiE5ORmDBw9WdRMvdxjLGkoKUVFRTp0XKLek3RBH/57kMC/UUZSwNB7ZT1Ay7devH7766it07twZ//Ef/1HnsQMHDmDv3r149NFHJQmwMeHh4UhISMC2bdtw69YtdOnSBcuWLaszxzQ3N9c8OMrDwwMGgwGpqam4ffs2goKCEBkZiUmTJjklXpKXhpJCQkKCU/rsOJm/PjnMC1Ubtk45lqBkOnPmTJw/fx7Lli3D5s2b0aZNG2g0GuTm5uLmzZsICgqq15/qSMOHD8fw4cOtPp6RkWH+u0mTJlixYoUzwiIFa6jP7tNPP0VVVRUASHJQUtNAGynJtY9RqZyVNN01aQtKpkFBQXjvvffw0UcfITMzE//7v/8LAGjVqhWioqIQHx9vU38nkVw11GdnMpmQlJQkSZ+dUgfaWCPlgVRJl9+iP6g9aVojuM9Up9Nh6tSpmDp1qpTxEMmCs/rslDrQxhpHHEjZx0hK4J4dMkSNqO2zu3/NZqn77GqTtiVKH2hD5E5sqplmZWUBqFko4d7bjandnkiJnNFn5y4Dbdy1H43ch03JdN68edBoNPj666/h5eVlvm1NdXU1NBoNDh48KFmgpA5Kmv7hrD47dxhow6RJamdTMl2wYEHNxp6edW4T2UrJ0z8c3WfHgTZEymdTMo2JiWnwNlFjOP2jcY5M2mxmJSH4u7GdpNczPX36NIqLixEeHs7LCSmEM3YWtU3/UCIe/EgI/m5sJyiZbtu2DVlZWXWuDPPaa68hMzMTQM0IxZSUFAQFBUkTJTmMM3YWtU3/ICK6n6DOqkOHDuGhhx4y3z5x4gSOHDmCkSNHYs6cObh9+zb+9a9/SRYkKRunfxCR2gmqmebn56Nt27bm25mZmWjZsiXmzZsHAPj9999x6NAhSQIk5XOX6R9E5L4EJdM7d+7UOQCePHkS4eHh5ttt2rTBzZs3xUdHquEO0z+IhOAgH3UQlExbtmyJixcvAqi5oPbly5cxceJE8+OFhYVo0qSJNBGSKnD6B5FlTJrqICiZDh48GB9//DFMJhOys7Ph5+eHAQMGmB83GAxo3bq1ZEGSenCdVSJSI0HJdOLEicjPz8e3334LnU6HJUuWmK8SU1JSgszMTIwdO1bSQImI5IJNs3Q/QcnU29sbCxcutPiYn58fdu7cWW+BcCIitWDSpPtJumgDUNM3xmuZEpErsMZIriI4mebl5WHr1q04fvw4CgsLsXz5cvTt2xeFhYXYtGkTRo0ahc6dO0sZq1vhQYHIftw/yFUEJdMrV65g5syZuHv3Lrp06YKTJ0+aHwsICMDPP/+MvXv3MpmKwIMCEZFyCEqm7733Hry8vPDee+/B29sbTz/9dJ3H+/fvjyNHjkgSIBERuScltdAJSqanTp3CuHHjEBQUhKKionqPP/jgg1y0gYiIRJFj0rRG0Nq85eXlCAgIsPp4RUUFqqqqhMZERESkKIKSaatWrfDzzz9bfTwrKwvt2rUTHBQREZGSCEqmTzzxBNLS0nD69GnzfRqNBgCwe/duHDt2DFFRUZIESEREJHeC+kzj4+Nx8uRJzJ8/Hw8//DA0Gg02bNgAo9GI/Px89OnTp96gJCIiIrUSVDP19vbGqlWr8NJLL0Gr1cLb2xuXL19G06ZNMXXqVLz11lvQarVSx0pERCRLdtdMKysr8cMPPyA4OBjjx4/H+PHjHREXEZGklDTNgpTH7mTq6emJv//975g5cyY6duzoiJiIiCTHpEmOZHczr4eHB4KCglBRUeGIeIiIiBRHUJ9pVFQUDhw4gDt37kgdDxERkeIIGs3bs2dPZGZmYurUqRg1ahRCQkLQpEmTetv16tVLdIBERGrAPlt1E5RMX331VfPf69evN88xrVVdXQ2NRoODBw+Ki46ISCWYNNVNUDJdsGCB1HEQEREplqBkGhMTI3UcREREiiX44uBERI7EPkZSEiZTIpIlJk1SEsUn06ysLHz00Uc4d+4cioqK8OGHHyIkJKTB51RVVWHjxo345ptvcOfOHTz++OOYO3cufH19nRQ1ERGpiaB5pnJSXl6OsLAw/OUvf7H5Odu2bcPBgwfx2muvYdWqVTh//jzWrFnjwCiJiEjNFF8z7d+/P/r374+8vDybtjeZTPjyyy/xl7/8BeHh4QCAV155Ba+++ipmzpwJf39/R4ZLCuKsPjv2DRIpn+KTqb2uXbuGoqIi9OnTx3xf7eIS58+fR0REhKtCI5lxVjJj0iRSPsU389qroKAAANC8eXPzfVqtFnq9HoWFhS6KioiIlEy2NdPVq1djz549Vh/v1asX1q5da/frVldXi4iKiIioPtkm02nTpuH555+3+riXl5eg1w0MDARQU0P18/MDUDO612g0IiAgwOrzFi9eDG9vbwBsliMiUrvasQxAzXW8GyPbZKrT6aDT6SR/3VatWsHf3x+nT582T6E5c+YMAKBTp05Wn7d06VLo9XrJ4yEiIvm5t9JkNBqRkpLS4PaK7zMtKyuDwWDAL7/8AgD49ddfYTAYYDQazdtMmjQJ3333HYCa67E+9dRT+OCDD3Dy5ElkZ2dj/fr1GDJkCEfyEhGRILKtmdrq/PnzmDdvnvn2kiVLAAALFy40ryGcm5uL0tJS8zaTJ09GWVkZXn/9ddy5cweDBg2q8xokPU7/ICI102RkZHBETgNKS0sxcuRIFBUVsZmXiMgNGY1G+Pv7Y+/evWjatKnFbRTfzEtERORqTKZEREQiMZkSERGJxGRKREQkEpMpERGRSEymREREIjGZEhERicRkSkREJBKTKRERkUhMpkRERCIxmRIREYnEZEpERCQSkykREZFITKZEREQiMZkSERGJxGRKREQkEpMpERGRSEymREREIjGZEhERicRkSkREJBKTKRERkUhMpkRERCIxmRIREYnEZEpERCQSkykREZFITKZEREQiMZkSERGJxGRKREQkEpMpERGRSEymREREIjGZEhERicRkSkREJBKTKRERkUhMpkRERCIxmRIREYnEZEpERCSSp6sDECsrKwsfffQRzp07h6KiInz44YcICQlp8DlvvfUW0tLS6twXFxeHWbNmOTJUIiJSKcUn0/LycoSFhWHQoEFYtWqVzc+LjIzEnDlzzLd9fHwcER4REbkBxTfz9u/fH1OmTEG/fv3sep63tzcCAwPN//z8/BwUoXj316Ldkbt/Biy/e5cf4Gcg9/IrPpkK9f3332P06NF44YUX8P7776OiosLVIVkl9x+RM7j7Z8Dyu3f5AX4Gci+/4pt5hYiIiMATTzyBli1b4tKlS9i0aROuX7+OxYsXuzo0IiJSINkm09WrV2PPnj1WH+/VqxfWrl0r6LWffPJJ89/t27dHQEAA5s+fj5dffhkBAQF1tq2urgYAGI1GQe8lhcrKSpe+vxy4+2fA8rt3+QF+Bq4sf+371uYDSzQZGRnWH3WhkpISlJWVWX3cy8urTuLLy8tDfHy8TaN5Lb1XbGwsNmzYgM6dO9d5LD8/H+PHj7fr9YiISH0+/fRTBAUFWXxMtjVTnU4HnU7nlPfKyckBAAQHB9d7rEWLFvj000/h6+sLjUbjlHiIiEg+qqurUVZWhhYtWljdRrbJ1FZlZWW4cuUKbty4AQD49ddfUVZWhpYtW0Kv1wMAJk2ahKlTp+Lxxx9HWVkZtm3bhsjISDRv3hyXLl1CcnIyBg8eXK+JFwA8PDysnokQEZF7aKxyp/hkev78ecybN898e8mSJQCAhQsXIiYmBgCQm5uL0tJSADXJ0WAwIDU1Fbdv30ZQUBAiIyMxadIk5wdPRESqINs+UyIiIqVw23mmREREUlF8M69a7NixA7t27UJJSQnCw8Mxf/58BAYGWty2rKwM77zzDg4fPgxPT09ER0dj+vTp0Gq1To5aOraW32g0YsuWLfjhhx9w48YNPPDAAxg6dCgmTJjgFuW/V2lpKf785z/j999/R3p6uqLLD9j/GRw8eBA7duxAbm4u9Ho9xo4di2effdaJEUvLnvJfunQJ7777LrKzs6HVatGrVy/MnDkTDz74oJOjlsbhw4fxxRdf4MKFCygtLW309yzHYyBrpjKQmpqK7du345VXXkFycjJKS0vx+uuvW91+7dq1OHv2LFasWIG//e1vyMjIwNatW50YsbTsKf/NmzdRWFiI2bNnY8uWLZg5cyZ2796NDz/80MlRS8fe77/WO++8g3bt2jkhQsez9zP45ptv8M4772DcuHH44IMPsHTpUnTq1MmJEUvL3vL/53/+J3Q6Hd59912sWrUKJSUl+Mc//uHEiKVVUVGBvn37Ij4+3qbt5XgMZDKVgd27dyMuLg6RkZEIDQ3FggULcObMGRgMhnrbFhcXIz09HbNnz0bXrl3Rt29fTJkyBV9++SWqqqpcEL149pS/ffv2+Pvf/44BAwYgJCQEjz76KMaOHYvMzEwXRC4Ne8pf67vvvsPly5fxzDPPODFSx7HnM7h79y42btyIGTNmICYmBiEhIQgLC0OfPn1cELk07Cl/YWEhrl69igkTJqBdu3YIDQ3F2LFjceHCBRdELo2oqChMnDgR3bp1a3RbuR4DmUxdrLKyEjk5OXUOBK1bt0ZwcDDOnj1bb/vaHaZ3797m+/r27Quj0YgrV644PF6p2Vt+S4qKitCsWTNHhehQQsp/69YtJCcnY9GiRYpv2gWE7QMFBQWoqqrCiy++iPHjx2PZsmUoKipyZtiSsbf8er0ebdq0wTfffIPKykqUlZXh4MGDdl/sQ6nkegxkMnUxo9EIk8mE5s2b17k/ICAAhYWF9bYvKCiATqeDp6dnnW0BWNxe7uwt//2uXr2K/fv3Y8SIEQ6K0LGElH/VqlUYM2YMHnroISdE6Hj2fgZ5eXkAavoYX3rpJbz22mu4fPmyYps57S2/h4cHVqxYgePHj2PYsGEYMWIErl69iqSkJCdF7FpyPQYymbpYQ2s92rq9kldmsrf89yooKMCiRYvw5JNP1llvWUnsLX9qaiqKioowbtw4B0XkfPZ+BiaTCQDw/PPPY8CAAejevTvmz5+P48eP4/r1644I0aGElH/t2rV46KGHkJKSgnXr1sHPz0+xJxP2kusxkKN5Xczf3x8eHh4oKCioc39hYaHFFZkCAwNRUlKCu3fvms/Map9raXu5s7f8tYqKipCYmIhOnTph7ty5jg3Sgewtf1ZWFrKzsxEVFVXn/ujoaMydOxexsbGODNch7P0Mamtw9w6+qv37+vXraNmypeOCdQB7y3/q1CmcOnUKe/bsgbe3NwAgKSkJ48aNw8WLF9GhQwdnhO0ycj0GsmbqYt7e3ujYsSNOnz5tvu/atWvIy8tD165d623/yCOPAKg5qNY6deoU9Hq93Qv8y4G95QdqBiC8+uqraNWqFRYtWgQPD+X+jO0t/5///Ge8//775n+JiYkAgE2bNmHw4MFOilpa9n4GnTp1gqenZ53+sdq/lTg1xN7yl5eXQ6PR1Pnd1/5dW2tXM7keA5V7FFKR0aNH4/PPP8d3330Hg8GAFStWoGfPnggNDUV+fj4mTZqE7OxsADWDD4YMGYL169cjOzsbp06dwpYtWzBq1CjFDkaxp/ylpaVYsGABtFotZs+ejaKiIty6dUuR/cW17Cl/UFAQ2rdvb/7XqlUrADWjnJU6CAuw7zPQ6XQYOnQoPvjgA5w5cwY5OTlYu3Yt+vfvr9h1tO0pf7du3eDl5YVVq1bh8uXLyMnJwcqVK9G6dWvF9qMbjUYYDAbzSZHBYIDBYEBZWZlijoFs5pWB4cOHo6CgAGvXrjVP2K6tcVRVVSE3NxcVFRXm7efNm4d169YhMTERWq0W0dHRmDx5sqvCF82e8v/88884d+4cANSZoP/ggw/i448/dn7wErD3+1cjez+D2bNnIyUlBUuWLIFWq0VERARmzZrlqvBFs6f8AQEBeOutt7B582a8/PLL0Gq16Nq1K5YtWwYvLy9XFkOwo0ePYvny5ebb06dPBwCsWbMGwcHBijgGcm1eIiIikdjMS0REJBKTKRERkUhMpkRERCIxmRIREYnEZEpERCQSkykREZFITKZEREQiMZkSERGJxGRKREQkEpMpERGRSEymREREIjGZEhERicRkSkREJBIvwUZEFn311VcwGo345ZdfEB0djfz8fNy6dQsGgwGzZs1S7LVDiRyBNVMiqic1NRUdO3bExIkTMWrUKLzxxhvw8vJCp06dcPjwYVy8eNHVIRLJCmumRFSP0WhEt27dAAA3b96ERqPBkCFDcOfOHaxZswa9e/d2bYBEMsNkSkT1PPPMM+a/z5w5g549e0Kr1UKr1TKRElnAZl4iatCpU6eYQIkawWRKRHVUVVXh+PHjqKqqQkFBAX755Rf07NnT/PiOHTtcGB2RPDGZElEd+/btw6uvvorff/8dGRkZaNKkCfz9/QEAR44cwcMPP+zaAIlkSJORkVHt6iCISD4uXbqEDz/8EO3atUPHjh1x+/ZtnDhxAq1bt0ZwcDCGDh3q6hCJZIfJlIiISCQ28xIREYnEZEpERCQSkykREZFITKZEREQiMZkSERGJxGRKREQkEpMpERGRSEymREREIjGZEhERicRkSkREJBKTKRERkUj/B9ZSKjGGC2G7AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.errorbar(x,r,u, fmt = 'ko')\n", "plt.xlabel('x')\n", "plt.ylabel('residuals ($y_i - f(x_i)$')\n", "plt.axhline(0)\n", "plt.xlabel('$x$')\n", "plt.ylabel('residuals ($y_i - f(x_i)$)')\n", "plt.xlim(-0.1, 1.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### (d) Determine the goodness-of-fit parameter $\\chi^2$ for this data set and model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Goodness of fit\n", "\n", "$$ \\chi^2 = \\sum_i\\left(\\frac{y_i - y(x_i)}{\\alpha_i}\\right)^2 $$ " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "36.61459756007561 0.7472366848995022 51\n" ] } ], "source": [ "chi2 = np.sum(r_norm**2) # Chi-square\n", "chi2_nu = chi2/(len(r)-2) # reduced chi-square\n", "print(chi2, chi2_nu, len(r))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###### Conclusions:\n", "\n", "#### (e) Given your analysis, does your fit to the data appear to be reasonable? Comment briefly. & (f) What are your resultant values for $a_1$ and $a_2$? (Include uncertainties.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The goodness of fit parameter $\\chi^2$ is less than the number of data points (but not much less), or, equivalently\n", "the reduced chi-square $\\chi^2_\\nu$ is less than 1 (but not much less than one), so the data is consistent with the \n", "assumed model. Also, in Fig.2 the residuals fluctuate around $0$ and show no systematic dependence on $x$.\n", "\n", "Therefore we can use the best-fit parameters found above:\n", "\n", "$$ a_1 = 1.96 \\pm 0.08\\quad\\quad \\mbox{and}\\quad\\quad a_2 = 0.88 \\pm 0.08 $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Version information\n", "`version_information` is from J.R. Johansson (jrjohansson at gmail.com); see Introduction to scientific computing with Python for more information and instructions for package installation.\n", "\n", "`version_information` is installed on the linux network at Bucknell" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "%load_ext version_information" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.7.8 64bit [GCC 7.5.0]" }, { "module": "IPython", "version": "7.17.0" }, { "module": "OS", "version": "Linux 3.10.0 1127.19.1.el7.x86_64 x86_64 with centos 7.9.2009 Core" }, { "module": "numpy", "version": "1.19.1" }, { "module": "scipy", "version": "1.5.0" }, { "module": "matplotlib", "version": "3.3.0" } ] }, "text/html": [ "
SoftwareVersion
Python3.7.8 64bit [GCC 7.5.0]
IPython7.17.0
OSLinux 3.10.0 1127.19.1.el7.x86_64 x86_64 with centos 7.9.2009 Core
numpy1.19.1
scipy1.5.0
matplotlib3.3.0
Thu Feb 10 11:15:41 2022 EST
" ], "text/latex": [ "\\begin{tabular}{|l|l|}\\hline\n", "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n", "Python & 3.7.8 64bit [GCC 7.5.0] \\\\ \\hline\n", "IPython & 7.17.0 \\\\ \\hline\n", "OS & Linux 3.10.0 1127.19.1.el7.x86\\_64 x86\\_64 with centos 7.9.2009 Core \\\\ \\hline\n", "numpy & 1.19.1 \\\\ \\hline\n", "scipy & 1.5.0 \\\\ \\hline\n", "matplotlib & 3.3.0 \\\\ \\hline\n", "\\hline \\multicolumn{2}{|l|}{Thu Feb 10 11:15:41 2022 EST} \\\\ \\hline\n", "\\end{tabular}\n" ], "text/plain": [ "Software versions\n", "Python 3.7.8 64bit [GCC 7.5.0]\n", "IPython 7.17.0\n", "OS Linux 3.10.0 1127.19.1.el7.x86_64 x86_64 with centos 7.9.2009 Core\n", "numpy 1.19.1\n", "scipy 1.5.0\n", "matplotlib 3.3.0\n", "Thu Feb 10 11:15:41 2022 EST" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "version_information numpy, scipy, matplotlib" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.8" } }, "nbformat": 4, "nbformat_minor": 2 }