gptide toeplitz example

test_toeplitz.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f063dffe",
   "metadata": {},
   "source": [
    "# 1D kernel basics\n",
    "\n",
    "This example will cover:\n",
    "\n",
    " - Initialising the GPtide class with a kernel and some 1D data\n",
    " - Sampling from the prior\n",
    " - Making a prediction at new points\n",
    " - Sampling from the conditional distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ad0745e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gptide import cov\n",
    "from gptide.gpscipy import GPtideScipy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef043045",
   "metadata": {},
   "source": [
    "We use an exponential-quadratic  kernel with a length scale $\\ell$ = 100 and variance $\\eta$ = $1.5^2$. The noise ($\\sigma$) is 0.5. The total length of the domain is 2500 and we sample 100 data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "32bb563a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change this to really test the Toeplitz\n",
    "N = 15000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6a64bbbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "####\n",
    "# These are our kernel input parameters\n",
    "noise = 0.5\n",
    "η = 1.5\n",
    "ℓ = 100\n",
    "covfunc = cov.expquad_1d\n",
    "\n",
    "###\n",
    "# Domain size parameters\n",
    "dx = 25.\n",
    "covparams = (η, ℓ)\n",
    "\n",
    "# Input data points\n",
    "xd = np.arange(0,dx*N,dx)[:,None]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "806fb022",
   "metadata": {},
   "source": [
    "## Initialise the GPtide object and sample from the prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "117e346b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'x')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAACBR0lEQVR4nO29eXwc1Z0v+j3VLcmWLEttrbZlSZZtjJENxPIKhC1AIANhCyGQN5kMQ4B5mTuTN5m52f24JJOXeTNzJ3fmcq9ZJjfvztiOAbOFgQQMBkPiVR4bS7blRbZkWbvcktrau+u8P6pO9alTp6qrW93qllXfz8eg7q6uOl11zm/9/n6HUErhwYMHDx5mHpR0D8CDBw8ePKQHngLw4MGDhxkKTwF48ODBwwyFpwA8ePDgYYbCUwAePHjwMEPhT/cA4kFxcTGtrq5O9zA8ePDgYVqhvr6+l1JaIr4/rRRAdXU1Dh48mO5hePDgwcO0AiGkRfa+FwLy4MGDhxkKTwF48ODBwwxFWhUAIaSQEPIKIeQEIeQ4IWRjOsfjwYMHDzMJ6c4B/DcAv6GUfokQkg0gN83j8eDBg4cZg7QpAEJIAYAbAXwdACil4wDG0zUeDx48eJhpSGcIaDGAHgD/ixDyH4SQFwkheeJBhJAnCCEHCSEHe3p6pn6UHjx48HCZIp0KwA9gNYD/SSn9DIAhAN8VD6KUPk8pXUMpXVNSYqGxephmqG8J4tldp1HfEkz3UDx4mPFIZw6gDUAbpXSf/voVSBSAh8sH9S1BfPXFvRgPq8j2K9jy+AbUVQXSPSwPHmYs0uYBUEo7AZwnhCzX3/ocgGPpGo+H1GNvcx/GwypUCkyEVext7kv3kDx4mNFINwvoPwHYojOAmgH8cZrH4yGF2FBThGy/gomwiiy/gg01RekekgeXqG8JYm9zHzbUFHle22WEtCoASulhAGvSOQYPU4e6qgC2PL7BEyTTDF7o7vJFuj0ADzMMdVUBT3hMM+xt7sPYhAoKYHxCC915z/DygNcKwoMHD44I5GaD7Ryu6q89XB7wFIALeNRFDzMZweFxKET7WyHaaw+XB7wQUAx48c/kwUskTk9sqCmC36cl7/0+L3l/OcFTADEgoy56wit+eIp0eiOiajmAiKpaPvMU+/SFFwKKAUZd9BF41MVJwKsBmL547qMziOhyP6JqrxmYYv+Hd5vw1Rf3emHSaQbPA4gBj7qYHHg1ALGRqZZ01+Co7WvPQ57e8BSADqfFx1MXM3WRZjqYIt1xqA0k3YPJQGRyiOzhtZU40nbU9JrBU+zTGzNOAcgEuNvFl8mLdLrglYPnMRGhePngeWx7YqN3/3RksiW9vDwffgUIq4Bf0V4zeB6yHNPFUJxRCsBOgLtdfPxx4xMqfr7zJL512xUZ/YAzCTsOtWE8ojHKxyMUOw61efdORyZb0nub+6DqhQCUwrI+Mqm4LxME73QyFGeUArAT9G4XHztufEKFCuB3p3tx4NzFtD3gTJjs8UAM/XihoCgy2ZLOZOXEI1MEbyZ7cyJmlALYUFMEv0IwEaHwKcSYyG4XHzvu5ztP4nene9P6gDNlsseDB1ZX4OX6NkOQPLC6It1DyihkkiXNI5OVE49MEbx2ciYTMaMUAACAEABU/38UbhdfXVUA37rtChw4dzGtFlGmTPZ4UFcVwNP31OKdhg7ctXJ+xo/XQxSZqpx4ZJSnYiNnMg0zSgHsbe5DOKIXtEQSF5p1VQFsuju9giyjJrtL1LcE8cxbjRgPqzhw7iKWl+dnvFDxEBuZEorMFE8lWXJmKjCjFECyhGYmCLJMmezxYDp6LVOJTBGk8SDTQpGT9VSS8Qymk3E2oxRAsoSmiQ2URkE2HdxyHtNpYUw1Mk2QukWmrIVkIJnP4IHVFSD6/zP5fswoBQAkR2gGcrMNWpxKvfa4bjEdvZapQqq9o1R5F5fTWkjGMxCVSKYTHWacAkgGgsPj0FM8UOC1x3WL+pagVwlsg1R6R6n0LhraBxxfTyck4xlMtzCnpwASwIaaIuRkJTZRpmOcNxmobwnikef3GIVgL9e3Yds3pkeYYyqQSu8olUIpU2o7krGukvEMArnZUHQG0HQIc85IBTDZyZLoRJmucd5kYG9zHyYi1Hg9HayjqUaqcjqp9C4yobYjmetqMs+AkUMiqsb/33R3bcbP7xmnAJglOhGhyPIRSz8at8ohkYky3dzDZGJDTRF8ej8ZQAufTed48XRCKr2LVNZ2uF2LTutqKj1uNg4KgFI6LULDM04B2PWjYfHpV+rbEI7YWxKTmVAzmQVTVxXAw2srsWVfKwAvdyLD1n2thiB9dH2l6bNkeK2pEICpokTHY9XzO5b5uB3Lptrj5te3TyFo7x9BfUswo428GacAZDFLNlHGJlRj82sZpW2yE2qms2AeWF2BHYfaZqQCjIWt+1rx/de0lssfn+oFAEMJZHLoMFVebdznpdT8f/0coxOayzk2kXqPm61vZkhu29+KHYfaMup5iZhxO4I9sLoC2X4FBDBoWrzrxiCjtLEJpdLohIoXdVUBfPOWpRk7IVIJtkD+8o7lGb0o0oF3GjpsX2fybmqp2jEvnvPube5DWKX6lpXUuD+hkQnjGCq8ThXqqgJYWDgbE1zX4Ex6XiJmnAdQVxXAtm9YrfBsv2LyAGQhinRMqMsJ6WZApfL67NyB3GwEh8fjvsZdK+cblj97zZDJocPJECKcvhPPee1CL40dg6bjxNepQiA325AjKpxzXeleE2lXAIQQH4CDAC5QSu+eimuK8VDRdYtE5AstkQmV7gec6nE4xa3F66czjJHK64shRIUg7ms8ur4SrX1D+E1jJ+6sLTfdy2SGDlMxD+LNL7h9FrHOy/8WWejlztpy0/G18+e6/1GTQHB4HArRoggKsc91pXtNABmgAAD8BYDjAKbm6diATbYH9ZCQbIHUzp9rstJiTaipesCxFnWqxvGzt49j8+5mANa4tQi7lgF8cZjbsvlEhNirh9oMAZ1sBpYYQkwkHl7fEsQvfn8OE2EVv/j9OdxeW24xUiY73lTMg0Sen5v4fiJzemHhbIQj0fP2DUUFLwGQPztrUr/VLTbUFJkUgJ3HlgmswLQqAEJIBYA/APA3AP5yKq7pxvW0m4yhsbDp/VgTin/AYxNqSnbAcrOoUzHR6luCeP7jZtN77zR02CoAsWVAaGQC33/tKF46eB7hOIrDEhFi9S1BvHzwvCGgeaZIMiBuFKQkEA9/9VAbxnWO7HhYxauSuTJZ631vc5+hBMfjTIrabaXKF/dtP3ge25/YaFzLbpyx+uXHO6eZQSGGyu5aOT/utu3J8JCaOkMG3Tmsaq/d3oepjhik2wP4OYD/DCDf7gBCyBMAngCAykr7EIMbJCo82HcUgUIUi8fOHvB4REtQvVLfhgeT3BzKjXBPRQyZ3yaQgY9bixDd4Bc+OQtVT9wxuFFOiSgzliQENEvwS3XJfQZ8iCbRHACN8ToZ1nuisWkA9lupcsV94QjF5o/O4ONTPbHH6dAv380zlvUgkoXKlpfnuxaofCjPpxA8c+9Kx7CmHbYfaDW9djKM+PvQ1BkyKLVTFRJKmwIghNwNoJtSWk8IudnuOErp8wCeB4A1a9aI6yIuJCo82HdEgRer70ldVQAPrVmErftaU9Yb3I1wTwX9dENNEbJ9mnIjAJ68scZxsYjCJiLeTLgrDktEmYnfeTDJ1arJsNoeXF2Blw60Ghuvi2NMhheXaGz6wdUVtlup+gjA6QB0D47GLMpq7x9x7Jfv9Iz5c+ii00TYkOX33N4n3kMKqxSb3miIu65h675WHG0zywU7w4jfNyAcUfFOQ8eUh4TS6QFcD+CLhJAvAJgFYC4h5N8opf9Hqi7oVnjwC/pUV8gi+BmstosVyeK+27nge5v78PWN1WjsGLRUYorfSfpkIgRE73lyu5BwEyFTlgSaIGICxE1xWCLKjE/yu3lmbsHi3y/rYawsvzK5/kYOVnEyvDi35xCVDQWk36urCuAbn60x8kAAsLGmCE1dIcuxvFLxKwR+n2JLtrB7xuI5svzRcwRys/HsrtOT7uPvU4jhLaqUxh0m2/RGA1TuvduvKnMdFq2dP3fKdxpMmwKglH4PwPcAQPcA/iqVwh9wFh5GnH9kAi9+chYqpVBItHWBDLULCuK6ZiA32+AET5Y1AcDCPOErMe1CBsmwVutbgvj5zpNx7XokijSFAI+sq0TtggI8/WYDJiIUfp+7/VNFZeYmYcgL6mQ0omP3lxUaAfaxezfYcajNyIWEuQp1HpPtMe9Weco8ppULCqTtHvJnZ5ks8fzZWdJrmKzrCMVX1i/CwsLZceXieMUUUSkeXqedI5CbbRs6iWe+11UF8My9K/GjNxqgqhT+OHJFxpoQrMVblpfafofvKsyS1FNdKJruHMCUwy7JK1YCA9aQj4jtB1pduYjs80RjuHaFQE7ME7vvTDaOLKM8urFWRGX5xWsW4G/uX4X6luCk9k+VKbr3GjsNOuXtteVJEdSiIGH3V0SiMcre0Jjt62T2mHfjCdZVmbc8BWDb7kHWGVd2DTH/sHJBQUzasCgIZYqpriqAZ3edNs31HYfaDIMr3pj68vJ8+Ig2Rr6q2Gls/JrgQeDs0fL3hAI4cr4fG2qK8M1bljqOMZnICAVAKf0QwIfpur6sEliGwtl+9I9EmUBH2gbw1Rf3uppYvEAenVDxzK8bsekea7dAWUGRnevuxDyRfceOihnvvWLCnwC4fmkxvnXbFTHPw1s7APDWpx34w43VJvZLWKCH2llC/GeiovvZO8dx4FwQALB5dzMOtQYnLahlSobdX9FoGBaYYuz7saiS/YKgONY+YPSRmWq6oNjfxy4HALj3KtzmH9j1ZYaKXShPLARj/bwUQhDRiQZu57usqlj0NsVmkrz8IND/Q2MbRuKaeO9YF3af6pnSeoCMUADpAi9sGVtHhALNGvArwBVl+divCxcGtz1G+HgfoCmPh5/fg+1cN9L6liAeeWGvIbD4giLZJvROzBNxYQLA4fP9k969SbRa3HaAFOOrEVULc7zEMSZUaPRQJ8aL+Nmmu2tNjcDOXxw2Xfd09yWLJ5ftI3ElgmUC+Ju3LMWmu2vx/O4zONcXvebrh9sxPB7BkzctMRSZm30QxgQl1dY/ahgXojJPJN4dT6Wy2xwAQ1NnyAif2imCeHIYsRQeMxr4Pjtsrrf3j2Db/laoVOvIyXvIsvkuGhqxxilrJvmg3l6GKSAVQCRCbT0I/p7kZEWNCL5Ghd2HVIeCZqwCEAXJzctL8e6xLtMx2X4FT99TayyYVw+1WRSA25YQMotHjPXy1jAAo5fIjkNtxme8C+42scvoZaOCi5rI7k28JUegUdz4kICd5V5XFcDjNyw2EoYUWpgjIhjnL35yFqGxsK0AEIVDQ/uAqRHYhpoivH643Tjf0tI5OHAuaFhnNyzTPBYAroWoTCgwK1l0+wHgXc6S29scex+E+pYgSufOAmB+Hryy4ZV9vGGNeCuVN9SYu2s+uLrCtkBSVgyYLWmzLoaVnMbsJITtlAP7x7ytCX0d8Tbdh03dprCTnaHhNE4xSElgNraOnO/He8e6bD0IHrxHw3cgCORmT1mF8IxUACxhw0+k4vwc+DkLFQBuuqLEEqd8ub7NElJw0xJCqwmwJpX5CSWzF1RogjLeEIC5foFIaZeJMGLEoqffne7FgXMXLYlpGY+aTxgCwKdt/SDEbCiputVmJwBE4UAAk8u+rCwfT91YY8kBsOOZ8I9ngcnCHCzubGfj8cVJWb6odymjNUaZLUB1UR5aLg5DVamFcQMAP995MlrMFSeV2S5fJIXQXZNXwOy1rBgQMLdZ53+nXR5BxlazCyuJ4R6x5bIojHmj7v0T3aZjZcqkqTOEH71+FCoF9jX3WXJ8D6yuwHadTOD3ESMfw475x50no/3EuCI3p5DmwsLZJkNzKkN+M04B2CUxH1xdgeGxsMl6/OBEl2VybfvGBmz+6Aze4yaWUwEUQ12VuR8+9GvzCb0HV1eYKmPZMcX5OVKB6DSp+EnEu8IMiSYT2QL7+c6T+N3pXkuSmedR/0jgUYsub+egOfFJ9HE5WZyy0JZIs62rCuC7X1hhfMdOeMezwERvixdEhBAL+4MvTtr2xEZTDgCIeh98iwoCYF1NEdbVFFnyBVv3tWLTGw2m67gN4/EFiQxO1dCyODhgVZp7m63FgAyicWEn1Jzi/XbsIFnfH54K/ej6StRVBfCzt4+bvqtyv4V5U2Jo7UevHzW8hnG9sO3aRYWmeajo9GdFIC28yjG5AM3Df/VQm22Bl1OoU7a/QSow4xSAXRIT0BKTPCIqLBOgriqAF762Jq4maEz4iEwYGemF789NoE0EGQ1P5EQ/tGaRSWDwZeY81x7QlMrTkgS0OF4nobhoXi58CgGNUNMk5a8VUamFbfPA6gp8cLxLKvyrinLxxI1LADjHP0XhECsJyR9f3xLEhf4RRx66G9RVBfD1jdX4TWMnrl1UiLcbOk2eoVicxH4TLwz8erzYeDQE2H7gPCilJgVd3xI0CSYebsJ4dVXmgkQC52roWASCCc67mZWlGKFFRffmZFtD2oV14u3RxEI8jRcGjL4/YxOqKQzV2jeE/NlZ2NNsbsNMCEzhFb9CcPPyUpTk5xht4cV7/MGJbrx/vMuk9PjiLX684uOhALbsa9WeM6WumHrGbxc8sFS1iJhxCsAuifnsrtMWKw7QJsDOY11GSIMvLY9VJs4ndX0KQV1loelzVYXpob96qM0UL6YAQCmaOkN4+teNmAir2Hf2ojEGg9ETodi6rxUvHzxvKAIABr2SKAQ+aq6+leUk3LQbkFqiqiYA6qoC+NyKMpPbzY7izy0LPVEA5/qG8aM3jsKnKMaubJvuro2ZtHSbCxGV5lfWVdqycmRUP/711n2thtA51zcMH6e5WYxdVgTFM1MmBGmjqgDz1fjwzo5DbVLhDwC/2t+KxgsDeHhtpeN8rF1QAJ9CjP1qVzrUsMhCME2dIctm5yxevumNBqiUGoZI7YICS70LrzDvrC3XFOGvG9HQPhCzRxOfwN70xlFTGJW1Z+HDiM/tbgYhsLRueeKzNQgOj5vWzXvHupCTpSmsDTXR6nZ2bioIblkLCga+mpsHu+cE1MLUUxQCNaIl1Nr7RwzDkvfA+BxgsnMCM04BNHIWE09H21BjLWsngNGvJqxS/Oj1o/D5FMctI3nwSd2ISqUJZDaBxIZlDBGVYvuBVkujsAdWV8DvU4z3KaKKgDETmKVCVYpbV5Rh14luQ3DLEtexYo92lmhY95Re+NoaPHnTEnzY1G3Q5Bjbhj+3U+4hogIRXaGMT6iGcHFSBqI3ZmctiYVECwpn2wp/kWkkuvBivxeWzFZg9iqf3XUa7f0jXDsRs4LP8hGojBnFzz0SjR+LNQI8VKoxyo60abuJyZRAfUsQT/+60cTAeuatRscaFtFrYt8XNzsPDo8b1i1TbLyHwwySps6QoTD5ymHj98LqldS3BLH5ozP44EQ3KNVajohz7/qlxaidP9d0TgpNIVCqVeKOTkRMc4On7/KexzdvWYqnv7gS2w+0omzuLNy8vNQoUmTN2vY290lbULB7tv3J6/CdHZ/idPcl431FNx7FudvUGTJCRhEVRniYnZ8AuuJAynICM0oBMCHLwFf6sQnPC+kry/NxqvsSVxoOqGH37qpdgpAHm0As7sqgEG0CZPkVlAkMEQpt8kRECo3+2fiElbr31E1LsKQ4D5t3N0Ol2iKsLMozCQw3FDg7S/T949F8CeNGi0U8diEpEYRoi4uQqOssKgOmfMWtFFv7hvDLPeek1lIgN9tixcogKkJZjxbxmTCl5vcRS6KZsaYsvxPAQ2sWAQB2N3WjrX/U+KyustDwPD5s6ra/WRzsmo4999EZU3jK7fxlEA0Z3ogSk7KNFwaiSWrOIFleZtvvEQAMr4SFeHpDY/igqdsUUxehEBg1KJ2Do/hNY6eJ6UahVeLK9laQ7f3BFJ02/0O4eXmpEaJjZ+UL3/gkNBDNK5zrjQp/QJvPMmUr7gLHj9v4PyGoXVCQsg2BZpQCeJXj8MosjmUCz391VQA3XVGC5z9u1mKbPgIQgkjE3abPTm42oAloxukWaxEUheBhLpzDW9UrFxTgR2802ApRVmn5oN6HiAknka0kCgwn9gW7Z3agNBrOkoVkmjpDiNDopCaSxDTDk5+tQf7sLIPyyBKtsjiquIh+09hpm2x85q1Gwx3nrVgRzDtiCljWVnhDTRF26c/Ep2iJwXCEGokdcxLe5p4BmJvjxy/3nLNQdJfqAlM0DJwgIyPUtwTx/nEzvZn9LreCpFvwQPjXokA9IjRCY8om22+/+yyBplieflPro2Mn9H1EF8ZcCKapM4SmzpCJvMGflxlYMqaRSDT4wWtHTZ72Lz5ptrTn+On9q4zf+/LB89i6rxXbD57Xn78KUEA0y6gNHVTcX0SGcFhFcHjcNYU2XswYBVDfEjS57RRWAf2AzsJhgrZ2QQGeeasRlGoWytNf1HIAOw614SX94b988LyF88zQ6JCgUwjw2HXVnIsJVBfPMVxHqocoAE0IfGHVfBw+3487a8vR2D4gpXXy52YTn1lvLx88j2sXFZqOkwkMp3j6A6srpDRYAKY+PrL4OR86UnUhTKn2/4gePiAA7r12AfJnZxnfZfkOXhnwwktcRNcuKkTn4KjlOJ4KqVJqW4nKx/YB4Osbq/Ho+kppW+Ff6Z4OKz6iiPZFEi1jZjgQLgegQFPK4v30K5rw2rqv1UhYh8OqRbDwWFGeL7X+RaaOAuC2q8pQkp/jcDYzSoVjxdd1VQFLDQsDUzZXlOUbFdoATF4RG56sEJPBpxD8+F4tPMMrmXcaOjAoCWcSADlZUcvejmkERKmt4tVHJiKWc/K/d4JTDk7+vp2yFWnR5XNzLOQIVhz5zx+cklJoJ4sZowD2NvdZio5kQoBw/xrbBwyhQXWhwZJyzDKQcZ4ZnOw2lQI7j3cZkz6swhD+CtEm/JHz/fjHnSdNFtHm3c1YV23/8AmiCUgxUXzgXBB+RUsIxkoaysBosDKO9UNrFtlS28TQEQWwurIQy8ryjWZwKgCfArzd0InwkXbpQuWVAVu0/KY8BJoXt25xkcVakm1IIysEEz0KpxoPNjZmXKi6A8DOKdJVZYqM9y58Oivlw5M9BmNHIYBfIbjtqjKcvziM450h6VgGR+XFiCxkMT6hQlEI7r56Pt76tAMqpaZKWifwit+nEAubbeu+VmzdZ86JaPFrra6hpmSOFsbwRUOAdmtDFi7zEeDHek3Jb4TnExqZQO8la46EEE15y2i/dr2CNt1da4zRpwCdg9GQnE+BiZXVcMFs3PkUYtnfgr0vept8UjuLC+1Uzsu1KAAA2JPCuoAZowC0VgTRQiyZxdreP2Jk3yciFN2hMUvs7WdvH8ebhy+Yzm0XGnlwdQW272+1DdWMChYGoFloqxYW4HhnyFKZzDAWVi0JayYoeDooY26wOgCWHLtD2HM2EdQU55leM2/KpHT0xdYoLBYA2H8uiMO6JcfueViFRoWBfMcq9jffi+XpL67ErCwzl5tV6O45oymJ5eX5eKehw5RcYx1fxTyBbHP2WJ1VT3WZd4D6zo5P8dj1iw0+ujh+AKakNe9d7G3uw87jXaairXCE4oMT3VqHWlhDDADQMTAqDUfyIZre0BjePNJuCFi3hWR1VQE8fU+U7cMnkFkLZHGK5+X4MDwewemeIZzuGcIHTd3GOiEKgY8QqVfDxsZbxkDUWGsVWn0cbrPOLXaeFz85i8qiPBPtl+8VxCeVJ/RQC8tfiQYOG5fY+I0AxjxsaB/AK/VtRvsMQPN2eaquyEQzfj+lCI1a+0gBQOncWciWtNhOBmaMAgAQ7ToJzaJv6gwZsctwRHsgjJZFAexq6satejvXkvwc/Ouec5ZYY6yCKuKQ8cyflYVs/7il98/KhQU4KhGaDBtrijTLlKOrrVpotur5mLeiuzSqCukWfDzEnaDEUA6bvKKV9j92ncLy8nyLpf2yUNjGgy0UWR8m2Y5V9S1BPPPrRlMvlob2AZOlzdd58Mwtdi2FaPF6WT4BiG7O/vrhC6icl2uh3IqdVcVmcIDmyX3/taPYf7YPP//KZyz3VKyI5SFWWovjtQOl1poVHq8earPkGXimEX+PxSI7AqAnNBalr3JW9JHz/dIcxaUxs3ET1jcOotDCm1/SWzk/99EZDEoEH39GvqL2ztpyKYtIBrapC6OofmWdtjZYryDTvdCvwQynT9v6zeOh2j1cUDjb1PiN1a6wdcfybqygk8K8E6BpLunzmEILgdqFJW9ZXoqnblri1QFMBloBR/SpR1QYfb+p8R7F4pJoHD4coYYVwIpceBTmZuFf/mit7QMRrylicHTCCKnwTboAvenUhNVC8ilAaCysNZvSwaiAjR0NaGwfMIpaojFvGOwXp5bLfKdDnwIoipnyyk9eEax52YN6z3q2QNgikEFRiHE8C3kwEJj7DIlN1fjjePB1HoCmIyPhqLV2/dJi3LVyvjSfwO7Bi580I6xqlcqPvLAXj11Xbdw/lvxnMW+nMN/rh9tRPneWiZUkdtbc/NEZ7NKte3afxSZ/gdxsPP1mg2OMnELrJskXLfExbll8/nNXllrotExgslCM7FlTRBW7+LHYTsX0mU559SkEvaExEGg5m92cx8XWGX8GlTvf7bXleOHjZkcWGQ+mtBjtd0NNEXZwxWfGb6LRe8BYZbLfLSro1ovDeOatRgAw1u9P719lmtORiDnsJMsN+RQiDQExqqnbWpd4MWMUgIznzydSCTTZ2NxzyfplyBfCzVeUOD6UWGX6HQNajFHW/5sJgdDIBF46eB4Xhyf0MWssDJ+PWJRLOEKxZV8rXq5vw2PXVZuSbEY9Q1jFz3eetLRwFq1rPhzDV35m+xWLJckwNqGi4cIAsvxa0pIQTcjLBAXR47psfPzvYQrk41O92Hf2IrZ9Y4OpCyOP/By/KTxz47IS83VgptXJwi78fdhxqM1UyDMeVvHiJ2cNT4pCsyDZrlYywcrj9cMXbDtr+hSCD050G/OQhczEjVJMeybEgMyrEYUWS8w+edMS43sslMNTnp2uodpIYDvhTxClvL508Lw0vMmEP1M+/HPjq2ddyn6DR0+puXiNhcS2Hzhv3Huq7/7FvDvxPIx9t+NQGxYX5eF8cBhDYxEpRXnT3bXoDo1B4YgOLMog1rMAMBhFB84FoRCgcHYW+kcmQClcb5KUKGaMAqirCuDH963S6JP6YvYpBJEIBSF6c6vWflOhDg82mfgJviwGtznW9oYqjTIQRC+Aj3fLrF57Up0mSHaeMHPHjZgkzIKVWdcPP/d7SwWjT1+R/OL5+sZqWxecAjh6YcAQuqzoa211wFIEd9uKMiwvz7ds1uL3EVQX5RleGCt8k/ktBNYEWReXuAOAhYWzTPx61hHSzqISr0NItEEdE3zMonx4nSbQtgneC49F83LRPzKBcf039obGDAHAKj958IKC33w9LKn5kEG2QQ+fkN7X3IdPLwxYjBeNJOEsWu1UkMxqF7+Xk6UYtMuIRHkQRM+hCufiw5axjBAeWX4Fj11n3S6VPfuVCwpMgntDTRECudmmHNBTN0YpyTIvjEA3cvR7Nzqh4gevHTXGTwhw8/JS7DzeZepgu1Bn+ImMIkphGHvGCVKIGaMAgGiVJHvoqkr1nX80C8jJraQArltShP0Oe3bKeovbJe0YQiMTJiHIFouslTCPWNxwvhJRBiZY66oC2PzRGYvwF1thM0Xx4idnbc9JYLUcwxGKi0NWRTgwPG50tuShqhR52T7TexRabFWkoFJoLB2+r8/i4jwTTTB/VhaAqAJgHSGbOkPSXk48FZgAuPeaBfhNY6fFZec3l3dK9N//mQrc/5kKgwb77rEufHiyB9u+oXVPffngeUOoUMD4faMTqlFdvaGmyLaYjMfa6gBmZfmkXPG6qgDea+w0wi0sPMWa5vFsITYW/nK369RRVjzl0wkH+Tl+NHYMojc0ZmEoEWiUUzEvwZMxGBTdW1RVa66jal6uSXhvedzakFGGysBsvPi7s4hEqKmzJ185vv1Jc9FiU2cIS0vyAELw2PWLDU+xsX1AuhbL5uZgQ425BbkYWuofHjd5fS8dPI8IR2LYJ/E6GNxstzoZzCgFAJhL13m4iSnuPtVrWARi6EDWnK12QQFIDA2wp7nPUr3IQgE9oTGpVcWHEURrKR4cagni2V2nceS82Tqfl5eNF762xjLpnKzE3GzFaMwlHiEWEgHQvC0JbU4hBBtrikxCfG6O30JBZT3XqUpxTVUhekJjuLO23NIArGtw1NRymlJqEh4fn+rFh03dxgYu7BhAu6//frQDf3L9Yuw83gUQgtuuLDU9/2d3nTbz7DlBTRD1Avn5xbce4Ju0ic7ne8e6sHVfKx5dX4nyuWZPRoaDLUGAQsoVr28J4ldcFTygFc0xBcCHRsSOtD6F4Cn9/vBNCZkHNzZh7e/Ewial+TmGofSN/30Q3YOjWF1p9giZ4UAoxedWlOH9E10mynbLxWFLV94XvrYGX978e4tnyeN0z5DxN+vsecvyUlPl+E/vX2WEYMX4f2vfkKlORxZ2HR6PoKHduR1868Vhw+vj2UXjEWryFmQgihcCSipCIxMxLSkn7Gnuwx215cZrnkLKBDkrgY/V8gCQd3OkFKb4pIjS/ByTS79bUk0osxjF9453hnCis8kyAZeU5Jniz8xCCo1MWPr3M1TNy8NP7l+FzR+dwfu6u8sgY3nYcaZZzxQ+3PDiJ2dxO3fPb15eit2nevQKYZi2gLxqvjksZ3KnobX/6BbCRO8e68IHJ7qNa/PW6USE4rndzcZYTndfwk/vX2WKr/M8+8VFuTilCx4KLQ/U2jdkuh4f131gdYVjCOn53WfQ2jcUU/gD0eci7lIn7jTHcGdtuYX1c/7isCVEE1Ep3mvsNPXlYVtF8j11eJTNzUF3aAxb9rXiV/tbAQKTUM/yaeFX/rthVdsjQhUMpnCEWvJW9S1BHD7fH/Oe8PjgRLfl2fPV8GJ/py37W0w5sbXVhWjvH8EF7lkMjoYxOOrsbXcOjuFHrx/F51aUWXasiyWKYmwqNmnMKAVQ3xLEC5INLOLBsY5BHL0wYGkSJj5JCndehTjZje86aKnaBQVGHLO9f8T02dKSPKzXW0+zfAcQ3bnogKQhnYgr9NyGuYul2W2vLso1bYWYP8uPps4Qdp3o1iqnieZJ9FyS50FYZ0r++gV6UZdWs0G4hKRWtPQKFy6pLsrF9UuL0XBhwOQtZPkUZOnFPDKsKM+3eBhAlDL4+A2LLd8Rz7T9QKshNHjL+T9agpYwyPYDrRYLkRXNAcC/7jnnGCJs6RvGcw60x4rCWQBgUhAU0UK3QG62qZmgCJGTbsdaE/M+o6yluo1BwLNZZIWyWfrzFb8rK4SiAD451Yvfn+5FXVXAyL25bZFhnIdSy85rrM5jb3MfcoR2FaFRM5W1viVoazwuLcnD/MLZqJ0/F//yu7OW+cfCf06Yl5uFgdEJk6JkyWkvBJQEPPfRGdf0MRHlc3NQOS/X2F5QbBIWCz6F4PolRSZr3adAa30co8xfBN9iQgyvLC6Zg7+5f5Xxmk9yuUknsbbGbFFEu1iajxMtqf3ngjjILRCVAr02wh8AArlZWFych0Pn+w2hc3FoHN9/7SieurEGt1xZalib2X5F2xWNe3jn+obRcrEV916zwCTMN9YUoXZhAXpDY0avHiCawNTosoNSD0mlVAv1xIAmRKL3KDQygZcOnJcKJFHR+PWmZ6wK+cOTPabPZ2UpRtvklr7hmBainWeweXezbdKWge+b5EQxtcPwmLxwydV3XSRxgehzYwbV/nNB7D8XhE+Bkfvh22uIWKE3dGRr4JblpQbzSiEaKeDpXzdG64Acci1O6/y2FWX47hdW4Nldp2Mm0+3w5TWLcHttueFFaywgb0OYpKC+JehqccuQ5SP4889dgaffjFY8+nwKaufP1atN5a1qeURUagnV+AjB0/fUYldTd8yEFg/+MgMC04h/zbu242HVNl7qU4DPXVkGQEuSbtnXipcOnscfrJpvO+llC5g/lsRIWvZcGkfPpXGskzCE+H7uK/UCN1lfJUq1TXyeurEGjR2DqJ0/1+Dc+xVi1CKIfZZkVi4BdHaJtTpbxC3LSy07y7mBj2gtMJjAyfYruLIsH4eHo7/tyrJ8DI6FcUVZPtoHRk1VpfEi1vc03vmodP9cNxD7DKUCdqePqEDBLB/WVJfg5uWlRl2HmBNbNC8XP7l/lRHmevVQmyGgVcEqD+t9+eMBO/yXe87h9tpy1yFmhQBrqsxzn7U1KcnPgaIzFFVVI2sASIkXMGMUwI5DbY4PhgC2MftwhOL1/2gzLDwCbb/gX+45h7BuSXzjszUAYOwQ9ZvGTmkhFw9W/ScTOj5Fa68wHlYtYQW+id2Y4Nrzr7fua7XdSco8Dm3SsWpPQIt9vyFUPbtjomvezjduWIxf/O5sTMvy4tC4Za9klhBlBW7HOxrw9BdXSr8fVinO9A5hQ02Rqe8+X2VJqbWxl2zMm+7WlHGsePuHTd1obB+IS/gvLclDa3DE8CABSIX7kbYBo72BQjRPqd+lUFlYOAsdA6OuhfLvzvThx3reIzQy4brCliHH7zNCeenAxeEJvHusCye7Qkb9x07dcmYozs8xUX5/aFPkxeicdr/FzjPgn+WOQ2146cB560Gwrh0CoECoE/qoqRv/7f1TpnBdWN8nYPuB85Y9tpMBJzr5ZQXZphoryvNxx1VlWptnwJZzSwHd7STwEY2mWZqfY1AYVQq8oPcdeWjNIvzhxmpseXwDVlUUSM8HmFvyil0511UH8NKT1+GNP7sBP9GrCnnw/eE3Cu4he80Ke9xadfua+ywMGtOEjcMyWl1ZiMaOQaMfzlfXV+KOq8qkx9aUzMH2J6/DivJ8+GyuMR6h2NXUjbU2TfDeO9aFf3i3yWjIxsbu92nPK8uhFTGDqivjp7jiKIbZWebvv3usC786YJ+4FZHtV7C+psjYoAeI1pWIzeb4c6pUE3Ju5euFfrPwJwDKHbp+st/8zVuWGp0p40GWT9vlKhHkClTfWFCIpuDEXb4ALRz47rEuvK/nn3ic7grhB68dRX1LEPUtQZywaaa3tjqAZ+5dadrZjYdTMpatZQL7vIRCgDuuKkO2XzHmZGl+jume7z8XtM3VsD222d4DycKM8QCKJQvhRGcI2X7FiB9SlTrGAK+aPxd31JYbMblfcfzviP6AWMx6y+MbsHJhgSUGzMAsTkBjAjELQSHAd+7SqHksiRfIzTKxWZp7Lhlc5tlZ5oXE3Mh4+sgDZsqcDFeW5eNs31BMrwaIsnI+5mizQ5J4sUKAJ29agvcaO227XDLwFbMyqNS8SAm0ZOvCwtkIjUxgy/4Wx/Oz3kNNknHk+BWMCCEvl3VZmJeXjRuXFesKSTFx6E91hRxpjJNFlo/g0rh9nJ51LgWibCY3BVYMNSVz0NzrPG/sMDweO9TGQ6WagvMr1tAJQ0SlllYULGew/eB5VBfl2aqrQ639+M5dK/Dw2kpjZy4edt/z+7R9O2oXFGBXU7dtUjxCNUry0/fUGuu9dkEBsnzWPlh2kO2xPVnMGAXwoIRux5KCQLSb5oSDkCHQCndCIxO4vbYci4vzTILTCJ/oPG+xbS4PlWqNzJ55q9G06FQK/Oyd4zhyvt92YpzpGbLtV8LaT4hbPs7J8ZkadF1TUYDFxXnSjTSk1+wdwmPXVWNPc5+tUpOBp1DyIAC+eM0C7DjUhm37rQtOhBpnmIEiGir7u982xTyeQFPEH0jyRP0jiSc7Lw6N4/XD7QYv/ivrKg1hIQoxN8VecYEQrK4MSGnCgNZygO9quunuWvzik2ZbY4APY9x+VRmevGkJdp/qcWUUJAthVQudyO4VgVas+fGpXsucC0eoY3FkRNXYNg+srjAV5/GoEKrKr6kowKZ7NCPu4ef3OPb9AjRZc/TCUS0/pmrz4eblpTHZQTySHWxLmwIghCwC8L8BlEH7Xc9TSv9bqq5XVxVA5bxctAg8XIbKebkozc9xtMhYbHbz7ma88MlZk1BijB5+izlZXxEGQqJ7fYo4030pZuMv6TmhCbFnd522tIK4NBZBll7IkuVXsLGmyJFeKCKs98SJ5VWIsU67oylgCEY3kzreY1gTrXcbO118M9rgzI4+OllQ6N1LLwzYbqojejBXluejY3AU/cPyXv+xEImoqJiXa6tYyubOMhrAsVzWleX27U2y/NH5fcvyUkNpBIfH8e+ftuNYh7MXlyzss0k+U8BW2bnBu42demtvOUmgpmSOSQFs5HYSiyX8GVQKY6KOR7QQnFMDPR4+AqP6PFlIpwcQBvBtSukhQkg+gHpCyHuU0mOpuFh9S9BW+ANatV5rn/lz1n9L9mjEcARj9LDujay3j93ii6gUQ2NhaV8TkY8swk5oMiEWkfCrAS0k0hMaw9meS7aWuQysJbGbZJ+iAHWVAYyFVfQOjeNCcMTxeDExNlnxy5L5bG+AC/3y67NYMl+1G6/wz4rhMcrg5D2JijNWWCwmCHHcUD44NI4fvn7URN091hGSPod5uVmYl5eNmpI5BuuGr3qfk+MsSornZDvSgm1/gmQsssLCyYJFA+yej0y5sPyNrNLdLQ619rtaV2xHtGQzgdKmACilHQA69L9DhJDjABYCSIkC2KFTqUSwRlayZ6AQ4HMryvDhyR6LxSYqh7CqJSq7B0dx9MIAVKp5BU7Jo9cPt+OpG2uw80S3yT3tkBTDGNcFsKR0jq07yzpvyhJ6w2PhuOim66oDuO8zFYZS47tF2iGianHXbB/BkpI5uOB4tAYWR31gdQWaOkNa8tqG181gx3j57LJiDIxMYDysOpbZ1wjhu0QgCn92z9PDibEiolJ8cKLLNqzUaSO4ZIdfHJ7AxeEJnO4ZQnPvkKl+YOu+1pgJ5HgTzE5jyRTUzp+LZ3edxqkuq6J2a8y4CW3eoYfbLlsaKCGkGsBnAOyTfPYEgCcAoLIycQqUzBJSCPDEZ2vwyz3npJQ+Ro9k7Yh/pW8k4SPAj+9bZeLvqxQW4eomUbinuQ9nbVpQMywMzEbnwKjWS91HcI5rLSBONL7n+oTQhtltvJ9hLKwa+wuwyRerdwnDeIS6smAJtAKYB/SNNHpDY7j1ylKU5Ofg/MVhqUvv9xGL8C/MzcLNV5Tg7aMdrpJqvUKDumQIGtk55gkJfBHJ8HqcEKNbdUJgxgdfpBULdhXh0xXl+TmOXrTbZ+r2uMu2EpgQMgfADgDfopRauipRSp8H8DwArFmzJqG1Ut8SNFEno+fWWDOslP8VSWyWwnzzCWAIxODwOHYe65rUAp6IqDGpmnwYhXKWMQFwdUUBPm0b0BhE0PrkjE5EcNfK+fin909KS+vdgrnEvzpwHj++dyWWl+dLOzlOBhRaszexHbWM7scQmJ1lEigKAf5glUaldcuokG0knmwoBHjhj9ZKd5JjIETeRXU6YIle26BV46ZG2Uwlqotyce2iQjRcGHD0DgnsvadYyM32xc2A2ne2T7p/dTKQVgVACMmCJvy3UEpfTdV19jb3SQWDohCD+nf+4jC+sLIcbx5uNxgNfp+2YxW/UxK/BaTYQjeR+R9v/FALLWn8a59CkONXjI0n/ArBRyd7EI5o2w2uq543KQXAEFGp0cwqFYt8T3Of5byxqohFsE1a3DJpZAySZMtgxnJ662iH7TEqTTw+PpWY7VcQodS0jmpK5uBvv3SNdA/d6QYCrTkeiwY4we/Qa8p0nMRYilf4A8DASBj/8G6TZae3ZCBthWCEEALgXwAcp5T+11Rey66XhqpSbHrjKL7/2lF8fKoXr3PCH9BK/ps6Q/jh60cRVqmx0Tlj97ANUkASE/5A/JYf65h564oyhFWK/eeCRhKpsijPKIcfnVDxyWlnRkQ8cdkIhSPnm0Dr2x4vCGBpaBcvWK+iiEoTXhz3XrsACwtnIdtF0VgsLC3Jw7rqAN5u6MTWfa0xGSKxhL+C+Aunko2RsGoxovjamragPcFiOoACeO7jZlfV3W7JAsk0lkTZkyyk0wO4HsAfAjhKCDmsv/d9SunbqbiYrC0/hfNDGhgex484lgSgsWE21BShviVobCs3Gff9xmXFeONwuyO1k//suiVFWF6er3sk0fdVat0EJta4br+qDCMTEXwi4U3LkOUQl/H7CG5YVowWSRGNEyiSEx8m0MIph1r7E/p+vPkROygEaA2O4EzPUNI8ChWJWY6pRLa+RWK8/ZAyGaluvTxZENgbs4kinSygT5A4OSAuxLOPKI+LQ+Om+DwhwOM3LJZumhEvcvwK/vi6avzi9+ccxyZ+tvtUL2Zl+Sbdf4VA6/3e0ueeCXNCwnZgmIhQnHT4PNVYUpKHM71DxmYu6QKz1NIFxk5jDLb+4fGkVxsTAjx2/WIEh8fTKvzXVgccWzRPB/h02vTFoXGtziA4bF9PkYLtIWdELyDWX94JfuFzn0IwL8/crKkykIsXPznryq2PhfGIikOt9r0/nLD7ZM+kFx3jmQ+7XMCMLuuEZPcpiQcD+ibaqYQ4g7LsmhelEdcsLMANy4rx4/tW4fmvrYm5b3UioFTrfRUamUir5V/fEsQTn63B0tI5aRzF5FAwKwv7zwVxumcI7x7rQo9DTlDVq5WTiRmhAOqqtEZPTsuV57cTAPdcPd9yTMvF4bg3obADpXBtmRXq/X0YRqfYwvQrJO7EaiyFm2y4DSMV5Cbu9Iq3IFVVw5PB4bYB/O50L555qxH1LUE8sLoCfgdF5fYxiU3SIiq1NA+caqgUeH53M3pCsXdLSwTzcrNiHzRJiBRhp3ks24d8spgRCgAAljuUuIug0GLCB9Jo0fLIy0lvAjARpZdol8hUY3DYWkW6rjowNbHIKYJKtW0hv/3SYTR1hrD9iY1YWx2wdFuNh34qO+7CwOSS98mACo0lk2wQAHNnp14BxIPHrqtOOg00pgIghJQQQv6eEPI2IeQD9i+po5gCPPfRmbhF0lSEk91YYBf6R+Ej8hDgjcuK4Ve0vkKZpM0zlRMuPtL7rl2ApWX5jj1weEwXRUGhtUn+/mtH8a97zuHAuaCl3iSe6U2p9bf3hjKbujoZKArBndw+1JkAsb9XMuBGZmwBcBzAYgD/BcA5AAeSPpIUor4liPcT3A0s1XCb14lQYEHBLCO0QqBZrn9x2xV45t6VGv89dcPMaLD+P4l8743D7di6r9V1353bbfY1iBd2fedTgUR3whMRmIKQSKZg/lz5PgqyeTZVuaARh9beicLNNCyilP4LgAlK6UeU0scA3Jr0kaQQU7F1XaJw21ce0DwBxv6h0HIIDz+/B7/4pDnhvY5lmJc3fRY6642UyPN128aAR6zaCrfIj9E8TcRkRAzfBnwycGppwWMqlVuq0NY/Kt0hTRYVmKpckFN7+UTh5lGxp95BCPkDQshnAMxL+khSCLE3fjIgJmbThXCETrqpmYhMTG7aQVEwpQTuZPHx4+1oOX2eSHxGzXRDup4D2zwp2XBjhvyEEFIA4NsA/hnAXADfSvpIUoSt+1otmjw/x4fQJK0iJ2ZFPHBqQeD3aRtDT/WkC02i3W7SNzWJgYga3c2MhYKmg/7K9KIjD5mFVJHq3HgAQUrpAKW0gVJ6C6W0DsDF1Awn+dh+wFqZOlnhD8Qu358MCmf7cU1FAZ754koEhFqETEc6Q20LCmelpFgmFaDQdpjykHnIxBBWWLVvaT8ZuPEA/hnAahfvZSRiba6SbsjkZf9IGP1tAzjSJt/20YMc7QOjGWlZEwI8+dkaNHYM4mOuvTW/u5SHzEGmhrBOp6DS3lYBEEI2ArgOQAkh5C+5j+YCSC8xPQ4sK8tP6cbbU43sODaRTibcbls3leB77RPEDqv4lPQs7iyfgttry3F7bTl+f6Zv0m08ZCic7cfgaFjqgVXNy0VYVdHePzqtcgkezNh/Loit+1rx6PrE90UR4WQeZwOYA01J5HP/BgF8KWkjSDHiZVtkOpIl/PNnxXdfMk34A2ZWiqvRpeknTITVqPueIhdlYCSMG5YWSz9ruTiMCzNA+E+P4N/k8E6DfWvxRGArBSilHwH4iBDyS0ppS1KvOoVg+3ZebpB1N40Hk0n0OoG46BmULqQrOUyh5aJ6Q2MpGwPF5DZEvxyQl+3DpQzrmuoGvjiIC0VJzgm6MQOHCSF/B6AWgJG1opROi1qA2vlzTXHXywUZGqbMWOGfbkRUoLHdfkN4D5PHdBT+QHyGyVmHPTkSgdtK4BOYppXA+RnC1/fg4YKX9PUwSSSb1HLZVwKnogiMIRPbAXvw4OHyxaWx5IZu3YSATJXAANoxjSqBk9UHRYbpVDHrwYOH6Y+LQ8mtP0q0Evj/SuooUolpUhjkwUOmIcdHMOYZOSlHPEngKd8SklL6lv7nAIBbknr1KcBtV5Za9sr14MFDbHjCf2oQz23OSzKt3akQ7J/hwJymlP55UkeSInhJYA8ePEwn5GX7MGTDaOp22DIyETglgQ8CqIdG/VwN4JT+71poRWLTAhtqilLWSGm6IjfbNyOKZjx4mI4YDUdsm02W5sv3KUgUToVg/x8AEEL+FMANlNKw/nozgI+TOooUoq4qgPmFs3EhmP7t6zIFyWpp7CG9mOrOqx6mBqoKLAjMksqsZHc2cEMDDUBL/DLM0d+bFqhvCaIrA/Yu9eAh2fCEfxSXk0dLAbT3y2VWsreFdKNOfgbgPwghu6Dd5xsBPJ3UUaQQmbwbmAi/krl76XrwkMmYJkvcNWwr6pNcau+GBfS/CCHvAFivv/UdSmlnUkeRQmyoKUK2X8HYhJrxk8QT/h48eHDCbSuSsyc1g6u6YkppJ6X0Df1f0oQ/IeROQkgTIeQ0IeS7yTovj7qqADbdXeuVA3jw4GHaI5TkSuC07ZZCCPEBeBbAXQCuAvAIIeSqVFyrsX1g2oSBZgKK86cNicyDh4zCqSRvCpPO7bLWAThNKW2mlI4D+BWAe1NxIU/2ZxZ6Q6nbTtODh8sZF2ySw4nClQIghNxACPlj/e8SQsjiJFx7IYDz3Os2/T3x2k8QQg4SQg729PQkdKHhJLtNHjykGwRAdoZvd+oh+ZidPcU0UELI/w3gOwC+p7+VBeDfkjoKB1BKn6eUrqGUrikpKYn7+/UtQbx+uD0FI/PgIX3Iy/Hh6XtqUTj78trxzoMzVi6YG/ugOOBm9twP4DMADgEApbSdEJKfhGtfALCIe12hv5dUvMq24vPg4TLC0FgEz7zV6OW2ZhjSsSHMOKWUQg+lE0LyknTtAwCWEUIWE0KyAXwFwJtJOrcBb314SAcKU9yDigIYnVARTscu9x7ShmSH/dyc7SVCyHMACgkh3wCwE8ALk72w3lrizwD8FsBxAC9RShsne14RD66ugLdvS/IxJ8eHddXTpiB8ytGfwo2IeHgewMzCFWXJCL5E4aYQ7O8JIbcDGASwHMAmSul7ybg4pfRtAG8n41x2qKsKYOXCAhxp8/ZjTSYujUVwcXhqhJwHDx40JLsXkKuzUUrfI4TsY8cTQuZRSi8mdSQpxMNrK3Gk7Wi6h3HZwdtnwYOHqUVjx2BSzxdTARBCnoS2GfwoABUaA40CqEnqSFKIR9dXorVvCK8fvgBQoDPJPbUBoKJwFtq8Tb89ePCQQtTOn3oW0F8BWEkp7U3qlacQ9S1B/HLPOYyHVfgVkpKma/1eOMSDBw8pRrJbQbhRAGcADCf1qlOMvc19GA+rUCkwnqJt7i55PfY9ePCQYiS7FYQbBfA9AL/XcwBG7GS6bAkJ6LuCKQSqt8epBw8epjHGkhy6cKMAngPwAYCj0HIA0xIRT/h78OBhmmNjTVFSz+dGAWRRSv8yqVedYmz+6IxXEObBg4dpj/wkFxi6KQR7R2/INp8QMo/9S+ooUozuQY+d48GDh+kNn6KFs5MJNx7AI/r/v8e9N61ooIuL87xCMA9TirxsH4Y8YoCHJCLJu0ECcFcJnIzWz2lF31Bm9Z9nhRQeLg9UzctFy0UzUc4T/h6SDZVqjMa6quS1YHHTDjqLEPLnhJBX9H9/RghJbaerJKMoLz07UPls7m7OFPdxz832Ten1MhGpvOWi8L9ckJ1gEy0C7X4vLcmD14YruQjkJleWuVkW/xNAHYD/of+r09+bNkiHB5DlI7Br1Dg6xbu/2ymiyWIq9dhkBckU3/LLApXzcuNWAj4CXFmeD5UCp3uGPE83yWhoT24o200OYC2l9Bru9QeEkCNJHUWKkQ4PYCKDaKd5WX6ERpMfkphKoZo5dzOzkMpw4pneIVQGrOEtJ0QocLwzucVKHqLoTXIbGzc2XIQQsoS9IITUAJg2Ac76liB+fWRm7wg2PDFtHpdrZF2muyHG6+mkUjFSmprw1uwsBcoUxIYUorUtv5xQnJ+T1PO58QD+GsAuQkgztPlZBeCPkzqKFGJvcx8yyBgHMPUMkcHRy29P5InLNKRDoYUPU+VBKkQT7OlcEiNT9PBUqrUtv5wwd6rbQVNK3yeELIO2FwAANFFKk99OM0XYUFMEH0FGKYHKebmem+zBFqkMH3obyExvJLsdtBsW0EMAsimlnwL4IoBthJDVSR1FClFXFcDikjmW9/PT6BoOj0eQ5W1T5uEyRaIzeyrCQtMdyW4H7SaS+iNKaYgQcgOAzwH4F0wzFtC8XCtrNTRFrqFPIZatE1suDmdUkjgZIA6L11vXyUN+ji/j72eiM7umOA9zZyUW4nCaf5cTkt0O2lUSWP//HwB4gVL67wDSQ6xPEIVJ5s7Gg4fXLsJ9n6lI2/WnCnZVigoB1np7BycNl8YiUxq/z0tBDYkdtfR0z1DC+SpKtblGAPh9BAsDsycxwsxFsp+9GwVwQd8U/mEAbxNCclx+zwOA010hNCaZuzudoFKgvjWY7mFcNkiGAIjHWJaRFSYbqnGzJ0cixZKqntymKnWVLE2n0zA7QRrbygUFSR2Hm1F8GcBvAXyeUtoPYB40ZtC0h28Kgo77zwWxdV9ryq/jhNlZCu64qgyFSe4k6BZ2BXGJYoZ4+ynDI+srJ+eVUbjOYd24rDghhTGZvvcRCpxwQbKoKspN+BqTRV6Coa5kF4LFVACU0mFK6auU0lP66w5K6btJHUWKIePOKgRQp4gSke5o/+dry/FhUzf6R6bvtpW8DJnu8V6/IBEJpk6pEQIMj4Vx4Fz8Xhkbpwr3TKWPT/emhXkU65IKAZ64cQlSxcWIVX3fG0qsO8HpJO8INiNCOQ+urrA86LqqgK0gmXuZFY/k5fjTlnReUZ4P/yRX2Y3LilFVlGsIyVR0RZwMFADz4qg2DwsScWlJXlJYYbNdhE0oBd5IsDDyyvL8uL+TzmfldEfXVAWwvDwf1y8tTsm1k+31MiR7R7AZoQDqqgK455oFhgWT7SM40jZga5kMXmbFIye7QpOK205GNJ3tG8LqRYWTOAOwp7kPLX3DoPpYfD4yaaWSDCgA/vrzy/Hyn16HL9clnujvGRrHzctL47qujC0z22WRUDxCmc0bvwKsrgq48mZL5mQGR8TJUzzYEsTDz/0eu0/1JnZuwMLumwo8vLYyqeebEQpg675WvH64XUsQAVhaOgcTKWpkk36xZEV9S3BShXCTMeLGJ1Sc7rk0iTNo4QY2Bgpte88rSq21HW5x1fz85DwnohUa1lUFsKe5L+HTXL2wAB+c6DJeO4UPvrq+Ej+5fxWGBDogIUC+SwUQz29ftbAAj66vxK1XlrnekHxhYWYwcFRq/1tVOrleVo+sr8R37loxpet9RXk+Hl2fXAWQ3LriDMU7DR2m18c6Qpp1kAL3VFEAgCAiuBfp3AMgHTFYRf/BKoCLQ9bcQ7z3gz+eQn+G3Ht21d4Ksf5+SjUCgBiKiRdU78/e1BnC0QQ3HFIIMCvLZxJGi4vn4HS3XGme7Aph39mLlt9KKRBR3Uk0Qtx7AQWzs+ImMRxO8+ZL/DNPxdTP9it4cHUF6qoC+Jv7V+H7rx1NwVWsWDQv+UnrtCgAQsjfAbgHwDiAMwD+WGcYpQR3rZyPjwVXL1WxSS32R6FAE36XO+wEuUqB0vxs9ITGpZ/He/uzfMRCH+Rf2Xk4Mhl/vDOEaysKJi2oKIC3jrTjRFfINJZsyVjtoFLgYMtF03sjDs37nJK3bu+p27m/rCQv4RBJIiic7Uf/SBIKnahc8YuIxwi5cVkxKuZpeagHdOFf3xJEQ/sAfErqYv48kt0IDkhfCOg9ACsppVcDOAnzdpNJx6PrK3HftQtM7xGk9sfz88GnICNi1vEgy0ds3dtsH8HS0jn46f2r8OSN9juDdifIdJDBrUB1i9M9l5LSeuB4Z8gkUBXEP1bRQ7oQHEloLLOz/bhxWeykptscytm+qd3oJinCH4Dfr+CJz9bEpHn7fMR18n1gZAI/vX8V/ub+VairCmDrvlZ8+bk92Lqv1ST8FaLNAb9CcN+1C3BNRXy8fbvh+H0ED65OfkFpWhQApfRdSil72nsBpLRUtr4liLc+NTMffD4CX5KFcrZP3uY2omqUUwLtAWcijVGs+OTj7iLGIxSnuy/h+68dxfMfNzueN92EHbYgRQyPR+BP8k45CgEWFM5K6jl5lOfnOG7Cc7r7kiuLfVnJHCjEXsEzuAmRJWsqJ7MkZ4XOVqIxxv+ZRYWgLt2hbO7G17cE8aPXj1rCvIDmXT2yvhLP3LsSbx/twKdtA5Z7tLBwlm0C2c52+PKaRUndCpIhE5LAjwF4x+5DQsgThJCDhJCDPT09CV1gx6E2S8Kndv5c6QOMBSedMR5RpW4nQbRKMUKnnhrnZm0NjUcSWoSZ1F2yeE62oVz9CsFTN9bg23csxyqJFUYp8KW6CuQnWJAjwxevWYCrklypyePqRYXY/uR1WDqJBDgAnOgKIRxRE1LOS0vyLO8xSzvW9HESNvMLZ7vauW5ddSAm++bTtgFs3t1s8sJl21O2D4y6TgQTAM/uOo36liBePdRmK6gpgNaLw/iwqRvjuhElHtodGsN9n6mIKyqQ7DbQDClTAISQnYSQBsm/e7ljfgAgDGCL3Xkopc9TStdQSteUlJQkNBZZ8cTi4ryEhFcikYilJXmmSeDmsc/2K1hakofyJMT97r12AXKzYz/qVAnzbL+CG5cVY06OT9qYD5i8JakQoPfSuKFcIypF/uwsfPOWpXLqHNEWVcih90xVnEm31w+3o6Y4L2FrdmlJHpaW5GFddUAqDIvzc1BXFcDfPng1ZmUpid8zCiiEQIFVAMQa+6iQn6CAYUjJpo9f0TwNhQBrHAT3heCI9AQ+JWrRA1plfavNJjVGnYjks7wcv+X92Vnu6332nwvi73/bhC9v/j1+dcA5Kf7xqV68d7zL9vOJCEVD+wBujYP6OxmWmRNSpgAopbdRSldK/r0BAISQrwO4G8BXqVs/LEHIiif6hsbjqi6djIASO/jJfmy1UJY+ElZxumcInXFsAUeIlrgT8danHRgen7qUNAHw0/tX4dH1lfjq+ko8dl01dp/qxaWxCC4Oy6uRJzsBLL1jdIomACwvt9I+FQLsPNFtez6FABcG4o/F72nuMyp9FZiFVyy0XBxGc+8QPr0wgCtKrd9jfWDqqgLY8vgGXM15NgTANRUFJsVhJ8yzfASP37AYit49za9o3/3p/avwk/tWGUJb9vWOgVHXv8enENx99XwQwpLdztXHVELb/NyVZZgQLJPOwfi2I1GIxp/nm9D5FOCx6xfHpayZB+8m4RtLom0/0IoPmuznn4iG9kHUx7h/iSAtISBCyJ0A/jOAL1JKU55p2qgLAh618+ciJ0uBj2hhnVgPLN5kDo8cf2xLY3QiMuncAKHAqZ4hy/uxQl3xWqzF+dmOCpH9DpY0S/YmFjKUzTXH3isDuUbMdMehNouCoSpwxoZqSaBVioYTcPdy/IoRO6cABkfdt98IRyhUCkyEVQSHrQl0sangce6+ZvkVbLqnFi89eR1uv6oM11QU4Cf3rcJP71+FayoKTM84QoHm3iFEVO16lAJ31Jbj0fWVCA6PQ6Xy0AX7rh0ItGSl5sEQUErx5uF2w7OMycohwGwhFxUcHrelxIqwOz2lmhHw9BdXajk4AD5CsLw8H7etKHN17mQjomr1LG6hqhR7U+AFpCsH8N8B5AN4jxBymBCyOZUXk/XQzp+dhS2Pb8Bf3rEcn3MxCToHR7GuOoCFMZJ8spCNm31VuwbHQCZrBttI5ViKpTQ/B7dfJb8HBMB91y4w4pU+BQgOTTha7CoFNr3RYFgsd62cH2Pg9hAtaIUw7rxiYrycvzhs+vktF4cN/rr48wl0LrzNNe+9dgGWlcXf9gAAQqNhwwOgAC7021vMvFD2+wiy/JpBkuVXDO+FRw/nDe5t7jMUDYGWz2AK7+NTPTh6YQBPv9mAxvYBlM2dZRK+EZVi5/Eu4/cThRjX29fc5zoUSKB5Dyy+TgFQSrGsLB+UasolHr8zQrXkPA87qzcem4UCeO6jMwgOjxuKbSJCseNQG568aUnShGA8Y1L05+z22nZzYrJISx0ApXTpVF2rviWIlw+et7wfyM1GXVXA4PN+cKLLMSHUOTiGzsExZ8sXwJzZWUCMsE1uts8y0an+fTf8ZTsUzs6ShlgWFMzChf5RW4G3oaYIy8rysfNYl9VSBvDmkXbcsLQYFFrc9L1j9vFNhohusdRVBYzqxXcaOlCUl423Pu2IyTBRiJY7+Mn9q9DUGcI7DR24a+V8LC/Px97mPgRys/FDrgAnQq339Z/eP4nl5fmoFRKzt19Vhv7hcey34dT/+tMO3HN1YkqrKzSGm5eX4t0Y92htdQBH2gYwrk+6u1fNxx9urMbe5j5sqCmSWnv8HdtQU4Rsv4KJsAqfHrKpbwlib3MfxsMaGWE8QrFlX6vUw+NvP6XAe42d+PNthxwVlgyEEJzmvM6ICpzqChljo3EqAadx8rCbPXZFbvvO9uHTtn7T91+pb8PKBQVYXp7vaovWqnm5MY05n2ItAr3v2gXIzfFj+4HziKgUCgF+ct8qtPYN4VcHz2NgeMJY/2Vzc9BzacwIMykEuG1FGZ68acllywJKKXhLiYEAFhdbUay3Yl5uliU2H0s2Z7mIp4iJUIVE6aGTScQGbeLrKxYUOPKdc3P82t7JNrNBpcDuU734+FQvPjjRBcmtssDHWZWA5oJvqCnCusVFuOXK2Mmv65cWY8vjGwwF8q9/st5UBv9hU7dFsIgLr3NwDF/e/Ht82NRtCEECLZm61MHCj6gUbxxpTyjvs7QkD/2S8A0Pv0JwRVm+IfwBLYH83EdnjNYSbC9rHqWcd8nyAA+vqwQIwbb9rfjqi3sRyM1Gtt+cIBZbIvgUcyv0iEqxeXdz3MKfWdIixsIqtjy+AbeuKIsp/H26opfda4U416PwxwHa77p9RRm+ur4SV803P9+BkTA6B8dM6yscVvHD147aCn9xzbRcHIYCrZWMWFcEmBPiPBouDBgNKTWviaC1bwibdzejfzjqTWf5CP78c1eYev4TANcsKkyJ8AdmQCsIZimNTURpbxSaB8Cwt7lP2huocl4uvry20nWpNwXQ5KJfSv6sLADRxaZSLdSRyEbxPkVLDh69oDW3UwCsqihAQ/sAIqrmot+yvBQfneyBnfqKsieI7TEMYRXInyVnzzDvhRDg8RsWG5O2viWIr7641/QMYoGFjZ7ddRqhkQk0dgyidv5c/HLPOdvzyJL9EQrsPN4FRSFQdVretn2tWFsdgN9HbOP8rF0EValFiM3OUjAyYb0WAXDfZyrww9ed58vjNyzG7bXl2Lqv1fQ73j3Whd2nerDl8Q0AgFtXlOH9411QqSYcHhAKgeqqApqBo9OPWe5gy+Mb8MyvG3GEq3S+7aoylOTnGJWs7zV2YvNu5xoOEbKKV5nFvbg4D3ub+9A9aFYo4n2bl5eNv7pjueHVhUYmsEf3YE52XwKlFJRSXF1RgIGRCZyzKUyLMr+0Z53tV/D1jdU41iFfT8zTJsS5HYhsbqgAzvUOYf3ieZbPmCIST3m6ZwjPfXQGYZUaSuI3jZ2W70coxdO/bjRkEQsTpSL0Y4w5ZWfOEDBL6YZlxYagI9A2VmC83tCIPKZ9uG0AH0oy9U7MDtl84jdiIQB6hqwhosSEP8E3bqjBpntqka3Hj7OzFDy8thI+RbOqFEVBY/sAwg7UhdoFBdjb3Od6f4RZkmqkbB8xcg2UAr/4/TkjfsvCEnZnz/IRi6e1q6kbj7ywF3/32yZs3t2Mj0/1YvPuZoxKhP+cHB8qHHIzKgWqi6LsKAqN1kcpddyZaXFRLh5ZX2lp9SwT/oCmGBvb7bvMAtrzz5+dhbqqgHRTlomwih2H2vDI83vw3rGuaE8bG5YCM3D43EFdVUCbE7r1nO0jeOqmJaZKVje7tIldPesqA5Zc0TULreSItz7twD+824SGC+aktZiovzg0jmfeagQAfPOWpfjuF1bgjT+7AXdfs8DIIYRV4EjbgC31EzCbLCoFxiZUR4bXleX5+Mo6rVjLrrBu7iwrbZQholI0XBiwVBp/ZV2lbVK5a3DU9JwqJRTjiApjnSgwe8GpwmWvAABNCdTOn2vyALYfOI9/eLcJX31xL3Y6cHYPnrtoSaIumpeLr9p05ZOFWviNWCi0DpmJQAzRMNe9qTNkJLS3PL4BweFxo9AnElFxsivkKJQa2gewoaYIOQK33KcQaWuBaysDxu/0KVqHyofWLALfi2wirBpxbCak7BCOUIt1d7x9wBQiccLweASdg87hC1n9QUS1F+aAxpTZcajNVXsFABgcDWP7wfOOSXeFAO39I6hvCeI7d60wCSBm8fWGxiztJMKqxmYCNI+KGS/MwGHPngmLuqoAtj2xEX/1+eXY9sRGixA5HyOWfftVZSgQdpC70D+CaxcVmubIEUk/JcYuEueczIIf5+YJA5svYhjLLShgYg5l+xXcd+0ClM/NgU8haOoKYcehNiwvz8f2J6+Ttq8W9z7gx0IBHL0wAAIKn0IMUsKDqyvw5E1LpKHUjTVFpuc0EGNzJr9fwbduuyKlwh+YIQqgviWIFz85a3qPTdKJsOpIk7k4PGFxcYvzc/DA6gppbHJpyRwU2hQ7Afp+BAnsd6qNWf7+Ow0dqKsK4Ju3LDXix8za8CkkJn+4NzSGvc192HR3LR5ZXwm/j/HAKfaevWg5fklxHn71xEb89eeX46Unr8Pf3L8KD6yuME18ny+aA2BCyo5KK7v9sdY7b3wxKuPtV5UhX7KZj0KApWX5jm0UZGDWZG6OH0/dWIPyuTkxi8MiEYoiQaCwdhQ+AihKNF4PANufvA5//fnl+On9q/BtXTiU2BT/9YbGjHAaM16YEmDPnkdTZ8joVirivmsXSq9BoAmzp25aYvF8OgZGLd6y7DnFk8ZSaTQcyxQbAEudgxvYreKbryjBz7/yGfzhxmrDs+ANlL4hc86mZE42Vlc6C1425x5euwhfWVdp7tMjuQHM62PPaUJYzLOyFFOeirG6eGWfClz2OQBAC0HIkjMEmsV125WlJouhcHaW7faJeu0M3mvslE50u1AOYwdQaBWrsRA7Gh/FXSvnGwwQFgLY8vgG7G3uQ3v/CLY4tPP1KcCHJ3u0ODkBiufkIBJhsUp5i+EXPzmL22vL8c1bBDIX12NbXIx1VQE8vLYSR9rM8XECrYXC20c7DKs3y0dQu6DAlJScnaVon1MKQsxMC/Ycb1leip0S9g1r3/vg6grsONSGXp2l9eHJHo2pYnt3mLfYih/fuwr9IxMYD4/BR7Q2yf0jExbL1Ocj6BOe7/KyfBTn52B2lg879Zg+C/UsLJxtPDMe2w+et8Sgi/NzTCwfJsTE79a3BLH5ozMGW4t1wuWT6N/9wgrsPtVjiZOvrQ6gMDcbrx5qQ2GuWQGoFGjsGHQ9Nyncsdoa2gcMxTYeVpHtV7Dl8Q1YXJwn9TB4lMzJRo9+v+0uw7po8swpFi579VCbZXw9l8bxwidnjd9JoK1fVaUguhJXVYosv4KVCwrwzFuNGNef54OrKyy1EtmSOH6W4CbkZvmgUs1jZ6yun719HC9+chYqpcY9SbZHMCMUQCA3Wzo5rq4owKZ7arG3uc/0sKuKctFvM/FUCkvyLhbWVgdw5Hw/4tlnzO35mdUgLh6e4vrywfOGcOUTVQTmBLJKo1WWBJogjlBq8TwiKsXPd540uah7m/tMhS3hCLUIp+DwuFR4/KaxE09/cSUa2geMJCWgCWgWBhoLq/ArBA+trUR+jt+UwFxbHcB37lqBvc19lnPfcZWZQsePZ+u+Vvzik2YTjVGGiAr84ndnDcELwES3ZYn3lXo8XFS4J7svoakrBL9C4PcpiERUEKKFIVVVCyM8c+9KQ0DXVQXw8JpFpvMo0J7VrqZuEEKggBpCbOu+VhNN9qsv7sWoENp6p6HDpADqW4LSAqtYmwfVzp+LfWcvavRTn27UOBw/v2AWrlpQgA9OdENV5Q0GmQfK7u/YhCZMD5/vtz8xNI9qcXGeoQCkxygwrHPeMGJKl4XVRERUCr9ezJalJ5WZMCbQ4v0PrK7A3uY+g5QwPqGiW6CAr9Pnpii4RWMoODyBLB/B51aU4f0TXZY5NG6j7CeLGRECYoJHxMNrK1FXFTApCAotXucULohH+PsIcEVZviPbgHkVbkAAS2x0+4HWKPdbiKmyWPDtV5VhaUke1lQF4Pdp4aEcPWEs64pZVZSLbU9sNFHSGCg0q/LLm3+Pn7193GDq8L/Q7zPTQOtbgrjQP2IUv7DfTBFlrzy4ugIL9N2k6qoC2PaNDfjssmLDigxHKFovDlsK+1iLgQ01RUa5P4HGv76G246Sd6frW4J45q1GW+EvPo+R8bBtB02/XoXLQmHi3FH1cGNEpfhSXQUeXlcJiqhHGFapqXAO0JQg37oABNj0ZgPeO9aFiKp5QZvurkVTZwjff+0oPj7Vi++/dhTPfXRGmjsRi/FEejSjJ8ay1gfHwgbtRgFw7zVWOiSP2dl+EMD4rTIU5+dgQ02RYZxQAC8fPI9rbbYSrZqXa8yJw+f7HZvIRVRIQ2AMD0qeF6BZ7Y/fsBjXLS3GprtrkT87CyqNPscFhbMtskNFlGEE/e+lZflSob28PN+UL2TsoOaeS9JQr0LI5VMINtVgCU6ePqiQaC2AWGK/83gXiB7OkFG7nFxg8TNF0cIZfNHOQ2sWofHCAE50hbCueh7+4rYrsONQG7bF8Cz8isZmOSMIrWy/Yiq3D+Rmm0JCAPCR3p0QPUPw+4hhwQCQVs4wzrksbMMQocDm3c0WBUYAPMS1r+VpoArRKIk3Ly/VXOcJFYQQhEYmpF7Mt267AgfOXcT4hAoVwO9O92Kf7iIbC49qlZ5P3rTECEMRAvz70Q5EjrQj269g0921eOatRoxNaM/g1itLpYJyRbkWrinKy8brh6MtxLtDY0YRj6IQRCJRemgkoqKpM2R4XQ+vrTQsOBY+YJbkg7rVKDKuVGr1mCqL8gwrXaWAypnmqkoRHB7HuwKdkLFNJsKal1G7oAAPr620bCUoFpI9tGYRahcUYNMbRx0LIntDYwadMRyheOvTDsMoKczLxkUhnr5ywVzTfRTB97nnb8lEhCIvx4/7rjWHBwGgNThsTNmJCNV6GjmsHOb91LcE8cgLe40Q0LZvbDCeF+/VLy2dg8euX2xQMvedvYindaYdHz4CNBnClJFCNGXmV4jRCfSV+jZj9zAeYliahTFlyWG/7iGmIiE8IxQAc/12HGrDK/VtiETMD1GcOrxVqFLNil9XHcBYWMXGmiLkz85CIDcb7zR0WHcaE87FFirvegIao2M8rGL/OS3J+uDqCrx0oNV28RFovOUzPUNmBQPNwzh4LmjQxxraB4zJm+VX8FBdhalgJxyh+N3pXtQuKEBweFxazHPgnCa0tzy+AXdcVWZUtsqWmmg1MkHHsLe5zwhJRCjwwYlu3Ly8FKsWFuBgSxARlRruNR8fZ/dry+Mb8POdJ/G7072GBZY/y49Brhaha3DUVM/BC8yJsIp3GjoMAyCsUrx/vAt+n2LJATR1hdCkV7I+dWMNGjsGjdg9O+6hNYtw/uIwPjnVazQI+9HrR9HYPoAHVlfgAT3XwO7/prtrERweN8X6c7IU454wYgCbG0xQObGgsvwKArnZaBCMl4fXVhq8ellugYEPhwRysw1jSFEUMDqXNueiz1dBVMBNRKj+GTVi/XeuLMcrXLjRpwBne53Da0yS75W0oDjZFdKqpSPRvBKLvxsg9hRZhtr5c7F1Xyv+6f2Txj0dD6t49VAb6qoCeGB1hSlMarRz5o5tbB/A1zdW4zeNnbiztty4r2Je4UGdHMIUCpvL4nOwU8A/FGqOWK1EsvcCZpgRCgCAYZ2tXFBgxEvZQ5GFOXioFDjU2g+VUjR1hQzr9DfCXsMimFZnC5HF5H++86QlkceuwyBLnoUlm7T49IRplj6Z/IxGyE3entCYZUvFc33D+P5rR3HftQtsmRxsbE/etAS7T/UYk1XMC/BKgWcwMASEZGJYpdqGGtyFI3osnFHrXqlvQzgS9QaYJ8DGMDRuDgMxwSfb69enENy1cj72nImGPVQKPFxXgcYLA6ZEI/vqeFhF/uws/OufrEd9S9D4/bxy488X0XNDOw61YcvjGyyxZh51VQFsursWm95oMH73prtrjeNe1Y0DGa7Rcw0PrK7Q+tJzh62rDhiWrhuw6zHPSyEEKhfeuX5ZMe5aOR/PvBU1JlYuKMArpA0ABVEIfIQYBpUo/EBjM94YvZVVyvJzIjg8YboPV1do3szTbzaYclo+RUE4rNpWHb/T0Clt4cAuVVcVwENrFhnjjkRUdAm04pNdIWM7zs27m1FZlGcoWpmCZwpF9AJ4z1ycI8/uOm1Zi/3DWq3E8nJ5KGmymDEKAIAR9x0Pqzhw7qJxU2WdF4FoLI8tDFFgi9Y/D7+P4OE1i4z9Q9n1+YpYsdKPVSwTojFjhsYj+OBENyjVElIgBGHdYmUTJRyh2NXUHQ3jSKyh4vwcbHtiI3YcasNvGjpNbvrh8/3RCl5oi+x4Z8jiJTHaKwsb7TjUBgItxNDQPmB4VnxfGva7Zclf0elgidDg8Dja+0ewbX+r6X5/85alJmbTtv3RJNnS0jnGs3zm3pWGYDUuoXd+fPyGxUbymEJT/L+zeYZ8zJUJbN5wqG8J4pYrS03JTV5pMrpffUsQ33/tqHHv+HvCrGeVUtMctLNns3wEm+6JKgoxgbm0LN/EpvHrliV7ZjKFxCdfQbVwClPClfNysbw83+QpvNPQYdSYUJXiS+sWWZhMvPfDe6d2IPo9/vF9q/DD148alc81xXmmRHXtwgI8ur4Sje0DJiVz0xUlGJ2IoCgvG78+0m6ZWzLhny1ssci8gImI9tsfXluJ4x0NmIhQzXgSFPL2A61o6gpZQpbQfwvfDyoSicoMMczJM+k21BQhy69gXA/fgcKR7ZUMzCgFYEehEy1UQBPOX1lXiYWFsxHIzTZZQaxZl9OkpipF7YIC6WJjoZrrlxabmDRf31iNzbubQanWG+an96/CUzctMRZuU2cIP3qjweQCU8BQEiyRVJyfg2wfMSbvygUFBv1xYMSs7O6sLccv95yLhivuqTXGygQgP2mZMBEXPaNYvlLfhm37o5Ywq0vgQx4iCIGJBbN1XysUfQXwSoj3onYcajMU6ZnuS0a4ip3j+d1n0NI3bFh0e5v7kC9UZDe0D0iFg0KsrSx4wwGAkU9QiMZCOtw2YFGa9S1BPPL8HsNafbm+zYg7hzgKKc+FZ/fyV/taLRbtQ8K2gA+ursArutDK0gWarBncrw60alZyxCysWGLe79MsaEUhuG5JEU52hdBzadx4jl/fWI09zX041jFoeKHMeFmphxEZRKYNEFUIhMAS4uTn1PLy/Gh4CUBNsXlvC+ap1y4o0GiZVBPWH53sMX7bj+9bhX98r8mRGcTCKhaByqSubjBse2Kjae3xubCyubMM9tz4hGpixdW3BE0dBHy+qMzgmU4ikw6ARnOGXjOiM8ZS2Q5ixigAfrKLN5VP5ADRmCyfvGntGzLF/0RmQVFuFirm5eLTtgFpXJgv0GLClp8we5v7LOXrv/ikGfevrjAErebyW9UOVaOWG3PFWairdv5ck8vM4/aryvDdL6zA7bXlFuuQ/f/ZXadNSnPHoTYjRCFSTvm+NDxtjQmFb26pl27mwTstTNjyoRE2Dj6UtunuWvzg9aMA1ZTg2ETUymLCWfSy+GdGAaMegAfz+n655xxu15+1aDjw+YSIHh5k3gt/D/c295nyK7z3yBcmMmXEfmNTZ8gi/H0Elk3BGcNLfHZ+n2KyWPl6Dn4M0dCP1hI6rFLLnsJjE6qlZ5ACYNXCApTNnYWnf91oUSz87w/kZuOB1RXoDY2huXfIZNFfo9OwxftFoSV39zT3mRKsweFxy/y4eXmpqbYiODyOBYWzLQqAX9+ysAqbuxRAOGL24ti9BmCsqcGxsKE4GTnhwLmLhvKTtepu6oxW5FMAn5yKfsdYP8yb5DwbPlydbMwIBSC6xYwBwydymOXhU7R4Mv/51n2txiJg8T8xrDEwGsa311aisb3BNi5sF0qwa5R2umcI//Buk7G4ZB4Hgdb/h49DAjAs1j1nrEVwBJqLXZqfY4Rq+KpDXpiISosAFoXABFAgN9vWqq2rCuDqikJ02rRJZkwN3kuiVNs6j/0Wv77gS/Jz0BMaMykOCu16Tl7WXkGgFOfnYJbumRAAK+bn44S+SHkPUbwHYj4homrj/On9q0y/aUNNkSn34vMRtPePWBS5QmDKeSyXdCpVbLrM8gKXQZUU7xGYlSGv1CiLpbiE30dwvDNkGDuA+X6Jc9qOo1M2d5YlV8SOo9CUomiw8bx7lXm7AjsnkJttstbvu3YB3m7oNJSiLKziNHcZHl1fadRZsPlYUzoHp7svmc7JN6AkJLqfrygzWMiQrSHWyZXluZhnw4erk40ZoQD4yc645BYw2qeimIQ/oMX7eGw/0IpN99SaEo4sjivGoPm4MBAVzPua+wz6qVOjNOZi7m3uszCFFKJZYiLNj7faQTXqIjNE/T6CW5eXYldTN7buazXCEoB9MZmdSy8ma29cZt6zWWSoPHnTEnzQ1C3tssh46k4KZzxCTWwkEUwByrwsdm6m6P0KMaqD+d/2yPN7jDhwIDfbUIjiPbjlylKjU6eY6GNgFjoLv314sgfb9rfCr2ibv7Cwy61XRq3Y8bCK0rmzAJjvHdsRKpYQ2HGoTcokW1sdwE3LS03KnRc2IMTCiFpaOgdZCjFVty8tnYP1i+dh2/4obZInOwBaEpsP99nN7WKh5YUoICMqcOuVJbh2UaEx7qbOkIl3v3JBgekZ8gqRGVrB4XGEj0SpqOJ4ASsVXHzNIMqSsz1Rj4aFeuqqAqZwLjMaWSiUUZq1BLZ5DX19Y7XBPHtP35+DrX9PASQIJhR4LrnF9RLcP/5miwuyVLdcmLBnpdrs4S8vz5dSTsX47NZ9rfApsUvlVQChkQnsbe7D4zfU4IWPmxHRE0RHLwygqcvszopCcNPdtaYq21cPtRmhifGwir995zhuWl5q22JAtDLFZCz7TrNA+ROFdF1VANuFuCpbpHwVrChsX663smLEW8b2HxDpjUzxGuPn4rzib6tvCRqfq4AlvPHNW5aavEn+90Uk84Y//7O7ThtCPqxSfG5FqeHeAzAUm0q1Xku7FJgUvds4sEwx+n3EUo0q3mf2LIrysnG2dwjHOgbR3HPJcr7brizF7bXlprg+qzVg1r9oMNnhdFcIP3jtqClEKrLVdjV1m/ZB4JUE0V+Lz5D9Jjan6luCBuXXJyFnANb5JFuSLIysKAQ0whL42md8qGfrvlb8274W03eZh8vPzYb2ATReGDDlEhgdWiEwKTqZR5IMzAgFwCY7zyUXk8BO7t9TNy3BLn3HML+ivQaiLqFd/ByAhf3B702gJW3lYxb7rEcnBjFNTjuB7URDFNkj+88FsboyIC10sbuffDKWWZFne6PWkKx/Pf9dhuDwOJYLnRctQlnCbPIpej/3CIWPAD/mCmXY/0WPhlf0MoHN6gi0zykioCYPTswHKICpXQCf/BXvvTjHPtQ9oT1n+nDrlaWGUFOg9dvhBYtIFnCCuPsZK2qSfZd/jswzZbkvJpRENHYM4rtfWGGqq+GNEK3AKeYwAWjzbv+5oCk5vu2Jjaa9DMK6ocTCqGKYiF+rsn5C7HczxhUAi/AH5Al1HuzcdkQGRfcot+5rle4fwhQ9u+db97XiJX2HMJar4tmG4r0XvelkYUYoAIZF83Jtk8D8AhRpoXVVAWx/8jqpQJXFYMWJ+MDqCkMoMGucp01C51ID0VCNKPMMWqNO1QOlhhspE9iycTE8uLrCUnXc2DHoqDTE3ydymY+c7zdtg3jL8lLTOUSh6LRY+ePb+0dMCTUGHyF47PrFaOwYNHkQDDLGl6wZGA9RuGT5ok2/2LFiAQ/LSTChYve7RMuVeWCsKC3LH52XtfPnGhRjCsSVBBQJDWe6L8XkkYv3iiIaHgLMlN3a+XMByDejYfeYt+L9PgJFn98sXCaCV7B1VVplbmO7fRjVbq3asfxePdRmhB3DEWoUgPGwS6jz5x5zaB1eMkfbYvbnO09aPrv9qjKTN7LjUJuxPST7PdcvNddciPddngGaPGaEAnCTBOYLqWTWr5NAFSFORBlzRow9723WdkSS7dJEYG4nwBK+rIIzlsCW/ZYnb6wxXYsJGZkyiyW4v3nLUovVQ4VziNx0ALYhJ/H4aAO1qIUUjkSrh2VJMpmwj+UZiVbWzctLTfFndu9465ftQMVz7WW/i8V/o8V00e0DKbTwAaPW7m3uczRInCDugCd6MHbfiYZIzCyysYmIae/k0FjYmBOB3GzT99h9YnkP5v2y+/Laf1yQNqDj+0Zt3deKTW9oRAo27yF4WOw+ikpcVM5s3wW34R0n42dDTZHWAsQmXsvaa9+1cr6pPihL34yHXYOnBTP4FGJ4ePye13w1v8ybTgZmhALgFyXfyMkESSEVP9FjCVq+I2Ms5oxIMQM0wSJzHVnC6ul7rNWGk8HtteU40zuE7sFRaa8Y9vtlYRSZgBMT1B+d7DEYRrLcR5aPSL0xwPq8HtaLjfh6DCIpznOKcdvlM8R7zaM0P8fa8hr21q+MMbShpsjk/bFn2NQZwo/0oieWkObHZSfkYqGuSmOabT/Qiob2AaiqdX9mKbj539QZMkJC4j3pDo2ZlLMqfh/ye8xCreIcJ4jWN9S3BA3hD8DweEWDTWTS8dfglTNfx8BDrPznWUtiZ1bpPdKxrCQPPUPjuPmKEnz3CysARFtubz/QirK5s0ydaHccarMIf7HPD3/v3LT0mCxmhAKI5frz/NuIzrYAYKnatevJzcf9Pj7VizuuKjOy+axFL18dabcYZW4eWwTJpIGJgl2MwTPEE0ZhrjtfTi8KRXPuIyrYxQku668iLgpZcZ6IeLw2Zilm6fsEx7K6nO6DmMTm7/Wmu2uN8ft8CtSwdUOiWJ5KrN/B10EAiLkxuzj/32noMOU5fHpIyacA3YOjxrlZGItfN05j5YXjsY5BI7z2IOcliBa2SmEy2Oyq+RlkynlPs3nHMZHhw1NLWWdW8bx7m629is72DUOl2v6+fOX7o+vlBpVYd7IwMBs3X1Fiu/7imb+JYkYogFgLSraYeT454FyS/Y7QE4jFwhUCg23kZkE/sLpCuhGIm8UVD0wWuUN4IN4witgETWylwPf/YaEGu+Sk3TVSYSEx15zVgTy6vlKaKExkjDwld3xCNVhjLOFnl5BOFOK8BbRwmawhGYOszoH1XeJZZK/Ut5n6JtnlSZzCKUw4yo6RVYyLiV47D9Tp9+QI/YjEnv0baooslG7xvGJug5BoXs5pYx7+N4q0186BUUvV/FRjRigAwMx4EIudmIAS3UqeOipLtrIHzCftePAcftmWfbIxbhd446koBXdT9MLGE08YxUkoBofHTbxxsWGc7Fxu7pfsGCcBJPuMd83DqiZ03CxGN2PkhREftuL77sjmlaxtsRuI3haDzLt0akzGK1dAM3JYcpJBgRa+4b24WMl9p3vH5s/Pd540Oq2KOZBY3jx/Hjb+Vw+1mfIYpYIgrquSU7rFY7Y9sRGbPzpj6v/EuPws3yDLY/F5P8Y00lIb9iHMqcKMUQCAPU3Mzq0U2+XyQoNv2atRuOR8fjccXlEosWv87O3jlvazTt9zi1isJx6yhWp3XafxyMI68cLN73USQHaficIxmYwLcR7xYStZF0nA3A10PKzimV83mlomuLkevyUkYKWH2jWNE5W8Hf2RQKtAF724WBZ6LOW841AbZmf5TKwoXhgnEh57YHUFXq6Peqay0J4dpVu8t9cuKsT7emtw1hLjeGfIYsnL7sM3b1lqMI3sQpiJrulEkVYFQAj5NoC/B1BCKbVvrZkk2E1Ou/edLDx+kfKCXwFQOjfH1PPGicNrZ+2J7ScAbR9XNkFCIxMJ7xfKM1Li9S6clKiT5ccn6BIRsG4tSycBZPeZGwExGYjzSCx+43/j3uY+9AghiiNtA0azO7dK4NpFhdh5rMtWycsS87JQBDtOxNVCHx8GJwvdybMRGTIKAW5bYd7Ok/99TvdBNle2fSO20ojXo8vyK1i5MLqlaiwyAA++y2q83lMykTYFQAhZBOAOAO7KBpMAu4fixq0UNbOMDMasosp5uSYFwAs88TyitffcR2dwzaJCy05Pz3+slZOLCT4g2gjN7WRJxIpiiFeJimC/N964p9vzOz1Lp8StGwExWTglMK3UV2LKBcUbJoil5GPRRXkGnCyk1DU4auyCxsNpbolznefj7202N85TKfDesS7cLNSTuIEdecEJbi1vPlxclJeNhgsDWgt01dq5NpaAZ40OncZ92SoAAP8I4D8DeGOqLugU03YSiHbxPJ72yPr/1y4owNNvNhjf5StiZecRFcn7J7rx3rEui5VMKQx2hozXLIaZYk1oN9aODJNRopOZ4G7Oz36XU3LWTeI2VRB/P99IT6S+3rqiDN2Doya2TLx00Fib0tjtkicTVA3tA0a/fADoHBwzmG+iJ2N3L534+GKSlX3+IwkjJxbEuRLIzY5ZdOjW8mZKXAyJ3XGV1VsR74PJ6+IIAeyabud4MpEWBUAIuRfABUrpEUJEUWc59gkATwBAZeXkt0Wzm5xOAsAunrf9yetMBS91VVrPF75yle/hLjuPmBjiNxdh5wA0XjhjZ8jK0XkX382ETjTWmKgSBdwL8Xiua3css7hkCf90JNsAa6ES3wRs0921ps9YJ0hZ4aJbNHWGDCveac6LzdTEFuDB4XH89P5VeHB1Bb790mGc64s2U2Q9btzAqd1CXVUAT39xJX7w2lFzq5MEGHDiXIlleMRjmNiFxEYmInGFj2R1LPymR9M+B0AI2QmgXPLRDwB8H1r4JyYopc8DeB4A1qxZIy/DSzGcQgexYqD8JLejVbLEUGhkAs9/3GyqN2Eb0zABsLw837IIFQKTMHWTiJtMrNFJibLr86/5zyczweMR3umIp8YC//vFRnr8vtH8Z7aFizEg1qYAVkudHxd//g01fHt0885oT9y4xFTMxXrcuP39Tu0WZGQEmaEgq06PxVRzMjziMUzsWFZu7gP//O2SwFNtoKRMAVBKb5O9TwhZBWAxAGb9VwA4RAhZRyntlH0n3XAjuJwodbHOw/7/1Rf3WphEFOZCmLqqAO6sLTe1cfjiNQssCzhW4VsqYo2sjN9tYjpeLySe49MRT3UD3jsRayacPosXshbmcW0sLnRNZWDnsEtkx4KTgGN5i3Hdw62RNLKThaf4RnayOec2HObWu2THhkYmbHtRufn9U1HpGwtTHgKilB4FUMpeE0LOAVgzFSygVEFmbTLXE4BFCcgeNu9aEkTXnYyTnD87y9RYbFmZtZtmvIVvk4VYxi8WmDGh9tLB84joBVeKZJtCp/PHY9FvqCmy9KnJJMTKR9i1O3CLMqGFufbaDDuFurfZuWuqXaWrGzgpcTEv0dxjbWQnKna+atlJ0ceyrOOxvJNlpaczHMkwo+oAnBCLn+wkfGTJPdm2ifx1xNoCUSjbccQBdzTOWJOL3+Q93kkou1d7m/uEXa6ioQMZlzysAhC2KXQaR0IWvaS/UybB7hnFanfgBk/etATvn+hCRNVaODypNyTjr2E3p+0MhHh6Y8ngRonXVdn3WZKNTaxaTqWin2qO/lRcL+0KgFJane4xxCvgReEjTkoCeadLcas8sb9QIm7oZBdhvJx3vm1Clo9g2xMbDQXG3HdFaHC1t1neSlfW6dEO8Xote5ut/Z3SbW25RbLCVz5Fgaqq8CmK5TMnRpLMA4k1d5P5u5yetWzui1XLz+46nbCSssNU55Sm6nppVwCZgHgFvCh8xEkJQBrDZdfh+wvx271NhUs4WeHCt00Y53rMOCmlDTXmXiuAJkB+LNlI3Q7xKr10UOqShWSMnQ/jyHa5i8VIEj0Q2dxNB5WXfS4LqyZDSdlBpHD+fOdJY7tJN3nBeK8/VTksTwEgfgHvJsYoO17GIEhku7fJWAeTFS7E4bWdAqurCuDxGxbjud3Nxu/2KSTu0Ea8cdqpptQlC8kYe6x+T/w1REaSLK7O5o1Tbyw3v8ttbiMRYygZSsoO4u//5FQvPj7Va6tkJmvBT5UB4ykAJCbg3ZzTLsllanZF4tvwA5icdTBZ4eLUNsGpR9Av95ybNL87XmRCki1RTHbsbvo98ZYz77HK4ur8vJlMDmCyuQ0nJENJ2YFfu2xbWcBeyUzWgp8qA8ZTAFOMuqoAvnXbFZNKXE3WOpiMcKmrkrdNcLJ4xOIZAhgVmmKhlofkwA1RgCFWXJ09m8kqJbdCMRlFisnOAbDzs7UbS8kkw4KfCgOG0AxlSMiwZs0aevDgwaSfNx1FQ5PN8E81IyEWnt11Gn//2ybD4vz255cbu2mx+8vizQ/pLTNi8bc9TA6ZNkf4eZDlUKGeaQV8ItyyoTLp/hNC6imla8T3PQ8A6Skamqx2z7TwBr+hupjXkFmYYruBVNzzTFqA6UCmzRE3YY1MLeDj4fa+Ztr9l8FTAJjejJFMQXB43NgTQZbXEBdDqu+5HV3VQ3oRSyh6a3Fq4SkATG/GSKYg3oUbDyMkEdjRVT1kNry1OLXwFICO6eCuZTLiXbipZoQ40VU9ZDa8tTh18BSAh6TBzcJlcfn2/pGUxnpTvcvX5YSZniuZyfAUgIcpg3XXK/m+r8mAHV3VgxnTgXXjIXXwFICHKQPP8IioFA+vW4SFhbNTJqC9UEJsTAfWjYfUwVMA0xjTzXUXE8UPJtCJ1ENy4bFuZja8QrBpiunquk83pTUT4D2Tyx9eIdgUI9aimuyim66uOx+W8QRPZsALlc1ceAogBYhlnSfDes/0Ha9iYbp6MB48XE6w7hThYdKQWefxfO4aGb7jlROSdg88eMhQ1LcE8eyu06hvCaZ7KLbwPIAUIFZiLWkbfkzTHa8AL/no4fLGdPFwPQWQAsSqik1Guft0F6Beyb+HyxnTJUfnKYAUIVZiLRndQKejABUTv9Nl3B48xIPpYqB5NFAPU4bp4hZ78JAMZBLLzaOBekg7potb7MFDMhDLw80EBeEpAA9ThuniFnvwkGpkijfsKQAPU4bpmrfw4CHZyBRv2FMAHqYUXuLXg4fM8YY9BeDBgwcPU4xM8YbTpgAIIf8JwDcBRAD8O6X0P6drLB48ePAw1cgEbzgtCoAQcguAewFcQykdI4SUpmMcHjx48DCTka5eQH8K4GeU0jEAoJR2p2kcHjx48DBjkS4FcAWAzxJC9hFCPiKErLU7kBDyBCHkICHkYE9PzxQO0YMHDx4ub6QsBEQI2QmgXPLRD/TrzgOwAcBaAC8RQmqopCyZUvo8gOcBrRI4VeP14MGDh5mGlCkASultdp8RQv4UwKu6wN9PCFEBFAPwTHwPHjx4mCKkKwT0OoBbAIAQcgWAbAC9aRqLBw8ePMxIpKUZHCEkG8AvAFwLYBzAX1FKP3DxvR4ALQlethiZr2QyfYyZPj7AG2MykOnjAzJ/jJk2vipKaYn45rTqBjoZEEIOyrrhZRIyfYyZPj7AG2MykOnjAzJ/jJk+PgZvS0gPHjx4mKHwFIAHDx48zFDMJAXwfLoH4AKZPsZMHx/gjTEZyPTxAZk/xkwfH4AZlAPw4MGDBw9mzCQPwIMHDx48cPAUgAcPHjzMVFBKL/t/AO4E0ATgNIDvTsH1zgE4CuAwgIP6e/MAvAfglP7/gP4+AfBP+tg+BbCaO88f6cefAvBH3Pt1+vlP698lLsb0CwDdABq491I+JrtruBzf0wAu6PfxMIAvcJ99T79WE4DPx3rWABYD2Ke/vx1Atv5+jv76tP55tc34FgHYBeAYgEYAf5GB99BujJl0H2cB2A/giD7G/5LoeZM1dpfj+yWAs9w9vDZdzzmpsirVF0j3PwA+AGcA1ECrOD4C4KoUX/McgGLhvf+XTUYA3wXwt/rfXwDwjj6RNgDYx02GZv3/Af1vJlz268cS/bt3uRjTjQBWwyxgUz4mu2u4HN/T0IoExWOv0p9jjr6oz+jP2fZZA3gJwFf0vzcD+FP97/8TwGb9768A2G4zvvlscQPIB3BSH0cm3UO7MWbSfSQA5uh/Z0ETyBviPW8yx+5yfL8E8CXJ8VP+nJMqq1J9gXT/A7ARwG+5198D8L0UX/McrAqgCcB8/e/5AJr0v58D8Ih4HIBHADzHvf+c/t58ACe4903HxRhXNcwCNuVjsruGy/E9DbngMj1DAL/Vn7P0WesLrReAX5wT7Lv63379ODce1RsAbs+0e2gzxoy8jwByARwCsD7e8yZz7C7H90vIFUDan/Nk/s2EHMBCAOe51236e6kEBfAuIaSeEPKE/l4ZpbRD/7sTQFmM8Tm93yZ5PxFMxZjsruEWf0YI+ZQQ8gtCCNs+Kd7xFQHop5SGJeMzvqN/PqAfbwtCSDWAz0CzDjPyHgpjBDLoPhJCfISQw9BCfu9Bs9jjPW8yx+44Pkopu4d/o9/DfySE5IjjczmOVK6VuDETFEA6cAOldDWAuwB8kxByI/8h1VQ8TcvIbDAVY0rgGv8TwBJoPaM6APxDCoYVFwghcwDsAPAtSukg/1mm3EPJGDPqPlJKI5TSawFUAFgH4Mp0jkeEOD5CyEpoXsSV0NrXzwPwnRSPYUpkxExQABegJccYKvT3UgZK6QX9/90AXoM2ybsIIfMBQP8/2wXNbnxO71dI3k8EUzEmu2vEBKW0S1+MKoAXoN3HRMbXB6CQEOIX3jedS/+8QD/eAkJIFjTBuoVS+mqM35eWeygbY6bdRwZKaT+0pPXGBM6bzLHHGt+dlNIOqmEMwP9C4vcwJWslUcwEBXAAwDJCyGK9C+lXALyZqosRQvIIIfnsbwB3AGjQr/lH+mF/BC0+C/39rxENGwAM6G7gbwHcQQgJ6C77HdBilh0ABgkhGwghBMDXuHPFi6kYk901YoItBh33Q7uP7JxfIYTkEEIWA1gGLbEmfda6NbULwJdsfisb35cAfKAfL46FAPgXAMcppf+V+yhj7qHdGDPsPpYQQgr1v2dDy1EcT+C8yRx7rPGd4AQzAXCfcA/TvlYSRqqTDJnwD1qm/iS0WOMPUnytGmjMA0Yj+4H+fhGA96FRvHYCmKe/TwA8q4/tKIA13Lkeg0YVOw3gj7n310CbgGcA/He4S1pug+b+T0CLO/7JVIzJ7houx/ev+vU/hbY45nPH/0C/VhM4FpTds9afy3593C8DyNHfn6W/Pq1/XmMzvhugueSfgqNTZtg9tBtjJt3HqwH8hz6WBgCbEj1vssbucnwf6PewAcC/IcoUmvLnnMx/XisIDx48eJihmAkhIA8ePHjwIIGnADx48OBhhsJTAB48ePAwQ+EpAA8ePHiYofAUgAcPHjzMUHgKwIMHDx5mKDwF4MGDBw8zFJ4C8OBhEiCErNUbhM3Sq8Ab9d4xHjxkPLxCMA8eJglCyE+gVazOBtBGKf1/0jwkDx5cwVMAHjxMEnrPmQMARgFcRymNpHlIHjy4ghcC8uBh8igCMAfaLlyz0jwWDx5cw/MAPHiYJAghbwL4FbStCedTSv8szUPy4MEV/LEP8eDBgx0IIV8DMEEp3UoI8QH4PSHkVkrpB+kemwcPseB5AB48ePAwQ+HlADx48OBhhsJTAB48ePAwQ+EpAA8ePHiYofAUgAcPHjzMUHgKwIMHDx5mKDwF4MGDBw8zFJ4C8ODBg4cZiv8fESxSRcSYdHwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "GP = GPtideScipy(xd, xd, noise, covfunc, covparams)\n",
    "\n",
    "# Use the .prior() method to obtain some samples\n",
    "yd = GP.prior(samples=1)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(xd, yd,'.')\n",
    "plt.ylabel('some data')\n",
    "plt.xlabel('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a51838a5",
   "metadata": {},
   "source": [
    "## Make a prediction at new points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7f56c41f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 40.8 s, sys: 5.63 s, total: 46.4 s\n",
      "Wall time: 25.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "del GP\n",
    "\n",
    "# Create a new object with the output points\n",
    "GP2 = GPtideScipy(xd, xd, noise, covfunc, covparams)\n",
    "\n",
    "# Predict the mean\n",
    "y_mu = GP2(yd)\n",
    "\n",
    "# plt.figure()\n",
    "# plt.plot(xd, yd,'.')\n",
    "# plt.plot(xo, y_mu,'k--')\n",
    "\n",
    "# plt.legend(('Input data','GP prediction'))\n",
    "# plt.ylabel('some data')\n",
    "# plt.xlabel('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2c8851e",
   "metadata": {},
   "source": [
    "## GPToeplitz test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "85210713",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gptide.gptoeplitz import GPtideToeplitz\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d7c175bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 388 ms, sys: 0 ns, total: 388 ms\n",
      "Wall time: 387 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "GP_toep = GPtideToeplitz(xd, xd, noise, covfunc, covparams)\n",
    "\n",
    "# Predict the mean\n",
    "y_mu2 = GP_toep(yd)\n",
    "\n",
    "# plt.figure()\n",
    "# plt.plot(xd, yd,'.')\n",
    "# plt.plot(xo, y_mu,'k--')\n",
    "# plt.plot(xo, y_mu2,'r--')\n",
    "\n",
    "\n",
    "# plt.legend(('Input data','GP prediction','GP toeplitz'))\n",
    "# plt.ylabel('some data')\n",
    "# plt.xlabel('x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bed0d68f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(55.885714285714286, 65.63307493540051)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "9.78/0.175, 25.4/0.387"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "4a341611",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.22 s, sys: 1.26 s, total: 2.48 s\n",
      "Wall time: 954 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-16553.80165318692"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "GP_toep.log_marg_likelihood(yd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "defa8fa1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 396 ms, sys: 20.5 ms, total: 416 ms\n",
      "Wall time: 415 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-16553.80165318762"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "GP2.log_marg_likelihood(yd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d9ee202a",
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_2049/2047172887.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mGP_toep\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprior\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/code/gptide/gptide/gp.py\u001b[0m in \u001b[0;36mprior\u001b[0;34m(self, samples)\u001b[0m\n\u001b[1;32m     60\u001b[0m         \u001b[0mPlaceholder\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     61\u001b[0m         \"\"\"\n\u001b[0;32m---> 62\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     63\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     64\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mconditional\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myd\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msamples\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "GP_toep.prior(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b2a2f49f",
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_2049/2210190309.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mGP_toep\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconditional\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/code/gptide/gptide/gp.py\u001b[0m in \u001b[0;36mconditional\u001b[0;34m(self, yd, samples)\u001b[0m\n\u001b[1;32m     66\u001b[0m         \u001b[0mPlaceholder\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m         \"\"\"\n\u001b[0;32m---> 68\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     69\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     70\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mlog_marg_likelihood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0myd\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "GP_toep.conditional(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1d663bb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}