{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Illustrating calculation of effective sample size (ESS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create two samplers of an underlying $\\mathcal{N}(0,\\sigma)$ process: an independent sampler, using numpy's built in pseudorandom number generator, and an AR1 process.\n",
    "\n",
    "The marginal distribution (i.e. across all time points) for the AR1 process is $\\mathcal{N}(0,\\sqrt{\\sigma^2/(1-\\rho^2)})$, so we correct its standard deviation (through judicidious choice of $\\sigma$) to ensure it is the same as that of the independent sampler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "import numpy as np\n",
    "\n",
    "def independentSampler(sigma,T):\n",
    "    return np.random.normal(0,sigma,T)\n",
    "\n",
    "# to make comparable with independent, sigma = sigma1^2 / (1-rho^2)\n",
    "def AR1(rho, sigma, T):\n",
    "    sigma1 = np.sqrt(sigma * (1 - rho**2))\n",
    "    vX = np.zeros(T)\n",
    "    for t in range(1,T):\n",
    "        vX[t] = rho * vX[t-1] + np.random.normal(0,sigma1)\n",
    "    return vX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the time course of each of these samplers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApMAAAGrCAYAAACc1xygAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAIABJREFUeJzsvXt4XGW5/n+vWXNMJskkTdqmObS0pVAQKIWiloKgIOhGUPCHCgVROchGYVu1/e6N/ABFZeu2eysKVlCUL3JQQEQ5uFFOLRRokXIolDalbZI2bdokk2Qyh7XWrPf7x5p3zZqZNTNrJjOZNZPnc129lMxk8maSrLnnfp7nfgTGGANBEARBEARBFIGj0gcgCIIgCIIgqhcSkwRBEARBEETRkJgkCIIgCIIgiobEJEEQBEEQBFE0JCYJgiAIgiCIoiExSRAEQRAEQRQNiUmCIAiCIAiiaEhMEkXxiU98At/73veK/vwVK1bgpptuKt2BykS1nJMgiPKjKAoEQcBzzz1X6aPkpNLn/P3vf49zzjlnyr+uqqpYsmSJ7X8+tQiJSaIonnzySdxwww2VPkbV8dxzz0EQBCiKUumjEARBWOamm27CihUr8t4vGo3i29/+Nm655ZaM25599lkIgoDPf/7zpo8viiL8fj8aGhowd+5cfOMb30AsFtPvs3fvXpx33nmYO3cuBEHAXXfdlfIYDocDN998M6699toivkNiMpCYJAiCIAiiJNx///3o7OzEkiVLMm67/fbbMWPGDDzyyCM4cOBAxu0f/vCHEQqFMD4+jr/85S+47777cOutt+q3OxwOfPzjH8d9992Hzs5O069/zjnn4NChQ3j66adL900ReSExSRTFaaedhu985zsAAEEQcNttt+Hkk0+G3+/HMcccgw0bNuj3VRQFq1evxuzZs9HW1oZ///d/z3i8vXv34qKLLkJHRwdmzpyJL3zhCzh48GDK1/va176G888/Hw0NDTj88MNxzz33pDzGK6+8gtNOOw0zZszA3LlzccMNN6Q4gFN1zuuuuw4XXXQRmpqa0NXVhTvuuAMA0Nvbi0984hMAgEAgAL/fjx/84AcFPe8EQUwdg4ODOP/88xEIBDB//nw8+OCDGfexct1Zu3Ytli1bBr/fj5NOOgmbN29OeYx77rkHxx13HJqamnD00UfjgQce0G/j1YyHHnoIixYtQkNDA84880zs3bu35OfMdn38/e9/jx/84AfYuHEj/H4//H4/1q9fb/qcPfzwwzjrrLMyPj4wMIBHH30Ut912G1paWnDnnXdme9oBAMceeyxOOeUUbNq0Sf9Ye3s7rrnmGpx88skQRdH080RRxBlnnIFHHnkk5+MTJYYRRBF85CMfYddffz1jjDEA7Nhjj2U7duxgsiyz6667jnV3d+v3veWWW9jcuXPZ1q1bWTQaZd/5zneY0+lkN954I2OMsWg0yo444gj2zW9+k4VCITY+Ps5WrlzJzjjjjJSv5/V62WOPPcZkWWaPP/44c7lcbMOGDYwxxrZt28b8fj+7//77mSzLbPfu3ezYY49lt9xyi/4YU3XOxsZG9o9//IPF43H20EMPMYfDwXbs2MEYY+zZZ59lAJgsy6X9gRAEUXLOPPNMdtZZZ7GhoSE2NDTE/uVf/oUBYM8++yxjzPp1Z8GCBfp15cYbb2Stra0sGAwyxhi7++67WVdXF9u0aROLx+Ns/fr1rKGhga1fv54xlrxmXHTRRSwYDLJgMMiWL1/OLr300pKfM9f18cYbb2Qnn3xy3uds1qxZ7IEHHsj4+E033cRmzpzJJElia9asYV1dXUxRFNPHV1WVvfbaa6y1tZV961vfMv06c+fOZXfeeafpbT/+8Y/ZCSeckPesROkgMUkURbqY/N3vfqff9vbbbzMAbP/+/YwxxhYuXMh+9rOf6bcrisLa2tp0kfbwww+zOXPmMFVV9fv09/czAKyvr0//eueff37KGS688EL25S9/mTHG2Ne//nX2hS98IeX2e++9ly1YsED/76k655e+9KWUc7S2tuoXVxKTBFEd8L/tN998U//Ym2++mSLSrF53jNeVeDzO2tvb2T333MMYY+yYY45h69atS3mMyy+/nH3lK19hjCWvGXv27NFv//nPf86OPPLIkp8z1/XRqph0uVzsySefTPmYoiiso6ODrV69mjHGWE9PDxMEgT366KP6fW688UYmiiJrampiHo+HAWCf//znWSgUMv06ucTkr371K3bYYYflPStROpxT54EStcycOXP0/19fXw8AGB8fx6xZs9Df34/58+frt4uiiO7ubv2/d+zYgQMHDqC5uTnlMT0eD3p7e/XeGONjAMBhhx2Gf/7zn/pjPPvsswgEAvrtqqpCVdUpP6fxa/CvMz4+DoIgqof+/n4Aqdedww47LOU+Vq87xsdwOBzo7u5GX1+f/hirVq3C6tWr9fsoioJTTz015THSr138mlLKc+a6PlqlpaUFo6OjKR/785//jH379uHyyy8HACxYsACnn346br/9dpx33nn6/T70oQ9hw4YNUBQFv/nNb/Dd734XIyMj+lmsMjo6ipaWloI+h5gcJCaJstPZ2Yldu3bp/x2Px/ULKQDMnj0bc+fOxc6dO3M+jvEx+H9zATd79mxcdNFF+M1vflPxc+bC4aA2ZYKoBvi1ZdeuXfjABz6g/38jVq87xs9TVTXlzefs2bNx880349JLL634OXNh9dp1wgknYOvWrSkfu/322wEAp5xyiv6xUCiEcDiMnp4eLFy4MOX+TqcTV155JZ5//nl8/etfx5/+9KeCzvrWW2/hxBNPLOhziMlBr2xE2fniF7+In/zkJ9i2bRtisRi++93vYnh4WL/9/PPPhyzLuOGGG/R3tIODgxlN5E888QQef/xxxONxPPXUU3j00UfxpS99CQDwr//6r3jooYfwxz/+EZIkIR6Po6enB0899dSUnzMXs2fPBgC89957lj+HIIipp6OjAx/72MewevVqjIyMYGRkBNdff33Kfaxed37605/i3XffhSRJ+P73vw9JknDuuecCAP7t3/4N3/ve97Bp0yaoqopYLIZNmzbhtddem/Jz5mL27Nno7e1FNBrNeb/zzz8ff/vb3/T/3r59O5555hn89re/xZYtW/R/27dvR2dnpz6gaMbNN9+Mv/71rymDktFoFNFoFIwxKIqCaDQKWZb12+PxOP7+97/jM5/5jOXvjZg8JCaJsrNmzRqcf/75+MhHPoLOzk5IkoQPfvCD+u0NDQ3YuHEjent7ccwxx6CxsRHLly/HCy+8kPI4X/7yl/HrX/8agUAA11xzDX75y1/q73SXLVuGp59+GnfeeSc6OjowY8YMfPazn8WePXum/Jy5WLRoEa699lqcfvrpCAQCKbEXBEHYi3vvvRdutxvz5s3D0qVL8bnPfS7ldqvXnauvvhqXXHIJWlpa8Nhjj+GJJ57QS87XXXcdbrrpJnz1q19FS0sLOjo68O1vfxsTExNTfs5cfO5zn8MRRxyBOXPmIBAIpAg8IxdffDF6e3uxZcsWAMAdd9yBww8/HCtXrsTs2bP1f3PmzMGqVavw29/+FpFIxPSxFi5ciEsvvRRr1qzRP+bz+eDz+dDb24urr74aPp8PV1xxhX77448/jpaWFnz84x+3/L0Rk0dgjLFKH4Ig8nHaaadhxYoVpkG4BEEQdkUQBDz99NM444wzKn2UKeP3v/897rvvPjz++ONT+nVVVcUJJ5yAn/zkJ/joRz86pV97ukM9kwRBEARBlIyLL74YF1988ZR/XYfDgddff33Kvy5BZW6CIAiCIAhiElCZmyAIgiAIgigaciYJgiAIgiCIoiExSRAEQRAEQRRNVQ/geDwetLW1VfoYBEFUKQcPHkQsFqv0McoCXR8JgpgMhVwfq1pMtrW16aukCIIgCoVvD6lF6PpIEMRkKOT6SGVugiAIgiAIomhITBIEQRAEQRBFQ2KSIAiCIAiCKBoSkwRBEARBEETRkJgkCIIgCIIgiobEJEEQBEEQBFE0JCYJgiAIgiCIoiExSRAEQRAEQRQNiUmCIAiCIAiiaEhMEgRBEARBEEVDYpIgCIIgCIIoGhKTBEEQBEEQRNGQmCQIgiAIgiCKhsQkQRAEQRAEUTQkJm2CHFfxmdtfxMOv9Vf6KARBEESB7B+N4pQfPYONO4cqfRSCmHJITNqE4QkJr/cGsWn3cKWPQhAEQRTI1n2j6BuO4OF/kiFATD9ITNoESVEBALHE/xIEQRDVw2hEBgC81HMIjLEKn4YgphYSkzZBinMxGa/wSQiCIIhCCYY1MblvNIrdQ+EKn4YgphYSkzZBdyZlciYJgiCqDe5MAsCLPYcqeBKCmHpITNoEKnMTxPRix44dWL58ORYtWoRly5Zh69atpvd76623cNppp2Hx4sVYvHgxHnnkkSk+KWEFo5h8aSeJSWJ64az0AQgNKnMTxPTiqquuwpVXXonLLrsMDz30EC677DJs2rQp5T7hcBjnnXce7rnnHqxYsQLxeBwjIyMVOjGRi7GEmFw404+NO4egqgwOh1DhUxHE1EDOpE0gZ5Igpg+Dg4PYvHkzVq5cCQC44IIL0NfXh56enpT73XffffjQhz6EFStWAABEUURra+uUn5fIz2hEhlt04GOLZ2IkLOOdgbFKH4kgpgwSkzaBeiYJYvrQ19eH9vZ2OJ1acUgQBHR3d6O3tzflfu+88w48Hg/OOeccLFmyBJdeeikOHjxYiSMTeQhGZDT6XDh5gSb2qdRNTCdITNqEmEJlboKYzpjFyciyjL///e9Yt24dXn/9dXR0dODqq682/fy1a9eis7NT/xcKhcp9ZMLAaERGk8+JZfNa4BYdeLGHwsuJ6QOJSZuQ7JkkZ5Igap2uri4MDAxAURQAmpDs6+tDd3d3yv3mzp2L008/HR0dHRAEAStXrsTLL79s+pirVq1Cf3+//s/v95f9+yCSjEZkBOrc8LlFHN8dwKu7hvWKE0HUOiQmbQL1TBLE9GHmzJlYunQp7r33XgDAww8/jM7OTixcuDDlfhdeeCE2bdqEsTGt/+6JJ57AcccdN+XnJXLDGEs4ky4AwMkLWxGR49jSF6zwyQhiaiAxaROSPZNU5iaI6cC6deuwbt06LFq0CLfeeivuvvtuAMDll1+Oxx57DADQ3d2N//iP/8Dy5ctx7LHH4plnnsEvf/nLSh6bMCEqq5AU1SAmZwDInje5NxihliaipqBoIJsgJS4s5EwSxPTgiCOOwMaNGzM+ftddd6X89yWXXIJLLrlkqo5FFAHPmORi8tjOAOrdIl7sOYRvnLko5b4v9hzCpb95FdectgCrPn7ElJ+VIMoBOZM2gfdMKiqDEidBSRAEUS2ki0mX6MCH5s/Alr4gJmKKfr/9o1Fc98DriKsM2/aPV+SsBFEOSEzaBGOjtkRikiAIompIF5MAsHxhKxSV4dVdwwAAOa7ia/f9E4dCEtxOB/YGIxU5K0GUAxKTNsEoJilrkiAIonoIhiUAqWIyvW/yR09tw+Y9I/jKisOwpDNAYpKoKUhM2oSYwY2kvsniWf3QG/jRU9sqfQyCIKYRZs7kEbMa0Op348WdQ3jq7QHcuX4XTpjbjP/ziSPR0exDMCynlMAJopohMWkTZCUZWExTfsXzzLaDeGEHbQghCGLq0MVkXVJMCoKADy9oxbsDY/jWH99ES70bP7/oeLhEBzoCPgAgd5KoGUhM2gQpnhSQ5EwWT0yJU1AwQRAlQ1JU/OPdA6YbijhjCTEZMDiTAHDyAq3UPSEp+Onnl6C9SRORHc0JMTlCYpKoDUhM2gTqmSwNMUUlMUkQRMl45J/9+MrvNuPl94ez3seszA0Apx85E00+F9acfSROObxN/zh3JvvJmSRqBMqZtAkpYpLK3EXBGIOkqOTsEgRRMnoGtR3nuUrSwYSYbEwTk7Mavdjy/58JQRBSPk7OJFFrkDNpEyQawJk0/HkjZ5IgiFLROxwGABwKxbLeZzQiw+N0wOsSM25LF5IAqGeSqDlITNoEciYnD4lJgiBKjS4mx3OLyUCdK+vt6XhdIlr9buwjMUnUCCQmbUKMeiYnTYxWUhIEUUIYY7qYPJjHmUzvl8zHnICPytxEzUBi0iakOpMkhoqBi3ApruacvCQIgrDC0ISEsKS9Sc1Z5g4XLiY7Aj4cGI9SJYWoCUhM2oTUnkkqcxeD8XmjlZQEQUwW7koCwKFxyfQ+jLGinMmOgA+Mafu6CaLaITFpE8iZnDxRQ3sAvdsnCGKy9BnFZBZnMizFoagMTT53QY/NJ7r7g+E89yQI+2MLMRmNRvHpT38aixYtwnHHHYczzzwTPT09lT7WlEI5k5MnRoKcIIgSsmdIE3pdLT4MhyUoJhWPbBmT+dAnuqlvkqgBbCEmAeDKK6/Ee++9hzfeeAPnnXceLr/88kofaUqR4ip4ggSVuYsjpcxNYpIgiEnCy9wndDeDMWB4IrPUHQwXKSabKR6IqB1sISa9Xi8++clP6nlcH/rQh7B79+7KHmqKkRQVfreWIU+uWnEYnzcSkwRBTJbe4TCafC4saPMDMJ/oTjqThe0A6QzUASBnkqgNbCEm0/nZz36G8847r9LHmFIkRUWDl8TkZDC2B9BzSBDEZOkbDqO7pQ6tDR4AwEGTrEldTBaQMwkAjT4n/B4nOZNETWA7MfmDH/wAO3bswA9/+MOM29auXYvOzk79XygUqsAJy4MmJrWLUUymMncxUJmbIIhSEZXj2D8W1cSkXxOTh0KZZe6xhJgMFDiAIwgCOgI+EpNETWArMflf//VfeOSRR/Dkk0+irq4u4/ZVq1ahv79f/+f3+ytwyvIQi5MzOVlSytxxEuQEQRTP3mAEjAFdLXVoa+BiMrszmb6X2wodzT4MBKNQVcrFJaob24jJtWvX4v7778fTTz+NQCBQ6eNMKYwxyCQmJw1tESIIolT0Jia5586oQ6tfcx3NVioGI5pbWegADqBNdEtxNed2HYKoBgrrGC4T/f39+OY3v4n58+fj9NNPBwB4PB688sorFT7Z1KCoDIwBfl7mpmnuojC2B8QotJwgiEnAJ7mNZe7cAzjFOZMA0D8SwaxGb7FHJYiKYwsx2dnZOa3X3/H+Po/TAbfTQa5akdA0N0EQpcIoJr0uEQ1eZ5YytwKgeGcS0ErqJ8xtnsRpCaKy2KbMPZ3hwsftdMDjdFCZu0hSnEl6DgmCmAS9w2GIDgHtTZpj2Ob3mK5UHI3IqHOLcDsLfznVsyYpHoiockhM2gC+R9otOuBxilTmLhJyJgmCKBV9w2F0BHxwitrLZKvfY+5MhqWiXEkg6Uzuo4luosohMWkDjGVuciaLh8QkQRClgDGG3uEw5s5Ipoq0NrhNVyqORuSixWSb3wO36KB4IKLqITFpA2LGMreLeiaLxejokrtLEESxHApJCEtxdLUkxWSb32O6UnE0IhcVCwQADoeA9oCXytxE1UNi0gboPZNU5p4URhFOziRBEMViHL7hmE10qyrDaERGoEgxCUAPLp/OQ6hE9UNi0gboPZNU5p4UVOYmCKIU9JmJSZOViiFJgcqKm+TmdAR8CMUUjCWmwgmiGiExaQO48HGJJCYnQ8o6RcqZJAiiSHI5k8aViqPh4jMmOXrWZDBc9GMQRKUhMWkDUqKBXCLt5i6SlA04JMgJgigSXUwaBnDMVipOJrCco2dNUt8kUcWQmLQBfI80lbknR9QgwqnMTRBEsfQOhxGoc6HRmxSJZisVxxJiMlA3eWeSJrqJaobEpA1IjwZSVJYRP0HkJ6ao8HtovzlBEJOjdyicUuIGzAdwggkxWew0NwB0BrSvQ84kUc2QmLQBsbRpboB6/oohJqto8HIxSa0CBEEUTlSOY/9YNCUWCIDpSsVSlLlnN3khCORMEtUNiUkbIKXlTAKgrMkiiClxXUxSmZsgiGLoTziEc9PEJJC5UrEUYtLtdGBWg5fEJFHVkJi0AXJcyxfjPZMAlWmLIaaoqHM7IToEEpMEQRSFWSwQJ32l4qjeM+me1NfsaPZRmZuoakhM2gApUZI1lrmpTFs4MUWFx+mAW6QhJoIgisMsFoiTvlKxFM4koE10D01IiEh03SeqExKTNiA9tBwgZ7IYYnIcHpcIt9NBziRBEEWxZ0gTk+k9k0DmSkWeM9mYaK8pFproJqodEpM2gHomSwN3Jj1OBw0wEQRRFL3DYTgdAtqbvBm3pU90j0Zk+D1OOMXJvZTOCZCYJKobEpM2IDUaiMrcxaDEVSgqgzfhTNLzRxBEMfQNh9HZ7DMViOkrFUcj8qRL3ADQmRCT+0hMElUKiUkbEONlblGkMneRcCfS43RQmZsgiKJgjKF3OGxa4gYyVyqWSkzqZW4awiGqFBKTNiClzK2LSXLWCoG3BXB3l8QkQRCFcigkISLHTYdvgMyVisGwVBoxSc4kUeWQmLQB6bu5AeqZLJRoQnx7nDSAQ1QHO3bswPLly7Fo0SIsW7YMW7duzXpfxhg++tGPIhAITOEJpx+5JrmB1JWKqsowHlNKIibrPU74XCIOTUj570wQNoTEpA0wdyZJDBWC7ky6HPBQNBBRBVx11VW48sorsX37dqxZswaXXXZZ1vv+93//NxYsWDB1h5um9A5PAMglJpMDOONRBYxNPhaI01LvxvBELP8dCcKGkJi0AbzfzyUKVOYukphxiMlFziRhbwYHB7F582asXLkSAHDBBRegr68PPT09GffdunUrHn30UaxZs2aqjznt2Dmoicm5M+pNbzeuVEwGlpdGTM7wuzEcImeSqE5ITNoAyWQ3NzlrhREzlrlFhz7URBB2pK+vD+3t7XA6tXxCQRDQ3d2N3t7elPvJsowrrrgC69at0+9LlI8tfUF4XQ4smuXPeh++UjEY0YRfY4mcyeY6N4YmJDDGSvJ4BDGVkJi0AZKiwi06IAgC5UwWidGZ5D2TdFEmqgmz39ebb74Z559/PhYvXpz389euXYvOzk79XygUKscxa5a4yrClL4hjOwM5cyP5SsVSbb/hzKh3I6aoiMhUlSKqDxKTNkCKq3AnyttU5i6OlJ7JxHNIweWEXenq6sLAwAAURQGgCcm+vj50d3en3O/555/Hbbfdhnnz5mHFihUYGxvDvHnzcPDgwYzHXLVqFfr7+/V/fn92d43IpGcwhFBMwfHduYec+EpFvgWnlD2TADBEpW6iCiExaQNiilFMUpm7GGJp09wAqG+SsC0zZ87E0qVLce+99wIAHn74YXR2dmLhwoUp91u/fj327NmD3bt3Y8OGDWhsbMTu3bvR1tZWiWPXNK/3jgAAju9qznk/vlJx1yGtv7JkYjIxKT5ME91EFUJi0gbwMjcAmuYuEv58eV0OXUzSc0jYmXXr1mHdunVYtGgRbr31Vtx9990AgMsvvxyPPfZYhU83/Xi9NwgAWJrPmUxMdPcMam0EJRvAqScxSVQv1NFtA2RjmVvvmZx8mfuZbQcwHlVw3pKOST+W3TE6k9zdLdSZ7BkMYfehCZxx1KySn48g0jniiCOwcePGjI/fddddpvefN28egsFguY81bXm9bwQdAR9mNmbu5DbCVypyMVkqZ7K5LlHmJjFJVCHkTNoAqUxl7v/5+w7c8vi7k36casC4AafYMvfP/rEDX733NcRVGtwhiOnEWFTGjsFQ3n5JIOlMlrrMPSNR5h4hMUlUISQmbYAUL0+ZOyzFMZaYOKx1ojJ3Jh36c1nocxiKKVBUBpkGdwhiWvFGXxCMAcd35+6XBJIrFfn1pcFbqgEc7XHJmSSqERKTNiDVmSzdNHdEiiOmqLrQqmX0aCBX8QM4/DlXyJkkiGkF75e05ky69f/f6HVCdAglOUOL3jNJW3CI6oPEpA0wiklBEOB2OkqSM8nzysajyqQfy+6kbMDRo4EKE9H8OVfImSSIacXrvSNwiw4cPacx7315mRsAmko0fANowtTpEGgAh6hKSEzaAElRdQEEaIKoFGXuiKSJqbFo7Ze6kwM4hv3mBQryaOIx5Dg5kwQxXWCM4fW+II7uaNR71nPBVyoCpeuXBDQjobneTWVuoiohMWkDYoaeSUAbwplsmZsxpjuT06FvMhlanixzF7pSkT8GDeAQxPRh16EJBMNy3nxJI20Jd7KUYhLQ4oFoAIeoRkhMVhjGWEqZGyiNM2n8fLuUufuGw2VbcZha5i4uGog/Bg3gEMT0oZB+SQ4vdQd87jz3LIwWciaJKoXEZIXhJVWX0Zl0Tb5nMiwlnU07lLnfHRjDKT96Fn/bur8sj28scxc7gMMHlWgAhyCmD6/3JTbfFCImGzQR2VhiZ7Kl3o3xqELbu4iqg8RkheH7o1OdycmXuSOGCe6xSOWdyT1DWibb7qFwWR4/6UyKRUcD8fvTAA5BTB9e7w2ircGDjoDP8ueUs8wNACNhcieJ6oLEZIXh70BLXeaO2MyZ5IJ2vExnickqHALgEgV9i1Cx0UA0gEMQ04OwpGDb/nEs7Q5AEKxH/LSWSUw2J8TkUIjEJFFdkJisMLqYFEsrJqMpzmTlxeRo4gzl6t+MKXF4nKIWrSRyMWnd3WWMIcqjgVRyJgliOvBm/yjiKrMUVm6Er1QkZ5IgNEhMVhjJMDjC8bjESe/mTilz28GZTJwhVCYxGZVV3ZF0F7FFSDKUtsmZJIjpgT5802W9XxIAlnQF4HOJOKajqaTnoS04RLXirPQBpjs8WLvUZW7jAI4dprm5OzpWVmcydb95IWVu4/NNPZMEMT14vXcEokPAMZ2FicLF7Y1493tnl/w8+hacEG3BIaoLciYrTCxLmVtR2aRETUrPpA3K3FxElq1nUlF1EalPcxfw/BnbAmiamyBqHx5WfuTsBtS57eGrJFcqkjNJVBckJisML6mmT3MDhYmhdFJ6Jm3kTJavZ1I1OJOFl7mNUUyUM0kQtU//SAQHx2NYWmC/ZDnhYpLK3ES1QWKywphOc7uKWwdoJGLXAZxYuZzJeEbPZPFlbnImCaLWeb2v8LDyctOc2PVNAzhEtUFissJkiwYCCs9JNMLL3IIwPQZwYnKyzF3M80dlboKYXmzdOwoAOK7A4Zty4hQdCNS5KBqIqDpITFYYfQBHzCxzTya4nDv4cwvuAAAgAElEQVSTM+o9NhnA4T2TSllWKhrL3MlpbuvPX4ozSdFABFHz9AcjAICu5roKnySVlno39UwSVQeJyQpTbmdyZoMHYSle8T5A7kwqajLPsZQYp7mTOZOFlLkNziSVuQmi5hkIRtDq96Rce+1ASx2JSaL6sNdf0TQkZpozWbqeyVmNWm5ZJd1JOa6mRRWVtuzOGENMUeF1aY6uU3RAdAhF90xWWngTBFF+9o9GMSfgrfQxMmipd2MkLEGldhuiiiAxWWHMncnSlblnNWoXy0oO4aR/7fFYaYWtHGdgLFWQu8XCsjpj1DNJENOGuMpwYDyG9ib7ickZfjdUlhxaJIhqgMRkheHxPy6xtGXuKC9zczFZwSEcHk3kdGi7b0vtkkYTopuLcEAT5xRaThCEGQfHY4irDO1NvkofJQOKByKqERKTFSbbbm5g8s6kIABtiR2yfACmEnBnsj1RUip1mZu3A/D2AEB7DgvJ6UzNmSRnkiBqmX2j2vCNHZ1JvlKx1vsmVZXhD5v7sC8xCEVUNyQmK4x5zmSizD2JnsmwFIfPJaLJp+WWlWvzjBW4K9oR8CXOUlphG9OdSUOZu0BnMmocwKFpboKoafaPRgEAs20pJrVr9vBEba9U/Mub+7D6oTex7vmdlT4KUQJITFaYsk1zy5qYbPBqa8IqWuZOuKIdAS2Co+TOpD7ElFrmLigaiJxJgpg2cDdsTsCOZW7NmazlMndUjuNHT70HANgxGKrwaYhSQGKywvBSrMdUTBZf5o7KcXhdIhq92rvcipa5uTPZXCZnUjZ7DkWKBiIIwhTdmWy0nzM5I9EzOVLDYvJ3L+3G3mAEggDsPEhishYgMVlhkj2TSVctOc09uZxJn1tEk6/yziSfSuwsd5nbNYkyt8GZjFOZmyBqmoHRKAQhmXZhJ2p9AGdkQsLPn+1BZ7MPn/xAOw6MxSrahkWUBhKTFSZWxt3cdW6jM1n5aKCyOZMmZW5PodFABmdSpmgggqhpBkbtGVgOJMVkrQ7g/OyZHRiPKlhz9pFY3N4AAHj/4ESFT0VMFvv9JU0zeJnbvGdyEtPcUqLMnRjAGatgaHnmAE65eiZTBTlFAxEEYcbAaBRzbDh8AwBel4g6t1iTYnLXoQn83417cFxXAOcc244FbX4AQA/1TVY9JCYrjJwztHzyAzgepwNu0VHZae5Ev+bsJi8EAQiVOLScB47zDThAIrS8AFEYNYSW0wAOQdQuSlzFgbGoLSe5OS31bgyFak9M/uipbVBUhus/uRiCIGDhTE1MUt9k9WMLMXnttddi3rx5EAQBW7ZsqfRxphTdmSxhaDljTBeTgiCgweus+ACO1+WA1yXC73GWscyd2TPJmDVhmOJMUs8kQdQsB0MxqAy2DCznzEisVKwlNu8expNv78dZR8/CSYe1AAC6Z9RBdAgkJmsAW4jJz372s9iwYQPmzp1b6aNMObwU6xIF/WPJnsniytwxRQVjgM+tOXWNPlfFB3B472aDx1lyl5S7iukDOAAsB5cb+1Npmrv8MMbw0s5DtH+YmHL2BbVJbjsGlnNa6t0YmpAsvxm2O4wx3PL4u3A6BKw5+0j94x6niO6WOuyknsmqxxZi8tRTT0VnZ2elj1ERJEWF2+mAIBjE5CTL3Fxc6WLS66z4AA4PT2/wuqZmAIeLSYvPYcoADonJsvPSziFcdOcr+N93DlT6KMQ0g8cCtdswY5LTUu+BpKiYkIrvm7cT/+wdwZa+IL5wUjfmJ/okOQva6rFnaAIy9apXNbYQk9MZKa7CI6b+GCZb5g4nLkA+l9GZrGSZW9EHgRq8zpKfJVuZ23hbPqKyqj9fVOYuP3y4YM8QORLE1DJg41WKHH0LTo30Tb78/jAA4Jxj2zNuWzDTDznO0DscnupjESWkqsTk2rVr0dnZqf8Lhaq/zyKWcCaNTHaaOyKniUmvC6GYgniFSopjERmNiU08DV4nQrHCXdKXdh7KusM1uU7ROICj/f9CnMl6j3ZGKnOXH+5CHBir7ZVxhP0YGK2GMjffglMbfx+v7hqGW3TguK5Axm18onsnTXRXNVUlJletWoX+/n79n9/vz/9JNkcyEZOCIGjrAIvMmYxIaWXuRHB5qALuZFSOI6aoBmfShaisFlTSCEsKLv31q/jvp7eb3q5vwHFlZnVaF5Mq6twiBAFUbpkC+M9lcDxa4ZMQ042B0YhtA8s5+hacGhjCUeIqXtszgiVdgZTEDY4uJqlvsqqpKjFZi0iKCpeY+WPwOAsL3TYSTYvKaeDB5RUYwuFfkw/g+BMOZSF9k0MhCYrKcChk/i7dtMwtFlrmjsPrcsDlcFTMwZ1OcME+OF4bzgtRPQyMRtHm95hed+2CvgWnBsrc7w6MIxRT9AnudBa01QOgeKBqxxZ/TVdddRU6OzvR39+Ps846CwsXLqz0kaYMKZ7pTAJaybZ0ZW5NwI1WYAiHRxI1GXomgcKCy/m782znNy1zFzyAo8LjFOEUBdqAMwVIiVaCwTFyJompZSAYtXWJGwBa/LWzBeeVXUMAkFVMBurcaPW7Kbi8ynFW+gAAsG7dukofoWJIigq/J/PHMBlnkg/g1BmigYAKO5OJUjt3KAtxJvkFNbuYNNmAo0cDWRPkmph0wOkQaAPOFEDOJFEJlLiKwfEolpj07tmJlrraEZOv7hqG6BCwdG5z1vssaPPjnYExMMZSkk2I6sEWzuR0xqxnEtB6/ortmUwvcyf3c099zySPJNJzJosoc3NnMtsUuGnPZMHT3Nr6SZfooAGcKYA7xmEpXvKNSASRjcFxLbDczttvgKQzOVTlYlJVGTbtHsYH5jSamiacBTP9GI8qOJillYmwPyQmK0xZytxZBnAqsVKRC0DujvILSiFnGZ7Q7ltMmduqmNSdSVGATNFAZcc45HSASt3EFMFjgeYE7C0mGzxOuEQBI1UuJnsOhjASlrOWuDnJiW4awqlWSExWGCkhYtKZTJnbLBoIyO7slRMuAI2h5UCBzmTigiopasoObU5UVuESBYiOzOB3Kz2TjDHt5+BywOkgZ3IyjEVl/On1/rybO4ybiQYpHoiYIngs0Gwbr1IEtEQPvgWnmnlll5YvedJhM3Lejw/h9NAQTtVCYrLCSIqaspebU0oxqU9zV2QAJ1uZuwBn0hCPYfY9xJR4iisJFDaAw59nr1OESxQoGmgS3P7sTnzjwTewbf94zvvJSlJsTtd4oB07dmD58uVYtGgRli1bhq1bt2bc55lnnsFJJ52Eo446CkcffTRWr14NlZzzohlIrFKcY/MyN6BlTVZ7z+Sru4YhCMCyedn7JQFg4UzKmqx2SExagDFm6oilcygUw+W/24S9WcK1zR43a5nbJRa9mztbmdsOAzhcTBbSJ2fcAmFW6o6ZuLuFRAMZey5FhwDFxtPcUTlu6+iiF7YfBABM5Pn5GgejpqszedVVV+HKK6/E9u3bsWbNGlx22WUZ92lubsYDDzyAd955B6+99hpeeukl3HPPPVN/2Boh6UxWg5h0VbWYZIzh1V1DOGJWAwKJgaJszGnywetyUDxQFUNi0gJPvr0fx970v+gfyb3u6aWdQ/j7u4NYn3hBzQffAW3eMzkJZzJdTJZxAEdSVDy//WDWsib/mroz6SlimjtcuJgsJLTc2HOpDeDY0/mJynGs+M9ncMdzPZU+iikHx2N4Z2AMQP7nfbo7k4ODg9i8eTNWrlwJALjgggvQ19eHnp7Un+3xxx+P+fPnAwC8Xi+WLFmC3bt3T/Vxa4ZqCCzntNR7EIopRffOV5re4TAOjMXwwTz9kgDgcAiY3+rH+xRcXrWQmLTAe/vHIcXVvAn9vLfPquvG+8aylbkVlRUlbNLL3HVuEaJDKIsz+eiWvfjib17Fq4nemHT412zwpjqThfRvGpvQzb6HmByHJ22zAn9OJQsX4qicjBZyioIu8u3GwfEYDoUk2+axvbTzkP7/Y3l+b1MHcKafM9nX14f29nY4ndrfgyAI6O7uRm9vb9bP2b9/Px566CGcc845U3XMmmNgNIqZDfYOLOfoW3Ampr6iVAqs9ktyFsz0Y28wgrBE6Q7ViP3/omwA/+XON1nHSxI85zEf3L3JNs0NpA4qWCVdTAqCgEavsyzT3PsTZaP3D5kL7bGIDL/HCWfi4u2fRGg5YO5Mmg0xFTLNrTuTLhFOG2/A4d97pMj2h3LzwvakmJTzPO/897rB65yWzqQZuYaWxsbG8KlPfQqrV6/GiSeeaHqftWvXorOzU/8XCtnzTUclGRiN2H74hqNvwanS/dzcYFh2WO5+SQ4fwiF3sjohMWmBUEx78c43WcdFT75+MU5OMZko0xaTNcn7O40Cq9HnKkuZm3/P+7L0iY5FZH0DDwC4RAd8LtGye6uqDCPh5GOMhi32TBYxgONxOuASBSg2HXDgz7XVNytTCWMM63ck2zvyvQmSFBWiQ8DsRu+0DC7v6urCwMAAFEX7O2CMoa+vD93d3Rn3HR8fx9lnn41zzz0Xq1atyvqYq1atQn9/v/7P7/eX7fzViBxXMTgeq4rhGwBorq/u4PJXdw1jfms9ZjZYe771IRzqm6xKSExagItDq86k5TJ3TmeysJxEIxEpDp9LhMMQldPgdZalzM3F3d6RLGIyqugZk8azWO2ZHI8qiKsMh7XW64+Xjtk0dyHObrJnUosGsmuZO5h4rq0Mg001OwZDGByPoatFc33yTcTLcS3OaVajd1oO4MycORNLly7FvffeCwB4+OGH0dnZmbFKNhQK4eyzz8ZZZ52FG264oRJHrRkGx2NgVRBYzplRxWJyYDSC3uFw3nxJI8msSRKT1QiJSQvwMrdxEMSMgp3JxESrx7RnUhNDxTRfh6W4PnzDafS6yhINFEw8Zn9OZzJVTPoLKLnz53zuDE1Mmg7gyGrK9hugMDHOeya9Lm03t10HcII2dib5FPfHjpwFwMIATpzBJTows0EbMrD6N1NLrFu3DuvWrcOiRYtw66234u677wYAXH755XjssccAAD/96U/x6quv4k9/+hOWLFmCJUuW4Pvf/34lj20rlLhq+Rq5nweWV1uZO1R9YvJVvV/Supg8rLUegoC8swmEPbHFbm67E7LsTMqJ+1vfBw3AtBl8Ms5kVI7r/ZKcRq8L4zEFqspSHMvJwgWOmTPJGMNYVNZjgTgNXhcGLW494e/K57XmEJOKOsmcyaQz6RIdkG3aM8mdSTv2TK7fcQhu0YFTDm/Fb1/aDSmPuyvFtdaEmYmp2sHxGA7LsW6tFjniiCOwcePGjI/fdddd+v+//vrrcf3110/lsaqK6x7cgld3DePnXzgeH5yfe9BjX7B6YoEAwwBOHhPDjhQjJr0uEZ3NPtsOGBK5IWfSAtwJyldu4GKz0j2TETkOb5pT1+hzgjFgosSTclzg7B+LZgyuRGUVcpxllLkbCyhz8+e0q9kHh5AZWq6qWlbnZJzJZM6kCKfDvs7kCBeTNnMmY0ocr+wawglzm/VNR/lEvKSoujMJwPKbC4Iwsm1gDAfHY7jorlfw6w27cg4xVcsqRU5yAKc6xWRHwIfO5rqCPm9hmx+7Dk3YdgiSyM70sgKKRHcmc7xDZIzpJVmrgs3KNHcxZe6IHNdf1DnGlYoNaWXnycDL3HGV4cBYFHMCyRLSaNr2G06D14mQRZeUP6cz/G40+lwZzqRxeMZIcprbQjRQ4j7ehDOpMpTcwS0FwYj2XNjNmXxt9wiisopTFrXqLnv+MndCTDYmxOQ0HMIhJs/whIT5bfUAA77313ewpS+I/7zgGNS5M1/aqmWVIidQ54YgpC5tsBujERnrnt8JRWXwOh3wuLQYuh2DIXzm+I6CH29Bmx/PvncQ/SNhvbWJqA5ITFqAO425nMmwFNdfQAvOmSz5AI6K9sbM0jKgOXsdgdJcTFWV6WVuANgbjKSIyeT2m7SeyUQ5MyQpGUIzHe5MNte50WQqJpOB40bcFkUNkOpM8v3esqrC4xBzfdqUE7SpM/nCDi0S6NTD25LPXwEDOABwgJxJokDiKkMwIuPDC2bg1guOxbf+8Ab+8sY+bN8/jl9ecoI+tMcZCEbhEKC74XaHpx30DudellFJ/vb2ftz+3E7T21YsbC348RYYJrpJTFYXJCYtEE70QI6E5ayOlVFohi32TOrOZM6eySKcSUmBN30Ah69ULOEQznhMgcqS09np8UDJvdzmwnY8ml9McmeypV4Tk+nN6NmcSaeorUYsNBrIKWo/WyXOYLcWPi7cY4qKuMp04VZpNvQcREu9G0e1N2LXkNY8b2UAx+sS9Rf2g+RMEgUyEpbAmHZtaPS6sO6SE3DH8zvxX397Dxfc8RKeuu4UvScXAAbGomirksByzpGzG/Biz5Du5NsNvhXuni+fhDkBH6JyHDFFhSAAx3UGCn48PtHdMxjCRxPDfER1YL/fTpvBGNPL1nGVZe31M5bAS9MzmShzF9gzyRhDRI7Dl94zaShzlwoeC3RUeyMAoD9tCIc7k+kl94YCgst1ZzLxgpHhTBr2aqfjFh2WooGMuZwuh/Y4ig3jgYKG790u8UBDoRje3juGkxe2wuEQ9DdG+ZxJSVHhFgU9g47K3ESh8DfwLfXaGxJBEPCvpy3E/3z+eAxPSFjz8JspPZQDwQjaq6TEzTmyvRFSXMWuLEshKg1P8Vg6txkLZ/rxgY4mnDC3GUu7m4t6s8uzJp/dln1FL2FPSEzmISLHYewFzhYPZHQmJyTF0h8CFzrprprxY4WWuaW4CpUhc5rblyxzlwrew3f0nCYAWpnbiL6XO0NMWt/PPTwhw+kQ0OBxosnnQiimpAzIZCtzA5pItyLG+XPMo4EArcxtN4KGwHa7xANt6NFK3KccrpW0rG4ekhJOi88tosHrpDI3UTCHQtobED71zDn3uDk4//gOPPveQTywqQ+A9ubmYCiG9iqZ5OYcObsBAPBuYue93dg7EkGTz6W3Lk2Wlno3PntCJza+P4T/+/KekjwmMTWQmMzDRFrJOlvfJHcmGzxOqMzakETMUmh5YaIhKmmPmZkzWfgaw3xwcbNgZj08TkdGmTvrAA7vmbQgJkfCEprr3dpKSF+mCM1W5uYfKzS0nJeS7DZNmN6fahdncsOOVDHpsuhMGst2Mxs85EwSBZN0Jt0Zt9147tGY0+TF9/76DvYMTeDAWBSMoeqcSV71eXdgvMInMWdvMFKyHnzOjZ86Cp3NPvzgiXcpJqiKIDGZB16y5r1d2cQkz5jsbNGiEKwM4SR7JjNdteQ0d2EOWXIvd+o7Rd2ZLGGZmwvo5jo3OgK+jKxJvWcyI2cy0b9psczdUudOeRxjqdsoBNNxOx2WeiajadFAQH4xNNXw/lSOHZxJbYXiISyc6ddfpK3me8qKqt93VqOXnEmiYPi1ON2ZBLTWmh//f8chLMXxzT+8oV+bqs2ZPKy1Hm7RgW377edMxlWG/aNRdDSXVkw2eF1Ye+ESxBQVq/6wxXbXYsIcEpN54KKwKyESswWXG/MQgUxH0ww51zR3kTmTfFuPz536mLqAK2GZm4u6QJ0LcwI+7A1GUsr7+jR3RjRQAWXusITmeu3+TbogNohJgxBMx+10WHJ2Y4ZoIKdoz55J3p/K+5DsEA+082AI+8eiuisJAC7RmhiX0pzJ8ahi6rbGlDi+8+hbti3zEZWDD+O1+DPFJACcvLAVXzp5HjbvGcEPn9wGAGivkoxJjlN04PBZflv+/h8Yi0JRWcmdSUALO//qRxbgzf5R3PaPHSV/fKL0kJjMA3eAOhMiMWvPZOLjXHRaGcKxtpu7MNGQdCaz9EyWocwd8GnOZFiKp/T18Z7JprpsAzi5nyMlrmI0IutlLC4mU53JXGVuseBoIC6GFJv1THIXeFbCIQ+XOHy+GNYbIoE4eiRTDjHJGIMcZ3A7teda34JjsqP7lfeHce/Lvfjrm/tKdm6iNkg6k9mjftacfSQWtNVjS18QQPU5kwCwuL0RB8ZittvRzXvkyyEmAeAbZyzC0XMa8fNne/DanpGyfA2idJCYzAMXhV3N+Z1JQUj+YRVW5s4VWl6YqOHujjdNTPrdTghCUuCVAi5wAnUuvdRhHMIZi8oQBO1rG+FiMhTLLWxHIzIY08roQNLhNCtzp3+/QAFlbmM0kIP3/NnLmeST3O2J3y879EzuGdJiQRYn+roAbaLWJeaOZOLPrdGZBIAD45mlbi4CrDj9xPRiWM+gzR4v5nWJ+O/PLdEd/WrrmQSSQzh2K3Xz1oFSl7k5bqcD//O5JXCJDqz6wxbLKSlEZSAxmYdkmTvhTGbtmZQQ8Ll0oWTJmSxDaHkkywCOIzERPZ5HwBUCL70217n1sHKjmByNyGjwODNyOf0WncmRcGqDvV7mNghivd/R7DkUHRbXKcYhOgS4xNScSTvBh2+4s2KHnkn+O+5PyxHVIpmyP3+8BK6LyRzOJBeTdnBiCXsxNBFDoM6lt6Zk49jOAG781FE4/Yg2PSS/muBv1rbZbAin3M4kABw+qwH/5xNHYs9QGD+jcretITGZB/4iNrPRC5coZF2pyKeO6xOTyhMWXuylXCVavWeyuDJ3nTvTqWv0uUrqTAYjMtxOB7wuh35BMQ7hjEXljFggIOkw5hOTfKiJO5PmZe7JD+DEFFX/fJdNo4F4+wAX7XbYgsPzV+vSXGHtec9+vvT2Dn0/d5ozyRhLOpM2+H4JezE8IZlOcptx6Yfn4e4vnWSboP9CsGs8UH+ZnUnOFz88D41eJ97sHy3r1yEmB4nJPIQS5TW/x4mWejeGckxzt9QZxKQFZ5K7ZmabDYotc+sDOCZl30avq8Q9k5obKwiC3lNqjAcai5hvuNEieIS8MUXp0R+NOXsmzSbiHYhZjAbiYlK0aWj5SJozaYcBnFAsjjq3mOE8u0RHzjYBffBMTE5zA8CBNGeydzis/w6EqcRFpDE8IZlOctcaM/wezGzwYNt++zmTXpej7D8Dh0NAq99ju55RIhUSk3ngorDe7URzndu0Z5IxpjuTfo+Y8nm5KEeZO1vPJKD1KpY0tDwsI5DoV5rV6IUgZPZMpm+/AbS+ugavK29MkR49lF7mNp3mzu5M5guQj8qq/nzpAzg2i6PgziTv+bKFMxlT9DdPRvI5wpJe5k4M4GRxJrkrCZAzSaSiqqwgZ7LaObK9EdsPjNvqurR3JIw5AR8Eofxub3MOI4ewByQm88BLefUeES31btN3R2NRBXGVpTiTBQ3glHKaWzKf5gYSZe6ote08VghGZAQSJWi304FZDV5dTDLGMBaRMzImOX6PM29oue5M6gM4hedMArkni7XHSJa59QEcm4WW8+95Mj2TIxNS1gGyYpiIKag3aadwi46c0UDpAzj1Hif8HmdGz+TrvUH98ahnkjASjMhQWXKVYq2zuL0BMUXF7iF7rFVkjJUlsDwbLfVujIQlqDa7LhNJSEzmQXcmPU4017sxFlUyXiiN+6Pr3QUM4OSY5hYEwfI6QCMR2XwAB9DK3HGVlWR4g29kCRicx47mZHB5KBGybVbmBjSXNN8wUPJ51R7DKTpQ7xZT3NV8ZW4gf4C2VuYWE19De5cdt1nP5EhYgtfl0PtHi5nmvuKezTj3FxtK5mpOSEU6kyZvorQtOJnOZKDOhcNa6xGmaW7CwPCE+SrFWmXxbHttwhmekBCVVb29qdzMqHcjrjJL2cREZSAxmYcJY89k4oXcmKUIJDMmW+qTO0pDFl78JCW13JeOx2ltGtlI7gEca1PUVuAbWQKGWI45AR+GJiRE5bhewjYbwAESYjKfM5k2zQ1opW5zMZndmcz3HMZkVS+TJ0O37fUOOBiWEfC59TcJxbwhODAeRd9wBOte2FmSM03E4qZi0iXmXmOZPs0NAG1pKxVjShzv7BvDkq4A6j2iXiEgCMAQWD5NxOSR7faKB5qKSW4j/Oc8NEFrV+0Kick8TMQUOARNrPBf6PRSt+6gFTiAI8W1lXLZek48TrGIMrf2dc16JrlLWIohHGMsEKfDEA/EBZ9ZzySgbcEZz1NyH5mQ4HE6Ukr2jT5Xapk7IZ5NeyYTayqtOJNe7kzadAAnGJYQqHPpYrKYARzuct/x3E70j4QnfaaJmKK/eTJitWfS6MjPavQiGJZ1x/XdgXFIcTUhJp3kTBIp6IHlWbbf1BoL2vxwiYJtnEk+aFnuSW5Ottdewj6QmMwDL+UJgpD1F9o4dex2OuAWHZbL3J4cGWmTcSaz9UwCpVmpGIxo37Nxu40eXD6SFJONXvOeyQaPE3GV6TmRZgyHte03RrGdISZzlbldFsvcJs6k3TbgaP2pLv3nWkypOqaoaPK5EFNU/PCJbZM6jxJXEVNUUwfcJQq5nUkTR54P4RxMuJNberWNF0u6AqhziwjL8ZL1+hLVz5CF7Te1hEt0YOHMBmyzSTwQjwWaM0Uh8PxNAw3h2BcSk3kIxeK6+8KnitOzJtOnjq2W5SRFNR2+4XhcRfRMZgktBwz7uUvgTBpXKXI6Entv9wUjlsrcAHLGA42YTGs2pQ0RcTHpNXUmrQ3gRA3RQHbcgBNXGUYjWplbdGi9tEU5k0ocS7sD+PhRs/D4WwPYuHOo6DMZ2z/ScTtFXTCakUwxSP6OzmzkE90JMZmY5F7SFUC9W3vjUegbK6J2SY8Nmw4snt2AfaNRvSpUSfZOsTPJK2DkTNoXEpN5CMcU3X1pyfILzcO1+e11bqel9W+xeB4xWUSZW48GMnHq9DJ3CYLLjasUOR0BbeXk3mBEdw+zD+Dwknv2s2QTk3GV6VEx/Ps1G2LSeyZzCPK4qu2J9rhSB3DsFMExHk2slUwMItW5xYKdScZYYmpdxHf+5Si4nQ7c/JetRX+fyZQDEzFpuWcy6UzO0rfgaEM4r/cFcVhrPQJ1btR5iu8TJWqT6VbmBuzVN7l3JALRIWD2FG0U4g40iUn7QmIyD8a+MP5inh6vYpzmBjS3xnAwqVAAACAASURBVGqZO7eYLK7M7XU5MoKkgeQATkl6JhNiMXUAR7uwpJS5i3QmJUXFeExJ6ckEMvdz81gfs75TfZo7nn8bS3IDjv2igUbCvP9Uey58Lq3sWwhynIExze3unlGHK0+Zj237x3H/q71FnSmZv2q2E11IfD3z51BSWOJ+qQM4gOZMDk9I2DMUxpKuQOJrWO9DJqYHh0Kag51+fahl+FpFO2zC2RuMYHajN+8qy1LRwsvcIRKTdoXEZB5ChmBm/u4ovW9jOCxBdAh6f2C9R7SUMynHVVNHjVOMmAxLimm/JGB9jaEVzMrcDV4XGr1ObQAnmn8AJ9dZ0vdyc5rS+j6N22vSsRL8zp1N3nPJ163ZyZnke7mbE8Ld5xYRLdClS8/j/NfTF2B2oxf/9b/bi8qeDMVyO5NA9vYCs2numQ18C04UbxhK3IDm9APkTBJJhickNHidOd+M1xpHJuKB7LAJZyozJoFkBFS2dcZE5Zk+f4lFwJiWychfzLgLl9EzOSGhuS45KFJfKmfSJRaxm1vNKibThdhkMCtzA0BHc11impv3TGYJLffmDncfNkzIG2nypQaXa8Mz5t+vlWig9J7L5AYc+ziTunDnYtIlIiwX9oYgfVCpzu3Ev3/ySIxGZKx9envBZ+JtHNmigYDsfadm+aqzDD2Tr6eJyXq+VYrigYgE02WVopG2Bg9a/W68W2ExORFTEAzLU9YvCWjpJHVukQZwbAyJyRzEFBWKyvQViV6XiHq3mNkzGZbQUp8UVX6PExNSPG9afznK3FEpDq9J6REoLBqIMZZzi4lZNBCgxQPtH43qblqu0HIge5l7RG+wT/389P3cxu016XDhlGuaO+nYpUYDyTaa5uaT8wG9J1fUB62sYpbHee5xc3DErAY88dZAwWfiwi5bNBCQ/XnXnUnDWfweJ3wuEYPjMWzpC8LtdOhlPd2ZpHggIsHQNFqlaGRxeyPe2z+GeAXbcKY6Y5KjbaCjnEm7QmIyBxMmpbzmxFonI9yZ5PD75+trkyyUuRWVFVRyjcjxrM4kdwOtDODc/eJunHjL37MKz2BEhtvpyJii7gh4oagMOwZDEB2CaXQMkIwMylbmHk6bkOcUUubOJ2oA6NFE+jQ334BjI2dyZIK3FGjfu9cl6nmiVjHL4xQEAZ3NvqLaHvjfhnk0EHcm85W5k32ugiBgVqMHB0a1MvfRcxr1n58VZ5Ixhm//8Y2ie0CJ6oExlhjOmx6xQEaOnN2AqKxiTwXXKvItZ1PpTAIJMUk9k7aFxGQOeI+WUUy21Lv1F3dAmwYORuSUd8l+i8Hl+Z3JhLNWoJjMJuBEhwC/x2nJmdy0exijERk9gyHT2/kqxfTBF36Bee/AOBq9zqyB7PmmuXVnMqPMbeZMZilzi1bK3Kkiy44DOMHE98qFdZ1bRKTA3MVseZwNXiekuFpwagD/3TZzJvOtsZQSQj39jdTMBi92DI5jNCLj+K5m/ePJnsnsf08xRcUfX+vHq7uGC/guiGpkLKJAUdm0K3MD9uib7K+gMzk0IVHerE0hMZkDfcjAnSomjWXu0QiPbTE6k2LK52dDUvI4k6780TbpRKW46fYbTqPXmTOOh9M7rG1I6Rs235QSDMumk5Q8HkhKBGRnQ187mc2ZnEgVUJz04HVj4Hg6VpxJvWdSL3MXN4ATlhQ8+dZAWS50o7w/1ZfsmVRZ/jWRRrKtnfTncYizMWHyRovjyiPiZT1nMvUsbY0ecA2/pDugf5y/OcoVt8XP35AlJJ+oHfhKvekUC8Sxw0R3JZ3JmKLSIJ5NITGZg2SZOynOWurciMhxPedv2MRB08vcOV78GGP6OsVsWJlGTn/McI4yN6CJsfE8AziMMfQO5RGTETll+w3HeIHJFgsEWOiZzDfNnRAPVqa5pRyuGxfq6c6kUqAz+aOn3sPVv/+nPjxSSvRoIH2aW3vuogUMZ+ll7rTnijvE2UR9Nsz+Njj8dzpbmTu5kz71LLMakpl1x3dlislcziT/PTJzSonaYjoGlnMWzKyH0yFg0+7hivVNVqpncgatVLQ1JCZzYOa+cKeM9/Slb78Bkk5mLmdSyuLOGOElSaslSDnOEFeZ6fYbTqPXlbfMPRqRMZ44e99wJON2VWV6mTsdnjXJv1Y26t1OCEKOnskJ82nxgsrcRUQD8Z7JXMNH6QyFYnhgk9arNxCMWv48qwQjMurcon5G/mahkHfoujOZ9kaDi69Cncmc0UBWB3DSy9yJie4Z9W50Gt6UJPfdZ/9++XkacvzOEbXB0DQMLOd4nCI+sqgNL78/jEt/84q+fnQq2TsSRqvfnbMCVg5aKLjc1pCYzIHZAA5/N8x7+oZNpo7rLfRMmsWjpFOoM5lrLzenud6FkbCcUyz1GtzIXhNnMiQpUFmm0AOA1nqPLiayxQIBgCPRvzkey+5M+j3ODKHoSew+19oLGKKJkHYzCipzu1LXKRYSDfS7jXv0QR4eplxK0oU7d+oKWamYrcytO8RZfg7ZMPvb4OQbwNHfSGX0TGovFku6Aim9ttacyUQPJ5W5ax4eXD0dB3AA4BcXL8XFH+zGiz1D+OTP1uPl94tfi1oMe4MRzJliVxIgZ9LukJjMQXLIIClo0neEjpjkIfotTJ/qYjLPbm7Aes8kd9lyOZPzZtQjrrKs5WsA2DOUW0wGJ8xjgQBNJPLyR66eSQBo8DhzOpPN9ZmfLwgCGn0ujEVkKCqDyjKHSjhWBpjSo4H0nEmL0UBhScE9G3frm2DKIyZlPRYISP58C1mpmP59chom2TNZZ/LGJa8zabIBB4D+ArV0bnPKx/UNODm+X71nksrcNQ+Ph5mOAziAlubw/c8cg59+fgkmYgouuvNl/OLZnrxRdKVAUlQMjsemvMQNJI0cypq0JyQmc5CMPzE6k6nB5cMmvX3crZnqMjcXF7mcyflt9QCA9w9mj5bgArK9yYuB0UiGw8RzD816JoFkqTtXmRvQSpLZevVGJqSMSW5Ok8+J0Yic1W3jWNuAk/oYfANOtsDtdB7c1IdgWMbVpy0AUB4xORKWUlxg/vMtyJmUswzgeIrvmfS6HKbr1NwJQZ5NxPP1lsZoIAA4aV4LfvTZY/HF5fNSPq7v5s7x95Qsc5OYrHWGpnHPpJHzlnTgsa+twMKZfvz4b+/h+kffKvvXHBiNgLGp75cEDC1mNsyaZIyhfyS7QTMdsCQm161bh9HRUQDANddcgxNPPBEvvPBCWQ9mB7gT4k8pcydWKoayO5MFlblLOIATtiQm/QCA9w+ZR/4AyaGb5QtaoTJgXzC1b9JslaIRfqHJNYADaC/82SbLh8NSxiQ3p9Hn0sSkSXaiEWtlbu0xeP+PPoBjoWdSjqu4a/0uzKh34ysr5sMlCjg4Xtp3zUpcxXg0dUe5zz2ZnknzMreV9Z9GjDvr07HqTLrSfvcdDgEXntiV8bhu0QGnQ8jjTJZ+AGe6XvfsznQewEln4Uw//nyNJiiffmew7F+vUpPcQNKJtqMz+fhbA1jxn8/qq2CnI5bE5C9+8Qs0NTXhxRdfxNtvv43vf//7+Na3vlXus1Ucs2DmDGcyUfI1y5kM5RgY4C+0Hks9kxadSQtl7vmt1pxJv8eJ47qaAGQO4Yyk7YpOh8cDNeZxiRq8TtNp7ogUR1RWcziT2hBRtuxEjhUxnu5u8mggKzmTf31zH/YGI7hs+Tz43CJm1HtK7kzyQaMmM2eyBGXu5ABOoT2TcdN+SSD/OkU5S89kNgRBC7/P1TMZKkPP5HS97tmd4QkJ9W5xygdA7IrPLWLejDqMRsqfwVipSW4AaEkMXNkxuPzpdw4AAPpHMgdWpwuWruZOp3aBfuaZZ3DppZfirLPOgqLU/p5cs2DmjJ7JsAS305EiOK04kzFLPZOJMneBPZO5LrIt9W40ep14/1BuMdnVUoeuljr9v42YCRwjfBK3KYsY5DR4XYgpaoaDlW37DafJ50JUVvU+ucltwEl1N3mZO98GHMYY1j3/PurcIi758FwAQGuDu+TTlXpguVFM6gM41v8G8w/gFNozqaS0fxjhIpGXs9ORskxz50Lbd59/mjtfa0UhTNfrnt0ZCkm6sCA0mnxuyHFW9gxGXUxWwJls8DjhEoWMDXSVhjGGl3ZqQ1D5FpXUMpau5g6HAw8++CAefPBBfOxjHwMASJK9fqDlgDuLRgemyeeCIBidSa23zzh96rewscNaz2SB09wWytyCIGB+mz+rMynHVewLRtDd4kN3FjGZr8z9yWPa8e2zjsCZi2flPK8/S4l1JE8ZiwsGLtzyb8DJkTOZFlouCAJcopB3AOe59w5i2/5xfOGkbn04ptWvOZOldAeCemC5ocytO5MFhJbL2crc2nNZ8ABOTEkZTDPiyiPiJUWF6BB04W6FfM4kb5coZZl7ul737M7whIQZ03SSOxu8pzqYJ0N4svAyd2ei+jSVCIKgb8GxEzsGQ/prUaHtQrWE5TL3/fffjyuuuALz5s3D9u3bcfrpp5f7bBWHv3gZxZlTdKDJ50pxJtMdtDp9A052ESMXFA1UWJk72zpFzvy2ehwKxUzzJvcFI1AZ0N1Sp5cy+kbMxaTZtDWgOWfXnL4wZ7kdyB5cPmzSh2qET4kPjmuZjtl6Jp2iA6JDyN0zaSKynA5H3gGcO57fCadDwFdWHKZ/rNXvQUxRS3pB0YV7XWY0UC5xlU6+MnfhAzjZy9we3ZnMXuZOH77Jh1VnstRl7ul43bMzjLGEmCRn0giPDguW2bXbG4zA73HmjH0rJy31HttFA73Yc0j//4Vck2sNS78RH/zgB/Hoo48iFtPU96JFi3DbbbeV9WDlJiLFccOf38ab/UE8fu0ppiW3UExBvVuEI81BMe7nHp6QUgKWAa1853Y6cg/g6M5kdsGVnOYuzJn05hFxC/gQzsEJLDFsGgGSLmT3jHp4XSJmN3ozYoTM3LJiaMziiiW335iL1aSY5M5kdkHuFh05o4GiJiLLmceZ/GfvCF7dNYwLlnam5K21+jW35FBIKll4dlJMJp9r3sZQ0AacLGVut9MBj9NRkACOqwwROZ6yZjT9MYFcoeWsoBI3oAnovpw9kzK8LkfBj5uLWrzuVTvjMQVSXKXhmzR0ZzJcZmcyGEFHwJdSiZtKWupd6M8Ra1cJeIkbyG0g1TqWrrxvvvkmPvCBD2DBAi3+5LXXXsPq1avLerBy0jccxgV3vISHXuvH9gMh7B8131oyEVNM3ZeWOs1ql00mbTl+jzN3NFAZciathJYDxiGczIluXUwmStxdLb5MMRmR4XY6soaFW4U7k+kOaT5nkr8rzlfmBrTnN9fzZxaZ4xJzO5PPbdOmJnmvJKc10cdVyiEcs2Gn4kLLzdcpAtkHobLB81PNVikC1kLLc70BMKPe7cybM8ljjkpFrV33agE+fEE9k6nwN5vlFJOqyjAQjFakX5LTUu/BeEyxXK0rN0pcxcvvD6GrRXtOprMzaemKfu211+KXv/wl2traAABLly7F448/XtaDlYsNOw7hUz/fgHcGxnDk7AYA2V/8w5J5Ka+53o2RsJSzt6/eI5YwGqjAae48YvKwRNbkLpMhnEwxWYeRsJwiNvhGlsm+O822yi9fz2S6M5lL1LqduZ1JM5HldAg5o4H489ya9oLWltjgcqiEQzh82ClgMoBTUDSQXs7P/N3w5wiPNyPX9hsgvzMpKWrhzqTHCUlRswrUUEwpecZkLV33agV9lSI5kykkeybLVwI+GIpBiqsVmeTmzNA30JXXgbXK2/vGMB5VcObi2QCoZzIvoVAIK1as0P9bEAS43dX1x6xN3+7Epb95BfE4w6+/eCIuP2U+AK0saUYoppi6Ly11bsRVhj0J4WXmoNW7ndbK3Dl6x4otc+frVZw3ox6CYB4P1DsUhiAkox+6mjVRaYwHCoblrK5hIfBScHq/Xr5pbp5feWAs0TOZs1XAkWeaW4XTIaSEb7tEB5Qc0UDZ3ggky9yFicnRiIzLf7cZb/WPZtzGnckmQ0tBncuZOPvky9yA9nMoTExm5q8ayedMaj2ThTqTuQX0eLT0YrLc170dO3Zg+fLlWLRoEZYtW4atW7ea3u/Xv/41Dj/8cCxYsABXXHEFZNkeL6SVIJkxSQM4RnjLUTmdSR57U4lVipzkFhx7BJfzfskzFs+EIABhKnPnxul0QpZl3Ynq6+uDKFZPxhdjDKv+8AZ++OQ2zG/z489fOxkfWzwrb1lyIqaY9oVxkbNzUCsTmzloWpk7/xRxKae5oxadSa9LREfAh51Zytxzmnz6ucwmuoMROWssUCFkG8Dh7zoDWULPea/lISs9k05HnmnueMbnO0UhZ5lbF2ZpfwNcTB4sMAft6XcO4O/vHsD9m3ozbjMbwPG6tfMWFloeh0NI5mgaydeSkY7ZZigjnrw9k4UP4NTlSUgYj8olneQGyn/du+qqq3DllVdi+/btWLNmDS677LKM++zatQs33HAD1q9fj56eHhw4cAC/+tWvSnaGamO6r1LMBr8+jJZxmruSsUCcFps5ky/tPASP04Glc5tR5xJzrlCudSyJya997Wv49Kc/jYMHD+I73/kOTj311KrqHRIEAUu6Ajj76Nl49JqT9S0wupOUpSw5kaXMzS9kXIyZOWja9KmFMneOF6dkz6Q10RC26EwC2iac3UMTKftcGWPoHQrr/R8A9KxJ3jepqkwvc08WXhb+05Z9GDW8ox6ekNDkc5mu6gNMBnBylbnF3M5kTFEzcjnFPGVuPXDele5MFtczuWHHQQDAq7uGM24LhjWRZHTy3Ikp9cJCy1V4nKJpa0KDVxOTViONzHbWG3GJud8EyXGWc/DMDF4hMJvoZowhlGMjT7GU87o3ODiIzZs3Y+XKlQCACy64AH19fejp6Um530MPPYRzzz0Xs2fPhiAI+OpXv4r777+/JGcwwv/2rWx+qiS0StGcprryT3O/mdjusiDRJlUJZtjImYzKcWzePYJl81rgdYl5X/NrHUticuXKlbj++utx0UUXQZIk3HvvvbjwwgvLfbaScumH5+KOlUtTXnD0HjeTF385roVpZ+uZBIAe7kxmGcCJyHHEs5RLrfVMFljmtuhMAtoQTlRWMTCWHD4ajcgYjym6GwkknUkeDxSSFKgs+3BMISxo8+OqU+fjjb4gPverjRhMnGUkLOV8seAXTu6m5Sxzu8S8YjLdmXQ5cpe5Y1linZrr3BAdQkE9k4wxbOjRpgF7BkMYSvtdDEZS93ID2psjn0sseDd3NtHt9zr1CW0r8EGYfD2TWQdwFDVne4cZuZzJiByHylCyCXpOOa97fX19aG9v14PRBUFAd3c3entT3ene3l7MnZsc9Jo3b17GfUrBbc/04NQfP4tt+8dL/tilRB/AITGZQoPHCdEhlK3MzRjDk2/vR0fAh6PaG8vyNayQ3M9d+Xigf+4ZQUxRsXzhDAD548tqHctv5ZcvX47ly5eX8yxlxcyR4Rcks57JXO4Lj6zZmeg5NMtb5E5KWFJMX+TKElpuYQMOZ35bcqKb90emD98AwMwGD9xOh35bcCKz7DoZ/v2TixGo+3/snXmYG2eV7t9atLSkXtybl17tdjvx3nbsxM7meEyAGxgnPEMwMxAcJgvxXMiA4T4sAzdw75BZyJiZQLiBTMjMxOAxhBAyhC1khcQQe2zHWxK3l97stt3d7k3q1l73j9JX2qpKJakklbrP73n8JN1SS19p+erUe855jx3/8Ku38cFH92H3Xdeo2i0l4rGL4DiACWl6aW6HwOuPUwxF0ppSRIHDlF4HeFhO06ZaRvE8hzq3PStl8p2Lkxj2BlDvsWPYG8T+nlG8d8U85fZRX0j181VhF7Iep6j1OlUmeE1qpa4TyZTmZilszQacXGomHdo1k6ze0+yaSaC4+54RZVjvPrt27cKuXbuUn73e9DIWLViAcLh/DCuaqg3/XbFhQUQddXMnwXEcqitsBTMtf3NgHOfGpnHPDQtLZgsExJVJKwSTzBLouo56APIeVWhrJitjaPfdvHmz6gfoxRdfNG0h3d3d2L59O4aHh1FdXY1/+7d/w/Lly017fDVsAo85LhuGVE7+TH1RO2EyVY6pdWpXyezvfIGIejBp0LTcLvKGx0f5gxHYRd7QZJFF9XGvyRs65W5VFjC2JASTPM+heU7cHoh1C5pRM8nYcVMHalw2fOmnR/HBR1/HZV8Qq5q1T2g8z6HKaVPqg3TT3BkacNSUSVHgEdatmYxovm/yFBzjG93vu+UC7vs2deBvn3sLb5y9nBRMjk+HlMA/kQqbkGXNZFRTwWWfzwl/GI0GRAfFIDxHZTKXBhw9ZbJQwWQh972WlhYMDg4iHA5DFEVIkoT+/n60trYm3a+1tRWnT59Wfu7t7U27D2Pnzp3YuXOn8nNzc7Ph9XS1yn6zh/vH8NENbRnunTssGM41IBn2BVFhEwxd9Mw2aipsSeVCZvLLo4MA5OlmpSTegFP6YPK108OodIrKxZfLLioTgmYjhr6Rn/vc55T/9/v9+OEPf4glS5aYuhBWjH7nnXfiqaeewp133on9+/eb+hxqsBF4qejZn7APNBMJtHwmAW2rACNpbo7j0FrrQt+IMZPW6VAk4/QbRqIyyWDBZFtdcvDSWuvC66dHYvWS+qMUc+XPr25FdYUNf/2fhxCOShnT6NUVCcFkhm7ugJ41UCiidIczbLx+A04wHFW12AGA+koHeka0556n8tqpYdgEDn9+dSu+8/Jp7O+J100GY9N0qlXqU112Ietubi1lUmuspRZTGXwm47O5NYLJcFT3c68G+1yrpZEyBbe5Ush9r7GxEWvXrsXu3btx55134ic/+Qmam5uxePHipPv92Z/9Ga6//np89atfxdy5c/Hoo4/iwx/+sClrSKTe40DznAocjtXFFYKpYBibvvEy7ry2Hf9z8+LMf6DCZV+AUtwaVLtsBQlmJEnCc0cHsaDamTbkotjUuOzguHi5Q6mY9IdwZGAcf3JloyLeeBz6XrgzHUO77/ve976kn2+99da03+UDK0b/zW9+A0DeQD/5yU/i1KlTaZur2dR7HDh+Pt2Sxaub5k6wabELqmllFoRqFeSyE20m8+b2OhdefmcI4UhUsyGFMRWMGKqXBIB5VU5U2AScSfCa7FdJcwOyPVAwPIQhb0DVRNssblk5H5VOEX/1g4MZU22J47wydXMHw1FIkqSqhvhVlUn9CTiBcFRHmbRjKhjRNLxPJBiO4o9nL2Nt6xy4HSLWt8/B8ycuYtIfQmWC8qoWWDtt2aVUAuGIalAKxBU9oyMVvRmsgdh882BYo144J2VSe4QkcwMwc5QiUPh977vf/S7uvPNOPPjgg6iqqsITTzwBALj77ruxdetWbN26FYsWLcLXvvY1XHfddQCAm266CZ/4xCdMW0MiXS01eO7oICb8IcUxwUzeuTCJockAXnr7Uu7BpDeI+kqyBVKjpsKG4+cnNPe6XDl2bgIDo9O46/rSprgBuTlyjsuu2MeVij+euYxIVMJ1HXXK71x2QfHCNXMSV7mQ0+4biURw9uxZ0xahV4xe8GCy0oEJfzhWUxYPxJhflFo6Re6uldUrLQXNo3Sf5q5MArInZDh6CefH/Gitc+ne1x8yHkzyPIf2eneS12Tf5Sl4HGJaoJhoD8QCHDPT3Inc0NmAQ1+5OWPgnBgYZQomATZ1Jf21CYTSawltGdLcsjKp/pwNCV6TmYLJQ32jmApGcEOnXHNz9cI6/Pr4RRzsG8OmJQ0Yj5UUqNWnuuw5NOBUqn824ubxxoJTpWZS5/i0xlhKkhTr5s5+Nrf83CrKpJLmLsxnkmH2vnfFFVdg3759ab//13/916Sf77nnHtxzzz2mPa8WXS01+PmRQRzpH8f1sc+kmXRflLMgJwYnEIlKhspxEpEkCSO+IJbEhk0Qycxx2REMR+EPRQ05ehjlOSXFPS/DPYtDrdte8prJ107L5UnXLY5/T9g+OhWIoNpFwaQqH/jAB5QrkkgkgiNHjuCWW24p6MLUCs3zKTDXgtm5jHiDSWaseqkzjpOvji5Naqdc3BnS3Fodwam0xUYf9oz4MgaTevOS1VjU4MYvjg7CH4rAaRPQd3kKLbWutKtPZhXUNzJVsDR3IpkCSSAlmNQJoBM9D1WDSRVrIJHnNOv9ADkw1auZBORgMrVcIJXfn0rekK5urwUAvHF2BJuWNGBUZS43o8KWbQOOdppb8fs0mOZWmtN0Pms2kUdIpVaVlQ+YqkzG1lNpcpq7FPteKVmj1E2OFiSYPHlR7hSfCkZwdtiLxY3ZBYVTwQgCYZrLrUV1whScCrs5XpByF/cg5lU5saZljimPmS+1LjtOqXgkF5PXT42gsdKBxY0e5XdMePIGwwUTW6yMod33tttui/+BKOKLX/wiNmzYYNoijBaj51NgrkXiyT8xmMw0Mq7WLQeTWlNaFCVFw8TUuDIpB5A9Iz7ciAbd+04HI8rxGKGj3g1JkscqLm704PyYH+9amt6B0ZJgDzQxLR+PWodxMUlMwzmNKJMqgU04EkU4Kqk34OhaA2mnjOsr5c/D0GTmK+ffn5ILuFfGUvpL51fC4xCx/+wogATDcpXnqrALCEaihsof2Jq1g0n58Y1OwVGa0zRqJgFtZTIUMXYRlQq7SFKrSWLrNjvNXeh9z2osX1ANkecKVjd58lI8ADh6bjzrYPIyjVLUJXEKzvxqc4LJE4MT6B2Zwp3Xtqe5V5SKWrcdo73BnNRtMxiaDOCdi5O4rWtBkvDCspFTs9Rr0tDuu3379oIuwmgxeiFo0BiBpzQZaKQLWHq7VuMKJN6Ao64ehQxYAwFymhsAeoYzN+Fk04ADQDFvPzPkg8suIBKV0uolgXgw2Xd5CojFWIVUJo3AgjkhZRRiKswUXs0eKD5iMPk1swkcIlFJs/YoqFszaWyk4oQ/hDf7x3DzsrnK+kWBx1Vtc7DvLzPSTAAAIABJREFU9Aj8oYhSn6qW5mblDNOhCCozBGaSJOl2c3sc2dVM+gJh2EVeV120aZjFs9/Zsm3A0dmovQXq5i70vmc1nDYBS+dX4XD/mOl1dwDQfXESNS4bxqZCODowgQ+sye7vR2iUoi7KfG4TO7p/YZEu7kRqPXZIkmzQXpeFeGIWr8dS3Nd2JKv3rgzZyJmO7u6bqAKqkZhyzhetYvRCw5Sk4RQliQWBmspkLD2uloJM/Dutq5RgOApOY7xdIgtqKmATOPQa6BBm6WqjLIyl0M8Oe5WGFrVgssppQ43LhoHL0/A4RdhFHk4dO55iwDqwMzUwsdpGtcCGBZOpxyLwzNpGvbZPLzAzGkzuOz2CqARcvzh5Q7p6YS1eOTmEIwPjis2H2meMXTRMh9StpxIJRSRIkraFktZYSy2MTJtxiLxqqUAhlElvINaAY1Kau5j7ntXoaqnB0XPjGBidTrIIy5fx6RAGx/24tWsBXnzrEo6dS296zASNUtQnPlLRnHpCSZLwi6MX0FjpwLo2a6S4gWSvyWIHk5Ik4fHfn4Vd5HHTlcmZQuWcP0s7unV33+rq4pnXahWjF5r4POXkk3/GNDdTJjU2towNOLG6u0xX/wLPoaXWhbMZgslQJIpQRDLcgAMk2gP5lHS91gmkZY5Lnttd48Qcl63kXX1Gg0k9mxpmrZOmTMYC/HA0CrvKkKigjrWN0WDytVPpBdyAHEwCct0k25TUlEknCyYNbFxsNrm2aXnyRKFMTAXDmrZADC1/T/Y+ZDubm32udX0mHeaUXhRz37MaXS01ePIPvTjcP2ZqMHnqklwvecW8Slyc8OPowDiiUSmr1OkITb/RhWVrzFIm374wibPDPnxsY5tlUtxAstdkZ5Gf+4W3LuHIwDjuvLYdjZXOpNtYFpOUSRUeeOCBYq2jZGid/H0ZvPRY8KVVM6kU42qkufVSpaksrHPj1e4h3RqRbEYpMiqdNjRUOnB62IeGKvl1UFMm2e+PnhuHwHMlT3ED8Y1Tz2MSiJcRBFQm2ihpbltqzaT8Gmt5Teo1s9S67eC5dKU7ld93D6OppkJRhxmrmqthF3n88exl5WSu5hiQmObOhFY6n8E+48YbcDI3etkEXnVTZa9ptj6TPM/BZRdUu7nZus2qmZwN+54WrAnnUN8Y/nT1AtMe92Ssk3tJYyVGfUH84cxlnB3xoaPBk+Ev4yhpbpp+owrLYJg1BccqRuWpsGBytMgd3ZIk4Zu/PQmHyGPHTR1pt8eVSQomdXnjjTdw+PBh+P3xWc73339/QRZVTNhYrtSpJZmVSVvsvxkacHTS3EZPqG11boTevoTzY9qpJ39MocrWEmJRvRsnBiewoNoJjgOaNMYYsuc9NzatqGelRAkmM6TblW7uSHoQwhQ7Z0qQxWoY1eaqhyNRRKKS5nsn8BxqM4xUPDc2jTPDPnxoXXOawusQBXS11OBg76gSMFapBEnx7uZsgkn1NYsCD5ddMNyA4w2EUevWV600lUlWM5mDD5vLLqpu1PIYSKEgxfgzdd/TYmG9G9UVNhzuHzX1cVkn95K5lcqF+rFz41kFk6wBp55qJlWpMVmZ/MWxC6j32LG+vfT7fSJ6U3CePjgAjgM+sCb/5txUfnPiIo6fn8Bd1y/E3Cpn2u3uDALSTMfQjv7ggw/ivvvuw5e+9CW88sor+PKXv4yXXnqp0GsrCg5RQJVTxPBkapo7VjOpocDcuKQB1yysxfqF6rUkiuStcZUSiBgPJhfWxzu6tVCUyWyDyQYPJv1hHO4fw4LqCk31itkDAerdxcWGBVgZayZ15pv7Q+rKpJLmVkmNx83mtV9nralKjNdiIxSv71Tvzr9mYS18wQje6LmMKqeo2mDEAk2/kWCSpfN1Am+PQ4Q3C5/JTB6asg+rds1kLsGk26GhTPpDpk+/AWb2vqcFx3FY3VKDY+cndMeQZkv3RS8qbAKa51QoAwmODmRXN6mkuUmZVCXegJO/Ynfy4iROXfLiPcvnlaRjWo9ajfnc/lAEX37mGB769UnTnzMalfDN50/CaeNx36Z0VRKIZ3hmaze3oR39hz/8IV5//XU0NzfjJz/5Cfbv3w+7feZ8oesrHRjxpae5nTbtOdeLGjzY+4mNaXUTDFGQm1TMUiYBoEdnrGIuaW5AViYBYHDcnxQwppKY/s406rAYZJ3mVmvACanXErLgLaSiTBqxdGqo1J/Pzfwlr02YnpAIUwLGpkKaDV4VsYscM9LcgNyEY0SZjEYlTAUjqpOhErGLgm7NZLYNOICOMhkIm97JDcz8fU+LrpYaBMNRvH1hwrTHfOfiJDrnesDzHBbWueFxiDimMnlMj8u+AOwir+mwMdupdNrAceYok8+fuAjAeiluAKiLKdOpweRrp4YxFYxgyBtQ9anOh18dv4C3L0xi+8Z2NGhMYMqUjZzpGNrRnU4nnE4nolF5LN0VV1yBnp6eAi+teMhKUnqaO1+1w+MQdYLJiOETatweSFuZZOnOrIPJhnjNnla9JCA34DDUGkKKTbXRBhwdn0mtIIvVTKopk5lSxoD8efIGwqqzs6NRCa+dGsbS+VWanqBr2+YoFzFaYyvjDSn5p7kBwOO0GSocnwppT4ZKxC6ozzcPGQjGtXDbBU2fSU8Bpt/M9H1PizUtzLzcHL/JsakghiYD6Iz5SvI8h2ULqnD83ASiOn6uqVz2BVHntpe8+c+qCDyH6gobxkzo5j7UNwaB53CVhbq4GczjODXN/evjFwDIe/2EwZIdI0RiqqTLLuDeGxdp3i/uLU1pbk1cLhdCoRC6urrw2c9+Ft/85jcxNZXZ97BcaPA4MDoVTAoevIFIxlReJtwOUbsBJ4s094IaZ0Z7IJbudOaQ5mboBZMLairARForuPtXGa6ZlF+PbKyBbAnWQKkEDQWTzLg8PdX9zsVJjPiCyghFNTwOESsWyObx1RrKZKI1UCa0FNhEqpyiIZ9Jn85kqETsomxanqoQBPNIc7scoqbPpNnTb4CZv+9psZoFk33mBJNK883c+F6zYkE1JgNh9F429npeGPfjrQuTyoU1oU5Nhc0UZfL4+XF0NnqysporFg5RQKVDVKyiAPnC/7dvXVJ+Vtt7c+W5o4PovuTF9mvbda2I2EQwUiZV+NnPfoZIJILvfOc7CAaD+Kd/+idMTk7itddew5NPPlmsNRac+pgJaqJsPhUMZ1RfMuHWSMsB2aW5RYFHyxxXQdLcLXMqFJsWPSsQu8grUxWskOa2CTwW1bsznlyY+qteM6luDaQok1G1AFT+G733TstuCgAO9FwGINdF6sFS3Vr1qWyTN2YNxGpDtT8bHocIbzCcUSmKN6bpf85sGpZMoRytgQBZmZwKRZLWGI1K8AbNTXPPln1Pi1q3HW11LtOUycTmG8bKZvli6ahBv8lHXzmNYDiKe25caMqaZirVLjvG8+zmHvYGMDjuVyZzWZFaj12poQWAA72juOwL6l7I50IkKuFffnsSbruAe2/QViWB+GAFUiZVeOCBB9Dc3Iwnn3wSAwMDaGhowGOPPYannnoKXV1dxVpjwVE7+ctp7vyuyvTT3MatgQCgrc6FvpEp1Q5jIB5MZjMBB5ADVaZI6imTQLwJxwoNOADw8/uvxwN/ulz3PsbS3CnKZOx9Casok0bT3ADSmroA4FDsBN0VU3+0YB3zWmnurJRJI2luhwhJ0h7/yfBlMPNnsM92qrobDOdmDQTIqXVJAvzh+DH7gmFIknmG5cDs2ff06GqpwZlhn2Kcnw/dsWCyM0GZZIGKEfPyixN+/PCNPqxqrsbmKxrzXs9MxgxlkgX4K5stHEy67cqEMCCe4t62vgWA+oV8Lvz8yHmcHvLhL69fqGkDyLAJPOyidp/ETEd3Rz98+DD+67/+C16vF9deey1uuOEG/Pu///uMS/XUVzKvyfiH0xeI5K1MuhyCZh1aNsokALTXuxGMRDE4Pq16+3SONZNAPNWdMZiM1U1aIc0NyMFFptdQsQYKpwddijVQymvG6hXVupEDBmr+1D5PDNkMuiLj5IaNHXXoaqlJMzVnVGRlDaSuwCZidD43+zxn8pnUCuJznYADxNXQxI5ur8kek8Ds2ff0YBc7hwfyVydPXvTCbRfQVBNv8FtY74HLLhgKJv/fy7Iqef+fdFK9ZAZqXDZMhyKq9dpGORbrsl9hZWXSZcdlXxCSJI+9/c3xi2iqqcCmJfLFhlnK5DOHzkHkOfzldcYUcbddoGBSi3Xr1uGRRx7B+fPnsWPHDuzevRtNTU249957i7G+opCqJEWiEqZDkbzVDrdDhD8U1bSYydQ8kghL5/ZqpLqZQpVLjcsnNy/GA3+6LGOAs3xBFXguuRnH6uh1c2tZA4nKBBy9mkk9ayDmXZq8oY1Ph3BmyIeulsxF7ZVOG575n9fh3cvnqd6uWAMZqpk00oBjbK6sMrM+ozUQUyaTX/d8fSYT1wAkzuU29wJnNux7enSZWDd58uIkOudWJgWCAs9h2fwqHDs3rtt5e2nCjz1v9GFFUxW2LCVVMhMsazSRR6qbDadYNr/KrGWZTq3bjlBEwmQgjOPnJ3BubBrvWT4PjbELeTOCyalgGK+dHsHVC2szqpIMt0PMmN2ZqRje0R0OBz70oQ9hx44d6OjowH/+538Wcl1FJfXkn2n6jVE8GvOEo1EpNvc5uzQ3AJzV6OieztG0HJAL7j9u4Mrroxva8NLnbjJ1zFqh0U9zqzemaAVCiY+jaw2kMVXpyICxFLcR4spk5o1La9JPIlXKfG5jymSmEhBHBmXSlmM3N5CsTE4ooxTNb8ABZva+p8eyBVWwC3ze5uUj3gBGfMGk5hvGiqZqTPjD6NNpwnn0lTMIhKP46y1LSJU0QLUJU3COnbNu8w2DeY1e9gbxm1iK+z3L5yq2PWYEk7/vHkYwHMWWpXMN/43bLmKKTMu1OXr0KD7zmc+gqakJDz30EO677z6cP3++0GsrGqkjFdmHIe8GHA3fqVy89tjYPa2O7lwbcLJBFHjF87JcULq51VLWIdbNrd6Ao1afaqT+sNZtB8elB5NM5TEjmFRqJoOZjaWNpLmZCj+ZwbjcZ/C7UYgGHJfKuDIW3BbCZ3Km73t6OEQBSxdU4XD/WF6effFO7sq021jdpFYTzqUJP37wx14sX1CFd5EqaQimTOY6anDEG8D5cT+WL7BuihsA6hKm4Pz6+EXUue1Y114Lt0OEyy6YUjP5Qqw7PJvPnluntG2mo7sDf+c738H3v/99DAwM4KMf/SheeuklLFu2rFhrKxoNKTVuXoP2J5nwKDVeyR8upW4sC3WmqaYCIs/h7LB+mjvbBpyZjpZCBsQbOdKUSV67AceIMikKPOa47GnzuQ/3j0HkOSxfkH/6iI2AnA5loUyakObONGaUoaUIB2OvaTYlHgxFmUxQ+lma28yaydmy72ViTUsN3uwfw7+80I2GSgc8DhFuu4g5bjvWtNSANzAZpfsSa75RCSabWRPOBN6/Kn0O+HdfZaok1UoaRZmCk6Myeey8bFS/ssm6KW4AqI0Zlx/qG8U7FyexbV2LUuveUOnIW5mMRiW88PYlLG70ZCWguB2ioTr2mYjuDvzzn/8cX/ziF7F161bYbNZouigETpsAj0OMp7kNnjAzwf4+9QRtJCBJRRR4NM+p0FQmFZ9JC6cmSoH+BBx903L1BpzMKh8gl04kKpOSJOFw/xiWzq8y5T3ieQ5OG2/MGkijNjQRow04rAQk04VWIWomldR+wveJKalmdnPPln0vE9ctrse/vd6Df/5td9pt/+fW5fjYxvaMj8Fsga5QCSYX1bvhtPGqTTiXJmVVctn8Kty8zHiacbbDgslcu/CPlUEnNxBXJve80QcAeM+K+GekwePQtdEzwpFz4xj2BvDBq7Kb8e22izGHCWnWXQDp7sC/+MUvirWOklPvsStXM2bVTMbT3Mkn/HiaO7vHb6934/XTI4hGpTRVYCqPmsmZjN5sbi3TcjZOUa8BJ9OFQL3HkZS+GxidxogvaOp4sgqbYFo3NwvGMhmXswutTAp4pm7unGZzq9Qgx9Pc5gV9s2nf0+PmZXPxyv+6CSO+IHyBMHyBMLyBCP7+l2/jkZdO4UPrWjJeGJ286EWlU8TcqvTmPlHgsWx+FY7GmnAST77/7+XT8Iei+Ot3kSqZDTVKzWRuae6jA+PgOWDZfGsHk6wh5vSQD267gGs74q4XDZUOHOwbRSQq5TxX/IW35HGS2ZZXuBwCJEnOFOZbJlduZL+jz1ASRyoqXnp5fhjYCTq1u0tRZ8TsPujtdW4Ew1EMTvjTblO6uXNIH85kdBtwNEzLbQasgTKlaRsqHZj0x0cqHjboL5kNLrtorJvbwJrjDTj6igab6JRxAk5M3TUzmGSmwIk1k5P+wtVMEkBbnRtrW+fghs4GvHfFfHzwqmbct2kRLk4EsHd/v+7fSpKEkxcnsSSlkzuRlU3VGJ8OYWBUtjw7dWkSdz7xBp54rQfLF1Th3aRKZgWrmdTymjw/No0/nhnR/Puj58axuNFjeVGiLqG7+qYrGpMuahoqHYhKwIgv91T3b9+6hFq3HWtasxsn6dHIRs4GKPKIUe9x4LIvgEhUMj3NnVozeSmmgFZnaf7dHuvo7lXp6PaH5FnfYg4n6ZlMfAKOms+kusoo6piWZ6NMAvH5sUow2WpeMOm08Vkqk5lrJidNsgZSgvjUNHfEWDCuhqJMJij9LJg0M81N6PORa9pQ73HE1EPtz9+QN4CxqZBqJzeDeRn+rnsYX332ON7zz7/Dy+8M4QNrmvD9O9eTKpklNRm6uR/8xVv48GN/QL9KB/2oL4hzY9OW9pdk1CYEk+9ennzBwdw0cq2bPDc2jbcGJ3DTFQ1ZK5uKfdks7OimyCNGfaUdUQkYnQqal+a2qzfgHOiR7TbWZnnV0xbr6FarB5kORSx/NVkKRIEHz2lbA9kELm3DMDJOMVMwlOpderh/DFVOEQtN7IZ32UWDs7mNjVMEjKW57bFJD3poprljE3BySnM7WAd7Yje3fNIkZbJ4VNgF3LdpES5M+PGjA9rqZHesk7uzMb1eksECly/99Cj+7fUerGquxk//6lp8c1sX5lY5zV34LIBlGLRqJt+5MAlJgqqqfOx8zKzc4p3cgFxm4xB52AQOm69MTkXnaw/0opLizl4VZ023pEzOYhLtgcxvwEk+4e/vuQybwGWd8mTG5T0qTTjTwUhBbYHKGYcoqFoD+UNRpSs6EdbNnToKEEg0Lc8UTMa9S0ORKI6dG8dqgx2wRqmwCdnN5tZZs9suguOM+Uy6DFxk2bTGKUYisdtzsAbSqZnMtySFyA5ZnbTjOy+dVlX9gYTmm3nawWRnowfzqpxYUO3Ev3y4C0/vuDbr1CIRRxR4VDpF1ZrJSFRShl786EB/2jCNchijyOA4DtcsqsNtXU2oSqmXzjeY/O1bl2ATONzQqT59TI/4YIXZp0zSDhwjriQFDdeFZUItzR2NSjjQcxkrm6qz7uptnlMBgefQo5LmngqSMqmFXeQVdS6RQDii2uHMlEq1yUVKajxD81R9wob29uAkAuGoqfWSAOC0CwZnc0fAc/HJPmrwPAePXTQwASdiKHCzKz6TyetTlMk80typNZMeh2hqkE5kRlYnO/C3z72FH+3vxx0qnd3MY7JTJ80tCjxe+OwmZa4xkT81LvX53OdGpxGMROGyC7g0GcCLb19KmrB17Nw4OA6WnnyTyH/85dWqv1eCyRy8Jn2BMPadHsGGRXU5NfV5NErbZgP07Y2RqExOGexYzYRaMW73JS8m/GGsb6/N+vFsij1Qcpp7cHwa50anKdWngV3k1U3Lw1HVDmebkubWMS3XsdkBkqfgHDZx8k0iLpscTGYylWbHman+rNIpGmjACRu6yGLBIgseGfnM5q5QmYAz6Q/T575EMHXyEQ118uTFSdS4bMp3QQu3Q6RA0kRqKuyqweTpYTm4v/v6heC5uK0O4+i5cXQ0ePLOyJWafJTJ33UPIxiJYsuVuZnks6zNbBypSN/gGA2V8bSkUS+9TLhVTMvf6LkMAFiXQzAJyN2VPSM+RGOBTiAcwX27D2IyEMaOTR15rXem4hB5jW7uqGpQqGcNFFcmDdZMeoOmTr5JpMIu21Co2R4lonWcqXicYsYGHJ/BNLeDNT5pNODkUjNpF3nYBT5tAg4135SGCruAT9zYIddOptTgKZ3cjdqd3ERhkJXJ9DT3mSE5o7Whow5/cuVcvHJyCOfG5C76sakg+i9PK1OJypk6d+7BJLMEymaEYiJaTbezAQomY7CT/5A3oKS58x6nqNR4xT9YB1gw2ZZbXdDCOhcC4SguTvohSRL+9zPH8Wb/GP7qpg78DxM9DGcSdpFXVU784YiqMinqpLmDBpXJuljN5JA3gMP9o2iprUBdBoUmWyqUkYr6qe5AOGKoe9rjEA004ESyUibTJuCEoxD49KYno7gcQooyGTJ1+g2RHR/Z0CrXTr58Gq+fGsbjvz+Lz/34Tdzy8O8x6Q/rpriJwlBdYYMvGEn77p0ZkpXJjgYP/vzqFkQlKBcBx87Jk2/KoZM7E3aRxxyXLetgMhKV8OLbl3DF3Eq01Lpyem41x4nZAgWTMRJrJqcMdqxmguc5uOzJJ78DPaPobPQopqvZwkY7nR324Ydv9GHvgX7cuKQBn333FXmtdSZTYRNUv9yBUFQ1yNJqHgHi3dyZlEmbwKPGZcPZIR9OD/mwutlcVRKIz2GfylA3qZXOT6XSadNtwJEkCb5gOKuayVSvzlAkmlPzDcNtF5OVSX/YVMNyIjtcdhGfuLEDg+N+/MW//hH/9+cn8NR/D2B8Koh3LW3EHRvbSr3EWYcyBSfFHuhMzOC7sdKBTUsaML/aiR8d6EckKinNNytMGPVqBRoqHVnXTB7uH8OIL4gtecyBV8tGzhbokj6G2yGiwiZg2BuAPxTJ2xYo8XHZB+vc2DTOjU3jz69uzfnx2uvlK6afHjyHZw6fQ2utCw9/uCtnpWc20NHgwbNvnseEP5TU+RcIR9Km3wD61kBMWTPi51nvceDEoHzFb3aKG4jX9GZWJtWD5lQ8TtlqKByJqh6fXJ9pzOWABeTppuVSTiluhssuKN3ckagEXzCCSkpzl5Q7NrZhfDqEGpcNy+ZXYen8qpwvlon8mRPzmhyfDir1gwBwZtiLRQ0ecBwHUeDwoXUt+JcXuvHKyUs4dl5uvlk+A5RJQA4mjwykj+nU41fHBgHknuIG1Kd0zRZImUygvtKu1EyaVYTsccQ7ZFmKe3177tYXzB7ox/89AJHn8d07rlKMagl1VjTJV9snzk8k/d4f0mjA4bVNywPhqOHmEWYPBABrTDQrZzA3AENpbgM1k8yjTqujW7HhMXChxVT9tNncEWOBrRYuh6g0yLH1UM1kaXHaBHzuPVfg7hsW4drF9RRIlphqlSk43kAYFycCWNQQ97n90PoWcByw541+HDs3joX17hnzXWrwJE8gy8Tg+DSe/EMvFjd68rrwp5pJAgAbqRjAVMCY/YkRZCVF/mDtV4LJ3JpvAKB5jgtMhPzHD67C0jKxcSglzIT3eEIwKUmSZi0hUya1fCaNBGZAvHRC5DksL4ARsKJMZkpzawTNqbATiVaqm011MHKh5dCpmcxHmXQnKJPxudwz4wRIEGagTMFJCCbPxppvFtXHa1ibaipw05IGvPj2JfSOTM2I5htGth3d//ird+APRfHl9y3NK8vnpm5uApBP/iPeICYDYZPT3PLJ70DPKOZVOdE8pyLnx7OLPO6+YRH+5pal+NPVC0xZ40yHBXLHz8XTHuGohKikPhVGN80dyUaZlDe0pfOrsvYUNYJSM5lh4zKc5nbIikYmZdJQA47iM6lWM5lPmjteM8mahagBhyDiKPO5E2omz8RsgRKVSQD48NWtiMRcK2ZkMGmgbvJw/xh+eugcNi1pwE1X5F4vCch7MseRMjnrqfc4EI5KGPYGTE9zj0+F8M7FSaxrn5O3VcaXblmKe25cZMr6ZgPVLhtaaiuUcWEAlPSHU60BR2cCjtGUMRDf0Fa3FGaTZt3cmVI5Rru5mcKnpUz6svBf1RynmG8DjkNAKCIhGI4qnpjUgEMQcVgDTqI90GmmTKYEk39yZSMaY/vUTOjkZhhVJiVJwv/9+QkIPIcvv29p3s/LcRzcdjHrbu53LkwaTslbFQomE2iI1bhJknnj2dwOEcFwFH84OwJJyi/FTeTO8vnVOHXJq9QX6pmPC4K+NZBRZZLNFl7TUpjxcHFlUnsTktP5BtPcSs2kunF5fGa9EWWSlQqkN+DYDaxFC1fCFBzmiUkNOAQRR62bm9kCLaxPDiZtAo9PbOpAa61rZimTHnnvzRRM/vzIIP67dxQfuaYVnXO1x35mg9shZJXmPnVpEv/jX17FY6+eMeX5SwXtwgnUJ3S+madMyifOl98ZAkDBZKlY0VSFXx2/gLcuTGBt65yEedU6DTgapuVGLzRuWTkPk/5QwcoRWGClVzMZikiQpMy+mEC8AUdbmTQ+ZlRLmQyGo8rz5IKbTcEJRijNTRAqVFfIoshogjJ5ZsiHBdVOVe/ku65fiLuuX1i09RUDI8qkPxTB3//ybVQ6RXz6XUtMe+5EBxcjvHpyGFEJ2N87atoaSgEpkwnUJ5hKe8yqmYx9eV955xIqHSKumGfO1Q+RHczygtVNKmluXWsg9QYco/6jLruIj1+3sGCj4irs8uPqdXMzX8xsaiYzpbkNKZO8es1kMN+aydhzTwXCyjqpAYcg4qR2c0ejEs4O+7CoYfYYyBupmXz892dxbmwaf72lE7UmOhBkm+bed2YEgHxuyjQa18pQMJlAYjDpMkmZZCfe8+N+rG2bQ36QJYJ1dLNJD4GQtjKpNwHHaDNLMaiwxZRJ3WBS+zhT8Ri1BjJQM8nzHGwCh6DKbO58u7mBmDIZS8caAowRAAAgAElEQVTPFDsTgjADu8jDbReUNPeFCT+mQ5G0esmZTE2FDSLPaSqTlyb9+M5Lp9Be58LHNrab+tzZpLkjUQl/jAWTI74gLkz4TV1LMbHGWdEiJPoCmnWCSnycfPwlifxoqHRgbpVDacLRU+w4joPIc5rWQIVSGrOFNeDoTcCJB5PZNOCo10yy2kyjJSA2gU/v5s7z9WPNP1OBsJLmJmWSIJKpcdkVZfKMYgs0e4JJnudQ73FoBpPfeek0fMEIvnTLUtP3c7ddVGzUMvHW4AQm/GHMrZKFrKNZGq1bCWucFS1CYs2kkY5VI7iTgkmqlywlKxZU4+TFSQTCEfhjyqSWZY8ocKrWQEY7o4uBkQk4gVigaaRmkjWyaM3n9mVpEm4XeYRSayZNsAYCZGVyQgkmqZubIBKpcdkwNi3XTMZtgWZPmhuIjVTUCCZfPTmEltoK3Lws92k3WrgdIoKRaFq9uBqvnx4GAGy/th0AcCxlsEY5YY2zokWodIjKVYpZDTjMr9ImcFhdgJF6hHGWN1UjFJHQfdGbsZbQxvNpE3CiUQmhiGQoZVwMjEzAySbNzYIyrZpJbxbWQIA8nztRmZQkKdbNnb8p8FQwTBNwCEKDGpctXZmcRWluID6fO7UO8bIviDPDPqxvq83bpk+NxD0qE/tOj8AmcPjI1W2wCRyOnSNlckbAcRwaYnWTZp2gWAPOyqbqghhXE8ZZsUCeFnT8/LiuNRAgK5NqowABWCbNbWQCTjZpbqeNh8BziuVOKrmkuRNfQ1Y2YIoyGYhg0h8Cz5mXRSCImUJNhR2T/jDCkShOD3nhtPFYUJ37sIxypMHjQDAcVTIYjIOxrum1bYUpO2PnfK3ac0Y4EsX+nlF0tdSg2mXDFfMqKZicSbC6SbNOUKypYf1CSnGXGmbKe+zchBJkOTUUO1Hg07q5swnMioFN4CHynK7PpJLmNrBmjuPgcYiaNZPeQBgizxk+fofIJ6V6WGBp1KdTjVRl0uMQC6IuFJJoNIpPfepT6OjowOLFi/Htb39b9X5+vx+33XYblixZgtWrV+Pmm2/GqVOnirxaohypjnlNTvjDODPkQ3udG/wsa/7Usgf6775YMNlamGBScZzQ2ZcB4Oi5cXgDYWxcVAdALsO6NBnApTJtwrHGWdFC1JusTK5tnYN7bliI7SZ3jBHZM7/aiVq3HcfOjyvWQFrKpI3n0rq5WWrcKsokIBuX601OiCuwxi6OKp2i5hW1LxCGO4vgLbUBhwWWtrwacOLKpNcfLst6yd27d+PEiRM4efIk3njjDXzjG9/A8ePHVe9777334p133sGbb76JW2+9FXfffXeRV0uUI2yk4oVxP86PT6NjltVLAjrBZO8o3HahYDZ9zFYwkzL5+mm5i3tDhxxMMvu6xElt5YR1zooWgQWTZtVM2kUef/O+ZVhQM7tSDFaE4zgsX1CFtwYnMBX7omvVEgpCejd30GLKJCB3dOvV5mSrpnocom4DjhFbIIa9EMpk4gQcf7gs6yX37t2Le+65B4IgoLa2Ftu2bcOePXvS7ud0OnHLLbcowfuGDRvQ09NT5NUS5QibgvPmwBgkafbVSwLqXpOhSBRv9o9hTWvhbPqUKV0ZOrr/cGYEdpFXFNKVTcn2deWGdc6KFmHp/EpU2ATMi43CI2YWK5qq4Q9F8dbgJIAMDTjRVGXSWjWTgBxMToe0uwbjjUbGgsAqp03btDwYyeoiS/aZTFAmY8FkPrO5XQ7mMymPUyxHW6C+vj60tbUpP7e3t6Ovry/j3z388MO49dZbC7k0YoZQE5uCw+oDZ3UwmaBMnjgvlzgVql4SiGc19ZTJYDiKAz2juKp1jtJLceW8Sgh8+TbhlN9OXGDu2NiOD6xpVmpOiJnF8lgTzoHeywAyWANpKpPWafjImOYOZalMOkXNBhxfIKzMGzeCXVRvwMknGHcnXPV7/WFLjlLcuHEjuru7VW87dOhQ2u+MTL148MEH0d3djRdeeEHzPrt27cKuXbuUn71er4HVEjMRdv5i9YGL6mdhmtuTHkz+dyy4vqqAwaTihauTMXpzYAzToQg2xlLcgHwu6mz0lG0waR2JxSIIPEeB5AyGTcI5HbPL0AqyRF67AcdqyqShNLcBn0lAvqoOhqOKopmINxBWGmCMYBOS09xKzWQeaW6njQfHAePTIUyHIpZMc+/btw/Dw8Oq/1paWtDa2ore3l7l/r29vWhtbdV8vIceeghPP/00fvnLX8Llcmneb+fOnRgYGFD+eTyzL4AgZFjN5Gy1BQLUlcmDfaPgOKCrgDZ9bE/Sm8/9+im5XvLahGASkDNn58f9GNEZA2lVrHNWJIgi0FrrUsy5AW2V0SakN+BYsWbSZRcMzuY23oADpBuXS5KEqWBEUQaN4BD5pLrTUCT/YJLjOLjtIi7FThDl2IBz++2347HHHkMkEsHo6Cj27t2Lbdu2qd53165d2LNnD55//nnU1JBPLWGMGld8mltDpaMsvyf54naIcNmFpJrJg72jWNJYqcwvLwSsm9unsy/vOzOMCpuAVc3J32lmX1eO5uXWOSsSRBHgeQ7LYl9YQFa61BAFHqGoeprbUsqkTTDNZxLQns8dCEcRiUpZKYF2Ue7mZmncoAkNOIAcQF+M2WeUY83kHXfcgSuvvBKdnZ1Yt24ddu7ciZUrVwIAnn32WaVje2BgAJ/97GcxNjaGzZs3o6urC9dcc00pl06UCXMSsmuzaYxiKolTcM6PTeP8uL+g9ZJAvJtbS5n0hyI42DeGde1z0s4lK5tZE075pbrLbycmiDxZ0VSNP56Vaya1FDtRxxrIUjWTdhGhiISQxphCpWbSYJq7SmMKDgsus2vAkZ+TTb0JmRSMux0iekfk9F2lBdPcmRAEAY888ojqbVu3bsXWrVsBAM3NzYbqKQkilaoE5W22jVFMpMHjQM/IFAA5xQ0Utl4SSLYvU+Ng3yiC4WhSvSRj6fwqcJw8WKPcsI7EQhBFYkVTXJnU9JkU0scpWlOZlNeipU5mGwAz5TE1mGQ2F64saiaZAskUyaAJaW5AViaZaGzFBhyCKDVOm4CKWHNhxyysl2Q0VDpw2RdAJCoVpfkGiF9waymT+06zesn6tNtcdhEdDR4cLUNl0jpnRYIoEqwJB9BpwBE4hDSsgaxVMylvXH6N+pxcfCYBpE3BUeZgZ1EzyczJmSIZMsEaCEBS3aYVG3AIwgowr8nZ2HzDaKh0ICoBI74ADvaOotZtR3uddhObGTAvXp9GY+S+0yPwOESlPjKVlU3V6L88jfEp9UlkVsU6Z0WCKBKLGjxw2njYRV5zmovIl4cyyayNtEZ3xZVJY2uu1KiZZBtjNmnuNGUynL81EJCsjs7GxgKCMAJrMpmNtkAMZg/Uf3kKx89PYG1rTcHHr4oCD4fIqyqTU8Ew3hwYw9ULayFqZGiYfV25pbqtc1YkiCIh8BxWNdUo9YFq2AQO4aiUVLOWbWBWDJinmWaaW6mZNJjm1ggm4zWT2U3AAeJBuBkTcIBkZbIcG3AIohjUeeywCzya58ze6WvMHuiFty4hHJUK3nzDcDtE1W7uIwPjCEUkXLOwVvNvV8Qm4ZRbqpt2YmJW8g8fXIWxqaDm7WzUVjgqKWlZK6a5KzIqk9mtWasBh9VM5qdMmlczyaA0N0Go87l3X4GLEwFNBWw2wILJXx2/AAC4qrVYwaSgqkyeuiQPEtCbC768TO2BaCcmZiUL690AtGuJWMATjkhgop5iWi5YqZtbXovWFJxAOAKek7vTjaDVgOPLo5s7VZm0mdDNzSBlkiDUWVOkwMnKsGDyzJAPIs+l+ToWCrddVL3AZybyHTod9pVOGxbWu8vOHmj2XrIQhA6iokymT3AxarNTDIwokw5RMFwnxNLcWg042ZiWszQ3CyLNasBJUiYpmCQIQgMWTAKy4ldhL44Q4HaIqrO5zwx7YRd5LKjRLz1Y0VSNs8O+tH3YyljnrEgQFkJMUCYZcWXSOl8bIzWT2QS/Wg04U8H8ayaDsdcy3zKBJGXSQQ04BEGoU+eOB5PFqpcE5H15SiWYPD3kxaJ6t1JGpQXr9D5RRqlu65wVCcJCMPUsZHFl0smCSQ0bikA4klXw5hAF2AU+bZyiN1YzmdUEnNhraHbNJFNjRZ7TnGBEEARhF3llGlCh/SUT8cQacKIJU9T8oQgGRqcNWTWVYxMO7cQEoYLIqymTckBlKWXSxoJJ/TR3NlQ6xbSayUux8YVZNeBodHPnG0wyddTjFAtu80EQRHnDUt3FDCbZPpmYMeodmYIkGbNqWjpfViZZw045QAVHBKECUyYTg8m4Mmm9BpwpzQacaNZpZY9TxGQsRXNh3I///bNj+M2Ji2ieU6F41xkhcZyi/F+zurnlbYs6uQmCyMTKphrYRR7zq4tnkaQYlwfCSmB5ekgODDsaMyuTc1w2iDyHYa+244jVoN2YIFQQVdLcVq6Z1J6AE8kqAATkIG1sKoQf/LEXf/+LtzEZCOP9q+bjgT9dnlUgmF4zaY61ElMmybCcIIhMPHT7KkSixZ1xr4xUTNiXz8SCSSPKJMdxqPPYMeILFGaBBYCCSYJQQS3NHQxHwXH5dyObScYJOKEoHJXZp7mPn5/A3/z0GOZXO/HPH+7ClqVzs16bYg0Ukddmns+kvG1VkjJJEEQGOI5TxIFioTaf+3TMFsjoeMs6twMjZaRMllxiee6557Bu3To4HA58+tOfLvVyCAJAQgNOJKEBJxKFXdAewVgKWDOKZjd3DmnuuVVOcBywfWMbnt+5KadAEkiwBgqnprnNmc1NtkAEQViRxDQ348yQF42VDsMZlTqPHSNeUiYN09nZiccffxw//vGP4fWWT7EpMbMRmDIZTW7AsdL0GyCu0mk34GS/5v+zdQXu39Kpa6xrBFYOEGA+k7GgMl/TcpeS5i759kUQBJGGS0lzy8GkJEk4M+TD8qYqw49R73HAF4xgOhgpmj9mPpT8zLhkyRKsXr0aokgnBsI6xBtwkq2B7Fl2Rhcah8iD49TT3JIk5dTNXe2y5R1IAonKpLmzuVkQmW0tKEEQRDHwKGlueV8emgxgMhDGoiz21Tq3HQDKpm6y5MFkNuzatQvNzc3KP1IyiUIhJszmZuSSMi40PM+husKGUZU546GIBEkqnS9m2mxuk7q5Gyud+Mc/W4W/vG5hfgskCIIoAK6UNPdpA2MUU6nzyJZG5VI3WXA5cOPGjeju7la97dChQ2hpaTH8WDt37sTOnTuVn5ubm/NeH0GooTYBJ2jBYBIAat12XPalbzjMF7NUa1asgcJx03KB5zJOfzDCh9Yb3zcIgiCKiSelm/vMcKyT22DzDSDXTALlo0wWPJjct29foZ+CIExHbQKOFZVJQE6HnIld+SbCrIyyTXObhWINlDCb20qd8ARBEIWA1bIryuSlmDJpwBaIUR8LJsvFa9J6Z0aCsABa1kBWDCZr3XaMTgWTRncBicFkidPcSs2kZCmPToIgiELgSWnAOTPshV3k0TTHuHE6myteLmnuku/sL7zwApqbm7Fr1y48/vjjaG5uxrPPPlvqZRGzHFGlAScQjihqm5WodTsQlYCx6VDS7wMxu6CS1UyK6TWTVnz9CIIgzIQ5TjBl8syQDwvr3FmV+Chp7jKxByp5C/WWLVswMDBQ6mUQRBJKvV9aA461urmBeDpkxBtAbawDECh9mpultIMJNZP5Nt8QBEFYHaZMTgUi8Ici6B+dwnuXz8vqMRRlUqUe3orQzk4QKijd3JHkmkkrKmu1ioVE8qZT8jQ3swZKqpm03utHEARhJg6RB88B3kAYvSNTkKTsOrkBoMIuwG0XMFwmyiTt7AShgi2lm1uSJEvXTAJI6+hW0twl7uYOJvhMWjEYJwiCMBOO4+B2iJgKRuIzubPo5GbUecpnpCLt7AShAqttYd3crO7PisGQVjqEKZNsfnexYc02oYik/JeUSYIgZgNuuwhvIIzTsWAyl0EQdR572VgD0c5OECrEG3DkQChY4pSxHooy6dVIc5eoAYfnOdgETllHMByFnayBCIKYBbgdAqaCYcW2LSdl0i0rk5IkZb5zibHemZEgLICS5o414LCAyIrKZL2GuW3ctLx0TUM2gVdqJoNUM0kQxCzB7RDhC0RwetiHhkoHKp3Zj3+t99gRjkqYmA4XYIXmQjs7QaiQ2oATLHFntB5ztBpwQqVXU+0iTzWTBEHMOlia+8wlLzpyUCWBuD3QcBmkumlnJwgVykmZtAk8qpyidpq7xMqkMgGHrIEIgpgluB0CxqdDmAyEsSiHekmgvIzLaWcnCBVYzWQoTZm05lemzuNI7+YOl9a0HJCbcCjNTRDEbMPtiNt4L6rPT5ksB+Ny2tkJQoXUcYpBCyuTgNyEk14zWfoAmKW5JUmSxymK1IBDEMTMh83nBoCOxtyUyXqPrEwOl4FxuTXPjARRYtj0FmYNZIVmFj3q3HaMToWS5nPHayZLt2Z7LM3N7IFoNjdBELMBjyO+73bU55jmJmWSIMobUSgvZbLOY0ckKmE8YT53PAAuvTLJUt2U5iYIYjbAlEm7yKNpTkVOj0E1kwRR5thSurmtkDLWQ22kYql9JgFZ4Q1FokowbrPo60cQBGEmbD73wjq3MgQjW+a4bOC4dNs3K0I7O0GoEJ+Ak9zNbd1gUr6CvZwUTJY+NZ+qTFKamyCI2YArlubOxaycIQo85rjsGCZlkiDKk3iaO7lm0qrBUJ07vbbGCj6Tsmm5ZOlxlARBEGbDlMl8gklA3tuHqWaSIMoT1oDDfCaDFkgZ66EUaquluUsYwDkUZVJ+HW00TpEgiFlAc6xOcnVzTV6PU+exl0XNpJj5LgQx+0i1BlJMywVrdnMr87lT0twCzykqaylgpuVM2aUGHIIgZgNXtdXitzs35Tz9hlHncWB8+jKCYWtPELPuygiihMSVyRTTcqsqk6o1k9GS13iyzW8qSMEkQRCzi8WNHnBcftmY+phQMDplbXWSdnaCUIHjOAg8p6Rn48qkNb8yc9w2AClp7lDpg0kWPPoCYQDWff0IgiCsSB0zLrd43STt7AShgchzSgOO1ZVJhyig0iEmN+CEIyU3WWfKpC8QSfqZIAiCyEzcuJyUSYIoS2wCrzTgWL2bG5A3nbQ0d4mDX3uKMklpboIgCOMoxuUW95qknZ0gNBBjhttAojJpzQYcgM3ntmbNpC/Igknq5o5Go/jUpz6Fjo4OLF68GN/+9rcz/s0TTzwBjuPwzDPPFGGFBEFYhYbK8lAmqZubIDQQeV6lm9u611+1bgeODIwjGpXA8xwC4QhqKuwlXRN7vbysZpLS3Ni9ezdOnDiBkydPYnx8HGvWrMHmzZuxfPly1fv39vbisccew4YNG4q8UoIgSg1TJq1uXE47O0FoYBM4ZQKO1WsmAdncNhyVMOGX53NbqQFnKmD9MoFisXfvXtxzzz0QBAG1tbXYtm0b9uzZo3rfaDSKu+66C9/61rfgcDiKvFKCIEpNvGaS0twEUZYICQ045VAzWZtiXG6JmkkxWZmkmkmgr68PbW1tys/t7e3o6+tTve+uXbtw3XXX4aqrrirW8giCsBAehwi7yCeVMFkRSnMThAY2IZ7mZuMAS6306VGXYFze0WCNbm5WI6kEkxZ+/cxi48aN6O7uVr3t0KFDab+TJEn1vsePH8dTTz2F3/3ud4aed9euXdi1a5fys9frNfR3BEFYF47jUO+2W16ZpGCSIDQQeS7JtNwu8Hkb0BaSRAsJSZIs0YDjUEzLZ08Dzr59+3Rvb21tRW9vLzZu3AhArolsbW1Nu9+rr76K3t5edHZ2AgAuXLiAe++9F4ODg9ixY0fa/Xfu3ImdO3cqPzc3N+dzGARBWIQ6j4NqJgmiXBGTrIFKH5hlojZhCk4oIkGSSq+k2pQGHLlMoNTrsQK33347HnvsMUQiEYyOjmLv3r3Ytm1b2v127NiBwcFB9PT0oKenBxs2bMD3vvc91UCSIIiZS53HjhFfQDOLYQVoZycIDWwCl9TNbfVOZJbmHvEGlBrPUqe546blVDPJuOOOO3DllVeis7MT69atw86dO7Fy5UoAwLPPPou77767xCskCMJK1Lkd8IeiylhaK0JpboLQQOTjPpPloUzGG3ACFuk+p2AyHUEQ8Mgjj6jetnXrVmzdulX1tpdffrmAqyIIwqrUJ5QwuR3WDNtoZycIDRLT3MEyUCZrExpwlGDSImnuuGm5tV9DgiAIq8Hq4YctPAWHdnaC0MAmJCqTEcsHk06bAI9DlIPJkNXS3FQzSRAEkQvKSEULN+HQzk4QGiROwAmGoyUPzIzARipaRZlMnYBDyiRBEER2lINxOe3sBKGBTYhbA5VDAw4QCya9AcvVTLIJQrPBGoggCMJM6j0xZdLCxuXWPzsSRIkQeC5msSPFlEnrf13q3HaMTgXht0iaO1WJnA2m5QRBEGai1EySMkkQ5YcYC4QiUaksaiYBedMJRSRl0yl1AJw6ftLK4ygJgiCsiOLUQTWTBFF+2Hg5JRuOlo8yyYzLB8f8AKyQ5k5Oa1PNJEEQRHY4RAGVThEj1M1NEOUHUyYDoSiiEmAvgwYcZlx+bmwaQOnT3HYh/vwCz0HgqWaSIAgiW+o9DlImCaIcYc0i3phHYnkok3IwOTjOgslSK5Px56fmG4IgiNyoc9stPZ/b+mdHgigRIi9/PaZitjblUDNZ62HBZCzNXfIGnHgASfWSBEEQuVHnseOyL4Bo1JrzuWl3JwgNxFgg5AuWj+F2faxm8rxlaiZ51f8nCIIgjFPncSAqAWPToVIvRRXa3QlCA2UUYBkqk1bp5k5suKHmG4IgiNyod1vbuJx2d4LQQIw1i7BgstQpYyOwBhxGqddsp2CSIAgib+pixuVWrZuk3Z0gNGDd3L4yasBx2gS47PEAstRr5nlOCcrLQdklCIKwIspIRYvaA9HuThAaxJXJ8qmZBOId3UDpayaBeBBJyiRBEERu1MXq4a1qD0S7O0FooDTglFHNJBBPhwClT3MD8SDSTtZABEEQOVFv8ZGK5XF2JIgSYONZmltWJsvF2iaxbtIKaiopkwRBEPlBNZMEUaYwZZL5TFohZWyEWqsFk0yZtMBaCIIgypGaChvsAo+LE/5SL0UV2t0JQoN4Aw5TJkufMjYCUyYFnlOOoZSQMkkQBJEfPM+hpbYCPcO+Ui9FFdrdCUIDW5o1UHl8XZgyaZX1MmWSgkmCIIjcWVjvRv/oFMKRaKmXkgbt7gShAVP1poLl2YBjlWDSJjJrIGrAIQiCyJW2OjdCEUmZcGYlrHG2IQgLwuZKe8tMmaxTlElrpOWVmklSJgmCIHKmvd4NADg7Yr1UN+3uBKGByDNlMlYzWSbBpJLmtkjDkI3S3ARBEHmzsE4OJq1YN0m7O0FoIKYpk9ZQ+jJhuZpJ1oBjkfUQBEGUI+31LgDAWQomCaJ8sCnWQOWlTLKxW1YJfinNTRAEkT8LqitgF3n0UJqbIMoHgZmWl1nNpMsuosImWCb4ZeuwynoIgiDKEZ7n0FbrsmSaWyz1AgjCqijWQMHyCiYB4BObFmFBTUWplwEg0WeSurkJgiDyob3ejRffvoRQJGqpOnQKJglCA2YNFJXkn62SNjbCp9+1pNRLUKAGHIIgCHNor3MhEpVwbnRa6e62ArS7E4QGYoqSRmna3KAJOARBEOZgVXugku/uDz/8MFasWIGVK1di1apV2L17d6mXRBAAABuf/PWgYDI3WONNOZUJEARBWBGr2gOVPM29fPlyvPbaa6iurkZ/fz/WrFmDjRs3oqOjo9RLI2Y5icqkyHMQeKr5ywVSJgmCIMyBKZNWCyZLvrtv2bIF1dXVAICWlhbMmzcP/f39JV4VQSQ3jJAqmTvsdaRgkiAIIj/mVTnhEHmcHZkq9VKSsNTu/tvf/hajo6NYv359qZdCEMoEHIBStPlgF+TGJermJgiCyA+e59Be5559yuTGjRtRX1+v+i9RgTx69Cg+/vGPY+/evXC71TuUdu3ahebmZuWf1+st9PKJWYxIyqQp2ET5daTXkCAIIn/a610YGJ1CMBwt9VIUCl4zuW/fvoz3OXHiBN7//vfj+9//Pq6//nrN++3cuRM7d+5Ufm5ubjZljQShRmJatpxsgawGTcAhCIIwj/Z6N6ISMDA6hUUNnlIvB4AF0txvvfUWbrnlFnzve9/DzTffXOrlEIRCYsMNqWq546AGHIIgCNNoZx3dFrIHKvnufv/992N8fByf//zn0dXVha6uLvz6178u9bIIIskaiGomc6eh0gEAqI/9lyAIgsgdFkyeHbZOE07JrYGef/75Ui+BIFShmklzePeyefivT16Plc3VpV4KQRBE2bMwgz3QL44Owi7weNeyuUVbU8mDSYKwKonBJCmTucPzHAWSBEEQJjG3yoEKm6Ca5vYGwtj5o8OorrBhy9JGcFxxXDToDEkQGiSmue3UgEMQBEFYAI7j0FbnwlkVZfI3xy/AH4ri4kQAg+P+oq2JgkmC0IDnObAeHFImCYIgCKuwsN6N82PTCIQjSb9/5vB55f8P9Y0VbT10hiQIHURma0PBJEEQBGERmD1Q/+Vp5XdDkwH8vnsIixtlu6BDfaNFWw+dIQlCB1tMmiRlkjCLaDSKT33qU+jo6MDixYvx7W9/W/O+gUAAn/zkJ9HZ2YmVK1fiox/9aBFXShCEVWmvcwFIbsL5+ZHziErAp/5kMaqcIg4WMZikBhyC0EFWJiMUTBKmsXv3bpw4cQInT57E+Pg41qxZg82bN2P58uVp9/3CF74AjuNw8uRJcByHCxculGDFBEFYDTWvyWcOn4fbLuDdy+bhJwfP4Q9nRhAIR4oydIPOkAShA5snTRNwCLPYu3cv7rnnHgiCgNraWmzbtg179uxJu5/P58Pjjz+Or3/960pH5rx584q9XIIgLAizB2JNOLy/TjQAAAxKSURBVGeHfXizfwzvWT4PFXYBa1trEAxH8dbgZFHWQ8EkQejApuBQzSRhFn19fWhra1N+bm9vR19fX9r9Tp8+jdraWjz44INYt24dbrjhBrzwwgvFXCpBEBalodIBtz1uD/Szw+cAALeuaQIArGmdA6B4dZN0hiQIHUSe5koT2bFx40bU19er/uvv70+7vyRJqo8TCoXQ29uLZcuW4cCBA3j44Yexbds2XLx4UfX+u3btQnNzs/LP6/WaelwEQVgH2R7IjZ7hKUiShJ8dPo96jx3XddQBALqaawAAB4vU0U1nSILQIZ7mpq8KYYx9+/ZheHhY9V9LSwtaW1vR29ur3L+3txetra1pj9PW1gae5/GRj3wEALBmzRosXLgQR48eVX3enTt3YmBgQPnn8XgKc4AEQViChfVunB+fxv6eUZwd9uH9qxYoDiTVLhsWN3pImSQIK0DWQITZ3H777XjssccQiUQwOjqKvXv3Ytu2bWn3q6+vx5YtW/DrX/8aAHD27FmcPXsWS5cuLfaSCYKwIO31LkgS8K0XuwEAt3YtSLp9TUsNBkancWmy8ObldIYkCB1EsgYiTOaOO+7AlVdeic7OTqxbtw47d+7EypUrAQDPPvss7r77buW+jz76KL7xjW9g5cqVuO222/Dd734XTU1NpVo6QRAWgnV0/657GG11LnS11CTdHq+bLHyqm6yBCEIHm6JMUjc3YQ6CIOCRRx5RvW3r1q3YunWr8vOiRYvw0ksvFWtpBEGUEe2xjm4AuLWrKW0O95pWObg81Cd3eRcSklsIQgeRaiYJgiAIC8KUSSA9xQ0AS+ZWwm0XilI3ScokQehg46lmkiAIgrAe9R476j0ONM2pQEdDesOdwHNY3VKDQ31jCEeiSg9AIaBgkiB0IGWSIAiCsCIcx+E/792ASqd2KLemtQavnx7B2xcmsaKpumBroTMkQehA3dwEQRCEVVnc6MHcKqfm7WtZE05/YZtw6AxJEDrYeBqnSBAEQZQnrMO70HWTFEwShA40TpEgCIIoV+o8DrTVuXC4wPZAdIYkCB2YNRDVTBIEQRDlyJqWGpwZ9mHUFyzYc9AZkiB0oAYcgiAIopxZ2ybXTR4uYN0knSEJQgeRZ8ok1UwSBEEQ5ceaFjYJp3B1kxRMEoQONoFqJgmCIIjy5cr5lXDaeBwsYN0knSEJQodFDW7Uue2ocdlKvRSCIAiCyBqbwGNVUw3eHBhDNCoV5DnItJwgdLj7+kX4+HULlUYcgiAIgig3/v7PVqLO7QDPc5nvnAMUTBKEDjzPgUdhvnwEQRAEUQwWqYxbNBOSWwiCIAiCIIicoWCSIAiCIAiCyBkKJgmCIAiCIIicoWCSIAiCIAiCyBkKJgmCIAiCIIicoWCSIAiCIAiCyBkKJgmCIAiCIIicoWCSIAiCIAiCyBkKJgmCIAiCIIicoWCSIAiCIAiCyBkKJgmCIAiCIIicoWCSIAiCIAiCyBkKJgmCIAiCIIicoWCSIAiCIAiCyBlOkiSp1IvIFYfDgYaGBsP393q98Hg8BVxR8aBjsR4z5TiA2XMsQ0NDCAQCRV5RcaD9kY7FasyUY5kpxwGYtz+WdTCZLc3NzRgYGCj1MkyBjsV6zJTjAOhYZiMz6XWiY7EmM+VYZspxAOYdC6W5CYIgCIIgiJyhYJIgCIIgCILIGeGrX/3qV0u9iGKycePGUi/BNOhYrMdMOQ6AjmU2MpNeJzoWazJTjmWmHAdgzrHMqppJgiAIgiAIwlwozU0QBEEQBEHkDAWTBEEQBEEQRM7MimCyu7sb1157LZYsWYL169fj+PHjpV6SYe6//360t7eD4zgcPnxY+X25HZPf78dtt92GJUuWYPXq1bj55ptx6tQpAMClS5fw3ve+F52dnVixYgVeffXVEq82M+9+97uxatUqdHV14YYbbsChQ4cAlN/7ksgTTzwBjuPwzDPPACjP96W9vR1XXHEFurq60NXVhb179wIo7/elGJTr6zNT9kdgZu2RtD9ak4Luj9IsYPPmzdITTzwhSZIk/fjHP5bWrVtX2gVlwSuvvCL19/dLbW1t0qFDh5Tfl9sxTU9PS88995wUjUYlSZKkb33rW9KmTZskSZKkj3/849IDDzwgSZIkvfHGG1JTU5MUDAZLtFJjjI6OKv//9NNPS6tWrZIkqfzeF0ZPT4+0ceNGacOGDdJPf/pTSZLK831J/Z4wyvV9KRbl+vrMlP1RkmbWHkn7ozUp5P4444PJixcvSpWVlVIoFJIkSZKi0ag0d+5cqbu7u8Qry47ED8FMOKb9+/dLbW1tkiRJktvtlgYHB5Xb1q9fLz3//PMlWln2PPHEE9Lq1avL9n2JRCLSli1bpAMHDkibNm1SNstyfF/UNstyfV+KxUx4fWba/ihJM2ePpP3ROhRyf5zxae7+/n7Mnz8foigCADiOQ2trK/r6+kq8styZCcf08MMP49Zbb8XIyAhCoRDmzZun3Nbe3l4Wx/Kxj30MLS0t+MpXvoInn3yybN+XXbt24brrrsNVV12l/K7c35eVK1firrvuwtDQUNm+L8Vipr0+M+V4yn2PpP3RmhRqf5zxwaQa0gx0QyqnY3rwwQfR3d2Nv/u7vwMgf3gTKZdj+Y//+A/09/fjb//2b/H5z39e9T5WP5bjx4/jqaeewpe//OW028rxfXn11Vdx5MgRHDx4EPX19di+fbvq/crhWErJTHt9yu14ZsIeSfuj9Sjk/jjjg8mWlhYMDg4iHA4DkF+k/v5+tLa2lnhluVPOx/TQQw/h6aefxi9/+Uu4XC7U1dVBEARcuHBBuU9vb29ZHAtj+/bteOmll9Dc3Fx278urr76K3t5edHZ2or29HX/4wx9w77334kc/+lFZvi9sfTabDZ/+9Kfxu9/9rqy/L8Vgpr0+5X48M22PpP3ROhRyf5zxwWRjYyPWrl2L3bt3AwB+8pOfoLm5GYsXLy7xynKnXI9p165d2LNnD55//nnU1NQov7/99tvx6KOPAgD279+Pc+fOYdOmTaVaZkbGxsZw/vx55ednnnkGdXV1Zfm+7NixA4ODg+jp6UFPTw82bNiA733ve9ixY0fZvS8+nw9jY2PKz3v27MGaNWvK8n0pJjPt9Snn45kJeyTtj9ak4PtjTlWcZcbbb78tbdiwQers7JSuuuoq6ciRI6VekmHuvfdeqampSRIEQWpsbJQ6OjokSSq/Y+rv75cASIsWLZJWr14trV69Wrr66qslSZKkCxcuSDfffLO0ePFiadmyZdKLL75Y4tXq09PTI61fv15asWKFtGrVKmnLli1KUXO5vS+pJBaYl9v7cvr0aamrq0tauXKltGLFCmnr1q3S2bNnJUkq//el0JTr6zNT9kdJmjl7JO2P1qTQ+yONUyQIgiAIgiByZsanuQmCIAiCIIjCQcEkQRAEQRAEkTMUTBIEQRAEQRA5Q8EkQRAEQRAEkTMUTBIEQRAEQRA5Q8EkURZ0dXWhq6sLy5YtgyAIys/btm3DgQMHsG3btlIvkSAIoiTQ/kiUGrIGIsqKnp4edHV1JZmvEgRBELQ/EqWDlEmi7Hn55ZfR1dUFQN5Ma2pq8JWvfAVr165FZ2cnXnvtNXzmM59BV1cXVqxYgWPHjil/++STT+Kaa67B2rVrceONN+LNN98s1WEQBEGYDu2PRDGgYJKYcYyPj+Oqq67CwYMH8YUvfAHvec97sHXrVhw+fBjbt2/H1772NQDAa6+9hj179uDVV1/FwYMH8fWvfx1/8Rd/UeLVEwRBFA7aH4lCIJZ6AQRhNk6nE7fddhsAYN26dfB4PNi8eTMA4Oqrr8YPfvADAMDPfvYzvPnmm7jmmmuUv718+TKmp6dRUVFR/IUTBEEUGNofiUJAwSQx43A4HMr/C4IAp9OZ9HM4HAYASJKE7du348EHHyz6GgmCIEoB7Y9EIaA0NzFrufXWW/GDH/wAfX19AIBoNIoDBw6UeFUEQRClh/ZHIhtImSRmLddffz3+4R/+AR/4wAcQDocRDAbxvve9D+vWrSv10giCIEoK7Y9ENpA1EEEQBEEQBJEzlOYmCIIgCIIgcoaCSYIgCIIgCCJnKJgkCIIgCIIgcoaCSYIgCIIgCCJnKJgkCIIgCIIgcoaCSYIgCIIgCCJnKJgkCIIgCIIgcoaCSYIgCIL4/+3WsQAAAADAIH/rOewuioBNJgEA2ALvYSSAl+OM2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f05a4e5aa90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "T = 50\n",
    "sigma = 1\n",
    "rho = 0.99 # for AR1 only\n",
    "\n",
    "lXIndependent = independentSampler(sigma,T)\n",
    "lXDependent = AR1(rho,sigma,T)\n",
    "\n",
    "plt.figure(num=None, figsize=(15, 6), dpi=80, facecolor='w', edgecolor='k')\n",
    "# independent\n",
    "ax1 = plt.subplot(131)\n",
    "ax1.plot(lXIndependent)\n",
    "ax1.set_xlabel('Time')\n",
    "ax1.set_ylabel('Values')\n",
    "ax1.set_title('independent')\n",
    "\n",
    "# dependent (AR1)\n",
    "ax2 = plt.subplot(132)\n",
    "ax2.plot(lXDependent)\n",
    "ax2.set_xlabel('Time')\n",
    "ax2.set_ylabel('Values')\n",
    "ax2.set_title('dependent (AR1)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The marginal distributions (the target $\\mathcal{N}(0,\\sigma)$) should look similar, although the dependent approximation of the target is poorer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApEAAAGrCAYAAACYIsydAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAIABJREFUeJzt3Xl0VGWexvGnQtqgJpJhXyoLu+yYgIZ9ERpFWmh0VCCiGIRhEenYgiM6YMvQtCM1PcARAyNRR0W7IWK3UfsYBIE+4gElgLGFQJPNDuLIGodAKnnnD8ca0kSoF1K5t8z3c07O4S5171PE+vF4a/MYY4wAAAAACxFOBwAAAED4oUQCAADAGiUSAAAA1iiRAAAAsEaJBAAAgDVKJAAAAKxRIgEAAGCNEgkrt956q55++unLvv3AgQO1aNGi2gsUIuGSE0Bo+f1+eTwebdmyxekoF+V0zldffVVjxoyp8/NWVVWpd+/erv/9/FhRImHl3Xff1ZNPPul0jLCzZcsWeTwe+f1+p6MAQFAWLVqkgQMHXnK/8vJyPfroo1q8ePEF2zZv3iyPx6N77rmnxuM3aNBA0dHRiomJUUJCgn7xi1/o7NmzgX2+/PJLjR07VgkJCfJ4PPrP//zPaseIiIjQU089pTlz5lzGPcSVokQCAIDLtm7dOnm9XvXu3fuCbc8995yaNGmirKwsffXVVxds79evn8rKynT69Gn98Y9/1GuvvaalS5cGtkdEROinP/2pXnvtNXm93hrPP2bMGP33f/+33n///dq7UwgKJRJWhg4dqieeeEKS5PF4tGLFCg0YMEDR0dHq0aOHtm/fHtjX7/dr3rx5atmypZo1a6Z//ud/vuB4X375pSZOnKg2bdqoefPmmjBhgr7++utq55s9e7bGjx+vmJgYdezYUS+//HK1Y3z88ccaOnSomjRpooSEBD355JPVrvjVVc6HH35YEydOVKNGjRQXF6dVq1ZJkoqKinTrrbdKkmJjYxUdHa0lS5ZY/b0DqBtHjx7V+PHjFRsbq3bt2umNN964YJ9gZo7P51Pfvn0VHR2tG2+8Ubt27ap2jJdfflm9evVSo0aN1K1bN73++uuBbd8/c7F+/Xp16tRJMTExGjlypL788staz/lDs/HVV1/VkiVL9NFHHyk6OlrR0dHatm1bjX9nGzZs0KhRoy5YX1paqo0bN2rFihVq3Lix1qxZ80N/7ZKknj17atCgQdq5c2dgXatWrTRr1iwNGDBADRo0qPF2DRo00IgRI5SVlXXR4yMEDGBhyJAhZsGCBcYYYySZnj17mvz8fFNRUWEefvhhEx8fH9h38eLFJiEhweTl5Zny8nLzxBNPmMjISLNw4UJjjDHl5eWmc+fO5pFHHjFlZWXm9OnTJjU11YwYMaLa+Ro2bGj+8Ic/mIqKCpOdnW1+8pOfmO3btxtjjPniiy9MdHS0WbdunamoqDAFBQWmZ8+eZvHixYFj1FXO6667zmzatMlUVlaa9evXm4iICJOfn2+MMWbz5s1GkqmoqKjdXwiAWjVy5EgzatQo880335hvvvnG3HbbbUaS2bx5szEm+JnTvn37wExZuHChadq0qTlx4oQxxpjMzEwTFxdndu7caSorK822bdtMTEyM2bZtmzHm/+fFxIkTzYkTJ8yJEydM//79zeTJk2s958Vm48KFC82AAQMu+XfWokUL8/rrr1+wftGiRaZ58+bm3LlzZv78+SYuLs74/f4aj19VVWU++eQT07RpU/PLX/6yxvMkJCSYNWvW1Ljt3/7t30xycvIls6J2USJh5e9L5EsvvRTY9tlnnxlJ5siRI8YYYzp06GCWL18e2O73+02zZs0C5WzDhg2mdevWpqqqKrBPSUmJkWSKi4sD5xs/fny1DHfddZd54IEHjDHGPPTQQ2bChAnVtr/yyiumffv2geW6yjllypRqOZo2bRoYrJRIwP2+f1zv3bs3sG7v3r3VylmwM+f8mVJZWWlatWplXn75ZWOMMT169DAZGRnVjjF16lSTlpZmjPn/eVFYWBjYvnLlSnP99dfXes6LzcZgS+RPfvIT8+6771Zb5/f7TZs2bcy8efOMMcYcPHjQeDwes3HjxsA+CxcuNA0aNDCNGjUyUVFRRpK55557TFlZWY3nuViJXL16tWnbtu0ls6J2RdbdNU/8GLVu3Trw52uvvVaSdPr0abVo0UIlJSVq165dYHuDBg0UHx8fWM7Pz9dXX32lf/iHf6h2zKioKBUVFQVe/3L+MSSpbdu2+vTTTwPH2Lx5s2JjYwPbq6qqVFVVVec5zz/H9+c5ffq0AISHkpISSdVnTtu2bavtE+zMOf8YERERio+PV3FxceAY6enpmjdvXmAfv9+vwYMHVzvG38+t7+dJbea82GwMVuPGjXXy5Mlq69566y397W9/09SpUyVJ7du317Bhw/Tcc89p7Nixgf1SUlK0fft2+f1+rV27Vr/61a90/PjxQJZgnTx5Uo0bN7a6Da4cJRIh4/V6dfjw4cByZWVlYIhKUsuWLZWQkKBDhw5d9DjnH+P75e+LW8uWLTVx4kStXbvW8ZwXExHBy48Bt/t+rhw+fFjdu3cP/Pl8wc6c829XVVVV7X84W7ZsqaeeekqTJ092POfFBDu3kpOTlZeXV23dc889J0kaNGhQYF1ZWZn+53/+RwcPHlSHDh2q7R8ZGalp06bpww8/1EMPPaQ333zTKuu+ffvUp08fq9vgyvEvG0Lmvvvu07Jly/TFF1/o7Nmz+tWvfqVjx44Fto8fP14VFRV68sknA/8Xe/To0QteIP7OO+8oOztblZWVeu+997Rx40ZNmTJFkjRz5kytX79ev//973Xu3DlVVlbq4MGDeu+99+o858W0bNlSkrR///6gbwOgbrVp00Y333yz5s2bp+PHj+v48eNasGBBtX2CnTn/8R//ob/85S86d+6c/vVf/1Xnzp3T7bffLkmaO3eunn76ae3cuVNVVVU6e/asdu7cqU8++aTOc15My5YtVVRUpPLy8ovuN378eP3pT38KLB84cEAffPCBXnzxReXm5gZ+Dhw4IK/XG3jTYU2eeuopvf3229Xe/FheXq7y8nIZY+T3+1VeXq6KiorA9srKSuXk5OjnP/950PcNtYMSiZCZP3++xo8fryFDhsjr9ercuXO66aabAttjYmL00UcfqaioSD169NB1112n/v37a+vWrdWO88ADD+iFF15QbGysZs2apeeffz7wf7d9+/bV+++/rzVr1qhNmzZq0qSJ7rzzThUWFtZ5zovp1KmT5syZo2HDhik2NrbaR1gAcI9XXnlFV111lRITE5WUlKS777672vZgZ86MGTN07733qnHjxvrDH/6gd955J/DU8sMPP6xFixbpn/7pn9S4cWO1adNGjz76qL799ts6z3kxd999tzp37qzWrVsrNja2WrE736RJk1RUVKTc3FxJ0qpVq9SxY0elpqaqZcuWgZ/WrVsrPT1dL774os6cOVPjsTp06KDJkydr/vz5gXVXX321rr76ahUVFWnGjBm6+uqr9eCDDwa2Z2dnq3HjxvrpT38a9H1D7fAYY4zTIYAfMnToUA0cOLDGD7EFADfyeDx6//33NWLECKej1JlXX31Vr732mrKzs+v0vFVVVUpOTtayZcs0fPjwOj03eE0kAAC4QpMmTdKkSZPq/LwRERHavXt3nZ8X3+HpbAAAAFjj6WwAAABY40okAAAArFEiAQAAYC3s3lgTFRWlZs2aOR0DQBj7+uuvdfbsWadjhAQzEsCVsJmPYVcimzVrFvjKJwC4HN9/48ePETMSwJWwmY88nQ0AAABrlEgAAABYo0QCAADAGiUSAAAA1iiRAAAAsEaJBAAAgDVKJAAAAKxRIgEAAGCNEgkAAABrlEgAAABYo0QCAADAGiUSAAAA1kJeIufMmaPExER5PB7l5uZesD0zM1Mej0cbN24MdRQAcBXmI4BwFvISeeedd2r79u1KSEi4YFthYaHWrFmjlJSUUMcAANdhPgIIZyEvkYMHD5bX671gfVVVldLS0rRixQpFRUWFOgYAuA7zEUA4c+w1kT6fTwMGDFBycrJTEQDAlZiPAMJBpBMnzcvL0/r167Vt27ZL7uvz+eTz+QLLZWVloYyGH4HEx7JDfo6CpbeF/Byon2zmo8SM/DFhdiHcOHIlcuvWrSosLFTHjh2VmJioHTt2aNq0aVq1atUF+6anp6ukpCTwEx0d7UBiAKgbNvNRYkYCcI4jJXLGjBkqLS1VQUGBCgoKlJKSotWrV2vGjBlOxAEA12A+AggXIS+R06dPl9frVUlJiUaNGqUOHTqE+pQAEBaYjwDCWchfE5mRkXHJfbZs2RLqGADgOsxHAOGMb6wBAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwFvISOWfOHCUmJsrj8Sg3N1eSVF5ernHjxqlTp07q1auXRo4cqYMHD4Y6CgC4CvMRQDgLeYm88847tX37diUkJFRbP23aNO3fv1979uzR2LFjNXXq1FBHAQBXYT4CCGchL5GDBw+W1+uttq5hw4YaPXq0PB6PJCklJUUFBQWhjgIArsJ8BBDOXPGayOXLl2vs2LFOxwAA12E+AnCrSKcDLFmyRPn5+dq0aVON230+n3w+X2C5rKysrqIBgKMuNR8lZiQA5zh6JfLZZ59VVlaW3n33XV1zzTU17pOenq6SkpLAT3R0dB2nBIC6F8x8lJiRAJzj2JVIn8+ndevWKScnR7GxsU7FAADXYT4CCAchvxI5ffp0eb1elZSUaNSoUerQoYNKSkr0yCOP6MSJExo2bJh69+6tm266KdRRAMBVmI8AwlnIr0RmZGTUuN4YE+pTA4CrMR8BhDNXvDsbAAAA4YUSCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYC3kJXLOnDlKTEyUx+NRbm5uYH1+fr769++vTp06qW/fvsrLywt1FABwFeYjgHAW8hJ55513avv27UpISKi2fvr06Zo2bZoOHDig+fPn6/777w91FABwFeYjgHAW8hI5ePBgeb3eauuOHj2qXbt2KTU1VZJ0xx13qLi4WAcPHgx1HABwDeYjgHDmyGsii4uL1apVK0VGRkqSPB6P4uPjVVRUdMG+Pp9PXq838FNWVlbXcQGgztjMR4kZCcA5rnljjTGmxvXp6ekqKSkJ/ERHR9dxMgBw1g/NR4kZCcA5jpTIuLg4lZaWyu/3S/puQBYXFys+Pt6JOADgGsxHAOHCkRLZvHlzJSUl6ZVXXpEkbdiwQV6vVx06dHAiDgC4BvMRQLgIeYmcPn26vF6vSkpKNGrUqMAgzMjIUEZGhjp16qSlS5cqMzMz1FEAwFWYjwDCWWSoT5CRkVHj+s6dO+ujjz4K9ekBwLWYjwDCmWveWAMAAIDwQYkEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArEU6HQD1S+Jj2U5HqBV1cT8Klt4W8nMAAHC5uBIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYokQAAALBGiQQAAIA1SiQAAACsUSIBAABgjRIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYokQAAALBGiQQAAIA1SiQAAACsUSIBAABgjRIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYcLZHvvPOOkpKS1Lt3b3Xv3l0vvfSSk3EAwDWYjwDcLtKpExtjlJqaqi1btqhnz54qKCjQ9ddfr/HjxysmJsapWADgOOYjgHDg6JVIj8ejEydOSJJOnTqlJk2aKCoqyslIAOAKzEcAbufYlUiPx6M33nhD48eP17XXXqvjx48rKytLV111lVORAMAVmI8AwoFjJdLv92vx4sXKysrS4MGDtXPnTt1+++3at2+fmjZtGtjP5/PJ5/MFlsvKypyICwB1Jtj5KDEj4T6Jj2WH/BwFS28L+TlwaUE/nT1u3Lig1gUrNzdXf/vb3zR48GBJUt++feX1erV79+5q+6Wnp6ukpCTwEx0dfdnnBIBQcGo+SsxIAM4JukQWFRVdsO6vf/3rZZ84Li5OpaWl+stf/iJJOnjwoA4dOqTOnTtf9jEBwAnMRwD10SWfzs7IyNDzzz+vAwcOKCkpKbD+5MmT6tat22WfuEWLFlq9erXuuusuRUREqKqqSitXrlR8fPxlHxMA6hLzEUB9dskSecstt6hz586aMWOG/v3f/z2w/rrrrlPPnj2v6OQTJkzQhAkTrugYAOAU5iOA+uySJTIhIUEJCQmBp1UAAN9hPgKoz4J+d3ZBQYF+85vf6NChQ/L7/YH1H3zwQUiCAUC4YD4CqI+CLpF33XWXbr75Zs2ePVsNGjQIZSYACCvMRwD1UdAlsry8XL/+9a9DmQUAwhLzEUB9FPRH/PTo0aPGj7EAgPqO+QigPgr6SuTRo0fVq1cv9evXTw0bNgysz8rKCkkwAAgXzEcA9VHQJTI1NVWpqamhzAIAYYn5CKA+CrpE3nfffaHMAQBhi/kIoD4KukQ+8MADNa5fu3ZtrYUBgHDEfARQHwVdIpOTkwN/Li8v14YNG6p9zRcA1FfMRwD1UdAlctasWdWWZ8yYoX/8x3+s9UAAEG6YjwDqo6A/4ufvNWzYUIcPH67NLADwo8B8BFAfBH0lMj09PfDnyspK7dq1S127dg1JKAAIJ8xHAPVR0CWyUaNG/3+jyEg99NBDuuOOO0ISCgDCCfMRQH0UdIlcuHBhKHMAQNhiPgKoj4J+TeTp06c1a9YsderUSZ06ddLs2bN1+vTpUGYDgLDAfARQHwVdImfOnCm/36/f/e53+t3vfqeqqirNnDkzlNkAICwwHwHUR0E/nb13717t2bMnsPzcc8+pV69eIQkFAOGE+QigPgr6SmRlZWW1p2dOnz6tysrKkIQCgHDCfARQH1l9d3ZKSoruvvtueTwevfHGG5oyZUoos6GOJT6W7XQEICwxH+sHZiRQ3SVL5KlTp3Ts2DE9+uij6t69uzZt2iRjjGbOnKnU1NS6yAgArsR8BFCfXfLp7Hnz5umTTz6RJN1666169tlntWzZMsXExGj+/PkhDwgAbsV8BFCfXbJEbt26tcYPzb333nu1devWkIQCgHDAfARQn12yRDZo0OCHbxxx2V+9DQBhj/kIoD675JTz+/06derUBetPnjypioqKkIQCgHDAfARQn12yRE6YMEH33nuvjh8/Hlh3/PhxTZkyRffcc09IwwGAmzEfAdRnlyyRCxYsUGxsrOLi4nTDDTfohhtuUFxcnGJiYvTkk0/WRUYAcCXmI4D67JIf8dOgQQO99NJL+pd/+Rd9+umnkqSkpCS1b98+5OEAwM2YjwDqs6A/bLx9+/YMRgCoAfMRQH3E2wcBAABgjRIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYokQAAALBGiQQAAIA1SiQAAACsUSIBAABgjRIJAAAAa5RIAAAAWHO0RJ49e1azZ89Wx44d1aNHD6WmpjoZBwBchRkJwM0inTz5Y489Jo/HowMHDsjj8ejIkSNOxgEAV2FGAnAzx0rkt99+qxdeeEElJSXyeDySpJYtWzoVBwBchRkJwO0cezr70KFDaty4sZYsWaI+ffpo0KBB2rRpk1NxAMBVmJEA3M6xK5EVFRUqLCxU165dtXTpUu3evVsjR45UXl6eWrRoEdjP5/PJ5/MFlsvKypyICwB1ihmJUEh8LNvpCPgRcexKZEJCgiIiIjRp0iRJ0g033KC2bdtq37591fZLT09XSUlJ4Cc6OtqJuABQp5iRANzOsRLZtGlT3XzzzfrTn/4kSTp8+LAOHz6sLl26OBUJAFyDGQnA7Rx9d/bzzz+vtLQ0zZ8/XxEREcrIyFCbNm2cjAQArsGMBOBmjpbIdu3aafPmzU5GAADXYkYCcDO+sQYAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMCaK0pkZmamPB6PNm7c6HQUAHAV5iMAt3K8RBYWFmrNmjVKSUlxOgoAuArzEYCbOVoiq6qqlJaWphUrVigqKsrJKADgKsxHAG7naIn0+XwaMGCAkpOTnYwBAK7DfATgdpFOnTgvL0/r16/Xtm3bLrqfz+eTz+cLLJeVlYU6GgA4Ktj5KDEjUT8lPpYd8nMULL0t5OcI9/vh2JXIrVu3qrCwUB07dlRiYqJ27NihadOmadWqVdX2S09PV0lJSeAnOjraocQAUDeCnY8SMxKAcxwrkTNmzFBpaakKCgpUUFCglJQUrV69WjNmzHAqEgC4AvMRQDhw/N3ZAAAACD+OvSby723ZssXpCADgSsxHAG7ElUgAAABYo0QCAADAGiUSAAAA1iiRAAAAsEaJBAAAgDVKJAAAAKxRIgEAAGCNEgkAAABrlEgAAABYo0QCAADAGiUSAAAA1iiRAAAAsEaJBAAAgDVKJAAAAKxRIgEAAGCNEgkAAABrlEgAAABYo0QCAADAGiUSAAAA1iiRAAAAsBbpdAAAgHMSH8sO+TkKlt4W8nMAqHtciQQAAIA1SiQAAACsUSIBAABgjRIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYokQAAALBGiQQAAIA1SiQAAACsUSIBAABgjRIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYokQAAALBGiQQAAIA1SiQAAACsUSIBAABgjRIJAAAAa46VyPLyco0bN06dOnVSr169NHLkSB08eNCpOADgKsxIAG7n6JXIadOmaf/+/dqzZ4/Gjh2rqVOnOhkHAFyFGQnAzRwrkQ0bNtTo0aPl8XgkSSkpKSooKHAqDgC4CjMSgNtFOh3ge8uXL9fYsWMvWO/z+eTz+QLLZWVl1sdOfCz7irIFo2DpbSE/R13cD7gHv+/g1MVjzw1COSMBXIgZfGmuKJFLlixRfn6+Nm3adMG29PR0paenB5a9Xm9dRgMAxzEjAbiR4yXy2WefVVZWlnJycnTNNdc4HQcAXIUZCcCtHC2RPp9P69atU05OjmJjY52MAgCuw4wE4GaOlciSkhI98sgjateunYYNGyZJioqK0scff+xUJABwDWYkALdzrER6vV4ZY5w6PQC4GjMSgNvxjTUAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArEU6HeDHIvGxbKcjAIArMR+BHyeuRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokAAAArFEiAQAAYI0SCQAAAGuUSAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUdLZH5+vvr3769OnTqpb9++ysvLczIOALgKMxKAmzlaIqdPn65p06bpwIEDmj9/vu6//34n4wCAqzAjAbiZYyXy6NGj2rVrl1JTUyVJd9xxh4qLi3Xw4EGnIgGAazAjAbidYyWyuLhYrVq1UmRkpCTJ4/EoPj5eRUVFTkUCANdgRgJwu0inA5zPGHPBOp/PJ5/PF1g+cuSIvF5vXcaqUVlZmaKjo52OUSu4L+7EfQmO9xX723z99de1H6QOuGFGuvW/S7fmktybza25JPdmc2su6Yez2c5Im/noMTVNpTpw9OhRdejQQceOHVNkZKSMMWrVqpW2b9+uDh06OBHJitfrVUlJidMxagX3xZ24L/WbW2ekW3+Xbs0luTebW3NJ7s3m1lySM9kcezq7efPmSkpK0iuvfFeRN2zYIK/XGxYFEgBCjRkJwO0cfTo7IyND999/v5YsWaLrrrtOmZmZTsYBAFdhRgJwM0dLZOfOnfXRRx85GeGypaenOx2h1nBf3In7AjfOSLf+Lt2aS3JvNrfmktybza25JGeyOfaaSAAAAIQvvvYQAAAA1iiRAAAAsEaJvALLly9X9+7d1aNHD/Xs2TPwLspwlJ2drT59+igqKkpz5851Oo61H9N3DM+ZM0eJiYkKrKWiAAAJ1klEQVTyeDzKzc11Os4VKS8v17hx49SpUyf16tVLI0eO5BtXwkhVVZUeeughtW/fXh06dNDKlSt/cN+zZ89q9uzZ6tixo3r06BH4ph03ZPteZmamPB6PNm7c6HiuunpsBDsbX3jhBXXs2FHt27fXgw8+qIqKilrPcjnZPvjgA914443q2rWrunXrpnnz5qmqqsrxXN8zxmj48OGKjY0NaSbbbPv27dPQoUPVpUsXdenSRVlZWaEJZHDZcnJyzIkTJ4wxxhQVFZkmTZqYgwcPOpzq8uzfv9/k5uaaBQsWmIcfftjpONaGDRtmMjMzjTHG/P73vzd9+vRxNtAV+PDDD01xcbFJSEgwu3fvdjrOFTlz5ozJzs42VVVVxhhjVqxYYYYMGeJsKATtpZdeMsOHDzd+v9988803Jj4+3nz22Wc17jt37lwze/bswO+6tLTUNdmMMaagoMD069fPpKSkmDfffNPxXHX12AhmNv71r381rVq1MqWlpaaqqsr87Gc/MytXrqz1LJeT7dNPPzWHDh0yxnz3dzZgwIDAbZzM9b1ly5aZqVOnmkaNGoU0k022b7/91rRt29Zs27bNGGOM3+83X3/9dUjyUCJrUbdu3czmzZudjnFFFi5cGHYl8quvvjIxMTGmoqLCGGNMVVWVadGihcnPz3c42ZX5MZTIv7dz506TkJDgdAwEafTo0WbdunWB5UcffdQsWLDggv3KyspMTEyMOXnypOuyGWNMZWWlufnmm82uXbvMkCFDQloibXKdLxSPjWBn4zPPPGOmT58eWM7OzjYDBgyo1SyXm+3vzZo1yyxcuNAVuT777DMzaNAgk5+fXyclMthsa9asMRMmTAh5HmOM4ensWpKTk6Pjx4+rb9++Tkepd/iO4fCxfPlyjR071ukYCFJRUZESEhICy4mJiTU+rg4dOqTGjRtryZIl6tOnjwYNGqRNmza5Ipv03VdDDhgwQMnJySHNZJvrfKF4bAQ7Gy83c11kO9+RI0e0fv16jRkzxvFcFRUVevDBB5WRkRHYN9SCzfb5558rKipKY8aMUe/evTV58uSQfdWrq74722369eun/Pz8Grft3r1bcXFxkr577cGUKVP0xhtv6Nprr63LiEEL9r78WBg+ucp1lixZovz8/JCXCwTvUnPh7/3Q46qiokKFhYXq2rWrli5dqt27d2vkyJHKy8tTixYtHM2Wl5en9evXa9u2bZeVI1S5zleXj41g8jg1Py923lOnTulnP/uZ5s2bpz59+tRhqppzPfXUUxo/fry6dOmigoKCOs1zvpqyVVRUKCcnRzt27FDr1q31+OOPa8aMGVq/fn2tn58SeRHBfMjv559/rjFjxmjt2rUaOHBgHaS6PG77wOLaFBcXp9LSUvn9/sB3DBcXFys+Pt7paPg/zz77rLKyspSTk6NrrrnG6Tj4P5eaC/Hx8SosLFS/fv0kSYWFhTU+rhISEhQREaFJkyZJkm644Qa1bdtW+/btu+wSWVvZtm7dqsLCQnXs2FHSd1ezpk2bptLSUs2YMcOxXN8L5WMj2NkYHx+vQ4cOBZYvlbkus0nS6dOndcstt+j2228P+QdqB5vrww8/VFFRkVauXCm/369Tp04pMTFRO3fuVLNmzRzNlpCQoGHDhqlNmzaSpNTUVI0aNSokmXhN5BX4/PPPTUJCgnnvvfecjlJrwvE1kcYYM2TIkGovNk5OTnY2UC34sbwmctmyZSYpKckcO3bM6SiwlJmZGXiTyLFjx0x8fLzZu3dvjfuOHDnSZGdnG2O+e6NGkyZNTElJiSuynS/Ur4m0yVUXj41gZuOhQ4cueGPNihUrQpbJJtvp06dN//79zaJFi0KexybX+Q4fPlxnb6wJJlthYaG5/vrrA69RfuaZZ8zo0aNDkocSeQVGjBhhYmNjTa9evQI/4Vooc3JyTJs2bUxMTIyJjo42bdq0MW+99ZbTsYL2xRdfmJSUFNOxY0eTnJwc1D8mbjVt2jTTpk0b06BBA9O8eXPTvn17pyNdtuLiYiPJtGvXLvAYufHGG52OhSD5/X4zc+ZM07ZtW9OuXTvz29/+NrDtrbfeMmlpaYHlQ4cOmaFDh5ru3bubnj17mvXr17sm2/lCXSKDzVVXj40fmo1paWnVZvzq1atNu3btTLt27cwDDzxgzp07V+tZLifb4sWLTWRkZLV/ZxcvXux4rvPVZYkMNtvLL79sunXrZnr06GFuueUWU1RUFJI8fO0hAAAArPHubAAAAFijRAIAAMAaJRIAAADWKJEAAACwRokEAACANUokXOPWW2/VypUrL1jfq1cvvfnmmz94u0WLFmnu3LmhjAYAjmI+wo0okXCNtLQ0ZWZmVlu3a9cuHTlyJKTflQoAbsd8hBtRIuEat99+u4qLi7Vnz57AurVr12ry5Mn64osvNGjQICUlJalr165avHhxjcd48cUXNW7cuMDy22+/raFDhwaW/+u//ks33XSTkpKSNHjw4MC5duzYoeTkZPXu3Vvdu3fXqlWrQnMnAeAyMB/hRnx3NlzjqquuUmpqqjIzM/Xb3/5W5eXlev311/XnP/9ZXq9XOTk5ioqK0pkzZ9S/f3+NGDFCKSkpQR//z3/+s9atW6etW7cqKipK27Zt08SJE5WXl6df//rX+uUvf6kJEyZIko4fPx6quwkA1piPcCNKJFwlLS1NQ4cO1TPPPKOsrCx16dJFXbp00dGjRzVz5kzl5uYqIiJCxcXFys3NtRqSb731lvbs2aObbropsO7YsWM6c+aMhg0bpqefflr5+fkaPny4Bg4cGIq7BwCXjfkIt+HpbLhKt27d1L59e/3xj3/U2rVrlZaWJkl6/PHH1aRJE+3evVt79uzR0KFDVV5efsHtIyMjVVlZGVg+fx9jjO677z7l5uYGfkpLS3X11Vdr7ty5evvtt9WqVSs9/vjjmjlzZujvLABYYD7CbSiRcJ20tDQtWbJEO3fu1F133SXpu6dP4uLiFBkZqf379+v999+v8bbt27fX3r17debMGfn9fr322muBbWPHjtWrr76qoqIiSVJVVZV27dolSdq/f7/atWunBx98UI8//rh27NgR4nsJAPaYj3ATns6G69xzzz36xS9+obvvvlvR0dGSpCeeeEL33nuvXn/9dSUmJmr48OE13rZfv34aPXq0evToobZt2yopKUkff/yxJGngwIH6zW9+o5///Ofy+/06d+6cbrvtNvXp00crVqzQ5s2bddVVV6lBgwZatmxZnd1fAAgW8xFu4jHGGKdDAAAAILzwdDYAAACsUSIBAABgjRIJAAAAa5RIAAAAWKNEAgAAwBolEgAAANYokQAAALBGiQQAAIA1SiQAAACs/S9I8hgZYiwYKwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f05a52ce278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(num=None, figsize=(15, 6), dpi=80, facecolor='w', edgecolor='k')\n",
    "# independent\n",
    "ax1 = plt.subplot(131)\n",
    "ax1.hist(lXIndependent)\n",
    "ax1.set_xlabel('Values')\n",
    "ax1.set_ylabel('Count')\n",
    "ax1.set_title('independent')\n",
    "\n",
    "# dependent (AR1)\n",
    "ax2 = plt.subplot(132)\n",
    "ax2.hist(lXDependent)\n",
    "ax2.set_xlabel('Values')\n",
    "ax2.set_ylabel('Count')\n",
    "ax2.set_title('dependent (AR1)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The reason for the poor fit is the dependence in the AR1 process with $\\rho = 0.99$\n",
    "\n",
    "Use spectral power density to calculate autocorrelation at all lags simultaneously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Autocorrelation(x) :\n",
    "    xp = x-np.mean(x)\n",
    "    f = np.fft.fft(xp)\n",
    "    p = np.array([np.real(v)**2+np.imag(v)**2 for v in f])\n",
    "    pi = np.fft.ifft(p)\n",
    "    return np.real(pi)[:x.size//2]/np.sum(xp**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqAAAAGrCAYAAAAfEyWbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAIABJREFUeJzs3X90VPWd//HXzSQCISAJAUEmyQgh2LUgYoIIKFLrqlWjK1pdipAiP3Rt7W7cSrcsR+laym638VipJYLBTanoFqzNatutIKX4FZvEH7WLVhIlyUSREIhKCD8mM5/vH5gpMSGZhJk7d2aej3PmHO7cT+a+7wx55z2f+/l8rmWMMQIAAABskhTtAAAAAJBYKEABAABgKwpQAAAA2IoCFAAAALaiAAUAAICtKEABAABgKwpQAAAA2IoCFBF37bXX6t/+7d/6/fMzZ87Ugw8+GL6AIiRW4gQQWe3t7bIsS7///e+jHUqPoh3nz3/+c11//fW2HzcQCGjy5MmO/3ziHQUoIu43v/mNVqxYEe0wYs7vf/97WZal9vb2aIcCACF58MEHNXPmzF7bHTt2TN/+9rf10EMPddm3fft2WZal22+/vdvXd7lcSktL05AhQ5STk6N/+qd/0vHjx4NtPvjgA914443KycmRZVlav359p9dISkrSypUrde+99/bjDBEuFKAAAMBWmzZtktvt1uTJk7vse+yxxzR8+HA9++yz2r9/f5f9l156qVpbW3X48GH9z//8j5566imtXr06uD8pKUl/+7d/q6eeekput7vb419//fVqbm7Wiy++GL6TQp9QgCLirrjiCv3rv/6rJMmyLD366KOaMWOG0tLSNHHiRL388svBtu3t7br//vs1atQojRgxQv/yL//S5fU++OADzZ07V2PGjNHIkSP193//9zpw4ECn433jG9/QzTffrCFDhmj8+PEqLy/v9Bp//OMfdcUVV2j48OHKycnRihUrOvU02hXnt771Lc2dO1dnn322srKy9NOf/lSS1NDQoGuvvVaSNGzYMKWlpWnVqlV9et8B2KOpqUk333yzhg0bprFjx+qZZ57p0iaUnFNSUqKCggKlpaVp6tSpqq6u7vQa5eXluvDCC3X22Wfrggsu0NNPPx3c13HFZPPmzcrLy9OQIUN01VVX6YMPPgh7nKfLjT//+c+1atUq7dq1S2lpaUpLS9POnTu7fc+2bNmiq6++usvz+/bt03PPPadHH31UGRkZWrdu3enedknSpEmTdNlll6mqqir43OjRo3XPPfdoxowZcrlc3f6cy+XSl7/8ZT377LM9vj4iyAARNmvWLLN8+XJjjDGSzKRJk0xNTY3x+XzmW9/6lsnOzg62feihh0xOTo7ZvXu3OXbsmPnXf/1Xk5ycbB544AFjjDHHjh0zEyZMMPfdd59pbW01hw8fNvPmzTNf/vKXOx1v4MCBpqKiwvh8PvPCCy+YlJQU8/LLLxtjjPnLX/5i0tLSzKZNm4zP5zN1dXVm0qRJ5qGHHgq+hl1xDh061Gzbts34/X6zefNmk5SUZGpqaowxxmzfvt1IMj6fL7wfCICwuuqqq8zVV19tDh48aA4ePGiuu+46I8ls377dGBN6zhk3blwwpzzwwAMmMzPTfPzxx8YYYzZs2GCysrJMVVWV8fv9ZufOnWbIkCFm586dxpi/5ou5c+eajz/+2Hz88cdm+vTpZv78+WGPs6fc+MADD5gZM2b0+p6dc8455umnn+7y/IMPPmhGjhxpTpw4YZYtW2aysrJMe3t7t68fCATMa6+9ZjIzM80///M/d3ucnJwcs27dum73/fCHPzQXX3xxr7EiMihAEXGfL0D/67/+K7jv//7v/4wk89FHHxljjMnNzTU//vGPg/vb29vNiBEjgoXdli1bzLnnnmsCgUCwTWNjo5FkvF5v8Hg333xzpxi++tWvmoULFxpjjPnmN79p/v7v/77T/o0bN5px48YFt+2K8+tf/3qnODIzM4NJmQIUcL6O3+u33nor+Nxbb73VqbALNeecmlP8fr8ZPXq0KS8vN8YYM3HiRFNaWtrpNRYtWmTuvPNOY8xf80V9fX1w/5o1a8z5558f9jh7yo2hFqApKSnmN7/5Tafn2tvbzZgxY8z9999vjDGmtrbWWJZlnnvuuWCbBx54wLhcLnP22WebAQMGGEnm9ttvN62trd0ep6cC9PHHHzfnnXder7EiMpLt62sFTjr33HOD/x48eLAk6fDhwzrnnHPU2NiosWPHBve7XC5lZ2cHt2tqarR//36lp6d3es0BAwaooaEhON7n1NeQpPPOO0+vv/568DW2b9+uYcOGBfcHAgEFAgHb4zz1GB3HOXz4sADEhsbGRkmdc855553XqU2oOefU10hKSlJ2dra8Xm/wNYqLi3X//fcH27S3t+vyyy/v9Bqfz1sd+SSccfaUG0OVkZGhTz75pNNzv/rVr/Thhx9q0aJFkqRx48Zp9uzZeuyxx3TjjTcG202bNk0vv/yy2tvbVVZWpu9973tqaWkJxhKqTz75RBkZGX36GYQPBSgcxe12a+/evcFtv98fTMCSNGrUKOXk5Oi9997r8XVOfY2O7Y6ib9SoUZo7d67KysqiHmdPkpIYog04XUde2bt3r774xS8G/32qUHPOqT8XCAQ6fVkdNWqUVq5cqfnz50c9zp6Emrcuvvhi7d69u9Nzjz32mCTpsssuCz7X2tqqtrY21dbWKjc3t1P75ORkLVmyRDt27NA3v/lN/fKXv+xTrH/+85+Vn5/fp59B+PAXDo6yYMEC/ehHP9Jf/vIXHT9+XN/73vd06NCh4P6bb75ZPp9PK1asCH57bmpq6jKY/te//rVeeOEF+f1+/fa3v9Vzzz2nr3/965Kkf/iHf9DmzZv1i1/8QidOnJDf71dtba1++9vf2h5nT0aNGiVJevfdd0P+GQD2GjNmjK688krdf//9amlpUUtLi5YvX96pTag555FHHtE777yjEydO6Pvf/75OnDihwsJCSdI//uM/6t/+7d9UVVWlQCCg48ePq6qqSq+99prtcfZk1KhRamho0LFjx3psd/PNN+t///d/g9t79uzRSy+9pCeffFJvvvlm8LFnzx653e7gBM3urFy5Us8//3yniaLHjh3TsWPHZIxRe3u7jh07Jp/PF9zv9/u1detW/d3f/V3I54bwogCFoyxbtkw333yzZs2aJbfbrRMnTuiSSy4J7h8yZIh27dqlhoYGTZw4UUOHDtX06dP1hz/8odPrLFy4UE888YSGDRume+65R2vXrg1+qy4oKNCLL76odevWacyYMRo+fLhuueUW1dfX2x5nT/Ly8nTvvfdq9uzZGjZsWKdlRgA4x8aNG3XWWWfJ4/FoypQpuu222zrtDzXn3H333brjjjuUkZGhiooK/frXvw5eDv/Wt76lBx98UHfddZcyMjI0ZswYffvb39aRI0dsj7Mnt912myZMmKBzzz1Xw4YN61QUnuprX/uaGhoa9Oabb0qSfvrTn2r8+PGaN2+eRo0aFXyce+65Ki4u1pNPPqmjR492+1q5ubmaP3++li1bFnxu0KBBGjRokBoaGnT33Xdr0KBBWrx4cXD/Cy+8oIyMDP3t3/5tyOeG8LKMMSbaQQDhdMUVV2jmzJndLnAMAE5kWZZefPFFffnLX452KLb5+c9/rqeeekovvPCCrccNBAK6+OKL9aMf/Uhf+tKXbD02/ooxoAAAwHZf+9rX9LWvfc324yYlJemNN96w/bjojEvwAAAAsBWX4AEAAGArekABAABgKwpQAAAA2CohJiENGDBAI0aMiHYYAGLYgQMHdPz48WiHERHkSABnoj/5MSEK0BEjRgRvQwYA/dFxJ5l4RI4EcCb6kx+5BA8AAABbUYACAADAVhSgAAAAsBUFKAAAAGxFAQoAAABbUYACAADAVhSgAAAAsBUFKAAAAGxFAQoAAABbUYACAADAVhSgAAAAsJUjCtB7771XHo9HlmXpzTffPG27J554QuPHj9e4ceO0ePFi+Xy+sMdijFFV3SH9otqrqrpDMsaE/RgAECryI4B45IgC9JZbbtHLL7+snJyc07bZu3evVqxYoZ07d6q2tlb79+/X448/HtY4GlvadGXJDs1d96oeqNituete1ZUlO9TY0hbW4wBAqMiPAOKRIwrQyy+/XG63u8c2mzdvVmFhoUaNGiXLsnTXXXdp06ZNYYvBGKP5ZZWqP9gmn9+o7YRfPr9R/cE2LSir5Js+gKggPwKIR44oQEPR0NDQqQfA4/GooaEhbK9fXd+ixkNH5Q90TqT+gFHDoTZV17eE7VgAEE7kRwCxJmYK0M/r6Rt3SUmJ3G538NHa2trr69U1H1Gyy+p2X4orSXXNR/odKwDYqbceyb7mSPIjgHCLmQI0Oztb9fX1we36+nplZ2d327a4uFiNjY3BR1paWq+v78kcLJ8/0O0+nz8gT+bg/gUOABHWl/wo9T1Hkh8BhFvMFKBz5sxRRUWFPvroIxljtHbtWt1+++1he/38nHRlZaTKldT5W74ryVJ2Rqryc9LDdiwACCfyI4BY44gCdOnSpXK73WpsbNTVV1+t3NxcSdKiRYtUUVEhSRo7dqxWrlypGTNmKDc3VyNGjNDSpUvDFoNlWSpfOFU5w1M/25ZSXJY8w1NVfuclsqzuLz8BQCSRHwHEI8skwPTFjuQdCmOMLvuP7fL5A1ozd4ryc9JJrgD6lEdiTajnRn4E0J3+5MfkCMUSsyzL0qAUlwaluFTgyYh2OADgGORHAOHiiEvwAAAASBwUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABslRztAAAA8cUYo+r6FtU1H5Enc7Dyc9JlWVa0wwLgIBSgAICwaWxp0/yySnkPtSnFlSSfP6CsjFSVL5wqd3pqtMMD4BBcggcAhIUxRvPLKlV/sE0+v1HbCb98fqP6g21aUFYpY0y0QwTgEBSgAICwqK5vUeOho/IHOhea/oBRw6E2Vde3RCkyAE5DAQoACIu65iNKdnU/1jPFlaS65iM2RwTAqShAAQBh4ckcLJ8/0O0+nz8gT+ZgmyMC4FQUoACAsMjPSVdWRqpcSZ17QV1JlrIzUpWfkx6lyAA4DQUoACAsLMtS+cKpyhme+tm2lOKy5BmeqvI7L2EpJgBBFKAA4HA1NTWaPn268vLyVFBQoN27d3dpEwgEVFxcrL/5m7/RpEmTNHv2bNXW1toeqzs9VduKZ8mdPkgjhwzQU4unaWvxLI0ZNsj2WAA4FwUoADjc0qVLtWTJEu3Zs0fLli1TUVFRlzYVFRX6f//v/+lPf/qT3nrrLV155ZX67ne/a3+wOtkTOijFpaEDU1TgyaDnE0AXFKAA4GBNTU2qrq7WvHnzJElz5syR1+vt0rtpWZaOHz+uY8eOyRijTz/9VG63OxohA0CvuBMSADiY1+vV6NGjlZx8Ml1blqXs7Gw1NDQoNzc32O6GG27Q9u3bNWrUKA0ZMkRjxozRjh07ohU2APSIHlAAiDHd3VGourpa//d//6cPPvhAH374oa688krddddd3f58SUmJ3G538NHa2hrpkAGgEwpQAHCwrKws7du3T+3t7ZJOFp9er1fZ2dmd2pWXl+tLX/qShg0bpqSkJC1YsEDbt2/v9jWLi4vV2NgYfKSlpUX8PADgVBSgAOBgI0eO1JQpU7Rx40ZJ0pYtW+R2uztdfpeksWPH6qWXXtKJEyckSc8//7y++MUv2h4vAITCMQVoLC0zAgB2Ki0tVWlpqfLy8rR69Wpt2LBBkrRo0SJVVFRIku655x6dd955mjRpkiZNmqRt27bppz/9aTTDBoDTcswkpI5lRoqKirR582YVFRWpqqqqU5tTlxlJSUnRQw89pO9+97v67//+7yhFDQCRN2HCBO3atavL8+vXrw/+e8CAAVq3bp2dYQFAvzmiB5RlRgAAABKHI3pAWWYEAAAgcTiiB7Q7Z7LMCEuMAAAAOJcjCtBwLzPCEiMAAADO5YgClGVGAAAAEocjxoBKJ5cZKSoq0qpVqzR06NBOy4wUFhaqsLBQ99xzj9555x1NmjRJZ511lkaNGqW1a9dGOXIAAAD0hWMKUJYZAQAASAyOuAQPAACAxOGYHlAAQOIxxqi6vkV1zUfkyRys/Jx0WZYV7bAARBgFKAAgKhpb2jS/rFLeQ21KcSXJ5w8oKyNV5Qunyp2eGu3wAEQQl+ABALYzxmh+WaXqD7bJ5zdqO+GXz29Uf7BNC8oqu10LGkD8oAAFANiuur5FjYeOyh/oXGj6A0YNh9pUXd8SpcgA2IECFABgu7rmI0p2dT/WM8WVpLrmIzZHBMBOFKAAANt5MgfL5w90u8/nD8iTOdjmiADYiQIUAGC7/Jx0ZWWkypXUuRfUlWQpOyNV+TnpUYoMgB0oQAEAtrMsS+ULpypneOpn21KKy5JneKrK77yEpZiAOMcyTACAqHCnp2pb8Sxd9h/b5fMHtGbuFNYBBRIEBSgAIGosy9KgFJcGpbhU4MmIdjgAbMIleAAAANiKAhQAAAC2ogAFAACArShAAQAAYCsKUAAAANiKAhQAAAC2ogAFAACArShAAQAAYCsKUAAAANiKAhQAAAC2ogAFAACArShAAQAAYCsKUAAAANiKAhQAAAC2ogAFAACArShAAQAAYCsKUAAAANiKAhQAAAC2ogAFAACArShAAcDhampqNH36dOXl5amgoEC7d+/utt2f//xnXXHFFfrCF76gL3zhC3r22WdtjhQAQpMc7QAAAD1bunSplixZoqKiIm3evFlFRUWqqqrq1KatrU033nijysvLNXPmTPn9frW0tEQpYgDoGT2gAOBgTU1Nqq6u1rx58yRJc+bMkdfrVW1tbad2Tz31lKZNm6aZM2dKklwulzIzM22PFwBCQQEKAA7m9Xo1evRoJSefvGBlWZays7PV0NDQqd3bb7+tAQMG6Prrr9fkyZM1f/58HThwoNvXLCkpkdvtDj5aW1sjfh4AcCoKUACIMcaYLs/5fD5t3bpVpaWleuONNzRmzBjdfffd3f58cXGxGhsbg4+0tLRIhwwAnVCAAoCDZWVlad++fWpvb5d0svj0er3Kzs7u1C4nJ0ezZ8/WmDFjZFmW5s2bp1dffTUaIQNAryhAAcDBRo4cqSlTpmjjxo2SpC1btsjtdis3N7dTu69+9auqqqrSp59+Kkn69a9/rQsvvND2eAEgFMyCBwCHKy0tVVFRkVatWqWhQ4dqw4YNkqRFixapsLBQhYWFys7O1ne/+11Nnz5dSUlJGjNmjB5//PEoRx4+xhhV17eorvmIPJmDlZ+TLsuyoh0WgH6iAAUAh5swYYJ27drV5fn169d32r7jjjt0xx132BWWbRpb2jS/rFLeQ21KcSXJ5w8oKyNV5Qunyp2eGu3wAPQDl+ABAI5ljNH8skrVH2yTz2/UdsIvn9+o/mCbFpRVdjshC4DzUYACAByrur5FjYeOyh/oXGj6A0YNh9pUXc9i+0AsogAFADhWXfMRJbu6H+uZ4kpSXfMRmyMCEA6OKUC51zEA4PM8mYPl8we63efzB+TJHGxzRADCwTGTkLjXMQDg8/Jz0pWVkar6g22dLsO7kixlZ6QqPyc9itEB6C9H9IByr2MAQHcsy1L5wqnKGZ762baU4rLkGZ6q8jsvYSkmIEY5ogCNxL2OAQDxwZ2eqm3Fs+ROH6SRQwboqcXTtLV4lsYMGxTt0AD0kyMK0O6cyb2OS0pK5Ha7g4/W1lY7QgYARIhlWRqU4tLQgSkq8GTQ8wnEOEcUoOG+13FxcbEaGxuDj7S0NFvOAwAAAL1zRAHKvY4BAAASh2NmwXOvYwAAgMTgmAI00e91DAAAkCgccQkeAAAAiYMCFAAAALaiAAUAAICtKEABAABgKwpQAAAA2IoCFAAAALaiAAUAAICtKEABAABgKwpQAAAA2IoCFAAAALaiAAUAAICtKEABAABgKwpQAAAA2IoCFAAAALaiAAUAAICtkqMdAAAA4WKMUXV9i+qaj8iTOVj5OemyLCvaYQH4nJAL0H379mnv3r1qb28PPnf55ZdHJCgAiCXkR2dobGnT/LJKeQ+1KcWVJJ8/oKyMVJUvnCp3emq0wwNwipAK0O9///v64Q9/qLFjx8rlckmSLMtSZWVlRIMDAKcjPzqDMUbzyypVf7BN/oCRz++XJNUfbNOCskptLZ5FTyjgICEVoGVlZXrvvfc0fPjwSMcDADGF/OgM1fUtajx0VP6A6fS8P2DUcKhN1fUtKvBkRCk6AJ8X0iSkc845h+QKAN0gPzpDXfMRJbu67+FMcSWprvmIzREB6ElIPaBXXXWV/vEf/1Fz587VwIEDg89PmjQpYoEBQCwgPzqDJ3OwfP5At/t8/oA8mYNtjghAT0IqQMvLyyVJv/rVr4LPWZal999/PzJRAUCMID86Q35OurIyUoNjQDu4kixlZ6QqPyc9itEB+LyQCtC9e/dGOg4AiEnkR2ewLEvlC6dqflml3j9wRJYlJX9WfJbfeQkTkACHCXkZpsrKSm3dulXSyUtOBQUFEQsKAGIJ+dEZ3Omp2lY8S5f9x3b5/AGtmTuFdUABhwppEtLjjz+uW265RU1NTWpqatItt9yi9evXRzo2AHA8O/JjTU2Npk+frry8PBUUFGj37t2nbWuM0Ze+9CUNGzYsrDHECsuyNCjFpaEDU1TgyaD4BBwqpB7QNWvW6LXXXtOIESMkScuXL9eVV16pRYsWRTQ4AHA6O/Lj0qVLtWTJEhUVFWnz5s0qKipSVVVVt20ffvhhjRs3Tq+//nrYjg8A4RbyveA7kuvn/w0AiS6S+bGpqUnV1dWaN2+eJGnOnDnyer2qra3t0nb37t167rnntGzZsrDGAADhFlIBOn78eC1fvlwNDQ3yer1asWKFxo8fH+nYAMDxIp0fvV6vRo8ereTkkxesLMtSdna2GhoaOrXz+XxavHixSktLg21Pp6SkRG63O/hobW0NW7wAEIqQCtC1a9fqvffe05QpU3TRRReptrZWP/3pTyMdGwA4XjTyozGmy3MrV67UzTffrC984Qu9/nxxcbEaGxuDj7S0tEiECQCnFdIY0BEjRujpp5+OdCwAEHMinR+zsrK0b98+tbe3Kzk5WcYYeb1eZWdnd2q3Y8cONTQ0aM2aNWpvb9enn34qj8ejqqoqhk0BcJweC9AdO3Zo1qxZqqio6HZ/YWFhRIICAKezKz+OHDlSU6ZM0caNG1VUVKQtW7bI7XYrNze3U7udO3cG/11XV6fJkyerrq4uLDEAQLj1WIBu3LhRs2bN0sMPP9xln2VZFKAAEpad+bG0tFRFRUVatWqVhg4dqg0bNkiSFi1apMLCQnIxgJjTYwG6bt06SdL27dttCQYAYoWd+XHChAnatWtXl+dPt96ox+PRxx9/HOmwAKDfQpqENHXq1JCeA4BEQ34EgL4LqQBtb2/vtO3z+XT48OGIBAQAsYT8CAB912MB+u///u9KT0/Xn//8Z2VkZAQfQ4cO1eWXX25XjADgOORHAOi/HseA3nXXXbrtttt09913a+3atcHnhw4dqvT09IgHBwBORX4EgP7rsQA9++yzdfbZZ+s3v/mNXfEAQEwgPwJA/4W0EH1TU5MeeOAB/elPf9KxY8eCz7/++usRCwwAYgH5EQD6LqRJSHfeeac8Ho+am5u1cuVKnXvuubruuusiHRsAOB75EQD6LqQC1Ov1atmyZRowYIBuuOEGPfvss3rllVciHRsAOB75EQD6LqQC9KyzzpIkDRw4UAcPHlRycrIaGxsjGhgAxALyIwD0XUhjQCdMmKBDhw5p3rx5uuSSSzR06FBddNFFkY4NAByP/AgAfRdSAfqzn/1MkvStb31LF198sT7++GNde+21YQ2kpqZGCxYsUHNzs84++2w9+eSTuuCCC7pta4zRlVdeqddff53bzQGIKjvyIwDEm5AK0FPNnDkzEnFo6dKlWrJkiYqKirR582YVFRWpqqqq27YPP/ywxo0bxyxTAI4SqfyIyDDGqLq+RXXNR+TJHKz8nHRZlhXtsICE0GMBmp7e/S+jMUaWZenQoUNhCaKpqUnV1dX63e9+J0maM2eOvvGNb6i2tla5ubmd2u7evVvPPfecysrK9Itf/CIsxweAvrIrPyIyGlvaNL+sUt5DbUpxJcnnDygrI1XlC6fKnZ4a7fCAuNdjAfrmm2/aEoTX69Xo0aOVnHwyHMuylJ2drYaGhk4FqM/n0+LFi/XEE08E2wJANNiVHxF+xhjNL6tU/cE2+QNGPr9fklR/sE0Lyiq1tXgWPaFAhPVYxeXk5AT/vW/fPr377ru64oor1N7erkAgENHAjDFdnlu5cqVuvvlmfeELX1BdXd1pf7akpEQlJSXB7dbW1kiECCCBRTM/4sxU17eo8dBR+QOd/874A0YNh9pUXd+iAk9GlKIDEkNIyzBt3rxZ06ZNU1FRkaSTl8FvuummsAWRlZWlffv2qb29XdLJ4tPr9So7O7tTux07dujRRx+Vx+PRzJkz9emnn8rj8ejAgQOd2hUXF6uxsTH4SEtLC1usAHCqSOdHhF9d8xElu7rv4UxxJamu+YjNEQGJJ6QC9Ac/+IFef/11paenS5IuvPBC1dfXhy2IkSNHasqUKdq4caMkacuWLXK73V3Gf+7cuVP19fWqq6vTyy+/rKFDh6qurk4jRowIWywA0BeRzo8IP0/mYPn83fdS+/wBeTIH2xwRkHhCKkBdLpeGDx/e6bmOxZfDpbS0VKWlpcrLy9Pq1au1YcMGSdKiRYtUUVER1mMBQLjYkR8RXvk56crKSJUrqXMvqCvJUnZGqvJz0qMUGZA4QprJM2TIEO3fvz84KHvbtm3KyAjv+JgJEyZo165dXZ5fv359t+09Hg9rgAJjIaGwAAAgAElEQVSIOjvyI8LLsiyVL5yq+WWVev/AEVmWlPxZ8Vl+5yVMQAJsEFIBunr1al177bV6//33NXPmTO3du1cvvPBCpGMDAMcjP8Ymd3qqthXP0mX/sV0+f0Br5k5hHVDARiEVoAUFBdq+fbteeeUVGWM0ffp0DRs2LNKxAYDjkR9jl2VZGpTi0qAUF7PeAZv1WoD6/X5ddNFFeuutt7i9HACcgvwIAP3T6yQkl8slt9uttrY2O+IBgJhBfgSA/gnpEvzYsWM1Y8YM3XrrrZ3W1Lz33nsjFhgAxALyIwD0XUgF6JEjRzR58mTV1NQEn2OgNgCQHwGgP0IaA3r99ddrzpw5dsQDADGD/AgA/RPSGNBT76sOADiJ/AgA/RPSnZAKCgr08ssvRzoWAIg55EcA6LuQxoBu27ZNa9as0dixYzsNsn/99dcjFhgAxALyIwD0XUgF6Jo1ayIdBwDEJPIjAPRdSAXorFmzJEkffvihJOncc8+NXEQAEEPIjwDQdyGNAX3nnXd0wQUXBB8TJ07UX/7yl0jHBgCOR34EgL4LqQD9h3/4By1fvlwtLS1qaWnR8uXLdffdd0c6NgBwPPIjAPRdSAVoS0uL5s6dG9y+/fbb1dLSErGgACBWkB8BoO9CKkBdLpfefvvt4Pbbb78tl8sVsaAAIFaQHwGg70KahLRq1SpdfvnlmjRpkizL0p///Gdt3Lgx0rEBgOORHwGg70IqQK+++mq98847+uMf/yhjjC699FJlZmZGOjYAcDw78mNNTY0WLFig5uZmnX322XryySd1wQUXdGrz0ksv6Tvf+Y5aW1tlWZauu+46rV69WklJIV3oAgBbhZSZqqqqNHDgQF1//fW64YYbdNZZZ6m6ujrSsQGA49mRH5cuXaolS5Zoz549WrZsmYqKirq0SU9P19NPP623335br732ml555RWVl5eHNQ4ACJeQCtClS5cqNTU1uJ2amqq77rorYkEBQKyIdH5sampSdXW15s2bJ0maM2eOvF6vamtrO7W76KKLNHbsWEnSwIEDNXnyZNXV1YUtjkRmjFFV3SH9otqrqrpDMsZEOyQg5oV0CT4QCHQaVJ+cnKz29vaIBQUAsSLS+dHr9Wr06NFKTj6Zri3LUnZ2thoaGpSbm9vtz3z00UfavHmznn/++bDFkagaW9o0v6xS3kNtSnElyecPKCsjVeULp8qdntr7CwDoVkg9oGeddZZqamqC23v27FFKSkrEggKAWBGN/NhTD9ynn36qG264Qffff7/y8/O7bVNSUiK32x18tLa2RirUmGaM0fyyStUfbJPPb9R2wi+f36j+YJsWlFXSEwqcgZB6QB944AHNnDlT1157rYwx+t3vfqcNGzZEOjYAcLxI58esrCzt27dP7e3tSk5OljFGXq9X2dnZXdoePnxY11xzjQoLC1VcXHza1ywuLu603+12hy3eeFJd36LGQ0flD3QuNP0Bo4ZDbaqub1GBJyNK0QGxLaQe0Ouuu047d+7U1KlTdckll+jll1/WNddcE+nYAMDxIp0fR44cqSlTpgSXdtqyZYvcbneXy++tra265pprdPXVV2vFihVhO34iq2s+omSX1e2+FFeS6pqP2BwRED9C6gGVTn5Dnjx5siRp9OjREQsIAGJNpPNjaWmpioqKtGrVKg0dOjTYw7po0SIVFhaqsLBQjzzyiCorK3XkyBH98pe/lCTdeuutWr58edjjSRSezMHy+QPd7vP5A/JkDrY5IiB+hFSAvvLKK5ozZ45GjRolY4yampq0ZcsWXXrppZGODwAczY78OGHCBO3atavL8+vXrw/+e/ny5RSbYZafk66sjFTVH2zrdBnelWQpOyNV+TnpUYwOiG0hXYIvLi7W5s2b9cYbb+jNN9/U5s2b9U//9E+Rjg0AHI/8GL8sy1L5wqnKGZ762baU4rLkGZ6q8jsvkWV1f3keQO9C6gE9evSoZsyYEdyePn26jh07FrGgACBWkB/jmzs9VduKZ+my/9gunz+gNXOnKD8nneITOEMh9YCmpaVp69atwe1t27Zp8GDGvgAA+TH+WZalQSkuDR2YogJPBsUnEAYh9YA+8sgjmjNnjlwulyzLkt/v15YtWyIdGwA4HvkRAPoupAI0Pz9ftbW1evfdd2WM0fnnn89C9AAg8iMA9EdIl+BvuukmpaSk6Itf/KImTpyolJQU3XTTTZGODQAcj/wIAH0XUg9oQ0NDl+fee++9sAcDoHvGGFXXt6iu+Yg8mYOZBOEg5EcA6LseC9DS0lKtXbtWe/bs0ZQpU4LPf/LJJ7rgggsiHhwAqbGlTfPLKuU91KYUV5J8/oCyMlJVvnCq3Omp0Q4vYZEfAaD/eixAr7nmGk2YMEF33323Hn744eDzQ4cO1aRJkyIeHJDojDGaX1YZXAjb5/dLkuoPtmlBWaW2Fs+iJzRKyI8A0H89FqA5OTnKycnRO++8Y1c8AE5RXd+ixkNHO92FRZL8AaOGQ22qrm9RgScjStElNvIjAPRfSGNAZ8+e3W0vy0svvRT2gAD8VV3zESW7LJ3wd92X4kpSXfMRCtAoIz8CQN+FVID+8z//c/Dfx44d01NPPaW8vLyIBQXgJE/mYPn8gW73+fwBeTJZ8DzayI8A0HchFaDXXXddp+0bb7yxy3MAwi8/J11ZGanBMaAdXEmWsjNSlZ+THsXoIJEf0RkrVgChCakA/Ty/36/3338/3LEA+BzLslS+cKrml1Xq/QNHZFlS8mfFZ/mdl/CHzYHIj4mLFSuA0IVUgP7d3/1d8A+d3+/XW2+9pa985SsRDQzASe70VG0rnqXL/mO7fP6A1sydQq+Kg5AfIbFiBdBXIRWgp97VIzk5Wd/+9rfV3NwcsaAAdGZZlgaluDQoxcWkI4chP0JixQqgr0IqQBcsWCBJevfdd1VWVqb77rtPbreb280BSHjkR0isWAH0Va8FaFtbm/77v/9bTzzxhN577z0dPXpUu3bt0vnnnx/WQGpqarRgwQI1Nzfr7LPP1pNPPtnlbiIvvfSSvvOd76i1tVWWZem6667T6tWrlZQU0i3tASCs7MqPcD5WrAD6psfKbcmSJcrKytKvfvUr3X///WpoaNCwYcMiklyXLl2qJUuWaM+ePVq2bJmKioq6tElPT9fTTz+tt99+W6+99ppeeeUVlZeXhz0WAOiNnfkRztexYoUrqfM4T1asALrXYwG6adMmTZo0SUuXLtX111+v5OTkiAyibmpqUnV1tebNmydJmjNnjrxer2prazu1u+iiizR27FhJ0sCBAzV58mTV1dWFPR4kDmOMquoO6RfVXlXVHZIxpvcfAmRffkRs6FixImd46mfbUorLkmc4K1YA3enxEvy+ffv09NNPa+XKlVqyZInmz58vn88X9iC8Xq9Gjx6t5OST4ViWpezsbDU0NCg3N7fbn/noo4+0efNmPf/882GPJ1ribf04p58PS6bgTNiVHxE7WLECCF2PPaBpaWlatGiRdu3apd/85jc6evSoTpw4oenTp+uxxx6LaGA99UR9+umnuuGGG3T//fcrPz+/y/6SkhK53e7go7W1NZKhhkVjS5uuLNmhuete1QMVuzV33au6smSHGlvaoh1avzj9fE5dMsXnN2o74ZfPb4JLptATit5EMz/CuTpWrBg6MEUFnowzLj65SoN4FfLsnQsuuEA/+tGP9MEHH6i4uDisPY9ZWVnat2+f2tvbJZ38hfN6vcrOzu7S9vDhw7rmmmtUWFio4uLibl+vuLhYjY2NwUdaWlrYYo2EeCuGYuF8QlkyBQhVJPMjEpfTv8gDZ6LP08eTk5N1yy236Ne//nXYghg5cqSmTJmijRs3SpK2bNkit9vd5fJ7a2urrrnmGl199dVasWJF2I4fbfFWDMXC+XQsmdKdjiVTgL6KRH5EYoqFL/LAmXDM+kWlpaUqLS1VXl6eVq9erQ0bNkiSFi1apIqKCknSI488osrKSv3yl7/U5MmTNXnyZH3/+9+PZthhEW/FUCycD0umAHCyWPgiD5yJft0LPhImTJigXbt2dXl+/fr1wX8vX75cy5cvtzMsW8RbMRQL59OxZErHbfM6sGQKADv0NkmThe0R7xxTgCayeCuGYuF8OpZMmV9WqfcPHJFlScmfxceSKQAiKZQVOPrzRd7pK48Ap3LMJfhEFm/rx8XK+XQsmeJOH6SRQwboqcXTtLV4lsYMGxTt0ADEqVDHdvZ1YXsmLCHWUIA6RLwVQ7FyPuFeMgUAehLq2M6+fJHv64QllnaCE3AJ3kE6iqFBKa64GNsTb+cDAGeqL2M7Q13YPpSituM1uQEHnIIeUAAAbNLXsZ2hXKUJdeURlnaCk1CAAgBgk76O7QxFqEUtSzvBSShAAQCwSSQmaYZa1MbCGs1IHBSgAADYKNyTNEMtavt6+Z/JSogkJiEBAGCzcE/SDGXCUl/WaGayEiKNHlAAAOJAbxOWQu0pZbIS7EABCgAOV1NTo+nTpysvL08FBQXavXt3t+2eeOIJjR8/XuPGjdPixYvl8/lsjhROF8rlfyYrwQ4UoDGIcTlAYlm6dKmWLFmiPXv2aNmyZSoqKurSZu/evVqxYoV27typ2tpa7d+/X48//rj9wcLxeuspZbIS7EABGmO43Rrsxhee6GpqalJ1dbXmzZsnSZozZ468Xq9qa2s7tdu8ebMKCws1atQoWZalu+66S5s2bYpGyIhx/bkPPdBXTEKKIaeOy/EHjHz+k7fS6BiXs7V4FreSRFgxESH6vF6vRo8ereTkk+nasixlZ2eroaFBubm5wXYNDQ3KyckJbns8HjU0NIQ9nsW/fUyZhw/ovZd6/vz/5eDJL8W9tetL23C3i5XXtPvYGUYqPdimE90UoWe5kpTxZqre6+ZPzYcfH5MknTtsYI/HDne7WHnNWDr2J8NG6vrnn+q17ZmgBzSG9GdcDr1X6C8mIjhXKO99T21KSkrkdruDj9bW1pCPnXqWSymu3v90pLiSQmrXl7bhbhcrr2n7sS1pTPognXXKfksni88x6YNObnTD5w+ctuc0ku1i5TVj6dht3d0rNszoAY0hfbmHsETvFc5MX+4vjcjJysrSvn371N7eruTkZBlj5PV6lZ2d3alddna23nvvveB2fX19lzYdiouLVVxcHNx2u90hxxNqr8i4kF8x9LbhbhcrrxmtY+cZo+r6FtU1H5Enc3C396E/1V0lOyRJLxbP6vF1w90uVl4z1o791V5bnhl6QGNIX8bl0HuFM8VEBGcYOXKkpkyZoo0bN0qStmzZIrfb3enyu3RybGhFRYU++ugjGWO0du1a3X777dEIGXHCsiwVeDJ0a37Wae9DD/QXBWgM6cs9hFlGA2eKiQjOUVpaqtLSUuXl5Wn16tXasGGDJGnRokWqqKiQJI0dO1YrV67UjBkzlJubqxEjRmjp0qXRDBsATotL8DGkYxHh+WWVev/AEVmWlPxZ8fn5ewj39XI98Hl9uWsKImvChAnatWtXl+fXr1/faXvx4sVavHixXWEBQcYYHfX55fMHVFV3qNfL9QA9oDEm1HsI03uFMxXqXVMAJLaO5QEbW46q6fBxlgdESChAY1BviwhLfbtcH4+Y/R8eoX7hAZCYTp1vcHJbzDdASLgEH6f6crk+3jD7P7w6vvAMSnExbANAJ6yWgf6iBzSOJWLvFbP/AcA+rJaB/qIAjXOhXK6PJ8z+BwD7MN8A/UUBirjCt3EAsE+izzdA/1GAIq7wbRwA7HPqahkpLuuzW7WyWgZ6xyQk9Inp463Z7MbalQBgr475Bk7+2wDnoQBFyGJhdnkiz/4HgGjpuG0nM94RKi7BIySxNLs8EWf/AwAQSyhAEZJYm12eaLP/AQCIJRSgCAmzy8OPuzWFD+8lAMQWxoAiJMwuD69YGE8bK3gvASD20AOKkLDWW/jE0nhap+O9BIDYRAGKkJy61tvJbbHWWz/F2nhaJ+O9BIDYRAGKkDG7PDwYTxs+fX0vGSsKAM7AGFD0Scfs8kEpLtZ76yfG04ZPX95LxooCgHPQAwrYjPG04RPqe8lYUQBwFgpQwGaMpw2fUN9LxooCgLNQgAJRwHja8AnlvWTcLQA4CwUoECXcrSl8ensvGXcLOIMxRkd9fn16zMdEwARHAQog7jHuFoi+xpY2XVmyQ40tR9V0+Ljmrnv1s+22aIeGKKAABRD3GHcLRNepEwFPbouJgAmOAhRAQmDcLRA9TATE5zmmAK2pqdH06dOVl5engoIC7d69u9t2TzzxhMaPH69x48Zp8eLF8vl8NkcKIFYx7haIDiYC4vMcU4AuXbpUS5Ys0Z49e7Rs2TIVFRV1abN3716tWLFCO3fuVG1trfbv36/HH3/c/mABAEDImAiIz3PEnZCamppUXV2t3/3ud5KkOXPm6Bvf+IZqa2uVm5sbbLd582YVFhZq1KhRkqS77rpLq1at0j333BPWeBb/9jFlHj6g917q+e4oH358TJJ07rCBvb5mqG3/5bPxMb0dO9R2fWkb7hgjcey+vOfhfi+j+XlHIs5IvJeROHYkfneah4yQimf1emwA4dExEbD+YFuny/BMBExcjugB9Xq9Gj16tJKTT9bDlmUpOztbDQ0Nndo1NDQoJycnuO3xeLq0kaSSkhK53e7go7W1tU/xpJ7lUoqr97fG5w+c9htdf9umuJJCOnao7frSNtwxRuLYfXnPw/1eRvPz7kvbaL6XkTh2JH53Us9yhXRsAOFx6kTAFJf12d9ZJgImMkf0gHYnlBlxp2tTXFys4uLi4Lbb7e7Tsa9//qmQ2t1VskOS9GIIPSmhth0X0pFDb9eXtuGOMRLH7st7Hu73Mpqfd1/aRvO9jMSxo/m7AyB8OiYCVte3qK75iDyZg5Wfk07xmaAcUYBmZWVp3759am9vV3Jysowx8nq9ys7O7tQuOztb7733XnC7vr6+SxsAAOBMlmWpwJOhAk9GtENBlDniEvzIkSM1ZcoUbdy4UZK0ZcsWud3uTuM/pZNjQysqKvTRRx/JGKO1a9fq9ttvj0bIAAAA6CdHFKCSVFpaqtLSUuXl5Wn16tXasGGDJGnRokWqqKiQJI0dO1YrV67UjBkzlJubqxEjRmjp0qXRDBsAAAB95IhL8JI0YcIE7dq1q8vz69ev77S9ePFiLV682K6wAAAAEGaO6QEFAABAYqAABQAAgK0oQAEAAGArClAAcKhAIKBvfvObGjdunHJzc7VmzZpu2x07dkw33XST8vLydOGFF+qqq65SbW2tzdECQOgoQAHAoTZu3Ki3335be/bsUWVlpX74wx9q9+7d3bZdsmSJ3n33Xf3pT3/SjTfeqEWLFtkcLQCEjgIUUWeM0VGfX58e86mq7lBId8ECEsEzzzyjxYsXy+VyKSMjQ7fddps2bdrUpd3AgQP1la98JXhHmWnTpqmurs7maAEgdBSgiKrGljZdWbJDjS1H1XT4uOaue/Wz7bZohwZEXUNDg3JycoLbHo9HDQ0Nvf7cj3/8Y914442RDA0AzggFKKLGGKP5ZZWqP9j22bbk8xvVH2zTgrJKekJjDD3ZfXfppZcqMzOz24fX6+3SPpT3dNWqVaqpqdEPfvCD07YpKSmR2+0OPlpbW8/oPIBwI5/EP8csRI/EU13fosZDR+UPdE4s/oBRw6E2Vde3cL/gGNHY0qb5ZZVqbDkqy5LmrntVWRmpKl84Ve701GiH51jd3XzjVNnZ2aqvr9ell14qSaqvr1d2dvZp2//nf/6nnn32WW3dulWpqad/34uLi1VcXBzcdrvdfYwciBzySWKgBxRRU9d8RMkuq9t9Ka4k1TUfsTki9Ac92ZFz6623at26dfL7/WppadEzzzyj2267rdu2JSUl2rRpk1588UUNGzbM5kiB8CCfJA4KUESNJ3OwfP5At/t8/oA8mYNtjgj9EUpPNvrnjjvu0Pnnn6/x48crPz9fxcXFmjhxoiSpoqIiONO9sbFR9913nz7++GPNnj1bkydP1iWXXBLN0IF+IZ8kDi7BI2ryc9KVlZGq+oNtnZKNK8lSdkaq8nPSoxgdQtXRk33C33VfR082Qyn6x+Vy6Sc/+Um3+woLC1VYWCjp5CV0eoYQD8gniYMeUESNZVkqXzhVOcNTleKylHqWSykuS57hqSq/85LgkjJwNnqyAYQL+SRx0AOKqHKnp2pb8SxV17eorvmIPJmDlZ+TTvEZQyLZk90xE9bnD6iq7hD/N4A4x5WxxEEPKKLOsiwVeDJ0a36WCjwZFBgxJlI92awRCyQerowlDnpAAZyxcPdk9zYTdmvxLP4QAXGKK2OJgQIUQFh09GSHY4IAa8QCiS2c+QTOxCV4AI7DGrEAEN8oQBEzonlrNm4LZy9mwgJAfKMARUyI5oQUJsPYr2MmrCupcy8oM2EBID5QgMLxonlrNm4LFx3MhAWA+MYkJDheNCekMBkmepgJCwD2snPtZXpA4XjRnJDCZJjoYo1YALCH3cPNKEDheNGckMJkGABAvIvGcDMKUDheNCekMBkGABDvQhluFm4UoDZgCZ8zE80JKUyGAQDEu2gMN2MSUoQ1trRpflmlGluOyrKkueteVVZGqsoXTpU7PTXa4cWMaE5IYTIMADiXnRNn4lU0hpvRAxpBLOETXtGckMJkGABwHtZp7l0oV2GjMdyMAjSCojGmAgCAREAnT+9CLdCjMdyMS/AR1DGm4oS/676OMRWsIYlEw+UyAOHAOs09661A31o8q1PutXu4GT2gEcQSPkBnXC4DEC6s09yz/lyFtXO4GQVoBLGED/BXXC4DEE508vTM6QU6BWgEsYQP8FeMiQYQTnTy9MzpBTpjQCOMJXyAkxgTDSCcOjp55pdVynuoTSmuJPn8AWVn0Mkj/bVArz/Y1umLv1MKdApQG3SMqeCPKxKZ07+NA4g9dPKcntMLdApQALZw+rdxALEpUTt5QllRxMkFOmNAIYnbhSLyGBMNAOHRlxVFnHojFXpAwe1CE1C01uJ08rdxAIgFfV3f06noAU1wLI0Tfk7vTY72WpxO/TYOALEgXlYUoQBNcPHyH9kpol3c9YYvHAAQ25y+vmeoKEATXLz8R3aCWCju+MIBALEtXlYUoQBNcPHyH9kJYqG44wsHAMS2eFmAP+oFaCAQ0De/+U2NGzdOubm5WrNmTbftjh07pptuukl5eXm68MILddVVV6m2ttbmaONPvPxHdoJYKO74wgEAsS1eVhSJ+iz4jRs36u2339aePXv0ySef6KKLLtLs2bN1wQUXdGm7ZMkSXXvttbIsS2vWrNGiRYv0+9//3v6g44jTF6qNJbFQ3LEWJwDEvnhYUSTqPaDPPPOMFi9eLJfLpYyMDN12223atGlTl3YDBw7UV77yleCbO23aNNXV1dkcbXzq+I/81OJpWll4gZ5aPE1bi2dpzLBB0Q4tpsRCb3K8fHMGgL5y+golfRXrK4pEvQe0oaFBOTk5wW2Px6NXX32115/78Y9/rBtvvDGSoSWURL2TRDjFSm9yPHxzBoC+YL1r54l4AXrppZeqpqam231vvPFGl+dC+UayatUq1dTUaNu2bd3uLykpUUlJSXC7tbU1xGiBMxMrxR1fOAAkinhZuD3eRLwA3bVrV4/7s7OzVV9fr0svvVSSVF9fr+zs7NO2/8///E89++yz2rp1q1JTu//WUlxcrOLi4uC22+3uR+RA/1DcAYBzhLJCCfnaflEfA3rrrbdq3bp18vv9amlp0TPPPKPbbrut27YlJSXatGmTXnzxRQ0bNszmSAHAXqGuEnKqDRs2yLIsPffcczZECDhfLKxQkoiiPgb0jjvuUFVVlcaPHy/LslRcXKyJEydKkioqKlRRUaH169ersbFR9913n8aOHavZs2dLkgYMGKA//vGP0QwfACKmL6uESCevIK1bt07Tpk2zOVLAuWJhhZJEFPUC1OVy6Sc/+Um3+woLC1VYWCjp5GX0WJ+xBgB9cbpVQh566KEubQOBgO688049+uijuu+++6IQLeBMsbb8XMdsfZ8/oKq6Q46cRxAOUb8EDwDoXnerhDQ0NHTbtqSkRDNmzNDFF19sV3hATIil5ecaW9p0ZckONbYcVdPh45q77tXPttuiHVrYRb0HFAASVbhWCdm9e7c2b96snTt3hnRcVgpBoomFFUoSbbY+BSgAREm4Vgn5wx/+oPr6eo0fP16S9NFHH2nJkiXat2+f7r777i7tWSkEicjpK5Qk2mx9LsEDgEOFukrI3XffrX379qmurk51dXWaNm2aHn/88W6LTwDhE867KyXabH0KUABwqDvuuEPnn3++xo8fr/z8/C6rhCxatCjKEQKJK9zjNRNttj6X4AHAoUJdJeTzfv/730cwKiD+9TYTPRLjNWNttv6ZogcUAADgM6H0bIYyXrOvYmm2fjjQAwoAAKDQezY7xmue8Hd9jY7xmv2ZMBQLs/XDhR5QII6Ec0A8ACSaUHs2Izles2O2/q35WSrwZMRl8SlRgAJxI5EWMAaASAh1JnrHeE1XUue28TpeMxIoQPuJniY4SW+Xjfj/CQC9C7Vns6/jNakZumIMaD80trRpflmlGluOyrKkueteVVZGqsoXTpU7PTXa4SEBJdoCxgAQCX2ZiR7qeE1qhu7RA9pH9DTBiRJtAWMAiIS+9mz2Nl6TmuH06AHtI3qa4ESJtoAxAERKOGeiUzOcHj2gfURPE5yIAfEAED7hmolOzXB6FKB9RE8TnCjRFjAGgFhAzXB6XILvo0S7VRZiRyItYAwAsYCa4fToAe0jeprgZImygDEAxAJqhtOjB7Qf6GkCAAChoGboHgVoP3X0NCXq7DUAABAaaoauuAQPAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALAVBSgAAABsRQEKAAAAW1GAAgAAwFYUoAAAALCVZYwx0Q4i0nsJN4AAAAbkSURBVAYMGKARI0b06WdaW1uVlpYWoYjsF0/nE0/nInE+TtdxPgcOHNDx48ejHU5E9DVHxutnHC/i6Xzi6Vyk+D2f/uTHhChA+8PtdquxsTHaYYRNPJ1PPJ2LxPk4XbydTzjE23vC+ThXPJ2LxPmcikvwAAAAsBUFKAAAAGzlevDBBx+MdhBOdemll0Y7hLCKp/OJp3OROB+ni7fzCYd4e084H+eKp3OROJ8OjAEFAACArbgEDwAAAFtRgAIAAMBWFKCfU1NTo+nTpysvL08FBQXavXt3tEM6Ix6PRxMmTNDkyZM1efJkPfPMM9EOqU/uvfdeeTweWZalN998M/h8LH5OpzuXWP2Mjh07pptuukl5eXm68MILddVVV6m2tlaS1NTUpGuuuUbjx4/XF7/4Rf3hD3+IcrS96+l8rrjiCp133nnBz+jhhx+OcrTREYu/dz2J1d+9DvGUH6X4ypHkxxAYdDJ79myzYcMGY4wxv/jFL0x+fn50AzpDOTk55o033oh2GP22Y8cO4/V6u5xHLH5OpzuXWP2Mjh49al544QUTCASMMcY8+uijZtasWcYYY77+9a+bBx54wBhjTGVlpRkzZow5ceJElCINTU/nM2vWLPPLX/4yitE5Qyz+3vUkVn/3OsRTfjQmvnIk+bF3FKCn2L9/vxkyZIjx+XzGGGMCgYA555xzTE1NTZQj679Y/MXtzqnnEeufUzwk1+5UVVWZnJwcY4wxgwcPNvv27QvuKygoMC+++GKUIuufU8+HAjT2f++6Ey+/e/GUH42JzxxJfuyKS/Cn8Hq9Gj16tJKTkyVJlmUpOztbDQ0NUY7szMyfP18TJ07UnXfeqQMHDkQ7nDMWj59TPHxGP/7xj3XjjTfq4MGD8vl8GjVqVHCfx+OJuc+n43w6fOc739HEiRN122236f33349iZNERj793Unz87p2Kz8mZyI9dUYD2wsT4KlV/+MMf9NZbb+n1119XZmamFixYEO2QIiKWP6d4+IxWrVqlmpoa/eAHP5B08o/eqWLt8/n8+fzsZz/TX/7yF7311lu67LLLdP3110c5QmeItc/18+Lhdy8UfE7RRX48jbD1x8aBeLh00ZMPP/zQpKWlRTuMfomnS0w9XU6Kxc/ohz/8obn44otNS0tL8LnU1NSYvcTU3fl83oABA0xzc7ONUUVfrP/e9SYWf/c6xFN+NCa+ciT58fToAT3FyJEjNWXKFG3cuFGStGXLFrndbuXm5kY5sv45cuSIPv744+D2pk2bdNFFF0UxovCIp88p1j+jkpISbdq0SS+++KKGDRsWfP7WW2/V2rVrJUlVVVX64IMPNGvWrGiFGbLuzqe9vV379+8PttmyZYvOOeccDR8+PFphRkU8/d5Jsf+7dzp8Ts5BfuwZd0L6nHfffVdFRUU6ePCghg4dqg0bNmjixInRDqtf3n//fc2ZM0d+v1/GGI0dO1aPPPKIPB5PtEML2dKlS/XCCy/oo48+0vDhwzVkyBDV1tbG5OfU3bn87ne/i9nPqLGxUVlZWRo7dqyGDBkiSRowYID++Mc/av/+/brjjju0d+9enXXWWVqzZo1mz54d5Yh7drrzeemllzRr1iwdP35cSUlJyszMVElJiS688MIoR2y/WPy9Ox3yo/PEU44kP/aeHylAAQAAYCsuwQMAAMBWFKAAAAD4/+3bIa7qUACE4SEpkJAQWAEhwWAImKYaDA6BQVWAZQtYbBEsAYNhDVADooIGFCFdAgYEUNPnEE/fd87N4/9ce8yoyeQ0NYoBCgAAAKMYoAAAADCKAQoAAACjGKD479XrdcVxbDsGAPw69CNsYYACAADAKAYovlIQBHJdV51OR67r6nA4fM72+706nY5arZYmk4na7bZ2u529sABgEP0IExig+Eq+7yuKIsVxrOVyqfF4LElK01Sj0UhBEOh8Psv3fZ1OJ8tpAcAc+hEmOLYDADYcj0fN53Pdbjc5jqPL5aLn86nr9SrHcdTr9SRJ3W5XjUbDcloAMId+hAkMUHydNE01HA613W7luq7u97sqlYre77eyLFMul7MdEQCsoB9hCp/g8XVer5fSNFWtVpMkLZfLz1mz2VSapgrDUJIUhqGSJLGSEwBMox9hCjeg+Ar9fl/5fP7zPJvN5HmearWaBoPB532xWNR6vdZ0OlWhUJDneWo2m6pWqzZiA8A/Rz/ChlyWZZntEMBv8ng8VC6XJUlRFGkwGChJEpVKJcvJAMAu+hE/hRtQ4C+bzUaLxUJZlslxHK1WK8oVAEQ/4udwAwoAAACj+AkJAAAARjFAAQAAYBQDFAAAAEYxQAEAAGAUAxQAAABGMUABAABgFAMUAAAARv0ByrubqZ/w4XUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f05a5195128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lAutoIndependent = Autocorrelation(lXIndependent)\n",
    "lAutoDependent = Autocorrelation(lXDependent)\n",
    "\n",
    "plt.figure(num=None, figsize=(15, 6), dpi=80, facecolor='w', edgecolor='k')\n",
    "# independent\n",
    "ax1 = plt.subplot(131)\n",
    "ax1.stem(lAutoIndependent)\n",
    "ax1.set_xlabel('Lag')\n",
    "ax1.set_ylabel('Autocorrelation')\n",
    "ax1.set_title('independent')\n",
    "\n",
    "# dependent (AR1)\n",
    "ax2 = plt.subplot(132)\n",
    "ax2.stem(lAutoDependent)\n",
    "ax2.set_xlabel('Lag')\n",
    "ax2.set_ylabel('Autocorrelation')\n",
    "ax2.set_title('dependent (AR1)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate ESS using,\n",
    "\n",
    "## $$ESS = \\frac{n}{1 + 2 \\sum_{t = 1}^{T} \\hat{\\rho}_t},$$ ##\n",
    "\n",
    "## where $\\hat{\\rho}_T$ is the last positive sample autocorrelation.## \n",
    "\n",
    "Creating function to find $T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AutocorrelateNegative(lAutocorrelation):\n",
    "    T = 1\n",
    "    for a in lAutocorrelation:\n",
    "        if a < 0:\n",
    "            return T - 1\n",
    "        T += 1\n",
    "    return -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And function to calculate ESS for a single chain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ESS_single(lX):\n",
    "    lRho = Autocorrelation(lX)\n",
    "    T = AutocorrelateNegative(lRho)\n",
    "    n = len(lX)\n",
    "    ESS = n / (1 + 2 * np.sum(lRho[0:(T-1)]))\n",
    "    return ESS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare ESS for independent and dependent (AR1) series - the equivalent number of independent samples for the dependent sampler is tiny relative to its sample size. We would need roughly 10X more samples from the dependent sampler to explore the posterior to the same degree as the same degree as the independent sampler. This sampler is very inefficient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "independent ESS = 50.0\n",
      "dependent ESS = 5.1472886755\n"
     ]
    }
   ],
   "source": [
    "ESSIndependent = ESS_single(lXIndependent)\n",
    "ESSDependent = ESS_single(lXDependent)\n",
    "\n",
    "print('independent ESS = ' + str(ESSIndependent))\n",
    "\n",
    "print('dependent ESS = ' + str(ESSDependent))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare ESS estimator for AR1 with true value (since we know $\\rho$ for the AR1).\n",
    "\n",
    "For large values of $\\rho$, I think the above estimator will underestimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "est. dependent (rho = 0.99) ESS = 5.1472886755\n",
      "true dependent (rho = 0.99) ESS = 0.2512562814070354\n"
     ]
    }
   ],
   "source": [
    "# geometric progression\n",
    "def ESSTrueAR1(nSamples, rho):\n",
    "    ESS = nSamples / (1 + ((2 * rho) / (1 - rho)))\n",
    "    return ESS\n",
    "\n",
    "print('est. dependent (rho = 0.99) ESS = ' + str(ESSDependent))\n",
    "print('true dependent (rho = 0.99) ESS = ' + str(ESSTrueAR1(len(lXDependent),rho)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For smaller values of $\\rho$ we get a better estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "est. dependent (rho = 0.8) ESS = 6.33865004864\n",
      "true dependent (rho = 0.8) ESS = 5.5555555555555545\n"
     ]
    }
   ],
   "source": [
    "rho = 0.8\n",
    "lXLowerRho = AR1(rho,sigma,T)\n",
    "\n",
    "print('est. dependent (rho = 0.8) ESS = ' + str(ESS_single(lXLowerRho)))\n",
    "print('true dependent (rho = 0.8) ESS = ' + str(ESSTrueAR1(len(lXLowerRho),rho)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use Pints function to do the same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "est. dependent (rho = 0.8) ESS = 6.33865004864\n"
     ]
    }
   ],
   "source": [
    "import pints._diagnostics as diagnostics\n",
    "\n",
    "print('est. dependent (rho = 0.8) ESS = ' + str(diagnostics.ess_single_param(lXLowerRho)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}