{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fully Bayesian QEPs - Sampling Hyperparameters with NUTS\n",
    "\n",
    "In this notebook, we'll demonstrate how to integrate GPyTorch and NUTS to sample QEP hyperparameters and perform QEP inference in a fully Bayesian way.\n",
    "\n",
    "The high level overview of sampling in QPyTorch is as follows:\n",
    "\n",
    "1. Define your model as normal, extending ExactQEP and defining a forward method.\n",
    "2. For each parameter your model defines, you'll need to register a QPyTorch prior with that parameter, or some function of the parameter. If you use something other than a default closure (e.g., by specifying a parameter or transformed parameter name), you'll need to also specify a setting_closure: see the docs for `qpytorch.Module.register_prior`.\n",
    "3. Define a pyro model that has a sample site for each QEP parameter. For your convenience, we define a `pyro_sample_from_prior` method on `qpytorch.Module` that returns a copy of the module where each parameter has been replaced by the result of a `pyro.sample` call.\n",
    "4. Run NUTS (or HMC etc) on the pyro model you just defined to generate samples. Note this can take quite a while or no time at all depending on the priors you've defined.\n",
    "5. Load the samples in to the model, converting the model from a simple QEP to a batch QEP (see our example notebook on simple batch QEPs), where each QEP in the batch corresponds to a different hyperparameter sample.\n",
    "6. Pass test data through the batch QEP to get predictions for each hyperparameter sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import os\n",
    "\n",
    "import qpytorch\n",
    "from qpytorch.priors import UniformPrior\n",
    "import matplotlib.pyplot as plt\n",
    "import pyro\n",
    "from pyro.infer.mcmc import NUTS, MCMC\n",
    "import torch\n",
    "\n",
    "# this is for running the notebook in our testing framework\n",
    "smoke_test = ('CI' in os.environ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training data is 4 points in [0,1] inclusive regularly spaced\n",
    "train_x = torch.linspace(0, 1, 4)\n",
    "# True function is sin(2*pi*x) with Gaussian noise\n",
    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use the simplest form of GP model, exact inference\n",
    "POWER = 1.0\n",
    "class ExactQEPModel(qpytorch.models.ExactQEP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super().__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.ScaleKernel(qpytorch.kernels.RBFKernel())\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running Sampling\n",
    "\n",
    "The next cell is the first piece of code that differs substantially from other work flows. In it, we create the model and likelihood as normal, and then register priors to each of the parameters of the model. Note that we directly can register priors to transformed parameters (e.g., \"lengthscale\") rather than raw ones (e.g., \"raw_lengthscale\"). This is useful, **however** you'll need to specify a prior whose support is fully contained in the domain of the parameter. For example, a lengthscale prior must have support only over the positive reals or a subset thereof."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sample: 100%|██| 200/200 [00:11, 16.96it/s, step size=5.38e-01, acc. prob=0.912]\n"
     ]
    }
   ],
   "source": [
    "num_samples = 2 if smoke_test else 100\n",
    "warmup_steps = 2 if smoke_test else 100\n",
    "\n",
    "likelihood = qpytorch.likelihoods.QExponentialLikelihood(power = torch.tensor(POWER))\n",
    "model = ExactQEPModel(train_x, train_y, likelihood)\n",
    "\n",
    "model.mean_module.register_prior(\"mean_prior\", UniformPrior(-1, 1), \"constant\")\n",
    "model.covar_module.base_kernel.register_prior(\"lengthscale_prior\", UniformPrior(0.01, 0.5), \"lengthscale\")\n",
    "model.covar_module.register_prior(\"outputscale_prior\", UniformPrior(1, 2), \"outputscale\")\n",
    "likelihood.register_prior(\"noise_prior\", UniformPrior(0.01, 0.5), \"noise\")\n",
    "\n",
    "def pyro_model(x, y):\n",
    "    with qpytorch.settings.fast_computations(False, False, False):\n",
    "        sampled_model = model.pyro_sample_from_prior()\n",
    "        output = sampled_model.likelihood(sampled_model(x))\n",
    "        pyro.sample(\"obs\", output, obs=y)\n",
    "    return y\n",
    "\n",
    "nuts_kernel = NUTS(pyro_model)\n",
    "mcmc_run = MCMC(nuts_kernel, num_samples=num_samples, warmup_steps=warmup_steps, disable_progbar=smoke_test)\n",
    "mcmc_run.run(train_x, train_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Samples\n",
    "\n",
    "In the next cell, we load the samples generated by NUTS in to the model. This converts `model` from a single GP to a batch of `num_samples` GPs, in this case 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.pyro_load_from_samples(mcmc_run.get_samples())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "test_x = torch.linspace(0, 1, 101).unsqueeze(-1)\n",
    "test_y = torch.sin(test_x * (2 * math.pi))\n",
    "expanded_test_x = test_x.unsqueeze(0).repeat(num_samples, 1, 1)\n",
    "output = model(expanded_test_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Mean Functions\n",
    "\n",
    "In the next cell, we plot the first 25 mean functions on the same plot. This particular example has a fairly large amount of data for only 1 dimension, so the hyperparameter posterior is quite tight and there is relatively little variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEYCAYAAABxx2wUAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX4hJREFUeJztnXlc1NX6xz/DvsiOyCIg7vu+l+VWbrm0L+ZSWlaalpVX65a32y3NbFVvWXnD1NK0NFMzt3BXBHFLRAFZFJFVkG1Y5vz++PyGAQVkkIGRnvfrNS+Y+W7nu33Oc57nOedolFIKgiAIgtlhUd8FEARBECpGBFoQBMFMEYEWBEEwU0SgBUEQzBQRaEEQBDNFBFoQBMFMEYEWBEEwU0SgBUEQzBQRaEEQBDNFBFoQBMFMMalAf/nll+jcuTOcnZ3h7OyMfv364ffffzflIQVBEBoMGlOOxfHbb7/B0tISrVq1glIKK1euxEcffYSIiAh06NDBVIcVBEFoEJhUoCvC3d0dH330EaZMmVKXhxUEQbjjsKqrA5WUlGD9+vXIzc1Fv379KlxHq9VCq9WWftfpdMjIyICHhwc0Gk1dFVUQBMFkKKVw/fp1+Pr6wsLiFl5mZWJOnTqlHB0dlaWlpXJxcVFbt26tdN358+crAPKRj3zk0+A/iYmJt9RPk7s4CgsLkZCQgKysLGzYsAHffvst9u7di/bt29+07o0WdFZWFgICApCYmAhnZ2dTFlMQBKFOyM7Ohr+/P65duwYXF5cq161zH/TQoUPRokULLF++/JbrZmdnw8XFBVlZWSLQgiA0CIzRtTrPg9bpdOWsZEEQBKFiTBoknDdvHkaMGIGAgABcv34dP/zwA0JCQvDHH3+Y8rCCIAgNApMKdEpKCiZOnIgrV67AxcUFnTt3xh9//IH77rvPlIcVBEFoEJhUoFesWGHK3QtCtSgpKUFRUVF9F0P4m2BtbQ1LS8ta2Ved5UELQl2jlEJycjKuXbtW30UR/ma4urrC29v7tvtviEALDRa9OHt5ecHBwUE6OwkmRymFvLw8pKSkAAB8fHxua38i0EKDpKSkpFScPTw86rs4wt8Ie3t7AIzBeXl53Za7Q4YbFRokep+zg4NDPZdE+Duif+5uN/YhAi00aMStIdQHtfXciUALgiCYKSLQgnCH0qxZM3z22Wf1XYxao6GdT20gAi0IZkZiYiKeffZZ+Pr6wsbGBoGBgZg1axbS09Pru2j1yr/+9S9oNBpoNBpYWVnB09MT99xzDz777DOjh48ICQmBRqMx+xRMEWhBuAVhYWEYPHgwwsLCTH6s2NhY9OzZExcuXMCPP/6I6OhofPXVV9i9ezf69euHjIwMk5ehMkpKSqDT6ert+ADQoUMHXLlyBQkJCfjzzz/x6KOPYsGCBejfvz+uX79er2UzBSLQgnALvv/+e/z5559YtWqVyY81ffp02NjYYMeOHbj33nsREBCAESNGYNeuXbh8+TLeeuutcutfv34dTz75JBwdHeHn54dly5aVLlNK4V//+hcCAgJga2sLX19fzJw5s3S5VqvF66+/Dj8/Pzg6OqJPnz4ICQkpXR4cHAxXV1ds3rwZ7du3h62tLb799lvY2dndZHnOmjULgwcPLv1+4MABDBgwAPb29vD398fMmTORm5tbujwlJQWjR4+Gvb09goKCsGbNmmpdHysrK3h7e8PX1xedOnXCyy+/jL179+LMmTP48MMPS9dbtWoVevbsCScnJ3h7e+Opp54qzU2Oi4vDoEGDAABubm7QaDSYPHkyAGD79u24++674erqCg8PDzzwwAOIiYmpVtlMQq2Mym8isrKyFACVlZVV30UR7jDy8/PV2bNnVX5+fo22j4uLU2FhYSo8PFx5eXkpAMrLy0uFh4ersLAwFRcXV8slVio9PV1pNBr1wQcfVLj8ueeeU25ubkqn0ymllAoMDFROTk5qwYIFKioqSn3xxRfK0tJS7dixQyml1Pr165Wzs7Patm2bio+PV0ePHlVff/116f6mTp2q+vfvr/bt26eio6PVRx99pGxtbdX58+eVUkp99913ytraWvXv318dPHhQnTt3TuXk5KgmTZqob7/9tnQ/xcXF5X6Ljo5Wjo6O6tNPP1Xnz59XBw8eVN26dVOTJ08u3WbEiBGqS5cu6vDhwyosLEz1799f2dvbq08//bTS6zN//nzVpUuXCpeNHTtWtWvXrvT7ihUr1LZt21RMTIw6fPiw6tevnxoxYkRpeX/++WcFQEVFRakrV66oa9euKaWU2rBhg/r555/VhQsXVEREhBo9erTq1KmTKikpqbRcFVHV82eMrolACw2S2xVolJn5QqPRlPur/9Q2R44cUQDUxo0bK1z+ySefKADq6tWrSikK9PDhw8ut8/jjj5cK0ccff6xat26tCgsLb9pXfHy8srS0VJcvXy73+5AhQ9S8efOUUhRoAOrEiRPl1pk1a5YaPHhw6fc//vhD2draqszMTKWUUlOmTFHPP/98uW3279+vLCwsVH5+voqKilIAVGhoaOnyyMhIBaDGAv2Pf/xD2dvbV7rtsWPHFAB1/fp1pZRSf/75pwJQWubKSE1NVQDU6dOnq1zvRmpLoMXFIQgVsHr1alhZsaOt+v85LfR/rayssHr1apMdWxkxh8aN83v269cPkZGRAIBHH30U+fn5aN68OZ577jls3LgRxcXFAIDTp0+jpKQErVu3RqNGjUo/e/fuLdekt7GxQefOncsdY/z48QgJCUFSUhIAYM2aNRg1ahRcXV0BACdPnkRwcHC5/Q4bNgw6nQ4XL15EZGQkrKys0KNHj9J9tm3btnT7mqCUKpd7HB4ejtGjRyMgIABOTk649957AQAJCQlV7ufChQt48skn0bx5czg7O6NZs2bV2s5USFdvQaiA8ePHo127duVERM/Ro0fRvXv3Wj9my5YtodFoEBkZiQcffPCm5ZGRkXBzc0Pjxo2rtT9/f39ERUVh165d2LlzJ1566SV89NFH2Lt3L3JycmBpaYnw8PCbuiI3atSo9H97e/ubOl306tULLVq0wNq1a/Hiiy9i48aNCA4OLl2ek5ODadOmlfN36wkICMD58+erVX5jiIyMRFBQEAAgNzcXw4YNw7Bhw7BmzRo0btwYCQkJGDZsGAoLC6vcz+jRoxEYGIhvvvkGvr6+0Ol06Nix4y23MxUi0IJwCywsLKDT6Ur/mgoPDw/cd999+O9//4tXX321dEwHgAM/rVmzBhMnTiwnmEeOHCm3jyNHjqBdu3al3+3t7TF69GiMHj0a06dPR9u2bXH69Gl069YNJSUlSElJwYABA4wu6/jx47FmzRo0bdoUFhYWGDVqVOmy7t274+zZs2jZsmWF27Zt2xbFxcUIDw9Hr169AABRUVE1Tnk7d+4ctm/fjnnz5pV+T09Px8KFC+Hv7w8AN2Xg2NjYAGBmip709HRERUXhm2++Kb0mBw4cqFGZagtxcQhCJXh5ecHb2xs9evTAV199hR49esDb2xteXl4mO+bSpUuh1WoxbNgw7Nu3D4mJidi+fTvuu+8++Pn54f333y+3/sGDB7Fo0SKcP38ey5Ytw/r16zFr1iwAzMJYsWIFzpw5g9jYWKxevRr29vYIDAxE69atMX78eEycOBG//PILLl68iNDQUCxYsABbt269ZTnHjx+P48eP4/3338cjjzwCW1vb0mX/+Mc/cOjQIcyYMQMnTpzAhQsX8Ouvv2LGjBkAgDZt2mD48OGYNm0ajh49ivDwcEydOrVchVQZxcXFSE5ORlJSEk6fPo0lS5bg3nvvRdeuXfHGG28AoJVuY2ODJUuWIDY2Fps3b8Z7771Xbj+BgYHQaDTYsmULUlNTkZOTAzc3N3h4eODrr79GdHQ09uzZg9mzZ9+yTCbFKM93HSNBQqGm3G6QUE9BQUFp1oROp1MFBQW1UbwqiYuLU5MmTVJNmjRR1tbWyt/fX7388ssqLS2t3HqBgYHq3XffVY8++qhycHBQ3t7e6vPPPy9dvnHjRtWnTx/l7OysHB0dVd++fdWuXbtKlxcWFqp33nlHNWvWTFlbWysfHx/14IMPqlOnTimlGCR0cXGptJy9e/dWANSePXtuWhYaGqruu+8+1ahRI+Xo6Kg6d+6s3n///dLlV65cUaNGjVK2trYqICBAff/99yowMPCWQUL8f4DW0tJSubu7q7vvvlt9+umnN92XH374QTVr1kzZ2tqqfv36qc2bNysAKiIionSdf//738rb21tpNBo1adIkpZRSO3fuVO3atVO2traqc+fOKiQkpMrAbWXUVpCwzmf1NgaZ1VuoKQUFBbh48SKCgoJgZ2dX38UR/mZU9fyZ9azegiAIQvUQgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEodpoNBps2rTptvYxefJkjBs3rlbK09AxqUAvWLAAvXr1gpOTE7y8vDBu3DhERUWZ8pCCcEeTmpqKF198sXSaKm9vbwwbNgwHDx6s76LVGfoJXd3c3FBQUFBu2bFjx0onjv07YFKB3rt3L6ZPn44jR45g586dKCoqwv33319ubjJBEAw8/PDDiIiIwMqVK3H+/Hls3rwZAwcO/FvO6O3k5ISNGzeW+23FihUICAiopxLVPSYV6O3bt2Py5Mno0KEDunTpguDgYCQkJCA8PNyUhxWEO5Jr165h//79+PDDDzFo0CAEBgaid+/emDdvHsaMGVO63ieffIJOnTrB0dER/v7+eOmll5CTk1O6XD/Z65YtW9CmTRs4ODjgkUceQV5eHlauXIlmzZrBzc0NM2fOLDcecrNmzfDee+9VOgltRSQmJuKxxx6Dq6sr3N3dMXbsWMTFxZUuLykpwezZs0snYZ0zZ061Z4yZNGkS/ve//5V+z8/Px9q1azFp0qSb1r3VJLVVTSILGKz23bt3o2fPnnBwcED//v3LtfhPnjyJQYMGwcnJCc7OzujRo4fJZ3qvUx90VlYWAMDd3b3C5VqtFtnZ2eU+gvB3QT891KZNm6DVaitdz8LCAl988QX++usvrFy5Env27MGcOXPKrZOXl4cvvvgCa9euxfbt2xESEoIHH3wQ27Ztw7Zt27Bq1SosX74cGzZsKLfdRx99hC5duiAiIgJz587FrFmzsHPnzgrLUVRUhGHDhsHJyQn79+/HwYMH0ahRIwwfPrx0BpKPP/4YwcHB+N///ocDBw4gIyPjJqu4MiZMmID9+/eXTjf1888/o1mzZjfNZhMTE4Phw4fj4YcfxqlTp7Bu3TocOHCgdPxpfVnfe+89nDx5Eps2bUJcXFzpTN5leeutt/Dxxx8jLCwMVlZWePbZZ0uXjR8/Hk2bNsWxY8cQHh6OuXPnwtraulrnUmOMGuT0NigpKVGjRo1Sd911V6XrlB3vtexHxoMWjKWi8XhzcpQKDa3bT06OceXesGGDcnNzU3Z2dqp///5q3rx56uTJk1Vus379euXh4VH6XT/Za3R0dOlv06ZNUw4ODqWTpiql1LBhw9S0adNKv99qElqlVLmxkVetWqXatGlTOl62UkpptVplb2+v/vjjD6WUUj4+PmrRokWly4uKilTTpk3V2LFjKz2fshO6jhs3Tr377rtKKaUGDRqkPv/8c7Vx48Zyk/beapLaiqhsEtmy42Vv3bpVASjdh5OTkwoODq603GW54yaNnT59Os6cOYO1a9dWus68efOQlZVV+klMTKyr4gmCWfDwww8jKSkJmzdvxvDhwxESEoLu3buXm/Nv165dGDJkCPz8/ODk5IQJEyYgPT0deXl5pes4ODigRYsWpd+bNGmCZs2alZtvsEmTJuWa+UDVk9DeyMmTJxEdHQ0nJ6dS69/d3R0FBQWIiYlBVlYWrly5gj59+pRuY2VlhZ49e1b7ejz77LMIDg5GbGwsDh8+jPHjx1dYjqomqQWqP4ls2QlyfXx8AKD0Gs2ePRtTp07F0KFDsXDhwnKT65qKOpmTcMaMGdiyZQv27duHpk2bVrqera1tualzBKE2cXQE/n8KPLPGzs4O9913H+677z68/fbbmDp1KubPn4/JkycjLi4ODzzwAF588UW8//77cHd3x4EDBzBlyhQUFhbCwcEBAG5qems0mgp/u505FnNyctCjRw+sWbPmpmXVndj2VowYMQLPP/88pkyZgtGjR8PDw6PCclQ1Sa0xk8iWvUb6TBH9NfrXv/6Fp556Clu3bsXvv/+O+fPnY+3atRVO8FtbmFSglVJ4+eWXsXHjRoSEhJTOuisIQvVp3759ae5xeHg4dDodPv74Y1hYsAH8008/1dqxbjUJbVm6d++OdevWwcvLq9KZQXx8fHD06FHcc889AFA6WWx1Z0W3srLCxIkTsWjRIvz++++VlqOqSWpPnz59y0lkq0vr1q3RunVrvPrqq3jyySfx3XffmVSgTerimD59OlavXo0ffvgBTk5OSE5ORnJyMvLz8015WEG4I0lPT8fgwYOxevVqnDp1ChcvXsT69euxaNEijB07FgDQsmVLFBUVlU6IumrVKnz11Ve1VoaqJqG9kfHjx8PT0xNjx47F/v37cfHiRYSEhGDmzJm4dOkSAGDWrFlYuHAhNm3ahHPnzuGll14yevbu9957D6mpqRg2bFiFy281SW11JpG9Ffn5+ZgxYwZCQkIQHx+PgwcP4tixY5VWXrWFSQX6yy+/RFZWFgYOHAgfH5/Sz7p160x5WEG4I2nUqBH69OmDTz/9FPfccw86duyIt99+G8899xyWLl0KAOjSpQs++eQTfPjhh+jYsSPWrFmDBQsW1FoZXnvtNYSFhaFbt274z3/+g08++aRSYXRwcMC+ffsQEBCAhx56CO3atcOUKVNQUFBQalG/9tprmDBhAiZNmoR+/frBycnJaIvTxsYGnp6elXZO6dy5M/bu3Yvz589jwIAB6NatG9555x34+voCoLslODgY69evR/v27bFw4UIsXrzYqDJYWloiPT0dEydOROvWrfHYY49hxIgRePfdd43aj7HIpLFCg0QmjTWeZs2a4ZVXXsErr7xS30W545FJYwVBEBo4ItCCIAhmSp2k2QmCYP6U7aItmAdiQQuCIJgpItCCIAhmigi00KC5nZ5yglBTauu5Ex+00CCxsbGBhYUFkpKS0LhxY9jY2PxtBnkX6g+lFAoLC5GamgoLCwvY2Njc1v5EoIUGiYWFBYKCgnDlyhUkJSXVd3GEvxkODg4ICAgo7Y5fU0SghQaLjY0NAgICUFxcXG5gekEwJZaWlrCysqqVFpsItNCg0Y/iZvKB1QXBBEiQUBAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwxqUDv27cPo0ePhq+vLzQaDTZt2mTKwwn1TFhYGAYPHoywsLD6LoogNAhMKtC5ubno0qULli1bZsrDCGbC999/jz///BOrVq2q76IIQoPApHMSjhgxAiNGjDDlIYR6Jj4+HmlpadBoNFi3bh0AYO3atZg0aRKUUvD09ERgYGA9l1IQ7kzMatJYrVYLrVZb+j07O7seSyNUh2bNmpX+r5/FODU1FT169Cj9XSlV18UShAaBWQUJFyxYABcXl9KPv79/fRdJuAGlgIwM4MIFYNMm4P77V0Ojsfr/ZarcXwsLK3zwwWrodPVVWkG4s9GoOjJvNBoNNm7ciHHjxlW6TkUWtL+/P7KysuDs7FwHpRQAIDsbiI4GfvsN+PNPICUF0GqBggLg+nWgsBCwsACsrQFHR6Co6DjS0nrctJ+goHBYWHRHURHQqBHQqhUwaxYwYABgZVZtN0GoO7Kzs+Hi4lItXTOr18TW1ha2trb1XYy/DUVFwOnTwMGDwOHDQFQUkJVFAbazA7y8AB8fwMEBSE+n9dy7N/Dww0BxMRAeDhw4wOXr1wMajQWU0pX+9fSkkHt4ADodkJYGTJnCfffvD3TqBHTpAvTrB8htF4SbMSuBFkzPuXPA778DJ04A588DGg0tY72Fa2kJ5OcbBPXKFcDeHrjrLuDuuynYPj5ARAQF/LffgMuXvbB/vzd8fPzRt+8UbNmyAteuJWLmTC/cey+wfDkweDCPefIk3SM7d7JiGDcO+OsvCn6jRsDQoYDEFAWBmNTFkZOTg+joaABAt27d8Mknn2DQoEFwd3dHQEDALbc3pikgVM7p08D33wNHjwLJybRgXVzoonByooXbpAng5kb3w4kTFGalgPvu42/5+dzP8ePAjz8CNjaAuzvX0WgADw8tmjWzQaNGGiil4OBQCHt7W8THswzx8cDrrwMtWlCM9+2jcO/bx8pg7lxg6lRg1y6ua2EB9OkD9OrF/wWhoWCMrplUoENCQjBo0KCbfp80aRKCg4Nvub0IdPXR6YCLF2mNXrrEQN7mzcDly1xubw80awY0bw4kJdFq9vMzWM5WVrSK09KAkhLA2ZkibGEB+PpSwAMD6dJ4+mmgc+fyxz5/HggNBa5do2g3bw40bswypabSp713L3DvvcBHHxlcGkoB//sf8NZbdK289Rbw2mssQ2gocOwY1wkIAIYMYbkE4U7GbAT6dhGBrpxLl4BDh2gRAwZfb0gIcOQIv3foQCENCqKlmpAAtGkDjBpFH7C1NbeNiKCYZ2UBcXFA376Avz8F+fRpuiXS0ijC3t4U7FatKJZ6CxqgUOv/t7QEcnIo/O3b03VRUgJMm8ayu7kBkycDw4cbLORNm4AXX+R+FiwAnn3WcL4JCbSus7NZ7r59gW7dxLoW7jxEoBsg169ToBIS+N3Pj35hBwdg40Zg3TqKda9egKsrreSWLYHYWIroK68APW5ItMjJAd55h+tev053hoMDrWuA23fvTut5yRJg2DAKfHIyK4Ljx4GYGJalRQuKsVIU2ORkuiosLSnGFhYU/FGjgF9/pU/6xx957F69GHj092e2yH/+A6xYQet+wQLgoYcMwg8wuHn0KCsWpWiN9+7NoKNkhwjmjgh0AyElBdi6lZZto0YU0MBACqA+yPbbb/TpdupEYbW3Bx59FMjMZObFo4/eHHRLSwMWLgQiIwFPTwpejx60igcP5rH0FBVRnJs3p4Dm5gLbt9Mn3bEjXRbp6XRFFBYCrVsD999P94j+WJs2AadO0fednk5LPzWVQgwAwcGsLNzdub/evYEzZ+jqiI1lBfDOOyxbRRQU0B1y5gyvBcBz7t2bAU1BMCdEoO9gsrNpYaalMUti1CgK8rlzFKFDh+hfTkqiGI8fT+u3RQuK8d69tGwffpiCq0cpYMcOpsMdP07LNjOTQubpyeNmZRmsZz2Zmcx19vCguBYW8lhNmtB10qQJXQ75+Vz/6lWW1dqalu+DD/JYOp0he8Tfn2XYtw+YMQOYOJHfd+6knzwxkf7mbt2Ad99ly6GwEOjZE3j5ZVYMVaEU9xEaymAnQAu8bVta6y4utXW3BMF4RKDvMJSisB45QqFs3JhNdX2z/to1Cra1NV0cdnbA9OnA/v0UyMcfZw7zunUUTVtb4OxZrpufTzHPyADy8mgB29kxe8PTk38dHelb9vPjdz2xsRTwQYPoRnniCYp+Sgqt79BQCmFKCl0YTZvSNzxoEM/pu+9YodjbM9f5/vspwCtXsozt2gGvvspjjBtHQf/hB1rtTZsyz/qJJ3j+r7zC/djZ8RjTptHiri46Ha9RaCgrIoDXbuBA/hWEukIE+g4hPh748ENaeS1a0D3QuTM/dnZ0b5w7x2b7hQt0G0ybRsuzoIDBukuX6I+1sqK4urtT4IOCgD17uG8bG4qtnR390/fcQ39yly6VB9l+/ZX+4KgoWrcODrRiAfqUmzfn8ZQyCHZ6Ol0X2dmsXNq1Y6pcjx48j+hoiqGnJ/3bx44BzzxDd4uVFfc/diy3W7WKLYKICFrlkyZRpNPSuG7btrxOTz5Z80DhlSv0pV+9yspv6FD62AXBlIhAmzHJycDXX1P0XF2ZG1y2yV5YCPz0E8XDzo6ZEwA7ifz+Oy3g3r2ZdvbHH/Td9unDYJ6fH90Bhw7RCm3WjAKXnU23wpAhrAQSErhdWTQafoqL6WpIT+d2TZtSMP38aME6OfFvRgYDhGfPUuhnzaJLRk9MDPDeeyy/tTVFvWtXZok4OHCbrCyDlXzyJM9jyxb6wJ9+mu4QS0vggQeAb75hOdLSaIE7OFCgbWyAp54q786pCYWFdAFFRbESeeghHkMwHUqx4o6I4POkbzGWzQwq+7+1NWMtXbvyGbxTEYE2M0pKmMa2Zw8FZuxYWoeWloZ18vLYvNcHBGNiKG5NmrB5rhT9tWfP0hWSmAjMmUNh+vVX+ooTE+muuHiR/+vvbJMmzPho2ZLrBwSUDwRev05r/epVBvNatqRr5NVXKaxlycqiNf7VVxTzLl147JQUWsZ+frTO9YHJ7Gzg44+5n7AwCmGfPhTCpCT6qMPCuF8PD+Dtt3k+J09yH/feS+F86SVa4Dt3sgXw6afcpm1bVnBWVkzLK3tNa8rVq6wk7ezo4xehrj2U4vMbGsrv7dsz1uDpeetttVqmfUZE8H2xsGBMoVev2rnvdYUItJmgF93QUFrLY8dSKMuSkwOsXk2XhZUVO4Jcu0af8IMP0oosm2+cnk6BKiqiyLm58cH95RcKnlIUwMuXud1jj9EdoR/cyNqa1oetLY8VF8f1Ro1iANHLi9u+8krFKWsFBXRJPPNM+eyQK1c4nkdsLN0dJSW06h97jBXE8uW03jt0oLj+9ZfhPA4eZKWRn8/9TJlCP/Py5Sx3+/ascAYPptW8fDndOKGhrLAaNaIVfvYsBbVFi9q5f+npvDfe3jyPsql+gnEoBWzbxvt+992MSdzu9SwpYeV+5Aj3P2gQnw9zv08i0PVMcTFzfA8c4Ms9bhythLLohTkvjyIZG8vsglat2LxfupTNPj8/CuzhwxSpoiIKsrU1v2dk0HpwceExcnMp3B07GsbW0Gj4AJeUUHRiYlhGLy+KdWEhrVMPD4py//70Zfv6UoRbtKCFo9UyL/nVV1nhVEVcHFMAjx5lC+CBB2iR2tkBI0fyt+++Y/CxQwceY/16lrN9e1rls2fTok9O5rG9vNhCePpp7nf3blYYubl0GY0bxwrL2pq+6dp6USMjWbaJE+k2EowjIoLiPHbsrTNwakpJCUdejIigW+6hh8x3AC4R6Hpk3z76SBs3pm+0bJdogIK8ejUFOS6O4tK9u6GX3t69FKemTWlVKmWwlH18KKIxMbSWu3fneiNG0Ep99VWKuT7opkefN336NIOHY8ZQxJQCNmxgYM/bm66O3r0ZeIyPpztG/8nO5m9Dh7LSCAriuQUE3FoIU1KAb7+lxZuby0DnI48Ylm/bBixbRv+iiwuwZg1dG5mZFMShQ5mx4ulJt4ZS9N1rtcAXX7DcbdpQ7C0sgDff5PnOmHHriqS66CsUb2+2NoRbo9Uy3uLvz8qzrkhMZIvSxoYtKnOTDhHoeiAri8Kg1bK78o299goKgM8+oz/12jU+PC1b0jpzdOR2+fls3j/wAC3Y2bO538BACsOFC/Qvd+7MdT09aZHu2kUr8557KLwlJfwUFzMgWFDAYwUG0h3g6kpXQnIy86WTk/nbyJEVZ0QoxU4lY8fy4Y+NpSV++TL33awZK6SmTdl0rapzyPHjTBG0s+MAScOGGZZt3kyh7tWLLQZfX1Y0xcWG7IqOHdmkLSjg9k2aUNAPHeKxT52iFfXII6xU7r2XLp/aYu9eulJeeMH8m9L1yZUrdEXNmFE9/7IpyMqii7GkhK2u2qqsbxcR6Dpm5UrmIM+aZRCczExajJGRfKkjI/m7pydFQz+gUatWfJgHDaKwnTsHrF1L69rDgw9Vbi735+VFwXJxoRheu2bo0ffww7QUHBxoXf/wAy3joCC6MYqLaclGR9PidHXlvi5epNA4OxuCkfonwsaGwpuTQ9Fs3pzHbdmSH2trnsPRo9xvRgYF3MWFVrZ+dLwWLW4Ws4ULuW5cHCuksoL35Ze0hi0t6bcuLOQyfRB0xAie98GD7G3YtSuzQT7/nHnNCQm8XpaWrOj8/emeqC1BvXiRlcI//mEYz0QwcPIkWzCVxTHqGr07USlgwoTyAfL6QAS6jkhIAGbOpEX7+OMU5JwcPgjFxewSfekSRcXLiyKckEBfmasrrWalmD7m4MC84CtX6AOeOJEBwqgoBszs7Cg69vbcrnVrrt+zJ63ZvDxajGFhfADvv5/7TE6m8BYX09IdMsTQGWXFCvp/+/at+PwyMiiWMTG0QnJzeRz9OerH4PDxYTZH9+4UrPBwWsrXrxuClY0bsxUwdCgrqeJiWuXz5nE0u99+o7i+8QaFtLCQWSo7djAHWj/Ako8P/ft2djzm3r2snCZP5nm+8QZ98dnZvOYWFnQfabXAokW193KmpvLavPmmeYiQuRASQlfYpEn1XZKbuXaN+fUODrSo68tHLQJtYvTujP37GYxwd+fNjotj8/rkSVp4bm4UBHd3ilpCAoXZy4uCN2CAoUNKejrFrlMnWmhZWYYg3rlzFDVHR34iIylwLVtSJB0dOQ5Fo0YUSU9PNv0DArhORQLy5ZeG/OnKOH+egc6yo8oBFOeEBIpjaCjT0jIzDefcogUrgvvvp2ju3ct9ZWfzOnh58Xq1b08xf+MN7vf772n5Dx1Ky1ijod/8ySd5rXx9WeH060eBDAkB/v1vjtNhY8PttFr61e3seN6NG/Ma6nSMDwwezEqp7LlUZlnrg6s3nrue69d5v2fMYKVQthfm34nCQrrfduzgc6h371WUz6y/fjf+XxEVXfuK8qNvXN/Bga1PPz/DUAZl101OZgvI15fZOXWdoicCXcvop3cKDeXfuDg+DC4utOyUojjm5vJh9fam1di1K4Xr+HGK1uzZzPEdMoTrfv45LbxWrbg8K4sCY2NDIQsPp6C3bs0HrqSEx3vhBW4TGwv8/DMfxEcfrV5zWylmiAwaVHVEvaCA4zb/85/Vdw0kJ1PQjx5lJXL1Kq9DYCDT8u66y+Dyyc/nOaalUegmT6YIW1oyoLh+PV+eKVNYlscf5/l17UpxtrHhUKUffMDc6YgI3p/Vq7mPd9+l9e7tzVbG9esU+pdf5v7mzWPldTvodPRHL1/Oii4nxzBw1aBBd1ZurrGcP8+OUgDvi5UVXWivv24YKKu+yMszxEji4viM6VXO0ZHPWcuWXLZ+Pd+lsWPrbuhaEejbRKdjsOnIEd7s1FSKTWwsBbK4mD7ZsWMprPPns1nn6UmrcOZMugX++IMPxmuv0d3h6koB37iR23l50aURGEgRDAvjw336NF/yu++mvzooiBkEDg6MSl+4QKEPCqIFX10hUIo5yCNHsoNHVSxeTHG8saOKMehHmdu7l62NtDTD+BfDhnHfO3YwXS4ggNfD0pIiPGAAy7p7N1/6e+6hyyM3l5WhkxN97ffdR/dIixa8Vh98YMhd/vZbXvd27VgpXrgAPPcce2Ru2sRtxo9nBXc7REYyqKlvaZw+zU5JjRqxqd+QXCB79tDgaNWKsQNLSwrdpk30OZs716/TiLhwge95u3a0pLdv5/ukH9yrKoqK2Eo+fpwVsbG9WEWgjUQpCuqePbQMkpLYNG7WjFbRjh20+GxsKG4+PvTP/v47xbtHDwrAAw8wE+HCBQrJuHEU9osXaRmHhlJ4PD3ZpbltW+4nKopNfn2njHfeMYy4duAAm+bPPkuB27qVGQ1jxhhX45eU0CLWd1ypis2bKZaV+aZrSl4eg0cbNvCaWFnR5eDjw8qwa1fDhLVaLSutPn0YgI2Ppw/5hx8ozsnJvB8pKaz43Nxozb7xBgV9/nxWcJs3053TowfvT2gorehz52g9+fjwWJMn355/evt2WpJDhhh+u3yZAeQhQ2o3k6Q+yM9ny6t///KdrQoKeF/++c87c/KEs2f53hcW8n2PjaVL8uGHDe9gerqhE5ZGw+dWH3OxszP+mCLQ1SApiZbWiRN8yLy8KJg9e7L5s3s3rdS4OH7PyjLclF27GKDr3t2Q67xuHQWosJCZDs2aMb83J4fHGjaMov/ZZ7zxTk7MfOjalTe8c2e6AfSBi8uX6ZPt25fr7t5Nl8TIkcZnIxQVMWvi2WdvbS3GxbFCev55Y6+o8YSFMY85PZ1WS0EBz8/amtfSyoovS0YGlx04wPL7+lIo9u/ndtnZbIn07m0Y22PlSvqkn3+eL9dbb9F3/cQTrCRmz6b76fPPGZDdsoXPwOOP11xovv66Yotq82aW/7HHbv+a1QeRkXSlVZRX/uGHzGs3RQrb9es0dpKTaQjdOH6MHltbVtDu7oaZ6N3cjHtPlGJW0NGjfC5iY1kp+fmxtTxgAC3ssvsMCwvDnDlzsGjRIvTs2bPaxxKBvgGdjk2SH3+kpQxQJB94gC+xszMfhm3b6MIIC6NAWllxeWoq3Rb//CcFfeBAvvBFRXxwGzXiQ5KZyRu5ahWtXZ2O2+flcbAf/RCfhYW07lxc+GBMnGhIps/NNQzHGRBA32qPHrQ0a5ImlpdHC+ell8oPZlQRxcXA++/Tp1uX1pBStECXLKE126ePIZ1QnwXj789g39GjbFEEBVEcNm3iixkdTavc35/nERTEl9vdnecTE0Pf/V13saLatImWdnExxwoZN47HXLeO/9ekx5tSHCBq3ryb4wG7d/OZePzx279edUlcHF1yr7xy8/P344+GadVul6Iitm5OnGBrT6Ph+9KqFStkT0++ZxW9A1qtYYKK1FQaRJmZ5YOMSjGO07y5oWesPlAZG8tnKjub63p7c1/5+YZZgjQavpMdO/KcXV2B556biRUrlmDmzJn4/PPPq32uItBgM/rAAb7gAK3aJ57gy6+/yXl5FOVz55iV4OREd0PXrhSHJ5+kT/PsWTaln33W0P14506KgZcXRWPAAEaG9U3url1ZKYSHG1wXxcXMw0xOZoXw9NN8aAAu++knLvPy4kM2cCCtwpqSnk4L9fXXq5dhsHQpO3h4e9f8mLfL++/zWqSlscVhZ0crWT9Qjrc3r/V33/GFnjePL3BWFi3ts2f5knbqRBdIdjZfsldeYcU4aRKt7xde4PWeM4fPw7ffsqU0cCAFKSmJlqGxec7JyXz2Zsy4eZl+sKw7xZLOzGTHoTffvLnCjoykmD75ZM33r9MxE+fkScM8k/oWZW2jFM8nJobv76FDhnHBPTzYem7VivEgf3++/xoNxXv7dqbopaYCV67EIysrDUVFGhw8OAIFBSnw8vLC77//DqUUPD09EXjjFEY3IAIN4L//pRh37Vre6a8fVnLfPnZ8aNOGlmpqKgNYHh4U1bw8+ig9PGgVDR9OUT55kq6Opk1phbVvzyBVSgpvro0NrbkrV1jj+vnxxR85kuKvzybQW7M6HbulRkUZ/FkPPHD74xLHxdGFM2dO9aLq+/bRYijbs68+KC5mpTh7tqEbup8fffkaDZuuERG0lv38DMHX0aN5Lfv3p7siIoLWTrt2HKDH35/uh6AgXv+77mKPRn2HEwsL3se8PFq5V68yO+PRR7kPY9i6lZZ7v343L/vlF1pxXbvWxtUyHVot78Obb96cL1xcTJfZW2/VrFWn1fLdSk6mf75rV9P2ylSKlfv+/bTU3d153LLuvqIiwxAH8fF8ZvRj3xQVsdJwcABWry5bUA0ABY1Gg7IyeitJFYG+AX1NvWEDH4q2bRnxHziQwrRlC0X6r7/4oFy4wNr1iy9Yq//8M10e991Hq+6HH+iS0DeLgoIoGAUFfHitrCjA3btTFPTz9T3xBLcHuN7GjRQSjYYPyyOP3NoNUR0OHeJ+X3yxeq6KtDS6ZV599faPXRskJrJlM20ar9OmTbRkhg9nxXf2LEXD3Z1xgiNHeI+HDOHyCRNYSe7ezXsxezY7wwwZwmby0KEU3oEDeY1WrTKI9NGjtMz1FvAPP3AfxromFi7kPioKPL7/Pls15jqYD8BYycSJFc9a8/XXNCJ8fY3bp1IU5kuXeD1vN3umKnQ6vgdhYfzeqRNbXnpjpaTE0JP17Fnmb6ek0EjJzeVynY7f9aNAMni9BocOTYZOV3zTMa2srBAcHIzx48dXWTYRaBjGnf36azZPAgIoynZ2rCEzMri8VStmRISHU3AvXzbk3/70E2vPhx7idqtW0cIKD+fNc3Y2+LQDAijWbm4U9YED2aS1saHFpo8I5+dzv4cO0Trv3JmpPbX1sq5fT0F58MHqra/TsUff3Ln1n79ali1bKA79+/O7TsffIiM5WFGLFrRUExO53s6dbBk5O7MC9vKin/HMGYr72rV0hzz1FO/xk08yxW7IEN7rH35ga8PCwtCV+403eF+OH2dg+JVXqn+NsrPZU7OiSi8tjf7bl1+utctVq5w5w2ymhx66edlff/FjrJvmxAner0ceMd2sNUqxcj18mEZP//6M3yQk8PiJiVyvpIT3+PJlQ7aQrS23d3JikLFzZyYMBAWVb4FfvQq8++5xfPllj5uOHx4eju5V9fz6f0SgwSasrS2F0sPDcOGbNqUQW1rSQggJYb6tVssHR98d292dD6FWS1eB3k9tZUXRTklhk6d9e95oJyc+0M2acZ/e3nwY9W6L5GS+sKdOsck8ZkzVvfiMRT+y28CBjDpXl+XLaZnewm1WL3zyCf30ZVsVSlGoz55lHnrbtnSD/PorX8Lz5ynIFhZ0a7RsSbeU3uWzbBkF/tIl3p8XXqD7afx4irjeJ52RwQyPWbP4LKSl0UdvzOA/mzbRnVFREE2fNXI7MQZToBQt/IrcFzod3R7GuDZ0OrZePDwYfDWFKyMhge4prZaCbGXFZ0Gn4/LAQLovf/2VrsSCApYjIMAwnG6/fjdb9CkpfF+jotgivXCBx7K3P47IyB4ALADoYGFhAZ1OJwJ9uxw4QD/U8OFsyiYk0Bru1YsXv21bvjAjRvAFXbGCqTfXrxsGvE9P5wvq40PR9fSkpXXhAoWhSxc2ofUR4p07aY0VF1P8n3ii9gdriY2l+MyYYdxEqvv3szk3fHjtlqe2KCqiIPzznzd3HtDp+MJduMDmsr6C0ec95+XxHgJs2oaHGyw/jYYtJ30+6/z5dHmMGUNRnT2b62i1zB0fP56WVGEhv0+YwJf7VujFrqIgG8CxtefONa9R8X76iRVKRR2ZVq+m4RMUVL19paTwXkycWP1tqktRESu52FgaSvr7ZW3N97lbN7o3/vtfQw/dli0pxkFBNGRuTA28eJHGVU4On6/kZLonMzPZCtfPs6nUJezZ0wsuLv6YPHkKdu9egcTERBw7dgxNmza9ZdlFoG/g6lW6OnQ6iqO9PS94bCwtI42GzfwOHfjCf/ABLWb9UJoWFoaOEa6ufDh69mQzWd/d9YEHDB1A4uOZuhUfT6v8hRdu3TmkJihF8S8sZA61MS96QgLFaObM2i9XbZKUxBhAZe6A4mKmxl29SiHw9KQwf/YZ79+OHXQ3jBzJ67R2LZ+Bli0ZROzbl8/FV19x+yFD2AFp1ixDJbt0KS2snj257mefMZhadkyPyoiNpTvr6advXnbqFC35kSNv5wrVHllZQHAwz/1GMjLolpk+vXr7OnuWVu3s2bU74t+VK4bgsX6i5MBAxoecnOiO+vxziqq1NStgDw9W0p06lX9HiorYmj5zht/t7FjuAwd4HP0wBdbWTFO96y5a0n/9BaSkaHH4sA2mT9dgzBiFwsJC2FbTTykC/f8UFLA5lpDAGzRoEGvIS5fYbLl2jdZT+/Zs+u7YQQvZ358vY1YWrUxra96oxo3ZTMvOpiA0a0ZhtrHhDV20iLWwmxtvaM+eprOO4uOZL/3EExyrwxhyc9mVu67znWvK/v0UiLFjK1+noIAde0pKmEpnY0OXT4cOtLR27DCkUP3xB19G/dCq06cbpg2bMYMWWEhIeTH64QfDmNlK0TK8557q5Ut/+SUt9IpcIwsXMkBpDlb0V1/RR1/Rq7Z4Ma9HdSZr3b+fsYLnnqu984qIYMV5+TJbL61asaVrb0+x/uQTCmejRlzm4cGWbPv25fdTUMBW7YULfKaKiqgP8fHUA3d3vlNPP833Sz8GTkYG8/E7dWJL4upVPkMvv2x8L1GzE+hly5bho48+QnJyMrp06YIlS5agdzWcb7cj0PPm8YbNmsVOIf/9L90aLVow4OPmxqaOviehszMvtI8Pb0BqKtfp35/NPVdXw00aPZo+5tBQ+nAzM7n99Omm9ynm5jIH2MmJD5GxA/LodGx2v/pq/Y+Lawzr1xsmBKiKjAy+WN7edH1s2WKYkGDVKr6Q3bvTUlSKAq5/2TQaBnHfeostn/BwYOpUw7537aJATJrEbZcsoSV9q6BXQQGfv9mzb1528iRbCSNGGHtFahetlgJdkfUcHs5rVB1Lf+tWnu/DD99+mfQuwpAQ3teOHRn89vPjst27mb9eUkIjxdWV7pmhQ8u/F0VFXDcigkF9Nzc+Ezk5XK9RI1ae48fTNbJtGysZ/TDBY8fSst61i60hV1caAJ6eLM/Qocadl1kJ9Lp16zBx4kR89dVX6NOnDz777DOsX78eUVFR8LpFTllNBTosLAyPPTYHY8cuwsmTPZGSwpuSl8eXoXdv3hj9bB1WVrQkjx3jze7YkcJsYcEbYWvLm1BSwub0X39R6Fq0YA1aHX/k7VJYSP9gSgrdGTUdxOjzz/nA3Ylz6y1bRn95dSaFjY3lverbl5bzqlVs7Sxdymj+Sy9RNCMieP/1EyZs3sxW1JQpfPkuXqTPWc+pU7TGX32Vz8cnn7BVdasybdjA56oi3645WNGrV7O36o2pc0rRV/7mm7fex/btfE7HjLm9suh0dI/o85b79jXMMVhQQANlxw4KddOmtKLHjSsf6NbpGD/au5etZR8fvr8WFjS+AgJoVBUXM2jcqRNb1osXc9t772UL+Kef+L5fvMiWc2AgLXQfHxoLffoY3wo1K4Hu06cPevXqhaVLlwIAdDod/P398fLLL2Pu3LlVbltTgZ45cyaWLFkCO7uZaNny89KBiAoK2MQF+FK2bcuXdf9+1qJ9+hi6lRYWsqaMjmYTqLiYN2XsWD7IdTVCWW4u/abZ2WwmVyMGUSnLl7Pi6dSp9spXlyhlmFG8ur0dDx3iSzp6NC0xf3++sH/8wWuRn8+KOSqKwaZOndis/e03+qn1HZnmzzdYZUlJ7Lr/2mu0shYs4JgfVWV3VCV0J07QoquvYK1Ox5iJflzusmzezPjJrVw5u3bxGa0oNc+Ycvz2G4VVp6PPVz8MaFYWLfyjR+lXtrc39PbVS4NSvJa7drFjiocHDZHCQgb6c3N5T4OCWDEPG8bWVEkJK+6DB1n+nBy6U/Tjr4eGsgLu3ZvPQ3ExWwpHj7KivnG89FthNgJdWFgIBwcHbNiwAePKzBo5adIkXLt2Db/++muV2xtzIvHx8UhLS4NGo8GIESOQkpICV1cvjBz5O6KjFdq2ZRfMS5d4s8LCaFF7e7OZq9Xywtvb84a5uvLFHDKEqTt12alAKd78gwf5gDz8sKGDS01ZuZL+OGNS8MwRfS+2KVOqnvuwLDodA42JiWwCh4XxJV+8mL3YrlxhxbxzJ+9/RgYr4MJCCnt8PJ+PoUMpGr168WX/5BO+nF5eDCzPnVv16GZHjnC7siPe6fn4Ywp+fbBlC0X4Rn9tURFTUOfMqXr7w4eZxljTbt9K0a0QGkojqnt3GiMWFnQffv0171nfvvytUSO69/T+8L/+ogsjJoZlbt+ez0ZcHL8XF9Pq7dyZmT+9ehnuwdGjdPk1b85KMj+f91mnoxvM3Z0tLwsL3rvdu9nSGTiQPnYLC+PfTbMR6KSkJPj5+eHQoUPoV8Z5OGfOHOzduxdHjx4tt75Wq4VWqy39np2dDX9//+qdSJn2oaHrJbti6mnaVCEjgxfVz49+3BYt+PH3NwyEcjtjINeU3Fxae3/9xe+9etHCq41mb3Awz3HAgNvflzlQE5EGDP57Gxv+36QJm7Dt2lEcxoxhRebvz2fgwAH6XidN4gt85gyF3NmZFfbddzMLQD8A/NKllafU6Vm4kJbqjbGDDRvYpK4P19OiRRWL8Hff0cqsqsdgTAxdG9XN7riRkBA+90rRKHrmGbZkcnPpXz5yhO+CflaciRP53uqHCc3JocvCy4vrXb5MofXxocuiQwfen+BgGmOPPWYQ2xkzeE/d3Hj/HnyQgebERC6//35uU1jIVradHTOyLl0yxKjGjjVeL4wR6DpqqFePBQsW4N13363RtqtXr8bkyZNRXFxcpi+8/q8VHByC4ehIf2KPHqwZHR35MDg4sFbWp+CZEqXYXDt3jg+ZfghFBwf6tPQ51LWBPiVsyBBajA0FKysGgRcvpoV7o+VXGY6OfCkTE5memJPDZyE5maK8cycDYX/+yVbUffcZBj/y9mal3qOHIZior0yXLOH1feopupFefLHyMjzyCPd3Y9fxMWPYhK/rtMeTJxkYu5GsLIpUVeKcmUn3W3X80xUdd+tWtlJ1Ol675s1Z+QYH8x60akWXw8WLvD4pKaw0dDq6N/Lz+fehhwxuje7dKdJubrTo166l6E6dagiKv/su+zgEBtI61miYnvfpp6w4g4IYW1q3jr7udu0oxNeusSyjRtXOkAzVwaxcHLdjQQPA8ePH0UM/IVoZPv00HJMmdUdREWvDggJ+8vPp5sjN5Scnh9+Biuej0/9Wdl40PVXNbVd2HYAPVdu2vPGmyqTIzWWz+ZlnDH73hoZSBsuoJlkQhw/zxXZwYKwhOtoQdNKnbOmHMi0upg/y9GlWEE89xeWXLlFYEhPp6/bx4b2dMaPy52HxYgrAjW6zL76guNflTOGffsrMjRut/i++4LNT2SiIRUV0DcybZ5z7LyGB+dRt21JIvbzowtNoeC82baKl+sADbMH4+NAlCRhGujt0iC3Cbt24vrMzg7wbN7IF8tBD9DFv386KsGVLivrSpRyiVj+psoUF75u9Pe97fj7jAMHBfH86d2b5XF1vFuXCQlrfXl7Gx4XMxsUBMEjYu3dvLFmyBACDhAEBAZgxY0atBwn1Aq3veqnRWEApHf7xj3D4+3eHUnyYWrem1XW7fl1zJTSUzb+XXqo4p7WhsWcPBXXaNOPHE9Hp6OaIiaHv+cABVmh9+jC1LyiIgnHoEC2uLl0oyKdPc/1WrfjJzGQFn5lJMXdxobh37sx86bKim5xM63HKlPJliY1lkOt2Am3GoNXSv3tjJ6ArVyhuzzxT+bYVdcOvCn0nGHd3unLWr6dV6+tLX/E333Cd++7jtU5IYIXYrx9dFzt28Jr368f3dssWtmhGjmSOuo0NW8d5eXSNtG/PZVotK6FPP2Ul8MQTvJ8XLlD8H3yQOduhoTxvnY5JAE2bUpQ9PWkIREfT3XLtGr/b2NAl2r07K3hjMCuBXrduHSZNmoTly5ejd+/e+Oyzz/DTTz/h3LlzaNKkSZXbGivQly5dQq9eveDv748pU6ZgxYqbu2AWFHC8Bv24wQBr0g4d+EAYe7HNiYICPpyBgWz6/51ISWGzdfDgmk0vlZPDMSPS0hjB13fT/te/6KvWi0LPnqz0hgyhFXj1KkX9+nVadQUFfJE1Gjav4+Pppy4qoqWvn4tx6dLyEzXoqctg4YYNFL8bx2FZvJiulsoqu02beC7VmRKtqMjgTpo0icKfl8dzLyzk8xoZyWNpNHT9vfSSoUOQvjv/8OF8T/VTvo0ebZjx6JlnaOVu3EiLeOpUVnZvv817AbASunKFVvvAgdzm0CGOKlhSQmt85EiKsqsry3TwIK1qjYb3tm9f44ZSqAyzEmgAWLp0aWlHla5du+KLL75An2q8RTVJs9NqtbCxsSkNFFanC2ZJCS2w0FCDi6NFC8OsJ+aOvrtzSgqtmobaMqgOv//Opucjj9Rs/Af9oEr797MinzuXL2pUFHNjv/uOTW1ra77QfftStCMjKR6XL/OZOXaM+5o2zdCZJTmZlqDeyk5NZfZHWffCpk1M9atOrvftUlFlEB1Nf2xlo9Xprc1Jk6ret1KGe/HUUxS2zz+n0HbrRnfQjz8axq8BaCA98gi33biRrZqxY2kFb9tGi3XUKO7z119p/XbowOu8YgUr0pwcwyiHbm687iUldNVMncoKadEi5tTb2zNN7pln+N7v2sV3SD/mxl13mcZgMzuBrin1Nau3fhqcgwcNTRq9xRAQYB7dcgFae7/8QuvtkUfMc0S6+qC4mNclPp5B15oMCH/gAAXl+HG+5Pffz2DTwIEUcH9/g+9z9GhDoC0vj+Jy6RKFLC6Ox/f0pMU9bhyFo7CQPRZdXQ2+7kGDKDLLl1c8I0ttkphIH/yNQvzhh8wyqSgTJS+ProI336z6euq7ZY8YweuiH8zrlVfYwnjtNVZk7dszK+KXXwwZOZs3s2J88EGuu3Ur3UQjRvA5X7GCFe+IESz/N9/w/QwKYo/HK1d4DXNzWVHefTcrk1272FkpNpbX/6efDDN8FxfTqBk6lNuaGhFoE5CczJzJ+Hh+12goiB07GgZUqgvi4+lzzcrig/bgg38PP3NN0OloqUVE0O84YADvlzFDZX7wAYXZ2Zlpj1FRfJn37GFHCDs7Wn5Tptw8fOxff7GZnZlJ14ulJa2/Jk0oHEOG0K0yezZ7L+7dS2H/809WuPfcY7oA77JlLHPZvO2ICD7nlQVcFy2ioFb2vCUm0h/cqRP3odHwOkVFMa9ZPzdk69YU6SNHaB1PnEjXx6lTtJj1HVY6dTJ0L1+3jq2aoCDGkS5e5LXs25dWfVwcr2tGhmHo36IiGln64YDd3Hj/9efctSut5LoMygIi0HWCUoapceLjDWPPAow0N2lCi8DLi9aRs3P1RTw/n83fxEQGS9LTDQP7+PvzZb8TXC/mREGBIc9cKYplmzYG67YqoqOZaWBvT9Fp1owv+Y4dvB+pqWymz5lzs19Wp2MvxLw83s/hwxnc0pclLY2Bq2efNfh8jx7lPgsLDXNq2tmxrF263H6nKaWYfnnjZAJVDX/62298livyTF67xvxxd3emtuknWl25kv/HxrIl0q4dreYHHmAFMW4cz//IEQq6nR1dPG3acN2jR5mOd+QIj922LSvHI0cYXPTwYAckHx9ue/kyRV2no/Db2THYl5xM18VddzEI2aVL/baCRaDrmeJiPhBJSfybmclusPorfasrbmdH0QgI4Mfd3XzcKg2F4mIGn06epEiUxc6OllpQEAXY2poC/8QT3C4mhvenTRt2tGjalFZcUBDwzjs3DyqUlESf6ciRHP/D1ZUi27w5t3n9dVbglpZ0oUyaxOyKsgJaUMCynjhhmAXEzo4WYceO1ZsUWM+hQzxWWbHds4f+1ooCfwkJdFlMm1b+9/x8wwiCEycaUkaLi5lrnJdnmANQp6OvNyWF/uSWLQ0TW6Slsaegfg5PjYbXNyWFhs7kycyU0Y+vkppK8fb2Zm77hQsU7K5dub1WyzJbWxvGUjGniXpFoAXhNsjLY6vo4kXDOCwARWbbNoppaip9pV26sPJt0sQwpsv06TcH3zZvpqXXowfdAPn5DDru3k3fadeuDKZt2sTmflQUBfSFFyqfeSc/n838s2cNucKAQfybN2cFc+O4MZ9+Sn9w2Xz+hQuZ03wjJSUcK73spAmFhQzwZWRQ/Mq2QC5dYhfobt1o0WZnsyJp25bXRx/T8fdnxRYfz2syc6YhQyIsjK2TyZN5zcLCWL6UFK6vN1ri4ynIw4fzGum783ftynsWEMDKztgRH02NCLQgmIj0dFrJR49SdH79le6mFi1oRerdGYMGMc2urH9zwQIKkaMjreqVK+mLvvtuCnrr1rSUfX1pne/fT6sxLo5i2rw5RbxLl6rdZQUFrFxiY1nBlJQYlpWU0C87YgRdBG5urCA6dqSwOThQ4PWi9s03zEIJCDBMeJCURDeFkxNjIVlZFMbly1mxNG3KyszRkefm6clz0c+mfffdzIP286NvWl+B5Oay5eDtzXM4epSVYHo6j92kiWFUSisrumPS0ngdHR05YJWVFY81caLxk9rWFSLQgmBi/viDnTU++IBWWkwMRWLbNgpHWBgFads2Q1f07GyKWNlR4/bupcuhbVtDfnVSEi3IX3+lOA4YYBgvZO1ag1/a25sdN/r2paukOjGOzZsNgwllZjKP+5tvmIly/TrFT6ulkF++zJZC+/b0IRcU0DL286M/3sqKfvTz52nN+vsbxk63sKBL5PBhThM2ejTPbeNGw1jd1ta0pk+fZvbI6dP077u48DjnzvHaWFhQbPVWt78/K5WEBLYy5swxjAPes2fFg1GZEyLQglAHREfTIn75ZQrtDz/Qorv/fgrNzz9TOKZMYcqeRkPrLj+f6+jRjzO+ejVnX9GnS65cydSvyEjDuB8eHvRT+/rytz/+YDnS0w0+XD8/WvStWlG4y1rxn3xSfuKAFStoTd9obRYUsDXg50eBfPxxQ379pUvs5KLT0WXh6krLuWVLWrpDhjBYpx/8/7XXWNlYWhrENTGRbgitlhXFk08yT/nECVZiO3YYLPARI3ic3bu5jY0N9927N4/3/fcUev1MOuaOCLQg1BGZmfQ5d+1Kod23j5aunx/F6Pp1io6tLQNtPXsyg2LCBIptWaKiOAZG8+YG0fzll/LjDaelsTK4fJnfHR3pw+3YkcJ/8iTT1VJSaP2mp3MdT0/+jYpiNoOHB8Xs2DGO/+HsTP+5RsM0v88+oxtj0iRDWlp8PPOH/fwo1seP09r99VeW4epV/v3zT5YhKIi/WVmx0vH0ZMV1zz28XqtXs6XRujUrkbg4DlgFUOSHD2cltHkzXS/33EMR13e/3raNPviJE40b1bC+EYEWhDqksJA5vvru3gcOGIasbNqUTXEnJ4pjUBCXf/EFO6rcmJ3z6acMsu3bx+DfmTO0qisbZTE3lxbpmTMGX7OjI90S7dpReNPSaGUHB/P4Dg6GnnPt2tGKTUtjOhpgyCIKCjLMln3qFCuZgAD+7+xMS1qrpdj7+tKivn6dn8BA9tqbOpXnHhNDF0ZyMiuNY8fopx83jmVbtYrpc/368Xe9WPv7M3tmzBhDS+DECXZgGTGi8gCqOSMCLQh1jFL05Z46Rety3z76UQ8coJhduULhataM4tSzJ320Tz1Vfj+ZmbQs9QMYrVxpmEfP35+W5a3GH87NpVskMpIBPD179rBjU9Om9IcfPGhwbbRqRas2OZl+7scfZ4rcrl2GCXsTEiiOI0dyiNVWrejecHKipZyQwJS3Xr0Y5NSP/mhtzVZBu3ZsRdjaMjj4yy/sHPPbb7TEhw3j+qGh9N/Pm8cy6UlIMHSEGTnyzk09FYEWhHpi1y76T/UdY556ikG/Y8fog42JoVDpMyxGjqSwl3V3rFvHTI22belzXbaMPtdLlyiy165RnPSWckAAMxwqSifT6Si6ISEUPl9fVgIhIbRKnZzKp9tt306/d1YWy9y1Kyucgwcpxvn5FH47O7pInJz4f0ICt5s+nZ1KbgxYhoby2tx/PwOHycl0XWRmcvwYV1eKv7U1U+r8/AzbXr3KSqtJE1rTdTXdnKkQgRaEeiQ2lsG3khJaiVOnUpR37KDfNjWVgtyzJ7t6W1tzEKDXX6fVrR9JT+8C+ewzWtQ3CnBODt0gly5RxEpKuP6NlqWPD0Vx7lxapj//TAv3xnkG16+nVXz2LN02EycyXW/1amZLHDvGFkJGBisQFxdWEqmpFOaKJk/OzGQFpJ+w9dQpujrS09mdPSCAfmZ3d1ZCZbu2p6VxJDwnJ1Z0VU0ndichAi0I9UxeHt0ATk5MIRs1igL53XcMxi1dSuv4jTcYZDt9mi6F7t0ZQOzYkeL46KP0L1++TBdATSgqoh975kyK6U8/3TxF1YULDLrpJydu04aW/NWr9DPrXThNm9IK7tmT4nzkCCuPG7MndDoec+dOWrzZ2XS95OWxwhg4kP/b2bEbfatWhm31bhZnZ2Z3mHqWo7pGBFoQzIRffqHg/Pe/DM7t30+rOSCAboZvvmHz/ptvKMRffEHL2NGRgbYXXqCALV1aPj3OGDZvpgXfqhU7y8yeXX48j8JCBiZ79GBGR3Y2K5fiYlrnZ84wK2TWLMPYy//7H10fY8eWP9a1a6yI9u2jVezlZXCB6HSsfFxcWOGMGcOKQM/58yxr48b0gTcUi/lG7tg5CQWhofHQQ/Q76zMX2rRh3q61NYUoNpai9cADHC1v40ZapVu30h3w4ovsuejlRYvyiSeMnybtwgWK4e+/c6CtsuKs1XKcimefpdj+/jtF/No1rufiQqv+o4/4f1oau34//TQDf1otxXjtWopuQgKt7P79WdFkZNCCf/FF+uWzstgSaN6cx1eKfvWICGbAvPqq+XXNrk/EghaEOqCkhO6NXbsoRn36UOjWrmVesosLxdPGhn7oDh0YKDtyhK6OwkKm7jk5UahbtaKwjxhRtaAlJdFSv+suiv5LLxmWnTlDX3e3blzvl1+4r7vvplj27MnA3RtvGPK4d+9mJRMTw8oFYFbJX3/RZdG1K8toYUEL2M6OLQdHR1YE+nE7MjN5vMxMptVVMJVog0VcHIJgppw7x4GH9DOvDB5sELxhwyiYx47RJWFtTTdDUhLFOjOTwm5lRUv15EnDsJqWlhTBVq0YwPP3p3Bu3kxx/vlnZlkUFNDVsGcPy6DPnbazo5sjJobWbU4OUwQHDGDmxrZtPG7LljzW9et0hSQmcts332RZYmIozhERTCns3p3dvPXnsmcP/e2urmxd3CplsCEiAi0IZoxSzNDYvp3jaOTkUMhyc7m8Vy8GDg8fphU7YAB90EVFdJX4+FCwmzfnbCPZ2RS/jAxWAPrxw/Py6DvWapmi5uBAQU5NpTCmpNDabdyYlrt+7PK//mKZAgO5TloaRb+4mBkiXl60qD086BY5cIDH1s/H2K4dg5BBQYbBmcLDKdxDhtycPfJ3QwRaEO4A9u2j+8Dfn5ZvcTGt23PnaM22acMc6MREWq9ZWRTQwYPZnfrKFQqjm5thSi1fX66Xns5u1IWFhs4eBw9SuJ2c6Jd2dKRw64fvdHen26JVK+5n507+bd+eLo/u3Tn2R3o6RfjkSVr3Gg23Gz6cvuasLLpyLl+mKN91F7e9UzuW1DYi0IJwh5CURJdHXJxhxhAnJ6a4tW3LAGF0NMXW0ZF5wQAt78BACnBYGC3drCxarPb2/HvlCn27mZkUypYtKaTx8fRd62cW0Vvb+/bRYleK27/+Ol0RV68yNS89ndZ0UhIDgfb2rFS6dmXZLl/mtq6utJTLdjYRDIhAC8IdRF4eR5nLzKR4urhQAM+do783L49+3MaNOZ7H+vUMrLVtS7eITkc/cVERsy/OnqV7xMWFAqof/tPdnW6Jhx9mHrKPDwX6+HF2ohkyhIHA/v25zW+/0WWRnc3tHB05glx0NMVfb/k3bUq3jKnmT2xoiEALwh3IunUGt0VGBi3WCxc4EFNiIgNs+hnmL1xgup2dnaFH37VrFHl9N+xevWhhN2rEEeNWr6bv2Nqa6+XlMUOkpIQfZ2cuS0ridx8fVgJt21LIL16kH3v8ePq+hZohAi0IdyhnzzLjont3+ngdHAy96jp2ZG/AmBj20rtyhQG74mKu5+NDS/foUQ5i378/fdclJRwlTz9hKsBMkYULuSw/3zDUaKtW9H83b07B37KFQUUfHwYEG1qvvvpABFoQ7mBKSmjtFhZSUPfuZRZE797MsMjMpCVcUsJ5+/r2pW/42DH6igcPNsxWkp/PgGK/fhTryEh+t7RkJdC/P7uht23L9c+eZW/HggJa5qNGsRIQag8RaEFoAFy5QuvZzY3W7+ef0wett4RjYhg0vHw5DImJc9Cs2SI89FBPeHjQ53z6NPcREMAu28nJzNr44guKbkwMu1frB/8HKNQDBjTcbtbmgAi0IDQgkpNpGStFq/fQIWZteHoy42PDhpnYu3cJevWaiS5dPi+daNXHx5CR4evLbfXCa2vLrI42bejTlhS4ukPG4hCEBoS3N33PJSXstm1nB9jbxyM0NA2WlhqEhq4DAJw5sxZNm05C69YKPj6eaNIkEH37Sg7ynYxY0IJwB6Ipp7gaAKrMX2LGr/bfGmN0rRoTtdeM999/H/3794eDgwNcXV1NdRhB+FuyevVqWJVOLaLK/bWyssLq1avrpVxC7WIygS4sLMSjjz6KF1980VSHEIS/LePHj8fRo0crXHb06FGMHz++jkskmAKT+aDfffddAEBwcLCpDiEIAgALCwvodLrSv0LDwayChFqtFlqttvR7dnZ2PZZGEMwbLy8veHt7w9/fH1OmTMGKFSuQmJgIL0lcbjCYlUAvWLCg1PIWBKFqmjZtiri4ONjY2ECj0eD5559HYWEhbMtOmSLc0Rjlg547dy40Gk2Vn3PnztW4MPPmzUNWVlbpJzExscb7EoS/A7a2tqUZHRqNRsS5gWGUBf3aa69h8uTJVa7TXD/ZWA2wtbWVB0wQBOH/MUqgGzdujMaNG5uqLIIgCEIZTOaDTkhIQEZGBhISElBSUoITJ04AAFq2bIlGxk5LLAiC8DfEZAL9zjvvYOXKlaXfu3XrBgD4888/MXDgQFMdVhAEocEgXb0FQRDqELPo6i0IgiDcHiLQgiAIZooItCAIgpkiAi0IgmCmiEALgiCYKSLQgiAIZooItCAIgpkiAi0IgmCmiEALgiCYKSLQgiAIZooItCAIgpkiAi0IgmCmiEALgiCYKSLQgiAIZooItCAIgpkiAi0IgmCmiEALgiCYKSLQgiAIZooItCAIgpkiAi0IgmCmiEALgiCYKSLQgiAIZooItCAIgpkiAi0IgmCmiEALgiCYKSLQgiAIZooItCAIgpkiAi0IgmCmmEyg4+LiMGXKFAQFBcHe3h4tWrTA/PnzUVhYaKpDCoIgNCisTLXjc+fOQafTYfny5WjZsiXOnDmD5557Drm5uVi8eLGpDisIgtBg0CilVF0d7KOPPsKXX36J2NjYaq2fnZ0NFxcXZGVlwdnZ2cSlEwRBMD3G6JrJLOiKyMrKgru7e6XLtVottFpt6ffs7Oy6KJYgCIJZUmdBwujoaCxZsgTTpk2rdJ0FCxbAxcWl9OPv719XxRMEQTA7jBbouXPnQqPRVPk5d+5cuW0uX76M4cOH49FHH8Vzzz1X6b7nzZuHrKys0k9iYqLxZyQIgtBAMNoHnZqaivT09CrXad68OWxsbAAASUlJGDhwIPr27Yvg4GBYWFS/ThAftCAIDQ2T+qAbN26Mxo0bV2vdy5cvY9CgQejRowe+++47o8RZEATh747JgoSXL1/GwIEDERgYiMWLFyM1NbV0mbe3t6kOKwiC0GAwmUDv3LkT0dHRiI6ORtOmTcstq8PMPkEQhDsWk/kcJk+eDKVUhR9BEATh1ohTWBAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFBFoQRAEM0UEWhAEwUwRgRYEQTBTRKAFQRDMFJMK9JgxYxAQEAA7Ozv4+PhgwoQJSEpKMuUhBUEQGgwmFehBgwbhp59+QlRUFH7++WfExMTgkUceMeUhBUEQGgwapZSqq4Nt3rwZ48aNg1arhbW19S3Xz87OhouLC7KysuDs7FwHJRQEQTAtxuiaVR2VCRkZGVizZg369+9fqThrtVpotdrS71lZWQB4QoIgCA0BvZ5VyzZWJmbOnDnKwcFBAVB9+/ZVaWlpla47f/58BUA+8pGPfBr8JzEx8Zb6abSLY+7cufjwww+rXCcyMhJt27YFAKSlpSEjIwPx8fF499134eLigi1btkCj0dy03Y0WtE6nQ0ZGBjw8PCpcvzKys7Ph7++PxMTEBusaaejnKOd359PQz7Gm56eUwvXr1+Hr6wsLi6rDgEYLdGpqKtLT06tcp3nz5rCxsbnp90uXLsHf3x+HDh1Cv379jDmsUfwdfNcN/Rzl/O58Gvo51sX5Ge2Dbty4MRo3blyjg+l0OgAoZyULgiAIFWOyIOHRo0dx7Ngx3H333XBzc0NMTAzefvtttGjRwqTWsyAIQkPBZHnQDg4O+OWXXzBkyBC0adMGU6ZMQefOnbF3717Y2tqa6rAAAFtbW8yfP9/kx6lPGvo5yvnd+TT0c6yL86vTPGhBEASh+shYHIIgCGaKCLQgCIKZIgItCIJgpohAC4IgmCl3rEAvW7YMzZo1g52dHfr06YPQ0NAq11+/fj3atm0LOzs7dOrUCdu2baujktYcY87xm2++wYABA+Dm5gY3NzcMHTr0ltekvjH2HupZu3YtNBoNxo0bZ9oC3ibGnt+1a9cwffp0+Pj4wNbWFq1btzb759TYc/zss8/Qpk0b2Nvbw9/fH6+++ioKCgrqqLTGsW/fPowePRq+vr7QaDTYtGnTLbcJCQlB9+7dYWtri5YtWyI4OPj2ClErA27UMWvXrlU2Njbqf//7n/rrr7/Uc889p1xdXdXVq1crXP/gwYPK0tJSLVq0SJ09e1b985//VNbW1ur06dN1XPLqY+w5PvXUU2rZsmUqIiJCRUZGqsmTJysXFxd16dKlOi559TD2/PRcvHhR+fn5qQEDBqixY8fWTWFrgLHnp9VqVc+ePdXIkSPVgQMH1MWLF1VISIg6ceJEHZe8+hh7jmvWrFG2trZqzZo16uLFi+qPP/5QPj4+6tVXX63jklePbdu2qbfeekv98ssvCoDauHFjlevHxsYqBwcHNXv2bHX27Fm1ZMkSZWlpqbZv317jMtyRAt27d281ffr00u8lJSXK19dXLViwoML1H3vsMTVq1Khyv/Xp00dNmzbNpOW8HYw9xxspLi5WTk5OauXKlaYq4m1Rk/MrLi5W/fv3V99++62aNGmSWQu0sef35ZdfqubNm6vCwsK6KuJtY+w5Tp8+XQ0ePLjcb7Nnz1Z33XWXSctZG1RHoOfMmaM6dOhQ7rfHH39cDRs2rMbHveNcHIWFhQgPD8fQoUNLf7OwsMDQoUNx+PDhCrc5fPhwufUBYNiwYZWuX9/U5BxvJC8vD0VFRXB3dzdVMWtMTc/v3//+N7y8vDBlypS6KGaNqcn5bd68Gf369cP06dPRpEkTdOzYER988AFKSkrqqthGUZNz7N+/P8LDw0vdILGxsdi2bRtGjhxZJ2U2NabQmTobD7q2SEtLQ0lJCZo0aVLu9yZNmuDcuXMVbpOcnFzh+snJySYr5+1Qk3O8kX/84x/w9fW96YExB2pyfgcOHMCKFStw4sSJOijh7VGT84uNjcWePXswfvx4bNu2DdHR0XjppZdQVFSE+fPn10WxjaIm5/jUU08hLS0Nd999N5RSKC4uxgsvvIA333yzLopscirTmezsbOTn58Pe3t7ofd5xFrRwaxYuXIi1a9di48aNsLOzq+/i3DbXr1/HhAkT8M0338DT07O+i2MSdDodvLy88PXXX6NHjx54/PHH8dZbb+Grr76q76LVGiEhIfjggw/w3//+F8ePH8cvv/yCrVu34r333qvvopktd5wF7enpCUtLS1y9erXc71evXoW3t3eF23h7exu1fn1Tk3PUs3jxYixcuBC7du1C586dTVnMGmPs+cXExCAuLg6jR48u/U0/MqKVlRWioqLQokUL0xbaCGpy/3x8fGBtbQ1LS8vS39q1a4fk5GQUFhZWOHxvfVKTc3z77bcxYcIETJ06FQDQqVMn5Obm4vnnn8dbb711y7GRzZ3KdMbZ2blG1jNwB1rQNjY26NGjB3bv3l36m06nw+7duysdJa9fv37l1geAnTt3mu2oejU5RwBYtGgR3nvvPWzfvh09e/asi6LWCGPPr23btjh9+jROnDhR+hkzZgwGDRqEEydOwN/fvy6Lf0tqcv/uuusuREdHl1Y8AHD+/Hn4+PiYnTgDNTvHvLy8m0RYXyGpBjAkkEl0psbhxXpk7dq1ytbWVgUHB6uzZ8+q559/Xrm6uqrk5GSllFITJkxQc+fOLV3/4MGDysrKSi1evFhFRkaq+fPn3xFpdsac48KFC5WNjY3asGGDunLlSunn+vXr9XUKVWLs+d2IuWdxGHt+CQkJysnJSc2YMUNFRUWpLVu2KC8vL/Wf//ynvk7hlhh7jvPnz1dOTk7qxx9/VLGxsWrHjh2qRYsW6rHHHquvU6iS69evq4iICBUREaEAqE8++URFRESo+Ph4pZRSc+fOVRMmTChdX59m98Ybb6jIyEi1bNmyv2eanVJKLVmyRAUEBCgbGxvVu3dvdeTIkdJl9957r5o0aVK59X/66SfVunVrZWNjozp06KC2bt1axyU2HmPOMTAwsMJ5z+bPn1/3Ba8mxt7Dspi7QCtl/PkdOnRI9enTR9na2qrmzZur999/XxUXF9dxqY3DmHMsKipS//rXv1SLFi2UnZ2d8vf3Vy+99JLKzMys+4JXgz///LPCd0p/TpMmTVL33nvvTdt07dpV2djYqObNm6vvvvvutsogw40KgiCYKXecD1oQBOHvggi0IAiCmSICLQiCYKaIQAuCIJgpItCCIAhmigi0IAiCmSICLQiCYKaIQAuCIJgpItCCIAhmigi0IAiCmSICLQiCYKaIQAuCIJgp/wczZdAXYNwUMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    # Initialize plot\n",
    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "    # Plot training data as black stars\n",
    "    ax.plot(train_x.numpy(), train_y.numpy(), 'k*', zorder=10)\n",
    "\n",
    "    for i in range(min(num_samples, 25)):\n",
    "        # Plot predictive means as blue line\n",
    "        ax.plot(test_x.numpy(), output.mean[i].detach().numpy(), 'b', linewidth=0.3)\n",
    "\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', 'Sampled Means'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate Loading Model from Disk\n",
    "\n",
    "Loading a fully Bayesian model from disk is slightly different from loading a standard model because the process of sampling changes the shapes of the model's parameters. To account for this, you'll need to call `load_strict_shapes(False)` on the model before loading the state dict. In the cell below, we demonstrate this by recreating the model and loading from the state dict.\n",
    "\n",
    "Note that without the `load_strict_shapes` call, this would fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_dict = model.state_dict()\n",
    "model = ExactQEPModel(train_x, train_y, likelihood)\n",
    "\n",
    "# Load parameters without standard shape checking.\n",
    "model.load_strict_shapes(False)\n",
    "\n",
    "model.load_state_dict(state_dict)"
   ]
  }
 ],
 "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
}