{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Class 5: Error estimates\n", "\n", "In this class we will review methods to determine statistical errors by re-sampling data, Monte Carlo simulations or error propagation. At the end of these activities you should be able to:\n", "\n", "- generate errors through re-sampling the data using bootstrap or jack-knife procedures\n", "- propagate errors in different quantities in linear or non-linear combinations\n", "- use Fisher matrices to forecast parameter errors\n", "- model errors using Monte Carlo simulations\n", "\n", "Often we have no analytic model for the statistics we estimate from our data. However, we can still determine errors in these statistics using **approximate sampling procedures**.\n", "\n", "The basic idea is to build up many _statistical realizations_ of the data by random re-sampling, and use the **scatter across these realizations** to estimate the error.\n", "\n", "## Jack-knife errors\n", "\n", "The **jack-knife procedure** allows us to estimate the error in a statistic by re-sampling the data (without replacement):\n", "\n", "- Given a dataset with $N$ entries, we create $N$ separate datasets, deleting 1 entry in turn\n", "- We measure the statistic for each of these datasets, creating $N$ measurements $S_i$\n", "- The jack-knife error $= \\sqrt{N-1} \\times \\sigma(S_i)$\n", "\n", "The factor $\\sqrt{N-1}$ is required since the $S_i$ are correlated with each other, given that the datasets $D_i$ all significantly overlap.\n", "\n", "Let's apply this procedure to the problem of fitting the straight line from Class 3. Here are the 10 jack-knife samples and best fits:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd5xU1f2/nw9L35UuShMEJAgoZRFsiIqCKAoG+RpiwWhiTFFj1yiCopEk+rPHYNRgRSxYUGkqYBchgooKKqKIFFFBFmmz8/79cWbY2TK7M7MzO3fZ87xe97Uzc9tn7zPnnHvOuXOOScLj8Xg8nqBRK9sBeDwej8dTFr6A8ng8Hk8g8QWUx+PxeAKJL6A8Ho/HE0h8AeXxeDyeQOILKI/H4/EEEl9AJYCZdTAzmVntbMfiiY/3FHy8o+ATJEcZKaDMbKWZbTWzAjNba2aTzSwvE+cKOmbWw8xmmdkGMwvUj868pyLMbIyZLTKzn8zsGzP7RxASqHdUhJn9ysyWmdkmM1tvZg+aWaMAxOUdlYGZvVLZgi6TNagTJeUBvYDewFUZPFeQ2Qk8AZyT7UDi4D05GgJ/AVoA/YFBwKVZjagI78jxJnCYpMZAR6A2cEN2Q9qFdxSDmZ0G1KnscTLexCdpLTALJw4AMzvYzN4ys41mtsTMjoy3v5l1NrP5kbumDWY2NWbd7Wa2KnLXu8jMBsSsG29mT5rZI2a22cw+NLMuZnZV5O5rlZkNjtl+npndZGYLIsd7zsyaxYmpsZndb2ZrzGy1md1gZjlx/v9lku4HliZz3aoa70n3SHpd0g5Jq4FHgcOSuIQZxzvSKkkbYj4qBDoncOmqjJruKLo9MA64PMHLFpeMF1Bm1hYYCnweed8GeBF359MMd5f6tJntGecQE4DZQFOgLXBnzLr3cF+EZsBjwJNmVj9m/YnAw5F938d9cWoBbYDrgUklznUmcDbQCggBd8SJaXJkfWfc3dJg4Ldxtq0WeE+lOIKA3VR4R2Bmh5vZJmAzMBK4Ld622cA7AuBvwD3A2nK2SQxJaV+AlUAB7ksk4BWgSWTdFcDDJbafBYyJc6yHgHuBtgmc90egZ+T1eGBOzLoTIzHlRN7vEYktGtc8YGLM9t2AHUAO0CGybW1gL2A70CBm29HA3Api6+wud/qvt/eUPk+R7c4GvgFaeEeBddQmElcX7yg4joC+wOLIvruOk+q1zWQNaoSkPYAjga64tn2A9sCoSHV3o5ltBA4HWpnZAHMdjQVmFr17vRwwYIGZLTWzs6MnMLNLzeyTSHV4I9A45jwA62JebwU2SCqMeQ8Q25m5Kub1V7g21NjjReOvA6yJiX8S0DKhqxI8vKcYzGwEcBMwVMWbk7KJd1QCuWbYmcDjFW1bRdR4R2ZWC/gXcKGkUJlXKUky/pSSpPlmNhm4GRiBuygPS/pdnF2KPf0i16b7O3DVe+BlM3sNVy29HNeZvVRS2Mx+xMlNlXYxr/fBPeCwocTnq3B3FC3SJSEIeE9gZscB/wFOkPRhJeLLCN5RKWoDnVKOMAPUcEeNcDWoqWYGrjYG8I2ZjZL0erIBVtXvoG4DjjWznsAjwIlmNsTMcsysvpkdGWm7LYWZjYpZ9yOuyhjGVVlDwHdAbTO7FneBKsPpZtbNzBri2myfirkDAUDSGlwb8S1m1sjMaplZJzMbGCd+i7QT1428r29m9SoZZ6aoyZ6Oxj0YMVLSgkrGl0lqsqPTzGyfyOv2wI245rSgUVMdbQJa4/rJegHHRz7PB95NJcAqKaAkfYdrW71W0ipgOPBX3MVeBVxWTiwHAe+aWQHwPK76uALXjjsTWI6rnm6jeJU1FR7GdQiuBeoDF8TZ7kxcgfMx7kv0FO4Opyza46rX0Sr8VmBZJePMCDXc01hck8lLMc0uMyoZZ9qp4Y66AW+Z2RbcI+fLiNQ2gkRNdSTH2uiC+38B1knakUqAFunYqvGY2TzgEUn3ZTsWT3y8p+DjHQWf6uLID3Xk8Xg8nkCS0EMSZrYS9whlIRCS1DeTQXmSxzuqHnhPwcc7Cg4JNfFFhPUN0GO3nhJ4R9UD7yn4eEfBwTfxeTwejyeQJFqD+pKiRx4nSbq3jG3OBc4FyM3Nze/atWuaQ939WLRo0QZJ8YY8SQrvKHNUpSfvKDW8o+CTkqNEhpsA2kT+tgSWAEeUt31+fr48FQMsVPqGW/GOMkS2PHlHieMdBZ9UHCXUxCc3rAiS1gPPAP2SKgU9Gcc7qh54T8HHOwoOFRZQZpZrZntEX+NGsv0oI9Fs3w6nngqdO0P//rByZelttm2Dfv2gZ0/o3h3GjStad9ddbl8z2BCnf3PyZPjzn1OL74cf4NhjYb/93N8ffyx7u+OOgyZNYNiw4p+/+ir06QM9esCYMWkbZypwjlatgqOOgm7dnKPbby+9zS23xPdUFY6+/hoGD4b993dxRv+Pc85x360DD4RTToGCgtTiKIMq81RZR6eeCr16uaVDB/e3JPPmlf5+pzO+KIWF0Lt38XMlks5TpFqlpaB7Ouss2HffXTEeAg2SDSGRGtRewBtmtgRYALwoaWayJ0qI+++Hpk3h88/hoovgiitKb1OvnsvolyyBxYth5kx45x237rDD4OWXoX37jITHxIkwaBB89pn7O3Fi2dtddhk8/HDxz8JhGDMGHn8cPvoI2rfnz6UHZkyVYDmqXdsVQB9/7Nzcfbd7HWXVKpg9G/bZJ/3xJerozDOdp08+gQULoGVk7Mtbb3XfrQ8+cPHddVc6o6saT5V1NHWqS1uLF8PIkfDLX1Z9fFFuv93dRMSS2XRefdJS0D0B/POfu2J8u2jA2sRJtk0wkWVu48ZSnz5St27SpEmJN1IOHiy99ZZ7vXOn1Ly5FA7H337LFql3b+mdd4p/3r699N13Ze/z3/9KJ50kDRwode4sjR+feHxdukjffutef/utex+PuXOlE04oer9+vdSxY9H7117TPNioDFz/RJYqcyS56z17dtH7kSOlxYvje8q0o6VLpcMOK/844bB03nnSxIlp7d9IZsmqo+g1aNtWWr689PZz50oDBkjHH++u8e9/LxUWpje+Vauko4+WXnmleFqKEvP9yZaj/Px8afhw76ksT2PGSE8+uettKo4yMpr59R06cOSiRbB1Kxx0kCvdmzd31cVlZQxDd/HF7o529WpoFxlIt3ZtaNwYvv8eWpSoaBQWQn6+K9n/9CdX/UyGBQtcLaZhQxffCSdA374wYABs3lx6+5tvhmOOgXXroFVkCKq993bvE6VFCwiFYOFCd66nnqJVZADZbJBxR1FWroT33y9y9Nxz0KaNa0Yrj0w6Wr7cNcH+8pfw5Zduv4kTIScy+PJvfgMvveSaVW65Ba68svxYM0TWHEV5/XXYay/XXFoWCxa4u/n27V2z9rRprlk0XfH95S/wj3+U7TtIPPAANGvmPZXl6eqr4frrYdAg6qcw8npGCqhfrV9flAGtWuWaW5o3d1XSdJCT46qNGzfCySe7jKxHj8T3P/ZYFw+4TOqNN1zm93oSo8GbuSWZ7R9/3FWVt2+HwYMJJ7532sm4I3D9NyNHwm23QaNG8PPP8Le/uea9isiko1DIHef9910z3qmnun6vc85x6//7X3cTdP756b0eSZIVR7FMmQKjR8fft18/6NjRvR492jk65ZT0xPfCC67ZNT/f9aMEmTvugGeeca+9pyJuusndJO7YAeeeywTYO9nDp7+AMjuyX26uq900bAhHHukebICKS+w2bZzgtm1dJrJpU1EmVRZNmrgOxJkzkyugSmZa0fcV3Z3vtResWePu0NesKeq3SJRDDinKYGfPZsWNN27Lyi8oqsLRzp0uQZ12WlHb+BdfuBpLNNP95hv30MiCBe6LXDzGst+nw1Hbtq7jNppoR4xw7fvRAgrcTdCvfuXuDLNBthxFCYXcnfaiReXFWPb7dMT35pvw/POuJrttG/z0E5x+OjzySPx4skD+5s2uP+ztt72nkp6iLRn16sFvfkP+Qw/lxg8yDsm2CVa4wPD5jRu7RsdPPpHq1XPtoIlw112ujVSSpkyRRo0qvc369dKPP7rXP/8sHX64NH168W0q6oNq1Ur6/nu3/wEHSO+9l1h8l14q3XSTe33TTdJll8XftmQflCStW+f+btsmHX20ToRlykYfVKYdhcPSGWdIF15Y/rHK64PKpKNQSDrwQPddkqSzznL/VzgsffZZ0f9wySXSJZdkp38j245mzJCOOCL+OebOlerXl1ascH0agwdLTz2VvvhKniugfVAXd+okDRvm4vGeinuK9gWHw9KFF+puWKMkr28mEla9Nxo1krp2dZ2HAwcmLmzrVumUU6ROnaSDDpK++MJ9vnq1NHSoe71kidSrl8u0uneXrruuaP/bb5fatJFyclwGd845pc/x3/+6uI48MvkO+A0bXGdg587SoEEuA5Vc5hl7rsMPl1q0cF+MNm2kmTPd55de6q5Lly7SrbdmLVFl3NHrr7uv1gEHSD17uuXFF0sfq7wCKtOOZs928fXo4Tpzt293CfjQQ91n3btLv/61tGlTtgqo7DoaM0a6557456hM53si8ZU8V2zGV0Y6z1ZaOrh3b+m447yn6LliPR11VFFaOu00NYL/Kcnrm5H5oPr27auFCxem/bi7G2a2SFkaKdk7SpxsefKOEsc7Cj6pOPKDxXo8Ho8nkPgCyuPxeDyBxBdQaWL8+PGY2a5l/PjxSa33ZB7vqHpQngfvKBhUWVpKttMqkWV3HeF33Lhxwg3BL0Djxo0rtn7gwIEaOHBgmfuGQiH16NFDXdq313tXXaXwiBHZe0jCOypz36ijDh06aPr06QqFQpKU3VEKdkMqciTF91Smo3DYO8oA6U5LqTjy0pIgXgZW4bqdO3Vhfr7+DfoeJND3dev6AioDpOwoFNKgQYNUq1YtAcrLy9OwgQMVeuopn/mlmfI8lLc+1pGBjmrQQFP32Ufhjh29owyQjrS0N+jCevX0XrNmvoDKJGVlYIMGDVIoFIq/7pNPpLFjVbDXXhKoAPQwaDCocW6ugM/kC6i0kZKjSMKaPn268nJz9QvQRaCXQdsjNxM+80sfFXkob/30Z5/VcfXr6w7Qqoib7aB1ffp4R2mmMmnplXvv1V/r1tUboMKIp8/MUsrvapy0RJoXymL69OnKy8srtm9eXp6mT59ebF1z0J9AC2rVcpe3Vi193qmTzgDlxuxrTthq+QKqFJl2FLvuxSeflF58UQv69dMXkcQk0Iegf4AeOvtsn/mVQSYclbW+DmhE/fpaOWSICnJzJdDPoGmg00BNQBMmTPCO4lBVaal/gwb69Ne/dj+Aj6ShRaBrQN1A5rZLOr+rcQ9JjB07lh49etChQwemT5/O2LFjE9rv/fffZ8uWLcU+27JlC4sXL+aDBQsYWlDAc8Aa4C6gTjjMy0OGwKpVfHLbbTyTl0fs3rm5uQA/p+e/2r3IhKPYdfsCfwKmFhRw7OjRcMIJ9F6yhE9zcjgPaA8cAFyfl0fTk09O57+225AJR9H1hQUFDAceBNYBz2zbxt5z57IxP58z6tdnT+CXwKNAKC+PXmXNg+QBMpeWfi4o4GDg78BnwDtbt9JlyhRo1Iil55xDj4YNyQduAD4GcvPyIJX8LtkSLZElqHcVFVVNy6PkXYOBhjRooK+OPVY7GjaUQN+A/g7qUeKuMN558X1QpUino+j+L06bprcnTNCdderok5ha0mdm+uKkk6RZsxQqKIh73mx5qkmOXpo6VXr8ca0+7DBtjvbTgv4LGlW/vl58+ulyz+sdlSbdnprk5urtCRP05dCh+tZM0SbWl0Dn16unWQ89VO55U3FUo6RV1LxQHtGLvr+ZbgCtjAgK5+Wp8MwzdVnv3qrtmu3K/CKk66mWdC27s6NatWqpbSTRvNGihcKRpqHttWppJuh8UM+GDRNyJMlnfiVIl6OmoHPr1dObLVooXK+eS08tW+r5Nm002Ey14/RReUeJkQ5PuWYaDnq0dm1tql3bOWrYUPP33FOnmalxhvO7GtXEV1HzQlzWryfnrruYs3EjH0tcCTTo3Zvwww9ja9dS68EHuem99+javTsdOnRgypQpzJo1i5zo/EJATk4OzZs3p3379gwbNqzYOk8RKTvauZOcN99kdp8+fFKnDquAO7Zv59DcXOzMM2H6dHJ+/JFLe/RgeocO3DB1qneUIik7AnI2bGD2yJG83qAB64BJ27dzSIMG2HnnwWuvYd9+y/FffcW33bvTtoy05B0lTsqeNm4k5/HHmdOoEeuBZ4FT6tcnb/RoePZZbMMGDluzhiXdu9M00/ldsiVaIktQ7yqSuqPYskV67DE3yGJOjgRSnz66s1MnnXzIIaU2r6gzsqz1+BpUKZJytGaNG1h21CipcWPnqHZtLWrSRP/q2NHNnBszA2gqjiR/d16SpO/Mv/5auu02N2hppOVhVYMGerRdO2nBglKztJbnyTtKnKQ8rV3rZgMeMkSqU8elpVatNK11a1184IHSjh3FNq+q/K5GSauwTbaw0E1bfNZZ0h57uMvTrp105ZXS0qUpPxETD19AlaZcR6GQ9Pbb0tixUn6+on1Ju0auf/pp/e3KK9PqSPKZX0kS6tv47DNp4kSpX78iTz16SOPG6e7zzvOOqoAKPa1YId1yi5t9IXLjoE6d3BQ1b7+t8ddem/X8rkZKK9WG/eGH0uWXS23bukvSqJF09tlu+PhEh6ZPAV9AlU2so5mPPqrChx+WTjtNat7c+alVSzrsMOnGG6X33y91B55ufOZXmlLpaOdOl46uu67Yo8bq29fNy7VsWUbj8Y7Kppin559XaMkS6frrpd69ixz17OmmtPngg4ympVQcZWTK9yATbRvdLy+PYcuXw9ixbvr42rXhuOPczKwnnQQNGmQ71JpJOEzOkiX8paCA/j/8QI8zzoBwGFq0gOOPd8vgwdCsWbYjrdHk5OTQvFkzDq1bl2FvvQWXXALLl7vZWg87DG69FU4+Gdq3z3aoNZocMw6vU4cBhYUMu+QSNx09wKGHurzu5JOLZpYOIDWrgNqyhWlnnMFV8+dzDMA777C6dWva3HGHm957zz2zHWHNZNMmmDPHTRs9YwasXcs5wAJgPNDmt7/ld//+t5uG3ZNdwmF4+23evvRSJr/zDh2A0P/+x1cdO9LpX/+CESOKpvr2ZIdQCF57DaZN46eHHuKezZvZCcwGdpxwAsP+85/q4yjZKlciS6arvcl0oo4fO1aaNctNmxx53FgdOkjXXCN9+mlG46wIduMmvnL768Jh3f2HP+gy0FzQjmhTQ5Mm0qmnSg8+6DptA0K2PGXVUcz6HNDRoAV9+0p77+1c1a3rZk994AE3i3GW2V0dVUTUUT3QMND/evWSmjVzjho0kE4+WXroIemHH7Iap5Sao2orraLRjoe2bq3PR4xQuFUr9282biz97nfSa6/F7VdK90MQFbE7F1BSCUcFBdLzz7spp9u1U7T9e2mdOlo+apRC8+ZJO3cmdNya4imb6Ujbtin03HOa1qSJfogM2xVu2FAaOdI93bppU7nH9Y6qgE2bVPjII5rZqJEKor/LbNxYOv106emnXZorh+rgKJDSKrpwcUc7/uorTercWUsimd920BstWij0xBPS1q0VnreiUZbTTXUvoMrzFAqFdPx++2l8s2Za17u3wnXruq9bXp7Cw4fr5q5d1a6cHzaXR03xlOnMr9R13LRJeuop6de/VrhRIwm0ETfA8ej69XX8kUd6R1Wcjkqxbp30n/9Ixx+/K02tAd0DOql+fQ056qjdylEgpSUzLcLeubma2K2bwkcfrXDkLuJN0HmgZhX9PqPEOVMdFiRVglBAVeYuqpSnggJp5kwVnn++VjVooGgt6VMzPbHPPgrNmiVt25a2kQh2d0/pcBSP6HVsYqbRoGdzcrQ1OsBx8+b66thjdXL9+qrrHVV9Oip5nb76yv2O7Igj3BOsIO27rz4fMULHNGigWruxo8AVUIlMi9A4N1fHgR4FbYlkggWtWmnekUeqc4wscKOGT5gwocLzVibTTJUgFFBS+ROPxSPqaV8znQd6KSaDC9Wtq5k5OfojaN8yruX1118fHc3de8qgo7hs2KDFF1ygmTk52hZJP6tB/65TR2/deKO0c6d3VMXpqFR+9+GH7mcUsb/369FDuvbaXT+tqAmOAjfU0YwZM3j33XcJh8MAFBQU8O677zLjpZdg0SJa/O1vfLplCzOAIcBk4FDg1j/8gc2XXMJaN2ruLnJzcxMa7bgyw7dUZwoLC/n+++/56quveOGFFygsLCx/hx07YO5cVo4axZ1z57JC4h7gF4WFTM7J4d1x4/jHFVcwNBzmX8CXkd1ir2Xv3r2jo7nvwnuKT9KOymLNGvjXv2DQINhrL3recQddCgu5E5d+2gJ/CIV4JRyG2rW9oyRJ1VFsfpcPXFVQwN2vvkrOAQfA1VdDnTrw97+7R/g//BCuuw569QKzmuEo2RItkaUyNaiSdwXtQFeB1u+5p7s7r11bz+bkaDhurhhiSv5MjLIctDuKdC35+fmJX6/Vq6X77pN++ctdI2yEcnI0G/QXUJcSd28VXUvvKQOOyuLLL91IAYceWjRSQNeu0tVXa/6ttyrPTZrpHWXLUSikyb/5jW4HfRWpJe0EzQG9NGyY9M03Fey++zsKXAE1ffp0tc7N1dm4R5Cj1dsN3bpJkyYp9N13Fc7ImUrHX3Vpk03Xkp+fH/dL+sKzz0pvvCH99a9Sr167HKhtW+ncc6Vnn9VLU6fG/YInci29p8o5ipuRfPqpaxrq06fIW69ebvSApUt3beYdZcnRtm3SCy+4oblatJBwEzA+CzqT5PrNpd3fUXAKqB07pOnTVThqlLZF+jI+BV1Xt65OO/TQpBJOqu311eGplnQt+fn5xWqrLUCngx4D/Rx9wCEnx3XMTpxYahiURKburuhaek/JOYouxfoZwmFp8WI3PmH37tpVKPXvL/3jH9Lnn8f9v7yjKnL000/S1KnuN37RMT4bNZJ+/WuFpk7VCQMHVqqg2J0dZa2Aij71chDoDlBBZNI/tWihwj/9Sb/q2FEd2rcv88KlMtpxoqS1M7oCsl1AvXbLLbqhbl29AyqMZGxrzfT1oEHSE09IP/5YbvzlfcFTHTk8UWqCp7h357m5ev3mm6XLL9f3TZtKoFCkxeHF446TVq1K6P/yjjLnqH1urhZfcIE0bJgUmetqHWgSaAjo+quv3hV7ZQqK3d1RdgqoFStUOH68VkSe4w/VqaPCUaOk6dOlHTuq/AdkUarDD9fSteRHJh8rBL0FugZ0eIMGOubooxNOIN5T5jO/aE21tpkGgO6uU0frIhmeatdWePBgjWvdWvnt2lXJXXCi1ERH+5jpz6D5OTkKRWuy7dtLf/mLQnPn6sDu3VO6mcsU1cFR1RVQP/zg5hs5/HBFmyFeBZ0Nap2bm/H2zyCS1QKqWTPpkUcUWrOmSqv51ZGsZX59+kizZqnwd7/ThsicZKG6dRU+6STpoYcq7I+tSWTNUffu0k03Kdy376587ad27VT4179KixZJ4XBW+nuCSPAKqO3bpWeeccOjREcS2H9/fXzmmerasGHiHb+7Kdlu4pOyd/dWncha5hcplLbVraspoFGg3BhH2XgSK6hkzVGkUPqmdWtdQdETrbHpyHtypOIood9BmdlxZrbMzD43sysr3GHLFvjjH92IuSefDK+/Dn/4AyxcCEuX8lTnzizburXELgF8Br8akbSjCOPHjy/2hRg/fnwGo/Qk5alJE3juOept2sSvJJ6QKIhxVG1+y1LNSMpRu3bw9de0Wb2aiRLLykhH3lPqVDjdhpnlAHcDxwLfAO+Z2fOSPo6706efwldfuaH3zzgDjj3WzbcUIfoDs4KCgl2fJfoDM09pUnLkqXKS9tShg5ubLA4+HaWfpB21bOkKqXLwnlInkRpUP+BzSSsk7QAeB4aXu0eHDrBuHTz2GAwdWqxwAhg6dCj9+/enVi13+ry8PPr378/QoUNT+R88qTjyZIO0evLpKCOkPS15T6mTSAHVBlgV8/6byGfxad4c9tgj7uqcnBxmzZpFt27d6NChA1OmTGHWrFnk+AnpUiV5R55skFZPPh1lhLSnJe8pddI2o66ZnQucC7DPPvtUuH106vXmzZszbNiwdIXhKYdkHXmqHp+Ogk8q6ch7So1EalCrgdhG1raRz4oh6V5JfSX13TOBqdPHjx/P/PnzmT9/PmbmO+crR0YcedJOhZ58Oso6aXcE3lPKVPSYH66WtQLYF6gLLAG6l7dPtqdBri6QpkdjvaPMki1P3lHieEfBJxVHFTbxSQqZ2Z+BWUAO8ICkpZUuGT1pwzuqHnhPwcc7ChbmCrY0H9RsM7As7QeuHC2ADdkOogS/kBT/aZIM4h0lRVY8BdQRBNOTd1Sc3cJR2h6SKMEySX0zdOyUMLOFQYwpi6f3jhIki54C5wiC6ck7Ks7u4ihwM+p6PB6PxwO+gPJ4PB5PQMlUAXVvho5bGXxMwTl3PIIYE2QvLn89Esc7Kk4Q40o6pow8JOHxeDweT2XxTXwej8fjCSS+gPJ4PB5PIElrAZXqnESZxsxWmtmHZrY4W4+jmtkDZrbezD6K+ayZmc0xs88if5tWQRzeUflxeE9x8I5KxRE4RxAMT+lylLYCKmYelaFAN2C0mXVL1/HTwFGSemXxtwGTgeNKfHYl8Iqk/YBXIu8zhneUEJPxnsrDOyLwjiD7niaTBkfprEH5OYnKQdJrwA8lPh4OPBh5/SAwIsNheEcV4D0FH+8o+KTLUToLqCDPSSRgtpktigyVHxT2krQm8notsFeGz+cdpYb35PCOigiqIwiup6QdZWqoo6BxuKTVZtYSmGNmn0ZK+MAgSWZWk5/5D7wjqPGevKPqQeA9JeoonTWohOYkygaSVkf+rgeewVXPg8A6M2sFEPm7PsPn845Sw3vCOypBIB1BoD0l7SidBdR7wH5mtq+Z1QV+BTyfxuOnhJnlmtke0dfAYOCj8veqMp4HxkRejwGey/D5vKPUqPGevKNSBM4RBN5T8o6SnUCqvAU4HlgOfAFcnc5jVyKmjrhJx5YAS+LS68QAACAASURBVLMVFzAFWAPsxLVXnwM0xz3N8hnwMtCsCuLwjrwn72g3dBQkT+ly5Ic68ng8Hk8g8SNJeDwejyeQ+ALK4/F4PIHEF1Aej8fjCSS+gPJ4PB5PIPEFlMfj8XgCiS+gPB6PxxNIfAHl8Xg8nkDiCyiPx+PxBBJfQHk8Ho8nkPgCyuPxeDyBxBdQHo/H4wkkvoDyeDweTyDxBVSCmFkHM5OZ1ZRJHqsd3lHw8Y6CT5AcZayAMrOVZrbVzArMbK2ZTTazvEydL8iY2VlmVhi5FtHlyADE5R3FYGYdzewFM9tsZhvM7B8BiMk7imBm/y6Rhrab2eYAxOUdRTDHDWa22sw2mdk8M+ue6vEyXYM6UVIe0AvoDVyV4fMFmbcl5cUs87IdUATvCIhMOjcHeBXYGzdD6iNZDaoI7wiQdF5sGsLNOfRktuOK4B05RgFnAwOAZsDbwMOpHqxKmvgkrQVm4eQBYGYHm9lbZrbRzJaUV6Mws85mNj9SIm8ws6kx6243s1Vm9pOZLTKzATHrxpvZk2b2SOSu+EMz62JmV5nZ+sh+g2O2n2dmN5nZgsjxnjOzZnFiamxm95vZmsjdwg1mllPJS5U1vCPOAr6V9P8kbZG0TdIHiV6/qsA7KrZfLjASeLCibasS74h9gTckrZBUiLvJ65bg5StFlRRQZtYWGAp8HnnfBngRuAFXyl4KPG1me8Y5xARgNtAUd2d7Z8y693BfhmbAY8CTZlY/Zv2JuBK8KfA+7stTC2gDXA9MKnGuM3F3AK2AEHBHnJgmR9Z3xt0xDQZ+G2dbgN6RL9xyMxtrAWjfjcU74mBgpZnNiHiaZ2YHxNk2K3hHxRgJfAe8lsC2VYZ3xONAp0jhWAc3tfvMONtWTAan/F0JFACbAeGm+m0SWXcF8HCJ7WcBY+Ic6yHgXqBtAuf9EegZeT0emBOz7sRITDmR93tEYovGNQ+YGLN9N2AHkAN0iGxbG9gL2A40iNl2NDC3nGmY98V9WQ4APgauysZUzN5RXEezcdNTDwXqApcBK4C63lEwHJWI7xVgfDbdeEdlxlQXuD2yfwj4Etg31Wub6RrUCEl7AEcCXYEWkc/bA6MiVd6NZrYROBxoZWYDrKgTdGlk+8sBAxaY2VIzOzt6AjO71Mw+iVSJNwKNY84DsC7m9VZgg1zVM/oeILZDc1XM66+AOiWOF42/DrAmJv5JQMuyLoJcdfdLSWFJH+LuZk4pa9ss4B0VnecNSTMk7QBuBpoD+8fZvirxjmIws30i1+Kh8rarYrwjx7XAQUA7oD5wHfCqmTWMs325VEkzk6T5ZjYZl+hH4C7Mw5J+F2eXYk/AyLXr/g7AzA4HXjaz13BV08uBQcBSSWEz+xEnOFXaxbzeB3dXvaHE56twdxUtJIVSOIeoXIxpxzviA+CwSsSUcbyjXZwBvClpRSXiywjeEb2AqZK+ibyfbGa34WpoC5MNsCp/B3UbcKyZ9cR1nJ1oZkPMLMfM6pvZkZH221KY2aiYdT/iMvgwrtoawrVF1zaza4FGlYzzdDPrFinxrweeirkLAUDSGlyT0C1m1sjMaplZJzMbGCf+oWa2V+R1V2As8Fwl48wENdZR5P892MyOMdcB/BdcYv2kkrGmm5rsKMqZuH6RoFKTHb2HqzHuFdn2DFwN7PNUAqyyAkrSd7gq+bWSVgHDgb/iLvgqXJt/vHgOAt41swLgeeDCyN3TLFwH3HJcFXUbxautqfAw7su/FldFvSDOdmfi2ls/xn2RnsLd5ZTFIOADM9sCvARMA/5WyTjTTk12JGkZcDrw78i2w4GTIs19gaEmOwIws0NwDw8E5fHyUtRwR38HlgCLgY3ARcBISRtTCdAiHVse3KOXwCOS7st2LJ6y8Y6Cj3cUfKqLIz/Ukcfj8XgCSUIPSZjZStwjlIVASFLfTAblSR7vqHrgPQUf7yg4JNTEFxHWV9KGjEfkSQnvqHrgPQUf7yg4+CY+j8fj8QSSRGtQX1L0yOMkSfeWsc25wLkAubm5+V27dk1zqLsfixYt2iAp3pAnSeEdZY6q9OQdpYZ3FHxScpTIcBNAm8jflrhHCI8ob/v8/Hx5KgZYqPQNt+IdZYhsefKOEsc7Cj6pOEqoiU/S6sjf9cAzQL+kSsFE2b4dTj0VOneG/v1h5crS26xaBUcdBd26QffucPvtRetOPRV69XJLhw7ub0nmzYNhwzIXX5TCQujdu/i5Bgwoiq91a+ZAp9QCKY13lER827ZBv37Qs6eLb9y4onWvvgp9+kCPHjBmDIRSGSgkPlXiqbKOliyBQw6BAw6AE0+En34qvX8209FZZ8G++xZ9jxYvTi2OOFSbtARw553Qtatbd/nlpffPpicJrr4aunSB/ffn6gqGsCqLCgsoM8s1sz2ir3Ej2X6U7IkS4v77oWlT+PxzuOgiuOKK0tvUrg233AIffwzvvAN33+1eA0yd6r6sixfDyJHwy19WfXxRbr8d9i8xjNvrrxfFd8ghPOt+yFZpvKMk46tXzxVES5a4OGbOdHGGw65Qevxx+OgjaN8eHkzfbA5V5qmyjn77W5g4ET78EE4+Gf75z6qPL0pZ6QhcTNHvUVk3OSlSrdLS3Lnw3HPue7x0KVx6adXHF6UsT5MnuwL200/hk0/4L/yQdAwVVbFwI3EviSxLgasr2mdu48ZSnz5St27SpEmJ1wEHD5beesu93rlTat5cCofL3+ekk6TZs4t/Fg5LbdtKy5eX3n7uXGnAAOn446UuXaTf/14qLExvfKtWSUcfLb3yinTCCaXXb9okNWmipvA/padJwjtKNb4tW6TevaV33pHWr5c6dixa99pr0tChaWs+StZT1hw1alS0/ddfS/vvX3r7bKajMWOkJ58stmm2HOXn50vDh2fH06hR0pw55W+fTU8HHSR99tmut6k4qrTQspajevZ0Ef38s9S9u7Rhg3v/f/8n9exZennwQbe+e3f3z0bp2FH67rv4F/DLL6V27VyGH8v8+VK8tuG5c6V69aQvvpBCIemYY4q+7OmKb+RIaeFCd66yCqgHH5RGjkxru3myS413FAq5/XJzpcsvd5+Fw9I++0jvvefeX3CB1KNH1jxlzdEhh0jPPONe33KLlJdXep9spqMxY1xme8AB0l/+Im3bljVH+fn50vffu7iq2lPPntK110r9+klHHCEtWFB6n2x6atZMuuEGl86PO0494MNkr29GRjP/1fr1rn0fXBXvs8+geXPXvJMuCgpcE9Ftt0GjEmMmTpkCo0fH37dfP+jY0b0ePRreeANOOSU98b3wArRsCfn5rv23LKZMcc0oTz9d+fOlSI12BJCT45qHNm50zVgffeT6nR5/3DVnbN8Ogwe77bJE1hw98ABccAFMmAAnnQR165a9b7bS0U03wd57w44dcO658Pe/V/58leGOO+CZZ9zrqvQUCsEPP7imv/feg//7P1ixAqzEAOfZ8rR9O9SvDwsXwrRpTJ45s0Oyh09/AWV2ZL/cXNdu2bAhHHmk65QG1+G2bFnpfS6+GM48E9q0cYLbtnUXf9MmJ7okO3c6WaedVroPIxSCadNg0aLyYiz7fTrie/NNeP55eOkl93//9BOcfjo88ohbv2EDLFhQ9IXOBjXdUSxNmrhO6JkzXQF1yCGurxBg9mxYvty18Vc12XTUtav738H9/y++GC/Gst9nOh21ioxTWq8e/OY3cPPNZcdXBeRv3gwvvwxvv131ntq2de/NXCFUq5bLX/Ys8SR3tjxF4wM4+WS6QIO4FzIeyVa5Klxg+PzGjV0V75NPXPVy7tz41dZY7rrLtZFK0pQpro21JOGwdMYZ0oUXln2MGTNcdTcec+dK9etLK1a4ttjBg6WnnkpffCXPVbKJ7557pDPPlKTsNfHVdEfr10s//uhe//yzdPjh0vTp7v26de7vtm272tWz4imbjqLXoLDQbXP//aW3yWY6+vbbov/hwgulK67IWlq6uFMnadgwF09Ve7rnHmnsWPd62TLXp1uyjyibnq64oui7M3euPoQtSvL6ZiJh1XujUSOpa1fXeThwYOLCtm6VTjlF6tTJdbB98YX7fPVqaehQ9/r1113YBxxQ1Gb64otFxxgzxokr7yKm2mmYSHwlz1WygBo40GXQymoBVbMdLVki9erl4uveXbruuqL9L73UXZcuXaRbb5WUJU/ZdHTbbdJ++7nliivK7hjPZjo66iipRw/n7rTTpM2bs5aWDu7dWzruuOx42r7d/f/du7sHfV55pexrly1PP/7oztujh3TwwTrYTbSY1PXNyHQbffv21cKFSU+eWOMws0XK0kCU3lHiZMuTd5Q43lHwScWRH4vP4/F4PIHEF1Aej8fjCSS+gKoixo8fj5ntWsaPH5/tkDwl8I6CTylH48a5J+g8gSJtaSnZTqtElt11AMVx48YJN8KxAI0bNy6hdVEGDhyogQMH7npPFn+o6x0l5kjKnqea6CiR9QOPOEK/7dNHuuwyqX17yW3nHVUhVZXfeWlJUlYGJkmhUEg9evRQhw4dNH36dIVCobjrX5k0SYUTJ/oCKkOkw1Hsep/5pZ94jqRyPH30kQr/+ld9WbeuBCrMyVF46FDpwQe9oyomlbTkC6gME09KKBTSoEGDVKtWLQHKy8vToEGDitbv3Klz+/XT9aAP3Bi/yuZdn7yj0o5KrN+/YUPdtd9+Cg8c6DO/NFNe5lbSQ8+GDfVAx44K9+ghgUKgOaBzQPvk5u5y6B1VHcmmpRa5ufprz56+gEqERKqmZVGelOnTpysvL6/YcRvl5urNiROliy7SlpYtdyWuuaALIhkg8Jl8AZU2knWUl5en6ZEf6E5//nn1a9BA14AWxdxEbGrf3md+ZZCJdCRJ06dP1/4NG+pS0HsxHjZ066YPzjtP+7p0U8qhd1Q2qXoqjwrT0vTp6tiwoc4BPQvaUnRDnnR+l5Gx+ILM+PHjmRcZM2pevLHyymDGjBm8++67hMNhAAoKCnj33XeZMWMG77//Plu2bKEucDRwMjB8yxb2uvJKqFuXbzp04O/r1/M8sCFyPNu6FaBhuv4vT2KOYtlaUMCG556DefM49L//5d2tWwkDbwOXAs8BY377Wxg7tqr/lcCT7nT06qOPcuzGjfS8+WY+/vlnABYAFwNPAeeOHo0kVrp0s4stW7awOM3zQe1OpOqpPMpKS1sKCvh25kxYsoTe99zDFxGHK4H7gRfcZknndzWugCosLOT777+noKCAF154gaFDh5KTwICgZUrZsoWl777L8Tt38oucHIaEQjQGfgJm167NvhdfTP4117B8/nyeGD2agoKCXfvm5uZSUFDwc5r/vRpNPEeLFy+md+/e5ObmsqOggEHACGCEGS3vuw/q1iV0wAFcsGULT2zfzrrIvnl5efRK41xDuxPpSEfNgZHAqQUFHHXWWSDRuEMHxtety8M7drAisk+sh0i62XW83Nxc76gcUvVUHtG0tL2ggIHAScBJZrS/+24AGnbpwnV16/L0jh18GNknLy8PUsnvkq1yJbIEtdpbUfNCecRWa5uDzgK9lJOjUJ06EuiHOnX0H9BQULOYtvHyzovvg0or8ZoeZkydqsJHH9Xcli31U7TpDvRqy5YqfOwxadOmcr8b2fIUVEeVSUczpkzR7+vV00ugnREXy8y0bPRo6eOPyz22d5QclfEUlw0bVDh5sua1bKlNEX9bQG+2aKHCSZOkb79Na35Xo6RV1HZaHqEvv9SdXbro1UhfkkBr6tdX4QUXSPPnK7R9e5U81ZKuJaiOKkNswtgLdH69enq3eXOFIzcR4ZYt9WTTphrTsqVenDbNP8WXIkmno82bpccek4YPVzjyBN4K0N9AhzRooEFHH13qQYnyHqLwjhKjMvldMT79VPrnP92YfrVqubTUqpWebNpUZ7dsqReffNI/xZcOrr/+eplZMWFmpgkTJpS9w8cfSzfeKPXtq2hn7fJ69XRH48aaf+utCu3cWWzz8h6dLauz0hdQaeazz1Q4caL+17ChCiO+wp06uQFg33hDCoWSdiT5zK8kCaWjn3+Wnn7aTYrXoIFLP61bSxddpNCbb6pH9+5xb+ak+GnJO0qcpPO7KDt3uglFL73UDTAbfVilZ083evqCBVJhYZXkdzVKWoV3FOGwu/hXXeVGJ46K6d9fmjhRd/z5z3GfiEnlaRlfQFWScNjN5HnNNW5E54ivhaBrQN1B4669dtfmqT7R5DO/4sRLRy9Mmya98IKbHmKPPZyPPfeU/vhHl+FFRtGu7A91y8I7Kk1SNahNm6QnnpBOP93NhAtSnTpueo677pJWriy2eVU5qlHSymobPfaooxSaM0c6/3w3nTJIOTnSoEFOzDffZCweX0ClwM6d0quvuunY99nH+apVSzrySDdNRImElA585lec2HSUAxpWv75ebN1a4aZNnY8mTaRzzpFmz3a+qgDvqDQV9kF9+aV0xx3Ssce6wgik5s3dfHVPPlk0tXyaSMVRjXqKLycnh1mzZtHvwAPpvWED1x54IO3efx879lg3NfGQIXDDDTBsGDRrlu1wPVF+/hnmzHGzEE+f7qa5rl/fTcl+3XXOV4sW2Y6yxpBjxqxrruGJ99/n2I0babFtG/rpJ2zECDdL6+DB8aeJ91QZ0fyuV69eFBQUcOfttzO0RQtyxo1zs+B+GHnGrmtXuOgiOPFEN6N0JZ/ySyc1p4DatAleeolPb7iB+R9/TB6w8eWX+eDAA+l5772ucMrNzXaUnig//AAvvADPPguzZrlCqkkTVxidfLL3VdVI8N57MHUqTJ1KzurVjACmA1OB3hdcwDU33pjlID0lmTh2LB0/+ogTgb7Dh5MDrgAaMABuucUVSvvtl+UoyyHZKlciS2CqvWvXSpMmuRkvo1XYvfeWzjtPmjXLzUiZRajBTXxltmGvWiXdeaebaj0nx/lq00b605+kOXOkHTuyEmu2PGXd0bXX6kDQjaAvov2xdepIJ50kPfqo9NNPWY0vlhrrqKx0tHq1y/eGDXPTvYPUqJF06qnSI49I33+flVhTcVQtpZXbQbdihXTLLdLhh0tm7l+MPsn15pvlTneciWFBymN3LqASGu34iCN0Zt++pZ6UVNeu7kGVyNNCqRw7neyumV/c6/jJJ9L48bseFAqBNGSI9MAD0g8/ZDSmVNldHUmJje5+Tn6+dP31xdPRvvu6vtqXXy7zZrw6pKNASqvowhV7xv755xV6/33puuukXr2K5PTs6RLZkiXuaa8EKe/RyXRT3Quo8jzF/S1LYaH09tsqvPRSrYj8JkagcL9+0k03ucwxAWqCp6p09PJ//qPCG2906QYkM4UHDtR1rVurT7t2cR8HDwq7qyMpTlratk2aOVOFf/iDvo3+zs9M4f79pb/9Tfrww4TyvaCno0BKq+iHesccfbQONdPfQZ9HaklhM+mww6Sbb5a++CLt580E1b2ASnTk8Ka5ubqiVy8V/v73UqtWEminmWaDzgPt17BhUr9wrymeMu3o1MMO00Vmeid6UwcKH3ywdNttCn39dfpHIcggu6Oj6Lqohxagc+vV0/w991Q4L08C/Vyrlp4B/Qa0726YjjImLe2jHW/dKs2erS+HDtW3kUJpO2gGbsSA2Q89VKmLl5FhQSogCAVUpkZ33ys3V6eAHgH9GMn8dtavL51yiv53ySVqk5ub2O8zkjhvpsh25pdORyMHDFDhnXdqQ4nfjV0G6tawYbERqdMyCkEVsTs52vV9Doc191//0jV16+p1ikawWW2mlccdp3fHjVOL3TwdZUxaqqVzbMJoCDoZ9Fjt2tqemyuBttepoydBo0GNSeLX0Umct6oSZBAKqHR4ii4dGjbU4vPP17IuXbQ1kpjWg+4DnQi6KfKj2ZR/4R7nvLurp3Q5ago6GzQ7JpNbv+eeGgvaL46DyjjKBtXdUfQa1wYdV7++vhg+XOrcWdGbiEWgcaA+IANNmDChRqSjWmSIIUOG8PHHH7Ny5UpGjx7NkCFDKCwsrHC/T958k18WFDAN+A6YBgwOhfikSxd47jlenjKF3+TlMQXYFNknHSMalzcS9u5KYWFhyp6i16sDcBEwH/j855/peeedtN24kfvq1OEIoBXwW2BuXh49DjoIKBoNOZZEHdY0Tyk7+uknQpMnM6WggLW4KQ86ABOBSX/+M+8+8AC35uXxWcwusQ4q46imUdl0VLuggFOBR4D1wIxt29jnhRdgv/348A9/oGvDhuQD1wH/A3Ijo7vXiHSUbImWyNK5c+fkSufVq6W775aOOUaFkcEIV4HuAB0Fapybu2vfTFVNq8sdRbqW/Pz81P7ncFhavFjLRo/WBxFXAi0G/a1OHc2//XaFdu5MasbNVEeV3909Je1oyxY3XM0vfynVqyeBvjbTPyJ33rH7Jjsrqu+DSnM6+uIL6bbb9N2BB2pHJA2tAz0AGl2/vl564glJ5XuoCekoI9Jat25dcdVz+XLp73+XDj54VyanX/xChZdfrj/27aucyP5lXfRMdO5VlzbZdC35+fmJNxGEQm4stYsuco+u4h5K+aBxY11spo5xMrhkR3dPhJrkKSFH27ZJzz0njR4tRZrBtffe0vnnK/Taazrm6KPLLYQy4SgbBNqR5NLQW29JV15ZbNzIcLdueqx9ex1mplop5He7ezrKiLQya1C5uZp/221uNNwYQe+B/gq6849/LHbxyrvomXp+vzo81ZKupcI7v61b9ejo0bov0o8k0M6cHOn446X//Edat67c65WJAUGj1BRP8Rw1yc3VO9ddJ511lrZGakrfgf4N+u+ZZ7rMMIFrlUlHVU3QHOXl5bla0LRp0m9+o4KGDSXQDtDLoBlDhkiffy6pcvnd7p6OMiZt0KBBqm2mw0F31qmjb6O/aK5VSxo4ULrtNv1f//5xn8Gvyufzo1SHH66la4k+yBJ7F9UmN1c3du+u8MiRu+7GN+fkaE7LltLUqWWOHOA9VY2j2mYaCLqvdm1tjI6K0qiRNGaMLjvgAB09YEDc+LPhqKoJgqNatWqpDeiCevX0TvPmCkduHtSkiTR6tK7bf3+dcNhhZcbv01FVFlCdO6vwnHO0ITJcTah2bYWPP1667z5p/XpJlbuz213IdgElSaGvv9b1rVppfv36KowOL7T33tLvf6/Qiy+qV7duVV6TDRrZzPz01lsqPP98ra9dWwLtrFdPhaeeKj37rLR1a9ZaG4JGVh0tWqTCsWO1NHoTTmQesosuciPv79jh8zul5igz0kDb6tbVFNAoUB6lf8FenTpgM0VWC6g2bdw8V5EEtRz0d9B/zj5bKiz0jmLIWuYXGWljZ06Ono6kpYYxack7KiJrjqKjOIBeB10O6krxeci8J0dwCqjOnV3nbRyq248AM0VWCyiQ8vOlG26QPvqo1LAo3lERWcv8GjWSHnoo7rw83lERWXPUpIk0efKulqGy8J4cqThK6HdQZnacmS0zs8/N7MoKd2jcGOrVi7u62jyDX41I2tEBB8DChXD11dC9O5gVW+0dZYakPO23H5xxBjRqVOZq7ygzJOWoUycYMwb23DPuJt5T6lRYQJlZDnA3MBToBow2s26VOan/EWB6SclRBRPKeUfpJ91pyTtKPz6/CxaJ1KD6AZ9LWiFpB/A4MLwyJx06dCj9+/enVi13+ry8PPr378/QoUMrc9iajHdUPUirJ+8oI/i0FCASKaDaAKti3n8T+awYZnaumS00s4XfffdduQeMTkXcrVs3OnTowJQpU5g1axY5AZpquJrhHVUPKvTkHWWdtDoC76kypG0sPkn3Suorqe+e5bTHRpkwYQIfffQRK1eu5MQTT2TChAnpCsUTB+8o+HhHwSdZR+A9pUrtBLZZDbSLed828lmlGD9+POPHj6/sYTwO76h6kHZP3lHa8WkpQCRSg3oP2M/M9jWzusCvgOczG5YnSbyj6oH3FHy8owBh7vH0CjYyOx64DcgBHpB0YwXbbwaWpSXC9NEC2JDtIErwC0l7pONA3lFGyYqngDqCYHryjoqzWzhKqIBKFjNbKKlv2g9cCXxMwTl3PIIYE2QvLn89Esc7Kk4Q40olpoxNWOjxeDweT2XwBZTH4/F4AkmmCqh7M3TcyuBjCs654xHEmCB7cfnrkTjeUXGCGFfSMWWkD8rj8Xg8nsrim/g8Ho/HE0h8AeXxeDyeQJLWAirpKR+qCDNbaWYfmtliM1uYpRgeMLP1ZvZRzGfNzGyOmX0W+du0CuLwjsqPw3uKg3dUKo7AOYJgeEqXo7QVUJkYpj7NHCWpVxZ/GzAZOK7EZ1cCr0jaD3gl8j5jeEcJMRnvqTy8IwLvCLLvaTJpcJTOGlTah6nfnZD0GvBDiY+HAw9GXj8IjMhwGN5RBXhPwcc7Cj7pcpTOAiqhKR+yhIDZZrbIzM7NdjAx7CVpTeT1WmCvDJ/PO0oN78nhHRURVEcQXE9JO0pkNPPdgcMlrTazlsAcM/s0UsIHBkkys5r8zH/gHUGN9+QdVQ8C7ylRR+msQWVkmPp0IGl15O964Blc9TwIrDOzVgCRv+szfD7vKDW8J7yjEgTSEQTaU9KO0llABXKYejPLNbM9oq+BwcBH5e9VZTwPjIm8HgM8l+HzeUepUeM9eUelCJwjCLyn5B1JStsCHA8sB74Ark7nsSsRU0dgSWRZmq24gCnAGmAnrr36HKA57mmWz4CXgWZVEId35D15R7uhoyB5SpcjP9SRx+PxeAKJH0nC4/F4PIHEF1Aej8fjCSS+gPJ4PB5PIPEFlMfj8XgCiS+gPB6PxxNIfAHl8Xg8nkDiCyiPx+PxBBJfQHk8Ho8nkPgCyuPxeDyBxBdQHo/H4wkkvoDyeDweTyDxBZTH4/F4AokvoBLAzDqYmcyspkzwWC3wXjye3ZuMFFBmttLMtppZgZmtNbPJZpaXiXMFHTOrZ2a3mtm3Zvajmf3LzOpkKRbvJYKZKefZVQAAFexJREFU9TCzWWa2oayZPc2smZk9Y2ZbzOwrM/t1NuL0eGoymaxBnSgpD+gF9AauyuC5gsyVQF+gB9AF6ANck8V4vBfHTuAJ3Dw1ZXE3sAPYCzgNuMfMuldRbB6Phypo4pO0FpiFyxABMLODzewtM9toZkvM7Mh4+5tZZzObb2abIne7U2PW3W5mq8zsJzNbZGYDYtaNN7MnzewRM9tsZh+aWRczu8rM1kf2Gxyz/Twzu8nMFkSO95yZNYsTU2Mzu9/M1pjZajO7wcxy4vwLJwJ3SPpB0nfAHcDZCV6+jFHTvUhaJul+3KRuJY+TC4wExkoqkPQGbjbQM+JfUY/Hk24yXkCZWVtgKPB55H0b4EXgBqAZcCnwtJntGecQE4DZQFOgLXBnzLr3cBlsM+Ax4Ekzqx+z/kTg4ci+7+My5FpAG+B6YFKJc52JKzxaASFcYVIWkyPrO+NqIYOB38bZFsBKvG5rZo3L2T7jeC/l0gUISVoe89kSwNegPJ6qJEPT/a4ECoDNgHDT/DaJrLsCeLjE9rOAMXGO9RBwL9A2gfP+CPSMvB4PzIlZd2IkppzI+z0isUXjmgdMjNm+G66JJwfoENm2Nq7JZzvQIGbb0cDcODHdALwJ7AnsDbwbOVarLEzD7L2Ujq2zSwbFPhsArC3x2e+AeVXtzC9+qclLJmtQIyTtARwJdAVaRD5vD4yKNCNtNLONwOFAKzMbEOnALzCzaNPL5bhaxwIzW2pmu5rHzOxSM/sk0sy0EWgccx6AdTGvtwIbJBXGvAeIfUhgVczrr4A6JY4Xjb8OsCYm/klAyzjX4UZcLWEx8BbwLK7/Y12c7TON91IxBUCjEp81whXsHo+nisj447mS5pvZZOBmYAQus3lY0u/i7FLsqTK5vpLfAZjZ4cDLZvYarrnncmAQsFRS2Mx+pHhzWrK0i3m9D64g2VDi81W4O/UWkkIVHVDSVuDPkQUzOxdYJClciTgrTU33UgHLgdpmtp+kzyKf9aSM/iqPx5M5qup3ULcBx5pZT+AR4EQzG2JmOWZW38yOjPSJlMLMRsWs+xHXpBPGNQWFgO9wmcm1lL7rTZbTzaybmTXE9YU8FXNnD4CkNbi+l1vMrJGZ1TKzTmY2ME78bcystTkOBsYC4yoZZ7qoyV4s0i9WN/K+vpnVixxrCzANuN7Mcs3sMGA4rt/M4/FUEVVSQMk9vfYQcK2kVbjE/ldcJrYKuKycWA4C3jWzAtyTVBdKWoHrH5mJu9v9CthG8aagVHgY19G+FqgPXBBnuzNxGdvHuMz5KVzNoSw64Zr2tgAPAldKml3JONNCDffSHtecGK0VbQWWxaz/I9AAWA9MAf4gydegPJ4qxKRSv1GskZjZPOARSfdlOxZPEd6Lx1Nz8UMdeTwejyeQJPSQhJmtxD3BVIj7fUjfTAblSR7vyOPx7G4k1MQXyfz6StqQ8Yg8KeEdeTye3Q3fxOfxeDyeQJJoDepLih4lniTp3jK2ORc4FyA3Nze/a9euaQ5192PRokUbJMUbSigpvKPMkU5PHo8ncRItoNpIWm1mLYE5wPmSXou3fd++fbVw4cI0hrl7YmaL0tVX5B1ljnR68ng8iZNQE5+k1ZG/64FngH6ZDMqTPN6Rx+PZ3aiwgIr8kn6P6GvcCNEfZSSa7dvh1FOhc2fo3x9Wriy9zapVcNRR0K0bdO8Ot99efP2dd0LXrm7d5ZeX3n/ePBg2LHPxRSkshN69i59rwADo1cstrVszx/2It9JUK0ennlp0DTp0cH9Lkk1H55wDPXvCgQfCKadAQUFqcXg8nkqTyGPmewHPmFl0+8ckzcxINPffD02bwuefw+OPwxVXwNSpxbepXRtuuQX69IHNmyE/H4491mWGc+fCc8/BkiVQrx6sX1/18UW5/XbYf3/46aeiz15/vej1yJE8O23axmPSE1n1cRS77SWXQOM0zzpSWUe33gqNIiMzXXwx3HVXeuPzeDyJk4kh0uc2biz16SN16yZNmqSEGTxYeust93rnTql5cykcLn+fk07S/2/v3IOjqvI8/j1pDJ10gxF5CETMjEYoCGUyjQSmIC2kMIkGsIBBIj7mpU5Zg+OMUzLrTCYhWQt0LcvVcZhidxRnpwyFPIY0ko2QQcUpSETwwTI85OGYlMNDdiMdE5Lu/u4fpwOhk37Q3Te3Sf8+Vb/K7Xv7nnNuf+H87j3n3N+Pb7+tt7/3PXL79tDf37mTnDmTvOsu8tZbyUcfJb3e+Lbviy/I2bPJhgby7rt7H29tJTMyeB2wjyaFsTdNo258PjIzkzxypPf3E0Ejn4/8yU/IVasIYC8TIPWAmFiymSHLzKuysoAPPwT27gVeegn46it9oOfwTk/705/08ZYW4EZ/gOpBg/Tddfe5fXHyJLB/vx7KAYAjR/RTSn4+4HQCH3zQ93lNTXoo8OBB4NgxYNOm+LbviSeA554DUoL8vH/5C1BYiP/VwVVNwTSNutm1Cxg1CsjO7vs8MzX6wQ+AG24ADh0Cli0Lfm2CIBiKIek2lpw+rcfxAT0fcfQocP31wYdaosHtBhYuBF588dKQjMcDnDsH7NmjndPixcDx44AKyPQwdSrw7W/r7bIy4P339XxDPNq3dSswcqQe1nrnnb6/U1MD/PjHwMaNsdcXJaZp1E1Njf7tg2GmRq+9puenli2L7+8hCMIVEX8HpdQdU202PQeQng7ccQfQ0aGP3XsvcPhw73N+8QvgwQeBsWN1Z5mZqZ1Na6vuNAPp6tId39KlwIIFl/ZnZurPSukOLiUFOHsWGBHwCkugw+r+HI/2/e1vQG0tsG2bvu6vvwbuvx/485/18bNn9dPB5s1Bf0LDMVMjQJ+3aZN+ggvexr4/94dGAGCxAEuW6KcsQRDMIe7jhsD8d6+9liTJv/+dHDxYzylEwu9+p+cbSLKmRs8pBeLzkQ88QP7sZ72PrV5Nlpfr7cOH9RxH4PzDzp2k1UoeP67nNe68k9ywIX7tC6wrcH5j9WrywQdJ0ry5DTM1Ism6OrKgIHgdZmnk85FHj17afvJJ8sknZQ5KTMwki3+hwOD3hw4lJ0wg588nnc7IO7/2dnLRIvLmm8nbbyePHdP7W1rIkhK9vWuXbvbkyeRtt2l76y197MIFculSctIkMi9PT4AHEssEfCTtC6wr0EE5nbqDJs10UOZpRJIPPaQddTDM0sjrJb/7XTInR/8buu8+srVVHJSYmElmSD4oiVIQGWZGKBCNIkciSQiCOUiwWEEQBCEhEQclCIIgJCTioPqJyspKKKUuWmVlpdlNSjpEA0G4uhAHdQXE0sGVl5cjJycHWVlZcLlcKC8vN66hSUwojSorK+F0OuF0OkEyvH4ksG+foe0VBCEERqy8cDgcHIh4PB7m5OQwKyuLLpeLHo8n4vMKCwuZkpJCABxqs/GnU6aYujosGTWKSD+PR69C/PnPyZtuInV+LVnFJyZmghlS6EDs/AKdjN1uZ2FhYUROyuVycYTNxrkA/wPgP/W9OQEcpUnCJ5tGIfXr7CTr6/Vy9lGj9H+L1FSytJR89VVxUGJiJpkhhSZy51dRUUFo50AArKioiOg8l8tFu91+2bl2u50ulyv4SadOkX/8Iw9NmMA2v1P6P4A1AMt0GS00Sfhk0yjwmBXgvVYr/zF7NpmRof8r2Gzk4sXkunU6oK8fcVBiYuaYIbH4Epny8nJs3LgRbrcbL7/8MkpKSiI6b//+/Whra7tsX1tbGz766COU9swndPiwTvmxZQuwezdAYtzw4Xj9mmuwsasL7wLoAmC32wG3+5v4XdnAwQiNSCLF7ca9ABYCKAFg7+hA++7dOmbjggU6JUhaWtyvRxCE6EiqRRJerxdFRUU4ePAgTp48ibKyMhQVFcHr9YY9Ny8vDzab7bJ9NpsNuZMn69huTz0FjB+vkyUuXw60twMVFcC+fUj98ktsKCjAX1NSLjqnfB3du9WQC72KibdG49LTMffMGTxSW4vTANYBmAngvwDMtVrx1zfeANauBebNE+ckCImGEY9liTp8FNUwnZ+ecxhpAJdYrawbPZq+ESNIgBw0iJwzR8eC+/zzPs8PnKCHLJLoRTw0ylSKjwHcabHQoxQJ0HfTTVw/bhxnKMWUK5xDNFMnMbFktqQa4ot4mK4PLGfP4u3Fi/FuUxOmu92wdnSA33wDddddwPz5QHFxyOyw1dXVOHBAZ2GfO3cuKioqYr+gAUjUGh07BsumTdjudkNRh+86P2YM1P33A4sWQeXlYYHPh6rcXIzrMXRosViMvBxBEGLBCK83IO7OfT7y4EFy5Upy+nTSfyf+5eDB3DB2LLljh179FQOQJ6heRKyRz0d++im5YoUORutfhMLvfIdrsrL4wJQpvcqOdvGFmTqJiSWzGVJoonZ+YZeKezzke+/pNAvZ2bzY6Tkc5IoV/P2jj0bVwQVDHFRvQmrk85FNTeTy5Zf0UYqcMYN84QXyxImonVAoxEGJiZljSRfN3Ov1Ijc399IKsYICWBoa9Kq7rVt1evBrrgFmz9YT5/Pm6eR3BiDRzPvmMo1efBElQ4bAsmWLTnLY3KxTuc+apVfe3XOPTs9uIBLNXBDMIanmoADAYrEg227H9M5OlP7hDzqN+IULQEYGcPfd2iEVF/dOUS70GxavFyVKYeaFCyh9+GHgzBnAagWKioBnngFKS4Fhw8xupiAIBpMcDooEDh4EamvR/Mor2NTSAgA4ceQITuXnY9rKlcCMGfrJSTCHtjagvh7YtAkdGzbguQsXcB5ADQDLokVY/NprgN1udisFQehHrsr3oMIFBFVKYZBScCqF3dOnA9nZQE4O8PTTyBw9GqiuBj75BN/y+TBtzx49XCTOKa6EC6xbWVmJDKWwVClsVApdGRnAwoVAXR2s990HbN2KIe3tKCOx+M03xTkJQjJixMRWrBPw4Sa6gwb9PH+envXruSUjg+dSUkiAvtRUnep79WqyuTmmeuMNrvJFEqF+r6AanTpFrllDX1ERO/0rI9uHDaP3scfIhgayqyumeo3ATJ3ExJLZDCk0Hp2f0+mk0+nstT9wldct6el8Yfx4+oqLtTMC+BXA1wEutVpZ6nReUdTxaKKVR8vV7qDIvnUK1Gh8ejpfvvVW+goKSP+NQ0taGv8NYD7AITZbxC/NdpefLDqJiSWzGVKow+GI6S43VAfkqq3l7WlpfBpgY/cycIDuG27gsfnzWWy10hJjpIju866k04yGRHBQRujkcrl4W3o6nwrQqHXcOPK3v+U7L71Eu80Wc0SPZNBJTCyZzZBCHQ5H3HIn2e12zpk1i57t28knnuC5665jd4e3B+C/AJwEsLqqilVVVVRKXdbxKaVYXV0dtt5YQuxESyI4qLjpZLPx4alT6f3Nb3hq5MiLGjUCXA5wPHBRB9FJTEwsEjNskUS0AT/r6urQ2NiIdJ8PiwD83u3Gup07YZkzB1i9Gp0334zHBw/GGADTAKwE8Lndjty8vOABXXNzw9YbKsTOQCWWwKx1dXVo2rMHt/t8eBbA/rY2rGlqgnrmGVhGjMBTqakYByAfwLMAWuz2izqIToIgRIQRXu+WW26J7i63uZlvlZZyG8AO/x34WYBrAa5fskQvgog2KV0Yku3O3OFwRHfNXV1kQwObpk5ls1+jCwC3AXwY4AvLl4fVQXQSExOLxAwpdMyYMZEN4fh85Mcfk1VVOpyQv8P7TCk+D3AmQEsfHVDMab37INnmNhwOR+RDbR0d5Nat5A9/SF5/PQnQk5rKLRYLlwK8tg9HEU4H0UlMTCycGVJoyCeozk69nPjxx8msLF6MpzZtGrlyJT2ffMLC2bPDdkDBVvkZtTjDCMx2UCGfRs6fJ998k1yyhBwyROs0dCi5dCm5cSM9ra1hHUUwjUjRSUxMLLwZUqjD4bis8xpts7EqJ4fesrKL6bU7Bw1iLcAfARyFCN+h8WPkezChOtV4Y7aDCnwaybTZuGriRPrmzSOtVhLgaYBrABYDrPr1ry9rfyidjH5XKVl0EhNLZjOkUIfDQc+JE/zX0aP5rtVK76BBuqrhw8nvf5/cvJme1lbTOrdgJNMLoBdX8TU3c8WYMXzPaqXXYtE6jR1LLltGz44dvG3SJFNuFEKRTDqJiSWzGVKoIz2d3fNJ/0hLI3/5S3LXLp3OgubMIyQipjqozEydpsIfzeGLtDSdxqKxkfR6RaMeiIMSEzPHjFlmrhS2FxZiAoBx7e1Qzz+Pyh07AH/20u6l5D6fDwDgdrvR2NiIuro6Q5oj9EFzM/D119hZUIDJAG5sb4d69llUbtsGpKSIRoIgmE5E0cyVUsUA/h2ABcB/klwV8oQJEzBnxw4cCnI4ltTrQt9csUY5OcDHH2MWgE/7OCwaCYJgNmGfoJRSFgCvACgBMBFAmVJqYiyVxvKiptCbqDQaPDjkYdFIEASziWSIbyqAz0geJ9kJYB2A+bFUWlJSgvz8fKSk6Ortdjvy8/NRUlISS7HJjGgkCMKAIxIHNRbAFz0+N/v3XYZS6hGl1F6l1N4zZ86ELNBisaC+vh4TJ05EVlYWampqUF9fD4t/jkq4YkQjQRAGHHFbJEFyDckpJKeMGDEi7Perq6tx4MABnDx5EnPnzkV1dXW8miIEQTQSBOFqIpJFEi0AbuzxOdO/LyYqKyt7ZVkVokY0EgRhwBHJE9QHALKVUt9SSqUCWAKg1thmCVeIaCQIwoAj7BMUSY9S6qcA6qGXML9K8n8Mb5kQMaKRIAgDEUUy/oUqdR7A4bgXHBvDAZw1uxEBjCc5xIyKRaMrwjSdBCGZiehF3Sg4THKKQWVHhVJqbyK2ycTqRaMIMVknQUhaDMuoKwiCIAixIA5KEARBSEiMclBrDCo3FqRNiVN3MBKxTUDitksQBjSGLJIQBEEQhFiRIT5BEAQhIREHJQiCICQkcXVQSqlipdRhpdRnSqlfxbPsWFBKnVRKfaqU+sisJcNKqVeVUqeVUgd67BumlNqulDrq/3tdP7RDNArdjoTQSRCEODooI/JGxZlZJHNNfM9mLYDigH2/AtBAMhtAg/+zYYhGEbEWJuskCIImnk9Qcc9JNJAg+R6AcwG75wN43b/9OoB7DG6GaBSGBNFJEATE10FFlJPIJAjgbaXUh0qpR8xuTA9GkfzSv/1PAKMMrk80io7+1kkQBBgX6ijRmEGyRSk1EsB2pdQh/51ywkCSSqlkXvOf8BoBopMg9CfxfIIyJCdRPCDZ4v97GsBm6KGuROCUUmo0APj/nja4PtEoOvpbJ0EQEF8HlZA5iZRSNqXUkO5tAHcCOBD6rH6jFsBD/u2HAGwxuD7RKDr6WydBEBDHIb4Ezkk0CsBmpRSgr/cNkv/d341QStUAuAPAcKVUM4AKAKsArFdK/QjA5wAWG9kG0Sg8iaCTIAgaCXUkCIIgJCQSSUIQBEFISMRBCYIgCAmJOChBEAQhIREHJQiCICQk4qAEQRCEhEQclCAIgpCQiIMSBEEQEpL/B/1FlfO2dRBTAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import minimize\n", "%matplotlib inline\n", "\n", "# Function which evaluates the chi-squared value of a model y=a*x+b for an (x,y) dataset\n", "def chisqstraightline2(p,x,y,yerr):\n", " a,b = p[0],p[1]\n", " return np.sum(((y - (a*x + b))/yerr)**2)\n", "\n", "# Read in the example dataset sampled from a straight line\n", "# You will need to change the file path to the location where you've saved the data\n", "data = np.loadtxt('../datasets/straightline.dat')\n", "x,y,yerr = data[:,0],data[:,1],data[:,2]\n", "nx = len(x)\n", "# The number of jack-knife samples is equal to the number of data points, we use separate variables for clarity\n", "nsamp = nx\n", "\n", "# For each jack-knife sample, we exclude one of the data points in turn\n", "# We loop over jack-knife samples, plotting the best-fit in each case\n", "fig = plt.figure()\n", "xmin,xmax = 0.,10.\n", "ymin,ymax = -0.5,5.\n", "asamp,bsamp = np.empty(nsamp),np.empty(nsamp)\n", "for isamp in range(nsamp):\n", "# Select jack-knife sample by removing each data point in turn\n", " cut = [np.arange(nx) != isamp]\n", " x1,y1,y1err = x[tuple(cut)],y[tuple(cut)],yerr[tuple(cut)]\n", "# Fit straight line to data by minimising the chi-squared value\n", " a0,b0 = 0.2,1. # initial guess (result should not depend on this)\n", " p0 = [a0,b0]\n", " result = minimize(chisqstraightline2,p0,args=(x1,y1,y1err))\n", "# Best-fitting values of straight line model parameters a,b for this sample\n", " a1,b1 = result['x']\n", " asamp[isamp],bsamp[isamp] = a1,b1\n", "# Plot each sample and best fit\n", " sub = fig.add_subplot(3,4,isamp+1)\n", " plt.plot(x1,y1,marker='o',markersize=5,color='black',linestyle='None')\n", " plt.errorbar(x1,y1,yerr=y1err,color='black',capsize=2.,linestyle='None')\n", " plt.plot([xmin,xmax],[a1*xmin+b1,a1*xmax+b1],color='red')\n", " plt.title('Re-sample ' + str(isamp+1))\n", " plt.xlim(xmin,xmax)\n", " plt.ylim(ymin,ymax)\n", " xc,yc = 0.5*(xmin+xmax),ymin+0.8*(ymax-ymin)\n", " caption = 'a=' + '{:4.2f}'.format(a1) + ' b=' + '{:4.2f}'.format(b1)\n", " plt.text(xc,yc,caption,horizontalalignment='center',color='red',fontsize=10)\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here are the best fits of $(a,b)$ for each of those samples, compared to the original $\\chi^2$ contours:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$b$')" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeVhUZf8G8Pth2HcVUEEQUQTEJRXXcl+zXs2y1Eotd3PJFn0z28w0zdTMNDWXXPK1TMu13MMlyX0Bc9/BNRWJnZn794c4kT8XPHOGYfD7ua65ZGD4nocTzc05z6ZIQgghhLCUg60bIIQQomiQQBFCCKELCRQhhBC6kEARQgihCwkUIYQQupBAEUIIoQubB4pSKlgptUkpdUgplaCUev0ur1FKqS+VUseVUgeUUjVs0VYhhBD35mjrBgDIAfAWyT1KKS8Au5VS60geyvOaJwGE5z7qAPg6918hhBCFhM2vUEheILkn9+MUAH8CCLrjZe0AzOMtcQB8lVKlC7ipQggh7qMwXKGYKaVCAVQH8McdXwoCcC7P8/O5n7twx/f3BtAbADw8PGpGRkZaq6lCCFEk7d69+ypJfy3fW2gCRSnlCWAJgMEkb2qpQXIGgBkAEBMTw127dunYQiGEKPqUUme0fq/Nb3kBgFLKCbfC5DuSS+/ykkQAwXmel8n9nBBCiELC5oGilFIAZgH4k+SEe7xsOYCuuaO96gJIJnnhHq8VQghhA4XhltfjALoAOKiU2pf7uXcBhAAAyWkAVgNoA+A4gDQAr9qgnUIIIe7D5oFCcisA9YDXEED/gmmREEIILWx+y0sIIUTRIIEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdOFo6wYAgFJqNoCnAVwmWfkuX28MYBmAU7mfWkry44Jrof37+++/cf36dWRlZSErKwtGoxFubm7w8PCAp6cnPDw8oJSydTOFEHasUAQKgG8BfAVg3n1es4Xk0wXTHPuVkZGBnTt3Ytu2bdi+fTtOnjyJc+fOITk5+b7f5+zsjFKlSqFUqVIIDAxEWFiY+REVFYWQkBA4OMgFrRDi3gpFoJDcrJQKtXU77FlsbCwmTZqEVatWISsrCwAQGRmJyMhING7cGMHBwShevDhcXFzg5OQEBwcHpKenIzU1Fampqbh69SouXryIixcv4siRI/j111+RkZFhru/u7o7IyEhUrVoVNWrUQExMDKpVqwZ3d3db/chCiEKmUARKPtVTSu0HkATgbZIJtm5QYbB//34MGTIE69atg7+/P/r164dmzZqhfv36KFGihOa6JHHx4kWcOHECf/75Jw4dOoSEhASsWrUK3377LQDA0dER1apVQ/369VG/fn088cQTKFOmjE4/mRDC3iiStm4DACD3CmXlPfpQvAGYSP6tlGoDYBLJ8Lu8rjeA3gAQEhJS88yZM9ZttI198803GDBgALy9vTFs2DD069cPbm5uVj0mSSQmJmL37t3YsWMHtm/fjj/++ANpaWkAgHLlyqFBgwZo0qQJmjdvLgEjhJ1RSu0mGaPpe+0hUO7y2tMAYkhevddrYmJiuGvXLt3aV9gsXLgQL730Elq1aoXvvvvOoqsRS+Xk5ODAgQPYsmULNm/ejC1btuDKlSsAbt12a9myJVq3bo3GjRtbPfCEEJYp8oGilCoF4BJJKqVqA/gRQFnep/FFOVD27duH+vXro1atWli/fj2cnJxs3aR/IYmDBw9i3bp1WLduHWJjY5GRkQFXV1c0atQITz/9NJ5++mmEhobauqlCiDvYfaAopf4HoDEAPwCXAHwIwAkASE5TSg0A0A9ADoB0AG+S/P1+NYtqoJBEtWrVcO3aNezevRslS5a0dZMeKD09HZs3b8Yvv/yC1atX49ixYwCAKlWqoF27dmjfvj2qV68uw5aFKATsPlCsoagGys6dO1G7dm1888036Nmzp62bo8nRo0exYsUKLF++HFu3boXJZELZsmXx7LPPokOHDqhbt64MURbCRiwJFPm/1s789NNPcHR0xHPPPWfrpmhWsWJFvPXWW4iNjcXFixcxa9YsVK5cGVOmTMHjjz+OsmXL4o033kBcXByK6h88QhRFEih25ujRo6hQoQKKFStm66bowt/fH927d8fKlStx+fJlzJ8/HzVr1sTXX3+NevXqISwsDMOGDcPBgwdt3VQhxANIoNiZ8+fPIygoyNbNsAofHx+8/PLL+Pnnn3Hp0iXMnTsXkZGRGDduHKpWrYpq1aph3LhxSExMtHVThRB3IYFiZ0jC0dGe5qNq4+Pjg65du+KXX37BhQsXMGXKFLi7u2Po0KEIDg5Gy5Yt8d133yE1NdXWTRVC5JJAsTPFixfHtWvXbN2MAuXv74/XXnsN27dvx7Fjx/D+++/j2LFjePnll1GqVCn07NkT27dvl/4WIWxMAsXOlClTBidOnHhk3zwrVKiAESNG4MSJE4iNjcXzzz+PRYsWoX79+oiOjsb48ePNkyqFEAVLAsXOxMTE4Nq1azh16tSDX1yEOTg4oGHDhpg9ezYuXLiAmTNnwsfHB2+//TaCgoLwwgsvYO3atTCZTLZuqhCPDAkUO9OgQQMAwNq1a23cksLDy8sLPXr0wPbt2xEfH48BAwZg48aNaNWqFSpWrIhx48bh6tV7rtIjhNCJBIqdiYqKQlhYGH7++WdbN6VQio6OxoQJE5CYmIiFCxciKCgIQ4cORVBQELp06SJ9LUJYkQSKnVFK4YUXXsD69euRlJRk6+YUWi4uLujcuTNiY2MRHx+P3r17Y9myZahfvz5q1qyJWbNmmVdIFkLoQwLFDvXs2RMmkwnTp0+3dVPsQnR0NCZPnozExER8/fXXyM7ORs+ePREcHIyhQ4c+8v1RQuhFAsUOlS9fHm3atMG0adOQnp5u6+bYDS8vL/Tt2xcHDhzAb7/9hqZNm2LChAkoX7482rdvj82bN8vtMCEsIIFip4YOHYrLly9jxowZtm6K3VFKoVGjRli8eDFOnz6NYcOGYfPmzWjUqBFiYmKwYMEC8zbKQoj8k9WG7ViTJk1w5MgRHDt2DB4eHrZuzl2lpqbizJkzSEtLQ2pqKjIyMuDm5gYvLy94eXkhICAA3t7etm4m0tLSsGDBAkycOBGHDx9GYGAgBg4ciD59+hSZddOEyA9Zvv4uHoVA2bZtG5544gmMGDECH3zwga2bAwC4dOkSFi5ciG3btuHAgQM4fvz4A28jeXt7IyQkBKGhoYiMjERkZCQqVaqEKlWqwNPTs4BafovJZMKvv/6KiRMnYv369fDw8ECPHj0wePBglCtXrkDbIoQtSKDcxaMQKADw/PPPY/Xq1Th8+DCCg4Nt1o6zZ89iyJAhWLJkCYxGI8qXL49q1aqhWrVqCA8Ph5eXF9zd3eHq6or09HSkpKTg5s2buHTpEs6dO4ezZ8/i1KlTOHLkCDIzMwHcujUVFRWFmJgY1K1bFw0bNkSlSpUKbCOu/fv3Y8KECVi4cCFMJhOee+45DB06FDExmv5fE8IuSKDcxaMSKKdOnUJ0dDRatGiBn3/+2Sa7Hq5duxbPPvssSKJ///7o3r07IiMjNdUyGo04ffo0EhISsHfvXuzatQs7d+7EpUuXAAB+fn5o2LAhmjZtimbNmiEiIsLqP3NiYiK+/PJLTJs2DTdv3kSzZs3wzjvvoFmzZrLLpChyLAkUkCySj5o1a/JR8dlnnxEAf/zxxwI/9ubNm+nq6sqqVavy9OnTVjmGyWTi8ePHOXv2bL7yyissW7YsARAAAwMD2b17dy5dupQpKSlWOf5tN27c4NixY1m6dGkCYM2aNbl48WIajUarHleIggRgFzW+78oVShGQk5OD2rVrIzExEQkJCfDz8yuQ45JETEwMrl+/jh07dhTocU+dOoUNGzZg3bp1WLt2LZKTk+Hs7IymTZvimWeeQdu2bVG6dGmrHD8zMxPz5s3DuHHjcOzYMURGRuKdd97Biy++CCcnJ6scU4iCIlcoj/gVCknu37+fTk5OfOGFFwrsmBs3biQAzpgxo8COeTdZWVnctGkT33zzTZYvX9589VK3bl2OHz+eZ8+etcpxc3Jy+P3337NatWoEwLJly3LKlClMT0+3yvGEKAiw4ArF5m/81no8aoFCkp988gkB8Ntvvy2Q4w0aNIgeHh6F6g3UZDIxPj6en3zyCatXr24Ol/r163Py5Mm8dOmSVY65cuVK1qtXjwBYqlQpfv755/z77791P5YQ1iaBIoFC8tZfzI0bN6a7uzsTEhKsfryOHTuyYsWKVj+OJY4ePcpRo0axSpUqBECDwcBWrVpx3rx5uve5mEwmbty4kU2aNCEA+vn5cfTo0bx586auxxHCmiRQJFDMEhMT6e/vz6ioKKt3Urdu3Zr2dJ4PHjzId999l6GhoQRAT09Pvvrqq9y8eTNNJpOux9q6dStbt25NACxevLgEi7AbEigSKP+ybt06KqX40ksv6f5GmVefPn1YokQJq9W3FpPJxM2bN7N79+709PQkAIaHh3P06NE8f/68rsfauXMnn376aQkWYTckUCRQ/p+RI0cSACdNmmS1Y0ycOJEAmJiYaLVjWNvff//Nb7/9lo0aNSIAOjg48D//+Q9XrVrFnJwc3Y6zY8cOtmnTxnwr7LPPPpM+FlEoSaBIoPw/RqOR7dq1o4ODA9esWWOVY+zbt69QjPLSy/Hjxzls2DCWLFnSPGpr9OjRunbkx8XFsVWrVgTAkiVL8osvvihUgxqEkECRQLmrlJQUVq1alT4+Pjx06JDu9U0mE8PCwtiiRQvda9tSZmYmf/jhB3PnupOTE1966SVu375dt1uIW7duNdcPCgri119/zczMTF1qC2EJCRQJlHs6ffo0S5YsybCwMF6+fFn3+h9++CEB8OjRo7rXLgwOHTrEgQMH0svLiwAYExPDuXPnMiMjQ5f6GzduZP369QmA5cqV47x583S91SbEw7L7QAEwG8BlAPH3+LoC8CWA4wAOAKjxoJoSKP+Ii4ujq6sr69evz7S0NF1rX7hwgU5OTuzZs6eudQubmzdvcsqUKYyKiiIA+vv78/333+fFixctrm0ymbhq1SrzvJlKlSpx6dKlVh1QIcS9FIVAaQigxn0CpQ2AX3KDpS6APx5UUwLl3xYvXkylFNu2bcvs7Gxda/fr14/Ozs48deqUrnULI5PJxHXr1vE///kPlVJ0cXFhz549dbmlaDQauXjxYkZGRhIAa9euzQ0bNujQaiHyz+4D5dbPgND7BMp0AJ3zPD8CoPT96kmg/H+TJ08mAPbs2VPXv37PnTtHNzc3du7cWbea9uDIkSPs27cvXV1dCYBPP/00Y2NjLT632dnZnD17NoODgwmAzZs35+7du3VqtRD39ygEykoAT+R5vgFAzF1e1xvALgC7QkJCdDq9Rcvw4cMJgO+8846udd977z0CeCT/or58+TJHjBhBPz8/AmCdOnW4dOlSi1chTk9P54QJE1iiRAkCYOfOnXny5EmdWi3E3UmgyBVKvplMJvbt25cAOHbsWN3qpqWlsUKFCixfvjxTU1N1q2tPUlNTOXXqVIaFhREAIyMjOWvWLIs78G/cuMF3332Xbm5udHJy4uuvv86rV6/q1Goh/u1RCBS55aWjnJwcdurUiQA4ZcoU3ereXn148ODButW0R9nZ2Vy0aBEfe+wxAmCZMmX45ZdfWjwg4vz58+zZsycdHBzo4+PDsWPHyhwWobtHIVCeuqNTfseD6kmg3F9WVhbbtm1LAJw9e7ZudQcMGEAAXL9+vW417ZXJZOKvv/7KBg0aEAADAgI4duxYi9dYi4+P51NPPUUADAkJ4YIFC2STL6Ebuw8UAP8DcAFANoDzAHoA6Augb+7XFYApAE4AOPig212UQMmX9PR0tmzZkkopzps3T5eaqampjIiIYOnSpXnhwgVdahYFsbGxbNmyJQGwRIkS/PTTTy1e02vDhg2sUaOGeX5MbGysTq0VjzK7DxRrPCRQ8ic1NZXNmjXTNVT2799PNzc3NmzYkFlZWbrUfBhZWVn86aef2KNHDzZt2pRhYWF0cnKiq6sr/fz8GBoaylq1arFz5858//33OW/ePB46dKhA/sqPi4vjk08+qVuwGI1Gzps3j2XKlCEAtm/fnseOHdOxxeJRI4EigWKRvKEyZ84cXWouWLCAAPj666/rUi+/Lly4wMqVK5tX961bty47derEIUOGcMiQIXzttdfYpUsXtmzZkmFhYXRwcDBvwuXt7c3mzZtz5MiRjIuLs+qM9bzBosdikampqRw5ciQ9PDzo5OTEN998k9evX9exxeJRIYEigWKx1NRUtmjRgkopfvPNN7rUHDx4MAFw1qxZutR7kMTEREZERNDDw4Pff/99viZwZmZmMj4+nnPmzGG/fv3M2/kCoK+vLzt06MA5c+boMiP+buLi4sz7puixWGRSUhJ79OhBpRT9/Pw4depU3SeyiqJNAkUCRRfp6enmv5q/+uori+tlZ2ezRYsWdHJyKpD7+0888QQ9PT25ZcsWi+pcvnyZ//vf//jqq6+ydOnS5oCpU6cOJ0yYYJXl+rdu3cqmTZuaF4ucPn26RbcL9+7da16Sv3Llyly3bp2OrRVFmQSKBIpuMjIyzKO/xowZY3G9a9euMTIykr6+voyPj9ehhffm4eHBQYMG6VrTZDJxz549HDlypLkDXCnFxo0b85tvvtH9ttKGDRvMe9NXqFCBCxcu1Ny3YzKZ+OOPP7JcuXIEwGeeeYbHjx/Xtb2i6JFAkUDRVVZWlnmeyvDhwy1eSuT06dMsXbo0y5Qpw7Nnz+rUyn/LyMggAL7//vtWqX/bkSNH+NFHHzE8PJwA6OLiwg4dOnDFihW63VoymUxcuXIlq1atSgCsWrUqV61apfm/Q3p6OkePHk0PDw86Oztz2LBhVt8eWtgvCRQJFN3l5OSwV69eBMC+ffta3EG9b98+ent7MzIy0irL6BuNRgYGBvLVV1/VvfbdmEwm7tixg4MGDaK/v795AuOHH36oW2gajUYuXLiQ5cuXJwA2btyYO3bs0FwvMTGRXbp0Md9WW7hwoaxoLP4fCRQJFKswmUx85513CIDPP/+8xUuIxMbG0tXVlTVq1OCNGzd0auU/WrZsyUqVKule90GysrK4ZMkStmrVikopOjg4sG3btly7dq0ub9iZmZmcPHmyObg6duzIEydOaK73+++/m2/fNWjQgPv377e4jaLokECRQLGqzz//nADYtGlTJicnW1Rr5cqVdHR05OOPP677bZdJkyYRAOPi4nSt+zBOnjzJYcOGmd/8IyIi+NVXX+nysyYnJ/P999+nu7u7eWjwtWvXNNXKycnhjBkz6OfnR4PBwEGDBskwY0FSAkUCpQDMnTuXBoOB1atXt3gG/OLFi+ng4MDGjRvrupBkSkoK/fz8+OSTT+pWU6uMjAzOmzePtWrVIgD6+PjwzTff1GW14MTERPPQ4GLFinHixImatw++du0aX3vtNTo4ODAgIIBz5syRZVwecRIoEigFYvXq1XR3d2e5cuV45MgRi2otWLCASik2adLEogl9dxozZgwBFJphsiaTib///js7depER0dHGgwGdu7cmfv27bO49v79+9miRQsCYHh4OH/++WfNt9j27NljHl32xBNP8MCBAxa3T9gnCRQJlAKzY8cO+vv7s0SJEvz9998tqjV//nw6ODiwQYMGFq9rdVt6ejorVKjA8PDwQrcS77lz5zhkyBB6enoSAFu3bs3NmzdbVPP29sG3d3ls1qyZ5j4Ro9HIWbNmsUSJEjQYDHzrrbd0++8i7IcEigRKgTp27BjLly9PV1dX/vjjjxbV+v7772kwGFi7dm3+9ddfurRv7dq1BMAhQ4boUk9v165d46hRo8z9LA0bNuS6dess6sDPysri5MmTWbx4cTo4OLBv3768cuWKplpXr141j/ALCgrikiVLZDTYI0QCRQKlwF2+fNl8i2Ts2LEWveEsW7aMzs7OrFq1qm5LnNzeRGzlypW61LOG1NRUfvHFFwwMDCQA1qtXj2vWrLHoXP711198/fXXaTAY6Ovry0mTJmmecb99+3bzUjRPPfUUT506pbldwn5IoEig2ERaWhpfeOEF8z71liwVsm7dOrq7uzM8PFyXN660tDRWq1aNxYsXt2iIbUFIT0/nlClTzHvI161b1+IhxwkJCeb+lejoaM1bM2dnZ3P8+PH08PCgm5sbx44da5MVpEXBkUCRQLEZo9HId999lwDYpEkTi25bbdu2jb6+vgwKCmJCQoLFbTt27BiLFSvGSpUqWWXei94yMjI4bdo0c7A0bNjQojXQTCYTf/75Z/PSKy+88ILmSZdnz55lu3btzDP3t2/frrldonCTQJFAsbl58+bR2dmZ4eHhPHz4sOY6Bw4cYOnSpVm8eHFu27bN4nZt3LiRjo6ObNWqleahtQUtIyODkydPZqlSpQiALVu25K5duzTXS0tL44gRI+jq6kp3d3d++umnms/FTz/9xKCgICql2L9/f4vnJYnCRwJFAqVQ2LJlC/38/Ojr68u1a9dqrnPixAlWqFBBl05/kpw5cyYBsFOnTlbd40RvqampHDduHEuUKGGeIW/J4o6nTp3iM888QwCMjIzkxo0bNdVJTk7moEGDqJRiUFAQly1bprlNovCRQJFAKTROnTrFKlWq0GAwcNKkSZr7Aa5cucJ69epRKcUJEyZYPMpo7NixBMDevXvb3YilGzdu8L333qO7uzsdHR05cOBAzSO4SHLVqlUMCwsjAHbu3FnzRNW4uDhWqVKFANihQwfZ8rmIkECRQClUbt68ab7f3qNHD81rgKWlpfG5554jAPbr18/i1XyHDRtGABw0aFDBhsrly+SOHbf+tUBSUhL79OlDBwcH+vj48PPPP7fo3H7wwQd0dnamj48Pp06dqmmGfFZWFkeNGkUXFxcWK1aMs2fPtrvAFv8mgSKBUugYjUa+9957BMD69etr/uvVaDRyyJAh5omAltyzN5lMfOONNwiAb731VsG88S1cSLq5kT4+t/5duNDikvHx8eaN0MqXL8+lS5dq/lmOHDnCZs2amTcQ27t3r6Y6hw8f5hNPPEEAbN68uQwxtmMSKBIohdYPP/xAd3d3BgYGWjQyaMaMGXR0dGR0dLRF62GZTCYOGDBAt2X57+vy5VshAvzzcHOz+Erltl9//ZXR0dHmGfJaNzAzmUxcsGABAwICaDAYOHToUE1rrBmNRk6ZMoWenp708PDg5MmTZV0wOySBIoFSqO3bt4/lypWjs7OzRfvVb9iwgb6+vvTz87NoyZK8y/J36tTJeqO/duy4dWWSN1C8vW99XifZ2dmcPHkyfX19LV41+K+//mLPnj0JgOXKleOaNWs01Tlz5gxbtWplXhfM0nXfRMGSQJFAKfT++usvtmzZ0jwJUus6W0eOHGHFihXp5OTEGTNmWNSm2wtJtmjRwjprVln5CiWvK1eusG/fvlRKMSAggHPnztV8G2zTpk2sWLEiAbBbt26a5haZTCZ+++239PX1paurK8ePH29XI+weZRIoEih2IScnx9wxXrNmTZ4+fVpTnevXr5v/Au7fv79FM7dnz55Ng8HAxx57jImJiZrr3NPtPhRvb936UO5n9+7drFu3rsWrBqenp/Pdd9+lo6MjAwIC+MMPP2gKqKSkJP7nP/8xLy1jyRwlUTAkUCRQ7MrPP/9Mb29vFi9enL/88oumGjk5OXzrrbcIgI0aNbJoDbBffvmFHh4eDA4Ots7uhTqN8sovo9HImTNnmlcNHjJkiOYtAvbu3Wve3bF9+/aaBlfc7qMpVqwYXV1dOWHCBOlbKcQkUCRQ7M7Ro0dZpUoVKqX4wQcfaL4dMn/+fLq6ujIoKMiinRr37NnDwMBAenl5cfXq1ZrrFCZXr14194mULVuWq1at0lQnOzubY8eOpYuLC4sXL84FCxZYfLXSoEEDiyZpCuuRQJFAsUupqans1q2beajppUuXNNXZs2cPQ0ND6ezszK+//lpz38G5c+f42GOP0cHBgV988UWRmU+xefNmVqpUyTwIQet5/vPPP80rTLdt25ZJSUkPXeN234qPjw/d3d05derUInOeiwq7DxQArQEcAXAcwDt3+forAK4A2Jf76PmgmhIo9sFkMnHmzJl0dXVlYGCg5sUQ//rrL/PcjJdeeknzLZ6UlBTz8iTdu3fXPHGwsMnMzOTHH39MZ2dnFi9enPPmzdP0Rp6Tk8Px48fT1dWVxYsX58KFCzXVOXfunHk15JYtW/LcuXMPXUNYh10HCgADgBMAwgA4A9gPoNIdr3kFwFcPU1cCxb7s27eP4eHhNBgMHD16tKZ77EajkSNHjqRSipUqVeKhQ4c0tcVoNPL99983T8q0Sme9jRw6dIj169cnALZp00bzG/nhw4dZp04d87IrWpaCMZlMnDp1Kt3d3enr68uFVh6wIPKnQAIFQEMAmwEkAFgIoLbWg95Rtx6ANXmeDwMw7I7XSKA8ApKTk9mxY0cCYKtWrXhZYyf2unXr6O/vTw8PD86fP19ze77//nt6eHgwICCAmzZt0lzHYjp36ufk5HDSpEl0d3ent7c3Z82apekqIzs7m59++imdnJxYsmRJzZuZHT161DwyrVOnTrx27ZqmOkIfBRUoxwG0AOCfe4vqdwAvaD1wnrodAMzM87zLneGRGygXABwA8COA4HvU6g1gF4BdISEhep9nUQBMJhOnTZtGFxcXli5dWvOKuOfPn2eDBg3M8160zPwmby1zEhERQYPBwLFjx1pndFLewLgzPKywdMttx48fZ6NGjcxXK1r6REhy//79rFq1KgGwV69eTElJeega2dnZ/OSTT+jo6MigoCCuX79eU1uE5QoqUOLueO4JIF7rgfPUyU+glADgkvtxHwAbH1RXrsGqK9cAACAASURBVFDs2759+xgREWEeBaZlYcjs7GzzvJfo6GjNS5PcvHmTzz//vHkrXEtW+v1/8gaGkxPp7PxPeEybZvnEyAdc3RiNRn755Zd0c3Nj8eLF+f3332v6MTIyMvjf//6XSimWL19e84i7Xbt2MTIy0rzeWlHpw7InVg0UAPMADAbwOYAPADjmft7JkgPnqf/AW153vN4AIPlBdSVQ7F9KSgpfeeUV8zBTrbsNrlmzhgEBAXRzc+OMGTM03d4xmUz86quv6OzszKCgIIuWfjG720z6vA8XF9LL6/9//p138lf/Ia5uDh8+zNq1a5sHNWjd4TI2NpYhISE0GAwcMWKEpj8EUlNT2a9fPwJgtWrVdNm9U+SftQOlEYBBAGYD2A3gNID1uR3p47QeOE99RwAnAZTL0ykffcdrSuf5uP2dV0t3e0igFB3z58+np6cnfX19NW+4lZSUxObNm5s7kbXep9+zZw8rVKhABwcHzVdOZndb6yvvw8vrVqjkPs8BuALgxwBXdOjAnPtNMtSw7Et2djZHjBhBg8HAsmXLcuvWrZp+rBs3bvCll14iAD7++OOaVx5evnw5/f396ebmxmnTpsnw4gJSILe8zN9wKwAqA3hZj0DJrdkGwNHckBqe+7mPAbTN/fjT3MEA+wFsAhD5oJoSKEXLsWPHWKtWLfN9ei3Dgo1GI8eMGUNHR0cGBwdrvsq4efMmu3TpYh4Fpnn14wddobi53boayQ2TZgA9Aarcf5spxZwFC+5e24KFKePi4hgWFmZxaH733Xf09vamj48PFy1apKnGhQsXzGvAtW/fnlevXtVUR+RfgQaKvTwkUIqezMxM8336iIgI7t69W1OdHTt2mK8yhg8frnktsNtvmF5eXtoXY8y71tftPpS8634dOkQ6OnJFboggz8MT4Apn57tfdVi4MOXNmzfZtWtXIndNMK3Di0+ePGkewdWjRw/NfwiMHz+eTk5ODAoK0jxXSeSPBIoEyiNlw4YNDAwMpJOTEz/77DNNI69u3rzJ7t27EwBjYmI0L7F+6tQp82iyZ599VluH/b1Ged0OG1dXfpx7ZZI3UBTAkS4u977q0GFhygULFtDT05PFixfn8uXLH/5n461dHYcPH06lFKOionjw4EFNdXbv3s3w8HA6ODjwo48+ktWLrUQCRQLlkXP16lU+++yzBMDGjRvzzJkzmur8+OOPLF68uEXLgOTk5HDs2LHm+Rha33j/5Y4rjIe+Qslbx8I5LEeOHGH16tUJgEOGDNF8C2z9+vUsWbIkXV1d+c0332g613mvnBo2bMjz589raou4NwkUCZRHkslk4qxZs+jp6UkfHx9+9913muqcP3/efJ++devWmmfG79u3zzwfo2vXrpZN0LujDyQHYDODgZ6Ojv/0oTg43LsPRWfp6ens27ev+Y1c65yVixcvmgdHvPjii5rmrJDkvHnz6OHhQT8/P80rVou7k0CRQHmknThxwrycSMeOHTV13N4eFnx7PobWTuTMzEy+//77NBgMDAwM5LJlyzTVuVsfSI6rK1csWMCR77zDFePH33+Ul5XMnz+f7u7uLFmyJLds2aKpRk5ODkeOHEkHBwdGRERonh90+PBhc4D/97//tWzEnTCTQJFAeeTl5OSYZ1qXLl1a8xL0eedjPP/885onMe7atYtVqlQhAHbu3FnbMjIFvDlXfsXHxzM8PJxOTk4Wbem8adMmlixZku7u7pw3b56mGmlpaezdu7fFgwfEPyRQJFBErj179jA6OpoA2Lt3b01b+2ZnZ3P06NF0cnJiQEAAlyxZoqktmZmZ/Oijj+jk5MQSJUpoW+G3gDfnyq9r166ZbxMOGDBA89VBUlISGzZsSADs06eP5pnx3333nfkW2Jo1azTVELdIoEigiDzS09P59ttvUynF0NBQzQs7HjhwwNwZ3bFjR81XK/Hx8eahsy1btiwyG0tlZ2fzzTffJHIX80xOTtZcZ+jQoQTA2rVra14R4c8//2R0dDSVUjIKzAISKBIo4i62bNnC8uXLEwAHDRqkaYHIrKwsjhw5kk5OTvT399e8t3pOTg6/+uorenl50dXVlaNGjWJmZuZD1ymMZsyYQUdHR0ZHR/P06dOa6yxZsoSenp708/PTvCjo33//bR4F1qpVK33XXXtESKBIoIh7+PvvvzlgwAACYPny5TXPjj9w4ABjYmIIgM8884zmUU7nz5/nc889RwCMjIzU/MZZ2Kxfv54+Pj4sVaqU5gmn5K2rjKioKBoMBk6ePFnzumszZsygs7MzQ0JCuCMfqwOIf0igSKCIB9i4cSPLlStHpRRff/11TTO2s7Oz+dlnn9HV1ZU+Pj6cOXOm5vWlVq1axXLlypk77YvCJl7x8fEMCQmhh4cHf/31V811kpOT2bZtW/Pseq39Kjt37mTZsmXp7OzM6dOna27Po0YCRQJF5ENKSsq/rla09q0cPXrU3JHcuHFjHj16VFOdtLQ0fvDBB3RxcaGnpyfHjh1r97fBkpKS+Nhjj9HR0dGiHRjz7pr5+OOP89KlS5rqXL16la1atSJwa18cWQ7/wSRQJFDEQ/jtt9/MfSv9+vXTNBLMaDRyxowZ9PHxoYuLCz/55BPNYXD8+HH+5z//IQBWrFiRq1at0lSnsLhx4wYbN25MAJw0aZJFtRYtWkRXV1eWLVuWBw4c0FQjJyeH7777rrnTX4YW358EigSKeEh///0333jjDSqlGBwcrHneSlJSEjt06EDg1iZeWif7keTq1atZsWJFAuCTTz7JQ4cOaa5la+np6Wzfvj0BcOzYsRbV2rFjB0uXLk0vLy/N/51IcunSpfT09GRAQIDmpfkfBRIoEihCo99//51RUVEEwJdfflnzqKDly5czJCTEfGvlr7/+0lQnMzOTn3/+Ob29vWkwGNi/f3+7HamUnZ3NTp06EQDHjx9vUa1z586xevXqdHBw4FdffaW5zqFDh1ihQgWLJ2UWZRIoEijCAhkZGXz//ffp6OhIf39/Lly4UFNne0pKCt9++20aDAb6+flxzpw5mjvtL1++zH79+tHBwYE+Pj4cM2YM09LSNNV6iIPqPokyOzvbvH2yJUFA3jq/tzvr33jjDU2rTJP6TcosqiRQJFCEDvbv32/exOvJJ5/UPKdi3759rFevnnk5EK33/slbI6fatGlDACxTpgxnz55tnQl7D7Fd8MPKyspiu3btCIBz5861qFZOTg4HDRpk3i5Aa8jmnZTZvHlzyxbyLGIkUCRQhE5ycnL4xRdf0MPDgx4eHpw4caKmN3Cj0chZs2axRIkSNBgMfOONNzTPJCdvrXt1O+yio6O5ZMkS/bbEtXAzrvxIT09ns2bNaDAYuGLFCovrTZw4kUop1qtXT/PtRZKcPXs2nZycGB4ernlPnKJGAkUCRejs9OnT5iuDmjVrap6sd/XqVfbq1YtKKZYqVYrz58/XHAQmk4k//PADIyIizO1avXq15cFiwXbBDyMlJYU1a9akh4cH9+zZY3G9H3/8kS4uLqxUqZJFI7e2bNlCPz8/+vr6cv369Ra3y95JoEigCCswmUxctGgRS5YsSQcHB77xxhua9+/YsWOH+Qrj8ccft+gNNTs7m3PmzGFoaChv72u/bt067cFSAFcotyUlJTE4OJiBgYG6TObctGkTvb29GRwczMOHD2uuc/LkSUZHR9PR0ZEzZsywuF32TAJFAkVY0bVr18xLpJcpU4Y//fSTpjpGo5EzZ86kv78/lVLs3bu3tmXtc2VmZnLatGksU6aMOajWrFlj+d72Vl4qf//+/XR3d2eDBg106RDfu3cvAwIC6O/vz71792quk5yczNatWxMA3377bc2d/vZOAkUCRRSA33//3bzHSdu2bTV32l+/fp2DBw+mwWCgr68vv/jiC2ZlZWluV0ZGBqdOnWoOljp16nDlypWFeqn8+fPnEwCHDx+uS70jR44wODiYPj4+/P333zXXyc7O5muvvUYAbN++vaYFRe2dBIoEiiggWVlZHDduHN3d3enu7s4xY8ZoniGfkJBgHr4aGRlp0aQ98lawTJs2jWXLliUAPvbYY1y0aFGhXca9e/fudHBw4B9//KFLvTNnzrBChQr08PDgb7/9prmOyWTiF198QaUU69ata7fzgLSSQJFAEQXszJkz5pngUVFRmtcFM5lMXLFiBcPDw81LrmvdEve2rKwszpkzx9x5X6FCBU6fPt3681geUnJyMoOCgli5cmXd1jBLSkpiVFQU3dzcuG7dOotqLV26lK6urgwPDy8ye9jkhwSKBIqwkZUrV+qyanBmZiYnTpxIX19fGgwG9u3bV/OCiLfl5ORwyZIl5mX3/fz8OHz48EK1svGyZcsIgF988YVuNS9dusSqVavS1dXV4lFb27ZtY/HixRkQEGDRsvz2RAJFAkXY0J2rBn/++eea+0SuXLnCgQMH0tHRkV5eXhw1apTFVxYmk4mbNm3iM888Q6UUHR0d2blzZ27fvl2/uSwWaNKkCQMCAjSPoLubK1eusHLlynRzc7Po9hdJHj58mGXLlqWnp6fFVz32QAJFAkUUAsePHzfPXYmKirLozefw4cPm2eVlypThnDlzdOkLOXHiBAcPHkxvb28CYI0aNThr1iybdj7//vvvuizNcqdLly4xKiqKHh4e3L59u0W1EhMTWaVKFTo5OXHRokU6tbBwkkCRQBGFhMlk4vLly83L4z/77LMWbYv722+/meevVKlShb/88osuVxU3b97klClTWKlSJQKgj48PBw4caHH/jVbh4eFs06aN7nWTkpJYvnx5FitWjAcPHrSo1vXr19mgQQMqpThlyhSdWlj42H2gAGgN4AiA4wDeucvXXQB8n/v1PwCEPqimBIqwpfT0dH7yySd0d3enq6srP/jgA81XASaTid9//z3DwsIIgE2aNGFcXJwu7TSZTNy8eTNffPFFOjs7EwBjYmI4adIki+bIPKwBAwbQ3d3dKiPSTp48ycDAQJYuXZqnTp2yqFZaWpp575qRI0cWiluGerPrQAFgAHACQBgAZwD7AVS64zWvAZiW+3EnAN8/qK4EiigMzp49a17CPTg4mIsWLdL8JpSZmckvv/yS/v7+BMB27dpZ/Fd3XleuXOH48eP52GOPEQAdHR351FNPceHChZq2TH4Yn3zyCQFYNB/nfuLj4+nr68uKFStaHJTZ2dns2rWreQJkUQsVew+UegDW5Hk+DMCwO16zBkC93I8dAVwFoO5XVwJFFCabN282v1E3aNDAohFDKSkpHDlyJL29vamU4ssvv8xjx47p2Fry4MGDHDp0qHmypIeHBzt16sTvvvvOKivz3t5R0Zpvzlu3bqWrqytr1aplcZ+R0Whk//79CYB9+vQpUrPq7T1QOgCYmed5FwBf3fGaeABl8jw/AcDvLrV6A9gFYFdISIhuJ1gIPeTk5HD69OnmpVe6d+/OCxcuaK539epVDhkyhG5ubjQYDOzZs6dF/TV3YzQa+dtvv7FXr17mKyODwcBGjRpxzJgx3Ldvn8UhkJ6ezpCQENauXVunVt/bsmXLqJRip06dLG63yWTiO++8QwB89dVXC+0E0oclgSJXKMKO3Lhxg0OGDKGTkxM9PT0tHhp84cIFDho0iM7OznRycuJrr71mlX3TjUYjt2/fzuHDh7NatWoEQAAsXbo0O3bsyAkTJnDbtm35/uvfZDIxLi6OLVq0IABu2LBB9zbfzZgxYwiAo0aNsriWyWTihx9+SADs0qVLkQgVew8UueUlHknHjh0zDw0OCQnRvFPkbWfPnmWfPn3o6OhIZ2dnDhgwgOfPn9exxf+WlJTEb7/9lp07dzZvfwyASimGhoaydevWHDhwIEeMGMHJkydz/vz5nD59Oj/++GP279+f0dHRBEB3d/cC3Y7XZDKxc+fOVEpZvNzNbaNGjSIAvvjii3a/A6S9B4ojgJMAyuXplI++4zX97+iU/+FBdSVQhL3YuHGjuX+ldu3a3Lp1q0X1Tp06xZ49e9LR0ZEuLi4cOHCgVYPltqSkJP70008cMWIEO3fuzOrVq5vnu9z5KFasGOvXr88ZM2ZYtPGYVqmpqaxWrRp9fX1163/69NNPzSsm2POVil0Hyq32ow2Ao7m3sobnfu5jAG1zP3YFsBi3hg3vABD2oJoSKMKe5OTkcM6cOQwMDCQAPvfccxa/0Z08edIcLM7OzuzXr5/Fw2a1yMrK4uXLl3nkyBGeO3dOt3W7LHXy5EkWK1aMtWrV0m102e3baV27drXbjnq7DxRrPCRQhD36+++/OWLECHp4eNDR0ZGDBg2yeLXbU6dOsU+fPnRycqKjoyNfffVV2e421+LFiwmAH3zwgW41P/74YwJgr1697DJUJFAkUEQRc+HCBfbu3ZsODg709vbm6NGjLR7qeu7cOQ4aNIiurq5USrFDhw7cuXOnTi22X127dqWDgwN36LTlsclkMg+DHjx4sN3NU5FAkUARRVRCQoJ5ZnZQUBBnzZpl8f35ixcv8t1336WPjw8BsGnTptp3eiwCrl+/zsDAQFauXJkZGRm61DSZTHz99dfNM+rtiQSKBIoo4mJjY1m7dm0CYKVKlfjzzz9bHADJycn87LPPzP021apV4/z58602W70wW7FiBQHwk08+0a2m0Whkt27dCMCu1v6SQJFAEY8Ak8nEH374gRUrViQA1q1b1+Kl2clbOz3OmjWLUVFR5tWNP//8c964cUOHVtuPZ555hl5eXhbvQ5NXdnY227ZtS6UUf/jhB93qWpMEigSKeIRkZ2fzm2++YVBQkHmXx127dllc12g0csWKFWzcuDEB0MvLi4MHD35kdis8fPgwDQYDBw0apGvdtLQ0PvHEE3RycuLGjRt1rW0NEigSKOIRlJaWxnHjxrFEiRLmocYJCQm61N61axc7d+5MR0dHKqX49NNPc82aNXY5aulhvPLKK3Rzc+PVq1d1rXv9+nVWqlSJvr6+/PPPP3WtrTcJFAkU8QhLTk7mRx99RC8vLzo4OLBLly66XVUkJibygw8+YEBAAAEwIiKCkyZNKrK3w+Lj463WkX7y5En6+/szLCzM4qHg1iSBIoEiBK9cucK33377X4tFnjlzRpfaGRkZnD9/PuvUqWNeLqVXr17cu3evLvULk2bNmjE8PNwqo962b99OFxcXNmzYsNBM8LyTBIoEihBmSUlJHDBgwL8Wi9Rz6ZVdu3axe/fudHNzMw8OmDNnjq57wtvSzJkzCUCXfqm7+e677wiAPXv2LJRDtSVQJFCE+H/OnDnD3r17/2tNr8TERN3qX7t2jRMnTmRERAQB0NPTkz169OC2bdsK5Rtlfl25coUAOGbMGKsd4/bExxkzZljtGFpJoEigCHFPJ0+eZI8ePWgwGOjq6srXX3+dSUlJutU3mUzcsmULX331VXp4eJj7Wj799NMCWZTSGkJCQti5c2er1c/JyWHLli3p7Oys2wx9vUigSKAI8UAnTpzgq6++SoPBQBcXF/bv359nz57V9Rg3b97krFmz+MQTTxAAHRwc2KxZM86ePduuOvKbNWvGOnXqWPUYV69eZUhICENDQ62yC6ZWlgSKA4QQj4SwsDDMnj0bR44cQZcuXTB9+nSUL18evXv3xsmTJ3U5hpeXF7p3744tW7bg6NGjGD58OE6fPo3u3bujZMmSeO6557B48WKkpaXpcjxrycrKgqurq1WPUaJECfzwww84f/48unfvfusvfHunNYkK+0OuUIS4vzNnzvC1116js7MzDQYDX375Zd3mseR1e2fGgQMHsmTJkuZRYi+88AL/97//2WQ/lPsxmUwsV64c27VrVyDHmzBhAgFw+vTpBXK8B4Hc8pJAEUKrxMREvvnmm3R3d6dSis8++6zVRjjl5ORww4YN7NOnj3lui7OzM9u0acNp06YVij6XTZs2EQDnzp1bIMczGo1s3rw5PTw8eOLEiQI55v1YEijq1vcXPTExMdy1a5etmyGE3bh69SomTZqEr776Cjdu3ECLFi0wbNgwNG7cGEop3Y9nNBqxfft2/PTTT/jpp59w6tQpAED16tXRunVrtGjRAvXr14eLi4vux76XtLQ0tGzZEgkJCUhKSoKbm1uBHPfcuXOoUqUKqlWrhk2bNsHBwXa9EUqp3SRjNH2z1iQq7A+5QhFCm+TkZI4ZM8Z8e6p27dpcunSpVZddMZlMjI+P55gxY9igQQM6OjoSAN3c3NiiRQuOHj2a27dvt+pKyKmpqWzSpAkdHBy4cOFCqx3nXmbPnk0AnDp1aoEfOy/IFcr/J1coQlgmPT0dc+fOxbhx43Dy5ElERETgrbfeQpcuXazeYZ2SkoLY2FisW7cOGzduRHx8PADAw8MDNWrUQK1atRATE4Nq1aohPDwcTk5Omo91/fp1LFq0CF9//TUSEhIwd+5cvPzyy3r9KPlGEi1btkRcXBwOHz6MoKCgAm8DYNkVigSKEOK+cnJysGTJEnz22WfYs2cPAgICMGDAAPTr1w9+fn4F0oYrV65g8+bNiI2Nxc6dO7Fv3z5kZGQAAJycnBAREYGoqCiUK1cOoaGhCA0NRYkSJVCsWDEUK1YMjo6OyMrKQlZWFq5fv44jR47g8OHD2LNnD1avXo3MzExUqVIFI0aMQPv27QvkZ7qbkydPIioqCi+++CLmzJljkzZIoNyFBIoQ+iKJ2NhYjBs3DqtXr4abmxu6deuGwYMHIyIiokDbkp2djYSEBBw8eBAJCQlISEjA4cOHcfbsWWRlZeW7TmhoKNq2bYtXXnkFjz32mFX6ih7Wf//7X4wbNw67d+9G9erVC/z4Eih3IYEihPUcOnQIEyZMwIIFC5CZmYmnnnoKb775Jpo0aWLTN2WTyYSLFy/izJkz+Ouvv3Djxg1cv34dRqMRzs7OcHZ2hqenJyIiIlCxYkV4eHjYrK33kpycjLCwMNSrVw8rV64s8ONLoNyFBIoQ1nf58mVMnToVU6dOxZUrV1C1alUMHjwYnTp1KrARUkXRmDFjMGzYMMTFxaFOnToFemxLAkVmygshNAsICMBHH32Es2fPYtasWTCZTOjevTuCg4MxbNgwnD171tZNtEsDBgxAiRIlMGbMGFs35aFIoAghLObq6oru3bvjwIED2LhxIxo2bIjPPvsM5cqVQ/v27bF+/XoU1bsh1uDp6Ym+ffti2bJlOHHihK2bk28SKEII3Sil0KRJEyxduhSnTp3C0KFDsXXrVrRo0QJRUVGYNGkSbty4Yetm2oX+/fvDYDBg5syZtm5KvkmgCCGsIiQkBJ9++inOnTuHefPmwdfXF4MHD0ZgYCB69OiBXbt2yVXLfZQuXRpNmzbFjz/+aDfnSQJFCGFVrq6u6NKlC+Li4rBnzx68/PLLWLRoEWrVqoUaNWpg6tSpctVyD8899xyOHz+OQ4cO2bop+WLTQFFKFVdKrVNKHcv9t9g9XmdUSu3LfSwv6HYKIfRRvXp1zJgxA0lJSZg6dSqUUujfvz8CAwPRtWtXxMbG2s1f4wWhadOmAIBt27bZuCX5Y+srlHcAbCAZDmBD7vO7SSf5WO6jbcE1TwhhDT4+PujXrx/27NmD3bt3o1u3bli2bBkaN26MihUrYvTo0Th37pytm2lz5cuXh7+/P/744w9bNyVfbB0o7QDMzf14LoBnbNgWIYQN1KhRA19//TUuXLiAb7/9FoGBgRg+fDjKli2LFi1aYMGCBUhNTbV1M21CKYUKFSrYzfBrWwdKSZIXcj++CKDkPV7nqpTapZSKU0pJ6AhRBLm7u6Nbt26IjY3FiRMn8OGHH+LEiRPo0qULSpYsia5du2Lt2rXIycmxdVMLVMmSJXHx4kVbNyNfrB4oSqn1Sqn4uzza5X1d7rLJ97p5WjZ35uaLAL5QSpW/x7F65wbPritXruj7gwghCkxYWBg+/PBDHD9+HJs3b8aLL76I5cuXo1WrVggKCsKgQYOwffv2R6K/xWQy2XR/lIdh06VXlFJHADQmeUEpVRrAbyTvu8qcUupbACtJ/ni/18nSK0IULRkZGVi9ejX+97//YcWKFcjMzETZsmXx/PPP4/nnn0etWrUKxeKOeqtfvz7c3d2xfv36AjmePS+9shxAt9yPuwFYducLlFLFlFIuuR/7AXgcgH2MoRNC6MbV1RXPPvssFi9ejMuXL2Pu3LmIjo7GpEmTUKdOHYSGhuKNN97A1q1bYTQabd1cXaSnp+PAgQOoWLGirZuSL7YOlDEAWiiljgFonvscSqkYpdTt6aFRAHYppfYD2ARgDEkJFCEeYd7e3ujatStWrVqFS5cu4dtvv0W1atXw9ddfo0GDBggKCkKvXr2wYsUKpKWl2bq5mm3YsAGpqalo167dg19cCMhqw0KIIiMlJQWrV6/G0qVL8csvvyAlJQWurq5o1qwZnnzySbRp0wblypWzdTPzxWg0olGjRjhy5AgSExPh7OxcIMe15JaXo96NEUIIW/Hy8kLHjh3RsWNHZGVlITY2FitWrMCqVauwatUqAEBERARatmyJli1bolGjRvDy8rJxq+/uvffew7Zt2zB37twCCxNLyRWKEKLII4ljx47hl19+wa+//orY2Fikp6fD0dERMTExaNSoERo2bIjHH38cPj4+Nm1rRkYGBg8ejOnTp6NXr16YMWNGgR5fNti6CwkUIcS9ZGZmYtu2bVi/fr15n/rs7GwopVClShU8/vjjqF+/PmrWrImKFSvCYDAUSLt2796NPn36YPfu3RgyZAhGjRoFJyenAjn2bRIodyGBIoTIr7S0NGzfvh1bt27Ftm3bEBcXh5SUFACAh4cHqlevjqpVq6Jy5cqIjo5GVFQU/Pz8dBmmfP78eSxbtgxz587Fzp07Ubx4ccyePdtmHfESKHchgSKE0MpoNOLQoUPYvXs3du/ejT179iA+Ph43b940v8bLywsVKlRAWFgYgoKCEBgYiMDADGOBvAAABoRJREFUQBQrVgze3t7w9vaGs7MzTCYTjEYjMjMzcfnyZVy8eBFJSUnYu3cvduzYgaSkJABA1apV0a1bN/To0cOmt90kUO5CAkUIoSeSOH/+PBISEnD06FEcP34cx48fx8mTJ5GUlGS+osmv8PBw1K5dG3Xq1EGjRo1QtWpVK7X84cgoLyGEsDKlFIKDgxEcHIzWrVv/v6+npKTgwoULSE5Oxs2bN5GSkoLMzEwYDAY4ODjAyckJJUuWRKlSpVCyZEm4uLjY4KewLgkUIYTQgZeXV6EdglxQbD1TXgghRBEhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXEihCCCF0IYEihBBCFxIoQgghdCGBIoQQQhcSKEIIIXQhgSKEEEIXNg0UpdTzSqkEpZRJKRVzn9e1VkodUUodV0q9U5BtFEIIkT+2vkKJB/AsgM33eoFSygBgCoAnAVQC0FkpValgmieEECK/HG15cJJ/AoBS6n4vqw3gOMmTua9dBKAdgENWb6AQQoh8s2mg5FMQgHN5np8HUOduL1RK9QbQO/dpplIq3sptsxd+AK7auhGFhJyLf8i5+Ieci39EaP1GqweKUmo9gFJ3+dJwksv0PBbJGQBm5B53F8l79ss8SuRc/EPOxT/kXPxDzsU/lFK7tH6v1QOFZHMLSyQCCM7zvEzu54QQQhQitu6Uz4+dAMKVUuWUUs4AOgFYbuM2CSGEuIOthw23V0qdB1APwCql1JrczwcqpVYDAMkcAAMArAHwJ4AfSCbko/wMKzXbHsm5+Ieci3/IufiHnIt/aD4XiqSeDRFCCPGIsodbXkIIIeyABIoQQghd2H2gPGhZFqWUi1Lq+9yv/6GUCi34VhaMfJyLhkqpPUqpHKVUB1u0saDk41y8qZQ6pJQ6oJTaoJQqa4t2FoR8nIu+SqmDSql9SqmtRXklivwu46SUek4pxfstCWXv8vF78YpS6kru78U+pVTPBxYlabcPAAYAJwCEAXAGsB9ApTte8xqAabkfdwLwva3bbcNzEQqgKoB5ADrYus02PhdNALjnftzvEf+98M7zcVsAv9q63bY6F7mv88Kt5aDiAMTYut02/L14BcBXD1PX3q9QzMuykMwCcHtZlrzaAZib+/GPAJqpB6z1YqceeC5IniZ5AIDJFg0sQPk5F5tIpuU+jcOt+U1FUX7Oxc08Tz0AFNWROvl5vwCAkQDGAsgoyMYVsPyei4di74Fyt2VZgu71Gt4agpwMoESBtK5g5edcPCoe9lz0APCLVVtkO/k6F0qp/kqpEwA+AzCogNpW0B54LpRSNQAEk1xVkA2zgfz+P/Jc7m3hH5VSwXf5+r/Ye6AIYRGl1MsAYgCMs3VbbInkFJLlAfwXwHu2bo8tKKUcAEwA8Jat21JIrAAQSrIqgHX4507PPdl7oORnWRbza5RSjgB8APxVIK0rWLJEzT/ydS6UUs0BDAfQlmRmAbWtoD3s78UiAM9YtUW286Bz4QWgMoDflFKnAdQFsLyIdsw/8PeC5F95/r+YCaDmg4rae6DkZ1mW5QC65X7cAcBG5vY4FTGyRM0/HngulFLVAUzHrTC5bIM2FpT8nIvwPE+fAnCsANtXkO57Lkgmk/QjGUoyFLf61tqS1LxYYiGWn9+L0nmetsWtlUruz9ajDXQYrdAGwFHcGrEwPPdzH+PWLwIA/F97d6waRRSFAfg/YJFefAELCQTLNJa+gsFS8iaW1hY+RSCFpQiKhUiqLJIuLxCIndbXYqZKo8Vhrm6+Dxa2GJbD4cLPvTM75yDJWZLrJBdJHs+ueWIvjrOclf7Ksku7ml3zxF58THKT5HL9vJ9d88RevE1ytfbhU5Kj2TXP6sWdaz9nT5/y+st18WZdF7t1XRz+6Te9egWAFv/7kRcA/wiBAkALgQJAC4ECQAuBAkALgQJAC4ECQAuBAhuqqpOq+lZVu3X2yKPZNUEXf2yEDVXVwzHGj/X76yS3Y4x3k8uCFnYosK3Tqrqoql2W4W/7PHODe+bB7ALgvqiqV1kGGz0fY/ysqi9Z3pUEe8EOBbbzNMnXNUxeJHmW5PvkmqCNeyiwkao6SnKeZWrohyQvxxhP5lYFfQQKAC0ceQHQQqAA0EKgANBCoADQQqAA0EKgANBCoADQ4jd8ivK/928x6wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Determine chi-squared across the 2D parameter grid (code is the same as Class 3)\n", "amin,amax,na = 0.,0.5,100\n", "bmin,bmax,nb = -1.,2.,100\n", "da,db = (amax-amin)/na,(bmax-bmin)/nb\n", "alst = np.linspace(amin+0.5*da,amax-0.5*da,na)\n", "blst = np.linspace(bmin+0.5*db,bmax-0.5*db,nb)\n", "chisq2 = np.empty((na,nb))\n", "for i in range(na):\n", " for j in range(nb):\n", " chisq2[i,j] = chisqstraightline2([alst[i],blst[j]],x,y,yerr)\n", "minchisq = np.amin(chisq2)\n", "i,j = np.unravel_index(chisq2.argmin(),chisq2.shape)\n", "a,b = alst[i],blst[j]\n", "alst2,blst2 = np.meshgrid(alst,blst,indexing='ij')\n", "\n", "# Plot the best-fitting values of the parameters from each jack-knife sample\n", "# Also plot the (1-sigma, 2-sigma, 3-sigma) contour values of chi-squared (code is the same as Class 3)\n", "clst = np.array([0.,2.30,6.17,11.8])\n", "clst += minchisq\n", "plt.contour(alst2,blst2,chisq2,levels=clst,colors='k')\n", "plt.plot([a],[b],marker='o',markersize=5,color='black',linestyle='None')\n", "plt.scatter(asamp,bsamp,marker='o',s=20,color='red')\n", "plt.xlim(amin,amax)\n", "plt.ylim(bmin,bmax)\n", "plt.xlabel(r'$a$')\n", "plt.ylabel(r'$b$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can compare the estimated jack-knife errors with the errors obtained from the $\\chi^2$ contours:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jack-knife error in a = 0.05193629364182536\n", "Jack-knife error in b = 0.3207081549826602\n", "Chi-squared error in a = 0.05292881652401213\n", "Chi-squared error in b = 0.33002843260768994\n" ] } ], "source": [ "from scipy import stats\n", "\n", "# Jack-knife error is the standard deviation of the best-fitting parameters to the samples, scaled by \\sqrt(nsamp-1)\n", "norm = np.sqrt(float(nsamp-1))\n", "s = stats.describe(asamp)\n", "print('Jack-knife error in a =',norm*np.sqrt(s.variance))\n", "s = stats.describe(bsamp)\n", "print('Jack-knife error in b =',norm*np.sqrt(s.variance))\n", "\n", "# For comparison, we determine the error in the parameters using the chi-squared method (code is the same as Class 3)\n", "# The different determinations of the error agree well\n", "chisqa = np.amin(chisq2,axis=0)\n", "chisqb = np.amin(chisq2,axis=1)\n", "i = chisqa.argmin()\n", "abot = np.interp(minchisq+1.,np.flip(chisqa[:i]),np.flip(alst[:i]))\n", "atop = np.interp(minchisq+1.,chisqa[i:],alst[i:])\n", "i = chisqb.argmin()\n", "bbot = np.interp(minchisq+1.,np.flip(chisqb[:i]),np.flip(blst[:i]))\n", "btop = np.interp(minchisq+1.,chisqb[i:],blst[i:])\n", "asig = 0.5*(atop-abot)\n", "bsig = 0.5*(btop-bbot)\n", "print('Chi-squared error in a =',asig)\n", "print('Chi-squared error in b =',bsig)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In some situations, the jack-knife samples could be created by deleting **regions** or **groups of points**, not individual points.\n", "\n", "## Bootstrap errors\n", "\n", "The **bootstrap procedure** is another method for estimating errors by re-sampling the data. If we have $N$ data points, the procedure is:\n", "\n", "- **Repeatedly draw samples of $N$ points at random, with replacement**\n", "- Hence, the same point can appear multiple times in each bootstrap sample\n", "- Measure the statistic for each of these bootstrap datasets (which can number $N_{samp} \\gg N$), creating measurements $S_i$\n", "- The bootstrap error is the standard deviation of $S_i$\n", "\n", "Here's a bootstrap analysis of the straight line dataset:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, '$b$')" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeXgT1f7G30nSfWVpEQqlUJYWKAiUUkpbNgVFRXEFdwFRlE1Er6h4LxdFEAX5IagooKBcRFRUQGUVkEVk32Rfu7AJFOiaZN7fH5NJJ2mapmm6wfk8zzwlycyZk+ThvPme7yaRhEAgEAgEZUVX2RMQCAQCwY2BEBSBQCAQeAQhKAKBQCDwCEJQBAKBQOARhKAIBAKBwCMIQREIBAKBR6h0QZEkqYEkSWslSTogSdJ+SZJGODhHkiTp/yRJOipJ0h5JktpVxlwFAoFAUDyGyp4AABOAl0nukCQpCMB2SZJWkjygOedOAE0tR0cAH1v+CgQCgaCKUOkWCslMkjss/74G4G8AEXan3QtgHhW2AAiVJKluBU9VIBAIBE6oChaKFUmSogC0BfCn3UsRAM5oHqdZnsu0u34wgMEAEBAQ0D4mJqa8pioQCAQ3JNu3b79IMsyda6uMoEiSFAjgOwAjSV51ZwySswDMAoD4+Hhu27bNgzMUCASCGx9Jkk65e22lb3kBgCRJXlDE5GuS3zs4JR1AA83j+pbnBAKBQFBFqHRBkSRJAjAbwN8kpxRz2k8AnrREeyUCyCKZWcy5AoFAIKgEqsKWV2cATwDYK0nSLstzrwOIBACSnwBYDqA3gKMAcgA8UwnzFAgEAoETKl1QSP4BQCrhHAJ4sWJmJBAIBAJ3qPQtL4FAIBDcGAhBEQgEAoFHEIIiEAgEAo8gBEUgEAgEHkEIikAgEAg8ghAUgUAgEHgEISgCgUAg8AhCUAQCgUDgEYSgCAQCgcAjCEERCAQCgUcQgiIQCAQCjyAERSAQCAQeQQiKQCAQCDyCEBSBQCAQeAQhKAKBQCDwCEJQBAKBQOARhKAIBAKBwCMIQREIBAKBRxCCIhAIBAKPIARFIBAIBB5BCIpAIBAIPIIQFIFAIBB4BCEoAoFAIPAIQlAEAoFA4BGEoAgEAoHAIwhBEQgEAoFHEIIiEAgEAo9gqOwJAIAkSXMA3A3gPMlWDl7vCuBHACcsT31P8r8VN8Pqz/Xr13H58mUUFBSgoKAAZrMZfn5+CAgIQGBgIAICAiBJUmVPUyAQVGOqhKAA+ALARwDmOTlnA8m7K2Y61Ze8vDz89ddf2LhxIzZv3ozjx4/jzJkzyMrKcnqdt7c3brnlFtxyyy2oV68eGjdubD1iY2MRGRkJnU4YtAKBoHiqhKCQXC9JUlRlz6M6s27dOkybNg3Lli1DQUEBACAmJgYxMTHo2rUrGjRogJo1a8LHxwdeXl7Q6XTIzc1FdnY2srOzcfHiRZw9exZnz57FoUOH8OuvvyIvL886vr+/P2JiYtC6dWu0a9cO8fHxaNOmDfz9/SvrLQsEgipGlRAUF+kkSdJuABkARpPcX9kTqgrs3r0br7zyClauXImwsDAMGTIEPXr0QFJSEmrVquX2uCRx9uxZHDt2DH///TcOHDiA/fv3Y9myZfjiiy8AAAaDAW3atEFSUhKSkpKQnJyM+vXre+idCQSC6oZEsrLnAACwWChLi/GhBAOQSV6XJKk3gGkkmzo4bzCAwQAQGRnZ/tSpU+U76Urms88+w9ChQxEcHIwxY8ZgyJAh8PPzK9d7kkR6ejq2b9+OrVu3YvPmzfjzzz+Rk5MDAGjUqBFSUlLQrVs33HbbbUJgBIJqhiRJ20nGu3VtdRAUB+eeBBBP8mJx58THx3Pbtm0em19VY8GCBXjsscfQq1cvfP3112WyRsqKyWTCnj17sGHDBqxfvx4bNmzAhQsXACjbbj179sQdd9yBrl272gqeLAMXLgDh4YAICBAIqgQ3vKBIknQLgHMkKUlSAoDFABrSyeRvZEHZtWsXkpKS0KFDB6xatQpeXl6VPSUbSGLv3r1YuXIlVq5ciXXr1iEvLw++vr7o0qUL7r77btzduzeinnkG2LQJSEoC1q4FhNNfIKh0qr2gSJL0PwBdAdQGcA7AvwF4AQDJTyRJGgpgCAATgFwAo0hucjbmjSooJNGmTRtcunQJ27dvR506dSp7SiWSm5uL9evX45dffsHy5ctx5MgRAEAcgHsB9NXr0TYtDdItt1TqPAUCwQ0gKOXBjSoof/31FxISEvDZZ59h0KBBlT0dtzh8+DB+/ukn/PT22/gjKwsygIYNG+L+++/Hgw8+iMTERBGiLBBUEmURFPG/tprxww8/wGAw4IEHHqjsqbhNs2bN8PLo0Vh36RLO7tuH2Z9/jlatWmHGjBno3LkzGjZsiJdeeglbtmzBjfqDRyC4EREWSjXjwQcfxP79+/H3339X9lQ8TlZWFn7++WcsXrwYv/76K/Lz8xEVFYV+/frh0UcfRVxcXGVPUSC44REWyk1EWloaIiIiKnsa5UJISAgef/xxLFmyBOfOncOXX36JmJgYTJ48Ga1bt0abNm0wefJkpKenV/ZUBQKBA4SgVDNIwmCoTvmo7hESEoInn3wSv/zyCzIzMzFjxgz4+/vj1VdfRYMGDdCzZ098/fXXyM7OruypCgQCC0JQqhk1a9bEpUuXKnsaFUpYWBheeOEFbN68GUeOHMHYsWNx5MgRPP7447jlllswaNAgbN68WfhbBIJKRghKNaN+/fo4duzYTbt4NmnSBOPGjcOxY8ewbt06PPTQQ1i4cCGSkpLQsmVLfPDBB9akSoFAULEIQSkJWQbOnQOqyAIeHx+PS5cu4cSJEyWffAOj0+mQmpqKOXPmIDMzE59//jlCQkIwevRoRERE4OGHH8aKFSsgy3JlT1UguGkQguIMWQa6dQPq1we6dlUeVzIpKSkAgBUrVlTyTKoOQUFBGDhwIDZv3ox9+/Zh6NChWLNmDXr16oVmzZph8uTJuHix2Co9AoHAQwhBccaFC0ppEJNJ+VsFtlJiY2PRuHFjLFmypLKnUiVp2bIlpkyZgvT0dCxYsAARERF49dVXERERgSeeeEL4WgSCckQIijPCw5U6UwaD8jc8vLJnBEmS8PDDD2PVqlXIyMio7OlUWXx8fNC/f3+sW7cO+/btw+DBg/Hjjz8iKSkJ7du3x+zZs60VkgUCgWcQguIMSVKKFqalAb//XmUq4g4aNAiyLOPTTz+t7KlUC1q2bInp06cjPT0dH3/8MYxGIwYNGoQGDRrg1Vdfven9UQKBpxCCUhI6HVCnTqGYVAEnfXR0NHr37o1PPvkEubm5lTaP6kZQUBCef/557NmzB7///ju6d++OKVOmIDo6Gn379sX69evFdphAUAaEoJSGKuSkf/XVV3H+/HnMmjWr0uZQXZEkCV26dMG3336LkydPYsyYMVi/fj26dOmC+Ph4fPXVV9Y2ygKBwHVELa/ScO6cIiYmk+JXSUtTrJdKolu3bjh06BCOHDmCgICASpuHM7Kzs3Hq1Cnk5OQgOzsbeXl58PPzQ1BQEIKCghAeHo7g4ODKniZycnLw1VdfYerUqTh48CDq1auHYcOG4bnnnkONGjUqe3oCQYUhytc7oFwEhVQsE7UpVCX7VTZu3Ijk5GSMGzcOb731VqXNQ8u5c+ewYMECbNy4EXv27MHRo0dL3EYKDg5GZGQkoqKiEBMTg5iYGLRo0QJxcXEIDAysoJkryLKMX3/9FVOnTsWqVasQEBCAgQMHYuTIkWjUqFGFzkUgqAyEoDig3KoNV7G2tQ899BCWL1+OgwcPokGDBpU2j9OnT+OVV17Bd999B7PZjOjoaLRp0wZt2rRB06ZNERQUBH9/f/j6+iI3NxfXrl3D1atXce7cOZw5cwanT5/GiRMncOjQIeTn5wNQtqZiY2MRHx+PxMREpKamokWLFpAq6HPfvXs3pkyZggULFkCWZTzwwAN49dVXER/v1v81gaBaIATFATdq+Xp7Tpw4gZYtW+L222/HkiVLKmyx1bJixQrcf//9IIkXX3wRAwYMQExMjFtjmc1mnDx5Evv378fOnTuxbds2/PXXXzh37hwAoHbt2khNTUX37t3Ro0cPNG/evNzfc3p6Ov7v//4Pn3zyCa5evYoePXrgtddeQ48ePSrl8xYIypOyCApI3pBH+/btebPw3nvvEQAXL15c4fdev349fX192bp1a548ebJc7iHLMo8ePco5c+bw6aefZsOGDQmAAFivXj0OGDCA33//Pa9du1Yu91e5cuUKJ02axLp16xIA27dvz2+//ZZms7lc7ysQVCQAttHNdVdYKDcAJpMJCQkJSE9Px/79+1G7du0KuS9JxMfH4/Lly9i6dWuF3vfEiRNYvXo1Vq5ciRUrViArKwve3t7o3r077rvvPvTp0wd169Ytl/vn5+dj3rx5mDx5Mo4cOYKYmBi89tprePTRR+Hl5VUu9xQIKgphodzkFgpJ7t69m15eXnz44Ycr7J5r1qwhAM6aNavC7umIgoICrl27lqNGjWJ0dLTVeklMTOQHH3zA06dPl8t9TSYTv/nmG7Zp04YA2LBhQ86YMYO5ubnlcj+BoCJAGSyUSl/4y+u42QSFJN9++20C4BdffFEh9xs+fDgDAgKq1AIqyzL37dvHt99+m23btrWKS1JSEqdPn85z586Vyz2XLl3KTp06EQBvueUWvv/++7x+/brH7yUQlDdCUISgkFR+MXft2pX+/v7cv39/ud/vkUceYbNmzcr9PmXh8OHDfOeddxgXF0cA1Ov17NWrF+fNm+dxn4ssy1yzZg27detGAKxduzYnTJjAq1evevQ+AkF5IgRFCIqV9PR0hoWFMTY2ttyd1HfccQer0+e8d+9evv7664yKiiIABgYG8plnnuH69espy7JH7/XHH3/wjjvuIADWrFlTCIug2iAERQiKDStXrqQkSXzsscc8vlBqee6551irVq1yG7+8kGWZ69ev54ABAxgYGEgAbNq0KSdMmMC0tDSP3uuvv/7i3XffLYRFUG0QgiIEpQjjx48nAE6bNq3c7jF16lQCYHp6erndo7y5fv06v/jiC3bp0oUAqNPpeM8993DZsmU0mUweu8/WrVvZu3dv61bYe++9J3wsgiqJEBQhKEUwm8289957qdPp+Ntvv5XLPXbt2lUlorw8xdGjRzlmzBjWqVPHGrU1YcIEjzryt2zZwl69ehEA69Spww8//LBKBTUIBEJQhKA45Nq1a2zdujVDQkJ44MABj48vyzIbN27M22+/3eNjVyb5+flctGiR1bnu5eXFxx57jJs3b/bYFuIff/xhHT8iIoIff/wx8/PzPTK2QFAWhKAIQSmWkydPsk6dOmzcuDHPnz/v8fH//e9/EwAPHz7s8bGrAgcOHOCwYcMYFBREAIyPj+eXX37JvLw8j4y/Zs0aJiUlEQAbNWrEefPmeXSrTSAoLdVeUADMAXAewL5iXpcA/B+AowD2AGhX0phCUArZsmULfX19mZSUxJycHI+OnZmZSS8vLw4aNMij41Y1rl69yhkzZjA2NpYAGBYWxrFjx/Ls2bNlHluWZS5btsyaN9OiRQt+//335RpQIRAUx40gKKkA2jkRlN4AfrEISyKAP0saUwiKLd9++y0lSWKfPn1oNBo9OvaQIUPo7e3NEydOeHTcSsdsJs+eJTULuyzLXLlyJe+55x5KkkQfHx8OGjTII1uKZrOZ3377LWNiYgiACQkJXL16dZnHFQhKQ7UXFOU9IMqJoHwKoL/m8SEAdZ2NJwSlKNOnTycADho0yKO/fs+cOUM/Pz/279/fY2NWOmYzmZpKGgzKXwcFIA8dOsTnn3+evr6+BMC7776b69atK/NnazQaOWfOHDZo0IAAeNttt3H79u1lGlMgcJWbQVCWAkjWPF4NIN7BeYMBbAOwLTIy0kMf743FG2+8QQB87bXXPDrum2++SQA3zi/qs2cVMQGUv062ts6fP89x48axdu3aBMCOHTvy+++/L3MV4tzcXE6ZMoW1atUiAPbv35/Hjx8v05gCQUkIQREWisvIssznn3+eADhp0iSPjZuTk8MmTZowOjqa2dnZHhu30pBlWwvFBasjOzubM2fOZOPGjQmAMTExnD17dpkd+FeuXOHrr79OPz8/enl5ccSIEbx48WKZxhQIiuNmEBSx5eVBTCYT+/XrRwCcMWOGx8ZVqw+PHDnSY2NWKg58KK5gNBq5cOFC3nrrrQTA+vXr8//efps5ZRTatLQ0Dho0iDqdjiEhIZw0aZLIYRF4nJtBUO6yc8pvLWk8ISjOKSgoYJ8+fQiAc+bM8di4Q4cOJQCuWrXKY2NWKG6KiKMxZLOZvy5fzpTgYAJguJcXJ02cWOYaa/v27eNdd91FAIyMjORXX30lmnwJPEa1FxQA/wOQCcAIIA3AQADPA3je8roEYAaAYwD2lrTdRSEoLpGbm8uePXtSkiTOmzev7AOazcw+fpzNmzdn3bp1mZmZWfYxPYGrIuGCI96le6WkkHq98jcjgzQYuA5gT0kiANaqVYvvvvtumWt6rV69mu3atbPmx6xbt65M4wkE5A0gKOVxCEFxjezsbPbo0aPsoqJZjHe3a0c/Pz+mpqayoKDAc5N1kYKCAv7www8cOHAgu3fvzsa+vvQC6KvTsXbt2oyKimKHDh3Yv39/jh07lvPmzeOBAwdotiz+JTrinQlURoZyvXqkp9uI1JZly3jnnXd6TFjMZjPnzZvH+vXrEwD79u3LI0eOuD2eQCAE5UYWFE9swZSAVlTmzp3r3iB2UVFfzZhBABwxYkTx76Ec3ltmZiZbtWplre6b2L49+0kSXwH4iiTxhaef5hNPPMGePXuycePG1Ol01iZcwcHBvC00lON1Om659VaaHOXrOLJizGZFSDIzlb+SpHwOkqQISkaG8jclpVBYNm2yCosnikVmZ2dz/PjxDAgIoJeXF0eNGsXLly+X4ZMU3KwIQblRBcWdLRg3F+ns7GzefvvtlCSJn332WennmZlps2BSljly5EgC4OxmzYq+B+17U7eGyigs6enpbN68OQMCAvjNN98oCZwlRGvl5+dz3759nDt3LocMGWJt5wuAoaGhfPDBBzl37tzCjHj7cOKMDGVcVUSSksiOHQvfl3rvTp2UbTA762fLli3WvimeKBaZkZHBgQMHUpIk1q5dmzNnzvR4IqvgxkYIyo0qKKXIhSBZZh9Abm6u9VfzRx995NpF9vfMyCBNJvLsWRoLCnh7ly70ArhOfQ+qX0X73tTX3PVbWEhOTmZgYCA3bNhQdI6lENnz58/zf//7H5955hnWrVvXKjAdO3bklA8+YLoqGKmpyvvRvg/VMmnXjkxLs/3+EhOLFbY//viD3bt3txaL/PTTT8u0Xbhz505rSf5WrVpx5cqVbo8luLkQgnKjCkppcyFKK0AOyMvLs0Z/TZw4seQLivvFbpnzpQsXGOPnx1CA+wDlV7vZXPje9PrCX/faORuN5N69pRKYgIAADh8+vNTv2RmyLHPHjh0cP3681QEuSRK7JiXxs1mzePniRbJ9e1tBUY/ERNvvzyK0zr7H1atXW3vTN2nShAsWLHA7gkuWZS5evJiNGjUiAN533308evSoux+F4CZBCMqNKihk6X5du5GM54iCggJrnsobb7zhvJSI/T21v9gtAnFy61bWBVgf4Gm9vlA01K0y+zkbjWRIiDJGSIjyuATy8vIIgGPHjnXrPbvKoUOH+J9//5tNLcmLPpLEBwH+7OdHo05HBgXRxurKyFDerwtioiLLMpcuXcrWrVsTAFu3bs1ly5a5XdIlNzeXEyZMYEBAAL29vTlmzJhybw8tqL4IQbmRBaW0eMjRbTKZ+OyzzxIAn3/+eecl1bX3dCRqssxd7doxGGCMnx/P2zessp/z3r20+aW/d2+J79FsNrNevXp85plnSjzXKSWdb9nik/V6bo2L43BJYphlS6x+nTr895tv8nRcXFGB7NSp1Nt6ZrOZCxYsYHR0NAGwa9eu3Lp1q2vvwwHp6el84oknrNtqCxYsEBWNBUUQgiIEpVyQZZmvvfYaAfChhx5yvYSIo0XZbOa6H36gr68v27VrxytXrji/XmuhaBdgJ36inj17skWLFoXnabffOnVSrITi5lfC2Fbst/gSE1mg1/O7Fi3Yq2dPSgB1APvUqMEVv/xC2WRStr60VksptyLz8/M5ffp0hoWFEQAfeeQRHjt2rFRjaNm0aZN1+y4lJYW7d+92eyzBjYcQFCEo5cr7779PAOzevTuzsrLKNNbSpUtpMBjYuXNn59suWh+KVgC0C7peX+jkJzlt2jQC4JZNm4r6Z1SfhtFYvGjYh/xmZBQVH3sLrKCgcJ5nz/K4Xs8xgNVqad6kCT+SJF7TzsFNqyArK4tjx46lv7+/NTT40qVLbo1lMpk4a9Ys1q5dm3q9nsOHDxdhxgKSQlCEoFQAX375JfV6Pdu2bVu2DHizmd9+9hl1Oh27du3K7OPHnS+w9laDyaQ49tVFXyMK165dY+3atXln9+5FI69UAdq7t/jAhcxM2/PtkhJpNNr6QwoKbLeyCgrIhAQSYB7AedHR7BAfTwAMATgK4PH4+DKHf6enp1tDg2vUqMGpU6e63T740qVLfOGFF6jT6RgeHs65c+eKMi43OUJQhKBUCMuXL6e/vz8bNWrEQ4cOOT7JmQ9CIw5fNW9OCWA3SeL1zp2LWiIqjiLXnGS0T5w4kQC4UvVjJCcrVoFeXxjWrM1/ycy0tT6Sk0mdjuzcmdy1yzZ3RBv2azTabmXp9baPLYeckMBNLVqwH0ADQD3A/vfdx12rVxcKk73DXh1bnXMxC/zu3bt5++23EwCbNm3KJUuWuO0T2bFjhzW6LDk5mXv27HFrHEH1RwiKEJQKY+vWrQwLC2OtWrW4adMm2xeLyyJ3tF1lMHC+JFEHMAXg1cOHHW9FybJtbSxHjn+NKOTm5rJJkyZs2rQpc0+eVJ5X/Slq4mVKipIjYp9YaTIVPhccbCsOgYG24rJnjyI86uv+/oWvaw+dzvqezwB8pW5dBur1BMA79Hqu1+sVP5EqHqoDXytU2qg4Bx0kly1bZu3y2KNHD7d9ImazmbNnz2atWrWo1+v58ssvl7nemKD6IQRFCEqFcuTIEUZHR9PX15eLFy8ufKGEnBSbBTs1lUxJ4Tc6HfUAE9q25T8OMsmLdZTbi4RqfcgyV6xYQQB85ZVXip+bdutLa4FoRcI+WbF9+0IfS0KCkhFvLzr213XooJyr0ynjp6Xxkk7HdzR+llSAKwHKer0iVNp5qT6XEgIGCgoKOH36dNasWZM6nY7PP/88L1y44Nb3e/HiRWuEX0REBL/77jsRDXYTIQRFCEqFc/78eesWyaRJk5QFx4WcFJtf2ZZ//7hkCb29vdk6IIBn1V/q6gLmxAlv85ok2WwRqU3Els6f79iq0ea6ODrsrY3AwKJik5Bg6/Qv6UhOthGhbIAfAqxnEZZOAH9r2ZKyKpKdOhUKh4tJq//88w9HjBhBvV7P0NBQTps2ze2M+82bN1tL0dx11108ceKEW+MIqhdCUISglB9OfCI5OTl8+OGHCSh96gsKCmyjs0wm2/BfJ7ksK1eupL+/P5s2bswT2ja36paXAye8jUhoLY2zZ5lz/TrbBASwJsBjHToo87IvK68VDXsBMRjItm2Ve/r4OBaI4qwZVw9/fxJgrp8fZwBsYBGWxPbtueKbbyhrrbH0dFu/SgkWw/79+63+lZYtW7rdmtloNPKDDz5gQEAA/fz8OGnSpEqpIC2oOISgCEEpH1zIyzCbzXz99dcJgN26deM/2qgn7aKt9QUUw8aNGxkaGsqIiAju37+/8AVnZeWLKUzJs2d5RK9nDYAtAF7ZtMl2DG2GviowycmFwtWhg2MR6NCBDAhwLELuHBqxygP4SVQUGzRooGyFtW/PdWvW2ApqYqJTYdYiyzKXLFliLb3y8MMP8/Tp0y5da8/p06d57733Us3c37x5s1vjCKo+QlCEoJQPpagNNm/ePHp7e7MpwIOOFm0XS8Hs2bOHdevWZc2aNblx40blSVdKythvpVlEZo1OR4MksVfPnsxPTrYdw/6atDQlsuvkyeKtj+bNyy4izo5du5h37Rqn+/ryFovF0hPgNq0lpd06dKGkS05ODseNG0dfX1/6+/vz3XffdTvM+IcffmBERAQlSeKLL75Y5rwkQdVDCIoQlPKhlLXBNqxfz9peXgwFuCIuruii7SLHjh1jkyZNbJ3+7nRd7NyZjIvj55ZOif0eeYSm9PTiM+RVf0irVu6JgU5Htmjh/JwOHci4OOfnWF7PBjgZYC2LsDwC8GhgIHnmTOH2XUhI0UoAxXDixAned999BMCYmBiuWbPG5e9ES1ZWFocPH05JkhgREcEff/zRrXEEVRMhKEJQyo9SCsKJY8cYFxtLvV7PadOmuR0ddOHCBXbq1ImSJHHKlClKCRO1iZVWqOx/oduXxbcckyyL8uDBg4vOqZhrXPKB6HRFQ4yLO7Zvdx4IYC9Oln9fAfgmQH8ouSzDJIkXHF2TmOhS0uSyZcvY2FLcsn///m4nqm7ZsoVxcXEEwAcffLDqtHwWlAkhKEJQqhRXr1617rcPHDjQ9RpgduTk5PCBBx4gAA6pW5dGQLEi1MZVjn6ha8vi2y24Y4YNIwAOHzaMsn1Co+qncHSoYwUFFR23bVty2zbn4qA60vfscV2sHAhcBsDnoNQKCwH4PhS/i8197CPpnHy2b731Fr29vRkSEsKZM2e6lSFfUFDAd955hz4+PqxRowbnzJkjQoyrOUJQhKBUOcxmM998800CYFJSktu/Xs1mM1954QUC4B0As9SF05FFof5CV3NUtJZDcjJls5kvWbpIvixJSniuWk4lPb14K0WvJ3fvLsyit883KSnaKy5OKcty+rRtmHG7duTx46TqeNfrFdFyllgZHMx9Oh3vDA0lAEYD/B6grEbA2ef6lCAShw4dYo8ePag2ENu5c6db39PBgweZnJxMALzttttEiHE1RgiKEJQqy6JFi+jv78969eq5Hxkky5zVpAkNAFvCUg9LXXDtQyfli80AACAASURBVH+1gQNGoyIEmvbCcmYmh1p8Ks9LEk0JCcp1ycm2UV87dyo5I6p1oc2pUTsyliZs2D4JUj1UgQkKUgRHFS1nW2G7d5NmM39t1YotoWzl9UhJ4b69ex0HUpRgsciyzK+++orh4eHU6/V89dVXmZ2dXeqvyWw2c8aMGQwMDGRAQACnT58u6oJVQ4SgCEGp0uzatYuNGjWit7d36fvVq5jNXL1oEUNDQli7Zk2uVxddnU5Z3NX8Eu1WlqOFVJYpp6TwNdVRDzBfXazT0gr9MvZl7+3zYVJSlIgwV8TE1eTHhQtdEyk/P8UxbzLRePo0p7/zDkNDQ5WqwcOG8bLavz4lpVQWyz///MNBgwYRABs1asTffvvNra/q1KlT7NWrF9W6YMXWfRNUSYSgCEGp8vzzzz/s2bMn1STI3Nzc0g9iNvPQxo1s1qwZvSSJs3S6wm0e+zIvapl6tWijNgLKsiU2ccQIAuDtAK8Cio/DUYn8vXuV5+3zYU6fLtz+Kk3GvKeOlBRr/s2FTp34/HPPUZIkhnt58UudTtnSS0uzrUHmQi+WtWvXslmzZgTAp556iv/880+pvypZlvnFF18wNDSUvr6+/OCDD5w3aRNUGYSgCEGpFphMJo4ZM4YA2L59e548edL1izXhwJeTktjLIk4vvvCCkrntqFaXdjssMbFw60sdS6/nHD8/6gHeqtMxXadTLBKtGKlOf61/QhUPrS+lrFnz7hx6feF71OnIjAxuX7GCiZZtsGSAe1q0KJyvWlzTBXJzc/n666/TYDAwPDycixYtcsvZnpGRwXvuuYcA2KlTJx48eLDUYwgqFiEoQlCqFUuWLGFwcDBr1qzJX375xbWL7ATDlJ7Ol19+mQDYpUsXnrVPojSbi5aT12bwayyQX6KiGACl9MluVXzUEjKOSuc7ypBPSCg+u768Dm3/ekDxK505Q3NyMj/X6VgLSrn8VwBet6+D5iI7d+60dnfs27evW8EVqo+mRo0a9PX15ZQpU4RvpQojBEUISrXj8OHDjIuLoyRJfOutt0reDikmyXL+/Pn09fVlRESE0qlR6zMxmRSLQ9u50T6Dv1MnUqfjDihFGoMALtfpCsexv6f6nL1F0r59YcXi7dtLLw6+vraPJ092T2QkSQkAWLmSFwEOslgrDQEuW7LEre/KaDRy0qRJ9PHxYc2aNfnVV1+V2VpJSUnh0aNH3ZqPoHwRgiIEpVqSnZ3Np556imqo6blz52xPsHeqFxOttGPHDkZFRdHb25sff/yx7WKnlmGxFwZ1LKNREQMo/UpuhZLn8eHUqco4ju6pFmtUI820/hO1bH9SUumEoE0bx+Lg7BoXt9nWQ6lnBoD9+vThOdUnVEr+/vtva4XpPn36MCMjo9RjqL6VkJAQ+vv7c+bMmSJvpYpR7QUFwB0ADgE4CuA1B68/DeACgF2WY1BJYwpBqR7IsszPP/+cvr6+rFevHtetW6e8UFJhSruF/p9//uGdd95JAHzsscd4/fp1p+cXuUdQEClJvNahA++zJGUOGDCAeTk5jkNuz54tvjhkcrIiVKdOWSsKu2SVlOORD/C/AL0B1gQ4r1kzpfpAKTGZTPzggw/o6+vLmjVrcsGCBW4JwpkzZ6zVkHv27MkzZ86UegxB+VCtBQWAHsAxAI0BeAPYDaCF3TlPA/ioNOMKQale7Nq1i02bNqVer+eECRNoLqnCsAOxMZvNHD9+PCVJYosWLXjgwAHnN3XkyJdlms1mjlWTMoODme6oFa/aLri4RTwjo6joeHnZnuOJasUlHQkJZEyM9fEBgEkWa6V3jx5uL+QHDx5kx44dqZZdcaeZlyzLnDlzJv39/RkaGsoFCxa4NReBZ6kQQQGQCmA9gP0AFgBIcPemduN2AvCb5vEYAGPszhGCchOQlZXFRx55hADYq1cvntf2cC+pz7yGlStXMiwsjAEBAZw/f37xNyyu+KVFsL7R6RgAMBzgWtWvomX37uIX8jNnlPEc9Jkv0+Go/EtJx4oVNo9NAKcB9PfzY3BQEGd//nnx23tOMBqNfPfdd+nl5cU6depw6dKlLl1nz+HDh5mYmEgA7NevHy9duuTWOALPUFGCchTA7QDCLFtUmwA87O6NNeM+COBzzeMn7MXDIiiZAPYAWAygQTFjDQawDcC2yMhIT3/OggpAlmV+8skn9PHxYd26dblm8eKiC5wLVZDT0tKYkpJCNe+l2MxvR4uoRrD2AWwOJVpq0rvv0mw0KtZHerqtH8X+2LFDGSs93fZ51S9SljDj1atdP1eSlDnYR4QBPHrrreyiWit33smMjh1t83hcFJfdu3ezdevWBMBnn32W165dc/XrtmI0Gvn222/TYDAwIiKCq1atKvUYAs9QUYKyxe5xIIB97t5YM44rglILgI/l388BWFPSuMJCqd7s2rWLzZs3t0aBGY1G2xNc+DVtNBqteS8tW7bkvn37XLu5tsCkJPEqwIcsmfV3hYbaVvotLlR4587CsifFlV1x91i1SilKWdzrHTrYNhwrKHAoKDQYaAb4fzod/Xx9WRPgN4DyvrWN0lwI8c3Ly+O//vUvSpLE6OhobtmyxbXP2o5t27YxJiaGAPjyyy+7XVhU4D7lKigA5gEYCeB9AG8BMFie9yrLjTXjl7jlZXe+HkBWSeMKQan+XLt2jU8//TTVMFN3uw3+9ttvDA8Pp5+fH2fNmuWaE9kuOkxOSeFHUVH0BhgBJXKqWEsjKKiw2GNiYukjvsp6+PkpAQGZmUro9Nq1Rc9p2dLGujrYpg0TAgIIgI/Vrs0rpcyuV1m3bh0jIyOp1+s5bty4oj8EXCA7O5tDhgwhALZp08a2e6eg3ClvQekCYDiAOQC2AzgJYJXFkT7Z3RtrxjcAOA6gkcYp39LunLqaf/e1t5YcHUJQbhzmz5/PwMBAhoaGFjbcKiUZGRm87bbbqDqRXd6nVy0hS3HIHQCbQAktfgugMTlZcc7r9UpDr127FN+Kuw73+Hhlu8zB9SaAP0OJ1vrZ8tj6erNmRcdq21axVuzDj4vZbjNKEsdFRFCv17Ohjw//KK5/fQnW4ZUrV/jYY48RADt37ux25eGffvqJYWFh9PPz4yeffCLCiyuICtnysl6gCEArAI97QlAsY/YGcNgiUm9YnvsvgD6Wf79rCQbYDWAtgJiSxhSCcmNx5MgRdujQgeo+fZGwYBcwm82cOHEiDQYDGzRowPXr17t+scZvc7VDBz5h6dOSlJTE40eP2haVVMvQq4u1q3W+2rZVtqfS0xXLxmCwlncxAewBMBCgZPnbw15USjpatHCpG+WWpUvZuHFj6nQ6vjV2rK2VUVI4t4avv/6awcHBDAkJ4cKFC0v9fZFkZmamtQZc3759efHiRbfGEbhOhQpKdTmEoNx45OfnW/fpmzdvzu3bt7s1ztatW9mkSRPqdDq+8cYbSi0wFWe/vu1eUxfMoKAgfvnll0rTLjX6TKcjW7cuXHjXrCl5wVdrh6kC1Lat1Zr42SIi0ByBluddFhRJKlncJInMzOTVK1f45EMPEVAqBlvDi0uIsLPn+PHj1giugQMHuv1D4IMPPqCXlxcjIiIKc5UE5YIQFCEoVQs3+siXhtWrV7NevXr08vLie++951ZdqKtXr3LAgAEEwPj4eKXEeil+faucOHHCGk12//3384LqzA4JUcQgMVGxOuwbfrly6HRWh/5/LZaJVlAkgONLM56zWmOS5LDk/VfNmzMwMJA1a9bkTz/95DjCroTvu6CggG+88QYlSWJsbCz37t1b6u+LJLdv386mTZtSp9PxP//5j6heXE4IQRGCUnVwY1F2h4sXL/L+++8nAHbt2pWnTp1ya5zFixezZs2aShmQiRMpu+GMNplMnDRpkjUf46fp022rAKvbVx07KhFazgTEkrFPQFncT52ixyyUkoTLUfFMg4GHNm5k27ZtCYCvjB5NY35+oYCU4vtetWoV69SpQ19fX3722Wdu+USuXr3KJ598kgCYmprKtLS0Uo8hcI4QFCEoVYdSbomUBVmWOXv2bAYGBjIkJIRff/21W+OkpaVZ9+nvqFGjMDO+NAue2cxdq1db8zGeDA/nJe2vf/XfiYm2GfZa53tgYGFxx/R0xbKxWDUu+1D8/EonJnq9kk2vFVL72mcFBcxNSODzltDp1NTUwjpepfy+z549aw2OePTRR93KWSHJefPmMSAggLVr13a9YrXAJYSgCEGpOriQdOhpjh07xqSkJALgI4884pbjVpZlfvTRR/Tz82PNGjW48H//c/1iza/0/ORkjn3zTer1etYD+KOjRTwtTbEC9uxx3MdeLTBpZ82oUV7j4SDKSysQpYkwa9/etq+LuuWlhj0nJ9tk+8+XJPr7+bFOnTrcsGFDyd+3g+0wk8nE8ePHU6fTsXnz5q7nB9lx8OBBq4D/61//citEWVAUIShCUKoW5exDcYTJZLJmWtetW5fLly93a5yDBw8yISGBAPjQQw+5VqPKwa/0bX/9xThLXkf/WrV4XruIq50hZdlx0mNKivM6YeV5tG1bpF+MjegFBXHf7t1s2rQpvby8lJbOxX3fJWyHrV27lnXq1KG/vz/nzZvn1veVk5PDwYMHs0jwgMBthKAIQRFY2LFjB1u2bEkAHDx4MK9evVrqMYxGIydMmEAvLy+Gh4fzu+++c35BMb/S83Nz+Z/Ro+nl5cVaBgPnSRLloKDC89LSlO0m7YLerp2y3WWfK7JiheuWR1naEet0hVte6nZYQIDt2BkZvHTwoHWbcOjQoYpfJSOjMJlSk7vjbDssIyODqampBMDnnnvO7cz4r7/+2roF9ttvv7k1hkBBCIoQFIGG3Nxcjh49mpIkMSoqimvXrnVrnD179lid0Y888ohza8WJVbZv3z5r6GxPSeJRdWHW6You/mpnRa2FEhJiE3lVru2G1TbBRqNj66lDB2uQgTElhaMs1kGv0FBmqecEByvz1PS8d7b9aTQa+eqrrxIAExIS3K6I8Pfff7Nly5aUJElEgZUBIShCUAQO2LBhA6OjowmAw4cPL75ApBMKCgo4fvx4enl5MSwszO3e6iaTiR9Nm8YgnY6+AN+B0qOkiGWRmlpYePLMGWV7TN0qMpuVGmHlJSaqmDm7jzYx0iKGswAaALYEeNL+fNXhn5KijOtEeL/77jsGBgaydu3aXLNmTak/Y5K8fv26NQqsV69ebpXVv9kRgiIERVAM169f59ChQwmA0dHRpcuO17Bnzx7Gx8cTAO+77z7n3QodLZpGIxkSwjSADxgMBMAYgGu0i++iRUr5FnWxTkxULBPtGLfe6jkBiYuzfZyaWuiQL+VYqwCGALwF4HZH56iBBiWEGP/999+MjY2lXq/n9OnT3RJvWZY5a9Ysent7MzIyklu3bi31GDczQlCEoAhKYM2aNWzUqBElSeKIESPcytg2Go1877336Ovry5CQEH6u9hHRUpwjeu9emwV22UcfsVFEBAGwP8B0zS9+myMxsfCXvb2/xdXDUWfIgABbZ7skKUKSkeH6lpq2vExQEPe1asVIgAEAf7UXk9RUl3wqpNIXp0+fPlSz6931q/z1119s2LAhvb29+emnn7o1xs2IEBQhKAIXuHbtmo214q5v5fDBg0y1+ES6du3Kw4cPF75YXF6G2az4QgDlr9nMnOPH+RZAHyh5JZPgYBtMr1fGyMgouqA7yzlxx88iScpWlzMrSCt62ntYMu0zAN4KZQtsAaD4gjIzC6PaXAwpN5vNHDt2LNUCk+fOnXPru7p48SJ79epFQOmLI8rhl4wQFCEoglLw+++/W30rQ4YMKV0kmMUCMev1nNW0KUNCQujj48O3336b+fn5zhdNo1GxVFSrxdLR8SjAe6Bkvzfz9eWyL79ULBM1wdJkKhpGrJao1zrOnfWv1wqUO1aOes/0dMWS0elsS8lIkrXczBWAXS3vZ9r48YVbgNr8FtWnUgILFy6kr68vGzZsyD179pTym1YwmUx8/fXXqTr9RWixc4SgCEERlJLr16/zpZdeoiRJbNCgget5K3YWSMbu3XzwwQcJKE28NmzYULo8HJNJEQ+djsubN2ezZs0IgHfecQcPrF+vjHH2bFGLo23bQmshPt42tLekw90oseBgMj+/MKTY/nWDQem9kpLCXL2efWvVIgBOatRIea1TJ9uMfPstL/vPzfJ4659/sm7dugwKCnI7v4gkv//+ewYGBjI8PJx//PGH2+Pc6AhBEYIicJNNmzYxNjaWAPj444+XHBVUjAXy008/MTIy0rq18s8//7g+Cc1Cmp+fz/fff5/BwcHU6/V88cUXeeH8eVsLRS3R4q6l4azbY0nH4sXFWzn+/opYpaaSp0/TmJDAfhZL5QNVRNS6Zo6sN22XSKPR5nM+c+oU27ZtS51Ox48++siNb1rhwIEDbNKkSWFSpqAIQlCEoAjKQF5eHseOHUuDwcCwsDAuWLDAeXRRMRbItWvXOHr0aOr1etauXZtz5851uynU+fPnOWTIEOp0OoaEhHBiVBRztIu26o8BnFsnHTvavh4UZFtqpbSHK6Vd9HoyNpYEaAT4kEVUPmrSpDDpUfu5mM025V1oMChbg3a+qGvXrlmd9S+99JJbVaZJ8tKlS7ZJmaJkiw1CUISgCDzA7t27rU287rzzTp48edKtcXbt2sVOnTpRLQfi7t4/qSRF9u7RgwBYH+AcWGp4GQxKZ8hdu2wXc60fJT7e8/3s3TgKAN5rEZUvv/yy6Js8e9ZWpDp1KhotZxEgk8nE4cOHU20XkJOT4/wDLEb8jUYjR40aRQC87bbbXO/geRMgBEUIisBDmEwmfvjhhwwICGBAQACnTp3qVsa12Wzm7NmzWatWLer1er700kvMyspyb1KyzLWtW7ODZVFuCfC7Fi0om82Kk7y4xbxdu6LWRGmtk9JWL7b3z/j7kwYDc5OT2aNHD+r1ev78889F3p/VL6MNk1ZLuTiw8qZOnUpJktipU6fitxddKK0/Z84cenl5sWnTpkpPHIEQFEeHEBRBWTh58iR79+5NAGzfvr3b3SEvXrzIZ599lpIk8ZZbbuH8+fPd2wYzmylnZnLRwoVs3qSJdV7L582j7GzrqSwWSmn8NI6EKiBAKcFvsRCuXbvG9u3bM8DfnzvsP0+tJWE02ka5FbO1tXjxYvr4+LBFixaOI7dcLK2/YcMG1q5dm6GhoVy1alXpv5sbDCEoQlAE5YAsy1y4cCHr1KlDnU7Hl156ye3+HVu3brVup3Xu3Jk7duxwe15Go5Fz585lVFQUATAJ4EpAERZtKG9yMrl6tfuCUpZDpyO1OR8WiyOjY0c2AFjP25vpjkTAbFa2vLSi6KTHytq1axkcHMwGDRrw4MGDti+WIu/l+PHjbNmyJQ0GA2fNmlXar+SGQgiKEBRBOXLp0iVrifT69evzhx9+cGscs9nMzz//nGFhYZQkiYMHD+b58+fdnlf+jh38BIpvBQA7t2rF35Yto6x2X1TbENsv9mr3SO05HTt6tqwLUJhzk5ZWaHEA3A3QH2BKYmJRh7jWqgCU60qw6Hbu3Mnw8HCGhYVx586d9h+6yyHcWVlZvOOOOwiAo0ePdtvpX90RgiIERVABbNq0iXFxcQTAPn36uO20v3z5MkeOHEm9Xs/Q0FB++OGHLCgoKP1Aluz7PIAzfX1Zv359AmBHgEtVi8XRNlhKSmHhSaNR8VWoNcS021VlFZTTpx3XBdPrOb95cwLgG2+8YfuetFaF6k9xgUOHDrFBgwYMCQnhpk2bSv9ZWjAajXzhhRcIgH379nWroGh1RwiKEBRBBVFQUMDJkyfT39+f/v7+nDhxopIh7wb79++3hq/GxMS4l7Snyb7Py8vjJx9/zIY+PgTAW3U6LpQkmlQrRf2r3RZLTHTcf2X7dqWbo05XsiO/OPFxZB0FBSn3k2UOGDCAOp2Of/75Z+H7MZuLz6Yvwdo4deoUmzRpwoCAAP7++++l/ywtyLLMDz/8kJIkMTEx8aarWCwERQiKoII5deoU+/btSwCMjY11uy6YLMv8+eef2bRpU6ol191tiatSkJfHuR9+yOYWK6BJkyb8dNQo5hTnZO/Y0dZCCQ4uLK8SFORauZb4+JLPUUXG4hPJyspiREQEW8XEMF/1tzirhVZCxBapNOyKjY2ln58fV65cWabP8fvvv6evry+bNm3Ko0ePlmms6oQQFCEogkpi6dKlbNSoEQGwf//+TE9Pd2uc/Px8Tp06laGhodTr9Xz++edLLojo7Be72UxTejq/W7zYWna/tpcX34ClsrH9kZamHHv22LYALsshSUqF5G3bCq2coCDFqrLM8UdLd80Po6OV91OcI93FiC2SPHfuHFu3bk1fX98yR21t3LiRNWvWZHh4uNuRftUNIShCUASVSE5ODt966y36+PgwMDCQ77//vns+EZIXLlzgsKFDaTAYGBQUxHfeecdx8p6zX+x2r8kmE9euXs377riDkiTRIEnsD3Cz6meRpMI+92Thou4sbNjRa/YZ+dpy+9rzExMV0UpPJw0GdgMYDvDasWOF87cXSlm23Qorwcl+4cIFtmrVin5+fmXa/iLJgwcPsmHDhgwMDCyz1VMdEIIiBEVQBTh69Kg1dyU2Nta9xcciBgf1et5rKa5Yv359zp071zbB0tkvdvvXMjKsi/Gx+HiOHDaMwYGBBMB2AGf7+jJbrcGl7Q65dq1jMenYkczNVV5PTFTEonNnxdGv+k30eudFKPV6q9WyCZbSLNOnl/i5lLTlpeXcuXOMjY1lQEAAN2/eXPrvQkN6ejrj4uLo5eXFhQsXlmmsqo4QFCEogiqCLMv86aefrOXx77///tJFg9mJwe/ff2/NX4mLi+Mvv/yiJEY6y7Gwf+3MmSKCcFWn44zoaLZo3JiA0m1xmCRx33ffFWanm822DnytoKSk2Foder2Sd6K9r71fpm1bx/4YnY5NGzdm7969Xf5cnG15acnIyGB0dDRr1KjBvXv3uv49OODy5ctMSUmhJEmcMWNGmcaqylR7QQFwB4BDAI4CeM3B6z4AvrG8/ieAqJLGFIIiqExyc3P59ttv09/fn76+vnzrrbdcC0F1IBSyLPOb//2PjRs2JAB269aNWzZtclqaxKYHibbwot0hp6VxfZs2fFSS6G2xFOIBTouO5vmzZxV/x2+/2V5XXIHI1attt6u0vhhJKvSp2Hee7NiRQ198kf7+/oVWmH3plVIkKdpz/Phx1qtXj3Xr1uWJEydcvs4ROTk5vOeeewiA48ePd7v4Z1WmWgsKAD2AYwAaA/AGsBtAC7tzXgDwieXf/QB8U9K4QlAEVYHTp0+zX79+BMAGDRpw4cKFJS9CjvqCpKYyX6/n/0VHMywsjAB4ryRxb/v2zrd/7BMF27UrtCwkSfFjpKWRa9fyApQy87dahMVgMPCu3r25wM+P17W+kNRUxyLVoYMiQOrcTabCkvTa87Zvt32ckcG3x48nABbk5xeGDmvvp9b3crXPjB379u1jaGgomzVrVqZkUlLJVXnyySepJkDeaKJS3QWlE4DfNI/HABhjd85vADpZ/m0AcBGA5GxcISiCqsT69et56623EgBTUlJKFzFkt91z7c8/OV6nYzBACeDjDz7II0eOOL5W+8u+UydlkVcfp6QU7QRpOfYGBPDV0aNZv04dAkqf+H4Avx4zhpf+/lsRjrg42+t0usJ+JykphQUf27WzPS8hQdkOkyTl/gUFfL1ePQKgrPa111pApdjicsYff/xBX19fdujQocwJi2azmS+++CIB8Lnnnruhsuqru6A8COBzzeMnAHxkd84+APU1j48BqO1grMEAtgHYFhkZ6bEPWCDwBCaTiZ9++qm19MqAAQOYmZlZ8oX22z0Wi+WiXs9X6tenn58f9Xo9Bw0a5NhfU0wnRGZmFp9jYlnEzWlp/B3gswDDLJaLHmCXkBBOHD+eu3S6woz8Dh0Kx9PrbX0e9u2JNc773Ph4RgJMUK/LzLTNsNducbljpWiu+fHHHylJEvv161dmy0KWZb722msEwGeeecatqtRVESEoDg5hoQiqKleuXOErr7xCLy8vBgYGFh8arKU4UZBlZmZmcvjw4fT29qaXlxdfeOEF1/qmq6G4jgRFXcRlWbEidDqaO3TgZr2ebwBsYxEXAKwL8BGAU7y9uRFgNlBooRgMSgixI/8NwC0Ab7eMsxoorN2lrTisZsy7Eenl6JqJEycSAN95552Sry/xI5T573//mwD4xBNP3BCiUt0FRWx5CW5Kjhw5wnvvvZcAGBkZWXKnSHvsROb06dN87rnnaDAY6O3tzaFDhzItLc35tSaT4kdRt6rUHBH7joqqg16zyGekpfGL6Gj2BxipERgJYFSDBryjWzcOe/RRjpMkTgc4H+CnAP8L8EUofV0A0F+n42c6nbIlp/pgMjOLRnW5E+nl4BpZltm/f39KklSmHvVa3nnnHQLgo48+Wu07QFZ3QTEAOA6gkcYp39LunBftnPKLShpXCIqgurBmzRqrfyUhIYF//PFHyRc5+bV+4sQJDho0iAaDgT4+Phw2dCjTdu4sFAlH/UZK2krSOspVCyQ/X+kYaRkrIyiIP+h0HBcZyf5hYWwLMFivtwqN9qih1zPJ35+zJIlZSUmKiNn7d+yjuuy3/hy1E7anmOiw7OxstmnThqGhocX7n0rJu+++S0CpmFCdLZVqLSjK/NEbwGHLVtYbluf+C6CP5d++AL6FEja8FUDjksYUgiKoTphMJs6dO5f1LM7pBx54wPlC58Kv9ePHj3PQwIE0WEKCh9StyxNHjrjeb0QrMvZterVbWYGBSq6LusBrrQuABTodz8fF8VC7djyj1zO/Y0fbkGK1h7y9VZKR4bj/vGpZubr9VYxYHj9+nDVq1GCHDh3crmxgj7qd9uSTT1ZbR321F5TyOISgCKoj169f57hx4xgQEECDs6zwfgAAFmVJREFUwcDhw4c7rnbral7G2bM8odfzOYBeUEKBn5EkHlKFobh+I/YWkMlkG8rbvn2huADKY9XSUZ3q9g7/xEQlRDkz03b8kJBCayklxfY9FWc5uZnoaM+3335LAHzrrbfcut4R//3vfwmAzz77bLUUFSEoQlAENxiZmZkcPHgwdTodg4ODOWHChKKhrq5EPGmE50zHjhw+bBh9dTpKAB+sWZN/aUvHa3G0YKvNsvbsUba7tLW7VItCu2XVoYOtoEhSYV5Kamph6X1VeOytEmdO+DIkOtrz5JNPUqfTcevWrW6PoUWWZb7++usEwJEjR1a7PBUhKEJQBDco+/fvt2ZmR0REcPbs2aXfn7cTnrMZGXx9xAiGhIQQALt3787ffvvNduFztGCrC7y2rH1AQGEk1p49tiJkX0CybVvb19Us+OKEwZkVohU3eyuglKHFly9fZr169diqVSvmnTpVJnEq/PhkjhgxgoCSUV+dEIIiBEVwg7Nu3TomJCQQAFu0aMElS5aU+ZdvVlYW33vvPavfpk2bNpw/f36hP8F+Ybb3o2j9MMHBitAEBxeKg32JlfR02yrGahKjI+e6duvMXmxUYdOOoy1qWdrQYpI///gjAfBt+yKZZcBsNvOpp54igGpV+0sIihAUwU2ALMtctGgRmzVrRgBMTEwsc2l2kszLy+Ps2bMZGxtLQKlu/P777/PKlSvKCaqwmM1O64LZCExKCllQUBiOrAqCNgtekmwjzVTsRUENY9YmZGrLuWgDC5w16HJmtZw9y/skiUEAzzkLVCglRqORffr0oSRJXLRokUfGLG+EoAhBEdxEGI1GfvbZZ4yIiCCgdHnctm1bmcc1m838+eef2bVrVwJgUFAQR44YwaMdOhQu7gUFih9Er1eqDmsd6zpdodWg9btoF3LtFpfW/+Ks/L46jtY/o612rO2PYr+FlpnpWkSYLPNgfDz1AIfXq+eRbS+VnJwcJicn08vLi2vWrPHYuOWFEBQhKIKbkJycHE6ePJm1LH1THnjgAe7fv98jY2/bto39+/enwWCgBPBugL/pdDRnZNiKhDaMNzPTNUe5s+0s0rFfxVGPl+KqLauVitVIs3btbB3/TsKkn37kEfr5+fHixYse+RxVLl++zBYtWjA0NJR///23R8f2NEJQhKAIbmKysrL4n//8h0FBQdTpdHziiSc81gM9PS2Nb0VGMtySkNi8eXNOmzatcDvMntI4xEtoYVysZeNKVJe9vycoSLGgSuj2uG/fvnJzpB8/fpxhYWFs3Lix41DwKoIQFCEoAgEvXLjA0aNH2xSLPHXqVNkHNpuZd+oU58+bx44dOyrlUvz9+eyzz3Lnzp0uj+Fu6Xm3xpBlW3+PM3+N3Zg9evRg06ZNyyXcd/PmzfTx8WFqairz8/M9Pr4nEIIiBEUgsJKRkcGhQ4faFIsstqaXG2zbto0DBgygn5+fNThg7ty5vHbtWtGTtdtPpYy8KjOO+rHYO+od+FY+//xzAvCIX8oRX3/9NQFw0KBBVTJHRQiKEBSBoAinTp3i4MGDC2t6DRvG9PR0j41/6dIlTp06lc2bNycABgYGcuDAgdy4caOyUGoX7OIc8OWNM39NMRFhFy5cIABOnDix3KalJj7OmjWr3O7hLkJQhKAIBMVy/PhxDhw4kHq9nr6+vhwxYgQzMjI8Nr4sy9ywYQOfeeYZBgQEWH0t777xBtMchQhXxq9yR9tlTvwykZGR7N+/f7lNx2QysWfPnvT29naeoe+JrcJSIgRFCIpAUCLHjh3jM888Q71eTx8fH7744os8ffq0R+9x9epVzp49m8nJyQRAHcAeksQ5zZrxyqFDrjngK5Ji7t2jRw927NixXG998eJFRkZGMioqipcuXXI8NzeSNMuKEBQhKAKByxw9etRa3t7Ly4vPPvssjx075vH7HD58mGPffJPRUVEEQB8fH95///1ctGhRYV2ySlo0SyIlJYVdunQp9/ts2bKFBoOB9913X1F/iocKYJYWIShCUASCUnPq1Cm+8MIL9Pb2pl6v5+OPP+6xPBYtsixzy5YtHDZsGOtYetT7+/vz4Ycf5v8++YRZruSIVCCyLLNRo0a89957K+R+U6ZMIQB++umn9hPxWAHM0iAERQiKQOA26enpHDVqFP39/SlJEu+///5yi3AymUxcvXo1n3vuOYaHhxMAvSWJvSWJnzRpwjS1dXElboOtXbuWAPjll19WyP3MZjNvu+02BgQEFLUUhQ+lahxCUASC0nHhwgW++eabDA0NJQDefvvtXLNmTbmFtppMJm7YsIGjXnqJjSIjqXZzbNu2Lcc0aMA1Oh3zkpMrdBssOzubnTt3ZmhoKHNycirsvqdPn2ZISAhTU1MrvYeKEBQhKAKBx8jKyuLEiROt21MJCQn8/vvvy3Whk2WZ+/bt48SJE5mSmEiDRVz8AN7epQsnTJjAzZs3e6yzoiOys7PZrVs36nQ6LliwoNzuUxxz5swhAM6cObPC761FCIoQFIHA4+Tk5PDjjz9m48aNraHAs2bNYm5ubvneWJZ5NSmJP+t0HG7pU6JaLwEBAUxJSeGoUaO4YMEC7t+/v8wic+nSJc6cOZNxcXHU6XScP3++h95I6ZBlmbfddhsDAwM9mohaWsoiKJJy/Y1HfHw8t23bVtnTEAiqPSaTCd999x3ee+897NixA+Hh4Rg6dCiGDBmC2rVrl89NZRm4cAEIDwckCRcuXMD69euxbt06/PXXX9i1axfy8vIAAF5eXmjevDliY2PRqFEjREVFISoqCrVq1UKNGjVQo0YNGAwGFBQUoKCgAJcvX8ahQ4dw8OBB7NixA8uXL0d+fj7i4uIwbtw49O3bt3zekwscP34csbGxePTRRzF37txKmYMkSdtJxrt1rRAUgUDgCiSxbt06TJ48GcuXL4efnx+eeuopjBw5Es2bNy/+Qjtx8ARGoxH79+/H3r17sX//fuzfvx8HDx7E6dOnUVBQ4PI4UVFR6NOnD55++mnceuutkDw0v7Lwr3/9C5MnT8b27dvRtm3bCr+/EBQHCEERCMqPAwcOYMqUKfjqq6+Qn5+Pu+66C6NGjUK3bt1sF2VZBrp1AzZtApKSgLVrAZ2u3OYlyzLOnj2LU6dO4Z9//sGVK1dw+fJlmM1meHt7w9vbG4GBgWjevDmaNWuGgICAcpuLu2RlZaFx48bo1KkTli5dWuH3F4LiACEoAkH5c/78ecycORMzZ87EhQsX0Lp1a4wcORL9+vWDn58fcO4cUL8+YDIBBgOQlgbUqVO6m5SDhVPVmThxIsaMGYMtW7agY8eOFXrvsghK+f1UEAgENzzh4eH4z3/+g9OnT2P27NmQZRkDBgxAgwYNMGbMGJzOy1MsE4NB+RseXvKgsqwIEVlo4dSvD3Ttqjy+CRg6dChq1ar1/+3df2xVZZ7H8feX8qNBW8RBKq0VKFAH0LogyG6JFFDwxyYwIyKzRlojiz+WjcHd6E6CCcluVNZNNkuywtoYQ40ZmZWQtJ2C4xShhEpFWoUBJkirE6nlR8dE6AKttn32j3uhFVp6gXPPw739vJKbe057eu+Xb0774Zzz3OewZs0a36VcEQWKiFyz1NRUnn76afbv38/HH3/MrFmzeOONNxibk8Mvhw+n8je/wW3f3vcRxsUBcuJE5HRZe3vkubk5lH+PbzfeeCPPPfccpaWlNDQ0+C4nZgoUEQmMmTFnzhw2b97M119/zcsvv8yu6mrmPf44EydNYu3atXz//fe9v0Bz808DxOzKj3CSxIoVK0hJSeHtt9/2XUrMFCgiEhe33347r7/+OkePHuXdd9/lpptuYuXKlWRmZrJs2TL27t3LJddwR478aYBkZEQu5Dc2wo4d/eYaCsCoUaOYO3cumzZturRP1ykFiojEVWpqKkuXLqWmpoa6ujqefPJJNm7cyPTp05k6dSrr1q3rOmoxuzRABgyIBEs/CpPzFi1aRH19PYcOHfJdSky8BoqZ3WxmfzCzI9Hn4b1s12FmX0QfZWHXKSLBmDJlCsXFxTQ1NbFu3TrMjBUrVpCZmUlhYSFVVVU4s34bIBebO3cuANXV1Z4riY3vI5RfA9uccxOAbdH1npxzzv1V9LEgvPJEJB6GDRvG888/T11dHbW1tRQVFVFaWsrs2bPJzc3ltdde4+jRo77L9G7cuHHccsstfPrpp75LiYnvQFkIlESXS4BfeKxFRDyYOnUq69ev59ixY2zYsIHMzExWrVrF6NGjmTdvHu+99x5nzpy5+jfoPgw5wZgZ48eP55tvvvFdSkx8B0qGc+5YdPk40NsnnlLNbK+Z1ZiZQkckCQ0dOpSioiKqqqpoaGhg9erVNDQ0sHTpUjIyMigsLOSjjz6ivb099hdNgs+xZGRkcPz4cd9lxCTugWJmlWZ2oIfHwu7bRWe57O2/EKOjn9x8AvgvMxvXy3s9Ew2evc39ZLy6SDLKyclh9erV1NfXs3PnTp544gnKysp48MEHycrK4oUXXmD37t19j366eBhyAv5d6OzsZEAcp6sJktepV8zsMDDbOXfMzEYBO5xzl5llDsxsA/A759ymy22nqVdEkktraytbtmzh/fffp7y8nLa2NkaPHs3ixYtZvHgx06dPv3RyR+ciRybn5xJLwKHH+fn5DB06lMrKylDeL5GnXikDiqLLRUDpxRuY2XAzGxJdHgHMBBJjDJ2IBCY1NZVHH32UDz74gJMnT1JSUsLkyZNZu3YtM2bMYMyYMbz44ovs2rWLjo6OyA/1NAw5gZw7d479+/eTm5vru5SY+A6UNcA8MzsCPBBdx8ymmdn5j4dOBPaa2T5gO7DGOadAEenH0tPTKSwspKKighMnTrBhwwbuvvtu1q9fz3333UdWVhbLly+nvLycs62tCTsMedu2bZw5c4aFCxf2vfF1QLMNi0jSaGlpYcuWLWzevJmtW7fS0tJCamoq999/Pw8//DCPPPIIY8eO9V1mTDo6OigoKODw4cN8++23DB48OJT3vZZTXgODLkZExJe0tDSWLFnCkiVL+OGHH6iqqqK8vJyKigoqKioAuOOOO5g/fz7z58+noKCAtLQ0z1X37JVXXqG6upqSkpLQwuRa6QhFRJKec44jR46wdetWPvzwQ6qqqjh37hwDBw5k2rRpFBQUMGvWLGbOnMmwYcO81tra2srKlSt56623WL58OcXFxaG+v26w1QMFioj0pq2tjerqaiorKy/cp/7HH3/EzLjrrruYOXMm+fn53HPPPeTm5pKSkhJKXbW1tTz77LPU1tby0ksv8eqrrzJo0KBQ3vs8BUoPFCgiEquzZ8+ye/dudu3aRXV1NTU1NbS0tABwww03MGXKFPLy8rjzzjuZPHkyEydOZMSIEYHcg76xsZHS0lJKSkr47LPPuPnmm3nnnXe8XYhXoPRAgSIiV6ujo4NDhw5RW1tLbW0tdXV1HDhwgNOnT1/YJi0tjfHjx5OTk0NWVhaZmZlkZmYyfPhw0tPTSU9PZ/DgwXR2dtLR0UFbWxsnT57k+PHjNDU18fnnn7Nnzx6ampoAyMvLo6ioiGXLlnk97aZA6YECRUSC5JyjsbGRgwcP8uWXX1JfX099fT1fffUVTU1NF45oYjVhwgTuvfdeZsyYQUFBAXl5eXGq/MpolJeISJyZGdnZ2WRnZ/PQQw9d8v2WlhaOHTvGqVOnOH36NC0tLbS1tZGSksKAAQMYNGgQGRkZ3HrrrWRkZDBkyBAP/4r4UqCIiAQgLS3tuh2CHBbfn5QXEZEkoUAREZFAKFBERCQQChQREQmEAkVERAKhQBERkUAoUEREJBAKFBERCYQCRUREAqFAERGRQChQRETipbMTTpyAJJ2E92IKFBGReOjshDlz4LbbYPbsyHqSU6CIiMRDczN88gm0t0eem5t9VxR3ChQRkXgYORLy82HgwMjzyJG+K4o7TV8vIhIPZrB9e+TIZOTIyHqSU6CIiMTLgAGQkeG7itDolJeIiARCgSIiIoFQoIiISCAUKCIiEgivgWJmi83soJl1mtm0y2z3kJkdNrN6M/t1mDWKiEhsfB+hHAAeBXb2toGZpQBvAg8Dk4C/M7NJ4ZQnIiKx8jps2Dn3JwC7/Pjse4F659xX0W03AguBQ3EvUEREYpYIn0PJAo52W28EZvS0oZk9AzwTXW0zswNxri1RjAD+4ruI64R60UW96KJedLnjan8w7oFiZpXArT18a5VzrjTI93LOFQPF0ffd65zr9bpMf6JedFEvuqgXXdSLLma292p/Nu6B4px74Bpf4lsgu9v6bdGviYjIdcT3RflYfAZMMLOxZjYY+BVQ5rkmERG5iO9hw780s0bgb4AKM/t99OuZZrYFwDnXDvwj8HvgT8D/OucOxvDyxXEqOxGpF13Uiy7qRRf1ostV98JcP7mTmIiIxFcinPISEZEEoEAREZFAJHyg9DUti5kNMbPfRr//qZmNCb/KcMTQi1lmVmdm7Wb2mI8awxJDL/7JzA6Z2X4z22Zmo33UGYYYevGcmf3RzL4ws13JPBNFrNM4mdkiM3OXmxIq0cWwXzxlZs3R/eILM/v7Pl/UOZewDyAFaABygMHAPmDSRdv8A/A/0eVfAb/1XbfHXowB8oB3gcd81+y5F3OAodHl5/v5fpHebXkB8KHvun31IrpdGpHpoGqAab7r9rhfPAX895W8bqIfoVyYlsU59wNwflqW7hYCJdHlTcD91sdcLwmqz1445/7snNsPdPooMESx9GK7c+5sdLWGyOebklEsvTjdbfUGIFlH6sTy9wLg34B/B1rDLC5ksfbiiiR6oPQ0LUtWb9u4yBDkU8DPQqkuXLH0or+40l4sA7bGtSJ/YuqFma0wswbgDeCFkGoLW5+9MLOpQLZzriLMwjyI9XdkUfS08CYzy+7h+z+R6IEick3M7ElgGvAfvmvxyTn3pnNuHPAvwCu+6/HBzAYA/wn8s+9arhPlwBjnXB7wB7rO9PQq0QMllmlZLmxjZgOBYcB3oVQXLk1R0yWmXpjZA8AqYIFzri2k2sJ2pfvFRuAXca3In756kQbcCewwsz8Dfw2UJemF+T73C+fcd91+L94G7unrRRM9UGKZlqUMKIouPwZ87KJXnJKMpqjp0mcvzGwK8BaRMDnpocawxNKLCd1W/xY4EmJ9YbpsL5xzp5xzI5xzY5xzY4hcW1vgnLvqyRKvY7HsF6O6rS4gMlPJ5fkebRDAaIVHgC+JjFhYFf3avxLZEQBSgQ+AemAPkOO7Zo+9mE7kXOkZIkdpB33X7LEXlcAJ4Ivoo8x3zR57sRY4GO3DdmCy75p99eKibXeQpKO8YtwvXo/uF/ui+8XP+3pNTb0iIiKBSPRTXiIicp1QoIiISCAUKCIiEggFioiIBEKBIiIigVCgiIhIIBQoIiISCAWKSIjM7DEzqzGzfdF7j9ziuyaRoOiDjSIhMrOfOee+iy6vBv7inHvTc1kigdARiki4njKzPWa2j8jN35L5nhvSzwz0XYBIf2FmhURubDTXOfd/ZraTyFxJIklBRygi4bkL+CQaJouAfOCPnmsSCYyuoYiExMwmA5uJ3DX0I+Bx51yu36pEgqNAERGRQOiUl4iIBEKBIiIigVCgiIhIIBQoIiISCAWKiIgEQoEiIiKBUKCIiEgg/h8d5vhaJzwr3AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Determine the errors in the model parameters by bootstrap re-sampling\n", "nsamp = 1000 # Number of bootstrap samples (can be much greater than the number of data points)\n", "asamp,bsamp = np.empty(nsamp),np.empty(nsamp)\n", "for isamp in range(nsamp):\n", "# Select a bootstrap sample using random sampling with replacement\n", " cut = np.random.choice(nx,nx,replace=True)\n", " x1,y1,y1err = x[cut],y[cut],yerr[cut]\n", "# Fit straight line to each bootstrap sample\n", " p0 = [a0,b0]\n", " result = minimize(chisqstraightline2,p0,args=(x1,y1,y1err))\n", " a1,b1 = result['x']\n", "# Save the best-fitting parameters for each bootstrap sample\n", " asamp[isamp],bsamp[isamp] = a1,b1\n", "\n", "# Plot the original chi-squared contours, as well as the best fits for each bootstrap sample\n", "# The two distributions should roughly agree\n", "plt.contour(alst2,blst2,chisq2,levels=clst,colors='k')\n", "plt.plot([a],[b],marker='o',markersize=5,color='black',linestyle='None')\n", "plt.scatter(asamp,bsamp,marker='o',s=5,color='red')\n", "plt.xlim(amin,amax)\n", "plt.ylim(bmin,bmax)\n", "plt.xlabel(r'$a$')\n", "plt.ylabel(r'$b$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimated bootstrap errors are then:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Bootstrap error in a = 0.04732010550152727\n", "Bootstrap error in b = 0.2631568374178819\n" ] } ], "source": [ "# Bootstrap error is the standard deviation of the best-fitting parameters to the samples (no re-scaling needed)\n", "s = stats.describe(asamp)\n", "print('Bootstrap error in a =',np.sqrt(s.variance))\n", "s = stats.describe(bsamp)\n", "print('Bootstrap error in b =',np.sqrt(s.variance))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity\n", "\n", "In this Activity we will return to our previous analyses of Hubble and Lemaitre's distance-velocity datasets and determine **bootstrap errors** on our measurements:\n", "\n", "- Find the bootstrap error in the correlation coefficient\n", "- Find the bootstrap error in the slope $H_0$\n", "- How do these compare to your previous measurements?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Combining and propagating errors\n", "\n", "A common situation in statistical analysis is, we have obtained measurements and errors of some variables, and need to determine the error in a function of those variables.\n", "\n", "Consider for example a **linear function of independent variables** $(x,y)$ with coefficients $(a,b)$:\n", "\n", "$z = a \\, x + b \\, y$\n", "\n", "The variances combine as:\n", "\n", "${\\rm Var}(z) = a^2 \\, {\\rm Var}(x) + b^2 \\, {\\rm Var}(y)$\n", "\n", "Now consider a **non-linear function of a single variable** $x$, $z = f(x)$. An approximation of the propagated error at $x = x_0$ is:\n", "\n", "$\\sigma_z = | \\frac{df}{dx} (x=x_0) | \\times \\sigma_x$\n", "\n", "Finally, a **non-linear function of 2 independent variables**, $z = f(x,y)$:\n", "\n", "${\\rm Var}(z) = \\left( \\frac{\\partial f}{\\partial x} \\right)^2 {\\rm Var}(x) + \\left( \\frac{\\partial f}{\\partial y} \\right)^2 {\\rm Var}(y)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Activity\n", "\n", "Here are a couple of examples to practice error propagation:\n", "\n", "_A galaxy of absolute magnitude $M = -20$ is observed to have apparent magnitude $m = 20.0 \\pm 0.2$. Assuming $m - M = 5 \\log_{10}(D_L) + 25$, what is the luminosity distance $D_L$ and its error?_\n", "\n", "_The total mass of a binary star system (in $M_\\odot$) is given by Kepler's law, $M = a^3/P^2$, where $a$ is the mean separation (in A.U.) and $P$ is the period (in years). The $\\alpha$ Centauri system has a period $P = 79.9 \\pm 1.0$ years and mean separation $a = 23.7 \\pm 1.0$ A.U. What is the total mass and error?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fisher matrix\n", "\n", "The **Fisher matrix** is a mathematical technique to propagate the statistical errors in a dataset to the errors in the model parameters describing the dataset. The matrix is given by:\n", "\n", "$F_{ij} = \\sum_k \\frac{\\partial m_k}{\\partial p_i} \\frac{1}{\\sigma_k^2} \\frac{\\partial m_k}{\\partial p_j}$\n", "\n", "where $i,j$ label the model parameters $p_i$, $k$ labels the $N$ data points, $\\partial m_k/\\partial p_i$ is the partial derivative of the model at point $k$, with respect to the parameter $p_i$, and $\\sigma_k$ is the error in data point $k$.\n", "\n", "After the Fisher matrix has been evaluated, the **covariance matrix of the parameters is then given by the inverse of the Fisher matrix**,\n", "\n", "$C = F^{-1}$\n", "\n", "The error in $p_i$ is then forecast as $\\sqrt{C_{ii}}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Monte Carlo simulations\n", "\n", "A **Monte Carlo simulation** is a computer model of an experiment in which many random realizations of the results are created and analysed like the real data. This allows us to determine the errors in our measurements as the **standard deviation of the fitted parameters over the realizations**.\n", "\n", "## Activity\n", "\n", "Run a Monte Carlo simulation of Hubble's distance-redshift investigation, and hence determine the error in $H_0$." ] }, { "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }