{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QEP Regression with a Spectral Mixture Kernel\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This example shows how to use a `SpectralMixtureKernel` module on an `ExactQEP` model. This module is designed for\n",
    "\n",
    "- When you want to use exact inference (e.g. for regression)\n",
    "- When you want to use a more sophisticated kernel than RBF\n",
    "\n",
    "The Spectral Mixture (SM) kernel was invented and discussed in [Wilson et al., 2013](https://arxiv.org/pdf/1302.4245.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import qpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we set up the training data for this example. We'll be using 15 regularly spaced points on [0,1] which we evaluate the function on and add Gaussian noise to get the training labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 15)\n",
    "train_y = torch.sin(train_x * (2 * math.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the model\n",
    "\n",
    "The model should be very similar to the `ExactQEP` model in the [simple regression example](./Simple_QEP_Regression.ipynb).\n",
    "\n",
    "The only difference is here, we're using a more complex kernel (the `SpectralMixtureKernel`). This kernel requires careful initialization to work properly. To that end, in the model `__init__` function, we call\n",
    "\n",
    "```\n",
    "self.covar_module = qpytorch.kernels.SpectralMixtureKernel(n_mixtures=4)\n",
    "self.covar_module.initialize_from_data(train_x, train_y)\n",
    "```\n",
    "\n",
    "This ensures that, when we perform optimization to learn kernel hyperparameters, we will be starting from a reasonable initialization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "POWER = 1.0\n",
    "class SpectralMixtureQEPModel(qpytorch.models.ExactQEP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(SpectralMixtureQEPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.power = torch.tensor(POWER)\n",
    "        self.mean_module = qpytorch.means.ConstantMean()\n",
    "        self.covar_module = qpytorch.kernels.SpectralMixtureKernel(num_mixtures=4)\n",
    "        self.covar_module.initialize_from_data(train_x, train_y)\n",
    "\n",
    "    def forward(self,x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return qpytorch.distributions.MultivariateQExponential(mean_x, covar_x, power=self.power)\n",
    "\n",
    "    \n",
    "likelihood = qpytorch.likelihoods.QExponentialLikelihood(power=torch.tensor(POWER))\n",
    "model = SpectralMixtureQEPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we handle using Type-II MLE to train the hyperparameters of the q-exponential process.\n",
    "The spectral mixture kernel's hyperparameters start from what was specified in `initialize_from_data`.\n",
    "\n",
    "See the [simple regression example](./Simple_QEP_Regression.ipynb) for more info on this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/100 - Loss: 1.536\n",
      "Iter 2/100 - Loss: 1.511\n",
      "Iter 3/100 - Loss: 1.479\n",
      "Iter 4/100 - Loss: 1.450\n",
      "Iter 5/100 - Loss: 1.418\n",
      "Iter 6/100 - Loss: 1.383\n",
      "Iter 7/100 - Loss: 1.348\n",
      "Iter 8/100 - Loss: 1.312\n",
      "Iter 9/100 - Loss: 1.277\n",
      "Iter 10/100 - Loss: 1.241\n",
      "Iter 11/100 - Loss: 1.206\n",
      "Iter 12/100 - Loss: 1.172\n",
      "Iter 13/100 - Loss: 1.139\n",
      "Iter 14/100 - Loss: 1.105\n",
      "Iter 15/100 - Loss: 1.069\n",
      "Iter 16/100 - Loss: 1.032\n",
      "Iter 17/100 - Loss: 0.994\n",
      "Iter 18/100 - Loss: 0.956\n",
      "Iter 19/100 - Loss: 0.919\n",
      "Iter 20/100 - Loss: 0.884\n",
      "Iter 21/100 - Loss: 0.850\n",
      "Iter 22/100 - Loss: 0.817\n",
      "Iter 23/100 - Loss: 0.783\n",
      "Iter 24/100 - Loss: 0.748\n",
      "Iter 25/100 - Loss: 0.712\n",
      "Iter 26/100 - Loss: 0.677\n",
      "Iter 27/100 - Loss: 0.641\n",
      "Iter 28/100 - Loss: 0.607\n",
      "Iter 29/100 - Loss: 0.573\n",
      "Iter 30/100 - Loss: 0.540\n",
      "Iter 31/100 - Loss: 0.507\n",
      "Iter 32/100 - Loss: 0.474\n",
      "Iter 33/100 - Loss: 0.441\n",
      "Iter 34/100 - Loss: 0.407\n",
      "Iter 35/100 - Loss: 0.374\n",
      "Iter 36/100 - Loss: 0.342\n",
      "Iter 37/100 - Loss: 0.310\n",
      "Iter 38/100 - Loss: 0.278\n",
      "Iter 39/100 - Loss: 0.247\n",
      "Iter 40/100 - Loss: 0.215\n",
      "Iter 41/100 - Loss: 0.183\n",
      "Iter 42/100 - Loss: 0.152\n",
      "Iter 43/100 - Loss: 0.121\n",
      "Iter 44/100 - Loss: 0.090\n",
      "Iter 45/100 - Loss: 0.059\n",
      "Iter 46/100 - Loss: 0.029\n",
      "Iter 47/100 - Loss: -0.002\n",
      "Iter 48/100 - Loss: -0.033\n",
      "Iter 49/100 - Loss: -0.064\n",
      "Iter 50/100 - Loss: -0.094\n",
      "Iter 51/100 - Loss: -0.124\n",
      "Iter 52/100 - Loss: -0.155\n",
      "Iter 53/100 - Loss: -0.185\n",
      "Iter 54/100 - Loss: -0.216\n",
      "Iter 55/100 - Loss: -0.246\n",
      "Iter 56/100 - Loss: -0.276\n",
      "Iter 57/100 - Loss: -0.306\n",
      "Iter 58/100 - Loss: -0.337\n",
      "Iter 59/100 - Loss: -0.367\n",
      "Iter 60/100 - Loss: -0.399\n",
      "Iter 61/100 - Loss: -0.431\n",
      "Iter 62/100 - Loss: -0.466\n",
      "Iter 63/100 - Loss: -0.504\n",
      "Iter 64/100 - Loss: -0.546\n",
      "Iter 65/100 - Loss: -0.594\n",
      "Iter 66/100 - Loss: -0.649\n",
      "Iter 67/100 - Loss: -0.708\n",
      "Iter 68/100 - Loss: -0.775\n",
      "Iter 69/100 - Loss: -0.852\n",
      "Iter 70/100 - Loss: -0.943\n",
      "Iter 71/100 - Loss: -1.048\n",
      "Iter 72/100 - Loss: -1.165\n",
      "Iter 73/100 - Loss: -1.279\n",
      "Iter 74/100 - Loss: -1.113\n",
      "Iter 75/100 - Loss: -1.400\n",
      "Iter 76/100 - Loss: -1.351\n",
      "Iter 77/100 - Loss: -1.625\n",
      "Iter 78/100 - Loss: -1.539\n",
      "Iter 79/100 - Loss: -1.628\n",
      "Iter 80/100 - Loss: -1.638\n",
      "Iter 81/100 - Loss: -1.833\n",
      "Iter 82/100 - Loss: -1.693\n",
      "Iter 83/100 - Loss: -1.681\n",
      "Iter 84/100 - Loss: -1.977\n",
      "Iter 85/100 - Loss: -1.805\n",
      "Iter 86/100 - Loss: -1.735\n",
      "Iter 87/100 - Loss: -1.903\n",
      "Iter 88/100 - Loss: -1.983\n",
      "Iter 89/100 - Loss: -2.066\n",
      "Iter 90/100 - Loss: -1.982\n",
      "Iter 91/100 - Loss: -1.971\n",
      "Iter 92/100 - Loss: -2.196\n",
      "Iter 93/100 - Loss: -2.079\n",
      "Iter 94/100 - Loss: -2.219\n",
      "Iter 95/100 - Loss: -2.051\n",
      "Iter 96/100 - Loss: -2.102\n",
      "Iter 97/100 - Loss: -2.262\n",
      "Iter 98/100 - Loss: -2.234\n",
      "Iter 99/100 - Loss: -2.318\n",
      "Iter 100/100 - Loss: -2.285\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iter = 2 if smoke_test else 100\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)\n",
    "\n",
    "# \"Loss\" for QEPs - the marginal log likelihood\n",
    "mll = qpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "for i in range(training_iter):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x)\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iter, loss.item()))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've learned good hyperparameters, it's time to use our model to make predictions. The spectral mixture kernel is especially good at extrapolation. To that end, we'll see how well the model extrapolates past the interval `[0, 1]`.\n",
    "\n",
    "In the next cell, we plot the mean and confidence region of the q-exponential process model. The `confidence_region` method is a helper method that returns 2 standard deviations above and below the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEYCAYAAABxx2wUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYSJJREFUeJztnXd8W+X1/99X0/KS99527GwncUhImIGwygpQVlMKlFLKL6wCpdB+25Tvt23aQoFCmaWFFkoZhTASCATIICRkb7Ls2I7jJLYT76V5f3/Iku14Sda9kuw879dLBEl3PFeWPvc855znHEmWZRmBQCAQhByaYA9AIBAIBP0jBFogEAhCFCHQAoFAEKIIgRYIBIIQRQi0QCAQhChCoAUCgSBEEQItEAgEIYoQaIFAIAhRhEALBAJBiCIEWiAQCEIUVQX6+eefZ/LkyURHRxMdHc2sWbP45JNP1DylQCAQjBokNWtxfPTRR2i1WsaMGYMsy/zzn//kscceY+vWrUyYMEGt0woEAsGoQFWB7o+4uDgee+wxbrvttkCeViAQCEYcukCdyOFw8M4779DW1sasWbP63cZisWCxWDzPnU4n9fX1xMfHI0lSoIYqEAgEqiHLMi0tLaSlpaHRDOFlllVmx44dckREhKzVamWz2SwvXbp0wG0XLlwoA+IhHuIhHqP+UVVVNaR+qu7isFqtHDp0iKamJv773//y8ssvs2rVKsaPH99n25Mt6KamJrKysqiqqiI6OlrNYQoEIcmGg/VsrKhX7fhXTkkjIy5c0WPWtXTy9sbDih4zVPn+6dmYw/U+7dPc3ExmZiaNjY2YzeZBt1XdxWEwGCgoKACgpKSEjRs38pe//IUXX3yxz7ZGoxGj0djndXcWiEBwqtGpaSUsIlK147c49Yr/tvbV21UdcygRHR1NtI8C7cYbt23A86CdTmcvK1kgEAzMiVarqsc/1typ+DErT7QrfsxTFVUt6EceeYRLLrmErKwsWlpaeOONN1i5ciWffvqpmqcVCEYFdoeTxnabquc42qSsQFvtTo40dih6zFMZVQW6traWH/zgBxw9ehSz2czkyZP59NNPueCCC9Q8rUAwKmhot+FUOQu2w+qgsd1KTLhBkeMdbmjH4VR3zKcSqgr03//+dzUPLxAMG6fTidWqrvvAX+oaWzFpHKqfp6KmkbGpyvihq443BWTMoYLV0knnSder1+vRarWKHD9gedACQahgtVopLy/H6XQGeyiDYrU5mBqr/hgdzTWUd55Q5FjhFjtTY08dC7r2SBXHNX2DfTExMaSkpPi9fkMItOCUQpZljh49ilarJTMzc+iFAkGkucOK1REYsYuLMKDxU0zsTvV95qFGjEmPTtv9HZJlmfb2dmprawFITU316/hCoAWnFHa7nfb2dtLS0ggPVzb/V2labRJ6bWAEWqPTE6b3b1rebrWjN4TuDU8NwsIMvQQawGQyAa4YXFJSkl/ujlPr0xScUsiyjMXe2z/ocLieGwzKBMXUwinLOAJYJsdi99+VYlXgGKMF983fZvNvRiEsaMGo5XBDB+9uOUx0mJ6EKCMJkQbijBJOWQ752i6BzoSw2B3Ism7Yn4ssy1gdQqDdKPX9EgItGLWU1rUiy9DUYaOpw0ZZLZg0DqbHh36WgT3AYifLYHM4MeiGNx23OWQCWxfz1EC4OASjlrLa1n5fd8hyyOfq2v0Y3/RJRbz03DM+7+ePm0NN98Zwr2c0IARaMCqpbe6kpdPe/5syWO2hbUX3J9DVh6u4b8EdFBflkpkQTcnEQv7n5w9QX69MipxfAu0I3uf52KLfkmI2kWI2kR4XyfjcDOZdMpeXnnvG57ISX3+1mhSziabGRnUG6yNCoAWjktK6/q1nN0r4Szdt2sR5553Hpk2b/D7WydhPSq+rLC/nonPPpLyslOf//i/Wbd3Fn558hq9WreSyuefSUO9/xTuHUx6Wa6XTasNiC+4Nr2jceHbsL2fz7v28u2QZl827mqefeJzLL5hDa0tLUMfmD0KgBaOSgdwbbmwO2e9l1P/6179YsWIFr732ml/HORmns+/YHn7wPgwGPW8uXsLsM88iIzOL8y+4iHc+WMrRo0dY9H+/6bV9a2srP/nhD8hNjWfK2Dz+8bcXPO/Jssxji35LyYQxZCWaKS7K5ZcP3Q+4rGiLxcKDDz5Ieno6ERERzJw5k5UrV3r2f/XVV4mJieHDDz9k/PjxRIab+Pc/XyE7KaaP5fk/P3+Aay672PN8/bqvufLi88lJjmXa+AJ++dD9tLW1ed6vq6vlpuuvISc5ltMmjeXdt//j1Wem0+lISk4hJTWNcRMm8qM7/h+LP/6MvXt289en/uzZ7p033+DCc84gPz2RSWNyuPO2m6mrc+UsH6qs5JrLLgKgKDuVFLOJe+68HYAvP/+MKy46j8KsFMblpPP9666m4uBBr8bmD0KgBaOOxnYrx72oAjccv2llZSWbN29my5YtvPXWWwC8+eabbNmyhc2bN1NZWenzMU/mZPdGQ309K79Yzi233eHJsXWTlJzCNdfewIfv/Zeepd2fe/pJJkycxOdffcNdP32QX/38QVZ9+QUASz5YzEvPPcOfnvora7fs5JU33mbc+ImAS6Dvuusu1q1bx5tvvsmOHTu49tprufjiizlw4IDn+O3t7fzxj3/khRdfYtX6LVx93Q1Em2NY+uH7nm0cDgcfvPcu11x3AwAVBw9y4zVXcukV8/hy7UZefOU1Nqxbxy9+9lPPPvfe+WOOVB/m3SXLePlfb/Dqyy9xoq5uWJ/jmMIizrvgQpZ+9IHnNZvNxs//59d8uWYDr77xFlWHKrn3zh8DkJ6Rwd9fc90Qvt68gx37y/ntHx53XW9bG3csuIdPV3zNOx9+jEaj4dbvX6/6alSRxSEYdZQN4d5wY7E5fF6ckZOT4/l/dypVXV0dJSUlntf97YFhP+lHX36wFFmWGVNU1O/2Y4qKaGxs4PjxOhITkwA4beYs7r7/ZwDkF4xh4zfrePG5ZzjnvPOpPlxFUlIyZ597Hnq9nozMLKaVnOY6V0UFr7zyCocOHSItLQ2ABx98kGXLlvHKK6/w+9//HnAJ3XPPPUdu4Xg6utwb8675Lu+98xbf+8EtAHy1cgXNTY1cesU8AJ5+4jGuufYGfvz/7gYgL7+A3/7pca76zoX88YmnqT5cxZfLP+WTL79iasl0AJ746wucddqUYX+WBWOKPDcmgO/ddLPn/7Nzc/ntH//MxXPOpK21lYjISGJi4wBISEjEHBPj2fayK6/qddwnn32BCXmZfPvtt0wpnjzs8Q2FsKAFo46y2rahN8JlLfoqpq+//jo6ncuuce/r/len0/H666/7dLz+GCiDw5exTp8xs8/zA/v2AnD5vKvp7OxgZvE4Hrj7//HxRx9gt7sCqnu+3Y3D4aCwsJDIyEjPY9WqVZSVlXmOZzAYmDBxkkecAa6+9gbWrlnNsaNHAHj3nTeZe+HFHqHbvWsHb73xGnlpCZ7HDVdfgdPp5FBlBQf27UWn01E8dZrnmGMKizCbY7y+7pORZRl65CRv37qFm66/hpIJY8hPT+SqSy8E4PDhqkGPc7CslJ/88AfMmDyOgowkTps0FoCqQ4eGPTZvEBa0YFTRbrVzpMm7esQyLjeH0Qcrev78+YwbN66Xxexm/fr1TJs2rZ+9fOPkAGFOXj6SJHFg3z64vO/2B/btIyYmloSERK+On56RyZpNO/hq5ZesWvEFDz9wH889/SSLP15OW2srWq2WzZs391miHBnZ3SXFZDL1EmeAqSXTycnN4/133+Hm237MJ0s+5C/PveR5v72tjZtuvY0f3bGg75gyMzlYeqDP6/5yYP8+srJzAGhra+PGq6/g3PPn8uzfXiE+IZHqqipuuPpybENUNvzB9deQkZnFn59+juTUVJxOJ+eeXoLVpm5FRCHQglFFWW2bTwsmLD4KdE80Gg1Op9Pzr1I4TjpWXFw858w5n1f//iI/XnB3Lz90bc0x3n3nTa69YX6v1WubN27odYzNGzcwpmis57nJZOLCSy7lwksu5dbbf8KZ04vZs3sXk4qn4HA4qKmp4eyzzx50nB3WvpkbV193A++9/RapaeloNBrmXnSJ571JxVPYv3cvufn5/R6voLAIu93O9q1bPC6O0gP7aWpqHHQcA3Fg/z5WfP6Zx9VTun8f9fUn+OVv/o/0jEzAZVH3xGBwta9yOLuvrb7+BKUH9vP4089y+uwzAVewMxAIF4dgVOGt/9nNcNwcSUlJpKSkUFJSwgsvvEBJSQkpKSkkJSX5dJz+cDhl+vNw/P7xJ7FYrNx49eWs+3qNy1/7+WdcN+8yUlPTeORXv+m1/cb16/jrU3+mrPQA//jbC3z0/nvc/hOX5frmv1/jjX+9yp5vd1NZXs67b/0Hk8lERlYW+QVjuOa6G/jBzTfz3nvvUV5ezoYNG1i0aBFLly7tdY7+PrVrrr2BHdu38pc//4nLrriqV4/Ru+57gE0bvuGRB+9j147tHCwrZdnSj3jkwfsAKBhTyJy5F/LQfXezZdMGtm/dwgN339knMNofdrud2ppjHDt6hD27d/Hyi89x1XcuZMKkySy4xxWETM/MxGAw8PcXn6eyvJxPP17Ck39a1Os4GZlZSJLE8mWfcPx4HW2trcTExBIXF8/rr/6D8rIy1qxaycJf/HzIMSmBEGjBqMFid1BV71s/PKcsY/OxpGdGRgYVFRWsX7+eO+64g/Xr11NRUUFGRoZPx+mPgfKQ8/IL+HTlGrJycvnxLd/n9CkT+Nm9CzjjrLNZ8vlKYuPiem3/k7vuZfvWLcw963SeeuyPPPr7PzJnrquTkdls5vV//oMrLjqPOWecxuqVX/KvN98lLi4egKeee4kbbpzPAw88QFFREfPmzWPjxo1kZWUBrs9soE8sNz+fqSXT+XbXTq7uyt5wM37iJN5b+hkHS0u58pK5zD3rdP70+/8jJSXNs81fnnuR5NRUrvrOhfzwphv4/i23EZ84tOtm355vmVyYS8mEQq6+9CI+Wvwu99z/IB8s+4KILtdMQkIif3n+JT56/z3OnjmVZ558nIW/7S3QqWnp/OwXv+J3v/kVkwqyeeRnP0Wj0fDCP/7F9m1bOXdWCb/+xUP8+v8W9TcMxZFkf0POKtLc3IzZbKapqUl09RYMyf6aFpbuODroNiaNg6mxNtKzstEbXNZduEFLVNjwOjMrTZvFTqtlgBWQAUQjQVSYHqNO06fwT6iMMRSIj+hbbhSgs7OT8vJycnNzCQsL6/WeL7omfNCCUcNQi1MGwmJ3EqXwWIaLLUQqwjm7ikxpJAmTXovJoEWrkVwF6fvxPQvUQQi0YFTgcMqUn/Auva6/fe0OZ7+WUKDxp0iSGjhlmTarnTarHaNOg1Yjqd7IVtCNEGjBqKCqvh2LbfjWp8UefIF2hniVPSWK+gt8I/gmg0CgAAePD8+94ebkzivB4OT8Z4FACLRgVHCksdOv/W2O4FuvJy/xFgiEQAtGPA6nTH2b/yu6gt2yydd0P8HoRwi0YMRT32ZVxPoNdgZFoNtcCUIfIdCCEU9di29dMwbCFsQgmCzLIZfBIQg+QqAFI57jrcoItN0p+10q1J9zCwQnIwRaMOJRSqAheH5g4d4Q9IeqAr1o0SJOO+00oqKiSEpKYt68eezbt0/NUwpOQZRycUDw/NBDWdD33Hk7KWYTD913d5/3Hn7gvl7tmQSjB1UFetWqVSxYsIBvvvmG5cuXY7PZuPDCC3v1IBMI/KHNYld06XGwBNobyz09I4P333uHjo7uetednZ0sfuct0jMz1RyeIEioKtDLli3jlltuYcKECRQXF/Pqq69y6NAhNm/erOZpBacQSlrPEBwXhytAOPSNYVLxFNLSM/j4o/c9r3380fukZ2YyaXKx5zWn08nTf36M0yaNJSc5lvPOmMFH77/ned/hcPDTBT/xvH9GyWT+9vxfe53rnjtv55bvXctzTz/J5MJcxuWk8/AD92Gz2fy/YIHXBHSpd1NTEwBxJ5VGdGOxWLBYun9wzc3NARmXYOTir/9ZlqG9V4VSGSNOtBr1wzPh4a5uTA6n7HWTgRu//wPefP01rrnuRgD+89q/uGH+Taxds9qzzdN/fox33/4Pf3ryGfLyC1i3dg13/fiHxCckMvvMs3A6naSmp/O3f/6b2Lh4Nm1Yx4P33kVScgpXXv1dz3G+/mo1ScmpvLtkGeUHy7jj1puYOGky37/lh4p+DoKBCZhAO51O7rvvPs444wwmTpzY7zaLFi3i0UcfDdSQBKMAfwW6vR3y08KG3lAFWlshIsK3DI5rrr+R3z/6a6oOubqHb1y/jhde+ZdHoC0WC3954k+888FSps84HXA1R92wbi2vvfIys888C71ez0O/+JXnmNk5OWzasJ4PF7/bS6DNMTEsevxJtFotYwqLmHvhxXy1aoUQ6AASMIFesGABu3btYs2aNQNu88gjj3D//fd7njc3N5MpfGuCQVDaxREMfMngSEhIZO6FF/PWG68jyzLnX3gx8fEJnvfLD5bR0d7OdfMu67WfzWplYg83yD/+9gJvvvYvDh+uorOzA5vVyoRJvbtTF40d36svYXJKCnt27/b18gR+EBCBvuuuu1iyZAmrV68etOuE0Wjs1SJHIBgMu8NJfZt/PtHwcCg70ruOh14rEReh/vcwPNz1r6850DfcdDO/eNDVxmnRn5/q9V57m6to1OtvLyY1Na3XewajAYD3//s2//s/j7Dwt39g+oyZREZG8dzTT7Jl88Ze2+v1veVBkiScskgHDCSqCrQsy9x9990sXryYlStXkpubq+bpBKcY9W1Wv2sTS5LLzdAbmfBwuU8nEbXwNTB53twLsdmsSJLEnPMv6PVeYdE4jEYj1YermH3mWf3uv2H9OqbPOJ1bb7/D81pF+UHfBy5QHVUFesGCBbzxxht88MEHREVFcezYMcDVE82bRpACwWDUKbhA5WRsDhmDTn2BdjWJ9U2gtVotX23Y5vn/nkRGRXHn3fex8JGHcDqdzDx9Ns3NTWxYv46oqGiu/973ycsv4J0332DF58vJysnhv2++wbatm8nKzlHoqgRKoapAP//88wCce+65vV5/5ZVXuOWWW9Q8teAU4Hir/xXsBsLmcGLQqZ/JMdwSo1GD9LL7+f8sJD4hgWeeeIwHK8qJNscwuXgK9zzwEAA33fojdu7Yzh0/vAkJiXnfvY5bbvsxX37+2bDGIlAP0TRWMGL57+bDPnfx7q9pbH+E6TSYww3+DnFIRAPWkY3aTWNFLY5TmGAXqPcXJWtwnIw1QAtWRA0OwWAIgT5FcTpl3t18eMSmqbV02uhQsbt0oPoD2kb4TVKgLkKgT1G2HW6kurGDj7YfodMW/H58vqKm/9mN2nU5Qr1JrCD4CIE+BWm12FlXdgKApg4bH+88inOECUUgLH+1BVo0iRUMhRDoU5DV++uw9ugeUnmina/LjgdxRL6jpv/ZjdqFk0STWMFQCIE+xTh0op19x1r6vL6poqHf10OVQAi03eFUtcOKaBIrGAoh0KcQDqfMin21A76//Ntj1LZ0Dvh+qGB3OGnwc4m3N8io24rKGsQeiIKRgRDoU4hNFfXUtw0cXLM5ZD7afjTkg4YnFFji7S1q+aHtDmfArkEwchECfYrQ1GFjY0X9kNs1d9hC3tURyNRAtQRaWM8CbxACfYqwcl+t1z7PihOh3ZJMzRocJ2Ozq2PlWkN4gYosyzx4zwLGZqeRYjaxa8d2rrr0Qn718IOD7jd9UhEvPfdMgEZ5ahDQjiqC4FDfZuVgnfeie7ihA4dTRqsJTDU3XzmuggX9t9UDV3MzGbRoFK1sJzN/Zvaw9qytOcZTj/+Rzz9dxrGjR0hITGTCpGJ+fOddnHXuHEVG9+Xnn/HWG6/x3tJPyc7JJS4+gX+89iZ6vV6R4wu8Rwj0KcDhBt/qVVjtTqobOsiKD1dpRP4RSAsaXMFVjVY5gXY4ZYZjlx+qrOSKi+YQbY7h1//3e8ZNmIjNZmPlF8t55MH7WLNpuyLjqyg/SFJKCqfNnOV5LXaANnUCdREujlOAww0dQ290EqHq5mjutGGxBdY9oPRqv+Ee7+EH7kWSJD758isuu/Iq8gvGMHbceH5y170s/XwVAIerDnHzjdeSl5ZAQUYSt988n7raGs8xHlv0W84/cybvvPkG0ycVMSYzmTtuvYnWFlfc4Z47b+eXP7uf6qoqUswmpk8qAujj4qirq+Wm668hJzmW0yaN5d23/9NnvE2Njdx/152Mz8ukICOJay67mN07d3g9FnC1yvvrU3/m9CkTyEo0UzJhDE899kfP+9WHq7j95vkUZqUwNjuNm2+8lkOVlcP6fEMRIdCnANWjSKCDUTvEJajKifRwBLqhvp4Vn3/GrT/6CRF9OwxgjonB6XRyy43X0dhQz+Kln/H2+0uorKjgx7fe1GvbivKDLFvyIa+99R6vvfUu675ewzNPPg7Ab//wOA/98tekpaezY385y1b036Lu3jt/zJHqw7y7ZBkv/+sNXn35JU7U1fXa5vab53P8eB1v/Pd9Plu1lknFU7j2iu/QUN8drB5sLAC/+82v+OuTf+anDz3C6vVbee7lV0lMSgLAZrNxw9VXEBkZxQeffM6Hn31JREQE37vmCqxW9UsBBALh4hjlNLRZh1XO8kSrleZOG9FhoeV3VMP/7A1K+eRlWWY4BnR5eRmyLFNQWDjgNl+tXMGeb3exYcce0jNcvTyfefFlzpk5ja2bNzG1ZDrgskr/8vzfiIyKAuC719/IV6tW8AiPEm02ExkZiUarJSk5pd/zlJUe4Mvln/LJl195jvnEX1/grNOmeLZZv+5rtm7ZxK7SQ542dr/53R9YtvQjlnywmJtuvW3IsbS2tPDyC8/y+8ee5PrvfR+AnLw8Zs46A4AP3vsvstPJE3993tP95qnnXqIoK4W1X63m3PPn+v5BhxhCoEc5w3FvuKk43sbkjBjlBqMAgSiS1B8ugVbgOMPNffZivwP795KWnuERZ4CiseMwm2M4sH+fR0wzs7I9ggiuZrDHT7J+Bz3Pvr3odDqKp07zvDamsAizOcbzfPeunbS1tjIuN73Xvp0dHb3aaw02lv3792KxWDjznP6Dn7t37qD8YBn56Ym9z9HZOWpaeAmBHuX4GiDsScWJ9hAU6OBZ0ME8Tm5eAZIkUbp/v99jODkbQ41msG2trSSnpPDekr5dWqJjzF6NxRQ2eFu8trY2Jk+ZynN/e7XPe/EJCX13GIEIH/Qop7px+BZ0VX17SJXDtDmcNLQHx4J2yihSl2O4n2dsXBznnn8Br7z8Am1tfeMDTY2NjCkcy5Hqw1QfrvK8vm/vHpqaGiksGjvsMZ9MQWERdrud7Vu3eF4rPbCfpqZGz/PJxVOpralBq9ORm5/f6xEf75145uYXYDKZWLNqRb/vTy6eQnlZGQmJiX3OEW0297vPSEMI9Cimsd1KS+fw2ylZ7U6O+CHwSnOi1erNTF81/L1ZOWXZr/H/4fGncDicXHLeWSz5YDEHy0rZv28vL7/wLJddcC5nzzmPceMnsuD2W9mxbStbNm/k7jt+xKwzz2LKtBK/xt6TgjGFzJl7IQ/ddzdbNm1g+9YtPHD3nb0aQZ895zymz5jJrfOvY+UXn3OospKN69ex6H8Xsm3LZq/OExYWxoL7HuD/fv1L3v7Pv6k4eJDNG9fzxr9eBeDq624gLj6em793Ld+sXUNlRQVff7WaXz50P0eqDyt2vcFECPQoxh//s5vy46GTzREs94abYfuP3fv7KfDZubksX72WM846m9/8z8Oce3oJ18+7lK9WreQPTzyNJEm8+p+3McfEMu87F3DdlZeSnZPDS6+85td5++Mvz71IcmoqV33nQn540w18/5bbiE/s9gVLksS/33mf02efwX0LfswZJZP4yQ9/wOGqQ54sDG+4/6FH+Mld9/Kn3/8vZ82Ywo9vvcnjow4PD+f9T5aTnpHJD79/I2fPmML9d/0ES6eFqKjR0cNUNI0dxSzbdZQ9R/2rq5EQaeCmWTnKDMhPVuyrZduhRr+O4W3T2P7QSJAQafRkDPhKQ7tV1OAYZYimsYJho4QFfbzVSkun+qU9vSHY/ROd8vDLj8qyjE2Is8BHhECPUprabX75n3tScXz4mSBKciJIKXY9Ga4FbHU4FVzqIjhVEAI9SqnyI73uZEJhVWFzpy0k6lRbhivQwnoWDAMh0KMUJdwbbg6FQLpdsFYQnoxtmIX2hUALhoNYqAKsLT3OvpoWJECjkZAkCQmYlR9PfmJksIc3LPzJfz4Zd7pdZlzwqtspt4KwK8DnR2zcancSptd6vb3DKavaOksQeiiVe3HKC7TN4WT74aZ+p89bDzWOSIFuarfR3KFsYK/iRFtQBVqpAKFddgX7HHY7et+SODz4KtAd1uC7ZgSBpb3d5WL0t4a2qgK9evVqHnvsMTZv3szRo0dZvHgx8+bNU/OUPrO/pmVA3+bhhnYa2qzERhgCPCr/ONyofFCv4kQ7Z41R/LBeo1QOtE2WqLdAeMNx4nRaJMl3L5/TJmGQvBNdu8NJo8I3S0Ho0Kl19kqzk2WZ9vZ2amtriYmJQav1/kbeH6oKdFtbG8XFxfzwhz/k6quvVvNUw2ZXdRNV+3fy0d8e4/Lbf0Zm4aRez3dmxXJ2YeLQBwohlPQ/uznRaqHT5vDJclQKm8NJY7tSIidR0W4gUmeho7OK4danqzfq0HhR3a7dasfuZasxwcijaYDvQUxMDCkp/VcD9AVVBfqSSy7hkksuUfMUflHXYuFIYycbl39A6fb1bPr8AzILJ/V6XjRxCmcUJIRs+6f+UEOgZRmONnWSm9C3FrHa1Cvcxdsqa9jaFIZRIyMNM/nttNw4JqQNXu+hqr6NjXtqh3V8wcjgmmkZRJl6uzH0er3flrObkPJBWywWLJbuqWxzc7Nq56qsrGTZpv0crmll26qlAGz6/EOyiiaz+YsPAdi68mNOu+Aq3u+sZvrYbLKzh9dHLpA0dSjvf3ZzpLEjKAKtxgIVGYlO5/BvuusqmkmKiRrQL293OFlTfpQOZ+BnHILAYTCGEaZizfSQEuhFixbx6KOPBuRcOTk5fV7raG3i33/8med5a+MJnlhwNU90PQ/hVfEehtM9xetjB6lwUqB7EHqDzSHzwbZqrpyS3q9Ib6psUNAtIzhVCSmBfuSRR7j//vs9z5ubm8nMzBxkj+Gz6OmX+OVP78TpmAC8DrjXxPe0qg4DV6DRNvL8S39XZRxKc3L9Z6cT3vlLMocPGF1XJoEkuf7V651cfPMJCoq9E97a5s6gdPseKge6vkbHjq8ikWUJSZLRaF3XKGlkcsZ3klGgjsDbHDIfbj/ClVPSyIjtFummDhubKuoH2VMg8I6QEmij0ehpj6M22TMv4r6n3+GJBQ5g0gBbZQF3cN/T0xl75pkBGZe/1DR39nq+b3M46z8Z2Ff6wYsa7n/2EN7U/7E5ZGpbOkk1D15IXWmGyoF+7fepVO7pf0xGk5Nf//sgpkh1FopY7U4+2HaEq6amkxbjGsOq/XXYFAwM1tfoeH1RKh2tGjRa0GhcNyGNViYqxsG199UQFStS+UYjISXQgaK2uZNjTZ20NhmBC7tenYfLYnYzB3gMuB2ncyvfHm1mdn58v5WrQgWbw0l9W+9p9bolMQBMndPM9PObkZFABrtN4vVFKVSXhlG130hWkXdW5pHGjoAKdMsQS7wPHzBSuceEVicz5ZwWZBlkJ8iyxMGdJprrdWxcHs3ZVzWqNkar3cnirdVcPS2dTpuTstpWRY//5VtxVHw78GeelGXlstuOK3rOQOJ0wiu/SaN0ezganYxOJ6PVyWj1MnqDzPk31FNynn9VGUcqqgp0a2srpaWlnufl5eVs27aNuLg4srKy1Dz1oOysbgLgwNZxgA5D2Fau+PE0vv5oH8cq9pGaW8jMi2P54IUmnM4cag61kFXkoLSulbEpoVv2tK7F0ivbobFOx7frXUG9C+fXk5zV2xItPruVzV9Es3ZJDFlFNV6do7qxk5IAxkqHChCuXeKaHUw+s4X5Pz/W6701H5p576/JrFtq5qx5jV7NEoaL1e7kvS3ViqchdrZp2PyF6zv33XtqSEi34nRIOJ0SR8qMfPxKAus/MXPRTSfQG0I/RtIfB7aEs/ubgReEffBCIsVntqIbodfnD6oK9KZNm5gzp7vho9u/fPPNN/Pqq6+qeeoBsdqdfLpqLe89/xh1h5cDcM3dKZx2wQ3MuvR6LO1tGMMjkCSJE8ccrH4P1i+LZ+NnV2O5/394+AeXBWXc3nCye2P9smicTom8Se19xBlg9mWNbP4imm0ro7jyJ3WYIoZ2AxwNcKBwMPdGR5uGLV+6xGv2ZU193p9+fgtLXk6k5pCRsp0mCiarO3ar3al4zY2Nn0dj6dCQlGlh1qVNvW4yRSVtrF1iprFOz/bVkUyfOzKtzNXvxwBw+iWNnHNNAw67hMMuYbdJvPa7VBqP69k2gq/PH1Sdr5977rnIstznESxxBth3rIV1yxZTtsNMc30EpigHxWe7pqSSJBEWEekpyD7rO64f/cGdKZRuP8zHi9+mvi34JS8HoraHtelw4PE9z7q0r3gB5IzvJCXbgtWiYfMXUf1uczLtVgcNAfwMBltBuOnzaKwWDSnZFvIm9RXfsAgn0+a4UjXXLRl5PepkGdZ+5Br3GVc09ZkBaLXdN6Y1H8YEeHTKUHtYz54NkUiSzJzrGkjOspGWZyWz0ELuhE5mX+66vq/ejw1qu7NgEboOVYWprKxk8+bNfPDFmq685zsAGDe9nJrKndTXVPfavr6mGmvnVjLG1OH6mG5j68qPefezr9i8eTOVlZUBv4ahqO1hQe/ZEEHjcT0RZjvFZ/bvE5WkbvFeuyTG6x9AINPtBhJoWe52b8y6rK94uXH/wHesiaKlYWTlJJduM1FzyIjR5OS0uf2vCZh5cRNancyhvSaq9gcmwK4kaz6IAWDcjDYS0/umJZ5+SRM6vZOq/WFU7g3r8/5o55QR6JycHKZPn86vbrmM1kYT4HJVbFlxOU/edQ2/vem8Xtv/9qbzePKuazh84J6uV35Ea2MjP756LtOnT+83jzqYnBwgXLfUJV6nXdA8qO9u+txm9EYnxyqMVHzr3Q/gaFPn0BspgN3hpKGt/1zi8l0maiqNGIxOpg8gXgAZBRayijpw2CU2fBq68YP+WPNRDAAl5zcTNoD7KSrWQfFZrqn/113bjxQ62jRs/Mz1PR0oiBsZ42BaV4Dwqy5XyKnEKSPQr7/+Olqd2+V+Gy73+ypgLxqtlvk/f6zX9vN//hgarRZ4D6gDMoHvAKDT6Xj99dcDNXSv6BkgrD+mY+9GV3BwIPeGG1Okk6nnun4A65bGeHWuQHX6PjHIEu+vu6znqee1DOk7d1vR65bG4Bwh2WgNtTp2rXUFzs68onHQbc/oen/LiijamkfOT3rDpy7/enK2hTFTBy7wdeaVjQBsXx1F04mRNQvyl5Hz1/ST+fPn8+fXlwBa4Eddr77o+u9/P+WJX97FTbOyuXl2DreekcOzC+/lqX8vBazAK13bu9wir33wOfPnzw/sBQxBT//zN5+YkWWJMVP7nzaejNvXvm1VpFc/8Po2a0BKaA6UwdHSoGXHVy6f+RmXNQ55nClnt2CKdFBfo2ff5uCVTPWFdUvNyE6J/MntpOQM7vPPGd9Jen4ndmu3RRrqOB3d7o2zrhw8wyajwELexHacDom1XWmjpwqnjEDLstzlo70ElzVch8s6hnGpUSRFhZEQaSQuwkBMuAGzSU9hijtw9reufy8BsqhpCcwU3xfcGRwOO6z/1PUjnT2E9ewma2zXD9ymYdPn3rkBAuGHHsj/vOHTaBx2iayiDjLGDJ2/bQiTOe0ClxtkJPzA7VaJbz5xBwcbh9xekrpnCV8vMeMcAc1b9myI4MRRA6YoByXnD11z58x5jYDrxmW3jpzCZf5yygj08VYruohYtPq7ASicdoDMwkJi4xPJTk/rd58J+VmY4xLJLDSRlHkQ0GAIu4dObeALBg2F24LetS6SlnodUbF2Js72bsFEz2DhOi+DhUebAiHQfS1Hp6PbFdNfat1AzLq0EYBvN0TQUBva67O2r4mktVFHdLydSV7+DafNaSYswsGJI4YRMUvwpNZd3ITRNPQXbtLsVmISbLQ26ti2euQ10Rgup4xAVzW0I8sZOB0XAHDN3anc98w7bNi5j4yMjH73ycjI4PMNO7nvmXe4+GZXAM0Yfg/aiKSQSrezOZzUd4mZOzjoiu57f4xp5zVjNDmpPWygbOfQKwUD4Yfuz4Letzmc+ho9pkgHU87xPi82OctGQXE7slPim49D2w3wdVfK3KxLG73+GxpNMjMubO61f6hyrMLAga0RSBrZqxkCgFYHsy93bbt68amTcnfqCHR9u8s365QYM6WdxHQbcREGClJjB91vUnYCRr2WibNaiYyx01KvZ/c3kX2KEgWT462uAGFdtZ79WyKQJJnTL/HeugQIC5eZ6s4ZXjq0gNU0W7A71JtLH2+19OvndrsoTrugGUOYb7/S2V3+6m+WmXHY/R2hOlSXGan41oRGKzPLx7+hW8D2bIig/ljozhLc1vOk2a3EJXv/hzj9O66Uu8MHwqjcc2qk3J0SAu10yhw60dFj4UYjAJMyhhYio05LUXIUOr3LKgWXgKlZ1tNXappdlqbbb1k0vZ24FN8VyO2z3vFVFK2Ng0fLHU6ZGhU7bVfV970BNtTq+HaDy70024vg4MlMnN1KVKydlvruDIlQw73gZPKZLUTH+xaITcqwUTitDVmWWOtlRk6gaWvuXrp+Vpdf2Vsizc5TLuXulBDo2hYLZbsNNNfrMEU6mDi7FZ1GYnyqd1PdyZmu7WZc5LIwD2wLp+xI6AQK3QtUdneJzsyLfLO83GSMsZBe0InDLrFn49B+TDXdHIf6EegdayI9mQ1Jmb7XWtbpYUbXZ7NlpXcrJwOJzSqx9UvXuM68Ynh/Q3dK3vpl0SGZUrh+mRmbRUN6fme/qz+Hwi3q27+KovF46M4SlOKUEOjDDe2U7nD5VfMnd6DTQ0FSJCaDdzmVSVFhpJjDSEy3EZtkw+mQ2LVVH9Alz4NR02Kh+YSW2sMGJEkeNKd0KIpKXPuW7QieQDudcr9tu8q2u8Y0bkbbsI894fQ2z7FCLduh4tswrBYN0XF2cicO77MdN7ONsAgHbU06qstCb2Whe+Zycl0Rb0nPt5A7oQOnQ2LHV6E5C1KSU0KgV3y9jhVvHwAgf7JLgLxxb/RkUrq51/5vP/k5n65cq+Aoh4e9K0DoDuyl5VsIjxq+8rivr3S7N4HCTlW6zNS0dPYpOuR04LnGguLh34AyCzsxhDlpb9FyrCK0urWXbnPdgAqK24ddeU+rxVMU6sC20Mrm6GjTcKhrubZfN9lZrsyW/VtC6/rUYNQLtMMps+Sdd2hvcRXlLyjuID7S0KsDhjcUpURh1Gs83Ucaasfw5n+Cv5qwritA6LZ484cxbexJ7oQONBqZ+mMG6msGn0J22hyqZLMcOtFXgI+UG+lo1WIMd5DuR4cUrQ6PdVq6PbR+4O7xFEzxLwDt3v/A1tC6vrIdJpxOicR0K7FJw4/SFk7rnuWFarBXKUatQLuLIy3/ah0bv6gGokBqwGnfgra+3OdiR0cOV6E9UU5E9PauV07ji6WfsWXLlqAWT6rtChC6LV5/rEtwZXNkFLp82mU7vLOilaY//3Ppti4X1aQO/G2Y7P6MSkPIwrR0SJ5iQGP8FGj3/gd3mbCHUFtE9w3DHxccQFqehQizHUuHZtQXUBq1XvbexYy6GsHKq3jy7qs9r/oyPe99vEogi9bGQkpKSoZ1PKWoae6kuV5LbZURSZKH7bvsSUFxB4f2mijbEc5pFwyea1xZ3+azu2gwbA4nx/opxuSZIShQ09k9CyrbacLpAE0IlHc4uMuE0yERl2wjPtU/szAlx0pkjJ3WRh2H9pqGFYxTg/1dAu22gIeLRuO6CW1bFc3+LRHkTQydgL3SjFoL+vXXX0fnKY50Tte/K4HhFTvqfbyVXf+eO+zjKUVti8Vj6ablWYiI9j/y5YsfuvJEu6L50EcaO7A7e9/olPI/u8kY04kx3EFHq5YjB0MjkObxP/tpPYNrZajbit4fIm6OphNaaipdRoT7++UPbpEf7X7oUSvQ8+fP5+u163AVRzqr69VVAKxfv97nYkfz589n/fr1Xc9Wdv17zrCPpwR2h5MTrVZFrUvo7Ycealm01e7sN+NiuFTV9z1WdZmRzjYtYREO0vP9z73Wart99aHihz7QI0CoBAVTuq5vW2Ab/A6E272RMUYZI8It0If2htHRNmplbPQKNLiroU0FooEGYKdCR17d9e9MIHg/gOOt1q4AoXLWJfjuhz54XLkmqf36n7d3B0CVcke43RzezBLUpr1F40mJcwurvxR2+Xkr95qwdAS/uJBS/mc3ccl2EtOtOJ2SV9/RkcqoFmibPgqj6VIAUnNrGTNhIikpKSQlJQ3reElJSaSkpJA7NgJTZBNgwBR5MWFRgy8XV4ua5k5aGrTUHBrc/5weYyLcy5xvN90CNrSFebBu+ClTPem0Oajtp1JgmSeHXbnl9e6bWdlOE44gL+g4uNOE7JRIyrASk6BMWkJcio3YZBsOu0T57uAKmCz39D8r812BbrEfzW6OUS3QlrAYcic8AMCMC2N5d9lKKioqBiyONBQZGRlUVFTwrw8/Z/xMl1Uy+9IXkCPiFBuzL/T0P6fm9j91zEuM4Kpp6ZxRkODTsd1i6I110tJp79Vua7gcbmjvUwTH0dP/rJB1CS5/vSnSgaVdS3VpcP3QBxT0P7vp6YcOdj50bZWepuN6dHonOeOVC+h5/NCbQ6+6pFKMWoG22p1U11uo+LZrelzcQUZsOEajfz9Go9FIZly4R8DKd0cq6oP1hZrmzu7pf3HfMYxLjebyyWnotRompEWTFuN9SpLbD33i6NB+aIAyBazo/twb1aVGLO1aTJEO0nKVq/2h0eLJbgi2H1oNgYYeAh3kQKH7+nIndmAwKpfpNGZKO5JGpvawd9/RkcioFejqxg4OlxrobHcFl8aMtRMboczKsVSzicIpLkugcl8YB48GPs3H7nBS39a9grDgpOn/1KwYLpqQjEbjsvQlSWLO2CQ0Xi5RC4Yfur8AoRr+ZzehkA/d0qDlWEWX/1mhIK8bt+BXlxppbwneT33/FpeFW6iQ/9mNKdJJZtd3NNg3IbUYtQJdVd/e68edGa/cH9Cg0zBhrAZzgg2HTcPOLVqaOwO7IuB4q5XGExpP6lLPXNfZ+fGcW5SEdJIYJ0WFUZzpfc6yL37o2mYLLX58Bi2dtn5XJXbPEJQv7+oW6IO7TEFbkeZxUeVZiIxR1hlujneQnGVBlqWgzRIcju5ArL/5z/1RNMrT7UatQDd22HoFl9JjlQ2UuNwc7gUP4Rzux/pTk2PNnR7ruaf/+cwxCczMix9wv1n58UQavZsO+uKHBv+Chf25Nxx2V/duUC5DpSepuVYioh1YOzVU7Q/OijT39H+MCtcHPdLtgpStcnh/GJ1tLheVEimSJzPGLdBbQ6/4lRKMWoF2OFzRcXBZgukxyn5BM2JN3VPk7aaAF/CvbujwVHdz3yj0WonijJhB9zPqtJxV6F3A0Fc/tD9ujv7cG4cPhGHp0BAe5SA1V/maHxoN5HkW5QTHAnO7V/xd3j0QwQ4U9vSvq7FiM2dcBwajk9ZGXcgVv1KCUSvQ5fv0Hv9zXpGdhEhl/3hpMSbGFHf5ofeGcfBYYP3Q1Y29S6gC5CREYNAN/ScdmxJNZtzQP1hf/dBV9R1Y7MObpvdXoN9t9eVPbkej0jfV48YJwoKOxjodddUGJI1MnsL+Zzf5k9uRJJmaSiPNJwK/pt3telDa/+xGp8fz2bl93aOJUSvQuza5BDlvYgcZ8WF9/LH+otdqmDhOQ3S8HYdNw44A+qEb2qwcOwY1la7gktsVUZjsfRH6OUWJaDVDfya++KEdTrnfSnRDUd9mpdXS1wlcqvAKyf5wuxbKvzUFvFu0+waUMcaCKUKd+XlEtNPjWjgQ4FmCtVOi/FuX60gN/7Oboq7c6tHohw6IQD/77LPk5OQQFhbGzJkz2bBhg+rn3LWpKzJe3EGGwv5nN5lx4Z7sibIdgfNDVzd2eNw3qV31N/RaidwE7y2I+EgjxZkxQ27nqx96OOl2Q/qfVZr+AyRnW4k027FZNBzaF1g/9AGV3Rtu3J9foLNVyneZcNg0xCTaSEhTz3hx+6HLdgb+Jqs2qgv0W2+9xf3338/ChQvZsmULxcXFXHTRRdTW1qp2zvXrN7H5a5dFlj+5nfQYdb6YGbEmT/7x6sWVfPlVYAr4H27oYNtq13mTs6oAl3tDr/Xtzzk2ZWiL21c/dMWJNpxO73NdZVlmz9HmPq9X7Q/D2qkhItpBSrZ6nWskqTuHPJCBNFnu4X9WKUDoxr3iLtCpaPu3dVevU3gC24vUHCtRsa6bbMUoayarukA/8cQT3H777dx6662MHz+eF154gfDwcP7xj3+ods4nn/wSpyMSra6d3CIbSVHqrBRLizFR2BUl72wbz7tvvavKeU6murGDA1td12RpXwr45t5wkxRlHLLtl69+6A6rgyNN3s8kdh9p7re8qCe9TkX/s5vuYG/gBOzEUT0NtXq0OpmcCerOvPImdqDRytTX6DlxNHALOg5sUbb+xkBIUnCWfQeiurCqX32r1crmzZuZO3du9wk1GubOncu6desUPZe7QP+WLVtYutSdTfAVlpoDbN26RfGC+pWVlezYtpUIaTeSdAwIY/3yGlav26BaAf/KykpWrV3Ppq+/pb3FtVz90N5nOVb2LQ2Ve30+pyRJZHrRWcYXPzR4n27XaXOwpvR4v+95AoT9rJBUGrdAV3wbhi1AU2S3eyNrbAdGk7q/dKNJJntsZ6/zqk1rU3cBKLUChD0JRj70S79I5/JLtGzZot45VL2dHj9+HIfDQXJycq/Xk5OT2bt3b5/tLRYLFkt3rmRzc9+p70D0Lqj/AQAO+3J+9v0/e15VsqB+7/P9G/geHa3TOGf2TFXO1/uc1wCXAjtoaz7An+68ij8N85zZ8eHsrxm8KH/+5Ha+fCvOaxfAwbpWzi5MHHK7r0uP02Htm/Vht0HF7v5XSKpBUqaN6Dg7zfU6KveEeW5IatLdAScwcYsxU9sp323iwNZwTr/E+9/VcCndHo4sS6TkWIiKVb8aldsPXXUgjPYWjV99Ob2hs921+GefXcKsXL+KPoRUFseiRYswm82eR2Zmptf79i6oXwVU40+Bft/Ot7Lr33NUO5/7nFqtju4GBKs87w33nN6k2/lSHxqgod1G+fHBreia5k52Vjf1+17Vfld36wiznZQc9Tunu/zQgXNz9PI/qxwgdOMJFG4PD8jUXOnyokMRk2B3rZp0SgGZJbj6IUpk58jk56t3HlUFOiEhAa1WS01NTa/Xa2pqSElJ6bP9I488QlNTk+dRVVXl9bl6F9S/C8gANgPqFNTvfT63UM4CjKoV8J8/fz6/+tti3J1cYIXnveGe02zSExuuH3QbX/3QAB9tPzKgZS7LMiv21g4oFJ7mqZM7VA0u9SSQ+dC1VQZaGnToDE6yxwUmfz5nbCd6o5OWhsAs6HALdCDcG24C2WXFfY7zzlfXUldVoA0GAyUlJXzxxRee15xOJ1988QWzZs3qs73RaCQ6OrrXYzhouqJKSuc+D8wB4CgQBszoN6dXCdqtdo7VSsCkrldWK3KNWV7UKfHVD+1wyny88yjbqxr7vLf7SDNH+wkMulGqAa4vuM9VudeEtVPd782BrptA7vhO9IbA9LHUGWTyuuqFq21h1tfoOH7EgEYj+91l3hfcAr1vs/oC7T7HnPPV/fup7uK4//77+dvf/sY///lP9uzZw5133klbWxu33nqr4udyF9QvKSnhzv/5A0UTi/0q0O/t+cZOKiaz8AQAhrDvYNf7nlHhDdUNHZw4mguA3nCA795zDxOLp/p9jVlxQ+dP+9Kn0I0sw5d7a1lXdsLz2mCBQQC7VaLi28AFCN0kpNmISQhMgXsl+w/6gifdTmWBdlvPWWM7CVNpAU5/FBS3o9W5XHHHjww+K/SHhlqdq0mzRubsc9QVaNVzbq6//nrq6ur49a9/zbFjx5gyZQrLli3rEzhUAndBfYPBwEc7jvLQvQtIjdL7XQN6qPMda7Hz4P+1UrUfMgvvwm5SrgVUT6obOzhakQrAjIsTOOfKG/n3n3+J7LD7dY0ZsSY0koRzEOdkTz90fY2OuGTvZwnfHDxBh83OnKIk1pb1Hxh0c2i/EZtFQ2SMneQs9f3PbiTJJZibPjdTut1EUYk64ul09nDhBHCGAN3+7rLtri4yWpVWfnv8zwG+ARlNMtnjOji4M5z9W8JJSOs/xuEvbus5q6iTmBh1JTQgQcK77rqLyspKLBYL69evZ+bMmUPvNEyMRiOSJKGRIC3G/wL93pwvLbY7H7pybzjlNepYftWNHT1+3B3kJkRi0Gn9vsYwvZYU8+DHGI4fuifbq5p4b0s1Ow4P/qPxWJfFgfM/u+mu/KaehXm03EB7ixZDmJOsosDWb0nPd3WR6WzXclil6n2y3L2kfIyKy7sHwp1ut09FP7Tb/6zWTbwnIZXFoSRJUWFeFQ5SAr1Ww+QJWqJi7ditGnZt0ylel8Nid1B+2OYp7p4/qYPC5EjFju9NNofbD+3uIu4rh+r7trQ6me4O5YH/cbst2qp9YXS2q3N3cN+A8iZ2oA1wExCNtvsa1XJz1FQaaKnXoTc6yRkb+EYWbj906bZwnCpk9zkd3UWZikqU6684EKNWoJWu/zwUmXGmXn7aaoXbYB1p7PSIV0qOhbgEmRwfam8MRXb80MfqWV5VDezW7uI6gcoP7klcsp24FFenaHcdEKXxzIACPP130+2HVuf63M1h8yZ2oAtQALQnmYWdmCIddLRqVanxXV1mpL1FS1i4IyAzoFEr0CnRgV2TnxHbo4D/jnDF+xRWN3T0KL/ZQXZ8uM+1NwYjNXroGcfJfmilqdwbht2qISrWTlJm4PzPPXHfGNSwMB2ObvdQoP2zbsZ0uXHKd5tUWTUZ6Pznk9Fou29+aqTbuf3PBcWBmQGNWoH2ppSmkqTGhHn80OXfhlFeo+wXtLqx3WNBF0xuJ8sLl4QvaDTSkFX/jCaZzC6roUwFP21pjw44gfY/u1GzLkd1qdFTo1yN7iLekJRpJTrO5Yqr2K2sERMKNyBQ1w+9L4DuDRjFAh1o9FoNkye6sg/sVg07t+r86tHXE5vDSVmVjaPlLv9z3qQOr2po+Io3bg73LKF0GIHCoSjrEQANFu5AYXWZkY5WZX8epT064KjRXcQbehYWUnqWULU/jM52LeFRwbsBQbcfunKPCUuHcnd6S4fkKUFQGIAAIQiBVpTMuPAegTSTYm6OY02dHOhyb6RkW0hLlRTrUN4Tb6xyt4WptAVts0pUePzPA3/5S7JjKcmOVW2GFJNgJzHDiuyUPD0flSJQ9Z+HwlP5TeHyo+7qdQXF6rS38paENBtxKVYcdmnYAe3+cC/vjkuxqlrfuidCoBUkI9ZE/iTXl3/XN05unPcdNm3a5PdxXfWfXcdNzj6sWgOCuAgDUWGDO9Y8fugaPfXHlHPCVe4Nw27TEBVnJzGj/y+/RpKYlh3L2YWJ3HR6tk8NCnzB4+ZQ0MK029RtgOsL7htE1f4wOtqUkwB3/edg+Z97UqiCm2Nfj/S6QLnghEArSKrZROFUl4+2utTM9g0bee211/w+bnVjB/u3uCxmS8fHZKjg3nAzlBXdyw+toHXiyX8exP+cnxTh6UgeG2Fg3tR05k1NH7KWiK+oka3iaUAQoAJQgxGbZCchzTVLOKiQq8ra2T0DCgWBducoKxko3L+5uwFBoBACrSBHq6uI1O7CFGnB6TACp/Gf//yHLVu2DKtGdGVlJd+s38jqFdtpb84C4NC+ZzlRuU+1mtPBSrdzB5cG8z/317E8NyGCm2blkJeonDXtHsORg2G0NinzE+m+AXWo3oDAGzxuDoVmCeW7u9pbJdhITA/M9H8wxkxpR9K4muU2Hvd/ptdQq6PmkGt5dyBdVAFOlR/ddNdrfgu4DjiHurq1lJSUeLbxpV5z9/Hm4ar/vJv25n2ce4Z6Nacz40xI0uDdIgomd/DFm91NXf3FZpWo7GpVlD/A9D8uwjDgYhqtRmJ8arTXjQKGIirWQXK2hZpKI2U7wik+y/+l+weCVH9jIMZMbWfd0hjF2mC5/dljVG5v5S3hUU4yx3RyaJ+JA1vCOe1C/2pguy3xrKJO1WtN9yQE7uWjh+4a0e7yo+cCLqUbTr3m119/HY1WR3d50ZWe99SqOR1u0JEQOfiy75wJrhZKDQr5oSu+dfmfo+PsA1pfkzIGr4qeHR+BTsHAoa/V+wbD2wBoIHGP41iFkZYG/yN6waq/MRhKVrcL5PLungiBVpD58+ez5ut1dAvpGYDLPzqces1zL7+G+55+m+4C/Ss976lVcxogPWbofGj3KiolrGhPfndx/9aXXuuykAfDoNN4VTbVW5R04/S8ASVlBn/6DxBpdpKW1/U39NPN0dasobrUdVMPBf+zG3cq3P6t4Tj9MHqdzu7l3YXTApP/7EYItMLodRrgW6AOCAdmDPtYOw830dGmByZ3vbLaU+taTVLMQy9gcC9rL1NAwIbqP1iUEk2YfmgrLz9Rudok7uurqfTfwuxZvS4Upv9ulPJDl+1wtbdKzrJgjle/vZW35IzrwBDmpLVRx9Hy4aelVpcaaWvWYgx3eHo7Bgoh0AqTlJREXEISpsitAJjjryMqNgFjVIxPx7HYHeyraemq/6xBZyjjlofupaSkRNUa1wBp5qFFt8CzYMW/H7fVIlG5d/Dpf/EQ7g03uQkRigmgkhampwFBCE3/AcZM7VrWvtW/m+z+AHXv9hWdvvtG67aAh8Peja59xwRoeXdPhEArTEZGBl9v28PFPygGIDn7Nn712grq8a2I/75jLVjtTo5WpAEw48I4fnDrj1i/fj0VFRVkZGQoPnY35nA94YbBrcaefugTR4f/ra38NgyHTUN0vL3f5P9UcxhJXtZViTDqSPXC+vcWT10OP2YJlg6JQ3vdy5+Dt0KyP/IntaPRumqrnDg6vFRFWe7ODw5keytvGTvdNaZNn0cPqxejwwHrPnYZCJPOUKfO+2AIgVaB3GQzhVNc1lfFbhNIBnYfacbu8N4R5q6bXNaj+7Mrw0JSvcY1QKovfmg/AmlDTf+HCg6ejJJujp5Fd4brw9y/xbX6LDbJRlxKaPif3RhNMrldbbA2Lh9ee7kD20ycOGLAaHKG3AwBoOT8ZgxGJ0fLjcOq4Lfz60ga6/REmu1MObf/HptqIgRaBYw6LZMmSUSa7VgtGvZviaDD6uBArXd34KNNHdS1WDhxVM+Rg93lN9VcoHIy3lii7int9tXDa/HldMKWla59i/pJ/g/TaylK9u3YBUnKCfSYKe0YTU7qjxmG7eZYtzQGgCnntISU/9nN7EsbAZeVaB/G/WPNB7EAnHZhE2HhgS8vOhThUU5Ou8iVYrf6vVif9/9qsWuf2Zc3Bax/ZE+EQKtEZpyJaee57rjuKdLOIbqJuHFv596vqKSNvCydV4EypfBGoGd05Zbu2xw+rCnyga3hnDhiICzcweSz+lonE9Ki0flYUjUm3EB8pDJ1SowmmelzXdf49Ue+WfIAx4/o2bspAkmSmXWpOu2X/GXyma1Ex9tpqdex/SvfboYNtTp2f+Pyz55xeaMKo1OGs69qQJJkvl0fSW2V99/Tqv1Gyneb0OpkZl/WqN4AB0EItEpkxJo8P8pv10fQWKejurGD462DV/nqtDnYX9OC3QYbPnVNO2dd2qRK9brBSI4OQzOEyRefaqNoehuyLHluJr6wbqlrn5K5LRhNfa2TyT66N9wo6eZw/zB3rYuksc43X7v7+opK2gNWXMdXtLpuK3rNBzE+7bt2iRnZKTFmahvJWaF5fQCJ6TbGz3Slx61e7L0V7d526rktRMcFJztFCLRKpMWYSM22kT+5Hdkp8c0n3lnRe4+1YHPI7FobSWujjug4OxNOb/WqJZWS6LUaEqOG9nW7f9wblkVj96EAfNMJLbvWuoR0VtcxehIbricmfHiWsJICnZprJW+S62/oy03IZpXY8Klr+9khbF2CywDQ6mQq95io2u9dfMNm7f5On3lFo4qjU4ZzrmkAXL72tuahZa/5hJZtq1wzirPmNag6tsEQAq0SYXotiVFGjxW9domJZx+8lU9Wfo3VPnDEaefhRgBWvOOy1sbNOIBeLw25eEQNUmOGdnOMP72N6Hg7rU06dn7tvTBuWGbG6ZTIGd9BWm7f4kHpfswYkqONQ1bl8wX39P2bT7z3025fHUlbs5aYRBvjZwR2cYOvRMU6KD7b5WL66v0Yr/bZtiqStiYdsUk2xp8e2tcHrhrc6fmd2Cwaz8xmML5eEoPDLpE7sZ3MwuDVthYCrSLpsSYmn9FKhNlOa2M4ZTvi+ebTxeyv6T8afKSxg+OtVuqq9VTtTwKcOJ0vkRxtDFgD3J5444fWauH0S7puQl588cHVeNNtfQ3k20vz4uYwEJIkKVo8adIZrUTFuvy0bqt/KNYuiQFc1mkwayN7y1lXNgKwdVWUVwtz1nzYFTy7rBHtCLg+SYKzr3ZZwms+jB30RmuzSqxb4vp+nj2vMQCjGxgh0CoitdZxtGInY6eXd71yB1tXfsw7y1bz4rvL+c8Xm1mxt5Yv99bw5d4a3lm5lar9u/jsdXvX9p+wZ8NrtB89oFr1usFI9WLBCrgEWtLIlO0Ip+bQ0G6JvZsiaKjVEx7lYPIAhYgyYvxz6Sjp5tDpu29CX38UM+T21WUGKr41odHKzLw4NIODJ5M9rpPMok4cNo3n5jkQlXvDqNoXhk7vZObF/hUhCiRTz20hKs5O8wndoJlHW1dE0do1O5gYhNznnohqdipyTsmErv8rAA4Al9DaGMH//ugKzzZPfLbP8//3X3UGYAQOd73yIq2NJ7ht3lzPNkpXrxsMs0lPhFFLm2XwAElMop0JM9vYtS6StUvNXHVn3aDbr+2yTk67oBmDse/1RIXpMPtZ4zkjNhyjXoPFpkzlsdO/08Tnb8ZRtiOcYxWGQWs6u63nyWe2Bi24NBzOurKBN/6UytolZs67rn7AVXNffxgDuAQvMmbkXJ9OD2de3sgn/0xg1XuxTDuvb+qjLMOqruDgGVcEf3YgLGgVef3119FqdUAp8Dmuj/tHAGi0Wub//LFe28//+WNI0neBBFwi/bHnPbWq1w2Ft1a0OxC2aXk01s6Bg4UNtTr2dC2d7S84CEMXa/IGrUYi14va1t4Sm2RnQpev9esuAe6PzjYNm79wZd8EKzVruEw5u5XIGDtNx/XsHMCV09KgZesq13tnjIDg4MnMuqwRncHJ4QNhng43PSnbYeLoQSMGo9MzawomQqBVZP78+by59IuuZy92/XsboOO+p9+h5Pwrem1fcv4VpOU/1fXsZaDbOlGzet1geLt0unBaO3EpVjpatWwbZPr4zSeu1KyC4vYBK7ulK9TSK1/BRSvQHSzc9HkUne3934Q2fxGFtVNDcpbF02B3pKAzyMz6jkuU1gwQLFy/LBqHTUNWUQdZRcELng2XSLPTk9u+8t1YLB0SrU0aGut01FXrWfGOy3qefkFzQOs+D4QQaJXpbiH1IVADpAGX9bttzSED1aUJuIT57wBIQW6/MdSSbzcaDT0yVvr3YTrssH6I4CAoY0EDZMeHD5nL7QtjpraTmG7F0q5ly5d9l0bLcrePevZlTSG5cnAoZl3WhEYrc3BXONVlveMJDke3++bMrqDiSOScqxsB2LU2kkeuHMOvry3gf+fnsejWXPZscN3Ug5la1xMh0CqTlppCbEISmYVFFJUcAECru4fm+jqe+9kPqNq/E4Cq/Tt54eF1AOj0n5FZGMuPHv490wNQvW4wkqOMXnfQnnFhM1qdzKG9Jg6X9s2n3f1NJM31OiJj7Eyc3X/wxWTQEqdQx3KjTkucQqsKwXUTct9Yvv4opk/xnfJdJo5VuqbH0y8YOcGznsQk2Jl8putv8+8/pvKPhWm89It0nn8og7/ck+WpS1F8dnCDZ/6QnGVl2nm9/z5anYwx3EGk2c45VzeEzMIb1YKEv/vd71i6dCnbtm3DYDDQ2Nio1qlCmoyMDHbuPcBbm49Sf8zA724Gh30O21YtoXT7ejZ9/gGZhZP4ZtnHNB13+aRv+sU4Js5+h4smpjD+9w9jtVoDUiCpP3RdC1aONQ1dBzcq1sHkM1vYujKadUvNXHtvba/33fmnMy5qRjdADDAtxlUQSilSosM43qLcVPy0C5v5+NUEjpYb+dfvUtHpZZwOcNgljlW4/kZTz2vBFBH86fFwOWteA9tWRXGswui5pp4Eqy6Fksz/+TGuuasWrU5Gp5dDNhVSNYG2Wq1ce+21zJo1i7///e9qnWZEkB4fTXpsIw21ZWSPDadybzJbVswHSli7xMChfQ6q9t8ExCFJVZgTtlJbLmEsMCBJ5qCJs5sUc5hXAg2uKfLWldF884mZbauicDpBdkrIMtgsGlddiu8MHHxRekFOSnQYu6qVC/aERzmZNqeF9cvM/aZqSZIc0nUpvCF3Qic/+r9qThzTo+sSMK1eRquTCQt3hkzbLn+QJDBFhv5NVDWBfvTRRwF49dVX1TrFiGJiupkbZpwHXAF8gNMxDZiGww6VewDGASDLz/PU3YsA+COBTasbiDSziW00erVt/qQOMos6qdoXRkdrX7Ok+OwW4lMHnj5mKBQgdJNsVv7mdvntdaRkW3A4JDQal/Wl1cpodJCYYSWjYOQFz07GXbtCEFxCKg/aYrFgsXR/uZubR6Yfrz8Kk6O4+ReP868//BzZeR2Q2s9WLYArlU6n04XMzc2bFlhuJAnu/nMVJ47pkDSu5xqNjKQBjRbM8fYB9zXoNCQO0bDWVxIijOi1EjaHcje68Cgn51zTqNjxBIKBCCmBXrRokcfyHm3otRpuvukm4tPzeGLB1UNuv379eqZNmxaAkQ2N2aQn0qij1TKwuPZEZ5CHFWRJiwlDo2BnbgCNRiIpKozqxpGV8iYQgI9ZHA8//DCSJA362Lt377AH88gjj9DU1OR5VFVVDftYocjE9O70sz6BsK7nSgbIlMSbwkn+4k0vxOGQrGAbLIEgkPhkQT/wwAPccsstg26Tl5c37MEYjcagB8TUJDHKSF5mGlGxCcQkpjLz4u/y9Uf/4VjFPlJzC5l92ffY+Ol/6WysDVpa3UCkmsM4UKNuapVSC1RORsk+hQJBIPFJoBMTE0lMTFRrLKcE50wby69eW4FWr0eSJGZdej2W9jaM4RFIksSC//cTZuUEP3PjZLxd8j1cdBqJFC+bw/pKskrHFQjURjUf9KFDh6ivr+fQoUM4HA62bdsGQEFBAZGRyi7BHUkUpUSxOiLMU8RHkiTCIro/j8kZMRiNyi2uUIrk6DAMOs2gtaz9Or45zOf2Vt5iNrm6lLdbR05hH4EAVFxJ+Otf/5qpU6eycOFCWltbmTp1KlOnTmXTpk1qnXJEoNdqGJfSfwfljFiTYqvolEarkRRPgetJhsoNCXzJRBEIQgXVLOhXX301ZNLEQo3pObH9TruTokPLrXEymXHhHKxTJz9WLf+zm+ToMNXGLhCoRUil2Z0qRIXpGZ/mX73jYJCtUl9EjSSp7uNWy78tEKiJKJYk8Jr4SGV7/blJjFK/pZdwcQhGIkKgBT6RpYIVrbZ7A1xNfGP87NIiEAQaIdACn8hWsEuJGzVEvz9EPrRgpCEEWuATWXHhihaiN+g0ZAbAgoaRkQ+t00hEm4SlL3AhBFrgEyaDlqQo5YQuKy5ctfznkxkJfuiCpEiumZZOuCFECxQLAooQaIHPKOmSyEtU3mUyEImR3neHCRYT083EhBu4amq66oFTQegjvgECn8mOV0agNZJEXkLgVpXqtBoSFC5nqiSx4Xoyu25+SdFhXFGchi7EbygCdRECLfCZtBiTItZdakwYpgBP5VNUKOCvFJMyejfbzYwL55JJqYo2vhWMLIRAC3xGq5EUaU2VH0D3hptQDRRqNRLjUvuWAChIiuT8caFV2VAQOIRAC4ZFlgJujkC6N9yE6orC/MRIwg39LwKamG5mZl5cgEckCAWEQAuGhb/LvuMiDMQGoTBUXIQBoz70vvaT0s2Dvj8hbfD3BaOT0PumCkYE/i77DmT2Rk8kSSJZwTRBJTCb9GTGDe4yMpv0Ij/6FEQItGDYZPphReclBq8meKjlQ09MN3vV6kwJv79gZCEEWjBshptuF27QkhZEkQylQKFGkpiQ1n998JNRsx63IDQRAi0YNsNd9p2TEBHU5rihVJMjLzGCCKN3riJhQZ96CIEWDJtwg47EKN/zivOD6N4AiDDqVCmbOhyGCg72JDbCQKSXYi4YHQiBFvhFboJvwT6dRlJsJaI/hIIfOtqk9/mzCERpVkHoIARa4Bcl2bE+WXVZ8eHoA1QcaTBCIR96fGq0z64e4eY4tQj+L0UwojHqtJxdmOj19sFYnNIfoWBB5yT4PpMQFvSphRBogd8UpUR5VeFOkiA3SPnPJ5MUFRbUGhcGnWZY+djxEYaA1y8RBA8h0AJFOG9s0pCV17Ljw0MmyGXQaYiLDPxKRjcp0WFohlGpTpKUqYMiGBkIgRYoQmyEgZLs2AHfz0kI57LJaQEc0dAE0w+d5ofICjfHqYMQaIFizMiNw9zPcuTC5CiuKE4PieBgT4KZD+3PopMMYUGfMoTWL0YwotFpNcwZ27s05sR0M9+ZlBKSnUyCtaJQq5H8ClImRhlDsuCTQHnEX1mgKLkJERQkuTI1SrJjuWB8clBXDQ5GfIQhKG2lkqKMfs0mJEkizSys6FMBIdACxTm3KJGzxiT4lH4XDDQaiaRhrIT0FyV8yKIux6mBEGiB4kSF6ZmeMzIKzAcjH1qJLAwRKDw1UE2gKyoquO2228jNzcVkMpGfn8/ChQuxWq1qnVIg8JlAZ3JIkn8ZHG6So8JE1+9TANWSUvfu3YvT6eTFF1+koKCAXbt2cfvtt9PW1sbjjz+u1mkFAp9IDrAFHR9pJEzv/0ITjUYiJTqMQ/XtCoxKEKqoJtAXX3wxF198sed5Xl4e+/bt4/nnnxcCLQgZosP0RBp1tFrsATmfkily6bEmIdCjnIAu62pqaiIubmDfpMViwWKxeJ43NzcHYliCU5xkcxitta0BOZeSvmMRKBz9BMyJVVpayjPPPMMdd9wx4DaLFi3CbDZ7HpmZmYEanuAUJpB+aCX8z25SosOGXF4vGNn4LNAPP/wwkiQN+ti7d2+vfaqrq7n44ou59tpruf322wc89iOPPEJTU5PnUVVV5fsVCQQ+EqgVhTHhekVrkei0moD70AWBxedvywMPPMAtt9wy6DZ5eXme/z9y5Ahz5sxh9uzZvPTSS4PuZzQaMRoDn5cqOLVJijYiSSDL6p5HjSJHydFhVDd0KH5cQWjgs0AnJiaSmOjdAoTq6mrmzJlDSUkJr7zyChqNSAsShB5GnZa4CAMnWtVNAVXSveEmGAttBIFDtSBhdXU15557LtnZ2Tz++OPU1dV53ktJSVHrtALBsEiODlNdoNUI6g2nJ6Rg5KCaQC9fvpzS0lJKS0vJyMjo9Z6s9lxSIPCRVHMY3x5RL2so0qgjJlz5+tNx4Qb0WgmbQ/ymRiOq+RxuueUWZFnu9yEQhBpqZ3Ko4d4A14KV+EhhRY9WhFNYIAASIo3oteqlrKlZOyNRCPSoRQi0QEBXZTsVrWg121QlRQuBHq0IgRYIushLUKehbZheS4KK/Q9FoHD0IgRaIOgiVyWBzog1qdq0ICHSlcctGH0IgRYIuoiPNBIT3renor/kJ0Yqfsye6LUa4iKC16FcoB5CoAWCHihtRWskibxEdSzznohA4ehECLRA0AOlrd30WJMi9Z+HQvihRydCoAWCHqTHmBTtmJ0fAOsZIClKFE0ajQiBFgh6oNFI5MQrJ6r5Ser6n90IC3p0IgRaIDgJpXzGSdFGosOUDzr2h8mgJSosoP03BAFACLRAcBI58RFoFMhbUzt742SEFT36EAItEJxEmF5LWoz/Pl0h0AJ/EQItEPSDv26OmHB9wAVTBApHH0KgBYJ+yEvwz/oNtPUMwoIejQiBFgj6ITbCQKwfqwoDlb3RE7NJH5Cca0HgEAItEAxA3jCt4HCDlrQgNXMVVvToQgi0QDAAw132nZcYqWpxpMEQAj26EAItEAxAeszwlmkHavVgf4gmsqMLIdACwQBoNBK5CeE+7WPQaciK820fJREW9OhCCLRAMAi5PmZzZMeHo9MG72cVF25ApxHFoUcLQqAFgkHIS4wgwQerNBjpdT3RaCSfxisIbYRACwSDoNdquGJyGibD0L7opGhj0AUaQr82tEEnZMdbxCclEAyBOVzPdyamDlqfIyHSwNVTM0JCfELZDy1JcN30TOJV7NE4mgj+t0kgGAFkxYdzVmFCv+/Fhuu5elqGV1Z2IAjlLt/Z8eEkRhm5cko6kUZRfW8ohEALBF4yLSuW8WnRvV6LNum5uiSDiBASm4RIoyLV+NRgQpoZcK16vHJqWkjMOEIZ8ekIBD5w/tgkUrtWCUYadVwzLT1gNZ+9Ra/VhKQLwWTQ9vLRJ0WFcdnkVLQi62RAhEALBD6g02q4rDiNhCgjV09LJyY89IQQICU69CrbFaVE9RHj7PgIzh+XFKQRhT6qCvQVV1xBVlYWYWFhpKamctNNN3HkyBE1TykQqE6kUcf3Z2YRH8LZEilBqgUyGBNOcg91v25mdn58gEczMlBVoOfMmcPbb7/Nvn37ePfddykrK+O73/2umqcUCAJCsGpteEuoCXRStHHQetUz8+Ixm0LLVRQKqBrZ+OlPf+r5/+zsbB5++GHmzZuHzWZDrxd/DIFALeIjDBh0Gqx2Z7CHAnQHBwcjLSaMpg5bAEYzcghY6Lm+vp5///vfzJ49e0BxtlgsWCwWz/OmpiYAmpubAzJGgWA0EaWxUd3WEexhoNVIpIXLQ/6Oo7U2OttaAzQqZWhubkay+2Zsuj8HWZaH3lhWmYceekgODw+XAfn000+Xjx8/PuC2CxculAHxEA/xEI9R/6iqqhpSPyVZ9kbGu3n44Yf54x//OOg2e/bsYezYsQAcP36c+vp6KisrefTRRzGbzSxZsqRfH97JFrTT6aS+vp74+HiffH7Nzc1kZmZSVVVFdHT/gYmRzmi/RnF9I5/Rfo3DvT5ZlmlpaSEtLQ2NZvAwoM8CXVdXx4kTJwbdJi8vD4Ohb/rR4cOHyczMZO3atcyaNcuX0/pEc3MzZrOZpqamUfnFgNF/jeL6Rj6j/RoDcX0++6ATExNJTEwc1smcTlfAoqeVLBAIBIL+US1IuH79ejZu3MiZZ55JbGwsZWVl/OpXvyI/P19V61kgEAhGC6rlQYeHh/Pee+9x/vnnU1RUxG233cbkyZNZtWoVRqO6Cf5Go5GFCxeqfp5gMtqvUVzfyGe0X2Mgrs9nH7RAIBAIAoOoxSEQCAQhihBogUAgCFGEQAsEAkGIIgRaIBAIQpRRKdDPPvssOTk5hIWFMXPmTDZs2BDsISnG6tWrufzyy0lLS0OSJN5///1gD0lRFi1axGmnnUZUVBRJSUnMmzePffv2BXtYivH8888zefJkoqOjiY6OZtasWXzyySfBHpZq/OEPf0CSJO67775gD0UxfvOb3yBJUq+He+W00ow6gX7rrbe4//77WbhwIVu2bKG4uJiLLrqI2traYA9NEdra2iguLubZZ58N9lBUYdWqVSxYsIBvvvmG5cuXY7PZuPDCC2lrawv20BQhIyODP/zhD2zevJlNmzZx3nnnceWVV7J79+5gD01xNm7cyIsvvsjkyZODPRTFmTBhAkePHvU81qxZo86JFKmIFELMmDFDXrBggee5w+GQ09LS5EWLFgVxVOoAyIsXLw72MFSltrZWBuRVq1YFeyiqERsbK7/88svBHoaitLS0yGPGjJGXL18un3POOfK9994b7CEpxsKFC+Xi4uKAnGtUWdBWq5XNmzczd+5cz2sajYa5c+eybt26II5MMFzcJWfj4uKCPBLlcTgcvPnmm7S1tY261bULFizg0ksv7fVbHE0cOHCAtLQ08vLymD9/PocOHVLlPKHTilgBjh8/jsPhIDk5udfrycnJ7N27N0ijEgwXp9PJfffdxxlnnMHEiRODPRzF2LlzJ7NmzaKzs5PIyEgWL17M+PHjgz0sxXjzzTfZsmULGzduDPZQVGHmzJm8+uqrFBUVcfToUR599FHOOussdu3aRVRUlKLnGlUCLRhdLFiwgF27dqnn3wsSRUVFbNu2jaamJv773/9y8803s2rVqlEh0lVVVdx7770sX76csLDQarulFJdcconn/ydPnszMmTPJzs7m7bff5rbbblP0XKNKoBMSEtBqtdTU1PR6vaamhpSUlCCNSjAc7rrrLpYsWcLq1avJyMgI9nAUxWAwUFBQAEBJSQkbN27kL3/5Cy+++GKQR+Y/mzdvpra2lmnTpnleczgcrF69mr/+9a9YLBa0Wm0QR6g8MTExFBYWUlpaqvixR5UP2mAwUFJSwhdffOF5zel08sUXX4w6H99oRZZl7rrrLhYvXsyXX35Jbm5usIekOk6nc9SU4D3//PPZuXMn27Zt8zymT5/O/Pnz2bZt26gTZ4DW1lbKyspITU1V/NijyoIGuP/++7n55puZPn06M2bM4KmnnqKtrY1bb7012ENThNbW1l536vLycrZt20ZcXBxZWVlBHJkyLFiwgDfeeIMPPviAqKgojh07BoDZbMZkMgV5dP7zyCOPcMkll5CVlUVLSwtvvPEGK1eu5NNPPw320BQhKiqqT7wgIiKC+Pj4URNHePDBB7n88svJzs7myJEjLFy4EK1Wy4033qj8yQKSKxJgnnnmGTkrK0s2GAzyjBkz5G+++SbYQ1KMFStW9Nvf7Oabbw720BShv2sD5FdeeSXYQ1OEH/7wh3J2drZsMBjkxMRE+fzzz5c/++yzYA9LVUZbmt31118vp6amygaDQU5PT5evv/56ubS0VJVziXKjAoFAEKKMKh+0QCAQjCaEQAsEAkGIIgRaIBAIQhQh0AKBQBCiCIEWCASCEEUItEAgEIQoQqAFAoEgRBECLRAIBCGKEGiBQCAIUYRACwQCQYgiBFogEAhCFCHQAoFAEKL8f0u1ddgRnjZ6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test points every 0.1 between 0 and 5\n",
    "test_x = torch.linspace(0, 5, 51)\n",
    "\n",
    "# Get into evaluation (predictive posterior) mode\n",
    "model.eval()\n",
    "likelihood.eval()\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(), qpytorch.settings.fast_pred_var():\n",
    "    # Make predictions\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "\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(rescale=True)\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'])"
   ]
  }
 ],
 "metadata": {
  "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.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}