{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hadamard Multitask GP Regression\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This notebook demonstrates how to perform \"Hadamard\" multitask regression. \n",
    "This differs from the [multitask gp regression example notebook](./Multitask_GP_Regression.ipynb) in one key way:\n",
    "\n",
    "- Here, we assume that we have observations for **one task per input**. For each input, we specify the task of the input that we observe. (The kernel that we learn is expressed as a Hadamard product of an input kernel and a task kernel)\n",
    "- In the other notebook, we assume that we observe all tasks per input. (The kernel in that notebook is the Kronecker product of an input kernel and a task kernel).\n",
    "\n",
    "Multitask regression, first introduced in [this paper](https://papers.nips.cc/paper/3189-multi-task-gaussian-process-prediction.pdf) learns similarities in the outputs simultaneously. It's useful when you are performing regression on multiple functions that share the same inputs, especially if they have similarities (such as being sinusodial).\n",
    "\n",
    "Given inputs $x$ and $x'$, and tasks $i$ and $j$, the covariance between two datapoints and two tasks is given by\n",
    "\n",
    "$$  k([x, i], [x', j]) = k_\\text{inputs}(x, x') * k_\\text{tasks}(i, j)\n",
    "$$\n",
    "\n",
    "where $k_\\text{inputs}$ is a standard kernel (e.g. RBF) that operates on the inputs.\n",
    "$k_\\text{task}$ is a special kernel - the `IndexKernel` - which is a lookup table containing inter-task covariance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "In the next cell, we set up the training data for this example. For each task we'll be using 50 random points on [0,1), which we evaluate the function on and add Gaussian noise to get the training labels. Note that different inputs are used for each task.\n",
    "\n",
    "We'll have two functions - a sine function (y1) and a cosine function (y2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x1 = torch.rand(2000)\n",
    "train_x2 = torch.rand(2000)\n",
    "\n",
    "train_y1 = torch.sin(train_x1 * (2 * math.pi)) + torch.randn(train_x1.size()) * 0.2\n",
    "train_y2 = torch.cos(train_x2 * (2 * math.pi)) + torch.randn(train_x2.size()) * 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up a Hadamard multitask model\n",
    "\n",
    "The model should be somewhat similar to the `ExactGP` model in the [simple regression example](../01_Exact_GPs/Simple_GP_Regression.ipynb).\n",
    "\n",
    "The differences:\n",
    "\n",
    "1. The model takes two input: the inputs (x) and indices. The indices indicate which task the observation is for.\n",
    "2. Rather than just using a RBFKernel, we're using that in conjunction with a IndexKernel.\n",
    "3. We don't use a ScaleKernel, since the IndexKernel will do some scaling for us. (This way we're not overparameterizing the kernel.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MultitaskGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(MultitaskGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.RBFKernel()\n",
    "        \n",
    "        # We learn an IndexKernel for 2 tasks\n",
    "        # (so we'll actually learn 2x2=4 tasks with correlations)\n",
    "        self.task_covar_module = gpytorch.kernels.IndexKernel(num_tasks=2, rank=1)\n",
    "\n",
    "    def forward(self,x,i):\n",
    "        mean_x = self.mean_module(x)\n",
    "        \n",
    "        # Get input-input covariance\n",
    "        covar_x = self.covar_module(x)\n",
    "        # Get task-task covariance\n",
    "        covar_i = self.task_covar_module(i)\n",
    "        # Multiply the two together to get the covariance we want\n",
    "        covar = covar_x.mul(covar_i)\n",
    "        \n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar)\n",
    "\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "\n",
    "train_i_task1 = torch.full((train_x1.shape[0],1), dtype=torch.long, fill_value=0)\n",
    "train_i_task2 = torch.full((train_x2.shape[0],1), dtype=torch.long, fill_value=1)\n",
    "\n",
    "full_train_x = torch.cat([train_x1, train_x2])\n",
    "full_train_i = torch.cat([train_i_task1, train_i_task2])\n",
    "full_train_y = torch.cat([train_y1, train_y2])\n",
    "\n",
    "# Here we have two iterms that we're passing in as train_inputs\n",
    "model = MultitaskGPModel((full_train_x, full_train_i), full_train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training the model\n",
    "\n",
    "In the next cell, we handle using Type-II MLE to train the hyperparameters of the Gaussian process.\n",
    "\n",
    "See the [simple regression example](../01_Exact_GPs/Simple_GP_Regression.ipynb) for more info on this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 1.000\n",
      "Iter 2/50 - Loss: 0.960\n",
      "Iter 3/50 - Loss: 0.918\n",
      "Iter 4/50 - Loss: 0.874\n",
      "Iter 5/50 - Loss: 0.831\n",
      "Iter 6/50 - Loss: 0.788\n",
      "Iter 7/50 - Loss: 0.745\n",
      "Iter 8/50 - Loss: 0.702\n",
      "Iter 9/50 - Loss: 0.659\n",
      "Iter 10/50 - Loss: 0.618\n",
      "Iter 11/50 - Loss: 0.580\n",
      "Iter 12/50 - Loss: 0.545\n",
      "Iter 13/50 - Loss: 0.511\n",
      "Iter 14/50 - Loss: 0.478\n",
      "Iter 15/50 - Loss: 0.445\n",
      "Iter 16/50 - Loss: 0.412\n",
      "Iter 17/50 - Loss: 0.378\n",
      "Iter 18/50 - Loss: 0.344\n",
      "Iter 19/50 - Loss: 0.309\n",
      "Iter 20/50 - Loss: 0.276\n",
      "Iter 21/50 - Loss: 0.243\n",
      "Iter 22/50 - Loss: 0.211\n",
      "Iter 23/50 - Loss: 0.182\n",
      "Iter 24/50 - Loss: 0.154\n",
      "Iter 25/50 - Loss: 0.129\n",
      "Iter 26/50 - Loss: 0.105\n",
      "Iter 27/50 - Loss: 0.085\n",
      "Iter 28/50 - Loss: 0.067\n",
      "Iter 29/50 - Loss: 0.052\n",
      "Iter 30/50 - Loss: 0.039\n",
      "Iter 31/50 - Loss: 0.029\n",
      "Iter 32/50 - Loss: 0.021\n",
      "Iter 33/50 - Loss: 0.015\n",
      "Iter 34/50 - Loss: 0.012\n",
      "Iter 35/50 - Loss: 0.010\n",
      "Iter 36/50 - Loss: 0.009\n",
      "Iter 37/50 - Loss: 0.009\n",
      "Iter 38/50 - Loss: 0.009\n",
      "Iter 39/50 - Loss: 0.009\n",
      "Iter 40/50 - Loss: 0.008\n",
      "Iter 41/50 - Loss: 0.008\n",
      "Iter 42/50 - Loss: 0.009\n",
      "Iter 43/50 - Loss: 0.009\n",
      "Iter 44/50 - Loss: 0.010\n",
      "Iter 45/50 - Loss: 0.010\n",
      "Iter 46/50 - Loss: 0.010\n",
      "Iter 47/50 - Loss: 0.010\n",
      "Iter 48/50 - Loss: 0.008\n",
      "Iter 49/50 - Loss: 0.007\n",
      "Iter 50/50 - Loss: 0.005\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iterations = 2 if smoke_test else 50\n",
    "\n",
    "\n",
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "for i in range(training_iterations):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(full_train_x, full_train_i)\n",
    "    loss = -mll(output, full_train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/50 - Loss: %.3f' % (i + 1, loss.item()))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make predictions with the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAADSCAYAAACW5MO6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABWg0lEQVR4nO2dd3xUVfbAv3dKMukhvdM7hCpFEBA7Iig2sKxtRbCs/iyoa8Oyu3a3iAV3V1FBsCC6IgqoSBOk9x4CpJDek0mm3N8fMwmTZCYZkkkmCff7+eSTmffe3Hfee/e8c8u55wgpJQqFQqFQKLyDxtsCKBQKhUJxLqMMsUKhUCgUXkQZYoVCoVAovIgyxAqFQqFQeBFliBUKhUKh8CLKECsUCoVC4UU6pCEWQswVQnzqbTnOBiHE7UKI9S1Q7gQhRJqny3XjvJcJIZa5cdzNQoiVDt+lEKJHE85X8zshxHtCiGfsn1vt+oUQqUKIi+2fHxBCvNIa5+3IKF2uVa7S5Q6qy+3SENsr+h4hRLkQ4rQQ4l0hRKi35WoJhBAGIUShEGKik31vCSG+9IZcbvAX4OXqL66UUkq5UEp5qSdPLKWcJaV80ZNlNoEPgJuFEFFelqNNo3S5Zp/SZSecK7rc7gyxEOIR4BXgMSAEGAV0BlYJIXxaUQ5da5xHSmkElgB/qHN+LTADWNAacpwNQojzgBAp5SZvy+It7M9tBXWem+IMSpdrzq90uQ3TGrrcrgyxECIYeB54QEr5g5TSJKVMBW4AugC3OBxuEEIsEUKUCCG2CyEGOZTzuBAi3b7vkBDiIvt2jRDiCSHEMSFEnhDicyFEmH1fF3tL8C4hxEngZyHECiHE/XVk3CWEmGb/3EcIsUoIkW8/zw0Ox4ULIb4VQhQLIX4Hujdw6QuAa4UQ/g7bLsP2/FYIIe4QQhywX0+KEOKeBu5hrdasEOIjIcRLDt8nCyF22lvuG4UQyY3dNydcAfzawPU4yuNyGE8IMVYIcUoIMcH+/U77dRYIIX4UQnR28bta12Tf9ogQIlsIkSmEuMNhe4gQ4mMhRI4Q4oQQ4mkhhMa+T2P/fsL+24+FECEOv73Vvi9PCPGUE1HWAFe6cx/ONZQuK11WuuyAlLLd/AGXA2ZA52TfAuAz++e5gAm4DtADjwLH7Z97A6eAOPuxXYDu9s8PApuABMAXeN+hzC6ABD4GAgA/bC2kDQ4y9AMK7b8NsJ/nDkAHDAFygX72YxcDn9uPGwCkA+sbuPbDwC0O3z8D/m7/fCU25RfAeKAcGGrfNwFIc/idBHo4fP8IeMn+eQiQDYwEtMBtQKr9elzeNyeyfgE8VmdbrfM6bL/d8bqrj7M/61PACPv2qcBRoK/9fj4NbHRWfp1rmoCtzrxgf/6T7Penk33/x8A3QJD9mg4Dd9n33Wk/ZzcgEFgKfOLwrEuBcfb786b9PBc7yDQUyPe23rTFP5QuK11Wunzm3nlbIc9SeW8BTrvY9zKwykF5Nzns0wCZwAX2ipENXAzo65RxALjI4XsstpeAjjPK281hfxBQBnS2f/8L8F/75xuBdXXKfx94zq4YJqCPw76/0rDyPg2stH8OtlfAIS6OXQY86FB53VXed4EX65R1CNsLweV9c3L+VcCsOtvORnmfBE4AAxy2r6hWKodnWu5w7xtS3gocXvj26xhlfw5V2F+o9n33AGvsn38C7nXY19uhPjwLLHbYF2Avy1F5ewIWb+tNW/xD6bLS5drP9JzW5XY1NI2tFRohnM/pxNr3V3Oq+oOU0gqkYWsBHgUewqbg2UKIxUKIOPuhnYGv7UM5hdiU2QJEuyi3BFgOTLdvmgEsdChrZHVZ9vJuBmKASGwVoKYsbJW1IT4BLrTLeh1wTEq5A0AIcYUQYpN92KwQW0sxopHynNEZeKSOzIk0ft/qUoDtxdZUHgI+l1LurSPbPxzkysfWa4h3o7w8KaXZ4Xs5tlZxBLaWteO9P+FQZpyTfTps9SGO2nWhDMirc94goMgN+c5FlC4rXVa6bKe9GeLfgEpgmuNGIUQgtrmMnxw2Jzrs12AbosoAkFIuklKOxVYhJDaHEbA9jCuklKEOfwYpZbpDubKOTJ8BM4QQowED8ItDWb/WKStQSjkbyME29JHoUE5SQxcupTwBrMPWk7gVu2OHEMIX+Ap4HYiWUoYC32Or2M4oBxznp2IcPp8C/lJHZn8p5Wd2GVzdt7rsBno1dD2NcD1wtRDiwTqy3VNHNj8p5cZmnCcXW6u4s8O2JGxDi2CrL3X3mYEsbL0yxzrmD4TXKb8vsKsZ8nVklC4rXVa6bKddGWIpZRE2B49/CSEuF0LohRBdsM3PpGFraVYzTAgxzd7ifgib0m8SQvQWQky0V3ojtqEOq/037wF/qXYcEEJECiGmNiLW99ge8AvAEnuLHeA7oJfdCUBv/ztPCNFXSmnBNkcxVwjhL4Toh20OpzEWAPcDYzjTWvfBNq+RA5iFEFcADS0h2AncJITQCiEuxzZUVc0HwCwhxEhhI0AIcaUQIqiR++bsnox3st1H2JZwVP9pXfw+A7gIeFAIMdu+7T3gSSFEf6hxzLi+getsFPtz+BzbMw+yP/eHgep1q58B/yeE6Go3EH/F9ozNwJfAZGFzQvHB9vzr6tN4bMNwijooXVa6rHS59gW0uz/gLmAvtgqUhW2+ppPD/rn2m7sEKAF2cMbhIRn43b49H5uSVTstaLA9vEP2/ceAv9r3dcHWcnTmXPIf+77z6mzvjW24KwfbUMfPwGD7vkj7uYvt8rxIA/NK9t8EYnMqWFFn+332+1CI7QW2mNrzKo7zSsOBffbr+wRbBX3JYf/lwBZ7WZnYnDWCGrpvLmTdAox0+C6d/P0RFw4e9s9dsQ0h/dH+/VZgj/2encI+h+fkdx+5un77tlTs8z9AJ2zKmmMv81lA41AfnrVvz7Ef51jPbgNO2p/tU3XKNWAzKNHe1pe2/IfSZaXLSpcR9hMpFB5FCHEpNueIq70tizcQQjwAJEop53hbFoWiOShdbnldVoZYoVAoFAov0uw5Yvv8wO/Ctvh9nxDieU8IplAoWh+lzwpF69PsHrEQQgABUspSIYQeWI9t3ds5GxJNoWivKH1WKFqfZsdYlTZLXmr/qrf/qfFuhaIdovRZoWh9PLJ8ye4+vxNblJNVUsrNnihXoVC0PkqfFYrWxSNZR6RtDddgYUtf9rUQYoCsHUkFIcRMYCZAQEDAsD59+nji1ApFh2bbtm25UsrI1jxnY/qsdFmhOHsa0mWPe00LIZ4FyqWUr7s6Zvjw4XLr1q0ePa9C0RERQmyTUg734vkb1GelywqFezSky57wmo60t5wRQvgBlwAHm1uuQqFofZQ+KxStjyeGpmOBBfYQZxpsAb6/80C5CoWi9VH6rFC0Mp7wmt6NLfelQqFo5yh9VihaH484aynaDiaTibS0NIxGo7dFUZwFBoOBhIQE9Hq9t0VpNqoOtk86Uh1sbyhD3MFIS0sjKCiILl26YIvNoGjrSCnJy8sjLS2Nrl27elucZqPqYPujo9XB9ka7SoOoaByj0Uh4eLh6AbYjhBCEh4d3mB6kqoPtj45WB9sbyhB3QNQLsP3R0Z5ZR7uecwH1zLyHMsQKj5OWlsbUqVPp2bMn3bt358EHH6SqqgqAjz76iPvvv9/LEtYnMDDQ6XatVsvgwYPp378/gwYN4o033sBqdZVD3UZqaiqLFi1qCTEVbqLqoKqD7QlliBVkZmYyfvx4Tp8+3eyypJRMmzaNq6++miNHjnD48GFKS0t56qmnPCCpc8xmc4uV7efnx86dO9m3bx+rVq1ixYoVPP98wwmJ1Evw7FF10DWqDp4DSClb/W/YsGFS0TLs37//rH8ze/ZsqdFo5OzZs5t9/tWrV8sLLrig1raioiIZFhYmy8rK5IcffiinTJkix48fL3v06CHnzp0rpZSytLRUTpo0SSYnJ8v+/fvLxYsXSyml3Lp1qxw3bpwcOnSovPTSS2VGRoaUUsrx48fLBx98UA4bNkzOnTtXJiUlSYvFUlNWQkKCrKqqkkePHpWXXXaZHDp0qBw7dqw8cOCAlFLKlJQUOWrUKDlgwAD51FNPyYCAAKfXU3f7sWPHZFhYmLRarfL48eNy7NixcsiQIXLIkCFyw4YNUkopR44cKYODg+WgQYPkm2++6fK4ujh7dsBW6QUddffPmS6rOtix6qDCMzSky21GeRWe4WwUyWAwSGyZdWr9GQyGJp//H//4h3zooYfqbR88eLDctWuX/PDDD2VMTIzMzc2V5eXlsn///nLLli3yyy+/lH/84x9rji8sLJRVVVVy9OjRMjs7W0op5eLFi+Udd9whpbS9BB1f2lOmTJE///xzzXF33XWXlFLKiRMnysOHD0sppdy0aZO88MILpZRSXnXVVXLBggVSSinffvttt1+CUkoZEhIiT58+LcvKymRFRYWUUsrDhw/L6nr9yy+/yCuvvLLmeFfH1eVcNMSqDrb9OqjwDA3pshqaPodJSUnhpptuwt/fHwB/f39uvvlmjh8/3qLnveSSSwgPD8fPz49p06axfv16Bg4cyKpVq3j88cdZt24dISEhHDp0iL1793LJJZcwePBgXnrpJdLS0mrKufHGG2t9XrJkCQCLFy/mxhtvpLS0lI0bN3L99dczePBg7rnnHjIzMwHYsGEDM2bMAODWW29t0nWYTCbuvvtuBg4cyPXXX8/+/fubddy5iKqDqg4q1Dric5rY2FiCg4MxGo0YDAaMRiPBwcHExMQ0ucx+/frx5Zdf1tpWXFzMyZMn6dGjB9u3b6/nnSmEoFevXmzfvp3vv/+ep59+mosuuohrrrmG/v3789tvvzk9V0BAQM3nKVOm8Oc//5n8/Hy2bdvGxIkTKSsrIzQ0lJ07dzr9fVO8RFNSUtBqtURFRfH8888THR3Nrl27sFqtGAwGp79566233DruXETVQVUHFcpZ65wnKyuLWbNmsWnTJmbNmtVsZ5mLLrqI8vJyPv74YwAsFguPPPIIt99+e02vZ9WqVeTn51NRUcGyZcsYM2YMGRkZ+Pv7c8stt/DYY4+xfft2evfuTU5OTs1L0GQysW/fPqfnDQwM5LzzzuPBBx9k8uTJaLVagoOD6dq1K1988QVgm4bZtWsXAGPGjGHx4sUALFy40K1ry8nJYdasWdx///0IISgqKiI2NhaNRsMnn3yCxWIBICgoiJKSkprfuTpOYUPVQVUHz3lcjVm35J+aI2452sIcz8mTJ+XkyZNljx49ZLdu3eT9998vjUajlFLKDz/8UE6dOlVOmDChlqPMDz/8IAcOHCgHDRokhw8fLrds2SKllHLHjh3yggsukMnJybJfv35y/vz5Ukrb/Fz1MdV88cUXEpBr1qyp2ZaSkiIvu+wymZycLPv27Suff/75mu3uOMpoNBo5aNAg2a9fP5mcnCxfe+21Goecw4cPy4EDB8rk5GQ5Z86cmjKqqqrkhRdeKJOTk+Wbb77p8ri6nItzxC2FqoOeq4MKz9CQLns8H7E7qBymLceBAwfo27evt8VQNAFnz87b+Ygbw5kuqzrYflHPruVo0XzECoVCoVAomo4yxAqFQqFQeBFliBUKhUKh8CLKECsUCoVC4UWUIVYoFAqFwos02xALIRKFEL8IIfYLIfYJIR70hGAKhaL1UfqsULQ+nugRm4FHpJT9gFHAfUKIfh4oV9FOEUJwyy231Hw3m81ERkYyefJkL0qlcJMOoc+qDiraE802xFLKTCnldvvnEuAAEN/cchXtl4CAAPbu3UtFRQVgi2IUH6+qRHugo+izqoOK9oRH54iFEF2AIcBmT5araH9MmjSJ5cuXA/DZZ5/VBLcHKCsr484772TEiBEMGTKEb775BrDlUL3gggsYOnQoQ4cOZePGjQCsWbOGCRMmcN1119GnTx9uvvlmvBGI5lyjveuzqoOK9oLHkj4IIQKBr4CHpJTFTvbPBGYCJCUleeq0igZ46CFwEWu+yQweDH//e+PHTZ8+nRdeeIHJkyeze/du7rzzTtatWwfAX/7yFyZOnMh///tfCgsLGTFiBBdffDFRUVGsWrUKg8HAkSNHmDFjBtVRm3bs2MG+ffuIi4tjzJgxbNiwgbFjx3r24hQ1NKTPZ6PLqg4qFI3jEUMshNBjU9qFUsqlzo6RUs4H5oMtLJ4nzqtouyQnJ5Oamspnn33GpEmTau1buXIl3377La+//joARqORkydPEhcXx/3338/OnTvRarUcPny45jcjRowgISEBgMGDB5Oamqpegi1EY/rcXnRZ1UFFe6HZhljY8nj9BzggpXyz+SIpPIU7vYaWZMqUKTz66KOsWbOGvLy8mu1SSr766it69+5d6/i5c+e6TNXm6+tb81mr1WI2m1v+As5BPK3Pqg4qFI3jiTniMcCtwEQhxE7736TGfqTo+Nx5550899xzDBw4sNb2yy67jH/96181c2w7duwAVKq2NkKH0mdVBxXtAU94Ta+XUgopZbKUcrD973tPCKdo3yQkJPCnP/2p3vZnnnkGk8lEcnIy/fv355lnngHg3nvvZcGCBQwaNIiDBw/WSrquaB06mj6rOqhoD6g0iB0Mlcas/aLSICq8jXp2LYdKg6hQKBQKRRtFGWKFQqFQKLyIMsQKhUKhUHgRZYgVCoVCofAiyhArFAqFQuFFlCFWKBQKhcKLKEOs8DinT59m+vTpdO/enWHDhjFp0qRaoQLdZd26dfTv35/BgweTnp7Odddd5/S4CRMmoJbDKeqi6qGiveCxpA+Ktslbq87+xdMQ/3dJrwb3Sym55ppruO2221i8eDEAu3btIisri169Gv5tXRYuXMiTTz5Zk1f2yy+/bJrQCq/S2nUQVD1UtC9Uj1jhUX755Rf0ej2zZs2q2TZo0CDGjh3LY489xoABAxg4cCBLliwBXKeX+/e//83nn3/OM888w80330xqaioDBgwAoKKigunTp9O3b1+uueaampyzYAvmP3r0aIYOHcr1119PaWkpAF26dOG5555j6NChDBw4kIMHDwJQWlrKHXfcwcCBA0lOTuarr75qsBxF+0DVQ0V7QhlihUfZu3cvw4YNq7d96dKl7Ny5k127drF69Woee+wxMjMzAVuc37///e/s37+flJQUNmzYwB//+EemTJnCa6+9xsKFC2uV9e677+Lv78+BAwd4/vnn2bZtGwC5ubm89NJLrF69mu3btzN8+HDefPNM3oKIiAi2b9/O7Nmza7LuvPjii4SEhLBnzx52797NxIkTGy1H0fZR9VDRnlBD04pWYf369cyYMQOtVkt0dDTjx49ny5YtBAcHn3V6ubVr19bED05OTiY5ORmATZs2sX//fsaMGQNAVVUVo0ePrvndtGnTABg2bBhLl9qy+61evbpm6BKgU6dOfPfddw2Wo2i/qHqocIbRZCG9sILiChOllWZKjGZMFisaIdBpBFqNwN9HR6BBR5D9L9TPBx+dZ/qyyhArPEr//v3Peg7NU+nlpJRccsklfPbZZw2ep7FzNFaOou2j6qGiMYoqTBzNLiElp4yMQiPWJuRdCPDVEurvw+TkWPx9mm5O1dC0wqNMnDiRyspK5s+fX7Nt9+7dhIaGsmTJEiwWCzk5Oaxdu5YRI0Y06Rzjxo1j0aJFgG0Icvfu3QCMGjWKDRs2cPToUQDKysoa9ZK95JJLmDdvXs33goKCJpWjcI2UEovVisUqsUrbn5QSScslnFH1UOGKonITK/ed5qMNqaw9nEtaQQVWKSnOy+btR26hOD/H7bLKKi2kF1RgsjSvLitD3ERMFivFRhO5pZVkFFZwMq+ctIJyThcZySmppMRowhuZrbyNEIKvv/6a1atX0717d/r378+TTz7JTTfdRHJyMoMGDWLixIm8+uqrxMTENOkcs2fPprS0lL59+/Lss8/WzAVGRkby0UcfMWPGDJKTkxk9enSNM4wrnn76aQoKChgwYACDBg3il19+aVI5CtdYpMRosmI0Waiosv2VV1kor7RQVmmmrMpMeZWZCpMFo8lCldmC2Wq191CapkOqHirqUlZpthngjansyyiu1wNeufAdju/dyspP57kooeVQaRAbQUrJvmMnuP2Wm5nzynuIgFCK7PMI1beuOC+bj//6MH946i2CwyJrfqvVCAJ9dQT76YkI9CEmxEBssB8h/voWk1elMWu/dNQ0iEaThaIKU5PKF4CPToOvTouvToNGI5ojrqIROur742h2KT8dyKK8ylJv35zJyZirKutt1/n48up3u90q/86xXQnxa/i93pAuqzliJxhNFo5ml3Isp5SMQiOfvvEM23//jX+8/jeu+9Pcesc7tqQc91uskqIKE0UVJk7ll9dsD/DV0iU8gG6RgXQO90evVQMTCoUzJFBptlJptgLgq9Pg76PFR6f1rmCKdkGV2cqvh3PYm17k8pinF6zm7UduITfjBAB6XwMDx1zClJmPt5aYyhBXY7ZYOZRVwqHTJZzKt80Z1G0pbfzuMzZ+91lNS6mx/a4oq7SwL6OYfRnF6LWCLhEBDIwPISnMHyFUi1+hcEW1UdZpzQT46PDVaZTOKJxSWF7FNzszyC+rcnmMs96wqdLIjjXLueUJ29IyVyOenuSc74oVG02sP5LLv9cfZ+W+LE7kldfMHTy9YDVDL5yM3tcA2FpKQydexdMf/+TWfncwWSRHskpZuj2dDzeksiU1H6Op/vCJQqE4g9liG20qKK/CZLF6WxxFG+N0kZElW041aIThzDu8GiE0RMR3pvewM8vWWmPu2CM9YiHEf4HJQLaUcoAnymxpCsur2JSSz6HTJS7d1oPDozD4B2KuqkTn44u5qhKDf2BNq6ix/WdLUYWtUfD78XyGJIYytHMnDPqzH4KTUqpeQjujrTj2eUqX3a2DFrOJvMw0wmMT0erO/nVkskjyy6rw99ES4KtDo+p9k2krdbC5HMspZcWeTLc8mV+67eJaPWIpreSmn6Aw53STRzybgqd6xB8Bl3uorBalxGhi9f4sFmw8wYHM+p5z9Y4vzOP8yTN48B+fc/7kGZQU5J7Vfmc05iZfZbay+Xg+b36zmSEjzudUeobb12cwGMjLy+swSnUuIKUkLy8Pg8HgbVHAA7rsrA5azCayTx3HUmfdbHFeDpUVZRTnZZ/VOeqWV15lIa+0CpNZ9Y6bQhurg01mb3oR3+1yzwiDrUccEhEDIgy4HMTd+Pq/wuBxxxl0wXFiu3yKRjsTuAStPvqsRzzdxWNe00KILsB37rSiveE1bbZY2XqigC3H8zFbvWekivOyefO+aRTn53D+5BlOnb+q+fKfc/lt+WLGTbmJ+e+/S6/ooEbLN5lMpKWlYTQaPSe0osUxGAwkJCSg19f2vPSG13RzdbluHTRZrOTl5VNZUYavXwD+QSEUZGfifGmSoFNUbKMylpcU1SrP4ecYdFqPRTw6l3BVB9sLe9OLWH0gC3dNWmGujn0bA1i9OIui3L7AmesOCDEjgNIix1EaM8FhRxg3LYIBo0uJSjyzEqC5XtOtZoiFEDOBmQBJSUnDTpw44ZHzukNqbhm/HMqmsLxpSyg8hbtu8q6O0/v4cjq/mLAAnxaVU9F2aIuG+Gx02c/Pz2mjUKv3YdDYS9mzcTWmSmMtT9WGpnbc1aEeUYFc2j8aX+VdfU6wP6OYlftPu2WEC7J1rFwYzpYfg7FaBT6GU0Qm7GDMVeGk7F5CRflR7nr+7/a6JoFYoAtwEXAVMBiAXkPLmHhDAT2HlHPXBe3EEDvS0j3izMxMpk+fzkefLOJAkYZDp0s8fo6yYg2lhToqSjVUlGqorNCg1Ul0etufX6CVTtEm/IOsCOH6BSI0Gp5btLbWy6c4L5tv57/i9CUVFhHFmJ4RDEkMVfPA5wBt0RA70pguZ2ZmMvP+B/lx+f9q6vKA8y/joul/ZtXCL9i1dgcaXTBWs5F+I8cy/tpbCQ4zEx5jRudT/93UkG7U1aHPXnmE5cu+onvnhCbdB0X74ODpYn7Y27gRLinQsvqzMDYut42gnH9lEedPLiQ6yXkHzVVdG3bRfXz1zyNUVc6ktNCHhB5GXnhWxx9u0qFtoN13zq0jfvHFF1m3fj23PzCHq+9/zuVx7rilSwlZJ3w4ttuPtKMGsk/5kH3Kh7Ji91rael8LcILuA49SVrSC9GMLkHIrYDPKwyZOqXfuhpzAzFbJr4dySMkp49L+0QQb2ucwkuLcIDg4loK8CZgq+yPEQEyV/dixpic7ftECttCSVvu08f7Ntj8AISQhEWaiEqroPqiCPsPLiO9R6baD5MqF73B41xbuevAJvvr0P4QH+qLoeBzNLuXHvc6Hox3f76n7u/L5W9EYyzScd2kxl96SR6eohmOJu6pr+35bQEH2YkZe/jtJfV7lly/CeOl5LX+4qenX0aEMcd1hsLXfLmLtt4tcerm5CsRhLNOw97cA9v0WyLHdfjXzBP7BVVgte+k/Ko747hoObFnCkR0rGDhmOJf/4U4sFoHZZPsrL9ZSkKXj95XbyDxexfF9Y6gy3g3cjc0I/0Rw2EZKi444vZZqJ7ABYy7hs1fnkJ+VXmv/qfxyPt10gkv6RtPTjbljhcIblJXBhl/vBayERpYjNAfQ6j5n7JSLCOpkJrCTBV+DFatFYLWC2SQoytWRl6knL9OHzOM+rPgoghUfReAXaEKr/4Hw2CDOnzyDUZNuZNP3S2o5PdYdefr1m4VEBC3E19eA0VjhREJFeyW9sIIVezJdOtyuXPgOKXv28+6cSrJOxpHQ08jNj58mOunMkqbGOmPV7+FRk27kzfuuYeN3ZxJwbFrxCZtWfIJWb+CnPcVotU3vFHlkaFoI8RkwAYgAsoDnpJT/cXV8Sw1N7ziYwuwH/o/t61Y2OGzlfJhYh0Z3Hf1HzefA5gDMJg3B4WZ6DSmne3I5PQZV8PPnT7Hp+8UgBNJa3zvT0eA7P0csiNEMvfBV9m8KxFgejRCSLv0rGD2piEHjStHXGY6rdtgafeV0l45dwzp3YmyPCBX+rwPS2kPTLaHLy1aXcagiAx9D0941JQVaDm/3Z+WnR8lJHwwE0blvBWOnFjJ4XAlah+6Eq+HEa2c/wc0XDqJzeECTZFC0LXJLK/lia5rTmAtn3r0DgC+BnsCraPV/4bXl22od6877tZqGpkUemjLC+yEupZQzPFFOM87P5uP5bE6zoPX1b3TY6ukFqx1uaBAa7b1otQ9gqgondb+ZUZOKGDKhhM59jWg0ToxqncaLs5Botc9R/dCGM2XmwwSHWZCyiNOpFezZEMiWVf4sejWWr981MfKyEsZOKeTlP/Zzew3bthMFnC42cuXAWAJ8O9Qgh6KVaQld7jfQyvHdTW/wv3hrfwddCAJu58SBBzhxoCerPg1n6uxs+p5nCyHrbDjxyM5NmCyS/+3K4KpBcU6NcbVfyZIlS5qcBELROhQbTSzbke4y8NHTC1bz6cs/cnTXY0ApWv0VDLrAlykzf6g5pilrhD0dN8KRdu/jX2I08eW2NH47lodVSrfW9QaHRwFdMFW+DZzCaplLQEgaf3wxnecWpjDtvhy69rcZYXAeQSsirjOAywfS2EMTAmK7VnHpLfn0GvpH4GL8A7ez9utO/PWOriSPTaX/6DvcjtqVXlDBos0nyS5Wy5YUHYva+leC3vcDhlz4KNMfPYCU8MFTCfzn2Thy0m09Esd3QFRSd0ryc1j56bwaY3wir6zeOV588UXWr1/PCy+80MpXpzgbjCYL3+xIp8R4Zn63blyGg9t6cGz3c0AqWv0FWM2r6r2fmxoVsSlxI9yhXXefjuWUsmp/FhUOGTXueO7tms/XPlDfUasoT8vqReHsWPMmQgPJY9PQaOZhNu2m30jbb+vOGzgzqlarhTFX3eR0nqoax/kFZ8fUbZXlZf4ExCE0T7Fr3Sys5vlIOQKt/q+Yq9IabX2VVpr5YlsaVwyIoVtkoNv3UaFoyzjTP7+AAEZcqmXohBOsXRbKqoXhvDqzM5PvyuX2Z9/m8auS2fC/RTVlOPZ4WLGHKYPiSQr3r+dX8u677/Luu+9iMBioqFBzym0Ji1Xy3e5Mcktrh62s9vX58ZN5BIf/kx8/jiAwZDv9R7/D2KlvOX33NrV325h9aSrt0hBbrJK1R3LYebLQ7d9UVgh+WhzGmq86YbUIRk0q4uIZ+XbPuZm1jnXmxOXMqFY/CFcPpLGH5nz4ehhTZg7CYk5l3qP7KMj6I1rtXST1XkRx/mf1yqhLldnK/3ZlMq5XBEOSOrl5dxSKto2rRq3ORzLxhgKGX1TMF/+M5pv3oji605+H317D6s/+4nQ+z2SRfLsrnWlDE0hJSeHRRx9l2bJllJeX4+/vzzXXXMPrr7/u5StW1GXV/qxaWezqdmR+Wz4YiEBoPuHZRSPR6ecArt/PDXWUtBqBlDQaedFTtDtDnFdayYq9p8kpqb8m1xlSwq51gXz7fiSFOXqGXljMFbfnER5bf+2Yu/MGDbWEziZTR8OtMjNPf9yb7FOnWDovisPbbyOu240c35dN1/4NDz9bpWTNoRyKjWbG9YxQ640V7Z7GGrXB4RbunJvBumWh/O+DSNKPDiOh5zDMVcud9nhMFsk3OzO4blgCwcHBGI1GDAYDRqOR4OBgNU/cxvjtWB4HMotrbavuyOzesBpz1VxgDhHxy7nv9W7o3HBgrq5TIX56pl40mi7hAQT46vD3seW+tlgluaVV5JRUkl1i5FBWCZWmlgmh2m7miKWU7DhZwKLNJ8kpqWw0XjNATrqe956I5+OX4vAPsnD/G6e45cnTTo0weCab0tlm6mhsziEq0cQ9f0vntqczKC/W8vbDiXw7PwJzVePGdfuJAlbuz8LqxZCeCkVzcUfXweZ3Me6aQv7095PofCR7Nz1Nj0HzXOqW0WRh2Y500jIymTVrFps2bWLWrFmcPn26JS9HcZbsyyhiU0peve3B4VH4+gXajfDjwLv0HLyQkPDGnae0GsHgpFBuHpnEnWO7ckHPSBLD/AkL8MGg1yKEQKfVEBNiYGBCCBf1jeausV0Z0yMCfx/PR2trFz3iEqOJ1QeySM09Myzhag0wgNUC674J5fsPI9DqJNPuz2L0lUUNRj2B5nnFNTVThztzDkLAoHGl9DmvjP99EMmaL8PYv9kXve/d3P3SfQ3Ktz+jGKPJwpUDY9Fp2027S6GooSFdd0Zir0r+b95JPno+jiM7Z7N/cy7T7u+Ds4Gh0koz1zz2Fjeel0iAr45581ou1Z3i7DmZV85PB1wnBDm6+yrgJgZdkEpAyIZGnaeEgD4xQYzuHtHocqO6+Oq0jOgaxpCkULadKGBzSr7Hhq49FuLybHB3HbGUkp2nCtl4LI8qe1aVxmLNZqfpWfx6DKn7/eg3spTrHswmNKLhCCqOfPj8/QSHRdaaN6g2loG+OsICfAgP9CE8wBd/Xy0GvRaDTkNBbjZ/fmIO3337DRUVFRj8/Lj4isnc+/jz6APDyCmtJKekEosHeqf7fw9gwYsBmCr96NLvcx54a6TTl4wjCZ38mDI4TsXebWd4I8Tl2eCOLh/OKmH57kyn+/RaQYifHh+dxvan1VJlsVBWaeGO8b3diivtCrMJlrwZw7afghl2USb5WdO47enXnTZcI4N8uX54gtKPNkROSSWfbz1V8+6vy5ovO/Ht/EhGXVHIdQ9m16xycUVMiIGL+kYRFeSZDFOn8sv5fk8m5VWWZid9aLM94uxiI6sPZJNVZzmOcwenS5hw3V389favKMydg94XbpqTybCLSho1UHVx7KHe+sgLJIX5k2j/a+hGhwcmEt4plMrKSgwGA1WVlSRGR3DFiL41x5gtVrJKKskorCAlp5TMIqPbmUKqOdMQCQc+IHX/rTxy2edo9ffy2vKNLn+XVlDBsh3pTB0c36QcxwqFJ/DRaegaEUBcqB+xIQYiA31dBqK5KPU4Dz/yCMuWLcNYUYGPr4GBYy/hqrsfd3p8XXR6mHzXbo7tPsK2n+4A5vDDx+9xw0PP1DrONvT9MMbX3uOOS4agVYFxvE6J0cQ3O9NdGuEtK4P5dn4kg8aVcN2fGjbCQsB5XcIY3S3co0GPEsP8mTEyyWUj82xos4b4x32n67mpV3Nk12bb8LHeB1OlEYulEx8+H09B1lRCI3fz4D8DCAl3vti7MQx6LT2iAukTE0RCJ7+zcnTKyspi1qxZzJw5k/nz55OZWfsB6bQa4kP9iA/147wuYZRVmjmWU8rB0yWkF7i3VKJ2Q2QaGu2TWC0vohHDOLaniO4DXS9byig08vWOdK4ZooyxonWJCvZlYHwIvWOC3O51xsbGEhoSQlV147aqiiHd45g0si970oooqmg8m9oLt4y3R8HbDvyLTd/Dpu8HoPPR1PSqq4e+P3z7dRLi3+Cy/jHKwdGLVJotLNuZUWutsCP7fgtgyZvR9BpSxs1zTqNpoDoFGXRc1j+GxDD/FpE12KDnhuGJzc7/3maHpj/5LdWpIf7yn3PZ+N1nRHfuQUznHuxaawI+AcKwTdj/E5BuD19VEx/qx5CkULpFBnqlRZxbWsnutEIOZJa4bAVW8+U/nuO375eg1ftgMVURGnkzBdlvoNEGccez+fQfXT9ggSNRwb5MG5KAXws4HSg8S0cYmrZYZZN1atq0acTGxtZq3C5duhQpJSfyyvktJY/TRfVXETifwroXmEdg6DrCov+P9JQ9WEz13zE+vr5UqnzeXsFssfL1jnTSHDomjitRcjMSee+JBGK7VDH71VMY/F3br7hQA1MGxbeZ91yrpEE8G5piiOsrlgZ4DngaOAjMAHa7ndcUbJ5zvaKDGJoUSlSwZ+YNmkul2cLOk4VsO1ng0lW+eh5704ovsJirewXxwNfAMDTaR3l9xSzA9XKqiCBfrhuqjHFbpyMY4pZESsmhrBI2HM2j2KGHXDcu8BlsxhiWMXTiQpBGp2uNrz6/P4MSQ1v5as5trFbJ8j2ZHM0urbW9Oh70oHEPc2j7KwSGWHjgzVMEhroe9ewZHcjl/WPalINqu5wjdqQ4L5v4bn0IDovk4Lb1mCo7AYuwxab/L0I8hJQlaPU+bnk6a4Sgb2wQI7uGE+LfttII+uq0jOwWzqDEULafKGDHqcJ6PeTqeexLbprt8LJJR+dzOQHByynKfZOl8wq4elaOS4/T3JJKvtqexnXDEtQwtaLdIoSgT0wwPSID2ergyeq4AqI279j/z2P7z+XALVSPoDm+O345lE2Ar44eUSpCXWvx08FsjmaX1nQeThzc5TBiEcPOX+8Hcqkyjicw9FuX5Qzt3KndxU9oF4Z45cJ3OHloN1GJ3TFVjgUWAgHAbeh8lmCuqiSmcw9ufuINl+EmwTZp3zs6iFHdwukU4NOal3DWGPRazu8RweCkUNYdyWV/RnG9Y+out7KYCug38l/4GHrx61dhrP9mHbAMkE6XU+XYjfG1Q5UxVrRvdFoNo7qF0yU8gBV7MyksN9Ws0R845hK+/Odc8jJPIaUVmzEOBF7Bx1BB8tjljL/29lrvDinhh72ZTBuaQFyonxev7Nxg/ZFc9qYXAWfm7IddNBWrxczuDZswVy0HIugz/BmmP+oyGRjjekUyrHP7iyjYpg1x7eFoQdbJ6cDzwCFCIqaj1R3hjuc+r1Gg+O59XK7FjQkxMKF3JLEh7Uup/H1szgb944L55WB2vXlzZ2Ha7ngul4DgAr7/cDJCrEbKSeh9K+tliALILlbGWNFxiAkxcPPIzvxyKLvWCoi809VGuJpXgUiqjI8SHjuV+O759d4dtlCYGdw4PLHNN9zbM78dy2NLan696cetq5dhM1H/A5KBKYTFhLoc7ZzQO7LdhvVt03PEKSfS7CHMdmCu+jdwOZ2i1zDrFV8i48LdOleAr5axPSLpGxvUroYqnGG1Srak5rP5eL5b65Hfe+J7Dm+/F8RBkJdx/uQLXQZEiA42MG2o8qZua6g54qazN72Inw5kY5WS+U/dTW76CfIyT9o8XIUgPDYJU+W7FOddwbT7sxg7pchpOcF+em48L5FAlWLU4/x2LK8malbdeX2djwGN+JCqyulcestOyopeqxXXoRohYELvKAa38Tn9hnS57cxkOyE4PApT5RDMVRuBC4HZ9Bk+3y0jLAQMjA/hD6O70C8uuN0bYQCNRjCyWzg3npdImBstdF+/7+k38k10+n4Y/HeQl+n6RZJVbGTZjnQqzU1b9qVQtDUGxIdw9ZA4fHQaZv7lA3oNGQ3YAoIIoPfQ83nm0+70H1XK1/Oi2L3e+XxwcYWJrxvIf6toGo5GGJxFNvwzVZXTuezWXC7/gz/XPvCcUyN8YTswwo3RZg2xlPDr0lD2/PYsvn4Gbn5iE2OuKqa0sPH8j5389Vw7NIGL+0V3yB5edLCBm0YmMSgxpMHj7njubf744jTufS0DoY0kM/W/ZKe5dk7LLDLyzY6MRpdPKZpOpdnC9pMFbDzmmTymiobpHB7ADcMTCTLonMZ112rh1j9nktTHyKcvx5C63/nqidySSr7dmYHJonTDE2w8lus0fnT1Mzr/yp+BZwiJ+I5Lb8l3Wc64XpEdwru9TQ5NFxXBxVeXsXVNAANGlzL90dP4BzWuAELA0KROnN89vE25rbckh7NKWLU/q1HjmXHch/fmJCC0MPuVNGI6Ow+WAhDfyY9rhsSjP0fuYWtQVGFi56lC9qYXUWW20jc2iMsHxDb6OzU07RlKjCa+2pZGQbnzICClhVr+8VAixlItf/rHSSLjnR/XLTKAq5LjPBqh6VzCapX8fDCbPenOpwEA9m0K4D/PxgIrGDbxvxTknHKazW5Ut3BGd3dvirIt0OJD00KIy4UQh4QQR4UQTzS3vPvug+3r/Lnq7hzumJuBf5C10QwsQQYd1w5NYFyvyHPGCAP0ig7ixvMS6dTIMqy4rlXc+3oaAO88mkDGcddD2+n2cJiq9d98Mgor+G53Bh9tSGX7iYJ2MdrgaX1uCwQZ9Fw33LWeBIZamPmXdBCSD56Kp7TQ+UhaSk4ZK/efbnYkpXMRk8XK/3ZnNGiEH5s0i/882wnYCtzAtp+XkrJnC8/fNK7WcYOTQtuVEW6MZlssIYQW2wr5K4B+wAwhRL/mlPnyy/Dk26e58PqCmljRDaUX7B0TxC2jOrdYGLO2TkSgLzNGJtEtMqDB42I6V3Hf66fQ6iXvPpZI+jFfl8emKWPcZKSUHMkqYfHvJ/ngh208eMvVFOa5ziDTlmgJfW4rBPrqGjTGkfEm7nohg8JcHf95Lo6qSue93gOZtlEoZYzdp6LKwlfb0th58LjLDlVGig86n5X4+uUDVwJnsu1Jq5WHL+3NnMnJ9IsLZkKvxjPitSc80XUcARyVUqZIKauAxcDU5hSYkAC9B9nc2OdMTubhS3uz8bvPkNK2Hrb6geg0gov7RjNpYGyHnAs+G3x1WqYMimt0DV1Ugon7Xk/Dx2Dl3TkJpB9z3TNOK6jgGzUv5jZmi5XdaYV8tDGV73ZnkllkPOv81G0Aj+tzW6LaGLtyduzS18gtT5zm5EEDi16Jweqi6u/LKOanA9nKGLtBZlEFCzefaFAfcjP0vP9kAn7+kgGjXwZyEeKMearODf/B8t+4pG90h3C+dcQThjgeOOXwPc2+rRZCiJlCiK1CiK05OQ0n+K6mOC+buG59GDDmEvS+NieK6gfy6udrufG8RAYmNOywdC4hhGBcr0gm9olC00BFjYgzce9rNmP83uOJDRrjU/nlfLR6BxeMG6cSprugymxla2o+/91wnJ8OZFNYbmqwAdnGaVSfm6LLbYlAXx3ThsYTZHC+iiB5bClTZuawe30Q3/07wmU5e9KL+OVQ+xjp8BY7ThbwxdY0Zl/cz6U+pB8t4NW7fTCbJfe8nE5V5WHGXHUTA8deCoAQGsxVlYSFhnDzhYM65Px8q02mSinnSymHSymHR0a6N6ywcuE7nDq0m5xTxx1c2iuJCg9l9qRhbSY+dFtjUGIoUwbblm24IjzWZoz1PjZjnJHi2hh/9PYbbFi/nufmzm0BadsvlWYLvx+3GeB1R3IpqzyzvOXpBasZeuHkeg3Ipz/+yVvieoym6HJbI8igZ5qTWOvVviiDJxxh7NQC1nwZxvpvXTf2d50qYtX+LKweyDPekTCaLCzfncmaQzlYrLJGH4Q9VZLQaBk68Soe/Oda3pkTj9kUQo9BzxOdVMUdz73NtQ88h7RaGHPVTTz8ztdcePXN+FtKO6wDqSdWqKcDiQ7fE+zbmoyfnx9Gh+wnWSePArZ5gqtvuh3KC8/5oejG6BoRwHXDEvh6RzoVVc7XP0bEmZj9WhrvPJbAu48ncO+racR2dZ1oY/777zP//fcxGAxUVLiXtrEjYrZY2ZVWyJbUApf3tv6ayMZjoLcRPK7PbZWwAB+uHhzPV9vTapzoqodOVy2cx7T75lKQpefrd6Lw8c3j95W3O/Xe3ZtehNFk4YoBbSvJQEuSmZnJ9OnTWbJkCTExMbX2Hc4qYc2h7FoN05duu7jWu0RaLWz/eTPbfw4EOgFXsGf9Oh6+9C81YXir1wwH+up45rYrGswH397xRK3ZAvQUQnQVQvgA0wHXEbndICUlhdGXTKnVmxh20VX8sm0fSz/9L0uXLm2+1OcA0cGGmjWUroiMt/WMtTp4d04Cp0+c6Rk769WNvmwq+w4daXHZ2yJWq2RPWhEfbUxl7eFcl0a4GmfrVtsBHtfntkxMiIHJybFOpxIevaI3B7clkdCzks/f6k7KHp3Luf6j2aUs25lxzgTEefHFF1m/fj0vvPBCzbYSo4lvdqazfHdmLSMMtndJSEQMGm11jzgWjWYtQtMVrX4qsM7pqJGPTsPUIXEd2giDBwyxlNIM3A/8CBwAPpdS7mtOmbGxsRgCgmr1JgZ0juGC5J7NFfecIyzAhxsaicQVGW/i3ldPIbQw79E43rzvKYrzc5z26oTen7WnTG4lZe9IHM8t49PNJ1h9IMtlwvK6VA+xVcdArxsVqC3SEvrc1ukcHsDqTbudTiUIyjh1KBGr9STwPzZ+t93lXP+p/HK+3JZGsbHj6oafnx9CCN59912sVivvvvsuQgh8DQY+/u0EKTnOc6EHh0fRf+QEpNWKVt8ZaV0NojP9RryG1fyz01EjnUYwZVAcUUEdfwrSI+MoUsrvpZS9pJTdpZR/8USZxQW5nD95Bs+8v5Q7755JcfvoTbRJgg16bhieSFSw6+VKUYk2Y1xlLCftyD/45r2vAOe9uoJyE59vOUV2ScdPnp5TUslX29JYtiOdPHvCjeK8bP7xpxv4+59udLmu3RV5mTo2rQjmlTmhbNzYEhI3n5bQ57bO+CG9SIqJqDeV8PTHPzH0wlHofKYAZuBHBoy5zeVcf3ZxJYs2n+RkXrnT/e2dlJQUbrrpJvz9bUtFfQ1+nHfxFP684KdG18iXFOYxeMLDBATvQ6PtRue+T6HRrnc6aiQEXDYg5pxZktpmo5g/+Lf30Gs1TB0cj9/tk7wtTrvHz0fLtUNtc8ani+ob0DlXDsRsqsK2dPQXdqx5mB1rJiBECs99to7gsMha2WlKK818sTWNq5LjSArveMpSUWVh47Fc9qYXY62zRGXlwnc4cXCX7XOdPM91sVggZbcfu9YFcXCrP/mnbSMTncItpKW1mPiKJiCMRUy64Q/0mTCtJpNZ9aiQxXQIre5qLOZVpOyZi97HDDg3PBVVFr7ekc6YHuEM7xLWuhfRwsTGxuLjF0BFhRGdjy9VlUb0hgC3fB8m3vBv/v1MPEJIHvxHOom97qnJPRzUKbzW+2V8r0h6RQe15KW0KdqsIe4aEcjIbmEd1kvOGxj0WqYNjWfZjnQyCmsb48HjJ7F19TI0mkNYrROBX9Bo1mG1XuDS2FSZrSzbmc7FfaPpFxfcOhfRwlitkl1phWxKya8X5L+u8xrgNM9zUW42Hzz9EbFdX+LQ1ghKi3T4+FrpObSc8dMK6Tm4nAmjfLliYOMhLhWtx9KlS7FYJUu3pxHfvU/NdsdUoys+epsDvz/C/KfM3PO3NAz+zr2lrVKy7kgu6YUVTOwTRZChfc9xGk0WjuWUciCzhG0HUzl/8vRaqVcbY//mAD5+KZaA0EoCQ24mJOIRILLWuuLqd8yIrmHtNp1hU2mTsaYVLUu1AU0vqHBqXGwMAH4GKoEJwLFaxqYuQ5JCGdczsl2v8TuZV86vh+vnfK6mOC+bL//1PPs2/Yy0R3oQGg0DRl/Epbfcx5f//Bd9z5vPr0t9qChNQKM1kjy2ikHjSul7Xhk+hjO6pmJNt12MJgtLtpwiv8x5Pdi9PpCPX4ql64AKbnx4J4vfeMipN3U1vnoNY3tEMDA+pF0FoiivMpOaW86R7BJO5JW7lXq1LlYr/LwkjBULwonvXklMl4fYtno+CFGjQ474+PpSaeyYU17tNg2iomXw0Wm4Zkg8CZ386nlGC42WHoNHERZdAOJiwBf4lX6jZja4BnbHycJ2myquqMLE/3Zl8NX2NJdGGGwOJ8GdImq9QKTVit6nH5++7M+JA//jh4+7UVGaDtyJ1RLOzl9DWPhKt1pGWNG2Mei1TB0cV2+NcTXJY0u56fHTpOzx453HAkjZs7vByGmVJis/Hcjmy21pLeZXkZmZyfjx45sVdMdksXIqv5wNR3NZuPkE89em8OO+06TklDXJCFeUavjw+Ti+/zAC+IK0IxFsXfU+Usp6Rljva2DCpGmkHj/eZPnbM8oQn6PotRquHhJPvx6da3lGI61EJXSlz/ALEOxBq7sC8OHIjpepKK0XMK0WJ/PLWbT5JJlF7WONcZXZysajuXy8MZWj2aVu/aakMA9fP1tM75CIy4DP2P7LG2SdGAcsBAYDo4AP0ftaO0wQj3ONUH8frhoUh87FCM/iN7oi5V0U5gwCvmDjd0sbjZyWVlDBos0n+XZXBtnFnjXIzpYTNUZZpZljOaWsO5LDki0neXfNMb7clsbvx/PJLq6kOYOlaUd9efP+JA78HsDVs7N59tM4hl54YS2v9Ii4zgghapzj+iRFERt7bk7XKEN8DlPtDGcuK6znuVg9L/bQv55lyITXkFbJvEcajk0Ntt7l51vS2JSS12ajDUkp2ZtexIKNqWw+no/ZDTmL87J5+LI+7NmwisqKJOALinJ/ACYBb6Dz6QvcjUazFwCt3qfecozGMogp2hbxoX5c3C/a6T7bSFIOGt0D2BIUrKDfyKmNNrqkhGPZpSzcfJJlO9I5klWCuRmx3F0tJ/Lz86t1XKXZwqn8cradyOf7PZn8Z/1x5q9N4dudGWxNLSCj0HhWvV5XddlcJVixIJy/P5CEqVLDfa+fYtw1hYRE1F8KabVaOH/yDF7+6Ftm3jOLrKysJt+H9k6bddZStA4+Og3rVi2vmTN29Fys5tY/9yE7rYj35gTxzpxEZv0tjcReZ+aVqz0fq+fJrFLy27E8TuSVcXn/WEIaSdHYmqTklLLxWB45Jc7mxV2zcuE7IDvj6/c6lRVXA+VotH9lwPl78PGpYNvPqTUvmJjOPbj5iTfqObI4OqaM/Nsbnr0wRYvQNzaYwnJTvST21d7U0jIPRAHIjzi8IxyN1h9X3tR1OZ5bxvHcMnx0GrpHBtAjKojYEAMBvu6/llNSUnj00UdZtmwZ5eXl+Pv7c+VVU/i/p19iU0oeuaWV5JZUUlhhalYPty7OnKxOHDCw+M1osk74MvziYqbOyiYg+My9cHR6q9aN/3vuFaYMjkM/41LPCdcOUc5aCsDuwLUjnfRC18PKeZk63p2TSHmJhjtfyKBHsu3YL/85l9+WL2b0ldPreVfrtYJhncMY3qWTVz3g0wrK2Xg0r8Hrc4bNmc0XeBJ4CJDYsgS+ghB5jL5yOiUFuQSHRdZ6wTgG73DlEOdOqFDlrNU2+HHfafZnFNfa9shlfZGy2tBcCXwBHEern8Jry79v8rkCfXVEBfsSHuCLn48WP70Wfx8tGiGwSIlVSqxWSXmVhbIqM3/78yMs/+ITdHofzKYqp3roKZzX5XiE5kXgdkLCzFz3UDb9RjgP7OFIUpi/zQifIytjGtJlZYgVNVSZrXy9I63e0iZHCrJ1vP9kPPmn9VjlTVjNS+od48y7OsigY5wX1gYezy1j24kCTuWffYAFiwXWfKHhx08iMJtCgQXAU0TG67n2gbns2bCyntF1RnFeNt/Of4U9G1djqjSi9zUw9uJJLPr3vHpxeuuiDHHbwGKVLNuRzkmHelT3uWr1FyOtyzAE6Ln92ayahmpL8+Hz9zfYEPQkta/ZD432aaS8DyH0jJlSxOV/yMMvoPERgXPNCEPDuqyGphU1+OhsDlxfb08n00nQD4BOUWYeeOsU/3k2ntQDn5HYaxynTzxSY2AGjrmEKTMfr/e7EqOZ5bsz2RZSwNCkTvSMCmyxpU5mi5VDWSVsP1HQoBd0Qxzb7cfSd6LITLF5jcMjwDYActLh38/e43IpV12chQr1Dwhs1Agr2g5ajWDyoFg+35pGrn1ao+5ztZh+YtC458g49lfem5PAlHtyuODqQlp6xZKj0XU2teRJgsOjQHbDVPk34A6sliAi4zdwz9/iCYtxL/Rrz+hALu9/7iTIcAd1JxS18NVpuXpIPDEhruO7BgRbmfVyGv1GlHHq8L2YKp9Gqze4lWHodJGR7/dk8uHGVLadKKC8yj3lbQwpJafyy1m1P4v561JYuS/rrIxwtfPJwa3Heeb648x7NBFjqYbrHzxAcPh0hGYnYFveFRIRc9ae0HVDhRbkKYet9oavTsvVg+NqJVGp+1wt5t089K+T9B1ZxrJ3o1j0agzGsqa/ZtuKg5/FDAd+9+c/z8axfc1bCHE/fYaXMHj8n4np8qLbRnhQYghXDoxVRrgOamha4RSjyeIyHGY1Fgv89bZjFGRfQc/BGYTHPk1ZcdpZDYsJAbEhBrpHBtI1IoCwAB+3gx4UG02cyi8nraCCU/nlbidjcMbnb73EphUJIJ4GqSGh53LufyMZH4Pky388x2/fL0Gr98FylnNwPjoNMcEGYkMMhPjr8dNr8fPREmzQu+WUo4am2x55pZV8sS2twexbViusWhTGyk/DCQ4zc90D2fQf3fi8aV0a8r9oaawWSNnrx441QexeF0RZsZbAEDOjJxdx/uRCQsLPLmbA6O7hjOoW3kLStn3UHLGiSbhjjKWEtUtD+faDSGI6V3Hn3HTCY5tuEHUaQWiAD+EBPoT46RECNEKg1QiqzFaKK0wUG00UV5gprWz6eao9vU8e3IXZdBHwD6AnsBTbMHSqTR4fX/qeN+6s5uAignzpHR1Et8gAws+iYeEMZYjbJtnFRr7YltZoooMTBwwseSua06m+DB5fwjX3ZhPUqXED5srBz1V0u7orF5pKUZ6WQ1sDOLg1gEPb/ako0eLja6X/+aUMmVBCn2Hl6HzOzmbotYKL+kbTN7ZjhMFtKsoQK5qMO8YY4NBWfz7+aywaDdz650x6DXXtHOWpl0Zz+PKfc9n43XZCIj6lKHc4cAh4AFgF2F54yWMvZcrMx92S0d9Hy8D4EPrEBjeYcvJsUYa47ZJWUM6yHemYLA2/Q80m+OWLMFYuDEOrhfMnFzLh2gKCG+hROnPwq/a/cFYfHXvOl958r1v6ZbVCTpqe1P1+pOz14/heP3IzbHU3sFMV/UZU0Oe8MvqeV4avX/1rdEePQ/z0TB4Ue06kMmwM5aylaDIGvZZrhsTzzc76iSIcie2aSkTcnzCWf8p7TyQw4bp8Jt2e57T17GwNYmvx2JUDsZj8gWeBBRTlVgCPAv9EaCxIqy1+tMVU1eh8N9h6v0MSQ+kTE6Tmvc4xEjr5Mzk5jm93ZTQYDEOnh0tuymfQuBJWLQzn16WdWP9NKCMvLyb5ghR++HgWtz1d25g5c/BzVh/r9pyrk5BA7cxgVgvkZuhJP+ZL2lEDpw4ZSDvii7HcFsYzIMRM1/5GAkI+58SB1xh4fn+uf3Bug9ffmB53ifDnigGxGPTOQ4UqzqB6xAq3qDJb+WZnOmkFzpdkVLfIR1z+B7Tav7Pxu1Diuxu55cnTRCfZnKbOdrjN05iqBG8/vJZTh6cCYcAH6Hz+SkBwFT2HjKYwJ5PSgjwCO4UTnditwSHoqGBfzu8eQdeIgBaVWfWI2z4n8sr4366MRnvG1eRm6Pl5SSe2rArBYhbAQeJ7HOSGh4YQ160Srb175M6ypGpnrtyME/YtAkgA+mJLadofxEB8fM6jqtLWUNTqJHHdKknsZSSpt5HOfY28PrsPFpN7umlrzNZ3hKw+Vq8VjO4ewdCk0HaV5KKlUUPTCo9gslj5364MTjgkPXdlXDW6a/AP+BxjhYZLZuQz4foCKkqyzmq4rSk4Gy6zWODxya9jtTwFdAZWY5sHtr1gzp88w+2eeUSgD6O7h9M9MrBVXjLKELcP0gsrWLYjvdE542psehMCXANcB1wIaAEjCT0grnslUYlV+AdZ8A+y4h9kobykgB8+eZ+J19+H3jeUT/76N6yWCCAam/HtAXTFlqjFhk5fSEJPQWJvK/HdK4nvXkl0YlW9kaqzGQpf9OrjbF29DKHRIq2WWsf26ZbIpf1i6OTB6ZmOghqaVngEvVbDlEFxLN+TSUqOzQP06QWrXSjwI8AJvp4XxYoFEWz+MZip9/jj69f4cNvZUNfwOg6XXT3rebb/EsTPn3fCapmPf9BhKo1XYDH9gEarpdfQ8QSFhlFSkNvoeQJ8tZzfPYL+ccGqla+oR3yoH9cNS+DrHekNelNXc0ZvFmCqfB+dTwIJPR4krtsMctND2bcpgN9/DKnzq0RgHoterf7+L2zhNHOATITYj0a7Aov5AFrdcSzmXYy47GK3GpnuDIXXbXRLq+06TZVGAgKDmDyqn+oFN5FmGWIhxPXAXGzjICOklKpp3MHRaTVclRzHqgNZ7M8obkSBLdz2TCZHdhTy9TtRfPh8PAEhc+g3ciCX/WEAm1e4l1S8Gme93WrD+/xN4xxSq4Ww8bt4Nn7nD8QQ393I7c9mcGjrc2xa8aM9+EIVYdFxjb6kdBrB0M6dOK9LGD66jj0HrPS5eUQHG7h+WALf7sqgsNzk8rjqehwWFecQDCSduG7bue5PU4AypITKCkFFqZa/3j4dizkAW9iHSqAKqESrL2LExWPY9MNnNUvrQiLi6XveOEZNeqBerPPGcBYL2pG6jW6NVkufYWNJjIvGV5YzrHOnptw2Bc3vEe8FpgHve0AWRTtBoxFc1j+GAB8dW1LzG1XgnkMqeOTdE2z4Xyg/LR7Ivk1DKMg2Mv7a8xgyocTt8zr2dn9fubT2kLjUAJcDfwCuBvwI6rSTq2cXMHi8D0LAtp8alrMuPaICGdcrkhC/tpO0ooVR+txMwgN9mTEiiR/2nuZ4rvN1w9X1uDA63mV9FAIM/hKDv5lnPnnT5bDxV/96vl4Z1dG1zjbKVmMRuupHEqtiVHJv/vOBqi7NxSNzxEKINcCj7rag1bxSx2H7yQLWHs5xO7OLqUqw/ecg1nzViawTvhj8LfQ5r4z+o23LJPyDas+xFedlM/emcTg/QRxa3aVYzGOxpSOMBfIQms+R1g84f3KfJnlld/LXc2GfKDqHt6wjljt4Y474bPRZ6bJzpJRsPJbH78fza7Y1x1mxOUFlPMmCF++nS2ICjz14L4sWfEhmZiZLly5tdTnaI21ijlgIMROYCZCUlNRap1W0MEOTOhFs0PHjviy3HFX0PpKRlxcz4rJiDm/zZ+faIPZtCmDnr8EIjSQi1kRUYhVRiVWERZvY+euvIO8kMLQT5SUarJZEhKY7Pr7JVFbEYzEDFAA/ERj6IwPPr+L8q65l7dIu7Fy7gktvuc/tOWgfnYYRXcMYmtQJbQvFwe4IKF1uHCEEY3pEEBti4OeD2ZQYzQ34U9SPzV6XxkadWhIhbEu1+sQEcd/PP9RM0YwZ0WZ9CNsdjRpiIcRqwFl0+qeklN+4eyIp5XxgPtha0W5LqGjz9IgKIsTPh293ZVBc4XpuzBEhoPfwcnoPL8dqhZOHDBzcEsDpVB+yT/mwb5Mem/fn3cDdlBZW/7IIaT1GQEg6naKXE514gonTR/L7D19TnJ/D9Q/Zhtf0vgYqSorcXqvcKzqIcb0iCDJ07GFoT+iz0mX36RYZSHwnPzYczWW3wK21wc5ozcQOYItmFxtqoHtkAL2igzq8XnibRg2xlPLi1hBE0b6JDPJlxohEvtudSbqLtcau0GigS18jXfqeCRhSmJPN1/Pmc2DrBsxVVSD0hMcEc92f/o+9G+unH0zsaXs5uQpw4Gr4LyLQh/G9okgK9z/bS26XKH1ufXx1Wib2iaZ3TDCLygparGfb3Ih1YQE+xIX60SXcn8QwfxWIoxVRy5cUHsPfR8e1QxP47VgeW0/ku5w3dueFERoZRVAnExbTcXQ+tnmx3sOm03vYaHoPG+1SBneH/wx6LaO7h5McH9Ji6RgVCkfiQ/34/ZcfOJlXzpbUfOK79/Fo+WcTsc6g1xId7EtUkIGYEAPxoX74+SjD6y2au3zpGmyL2SKB5UKInVLKyzwimaJdotUIxvaMoHO4Pz/uO+00I5K7L4ymzIs1th5SqxEMjA9hdPdw1eKvg9Ln1iEp3J+kcH+yio0cPF3CsexSityc0nFGQ6NAb36/pyaJSkSgL+GBtv/n0EqAdoGKrKVoMYwmCz8dyOZwlm2JUmuFuHQWGvDOuW/TOzqI0d3DCfVvP1F/VGStc4OckkpS88rIKakkr7SSgnJTg/GrqxECKovz+Ob9l9m+dhWVxgoMBj8umzyFV195jZ5dE1SAjTZCm/CaVpx7GPRarkyOZUBeMGsO5TTLa/RscJw7vu5Pz9E1IoDR3cKJClYZYBRtk8ggXyKDzoSmtFglpZVmqsxWqizWmhUJOo0tJahOK/D30eGv16LR9CJ1VRybV1diMBioqqokLjKMXt0SvXU5irNEGWJFi9M5PIBbRvmzMz6Y1YuCPBri0hUaIegdE8iwzmG1XnAKRXtAqxFnNXyclZXFrFmzmDlzJvPnzyczM7MFpVN4GmWIFa2CViMY1jmMcE0F1958B/0mTuOnrxd5fD1kkEFHn5hgBiaEqHkwxTmDY1CNefPmeVESRVNQhljRqixb9jVgG3qbfvkFHMgsJq2gwq35MFcY9Fo6h/vTLzaYpDB/5QWtUCjaFcoQK7yCViPoGxtM39hgKs0WTuaVk5JbRm5pJYXlJpdRujRCEGTQ0SlAT1yIH53DA4gO9lUOKQqFot2iDLHC6/jqtPSMDqJndFDNttJKMyVGU81aZCHAoNMS7KdX4ScVCkWHQhliRZsk0FdHoK+qngqFouPTsROsKhQKhULRxlGGWKFQKBQKL6IMsUKhUCgUXkQZYoVCoVAovIgyxAqFQqFQeBFliBUKhUKh8CLKECsUCoVC4UWUIVYoFAqFwosoQ6xQKBQKhRdRhlihUCgUCi+iDLFCoVAoFF6kWYZYCPGaEOKgEGK3EOJrIUSoh+RSKBStjNJnhcI7NLdHvAoYIKVMBg4DTzZfJIVC4SWUPisUXqBZhlhKuVJKabZ/3QQkNF8khULhDZQ+KxTewZNzxHcCKzxYnkKh8B5KnxWKVqLRhK9CiNVAjJNdT0kpv7Ef8xRgBhY2UM5MYCZAUlJSk4RVKBTNwxP6rHRZofAsjRpiKeXFDe0XQtwOTAYuklLKBsqZD8wHGD58uMvjFApFy+EJfVa6rFB4lkYNcUMIIS4H5gDjpZTlnhFJoVB4A6XPCoV3aO4c8dtAELBKCLFTCPGeB2RSKBTeQemzQuEFmtUjllL28JQgCoXCuyh9Vii8g4qspVAoFAqFF1GGWKFQKBQKL6IMsUKhUCgUXkQZYoVCoVAovIgyxAqFQqFQeBFliBUKhUKh8CLKECsUCoVC4UWUIVYoFAqFwosoQ6xQKBQKhRdRhlihUCgUCi+iDLFCoVAoFF5EGWKFQqFQKLyIMsQKhUKhUHgRZYgVCoVCofAiyhArFAqFQuFFlCFWKBQKhcKLKEOsUCgUCoUXUYZYoVAoFAov0ixDLIR4UQixWwixUwixUggR5ynBFApF66L0WaHwDs3tEb8mpUyWUg4GvgOebb5ICoXCSyh9Vii8QLMMsZSy2OFrACCbJ45CofAWSp8VCu+ga24BQoi/AH8AioALmy2RQqHwGkqfFYrWR0jZcKNXCLEaiHGy6ykp5TcOxz0JGKSUz7koZyYw0/61N3CoEdkigNxGjmltlEzuoWRyD3dk6iyljPTUCT2hz03QZWi/9781aWvygJLJXZqly40aYncRQiQB30spB3iovK1SyuGeKMtTKJncQ8nkHm1RpmqUPrc+bU0eUDK5S3Nlaq7XdE+Hr1OBg80pT6FQeA+lzwqFd2juHPHLQojegBU4AcxqvkgKhcJLKH1WKLxAswyxlPJaTwnihPktWHZTUTK5h5LJPdqUTEqfvU5bkweUTO7SLJk8NkesUCgUCoXi7FEhLhUKhUKh8CJeN8RCiMuFEIeEEEeFEE842e8rhFhi379ZCNGlDcj0sBBivz0c4E9CiM7elsnhuGuFEFII0eJehe7IJIS4wX6v9gkhFnlTHiFEkhDiFyHEDvuzm9SS8tjP+V8hRLYQYq+L/UII8U+7zLuFEENbWqaWQumyZ2RyOO6c1WV3ZGptfW5RXZZSeu0P0ALHgG6AD7AL6FfnmHuB9+yfpwNL2oBMFwL+9s+z24JM9uOCgLXAJmC4t2UCegI7gE7271Felmc+MNv+uR+Q2pL3yH6eccBQYK+L/ZOAFYAARgGbW1omL95/pctKlz0pU6vqc0vqsrd7xCOAo1LKFCllFbAY27IJR6YCC+yfvwQuEkIIb8okpfxFSllu/7oJSGhBedySyc6LwCuAsYXlcVemu4F5UsoCAClltpflkUCw/XMIkNGC8thOKOVaIL+BQ6YCH0sbm4BQIURsS8vVAihd9pBMds5lXXZXplbV55bUZW8b4njglMP3NPs2p8dIKc3YQu+Fe1kmR+7C1gpqSRqVyT4MkiilXN7CsrgtE9AL6CWE2CCE2CSEuNzL8swFbhFCpAHfAw+0oDzucrb1ra2idNk9lC57Tqa5tC19brIuNzvW9LmMEOIWYDgw3styaIA3gdu9KYcTdNiGtCZg62msFUIMlFIWekmeGcBHUso3hBCjgU+EEAOklFYvyaNoIyhdbpS2psvQgfTZ2z3idCDR4XuCfZvTY4QQOmxDEHlelgkhxMXAU8AUKWVlC8rjjkxBwABgjRAiFdv8xLct7OThzn1KA76VUpqklMeBw9iU2Vvy3AV8DiCl/A0wYIsR603cqm/tAKXLnpFJ6bL7MrU1fW66Lrfk5LYbk986IAXoypkJ+f51jrmP2g4en7cBmYZgcyTo2VbuU53j19DyDh7u3KfLgQX2zxHYhm3CvSjPCuB2++e+2OaURCs8vy64dvC4ktoOHr+3Rp3y0v1Xuqx02ZMytbo+t5Qut3jFc+PCJmFrXR3DlgEG4AVsrVOwtXK+AI4CvwPd2oBMq4EsYKf971tvy1Tn2BZXXjfvk8A2zLYf2ANM97I8/YANdqXeCVzaCvfoMyATMGHrVdyFLXTkLId7NM8u857WeG5evP9Kl92Qqc6x56QuuylTq+pzS+qyiqylUCgUCoUX8fYcsUKhUCgU5zTKECsUCoVC4UWUIVYoFAqFwosoQ6xQKBQKhRdRhlihUCgUCi+iDLFCoVAoFF5EGWKFQqFQKLyIMsQKhUKhUHiR/weXjBCog95gEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initialize plots\n",
    "f, (y1_ax, y2_ax) = plt.subplots(1, 2, figsize=(8, 3))\n",
    "\n",
    "# Test points every 0.02 in [0,1]\n",
    "test_x = torch.linspace(0, 1, 51)\n",
    "test_i_task1 = torch.full((test_x.shape[0],1), dtype=torch.long, fill_value=0)\n",
    "test_i_task2 = torch.full((test_x.shape[0],1), dtype=torch.long, fill_value=1)\n",
    "\n",
    "# Make predictions - one task at a time\n",
    "# We control the task we cae about using the indices\n",
    "\n",
    "# The gpytorch.settings.fast_pred_var flag activates LOVE (for fast variances)\n",
    "# See https://arxiv.org/abs/1803.06058\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    observed_pred_y1 = likelihood(model(test_x, test_i_task1))\n",
    "    observed_pred_y2 = likelihood(model(test_x, test_i_task2))\n",
    "\n",
    "# Define plotting function\n",
    "def ax_plot(ax, train_y, train_x, rand_var, title):\n",
    "    # Get lower and upper confidence bounds\n",
    "    lower, upper = rand_var.confidence_region()\n",
    "    # Plot training data as black stars\n",
    "    ax.plot(train_x.detach().numpy(), train_y.detach().numpy(), 'k*')\n",
    "    # Predictive mean as blue line\n",
    "    ax.plot(test_x.detach().numpy(), rand_var.mean.detach().numpy(), 'b')\n",
    "    # Shade in confidence \n",
    "    ax.fill_between(test_x.detach().numpy(), lower.detach().numpy(), upper.detach().numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "    ax.set_title(title)\n",
    "\n",
    "# Plot both tasks\n",
    "ax_plot(y1_ax, train_y1, train_x1, observed_pred_y1, 'Observed Values (Likelihood)')\n",
    "ax_plot(y2_ax, train_y2, train_x2, observed_pred_y2, 'Observed Values (Likelihood)')"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "anaconda-cloud": {},
  "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}