{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Estimating the Error Variance\n", "\n", "Author(s): Paul Miles | Date Created: July 19, 2019\n", "\n", "Included in the [pymcmcstat](https://github.com/prmiles/pymcmcstat/wiki) package is the ability to estimate the error variance as part of the sampling process. Furthermore, when using multiple data sets to inform parameter values, you can estimate error variances for each data set separately. \n", "\n", "For more details regarding how to estimate the error variance please refer to:\n", "- Smith, R. C. (2013). Uncertainty quantification: theory, implementation, and applications (Vol. 12). SIAM.\n", "- Marsaglia, G., & Tsang, W. W. (2000). A simple method for generating gamma variables. ACM Transactions on Mathematical Software (TOMS), 26(3), 363-372. [https://doi.org/10.1145/358407.358414](https://doi.org/10.1145/358407.358414)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9.0\n" ] } ], "source": [ "# import required packages\n", "import numpy as np\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt\n", "import pymcmcstat\n", "print(pymcmcstat.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Model and Sum-of-Squares Functions\n", "- Note, the sum-of-squares function is designed to loop through the data sets." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "# define model function\n", "def modelfun(xdata, theta):\n", " m = theta[0]\n", " b = theta[1]\n", " nrow, ncol = xdata.shape\n", " y = np.zeros([nrow, 1])\n", " y = m*xdata.reshape(nrow,) + b\n", " return y\n", "\n", "# define sum-of-squares function\n", "def ssfun(theta, data):\n", " n = len(data.xdata)\n", " ss = np.zeros([n])\n", " for ii in range(n):\n", " xdata = data.xdata[ii]\n", " ydata = data.ydata[ii]\n", " # eval model\n", " ymodel = modelfun(xdata, theta).reshape(ydata.shape)\n", " # calc sos\n", " ss[ii] = ((ymodel - ydata)**2).sum()\n", " return ss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Data Set - Plot\n", "We consider two simulated data sets. Both are linear, but each one has a different level of observation errors. The first data set has $\\varepsilon_i \\sim N(0, 0.1)$, whereas the second data set has $\\varepsilon_i \\sim N(0, 0.5)$." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEVCAYAAADtr5+DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXxU5fX48c9hICiyKIhikQhudalLNSzpV2sUbN1YrKi4gYqi1g2tVqkbBRW3WuxPRRFFcAEXUEJrRUVTaRsUcMel4gqIioAKAoEk5/fHvZOZDDPJTHLnzp075/165ZWZe587c5KBmZPnOc/ziKpijDHGGBMmLXIdgDHGGGOM1yzBMcYYY0zoWIJjjDHGmNCxBMcYY4wxoWMJjjHGGGNCxxIcY4wxxoSOJTjGGGOMCR1LcIwpMCLyuYhsEJG1IvK9iPxXRM4XkbTeD0Sku4ioiLRsRgwDReQtEflRRL4TkZdFpIcXzy0ivxCROe7j2kJfxhQoS3CMKUz9VbUdsAtwC3AV8KAfTywiuwNTgT8AHYAewD1AjUdPsRl4Ehju0eMZY/KQJTjGFDBV/UFVy4GTgWEi8gsAETlWRN50e1iWisjouMtedb9/LyLrRKRURHZze2FWuT0nj4nItime9kDgM1Wdq461qjpDVb90n7uFiFwtIp+4j/ekiHRM9dxJfqaPVPVBYHHzfjvGmHxmCY4xBlV9HVgGHOoe+gkYCmwLHAtcICKD3HO/dr9vq6ptVbUSEGAc8DNgb6AbMDrF070B7CUifxWRw0WkbcL5i4FBwGHu463B6eFJ9dzGGLMFS3CMMVFfAR0BVLVCVd9V1VpVfQeYhpNwJKWqS1T1RVWtUtWVwJ2p2qvqp0AZ0BVnKOk7EXk4LtE5H7hGVZepahVOojS4OTU/xpjCYwmOMSaqK7AaQER6i8grIrJSRH7ASTq2T3WhiOwoItNFZLmI/Ag82lB7VZ2vqiepamecXqNfA9e4p3cBnnELoL8HPsCpz9nRg5/RGFMgLMExxiAiPXESnH+7hx4HyoFuqtoBuA9nGAog2cykm93j+6lqe+D0uPYNUtUFwEzgF+6hpcDRqrpt3NdWqro8xXMbY8wWLMExpoCJSHsROQ6YDjyqqu+6p9oBq1V1o4j0Ak6Nu2wlUAvsGnesHbAO+EFEugJXNvCch4jIuSKyg3t/L2AAMN9tch9wk4js4p7vLCIDG3juxMcXEdkKKHLvbyUirRv7XRhjwsUSHGMK02wRWYvTW3INTs3MWXHnfw+Mcdtcj1MrA4CqrgduAv7jDiP1Af4MHAT8APwDp0cmle9xEpp3RWQd8DzwDHCbe/4unN6jF9znnw/0buC5E+0CbCA2i2oD8FGjvxFjTKiIqvX4GmOMMSZcrAfHGGOMMaFjCY4xxhhjQscSHGOMMcaEjiU4xhhjjAkdS3CMMcYYEzqW4BhjjDEmdCzBMcYYY0zoWIJjjDHGmNDJ2915t99+e+3evXuuwzCmYC1atOg7d7PMgmbvRcbkVqr3orxNcLp3787ChQtzHYYxBUtEvsh1DEFg70XG5Faq9yIbojLGGGNM6FiCY4wxxpjQsQTHGGOMMaFjCY4xxhhjQscSHGOMMcaEjiU4xhhjjAkdS3CMMQBUVsK4cc53Uzgql1Yybt44KpfaC2/CJW/XwTHGeKeyEvr2hU2boKgI5s6F0tJcR2WyrXJpJX2n9mVTzSaKIkXMHTqX0m72wptwsB4cYwwVFU5yU1PjfK+oyHVExg8Vn1ewqWYTNVrDpppNVHxekeuQjPGMJTjGGMrKnJ6bSMT5XlaW64iMH8q6l1EUKSIiEYoiRZR1L8t1SMZ4xhIcYwylpc6w1Nix4RieEpGjROQjEVkiIlcnOX+5iLwvIu+IyFwR2SXu3DAR+dj9GuZv5P4q7VbK3KFzGXv4WBueMqFjNTjGFLjKSmdIqqwMRo3KdTTNJyIR4B7gSGAZsEBEylX1/bhmbwIlqrpeRC4AbgNOFpGOwA1ACaDAIvfaNf7+FP4p7VZqiY0JJUtwjClgIS0u7gUsUdVPAURkOjAQqEtwVPWVuPbzgdPd278FXlTV1e61LwJHAdN8iNsY4yEbojKmgIW0uLgrsDTu/jL3WCrDgX828VpjTEBZD44xBSxaXBztwSm04mIROR1nOOqwDK8bAYwAKC4uzkJkxpjmsh4cYwpY2IqLXcuBbnH3d3aP1SMi/YBrgAGqWpXJtao6UVVLVLWkc+fOngVujPGOLwmOiFwqIu+JyGIRGZnkvIjI39wZD++IyEF+xGVMIYuuXAxOcXFIkhuABcAeItJDRIqAIUB5fAMR+SVwP05y823cqTnAb0RkOxHZDviNe8wYk2eyPkQlIr8AzsUp/NsEPC8if1fVJXHNjgb2cL96AxPc78aYLEhZXLxsGey8c67DaxZVrRaRi3ASkwjwkKouFpExwEJVLQduB9oCT4kIwJeqOkBVV4vIWJwkCWBMtODYGJNf/KjB2Rt4TVXXA4jIv4Df4UzLjBoITFVVBeaLyLYispOqrvAhPmMKTrLi4tLPHoezz4aJE2Ho0FyH2Cyq+hzwXMKx6+Nu92vg2oeAh7IXnTHGD34MUb0HHCoinUSkDXAM9ce4wWYuGJN18Ztpxq9c3LpVLWd8eA2cdhpUVcG558J//5vrcI0xplmy3oOjqh+IyK3AC8BPwFtATVMey2YuGNM0yYak5s6F/8xZx1kVQ+k09ZlY4912gx13zF2wxphAqlxaScXnFZR1L8uLxSF9mSauqg8CDwKIyM04PTTx0p65AEwEKCkp0awEa0wIJRuSGnXal5Q+OwDefjvW8OijYdo06NAhV6EaYwIoH3ee92sW1Q7u92Kc+pvHE5qUA0Pd2VR9gB+s/sYY7yRupnlcp0ro1at+cjNyJMyebcmNMWYL+bjzvF8L/c0QkU7AZuBCVf1eRM4HUNX7cIoBjwGWAOuBs3yKy5iCEF3vpqICBm98lD0uHu505QC0bAkTJsA55+Q0RmNMcEV3no/24OTDzvN+DVEdmuTYfXG3FbjQj1iMyTfxm2E2Z62a0t61lJZfA7fcEjvYqRPMmAGHZbSQrzGmwER3nrcaHGOMJ9LZDDOtBGjdOjj9dJg1K3Zsn32cIaldd81S9MaYMMm3nectwTEmwJKuVxP3/hKfAEUizjI2Q4cmJDpffAEDBsA779QdWvLzY1h91zR67drerx/FGGN8ZXtRGRNgicXBiZthJiZA99/vJDyVlW6D//yHzQf2rJfcjG/5B/b5uJyyAe1j7YwxJmQswTEmwBrbDDOaADm7DYCqs1bf6NGw5Lop1B5+BK2+XwnAJloxtewhrtA72FwbqesRMsaYMLIhKmMCrrQ0dW1NNAGaOhUmT4bNm4HaGvq9MIrdX7i9rt1KtufEFjPZe69DKXotVtOT2CNkjDFhYT04xuS50lJnlvcrr0D/srU8y/FcSSy5eU9+QZ8WC3i99aEMHdpwj5AxxoSF9eAYk6cSZ0+VdvmMxz8fQBveq2uz+pD+/HTDY5yzoF29WVaW2BgTXPm2JUJQWYJjTB5KnD7++l/m8Yvrf0eb776ra7P8tD/SdcrN9I5E6J1y72xjTJDk45YIQWVDVMbkofjZU6dufIi9L+4L0eSmqAgmT6bro7c606+MMXkjH7dECCpLcIzJQ2VlsFWrGu6UPzBJhxOp2eyc6NwZXn4Zzjwzl+HllIgcJSIficgSEbk6yflfi8gbIlItIoMTztWIyFvuV7l/URvjiG6JEJFI3myJEFQ2RGVMHird90eW/vIUtqt8LnZw//2hvBx22QXwbouHfCIiEeAe4EhgGbBARMpV9f24Zl8CZwJXJHmIDap6YNYDNSaFfNwSIagswTEm33z6KfTvz3bvx31mDxwIjz4KbdsC6W3xEFK9gCWq+imAiEwHBgJ1vyxV/dw9V5uLAI1pTL5tiRBUNkRlTD559VXo1Qvik5tRo2DmzLrkBpJv8RAkIrKN29vita7A0rj7y9xj6dpKRBaKyHwRGeRtaMYYP1kPjjH5YtIkuOACqK527rdu7Rw7/fQtmkZXOA7Kgn4i0gIYApwG9ASqgNYi8h3wD+B+VV2SwxCjdlHV5SKyK/CyiLyrqp8kNhKREcAIgOLiYr9jNMakwXpwjAm66mq47DI499xYcrPDDs7KfkmSG2h8i4cceAXYDRgFdFHVbqq6A3AIMB+4VUSS/zCZWQ50i7u/s3ssLaq63P3+KVAB/DJFu4mqWqKqJZ07d844yMqllYybN47KpcHeDCxf4jQmGevBMSbIfvgBhgyB55+vO/TdzgfwxKnlHEQxDeUtDW3xkAP9VHVz4kFVXQ3MAGaISCsPnmcBsIeI9MBJbIYAp6ZzoYhsB6xX1SoR2R74P+A2D2KqJ1/WOcmXOI1JxXpwjAmqJUugT596yc0zHE/3Zf/mkjuK6+8aHnCqullE9hKRviLSNv6ciBwVbePB81QDFwFzgA+AJ1V1sYiMEZEB7vP1FJFlwInA/SKy2L18b2ChiLyN0+N0S8LsK0/kyzon+RKnCR+veg6tB8eYgKg3rbuqAk44AVavrjt/s1zDtToGpQXUxoqHA9RLk5KIXAJciJN0PCgil6rqLPf0zcDzKS/OkKo+BzyXcOz6uNsLcIauEq/7L7CfV3GkEl3nJNozEtR1TryI07YcSI/9nmK87Dm0BMeYAIif1n1+i4n0rr2QFjWxYuKPRz3EjbeeilSB1kKLFsEoHs7AucDBqrpORLoDT4tId1W9C5CcRuazfFnnpLlx5ssQV66Ti0x/T7mON9uS9RxagmNMAKW72F5FBdRUVfOX2j9wac3fYie6dIFnn2WP3r2Z+xunXadOsGpV3i3g10JV14GzDo2IlOEkObtQYAkO5M86J82J08sPqmwJQhKWye8pCPFmm5c9nJbgGJMlmSy21/fg7ynhZI7khdjBX/4SZs2Cbs6koIAVDWfqGxE5UFXfAnB7co4DHsKHYSHjv3wYimssufCjtyST31Oquqgw9eh42cNpCY4xWZJssb2kCcrHH9Prkv5Q+1Hs2AknwJQpsM02PkWbdUOB6vgDbkHwUBG5PzchmWzK1lCcl0lHQ8lFOr0lXsSSye8pMd5ObTqFskfHqx5OS3CMyZK0Ftt7+WUYPBjWrIkdu+46GD3aKbQJCVVd1sC5//gZi3H40Tvh9VCc10M0DSUX6fTueBVLur+nxHjzYRgwlyzBMSZLoovtRetmotsl1PXiTJgAF1/sdPEAbLUVTJ7srHsTQiKyF86+UNGtE5YD5ar6Qe6iKkz5WsuRjQ/0VMlFY0NHuUouEuMN+jBgLlmCY0wWRZOZerU4c6opfWIk3HNPXbtN2/+Mon/OgpKSHEWaXSJyFXAKMB143T28MzBNRKar6i05C64A5etf/n7W9TQ2dBSEGqN8mZGXK5bgGJNl8bU4barW0OXsk2DJS3XnF3IwQ9bN4pHNXeutTJzuDKw8MRzYN3ExPxG5E1gMWILjoyB8ODeF3x/oDQ0dBSW5SHd4K+zTy5OxBMeYJsgk+YjW4uxS9T9maX96LPlf3bmn5ESG6cNs2tymXhFyJjOw8kQt8DPgi4TjO7nnjI+C8uHcFEGaYh+kWBqSr0OSzWUJjjEZyjT5KC2Fhbe8RI+rTmTrjd/XHf97yWjOeOd6qmukXhFyZaVTY1xVBbV5tmJxA0YCc0XkY2Cpe6wY2B1nawXjs4Y+nL3+a78Qew+CJF+HJJvLEhxjMtTQ9O+kPTv33MM+l19aV0xc03prhtU+zPQ3TyIScTYJHzrUaR9NnqLJTR6uWJyUqj4vInsCvYgVGX8FvK6qNbmLzCTy+q/9Qu09CJJ8HZJsLktwjElDfOKSavp3Ys/Oy3M202fapc5sqaif/Yypx5cz/b6D6yZPFRfHkqFo8hRNbvr1c3pz8rz3BgBVrQXmR++LyBuqelAOQzJJeP3XfkOPZz07/sjnIcnmsATHmEYkG5KKTv+O76mpV0y8cTWtB50Iq1+OPVDPnvDss+z1xc8oeij5+jiJyVNYkpsUCm6Lhnzg9V/7qR7Penb8lS/1Ql6yBMeYRiQbkho1asvEI5qcdN/4IbO0P3usXhI7OWQIPPQQlW9tTUUFjB+ffD+p+LVzQjJ7qiEP5DoAsyWv/9pP9XhB2CbBhJslOMbESVZD09CKxIntF978AsVXnkTb6h/q2vzryBs57PE/UTlf0ipOzvM9pxokIgKcBuyqqmNEpBjooqqvN3Jpps9zFHAXEAEmJa6zIyK/BsYD+wNDVPXpuHPDgGvduzeq6hQvY8sHXv+1n+zxmrtNgjGN8SXBEZHLgHMABd4FzlLVjXHnzwRux1nZFOBuVZ3kR2zGRKWaHZWqV6Ve+1bK4gvuZp+7RjoFNMB6tmZ40SNc8ucTQDLYmyrc7sWZFn4EMAZYC8wAenr1BCISAe4BjgSWAQtEpFxV349r9iVwJnBFwrUdgRuAEpz3q0XutWswnmrONgnpsB4gk/UER0S6ApcA+6jqBhF5EhgCPJzQ9AlVtemiJmcaSkCiiU5lJYwb5yQ60fZSs5m/1l5Mj7/G9oys2mFnpp9UziWn/jKtnqAC0ltVDxKRNwFUdY2IFHn8HL2AJar6KYCITMfZIqIuwVHVz91ziWvw/BZ4UVVXu+dfBI4Cpnkco6Hp2yQ0pHJpJVPfnsrktyZTXVtd1wMEud9125Iuf/k1RNUS2FpENgNtcKaHGhMojSUgiT0848dDl1areLRmMGVaEWvYuzetn3mGs3faqd71BVZfk8pmt4dFAUSkM94v9NeV2Fo74PTi9G7GtV1TtDVZ0tQ6oOjQ1sbqjajzT4xNNZuY+vZUprw9JadDXn7tTp7PvP75s57gqOpyEbkDp0t4A/CCqr6QpOkJ7rj4/4DLVHVpkjbGZE1jCUhiDw/vv8+Sjv3Z6qtPY41OOw0mTXI2zkzxHAWa2ET9DXgG2EFEbgIGE6t3yRsiMgIYAVBcXJzjaMKpKXVA0aGtaHIjCEURp4Mw1wvd+bk7eT7Kxs/fwqPYUhKR7XC6h3vgLNW+jYicntBsNtBdVfcHXgSSFvWJyAgRWSgiC1euXJnNsE2BKi1NPkMKYj08kQgcF/knZ08qrZ/c3HQTPPJIyuTGgKo+BvwRGAesAAap6lMeP81yoFvc/Z2J1fd5cq2qTlTVElUt6dy5c5MDDYvKpZWMmzeOyqWVOY0jOrQVkQhFkSLOO/g85g6dy9ADhtY7nouF7hJjS2d38nzTnH8H2fj5/Rii6gd8pqorAURkJvAr4NFoA1VdFdd+EnBbsgdS1YnARICSkhLNVsCmMGS6mWVpKcx9Sfnp5rvo+88/IJvckZU2beDRR+H447MZbmio6ofAh1l8igXAHiLSAyc5GQKcmua1c4Cb3T/MAH4DjPI+xPAIUs9DQ0NbuV7oLh92J2+O5v47yMbP70eC8yXQR0Ta4AxR9QUWxjcQkZ1UdYV7dwDwgQ9xmQKWrJ4m2bo09WzaROlDv4d/PBg71q0bzJ4NBxzgR9ihJCJnqepkrx5PVatF5CKcZCUCPKSqi0VkDLBQVctFpCfOUNl2QH8R+bOq7quqq0VkLE6SBDAmWnBskgvaPkephraCsNBdPuxO3lTN/XeQjZ/fjxqc10TkaeANoBp4E5gY/2YDXCIiA9zzq3GmbxqTNfH1NFVVcNFFzuzulOvTfPcdnHACvPpq7FifPvDMM9Cli5+hh9GfAc8SHABVfQ54LuHY9XG3F+AMPyW79iHgIS/jCbN873kIkiAkYU3lxb8Dr39+X2ZRqeoNOGtLxIt/sxmFdQMbH8XPmBJxEp2UO3cvXgz9+8Nnn8WOnXEGTJxo9TZpEpF3Up0CdvQzFuOtfO95MN4I4r8DW8nYFKT4GVOdOsHIkSmmh//jH3DKKbB2rXNfBG6+Ga66yrlt0rUjzhoziQvmCfBf/8MxXsrnngfjnaD9O7AExxSs+Cnb++2XUHCsCnfeCVde6dwG2GYbePxxGDBgi8fKtGC5AP0daKuqbyWeEJEK/8MxxoSdJTjGkLA+TVUVXHABTI4rCykudoqJ999/i2tTbfFgYlR1eAPn0p3hZIwxabMEx5h4K1fC734H//537NivfuUUE++wQ9JLbI+p9IjIXjhrYkVXBl4OlKuqzZo0xngu6wv9GZNr0f2jKhtbe+rdd6Fnz/rJzZlnwssvp0xuoP4CgAW8x1SDROQqYDpOzc3r7pcA00Tk6lzGZowJJ+vBMaGW9vDR7Nlw6qmwbp1zXwRuuw3+8IdGi4ltj6m0DAf2VdXN8QdF5E5gMXBLTqIyxoSW9eCYUEs2fFSPqpPIDBwYS27atuWD28oZt/kKKuenN1OqoS0eDOBsqPmzJMd3wvvNNk0DvNhWIShbMxjTEOvBMaHW4A7hGzfCeefB1KmxY9278/bYckpH7MemTc6w09lnw9Chlrw000hgroh8TGy37mJgd+CinEXls1zvFu3FtgpB2prBmIZYgmNCLdnwUWUlvD77G86efTzt3ov7C/SQQ2DmTJ6b1Lmu16emBu6/H6ZMsdlRzaGqz4vInkAv6hcZL1DVmtxF5p8gJAZebKsQtK0ZTPA0lsj7lehbgmNCL34KeGUlXHr4OzxV1Z92fBlrdPbZMGECFBXV9fps3OiMYKna7CgvqGotMD/XceRKthKDTD4svFhO37Zm2FKue+aCpLFE3s9E3xIcU1CW3zuLl6tOoy0/AVArLWhxx+1w2WV1xcTRXp+pU52lcKqrbXZUNnm92WZQZSMxyPTDwovl9NN5jHz4wPcqxiD0zAVJY4m8nz2AluCYwqAKt97KCY/9CcFZmfhH2rH8jifY+/Kjt2ge7fUZOtRmR/nA8802g6ipyUVDH8RN+bDwYjn9hh4jWx/4XiZNXsYY/xpUVVcxumI0o8tGF2yS01gi72cPoCU4Ju9kvC3Cxo1w7rnw6KNE50St6bgrX/ytnANP27fBS+utcGyazDbbdGSaXDT2QRzE4aJs/IXuddLkZYzR16Cquopaannps5eY9+W8gu3JaSyR93NTTktwTF7JZFuEykpYMPtrzio/nnaL40o/DjuM7Z5+mu22396foA3YZptN0tgHcRB3cM5G0uV10tTUGJP1IkVfg9EVo3nps5eo1dqCL75uLJFPPJ+tIU1LcExeSbWuTWKPTmUlXHb4WzxZNYB2dbOSgXPOgXvucbIj4ydfN9sUkaOAu4AIMElVb0k43xqYChwMrAJOVtXPRaQ78AHwkdt0vqqe73V86UrngzhoOzhnI+nyOmlqSowN9SKVditldNlo5n05z7fetHyoc0pHNmuYLMExeSVxXZtOnZL36Hx190zmVp3BNqwH3GLiO/8Cl17a6MrExnt+brYpIhHgHuBIYBmwQETKVfX9uGbDgTWquruIDAFuBU52z32iqgd6GVNTBbGHJh2NJV2Zfjhn4/eQaWIYpN60MBU2Z7Po2BIck1fi17Xp1AlmzHA2/66tdXt0XlFKX76ZEx6/tu6aH2jPV395gr1HHpW7wAuciFwCzFTVZT48XS9giap+6j73dJxNPuMTnIHAaPf208DdIsHMfIPWQ9NcTf1wzvXvIUi9aekkBfnSw5PNOjJLcEzeiQ5D9e0bS25atID2rTZwbsVweHFaXds1HXfji7tnc+Ape+coWuMaC1wtIp8A04CnVHVllp6rK8SPS7IM6J2qjapWi8gPQCf3XA8ReRP4EbhWVedlKc5myZcPsET5ulBgkHrTGksK/O7hac6/xWz+Xi3BMXkpWosTTW5OOnQFE1cOot2Lr8calZU5xcSdOqV8HOObT3HqXfrhDAX9WUQW4SQ7M1V1bS6Di7MCKFbVVSJyMPCsiOyrqj/GNxKREcAIgOLiYt+DzOchiiDO/EpXrnuR4uNoKCnwM4n04t9itn6vluCYwElnGnh8LU7vlot4+P2BtF65PNZgxAi4+25o1cqHiE0a1F3J+AXgBRFpBRwNnALcAXT28LmWA93i7u/sHkvWZpmItAQ6AKtUVYEqN+BFbo/TnsDChB9mIjARoKSkRD2MPS352gsCweoJyWcNzURqKIn0uucvyP8WLcExgZLuNPBoLc6K//c0A2cOJbJyAwC1LSJUnvhXXt3lIsoWiq1hExz16ltUdTNQDpSLSBuPn2sBsIeI9MBJZIYAiYXM5cAwoBIYDLysqioinYHVqlojIrsCe+D0PgVKPveCQHB6QsIiWS9KsiQyGz1/Qf63aAmOCZR0p4GjSumLY2HaDXXXfk8HTtaneOGJI2nxFLRubRtkBsjJqU6o6novn8itqbkImIMzTfwhVV0sImOAhapaDjwIPCIiS4DVOEkQwK+BMSKyGagFzlfV1V7G5wXrBTHxkvWijDp0lC9DV0H+t2gJjgmUtKaBH7gBzjoLnnii7rrVnfbgkDWz+aD250DcrKoKS3CCQFX/l+qciHRR1a89fr7ngOcSjl0fd3sjcGKS62YAM7yMJVuy0QuST4XL+RRrtqXbi5JOu6b8XoPaI2cJjsmJVHU28dPAy8q27NFZWP4VpZcOggULYhf17cuSK5/i8+O3o0XcrKrEDTIz3uLB+OVB4NhcB5HvmvuBH5TC5XR+jqDEGhTp9qI01i5sv1dLcIzvGquzSdz/Kdqj06flQs6bNBC++yp28oIL4K676NWqVb31cVat2nJl43S3eDD+UlVLbprJiw+mIBSLpvtzBCHWoEm3F6WhdmH7vbbIdQCm8KSqs0km2qPz1IlP8qoeSlE0uYlEnC0X7r23bqZUaSmMGuVMoBo1qn4Ck8lzmuwQkb1EpK+ItE04biswNlOyD6ZMRYcvIhLJWbFouj9HEGINo7D9Xq0Hx/gusc4mfhhpC7W1lM4ZA9P/HDu27bbw1FPQr192ntN4zl3J+EKcfZ4eFJFLVXWWe/pm4PmcBZcnGhq68WImSxCKRdP9OYIQaxiF7fcqzrIP+aekpEQXLlzYeEMTSGnVw6xfD2ee6SQzUXvuCbNnO9+z8ZwmbSKySFVL0mz7LlCqquvcDS2fBh5R1btE5E1V/WUWQ80qP96L0hm6CUvRbbZ/jrD8nkxMqvci68ExOZFYZ7OF5cth4EBYtKju0DH4MkQAACAASURBVPc9+zH5qCfps2o7mvK21OhzmmxqoarrANxdu8uAp0VkFxLWyAkjL7e52sAGfnXtrzx7PGOCrDmdMFaDY4Ln9dehZ896yc3XJ1xIt3f/yZU3b0ffvk5vjMkr34hI3Q7dbrJzHLA9sF/OojLGhJYlOCZYpk+Hww6DFSuc+y1bwoQJTD74bjZsbmlFwvlrKFBvrRtVrVbVoTiL6xljjKdsiMpkVWLdS8o6mNpauOEGuPHGukPV7bbjyZOepscBR1CGFQnnM1Vd1sC5//gZSy6k282ejXVImltzMm7eOK575TpqtIaIRBh7+FhGHTqqWTF5LZPfm9Xg+CMIa+pYgmOyJnHtmfHjYeTIJGvR/PQTDBsGM2ILyK7fZS96fT2bDx/enaLHnbbxCwBaLY0JI6/XIfHiQybIew1FZTL7J6ir7oZNEGZkNTnBEZGrVPVWL4Mx4ZK49syMGVuuRVO681IYMADeeit24W9/y329n+DDmzrUa5u4to0xYeN1MuFFwhSED6p05GPiEvbepFy/JmknOCLyZPxd4EAgrQRHRC4DzgEUeBc4y90LJnq+NTAVOBhYBZysqp+nG5sJpsS1Z044AebNi90/rvNr0HMgfPNN7KJLL4U77qB0QUuKbrchKVNYvE4mvEqYcv1BFUZBGMIJu0x6cH5U1XOid0RkQjoXiUhX4BJgH1Xd4CZKQ4CH45oNB9ao6u4iMgQncUq5+7AJrsQam8Rhpf32c+4P3vQ4e1x0NlRVORe2bOmsSnzuuUDya004iMgRqvpy9Huu4wkaL5OJfOl9CbtkPTVh2xYhiDJJcG5KuH9Nhs+ztYhsBtoAXyWcHwiMdm8/DdwtIqL5ugphgUq131O0uHjcOCj7dS2j1l7r3Inq2NEZv0roprF1a0LrDuCguO8mi6z3JbdS9dTkQ21Tvms0wRGRF4ErVPXt+OOqujqdJ1DV5SJyB/AlsAF4QVVfSGjWFVjqtq8WkR+ATsB36TyHyb3KShg92umQqa2Nq7EpjSU+rarWsS9nQO2zsQv33ttZmXi33Wyl4cIT+gX+siHsdRthk6qnxq/etUL+95JOD85VwHgR+Rz4k6quyOQJRGQ7nB6aHsD3wFMicrqqPpppsCIyAhgBUFxcnOnlJkuiCUw0uWnRon7dTEUF7Fj1Jc/UDuBA4vLko4+GadOgQwfb7dt4yt3A8y4gAkxS1VsSzqes+xORUTjD5jXAJao6x8fQG1RIWzaERUM9NdnuXSv0Op9GF/pT1TdU9XDg78DzInKDiGydwXP0Az5T1ZWquhmYCSSuM74c6AYgIi2BDjhvOomxTFTVElUt6dy5cwYhmGyKzpaKJjf9+tVPUI7rVMn82p71k5vLL3d6bjp0qPcYtpCfaS4RiQD3AEcD+wCniMg+Cc3q6v6Av+JOmHDbDQH2BY4C7nUfLxAa2207+oF23SvX0XdqXyqX2pLfuRbtqRl7+FjfEwwvdpnPZ2mtZCzORiofAROAi4GPReSMNJ/jS6CPiLRxH6cvzo7C8cqBYe7twcDLVn+TP6KzpSIRaN3aGaqq63155BH2u7iMHfkWgNqWrWDSJPjLX5wLkjyGzZoyzdQLWKKqn6rqJmA6Ti9yvIHAFPf200Bf9/1pIDBdVatU9TNgift4gRDtDYhIJGndRqF/oAVVabdSRh06yvfek8b+vYRdOjU4/8EZXloMzAfOBD4ELhWRQ1V1REPXq+prIvI08AZQDbwJTBSRMcBCVS0HHgQeEZElwGqcv6BMgDRUH5N0xlNNDVxzDdwat5JAp060mDkTfr3lyvw2a6qgrHO/r83S49fV9LmWAb1TtUmo++uK8z4Xf23XxCfI1XB5Y3UbVrhq4hX6LDpprKNERPYF3k/WoyIiH6jq3tkKriElJSW6cOHCXDx1wcm4PmbtWjj9dCgvjx3bd19nSKpHj6zHa/whIotUtSTXcSQSkcHAUdFlLdze5t6qelFcm/fcNsvc+5/gJEGjgfnRGkEReRD4p6o+ner5gvZeZDU4ptCkei9qtAdHVRc3cPrYZkVl8kKy+piUCc7nnzsrE7/7buzYscfC449D+/bZD9aYuJo+187usWRtliXU/aVzbaDZtHBjHM3aTVxVP/UqEBNcadfH/Oc/0KtX/eTmiitg1ixLboyfFgB7iEgPESnCGfIuT2iTqu6vHBgiIq1FpAewB/C6T3EbE2qVSysZN2+cb8XvttmmaVRD9THR2pwT109h99tGOF08AK1awcSJcOaZ/gdsCppbU3MRMAdnmvhDqro4nbo/t92TwPs4NYMXqmpNTn4QY0IkF1PWLcExaUm2qnBlJRx5RA2jq0axu94eO9G5M++NmcnsFYdQVmkFw8Z/qvoc8FzCsevjbm8ETkxx7U1suXK7MaYZcrE1hSU4JiPxs6n++/yPTN94Gsfx91iD/fbjjdHlHHJ6d1u0z2xBRN4F3on7ehcY5iYVxpiQysUMP0twTD0NTQePn021Z8tP+XenAXQkVoO++pD+dHzuMebc3S79omRTaA4D9ne/hgDTcJagsATHmBDLxZR1S3BMnWTTwSGW8ERnU5XWzGNmze/o+FVsq7Dlp19F14dvgkikrig5+jiNLdpne1AVDncPuwr3CxHZA7g2hyEZY3zi9ww/S3BMncTp4FOnwpQpsURl/Hg4p8VD/K3mfIrY7FxUVAQPPEDXoUPrHieTRftsD6rCIiJ7qur/ovdV9WMR2T+XMRljwrl+kiU4pk5izwvEEp7qqhoOeOSPjNh8Z137n7bZgU/vfIb9hiZuLZa8KDmZjNbYMWFwv4jshrO2zDvAVsB7IrKNqv6U29CMKUyZznDKl2TIEhxTJ7HnBZwenK2qfmAap9D73/+sa/uOHMCgDbP4euQuzN2v6UlJpsNZJj+JiKjjcPd+MXAAcKD7/Q0RQVV/nss4jSlEmcxwyqcdyi3BKUCN7SsVf+zfUz6h+ML+bL8ytj/qR/sM4tAPH+HH2rZEmtnrYntQFYxXRGQGMEtVv1TVL4EvRWQOcCjOonvB2e/AR/ny17AJr0xmOOViundTWYJTYDKqefnXvzjo/N/B6tWxY6NGsfrYG9l8ZAsiHvW6pDucZfLaUcDZwDR3heDvcYanIsALwJ2q+lYO48uJfPpr2KQv35LWTGY45dOGrpbgFJi0a14eeAB+/3uornbut24NkybB6adTivW6mMy4C+vdC9wrIq2A7YENqvp9biPLrXz6aziogpZM5GvSmu4MJy+ne2f7tbMEp8CkqnmpG7Y6pJrSGVfAXXfFLtpxR3j2WejTp+6Q9bqYplLVzcCKXMcRBPn013A6/E42gphM5FvS2pTXLJoMRfeWasrr7cdrZwlOCDVWY5PY+xIdttqq6gdKOBlq58QuOPBAKC+Hbt0wxngrF4ufZUsuko0gJhP5lLQ25zVr7uvtx2tnCU7IxNfYRCJw9tkwdGj9RCex96WiArpVLeHZ2v7szYexE7/7nbMYzjbb+Ba/MYXG78XPsiUXyUYQk4l8Slqb85o19/X247WzBCdk4mtsamrg/vudqd4NFRMPaPsy59UOpiNrYgevuw5Gj4YWLfwI2xiT53KRbPiVTGQ6jJMvSWtzXrPmvt5+vHaiqp4/qB9KSkp04cKCnFXaoGgPzsaNEH1pIxEYOza23UK9oav77oOLL64rJq4tak2LhyfDKafkIHqTT0RkkaqWNPMxdgJWq2qVR2H5zt6LYoJW8OuFZEMxQGh+zua8ZkF5vVO9F1kPTshEa2ymToXJk528pagIOnVKmB4+p5rSJy+Du++OXdylCy1mzYJevWx/KOOXR4DdRGSGql7R3AcTkY7AE0B34HPgJFVdk6TdMGJ7YN2oqlPc4xXATsAG99xvVPXb5sZVKPKl5yITiUMxU9+eypS3pwSqsLk5mvOaBf31tgQnhKI1NkOHbrlRZk0NtKlaQ5ezT4YlL8YuOuggmDULdt7Z9ocyvlHVfiIiwD4ePeTVwFxVvUVErnbvXxXfwE2CbgBKAAUWiUh5XCJ0mqpal0wagvIXfDYlDsUAvtYaFcLvOFsswQmxxGLioiLYpep/zNL+9Fjyv9iJwYOdQp02bQDbH8r4S51x8sUePdxAoMy9PQVn1/KrEtr8FnjR3dkcEXkRZyHCaR7FEAjZ/mAM4hTtbEisFQHq9eA0VHvS3NegUH7H2WIJToEoLYUFt85l16tOZOsNcT32118PN9xQr5jY9ocyXhORtTi9JVucwslx2nv0VDuqanSNna+BHZO06Qosjbu/zD0WNVlEaoAZOMNXeVeo6McHYxCnaGdL4lBMOsWxXrwGhfQ7zgZLcArFvfey72WXON0yAFttBQ8/DCefvEVT2x/KeE1V23n1WCLyEtAlyalrEp5TRSTT5OQ0VV0uIu1wEpwzgKlJYhgBjAAoLi7O8Cmyz48PxiBO0fZLOrUnXrwGhfw79oIlOHkssRA4aWHw5s0wciTce2/swp/9zKm3KUk9AcZWKjbZIiLbAXvg7EUFgKq+mu71qtqvgcf+RkR2UtUV7gytZAXCy4kNYwHsjDOUhaoud7+vFZHHgV4kSXBUdSIwEZxZVOnG7hc/PhjzZb2XXNWwePEa5MvvOKhsmnieSiwEHj/eyWPqFQbvtQZOPNG5E1VS4my70LVr6gc3Jg1NmSYuIucAl+IkFW8BfYBKVT3Co5huB1bFFRl3VNU/JrTpCCwCDnIPvQEcDPwIbKuq37n7ZU0DXlLV+xp6zqC+F1lxau5rWOw18IdNEw+ZxELgGTPq33/nqY8o/Xt/+Pjj2EUnneTMHXeLiY3JgUuBnsB8VT1cRPYCbvbw8W8BnhSR4cAXwEkAIlICnK+q56jqahEZCyxwrxnjHtsGmOMmNxHgJeABD2PLSHM/HIM+hdcPua5hsdcgtyzByVOJhcAnnADz5jn3j468wPAHToJ1P8QuGDMGrr0WRHIWszHARlXdKCKISGtV/VBEfu7Vg6vqKqBvkuMLgXPi7j8EPJTQ5iecnpycy3XPQ1hYDUtqhdC7ZAlODjVnMb1khcD7/UJZO+5ujvznZcg6t5h4662dVf8GD/Y2eGOaZpmIbAs8C7woImtwelpMnFz3PIRFU2tYwv7hXygJtCU4OdLQYnrpJj71CoE3b6Z06sXw9/tjDbp2dYqJDw7EH6XGoKrHuzdHi8grQAfgnzkMKTDiP1St58E7mQ4TFcKHf6Ek0Jbg5EiqxfSalPisXu300LzySuxYz55OcrPTTr79TMY0RkSuT3L4QGCM37EESbIPVZs94/C7N6UQPvwLJYG2BCdHUi2ml3Hi88EH0L8/fPJJ7MFPOQUefNAZnjImWH6Ku70VcBzwQY5iCYxkH6qjDh0Vug/WTOWiN6UQPvwLZfq5JTg5kmoxvYwSnx/nODOjfvwx9sA33gh/+pMVE5tAUtW/xN8XkTuAOTkKJzAK4UO1KXLRm1IoH/6FMMPLEpwcSlxMLzoENX48rFrVQOLTShnyzd/gmMuhttZp0KYNPPII/O53GcVgu4abHGuDsyZOQSuUD9VM5SrxK4QP/0KQ9QTHnQL6RNyhXYHrVXV8XJsyYBbwmXtopqoW1Jh8Yzt4R3t8Xn1pE2cuvIgd74otz1G1YzcePbGcfXY6kEz+S9qu4cZvIvIusT2pIkBnYGzuIgoO+1DdkiV+MWGf2ZUNWU9wVPUjnCJCRCSCs0z6M0mazlPV47IdT1Cls4N36Z6rKB11AvzrX3XH1u7bm/0/eZalE7pQ9GBmSYrtGm5yIP7/eDXwjapW5yoYE3yW+BXGzK5s8HuIqi/wiarauhcJGt3B+/33nWLiTz+NHTvtNCbsOYmlY7ZqUpJiu4Ybv4jI5Q2cQ1Xv9DMeY/JJIczsyga/E5whOPu7JFMqIm8DXwFXqOpi/8LyV7K6l2RFx9F2g4qeY+8/D4G1a53GInDzzXDVVRw6Xyi6pWlJiu0abnwU3U385zhbNZS79/sDr+ckImPyhBWhN41vm22KSBFO8rKvqn6TcK49UKuq60TkGOAuVd0jyWOMAEYAFBcXH/zFF/nXEZRu3UtlJfQ9Qvl91V+5Va8kgltMvM02fHjdYzxTO7DhXcQTHsuSGOO1Jm62+SpwrKqude+3A/6hqr/ORox+COpmmyZc8rkGJ9uxB2GzzaOBNxKTGwBV/THu9nMicq+IbK+q3yW0mwhMBOdNJdsBZ0O6dS+vvrSJu6su4GyN2y6nuJi3x5ZTev4BWyRIqRIXKyQ2AbMjsCnu/ib3mDGmAflai5TL+qEWvjyL4xRSDE+JSBcRZ+EWEenlxrXKx9h8E617iUQaGFJauZILZvSrl9z8uN+v4PXXeW75AVskSA1JllAZk0NTgddFZLSIjAZeAx7OaUTGmKxJVj/kF196cERkG+BI4Ly4Y+cDqOp9wGDgAhGpBjYAQ9SvsTOfNVRrU1YGpe3eg/79af/553XXfHv0UHZ4ZiK0bp1xYbAVEpsgUdWbROR54BD30Fmq+mYuYzLGZE8u64d8q8HxWljGveOHkAZG/s6TLU8hsn6dc1IEbrkFrryy3srEmdbUWA2OyYam1OCEUVjei9KRz3Ugpuma+7oXQg1OQUuVZFRUwKYqZWTtX7it5o+02OQmnG3bwmOPwYABWzxWQzU3yWTa3hivici/VfUQEVlLbKE/AAFUVdt79DwdcRYW7Q58DpykqmuStHse6AP8O379LRHpAUwHOgGLgDNUdVPi9YXI1mIpTF687rmqH/KzBqdgRXtprrvO+V5ZGTt3+K+qmCxncwdX0iL6vt+9O/z3v0mTG2Pykaoe4n5vp6rt477aeZXcuK4G5rqzMOe695O5HTgjyfFbgb+q6u7AGmC4h7HltVzWUpjcyefX3RIcH6Qs9P32W/pc05czah6ONf6//4PXXoP99qs7VFkJ48bVT4yMyUcicqI7NRwRuVZEZorILz18ioHAFPf2FGBQskaqOhdYmxCbAEcATzd2fSGK1lJEJGJrsRSAyqWVjJs3jk5tOuXt625DVD5IWuj7zjtOD038Wj5nnQUTJkDr1nWHbJq3CZnrVPUpETkE6IfTk3If0Nujx99RVVe4t78msynonYDv47aOWAZ09SiuvGf7QhWOxGGp8UeNZ9X6VXn3uluC44MtZk6tLIcjT4WffnIaiMDtt8Pll9crJgbbL8qETo37/Vhgoqr+Q0RuzOQBROQloEuSU9fE31FVFZGszKJIWHQ0G08RSPm6FovJTOKw1Kr1qxh16Khch5UxS3B8UloKpX0UbrsNRo2C6Oy1du1g2jQ49tik19k0bxMyy0XkfpxlI24VkdZkOFSuqv1SnRORb0RkJ1VdISI7Ad9m8NCrgG1FpKXbi7MzzubAyWLI+0VHTTh5MWMpLFtDWILjl40b4bzzYOrU2LEePWD2bNh335SX2X5RJmROAo4C7lDV790k5EoPH78cGAbc4n6fle6Fbo/PKzjrck3P9Hpjcs2rmW5hGY60BMcP33wDxx9fv0r40ENh5kzYfvtGL7dp3iYsVHU9MDPu/gpgReorMnYL8KSIDAe+wEmoEJES4HxVPce9Pw/YC2grIsuA4ao6B7gKmO4Om70JPOhhbMZklZe7jodhONISnCyJrntzTNe3OeDa/rB0aezk8OFw773OmJMxBcSdqXQasKuqjhGRYqCLqnqyo7iqrgL6Jjm+EDgn7v6hKa7/FOjlRSzG+C0sQ0tesQQnC6Izn46ueoZLak8H1jsnWrSAO+6AkSO3KCY2pkDcC9TiTMcegzNVewbQM5dBGRMGYRla8oolOFlQ8Ypy+cZx3Khxkzrat4fp0+Hoo3MXmDG511tVDxKRNwFUdY2IhL4r07Y4MH4Jw9CSVyzBaaKU+ztt3Mg5/zqHzvpY7NDPdmWrF2fDPvv4HaYxQbNZRCK42zWISGecHp3Qsi0OjMkNS3CaIH7xvUgEzj4bhg6F0h5fw6BBdH7ttbq2PxxURocXnoZOnXIYsTGB8TfgGWAHEbkJZ8bStbkNKbu8LPw0xqTPEpwmiF98r6YG7r8fFk16k+dbD6DjT8tiDc89lw53323FxMa4VPUxEVlErBB4MLB/DkPKOiv8NCY3LMHJUGUlfPkltGwJtbXOen2DdCaPVJ/BNtVOMbG2aIHceSdccokVExsDiEh74EKcrQ/KcYqNLwJmA28Dj6W+Or8FpfDT6oBMobEEJwOJQ1MDBygH/P0mRtdcV9fmezrw/NAnGHLpb3MYqTGB8wjO7tyVONO1/wQIMEhV38plYH7IdeGn1QGZQmQJTgbih6a20g3c+Nlw9q2ZVnd+CbtxQqvZ/Gqrvdml0hbnMybOrqq6H4CITMJZ3K9YVTfmNqzCYHVAphBltAdMoYvuC9W1xQoq9DD2fSeW3Pxw0OFMOPN1PmqxNw884PT0xC9cbEyB2xy9oao1wDJLbvwTrQOKSMTqgEzBsAQnhcpKGDeufpJSWgrz71nE+217UqILYifOP58O8+ew/Z4dqa6uv/O3MQaAA0TkR/drLbB/9LaI/Jjr4MIuWgc09vCxNjxlCoYNUSWRchr4sqfY/8JhsGGD0zASgfHj4cILQSTpzt8p18sxpoCoaiTXMRS6XNcBGeO3gktw0kk4tpgGfp/S5YGxlNbcEGvUoQM89RQceWTdocSdvyGWKBUVOecsyTHGGGOyr6ASnPiemYYSjmhPzMaNsJWuZzJncXLNk3XnN3Tbg61fnA0///kW18bv/D1uXCxRig5ZWYJjjDHGZF9B1eDE98w0VCMT7Yn54+lf8aocxsnEkpu59OW+YfOTJjeJoolSJBIbsjLGGGNM9hVUD06yGplUSlsuoPSlgaAr6o7dIxcyqvVfmXNMq7SeL3HIynpvjDHG2KKL/iioBCfthOOJJ+DMM50xKoBIhE8v+3/82PEC5jR0XYrntMTGGH+ISEfgCaA78DlwkqquSdLueaAP8G9VPS7u+MPAYcAP7qEzC2EhQuMfW3TRPwWV4EDqhKOyEiperuWMT0az8+SxsRPbbQdPPcWuffsyyr8wjTFNczUwV1VvEZGr3ftXJWl3O9AGOC/JuStV9eksxmgKmC266J+CS3CSqayE/kf8xP0bh7EzM2In9twTZs92vhtj8sFAoMy9PQWoIEmCo6pzRaQs8bgx2Wabr/rHEhxg0axlzNk4kIN5I3bwN79xhqq23TZ3gRljMrWjal3h3NfAjk14jJtE5HpgLnC1qlYlNhCREcAIgOLi4qbGagpQUDZfLQQFkeA0uPbNa68xYtIgivi67tCKEy9hp8f/4mwZbowJFBF5CeiS5NQ18XdUVUVEM3z4UTiJUREwEaf3Z0xiI1Wd6J6npKQk0+cwBc4WXfRH6D/BG1z7Zto0OOssiqqcP9BqWrTk8yvuZrdbkw3LG2OCQFX7pTonIt+IyE6qukJEdgK+zfCxo70/VSIyGbiiGaEaY3Io9OvgJF37prYWrr0WTj0V3OSGjh2JvPQC3w46b4s9qIwxeaMcGObeHgbMyuRiNylCRAQYBLznaXTGGN+Evgcnce2bI3r/BCcOhZkzY4322gtmz6Zy5e62tYIx+e0W4EkRGQ58AZwEICIlwPmqeo57fx6wF9BWRJYBw1V1DvCYiHQGBHgLOD8HP4MxxgOhT3Di1775zd5LOfgPA+CtuGUtjj7aGarq0IGKp2xrBWPymaquAvomOb4QOCfu/qEprj8ie9EZW+DO+Cn0CQ64a9/IfBg0CL75pu746/83kppRt1Pawfk1ZLLSsTHGmPTZAnfGb1mvwRGRn4vIW3FfP4rIyIQ2IiJ/E5ElIvKOiBzkaRCPPeZkK25yUxtpye9bPcCv5v+Vsn4tueACp+Ym2tszdqwNTxljjJeSLXBnTDZlvQdHVT8CDgQQkQiwHHgmodnRwB7uV29ggvu9eWpr4brr4OabY8c6deLx42cwcfJh1NQ4w1H33w9TpsSSGktsjDHGW7bAnfGb30NUfYFPVPWLhOMDgamqqsB8Edk2OtWzyc+0bh2ccQY8+2zs2D77QHk5u327G0WPOVtNqTpfVnNjjDHZ09ACd1abY7LB7wRnCDAtyfGuwNK4+8vcY01PcCZMqJ/cHHOMU0zcvj2luzm9NVOnwuTJUF1tNTfGGJNtyRa4s9ocky2+rYMjIkXAAOCpZjzGCBFZKCILV65c2XDjyy6Do45ybv/hD1BeDu3b150uLXVyoFdesZobY4zJFavNMdniZw/O0cAbqvpNknPLgW5x93d2j9WT0fLoLVvC9OkwZw6cdFLKZlZzY4wxuWO1OSZb/ExwTiH58BQ4q49eJCLTcYqLf2hW/U1Uhw4NJjfGGGNyyzafNNniS4IjItsARwLnxR07H0BV7wOeA44BlgDrgbP8iMsYY0zu2eaTJht8SXBU9SegU8Kx++JuK3ChH7EYY4wxJvxCv9mmMcYYYwqPJTjGGGOMCZ1QJziVlTBunPPdGGOMSaVyaSXj5o2jcql9YIRFaDfbrKyEvn1jG2faOjfGGGOSscUGwym0PTgVFU5yU1MT24bBGBNuItJRRF4UkY/d79slaXOgiFSKyGJ3c9+T4871EJHX3I1/n3AXKDUhZ4sNhlNoE5yyMqfnJhKxbRiMKSBXA3NVdQ9grns/0XpgqKruCxwFjBeRbd1ztwJ/VdXdgTXAcB9iNjkWXWwwIhFbbDBEQjtEVVrqDEtVVDjJjQ1PGVMQBgJl7u0pQAVwVXwDVf1f3O2vRORboLOI/AAcAZwad/1oYEJWIzY5Z4sNhlNoExywbRiMKUA7xq2C/jWwY0ONRaQXUAR8grNW1/eqWu2ejm76awqALTYYPqFOcIwx4SMiLwFdkpy6Jv6OqqqIpNyzTkR2Ah4BhqlqrYhkEsMIYARAcXFx2tcZY/xjCY4xJq+oar9U50TkGxHZSVVXuAnMtynatQf+AVyjqvPdw6uAbUWkpduLk3TTXzeG6sdQ8wAABlRJREFU9Df+NcbkRGiLjI0xBakcGObeHgbMSmzgzox6Bpiqqk9Hj7tbxrwCDG7oemNMfrAExxgTJrcAR4rIx0A/9z4iUiIik9w2JwG/Bs4UkbfcrwPdc1cBl4vIEpyanAf9Dd8Y4xUbojLGhIaqrgL6Jjm+EDjHvf0o8GiK6z8FemUzRmOMP6wHxxhjjDGhI86wc/4RkZXAF2k03R74LsvhNJfF6A2LsfkyiW8XVe2czWDygb0X+c5i9EbQY2z2e1HeJjjpEpGFqlqS6zgaYjF6w2JsvqDHl8/y4XdrMXrDYmw+L+KzISpjjDHGhI4lOMYYY4wJnUJIcCbmOoA0WIzesBibL+jx5bN8+N1ajN6wGJuv2fGFvgbHGGOMMYWnEHpwjDHGGFNgQpHgiEg3EXlFRN4XkcUicmmSNiIifxORJSLyjogcFMAYT3Nje1dE/isiBwQtxri2PUWkWkQGp2qTyxhFpMxdoXaxiPwrSPGJSAcRmS0ib7ttzvIrPvf5txKR1+Oe/89J2rQWkSfc/y+viUh3P2PMV/Ze5F+McW3tvaiJ8YX+vUhV8/4L2Ak4yL3dDvgfsE9Cm2OAfwIC9AFeC2CMvwK2c28fHcQY3XMR4GXgOWBw0GIEtgXeB4rd+zsELL4/Abe6tzsDq4EiH2MUoK17uxXwGtAnoc3vgfvc20OAJ/x8nfP1y96L/IvRPWfvRc2LL9TvRaHowVHVFar6hnt7LfAB0DWh2UCczfVUnd2DtxVnt+HAxKiq/1XVNe7d+Ti7Gfsmzd8jwMXADFLs1JxNacZ4KjBTVb902/kWZ5rxKdBORARoi/OmUu1jjKqq69y7rdyvxGK8gcAU9/bTQF83XtMAey/yL0aXvRc1L75QvxeFIsGJ53Zf/RInE4zXFVgad38Zyf/DZF0DMcYbjvNXXk6kilFEugLHAxP8j6q+Bn6PewLbiUiFiCwSkaF+xwYNxnc3sDfwFfAucKmq1vocW0RE3sL5YHhRVVP+f1HVauAHnM0nTZrsvcgb9l7UfIX6XhSqzTZFpC1ONj9SVX/MdTzJpBOjiByO86ZyiJ+xxT1/QzGOB65S1dpc/kHfSIwtgYNxNl3cGqgUkfmq+r+AxPdb4C3gCGA34EURmefnv1lVrQEOFJFtgWdE5Beq+p5fzx929l7kDXsvynp8oX4vCk0Pjoi0wnkRH1PVmUmaLAe6xd3f2T3mmzRiRET2ByYBA9XZGdlXacRYAkwXkc+BwcC9IjLIxxDTiXEZMEdVf1LV74BXAd+KJNOI7yycbmtV1SXAZ8BefsUXT1W/B14Bjko4Vff/RURaAh0A3/895iN7L/KGvRf5El+o34tCkeC443EPAh+o6p0pmpUDQ50JDNIH+EFVVwQpRhEpBmYCZ/iZ4cc9f6MxqmoPVe2uqt1xxkN/r6rPBilGYBZwiIi0FJE2QG+c8eegxPclzl90iMiOwM+BT/2Iz33Ozu5fS4jI1sCRwIcJzcqBYe7twcDLqmqLZjXC3ou8Ye9FvsUX6veisAxR/R9wBvCuO5YHTnV4MYCq3odTZX8MsARYj5O5Bi3G63HGFu91u1yr1d/N0NKJMdcajVFVPxCR54F3gFpgko/DL+n8DscCD4vIuzizCK5y/7rzy07AFBGJ4PyR86Sq/l1ExgALVbUc543xERFZglN4OMTH+PKZvRf5F2Ou2XtR82X1vchWMjbGGGNM6IRiiMoYY4wxJp4lOMYYY4wJHUtwjDHGGBM6luAYY4wxJnQswTHGGGNM6FiC8//bu0PcrIIwCqD3WtIEQUKCre8eME1QYNgBS0BWNCElLACLxbIFwipw3USrMB/iZwv9H5k5Rz151c1N5mUGAFiOgQMALMfA4Sza/mx7/e/7ru3XozMB+9FF+1jlJmP+f7dJPrV9mdOrtm8PzgPsSRdtwk3GnE3bX0kukryemYe2l0lukjyfmffHpgN2oYv24IiKs2h7ldO7I39m5iFJZuZ+Zj4cmwzYiS7ah4HDk2v7Ksn3JO+SPLZ9c3AkYEO6aC8GDk+q7bMkP5J8nJnfOb1ee3tsKmA3umg//sHhMG1fJPmc5DrJt5n5cnAkYEO6aE0GDgCwHEdUAMByDBwAYDkGDgCwHAMHAFiOgQMALMfAAQCWY+AAAMsxcACA5Rg4AMBy/gJrrShfQuV8WwAAAABJRU5ErkJggg==\n", "text/plain": [ "