tutorials/saasbo.ipynb

{
  "metadata": {
    "kernelspec": {
      "display_name": "python3",
      "language": "python",
      "name": "python3",
      "metadata": {
        "is_prebuilt": false,
        "kernel_name": "ae_local",
        "cinder_runtime": false
      }
    },
    "last_base_url": "https://devbig817.ftw3.facebook.com:8090/",
    "last_kernel_id": "a7bbe73b-ad66-4be3-9bac-27063ad0e3e4",
    "last_msg_id": "54931601-325c365835eb3d8a798d24d6_98",
    "last_server_session_id": "a4cb3c34-3d75-4f47-ad5e-0609188502cd",
    "outputWidgetContext": {}
  },
  "nbformat": 4,
  "nbformat_minor": 2,
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "501041cc-0473-4971-bff9-fd6e92e1eae4",
        "showInput": false
      },
      "source": [
        "## High-Dimensional sample-efficient Bayesian Optimization with SAASBO\n",
        "\n",
        "This tutorial shows how to use the Sparse Axis-Aligned Subspace Bayesian Optimization (SAASBO) \n",
        "method for high-dimensional Bayesian optimization [1]. SAASBO places strong priors on the \n",
        "inverse lengthscales to avoid overfitting in high-dimensional spaces. Specifically, SAASBO \n",
        "uses a hierarchical sparsity prior consisting of a global shrinkage parameter \n",
        "$\\tau \\sim \\mathcal{HC}(\\beta)$ and inverse lengthscales $\\rho_d \\sim \\mathcal{HC}(\\tau)$ \n",
        "for $d=1, \\ldots, D$, where $\\mathcal{HC}$ is the half-Cauchy distribution. \n",
        "While half-Cauchy priors favor values near zero they also have heavy tails, which allows the \n",
        "inverse lengthscales of the most important parameters to escape zero. To perform inference in the \n",
        "SAAS model we use Hamiltonian Monte Carlo (HMC) as we found that to outperform MAP inference.\n",
        "\n",
        "We find that SAASBO performs well on problems with hundreds of dimensions. As we rely on HMC \n",
        "and in particular the No-U-Turn-Sampler (NUTS) for inference, the overhead of SAASBO scales \n",
        "cubically with the number of datapoints. Depending on the problem, using more than a few hundred\n",
        "evaluations may not be feasible as SAASBO is designed for problems with a limited evaluation budget.\n",
        "\n",
        "In general, we recommend using [Ax](https://ax.dev) for a simple BO setup like this one. See [here](https://ax.dev/tutorials/saasbo.html) for a SAASBO tutorial in Ax, which uses the Noisy Expected Improvement acquisition function. To customize the acquisition function used with SAASBO in Ax, see the [custom acquisition tutorial](./custom_acquisition), where adding `\\\"surrogate\\\": Surrogate(SaasFullyBayesianSingleTaskGP),` to the `model_kwargs` of `BOTORCH_MODULAR` step is sufficient to enable the SAAS model.\n",
        "\n",
        "[1]: [D. Eriksson, M. Jankowiak. High-Dimensional Bayesian Optimization with Sparse Axis-Aligned Subspaces. Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence, 2021.](https://proceedings.mlr.press/v161/eriksson21a.html)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "executionStartTime": 1650234170354,
        "executionStopTime": 1650234172777,
        "hidden_ranges": [],
        "originalKey": "a6c53717-9ad1-4756-81e0-6d8756a23e11",
        "requestMsgId": "26933c08-82d6-439d-9fcb-6e358b080ab6",
        "collapsed": false
      },
      "source": [
        "import os\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import torch\n",
        "from botorch import fit_fully_bayesian_model_nuts\n",
        "from botorch.acquisition import qExpectedImprovement\n",
        "from botorch.models.fully_bayesian import SaasFullyBayesianSingleTaskGP\n",
        "from botorch.models.transforms import Standardize\n",
        "from botorch.optim import optimize_acqf\n",
        "from botorch.test_functions import Branin\n",
        "from torch.quasirandom import SobolEngine\n",
        "\n",
        "\n",
        "%matplotlib inline\n",
        "SMOKE_TEST = os.environ.get(\"SMOKE_TEST\")"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234172872,
        "executionStopTime": 1650234172878,
        "hidden_ranges": [],
        "originalKey": "a19ae9f3-2b63-411c-94d0-ffc2fdf2a45a",
        "requestMsgId": "f1e3c7f0-1afc-42e2-af59-5f5fae755ce5",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "tkwargs = {\"device\": torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\"), \"dtype\": torch.double}"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "08a3d790-52a5-4821-af21-1040f1a037f0",
        "showInput": false
      },
      "source": [
        "The time to fit the SAAS model can be decreased by lowering\n",
        "`WARMUP_STEPS` and `NUM_SAMPLES`. We recommend using 512 warmup steps and 256 samples when\n",
        "possible and to not use fewer than 256 warmup steps and 128 samples. By default, we only\n",
        "keep each 16th sample which with 256 samples results in 32 hyperparameter samples."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "customInput": null,
        "executionStartTime": 1650234172890,
        "executionStopTime": 1650234172961,
        "originalKey": "5cdef3eb-e589-4df0-b0d9-78919e0c6a7e",
        "requestMsgId": "363224de-347c-46a7-9c84-970cbb8e825d",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "WARMUP_STEPS = 512 if not SMOKE_TEST else 32\n",
        "NUM_SAMPLES = 256 if not SMOKE_TEST else 16\n",
        "THINNING = 16"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "af8beafd-352c-421d-8797-7660ddfa39f3",
        "showInput": false
      },
      "source": [
        "## Simple model fitting\n",
        "We generate a simple function that only depends on the first parameter and show that the SAAS\n",
        "model sets all other lengthscales to large values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234173003,
        "executionStopTime": 1650234173056,
        "hidden_ranges": [],
        "originalKey": "74eb6de2-5ad9-46a8-bb89-289ea7e0795a",
        "requestMsgId": "f506aa6b-904c-4a7e-8a38-0443e983df06",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "train_X = torch.rand(10, 4, **tkwargs)\n",
        "test_X = torch.rand(5, 4, **tkwargs)\n",
        "train_Y = torch.sin(train_X[:, :1])\n",
        "test_Y = torch.sin(test_X[:, :1])"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "cc9314b1-f255-4f7d-9f6d-eb349b34805e",
        "showInput": false
      },
      "source": [
        "By default, we infer the unknown noise variance in the data. You can also pass in a known \n",
        "noise variance (`train_Yvar`) for each observation, which may be useful in cases where you for example\n",
        "know that the problem is noise-free and can then set the noise variance to a small value such as `1e-6`.\n",
        "\n",
        "In this case you can construct a model as follows:\n",
        "```\n",
        "gp = SaasFullyBayesianSingleTaskGP(train_X=train_X, train_Y=train_Y, train_Yvar=torch.full_like(train_Y, 1e-6))\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234173089,
        "executionStopTime": 1650234188812,
        "hidden_ranges": [],
        "originalKey": "69db96ac-4974-4934-acc8-5a2662ec9dc4",
        "requestMsgId": "148855fb-cf0c-4fc5-8431-06a5e61c5da5",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "gp = SaasFullyBayesianSingleTaskGP(train_X=train_X, train_Y=train_Y)\n",
        "fit_fully_bayesian_model_nuts(\n",
        "    gp, warmup_steps=WARMUP_STEPS, num_samples=NUM_SAMPLES, thinning=THINNING, disable_progbar=True\n",
        ")\n",
        "with torch.no_grad():\n",
        "    posterior = gp.posterior(test_X)"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "5f4fa168-2662-499b-ac82-3ab122dfe2ad",
        "showInput": false
      },
      "source": [
        "Computing the median lengthscales over the MCMC dimensions makes it clear that the first feature has the smallest lengthscale\n",
        ""
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234188822,
        "executionStopTime": 1650234189056,
        "hidden_ranges": [],
        "originalKey": "02b17ef6-53b6-4b36-a805-66a1fdaa791c",
        "requestMsgId": "44a1f7c0-9649-4d89-8226-0405fdf88518",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "print(gp.median_lengthscale.detach())"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor([  4.4017, 207.1334,  96.3712, 190.6768], dtype=torch.float64)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "cf15a6ca-3377-40d1-9821-fad8f6600657",
        "showInput": false
      },
      "source": [
        "### Make predictions with the model\n",
        "\n",
        "In the next cell we show how to make predictions with the SAAS model. You compute the mean\n",
        "and variance for test points just like for any other BoTorch posteriors. Note that the mean \n",
        "and posterior will have an extra batch dimension at -3 that corresponds to the number of MCMC\n",
        "samples (which is 16 in this tutorial)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234189064,
        "executionStopTime": 1650234189076,
        "hidden_ranges": [],
        "originalKey": "cfe56132-ebf0-4897-a978-00c5d73a91b4",
        "requestMsgId": "898039a4-6ec8-46bd-a583-5a1614a3ccf6",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "print(posterior.mean.shape)\n",
        "print(posterior.variance.shape)"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "torch.Size([16, 5, 1])\ntorch.Size([16, 5, 1])\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "02b33f7f-4f31-432a-bac7-8cad1831e9a1",
        "showInput": false
      },
      "source": [
        "We also provide several convenience methods for computing different statistics over the MCMC samples:\n",
        "```\n",
        "mixture_mean = posterior.mixture_mean\n",
        "mixture_variance = posterior.mixture_variance\n",
        "mixture_median = posterior.mixture_median\n",
        "mixture_quantile = posterior.quantile(q=0.95)\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234189089,
        "executionStopTime": 1650234189100,
        "hidden_ranges": [],
        "originalKey": "5112b6e2-78af-4018-b7fd-6c0dfc9857e6",
        "requestMsgId": "b387d057-a497-401b-bfc2-ab427669c451",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "print(f\"Ground truth:     {test_Y.squeeze(-1)}\")\n",
        "print(f\"Mixture mean:     {posterior.mixture_mean.squeeze(-1)}\")\n",
        "print(f\"Mixture median:   {posterior.mixture_median.squeeze(-1)}\")"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Ground truth:     tensor([0.3282, 0.8348, 0.1427, 0.0913, 0.6747], dtype=torch.float64)\nMixture mean:     tensor([0.3274, 0.8415, 0.1443, 0.0924, 0.6733], dtype=torch.float64)\nMixture median:   tensor([0.3277, 0.8422, 0.1436, 0.0919, 0.6738], dtype=torch.float64)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1644277314184,
        "executionStopTime": 1644277314189,
        "hidden_ranges": [],
        "originalKey": "d9bec8be-2acd-40b8-aebb-612b62bbdfc3",
        "requestMsgId": "d9bec8be-2acd-40b8-aebb-612b62bbdfc3",
        "showInput": false
      },
      "source": [
        "## Optimize Branin embedded in a 50D space\n",
        "We take the standard 2D Branin problem and embed it in a 50D space. In particular,\n",
        "we let dimensions 0 and 1 correspond to the true dimensions. We will show that\n",
        "SAASBO is able to identify the important dimensions and efficiently optimize this function.\n",
        "We work with the domain $[0, 1]^d$ and unnormalize the inputs to the true domain of Branin \n",
        "before evaluating the function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234189112,
        "executionStopTime": 1650234189147,
        "hidden_ranges": [],
        "originalKey": "076b953c-6027-40c6-a807-b2bb82be770a",
        "requestMsgId": "15baa08e-ca35-4da7-a495-c63fe5d5779d",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "branin = Branin().to(**tkwargs)\n",
        "\n",
        "\n",
        "def branin_emb(x):\n",
        "    \"\"\"x is assumed to be in [0, 1]^d\"\"\"\n",
        "    lb, ub = branin.bounds\n",
        "    return branin(lb + (ub - lb) * x[..., :2])"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234189178,
        "executionStopTime": 1650234189438,
        "hidden_ranges": [],
        "originalKey": "10202f24-e344-47e2-8198-1737fccf7ac8",
        "requestMsgId": "98b6936b-2f06-4d1f-82c0-2f1bd660d0b2",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "DIM = 50 if not SMOKE_TEST else 10\n",
        "\n",
        "# Evaluation budget\n",
        "N_INIT = 10\n",
        "N_ITERATIONS = 10 if not SMOKE_TEST else 1\n",
        "BATCH_SIZE = 4 if not SMOKE_TEST else 1\n",
        "print(f\"Using a total of {N_INIT + BATCH_SIZE * N_ITERATIONS} function evaluations\")"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Using a total of 50 function evaluations\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "27fd793f-18ee-49cb-9aa5-c8cd78b0b807",
        "showInput": false
      },
      "source": [
        "### Run the optimization\n",
        "We use 10 initial Sobol points followed by 10 iterations of BO using a batch size of 4, \n",
        "which results in a total of 50 function evaluations. As our goal is to minimize Branin, we flip\n",
        "the sign of the function values before fitting the SAAS model as the BoTorch acquisition\n",
        "functions assume maximization."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234189470,
        "executionStopTime": 1650234891471,
        "hidden_ranges": [],
        "originalKey": "7342e71a-3fa9-4389-a98d-88d13beadc99",
        "requestMsgId": "269287e0-500f-474d-891a-5439487e9a77",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "X = SobolEngine(dimension=DIM, scramble=True, seed=0).draw(N_INIT).to(**tkwargs)\n",
        "Y = branin_emb(X).unsqueeze(-1)\n",
        "print(f\"Best initial point: {Y.min().item():.3f}\")\n",
        "\n",
        "for i in range(N_ITERATIONS):\n",
        "    train_Y = -1 * Y  # Flip the sign since we want to minimize f(x)\n",
        "    gp = SaasFullyBayesianSingleTaskGP(\n",
        "        train_X=X, train_Y=train_Y, train_Yvar=torch.full_like(train_Y, 1e-6), outcome_transform=Standardize(m=1)\n",
        "    )\n",
        "    fit_fully_bayesian_model_nuts(\n",
        "        gp, warmup_steps=WARMUP_STEPS, num_samples=NUM_SAMPLES, thinning=THINNING, disable_progbar=True\n",
        "    )\n",
        "\n",
        "    EI = qExpectedImprovement(model=gp, best_f=train_Y.max())\n",
        "    candidates, acq_values = optimize_acqf(\n",
        "        EI,\n",
        "        bounds=torch.cat((torch.zeros(1, DIM), torch.ones(1, DIM))).to(**tkwargs),\n",
        "        q=BATCH_SIZE,\n",
        "        num_restarts=10,\n",
        "        raw_samples=1024,\n",
        "    )\n",
        "\n",
        "    Y_next = torch.cat([branin_emb(x).unsqueeze(-1) for x in candidates]).unsqueeze(-1)\n",
        "    if Y_next.min() < Y.min():\n",
        "        ind_best = Y_next.argmin()\n",
        "        x0, x1 = candidates[ind_best, :2].tolist()\n",
        "        print(f\"{i + 1}) New best: {Y_next[ind_best].item():.3f} @ \" f\"[{x0:.3f}, {x1:.3f}]\")\n",
        "    X = torch.cat((X, candidates))\n",
        "    Y = torch.cat((Y, Y_next))"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Best initial point: 13.121\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "4) New best: 4.679 @ [0.479, 0.238]\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "5) New best: 0.464 @ [0.131, 0.795]\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "6) New best: 0.457 @ [0.536, 0.149]\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "7) New best: 0.400 @ [0.542, 0.155]\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "10) New best: 0.399 @ [0.543, 0.150]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "a9704a99-0712-40bb-a263-6798e0925291",
        "showInput": false
      },
      "source": [
        "## Plot the results\n",
        "\n",
        "We can see that we were able to get close to the global optimium of $\\approx 0.398$ after 50 function evaluations.\n",
        ""
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234891537,
        "executionStopTime": 1650234892438,
        "hidden_ranges": [],
        "originalKey": "22f74895-3c95-410c-a49f-bd6a0094695d",
        "requestMsgId": "fd0d7aa7-8d55-4942-adc2-de356666ac84",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "Y_np = Y.cpu().numpy()\n",
        "fig, ax = plt.subplots(figsize=(8, 6))\n",
        "ax.plot(np.minimum.accumulate(Y_np), color=\"b\", label=\"SAASBO\")\n",
        "ax.plot([0, len(Y_np)], [0.398, 0.398], \"--\", c=\"g\", lw=3, label=\"Optimal value\")\n",
        "ax.grid(True)\n",
        "ax.set_title(f\"Branin, D = {DIM}\", fontsize=20)\n",
        "ax.set_xlabel(\"Number of evaluations\", fontsize=20)\n",
        "ax.set_xlim([0, len(Y_np)])\n",
        "ax.set_ylabel(\"Best value found\", fontsize=20)\n",
        "ax.set_ylim([0, 8])\n",
        "ax.legend(fontsize=18)\n",
        "plt.show()"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<Figure size 576x432 with 1 Axes>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfwAAAGUCAYAAADQ9ZAkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3gU1f7H8feSBBJKKIYOCtgFxXKxK6COZQR1vHr9Wa7l2q+9d0WvBRV799oLYh0VHJVRAZUiYuNiQRGsECEC0gIkkN8fOYHNkrKb7Oxudj6v58mz7MyZmW+Oa757zpxzJlJRUYGIiIhkt2bpDkBERESCp4QvIiISAkr4IiIiIaCELyIiEgJK+CIiIiGghC8iIhICSvgiIiIhkJvuAEQkOJbtnAg8CZzke+5T6Y5HRNJHCV+kAaISaazlwGxgDDDC99yFaQgv2jjgSODToC9k2c5TwAkxm1cAC4BpwEvAq77nrgk6lqiYxgMD6yji+J77elT5dsAw4DCgK/An4AFX+55bnJqoRYKhhC/SOM8DVQkjAnQEDgCuABzLdnbyPXdFuoLzPfdn4OcUX/ZK4Afz79ZAb+AQ4EVgmmU7R5i4UunIWrZPrfqHZTutgAnANsD95kvKFsDFwD6W7QzwPffP1IUsklxK+CKNM8P33Fditj1o2c5I4GiT6EbVdrBlO/m+564MPsyUGud77pSYbddZtvMv4BFgrGU7O/ue+1eqAqrhv1FNzgW2A/7te+5DVRst2/kSeM18kbko2EhFgqOELxKMiSbhd6zaEHUb4HhgAHCiSSQnmv1/Ay4DdgOKgD+AT4BrfM+dGXWeM4CHAAvoDlwCbAoUA6+a8qXUcg/fsp1iYBZwBHC3OU8u8D/gMt9zJwZRIb7nPmHZTlfgRhPz1bWVtWynOdAyjtOu8D13dZJCPBZYBjwes/114DfgWMt2LvY9Vw8gkSZJo/RFgrGdef2ihn1HmqR+MVCVhLcz3cm7AfcA/wIeM8n4E8t2ekYdv8q8ngBcZRL66cBM0wK9qZ7YVgEFwHsmwV0I3GJiHmvZTpfkVEGN7jbjHI6vp9wxwKI4fo5J5OKW7TS3bGeDv3uW7bQG+gJfxn6BMAn+E6CzuT0h0iSphS/SOPlmoBfmHn5n4O/AKcCTvud+XMMxewG9Yrq0twE+Bm73Pfe9qo2W7SwwrfkTTMsYoKqFeQCwhe+5i03Zl0wr/3CTxGtTAewIXOV77s1R14oANwMHVn0RSTbfc5dbtjMNGGjZTjffc+fWUvQdYHAcp/wunutatjMMOAnoCayxbGcScIXvuZNMkT7m9ddaTvFrVLnZ8VxTJNMo4Ys0znXmJ1oF8LDpnq+JF3v/2vfcUVX3+i3byQFamR64OaZIrxrO83RVsjfnWGnZzkyTzOtTAdwbs+1L89o9juMbY5557QjUmPDNiPhkjorfz/R8zDW3Uy4Gxlu2s4/5UtbGlFtey/FV2wuTGJNISinhizTOf4GRUe/bAlsDZwJDzYj0T2KO2WCEuulmPt/0DGxZw+22mv5fnVXDttI4/7/+w/fcZTUcC5AXx/GNUWBeUzFY8VKTzMf5nrvWbBtj2c7HwFjgDmCXOM4TMa+6fy9NlhK+SOPM9j13fMy2N8yc9BnASMt2tvQ9tzxq/5IaznOLSU5fAKeZLuTVpqv/gVquvaqW7fFozLGNtbF5/b22ApbttIhqdddlqe+5tf4uvudOrWW7b9nO98AAy3baAFU9LrVds7V5TdnMApFkU8IXCYDvucWW7YwzI+G3AL6praxlO3nAWWYQ2iDfc5dE7cuqgbWW7XQHtgc+raGHIdrRtSxsFOukRow3mG/+27Q29+Uror6MxKq6pVJTr4pIk6CELxKc5ua1oJ5yReae/afRyd6oa5W4pug60z1eX5Ju9KA9M9tgf+Bn33MnxOyLmKmMpUCJ77llZr799pbtFFRNa2T9mIrdzHl+ieeXFMlESvgiATDT6AaZVvuMeoovAMqBjS3biVTN87ZsZ6uqOfpxfGkIhGU7bc0Ss/Mbu0ywZTuXmjEKk80CPLVK0qC9CPAEMMeynR1iehRONL/XS77nlpltT5kpkWcAd0WV/acZYBg7OFOkSVHCF2mcfpbtHBH1vpXpJj7VdBUfV9c9ZiqTW7llO68CRwHPWbbzDrCVOccJwJvAvpbtnAS8EfyvVI1jutavAIbHecxgy3Z6mH+3MN3hRwHbAh8C/0jFevq+586zbOdGk6inWLbzGLDYtNZPMWMILo465BEzr/92y3Y2Mc8f6GemOH5lBviJNFlZdX9QJA2OBV6O+nnMtB7HAbv5nvtCnOc5E3jaTB97ENjTJMa3zfz75mZgX4eAf5/arI2jTJWbo+rjadNinmvqal/fc/8IMM5qfM8dZr5slJilcR8xixk9AOzke+6vUWVXmX13mS86TwDHmSmWg3zPrW3KnkiTEKmo0CwTEamdZTvPAeN9z30s3bGISMOlvUvf3Ke8DdjdzP/9DrjB99y30h2bSNiZ6XEDgRvSHYuINE5au/TNSNl3gaVmxGxR1eNGzRcBEUmvXsBw33O/T3cgItI46W7hdzLzXkdWLTVq2c7jZqRs/3jXyRaRYJin9M2Mo6iIZLi038O3bGc88Kd52tdSswDJlUDfVA7uERERyWbpbuED/AN428xFrjCjaY+oL9kX//5LRW5OJoSfvcrLy8jNDXpZdUl1Pb8xOp9TTls/2H/IwaU8+diilF0/HfRZDp7qOHhFXbpF4ihWq7RmTMt2mpt7+N8BB5lnc58AjLZsZzffc2tdsCQ3J5eiLt1SG3DIlBTPVR2nQKrruVPM0+7XVhRQ1CUt6/qkjD7LwVMdZ750N5H3M+tqH+R77nyz7SHLds4E/lXPM71FpAEKYnJ7aWltJUUkm2TKwjux3RS5egylSDCU8EXCKd0t/Elmvezhlu2cDyw3K1ttaQbxiUiSKeGLhFNaW/i+5y4GDgA2Ar43A/fOBI70PfejdMYmkq2U8EXCKd0tfHzPnQ4MSXccImGhhC8STplyD19EUkQJXySclPBFQkYJXySc0t6lLyKpVVPCr6iASKOW9MgcFRUVLF68mPLy8nXbVq4uZ8GCBWmNK9upjhuuXbt25OUFv2iREr5IyOTlQU4OrFlT+X7tWigrg+bN0x1Z461evZqFCxfSoUMHmkf9QmVlq8nLy4JfMIOpjhtm7dq1zJ8/n4022ijwpK8ufZEQytZu/b/++otOnTpVS/YimaxZs2Z06tSJRYuCX95aCV8khGIT/sqV6Yok+Zo10581aVqaNWtGJAX31PR/hkgIZWsLX0Rqp4QvEkL5+dXfK+GLZD8lfJEQUgtfJHw0Sl8khJTws8efCxcy6qXXmPbZF8xfUEIkAu3atqVf32044vBD2GzTPhscs6CkhONOPJ21a9fy0H13stmmveu8RrzlE4nlmedG8ezIF6sdn5ebS8eORew8YCeOPfpI2rVtu8E15s9fwCvum3z2+ZfMn7/ADHorYtedB3DI0IPoWFSUQO2FixK+SAgp4WeHRYsWc9Z5l1CxtoIj/n4om/bpTU5OM3766Rdeee0NJk7+hBHDb2DLLTavdtw7Y98np1kzmufl8a7/Ppttekqd14mn/KLFizn3wisSjuXaKy+lU6eOAKxYsYIZX3/LCy+9ypdf/Y+H7ruD3Nz1aWrqp59x4/ARtGvblsMOOZhNzRePH3+cg/vmW7z1zliuvvxidtyhf6PqNVsp4YuEkBJ+dvDe8fnzz4XcdfvN9Ou79brt/bftx267DODEU/6N+8YYLr/kgnX7KioqGOt/wPbbb0dBfj7vj5vAqSefQPNa5oDHW/6dse8nHAvAJpv0ZOOePda932H77aioqODZkS8yfcY37Lj9dgDMm1fMjcNH0LtXL4bfeC0FUR/i/tv248D99+Oyq4Zx4y0jePTBuykq2qhRdZuNdA9fJISU8LPDnwsXAtC9e9cN9nXq1JFRzz+xQYL9/IuvKP5jPoP23oPBg/Zi6dJlTJ48tdZrxFt+4cJFCcdSm969NwFgyV9L1m0b9fJrrFq1mksuPKdasq/SsmUBl150LsuWL+fl196I6zpho4QvEkJK+Nlh0z6VXdoPPvw4S5Ys2WB/YZs2G2x7+933KCjIZ689dmOXATtRWNiGd/z3a71GvOX79O6VcCy1+e23uWBa/1WmTJ3GVltuQY/u3Wo9rmeP7my15eZ89PGkuK8VJurSFwmhsCX8TF95r6KiYccdYO2D//44xn/4MRMnTaHvNlux3Xb92K5fX/r13ZqcnJxq5ZcsWcKkKVPZZ9Be61rJ+wzamzfHvM2CkpINBrwlUt7adxAfjP8o7lhqsnz5cj7/cjovv/o61r6D6N1rk3XbFy5cxO677lzvOfr07sW3333P6tWrM/6/e6op4YuEUNgSfrbKzc1lxPD/MPa9cbw/bgIzvvmOL6fPAKBNm9Y4hwzhmP87Yl2yfe+DCZSVlXGAte+6c+y/32Bef/Mtxr43jmP/78hq50+kfKKxVDn59HM2+L322H0XTj/lpHXvS1euAqCgIH+DsrGqvpgsXbaMjTp0qLd8mCjhi4SQEn72yM3NxT7Qwj7QonTlSr79diafffEV730wnmeeH8X8BSVcdP5ZYLrnO3YsYuOePfjL3B/v1LEjPXt0Z6z/AcccdUS1JV4TLZ9ILFWGXX05nTt3AqCsrIz58xcw2nuHU844hxuuu4qtt9qCwsLK2wFLly6vtz6WL19OJBKhbWFhUuo3myjhi4SQEn52KsjPZ8cd+rPjDv057ph/cP7FVzD2vQ84+8xT+HHOT/z08y8AHHH0CTUeP33G1/Tfth8A33w3M6HyicTSokWLdeV69uxebZT+1lttwe677cwpZ57HvQ88wkP33UHzvDy6dOnMrB9n11sHc+b8TLeuXapN55NKqhGREMrmh+fUZPXq7Ht0a1lZGV9/8x3NWzRnm6223GB/QX4+O+3Qn9mzf2Lx4r945933iEQiXHvVpbSM+QCsXr2aYTfeyrtj31+XwBMpX1ZWxlfTZ9CyVau4Yqlq0dcmLy+PTfv0YuKkT6ioqCASibDXHrvy8qtvMOvH2TUuJgTw62+/M/OHWRxz1BFx1WHYKOGLhJBa+E1fBXDzbXfSsqCAB++9g5Ytq/9HXbNmDdP/9zWFhW1o3boVEz6cyA79t2PP3Xet8Xy7DNiJjz6ezNlnnkakWSSh8rl5udx2x720bNmy3ljimR9fXl7Oj7N/omijjdbdMjjy8MMY64/j9jvvY8Tw/9CmTetqx5SuXMmIu+6nQ4f2HO4MrfcaYaSELxJCSvhNX/O8PM464xSG334351xwKUOHHMQmPXuQm5vLvOJivHd8fpg1m8suOo8PP57MitJSDth/31rPd8D++zJpylTGf/gRkUizhMrbB+7PGaedxO133l9vLLGD9n7++VdKSyu7mNasXUNJyZ+89fZY5s6dx0Xnn72uXPv27bj2qku54abbOP2s83EOG8oWm29KhAg/zp7D62++xcpVq7j+mssTmgIYJkr4IiGkp+Vlh4F77UHXLl0Y/dbbvP7GGBYuWkRZWTnt27ej7zZbcfopJ7H1Vltw3kWX07p1K/bcfZdaz7XLgJ1o374d74x9n0gkklB5+8D92WuP3ejRvUe9scS64ebb1v27WbNmbLRRBzbt3YsRt/5ng/EB223bl4cfuIvRY95mrP8Bzz3/IkQidO7UkX0G780hQw6ifbt2Da7PbBepaOgE0DQrKZ5bUdSl9gUYpPFKiueiOg5eOur51VfhiKjbnIcdBq6b0hACsWDBAjp27LjB9rKy7LuHn2lUx41T22c3RqS+AnXRSnsiIaQufZHwUcIXCSElfJHwUcIXCSElfJHwUcIXCSElfJHwUcIXCSElfJHwUcIXCSElfJHwSes8fMt29gbG1rArD3jG99yTatgnIo2khC8SPmlN+L7nfghUWwLEsp1uwHTgqfRFJpLdlPBFwicTu/QfBV72PXdCugMRyVY1PTynia7BJSJxyqildS3bcYCdgePSHYtINmvWDJo3h9Wr129btWrDJXdFJHtkTMK3bCcXGA7c4Hvu4vrKl5eXUVI8NzXBhZTqODXSVc/5+V1YvXp9J99vP82jXbum3cxfubqcsrLVG+6oqKh5uySP6rhRVpaW1vt3oLFLcGdMwgf+DhQBT8RTODc3T+u8B0xr6adGuuq5oACWLFn/vmVhV4q6pDyMpFqwYEGN67lrnffgqY4bJ7+ggKL619JvlExK+McDL/meuyLdgYiEgQbuZYf58xfwivsmn33+JfPnL6BZs2Z06lTErjsP4JChB9GxqChp17rosqv544/5PPfUo0k7Z31uu/Ne/PfG4XupebrTM8+N4tmRL/Lsk4/QpXOnlFwzVTIi4Vu20wrYB/i/dMciEhZK+E3f1E8/48bhI2jXti2HHXIwm27aG4Aff5yD++ZbvPXOWK6+/GJ23KF/wuee+f0szj7/kmqJ9vxzzqS8rDypv4OkTkYkfKCfmZ73VboDEQkLJfymbd68Ym4cPoLevXox/MZrKYj6D9p/234cuP9+XHbVMG68ZQSPPng3RUUbJXT+GV9/s8G2nj26JyV2SY9MSfhVn6IFaY5DJDSU8Ju2US+/xqpVq7nkwnOqJfsqLVsWcOlF53Ly6efw8mtvcOZp/wLgvIsuZ+XKVVx8wdnc/9B/+XH2HFq0aMHAvfbgjFNPonnz5lx02dVM/9/XUNkDy3bb9uWOW2/coEv/ltvu4ouvpnPnbTdxx93388OsH2ndujX/+PthOIcO4YWXXmX0mHdYsWIFW2+9JReedxadOq6/xTD5k0956RWXn3/+lTVr1tC9ezeOOsJh4N57xF0PTz79PCNffIXHHr6XTTbuWW3fTbfewZQpn/LSyCcpKChgxtff8vyol5k160dWrlxFly6dGWIfwKFD7VrPX9XF//gj97Fxzx4bbI/u+l+4cBFPPP0cU6d9ztKlyygq2ojBA/fkuGOOonleXty/U1AyIuH7nvsaEEl3HCJhooTftE2ZOo2tttyCHt1rH/DZs0d3ttpycz76eNK6hJ+Tk8PCRYsYcdf9HH7YULp07sTHk6bw+ptvAXDuWadz/jln8uhjTzFl6jTuv/t2Wrbc8AtF5bmaUV5ezoi77mP//Qbzz2OO4unnRvHQo08we85PlJev4dKLzmXOTz/zyGNPcf9Dj3LDtVcC8OVX/2PYf4azy847cdXlFxGJRHDfGMONw0fQsmUBA/62Y1z1MHjQXox88RU+njSlWsJfvXo1n0ydxq67DqCgoIBffv2Ny68exmab9uHiC88lv0UL/PfHcf9D/yUnJ4ch9gEJ1X+sFStKueCSKykrK+NfJxxH166dmfH1t4wc9Qq//z6Pa668pFHnT4aMSPgiknpK+E3X8uXLWbhwEbvvunO9Zfv07sW3333P6tWrad68OZFIhMWL/+K8s89gz913BaD/dv345Zff8N8bxxmnnkTPHt1pU9gGgC232Kz2k0ciLF26jCH2AQzca/fKUfqRCJddeR3ffDeTxx++j0gkwg7bb8dHEyczY8a36w5dUFLCgJ124NILz6N161YAbLH5Zhx+1D8ZN+GjuBN+r002ptcmGzNx0hSO/b8j122fOu1zSktXsu+gvQGYV/wHO2zfn1P/dfy6lnq/vlvzydRpfDD+w0Yn/Dffepu584p58N4RbL7ZpmBurVABTz07kjlzfqZ3700adY3GUsIXCakwJvxh44dx/YTr4yp76o6n8ujQ6qPRTxt9Gv/9/L9xHX/dwOsYNmhYtW1DXxjKmO/H1FkmHqUrVwFQUFD/SklV3f1Lly1jow4dwLTyd45JqP379+PzL7/it9/n0qd3r4Ti+dtOO6z7d1WX/fbbbUskEona3pGvv/lu3Xtr38FY+w6udp7WrVtR2KYNCxaUJHT9QQP35KlnRlL8x/x13esffjSJwsI262LbZcBO7DJgp2rH5eTk0KlTJxaU/JnQ9Wry2edf0qFD+3XJvsoeu+/CU8+O5Mvp/1PCF5H0CGPCzxaFpvW9dOnyessuX76cSCRC28LCddvaFhbSvHn1OfPt2rYF4K+/lmxwjvq0a9t23aI7ubm51c5XJSc3h4qo9ZtLS0t5+rlRTJr8CSUlf1JWvn70/9oE13kebBL+xxMnc8Thh67rzt93n4Hr4lmzZg0jX3yF8RM+pviP+ayOWmayc6fGz38vKfmThQsXYdlOjfv/XLiw0ddoLCV8kZBSwm+6mufl0aVLZ2b9OLvesnPm/Ey3rl3WJT6AZs1qeoxKZZKNbpUH6c57HuTDjydxxOGHsvPfdqRVq5ZAhMuuui7hc3Xr2pUtt9iMjydN4YjDD+XTaV+worSUfQbtta7Mk888z4svuxx0wH78+/STKSxsQyTSjNvuuIcVKxq//EskEqFbt65cc0XN9+rbtS2scXsqKeGLhFRND9DJdsMGDWtQF3qVR4c+ukE3fyJGHz26wcfG2muPXXn51TeY9eNsNtu0T41lfv3td2b+MItjjjqi2va//vqLNWvWkJOTs27bosV/AdC+fbukxVib1WVlfDRxMnvtsRun/uv4ddvLysooLW3YB3HQ3nvy3yeeYdGixXw0cRKdO3Wk7zZbr9v//gcfstUWm3PheWdVO660nm+6kWaVX4AqYnodFi2uvgJ8x6Iifpj1I5v26ZWyL02JysSn5YlICqiF37QdefhhtC0s5PY772Pp0mUb7C9duZIRd91Phw7tOdwZWm1fWXk5n39RfdmTr6bPoHXrVutG/UfMxKk1a9YkPfbVq1azZs2adbcmqrz19ljKyspYu3ZtwucctPeeYGYvTJk6jcGD9q6WeEtXlq4biFhl8iefUvzH/Dqv16Z1awD+/HN9l/zatWv5avqMauV23KE/S5ct45Op06ptn/n9Dzzw8GMbfEFIB7XwRUJKCb9pa9++HddedSk33HQbp591Ps5hQ9li802JEOHH2XN4/c23WLlqFddfczmFbaonujZtWvPgI49zZEkJ3bp2ZeLkT/jiy+n835GHr2v1d+jQHsx8/14bb8weu++StNhbt25Fnz69GD/hY7bttw0dOrTnk6nT+PmX39hm6y2Z89PPfPb5l9Va6PUpKtqIvttsxQsvvsLy5SvYd/De1fb337YfU6d9jvfOWHr26MH0GV8z5ZNp7LXHbkyc/AkTJ33Cdtv23eC8/bftRyQS4cWXXTp1LGLN2rW8/Orr5Oe3qFZuyMEHMObtdxk+4m5OOv5Y+vTuxa+//cbTz42ifbt2G/w3SAclfJGQUsJv+rbbti8PP3AXo8e8zVj/A557/kWIROjcqSP7DN6bQ4YcRPt2G3bRtywo4LKLz+OhR59k1o+zKcjP5++HDeXE449ZV2bIQfszddpnjHrxVTbbrE9SEz7AlZdeyL0PPMKd9zxIfn4L9thtF66+/CKmff4ld93zICPuup+7R9yc0DkHD9yLex94hD69e9Frk42r7Tvn36dxz/0P88hjT5PTrBk77bg9N1x7BXPnFfP1t99x9/0Pcf01l29wzt69N+GSC8/lhRdf4bSzLqBD+3YcMuQgttl6K+6690HKzWDDVi1bctftN/PUM8/z/KiXWbp0Ge3bt2PP3XflxH8eU+32SbpEYu9LNBUlxXMr9CS3YOlpeamRrnp+4AE4++z17884Ax56KOVhJNWCBQvoWMMTx/Qkt/WCegCO6rhxavvsxmjU4ADdwxcJKbXwRcJFCV8kpJTwRcJFCV8kpJTwRcJFg/ZEQkoJP5zuuPXGdIcgaaIWvkhIKeGLhIsSvkhIKeGLhIsSvkhIKeGLhIsSvkhIZWPCr6io2GDNcxGppIQvElLZ+PCcFi1asGrVqnSHIZKQ1atXV3uaYVCU8EVCKhtb+IWFhSxevJiVK1eqpS9Nwtq1a1m4cCHtalgCOdk0LU8kpLIx4UciETp37sySJUtYsmTJuqelrSwtJT/2F5akUh03XIcOHVLySF0lfJGQalH9YV+sXg1r1kAGPOOjUSKRCG3btq22raR4LkX1r1MujaA6znzq0hcJqUgE8vOrb8uG+/giUjMlfJEQy8ZufRGpmRK+SIgp4YuEhxK+SIgp4YuEhxK+SIgp4YuEhxK+SIgp4YuEhxK+SIgp4YuEhxK+SIgp4YuEhxK+SIgp4YuER0astGfZzr+Ay4BNgLnAfb7n3pXuuESyXTY+QEdEapb2Fr5lO0cB1wLHA22Bc4AzLNsZkO7YRLKdWvgi4ZEJLfzrgEt8z/3EvH/L/IhIwJTwRcIjrQnfsp2uwNZAjmU704AtgdnATb7nvpTO2ETCQAlfJDzS3cLf2LyeARxt7t+fArxo2c4fvudOqO3A8vIySornpi7SEFIdp0Y667liTRugzbr3fy5YQknxsrTEEiR9loOnOg5eUZdujTo+3Qm/6vrX+577g/n3PZbtHAucBNSa8HNz8xr9y0vdSornqo5TIJ313KGo+vtITiFFXQrTEkuQ9FkOnuo486V70F6JeV0Us3020DUN8YiEirr0RcIj3Ql/lkn6u8Rs3wyYk6aYREJDCV8kPNLape977hrLdu4ErrNs5wtgOnAqsIN5FZEAKeGLhEe67+EDDDc9Da8AGwHfAUN9z/0i3YGJZDslfJHwSHvC9z23ArjJ/IhICinhi4RHuu/hi0gaKeGLhIcSvkiIKeGLhIcSvkiI6eE5IuFR6z18y3b2bsD5IkAL33PHNi4sEUkFtfBFwqOuQXvjgYoGnjengceJSAop4YuER10J/4YaEv5AYFfgfTN9rtQsxN3X7HsXeDvgmEUkSZTwRcKj1oTve+6w6PeW7RxiFsPp63vuBqvgWbazOeADjwcVrIgklxK+SHgkMmjvGuCBmpI9lV8QfgDuBa5NXngiEiQlfJHwSCTh9wX+qKfMQvN8exFpAnJzoVnUX4E1a6CsLJ0RiUhQEkn4S4FD6ylzCJB9D9MWyVKRiFr5ImGRyNK6rwJnWLbjAy8CP5pBe/lAH+AI4ADgqQDjFZEkKyiA5cvXvy8thcLCdEYkIkFIJOFfCvQAhgD71LA/YqbyXZTE+EQkYGrhi4RD3Anf97Uo4mYAACAASURBVNxlwCGW7fQ3Cb830NK08n8BJvieOzXYcEUk2ZTwRcIh4afl+Z77FfBVMOGISKop4YuEQ8IJ37KdjYHuQF5tZXzP/bDRkYlISijhi4RD3Anfsp0ewOvADnEU19K6Ik2EHqAjEg6JtPDvAnYExgJfmnv3ItLEqYUvEg6JJPy9gPt8zz0vwHhEJMWU8EXCIZGFd1oBnwQYi4ikgRK+SDgkkvA/A7YIMBYRSQMlfJFwSCThXwacbtnOrgHGIyIppoQvEg6J3MM/GJgITLRsZxowC1hVQ7kK33NPTmKMIhKg/Pzq75XwRbJTIgn/6qh/DzA/NakAlPBFmgi18EXCIZGEv7Vp0VcEGI+IpJgSvkg4JLKW/sxgQxGRdFDCFwmHRFba2ziOYhGgue+5PzQuLBFJFSV8kXBIpEv/pwS687W0rkgToYQvEg6JJPxnakn4zc2jcncG3gX0iFyRJkQJXyQcErmHf2Jd+y3b2QYYDdyalMhEJCX08ByRcEhk4Z06+Z77DXA9cHOyzikiwVMLXyQcEunSj8d3QP9EDrBs5yegO7AmZtd2vud+n9zwRCSWEr5IOCQ74e8CLG/Acaf6nvtUkmMRkTgo4YuEQyLT8q6tY3dzYEvgEGBMckITkVRQwhcJh0Ra+MPiKPM/4IIGxHGkZTuXAd2AH4Bhvufqi4NICijhi4RDIgl/cB37yoB5vufOaUAM04EfgdOApcD5wJuW7ezhe+7k2g4qLy+jpHhuAy4n8VIdp0a663nFkmZAl/XvV6ylpLg4bfEEId11HAaq4+AVdenWqOMjFRWZtzS+ZTufAdN9zz2ptjIlxXMrGvvLS91Kiuc2+gMm9Ut3PS9bBm3arH9fUAArVqQtnECku47DQHWcEpHGHJzwoD3LdnYGhgCbAa1Nq/xb4DUzNS8ZZgFdk3QuEalDTV36FRUQadSfFhHJNIkM2ssBngaOruVbxvWW7dzhe+6lCZyzN3ApcIXvuYujdvUFxsV7HhFpuJwcyMuDsrL121atgvz8dEYlIsmWSAv/AuAYs3zu48BMMwWvlUnQpwEXWbbzTQJT7IqBoUChZTvnAKuBi4HNgb838HcSkQQVFFRP+KWlSvgi2SaRhH8c8J7vuQfVsO9/wCjLdnzgDCCuhO97bqllO/uZ5Xh/MD0H04FBehyvSOoUFMCSJevfl5ZC+/bpjEhEki2RhL858HA9ZV4A7kkkAN9zvwMOTeQYEUkuracvkv0SWUu/AlhVT5nyZK7PLyKpobn4ItkvkeT8PXBEPWX+YdbTF5EmRAlfJPsl0qX/FHC3ZTtjgEeBr82gvdZm0N7pwAHA2QHGKyIBUMIXyX6JJPz7zJPwTgJqGrgXAR7wPfehJMYnIimghC+S/eJO+L7nVgAnW7bzCOAAW5nW/TLT2n/F99wvgw1XRIIQm/CHDKmcnx+kzp3hppvg2GODvY6IVKo14Vu2cyTwue+5P5r3xwPjfM+dCkxNaZQiEqjYhF9WVn1efhB+/hlOOQUOOwxatQr2WiJS96C9p4Bdo94/CQxIQUwikmKbbZae665cCXMa8sgtEUlYXV36fwFXWLbT1nTbR4BBlu20ru+kvuc+k9wwRSRI550HH31U+ZNq0Qv+iEhw6kr49wC3APebOfgVwFn1nC9iyinhizQhnTrBhx9WrqG/dm2w1zr4YBgX9aSMpUuDvZ6IVKo14fuee6tlO28C/cx6+U8ADwDTUhuiiKRKixbBX6Ndu+rvlfBFUqPOUfq+535rHn2LZTvXAWN8z303ZdGJSNZp06b6e3Xpi6RGItPyegcbioiEQWzCVwtfJDW07r2IpFRhYfX3SvgiqaGELyIppS59kfRQwheRlFKXvkh6KOGLSEop4YukhxK+iKSU7uGLpEciT8uDyul5RcC+wMbAKN9zfzXbW/meuzyQKEUka+gevkh6JNTCt2znauBXYCQwHOhttucAMy3buSawSEUkK6hLXyQ94k74lu0cA9wAfAlcbpbRrdLabB9myomI1Ehd+iLpkUgL/yxgErA78Gj0Dt9z//I9dwgwIY719kUkxNSlL5IeiST8fsBI33Mr6ijzArB1EuISkSylLn2R9Egk4ecB9Q3KKwNS8PgNEWmqWsc8YHv5clizJl3RiIRHIgl/FrBXPWWOMeVERGqUkwOtWlXftmxZuqIRCY9EpuU9B9xk2c4MYIzZFrFspyOwC3A+MBi4NKBYRSRLtGlT2bKvsnQptG2bzohEsl8iCf8OYFfgTvNvgPejRutHABe4K4A4RSSLFBZCcfH697qPLxK8RB6PuwY43LKdw4B/ANsArYClwNfAS77njg42XBHJBhq4J5J6Ca+053vu68DrwYQjImGgqXkiqae19EUk5dTCF0m9uFv4lu3MjqNYBMjzPbdH48ISkWym1fZEUi+RLv1mQE2L7jQHuphkPx1Y1NBgLNvZA/gQ+I/vucMaeh4RyWzq0hdJvUQG7fWqbZ9lOy2Ak4BzgeMbEohlOwXAk4Bm5IpkOXXpi6RewoP2auJ77irgYct22pppeUc04DQ3A98Bc5MRk4hkLnXpi6ResgftfWgW30mIZTt7Av8ETk9yPCKSgdSlL5J6SWnhR+kK5CRygGU7LU1X/vm+586zbCeu48rLyygpVmdAkFTHqRHKel5bALRf93bB/BWUFC8O7HKhrOMUUx0Hr6hLt0Ydn8go/b3r2N0c2BK43CzCk4ibgW98z30ukYNyc/Ma/ctL3UqK56qOUyCM9dy9Z/X3ZeUtKerSMrDrhbGOU011nPkSaeGPr2WUfpUIsNI8QCcupiv/GGC7BOIQkSZOg/ZEUi+RhH9DHQm/DJgHvON77rwEznkyUAjMiOrKbwvsbNnOIb7n7pjAuUSkidA9fJHUS2RaXhDz4i8EronZ9jIwGbgtgOuJSAZQC18k9ZI9aC8hvucuil2ox7KdVcAS33OLaz9SRJoyTcsTSb1aE36cS+nWpML33E0bGpDvuYMaeqyINA3q0hdJvbpa+LUtpVufSCPiEZEQiE34y5ZBRQVE9NdDJDC1Jvy6ltIVEWmMvDzIz4eVKyvfr10LK1ZAq1bpjkwkeyV1pT3LdvaybOfuZJ5TRLKTBu6JpFbCg/Ys28kFOtVwbIGZZvcP4PzkhSgi2ahNG1iwYP37JUugS5d0RiSS3RJZaa+ZmSp3hknuNal6RK6ISJ00Ul8ktRLp0j/dzJtfCEwwyX0a8DGwAlhslsk9LMB4RSRLqEtfJLUSSfinAm8DvYDDzbZLfM8daB6a8wrwNz3eVkTioal5IqmVSMLvA7zse+7aqOl6ESpH9C/zPfd0c4vgP8GEKiLZRF36IqmVSMKvANaaf68yrzH/yzIKOCpJsYlIFlOXvkhqJZLwfwb2p7JFvxL4C9grpkwLoCi5IYpINlKXvkhqJTIt70XgJst2yn3PPREYB5xr2c4iM3CvF3Al8EOA8YpIllALXyS1Ekn4twN9o6bkXQ0MAm407yPAGuDMAOIUkSyje/giqZXI43HLgeMs28kz77+xbKcfcAKwCTDPDOr7JtCIRSQrqIUvklqJLLwz0PfcCb7nllVt8z13LnBLYNGJSNbSPXyR1EqkS3+cZTu/ACOB533P/TrAuEQky6lLXyS1Ekn4jwOHApcDl1m2Mx14FnjB99x5AcYoIllIXfoiqRX3tDzfc081K+odCDwBdAdGAL9YtjPWsp3jLdtpHWy4IpIt1KUvkloJPR7X99w1vueOjUn+TwHbm9c/LNsZGVy4IpIt1KUvkloJJfxoMcm/C3AasEQr7YlIPNSlL5JaidzD34BlO7sAfzdPyNvUbJ6cnNBEJJvV1KVfUQGRSLoiEsluCSd8y3b2Ao4wT8zrZhbcmW5W2XvB99xfgglVRLJJixaQmwvl5ZXvy8th1SrIz093ZCLZKZF5+A+Zlnwnk+RnAzcDI33P/TbYMEUk20QilffxFy5cv23pUiV8kaAk0sI/HSgG7jNJfmqAcYlICLRps2HC79gxnRGJZK9EEr4FfOB7bkWA8YhIiGhqnkjqJLKW/vvBhiIiYaOpeSKp0+BpeSIijaWpeSKpo4QvImmjLn2R1FHCF5G0UZe+SOoo4YtI2qhLXyR14k74lu08YdnOjvWUOdmyneeTEpmIZD116YukTiLT8k4ERgOf11GmB3BQIgFYtrOtWcBnN6A5MBO4yffc1xM5j4g0PerSF0mdOhO+ZTsnACdEbbrBsp1zaimeD+wELIj34uZxuuOAZ4FjgNXAxcArlu1s53vuN3H/JiLS5KhLXyR16mvhzwJ+AQYAFUDfOsquBeYAFyVw/QLgcrNy3woqvwTcB9wI9AOU8EWymBK+SOpEKiriWzjPsp21wBG+574WVDCW7RQBVwD/AHbwPbektrLFv/1ckZubF1QoApSXl6E6Dl6Y6/mdd1vwzxM3Wvd+v31X8sJzC+s8piHCXMepojoOXlGXbo16lmQi9/AHAzMac7G6WLazytzD/wyw6kr2ALm5eRR16RZUOAKUFM9VHadAmOu5xybV369anR9IXYS5jlNFdZz54h6l73vuBKDAsp39o7dbtnO0ZTuvWLbzgmU7AxsaiO+5LYCOgAdMsmxni4aeS0SaBnXpi6ROItPy+gJfRd+jt2znXOA54HDgKMC3bGdAQ4PxPbfE99xrgT+AMxp6HhFpGjQtTyR1Ell452pgOXAelck+B7gSmAtsBfQ2g/wuifeElu0cbNnOr5bttIrZFQHKE4hNRJogTcsTSZ1E7uHvCdzve+535v0eQCfgUt9zv6cygf8XODeBc04x9+3vs2znIqAUOB3YDAhscKCIZAZ16YukTiIt/CLg56j3+5mpel7UtnlAl3hP6Hvun4Blvjh8b7ryjwEc33OnJBCbiDRBLVtCs6i/QitXQllZOiMSyV6JtPD/NIPqqhwE/BazOE4HIKG7cL7nTgeGJHKMiGSHSKSylf/XX+u3LV0KHTqkMyqR7JRIwv8SON2ynYnAQGBH4O6qnZbtRMz8+e/qPo2IyHpK+CKpkUjCvx14D/jUDKr7E7gjav+7wF7ASQHEKSJZSvfxRVIj0Xn4ewP3AXcCu/ueOzeqSBlwpe+5zwQTqohko9iR+pqaJxKMRFr4+J47GZhcy76DkxaViISGWvgiqZFQwmf9evf7AhsDo3zP/dVsb+V77vJAohSRrKWEL5IaiUzLw7Kdq4FfgReA4WaxnapFeGZatnNNYJGKSFbSansiqZHI0rrHADeY0fqXmYF7VVqb7cNMORGRuGi1PZHUSKSFfxYwCdgdeDR6h++5f/meOwSYYMqJiMRFXfoiqZFIwu8HjPQ9t6KOMi8AWychLhEJCXXpi6RGIgk/zzw8py5lQItGxiQiIaIufZHUSCThzzIL69TlGFNORCQu6tIXSY1EpuU9B9xk2c4MYIzZFrFspyOwC3A+MBi4NKBYRSQLKeGLpEYiCf8OYFezyl7VkrrvR43WjwAucFcAcYpIltJKeyKpEXfC9z13DXC4ZTuHmYfkbAO0ApYCXwMv+Z47OthwRSTbqIUvkhoJr7Tne+7rwOvBhCMiYaOEL5IaCSf8aJbtNAO2AAqAGb7nliUvNBEJA03LE0mNehO+ZTv5wCFAR+AL33Mnme07A88DfUzRhZbtnOV77kuBRy0iWUPT8kRSo86Eb9lOe+BjYCszKK/Csp0nzIj8V0zL/lWgOTAIeN6yndm+505L3a8gIk1Z69bV3y9fDmvWQE5OuiISyU71tfAvNivnvQBMBPoDJwOrgLVAX99z51P55aAnMA240MzHFxGpV7Nm0KpVZaKvsmwZtG2bzqhEsk99CX8I8IbvucdWbbBs5xfgeuDKqmRP5WC+Xy3beQw4PtCIRSTrFBZWT/hLlyrhiyRbfSvt9QHejdn2ijnu2xrKzwQ6JTE+EQkBjdQXCV59Cb8VsCBm2x/m9a8ayi9v7Mh/EQkfJXyR4MWzln7s0/HqelqeiEjCtNqeSPASeXiOiEgg1MIXCV483e+HWbazWdT7fNPKP8aynV1jym6b5PhEJASU8EWCF0/CP66W7afVsl1d/iKSEHXpiwSvvoR/UoriEJEQUwtfJHh1Jnzfc59OXSgiElZK+CLB06A9EUk7PUBHJHhpnzNv2U5n4FbgQLM2/9dmFb/x6Y5NRFJDD9ARCV4mtPDfADoD25vXCcAYy3a6pTswEUkNdemLBC+tCd+ynULTor/A99xi33NXmtZ+KyB2yp+IZCklfJHgpbVL3/fcJebpe9H6mNe5aQhJRNJA0/JEgpf2e/jRTIv/SeAt33On1FW2vLyMkmJ9JwiS6jg1VM9Qviq32nO3Fi8qo6Q49jEejTi/6jhwquPgFXVp3J3ujEn4lu1sAowB5gNH11c+Nzev0b+81K2keK7qOAVUz7DxqurvV5Qm9/9v1XHwVMeZLxMG7WHZzgBgKjAROND3XN3BEwkRdemLBC/tLXzLdvoB7wA3+J57T7rjEZHUix20t2wZVFRAJJKuiESyT7pH6ecATwMPKtmLhFduLuTnr3+/di2sWJHOiESyT7pb+LsBOwL9LNu5JGbfs77nnpqmuEQkxdq0gZUr179fuhRatUpnRCLZJd3T8j4G1GknIhQWwoKogflLlkCXLumMSCS7ZMSgPRERLb4jEiwlfBHJCEr4IsFSwheRjKCpeSLBUsIXkYygFr5IsJTwRSQjKOGLBEsJX0Qygrr0RYKlhC8iGUEtfJFgKeGLSEZQwhcJlhK+iGSE2C59JXyR5FLCF5GMENvC1z18keRSwheRjKAufZFgKeGLSEZQwhcJlhK+iGQETcsTCZYSvohkBLXwRYKlhC8iGUEJXyRYSvgikhFq6tKvqEhXNCLZRwlfRDJCixaQl7f+fXk5rFqVzohEsosSvohkDHXriwRHCV9EMoZW2xMJjhK+iGQMrbYnEhwlfBHJGOrSFwmOEr6IZAwlfJHgKOGLSMbQansiwVHCF5GMoRa+SHCU8EUkYyjhiwRHCV9EMoa69EWCo4QvIhlDLXyR4Cjhi0jGUMIXCY4SvohkDK20JxIcJXwRyRhaaU8kOLnpDgDAsp3ewJPAQKC377k/pTsmEUk9demLBCftLXzLdhxgCvBzumMRkfRSl75IcNKe8IEOwN7As+kORETSS136IsFJe5e+77mPU9nS75nuWEQkvdSlLxKctCf8hiovL6OkeG66w8hqquPUUD2vt2pFBOi67v3SpRWUFM9r9HlVx8FTHQevqEu3Rh3fZBN+bm5eo395qVtJ8VzVcQqontfbqAKaNYO1ayvfr1wZoe1G3cjLa9x5VcfBUx1nvky4hy8iAkAkom59kaAo4YtIRlHCFwmGEr6IZBRNzRMJRtrv4Vu2MxPYJOrLx0zLdiqAZ33PPTXN4YlIimlqnkgw0p7wfc/dMt0xiEjmUJe+SDDUpS8iGUVd+iLBUMIXkYyiLn2RYCjhi0hGUZe+SDCU8EUko6hLXyQYSvgiklHUpS8SDCV8Ecko6tIXCYYSvohkFCV8kWAo4YtIRtE9fJFgKOGLSEbRPXyRYCjhi0hGUZe+SDCU8EUko6hLXyQYSvgiklHUpS8SjLQ/PEdEJFpNCX/OnMadc+GCHJaW1l0mNxcKCip/8vMhJ6dx1xTJNEr4IpJRWreu/n7FCujTp7Fn7ZzwEXl5678AVH0JyE3RX8xmzWCbbWCffSp/+vSBSCQ115bspYQvIhmlWbPKpL9sWXrjKCur/EnXLYXp02HUqMp/b7zx+uQ/eDD06JGemKRpU8IXkYyzww7w0UfpjiJz/PILPPVU5Q/AFlvAoEGw0Ubpjmy9Fcva0LJ1HAVDKhKp/DJb388OO1R+sQuCEr6IZJzHH4dzzoGZM5NzvrVrymmWU/ufu4qKytZ8aSmsXFn5msm+/77yJ7O0iaOM1Oecc4JL+JGKiopgzhywS8dcVHH7Z3fGVfbUHU/l0aGPVtt22ujT+O/n/43r+OsGXsewQcOqbRv6wlDGfD8mruMfGfIIp+10WrVtOz26E5/P+zyu49/8vzcZuuXQatu63dGNecvmxXX8tFOnsVO3napti1wf/w3B3y/8nW5tuq17P3fpXLrf2T3u4yuuq/4Z+2zuZ/ztv3+L69iurbsy96K51baNnjmaQ0YdEtfxO3bdkc9O+6zatkc/e5TTx5we1/FDthjC6KNHV9s2bPwwrp9wfVzH67OXGZ+9kuK5FHXpltBnr0vrrnx78lxKS1n388Hvozl3cnyfva3b7ciL+1T/7L0y51Fu+CK+z96WkSF0nzCaiRNh1SqzcdAwGBTfZ4/PToXR1T97DD0Ndorvs8f462B89c8eRw+FLeP77DH6Efis+meP03aCbvF99hj5Jnxf/bPHRd2gTXyfPR6ZBvOqf/YYlsBAiDt+h6XrP3u0mQsXxf/ZY1hMbu36GZzeuL97Q7cc2qiRHGrhi4jUIAK0a1f5U2V2LjA5vuMLCmDbbatvm7wa+CK+4zffHEZfW9njMHkyfPABPPUz/JbA7yASTQlfRCSD5edXDtQbPBhyxsP1E+I77m9/A2fX6tvccpgWZ6fuPvvAvlb1bc+Uw8w4jz/0MNj58OrbHiiDubUdEOO4f8LWMSvFDC+DeNdh+vdZ0D2mPXxVWZwHA5ddBoVRxy+pgFvL4z/+ppuqv/9tLTy0Jv7jg9Bku/RLiudWFHXpFkdJaaiqblAJluo5eKrj4KmOU6JRXfpaaU9ERCQElPBFRERCQAlfREQkBJTwRUREQkAJX0REJASU8EVEREJACV9ERCQElPBFRERCQAlfREQkBNK+tK5lOwXA7cCRQCHwNXCp77kfpDs2ERGRbJEJLfwHgf2AwUAR8ALwlmU7m6c7MBERkWyR1ha+ZTsdgOOAI3zP/cZsvsOynWOBM4CL0hmfiIhItkh3C39H86Xj05jtU4FdazlGREREEpTue/idzOufMdtLgC51HVjUpVujnhok9dOTr1JD9Rw81XHwVMeZL90t/NpEgKb53F4REZEMlO6EX2xeO8Zs7xS1T0RERBop3Qn/M2A1sFvM9t2BiWmKSUREJOtEKirS23Nu2c6DwP7AocDPZmT+JUBf33N/TmtwIiIiWSLdLXyAC4AxwDhgAXAgsL+SvYiISPKkvYUvIiIiwUv3tLyEaBneYFi20xt4EhgI9PY996eofarzJLBspzNwq+nBKjD1eKXvueNRPSeFZTvbAjebMUHNgZnATb7nvo7qOOks29kD+BD4j++5w1AdJ4VlOz8B3YE1Mbu28z33+8bUcSZ06SdCy/AmmWU7DjDFjJ+oieo8Od4AOgPbm9cJwBjLdqomL6ueG8GyndbmtuAsoLeZ+fM68IplO9uYYqrjJDFJ50lgWcwu1XFynOp7bn7Mz/dmX4PruMm08LUMb2A6AHsDPYHjo3eozpPDsp2qb+G3+55bbLbdClwO7GrZznjVc6MVmPoc6XvuCirr+D7gRqCfZTvFquOkuhn4DphbtUF/L4LX2DpuMglfy/AGw/fcx6n8IPWsYbfqPAl8z10CnByzuY95nat6bjzfcxcAj1W9t2ynCLgC+A34QHWcPJbt7An8E9jWtC6rqI6T50jLdi4DugE/AMN8zx3T2DpuSl36DV6GVxpMdR4A0+J/EnjL99wpqufksmxnlZnxMxCwfM8tUR0nh2U7Lc1n93zfc+fF7FYdJ8d04HvTbd8TeBN407Kd3Rpbx02phV8bLcObeqrzBrJsZxMzDXU+cHQ9xVXPDeB7bgvTwj8XmGTZTl0tH9VxYm4GvvE997kEjlEdJ8D33ENiNt1g2c6hwGmAX8thcdVxU2rhaxne1FOdJ5FlOwNM19tE4EDfc5eaXarnJPM9t8T33GuBP8y9TdVxI5mu/GOA02spojoOziyga2PruCm18KOX4X05avvupsUkyac6TxLLdvoB7wA3+J57T8xu1XMjWbZzMPAwsJXvucujdkWActVxUpxspoHNsGynaltbYGfLdg4xo8ZVx41gpkhfClzhe+7iqF19zSyURn2Om0zC9z33L8t2HgdusWznm6hleDcx0xQkyVTnyWHZTg7wNPBgDcle9ZwcU8zc+/ss27kIKDUt0c2A11THSXEhcE3MtpeBycBtquOkKAaGAoWW7ZxjkvvFwObA3xtbx00m4RsXmMVLxgFtgC+1DG/jWLYz03xYqm7vzLRspwJ41vfcU1XnSbGbGV3bz7KdS2L2qZ6TwPfcPy3bscw95u9N8v8OcMzASFTHjeN77iJgUfQ2M0BySdV0U9Vx4/ieW2rZzn6mDn8wPVTTgUG+5840xRpcx1paV0REJASa0qA9ERERaSAlfBERkRBQwhcREQkBJXwREZEQUMIXEREJASV8ERGREGhq8/BFMpplO8OA64DBvueOT3c8ibBs527gWPOoWcv33Mnpjqkm6a5jy3YGmTnQ1/ueOyzV1xdpKCV8aRIs2znRPKVrJdDX99zZtZT7CfjJ99xBqY+y6bJs5yDgPLPO/8PAT+mOKRNYttParHR2d9RSp18DRwLf1HO4SEZRwpemJh+4H7DTHUiW2da83up77ug0x5JJBpjehKeAxVSuhrYAeCXdgYkkSvfwpakZBxxk2c7h6Q4k3SzbyU/i6Vqa16X1lAubAekOQCRZ1MKXpuZWoCdwt2U778Y8GW0DUbcCTvI996mYfeOBgb7nRsz7A4G3gVOBZcC1QB9grnk4yMOmzI3ANuZBF4+aVnHsGtURs27+6Sbe+cCzwHW+55ZFxZBrHkryT/OglzJgBvCA77nPR5Wrum98jXn62yVmje3B9fz+B5qyO5l7878Db5mn9i0wZaJjH2eehFbn/XHLdvqb+tnbPDGtGHgXGOZ77u+W7UTMgz3aAp18z10Vc/wewMfAM77nnmC27WceBLI90N7EOg641vfcuXXE0guYAzzte+6JMftqvN9v2c4ZwEmmzvPNQQyLNgAACRxJREFULYyRwJ2+55ay/vbQJuaQOZbt4HtupLZ7+JbtbFV1LWAjYCHwkanr6VHlRgFHmS9Z1wLHmfJzzPUfjyqbB5wNHA/0Nn+zfzEPrbk5tl5F6qIWvjQ1q8wfwJ5AsgdMVf3x3Nec+yGTLJsDD1m2c4FJ8C8A55snst1i/njHuhj4BzDClJ0LXGHeR3vJnGMa8G/gSrP9Oct2rq3hvDsD/zKJIvZc1Vi2cwLgAd1MIjrDJKqzgYnm/jTmfnTVozaHmfdf13He3cwT0vqbL2CnmN/jaOBTy3Y6my9Ao8zjVK0aTnNk1e9pznmAeXxwL1MfJ5uYjgUmR8XaaJbt3Gj+2y42jyL9t3kAyY3A81FFzwSqviT8Oyrmms7ZD/gEONCc+yTgv8ABJv4doopXfc4eNT0I/zFf+nKAx8yjZqvcC9xpHgR0gYnjQ+Bq8zkUiZta+NLk+J77rmU7rwLnW7bztO+5M5J06qqW7lCgj++586n8Y77A/HEdAfTzPfdbs326SXwHm+QWbXOgf1Rr8SnTcj/Tsp1hvucusmznMMABLvY9946qAy3becgMnrvKsp2Hq+IwbGBL33N/rOsXsWynALjL9CzsFjXg7EnLdmab5HY2MNz33FdMwgKYEMfI94eAP4GdfM/9K+qa44HR5ovN+abOLgH+Hv2sbtP6/zswD/jAbN7SfBk53/fcr6PKAlwGHFb15SAJOpuenEN8zy032562bKcP4Fi208P33N98z33bsp2qL3Nv+55b10DG282Xm92jZzdYtvOOaeXfYr4MEPU562iecrbWlP0KmAQcDrxpyhwDfO177tFR13rGsp1ZwO6W7bSqr5dLpIoSvjRV55vW04OW7QysoUu9Md6MSbJVCWhKVbKP2d6lhnM8VZXsWf/YyzdMK20PkwD/YXaPtmynXczxrwG7mN6G6JbcJ/Ule2Og6Ra/JyrZV3nUJPwhwPA4zrWOZTtbmJb9k+a2RXTcH5tu7INM4v7Csp3vgEMt28mLupWxO9DDdF+vMfVzr2nNVt3maGUeDVr1u/ZKJM66mMcRV/0+LUyXfsQ8VndXc63f4j2fZTutgP2Br2KnMvqe+7F5bvm+lu3k+567Mmr3fVXJ3vjSvHaP2lYO9LRsp0/0zBTfc+vs3RGpibr0pUnyPfc34HpgL+CEJJ8+9rnSVQPZfomJoWp7Xg3nqKlLvCp5Vd0X7mteZ5rnjEf/3Gr2bVxPbLXZqrY4zL37Reb+daKqYj6phpgXAR1iYn7BfPGIHmtQrTsfMwDRsp3/mN6HMtPdvsh8OSGZjRPLdnpYtvO06blZGXWt4xt4rc3N39LaboPMNOeM/dIyK/pN1BfE6M/TLabnYIZlO69ZtvNvy3Y2TTA+EVALX5q4u02yv82ynTd8z12UpPPWNhAqkQFSNY12r/qDXmBeW5nu3f2AtTWUxwzkirYkzutX3fOurbt3uUnEiWplXp8xrfyaRPe2jDRfzP4OjI3qzv/G99wvoso9bXo8xpnxBvNM63YfM1AxKSzbaWnugfcGHjPjBhaZ+r/Q3M5JVDx1HV2uSr2fJ99zR5iu/nNMz4ljfo/JwNm+537egHglpJTwpcnyPbfcsp2qQUy3mEFpiWgeUGhETXOradsK87rUdCV/7XvuH0m+/jLzWttgt1YNnIJXdcyKeFa58z13lmU704DDzH+rXU13ftXgRCzb6W6S/UxzT7s8al+fBsQYLfa/8WEm2T8e3bVvrvWvBl4jnrr+//buJrSuIgwD8JsuKhUKXag1i24btRBQpGKxbRYZsCM1DJUQFEIgpAV/ilAKhVIk0FoIqNDQPygBW4gYKBMTHFvHH2pqW9SNizYYXEjBv5hNN9akSLp5j0xOJrb39go25302B85J7j05J9zvzDffNxf1tjzG4COAyLqM5wB0AegBcM5Y19LAB11Z5pTSl/taDH6Co80+Y90zmR8p5o1zwb1h88IZj2f2FSn0Yi62SAEvOm9j3RrOZdereO0N5QPGumaO7icX/9pdv27uWsNY93Bm9zCARwA8C+AlZgCGk+PFFMDlNNjT1rs4p1rucTGdMlE67xWcHqrHFLMRi641beBovpytqUkM/mYMPsbgewEMsOiv3nOWClLAl+VgL1Pdx9nalPqV25Z0J/vTm//Dc+pm+rh4vwcBdHDO+GvuHuH2Taa6i59t4vz2dWPd6jrf/wKAPwB0ZQoCd3F7ttYXjcH/yOKyJ9mP/g8+cP3G9QdSHzJl/jxT0hMx+LQWobhHC4KzsW4L09hIpkFyZhhwy/f4oUyKPvteAPZxzYDye/3N7ZKLHMXg/2T7Y6uxblPpHNoArAcwHoOf+5e/YRFj3dPGuiljXV/mcDEdcCtzTCRLKX2578Xgp411+wEc5a60iv0SA18f2+h+APAUHxK+ArDFWNfUwCr/InBPA7horDvDD+ceBpm3i2K/GPyYsc4zCH5qrBtmwVYXi9z6k8LAmsTg54x1r7No7qKx7hiAGxwx9wK4wj7xerwK4HMAo8a691h89gSA17gATzp6Rwz+F2PdBfbrrwVwqHT8J2PdFQBtxrojAL7lQkEv81zHAOww1l3NtD8iBj9rrPsEwHZjXT9b7taxFuALAC8m9yUwtb7HWPcX79MLfFh4i50CbxjrVsTgP05G5QPGuiKblLOXo+1xY927zOI8xu8nmGFrYa2+5//OUS509B0fnFp5rSeTtkaRO9IIX5aLE/xAXIBtUAbAN+xLjwywHcmHeSPn8h/g9gCAjwDsZnHhWi6WUy5A6+TiL48COMZU7SoA3ff6TWwx+BGOqn9n+90QgDbWO7TXOuJMXvcy5+I/Yy//EFeLG2HP/8+ZX/uA12A2WeQn1cnr9QqAQY6K27mu/ylmYw4u0REBLkb0Ph8qvgSwn4sIeR5fyXOf5loG13g/DrNKfytXQrzEQsFiVH2S7YbtAHYm8/HlazLFaY7zDPKn+RqjADYu9WVPd7jOt3i/Bvk/PMiH2m0A3gGwWSvtSS2a5ucb2b4sIiIi/0ca4YuIiFSAAr6IiEgFKOCLiIhUgAK+iIhIBSjgi4iIVIACvoiISAUo4IuIiFSAAr6IiEgFKOCLiIhUgAK+iIhIBdwG681QH4edlWwAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "bento_obj_id": "140494428153504",
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "d81134ff-cec6-45cb-92bf-4170f428af40",
        "showInput": false
      },
      "source": [
        "## Predict on some test points\n",
        "We fit a model using the 50 datapoints collected by SAASBO and predict on 50 test \n",
        "points in order to see how well the SAAS model predicts out-of-sample.\n",
        "The plot shows the mean and a 95% confidence interval for each test point."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234892466,
        "executionStopTime": 1650234951004,
        "hidden_ranges": [],
        "originalKey": "5b5b8b72-77fd-4a4d-b46d-49cc842f2c62",
        "requestMsgId": "970977ea-ee5e-46eb-b500-683673ce723e",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "train_X = SobolEngine(dimension=DIM, scramble=True, seed=0).draw(50).to(**tkwargs)\n",
        "test_X = SobolEngine(dimension=DIM, scramble=True, seed=1).draw(50).to(**tkwargs)\n",
        "train_Y = branin_emb(train_X).unsqueeze(-1)\n",
        "test_Y = branin_emb(test_X).unsqueeze(-1)\n",
        "\n",
        "gp = SaasFullyBayesianSingleTaskGP(\n",
        "    train_X=train_X, train_Y=train_Y, train_Yvar=torch.full_like(train_Y, 1e-6), outcome_transform=Standardize(m=1)\n",
        ")\n",
        "fit_fully_bayesian_model_nuts(\n",
        "    gp, warmup_steps=WARMUP_STEPS, num_samples=NUM_SAMPLES, thinning=THINNING, disable_progbar=True\n",
        ")"
      ],
      "execution_count": 13,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234951111,
        "executionStopTime": 1650234951168,
        "hidden_ranges": [],
        "originalKey": "b18a9ebf-7fe7-4c0e-9d4e-882c149128bd",
        "requestMsgId": "25139c91-a34c-4fa8-808f-70c1cf6952fd",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "with torch.no_grad():\n",
        "    posterior = gp.posterior(test_X)\n",
        "median = posterior.mixture_median\n",
        "q1 = posterior.mixture_quantile(q=0.025)\n",
        "q2 = posterior.mixture_quantile(q=0.975)"
      ],
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234951208,
        "executionStopTime": 1650234951511,
        "hidden_ranges": [],
        "originalKey": "b85763a7-759e-4b04-b73d-52db5ce9537e",
        "requestMsgId": "39163b27-e252-4244-9712-f52503e00f74",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
        "ax.plot([0, 80], [0, 80], \"b--\", lw=2)\n",
        "\n",
        "yerr1, yerr2 = median - q1, q2 - median\n",
        "yerr = torch.cat((yerr1.unsqueeze(0), yerr2.unsqueeze(0)), dim=0).squeeze(-1)\n",
        "markers, caps, bars = ax.errorbar(\n",
        "    test_Y.squeeze(-1).cpu(),\n",
        "    median.squeeze(-1).cpu(),\n",
        "    yerr=yerr.cpu(),\n",
        "    fmt=\".\",\n",
        "    capsize=4,\n",
        "    elinewidth=2.0,\n",
        "    ms=14,\n",
        "    c=\"k\",\n",
        "    ecolor=\"gray\",\n",
        ")\n",
        "ax.set_xlim([0, 80])\n",
        "ax.set_ylim([0, 80])\n",
        "[bar.set_alpha(0.8) for bar in bars]\n",
        "[cap.set_alpha(0.8) for cap in caps]\n",
        "ax.set_xlabel(\"True value\", fontsize=20)\n",
        "ax.set_ylabel(\"Predicted value\", fontsize=20)\n",
        "ax.set_aspect(\"equal\")\n",
        "ax.grid(True)"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<Figure size 576x432 with 1 Axes>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAGFCAYAAADuNsSCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3xU9Z3/8VfI5B4DakBIu4qgrbWtFSNUxWvtSfVY256oZRvv1iKobGFRu7XrtrVbbb1Vtq5u2VqVWv3VC8cqnGqOAsqurcRUdFeqtlBZLaAJCpYwSWYgvz/ynTgMk8lMMpnr+/l48IhzLjOffCX58D3f7/fzLenr60NERGQoY7IdgIiI5AclDBERSYoShoiIJEUJQ0REkqKEISIiSVHCEBGRpChhiIhIUgLZDsCyncOAm4DjgDLgNeB633OXm/NVwM3AOUAd8Cpwje+5K7Idu4hIMclqD8OynRLgKeBvwFSgHvgV8JhJJAB3Ap8HTjHnHwSWW7ZzaDZjFxEpNiXZXOlt2c4BwBbgi1E9ihpgB/D3gA+8A5zte+5vou77A7DS99yFWQteRKTIZLWH4XvuO8CzwCWW7dRbtlMBXAZsBVYBR5nHZm0xt64BjslS2CIiRSkXBr2/CkwGOoAg8E+mR/EOMMFcszXmnk5gYhZiFREpWlkd9LZsp9yMYbwGnG4eRV0IPGHZzrEJbi0BEj5L2/LX/+sLlGZ9TD/vhcMhAoGybIeR99SOI6c2HJ5QCC67fF+eWFZFXd1utm8fUzLc98r2b9TPA0cCp/ue+645dpdlO3OBS4Bl5th44O2o+yaYsY9BBUoD1E9sGL3Ii0Tnlk1qxzRQO46c2jB1oRC0tMATy2DsWGhtHdlDpWwnjIjYjBcwPYh2oBc4Fng46vxxUclERKSotbe3097eTm9vL6tXr2br1q1MnjyZffedxn77lXH88Y3cemsjM2aM7HOynTCeNz2FH1m2Mx/oAs4DPg5c5nvudst27gZutGxnHbARWAgcZKbbiogUvcbGRnbv3k1TUxPBYJCenh5qa2vp63uUsrJWnnpq5MmCbA96+567DfgCsD/whhn4nguc43vuanPZAtObWGnOnwY0+Z67MZuxi4jkimAwSFNTE9u2baOnpweAHTt20NW1DWjiiCO60/I52e5h4HvuK8AXE5zvAeabPyIiEsN1XcLhcNxz4XCYpUuX0tLSMuLPyYVptSIiMgLr16+nq6sr7rmuri42bNiQls9RwhARyXMHHTSV0tKauOdqamqYMmVKWj5HCUNEJI+FQvDYY82Ew/FHGAKBAM3NzWn5LCUMEZE8FVln4bqV1NS0ss8+46ioqACgtraWcePG0draSmVlZVo+L+uD3iIiMjxr18Ljj0cW5U3n05/exLx58+jo6GDWrFk0NzenLVmQ7Wq1o6lzy6Y+rQodOa2uTQ+148ipDeP77W+hq6ud995rH/SaxsZGGhsbIy/ztjSIiIikIBSCV1+FI4/sf3366QCN5s/o0hiGiEieiIxZHHccrMjCnqNKGCIieSCSLB55BMrLobY28zEoYYiI5LjoZNE/wE1aakOlSglDRCSH5UqyQAlDRCS3XXJJbiQLlDBERHLbuefCAQdkP1mgabUiIrnttNNg/XqoiV8qKqPUwxARySGhEFx4ITz99IfHciFZoIQhIpI7QqH+R1BLlsD550MwmO2I9qSEISKSAyLJ4uGHoa4OfvMbqKrKdlR7UsIQEcmy2GTh+9kf4I5HCUNEJIvyJVmghCEikl2vvNJfojzXkwWaVisikl2Njf0JY9y43E4WKGGIiGReKAT/8z9w1FH9r5uash1RcvRISkQkgyK1oWbO7H8ElU+UMEREMiS6kGBFRX99qHyihCEikgG5VHV2uJQwRERGWSEkC5QwRERG39e/nv/JAiUMEZHRd/75MHFificLNK1WRGT0WRZs2JB7taFSpR6GiEiahUL9vYqnnvrwWL4nC7Ldw7Bs50SgNc6pMmCJ77kXW7ZTBdwMnAPUAa8C1/ieuyILIYuIJBQ9wO37/T2L6upsR5UeWU0Yvuc+B1RGH7NspwF4BbjXHLoTOBY4BdgIzAGWW7ZzhO+5f8pO5CIie4udDfX444WTLMh2whjEYuBh33OftWxnP+A84Gzfc9eZ87datnOuSRwLsxyriAgU0NTZRHIqYVi24wAzTJIAOMrE2BZz6RrgmCyEKCKyl2JIFuRSwrBsJwD8CLje99xt5vAE83VrzOWdwMRE7xcOh+jcsml0gi0iasf0UDuOXC634Sv/E2DZsvHU1fXx0INbmXJgiM4t2Y4qvvqJDcO+N2cSBnAWUA/8IolrS4C+RBcEAmUjahjp17llk9oxDdSOI5fLbfi5ibBsGeyzTwkzZozPdjijJpem1V4APOR77s6oY5EcHft/YELUORGRjAuF4MUXP3x96qmF+RgqWk4kDMt2aoDPAU/GnGoHes0sqWjHAf+dwRBFRAZElyh/Mva3VgHLlUdSnzLTa1+OPuh77nbLdu4GbrRsZ52ZVrsQOMhMtxURyajYAe799st2RJmTEz0M4CPma0eccwuAZcBKc/40oMn33I0ZjlFEilyxzIYaTElfX8Kx47zVuWVTX64OkOWTXB5ozCdqx5HLdhsWULIoGe6NudLDEBHJaZdeWhDJYkSUMEREknDRRTBpUvEmC3Jo0FtEJKedcgqsX18YVWeHSz0MEZE4QiE47zxYvvzDY8WcLFAPQ0Rkb/FKlNfUZDuq7FMPQ0QkSuxsqCeeULKIUMIQETEKaOrsqFDCEBFRskiKEoaICPDHP4LnKVkkokFvERHgiCP6Z0RVVytZDEYJQ0QKXnt7O+3t7QB0dnYCUF9fT18f7NwJJ57YSGNjIyefnOVAc5wShogUrPb2dlavXs3OnTsJh8NEaueVlJTQ0dHJzp3VPP30CRx6aGO2Q80LShgiUrAaG/t7DgB33HEHK1eu5PXXX6e8vIJAYB5tbX9PXV0lEyYM+VaihCEixaCtrY2rrrqKnp6eqKMXA/NYtGgFM2ZMz2J0+UOzpESkoAWDQSzLikkWETuYP9+iu7s7C5HlHyUMESlot9xyC9u3bx/0fHd3N0uXLs1oTPlKCUNEClYwGOSGG25IeE1PTw8bNmzIWEz5TAlDRAqW67oMtatoRUUFU6ZMyVhM+UwJQ0QK1htvrB9k7OJDlZWVNDc3ZyymfKaEISIFKRSC3/52KlA76DUVFRX4vk9lZWVGY8tXmlYrIgUnUkhwzZpm4Iq411RVVbFp0ybGjRuX8fjylXoYIlJQ9qw6W8mNNy6iurqa8vJyAMrLy6msrGT+/Pk89NBDAyVDZGjqYYhIwYhfovwCvvnNc5g3bx4dHR2MHz+eadOmccUV8XseMjglDBEpGK+/Dr/97d4lytetW8eUKVOYMmUK9fX1ACxevHjgvugSIjI4JQwRKRif+hQ8+SSUl+9ZoryxsXHg0dPs2bOzF2CeU8IQkbwWCsGLL8Kxx/a/Pv74bEdUuJQwRCRvRcYsHnusf9ziy1/e83z0PhgRehQ1fEoYIpKXYge4J03a+xolhPTStFoRyTvxZ0NlO6rCp4QhInlFySJ7cuKRlGU7lwDfAg4CNgE/9T33J+ZcFXAzcA5QB7wKXON77opsxy0imTd7tpJFtmS9h2HZzizgX4ALgLHAPGCOZTuRLbDuBD4PnALUAw8Cyy3bOTTLoYtIFlx6KXz0o0oW2ZALPYzvAlf7nvuCeb3c/MGynf2A84Czfc9dZ87fatnOucAcYGH2whaRbJg5E/78Z6ioyHYkxSerCcOynUnAJ4BSy3ZeBD4ObAB+6HvuQ8BRJsa2mFvXAMdkKWwRyaBQCC68EJpOreSir/cfU7LIjmz3MA40X+cAXzPjF5cCv7Zs5x1ggjm/Nea+TmBiojcOh0N0btk0OlEXEbVjeqgdhycUgssu35cnllXx1FNjOfmkzdTWJt4QSRKrn9gw7HuznTAin/9933P/ZP57kXnkdDHQOsh9JUDCvzWBQNmIGkb6dW7ZpHZMA7Vj6iKzoZ5YBnV18OD97zH5kDiLLSRjsp0wOs3X92OObwAmAVvM6/HA21HnJ0SdE5ECEz11tq4OfB+mHBjKdlhFL9uzpP5sksZnY44fAvwFaAd6gWNjzh8H/HcG4xSRDImXLDQbKjdktYfhe+4uy3ZuA75r2c5LwCvAN4BpwDd8z91u2c7dwI2W7awDNpqZUQeZ6bYiUmD+9Cd46ikli1yU7UdSAD8yPZ1HgP2B14Azfc99yZxfAPwYWAnsA6wFmnzP3ZjluEVkFBx+eH/CKC1Vssg1JX19hTnjoHPLpj4NMo6cBmvTQ+2YWCgEa9b0r7EYjNowbUqGe2O2xzBEpMhFxixOOgmWLs12NJKIEoaIZE30AHdtbX/JD8ldShgikhWqOpt/lDBEJOOULPKTEoaIZNycOUoW+UgJQ0Qy7hvfgAMPVLLIN7mwDkNEiswxx/Qv0Csvz3YkkgolDBFJm/b2dtrb2wHo7OwvFVdfX09fH/zlL3DwwY1cdlkjoGSRj1JKGJbtjAe+DZwO/B1wuu+5q825hcADvuduHrVoRSSnNTY20tjYnxAWL14MwMUXzx4Y4N5/f/j7v+8fu5D8k/QYhmU7DcBLwHxTPbYqsmLQsp16U77jBct2NJNapMgFg0HWrFnDE08s5/jjH+CRR7qpqwPPU7LIZ6kMen/P7Ln9RVNNdmB5ue+5ncBMU+vpO6MTqojkg7a2NhoaGrj//vtZtuwJ1qy5DJjEokVtGuDOc6kkjC8Ad/me68XbvMjsyX0ncFp6QxSRfBEMBmlqamLbtm309PSYXxU7gG0sWNBEd3d3tkOUEUglYUw05ccTeXWorVNFpHC5rks4HI57LhwOs1TFovJaKglje9Qe24M5BPhghDGJSJ5av349XV1dcc91dXWxYcOGjMck6ZNKwlgNzLFspzbeSct2TjQD4s+mLzwRyRehEPT2TqWmpibu+ZqaGqZMmZLxuCR9UplWez3wvHns5JuHkxdbtjPLbLE6DQia60SkAAy2riIiMo02Uhvq0Uebqa6+Iu57BQIBmpubMxS5jIaUNlCybOdY4C7giDin1wKX+577+7RGOEzaQCk9tGlNehRCO0bWVcyePXuP47GFBG+/vY0FC5oIBoP09PRQW1tLIBCgtbWV6dOnD/vzC6ENc8SwN1BKaeGe77m/A460bOdQ4JNADfA34FXfc9cPNwgRyU/xq85OZ9asTcybN4+Ojg5mzZpFc3MzlZWV2Q5XRmhYpUF8z/0T8Kf0hyMi+SJesigtbWfx4v5HWDPMoosdO3awZMkSiFkJLvkn6YRh2c4FSVxWApT7nvufIwtLRHLd3LnxSpQrIRSyVHoY90Yt2It+BhY9CFJiXithiBSQSKmPjo4OamtraW5uZs6cSp55Bn79a5UoLxapJIyLBzleDhwMtAB/ABalKTYRyQFtbW00NX04iL1ixQquuOIKWltbef316ao6W0SSThi+596X6LxlO983020P01oMkcIQXeojYseOHQA0NTWxefNmQIPZxSJtO+75ntsD3AJ8M13vKSLZpVIfEi3dGyh1AJPT/J4ikgHxFum98MILAz2KWCr1UXzSljAs2ykBLgDeSdd7isjoipck6urqePnll9m6dSv77rsfpaU17Nq1d30olfooPqlMq12R4HS5KTw4Hrg9PaGJyGiL3SHvzTff5Oabbx4Y4C4trWXXrp1x71Wpj+KTSg/j5ATnwsBm4A7gn9MQl4hkWG9vL4sWLWLnzg8TxK5dkcdRJZSXl9Hb27tHqQ+t3i4uqcySStsAuYjklmAwyKOPPmo2PdpbbW0NRx99NHV1dSr1UcTSPegtIjkudtzi7bff5r777qOnp4ddu3bFvaerq4u6ujrOOOMMWlpaMhyx5IpBE0aSpUDi8j13SbLXWrbzJvARIPZv6hG+575h2U4VcDNwDlBnyqtf43tuojEVEYkRb4A7FApxzz33DLl1ak1NDePHj89InJK7EvUw7o23d/cQIqVBkk4Yxjd8z713kHN3AscCpwAbgTnAcst2jjBFEEUkCbED3ABr1qwhmS0OAoEA06ZNG/UYJbclShiDlQLJGMt29gPOA872PXedOXyrZTvnmsSxMMshiuS1jo6OQcctAMrKyigrK2Pu3LmUlZVBVLJB1WeLzqAJY6hSIPFYtrMvsO8w4jjHsp1vAQ2mbPr3fM9dBhxlYmyLuX4NcMwwPkdEouy//3jGjKlm9+69p86WlpYyc+ZMzj77bMrKypQcJO2D3l8GbjC/+JP1CrAemG02Y5oPPG7Zzkxggrlma8w9ncDERG8aDofo3LIp9e9A9qB2TI9ca8dwKATAa3+8nN27HwX2Thi1tbX88p6795gNlc3vIdfaMF+NZNfClBKGZTsfAS405T9i760CTjWL+JLme+6XYg5db9nOl00C8Qe5rWSo8ZVAoEzbOaaBtsVMj2y0Y6L9uLdt304gEGDCAUuoqHAoKfk1u3eH91pn8dHJubOSW38Xsy+Vld4fA34PjDOH+uLsi7EL+EEa4vozMAnYYl6PB96OOj8h6pyIxBFvkHv27Nn09cF3vnMtixYtYteuXQP7bvf2wgknnMCcOXO0zkLiSqWHEUkEc4E3gd8ClwOvm1XgFwKX+Z77VLJvaNnOwcA1wLd9z90WdeqTwEqgHeg1s6Qejjp/HLAshdhFilb05keVlbW47ul43k/p7f3wMVSkwGB7e7uShQwqlYTxWeBW33N/ZtnOWHNsne+5zwErLdt5EFhl2c45vueuTvI9twBnAnWW7cwzyeEq4FDgLN9zt1u2czdwo2U768y02oXAQWa6rYgkELv5keetIBzeZar57G3Xrl0sXbpUi/MkrlTKfRxgBqcBdpuvA+MVvue+BtwNXJ/sG/qeGwQ+D9Sa2VH/Z3orJ/ue+7q5bIHpTaw05dNPA5p8z92YQuwiRSd686PI1NlweAcQBEJx7+np6eHXv/41ixcvZvHixQNjICKk2MPoNCuyAbpM0vhozDV/BK5MJQCTaL6c4HyPmTk1P5X3FSl2iTY/GkxtbS2zZs1SD0PiSqWH8Xvgm5btnOJ77m4zdnGFZTs1UdecFnd+nohk3Pr16+nq2nsfi0RUslwSSaWH8WNglRmkXmkeP90CbLRs52XgQGAKkPKCPxFJv6lTp1JVVcPOnfF3zItVWVmpkuWSUNI9DN9zXzSzlSI1n35iigJWmDpPBwEPmTEHEcmy5uZmysuT+zdhZWUld911F9OnTx/1uCR/lSRTeCwRy3YCQD3Q4Xtu/NrIWdC5ZVOfFvmMnBZLpUcm2zEUgtWr4XOf23uW1GDGjRvH5s2bc7p3ob+LaVOSxDVxpbJwbyXwS+AR33M/iBz3PTesRXQiuSEUgpYWePRRuPdeuOCC6WzatIl58+bR0dHB0UcfzW233TaQQCoqKigtLdWjKElKKmMYJwEnAndYtrMMuB/wTMIQkSyLJItHHoG6OjjssP7jVVVVzJgxA8xK76uuumoggYwfP55p06bpUZQkJZWE8XdmE6OvAmeZP+9btvMQcL/vuc+PYpwikkBssvB9MDliL9EJRCQVwxrDMEUII8njs+bwRtPr+FXUorus0RhGeui5cXqMZjsOliyiiw/uFU9nJ9XV1VRXV4PpeeQ6/V1Mm9Efw4jme+5fgduB203yOAtwgH8CrtVe4SKZc+WV8XsWifaviN4ESSRZ6fjFvsOsAn8b+GCYGyiJyDBdfjk88ww88MDgj6Gitbe371XuXLvoSTKGlTAs26kHvmJ6Fp8z7/MBsBR4MP1hiki0vj4oMQ8WPvMZeO01CCT509zY2DjwqCofHkVJ7khlWu1EoNkkiRPMvUHgN8ADZsZU7+iGKyKhEJx3HnzhC3DJJf3Hkk0WIiORyl+zv5qvu4CnTU/C9T03uboDIjJi0QPcra3wla/AfvtlOyopFqkkjN8BvwIe9j23cxRjEpE4YmdDPfWUkoVk1ohLg+QqTatND01lTI+h2jHR/tsARx7ZyM03Nya1ziKZz4gn1we79XcxbYY9rVYJQxLSD2l6pNKO0ftvk+KivEKmv4tpM+yEkcp+GCKSBW+9BatWFXeykNyguRUiOW7KlP51Ft3dShaSXUoYIjkkGAyyZs0a3nmng7/8pZbvfreZyspKjjgi25GJ6JGUSM5oa2ujoaGB+++/n2XLnuBHP7qM/fefRFtbW7ZDE4FEPQzLdg4cxvuVABW+574xsrBEClvsrKhQKMRNN93Ezp07o67awc6d0NTUlPObG0lxSPRI6k1guFOoSod5n0hRiJ7CunjxYtasWcOYMfE7/OFwmKVLl9LS0pLhKEX2lChhLImTMD5m9vV+FXjNlAbZBzgcOBTwgbWjHLNIwYiMWaxe/V/s2BG/aEJXVxcbNmzIeGwisQZNGL7nXhT92rKdmaa44PHxNkuybOdksx/G9aMVrEghSXa/7ZqaGqZMmZLR2ETiSWXQ+0bgZ4PtrOd77irgP4AfpS88kcIUDAZpampi27ZtCZMFQCAQoLm5OWOxiQwmlYRxFDBUv3gjcOQIYxIpeK7rEg6HE15TW1vLuHHjaG1t1YC35IRU1mH0AicD9ya45mQg8U+BSBGJzIYKh0Js274dTI2oZcuWDzpmAXDYYYdx3XXX0dzcrGQhOSOVhPEkcL5lO7tMafP1ZtC7Ephi9vi+yIxziEjUbKi3/rKe7//wRjo6OjjrrFls2eIAK8yGlXuqqKhg5syZmhUlOSeVhDEfOAy42CSGWCXAG8A/pjE+kbzX1taGZVl0d3fT09OD560gHB78R6+0tJRp06ZlNEaRZCSdMHzPfdeynaOBLwGnAAcD1aaX8X/As2ZDpdBwgzEzsZ4DfuB77vfMsSrgZtODqTNTeq/xPXfFcD9HJFMig9vbzeMogHC4v1dRXb0PZWVjBxJJbW0tgUCAuXPnUlZWlsWoReJLqZaU77m7gcfMn7QyieGeOH30O83aj1PMoPocYLllO0f4nvundMchkk6JBrfHjOnj9ttv57/+67/YtGkTDQ0NTJs2bSBZRMqckwd7VUhxSLn4oGU79cCpwIHA//M99y1zvMb33K4RxHKDWQy4Keqz9gPOA872PXedOXyrZTvnmsSxcASfJzLq1q9fT1dX/B+Lrq4u3n77bWaYErSR/S9EclVKCcOynX8GvgNUmFXgLwBvWbZTCrxu2c7PfM/9QapBWLZzPHA+8GkzoB5xlIkxtvraGuCYVD9HJNOmTp1KTU1N3BlR5eXlvPnmm0yYMAHUo5A8kHTCsGynxazifsHMhPpx1OlaUxLke5btrPc994EU3rfaPIqa73vuZst2ok9PMF+3xtzWCUxM9jNEsiEUgtraZgKBK+Ker6qq4o477tC0WckbqfQwrgCeB04wg88DCcP33O3AFy3bWWGuSzphmEdR63zPvT+Fe0qGKowYDofo3LIp0SWSBLXj8IRCcNnl+/LEsirm/8Oj3HOvQ3dPT//gdk0NpYEADz34K3Zsey/OxFqJR38X02Mk29ymkjA+BXzb99y+mF5AtAdjeh4JmUdRLcBg28NsMV/HA29HHZ8QdS6uQKBM+/+mgfZRTl1kD+4nlsHYsfC1cz/HP87/w8A6jFmzZmlB3jDo72L2pZIwyoChBrVDZnwjWV83vZX/jUpCY4EZlu1Epu/2mllSD0fddxywLIXPEcmISLJ45JH+ZPHzn7ezdm07L4ZCA4PbO3bsYMmSJaCxCskzqSSMP5vHUfcluKbFXJesfwSuizn2MPA74Cbfc7dbtnM3cKNlO+vMtNqFwEFmuq1IzohNFq2tMGNGI9Cofx1LQUglYdwP/NCynf+N+td9iWU744HPmpXgpwDXJPuGvue+D7wffcyynR7gA99zI4+cFpjHXCvN3htrgSbfczemELvIqPvmN2OTRbYjEkmvkr6+5DbVM1NnHwa+YgacS4Dd5ivmqwucYxb4ZVXnlk19+hfdyOlfxh+K3VYVU0gwoqGhkauuamTJkr2Thdpx5NSGaVOSxDXxb0w2YURYtvMV4Ktml70a4G+mXMdDvuc+MdxA0k0JIz30Q7q3YDDIvHnz6Ojo4KtfncVZZ304gL1rF5TG2aBY7ThyasO0yVzCSMSynQBQ6ntu4h1hMkAJIz2K9Yd0sN7Em2++yaJFi9i1axc9PT0EArWUlQV49tlWpk+fPuj7FWs7ppPaMG2GnTBSWbi3AZjre+5TCS5bCFwAfHK4AYnkgujZS5EV2Oeffz4NDQ3s3Llz4LpweAfhMFhWE1u2bNZUWSloqQx6TzbVaeOybKfclD8/KD2hiWRPe3s7q1evZvv27bzyyiu8//77/PznPx+0LtSuXWGWLl2qPSykoCVMGJbtfBf4F/OyD3gkwaK9iFfSFp1IljQ2NrJ7926ampro6uoiFApRUlLCYI9wu7q62LBhqB2MRfLbUD2Mn5u9Lqab6rDrgI5Brt0F/AW4aRTiFMmoyD4W27ZtGziWaLyvpqaGKVOmZCg6kexImDB8z/2rKQx4j2U7c4B/8T1XW7BKwXNdl+7u7qSvDwQCNDc3j2pMItmWyo57Y+h/TFUG9PmeO7ArjGU7+wJdvuf2jlagIpn02muvJZUwIrvktba2asBbCt6YVC62bOdyU/TvUzGnvgR0WrZzZXrDE8mOrVtjK+rv7bDDDuNnP/sZmzdvTjilVqRQJJ0wLNs5C7gD2AbErrNYB2wAFlm2c3b6wxTJrOgV3PGUlpYyc+ZMWlpa1LOQopHKtNqrgGeBL8Q+evI9t82ynUbgKXPdI+kPVSRzpk79OKWlVezaFYx7vrS0lGnTpmmXPCkqqSSMw82ueHHHKXzP3WXZzv3AovSFJ5J5oRA89lgzu3ZdAeydMEpKSnjmmWc4/vjjsxKfSLakMoaxGxg3xDU15jqRvBQpUe66ldTUtLLPPuMoLy8H06uorq7mueeeU7KQopRKD+P3wGWW7dxrypLvwbKdj5hS5L9Lb4gimTN/fn+J8ro6uPvuMWzZ8gNWrVrF+++/z8EHH8y0adNYt24d69at0yMoKTqpJIzvA6uAP1u281tgvemv15rHVdCGlm4AACAASURBVE1AKXDeKMYrMqr+4R9g5Uq4/vp23nuvnfLyco466iiIGQhXspBilFK1Wst2TjIzpeIVF3wVmOd77qq0RjhMqlabHsVQIbSvD0qi6nfGliiPDGzPnj172J9RDO042tSGaZPZ8uaW7XwM+ITpXewA/uh77hvDDWI0KGGkR6H/kEbGLI4/vn/HvIjo8ubxpNrDKPR2zAS1Ydrkxn4YuUQJIz0K+Yc0dg/uN96ACRNG57MKuR0zRW2YNunfD8OynQuAlb7nvhX1Oim+5y4ZbkAimRCbLFpbRy9ZiBSKRIPe9wJnA29FvY50R+JlqMg+332AEobkrHjJInYPbhHZW6KEcTHQFvNaJK8pWYgMn8YwJKFCe2781ltwzDHQ1ZXZZFFo7ZgNasO0Gf09vUUKwd/9HaxaBdu2gQrMiqQm0aD3cPabLAHKfc/9yMjCEkmfUAieegq++MX+14cemu2IRPJTolpSY0wCiP5TAUw2f/YxCWffqGNh4O0Mxi+SUGTM4swz4Y47sh2NSH4btIfhe+7k6NeW7TQAPvDvwM98z90adW4ScBnwFcAZ7aBFkhE9wF1Xp8FtkZFKZQzjduC/fc+9IfaE77mbge9ZtlNvrvtqesMUSU1ssvB9JQyRkUqlvPmpwPNDXPO8uU4ka5QsREZHKgmjEpg6xDUfN9eJZM2CBUoWIqMhlUdSLwILLdt5A3jY99zuyAnLdsrN2MVC4KXRCVUkOfPn90+d/cUvlCxE0inVPb2fMSVC7rZsZ7PZD6MSmGTeqwf4dioBWLbzaeAG4FigHHgd+KHvuY+Z81XAzcA5QJ0po36N77krUv5upWDt3g1jTH/5kEPglVc+fC0i6ZH0j5TvuW3Ap4GfAn8002kPAfYH/gz8DJjme+7qZN/Tsp1aYKW5/2BgPPAY8IhlO4eby+4EPg+cAtQDDwLLLdvRbHoBM2YxaxbccsuHx5QsRNIvq6VBLNsZD3wZeMD33J3mWB2wHZgFPA28A5zte+5vou77g6mku3Cw91ZpkPTI9XIMsbWhXn8dDjgg21HtLdfbMR+oDdMms6VBLNs5ADgQeM333L8N98N9z+0Afh71vvXmkdbbwArgKBNjW8yta4Bjhvu5UhjiFRLMxWQhUihSShiW7ZwC3Ap8xhw6BXjOnHsc+Dffc58eTiCW7fSYMYx2wPI9t9OyncgOBVtjLu8EJiZ6v3A4ROeWTcMJRaLkajuGQnDZ5fvyxLIq6up289CDW5lyYIjOLdmOLL5cbcd8ojZMj5H00pJOGJbtfBZ4EtgJLAO+GHVuf9MbeMKynZN9z30h1UB8z60wPYx/AJ63bCdRD6Ikam+OuAKBMnVf0yAXHwNEehZPLIv0LMYwY8b4bIeVUC62Y75RG2ZfKkOD1wKbgMOBC6Ofg5kyIZ82my2lNEsqmu+5nb7n/osZt5gDRP69GPvbYELUOSky774LL7yg/SxEMi2VR1LHAbf4nrvZsp2xsSd9z33fsp2fmem3SbFs5wzgP4DDfM/tijpVYgoZtgO9ZsrtwzGxLEshdslh7e3ttLe3A9DZ2QlAfX39wPnGxkYaGxsHXn/kI7ByJbz3nkqUi2RSKgmjLolKtJvNdNtk/d6MW/zUsp2FZl3HZWa67lLfc7dbtnM3cKNlO+uAjWZx4EFmuq0UgOiEsHjxYgBmz569xzWhEHgefPnL/a+nTu3/IyKZk8ojqc3AYUNcc6J5bJUU8yjLMo+Y3jCPoloAx/fc35vLFpjexEqgAzgNaPI9d2MKsUsei4xZfOUrsGhRtqMRKV6p9DCeAK60bOcxs9COyMCzZTv7AVcCX0/1X/6+574SPYAe53wPMN/8kSITO3X22GOzHZFI8UolYXwfsM1jpPUmWfybqSM1xWyutAG4fhTjlQIWDAZZs2YNHR0d1NbWcuaZzVxySeUe6yw0wC2SPSmt9LZsZ1/gu6au06SoU38FHgL+1ffc90cl0hRppXd6ZGoqY1tbG01NTQSDQXp6eqitraWnJ0Ao1MrYsdPzPlloSujIqQ3TZtgrvYddGsSynXFADfA333M/GG4Ao0UJIz0y8UMaDAZpaGhg27Ztcc6O47nnNnPCCfldNV+/7EZObZg2o1saxLKdCuAu4F7fc5+jf2xhGxDvJ1wkJa7rEg6H456rrg7z1ltLzVwIEcmmpGZJmYHnswBViJW0W79+PV1dXXHPBYNdbNiwIeMxicjeUplWex8w25QkF0mbqVOnUlNTE/dcTU0NU6ZMyXhMIrK3VGZJvWAWzG20bMczi+jijl34nntT+kKUQnfmmc1cdNEVcc8FAgGam5szHpOI7C2VhPFLM5W2BDg3wXV9gBKGJCUUgq9/vZJQqBVooqysm1Com9raWgKBAK2trVRW5veAt0ihSCVhXDyKcUiBSlQnqq8P/vCHRh5+uJG6uuk88cQmliyZR0dHB7NmzaK5uVnJQiSHJJ0wfM+9b3RDkUI0WJ2oUAjOPRcefhjq6sD3YcaMKl57rX+xRUuLZkWJ5Jpkp9UGgE8CZcCrvucGRz80KSSxq7iPOaaZNWsqo5JFtiMUkaEMuXDPsp1zgUVRVWi7gRt8z/1hRiIcJi3cS4/hLJaKfQz19ttvc9999xEOh+nt7aW6upry8nLuvbeVhobpjBnz4fXxxJY3z0dadDZyasO0GZ2V3pbtHA88C/wNWGWSxUygAfhH33NztnaoEkZ6jPSH9I477uBb3/oWO3fu3OvcuHHj2Lx5c1GMU+iX3cipDdNm2AljqHUYC01J8U/5nvsV33P/HviYSR7XWrYz7A+WwhcMBnn00Ufp6emJez4cDrN06dKMxyUiwzNUwjgW+E/fcwc2TjLjF9cD9cDHRz9EyUdtbW00NDSwevVqdu3aFfeari6t4hbJJ0MNetcDr8c5/prp1tTHOSdFLhgM0tTUNEgxwQ9pFbdIfhmqhzEG2Pvhc/9YRjL3SxFKVEwwmlZxi+SXVBbuicQVmRXV29vL6tWrefHFF9mxY8eg15eXl1NdXa1V3CJ5RglDRqyxsZHdu3fvsQHSYEpLS7n00ku59dZblSxE8kwyCeMrlu0cEnOs0tSMarFs55jYG1R8sLgkO2YBUFFRoWQhkqeSSRjnJTg3O84xFR8sMsmMWVRUVFBaWso3v/lNlixZMnC8EBbliRSLoRKGCg7KkBJtgARw2GGHcd1116mYoEieS5gwVHBQknHQQVMpLa0hHN57oLuiooKZM2eqmKBIAdCgt4xIKASPPdZMOBx/A6TS0lKmTZs2UKkWPYYSyVtKGDJsoRC0tIDrVlJT08qYMU309vbPkoreAGn69OnZDlVE0kAJQ4Zt61Zob4eZM9s566yXKCv7AWvXrqWjo4Px48czbdo0XnrpJcaMGaMehUgBUMKQYZs4EVatgnffbeToo5UQRAqdSntISkKh/l3yIg48EI4+OpsRiUimqIchcUXKfYRDIbZt3w7A/vvX85e/wPvvwxtvNPKd76hXIVJMlDBkQOxOeQDVVVW8/PLLbN26ld27T2HlyoWMHVtJa2uWgxWRjBtyi9bRZtnOAcCPgdOAKuBV4Frfc1eZ81XAzcA5QJ05f43vuSsSva923BuZxYsX8+abb3LnnXfS3d1NT08vUAMEuOeeVi66SDOfUqHd4kZObZg2o7bjXib8BjgAONJ8fRZYZtlO5G/GncDngVPM/hsPAsst2zk0y3EXnPb2dhYvXsxtt93G4sWL+fGPf8z27dtNMcE+YAewjQULmuju7k7iHUWkkGQ1YVi2E+kxLPA9d4vvud2mt1EDHGPZzn6mltW3fM9d53tul++5twJ/BOZkM/ZC1NjYSEVFBVdddRXt7e3s3r077nXaWlWkOGV1DMP33A+Ar8ccjmzBtgk4ysTYFnPNGmCvKrkyMu+99x4XX3wxQz2m1NaqIsUppwa9TY/jHmC577m/t2wnUoBoa8ylncDERO8VDofo3LJp9IItQN/69rVDJguAmupqxu+3r9o3Bfr7OHJqw/QYyThQziQMy3YOApYB7wJfG+LyEvNQfVCBQJkGyFL05v+9ldR1gbIyLrzk66o8mwIN2I6c2jD7ciJhWLYz3SQLF5jne27InNpivo4H3o66ZULUORlCZLrszp07+eCDD/Y4FwgEqK6uprq6mv322z/h+wQCAWpra7W1qkiRynrCsGznU8CTwPW+5y6KOd0O9ALHAlHriznOJBhJQnR12OiqsQCzZ/fvgRUKwfLlXwUeGbTzduedd3L++ecrWYgUqawmDMt2SoH7gDvjJAt8z91u2c7dwI2W7awDNgILgYPMdFtJg0jV2ccfH0dl5T10d1+8V9K48MIL6evrG9gtTyXKRYpPtnsYx5qZUJ+ybOfqmHO/9D33G8ACM9V2JbAPsBZo8j13Y5Ziznu9vb0DVWUrK2t57LFmXLeSsWOhtfVCDjnkTGzbpqOjg5NPPIFbf3I748aNy3bYIpJlWV/pPVq00ju+a6+9lkWLFrFr1y56enqoqaklGAxQVdXKbbe9RWfna3tc39fXR0lJ/8LQww47jObm5ixFnt80YDtyasO0GfZKbyWMIvLee+8xadIkent79zq3zz7jePfdzXuNT+iHND3UjiOnNkybvC4NIqMoUu7j2muvZeLEiXGTBUBfn1Zvi0hi2R7DkFHW2NjI4YcfTkNDA6FQaNDrtHpbRIaiHkYRuOWWW9i2bVvCa2pqapgyZUrCa0SkuClhFLhgMMgNN9ww5HWBQEAD2iKSkBJGgXNdd8j6UNXV1Vq9LSJD0hhGnou3S159ff3A+ZdeesnsZxFfWVkZf/3rX7XOQkSGpISR5+KV/YiU+wB44IEH+OUv76era0fc+6+77jolCxFJih5JFYhgMMiaNWtYvnw5DzzwAF1d3fT1QXNzM2Vl8f9dUFVVxdVXxy6wFxGJTwmjALS1tdHQ0MD999/P448/zmWXXca++07i4ovbqKiopLW1lXHjxlFRUQFARUUF1dXVzJ8/X+MWIpI0rfTOc8FgkIaGhkGmzY7jjTc2c+ihlQSDQWbPns2mTZuor6/nxBNPpKysbI+r4xUU1Ora9FA7jpzaMG2GvdJbYxh5znVdwuFw3HPV1WHa2pZy6KEtVFVVccIJJ9DZ2Ul9ff0e4xwiIslQwsgTg82GWr58OTt2xB/QDga7WL16Ne+++y47d+4cON7Z2TmwNqO6upoTTjhBpcpFZEhKGHlisNlQtbW1rFixIm7SqKmp4YQTTqClpWWvcyIiqdKgd55rbm6mtDR+3tfqbRFJJ/UwctBQi/F27txJdXU1AJWVlfh+K5bVRDDYTW9vN7W1tQQCAa3eFpG0UsLIQYkW4wWDQebNm8c773Tw8su13HJLM9OnT2fz5k3MmzePjo4OZs2aRXNzs5KFiKSVEkYeaWtro6mpiWAwSE9PL7CCX/ziCp57rpXp06czY8YMAI1ZiMio0BhGnggGgzQ1NbFt2zZTG6oP2EF39zaampro7u7OdogiUuDUw8hhkXIfHR0dvPzyy4OutwgGg1x55ZUDPYzIYywGWYwnIjIcShg5as/HTz2Ul5cPur1qb28vkydP1mI8ERlVShg5KPrxU8RgyQLtliciGaIxjByUqNxHPFpvISKZoB5GDlq/fj1dXV2Dni8tLWXXrl1abyEiGaWEkYOmTp1KTU3NoOU+pk+fTl1dndZbiEhG6ZFUDjrzzGZ6euLn8rKyMs4++2zOOOMMWlpalCxEJGPUw8gxoRBcckkloVAr0ERZWTehUPdAKZC5c+cO7GOh6bMikklKGBkyVH2oxsZGjjiikZYWeOQRGDt2Oo8/voklS1TuQ0RygxJGhiSqD4XpWXyYLKC1FWbMqOK111TuQ0RygxLGMCTTW0j18dD27fDqq9HJIs1Bi4iMUE4kDMt2DgbuAU4CDvY9982oc1XAzcA5QB3wKnCN77krshXvUL2FRKLLfdTW1g48Zqqvh5Ur4e23AdpZvLh9j/s0XiEi2Zb1hGHZjgP8B/DkIJfcCRwLnAJsBOYAyy3bOcL33D9lONyURfdG1q5dyz333MPu3bvp7e3lqaee4tJL57Bq1TPMmDGdAw6AAw4AUEIQkdyT9YQB7AecCPwdcEH0Cct29gPOA872PXedOXyrZTvnmsSxMDshJy/SGwgGgyxYsGCPqrL9VWd7OOmkJt5/f7MGtEUkp2V9HYbvuXf7nvv6IKePMkmtLeb4GuCYDISXNq7rsmvXrrjnxowJs3Tp0ozHJCKSilzoYSQywXzdGnO8E5iY6MZwOETnlk2jF1nkc0Kh/oBiPuvVdX/k1T++BsD727bx9NNPmx7F3oLBLv7n5bV0fu7kUY83VZlqx0Kndhw5tWF61E9sGPa9uZ4wBlNidhAaVCBQNqKGSVbALKKL/ayTJjZw0udOBeDaa6/l+eefH/Q9ampq+PRnjsxIvKnq3LIpJ+PKN2rHkVMbZl/WH0kNYYv5Oj7m+ISoczktGAyyaNGihOXJVW1WRPJBrieMdqDXzJKKdhzw31mKKSW//vWvEyaL6upqVZsVkbyQ04+kfM/dbtnO3cCNlu2sM9NqFwIHmem2Oa2trY25c+cOurdFSUkJV199NdOnT894bCIiqcp6wrBs53WTACK9ndct2+kDful77jeABcCPgZXAPsBaoMn33I1ZDj2hyK550dNoY9XU1PCxj30so3GJiAxX1hOG77kfH+J8DzDf/MkbDz/ssmNH4l3zdu/ezfvvv8/ixYu1eltEcl7WE0YhCoXgppvWEw4PvmteZWUlq1at0uMoEckbShgjEK8uVGlpJeeeC6++OhWoAfbeNa+iooK77rpLyUJE8kpJX1/C5Qx5q3PLpr7RnLPd1tZGU1MTwWCQnp6egf21GxtbeeaZ6eyzTzcwib/9bdte944bN47Nm/OjFIjmvqeH2nHk1IZpUzLcG9XDGIbIgPa2bR8mg8j+2ytXNrHPPpt5+ulKSkpaOfXUU+np6aG3t5eKigpKS0uZO3cuS5YsAVWeFZE8ooQxDK7rDjpVtqoqzLe+tZQZM1qA6bzzzjvMm6dd80Qk/ylhDMP69evp6oo/oL1zZxd9fRsGXldVVTFjhnbNE5H8V/QJYzi7502dOpWampqBx1DRampqmDJlyqjHLSKSaUWfMIaze15zczNXXHFF3HOqCyUiharoE8ZwlJZWctRRraxY0QSEgS5qa2sIBAIDdaGiey4R2mZVRPKZEkaKQiFoaYEVK6ZTV7eJk066lpKSDXsNaCshiEihUcJI0XXXwSOPwNix0Npaxdq1nwA+oQFtESl4uV7ePOdcdRWceCK0toKZ/CQiUhTUwzDilfmIPF4Kh6G0FEpKoL4eVq3q/28RkWKiHoYp89HQ0MD999/P448/zmWXXcakSZNoa2sjFIKvfQ3mzYNIFRUlCxEpRgXbw3h13R95/fFlMMT6ikRlPpqamjjllM24biVjx8LChXDwwRn/VkREckLBJoxPHv4JTvrcqTDE+opEZT527AjjuksZO7aF1lYlCxEpbkX/SCpRmY9wuIuKig0a4BYRUcL4sMxHfDV85ztTlCxERAr5kVSsnTt3Djyaih7TCIVC7N69O+49++wT4OqrVeZDRIRiShjV1dUDYxixYxozZszYazOkkpIAzzzTqlLkIiJG0SSMRI48cjo33LCJtrZ5bN2qfStEROIp+jGMSG2oyy+v4u2353HGGWfQ0tKiZCEiEqOoexiRZBGpDdXYuCGJu1AlWhEpSgWfMGJLfpx++umsWbOGd97p4O67a1mzppmxYytpbYW1azuSek8lBBEpRgWdMNra2vYYzPZ9n+7ubgKBAKFQGFgBXMHtt7cyY8Z01q7NdsQiIrmrYMcwokt+9PT0DBzr6+sjFAoBfcAOYBsLFjTR3d2d7ZBFRHJawfYwvN8+OWjJj1jBYJArr7ySGWaFnsYjRET2VrAJ48033xy05Ees3t5eJk+ePORe3iIixaxgH0lNnjw5QcmPPdXU1DBlypRRj0lEJJ/lRQ/Dsp0q4GbgHKAOeBW4xvfcFYPdc4Z9Ov/0nX9O6v0DgQDNzSoBIiKSSL70MO4EPg+cAtQDDwLLLds5dLAbKisr8bxWysrGAbVACRUVVZSUlFBeXg5AbW0t48aNo7VVJUBERIaS8z0My3b2A84DzvY9d505fKtlO+cCc4CF8e4LheC226YTCm2irOwhjjnmKebM+SKnn346V199NR0dKgEiIpKKnE8YwFEmzraY42uAYwa76bLL9+WJZTB2bBVz59Zx8MEn09LSAqbYIDDwWkREhpYPj6QmmK9bY453AhMHu2njxlLq6jA75SW3gltERAaXDz2MwZSY1XdxvfxKeUnkv2fM2HO6rKbPJq9+YkO2QygIaseRUxtmXz70MLaYr+Njjk+IOiciIqMsHxJGO9ALHBtz/Djgv7MUk4hI0Snp6xv0qU7OsGznTqAJ+DKw0cyMuhr4pO+5G7Mdn4hIMciHHgbAAmAZsBLoAE4DmpQsREQyJy96GCIikn350sMQEZEsy+dptXsZTs0pAct2DgB+bB71VZl2u9b33FWoXVNm2c5M4DngB77nfg+1YUos27kE+BZwELAJ+KnvuT9B7ZgUy3YOA24yE4PKgNeA633PXc4I27DQehgp15wSAH4DHAAcab4+CyyzbCcy8V3tmiTzw3iP2Z0rmtowCZbtzAL+BbgAGAvMA+ZYtjPdXKJ2TMCynRLgKeBvwFTTRr8CHjOJhJG0YcH0MIZbc6rYWbYT+RfGzb7nbjHHfgz8E3CMZTur1K4pucH8i25T5ID+bqbku8DVvue+YF4vN3/UjsmZABwIPOB77nb62+1uYBHwGct23h1JGxZMwhhuzali53vuB8DXYw5HNgfZpHZNnmU7xwPnA582/2qLUBsmwbKdScAngFLLdl4EPg5sAH7oe+5Daseh+Z77jmU7zwKXWLbzgulpXGZKK60aaRsW0iOpYdWckj2ZHsc9wHLfc3+vdk2OZTvVpt3m+567Oea02jA5B5qvc4Cvmbb5BfBry3ZOUjsm7avAZLMEIWieFpzte+47I23DQuphDCZhzSn5kGU7B5n1Lu+aH9hE1K57ugFY53vu/SncozbcU+T30fd9z/2T+e9F5nHJxUDrIPepHQ3LdsrNGMZrwOlmLO1C4AnLdmKrZURLqg0LqYehmlMjYAYV15hyK6f5nvs3c0rtOgTzKKrFdP3jURsmp9N8fT/m+AZgktoxKZ83k1cW+J77ru+5O33PvQv4C3DJSNuwkBKGak4Nk2U7nwKeBG7wPXeO77mhqNNq16F93UxP/F/Ldjot2+kEZgLXWLbzB7Vh0v5sksZnY44fYn7hqR2TVxLzOmB6ECNqw4Ja6a2aU6mzbKfU9Cw833OvG+QatWsClu3sC9TEHH4Y+B1wk++5W9SGybFs59tmKu1XgFeAbwC3A0f7nvuS2jExy3bGAX80j+/mA11mVtR/Aif7nrt6JG1YaGMYC8wCtJXAPsBa1Zwa0rFm5sSnLNu5OubcL33P/YbaNTHfc9+PfYxi2U4P8EFkqrLaMGk/Mk8+HgH2N8/iz/Q99yVzXu2YgO+52yzb+YIZU3sDKDdfz/E9d7W5bNhtWFA9DBERGT2FNIYhIiKjSAlDRESSooQhIiJJUcIQEZGkKGGIiEhSlDBERCQpShgiecqynVWW7WhevGRMoS3ckwJl2c73zF4JyXjW99yTRzkkkaKjhCH54iHgf2OOfR843JSP2BZ1vCPDsYkUBSUMyQtmd7B10ccs27nS/OeyqBIcCVm2U+l7bveoBClS4JQwpGBZtnOR2dToAmA6cBGwFLjIbD17ku+5JTH3TDaVUe/zPfeiqOPjzV7TXzKltj8AnjcVfn+fIIZDgD8Brb7nfiHO+e8A/wpc4nvuPebYHLP/wyFAJfAm8ABwm++5wSS+34t9z7035txg3++lpiz74ebQ68C9wB2+5+5OsqmlSGjQW4rBOabI4lXml2FKzF7SvzeJ51emnPmtwGeA5yzbOXWwe33P/bPZDvMUU9U2XmzdwKPms/4VuMs8YrsGuNwUh/tX89lpY9nOT0wV07dMZdOFwGaz//PP0/lZUhjUw5BicAIw2ffc7cO8/1/MlpfH+p67JnLQsp0lwKsmeRyZ4P4HTQ/nS8B9UfcfapLOQ2ZvdYADgN8CX/I9N2yO3WfZzhTAsWzno77nvj3M72OAZTvTTJK4w/fceVGn/sOynUeAiy3b+Tffc9eO9LOkcChhSDHwRpAsMHskrwfeMPsNRHQBzwFnWrYzKc5e3hH/D7gFODs6YZj3BRjY1tWUk4f+X+oV5pFUiSlRfYxJXCNOGKZnA/CbmO8JU1r8LLPFpxKGDFDCkGIw7L0SzGOkSeZl7Nah0Q40j3P24nvuZst2ngUsy3b2idr+9hxgq9ntMPJ5HwV+CNhAfZy3S9fP7Ccj4SW45sA0fZYUCCUMKQYfJHHNYCI76b0CfDPBda8N8T4PAKcAXwQejHocdWdkS1zLdqpNj+VgM4bwpElSu4F/BM4cwfcRK/J9tQyW6BIclyKlhCGyp/KY15HeQLnvuatG8L6PAv9uHvU8GO9xlNmW9GDg7uhHU/Qnk0tG8Nkk+L42+J77wgjfW4qEZklJsYr8qz72F+nk6Bdm7OOvwBTLdvZ6RBTvWDxmG9cngdPN2MTZwHrfc38XddlB5uvq6Hst2xljBu6T+p7iJIe9vi8zWA/w2dgLLduptWynMonPkyKjhCHFKvK45eMxx+P9S/4h80v4yuiDZnxjrWU7y5L8zAeBajM998g402QjMcX+cv8nYKz576oE7x/3e7Js57SocZiIh83Xy+Mkh5uAdy3bOXjob0mKiR5JSbF6DDgf+HfLdr4P9AFfM7OSMDOTIv4V+DJwnWU7k8w4w0RgDjAe+GmSn/m4mVn1A/P6/pjznnlUtNCynW7g0WrP8AAAAQtJREFUXeAMkwC+C/wbMM+ynTG+5y6P8/7Pm7Io37Bs5xWzCO8o4GoT84mW7ZT4ntvne+7Llu3cbqbWrrZs5z+BsPm8ZrNw8S/DaFcpYOphSFHyPXcpMNfMRFoG/BLYYZIG0Y91fM99z0xpvRP4gllNfY0Z6D7F99ynkvzMncBvzFqLNb7n/inm/LtmdtQ6s/bjRrOA7yQT3/PA50ztrHjv3w1YwBrgJ2YGlGOSXeSXf/T3tQC41CTL282CvakmiYx0zEQKUElfn6oji4jI0NTDEBGRpChhiIhIUpQwREQkKUoYIiKSFCUMERFJihKGiIgkRQlDRESSooQhIiJJUcIQEZGkKGGIiEhS/j94u2+fzq2pfwAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "bento_obj_id": "140494171138128",
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "hidden_ranges": [],
        "originalKey": "34e976cd-7d09-40d2-8987-aecdefa7c0fd",
        "requestMsgId": "34e976cd-7d09-40d2-8987-aecdefa7c0fd",
        "showInput": false
      },
      "source": [
        "## Look a the lengthscales from the final model\n",
        "\n",
        "As SAASBO places strong priors on the inverse lengthscales, we only expect parameters \n",
        "0 and 1 to be identified as important by the model since the other parameters have no effect.\n",
        "We can confirm that this is the case below as the lengthscales of parameters 0 and 1 are \n",
        "small with all other lengthscales being large."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "code_folding": [],
        "customInput": null,
        "executionStartTime": 1650234951551,
        "executionStopTime": 1650234951627,
        "hidden_ranges": [],
        "originalKey": "f08ad85a-4d54-4958-a30c-62c2e77d39ff",
        "requestMsgId": "33147b57-ea6b-4c67-9c7d-796bb54d5c84",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        "median_lengthscales = gp.median_lengthscale\n",
        "for i in median_lengthscales.argsort()[:10]:\n",
        "    print(f\"Parameter {i:2}) Median lengthscale = {median_lengthscales[i].item():.2e}\")"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Parameter  0) Median lengthscale = 7.63e-01\nParameter  1) Median lengthscale = 2.46e+00\nParameter 14) Median lengthscale = 1.52e+03\nParameter 44) Median lengthscale = 1.54e+03\nParameter 46) Median lengthscale = 1.59e+03\nParameter 48) Median lengthscale = 1.65e+03\nParameter 17) Median lengthscale = 1.70e+03\nParameter 25) Median lengthscale = 1.71e+03\nParameter 42) Median lengthscale = 1.72e+03\nParameter 45) Median lengthscale = 1.74e+03\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "customInput": null,
        "executionStartTime": 1650234951631,
        "executionStopTime": 1650234951646,
        "originalKey": "1903180f-c926-4f00-98a6-dc1b38d8d8fc",
        "requestMsgId": "6b419c9d-6e7f-454b-939d-e099c9715781",
        "showInput": true,
        "collapsed": false
      },
      "source": [
        ""
      ],
      "execution_count": 16,
      "outputs": []
    }
  ]
}