{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Least squares fitting using ``linfit.py``" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook demonstrates the function ``linfit``, which I propose adding to the SciPy library. \n", "\n", "``linfit`` is designed to be a fast, lightweight function, written entirely in Python, that only calculates only as much as the user desires, and no more. It can handle arbitrarily large data sets." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What is ``linfit`` and what does it do?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``linfit`` is a function that performs least squares fitting of a straight line to an $(x,y)$ data set. It has the following features:\n", "\n", "* ``linfit`` can fit a straight line $ax+b$ to unweighted $(x,y)$ data where the output is the slope $a$ and y-intercept $b$ that minimizes the square of the residuals $$E = \\sum_{i=1}^n \\left[ y_i-(ax_i+b) \\right]^2$$ where $n$ is the number of data points.\n", "\n", " - Optional outputs:\n", " \n", " - $\\Delta a$ and $\\Delta b$, the respective uncertainties in the fitted values of $a$ and $b$. By default, these estimates of $\\Delta a$ and $\\Delta b$ use the residuals $|y_i-(ax_i+b)|$ as estimates of the uncertainties $\\sigma_i$ of the data (\"relative weighting\").\n", " - $y_i-(ax_i+b)$, the residuals.\n", "\n", "* Alternatively, ``linfit`` can fit a straight line $ax+b$ to $(x,y)$ data that is weighted using estimates of the uncertainties of the data provided as an optional keyword argument. The uncertanties can be expressed as either as (1) a single number $\\sigma$ or (2) as an array of uncertainties $\\sigma_i$. In the first case, the output is the slope $a$ and y-intercept $b$ that minimizes $\\chi^2$, which is defined as $$\\chi^2 = \\frac{1}{\\sigma^2}\\sum_{i=0}^n \\left[ y_i-(a x_i + b) \\right]^2$$In the second case, the output is the slope $a$ and y-intercept $b$ that minimizes $\\chi^2$, which is defined as $$\\chi^2 = \\sum_{i=0}^n \\left[ \\frac{y_i-(a x_i + b)}{\\sigma_i} \\right]^2$$\n", "\n", " - Optional outputs:\n", "\n", " - $\\chi^2/(n-2)$, the reduced value of $\\chi^2$, which should be approximately equal to 1 if a straight line is a good model of the data and good error estimates $\\sigma_i$ are provided.\n", " \n", " - $\\Delta a$ and $\\Delta b$, the uncertainties $\\Delta a$ and $\\Delta b$ in the fitted values of $a$ and $b$. By default, these estimates of $\\Delta a$ and $\\Delta b$ use the residuals $|y_i-(ax_i+b)|$ as estimates of the uncertainties $\\sigma_i$ of the data.\n", " - $y_i-(ax_i+b)$, the residuals.\n", "\n", "``linfit`` can also be used for nonlinear fitting for functions that can be transformed to be linear in the fitting parameters. In this case, using weighting is essential to get the fit right." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Demonstrations of ``linfit.py``" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import libraries we will need." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.gridspec as gridspec # for unequal plot boxes\n", "from linfit import linfit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Very simple demonstration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perform a simple linear fit without any weighting." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slope = 1.00, y-intercept = -0.95\n" ] } ], "source": [ "x = np.array([0, 1, 2, 3])\n", "y = np.array([-1, 0.2, 0.9, 2.1])\n", "fit, cvm = linfit(x, y)\n", "print(\"slope = {0:0.2f}, y-intercept = {1:0.2f}\".format(fit[0], fit[1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the above data and fit." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deVyU9fr/8ddHBBFB3BVFQM0N3EVxLa0srdzrpMeTZSfNOp3Szi+XXLKyo3Y6LR47mZ3ULMtK3NM0K81SEy1lUVFcQUgUlR1kuX5/gHzJQEEGZoa5no/HPJy578/Mfc0dvbm5789cY0QEpZRSlV8VaxeglFKqYmjgK6WUg9DAV0opB6GBr5RSDkIDXymlHERVaxdwI/Xq1RM/Pz9rl6GUUnbjwIEDF0WkflHrbDrw/fz82L9/v7XLUEopu2GMOVPcOj2lo5RSDkIDXymlHIQGvlJKOQibPodflKysLGJiYsjIyLB2KRXG1dUVb29vnJ2drV2KUsqO2V3gx8TE4OHhgZ+fH8YYa5dT7kSEhIQEYmJiaNasmbXLUUrZMbs7pZORkUHdunUdIuwBjDHUrVvXof6iUUqVD7sLfMBhwv4aR3u/SqnyYZeBr5RSlVXI6Uss3nmiXF5bA7+M5syZwxtvvFHs+nXr1nH48OEKrEgpZY9SMrOZvT6chxbv4dOfz5J2Ndvi26j0gZ+4cSPH77yLI239OX7nXSRu3Fih29fAV0rdzM5jF7j3rR/4eO8ZxvX2Y8tzfXFzsfycmkod+IkbNxI3azbZsbEgQnZsLHGzZpc59F977TVat27N3XffTWRkJAAffPAB3bp1o2PHjowcOZK0tDR2797Nhg0beOGFF+jUqRMnTpwocpxSyjFdTr3K818c5NGl+6ju4sTqib14aXAANaqVzwTKSh348W+9jVw3u0UyMoh/6+1bfs0DBw6watUqfv31V9asWUNISAgAI0aMICQkhEOHDtG2bVs+/PBDevXqxZAhQ/jXv/7FwYMHadGiRZHjlFKORUTYHBbHgLd2suFgLH+/8za+erYPXX1rl+t27W4efmlkx8WVanlJ7Nq1i+HDh+Pm5gbAkCFDAAgPD2fmzJlcuXKFlJQU7r333iKfX9JxSqnKKT4pg1nrw9kacZ72TTxZ8XgQ/o1rVsi2K3XgV/XyyjudU8TysihqmuRjjz3GunXr6NixI8uXL2fHjh1FPrek45RSlYuI8OWBGOZuOkxmdi7TBrXhiT7NqOpUcSdaKvUpnQaTJ2FcXX+3zLi60mDypFt+zdtvv521a9eSnp5OcnIyG/OvByQnJ+Pl5UVWVhYrV64sGO/h4UFycnLB4+LGKaXs180mh0RfSuORD/cxZXUobRrVZMtzfZl4R4sKDXuo5Ef4noMHA3nn8rPj4qjq5UWDyZMKlt+KLl268PDDD9OpUyd8fX3p27cvAK+++ipBQUH4+vrSvn37gpAfNWoU48ePZ+HChaxevbrYcUop+3Rtcsi164XXJocAuN//AB/tPs2/tkbiVMXw6rB2jOnuQ5Uq1vkwpRERq2y4JAIDA+X6L0A5cuQIbdu2tVJF1uOo71spW3f8zruKPHVMw0bMeHguv5y9Qr/W9fnn8PY0rlW93OsxxhwQkcCi1lXqI3yllCpvxU0CyT3/G6cupvL2w50Y2qmxTbRI0cBXSqkyKG5ySIpnPb55/g7quVezQlVFq9QXbZVSqrwVNTkkt1o12sycYlNhDxYIfGNMU2PM98aYI8aYCGPMc0WMMcaYhcaYKGNMqDGmS1m3q5RStuBoQC8+7DGa89VrIRicvLzwnvtqmSaHlBdLnNLJBv4hIr8YYzyAA8aYb0SkcAOZQUDL/FsQ8F7+v0opZZeSM7JY8PVRPtl7Fp/m3Rj+/x7H/7Z61i7rhsoc+CISB8Tl3082xhwBmgCFA38osELypgTtNcbUMsZ45T9XKaXsyvdH43lxbRjnkzJ4ok8znr+nVbk0O7M0i57DN8b4AZ2Bn69b1QSILvQ4Jn9ZUa8xwRiz3xiz/8KFC5Yszyb169eP66eeKqVs06XUq0xa9SvjlofgXq0qwU/1YuYD/nYR9mDBwDfGuAPBwCQRSbp+dRFPKfIDACKyREQCRSSwfv36lirP6nbs2MFjjz1m7TKUUrdARNh4KJYBb+5kU2gcz93Vkk3P9qGzT/k2O7M0iwS+McaZvLBfKSJrihgSAzQt9NgbKOKTCvZh2LBhdO3alYCAAJYsWcKZM2do2bIlFy9eJDc3l759+7Jt27YSv94nn3xCr169aNeuHfv27SvHypVSpXU+KYPxKw7w989+pUnt6mx6tg+TB7SiWlUna5dWamX+O8TkfZrgQ+CIiLxZzLANwDPGmFXkXaxNtMT5+5c3RnA49vo/JsrGv3FNXhoccMMxS5cupU6dOqSnp9OtWzdGjhzJ1KlTmThxIkFBQfj7+3PPPfeUeJupqans3r2bH374gccff5zw8PCyvg2lVBmJCJ+HRPPa5iNk5eQy4762jOvtV+H9byzJEieeegOPAGHGmIP5y14EfABEZDGwGbgPiALSgHEW2K7VLFy4kLVr1wIQHR3N8ePHeeKJJ/jyyy9ZvHgxBw8eLBgbFBREZmYmKSkpXLp0iU6dOgGwYMGCgtbIo0ePBvIasyUlJXHlyhVq1apVwe9KKXXNmYRUpq8JY/eJBHo0r8P8ER3wq1fD2mWVmSVm6fxI0efoC48R4G9l3db1bnYkXh527NjB9u3b2bNnD25ubvTr14+MjAzS0tKIiYkBICUlBQ8PDwB+/vnnguctX76c5cuX/+E1r//ItS18BFspR5STKyz76RRvbIvEuUoV/jm8PaO6NbVaszNLs49LyzYkMTGR2rVr4+bmxtGjR9m7dy8AU6dOZcyYMfj6+jJ+/Hg2bdpU4tf8/PPP6d+/Pz/++COenp54enqWV/lKqWJE/pbMlOBQDkVf4a42DZg7vB1enuXf7KwiaeCX0sCBA1m8eDEdOnSgdevW9OjRg507dxISEsJPP/2Ek5MTwcHBLFu2jHHjSnbmqnbt2vTq1YukpCSWLl1azu9AKVXY1exc/rsjine/j8LD1Zl3RnViSEfbaHZmadoe2U446vtWqjwdir7ClNWhRJ5PZminxsx+wJ+6Ntb/prS0PbJSShWSfjWHN7+J5MMfT9HAw5UPHw3krrYNrV1WudPAV0o5lN0nLjItOIyzl9L4c5AP0wa1oaars7XLqhB2GfgiUinPrxXHlk+7KWUvkjKymLf5KJ/tO4tvXTc+G9+Dni3qWrusCmV3ge/q6kpCQgJ169Z1iNAXERISEnC9rt+2Uqrkth8+z4x1YVxIzmTC7c2ZfHcrqrvY3ydly8ruAt/b25uYmBgcobHaNa6urnh7e1u7DKXsTkJKJi9vPMyGQ7G0aeTBkkcC6djUcT/UaHeB7+zsTLNmzaxdhlLKhokIGw7FMmdDBCmZ2Tw/oBUT72iBS1X7bYtgCXYX+EopdSNxienMXBvOt0fj6dS0Fq8/2IFWDT2sXZZN0MBXSlUKubnCZyFnmbf5KDm5wqwH/Hmslx9OlaQtgiVo4Cul7N6pi6lMCw7l51OX6H1bXeYN74BPXTdrl2VzNPCVUnYrOyeXpT+d4t/bjuFStQoLRrbnT4FNHWIG363QwFdK2aUjcUlMDQ4lNCaRAf4NmTusHQ1r6vTlG9HAV0rZlczsHN79Lor/7jiBZ3VnFv25M/e399Kj+hLQwFdK2Y1fzl5m6upQjsenMKJzE2Y94E/tGi7WLstuaOArpWxe2tVs3th6jGW7T+FV05Vl47rRv3UDa5dldzTwlVI27aeoi0xbE0r0pXQe6eHLlIGt8XCQZmeWpoGvlLJJielZ/POrI3y+P5pm9Wrw+YQeBDV3rGZnlqaBr5SyOdsifmPmunASUq8y8Y4WTLq7Ja7OjtfszNIsEvjGmKXAA0C8iLQrYn0/YD1wKn/RGhF5xRLbVkpVHheSM5mzMYKvQuNo61WTDx/tRntv/Y5nS7HUEf5yYBGw4gZjdonIAxbanlKqEhER1v56jlc2HSYtM4cX7m3NhNub4+zk2M3OLM0igS8iPxhj/CzxWkopx3LuSjoz1oaxI/ICXXzymp3d1kCbnZWHivz12dMYc8gYs8UYE1DcIGPMBGPMfmPM/vLueR8XF8eQIUNo3DjvG+pPnz59w/GnT5+mf//+uLm50aZNG7Zv3/679Z9++im+vr7UqFGDYcOGcenSpXKsXin7lpsrfLznNPe8uZN9py4xZ7A/X07spWFfjioq8H8BfEWkI/AfYF1xA0VkiYgEikhg/fr1y7WoKlWqMHDgQIKDg0s0fvTo0XTu3JmEhARee+01HnzwwYIvYomIiODJJ5/k448/5vz587i5ufH000+XZ/lK2a2TF1IYtWQvs9ZH0MW3Nlsn3c5jvZtpZ8vyJiIWuQF+QHgJx54G6t1sXNeuXeVG5s2bJ82bNxd3d3dp27atrFmz5obji5OVlSWAnDp1qtgxkZGR4uLiIklJSQXL+vTpI++9956IiEyfPl1Gjx5dsC4qKkqcnZ1/N14pR5eVnSP//T5KWs7YLO1f+lq+CDkrubm51i6rUgH2SzGZWiHTMo0xjYDzIiLGmO7k/WWRUNbXbdGiBbt27aJRo0Z8+eWX/OUvfyEqKooTJ07wwAPFXx/etGkTffr0KdW2IiIiaN68OR4e//fnZseOHYmIiChY36tXr9/V5uLiwrFjx+jatWsp35lSlU9EbCJTg0MJP5fEwIBGvDI0gAba7KxCWWpa5mdAP6CeMSYGeAlwBhCRxcCDwFPGmGwgHRiV/5uoTB566KGC+w8//DDz5s1j3759DB06lCtXrpT15X8nJSUFT8/fTw/z9PTk3LlzN1yfnJxs0TqUsjcZWTn857vjLN55ktpuLrw3pguD2ntZuyyHZKlZOqNvsn4RedM2LWrFihW8+eabBRdbU1JSuHjxoqU3A4C7uztJSUm/W5aUlFRwxH+z9Uo5ogNnLjFldSgnLqQysos3sx5oSy03bXZmLXY7yfXMmTOMHz+eRYsWkZCQwJUrV2jXrh0iwq5du3B3dy/2tmvXrlJvLyAggJMnT/7uiP3QoUMEBAQUrD906FDBupMnT5KZmUmrVq3K/maVsjOpmdnM2RDBg4v3kJGVy0ePd+fff+qoYW9ldttaITU1FWMM12byLFu2jPDwcAD69u1LSkpKiV4nIyODnJwcADIzM8nIyMDV9Y/nFVu1akWnTp14+eWXmTt3Llu2bCE0NLRghs+YMWPo2bMnu3btokuXLsyePZsRI0boEb5yOD8cu8D0NWHEJqYztocvLwxsg3s1u42aSsVu/yv4+/vzj3/8g549e1KlShXGjh1L7969S/061atXL7jfpk0bgGsziZg4cSIAixcvBmDVqlU89thj1K5dGx8fH1avXl3wCycgIIDFixczZswYEhISuPvuu1m2bFmZ3qNS9uRK2lXmfnWE1QdiaF6/Bl8+2ZNAvzrWLksVYixw7bTcBAYGyv79+61dhlLqJraExTFrfQSX064y8Y7m/P1ObXZmLcaYAyISWNQ6uz3CV0pZX3xyBi+tj2BL+G8ENK7JR493I6CxNjuzVRr4SqlSExFWH4hh7ldHSM/KYcrA1ozvq83ObJ0GvlKqVKIvpfHi2jB2Hb9IN7/azB/ZgRb13a1dlioBDXylVInk5gor9pzm9a2RGOCVoQH8JciXKtr/xm7Y7d9ffn5+f+hWqZQqH1HxKfzp/T3M2XiYbn512Dr5dsb29NOwtzN2G/jl4dKlSwwfPpwaNWrg6+vLp59+WuzYzMxMJk+eTOPGjalduzZPP/00WVlZBev79euHq6trwYe9WrduXRFvQSmLysrJ5d3vo7jvnV1EXUjhzT91ZPm4bnjXdrN2aeoWaOAX8re//Q0XFxfOnz/PypUreeqppwqao11v/vz57N+/n/DwcI4dO8Yvv/zC3Llzfzdm0aJFpKSkkJKSQmRkZEW8BaUsJvxcIkMX/cS/tkYywL8h30y+gxFdvDFGj+rtlV0HfkhICP7+/tSuXZtx48aRkZFxy6+VmppKcHAwr776Ku7u7vTp04chQ4bw8ccfFzl+48aNPPvss9SpU4f69evz7LPPsnTp0lvevlK2IiMrhwVfH2Xouz9xISWTxX/pyrtjulDfo5q1S1NlZNeBv3LlSrZu3cqJEyc4duwYc+fO5ezZs9SqVavYW3GnaY4dO4aTk9Pvet8Ubn98Pfm/3v4Fj2NiYkhMTCxYNn36dOrVq0fv3r3ZsWOHZd60UuUo5PQl7ntnF+/tOMHILk3YPvkOBrZrZO2ylIXY9SydZ555hqZNmwIwY8YM/v73vzN37txbao1c2vbGgwYN4p133qF///7k5OSwcOFCANLS0vD09GTBggX4+/vj4uLCqlWrGDx4MAcPHqRFixalrk2p8paSmc3rXx9lxZ4zeNeuzid/DaJPy3rWLktZmF0f4V8LewBfX19iY2NL/NxBgwYVXFBduXJlqdsbz5gxg86dO9OpUyd69erFsGHDcHZ2pkGDBgAEBQXh4eFBtWrVePTRR+nduzebN2++hXepVPn6PjKee97cycd7z/B472ZsnXS7hn0lZdeBHx0dXXD/7NmzNG7cmLNnz96wNfLKlSsB2LJlS8EF1TFjxtCqVSuys7M5fvx4wWsWbn98verVq7No0SLOnTvHyZMnqVu3Ll27dsXJqej+IcYYbLlvkXI8l1Ov8vznBxm3LAS3alVZPbEXswf7U0M7W1ZexX33oS3cbvSdtr6+vtKuXTuJjo6WhIQE6dOnj0yfPr10X/54nYcfflhGjRolKSkp8uOPP0rNmjUlPDy8yLExMTFy7tw5yc3NlT179oi3t7ds3bpVREQuX74sX3/9taSnp0tWVpZ88skn4ubmJkePHi1TfUpZQm5urmw6FCtdX90mLaZ/Jf/eelQysrKtXZayEG7wnbZWD/Ub3W4W+P/85z+lbdu24unpKWPHjpXU1NQy7aiEhAQZOnSouLm5SdOmTWXlypUF686cOSM1atSQM2fOiIjIzp07xdfXV6pXry6tWrWSTz75pGBsfHy8BAYGiru7u3h6ekpQUJBs27atTLUpZQnnE9Nl/Ech4jt1kzywcJdEnEu0dknKwm4U+NoeWSkHICJ8uT+GV786zNXsXJ4f0Iq/9mlGVW12Vuloe2SlHNjZhLxmZz9GXaR7szrMH9Ge5trszCFZJPCNMUuBB4B4EWlXxHoDvAPcB6QBj4nIL5bYtlKqaDm5wvLdp3ljayROVQxzh7Xjz919tP+NA7PUEf5yYBGwopj1g4CW+bcg4L38f5VS5eD4+WSmBIfy69kr9G9dn9eGt6dxreo3f6Kq1CwS+CLygzHG7wZDhgIr8i8o7DXG1DLGeIlInCW2r5TKczU7l8U7T7DouyhqVHPi7Yc7MbRTY+1/o4CKO4ffBIgu9Dgmf9kfAt8YMwGYAODj41MhxSlVGYTGXGHK6lCO/pbM4I6NeWmwP/Xctf+N+j8VFfhFHV4UOT1IRJYASyBvlk55FqVUZZB+NYe3tx/jg10nqe9RjQ/GBjLAv6G1y1I2qKICPwZoWuixN1DyPghKqSLtPZnAtOBQTiekMbp7U6YNaotndWdrl6VsVEVNwt0AjDV5egCJev5eqVuXnJHFjLVhjFqyl1yBT58IYt6IDkWGfeLGjRy/8y6OtPXn+J13kbhxoxUqVrbAUtMyPwP6AfWMMTHAS4AzgIgsBjaTNyUzirxpmeMssV2lHNF3R88zY20455MyeKJPM/5xT2uquxTdwylx40biZs1G8r8rIjs2lrhZswHwHDy4wmpWtkE/aauUnbiUepVXNkaw7mAsrRq6s2BkBzr71L7hc47feRfZRXSRrdq4MS2/+7a8SlVWpJ+0VcqOiQgbQ+OYsyGC5IwsnrurJX/rfxsuVW9+RjY7rugzp8UtV5WbBr5SNuy3xAxmrgtn+5HzdPT2ZMGDQbRpVLPEz6/q5VX0Eb6XlyXLVHZCOycpZYNEhM/2nWXAmzv5MeoCM+9vy5qne5cq7AEaTJ6EcXX93TLj6kqDyZMsWa6yE3qEr5SNOZOQyrTgMPacTKBn87rMH9ke37o1bum1rl2YjX/rbbLj4qjq5UWDyZP0gq2D0sBXykbk5ArLfjrFG9sica5ShXkj2jOqW9Myt0XwHDxYA14BGvhK2YTI3/KanR2KvsLdbRswd1h7Gnm63vyJSpWCBr5SVnQ1O5d3v4/ivzui8HB1ZuHozgzu4KXNzlS50MBXykoORl9hyupDHDufwtBOjXlpcAB1arhYuyxViWngK1XB0q/m8O9tkSz96RQNPFz58NFA7mqrzc5U+dPAV6oC7T5xkWnBYZy9lMaYIB+mDWqDh6s2O1MVQwNfqQqQlJHFvM1H+GxfNH513Vg1oQc9mte1dlnKwWjgK1XOth8+z4x1YVxIzuTJ25sz6e5WxTY7U6o8aeArVU4upmTy8sbDbDwUS5tGHnwwNpAO3rWsXZZyYBr4SlmYiLD+YCwvb4wgJTOb5we0YuIdLUrU7Eyp8qSBr5QFxV5JZ+a6cL47Gk9nn1osGNmBVg09rF2WUoAGvlIWkZsrfLrvLPO3HCUnV5j9gD+P9vLDqYp+gErZDg18pcro1MVUpgWH8vOpS/S+rS7zhnfAp66btctS6g808JW6Rdk5uXz44yne/OYYLlWr8PrIDjwU6K1tEZTN0sBX6hYcjk1ianAoYecSGeDfkLnD2tGwpjY7U7ZNA1+pUsjMzmHRd1G8t+MEtdyceffPXbivfSM9qld2wSKBb4wZCLwDOAH/E5H5163vB6wHTuUvWiMir1hi20pVlANnLjM1OJSo+BRGdGnCrPv9qa3NzpQdKXPgG2OcgHeBAUAMEGKM2SAih68buktEHijr9pSqaGlXs/nX1kiW7z6NV01Xlo3rRv/WDaxdllKlZokj/O5AlIicBDDGrAKGAtcHvlJ258fjF5m2JpSYy+mM7enLlIFtcK+mZ0KVfbLET24TILrQ4xggqIhxPY0xh4BY4P+JSERRL2aMmQBMAPDx8bFAeUqVXmJaFq9tPswX+2NoVq8GXzzZk+7N6li7LKXKxBKBX9TVKrnu8S+Ar4ikGGPuA9YBLYt6MRFZAiwBCAwMvP51lCp3X4f/xqz14VxKvcpT/Vrw3F0tcXXWZmfK/lki8GOApoUee5N3FF9ARJIK3d9sjPmvMaaeiFy0wPaVsogLyZnM2RDBV2Fx+HvVZNlj3WjXxNPaZSllMZYI/BCgpTGmGXAOGAX8ufAAY0wj4LyIiDGmO1AFSLDAtpUqMxFhzS/neGXTYdKv5vDCva2ZcHtznJ202ZmqXMoc+CKSbYx5BthK3rTMpSISYYyZmL9+MfAg8JQxJhtIB0aJiJ6uUVZ37ko6L64JY+exC3T1rc2CkR24rYG7tctSqlwYW87dwMBA2b9/v7XLUJVQbq7wyc9nWLDlKAJMubc1Y3v6UUWbnSk7Z4w5ICKBRa3T+WXK4Zy4kMK04FBCTl+mb8t6/HN4e5rW0WZnqvLTwFcOIysnlw92neTt7cdxrVqFfz3YgQe7arMz5Tg08JVDCD+XyNTgUCJikxjUrhEvDw2ggYc2O1OORQNfVWoZWTn857vjLN55ktpuLrw3pguD2ntZuyylrEIDX1Va+09fYkpwKCcvpPJgV29m3t+WWm7a7Ew5Lg18VemkZuY1O/toz2kae1ZnxePdub1VfWuXpZTVaeCrSmXnsQu8uCaM2MR0Hu3pxwv3tqaGNjtTCtDAV5XElbSrvLrpCMG/xNCifg2+fLIngX7a7EypwjTwldUkbtxI/Ftvkx0XR1UvLxpMnoTn4MGlfp0tYXHMWh/B5bSrPNP/Np658zZtdqZUETTwlVUkbtxI3KzZSEYGANmxscTNmg1Q4tCPT8pg9voIvo74jYDGNfno8W4ENNZmZ0oVRwNfWUX8W28XhP01kpFB/Ftv3zTwRYTVB2J4ddNhMrJzmTqwDeP7NqOqNjtT6oY08JVVZMfFlWr5NdGX0nhxbRi7jl+km19t5o/sQIv62uxMqZLQwFdWUdXLi+zY2CKXFyUnV1ix5zT/2hqJAV4dGsCYIF9tdqZUKejfwMoqGkyehHH9fWsD4+pKg8mT/jA2Kj6ZP72/h5c3HqabXx22PX8Hj2hnS6VKTY/wlVVcO09/o1k6WTm5vL/zBAu/jcKtmhNv/qkjwzs30WZnSt0iDXxlNZ6DBxd7gTb8XCIvrA7lSFwS93fwYs7gAOp7VKvgCpWqXDTwlU3JyMrh7e3H+WDXSerUcOH9R7pyb0Aja5elVKWgga9sxs8nE5i2JoxTF1N5OLApL97XFk83Z2uXpVSloYGvrC45I4vXv47k471n8K5dnU/+GkSflvWsXZZSlY4GvrKq7yPjmbEmjLikDB7v3Yz/d28r3Fz0x1Kp8mCRaZnGmIHGmEhjTJQxZloR640xZmH++lBjTBdLbFfZr8upV3n+84OMWxZCjWpVCX6qF7MH+2vYK1WOyvx/lzHGCXgXGADEACHGmA0icrjQsEFAy/xbEPBe/r/KwYgIX4XF8dL6CBLTs3j2ztv42523Ua2qNjtTqrxZ4nCqOxAlIicBjDGrgKFA4cAfCqwQEQH2GmNqGWO8ROTGn6NXlcr5pAxmrgvnm8Pnad/Ek0+eCKKtV01rl6WUw7BE4DcBogs9juGPR+9FjWkC/CHwjTETgAkAPj4+FihPWZuI8MX+aOZ+dYSr2blMH9SGv/bRZmdKVTRLBH5RH3uUWxiTt1BkCbAEIDAwsMgxyn6cTUhj2ppQdp9IoHuzOiwY2YFm9WpYuyylHJIlAj8GaFrosTdwfVeskoxRlUhOrrB892ne2BqJUxXDa8PbMbqbj/a/UcqKLBH4IUBLY0wz4BwwCvjzdWM2AM/kn98PAhL1/H3ldex8MlNWh3Iw+gp3tmnAa8Pb4eVZ3dplKeXwyhz4IpJtjHkG2Ao4AUtFJMIYMzF//WJgM3AfEAWkAePKul1le65m57J45wn+891x3KtV5Z1RnRjSsbE2O1PKRlhk0rOIbCYv1AsvW1zovgB/s/hHhk8AAA4iSURBVMS2lG06FH2FqcGhHP0tmcEdGzNnsD913bXZmVK2RD/losok/WoOb20/xv92naS+RzU+GBvIAP+G1i5LKVUEDXx1y/acSGD6mlBOJ6QxursP0+9rQ01XbXamlK3SwFellpSRxfwtR/n057P41nXj0/FB9Gqhzc6UsnUa+KpUvjt6nhfXhBOfnMH4vs14fkBrqrtoWwSl7IEGviqRhJRMXtl0mPUHY2nd0IPFj3SlU9Na1i5LKVUKGvjqhkSEDYdieXnjYZIzsph0d0ue7ncbLlW1LYJS9kYDXxUrLjGdmWvD+fZoPB2b1uL1kR1o3cjD2mUppW6RBr76g9xcYVVINPM2HyErN5eZ97dlXO9mOGlbBKXsmga++p3TF1OZtiaUvScv0bN5XeaPbI9vXW12plRloIGvgLxmZ0t/PMW/v4nEuUoV5o9oz8PdmmpbBKUqEQ18xdHfkpi6OpRDMYnc3bYBc4e1p5Gnq7XLUkpZmAa+A8vMzuHd70/w3++j8KzuzH9Gd+aBDl56VK9UJaWB76B+PXuZqcGhHDufwrBOjZk9OIA6NVysXZZSqhxp4DuYtKvZ/HvbMZb+dIpGNV1Z+lggd7bRZmdKOQINfAeyO+oi09aEcfZSGn/p4cPUgW3w0GZnSjkMDXwHkJiexbzNR1gVEo1fXTdWTehBj+Z1rV2WUqqCaeBXctsifmPmunAupmTy5B3NmXx3K1ydtdmZUo5IA7+SupiSyZwNEWwKjaNNIw/+92ggHby12ZlSjkwDv5IREdYdPMfLGw+TlpnDPwa04sk7WmizM6WUBn5lEnslnRlrw/g+8gKdffKanbVsqM3OlFJ5yhT4xpg6wOeAH3Aa+JOIXC5i3GkgGcgBskUksCzbVb+Xmyus3HeWBVuOkpMrzH7An0d7+WmzM6XU75T1CH8a8K2IzDfGTMt/PLWYsf1F5GIZt6euc/JCCtPWhLHv1CX63FaPeSPa07SOm7XLUkrZoLIG/lCgX/79j4AdFB/4yoKyc3L534+neOubY7hUrcLrIzvwUKC3tkVQShWrrIHfUETiAEQkzhjToJhxAmwzxgjwvogsKe4FjTETgAkAPj4+ZSyvcjocm8SU4EOEn0viHv+GvDqsHQ1rarMzpdSN3TTwjTHbgUZFrJpRiu30FpHY/F8I3xhjjorID0UNzP9lsAQgMDBQSrGNSi8zO4dF30Xx3o4T1HJz5r9jujCoXSM9qldKlchNA19E7i5unTHmvDHGK//o3guIL+Y1YvP/jTfGrAW6A0UGviragTN5zc6i4lMY0aUJs+73p7Y2O1NKlUJZT+lsAB4F5uf/u/76AcaYGkAVEUnOv38P8EoZt+swUjOzeWNbJMt3n6axZ3WWj+tGv9bFnTlTSqnilTXw5wNfGGP+CpwFHgIwxjQG/ici9wENgbX5px2qAp+KyNdl3K5D2HX8AtPXhBFzOZ2xPX2ZMrAN7tX0oxNKqVtTpvQQkQTgriKWxwL35d8/CXQsy3YcTWJaFnO/OsyXB2JoXq8GXzzZk+7N6li7LKWUndPDRRvzdfhvzFofzqXUqzzVrwXP3dVSm50ppSxCA99GxCdnMGdDBJvDfsPfqybLHutGuyae1i5LKVWJaOBbmYiw5pdzvLLpMOlZObxwb2sm3N4cZydtdqaUsiwNfCuKuZzGi2vD+eHYBbr61mbByA7c1sDd2mUppSopDXwryM0VPt57hgVfHwXg5SEBPNLDlyra7EwpVY408CvYiQspTF0dyv4zl+nbsh7/HK7NzpRSFUMDv4Jk5eSy5IeTvPPtcao7O/HGQx0Z2aWJtkVQSlUYDfwKEH4ukanBoUTEJnFf+0bMGRJAAw9tdqaUqlga+OUoIyuHhd8e5/0fTlLbzYXFf+nCwHZe1i5LKeWgNPDLScjpS0wNDuXkhVQe6urNzPv98XRztnZZSikHpoFvYSmZ2bz+9VFW7DlDk1rVWfF4d25vVd/aZSmllAa+Je08doEX14QRm5jOY738eOHe1tTQZmdKKRuhaWQBV9Ku8sqmw6z55Rwt6tfgyyd7Euinzc6UUrZFA7+MNofFMXt9OFfSsnim/208c+dt2uxMKWWTNPBvUXxSBrPWh7M14jztmtTko8e7E9BYm50ppWyXBn4piQhfHohh7qbDZGTnMnVgG8b3bUZVbXamlLJxGvilEH0pjelrwvgx6iLd/eowf2R7mtfXZmdKKfuggV8CObnCij2nef3rSKoYeHVoAGOCtNmZUsq+aODfRFR8MlNWh/LL2Sv0a12f14a3p0mt6tYuSymlSk0DvxhZObm8v/MEC7+Nwq2aE2893JFhnbTZmVLKfpUp8I0xDwFzgLZAdxHZX8y4gcA7gBPwPxGZX5btlrewmEQ+m7eEe3avYV36FZwaedGo82RMZ29rl6aUUresrEf44cAI4P3iBhhjnIB3gQFADBBijNkgIofLuG2Ly8jK4a3txzi28kue+3U11XKyAMj9LY64WbMB8Bw82JolKqXULSvTXEIROSIikTcZ1h2IEpGTInIVWAUMLct2y8PPJxMY9M4u3t95kqePf1MQ9tdIRgbxb71tpeqUUqrsKmLyeBMgutDjmPxlRTLGTDDG7DfG7L9w4UK5F5eckcXMdWE8vGQv2bm5rHwiCPekhCLHZsfFlXs9SilVXm56SscYsx1oVMSqGSKyvgTbKOoqpxQ3WESWAEsAAgMDix1nCd8fjWfG2jDikjL4a59m/OOeVri5VOW4lxfZsbF/GF/VS3vZK6Xs100DX0TuLuM2YoCmhR57A39M0wp0KfUqr246zNpfz9GygTvBT/Wii0/tgvUNJk8ibtZsJCOjYJlxdaXB5EnWKFcppSyiIqZlhgAtjTHNgHPAKODPFbDdPxARNoXGMWdDBInpWTx7V0v+1r8F1ar+vtnZtQuz8W+9TXZcHFW9vGgweZJesFVK2bWyTsscDvwHqA98ZYw5KCL3GmMakzf98j4RyTbGPANsJW9a5lIRiShz5aV0PimDGWvD2X7kPB28PfnkiSDaetUsdrzn4MEa8EqpSsWIlOtp8jIJDAyU/fuLnNpfYiLC5yHRvLb5CFezc/nHPa14vLc2O1NKVU7GmAMiEljUukr9SduzCWlMWxPK7hMJBDWrw4KRHfCrV8PaZSmllFVUysDPyRWW/XSKN7ZFUrVKFV4b3o7R3Xy02ZlSyqFVusBPTMvi0WX7OBh9hTvbNOC14e3w8tRmZ0opVekCv2b1qvjWdWNcbz+GdGyszc6UUipfpQt8YwzvjOps7TKUUsrm6FQVpZRyEBr4SinlIDTwlVLKQWjgK6WUg9DAV0opB6GBr5RSDkIDXymlHIQGvlJKOQib7pZpjLkAnLFyGfWAi1auoTTsrV7QmiuKvdVsb/WCbdTsKyL1i1ph04FvC4wx+4trNWqL7K1e0Jorir3VbG/1gu3XrKd0lFLKQWjgK6WUg9DAv7kl1i6glOytXtCaK4q91Wxv9YKN16zn8JVSykHoEb5SSjkIDXyllHIQGviAMWagMSbSGBNljJlWxHpjjFmYvz7UGNPFGnVeV9PNau5njEk0xhzMv822Rp2F6llqjIk3xoQXs94W9/HNara1fdzUGPO9MeaIMSbCGPNcEWNsaj+XsGZb28+uxph9xphD+TW/XMQYm9rPBUTEoW+AE3ACaA64AIcA/+vG3AdsAQzQA/jZDmruB2yy9v4tVM/tQBcgvJj1NrWPS1izre1jL6BL/n0P4Jgd/CyXpGZb288GcM+/7wz8DPSw5f187aZH+NAdiBKRkyJyFVgFDL1uzFBgheTZC9QyxnhVdKGFlKRmmyIiPwCXbjDE1vZxSWq2KSISJyK/5N9PBo4ATa4bZlP7uYQ125T8fZeS/9A5/3b97Beb2s/XaODn/XBFF3ocwx9/4EoypiKVtJ6e+X92bjHGBFRMabfM1vZxSdnkPjbG+AGdyTv6LMxm9/MNagYb28/GGCdjzEEgHvhGROxiP1e6LzG/BaaIZdf/ti7JmIpUknp+Ia+nRoox5j5gHdCy3Cu7dba2j0vCJvexMcYdCAYmiUjS9auLeIrV9/NNara5/SwiOUAnY0wtYK0xpp2IFL7WY5P7WY/w837zNi302BuIvYUxFemm9YhI0rU/O0VkM+BsjKlXcSWWmq3t45uyxX1sjHEmLzhXisiaIobY3H6+Wc22uJ+vEZErwA5g4HWrbG4/gwY+QAjQ0hjTzBjjAowCNlw3ZgMwNv/Kew8gUUTiKrrQQm5aszGmkTHG5N/vTt5/64QKr7TkbG0f35St7eP8Wj4EjojIm8UMs6n9XJKabXA/188/sscYUx24Gzh63TCb2s/XOPwpHRHJNsY8A2wlb/bLUhGJMMZMzF+/GNhM3lX3KCANGGetevNrKknNDwJPGWOygXRglORPH7AGY8xn5M22qGeMiQFeIu9il03uYyhRzTa1j4HewCNAWP75ZYAXAR+w2f1ckpptbT97AR8ZY5zI++XzhYhssuXMuEZbKyillIPQUzpKKeUgNPCVUspBaOArpZSD0MBXSikHoYGvlFIOQgNfKaUchAa+Uko5CA18pUrIGNMtv7e5qzGmRn4v9HbWrkupktIPXilVCsaYuYArUB2IEZF5Vi5JqRLTwFeqFPJ7F4UAGUCv/K6JStkFPaWjVOnUAdzJ+3YmVyvXolSp6BG+UqVgjNlA3jeMNQO8ROQZK5ekVIk5fLdMpUrKGDMWyBaRT/M7Je42xtwpIt9ZuzalSkKP8JVSykHoOXyllHIQGvhKKeUgNPCVUspBaOArpZSD0MBXSikHoYGvlFIOQgNfKaUcxP8HP+t1b2alfikAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xfit = np.array([-0.2, 3.2])\n", "yfit = fit[0]*xfit + fit[1]\n", "plt.plot(x, y, 'oC3', label=\"data\")\n", "plt.plot(xfit, yfit, zorder=-1, label=\"ax+b\")\n", "plt.text(-0.3, 1.1, \"a={0:0.2f}\\nb={1:0.2f}\"\n", " .format(fit[0], fit[1]), fontsize=12)\n", "plt.legend(loc=\"upper left\")\n", "plt.xlabel('x')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perform same fit without weighting but get estimates for uncertainties in slope and y-intercept from covariance matrix." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slope = 1.00 ± 0.07\n", "y-intercept = -0.95 ± 0.13\n" ] } ], "source": [ "fit, cvm = linfit(x, y, relsigma=True)\n", "dfit = [np.sqrt(cvm[i,i]) for i in range(2)]\n", "print(u\"slope = {0:0.2f} \\xb1 {1:0.2f}\".format(fit[0], dfit[0]))\n", "print(u\"y-intercept = {0:0.2f} \\xb1 {1:0.2f}\".format(fit[1], dfit[1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Demonstration of a fit to data with error estimates for each data point" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slope = 21.7 ± 1.2\n", "y-intercept = -72 ± 15\n" ] } ], "source": [ "# data set for linear fitting \n", "x = np.array([2.3, 4.7, 7.1, 9.6, 11.7, 14.1, 16.4, 18.8, 21.1, 23.0])\n", "y = np.array([-25., 3., 110., 110., 230., 300., 270., 320., 450., 400.])\n", "sigmay = np.array([15., 30., 30., 40., 40., 50., 40., 30., 50., 30.])\n", "\n", "# Fit linear data set with weighting\n", "fit, cvm, info = linfit(x, y, sigmay=sigmay, relsigma=False, return_all=True)\n", "dfit = [np.sqrt(cvm[i,i]) for i in range(2)]\n", "print(u\"slope = {0:0.1f} \\xb1 {1:0.1f}\".format(fit[0], dfit[0]))\n", "print(u\"y-intercept = {0:0.0f} \\xb1 {1:0.0f}\".format(fit[1], dfit[1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the data with error bars together with the fit. Plot the residuals in a separate graph above the data with fit." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAHgCAYAAADT1NXlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeVzVVf7H8dcBV9xNLRUFzRUXEFzLErPUUbN0cikzbdPRppnpVxZNZBuWk61W1tiiprjlTGW55u6UqahYCiKmuCeaoii4AOf3BxfCREG5lwvc9/PxuI97v+d7vt/vh3uB7+ee7/meY6y1iIiIiGfxcncAIiIiUviUAIiIiHggJQAiIiIeSAmAiIiIB1ICICIi4oGUAIiIiHigUu4OwNVq1Khh/f393R2GiIhIodi0adMxa23NvOqV+ATA39+fqKgod4chIiJSKIwxe/NTT5cAREREPJASABEREQ+kBEBERMQDKQEQERHxQEoAREREPJASABEREQ+kBEBERMQDKQEQERHxQEoAREREPJASABEREQ+kBEBERMQDKQEQERHxQEoAREREPJASABEREQ+kBEBEir3Q0FBCQ0PdHYZIsaIEQEREPIISxYspARAREfFAbk0AjDGfGWMSjTHbcpRVN8Z8Z4yJdzxXy7HuWWPMLmNMnDGmh3uiFhERKf7c3QIwFej5h7IwYLm1tjGw3LGMMSYAGAy0cGwzyRjjXXihioiIlBxuTQCstWuA438ovguY5ng9Dbg7R/lsa+05a+0eYBfQvlACFRERKWHc3QKQm+uttYcBHM+1HOV1gf056h1wlImIiMhVKooJwOWYXMpsrhWNGWGMiTLGRB09etTFYYmIiBQ/RTEBOGKMqQ3geE50lB8A6uWo5wscym0H1trJ1tq21tq2NWvWdGmwIiIixVFRTADmA8Mcr4cBX+coH2yMKWuMaQA0Bja4IT4REZFir5Q7D26MmQWEAjWMMQeAF4DxwFxjzMPAPmAAgLV2uzFmLhADpAGPWWvT3RK4iIhIMefWBMBae+9lVnW7TP1xwDjXRSQiIuIZiuIlABEREXExJQAiAmicdBFPowRARETEAykBEBER8UBKAERERDyQEgAREREPpARARETEAykBEBER8UBKAESkWJsZGcnPGzawZvVqAho0YGZkpLtDEikWlACISLE1MzKSsFGjebNGDaKbNOXp9HTCRo1WEiCXUKJ4qXwlAMaYm40xFRyv7zfGvGWM8XNtaCIiVxYRHs7LVavQwacCpY2hg08FXq5ahYjwcHeHJkWIEsXc5bcF4EMgxRgTCDwN7AU+d1lUIiL5ELd3L8HlfS4qCy7vQ9zevW6KSIoiJYq5y28CkGattcBdwLvW2neBSq4LS0Qkb039/NicmnJR2ebUFJr6qYFSfqdEMXf5TQCSjTHPAvcDC4wx3kBp14UlIpK38IgIxiadZH3KGS5Yy/qUM4xNOkl4RIS7Q5MiRIli7vKbAAwCzgEPW2t/BeoCE1wWlYhIPtw3ZAjjP5zEk8eOEbQzjte9vRn/4STuGzLE3aFJEaJEMXel8lPJcdJ/K8fyPtQHQApJ1gx1q1atcmscUjTdN2QIkz/+GNDviOQuKyF8/NFHOZGaSjN/fyWK5JEAGGOSAZvbKsBaayu7JCoRKVRZt0idSE0loEEDwiMiPP6fo5QsShQvdcUEwFqrjn4iJVzOW6SCy/uwOTWFsFGjAZQEiJRgVzUQkDGmljGmftbDVUGJSOHRLVIinim/AwH1NcbEA3uA1UACsMiFcYlIIdEtUiKeKb8tAK8AHYGd1toGQDfge5dFJSKFRrdIiXim/CYAF6y1vwFexhgva+1KIMiFcYlIIdEtUiKeKV+3AQJJxpiKwBog0hiTCKS5LiwRKSy6RUrEM+U3AbgLOAs8AQwBqgAvuyooESlcukVKxPPkdyCgMzkWp7koFhERESkk+UoA/jAgUBky5wE4o4GAREREiqd8dQK01lay1lZ2PMoBfwbed21oIr+PULdm9WoCGjTw+Pm7RUScJb99AC5irf3KGBPm7GBEctIIdSIirpPfSwD9cyx6AW3JfY4AEafJOUIdkDlCnaNcCYCISMHktwXgzhyv08gcCfAup0cjkkPc3r0EN25yUVlweR/i4ne6KSIRkZIjv3cBPOjqQET+KGuEuqwWANAIdVJyaJprcbe8pgN+jys09Vtr/+b0iEQcwiMiCBs1mpchuw/A2KSTjP9wkrtDExEp9vJqAYhyPN8MBABzHMsDgE2uCkoENEKdiIgrXTEBsNZOAzDGDAe6WmsvOJY/Apa6PDrxeBqhTkTENfLbCbAOUAk47liu6CgTEXE7JYciVy+/CcB4YIsxZqVjuQvwoksiEimB1OFLSpLi+vtc3OJ1tfzeBTDFGLMI6OAoCrPW/uq6sIqm4vpLLyIi8kdXHArYGNPM8RxMZpP/fsejjqNMREREiqG8WgD+DxgBvJnLOgvc5vSI8mCM6Qm8C3gDn1hrxxd2DCIlkVq2RDxLXncBjHA8dy2ccK7MGOMNfADcARwANhpj5ltrY9wbmYiISPGSr9kAjTEDjDGVHK/DjTH/Nca0cW1ouWoP7LLW7rbWngdmU0hDEhf3WelCQ0Oz+zCIiEjR4M7/zflKAIDnrbXJxpjOQA9gGvCR68K6rLpk9kHIcsBR5lI5Z6WLbtKUp9PTCRs1utglASIiIlnyextguuO5N/ChtfZrY8yLrgnpikwuZZcMVWyMGUFm3wVq165d4Gub4WPG5DorXfiYMdSp6/L8wymSkpKA4nmdtzjHnqUk/AziXMX5d6I4x17UuPO9zG8CcNAY82/gduBfxpiy5L/1wJkOAPVyLPsCh/5YyVo7GZgM0LZtW1vQ5pW9v/6a66x0e+N3Fptm9apVqwIUm3hzio6OdncIBVac339xjeL8O1GcYy9q3Ple5vckPhBYAvS01iYB1YExLovq8jYCjY0xDYwxZYDBwHxXHzRrVrqcNCudiIgUZ/lKAKy1KUAi0NlRlAbEuyqoK8SRBvyVzGQkFphrrd3u6uOGR0QwNukk61POcMFa1qecYWzSScIjIlx9aBEpgYp7p2IpGfJ1CcAY8wLQFmgKTAFKAzPInCWwUFlrFwILC/OYmpVORJwlZ6firGmuw0aNBigW/1OykpcTqakENGhAeEREsYhbLpXfPgD9gDbAZgBr7aGs2wI9hWalExFniAgPz7VTcUR4eJE/kRb35EUult8+AOettRZHj3tjTAXXhSQiUnLF7d1LcHmfi8qCy/sQt3evmyLKv5zJS2ljMpOXqlWICA93d2hyDfJMAIwxBvjWcRdAVWPMo8Ay4GNXByciUtIU507FxTl5kUvlmQA4vvnfDcwD/kNmP4Cx1tr3XBxbkbNq1So1/4tIgRTnTsXFOXmRXFhr83yQOf5+u/zULWqPkJAQ6+kiZ8yw1cuXtwZsc39/GzljhrtD8ih6/+WPiuvvROSMGbZepcp2Sr16dmuTpnZKvXq2XqXKxSb+osZVvwdAlM3PuT1flSCGzFv/fgF+ynrkZ1t3Pzw9AdAfrHvp/ZfL6dKli+3SpYu7w7hqxTV5KWpc+b8hvwmAyax7ZcaYXNt3rLVF/sJP27ZtbVRUlLvDcJuABg14Oj09u8cxwPqUM7zu7U3Mnj1ujMwz6P2Xy8ka+a04XlYszrEXFa7832CM2WStbZtXvfwOBLQ3t0eBIpRCoU477qX3X0RyUxT+N7hjPH8pROq04156/0UkN0Xhf4MSgBKuOPc4Lgn0/otIborC/4b8jgQoxZSGMXYvvf8ikpui8L9BCYAH0DDG7qX3X0Ry4+7/DboEICIi4oGUAIiIiHggJQAiIiIeSAmAiIiIB1InQBERuSrqzFoyqAVARETEAykBEBER8UBKAERERDyQ+gCIiLiBrqOLuykB8BD6ZyMiUvS483+zLgGIiIh4ICUAIiIiHkgJgIiIiAdSHwCRQqA+GCJS1KgFQERExAMpARAREfFASgBEREQ8kBIAERERD6QEQERExAMpARAREfFASgBEREQ8kBIAERERD6QEQERExAMpARAREfFASgBEREQ8kFsSAGPMAGPMdmNMhjGm7R/WPWuM2WWMiTPG9MhRHmKM+dmxbqIxxhR+5CIiIiWDu1oAtgH9gTU5C40xAcBgoAXQE5hkjPF2rP4QGAE0djx6Flq0IiIiJYxbEgBrbay1Ni6XVXcBs62156y1e4BdQHtjTG2gsrV2nbXWAp8DdxdiyCIiIiVKUesDUBfYn2P5gKOsruP1H8tFRETkGpRy1Y6NMcuAG3JZ9Zy19uvLbZZLmb1C+eWOPYLMywUAp40xubU2yO9qAMfcHYTocygi9DkUDfocrp1ffiq5LAGw1t5+DZsdAOrlWPYFDjnKfXMpv9yxJwOTr+H4HskYE2WtbZt3TXElfQ5Fgz6HokGfg+sVtUsA84HBxpiyxpgGZHb222CtPQwkG2M6Onr/PwBcrhVBRERE8uCu2wD7GWMOAJ2ABcaYJQDW2u3AXCAGWAw8Zq1Nd2w2CviEzI6BvwCLCj1wERGREsJkdqoXT2aMGeG4bCJupM+haNDnUDToc3A9JQAiIiIeqKj1ARAREZFCoATAwxljEhxDLEcbY6LcHY+nMMZ8ZoxJNMZsy1FW3RjznTEm3vFczZ0xeoLLfA4vGmMOOv4moo0xvdwZY0lnjKlnjFlpjIl1DBH/d0e5/h5cTAmAAHS11gbplptCNZVLh7MOA5ZbaxsDyx3L4lpTyX1Y8bcdfxNB1tqFhRyTp0kDnrTWNgc6Ao85hoXX34OLKQEQcQNr7Rrg+B+K7wKmOV5PQ8Ndu9xlPgcpRNbaw9bazY7XyUAsmSO96u/BxZQAiAWWGmM2OUZQFPe53jHmBY7nWm6Ox5P91Rjzk+MSgZqeC4kxxh9oA6xHfw8upwRAbrbWBgN/IrPp7VZ3ByTiZh8CNwJBwGHgTfeG4xmMMRWB/wD/sNaecnc8nkAJgIez1h5yPCcCXwLt3RuRRzvimPkSx3Oim+PxSNbaI9badGttBvAx+ptwOWNMaTJP/pHW2v86ivX34GJKADyYMaaCMaZS1mugO7DtyluJC80HhjleD0PDXbtF1knHoR/6m3Apx/DunwKx1tq3cqzS34OLaSAgD2aMaUjmt37InBhqprV2nBtD8hjGmFlAKJkznh0BXgC+InMo7PrAPmCAtVYd1FzoMp9DKJnN/xZIAEZmXYsW5zPGdAbWAj8DGY7if5LZD0B/Dy6kBEBERMQD6RKAiIiIB1ICICIi4oGUAIiIiHggJQAiIiIeSAmAiIiIB1ICICIi4oGUAIiIiHggJQAiIiIeSAmAiIiIB1ICICIi4oGUAIiIiHggJQAiIiIeSAmAiIiIB1ICICIi4oGUAIiIiHggJQAiIiIeSAmAiIiIB1ICICIi4oGUAIiIiHggJQAiIiIeSAmAiIiIB1ICICIi4oGUAIiIiHggJQAiIiIeqJS7A3C1GjVqWH9/f3eHISIiUig2bdp0zFpbM696JT4B8Pf3Jyoqyt1hiIiIFApjzN781NMlABEREQ+kBEBERMQDKQEQERHxQEoAREREPJASABEREQ+kBEBERMQDKQEQERHxQG5NAIwxCcaYn40x0caYKEdZdWPMd8aYeMdztRz1nzXG7DLGxBljergvchERkeKtKLQAdLXWBllr2zqWw4Dl1trGwHLHMsaYAGAw0ALoCUwyxni7I2AREZHirigkAH90FzDN8XoacHeO8tnW2nPW2j3ALqC9G+ITEREp9tydAFhgqTFmkzFmhKPsemvtYQDHcy1HeV1gf45tDzjKLmGMGWGMiTLGRB09etRFoYuISGELDQ0lNDTU3WGUCO6eC+Bma+0hY0wt4DtjzI4r1DW5lNncKlprJwOTAdq2bZtrHREREU/m1hYAa+0hx3Mi8CWZTfpHjDG1ARzPiY7qB4B6OTb3BQ4VXrQiIiIlh9sSAGNMBWNMpazXQHdgGzAfGOaoNgz42vF6PjDYGFPWGNMAaAxsKNyoRURESgZ3XgK4HvjSGJMVx0xr7WJjzEZgrjHmYWAfMADAWrvdGDMXiAHSgMestenuCV1ERKR4c1sCYK3dDQTmUv4b0O0y24wDxrk4NBERkRLP3XcBiIiIiBsoARAREfFASgBEREQ8kBIAERERD6QEQERExAMpAXCy0NBQoqKi3B0GACkpKfTu3ZtmzZrRokULwsLCstetWbOG4OBgSpUqxbx583LdPjk5maCgoOxHjRo1+Mc//lHguJ577jnq1atHxYoVL1vnu+++IyQkhFatWhESEsKKFSsKfFwREfmdEoAS7qmnnmLHjh1s2bKF77//nkWLFgFQv359pk6dyn333XfZbStVqkR0dHT2w8/Pj/79+1/xeMOHD2fVqlVXrHPnnXeyYcOVx3CqUaMG33zzDT///DPTpk1j6NChV6wvIiJXRwnANTpz5gy9e/cmMDCQli1bMmfOnEvqzJo1i1atWtGyZUueeeaZ7PKKFSvy5JNPEhwcTLdu3ciasOiXX36hZ8+ehISEcMstt7Bjx5WmRsibj48PXbt2BaBMmTIEBwdz4MABAPz9/WndujVeXvn7FYiPjycxMZFbbrmlQDEBdOzYkdq1a1+xTps2bahTpw4ALVq04OzZs5w7d67AxxYRkUxKAK7R4sWLqVOnDlu3bmXbtm307NnzovWHDh3imWeeYcWKFURHR7Nx40a++uorIDN5CA4OZvPmzXTp0oWXXnoJgBEjRvDee++xadMm3njjDUaPHn3JcVeuXHlRs3zW46abbrpivElJSXzzzTd065brGEt5mjVrFoMGDcIxcuNFlixZkh3H/PnzeeSRRwgKCqJDhw7XdKw/+s9//kObNm0oW7asU/YnIiLunw2w2GrVqhVPPfUUzzzzDH369Lnkm/HGjRsJDQ2lZs2aAAwZMoQ1a9Zw99134+XlxaBBgwC4//776d+/P6dPn+aHH35gwIAB2fvI7Rtv165diY6OvqpY09LSuPfee/nb3/5Gw4YNr/ZHBWD27NlMnz4913U9evSgR48eQOYlgOHDhzttus7t27fzzDPPsHTpUqfsT0REMikBuEZNmjRh06ZNLFy4kGeffZbu3bszduzY7PXW5n8WYmMMGRkZVK1aNc+T+8qVK3niiScuKffx8eGHH37IdZsRI0bQuHHja+7At3XrVtLS0ggJCbmm7a/VgQMH6NevH59//jk33nhjoR5bRKSkUwJwjQ4dOkT16tW5//77qVixIlOnTr1ofYcOHfj73//OsWPHqFatGrNmzeLxxx8HICMjg3nz5jF48GBmzpxJ586dqVy5Mg0aNOCLL75gwIABWGv56aefCAy8eLqEq20BCA8P5+TJk3zyySfX/LPOmjWLe++9N191//g+XKukpCR69+7Na6+9xs033+yUfYqIyO/UB+Aa/fzzz7Rv356goCDGjRtHeHj4Retr167Na6+9RteuXQkMDCQ4OJi77roLgAoVKrB9+/bs29uyWg4iIyP59NNPCQwMpEWLFnz99deXHPdqHDhwgHHjxhETE0NwcDBBQUHZicDGjRvx9fXliy++YOTIkbRo0SJ7u6CgoIv2M3fu3CsmADn7AOR8XK4PwNNPP42vry8pKSn4+vry4osvAjB//vzs9+L9999n165dvPLKK9n7S0xMLMjbISIiOZiraaoujtq2bWuLyn35WSpWrMjp06fdHYaISLGT1b8or9uNPZkxZpO1tm1e9dQCICIi4oGUALiBvv2LiFy9mZGR/LxhA2tWryagQQNmRka6O6RiTZ0ARUSkyJsZGUnYqNG8WaMGweV92JyaQtiozLFS7hsyxM3RXbu09AwyLJQpVfjfx9UHQEREiryABg14Oj2dDj4VssvWp5zhdW9vYvbscWNk125t/FEivo3l7jZ1GRXqvFud89sHQC0AIiJS5MXt3Utw4yYXlQWX9yEufqebIrp2vxw9zasLYlm+I5F61ctzY80KeW/kAkoARESkyGvq58fm1JSLWgA2p6bQ1M/PjVFdnaSU87y7PJ7p6/ZSrrQ3YX9qxvCb/ClX2tst8SgBEBGRIi88IoKwUaN5GbL7AIxNOsn4Dye5O7Q8XUjPYMaPe3lnWTzJZy8wqF19/u+OJtSs5N75TZQAiIhIkZfV0e/xRx/lRGoqzfz9Gf/hpCLdAdBay8q4RCIWxLL76Bk6N6pBeJ/mNLuhsrtDA5QAiIhIMXHfkCFM/vhjoOgPBBT3azIRC2JYG3+MhjUq8OmwttzWrFauM6q6ixIAERERJzl2+hxvfbeT2Rv2Ualcacb2CeD+jn5uuc0vL0UvohJm7NixLFu27Ip1Vq1addmZ/ArD1KlTOXTokNP3+8QTT2SP49+kSROqVq0KQHR0NJ06daJFixa0bt2aOXPmXHE/8+bNwxiDM27nXLNmDcHBwZQqVYp58+Zdsv7UqVPUrVuXv/71r7luv2/fPrp27UqbNm1o3bo1CxcuLHBMIlL8nUtL59+rf6HrhFXM2bifBzr5s+qpUB7q3KBInvxBLQAu9/LLL+dZZ9WqVVSsWJGbbrop3/tNS0ujVCnnfHxTp06lZcuW1KlTxyn7y/L2229nv37vvffYsmULkDl18eeff07jxo05dOgQISEh9OjRIztByCk5OZmJEydedmKhnFatWsXUqVOvOCNh/fr1mTp1Km+88Uau659//nm6dOly2e0jIiIYOHAgo0aNIiYmhl69epGQkJBnbCJSMllrWbztV15btIN9x1O4rVkt/tmrOY1qVXR3aHkqmmlJMfD888/z7rvvZi8/99xzTJw48ZJ6w4cPz/6m6e/vzwsvvEBwcDCtWrVix44dJCQk8NFHH/H2228TFBTE2rVrOXr0KH/+859p164d7dq14/vvvwfgxRdfZMSIEXTv3p0HHniAI0eO0K9fPwIDAwkMDMxuRZgxY0b2TIUjR44kPT0dyJyE6MknnyQ4OJhu3bpx9OhR5s2bR1RUFEOGDCEoKIjU1FSXvF85pxRu0qQJjRs3BqBOnTrUqlWLo0eP5rrd888/z9NPP025cuWcEoe/vz+tW7fGy+vSX/1NmzZx5MgRunfvftntjTGcOnUKgJMnTzo9aRKR4mPbwZMMmvwjoyI3U760N9Mfbs9nw9sVi5M/kJm9lORHSEiIdYU9e/bYNm3aWGutTU9Ptw0bNrTHjh27pN6wYcPsF198Ya211s/Pz06cONFaa+0HH3xgH374YWuttS+88IKdMGFC9jb33nuvXbt2rbXW2r1799pmzZpl1wsODrYpKSnWWmsHDhxo3377bWuttWlpaTYpKcnGxMTYPn362PPnz1trrR01apSdNm2atdZawM6YMcNaa+1LL71kH3vsMWuttV26dLEbN27M9ed8/fXXbWBg4CWPxx9/PN/vVUJCgr3hhhtsWlraJevWr19vmzVrZtPT0y9Zt3nzZtu/f/88Y2zfvr0NDAy0N954o61WrVp2jIsXL75sTDk/F2szP8MuXbrYffv22SlTpmS/N3906NAh27JlS1u3bl1btWpVGxUVdcWfXUScq0uXLrZLly5ujaFT51usb0CI9Q/71rZ5eamdvi7BXki79H+YuwBRNh/nR10CuEb+/v5cd911bNmyhSNHjtCmTRuuu+66PLfr378/ACEhIfz3v//Ntc6yZcuIiYnJXj516hTJyckA9O3bl/LlywOwYsUKPv/8cwC8vb2pUqUK06dPZ9OmTbRr1w6A1NRUatWqBYCXlxeDBg0C4P7778+O5UrGjBnDmDFj8qx3JbNnz+aee+7B2/viwS4OHz7M0KFDmTZt2iXfyDMyMnjiiSeu2JyfZf369UD+LgFczqRJk+jVqxf16tW7Yr1Zs2YxfPhwnnzySdatW8fQoUPZtm1bri0KIlKypJ5P5+O1u9m6P4kMC2NvachjtzWicrnS7g7tmigBKIBHHnmEqVOn8uuvv/LQQw8B8OCDD7Jlyxbq1KmTawexsmUzB37w9vYmLS0t1/1mZGSwbt267BN9ThUqXHnISGstw4YN47XXXssz/vzcjjJhwgQic5lx69Zbb73kksdzzz3HggULgMyOfllmz57NBx98cFHdU6dO0bt3byIiIujYseMl+09OTmbbtm3Zc3//+uuv9O3bl/nz59O2bZ5DXF+1devWsXbtWiZNmsTp06c5f/48FStWZPz48RfV+/TTT1m8eDEAnTp14uzZsxw7diw7yRKRksday/yth/jXoh0cOnmWKuVLU7+6D8/2au7u0ApEX1sKoF+/fixevJiNGzfSo0cPAKZMmUJ0dPRV9Q6vVKlS9jd8gO7du/P+++9nL+c8mebUrVs3PvzwQwDS09M5deoU3bp1Y968eSQmJgJw/Phx9u7dC2QmFln9EWbOnEnnzp1zPX5OY8aMITo6+pJHbv0dxo0bl70+S1xcHCdOnKBTp07ZZefPn6dfv3488MADDBgwINfjVqlShWPHjpGQkEBCQgIdO3bM8+QfGhp6Td/+ASIjI9m3bx8JCQm88cYbPPDAA5ec/CGzE+Hy5csBiI2N5ezZs9SsWfOajikiRd/mfSfoN+kH/j47mmoVyjB7REeaXF/JbcP3OpMSgAIoU6YMXbt2ZeDAgZc0b1+NO++8ky+//DK7E+DEiROJioqidevWBAQE8NFHH+W63bvvvsvKlStp1aoVISEhbN++nYCAACIiIujevTutW7fmjjvu4PDhw0Bm68H27dsJCQlhxYoVjB07FsjsqPiXv/zFJZ0AZ82axeDBgy9qbZg7dy5r1qxh6tSp2bcJZiUNY8eOZf78+Vd1jA4dOmTvJ+djyZIll9TduHEjvr6+fPHFF4wcOZIWLVrkuf+cMb355pt8/PHHBAYGcu+99zJ16tQiNbCHiDjHwaRU/jZrC/0n/cDBpFQm3NOab/7amY4N877UW1xoOuACyMjIIDg4mC+++CK7V3tRVrFiRU6fPu3uMERErlnWZUFXjQR45lwaH676hY/X7gZgxK0N+UuXG6lQ9vcr5q6OoaA0HbCLxcTE0KdPH/r161csTv4iInJ5GRmWeZsPMGFJHEeTz9E3sA7P/KkZdate2herpI+9yAYAACAASURBVFACcI0CAgLYvXu3u8O4Kvr2LyJyqR93/8Yr38aw/dApgupV5d9DQwiuX83dYbmcEgAREfFIe387w2sLd7B4+6/UqVKOdwcH0Tewjsf061ECICIiHuXU2Qu8v2IXU79PoJS34ck7mvDorQ1LRM/+q+H2BMAY4w1EAQettX2MMdWBOYA/kAAMtNaecNR9FngYSAf+Zq29tJu3iIhILtLSM5i1cT9vf7eTEynnuSfYl6d6NOX6ys4Zary4KQq3Af4diM2xHAYst9Y2BpY7ljHGBACDgRZAT2CSI3kQERG5ojU7j9Jr4lqe/2objWpV5Ju/dmbCgMCrPvnPjIzk5w0bWLN6NQENGjAzl4HSigu3tgAYY3yB3sA44P8cxXcBoY7X04BVwDOO8tnW2nPAHmPMLqA9sK4QQxYRkWJkV+Jpxi2IYWXcUepX9+Gj+4Pp0eKGa7rOPzMykrBRo3mzRg2Cy/uwOTWFsFGjAbhvyBBnh+5y7r4E8A7wNFApR9n11trDANbaw8aYrDFW6wI/5qh3wFEmIiJykRNnzvPu8nim/7gXn9LePPunZgy/2Z+ypa694TgiPJyXq1ahg0/mkOwdfCrwsqNcCcBVMMb0ARKttZuMMaH52SSXslxHMTLGjABGQObQrSIiUjLkNfjO+bQMpv+4l4nL40k+e4F729fniTuaUKNi2QIfO27vXoIbN7moLLi8D3HxOwu8b3dwZwvAzUBfY0wvoBxQ2RgzAzhijKnt+PZfG0h01D8A5JyqzRc4lNuOrbWTgcmQORKgq34AEREpGqy1LI9N5NWFsew+doZbGtcgvHcATW+olPfG+dTUz4/NqSnZLQAAm1NTaOrn57RjFCa3dQK01j5rrfW11vqT2blvhbX2fmA+MMxRbRjwteP1fGCwMaasMaYB0BjYUMhhi4hIERN7+BT3f7qeRz6PAgOfDW/L5w+1d+rJHyA8IoKxSSdZn3KGC9ayPuUMY5NOEh4R4dTjFJaicBfAH40H7jDGxAN3OJax1m4H5gIxwGLgMWttutuivEYvvvgib7zxxmXXr1q1ij59+uS67pFHHiEmJsZVoeXqiy++oEWLFnh5eXGlORUWL15M06ZNadSo0UWz6D3//PO0bt2aoKAgunfvzqFDuTbaiIhctaPJ53j2vz/Re+Jath08xYt3BrDkH7dyW7PrXTKYz31DhjD+w0k8eewYQTvjeN3bm/EfTiqW1/9BkwE5hbUWay1eXnnnUy+++CIVK1bkqaeeynX9qlWreOONN/j222+dHeY1iY2NxcvLi5EjR/LGG2/kOh1veno6TZo04bvvvsPX15d27doxa9YsAgICOHXqFJUrVwZg4sSJxMTEXHZ2QxGR/Dh7IZ0p3yfwwcpdnL2QztBOfvy9W2Oq+pQplOOXlMmAimILQLGQkJBA8+bNGT16NMHBwezfv58JEybQrl07WrduzQsvvJBdd9y4cTRt2pTbb7+duLi47PJdu3Zx++23ExgYSHBwML/88guQOWb/PffcQ7NmzRgyZAhZSVpoaChRUVGkp6czfPhwWrZsSatWrXj77bcB2LRpE4GBgXTq1IkxY8bQsmXLAv+czZs3p2nTpless2HDBho1akTDhg0pU6YMgwcP5uuvM6/cZJ38Ac6cOeMxQ2yKiPNZa1n482HueHs1/1q8g44Nq7PkiVt54c4WhXbyL0ncfRtgsRYXF8eUKVOYNGkSS5cuJT4+ng0bNmCtpW/fvqxZs4YKFSowe/ZstmzZQlpaGsHBwYSEhAAwZMgQwsLC6NevH2fPniUjI4P9+/ezZcsWtm/fTp06dbj55pv5/vvv6dy5c/Zxo6OjOXjwINu2bQMgKSkJgAcffJD33nuPLl26MGbMmFxjTk5O5pZbbsl13cyZMwkICLjq9+HgwYPUq/d7/0xfX1/Wr1+fvfzcc8/x+eefU6VKFVauXHnV+xcR+fnASV75NoYNCcdpdkMlZjzcgc6Na7g7rGJNCUAB+Pn50bFjRwCWLl3K0qVLadOmDZD5LT4+Pp7k5GT69euHj48PAH379gUyT8QHDx6kX79+AJQr9/toVO3bt8fX1xeAoKAgEhISLkoAGjZsyO7du3n88cfp3bs33bt35+TJkyQlJdGlSxcAhg4dyqJFiy6JuVKlSkRHRzv1fcjtMlLOb/rjxo1j3LhxvPbaa7z//vu89NJLTj2+iJRcv548y+tLdvDfzQe5rkIZXu3XikHt6uHtpdbEglICUAAVKvx+K4i1lmeffZaRI0deVOedd97Jtdn7Sn0vypb9/X5Vb29v0tLSLlpfrVo1tm7dypIlS/jggw+YO3cub731Vr6a113RAuDr68v+/fuzlw8cOECdOnUuqXfffffRu3dvJQAikqfU8+lMXrObj1b/QnqGZWSXhjzWtRGVy5V2d2glhvoAOEmPHj347LPPOH36NJDZLJ6YmMitt97Kl19+SWpqKsnJyXzzzTdA5rVxX19fvvrqKwDOnTtHSkpKvo517NgxMjIy+POf/8wrr7zC5s2bqVq1KlWqVOF///sfAJGXGZ86qwUgt8e1nPwB2rVrR3x8PHv27OH8+fPMnj07u6UjPj4+u978+fNp1qzZNR1DRDxDRoblyy0HuO3NVby9bCddm9Vk2f914dk/NdfJ38nUAuAk3bt3JzY2lk6dOgFQsWJFZsyYQXBwMIMGDSIoKAg/P7+Lvn1Pnz6dkSNHMnbsWEqXLs0XX3yRr2MdPHiQBx98kIyMDABee+01AKZMmcJDDz2Ej48PPXr0cMrP9eWXX/L4449z9OhRevfuTVBQEEuWLOHQoUM88sgjLFy4kFKlSvH+++/To0cP0tPTeeihh2jRogUAYWFhxMXF4eXlhZ+fn+4AEJHL2rT3OC9/G8vW/Um0qluFdwYF0aHhde4Oq8TSbYAlVEJCAn369MnuKCgiUlQdOJHC+EU7+Panw1xfuSxjejSjf5u6eBXR6/wl5TZAtQCIiIhbnD6XxqSVu/jkf3vwMvC3bo35S5eG+JTRqakw6F0uofz9/fXtX0SKpPQMy7xN+5mwZCfHTp/j7qA6PN2zGXWqlnd3aB5FCYCIiBSaH345RsS3scQcPkVw/ap8/EAIbepXc3dYV6WoNv1fLSUAIiLicgnHzvDqwliWxhyhbtXyTLy3DXe2rq3RQd1ICYCIiLjMydQLvL8inqk/JFDG24sxPZrycOcGlCvt7e7QPJ4SABERcbq09AxmbdjH28viOZFyngEhvjzVvSm1KpfLe2MpFEoARETEqVbFJTJuQSzxiafp2LA64b0DaFm3irvDkj9QAiAiIk6xKzGZiAWxrIo7it91Pvx7aAjdA67Xdf4iSgmAiIgUyPEz53ln2U4i1+/Dp4w3z/VqzgM3+VG2lK7zF2VKAERE5JqcT8vg83UJTFwez5nz6dzXvj7/uL0x11Usm+e24n5KAERE5KpYa/ku5givLdrBnmNnuLVJTcJ7N6fJ9ZXcHZpcBc0GWEKdO3eOQYMG0ahRIzp06EBCQsIldZKTkwkKCsp+1KhRg3/84x8A7Nu3j65du9KmTRtat27NwoULC/knEMlbaGho9rjsUnD5eT9jDp1iyCfrGTF9E95ehikPtuPzh9rr5F8MqQXACay1WGvx8io6+dSnn35KtWrV2LVrF7Nnz+aZZ55hzpw5F9XJmho4S0hICP379wcgIiKCgQMHMmrUKGJiYujVq1euSYSIeIbE5LO8tXQnc6L2U6V8aV7q24L7OtSntHfR+b8nV0ef3DVKSEigefPmjB49muDgYPbv38/rr79Oq1atCAwMJCwsDIDo6Gg6duxI69at6devHydOnCA2Npb27dtftK/WrVs7Nb6vv/6aYcOGAXDPPfewfPlyrjTzY3x8PImJidnTFRtjOHXqFAAnT56kTp06To1PRIqHsxfS+WDlLrpOWMW8TQd46OYGrH6qK8Nu8tfJv5hTC0ABxMXFMWXKFCZNmsSiRYv46quvWL9+PT4+Phw/fhyABx54gPfee48uXbowduxYXnrpJd555x3Onz/P7t27adiwIXPmzGHgwIGX7H/ChAlERkZeUn7rrbcyceLEK8Z28OBB6tWrB0CpUqWoUqUKv/32GzVq1Mi1/qxZsxg0aFD27Tovvvgi3bt357333uPMmTMsW7bsqt4bESnerLUs+Pkw4xft4MCJVO4IuJ5/9mpOgxoV3B2aOIkSgALw8/OjY8eOACxbtowHH3wQHx8fAKpXr87JkydJSkqiS5cuAAwbNowBAwYAMHDgQObOnUtYWBhz5sy5pHkeYMyYMYwZM+aaYsvt2/6V7sWdPXs206dPz16eNWsWw4cP58knn2TdunUMHTqUbdu2FanLHCLiGlv3J/HKtzFE7T1BsxsqMfORDtzUKPcvD1J8KQEogAoVfs+ErbVXNdjFoEGDGDBgAP3798cYQ+PGjS+pczUtAM899xwLFiwAMi87+Pr6sn//fnx9fUlLS+PkyZNUr14911i2bt1KWloaISEh2WWffvopixcvBqBTp06cPXuWY8eOUatWrXz/jCJSvJxPy2D/8RTu+uB7alQsw/j+rRjQth7eXhrIpyTS1zkn6d69O5999hkpKSkAHD9+nCpVqlCtWjXWrl0LwPTp07NbA2688Ua8vb155ZVXGDRoUK77HDNmDNHR0Zc8cmv+HzduXPZ6gL59+zJt2jQA5s2bx2233XbZBGXWrFnce++9F5XVr1+f5cuXAxAbG8vZs2epWbPm1b4tIlIMpJxP4+3vdhK9P4nfzpxnVOiNrHwqlMHt6+vkX4KpBcBJevbsSXR0NG3btqVMmTL06tWLV199lWnTpvGXv/yFlJQUGjZsyJQpU7K3GTRoEGPGjGHPnj1Oj+fhhx9m6NChNGrUiOrVqzN79uzsdUFBQRf1/p87d+4lt/m9+eabPProo7z99tsYY5g6daqG8xQpYTIyLF9FH+T1xXH8euos1SqUoX51H57p2czdoUkhMFfqGV4StG3b1kZFRbk7DBFxgax71letWuXWOIqjqITjvPxtDD8dOElr3yo83yeAMcMzbwPW+1m8GWM2WWvb5lVPLQAiIh5k//EUxi/ewYKfDnND5XK8NTCQu4Pq4qWmfo+jBEBExAMkn73ApFW/8On/9uBl4O/dGjOyS0N8yug04Kn0yYuIlGDpGZa5Uft5c2kcx06fp3+buozp2ZTaVcpfVG9mZCQ/b9jAidRUAho0IDwigvuGDHFT1FIYlACIiJRQP+w6xsvfxrDj12RC/KrxybB2BNWrekm9mZGRhI0azZs1ahBc3ofNqSmEjRoNoCSgBFMnQBEpttQJMHd7jp1h3IJYlsUeoW7V8jzbqxm9W9W+7J08AQ0a8HR6Oh18fh/bZH3KGV739ibGBXcpiWupE6CIiIc5mXKBiSvi+XxdAmW8vXi6Z1MeurkB5Up7X3G7uL17CW7c5KKy4PI+xMXvdGG04m5KAEREirkL6RnMXL+Pd5btJCn1AoPa1uP/ujehVqVy+dq+qZ8fm1NTLmoB2JyaQlM/P1eFLEWAEgARkWJsZVwi4xbEsivxNJ0aXsfzfQIIqFP5qvYRHhFB2KjRvAzZfQDGJp1k/IeTXBO0FAlKAEREiqGdR5KJWBDLmp1H8b/Oh8lDQ7gj4PprGrEzq6Pf448+yonUVJr5+zP+w0nqAFjCKQEQESlGfjt9jreX7WTWhv1UKONNeO/mPNDJnzKlCja1y31DhjD5448Bdar0FEoARERcyFl3KpxPy2DaDwlMXBFPyvl07u9Qn7/f3oTqFcoUPEjxSG5LAIwx5YA1QFlHHPOstS8YY6oDcwB/IAEYaK094djmWeBhIB34m7V2iRtCFxEpNNZalsYc4dWFsez9LYXQpjV5rldzGl9fyd2hSTHnzumAzwG3WWsDgSCgpzGmIxAGLLfWNgaWO5YxxgQAg4EWQE9gkjHmyve2iEiJlTVy3ZrVqwlo0ICZkZHuDsnpth86yb0f/8jI6Zso4+3F1AfbMfXB9jr5i1O4rQXAZo5AdNqxWNrxsMBdQKijfBqwCnjGUT7bWnsO2GOM2QW0B9YVXtQiUhSU9JHrEk+d5Y2lcXyx6QBVy5fmlbtacG/7+pTydud3Nilp3NoHwPENfhPQCPjAWrveGHO9tfYwgLX2sDGmlqN6XeDHHJsfcJSJiIeJCA/n5apVsu9b7+BTgZcd5cU5ATh7IZ1P/7eHD1bu4kJ6Bo90bsBfb2tMlfKl3R2alEBuTQCstelAkDGmKvClMablFarndm9LruMYG2NGACMA6tevX+A4RaRoKWkj11lr+eanw/xr0Q4OJqXSPeB6/tmrOf41KuS9scg1KhJ3AVhrk4wxq8i8tn/EGFPb8e2/NpDoqHYAqJdjM1/g0GX2NxmYDJlzAbgscBFxi5I0cl30/iRe+TaGTXtPEFC7MhMGtOamG2u4OyzxAG67oGSMqen45o8xpjxwO7ADmA8Mc1QbBnzteD0fGGyMKWuMaQA0BjYUbtQiUhSER0QwNukk61POcMFa1qecYWzSScIjItwdWr4dSkrlH7O3cPcH37P3txT+9edWfPN4Z538pdC4swWgNjDN0Q/AC5hrrf3WGLMOmGuMeRjYBwwAsNZuN8bMBWKANOAxxyUEEfEwxXnkujPn0vj36l+YvHY3GRYe63ojo0IbUbFskWiQFQ+i6YBFpNgqDtMBZ8W4YsVK/rvlIBOW7ODIqXP0aV2bZ3o2o151H/cGKCWOpgMWESkiks+mcdcH3/PzwZME1qvKpCHBhPhVd3dY4uGUAIiIuMj+4ynEH0nmtzPnKXP6HO8MCqJvYB28vK5+wh4RZ1MCICLiZMlnL/D+yl1M+V8CJ1Iu4FvNhxVPhlK+jAYvlaJDCYCIiJOkZ1jmbNzPW9/Fcez0ef4c7Mu61VUpU8pLJ38pcpQAiIg4wf/ijxGxIIYdvybTzr8anw1vR2vfqoRO0vC9UjQpARARKYBfjp7m1QWxLN+RSL3q5Zk0JJg/tbwBY0z2hEUnUlMJaNCA8IiIYnGrongGJQAiItcgKeU87y6PZ/q6vZQr7U3Yn5ox/CZ/ypXObOov6RMWSfGncQBEpNhyxzgAF9IziPxxL+8sj+dU6gUGtavP/93RhJqVyl5UL6BBA55OT79ouOL1KWd43dubmD17Ci1e8TwaB0BErklxGFzHHay1rIxLZNyCWH45eoabG11HeO8AmteunGv9kjZhkZQ8SgBERPIQ92syEQtiWBt/jIY1KvDJA23p1rwWxlz+fv6SNGGRlExKAERELuO30+d467udzNqwj0rlSjO2TwD3d/SjTKm8e/aHR0QQNmo0L0N2H4CxSScZ/+Ek1wcukg9OSQCMMW8AU6y1252xPxERdzqXls7U7xN4f8UuUi6k80Anf/7erTHVKpTJ9z6K84RF4hmc1QKwA5hsjCkFTAFmWWtPOmnfIiKFwlrLku2/8urCHew7nsJtzWrxz17NaVSr4jXt774hQ5j88ceA+lRI0eOUBMBa+wnwiTGmKfAg8JMx5nvgY2vtSmccQ0TElbYdPMkr38awfs9xmlxfkc8fas+tTWq6OywRl3FaHwBjjDfQzPE4BmwF/s8YM9JaO9hZxxERcaYjp84yYUkc/9l8gGo+ZYi4uyWD29WjlLdG8JOSzVl9AN4C7gRWAK9aazc4Vv3LGBPnjGOIiDjT2QvpfLxmNx+u/oUL6Rk8ektDHuvaiCrlS7s7NJFC4awWgG1AuLU2JZd17Z10DBGRArPWMn/rIf61aAeHTp6lZ4sbeLZXM/yuq5D3xiIliLMSgCHW2s9yFhhjlltru6kzoIgUFZv3neCVb2PYsi+JFnUq89agIDo2vM7dYYm4RYESAGNMOcAHqGGMqQZkjYpRGahTwNhERK4ovz3rDyal8q9FO5i/9RA1K5Xl9Xta8+dgX7y9Lj+Qj0hJV9AWgJHAP8g82W/OUX4K+KCA+xYRKZAz59L4aPUvTF6zG4C/dm3EqNAbqVBWY6CJFOivwFr7LvCuMeZxa+17TopJRKRAMjIs8zYf4I0lcSQmn6NvYB2e+VMz6lYt7+7QRIqMgl4CuM1auwI4aIzp/8f11tr/FmT/IiJXa/3u33hlQQzbDp4iqF5VPhoaQnD9au4OS6TIKWg7WBcyb/27M5d1FlACICKFYu9vZ3ht4Q4Wb/+VOlXK8e7gIPoG1rnihD0inqyglwBecDw/6JxwRESuzqmzF/hgxS6mfJ9AKW/Dk3c04ZFbGlK+jLe7QxMp0pwy1JUx5lVjTNUcy9WMMRHO2LeIFJ6ZkZH8vGEDa1avJqBBA2ZGRro7pMtKS89gxo976TphFZPX7qZvUB1WPhXK490a6+Qvkg/O6gr7J2vtP7MWrLUnjDG9gHAn7V9EXGxmZCRho0bzZo0a2dPXho0aDVDkZrBbs/MoEQti2HnkNO0bVGdanwBa1q3i7rBEihVnJQDexpiy1tpzAMaY8kBZJ+1bRApBRHg4L1etQgefzBHxOvhU4GVHeVFJAHYlnubVhbGs2JFI/eo+fHR/MD1a3KDr/CLXwFkJwAxguTFmCpmd/x4Cpjlp3yJSCOL27iW4cZOLyoLL+xAXv9NNEf3uxJnzvLs8nhk/7qV8aW+e/VMzht/sT9lSRb+pX9MAS1HlrOmAXzfG/ATc7ih6xVq7xBn7FpHC0dTPj82pKdktAACbU1No6ufntpgupGcwfd1e3l0eT/LZC9zbvj5P3NGEGhXVwChSUM4cDmsLUJrMFoAtTtyviBSC8IgIwkaN5mXI7gMwNukk4z+cVOixWGtZHpvIqwtj2X3sDLc0rsFzvZvT7IbKhR6LSEnlrOmABwITgFVkzgfwnjFmjLV2njP2LyKul3Wd//FHH+VEairN/P0Z/+GkQr/+H3v4FBELYvh+1280rFmBz4a3pWvTWrrOL+Jkxlpb8J0YsxW4w1qb6FiuCSyz1gYWeOcF1LZtWxsVFeXuMESKjdDQUKDwr10fO32ON5fuZM7GfVQqV5onbm/MkI5+lPZ2yt3KIh7DGLPJWts2r3rOugTglXXyd/gNJ40xICIl27m0dKZ8n8D7K3Zx9kI6w27y5+/dGlPVp4y7QxMp0ZyVACw2xiwBZjmWBwELnbRvESmBrLUs2vYrry2KZf/xVLo1q8U/ezfnxpoV3R2aiEdw1l0AY4wxfwZuJrMPwGRr7ZfO2LeIlDw/HzjJK9/GsCHhOE2vr8SMhzvQuXENd4cl4lGcdheAtfY/wH+ctT8RKXmOnDrL64vj+M/mA1xXoQzj+rVkUNt6lNJ1fpFCV9DpgJPJvO3vklWAtdbqnh0RIfV8OpPX7Oaj1b+QnmEZ2aUhj3VtROVypd0dmojHKuhsgJWudVtjTD3gc+AGIIPMywbvGmOqA3MAfyABGGitPeHY5lngYSAd+JsGGxIp2jIyLPO3HuJfi3dw+ORZerW6gbCezal/nY+7QxPxeE67BGCM6Qw0ttZOMcbUACpZa/dcYZM04Elr7WZjTCVgkzHmO2A4sNxaO94YEwaEAc8YYwKAwUALoA6wzBjTxFqb7qyfQUScZ9PeE7z8bQxb9yfRsm5l3hkURIeG17k7LBFxcNZAQC8AbYGmwBSgDJnzA9x8uW2stYeBw47XycaYWKAucBcQ6qg2jczBhZ5xlM92TDi0xxizC2gPrHPGzyAiznHgRAr/WhzHN1sPUatSWd4YEEj/NnXx8tJAPiJFibNaAPoBbYDNANbaQ45v9flijPF3bL8euN6RHGCtPWyMqeWoVhf4McdmBxxlIlIEnD6XxoerdvHJ2j0YA3+7rREju9xIhbLOHHFcRJzFWX+Z56211hhjAYwxFfLaIIsxpiKZdw/8w1p76grDfea2ItdhDI0xI4ARAPXr189vKCJyDdIzLP/ZdIAJS+M4mnyOu4Pq8HTPZtSpWt7doYnIFTgrAZhrjPk3UNUY8yiZ0wF/nNdGxpjSZJ78I621/3UUHzHG1HZ8+68NZI0weACol2NzX+BQbvu11k4GJkPmUMDX8gOJSN7W/fIbr3wbQ8zhUwTXr8rkoSG0qV/N3WGJSD44KwHIANYCp4AmwFhr7XdX2sBkftX/FIi11r6VY9V8YBgw3vH8dY7ymcaYt8jsBNgY2OCk+EXkKiQcO8OrC2NZGnOEulXLM/HeNtzZurYm7BEpRpyVAFQi8/a848Bs4Kd8bHMzMBT42RgT7Sj7J5kn/rnGmIeBfcAAAGvtdmPMXCCGzDsIHtMdACKF62TqBd5fEc/UHxIo7e3FmB5NebhzA8qV9nZ3aCJylZw1FPBLwEvGmNZkzgOw2hhzwFp7+xW2+R+5X9cH6HaZbcYB4woar4hcnbT0DGZt2Mfby+I5kXKeASG+PNW9KbUql3N3aCJyjZzdPTcR+JXM2QBr5VFXRIqgP04DvHrnUSK+jSE+8TQdG1YnvHcALetWcU9wIuI0zhoHYBSZ3/xrAvOAR621Mc7Yt4i4x67EZCIWxLIq7ih+1/nw76EhdA+4Xtf5RUoIZ7UA+JF5G190njVFpEg7ceY87yzbyYz1+/Ap481zvZrzwE1+lC2l6/wiJYmz+gCEOWM/IuI+59My+HxdAhOXx3P6XBr3dajPE7c34bqKZd0dmoi4gIboEvFw1lq+iznCa4t2sOfYGW5pXIPn+wTQ5PprnutLRIoBJQAiHizm0CkiFsTwwy+/cWPNCkwZ3o7QpjV1nV/EAygBEPFAR5PP8ebSOOZE7adK+dK81LcF93WoT2lvL3eHJiKFRAmAiAc5eyGdz77fw6SVv3D2QjoP3dyAv93WmCo+pd0dmogUMiUAIh7AWsvCn3/ltUWxHDiRyu3Nr+efvZrR8P/bu+/4qsrDj+OfJyEJSZghYQgJCSsQkRmliCMWF2Jr1eKAOvpzFcdLbUtFf2j7kyAUxdUKitbVBluttlqmiAwHMgWBdDjMJwAAFVRJREFUhAwgEGZYCQk3hOTe5/dHrhYVMOvm3PF9/5Obk3tvvjyc5HxznjMSWjgdTUQcogIgEuTWF5WQOSebVYWH6d2xJVm3D2FYj3inY4mIw1QARILUntIKnpyfy3tf7iK+RSSTrzmL69ITCQ/TAX4iogIgEnRcx6uZuWwrLy7dgsfC2Izu3J3RnZbNNc8vIv+lAiASJDwey7/X7WLq/Fz2HjnGyLM6MX5EbxLjYpyOJiJ+SAVAJAisLjzExNnZrN9ZSr8urfnT6IGcnRzndCwR8WMqACIBrOiQiynzNzPnqz10aBXFtFH9uXpgZ8I0zy8iP0AFQCQAlR2rYvqSLfzl022EGbh/eE/uurAbMZH6kRaR2tFvC5EA4vZY3lldxFMf5nGgvJJrBnZm3OWpdGod7XQ0EQkwKgAiAeLzggNMnJNDzp4jDO7allduSWdAYhunY4lIgFIBEHFQRkYGAEuWLDnlc7YdOMoTc3NYmL2Pzm2i+fPogYw8q5Nu2CMiDaICIOKnSl1VPP9xPm8uLyQyPIxxl6Vy23kpNI8IdzqaiAQBFQARP1Pt9jBr5Q6eWZhHSUUV16cn8utLe9G+ZXOno4lIEFEBEPEji3OLmTQnh4LicoZ2a8ejV6aRdkYrp2OJSBBSARDxA3n7ysick8OyvP0kt4th5k2DuSStg+b5RcRnVAAkaNTmgDp/U+22FB12MeK5T4iJDGfCyD7cPDSZyGZhTkcTkSCnAiDigOPVHt5cXsiXRYfxWBg3JIkHLu5FXGyk09FEJESoAIg0IWstH2bvY/LcHAoPumjZPIKucTE8flVfp6OJSIjRfkaRJrJpdyk3vvwFd/11DRHhYbz+y7Pp3bEl0ZE6rU9Emp72AIj4WHHZMaYtyOPtNUW0iY5g4lVncuM5STQLV/8WEeeoAIj4yLEqN3/5dBvTFxdw3O3htmEp3De8J62jI5yOJiKiAiDS2Ky1zP5qD1PmbWZXSQWXpnXg4Sv6kBIf63Q0EZFvqACINKJ1RSVMnJ3Nmu2H6dOpFU+O6se53eOdjiUi8j0qACKNYHdJBVPnb+bf63YT3yKKP157Fj8fnEh4mC7kIyL+SQVApAFcx6t5celWZi7bgsfCPRd1Z2xGD1pE6UdLRPybfkuJ1IPHY3nvy108uWAz+45UcmW/Tjx0eW8S42KcjiYiUis6D0mCwqysLDasXMmypUtJS0lhVlaWz77XqsJD/Gz6Z/z2nfV0bB3Nu2OH8ufRg+q88W/KzCIi36U9ABLwZmVlMX7s3UyLj2dQdAxrK1yMH3s3AKPHjGm071N0yMXkeTnM3bCXjq2a88z1/bmqf2fC6jHP31SZRUROxVhrnc7gU+np6Xb16tVOxxAfSktJ4XduN0Ni/nua3QrXUaaGh5O9bVuD37/sWBUvLN7Cq59uIzzM8KsLu3PnBd0adAU/X2cWkdBljFljrU3/oedpD4AEvNzt2xnUs9e3lg2KjiE3P69B7+v2WP6xqoinF+ZyoPw41w7qwrjLUunYunmD3hd8l1lEpLYcPQbAGPOqMabYGLPxhGVxxpiFxph878e2J3ztYWNMgTEm1xhzmTOpxd+kdu3K2grXt5atrXCR2rVrvd/zs4IDjHz+Ex751wZS4mP54N5hTLuuf6Ns/ME3mUVE6sLpgwBfBy7/zrLxwCJrbU9gkfdzjDFpwA3Amd7XTDfG6C4qwoTMTB4rKWWF6yhV1rLCdZTHSkqZkJlZ5/faur+c299YxZhXVlBeWc30MYN4+66h9OvSxm8zi4jUh6NTANbaZcaY5O8svgrI8D5+A1gCPORd/ndrbSWwzRhTAJwDLG+KrOK/vj5o7r477uBwRQW9k5OZMmN6nQ6mK3VV8dyifN5cXkjziHAeurw3vxyWTPMI33TMxsgsItIQ/ngMQAdr7R4Aa+0eY0x77/LOwBcnPG+nd5kIo8eMYebLLwOwZMmSWr+uyu0h64vtPLsonyMVVVx/dhK/vqQXCS2jfJT0v+qbWUSkMfhjATiVk51rddJTGIwxdwJ3AiQlJfkykwQoay2Lc4uZNCeHLfuPMqxHOyaMTKNPp1ZORxMRaRL+WAD2GWM6ef/67wQUe5fvBBJPeF4XYPfJ3sBaOxOYCTWnAfoyrASe3L1lZM7J5pP8A3SLj+WVm9MZ3qc9xui6/SISOvyxAHwA3AJM8X58/4Tls4wxTwNnAD2BlY4klIB0sLySpxfm8dbKHbSIasajV6Zx04+6EtnM6WNhRUSanqMFwBjzFjUH/MUbY3YCv6dmw/+2MeY2YAcwCsBau8kY8zaQDVQD91hr3Y4El4BSWe3mjc8L+dOiAlxVbm4emsz9w3vSNjbS6WgiIo5x+iyAG0/xpeGneP4kYJLvEkkwsdayYNNeJs/bzPaDLi5KTeB/R/ahR/uWTkcTEXGcP04BiDTYxl2lTJydzYpth+jVoQVv/s85XNArwelYIiJ+QwVAgkqV28O4d9bzz7U7aRsTSebP+nLD2Yk0C9c8v4jIiVQAJCgcq3Kzq6SC3SUVFK/bxR3nd+Oei3rQOjrC6WgiIn5JBUACmrWWD9bvZur8XIoOuYiLjWThgxeSHB/7wy8WEQlhKgASsNbuOMzE2dl8uaOEM89oxeLFSxjavZ3TsUREAoIKgAScXSUVTJ2/mffX7SahZRRTf96Pawd1ITxMF/IREaktFQAJGEcrq3lx6RZmLtsKwL0X9WBsRndio7Qai4jUlX5zit/zeCzvrt3JkwtyKS6r5Kf9z+ChEb3p3Cba6WgiIgFLBUD82oqtB5k4J5uNu44wILENM34xmMFd2zodS0Qk4KkAiF/acdDF5Hk5zNu4lzNaN+e5Gwbwk35nEBZk8/y6DbCIOEUFQPzKkWNVvPBxAa99VkizcMNvLunF7ed3Izoy3OloIiJBRQVA/EK128PfVxXxzMI8DrmOc+2gLoy7LJUOrZo7HU1EJCipAIjjPsnfT+bsHHL3lXFOShxvXJlG386tnY4lIhLUVADEMQXF5TwxN4ePNxeTGBfNjDGDuLxvR4wJrnl+ERF/pAIgTa7EdZxnP8rnb19sJzoinIdH9ObWYclENdM8v4hIU1EBkCZT5fbw1+XbeW5RPmXHqrjxnCQevKQX8S2inI4mIhJyVADE56y1fLy5mElzcth64Cjn94znf0f2oXfHVk5HExEJWSoA4lOb9x4hc3YOnxYcoFtCLK/ems5Fqe01zy8i4jAVAPGJA+WVPL0wj7+v3EHL5hH8/idp/OJHXYkID3M6moiIoAIgjayy2s1rnxXywscFVFS5uXloMg9c3JM2MZFORxMRkROoAEijsNYyb+NeJs/LoehQBcN7t+eRkX3ontDC6WgiInISKgDSYBt2ljJxdjYrCw+R2qElf7ttCOf1jHc6loiInIYKgNTbviPHmDo/l/e+3ElcTCSTru7L9emJNNM8v4iI31MBkDqrOO7m5U+2MmPJFtwey50XdOOei3rQqnmE09FERKSWVACk1jweywfrd/PH+ZvZU3qMEX078vCIPiS1i3E6moiI1JEKgNTKmu2HmTg7m3VFJfTt3Ipnrx/AkG7tnI4lIiL1pAIgpzX0vAvYcchFxE//j/Yto3hqVH+uGdiZsDBdyEdEJJCpAMhJlVdWM2NJAeuLSgDI/HEP7rqwO7FRWmVERIKBfpvLt7g9lnfX7OTJD3PZX1ZJXGwkSXEx/PrSVKejiYhII1IBkG8s33KQibOzyd5zhIFJbZh502Ae/EIX8hERCUYqAML2g0d5Ym4OCzbto3ObaJ6/cSA/6ddJN+wREQliKgAhrLSiij9/nM/rnxcSER7GuMtSue28FJpHhDsdTUREfEwFIARVuz28taqIZxbmcdh1nFGDu/DbS1Np36q509FERKSJqACEmKV5+5k0J5u8feUMSYnj0SvT6Nu5tdOxRESkiakAhIiC4jImzclhce5+uraL4cVfDOayMztonl9EJESpAAS5w0eP8+xHefxtxQ5iIsJ55Ire3HJuMlHNNM8vIhLKVACC1PFqD28uL+T5RfmUV1YzekgSD17ci3YtopyOJiIifiDgCoAx5nLgOSAceMVaO8XhSH7FWstHOcU8MTeHbQeOcn7PeB69Mo1eHVrW+b1mZWWxYeVKDldUkJaSwoTMTEaPGeOD1CIi0tQCqgAYY8KBF4BLgJ3AKmPMB9babGeT+YecPUeYODubz7ccpHtCLK/dejYZqQn1mueflZXF+LF3My0+nkHRMaytcDF+7N0AKgEiIkHAWGudzlBrxpihwB+stZd5P38YwFo7+VSvSU9Pt6tXr26ihM7YX1bJ0wtz+ceqIlpFR/Dgxb0YPSSJiPCwer9nWkoKv3O7GRIT+82yFa6jTA0PJ3vbtsaILSIiPmCMWWOtTf+h5wXUHgCgM1B0wuc7gSEOZXHcsSo3r31WyAuLCzhW5ebWc1O4f3hPWsdENPi9c7dvZ1DPXt9aNig6htz8vAa/t4iIOC/QCsDJ9mV/bxeGMeZO4E6ApKQkX2dqctZa5m7Yy+R5Oew8XMHFfTrwyBW96ZbQeNftT+3albUVrm/tAVhb4SK1a9dG+x4iIuKc+u8jdsZOIPGEz7sAu7/7JGvtTGtturU2PSEhocnCNYWvdpZw3UvLuWfWWlpENSPr9iG8ckt6o278ASZkZvJYSSkrXEepspYVrqM8VlLKhMzMRv0+IiLijEDbA7AK6GmMSQF2ATcAo52N1DT2lh5j6oLNvLd2F/EtIpl8zVlcl55IeJhvLuTz9YF+991xB4crKuidnMyUGdN1AKCISJAIqAJgra02xtwLLKDmNMBXrbWbHI7lUxXH3by0bAsvLd2K22P51YXdueei7rRs3vB5/h8yeswYZr78MgBLlizx+fcTEZGmE1AFAMBaOxeY63QOX/N4LO+v38Uf5+Wy98gxRp7VifEjepMYF+N0NBERCQIBVwCCQUZGBnDqv6rXbD/E4//JZv3OUvp1ac2fRg/k7OS4pgsoIiJBTwXAjxQdcjFl/mbmfLWHDq2imDaqP1cP7EyYj+b5RUQkdKkA+IHyymqmLy7glU+3EWbg/uE9uevCbsRE6r9HRER8Q1sYB7k9lndWF/HUh3kcKK/k6oGd+d3lqXRqHe10NBERCXIqAE3sxBvsdGjfGTN0NBdccQ2v3JLOgMQ2TscTEZEQEWgXAgpos7KyeGjsWKbFx7OuVypPxTQjcslLXNNqmzb+IiLSpFQAmkhpRRXjfvMQE9u0YUhMLBHGMCQmlklt2zJpwgSn44mISIhRAfCxareHN5cXkvHkYvbs282g6G+fxz8oOobc7dudCSciIiFLxwD40JLcYjLn5FBQXM7Qbu0oT0wKuBvs6AqAIiLBSQXAB/L3lZE5J4eleftJbhfDzJsGc0laB96KncT4sXfzODV/+a+tcPFYSSlTZkx3OrKIiIQYFYBGdOjocZ79KI+sFTuIiQxnwsg+3Dw0mchmNTMtusGOiIj4CxWARnC8umae/7lF+biOuxkzJIkHLu5FXGzk956rG+yIiIg/UAFoAGstH2bvY/LcHAoPuriwVwITRvahZ4eWTkcTERE5LRWAetq0u5TM2Tks33qQHu1b8PovzyYjtb3TsURERGpFBaCOisuOMW1BHm+vKaJNdASPX3Umo89Joll47c+o1K5/ERFxmgpALVlrmbF0Cy98XMBxt4fbhqVw34970jomwuloIiIidaYCUEvGGLJ3H2FYj3gevqIPKfGxP/wiERERP6UCUAdPXzfgm1P6REREApm2ZnWgjb+IiAQLbdFERERCkAqAiIhICFIBEBERCUEqACIiIiFIBUBERCQEqQCIiIiEIBUAERGREKQCICIiEoJUAEREREKQCoCIiEgIMtZapzP4lDFmP7Dd6RxBKh444HSIIKcx9j2Nse9pjH3vxDHuaq1N+KEXBH0BEN8xxqy21qY7nSOYaYx9T2Psexpj36vPGGsKQEREJASpAIiIiIQgFQBpiJlOBwgBGmPf0xj7nsbY9+o8xjoGQEREJARpD4CIiEgIUgGQOjPGFBpjNhhj1hljVjudJxgYY141xhQbYzaesCzOGLPQGJPv/djWyYzB4BTj/AdjzC7v+rzOGHOFkxkDmTEm0Riz2BiTY4zZZIy537tc63IjOc0Y13k91hSA1JkxphBIt9bqvN5GYoy5ACgH3rTW9vUumwocstZOMcaMB9paax9yMmegO8U4/wEot9Y+5WS2YGCM6QR0stauNca0BNYAPwNuRetyozjNGF9HHddj7QEQ8QPW2mXAoe8svgp4w/v4DWp+yKUBTjHO0kistXustWu9j8uAHKAzWpcbzWnGuM5UAKQ+LPChMWaNMeZOp8MEsQ7W2j1Q80MPtHc4TzC71xjzlXeKQLunG4ExJhkYCKxA67JPfGeMoY7rsQqA1Mcwa+0gYARwj3e3qkigmgF0BwYAe4BpzsYJfMaYFsC7wAPW2iNO5wlGJxnjOq/HKgBSZ9ba3d6PxcC/gHOcTRS09nnn+76e9yt2OE9Qstbus9a6rbUe4GW0PjeIMSaCmg1TlrX2Pe9ircuN6GRjXJ/1WAVA6sQYE+s98ARjTCxwKbDx9K+SevoAuMX7+BbgfQezBK2vN0xeV6P1ud6MMQb4C5BjrX36hC9pXW4kpxrj+qzHOgtA6sQY042av/oBmgGzrLWTHIwUFIwxbwEZ1NzRax/we+DfwNtAErADGGWt1QFsDXCKcc6gZrepBQqBu76er5a6McacB3wCbAA83sWPUDNHrXW5EZxmjG+kjuuxCoCIiEgI0hSAiIhICFIBEBERCUEqACIiIiFIBUBERCQEqQCIiIiEIBUAEak1Y0wbY8zd3sdnGGP+6XQmEakfnQYoIrXmvfb47K/vpCcigauZ0wFEJKBMAbobY9YB+UAfa21fY8yt1NzhLRzoS811yCOBm4BK4Apr7SFjTHfgBSABcAF3WGs3N/0/Q0Q0BSAidTEe2GKtHQCM+87X+gKjqbkG+STAZa0dCCwHbvY+ZyZwn7V2MPBbYHqTpBaR79EeABFpLIu99ycvM8aUAv/xLt8A9PPevexc4J2ay5kDENX0MUUEVABEpPFUnvDYc8LnHmp+14QBJd69ByLiME0BiEhdlAEt6/NC7z3LtxljRkHNXc2MMf0bM5yI1J4KgIjUmrX2IPCZMWYj8GQ93mIMcJsxZj2wCbiqMfOJSO3pNEAREZEQpD0AIiIiIUgFQEREJASpAIiIiIQgFQAREZEQpAIgIiISglQAREREQpAKgIiISAhSARAREQlB/w99kH+NFQuaOgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Open figure window for plotting data with linear fit\n", "fig = plt.figure(1, figsize=(8, 8))\n", "gs = gridspec.GridSpec(2, 1, height_ratios=[2.5, 6])\n", "\n", "# Bottom plot: data and fit\n", "ax1 = fig.add_subplot(gs[1])\n", "\n", "# Plot data with error bars\n", "ax1.errorbar(x, y, yerr=sigmay, ecolor='k', mec='k', fmt='oC3', ms=6)\n", "\n", "# Plot fit (behind data)\n", "endx = 0.05 * (x.max() - x.min())\n", "xFit = np.array([x.min() - endx, x.max() + endx])\n", "yFit = fit[0] * xFit + fit[1]\n", "ax1.plot(xFit, yFit, '-', zorder=-1)\n", "\n", "# Print out results of fit on plot\n", "ax1.text(0.05, 0.9, # slope of fit\n", " u'slope = {0:0.1f} \\xb1 {1:0.1f}'.format(fit[0], dfit[0]),\n", " ha='left', va='center', transform=ax1.transAxes)\n", "ax1.text(0.05, 0.83, # y-intercept of fit\n", " u'y-intercept = {0:0.1f} \\xb1 {1:0.1f}'.format(fit[1], dfit[1]),\n", " ha='left', va='center', transform=ax1.transAxes)\n", "ax1.text(0.05, 0.76, # reduced chi-squared of fit\n", " 'redchisq = {0:0.2f}'.format(info.rchisq),\n", " ha='left', va='center', transform=ax1.transAxes)\n", "ax1.text(0.05, 0.69, # correlation coefficient of fitted slope & y-intercept\n", " 'rcov = {0:0.2f}'.format(cvm[0, 1] / (dfit[0] * dfit[1])),\n", " ha='left', va='center', transform=ax1.transAxes)\n", "\n", "# Label axes\n", "ax1.set_xlabel('time')\n", "ax1.set_ylabel('velocity')\n", "\n", "# Top plot: residuals\n", "ax2 = fig.add_subplot(gs[0])\n", "ax2.axhline(color='gray', lw=0.5, zorder=-1)\n", "ax2.errorbar(x, info.resids, yerr=sigmay, ecolor='k', mec='k', fmt='oC3', ms=6)\n", "ax2.set_ylabel('residuals')\n", "ax2.set_ylim(-100, 150)\n", "ax2.set_yticks((-100, 0, 100))\n", "\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the same $(x,y)$ data set but with a single value of $\\sigma$ for the entire data set." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "slope = 22.4 ± 1.7\n", "y-intercept = -72 ± 24\n", "redchisq = 1.34\n" ] } ], "source": [ "sigmay0 = 34.9\n", "fit, cvm, info = linfit(x, y, sigmay=sigmay0, relsigma=False, return_all=True)\n", "dfit = [np.sqrt(cvm[i,i]) for i in range(2)]\n", "print(u\"slope = {0:0.1f} \\xb1 {1:0.1f}\".format(fit[0], dfit[0]))\n", "print(u\"y-intercept = {0:0.0f} \\xb1 {1:0.0f}\".format(fit[1], dfit[1]))\n", "print(u\"redchisq = {0:0.2f}\".format(info.rchisq))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fit a huge data set." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create data set with 100000 data points." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def randomData(xmax, npts):\n", " x = np.random.uniform(-xmax, xmax, npts)\n", " scale = np.sqrt(xmax)\n", " a, b = scale * (np.random.rand(2)-0.5)\n", " y = a*x + b + a * scale * np.random.randn(npts)\n", " dy = a * scale * (1.0 + np.random.rand(npts))\n", " return x, y, dy\n", "\n", "npts = 100000\n", "x, y, dy = randomData(100., npts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit straight line to the data." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "fit, cvm = linfit(x, y)\n", "slope, yint = fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the data together with the fit." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeXiU5dWH72f27OsESAIkQUiCVbRitdCKOwnaWm1VQG2tIAoqixZcqHXfC4RFURDcWax+1aokuFcUd6moJBEJIGHLTJJJZpLM+j7fH5MZkhAgKDMJ8NzXlWtm3nmXM6P83jPnOYuQUqJQKBSKowtddxugUCgUiuijxF+hUCiOQpT4KxQKxVGIEn+FQqE4ClHir1AoFEchSvwVCoXiKOSQiL8QYqkQokYI8W2bbalCiLeEEBtbH1PavHerEOIHIUSlEGLkobBBoVAoFF3nUHn+TwNFHbbdArwjpRwIvNP6GiHEYGA0cGzrMY8JIfSHyA6FQqFQdIFDIv5Syg+Aug6bLwCeaX3+DPCHNttXSCk9UsrNwA/Arw6FHQqFQqHoGoYInruXlHIngJRypxAio3V7FvBJm/2qW7ftl/T0dJmTk3PIjVQoFIojmS+//NIupbR23B5J8d8XopNtnfaYEEJMACYA9OvXjy+++CKSdikUCsURhxBia2fbI5nts1sI0af14n2Amtbt1UDfNvtlAzs6O4GUcpGUcqiUcqjVuteNS6FQKBQ/kUiK/3+Av7Q+/wvwapvto4UQZiFELjAQ+CyCdigUCoWiA4ck7COEWA6cDqQLIaqBO4AHgReFEOOAH4GLAaSU3wkhXgQ2AH7gOill4FDYoVAoFIqucUjEX0o5Zh9vnbWP/e8D7jsU11YoFArFwaMqfBUKheIoRIm/QqFQHIUo8VcoFIqjECX+CoVCEWGklLjLy+lJY3OV+CsUCkWEkFLSsmEDjjfeYNs11+IuL+9uk8J0R4WvQqFQHBV4KiqonjgJzetF6DprbtB9KPFXKBSKCGEuKCB74WNomoZv61bMBQXdbVIYFfZRKBSKCCGEwFJYiG/LFnbfd1+PCvso8VcoFIqfQcfF3I6vPRUV1Nz/AFq9A++WLd1oaXtU2EehUCgOAiklnooKzAUFCCGCcf0bJpM1by4Ans2bsc8pCb/WNA3rrbeCgMTi4u40vR1K/BUKheIgaCv2QghM+flkz5+HlJKtE66BlhZ63XkHns2bqXnwIaTfj+b3kz5+XHeb3g4V9lEoFIouIqVEShn26qtvmIynogIAU34+cSefjGxqomXdOmwPP0L69OmYf3UyOJ3YS+bS8MYb3Wl+O5T4KxQKBXty8ls2bNhnMZanooLtk6cghMBcUID1xmlIKam+YTK1Tz5J80cfQWwsIj0d6fXS8r91tKx+M3QB3D/+GMVPtH+U+CsUCgWt4ZxrJ1I9cVLYm++IuaCA7PnzMBcU4K2sxDZ7DkII0qdNpXbRYrTGRgwDBuBYtJiAy0Xjv15qd7zzvfej8Em6hor5KxQKBa3C/vjC8PMDoT/mGCxDjscwcGC7Xwr+b77Z5zFxp4/4+YYeIpT4KxQKBcGc/JjBg4E96Zqhm4CnogJTfj6NpaXYHnoY680zaFm3Ducbq/A3NJD4+99DoHUmlcUCCQlgs+11Dcupp0bt8xwIFfZRKBRHPAfbWC2U0eOpqMBTUcG262+gbskSah58CM3jYefd90BqKsYhx9Py0Vp23347eDzBg/3+ToUfwH7LLYfqI/1sIir+Qoh8IcT/2vw1CiGmCiHuFEJsb7N9VCTtUCgURzdtxXxftL1BmPLzSZ82FU3Tgs+nTsG+aDFxZ5+FdeZM8PtxzJuPb+MPICV4vMFHCIr/PtDi4w/1R/vJRFT8pZSVUsoTpJQnACcBzcC/W9+eE3pPSrkqknYoFIqjm7YLtZ0hpaSxtJRt198Q9vZr7ruf7ZOuw1NRgeurdUink8blKwBJ3GmngRDo8nIPyg7T4GMPwac5NEQz7HMWsElKuTWK11QoFIpwjx0h9nTWbOvpeyoqqJk1m7gRI3C8/z6u778n0NxM8rSpOP7v/3C98krwIIMB744dNJWWgk6gffvdQdnh/e9/u7yvpmmsXLmSiy++GE3TDuo6XSGa4j8aWN7m9fVCiPVCiKVCiJTODhBCTBBCfCGE+MK2jxiaQqFQ/BRCsfyGVavQNI3kSy+hYfly6ufOw37HneB241j6FI7nX0Cflxc8SKfDt2t38HngJwhyTMwBd5FSUlZWxtChQxk9ejQVFRXs2rXr4K91AEQ0JssIIUzADuBYKeVuIUQvwA5I4B6gj5Tyqv2dY+jQofKLL76IuK0KheLIo2M/npDX7960iZp770MDes28De+OHQT8ftxffIH308/2xPFDCAE63Z7Mnp9AYcW+O3uuXbuWW2+9lQ8++IDc3FzuvvtuxowZg16v/8nXE0J8KaUc2nF7tDz/YuArKeVuACnlbillQEqpAYuBX0XJDoVCcRQRFvny8nYLvu7ycqonTsK7fTtaYyN4PNjuvof6J5fg/Wod3k8+3Vv4gyf8WcJPbGynm9evX8/vfvc7hg8fTmVlJQsWLKCiooLLL7/8Zwn//oiW+I+hTchHCNGnzXsXAt9GyQ6FQnEUEcryAciaN5dAIIDjjTcIBAJIrxdDZia6xERiTh+B5nKBy4X747XBg0UEJm81N7d7WVVVxeWXX84JJ5zAmjVruP/++9m0aRPXXXcdJpPp0F+/DREv8hJCxALnANe02fywEOIEgmGfLR3eUygUioMmFNox5efjrazEXFCAuaCAzLkluKuqQELNPfeguVyk3nA9MhAI5v1LScs77wbFXsqgKgVPGDFbd+7cyb333suiRYswGo3MmDGDGTNmkJqaGrFrdiTi4i+lbAbSOmy7ItLXVSgURxchL9964zRss+eQPX9ecIrW1q3smnFzsPJW08BoRAqB1tRE0/v/RfN4wOeLio0NgQC33XYbJSUl+Hw+xo8fz+23305mZmZUrt8W1d5BoVAcEZjy87HeOI24c89FahrNGzfS8sMPGHNySBozGueqUvD50Jqa8O3cheXUU2guLY2KbS2axvP19Sypq8X54IOMGTOGu+66i2OOOSYq1+8MJf4KheKwJhTukVJSM2s2ydXV1C1ZitbQECzESkxAmMxYb7uVpnXrcC1bjmvlyqjY5pOSlx0OFtbWYgv4GREXx9yPPmLIkCFRuf7+UOKvUCh6JB3TMzt7311ejmfzZmyz55A1by4poy+l7rnn0FpaAIi/9BIM1gwCtXaav/wK14oVUbFdk5I3nI0ssNvZ5vPxy5gYZqdnclJsLIU9QPhBib9CoeihhNIxsxc+tle3zRDVEychfT4yZt4GQN1zz2M58yyaXnwRAK/djmtFGy9fr/95qZoHQErJf5uamGu3UenxkG82szArm9Pi4jq9gXUnqqunQqHoubTJuAn337nmWqonTgLAOmM66PWYc4M9dqTbHRT+1nYI3rffaX++CAr/F83NXLHtRyZtr6ZZ0/hnn0xe7p/DiPj4PcLfg24AyvNXKBQ9irY9d7IfX4ilsBAIZvPYZs8h6bLL0IngAq8pPx9fdTXNGzfi27EDqWlRF9hyt5sSu401TU1kGAzc0asXFyUlY+zMjih0VOgqSvwVCkWPIjROESFIn/43PFVVGHNyAEi65GLsS5aC04kk2LDNPqek/QkMBujfD7ZGdl7uFq+XBXYbq5xOEnU6bkq3MjYlhRjdfgIqPcjzj0pvn0OB6u2jUBwdhAapuz78kNrFT4LLFWyIptMhjEZSrvwLjqVPoWkaQkp0vXsT2LQpavbt9vl4rLaW/2twYBKCP6ek8tfUVBK70obBZKJw/deRN7IN3d3bR6FQKNrRcbpW29fezZupX7IUnU4X9OR1OmhuRrrdCATpM2eCx4Nsaoqa8DsCAf5ZU0PR5ir+3eDg0uRkyvIGMMVq7ZrwA0S4ZcPBoMI+CoUiqrTNy98+eQrZ8+dhys/HvngxdU8/Q+qVf6Hu6aeRjY3EDB9Oy0cfIfr0QTY1gc9H7dy56AcMCFblpqVBbW1E7W3SNJ6rr+OpujpcmsbvEhO5Pi2d7J8i5D7voTfwJ6LEX6FQRJVQG4b0aVPJmjcXU34+dUuWUFsyFwDHM8/S67bb8O3YgfHkX+Hdto3Aj63x+9ZsncAPPwRfR1D4vZrGiw0OnqitpTYQ4Mz4eKakWxloNv/0k/agSV5K/BUKRVQxFxRgvXEaNbNm03fBfLyVldSvWEna1CkYMzMx5ebi27KFuqeeRi54dE/fHbN5z5D0CBKQktcbG1lQa2e7z8evYmKZb7VyQhcGsRyQr6Mb798fSvwVCkVE6VipK4QgoagIKSWapmEuKCDjphuJHzkS3/ffo2kau++7H+lytR+GHmHhl1LyrsvFXLuNH7xeBpvN3Jndl2GxsYeuQGsf/fy7AyX+CoXikNBZS2UhRDjMkzVvLkIITPn5OMvK2P3gQ+B2Y/3H7dQ+8k80TcP20MMkXTYWze1uL/wR5rPmJubYbHztdpNjNDE7M5Nz4xPQHeLUTNO55xzS8/0clPgrFIpDQkjkky+9hPoVK+m7YD6WwkLMBQVkz5+HpmlUT5xE8mVjqVu+gtihQ3GVleFe9z8QAm91NQGXi7r5C6JWDPWd202JzcZHzU30Nhi4u1dv/pCUhCFC+fixWdkROe9PQaV6KhSKQ0Ioll+/YiUZN92IuaAACBZiWQoLEUKgeb3UPfU0SZdcTNO774LJRMOrrxJz5pnULXwcWlqCi7raTxiOfhBs9nqYun07F2/dwnceNzOsGZTm5vGn5OSICT+AXt9zJDcak7y2AE4gAPillEOFEKnASiCH4CSvS6SU9ZG2RaFQ/HQO1GUzFMsHSCgqCg9K91RUYBw0CE9VFUl/voLG557D9d77SI8HUlLA4cD50ktRGaiy0+fjsVo7rzQ0YBY6Jqal8deUVOIjNCe3Lbq8PFLGjYv4dbpKtMI+Z0gp7W1e3wK8I6V8UAhxS+vrm6Nki0Kh+AmEwjqhCVkh2t4UvJWV2GbPwZSTExb/7ZOnYD7+OFxlq4MFW1IS+P774MH19cGWBxEW/nq/n0V1tSx3OJDA2JQUJqSmkWaIXuRbq6rCVVZGyu9+F7Vr7o/u+g1yAfBM6/NngD90kx0KhaKLhGL3oXBOCHd5OduuuZamb7/FuWYNfUrmALDt+hvwbN5M4p/+hGtV68Qsrzco9G73nhNEML7fpAV41G7n3M1VPFdfz3kJiazKzePWjF7REf7QL4oMKwC+HTsif80uEo3bngTeFEJI4Akp5SKgl5RyJ4CUcqcQIqOzA4UQE4AJAP369YuCqQrF0cOBwjid7dfW42+HEDhffRXH8y+QYrORfvPNxJ9xOrsffAjTiScE9wl13IzCYq5H01jpcPBEXS31gQDnxMczOd3KgJ9ToNUV0tPB3hrk0OkwnjwU3/pvSB49Gr1eT9r48ZG9/kEQDfEfLqXc0SrwbwkhKrp6YOuNYhEEG7tFykCF4mhkX2Gc/e1nLijAXV6Opml4t2zBFOq2OXYM5lNOxVtbR/0rr+K122kqLQODAfebb+05WYSF3y8l/2ls4FG7nZ1+P6fGxjI13crxh6JAqyu0Cr9+4DEENv6A75NPMQ8ZgmPBo2T+8xH0UVhb6CpR7eophLgTcAFXA6e3ev19gPellPn7O1Z19VQoDi0/xfP3VFSw7dqJaO4WNKcLXVwc6HRoDQ2IuDhkKKwD6Pr2Rdu2LWqf5S2Xk3l2O1VeL8dZLExNt/LruLiIX9t41plou3cT+Pa74Aa9Pvgrx2AAnw9daippf72S1HHjgo3qosy+unpGVPyFEHGATkrpbH3+FnA3cBZQ22bBN1VKOWN/51Lir1B0P6HOm5qm4amqwrdjByIjg9rb/4GIjUWG4vlRCu8ArG0dm/iN202eycSUdCtnt52eFWkMBvD7MQ0ZgjE/PzhJLCaGhAt+j/O11+l95x0kn39+t41x3Jf4Rzrs0wv4d+uHNgDLpJRlQojPgReFEOOAH4GLI2yHQqH4CXRWtRszeDCaptH4yis4nn8B0/HHQyCA+fjjID4e9+o3oyL861tamGO38WlzM30MBu7r3ZvfJyahj7bI+v1gMKDt2kniFZfT8tZbJIwqxvXmW+gtZiwDBvS4+b0QYfGXUlYBe42ql1LWEvT+FQpFBNlXy4X97WscNAjX6tUkFBUF4/0TJ5F2003U3H03iVf9FYNej6FPJo4XXgCTCe/mzaDX4/78i2A2T4T5weNhnt3G2y4XqXo9t2ZkcGlSMqZohlRyc2DzlnDqKrGxJF92GQnFxQidjt3/nEXGLTdjzs3d73pKd6LaOygURzChxVrrjdOwzZ6z38Xd0L5Jl1yMvWQufaTE0K8fAacT15dfIF0uGubNByDpsrFgau2yKWWwKjeCw9EBtvt8PGq385/GBmJ0Oq5PS+cvqSnE6SK8iGqxtEtNNRx/HLK2jrixY5AInKtXk1xchGPliyT89rckFheHF8JD6yQHWlfpDpT4KxRHMKHcfFN+PubcXMwFBftc6DXl55M+bSpaIIA+ORlj//44Xv4/ZHMz/rq64E5GI5aTTybhggvwaxpNy1dEvECr1u/nidpaVjY4EMCfU1K5OjWVlGgVaLWtSQAshYMx9+lN7ROLwOMheewYmt7/b7ilReg7rb52ItabZ2CfU3LAjKruQIm/QnGY0pVsnVBfHQBLYSFSShpLS8O99EOeaajT5q77H0B6PKRdPT6470svAeBpHZUo+vXFvXYt29etQ6amRvTzOQMBnq6v4+m6OrxScmFSEhPT0uljNEb0uvtECDAYcL3xBk1CQEsLIi6OhAsuIObEEzH277/X/qacnE4L43oCSvwVisOUrubpQ/vRibbZc7DeOC2cubN98pTwcJW4X/0KV2kp9YufxLu7Zo9Xv2MnGI3ITVXB87W0wPbtEflcbk1juaOexXV1OAIBRiYkMDk9nVxThAu0OsF09tnEFBbg/2ETTZ9+Stpf/kzDsuWk/+0mhBAYc3Lwbd2K7eFHkJpGr1tvIbG4GEthIX0fX9gjwz0hek6LOYVCcVDsq91CZ4RuFEDwmNxctk+eAkDWvLkY+vUn6ZKLaf7002DevssVbLYWwu2OeHjHLyX/cjgo3lzFIzYbvzBb+Ff/HOZkZkVP+FsXjXW9ewcXs995B6PRiOerr0geVUzSlVcSf87ZmHJzSTrvPPR6ffBmOmM6GbfcjG32HDwVFe06mfZUolrk9XNQef4KxU+nY4hI0zScZWXEnXsu9UuXUvfU02heb3AB1+9H168fWmhuboTRpGS108l8u50tPi9DLBamWa38KjbyBVodEXFxmIcMwb12LZZhv8aUk0PShRdSfdU4pNNJfHExrtJSdPHx9Hv2GSyFheHvFeiRi7vdleevUCh6AG1j/wDeykpqZs3G8vbbwW6bJlMwayeUsx4F4ZdS8lHrBK1yj4eBJjMLsrI4Iy6KBVqdYBk6FPeXX+Je+zHutR8T+8tfkjZ+HPa584g9/XQMaak0lK0G9nyvXa2W7kko8VcojkBC8fx9/bLXNInpuOOC3TZDueqhsYkDBkBlZUTt+19LC3NsNXze0kKW0ciDvftwXmJi9Au02qA/9lgCP2zEsXQpGI0k/PGPuFatwpybS9KoUZj69g3OKfj970j54x/bZU6FWlf3xKyefaHCPgrFEUjLhg1sG381MuBH8wVFXWc2Iwi22ZVeL7Kpac8BMTGQkAA1NRG1a6PHQ4ndxnsuF2l6PdempXNxcjKm7hB9kwlDfj7+bduwDB6Me906aGnBMnw4vvJyshYtQq/X7debd5eXt5tP3BM9fxX2USgOQ9qGE6BrMWUpJZ7Nm5FSonl9wdGIQOKFfyDxwgvxVFWx+867wGhEDDwGWVER3Kd1v0iwzetlQa2d1xsbidfpmJKezhUpqcRGsyo3IQGczuBzkwl8PtKuuAK9Qc/uBx9CmM3EnXE6TZ99jk5KdDqxz6E1oe+/7aJ7TxP9A6HEX6HowbRN5wTapXZ2tojbWFoaTufMmHkbAX8A1/vvISQ43/8vAF6/H5qbAZAbyiNqv83v54laOy86HOiF4KrUVMalppHcHa2NQ8JvNBJ/0UW4VqxA6ESnFbnAXuGbzlJrO66lHE6osI9C0QPprCcPBMMMQDjLJNS6IbG4mIZVq9j5t+noEhLI+MftICW7bv9HuELVPGQInq+/jor9jYEAS+vqeK6+Dp+U/DEpmYnpaWQYolCg1dplM4wQJI4Zgy49DceixWTcfRfJ552Hs6wMU04OMYMHt5s3vC8v/nBc1IV9h31Unr9C0QPouEAbEnZvZWU4Xzz0t33ylLAIhYqzHG+8gXdbNcTGIgwGhBDYHniwXWsCz7ffQmwsRLBCtkXTeLK2lnOrNrGorpYz4xN4LTePO3r3jo7wQ3vhBzAaSbrwDxgMBvSJicQMGIDv++8x5+ayY8rUsKcf+s5DrzsSiumHFngPd5Tnr1D0AFo2bKD62olkP75wT8vk0tJ2nikQzs9PKCpCCEHzd98FWyu//gY4HCSOHUvyRRcC0LhmDY75C4KDRSKMT0pebnCw0F6LLeDntLg4pqRbKbRYIn7ttpiGDcP72WegE2AyY8zLI3bEaRj0emrnLyB9ymSMWVnY55TstUjbFc8+tMB7WGX1qAVfhaKHIwSezZuxFBbirazE9vAjICV9n3gcS2EhmqZRt2QJtcuWI6XE7/dju/0f4PWiGzAAzeGg8aWXaNm0Cd9nn0Wlp74mJaVOJ/PsNrb5fPwyJoZZ6ZkMjY2N+LU7ou/Xj/jzzqNu7drwVDHf+vU0fPMNuuRk0qdMxpCZSc2s2fT62017VeB2JX5/MFXVPR0l/gpFD8CUn0/CuecEC6/y8jAXFJD12KN4Nm9G0zSklDjLyrCVzCWhaCQ1Dz1EoK4+GOIQAq218RpeL75PP424vVJKPmhqosRuo9LjId9sZmFWNqfFxUU3Hh6K76emEqiupnnVKgCMxxyDd/16EsaMwdArA3NWFqbcXHZcdz0ZM6aTWFz8k+w8nBd4OxLRmL8Qoq8Q4j0hRLkQ4jshxJTW7XcKIbYLIf7X+jcqknYoFD0dZ1kZ9S8sI2X0peGwgxCCmvvuZ9v4q3GXl5NQVET6lMk4P/4EkZwSFn5y+h/4AoeQL5ubuWLbj0zcXk2zpvFIn0xe7p/DiGiMTmz7i0KnC/baNxigrg50OjyVlcSPKsa7dSsiLo7YX55I4zPPUvPAA3i3bAEpMefmHlYLtpEi0p6/H7hJSvmVECIB+FII8Vbre3OklP+M8PUVisMCU04O+uRk4oYP3ysUAeDetAn3pio0TYPmZvzffx/cQcrgRKkoUOF2U2K38UFTE1a9gX/06sUfk5IxRklIzeecg3Q68X7ySXCDphH329/Q9MEHwRuh2UzvW2/BlJOD+8uvsN7+92AKpxBInw9TTg59n3i805DN4ZrJ83OI9BjHncDO1udOIUQ5kBXJayoUPZl9FW1ZCgvp9feZ4ZBCKBEja9ETNPz73+y6865gbr6UmM4+G39VFVpVVdD7jfCC7lavlwV2O284G0nU6bgx3cplKSnERLNASwg8a9dCUxOYzegGHoO2ZSstn3xK7zvvAgHm3FxiBg8GCLdTBsiYeRs1Dz2MTqc74BSzw2kh9+cStZi/ECIHOBH4FBgOXC+E+DPwBcFfB/XRskWh6A46DlKRUlI9cRLZCx8DYPeDD6FpGubcXLxbtmCfU0L8GafTsHwFmEzEjx6N61//wvv223tOGkHhr/H7WGiv5eUGB0YhmJCaxlWpqSRGq0ArPh7jqafie/vt8KjI2KKR6FJScC1fQVxREe4vv8Scl4tO174NQ1sBTxo1KryOsi+OpIXcrhKVW7cQIh54GZgqpWwEFgIDgBMI/jKYtY/jJgghvhBCfGGz2aJhqkIRMTwVFcHK25tuxJSfj7uqCs3rxV1VFfT0fT5233U3P/75L+y84070gwup//cr6EeMALebpq+/3juHPQI4AgFm2Wooqqri5QYHlyQnU5Y3gKlWa3SEPyTcLheWjAyIsYS7jja/+RbS4QAhiD/zTPotegIhBNuuvyFc3dyRrvTWPxz67x9qIp7nL4QwAq8Dq6WUszt5Pwd4XUr5i/2dR+X5Kw532oZ8PBUVbLvmWjSPB3Q6shcvoumjj6hd+PheM2P3qliNEM2axnP1dSytq8OlaZyfmMj1aen0NZkifu12fXd0umAhmhAkXHghrlWrkD4fNDdjGT4M9w+bsF5+GanjxqHT6dr9osq46cafnMlzpNItef4i+F9gCVDeVviFEH1a1wMALgS+jaQdCkV30lb0Q+EIc0EB2Y8vxL1pEzUPPIBvyxbqn3seAoFgBk9bpyzCwu+VkpccDhbW2qkNBDgjLp4p1nQGmaNYoOVyoR84kMDGjWCx0OvOO3F//T8aVqwg9YbJBGpqaPjPf/CVV9Bn5m0kjRoVFnghgv15pJTUPPRwuDBOsX8iHfYZDlwBnNkhrfNhIcQ3Qoj1wBnAtAjboVB0G521DQhVlvp27EALaPj9fgyZmcFRiZ39Go9AS4aAlPynoYHzN1dxb81u8kwmXujXj0ezs6Mr/AYD5uOPDwo/QMCP0Akay1YjLDH4a3bTsHIlyX+4APT6TlM1hRCYc3OjUth2pKDaOygUh4D9pQqGWjLEjxyJ7/vvMeXn07KhnNolT9JUtjoY5pAyKm0YQra+1+Rirs3ORq+HQrOZaVYrw2OjVKAlBCQnQX1r7H70aFzLlyMGDEBu2gR6PdnLl+H+5BPqnnoaDAbS/nwFSVdeScPTT5Ny1VXo9fq9vvOjMV2zK6j2DgpFBAgJjqZpbJ84iezHF7ab6yqECLZqmD0Hz7ZtOJ55luS//Jn6pU+hNTQETxIIRM3ez5ubmWOz8T93C/2NRmb1yWRkQgK6aIqllJhOOgnvu++Rdt99aDuDEWC5dWvw/UCAlrUf43j2Way33UrMgAFYCgtpLC3FVjIXY1YWSaNG7ZWeeSRV30YD1dVToThI2nbgDId0Nm9G+v3tt7V2f9Q0jeIkh7kAACAASURBVMSL/0Tt4icJ1NdTu2QpxijHpDe43UzYto2/bPuRnX4fd/XqzX9y8yhOTIyu8Lfifecd0qdOIf2C3+O3BaeHxV50IWlTp9LrkYcxZGUScDjQ6XThxnYJRUVkPvIwCUVFQNfSMw80zvJoRnn+CsVBEMossc2eQ9a8uUgpSZ82FUO/YIsFT1UV5rw8MueWIKUMduu8egIBjye4cBsTQ8yvT6XlzbcOcKVDwxavl3l2G2VOJ0k6HdOtVsYkp2CJZoGWXo8uNwetejvJE64mYLfTXLaauOHD8VZW4nrrbXTx8SSefDLJ558fHkyj1+vDQg+g0+lIGrWnE0xXPP2jsXirqyjxVygOAk9FRTilUEpJ9firESYT1hnTAai5/wF0RiPpM6az6557sRw7mIDbDVoAPB5ETAwt77wbjPNHMNyzy+fjsVo7/25owCQEE9PSuDIllYTumKCl16P9sKn1qZ6MmTNxnnQSAMZBg+h96y0gJbY5JcQccwyWwsK9hP6ncjQWb3UVJf4KxQFou5AYEpPQdoxGrDOmhz1UXXY2Ta+9hhbQoKEB99qPwWAg9o8X0fzaa0iPN6Ki7wgEWFRbyzJHPRIYk5zCNWlppBmi/E+9tTYhrrgIqUma33wTzGYcz7+ATqcLprUKQcbNM8K99fsumH/IRVqtA+wble2jUByAtgM8zAUF4YKipEsuxpiZiT8QwP3fD3B99CHxw3+Dq7QUy8hz8TkcBKo2g80WFEOIWM5+kxbgmbp6nqqvo0XT+H1iEtelp5MVwaldHTGccgrGlGRaylYHXw8ciLZ7N5rTSczwYaRNnhzM4Fm+gl5/uwlzbi7mgoLwmEqVoRMZ9pXto8RfcVTT2azcjiK0V2Xu9TcQf8bpOF5YFmw74PHs2dliCebqRymDx6tprGxw8ERtLXWBAGfHxzM53coxZnPkL962qZzZTL8XnkfTNOwlJbg/+zx4o4uNJf70EbjKVpNy2Vic776nqnCjjJrhqzgqOVC2R2hB0FlWts/5rUIITPn5NJaWYhw0iKx5c9GlWxEWC3g8GE45JSiEEGzNEAXhD0jJvxscFG+u4oGaGgaazSzv1595WdnREX5oV5cQd+aZAOwYfzXujz8BIYgtKqLvM8+Q+cgjWKdNxfnue1hvnBZss6zodpT4K45oDjSUOxTDTygqCod1OrthOMvK2DF9Bq7Vq/Fu3kxdSUlQ8HU6At9+G9UCrbedTi7YspmZu3aRZjDwZHZfnurbjyExMVGxAYDWPvkACEHT2rW41nyIBGKLRoLfj/vTTwls+xGdTkfa+PHBmH5ubngAfcfPpVIyo4ta8FUc0Rwo2yO0INhWdDpLD4wfOZK0H7fRsnXrng4CKSnQ1IRsaor0xwDg49axid+43eSZTJRkZnJOfEL3hE+2bAGLBWP+IHybqhB+P3ULF4LPR8IZZ2BKS8N8/BBss+dgzs3FUlgY/p47+++hUjKjjxJ/hYLWRd1rJ5K18DGEEGTOLUHTNJq+/Rbfli3B4elPPIFsaUHEx4PJhNy+PSq2fdPSwhy7jU+am+ltMHBv7978PjEJQ7RFX6cLevuBAMZTT0H7fiOBrT8idALpC4DPR/JlYxE6HY7lK+hz4ol7Cf2+sm9USmb0UeKvOKI5KI9SCJo++gjHi//CfNxxNL3/PjqjEa2xEczmYHz/F7/AUJCP+6WXg4cMHozcsCEitm/yeJhnt/GWy0WKXs8t1gwuTU7GHMUCLfOwYST8/ncEdu0GqaHvk0ntXXcFu2ve/ncM/fvTsnYtll//msC2bSQWFwOt6yQ5OV3O4lEpmdFHib/iiKYzjzIUXw69H8rysc6Yzq77H8BSWEBTaSkAlhEjaH7vvXBM3//tt/i/3dOBPBLCv8Pn41G7nVcbG4jR6bguLZ0rU1OI00W3QEvXu3d4dKK3qgrpdJI4dgzS7Qa9HoTA/+OP2OfNJ7NvX5LPOy98rDk3l23XXEvGLTeH2y+rxms9CyX+iiOajh6llBLH629Qc9dd6GJjybjl5nDFrrF/Dvh8tHzzLZZzz0FraEQLBNqnckaQOr+fRXW1LHc4EMAVKSlMSE0jJUoFWrrBg9EqK4PZSunpaLt2oevbF8/XX4eHq1hOOAFTr15oUmKfU0Lm3JJ2/XbaEQhge/gRLHl54WZ3Kq7fc1DirzhqCAr/6+z6+9/B4yX2tN8SP3IkALbZc0idfAPS4wFNw/3ue1GZngXgCgR4qr6OZ+rqcUuNC5OSmJSWTp9oFWjpdGAwoG3dGgxv+XxgtwNgGlyIp64OKQRp119HzIABbJ81G6lp9Lr1lnArhpAn39a77/vkYoDwry4V1+9ZKPFXHNG0DfEA7L7rbvB4MQ05npb/fY23shJj//6kTL4Bxyuvtvfyk5PB5YrYTcCjaSx3OFhUV4sjEODc+AQmp6eTF608/RCaBl4vIEFvAIsFy8lDca/7H+6P1pJ69Xgcz7+AKTsbS2FheOB8Z95829cdp2mpuH7PQlX4Ko5o3OXlbLvmWhCC5MvGIq1Wau+8i7Q77kDbvYuWigrc770fFEC/v33VaoTwS8mrDQ08Wmtnl9/P8Ng4plqtHGuJ4vSsNujy8tCqqsBiRme2kHLVX2l48V/BKublK+jz0IMInQ7b7Dl7hWzUQJWeT48b5iKEKALmAnrgSSnlg91li+Lwpa3YAHs9Nw4ahHXGdJrXrcNeMpf4opHg81F7zz3Q0rL3CSMo/JqUvOlyMt9uZ7PXy/EWCw/06cMpsXERu2aX7KqqQpjNJF9zDY3LlhP/m9+QeNppGAcNIvaXvyShqCicvSOlRErZbn5u25uB8u4PH7pF/IUQeuBR4BygGvhcCPEfKWVkcuYURyxtwwwA266/gYybbsSUk8P2yVOC3mtpGdLtJm7kubg+/yI4TtHvx/NWdHrqSylZ29xMic3Gdx43A0wm5mdmcWZ8fOS947g4aGqCzEzYsQMS4jHm5uFbvx4xIA+5qSq4n8FAwm9/S9KIEeEBNN7KynY9eIQQasH2CKK7PP9fAT9IKasAhBArgAsAJf6KLhPyQrPmzQ17+9ZpU9l9731kLXqCuNNHUP/c8+h79ybgctH07nvo+/XF1ZrGGQ2+bmlhjs3GZy3NZBmNPNC7D+cnJqKPRkgkPR1cLkRMDEl/+hOOefOIHT6c5g/WBGfnnnIqlgkT8H69Huebb4YXbquvnYgM+BF6A32feDws9GrB9siiu8Q/C9jW5nU1cErHnYQQE4AJAP369YuOZYoeTyjUI6Vk++QpZM2bG17UlVISqKujdskSmj77PNh7Z9eu4IEeD4GNP0TFxo0eD3PtNt51uUjT67ktI4NLkpIxRXOCVmvGjjQYMGb2IW3aVBzPPIMwGLAMG4ZzxQpkfT1Nn3xCxm237rmZti7oCiG6VJ2rODzprsZunbk9e608SykXSSmHSimHWq3WKJil6OmExihuu/4GADLnluBcs4Zt469m27jxaJoGJiNNZauJP+ecPd02IdiLJ8JUe73csnMHf9iymc+am5mcnk5Z3gAuT0mNrvALAWYzhoEDwe/Hft/9reMRJWgaLRUVxA4bhqu0FOlwENi5k+0TJ7F90nXodDpijz02PDtXcWTSXZ5/NdC3zetsYEc32aLoYXSWQdIuXfOfs0g484xgXHrLFmrnlIDBgIiPR9u1C7w+AFyfftI+TbO+PmI22/1+Hq+186LDgV4I/pqSyvi0NJKjPTZRCJAy+OfxELDbSL78MixDhlA7dx69b/87ppwcvFu2YJs9h7RpUzFmZmLKzSV22DB8W7diHDQId3m5ytg5wuku8f8cGCiEyAW2A6OBsd1ki6KHEUrPDLUG8FRUUD1xEkhJ9uMLSR0zGtucEhpfe430mX8PDlTxejEXFiJ690b/q18R+Owz2Lwl4rY6AwGW1tXxbH0dXin5Y1IyE9PS6BXFCVohTKeeinfduna1CnGn/preM2cCEHPMMWFBjxk8GEteXnhATfUNk7HeOA37nBKEEJ2mdSqOLLpF/KWUfiHE9cBqgqmeS6WU33WHLYruRdM0nGVlJBQVoWsbFmnTGsBcUED2wsfCbZdTrroKKSW1zz6H6913WguUwP3xx7g//jgqdrdoGsvq61lcV0ujpjEqIYHr063kmExRuX5n+DZsQJhMSMBy+ggMKSmk/ulPwP5TMkMLuab8fMy5ueFHtbB7ZNNtef5SylXAqu66vqJnEBqSkgkkjRoVFvjsxYvCC46apuFcswYpwbFyJb3+dhNxw4cHbwBzSoJNxkLTs5KTweGImL0+Kfm/BgcLa2up8fv5bVwcU9KtDO6OAq3WEI9l+HAsJ/2SxNNOw7tlCzX334/vq3Wk3noLO6dOw3rjtL3GJnYMrYVuBB0fFUcuqr2DoltJKCois/URgnn7oQyeULzf9uCD1D/3fDA9sbiImvvuAwSpN8/YI/yhytwICb8mJWVOJ/PsNn70+TjREsM/+2QyNDY2ItfrEmYz8WecgWv1atzr1mHp14+kUaOw5OWFb6Lp06ZSM2t2eKBKCNVkTaHaOyi6nbZeqJQSZ1kZhn79qL52Igkjz8X11tuQkYF01KPVO4IxbaMRw8Bj8K//JuK2rWmdoFXh8ZBvNjMl3cqIuLjoL4YajRATQ9L55xHQJK4VK7A+9CC+b77B8fwL6FNS6Ld0CZbCwuBwmhsmh2+iHRdvVRuGo4ce195BcfTSUXhCU7SyH1+IEIKaWbOxHH88WlMTDcuWgcEANTXtTxIIRFz4v2puZo7dxpctLfQ1Gnm4Tx9GJSSi6y6xlBLLL36B6+13iDv7bNDpMOj1pN52GzEnnhgengLtC7I6E3eVs69QA9wVEWF/A7lDIYeGVato2bABTdOQfj+apqFpGpYhxwercKXE8uthYOgkcyaCv1gr3W4mVm/j8m0/stXr5R8ZvXgtN4/zE5OiK/xCYDz+OLBYgqmsRiPutWsJNDTgfPNNrNOmEl9URGNpKaacnHAXzVBarKWwUHn1in2iPH9FRPBUVPDjddeTOmY0qePGtcvkMRcUkDZ1Cjv/cQdoGtY770CTkvqXX8b1+htIgwHzr38dvIGsXRs1m3/0eplvt7PK2UiCTse0dCuXpaQQG+nirIQEcDohNxc2bw5vjr34Yjwffoihb1/8GzdiOvlkPB9/TMq115I04jTMBQXULVmCbU4J+uRk+i15EkDF8hVdQom/4mfTWfzYXFAQzMcvmYsxK4vE4uJ2+wgINhyDYJMxTcO5fAUAIi8PT5RSNgFsfj8La+285HBgEIJxqamMS00jKVoFWk4n5iFDiLv0Eupum4nIzkZWV2Pu0xv9iSfgfGMVmExYp01Fu/hP4S6bDW+8gW3xkySOHk3yRReGQz6q/46iK6iwj+JnEwrjeCoqwtuEEKSOGxce8ddxH0P/HIiNxTJsGGk334x15m3hY2VVVVTsbggEmG2rYWTVJl5yOLg4OZnVeQO40ZoRPeFvxfPNN3jXrIG4OJIvugiRmIA5Oxv31+uJH1WMSExsNyC99skn2X3f/eB00lRWhl6vD95UW2P5KtyjOBAq20fxkwh5+6b8fNzl5Xi3bCGxuDgc3unYZz8U/9c0jaYPP8JdWUFz2WoADAMHYj77LJoWPh4V25s1jefr61laV4tT0zgvIZHr09Pp100FWqZfn4olN4/G5csRsbHo4+OxzpiOKScnnKnjLCsLV916Nm9mx/QZpE++AUNWFubcXNWHR7FPVLaP4pAS8uSTL72E+udfCC7Otg7qbvt+1ry5eDZvpubBh8i45WZ23XkX0ukMniQmBlpa8G/ciH/jxojb7JWSlx0OFtbasQcCnBEXz+T0dPK7aYIWAEYjgU1VpPztb8Sd9EuM/fuHb6DbJ08he/48dDodicXF4apbU35+uDZCF81mcYojCuX5Kw4aTdNoLC1FSknNrNmkjhlN3PDhYe9TSknLhg14qqrwbd9O3dPPID0esp95mqYPP6T+sYVRG44OwQKtN5yNLLDb2ebzMTQmhmlWKyfGdEOBVusNDyEgJobEq/4KdfXEnDCE5PPPR6fTtWtkp0I4ip+L8vwVhwxnWRk7Z9xMn4cfovf0v2GbPYeE3/42LFLu8nKqr55AwOsNLuoKAZpG3Usv0/zyy3taMUQYKSXvN7mYa7PzvddDgdnME9nZ/CY2igVavXoFaxSMRiwnnRTuPRQ/ejTNq96gcdFi8HppXL4cvcFAYnExjaWl1MyaTd8F85XwKyKGEn9FlwnF8ePOPRfr9u3hsENotmsgEMBZVobf70fzeLD8Zjju/34AbjcAze+/HzXh/6K5mTk2G+vcLfQ3GpnVJ5ORCQnRL9DavTv46PXiXr+e1KlTMGZm4tuxE5c/AF4vccVFJJx5Znhh3DZ7DimjL8WUnx9dWxVHFUr8FV2mbetfx8oXif/Nb8Jhia2TrsM8II+WDz8KNxxzr36z/QlqaiA+HlyuiNm4we2mxG7jw6YmMgwG7uzVmwuTkjBGUfRNZ59NYOtWAjYbicXF+OrraClbja51Tm7z2rXUzZ9P+uQbMGZnt1soNxcUYL1xGjWzZoe/X4UiEijxV+xF20web2VlODfflJ+P9cZpxJ17LkDYMzUOGrSX8O+TCAn/Fq+X+XYbpU4nSTodf7NaGZucgqUbFkS9778PmkbKZWNxvfc+fUpKaB48mNhhwxFCULd8BdapU/YqfoNgimzbxV2FIlIo8T/K6axAq62HXzNrNhk33Uj8yJHULVlC/fMvkCEl9jklmHJyAHCuWdNe+FNSIjo1qy27fD4W1tr5v4YGTEJwbVoaf01JJSHKefrGU07Bt25dcLaAXk/aDdeTNn48yRddhJQSx8oXMWVnk1BURN8F8wH2Gc9XfXcU0UCJ/1FOZ6192w73kFKy+4EHaf7qKxzLloPJhN/vJ3XKZBo/WEPdokUIgwERE4N+4ED869cHhf9AvwB+Jo5AgMW1tSxz1BOQktHJKVyTlka6oXv+l/Z9/TXodIj4eIRejykrC51Oh6WwECll8Eb64EPh3H3VgkHR3UTsX4oQ4hHgd4AX2AT8VUrpEELkAOVAZeuun0gpr42UHYr9EwrltF1cbDflKTcXtACNZauJLxqJa1Uptn/cEdzR6wUpkQYD+P1B4Q8RIeFv0jSeq69jaV0dTZrGBYmJXJeeTpYx+gVa4tjB6L2+YI1C66J2/JjRuF5/g90PPYw5Ly9cpGXOzQ3eEGnfcVOh6C4iGRB9C/iFlPJ44Hvg1jbvbZJSntD6p4S/G/FWVmKbPQdvZWV4W9uOnJbCQrKeeILUP19B74ceInHs2GA/fY9nj8B3zNmPiTn0draKflHVJubZ7ZwSG8srObnc3ycz8sIfKgJLSAg+toaULFnZaC4XSWPHkjBmDLqUFOJ++Ut0RiOpl1+Gpmlsu+Za3OXlWAoL6fv4wvACucrfV3Q3EfP8pZRtUz0+Af4UqWspfjptQzzu8nKMgwZRv3Qp9StW0nfBfIyDBlG3dCmustV4d+3Ca7cf+KQtLYfMvoCUvNbYyAK7jR1+P6fExrIg3cqQCNxg9onbHWyprNeTeNlYdGnp1M+bh+fzz0m76q84Vr5I1ry5pF9ycbA9tZTUPfU0xqyssLev4viKnka0AqRXASvbvM4VQqwDGoG/SynXdHaQEGICMAGgX79+ETfySOBgJzSFRKllwwaqJ04i7qwzaVi2nMSxYwgEAtQ/8ACuVaWg19G4bHkUPkEQKSXvuFzMtdvY5PVyrNnC3b378OvY2Kh7zIbjjiP+jNNxLl9BzIknYuyfg9DpMPTpg04nyJxbgrmgAG9lJZbCQnr9fSa2hx/BlJtL38cXqvCOokfys9o7CCHeBnp38tZMKeWrrfvMBIYCF0kppRDCDMRLKWuFECcBrwDHSikb93ct1d6ha4TG9x1oMbFj47WGVavYfe99aD4fNDUhLBZ0MRY0TWIqLAy2WG47KD2CfNrcxBybjfVuN7kmE5PT0zk3PqHbwiSJY8fQ/N8PSBl9KXXPPgdaABnQkH4fsrmFzH8+gjk3N9zLKIQK7Sh6AhFp7yClPPsAF/0LcD5wlmy9y0gpPYCn9fmXQohNwCBAKfshoCuLiVJKGktL2f3PWaSMvhRd797U3P8AsSPPRXM0EHDU4/v0MwwnDMHzzbf4Q6GeCAv/t+4WSmw21jY309tg4J7evbkgMQlDdwmoXo9hcCGN/36FXnffRfJ55xE3fDhSSrybN1Pz8COkTJhA/MiReCsrw8IfasimhF/Rk4lYYzchRBEwGxghpbS12W4F6qSUASFEHrAGOE5KWbe/8ynP/8BomoazrGyf3R5D3r6maWyfOAnLSSfhWr066NH7fO13NhiCYm8w7P3eIabK42Ge3c6bLicpej0TUtMYnZyMOdoFWq1ZS8RYiDv9DKSm0bw62HY6deoULP36hb/btr+c2qbLhl6rweiKnkJ3NHZbAJiBt1r/EYRSOk8D7hZC+IEAcO2BhF/RNZxlZeyYPoNMIGnUqL3ed5eXUz1xErFnnUXA7cb13nvthb9tbn4og0fTImbvDp+Px2rtvNLQgEXomJSWxpUpqcRHq0BLr4NA6+cTgvR77kbW1GDMzMT+z1n4PR4SRo9G6HQYM7Pafbft0mE7DEs/2IXdg12nUSgOBZHM9jlmH9tfBl6O1HWPZhKKisJ93jVNo2HVKpASc2uffXdVFf7GRhqXLQOzOZiuGSqKCgl/ejrY7ZCcDA5HREI9dX4/i+tqWe5wIIErUlK4OjWN1CgXaFnOPIuA04nvk0/AbCawcyeNK1YQd/bZxJ55Js433yT+5KHY55SQ8seLwlPJOvJzM3k6K7RTKCKN6ud/GNLRU+zY/x2CY/5ss2YDoEtJIen883CUloHTifR4MA4ZEqxKNRhACGIv/AMtr72OdLsjVqDlCgR4ur6Op+vqcUuNPyQlMSktnUyjMSLX2wudbs8vGSEgNjY8RxgI1idoWvCmCFhvupHUcePa9TeKhJeuPH9FJNlX2EeNAToMCXmK7vJyWjZsoGHVKrZdcy3VEyfhqajAU1ER7LpZXIw+NZW44iLqn3se08CByFZv31dVhejVKxje8flofullZHx8RITfo2k8U1fHyM1VPFZby2/i4ng1J5d7e/eJnvAnJCDa1gaYjAidDkwmaLUhZcLV9LrnbnSJiaRNnULKVVe1E37ofF7xz0UVfSm6A+X5H2a09fIBqidOAilJm/43/Nu3Y8zKIqG4GMdTT2F/YRnpl41F16sXu2++JejZut0RjeO3xS8lrzY28Kjdzi6/n2GxsUy1WvmFJYoFWhAUeK8XwzHH4K+uBilJnTgRQ5/e1NxzLynjxmHumx2O5Xe2kBv6RaW8dMXhxr48fyX+hxlt8/jNBQU0f7cB14dr8NXU4Fy2HIQgZtgwvBs3En/2WTSuKiW+uJiWTZvwf/llVPL0pZS86XIy326nyuvlOIuFaelWTo2Li/i194vZjC42loRRxTS9/Q5JY8ZQO3cuxMWR8+wzxAwe3G73A2VPKRSHA2qM42FOSIjiR45sl1ni/3ErdfPmg5SI+HgsQ4bQ8tFHiJgYGv79CrS04FwencpcKSUfNzdTYrfxrdvNAJOJeZlZnBUf371eckwMcSNG4P7qK1IuvyxYqOX1IqWGLiEBuQ9hD/U9MufmqoVYxRGHEv8eSGehhY5pnKEh6YZ+/ej14AP4d+wkdvgwAFxDh+KvqcH50kvBE7Zd6IwQX7e0UGK38WlzM5kGA/f37sPvEhPRR0H0Dccfh3/9N+03tqlGtpx4Iu7168m4eQbG/v3DvXccLyyj1x3/CA9OcZeXt/vOO/Y9UqEexZGEEv8eSCjWnDm3BO+WLZhycsJzc+NHjgSC4Z8fx41HSknKlVfS8OyzaFKjftFizCecEGzHECKCgrXR42Ge3cY7Lhepej23ZmRwaVIypiiGSfQDB+Ivr8B0+gi8b70Nej3W++8jsHMn9YufxP3JJ6RPnQJSUj3+aoTJRO+/z4TWqVk6na7TthihhdiutsxQKA4nVMy/hxHy6L1btoCU7Jw+AxEfT+r4cTheWIZ1+nQkEt/2HdQ9+SQiEAiGLZqaIj5ApS3bfV4etdv5T2MjsTodV6WmckVKKnHdVZULEBtL8kUX0lC2ml633Iw5NxdN02j5+GNihw1nx3XXIX0+Mmbehjk3N9yGIeT1Q+f9eNQir+JwRi34/kwOtQB0dj5N06hbsoS6Z59D6HT0WTCf+qeeouWTT5FSojU3BwW+dXAIRiP4fJjPORvPW2//bJu6gt3v54naWlY66tELwdjkFK5OSyM5mmMTLRZwu9Hl5aFVVYHJCEJHxt13k/K782ksLaXm3vvAaCTl8suof/4FrDOmY87NbVeQtb+sHoXiSEHl+f9MDja/u+1AlM5wl5fz44RraFi1Ck3TcJeX01haiq1kLnHnnE3SmDE0ffQRrtIykv/yZxJGFe8ZopKSAoB+0KCgbZ99fmg+5H5wBgLMs9sYWbWJFY56LkxKojQ3j+kZGVERfuPxx+8ZqtJ68wsKv4mUaydiSEpCb2hjR6vwO15YhnS7sT30MDqdjpjBgxFCtMutV5O1FEcjKubfRQ5WIA6UIw6AFsD20MNIKbHddz+Zi57AOnUKtUufQmtoQMTFkTR2DIY+fah79jmMQ47H9/X68HD0wHffBc/T0HDIP28It6axzFHP4tpaGjSN4oQEbki3kmOK7thE3/r16AYMQNu0CQDzqaeiS0qi5c03MfXNJuWKy6mZNRshBLbZc+h1y80kFBUR/5vfIKUMi3xnqEEriqMRFfaJAB3bLYTCOh1z9N3l5QQCGvVLl+JavZr0KZOJaxWrpg8/onbxYoROh/T7wOuDQACRl4esqvr/9u48POrqXvz4+zMzmSX7CglJMICIBHd4lFXRq7JYRexPRfFqFdeiCPbX5Q50egAAG0xJREFUKioWa12uKKuKl2qLFRC1tpW2srhSFdtCq+ICqGwSwpKF7MkkM3PuH7MQkgwgMJOQ+byeJ08m5zvLyXkyn5w52yfiv0OTMfypspL5ZaXs8XgYmpDA5MwsCoO97/ZiteK66CLqV6zAmp5O8iWjSRpzObsmTybrnikkjRzZaleuUrFMh32iyL1xIzsn3R0aXghynHwyuXPnhIaCnH37UrPsTWqWL8d52mmUL1pM0e130LR9O9acHPB6MbW1UN8QWrYY6cDvM4blVVVctnUL0/fsJscWx0v53VmQl98ugd964v7zAS15eeD1YktPQ5KSSLn+eqrfex+LRcibNze0ckePSlDq0HTYJwJaDhG1nNwtumsSmVMmA1D9/gfEDx1Kw5YtpF03HmMMex99FG91TeioZcnuitm9J6J1NsbwUV0ts0tK2OB2c5LdwbO5uQxPaN8NWnFnnol3xw5wu/FVVpJ63XiSxlxO/bvvkTRsKEnDhuLeupXkUaM04Cv1A2jP/yi1NbHb8qCu5pPFcSedRPx557J7+sPsefwJUq+6kvoNG0i+4Hz2LVpM+Qsv4q2rPyCBitm9x79RK0I+ra/jhh3fc1tREdU+H/+Tk8MbBQWcH63UiRaLP0F6Xm6oyHrKKf5LNdVYkpNJGX8tVpcLW5cuuAr7kv/8fFyFhTRt386uX9xL9YoVka+nUp2I9vyP0uEsE2w+3FO9YgWVi5eACOk3T8CWk4M0NVH11nLizzmHmr+vhsbG1k8SgR2637gbmFNSyvu1NWRYrUzr0pUfp6Zij3YP2ucDnw9TtBOsVhIuugjXeedR+vB0Gtb9m+yp95E0ciTl2dmUzJqNPS8vlKymeQ4DpdThi2Qax+nALUAwheP9xpi3AtemAhPwZ/KaZIxZeajn66gTvm2t12+5qqdhwwbcW7ZQMuMpMv//z/B6fdSvXk3yDddTPOFm//LNQACMxkatHY2NzCst5W/VVSRaLExIz+C6tDTio7FBq60k8MGlol4vxMVhcblABCNC9rQHQ6dttjxoTTdfKXVo7XWw2yxjzFMtKlIIjAP6Ad2Ad0TkJGNM5I+bjIC2lgk2/zQAUHT7HXiqqxBbHLunP4zj9NNxr1lDzfr1ByYTgYgG/hKPh+fLSnm9ogKbCBPS07kpPcobtLxeSE/H0qsnvrXrEIcD43bjGDyIxi+/wni9/qWZVivZDz4QCvzgb+vgRi3QDFhKHY32GPYZAyw1xriBrSLyHXA28MnBH9YxHE5vs/kwj71PHxz9z8KzfAX2wWfiXvMJ7jVr/HfcuTMqda7yenmxvJxF+8ppMob/l5rKHRmZZEU5bWLwXH3Ky/GVl4PTQddHHqHh88+pWPIKqddeQ9KYMXi2b8feo0doQ1ZQy2Cvm7OUOnKR/px/p4isF5HfikhaoCwX2NHsPkWBsuNCyyxa9V9/HZrsDU7++nw+aj76iO9vu52SBQuoff8DcDiI696dlDvv9OfPjcIQS73Px2/Kyrh4y2ZeKC/jgsQk/tqjJw91zY564Je+fUNzGfaBAxGXi66/eoTUH/2IrvffT+bku6leuQrP999TOnsOTdu3t3qOlsFeM2ApdeSOKgKIyDtAdhuXHgDmA48AJvD9aeAmoK13aptjHSJyK3ArQPfu3Y+mqkek5dh9cLVO1j1TMMaw846fggj5z88Pnf644/Y7sJ/Um/o1n4DVyr5nng2NcdcsfRUSE0M5YiOlyRjeqKhgflkZJV4PwxMSuDsziz7tuEHLbN0KLhfx551L3Sf/gIYGvLt3AWCxWEgcOpSKxUuwFxSQdc8U9j49s9U5+roTV6lj56iCvzHmwsO5n4j8Bvhr4MciIL/Z5TygOMzzLwAWgH/C98hremRajd3fNYmse6ZQMnMW3ebMDqROLMbr9eL1eqn56CO8NTXUf/QxtpN64/nmW/9Qh9cLffvChg1QUxOx+vqM4a3qKuaVlrKjqYn+LhezMrtxVnx8xF4zHMnLw+zahTidGLcbiYsj5bJLqX7/A7KnPYinuJjK114ncejQ0NEL+c/Px3HyyaEAb+/TJ+r1VipWRHK1T44xZlfg9hTgHGPMOBHpByzBP87fDXgX6H2oCd9or/Zpvn4/GJwaN23C3qcP7o0bqd+8mb2/nI6pq0Pi40keezmVS18lbeJE6teto+E//wFjcA4fTsMHH0S0t2+MYXVtLXNKS9jkdtPH4WBKZhbDEhLab0hEBBIScJ5+Gg0frwG7HVtaGl3u/QXJo0YB/n+uxpjQ0crBoK/n5yt17ET9SGcReRk4A/+Qzjbgtmb/DB7APwTkASYbY5Yf6vmiHfyDASjlqivZ9/Ii0v77OjJuvjmU+GPrjTdBZSX2006jcf16sNtx9u+Pc0B/KubOi1o91wXSJv6nvp78uDjuzsxiZFISlvYI+nFxkJ8PgSMobKeeimfbNizGkDntQaw2W+gIhqBDLZXV8Xyljo6e5/8DGWOoWr6c4kd+DdXV4PWS89QMUi+5xH8K54IFlM2eA04ntvx8PN98E7W6AWxoaGB2aQkf1tbSxWbjjowMrkhJJS7KwdI5eDAJZ5/Nvt//Hl9NTesNag4HOY/+GkfPnhTdNYm0cVeTPmGCJkRXKko0gfsRsBcUkDp6FBWLl4DTib2gIDQcFD94CO7du6l57XV/4I9CnlyA7Y2NzCst4a3qapItFn6WlcX41DSc7RBMLb16kXTZpbhOPBFbt27sefhh4s8/n8QLLsBgqHv/A2rXrMEemLhNvfoqSmbNxtatG6mXXBL1+iql9tPgH4Z740a+n3gniecPJ2PyZOKHDEZEqP/6a7ZPuBnq6/eP44tEPPjvaWpiflkZb1RWYBfhtvQMbkxPJzmaG7QCbAMH4vnHPzDFxZTcNxVJTiZ+yBCorcW9di2pIy7GXlDAvpmzsNpsodNNE4YMYd/Cl7AXFES9zkqpA+ln7zDsffqQcNaZVC15hbq1a2nats1/ONuWLf4VO80ncI3Zn0f2GKvwenlq715Gbt3CHysruDo1lRU9e3F3Vlb0An9m5v7bVivOEwLLbuPiSBl/LRhD7fLl4HKResP17H16JgC5zz3rz5cbWCrrKiyk+4sv4CosjE69lVJhac8/DPfGjdStXYu1e3fqP/6Yhs8/I/2WW6j//nt/oI9wT7/W5+PlfeX8rrycGp+Py5KTmZiRSV6UM2gBUFq6/7bxUbNiJRlTJpM4dGgosFcuXkLG7beRecstJJ97big3bums2Th79gxtxtLVO0p1DBr8wzDG4Kutw5SUBtbq+/wTvMEJ8ggF/kZjeL2igufLSinzerkgMZG7M7Po7XBE5PV+MLGAx0Pi0KHE9+tHw4YN1H6wmqyf3ROayA0GeD1+QamOS4N/GD6fD+N2gzFYTz8d79rIJkn3GsNfq6p4pqyUnU1NnO2K55msLE53uSL6uocj4eqrcOTkUL9xI/UrV4EIjdu24Sos9G/OemZeaFlmy2Wa2tNXqmPS4B9G3Zo1oXH8SAZ+Ywzv1dQwp7SE7xobKXQ4mJ6Xz+D4+KiucZfzzsXZNRtf0Q68Thee997zX3A4kKpqyl97HUlMRBITSZtwEyVPzjjg+IVg0G9r05ZSquPR4B+G45xzIv4a/6qrZVZJCZ83NFAQZ2dmt25cnNg+G7TMmk+o93j8w1p2OyQkkDxmDHFdulA2bx6p468l+fLLQ+foVy15Bdh/BEbw2IvcuXN0qEep44AG/zCq/vzniD33Vw0NzC4p4eO6WrJtNh7pms2YlBRs7bSb1ZKdja+uFtfgwYiBupUrSb3qSrrefz8Aju75ByRQadiwgbzAYXYAefPm+o+u7tFDd+UqdZzQpZ5hePdVHPPn3NroZkrxTq7cvo2v3A38IqsLy3v05MepqZEL/Ic6utlqxVdbC00e6leuwtn3ZFKvG0/1e++H8uK2TKCyc9LdobX7wXH94ESvBn6ljg/a8w/D3vtE6lceMrvkYdnV1MRzZaX8ubISh1i4IyODG9PSSYzGOv229h/YbKG0kSlXX038WWdSMmMGiT++gsrFSxCLhfTrxlMycxYAJTNnhYZyjDHkzp2jwzpKHec0+EfQPo+HBeVlvFJRgQHGp6Vxa3oG6dFMpGKzHfgPIPgPx2pFXC5SrhhLfL9+OHv1wt6nDyljx9K4bRtJI0eSOHToAcM5wV5/3ry52sNX6jinwT+MyuBqlyNQ6/PyUvk+frevnHqfjzHJKfw0M5PcuLhjWMPDkJRE2k030rRnLzV//OP+5PCNjeCwk/2rh4nv1y80dGOMoWn7dkpmzgptzAJ03b5SnZAG/3C2f/+DH+L2+Xi1ooIF5WWUe71clJjIpMwserXDBq2kiT/FUlFBxeIlpP/3ddR4POBwkHDppdQuW0bW9Ok4evQ4sP4bN1IycxZdfnZPmwFe1+0r1XnohG84uYefVthjDH+srGD01i08UbKXPg4Hr3Y/gTm5ee0S+HE6seyroPKVpUhTI65Bg0kcOQKL00n9229jdbmwxdnYOelu3Bs3hh4W7NknjxqlwzpKdXLa8w9n8+ZD3sUYwzuBDVpbGhs51enk19k5DEpIiEIFAacDGtwggrWwEO9XX/nLfT4q//IXMiZNIi4vF8/2bdQsX0HG5LsPSJvo7NnzgB6+9uyVih0a/MM5xNk9awJpE79oaKCn3c7cbrn8V2Ji5HvMDsf+E0UbAt+7dcPs2LH/Po2NYLXi2buHiiVLSL32WqxpaaHzeII00CsVuyIW/EXkVSCYgTsVqDDGnCEiBcAGYFPg2j+MMbdHqh5HLCsL9u5tVby+vp5ZpSX8s66OHJuNR7OzuSw5BWu0hkla5gIWYOdObKefhu3kviT2PwtvcTH7freQyleWgstF+cKFZNx0owZ7pVRIxIK/Mebq4G0ReRqobHZ5szHmjEi99jFRUnLAj9+53cwtLeGdmhrSrVamdunC1Smp2KOQQUvy8jBFRYEfZP/JogAWK7ZePWn8fD1Nm7fQddzVOC69lMRhw3Bv3YoxhtLHn2Dfy4tIGjYs7D8AzZurVGyJ+LCP+CPJVcAFkX6tYyoQYHc2NfFsaSnLqipxWSzclZHJ9elpJFiil0HL7N0bWq/vHDSIhjVrQtnDJCEBU1oGTidpEyZg79MHESG+Xz/i+/XDGIOrVy+Agy7RDJ7RoweyKRUbojHmPwzYY4z5tllZDxH5FKgCHjTGfNjWA0XkVuBWgO7du0e8os2VeTz8b1kZr1ZWIMD1aenckp5OWjQ3aAX5vJCQgOuUU0m9cyLVKSnUBDJnpd10I/bcXBChbPYcUoafd0DwFpHDypyla/iVii1HFclE5B0gu41LDxhj3gzcvgZ4pdm1XUB3Y0yZiPQH/iwi/YwxVS2fxBizAFgAMGDAANPyeiRUVlby9NNP89SWzTQaw9iUFO7IyCQn2hu0mpGkZJyFhf6MYuvX0/WXD1H/z3+SesP1VCxajFgs5M1/7qiCt670USq2HFXwN8ZceLDrImIDrgD6N3uMG3AHbv9bRDYDJwHrjqYuR6u+vp7nnnuOxx9/nLKyMkYmJzMpI5OCKKVNlJ49sffqhfvDD/3HMQSOZEi8ZhxJ/fuz98kZSEICXR6aRuKIETR89pl//X5g6aYeqqaU+iEiPYZxIbDRGFMULBCRLKDcGOMVkZ5Ab2BLhOsRlsfjYeHChUyfPp2dO3cyYsQIHnvsMVzXjo9qPcyWLTSVlWGNjyfzvnvxFBfjAyoWLcZ1+uk4+59FwvDzSb3kEspffJGKRYup+stfOGHhQu2xK6V+sEgH/3EcOOQDcC7wKxHxAF7gdmNMeYTr0YrP5+MPf/gD06ZN45tvvmHgwIEsWrSI4cOHA7AhIQFqa6NXIYuFzAfux7d7N46evSifO4/MKZOxGB8lD/8K6uup+/AjEk7qTdpNN+Hz+Yjr1g17nz6Hfm6llGohosHfGPOTNsreAN6I5OsejDGGVatWMXXqVD799FNOOeUU3nzzTS699NIDh03s9qgGf+eF/4Vv1y5K58wlJzeX3Llz8Pl8dJk6FUt+d2qWvUny5WNDSzGTzz2Xorsm4ejZk6bt20PJVpRS6nDE1A7fTz75hKlTp7J69WoKCgp4+eWXueaaa7C2da6+13vsK2Cx7N85bLH4v0QQmw33v9bS9NnnpF57DXEnnABA8U8nggj5z8+n27RpBzxVcHVOw5Yt7PrFvXQDUkaPPvZ1Vkp1SjER/L/88kseeOABli1bRteuXXnmmWe45ZZbsB9sMtcc/HiHI9I8+AeSqTgGDyJt7BXYexT4j1N+cgY1q94m7/n55D0/HzhwfX7zzVjOvn1D6/qTRo489vVVSnVanXqcwOv1csMNN3DaaaexevVqHn30UTZv3szEiRMPHvgBUlKP7sVbDsFYLFgCPXq6dkFSU3ENGULjpm8ofeoprFYrKaNH+5dsBvLjtrWCJ7gZK3gap8ViIWX0aB3yUUr9IJ265x8czvn5z3/OvffeS3p6+uE/uO4wx/sTE6GmpnV5sIfvcuE6dxj1K1fh27wZ+0kn0fjtt5AQT9aUyaGgHRzLb74hq2HDhla7bnUzllLqWOjUwR9g4cKFR7b+vbqNgN6WtgJ/QNK4cST0P4uGHUXUr1yFY9Agcp+fT/mMGVSvXIXFYjno7tu2Ar1uxlJKHQudPvgf8can5CQoO7IVqLbeJ+IrK8eek03pjKfw1NWBCAnnnINv82Zq3v+ArlPvO2QQ10CvlIoUHSgOwzbs3EPfyWr1f1ks4HBgOdF/gJqvtIyUS0ZT/spSUq8bT9eHpmFJSqJi0SIA8p+ZF8qWZYyhYcMGjInK6RVKKQVo8A/L8/nnB7+DCLbzzoO4OGyFhUh8PLmPP07GlMmk3nA9VStWkn7NOCoWL8FisZD/u9/SfcH/tprIbTmBq5RS0aDBPwxbly5tXwimaDQGz3vvQUMDni+/xFRVUfPxx5QvXEj5Cy9iGhuJ69YN09jInsceC03mthyG0glcpVR70OAfhjU/v3VhQvyBu37tduJHjvAP/fh8VLzwIsbdiFgsGBGaiovJun8qYg0/taKHsiml2oMG/zAcLfcBZGZiGzhw//VBg8Dp9AftlBQSRo7AWK0QZyP50h+R9pOfUDZ3HmKxhIZ7lFKqo9DgH0bSmMsO3KhVWorn3fdwDhkCiYm4v/gCGhupXb4CV2FfGj5fT85D08icMIHKpa9iz+1G1uS7SRo5Unv2SqkOp9Mv9TxSvh07/Kkc4+IQux37qafi/uILsiZPxvP9dnY//gRJF19M1Ztv0rBxEzn3TyV51CiMMdjz8og74QSK58wlcehQrNrrV0p1MBr8w4grKACXC6xWMm6ewL7FS7Da7VitFhJGj0ZEKJk5i8zbbiVhyJDQZK6IkDJ6NMYYnchVSnVYGvzDsFgsWOx2ECF+yBDi8vKwFxSEhnCSR40CYO/TM0kaNqzVsI5u0FJKdWQa/MNw9u1L/osv0LhtGyJC6azZ5M2bGwrywX8Ajh49tHevlDruHNWEr4hcKSJfiYhPRAa0uDZVRL4TkU0iMqJZeX8R+SJwba500JnQ4BBOyf88CUDevLnY+/Q5YDeuLtNUSh2vjna1z5f4E7T/vXmhiBTiT+HYDxgJPCciwYwp84Fb8efu7R243nEF/gk4+/alcdMm3Y2rlOoUjir4G2M2GGM2tXFpDLDUGOM2xmwFvgPOFpEcINkY84nxd59/D1x+NHWIJGffvuQHztYH3Y2rlOo8IrXOPxfY0eznokBZbuB2y/I2icitIrJORNaVlJREpKIH03JYR4d5lFKdxSEnfEXkHSC7jUsPGGPeDPewNsrMQcrbZIxZACwAGDBggB57qZRSx8ghg78x5sIjeN4ioPnhOHlAcaA8r41ypZRSURSpYZ9lwDgRcYhID/wTu/8yxuwCqkVkYGCVz/VAuE8PSimlIuRol3qOFZEiYBDwNxFZCWCM+Qp4DfgaWAFMNMZ4Aw+7A3gB/yTwZmD50dQh0jTZilKqM5LjJagNGDDArFu3Luqv21YSdaWUOl6IyL+NMQNaluupnoegyzuVUp2RHu9wCHpGj1KqM9Kev1JKxSAN/kopFYM0+CulVAzS4K+UUjFIg79SSsUgDf5KKRWDNPgrpVQMOm52+IpICbA9yi+bCZRG+TU7Om2T1rRNWtM2aa292uQEY0xWy8LjJvi3BxFZ19a26FimbdKatklr2iatdbQ20WEfpZSKQRr8lVIqBmnwP7gF7V2BDkjbpDVtk9a0TVrrUG2iY/5KKRWDtOevlFIxSIM/ICJXishXIuITkQEtrk0Vke9EZJOIjGhW3l9EvghcmxtIS9kpich0EdkpIp8FvkY3u9Zm+8QKERkZ+N2/E5H72rs+7UVEtgXeD5+JyLpAWbqIvC0i3wa+p7V3PSNJRH4rIntF5MtmZWHboL3fOxr8/b4ErgD+3rxQRAqBcUA/YCTwnIhYA5fnA7fiz0/cO3C9M5tljDkj8PUWHLJ9Or3A7/osMAooBK4JtEmsOj/w9xHsQN0HvGuM6Q28G/i5M1tI6zjQZht0hPeOBn/AGLPBGLOpjUtjgKXGGLcxZiv+vMNni0gOkGyM+cT4J01+D1wexSp3FG22TzvXKZrOBr4zxmwxxjQCS/G3ifIbA7wUuP0Snfw9Yoz5O1DeojhcG7T7e0eD/8HlAjua/VwUKMsN3G5Z3pndKSLrAx9tgx9dw7VPrIj13785A6wSkX+LyK2Bsq7GmF0Age9d2q127SdcG7T7307MpHEUkXeA7DYuPWCMeTPcw9ooMwcpP24drH3wD3E9gv93fAR4GriJTtgOP1Cs//7NDTHGFItIF+BtEdnY3hXq4Nr9bydmgr8x5sIjeFgRkN/s5zygOFCe10b5cetw20dEfgP8NfBjuPaJFbH++4cYY4oD3/eKyJ/wD2HsEZEcY8yuwFDp3natZPsI1wbt/rejwz4HtwwYJyIOEemBf2L3X4GPb9UiMjCwyud6INynh+Ne4I82aCz+CXII0z7Rrl87Wgv0FpEeImLHP4G3rJ3rFHUikiAiScHbwMX4/0aWATcE7nYDnfg9chDh2qDd3zsx0/M/GBEZC8wDsoC/ichnxpgRxpivROQ14GvAA0w0xngDD7sD/+y+C1ge+OqsnhSRM/B/LN0G3AZwiPbp9IwxHhG5E1gJWIHfGmO+audqtYeuwJ8Cq51twBJjzAoRWQu8JiITgO+BK9uxjhEnIq8Aw4FMESkCfgk8QRtt0BHeO7rDVymlYpAO+yilVAzS4K+UUjFIg79SSsUgDf5KKRWDNPgrpVQM0uCvlFIxSIO/UkrFIA3+SikVg/4Pro6TIHNNNYMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "xm = 0.05*(x.max()-x.min())\n", "xfit = np.array([x.min()-xm, x.max()+xm])\n", "yfit = slope*xfit + yint\n", "plt.plot(xfit, yfit, '-k')\n", "plt.plot(x, y, \".C3\", ms=1, zorder=-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare execution times of linfit and polyfit to fit a randomly generated data set of 10000 $(x,y)$ data points. On my computers, linfit is about 6 times faster than polyfit for unweighted data and about 3 times faster for weighted data." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "TIME COMPARISON WITH NO WEIGHTING OF DATA\n", " linfit time = 0.04812828300000005\n", " polyfit time = 0.3567316079999996\n", " ratio = 7.412099201627434\n", "TIME COMPARISON WITH WEIGHTING OF DATA\n", " linfit time = 0.08441347999999937\n", " polyfit time = 0.4123130230000003\n", " ratio = 4.884445268694091\n" ] } ], "source": [ "import timeit\n", "setup = '''\n", "from linfit import linfit\n", "import numpy as np\n", "\n", "def randomData(xmax, npts):\n", " x = np.random.uniform(-xmax, xmax, npts)\n", " scale = np.sqrt(xmax)\n", " a, b = scale * (np.random.rand(2)-0.5)\n", " y = a*x + b + a * scale * np.random.randn(npts)\n", " dy = a * scale * (1.0 + np.random.rand(npts))\n", " return x, y, dy\n", "\n", "npts = 100000\n", "x, y, dy = randomData(100., npts)\n", "'''\n", "nreps = 7\n", "nruns = 100\n", "\n", "linfitNOwt = min(timeit.Timer('fit, cvm = linfit(x, y)', setup=setup).repeat(nreps, nruns))\n", "polyfitNOwt = min(timeit.Timer('slope, yint = np.polyfit(x, y, 1)', setup=setup).repeat(nreps, nruns))\n", "print(\"TIME COMPARISON WITH NO WEIGHTING OF DATA\")\n", "print(\" linfit time = {}\\n polyfit time = {}\\n ratio = {}\"\n", " .format(linfitNOwt, polyfitNOwt, polyfitNOwt/linfitNOwt))\n", "linfitWT = min(timeit.Timer('slope, yint = linfit(x, y, sigmay=dy)', setup=setup).repeat(nreps, nruns))\n", "polyfitWT = min(timeit.Timer('slope, yint = np.polyfit(x, y, 1, w=dy)', setup=setup).repeat(nreps, nruns))\n", "print(\"TIME COMPARISON WITH WEIGHTING OF DATA\")\n", "print(\" linfit time = {}\\n polyfit time = {}\\n ratio = {}\"\n", " .format(linfitWT, polyfitWT, polyfitWT/linfitWT))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using linear fitting routine for non-linear fitting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Linear fitting with weighting can be used to fit functions that are nonlinear in the fitting parameters, provided the fitting function can be transformed into one that is linear in the fitting paramters. This can be done for exponential functions and power-law functions. This approach is illusutrated in the next two examples." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Using an exponential fitting function with linfit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nuclear decay provides a convenient example of an exponential fitting function: $N(t) = N_0 e^{-t/\\tau}$.\n", "\n", "Here are the $N$ vs $t$ data together with the uncertainties $\\Delta N$." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "t = np.array([0., 32.8, 65.6, 98.4, 131.2, 164., 196.8, 229.6, 262.4, \n", " 295.2, 328., 360.8, 393.6, 426.4, 459.2, 492.])\n", "N = np.array([5.08, 3.29, 2.23, 1.48, 1.11, 0.644, 0.476, 0.273, 0.188, \n", " 0.141, 0.0942, 0.0768, 0.0322, 0.0322, 0.0198, 0.0198])\n", "dN = np.array([0.11, 0.09, 0.07, 0.06, 0.05, 0.04, 0.03, 0.03,\n", " 0.02, 0.02, 0.015, 0.014, 0.009, 0.009, 0.007, 0.007])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Linear and semi-log plots of the data with error bars:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAADrCAYAAAAhQUHcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAXx0lEQVR4nO3df2zcd33H8derToy5pgstNYg1jRJIU61iG4fOhrUT8hijBRqKENJIYWJJVYMGoxVIrBUIDU2obAjUIZBWA0mLoKnQgK2pyko1sCoGw3a4wtKF0tJ2JStaElpC3RM4Td/7487gumen9vc+/v6450M6+e7ru+/3/VGST973+Xw/n7cjQgAAAFhbp+UdAAAAQD8iCQMAAMgBSRgAAEAOSMIAAAByQBIGAACQA5IwAACAHKzLO4CFzj777NiyZUveYQBYIwcOHDgWEcN5x5GF7R2SdpxxxhlXbt++Pe9wAKyhrH2Yi7RPWKPRiJmZmbzDALBGbB+IiEbecfQC/RfQf7L2YUxHAgAA5IAkDAAAIAckYQCQge0dtieOHz+edygASoYkDAAyiIj9ETG+cePGvEMBUDIkYQAAADkgCQOABFrNpo7dMKFWs5l3KAAKqlD7hC02NjYmSZqcnMw1DgBYiVazqft2Xq4BSeuGhrR57x7V6vW8wwJQMIyEAUAG3W7Mb01Na/C00zQgKU6cUGtqOr8AARRWoZOwrXNzunj2CYbzARRWtxvza6Mj8uCgNDAgr1+v2uhIjhECKKqk05G2H5L0uKSTkp5cya6yrWZTVz/6mAYkPbxrN8P5AEqjVq9r8949ak1NqzY6Qt8FoKu1uCfsTyLi2Eo/1Jqa1oD0tOF8OjIAZVGr1+mzACyrsNORtdERnVR7CI3hfAAAUDWpk7CQ9A3bB2yPr+SDtXpd1591pvZv2MBUJAAAqJzUSdhFEfFySa+T9G7br1r8Btvjtmdszxw9evRpv3twcFB3bDidBAwAAFRO0nvCIuKRzs8jtr8maVTSXYveMyFpQpIajUYs/B37gwEoOts7JO3Ytm1b3qEAKJlkI2G2T7d9xvxzSa+VdDDV9QAgD9SOBLBaKUfCXijpa7bnr3NzRPxbwusBAACURrIkLCIekPSHqc4PAFXWajbZZwyouELXjgSAftRqNvXwrt2KuTl5cJAV4kBFFXafMADoV62paT35q19JTz1F7UmgwkjCAKBgaqMjWjc0RO1JoOKYjgSAgqH2JNAfSMIAIINU+4RRexKoPqYjASAD9gkDsFokYQAAADkgCQMAAMgBSRgAAEAOSMIAAAByQBIGAACQA5IwAKioVrOpYzdMqNVs5h0KgC7YJwwAKqjVbOq+nZdrQNK6oSHqTwIFxEgYAFRQa2pag6edpgGJ+pNAQZGEAUAF1UZH5MFB6k8CBcZ0JABUEPUngeIjCQOADFLVjuwF6k8CxcZ0JABkQO1IAKtFEgYAAJADkjAAAIAckIQBAADkgCQMAAAgByRhAICuKHsEpMUWFQCAZ6DsEZAeI2EAgGdoTU1rQKLsEZBQ8iTM9oDtpu3bUl8LANAbtdERnZR0UqLsEZDIWkxHXiXpkKTfWYNrAQB6oFav67x9N1P2CEgoaRJme5OkN0j6qKT3pbwWAKC3KHsEpJV6OvJ6SR+Q9NRSb7A9bnvG9szRo0cThwMAAFAMyZIw25dKOhIRB5Z7X0RMREQjIhrDw8OpwgEAACiUlCNhF0l6o+2HJN0i6dW2v5jwegAAAKWRLAmLiGsjYlNEbJH0VknfjIi3p7oeAKBYdl14oa5rjLDZK7AE9gkDgEVsn277Jtuftf22vOMpo1azqasffUw7Zmf18K7dJGJAF2uShEXEZERcuhbXAoBubO+xfcT2wUXHL7F9r+37bV/TOfxmSf8cEVdKeuOaB1sBbPYKnBojYQD6xY2SLll4wPaApM9Iep2kCyTttH2BpE2Sftp528k1jLEyaqMjWjc0JA0MsNkrsARqRwLoCxFxl+0tiw6PSro/Ih6QJNu3SLpM0mG1E7G7xZfVVanV69q8dw+bvQLLIAkD0M/O0W9HvKR28vUKSZ+S9Gnbb5C0f6kP2x6XNC5JmzdvThhmObHZK7A8kjAA/cxdjkVEPCFp16k+HBETkiYkqdFoRI9jA1BxDLMD6GeHJZ274PUmSY/kFAuAPkMSBqCfTUs6z/ZW24Nq72l460pOYHuH7Ynjx48nCRBAdZGEAegLtvdJ+q6k820ftn1FRDwp6T2S7pB0SNKXI+KelZw3IvZHxPjGjRt7HzSASuOeMAB9ISJ2LnH8dkm3r3E4AMBIGABkwXQkgNUiCQOADJiOBLBaJGEAAAA5IAkDAADIAUkYAGTAPWHptJpNHbthQq1mM+9QgCRYHQkAGUTEfkn7G43GlXnHUiWtZlP37bxcA5LWDQ1p8949lEBC5TASBgAonNbUtAYkDUiKEyfUmprOOySg5yqdhI2NjWlsbCzvMAAAK1QbHdFJSSclef161UZH8g4J6LlKT0dunZvT9rkTajWbDGMDSML2Dkk7tm3blncolVKr13XevpvVmppWbXSEPhyVVNmRsFazqasffUw7Zmf18K7d3NgJIAn2CUunVq/r7HeOk4ChsqqbhHE/AQAAKLDKJmHcTwAAAIqsuklYva7rzzpT+zdsYGkzAAAonErfmL/3O9/JOwQAQI52XXihts+d0FWfneDLOAqnsiNhALAW2DG/uFighaIjCQOADFgdWVws0ELRJUvCbA/ZnrL9A9v32P5IqmsBALAYC7RQdCnvCfu1pFdHxKzt9ZK+bfvrEfGfCa8JAIAkNnxF8SVLwiIiJM12Xq7vPCLV9QAAWKxWr5N8obCS3hNme8D23ZKOSLozIr7X5T3jtmdszxw9ejRlOAAAAIWRNAmLiJMR8TJJmySN2n5pl/dMREQjIhrDw8MpwwGAnmN1JIDVWpPVkRHxC0mTki5Zi+sBwFphdSSA1Uq5OnLY9vM6z58r6TWSfpTqegAAAGWScnXkiyTdZHtA7WTvyxFxW8LrAQAAlEbK1ZE/lMSSFAAAgC7YMR8AACAHJGEAAAA5IAkDAADIAUkYAGTAPmHVNjY2prGxsbzDQEWRhAFABuwTVm1b5+Z08ewTajWbeYeCCiIJAwCgi1azqasffUw7Zmf18K7dJGLoOZIwAAC6aE1Na0DSgKQ4cUKtqem8Q0LFLLtPmO0PL/PriIi/63E8APoc/Q6KojY6opOd5+vWr1dtdCTXeFA9pxoJe6LLIyRdIelv0oYGoE/R76AQavW6rj/rTO3fsEGb9+5Rrc7+4+itZUfCIuIT889tnyHpKkm7Jd0i6RNLfQ4AVot+B0Wy9zvfyTsEVNgpyxbZPkvS+yS9TdJNkl4eEY+lDgxA/6LfAdAPTnVP2MclvVnShKTfj4jZNYkKQN+i3wHQL051T9j7Jf2upA9JesT2LzuPx23/Mn14APoQ/Q6AvnCqe8LYwgLAmqLfAdAv6OwAAEiEskdYDkkYAGRA7Ugsh7JHWA5JGABkQO1ILIWyRzgVkjAAABKg7BFOhSRsGbsuvFDXNUb49gIAWLH5skcnJZmyR+iCJGwJDCMDALKg7BFO5ZQ75verbsPI/AMCAKwEZY+wHEbClsAwMgCgCNjmoroYCVtCrV7Xeftubo+AjY4wCgYAAHqKJGwZtXqd5AsAkKutc3PaPndCrWaT/5MqJtl0pO1zbX/L9iHb99i+KtW1AACoIhaJVVvKe8KelPT+iPg9Sa+U9G7bFyS8HgAAlcJeY9WWLAmLiJ9FxPc7zx+XdEjSOamuBwBA1bBIrNrWZHWk7S2S6pK+txbXAwCgCthrrNqSJ2G2N0j6iqSrI+KXXX4/bnvG9szRo0dThwMAQKk8ODioOzacTgJWQUlXR9per3YC9qWI+Gq390TEhKQJSWo0GpEyHgAAymZycjLvEJBIytWRlvR5SYci4pOprgMAAFBGKacjL5L0F5JebfvuzuP1Ca8HAABQGsmmIyPi25Kc6vwAkIrtF0v6oKSNEfGWvOMBspgvecS0ZvFQOxJApdjeY/uI7YOLjl9i+17b99u+ZrlzRMQDEXFF2kgB9DvKFgGomhslfVrSF+YP2B6Q9BlJfybpsKRp27eqvQfmdYs+vzsijqxNqAD6GUkYgEqJiLs6exMuNCrp/oh4QJJs3yLpsoi4TtKlq72W7XFJ45K0efPm1Z4GSIrak8XFdCSAfnCOpJ8ueH1Yy1TwsP182/8kqW772qXeFxETEdGIiMbw8HDvogV6hNqTxcZIGIB+0G2R0JL7EkbEzyW9K104wNroVnuS0bDiYCQMQD84LOncBa83SXqkFye2vcP2xPHjx3txOqCnelF7cmxs7DcrLNFbJGEA+sG0pPNsb7U9KOmtkm7txYkjYn9EjG/cuLEXpwN6itqTxUYSBqBSbO+T9F1J59s+bPuKiHhS0nsk3SHpkKQvR8Q9ecYJrBVqTxYX94QBqJSI2LnE8dsl3d7r69neIWnHtm3ben1qoBBYXZkOI2EAkAHTkSi6ycnJVe+Wz+rKtEjCEtp14YW6rjHCX1oAQCl1W12J3iEJS4RvDwCAsuvF6kosjSQsEb49AP2BLSpQZayuTIskLBG+PQD9gXvCUHWsrkyH1ZGJ1Op1nbfv5vbuxKMj/OUFAABPQxKWUK1eJ/kCAJTaaldW4tSYjgQAAMgBSRgAZMCN+QBWiyQMADLgxnwAq0USBgAAkAOSMAAAgByQhAEAAOSAJKzgWs2mjt0wQdkjoKC4MR/AarFPWIG1mk3dt/NyDUhaNzREyQiggCJiv6T9jUbjyrxjAVAujIQVGPUnAQCormRJmO09to/YPpjqGlVH/UkAAKor5XTkjZI+LekLCa9RadSfBACgupIlYRFxl+0tqc7fL6g/CQAou7GxMUmrr0OZ9fNFlfs9YbbHbc/Ynjl69Gje4QDAirA6EsBq5Z6ERcRERDQiojE8PJx3OACwIpQtArBauSdhAAAA/YgkDAAAIAcpt6jYJ+m7ks63fdj2FamuBQAAUDYpV0fuTHVuAABQHlvn5rR97oRazeaqVvxn/XxRMR0JAACSaTWbuvrRx7RjdlYP79q94lrIWT9fZCRhAAAgmawl+Kpcwo8kDAAAJJO1BF+VS/iRhAFABmzWCiyvVq/r+rPO1P4NG7R5754V39OV9fNFRhJWca1mU8dumKjUHDpQJGzWCpzag4ODumPD6atOoLJ+vqhSFvBGzlrNpu7bebkGJK0bGqrcNwgAAMqMkbAKq/LNjAAAPFtjY2O/KQKe5zkWIwmrsNroiNYNDUkDA5W7mREAgLJjOrLCavW6Nu/do9bUtGqjI0xFAgBQICRhFVer10m+AAAoIJIwLKvVbDKSBgDIZHJyMtfPFxVJGJbE6koAANLhxnwsidWVAACkQxKGJbG6EgBQBVvn5nTx7BOZNi7vxTkWIwnDkuZXVw6/971MRQJLoGwRUGytZlNXP/qYdszO6uFdu1eVRPXiHN2QhGFZtXpdZ79znAQMWAJli4Bi68WtNaluzyEJQ3LUrwQA5KU2OqKTkk5Kq761phfn6IbVkUiKFZYAgDzV6nVdf9aZ2j53Qld9dmJV/wf14hzdMBKGpFhhCQDI24ODg7pjw+mZkqdenGMxkjAk1YsVlkxnAgCqiOlIJJW1fmWr2dTDu3Yr5ubkwUGmMwEAlUEShuSy1K9sTU3ryV/96mnTmatJ5Ci9BAAoGpIwFNr8dGacOLGq6cxeLQwgkQMA9BpJGAot83Rml4UBeUyJksQBABZLmoTZvkTSP6r9f+DnIuJjKa+HasoynZl1JE3KPiXKfW0AgG6SJWG2ByR9RtKfSTosadr2rRHx36muCSyWdSRN6sGU6NS0Ym5OeuqpXO9ry3qOvD8PAFWTciRsVNL9EfGAJNm+RdJlkkjCsKayjKTNfz5LIlcbHZEHB3O9ry3rObKO5jEaCADPlDIJO0fSTxe8PizpFQmvBySTaUq0B/e1DZ52WraRtIznyDqa14vRQABYrcnJyUKcY7GUSZi7HItnvMkelzQuSZs3b04YDpCfrPe1ZRlJ68U58v48AFSRI56RF/XmxPYfSfrbiLi48/paSYqI65b6TKPRiJmZmSTxAGVW1XvCbB+IiMaKT1ZA9F9A/8nah6VMwtZJ+rGkP5X0v5KmJV0eEfcs9Rk6MaC/kIQBKLOsfViy6ciIeNL2eyTdofYWFXuWS8AAoEhsv0nSGyS9QNJnIuIbOYcEoGKSFvCOiNsjYntEvCQiPpryWgAwz/Ye20dsH1x0/BLb99q+3/Y1y50jIv4lIq6U9JeS/jxhuAD6FDvmA6iiGyV9WtIX5g8stXeh2iP1i+9V3R0RRzrPP9T5HAD0FEkYgMqJiLtsb1l0uOvehZ3FQpcuPodtS/qYpK9HxPe7XYfV3QCySDodCQAF0m3vwnOWef9fS3qNpLfYfle3N0TEREQ0IqIxPDzcu0gB9IVCjYQdOHDgmO3/WXT4bEnH8oinh6rQBol2FEkV2iBJ56/htZ7V3oW/+UXEpyR96tmevML9l1SNdlShDRLtKJpMfVihkrCIeMZXSdszZV/CXoU2SLSjSKrQBqndjjW83GFJ5y54vUnSI706eVX7L6ka7ahCGyTaUTRZ+zCmIwH0i2lJ59neantQ0lsl3ZpzTAD6GEkYgMqxvU/SdyWdb/uw7Ssi4klJ83sXHpL0ZfYuBJCnQk1HLmEi7wB6oAptkGhHkVShDVKidkTEziWO3y7p9hTXXAJ/TsVRhTZItKNoMrUjWdkiAAAALI3pSAAAgBwUNglbSXmRvHUrkWL7LNt32r6v8/PMBb+7ttOue21fnE/UT2f7XNvfsn3I9j22r+ocL1s7hmxP2f5Bpx0f6RwvVTuk9g7vtpu2b+u8LmMbHrL9X7bvnl9FVMZ2rBT919qrQh9Wpf5Log97Vu2IiMI91C4j8hNJL5Y0KOkHki7IO65l4n2VpJdLOrjg2D9Iuqbz/BpJf995fkGnPc+RtLXTzoECtOFFkl7eeX6GpB93Yi1bOyxpQ+f5eknfk/TKsrWjE9v7JN0s6bYy/p3qxPaQpLMXHStdO1bYZvqvfNpR+j6sSv1XJz76sFO0o6gjYb8pLxIRc5JukXRZzjEtKSLukvToosOXSbqp8/wmSW9acPyWiPh1RDwo6X6125uriPhZdEqzRMTjaq8eO0fla0dExGzn5frOI1SydtjeJOkNkj634HCp2rCMqrRjKfRfOahCH1aV/kuiD3u27ShqErbS8iJF9MKI+JnU7hwkvaBzvPBtc7vmXl3tb2Gla0dnCPxuSUck3RkRZWzH9ZI+IOmpBcfK1gap/R/IN2wfcLvOolTOdqxEFdpR6j+jMvdhFem/JPow6Vm0o6hbVKyovEjJFLpttjdI+oqkqyPil3a3cNtv7XKsEO2IiJOSXmb7eZK+Zvuly7y9cO2wfamkIxFxwPbYs/lIl2OF+LOQdFFEPGL7BZLutP2jZd5b5HasRFXa0U3h21b2Pqzs/ZdEH7bIsu0o6khY0vIia+T/bL9Ikjo/j3SOF7Zttter3Xl9KSK+2jlcunbMi4hfSJqUdInK1Y6LJL3R9kNqT2W92vYXVa42SJIi4pHOzyOSvqb20Hzp2rFCVWhHKf+MqtSHlbj/kujD5p2yHUVNwqpQXuRWSe/oPH+HpH9dcPyttp9je6uk8yRN5RDf07j9dfHzkg5FxCcX/Kps7RjufIOU7edKeo2kH6lE7YiIayNiU0RsUfvv/jcj4u0qURskyfbpts+Yfy7ptZIOqmTtWAX6rxxUoQ+rQv8l0YetqB1rtcJgpQ9Jr1d7dctPJH0w73hOEes+ST+TdELtTPgKSc+X9O+S7uv8PGvB+z/Yade9kl6Xd/ydmP5Y7WHTH0q6u/N4fQnb8QeSmp12HJT04c7xUrVjQWxj+u3KolK1Qe3VgT/oPO6Z/3dctnassu30X2vfjtL3YVXrvzrx0Yct82DHfAAAgBwUdToSAACg0kjCAAAAckASBgAAkAOSMAAAgByQhAEAAOSAJAy5sv0823+VdxwAsFL0X8iKJAx5e54kOjEAZUT/hUxIwpC3j0l6ie27bX8872AAYAXov5AJm7UiV7a3qL2b8nJFagGgcOi/kBUjYQAAADkgCQMAAMgBSRjy9rikM/IOAgBWgf4LmZCEIVcR8XNJ/2H7IDe2AigT+i9kxY35AAAAOWAkDAAAIAckYQAAADkgCQMAAMgBSRgAAEAOSMIAAAByQBIGAACQA5IwAACAHJCEAQAA5OD/AeFvPJKWmQ2gAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(1, figsize=(10, 3.5))\n", "ax1 = fig.add_subplot(1,2,1)\n", "ax2 = fig.add_subplot(1,2,2)\n", "ax2.set_yscale(\"log\")\n", "for ax in [ax1, ax2]:\n", " ax.errorbar(t, N, yerr=dN, xerr=None, fmt='oC3', ecolor='k', ms=3)\n", " ax.set_xlim(-10, 500)\n", " ax.set_xlabel('t')\n", " ax.set_ylabel('N')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The semi-log plot suggests we can use linfit to fit the data by taking the logarithm of the $y$ data. Taking the logarithm of the exponential fitting function gives $$\\ln N = -\\frac{t}{\\tau} + \\ln N_0\\;.$$ Defining $y=\\ln N$, $a=-1/\\tau$, and $b=\\ln N_0$, the equation takes the form $y = at+b$ and can be fit using linfit.\n", "\n", "The uncertainties $\\Delta y$ are related to $\\Delta N$ by taking the differential of the tranformation $y=\\ln N$:\n", "$$\n", "\\begin{align}\n", " \\Delta y &= \\left(\\frac{\\partial y}{\\partial N}\\right)\\Delta N \\\\\n", " &= \\frac{\\Delta N}{N}\n", "\\end{align}\n", "$$\n", "To fit the data, we tranform the $N$ and $\\Delta N$ data:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "y = np.log(N)\n", "dy = dN/N" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we perform the fit on the tranformed data:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "fit, cvm, info = linfit(t, y, sigmay=dy, relsigma=False, return_all=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extract $\\tau$ and $N_0$ from fit of transformed data" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "a, b = fit[0], fit[1]\n", "tau = -1.0/a\n", "N0 = np.exp(b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extract the uncertainties in the fitting parameters $a$ and $b$." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "dfit = [np.sqrt(cvm[i,i]) for i in range(2)]\n", "da, db = dfit[0], dfit[1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the uncertainties in $\\tau$ and $N_0$ from the transformation equations:\n", "$$\n", "\\begin{align}\n", " \\tau=−1/a\n", " \\quad &\\Rightarrow \\Delta \\tau = \\left|\\frac{\\partial \\tau}{\\partial a}\\right|\\Delta a\n", " \\quad \\Rightarrow \\Delta \\tau = \\frac{\\Delta a}{a^2} \\\\\n", " N_0=e^b\n", " \\quad &\\Rightarrow \\Delta N_0 = \\left|\\frac{\\partial N_0}{\\partial b}\\right|\\Delta b\n", " \\quad \\Rightarrow \\Delta N_0 = e^b\\Delta b \\\\\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "dtau = da/(a*a)\n", "dN0 = np.exp(b)*db" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the data and fit on linear and semi-log plots" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAADxCAYAAACDDGAhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxVVZ7v/c8vw0mADBAIkwSSMMogQSNooUIhKggIWk5oKSKCCVq3qrvrdunt7qt966mnurrv7af6tpUgkyOCMwIyOYHggGATxoAaEgmgDEIIYQpJ1vNHQipAGEJyss85+b5fr/3inH3W3vu3CKz8ztprr2XOOURERESkcYV5HYCIiIhIU6QkTERERMQDSsJEREREPKAkTERERMQDSsJEREREPKAkTERERMQDSsJEREREPKAkTERERMQDSsJERM5iZqlmNsvM3vI6FhEJXUrCRKRJMLPZZrbPzDaftX+EmW03s+/M7CkA59wO59wkbyIVkabCb0mYmfU0s5waW7GZ/cZf1xMRuYgXgRE1d5hZOPAXYCTQGxhvZr0bPzQRaYoi/HVi59x2IA2qG7rdwLv+up6IyIU45z41s+Szdg8EvnPO7QAws3nAWGBr40YnIk2R35Kws9wM5Dnnvr9QoTZt2rjk5OTGiUhEPPf1118fcM4lehjCFUBhjfe7gEFm1hr4AzDAzJ52zv2xtoPNbAowBaBFixbX9OrVy9/xikgAqW8b1lhJ2P3A3No+qNmIde7cmXXr1jVSSCLiNTO74Bezxgihln3OOfcTkHGxg51z04HpAOnp6U7tl0jTUt82zO8D883MB9wBvFnb58656c65dOdcemKil1+IRaQJ2gUk1XjfCdhTlxOY2Rgzm3748OEGDUxEQl9jPB05Evgv59zeRriWiEhdrAW6m1lK1RfG+4EFdTmBc26hc25KfHy8XwIUkdDVGEnYeM5zK1JEpLGY2VzgC6Cnme0ys0nOuTLgSWAZkAu84Zzb4mWcItJ0+HVMmJk1B24BHvfndURELsY5N/48+xcDixs5HBER//aEOeeOOedaO+c0WEJEQpLGhInI5dKM+SIi9aAxYSJyuZSEiYiIiHhASZiISD3odqSIXC4lYSIi9aDbkSJyuZSEifhJQUEBt99+O61ataJ9+/Y8+eSTlJWVAXDy5EkmTZpEly5diI2NZcCAASxZsuSC5xs6dCjR0dHExMQQExNDz549GyTO5557jvT0dKKionjkkUcarKyIBDe1Yf6nJEzET6ZOnUrbtm354YcfyMnJYeXKlWRlZQFQVlZGUlISK1eu5PDhw/z+97/n3nvvpaCg4ILnfO655ygpKaGkpITt27dfNIZnn32WZ5999oJlOnbsyD/+4z/y6KOPXvR8dSkrIsFNbZj/KQkTAa688kpiYmLw+Xz4fL7qb2q5ubmXfc78/HzuvfdeoqOjad++PSNGjGDLlsp5QFu0aMGzzz5LcnIyYWFhjB49mpSUFL7++uuGqtIlu+uuuxg3bhytW7du0LIi0njUhgVnG6YkTATIzc2lpKSEhx56iGeeeab6m9qVV15ZXWb06NG0bNmy1m306NHnnPPXv/418+bN49ixY+zevZslS5YwYsSIWq+/d+9evvnmG/r06XPBOJ9++mnatGnD4MGDWbFiRb3qLA3j9MD8AwcOeB2KNGFqw4KTX2fMFwk2Gzdu5I477qj1s0WLFtXpXEOGDGHGjBnExcVRXl7OhAkTGDdu3DnlTp06xYMPPsiECRPo1avXec/3pz/9id69e+Pz+Zg3bx5jxowhJyeHrl27nlFu9OjRrF69GoATJ04A8Oc//xmAG264oc71kAtzzi0EFnbs2HHy4sWLGT58OD6fz+uwpIlSGxZc1BMmUqWiooKtW7fSt2/fBjnXbbfdxl133cXRo0c5cOAAhw4d4ne/+9055R566CF8Ph/PPffcBc85aNAgYmNjiYqKYsKECQwePJjFi89dbWfRokUUFRVRVFTEU089xVNPPVX9PpQar0ATExPD2rVrmTZtGvlffkneqNHk9ulL3qjRlBYWeh2eNAFqw4KPkjCRKjt37qSiooLU1NRaPx85cmT1OIuzt5EjR55R9uDBgxQWFvLkk08SFRVF69atmThx4hkNjnOOSZMmsXfvXt5++20iIyPrFK+Z4Zyre0XFL+Li4pgwYQIAe//brzm5YweUl1Oan09hRqbH0UlToDYs+Oh2pEiV4uJiWrRoQWlpKVFRUed8frHHr2tq06YNKSkpZGdn89vf/paSkhJeeukl+vfvX10mMzOT3NxcPvzwQ5o1a3bB8xUVFbFmzRqGDBlCREQEr7/+Op9++ml1F/35XOypIqh8yqmsrIzy8nLKy8s5ceIEERERRESc2zzUpWxTlJycTEZGBt9Nex47/culooLSizwxJtIQ1IYFYRvmnAuY7ZprrnEiXiktLXXDhg1zzZo1c7m5ufU+3/r1692QIUNcy5YtXevWrd3dd9/t9u7d65xzrqCgwAEuKirKtWjRonp79dVXq48fMWKE+8Mf/uCcc27fvn0uPT3dxcTEuPj4eDdo0CC3fPnyWq87YsSIM85ZcxsxYsQ55Z955hkHnLE988wztcZxsbJ1BaxzAdD21GcDxgDTu3XrVl2v724f5bb2utJt7dnLbe7Zy20YMtSVlpZe9t+TyKVQGxZ8bZi5AOoKTE9Pd+vWrfM6DBFpJGb2tXMu3es4GkLN9qu0sJDCjExKCwoobdOGDwZeS3SXLowbN45OnTp5HKmINJT6tmG6hyAi0sB8SUl0ff+vA4hj8vJYsGABs2fP5mc/+xlDhw7VLVwR0cB8ERF/69q1K1OnTiUtLY3PPvuM6dOns3v3bq/DEhGPKQkTEWkEUVFR3HHHHTz44IOcOHGCWbNm8dFHH1WvxSciTY9fkzAza2lmb5nZNjPLNbPr/Xk9EZFA161bN6ZOnUr//v1ZvXo1M2bM4IcffvA6LBHxgL97wv4DWOqc6wX0By5/ESsRkQB0etmiw4cPX/Ix0dHRjB07lvHjx3Ps2DFmzJjBJ598Qnl5uR8jFZFA47ckzMzigJuAWQDOuVLnXJG/ridy6NAhzIzrrz+zw/Xxxx/nb/7mb/xyzaFDhxIdHV094WHPnj3PW/bgwYPceeedtGjRgi5duvDaa6/V6fPLVdfzXqh8QUEBt99+O61ataJ9+/Y8+eSTTf52mnNuoXNuSnx8fJ2P7dGjB1OnTqVfv358+umnzJgxgx9//NEPUUqgU/t1edetS/lAbL/82ROWCuwHXjCz9WY208xanF3IzKaY2TozW7d//34/hiOhLicnh/bt27N169Yzbu/k5OSQlpbmt+s+99xz1Yvlbt++/bzlnnjiCXw+H3v37mXOnDlkZmayZcuWS/68Ns8+++xFJzOs63kvVH7q1Km0bduWH374gZycHFauXElWVtYFry8X1qxZM+68807uu+8+SkpKmDFjBitXrlSvWBOj9uvyrluX8oHYfvkzCYsArgaynXMDgKPAU2cXcs5Nd86lO+fSExMT/RiOhLqcnBzS09O55ZZbWLBgAQDl5eVs2rSJAQMGeBrb0aNHefvtt/n9739PTEwMN9xwA3fccQevvPLKJX3ur+vWtXx+fj733nsv0dHRtG/fnhEjRly0oZVL06tXL6ZOnUqfPn1Yu2gR62+8idzefbT2ZBOh9qvu161r+UBsv/yZhO0Cdjnn1lS9f4vKpEzEL9avX09aWhrjxo1j/vz5AGzbto3y8nKuvPLKSzrH6NGjadmyZa3b6NGjaz3m6aefpk2bNgwePJgVK1bUWuabb74hPDycHj16VO/r379/dQNwsc8vV13Pe7Hyv/71r5k3bx7Hjh1j9+7dLFmyhBEjRtQrRvmr5s2bc9dddzEqJ4fon36CigpO7tihtSebALVfdb9uXcsHYvvlt9kCnXM/mlmhmfV0zm0Hbga2+ut6Ijk5OYwdO5Zhw4aRkZHBkSNHyMnJoW/fvkRGRvIP//APfPrpp7Rr146XX36Z5s2bn3OORYsW1XLm8/vTn/5E79698fl8zJs3jzFjxpCTk0PXrl3PKFdSUsLZY4bi4+M5cuTIJX1e0+jRo1m9ejUAJ06cAKhef+2GG244ow51Oe+llB8yZAgzZswgLi6O8vJyJkyYwLhx42o9l9TDnh8INwPAnOPkjh3s378f3S0IXWq/mmb75e+nI38FzDGzjUAa8P/6+XrSRJ08eZLc3FzS0tJo1aoVAwcOZMmSJdXfLjdv3kxeXh6rVq1i+PDhzJ49u0GuO2jQIGJjY4mKimLChAkMHjyYxYsXn1MuJiaG4uLiM/YVFxcTGxt7SZ/XtGjRIoqKiigqKuKpp57iqaeeqn5/diNcl/NerHxFRQW33XYbd911F0ePHuXAgQMcOnSI3/3ud+f525HL5UtOhrDK5tmZcSQ+nueff57PPvuMiooKb4OTBqf2q+m2X35NwpxzOVXjva5yzo1zzh3y5/Wk6dq8eTPNmjUjNTUVoLpLf/369QwYMIBVq1YxcuRIAEaOHFn9TexsI0eOrH5S6Ozt9PEXYmbUth5rjx49KCsr49tvv63et2HDBvr06XNJn1+uup73QuUPHjxIYWEhTz75JFFRUbRu3ZqJEyfW2mhL/SRNy8aXkgLh4USlptLzxRfo0aMHH374IS+88AIHDhzwOkRpQGq/atck2q/6rP7d0Ns111xz2SuZS9M2Y8YMd8MNN1S/LygocPHx8a5Vq1Zu1apV7g9/+IN79913nXPOHTp0yN1yyy31vuahQ4fc0qVL3fHjx92pU6fcq6++6po3b+62bdtWa/n77rvP3X///a6kpMStXr3axcXFuc2bN1/y55errue9UPmUlBT3xz/+0Z06dcodOnTIjRs3zj3wwAOXHRuwzgVA21OfDRgDTO/Wrdtl/z1cioqKCrdx40b3pz/9yf3+9793n332mSsvL/frNaVxqP06v0Buv5yrfxvmeQNWc1MSJpfriSeecE8++eQZ+/r37+/MzBUXF7usrCz34osvOucqG7j77ruv3tfct2+fS09PdzExMS4+Pt4NGjTILV++/IwyI0aMcH/4wx+cc8799NNPbuzYsa558+YuKSnJzZkz54yyF/u85jlbtGhR6zZixIhzyl/svDVjvFj59evXuyFDhriWLVu61q1bu7vvvtvt3bu3bn9xNYRCEnZ6a6z268iRI27u3Lnu2WefdbNmzXIHDhxolOuK/6j9Cs72y7n6t2FWeY7A0L9/f7dhwwavw5AQtGnTJv74xz/y2muvMX36dE6ePMmvfvUrr8Nq8szsa+dcutdxNIT09HS3bt26RrmWc45NmzaxZMkSysrKuPnmmxk0aBBWNZhfQovar8BV3zbMb09HXo7TT0qINLR+/frRpUsXbrzxRtq2bcvLL7/sdUgil83MuOqqq0hJSWHhwoUsW7aM3Nxcxo4dS0JCgtfhSQNT+xW6AqonLDk52RUUFHgdhog0EvWE1Z9zjg0bNrB06VIqKioYPnw41157rXrFRBpBfdswf09RUSelpaUcPXrU6zBERIKGmZGWlsbUqVPp0qULS5Ys4eWXX+bQIT2MLhLoAioJA/juu++8DkFEJOjExcXxwAMPMGbMGPbs2cNL//qvbLl5OLl9+mrpI5EAFVBJWHh4+Bnze4iIyKUzM66++mqmTp3KTatW43btgvJySrX0kUhACqgkLCoqiry8PM0ILSJSD/Hx8TQ7dKh66SOc42R+PoE0BlhEAiwJi46O5sSJE+zatcvrUEREglrNpY8qzCiOiWHOnDkcPnzY28BEpFpAJWFRUVGYmW5JiojUU82lj6JTU4l45n+yc+dOsrOzWb9+vXrFRAJAQM0TZmZ07tyZb7/9lptvvtnrcERELsrMxgBjunXr5nUoZ/AlJdH1/TMXRE69/noWLFjAggULyM3NZcyYMeddDFlE/C+gesIAunfvzt69e89ZCV1EJBA55xY656bEx8d7HcpFJSQkMGHCBEaMGEF+fj5ZWVls2LBBvWIiHgnIJAw0VYWIiD+YGYMGDSIjI4PExETmz5/P66+/TklJidehiTQ5AZeEJSYmEhcXp3FhIiJ+1Lp1ax555BFuvfVW8vLyyMrKYtOmTeoVE2lEAZeEmRndu3dnx44dDBkyhKFDh3odkohISAoLC+P666/n8ccfp3Xr1rzzzju8+eabWrlEpJH4NQkzswIz22RmOWZ2yYuq9ezZk8iDB/mn/Qd47se9mu1ZRMSP2rRpw8SJExk+fDjffPMNf/nLX9iyZYvXYYmEvMboCfu5cy6tLgtcpqamcuOnq2hfUUE4UJqfr9meRUT8KCwsjMGDB/P444/TqlUrlrz4IutvvInc3n30RVjETwLudiRULl8UW1xM+OkdFRWUFhR4GJGISNOQmJjIpEmTuG3d10Tu2wcVFZzUskcifuHvJMwBy83sazObUpcDrVMnyk8PEA0Lq5z9WURE/C4sLIyIffuqlz0y5zi5YwfHjh3zODKR0OLvJGywc+5qYCTwhJnddHYBM5tiZuvMbN3+/fur9ydPf55dzlEO+FJSSJqW7edQRUTktJrLHjkziuNiycrKYtu2bd4GJhJC/JqEOef2VP25D3gXGFhLmenOuXTnXHpiYmL1/mYpKRz/07/w/sRHSF7wHr6kJH+GKiIiNdRc9igqNZWk7GxiY2N5/fXXeeeddzh+/LjXIYoEPb8tW2RmLYAw59yRqte3Av+rLufo3bs3GzduJD8/n0BbEkREJJTVtuzRY/37s2rVKlatWkV+fj5jxoyhR48eHkUoEvz82RPWDlhtZhuAr4D3nXNL63KCrl274vP52Lp1q18CFBGRSxceHs7QoUN57LHHaN68OXPnzmX+/PmcOHHC69BEgpLfkjDn3A7nXP+qrY9z7g91PUdERAQ9e/Zk27ZtlJeX+yNMERGpow4dOjBlyhRuvPFGNm7cSFZWllY5EbkMATlFRU29e/fm+PHjFGiKChGRgBEeHs6wYcN47LHHiI6O5rXXXmPBggXVvWKlhYXkjRpNbp++mmdM5DwCPgnTLUkRaWxm1sLMXjKzGWb2oNfxBLKOHTsyZcoUBg8eTE5ODtnZ2eTl5VGYkcmJvDwoL9eE2yLnEfBJWGRkJD169GDbtm1UVFR4HY6IBCkzm21m+8xs81n7R5jZdjP7zsyeqtp9F/CWc24ycEejBxtkIiIiGD58OI8++ig+n49XX32Vk/n5f/0Fowm3RWoV8EkYVN6SPHbsmG5Jikh9vAiMqLnDzMKBv1A5l2FvYLyZ9QY6Aafvn2lA6iXq1KkTU6ZM4frrr6c4JoaKqsleNeG2SO2CIgnr1q0bkZGRuiUpIpfNOfcpcPCs3QOB76oeJCoF5gFjgV1UJmIQJO1koIiMjOTWW2+lw3P/ybGWLakwozQxkfb/+X+9Dk0k4ARF43L6lmRubq5uSYpIQ7qCv/Z4QWXydQXwDvALM8sGFp7v4POt+CHQZeBA0j5dya4//3/MH/ZzZi5cqLsZImcJiiQM/npL8vvvv/c6FBEJHVbLPuecO+qcm+icy3TOzTnfwedb8UMqRUZGctttt/HII49gZrz00kssWbKE0tJSr0MTCQhBk4R1796dyMhINm/efPHCIiKXZhdQc020TsAej2IJWV26dCEjI4OBAwfy1VdfMW3aNHbu3Ol1WCKeC5okLDIykt69e7NlyxZOnTrldTgiEhrWAt3NLMXMfMD9wIK6nMDMxpjZ9MOHD/slwFDh8/kYOXIkEyZMwDnHCy+8wLJly9SeS5MWNEkYQFpaGidPniQ3N9frUEQkyJjZXOALoKeZ7TKzSc65MuBJYBmQC7zhnNtSl/M65xY656bEx8c3fNAhKDk5mczMTNLT0/nyyy95/vnnKdRErtJEBVUS1qVLF1q2bElOTo7XoYhIkHHOjXfOdXDORTrnOjnnZlXtX+yc6+Gc63o5y6tJ3fl8PkaNGsVDDz1EWVkZL7zwAh988AFlZWVehybSqIIqCTMz0tLSyM/Pp6ioyOtwRER0O7IeUlNTyczMZMCAAXz++ec8//zz7N692+uwRBpNUCVhAP379wdgw4YNHkciIqLbkfUVFRXFmDFj+OUvf0lpaSmzZs3iww8/VK+YNAlBl4S1bNmSlJQUcnJycM55HY6IiDSArl27kpmZSf/+/fnss8+YPn06e/boQVUJbUGXhEHlAP2ioiLNGSYiEkKio6MZO3YsDzzwACdOnGDmzJl8/PHHlJdr5SgJTUGZhF155ZVERUVpgL6IeE5jwhpe9+7deWzMGK57bS5tM6ey/sab2PX1116HJdLggjIJi4yMpE+fPmzdupWTJ096HY6INGEaE+Yf+3/zN1wBhJvR7OBBdk2dyooVK9QrJiHF70mYmYWb2XozW9SQ501LS+PUqVNa1FtEJASVFhRU/4IKA+KKj7By5UpmzpzJ3r17vQxNpME0Rk/Yr6mcBLFBderUidatW+uWpIhICPIlJ1Nx+k1YGFGpqdx3330cOXKE6dOn8+mnn6pXTIKeX5MwM+sEjAJm+uHcpKWlsXPnTg4ePNjQpxcRuSQaE+YfSdOyie7aFcLD8aWkkDQtm169ejF16lSuvPJKPvnkE2bNmsW+ffu8DlXksvm7J+zPwN/DX7/QNKT+/ftjZuoNExHPaEyYf/iSkuj6/iKu3LKZru8vwpdUuc568+bNufvuu7nnnns4fPgw06dPZ/Xq1VRU+OXXjIhf+S0JM7PRwD7n3AUfaTGzKWa2zszW7d+/v07XiI2NpWvXrmzYsEH/AUVEmpDevXszdepUevTowUcffcTs2bOp6+8QEa/5sydsMHCHmRUA84BhZvbq2YWcc9Odc+nOufTExMQ6X+Tqq6+muLiYb775pt4Bi4hI8GjRogX33HMPv/jFLzh48CDPP/88n332mb6US9DwWxLmnHu6apHcZOB+4GPn3C8b+jo9e/YkPj6eNWvWNPSpRUQkwJkZffv2ZerUqXTv3p0PP/yQF154gZ9++onSwkLyRo0mt09f8kaNprSw0OtwRc4QlPOE1RQWFsa1115LQUGBHlsWkUangfmBISYmhnvvvZe77rqLAwcOMG3aNL6Z8Agn8vKgvJzS/HwKMzK9DlPkDI2ShDnnVjjnRvvr/FdffTWRkZF8+eWX/rqEiEitNDA/cJgZ/fr1Y+rUqaSmpmI//PDXX3IVFZQWFHgYnci5gr4nDKBZs2ZcddVVbNq0iWPHjnkdjoiIeCg2Npb7778frriCcucAcGb4kpO9DUzkLCGRhAEMGjSI8vJyvtb6YiIiTZ6Z0f2F2USlpuLMKI6N5atbhnPo0CGvQxOpFjJJWGJiIqmpqaxdu1azKIuICL6kJLovWcyVW7cQ9fw0vj9xguzsbNauXYur6iET8VLIJGFQ2Rt25MgRcnMbfJUkEREJUmbGgAEDyMzMpHPnzixevJhXXnmFoqIir0OTJi6kkrDu3buTkJCg6SpEpNHo6cjgER8fz4MPPsjo0aPZvXs32dnZfP311+oVE8+EVBJmZgwcOJBdu3axe/dur8MRkSZAT0cGFzPjmmuuITMzkyuuuIJFixYxZ84clESLF0IqCQNIS0vD5/OpN0xERM6rZcuWPPTQQ4waNYqdO3eSnZ3N+vXr1SsmjSrkkrCoqCjS0tLYsmULR44c8TocEREJUGZGeno6mZmZdOjQgQULFvDaa69RXFzsdWjSRIRcEgaVA/QrKipYu3at16GIiEiAa9WqFQ8//DAjR47k+++/Jysri5ycHPWKid+FZBKWkJDAlVdeyVdffcXx48e9DkdERALc6THFGRkZtGvXjvfee4958+bpjor4VcSFPjSz/3mBj51z7vcNHE+Duemmm8jNzWXNmjUMHTrU63BE5BIFc7sjwS8hIYFHHnmENWvW8NFHH5GVlcXIkSPp168fZuZ1eBJiLtYTdrSWzQGTgN/5N7T6STh1iuvmzqPN4xl8d/soSgsLvQ5JRC5N0LY7EhrMjOuuu46MjAwSExN59913eeONNygpKfE6NAkxF+wJc879n9OvzSwW+DXwKDAP+D/nOy4QFGZkcgUQbkZpfj6FGZl0fX+R12GJyEUEW7tjZmOAMd26dfM6FGlgrVu35uWXX6a4uJjw8HCysrK4/fbb6dOnj3rFpEFcdEyYmSWY2f8DbKQyabvaOfc759w+v0dXD6UFBYRXvTbnKC0o8DIcEamDYGp3NE9Y6CotLOSf9h/g5WPHuXf1Z3QIC+Ptt9/mrbfe4ujRo16HJyHggkmYmf0bsBY4AvRzzj3rnAuK1U99yclUVL2uAMratvUyHBG5RMHc7khoKczIpH15OeFA+c6d/OyTFQwfPpzt27eTlZXF1q1bvQ5RgtzFesL+DugI/COwx8yKq7YjZhbQE6kkTcsmumtXCA/nZJs2rBw8mNLSUq/DEpGLC9p2R0JLaUHBX39JVlRQWlDA4MGDmTJlCi1btuTNN9/krbfe4tixY16GKUHsgkmYcy7MOdfMORfrnIurscU65+IaK8jL4UtKouv7i7hyy2bavT6PnyIjNG+YSBAI5nZHQkvNOyqEheFLTgagbdu2TJo0iWHDhpGbm0tWVhbbtm3zKkwJYn6bJ8zMos3sKzPbYGZbzOyf/XWti7niiivo1q0bn3/+uXrDRETkkiRNy+bH8HDKAV9KCknTsqs/CwsL48Ybb2TKlCnExsby+uuv884772huSqkTf07WehIY5pzrD6QBI8zsOj9e74Juuukmjh07xrp167wKQUREgogvKYmbt2ym77Zcur6/CF9S0jll2rVrx2OPPcbQoUPZsmULWVlZbN++3YNoJRj5LQlzlU5PqhJZtXm2BkRSUhKpqanqDRMRkQYVHh7OkCFDmDx5Mi1atGDevHnMnz+fEydOeB2aBDi/LltkZuFmlgPsAz5wzq2ppcwUM1tnZuv279/vz3D4+c9/ztGjR1m9erVfryMiIk1P+/btmTx5MjfddBMbN24kKyuLb7/91uuwJID5NQlzzpU759KATsBAM+tbS5npzrl051x6YmKiP8OhU6dO9OvXjy+++ILDhw/79VoiItL0hIeH8/Of/zK088AAAB9ySURBVJzHHnuMZs2a8dprr/Hee++pV0xq1SgLeDvnioAVwIjGuN6F3HzzzQB89NFHHkciIiKhqmPHjkyePJl33nmH3/zmN2RnZ5OXl+d1WBJg/Pl0ZKKZtax63QwYDnj+DG98fDzXX389mzZtYvfu3V6HIyJBzszGmNl09a7L2Sp++IF/LytnYVQ0N7z5Ju9Om8bChQs5efKk16FJgPBnT1gH4BMz20jl7NcfOOcCYvHGwYMHExMTw7Jly3DOs2cFRCQEaNkiOZ+aM+7HFB3m1rXrWL9+PdnZ2ezYscPr8CQA+PPpyI3OuQHOuaucc32dc//LX9eqq6ioKH7+859TWFioZSdERMQvzp5xP3L/fiZOnEhERASvvPIK77//vp7Wb+IaZUxYIEpLS6Ndu3Z8+OGHlJWVeR2OiIiEmNpm3E9KSuLxxx/nuuuuY926dWRnZ1NQUOBhlOKlJpuEhYWFcdttt1FUVMSXX37pdTgiIhJizjfjfmRkJLfddhsTJ04kLCyMl156iSVLlqhXrAmK8DoAL6WkpNCzZ09WrVpFWloaMTExXockIiIh4vSM++fTuXNnHn/8cT766CO++uorvv32W8aOHUuXLl0aMUrxUpPtCTvtlltuoaysjA8//NDrUEREpInx+XyMHDmSCRMmAPDiiy+ydOlSTp06VV1m6NChDB061KMIxZ+afBLWunVrBg8ezIYNGzSHi4iIeCI5OZmMjAzS09NZs2YN06ZNo7Cw0OuwxM+afBIGlYt7t27dmoULF+qevIiIeMLn8zFq1CgefvhhysvLmT17Nh+99hr/tP8Az/24l7xRoylVYhZSlIQBERER3HHHHRw+fJiPP/7Y63BERKQJS0lJITMzk2uuuYYWf/4P2pWVEQ6U5udTmJHpdXjSgJSEVencuTPXXnsta9asURewiIh4KioqitGjRxNfUkK4WeXOigpKNZ1FSFESVsPNN99MXFwcCxcu1NxhIiLiOV9KSvVcYxXA0fh49uzZ42VI0oCUhNVw+pvH/v37WbVqldfhiIhIE1dzrrGwpCTW3noLM2fO5OOPP1ZnQQho0vOE1aZ79+7069eP1atX07t3b9q1a+d1SCIi0kT5kpL4fWIbAFZ8sJyUEydYtmwZq1atYvv27YwbN44OHTp4HKVcLvWE1WLEiBFER0ezYMECysvLvQ5HRESasBUrVrBixQoAoqOjGTt2LOPHj+fYsWPMnDmTFStW6HdVkFISVovmzZtz+9XXkPRv/5vcPn31WLCIiASUHj16MHXqVPr27cvKlSuZOXMmP/74o9dhSR0pCTuP6H//dzqFhVU+Frxjhx4LFmlCzCzVzGaZ2VtexyJyPs2aNePOO+/kvvvu48iRI8yYMYOVK1eqVyyIKAk7j9KCAsJPv3GO0vx8L8MRkUtkZrPNbJ+ZbT5r/wgz225m35nZUxc6h3Nuh3Nukn8jFWkYvXr1YurUqfTu3ZsVK1Ywa9Ys9u3bV/25lj0KXErCzsOXnHzGY8HHExKoqKi40CEiEhheBEbU3GFm4cBfgJFAb2C8mfU2s35mtuisrW3jhyxSP82bN+cXv/gF99xzD4cPH2b69OmsWrVKv7cCnJ6OPI+kadkUZmRWTozXvj0fXz2AktWruemmm7wOTUQuwDn3qZkln7V7IPCdc24HgJnNA8Y65/4IjL7ca5nZFGAKVE74LOK13r1706VLFxYvXszHH39M/pdf8o/79tO+ooK8UaNJmpaNLynJ6zClit96wswsycw+MbNcM9tiZr/217X8wZeURNf3F3Hlls30/vADkgcNYsWKFezcudPr0ESk7q4Aaj5ds6tqX63MrLWZTQMGmNnT5yvnnJvunEt3zqUnJiY2XLQi9dCiRQvuuece7r77bnq+/Q7ty8u17FGA8uftyDLg75xzVwLXAU+YWW8/Xs9vzIwxY8bQsmVL3n77bY4dO+Z1SCJSN1bLPne+ws65n5xzGc65rlW9ZSJBp0+fPsQVF5+57JHGNwcUvyVhzrkfnHP/VfX6CJDLBb55BrqoqCjuvvtuSkpKePvtt/X0iUhw2QXUvAfTCWiQtV/MbIyZTT98+HBDnE6kQZ297FFxbCxffPFFncaKaWC//zTKwPyq8RkDgDW1fDbFzNaZ2br9+/c3RjiXrWPHjowePZodO3awZMkSnDvvF2kRCSxrge5mlmJmPuB+YEFDnNg5t9A5NyU+Pr4hTifSoGoue+RLSWHXhIdZvnw5L774Ij/99JPX4TV5fh+Yb2YxwNvAb5xzxWd/7pybDkwHSE9PD/isZsCAARw4cIDPP/+cNm3acN1113kdkojUYGZzgaFAGzPbBTzjnJtlZk8Cy4BwYLZzbouHYYo0ijOWPVqymG7O0W3jRpYuXcq0adO4+eabGTRoEGa13bEXf/NrEmZmkVQmYHOcc+/481qNafjw4Rw8eJDly5eTkJBAjx49vA5JRKo458afZ/9iYHEjhyMSUMyM/v37k5qaysKFC1m2bBnbtm3jjjvuICEh4ZzypYWF/NP+A7QrL9fTlX7gz6cjDZgF5Drn/t1f1/GCmXHnnXfSrl073n77bfbu3et1SCLiEY0Jk0BXc+3J02JjYxk/fjxjx47lxx9/ZNq0aXz11VfnDLMpzMjU05V+5M8xYYOBh4BhZpZTtd3ux+s1Kp/Px/jx44mKimLu3LmUlJR4HZKIeEBjwiRYmRlpaWlMnTqVzp07s2TJEl5++WWKioqqy5QWFPw1UaioqJw7UxqMP5+OXO2cM+fcVc65tKotpG4FxMXFVa9kP2/ePEpLS70OSUREpE7i4uJ48MEHGTNmDHv27CErK4t169bhnDtj9RjCwvAlJ3sYaejRskX11KFDB+666y727NnD3LlzOXXqlNchiUgj0u1ICQVmxtVXX83UqVNJSkri/fff59VXXyX+X/90xtOVSdOyvQ41pCgJawC9evXizjvvpKCggNdff52ysjKvQxKRRqLbkRJK4uPj+eUvf8moUaMoLCxkxoIF/K5ZNE+2b0fX9xdpUH4DUxLWQPr168cdd9xBXl4eb731liZzFRGRoGRmpKenk5mZSceOHTlw4AB79+6luPicWaaknpSENaABAwYwcuRItm/fzjvvvKPV60VEJGi1atWKhx9+mDfeeIOHHnqIrKwscnJyNFF5A/L7ZK1NzcCBAykrK+ODDz4gIiKCcePGaRI8kRBmZmOAMd26dfM6FJEGZ2YMHDiQbt26sWDBAt577z22bt3KmDFjiI2N9Tq8oKeeMD/42c9+xs19+9Lhj/9Cbu8+5I0aTWlhoddhiYgfaEyYNAUJCQlMmDCBESNGkJ+fT1ZWFhs2bFCvWD0pCfOTDrNfIKa4GHOOkzt2sPPxDK9DEhERuWxmxqBBg8jIyCAxMZH58+fz+uuva57MelAS5ielBQWEV7025yjNz+fYsWOexiQiIlJfrVu35pFHHuHWW2/lu+++Iysri02bNqlX7DIoCfOTMya4M+NIXByzZ88+YyZiERGRYBQWFsb1119PRkYGCQkJvPPOO7z55pscPXrU69CCipIwP0malk10164QHo4vNZW2//Fnjh49yuzZs7XWpEgI0WSt0pS1adOGRx99lOHDh/PNN9+QlZXFli1bvA4raFggdR+mp6e7devWeR2G3+zdu5c5c+ZQWlrKnXfeSc+ePb0OScRTZva1cy7d6zgaQqi3XyIXs3//fubPn8+ePXvo3bs3t99+Oy1atPA6LL+qbxumnrBG1K5dOyZNmkRCQgLz5s3jk08+0T10EREJCYmJiUyaNIlhw4axbds2srKyyM3N9TqsgKYkrJHFx8czceJE+vfvz6effsrcuXM5ceKE12GJiIjUW1hYGDfeeCNTpkwhPj6eN954g7ffflsPpp2HkjAPREZGMnbsWEaOHEleXh4zZsxg3759XoclIiLSIE7f+Rk6dChbt24lKyuLbdu2eR1WwFES5pHTsxBPmDCB0tJSZs6cqYnvREQkZISHhzNkyBAmT55MTEwMr7/+Ou+++y7Hjx/3OrSAoSTMY507d2bKlCl06NCB+fPn8+abb1Z325YWFpI3ajS5ffpq1n0REQlK7du3Z+7cubz77rts3ryZrKwsvvnmmzqdY+jQoQwdOtQ/AXrIb0mYmc02s31mttlf1wgVsbGxTJgwgZtvvpnt27eTnZ3Nt99+S2FGJify8qC8nNL8fAozMr0OVUTOoikqRC7OzGjZsiWPPfYYzZs3Z+7cubz33ntNfky0P3vCXgRG+PH8ISUsLIwbbriByZMn07x5c1577TVO5uf/9QdUUUFpQYGHEYpIbbR2pMil69ChA5MnT+bGG29kw4YNZGVl8d1333kdlmf8loQ55z4FDvrr/KGqffv2TJ48meuvv57imBgqzCo/CAvDl5zsaWwiIiL1FRERwbBhw3jssceIjo5mzpw5LFiwoEn2imlMWACKiIjg1ltvpe1//JljLVtSYcaJ1q1p+W//6nVoIiIiDaJjx45MmTKFwYMHk5OTQ3Z2Nnl5eV6H1ag8T8LMbIqZrTOzdfv37/c6nICSev31DFj1KQenP8/iEbfx/Pz5rF69mvLycq9DExERuSSlhYX80/4DPPfj3nMeMouIiGD48OE8+uijREZG8uqrr7Jo0SJOnjx5SccHO78uW2RmycAi51zfSymvZT/Or6ioiKVLl7J9+3batGnDbbfdRteuXbHTtytFgpCWLRIJfXmjRnMiL6+y1ycsDF9KCl3fX3ROuVOnTvHJJ5/wxRdfEB8fz9ixY0lJSbnk472gZYuaiJYtW3L//fczfvx4ysvLmTNnDi+99BKFIfSNQEREQk9pQcElPWQWGRnJrbfeysSJEwkPD+fll19m8eLFl3x8MPLnFBVzgS+Anma2y8wm+etaTUmPHj144oknGDlyJAcOHGD27NnMmzdPM+6LiEhA8iUnU3H6zSU8ZNa5c2cyMjIYNGgQa9eupSQujupBOCH2kJpfb0fWlbrz66a0tJQvv/ySzz//nJMnT3LVVVcxePBg2rZt63VoIpckFG5HmtkYYEy3bt0mf/vtt16HIxJwSgsLWTViJO3Ky2nWtStJ07LxJSVd0rHff/89y19+haRXXqFTWBjRXVPpPG3aJR/vb/Vtw5SEhYBjx47x2WefsXbtWk6dOkWPHj342c9+RufOnTVmTAJaKCRhp6n9Ejm/07Pdr1ixos7HlpaWkp6eTnFxMX/3d3/H2LFj6dy5c8MGeJnq24ZFNGQw4o3mzZtzyy23MHjwYNauXctXX33Fiy++SKdOnRg8eDA9e/ZUMiYiIkHJ5/ORkJBA8+bNqaio4IUXXuC6665j2LBhREZGeh1evWhgfghp3rw5Q4YM4Te/+Q0jR46kfM8eih6ewNbefdj082EUaQV7EREJUtHR0WRmZpKens6XX37J888/z65duy7p2IZYe9If61cqCQtBkZGRDBw4kJu/XENMcTFhzhH2ww9smzCBt956i4KCAgLpNrSIiMil8Pl8jBo1ioceeoiysjJmz57NBx98QFlZmdehXRbdjgxhpQUFhFe9DgPijpSQl5fHli1bSEhIoF+/flx11VUkJCR4GaaIiEidpKamkpmZyfLly/n888/55ptvGDduHFdccYXXodWJesJCmC85GcKqfsRhYUSlpPC3f/u3jB07lri4OFauXMl//ud/MnPmTNasWUNJSYmn8YqIiFyqqKgoxowZw4MPPkhpaSmzZs3io48+CqpeMT0dGcJKCwspzMiktKAAX3LyOY8FFxcXs2nTJjZt2sTevXsxM7p06ULPnj3p1asXzY8cueDxIvWlpyNFpCGcOHGCZcuWkZOTQ2JiIuPGjaNjx47Vn9fn6cwLnUNTVEiD2LdvH5s3b2bbtm2cXsNz1LLlRB88WHlLM8CWipDQoCRMRBrSt99+y8KFCykpKeGGG25gyJAhhIeHKwm7FGrEAsPBgwfZtm0bLSc9htX49+HMOPnaHFJTU4mPj/cwQgkVSsJEpKEdP36cpUuXsnHjRtq1a8foa69l+8MTLmuy2NPON+GskjDxm7xRoyndsQOcw5lREh/PkpEjAGjTpg1dunShS5cudO7cWUmZXBYlYSLiL9u3b2fRokUMfuMNYg4X1+uuzvkWEddkreI3SdOyq8eERSUn0zU7iy5RUezYsYMdO3awadMmvv76awDi4+Pp3LkznTt3pmPHjrRr147w8PCLXEEk+NVYtsjrUESkhp49e9K5c2fyZ86ierryy1wA3F+LiCsJk/PyJSWd822hHdCuXTuuv/56Kioq2Lt3Lzt37mTnzp3ViRlAeHg47du3p3NUFFe89BJhe37Al5Kiwf0ScpxzC4GF6enpk72ORUTO1KxZM6JSUzmel0c4UAFUtG1LRUUFYWGXPkGELzn5zJ6wBlpEXLcjpcE45ygqKmLPnj3s3r2bPXv2cGX2NGIOHyaMyn/8xxMSKPzvv6Vdu3a0bduWxMRE4uLitKxSE6XbkSLibzXHc51q25aPBg0krnt3xo0bR9u2bet8joYcE6aeMGkwZkarVq1o1aoVffr0ASD33/539edhQPNDh8jPz2fjxo3V+yMjI2nTpk311rp1axISEmjVqhXR0dGNXQ0REQkhvqQkfp/YBqh8sjFq61bef/99pk+fztChQ/nZz3520V6xM87RgLMEKAkTv/IlJ1Oanw8VFWdMGHvs2DH27dvHgQMHqredO3dW3848rU1ZGYM+/oTmBw9S1q4dp/7+vxPfowfx8fHEx8fj8/k8qpmIiASj3r1706VLF95//30++ugjtm3bxtixY0lMTGz0WJSEiV/VHNx/esJXqFxsPDk5meSz7quXlpZy6NAhDh48yMGDB4n5H/9A+E8/YUD4jz9y7Nl/Zs6o26vLR0dHExcXR2xsLLGxsdWvY2JiiImJIfrwYYr+/nec0oSzIiJSpUWLFtxzzz1s2bKFxYsX8/zzzzNs2DCuu+66Oo0Vqy8lYeJXtQ3uv2B5n4927drRrl07AHIPHKj+LAyILylh4sSJHD58mKKiIoqLizly5AjFxcXs3bv3nKWXbnt/MTGHDxNuxokdO9h0/3i+/dWTNG/evHpr1qxZ9Xb6fXR0NBERlf89LrbywMXU93gREWl4Zkbfvn1JTk5m0aJFfPDBB+Tm5jJu3Dhat27dKDH4NQkzsxHAfwDhwEzn3L/483oSes6+nelLSaFz587nLV9eXs7Ro0cpKSmhpKSEstffqB70H+YczQ4epKioiN27d3Ps2DEqKirOe67w8HCio6MZ8vY7NC8qIsw5Tu7YwdYHf8mP/+NpoqKiiIqKwufz4fP5zngdGRlZ/ecPGRmcyi+ofKw5P5/CjEytPCAiEiBiYmK477772LRpE0uWLGHatGkMGzaMQYMG+b1XzG9JmJmFA38BbgF2AWvNbIFzbqu/rimh53y3M88nPDycuLg44uLiAMhLSTlnTFpGRgZQ+TTnyZMnOXHiBMePH+f48eMcO3aM48ePc+LEieqtRVFR9coB5hyR+/ezfv16SktLL6kOd+/IJ+z0U8gVFZzYsYPnnnuOiIgIIiMjiYyMJCIi4pwtPDyciIgIfAcP0jIrm/C9e6no0J5Tv/3vhHXsQHh4+DlbWFjYOa/DwsKo+OEHDv72t5QVfE9kcjIdn3uOqC6dCQsLu6QnU9UbKCKhzMy46qqrSElJYdGiRSxfvrx6rFhCQoLfruvPnrCBwHfOuR0AZjYPGAsoCZNLVtfbmWe7UBJnZkRHRxMdHU3Lli3Pe468ac+fkchFp6Tw9NNP45yjtLS0ejt58uQZ70+dOsWpU6eoWLkS++FHrGrlgVOJibRv356ysrLqMidOnKCsrOycrby8nFsWLsIOH8bMsN17OPnMMyyrMS7uUtS8LXtyxw423n//GecICwurdTMzwsLCGPz6G9W9gSd27GDjfffz9YSHK2OqKnOh172ys4k+8BNW1Zu45YEHyfv1f6v7D1RExI9iY2O5//772bBhA0uXLiU7O5vhw4czcOBAv1zPn0nYFUBhjfe7gEFnFzKzKcAU4IK3mUQuR32TODh/Imdm1bckL6T0pZfOWHmg27Rs0urQC5Q773U4fUuVynFxv/rVrygvLz9nq6ioOOd1RUUFvpq3ZYG4khKGDRtGRUVFdTnnXPX705tzDufcGb2BYc7R/NAh4uLiqj+vWdY5V32+09vpBAwqexOjDhxg7969dfxJiIj4n5mRlpZGamoqCxcuZOnSpeTm5lJWVlY9Vrih+DMJq+0exzkzwzrnpgPToXKyQz/GI3JZ6pvI1fv4WsbF1bV7vLbbsjfeeOOlHz99xjnHjx8//tKPn//eOcc/8cQTPPnkk3Wqh4hIY4mLi+OBBx4gJyeHZcuWsWfPHlq1aoVzrsEmGPfniLNdQM2v+52APX68nkhISpqWjS8lBcLDq5d+auxzeH28iEh9rFixghUrVtT5ODNjwIABZGZm8s///M/ccccdvPLKKxQVFTVIXH5btsjMIoBvgJuB3cBa4AHn3JbzHaNlP0SaFi1bJCLBwjnHf/3Xf7F8+XIAbr31VtLT0wNz2SLnXJmZPQkso3KKitkXSsBEREREApWZcc0119C1a1cWLFjAokX1n2rIrxNgOOcWO+d6OOe6Ouf+4M9riYg0JDMbZ2YzzOw9M7vV63hEJDC0bNmShx56iNtvr9tT6rVpvLn5RUQaiZnNNrN9Zrb5rP0jzGy7mX1nZk9d6BzOufnOucnAI8B9fgxXRIKMmXHttdfW+zxatkhEQtGLwHPAy6d3nG8CaSqHS/zxrOMfdc7tq3r9j1XHiYg0KCVhIhJynHOfmlnyWbtrnUDaOfdHYPTZ57DKZ9D/BVjinPsv/0YsIk1RQCVhX3/99QEz+/6s3W2AA7WVDyKhUAdQPQJJKNQBoGcjXuuSJpCu4VfAcCDezLo556adXaDmZNNAiZltP6tIqPycQqEeoVAHUD0CTb3asIBKwpxziWfvM7N1wf4IeyjUAVSPQBIKdYDKejTm5WrZd945epxz/xf4vxc6Yc3Jpmu9YAj9nIK9HqFQB1A9Ak192zANzBeRpkITSItIQFESJiJNxVqgu5mlmJkPuB9Y4HFMItKEBUMSdt6u/iASCnUA1SOQhEIdwE/1MLO5wBdATzPbZWaTnHNlwOkJpHOBNxphAmn9nAJHKNQBVI9AU696+G3ZIhERERE5v2DoCRMREREJOQGbhNVlZmuv1TY7t5klmNkHZvZt1Z+tanz2dFW9tpvZbd5EfSYzSzKzT8ws18y2mNmvq/YHWz2izewrM9tQVY9/rtofVPWAyslFzWy9mS2qeh+MdSgws01mlnP6KaJgrEddqf1qfKHQhoVS+wVqwy6pHs65gNuonME6D0gFfMAGoLfXcV0g3puAq4HNNfb9K/BU1eungD9Vve5dVZ8oIKWqnuEBUIcOwNVVr2OBb6piDbZ6GBBT9ToSWANcF2z1qIrtb4HXgEXB+G+qKrYCoM1Z+4KuHnWss9ovb+oR9G1YKLVfVfGpDbtIPQK1J6x6ZmvnXCkwDxjrcUzn5Zz7FDh41u6xwEtVr18CxtXYP885d9I5lw98R2V9PeWc+8FVzQrunDtC5cDlKwi+ejjnXEnV28iqzRFk9TCzTsAoYGaN3UFVhwsIlXqcj9ovD4RCGxYq7ReoDbvUegRqElbbzNZXeBTL5WrnnPsBKhsHoG3V/oCvm1Uu9zKAym9hQVePqi7wHGAf8IFzLhjr8Wfg74GKGvuCrQ5Q+QtkuZl9bZWzy0Nw1qMuQqEeQf0zCuY2LETaL1AbBpdQj4CaMb+GOs1sHWQCum5mFgO8DfzGOVdsVlu4lUVr2RcQ9XDOlQNpZtYSeNfM+l6geMDVw8xGA/ucc1+b2dBLOaSWfQHxswAGO+f2mFlb4AMz23aBsoFcj7oIlXrUJuDrFuxtWLC3X6A27CwXrEeg9oSFwszWe82sA0DVn/uq9gds3cwsksrGa45z7p2q3UFXj9Occ0XACmAEwVWPwcAdZlZA5a2sYWb2KsFVBwCcc3uq/twHvEtl13zQ1aOOQqEeQfkzCqU2LIjbL1AbdtpF6xGoSVgozGy9AJhQ9XoC8F6N/febWZSZpQDdga88iO8MVvl1cRaQ65z79xofBVs9Equ+QWJmzahcgHkbQVQP59zTzrlOzrlkKv/tf+yc+yVBVAcAM2thZrGnXwO3ApsJsnpcBrVfHgiFNiwU2i9QG1anejTWEwZ13YDbqXy6JQ/4B6/juUisc4EfgFNUZsKTgNbAR8C3VX8m1Cj/D1X12g6M9Dr+qphuoLLbdCOQU7XdHoT1uApYX1WPzcD/rNofVPWoEdtQ/vpkUVDVgcqnAzdUbVtO/z8OtnpcZt3VfjV+PYK+DQu19qsqPrVhF9g0Y76IiIiIBwL1dqSIiIhISFMSJiIiIuIBJWEiIiIiHlASJiIiIuIBJWEiIiIiHlASJp4ys5ZmNtXrOERE6krtl9SXkjDxWktAjZiIBCO1X1IvSsLEa/8CdDWzHDP7N6+DERGpA7VfUi+arFU8ZWbJVM6mfKFFakVEAo7aL6kv9YSJiIiIeEBJmIiIiIgHlISJ144AsV4HISJyGdR+Sb0oCRNPOed+Aj4zs80a2CoiwUTtl9SXBuaLiIiIeEA9YSIiIiIeUBImIiIi4gElYSIiIiIeUBImIiIi4gElYSIiIiIeUBImIiIi4gElYSIiIiIeUBImIiIi4oH/Hxx50FDcWxUUAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "tm = 0.05*(t.max()-t.min())\n", "tfit = np.linspace(t.min()-tm, t.max()+tm, 50) \n", "Nfit = N0*np.exp(-tfit/tau)\n", "fig = plt.figure(1, figsize=(10, 3.5))\n", "ax1 = fig.add_subplot(1,2,1)\n", "ax2 = fig.add_subplot(1,2,2)\n", "ax2.set_yscale(\"log\")\n", "ax2.set_ylim(0.01, 10.)\n", "for ax in [ax1, ax2]:\n", " ax.errorbar(t, N, yerr=dN, xerr=None, fmt='oC3', ecolor='k', ms=4)\n", " ax.plot(tfit, Nfit, '-', color=\"gray\", zorder=-1)\n", " ax.set_xlim(-10, 500)\n", " ax.set_xlabel('t')\n", " ax.set_ylabel('N')\n", " ax.text(0.95, 0.95,\"$\\\\tau = {0:0.1f}\\pm{1:0.1f}$\\n$N_0 = {2:0.2f}\\pm{3:0.2f}$\"\n", " .format(tau, dtau, N0, dN0), fontsize=12, \n", " ha='right', va='top', transform=ax.transAxes)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A similar procedure can be used to fit data to a power-law $y=Ax^p$ with $A$ and $p$ as the fitting parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Need for linfit.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Currently, there is no function available in numpy or scipy expressly designed to perform chi-squared least square fitting of a straight line to a single data set with weighting (error bars). There are two functions available in numpy and scipy that can be adapted to fit single straight lines to data with weighting: ``numpy.polyfit`` and ``scipy.linalg.lstsq``. Both of these have drawbacks if the desired task is to fit a straight line to a single data set. In addition to these two functions, ``scipy.stats.linregress`` fits a straight line but without any provision for weighting. \n", "\n", "* numpy.polyfit and scipy.linalg.lstsq are both slower than ``linfit``: ``numpy.polyfit`` is generally a few times slower; scipy.linalg.lstsq is also a few times slower but can be several hundreds of times slower depending on the weighting used and the size of the data set. This is because, in part, both ``numpy.polyfit`` and ``scipy.linalg.lstsq`` involve matrix inversion while linfit does not. More generally, both ``numpy.polyfit`` and ``scipy.linalg.lstsq`` have significant overhead associated with the more general fitting problems they are designed to address.\n", "\n", "* In its current configuration, ``numpy.polyfit`` does not allow absolute weighting of the data; only relative weighting is implemented. ``linfit`` allows either relative (the default) or absolute weighting (by setting the keyword argument relsigma=False)" ] }, { "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.7" } }, "nbformat": 4, "nbformat_minor": 1 }