{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GP Regression with Grid Structured Training Data\n",
    "\n",
    "In this notebook, we demonstrate how to perform GP regression when your training data lies on a regularly spaced grid. For this example, we'll be modeling a 2D function where the training data is on an evenly spaced grid on (0,1)x(0, 2) with 100 grid points in each dimension. \n",
    "\n",
    "In other words, we have 10000 training examples. However, the grid structure of the training data will allow us to perform inference very quickly anyways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gpytorch\n",
    "import torch\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make the grid and training data\n",
    "\n",
    "In the next cell, we create the grid, along with the 10000 training examples and labels. After running this cell, we create three important tensors:\n",
    "\n",
    "- `grid` is a tensor that is `grid_size x 2` and contains the 1D grid for each dimension.\n",
    "- `train_x` is a tensor containing the full 10000 training examples.\n",
    "- `train_y` are the labels. For this, we're just using a simple sine function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_bounds = [(0, 1), (0, 2)]\n",
    "grid_size = 25\n",
    "grid = torch.zeros(grid_size, len(grid_bounds))\n",
    "for i in range(len(grid_bounds)):\n",
    "    grid_diff = float(grid_bounds[i][1] - grid_bounds[i][0]) / (grid_size - 2)\n",
    "    grid[:, i] = torch.linspace(grid_bounds[i][0] - grid_diff, grid_bounds[i][1] + grid_diff, grid_size)\n",
    "\n",
    "train_x = gpytorch.utils.grid.create_data_from_grid(grid)\n",
    "train_y = torch.sin((train_x[:, 0] + train_x[:, 1]) * (2 * math.pi)) + torch.randn_like(train_x[:, 0]).mul(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the Grid GP Model\n",
    "\n",
    "In the next cell we create our GP model. Like other scalable GP methods, we'll use a scalable kernel that wraps a base kernel. In this case, we create a `GridKernel` that wraps an `RBFKernel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GridGPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, grid, train_x, train_y, likelihood):\n",
    "        super(GridGPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        num_dims = train_x.size(-1)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.GridKernel(gpytorch.kernels.RBFKernel(), grid=grid)\n",
    "    \n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GridGPRegressionModel(grid, train_x, train_y, likelihood)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 0.978   lengthscale: 0.693   noise: 0.693\n",
      "Iter 2/50 - Loss: 0.933   lengthscale: 0.644   noise: 0.644\n",
      "Iter 3/50 - Loss: 0.880   lengthscale: 0.598   noise: 0.598\n",
      "Iter 4/50 - Loss: 0.816   lengthscale: 0.554   noise: 0.554\n",
      "Iter 5/50 - Loss: 0.738   lengthscale: 0.513   noise: 0.513\n",
      "Iter 6/50 - Loss: 0.655   lengthscale: 0.474   noise: 0.473\n",
      "Iter 7/50 - Loss: 0.575   lengthscale: 0.438   noise: 0.436\n",
      "Iter 8/50 - Loss: 0.507   lengthscale: 0.403   noise: 0.401\n",
      "Iter 9/50 - Loss: 0.450   lengthscale: 0.372   noise: 0.368\n",
      "Iter 10/50 - Loss: 0.398   lengthscale: 0.345   noise: 0.338\n",
      "Iter 11/50 - Loss: 0.351   lengthscale: 0.321   noise: 0.309\n",
      "Iter 12/50 - Loss: 0.306   lengthscale: 0.301   noise: 0.282\n",
      "Iter 13/50 - Loss: 0.259   lengthscale: 0.283   noise: 0.258\n",
      "Iter 14/50 - Loss: 0.214   lengthscale: 0.269   noise: 0.235\n",
      "Iter 15/50 - Loss: 0.167   lengthscale: 0.256   noise: 0.213\n",
      "Iter 16/50 - Loss: 0.120   lengthscale: 0.245   noise: 0.194\n",
      "Iter 17/50 - Loss: 0.075   lengthscale: 0.236   noise: 0.176\n",
      "Iter 18/50 - Loss: 0.027   lengthscale: 0.228   noise: 0.160\n",
      "Iter 19/50 - Loss: -0.021   lengthscale: 0.222   noise: 0.145\n",
      "Iter 20/50 - Loss: -0.067   lengthscale: 0.217   noise: 0.131\n",
      "Iter 21/50 - Loss: -0.118   lengthscale: 0.213   noise: 0.118\n",
      "Iter 22/50 - Loss: -0.167   lengthscale: 0.210   noise: 0.107\n",
      "Iter 23/50 - Loss: -0.213   lengthscale: 0.208   noise: 0.097\n",
      "Iter 24/50 - Loss: -0.267   lengthscale: 0.207   noise: 0.087\n",
      "Iter 25/50 - Loss: -0.311   lengthscale: 0.207   noise: 0.079\n",
      "Iter 26/50 - Loss: -0.363   lengthscale: 0.208   noise: 0.071\n",
      "Iter 27/50 - Loss: -0.414   lengthscale: 0.209   noise: 0.064\n",
      "Iter 28/50 - Loss: -0.460   lengthscale: 0.212   noise: 0.057\n",
      "Iter 29/50 - Loss: -0.514   lengthscale: 0.216   noise: 0.052\n",
      "Iter 30/50 - Loss: -0.563   lengthscale: 0.221   noise: 0.047\n",
      "Iter 31/50 - Loss: -0.616   lengthscale: 0.227   noise: 0.042\n",
      "Iter 32/50 - Loss: -0.668   lengthscale: 0.234   noise: 0.038\n",
      "Iter 33/50 - Loss: -0.717   lengthscale: 0.243   noise: 0.034\n",
      "Iter 34/50 - Loss: -0.773   lengthscale: 0.253   noise: 0.031\n",
      "Iter 35/50 - Loss: -0.821   lengthscale: 0.264   noise: 0.027\n",
      "Iter 36/50 - Loss: -0.868   lengthscale: 0.276   noise: 0.025\n",
      "Iter 37/50 - Loss: -0.927   lengthscale: 0.290   noise: 0.022\n",
      "Iter 38/50 - Loss: -0.981   lengthscale: 0.305   noise: 0.020\n",
      "Iter 39/50 - Loss: -1.026   lengthscale: 0.320   noise: 0.018\n",
      "Iter 40/50 - Loss: -1.084   lengthscale: 0.338   noise: 0.016\n",
      "Iter 41/50 - Loss: -1.133   lengthscale: 0.355   noise: 0.014\n",
      "Iter 42/50 - Loss: -1.183   lengthscale: 0.371   noise: 0.013\n",
      "Iter 43/50 - Loss: -1.231   lengthscale: 0.385   noise: 0.012\n",
      "Iter 44/50 - Loss: -1.277   lengthscale: 0.397   noise: 0.011\n",
      "Iter 45/50 - Loss: -1.323   lengthscale: 0.402   noise: 0.009\n",
      "Iter 46/50 - Loss: -1.374   lengthscale: 0.403   noise: 0.008\n",
      "Iter 47/50 - Loss: -1.426   lengthscale: 0.401   noise: 0.008\n",
      "Iter 48/50 - Loss: -1.481   lengthscale: 0.395   noise: 0.007\n",
      "Iter 49/50 - Loss: -1.529   lengthscale: 0.388   noise: 0.006\n",
      "Iter 50/50 - Loss: -1.576   lengthscale: 0.379   noise: 0.006\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iter = 2 if smoke_test else 50\n",
    "\n",
    "\n",
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "for i in range(training_iter):\n",
    "    # Zero gradients from previous iteration\n",
    "    optimizer.zero_grad()\n",
    "    # Output from model\n",
    "    output = model(train_x)\n",
    "    # Calc loss and backprop gradients\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (\n",
    "        i + 1, training_iter, loss.item(),\n",
    "        model.covar_module.base_kernel.lengthscale.item(),\n",
    "        model.likelihood.noise.item()\n",
    "    ))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we create a set of 400 test examples and make predictions. Note that unlike other scalable GP methods, testing is more complicated. Because our test data can be different from the training data, in general we may not be able to avoid creating a `num_train x num_test` (e.g., `10000 x 400`) kernel matrix between the training and test data.\n",
    "\n",
    "For this reason, if you have large numbers of test points, memory may become a concern. The time complexity should still be reasonable, however, because we will still exploit structure in the train-train covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "likelihood.eval()\n",
    "n = 20\n",
    "test_x = torch.zeros(int(pow(n, 2)), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        test_x[i * n + j][0] = float(i) / (n-1)\n",
    "        test_x[i * n + j][1] = float(j) / (n-1)\n",
    "\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    observed_pred = likelihood(model(test_x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAADNCAYAAACLr42mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHehJREFUeJztnUmPHEd2x/9R+9JL9cImKWooTXOsGcxiw82cgwEDNuDm2RdJ/ASiToZPFscnwwfDom4+ildfrOUbiDrYPhjw5DRgQLNJYGs0EqXm1l29d20ZPnRUZ0Z25ousqKqurOT7AYmqyMjMiorMl7G8iH8IKSUYhskWuUkngGGY0cOGzTAZhA2bYTIIGzbDZBA2bIbJIFNv2EKINSHEQyHEPSHE6+pz3eI694QQ76jvq0KIj0aQtrNrhvavCyF2hBB3AvveEUJ8IoRoxBz/ybDpibluI5CH7wshVgPxZ/mg8jY2T4JpDJ03srRH3KO1UVw3i0y9YUspNwBsAPhASvmxlPIuABuj/CBwzU0Ab8UdKIR4fdBrBpFSPgDwIYDtwO4NAG9IKZsxx5/bPwzqBdKQUjYDefi++u/93z3LBynlx9T1gmkMnTfKtIfvkTOi62aOqTfsGLaDJc+gqHMjz1cGccv22gE+AnA7EG5EGfUYuWMyViofxnGeBQ8GeMm+UBQmnYBRowyvKaXcVDf9bQDv4/RBuw/gDk5Lp1Up5X1VtdsAEK7W3YMy4MAxqwA2AThCiNellB+r30t6zTOklA/iqrYq3YvquPuhuHUAd6WUt4QQ99Qxd8PpAOD2r9H/vdDP3IhLW4izfAik4VcA/h1ACYH/TZ2n0r0G4IGqIUA1RTaD58fsi8xPdY/vAiBfUC8iWSqx19XD8yaAvwHOqo+r6vM+gH8EsKEe8hvqIeqHzx58Vc1rAv6Dpo5pqM/tQGmX6JoxfKjaruv9Y1W7sf9Qvx0+IWSg7we+a+mAqg2o8CYsCOZDH5XWNwBcDf0edd6qOu7jfrqUsfbzdVP1MUTtM+XnYsS+F54sGfaGlPKBlPJ+qEq7AQBq3yqAhjKe5wBuQm/nRnETyjCklO9FxNtcs0+/On5WDe+3d5UBJb1OVDr+FcAtVbqe65AzEVPFXcXpy6YR8XsUUf/j5/BfOJsqHLXPlJ+D5NELQ5YMOwm/xGmJsIHTEvxXML/xH0K1FwM91v3SfM3ymgDOSlOtB1+VUP0Srt9eDdN/cQXjwulYl1LelVLeDP9GQqJ+t98pdi/i9wZlM/Abq/3rRexLnJ+Mz9S3sftVVwBvCCHcYGndb9cJIdaklBtSyvdU9e6s/arCwGn77ZYQ4r663poQYlWdc08dA5xWJ/vt9wdJr0l0jN2HXsXcVL+9DtWuVC+Us/8B4JcqvoHTJshqOB0Afh5KM5WH/fbvbfUiuY3TDsi1fj708wS+kd1Q5y6qz+1+GlV88Ly+a6p/PxqqX+Asn/q1IdO+BPnJABA8u+vFRAjxTkzTYmpQL441U+/+i8iLVhVnfO5nwFXERh0DG/YLiqrKNqNGuk0DqrS26u1/EeCqOMNkEC6xGSaDsGEzTAaxcnc5jvM6Tn2pa67rnutZNcX3+cUvfsHtACY1vPvuu8J8FPDuv/ydbO7XqUO+evfdd18dSaIsGdiwHcdZAwDXdR84jrPqOM6a67obSePD/Mcry2fff1Ku49etw7PwK69tkWn525f+z5jeW/XfkvGvFUtkfFHktfDWzo9xZeE3Z+EnvcPwKRobLfPYik92f0LG/9d3PyDjn/1R7//6abmOzwL5WP8jfZvrj8zv19qTDhlf2mmR8blW9+z7K6/N46vPd/UDBG1TvVqRTiCAzix9THs+Ph/+6tWnxuv3ae7X8c9//2Fs/D/925uvJL7YmLCpit+GP/JpE+dHNZniGWbq8SBjtzRgUxVvQB+fuzRgPMNMPR3Zm3QSSCY+pPQnZb+tcr1Y1uJWunQ1tnbwqvH6h908Gf8kT1ff8kKv1DQPvqeFd70T8nyvM0vGA8DS0TUy/jWPdjVfKevtvevFihYuN+g8qPTMpUy5Sj/IhaUyGS86/vkr12oRB9BVca9M/wcA6FXpx7lbG11fcVpK5jhsDLsJf7xwA+dn9pjiNYJt6nD4oEBP3PnZzB8MSQXqhjb2yoBtbABaGztnaGPnErSxn0v6of18t0rGP4to3mpt7Kahjb11sW1sAKlrY1971Xh5/bfgDXaCwnGcvhzWDdd170bEJ+p4NmHzCvsA+gycBypBDSqeYbJET8rYLQ7HcdYBPHBd9z6AVRUOxp91PANo9sM2DFxiu6674ZyyDqAZ6PH+FMBNIj6S+Zf9N3ddAvPCD//Z4iMyLX9SpnvNAaAm6CrkgUeXNNue/mbe7x7jsHNwFv51e4U8/z/3fmRIIfDfW7SYybNv6Kp47Rv9Npbn86jt+vtqhhK5ut0l4wGgcGg4Jmcocet+zcirFLQwAHRr9KNIlbZ9WvN0OdWeS+TNSkTHriren+12H/oU1T63AfSFH/sdz6T9xGHVxlZvnPC+m1Q8w2SJjoVdh+xiDefFLkfW8TzxzjOGmUZ6sC/9VRV7w1SbHQY2bIaxoCOHqtavR3WcYcCOZwoeK84wFvQgYjcKx3Hu9Hu7+51n4+h4ZsNmGAs6Mhe7xaEM+Z7jOA8dx9kJRH0KnHZMB44zdjxTcFWcYSzoWZSJyo21ELF/5B3PbNgMYwFVMqeBiRv2D5efnH2/1gJQ9sM3KvSMm9kcPZwTAI4Mo7qeduks+KJ9RQsXOivonviTd/5nj5559b9PrhtSCDz9+txLXKP6DZ3GsJ+6LCVqj/191ee0L794YPZjC4/275iGfHZm/P/QnSmgvaD7sVvz9PknC2ZDas/T8Z1Z4j/sxEdF0WPDZpjs0YF57PokYcNmGAs6hprgpGHDZhgLbDrPLhI2bIaxoCPTbTrpTh3DpJTecCPPxo5VfcJxnDtquxcTf69/3DCJY5i00pGF2C0NDGzYpjmlijuO4zwEr9TAZJQecrFbGrB5vZjmlALAW67r8ppKTGbJXK94gjmlgF+SG+Vd/hy+7M8MSlgJhBePXiLT0jwxyw49M/gbn3VnyPgvj5e18EK7gZ2A7tmzXToNl5+T+tMAgO/t02ksGzTJylVdDOJKqYRSYF/hEj1ARSzQemUAIAt0m7JbMQ1Q8Uuy5YUy2rP68e1Z+vqdGfME6F6VPkaWRjdAJXOG3YeaUxqYvXLLcZx1NUY2ku3ao9jw5eq3dCLyB3Q8AGlo8xx1aMPcilCj3Ko8Pvv++R59g3faZDQAoLprGFn2hH5ga8/Op/HhN77mWWmPHlkmOmb9Llmkq5jBkWVRnCzo+fTF1pEe36ENu5VgjTmTqolH6JT9zHh1nSyPPIucU6o6zLZVVfw5oqvqDDPVpL3Etu4VJ+aUuvDnkd5QYYbJFMN0nlEihaPyKNn2ipvmlL6pZFQfjlP+hWEmRUfmYzcKZT8fEYeMxKNk03l2YXNKGSateJZtbLWmHWW0I/EopcObzjBTxhjb2Ik9ShTp7tpjmJRiWxU34brue6pWvBQz+CsREy+xKzl/6Zii6GnhbYOP+WnXvC7WVoueff/FAS34v/lcl3b+EebwO/h+36Pv6DRWvzXf6KAoQmS8SShhX4/PH3vaPtE1uIFK5jR25gxCCIb1wY6XfXdWa17gyNPdW61FOo2dhnkRPFGjj8kX7ZbliWIc7q5RepS4xGYYC0ZZYo/Do8SGzTAWeFLEbhTKW+Sozz4j9yhNvCrOMNMIVTJTA3RVNfvj0D5WKWWYNGAqmScNGzbDWECV2PRq5hcDGzbDWND10j1WnA2bYSzwhlht8yKYuGE/b/t+4Fy3ooW/9egF33fa5krPowPaj/10e46Ml4/1rpB2tYrj44AY/1PasVA1TLkEgOo27V8tHBjmU4fE/IWU2r5e2TDlctb8GJgE+0+W6Qf9JOCnbtclTvJ6mrsGP3V+tkPGA0ChSE9P9bwRLnyf8hLbdnYXOQPFcZzXHcdZdxznnWESxzBpxdbddVHY+rFjZ6D0p6SpYXFNaooaw0wrXZmL3dKAbVWcmoFyG8An6vsmgHUAPHWTyRS2s7suCtvUrRJV7QaA7UB4KeIYhplqMlliD6JpZuLS8dWz7/NtXX/MNO623imR8QBQ7dFigit5ugNOVota+HpZ70wrNegbWTLocAFAqUYfUzg2dJ6FxA4vX9b/k8wbhAirCSaBzNDXaM8Z4ut+Gq9Xzo/N8kp0B2Iub14RNCfoa9A2d0hFnqPrpcOA4xjYsBPMQGkC6FtoQx0Ty9Pqd7HhlkcnbyeXoFe8Q6sJPu3RRiOPzz+Evz72hfgqzfH3ihf36Yc6FzF768svfaFHz6AwOpJe8Z6hV1zoafzNoS5m2C0aesXLCXrFc/a94j82Xj10rZR0ksVh89qJnIESmKHyAXxjXw0cyzCZIXNVcdd1N9TyPtvQZ6B8CuCminfUJPGmaYbKF3uXzr53OvP4Q8d/6+6c0CVyc99cYneaFTK+9Iyuhlae62/mymIOM9u5QLxB83vXPI+4cGjyU9Pn90Ka3l4pp+0LanpHkWRR+ZMlgzywYT51t+HfV6/UQ7esl67FWbpmVSyZq+ImugmabklJe4lt28Y+NwOFNc+YF4m0t7HTnTqGSSnDDFAxyA+PZHAXGzbDWNCTudiNgpIfHuXgLjZshrGg5+ViNwpltHHyw7dx6lUC/MFdVkx8EgjDTCNj6jwb2eAuNmyGsYAumUenhmoLGzbDWJBg8U8bBhrcRTFxw9587K8/XcnNYjMwOqi7R/sdC03zUMj6Dl1lqmwb/NBN3cdclR5mHvn7igf02znfGv7t3avQ7bZ2SPO7U8+hFdD5PlkwzZU2Vyvbi/T/6DVoP3NpxvdTF0QJpWpLiy8bRpbJBFXfk2P6eekdjO5xpzvJBrvnjuM0XNdt4nRwl6N2DzW4izvPGMaCMcoPI+ngLoqJl9gMM43YqrGw/DDDpBiTW2vSsGEzjAVj6jwbGTYL3685jiPVwvcPHcd5P+IYUhONYaYdz8vFbmnApsRedF1XAGdD4JoRx9xRnQNvD5M4hkkrKS+wraZtBrvgnZjGPqWJxjBTjxyhlPE4sK43qC75D2OiKU00hpl6pBSxWxoYpvPsVpzW2SCaaD/e8YfDXi+XgZY/yKB4QGdS8cBcISq26WOKOcMAk9Bi6pfLBZRrgYEQtI4DkEBXvluhD+rUTYL/evyluTLac/6+Nr0mAjozCQZU1OkBKOUKfa8KBT8912QJZaFnnOfVyPPbHXNG9lpFMj7fpq6xTcSdZ5SLD4yDYQw7ckpZAk00jbD2VTBcag43agwAys3hFE4Kh+cf6C83fT2x8CocYaRBbwwAOnX6NgRHkUVx0jpv+L9/5ufjieH3WwWzygsq9MiwKlpkfFno+fi52NXCPUG/vI492mgBoNOmFrAFisTIs0E1zzJZFXcc55yxBjTPIjXRGCZTSGJLAcP0zYfnlAaHxb2pesUfDjMsjmHSivRE7JYGbDXPNhFyZbHmGfMikZZOsjh45BnD2MCGzTAZJCVt6TjYsBnGhpS0peOYuGHXv/YzqLIkUA8I9Fef0/7VyrbZTVPcpYXo88eGpWNCo/1zR2UU9nwHklems7Bbp10wANCap91Zx8sGV9AlPY2tWYmjwMLy7Ut0PhXnaVcVAMzUaadZ2bDofLfn/8e8B+RzepqPTuh87OyZ87G4TV+jukUYo/nyGnLy6kck6RixzjDThhTxG4FJN3xUE6jYsBnGAuHFb3Ek1A2/4zjOQ8RLFCdi4lVxhplK7HrFbwP4RH3v64aHx3mMZAIVl9gMY4NHbPEk0Q0fyQQqNmyGsWFMQ0pd131PVdWX1AxKK7gqzjAWCDt3F6kbPugEKgousRnGBrsS+wP4xnqmGz6OCVSJSmzHcdaCkznUBI8mgLX+3OvQ8WR8kNlvfB9rFZ4Wrjyl/avFnWNj2sWx2UdL4dVCE649CfT8u2fyY4fF/KM4WTL4qS/TT0v7ku5D7pZ6aNf8fbWlo/ApGgt1cz5Wi7S/v9Wl8+EwMM++I3NoeXq+tA4MUy4NPmrA4KcGMPMt4c//vvHyGsKiyu267oZzSlg3/FMAN1X8HcdxtjHkBCpjbqlEvI/TN4jWZe84zmqE0ZPxDJMJ7HXFz02QGscEKmNVPGLZT9NSnyNbCpRhUkvK52PbdJ6ZuuxHthQow6QVaiBKGph4r/iNl+tn368u6e2s0iy9yFr+gI4HANGmxzCb8Mq6JM/KtboeP0OnodUwS/q0FuiKU2uGLga6Jb3teL2g9wuUDbe53jMJtwFl0OPNOz36PxxKP59eRu1cyXZcpNOQr5v7KkpLdPW4SlrjM+P1NVJSMsdhY9impT4HWgr04TeHseFUdp4B+Or3vl5Xd6FKnn982fzyOezRD+1R0dB5Vj3/8vqs7edjDYbOs/wIOs8E/Sg1O4F8lMDvQppn+x36BVw8NL8gq8/tO8+uDdp5lvLZXTbuLlOXfWQ8w2QJm7HiF4nRsMPLfhJLfY58KVCGSS3T3nkWs+znyLrsa4/8amK5WtDC+e2DqFPOEIfmKqRx9bQq3baTlVAWlfLavrahDX1kmEsNAMcrhqr2Cl1Nnb2k51NV5jEr/H2XZ+l8nC2aBIqBtkc/KgdtuslxdOjncytfxFGoXZ83zaV+bK76kn5qADNfE8/LwFXxwY6/aCbeecYwU0lKSuY42LAZxgKbkWcXCRs2w9jAhs0w2YPb2AyTQdiwGSaLcFWcYbIHl9gGgr7q3H5R913vH0ac4SO75nHgomLwU0cMGQ3SXtTjuzMlbd/RsmGJ20tm/2trQD91mJfn9eGZS50CXi76+5Yr9PlJeHJM+6n3juh89Jr++V65CK+lX6/2jPb317bMllR/RPvjC4+oNbAHFRYf7PA+o9QyoGAFFYaxYBzywwnliRPBhs0wNtgNKb0wLYOJV8UZZhoh29jxrbML0zJIVGJHVBnuqO1ezPEjWaaEYdKKkPFbGkgyu2sdwEeh8AM10WM1Rvt4JMuUMExqsVswYKRaBhQ2mmer8Ov+m4jWPn7Ldd0b6lyGyRyWJfaFaRkM3Hnmuu79wLTMNURrH49kmRKGSSs2hn2RWgbWnWeq3b0R9eN9/5vjOLccx1nnkpvJHJYDVC5KfniYXvF113XvhncOukzJqzdmzr6vXK1CtH1tLdml9cRE3lzhkHX6Gt0GHX+ypGfR4kIZrXn/d02L0reWzU9AeYEWjLhUo2/TUk7/Dw1R1vYVO/Tgkv2ueXCGPJol44tHdTK+0Pb/w/fyZRTa+sCdmmEcT61Ka64BQHGB1rfzYNZNS0paOsnisDJsx3HuBErldbU4QMN13SZOq+b9NvkNnC42EMtXnz2PDcs2fTNFwaxcKRv0A9m+TKtuRCmDfPGdr/JyUKAN+2iGvj4AVEGPDDsq7pDx3cruuX2PK1tn38s5emTbdrtGxgPAly36Gltd+kkv7+v5+Nt9XWBx9jH9+/NfmkUpy3+gRpYB3nfxP7L4l9eM1w+S9iGlA2ueqfr/PcdxHjqOE3zigu2EN9XxQy1TwjCpJWuaZ6q9vBBx3MjbCQyTVtJeYvPIM4axgA2bYbJISqrccbBhM4wFwku3ZbNhM4wFmXR3jTYFAZdVPq+FRZFOnslHDQCdZdqVc3TFIHR/RXcctBZyOJIBP/YK3dgqLJsXNXipsUfGX63S8XMF/Tdquba276BH+6mbbXM+7uzT+ZjboX3E5cC6WkUILQwA1eeGtbue0uuPAYDcaZLx3ol5YYSkcBubYbIIl9gMkz24xGaYDMKdZwyTQbjzjGEyiDBPAZgoLGbIMDaMaay44zivU1oGSWXH2LAZxgLhydjNloTyw4lkxxJVxR3HWQvO0nIc557runfV9M1zEz4GET0P+qJlpaSFvSrtG+00aJF6ADhaoa9xeJV+tx1f1m9Uqy5xXAzsW6GnE640zGL9l6v7ZLxpYfpe6P3sQWj7dju0n/rp4QwZDwCtXTqvqzt0PlZ2/DwrlaUWBoDyNj1FN3dg9mN7hmm+yJmn+SZlTG3s2wA+Ud/78sPh2ZFvqYlZJAOLGSpi3xqjFD1nmLRis2BAApLIDyeSHbMRMwRoscKRiZ4zTFoZR1U8Ca7rvqfsbilGIRiAfa94X3Y4qqo9MtFzhkkt9mt3RXV6bfZruCDkhweRHbMy7FGKFb7yo8bZ95WXdd0sWaLbRN06Pc4bAFrz9DVaiwbNsrp+B6+X9XHXMk9n4ULP3D95uUVfY643Z7yGdnxrUQvnO/Q477xnvv6VEn2N0jz9H6pX/Hx8aeH82PWaoM8vzJktybtmmFvQMeumJUX07CzbIELyAQBHfdfkiQeVHRvYsBO8NQYSPf/qd83Y8IV0nhlaI8el8zfwN4d+R443T3dsXc2f1yML0yvTWl3t8uCrZT6rfXv2feuENtyHR+aG4bM2/SBXd+l8nt3Sz/9iS+8Mm/uSzsfS13QeAYD35BkdfxLf0bn0168ar68xhhq367obzilR8sQ3Vfwdx3G2YZAdsymxI98agbdK5FuHYbLEuNrSo5InHljMkBArHLnoOcOklbSv3TWwmKHaNzLR886S367uzZa1cGeOTt7xornCcbxCC1Yfr9B3ondZr755hSJ6M/6+5QW6mnylTs+lBoBG0TxnmyIsH1zulrV9W4d0VXxnzyw/nG/SeV0ytDhK+351vzAvtTAA5I8N7d+eubkgyvS883zV3HRLCk8CYZgMYtt5dlGwYTOMDem2azZshrGBq+IMk0UkGzbDZA5uYzNMFkm3XbNhM4wNwku3muHEDft4xR/v3ZovhsL0+JmTZcOiygBOlulXa/dSm4xfWDjUwnVZQEP4+67M0HOpF0rD+agB83zqJ8f6UsG1ThVPpL/v6T69dnVv1zzmvrpL53XxgM7n/IlvCLm2p4WTIKvmNbyNOvQV6n8OpnXEVXGGySLcecYwGYRLbIbJHiLlJXYiMcOgvJHjOGuO40jHcR6q7dyc0KRKigwztXhe/DYklJyYScW0j43m2aLrusJ13RsA3gBwL+K0REqKDDO1SBm/DUGMxmA/LrGe4MCaZyG1FMd13SjjpTTRGGbqET0Zuw1DjMZgn8R6gta64urN8mFMdCIlRYaZWnpe/DY+EusJDtN5diuuRB5EE+3GNd/HenVJ91V267TvtD1n9mO36/Qb1CvQPty61DXTrkndJ7zcpiWBGqB9yABQNKwXM+vRt6kc0jRb6ep+7RlDGlpl8zzlkkE7rmIoqcpF/z9euXT+90oNOh/FkTmNwmBUskTlI6ngFXEx6v+an8txM4xhR9bvB1FSBIDPWv4AjnZH4PNA+MRg2K2CudrTq9EDUOZmaSH6q7O6UEKjs4Rm0X8ICpXD8CkanQK9oAAAtA2G+zg0ACXMdyf6+T+UXfy+6//vpmEASvGJIYEAqo/pe1F9QhtVZVt/eW3+Uc+34r5hwYCuuSSUeboC2pmJv8ZfGK8eguwki38JGlRKTSTWE7QybMdxzhmrjZIiw0wtltM2B1EX6mOjJziw5lmAcAM/qHkWpYnGMNnB68VvQxBjbwPrCdpqnm0CeDu0z0rzjGGmkvGplEbZ28C2xSPPGMYGnt3FMBmkN1yVe9ywYTOMDSkfK86GzTA2jHcgytBM3LBPFn3/aGdW4KTnh1uLBjH/RfMiazMNWuhgeYb2Q8+V9DWlarKj7Svm6CrZYdcsENBs00IKjw/ohen39vTzj/Il7PX8fbldw4J3++YBFfkT+l6YVsDwyv5vyKLQwgDQydEDVLyCOY3dCu3kac+MbuCIlGzYDJM9uMRmmAzCveIMkz0k94ozTAbhqjjDZBDuPGOY7MFVcYbJIJIX5aP5295jP/AceC0Y+dRwsik+EbT/9HTaa5A2fnBu33BcN8T/qVHMfu9c+KdDpCcSekq4OV7jOa6+Zj4qxXx18x++/woVf2EpiUHIlA+NYxhmcKw1zxiGSS9s2AyTQdiwGSaDsGFnhLB4fNIVIy6KiPTxajFjJBWGnbaHMIo0P4jh1SMGWTHiIohZ3YJXixkjEzfstD2EBKl9ECNWj0i8YsRFELO6Ba8WM0Ym7sfG6UP4ifrefwjTqGz6lhKamwYSrxgxQVZVSb7WX2CCGR0TL7ExHQ8hwMsWjRTXdd9TpfWSMnBmhKTBsKeCKXsQE68YMQkcx7kT0M02rhbDDE4aDDvVDyEwlQ/iB/DTSK4YMSFc+Gm6ocLMCEmDYaf9IQRS/iCGV48YZMWICaaPV4sZI6kYK65cSJsAVtO6ikh/sUGcppE7e5hUkwrDZhhmtKShKs4wzIhhw2aYDMKGzTAZhA2bYTIIGzbDZBA2bIbJIGzYDJNB/h+DTFk57iX0fgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAADNCAYAAAB+UQYWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGsNJREFUeJztncFvHcd9x79LSpQoinmkCDKRKDkKabh25B5CjYBeilwotEB7lKO/QNIp6KlScip6KCr61qP0H9hxgF56MnPqKc5KMJAoRhyIsNtadmRYpCI+mgr1dnvgLN/MvN2Z3Zn3+PZxvx9goTc7b3b3Ufub+c38Zr4TpWkKQsjRZ2zYD0AIORxo7IQ0BBo7IQ2Bxk5IQ6CxE9IQGmHsURStRVE0M6Brr0ZR9GHB+c0oim4q525HUfRh3rMUXadPzzcTRdFKFEWPoyi6F0XRkpK/FEXRL+Tna9lnx/U+zCnbt+eX/1+3lXus9OO6TacRxg7gGoBV2xeiKLrmc+E0TdcBbBWcfx/AM+X0QwDvpGla9P2e8yHISmUmTdOtNE0fyvvfS9N0Q7nvBoAb8vMHrmuqz2mU7efzv2c8n+jTdRvNkTd22Sp8AOC65TszAK4O4Pa/MO47k2foA+Smy4BlK79k+84gylZk3bcyJl2OvLEDEGma3oHRskuXelW62QKAyF4ow1Vdi6JoTSl3LYqim6p7XoRs7XI9Ctt1iu4v3fGD55au+Wp25Nxm2fWMkjXzRBRFD+Q9Zoy/lbWs/N5t1fWW19HKF5y7LX+H9ltk6z6IyrhRNMHYM+LMIOQLtiGNcUb++yxrBWU64172Qb7AS2ma3gdwq+R935eGvQpgvcx1iu4P4OcAHsr8ZUivQaY34IE0JM3bkM/6jnw+8562skvyeweelOx7Z3/rDWnQeeduKvdRf3/GGZ/fR7ocaWOXbuaybLE3ALwjsy7LNNI0fbfs9bJ+rzSGZ67vSzJX/sCF97wOsO8yz8jK4hsA/w7gahRFDwBUHoAscI2XsF8BzShp9Z428n7LFXQrog2Zzjt3uaC87dqkAkfa2AGspGl6R7bYqiv/GLKvqYyMb8n0ipqG0ieVrU/WemWViZU8V77kdXruD+A32G8RHwK4D2BV/r7L5j1KknffbNAtc8/Ne1ZlQ7nPUna9nHMPwNZ7oBwb9gMMCtlq/jyKouxFXQJwJoqi22maviv7wtnXP8C+O3kNXRfyN/IaMwBWpUFuAFiR5x/KzzPy3xV5nzzuQ3dNy1yn5/7yuW9HUZQZxRXjN7j+HisArsvfch3AM1m5rSiDbSuQRhdF0b00TW+p94yi6Fn2nPLSatksTLaa/aY0Te/I8sB+5fuuvI71HPY9lvuHPKB5pIm46u3oklVsw36OUGRlslImNEiKOepufNO5f0RCVjT0PkBjP8JIF3grGtDswcNA6T6RQOjGE9IQ2LIT0hBo7IQ0BK/QmxDiGvbjwCtxHPeM9rryM372s5+xD0Fqw927dyP3t4C7//bTdOvFlO0rn9+9e/diXx6qj1Q2diHECgDEcbwuhFgSQqzEcfywbL7Je2fnDz5fmjyFR9/uHKTn3rBP2PqH84+cz/uP3/nYmv/2hP3/90R0XEt/tflDfG/29wfpL19tW8v/9+6i4wmB//zavoLz1xsXrfknP5nU0m+1TuGT592/4+wfOtby0394bn9AANH/PLHmd/78Z3v5Y91X7eLlc/jsgX698e9911r+5ev2fADY/KsT1vytN4vblutffu28/sF1XkzhX//p/cL8f/mPn3y/9MUOER83/jq6s7s20Dtzy5VPyMiTIC086oqPGz8DfZ7yXMV8QkaevdTuLdWRoU+XvTR56uDzaxO6Gzad2P+grfYF5/V301fW/KeOv8DxSP/C1rZ+z63Ot9byY3vuuu7cywVr/qXx09b8iZb+dzs/paenFhNr+Ulr7j7RGfv/RbIzbS8/Pn7weeH13inwY7P2qQB7Z+1/AwCYnz9uzW9P9q/VrXMLXoSPsW+hu2BhBr0roVz5Gmof3UzPjdmnRb8x9b+uZ8XJ7/zOmr9Qsc8OQOuzp44+e1Kiz/5k124ojzqnrPknn/caotZn/6IfffY/WfOr9NkB5PTZ7c/4ctf9qm5OOPrsc8UG+rbz6jp7sFegdcSnz/4e9BVL6wAghJix5RNylOikaeFRVyq37HEcPxT7rALYUkbafwXgsiU/F3XEfTrpaK25a7TdNdIOVB9tNzFH27c632qtuWu03TXSDlQfbTcxR9unFhOtNXe13K6RdqB6y22ijraPzc70tOSu0XbXSDtgH20HHNGdL52X19hriBuPOI571jXHcXzZlk/IUWJv9Gx9+AN0hIwiHZSaf1MraOyEeLCX0tgJaQRs2QlpCHvp6K0ho7ET4kFnBBeM0tgJ8YAtuwdqLL3VvqDNigtdsQZUj6ObmHH0sb05bVZc6Io1oHoc3cSMo08a50JXrAHV4uh5qHH0vbOne2bEhaxYywhZJfkSbzmvr9KhsRPSDPYw7v5SDjatB7k8/AG6mnvrcRzfEkKsxXF8RwhxM2QOy+hVT4TUgL10vPAoQtV6ALCVpRXOxHEcxXG8jP3di7KNOm4KIR4jUHiTLTshHngO0F0HkO1hn2k9HEwnl5VAhlBa8RtxHAdLadPYCfFgL/UynVJaD3JdiSqFsyTPWWXeXNCNJ8SDThoVHn3gahzHByvC4jh+V7b6c9LovfAydiHETXn07Ost89ey7/k+GCF1Zi89VnhYKKv1cNCXl3aW7erzDfI34yxFZWOXNcu67E8sFdQ0fRlQIKSudDBWeFhwaUFACGEac4yuJsSyTHvh07IvoSsiqW69q3IjjuNlY8CBkCODz2h8pu1QoAWhsmGU+Yls3R+79CFs+IhXqHG+FezXVialBxRWFY24v6CDCSW90LaLBW3ulJAqcmjEPXJoxH30fFlLz+218I1ybvsru37cj564RRemtl0acW1rvqkPtzCdaOeSky59OPdWcG6NOHv+9mJ3clOrdQLzi/pkp1Pn7ZNm3lx84XhC4MqsfXLQFYse4X85r65jM2obJbQgNgDccpXxwXs0XsYIH+bVNJmBCyGuCiFWbS38ydbvCtM++nAmoRpxefpwT048Pfjsow9nEqoRl6cP9/nH3XOhKjNAuEacqQ/36de69qBNHw5w6xECbk1Cux7hXzuvr9K0GXSrcRzfMU/KQblnMi4YNKBASF3xbdmHifdovNJ6r8p/Mz+ubwMKhNQVzwG6oeI7Gr8mhHgshNhUsn4F9HdAgZC64jNAN2x8BujWAczmnKfgJGkMScP67IQ0ljq34EXQ2AnxgMbugSpA8fSYHm4LFZ4AwjdxMMUnLo2f1sJtocITQPgmDmZoLdmZ1s6FCk8A4Zs4qOIT7cm0J9Q27O25f+m8uk7TQm+ENBa27IQ0hMSyuq2ubT6NnRAPbC27e4L0cKCxE+KBrWWvK3X1OAipNb6TaoQQ14QQq0KI2wX5PVoQrjJlobET4sGrZLzwKKKE4CRgaEGULFMKGjshHiSICg8L17GvVgN0BSdNTC2IMmVKMfQ+uxpLPx4d09JVN3DII3QTBzOOPtE6oS1brbqBQx6hmziYcfRofFw7V2UDhyJCN3FQ4+jTSadnyaorju6KoQPuTUPKLIkuy56lBbdQRnDS1IIoJVJZBt9Vb1aNuX71MQipK0kaFR4h9EtcMg9fN75QY66ffQxC6sqrdKzwsGAVnCwQlywrUunE19htGnN962MQUleSdKzwsOASnMzTgsgt44OvsS9Z3PS+9TEIqSs+LbtLcDJPC8JSpjJeA3RVNOZcfLX5w4PPW9sXtLwth1jkmEMsEgDOvbQLQl4aP23Nn2jpA1Pnp/T01KJLLNKNKRhpkuy4BCP1waKF189oabdYpP1vAADz8/bBrfakfYBuOun+xnPJVE9+q32h55zKrkUsMuOp420+HvVvPPpV4tdOlhCczMsfjuBkCY25Sn2M783+vjAdKhYJ5AtGqvgIRn7yvCuWGCoWCeQLRqr4CEZ+9qA7wh8qFgn0CkaaVBWM/KORDhOL3CdMoLTals1NmUGXqzGn9Dv61scgpK54DtANFR9Zqody1PAZdI25XwG4LPNF2T6GGkvf6nyrteaha9GB6nF0EzOOPrWYaK156Fp0IFzq2Yyjj83OaK156Fp0oFocPQ81jt5qX+hpyUPWomf0Q/+gLKPYsvv22Sv3Owg5Svj22YfJ0GfQETKKNKZlJ6TpUJaKkIbQoRtPSDOgG09IQ7C37PaJVsOCxk6IB6k9EllLhm7saix9bG9OmxUXuhYdCNd1N+Pok8a50LXoQLiuuxlH3zt7WpsVF7oWHQjXdVfj6Lvpq54Zcf1Yi94P/YOy2Afo6tmyj94oAyE1wHc9ewkNupvyWFPOWfUjykJjJ8SDJIkKjyJcWg9y1um6nJSWKdYAFv2IKgzdjSdkFPEMvV0H8KH8nGk9qNPJl+RxX+Zna0xuyIVnQdDYCfHAc4DOqvVgTDNfwf6iMqBXl86LytWTEGJFCJEKIR7L417Od/rSxyCkriTJWOERinTvD4Qr+qVL59Oyn4njOFIeaivnO5mW1i3fByOkznhG3spqPazGcXwHKKUfUZrK1ZChSiPiOM4bNLBp1BEy8qRJVHhYcGnQQQhxU1GCWkWBfoQP3j6HfJD3C7JtGnWEjDxpGhUeRbg06OT5Ndk93lTKaLp0vs8cMkB3tajlrqJR99H//f3B57m9Fr55vnyQ3v7Krh/3oydu0YWpbZdGXNuab+rDLUwn2rnkpEsfzq7/BpTRiLPnby/qE05arROYV86dOm93Ot9cfOF4QuDKrH1y0BWHRtxC++2Dz893X0PLePM2d+yv4qZDjxAAHjk0CT9S3i2TCbg17lRsITYbNi0IaSezZcr4EGLsudPbqvYxnpx4Wpj20YczCdWIy9OH+/zj7rlQlRkgXCMuTx/u06+7OnlV9eHyCNWIM/XhFmb1GXeuGXIuPULArUlo0yO8iDOFebnP42nsw8R3R5geA3ZoXxNytEgtR00JiROYA3OF2tcB9yCklngO0A0VXw26DRhhNWrQkSZhG4irK5xBR4gPNHZCGkKN++ZF0NgJ8aHGffMihm7sqgDFpfHTWrgtVHgCCN/EwQytJTvT2rlQ4QkgfBMHU3yiPZlq4bZQ4QkgfBMHNbR2PDrWE2rrh/BEiNjJRbhFRlTSeupTWBm6sRMykrDPTkgziNiyE9IQ2LIT0hDYshPSEDxDb3Jm6RYKVGfy8l1lykLBSUI8iJKo8CiihOBkT76rTBVo7IT44LcQ5jq6yk6Z4KQr31WmNKXceCHEirqgxccVKUKNpU+0TmjLVqtu4JBH6CYOZhw9Gh/XzlXdwCGP0E0czDj6dNLRlq1W2cChiNBNHNQ4+lbn254lq644uiuGDrg3DbHP26gWZ48GIDhZkO8qUxpnyy7VM36hpCu7Ir4PR0htSaLio6Y4jV0arbqc1ccVIeRo4efGuwQn8/LLilQ68emz+7gihBwpoqT4sOASnMzLzy3jw9BDb2+1unPhz0/pfdepRZd+nBtTQ84k2XFpyI1r6YXXdfkit37caWs+AMzP2/u77Ul7czGd6L/xXDKlpVvtC9byuw79OAB46nhTjkf2L2wpGnLt9sWe/DGHfty5l3Y9QmB/bYWNiZZlbKRqe+nRZ4/j+KHYJ09w8nJRfkGZyvgYu48rUsgnz3cK06H6cUC+hpyKj4bcZw+6g36h+nFAvoacio+G3B+Vc6H6cUCvhpxJVQ25mZnfaukQ/biMEM3CN5xX17GF2GzYBCct+X0Rg/Fx431cEUKOFJ5u/FApMxp/bf8fcQ1wa19b8gk5Ooyg4KTTx5SS0B8Y5yq7IkWosfSpxURz3UPXogPhUs9mHH1sdkZz3UPXogPV4+gmZhy91b6gue6ha9GBanH0PNQ4+tjeXI/bHrIWPSNI/6BYUj6XOrfgRQx9gI6QkaTGLXgRNHZCPPCcQTdUaOyE+EBjJ6QZsM9OSEOgsRPSFOjGE9IM2LJ7oMbSJ4106Fp0IFzX3Yyj7509rU2BDV2LDoTruptx9N30lTYFNnQtOhCu667G0c+9XOiZ/hq2Fn2fIP2DinF2tuyENIRRbNkpS0WIDwOaLiuEuCaEWBVC3C7IvymPNeXcWpZnuzaNnRAPBrEQpoQK1CqAdTkdfUmmAeCmEOIxdJGZHkoZe85Ne2oXI79UTUPIqBKlxUcALpWnJeXcBrqrS2/EcbwsK4lCfDToimoXlVI1DSEjS2I5/LGqPMVxfF9ZZLYCIJafl2yuf4aPBl1R7aJSqqYhZFQZUMteCulpP8yWj8dx/K60tbmCxheAx2i8sXx1BftiFSZZix+0gwUhdcXXqAu6thtZPx3lVJ5W4zi+o1zvmVyK/g3yG18AAaE3s3ZRUbatuSqEWGULT44cnu66Q+vhPQBCftZUoOI43pKfbyr2tYp9Vz7zvJcB3Cu6eEic/aB2UalS0wDARUUQcmE60QQik5MuMUi72CNQRhDSnr+9qE84abVOYF45d+q8vYp/c/GF4wmBK7P2yUFXHIKQC+23tfTz3dfQUv5nN3fs/82bihhkEY8cgpAfPbfPStn+qisYuTs2ie1EF5D80RP75KSpbbd1TaJtzXeJj1ZhEO66S5BSnl8TQtzBvgfwjixzUwjxDMBjmzKUl7GbtUscx+tK7VO6pgGAzz/+U2E6VGUGCBeEzBOD/PTrriimjxikSaggZJ4Y5MJsd9ZdVTHIPEIFIU0xyEcd/Z42MUjALT4KlFE2KhYfnf2beef1tWsNaFKNTQVKesizZcrkUVmDTqldHgshNpWvqhp0P5Hft9Y0hIwsTdCgs9QuXhp0hIwiozhdlnPjCfGAxk5IU6ixu14EjZ0QD6Jk9Kydxk6IB1SX9UANryU701o6VHgCCN/EwRSfaE+mWrgtVHgCCN/EwQytHY+OaedChSeA8E0cVPGJidaJnlBbkPCEJEzspB6ht0EydGMnZCRhy05IM2DLTkhD4AAdIQ1hFAfoKEtFiAdRp/gIoYQGXY8KlKtMBo2dEB8GMDfepUEn0VSgSpYBQGMnxIsoSQuPAFwadECvClSZMgBK9tmFECvq6jUhxFocx3fkUteeRS9yxdsWSijVqLH0aHxcS1fdwCGP0E0czDj6dNLRlq1W3cAhj9BNHMw4+lbnW23ZapUNHIoI3cRBjaNPLSY9S1bdy1PtMXSgP0uiyzKgPrtVg05iqkCVKQPAQ3BSUigoWcWtIGRUGYSUdBnK6s3lUWaJ67oQwjTqG3Lpax7XAXwoP2duBde0kyOFr7seokFXoAJVVrfOO/RmE5Qs7VYQMrJ4uvGBGnR5KlBxXpk8vIy9n4KSFy+fO/i88PoZLc+tH3faef35eXt/tz1p/1+bTvS+5blkSku32hes5Xcd+nEA8NTxv3A8sn9hy9CQa7cvaukxh37cuZcL1nwAuDRu/1tPtBwacotd//bsXO93Xds2ltGPS3ZcmoXjzmuUJer0v9Pu0qAr0psrKNNDZWMvIShZ2q0AgM8ePClMh+rHAfkacio+GnJ/VM6F6scB+RpyKj4acjMzvz34HKofB/RqyJlU1ZB7/IUuDhmiH5cRMkB35sffd15fY0CTamwadJb8/mjQ5T0Puq7CskxDCJE1w++hWwFY3QpCRpUBhd4GSmXBSYugpCo4iTJuBSGjyjB3hPGlsuCkPFfZ1ShCjaWPzc5ornvoWnSgehzdxIyjt9oXNNc9dC06UD2ObmLG0cf25jTXPXQtOlAtjp6H6qZPotdtD1uLLq/RB/2DstS5BS+CC2EI8WAQA3SDhsZOiA+jZ+s0dkJ8oBtPSFNIaeyENAL22QlpCqNn6zR2QnyIktFTnBy6saux9L2zp7UpsKFr0YFwXXczjr6bvtKmwIauRQfCdd3NOPq5lwvaFNjQtehAuK67GkePznR6pr/2Z3vuEP0Dty699jx04wlpCBygI6QhDKhlt6k8SSGYB+guc12P4/iWSzkqgxp0hHgQpWnh4UsJlaczcRxHcRwvA3gHwJo8X6gcpVLK2NWbCiFWhBCpEOKxPO7lfL9H7paQI0WSFB/+WMUjDd0IEcdxZtymCGUuPhp0RbWLSqmahpCRJU2LD39KqTxJm3xfObXUF914WVtsGOmD+yq1i0qpmoaQUSXqpIXHIXBVylQBKC9C6T1Al1O7qNg06ggZfTp+7nqI4KSC2q12KUcdEDIaf7Wo5a6iUfeDH3S1zb77XT3eG6ofB/RqyJlU1ZD7y85rWjpUPw7o1ZAzqaohN7fX0tKh+nGAriGXRxUNuYWlXm3BfujHhWkWVouz29314rkXgYKTEEKYxpwnQplLiLHnKiJUqWkA4OHE3sHnN44dx6dKemvuL9YHmLtgnzADAJcdk2b+tqL4xNNNYGG2e81Q4QkA+K1r0syLi9b8X3+h68NdGj+hacad/MShD/eHF/YHRPgmDuqkmWRnEZ/FX+jl+yA88XLSfo1N5d0y+Tvn1Q2sA3F+wpYuwUnlqxtGmR4Ryjy8jD2ndnHJ3RJytBjQEtcSKlAbAG65yuRRWYNOwRyYUzXo8jTqCDk6JJ3io6b4atDl1S5eGnSEjCQUryCkIXDVGyENoVNfd70IGjshPnDVGyENwXNSzTAZurGrAhTtyVTbey1UeAII38TBjKMfj45p50KFJ4DwTRxM8YmJ1glt77VQ4QkgfBMHNY4ejY/3xNXDhCf2CRI7+dJ5eY00pbET0gzYshPSEDgaT0gzSDkaT0hDoBtPSEPgAB0hzYBuPCENIeXc+Opc//JrLf22mnDEPl/iLef1f+n4zi+dV8jDfd8qXHTm22PYMPO/Ad5Q08uO4q58AMB8YL7OmR9/v9L3K4tL5FExlm7h88v//APbD/i8b3fqI1E6gtP+CCHVoW48IQ2Bxk5IQ6CxE9IQaOxHBHOrICHEtTIbBxwWOc/HXYMOmVoYe91ezDzq/HKau/aU2DPsUMnZVQjgrkGHztCNvW4vpoXavpzmrj1w7Bl22OQ8H8Bdgw6docfZsf9ifig/Zy9mHRVpb0jxzVGg1J5hQ4a7Bh0yQ2/ZMRovJlBy8zxSjrL7k5H+UQdjHwlG7OUsu2fYUJA7mGT7EDh3DSL9oQ7GXusXExjJl/M9dJ/xYM+wGhGj+0zLMk0GTB2Mve4vJlDzl9PctSfbhSdnz7A6PR93DTpkajE3XoazNgAs1XU3mWzDSuw/IweUyMhRC2MnhAyeOrjxhJBDgMZOSEOgsRPSEGjshDQEGjshDYHGTkhDoLET0hD+H+CaBd+uxWczAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO8AAADNCAYAAABQFCeEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHONJREFUeJztnUmTHMd1x/9Zvc707EMMAIIkxIGoEGV6iWHyGwzDN59A4RNwdHL4ZFA6OXxwkODNR9KfgCIOvhMKH31gEZLDlERKmiFBiyQAAZh96a3Sh8meqqzueq87p9HTNXi/iIruqqwlu7pe5fbeP5UxBoIg5I/grDMgCIIfYryCkFPEeAUhp4jxCkJOEeMFoJS6pZSaS6xfV0p97HmuVaXUJ0PM26pSat3m8bpSak0ptT6s8/e43k17netKqZsex3fyeP1p5E+IEeM95jqA1c6KMea274mMMXcAbHH79ftw2/PdBfCRMea2MeZDAO8kXzbDQim1CuCuvc5tANc8TnPD5vHOcHMnpHnmjVcptQLgNoAbI7zmHIA3T3HsHQALQ83UMel8feB7ImMM+wITTsczb7wAtDHmHSRKXsuyrbJetyUSlFIrdttqYtuaXV9LnzhZhbbV3ludawLQndJXKTVnq6s9z5Ni1RizZYzZsHn7pFPFTa/3yl+vfTrY0nZVKfWZUuqmMeYu9Tt6XH/F3rfrifS15G9K/s4Bf7eQxhjzTC8A1uznJzg2jM72TxLfP7Oftzr7AFgGcDOxvgrgpv3+cfo8dv9bGedPnvdWjzx+bPe51blGIm3dfs6l14n8Ocf0uN6qveYHffyO9PU/tp8riet17t8agOv2+03ud8tCL890yauUWgZwzZYUGwDeSiQnq30bdt93AbyplPoMx8bxhj0O9vMNz6wsA5izJdfjjH0+Msc1hE5puGy33wWcampyPSt/6WOQPKcx5o4x5i2bL46e5zLHpfZdW0N5Yje/3smPMeb9Pn+3kMEzbbwAVowx75jj6mK66pzsEJozxmzguJR4xxjzut13A/EDvgzg0x7X2Eqkd223D+6nADbsA/8hlWFjzB3b7l1hf11/+UuymngpJPOe/N6PQcNWg5fNcYdb58Ww3jne/oa+f7fQTfGsM3BW2BLhF0qpzsOzDGDBtvXex3Fpu4pjI+60Vd9QSnVOcdsctztv2m0rxpj3rTGuKKVW7Hk/TZxnVSm1bF8EG7bEv2OPu6mU6nRCnfTU2mNXANywBrAA4GcA3uqkda6VXjfGvNMjf6up/CV50kmz13k3kdb1O+w9S16/89s7L7bO97uJ69/q3EPqdws8yrY3BEHIGc96tVkQcosYryDkFDFeQcgpYryCkFPEeAUhp4jxCkJO8Rrn1Vpfx/Gg/UoYhu8Pmt7h5z//uYxTCWPDe++9p/i9gPf+7R/N1m6N2uXee++994OhZIpgYOPVWq8AQBiGd7TWy1rrlTAM7/abnuY/yxdPvr86O4nfbx+crJf0JpmXt5Z/zeb3H2Z+Q6a/WiqR6SVVcNbvb/4El+Z/d7L+sL1PHn+3zgf//Nfuq2T6fz98mUz/v+/ca7xWnMLnrb2T9dL3ZfL4ye/5Z7Z2v02nf3tEppe+eXTy/epri7j3uesN2b7/gDw+eIW+BwBw8AM6SnL3xezH/e/nHrLn77C1W8O//tMvM9P/5d9/erXvk50Cn2rzDcSuchvojsbh0gUh90Qwmcuo8Kk2zyF2NAeAxQHTBSH3NA1dExkFZ+7b/Ors5Mn3FyYrTlrBROSxM/svsec/MC0y/WGBvgUF5VZOtvZedNa3I7q6GDWnyXQAWDi8QqZfa82T6bPFKWf9pULVWS9M0b+xsshXm6sF+r+oVOnmR3EmLpGWrnbfk+hKoWtbEvU8LxxSXyLboTi8QFyjwZ7eYZQlbBY+xruFWMVhDt2hXFy6Q7KNm14vqW0yI39T+4bLKyZnPifTlwZs8wJw2rwB0+YN+mjzPlG08azv0W3W/2t1H++0eZnjJx+Pos3rPgYDt3l3KmQ6ABzQWcCuyX7cXx5QVKgJ+mU2CnzavB/BDTO7AwBa6zkqXRDOE21jMpdRMXDJG4bhXX3MKoCtRE/yrwC8TqT3pL4Yv8GaExHqxXh9qkRXebdbE2x+v2hcJNN3I7pHu2nc91urofDFYbztmybdsfi/By+S6QBw9wm9z3ePZ8n0whO39hBMFlA4iLeVt+mStbTHP3ClA7qkCRpMGzD5UBvjrgMwEZOHJv0sAEChTuexyJTMg9DMabUZYRh2BU6HYfg6lS4I54nm2dvu2XdYCUIeaaMvf46nihivIHjQNGK8gpBLpOQVhJyS7sg8C8R4BcGD9hgE5InxCoIHUvICiKrx2JwpG0QJl8h6k87evQPeeylQdJ/+nwqXyPR6yitn9vBFbCf+t++OaNecP+08R2cQwLeP6HNEf6mS6dUn7oNUbAcoJzJZ2aTvQWWH9xYq7dHjrMEB7V9o6vX4e7PlrAOACug2pGrzeeTGmotHwxvfaYvxCkI+aYL2xR4FYryC4EHTiPEKQi6RDitByClNIkJpVHi9PrTWa3a5lZF+q7PfaTInCONK26jMZVQMbLw2WuiODT5Ytutp1rTW64inlxSEc0XTFDMXjo7OW0bada31qtb6ZtY+HXxK3mXEulTJKSSTvB2G4bUwDCWWVziXtBFkLhS2sPs4I+1EvBHAFmXkgIfxhmH4YSLkbwVA2GO35X7fHoKQR5qmkLlQWMPMqpEOJN7o3eq2b4W7vYLtO1rNWus3tdarVAn810Gsv/SSqriD9Qe09EnwaIbN5+M9OtB9J6AH9tuR+2dE7Wk82IkdOw4atMTMDPMbAGBqj5HiadLv2HLqebkSVFBKbKtM0g4OpTYvplYs0HkIarQDRLQQX+PCCzWYRsoxJaL/y2COFiQAgGiW1rBqzNH3eRCe0lDRQOKNp+kyWw3D8J30RttJ9SQMw9s41q8iZ1L/7eFB9vpkHRQz6pDN5ELxgEyvFGjPoWbU/Sd9VYr1l7YYNY/NaJJMB4DoiDbwyg79oFQfdW/78lH8uycf0cZbfdwk0wGg/Ji+18HmLpkePXrirN/79X13B+YFElzklTTaS/Q5jpayPdUu/RV7evdaY+Bh5d3bnChdV+1n51UaItatuobe1WpByDW+1WaGgcQbfXubb2mt17V2pjT4FXCscQXgp3bKk3VOw0oQ8ohvh1UvfMUbfQTo7gDoEhIWDSvhWcK3hLWFmtZaX7dNS8BTvPHs3UQEIYdEnm1ea7C3U9u8Cj4xXkHwQAITBCGniPECqDyOqx+l2QCVRBB5s0EHoW/t8+N22xP0UE1QYoK8U76qk4UpfNWOt7WPmDl29vhbXNmkq2CVLTIZlW33N5QmI2dbeYceQint8BP1qD16qMjs0dO+RIex4rlpNJx1AAhq/JAah2rQw0nF/eFNDjYOQ0VnbryCkEek5BWEnBIR0UOjKpPFeAXBA6rk5R1ih4MYryB4QJW8o0KMVxA8oEpefu7K4SDGKwgetHoErIwaMV5B8CCSuYqA2p/jONBq06D2IF5v1egb1Jzis9+apN+QbToct7vrsFaG2o/Hn5n5v1E85P/kEh1Nh/IOHStb3nPHeYuHxtlW2qVD/oJ9OvQSANQ+Pc6bHrftIiGmfzy5tptnxcQLg5t8G4Bq0uO4hSM+rLBfeoWKjhrfkEBSYG4QHR5ByCORUZnLqPAdksoUmBtUh0cQ8kjLBJnLqPCtNr+dCGdKcwPAJ/Z7R4dHYnqFc4VvVNEw8c0BJTA3kA6PIOSR3Ja8gwjMcfzoYuyQ/vyc65vSrtLthxYdt3B8jgojjMbdgdR/8VLFzaNifN0LZb4NVGA6zUrMbyjX3M6fS8+5N6Y0Q/v8FPa5XjtA7dH7RPt0YIFpxMEPSz/sfp8Hk/TxapIfPTUT9O+MJqhAFjqwIk0rOvuSd2Dj7UNgbiAdnj88OMhc53ub+fy2JukHf+DeZgC/24/zqEbQ21zZon9DdbP7DbL+5/hhrD6io4aKm7RIHwCozR0yPdpm0g/d3uqvw2+d9cIsrR6ppqfJdAAw0/QLoD2T/bZfvMKe3mEcPKx8Xh89BeZ8dXgEIY/kstpsdXbWtNZP4ArMeenwzH8Zv/VrrQLm1xMl7xQdr9uY4cfa6tP0zWxN0G/QKFUyVxYVJh/3/9YN+FBZFA/pkrW0T8ccFw/d9ELDONu4SafR4KVfTZPexzTpKogqxv+lKhScdQBAgfkvmcm3jzNB30dVH1487ziUvL5t3i6dHRGgE54lfNu8VoBuC8BKp+8oI32Zs6Ozb3ULQg7xcdLgfCDs+kZnSpShz1UkCMKxDE7WQtDPXESdaXOXuSanGK8geNCOgsyFgPSBsMa6YSczcOeH6YEYryB48DR8m+2IzRaAdwH8h9aanOdLjFcQPPAseTkfiDUA79qOrLcBXKdOJsYrCB4Yk70Q9PSBSPhInGCdoEjR3zOP5y3c/fLke4DLKNz9/mS9xMzJWl3g52xtLtBeN81Zeiy5WXPHHycQYeq7eAyVda/sg6BJ/+PcOHDQTMXGtiN3W8RoUw8DQ19DFRM3Kgi643cDZpxXnX5cVTGWNQh0x1Tve0H4QHR8JN7XWt/UWm8AWOCGis7ceAUhj/i2bfvwkega+81CjFcQPIiinHpYCcKzDtMxNRLEeAXBgyE2n70Z+PWhtV7RWhut9bpdPuixD6lxJQh5J4qCzGVU+JS8C2EYKuDEF7NXd/aadbD+2WkyJwjjyhgUvF4hgcn4XJ3RnU1pXAlC7jFj0GHlXcbbsapfZiRTGleCkHuMUZnLqDhNh9WbWdpVg2hc/eDNn5x8X3pxCmpxPk7khL53aOkVACiQukVAVGV0jyru8RdLRVQS26Iq7VwwDCcOxfhYBA13h4tLE069rnBE6z+pQ14fSjEaVdEeI1OTcLJYujbflcxNrm2Y/wkATIWZbL1E/VcP2PMnyftQUc9Ywz40rhzu/X4ze32LNs72I1IeCwBQvEQrPETP0eoK7dlu3aOvN2LRqVaNvoVRaQieQZzI3VH3Dl99tXfyvbjPzJiwQ78kAUBt00Jb7U3Skw8q5SF1L+FJBwAB401npmjjBgAzSRu4KWf/Vwsvsqd3zzUGxus7Y0KXQSb8M3tqXAnCucIQy4g4Tb92eraEXwEnMYk/tb3N61xAsSDkEROpzGVU+GpYbSA1DCQaVsKzxCg7prIQDytB8EGMVxByyhh4aYjxCoIPY9DbfObGW7+6cPK9uTSF+kF8U0pzNfLYIjM2CADt+w/J9OCInhS6eLTgrBcWAhQfbMfHM3ls15ixRwDtEtNvyAiOm4KbbgLlbIvK9Fi0qvJ5VEf0vDBOsH0PTCsxZGcMTNqznxNdL/N55MZ5ufswCIz2wEg4c+MVhFwibV5ByCec19soEOMVBB+k5BWEnOJZ8vYxV9EKrEsxF5l39loegpBHPNwjubmKLL+wRrvMzVUkJa8geKD8hopuAPjEfu/MVXTiPmxL5U+B/lQk+zJerfVK0kd5gGkKe6YLQu7xc9Ig5yoC8AZwUkKvcrbDGq8Nuv8AxxFCTtGvtV7uYdhkepr6XJyFZq3grHN9AuX2DJd9FNp046R9n47jDFSqZbFbBh7H4W9BwLQ8+hALN8wk4u0CM86bGic2RYUouY150FS7n3FeOtxOMfG2Zrser0QR0E6FMTL3yZCxuMe0q/Tj3J4Y3jivenoeVo+tOPuq1vo61e5l27yduUITm7hpCvuZxlAQ8k2kspdsuLmKHiO2tS3YkjgLnw4rrujn0gUh//jF83JzFd1OpM/Btn+zkN5mQfBARdlLFp3mY8ZcRZ1Q2y3bZ7TIDRX59DZzRT+X7nDtSuwbfHnRbTeV5mn9l+Ii7W8LAMFFRmPqCuP3O+HqOy29nJJrmaW1mzhpFoBvi0VFrl3trl664Er3KGYis0KdlgoCgGCfudc7jG/z/tTJ96UfLnSlB/NdE+U5RNO8H3s0eRrf5k0irQeebd4+5irqpLPqqz7G+xEAbb87RX8YhltZ6Vmsf7ufuV7eph+q8mM6qAAACn+htZXYDqvpbuO89z+JYIcLjEbWLP/QtbgOq/LgnWIb9+L7mBaoS1M8oDWuAKCweUDv8IieyD3advXIvv7sO/f8l+gXdbTI57E10603loR6SS51a+KReA4VDRW22myLcG0/+yn6s9IF4dzgU20eNmzJa+vdt1Pb+i36BeF8IsH4cNtryl03TFvPVPhxOzNJV6V6VYuTRLuu5Kk5nHa2FZiYYjXBt8sVI+6cjtdNk5aXjYoK7Uoinrd4+vFNVafb7oUqfZ/N47habdptN74XGM44L9N30JocXv+sRBUJQl6RklcQ8slT9LDqGzFeQfBBjFcQ8om0eQUhp4jxCkJekWqzIOQTKXkBFPfju1A4ipz14iHteqhaQ7iDg+gNo3uM0jQZt702/4o2zBgnN01ou+KOX0Yl1bWNvH4fw8BBkx6vDvbocV5Vio9XxaKzDuA4xpfAMNrVAH+fWtUhxuFIySsI+URKXkHIK1TJO6KYBTFeQfCALHmHp7ZD0lcjIC1BqbVes8utjP1vdfY7fRYFYfxQJnsZFf2EBK4C+Di1fsdGDi3b9TRrWut1uNpXgnB+iIhlRPgI0C0jFpXbQKy5k+TtMAyv2WMF4dzhW/Jqra9bZcibzH5kOuChYRWG4YeJeN0VAGGP3Zb7yaAg5BUf4+1zxoRO7fZNLg/eA1/2wnd7KWWEYfi+zeBiRrVaEPKNX7V5qLLIp+ltXg3D8J30RttJ9cQqcDxG72r1CT+eibOwNFlAObEOJnhazfPdeqrJOGEcMqJlO64zwdIrrpJtMJsSpEuff54O9geA1izt4NCcpn9nu+yOTVycrzjbDDd3d4vvZSnu0/tUrjBOGlfi+7z08ixUOnifEa83F3mB/foC7UjSJJ+nv7DnT+LZMcXKIttJCu5orbtsK42X8Wqt1zpTMWitV+3FOgJ0IeI28jUcz7aQydd/2MlcV01alEw1edVDVW+Q6WZ3j0xvP+4WVvs6/Pbke+ECc/7LfA9GgxGxO1qg/6bWRPdD+cf7sWBcxLzjCg3+Saxs079j4vt9Mj3YeOisOyJ+AFCgjTd6mX9UD56fINPrM9k34qXL7OkdnqKTRre0ZgYDC9DZavAtrfW61jqpl5kUoPup3X9dBOiEc4mf6Dopi9wpdfvNwsACdPbkXUKZIkAnPEt4lrycbPKy1noZxwa+wM3zJTMmCIIHT2nGhNuJWRJoFXqIe6Qg+PGUZkxI7MPWXsV4BcEDFZ19TKAYryB4IOqRAIJvvj/5rubbCL5JDCEEzBhHHwHahnlDmgYTTJ8OlFfK2aa4IY5CH2LhTOA8F0TemHLz2KoqZ5th/uU2Pw0QuO6R4iE9zludi8dpVW0Sas4dt21/d588XjX6GBZkOpH6ER3oF4nnFYS8IiWvIOQTKXkFIadIh5Ug5BTpsBKEnKJot/uRIMYrCD6MQcnrq2FFalT1qxYgCHlFRSZzGRVsyWv9MD/AcXhfhzUbNfSzHvufqAVorZc55+rkOKtptpz1aG+z1yEnqCIdiwsAqsTE8zbokL6g4k4qrUold1uVnnQ6muArN9ykz80p+vhmKtS1Nelui5jbpNp9CJoX6X0KTfoipe04rjmamUB70Y1zVt8/II9X+0dMDoGgUaN3GOJA7zi0eX00rABao2qoagGCMI74BCYMG9+oIkqjilULEIS8Mw7VZi/jFY0q4ZnHLxh/qAzc29yHRhWpFpDm6t9dPPm+tOyGMJojWh9K9eE3DG6fFuMzmzp+6ZqrQxDM03lsLTANVvDaS/VZZgKtSTf9hZrbDud8m/sZ9ijQWUS1wMjkJHzAly5Pdu9QpnVo1GKX/kMXjcv0vT6iNM8GNDrVxwRyTxufoaKeGlUJNYCeagFZ3PvNg8z1aI/WlxpFh5XqMYvgvbtxMEVwiX7ym8/zf/IB49R/0GACE3oY9xdbCQ0rtsOKTgeA0i6dXntAG+/0hvtffv1HV7tMff49KNQL9PUB4JC5T/v17Btx9WJmUm/O3nYH17AiNKqSGla91AIE4dwwDm3egTWs7DZSDUA0rITzzjgMFZ25h5W6dCH+Pj8LdSkx2fYm3dCKtrbZ80cHB/QOzMTWwaw7iKqqFaha3GYzk5zmMl+1r08PXi120ufcJ6lVM2gk6nXtCv2k9TO8EZXpPAYtOo+l3ViWtTFTxtEFV6a1duE58niztUOmA0DxgJZ9ClrDe9wlMEEQckpeO6wEQfC0XdtXtAVgpTNxQSq943J8rdeMJElE+lUQPPDpsOImGutz+twTxHgFwQdjspdsONfhfqbPPUGqzYLggWebl3QdTo3SrODYZyITMV5B8OEp9ldR0+cmEeMVBA9U5BU+1K/rcM/pc9OcufGaWjxOaqplZz04rJPHqgof48mhKnQ8rppyY0RVtepsa84x00rO87e4vkCPkdYX6Nd8a971z26X22hW423BBO2/HfURz1uv0OPVqkX7kJd34/TmTICjRXf/6nPMPMe//jOTQ6Cwd4lML9YZB+0B8Kw2cxON9Zw+N+tk0mElCD54dFhxE40R0+f25MxLXkHIJZ5OGpRrcdb0uVkMrGGltV7RWhv7dljXWn/QY39S40oQ8o4yJnMZFT4aVgthGCqbtoJ43CpJpsaVIJwL/DqshsrAGlapBrQOwzCtbwXQGleCkH/8nDSGineHlS2Rf5mRTGlcCULuUW2TuYyK0/Q2v9np3k4jGlfCuacdZS8j4jS9zSu9NvahceVw9Ufx+N7S8662UXCJzl60R8fSAgCazOSzJXr8Mphx9YWXXnLXW4uMbtIFfmzxaIF+h9Zr9Ns8Krs6Ni8V3fuiCrTOjWFimgEAVfq/KM/Sv2HycvxQX57vHluvNek8mMIVMh0A1LUZMv3wEjUm/4g9v5sh6j/p434OAS/j1Vp3GWRioLmnxlUW975IDGdFEe79IQ6wDx7SQ13RDiOsZM9JoWbpP9xcco0vqlVx79vYeeSAUWbbmeXfxHsF2jjbNVpna3rWFRyYMLP4RsX3capKO7scNfnHYFv1EI1zdqBfpFNHsXE3GgG+PHIdbBYf0b+xss6IKgBQT2gdrP0fX8hMu/Iae3oX8rka4izeBANrWCVId1QlNax6aVwJwvkhMtnLiPDVsNpAahhINKyEZ4qIaoqMxvdJPKwEwQfRsBKEnDIGThpivILgQ/vsZ9cW4xUEH0boSZWFGK8g+DBCZ4wsztx4g+39+PtB1VmPdum5irh5hgBATdDB8pimJ2RuPOeOb7ZmKs62/Yv0mN7BJf4N3b5Mj8NeXKLF5V+cdh3dXqgDqDw8Wa8WaEeVnQZzjwAYQzsebNbo8e5WLR6VbFeBVuq2N2boR7E6R4/HA4DZpoXZSzu0KPsgGCPGKwj5REpeQcgp0tssCPnESG+zIOQUqTYLQk6RDitByCdSbRaEnGLEtxn42xvJ8cFdzL2WXKdFtEdDWixkC0uLPXfsDT8nNPDbU6YjPX7ZwA+7tp2On4AWbs8W/+/BE+BH6W2vMMe80s+jyqmm9hH/3R/3Xv/nl69S6cO6EIUyY+DmJQjC4MiMCYKQU8R4BSGniPGeE3rMsn59nOR3e+RPZtU4JWNhvOP2oPVinB82K6/7cWJ9BTgRyN9KG86oSefPsqa1Xke3FprQJ2duvOP2oBGM7cOWntUCwA3E3eQbAM5UO7tH/gCZVePUnPlQEY4ftE/s986DNo6Kk29bMb48MAfgSWJ9kMGtUbFsS+SVzny0wmCcecmLfDxogEzhMlRkVo3TMw7Gmwty9rBtAViw3+cwkAfF00drvZbQAWdn1RB6Mw7GO9YPGpDLh+0jxHlcBjBu7coQcZ6u2XVhQMbBeMf9QQPG/GFLz2rRmaXC1hC2znrWioz8yawap2Qs3CPt8MsGgOVxnW2hM4EajvMoHSzCmTMWxisIwuCMQ7VZEAQPxHgFIaeI8QpCThHjFYScIsYrCDlFjFcQcooYryDkFDFeQcgp/w+S+wRUgkp9IAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "pred_labels = observed_pred.mean.view(n, n)\n",
    "\n",
    "# Calc abosolute error\n",
    "test_y_actual = torch.sin(((test_x[:, 0] + test_x[:, 1]) * (2 * math.pi))).view(n, n)\n",
    "delta_y = torch.abs(pred_labels - test_y_actual).detach().numpy()\n",
    "\n",
    "# Define a plotting function\n",
    "def ax_plot(f, ax, y_labels, title):\n",
    "    if smoke_test: return  # this is for running the notebook in our testing framework\n",
    "    im = ax.imshow(y_labels)\n",
    "    ax.set_title(title)\n",
    "    f.colorbar(im)\n",
    "\n",
    "# Plot our predictive means\n",
    "f, observed_ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax, pred_labels, 'Predicted Values (Likelihood)')\n",
    "\n",
    "# Plot the true values\n",
    "f, observed_ax2 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax2, test_y_actual, 'Actual Values (Likelihood)')\n",
    "\n",
    "# Plot the absolute errors\n",
    "f, observed_ax3 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax3, delta_y, 'Absolute Error Surface')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}