tutorials/fit_model_with_torch_optimizer.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Fitting models in BoTorch with a torch.optim.Optimizer\n",
        "\n",
        "BoTorch provides a convenient `botorch.fit.fit_gpytorch_mll` function with sensible defaults that work on most basic models, including those that botorch ships with. Internally, this function uses L-BFGS-B to fit the parameters. However, in more advanced use cases you may need or want to implement your own model fitting logic.\n",
        "\n",
        "This tutorial allows you to customize model fitting to your needs using the familiar PyTorch-style model fitting loop.\n",
        "\n",
        "This tutorial is adapted from GPyTorch's [Simple GP Regression Tutorial](https://github.com/cornellius-gp/gpytorch/blob/master/examples/01_Exact_GPs/Simple_GP_Regression.ipynb) and has very few changes because the out-of-the box models that BoTorch provides are GPyTorch models; in fact, they are proper subclasses that add the `botorch.models.Model` API functions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [],
      "source": [
        "import math\n",
        "\n",
        "import torch\n",
        "\n",
        "# use a GPU if available\n",
        "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
        "dtype = torch.float"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Set up function to model\n",
        "In this tutorial we will model a simple sinusoidal function with i.i.d. Gaussian noise:\n",
        "\n",
        "$$y = \\sin(2\\pi x) + \\epsilon, ~\\epsilon \\sim \\mathcal N(0, 0.15)$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Initialize training data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "# use regular spaced points on the interval [0, 1]\n",
        "train_X = torch.linspace(0, 1, 15, dtype=dtype, device=device)\n",
        "# training data needs to be explicitly multi-dimensional\n",
        "train_X = train_X.unsqueeze(1)\n",
        "\n",
        "# sample observed values and add some synthetic noise\n",
        "train_Y = torch.sin(train_X * (2 * math.pi)) + 0.15 * torch.randn_like(train_X)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Initialize the model\n",
        "We will model the function using a `SingleTaskGP`, which by default uses a `GaussianLikelihood` and infers the unknown noise level.\n",
        "\n",
        "The default optimizer for the `SingleTaskGP` is L-BFGS-B, which takes as input explicit bounds on the noise parameter. However, the `torch` optimizers don't support parameter bounds as input. To use the `torch` optimizers, then, we'll need to manually register a constraint on the noise level. When registering a constraint, the `softplus` transform is applied by default, enabling us to enforce a lower bound on the noise.\n",
        "\n",
        "**Note**: Without manual registration, the model itself does not apply any constraints, due to the interaction between constraints and transforms. Although the `SingleTaskGP` constructor does in fact define a constraint, the constructor sets `transform=None`, which means that the constraint is not enforced. See the [GPyTorch constraints module](https://github.com/cornellius-gp/gpytorch/blob/master/gpytorch/constraints/constraints.py) for additional information.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "from botorch.models import SingleTaskGP\n",
        "from gpytorch.constraints import GreaterThan\n",
        "\n",
        "\n",
        "model = SingleTaskGP(train_X=train_X, train_Y=train_Y)\n",
        "model.likelihood.noise_covar.register_constraint(\"raw_noise\", GreaterThan(1e-5))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Define marginal log likelihood \n",
        "We will jointly optimize the kernel hyperparameters and the likelihood's noise parameter, by minimizing the negative `gpytorch.mlls.ExactMarginalLogLikelihood` (our loss function)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "from gpytorch.mlls import ExactMarginalLogLikelihood\n",
        "\n",
        "mll = ExactMarginalLogLikelihood(likelihood=model.likelihood, model=model)\n",
        "# set mll and all submodules to the specified dtype and device\n",
        "mll = mll.to(train_X)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Define optimizer and specify parameters to optimize\n",
        "We will use stochastic gradient descent (`torch.optim.SGD`) to optimize the kernel hyperparameters and the noise level. In this example, we will use a simple fixed learning rate of 0.1, but in practice the learning rate may need to be adjusted.\n",
        "\n",
        "Notes:\n",
        "- As the `GaussianLikelihood` module is a of child (submodule) of the `SingleTaskGP` module, `model.parameters()` will also include the noise level of the `GaussianLikelihood`. \n",
        "- A subset of the parameters could be passed to the optimizer to tune those parameters, while leaving the other parameters fixed."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "from torch.optim import SGD\n",
        "\n",
        "optimizer = SGD([{\"params\": model.parameters()}], lr=0.025)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Fit model hyperparameters and noise level\n",
        "Now we are ready to write our optimization loop. We will perform 150 epochs of stochastic gradient descent using our entire training set."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Epoch  10/150 - Loss: 1.966 lengthscale: 0.645 noise: 2.005\n",
            "Epoch  20/150 - Loss: 1.930 lengthscale: 0.599 noise: 1.868\n",
            "Epoch  30/150 - Loss: 1.894 lengthscale: 0.560 noise: 1.730\n",
            "Epoch  40/150 - Loss: 1.857 lengthscale: 0.527 noise: 1.590\n",
            "Epoch  50/150 - Loss: 1.819 lengthscale: 0.497 noise: 1.449\n",
            "Epoch  60/150 - Loss: 1.779 lengthscale: 0.471 noise: 1.310\n",
            "Epoch  70/150 - Loss: 1.737 lengthscale: 0.448 noise: 1.172\n",
            "Epoch  80/150 - Loss: 1.692 lengthscale: 0.427 noise: 1.038\n",
            "Epoch  90/150 - Loss: 1.645 lengthscale: 0.407 noise: 0.908\n",
            "Epoch 100/150 - Loss: 1.595 lengthscale: 0.389 noise: 0.785\n",
            "Epoch 110/150 - Loss: 1.542 lengthscale: 0.372 noise: 0.671\n",
            "Epoch 120/150 - Loss: 1.487 lengthscale: 0.355 noise: 0.566\n",
            "Epoch 130/150 - Loss: 1.429 lengthscale: 0.341 noise: 0.471\n",
            "Epoch 140/150 - Loss: 1.370 lengthscale: 0.328 noise: 0.389\n",
            "Epoch 150/150 - Loss: 1.311 lengthscale: 0.317 noise: 0.318\n"
          ]
        }
      ],
      "source": [
        "NUM_EPOCHS = 150\n",
        "\n",
        "model.train()\n",
        "\n",
        "for epoch in range(NUM_EPOCHS):\n",
        "    # clear gradients\n",
        "    optimizer.zero_grad()\n",
        "    # forward pass through the model to obtain the output MultivariateNormal\n",
        "    output = model(train_X)\n",
        "    # Compute negative marginal log likelihood\n",
        "    loss = -mll(output, model.train_targets)\n",
        "    # back prop gradients\n",
        "    loss.backward()\n",
        "    # print every 10 iterations\n",
        "    if (epoch + 1) % 10 == 0:\n",
        "        print(\n",
        "            f\"Epoch {epoch+1:>3}/{NUM_EPOCHS} - Loss: {loss.item():>4.3f} \"\n",
        "            f\"lengthscale: {model.covar_module.lengthscale.item():>4.3f} \"\n",
        "            f\"noise: {model.likelihood.noise.item():>4.3f}\"\n",
        "        )\n",
        "    optimizer.step()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "#### Compute posterior over test points and plot fit\n",
        "We plot the posterior mean and the 2 standard deviations from the mean.\n",
        "\n",
        "Note: The posterior below is the posterior prediction for the underlying sinusoidal function, i.e., it does not include the observation noise. If we wanted to get the posterior prediction for the observations (including the predicted observation noise), we would instead use `posterior = posterior = model.posterior(test_X, observation_noise=True)`. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "# set model (and likelihood)\n",
        "model.eval()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XdgU+X6wPFvRpN07z0oggyFspEhsgQVC1cEoqj3Il7FhaKouAXxqteJiID6c+IAIigiQ0YRARFkT9l07502zc7vjyIXsGhp056T9P38RZP0nKeHJM8573ne51W4XC4EQRAEQW6UUgcgCIIgCHURCUoQBEGQJZGgBEEQBFkSCUoQBEGQJZGgBEEQBFlSSx0AQFpamiglFARBaMGGDh2quPAxWSQoaoNr9DY2btzIoEGD3BKPtxDHpG7iuNRNHJe6ieNSN3cdl7S0tDofF0N8giAIgiyJBCUIgiDIkkhQgiAIgizJ5h6UIAiCO7hcLkpLS3E6nW7bZkREBEVFRW7bnre41OOiVCoJCwtDofhTPUSdRIISBMGrlJaW4u/vj06nc9s2dTodgYGBbtuet7jU42I2myktLSU8PLxerxdDfIIgeBWn0+nW5CS4j06nu6QrW5GgBEEQBFkSCUoQBEGQpUbdg9Lr9Z2A74FZBoPhvQueSweyAMeZh243GAw5jYpWEATBA5w8eZJHHnmE/Px8HA4H/fv35/XXX8fX15c777yTsWPHkpqaKll8K1asYMmSJXz22WdnH0tPT6dz58706NEDl8uFWq3mmWee+csmCllZWVRVVdG7d+8mibPBCUqv1/sDc4C6pwDXusFgMFQ1dB+CdMw2BzaHE5vDhc1RO2asUipQKRT4qJX4+ahQKutXiSMIcpeXl8ett97K4sWLiYmJadS2nE4nY8aM4a233jr75f7WW28xadIkvvjiCzdF3DTat2/Pxo0b4UySHTlyJIsWLSIlJaXO12/atAm73S6/BAVYgBHAk26MR2hGTqeLMpOVQqOFQqOFsmorRrONSrMdq/2vb2QqFOCvUeOvVRPq50NEoJaIAC2RgVoCtKI4VPAsL730Elu2bGHmzJnMmzevUdtau3Yt7dq1O+/KY+rUqbRv357CwkIAfvjhB2bNmkVxcTGffvopnTt35o477iAvLw+LxcKLL77I9ddfz7x58/jqq69QKpXcdNNNPPbYY8yYMYNTp05x+vRpQkJCeOKJJ7jmmmuoqamhY8eOnDx5khdeeIHNmzfjcDiYPHky48eP58CBA/zrX/8iLCyMNm3a/O3f0aZNG5599lnmzp3LBx98wNSpU/ntt98wm83cd999/OMf/+DVV19Fo9GQlJSEn58fzz//PBqNhtDQUAwGAxqNplHHssHfJAaDwQ7Y9Xr9X73sfb1enwxsAZ42GAwXbQr7R9ZujKqqKrdsx5tceExO5xbz6sszGfPgc1i1IfxNHmqQAB8FEX4KInyVRPkp0Krkd6Ul3it184bjEhERUe8qvsjISCwWy9mf58+fz/z589FqtefN73E4HBiNxnptc+/evVxxxRV/en2HDh3Yt28fNpsNlUrFsmXLWL16NS+++CKPP/44BQUFrFy5kvLyctauXcuBAwdYtGgRq1evBmDYsGHccMMNWCwWqqurWbVqFV9//TVLly6lW7durFq1iiFDhrB+/XpOnDjBypUrsVgsDBgwgKFDhzJ9+nSefPJJbrzxRh599FFsNtt5MVZVVeF0Os97rGPHjsydO5eioiJiYmJYvXo1NTU1dOnShVtuuYXx48cTERHB4MGD+e677/jwww9JTk5m0qRJLFu2jBtuuOFPx6egoIBDhw7V61g25anuC8CPQCmwDBgDLLnYi93RcFA0dPyzjRs30qFbH44VGDleWMUnCz/kxO8H2bR6GWMfntGk+y4GSl0K4gN9aRsVwOVRAfjL5OpKvFfq5g3HpaioqN5zc06fPs3jjz/OsmXLMJlM+Pn5MXr0aN58883ztmE0Guu9TZ1Oh91u/9PrVSoVgYGB+Pj4MHz4cAIDAxk0aBAzZ86kR48emEwmHnjgAUaPHs3EiRP55ptvOHXqFKNGjQLAZDJRXFyMVqulf//+BAYGcssttzBgwAACAwNZu3Yt48ePZ/fu3ezatYuRI0ee3XdVVRXHjh1j6NChBAYGMmzYMFavXn1ejAEBASiVyvMeczqdaDQaIiMjMZlMXHfddWg0GoqLiwkMDEShUJydC5WUlMQjjzyC3W7n1KlTXHfddXUes+joaDp16nTeYxdrFttk3xYGg2HBH//W6/WrgM5/laAE9zJZ7RzKrWRduo09tkympaZgt/7vTHHrioVsXbEQtUbL6yv2N1kcTpeLrFITWaUmfj5aRHKEH53ig2kd7i/uYQmSi42NJSgoCLPZjE6nw2w2ExQU1Kj7UB07dmT+/PnnPeZyuTh06BDt2rUDOK+TgkKhwM/Pj23btrF161Y+++wzVqxYwciRI7nxxhv54IMPztvWhg0bzg6dhYSEEBcXx5EjR/j111/54IMPOHjwIP/+9795+umn/xSDUllbuF3fuUg7d+6kW7du/Pzzz2zYsIGff/4ZHx+fOhPPXXfdxcqVK+nYsSOTJ0+u9/H6K01SZq7X64P1ev0avV7/xwDkQOBgU+xLOF9ueQ0r9+fx0ebTbDleTJWtdlT1uc/X031wKj7a2qEPH62O7kNG8tyCv6pxcS+ny8WpomqW783lk19O89vpUsw2Rz1+UxCaTkFBAffddx/btm3jvvvuIz8/v1HbGzZsGKdPn2bVqlVnH5s1axYDBgwgLCwMgC1btgCwbds2OnbsyO7du/n666+5+uqrmT9/PocPH6ZHjx789NNPmEwmXC4XU6ZMoaam5k/7Gz16NK+88gp9+/ZFrVZz1VVX8cMPP+B0OjGbzTz00ENwpgBi586dAPz0009/+3ecPHmSt99+m0cffZTi4mISExPx8fFh+fLl2O12rFYrSqUSu90OQEVFBUlJSZSXl/PTTz9htVobdRxpZBVfD+AtIBmw6fX6scBy4LTBYPjuzFXTNr1eXwPsAZY2OlqhTi6Xi1PF1exKLyOn/M9vYICg8Ch0fgHYrRbUGi12qwWdXwBBYZHNHi+A0WznlxPF7EgvpXN8MN2SQgjU+UgSi9Cyffvtt2f/PXfu3EZvT6lUsmbNGu677z5eeOEFnE4nPXv25N133z37GpfLxciRI8nKyuKLL74gISGBZ555hg8++ACVSsUTTzxxdsjsmmuuQaVScdNNN+Hr6/un/Y0ePZqHH36Y77//HoB+/foxePBg+vbti8vl4oEHHgDgueee46677uLdd9+ldevWdSaQo0ePMmjQICwWCw6Hg7lz55KUlERwcDCvvfYaAwcO5KabbiI1NZX777+fUaNGcf/99xMZGcmDDz5I//79adeuHdOmTWPGjBmMHDmS2NjYBh9Lhcsl/WK2aWlpLrFgYcOcKDTy68kSiqvqPltJT08nOTkZgE9fnExQWCR9RtzCtlWLqSwtYuL09+r8veamVipISQyhd3IYvhpVk++vJb5X6sMbjktRURGRke498bqUe1AtSUOOS13/P2lpafJeUVe4NBkl1Ww9WUJ+hbnev3NuMhrz0PQmiqxh7E4XuzPKOJhTQY9WofRoFYqPSjQ6EYSWTCQoD1NcZWHTsSIySkxSh9IkrHYnv54s4WBOBQPbRXJ5tDhrFYSWSiQoD1FjdfDrqWIOZFfilMGwbFMzmu2s2J9Hq/AKBrWPIsy/cRP+BEHwPCJByZzL5eJgTiVbThS3yIq3jBITX23LoE+bcHokhYrSdEFoQUSCkrFCo5mfjhSSW17/+0zeyO50seV4MccLqhh+ZTQRAVqpQxIEoRmIBCVDdoeTbadK2ZVR1iKG8+qroNLM19sz6d82nO5JofVeNloQBM8kyqRkJrvMxJfbMtiRXiqSUx0cThebjhXz3Z4cqi12qcMRhDqlp6ejUCjYtm3beY/37NmTO++8U7K4PI1IUDJhtTvZcKSAJbuyKTPZpA5H9jJKahN5enG11KEIQp0uu+wyFi5cePbnEydOUF5eLmlMnkYkKBnIKq39st2XVYG4aKo/k9XBsr05/HqyBDlMOBeEc/Xp04d169bhcNQWNy1atIjhw4cDsHnzZgYMGMCQIUOYMGECVqsVu93O7bffzsCBA+nZsycrVqyAM420X375ZYYOHUqXLl3IzMyU9O9qTiJBScjmcLLxaCFLd2dTUSOumhrC5YJtp0pYvi+3RVY5CvLl4+PDVVdddbbv3ffff8+IESMAzrYm2rBhA9HR0XzzzTeUlpYyfPhwfv75ZwwGA9On/28yfVBQEGlpadxwww3ntWbydqJIQiIFlWZ+PJhPaXXjGyoKcKqomq+3ZzKqa5yo8hPO455amvMnjNf3gn3cuHEsXLiQ2NhY4uPjCQgIoKCggOPHj3PzzTcDUF1dTUREBKGhoezYsYMPP/wQpVJJSUnJ2e0MGDAAgISEhPMe93YiQTUzp9PF9tOl/HZaFEG4W0WNjcU7shjROZbWEf5ShyPIhDs+Zg3txTds2DAmT55MbGwsY8eOBUCj0RAfH/+nhSE///xzSktL2bx5M6WlpfTs2fPsc2r1/76qW9Jwthjia0Zl1VYW78xi26kSkZyaiNXuZPneXHZnlkkdiiDg4+PDNddcw8cff3x2AcHQ0FAADh8+DMCcOXPYv38/xcXFtG7dGqVSybfffuuW5So8nbiCaib7ssrZfLwIm0MkpqbmdLn4+WgR5SYrg9pFie4TgqTGjRtHUVERwcHBZx/7+OOPmThxIhqNhri4OCZNmkRQUBCjRo1i27Zt3HXXXSQkJPDSSy9JGrvUxHIbTazaYmfd4QJOS1QOfe5yGy1Rm6gAbugU86fO6HJ8r8iBNxwXsdxG82nq5TbEEF8TOlFo5IttGZIlJwFOFlbx7e5sUeEnCB5IJKgmYLE7+PFgPj/sy6PGKr4YpZZbbsawMwujWZTyC4InEQnKzTJLTHzxawa/51VKHYpwjpIqK4ad2VSILh2C4DFEgnITm8PJT0cK+XZPNkaz6BEnR5U1Ngw7syipskgdiiAI9SCq+Nwgu8zEusMFlIuzc9mrstj5Zlc2MWan1KEIgvA3RIJqBKvdyZYTRezPFj30PEmN1cEvuXb6VpqJDtJJHY4gCBchElQDpRdXk3akkErRQ88jWR2wdHc2Y7oniCTl5WatO9bobVitFjSa2hZajw5rV6/fOX78OI888ghFRUU4HA769evHm2++iVZb/1ZckydPZuvWrbzzzjukpaXx4osvnvf82LFjmTx5ssdPDbgYcQ/qElVZ7Kzcn8d3e3JEcvJwFpuTpbuzya9o2SsWC+7ncDgYM2YM06ZN47fffmPnzp0AzJw585K2s2rVKjZs2MA111zzp+TUEogrqHpyOl3sz6nglxPFWO3i/oW3sNicfLsnm7HdE4gSV1KCm6xbt44OHTowcOBAABQKBa+//jpKpZLZs2ezaNEiAG666SaefPJJ7rzzTuLi4ti1axeZmZl89dVXpKWlkZuby8iRI3n88cf54osvWLJkCa+//joLFy6kVatWVFbWVgsbjUYmTpxIWVkZdrudOXPmkJKSQtu2bbn33nv54YcfsFgsrF+/Hp1Ox4QJE8jIyECn07FgwQJiYmKYNGkSp06dwmazMXPmTIYMGSLpMUQkqPrJLDHx8/Eiio2i+ssb1SapHMb2SBCd0AW3OHLkCF27dj3vMV9fX06fPs1nn33Gjh07AOjdu/fZJrIWi4U1a9bw/vvvs2DBAt555x3mzp3L6tWrz16BlZeXM2/ePI4cOYLNZqNNmzYAvPPOO1x//fXcfffdHD58mClTprBu3TrsdjsdOnTgiSee4NZbbyUtLY3i4mJiYmL4+uuvWbRoEcuXLycgIIDY2Fg+/vhjiouLGTJkCPv372/243YhkaD+Qlm1lU3HizhVJDpBeLsaq4Nvd2czrkciof4aqcMRvMAfCxWea8+ePfTp0+dsd/L+/fuzb98+uGBJje3bt9e5zRMnTnDllVei0+nQ6XT06NEDgK1bt1JUVMSXX34JgMlkOvs75263oqKC3bt380druVtvvRWA+++/n82bN7NlyxYAampqsFqtaDTSfhZEgqpDpdnG9lOlHM6tFF3HW5Bqi4Olu7MZ1zORYF8fqcMRPFjHjh157733znvMYrFw6NCh85bLsFqtKJW1pQD1WVLD5XKdfT2A01l7u0Gj0TBnzhz69u37p9+5cLsqlers7/1Bo9Hw7LPPMn78+Ab8tU1HFEmcw2i28dORQj77JZ2DORUiObVARrOdb3dnU20Rk62Fhhs2bBgZGRn88MMPcCaRPPnkkxw7doxff/0Vu92O3W5n+/btdOvWrd7bbdOmDb///js2m43Kykp27doFwFVXXcWyZcvgzDIeb7/99kW30atXLzZs2ADAihUreOWVV877/cLCQp555plG/f3uIq6ggEKjmd0ZZRwrqMLhFEmppSs32fjuzD0pnY9K6nCERqpvWfhfudSu3UqlkjVr1nDvvffy4osvotFoGDZsGG+//Tbz589n4MCBOJ1O7r77blq1alXv7YaFhTFhwgT69OnDZZddRq9evQB46KGHuPPOOxkwYAAOh4N33333otu49dZbWb9+PQMHDkStVrNgwQKio6PZsGED/fr1w+FwMGPGjHrH1JQatdyGXq/vBHwPzDIYDO9d8Ny1wCuAA1hlMBguurCJFMtt2B1OjhdWcSi3kqxSUz1+wzO19OU2LqY+xyU+1Jebu8WjVrWcgQax3EbdxHIbdZPtcht6vd4fmAOkXeQl7wJjgP7AcL1ef0VD9+UuTqeL7DITG44U8OHmU/x4MN+rk5PQODllNaw8kIdTXFULgiQaM8RnAUYAT174hF6vvwwoNRgMWWd+XgUMBQ43KtoGMFnt5JbXcLKomtPF1WL5C+GSnCqqZsORQq69IlrqUAShxWlwgjIYDHbArtfr63o6Big65+dCoE1D91UfLpcLk81FRkk1ZSYbhZVmcstrKBMNXIVGOpBTgb9WTd824VKHIggtSlMVSVw4lqgA/nKcZOPGjY3a4YlyB7tyaliT/kujtuNtrFYr6enpUochO5d6XNLT0zl6UEVysHcXTVRVVTX6syg1f39/AHQ693UGcTgcGI1Gt23PW1zqcTGbzaSnp3Po0KF6vb6pElTOmauoP8QDeX/1C429MRuYUcaBoh2iIOACl1ok4XKByaikokSNxaTEWqPEYj4zT8PHhVrjRKN1ERRmJyjMgVojzf2ZypJCFrwylX89O4ugsEu/Id6Q4pEKhYJWneJoHeF/yfvzFN5QJOFyuSgtLcVsdl+PxYKCAqKjxTDvhS71uCiVSnr27IlCcf41TFpa3aUMTZKgDAZDul6vD9Lr9clANpAK3N4U+xIarqJERfYxHdknteSc0FGU40N5oQ8KpYvgCDu+/k40WicanQsUYLcpsFsVWM0KKkvVGMvU6PwdRMTZiE22ENvaSnwbC4ntzfg0ceJa+9U8Th/cydov5zL24eYpiXW6XKw6kMe4HqJvn5wpFArCw907HHvo0CE6derk1m16g6Y+Lg1OUHq9vgfwFpAM2PR6/VhgOXDaYDB8B9wPLDzz8sUGg6HxPe+FRjFXKzm2x49ju/04vscPk1FFwuVm4tta6D64kqhEK6HRtYmpPpxOqK5QUZilIS9dQ95pLTvTgijI0JDYzkyblBo69Komqb0ZpZsqtaelpmC3/q8n4tYVC9m6YiFqjZbXVzR97zCr3cn3e3O5tXcigTrRbUIQmlJjiiR2ARcdCzAYDJuAP/fdEJqVpUbFrrRA9m4K5OQ+X1p1NNOuu4m+N+YS29raqMShVEJgqIPA0BrapNScfdxsUnD6kC8n9vmx+O1oaqpUdO5fRcrVRtqk1DRqn899vp7lH77Gga3rsVnM+Gh1dO4/jFGT/lRM2mSqLHaW7c1F3zMBrdq770kJgpREJwkv5HJBxu86fl0ZzL4trWmbYqHrQCO3Tcuv99VRY+j8XHTsZaJjLxMj7y6mINOHA78EsOz9KCwmBVddX0mv4ZWERFx6O6Gg8Ch0fgHYrRbUGi12qwWdX0CD7kM1RrHRwuoD+YzqEodS+af5hYIguIFIUF7EblWwMy2QLd+HYDUr6XtjBV2v303HzvGSxhWdZCM6qYyht5aRfVzLttXBvHFvKy7vYmLwuDJadby0m9nG8hL6pY6nz4hb2LZqMZWlRfX4Lfc7XVzNz8eLGNw+SpL9C4K3EwlKRhpamWY2Kfh1ZQg/fxtK3GUWRk0qpm1XE0olpKc3/zywi/0dCgUktrOQ2K6QUZOK2P5jMAteiSU0ysaQW8ro2KsaRT0uRiZO/19XrTEPTW+qP6Ne9maWE+anoUtiiKRxCII3ajlNxjzAuZVp9WE1K0hbHMrLE1qTfVzLPf/JYdLLObTrbnJbUUJD1Ofv0Pq6uGZ0Oc98dpp+qRWs+CiCOY8mcmKfb7PG6g4bjxaRUSLWDBMEdxNXUDJwqZVpDjts/zGYtV+F0/qKGia/lUV0kvQdMxpSYadSQffBRrpeY2TPxkAWvx1NRJyN1LuLiG9jbcboG87pcrHyQB639koiTCx2KAhuI66gZOC5z9fTfXAqPtrauTU+Wh3dh4zkuQV/nrx2dJcfb9zbin2bA7hrRg4Tns+TRXLiEv+OCylV0GOokSc/SqdT3yo+eDqBpe9FUV3pGW9Ri83J8r05mG2i16MguIu4gpKB+lSmleSpWf5hJLmntPzjviKu7FO/+zXNyR0Vdmof6D+qgq6DjKz+LILX70nmhjuL6X1dpaTDlvVRZrKx6kAeN3WNF5V9guAGIkHJxMUq0xx22LgklI1Lwrjm5jLueDq/ybs0NIa7Kuz8g5yMfbiQPiMqWDI7it0bgrhlagHhsfK4WryYjBITm44XMUhU9glCo4kEJRN1VaZlHdOy+O1oAkMdPDIng/BY+S9D7u4Ku4S2Fh56J4tN34byzkNJDLu9hKtHlaOU8fzYPZnlRARo6RQfLHUoguDRRIKSIZtVwY8Lwtm5LoiR9xTRY6hRdsN5zUmlgsHjyriybxWL347mwNYAbp+WT0ikfBP2hiOFhPlriAvxvKpEQZALmY/qtzw5J7XMmpxEcY4Pj7+fQc9rW3ZyOldUgo0H38imXTcTbz+YxP4tAVKHdFEOp4uV+/Oossg3iQqC3IkrKJlwOmCDIYxN34Yw6l5x1XQxShUMu62Uy7uZ+PLVGI7s9GP0A0WyvC9XZbGzYl8uY3skoFaJc0FBuFTiUyMD5UVq5j+ZwNHdfjw6N1NcNdVDckczj8/PpKZKxZxHEynNl+e5Vl6FmQ1HCqUOQxA8kkhQEju41Z+3H0yifQ8T9/83m9AoMSRUXzp/J/96No8eQyqZPSWJIzv9pA6pTodyK9mbVS51GILgceR52tkC2G2w/MNIDm8L4K4ZuSRf4b7VP1sShQIGjikn4XILX7way4B/lDHkljLZXYFuOlZERICGhFB5JlFBkCNxBSWB0gI1701NpLzIh8fmZ4jk5AZtUmp4ZE4m+zYHsvCNGOxWeWUoh7N2NV6jWd7zuARBTkSCamaHf/PnnYeS6DqwionTc/ENaPr1mVqKkAg7k9/KwmZRMG9aAsYyeU2WqrY4WLk/D4dTfgUdgiBHIkE1E6cTflwQzjfvRHHnC7kMGiu/YShvoNG5+OezeVzezcTsKYkUZslrWfa8CjM/iaIJQagXkaCaQU2Vkk+mx3Firy+PvpfJZZ3EkF5TUirhhgklDL+jlLmPJ5J+WCd1SOc5kFPBwZwKqcMQBNkTCaqJ5adrmDU5ibAYG/e/nk1QmOh23Vx6D6/k1sfy+WR6HAd/9Zc6nPP8dKSQ/ApxoiIIf0UkqCZ0cKs/855IYNjtJdz8YBEqUTPZ7Dr2NnH3f3JYMjuabauDpA7nLLvTxYr9udRYxQmLIFyMSFBNwOWCtV+GsfS9KP79Ui69hhmlDqlFS2pv4cG3slj3dTg/fyufpdmNZjurDuThFEUTglAnkaDczFKjYMHLsRz+zZ9H5mTSqoMYxpGDyHgbk9/KYusPIaz5MgyXTHJCZqmJX04WSx2GIMiSSFBuVF6k5r3HEvHRuHjwzWyCw8XwjZyERtmZ/HYW+zcHsmVJO9kkqZ3pZZwoFFfZgnAhkaDcJOOIjtlTEuk60Mj4J+S9qGBLFhjq4IE3ssg5Fsqy9yP/MklVlhTy3mN3NHjRxUux5lABpdXWJt+PIHgSkaDcYM/GAD56Po4xkwsZKsM2O8L5/IOcjH50F+mHfVk2/+JJau1X8zh9cCdrv5zb5DFZ7U5W7M/FahcTtwXhDyJBNYLLBWu+DGPFR5Hc999sOvWrljokoZ60fnbuezWbjN91fDfv/CQ1LTWFqcPbs3XFQlwuF1tXLGTq8PZMS01p0phKqqysO1zQpPsQBE8iElQD2awKvvpvDL9v92fKu5nEtxHDM57GN8DJvf/NIfOo7rwrqec+X0/3wan4aGsn+PpodXQfMpLnFqQ1eUzHCozsyihr8v0IgicQCaoBjGUq5k9LwOFQ8MCbYvKtJ/P1d3LvKzmcOujLyk8icLkgKDwKnV8AdqsFtUaL3WpB5xdAUFhks8S05Xgx2WWmZtmXIMiZSFCXKD9Dw+wpibTtYuKfz+Sh0YpiCE/nG+Dk3lezObzNn/VfhwFgLC+hX+p4psw20C91PMay5isFd7pqO5+L5eKFlq7BvQ30ev0soA/gAqYYDIYd5zyXDmQBf1xa3G4wGHLcErGE9vxk5avXIxh1bzrX3CR1NII7BQQ7ue+1bOY+loiP1sXE6e+dfW7MQ9ObPZ5qi4NV+/MY2yMBpVJU3QgtU4MSlF6vHwhcbjAY+ur1+o7AJ0DfC152g8FgqHJPmNL7dWUwy97X4nT8g8LMeGCG1CEJbhYU5uC+17J577FEdP4O+txQKWk8OeU1bDpexKD2UZLGIQhSaegQ31BgGYDBYPgdCNXr9fJpdOZGTgc8dv0qvpldjs3SC9jUbFVdQvMLjbJz76vZ/Ph5BPu3BEgdDnsyyzmSL22iFASpNHT6FGiLAAAgAElEQVSILwbYdc7PRWceO/eT9L5er08GtgBPGwyGv7xZs3HjxgaGUutEmQOr1Up6enqjtnMuq1nFmo87E9N6PH7Bk8g4mIXdCmqNljZd+3P1uHvcur+m4O5j4i3+7riMuD+fxbN6UFGVQ2IHaavq/i8znYEJPgRpm36or6qqqtGfRW8kjkvdmvq4NDRBXfhJUZy5F/WHF4AfgdIzV1pjgCV/tcFBgwY1MJRagRllHCjeQXJycqO284fyIjUfvx5HfBsLYx8uYNk8P07tsaLWaHHYrIRHRXNll+5u2VdTSk9Pd9sx8SZ/d1ySkyEspIAF/+nGPS/nkNjO0qzxXajUz4frrkpCq27aVYI3btzY6M+iNxLHpW7uOi5paXVP4WjoEF/OmSumP8QB+X/8YDAYFhgMhkKDwWAHVgGdG7gfSWQd055tW3TL1ALUPtJWdQnSaNulhnGPFPDxC/GU5Em7Mm+ZycaaQwW45NJAUBCaQUOvoNYCLwIf6PX6bkCuwWAwUltAEQwYgJEGg8EKDPy7qyc52bc5gCXvRjFuSiEpV/+vxkPqqi5BGp37V1NZWsKHz8Tz0KwsAkKkm/N2srCKHell9G4dJlkMgtCcGnQFZTAYtgK79Hr9VmAO8KBer79Tr9ePNhgMFWeumrbp9fpfztyfWur+0N3L5YL1C0P5fn4k976Sc15yElq2/iMr6HKNkY9eiMNqlrbke+vJYjJKREstoWVo8Dwog8Hw1AUP7TvnudnA7EZF1oxsVgWGWdEUZGp4+N0sQiLEBEnhfDfcWUJ5kZovXoll4vRclE17K+iiXC5YfTCf8b2TCPaVdthREJpai+8kYSxTMe+JBBx2BZPfEslJ66MkIlBL6wh/2scEcmVcEF2TQuiWFEJKQjBXxAXRLjqQ+FBfgn198FG1jEmkCgXcMrUAm0Xxt8t0NLUaq4MV+3OxO0Tnc0E6VRY7VkfTfhAafAXlDXJOavlkRhy9h1cw/I7SFrVMhlKhICpIS1SglqhAHVFBWkL8fBpUJWY02yipslJcZaHIaCGnvAaj2fsSvUoNE17IY84jiWz6LoSBN5dLFkthpYW0I4Vcd2VMPV4tCO5ldzj5YV8uQU18jtRiE9TeTQEsnRPFmIcK6XpNy7jfFOTrw2UR/iSG+ZEQ6ovOxz3jVIE6HwJ1PiRH+J99rMJkI6vMRHpJNRklJq9Z58jX38nd/8nh3UeSCI+xSbrEyuHcSqKDdHRNDJEsBqFlSjtSSH6FmabuztDiEpTTCWsWhLMzLYh7X80hoa2081uamq9aQfdWobSPDiQmWNds+w328yHYL5hO8cHYHE4ySkycKDRyorAKWxMPCzS1sGg7d83I4f+eiyc4Qto5UpuOFREZqCU+xFeyGISWZXdmGYdzm6e7SYtKUDXVSr5+PQaTUckj72YSGOqdy2QoFNA6wp/O8cFkqHMY2K55lom4GB+VkrZRAbSNCsBsc3A038jB3AoKKz335CCpvYVxUwr5ZEYcU2ZnERIpzZCmw+li5f5cbruqFQHaFvVxFiSQWWJi87HmmwPaYt7RBZk+fDojnrZdTUx4rhC1FxZAadRKOsUH0zUx5GyFV6bMbqzpfFR0SQyhS2IIueU17Mwo41RRlaRFBw2VcnUVRTk+fDw9jslvZaH1leaPqLY4WLEvl3E9E1GJzudCEymrtrLyQB7OZvywtogEdeAXfwzvRHPjXcWSd6huCrVf+sF0Twp1232l5hAX4suoEF/Kqq3szKgdNmjON787DNGXUZip4es3YpjwXB5Kiepi8yrMbDhSyLAroqUJQPBqZpuD5ftyMduad9TJq8vMHQ5Y+Uk4382N4u6ZuV6XnHxUCq5qHcZdVyfTr02ERyWnc4X6axh2RTT/6tuK9jGBHlVNqVDAuCmFVJerWP1ZuKSxHMypYG+WdJWFgndyOmsX0Cyttjb7vr32CspYpuLLV2NB4eLRud51v0mpUHBlXBB92oR71X2HUH8NIzrH0jM5lE3Hiskq9Yxlz9UaF3dOz2X2w0nEtLLSY6hRslg2HSsiIkBDQqifZDEI3uXnY0VklEjzWfTKK6jTh3TMejCJVh1ruPeVHK9KTnEhOsZflci1V0R7VXI6V1SgjrE9ErgxJZZAnWf8jQHBTu6akcuy9yPJPKqVLI7aook8Ks02yWIQvMeezDJJr8q9KkG5nLDBEMqnL8Yx5uFCRkwskawljbv5alQMuyIafc9EogKbr1xcSu2iA5nQL5nercNQesC4X2xrK7c8WsCnL8ZRUSLdG89kdfDDvlxsotOE0AiniqrY1IwVe3XxmgRVVqJg+XvdOLg1gEffy+TKPt7TULN9TCAT+ibTKT4YhQd8UbuTj0pJ/7YRjO+dSGSgdFcm9dWpXzX9Uiv4ZEY8Vot0/1eFlRbWHiqQbP+CZys0mll9MF/yoiWvSFBWK9w1JpCI+CoefDOL0CjvaLPjr1UxskssIzrH4qvxkkvBBooK0jG+dxL92oTLvpT62vGlhMdYWTI7WtLy+WMFRradKpEuAMEjGc02lu/NlUX3F69IUBoNfLDISP8xx1F5xi2Lv9U2KoB/9kmmbVSg1KHIhkqp4KrLwrmlVyKhfvKdyPZHY9ncUxo2fSdtG6Jtp0o4UShd0YbgWSx2B8v25sqml6ZXJCiAqBjPmj9zMT4qBdd2jGZkl7gWf9V0MdFBOm67qhWd4oOlDuWitL4u7pqRywZDGMf2SNeGyOWCNYcKKDSaJYtB8AwOp4sf9uVRbJRPhxevSVDeICJQy/jeSXROkO8Xr1xo1EqGXRHNjSmxaNTyfBuHxdi546k8vvpv7EWXjK8sKeS9x+6gsrSoyeKw2p0s35tLlUUeZ8WC/LhcLtYeypfd1A55frJboI6xQdzaK5HwAPkXAshJu+hAbu2VSJi/RupQ6nR51xqG3lrKJzPisNT8+d7Z2q/mcfrgTtZ+ObdJ4zCa7aKyT7ioLSeKOZIvv6FgkaAkplIqGNIhius7xeCjEv8dDREeUHvl2S5anvfrBtxUTnwbM4ZZ/yuamJaawtTh7dm6YiEul4utKxYydXh7pqWmNFkc+RVm1h4qwOVh7aSEprUro5Sd6WVSh1En8Y0ooQCtmnE9E+gi1vNpNI1ayY0psfRvGyG7VkkKBYydUkhRtoafl9b+Xz/3+Xq6D07FR1s7p81Hq6P7kJE8tyCtSWM5VmBk60lR2SfUOphTIflcp78iEpREooN0jL8qidhgsY6PO/VuHUZqSqzslqLXaGvbIf30TRjH9/gSFB6Fzi8Au9WCWqPFbrWg8wsgKKzpl0b57XQpB3Mqmnw/grydKKwi7fdCqcP4SyJBSaBddCDjeiZ4basiqbWNCkTfM1F2bZLComuLJr58LZbSAjXG8hL6pY5nymwD/VLHYyxrvjPZtN8LySjxnsnswqXJKjWxupmXzmgIeX2CW4CrLguj72XhLa4jRHOLCtJxS69Elu3NlVXZ7OXdahg8tozPZsbx0Ky5+GhqvyDGPDS9WeNwulys2J+HvqdndOgQ3CenvIbl+3KxO+WdnBBXUM1HqVAw7Ipo+rWJEMmpmQTqfBjXI4GEUHkNow4cU0ZEnI2lc6Ik7TRhtTv5fm8ORtFYtsXIrzCzbE+OLLpE1IdIUM1Ao1YyqmucrCeWeiudj4rR3eJlVeFX22kin8wjOratkvY9YTTbWbYnp9kXohOaX6HRzHcelJwQCarp+WpUjOmeQOsIf6lDabHUKiUjOsfQJVE+Jwha39qiidWfhZPxu7Td6YurrLVDPmKOlNcqMlr4brfnnYiIBNWEAnVqxvVIICa4ZSyPIWcKhYIhHaLplRwmdShnRSXY0E8t4PP/xGIsk7atVU5ZDasP5os5Ul6o0Ghm6e5sTFbPSk6IBNV0Qvx8GNdTdIaQm6svj6B/2wipwzirU99qel5byRevxOKQ+PvjRGEVe4s870tMuLjCSjNLd+VQ44HJCZGgmkZEgAZ9z0SCfeXbcbsl6906jM6R8mnEe/2/SlCqXKz6VPrEmV7hZPPxpusLKDSfgkozSz1wWO9cIkG5WWSglrE9EvEXc5xkrW2IikHtm35SbH0oVXDH03ns3RjIvs0BUofDzvQyfjtdKnUYQiNkl5lYsivbo5MTjZkHpdfrZwF9ABcwxWAw7DjnuWuBVwAHsMpgMLzktohlLCpIy5juCeh85HN2Llxct6RQFAoFPx2RfjZ9QLCTCc/n8n/PxRPTykJ0krSl37+cKEarVoo2XB7odHE1K/fnYnN4/v3EBl1B6fX6gcDlBoOhL/Bv4N0LXvIuMAboDwzX6/VXuCdc+YoO0onk5IG6JoYwpEOU1GEAkNTewoiJJXw2s+7O583tp6OFoiWShzlWYDzTtd7zkxONGOIbCiwDMBgMvwOher0+iNrkdRlQajAYsgwGgxNYdeb1XisqSMvN3eNFcvJQXRJDGCiT4b4+N1TQqqOZxW9Lu1w8ZxY7XP97AYdzK6UNRKiXvVnlrDqQh6OZOkRUlhTy1GOPkJ+f32T7aGiCigHOvZNadOaxup4rBGIbEaOsRQSKYT1v0D0pVBbVfQoF3Dy5kOIc6ZeL50ySWne4gCP5IknJlcvlYvPxIn46UtisJzVrv5rHoYMHmDlzZpPto6H3oC4cf1CcuRf1d89d1MaNGxsYSq0TZQ6sVivp6emN2s6lCNQouCJezbZfTjfbPi9VVVVVo4+tN7rYcdFWOjhaKv2N5aF35WN4tTfqwHTiLy9vtv1e7DP0QfppekSrSQhsmXVVcv0cOV0udhc4yDI23yTruQ+OxGGznv15/vz5zJ8/H41Gw5o1a9y6r4YmqJxzrpgA4oD8izwXD+T93QYHDRrUwFBqBWaUcaB4B8nJyY3aTn0F+/owrmcCgTp5l5Jv3Lix0cfWG13suAwCfj5WxO4MiRdwSwbNU0Usfrs7U9/LICi8eZJmenr6RT9DJQoF3dtFc0VcULPEIidy/ByZbQ5W7M9DZTeRHN58+31+QRrLP3yNA1vXY7OY8fPzY/To0bz55pvExMTUYwt/lpZW9zpoDT0dWguMpfaeUzcg12AwGKm9J5UOBOn1+mS9Xq8GUs+83msEaNWM6SH/5CQ0zDWXR8jiS7hjLxN9byzn85djcdiljqb2bH3t4XwOZIvCCamVVltZ+FsmWaWmZt/3uWuZ+Wg0mM1mgoKCGpyc/kqDEpTBYNgK7NLr9VuBOcCDer3+Tr1eP/rMS+4HFgKbgcUGg+GYe8OWjs5Hxeju8WISrhdTKBQM6xhNmyjp5yQNu60UnZ+TH/5PHkUcfxRO7M6U5xLhLUFmiYlFOzIpN0k3FeGPtczemj2X++67r8kKJRo8D8pgMDx1wUP7znluE9C3UZFdouLCfJa88TiTXprfZKuS+qgU/KNrHBGifZHXUyoVjOgUw3d7csguq5EwDrj9yXxmTU6iVQcz3QYbJYvlXD8fLcJsc9CvjfSFJS3JrowythwvlnyhwYnT3wPgMnJ46L67m2w/XnPH86N33yT3xEHWfjm3SbavUipITYkjLkReawsJTUetql0mReoF/fwCndz5Qh7fzosk77RG0ljOtf1UKesPF+D0gIXvPJ3V7mTVgTw2HSuSPDk1J49PUL6+vigUCpZ8+Qm4XGxdsZCpw9szLTXFrfu5tmM0yWLJjBZHq1ZxU7d4giQe0o1vY2HUpCI+mxlHTbV8PrYHcipYeSBPLNXRhMpNVhbvzOJovjyunpuTfN7pDXTq1Cluu+02tLraKxsfrY7uQ0by3IK6q0Iaon9bedw0F6QRoFUzupv0E7F7DTPSrruJr1+PwSmjfHCisOrMcg4yqOTwMkfzjXy1PZNio0XqUCTh8QkqNjaWoKAgrBYLKh8NdqsFnV+A2+5DpSQE07u1fNYQEqQR5q9hVNc41EppWxD9475CqitUpC2S13syt9zMot+yKKlqmV+k7mZzOFl3uIBVB/I8agVcd/P4BAVQUFDAmNsnon/qHfqljsdYVuyW7V4W6c/g9vLo0yZILz7El+s7xaCQMEepfWDC87n88kMIR3b6SRdIHSpqbCzakUV6cbXUoXi0IqOFRb9lij6Ijanik5Nvv/2WXRllLFi7g14D3NP2LypIyw2dYlFKfMYsyMvl0YEMMNvYdMw9J0ENERzu4F/P5PHZS7FMmZ1JeKx8htasdiff782lX9twWa1e7AlcLhc7M8r49WRJs/XTkzuvuIJyt0CdmlFd4tCoxeER/qxHqzBSEoIljeGyzjVcO76UT1+Mw2qW10mU0+Viy/FiftiXi8UufdsoT1BhsvHNzmy2HC8Wyekc4hv4Ahq1klFd4kSXCOEvDW4fRWuJqzoH3FRO3GVWDLOk73xelxOFVSzcnkmxuC91UU6ni10ZpXyxLZ2ccunm28mVSFDnUCjg+k4xRAXppA5FkDmlUsGIzrFESDhHSqGAsVMKKMiSR+fzupSZbCzcnsmezDJccsyiEio0mlm0I4tNx4q9Zv0mdxMJ6hxXt42gTaT07W0Ez6BRK/lH1zj8tdKVn2u0Lia+kMuGxWEc3yvPSeR2p4uNR4tYtjeHaot87pdJxWxzsPFoIQu3Z1FQaZY6HFkTCeqMjrFB9BQ3dYVLFKTzYWQXacvPw2Ls3P5UHl++GktpvnzrntKLTXy5LaPFri3lcrk4mFPB51vT2ZNZ3qI6QjSUSFBAbLCOazuKcnKhYWKDfRl2ZbSkMbTrVsMQvTyLJs5lsjpYfSCfZXtyqDRL1+y0uWWWmFj4WxbrDhdgsorCkfpq8QkqUKcmtUscalWLPxRCI3SICeIqiSd0X3NzOTHJVhbLtGjiXKeLq/ni1wx2ZZR6ddVaYaWZb3dns3R3thjOa4AW/a2sPtMANkAr32ERwXP0bRNOWwmX6FAoYNwjBRRla/jpm1DJ4qgvq93JpmPFLPg1nROF3tVnLr/CzPd7c/j6t0wySpp/zSZv0aK/mQd3iCImWFTsCe6hUCi47soYymuyJOudptG6mDg9l9kPJxGbbKFjb/l/OZabbPywL4/40HL6tQknIVReHTIuRVapiR3ppSIpuUmLvYLqkhhMp3hpJ1sK3uePeXS+Gukq+0Kj7Ex4PpeFb8RQkOk58/lyymr4Zmc2hp1ZZHrQF7zN4eRgTgVfbMtgya5skZzcqEUmqPgQXwa2E0URQtMI9vUhNSUWlYSVfa2vNHPjv4v5ZHo8JqNnfcxzympYujubr7fX9qOzyXQpj4JKMz8dLeTjLadZd7igxXYcb0otbogvQKvmRom/PATvlxDqx8B2kWw4UihZDFddX0neaS0LXo7lnpdzUEm7WsglK6g0s+6wmU3Hi7giNogOMUGSD8mXVls5WVTFkbxKiqusksbSErSoBKVSKhiREou/KIoQmkGXxBAKjRZJu1KPnFTER8/H8/37kdz8YJFkcTSGxeZkT2Y5ezLLCfL14fKoANpEBRATpGvyE02bw0leuZlDxQ4ytqZTWi2S0h+O7vIjqySGG4c13T5a1Df1gMsjiBdLtgvNaEiHKEqrLeSWS1NirFLBv57JY/aURH5ZHkz/UZ69hENljY1dGWXsyihDo1YSG6wjMcyPqEAtEQHaRp18ulwuykw2iowWiowWcstryK8043C6SC9zkBwskhNAdaWS5R9EcmK/H49MKWnSfbWYBNU+JpBuSfIvvRW8i0qp4MaUOBZuz6RKojY/vgFO7n4plzmPJhIRZ6N9T++4iW+1O8koMZ1XlOCrURHmp8Ffq8ZfqyJAq0atUqJSKM5ebdkcTmwOJ1aHk2qLgyqLDaPZTmWNTfTE+wsuF+zdGMiyDyLpOtDItA/T6exb2qT7bBEJKjxAw7UdpZ3pL7RcAVo1qV1i+WZntmSTUiPibEx4Lo/PZsby4FvZRCd559VAjdVBjlV0BXe3yhIVS96NpijXh4nTc0nu2DwjAp5V3tMAGrWS1BSxtpMgrdhgX8lWZ64sKeS9x+4gIj6T1HuK+ej5OIxlHlYxIUjC5YKd6wN58/5WxCRbeGxuZrMlJ1pCghraMYowf43UYQgCnROC6SzB3Lu1X83j9MGdrP1yLr2HV9JtkJGPp8dhtYhKVuHiKktVfDI9jp++CeOe/+QwYmIJak3zjgB49RBfSkIwHWKCpA5DEM4a3CGK4ioLeRVNfxY6LTUFu/V/c3O2rljI1hULUflo6XJ1IV+/FsO/nstD6fWnqcKl2rMxgO/mRdF3RAUTns9FLdF8b699a0YH6RjYLlLqMAThPCqlgtQuzbOG1HOfr6f74FR8tLVzh3y0OroPGcnzX6Rx62MFVJWrWPFRRJPHIXiO6kolC16OZc0X4dz9Ug433FkiWXLCWxOU1kfJjZ1jRYdyQZYCtGpGdI5FqWjaIbag8Ch0fgHYrRbUGi12qwWdXwBBYZGoNS4mzsjl8LYANn8vz9V4heb1+29+vHFvMsHhdqbOyySpvfSdMbxyiG/4FdEE+3lODzKh5UkI9ePqyyPYdKxpJ88ay0volzqePiNuYduqxVSW/m9//kFO7nk5mzlTkwgKtdPlmqomjUWQJ0uNguUfRnJkhz+3P5nH5V3lUwXpdQmqa2IIbaMCpQ5DEP5Wj1ahFFSaOZrfdEtNTJz+3tl/j3lo+p+eD4+1c89/cvjgqXj8gx2oRf/kFiXjdx1fvRZD8hU1PP5BBr7+8up76FVjYKE6BdeI+06CB7m2YzQRAdJWmca3sfDPZ/JY8J9YirOlW89KaD4OB/y4IJyPp8dx413F3DatQHbJiYZeQen1eh/gM6AV4AAmGgyGUxe8xgb8cs5DQw0GQ5OtdaxVK+kVoxZNYAWP8sc8vYU7MrHYpPuCuLxbDaMfLOTbd7vTuk0e4bEtZzn2lqYox4evXotB5+fksfkZBIfLdwn6hg7x3QaUGwyG2/V6/XDgVeCWC15TYTAYBrkhxnq5IjaI4uMiOQmeJ9Rfw3VXxvDDvlxJl2rvNqiKrPRy3n+qLZPfzpL1F5dw6Vwu2P5jECs/iWDYbaVc/Y9y2U8xaGh4Q4Hvzvx7PdDfjTE1iFJcOQkerE1kAL2Sw6QOg5RB2Vx1fSUfPJVAdaXMv72EequuVPLZzFg2fx/KA69nc81o+ScnAIWrAadser1+LfCEwWDYd+bnLKCNwWCwnvOaKuB7IBlYajAY3r7Y9tLS0lwqNyxWU1VVRUCAGEM/lzgmdZPjcXG5XGzNtVNoku4yymq14uOj4Zel7cg5FsroqTvR6MSVlNVqRaPxzI40mYfDWPdZJ9r1yqfvTcdR+7jv/dU/0kpUSOM/Rw6Hg6FDh/7pKuNvh/j0ev3dwN0XPHzVBT8rgAv/6seBL888vkmv128yGAw7L7afQYMaPxq4ceNGt2zHm4hjUje5Hpc+Vgdfbc/AaJam83l6ejrJyckkP2blm9ku1v5fX+75Tw5a35bd5fuP4+JJbFYFKz+OYP/mAP75VAHtulvPlA24jz85bvkcpaWl1fn43yYog8HwEfDRuY/p9frPgBhg35mCCYXBYLBd8Hvvn/P6NKAzcNEEJQhC7XIRqSlxfLMzC7tEnc8BFAoY+3Ahi9+O5pPp8fx7Zg4aXctOUp4k97SGr/4bS2S8lcfez8A/SH4VevXR0CKJtcA4YA0wEvjp3Cf1en17YDpwO6ACrgaWuCdkQfBuMcE6BrWPYv3vBZLGoVTCLY8W8PUbMXwyI45/z8zFp5mbhQqXxumEzd+FsH5RGKn/Lqb3dZU0ccOSJtXQ22SLAZVer98CPAg8TW1iekqv1/c1GAxHgUzgtzOl5isNBsNv7g1dELxX54RgroyTvtGxUgXjn8jHL9DBpy/GYbN68LedlysvVvPh0/Hs2xzIlNlZXHW9ZycnGnoFdWY+08Q6Hv/vOf9+qrHBCUJLNqRDFEVVFgorpe2JplLB7U/m8/XrMXz0fBx3zcht8fek5OaP7uP9R5Vz7fhSqssLee+xqfzr2VkEhXlu8wIPKDQUhJZJrVKS2jkOnY/0iwuq1LVJKiTCzv89G4+5Wnx1yEFNlZIvX41hzYII7n4ph+vuKEWlOn8NME8m3mWCIGPBfj5c3ylGFkM1ShXc8lgBUUlW3n86HpNRfH1I6dhuP964txW+AU6mzssgqb2FaakpTB3enq0rFtZOW1ixkKnD2zMtNUXqcBtEvMMEQeZaR/hzVetwqcOAM4UT46YUktzRzHuPJVJRIv3VXUtjNSv4dm4ki96MRv9oAWMeKjxbYXmxNcCeW1B3GbfciQQlCB6gz2VhtI7wlzoMOFOC/o/7iug+uJI5jyRRmC2Wtmku6Yd1vHV/K0xGFY9/kEGHnqbznv+rNcA8kdcttyEI3kihUHB9pxi+3p5JRY30jVwVCrh2fBkBIQ7mPp7I3TNzSGwn/QJ33spqUfDj5+Hs2hDEzQ8W0mXAxdfu+qs1wDyNSFCC4CF0PipSu8Ri2JGFzSGPKro+N1QSEOzg/56NRz+1gE59q6UOyeukH9ax6M0Y4tpYeOL9DAJC/rr11N+tAeZJRIISBA8SFahjSIdo1hzKlzqUszr1qyYoPIdPX4yjOKecgWPKZFHU4eksNQpWfRrB3k2B3PxAYYtc8VjcgxIED3NFXBBdE0OkDuM8Se0tPPxOFjvXB/LNO1E4pGkl6DWO7PTj9UnJ1FQpmfZheotMTogEJQie6Zp2kcSH+EodxnlCo+xMfjsLY5maeU8kUikq/C6ZsUzFV6/FsGR2NOOmFHDbtAKP7aPnDiJBCYIHUikVjEiJJUArj1H6ypJC3nvsDqzmQibOyKV9j2rentyKUwd1UofmEZxO+HVlMG9MakVQmJ0nPkz/U4VeSyQSlCB4qACtmhtTYlHJYLHOczsXKJUw/I5Sbpmaz2cz4/h5aYikKwXLXdYxLXMeTWTHuiDuey2bkfcUi1ZSZ8jj9EsQhAaJC/FlYLtINhwplGT/01JTsFv/V16+dfchGCUAABGASURBVMVCtq5YiFqj5fUV+5kyO5MvXonl6G5/xj+eT2CoWPzwD1XlKlZ9Gs6h7QHcMKG287gnrHLbnMThEAQP1yUxRLLO53/XuSA81s5Ds7JIaGvmzftbcXi7PCYbS8luVfDz0hBeu6cVGp2Lpz5Kp88NIjnVRVxBCYIXGNIhipJqK/kV5mbdb306F6jUMGJiCe17mPj6jRgO/GJi5D1F+AW2rJv/Lhfs2xTAio8jiGll5cE3s4lpZZU6LFkTOVsQvIBapSQ1JRZ/bfNXzv3RuWDKbAP9UsdjLCuu83VtUmp44v0M1D4uXp+UzL7NAc0eqxRcLji6y4/ZDyeStjiMWx4t4O6XckVyqgdxBSUIXiJQ50NqShxLdmXjaMbl4i+lc4HO38mYhwrpNrgSw6xodq4PYtSkIiLjpW/f1BRO7PPlx8/DqSpXc90/S+gy0CiG8i6BSFCC4EXiQnwZLIPl4v/OZZ3MPDY/k5+XhvDulCR6Da9g2G2l+AZ4/rCfywW/7/Bnw+JQKorVDL+jhO5DjKjEtLBLJnK5IHiZzgnBpCQESx3G3/LRuLh2fBlPfJhOTbWK//47mY1LQrHUSF823xB2G+xKC+TN+1ux8pMI+qVW8NQn6fQadn5y+mPOmCc3cW0u4gpKELzQoPZRlFZbyS6rkTqUvxUU5uCWRwvIPa1h3VfhbDC0ZuDNZfQfVY7OT/7zgcqL1Py6MphtPwYTk2Ql9a4iOvQyXbQf4blzxsY+PKO5w/UoIkEJghdSKRWkpsSx8Dd5LM9RH3GtrUx4Lo/8dA3rFobxn39eRo+hlfRLLSc6SV5/g82q4NCv/uxcH8Tpw770GFLJA69nE5108cKHv5szJvyZSFCC4KV8NSpGdonDsDMLq91z7u3EJFv559P5lBao2bYqmHlPJBLTykrPayvp1K9KsvtUNquC43v8OPBLAPt/CSCxrZmewyr55zN59er88Nzn61n+4Wsc2Loem8WMj1ZH5/7DGDXpyWaJ3xOJBCUIXiwyUMt1V8awYn+ux7UbCou2M2JiCcPvKOHALwHs2RjEd/MiuaxzDSlXV9G2q4mw6KZrm+5yQVGODyf3+bF7UxdyjkYQ18ZCp75VPPHPEkIiL23f3rbabXMQCUoQvFzbqAD6t41gy/G65yfJndoHug2qotugKszVSg5t8+fQtgBWfBSB1s9J2y41JLYzE9vaQmyyFZ3/pV9huVxQUaIm95SW3JNack5qOXXQF5XKRdsuJi7r8v/t3X1sXfV9x/H3fX6+tq8dPyaO4zxAUmiewRAKKYFqooFRyn7TBlrXFmmb0DbG1gq13VA1qbBpDK1j01qhilXd0H6kg/DUgmbVtJRSnil0VWgBQ7Ahz4njBD9c2/vDJ8yE62vnPp1z7/28pCjX5xyf+/VX55zv/f3u7/zOXj73leMkG4prvdXS024rQQVKpA5s7clw+MQE/zs84nYoRYkmptm84zibdxxnZgbeGwzzm1/EeXtPlJ//sIF9b4WJpaZoaM6Sbp4inckSiU0TDM0QCM7g88H4mJ+JMR/jJ/0cOxTk6P4QRw4EiUSn6ewdp3PlOB/rG2XnFw+Qac/i88Hg4DDJhp6i46+lp91WggqUSJ24bG0bx96fZKgKRvYths8HHSsm6Fjx/wMTpqdmR9WNHA5y7NDs/xNjPqYmfWQnZ4fVReNTNDTPEIlNk27O0tSapXHJpGYQ9yAVKJE6EfD7uLLKRvadKX8AMu1ZMu2FfTc1cmg/3/3GzfzBV+/Ud0MeoBt1RepILBzg6o1dREOa1iCXufcoifvUghKpM5lEmCvXd/DfLwxVdM4+L9M9St6kFpRIHVraFOfydW1uh+EZCz3XStyhAiVSp9Z2pLlgZbPbYXiC7lHypoK7+IwxlwD3AV+w1j6cY/11wE3ANPAta+13io5WREqqr7eZ0bEsrwwdczsU1+keJe8pqEAZY1YCNwNPzrM+AfwNcB4wATxrjHnAWnu46IhFpKQuPbuVExNZ3jhwwu1QXKV7lLyn0C6+d4FrgPnu+jsfeNZae8xa+z7wU2BbEXGKSJn4/T6uOLeDzsao26GIfEhBLShr7UlmW0rzbdIOzG0f7wc68u1zYGCgkFA+ZHR0tCT7qSXKSW7Ky0c1Tc0QnplgcHDQ7VA8Z2JCecmla8lEWc+jBQuUMeYG4IbTFt9qrX0sz6+d/iQUH5B3POv27dsXCmVBAwMDJdlPLVFOclNecpua+RHD0R5GavRG3kINDg7S01P8VEe1JsFQSc6j/v7coyUXLFDW2ruBu8/w/YaAnXN+7gKePsN9iEiFxYI+rtnYxX3P7+XE+JTb4UidK9eNuj8H7jbGNAJZ5/unm8r0XiJSQk2JMFdv7GLX8+8wPlk9z5GS8gkH/aSjQVLREJGgn3DQTyjgZ2zvUFnft9BRfJ8GvgScDWw2xvyZtfZTxphbgCestT9zXj/mdO193VqrcawiVaI1FeXqDV3c/+JQVT3sUIrj80FLMkJbOsqSVITWVIRMIjzv1FgDw/M8175ECh0k8QjwSI7lt895vQvYVWyAIuKOzsYYV63vZPdLQ0xOaUqkWpWOhehtSbAsE2NpU9xT8zRqLj4RmdeyTJyr1nex+6Uhspq3r2Y0xUOsaUuxqjVJa9q7txeoQIlIXt3Nca5c38lDLw+rSFWxSMjPmtYU6zrTdDbG3A5nUVSgRGRBPS0JFakqtSQVYcOyRs5qTxEKVNf0qypQIrIoPS0JfntDFw++rO+kvM7ng1WtSTZ2N9FVJa2lXFSgRGTRupvjXL2xi90vDWt0nwcF/D7WdqTZsryJpkTY7XCKpgIlImdkaVOcazZ18cCLw4xN6mZeLwj6fZzT1cCWniZS0ZDb4ZSMCpSInLGOhhi/s2UpD7w4xPGxrNvh1K2A38c5XWm29mRqqjCdogIlIgVpSUYwW5dx/wtDHD4x4XY4dcXv83F2R4q+3mYaYrVXmE5RgRKRgqWjIcyWZex+aYh3j425HU5dWNmaZNvKZpqTEbdDKTsVKBEpSiwc4NrNS3nsl/t4bd9xt8OpWV2NMS5a3VI19zCVggqUiBQtGPBzxbntNMZDPPOmHpxdSplEmG2rWljVmnQ7lIpTgRKRkvD5fGxb1UJTPEz/r/bpht4iJSIB+nqbOaezAb+/vJOyepUKlIiU1LrONM3JMA+9PKwRfgUIB/1s7G5ky/IM4WB1zfxQavX914tIWbSlo1x3/nK6M/GS7XPk0H7u+svrGTl8oGT79BK/z8e5XQ187sIeLlzZUvfFCRUoESmXWDjAZzZ2cX5vBl8Jeqge/49/5c1Xn+Px7/1LKcLzlN4lCa7v6+aydW0kI+rYOkWZEJGy8ft9XLiyhe5MnB+++l5BXX5f3vlxshPjH/z81MP38tTD9xIMR/j7h39R4ogrq6sxxrbVLVU9X145qQUlImW3tCnO9X3LOas9dca/+7V//x82fXInocjsc4tCkSibLr2Sr323vwyRVkZLKsJVGzoxW5epOOWhFpSIVEQ0FOCKcztY3ZrkR3v2c2J8cfP4pZtbicaTZCfGCYYjZCfGicaTpDNLyh5zqWUSYfp6m1nTlsRXin7PGqcCJSIVtbotxbJMnJ/8+iCvDh1b1O8cP3qIC3f+Hn1X/C5PP/pfVTdQIpMIs7Unw9ntqbodMl4IFSgRqbhoKMDl69pY25HiidcOsH9kPO/2n7/1rg9ef/ZPb61AhKXRkopw/ooMq1vVYiqECpSIuGZpU5zfP6+bXw6P8NTrBxfd7ed13Zk4m5Y30dMcV2EqggqUiLjK55t9ltGathQvv3OUF946wsmJ6itUQb+P1W0pNi1vpDUVdTucmqACJSKeEA762dqTYf3SRl4ZOsrzbx2pihZVYzzEx5c2sK6jgVg44HY4NUUFSkQ8JRz0s3l5hg3LmvjN/lFefucoQ0fedzusD4mGAqxpS3J2R5rOhqi68cpEBUpEPCng93FWe4qz2lMcHB1nz3vHeW3fcY6enHQlnng4wIqWBCtbk/Q0JwhoNF7ZqUCJiOe1JCO0rIqwbVUL+0bGeH3/KHuPnOS9Y+NMz5Rn1vRQwEd7Q4ylTTG6s+9wzcW9ailVmAqUiFSVtnSUtvTsIISJ7DRDR99n/8gYB0cnODg6ztGTk2dctGLhAI2xEEtSEVqSEVrTEVpT0Q9aSQNv+1WcXKACJSJVKxz0s6IlwYqWxAfLpqdnODk5xYnxLKPjWbJTM2Snp5menl0fDPgIB/2EA34SkSCpaJBQQLO+eZEKlIjUFL/fRzISJBkJ0uZ2MFKUgguUMeYS4D7gC9bah3OsnwR+OmfRDmut98eMioiIJxRUoIwxK4GbgSfzbHbMWru98NBERKSeFdrx+i5wDTBS4nhEREQA8M0UMUTTGHMPsGueLr5RYDfQA3zfWvuP8+2nv79/JhAo/g7s0dFRkslk0fupJcpJbspLbspLbspLbqXKy9TUFDt27PjIMMkFu/iMMTcAN5y2+FZr7WML/OpfAd8DZoAfG2N+bK19br6Nt28vvjdwYGCgJPupJcpJbspLbspLbspLbqXKS39/7odPLligrLV3A3ef6Rtaa//t1GtjTD9wLjBvgRIREZmrLMPMjTFnAbcC1wEB4CJgVzneS0REalNBgySMMZ82xgwAvwXcZox53Fl+izHmAmvtHuBt4BlnqPkj1tpnSh69iIjUrIJaUNbaR4BHciy/fc7rW4oNTkRE6pfm9xAREU8qaph5qfT397sfhIiIuCbXMHNPFCgREZHTqYtPREQ8SQVKREQ8SQVKREQ8SQVKREQ8SQVKREQ8qSqfqGuMuRPocyai/XNr7bNz1l0GfAOYAh611v6tu9FWzgJ5+SRwm5OXPcAN1tppdyOujHx5mbPNbcAF9fQMswWOl2XAvUAYeMFa+8fuRls5C+TlRuB65zx6zlp7k7vRVo4x5hznCRV3WmvvOm1dWa67VdeCcp7ku9paewHwReCbp23yTeCzwDbgU8aYdS6FWlGLyMu3gWuttduAlDNNVc1bRF5wjpGL3YnQHYvIyx3AHdba84ApY0y3S6FWVL68GGPSwJeAT1hrLwLWGWP63I24MowxCeCfgdzTjpfpult1BQrYATzA7HRKvwKanAMHY0wvcNhau9dpHTzqbF8P5s2LY7O19h3n9QGg2Z0wK26hvOBcjL/qTniuyXce+YFPAA8662+01r7tdsAVku94mXD+JY0xQSAOHHY33IoZB64Ahk9fUc7rbjUWqHbnAnvKAWdZrnX7gY4Kx+eWfHnBWjvC7MHUAVzuHET1IG9ejDF/CDwBDLoTnmvy5WUJcBy40xjzpDHmNmPMR+7yr1Hz5sVaOwZ8HXjDOV6etta+5l6olWOtzVpr359nddmuu9VYoE4/UXxOX/FC62rdgn+7MaYVeAi40Vp7qLLhuWbevBhjMsDnnRZUvVnoPOoC/gm4BNjofHquB/mOlzTwFWAN0Av0GWPWuxOmp5TtuluNBWpo7idgoBN4b551XcC7FY7PLfnycurk+gHw19bax90J0RX58nKp01r4CXA/sMn5grwe5MvLQeAta+3r1top53uHj7kUZ6Xly8ta4A1r7UFr7YRz3Gx2KU4vKdt1txoL1OPAtcxedDcCw9ba48w2QweBtDGmx+kj3ulsXw/mzYvjDmf0zQ/cC9EV+Y6XXdbaddbaPuAzzmi1v3A74ArJl5cs8IYxZrWz7WZn5Gc9yHceDQJrjTExp8tzC/Brd8N1Xzmvu1U5Wawx5nZn1NU0cKPTBXHMWnu/MeZi4O+cTb9vrf0Hl8OtmPnyAjwGHAF+Nmfz/7TWftvFcCsm3/EyZ5se4J46G2ae7zxaBdzjfIh9BfiTOrotIV9e/sjpFs4CT1lrv+x2vJVgjNnsfMjtASadVtODwJvlvO5WZYESEZHaV41dfCIiUgdUoERExJNUoERExJNUoERExJNUoERExJNUoERExJNUoERExJP+D1bC4aQjqAF4AAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "bento_obj_id": "139942981246032"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "from matplotlib import pyplot as plt\n",
        "\n",
        "%matplotlib inline\n",
        "\n",
        "# Initialize plot\n",
        "f, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
        "# test model on 101 regular spaced points on the interval [0, 1]\n",
        "test_X = torch.linspace(0, 1, 101, dtype=dtype, device=device)\n",
        "# no need for gradients\n",
        "with torch.no_grad():\n",
        "    # compute posterior\n",
        "    posterior = model.posterior(test_X)\n",
        "    # Get upper and lower confidence bounds (2 standard deviations from the mean)\n",
        "    lower, upper = posterior.mvn.confidence_region()\n",
        "    # Plot training points as black stars\n",
        "    ax.plot(train_X.cpu().numpy(), train_Y.cpu().numpy(), \"k*\")\n",
        "    # Plot posterior means as blue line\n",
        "    ax.plot(test_X.cpu().numpy(), posterior.mean.cpu().numpy(), \"b\")\n",
        "    # Shade between the lower and upper confidence bounds\n",
        "    ax.fill_between(\n",
        "        test_X.cpu().numpy(), lower.cpu().numpy(), upper.cpu().numpy(), alpha=0.5\n",
        "    )\n",
        "ax.legend([\"Observed Data\", \"Mean\", \"Confidence\"])\n",
        "plt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Interfacing with Ax\n",
        "\n",
        "It is simple to package up a custom optimizer loop like the one above and use it within Ax. As described in the [Using BoTorch with Ax tutorial](./custom_botorch_model_in_ax), this requires defining a custom `model_constructor` callable that can then be passed to the `get_botorch` factory function."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "def _get_and_fit_model(Xs, Ys, **kwargs):\n",
        "\n",
        "    train_X, train_Y = Xs[0], Ys[0]\n",
        "    model = SingleTaskGP(train_X=train_X, train_Y=train_Y)\n",
        "    mll = ExactMarginalLogLikelihood(model.likelihood, model).to(train_X)\n",
        "    model.train()\n",
        "\n",
        "    optimizer = SGD([{\"params\": model.parameters()}], lr=kwargs.get(\"lr\"))\n",
        "    for epoch in range(kwargs.get(\"epochs\")):\n",
        "        optimizer.zero_grad()\n",
        "        output = model(train_X)\n",
        "        loss = -mll(output, model.train_targets)\n",
        "        loss.backward()\n",
        "        optimizer.step()\n",
        "\n",
        "    return model"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "python3",
      "language": "python",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}