07_MTGPs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9f05ccbc-acde-45b9-b9e7-e47d80b9e402",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "plt.rcParams.update({'figure.figsize':(15,12), 'figure.dpi':120})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0641a9c0-1de3-4c4c-bca3-c1b330f0e25a",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_data = pd.read_csv('data.csv', index_col=0)\n",
    "base_data.index = pd.to_datetime(base_data.index)\n",
    "base_data = base_data.asfreq('H')\n",
    "\n",
    "imputed_data = base_data.copy()\n",
    "for idx, row in imputed_data.iterrows():\n",
    "    if np.isnan(row.tmp_f_KDFW):\n",
    "        imputed_data.loc[idx].tmp_f_KDFW = base_data[\n",
    "                                                    (base_data.index.month == idx.month) & \n",
    "                                                    (base_data.index.day == idx.day) & \n",
    "                                                    (base_data.index.hour == idx.hour)\n",
    "                                                ].tmp_f_KDFW.mean()\n",
    "        \n",
    "obs_data = imputed_data[imputed_data.index < '2021-07-01']\n",
    "obs_data = obs_data.asfreq('H')\n",
    "\n",
    "scaler = MinMaxScaler()\n",
    "scaler.fit(obs_data[['tmp_f_KDFW', 'north_load']])\n",
    "imputed_data[['tmp_f_KDFW', 'north_load']] = scaler.transform(imputed_data)\n",
    "\n",
    "imputed_data.fillna(0, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "9c679b8d-a5de-4314-9f81-325fa4ec52c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_X_time_delta = imputed_data[:'2021-06-30'].index.to_pydatetime() - min(imputed_data.index.to_pydatetime())\n",
    "train_x = torch.Tensor([train_X_time_delta[i].total_seconds() for i in range(train_X_time_delta.shape[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "5ba2f636-0207-4c1b-a4b2-14bd57f64020",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_y = torch.Tensor(imputed_data[:'2021-06-30'].north_load)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "cc03a54c-f6bf-4f43-936c-95d0d49b0cd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use the simplest form of GP model, exact inference\n",
    "class ExactGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "# initialize likelihood and model\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = ExactGPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "a6ec2650-a13c-43d3-a25e-8f3ae5abb5c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ExactGPModel(\n",
       "  (likelihood): GaussianLikelihood(\n",
       "    (noise_covar): HomoskedasticNoise(\n",
       "      (raw_noise_constraint): GreaterThan(1.000E-04)\n",
       "    )\n",
       "  )\n",
       "  (mean_module): ConstantMean()\n",
       "  (covar_module): ScaleKernel(\n",
       "    (base_kernel): RBFKernel(\n",
       "      (raw_lengthscale_constraint): Positive()\n",
       "    )\n",
       "    (raw_outputscale_constraint): Positive()\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "a1a07070-bd5d-452b-a4d0-e3a97ec3c47d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussianLikelihood(\n",
       "  (noise_covar): HomoskedasticNoise(\n",
       "    (raw_noise_constraint): GreaterThan(1.000E-04)\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "likelihood.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "e7813f35-11cc-42e7-8941-154aaf9be097",
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "88de7f70-3750-4027-85a5-4742f3004be9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e5c123f-0e6f-4c0b-aff4-5204f6f83cc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "training_iter = 50\n",
    "for i in range(training_iter):\n",
    "    # Zero gradients from previous iteration\n",
    "    optimizer.zero_grad()\n",
    "    # Output from model\n",
    "    output = model(train_x)\n",
    "    # Calc loss and backprop gradients\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (\n",
    "        i + 1, training_iter, loss.item(),\n",
    "        model.covar_module.base_kernel.lengthscale.item(),\n",
    "        model.likelihood.noise.item()\n",
    "    ))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "97ee23e9-8ea0-4ec9-8808-6404ca6cdee8",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'model' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [1]\u001b[0m, in \u001b[0;36m<cell line: 2>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Get into evaluation (predictive posterior) mode\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[43mmodel\u001b[49m\u001b[38;5;241m.\u001b[39meval()\n\u001b[1;32m      3\u001b[0m likelihood\u001b[38;5;241m.\u001b[39meval()\n\u001b[1;32m      5\u001b[0m \u001b[38;5;66;03m# Test points are regularly spaced along [0,1]\u001b[39;00m\n\u001b[1;32m      6\u001b[0m \u001b[38;5;66;03m# Make predictions by feeding model through likelihood\u001b[39;00m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'model' is not defined"
     ]
    }
   ],
   "source": [
    "# Get into evaluation (predictive posterior) mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Test points are regularly spaced along [0,1]\n",
    "# Make predictions by feeding model through likelihood\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    observed_pred = likelihood(model(test_x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f441daf9-9e19-4f8f-8712-382146e36fcd",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/esherma9/miniconda3/envs/NT_Demand/lib/python3.9/site-packages/gpytorch/utils/cholesky.py:38: NumericalWarning: A not p.d., added jitter of 1.0e-06 to the diagonal\n",
      "  warnings.warn(\n",
      "/home/esherma9/miniconda3/envs/NT_Demand/lib/python3.9/site-packages/gpytorch/utils/cholesky.py:38: NumericalWarning: A not p.d., added jitter of 1.0e-05 to the diagonal\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "f_preds = model(test_x)\n",
    "y_preds = likelihood(model(test_x))\n",
    "\n",
    "f_mean = f_preds.mean\n",
    "f_var = f_preds.variance\n",
    "f_covar = f_preds.covariance_matrix\n",
    "f_samples = f_preds.sample(sample_shape=torch.Size((1000,)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "62382eb9-27c7-4252-9109-dfae22d638ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAFICAYAAAAxj6teAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABJ0AAASdAHeZh94AABsjUlEQVR4nO3dd3hTZfvA8e/TdO/dssvee28QVEQQQVEEXxU36E9A3AtUVPTVV1EU90ZxQVFBUGSIshVQNsge3Xu3yfn9cZrQkaRpm7ak3J/rykVz5pPTcu6cZ9yP0jQNIYQQwtW41XUBhBBCiKqQACaEEMIlSQATQgjhkiSACSGEcEkSwIQQQrgkCWBCCCFckgQwIYQQLkkCmBBCCJckAUwIIYRLkgAmhBDCJUkAE0II4ZJqNYAppboppVYopU4qpXKVUilKqc1KqRtrsxxCCCFcn3stny8YOAV8CZwB/IApwGdKqRhN0+bVcnmEEEK4KHUhZKNXSm0BGmqa1rSuyyKEEMI1XChtYElAUV0XQgghhOuo7SpEAJRSbujBMwSYCFwO3FsXZRFCCOGa6iSAAW8BdxX/XADcp2naO/Z2UEpFAhFlFvsDbYA9xccRQgjhmjyBJsAGTdPSHdmhTtrAlFJNgcji11jgTuBhTdNetrPPXGBOrRRQCCFEXRmnadr3jmx4oXTiWATcjt6RI9HGNtaewNoB38bGxtKqVasaLqUQQoiacuTIEa6++mqAHpqm7XRkn7qqQixrG3A30AKwGsA0TUsAEkouU0oB0KpVKzp27FjDRRRCCFELHG4OulB6IQ4HTMDRui6IEEII11CrT2BKqXeBDPQnrnggHL0X4vXAf21VHwohhBBl1XYV4mZgKnAzelaOLGA38B9N0z6v5bIIIYRwYbUawDRN+wj4qDbPKYQQon66UDpxCHHRM5lMxMfHk5+fj8lkquviCOEUbm5ueHl5ERUVhZubc7tdSAAT4gJgMpk4efIkubm5GAwGDAaDpZetEK5K0zQKCgrIzc0lPz+fpk2bOjWISQAT4gIQHx9Pbm4uoaGhREZGSvAS9YamaSQkJJCSkkJ8fDwNGjRw2rEvlG70QlzU8vPzMRgMErxEvaOUIjIyEoPBQH5+vlOPLQFMiAuAyWSSakNRbymlMBgMTm/blQAmxAVCgpeoz2ri71sCmBBCCJckAUwIIYRLkgAmRD23ZcsWtmzZUtfFYMuWLUycOJEGDRrg6elJdHQ01157LZs3by637dy5c1FKkZSUVAclrVsxMTHccsstFW6nlLK8DAYDISEhdO3albvuuqvav+/nn3+e2NjYah2jNkgAE6IeS0lJYfDgwQwePJjU1NQ6K8cbb7zBwIEDOX36NC+99BJr1qzh5Zdf5syZMwwaNIiFCxfWWdlcmfkLwO+//86SJUu46aab2LJlC/3792fGjBlVPq6rBDAZByZEPRYbG0tRUZHl56lTp9Z6Gf744w9mzpzJ6NGjWbZsGe7u5287kyZNYvz48cyYMYPu3bszcODAWi9fRXJzc/Hx8anrYlgVFRVFv379LO8vv/xyZs6cyZ133snrr79Ou3btmDZtWh2WsGbJE5gQ9cihQ4fo2LGj5fXAAw9Y1s2ePbvUusOHD9dKmV544QWUUixatKhU8AJwd3fnrbfeQinF/Pnzy+176tQpJkyYQGBgIEFBQdx4440kJpaetGLt2rUMGzaMsLAwfHx8aNq0Kddccw05OTmWbQoKCpg3bx7t2rXDy8uLiIgIpk6dWu5YMTExjBkzhqVLl9K9e3e8vb15+umn6d69O4MHDy5XPqPRSKNGjZgwYUKlz1VYWMhDDz1EdHQ0vr6+DBo0iG3btjl+YW0wGAwsXLiQ8PBw/vvf/1qW5+XlMXv2bLp160ZQUBChoaH079+f5cuXl9pfKUV2djaffPKJpYpy2LBhACQmJjJ9+nQ6dOiAv78/kZGRXHLJJWzcuLHa5a4KeQIToh5ZuXIl+/bts7ouNTW1VDXiypUrq1XN5Aij0ci6devo1asXjRs3trpNkyZN6NmzJ2vXrsVoNGIwGCzrxo8fz3XXXcfdd9/N3r17efLJJ9m3bx9bt27Fw8OD48ePc+WVVzJ48GA+/PBDgoODOXPmDKtWraKgoABfX19MJhPjxo1j48aNPPTQQwwYMIATJ04wZ84chg0bxo4dO0o9Yf3111/s37+fJ554gubNm+Pn50fDhg2ZMWMGhw8fpnXr1pZtf/75Z86ePWt5sq3Mue644w4+/fRTHnjgAS699FL27NnDhAkTyMzMrPZ19/HxYeTIkSxZsoTTp0/TuHFj8vPzSUlJ4YEHHqBRo0YUFBSwZs0aJkyYwEcffcRNN90EwObNm7nkkksYPnw4Tz75JACBgYGAXiUNMGfOHKKjo8nKymLZsmUMGzaMX3/91RLoao2maS77AjoC2p49ezQhXNm///6r/fvvv9U+TkFBgTZx4kQNsPuaOHGiVlBQ4ISS2xcXF6cB2qRJk+xud/3112uAFh8fr2maps2ZM0cDtFmzZpXabvHixRqgff7555qmadq3336rAdquXbtsHvvLL7/UAO27774rtXz79u0aoL311luWZc2aNdMMBoN28ODBUtsmJSVpnp6e2mOPPVZq+XXXXadFRUVphYWFlTrX/v377X6+m2++2ebnMQO0e+65x+b6hx9+WAO0rVu3Wl1fVFSkFRYWarfddpvWvXv3Uuv8/PwcKoP5GCNGjNDGjx9vd9uK/sb37Nlj/vvsqDkYA6QKUYh6xMPDg8WLFzNx4kSb20ycOJHFixfj4eFRiyWzT9O/kJYb7DplypRS76+77jrc3d1Zt24dAN26dcPT05M777yTTz75hKNHy0/q/uOPPxIcHMzYsWMpKiqyvLp160Z0dDTr168vtX2XLl1o06ZNqWVhYWGMHTuWTz75xJJNIjU1leXLl3PTTTdZqkYdPZe5/LY+nzOYr2lJ33zzDQMHDsTf3x93d3c8PDz44IMP2L9/v8PHffvtt+nRowfe3t6WY/z666+VOoazSAATop7x8PCw26vvzTffrLXgFR4ejq+vL8eOHbO73fHjx/H19SU0NLTU8ujo6FLv3d3dCQsLIzk5GYCWLVuyZs0aIiMjueeee2jZsiUtW7ZkwYIFln3i4+NJS0vD09MTDw+PUq+4uLhyXfVtJZu99dZbOXPmDL/88gsAX375Jfn5+aW6vDt6LnP5bX0+Zzhx4gQADRs2BGDp0qVcd911NGrUiM8//5zNmzezfft2br31VvLy8hw65v/+9z+mTZtG3759+e6779iyZQvbt29n1KhR5ObmOqXclSFtYELUQytWrLC7zpFxRs5gMBgYPnw4q1atsrTFlHX69Gn+/PNPrrjiilLtXwBxcXE0atTI8r6oqIjk5ORSN3nzMAGj0ciOHTt44403mDlzJlFRUUyaNInw8HDCwsJYtWqV1TIGBASUem8r5dHll19Ow4YN+eijj7j88sv56KOP6Nu3Lx06dLBs4+i5zOW39fmqKzc3lzVr1tCyZUvLNf/8889p3rw5X331VanPWJkEu59//jnDhg1j0aJFpZY7o92uKuQJTIh6aPXq1YB+o/z666/56quvLDdNWzfXmvLoo4+iaRrTp0/HaDSWWmc0Gpk2bRqapvHoo4+W23fx4sWl3n/99dcUFRVZ7SxgMBjo27cvb775JqB3xgAYM2YMycnJGI1GevXqVe7Vtm1bhz6HwWDgP//5D7GxsWzcuJEdO3Zw6623ltrG0XOZy2/r81WH0Wjk3nvvJTk5mYcfftiyXCmFp6dnqeAVFxdXrhcigJeXl9UnKqUUXl5epZb9/fffVgej1wZ5AhOiHlq4cCGtWrXivvvuIzIyEtBvmgsWLGDWrFm1WpaBAwfy2muvMXPmTAYNGsS9995L06ZNOXnyJG+++SZbt27ltddeY8CAAeX2Xbp0Ke7u7lx66aWWXohdu3bluuuuA/T2mLVr13LllVfStGlT8vLy+PDDDwEYOXIkoI81W7x4MaNHj2bGjBn06dMHDw8PTp8+zbp16xg3bhzjx4936LPceuutvPjii0yePBkfHx+uv/76UusdPVf79u258cYbee211/Dw8GDkyJHs2bOHl19+2dLjzxHx8fFs2bIFTdPIzMxkz549fPrpp+zevZtZs2Zxxx13WLY1Dw+YPn061157LadOneLZZ5+lQYMG5YZUdO7cmfXr1/PDDz/QoEEDAgICaNu2LWPGjOHZZ59lzpw5DB06lIMHD/LMM8/QvHnzagfeKnG0t8eF+EJ6IYp6wlm9EC9kmzdv1q699lotKipKc3d31yIjI7UJEyZomzZtKretuRfin3/+qY0dO1bz9/fXAgICtBtuuMHSU9F8zPHjx2vNmjXTvLy8tLCwMG3o0KHa999/X+p4hYWF2ssvv6x17dpV8/b21vz9/bV27dppd911l3b48GHLds2aNdOuvPJKu59jwIABGqBNmTLF6npHz5Wfn6/Nnj1bi4yM1Ly9vbV+/fppmzdv1po1a+ZwL0Tzy83NTQsMDNQ6d+6s3XnnndrmzZut7jN//nwtJiZG8/Ly0tq3b6+99957lmtd0q5du7SBAwdqvr6+GqANHTrUUuYHHnhAa9Sokebt7a316NFDi42N1W6++WatWbNmdstbE70QlWalp4qrUEp1BPbs2bOHjh071nVxhKgyc++5Fi1a1HFJhKgZFf2N7927l06dOgF00jRtryPHlDYwIYQQLkkCmBBCCJckAUwIIYRLkgAmhBDCJUkAE0II4ZIkgAkhhHBJEsCEEEK4JAlgQgghXJIEMCGEEC5JApgQQgiXJAFMCFHjPv74Y5RSKKXKTSAJek7WVq1aoZSq/Wnphcuq1QCmlLpEKfWhUuqAUipbKXVGKbVcKdWzNsshhKgbAQEBfPDBB+WWb9iwgX///bfc3FxC2FPbT2DTgBhgATAamAFEAluUUpfUclmEELXs+uuv57vvviMjI6PU8g8++ID+/fvTtGnTOiqZcEW1HcDu0TTtEk3TFmmatkHTtG+BS4Fk4LFaLosQopbdcMMNAHz55ZeWZenp6Xz33XflJocEKCgoYN68ebRr1w4vLy8iIiKYOnUqiYmJpbb76quvuOyyy2jQoAE+Pj60b9+eRx55hOzs7FLb3XLLLfj7+3PkyBFGjx6Nv78/TZo0Yfbs2ZWamVhcGGo1gGmalmBlWRawD2hSm2URQtS+wMBArr32Wsukk6AHMzc3t3KTQ5pMJsaNG8f8+fOZPHkyK1asYP78+fzyyy8MGzas1IzBhw8fZvTo0XzwwQesWrWKmTNn8vXXXzN27NhyZSgsLOSqq65ixIgRLF++nFtvvZVXX32VF198seY+uKgRdT4js1IqCOgBrK3rsghxoZk5E3btqutSnNetG7z2WvWOceuttzJ8+HD27t1Lx44d+fDDD5k4cWK59q+vv/6aVatW8d133zFhwgTL8q5du9K7d28+/vhjpk2bBsATTzxhWa9pGgMHDqR9+/YMHTqUv//+my5duljWFxQU8PTTTzNx4kQARowYwY4dO/jiiy946qmnqvfhRK2q8wAGvAn4Ac/Z20gpFQlElFncsqYKJcSFYNcu2LChrkvhXEOHDqVly5Z8+OGH3HLLLWzfvp1XXnml3HY//vgjwcHBjB07ttR09d26dSM6Opr169dbAtjRo0d54oknWLt2LQkJCZScqHf//v2lAphSqtyTWZcuXVi7Vr5Du5o6DWBKqWeBKcD/aZr2ZwWbTwfm1HyphLhwdOtW1yUozRnlUUoxdepUXn/9dfLy8mjTpg2DBw8ut118fDxpaWl4enpaPU5SUhIAWVlZDB48GG9vb+bNm0ebNm3w9fXl1KlTTJgwoVRVI4Cvry/e3t6llnl5eZGXl1f9DydqVZ0FMKXUHOAJ4HFN0xY6sMtbwDdllrUElju7bEJcKKpbXXehuuWWW3jqqad4++23ee4565Uv4eHhhIWFsWrVKqvrzVWOa9eu5ezZs6xfv56hQ4da1qelpTm93OLCUicBrDh4zQXmapr2vCP7FHcAKdUJRCnl/MIJIWpco0aNePDBBzlw4AA333yz1W3GjBnDkiVLMBqN9O3b1+axzPcBLy+vUsvfeecd5xVYXJBqPYAppZ5ED17zNE17urbPL4S4MMyfP9/u+kmTJrF48WJGjx7NjBkz6NOnDx4eHpw+fZp169Yxbtw4xo8fz4ABAwgJCeHuu+9mzpw5eHh4sHjxYnbv3l1Ln0TUldrOxDEbeAZYBaxQSvUr+arNsgghLmwGg4Hvv/+exx57jKVLlzJ+/Hiuvvpq5s+fj7e3N507dwYgLCyMFStW4Ovry4033sitt96Kv78/X331VR1/AlHTVMneOjV+MqXWA0Ntrdc0rVJ1gkqpjsCePXv20LFjx2qWToi6c/ToUQBatGhRxyURomZU9De+d+9eOnXqBNBJ07S9jhyzVqsQNU0bVpvnE0IIUX9JNnohhBAuSQKYEEIIlyQBTAghhEuSACaEEMIlSQATQgjhkiSACSGEcEkSwIQQQrgkCWBCCCFckgQwIYQQLkkCmBBCCJd0IczILISw49VfDtV1EQCYdWmbah/j77//5tVXX2X9+vWcO3cOd3d32rRpw6RJk7j99tsJDQ11QknL27lzJzNnzmTXrl1kZGTw6quv0q1bN4YPH866desYNmyY3f1vueUW1q9fz/Hjx2ukfKJqJIAJIWrFe++9x/Tp02nbti0PPvggHTp0oLCwkB07dvD222+zefNmli1bViPnvvXWW8nOzmbJkiWEhIQQExODr68vmzdvpkOHDjVyTlHzJIAJIWrc5s2bmTZtGpdeeimxsbGlJp+89NJLmT17ts2Zl51hz5493HHHHVxxxRWllvfrJ7M4uTJpAxNC1Ljnn38epRTvvvtuuZmTATw9PbnqqqsAMJlMvPTSS7Rr1w4vLy8iIyO56aabOH36dKl9hg0bRqdOndi+fTuDBw/G19eXFi1aMH/+fEwmEwAff/wxSimKiopYtGgRSinLDM7r169HKcX69etLHffjjz+mbdu2eHl50b59ez799FOrn6mgoIB58+ZZyhkREcHUqVNJTEwstV1MTAxjxoxh1apV9OjRAx8fH9q1a8eHH35Y7phnzpzhzjvvpEmTJnh6etKwYUOuvfZa4uPjLdtkZGTwwAMP0Lx5czw9PWnUqBEzZ84kOzu7gt9C/SNPYEKIGmU0Glm7di09e/akSZMmFW4/bdo03n33Xe69917GjBnD8ePHefLJJ1m/fj1//fUX4eHhlm3j4uKYMmUKs2fPZs6cOSxbtoxHH32Uhg0bctNNN3HllVeyefNm+vfvz7XXXsvs2bPtnvvjjz9m6tSpjBs3jldeeYX09HTmzp1Lfn4+bm7nv++bTCbGjRvHxo0beeihhxgwYAAnTpxgzpw5DBs2jB07duDj42PZfvfu3cyePZtHHnmEqKgo3n//fW677TZatWrFkCFDAD149e7dm8LCQh577DG6dOlCcnIyq1evJjU1laioKHJychg6dCinT5+2bLN3716eeuop/vnnH9asWWMJ0BcDCWCiRmmaRpFJw91Nucx/rC1btgBSveQsSUlJ5OTk0Lx58wq3PXDgAO+++y7Tp0/njTfesCzv3r07ffv25dVXX+W5556zLE9OTmblypX06dMHgJEjR7J+/Xq++OILbrrpJiIiIoiIiAAgKirK7u/UZDLx+OOP06NHD5YtW2b5ex00aBCtW7emYcOGlm2//vprVq1axXfffceECRMsy7t27Urv3r35+OOPmTZtWqlr8Mcff9C0aVMAhgwZwq+//soXX3xhCWBPPfUUSUlJ7N69m/bt21v2ve666yw/v/766/z9999s3bqVXr16ATBixAgaNWrEtddey6pVq8pVk9ZnEsBEtWTlF5GQkUdydgFZ+UXkFhjJKTCSW1BEToGRvEITpuJZvw1uyvJyL/7Xz8udEF9PQnw9CC7xr8GtboJdSkoKgwcPBiAhIYGQkJA6KcfFat26dYDe66+kPn360L59e3799ddSASw6OtoSvMy6dOnCrl27Kn3ugwcPcvbsWe6///5SX7aaNWvGgAEDSvVA/PHHHwkODmbs2LEUFRVZlnfr1o3o6GjWr19fKoB169bNErwAvL29adOmDSdOnLAs++mnnxg+fHip4FXWjz/+SKdOnejWrVup815++eWW6lAJYEJYkZlXSHxGHgkZ+SRk5pOQmUd2vtHh/Y0mDaNJK7UsLaeQM6m5pZa5KUWQjzuNQnxpFuZLkxBfdv+1Haj5p6LY2FjLjSE2NpapU6fW6PkuBuHh4fj6+nLs2LEKt01OTgagQYMG5dY1bNiw1A0fICwsrNx2Xl5e5Obmllvu6Lmjo6PLrYuOji4VwOLj40lLS8PT09PqsZKSkipdzsTERBo3bmy3jPHx8Rw5cgQPDw+HzlvfSQATNplMGucy8jiWmM2x5GySMvNr57yaRmpOIak56ew5k05OZhpPXTcIBew8dIIOMQ2d9oR26NAhxo8fb3l/7tw5y8+zZ8/m5ZdftryPjY2ldevWTjnvxcRgMDBixAh++uknTp8+bfcmbb7Rnzt3rtx2Z8+eLdX+5Wzmc8fFxZVbV3ZZeHg4YWFhNntOBgQEVPr8ERER5TqqlBUeHo6Pj4/VDiDm9RcTCWCilLxCI8eTszmWmM3x5BzyCh1/wqop//yxBpNRL8fzb33K4DETaR8dSMdGgUQGeFfr2CtXrmTfvn1W16WmppKamlpq2xkzZlTrfBerRx99lJUrV3LHHXewfPnyck8uhYWFrFq1iksuuQSAzz//nN69e1vWb9++nf379/P444/XWBnbtm1LgwYN+PLLL0tVI544cYJNmzaVagMbM2YMS5YswWg00rdvX6ec/4orruCzzz7j4MGDtG3b1uo2Y8aM4fnnnycsLMyhNsX6TgKYQNM0TqfmsvdsOkcSsig0ahXvVIMSTh/jo6fvtbzPSDnfLXn5uy+y7tvz3z4f+u/7jOzflbbRAXi5Gyp9rnvuuYdNmzbxzTff2N1u4sSJTJ8+vdLHF7r+/fuzaNEipk+fTs+ePZk2bRodO3aksLCQnTt38u6779KpUyeWLVvGnXfeyRtvvIGbmxtXXHGFpRdikyZNmDVrVo2V0c3NjWeffZbbb7+d8ePHc8cdd5CWlsbcuXPLVStOmjSJxYsXM3r0aGbMmEGfPn3w8PDg9OnTrFu3jnHjxpV6snfEM888w08//cSQIUN47LHH6Ny5M2lpaaxatYr777+fdu3aMXPmTL777juGDBnCrFmz6NKlCyaTiZMnT/Lzzz8ze/ZspwVUVyAB7CKWkVfIvrMZ7DubQXpuYY2d5/j+XQDEtO/m0Pb7t20g/sQRq+tyM9PJzUy3vN+47mdUcAN+O5RI66gAujcJJjLQ8acyDw8PFi9eDGAziE2cOJHFixfbbHeoac5I4XQhuOOOO+jTpw+vvvoqL774InFxcXh4eNCmTRsmT57MvffqX1oWLVpEy5Yt+eCDD3jzzTcJCgpi1KhRvPDCC1bbkpzptttuA+DFF19kwoQJxMTE8Nhjj7Fhw4ZS48UMBgPff/89CxYs4LPPPuOFF17A3d2dxo0bM3ToUDp37lzpczdq1Iht27YxZ84c5s+fT3JyMhEREQwaNMiSYsvPz4+NGzcyf/583n33XY4dO4aPjw9NmzZl5MiRxMTEOOMyuAylaXX7bbs6lFIdgT179uyhY8eOdV0cl6BpGseSstl1Ko2TKTnU9K8/OyONOdcNAAXPfL0Z34CgCvcxFhXy+fwH2P2b/cwMXYeM4sZHXsbgXjqwtIjwo3/LsEpVLyYkJBAVFWVznbkrdk05evQoAC1atKjR8whRVyr6G9+7dy+dOnUC6KRp2l5HjimZOC4SRpPGnjPpfLblBMt3neVEcs0HL4A9m9ZgMhkxGY3s2bTGoX0M7h7c+MjLdB0yyuY2toIXwNHEbL7YepIfdp8l0cGOJytWrKjSOiFE3ZEqxHour9DIP2fS2XUyjaz8oop3qKbKtF/dOvdNIhrFWD2Owd2DCfc8ZfMp7Jp751gNXmaaBkcSsvg3MYtWkf70axFGuH/5FEZmq1evBvSeaIsWLULTNKZPn05ycjKrVq0qNy5JCFH3JIDVUzkFRew4nso/Z9IpKDLV2nkr0361f9sGIsbH2DnWepvr9m1bT5/LJthcb6ZpcDg+iyMJWXRqGMSg1uF4e5Tv7LFw4UJatWrFfffdR2RkJKDn2luwYEGNdhwQQlSdBLB6Jq/QyF8nUtl5Kq1WA5fZoKumcHzfTofarwaOnWx3mwM7fgfANyCYa2fMBU3j29efJiczjQPbNzoUwMw0Df45k86/iVkMah1Ox4al2+LCw8OZN29eqWWRkZGlsj4IIS4sEsDqiYIiE7tOpfHnidQ6Hbtlbr8CbAYxe+1XJU2490nCGzZl8NU3ERCi9z5r2aUvG2M/ZciEm6tUvpwCIz/vjWfv2QwuaRdpt1qxtrlyhyohKqJpmtPzoUoAc3FFRhN/n0lnx/GUSqV1qknVbb8y8w8KZfTU0tV3ASFh5ZZVxZnUXL7YepLuTYPp2zwMT/e67c/k5uZGQUFBjfwnF6KuaZqG0Wi0mXqrqqQXogs7GJfJJ5tPsOFgYrWC1/H9uyxjtay9r4qK2q+cea6qMpo0dhxP5bMtJziVklMnZTDz8vLCaDSSkJAgT2KiXtE0jYSEBIxGo9W54KpDnsBcUFx6HhsOJXA2La/ax8rOSGPhrMmWcVqappV678i4LWvstV/9+ev39LlsQrlzV/Vc1ZWRW8h3f52mV7NQ+rcMq5NM+FFRUeTn55OSkkJ6ejoGg0GexITLMz95GY1GfHx8bI61rCoJYC4kM6+QP44kcSAu02ljuMzjtMw/a5pW6n2fy6+p0nGttV81aNGel26/gn93byMnM73cuat6LmfQNNh+PIWTKTmM6hRNqJ9zqzoq4ubmRtOmTYmPjyc/P98yo7AQrkwphaenJ15eXkRFRZWaFNQpx3fl6oqLJRNHQZGJHSdS+OtEarXzFJYdp5WWFEd+dhYAqviPSyu+eXr7+RMUfj4HXNlxW5VNEbXmy7dZ+dGrAASGRVJYkG/pVu8TEERg6PlsF/bGiNU0D4NiSJsIujQOrpPzC3Exqkomjlp/AlNKBQBPAt2A7kA48LSmaXNruyyu4FB8Jr8dSiQzzzmDkO2N09LKfOvPy84iL/v8tiXHbTlS/Vc2WKbEnZ8qIiM5odS2lR0jVpMKjRq/7k/gWFI2l3aIwtdTKiqEuBDVRSeOMOBOwAuIrYPzu4TkrHy++/M0K/4+57TgBfo4LXspmmzp2H9EqXFbjqSIMgdL86sw37E2u7LnqivmlFTn0is/OaIQoubVxVfLE0CIpmmaUiocuL0OynDB2bJlCwA9evVhy9Fkdp1KKzd7sTM4Mk7LmqZtOvHytKst78umiFq9+C0APL28LdV/jg5qLqtj32EOdbO3prLVmhXJzCvimx2nGdZWqhSFuNDUegDTXLnRrYakpKQwePBgNE3j5eXb0Tz9nHJcWzfzisZpWbNmyTs2n6BsVf+5GTy45v/+R2F+CPu2rgLygTzAfgeFHz54mQ3LPrG8d7Q9rKZ6NRpNepViXHoel7SLxN0go0+EuBC4TOW+UioSKDunRcu6KIuzfb7kG4qK9GrCretWOaU3XkU3842xn1bqeNaDlwFoD/QAOgCRePu2YscvvVj/rTtZaQaKCt2AL8rsV4QeyPKBeOBQiddBcjMPkZtpve3Nnpru1bj3bAZJWQWM6dqAQO+6mRtMCHGeywQwYDowp64L4QyHDh1i/PjxaEBugZGE+DjLuspkbLfH3s08OyONNV++A4CPfxATZz7Nb0s/4fi+nXaOqIAuQM/iVw+gK+BTaqu8HDhtvY9ICe6Af/ErDD34lZUGbKVBi5M07zQcTTNRdliUI5nvC4oD790vfOCUXo3xGXl8sfUkozs1oGmYb7WPJ4SoOlcKYG8BZafMbQksr4OyVMvKlSvZt2+f1XVV7Y2XcPoY7z6mNyd6eHnbncak5yVXAXpN7uX/uYduQ66gVde+rPnibUwmE6kJZ9i7eS16gLkUGANcCdgbhFiA/jSVUPzSf+7YvzPpSf9y+vBhPDyD6dBvFAV5Rg7s2I5m8sTdsyVFBc3Qf5Ulx14FA5dz7ii8eg8EhBbRpnsObXtm06Z7DoFhxkplvt+9cRWtuvYDqt8+lltgZNnOMwxuE06PpiHVOpYQoupcJoBpmma+M1q4YqaCtJwCGg+4mq5DVjslY7vZ7t9WleqmXlLZm/kvX7xl+Xn152+yeeXXlvcT7vmQ4/sas3fzYWAYpYOKWTawE/ir+PUnsB8onc5KT9p7JbnZIfy2dC+Drx5lGdScmdqPjbGf0mNEc166vQN6dWQM0Kb41ZngiBtJS9RTz2SmuPPnr4H8+WsgAK275dD7sml0HrSHf37/ocLr4+Hp7dT2MZOmseFgImk5BQxrE4lbHWTvEOJiV6cDmYt7ISZSxXFgrjSQuchoYtvxFP48nkqRScNYVMjn8x+odsZ2c0eNc0cP8s2Cp6pYOl/gGuA2YKiV9TnAGmAFDWLiOXf8B3z8A4hs0pwTNvIYOlr+bau/Y8krj1ldd/3sF2jR6XoO/enLwT/9OLLbh7yc0nN5+QYU4Rv4I0lnngT22DyPh7cPhXl6d/jAsEjuefkzpw2Ubh7uxxWdo/FyLz/PmBDCMS4xkPli9G9iFhsOJpKeW2hZVt2M7Qmnj/HBU9NIPH0MAB//wCqUrA9wK3ADUHb/k7h7/MKQa8KJbnqS2EVPkZOZRljDkXTsfyeDr74JgDnXD6hS+c3s5Uw8uGMjfS+fQESjdAZelY7RCCcPeLP7twD+/DWQ7AwDOZnu5GReDVwNbAHeAz5Hr9I8zxy8QB9EvfD+KfgGBgPVz/pxLCmbr3ec5upuDQmQzh1C1Jo6CWBKqSsAPyCgeFEHpdS1xT+v1DStblODO0l6TiHrDyVwNDHb6vrqzDi8f9sGS/ACyM3KcLBUfuhB606gU6k1/sFF9BqZQYe+pzn45yKGjDfnMWxP2569LPNw+QeFAvrTU1XLb1aZOb8MBmjeMY/mHfMYc1sSezb78ccPfvz7t7k6sF/xaw7wHPARUIg1malJZKYmAc7J+pGUmc+SbacY160hkYHe1TqWEMIxdfUEtghoVuL9xOIXQHPgeG0XyJkKjSa2l6gutKWyMw6X7HWnaRoGdw+MRdZv0OVFAv+H3pkztMRyI7CSGx5oQ49LFAZ3AG9ada14Hi5nzJhc1Tm/3D01ug3NIj/3Y/79+z1gKnoVaCOgKfAO8Ch6IPsYvet+eeZ2RmcMgM7KL+KbP08zqlM0LSP8q3wcIYRj6iSAaZoWUxfnrQ2H4jPZeDiJjNyKA0tlZxy21+vOtlbAbOAWoOSTwRHgfeBT4BwaL2BwrzjgOFr+5p16cnz/LqdlxDCz3nU+HZiDl+9reHnfRVb6PZiMjdE7hbyHHsjmAZ9RMpCZ2+nycrKd1sGjoMjED7vPMqxtJHlnDgDQr1+/Kh9PCGGbpBRwkviMPL7efooVf59zKHjB+acP880fzj99mKvp4Pykj5XLY9gDfdTBQeBuzgevLcB4oC3N2q/CN0CvrT2wfaODx624/EOvmcoHT97Nwvsnk1Oi92NlWZvssmx+xZK9K/NzUslImY/J2AK9ivRk8ZoWwIfAXvRhATpzO50jeR0rQ9Pgx22HGDRoEIMHDyY1NbXaxxRClCedOKopK7+I3w8ncSAuw2lzdJVUNqPG+OlPsvu31ZjHcZXXHngWvVdhST8ALwF6tZ/56SMnM8PmE19VOSMjhq1MIo7lVywE3sPHP5YGMfM5uucyoDF69/yfga+AWbw8bRw+/oF2x8xVZyC50ahfg6VLl3HbbbdW+hhCCPskgFVRodHEnydS+fNEKgVFNTf5oLUJJ60Hr+bAXOA/6FkzAFUE2mfAy0DpgdPmpw9H2psq4khGDDNHA4KtIOhIMuKO/S8humkrhky4hWVvzQOm4+45G1PR45hMvsD1wBVkJD9JRvKblBy/Vp2B5Lauwcz7Z/PK/14x/1aIjY2ldevWFR5TCGGfBLBKMpo09pxJZ9uxFLLynTfNiZm1CSfNvnr1CSt7NASeQE/qr3fhVkojOmYLnQZu4pfPH7B6Hkd7CTqiMhkxbAWEygbB3pdNsBnArp/1HP7BehXs+Xa6MRiNcSx/O4LdvwWgDxtYgN42OA3YWu44lRlIbu8aZGWksX9fmuX9ypUrmTFjhkPHFULYJgHMQSaTxr5zGWw9luJwG1dVOD7hZAB654SZlMxHGBT+O3fPb0xU0zA+fW4tUL1egiXZ6qnn6LQp9gJCpdJC/baKnz5+zeZ5Sgbn0r0ci7j5iXPs357Ot6+HkRrvgz6n6ib0TGUPAbmWsjoyENvM0WswfsI1TJ8+3aFjCiHskwBWAU3TOBSfxZajyaRkF1S8QzVVfCN0Qx/HNY+SuQn9g7fSvvcvjL1zJP5Bejkr28vRHnvZ7R2p1qsoIFQmCPoFhWDOIOPp7cOkB16oVHBu3zuHEdcv5tvXc9C/BHgD9wIjgCnATocHYps5eg0uu/c5sgo0QmS8sxDVVqeppKqrJlNJmUwaRxKz2HYshcTMfKceuyK200xdAvwPPQu82RZadv2Wu1+YWuVJIB2xddW3fPW/xwGYNPt5qx0zMlOTbWbmeObrzZZqPVsqSq/l5etPcHgUGalJlicyb78AgsL1QG4yFtG6Wz9G3TyjVC9OsF5FqR+jFfABMKR4TQHwJNfdH0q/UePtltcaR66Br6eB8T0aERkgA56FMJNUUk6QV2hk79l0dp1Kr9GqQnvKp5lqjd4R46oSW50EHgaWcPPjm50evKrSMaM6mUWg4vRa+TlZxJ/MKrUsLzuTvOxMy/uBYyeXC15gr4ryCDAc/Vo+jZ68+EVWfPA3bXu4ExJZuXZOR65BToGRb/88zbhujWgU7GNzeyGEfTIOrFh6TiHrDybwwe/H+O2QYwORa5J+I/QHXkQfv2QOXlnA40BbYAmg3xidf37b461yM9NLrdu/bQNQOjPHTU+8xk2Pv4pvQLC+zsFxZvYCQEXstbHZH0NnAl4ABuDpfQqA7PQuvHx3M3aur1xGDUevQX6hiWV/neZ4kvU0Y0KIiskTGLDuQAJ/n07HZKc61RmphhylabD+2zzgAHpqJNBvsh8BTxDToQFDxr9QrQ4ZFalKxwxntLnZS08VHBFNWmKc1f0qamNzrI0qnIkzsvjx/TS2/BRMbpaBz55vyKG/0plwbwIenhVXtztyDUr+LX2/+yyjOkXTJirA1iGFEDZIGxiweOsJEjJst3NlZ6Qx57oBNlMNOTO4xR335JsFoRzbez47fEyHXEbdfJh9W/4HwMjJd+MfFEpmanK5BLvO5KwpXyojKz2F35Z+UioAmD9njxFjeen2K63u50gbm/lYjrTT/fOHH1+/Gk12hj5FSpM2edzy1FnSk3cAVf9dW/tbUgpGto+iU6PqzVEmhCurShuYVCE6wF6qIXPvvOqmTcrLdmP52xG8fHezEsErgT6X/8i9/ztFm+7eXD3tMa6e9pglWFlLO+VM5jYpWyrbU88R9tJrndy/2+Z+jlajVtRGZdZ5YDYPvHOcmA56t/pTh7x5ZXoT3pj5brV+19b+ljQNftkXz58nUqp0TCEuVlKFaEVlOjD0vOSqaqVN0jTYsDSXH99vgMlonp7eCLwJPMU/m+DEgQjL9tWdu6qyqtsxw5mckf2+MscICjMy/b+nWP5OJH98H0xOhgewGs34EP/8sYa+oyr+XVfmbylx7pvkD+7BgJbhFR5XCCEBzKrKDKr95Yu3LD9XNm1SwmkPlr4RyaGdfiWWbkQfk/R38fmoUmojZ3FG0HAWZ7SxVeYYJYNPYNh4MpKfRh80/grfvr6Utd9MxM1NT4Zs63dd6SwljWIoNGoMbRNhdR8hxHkSwKxwtAMDQGF+nuVnR9MmFeQrfl0SytqvQzAW6rW47h5pFBXOQJ/exLrKpDZyFkdv+N4eBvy9DPh4uuPraSh+6T97GNwwmjSMJo0ikwmTplFk1N9n5ReRllNIak4BOQVGW8UAqj53WFWPUTr4/BdYAywFYjAWTSDxVBv03qHHbP6uq9IZ5q/i/Joj20eilLK7nxAXMwlgVjjSY60itoLN/u2+LF0YSfI5T0DPW9hvdDqjbopj6ZsJ7P7N9vGc3WHCEdZu+A2io7j3oSeIDPAmMtCLqEBv/L2q/6eUV2i0BLPU7AJOp+USl56H0c6koDWpfPDZCfQCvkSflqUTsJVGrR+3+cWiqllK9pxJp9Bo4vKO0RjcJIgJYY0EMBsqGlTbsd8l7N2y1uo6azektCR3li+KYPfG892lG7XM49r7EmjWPg9ws3u+mugw4Sh/L3diwv1oHu5HgyBv/JwQrKzx9jAQHWQgOuh8hor8IiOnU3M5mZzDyZScWknnZWY9+CQDo9DHjT0ERHDm8Ots/+UM/a6w/gRZ0d+Srd/twbhMCo0mruzcAHeD9LcSoiwJYGWU7BJvrwND6+79bQawkjckoxF+jw1m1afh5OfqNyEvXyNX3JzMwKvSMBjO73ehdJhQCiICvGgR7k+LCD8iA7zqrCrLy91Aywh/WkboA4oz8go5npTN3rMZxKXnVbB39VkPPib0zB2HgEWAN1+/2pLM1CRG3pCCtUu1MdZ21bC93+3RxGxid53lqq4N8XSXICZESRLASiibsNZeB4Y/135v8zjmG9Lx/d58uyCSs0fPP1F0G5rBuLsTCQor/229rjtMRAR40bFhIK2jApxSJVgTAr096NI4mC6Ng0nMzGfP2XQOnMskr9B++1l12P5i8QFwHPgWCOanj8NZ981vBIY9wW1Pv2bp1JGdkcaaL98BwMc/iIkzn67U7/ZUSg5L/zrN1d0b4e1hsLmdEBebC/MuVUfKTqJorwND3IkjnOIfq8Fmz6adHN93N1t/CkLT9K/j4Q0LuObeBNr2yrGcr+wAaGdmj3eUl4cbbaMC6NQoiKhA10ouGxHgxfC2kQxuFc6RxCz2nMngdGqO02fGNn+x8PDyLtVpR/crMABYATQnL3sCedmh7P7td0beEANQPN5LL9Tl/7mHbkOuACr3uz2Xnsc3f55mQvdGNVaFK4SruagzcRw6dIjx48eTnltIkUkrkaEcfAKCCAy1Pf7KWsaI9ORkvnntECf2TyY7Q69CNHiYGHF9KiMmpZRKRVRRdo+a1jjEh44Ng2gd5Y9HPWpfSczMZ8vRZP5NzHJaIDP/rgeMmczyd5630ZYVASwH+gNgcD9KSORtGDxOV+rvqiLBvh5M6NGYIB+Zj0VcWPIKjdWqIZBs9JW0cuVK9u3bZ3VdRV3iy/bOO/OvJ0sXduPY3v6WZW16ZNPvig0ER2bh4dmt1PHLPu1VdgB0VbgpRdtof3o2CyUiwKvGz1cXIgK8GNu1IQmZeWw5msK/CVkV71SBkr9r2z0KE4FLCIpYRXriUIxFLUg6+wV6h4/SiZDt/V1VlJYsLaeQb3acYnz3RoT518/foXAtKdkFrD2QQLvogFpPh3ZRP4EVFhYyZcoUvvnmG7vb2evCnpvlxqpPw/j9+2A0k15dGBhaxLi7Emnd4xRzr9efsu57bQlfvvyoZT9nfiuviKe7Gx0aBtKjachF9809ISOPzUeTOZrovKzv9vIpzl2ymXXftGLDd+ZsGunoY8XKj48o+3dVmadyH08D47s3crlqX1F/FBlNbDuWwo4TqRhNGpd2qF4+T3kCqyQPDw8WL17MiZQctv26wuo2toLXsX27OLCjGVtW9CQzVb+Mbm4aQyakctmNyXj7amxddf4pa8N3H1cuI4MTsm34ehro2iSYbk2CL9rG/8hAb8Z1a0Rceh6/HUrkTFputY9pr7fogR3rGXdXKAEhRfz4fjQQBKwGJgPLLNtZHftViafyXMucYg1pHOJbvQ8kRCUdT8pm3cEE0nL0aaeO79/F7txQOjUaWavluKifwMze+ulP7hndy+o6a1nOj+wu4K0HU4HBlmUtu+Qw7Nqd/PjBLZZlJZ+yvP0D0Uwm8nPsV2k5Y8Cyt4eB3jEhdG0SXK/at6pL0zT2ns3g9yNJ5FaQ9cOeT5+bxa4NK6124Ok2dDQ3Pf4qAH98r/huYQzggd71fjqg90Z85uvN5GSl25gl2vGncg+DYnTnBrSIqNy8ZUJURVZ+ERsOJnIo/vwksuaaA6UgMTGRkJAQO0ewTZ7AqmjXJuvjuaD0GJ2sNAMrPw5jy8ogQK8u9PbL4pr/y6TH8Ex+W7bC5lNWXlZGheWobvDydHeje9NgejYLwcv94nziskcpRadGQbSK9Of3w0nsOZvuUEePqvYW9fD6FlgPfIc+OenbQDTwNPu2rSc3K6PaT+WFRo0fdp9jZIdIOjaU6VhEzTCZNHaeSmPL0WQKikyl1pWsOYiNjWXq1Km1Vi55AgP6XzqWLWt+tPqNunX3/lw9bS6LHt5Bdtq9aJp5qpNCYCHe/v8jKExvh9A0jdCohg7PPlyWo3NaleXupujSJJg+MaH4eErgctTZtFzWHkggMbPqc8HZY35S8/IdDvxIfo65qu9tug5ZwY2PzLc735qZo19sBrcOp1dMzUytIy5eZ4r/nyQV/z+xNsOC+ctWSEgIDRo0sKyLjY2ldevWDp1HnsCq6OYHnsE/onGpb9QNWrTnpduv4MiuMN56oBvZGZeV2OMnYBZwkLwsyCtRK9hj+JUc2P475nE/ZTVp04lTh/ZYXVfZbBtuStGxYSB9W4QS4H1xdc5whobBPkzu05Sdp9LY/G8Shcbyv7Pq9BYt+aSWmx3PO482JjXeA7gbY9G1aFpKlfIk2rLxcBLZBUaGtA6XJMCi2nIKith4OIn95zJK1VTYm2EhNTWV1NRUy/uVK1cyY8aMGiujPIFhfUbm1Z9tYPVnjYErSiw9iB64frJ6nK5DRtGmxwC+ec32JJCdBoxkz6Y1FbafVKRpqC9D2kTU2+7wtS05K5+f9sSxd/+BardL2ZKebOCdRxsTd1z/nbXrnc0tT54lPzfJoVmiHdW+QSCXdYjCTZIAiyrQNI2/T6fzx79J5Beayq2vaKZ2s4kTJ7J48WI8PBz7ci1PYNVgfiw2FkWTlTaTvOzbOD9hdTr+wW9SkP8yBbmpVvf38vUn7vhhDu3cbPc8odGN6X3ZeMbc9mCVsm2E+nkyuHW4NNo7WZi/Fzf0acrM7z+rsd6iQWFG7nn5FO890YiTB3w4sN2Pdx9rRLdhP9jcp+RTeUVjxCzlOpdBXqGRK7s0kE48olLOpevVhWW/0JfkyAwLlQ1eVSVPYOhPYIvf+YSfPvIG7gPMY2uK0PPdPYk+UNUJlMLNza3S7SneHgb6tQilS+NgmV6jBhUWFjJ+4iRWLF9qdztrVXuOBpj8XMWHcxtyuHgiUx//I+Rm9cc3oMjmU3lV2uIaBHlzVbeG+HrK91RhX05BEb8fTmJfmepCe+yNh0xISCAionKTslblCeyi/3qWlwcrFgey/pvH0KfHMAevZejzPd2NM4JXh37Dadd7CGgaJqOxOD9exdyUoluTYKYOjKF70xAJXjXMw8ODZd8sYcI119rcxlrwMieCXnj/ZHJKPKVZ4+WjcfuzZ+k0QG88zc1qhW/g39w9fw3dhlxBt6Gjefj9lYy84W4m3PskcL4trjJ/O+fS8/hq+ylSa3EKGuFaTCaNnSdT+XjTcfaedTx4gf3xkCtWWB9X62y1/gSmlPIH5gHXAaHAAWC+pmlLqnCsaj2Bbd0K110HJ0+eX+YbuJecjDuBTeW2D46IJi0xzmr7lXmdNV6+/gSHR5GRmlSp9pRGwT4MaxdBZIBkW6htCQkJREVFWV1nrV1q66pv+ep/jwMwafbzDnX2MBrhq1ei2LFGf5oKiSzk7hdPE9Go0G5Pr8q2xXl7GBjbtYEMeBalnEnLZV0FvXDtKTsesnOjIOY/8QDJyclcf/31LFlSuVu6q7SBLQV6A4+gT6g0GfhSKeWmadoXtVmQ5s0hJUX/ObJxAVfelkiz9hnMnVQ+eAHc9cJH/Pnrcqvjf3qMGMtLt19pdb/8nCziT5YewGyvPcXX08Cg1uF0aBAovcnqiL1vkPu2rSemQ/dyAcZs+bsvsu7bDy3vbQUYgwEmPRCPj7+JjbEhpCZ4sPD+Jtz1wmkO77Ld06uybXF5hUaW/XWGSztG0S460OZ24uKQlV/E74cT2X8us+KNyyhZTV52POSlHaK4acJoFixYwKxZs+wfyElqNYAppUajz8U+WdO0L4sXr1NKNQP+q5T6StO0mpvYqYzISLjrrmMcTEhj6GQ/DAbYtnq9ze1PHNhVKoEvQEBIGKOnzmLb6u+qXI6uQ0YxcOxk3JSiS5Mg+rcIu2hTP10oVq9eDUBYWBiLFi0iv9DItOnTyUpP5cD2jU4ZhAzg5gZXT0vEx9/Ez5+HkZnqzpsPNuH2Z26n65CdDo0RGzh2coWfp8ik8dM/caTlFNKvRViF24v6x1hcXbj1WEq5wciOKDtfYtmE5gCRkZE899xzzipyhWr7CWw8kAWUzZ77EfAF0BdrdXc1JCUlhQUL2mDSNPpdrTeMV3VSSXv72ateNLenNAkP4JJ2UdIt/gKxcOFCWrVqxX333UdkZCQAIy8ZzkNPz6f1iOvx8Qvg+D7nBBilYNRNyXj7mvj+3QhyMw28+1gMU+e8RXbGZI7s2mLz2JXN3LL532QycgsZ0T5K2lMvIseTsll/MIHU4tyFVVEXM2hUpFbbwJRSmwGDpml9yizvCOwB7tI07V0b+0aiT7pUUktgeVXbwD788ENuu+024Hy7hbV5vjJTky3d3P2DrI/JsbefverFF5dt5YpebenUSKoLXcXh+Ex+3hdPbl6+3fEwjgaYktUym1cE8e3rkWiawt3TRFHBVeiTZZZX1cwtAE1CfbmycwPJ3FLPpecUsv5QQpVmY6hsO+zkS/vW+2z0YcBRK8tTSqy3ZTowpzonN09gaXbu3DnLz2XbLXpfNt4SiMzVhPZYe5x2pHoxKOkfOjfuY3O9uPC0jgog1M+TH3afZcI9T9kMYNfcO6fC4FW2Wqb/leDla+KLF6MpKnBD7w07Gfi23L6VzdxS0qmUHL7YdpKxXRtIJ6F6qKBIn+pk58lUikxVe0ixl3HDWjU5l/at0nmqoy660du7mvbWvYXer73ka1xlTmyewNL8KpnyJDcznfgTRyyv/ds2VObQdpWsXrzpideY/swbBIfo35zXrfnFaecRtSfM34tJfZoSt8d2jfc+O92MzUp2j9/04xJevONKfvliCAFhdwN56JnslwC3EhQeTVB4FMpN/29b1ZybZhm5hXy9/RQH4ipONC1cgz7jQjqfbDrO9uMpVQ5eAIOumkLXIaMq3M7RdtiaUNtPYMlYf8oy14OkWFkHgKZpCUBCyWWVrXK755572LRpk0MTWDrzF2LurTN8ws2M6NGaXjGhzLnrulrtrSOcz9vDQPLB7YDjbabWqmXMfvniLQrz84rfHSl+fQ/4AR+QnjQDeB2AJm070/uy87UJVVVo1Dt3JGTkM6hVuKSfcmFn03LZcCiRuPS8ijd2gCMZN5wx/VN11HYb2LvADUCIpmlFJZZPAr4EBmqa5nAnjqqMA6toFuaa+oU0DfVlRPtIgn09nXpcUbeSkpJ47bXXGHjVjRzMcEPT7LeZblj6McvffqESZ+gPrASCi9/PoUO/3zmwbUOVMuTb0zTUlyu7NJAesC4mM6+Q3w8ncSCu8t3iHTq+nYwbJdth62JG5tquQlyGPjFS2e4rNwNnga01XQAPDw8WLlxoc70j7RaV4eNp4PKO0VzTs7EEr3ooPDycefPmcUWfdlzRqQHubsrS9mmtw4+j1TLnbQaGAfHF75+mMP85TCZTpbJyOOJkSg5fbD1JQqZzvsGLmlVQZGLTv0l8sum4U4PX8f27LB2LwH7GDUeqyWtSrQYwTdN+An4BFiml7lBKDS9+KhsFPFRbY8AqGqTqLB0aBnJz/xg6NJTBoxeDttEBXNOzsd2efeZqGXtBrEO/4WWW7Eaf/VtPGXN4Zz/gfcDA8ndf5MU7rrS8Es8cr9ZnSM8t5Kttp9h50nrSalH3NE1jz5l0Pt50jK1HU6xOA1RV1lKilW3Dv+nxV/ENCNbXlWiHTYhTmCo/vKxa6qITxwTgM+AZYBX62K8bNE1bXFsFKDlI9euvv+beeQut/kKqKsTXg2t7NubyjtHSTfki0zDYhxt6NyXE1/ZTvMHdgwn32J5yp013a9U1h4FB6FP6ANwKLCE3M9fpnY+KTBrrDyayfNcZcgsq951yy5YtbNlifdyaqL4Tydl8vvUkv+yLJzvf+d/3reXcnHDvk4y84W4efn+l1VydmakGli2K4IqBASy1nwPb6S7KbPTmdgvzINXFW0/w74mzFY71qojBTdGrWQh9mofiLtNYXNRyC4zE7jpjs0F92+rvWPLKY1bXNW3bhZMH/7Z0DDEVFfHly49iLCoEIoHVQLfirVejfyfMARRPf/U7ASHhTvsc/l7uXN4xmqZhFedRTElJseSPTEhIICQkxGnluNglZeXz++EkjiVVfjyXPdXJuZmd4ca6b0L5PTaYgnz9fte+PezZo2eYqSxXGAd2QTC3W5TkyFgvexoGezOifRTh/pJJQ+htn9f0aMyPf5/lRHJOufX2MrcEhkUy8oa7Sw2Kb9G5D28/cgsJp46it4mtAAYClwM/A2OANPZv2+DUDAlZ+UUs3XmaXs1CGdAyzG4vxdjYWIqKiiw/T5061WnluFil5xay+d9kDsRVLlO8oyo71itifAx52W5sWBbMhm9DyMs5X8PUd2AR/3vJvUrBq6ouyiewsqzNyOwoT3c3BrUKp0vjIMmkIcoxmjR+2RdXLnFqVTK+lO4N5oueF/vy4vf/AJfjE5BTrVmj7YkO8mZUx2hC/PTOSNYSA5jHVoaEhNCgQQPLutjYWFq3bu2UclwMcgqK2HYshX9Op1drLFdFjEWFvPPY7TbTlZl1HTKK62e9wuaVEaz9KpTsjPOBq1n7XEbfksT0KYH1PhNHvdIq0p/h7SLx95LLKKwzuKnitlB3/jpxvmOEvcwttpTuDZYDXIXenHwd0BnYRG7m5eRmHiqxT+lEwiXTVjk6AadZXHoen285QZ/mofSKCbUkBrAmNTW1VKKAlStXMmPGDIfOczHLLzLy14k0/jqZWqWEu5WVl5PNv7u36Qk5bT7MGGje8XVeuiOGtKTzbbuNWuZxxS3JtO+Tjf7dvfY7q8mdtwr8vdwZ3i6SVpH+dV0U4QKUUgxtE4Gfp4HfjyRVuSrIXO3o4x9EZJPmnNi/C31YZSJwDxAD/A6MBnaUG5BfMm3VIx+sKpXCytGxZEUmjU3/JnMoIYsJN97qUGKAiRMnMn369Ep/3otJQZGJf86kseN4KjkVdJypzpeQsvZsWoOm2QuUo4AXiV10/uk5qmk+o25KpvOgrFqtLrRGAlglKAVdGgcxsFU4Xu7Su1BUTq+YUHw93fllXzymKkSxkvMv+QYElkgkfC96kpqn0fNdr6NFpxe48ZHJpcY0lswmvuaLRdXKLJ6Umc/SXXHc8dSrmDSN774tn6sR9OC1ePFiPDzqJlPDha6gyMTu02n8daLiwAXV/xJi7rRhzviSm21r/FgP4CVghGWJmyEO/+AFaCylYcs3cHOLqfgD1jAJYA4K8/dkRPsoGgX71HVRhAvr0DAQLw83Vv59rtJtG2WrHUsnEn4GfbDzW4A/J/Y/y8bYf9i6+moAigoLSIk7bdl328/LLD9/9eoTrPnybQweetuWo+1mmgZ7zmUz4D8P2gxgb775pgQvK/KLjOw+lc5fJ1MrNVShul9C7HXa0DUFXkBPIG2WAbyAybiAjORcMpIrnkS1tkgAq4DBTdE7JpQ+zUNl/iThFC0j/Lm6eyO+3322Wu0c5TMkvAMkAYsxGr34/t2u6NWJr1vZ+3zw1Ewmks6eLHHcyt2ctv/2q811K1as4JZbbnH4WPVdXqGRXafS2HkyjbzCigNXyW7ujnwJMRaPJL77hQ+sfgkZdNUU/lr7A6cO7SmzJgB4BLgfMM9OUAAsAuah/13p6jJ5b1kSwOxoGOzNyPZRhEnXeOFkTUJ9uaZHY2KrMFjYzHZX/FG4GVZgMvoCC9Dbxh4AKg6WVbk5lSzHxBlzaRDkzTsvPEZqSgqrVq2SAAakZBew82Qq+89lVCpzhv0nJttfQnZvXMXISXcD5cd6pSeXzInuhj4ofh4QVWL5V8BjlJ39qq6T95YlAcwKT3c3BrYKp6t0jRc1KDrIm2t7NmbZX2fIyi+qeIcySraJmbvit+zSlx/f/y+N26xmzZdjyUxxB2ahB7Eb0XsvWlfVm5O1csx+ZwX71yzhuScfqfTnqk9OpeTw18lUjiVlV6nzzqCrpjg083dZ/oHnh2HYDoIjgP8BXUos28KUR0wsnj/J6nGdnSu2umQcGKXHgbWI8OOSdpEEeF84vyRRv6XnFPLdX6dJz636dO9m2RlpzLluACiYtfBPFs9vR9wJcw3CdmAs5xMDl1adGZ5t8XR3o3uTYLo1DcbX8+L4vlxQZOJQfCa7TqWRmFm18aUlGYsK7c78bU3JLBqaphEa1bBEmrzWwMvowzDMTqJXIS6hz+UTbE7CO+mBF2xOonoxZKO/YPl6GhjduQHjujWS4CVqVZCvB9f1bkK4f/VnKyiZy+7MkZX832uniGp6rHhtb2AL0N7qvjWRWbygyMTWYym8v/EYq/bEEZ9hPbVWfciheCYtl5/3xvHexqP8si/eKcELKs6daU3JCXoTTv5Lm+4D6DTgGuAVYC/ng1cW8DhN2k7EN2AloHHoL31Gq4qS914ILo6vRBVoFx1Ix4aBMg+SqDP+Xu5M7NWEZTtt50+0xt4EmcvffZF1335IamIS+jfu29CrEjehz2i0lmbtu5N4+pjNCTidxWjS2H8ug/3nMmgY7E23JiG0jvTHzU2RkpLC4MGD9c/jYjkUM/MK2X8uk31n00nNqf4TtC32pjSpSJfBo3H3uI+j/0Sgz/ANenvoR8ATQBx3PLsZTdPYGPspPUeO4881y8tVTZuzxFxIJIABPZu5zn8YUX95exiY0KMRy3ed5UxqrkP7OJ7L7nZad4/m8M4rgWDc3H6hba/F3PBgazSTVqs3p7NpeZxNO0eAtzsdGwaxbfW3LpVDMT23kKOJWRxNzOZUak6N5Cgsq2RHmcgmzTm+byegAA3fgGByMtOs7tei88MknJrL3xu9SyzdCMwAdlqW7Nu2nj6XTbAM06hslpi6Im1gQlxgCo0mVvx9zqHM4462j5g7aOzeGMqXL0dhLNRbDwZelcbVdydgqMWvsvYyoAcHh9Cw4YWVQ1HTNOIy8jiamM3RpGySnFQ1WBklc2cqN8VvSz+h+/Ax7Fz3Iz1GjOWl268ss0cL9Kfu87kq3T3PUVRwHz7+a5g4c26pBNLdho7mpsdfrVYZ66INTJ7AhLjAeBjcGNu1IT/tOcfh+Cy725onyARsBrGSvQt7DM8kOLyIj59pQFa6O398H8y5o57c/OQ5AkJsd+evSsoiW/vYe2pMS0slLe18DsVl3//AQ7Pvd/iczmA0aSRm5nMuPZe49DxOpebUyNxblVF2ELv55wZTZ5XpcBGAXi04A9A773h6GynMn0tR4f8Yes31XHLdygu+atBREsCEuAAZ3BSjOzXgF0M8+85m2N+2uJHfVgAr2/W5RedcZr15ko+ebsjpw94c3ePLq/c05ZY5Z2natvzTRcn0RY6mLLK3z6CrprB3yzqHMqC7dbic9zceJTLQm6gAL6KDvIkK9HZae3VugZHMvEJScwotASsxM99pGeDLBvHq5i60Rq9edMPT617c3F8gL7vk3G2fcMn18az6RJ8+qkFMa0vwggu3atBREsCEuEC5uSku6xCFp7sbu06m2d3WXiP/b7Gf0qHvsFI3zZDIIu793ym+XRDJjjVBpCV5sPD+JkyckUDvy0oHzJLpixxNWWRvH0cyoJd8aszMKyIzL4t/E84/jXp7GPDzMuDr6Y6fpwEfTwN+Xu74ehpQKEyahqaBUdOKf9YoMmpk5ReRmVdERl4hmXlFNZrxvWwQ1zStSgmUKwp63Yf/lyO73yYrrQkUf/9wM2zDZLwH2MHar86nvzN37DFz5nQ7dUECmBAXMKUUw9tG4mVwY+uxFJvb2Zsgc82Xb7P2q3fL3TQ9vTRueDCexq3z+f6dCIoK3fjy5WgO7czn1OFRKKV3rCjbs3H14reK9/e23AAd6Q1p1vOSqyrIgF7xgNm8QiN5hUaSKbB7nLpUNohrmlbpLwL2nmQTz3jww3sR7Nl0flYMN8NZPLznkJ/9vmVZQd75DkG2Jql0VRLAhHABA1qF4+nuxsbDSVbX28rKseSVx9i/bT0mo9HqTVMpGDI+jQbN8/n0uQZkp7vz56/N0LtY3wAcL7W9rRtgZWb2/eWLtyr8vOZeca6kbBBPS4qz/PzVq0+U2jb27ecdehKy9iSbneHGz5+H8ccPwZiMeqYgg3sBxqJnMRlfIT/bsR6sF1JOw6qSACaEi+gVE4qXu4FfD8SXq3kzN/InnD7GWw/dZFlu70mo5E2zdbdcZi3U28XOHPEG+gG7gLuBJVbL06RNJwaOnczx/bto3KYTXYeMcihbhHkqD2s8vHwozM+t0TFpNcVeENdMpZ8487KzyMs+v635i4C9J9nYd17lxw8hO306mkl/ElNKo9fIDC77Txw/vr+D3b85Hryqm9OwJtrzKksCmBAupHPjIDzd3Vi9Nw6jlY4GlXkSKlt9FBpVxH2vneLH98PZGBsCBAFfApcB/weU7tafdPYkL905lsTTeqaP/1ugB7rK5u0z6zpkFFff/Sh//PClS/aKq2rewpJPQrZ/f9eSlzUfaGlZEt7wX256wp3GrfIBVWFv1JKqm9OwKh17aoKkkhLCxbSNDmBs14Z4GMonmh501RS6DhlV4TFsVR95eGqMn57IbU+fwTfAnFliKvAX0L3UtrlZGZbgBfB77Gd2Ux517HeJ3fLc+MjLBIVHM3rqLPyDnJuTsTaYhzQ4cv3Nyj4Jlf/9DQe2At9wPnjtp3nHOTz0Xl5x8Dp/fkdTTlU3bVjJlGV7Nq2p1rGqQwKYEC6oebgfV3dvhKd76f/CjtxEHak+6tg/m0snvw2Y5/pqg55HcRZ6BojyDmzfyCvTr7Z5zNbd+9tcV5UnguP7d1mqsS4UBncPelei6rPs5zb//lp3vxdYDawF+hSvTQSm02XQg0z/73W4W5ko1F5v1AFjbqhyTsOE08d48Y4rLa/v33vJsm75uy/y4h1XcvUl/ejYsSOHDx+u1LGrQwKYEC6qcYgv1/ZsjI9n6TFRFX0TdzRYnDiwCrgMd8+5QBHgiT79xu9Ah3Lb52ZlkFFirqnel40vlQj2z7Xf2zxXZZ8IzFVYC++fTE6JatG6lp2RxgdP3uXw9mU/d8JpDxa/2JTDO99Ar7oFPeHuM0ArYBHX3veEzd9fyd6o190/j84DL8XHLxCAnMx0Hn5/JSNvuJsJ9z5Zqc9lrto0v0pWRZsTB/976AD79u1j5cqVlTp2dUgAE8KFRQV6c32vJgT6lL6h2fsm7miwmHDvk4y84U6e/Oxy7pi3H/i3eM0A9Dx6T2PO9lBWh37DuW7ms3QbOtpy0wwK0ydMdEaW8wulCqssc3d5AE9vH2564jViOpyveo3p0N3q506Jc+frVyN56fYYdm0IKN66AHgDvepwDqCPz7P3+9N/Z3fz8Psr6TdqIlPnLOSRD1dZgpZ54HJlq2gdrZqeOHEi06dPr9Sxq0NyIQpRD2TlF7Fs5xlLnr5Pn5vFrg0rrY4Lq0reu22rv2PJK/PQg9YswPzUdwC4A/2pTGerirJkPj9zV//M1GRLKiN7N1V7+RNLzn0Feu/K7Iw0oOZ7yNkrl7dfAEHhUZiMRrLTU+jQbzhX3fkw/kGhls/daeCdbP6xBdt/CbR0iVdKwz/4ZzJTp+EbkOqU358zGIsK+fKlB/lr/U9W10+cOJHFixfjYaVq0xGSC1GIi5S/lzsTezbm+11nOZOWa3NcWFXz3ulVUzn4BjzH8Osa8cf3w0hLbAC0Q89u/g7wMJBus4qybD4/cDyVUWV6V+7+bRWrPnm9VnrI2StXXnYmedmZlveNW3WwBOm8nGjSk+fz+n2BmEzn2xQ79M1i2MRjLHpwNCiN+xb8RGTj5kDd5i308zLQs00EN3zxIY1LJFsu6c0336xy8KoqqUIUop4wT8fSIsLPEiys5b1zpPqobAeJklVTI67vyWU3fgI8COQUb3EXcAi4l382/+bET6WrTO9Kv6CQWqteHHTVFFp16+dQuQaOnUz8SU++eCma+bfFsP3nIEvw6tQ/i1kLT3D7s2dJPvuDnqlE0zi+9y/LMapa/VcdAd7uDG8Xya0Dm9OzWQg/r7L+9AWwYsWKWiuXmVQhClHPmEwavx5IYM+ZqnVuyM5IY851A+w+wZirKL39uhIS+T3njjW1rPP0PsfEmdB9WCZuTvyKXNHUMV6+/gSHR5GRmmS3etHR3H+ODNTNzkjjqYn90dBs5nXsMngUfS9fxMbYCA7s8Cu1rvOgTHqN/IeVH51/qqqoerQ2chcG+XjQp3ko7RsEYnA7/4Q4adIkvvrqK8LCwli0aBGapjF9+nSSk5O5/vrrWbLE+qB3R0gVohACNzfFpR2i8PU0sM1O/kRbHEneW7KK0j84j92/neXH90NJifemIK8Bi+fD+m9CuPK2JNr2zEFZ73lfKRVl3c/PySL+ZOnpZ6qa+8/Rgbp6pw1beR29gCnEHX+D9544nyFeKY0ug7O4dEoyDZsXsGHpj1UefO5sAd7u9GsRRocGgbi5lf+lLVy4kFatWnHfffcRGRkJwLBhw1iwYAGzZtV+Vnt5AhOiHvv7dBrrDiRisvP/vLIdJGw9ARiLYMvKIH5eHEZm6vnvxq265nD5f5Jp0Tm32oFM70zyWJX2bdWtH3c9/75DQwi2rvqWr/73OACTZj9vCeL2rtV5DYE7gWlApGWpp5eJ3pdlMHh8KpGNCy3LKzspaXUyaNji42mgd0woXRsH4W6om5YleQITQpTSpXEwgd4erPjnnM2pQ6qTfqokgzsMvCqdXpdm8NvSENZ+E0J+joEju305stuXRi3zGHx1Gt2HZ+LhWbUvzvay7gdHRJOWGGdjT8XRv7eTn5tj9WnK0Wz62RlpZKVaS6jsBYxDz1pyKed7aYKHZxKX/Uej3xXp+AWW/x1UdlJSZ/J0d6N7k2B6xoTg5e6cOdZqU60GMKVUAPAk0A09L0048LSmaXNrsxxCXExiwv2Y2EvvoZiZV1RuvaM5/BzNXu7lo3HplBT6X5nOmi9D+ePHIIyFbpz515slr0Tzw/vh9B+dzoAx6QRHlC+PPfZ6V/YYMZaXbr/Sxp76VCa2qkQrE8SjY1oTd9ycbaIXetC6AQgptV/DFpmERn3DNfd1JCis9LqyzBk8HJ2UtLoMborOjYLo0zwUPy/XfY6p1SpEpVQMeorr3ehdlm6nGgFMqhCFcFxmXiHLd50lMbP8rMsVVWNV5wkgM9XAlp+C+OOHYDKSz98s3dw0Og/KoueIDNr0yMHTq/y9qDIZzx2pXrRVJepoNR64MWjcAn5fngaMB8red1IIa/g7Ux7qSLP2RoerTPXOIP2wdT+e9MALTsvO3zLSn8Gtwgnx83TK8ZzFFaoQTwAhmqZpSqlw9AAmhKgFAd4eTOzVmJ/+ieNYUunM8hV1kKjOE0BAiJFLJ6dwyXUp/L0xgI3Lgzm+zweTSbH7twB2/xaAp7eJ9r2z6TQwiw59s/HxM1U647m5etHDy9vmlC22qkTtV+P5AiOBscAYfl8eXWa9ET1v4cfA98x4bT3+wUbHLk6xshk8Jj3wQqnqUWdMLxMV6M2QNuE0DvGteGMXUasBTHPlHiNC1ANe7gau6tqQ9YcS2H2qdOeDitJPVfcGanCH7sMz6T48k5MHvdgYG8Lujf4UFbhRkOfG7o0B7N4YgMFdo3W3HPyCjmIydQb2ODSDsbl6ccCYySx/5/kKn6Y69h9Rqkr0fBDfhp5At1/xawjgY+UIO4Bvgc+As5aljlyrhNPHePcx/fu7h5d3qTY3N4M7qz9bCIBfUDDdho5i1M0z7B7PngBvdwa0DKd9gwCUM7qDXkBct/JTCFElbm6KS9pFEeHvzbqDCZZ5xex1kHD2BJNN2+Yz5eE4rr1PcfBPP/7+3Z+9m33Iz/XAWKSKx0uNKn7l8dWr/7D83S14eO7Bw/Mfbn/2fqKaNit1zJKZPhyZG6tdzxGkxPuSluBO3AkvThzw5tBf4UCyjT1y0bPz/wD8SMmgBXoANBYVOnStdv+2ipS401bXlc3gMXDs5CoNXvZ0d6NXsxB6NAvBo456FtY0lwlgSqlIIKLM4pbWthVCVKxz4yDC/D1Z8fc5svKLnJ5+yhFePhpdBmXRZVAW6775jB/e24betnQ1YK6q80Yz9SY3szfm+YZfvMNEQLAJ/xAjgSFFBIQUERBqJCC4CDcDmIyKiMZvA++i3+bc0asCGwPNgGZ8t9B6SqTSDqDnefwBWMP5zCOlVXYyTkcDkqMdZ0pSCjo0CGRAq3D8XbiDhiOq3IlDKTUMWOfg5t01TdtVZv9w9AluHOrEoZSai56SuRzpxCFE1WXlF/Hj7rOcS7feblRbSnekUEBb9F5+vYtf3bBeleccnt7ZxHQooKhwIwPHtqJVV8XSN+93aIbjZ77ejH9w1ZIR21KVjjONQ3wY2iaCyEBvh/e5UNR2J46D6GmoHXGyGucxewt9WtKSWgLLnXBsIS5a/l7uTOzVhLXVSD/lDOU7Uhwofn1evIU70AmlejP0mrnkZvuRmWIgM82djGRFepJCn7PMGiP69CRn0G9HJ0r8e4J2vZtw69wHiyeJbGfZy17HlpIqavey103flsp0nAn29WBw63BaRQZUvHE9UuUApmnaOeB9J5alovMlAAkll9W3Bkkh6orBTRGQcYzwnExS/WMs7WK1Xg67vSGLgF1o2i6imzUv1anjfOaMYCbc8xS9Lh1DblYKf/zwKQPHTOL792x36tCfdB7E4O5Rrtt+RTMc79rwk0NthI6OtSvJkc4gXh5u9G0eSrcmIaVyFl4s6ncFqRDCISkpKQwePBiAv4+cYNOpfDJyCyvYq2bYCxpmy95+nm9en2t57+FpfvJK46dPn+aPH9+yrOs7apxDQwSsddsv2bFlzB0PsH/rBo7s2kpudoZlhmNH2ggdybbRrH03Bo+7kaVvzqswKBrcFJ0bB9GveVi5GbkvJrUewJRSVwB+gPlZt4NS6trin1dqmma9lVQIUWNiY2MpKtKzYmxZu5opN97E2gMJHIzLrGBP53NkPFd+dumkvcbCAsvP1sZ6efv62zyf+UnHWhLjsh1b+o2aWG4STkfmM4OKx9rd9vQi/INDad19gN2g2CrSn0EX4EDkulDryXyVUsfRuwFZ01zTtOOVOJZk4hCiCg4dOsT48eMt78+dO0dqaioAISEhNGig99DLLzIx5Yk3CI5uavU4NWH/9t/YtWElV9w8y6HxXPaYO0IsfvGhcjNUf/3aU+RlZ+Lp7UNIVKNamcbEXraQirJt1MeByCW5QiYONE2Lqe1zCiFKW7lyJfv27bO6LjU11RLMALzP7SayTWsSMsqnoHK27Iw0PnjyblAw7q5HK6x2CwqPJj3JegLfkr34rA0ROP3vAdYueYeCvNxyHSxqahqTqoy1C/b1oH/LMNpG1b+ByNUlbWBCXITuueceNm3axDfflO3YW9rEiRN5YNZ9uBnc+f1IEjtPptqat9EprFbj2al2u+uFD3jpDusJfEv24rNW1RcW1cihMlVlLJYtlRlr5+dloG/zMDo1CrooO2g4QuYDE+IiVVhYyJQpU2wGsYkTJ7J48WI8PM535T6VksOv++NJzXFOBw9H5iLLzcogIznB6v59Lr+Gbau/s7qubJVcbY3Fqi4vDzd6NQule9PgeptBwxqXqEIUQlwYPDw8WLhwoc0A9uabb5YKXgBNQn25sV8zth5L4c8TqdXubl+ZaUw8vHy44cHSSW4P/bUJcKxKrqbHYlWXh0HRtUkwvWNC8fa4eHsWVoYEMCEuYitWrLC77pZbbim33N3gxsBW4bSJCuDX/fHVyuDh6PioyCYtmPbiRwSF6+mlzNVuPUeO4881yx2qkqupsVjV5enuRudGQfRsFuLSc3PVhYvn+VQIUc7q1asBCAsL4+uvv+arr74iLEwPBKtW2b/RRwR4cX3vJgxrG4Gne9VuJebxUV2HjLK5Tdcho3jwne8twQsgICSM0VNnEdWkBaOnzrIEr5LryuYbdORczdp348ZHXsY3IBiAA9s3VulzOcLT3Y3eMaFMHRjDkDYREryqQNrAhLiIJSUl8dprr3HfffcRGRkJQEJCAgsWLGDWrFmEh4c7dJyMvEI2HUniQFxmlTp5ZKYmM+f6AVbXVZRnsCbOVXaslzN5urvRvUkw3ZuGXNSDkMuqShuYPIEJcRELDw9n3rx5luAFEBkZyXPPPedw8AII9PZgVKcGTO7blObhfpUuR0VzkTmTI+ey9RRXHQHe7gxuHc5tg5ozoFW4BC8nkAAmhHCayABvru7eiGt7NqZBkOMZ0UuOj7rpide46fFXa6warzbPBdAgyJvRnRtw68Dm9JIOGk4lla5CCKdrEurLpD5NOZKQyR9HkknJLrC7fW3ORVYb53JTilaR/nRvGkzD4JqbAuZiJ21gQogapWkaR5Oy2XUyjZMp9TvVaYivB+0bBNKhYSAB3rU3dqw+kHFgQogLjlKKlhH+tIzwJykrn92n0th/LoNCo+t+eS7J092NNlEBdGgYSCN52qpVEsCEELUm3N+LEe2jGNgqnL1n09l9Kp30Opq2pTrc3RSNQ31oFx1Iq0j/iypjxoVEApgQolZs2bIFgH79+uHtYaBns1B6NgslISOPI4lZ/JuQRVKW/bayuuTnZSAmzI8WEf40DfWt8tg34TwSwIQQNa7khJkJCQmEhIRY1kUGehMZ6M2AluGkZhdYgllcRl6NJg6uiLubIjLQiyYhvjSP8CM60FuywV9gJIAJIWpcyQkzY2NjmTp1qtXtQvw86e0XSu+YUPIKjcRn5BGfkV/8bx6ZeUU1Uj6lIMTXk+ggb6IDvYkO8ibc30uywF/gJIAJIZzO2oSZZrNnz+bll1+2vI+NjaV169bljuHtYaBZmB/Nws4PjM7OLyI+I4+03EJy8o1kFxSRU1BEdr6RnIIicgqMaJoekNyUwk3pnUjclMLdTeHn5U6gjzsB3h4EeLsT6K3/HOTjIeOzXJAEMCGE01VmwsyVK1cyY8YMh47r5+VOiwh/m+vNw4Kkqu/iIK2QQginu+eee5g4cWKF202cOJHp06c77bxKKQleFxEJYEIIp/Pw8GDx4sV2g5i1CTOFqAwJYEKIGmGeMNMWaxNmClEZEsCEEDWmogkzhagOCWBCiBpTnQkzhaiI9EIUQtSYhQsX0qpVq1ITZg4bNswyYaYQ1SHZ6IUQQtQ5mZFZCCHERUMCmBBCCJckAUwIIYRLkgAmhBDCJUkAE0II4ZIkgAkhhHBJEsCEEEK4JAlgQgghXFKtBjCl1CVKqQ+VUgeUUtlKqTNKqeVKqZ61WQ4hhBCur7afwKYBMcACYDQwA4gEtiilLqnlsgghhHBhtZ0L8R5N0xJKLlBKrQKOAI8Ba2u5PEIIIVxUrT6BlQ1excuygH1Ak9osixBCCNdW5504lFJBQA/AoeSNQgghBFwY06m8CfgBz9nbSCkVCUSUWdyypgolhBDiwlblAKaUGgasc3Dz7pqm7bJyjGeBKcD/aZr2ZwXHmA7MqUQRhRBC1GPVeQI7CNzh4LYnyy5QSs0BngAe1zRtoQPHeAv4psyylsByB8sghBCiHqlyANM07RzwflX2LQ5ec4G5mqY97+D5EoCyPRircnohhBD1QK134lBKPYkevOZpmvZ0bZ9fCCFE/VCrnTiUUrOBZ4BVwAqlVL+S6zVN21Kb5RFCCOG6arsX4tjif0cVv8qSOkEhhBAOqdUApmnasNo8nxBCiPqrzgcyCyGEEFUhAUwIIYRLkgAmhBDCJUkAE0II4ZIkgAkhhHBJEsCEEEK4JAlgQgghXJIEMCGEEC5JApgQQgiXJAFMCCGES5IAJoQQwiVJABNCCOGSJIAJIYRwSRLAhBBCuCQJYEIIIVySBDAhhBAuSQKYEEIIlyQBTAghhEuSACaEEMIlSQATQgjhkiSACSGEcEkSwIQQQrgkCWBCCCFckgQwIYQQLkkCmBBCCJckAUwIIYRLkgAmhBDCJUkAE0II4ZIkgAkhhHBJEsCEEEK4JAlgQgghXJIEMCGEEC5JApgQQgiXVKsBTCnVTSm1Qil1UimVq5RKUUptVkrdWJvlEEII4frca/l8wcAp4EvgDOAHTAE+U0rFaJo2r5bLI4QQwkXVagDTNG09sL7M4h+VUs2BOwEJYEIIIRxyobSBJQFFdV0IIYQQrqO2qxABUEq5oQfPEGAicDlwbwX7RAIRZRa3Azhy5EgNlFIIIURtKXEf93R0H6VpWs2Uxt5JlXobuKv4bQEwU9O0RRXsMxeYU8NFE0IIUbfGaZr2vSMbVjmAKaWGAesc3Ly7pmm7SuzbFIgsfo1Fb/96WNO0l+2cz9oTmD/QBtiDHgiroiWwHBgH/FvFY9Qncj3Ok2tRmlyP0uR6nOeMa+EJNAE2aJqW7sgO1alCPAjc4eC2J0u+0TTtZIllK5VSAC8opT7RNC3R2gE0TUsAEqys2upgGawqPjfAv5qm7a3OseoDuR7nybUoTa5HaXI9znPitdhZmY2rHMA0TTsHvF/V/cvYBtwNtACsBjAhhBCipAulF+JwwAQcreuCCCGEcA212gtRKfUukIH+xBUPhKP3Qrwe+K+t6kMhhBCirNruRr8ZmArcjJ6VIwvYDfxH07TPa7ksZonA00jVpZlcj/PkWpQm16M0uR7n1cm1qJNu9EIIIUR1XShtYEIIIUSlSAATQgjhkiSACSGEcEkSwIQQQrikehvAlFL+SqnXlFJnlVJ5SqldSqlJDu4bqZT6WCmVpJTKKZ50c0RNl7kmVfV6KKUmKKW+VEodKZ6E9LhSarFSqnVtlLsmVOdvo8xx5imlNKXUnpooZ22p7vVQSo1TSm1QSmUopbKVUnuVUnfWZJlrUjXvHcOVUr8opRKUUllKqb+VUvcppQw1Xe6aoJQKUEq9pJT6WSmVWPz3PrcS+9fsvVTTtHr5An4GUtGTBg8H3gM0YHIF+3kB/6BPvDkFuBSIBQqBoXX9uergemxFz3E2FRgK3AjsAzKBjnX9uWrzWpQ5RjcgD4gD9tT1Z6qr6wE8AhiBN4FRwAjgHuDeuv5ctX09gJHF12Idek7AkcDrxfsuqOvPVcVrEQOkARtKXIe5Du5b4/fSOr9ANXTRRxdf6BvKLP8ZfSZog519pxfv27/EMndgL7C1rj9bHVyPSCvLGqInT36/rj9bbV6LMn8PO4EF6BO0umwAq+bfRs/iG/ZDdf05LpDr8Xnxlxq/MstXA+l1/dmqeD0U54dbhVcygNX4vbS+ViGORx8k/U2Z5R+h33z7VrDvQU3TNpsXaJpWhP7H2Ucp1cjJZa0NVb4emp5Eueyys8Bp9MzRrqY6fxtmjwChwOPOLVqdqM71uBfIB96omaLViepcj0L0L3a5ZZanoQc2l6MVq+LuNX4vra8BrBOwv/hilfR3ifX29v3bynLzso7VLFtdqM71KEcp1QJohv5NytVU61oopToATwDTNE3LqoHy1bbqXI8hwH7gGqXUQaWUUSl1Wik1Xynl8KSEF5jqXI+30acEeV0p1VApFayU+g/6jfwl5xf1glfj99L6GsDCgBQry1NKrK+JfS9UTvtMSil34AP0b6mvVr9ota7K16J4JvEPgaWapq2sgbLVher8bTQCWqO387yO3ubzMfAA+hOLK6ry9dA0bStwCXrAOoPejvYR8Limaa84uZyuoMbvpbWdC7E22XvsreiRuDr7Xqiq/ZmUPunPB8Bg4BpN0045o2B1oKrX4n70G/ZVzi1Onavq9XADAtDbi5YUL1unlPIDZiql5miadsT27hesKl0PpVRPYBl6x6e7gGz0gDZPKeWtadqzTi2la6jRe2l9DWDJWI/uocX/WvtW4Ix9L1TV/kzFwet99F6IN2uattx5xatVVboWxbOIP4Pe/lWglAouXuUOuBW/z9c0rWz7x4Wuuv9XotE7KZT0EzAT6AG4WgCrzvV4E32WjfGaphmLl61TSpmAuUqpxZqmXUxTRtX4vbS+ViH+A7Qvru4qqXPxv/bG7fxTYrvK7nuhqs71KBm8pgK3a3U3c4AzVPVatAB80HseppZ4DQTaF//8gtNLW/Oq87dhrX0D9J5roM/x52qqcz26AX+WCF5m29Hvte2dUkLXUeP30voawJYB/sA1ZZbfDJxFf8S3t287pZSlt1HxH/ON6F0/zzq5rLWhytejOHi9hx687tI0zVXbNsyqei12oY8JKvvaDRwv/nmh00tb86rzf+W74n+vKLN8NHrw2u6MAtay6lyPs0AvK4OW+xf/e9opJXQdNX8vretxBjU4fuFn9EfUO9BvLu+i17lOKbHNB0AR0KzEMi/0bwYngcnoDdNLqR8DmatyPd4o3u4DoF+ZV/e6/ly1eS1sHGs9LjwOrJp/Gx7An+jdxO8r/r8yv3i7N+r6c9XB9fi/4u1Wog9kvrT4ehQCv9T156rG9bgCuBb9S6wGfF38/lrA1871qPF7aZ1fnBq86P7o1T3n0Meq7AYmldnm4+JfSEyZ5VHAJ+h1uLnoE3GOrOvPVBfXA/3pQrPxOl7Xn6u2/zasHKs+BLDq/F8JRe8+Hoc+Buogei9Et7r+XHV0PSYAG9EndswqvoE/QZnBza70quAeEFPB9ajRe6lMaCmEEMIl1dc2MCGEEPWcBDAhhBAuSQKYEEIIlyQBTAghhEuSACaEEMIlSQATQgjhkiSACSGEcEkSwIQQQrgkCWBCCCFckgQwIYQQLkkCmBBCCJckAUwIIYRLkgAmhBDCJUkAE0II4ZL+H2KZvaxoZDCCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 480x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    # Initialize plot\n",
    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "    # Get upper and lower confidence bounds\n",
    "    lower, upper = observed_pred.confidence_region()\n",
    "    # Plot training data as black stars\n",
    "    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')\n",
    "    # Plot predictive means as blue line\n",
    "    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')\n",
    "    # Shade between the lower and upper confidence bounds\n",
    "    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c31ad53f-c573-46ce-9d6b-d9b8de49a874",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}