{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Estimating the Error Variance\n", "\n", "Author(s): Paul Miles | Date Created: August 21, 2018\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": [], "source": [ "# import required packages\n", "import numpy as np\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt" ] }, { "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": 2, "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[:,0] = m*xdata.reshape(nrow,) + b\n", " return y\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)\n", " # calc sos\n", " ss[ii] = sum((ymodel[:, 0] - ydata[:, 0])**2)\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": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEVCAYAAADtr5+DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPXVx/HPSSDghguiUoXiVn3cFyrGNQg+VVxwoUq1RcW1aqu1rUoVRamiPq11r2tV1KrgikrdKGlpHVTcatXaokUBdxTFBUKS8/xxZ5LJMJNMkpl7Z+58377ympk7d+aezMi9J7/tmLsjIiIiEidVUQcgIiIiUmhKcERERCR2lOCIiIhI7CjBERERkdhRgiMiIiKxowRHREREYkcJjoiIiMSOEhyRCmNm88zsGzNbYmaLzewZMzvJzPI6H5jZIDNzM+vRjRhGmtnLZvaFmX1iZn82sw0LcWwz28rMnki+rxb6EqlQSnBEKtMB7r4a8G3gEuAs4JYwDmxmmwCTgZ8DqwMbAtcCTQU6xHJgCnBsgd5PRMqQEhyRCubun7v7NOBw4Cgz2wrAzPYzs5eSLSzzzWxC2sv+mrxdbGZfmlmtmW2cbIVZlGw5ucvM1shx2O2A/7r7DA8scff73f3d5LGrzOxsM3sr+X5TzGytXMfO8ju96e63AK9179MRkXKmBEdEcPfngAXA7slNXwFjgDWA/YAfm9lByef2SN6u4e6runsCMGAS8C3gf4ABwIQch3sR2NzMfmdmQ81s1YznfwIcBOyZfL/PCFp4ch1bRGQFSnBEJOU9YC0Ad69391fdvdnd/wHcTZBwZOXuc939KXdf5u4fA5fn2t/d3wbqgPUJupI+MbPb0hKdk4Bz3H2Buy8jSJRGdWfMj4hUHiU4IpKyPvApgJkNMbOZZvaxmX1OkHSsneuFZraumd1jZgvN7Avgzvb2d/fZ7n6Yu/cjaDXaAzgn+fS3gQeTA6AXA28QjM9ZtwC/o4hUCCU4IoKZfZcgwflbctMfgWnAAHdfHbieoBsKINvMpIuT27d29z7AD9P2b5e7Pw88AGyV3DQf2Nfd10j76e3uC3McW0RkBUpwRCqYmfUxs/2Be4A73f3V5FOrAZ+6+1Iz2wk4Iu1lHwPNwEZp21YDvgQ+N7P1gV+2c8zdzOx4M1sn+Xhz4EBgdnKX64GLzOzbyef7mdnIdo6d+f5mZr2BmuTj3mbWq6PPQkTiRQmOSGV6xMyWELSWnEMwZuaYtOdPBi5M7nMewVgZANz9a+Ai4O/JbqSdgQuAHYDPgccIWmRyWUyQ0LxqZl8CjwMPApcln7+SoPXoyeTxZwND2jl2pm8D39A6i+ob4M0OPxERiRVzV4uviIiIxItacERERCR2lOCIiIhI7CjBERERkdhRgiMiIiKxowRHREREYkcJjoiIiMSOEhwRERGJHSU4IiIiEjtKcERERCR2lOCIiIhI7CjBERERkdhRgiMiIiKx0yPqALpq7bXX9kGDBkUdhkhFeuGFFz5x935RxxEVnX9EopPv+adsE5xBgwYxZ86cqMMQqUhm9k7UMeTLzPYBrgSqgZvd/ZKM508CTgGagC+BE9z99fbeU+cfkejke/5RF5WIxJaZVQPXAvsCWwA/MLMtMnb7o7tv7e7bAZcBl4ccpogUgRIcEYmznYC57v62uzcA9wAj03dw9y/SHq4CeIjxiUiRlG0XlYhIHtYH5qc9XgAMydzJzE4BzgBqgL3CCU1EikktOCJS8dz9WnffGDgLODfbPmZ2gpnNMbM5H3/8cbgBikinKcERkThbCAxIe7xBclsu9wAHZXvC3W9098HuPrhfv4qdQCZSNpTgiFSoRAImTQpuY+x5YFMz29DMaoDRwLT0Hcxs07SH+wH/KcSBE/MTTJo1icT8eH/AIqUqlDE4ZnYacDxgwE3ufkXG80YwjXME8DVwtLu/GEZsIpUokYBhw6ChAWpqYMYMqK2NOqrCc/dGMzsVeIJgmvgf3P01M7sQmOPu04BTzWw4sBz4DDiqu8dNzE8wbPIwGpoaqKmuYcaYGdQOiOEHLFLCip7gmNlWBMnNTkAD8LiZPeruc9N22xfYNPkzBPg9WQYCikjXJRJQXw91dcFtQwM0NQW39fXJBOeTT2DttSONs9DcfTowPWPbeWn3Tyv0Mevn1dPQ1ECTN9HQ1ED9vHolOCIhC6OL6n+AZ939a3dvBP4CHJKxz0hgsgdmA2uYWf8QYhOJvUQCfvxjGDoUxo8PWm769g1abqqrg9u63ZvgjDNg221h/vyO31TaVTeojprqGqqtmprqGuoG1UUdkkjFCaOL6p/ARWbWF/iGoBsqcwnQbFM51wfeDyE+kdhKdUUtXQqeXN2loQEWLQq6perrYa+dv2bI5T+EBx8MdthvP/jb36BPn8jiLne1A2qZMWYG9fPqqRtUp9YbkQgUPcFx9zfM7FLgSeAr4GWCJdE7zcxOAE4AGDhwYMFiFImrVFdUKrkxC1pt3n03eDzu2I9YMvQAeP251hdtsgn00BJZ3VU7oFaJjUiEQplF5e63uPuO7r4HwSC+f2fsktdUTk3TFOmcurq2XVEjRwZJzk03wY+H/osvttyZ1dKSm/dGnwFTp8LKK0cXtIhIAYSS4JjZOsnbgQTjb/6Yscs0YIwFdgY+d3d1T4l0U21t0BU1cWLQmrPTTtDYCLs2/YWZy2rp88l/AWiiip/a1dy+zW+DbEhEpMyF1Q59f3IMznLgFHdfnKzgi7tfTzDDYQQwl2Ca+DEhxSUSe7W1baeAj6m6k+ubxlLDcgC+YmWOqLqXp3rtz4y6aGIUESm0UBIcd989y7br0+47cEoYsYhULHdqn/41tcvPa9223nq8dfGj7PzBjpxdF8+1cESkMmkkoUgZSF/DpktJyPLlcOKJcOutrdu23BKmT2ebgQPZhtaVjbt8DBGREqIER6TEdXvV4cWLYdSo4IUpw4bBfffBGmsU5hgiIiVGtahESly2VYfz9s47sOuubZObY46B6dNbkptuH0NEpAQpwREpcZlTvfv2zV4kc4XimXPmwM47w+uvt+40cSLcckvwRu0co66uiL+QiEgI1EUlUuJSU73r64Pk5vTTV+xKyuxievGCR9h8wmj4+uvgTWpq4A9/gCOP7PAYGoMjInGgBEekDKSmek+alL1IZnoX03FLr+E7Z50G3hy8eM01gzIMe+6Z1zFEROJAXVQiZSRXV1JdHfTu2cQV9jOu8p9QlUpuNtwQnnmmw+RGJA4S8xNMmjWJxPxExztL7KkFR6SM5OpKqt32a97d6UjW+utDrTsPGQLTpsE660QRqkioEvMTDJs8jIamBmqqa5gxZoZqgVU4JTgiZWaFrqQPP4QDD2St59IKZh5yCNx5J6y0UujxiUShfl49DU0NNHkTDU0N1M+rV4JT4dRFJVLO3ngjmCmVntz8/OcwdSqJl1fKOttKJI7qBtVRU11DtVVTU11D3aC6qEOSiKkFR6Rc1dfDwQcHC/kBVFXBVVeR2OEUJp8SLFrc2KiF+6Qy1A6oZcaYGdTPq6duUJ1ab0QJjkhZuvNOGDs2KMEAsMoqcM89JPruz7BhsHQpuAdPpc+2Eomz2gG1SmykhbqoRMqJe7BY349+1Jrc9O8Pf/0r7L9/y3TxVHJjpoX7RKQyKcERKUErrEoMQeYydiycl1YNfKutYPZs2GEHYMVp5CeeqO4pEalM6qISKTFZC1/+z2I49FD4859bdxw+PCiYufrqLZu0IrGISEAJjkiJySx8+eKD71B73Ii2NaXGjoXrr4eePVd4vVYkFhFRF5VIyUnvZtq5xxyOv2VI2+Tm17+Gm2/OmtyIiEhACY5IiUl1M919xCP8xfak5tMPgydqauCuu+Ccc4LRw0lZx+uIiFQ4dVGJlKDaOVfDXadDc1rBzIcegj32aNknkYDJk7XejYhINkpwREpJUxP84hdwxRUtmz5bayPmXTud7ffYrGVbaiCy1rsREclOXVQipeLrr2HUqDbJzXM2hC0WJ9h17GZtuqDitt6Nma1iZtVRxyEi8aEER6QUfPhhkJ081FoN/F9bHspeNpMPmtdpaZ1JKff1bsysysyOMLPHzOwj4F/A+2b2upn9n5ltEnWMIlLe1EUlErU33oARI2DevNZtP/85nx18Gc17V1HdsGLrTAzWu5kJPA2MA/7p7s0AZrYWMBS41MwedPc7u3sgM9sHuBKoBm5290synj8DOA5oBD4Gxrr7O909rohESwmOSIReu2YmG/3yEFZamlYw8+qr4eSTqaX9JKbM17sZ7u7LMze6+6fA/cD9ZtbtefDJbq9rgb2BBcDzZjbN3dPm3fMSMNjdvzazHwOXAYd399giEi0lOCIR+c95k9l04nHUEFznm1Zaheqp98J++7XsU+ZJTE7uvtzMdgru+vNmtgWwD/Avd5+e2qcAh9oJmOvubwOY2T3ASKAlwXH3mWn7zwZ+WIDjikjElOCIhC1ZMHPTiee3bHqP/jw+9lHG7rdDhIGFx8zOB/YFepjZU8AQgm6rs81se3e/qECHWh+Yn/Z4QfJYuRwL/KlAxxaRCCnBEQlTQwOccALcfnvLplfZikN7PcbtRw6MMLDQjQK2A3oBHwAbuPsXZvYb4FmgUAlO3szsh8BgYM8cz58AnAAwcGBFfVciZUkJjkhYFi+GQw6Bma09Iou/uzdP7TuV2/dZPZZdUe1odPcm4Gsze8vdvwBw92/MrLmAx1kIDEh7vEFyWxtmNhw4B9jT3ZdleyN3vxG4EWDw4MFewBhFpAg0TVykADoslzBvHuyyS5vkhrFjWePvj3HGBUE18Aort9BgZisn7++Y2mhmqwOFTHCeBzY1sw3NrAYYDUxL38HMtgduAA50948KeGwRiZBacES6KbWqcENDjnIJzz8P++8PH6VdO3/9a/jVr8Cs49fH0x6plpLUFPGknsBRhTqIuzea2anAEwTTxP/g7q+Z2YXAHHefBvwfsCow1YIaX++6+4GFikFEoqEER6SbUqsKNzVlKZfw8MPwgx/AN98Ej2tq4Lbbgm35vD6m2ukG+gT4pMDHmg5Mz9h2Xtr94YU8noiUBiU4It2UWlW4IXNBviuvhJ/9rLWewlprBSsV7747iUTr+jY5Xx9zZrY5wZTt9ZObFgLT3P2N6KISkbhQgiPSTSusKrxTE+8fdgb9p17VutNGG8H06bDZZlm7pFKv79u3tSRDnFtxzOws4AfAPcBzyc0bAHeb2T2Zqw2LiHSWEhyRAmhZkO+rr/h06JH0n/Vwy3NLttyZ1WZOg379gOxdUuPGBftW0FicY4EtMxfzM7PLgdcAJThAYn6C+nn11A2qo3ZAfP9nECkGJTgihfLBB3DAAaw1Z07LpvsYxX8Pm8wv+63Usi1Xl1SFjcVpBr4FZNZ86k9hZ1GVrcT8BMMmD6OhqYGa6hpmjJkRWpKjxEriIJQEx8x+RlDMzoFXgWPcfWna80cTzGRIrU9xjbvfHEZsIp2RPnamTfLx+utBwcx3Wq/Xv7VfcF6vS3l676oVXpetxlSFjcU5HZhhZv+hdaXhgcAmwKmRRVVC6ufV09DUQJM30dDUQP28+lCSjSgTK5FCKnqCY2brAz8Ftkgu4jWFYC2K2zJ2vdfddWKTkpVzOvfMmXDwwfD558GOVVW8fcY1NKz1Y56uCzalXlddDWPHwpgxrd1SKTGoEJ43d3/czL5DUCsqNcj4PeC55AKAsdPZVpG6QXXUVNe0JBp1g+qKHyTRJVYihRZWF1UPYCUzWw6sTHAiEykrWbuQ/jMZjjsOlieHkqyyCkyZwkYjRpDKXyZNan1dUxPccENQqSHbGJu4FtfMJrn+zezUYzN70d1jWYyrK60itQNqmTFmRuhdRVElViKFVvQEx90XJuvLvAt8Azzp7k9m2fVQM9sD+DfwM3efn2Ufkci06ULq6fzwrQvhVxNad+jfHx57DLbfPuvrli4NZoy7V8QYm66wqAMolq62itQOqA299SSqxEpKV7mOyQqji2pNgrUuNgQWE6wW+kN3vzNtt0eAu919mZmdCNwO7JXlvVTsTiKT6kL669MNHJM4nnVumdz65NZbB8nNgAFtXpMae3PFFfDSS3DrrdDYWBFjbLripqgDKJZyaxWJIrGS0lTOY7LC6KIaDvzX3T8GMLMHgF2AlgTH3Rel7X8zcFm2N1KxO4lKKlEZtuNizprZtmAm//u/MHUq9Omzwmsyx+yMGVMZY2zyZUFthCOBjdz9QjMbCKzn7s918NKyolaRcJVri0MpKucxWWEkOO8COycL630DDAPmpO9gZv3d/f3kwwMBrWQqJSOVqPRfNo+DfASkL7R73HFw3XXQs+cKr8u13o0SmzauI5gWvhdwIbAEuB/4bpRBFYNaRcKRrcUBUMLTReXW+pgujDE4z5rZfcCLQCPwEnBjRrG7n5rZgcnnPwWOLnZcUnlyTvHuQH09bLvsOR5qPoB1SSuYefHFcPbZYNmHjlTYtO+uGuLuO5jZSwDu/lmy6rdIl2S2OEx+ZTK3v3J75F0s5dqqVM6tj6HMonL384HzMzanF7sbBy2TTkQKrjsVuw+2hzit+QhWJiiY2dyzhrnn3s79jKZudu73qaRp392w3MyqCdbIwsz6oYX+pBsyWxyAyLtYynkcC5Rv66NWMpaK0OVVgq+8ks1/9TOS11+W91mLNy95iJ1+vnteyVIlTfvuoquAB4F1zOwiYBRwbrQhSbrutjyE3XKR2eIAtGnBiaKLpZzHsZQzJThSETrdXdTUBGecAVelFczceGN6Tp/OI/d/p5JKKhSVu99lZi8QjM0z4CBVEy8d3W15iKrlIrPFIeoulnIex1LOlOBI7KVP1V60KI/uoq++giOOgGnTWrfV1sLDD0O/fhpbU2Du/i/gX1HHISvqbstDqbRcRN3FUs7jWMqZEhyJtU6PvfngA74cegCr/ittot+oUTB5MqwUFMzU2JruM7MhwBvu/oWZrQScDewAvA5c7O6fRxqgAN1veVDLRauok6xKVBV1ACLFlG3sTUoiEZRRSCSSG157jaXb79wmuVl45Jlw770tyU1KbW1rLak27yH5+gPwdfL+lcDqwKXJbbdGFZS0lWp5mDh0Ype6l7r7epHuUAuOxFqu7qTMlp3nL53BluMPpXeyYGYTVfzErmXAlicxLsefAd2ZmSVUuXtj8v7gtBpUfzOzl6MKqtLkMwC4uy0ParmQqCjBkVjL1p2USMCECbBsGTQ3w+ilt7H56cdDc3C9XcKq/KBqCn/utS8z6nK/d67WIXVd5eWfZnaMu98KvGJmg919TrLC+PKog6sE5T51uVyU6/o3caAER2Ivfap2qtUlSG6cC5jAeX5hahY4S/p8iyd/8ii7rrI959S1n6Rktg717asWnU44DrjSzM4FPgESZjYfmJ98ToqsVAYAx5mSyGgpwZGKkmp16dG8jJs5nh9xR8tz/7BtOODLx/j48g3ySk4yW4e6vNZOBUoOIj7azPoQFOLtASxw9w+jjaxyaABw8SmJjJYSHKkodXWwTs/P+GPTIdRR37L97U2/x55vTWFxcx+qO5GcZC7kp+njnePuXwCvpG9L67qSIop66nIldN0oiYyWEhypKLXr/Ze5645g5XfSll05/ng++tG1LPteT6q7kZxo+njBXIBmUoUiqgHAUXfdhJVcRZ1EVjolOFI5nnsODjiAlT9KK5h5ySVw5pnsbFaQ5ESlGfJjZv/I9RSwbpixSH4KmRRE2XUTdnKlWWTRUYIjleHBB+HII+GboGAmvXrBbbfB6NEtuyg5CdW6wPeAzzK2G/BM+OFIewqdFETZdaNxMZVDCY7EmztceWVQV8qTU6XWWisou7DbbkBrKQd1K4XqUWBVd19hzRszqw8/HGlPoZOCKLtuNC6mcijBkdhYIVFpaoKf/Qyuvrp1p403hunT4TvfaXmNpnaHz92Pbee5I8KMRTpWjKQgqq4bjYupHEpwJBYyE5WZj37FkCszCmbusgs89BD069eySVO7RTpWiKSgq2N4ijEgWONiKoMSHCkbubqSMlcmXmvZ+wwccwAsfKF1p8MOC8bcZNSUUmXw0mBmuwE7Af909ycL/N77ENS7qgZudvdLMp7fA7gC2AYY7e73FfL4paQ7yUJ3koKujuGJerZVehxq8Sk/SnCkLKS30FRXw9ixMGZM8FzrysSwlb3GYz6C/gvfbX3xWWfBxRdD1YpFpTS1Oxpm9py775S8fzxwCvAgcL6Z7ZCZhHTjONXAtcDewALgeTOb5u6vp+32LnA08ItCHLNURZksdHUMTykMCC6VJKsUlFuipwRHykJ6V1JTE9xwA9x+Oxx1VLC9uRmG2wweqj6EVRq/CF5UXQ3XXQcnnNDue2v2VCR6pt0/Adjb3T82s98As4GCJDgErUJz3f1tADO7BxgJtCQ47j4v+VxzgY5ZkqJMFro6hqe7Y38KcUEuhSSrFJRjoqcER0peIgHvvgs9egSJjHvw09AQPF9TExTMvMGPp2djskD1qqvClCmw777dOq5adoqmyszWBKoAc/ePAdz9KzNrbP+lnbI+QX2rlAXAkAK+f9mIcvZQV8fwdGfsT6EuyJp1Fehqohdlq48SHClpmV1TI0fCn/4EjY1BYjPmR845y85ng1sntr5o/fXhscdIfL0t9ZO6lqBodlXRrQ68QLDujZtZf3d/38xWTW4rOWZ2AkFrEwMHDow4ms6LevZQV8fwdPV1hWp5ifpzKxVdSfTSk8zqqmrGbjeWMduOCe0zVIIjJS29awpgp53gzDOD7UN3WcbOvz8O7ryz9QXbbBMkN/M36FaCotlVxeXug3I81QwcXMBDLQQGpD3eILmt09z9RuBGgMGDB3v3QwtfJc0eyueC3FHrQvrz43YfF0LUpasriV56ktnU1MQNL9zA7a/cHlr3lhIcKWnZZjnV1kLt5p/x+fBD4MX61p2/972gW6pPH+rv6F6CotlVxWVmPwUecPcF6dvd/WvgvwU81PPApma2IUFiMxrQOjsVoKMLckddWOU45qTYOpsgp5LMpY1L8eR/YY5jWnFaiUg3JBIwaVJwWwipWU4TJ6a1wvz3v3y9/S6snpbcfHjg8fDII9CnD9CaoFRXdy1ByXpcKaSJwHNmNsvMTjazfh2+ogvcvRE4FXgCeAOY4u6vmdmFZnYggJl918wWAN8HbjCz14oRi4SvdkAt43Yfl/Vimq0LqzPPl5LE/ASTZk0iMb9AJ94CSSWZJ+54Ir2qe1Ft1aGOY1ILjhRMscattJnllKVg5ji7hD5DzmRcT2vzmu5O/9bsqqJ6G9gRGA4cDlxgZi8AdxO07Cwp1IHcfTowPWPbeWn3nyfoupIK0lEXVrkMLi71lqZUq8+YbceEPo5JCY4UTNHHrTz4IBxxBCxdCsBSenGM3c7DvQ9nxtAVd1eCUtLc3ZuBJ4EnzawnsC/wA+A3QFFadERSOurCKpfBxeUyjT2K8V9KcKRg2hu3kjnlOvW4b19YtKiDVhZ3uOIK+PnP2xTMnHvRw2zz2W78tL3XSqlqM1PK3ZcD04BpZrZyNCFJpenoolsOg7LLpaUpCkpwpGBydQtldl1dcQWcfnrr6sNVVdCrV44uraamYOdrrmndtskmMH06W226KVuF9LtJwR2e64nkQGMRyUO5tDRFQQmOFFS2bqHMrqv7729dfRiC26xdWl9+GXRJPfJI67ZddoGHH4a1184ZgxboK33u/u9cz5nZeu7+QZjxRKnclr+vZKX6XZVDS1MUlOBI0WV2XR16KMya1bYFZ4WZTu+/D/vvDy++2LrtsMOC+gy9e+c8lhboi4VbgP2iDiIMpT5AVFrpuyqsMJJFJThSFJmtKJldV1tv3c4YnNdegxEjgvoMKe0UzEynBfrKn7tXRHID5TNAVKL/rkq19agrwkoWleBIweVqRUlPNHLOcJoxAw45BL7oXMHMFC3QV77MbLK7j4k6jjDFcYBonC7E6aL8ruLWehRWsqgERwquy60ot90Gxx8fFJqCLhXMLMT6N1J8ZjYtcxMw1MzWAHD3A8OPKnyFHiAadXIRtwtxuigH80bdelRoYSWLSnCk4NJbUaqrg56mRKKDaeDnnx8sG5ySLJjJttt2+vha/6YsbAC8DtwMOEGCMxj4bZRBRaFQA0RLIbkoxIU46iStPVEN5g279ajY30FYyaISHCm4VCvK5Mlw661w003B2OCsA36XLYNjj4W77mrZ9OF627Dg2sfYcVstLhtjg4HTgHOAX7r7y2b2jbv/JeK4ylYp/JXf3QtxKSRppSjM1qOwvoMwksVQalGZ2c/M7DUz+6eZ3W1mvTOe72Vm95rZXDN71swGhRGXFE9tLQwcGPQ2pXdVtfHpp/C//9smuXmq6nts/tEsdv/BBgWrZyWlx92b3f13wDHAOWZ2DfqDq1tSyUXY9X7SpS7EE4dO7NKFsZzqP4WtvbpahRSn76DoJxQzWx/4KbCFu39jZlMIKvrelrbbscBn7r6JmY0GLqWdhcCkPLQ74Pftt4OZUm++2bLppe+ewAEvXMOy5p5UawZURUhWE/++me0HfBF1POWsVBZ8685f5nEcdF1uuvIdlGq3Ylh/MfUAVjKz5cDKwHsZz48EJiTv3wdcY2bmnlqXX8pRzgG/zz4LBxwAH3/cuvMll7B09zOpGm5UawZUxXH3x4DHoo6j3JX7gm+lkqR1Rale5Durs99BKXcrFj3BcfeFZvYb4F3gG+BJd38yY7f1gfnJ/RvN7HOgL/BJseOT4lphwO8DD8CRR7YUzKRXL5g8mcSAw6ivD8o4dFibSkRiqxyTtM5c5MshEerMd1AKY79y6XKCY2Znufuleey3JkELzYbAYmCqmf3Q3e/swjFPAE4AGDhwYGdfLlHKVjCzb194+GESVbtq9WERKUuJ+Qkm1E9gWdMymr253Yt8eiJUXVXN2O3GMmbbMSWTEHRFKXcr5j3I2MympP1MBY7L86XDgf+6+8fJisEPALtk7LMQGJA8Tg9gdWBR5hu5+43uPtjdB/fr1y/f0CVqTU3w05/CGWe0JjebbBLMHd9116zr5oiIlLpUwvL020/T7M1UWVW7F/mwyB7IAAAgAElEQVTM1o4bXriBYZOHkZhfvjMqujuwvJg6M4vqC3c/LPnzfeDpPF/3LrCzma1sZgYMA97I2GcacFTy/ijgzxp/ExNffgkHHdS2GviuuwbJzaabAq2DkaurNfamEpnZXum3IuUilbA000wVVQzfcHi7F/lUa4dhADhe9jOVILwZXp3VmQTnoozH5+TzInd/lmDg8IvAq8lj3mhmF5pZarXSW4C+ZjYXOAM4uxNxSal67z3Yc0949NHWbYcfDk8/3aYaeGow8sSJ6p6qUL/JuJUCS8xPMGnWpLJuKShF6VPze/XoxYS6Ce1e5FOtHSfueCK9qntFOqW/ElhHDSVm9hTwC3d/JZyQ8jN48GCfM2dO1GFILq++CvvtB/Pnt247+2y46KIOC2amZBbslNJhZi+4++ACvdeL7r6Dmb3k7tsX4j2LrZzOP6U8yyUOujpouBwGGxdDIX7vfM8/+QwyPgu4wszmAb9y9/e7FJFUjqefhkMP7XLBTMhdsFNEOqeUZ7nEQVdnfeX7ujglQmEn2x3+Ke3uL7r7UOBR4HEzO9/MVipaRFKSEgmYNImOVxf+wx+C4pip5Ga11YKaUp1IbiB7wU6ROCtWN1IprHAsXZNKCMbPHF/2g5Eh/FWS85omnhwc/Cbwe+DXwPFmNs7d7yhmcFIa8mpNcYfx44MuqJQsBTPz7XZqdxVkkZgp5l+25bx4XqWLW+tb2FPKO0xwzOzvBGvYvAbMBo4G/gWcZma7u3vn/jSXspOtNaVNcrJsGYwdC3/8Y+u2bbcNkpv112/Z1Jlup5yrIEscfZm8XRJpFBEq9oWsHBfPk9JeY6Yrwk6282nBOQF4Pcu07Z+YWeZ0b4mh9lpTnn/iU9Y96WAGzvtr68Z99oEpU4LuqTQdJkoZVlgFWWLJ3fdIv61EcbuQSWHEsfUtzGS7wwTH3V9r5+n9ChiLlKhcrSkv3vc2fQ4bwUBvLZjJSSfB1VdDjxX/11K3k0h2HV3I4jTQVDpHrW9d161aVO7+dqECkdK2QmvK7NlsdvSBrOKtBTP/vM+l7HXdL8Es53uo20kku1wXMk3zFumaziz0JxK4/34YOpRVvgqSm6X04siaKax03pk5k5uU2loYN07JjUi+wp55UmyFnC2mBQylPUWvJi6lq9ML6bnD5ZfDL3/ZUlNq+RprM+Xwhzn1qF2UtIgUQT7jc8qlC6uQrVFq2ZKOKMGpUB3NaFoh+WlshNNOCxbsS/q076a88/vpjPn+JmGHLzFhZq8C/0j7eRU4yt0zS8NUrHzG55TLhb6Qs8XiNoVaCk9dVBWqvYX0UsnP+PHB7bMzkgUz05Kbv1ftxuafJdj1qE06XvxPJLc9gZuAb4DRwD+BEYU8gJntY2ZvmtlcM1uhzp2Z9TKze5PPP2tmgwp5/EJor5hhOXVhFXLRQS1gKB1y97L82XHHHV3y88wz7hdfHNymb1tpJffq6uA2/bmLLw62g/sGVQv9/W9tHzxI/ry29eG+ctU3DsF+F18c/u8k0QLmeBH+XQObArcX8P2qgbeAjYAa4BVgi4x9TgauT94fDdzb0fvmc/4B9KMf/XThp1DnH3VRxVyurqj2ZjSlpnNvuuyfPOojWO+9tIKZ48bx+X6/xveuolrTvaWbzOw77v7v1GN3/4+ZbVPAQ+wEzPXkjE8zuwcYCbyets9IYELy/n3ANWZmyROpiJQpJTgxl94VtWwZTJgQ/KSSnGwDg2trYc6kp9jorFH0XpZWMPP66+G446hlxeRIlb+li24ws42BhQRjcHoD/zSzVdz9qwK8//pAWobOAmBIrn3cvdHMPgf6Ap+k72RmJxAsfMrAgQMLEJqIFJMSnJhLtcYsWwbNzUGh71mzOqjOfcstbPGLk4KBxRCsSDx1Knzvey27pCdHqvwtnZVqIfGgkC9mNhDYFtguefuimeHum0UZZzp3vxG4EWDw4MEdtu7EsQGonAY0F9OkWZMYP3M8Td5EtVUzcehExu0+LuqwgPy/o1L5HYo5A1CDjGMu1RU1fDhUVQVJTs7q3O5w7rlw3HGtyc0GG8Df/tYmucmkyt/SBTPN7CfJxAZ3f9fdHwEuJSjq+yxwbQGOsxAYkPZ4g+S2rPuYWQ9gdWBRAY4duUKvE1NOA5qLqZgDnLv7neX7HZXCIO1iV0tXC04FqK0NuqVmzWqnTMKyZXDMMXD33a3bttsOHn20TcHMbFSCQbpgH2AscLeZbQgsJuieqgaeBC5395cLcJzngU2Tx1hIMIj4iIx9pgFHAQlgFPDnOIy/KUZri2pmBYpVI6oQ31m+31Ep1Lkq9lR/JTgVot0yCYsWwcEHBxlQyr77wr33rlAwM9tYG5VgkM5y96XAdcB1ZtYTWBv4xt0XF/g4jWZ2KvAEQfL0B3d/zcwuJJiJMQ24BbjDzOYCnxIkQWWvGBePUrgolopi1IgqxHfWme8o6jpXxU6YrVz/UBk8eLDPmTMn6jDK31tvwYgR8O9/t2478US45hro0aNNQgMaayMBM3vB3QdHHUdUyuH8o/Ey5acSv7OujMHJ9/yjFpxKNns2HHAAfJI2WeSyy+AXvwCzFQYPH3XUimNtlOCIlCa1tpSfSvzOitmKpASnAiUS8P7V93HQ/T+iqmFpsLF3b7jjDhg1qmW/zMHDoLE2IuUkjC6IcqmDVS6i7jaKEyU4FSbxjDNtz98yqfGXrRvXXhumTVuhOSZz8PCYMcGPxtqICFRml4qUDyU4MZRz0b3GRnqd8VMmNf6+ZdOitb9D39nTYeONV3ifXIOHldhIsZhZf+BTd18WdSzSsfRBscsalzGhfgIT6iYoyZGSoAQnZnIuuvfllzB6NDs8+1jLvn+v2o2aOx6i78Z9c75frtWORYrkDmBjM7vf3X8RdTDSvtQsmGWNy2immaf/+zSz3p2llhwpCVroL2ayLrr33nuwxx7wWGty8/o2o6me8RTf3Sd3ciMSNncfTlAY89aoY4lSoRfoK5bUoNjhGw2nyqpo9uZQFgAsl89HoqUWnBLWlfpOmeNmRgx4FYaMgAULWncaN44tfv3rYGljkRKTXGTvtajjiEohxrWEOfC3dkAtE+omMOvdWaEsABjXcT8arF14SnBKVHpXU3U1jB0L228frMnXXsKTPm5m5EpPssXJo2DJkuDJtIKZIlEysyVAtkW4jCDH6RNySCWju4u9RZEAhDm9udir30ahM99ZmIlQuSddSnBKVHpXU1NTkJdA0OjSq1f7i+zV1kLtazfDSScFL4asBTNTVAlcwubuq3W8V2XqzOqu2S5AUSUAYU1vLubqt1Fd0PP9zsJMXuPQUqYEp0SlupqWLg1qYKakF8vMmow0N8P48XDxxa3bNtggGH+zzTYr7K5K4BI1M1sT2JSgFhUA7v7X6CKKVr6tIbkuQHGvF1XKdaC6Kt/vLMzkNQ4tZUpwSlSqq2nyZLj1Vli+PMhdqqraWWQvW8HM7beHRx7JWTAz26BkJTgSFjM7DjiNoMr3y8DOBEUv94oyrqjl0xqS6wJUCavhlmodqK7K9zsLM3mNQ6KsBKeEpaZopxbX69u3nTE42QpmjhgRFMxcddWcx1AlcInYacB3gdnuPtTMNgcu7uA1QvsXIK2G23lRX9Dz+c7CTF7jkCir2GYcvPVWUP37P/9p3XbSSXD11dCj4xxWY3CkswpVbNPMnnf375rZy8AQd19mZq+5+5YFCLNoSuX8U+6DQEuNPs/uC+MzVLHNSpFIwIEHtimYOe/U/2PQVT8Hs7zeIn0xPyU7ErIFZrYG8BDwlJl9BrwTcUxlQy01haXPs3tKbWCyFkIpZ/fdB0OHtiQ339Cbw6qmssUtvyAxO7/kJl1qwPH48cFtQmtoSZG5+8HuvtjdJwDjgVuAkdFGJdJ9lbgYYbbSHVH+/mrBKQGdbjVxh9/8Bs48s2XT4p5rM6LxERLNO1PdxcHCGnAsYTOz87Js3g64MOxYpHIVulul1FoywlJqpTuK3oJjZpuZ2ctpP1+Y2ekZ+9SZ2edp+2Q76cVSp1tNGhv54NCT2yQ3/7bv8N3G2SR85zazrBIJmDQp/5aY1IDj6moNOJbQfJX20wTsCwyKMiCpLKlkZPzM8QybPKwgLQ7ZZmRVgqhKd+RS9BYcd3+T4C8yzKwaWAg8mGXXWe6+f7HjKTWdajVZsoTPvjea9RLTWzbNXW83dvvoIT5u7ktVFQwfDhMmBM91dn2bXNXDRYrF3X+b/tjMfgM8EVE4oYjLQNZi/x5hfU7FmB4e9YysjhTzsw27dEd7wu6iGga85e4aRJiU9zTthQth//1Z8+WXWzbdbT/gmf3/wJd39aY6+foJE4LEZNKkrnU3qXq4RGxlgjVxYqlUui4yL3CdveAV+/cI83MqRjJSylOsw/hsS+X3DzvBGQ3cneO5WjN7BXgP+IW7V0SxvbxaTf7xj2BNm4ULWzZdbOdwUa8LeXpsFUeMXfH1Wt9GyoGZvUprTapqoB8wMbqIiqsUVofNvMBdsc8VnP746Z264BX79wjzcyrWxbhUZ2SF9dmWwu8fWoJjZjXAgcC4LE+/CHzb3b80sxEEU0Y3zfIeJwAnAAwcOLCI0Yar3VaTJ5+EUW0LZr515g3Yasfyu76tic24jE9V3U1SJtK7pRuBD929Mapgiq0Uui4yL3D3v35/py94xf49wv6cSuFiHJZsn21cuk0zhbbQn5mNBE5x9//NY995wGB3/yTXPqWy0FZndWrG1E03wY9/3Fows0+fYGr43nurhpREqrsL/ZnZGe097+6Xd/W9w9Cd80/UF5NCtOCk3icOY3CiFsXvmX5MoCS7TdtTigv9/YAc3VNmth7BX25uZjsRzO5aFGJsocg7KWluhnPPDQbSpAwYEBTM3HprQFO6peylqolvRlCqYVry8QHAc5FEFJKoWwuydclsvc7Wnb7IFvv3iPpzCkNYY40yk4f0z3bSrEkl121aqM8hlATHzFYB9gZOTNt2EoC7Xw+MAn5sZo3AN8BoL9caEu3IKylZujQomHnPPa3bdtghKJj5rW+1bNIYGyln7n4BgJn9FdjB3ZckH08AHoswtIqQmTxUQjJRisIYD9NR8lCK3aaF+hxCSXDc/Sugb8a269PuXwNcE0YsUeowKVm0CA46CP72t9ZtOQpmaoyNxMS6QEPa44bkNpHYCyO56Ch56Mog60J3qxXrc9BKxiFqNymZOzdIZtILZv74x3DVVTkLZmpKt8TAZOA5M0utjXUQcFt04RRepYwlKQel9l2EMZ06n+ShMy146S1C1VXVjN1uLGO2HdOt2Iv1OaiaeCl45hkYObK1YKYZ/N//wRln5F0wUyRMhaomnnyvHYHdkg//6u4vFeJ9iynf80+prHsjlf1dFDKxmzRrEuNnjqfJg8kvhtG7R+9QP89SHGQs2UydCj/6ESxbFjzu3RvuvBMOPTTauERC4u4vAC9EHUcxlMK6NxKo5O+ikGOsUi1CSxuX4sn/SvXzVDXxqLjDZZfBYYe1Jjf9+sHMmUpuJPbM7G/J2yXJ+nSpnyVm9kWBjrGWmT1lZv9J3q6ZY7/HzWyxmT1aiOOmS10Mqq26JJfsryT6Lgoj1Z104o4n0qu6V0l/nuqiCkmb9W++2winngo33NC6w2abwfTpsNFGUYUokrdCdlEVi5ldBnzq7peY2dnAmu5+Vpb9hhGUiDgx33p4nTn/lNq4j3LW3c9S30VhRfV55nv+UYITgvT1b9bquYQ3tzuMNWc/3rrDHnvAgw/CWmtFF6RIJxQqwTGz7wOPu/sSMzsX2AGYWIhxOGb2JlDn7u+bWX+g3t03y7FvHUGJmIInOFIYlTyGRtrK9/yjLqoQpNa/WbdpIU8u3b1tcnPEEUE5hrXWIpEI1vZLJCILVSRs45PJzW7AcOAW4PoOXpOvdd39/eT9D9D087KWbQyNtC8xP8GkWZNIzE9kfRx3GmRcYNlKMdTVwY49XuH+pv3YgNaCmZx7Llx4IZip9IJUqmQdEvYDbnT3x8zs1/m+2MyeBtbL8tQ56Q+Sq6R3q7k6rrXwiqVc1kqJq0KV5ChnSnA60JnaUbmSlNovnuDv1d+nB8mCmT16BONvxo5tea1KL0iFWmhmNxCsdH6pmfWiEy3L7j4813Nm9qGZ9U/rovqoO4G6+43AjRB0UXXnveKuGN1JYawZEyeFKKpa7pTgtKOzrSpZk5R/BgUze6QXzLz/fhje9rys0gtSoQ4D9gF+4+6Lk4nILwv03tOAo4BLkrcPF+h9pQPFmpKda7qzBg+vKLPF69AtDmXWu7MqqgVMCU47Otuqkp6k9OrZzJjXfwV3Xtq6Q0bBzMzWodQqx337BregVhyJN3f/Gngg7fH7wPu5X9EplwBTzOxY4B2CZAozGwyc5O7HJR/PAjYHVjWzBcCx7v5EgWKoSJ3pTirEzCgNPl5RoYqqljMlOO3obKtKKkmZ9dRSxv71aNa+896W516yHWi65hEGbx0UzGyvdUhjcaRSmJkBRwIbufuFZjYQWM/du11R3N0XAcOybJ8DHJf2ePfuHkvayrc7qRDJSSUv4NeRSi+qGutZVN2dlZRKWCZOzD/RqN30E858cjhrz2hNbh5lP+rsLzz1Wms18GytQ+1tF4mp64Ba4AfJx0uAa6MLRwqldkAt43Yf1+4FtRAzo7SAn+QS2xacQs1KyregZSIBr9z3H46aMoKVFsxt2X5D9cn81K+kulePNi1AuVqHNBZHKswQd9/BzF4CcPfPzKwm6qAkHIWYGaXBx5JLbBOcMGclJRLwq7pnmNpwICuxCAA3Y8a+v8EO/BkTPrUVZmHlqizebsVxkfhZbmbVgAOYWT+gOdqQJCyFSk4qresFNLA6H7FNcMJsCfngqin8qWEMvQlqSjVU9+ao6ruY+sQh1MzM3XqUq3Uo31YjkRi4CngQWMfMLgJGAedGG5KEqRKTk+7SwOr8xHYMTlfGz3RasmDmwfcc3pLcfEQ/fndgPVObDtE4GpEOuPtdwJnAJOA9ggRHXVQi7dCqzvmJbQsOFLklpHHFgpmL1t6MBb+fzh7rb0TN4xpHI5KLmfUBTgHWJ1iv5jrgVOAR4BXgruiiEyltWtU5P7FOcIpmyRI47DB4PK2m1J570veBB+ibLJipcTQi7boD+AxIEEzZ/hVgwEHu/nKUgYl0RxhjYzSwOj8Vl+B0pvRCVgsXwn77wSuvtG478ki45Rbo1atlk8bRiLRrI3ffGsDMbiZY3G+guy+NNiyRrgtzbEwxxi7FbeByRSU43Zk6nkjA63e/wg/v2Y9eH7cWzFxw9LncsemF1L1oSmhE8rc8dcfdm8xsgZKb6MXtAhe2cl50MI4Dlysqwenq1PFEAi6t+xN3NBxGL74MNvbowdwzb2Sb3x2jVYdFOm9bM/sied+AlZKPjaD4d5/oQqtMcbzAha2cx8aUc3KWS0UlOO1NHW+v62rxZTdyX8PJ9CAomLm0Vx96P3o/U58frgrgIl3g7tVRxyBtxfECF7ZyHhtTzslZLhWV4ORaRC9n11VzM5xzDvs+dEnLe8y3ASy6eTrbDd+KulW06rCIxEMcL3BRKNd1fco5OculohIcyD74N73ratkymDABLhi3lJ1/fxRMmdKy3/vf2oEPr3+UwQf0b3kvzZYSkTiI4wVOOqdck7NcKi7BySbVdbVsWdBo89JTn+BPj4TmZ1p32n9/+t99N/1XXbXNazVbSkTiIm4XOKlssV3JuDNSLTHDh8N37D/83WupTU9uTjkFHnoIMpIbERGJn8T8BJNmTSIxPxF1KNINasFJqq2F3x7yd9Z7ciRrpxXMtMsvh9NOA7OIIxQRkWLTbLL4UAtOyr33stVpw1qSm8aeK/HAEfeTGHK6khsRkQqhOk/xoQTHHS69FEaPDgbhAA1rrsNeVfUcfs/BDBsWzLISEZH4S80mq7ZqzSYrc5XdRbV8eTC+5qabWrdtvjk37z+dZ363oda3ERGpMJpNFh+Vm+B88UVQMPOJJ1q31dXBAw+w/b/WpOZarW8jIlKJNJssHiozwVmwICiY+Y9/tG774Q/h5puhVy+tbyMiIlLmKi/BefnlILl5773WbePHwwUXtBlMrPVtREREyldlJTiPPw7f/z582Vowk5tugqOPjjQsERERKazKmUV1ww2w//6tyU2fPkHCo+RGRERiQAsUtlX0Fhwz2wy4N23TRsB57n5F2j4GXAmMAL4Gjnb3FwsSQHMzjBsHl13Wum3gQJg+HbbcstNvl6o63rcvLFqkMToiIhI9LVC4oqInOO7+JrAdgJlVAwuBBzN22xfYNPkzBPh98rZ7vvkGjjoKpk5t3bbjjvDII9C/f6ffLlV1PFWzqqoKevVKqz4uIpKnxPxEp6cid+U15SLOv1sYsi1QWOmfY9hjcIYBb7n7OxnbRwKT3d2B2Wa2hpn1d/f3u3ykTz6BkSPhmbSaUgccAHffDaus0tIS05kWmFTV8ebm4HFzs9bJEZHO68pf23H+Cz3Ov1tYUgsUpj5DLVAY/hic0cDdWbavD8xPe7wgua0NMzvBzOaY2ZyPP/64/SPddlvb5OYnP4EHH2xJboYNCyZPdWal4lTV8arkp1ZVpXVyRKTzulIOIM4lBOL8u4UltUDhxKETlSAmhdaCY2Y1wIHAuK6+h7vfCNwIMHjwYG935zPOgBdegHvvhd/9LiiYmZRqiensSsXp6+NoDI6IdFVX/tqO81/ocf7dwqQFCtsKs4tqX+BFd/8wy3MLgQFpjzdIbuu6qiq49VY4/njYa682T6VaYrqyUrHWxxEpD2a2FsEEh0HAPOAwd/8sY5/tCMb89QGagIvc/V6KrCvlAOJcQiDOv5tEx4JhLyEcyOwe4Al3vzXLc/sBpxLMohoCXOXuO7X3foMHD/Y5c+Z0OZ6ujMERkYCZveDug6OOoz1mdhnwqbtfYmZnA2u6+1kZ+3wHcHf/j5l9C3gB+B93X9zee3f3/CMiXZfv+SeUFhwzWwXYGzgxbdtJAO5+PTCdILmZSzBN/JhCHj9bMqOWGJHYGwnUJe/fDtQDbRIcd/932v33zOwjoB/QboIjIqUvlATH3b8C+mZsuz7tvgOnFOPYqQHFqe4oTekWqRjrps3E/ABYt72dzWwnoAZ4q9iBiUjxxb5UQ1cHFItI6TOzp4H1sjx1TvoDd3czy9kfb2b9gTuAo9y9Occ+JwAnAAwcOLDLMYtIOGKf4HRnQLGIlDZ3H57rOTP7MLWeVjKB+SjHfn2Ax4Bz3H12O8fKfxZnkWlRPJGOxT7BSZ/anW1AsQYbi8TWNOAo4JLk7cOZOySXr3iQYKHR+8INr2u0KF5lUBLbfbFPcCD3gGKNzxGJtUuAKWZ2LPAOcBiAmQ0GTnL345Lb9gD6mtnRydcd7e4vRxBvXrQkf/wpiS2MyqkmnkW28TkiEg/uvsjdh7n7pu4+3N0/TW6fk0xucPc73b2nu2+X9lOyyQ20LopXbdVaFC+mtLJzYVREC04uGp8jIuVGi+LFn1Z2LoyKTnA6Gp8jIlKKtCR/vCmJLYyKTnBAC/6JiEjpURLbfRU9BkdERETiSQmOiIiIxI4SHBEREYkdJTgiIiISO0pwREREJHYsKORdfszsY4LVSTuyNvBJkcPpLsVYGKUeY6nHB/nH+G1371fsYEqVzj+hK/UYSz0+iFeMeZ1/yjbByZeZzXH3wVHH0R7FWBilHmOpxwflEWM5KYfPUzF2X6nHB5UZo7qoREREJHaU4IiIiEjsVEKCc2PUAeRBMRZGqcdY6vFBecRYTsrh81SM3Vfq8UEFxhj7MTgiIiJSeSqhBUdEREQqTCwSHDMbYGYzzex1M3vNzE7Lso+Z2VVmNtfM/mFmO5RgjEcmY3vVzJ4xs21LLca0fb9rZo1mNqrU4jOzOjN7ObnPX8KKL98YzWx1M3vEzF5J7nNMyDH2NrPn0o5/QZZ9epnZvcl/L8+a2aAwYywnOv+EE1/avqGfezoTo84/HcYY3vnH3cv+B+gP7JC8vxrwb2CLjH1GAH8CDNgZeLYEY9wFWDN5f99SjDH5XDXwZ2A6MKqU4gPWAF4HBiYfr1NqnyHwK+DS5P1+wKdATYgxGrBq8n5P4Flg54x9TgauT94fDdwb5udYTj86/4QTX/K5SM49nfgMdf7pOMbQzj+xaMFx9/fd/cXk/SXAG8D6GbuNBCZ7YDawhpn1L6UY3f0Zd/8s+XA2sEFY8eUbY9JPgPuBj0IML9/4jgAecPd3k/uVYowOrGZmBqxKcIJpDDFGd/cvkw97Jn8yB+ONBG5P3r8PGJaMVzLo/BNOfEmRnHtA558Cxhja+ScWCU66ZFPW9gRZYbr1gflpjxeQ/R9Q0bUTY7pjCf7ii0SuGM1sfeBg4PfhR9UmjkFk/wy/A6xpZvVm9oKZjQk7tpR2YrwG+B/gPeBV4DR3bw45tmoze5ngQvGUu+f89+LujcDnQN8wYyxHOv90X6mfe5KxDELnn+7EFsr5p0d3Ay0lZrYqQXZ/urt/EXU82eQTo5kNJTjB7BZmbGnHby/GK4Cz3L05qj/oO4ivB7AjMAxYCUiY2Wx3/3cJxfg94GVgL2Bj4CkzmxXm/7Pu3gRsZ2ZrAA+a2Vbu/s+wjh9HOv90X6mfe0Dnn0II6/wTmxYcM+tJ8IXe5e4PZNllITAg7fEGyW2hySNGzGwb4GZgpLsvCjO+5PE7inEwcI+ZzQNGAdeZ2UElFN8C4Al3/8rdPwH+CoQ9WLujGI8haMZ2d58L/BfYPMwYU9x9MTAT2CfjqZZ/L2bWA1gdCP3/x3Kh808o8UV67gGdfwqt2OefWCQ4yb65W4A33P3yHLtNA8YEkxlsZ+Bzd3+/lGI0s4HAA8CPws74k8fvMEZ339DdB7n7IIK+0ZPd/aFSiQ94GNjNzHqY2crAEIJ+6FDkGeO7BH/hYWbrApsBb4cTIZhZv+RfTpjZSggk7qMAAAHySURBVMDewL8ydpsGHJW8Pwr4s7tr0awsdP4JJ74ozz35xojOPx0K8/wTly6qXYEfAa8m+/UgGCk+EMDdrycYdT8CmAt8TZDFllqM5xH0M16XbIJt9HCLo+UTY5Q6jM/d3zCzx4F/AM3AzSF3veTzGU4EbjOzVwlmFJyV/GsvLP2B282smuCPnCnu/qiZXQjMcfdpBCfJO8xsLsEgxNEhxldudP4JJ76o6fxTGKGdf7SSsYiIiMROLLqoRERERNIpwREREZHYUYIjIiIisaMER0RERGJHCY6IiIjEjhIcERERiR0lOCIiIhI7SnAkFGY208z2Tt7/tZldHXVMIlIZdP6pTHFZyVhK3/nAhWa2DkGF2wMjjkdEKofOPxVIKxlLaMzsL8CqQJ27LzGzjYBzgNXdfVS00YlInOn8U3nURSWhMLOtCWqQNLj7EgB3f9vdj402MhGJO51/KpMSHCk6M+sP3AWMBL40s30iDklEKoTOP5VLCY4UlZmtDDwA/Nzd3yCoZHt+tFGJSCXQ+aeyaQyORMbM+gIXAXsDN7v7pIhDEpEKofNP/CnBERERkdhRF5WIiIjEjhIcERERiR0lOCIiIhI7SnBEREQkdpTgiIiISOwowREREZHYUYIjIiIisaMER0RERGJHCY6IiIjEzv8D7ALi7T5pBuAAAAAASUVORK5CYII=\n", "text/plain": [ "