ML_course/assignment 6/iml_assgnment6_unsolved.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Solution for Assignment 6 of the course “Introduction to Machine Learning” at the University of Leoben.\n",
    "Author: Fotios Lygerakis\n",
    "Semester: SS 2022/2023"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Implement the Gaussian Processes algorithm as a class in Python. The class should be able to perform the following tasks:\n",
    "1. Fit the model to a given training data set (X, y) with X ∈ Rn×d and y ∈ Rn.\n",
    "2. Predict the output for a given test data set X∗ ∈ Rm×d.\n",
    "3. Compute the log marginal likelihood of the training data set.\n",
    "4. Compute the gradients of the log marginal likelihood with respect to the hyperparameters.\n",
    "5. Compute the predictive mean and variance for a given test data set X∗ ∈ Rm×d.\n",
    "6. Compute the gradients of the predictive mean and variance with respect to the hyperparameters.\n",
    "7. Use the gradients of the log marginal likelihood and the gradients of the predictive mean and variance to optimize the hyperparameters of the kernel function.\n",
    "\n",
    "Implement the following kernel functions:\n",
    "1. linear kernel: k(x, x′) = xT x′\n",
    "2. polynomial kernel: k(x, x′) = (xT x′ + 1)d\n",
    "3. Periodic kernel: k(x, x′) = exp(−2sin2(π||x − x′||/p)/σ2)\n",
    "4. Gaussian kernel: k(x, x′) = exp(−||x − x′||2/2σ2)\n",
    "\n",
    "Test your implementation on 1D, 2D and 3D data sets. Compare the results with the results of the scikit-learn implementation of Gaussian Processes."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Import the necessary libraries"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Datasets"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "outputs": [],
   "source": [
    "step_size = 0.5\n",
    "#1D data set\n",
    "X_1D = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "y_1D = np.sin(X_1D) + np.random.normal(0, 0.1, size=(len(X_1D), 1))\n",
    "X_1D_star = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "\n",
    "#2D data set\n",
    "X_2D_1 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_2D_2 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_2D = np.hstack((X_2D_1, X_2D_2))\n",
    "y_2D = np.sin(X_2D_1) + np.exp(X_2D_2) + np.random.normal(0, 0.1, size=(len(X_2D), 1))\n",
    "X_2D_star_1 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_2D_star_2 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_2D_star = np.hstack((X_2D_star_1, X_2D_star_2))\n",
    "\n",
    "#3D data set\n",
    "X_3D_1 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_3D_2 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_3D_3 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_3D = np.hstack((X_3D_1, X_3D_2, X_3D_3))\n",
    "y_3D = (\n",
    "    np.sin(X_3D_1)\n",
    "    + np.cos(X_3D_2)\n",
    "    + X_3D_3 ** 2\n",
    "    + np.random.normal(0, 0.1, size=(len(X_3D), 1))\n",
    ")\n",
    "X_3D_star_1 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_3D_star_2 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_3D_star_3 = np.arange(-5, 5, step_size).reshape(-1, 1)\n",
    "X_3D_star = np.hstack((X_3D_star_1, X_3D_star_2, X_3D_star_3))\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Implement the Gaussian Processes algorithm as a class in Python"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "outputs": [],
   "source": [
    "class GaussianProcess:\n",
    "    def __init__(self, kernel, sigma_n=1e-8):\n",
    "        pass\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        # Compute the kernel matrix. Add a small value to the diagonal for numerical stability.\n",
    "        # Compute the inverse of the kernel matrix.\n",
    "        # Compute the coefficients alpha.\n",
    "        pass\n",
    "\n",
    "    def predict(self, X_star):\n",
    "        # Compute the kernel matrix between the training data and the test data.\n",
    "        # Compute the predictive mean.\n",
    "        # Compute the predictive variance.\n",
    "        y_star, var = None, None\n",
    "        return y_star, var\n",
    "\n",
    "    def log_marginal_likelihood(self):\n",
    "        # Compute the log marginal likelihood.\n",
    "        lml=None\n",
    "        return lml\n",
    "\n",
    "    def log_marginal_likelihood_gradient(self):\n",
    "        # Compute the gradient of the log marginal likelihood with respect to the kernel hyperparameters.\n",
    "        lml_grad=None\n",
    "        return lml_grad\n",
    "\n",
    "    def predict_gradient(self, X_star):\n",
    "        # K_star is the gradient of the kernel function\n",
    "        # y_star is the gradient of the mean function\n",
    "        # v is the gradient of the covariance function\n",
    "        # var is the gradient of the variance function\n",
    "        # K_star_grad is the gradient of the kernel function\n",
    "        # y_star_grad is the gradient of the mean function\n",
    "        # v_grad is the gradient of the covariance function\n",
    "        # var_grad is the gradient of the variance function\n",
    "        y_star, var, K_star_grad, v_grad = None, None, None, None\n",
    "        y_star_grad, var_grad = None, None\n",
    "        return y_star, var, y_star_grad, var_grad, K_star_grad, v_grad\n",
    "\n",
    "    def plot(self, X_star, y_star, var):\n",
    "        plt.figure(figsize=(10, 10))\n",
    "        plt.plot(self.X, self.y, \"r.\", markersize=10, label=\"Observations\")\n",
    "        plt.plot(X_star, y_star, \"b-\", label=\"Prediction\")\n",
    "        plt.fill_between(\n",
    "            X_star[:, 0],\n",
    "            y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
    "            y_star[:, 0] + 1.96 * np.sqrt(np.diag(var)),\n",
    "            alpha=0.2,\n",
    "            color=\"k\",\n",
    "            label=\"95% confidence interval\",\n",
    "        )\n",
    "        plt.xlabel(\"$x$\")\n",
    "        plt.ylabel(\"$y$\")\n",
    "        plt.title(\n",
    "            f\"{self.kernel.__class__.__name__} \\n Log-Marginal-Likelihood: {self.log_marginal_likelihood()[0][0]:.3e}\"\n",
    "        )\n",
    "        plt.legend(loc=\"upper left\")\n",
    "        plt.show()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "outputs": [],
   "source": [
    "#Implement the linear kernel function\n",
    "class LinearKernel:\n",
    "    def __init__(self, theta=1.0):\n",
    "        self.theta = theta\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def gradient(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def set_params(self, params):\n",
    "        pass\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "outputs": [],
   "source": [
    "#Implement the polynomial kernel function\n",
    "class PolynomialKernel:\n",
    "    def __init__(self, theta=1.0, d=3):\n",
    "        self.theta = theta\n",
    "        self.d = d\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def gradient(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def set_params(self, params):\n",
    "        pass"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "outputs": [],
   "source": [
    "#Implement the periodic kernel function\n",
    "class PeriodicKernel:\n",
    "    def __init__(self, theta=1.0):\n",
    "        self.theta = theta\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def gradient(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def set_params(self, params):\n",
    "        pass\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "outputs": [],
   "source": [
    "#Implement the radial basis function kernel function\n",
    "class RBFKernel:\n",
    "    def __init__(self, theta=1.0):\n",
    "        self.theta = theta\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def gradient(self, X1, X2):\n",
    "        pass\n",
    "\n",
    "    def set_params(self, params):\n",
    "        pass\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Test your implementation. Compare the results with the results of the scikit-learn implementation of Gaussian Processes."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Initialize the kernels"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "outputs": [],
   "source": [
    "#Linear kernel\n",
    "kernel_linear = None\n",
    "#Polynomial kernel\n",
    "kernel_poly = None\n",
    "#Periodic kernel\n",
    "kernel_periodic = None\n",
    "#Radial basis function kernel\n",
    "kernel_rbf = None"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Initialize the Gaussian Process with the kernels and fit the data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "outputs": [],
   "source": [
    "gp_linear = None\n",
    "gp_poly = None\n",
    "gp_periodic = None\n",
    "gp_rbf = None"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Fit the 1D data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "outputs": [],
   "source": [
    "gp_linear.fit(X_1D, y_1D)\n",
    "gp_poly.fit(X_1D, y_1D)\n",
    "gp_periodic.fit(X_1D, y_1D)\n",
    "gp_rbf.fit(X_1D, y_1D)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict the values of the test data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "outputs": [],
   "source": [
    "y_star_linear, var_linear = gp_linear.predict(X_1D_star)\n",
    "y_star_poly, var_poly = gp_poly.predict(X_1D_star)\n",
    "y_star_periodic, var_periodic = gp_periodic.predict(X_1D_star)\n",
    "y_star_rbf, var_rbf = gp_rbf.predict(X_1D_star)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Plot the results"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_37332/3048407143.py:67: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
      "/tmp/ipykernel_37332/3048407143.py:68: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] + 1.96 * np.sqrt(np.diag(var)),\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACLsElEQVR4nOzdd3hUVf7H8c+kTHohJBBKIEACofcSWihxUdS1sYiiEAno2l1xV9mfuq6oKPbVta2UXcFeWNuKEAkd6UUETIAISE8gvc7c3x8jA0MSDHCTSXm/nmeece49997vTAbJh3PuORbDMAwBAAAAAEzj4e4CAAAAAKC+IWgBAAAAgMkIWgAAAABgMoIWAAAAAJiMoAUAAAAAJiNoAQAAAIDJCFoAAAAAYDKCFgAAAACYjKAFAAAAACYjaAEA3CYjI0MWi0Vz5851dyn1Gp8zANQ8ghYAoFrMnTtXFotF69evd3cpphg2bJi6dOlSbntKSor8/f3Vq1cvZWVluaEyAEBt5OXuAgAADVfr1q1VWFgob29vd5dyQb777jtdeeWV6tChgxYvXqywsDB3lwQAqCXo0QIAuI3FYpGvr688PT3dXUqlCgoKKty+dOlSXXnllWrfvr1pISs/P/+izwEAqB0IWgAAt6no3qGkpCQFBgbql19+0dVXX63AwEBFRETogQcekM1mcznebrfrpZdeUufOneXr66umTZvqtttu04kTJ1za/fe//9Xll1+u5s2by8fHR+3atdP06dPLne/U8MANGzZo6NCh8vf311//+tdydS9fvlyXX365YmJitHjxYjVu3Nhl///+9z8NGTJEAQEBCgoK0uWXX67t27e7tDn1Pnfv3q3Ro0crKChI48ePl+QIoHfddZcWLFigLl26yMfHR507d9Y333xTrpZffvlFkyZNUtOmTZ3tZs+e/dsfPgCgWjF0EABQ69hsNo0aNUr9+/fXc889p8WLF+v5559Xu3btdPvttzvb3XbbbZo7d65uueUW3XPPPdq7d69effVVbdq0SStXrnQOSZw7d64CAwN1//33KzAwUN99950effRR5eTk6Nlnn3W5dmZmpi677DKNGzdON910k5o2beqyf+XKlRo9erTatGmjlJQUhYeHu+x/5513NHHiRI0aNUrPPPOMCgoK9Prrr2vw4MHatGmToqOjnW3Lyso0atQoDR48WM8995z8/f2d+1asWKFPP/1Ud9xxh4KCgvSPf/xD1113nfbt2+cMdkeOHNGAAQOcwSwiIkL/+9//lJycrJycHN13331m/DgAABfCAACgGsyZM8eQZKxbt67SNnv37jUkGXPmzHFumzhxoiHJePzxx13a9uzZ0+jdu7fz9fLlyw1Jxvz5813affPNN+W2FxQUlLv2bbfdZvj7+xtFRUXObQkJCYYk44033ijXPiEhwQgLCzOCgoKMzp07G0ePHi3XJjc31wgNDTWmTJnisv3w4cNGSEiIy/ZT7/Ohhx4qdx5JhtVqNdLT053btmzZYkgyXnnlFee25ORko1mzZsbx48ddjh83bpwREhLifN8Vfc4AgOrF0EEAQK30xz/+0eX1kCFDtGfPHufrjz76SCEhIbrkkkt0/Phx56N3794KDAzUkiVLnG39/Pyc/52bm6vjx49ryJAhKigo0M6dO12u4+Pjo1tuuaXCmvLz85Wbm6umTZsqODi43P5Fixbp5MmTuuGGG1xq8vT0VP/+/V1qOuXMHrozJSYmql27ds7X3bp1U3BwsPMzMAxDn3zyia688koZhuFyvVGjRik7O1sbN26s8NwAgOrH0EEAQK3j6+uriIgIl22NGjVyufcqLS1N2dnZatKkSYXnOHr0qPO/t2/frocffljfffedcnJyXNplZ2e7vG7RooWsVmuF54yJidGECRP04IMP6oYbbtBHH33kMpFHWlqaJGnEiBEVHn92OPPy8lLLli0rbNuqVaty2878DI4dO6aTJ0/qrbfe0ltvvVXhOc78DAAANYugBQCodaoyC6HdbleTJk00f/78CvefCmonT55UQkKCgoOD9fjjj6tdu3by9fXVxo0b9eCDD8put7scd2bvV0X+8pe/KDMzUzNnztSUKVM0a9YsWSwWZ02S4z6tyMjIcsd6ebn+tevj4yMPj4oHl1T2GRiG4XKtm266SRMnTqywbbdu3c75XgAA1YegBQCok9q1a6fFixdr0KBB5wxHqampyszM1KeffqqhQ4c6t+/du/eCr/3MM88oKytLb7/9tho1aqTnn3/eWZMkNWnSRImJiRd8/qqIiIhQUFCQbDZbtV8LAHD+uEcLAFAnjR07VjabTdOnTy+3r6ysTCdPnpR0umfoVE+QJJWUlOi11167qOu/+eabGjNmjF544QU98cQTkqRRo0YpODhYTz31lEpLS8sdc+zYsYu65pk8PT113XXX6ZNPPtEPP/xQrdcCAJw/erQAANVq9uzZFa7/dO+9917UeRMSEnTbbbdpxowZ2rx5s373u9/J29tbaWlp+uijj/Tyyy9rzJgxGjhwoBo1aqSJEyfqnnvukcVi0TvvvOMSvC6Eh4eH5s+fr+zsbD3yyCMKCwvTHXfcoddff10333yzevXqpXHjxikiIkL79u3TV199pUGDBunVV1+9qOue6emnn9aSJUvUv39/TZkyRZ06dVJWVpY2btyoxYsXKysry7RrAQDOD0ELAFCtXn/99Qq3JyUlXfS533jjDfXu3Vtvvvmm/vrXv8rLy0vR0dG66aabNGjQIElS48aN9eWXX2rq1Kl6+OGH1ahRI910000aOXKkRo0adVHXt1qt+uyzz5SYmKi7775boaGhuvHGG9W8eXM9/fTTevbZZ1VcXKwWLVpoyJAhlc5meKGaNm2qtWvX6vHHH9enn36q1157TY0bN1bnzp31zDPPmHotAMD5sRgX+096AAAAAAAX3KMFAAAAACYjaAEAAACAyQhaAAAAAGAyghYAAAAAmIygBQAAAAAmI2gBAAAAgMkIWgAAFxaLRY899li1XmPYsGEaNmyYaec7u+bHHntMFotFx48fr5bzz507VxaLRRkZGc5t0dHRuuKKK0y5nhkyMjJksVg0d+5cd5cCAA0SQQtAvXXqF83nnnvO3aVIcizQa7FYFBwcrMLCwnL709LSZLFYalXN7paamiqLxaKPP/7Y3aXgV6dC7NkPX1/fKh3/1FNPacCAAYqIiJCvr69iY2N133336dixY+c8bv78+bJYLAoMDKxwv91u1+uvv64ePXrIz89PjRs31ogRI7Rly5bzfo8Xa9asWerYsaPz/b3yyisVtlu8eLGGDx+u8PBwhYaGql+/fnrnnXdquFoA1cXL3QUAQEPi5eWlgoICffHFFxo7dqzLvvnz58vX11dFRUVuqs6hsLBQXl5166+Hmq755ptv1rhx4+Tj41Nj16xtXn/9dZfQ4+npWaXjNmzYoB49emjcuHEKCgrSjh079K9//UtfffWVNm/erICAgHLH5OXl6S9/+UuF+06ZNGmS5s+frwkTJuiuu+5Sfn6+Nm3apKNHj57/m7sIb775pv74xz/quuuu0/3336/ly5frnnvuUUFBgR588EFnu88//1xXX3214uPjneH1ww8/1IQJE3T8+HH96U9/qtG6AZivbv1NCgB1nI+PjwYNGqT33nuvXNB69913dfnll+uTTz4x7XpFRUWyWq3y8Kj6AIaq9kzUJjVds6enZ5WDRX01ZswYhYeHn/dxFX2/4+PjNWbMGH3xxRcaN25cuf1PPPGEgoKCNHz4cC1YsKDc/g8//FD//ve/9emnn+qaa64575qqYu7cubrllltkGEalbQoLC/V///d/uvzyy529sFOmTJHdbtf06dN16623qlGjRpKkV199Vc2aNdN3333nDOy33Xab4uLiNHfuXIIWUA8wdBBAg3f06FElJyeradOm8vX1Vffu3fXvf/+7XLvMzEzdfPPNCg4OVmhoqCZOnKgtW7ac930wN954o/73v//p5MmTzm3r1q1TWlqabrzxxnLts7Ky9MADD6hr164KDAxUcHCwLrvssnJDok4Ns3v//ff
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACMyUlEQVR4nOzdd3gU9drG8XvTE9IogQQIBEiA0JsgzaggzQIiIp1AKBYsB2x4LCgH8GADRawUD6IgKoiIICKgIIoiTUVMIBFQSmgJIaTP+8e8WVhSCDDJJuH7ua65lp2ZnXl2WTR3fr95xmYYhiEAAAAAgGVcnF0AAAAAAJQ3BC0AAAAAsBhBCwAAAAAsRtACAAAAAIsRtAAAAADAYgQtAAAAALAYQQsAAAAALEbQAgAAAACLEbQAAAAAwGIELQBAsbj++ut1/fXXO7sMS8yfP182m00JCQmX/Nro6GiFhYVZXpOzXMlnAQBXE4IWAEDSuR+gcxcvLy/Vr19f48aN05EjR5xdXrlis9k0bty4POunTp0qm82mkSNHKicnxwmVAQCs4ubsAgAApctzzz2nOnXqKC0tTRs3btQbb7yhlStX6tdff5WPj4+zy3OKoUOHasCAAfL09Cy2czz//PP697//reHDh+vdd9+Viwu/CwWAsoygBQBw0LNnT7Vp00aSNGrUKFWuXFkvv/yyPvvsMw0cONDJ1TmHq6urXF1di+34L7zwgiZOnKhhw4Zp7ty5VxyyDMNQWlqavL29LaoQAHCp+HUZAKBQN954oyQpPj5ekpSVlaXJkyerXr168vT0VFhYmJ544gmlp6cXeIyUlBRVqFBBDz74YJ5tBw8elKurq6ZNmybp3BTGTZs2afz48QoKClKFChV0++23KzExMc/rZ8+ercaNG8vT01PVq1fXfffdp1OnTjnsc/3116tJkybauXOnoqKi5OPjo/DwcH388ceSpA0bNqhdu3by9vZWgwYN9PXXXzu8Pr/rkj777DPdfPPNql69ujw9PVWvXj1NnjxZ2dnZF/9Qz/Pyyy/r0Ucf1ZAhQzRv3jyHkJWTk6MZM2aocePG8vLyUrVq1TR27FidPHnS4RhhYWG65ZZbtHr1arVp00be3t566623tH79etlsNn300UeaMmWKatasKS8vL3Xp0kVxcXF5avnxxx/Vo0cPBQQEyMfHR1FRUdq0adMlvR8AgImgBQAo1N69eyVJlStXlmSOcj399NNq1aqVXnnlFUVFRWnatGkaMGBAgcfw9fXV7bffrsWLF+cJIh9++KEMw9DgwYMd1t9///3asWOHnnnmGd1zzz36/PPP81zXNGnSJN13332qXr26XnrpJd1xxx1666231K1bN2VmZjrse/LkSd1yyy1q166dpk+fLk9PTw0YMECLFy/WgAED1KtXLz3//PM6c+aM+vXrp9OnTxf6ucyfP1++vr4aP368Zs6cqdatW+vpp5/W448/XvgHep6ZM2dqwoQJGjRokObPn59nJGvs2LF65JFH1LFjR82cOVMjRozQwoUL1b179zzvb8+ePRo4cKBuuukmzZw5Uy1atLBve/7557V06VI9/PDDmjhxon744Yc8n/c333yj6667TsnJyXrmmWc0depUnTp1SjfeeKO2bNlS5PcEAPh/BgAAhmHMmzfPkGR8/fXXRmJionHgwAFj0aJFRuXKlQ1vb2/j4MGDxvbt2w1JxqhRoxxe+/DDDxuSjG+++ca+LioqyoiKirI/X716tSHJ+PLLLx1e26xZM4f9cuvo2rWrkZOTY1//r3/9y3B1dTVOnTplGIZhHD161PDw8DC6detmZGdn2/ebNWuWIcmYO3euQy2SjA8++MC+7o8//jAkGS4uLsYPP/yQp8558+blqSk+Pt6+LjU1Nc9nOHbsWMPHx8dIS0uzrxs+fLhRu3Zth/0kGbVr1zYkGQMHDjSysrLyHOu7774zJBkLFy50WL9q1ao863OPtWrVKod9161bZ0gyIiMjjfT0dPv6mTNnGpKMXbt2GYZhGDk5OUZERITRvXt3h888NTXVqFOnjnHTTTcV+lkAAPJiRAsA4KBr164KCgpSaGioBgwYIF9fXy1dulQ1atTQypUrJUnjx493eM2ECRMkSV988UWhx61evboWLlxoX/frr79q586dGjJkSJ79x4wZI5vNZn/euXNnZWdn66+//pIkff3118rIyNBDDz3kMBI0evRo+fv756nF19fXYdStQYMGCgwMVGRkpNq1a2dfn/vnffv2FfheJDlc/3T69GkdO3ZMnTt3Vmpqqv74449CXyvJ3smxTp06+V7/tWTJEgUEBOimm27SsWPH7Evr1q3l6+urdevWOexfp04dde/ePd9zjRgxQh4eHvbnnTt3dniP27dvV2xsrAYNGqTjx4/bz3XmzBl16dJF3377LV0QAeAS0QwDAODg9ddfV/369eXm5qZq1aqpQYMG9iDz119/ycXFReHh4Q6vCQ4OVmBgoD0E5cfFxUWDBw/WG2+8odTUVPn4+GjhwoXy8vLSnXfemWf/WrVqOTyvWLGiJNmvT8o9V4MGDRz28/DwUN26dfPUUrNmTYfgJkkBAQEKDQ3Ns+788xTkt99+05NPPqlvvvlGycnJDtuSkpIKfa0kDR8+XP/884+mTp2qKlWq6F//+pfD9tjYWCUlJalq1ar5vv7o0aMOz+vUqVPguS72WcbGxtprKkhSUpL9dQCAiyNoAQActG3b1t51sCAXBpaiGjZsmF544QUtW7ZMAwcO1AcffKBbbrnFHm7OV1CXP8MwLuvcBR3vcs5z6tQpRUVFyd/fX88995zq1asnLy8v/fLLL3rssceKNPrj5uamjz76SD169NCECRMUGBioESNG2Lfn5OSoatWqDiOA5wsKCnJ4XliHwYu9x9x6X3jhBYdru87n6+tb4PEBAHkRtAAARVa7dm3l5OQoNjZWkZGR9vVHjhzRqVOnVLt27UJf36RJE7Vs2VILFy5UzZo1tX//fr322muXXYtkNoGoW7eufX1GRobi4+PVtWvXyzpuUaxfv17Hjx/Xp59+quuuu86+PrczY1F5eXlp+fLluuGGGzR69GgFBgbq9ttvlyTVq1dPX3/9tTp27Fjsbdrr1asnSfL39y/Wzw0AriZcowUAKLJevXpJkmbMmOGw/uWXX5Yk3XzzzRc9xtChQ/XVV19pxowZqly5snr27HlZtXTt2lUeHh569dVXHUaf5syZo6SkpCLVcrlyR4jOP29GRoZmz559ycfy9/fXqlWrFB4eroEDB2rt2rWSpP79+ys7O1uTJ0/O85qsrKw8LeyvROvWrVWvXj29+OKLSklJybM9v7b6AIDCMaIFACiy5s2ba/jw4Xr77bft0+e2bNmi9957T3369NENN9xw0WMMGjRIjz76qJYuXap77rlH7u7ul1VLUFCQJk6cqGeffVY9evTQbbfdpj179mj27Nm65ppr8m2wYZUOHTqoYsWKGj58uB544AHZbDYtWLDgsqc1BgUFac2aNerYsaP69OmjtWvXKioqSmPHjtW0adO0fft2devWTe7u7oqNjdWSJUs0c+ZM9evXz5L34+LionfffVc9e/ZU48aNNWLECNWoUUN///231q1bJ39/f33++eeWnAsArhYELQDAJXn33XdVt25dzZ8/X0uXLlVwcLAmTpyoZ555pkivr1atmrp166aVK1dq6NChV1TLpEmTFBQUpFmzZulf//qXKlWqpDFjxmjq1KmXHeCKonLlylqxYoUmTJigJ598UhUrVtSQIUPUpUuXAjv/XUxoaKi++uorde7cWT179tS3336rN998U61bt9Zbb72lJ554Qm5ubgoLC9OQIUPUsWNHS9/T9ddfr82bN2vy5MmaNWuWUlJSFBwcrHbt2mns2LGWngsArgY243J//QYAwGW6/fbbtWvXLsXFxTm7FAAAigXXaAEAStShQ4f0xRdfXPFoFgAApRlTBwEAJSI+Pl6bNm3Su+++K3d3d6ajAQDKNUa0AAAlYsOGDRo6dKji4+P13nvvKTg42NklAQBQbLhGCwAAAAAsxogWAAAAAFi
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACUJUlEQVR4nOzdd3RU1d7G8WeSkFBCEnoTCJLQBBVUEERB5YpiARURRSEmgAoIXrGAildBRbHhFcVGgldQBCli5SqKAmIBuyI3IBFUmvQeyOz3j/1mMkMKKWda8v2slTVk5sw5v5mEzH7OLsdljDECAAAAADgmItgFAAAAAEB5Q9ACAAAAAIcRtAAAAADAYQQtAAAAAHAYQQsAAAAAHEbQAgAAAACHEbQAAAAAwGEELQAAAABwGEELAAAAABxG0AIAhIXu3bure/fuju0vKytLLpdL06dP99x3//33y+VyOXaMcFDQ+wAAKDuCFgCg1KZPny6Xy+X5qly5slq0aKERI0Zoy5YtwS4vYLp37662bdvmu3/x4sWqWrWqOnTooB07dgShMgBAsEQFuwAAQPgbP368mjVrpkOHDmnZsmWaOnWq3nvvPf3000+qWrWqI8f473//68h+inLvvfdqzJgxjuzr448/1qWXXqqWLVvqo48+Us2aNR3ZLwAgPBC0AABldtFFF+n000+XJA0ePFi1atXSk08+qbfeekvXXHNNmfZ94MABVa1aVdHR0U6UWqSoqChFRZX9o/HTTz/VpZdeqhYtWjgWsvbv369q1aqVeT8AgMBg6CAAwHHnnXeeJGn9+vWe+2bMmKHTTjtNVapUUc2aNdW/f39t3LjR53m5Q/BWrVqlc845R1WrVtXdd9/teezYOVpbt25VWlqa6tWrp8qVK+uUU07RK6+8kq+eXbt2KSUlRfHx8UpISNCgQYO0a9eufNsVNkdrxowZ6tixo6pWraoaNWronHPOKbSHbenSpbr44ouVlJSkjz76SLVq1fJ5/P3339fZZ5+tatWqqXr16rr44ov1888/+2yTkpKi2NhYrVu3Tr169VL16tU1YMAASZLL5dKIESO0YMECtW3bVjExMTrppJP0wQcf5Kvlzz//VGpqqurVq+fZLj09vcC6AQDOokcLAOC4devWSZInZDz00EMaN26c+vXrp8GDB2vbtm165plndM455+jbb79VQkKC57nbt2/XRRddpP79++u6665TvXr1CjzGwYMH1b17d61du1YjRoxQs2bNNGfOHKWkpGjXrl0aNWqUJMkYo969e2vZsmW66aab1Lp1a82fP1+DBg0q1mt54IEHdP/996tLly4aP368oqOj9eWXX+rjjz/WBRdc4LPt8uXL1atXLzVr1kyLFy9W7dq1fR5/9dVXNWjQIPXs2VOPPvqoDhw4oKlTp6pr16769ttvlZiY6Nn26NGj6tmzp7p27arHH3/cZwjmsmXLNG/ePA0bNkzVq1fXv//9b1155ZXasGGD5z3fsmWLzjzzTE8wq1Onjt5//32lpaVpz549uvXWW4v1+gEApWQAACiljIwMI8l89NFHZtu2bWbjxo1m1qxZplatWqZKlSrmjz/+MFlZWSYyMtI89NBDPs/98ccfTVRUlM/93bp1M5LM888/n+9Y3bp1M926dfN8P3nyZCPJzJgxw3Nfdna26dy5s4mNjTV79uwxxhizYMECI8lMmjTJs93Ro0fN2WefbSSZjIwMz/3/+te/jPdHY2ZmpomIiDCXX365ycnJ8anH7Xb71FazZk1TvXp1c9JJJ5mtW7fmq3/v3r0mISHBDBkyxOf+zZs3m/j4eJ/7Bw0aZCSZMWPG5NuPJBMdHW3Wrl3rue/77783kswzzzzjuS8tLc00aNDA/P333z7P79+/v4mPjzcHDhwwxhizfv36fO8DAKDsGDoIACizHj16qE6dOmrcuLH69++v2NhYzZ8/X40aNdK8efPkdrvVr18//f33356v+vXrKzk5WZ988onPvmJiYnTDDTcc95jvvfee6tev7zMHrFKlSho5cqT27dunTz/91LNdVFSUbr75Zs92kZGRuuWWW457jAULFsjtduu+++5TRITvR+axQwz379+vvXv3ql69eoqLi8u3rw8//FC7du3SNddc4/M+REZGqlOnTvneB0k+NXvr0aOHmjdv7vn+5JNPVlxcnH777TdJthdv7ty5uvTSS2WM8Tlez549tXv3bn3zzTfHff0AgNJj6CAAoMyeffZZtWjRQlFRUapXr55atmzpCSaZmZkyxig5ObnA51aqVMnn+0aNGhVr4Yvff/9dycnJ+QJQ69atPY/n3jZo0ECxsbE+27Vs2fK4x1i3bp0iIiLUpk2b426blJSkgQMH6q677tI111yjOXPmKDIy0vN4ZmampLz5a8c6NpxFRUXphBNOKHDbJk2a5LuvRo0a2rlzpyRp27Zt2rVrl1588UW9+OKLBe5j69atx31NAIDSI2gBAMqsY8eOnlUHj+V2u+VyufT+++/7BI9cxwagKlWq+KXGQLjzzju1fft2TZo0SUOGDNG0adM8PV9ut1uSnadVv379fM89drXDmJiYfCEyV0Hvo2R7sryPdd111xU6F+3kk08uxisCAJQWQQsA4FfNmzeXMUbNmjVTixYtHNtv06ZN9cMPP8jtdvsEkl9//dXzeO7t4sWLtW/fPp9Qt2bNmmLV7na79csvv+jUU08tVl2PPvqoduzYoZdfflk1atTQE0884dmXJNWtW1c9evQo1r5Kq06dOqpevbpycnL8fiwAQMGYowUA8KsrrrhCkZGReuCBBzw9LrmMMdq+fXup9turVy9t3rxZb7zxhue+o0eP6plnnlFsbKy6devm2e7o0aOaOnWqZ7ucnBw988wzxz1Gnz59FBERofHjx3t6ibxrL8wLL7ygvn376sknn9SDDz4oSerZs6fi4uL08MMP68iRI/mes23btuPWU1yRkZG68sorNXfuXP30009+PRYAoGD0aAEA/Kp58+Z68MEHNXbsWGVlZalPnz6qXr261q9fr/nz52vo0KG6/fbbS7zfoUOH6oUXXlBKSopWrVqlxMREvfnmm1q+fLkmT56s6tWrS5IuvfRSnXXWWRozZoyysrLUpk0bzZs3T7t37z7uMZKSknTPPfdowoQJOvvss3XFFVcoJiZGX3/9tRo2bKiJEycW+LyIiAjNnDlTu3fv1rhx41SzZk0NGzZMU6dO1fXXX68OHTqof//+qlOnjjZs2KB3331XZ511lqZMmVLi96EwjzzyiD755BN16tRJQ4YMUZs2bbRjxw598803+uijj7Rjxw7HjgUAyI+gBQDwuzFjxqhFixZ66qmn9MADD0iSGjdurAsuuECXXXZZqfZZpUoVLVmyRGPGjNErr7yiPXv2qGXLlsrIyFBKSopnu4iICC1cuFC33nqrZsyYIZfLpcsuu0xPPPGE2rdvf9zjjB8/Xs2aNdMzzzyje+65R1WrVtXJJ5+s66+/vsjnRUdHa/78+erRo4duueUWJSQk6Nprr1XDhg31yCOP6LHHHtPhw4fVqFEjnX322cVaabEk6tWrp6+++krjx4/XvHnz9Nxzz6lWrVo66aST9Oijjzp6LABAfi5T1NgHAAAAAECJMUcLAAAAABxG0AIAAAAAhxG0AAAAAMBhBC0AAAAAcBhBCwAAAAAcRtACAAAAAIcRtAAAPlwul+6//36/HqN79+7q3r27Y/s7tub7779fLpdLf//9t1/2P336dLlcLmVlZXnuS0xM1CWXXOLI8ZyQlZUll8ul6dOnB7sUAKiQCFoAyq3chubjjz8e7FIkSSkpKXK5XIqLi9PBgwfzPZ6ZmSmXyxVSNQfbkiVL5HK59Oabbwa7FPy/3BB77FflypVLvK9du3apbt26Bf6Mc/+/FPb1559/SpIOHDigZ599VhdccIEaNGig6tWrq3379po6dapycnIcec0lNW3aNLVu3VqVK1dWcnKynnnmmQK3++ijj3Tuueeqdu3aSkhIUMeOHfXqq68GuFoA/hIV7AIAoCKJiorSgQMH9Pbbb6tfv34+j82cOVOVK1f
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADTjklEQVR4nOzdd3gU9drG8e9k0zudUAOEDkIARUAEFAQRARWpIqFZQNH3qEfxHPtRFOzHdpQqooBSRRQQRaqAFEFFCCWASCcJhPTdef8YNxAJECDJZDf357pybdiZ3b0Tkuw++3vmGcM0TRMREREREREpMD52BxAREREREfE2KrREREREREQKmAotERERERGRAqZCS0REREREpICp0BIRERERESlgKrREREREREQKmAotERERERGRAqZCS0REREREpICp0BIRERERESlgKrREREQ81LJlyzAMg2XLltkdRURE/kaFloiI2Gby5MkYhpHz4evrS+XKlYmLi+PAgQPn7N++fftc+/v7+1OjRg3uuece9u/ff8H7PvvjiSeeyNkvOjqabt26nfNYU6dOxeFw0KVLF9LT0wv+ixcREa/ma3cAERGR559/nho1apCens6PP/7I5MmTWblyJb/88guBgYG59q1SpQpjxowBIDMzk99++40PPviARYsWsW3bNoKDg/O877M1atTognmmTZtGXFwcHTt2ZO7cuedkEBERuRgVWiIiYrubb76ZFi1aADBs2DDKli3LK6+8wvz58+ndu3eufSMiIrjrrrtyXVejRg0eeOABVq1aRadOnc573/kxffp0Bg0axA033MC8efMKpMhKTU09pwAUERHvptZBEREpdtq2bQvArl278rV/xYoVAfD1vbL3D2fOnMldd91F+/btmT9//jlF1ieffELz5s0JCgqidOnS9O3b95yWxfbt29OoUSM2bNjA9ddfT3BwME8++SQJCQkYhsGrr77Khx9+SK1atQgICODqq69m/fr152T5/fff6dWrF6VLlyYwMJAWLVowf/78K/r6RESk6GhFS0REip2EhAQASpUqdc42p9PJsWPHAMjKymLbtm0888wzxMTE0KZNm3P2T05OztnfrWzZsufsN2vWLAYMGMD111/Pl19+SVBQUK7tL774Ik899RS9e/dm2LBhHD16lP/+979cf/31bNq0icjIyJx9jx8/zs0330zfvn256667qFChQs62Tz/9lFOnTnHvvfdiGAZjx47l9ttvZ/fu3fj5+QHw66+/0qZNGypXrswTTzxBSEgIM2fOpGfPnsyaNYvbbrstf99IERGxjQotERGxnbsYSk9PZ+3atTz33HMEBATkOaTi999/p1y5crmuq1+/PosXL8bf3/+c/Tt27HjOdaZp5vr3pk2bWLRoEddddx0LFiw4p8jau3cvzzzzDP/5z3948sknc66//fbbiY2N5b333st1/aFDh/jggw+49957c65zF4/79u0jPj4+p4isW7cuPXr0YNGiRTlf70MPPUS1atVYv349AQEBAIwYMYLrrruOxx9/XIWWiIgHUKElIiK2+3sxFB0dzSeffEKVKlXO2Tc6OpqPPvoIgOzsbLZv387YsWO5+eabWbFixTlF2LvvvkudOnUu+PgnTpwgOzubKlWqnFNkAcyePRuXy0Xv3r1zrY5VrFiR2rVr8/333+cqtAICAhg8eHCej9WnT59cK3XuNsndu3fnZPnuu+94/vnnOXXqFKdOncrZt3PnzjzzzDMcOHCAypUrX/BrEhERe6nQEhER27mLoeTkZCZOnMjy5ctzVnL+LiQkJFdh1qVLF6677jpatGjByy+/zGuvvZZr/2uuueaiwzBuvPFGqlWrxvvvv0/p0qV56623cm2Pj4/HNE1q166d5+3dLX9ulStXznN1DaBatWq5/u0uuhITEwHYuXMnpmny1FNP8dRTT+V5H0eOHFGhJSJSzKnQEhER251dDPXs2ZPrrruO/v37s337dkJDQy96++bNmxMREcHy5csvO8M777xDYmIib7/9NqVKleLZZ5/N2eZyuTAMg6+//hqHw3HObf+eMa9VMbe8bg9n2hldLhcAjz76KJ07d85z35iYmAt+LSIiYj8VWiIiUqw4HA7GjBlDhw4deOedd3KdXPhCnE4nKSkpl/24Pj4+fPzxxyQnJ/Pcc89RunRpRo0aBUCtWrUwTZMaNWpctA3xStWsWROwVsnyOr5MREQ8g8a7i4hIsdO+fXuuueYa3nzzTdLT0y+6//fff09KSgpNmjS5osf18/Pjiy++oE2bNjz88MNMnToVsIZeOBwOnnvuuXMGaZimyfHjx6/occ9Wvnx52rdvz//+9z8OHjx4zvajR48W2GOJiEjh0YqWiIgUS4899hh33nknkydP5r777su5Pjk5mU8++QQ4Mwzj/fffJygoKN+rXxcSHBzMV199Rbt27RgyZAgRERF0796d//znP4wePZqEhAR69uxJWFgYe/bsYc6cOdxzzz08+uijV/zYbu+++y7XXXcdjRs3Zvjw4dSsWZPDhw+zZs0a/vjjD37++ecCeywRESkcKrRERKRYuv3226lVqxavvvoqw4cPzzm26Y8//mDgwIEAGIZBqVKlaNeuHc888wxNmzYtkMeOiIjIGffep08fvv76a5544gnq1KnDG2+8wXPPPQdA1apVuemmm+jevXuBPK5bgwYN+Omnn3juueeYPHkyx48fp3z58sTGxvL0008X6GOJiEjhMMy/90CIiIiIiIjIFdExWiIiIiIiIgVMhZaIiIiIiEgBU6ElIiIiIiJSwFRoiYiIiIiIFDAVWiIiIiIiIgVMhZaIiIiIiEgBU6ElIiK5GIbBs88+W6iP0b59e9q3b19g9/f3zM8++yyGYXDs2LFCuf/JkydjGAYJCQk510VHR9OtW7cCebyCkJCQgGEYTJ482e4oIiIlkgotEfFa7hear776qt1RAIiLi8MwDMLDw0lLSztne3x8PIZhFKvMdlu2bBmGYfDFF1/YHUXOo1OnThiGwQMPPJCv/V0uFx988AFNmzYlNDSUChUqcPPNN7N69epc+7l/X873ceDAgTzvPykpifLly9v2c+NyuRg7diw1atQgMDCQq666is8+++yc/datW8eIESNo3rw5fn5+GIZR5FlFpHD52h1ARKQk8fX1JTU1lS+//JLevXvn2jZt2jQCAwNJT0+3KZ0lLS0NX1/Penoo6swDBw6kb9++BAQEFNljFkezZ89mzZo1l3Sbxx57jNdff5277rqLESNGkJSUxP/+9z/atWvHqlWruOaaawC499576dixY67bmqbJfffdR3R0NJUrV87z/p9++mlSU1Mv7wsqAP/61794+eWXGT58OFdffTXz5s2jf//+GIZB3759c/ZbuHAh48eP56qrrqJmzZrs2LHDtswiUji0oiUiUoQCAgK48cYb83yH+9NPP+WWW24p0MdLT0/H5XJd0m0CAwM9rtAq6swOh4PAwMASvQqRnp7OI488wuOPP57v22RnZ/P+++/Tq1cvpk6dyj333MM///lPvv32W7Kzs5k2bVrOvq1ateKuu+7K9VGjRg1SU1MZMGBAnvf/yy+/8P77719Spvx69tlniY6OvuA+Bw4c4LXXXmPkyJF8+OGHDB8+nC+//JK2bdvy2GOP4XQ6c/a9//77SU5O5qeffqJTp04FnldE7KdCS0RKvCNHjjB06FAqVKhAYGAgTZo0YcqUKefsd/z4cQYOHEh4eDiRkZEMGjSIn3/++ZKPg+nfvz9ff/01SUlJOdetX7+e+Ph4+vfvf87+J06c4NFHH6Vx48aEhoYSHh7OzTffzM8//5xrP3eb3fTp0/n3v/9N5cqVCQ4O5uTJkwB8/vnnNGjQgMDAQBo1asScOXOIi4s758Xj+Y532rlzJ3FxcURGRhIREcHgwYPPWTmYNGkSN9xwA+XLlycgIIAGDRrw/vvv5/t7c7nyc1zZ3r17iYmJoVGjRhw+fBiw2swefvhhqlatSkBAADExMbzyyisXLU7zOkbLbeXKlVxzzTUEBgZSs2ZNPv7443P22b17N3feeSe
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gp_linear.plot(X_1D_star, y_star_linear, var_linear)\n",
    "gp_poly.plot(X_1D_star, y_star_poly, var_poly)\n",
    "gp_periodic.plot(X_1D_star, y_star_periodic, var_periodic)\n",
    "gp_rbf.plot(X_1D_star, y_star_rbf, var_rbf)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Fit the 2D data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "outputs": [],
   "source": [
    "gp_linear.fit(X_2D, y_2D)\n",
    "gp_poly.fit(X_2D, y_2D)\n",
    "gp_periodic.fit(X_2D, y_2D)\n",
    "gp_rbf.fit(X_2D, y_2D)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict the values of the test data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "outputs": [],
   "source": [
    "y_star_linear, var_linear = gp_linear.predict(X_2D_star)\n",
    "y_star_poly, var_poly = gp_poly.predict(X_2D_star)\n",
    "y_star_periodic, var_periodic = gp_periodic.predict(X_2D_star)\n",
    "y_star_rbf, var_rbf = gp_rbf.predict(X_2D_star)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Plot the results"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_37332/3048407143.py:67: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
      "/tmp/ipykernel_37332/3048407143.py:68: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] + 1.96 * np.sqrt(np.diag(var)),\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1YAAANsCAYAAABYvTmgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACU10lEQVR4nOzde3zO9f/H8ee1I7OT80y2CTkrZ6NSOaxIfEuiFUP0y7koFJJjSCFfVDQdFpFIUpp9myKnSCgJjSnn0xxmZtf1+f1x5eKy0eaz7drhcb/drttcn8/n+nxe12Gu67n3+/O6LIZhGAIAAAAA3DI3VxcAAAAAAPkdwQoAAAAATCJYAQAAAIBJBCsAAAAAMIlgBQAAAAAmEawAAAAAwCSCFQAAAACYRLACAAAAAJMIVgAAAABgEsEKAJAr9u/fL4vFovnz57u6lAKNxxkAXINgBQAwbf78+bJYLPrpp59cXUq2uO+++1SrVq10y+Pi4uTj46N69erp1KlTLqgMAJBXebi6AABA4RAaGqqLFy/K09PT1aXckv/9739q166dqlatqtWrV6tEiRKuLgkAkIcwYgUAyBUWi0VFihSRu7u7q0u5oeTk5AyXr1mzRu3atdMdd9yRbaHqwoULpvcBAMg7CFYAgFyR0bk/UVFR8vX11d9//60OHTrI19dXpUuX1pAhQ2S1Wp1ub7PZNG3aNNWsWVNFihRR2bJl9eyzz+r06dNO233xxRdq27atgoOD5e3trUqVKmns2LHp9ndlut+WLVt07733ysfHRy+//HK6un/44Qe1bdtWlStX1urVq1WyZEmn9V9//bXuueceFStWTH5+fmrbtq1+/fVXp22u3M99+/apTZs28vPzU2RkpCR74OzXr5+WLVumWrVqydvbWzVr1tQ333yTrpa///5bPXr0UNmyZR3bvf/++//+4AMAchxTAQEALmW1WhUREaHGjRvrjTfe0OrVqzV16lRVqlRJzz33nGO7Z599VvPnz1f37t01YMAAJSQkaObMmfr555+1bt06xxTD+fPny9fXVy+88IJ8fX31v//9T6NGjdLZs2c1ZcoUp2OfPHlSDz30kDp37qynnnpKZcuWdVq/bt06tWnTRhUrVlRcXJxKlSrltP6jjz5St27dFBERoUmTJik5OVmzZ8/W3XffrZ9//llhYWGObdPS0hQREaG7775bb7zxhnx8fBzr1q5dq88//1x9+vSRn5+fZsyYoccee0yJiYmOIHf06FE1adLEEcRKly6tr7/+Wj179tTZs2c1aNCg7Hg6AAC3ygAAwKTo6GhDkrF58+YbbpOQkGBIMqKjox3LunXrZkgyxowZ47Rt3bp1jfr16zuu//DDD4YkIyYmxmm7b775Jt3y5OTkdMd+9tlnDR8fHyMlJcWxrHnz5oYkY86cOem2b968uVGiRAnDz8/PqFmzpnHs2LF025w7d84IDAw0evXq5bT8yJEjRkBAgNPyK/dz2LBh6fYjyfDy8jL27t3rWPbLL78Ykoy3337bsaxnz55GuXLljBMnTjjdvnPnzkZAQIDjfmf0OAMAch5TAQEALvd///d/Ttfvuece/fnnn47rixcvVkBAgFq1aqUTJ044LvXr15evr6++++47x7ZFixZ1/PvcuXM6ceKE7rnnHiUnJ+v33393Oo63t7e6d++eYU0XLlzQuXPnVLZsWfn7+6dbHxsbqzNnzqhLly5ONbm7u6tx48ZONV1x7QjctVq2bKlKlSo5rtepU0f+/v6Ox8AwDC1ZskTt2rWTYRhOx4uIiFBSUpK2bt2a4b4BALmDqYAAAJcqUqSISpcu7bSsePHiTudO7dmzR0lJSSpTpkyG+zh27Jjj37/++qtGjBih//3vfzp79qzTdklJSU7Xy5cvLy8vrwz3WblyZXXt2lVDhw5Vly5dtHjxYqfGG3v27JEkPfDAAxne/vow5uHhodtuuy3DbUNCQtItu/YxOH78uM6cOaN3331X7777bob7uPYxAADkPoIVAMClMtMl0GazqUyZMoqJiclw/ZVgdubMGTVv3lz+/v4aM2aMKlWqpCJFimjr1q0aOnSobDab0+2uHd3KyEsvvaSTJ09q8uTJ6tWrl+bNmyeLxeKoSbKfZxUUFJTuth4ezm+x3t7ecnPLeKLIjR4DwzCcjvXUU0+pW7duGW5bp06dm94XAEDOIlgBAPK8SpUqafXq1WrWrNlNw1B8fLxOnjypzz//XPfee69jeUJCwi0fe9KkSTp16pTmzp2r4sWLa+rUqY6aJKlMmTJq2bLlLe8/M0qXLi0/Pz9ZrdYcPxYA4NZwjhUAIM/r1KmTrFarxo4dm25dWlqazpw5I+nqyM+VkR5JSk1N1axZs0wd/5133lHHjh315ptvaty4cZKkiIgI+fv7a8KECbp8+XK62xw/ftzUMa/l7u6uxx57TEuWLNHOnTtz9FgAgFvDiBUAINu8//77GX7/0sCBA03tt3nz5nr22Wc1ceJEbdu2Ta1bt5anp6f27NmjxYsXa/r06erYsaOaNm2q4sWLq1u3bhowYIAsFos++ugjp6B1K9zc3BQTE6OkpCSNHDlSJUqUUJ8+fTR79mw9/fTTqlevnjp37qzSpUsrMTFRX331lZo1a6aZM2eaOu61Xn/9dX333Xdq3LixevXqpRo1aujUqVPaunWrVq9erVOnTmXbsQAAWUewAgBkm9mzZ2e4PCoqyvS+58yZo/r16+udd97Ryy+/LA8PD4WFhempp55Ss2bNJEklS5bUihUrNHjwYI0YMULFixfXU089pRYtWigiIsLU8b28vLR06VK1bNlS/fv3V2BgoJ588kkFBwfr9ddf15QpU3Tp0iWVL19e99xzzw27Dd6qsmXLatOmTRozZow+//xzzZo1SyVLllTNmjU1adKkbD0WACDrLIbZP+MBAAAAQCHHOVYAAAAAYBLBCgAAAABMIlgBAAAAgEkEKwAAAAAwiWAFAAAAACYRrAAAAADAJIIVABRiFotFo0ePztFj3HfffbrvvvuybX/X1zx69GhZLBadOHEiR/Y/f/58WSwW7d+/37EsLCxMDz/8cLYcLzvs379fFotF8+fPd3UpAFBoEawAFAhXPli+8cYbri5Fkv0LcS0Wi/z9/XXx4sV06/fs2SOLxZKnana1+Ph4WSwWffbZZ64uBTfQqlUrWSwW9evXL1Pbf/vtt+rZs6dq1aold3d3hYWFZbjdlXB8o8u6deuctt+1a5cefPBB+fr6qkSJEnr66ad1/Phxs3cvy2bPnq3HH39cISEhslgsN/wi7MOHD2vYsGG6//775efnJ4vFovj4+FytFUDO83B1AQBQUHl4eCg5OVlffvmlOnXq5LQuJiZGRYoUUUpKiouqs7t48aI8PPLXW0Fu1/z000+rc+fO8vb2zrVj5kWff/651q9fn6XbfPLJJ/r0009Vr149BQcH33C7Rx99VJUrV063/OWXX9b58+fVsGFDx7K//vpL9957rwICAjRhwgSdP39eb7zxhnbs2KFNmzbJy8srSzWaMWnSJJ07d06NGjXS4cOHb7jd7t27NWnSJFWpUkW1a9fO8uMIIH/IX++mAJCPeHt7q1mzZlqwYEG6YPXJJ5+obdu2WrJkSbYdLyUlRV5eXnJzy/xkhCJFimTb8XNLbtfs7u4ud3f3XD1mXpOSkqLBgwdr6NChGjVqVKZvN2HCBL333nvy9PTUww8/rJ07d2a4XZ06dVSnTh2nZQcPHtRff/2lZ555xiksTZgwQRcuXNCWLVsUEhIiSWrUqJFatWql+fPnq3fv3rdwD53Nnz9f3bt3l2EYN91uzZo1jtEqX1/fG25Xv359nTx5UiVKlNBnn32mxx9/3HSNAPIepgICKFSOHTumnj17qmzZsipSpIjuvPNOffDBB+m2O3nypJ5++mn5+/srMDBQ3bp10y+//JLl81iefPJJff311zpz5oxj2ebNm7Vnzx49+eST6bY/deqUhgwZotq1a8vX11f+/v566KGH9Msvvzhtd2Xa3MKFCzVixAiVL19ePj4+Onv2rCRp8eLFqlGjhooUKaJatWpp6dK
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACkZElEQVR4nOzdeZyN5f/H8feZfZjNPghji7EVQqhG9qWiLNmNtSyhaPENyZJSRAmJhmqsiVRKyBKRUqKShmzJvsw0xqzn/v1x/5ycM4PBzNyzvJ6Px3mc69z3fe77c86cU+ftuu7rthmGYQgAAAAA4OBmdQEAAAAAkN0QlAAAAADABUEJAAAAAFwQlAAAAADABUEJAAAAAFwQlAAAAADABUEJAAAAAFwQlAAAAADABUEJAAAAAFwQlAAAaWrUqJEaNWpkdRkZYsGCBbLZbDp8+PBNPzc8PFwhISEZXpNVbue9AIC8hKAEALnElR/AV24+Pj668847NWTIEJ06dcrq8nIVm82mIUOGpFr+yiuvyGazqU+fPrLb7RZUBgDIKB5WFwAAyFjjx49X2bJlFR8fr61bt2r27Nlas2aNfv31V+XLl8/q8izRo0cPde7cWd7e3pl2jFdffVUvvviievXqpXnz5snNjX+LBICcjKAEALlMq1atdM8990iS+vXrp0KFCmnatGn69NNP1aVLF4urs4a7u7vc3d0zbf+vv/66Ro0apZ49e+r999+/7ZBkGIbi4+Pl6+ubQRUCAG4W/9wFALlc48aNJUmHDh2SJCUnJ2vChAkqX768vL29FRISov/9739KSEi45j5iY2OVP39+DRs2LNW6v//+W+7u7po8ebKk/4YAbtu2Tc8884yKFCmi/Pnz69FHH9WZM2dSPX/WrFmqWrWqvL29VaJECQ0ePFgXL1502qZRo0aqVq2a9uzZo7CwMOXLl08VKlTQxx9/LEnavHmz6tWrJ19fX1WqVEnr1693en5a5+V8+umnatOmjUqUKCFvb2+VL19eEyZMUEpKyo3f1KtMmzZNzz33nLp3766IiAinkGS32zV9+nRVrVpVPj4+KlasmJ544glduHDBaR8hISF66KGHtHbtWt1zzz3y9fXVu+++q02bNslms2nZsmWaNGmS7rjjDvn4+KhJkyY6cOBAqlq+//57tWzZUoGBgcqXL5/CwsK0bdu2m3o9AAATQQkAcrmDBw9KkgoVKiTJ7GUaO3asatWqpTfffFNhYWGaPHmyOnfufM19+Pn56dFHH9XSpUtTBYnFixfLMAx169bNaflTTz2lX375RS+99JIGDhyozz77LNV5PePGjdPgwYNVokQJTZ06Ve3bt9e7776r5s2bKykpyWnbCxcu6KGHHlK9evU0ZcoUeXt7q3Pnzlq6dKk6d+6s1q1b69VXX9WlS5fUoUMH/fvvv9d9XxYsWCA/Pz8988wzmjFjhmrXrq2xY8fqhRdeuP4bepUZM2ZoxIgR6tq1qxYsWJCqJ+mJJ57Qs88+q4YNG2rGjBnq3bu3IiMj1aJFi1Svb//+/erSpYuaNWumGTNm6O6773ase/XVV7Vy5UqNHDlSo0aN0o4dO1K93998840eeOABxcTE6KWXXtIrr7yiixcvqnHjxtq5c2e6XxMA4P8ZAIBcISIiwpBkrF+/3jhz5oxx7NgxY8mSJUahQoUMX19f4++//zZ2795tSDL69evn9NyRI0cakoxvvvnGsSwsLMwICwtzPF67dq0hyfjyyy+dnlujRg2n7a7U0bRpU8NutzuWP/3004a7u7tx8eJFwzAM4/Tp04aXl5fRvHlzIyUlxbHdzJkzDUnG+++/71SLJGPRokWOZX/88YchyXBzczN27NiRqs6IiIhUNR06dMixLC4uLtV7+MQTTxj58uUz4uPjHct69epllClTxmk7SUaZMmUMSUaXLl2M5OTkVPv69ttvDUlGZGSk0/Kvvvoq1fIr+/rqq6+ctt24caMhyQgNDTUSEhIcy2fMmGFIMvbu3WsYhmHY7XajYsWKRosWLZze87i4OKNs2bJGs2bNrvteAABSo0cJAHKZpk2bqkiRIipVqpQ6d+4sPz8/rVy5UiVLltSaNWskSc8884zTc0aMGCFJ+uKLL6673xIlSigyMtKx7Ndff9WePXvUvXv3VNsPGDBANpvN8fj+++9XSkqKjhw5Iklav369EhMTNXz4cKeemP79+ysgICBVLX5+fk69XpUqVVJQUJBCQ0NVr149x/Ir7b/++uuar0WS0/k///77r86ePav7779fcXFx+uOPP677XEmOmQTLli2b5vlPy5cvV2BgoJo1a6azZ886brVr15afn582btzotH3ZsmXVokWLNI/Vu3dveXl5OR7ff//9Tq9x9+7dioqKUteuXXXu3DnHsS5duqQmTZpoy5YtzMIHADeJyRwAIJd55513dOedd8rDw0PFihVTpUqVHEHkyJEjcnNzU4UKFZyeExwcrKCgIEeISYubm5u6deum2bNnKy4uTvny5VNkZKR8fHzUsWPHVNuXLl3a6XGBAgUkyXF+zpVjVapUyWk7Ly8vlStXLlUtd9xxh1PwkqTAwECVKlUq1bKrj3Mtv/32m0aPHq1vvvlGMTExTuuio6Ov+1xJ6tWrl/755x+98sorKly4sJ5++mmn9VFRUYqOjlbRokXTfP7p06edHpctW/aax7rRexkVFeWo6Vqio6MdzwMA3BhBCQBymbp16zpmvbsW18CRXj179tTrr7+uVatWqUuXLlq0aJEeeughRzi52rVmmTMM45aOfa393cpxLl68qLCwMAUEBGj8+PEqX768fHx89NNPP+n5559PV++Lh4eHli1bppYtW2rEiBEKCgpS7969HevtdruKFi3q1AN3tSJFijg9vt4Mdzd6jVfqff31153Obbqan5/fNfcPAEiNoAQAeUiZMmVkt9sVFRWl0NBQx/JTp07p4sWLKlOmzHWfX61aNdWsWVORkZG64447dPToUb399tu3XItkTmJQrlw5x/LExEQdOnRITZs2vaX9psemTZt07tw5ffLJJ3rggQccy6/MDJhePj4+Wr16tR588EH1799fQUFBevTRRyVJ5cuX1/r169WwYcNMn+a7fPnykqSAgIBMfd8AIC/hHCUAyENat24tSZo+fbrT8mnTpkmS2rRpc8N99OjRQ19//bWmT5+uQoUKqVWrVrdUS9OmTeXl5aW33nrLqfdn/vz5io6OTlctt+pKD83Vx01MTNSsWbNuel8BAQH66quvVKFCBXXp0kUbNmyQJHXq1EkpKSmaMGFCquckJyenmgL9dtSuXVvly5fXG2+8odjY2FTr05qWHQBwffQoAUAectddd6lXr16aO3euY/jZzp07tXDhQrVr104PPvjgDffRtWtXPffcc1q5cqUGDhwoT0/PW6qlSJEiGjVqlF5++WW1bNlSjzzyiPbv369Zs2apTp06aU4QkVEaNGigAgUKqFevXho6dKhsNps+/PDDWx4WWKRIEa1bt04NGzZUu3bttGHDBoWFhemJJ57Q5MmTtXv3bjVv3lyenp6KiorS8uXLNWPGDHXo0CFDXo+bm5vmzZunVq1aqWrVqurdu7dKliyp48ePa+PGjQoICNBnn32WIccCgLyCoAQAecy8efNUrlw5LViwQCtXrlRwcLBGjRqll156KV3PL1asmJo3b641a9aoR48et1XLuHHjVKRIEc2cOVNPP/20ChYsqAEDBuiVV1655QCWHoUKFdLnn3+uESNGaPTo0SpQoIC6d++uJk2aXHPmuRspVaqUvv76a91///1q1aqVtmzZojlz5qh27dp699139b///U8eHh4KCQlR9+7d1bBhwwx9TY0aNdL27ds1YcIEzZw5U7GxsQoODla9evX0xBNPZOixACAvsBm3+s9nAIA869FHH9XevXt14MABq0sBACBTcI4SAOCmnDhxQl988cVt9yYBAJCdMfQOAJAuhw4d0rZt2zRv3jx5enoynAsAkKvRowQASJfNmzerR48eOnTokBYuXKjg4GCrSwIAINNwjhIAAAAAuKBHCQAAAABcEJQAAAA
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC1+klEQVR4nOzde3yP9f/H8ednm53sZE4bmS1zmFM555SzOSQSUsKifL8hKfrGr0gOiaKkLx3UVJRCpJBYJmcip5JGRDmGbTnMbJ/r98f19WGfDRvbrh0e99vtc3tfn+u6Ptf1+nw29Xnuel/vt80wDEMAAAAAAAcXqwsAAAAAgLyGoAQAAAAATghKAAAAAOCEoAQAAAAATghKAAAAAOCEoAQAAAAATghKAAAAAOCEoAQAAAAATghKAAAAAOCEoAQAyBXNmzdX8+bNs+14hw4dks1m0+zZsx3rxowZI5vNlm3nyA8y+hwAALePoAQAhdjs2bNls9kcD09PT1WqVEmDBw/WiRMnrC4v1zRv3lzVq1dPtz4mJkbe3t6qXbu2zpw5Y0FlAACruFldAADAemPHjlVYWJiSkpK0bt06zZw5U8uWLdOePXvk7e2dLef47rvvsuU4N/Liiy9qxIgR2XKs77//Xp06dVLlypW1atUqBQYGZstxAQD5A0EJAKD27durbt26kqTHH39cxYsX19SpU/XVV1/p4Ycfvq1jX7hwQd7e3nJ3d8+OUm/Izc1Nbm63/7+2NWvWqFOnTqpUqVK2haTz58+raNGit30cAEDuoOsdACCdli1bSpIOHjzoWDdnzhzVqVNHXl5eCgwMVM+ePXXkyJE0r7vShW3btm2699575e3trf/7v/9zbHO+R+nkyZPq37+/SpcuLU9PT91111366KOP0tUTHx+vqKgo+fv7KyAgQH379lV8fHy6/a53j9KcOXNUv359eXt7q1ixYrr33nuve4Vr7dq16tixo8LDw7Vq1SoVL148zfbly5eradOmKlq0qHx9fdWxY0f9/PPPafaJioqSj4+PDhw4oA4dOsjX11e9evWSJNlsNg0ePFiLFy9W9erV5eHhoWrVqunbb79NV8tff/2lfv36qXTp0o79PvzwwwzrBgBkL64oAQDSOXDggCQ5QsKECRM0atQo9ejRQ48//rhOnTql6dOn695779VPP/2kgIAAx2tPnz6t9u3bq2fPnnr00UdVunTpDM9x8eJFNW/eXPv379fgwYMVFham+fPnKyoqSvHx8Xr66aclSYZhqHPnzlq3bp3+/e9/KyIiQosWLVLfvn0z9V5efvlljRkzRo0aNdLYsWPl7u6uzZs36/vvv1fbtm3T7Lt+/Xp16NBBYWFhiomJUYkSJdJs/+STT9S3b19FRkZq0qRJunDhgmbOnKkmTZrop59+UmhoqGPflJQURUZGqkmTJnr99dfTdGFct26dvvzySw0cOFC+vr5666239OCDD+rw4cOOz/zEiRO65557HMGqZMmSWr58ufr376/ExEQNHTo0U+8fAHCLDABAoRUdHW1IMlatWmWcOnXKOHLkiDFv3jyjePHihpeXl/Hnn38ahw4dMlxdXY0JEyakee3u3bsNNze3NOubNWtmSDLeeeeddOdq1qyZ0axZM8fzN99805BkzJkzx7EuOTnZaNiwoeHj42MkJiYahmEYixcvNiQZkydPduyXkpJiNG3a1JBkREdHO9a/9NJLxrX/a4uLizNcXFyMBx54wEhNTU1Tj91uT1NbYGCg4evra1SrVs04efJkuvr/+ecfIyAgwHjiiSfSrD9+/Ljh7++fZn3fvn0NScaIESPSHUeS4e7ubuzfv9+xbufOnYYkY/r06Y51/fv3N4KDg42///47zet79uxp+Pv7GxcuXDAMwzAOHjyY7nMAANw+ut4BANS6dWuVLFlS5cqVU8+ePeXj46NFixapbNmy+vLLL2W329WjRw/9/fffjkdQUJAqVqyo1atXpzmWh4eHHnvssZuec9myZQoKCkpzD1SRIkU0ZMgQnTt3TmvWrHHs5+bmpieffNKxn6urq5566qmbnmPx4sWy2+0aPXq0XFzS/i/PuYve+fPn9c8//6h06dLy8/NLd6yVK1cqPj5eDz/8cJrPwdXVVQ0aNEj3OUhKU/O1WrdurQoVKjie16xZU35+fvr9998lmVfRFi5cqE6dOskwjDTni4yMVEJCgrZv337T9w8AuHV0vQMA6L///a8qVaokNzc3lS5dWpUrV3YEi7i4OBmGoYoVK2b42iJFiqR5XrZs2UwN3PDHH3+oYsWK6QJMRESEY/uVNjg4WD4+Pmn2q1y58k3PceDAAbm4uKhq1ao33Tc8PFx9+vTR888/r4cffljz58+Xq6urY3tcXJykq/dvOXMOV25ubrrjjjsy3DckJCTdumLFiuns2bOSpFOnTik+Pl7vvfee3nvvvQyPcfLkyZu+JwDArSMoAQBUv359x6h3zux2u2w2m5YvX54mOFzhHGC8vLxypMbc8J///EenT5/W5MmT9cQTT+iDDz5wXHmy2+2SzPuUgoKC0r3WebQ9Dw+PdCHwiow+R8m8knTtuR599NHr3otVs2bNTLwjAMCtIigBAG6oQoUKMgxDYWFhqlSpUrYdt3z58tq1a5fsdnuaQPHrr786tl9pY2JidO7cuTShbN++fZmq3W6365dfftHdd9+dqbomTZqkM2fOaNasWSpWrJimTJniOJYklSpVSq1bt87UsW5VyZIl5evrq9TU1Bw/FwAgY9yjBAC4oa5du8rV1VUvv/yy44rHFYZh6PTp07d03A4dOuj48eP6/PPPHetSUlI0ffp0+fj4qFmzZo79UlJSNHPmTMd+qampmj59+k3P0aVLF7m4uGjs2LGOqzTX1n497777rrp166apU6dq/PjxkqTIyEj5+fnplVde0eXLl9O95tSpUzetJ7NcXV314IMPauHChdqzZ0+OngsAkDGuKAEAbqhChQoaP368Ro4cqUOHDqlLly7y9fXVwYMHtWjRIg0YMEDDhw/P8nEHDBigd999V1FRUdq2bZtCQ0O1YMECrV+/Xm+++aZ8fX0lSZ06dVLjxo01YsQIHTp0SFWrVtWXX36phISEm54jPDxcL7zwgsaNG6emTZuqa9eu8vDw0NatW1WmTBlNnDgxw9e5uLho7ty5SkhI0KhRoxQYGKiBAwdq5syZ6t27t2rXrq2ePXuqZMmSOnz4sJYuXarGjRvr7bffzvLncD2vvvqqVq9erQYNGuiJJ55Q1apVdebMGW3fvl2rVq3SmTNnsu1cAID0CEoAgJsaMWKEKlWqpDfeeEMvv/yyJKlcuXJq27at7r///ls6ppeXl2JjYzVixAh99NFHSkxMVOXKlRUdHa2oqCjHfi4uLlqyZImGDh2qOXPmyGaz6f7779eUKVNUq1atm55n7NixCgsL0/Tp0/XCCy/I29tbNWvWVO/evW/4Ond3dy1atEitW7fWU089pYCAAD3yyCMqU6aMXn31Vb322mu6dOmSypYtq6ZNm2ZqpL+sKF26tLZs2aKxY8fqyy+/1IwZM1S8eHFVq1ZNkyZNytZzAQDSsxk36nsAAAAAAIUQ9ygBAAAAgBOCEgAAAAA4ISgBAAAAgBOCEgAAAAA4ISgBAAAAgBOCEgAAAAA4ISgBQAFjs9k0ZsyYHD1H8+bN1bx582w7nnPNY8aMkc1m099//50jx589e7ZsNpsOHTrkWBcaGqr77rsvW86XHQ4dOiSbzabZs2dbXQoAFEoEJQB51pUviq+//rrVpUiSoqKiZLPZ5Ofnp4sXL6bbHhcXJ5vNlqdqtlpsbKxsNpsWLFhgdSn4n3379umZZ55Ro0aN5OnpmS4w3syV3/GMHm3atEm3/4EDB/TII4+oVKlS8vLyUsWKFfXCCy+k2++LL77QPffco4CAABUvXlzNmjXT0qVLb+et3pKZM2eqe/fuCgkJkc1mSzP58bWOHTumESNGqEWLFvL19ZXNZlNsbGyu1gogZ7lZXQAA5Cdubm66cOGCvv76a/Xo0SPNtrlz58rT01NJSUkWVWe6ePGi3Nzy13/ec7v
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAClIklEQVR4nOzdeZxOdf/H8fc1M2bfbbOYMRjLWCp7liwlSyqSEMVEdIdw31RUJEtKG93dqBQVkRJJSnJHCCGikpAJMcgyY4xZr/P74/zmul0zg8GMM9fM6/l4XI/PdZ1zrnM+16Ku95xzvsdmGIYhAAAAAICDm9UNAAAAAEBxQ1ACAAAAgFwISgAAAACQC0EJAAAAAHIhKAEAAABALgQlAAAAAMiFoAQAAAAAuRCUAAAAACAXghIAAAAA5EJQAgDAImvWrJHNZtOaNWusbgUAkAtBCQBw1ebOnSubzea4eXh4KDIyUvHx8frrr7/yLN+mTRun5T09PVWlShUNGjRIhw4duuS6L7yNHj3asVxMTIzuvPPOPNv64IMP5O7uro4dOyotLa3wXzwAoETzsLoBAIDrmzBhgqpUqaK0tDRt2rRJc+fO1fr16/Xzzz/L29vbadlKlSppypQpkqSMjAz9+uuvmjVrllauXKndu3fL19c333VfqG7dupfsZ/78+YqPj1e7du20dOnSPD0AAHA5BCUAwDXr1KmTGjVqJEl6+OGHVa5cOb344otatmyZevTo4bRsUFCQHnjgAadpVapU0dChQ7VhwwbdfvvtF113QSxcuFD9+vXTrbfeqs8++6xQQlJqamqeAAcAKNk49A4AUOhuueUWSdL+/fsLtHxYWJgkycPj2v5+t2jRIj3wwANq06aNli1blickzZs3Tw0bNpSPj49CQ0PVq1evPIf8tWnTRnXr1tW2bdvUqlUr+fr66qmnnlJCQoJsNptefvllvfXWW6pWrZq8vLzUuHFjbdmyJU8vv/32m7p3767Q0FB5e3urUaNGWrZs2TW9PgDA9cMeJQBAoUtISJAkhYSE5JmXnZ2tv//+W5KUmZmp3bt369lnn1VsbKxatGiRZ/mkpCTH8jnKlSuXZ7nFixerT58+atWqlT7//HP5+Pg4zZ88ebLGjh2rHj166OGHH9aJEyf073//W61atdL27dsVHBzsWPbkyZPq1KmTevXqpQceeEAVK1Z0zPvwww919uxZPfLII7LZbJo6daq6deumP/74Q2XKlJEk/fLLL2rRooUiIyM1evRo+fn5adGiReratasWL16se+65p2BvJADAMgQlAMA1ywkzaWlp2rx5s5577jl5eXnlO8jCb7/9pvLlyztNi4uL09dffy1PT888y7dr1y7PNMMwnB5v375dK1euVMuWLbV8+fI8IenPP//Us88+q0mTJumpp55yTO/WrZvq16+vGTNmOE1PTEzUrFmz9Mgjjzim5YS/gwcPau/evY4QWLNmTXXp0kUrV650vN7hw4crOjpaW7ZskZeXlyRp8ODBatmypZ588kmCEgC4AIISAOCa5Q4zMTExmjdvnipVqpRn2ZiYGL399tuSpKysLO3Zs0dTp05Vp06dtG7dujwh6j//+Y9q1Khxye2fOnVKWVlZqlSpUp6QJEmffvqp7Ha7evTo4bR3KiwsTNWrV9e3337rFJS8vLz00EMP5butnj17Ou0pyznM8I8//nD08t///lcTJkzQ2bNndfbsWceyHTp00LPPPqu//vpLkZGRl3xNAABrEZQAANcsJ8wkJSXp3Xff1XfffefYk5Kbn5+fU7Dq2LGjWrZsqUaNGumFF17QK6+84rR8kyZNLjuYw2233abo6GjNnDlToaGhmj59utP8vXv3yjAMVa9ePd/n5xwylyMyMjLfvVuSFB0d7fQ4JzSdPn1akrRv3z4ZhqGxY8dq7Nix+a7j+PHjBCUAKOYISgCAa3ZhmOnatatatmyp3r17a8+ePfL397/s8xs2bKigoCB99913V93DG2+8odOnT+v1119XSEiIxo8f75hnt9tls9n05Zdfyt3dPc9zc/eY316pHPk9X/rf4YB2u12SNGrUKHXo0CHfZWNjYy/5WgAA1iMoAQAKlbu7u6ZMmaK2bdvqjTfecLo47KVkZ2crJSXlqrfr5uam999/X0lJSXruuecUGhqqYcOGSZKqVasmwzBUpUqVyx7Gd62qVq0qydxLld/5VQAA18Dw4ACAQtemTRs1adJE06ZNU1pa2mWX//bbb5WSkqIbb7zxmrZbpkwZffLJJ2rRooVGjBihDz74QJI5aIO7u7uee+65PANBGIahkydPXtN2L1ShQgW1adNGb775po4ePZpn/okTJwptWwCAosMeJQBAkXj88cd13333ae7cufrHP/7hmJ6UlKR58+ZJ+t9gDjNnzpSPj0+B9z5diq+vr7744gu1bt1a/fv3V1BQkO6++25NmjRJY8aMUUJCgrp27aqAgAAdOHBAS5Ys0aBBgzRq1Khr3naO//znP2rZsqXq1aungQMHqmrVqjp27Jg2btyow4cP66effiq0bQEAigZBCQBQJLp166Zq1arp5Zdf1sCBAx3n9hw+fFgPPvigJMlmsykkJEStW7fWs88+q5tuuqlQth0UFOQYLrxnz5768ssvNXr0aNWoUUOvvfaannvuOUlSVFSU2rdvr7vvvrtQtpujdu3a2rp1q5577jnNnTtXJ0+eVIUKFVS/fn2NGzeuULcFACgaNiP3MQgAAAAAUMpxjhIAAAAA5EJQAgAAAIBcCEoAAAAAkAtBCQAAAAByISgBAAAAQC4EJQAAAADIhaAEACWMzWbT+PHji3Qbbdq0UZs2bQptfbl7Hj9+vGw2m/7+++8iWf/cuXNls9mUkJDgmBYTE6M777yzULZXGBISEmSz2TR37lyrWwGAUomgBKDYyvmh+PLLL1vdiiQpPj5eNptNgYGBOn/+fJ75e/fulc1mK1Y9W23NmjWy2Wz65JNPrG4FF3H77bfLZrNp6NChV/zcM2fOqEKFCvl+xr/88ovuu+8+Va1aVb6+vipXrpxatWqlzz//PM96cv5t5b7VqlXrql/XtXjnnXcUFxcnb29vVa9eXf/+97/zLLNkyRJ16NBBERER8vLyUqVKldS9e3f9/PPPFnQMoCh4WN0AALgSDw8Ppaam6vPPP1ePHj2c5s2fP1/e3t5KS0uzqDvT+fPn5eHhWv95v949P/jgg+rVq5e8vLyu2zaLo08//VQbN2686uePGzdOqamp+c77888/dfbsWfXr108RERFKTU3V4sWLdffdd+vNN9/UoEGDnJb38vLS7NmznaYFBQVddW9X680339Q//vEP3XvvvfrXv/6ldevWadiwYUpNTdWTTz7pWG7Xrl0KCQnR8OHDVa5cOSUmJurdd99VkyZNtHHjRt14443XvXcAhcu1/k8KABbz8vJSixYttGDBgjxB6cMPP1Tnzp21ePHiQtteWlqaPD095eZW8AMAvL29C23718v17tnd3V3u7u7XdZvFTVpamkaOHKknn3xS48aNu+Ln//zzz5o5c6bGjRuX7/PvuOMO3XHHHU7Thg4dqoYNG+rVV1/NE5Q8PDz0wAMPXHEfBTV37lw99NBDMgzjosucP39eTz/9tDp37uzYQzZw4EDZ7XZNnDhRgwYNUkhIiCTl+5offvhhVapUSTNnztSsWbOK5oUAuG449A6Ayzt+/LgGDBigihUrytvbWzfeeKPee++9PMudPHlSDz74oAIDAxUcHKx+/frpp59+uuLzQHr37q0vv/xSZ86ccUzbsmWL9u7dq969e+dZ/tSpUxo1apTq1asnf39/BQYGqlOnTvrpp5+clss5TG3hwoV65plnFBkZKV9fXyUnJ0uSPv74Y9WuXVve3t6qW7eulixZovj4eMXExDit52Ln++zbt0/x8fEKDg5WUFCQHnrooTx7A+bMmaNbb71VFSpUkJeXl2rXrq2ZM2cW+L25WgU5r+rPP/9UbGys6tatq2PHjkkyD/0aMWKEoqKi5OXlpdjYWL344ouy2+2XXFd+5yjlWL9+vZo0aSJvb29VrVpV77/
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gp_linear.plot(X_2D_star, y_star_linear, var_linear)\n",
    "gp_poly.plot(X_2D_star, y_star_poly, var_poly)\n",
    "gp_periodic.plot(X_2D_star, y_star_periodic, var_periodic)\n",
    "gp_rbf.plot(X_2D_star, y_star_rbf, var_rbf)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Fit the 3D data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "outputs": [],
   "source": [
    "gp_linear.fit(X_3D, y_3D)\n",
    "gp_poly.fit(X_3D, y_3D)\n",
    "gp_periodic.fit(X_3D, y_3D)\n",
    "gp_rbf.fit(X_3D, y_3D)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict the values of the test data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "outputs": [],
   "source": [
    "y_star_linear, var_linear = gp_linear.predict(X_3D_star)\n",
    "y_star_poly, var_poly = gp_poly.predict(X_3D_star)\n",
    "y_star_periodic, var_periodic = gp_periodic.predict(X_3D_star)\n",
    "y_star_rbf, var_rbf = gp_rbf.predict(X_3D_star)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Plot the results"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_37332/3048407143.py:67: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
      "/tmp/ipykernel_37332/3048407143.py:68: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] + 1.96 * np.sqrt(np.diag(var)),\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA00AAANsCAYAAACH80O/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACed0lEQVR4nOzdeVxVdf7H8fcV2VdxYwlBxQWXHJc0Ndcst0wrNY1U1LRxSZ3RqZyWscyszNLGUWsssMI1l8xMUxMXxrQsl4oMDdJyTQNEBBHO74/78+YVOKICF/D1fDzu43bP93vP+dzDZYa33+/5HothGIYAAAAAAPmq4OgCAAAAAKA0IzQBAAAAgAlCEwAAAACYIDQBAAAAgAlCEwAAAACYIDQBAAAAgAlCEwAAAACYIDQBAAAAgAlCEwAAAACYIDQBAIpEcnKyLBaLYmJiHF1KucZ5BoCSR2gCAFxTTEyMLBaLvv76a0eXUiQ6duyoRo0a5dm+efNmeXh4qFmzZjp79qwDKgMAlEYVHV0AAKB8CA0N1YULF+Ts7OzoUm7IF198oV69eqlevXratGmT/P39HV0SAKCUYKQJAFAkLBaL3Nzc5OTk5OhSCpSRkZHv9q1bt6pXr16qW7dukQWm8+fP3/Q+AAClA6EJAFAk8rvWJioqSl5eXvrtt9/Up08feXl5qWrVqpo0aZJycnLs3p+bm6tZs2apYcOGcnNzU/Xq1fX444/rjz/+sOv38ccfq2fPngoKCpKrq6tq166tqVOn5tnf5Sl4e/bsUfv27eXh4aF//vOfeerevn27evbsqfDwcG3atEmVK1e2a//ss8/Url07eXp6ytvbWz179tT3339v1+fy5zx8+LB69Oghb29vRUZGSrKGybFjx2r16tVq1KiRXF1d1bBhQ61fvz5PLb/99puGDRum6tWr2/q999571z75AIBixfQ8AECxysnJUdeuXdWqVSu9/vrr2rRpk2bOnKnatWtr1KhRtn6PP/64YmJiNHToUI0bN05JSUmaM2eOvv32W8XHx9um/cXExMjLy0t///vf5eXlpS+++ELPP/+80tLSNGPGDLtjnzlzRt27d9eAAQP06KOPqnr16nbt8fHx6tGjh2rWrKnNmzerSpUqdu0ffPCBhgwZoq5du+rVV19VRkaG5s2bp7vuukvffvutwsLCbH0vXbqkrl276q677tLrr78uDw8PW9uOHTu0cuVKjR49Wt7e3nrrrbf00EMP6ciRI7aQdvLkSd155522kFW1alV99tlnGj58uNLS0jRhwoSi+HEAAG6EAQDANURHRxuSjK+++qrAPklJSYYkIzo62rZtyJAhhiTjxRdftOvbtGlTo3nz5rbX27dvNyQZsbGxdv3Wr1+fZ3tGRkaeYz/++OOGh4eHkZmZadvWoUMHQ5Ixf/78PP07dOhg+Pv7G97e3kbDhg2NU6dO5elz7tw5w8/PzxgxYoTd9hMnThi+vr522y9/zqeffjrPfiQZLi4uxqFDh2zb9u3bZ0gy/v3vf9u2DR8+3AgMDDR+//13u/cPGDDA8PX1tX3u/M4zAKB4MT0PAFDs/vrXv9q9bteunX7++Wfb6+XLl8vX11f33HOPfv/9d9ujefPm8vLy0pYtW2x93d3dbf997tw5/f7772rXrp0yMjL0448/2h3H1dVVQ4cOzbem8+fP69y5c6pevbp8fHzytG/cuFEpKSkaOHCgXU1OTk5q1aqVXU2XXTlydqUuXbqodu3atte33367fHx8bOfAMAytWLFCvXr1kmEYdsfr2rWrUlNT9c033+S7bwBA8WN6HgCgWLm5ualq1ap22ypVqmR3rVJiYqJSU1NVrVq1fPdx6tQp239///33evbZZ/XFF18oLS3Nrl9qaqrd6+DgYLm4uOS7z/DwcA0ePFhPPfWUBg4cqOXLl9stYpGYmChJ6ty5c77vvzpoVaxYUbfddlu+fWvUqJFn25Xn4PTp00pJSdE777yjd955J999XHkOAAAli9AEAChWhVlNLzc3V9WqVVNsbGy+7ZdDV0pKijp06CAfHx+9+OKLql27ttzc3PTNN9/oqaeeUm5urt37rhyVys+TTz6pM2fO6LXXXtOIESP07rvvymKx2GqSrNc1BQQE5HlvxYr2/xfq6uqqChXyn8BR0DkwDMPuWI8++qiGDBmSb9/bb7/d9LMAAIoPoQkA4HC1a9fWpk2b1LZtW9OgExcXpzNnzmjlypVq3769bXtSUtINH/vVV1/V2bNntWDBAlWqVEkzZ8601SRJ1apVU5cuXW54/4VRtWpVeXt7Kycnp9iPBQC4flzTBABwuP79+ysnJ0dTp07N03bp0iWlpKRI+nPE5vIIjSRdvHhRc+fOvanjv/322+rbt6/eeOMNvfTSS5Kkrl27ysfHRy+//LKys7PzvOf06dM3dcwrOTk56aGHHtKKFSv03XffFeuxAADXj5EmAEChvffee/neX2j8+PE3td8OHTro8ccf1/Tp07V3717de++9cnZ2VmJiopYvX67Zs2erb9++atOmjSpVqqQhQ4Zo3Lhxslgs+uCDD+xC1I2oUKGCYmNjlZqaqueee07+/v4aPXq05s2bp0GDBqlZs2YaMGCAqlatqiNHjujTTz9V27ZtNWfOnJs67pVeeeUVbdmyRa1atdKIESPUoEEDnT17Vt988402bdqks2fPFtmxAADXh9AEACi0efPm5bs9Kirqpvc9f/58NW/eXG+//bb++c9/qmLFigoLC9Ojjz6qtm3bSpIqV66stWvXauLEiXr22WdVqVIlPfroo7r77rvVtWvXmzq+i4uLVq1apS5duuiJJ56Qn5+fHnnkEQUFBemVV17RjBkzlJWVpeDgYLVr167AVfluVPXq1bV79269+OKLWrlypebOnavKlSurYcOGevXVV4v0WACA62Mxbvaf5wAAAACgHOOaJgAAAAAwQWgCAAAAABOEJgAAAAAwQWgCAAAAABOEJgAAAAAwQWgCAAAAABOEJgAopywWi6ZMmVKsx+jYsaM6duxYZPu7uuYpU6bIYrHo999/L5b9x8TEyGKxKDk52bYtLCxM9913X5EcrygkJyfLYrEoJibG0aUAwC2L0ASg1Lv8R+Prr7/u6FIkWW/karFY5OPjowsXLuRpT0xMlMViKVU1O1pcXJwsFos++ugjR5eC/3fw4EH97W9/U5s2beTm5pYnPF7L7t27NXr0aDVv3lzOzs6yWCz59rscTAt6xMbG2vquWrVKXbt2VVBQkFxdXXXbbbepb9+++u6772724163zz//XMOHD1ejRo3k5OSksLCwAvtOmzZN999/v6pXr14i/1gBoORVdHQBAFAWVaxYURkZGfrkk0/Uv39/u7bY2Fi5ubkpMzPTQdVZXbhwQRUrlq3/mS/pmgcNGqQBAwbI1dW1xI5ZWuzcuVNvvfWWGjRooIiICO3du/e63r9u3TotWLBAt99+u2rVqqWffvop337t27fXBx98kGf7m2++qX379unuu++2bTtw4IAqVaqk8ePHq0qVKjpx4oTee+89tWzZUjt37lSTJk2uq8absWjRIi1dulTNmjVTUFCQad9nn31WAQEBatq0qTZs2FBCFQIoSWXr/00BoJRwdXVV27ZttXjx4jyhadGiRerZs6dWrFhRZMfLzMyUi4uLKlQo/AQBNze3Ijt+SSnpmp2cnOTk5FSixywt7r//fqWkpMjb21uvv/76dYemUaNG6amnnpK7u7vGjh1bYGiqVauWatWqZbftwoULGj16tDp37qyAgADb9ueffz7P+x977DHddtttmjdvnubPn39dNeYnOTlZNWvW1JYtW0ynlr788sv673//K2dnZ913332mo11JSUkKCwvT77//rqpVq950jQBKH6bnASg3Tp06peHDh6t69epyc3NTkyZNtHDhwjz9zpw5o0GDBsnHx0d+fn4aMmSI9u3bd93XjTzyyCP67LPPlJKSYtv21VdfKTExUY888kie/mfPntWkSZPUuHFjeXl5ycfHR927d9e+ffvs+l2eyrZkyRI9++yzCg4OloeHh9LS0iRJy5cvV4M
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADGXUlEQVR4nOzdd1xV5R8H8M9lXPYQRBmyBAWcuVNSNDX3z5EbBZylkFpYaqWZ5EhzUOYoFSxx5sjMIiVxkCPNWYioEA5cKCCyL+f3x4kLlyX7XODzfr3u63k483MvWHw5z3mOTBAEAURERERERKSkIXUAIiIiIiIidcNCiYiIiIiIqAAWSkRERERERAWwUCIiIiIiIiqAhRIREREREVEBLJSIiIiIiIgKYKFERERERERUAAslIiIiIiKiAlgoERERERERFcBCiYiIitS9e3d0795d6hiVIjg4GDKZDLGxsWXe18fHBw4ODpWeSSoV+SyIiOoSFkpERLVE7i/AuS9dXV00bdoUfn5+ePjwodTxahWZTAY/P79Cy5csWQKZTIaJEyciJydHgmRERFRZtKQOQERElWvRokVwdHREeno6Tp06hfXr1+Pw4cO4du0a9PX1pY4nifHjx2P06NHQ0dGpsnMsW7YMH330Eby9vbFp0yZoaPBvkURENRkLJSKiWqZfv35o3749AGDy5MkwNzfHqlWr8OOPP2LMmDESp5OGpqYmNDU1q+z4K1aswLx58+Dl5YUtW7ZUuEgSBAHp6enQ09OrpIRERFRW/HMXEVEt9/rrrwMAYmJiAADZ2dkICAiAk5MTdHR04ODggA8//BAZGRnFHiMlJQUGBgaYOXNmoXV3796FpqYmli5dCiBvCGBERATee+89WFhYwMDAAEOHDsXjx48L7b9u3To0b94cOjo6sLa2hq+vLxITE1W26d69O1q0aIErV67Aw8MD+vr6cHZ2xg8//AAAOH78ODp16gQ9PT24uLjg6NGjKvsXdV/Ojz/+iAEDBsDa2ho6OjpwcnJCQEAAFArFyz/UfFatWoUPPvgA48aNQ1BQkEqRlJOTgzVr1qB58+bQ1dVFw4YN8dZbb+HZs2cqx3BwcMDAgQMRGhqK9u3bQ09PDxs3bkR4eDhkMhl2796NxYsXo1GjRtDV1UXPnj1x8+bNQlnOnj2Lvn37wsTEBPr6+vDw8EBERESZ3g8REYlYKBER1XK3bt0CAJibmwMQrzItWLAAbdu2xerVq+Hh4YGlS5di9OjRxR7D0NAQQ4cOxa5duwoVEjt27IAgCPD09FRZ/s477+Dy5cv45JNPMG3aNPz000+F7utZuHAhfH19YW1tjZUrV+LNN9/Exo0b8cYbbyArK0tl22fPnmHgwIHo1KkTli9fDh0dHYwePRq7du3C6NGj0b9/fyxbtgwvXrzA8OHD8fz58xI/l+DgYBgaGuK9995DYGAg2rVrhwULFmDu3Lklf6D5BAYGwt/fH2PHjkVwcHChK0lvvfUW3n//fbi7uyMwMBATJkxASEgI+vTpU+j9RUVFYcyYMejduzcCAwPxyiuvKNctW7YM+/fvx+zZszFv3jycOXOm0Of9+++/o1u3bkhOTsYnn3yCJUuWIDExEa+//jrOnTtX6vdERET/EYiIqFYICgoSAAhHjx4VHj9+LNy5c0fYuXOnYG5uLujp6Ql3794VLl26JAAQJk+erLLv7NmzBQDC77//rlzm4eEheHh4KL8ODQ0VAAi//PKLyr6tWrVS2S43R69evYScnBzl8nfffVfQ1NQUEhMTBUEQhEePHglyuVx44403BIVCodxu7dq1AgBhy5YtKlkACNu3b1cuu379ugBA0NDQEM6cOVMoZ1BQUKFMMTExymWpqamFPsO33npL0NfXF9LT05XLvL29BXt7e5XtAAj29vYCAGHMmDFCdnZ2oWOdPHlSACCEhISoLP/1118LLc891q+//qqy7bFjxwQAgpubm5CRkaFcHhgYKAAQrl69KgiCIOTk5AhNmjQR+vTpo/KZp6amCo6OjkLv3r1L/CyIiKgwXlEiIqplevXqBQsLC9ja2mL06NEwNDTE/v37YWNjg8OHDwMA3nvvPZV9/P39AQA///xzice1trZGSEiIctm1a9dw5coVjBs3rtD2U6dOhUwmU37dtWtXKBQK/PvvvwCAo0ePIjMzE7NmzVK5EjNlyhQYGxsXymJoaKhy1cvFxQWmpqZwc3NDp06dlMtz+7dv3y72vQBQuf/n+fPnePLkCbp27YrU1FRcv369xH0BKGcSdHR0LPL+pz179sDExAS9e/fGkydPlK927drB0NAQx44dU9ne0dERffr0KfJcEyZMgFwuV37dtWtXlfd46dIlREdHY+zYsUhISFCe68WLF+jZsydOnDjBWfiIiMqIkzkQEdUyX3/9NZo2bQotLS00bNgQLi4uykLk33//hYaGBpydnVX2sbS0hKmpqbKIKYqGhgY8PT2xfv16pKamQl9fHyEhIdDV1cWIESMKbW9nZ6fydb169QBAeX9O7rlcXFxUtpPL5WjcuHGhLI0aNVIpvADAxMQEtra2hZblP09x/v77b3z88cf4/fffkZycrLIuKSmpxH0BwNvbG/fv38eSJUtQv359vPvuuyrro6OjkZSUhAYNGhS5/6NHj1S+dnR0LPZcL/sso6OjlZmKk5SUpNyPiIhejoUSEVEt07FjR+Wsd8UpWHCUlpeXF1asWIEDBw5gzJgx2L59OwYOHKgsTvIrbpY5QRDKde7ijlee8yQmJsLDwwPGxsZYtGgRnJycoKuri7/++gtz5swp1dUXLS0t7N69G3379oW/vz9MTU0xYcIE5fqcnBw0aNBA5QpcfhYWFipflzTD3cveY27eFStWqNzblJ+hoWGxxyciosJYKBER1SH29vbIyclBdHQ03NzclMsfPnyIxMRE2Nvbl7h/ixYt0KZNG4SEhKBRo0aIi4vDV199Ve4sgDiJQePGjZXLMzMzERMTg169epXruKURHh6OhIQE7Nu3D926dVMuz50ZsLR0dXVx8OBB9OjRA1OmTIGpqSmGDh0KAHBycsLRo0fh7u5e5dN8Ozk5AQCMjY2r9HMjIqpLeI8SEVEd0r9/fwDAmjVrVJavWrUKADBgwICXHmP8+PH47bffsGbNGpibm6Nfv37lytKrVy/I5XJ8+eWXKld/Nm/ejKSkpFJlKa/cKzT5z5uZmYl169aV+VjGxsb49ddf4ezsjDFjxiAsLAwAMHLkSCgUCgQEBBTaJzs7u9AU6BXRrl07ODk54YsvvkBKSkqh9UVNy05ERCXjFSUiojqkdevW8Pb2xjfffKMcfnbu3Dls3boVQ4YMQY8ePV56jLFjx+KDDz7A/v37MW3aNGhra5cri4WFBebNm4dPP/0Uffv2xf/+9z9ERUVh3bp16NChQ5ETRFSWLl26oF69evD29saMGTMgk8nw/fffl3tYoIWFBY4cOQJ3d3cMGTIEYWFh8PDwwFtvvYWlS5fi0qVLeOONN6CtrY3o6Gjs2bMHgYGBGD58eKW8Hw0NDWzatAn9+vVD8+bNMWHCBNjY2ODevXs4duwYjI2N8dNPP1XKuYiI6goWSkREdcymTZvQuHFjBAcHY//+/bC0tMS8efPwySeflGr/hg0b4o033sDhw4cxfvz4CmVZuHAhLCwssHbtWrz77rswMzPD1KlTsWTJknIXYKVhbm6OQ4cOwd/fHx9//DHq1auHcePGoWfPnsXOPPcytra2+O2339C1a1f069cPJ06cwIYNG9CuXTts3LgRH374IbS0tODg4IBx48bB3d29Ut9T9+7dcfr0aQQEBGDt2rVISUmBpaUlOnXqhLfeeqtSz0VEVBfIhPL++YyIiOqsoUOH4urVq7h586bUUYiIiKoE71EiIqIyiY+Px88//1zhq0lERETqjEPviIioVGJiYhAREYFNmzZBW1ubw7mIiKhW4xUlIiIqlePHj2P8+PGIiYnB1q1bYWlpKXUkIiKiKsN7lIiIiIiIiArgFSUiIiIiIqICWCgREREREREVwEKJiKiWkclkWLhwYZWeo3v37ujevXu
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADMzElEQVR4nOzdeVxV1frH8c8GZJJJFEUMQcUxret81RwqyynTzNSilDTtOvtLK60s09TMNC2vWllYSQ7lkNcsMxOna5PlUJGhQVjikAYKMrN/f5w4VwYVFdgM3/frdV77nL3XWfs5g3Ue1rPWNkzTNBERERERERE7B6sDEBERERERKW2UKImIiIiIiOShRElERERERCQPJUoiIiIiIiJ5KFESERERERHJQ4mSiIiIiIhIHkqURERERERE8lCiJCIiIiIikocSJRERERERkTyUKImISIno0qULXbp0KbL+YmNjMQyD5cuX2/dNmzYNwzCK7BxlQUHvg4iIXD8lSiIiFdjy5csxDMN+c3V1pUGDBowZM4aTJ09aHV6J6dKlC02bNs23f9u2bbi7u9OiRQvOnj1rQWQiImIVJ6sDEBER602fPp06deqQmprK7t27WbJkCZs3b+aHH37A3d29SM7x2WefFUk/l/PMM88wefLkIunriy++oHfv3jRs2JDPP/8cX1/fIulXRETKBiVKIiJCjx49aNWqFQCPPPIIVatWZf78+Xz00Ufcf//919X3hQsXcHd3x9nZuShCvSwnJyecnK7/f207duygd+/eNGjQoMiSpOTkZCpXrnzd/YiISMlQ6Z2IiORz2223ARATE2Pft2LFClq2bImbmxu+vr4MGjSIY8eO5XpeTgnbvn376NSpE+7u7jz11FP2Y3nnKJ06dYphw4ZRo0YNXF1dufnmm3nnnXfyxZOQkEBYWBje3t74+PgwZMgQEhIS8rW71BylFStW0KZNG9zd3alSpQqdOnW65AjXrl276NWrFyEhIXz++edUrVo11/FPPvmEjh07UrlyZTw9PenVqxc//vhjrjZhYWF4eHhw9OhRevbsiaenJ6GhoQAYhsGYMWPYsGEDTZs2xcXFhRtvvJFPP/00Xyx//PEHQ4cOpUaNGvZ2b7/9doFxi4hI0dKIkoiI5HP06FEAe5Iwc+ZMpk6dyoABA3jkkUc4ffo0r732Gp06deL777/Hx8fH/twzZ87Qo0cPBg0axIMPPkiNGjUKPEdKSgpdunThyJEjjBkzhjp16vDBBx8QFhZGQkIC48ePB8A0Tfr06cPu3bv517/+RePGjVm/fj1Dhgwp1Gt5/vnnmTZtGu3bt2f69Ok4Ozvz1Vdf8cUXX3DnnXfmartnzx569uxJnTp12LZtG9WqVct1/L333mPIkCF069aNOXPmcOHCBZYsWcItt9zC999/T3BwsL1tZmYm3bp145ZbbuHll1/OVcK4e/du1q1bx6hRo/D09OTVV1/l3nvvJS4uzv6enzx5kn/+85/2xMrPz49PPvmEYcOGce7cOSZMmFCo1y8iItfIFBGRCis8PNwEzM8//9w8ffq0eezYMXPVqlVm1apVTTc3N/P33383Y2NjTUdHR3PmzJm5nnvo0CHTyckp1/7OnTubgLl06dJ85+rcubPZuXNn++MFCxaYgLlixQr7vvT0dLNdu3amh4eHee7cOdM0TXPDhg0mYL700kv2dpmZmWbHjh1NwAwPD7fvf+6558yL/9cWHR1tOjg4mPfcc4+ZlZWVK57s7Oxcsfn6+pqenp7mjTfeaJ46dSpf/OfPnzd9fHzM4cOH59p/4sQJ09vbO9f+IUOGmIA5efLkfP0AprOzs3nkyBH7vgMHDpiA+dprr9n3DRs2zKxZs6b5559/5nr+oEGDTG9vb/PChQumaZpmTExMvvdBRESun0rvRESErl274ufnR2BgIIMGDcLDw4P169dTq1Yt1q1bR3Z2NgMGDODPP/+03/z9/alfvz7bt2/P1ZeLiwsPP/zwFc+5efNm/P39c82BqlSpEuPGjSMpKYkdO3bY2zk5OTFy5Eh7O0dHR8aOHXvFc2zYsIHs7GyeffZZHBxy/y8vb4lecnIy58+fp0aNGnh5eeXra+vWrSQkJHD//ffneh8cHR1p27ZtvvcByBXzxbp27Uq9evXsj2+66Sa8vLz49ddfAdso2tq1a+nduzemaeY6X7du3UhMTOS777674usXEZFrp9I7ERHh3//+Nw0aNMDJyYkaNWrQsGFDe2IRHR2NaZrUr1+/wOdWqlQp1+NatWoVauGG3377jfr16+dLYBo3bmw/nrOtWbMmHh4eudo1bNjwiuc4evQoDg4ONGnS5IptQ0JCGDx4ME8++ST3338/H3zwAY6Ojvbj0dHRwP/mb+WVN7lycnLihhtuKLBt7dq18+2rUqUKf/31FwCnT58mISGBN954gzfeeKPAPk6dOnXF1yQiItdOiZKIiNCmTRv7qnd5ZWdnYxgGn3zySa7EIUfeBMbNza1YYiwJTzzxBGfOnOGll15i+PDhvPXWW/aRp+zsbMA2T8nf3z/fc/Outufi4pIvCcxR0PsItpGki8/14IMPXnIu1k033VSIVyQiItdKiZKIiFxWvXr1ME2TOnXq0KBBgyLrNygoiIMHD5KdnZ0rofj555/tx3O227ZtIykpKVdSdvjw4ULFnp2dzU8//cQ//vGPQsU1Z84czp49y7Jly6hSpQrz5s2z9wVQvXp1unbtWqi+rpWfnx+enp5kZWUV+7lERKRgmqMkIiKX1a9fPxwdHXn++eftIx45TNPkzJkz19Rvz549OXHiBKtXr7bvy8zM5LXXXsPDw4POnTvb22VmZrJkyRJ7u6ysLF577bUrnqNv3744ODgwffp0+yjNxbFfyuuvv07//v2ZP38+L7zwAgDdunXDy8uLWbNmkZGRke85p0+fvmI8heXo6Mi9997L2rVr+eGHH4r1XCIiUjCNKImIyGXVq1ePF154gSlTphAbG0vfvn3x9PQkJiaG9evXM2LECCZNmnTV/Y4YMYLXX3+dsLAw9u3bR3BwMB9++CF79uxhwYIFeHp6AtC7d286dOjA5MmTiY2NpUmTJqxbt47ExMQrniMkJISnn36aGTNm0LFjR/r164eLiwvffPMNAQEBzJ49u8DnOTg4EBERQWJiIlOnTsXX15dRo0axZMkSHnroIVq0aMGgQYPw8/MjLi6Ojz/+mA4dOrBo0aKrfh8u5cUXX2T79u20bduW4cOH06RJE86ePct3333H559/ztmzZ4vsXCIikp8SJRERuaLJkyfToEEDXnnlFZ5//nkAAgMDufPOO7n77ruvqU83NzciIyOZPHky77zzDufOnaNhw4aEh4cTFhZmb+fg4MDGjRuZMGECK1aswDAM7r77bubNm0fz5s2veJ7p06dTp04dXnvtNZ5++mnc3d256aabeOihhy77PGdnZ9avX0/Xrl0ZO3YsPj4+PPDAAwQEBPDiiy8yd+5c0tLSqFWrFh07dizUSn9Xo0aNGnz99ddMnz6ddevWsXjxYqpWrcqNN97InDlzivRcIiKSn2FervZARERERESkAtIcJRERERERkTyUKImIiIiIiOShRElERERERCQPJUoiIiIiIiJ5KFESERERERHJQ4mSiIiIiIhIHkqURETKGcMwmDZtWrGeo0uXLnTp0qXI+ssb87Rp0zAMgz///LNY+l++fDmGYRAbG2vfFxwczF133VUk5ysKsbGxGIbB8uXLrQ5FRKRCUqIkIqVWzg/Fl19+2epQAAgLC8MwDLy8vEhJScl3PDo6GsMwSlXMVouMjMQwDD788EOrQ5G/rVu3joEDB1K3bl3c3d1p2LAhEydOJCEhodB9ZGdns2TJEv7xj3/g5uZG1apVue222zhw4ECudvHx8YwYMYI6derg5uZGvXr1eOyxxzhz5ky+PqOioujevTseHh74+vry0EMPcfr06et9udfkrbfeonHjxri6ulK/fn1ee+21AtutWrWKFi1a4Orqip+fH8OGDSuy5F5ErOdkdQAiImWJk5MTFy5c4D//+Q8DBgzIdSwiIgJXV1d
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADJB0lEQVR4nOzde3zO9f/H8cdl5/PMaTabYQ4T+jrnECrl0IEkqoVFfH8hKSoq5RBKKfpKKmWVc4kkJckcllQiqqWJRc45jNn52uf3x6frymwYNp8dnvfbbbf3tc/nc30+r+vapfba+/16fWyGYRiIiIiIiIiIUzmrAxARERERESlulCiJiIiIiIicQ4mSiIiIiIjIOZQoiYiIiIiInEOJkoiIiIiIyDmUKImIiIiIiJxDiZKIiIiIiMg5lCiJiIiIiIicQ4mSiIiIiIjIOZQoiYiIWCQuLg6bzUZcXJzVoYiIyDmUKImIyGWLjY3FZrM5v1xdXQkNDSUmJob9+/fnOb5Dhw65jnd3d6dGjRoMGjSIffv2XfDcZ3+NGjXKeVxERAS33XZbnmt98MEHuLi40LlzZ9LT0wv/xYuISKnmanUAIiJS8o0fP54aNWqQnp7Ot99+S2xsLBs3buTnn3/G09Mz17HVqlVj8uTJAGRmZvLrr78ya9YsVq1aRUJCAt7e3vme+2wNGjS4YDzz5s0jJiaGjh07smzZsjwxiIiIXIwSJRERuWJdunShWbNmADz44INUrFiRF198keXLl9OrV69cxwYEBHD//ffn2lajRg2GDh1KfHw8N99883nPXRALFy6kX79+3HjjjXzyySeFkiSlpqbmSeBERKR009I7EREpdNdffz0Af/zxR4GODw4OBsDV9cr+frd48WLuv/9+OnTowPLly/MkSXPnzqVp06Z4eXkRFBTEPffck2fJX4cOHWjQoAFbtmyhXbt2eHt789RTT5GUlITNZuPll1/mrbfeolatWnh4eNC8eXO+//77PLH89ttv9OzZk6CgIDw9PWnWrBnLly+/otcnIiJXj2aURESk0CUlJQFQvnz5PPvsdjt///03AFlZWSQkJPDcc88RGRlJmzZt8hyfnJzsPN6hYsWKeY5bsmQJ0dHRtGvXjk8//RQvL69c+ydOnMiYMWPo1asXDz74IEePHuV///sf7dq1Y+vWrQQGBjqPPXbsGF26dOGee+7h/vvvp0qVKs598+fP5/Tp0/z3v//FZrMxZcoUevTowe7du3FzcwPgl19+oU2bNoSGhjJq1Ch8fHxYvHgx3bt3Z8mSJdx5550FeyNFRMQySpREROSKOZKZ9PR0Nm/ezLhx4/Dw8Mi3ycJvv/1GpUqVcm2Lioriyy+/xN3dPc/xHTt2zLPNMIxc32/dupVVq1bRtm1bVqxYkSdJ+vPPP3nuued4/vnneeqpp5zbe/ToQePGjZk5c2au7YcOHWLWrFn897//dW5zJH979+4lMTHRmQTWrVuXbt26sWrVKufrfeSRRwgPD+f777/Hw8MDgMGDB9O2bVuefPJJJUoiIiWAEiUREbli5yYzERERzJ07l2rVquU5NiIigrfffhuA7Oxsdu7cyZQpU+jSpQsbNmzIk0S9/vrr1KlT54LXP378ONnZ2VSrVi1PkgTw8ccfk5OTQ69evXLNTgUHB1O7dm3Wrl2bK1Hy8PDggQceyPdavXv3zjVT5lhmuHv3bmcsX3/9NePHj+f06dOcPn3aeWynTp147rnn2L9/P6GhoRd8TSIiYi0lSiIicsUcyUxycjLvvvsu69evd86knMvHxydXYtW5c2fatm1Ls2bNeOGFF5g6dWqu41u0aHHRZg433XQT4eHhvPHGGwQFBTF9+vRc+xMTEzEMg9q1a+f7fMeSOYfQ0NB8Z7cAwsPDc33vSJpOnDgBwK5duzAMgzFjxjBmzJh8z3HkyBElSiIixZwSJRERuWJnJzPdu3enbdu23HfffezcuRNfX9+LPr9p06YEBASwfv36y45hxowZnDhxgtdee43y5cszduxY576cnBxsNhuff/45Li4ueZ57boz5zUo55Pd8+Hc5YE5ODgAjR46kU6dO+R4bGRl5wdciIiLWU6IkIiKFysXFhcmTJ3PDDTcwY8aMXDeHvRC73U5KSsplX7dcuXK8//77JCcnM27cOIKCghg2bBgAtWrVwjAMatSocdFlfFeqZs2agDlLlV99lYiIlAxqDy4iIoWuQ4cOtGjRgmnTppGenn7R49euXUtKSgrXXnvtFV3Xzc2Njz76iDZt2jB8+HA++OADwGza4OLiwrhx4/I0gjAMg2PHjl3Rdc9WuXJlOnTowJtvvsnBgwfz7D969GihXUtERIqOZpRERKRIPP7449x9993Exsbyf//3f87tycnJzJ07F/i3mcMbb7yBl5dXgWefLsTb25vPPvuM9u3b079/fwICArjjjjt4/vnnGT16NElJSXTv3h0/Pz/27NnD0qVLGTRoECNHjrziazu8/vrrtG3bloYNGzJw4EBq1qzJ4cOH2bRpE3/99Rc//fRToV1LRESKhhIlEREpEj169KBWrVq8/PLLDBw40Fnb89dff9GnTx8AbDYb5cuXp3379jz33HP85z//KZRrBwQEONuF9+7dm88//5xRo0ZRp04dXn31VcaNGwdAWFgYt9xyC3fccUehXNehfv36/PDDD4wbN47Y2FiOHTtG5cqVady4Mc8++2yhXktERIqGzTh3DYKIiIiIiEgZpxolERERERGRcyhREhEREREROYcSJRERERERkXMoURIRERERETmHEiUREREREZFzKFESERERERE5hxIlEZFSxmazMXbs2CK9RocOHejQoUOhne/cmMeOHYvNZuPvv/8ukvPHxsZis9lISkpybouIiOC2224rlOsVhqSkJGw2G7GxsVaHIiJSJilREpFiy/GL4ssvv2x1KADExMRgs9nw9/cnLS0tz/7ExERsNluxitlqcXFx2Gw2PvroI6tDkX9EREQ4P6fnftWuXbtA58jMzGTSpEnUq1cPT09PqlSpwq233spff/113udMnDgRm81GgwYN8uzr0KFDvvF07tz5sl/n5crJyWHKlCnUqFEDT09PGjVqxIIFC/IcExsbyx133EFYWBg+Pj40aNCA559/nvT09Kses4gUDVerAxARKUlcXV1JTU3l008/pVevXrn2zZs3D09PT8t/UUpLS8PVtWT95/1qx9ynTx/uuecePDw8rto1i4tp06aRkpKSa9uff/7JM888wy233HLR52dlZXHrrbfyzTffMHDgQBo1asSJEyfYvHkzycnJVKtWLc9z/vrrLyZNmoSPj895z1utWjUmT56ca1tISEgBX1Xhefrpp3nhhRcYOHAgzZs355NPPuG+++7DZrNxzz33AJCamsoDDzzAddddx//93/9RuXJlNm3axHPPPceaNWv4+uuvsdlsVz12ESlcJev/pCIiFvPw8KBNmzYsWLAgT6I0f/58br31VpYsWVJo10tPT8fd3Z1y5Qq+AMDT07PQrn+1XO2YXVxccHFxuarXLC66d++eZ9vzzz8PQHR09EWf/+qrr7Ju3To2btxIixYtCnTNkSNHct1112G328+7nDIgIID777+/QOe7HGPHjiU2NjbXcstz7d+/n6lTpzJkyBBmzJgBwIMPPkj79u15/PHHufvuu3FxccHd3Z34+Hhat27tfO7AgQOJiIhwJksdO3YsstciIleHlt6JSIl35MgRBgwYQJUqVfD09OTaa6/lvffey3PcsWPH6NOnD/7+/gQGBtKvXz9++umnS64Due+++/j88885efKkc9v3339PYmIi9913X57jjx8/zsiRI2nYsCG+vr74+/vTpUsXfvrpp1zHOZapLVy4kGeeeYbQ0FC8vb05deoUAB9++CH169fH09OTBg0asHTpUmJiYoiIiMh1nvPV++zatYuYmBgCAwMJCAjggQceIDU1Nddz58yZw4033kjlypXx8PCgfv36vPHGGwV+by5XQeqq/vzzTyIjI2nQoAGHDx8G4OTJkwwfPpywsDA8PDyIjIzkxRdfJCcn54Lnyq9
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gp_linear.plot(X_3D_star, y_star_linear, var_linear)\n",
    "gp_poly.plot(X_3D_star, y_star_poly, var_poly)\n",
    "gp_periodic.plot(X_3D_star, y_star_periodic, var_periodic)\n",
    "gp_rbf.plot(X_3D_star, y_star_rbf, var_rbf)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compare the results with the results of the scikit-learn implementation of Gaussian Processes"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "outputs": [],
   "source": [
    "#import the scikit-learn implementation of Gaussian Processes\n",
    "# TODO\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/liger/PycharmProjects/ML_course/venv/lib/python3.8/site-packages/sklearn/gaussian_process/kernels.py:420: ConvergenceWarning: The optimal value found for dimension 0 of parameter length_scale is close to the specified lower bound 1e-05. Decreasing the bound and calling fit again may find a better value.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# Fit the sklearn model\n",
    "# TODO\n",
    "\n",
    "# Fit the custom model\n",
    "gp_rbf.fit(X_1D, y_1D)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Predict the values of the test data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "outputs": [],
   "source": [
    "# Predict the mean and variance of the test data using the scikit-learn implementation\n",
    "y_star, var = None, None # TODO\n",
    "# Predict the mean and variance of the test data using the custom implementation\n",
    "y_star_rbf, var_rbf = gp_rbf.predict(X_1D_star)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plot the results comparing with the results of the scikit-learn with the ccustom implementation of Gaussian Processes using the RBF kernel"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "outputs": [
    {
     "data": {
      "text/plain": "<matplotlib.legend.Legend at 0x7ff0dc579580>"
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x600 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAAINCAYAAADMabVmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACa7UlEQVR4nOzdd3gU5d6H8Xt20yEJJKRSA4gQei+CdMGCYAcLVlDsYi9HRN9jL6jH3lCxF0BQUIpgoYTeO6EnBAikkrYz7x/RHPWAUpI8u5vv57pyZXczu7nDpUl+mZlnLMdxHERERERERKRcuEwHiIiIiIiI+BMNWSIiIiIiIuVIQ5aIiIiIiEg50pAlIiIiIiJSjjRkiYiIiIiIlCMNWSIiIiIiIuVIQ5aIiIiIiEg50pAlIiIiIiJSjgJMB3g727bZs2cP4eHhWJZlOkdERERERAxxHIecnBwSExNxuY6+v0pD1j/Ys2cPdevWNZ0hIiIiIiJeYufOndSpU+eoH9eQ9Q/Cw8OB0n/IiIgIwzUiIiIiImJKdnY2devWLZsRjkZD1j/4/RDBiIgIDVkiIiIiIvKPpxFp4QsREREREZFypCFLRERERESkHGnIEhERERERKUcaskRERERERMqRhiwREREREZFypCFLRERERESkHGnIEhERERERKUcaskRERERERMqRhiwREREREZFypCFLRERERESkHGnIEhERERERKUcaskRERERERMqRhiwREREREZFyFGA6QERERESOn8d2SEnNJCOngNjwEDolReF2WaazRAQNWSIiIiI+Z/rqNMZOWUtaVkHZYwmRIYwZlMzAFgkGy0QEdLigiIiIiE+ZvjqNUROW/mnAAkjPKmDUhKVMX51mqExEfqchS0RERMRHeGyHsVPW4vx2PzFoA6eEpACUPTZ2ylo8tnPE54tI5dDhgiIiIiI+YsGWDGoUzeKUuIXsr76XnUEWOUDLHbmsyuuDA6RlFZCSmknXRtGmc0WqLA1ZIiIiIt6sKB+2zoEN3zE+fRa7GgSwC4D/LnJRI3I+5PUpu5+RU/DXVxGRSqQhS0RERMTL7N+3jrnL32Zu2nwe37mV6sWHAWhXI4JVQRE0zKuGJ6cZthPA+top7KiWhYsS7N9+tYsNDzGZL1LlacgSERERMcyxbTanzmDOmo+Zc2AlK10lZR+bFwhnhNWDpmdxWVJvPvq6hPnZFg4QQCH/CprNwMNZ3GhtYalzKvGRpcu5i4g5GrJERERETPAUw475LFk1gYf2z2eX+7fHf1uWrIUdQK+oFjTrcRUk9QHLIgJ46NzS1QUtoIRgIrKaEef+lQHuxSwtOZUxg5J1vSwRwzRkiYiIiFSSnOzd/Lr8baLSVtMpdSEUZBEf4GZX3doEOQ5dXOH0iu9Cz9bXEhvX4oivMbBFAq9d3q7sOlnfezpynvtXzgpYQr1LntV1skS8gIYsERERkQq0e3cKc1aOZ07GYhY7+ZRYFj3zD9OpIAvCalG7yUDeiG9Em1bDCQurdUyvObBFAv2T40lJzeTAwVP4au57/BDm4e6w1UBixX5BIvKPNGSJiIiI/AOP7ZCSmklGTgGx4aXnPB31kDzbxtmzlNcXPcfMQ2vZ6LL/+zHLoqHHIjm2LZx9F9TpAC433U6gye2yflumPZqb19ZmnpPL7NUf0rjRGSfyJYpIOdKQJSIiIvI3pq9OKzs073cJkSGMGZRcdmheYUEWa1d/Stu0tbBhOlZuOinxsWwMDcHlOLQlhN4xbenV4grq1z+93Bv7JJ7G3N3fMytzFSPL/dVF5HhpyBIRERE5iumrSxeZcP7yeHpWAQ988i2b265g8+GlzCvJosCCOTt2E2XbEFSdqyKTOT++BT3aXEeNmkkV2tmr3Q24dk1nrctDetoy4hPaVujnE5G/pyFLRERE5Ag8tsPYKWv/MGA5NAhaT92IOeRW38GWEJs38n47ZNBlEedx2NHqfKJaDIMGPegZEFxprVFRjWlDMEspYvaKt7k04ZVK+9wi8r80ZImIiIgcQUpq5p8OEXw64E0CIxbzcEz0b49YJBVCt/CmDG53OU1PGYTlcpmJBfrEtGfp/vnM3pvCpcYqRATKrsQgIiIiIn+UkfPfAauHayUXB8yle14ByXkBtNp7KmGbrmfl1idp1mwczU4dbHTAAujT6ioAFjuHyTq0zWiLSFWnPVkiIiIiRxAbHgKAixIeCPgIgG+Kz2DhjiuOuJ1pdet2o5ntIrIwj8wN3xLZ+SbTSSJVlvZk+QiP7TB/ywEmL9/N/C0H8Nh/PQVXREREylOnpCgSIkPoEvUxTyUWsjAwgpdKziv7uEXpKoOdkqLMRf7Fx/Uv4q30fSRtm286RaRK054sH3AsS8eKiIhI+XK7LB7oH8MzK1ZzMCCEl4Oak0V1oHTAAhgzKPno18syIKDZIPjpGdg8C4oPQ2Co6SSRKkl7srzc70vH/nHAgtKlY0dNWMr01WmGykRERPxf6o7/42CAi8Rih4UH/nuYYHxkCK9d3s77/tgZ3woi67LPLmDXmq9M14hUWdqT5cX+unRsKAXYuCgkCIfSv6KNnbKW/snxXvVXNBEREX+QnraMD7LWgsvirmaXUb13LzJyCogNLz1E0Ct/9loW79drxrO5FkPWjuexNpebLhKpkjRkebE/Lx3r8FTgW2RX282r+Vezp+hUHCAtq4CU1Ey6Nor+u5cSERGR4/TS7NEUuizaO8H063qv8dUDj1WzhgNg5XrmFKZTUlxAQKB3LMwhUpX4xneLKuqPS8cmkEnt4HW8FGfjNHiHduFTj7idiIiInLw1675iSsl+AO7ufL/PDFgA7VpcTqTtcMhlsWzNx6ZzRKok3/mOUQX9cUnYNKIZVXAPiUVuct0uNtX5hdPjxuGixGuWjhUREfELjsOnC54GYFBADM2bXWA46PgEBIbQMzgOgNmbJpmNEamiNGR5sd+Xjv39iO+MkgasSR1D20M1AVgWlU7HBo/QODLTXKSIiIi/WT+VMdvX81BmDrf2ec50zQnp0+AMAH7MScWxbcM1IlWPhiwv5nZZjBmUDPx3qdhiQvgp7V6a72lLiO2wNrSES785j5VrPjMXKiIi4i9KimDGwwQAl7S+jviEtqaLTki3NiMIsR12u2Hj5u9M54hUORqyvNzAFgm8dnk74iP/fEjgdq7krvoP0sAD6W6LmbPug5S3wNFFikVERE7Ump+foChzK1SLhdNuM51zwkLDougaEAnArLU6L0uksml1QR8wsEUC/ZPjSUnN/J+lY8/O6cWH31zBiMwd8N1dsGsRnDMOgsJMZ4uIiPiUrEPbGLntC8LrJPJW8xuoGxxuOumkXJ50Dv0XvMTp+dtMp4hUOdqT5SPcLouujaIZ3KY2XRtFl12bo3p4AqMunUFA/8fAclO88jPGfNCd7dt/NlwsIiLiW16fcSvZLoswdxCJHUeZzjlpndpdz6C8QiL3roGD20zniFQpGrL8gWXBabfCld/wWmwiXwcWM3T2KGbNe8p0mYiIiE/Yvv1nPs3bCsBdLa/HHRBkuKgchEVB/W6lt9frvCyRyqQhy5806M6wC7+gnRNErsvi9k0TeP6r8ykp1nW0RERE/s4LP91HiWXRnWp063Cj6Zxyk9m4L+9GhvP4hg9Np4hUKRqy/ExMbAvevuxnhoc1BOC93E2M/Kg7+/evN1wmIiLinRYvH88sOxu343BX98dM55Srww2780JUTT4jh8zMzaZzRKoMDVl+KDAwjLsvmsyzDS8mzHZYZBVy8TcXsmr1J6bTREREvIrtKeGZZeMAuCCkDo0a9TcbVM5qJ3akme3GtizmLnvTdI5IlaEhy48N6PEvPun1Eg09Focth/BvboMFr2mZdxERkd/kLv+IWgW5VLMdbuz3kumcCtE7qgUAs3f/YrhEpOrQkOXnGib14ZOLZ/J6aDINigph+n3w5TV4DmeZThMRETGrKJ+IuU/xyt59TGpwCdG1mpguqhB9ml8KwHxPNvn5+w3XiFQNPjVk/fTTTwwaNIjExEQsy2LSpEn/+Jw5c+bQrl07goODady4MePHj6/wTm8
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(X_1D, y_1D, \"o\", label=\"data\")\n",
    "\n",
    "plt.plot(X_1D_star, y_star, \"-\", label=\"prediction(scikit-learn)\")\n",
    "plt.plot(X_1D_star, y_star_rbf, \"--\", label=\"prediction(custom)\")\n",
    "\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.legend(loc=\"best\")\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Bonus: Use the Gaussian Process to fit any dataset of your choice"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}