ML_course/assignment 6/iml_assgnment6_solved.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Solution for Assignment 5 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": 5,
   "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": 6,
   "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": 7,
   "outputs": [],
   "source": [
    "class GaussianProcess:\n",
    "    def __init__(self, kernel, sigma_n=1e-8):\n",
    "        self.kernel = kernel\n",
    "        self.sigma_n = sigma_n\n",
    "        self.X = None\n",
    "        self.y = None\n",
    "        self.K = None\n",
    "        self.K_inv = None\n",
    "        self.alpha = None # coefficients alpha of the Gaussian process for\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "        # Compute the kernel matrix. Add a small value to the diagonal for numerical stability.\n",
    "        self.K = self.kernel(self.X, self.X) + self.sigma_n * np.eye(len(self.X))\n",
    "        # Compute the inverse of the kernel matrix.\n",
    "        self.K_inv = np.linalg.inv(self.K)\n",
    "        # Compute the coefficients alpha.\n",
    "        self.alpha = np.dot(self.K_inv, self.y)\n",
    "\n",
    "    def predict(self, X_star):\n",
    "        # Compute the kernel matrix between the training data and the test data.\n",
    "        K_star = self.kernel(self.X, X_star)\n",
    "        # Compute the predictive mean.\n",
    "        y_star = np.dot(K_star.T, self.alpha)\n",
    "        # Compute the predictive variance.\n",
    "        v = np.dot(self.K_inv, K_star)\n",
    "        var = self.kernel(X_star, X_star) - np.dot(K_star.T, v)\n",
    "        return y_star, var\n",
    "\n",
    "    def log_marginal_likelihood(self):\n",
    "        # Compute the log marginal likelihood.\n",
    "        return (\n",
    "            -0.5 * np.dot(self.y.T, self.alpha)\n",
    "            - 0.5 * np.linalg.slogdet(self.K)[1]\n",
    "            - 0.5 * len(self.X) * np.log(2 * np.pi)\n",
    "        )\n",
    "\n",
    "    def gradient(self):\n",
    "        # Compute the gradient of the log marginal likelihood with respect to the kernel hyperparameters.\n",
    "        # Compute the gradient of the kernel matrix with respect to the kernel hyperparameters.\n",
    "        K_grad = self.kernel.gradient(self.X)\n",
    "        # Compute the gradient of the log marginal likelihood.\n",
    "        return (\n",
    "            0.5\n",
    "            * np.einsum(\"ijk,ij->k\", K_grad, np.dot(self.K_inv, np.dot(self.y, self.y.T)))\n",
    "            - 0.5 * np.einsum(\"ij,ijk->k\", self.K_inv, K_grad)\n",
    "        )\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": 8,
   "outputs": [],
   "source": [
    "#Implement the linear kernel function\n",
    "class LinearKernel:\n",
    "    def __init__(self, theta=1.0):\n",
    "        self.theta = theta\n",
    "        self.bounds = ((1e-5, 1e5),)\n",
    "        self.num_params = 1\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        return self.theta * np.dot(X1, X2.T)\n",
    "\n",
    "    def set_params(self, params):\n",
    "        self.theta = params[0]\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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",
    "        self.bounds = ((1e-5, 1e5), (1, 10))\n",
    "        self.num_params = 2\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        return (self.theta * np.dot(X1, X2.T) + 1) ** self.d\n",
    "\n",
    "    def set_params(self, params):\n",
    "        self.theta = params[0]\n",
    "        self.d = params[1]\n",
    "\n",
    "    def gradient(self, X):\n",
    "        return np.array(\n",
    "            [\n",
    "                self.d * self.theta * np.dot(X, X.T) ** (self.d - 1),\n",
    "                self.theta * np.dot(X, X.T) ** self.d * np.log(np.dot(X, X.T)),\n",
    "            ]\n",
    "        )"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [],
   "source": [
    "#Implement the periodic kernel function\n",
    "class PeriodicKernel:\n",
    "    def __init__(self, theta=1.0):\n",
    "        self.theta = theta\n",
    "        self.bounds = ((1e-5, 1e5),)\n",
    "        self.num_params = 1\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        return np.exp(-2 * np.sin(np.pi * cdist(X1, X2) / self.theta) ** 2)\n",
    "\n",
    "    def set_params(self, params):\n",
    "        self.theta = params[0]\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [],
   "source": [
    "#Implement the radial basis function kernel function\n",
    "class RBFKernel:\n",
    "    def __init__(self, theta=1.0):\n",
    "        self.theta = theta\n",
    "        self.bounds = ((1e-5, 1e5),)\n",
    "        self.num_params = 1\n",
    "\n",
    "    def __call__(self, X1, X2):\n",
    "        return np.exp(-self.theta * cdist(X1, X2) ** 2)\n",
    "\n",
    "    def set_params(self, params):\n",
    "        self.theta = params[0]\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": 12,
   "outputs": [],
   "source": [
    "#Linear kernel\n",
    "kernel_linear = LinearKernel(theta=1.0)\n",
    "#Polynomial kernel\n",
    "kernel_poly = PolynomialKernel(theta=1.0, d=2)\n",
    "#Periodic kernel\n",
    "kernel_periodic = PeriodicKernel(theta=1.0)\n",
    "#Radial basis function kernel\n",
    "kernel_rbf = RBFKernel(theta=1.0)"
   ],
   "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": 13,
   "outputs": [],
   "source": [
    "gp_linear = GaussianProcess(kernel=kernel_linear, sigma_n=1e-8)\n",
    "gp_poly = GaussianProcess(kernel=kernel_poly, sigma_n=1e-8)\n",
    "gp_periodic = GaussianProcess(kernel=kernel_periodic, sigma_n=1e-8)\n",
    "gp_rbf = GaussianProcess(kernel=kernel_rbf, sigma_n=1e-8)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Fit the 1D data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 15,
   "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": 16,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_12816/3712199621.py:45: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
      "/tmp/ipykernel_12816/3712199621.py:46: 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/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACKVElEQVR4nOzdd3hUZd7G8XsmZdILNQGBAAmE0FtCKGJhRVEXV1ksKASCunbXjq/uuqKioi6urqIrZVexoCLr2lhlRem9l5hAIiA9kN5nzvvHmJFhEgxwkkn5fq5rriHnPOec30xG4Z6nHIthGIYAAAAAAKaxersAAAAAAGhsCFoAAAAAYDKCFgAAAACYjKAFAAAAACYjaAEAAACAyQhaAAAAAGAyghYAAAAAmIygBQAAAAAmI2gBAAAAgMkIWgAAr8nKypLFYtHcuXO9XUqjxvsMAHWPoAUAqBVz586VxWLRunXrvF2KKS644AL16NHDY/vixYsVFBSkfv366fjx416oDABQH/l6uwAAQNPVoUMHFRcXy8/Pz9ulnJX//e9/uvLKK9W1a1d98803atasmbdLAgDUE/RoAQC8xmKxKCAgQD4+Pt4upVpFRUVVbv/uu+905ZVXqkuXLqaFrMLCwnM+BwCgfiBoAQC8pqq5QykpKQoJCdFPP/2kq666SiEhIWrZsqUeeOAB2e12t+MdDodmzJih7t27KyAgQK1bt9att96qEydOuLX797//rcsvv1xt2rSRzWZT586dNXXqVI/zVQ4PXL9+vc4//3wFBQXp0Ucf9ah76dKluvzyyxUbG6tvvvlGzZs3d9v/5ZdfatiwYQoODlZoaKguv/xybd++3a1N5evcvXu3Ro0apdDQUI0bN06SM4DeeeedWrhwoXr06CGbzabu3bvrq6++8qjlp59+0qRJk9S6dWtXu9mzZ//6mw8AqFUMHQQA1Dt2u10jR45UUlKSXnjhBX3zzTd68cUX1blzZ912222udrfeeqvmzp2riRMn6u6771ZmZqZeffVVbdy4UcuXL3cNSZw7d65CQkJ03333KSQkRP/73//0pz/9SXl5eZo+fbrbtbOzs3XZZZfpuuuu04033qjWrVu77V++fLlGjRqljh07avHixWrRooXb/rffflsTJkzQyJEj9dxzz6moqEivv/66hg4dqo0bNyomJsbVtqKiQiNHjtTQoUP1wgsvKCgoyLVv2bJlWrBggW6//XaFhobqb3/7m6655hrt3bvXFewOHz6sQYMGuYJZy5Yt9eWXXyo1NVV5eXm69957zfh1AADOhgEAQC2YM2eOIclYu3ZttW0yMzMNScacOXNc2yZMmGBIMp588km3tn379jX69+/v+nnp0qWGJGPevHlu7b766iuP7UVFRR7XvvXWW42goCCjpKTEtW348OGGJGPmzJke7YcPH240a9bMCA0NNbp3724cOXLEo01+fr4RERFh3HzzzW7bDx06ZISHh7ttr3ydjzzyiMd5JBn+/v5GRkaGa9vmzZsNScYrr7zi2paammpER0cbx44dczv+uuuuM8LDw12vu6r3GQBQuxg6CACol/7whz+4/Txs2DDt2bPH9fOHH36o8PBw/eY3v9GxY8dcj/79+yskJETffvutq21gYKDrz/n5+Tp27JiGDRumoqIi7dq1y+06NptNEydOrLKmwsJC5efnq3Xr1goLC/PY//XXXysnJ0fXX3+9W00+Pj5KSkpyq6nSyT10JxsxYoQ6d+7s+rlXr14KCwtzvQeGYejjjz/WlVdeKcMw3K43cuRI5ebmasOGDVWeGwBQ+xg6CACodwICAtSyZUu3bZGRkW5zr9LT05Wbm6tWrVpVeY4jR464/rx9+3Y99thj+t///qe8vDy3drm5uW4/t23bVv7+/lWeMzY2VuPHj9fDDz+s66+/Xh9++KHbQh7p6emSpIsuuqjK408NZ76+vjrvvPOqbNu+fXuPbSe/B0ePHlVOTo7efPNNvfnmm1We4+T3AABQtwhaAIB6pyarEDocDrVq1Urz5s2rcn9lUMvJydHw4cMVFhamJ598Up07d1ZAQIA2bNighx9+WA6Hw+24k3u/qvLQQw8pOztbzz//vG6++WbNmjVLFovFVZPknKcVFRXlcayvr/tfuzabTVZr1YNLqnsPDMNwu9aNN96oCRMmVNm2V69ep30tAIDaQ9ACADRInTt31jfffKMhQ4acNhwtWbJE2dnZWrBggc4//3zX9szMzLO+9nPPPafjx4/rrbfeUmRkpF588UVXTZLUqlUrjRgx4qzPXxMtW7ZUaGio7HZ7rV8LAHDmmKMFAGiQxo4dK7vdrqlTp3rsq6ioUE5OjqRfeoYqe4IkqaysTK+99to5Xf+NN97QmDFj9NJLL+mpp56SJI0cOVJhYWF65plnVF5e7nHM0aNHz+maJ/Px8dE111yjjz/+WNu2bavVawEAzhw9WgCAWjV79uwq7/90zz33nNN5hw8frltvvVXTpk3Tpk2bdMkll8jPz0/p6en68MMP9fLLL2vMmDEaPHiwIiMjNWHCBN19992yWCx6++233YLX2bBarZo3b55yc3P1+OOPq1mzZrr99tv1+uuv66abblK/fv103XXXqWXLltq7d68+//xzDRkyRK+++uo5Xfdkzz77rL799lslJSXp5ptvVkJCgo4fP64NGzbom2++0fHjx027FgDgzBC0AAC16vXXX69ye0pKyjmfe+bMmerfv7/eeOMNPfroo/L19VVMTIxuvPFGDRkyRJLUvHlzffbZZ7r//vv12GOPKTIyUjfeeKMuvvhijRw58pyu7+/vr08++UQjRozQXXfdpYiICN1www1q06aNnn32WU2fPl2lpaVq27athg0bVu1qhmerdevWWrNmjZ588kktWLBAr732mpo3b67u3bvrueeeM/VaAIAzYzHO9Ss9AAAAAIAb5mgBAAAAgMkIWgAAAABgMoIWAAAAAJiMoAUAAAAAJiNoAQAAAIDJCFoAAAAAYDKCFgDAjcVi0RNPPFGr17jgggt0wQUXmHa+U2t+4oknZLFYdOzYsVo5/9y5c2WxWJSVleXaFhMToyuuuMKU65khKytLFotFc+fO9XYpANAkEbQANFqV/9B84YUXvF2KJOcNei0Wi8LCwlRcXOyxPz09XRaLpV7V7G1LliyRxWLRRx995O1S8LPKEHvqIyAg4FePLSoq0t///nddcsklio6OVmhoqPr27avXX39ddrvdo31GRobGjBmjyMhIBQUFaejQofr222+rPLfD4dDrr7+uPn36KDAwUM2bN9dFF12kzZs3n/NrPlOzZs1St27dFBAQoLi4OL3yyitVtvvmm2904YUXqkWLFoqIiFBiYqLefvvtOq4WQG3x9XYBANCU+Pr6qqioSP/5z380duxYt33z5s1TQECASkpKvFSdU3FxsXx9G9ZfD3Vd80033aTrrrtONputzq5Z37z++usKCQlx/ezj4/Orx+zZs0d33XWXLr74Yt13330KCwvTokWLdPvtt2vVqlX65z//6Wq7b98+JScny8fHRw8++KCCg4M1Z84cXXLJJVq8eLHOP/98t3NPmjRJ8+bN0/jx43XnnXeqsLBQGzdu1JEjR8x70TXwxhtv6A9/+IOuueYa3XfffVq6dKnuvvtuFRUV6eGHH3a1+/TTT3XVVVcpOTnZFV7nz5+v8ePH69ixY/rjH/9Yp3UDMF/D+psUABo4m82mIUOG6L333vMIWu+++64uv/xyffzxx6Zdr6SkRP7+/rJaaz6AoSY9E/VNXdfs4+NTo2DRmI0ZM0YtWrQ4o2OioqK0detWde/e3bXt1ltv1aRJkzRnzhw9/vjjio2NlSQ9++yzysnJ0bZt29S1a1dJ0s0336z4+Hj98Y9/1Pr1613nmD9/vv75z39qwYIF+t3vfmfCq/M0d+5cTZw4UYZhVNumuLhY//d//6fLL7/c1Qt78803y+FwaOrUqbrlllsUGRkpSXr11VcVHR2t//3vf67Afuuttyo+Pl5z584laAGNAEMHATR5R44cUWpqqlq3bq2AgAD17t3b7Zv1StnZ2br
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACL1ElEQVR4nOzdd3yNd//H8ffJTmQZIUYIEsRepVYpalWLVtVqhRgdqkOX3h3u+rVaXbSq2iI6dGlLN1U1SrVaNVs0USlaI0YSEZnn+v1x3TkcSQiu5GS8no/H9TjONT/nOCTv8x2XzTAMQwAAAAAAy7i5ugAAAAAAKGsIWgAAAABgMYIWAAAAAFiMoAUAAAAAFiNoAQAAAIDFCFoAAAAAYDGCFgAAAABYjKAFAAAAABYjaAEAAACAxQhaAIAi0a1bN3Xr1s3VZVhi4cKFstlsSkhIuOhjo6OjFR4ebnlNrnI57wUAlCcELQCApDO/QOcuPj4+atCggSZOnKjDhw+7urwyxWazaeLEiXnWP/3007LZbBozZozsdrsLKgMAWMXD1QUAAEqWJ598UnXr1lV6errWrVun1157TV9//bV27NghPz8/V5fnErfccouGDh0qb2/vIrvGM888o//85z8aNWqU5s2bJzc3vgsFgNKMoAUAcNK3b1+1bdtWkjR27FhVrlxZL774oj777DMNGzbMxdW5hru7u9zd3Yvs/M8995ymTJmiW2+9VQsWLLjskGUYhtLT0+Xr62tRhQCAi8XXZQCA8+revbskae/evZKk7OxsTZs2TfXr15e3t7fCw8P1yCOPKCMjo8BzpKamqkKFCrr77rvzbDtw4IDc3d01ffp0SWe6MK5fv1733XefQkJCVKFCBQ0aNEiJiYl5jp8zZ46aNGkib29v1ahRQ3feeaeSkpKc9unWrZuaNm2qbdu2qWvXrvLz81NERIQ+/vhjSdKaNWvUvn17+fr6qmHDhvruu++cjs9vXNJnn32ma6+9VjVq1JC3t7fq16+vadOmKScn58Jv6llefPFFPfjggxo5cqRiY2OdQpbdbtfMmTPVpEkT+fj4qFq1apowYYJOnDjhdI7w8HD1799fy5cvV9u2beXr66vXX39dq1evls1m00cffaSnnnpKtWrVko+Pj3r06KH4+Pg8tfz888/q06ePgoKC5Ofnp65du2r9+vUX9XoAACaCFgDgvPbs2SNJqly5siSzlevxxx9X69at9dJLL6lr166aPn26hg4dWuA5/P39NWjQIH344Yd5gsj7778vwzA0YsQIp/V33XWXtm7dqieeeEK33367vvjiizzjmqZOnao777xTNWrU0AsvvKAbb7xRr7/+unr16qWsrCynfU+cOKH+/furffv2mjFjhry9vTV06FB9+OGHGjp0qPr166dnnnlGp06d0uDBg3Xy5Mnzvi8LFy6Uv7+/7rvvPs2aNUtt2rTR448/rocffvj8b+hZZs2apcmTJ2v48OFauHBhnpasCRMm6IEHHlCnTp00a9YsjR49WosWLVLv3r3zvL7du3dr2LBhuuaaazRr1iy1bNnSse2ZZ57RkiVLdP/992vKlCn66aef8rzf33//va666iqlpKToiSee0NNPP62kpCR1795dGzduLPRrAgD8jwEAgGEYsbGxhiTju+++MxITE439+/cbH3zwgVG5cmXD19fXOHDggLFlyxZDkjF27FinY++//35DkvH999871nXt2tXo2rWr4/ny5csNScY333zjdGzz5s2d9suto2fPnobdbnesv/feew13d3cjKSnJMAzDOHLkiOHl5WX06tXLyMnJcew3e/ZsQ5KxYMECp1okGe+9955j3a5duwxJhpubm/HTTz/lqTM2NjZPTXv37nWsS0tLy/MeTpgwwfDz8zPS09Md60aNGmXUqVPHaT9JRp06dQxJxrBhw4zs7Ow85/rhhx8MScaiRYuc1i9btizP+txzLVu2zGnfVatWGZKMqKgoIyMjw7F+1qxZhiRj+/bthmEYht1uNyIjI43evXs7vedpaWlG3bp1jWuuuea87wUAIC9atAAATnr27KmQkBCFhYVp6NCh8vf315IlS1SzZk19/fXXkqT77rvP6ZjJkydLkr766qvznrdGjRpatGiRY92OHTu0bds2jRw5Ms/+48ePl81mczzv0qWLcnJy9Pfff0uSvvvuO2VmZuqee+5xagkaN26cAgMD89Ti7+/v1OrWsGFDBQcHKyoqSu3bt3esz/3zX3/9VeBrkeQ0/unkyZM6evSounTporS0NO3ateu8x0pyzORYt27dfMd/LV68WEFBQbrmmmt09OhRx9KmTRv5+/tr1apVTvvXrVtXvXv3zvdao0ePlpeXl+N5ly5dnF7jli1bFBcXp+HDh+vYsWOOa506dUo9evTQ2rVrmQURAC4Sk2EAAJy8+uqratCggTw8PFStWjU1bNjQEWT+/vtvubm5KSIiwumY0NBQBQcHO0JQftzc3DRixAi99tprSktLk5+fnxYtWiQfHx/ddNNNefavXbu20/OKFStKkmN8Uu61GjZs6LSfl5eX6tWrl6eWWrVqOQU3SQoKClJYWFiedWdfpyC///67Hn30UX3//fdKSUlx2pacnHzeYyVp1KhR+vfff/X000+rSpUquvfee522x8XFKTk5WVWrVs33+CNHjjg9r1u3boHXutB7GRcX56ipIMnJyY7jAAAXRtACADhp166dY9bBgpwbWArr1ltv1XPPPaelS5dq2LBheu+999S/f39HuDlbQbP8GYZxSdcu6HyXcp2kpCR17dpVgYGBevLJJ1W/fn35+Pjot99+00MPPVSo1h8PDw999NFH6tOnjyZPnqzg4GCNHj3asd1ut6tq1apOLYBnCwkJcXp+vhkGL/Qac+t97rnnnMZ2nc3f37/A8wMA8iJoAQAKrU6dOrLb7YqLi1NUVJRj/eHDh5WUlKQ6deqc9/imTZuqVatWWrRokWrVqqV9+/bplVdeueRaJHMSiHr16jnWZ2Zmau/everZs+clnbcwVq9erWPHjunTTz/VVVdd5VifOzNjYfn4+Ojzzz/X1VdfrXHjxik4OFiDBg2SJNWvX1/fffedOnXqVOTTtNevX1+SFBgYWKTvGwCUJ4zRAgAUWr9+/SRJM2fOdFr/4osvSpKuvfbaC57jlltu0bfffquZM2eqcuXK6tu37yXV0rNnT3l5eenll192an2aP3++kpOTC1XLpcptITr7upmZmZozZ85FnyswMFDLli1TRESEhg0bppUrV0qShgwZopycHE2bNi3PMdnZ2XmmsL8cbdq0Uf369fX8888rNTU1z/b8ptUHAJwfLVoAgEJr0aKFRo0apTfeeMPRfW7jxo166623NHDgQF199dUXPMfw4cP14IMPasmSJbr99tvl6el5SbWEhIRoypQp+u9//6s+ffro+uuv1+7duzVnzhxdccUV+U6wYZWOHTuqYsWKGjVqlCZNmiSbzaZ33nnnkrs1hoSEaMWKFerUqZMGDhyolStXqmvXrpowYYKmT5+uLVu2qFevXvL09FRcXJwWL16sWbNmafDgwZa8Hjc3N82bN099+/ZVkyZNNHr0aNWsWVP//POPVq1apcDAQH3xxReWXAsAyguCFgDgosybN0/16tXTwoULtWTJEoWGhmrKlCl64oknCnV8tWrV1KtXL3399de65ZZbLquWqVOnKiQkRLNnz9a9996rSpUqafz48Xr66acvOcAVRuXKlfXll19q8uTJevTRR1WxYkWNHDlSPXr0KHDmvwsJCwvTt99+qy5duqhv375au3at5s6dqzZt2uj111/XI488Ig8PD4WHh2vkyJHq1KmTpa+pW7du2rBhg6ZNm6bZs2crNTVVoaGhat++vSZMmGDptQCgPLAZl/r1GwAAl2jQoEHavn274uPjXV0KAABFgjFaAIBidfDgQX311VeX3ZoFAEBJRtdBAECx2Lt3r9avX6958+bJ09OT7mgAgDKNFi0AQLFYs2aNbrnlFu3du1dvvfWWQkNDXV0SAABFhjFaAAAAAGAxWrQAAAAAwGIELQAAAAC
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACO00lEQVR4nOzdeVhUZf/H8c8AggsC7qipaOBui5WmWVr6ZFlmi5lpKoH6VJr5ZItWlmllWZZlWZmCllZuaXumlqZmm2VamaFJWrkv4L4w9++P82OYkR3OMDPwfl0XFzJz5pwvAzLzOfd9f4/DGGMEAAAAALBNkK8LAAAAAIDShqAFAAAAADYjaAEAAACAzQhaAAAAAGAzghYAAAAA2IygBQAAAAA2I2gBAAAAgM0IWgAAAABgM4IWAAAAANiMoAUACAidOnVSp06dbNtfamqqHA6HZsyY4bptzJgxcjgcth0jEOT0PAAAio+gBQAoshkzZsjhcLg+ypcvr8aNG2vo0KHatWuXr8srMZ06dVLLli2z3b5s2TJVrFhRrVu31v79+31QGQDAV0J8XQAAIPCNHTtWDRs21PHjx7Vq1Sq9+uqr+uSTT/TLL7+oYsWKthzj888/t2U/eXnkkUc0cuRIW/b1xRdfqHv37mrSpImWLl2qqlWr2rJfAEBgIGgBAIrt6quv1oUXXihJGjhwoKpVq6bnn39e77//vm699dZi7fvo0aOqWLGiQkND7Sg1TyEhIQoJKf5L44oVK9S9e3c1btzYtpB15MgRVapUqdj7AQCUDKYOAgBsd8UVV0iStm7d6rpt1qxZuuCCC1ShQgVVrVpVvXv31vbt2z0elzkFb+3atbrssstUsWJFPfTQQ677zlyjtXv3biUmJqpWrVoqX768zj33XM2cOTNbPQcPHlR8fLwiIyMVFRWlAQMG6ODBg9m2y22N1qxZs9SmTRtVrFhRVapU0WWXXZbrCNvKlSt1zTXXKDY2VkuXLlW1atU87v/000916aWXqlKlSqpcubKuueYa/frrrx7bxMfHKzw8XFu2bFG3bt1UuXJl9e3bV5LkcDg0dOhQLVq0SC1btlRYWJhatGihzz77LFst//zzjxISElSrVi3XdklJSTnWDQCwFyNaAADbbdmyRZJcIePJJ5/U6NGj1atXLw0cOFB79uzR5MmTddlll+mnn35SVFSU67H79u3T1Vdfrd69e+u2225TrVq1cjzGsWPH1KlTJ23evFlDhw5Vw4YNNW/ePMXHx+vgwYO65557JEnGGPXo0UOrVq3SHXfcoWbNmmnhwoUaMGBAgb6Xxx9/XGPGjFH79u01duxYhYaG6ttvv9UXX3yhK6+80mPb1atXq1u3bmrYsKGWLVum6tWre9z/1ltvacCAAerataueeeYZHT16VK+++qo6dOign376STExMa5tT58+ra5du6pDhw567rnnPKZgrlq1Su+9957uuusuVa5cWS+99JJuuukmbdu2zfWc79q1SxdffLErmNWoUUOffvqpEhMTlZ6eruHDhxfo+wcAFJEBAKCIkpOTjSSzdOlSs2fPHrN9+3bz7rvvmmrVqpkKFSqYv//+26Smpprg4GDz5JNPejx2w4YNJiQkxOP2jh07Gknmtddey3asjh07mo4dO7q+njRpkpFkZs2a5brt5MmTpl27diY8PNykp6cbY4xZtGiRkWQmTJjg2u706dPm0ksvNZJMcnKy6/bHHnvMuL80pqSkmKCgIHPDDTeYjIwMj3qcTqdHbVWrVjWVK1c2LVq0MLt3785W/6FDh0xUVJQZNGiQx+07d+40kZGRHrcPGDDASDIjR47Mth9JJjQ01GzevNl1288//2wkmcmTJ7tuS0xMNLVr1zZ79+71eHzv3r1NZGSkOXr0qDHGmK1bt2Z7HgAAxcfUQQBAsXXp0kU1atRQvXr11Lt3b4WHh2vhwoWqW7eu3nvvPTmdTvXq1Ut79+51fURHRysuLk5ffvmlx77CwsJ0++2353vMTz75RNHR0R5rwMqVK6dhw4bp8OHDWrFihWu7kJAQ3Xnnna7tgoODdffdd+d7jEWLFsnpdOrRRx9VUJDnS+aZUwyPHDmiQ4cOqVatWoqIiMi2ryVLlujgwYO69dZbPZ6H4OBgtW3bNtvzIMmjZnddunTR2Wef7fr6nHPOUUREhP78809J1ijeggUL1L17dxljPI7XtWtXpaWl6ccff8z3+wcAFB1TBwEAxfbKK6+ocePGCgkJUa1atdSkSRNXMElJSZExRnFxcTk+tly5ch5f161bt0CNL/766y/FxcVlC0DNmjVz3Z/5uXbt2goPD/fYrkmTJvkeY8uWLQoKClLz5s3z3TY2Nlb9+/fXgw8+qFtvvVXz5s1TcHCw6/6UlBRJWevXznRmOAsJCdFZZ52V47b169fPdluVKlV04MABSdKePXt08OBBTZ06VVOnTs1xH7t37873ewIAFB1BCwBQbG3atHF1HTyT0+mUw+HQp59+6hE8Mp0ZgCpUqOCVGkvCAw88oH379mnChAkaNGiQpk+f7hr5cjqdkqx1WtHR0dkee2a3w7CwsGwhMlNOz6NkjWS5H+u2227LdS3aOeecU4DvCABQVAQtAIBXnX322TLGqGHDhmrcuLFt+23QoIHWr18vp9PpEUh+//131/2Zn5ctW6bDhw97hLpNmzYVqHan06nffvtN5513XoHqeuaZZ7R//35NmzZNVapU0cSJE137kqSaNWuqS5cuBdpXUdWoUUOVK1dWRkaG148FAMgZa7QAAF514403Kjg4WI8//rhrxCWTMUb79u0r0n67deumnTt3as6cOa7bTp8+rcmTJys8PFwdO3Z0bXf69Gm9+uqrru0yMjI0efLkfI9x/fXXKygoSGPHjnWNErnXnpvXX39dPXv21PPPP68nnnhCktS1a1dFREToqaee0qlTp7I9Zs+ePfnWU1DBwcG66aabtGDBAv3yyy9ePRYAIGeMaAEAvOrss8/WE088oVGjRik1NVXXX3+9KleurK1bt2rhwoUaPHiw7rvvvkLvd/DgwXr99dcVHx+vtWvXKiYmRvPnz9fq1as1adIkVa5cWZLUvXt3XXLJJRo5cqRSU1PVvHlzvffee0pLS8v3GLGxsXr44Yc1btw4XXrppbrxxhsVFham77//XnXq1NH48eNzfFxQUJBmz56ttLQ0jR49WlWrVtVdd92lV199Vf369VPr1q3Vu3dv1ahRQ9u2bdPHH3+sSy65RC+//HKhn4fcPP300/ryyy/Vtm1bDRo0SM2bN9f+/fv1448/aunSpdq/f79txwIAZEfQAgB43ciRI9W4cWO98MILevzxxyVJ9erV05VXXqnrrruuSPusUKGCli9frpEjR2rmzJlKT09XkyZNlJycrPj4eNd2QUFB+uCDDzR8+HDNmjVLDodD1113nSZOnKjzzz8/3+OMHTtWDRs21OTJk/Xwww+rYsWKOuecc9SvX788HxcaGqqFCxeqS5cuuvvuuxUVFaU+ffqoTp06evrpp/Xss8/qxIkTqlu3ri699NICdVosjFq1aum7777T2LFj9d5772nKlCmqVq2aWrRooWeeecbWYwEAsnOYvOY+AAAAAAAKjTVaAAAAAGAzghYAAAAA2IygBQAAAAA2I2gBAAAAgM0IWgAAAABgM4IWAAAAANiMoAUA8OBwODRmzBivHqNTp07q1KmTbfs7s+YxY8bI4XBo7969Xtn/jBkz5HA4lJqa6rotJiZG1157rS3Hs0NqaqocDodmzJjh61IAoEwiaAEotTLfaD733HO+LkWSFB8fL4fDoYiICB07dizb/SkpKXI4HH5Vs68tX75cDodD8+fP93Up+H+ZIfbMj/Llyxd6XwcPHlTNmjUL9DN+8skn5XA41LJlyxzv//rrr9WhQwdVrFhR0dHRGjZsmA4fPlzomorL6XRqwoQJatiwocqXL69zzjlH77zzTo7bzp07VxdffLGioqJUrVo1dezYUR9//HEJVwzAW0J8XQAAlCUhISE6evSoPvzwQ/Xq1cvjvtmzZ6t8+fI6fvy4j6qzHDt2TCEhgfX
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAANsCAYAAABCgdkuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADW00lEQVR4nOzdd3gU5f7+8fdk0zuhJbQECB2kd6Q3FQULoKLS8ViO+rMcxfO1F+yKxy5NERUUQcSCiPTepRtKaBJ6EkJ6dn5/jLsQEyBAwmST+3VduTbMzO7eGzbZ/ew8n+cxTNM0ERERERERkULjZXcAERERERGRkkaFloiIiIiISCFToSUiIiIiIlLIVGiJiIiIiIgUMhVaIiIiIiIihUyFloiIiIiISCFToSUiIiIiIlLIVGiJiIiIiIgUMhVaIiIiIiIihUyFloiIiIdasGABhmGwYMECu6OIiMg/qNASERHbTJo0CcMw3F/e3t5UrlyZIUOGcPDgwTzHd+7cOdfxvr6+VK9enVGjRrF///7z3vbZX0888YT7uJiYGPr06ZPnviZPnozD4aB3796kp6cX/oMXEZESzdvuACIiIs8//zzVq1cnPT2dFStWMGnSJJYsWcLmzZvx9/fPdWyVKlUYM2YMAJmZmWzdupWPPvqIOXPmsG3bNgIDA/O97bM1bNjwvHmmTJnCkCFD6N69OzNnzsyTQURE5EJUaImIiO2uueYaWrRoAcCIESMoV64cr776KrNmzWLAgAG5jg0LC+OOO+7Ita169ercf//9LF26lB49epzztgvi66+/ZvDgwXTt2pXvv/++UIqs1NTUPAWgiIiUbBo6KCIixc7VV18NwK5duwp0fGRkJADe3pf3+eG0adO444476Ny5M7NmzcpTZH3xxRc0b96cgIAAIiIiuPXWW/MMWezcuTMNGzZk7dq1dOzYkcDAQJ588kni4+MxDIM33niDTz75hJo1a+Ln50fLli1ZvXp1nizbt2/nlltuISIiAn9/f1q0aMGsWbMu6/GJiMiVozNaIiJS7MTHxwNQpkyZPPtycnI4duwYAFlZWWzbto1nnnmG2NhY2rdvn+f4pKQk9/Eu5cqVy3Pc9OnTGTRoEB07duSHH34gICAg1/6XXnqJp556igEDBjBixAiOHj3K//73Pzp27Mj69esJDw93H3v8+HGuueYabr31Vu644w4qVqzo3vfll19y6tQp7r77bgzD4LXXXuOmm25i9+7d+Pj4ALBlyxbat29P5cqVeeKJJwgKCmLatGn069eP6dOnc+ONNxbsBykiIrZRoSUiIrZzFUPp6emsXLmS5557Dj8/v3wnqdi+fTvly5fPta1evXr8+uuv+Pr65jm+e/fuebaZppnr3+vXr2fOnDl06NCB2bNn5ymy9u7dyzPPPMOLL77Ik08+6d5+00030bRpUz744INc2xMSEvjoo4+4++673dtcxeO+ffuIi4tzF5F16tShb9++zJkzx/14H3zwQapVq8bq1avx8/MD4N5776VDhw48/vjjKrRERDyACi0REbHdP4uhmJgYvvjiC6pUqZLn2JiYGD799FMAsrOz2bFjB6+99hrXXHMNixcvzlOEvf/++9SuXfu893/ixAmys7OpUqVKniIL4LvvvsPpdDJgwIBcZ8ciIyOpVasW8+fPz1Vo+fn5MXTo0Hzva+DAgbnO1LmGSe7evdud5ffff+f555/n1KlTnDp1yn1sr169eOaZZzh48CCVK1c+72MSERF7qdASERHbuYqhpKQkJkyYwKJFi9xncv4pKCgoV2HWu3dvOnToQIsWLXjllVd48803cx3fqlWrC06G0a1bN6pVq8aHH35IREQEY8eOzbU/Li4O0zSpVatWvtd3DflzqVy5cr5n1wCqVauW69+uouvkyZMA7Ny5E9M0eeqpp3jqqafyvY0jR46o0BIRKeZUaImIiO3OLob69etHhw4duP3229mxYwfBwcEXvH7z5s0JCwtj0aJFl5zhvffe4+TJk7z77ruUKVOGZ5991r3P6XRiGAY///wzDocjz3X/mTG/s2Iu+V0fzgxndDqdADz66KP06tUr32NjY2PP+1hERMR+KrRERKRYcTgcjBkzhi5duvDee+/lWlz4fHJyckhJSbnk+/Xy8uLzzz8nKSmJ5557joiICB544AEAatasiWmaVK9e/YLDEC9XjRo1AOssWX79ZSIi4hk0vbuIiBQ7nTt3plWrVrzzzjukp6df8Pj58+eTkpJC48aNL+t+fXx8+Pbbb2nfvj0PPfQQkydPBqxJLxwOB88991yeiTRM0+T48eOXdb9nq1ChAp07d+bjjz/m0KFDefYfPXq00O5LRESKjs5oiYhIsfTYY4/Rv39/Jk2axL/+9S/39qSkJL744gvgzGQYH374IQEBAQU++3U+gYGB/Pjjj3Tq1Ilhw4YRFhbGDTfcwIsvvsjo0aOJj4+nX79+hISEsGfPHmbMmMGoUaN49NFHL/u+Xd5//306dOhAo0aNGDlyJDVq1ODw4cMsX76cAwcOsHHjxkK7LxERKRoqtEREpFi66aabqFmzJm+88QYjR4509zYdOHCAO++8EwDDMChTpgydOnXimWeeoUmTJoVy32FhYe7p3gcOHMjPP//ME088Qe3atXn77bd57rnnAKhatSo9e/bkhhtuKJT7dalfvz5r1qzhueeeY9KkSRw/fpwKFSrQtGlTnn766UK9LxERKRqG+c8xECIiIiIiInJZ1KMlIiIiIiJSyFRoiYiIiIiIFDIVWiIiIiIiIoVMhZaIiIiIiEghU6ElIiIiIiJSyFRoiYiIiIiIFDIVWiIikothGDz77LNFeh+dO3emc+fOhXZ7/8z87LPPYhgGx44dK5LbnzRpEoZhEB8f794WExNDnz59CuX+CkN8fDyGYTBp0iS7o4iIlEoqtESkxHK90XzjjTfsjgLAkCFDMAyD0NBQ0tLS8uyPi4vDMIxildluCxYswDAMvv32W7ujyFm+/vprmjVrhr+/P+XLl2f48OEXVdQuW7aMDh06EBgYSGRkJA888AApKSl5jouLi+PWW2+lSpUqBAYGUrduXZ5//nlSU1PPeduJiYlUqFDBtueN0+nktddeo3r16vj7+3PVVVfx1Vdf5Xvstm3b6N27N8HBwURERHDnnXdy9OjRK5xYRIqKt90BRERKE29vb1JTU/nhhx8YMGBArn1TpkzB39+f9PR0m9JZ0tLS8Pb2rJeHK535zjvv5NZbb8XPz++K3Wdx8eGHH3LvvffSrVs33nrrLQ4cOMDYsWNZs2YNK1euxN/f/7zX37BhA926daNevXru67/xxhvExcXx888/u4/bv38/rVq1IiwsjPvvv5+IiAiWL1/OM888w9q1a/n+++/zvf2nn376vIVYUfvvf//LK6+8wsiRI2nZsiXff/89t99+O4ZhcOutt7qPO3DgAB07diQsLIyXX36ZlJQU3njjDTZt2sSqVavw9fW17TGISOHwrFdSEREP5+fnR/v27fnqq6/yFFpffvkl1113HdOnTy+0+0tPT8fX1xcvr4IPYLjQG+Xi6EpndjgcOByOK3qfxUFmZiZPPvkkHTt2ZO7cuRiGAUC7du24/vrr+fTTT/n3v/993tt48sknKVOmDAsWLCA0NBSwhl2OHDmSX3/9lZ49ewIwefJkEhMTWbJkCQ0aNABg1KhROJ1OPv/8c06ePEmZMmVy3fbmzZv58MMPefrpp3n66acL9bE/++yzTJo0Kddw0X86ePAgb775Jvfddx/vvfceACNGjKBTp0489thj9O/f3/28efnllzl9+jRr166lWrVqALRq1YoePXowadIkRo0aVaj5ReTK09BBESn1jhw5wvDhw6lYsSL+/v40btyYzz77LM9xx48f58477yQ0NJTw8HAGDx7Mxo0bL7oP5vbbb+fnn38mMTHRvW316tXExcVx++235zn+xIkTPProozRq1Ijg4GBCQ0O55ppr2LhxY67jXMPsvv76a/7v//6PypUrExgYSHJyMgDffPMN9evXx9/fn4YNGzJjxgyGDBlCTExMrts5V7/Tzp0
     },
     "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": 17,
   "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": 18,
   "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": 19,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_12816/3712199621.py:45: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
      "/tmp/ipykernel_12816/3712199621.py:46: 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/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACOz0lEQVR4nOzdfXzN9f/H8efZ9WZXrmdik8lFVEgalcplVFRCLSyicl0U+iEhchkSKpouFpFIutAsFLlKSUpC05TLXGyY2XbO5/fH+e5wbNh8tp1dPO6327kd+3w+5/N5nYvtnKf3+/M6FsMwDAEAAAAArpmbqwsAAAAAgKKOYAUAAAAAJhGsAAAAAMAkghUAAAAAmESwAgAAAACTCFYAAAAAYBLBCgAAAABMIlgBAAAAgEkEKwAAAAAwiWAFACgQ+/fvl8Vi0YIFC1xdSrHG4wwArkGwAgCYtmDBAlksFv3444+uLiVP3H333apbt26W5fHx8fLz81ODBg104sQJF1QGACisPFxdAACgZAgLC9O5c+fk6enp6lKuybfffqsHHnhANWvW1OrVq1WmTBlXlwQAKEQYsQIAFAiLxSIfHx+5u7u7upTLSklJyXb5unXr9MADD+iGG27Is1B19uxZ0/sAABQeBCsAQIHI7tyf6Oho+fv7699//1WHDh3k7++v8uXLa8iQIbJarU63t9lsmj59um688Ub5+PioYsWKevrpp3Xy5Emn7T777DO1a9dOoaGh8vb2VvXq1TV27Ngs+8uc7rdt2zbddddd8vPz00svvZSl7u+//17t2rVTRESEVq9erbJlyzqt/+qrr3TnnXeqVKlSCggIULt27fTbb785bZN5P/ft26e2bdsqICBAUVFRkuyBs1+/flq+fLnq1q0rb29v3Xjjjfr666+z1PLvv/+qR48eqlixomO7d9999+oPPgAg3zEVEADgUlarVa1bt1bjxo01ZcoUrV69WlOnTlX16tX17LPPOrZ7+umntWDBAj355JMaMGCAEhISNGvWLP3888/asGGDY4rhggUL5O/vr+eff17+/v769ttvNWrUKCUnJ2vy5MlOxz5+/Ljuu+8+denSRU888YQqVqzotH7Dhg1q27atqlWrpvj4eJUrV85p/QcffKDu3burdevWmjhxolJSUjRnzhzdcccd+vnnnxUeHu7YNiMjQ61bt9Ydd9yhKVOmyM/Pz7Fu/fr1+vTTT9WnTx8FBARo5syZeuSRR5SYmOgIckeOHNHtt9/uCGLly5fXV199pZ49eyo5OVmDBg3Ki6cDAHCtDAAATIqJiTEkGVu3br3sNgkJCYYkIyYmxrGse/fuhiRjzJgxTtvWr1/faNiwoePn77//3pBkxMbGOm339ddfZ1mekpKS5dhPP/204efnZ6SmpjqWNWvWzJBkzJ07N8v2zZo1M8qUKWMEBAQYN954o3H06NEs25w+fdoIDg42evXq5bT88OHDRlBQkNPyzPs5bNiwLPuRZHh5eRl79+51LPvll18MScYbb7zhWNazZ0+jUqVKxn///ed0+y5duhhBQUGO+53d4wwAyH9MBQQAuNwzzzzj9POdd96pv/76y/HzkiVLFBQUpJYtW+q///5zXBo2bCh/f3+tWbPGsa2vr6/j36dPn9Z///2nO++8UykpKfrjjz+cjuPt7a0nn3wy25rOnj2r06dPq2LFigoMDMyyPi4uTqdOndJjjz3mVJO7u7saN27sVFOmi0fgLtaiRQtVr17d8fNNN92kwMBAx2NgGIaWLl2qBx54QIZhOB2vdevWSkpK0k8//ZTtvgEABYOpgAAAl/Lx8VH58uWdlpUuXdrp3Kk9e/YoKSlJFSpUyHYfR48edfz7t99+04gRI/Ttt98qOTnZabukpCSnnytXriwvL69s9xkREaFu3bpp6NCheuyxx7RkyRKnxht79uyRJN17773Z3v7SMObh4aHrrrsu222rVq2aZdnFj8GxY8d06tQpvf3223r77bez3cfFjwEAoOARrAAALpWTLoE2m00VKlRQbGxstuszg9mpU6fUrFkzBQYGasyYMapevbp8fHz0008/aejQobLZbE63u3h0Kzsvvviijh8/rkmTJqlXr16aP3++LBaLoybJfp5VSEhIltt6eDi/xXp7e8vNLfuJIpd7DAzDcDrWE088oe7du2e77U033XTF+wIAyF8EKwBAoVe9enWtXr1aTZs2vWIYWrt2rY4fP65PP/1Ud911l2N5QkLCNR974sSJOnHihObNm6fSpUtr6tSpjpokqUKFCmrRosU17z8nypcvr4CAAFmt1nw/FgDg2nCOFQCg0OvUqZOsVqvGjh2bZV1GRoZOnTol6cLIT+ZIjySlpaVp9uzZpo7/1ltvqWPHjpo2bZrGjRsnSWrdurUCAwM1fvx4paenZ7nNsWPHTB3zYu7u7nrkkUe0dOlS7dy5M1+PBQC4NoxYAQDyzLvvvpvt9y8NHDjQ1H6bNWump59+WhMmTND27dvVqlUreXp6as+ePVqyZIlmzJihjh07qkmTJipdurS6d++uAQMGyGKx6IMPPnAKWtfCzc1NsbGxSkpK0siRI1WmTBn16dNHc+bMUdeuXdWgQQN16dJF5cuXV2Jior744gs1bdpUs2bNMnXci7322mtas2aNGjdurF69eqlOnTo6ceKEfvrpJ61evVonTpzIs2MBAHKPYAUAyDNz5szJdnl0dLTpfc+dO1cNGzbUW2+9pZdeekkeHh4KDw/XE088oaZNm0qSypYtq5UrV2rw4MEaMWKESpcurSeeeELNmzdX69atTR3fy8tLy5YtU4sWLdS/f38FBwfr8ccfV2hoqF577TVNnjxZ58+fV+XKlXXnnXdettvgtapYsaK2bNmiMWPG6NNPP9Xs2bNVtmxZ3XjjjZo4cWKeHgsAkHsWw+x/4wEAAABACcc5VgAAAABgEsEKAAAAAEwiWAEAAACASQQrAAAAADCJYAUAAAAAJhGsAAAAAMAkghUAlGAWi0WjR4/O12Pcfffduvvuu/Nsf5fWPHr0aFksFv3333/5sv8FCxbIYrFo//79jmXh4eG6//778+R4eWH//v2yWCxasGCBq0sBgBKLYAWgWMj8YDllyhRXlyLJ/oW4FotFgYGBOnfuXJb1e/bskcViKVQ1u9ratWtlsVj0ySefuLoUXEbLli1lsVjUr1+/HG3/zTffqGfPnqpbt67c3d0VHh6e7XaZ4fhylw0bNjhtP2vWLNWuXVve3t6qXLmynn/+eZ09e9bs3cu1OXPm6NFHH1XVqlVlsVgu+0XYhw4d0rBhw3TPPfcoICBAFotFa9euLdBaAeQ/D1cXAADFlYeHh1JSUvT555+rU6dOTutiY2Pl4+Oj1NRUF1Vnd+7cOXl4FK23goKuuWvXrurSpYu8vb0L7JiF0aeffqqNGzfm6jYfffSRPv74YzVo0EChoaGX3e7hhx9WREREluUvvfSSzpw5o0aNGjmWDR06VJMmTVLHjh01cOBA/f7773rjjTf022+/adWqVbmqz6yJEyfq9OnTuu2223To0KHLbrd7925NnDhRNWrUUL169XL9OAIoGorWuykAFCHe3t5q2rSpFi5cmCVYffTRR2rXrp2WLl2aZ8dLTU2Vl5eX3NxyPhnBx8cnz45fUAq6Znd3d7m7uxfoMQub1NRUDR48WEOHDtWoUaNyfLvx48frnXfekaenp+6//37t3Lkz2+1uuukm3XTTTU7LDhw4oH/++UdPPfWUvLy8JNlHfqZNm6auXbvq/fffd2x7ww03qH///vr888/1wAMPXMM9dLZgwQI9+eSTMgzjitutW7fOMVrl7+9/2e0aNmyo48ePq0yZMvrkk0/06KOPmq4RQOHDVEAAJcrRo0fVs2dPVaxYUT4+Prr55pv13nvvZdnu+PHj6tq1qwIDAxUcHKzu3bvrl19+yfV5LI8//ri++uornTp1yrFs69at2rNnjx5//PEs2584cUJDhgxRvXr15O/vr8DAQN1333365ZdfnLbLnDa3aNEijRgxQpUrV5afn5+Sk5MlSUuWLFGdOnXk4+OjunXratm
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACkLElEQVR4nOzdeZyN5f/H8feZGbMZM2MfI8xYsleWCNWQnYpKshtr2StafKNECCVKSBiqsYuEEgrxRSilkkaNpaxhRmPMeu7fH/fP+TpnBoOZuWd5PR+P8zjXue/73PfnLOq857ru67YZhmEIAAAAAODgZnUBAAAAAJDTEJQAAAAAwAVBCQAAAABcEJQAAAAAwAVBCQAAAABcEJQAAAAAwAVBCQAAAABcEJQAAAAAwAVBCQAAAABcEJQAAOlq3LixGjdubHUZmWLBggWy2Ww6cuTITT83PDxcISEhmV6TVW7nvQCA/ISgBAB5xJUfwFdu3t7euvPOOzV48GCdPn3a6vLyFJvNpsGDB6dZPmHCBNlsNvXu3Vt2u92CygAAmcXD6gIAAJlr7NixCg0NVUJCgrZv365Zs2Zp/fr1+vnnn+Xr62t1eZbo3r27OnXqJC8vryw7xptvvqlXXnlFPXv21Ny5c+Xmxt8iASA3IygBQB7TunVr1a1bV5LUt29fFS1aVFOnTtVnn32mzp07W1ydNdzd3eXu7p5l+58yZYpGjhypHj16aP78+bcdkgzDUEJCgnx8fDKpQgDAzeLPXQCQxz300EOSpOjoaElSSkqKxo0bpwoVKsjLy0shISH6z3/+o8TExGvuIy4uTgULFtSwYcPSrPvrr7/k7u6uiRMnSvrfEMAdO3bo+eefV/HixVWwYEE99thjOnv2bJrnz5w5U9WrV5eXl5eCg4M1aNAgxcTEOG3TuHFj1ahRQz/99JPCwsLk6+urihUrasWKFZKkrVu3qn79+vLx8VHlypW1adMmp+end17OZ599prZt2yo4OFheXl6qUKGCxo0bp9TU1Bu/qVeZOnWqXnzxRXXr1k0RERFOIclut2vatGmqXr26vL29VbJkST399NO6cOGC0z5CQkL08MMPa8OGDapbt658fHz0wQcfaMuWLbLZbFq2bJnGjx+vO+64Q97e3mratKkOHz6cppbdu3erVatWCggIkK+vr8LCwrRjx46bej0AABNBCQDyuD/++EOSVLRoUUlmL9Orr76q2rVr65133lFYWJgmTpyoTp06XXMffn5+euyxx7R06dI0QWLx4sUyDENdu3Z1Wj5kyBD9+OOPeu211zRgwAB9/vnnac7rGTNmjAYNGqTg4GC9/fbbeuKJJ/TBBx+oRYsWSk5Odtr2woULevjhh1W/fn1NnjxZXl5e6tSpk5YuXapOnTqpTZs2evPNN3Xp0iV16NBB//7773XflwULFsjPz0/PP/+8pk+frjp16ujVV1/Vyy+/fP039CrTp0/X8OHD1aVLFy1YsCBNT9LTTz+tF154QY0aNdL06dPVq1cvRUZGqmXLlmle36FDh9S5c2c1b95c06dP1z333ONY9+abb2rVqlUaMWKERo4cqV27dqV5v7/++ms9+OCDunjxol577TVNmDBBMTExeuihh/Tdd99l+DUBAP6fAQDIEyIiIgxJxqZNm4yzZ88ax48fN5YsWWIULVrU8PHxMf766y9j//79hiSjb9++Ts8dMWKEIcn4+uuvHcvCwsKMsLAwx+MNGzYYkowvvvjC6bl33XWX03ZX6mjWrJlht9sdy5977jnD3d3diImJMQzDMM6cOWN4enoaLVq0MFJTUx3bzZgxw5BkzJ8/36kWScaiRYscy3777TdDkuHm5mbs2rUrTZ0RERFpaoqOjnYsi4+PT/MePv3004avr6+RkJDgWNazZ0+jXLlyTttJMsqVK2dIMjp37mykpKSk2de3335rSDIiIyOdln/55Zdpll/Z15dffum07TfffGNIMqpWrWokJiY6lk+fPt2QZBw4cMAwDMOw2+1GpUqVjJYtWzq95/Hx8UZoaKjRvHnz674XAIC06FECgDymWbNmKl68uMqUKaNOnTrJz89Pq1atUunSpbV+/XpJ0vPPP+/0nOHDh0uS1q1bd939BgcHKzIy0rHs559/1k8//aRu3bql2b5///6y2WyOxw888IBSU1N19OhRSdKmTZuUlJSkZ5991qknpl+/fvL3909Ti5+fn1OvV+XKlRUYGKiqVauqfv36juVX2n/++ec1X4skp/N//v33X/3zzz964IEHFB8fr99+++26z5XkmEkwNDQ03fOfli9froCAADVv3lz//POP41anTh35+fnpm2++cdo+NDRULVu2TPdYvXr1kqenp+PxAw884PQa9+/fr6ioKHXp0kXnzp1zHOvSpUtq2rSptm3bxix8AHCTmMwBAPKY999/X3feeac8PDxUsmRJVa5c2RFEjh49Kjc3N1WsWNHpOUFBQQoMDHSEmPS4ubmpa9eumjVrluLj4+Xr66vIyEh5e3vrySefTLN92bJlnR4XLlxYkhzn51w5VuXKlZ228/T0VPny5dPUcscddzgFL0kKCAhQmTJl0iy7+jjX8ssvv2jUqFH6+uuvdfHiRad1sbGx132uJPXs2VMnTpzQhAkTVKxYMT333HNO66OiohQbG6sSJUqk+/wzZ844PQ4NDb3msW70XkZFRTlqupbY2FjH8wAAN0ZQAoA8pl69eo5Z767FNXBkVI8ePTRlyhStXr1anTt31qJFi/Twww87wsnVrjXLnGEYt3Tsa+3vVo4TExOjsLAw+fv7a+zYsapQoYK8vb31/fff66WXXspQ74uHh4eWLVumVq1aafjw4QoMDFSvXr0c6+12u0qUKOHUA3e14sWLOz2+3gx3N3qNV+qdMmWK07lNV/Pz87vm/gEAaRGUACAfKVeunOx2u6KiolS1alXH8tOnTysmJkblypW77vNr1KihWrVqKTIyUnfccYeOHTum995775ZrkcxJDMqXL+9YnpSUpOjoaDVr1uyW9psRW7Zs0blz5/Tpp5/qwQcfdCy/MjNgRnl7e2vNmjVq0qSJ+vXrp8DAQD322GOSpAoVKmjTpk1q1KhRlk/zXaFCBUmSv79/lr5vAJCfcI4SAOQjbdq0kSRNmzbNafnUqVMlSW3btr3hPrp3766vvvpK06ZNU9GiRdW6detbqqVZs2by9PTUu+++69T7M2/ePMXGxmaollt1pYfm6uMmJSVp5syZN70vf39/ffnll6pYsaI6d+6szZs3S5I6duyo1NRUjRs3Ls1zUlJS0kyBfjvq1KmjChUq6K233lJcXFya9elNyw4AuD56lAAgH7n77rvVs2dPzZkzxzH87LvvvtPChQvVvn17NWnS5Ib76NKli1588UWtWrVKAwYMUIECBW6pluLFi2vkyJF6/fXX1apVKz366KM6dOiQZs6cqXvvvTfdCSIyS8OGDVW4cGH17NlTQ4cOlc1m08cff3zLwwKLFy+ujRs3qlGjRmrfvr02b96ssLAwPf3005o4caL279+vFi1aqECBAoqKitLy5cs1ffp0dejQIVNej5ubm+bOnavWrVurevXq6tWrl0qXLq2///5b33zzjfz9/fX5559nyrEAIL8gKAFAPjN37lyVL19eCxYs0KpVqxQUFKSRI0fqtddey9DzS5YsqRYtWmj9+vXq3r37bdUyZswYFS9eXDNmzNBzzz2nIkWKqH///powYcItB7CMKFq0qNauXavhw4dr1KhRKly4sLp166amTZtec+a5GylTpoy++uorPfDAA2rdurW2bdum2bNnq06dOvrggw/0n//8Rx4eHgoJCVG3bt3UqFGjTH1NjRs31s6dOzVu3DjNmDFDcXFxCgoKUv369fX0009n6rEAID+wGbf65zMAQL712GOP6cCBAzp8+LDVpQAAkCU4RwkAcFNOnjypdevW3XZvEgAAORlD7wAAGRIdHa0dO3Zo7ty5KlCgAMO5AAB5Gj1KAIAM2bp1q7p3767o6GgtXLhQQUFBVpcEAECW4RwlAAAAAHBBjxIAAAAAuCAoAQAAAIALghI
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC280lEQVR4nOzde3zO9f/H8ee1zU52MqchbMxhTuWUkEM5k0hChEX0y7moqJQckoooXzqo6aBEDkkklilyLlGk0UTlFDYxY9v1+f3xaVeua8PGts8Oj/vtdt3en+vz+Vyfz+u6NnU9935/3h+bYRiGAAAAAAAOblYXAAAAAAB5DUEJAAAAAFwQlAAAAADABUEJAAAAAFwQlAAAAADABUEJAAAAAFwQlAAAAADABUEJAAAAAFwQlAAAAADABUEJAJArWrZsqZYtW2bb8Q4dOiSbzab58+c71k2YMEE2my3bzpEfZPQ5AABuHEEJAAqx+fPny2azOR7e3t6qWrWqhg0bpuPHj1tdXq5p2bKlatWqlW59dHS0fH19Va9ePZ0+fdqCygAAVvGwugAAgPUmTpyosLAwJSUlaePGjZo7d65WrVqln376Sb6+vtlyjq+++ipbjnM1zzzzjMaOHZstx/r666/VuXNnVatWTevWrVNwcHC2HBcAkD8QlAAA6tChgxo0aCBJeuihh1S8eHHNmDFDn332me6///4bOnZiYqJ8fX3l6emZHaVelYeHhzw8bvx/bRs2bFDnzp1VtWrVbAtJ58+fV9GiRW/4OACA3MHQOwBAOnfeeackKS4uzrHuww8/VP369eXj46Pg4GD16tVLR44ccXpd2hC2nTt3qnnz5vL19dVTTz3l2OZ6jdKJEyc0cOBAlS5dWt7e3rr55pv13nvvpasnPj5ekZGRCgwMVFBQkPr376/4+Ph0+13pGqUPP/xQt956q3x9fVWsWDE1b978ij1c3377rTp16qTw8HCtW7dOxYsXd9q+evVqNWvWTEWLFpW/v786deqkn3/+2WmfyMhI+fn56eDBg+rYsaP8/f3Vp08fSZLNZtOwYcO0fPly1apVS15eXqpZs6a+/PLLdLX8+eefGjBggEqXLu3Y7913382wbgBA9qJHCQCQzsGDByXJERKmTJmi8ePHq0ePHnrooYd08uRJvf7662revLl++OEHBQUFOV576tQpdejQQb169dIDDzyg0qVLZ3iOCxcuqGXLljpw4ICGDRumsLAwLV68WJGRkYqPj9fIkSMlSYZhqEuXLtq4caP+7//+TxEREVq2bJn69++fqffy/PPPa8KECWrSpIkmTpwoT09Pbd26VV9//bXatm3rtO+mTZvUsWNHhYWFKTo6WiVKlHDa/sEHH6h///5q166dpk2bpsTERM2dO1e33367fvjhB4WGhjr2TUlJUbt27XT77bfrlVdecRrCuHHjRi1dulRDhgyRv7+/XnvtNd177706fPiw4zM/fvy4brvtNkewKlmypFavXq2BAwfq7NmzGjVqVKbePwDgOhkAgEIrKirKkGSsW7fOOHnypHHkyBFj4cKFRvHixQ0fHx/jjz/+MA4dOmS4u7sbU6ZMcXrtnj17DA8PD6f1LVq0MCQZb7zxRrpztWjRwmjRooXj+cyZMw1JxocffuhYd+nSJaNx48aGn5+fcfbsWcMwDGP58uWGJOOll15y7JeSkmI0a9bMkGRERUU51j/33HPG5f9ri42NNdzc3Ix77rnHSE1NdarHbrc71RYcHGz4+/sbNWvWNE6cOJGu/n/++ccICgoyBg0a5LT+2LFjRmBgoNP6/v37G5KMsWPHpjuOJMPT09M4cOCAY92PP/5oSDJef/11x7qBAwcaZcqUMf7++2+n1/fq1csIDAw0EhMTDcMwjLi4uHSfAwDgxjH0DgCg1q1bq2TJkipfvrx69eolPz8/LVu2TOXKldPSpUtlt9vVo0cP/f33345HSEiIqlSpovXr1zsdy8vLSw8++OA1z7lq1SqFhIQ4XQNVpEgRjRgxQufOndOGDRsc+3l4eOiRRx5x7Ofu7q7hw4df8xzLly+X3W7Xs88+Kzc35//luQ7RO3/+vP755x+VLl1aAQEB6Y61du1axcfH6/7773f6HNzd3dWoUaN0n4Mkp5ov17p1a1WuXNnxvE6dOgoICNBvv/0myexFW7JkiTp37izDMJzO165dOyUkJOj777+/5vsHAFw/ht4BAPS///1PVatWlYeHh0qXLq1q1ao5gkVsbKwMw1CVKlUyfG2RIkWcnpcrVy5TEzf8/vvvqlKlSroAExER4die1pYpU0Z+fn5O+1WrVu2a5zh48KDc3NxUo0aNa+4bHh6ufv366cknn9T999+vxYsXy93d3bE9NjZW0n/Xb7lyDVceHh666aabMty3QoUK6dYVK1ZMZ86ckSSdPHlS8fHxeuutt/TWW29leIwTJ05c8z0BAK4fQQkAoFtvvdUx650ru90um82m1atXOwWHNK4BxsfHJ0dqzA1PPPGETp06pZdeekmDBg3SO++84+h5stvtkszrlEJCQtK91nW2PS8vr3QhME1Gn6Nk9iRdfq4HHnjgitdi1alTJxPvCABwvQhKAICrqly5sgzDUFhYmKpWrZptx61YsaJ2794tu93uFCh++eUXx/a0Njo6WufOnXMKZfv3789U7Xa7XXv37tUtt9ySqbqmTZum06dPa968eSpWrJimT5/uOJYklSpVSq1bt87Usa5XyZIl5e/vr9TU1Bw/FwAgY1yjBAC4qm7dusnd3V3PP/+8o8cjjWEYOnXq1HUdt2PHjjp27Jg++eQTx7qUlBS9/vrr8vPzU4sWLRz7paSkaO7cuY79UlNT9frrr1/zHF27dpWbm5smTpzo6KW5vPYrefPNN9W9e3fNmDFDkydPliS1a9dOAQEBeuGFF5ScnJzuNSdPnrxmPZnl7u6ue++9V0uWLNFPP/2Uo+cCAGSMHiUAwFVVrlxZkydP1rhx43To0CF17dpV/v7+iouL07JlyzR48GCNGTMmy8cdPHiw3nzzTUVGRmrnzp0KDQ3Vp59+qk2bNmnmzJny9/eXJHXu3FlNmzbV2LFjdejQIdWoUUNLly5VQkLCNc8RHh6up59+WpMmTVKzZs3UrVs3eXl5afv27SpbtqymTp2a4evc3Ny0YMECJSQkaPz48QoODtaQIUM0d+5c9e3bV/Xq1VOvXr1UsmRJHT58WF988YWaNm2q2bNnZ/lzuJIXX3xR69evV6NGjTRo0CDVqFFDp0+f1vfff69169bp9OnT2XYuAEB6BCUAwDWNHTtWVatW1auvvqrnn39eklS+fHm1bdtWd99993Ud08fHRzExMRo7dqzee+89nT17VtWqVVNUVJQiIyMd+7m5uWnFihUaNWqUPvzwQ9lsNt19992aPn266tate83zTJw4UWFhYXr99df19NNPy9fXV3Xq1FHfvn2v+jpPT08tW7ZMrVu31vDhwxUUFKTevXurbNmyevHFF/Xyyy/r4sWLKleunJo1a5apmf6yonTp0tq2bZsmTpyopUuXas6cOSpevLhq1qypadOmZeu5AADp2YyrjT0AAAAAgEKIa5QAAAAAwAVBCQAAAABcEJQAAAAAwAVBCQAAAABcEJQAAAAAwAVBCQAAAABcEJQAoICx2WyaMGFCjp6jZcuWatmyZbYdz7XmCRMmyGaz6e+//86R48+fP182m02HDh1yrAsNDdVdd92VLefLDocOHZLNZtP8+fOtLgUACiWCEoA8K+2L4iuvvGJ1KZKkyMhI2Ww2BQQE6MKFC+m2x8bGymaz5amarRYTEyObzaZPP/3U6lLwr/379+vRRx9VkyZN5O3tnS4wXkva73hGjzZt2qTb/+DBg+rdu7dKlSolHx8fValSRU8//XS6/WbPnq2IiAh5eXmpXLlyeuyxx3T+/PkbeavXZe7cubrvvvtUoUIF2Ww2p5sfX+7o0aMaO3as7rjjDvn7+8tmsykmJiZXawWQszysLgAA8hMPDw8lJibq888/V48ePZy2LViwQN7e3kpKSrKoOtOFCxfk4ZG//vOe2zX37dt
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACkiklEQVR4nOzde3yP9f/H8edn5/Nmc9jM2Jgz9XXOIVMUKpGEkiyibxT6UdEXCaWUom9fdKRyKIdIR8nXFAlJoiRkJufTxuy8z/X742qfr48Nw+baZ3vcb7fdXtt1XZ/ren0++3D7PHdd1/ttMwzDEAAAAADAwc3qBgAAAACgpCEoAQAAAMB5CEoAAAAAcB6CEgAAAACch6AEAAAAAOchKAEAAADAeQhKAAAAAHAeghIAAAAAnIegBAAAAADnISgBAGCRhIQE2Ww2JSQkWN0KAOA8BCUAwBWbM2eObDab48vDw0ORkZGKj4/XgQMH8m3frl07p+29vLwUExOjQYMGaf/+/Rfd97lfo0aNcmwXHR2tO+64I9+xPvjgA7m7u6tTp07KyMgo+icPACjVPKxuAADg+iZMmKCYmBhlZGTohx9+0Jw5c7R27Vpt375dPj4+TttWqVJFkydPliRlZWXpt99+06xZs7RixQrt2LFDfn5+Be77XA0aNLhoP/PmzVN8fLw6dOigZcuW5esBAIBLISgBAK5a586d1bRpU0nSQw89pPLly+vFF1/U8uXL1bNnT6dtg4ODdf/99zsti4mJ0aOPPqp169bplltuueC+C+PDDz9Uv379dPPNN+uTTz4pkpCUlpaWL8ABAEo3Lr0DABS5G2+8UZK0Z8+eQm0fHh4uSfLwuLq/3y1cuFD333+/2rVrp+XLl+cLSXPnzlWTJk3k6+ur0NBQ9e7dO98lf+3atVODBg20efNmtW3bVn5+fnr66aeVmJgom82ml19+WW+++aZq1Kghb29vNWvWTJs2bcrXy++//64ePXooNDRUPj4+atq0qZYvX35Vzw8AcO1wRgkAUOQSExMlSeXKlcu3Ljc3V8ePH5ckZWdna8eOHXrmmWcUGxur1q1b59s+JSXFsX2e8uXL59tuyZIl6tOnj9q2batPP/1Uvr6+Tuufe+45jR07Vj179tRDDz2kY8eO6d///rfatm2rLVu2KCQkxLHtiRMn1LlzZ/Xu3Vv333+/KlWq5Fg3f/58nTlzRg8//LBsNpumTJmi7t27688//5Snp6ck6ddff1Xr1q0VGRmpUaNGyd/fXwsXLlS3bt20ZMkS3XXXXYV7IQEAliEoAQCuWl6YycjI0IYNG/Tss8/K29u7wEEWfv/9d1WoUMFpWd26dfX111/Ly8sr3/YdOnTIt8wwDKeft2zZohUrVqhNmzb67LPP8oWkffv26ZlnntGkSZP09NNPO5Z3795djRo10owZM5yWHz58WLNmzdLDDz/sWJYX/pKSkrRr1y5HCKxdu7a6du2qFStWOJ7vsGHDVLVqVW3atEne3t6SpMGDB6tNmzZ66qmnCEoA4AIISgCAq3Z+mImOjtbcuXNVpUqVfNtGR0frrbfekiTl5ORo586dmjJlijp37qzvvvsuX4j6z3/+o1q1al30+CdPnlROTo6qVKmSLyRJ0scffyy73a6ePXs6nZ0KDw9XzZo1tXr1aqeg5O3trQcffLDAY/Xq1cvpTFneZYZ//vmno5f//ve/mjBhgs6cOaMzZ844tu3YsaOeeeYZHThwQJGRkRd9TgAAaxGUAABXLS/MpKSk6N1339W3337rOJNyPn9/f6dg1alTJ7Vp00ZNmzbVCy+8oKlTpzpt37x580sO5tC+fXtVrVpVM2fOVGhoqKZPn+60fteuXTIMQzVr1izw8XmXzOWJjIws8OyWJFWtWtXp57zQdOrUKUnS7t27ZRiGxo4dq7Fjxxa4j6NHjxKUAKCEIygBAK7auWGmW7duatOmje677z7t3LlTAQEBl3x8kyZNFBwcrG+//faKe3j99dd16tQpvfbaaypXrpzGjx/vWGe322Wz2fTll1/K3d0932PP77Ggs1J5Cnq89L/LAe12uyRp5MiR6tixY4HbxsbGXvS5AACsR1ACABQpd3d3TZ48WTfddJNef/11p8lhLyY3N1epqalXfFw3Nze9//77SklJ0bPPPqvQ0FANHTpUklSjRg0ZhqGYmJhLXsZ3tapXry7JPEtV0P1VAADXwPDgAIAi165dOzVv3lzTpk1TRkbGJbdfvXq1UlNTdf3111/VcT09PbV48WK1bt1aw4cP1wcffCDJHLTB3d1dzz77bL6BIAzD0IkTJ67quOeqWLGi2rVrpzfeeEOHDh3Kt/7YsWNFdiwAQPHhjBIAoFg88cQTuueeezRnzhz985//dCxPSUnR3LlzJf1vMIeZM2fK19e30GefLsbPz0+ff/654uLi1L9/fwUHB+vOO+/UpEmTNHr0aCUmJqpbt24KDAzU3r17tXTpUg0aNEgjR4686mPn+c9//qM2bdqoYcOGGjhwoKpXr64jR45o/fr1+uuvv7R169YiOxYAoHgQlAAAxaJ79+6qUaOGXn75ZQ0cONBxb89ff/2lvn37SpJsNpvKlSunuLg4PfPMM/rHP/5RJMcODg52DBfeq1cvffnllxo1apRq1aqlV199Vc8++6wkKSoqSrfeeqvuvPPOIjlunnr16unHH3/Us88+qzlz5ujEiROqWLGiGjVqpHHjxhXpsQAAxcNmnH8NAgAAAACUcdyjBAAAAADnISgBAAAAwHkISgAAAABwHoISAAAAAJyHoAQAAAAA5yEoAQAAAMB5CEoAUMrYbDaNHz++WI/Rrl07tWvXrsj2d37P48ePl81m0/Hjx4tl/3PmzJHNZlNiYqJjWXR0tO64444iOV5RSExMlM1m05w5c6xuBQDKJIISgBIr74Piyy+/bHUrkqT4+HjZbDYFBQUpPT093/pdu3bJZrOVqJ6tlpCQIJvNpsWLF1vdCi7glltukc1m06OPPnrZj01OTlbFihUL9Tt+7rnnZLPZ1KBBgyLbZ3Gw2+2aMmWKYmJi5OPjo+uuu04LFizIt91bb72luLg4VapUSd7e3oqJidGDDz7oFL4BuDYPqxsAAFfi4eGhtLQ0ffrpp+rZs6fTunnz5snHx0cZGRkWdWdKT0+Xh4dr/fd+rXvu27evevfuLW9v72t2zJLo448/1vr166/48ePGjVNaWtolt/vrr7/0/PPPy9/fv8j2WVz+9a9/6YUXXtDAgQPVrFkzffLJJ7rvvvtks9nUu3dvx3ZbtmxRTEyM7rzzTpUrV0579+7VW2+9pc8++0xbt25V5cqVLXsOAIoGZ5QA4DJ4e3urffv2Bf6Fef78+br99tuL9HgZGRmy2+2X9RgfHx+XC0rXumd3d3f5+PjIZrNds2OWNBkZGRoxYoSeeuqpK3r89u3bNXPmzEI9fuTIkbrhhhvUtGnTItvn5Ro/fryio6Mvus2BAwc0depUDRkyRG+++aYGDhyoTz/9VDfeeKOeeOIJ5ebmOradMWOG5syZoxEjRqh///6aOHGiPv/8cx0/flzvv/9+kfcP4NojKAFweUePHtWAAQNUqVIl+fj46Prrr9d7772Xb7sTJ06ob9++CgoKUkhIiPr166etW7de9n0g9913n7788kslJyc7lm3atEm7du3Sfffdl2/7kydPauTIkWrYsKECAgIUFBSkzp07a+vWrU7b5V2m9uGHH2rMmDGKjIyUn5+fTp8+LUlatGiR6tWrJx8fHzVo0EBLly5VfHx8vg9/F7rfZ/fu3YqPj1dISIiCg4P14IMP5vvL/ezZs3XzzTerYsWK8vb2Vr169TRz5sxCvzZXqjD3Ve3bt0+xsbFq0KCBjhw5Ism8TGv48OGKioqSt7e3YmNj9eKLL14yXBZ0j1KetWvXqnnz5vLx8VH16tUL/ND7559/6p577lFoaKj8/Px0ww036PPPP8+3XWHfm8nJyYqPj1dwcLDjvXnu+ytPdna2fv/9dx06dOiiz68wpkyZIrvdrpEjR17R44cNG6a77rpLN95440W3+/bbb7V48WJ
     },
     "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": 20,
   "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": 21,
   "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": 22,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_12816/3712199621.py:45: RuntimeWarning: invalid value encountered in sqrt\n",
      "  y_star[:, 0] - 1.96 * np.sqrt(np.diag(var)),\n",
      "/tmp/ipykernel_12816/3712199621.py:46: 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/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACd2klEQVR4nOzdeVxU9f7H8feIsm/ixhKCigsued1Tc81yy7QyN1JR066792qLtyzTzMo07Zpa1wIzcss0M8vUwoW0xXKpyNAgzT0VEBVEOL8/5ufkCBxRgUF8PR+Pecyd8/3OOZ85jDfefL/neyyGYRgCAAAAAOSqlKMLAAAAAIDijNAEAAAAACYITQAAAABggtAEAAAAACYITQAAAABggtAEAAAAACYITQAAAABggtAEAAAAACYITQAAAABggtAEACgQSUlJslgsio6OdnQpJRrnGQCKHqEJAHBN0dHRslgs+v777x1dSoFo27at6tatm2P7pk2b5O7uroYNG+r06dMOqAwAUByVdnQBAICSISQkRBcuXFCZMmUcXcoN+fLLL9WtWzfVrFlTGzdulJ+fn6NLAgAUE4w0AQAKhMVikaurq5ycnBxdSp7Onz+f6/bNmzerW7duqlGjRoEFpnPnzt30PgAAxQOhCQBQIHK71iYyMlKenp46fPiwevToIU9PT1WoUEETJkxQVlaW3fuzs7M1e/Zs1alTR66urqpUqZIef/xxnTlzxq7fxx9/rK5duyowMFAuLi6qVq2apk6dmmN/l6fg7dy5U61bt5a7u7v+85//5Kh769at6tq1q8LCwrRx40aVK1fOrv2zzz5Tq1at5OHhIS8vL3Xt2lU///yzXZ/Ln/PAgQPq0qWLvLy8FBERIckaJkeNGqXVq1erbt26cnFxUZ06dfT555/nqOXw4cMaPHiwKlWqZOv37rvvXvvkAwAKFdPzAACFKisrSx07dlSzZs302muvaePGjZo5c6aqVaum4cOH2/o9/vjjio6O1qBBgzRmzBglJiZq7ty5+vHHHxUXF2eb9hcdHS1PT0/9+9//lqenp7788ks999xzSk1N1YwZM+yOferUKXXu3Fl9+vTRo48+qkqVKtm1x8XFqUuXLqpSpYo2bdqk8uXL27UvXrxYAwcOVMeOHfXKK6/o/Pnzmj9/vu6++279+OOPCg0NtfW9dOmSOnbsqLvvvluvvfaa3N3dbW3btm3TRx99pBEjRsjLy0tvvPGGHn74YR08eNAW0o4fP6677rrLFrIqVKigzz77TEOGDFFqaqrGjRtXED8OAMCNMAAAuIaoqChDkvHdd9/l2ScxMdGQZERFRdm2DRw40JBkTJkyxa5vgwYNjEaNGtleb9261ZBkxMTE2PX7/PPPc2w/f/58jmM//vjjhru7u5Genm7b1qZNG0OSsWDBghz927RpY/j5+RleXl5GnTp1jBMnTuToc/bsWcPX19cYOnSo3fZjx44ZPj4+dtsvf86nn346x34kGc7Ozsb+/ftt23bv3m1IMv773//atg0ZMsQICAgw/vrrL7v39+nTx/Dx8bF97tzOMwCgcDE9DwBQ6P75z3/avW7VqpV+//132+sVK1bIx8dH9957r/766y/bo1GjRvL09NRXX31l6+vm5mb732fPntVff/2lVq1a6fz58/r111/tjuPi4qJBgwblWtO5c+d09uxZVapUSd7e3jnaN2zYoOTkZPXt29euJicnJzVr1syupsuuHDm7UocOHVStWjXb6zvvvFPe3t62c2AYhlauXKlu3brJMAy743Xs2FEpKSn64Ycfct03AKDwMT0PAFCoXF1dVaFCBbttZcuWtbtWKSEhQSkpKapYsWKu+zhx4oTtf//888969tln9eWXXyo1NdWuX0pKit3roKAgOTs757rPsLAwDRgwQE899ZT69u2rFStW2C1ikZCQIElq3759ru+/OmiVLl1ad9xxR659K1eunGPblefg5MmTSk5O1ttvv6233347131ceQ4AAEWL0AQAKFT5WU0vOztbFStWVExMTK7tl0NXcnKy2rRpI29vb02ZMkXVqlWTq6urfvjhBz311FPKzs62e9+Vo1K5efLJJ3Xq1Cm9+uqrGjp0qN555x1ZLBZbTZL1uiZ/f/8c7y1d2v4/oS4uLipVKvcJHHmdA8Mw7I716KOPauDAgbn2vfPOO00/CwCg8BCaAAAOV61aNW3cuFEtW7Y0DTqxsbE6deqUPvroI7Vu3dq2PTEx8YaP/corr+j06dNauHChypYtq5kzZ9pqkqSKFSuqQ4cON7z//KhQoYK8vLyUlZVV6McCAFw/rmkCADhcr169lJWVpalTp+Zou3TpkpKTkyX9PWJzeYRGki5evKh58+bd1PHfeust9ezZU7NmzdKLL74oSerYsaO8vb310ksvKTMzM8d7Tp48eVPHvJKTk5MefvhhrVy5Uj/99FOhHgsAcP0YaQIA5Nu7776b6/2Fxo4de1P7bdOmjR5//HFNnz5du3bt0n333acyZcooISFBK1as0Jw5c9SzZ0+1aNFCZcuW1cCBAzVmzBhZLBYtXrzYLkTdiFKlSikmJkYpKSmaNGmS/Pz8NGLECM2fP1/9+/dXw4YN1adPH1WoUEEHDx7Up59+qpYtW2ru3Lk3ddwrvfzyy/rqq6/UrFkzDR06VLVr19bp06f1ww8/aOPGjTp9+nSBHQsAcH0ITQCAfJs/f36u2yMjI2963wsWLFCjRo301ltv6T//+Y9Kly6t0NBQPfroo2rZsqUkqVy5clq7dq3Gjx+vZ599VmXLltWjjz6qe+65Rx07dryp4zs7O2vVqlXq0KGDRo8eLV9fX/Xr10+BgYF6+eWXNWPGDGVkZCgoKEitWrXKc1W+G1WpUiV9++23mjJlij766CPNmzdP5cqVU506dfTKK68U6LEAANfHYtzsn+cAAAAAoATjmiYAAAAAMEFoAgAAAAAThCYAAAAAMEFoAgAAAAAThCYAAAAAMEFoAgAAAAAThCYAKKEsFosmT55cqMdo27at2rZtW2D7u7rmyZMny2Kx6K+//iqU/UdHR8tisSgpKcm2LTQ0VPfff3+BHK8gJCUlyWKxKDo62tGlAMBti9AEoNi7/Evja6+95uhSJFlv5GqxWOTt7a0LFy7kaE9ISJDFYilWNTtabGysLBaLPvzwQ0eXgv+3b98+/etf/1KLFi3k6uqaIzxey7fffqsRI0aoUaNGKlOmjCwWS559U1JS9OSTT6p69epyc3NTSEiIhgwZooMHD9r1uxySr364urre6Me8YcuWLdOjjz6q6tWry2Kx5PnHgbS0ND3//PPq1KmT/Pz8CLhACVXa0QUAwK2odOnSOn/+vD755BP16tXLri0mJkaurq5KT093UHVWFy5cUOnSt9b/zRd1zf3791efPn3k4uJSZMcsLrZv36433nhDtWvXVnh4uHbt2nVd71+3bp0WLlyoO++8U1WrVtVvv/2Wa7/s7Gzde++9+uWXXzRixAjVqFFD+/fv17x587R+/XrFx8fLy8vL7j3z58+Xp6en7bWTk9N1f76bNX/+fO3cuVNNmjTRqVOn8uz3119/acqUKapcubLq16+v2NjYoisSQJG5tf5rCgDFhIuLi1q2bKklS5bkCE0ffPCBunbtqpUrVxbY8dLT0+Xs7KxSpfI/QcARf52/WUVds5OTk0N+IS8OHnjgASUnJ8vLy0uvvfbadYem4cOH66mnnpKbm5tGjRqVZ2jasWOHvvvuO82dO1cjR460ba9Zs6YGDx6sjRs36sEHH7R7T8+ePVW+fPnr/kz5ERsbq3bt2ikxMVGhoaF59lu8eLGCgoJUqlQp1a1bN89+AQEBOnr0qPz9/fX999+rSZMmhVA1AEdjeh6AEuPEiRMaMmSIKlWqJFdXV9WvX1+LFi3K0e/UqVPq37+/vL295evrq4EDB2r37t3XPa2mX79++uyzz5ScnGzb9t133ykhIUH9+vXL0f/06dOaMGGC6tWrJ09PT3l7e6tz587avXu3Xb/LU9mWLl2qZ599VkFBQXJ3d1dqaqokacWKFapdu7ZcXV1Vt25drVq
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADJ00lEQVR4nOzdd1xV5R8H8M9lXPYQBBkyFBBRsXKmpLhx5siNAs5ScRSWWlamOdIclDlKBSvcppn5k5TEQY60XIWIBuFAxQGI7HvP748TFy5LQOBc4PN+ve7reTjzcy/oiy/nOc+RCYIggIiIiIiIiFS0pA5ARERERESkaVgoERERERERFcJCiYiIiIiIqBAWSkRERERERIWwUCIiIiIiIiqEhRIREREREVEhLJSIiIiIiIgKYaFERERERERUCAslIiIiIiKiQlgoERFRsbp06YIuXbpIHaNShIaGQiaTIT4+vtz7BgQEwNnZudIzSeVFPgsiorqEhRIRUS2R9wtw3ktfXx9NmjRBYGAg7t+/L3W8WkUmkyEwMLDI8iVLlkAmk2H8+PFQKpUSJCMiosqiI3UAIiKqXAsXLkSjRo2QmZmJU6dOYf369Th06BCuXr0KQ0NDqeNJYuzYsRg5ciT09PSq7BzLli3DBx98AH9/f2zatAlaWvxbJBFRTcZCiYiolunTpw/atGkDAJg4cSIsLS2xatUq/Pjjjxg1apTE6aShra0NbW3tKjv+ihUrMG/ePPj5+WHLli0vXCQJgoDMzEwYGBhUUkIiIiov/rmLiKiW69atGwAgLi4OAJCbm4tFixbBxcUFenp6cHZ2xvvvv4+srKwSj5GWlgYjIyPMnDmzyLrbt29DW1sbS5cuBZA/BDAqKgrvvPMOrKysYGRkhMGDByMpKanI/uvWrUPz5s2hp6cHOzs7TJs2DcnJyWrbdOnSBS1atMDly5fh7e0NQ0NDuLq6Ys+ePQCA48ePo3379jAwMIC7uzuOHj2qtn9x9+X8+OOP6NevH+zs7KCnpwcXFxcsWrQICoXi+R9qAatWrcJ7772HMWPGICQkRK1IUiqVWLNmDZo3bw59fX00aNAAb775Jp48eaJ2DGdnZ/Tv3x/h4eFo06YNDAwMsHHjRkRGRkImk2HXrl1YvHgxGjZsCH19fXTv3h03btwokuXs2bPo3bs3zMzMYGhoCG9vb0RFRZXr/RARkYiFEhFRLXfz5k0AgKWlJQDxKtNHH32EVq1aYfXq1fD29sbSpUsxcuTIEo9hbGyMwYMHY+fOnUUKie3bt0MQBPj6+qotnz59Oi5duoSPP/4YU6ZMwU8//VTkvp4FCxZg2rRpsLOzw8qVK/HGG29g48aN6NWrF3JyctS2ffLkCfr374/27dtj+fLl0NPTw8iRI7Fz506MHDkSffv2xbJly/Ds2TMMHToUT58+LfVzCQ0NhbGxMd555x0EBwejdevW+OijjzB37tzSP9ACgoODERQUhNGjRyM0NLTIlaQ333wT7777Lry8vBAcHIxx48YhLCwMPj4+Rd5fTEwMRo0ahZ49eyI4OBgvv/yyat2yZcuwb98+zJ49G/PmzcOZM2eKfN6//vorOnfujNTUVHz88cdYsmQJkpOT0a1bN5w7d67M74mIiP4jEBFRrRASEiIAEI4ePSokJSUJt27dEnbs2CFYWloKBgYGwu3bt4WLFy8KAISJEyeq7Tt79mwBgPDrr7+qlnl7ewve3t6qr8PDwwUAwv/+9z+1fVu2bKm2XV6OHj16CEqlUrX87bffFrS1tYXk5GRBEAThwYMHglwuF3r16iUoFArVdmvXrhUACFu2bFHLAkDYtm2batm1a9cEAIKWlpZw5syZIjlDQkKKZIqLi1MtS09PL/IZvvnmm4KhoaGQmZmpWubv7y84OTmpbQdAcHJyEgAIo0aNEnJzc4sc6+TJkwIAISwsTG354cOHiyzPO9bhw4fVtj127JgAQPDw8BCysrJUy4ODgwUAwpUrVwRBEASlUim4ubkJPj4+ap95enq60KhRI6Fnz56lfhZERFQUrygREdUyPXr0gJWVFRwcHDBy5EgYGxtj3759sLe3x6FDhwAA77zzjto+QUFBAICff/651OPa2dkhLCxMtezq1au4fPkyxowZU2T7yZMnQyaTqb7u1KkTFAoF/v33XwDA0aNHkZ2djVmzZqldiZk0aRJMTU2LZDE2Nla76uXu7g5zc3N4eHigffv2quV5/X/++afE9wJA7f6fp0+f4uHDh+jUqRPS09Nx7dq1UvcFoJpJsFGjRsXe/7R7926YmZmhZ8+eePjwoerVunVrGBsb49ixY2rbN2rUCD4+PsWea9y4cZDL5aqvO3XqpPYeL168iNjYWIwePRqPHj1SnevZs2fo3r07Tpw4wVn4iIjKiZM5EBHVMl999RWaNGkCHR0dNGjQAO7u7qpC5N9//4WWlhZcXV3V9rGxsYG5ubmqiCmOlpYWfH19sX79eqSnp8PQ0BBhYWHQ19fHsGHDimzv6Oio9nW9evUAQHV/Tt653N3d1baTy+Vo3LhxkSwNGzZUK7wAwMzMDA4ODkWWFTxPSf766y/Mnz8fv/76K1JTU9XWpaSklLovAPj7++Pu3btYsmQJ6tevj7ffflttfWxsLFJSUmBtbV3s/g8ePFD7ulGjRiWe63mfZWxsrCpTSVJSUlT7ERHR87FQIiKqZdq1a6ea9a4khQuOsvLz88OKFSuwf/9+jBo1Ctu2bUP//v1VxUlBJc0yJwhChc5d0vEqcp7k5GR4e3vD1NQUCxcuhIuLC/T19fHHH39gzpw5Zbr6oqOjg127dqF3794ICgqCubk5xo0bp1qvVCphbW2tdgWuICsrK7WvS5vh7nnvMS/vihUr1O5tKsjY2LjE4xMRUVEslIiI6hAnJycolUrExsbCw8NDtfz+/ftITk6Gk5NTqfu3aNECr7zyCsLCwtCwYUMkJCTgyy+/rHAWQJzEoHHjxqrl2dnZiIuLQ48ePSp03LKIjIzEo0eP8MMPP6Bz586q5XkzA5aVvr4+Dhw4gK5du2LSpEkwNzfH4MGDAQAuLi44evQovLy8qnyabxcXFwCAqalplX5uRER1Ce9RIiKqQ/r27QsAWLNmjdryVatWAQD69ev33GOMHTsWv/zyC9asWQNLS0v06dOnQll69OgBuVyOL774Qu3qz+bNm5GSklKmLBWVd4Wm4Hmzs7Oxbt26ch/L1NQUhw8fhqurK0aNGoWIiAgAwPDhw6FQKLBo0aIi++Tm5haZAv1FtG7dGi4uLvj888+RlpZWZH1x07ITEVHpeEWJiKgOeemll+Dv74+vv/5aNfzs3Llz2Lp1KwYNGoSuXbs+9xijR4/Ge++9h3379mHKlCnQ1dWtUBYrKyvMmzcPn3zyCXr37o3XX38dMTExWLduHdq2bVvsBBGVpWPHjqhXrx78/f0xY8YMyGQyfPfddxUeFmhlZYUjR47Ay8sLgwYNQkREBLy9vfHmm29i6dKluHjxInr16gVdXV3ExsZi9+7dCA4OxtChQyvl/WhpaWHTpk3o06cPmjdvjnHjxsHe3h537tzBsWPHYGpqip9++qlSzkVEVFewUCIiqmM2bdqExo0bIzQ0FPv27YONjQ3mzZuHjz/+uEz7N2jQAL169cKhQ4cwduzYF8qyYMECWFlZYe3atXj77bdhYWGByZMnY8mSJRUuwMrC0tISBw8eRFBQEObPn4969ephzJgx6N69e4kzzz2Pg4MDfvnlF3Tq1Al9+vTBiRMnsGHDBrRu3RobN27E+++/Dx0dHTg7O2PMmDHw8vKq1PfUpUsXnD59GosWLcLatWuRlpYGGxsbtG/fHm+++WalnouIqC6QCRX98xkREdVZgwcPxpUrV3Djxg2poxAREVUJ3qNERETlkpiYiJ9//vmFryYRERFpMg69IyKiMomLi0NUVBQ2bdoEXV1dDuciIqJajVeUiIioTI4fP46xY8ciLi4OW7duhY2NjdSRiIiIqgzvUSIiIiIiIiqEV5SIiIiIiIgKYaFERERERERUCAslIqJaRiaTYcGCBVV6ji5duqBLly6
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADMeUlEQVR4nOzdeVxU9f7H8dcBZJNNcUEQRMUFU8v9qrlU5lamlZlFKml6c0m9aand7JpmlqZpmVrX0sq9XDKzzLjhwk+zLJeKDBWuG+4CIjuc3x8Tc2VRUYFheT8fj3mcmTlnzvnMMBZvvp/zPYZpmiYiIiIiIiJiZWfrAkREREREREoaBSUREREREZFcFJRERERERERyUVASERERERHJRUFJREREREQkFwUlERERERGRXBSUREREREREclFQEhERERERyUVBSUREREREJBcFJRERKRadO3emc+fOhba/mJgYDMNg6dKl1uemTJmCYRiFdozSIL/PQUREbp+CkohIObZ06VIMw7DenJ2dqV+/PqNGjeLMmTO2Lq/YdO7cmcaNG+d5PiwsDFdXV5o3b87FixdtUJmIiNiKg60LEBER25s6dSq1a9cmJSWFnTt3snDhQjZv3syvv/6Kq6troRzj22+/LZT9XM/LL7/MxIkTC2Vf//nPf+jVqxcNGjTgu+++o3LlyoWyXxERKR0UlEREhB49etCyZUsAnnnmGby9vZkzZw5ffPEFTzzxxG3tOykpCVdXVxwdHQuj1OtycHDAweH2/9e2bds2evXqRf369QstJF25coWKFSve9n5ERKR4qPVORETyuPfeewGIjo62Prds2TJatGiBi4sLlStXpn///hw/fjzH67Jb2Pbu3UvHjh1xdXXlpZdesq7LfY7S2bNnGTJkCNWrV8fZ2Zk777yTjz/+OE89cXFxhIaG4unpiZeXF4MGDSIuLi7Pdtc6R2nZsmW0bt0aV1dXKlWqRMeOHa85wrVjxw4eeOABgoKC+O677/D29s6x/uuvv6ZDhw5UrFgRd3d3HnjgAX777bcc24SGhuLm5saRI0fo2bMn7u7uhISEAGAYBqNGjWLDhg00btwYJycn7rjjDr755ps8tZw8eZLBgwdTvXp163YfffRRvnWLiEjh0oiSiIjkceTIEQBrSJg+fTqTJ0+mX79+PPPMM5w7d453332Xjh078ssvv+Dl5WV97YULF+jRowf9+/fnqaeeonr16vkeIzk5mc6dO3P48GFGjRpF7dq1+eyzzwgNDSUuLo4xY8YAYJomvXv3ZufOnTz77LMEBwezfv16Bg0aVKD38uqrrzJlyhTatWvH1KlTcXR05IcffuA///kPXbt2zbFtREQEPXv2pHbt2oSFhVGlSpUc6z/99FMGDRpEt27dePPNN0lKSmLhwoXcfffd/PLLLwQGBlq3zcjIoFu3btx999289dZbOVoYd+7cybp16xgxYgTu7u688847PProoxw7dsz6mZ85c4a//e1v1mBVtWpVvv76a4YMGUJCQgJjx44t0PsXEZFbZIqISLm1ZMkSEzC/++4789y5c+bx48fNVatWmd7e3qaLi4t54sQJMyYmxrS3tzenT5+e47UHDx40HRwccjzfqVMnEzAXLVqU51idOnUyO3XqZH08d+5cEzCXLVtmfS4tLc1s27at6ebmZiYkJJimaZobNmwwAXPmzJnW7TIyMswOHTqYgLlkyRLr8//617/Mq//XFhUVZdrZ2ZkPP/ywmZmZmaOerKysHLVVrlzZdHd3N++44w7z7Nmzeeq/fPmy6eXlZQ4dOjTH86dPnzY9PT1zPD9o0CATMCdOnJhnP4Dp6OhoHj582Prc/v37TcB89913rc8NGTLErFGjhnn+/Pkcr+/fv7/p6elpJiUlmaZpmtHR0Xk+BxERuX1qvRMREbp06ULVqlXx9/enf//+uLm5sX79evz8/Fi3bh1ZWVn069eP8+fPW28+Pj7Uq1eP77//Pse+nJycePrpp294zM2bN+Pj45PjHKgKFSowevRoEhMT2bZtm3U7BwcHhg8fbt3O3t6e55577obH2LBhA1lZWbzyyivY2eX8X17uFr0rV65w+fJlqlevjoeHR559bd26lbi4OJ544okcn4O9vT1t2rTJ8zkAOWq+WpcuXahbt671cdOmTfHw8ODo0aOAZRRt7dq19OrVC9M0cxyvW7duxMfH8/PPP9/w/YuIyK1T652IiPDee+9Rv359HBwcqF69Og0aNLAGi6ioKEzTpF69evm+tkKFCjke+/n5FWjihv/+97/Uq1cvT4AJDg62rs9e1qhRAzc3txzbNWjQ4IbHOHLkCHZ2djRq1OiG2wYFBTFw4EAmTJjAE088wWeffYa9vb11fVRUFPC/87dyyx2uHBwcqFmzZr7bBgQE5HmuUqVKXLp0CYBz584RFxfHBx98wAcffJDvPs6ePXvD9yQiIrdOQUlERGjdurV11rvcsrKyMAyDr7/+OkdwyJY7wLi4uBRJjcXhxRdf5MKFC8ycOZOhQ4fy4YcfWkeesrKyAMt5Sj4+Pnlem3u2PScnpzwhMFt+nyNYRpKuPtZTTz11zXOxmjZtWoB3JCIit0pBSURErqtu3bqYpknt2rWpX79+oe23Vq1aHDhwgKysrByB4o8//rCuz16GhYWRmJiYI5QdOnSoQLVnZWXx+++/c9dddxWorjfffJOLFy+yePFiKlWqxOzZs637AqhWrRpdunQp0L5uVdWqVXF3dyczM7PIjyUiIvnTOUoiInJdjzzyCPb29rz66qvWEY9spmly4cKFW9pvz549OX36NKtXr7Y+l5GRwbvvvoubmxudOnWybpeRkcHChQut22VmZvLuu+/e8Bh9+vTBzs6OqVOnWkdprq79Wt5//3369u3LnDlzeO211wDo1q0bHh4evP7666Snp+d5zblz525YT0HZ29vz6KOPsnbtWn799dciPZaIiORPI0oiInJddevW5bXXXmPSpEnExMTQp08f3N3diY6OZv369QwbNozx48ff9H6HDRvG+++/T2hoKHv37iUwMJDPP/+ciIgI5s6di7u7OwC9evWiffv2TJw4kZiYGBo1asS6deuIj4+/4TGCgoL45z//ybRp0+jQoQOPPPIITk5O/Pjjj/j6+jJjxox8X2dnZ8fy5cuJj49n8uTJVK5cmREjRrBw4UIGDBhA8+bN6d+/P1WrVuXYsWN89dVXtG/fnvnz59/053Atb7zxBt9//z1t2rRh6NChNGrUiIsXL/Lzzz/z3XffcfHixUI7loiI5KWgJCIiNzRx4kTq16/P22+/zauvvgqAv78/Xbt25aGHHrqlfbq4uBAeHs7EiRP5+OOPSUhIoEGDBixZsoTQ0FDrdnZ2dmzcuJGxY8eybNkyDMPgoYceYvbs2TRr1uyGx5k6dSq1a9fm3Xff5Z///Ceurq40bdqUAQMGXPd1jo6OrF+/ni5duvDcc8/h5eXFk08+ia+vL2+88QazZs0iNTUVPz8/OnToUKCZ/m5G9erV2bNnD1OnTmXdunUsWLAAb29v7rjjDt58881CPZaIiORlmNfrPRARERERESmHdI6SiIiIiIhILgpKIiIiIiIiuSgoiYiIiIiI5KKgJCIiIiIikouCkoiIiIiISC4KSiIiIiIiIrkoKImIlDGGYTBlypQiPUbnzp3p3Llzoe0vd81TpkzBMAzOnz9fJPtfunQphmEQExNjfS4wMJAHH3ywUI5XGGJiYjAMg6VLl9q6FBGRcklBSURKrOxfFN966y1blwJAaGgohmHg4eFBcnJynvVRUVEYhlGiara18PBwDMPg888/t3Up8pd169bx+OOPU6dOHVxdXWnQoAHjxo0jLi6uwPvIyspi4cKF3HXXXbi4uODt7c29997L/v37rdtk//vN77Zq1aoc+9uzZw8jRoygRYsWVKhQAcMwCuvt3pKNGzfSvHlznJ2dCQgI4F//+hcZGRl5ttu6dSt33303rq6uVKpUib59++YI3yJSujnYugARkdLEwcGBpKQkvvzyS/r165dj3fLly3F2diYlJcVG1VkkJyfj4FC6/vN
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x1000 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAANsCAYAAABlL1jGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADLrklEQVR4nOzdeZyNdf/H8dc1+z5jbLMYBmPNclsLWSqFNpKoppiI+w7hjkKlZC1S9CupFC1jK5HkTpIhk1SyVZOGSNnXYcx+5vr9cTUnYwaDGdcs7+fjcR7fM+e6znV9zpkzdd6+y2WYpmkiIiIiIiIiTi52FyAiIiIiIlLcKCiJiIiIiIicQ0FJRERERETkHApKIiIiIiIi51BQEhEREREROYeCkoiIiIiIyDkUlERERERERM6hoCQiIiIiInIOBSUREREREZFzKCiJiIjYJC4uDsMwiIuLs7sUERE5h4KSiIhctrlz52IYhvPm5uZGeHg4MTEx7Nu3L8/+HTp0yLW/h4cH1atXZ8CAAfz5558XPPbZt1GjRjn3i4yM5Pbbb89zrvfffx9XV1c6d+5MWlpa4b94EREp1dzsLkBEREq+cePGUb16ddLS0vj222+ZO3cu69ev56effsLLyyvXvlWqVGHy5MkAZGRk8MsvvzBr1ixWrlxJQkICPj4++R77bA0aNLhgPbGxscTExNCxY0eWLl2apwYREZGLUVASEZEr1qVLF5o3bw7Aww8/TIUKFXjhhRdYtmwZPXv2zLVvYGAgDzzwQK7HqlevzuDBg4mPj+fmm28+77ELYsGCBfTp04cbb7yRTz75pFBCUkpKSp4AJyIipZuG3omISKFr27YtALt27SrQ/iEhIQC4uV3Zv98tWrSIBx54gA4dOrBs2bI8IemDDz6gWbNmeHt7ExwczL333ptnyF+HDh1o0KABmzZtol27dvj4+PDkk0+yZ88eDMPgxRdf5M0336RmzZp4enrSokULvv/++zy1/Prrr/To0YPg4GC8vLxo3rw5y5Ytu6LXJyIiV496lEREpNDt2bMHgHLlyuXZ5nA4OHr0KACZmZkkJCTw7LPPEhUVRZs2bfLsn5SU5Nw/R4UKFfLst3jxYqKjo2nXrh2ffvop3t7eubZPnDiRMWPG0LNnTx5++GGOHDnC//3f/9GuXTs2b95MUFCQc99jx47RpUsX7r33Xh544AEqV67s3DZv3jxOnz7Nv//9bwzDYMqUKXTv3p3ff/8dd3d3AH7++WfatGlDeHg4o0aNwtfXl0WLFtGtWzcWL17MXXfdVbA3UkREbKOgJCIiVywnzKSlpbFx40aee+45PD09811k4ddff6VixYq5HqtXrx5ffPEFHh4eefbv2LFjnsdM08z18+bNm1m5ciXXX389y5cvzxOS/vjjD5599lkmTJjAk08+6Xy8e/fuNGnShJkzZ+Z6/ODBg8yaNYt///vfzsdywt/evXtJTEx0hsA6derQtWtXVq5c6Xy9Q4cOpWrVqnz//fd4enoCMHDgQK6//npGjhypoCQiUgIoKImIyBU7N8xERkbywQcfUKVKlTz7RkZG8tZbbwGQlZXFjh07mDJlCl26dOHrr7/OE6Jee+01ateufcHzHz9+nKysLKpUqZInJAF8/PHHZGdn07Nnz1y9UyEhIdSqVYs1a9bkCkqenp489NBD+Z6rV69euXrKcoYZ/v77785avvrqK8aNG8fp06c5ffq0c99OnTrx7LPPsm/fPsLDwy/4mkRExF4KSiIicsVywkxSUhLvvPMO69atc/aknMvX1zdXsOrcuTPXX389zZs35/nnn2fatGm59m/ZsuVFF3O46aabqFq1Kq+//jrBwcHMmDEj1/bExERM06RWrVr5Pj9nyFyO8PDwfHu3AKpWrZrr55zQdOLECQB27tyJaZqMGTOGMWPG5HuMw4cPKyiJiBRzCkoiInLFzg4z3bp14/rrr+f+++9nx44d+Pn5XfT5zZo1IzAwkHXr1l12Da+++ionTpzglVdeoVy5cowdO9a5LTs7G8Mw+N///oerq2ue555bY369Ujnyez78MxwwOzsbgBEjRtCpU6d8942KirrgaxEREfspKImISKFydXVl8uTJ3HDDDbz66qu5Lg57IQ6Hg+Tk5Ms+r4uLC++99x5JSUk899xzBAcHM2TIEABq1qyJaZpUr179osP4rlSNGjUAq5cqv/lVIiJSMmh5cBERKXQdOnSgZcuWTJ8+nbS0tIvuv2bNGpKTk2ncuPEVndfd3Z2PPvqINm3aMGzYMN5//33AWrTB1dWV5557Ls9CEKZpcuzYsSs679kqVapEhw4deOONNzhw4ECe7UeOHCm0c4mISNFRj5KIiBSJxx9/nHvuuYe5c+fyn//8x/l4UlISH3zwAfDPYg6vv/463t7eBe59uhAfHx8+++wz2rdvT9++fQkMDOTOO+9kwoQJjB49mj179tCtWzf8/f3ZvXs3S5YsYcCAAYwYMeKKz53jtdde4/rrr6dhw4b079+fGjVqcOjQITZs2MBff/3F1q1bC+1cIiJSNBSURESkSHTv3p2aNWvy4osv0r9/f+fcnr/++osHH3wQAMMwKFeuHO3bt+fZZ5/lX//6V6GcOzAw0LlceK9evfjf//7HqFGjqF27Ni+//DLPPfccABEREdxyyy3ceeedhXLeHPXr1+eHH37gueeeY+7cuRw7doxKlSrRpEkTnnnmmUI9l4iIFA3DPHcMgoiIiIiISBmnOUoiIiIiIiLnUFASERERERE5h4KSiIiIiIjIORSUREREREREzqGgJCIiIiIicg4FJRERERERkXMoKImIlDKGYTB27NgiPUeHDh3o0KFDoR3v3JrHjh2LYRgcPXq0SI4/d+5cDMNgz549zsciIyO5/fbbC+V8hWHPnj0YhsHcuXPtLkVEpExSUBKRYivni+KLL75odykAxMTEYBgGAQEBpKam5tmemJiIYRjFqma7xcXFYRgGH330kd2lyN8iIyOdn9Nzb7Vq1SrQMTIyMpg0aRJ169bFy8uLypUrc9ttt/HXX38598n5eznfbd++fc59O3TokO8+nTt3LvTXfzHZ2dlMmTKF6tWr4+XlRaNGjZg/f36efebOncudd95JREQEvr6+NGjQgAkTJpCWlnbVaxaRouFmdwEiIiWJm5sbKSkpfPrpp/Ts2TPXttjYWLy8vGz/opSamoqbW8n6z/vVrvnBBx/k3nvvxdPT86qds7iYPn06ycnJuR77448/ePrpp7nlllsu+vzMzExuu+02vvnmG/r370+jRo04ceIEGzduJCkpiSpVqgDw73//m44dO+Z6rmma/Oc//yEyMpLw8PBc26pUqcLkyZNzPRYWFnY5L/GKPPXUUzz//PP079+fFi1a8Mknn3D//fdjGAb33nsvACkpKTz00ENcd911/Oc//6FSpUps2LCBZ599ltWrV/PVV19hGMZVr11EClfJ+j+piIjNPD09adOmDfPnz88TlObNm8dtt93G4sWLC+18aWlpeHh44OJS8AEAXl5ehXb+q+Vq1+zq6oqrq+tVPWdx0a1btzyPTZgwAYDo6OiLPv/ll19m7dq1rF+/npYtW553v1atWtGqVatcj61fv56UlJR8zxMYGMgDDzxw0fNfrrFjxzJ37txcwy3PtW/fPqZNm8agQYN49dVXAXj44Ydp3749jz/+OPfccw+urq54eHgQHx9P69atnc/t378/kZGRzrB0bkgUkZJHQ+9EpMQ7fPgw/fr1o3Llynh5edG4cWPefffdPPsdO3aMBx98kICAAIKCgujTpw9bt2695Hkg999/P//73/84efKk87Hvv/+exMRE7r///jz7Hz9+nBEjRtCwYUP8/PwICAigS5cubN26Ndd+OcPUFixYwNNPP014eDg+Pj6cOnUKgA8//JD69evj5eVFgwYNWLJkCTExMURGRuY6zvnm++zcuZOYmBiCgoIIDAzkoYceIiUlJddz58yZw4033kilSpXw9PSkfv36vP766wV+by5XQeZV/fHHH0RFRdGgQQMOHToEwMmTJxk2bBgRERF4enoSFRXFCy+
     },
     "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": 23,
   "outputs": [],
   "source": [
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "from sklearn.gaussian_process.kernels import RBF\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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",
    "gp = GaussianProcessRegressor(kernel=RBF(length_scale=1.0), alpha=1e-8)\n",
    "gp.fit(X_1D, y_1D)\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": 25,
   "outputs": [],
   "source": [
    "# Predict the mean and variance of the test data using the scikit-learn implementation\n",
    "y_star, var = gp.predict(X_1D_star, return_std=True)\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": 26,
   "outputs": [
    {
     "data": {
      "text/plain": "<matplotlib.legend.Legend at 0x7fdfb0cd7100>"
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 1000x600 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAAINCAYAAADMabVmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACdTUlEQVR4nOzdd3gU5d6H8Xt20whpBAgJELr0EnpHBZQmgo2iiBUVu1jxqIh67L1hL4iKFZBipNmokRJ6EQgQICFAIJW03Xn/iOaVIyglybO7+X6uay82m9nNnVMYfjuzz1i2bduIiIiIiIhIqXCYDhAREREREfElGrJERERERERKkYYsERERERGRUqQhS0REREREpBRpyBIRERERESlFGrJERERERERKkYYsERERERGRUqQhS0REREREpBT5mQ7wdG63m3379hEaGoplWaZzRERERETEENu2ycrKombNmjgcJz5epSHrX+zbt4/Y2FjTGSIiIiIi4iGSk5OpXbv2Cb+vIetfhIaGAsX/QYaFhRmuERERERERUzIzM4mNjS2ZEU5EQ9a/+PMUwbCwMA1ZIiIiIiLyrx8j0sIXIiIiIiIipUhDloiIiIiISCnSkCUiIiIiIlKKNGSJiIiIiIiUIg1ZIiIiIiIipUhDloiIiIiISCnSkCUiIiIiIlKKNGSJiIiIiIiUIg1ZIiIiIiIipUhDloiIiIiISCnSkCUiIiIiIlKKNGSJiIiIiIiUIg1ZIiIiIiIipcjPdICcHJfbJiEpnbSsPKJCg+hUPxKnwzKdJSIiUiFoPywip0JDlheIX5/CxJkbScnIK3ksJjyICYOb079ljMEyERER36f9sIicKp0u6OHi16cwdsqqY/5iB0jNyGPslFXEr08xVCYiIuL7tB8WkdOhIcuDudw2E2duxC55xCaQgj/uFZs4cyMut32cZ4uIiMiZ+N/9cKjjAPUC1wDaD4vIP9OQ5cESktL/8s6ZzUN+U3im8pNUdhz54xFIycgjISndVKKIiIjP+ut+ONjKIKbB8xxq8Dm9qr8GuLUfFpET0pDlwdKy/v/UhFocpFPQr7xUK4ez6j5DqOPQcbcTERGR0vHX/Wv76LdI8S9e6GJ1tb30jHkWB0V/205EBDRkebSo0KCS+3upzj1F11OExfYgm3r1niPcuf9v24mIiEjp+HP/2jz4F9aEFx+tissIxbJtEiOO0KX2EwSQp/2wiPyNhiwP1ql+JDHhQfy5QOzWvE4E7R5BhMvNzkCoXfclmkQcpFP9SKOdIiIivqhT/Ujqhjl4yj2Di7NziMsI4dd9/6FpSif8bJtW7gN8WulFOtXUYs0iciwNWR7M6bCYMLg5QMmgtT2vHc5do6ha5GZ3IPhXf55DBzeYixQREfFRTofFh/UX0IYUbj7gJnHf7QAkZFxCg12DuOXQUTra63FOHgzZaYZrRcSTaMjycP1bxjBpVDuiw///VISd+a0JP3gDUS6bXX5wzayR7E9NNBcpIiLig47s/IX6Wz8A4Hn/m8girOR7+wLOZ33fKRBcjcKUNUz8rC/JyUtNpYqIh7Fs29a6o/8gMzOT8PBwMjIyCAsL+/cnlJHjXWk+Zd9yrp87huqF+bx9NJDgq2ZClXrGGkVERHxFYX4OIz7rRvW8HB6L7EzVy6b8bT/sdFhwaDvPfT2UyUFQzWUzqcdTNG082HS+iJSRk50NNGT9C08Zsk4kNXU1lb+8htD0JAirBVfNhKoNTWeJiIh4tbdmjOKNI2uo4raZPmQ6kZGNTrjtgbQNjJ19OVscbkLcNq+2HUfHuGvLsVZEysvJzgY6XdDLRUe3JfSaeKjWBDL3MmXqBWzfPs90loiIiNfatn0ubx9OBOCB+pf844AFUD2qBR9ePJv2diDZDoubEl9k/uKnyqFURDyVVw1Zv/zyC4MHD6ZmzZpYlsX06dP/9Tk//fQT7dq1IzAwkEaNGvHRRx+VeWe5C42Gq2czvWZjngnx45pf7mLL1lmmq0RERLyOq6iAR365nyLL4hwrlAG9JpzU80LDa/P2iAX0doRRYFnc/funfD1vXBnXioin8qohKycnhzZt2vDGG2+c1PZJSUkMGjSIc889l8TERO68806uv/56fvjhhzIuNSCkOudc9iXN3U4OOyyuXfwAGzZ9Y7pKRETEq0z54RbWOYoIcds81O8tLMfJ/1MpMCicF0Ys4JKAGNyWxQt7fuDQwsdAn8wQqXC89jNZlmUxbdo0hg4desJt7r//fmbPns369etLHhsxYgRHjhwhPj7+pH6Op38m639lZiQz9tshrHUUEuK2mdTxQeJaXm46S0RExOPt3r2IixfcRL7D4tGa53PJeS+c1uvYbjeTZlxOxy0L6ZiXDx3HwIBn4RQGNhHxTPpMFrB06VL69u17zGP9+vVj6dITL7Gan59PZmbmMTdvEhYeyzuXzaGdHUC2w+LG355kReJHprNEREQ8m9tN9oKJxBQV0ZkgLu7z3Gm/lOVwcPNFU+l47uOABb+9y7avRlKQn1V6vSLi0Xx6yEpNTaVGjRrHPFajRg0yMzM5evTocZ/z1FNPER4eXnKLjY0tj9RSVTkkmknDfqALlch1WNy8+nn2b5xmOktERMRzrfyA5rsS+Dotg6fPf+eUThM8oc43wqXvsy2wElflrOPmqb3JyU4989cVEY/n00PW6Rg/fjwZGRklt+TkZNNJpyU4uBqvD59HT6sytx4+Qo1vboStc01niYiIeBz78G6YV7zARWCfCVSLaVt6L97yEtLPn0iR5WA5eVzzVX8OHdxaeq8vIh7Jp4es6Oho9u/ff8xj+/fvJywsjEqVKh33OYGBgYSFhR1z81aBQeG8NmIho2ueDa58mHo57o0zTGeJiIh4DNvtZtx3w/k4yMIV2wk63VDqP6NT+xv5oMtEIt02mxwuRn93CXv2LCv1nyMinsOnh6yuXbuyYMGCYx6bN28eXbt2NVRU/pwBwXDZR9DiYjIp4spf7yP+l4mms0RERDzCzJ8eYj7ZvBoZwd4+D5XZ4hQtml3C5N5vUssFu51w5dzr2fL77DL5WSJinlcNWdnZ2SQmJpKYmAgUL9GemJjI7t27geJT/UaPHl2y/U033cSOHTu477772Lx5M2+++SZffvkld911l4l8c5z+cMl7fNG4O2uDArh/x1d8t/BB01UiIiJGHTywiWd2fQfA2KodqFPv7DL9eXXr9mLyoM9p7HJw0GlxzaL72bL+izL9mSJihlcNWStWrKBt27a0bVt8rvS4ceNo27YtjzzyCAApKSklAxdA/fr1mT17NvPmzaNNmza88MILvPfee/Tr189Iv1EOJ9ddNqPk2h0P7f5OF0kUEZEK7b/xY8h0WDRzO7m6/1vl8jOjarTkw0tm0d4OoEl+PvWm3QabZpXLzxaR8uO118kqL952nax/43YV8fTXQ/g874+jfzV6cXn/k7u4s4iIiK+Y++sT3L3jC/xsm6ndnqZJ4wvK9efnHT1M0fSxhGz5HiwH9qCXsDpcXa4NInLqdJ0sOS6H04/xl83k6sqNAHhq/y98NOt6w1UiIiLl58jhJP67bSoA14a3KPcBCyCoUhVChk2BdqPBdvPS4gm8PeNKbLe73FtEpPRpyKqALIeDcRd/ww1hzQGYvH8xGQufMFwlIiJSPlbOvZdMCxq6LG4c+J65EKcfDH6V1Z2u5sOIMF4/kshTXw3G7Soy1yQipUJDVgVlORzcdtEX3FulA++lphH+y3Ow4HHQ2aMiIuLLts6lz6Z5fLlvP091nUBAYKjZHsui7cBXeCCqJ5Zt83nebh74vA+F+Tlmu0TkjGjIquBGX/ghDc59tPiLX59n+5w7dKqCiIj4prxMmHUnAGe1v4FmzS4x2/MXVwx4k2fqX4qfbfO9K52bp55LTnaq6SwROU0asgS63QYDn2dZUCDD0hby5FcX6FQFERHxOe9/N5oNeWkQ2QDO/Y/pnL8ZcPajvNHyZiq5bZZxlOu+6k96+jbTWSJyGjRkSbFOY0jpcCWFFkzNS+axLwdq0BIREZ+RsPo9Xj66nVE1o9l3/qMQEGw66bi6dbiZD7s
     },
     "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
}