{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lab exercise #4: Convolution likelihoods, convolution kernels and SM kernels \n", "\n", "CS-E4075 2021" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Learning objectives\n", "\n", "After completing this lab exercise, you should be able to:\n", "\n", "- Apply a GP model to a deconvolution problem.\n", "- Using GPflow implement a convolutional kernel within a sparse variational gaussian process model (SVGP)\n", "- Implement spectral mixture kernels in GPflow" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import time\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "import gpflow\n", "import tensorflow as tf\n", "import tensorflow_probability as tfp\n", "\n", "from gpflow import set_trainable\n", "from gpflow.utilities import to_default_float, positive\n", "\n", "\n", "gpflow.config.set_default_float(np.float64)\n", "gpflow.config.set_default_jitter(1e-4)\n", "gpflow.config.set_default_summary_fmt(\"notebook\")\n", "\n", "# for reproducibility of this notebook:\n", "np.random.seed(123)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Convolutional likelihood models." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import scipy.integrate as integrate\n", "##https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html##\n", "\n", "np.random.seed(111)\n", "\n", "n = 100\n", "magnSigma2 = .5\n", "ell = .5\n", "sigma2 = 0.001\n", "x = np.linspace(-2,2,n)\n", "\n", "f = lambda t: np.exp(-(2*t)**2)\n", "\n", "def g(t):\n", " bool1 = t > -0.5\n", " bool2 = t < 0.5 \n", " return bool1*bool2*1\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 1a: Use the scipy numerical integration function to compute the convolution.\n", "- plot g, f and f_conv.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'f_conv' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m# we now sample noisy data from the convolution function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf_conv\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msigma2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mf_conv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'+'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'f_conv' is not defined" ] } ], "source": [ "# we now sample noisy data from the convolution function\n", "\n", "y = f_conv + np.sqrt(sigma2)*np.random.normal(size=100);\n", "plt.plot(x,f_conv)\n", "plt.plot(x,y,'+');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 1b: Compute a discrete convolution.\n", "- Use the discrete grid (x) as above and turn the integral in to a sum. (writing this out on paper may help)\n", "- Output a Convolution matrix which acts as the linear operator on the vector f, $f_{conv} = Cf $. \n", "- Plot the convolved f together with the discrete version." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 1c: Now solve the deconvolution problem.\n", "- Follow the adjusted GP regression formula's from the lecture slides using an SE kernel and the C you have calculated.\n", "- Plot the true underlying function alongside with the mean and 95% confidence interval, as shown on slide 30." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Convolutional kernel.\n", "\n", "We follow the example https://gpflow.readthedocs.io/en/master/notebooks/advanced/convolutional.html" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# for reproducibility of this notebook:\n", "np.random.seed(123)\n", "tf.random.set_seed(42)\n", "\n", "def make_rectangle(arr, x0, y0, x1, y1):\n", " # Draw the edges\n", " arr[y0:y1, x0] = 1\n", " arr[y0:y1, x1] = 1\n", " arr[y0, x0:x1] = 1\n", " arr[y1, x0 : x1 + 1] = 1\n", "\n", "\n", "def make_random_rectangle(arr):\n", " # Randomnly initialise the corners\n", " x0 = np.random.randint(1, arr.shape[1] - 3)\n", " y0 = np.random.randint(1, arr.shape[0] - 3)\n", " x1 = np.random.randint(x0 + 2, arr.shape[1] - 1)\n", " y1 = np.random.randint(y0 + 2, arr.shape[0] - 1)\n", " make_rectangle(arr, x0, y0, x1, y1)\n", " return x0, y0, x1, y1\n", "\n", "\n", "def make_rectangles_dataset(num, w, h):\n", " # Generate the rectangle dataset\n", " d, Y = np.zeros((num, h, w)), np.zeros((num, 1))\n", " for i, img in enumerate(d):\n", " for j in range(1000): # Finite number of tries to ensure \n", " x0, y0, x1, y1 = make_random_rectangle(img)\n", " rw, rh = y1 - y0, x1 - x0\n", " if rw == rh:\n", " img[:, :] = 0\n", " continue\n", " Y[i, 0] = rw > rh\n", " break\n", " return (\n", " d.reshape(num, w * h).astype(gpflow.config.default_float()),\n", " Y.astype(gpflow.config.default_float()),\n", " )" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "MAXITER = 100\n", "NUM_TRAIN_DATA = 100\n", "NUM_TEST_DATA = 300\n", "\n", "H = W = 14 # width and height of the images.\n", "IMAGE_SHAPE = [H, W]\n", "\n", "X, Y = data = make_rectangles_dataset(NUM_TRAIN_DATA, *IMAGE_SHAPE)\n", "Xt, Yt = test_data = make_rectangles_dataset(NUM_TEST_DATA, *IMAGE_SHAPE)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAACQCAYAAADQgbjgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAKh0lEQVR4nO3dTYxddRnH8e8jhTYUSCxI5VUgqQsg0kWDJsRE0ghoSGCjAVzUpLFu2IibRk1EwwLZuSAmjalTiIAuRFgQgXTDwsRQEohgEJDwUgZbeUl5SeQtj4s5xelw56Vz79zn3P/9fpKTO/fM6Zwn/c3pr+fMPXciM5EkSeP1ueoBJEmaRhawJEkFLGBJkgpYwJIkFbCAJUkqYAFLklTAApYkqYAFvEBEbIqI+yPi/Yh4OSJuWmS7iIhfRcSb3XJHRMS459XSzLMtEXFzRByIiA8iYmaZbX8UEf+OiCMRsTci1o9pTB2HaT5GLeDPuhP4ENgMfA/4TURcMmC7XcD1wGXAV4BrgR+Oa0itmHm2ZRa4Ddi71EYRcTWwG9gOXABcBPxirYfTqkztMRq+E9b/RcRG4G3g0sx8rlt3N/BaZu5esO1fgZnM3NM93wn8IDO/NuaxtQjzbFdE3Aacm5nfX+Tz9wAvZeZPuufbgd9n5hfHN6WWM+3HqGfAx/oy8MnRb4TOU8Cg/41d0n1uue1Uxzyn16A8N0fE6UXzaLCpPkYt4GOdAhxZsO4IcOoKtj0CnDLpP5NojHlOr0F5wuDsVWeqj1EL+FjvAactWHca8O4Ktj0NeC+9pt8n5jm9BuUJg7NXnak+Ri3gYz0HrIuILfPWXQY8M2DbZ7rPLbed6pjn9BqU56HMfLNoHg021ceoBTxPZr4P/An4ZURsjIgrgOuAuwdsfhdwS0ScExFnAz8GZsY2rJZlnu2JiHURsQE4ATghIjZExLoBm94F7IyIiyPi88DPMM/emfpjNDNd5i3AJuDPwPvAK8BN3fqvM3e54+h2AdwBvNUtd9C9qtylP4t5trUAtwK5YLkVOJ+5S5Tnz9v2FuAQ8A7wO2B99fwuAzOd2mPU25AkSSrgJWhJkgpYwJIkFbCAJUkqYAFLklTAApYkqcCg++dWLCKuAX7N3D15v83M25fa/qRYnxvYOMwuNaT/8j4f5geLvnXb8WRqnvVGmSeYaR94jLZlqTxXfRtSRJzA3LuYfBM4CDwO3JiZ/1jsz5wWm/KrsX1V+9No/C33806+Nfib4TgzNc96o8wTzLQPPEbbslSew1yCvhx4ITNfzMwPgfuYewcTTS4zbYt5tsdMGzJMAZ8DvDrv+cFu3TEiYldEHIiIAx/xwRC70xgsm6l5ThSP0fZ4jDZkmAIedEr9mevZmbknM7dl5rYTWT/E7jQGy2ZqnhPFY7Q9HqMNGaaADwLnzXt+LjA73DgqZqZtMc/2mGlDhnkV9OPAloi4EHgNuAG4aSRTqYqZtmVseT48++RafNlFXX321rHur0dGlum4M2vdar4nV13AmflxRNwMPMzcy+H3ZuZE/27GaWembTHP9phpW4a6DzgzHwIeGtEs6gEzbYt5tsdM2+E7YUmSVMACliSpgAUsSVIBC1iSpAJDvQhLkpaz2luGvE2mzhTf5rWkUX9PegYsSVIBC1iSpAIWsCRJBSxgSZIKWMCSJBWwgCVJKjDxtyFN860K3iogSZPLM2BJkgpYwJIkFbCAJUkqYAFLklTAApYkqYAFLElSgYm/DWkpLdymM823WUlSyzwDliSpgAUsSVIBC1iSpAIWsCRJBSxgSZIKWMCSJBUY6jakiHgJeBf4BPg4M7eNYijVMdO2mGd7zLQdo7gP+MrMfGMEX0f9YaZtMc/2mGkDvAQtSVKBYQs4gUci4omI2DWKgVTOTNtinu0x00YMewn6isycjYgzgUcj4tnMfGz+Bt03yC6ADZw85O40Bktmap4Tx2O0PR6jjRjqDDgzZ7vHw8D9wOUDttmTmdsyc9uJrB9mdxqD5TI1z8niMdoej9F2rLqAI2JjRJx69GPgKuDpUQ2m8TPTtphne8y0LcNcgt4M3B8RR7/OPZn5l5FMNSL+JqHj1vtMdVzMsz1jydR/O8dj1QWcmS8Cl41wFhUz07aYZ3vMtC3ehiRJUgELWJKkAhawJEkFLGBJkgpYwJIkFRjFL2ModfXZW6tHkKSJ47+d9TwDliSpgAUsSVIBC1iSpAIWsCRJBSxgSZIKWMCSJBWwgCVJKmABS5JUwAKWJKmABSxJUgELWJKkAhawJEkFLGBJkgpYwJIkFbCAJUkqYAFLklTAApYkqYAFLElSAQtYkqQCFrAkSQXWLbdBROwFrgUOZ+al3bpNwB+AC4CXgO9m5ttrN6ZGaVIzfXj2yeoRylx99tZFPzepeWpxZjodVnIGPANcs2DdbmB/Zm4B9nfPNTlmMNOWzGCerZnBTJu3bAFn5mPAWwtWXwfs6z7eB1w/4rm0hsy0LebZHjOdDqv9GfDmzHwdoHs8c3QjqYiZtsU822OmjVn2Z8DDiohdwC6ADZy81rvTGjPP9phpW8xzcqz2DPhQRJwF0D0eXmzDzNyTmdsyc9uJrF/l7jQGK8rUPCeGx2h7PEYbs9oCfhDY0X28A3hgNOOokJm2xTzbY6aNWcltSPcC3wDOiIiDwM+B24E/RsRO4BXgO2s5pEarxUyXuk1nUqz2NqsW85x2Zjodli3gzLxxkU9tH/EsGhMzbYt5tsdMp4PvhCVJUgELWJKkAhawJEkFLGBJkgpYwJIkFVjzd8KSxmGaf1NS35mNNJhnwJIkFbCAJUkqYAFLklTAApYkqYAFLElSAQtYkqQC3oakidHCbzxqldlIx88zYEmSCljAkiQVsIAlSSpgAUuSVMACliSpgAUsSVKByMzx7SziP8DL3dMzgDfGtvOlTdMsX8rML4ziCy3IE/rz99iXOWCC8gSP0RWamEw9RlekLM+xFvAxO444kJnbSna+gLOMRl9m78sc0K9ZjlefZneW0ejL7H2ZA2pn8RK0JEkFLGBJkgpUFvCewn0v5Cyj0ZfZ+zIH9GuW49Wn2Z1lNPoye1/mgMJZyn4GLEnSNPMStCRJBUoKOCKuiYh/RsQLEbG7YoZ5s7wUEX+PiCcj4sAY97s3Ig5HxNPz1m2KiEcj4vnu8fPjmmcY5vnpvs109HOY5wj0Jc9uFjPtjL2AI+IE4E7gW8DFwI0RcfG451jgyszcOuaXos8A1yxYtxvYn5lbgP3d814zz2PMYKZrwTyH0MM8wUyBmjPgy4EXMvPFzPwQuA+4rmCOUpn5GPDWgtXXAfu6j/cB1491qNUxz46ZtsU829O3TCsK+Bzg1XnPD3brqiTwSEQ8ERG7CucA2JyZrwN0j2cWz7MS5rk0Mx2OeQ6vT3mCmX5q3bh2NE8MWFf5UuwrMnM2Is4EHo2IZ7v/JWllzLM9fcrUPIfXpzzBTD9VcQZ8EDhv3vNzgdmCOQDIzNnu8TBwP3OXa6ocioizALrHw4WzrJR5Ls1Mh2CeI9GbPMFM56so4MeBLRFxYUScBNwAPFgwBxGxMSJOPfoxcBXw9NJ/ak09COzoPt4BPFA4y0qZ59LMdJXMc2R6kSeY6Wdk5tgX4NvAc8C/gJ9WzNDNcRHwVLc8M85ZgHuB14GPmPsf6k7gdOZehfd897ip6u/GPM20OlPzbCtPM/3s4jthSZJUwHfCkiSpgAUsSVIBC1iSpAIWsCRJBSxgSZIKWMCSJBWwgCVJKmABS5JU4H/7opTQDMYLkwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(8, 3))\n", "for i in range(4):\n", " plt.subplot(1, 4, i + 1)\n", " plt.imshow(X[i, :].reshape(*IMAGE_SHAPE))\n", " plt.title(Y[i, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each image contains the outline of a randomly generated rectangle, and is labelled according\n", "to whether the rectangle has larger width or height. $y_i = 1$ if height is greater than width and $y_i = 0$ if width is greater than height." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 2a: Comment on why a standard kernel would struggle with the rectangle dataset?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now try to fit a squared exponential kernel to the rectangle dataset, using the model class SVGP. SVGP is the sparse version of VGP discussed in lecture 4." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# This is a standard initialisation of the SVGP model class we need a kernel, \n", "# likelihood and inducing point locations.\n", "rbf_m = gpflow.models.SVGP(\n", " gpflow.kernels.SquaredExponential(),\n", " gpflow.likelihoods.Bernoulli(),\n", " gpflow.inducing_variables.InducingPoints(X.copy()),\n", ")# We use our data points as inducing points but usually you would want to use a sparse set." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# We turn off training for inducing point locations\n", "set_trainable(rbf_m.inducing_variable, False)\n", "\n", "# Similar to how we optimised in lab_3 but now using Adam and gradients are called inside optimizer.\n", "def run_adam(model, iterations, data): \n", " \"\"\"\n", " Utility function running the Adam optimizer\n", "\n", " :param model: GPflow model\n", " :param interations: number of iterations\n", " \"\"\"\n", " # Create an Adam Optimizer action\n", " logf = []\n", " training_loss = model.training_loss_closure(data, compile=True)\n", " optimizer = tf.optimizers.Adam(lr=0.01)\n", "\n", " @tf.function # This creates a static graph so we can reuse it every iteration.\n", " def optimization_step():\n", " optimizer.minimize(training_loss, model.trainable_variables)\n", "\n", " for step in range(iterations):\n", " optimization_step()\n", " if step % 10 == 0:\n", " elbo = -training_loss().numpy()\n", " logf.append(elbo)\n", " return logf\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXhW9Z338fc3hEAIgQQSCAkkAQn7TmQQ1yJ1rEvt4gIdtdVpGa1V28e2Y2tnetm57NPpOvZxbIdRsE5d61LF1mp13KBubEFQQCBhC4FAIAkhe77PH/cBAwmEJblPkvvzuq5cOef87uWbc53cn/v8zvIzd0dERKS5uLALEBGRzkfhICIiLSgcRESkBYWDiIi0oHAQEZEW4sMuoD2kpaV5bm5u2GWIiHQpy5cv3+Pu6a21dYtwyM3NZdmyZWGXISLSpZjZlmO1qVtJRERaUDiIiEgLCgcREWlB4SAiIi0oHEREpAWFg4iItBBaOJjZrWa23szWmtlPg2W5ZlZtZquCn9+GVZ+ISCwL5ToHM/sUcAUwyd1rzWxQs+ZN7j4ljLpERDqzxiZnd2UNO/ZVs2N/NcX7axg2IJHLJmW2+3uFdRHczcBP3L0WwN13h1SHiEincbCugeL9NcEHfzU79kV+bw/mS8praGg6cgyez07O7FbhMAo418zuAWqAb7v7+0HbcDNbCVQAP3D3t1p7ATObD8wHyM7OjkLJIiKnzt3ZW1V3+AN/x6GffdUUl0f2Asqq6o54To84I6Nfb7JSEsnPSSUrNZHMlESygp/MlESSenXMx3iHhYOZvQJktNJ0V/C+qcBM4EzgSTMbAewEst19r5lNB/5oZuPdveLoF3H3BcACgPz8fA1nJyKhqmtoYmf5J909zUPg0O/ahqYjnpOU0OPwB/7koSlkpiQytFkADEruRXyPcA4Nd1g4uPucY7WZ2c3AMx4Zo/Q9M2sC0ty9FDjU1bTczDYR2cvQjZNEJFRNTU5xeTVFew5SuLeK7WUHD3f37NhXTemBWo4edXlQci8yUxIZm9mPOeMGH/62f+ibf7/EeMwsnD+oDWF1K/0RmA28bmajgARgj5mlA2Xu3hjsSeQBm0OqUURiTFOTs7Oihi17qijcW0XRnioK9xykaG8VW8sOUtfsm39CfNzhD/kLRqeTldKHzJTeZKVGlmX0702v+B4h/jWnJ6xwWAgsNLM1QB3wZXd3MzsP+JGZNQCNwE3uXhZSjSLSDTU1ObsqayjcU8WWvQeDAKiiaG9kvnnXT6/4OHIHJjEiLYkLxwwiNy2J3IFJDE9LYnC/Xp32W397CCUc3L0OuLaV5U8DT0e/IhHpTtyd3ZW1kQ/9PVUUBSFQtDfyU1N/5B5AzoA+5KYlcf6odHLTkhg+MInctCQy+vUmLq77BsDxdIvxHEQk9rg7pQdqKdoTfPvfW8WWvZFuoC17qzhY13j4sQk94hg2IJHhaUmcMzKNnMMB0Ich/RPpEaMBcDwKBxHp1JqanK1lB1lXUsG6kko+3nUg6BKqoqpZAMTHGdnBHsBZIwYyPK0POUEXUGaKAuBkKRxEpNPYV1XHupLKSBDsrGTdrko2lFRSXR8JATPIHtCH4WlJzBg+gOFpScFxgD5kpSSGdtpnd6RwEJGoq21oZNPuKtaVVLC+pJKPSipZX1LBroraw48ZkJTAmIxk5s4YxtiMfozOSGbU4GQSE7ruGUBdicJBRDqMu1NcXsO6nRXBHkEkBDaXVh2+DURCjzhGDurL2SPTGJORzJiMfowZkkx63+59NlBnp3AQkXZRWVPP+iAADu0RrCuppLKm4fBjslISGTskmU+PGxwJgYxkctOS6KnuoE5H4SAiJ6WhsYmivVV8tLMyCIDIXsH2fdWHH5PcK57RGclcMSWT0Rn9GJuRzKiMZPr17hli5XIyFA4ickxNTc6G3ZUsK9rHyq37WVdSwce7Dxy+UrhHnDEiLYmp2anMm5HNmIxkRmckk5WSqC6hLk7hICKHVdc1smrbfpYVlbFsyz5WbN13uFsorW8CY4f04yuzchk9OJkxQ5I5I70vvXvqAHF3pHAQiWG7K2tYXrSP94v2sXxLGWuLKw4fKM4b1JfLJg0hP2cA+bmpZA/oo72BGKJwEIkRTU3OxtIDLCvad3jPYGvZQSByD6HJQ1OYf94I8nNTmZadSkqfhJArljApHES6qZr6Rgq27WfZln0sD37Kq+sBGJiUwPScVK6bmcP03FQmZPYnIV5nDMknFA4i3cSeA7UsC7qHlm3Zx5od5dQ3RrqIzkhP4uLxGeTnppKfO4DcgeoikuNTOIh0Qe7OpkNdRMFeQeGeKiByUdmkof258ZzhnJkzgGk5qQxIUheRnByFg0gXUFPfyAc7yg/vGSzfso99ByNdRKl9ejI9ZwBzzxxGfm4qE7L6d+lBZqRzUDiIdFK7K2r40wc7efGDElZt209dY+TaghFpScwZO5gzcwcwPTeVEWlJ6iKSdqdwEOlEyqrqeHHNThYXFPNuYRnuRK4tODuX/JxUpuekMrBvr7DLlBigcBAJWUVNPS+v3cXigmKWbNxDY5MzIj2J2y/M47JJmYwc1DfsEiUGKRxEQnCwroFXP9rN4oJiXl9fSl1jE0NTE5l/3ggun5TJ2CHJ6iqSUCkcRKKkpr6RNzaUsrigmFc/2k11fSOD+/XiurNyuHxyJpOH9lcgSKehcBDpQPWNTSzduIfFBTt5eW0JlbUNDEhK4IvTs7h8UiZn5g6I2QHspXMLJRzM7AlgdDCbAux39ylB2/eAfwQagdvc/aUwahQ5VY1NznuFZSxeXcyLH+xk38F6knvHc/GEDC6fnMmsMwZqOEvp9EIJB3e/5tC0mf0CKA+mxwFzgfFAJvCKmY1y98ZWX0ikk3B3Vmzdz+KCYv70wU5KK2vpk9CDOWMHc/nkTM4blaZrD6RLCbVbySIdrFcDs4NFVwCPu3stUGhmG4EZwNshlShyTO7O2uIKFq8u5oWCnezYX01CfByzRw/i8smZzB4zSOMdS5cV9jGHc4Fd7v5xMJ8FvNOsfXuwrAUzmw/MB8jOzu7IGkWO8PGuShav3skLBcVs3lNFfJxxbl4ad1w0ik+PG0yyRjuTbqDDwsHMXgEyWmm6y92fC6bnAY81f1orj/fWXt/dFwALAPLz81t9jEh72bK3ihdWRy5OW1dSSZzBzBED+dp5I7h4fAapuneRdDMdFg7uPud47WYWD3wBmN5s8XZgWLP5oUBx+1cn0raS8hpeWF3M4oJiCraXA5Cfk8rdnx3PZyZmMCi5d8gVinScMLuV5gDr3H17s2XPA4+a2S+JHJDOA94LoziJXTv2V3Pf/27kqeXbqG90Jmb15/uXjOHSSZlkpSSGXZ5IVIQZDnM5sksJd19rZk8CHwINwC06U0miZcf+au5/bSNPLtuGYcw9M5sbzs5lRLpuXyGxJ7RwcPevHGP5PcA90a1GYtnO8mr+87WNPPH+NgCuzh/GLZ8aSab2EiSGhX22kkhoSspruP/1jTz+3jYc56r8YXz9gjMYmton7NJEQqdwkJizq6KG37y+iUff20pTk3NV/lBu+dRIhYJIMwoHiRm7K2q4PwiFxibnymlD+cbskQwboFAQOZrCQbq93ZU1/Pb1zTzy7hYampwvTsviG5/KI3ugQkHkWBQO0m2VVtbyX29s4vfvbqG+0fn81CxunT2SnIFJYZcm0ukpHKTb2XMgEgr/884W6hqa+NzULG6bnUdumkJB5EQpHKTb2HuglgVvbubht7dQ29DI56Zk8Y3ZI3WdgsgpUDhIl1dWVReEQhHV9Y1cMTmTWy/M4wyFgsgpUzhIl7Wvqo4Fb23md3+LhMLlkzK57cI8Rg5SKIicLoWDdDn7D9bx329t5qGlRRysb+SySZncNnskeYOTwy5NpNtQOEiXUX6wngeWbGbR0iKq6hq4ZOIQbr8wj1EKBZF2p3CQTq/8YD0PBqFQWdvAJRMzuP3CUYzOUCiIdBSFg3Ra5dX1LFxSyMKlhVTWNPCZCRncdmEeY4f0C7s0kW5P4SCdTkVNPYuWFPHAks1U1jTw9+MHc/uFoxiXqVAQiRaFg3Qqb2/ay7eeWEVJRQ2fHjeYb87JY3xm/7DLEok5CgfpFOobm/iPVzZw/+ubGD4wiWe/Poup2alhlyUSsxQOErptZQe57fGVrNy6n6vzh/LDy8eT1EubpkiY9B8ooXq+oJi7nvkAgP83byqXT84MuSIRAYWDhKSqtoEfPr+Wp5ZvZ1p2CvfOnapxFUQ6EYWDRN2aHeXc+thKivZWcevskdx+YR7xPeLCLktEmlE4SNQ0NTkLlxby739Zx8CkXjz61ZmcdcbAsMsSkVaEEg5m9gQwOphNAfa7+xQzywU+AtYHbe+4+03Rr1DaW2llLXf8oYA3N5Ry0bjB/PsXJ5GalBB2WSJyDKGEg7tfc2jazH4BlDdr3uTuU6JflXSUNzaUcseTq6isaeDfPjeBa/8uGzMLuywROY5Qu5Us8glxNTA7zDqkY9Q2NPKzv6zngSWFjBrcl0e+OlP3QxLpIsI+5nAusMvdP262bLiZrQQqgB+4+1utPdHM5gPzAbKzszu8UDk5m0sPcNvjK1mzo4LrZuZw16Vj6d2zR9hlicgJ6rBwMLNXgIxWmu5y9+eC6XnAY83adgLZ7r7XzKYDfzSz8e5ecfSLuPsCYAFAfn6+t2/1cqrcnaeWb+eHz68lIT6O/7puOn8/vrXNQEQ6sw4LB3efc7x2M4sHvgBMb/acWqA2mF5uZpuAUcCyjqpT2k9FTT13PbuGxQXFzBwxgF9dM4Uh/RPDLktETkGY3UpzgHXuvv3QAjNLB8rcvdHMRgB5wOawCpQTt2LrPm57bCU7y2v49kWjuPmCkfSI00Fnka4qzHCYy5FdSgDnAT8yswagEbjJ3cuiXpmcsMYm57dvbOKXf93AkP69efKfzmJ6jm6YJ9LVhRYO7v6VVpY9DTwd/WrkVJSU1/DNJ1byzuYyLps0hB9/YSL9evcMuywRaQdhn60kXdRfP9zFd54qoLa+iZ9eOYmrpg/VtQsi3YjCQU5KTX0jP/7zRzz89hbGZ/bj1/OmckZ637DLEpF2pnCQE7ZhVyW3PrqS9bsq+eo5w/nOxaPpFa9rF0S6I4WDtMndeeTdrfzbCx+S3Dueh244kwtGDwq7LBHpQAoHOa79B+v47lOrefnDXZw3Kp1fXDWZ9OReYZclIh1M4SDH9M7mvXzriVXsOVDLXZeM5R/PGU6crl0QiQkKB2mhobGJX7/6Mfe9tpGcgUk8c/PZTBzaP+yyRCSKFA5yhG1lB/nmE6tYvmUfV04fyt2fHU9SL20mIrFG//Vy2JsbSrnl0RW4w71zp3DFlKywSxKRkCgcBIAte6u45dEVZKUksuC6fLIH9gm7JBEJkcJBqKlv5OuPrCDOjP++Pp9hAxQMIrFO4SDcvfhD1hZX8OCXFQwiEhEXdgESrmdWbOex97Zy0/lncOHYwWGXIyKdhMIhhm3YVcldz65hxvABfPuiUWGXIyKdiMIhRlXVNnDz75eT1Cue++ZNJb6HNgUR+YQ+EWKQu3PnMx9QuKeKX8+bwqB+vcMuSUQ6GYVDDPr9O1tYXFDMHReNZtYZaWGXIyKdkMIhxqzevp9/e+EjPjU6nZvPPyPsckSkk1I4xJD9B+u4+fcrSE/uxS+vnqKb6InIMek6hxjR1OTc8WQBuytrePKfziI1KSHskkSkE9OeQ4z4rzc38+q63dx1yVimZqeGXY6IdHKhhYOZTTGzd8xslZktM7MZwXIzs1+b2UYzW21m08Kqsbt4Z/Nefv7yei6dOIQvz8oNuxwR6QLaDAczm2BmDwcf4O+b2e/MbFI7vPdPgbvdfQrwr8E8wGeAvOBnPvCbdnivmLW7soZbH1tJzoA+/OSLEzHTcQYRadtxw8HMrgCeBV4HbgS+CrwBPB20nQ4H+gXT/YHiYPoK4GGPeAdIMbMhp/leMamxybn9sVVU1tRz/7XTSO7dM+ySRKSLaOuA9I+AT7t7UbNlBWb2v8Bzwc+p+ibwkpn9nEhIzQqWZwHbmj1ue7BsZ/Mnm9l8InsWZGdnn0YZ3dev/rqBtzfv5WdXTmJMRr+2nyAiEmgrHHoeFQwAuHuRmbX5NdTMXgEyWmm6C7gQ+Ja7P21mVwMPAnOA1vo9vJUaFgALAPLz81u0x7rX1u3mvtc2ck3+MK7KHxZ2OSLSxbQVDvVmlu3uW5svNLMcoKGtF3f3OcdqM7OHgduD2T8ADwTT24Hmn2ZD+aTLSU7Ajv3VfOvJVYwd0o+7rxgfdjki0gW1dUD6h8ArZvYVM5sYHJy+AXiZyEHk01EMnB9MzwY+DqafB64PzlqaCZS7+87WXkBaqmto4uuPrKCx0bn/H6bRu2ePsEsSkS7ouHsO7v5HMysE7gBuJdLlsxa42t0LTvO9vwbca2bxQA3B8QPgz8AlwEbgIHDDab5PTPnxnz+iYNt+fnvtNIanJYVdjoh0UW1eIR2EwPXt/cbuvgSY3spyB25p7/eLBX9avZOH/lbEjWcP5+IJOsFLRE5dW6eyppnZD83sNjPra2a/MbM1ZvacmY2MVpHStk2lB/juUwVMy07hzs+MCbscEeni2jrm8CjQi8gFae8BhcCVwAt8cgBZQlZd18jXf7+ChPg47vvSNBLidVcUETk9bXUrDXb371vkstot7n7oKuZ1Zqaun07iX55bw4bdlTx0wwwyUxLDLkdEuoG2vmI2wuHjAHuOamvqkIrkpDz5/jaeWr6dW2fncf6o9LDLEZFuoq09hxFm9jyRs5QOTRPMD+/QyqRNHxZX8C/PreGckWncfmFe2OWISDfSVjg0v3/Sz49qO3peoqiipp6vP7KclD49+Y+5U+ihgXtEpB21dZ3DG8dqM7MniNyET6LM3fnnp1azbV81j8+fSVrfXmGXJCLdzOmc1nJWu1UhJ2Xh0iJeXFPCP188mjNzB4Rdjoh0QzrnsYtZvmUf//fPH/HpcYP52rkjwi5HRLqp43YrHWcUNgM0OECUlVXV8Y1HVzAkpTc/v2qyBu4RkQ7T1gHpXxynbV17FiLH19TkfPOJVeytquOZm2fRP1HZLCIdp60D0p+KViFyfPe9tpE3N5Ryz+cnMCGrf9jliEg319a9lb7bbPqqo9p+3FFFyZGWfLyHX72ygc9PzeJLMzTqnYh0vLYOSM9tNv29o9oubudapBUl5TXc/vhKRqb35Z7PT9BxBhGJirbCwY4x3dq8tLP6xiZufWwF1fWN/ObaafRJaPMO6yIi7aKtTxs/xnRr89LOfvbSet4v2se9c6cwclBy2OWISAxpKxwmm1kFkb2ExGCaYL53h1YW415eW8KCNzdz7cxsrpiSFXY5IhJj2jpbSQMQh2Dr3oPc8YcCJg3tz79cNi7sckQkBukK6U6mpr6Rmx9ZjgH/+aVp9IpXPotI9OkIZyfzoxc+ZG1xBQ9cn8+wAX3CLkdEYpT2HDqRZ1du59F3t3LT+WcwZ9zgsMsRkRgWSjiY2RQze8fMVpnZMjObESy/wMzKg+WrzOxfw6gvDBt2VfL9Z9YwY/gAvn3RqLDLEZEYF1a30k+Bu939RTO7JJi/IGh7y90vC6muUFTVNvD1R1aQ1KsH982bSnwP7dCJSLjCCgcH+gXT/YHikOroFB54q5CNuw/w6Ff/jkH9dIawiIQvrHD4JvCSmf2cSNfWrGZtZ5lZAZHA+La7r23tBcxsPjAfIDu7695vqLahkf95ZwufGp3OrJFpYZcjIgJ0YDiY2StARitNdwEXAt9y96fN7GrgQWAOsALIcfcDQXfTH4G81l7f3RcACwDy8/O77NXaf1q9kz0HarnxnOFhlyIicliHhYO7zzlWm5k9DNwezP4BeCB4zqErsHH3P5vZ/WaW5u57OqrOMLk7Dy4pJG9QX87RXoOIdCJhHfksBs4PpmcDHwOYWYYFtx0NzmCKA/aGUmEUvF+0j7XFFdxw9nDdbVVEOpWwjjl8DbjXzOKBGoJjB8CVwM1m1gBUA3Pdvct2GbVl4ZJCUvr05PNTde8kEelcQgkHd18CTG9l+X3AfdGvKPq2lR3k5Q9LuOn8M0hM0C0yRKRz0Qn1IXn47SLMjOvOygm7FBGRFhQOIaiqbeDx97dxycQhDOmfGHY5IiItKBxC8PSK7VTWNHDj2blhlyIi0iqFQ5Q1NTmLlhYxNTuFqdmpYZcjItIqhUOUvb5hN4V7qrjhbF30JiKdl8IhyhYuKSKjX28+M6G1i8dFRDoHhUMUbdhVyZKNe7h+Vg49dedVEenE9AkVRYuWFtK7Zxzzzuy6NwoUkdigcIiSsqo6nlmxg89PHUpqUkLY5YiIHJfCIUoee28rtQ1NOn1VRLoEhUMU1Dc28fDbRZybl0be4OSwyxERaZPCIQr+/MFOdlXUcqNOXxWRLkLhEAULlxYxIi2J80elh12KiMgJUTh0sBVb91GwbT83nJ1LXJzGbBCRrkHh0MEWLikkuXc8X5g2NOxSREROmMKhAxXvr+bFNSXMm5FNUq+wxlUSETl5CocO9PDbW3B3rteYDSLSxSgcOkh1XSOPvbeViydkMDS1T9jliIicFIVDB3lm5XbKq+t191UR6ZIUDh2gqclZuKSQiVn9yc/RmA0i0vUoHDrAWxv3sKm0ihvPycVMp6+KSNcTSjiY2WQze9vMPjCzxWbWr1nb98xso5mtN7O/D6O+07VoaSHpyb24dGJm2KWIiJySsPYcHgDudPeJwLPAdwDMbBwwFxgPXAzcb2Y9QqrxlGzcfYDX15dy3cwcEuK1YyYiXVNYn16jgTeD6b8CXwymrwAed/dady8ENgIzQqjvlD30t0IS4uP40t9pzAYR6brCCoc1wGeD6auAYcF0FrCt2eO2B8taMLP5ZrbMzJaVlpZ2WKEno/xgPU8v38HnpmSS1rdX2OWIiJyyDgsHM3vFzNa08nMFcCNwi5ktB5KBukNPa+WlvLXXd/cF7p7v7vnp6Z3jhnaPv7+V6vpGnb4qIl1eh93Twd3ntPGQiwDMbBRwabBsO5/sRQAMBYrbv7r219DYxO/+VsRZIwYydki/tp8gItKJhXW20qDgdxzwA+C3QdPzwFwz62Vmw4E84L0wajxZL63dRXF5DTeeo70GEen6wjrmMM/MNgDriOwZLAJw97XAk8CHwF+AW9y9MaQaT8qipYVkD+jD7DGDwi5FROS0hXKrUHe/F7j3GG33APdEt6LTU7BtP8u27ONfLxtHD43ZICLdgE7EbweLlhbSt1c8V+VrzAYR6R4UDqdpV0UNf/pgJ1fnDyO5d8+wyxERaRcKh9P0+3e20NDkfGVWbtiliIi0G4XDaaipb+SRd7cyZ+xgsgdqzAYR6T4UDqfhuVU7KKuq40Zd9CYi3YzC4RS5O4uWFjF2SD9mjhgQdjkiIu1K4XCK3t60l3UlldxwtsZsEJHuR+FwihYuLWRgUgKfnawxG0Sk+1E4nIKiPVW8um43/zAzh949u9RwEyIiJ0ThcAoe+lsR8XHGtTM1ZoOIdE8Kh5NUUVPPH5Zt4/JJmQxK7h12OSIiHULhcJKefH8bVXUas0FEujeFw0lobHJ+93YRM3IHMHFo/7DLERHpMAqHk/DKR7vYVlbNDWfnhl2KiEiHUjichIVLCslKSeTT4waHXYqISIdSOJygtcXlvFtYxldm5RLfQ6tNRLo3fcqdoEVLi+iT0IOrzxzW9oNFRLo4hcMJKK2s5flVxVw5fSj9EzVmg4h0fwqHE/DIu1uoa2zSmA0iEjMUDm2obWjk9+9sYfaYQYxI7xt2OSIiUaFwaMMLBTvZc6BOp6+KSEwJJRzMbLKZvW1mH5jZYjPrFyzPNbNqM1sV/Pw2jPoOcXcWLi0kb1BfzhmZFmYpIiJRFdaewwPAne4+EXgW+E6ztk3uPiX4uSmc8iLeKyxjbXEFN54zXGM2iEhMCSscRgNvBtN/Bb4YUh3HtWhpESl9evK5KVlhlyIiElVhhcMa4LPB9FVA84sHhpvZSjN7w8zOPdYLmNl8M1tmZstKS0vbvcBtZQd5+cMSvjQjm8QEjdkgIrGlw8LBzF4xszWt/FwB3AjcYmbLgWSgLnjaTiDb3acC/wd49NDxiKO5+wJ3z3f3/PT09Hav/3d/KyLOjOvOymn31xYR6eziO+qF3X1OGw+5CMDMRgGXBs+pBWqD6eVmtgkYBSzrqDpbc6C2gSfe38YlE4cwpH9iNN9aRKRTCOtspUHB7zjgB8Bvg/l0M+sRTI8A8oDN0a7v6eXbqaxt0OmrIhKzwjrmMM/MNgDrgGJgUbD8PGC1mRUATwE3uXtZNAtranIWLS1kanYKU7NTo/nWIiKdRod1Kx2Pu98L3NvK8qeBp6Nf0SdeW7+bor0HueOi0WGWISISKl0hfZRFS4vI6NebiydkhF2KiEhoFA7NrC+pZMnGPVw/K4eeGrNBRGKYPgGbWbS0kN4945h3ZnbYpYiIhErhECirquPZlTv4wrShpCYlhF2OiEioFA6Bx97bSm1DEzdozAYREYUDQF1DEw+/XcS5eWnkDU4OuxwRkdApHIAX1+xkV0UtN54zPOxSREQ6hZgPB3dn4ZJCRqQncX5e+9+jSUSkK4r5cFixdT8F28u5YVYucXEas0FEBBQOLFxaSL/e8Xxh2tCwSxER6TRiOhyK91fzlzUlzJuRTVKvUO4kIiLSKcV0OBysa+DcvDSN2SAicpSY/ro8clAyD90wI+wyREQ6nZjecxARkdYpHEREpAWFg4iItKBwEBGRFhQOIiLSgsJBRERaUDiIiEgLCgcREWnB3D3sGk6bmZUCW07jJdKAPe1UTlendXEkrY9PaF0cqTusjxx3b/V21N0iHE6XmS1z9/yw6+gMtC6OpPXxCa2LI3X39aFuJRERaUHhICIiLSgcIhaEXUAnonVxJK2PT2hdHKlbrw8dcxARkRa05yAiIi0oHEREpIWYDgczu9jM1pvZRjO7M+x6os3MhpnZa2b2kZmtNbPbg+UDzOyvZvZx8Ds17Fqjxcx6mNlKM3shmB9uZu8G6+IJM0sIu8ZoMbMUM3vKzNYF28hZsbptmNm3gqXhqJYAAASiSURBVP+RNWb2mJn17u7bRsyGg5n1AP4T+AwwDphnZuPCrSrqGoA73H0sMBO4JVgHdwKvunse8GowHytuBz5qNv/vwK+CdbEP+MdQqgrHvcBf3H0MMJnIeom5bcPMsoDbgHx3nwD0AObSzbeNmA0HYAaw0d03u3sd8DhwRcg1RZW773T3FcF0JZF//iwi6+F3wcN+B3wunAqjy8yGApcCDwTzBswGngoeEkvroh9wHvAggLvXuft+YnTbIDKkcqKZxQN9gJ10820jlsMhC9jWbH57sCwmmVkuMBV4Fxjs7jshEiDAoPAqi6r/AL4LNAXzA4H97t4QzMfSNjICKAUWBd1sD5hZEjG4bbj7DuDnwFYioVAOLKebbxuxHA7WyrKYPK/XzPoCTwPfdPeKsOsJg5ldBux29+XNF7fy0FjZRuKBacBv3H0qUEUMdCG1JjiucgUwHMgEkoh0Rx+tW20bsRwO24FhzeaHAsUh1RIaM+tJJBgecfdngsW7zGxI0D4E2B1WfVF0NvBZMysi0sU4m8ieRErQlQCxtY1sB7a7+7vB/FNEwiIWt405QKG7l7p7PfAMMItuvm3Ecji8D+QFZxwkEDnA9HzINUVV0Kf+IPCRu/+yWdPzwJeD6S8Dz0W7tmhz9++5+1B3zyWyLfyvu/8D8BpwZfCwmFgXAO5eAmwzs9HBoguBD4nBbYNId9JMM+sT/M8cWhfdetuI6SukzewSIt8OewAL3f2ekEuKKjM7B3gL+IBP+tm/T+S4w5NANpF/jKvcvSyUIkNgZhcA33b3y8xsBJE9iQHASuBad68Ns75oMbMpRA7OJwCbgRuIfKGMuW3DzO4GriFyht9K4KtEjjF0220jpsNBRERaF8vdSiIicgwKBxERaUHhICIiLSgcRESkBYWDiIi0oHAQOYqZ/S34nWtmX2rn1/5+a+8l0tnoVFaRY2h+vcNJPKeHuzcep/2Au/dtj/pEOpL2HESOYmYHgsmfAOea2argfv49zOxnZva+ma02s38KHn9BMC7Go0QuKMTM/mhmy4MxAOYHy35C5M6eq8zskebvZRE/C8YL+MDMrmn22q83G1fhkeAqXZEOFd/2Q0Ri1p0023MIPuTL3f1MM+sFLDWzl4PHzgAmuHthMH+ju5eZWSLwvpk97e53mtk33H1KK+/1BWAKkXET0oLnvBm0TQXGE7l3z1Ii94Fa0v5/rsgntOcgcuIuAq43s1VEbjEyEMgL2t5rFgwAt5lZAfAOkRs85nF85wCPuXuju+8C3gDObPba2929CVgF5LbLXyNyHNpzEDlxBtzq7i8dsTBybKLqqPk5wFnuftDMXgd6n8BrH0vz+/U0ov9biQLtOYgcWyWQ3Gz+JeDm4DbnmNmoYACco/UH9gXBMIbIEKyH1B96/lHeBK4JjmukExmF7b12+StEToG+gYgc22qgIegeeojImMq5wIrgoHAprQ8N+RfgJjNbDawn0rV0yAJgtZmtCG4JfsizwFlAAZFBY77r7iVBuIhEnU5lFRGRFtStJCIiLSgcRESkBYWDiIi0oHAQEZEWFA4iItKCwkFERFpQOIiISAv/HyD81sD1QUK3AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "maxiter = 100\n", "\n", "logf = run_adam(rbf_m, maxiter,(X,Y))\n", "plt.plot(np.arange(maxiter)[::10], logf)\n", "plt.xlabel(\"iteration\")\n", "_ = plt.ylabel(\"ELBO\")" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Train acc: 100.0%\n", "Test acc : 70.33333333333334%\n" ] } ], "source": [ "train_acc = np.mean((rbf_m.predict_y(X)[0] > 0.5).numpy().astype(\"float\") == Y)\n", "test_acc = np.mean((rbf_m.predict_y(Xt)[0] > 0.5).numpy().astype(\"float\") == Yt)\n", "print(f\"Train acc: {train_acc * 100}%\\nTest acc : {test_acc*100}%\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
name class transform prior trainable shape dtype value
SVGP.kernel.variance ParameterSoftplus True () float640.7859763281946253
SVGP.kernel.lengthscalesParameterSoftplus True () float641.9783777154062827
SVGP.inducing_variable.ZParameterIdentity False (100, 196) float64[[0., 0., 0....
SVGP.q_mu ParameterIdentity True (100, 1) float64[[-0.49359234...
SVGP.q_sqrt ParameterFillTriangular True (1, 100, 100)float64[[[8.44783537e-01, 0.00000000e+00, 0.00000000e+00...
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "gpflow.utilities.print_summary(rbf_m) # We can here see all the parameters that we optimise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 2b: Change the number of inducing points in this model?\n", "\n", "Currently we are not using a reduced number inducing points $M =N$ and SVGP is the same as VGP. Try picking less but keep them randomnly selected from the data and see the effect on the model. Try $M=75,50,25$ and state the train and test accuracy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Convolutional Kernel\n", "To construct a convolutional GP we use a patch response function $g: \\mathbb{R}^E \\to \\mathbb{R}$ mapping form patch size $E$ to a real value. Then we make the overall function of the image $f: \\mathbb{R}^D \\to \\mathbb{R}$ a sum of the patch functions. If $g(.)$ is given a GP prior then a GP prior is induced on $f(.)$.\n", "\n", "$$ g \\sim \\mathcal{GP}(0,k_g(\\bf{z},\\bf{z}')), \\quad f ({x}) = \\sum_{p} g(\\bf{x}^{[p]}) $$\n", "\n", "$$ f ({x}) \\sim \\mathcal{GP} \\big(0,\\sum_{p} \\sum_{p'} k_g(\\bf{x}^{[p]},\\bf{x}'^{[p']}) \\big)$$\n", "\n", "- The kernel now works on the patches.\n", "- $k_g$ can be any base kernel but we will use the squared exponential" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [], "source": [ "def get_patches(X): # A function that will give us vector patches\n", " \"\"\"\n", " Extracts patches from the images X.\n", " :param X: (N x input_dim)\n", " :return: Patches (N, num_patches, patch_shape)\n", " \"\"\"\n", "\n", " num_data = tf.shape(X)[0]\n", " castX = tf.transpose(tf.reshape(X, [num_data, -1, 1]), [0, 2, 1])\n", " patches = tf.image.extract_patches(\n", " tf.reshape(castX, [-1, 14, 14, 1], name=\"rX\"),\n", " [1, 3, 3, 1],\n", " [1, 1, 1, 1],\n", " [1, 1, 1, 1],\n", " \"VALID\",\n", " )\n", " shp = tf.shape(patches) # img x out_rows x out_cols\n", " reshaped_patches = tf.reshape(\n", " patches, [num_data, 1 * shp[1] * shp[2], shp[3]]\n", " )\n", " return to_default_float(reshaped_patches)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 2c: Compute a value of the kernel and compare it against the gpflow implementation.\n", "\n", "- Do not compute the whole kernel simply K[0,0] will do. Remember to normalise by the total number of sums.\n", "- You can check your solution against the evaluated GPflow kernel.\n", "- Comment on the computation of the convolutional kernel against a standard kernel." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "base = gpflow.kernels.SquaredExponential()\n", "Xp = get_patches(X)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we build the SVGP model with the convolution kernel." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "WARNING:tensorflow:From :3: AffineScalar.__init__ (from tensorflow_probability.python.bijectors.affine_scalar) is deprecated and will be removed after 2020-01-01.\n", "Instructions for updating:\n", "`AffineScalar` bijector is deprecated; please use `tfb.Shift(loc)(tfb.Scale(...))` instead.\n" ] } ], "source": [ "## We constrain our parameters to help with optimisation.\n", "f64 = lambda x: np.array(x, dtype=np.float64)\n", "positive_with_min = lambda: tfp.bijectors.AffineScalar(shift=f64(1e-4))(tfp.bijectors.Softplus())\n", "constrained = lambda: tfp.bijectors.AffineScalar(shift=f64(1e-4), scale=f64(100.0))(\n", " tfp.bijectors.Sigmoid()\n", ")\n", "max_abs_1 = lambda: tfp.bijectors.AffineScalar(shift=f64(-2.0), scale=f64(4.0))(\n", " tfp.bijectors.Sigmoid()\n", ")\n", "\n", "# Build the convolutional kernel given specific image data.\n", "conv_k = gpflow.kernels.Convolutional(gpflow.kernels.SquaredExponential(), IMAGE_SHAPE, patch_shape)\n", "conv_k.base_kernel.lengthscales = gpflow.Parameter(1.0, transform=positive_with_min())\n", "\n", "# Weight scale and variance are non-identifiable. We also need to prevent variance from shooting off crazily.\n", "conv_k.base_kernel.variance = gpflow.Parameter(1.0, transform=constrained())\n", "conv_k.weights = gpflow.Parameter(conv_k.weights.numpy(), transform=max_abs_1())\n", "conv_f = gpflow.inducing_variables.InducingPatches(\n", " np.unique(conv_k.get_patches(X).numpy().reshape(-1, 9), axis=0)\n", ") # We place the inducing variables in the patch space and remove duplicates.\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU5b3H8c8PQgIkgRAStoSdgILIFha3umFVXFDrrrXVKtVq1fbW/V57W2trq9ZWvVeL1t5iXaDudakbVevCErawKoGAZCMJkH1PnvvHHHRMBiKS5Ewy3/frNa+c85wzM9/AJL+c53nOOeacQ0REJFg3vwOIiEj4UXEQEZEWVBxERKQFFQcREWlBxUFERFqI8jtAW0hKSnIjRozwO4aISKeyYsWKYudccqhtXaI4jBgxgoyMDL9jiIh0Kma2fV/b1K0kIiItqDiIiEgLKg4iItKCioOIiLSg4iAiIi2oOIiISAu+FAczu9fMNplZppm9aGYJQdtuM7MsM/vUzE72I5+ISKTz6zyHt4HbnHMNZvZb4DbgFjMbD1wITACGAO+Y2VjnXKNPOUVEfFXX0ERJdR1l1fWUVNVTWv3lo6y6ganDEzgmLeR5bAfFl+LgnHsraHUJcK63PBd41jlXC2SbWRYwA/ikgyOKiLS56rpGdlXWsruyjt2VdeypqmNPZT17qgLrJVWB5ZKqekqq6iiprqeqbv9/G19z3OiuUxyauQJY6C2nECgWe+V4bSIiYcc5R2l1PcUVtRSW11JcUUdxeS3FFbXsqqijuKKW4so6dnnr1fWhf9F3M0joHU1C7x706x3NkISeHDq4j7feg769etC3d3Tga9CjT88oorq3z+hAuxUHM3sHGBRi0x3OuZe9fe4AGoCn9j4txP4hb1VnZvOAeQDDhg076LwiIsEqaxvIL62hoLSGnWU1FJTVULj3a3kthWW1FFXUUtfQ1OK5Ud2M/nHR9I+NISk+hlFJsfSPjSYxLjrwNTaGxNhAIUiMjaZPzx506xbq159/2q04OOdm72+7mX0POB040X15r9IcYGjQbqlA3j5efz4wHyA9PV33OhWRr62+sYmC0hpy9lSTV+I9SqvJLamhoLSa/NIaymsaWjyvT88oBvbpyYA+McwYmciA+BiSgx9xga99e/XALLx+2R8oX7qVzOwU4BbgWOdcVdCmV4Cnzez3BAak04BlPkQUkU7MOUdReS2f7676yiNnTzW5e6rJL62mqdmflElxMQxJ6MnIpFiOHJ3EoL49Gdy3JwP79GSQVxB6R4dDT3zH8Os7fRiIAd72qusS59zVzrn1ZrYI2ECgu+lazVQSkVCcc+yqrGNrUSVbiyrI3lXJ9uIqtu2qZPuuqq/075vBwPieDEvszcyRiaT260Vqv96k9OtFSkIvBvXtSc8e3X38bsKPX7OVxuxn293A3R0YR0TCWFOTI2dPNZsLy8kqrAg8iirYUlhBWVDXT3T3bgxN7MWI/oG//If3782w/r0ZltiblIRe+uV/gCLnGElEwl5xRS0b88v4tKA88NhZzuadFV85CkiOj2FMchxnTh7CqKQ4RiXHMiopjpR+vegeZoO6nZmKg4h0OOccO3ZXsza3lHV5pWzIK2NjfhmF5bVf7JMUF8O4QXFcNGMYYwfGkTYwjjHJ8fTt3cPH5JFDxUFE2t3OshpWfV7CmpwSMnNKWJdbRml1PRCY9pk2MJ6j05IYP7gP4wf3YdygePrHxficOrKpOIhIm6praGJ9Xikrtu9hxfY9rN5RQn5pDRAoBIcMjmfOxEFMTElgYkpfxg6KIyZK4wHhRsVBRA5KZW0DK7bvYWn2LpZn72FNTgm13olhQxN7MX1EIpOGJjB5aAIThvTRwHAnoeIgIgekpr6RjG17+DCrmCVbd7Eut5SGJkf3bsZhKX25dNZw0of3Y9rwfgzo09PvuPINqTiIyH41NTk25Jfx/mdFfJRVTMb2PdQ1NBHVzZg0NIEfHjuKmSP7M214P2Jj9Culq9D/pIi0UFJVxwebi3nv00I++KyY4orALKJDBsVz2azhHJWWxIwRiSoGXZj+Z0UEgO27Knl7w07e3rCTjO17aGxyJPTuwbfSkjl2bDLfGptMcrxmEEUKFQeRCOWcY1NBOW+szeeNdQVsLqwAAkcH1xw7mhMOHcCk1ASdWBahVBxEIohzjo355by2No/X1xaQXVxJN4MZIxO5eOZ4Zh86kKGJvf2OKWFAxUEkAuzYXcUra/J4eXUun+2soHs344hR/bnymJGcPGEQSTrhTJpRcRDpoipqG3gtM4/nVuSwfNseANKH9+Ousw5jzmGDdAay7JeKg0gX4pxjWfZu/r4ih9cy86mub2R0ciw3nTyOMycNUZeRfG0qDiJdQGlVPc+tzOGppdvZWlRJXEwUZ00ZwnnpQ5kyNKHT35VMOp6Kg0gntjanlAWfbOMfmXnU1DcxZVgC9503iTkTB0XUXcuk7enTI9LJNDY53t5QwJ8/zGb5tj30ju7O2VNSuXTWMCYM6et3POkiVBxEOonK2gYWLt/BXz7OZsfualL79eI/TzuU86cPpU9P3eNA2paKg0iYK62q5/8+3sZfPs6mpKqe9OH9uP3UQzlp/ECiunfzO550USoOImGqqLyWxz/cyt8+2U5lXSOzDx3ANceNYdrwfn5Hkwig4iASZvZU1vHoB1tY8PF2ahsaOf3wIVxz3GgOHdzH72gSQXwpDmZ2FzAXaAIKge875/IsMN/uj8AcoMprX+lHRpGOVlZTz+P/zuaJD7OprGvgrMkpXH9iGiOTYv2OJhHIryOHe51z/wVgZtcDdwJXA6cCad5jJvCI91Wky6praGLBJ9t4aHEWpdX1zJk4iBtnj2XswHi/o0kE86U4OOfKglZjAectzwUWOOccsMTMEsxssHMuv8NDirQz5xxvrCvgt//cxPZdVRyTlsQtpxzCYSmajir+823MwczuBi4DSoHjveYUYEfQbjleW4viYGbzgHkAw4YNa9esIm1tzY4S7np1Axnb9zBuYDx/vWIGx45N9juWyBfarTiY2TvAoBCb7nDOveycuwO4w8xuA64Dfg6EOsffhWjDOTcfmA+Qnp4ech+RcLOnso573/qUZ5Z9Tv/YGO45ZyLnpQ/VPRMk7LRbcXDOzf6auz4NvEagOOQAQ4O2pQJ5bRxNpMM1NTkWZezgt//cRFlNA1ccNZIbZ6cRr5PXJEz5NVspzTm32Vs9E9jkLb8CXGdmzxIYiC7VeIN0dp/tLOeW5zNZ9XkJM0Yk8suzJnDIIE1LlfDm15jDPWY2jsBU1u0EZioBvE5gGmsWgamsl/sTT+Tg1TU08ej7W3ho8Wbie/bg9+dP4uwpKbpCqnQKfs1W+s4+2h1wbQfHEWlzmTkl3PxcJpsKyjlj0hD++4zxurmOdCo6Q1qkDdU3NvHQu5t5+F9ZJMXF8Nhl6Zw0fqDfsUQOmIqDSBvJLq7kxoWrWbOjhO9MTeXOM8bTt5cGnKVzUnEQOUjOOZ5dvoNf/mMD0VHd+N9LpjJn4mC/Y4kcFBUHkYNQWl3PLc9l8s/1BRw1pj/3nzeZQX17+h1L5KCpOIh8Q+tyS/nRUyvJK6nmjjmH8oOjR9JNJ7NJF6HiIHKA9nYj/fyV9fSPjWbhD2cxbXii37FE2pSKg8gBqKlv5PYX1/LCylyOSUviDxdM1hRV6ZJUHES+poLSGuY9mcHa3FJuODGN609M0zWRpMtScRD5GlbvKGHeggwqaxuY/12duyBdn4qDSCteXp3LTc9lMiA+hgU/OFLXRZKIoOIgsg/OOR54+zMeXJzFjJGJPHLJVI0vSMRQcRAJoaGxidtfXMuijBzOm5bK3WdPJDqqm9+xRDqMioNIM1V1DVz39CoWbyrk+hPG8JOTxupKqhJxVBxEguyqqOWKv2awNqeEu88+jEtmDvc7kogvVBxEPHkl1Vz6+FJyS6p55NJpnDwh1F1uRSKDioMIsGN3FRc9toTSqnr+duVMpo/QGc8S2VQcJOJtLargkseXUlXXyFNXzeTw1AS/I4n4TsVBItrmneVc/PhSmpocz1w1i/FDdA6DCKg4SATbVFDGxY8tJaqb8ey8WaQNjPc7kkjYUHGQiJRVWMGljy8luns3npk3i5FJsX5HEgkrOqtHIs6O3VVc+vhSAJ66aqYKg0gIvhYHM/uZmTkzS/LWzcweNLMsM8s0s6l+5pOuJ7+0moseW0JNQyN/u3Imo5Pj/I4kEpZ8Kw5mNhQ4Cfg8qPlUIM17zAMe8SGadFFF5bVc8thSSqvqWXDFDF1AT2Q//DxyeAC4GXBBbXOBBS5gCZBgZrpTuxy0spp6vvvnpeSX1vCXy6druqpIK3wpDmZ2JpDrnFvTbFMKsCNoPcdrC/Ua88wsw8wyioqK2impdAW1DY3MW5BBVmEF8y+bRrpOcBNpVbvNVjKzd4BQ1x+4A7gd+Haop4VocyHacM7NB+YDpKenh9xHpKnJ8R+L1rBk627+cMFkjklL9juSSKfQbsXBOTc7VLuZTQRGAmu8K12mAivNbAaBI4WhQbunAnntlVG6vrtf38irmfncduohnDUl5EGoiITQ4d1Kzrm1zrkBzrkRzrkRBArCVOdcAfAKcJk3a2kWUOqcy+/ojNI1PP7vrfz5w2y+f+QI5n1rlN9xRDqVcDsJ7nVgDpAFVAGX+xtHOqvX1+bzq9c2ctrEwdx5+njdj0HkAPleHLyjh73LDrjWvzTSFazNKeWni1YzdVgC958/iW7dVBhEDpTOkJYuZWdZDVcuWE7/2Bj+9N10evbo7nckkU7J9yMHkbZSXdfIVQsyKK9p4PlrjiQ5PsbvSCKdloqDdAlNTY6f/X0Na3NLmf/ddA4drLOfRQ6GupWkS3hocRavrc3n1lMO4aTxA/2OI9LpqThIp/fuxp088M5nnDM1RVNWRdqIioN0atnFldy4cDWHpfTh12dP1JRVkTai4iCdVmVtA1c/uYKobsajl07TzCSRNqQBaemUnHPc8nwmmwvL+esVM0jt19vvSCJdio4cpFP684fZvJqZz00nH6KL6Ym0AxUH6XSWbt3Fb97YxCkTBnH1sRqAFmkPKg7SqRRX1PLjZ1YxLLE39553uAagRdqJioN0Go1NjhufXU1pdT3/c/FU4nv28DuSSJelAWnpNB5avJkPs4r57XcmMn6IzoAWaU86cpBO4cPNxfzx3c18Z2oq56cPbf0JInJQVBwk7O0sq+GGZ1eRNiCOu86aoHEGkQ6g4iBhraGxiR8/s4rq+kb+95Kp9I5WT6hIR9BPmoS1BxdnsSx7N78/fxJjBsT7HUckYujIQcLWx1nFPLQ4MM5wztRUv+OIRBQVBwlLxRW13LBwNaOSYvnl3Al+xxGJOOpWkrDT1OT46aI1lFbXs+CKGcTG6GMq0tF05CBh508fbOWDz4r4+RnjdUc3EZ/4UhzM7L/NLNfMVnuPOUHbbjOzLDP71MxO9iOf+GfF9t3c99annDZxMBfPGOZ3HJGI1erxupkdBtwMjAccsAG43zmXeZDv/YBz7r5m7zUeuBCYAAwB3jGzsc65xoN8L+kESqrquP6Z1QxJ6MlvvqMb94j4ab9HDmY2F3gReA+4ArgSeB943tvW1uYCzzrnap1z2UAWMKMd3kfCjHOOm5/LpLC8hocvmkofXTdJxFetHTn8EjjJObctqG2NmS0GXvYe39R1ZnYZkAH8h3NuD5ACLAnaJ8dra8HM5gHzAIYNU/dDZ/fXj7fx1oad/OdphzJpaILfcUQiXmtjDj2aFQYAvLb9/mlnZu+Y2boQj7nAI8BoYDKQD9y/92khXsqFen3n3HznXLpzLj05WTd76czW5Zby69c3ccIhA/jB0SP9jiMitH7kUG9mw5xznwc3mtlwoGF/T3TOzf46AczsMeBVbzUHCL6qWiqQ93VeRzqnitoGrnt6JYmx0dx33iSNM4iEidaOHH5OYFD4+2Y20cwOM7PLgbeAO7/pm5rZ4KDVs4F13vIrwIVmFmNmI4E0YNk3fR8Jb845bn0+k893V/HHCyeTGBvtdyQR8ez3yME595KZZQP/AfyYQLfPeuB859yag3jf35nZZAJdRtuAH3rvt97MFhGYEdUAXKuZSl3X35Zs9+4DPY6Zo/r7HUdEgphzIbv0O5X09HSXkZHhdww5AJk5JZz7yCccOaY/T3xvOt26qTtJpKOZ2QrnXHqoba1NZU0ys5+b2fVmFmdmj3iDyi+b2Zj2iStdXWlVPT96aiVJcdE8cP5kFQaRMNTamMPTQAxf9v1nA+cSGEB+vH2jSVfknONnz62hoLSGhy6eSj+NM4iEpdZmKw10zt1ugSkk251zv/PaN5nZte2cTbqgx/+dzdve+QzThvfzO46I7ENrRw6NAC4wMFHcbFtTuySSLuvjLcXc889NnDxhoM5nEAlzrR05jDKzVwjMUtq7jLeun2752nJLqrnu6VWM6N9b5zOIdAKtFYfg6yfd12xb83WRkGrqG7n6yRXUNTQx/7J04nXdJJGw19p5Du/va5uZLSRwET6RfXLOcfuLa1mbW8pjl6UzOjnO70gi8jUczP0cjmizFNJlLfhkOy+szOWGE9M4afxAv+OIyNekO8FJu/koq5i7Xt3AiYcM4IYT0/yOIyIHYL/dSmY2dV+baOWqrBLZNu8s5+q/rWBUciwPXKgT3UQ6m9YGpO/fz7ZNbRlEuo6i8lou/7/lxER154nvT9eNe0Q6odYGpI/vqCDSNVTXNXLlggyKK2pZOO8IUvv19juSiHwDrV1b6eag5fOabft1e4WSzqmpyfHTRavJzCnhDxdM0R3dRDqx1gakLwxavq3ZtlPaOIt0Ys45fvXaRt5YV8Adcw7llMMG+R1JRA5Ca8XB9rEcal0i2APvbOaJj7K5/KgRujSGSBfQWnFw+1gOtS4Rav4HW3jw3c2cn57Kf502XpfGEOkCWputNMnMyggcJfTylvHWe7ZrMukUnlq6nV+/vonTDh/Mb845XFNWRbqI1mYrde+oINL5vLQql/98aR0nHDKAB86fTHcVBpEuQ2dIyzfy94wd/HTRamaOTOR/L5lKdJQ+SiJdiX6i5YD95aNsbnouk6PGJPHE96fTs4cOMEW6mtbGHES+4Jzj4cVZ3P/2Z5w8YSAPXjSFmCgVBpGuyLcjBzP7sZl9ambrzex3Qe23mVmWt+1kv/LJVznn+M0bm7j/7c84Z2oK/3PxVBUGkS7MlyMHMzuewI2EDnfO1ZrZAK99PIET7yYAQ4B3zGysc67Rj5wSUFPfyC3PZ/Ly6jy+d8Rwfn7GBM1KEuni/OpWuga4xzlXC+CcK/Ta5wLPeu3ZZpYFzAA+8SemFJbVcNWTK1izo4SbTh7Hj44brfMYRCKAX91KY4FjzGypmb1vZtO99hRgR9B+OV5bC2Y2z8wyzCyjqKioneNGprU5pZz58Eds3lnOn747jWuPH6PCIBIh2u3IwczeAUJdYOcO7337AbOA6cAiMxtF6EtyhDwT2zk3H5gPkJ6errO129g/1uRx03Nr6B8bw3NXH8n4IX38jiQiHajdioNzbva+tpnZNcALzjkHLDOzJiCJwJHC0KBdU4G89sooLVXXNfLLV9fzzLIdTBvej0cvnUZyfIzfsUSkg/nVrfQScAKAmY0FooFi4BXgQjOLMbORQBqwzKeMEWdjfhlnPPwhzy7fwTXHjebZebNUGEQilF8D0k8AT5jZOqAO+J53FLHezBYBG4AG4FrNVGp/zjmeXLKdX722kb69evDkFTM5Oi3J71gi4iNfioNzrg64dB/b7gbu7thEkSu7uJLbX1jLJ1t3cdy4ZO47bxJJcTpaEIl0OkM6QtU3NvHYv7fyx3c2Ex3VjV+fPZELpw/V+QsiAqg4RKRVn+/h9hfXsTG/jFMmDOIXcycwsI+uwC4iX1JxiCC5JdX87p+beHl1HgPiY3j00qmccthgv2OJSBhScYgAlbUNPPr+FuZ/sBUHXHv8aK45bgxxMfrvF5HQ9NuhC6upb+SppZ/zyHtbKK6oZe7kIdx08jhS+/X2O5qIhDkVhy6opr6RZ5YFikJheS1Hju7P/MumMXVYP7+jiUgnoeLQhZTX1LNw+Q4e+/dWdpbVMmtUIg9eNIVZo/r7HU1EOhkVhy6goLSGv3yUzdNLP6e8toFZoxJ54ILJHDlaJ7KJyDej4tBJOedYvaOEJz/Zzitr8mhyjjkTBzPvW6M4PDXB73gi0smpOHQyNfWNvLI6jyeXbGdtbilxMVFcOms4Pzh6JEMTNdAsIm1DxaGT2JhfxsLlO3hxVS6l1fWMHRjHXWcdxtlTUjQlVUTanH6rhLHymnr+sSafhcs/Z01OKdHdu/HtCQO5dNZwZo5M1I13RKTdqDiEmcYmx783F/HCylzeXF9AbUMT4wbGc+fp4zl7Sgr9YqP9jigiEUDFIUxszC/jxVW5vLQql8LyWvr26sH56UP5zrRUJqX21VGCiHQoFQcfFZTW8PLqXF5clcumgnKiuhnHjRvAudNSOP6QAcREdfc7oohEKBWHDlZV18Cb6wt4YWUuH2YV4xxMGZbAL+dO4LSJg+mveymISBhQcegATU2OJVt38fzKXN5Yl09VXSNDE3vx4xPSOHtKCiOTYv2OKCLyFSoO7aiwrIa/r8hhUcYOtu+qIj4mirmTh3DO1FTSh/fTOIKIhC0VhzbmnGNZ9m4e/zCbxZsKaWxyzBqVyE9mj+WUwwbRs4fGEUQk/Kk4tJGmJse7mwp55L0sVn5eQv/YaK46ZhQXTB+qbiMR6XR8KQ5mthAY560mACXOucnettuAHwCNwPXOuTf9yHgg3lpfwL1vfsrmwgpS+/XirrkTOC99qI4SRKTT8qU4OOcu2LtsZvcDpd7yeOBCYAIwBHjHzMY65xr9yNma8pp6fvGPDTy3Ioe0AXH84YLJnH74YKK6d/M7mojIQfG1W8kCI7LnAyd4TXOBZ51ztUC2mWUBM4BPfIq4T8u37eYnC1eTV1LNdceP4YbZafRQURCRLsLvMYdjgJ3Ouc3eegqwJGh7jtcWNpxzPPD2Zzz8ryxS+vVi0Q+PIH1Eot+xRETaVLsVBzN7BxgUYtMdzrmXveWLgGeCnxZif7eP158HzAMYNmzYQSQ9ME98tI0HF2dxztQUfnHmBOJ79uiw9xYR6SjtVhycc7P3t93MooBzgGlBzTnA0KD1VCBvH68/H5gPkJ6eHrKAtLU1O0q4542NnDR+IPefN0nnKYhIl+VnJ/lsYJNzLieo7RXgQjOLMbORQBqwzJd0zZRW13PdMysZEN+Te889XIVBRLo0P8ccLuSrXUo459ab2SJgA9AAXBsOM5Wcc9z6fCb5JTUsuvoIEnrrstki0rX5Vhycc9/fR/vdwN0dm2b/nlyynTfWFXDbqYcwdVg/v+OIiLQ7zb1sxfq8Un716kaOH5fMVceM8juOiEiHUHFoxaPvb6VXdHfuP38y3bppnEFEIoOKw35U1zXy7sadzJk4mETdnlNEIoiKw34s3lRIVV0jZxw+2O8oIiIdSsVhP17NzCMpLoaZo/r7HUVEpEOpOOxDRW0DizcVMmfiILprrEFEIoyKwz68u3EntQ1NnH74EL+jiIh0OBWHffjHmnwG9elJ+nCd1yAikUfFIYTS6no++KyIORMHa/qqiEQkFYcQ3t6wk7rGJk6fpFlKIhKZVBxCeDUzj5SEXkwZmuB3FBERX6g4NLOnso4PNxdz+uGDdeVVEYlYKg7NvLm+gIYmp1lKIhLRVByaeTUzn+H9e3NYSh+/o4iI+EbFIUhjk2PZtt3MPnSgupREJKKpOATZsbuKuoYmxg2M9zuKiIivVByCbCmqAGD0gDifk4iI+EvFIUhWYaA4jFFxEJEIp+IQJKuwguT4GPr26uF3FBERX6k4BMkqqmBMso4aRERUHDzOObYUVjB6QKzfUUREfOdLcTCzyWa2xMxWm1mGmc3w2s3MHjSzLDPLNLOpHZWpqKKWspoGHTmIiODfkcPvgF845yYDd3rrAKcCad5jHvBIRwX6cjBa01hFRPwqDg7YewpyXyDPW54LLHABS4AEM+uQS6Nu0UwlEZEvRPn0vjcCb5rZfQQK1JFeewqwI2i/HK8tv/kLmNk8AkcXDBs27KADbSmqJC4mioF9Yg76tUREOrt2Kw5m9g4wKMSmO4ATgZ845543s/OBPwOzgVDXrHChXt85Nx+YD5Cenh5ynwORVVjB6ORYXTZDRIR2LA7Oudn72mZmC4AbvNW/A497yznA0KBdU/myy6ldZRVWcOSY/h3xViIiYc+vMYc84Fhv+QRgs7f8CnCZN2tpFlDqnGvRpdTWymvqKSir0XiDiIjHrzGHq4A/mlkUUIM3dgC8DswBsoAq4PKOCLOlqBKA0ZrGKiIC+FQcnHMfAtNCtDvg2o7Oo5lKIiJfpTOkCVw2o0d3Y3hib7+jiIiEBRUHAoPRI/rHEtVd/xwiIqDiAAS6lTTeICLypYgvDnUNTWzfXaXxBhGRIBFfHLbvqqSxyak4iIgEifjioLu/iYi0pOLgFYdRybqPg4jIXhFfHLYUVZCS0Ive0X6dDygiEn4ivjhkFVUwWl1KIiJfEdHFoanJsaWwUnd/ExFpJqKLQ15pNdX1jbpvtIhIMxFdHL6YqaQjBxGRr4jo4hAbE8VJ4weSNlD3jRYRCRbRU3Smj0hk+ohEv2OIiISdiD5yEBGR0FQcRESkBRUHERFpQcVBRERaUHEQEZEWVBxERKQFFQcREWlBxUFERFow55zfGQ6amRUB27/h05OA4jaM0xE6W2blbV/K2766ct7hzrnkUBu6RHE4GGaW4ZxL9zvHgehsmZW3fSlv+4rUvOpWEhGRFlQcRESkBRUHmO93gG+gs2VW3valvO0rIvNG/JiDiIi0pCMHERFpQcVBRERaiOjiYGanmNmnZpZlZrf6nQfAzJ4ws0IzWxfUlmhmb5vZZu9rP6/dzOxBL3+mmU31Ie9QM/uXmW00s/VmdkM4Zzaznma2zMzWeHl/4bWPNLOlXt6FZhbttcd461ne9hEdmTcod3czW2Vmr4Z7XjPbZmZrzWy1mWV4bWH5efAyJJjZc2a2yWcoPi4AAAWVSURBVPscHxHmecd5/7Z7H2VmdmObZ3bOReQD6A5sAUYB0cAaYHwY5PoWMBVYF9T2O+BWb/lW4Lfe8hzgDcCAWcBSH/IOBqZ6y/HAZ8D4cM3svW+ct9wDWOrlWARc6LU/ClzjLf8IeNRbvhBY6NPn4qfA08Cr3nrY5gW2AUnN2sLy8+Bl+CtwpbccDSSEc95m2bsDBcDwts7s2zfl9wM4AngzaP024Da/c3lZRjQrDp8Cg73lwcCn3vKfgItC7edj9peBkzpDZqA3sBKYSeCM0qjmnw3gTeAIbznK2886OGcq8C5wAvCq90MeznlDFYew/DwAfYDs5v9G4Zo3RP5vAx+1R+ZI7lZKAXYEred4beFooHMuH8D7OsBrD6vvwevCmELgr/Gwzex10awGCoG3CRxBljjnGkJk+iKvt70U6N+ReYE/ADcDTd56f8I7rwPeMrMVZjbPawvXz8MooAj4i9dt97iZxYZx3uYuBJ7xlts0cyQXBwvR1tnm9YbN92BmccDzwI3OubL97RqirUMzO+canXOTCfxFPgM4dD+ZfM1rZqcDhc65FcHNIXYNi7yeo5xzU4FTgWvN7Fv72dfvvFEEunEfcc5NASoJdMnsi995v+CNM50J/L21XUO0tZo5kotDDjA0aD0VyPMpS2t2mtlgAO9rodceFt+DmfUgUBiecs694DWHdWYA51wJ8B6BftgEM4sKkemLvN72vsDuDox5FHCmmW0DniXQtfSHMM6Lcy7P+1oIvEigAIfr5yEHyHHOLfXWnyNQLMI1b7BTgZXOuZ3eeptmjuTisBxI82Z9RBM4PHvF50z78grwPW/5ewT69fe2X+bNRpgFlO49rOwoZmbAn4GNzrnfB20Ky8xmlmxmCd5yL2A2sBH4F3DuPvLu/T7OBRY7r+O2IzjnbnPOpTrnRhD4jC52zl0SrnnNLNbM4vcuE+gTX0eYfh6ccwXADjMb5zWdCGwI17zNXMSXXUrQ1pn9GkgJhweBUfzPCPQ53+F3Hi/TM0A+UE+g4v+AQJ/xu8Bm72uit68B/+PlXwuk+5D3aAKHqJnAau8xJ1wzA4cDq7y864A7vfZRwDIgi8BheozX3tNbz/K2j/Lxs3EcX85WCsu8Xq413mP93p+rcP08eBkmAxneZ+IloF845/Vy9AZ2AX2D2to0sy6fISIiLURyt5KIiOyDioOIiLSg4iAiIi2oOIiISAsqDiIi0oKKg0gQM/vY+zrCzC5u49e+PdR7iYQjTWUVCcHMjgN+5pw7/QCe090517if7RXOubi2yCfS3nTkIBLEzCq8xXuAY7zr5f/Eu1jfvWa23Lsm/g+9/Y+zwP0sniZwghFm9pJ30bn1ey88Z2b3AL2813sq+L28M1fvNbN1FrgPwgVBr/2efXmvgae8M9JF2l1U67uIRKRbCTpy8H7JlzrnpptZDPCRmb3l7TsDOMw5l+2tX+Gc2+1dnmO5mT3vnLvVzK5zgQv+NXcOgbN0JwFJ3nM+8LZNASYQuBbORwSutfRh23+7Il+lIweRr+fbBK5Ps5rAJcn7A2netmVBhQHgejNbAywhcMGzNPbvaOAZF7ha7E7gfWB60GvnOOeaCFyaZESbfDcirdCRg8jXY8CPnXNvfqUxMDZR2Wx9NoEb7lSZ2XsErnfU2mvvS23QciP6mZUOoiMHkdDKCdz2dK83gWu8y5NjZmO9q4421xfY4xWGQwhcDnyv+r3Pb+YD4AJvXCOZwK1il7XJdyHyDemvEJHQMoEGr3vo/4A/EujSWekNChcBZ4V43j+Bq80sk8DtGJcEbZsPZJrZShe47PZeLxK41ecaAle4vdk5V+AVFxFfaCqriIi0oG4lERFpQcVBRERaUHEQEZEWVBxERKQFFQcREWlBxUFERFpQcRARkRb+HwrNCFCX09V0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "maxiter = 700\n", "\n", "# we now define a new model with the special convolutional kernel structure.\n", "conv_m = gpflow.models.SVGP(conv_k, gpflow.likelihoods.Bernoulli(), conv_f)\n", "\n", "set_trainable(conv_m.inducing_variable, False) # for not do not train inducing variables.\n", "set_trainable(conv_m.kernel.base_kernel.variance, True)\n", "set_trainable(conv_m.kernel.weights, True)\n", "\n", "\n", "logf = run_adam(conv_m, maxiter,(X,Y))\n", "plt.plot(np.arange(maxiter)[::10], logf)\n", "plt.xlabel(\"iteration\")\n", "_ = plt.ylabel(\"ELBO\")" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Train acc: 100.0%\n", "Test acc : 98.33333333333333%\n" ] } ], "source": [ "train_acc = np.mean((conv_m.predict_y(X)[0] > 0.5).numpy().astype(\"float\") == Y)\n", "test_acc = np.mean((conv_m.predict_y(Xt)[0] > 0.5).numpy().astype(\"float\") == Yt)\n", "print(f\"Train acc: {train_acc * 100}%\\nTest acc : {test_acc*100}%\")\n" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
name class transform prior trainable shape dtype value
SVGP.kernel.base_kernel.variance ParameterSigmoid + AffineScalar True () float6462.711377272119876
SVGP.kernel.base_kernel.lengthscalesParameterSoftplus + AffineScalar True () float640.6426671133394776
SVGP.kernel.weights ParameterSigmoid + AffineScalar True (144,) float64[-0.77944789, -0.37365567, -0.2530622...
SVGP.inducing_variable.Z ParameterIdentity False (45, 9) float64[[0., 0., 0....
SVGP.q_mu ParameterIdentity True (45, 1) float64[[0.03312427...
SVGP.q_sqrt ParameterFillTriangular True (1, 45, 45)float64[[[0.07696406, 0., 0....
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "gpflow.utilities.print_summary(conv_m) # can you spot the difference with the standard SVGP summary." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 2d: Now turn training on for the patches and plot them" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sprectal mixture kernel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The SM kernel implementation compatible with gpflow. There are several kernels already defined in gpflow.kernels, and they all share the base 'Kernel' class. The kernel computation uses tensorflow functions (not numpy) so that we can utilize autodiff for later optimization.\n", "\n", "__Task 3a__: Complete the below SM kernel function to implement 1D Spectral Mixture kernel within the gpflow framework\n", "\n", "$$k(x,x') = \\sum_{q=1}^Q a_q \\exp( -2 \\pi^2 \\sigma_q^2 \\tau^2) \\cos( 2 \\pi \\mu_q \\tau), \\quad \\tau = x-x'$$\n", "\n", "that corresponds to a spectral mixture density\n", "\n", "$$S(\\omega) = \\sum_{q=1}^Q a_q \\mathcal{N}(\\omega | \\mu_q, \\sigma_q^2)$$\n", "\n", "where $a_q$ is the amplitude, $\\mu_q$ is the mean frequency and $\\sigma_q^2$ is the frequency variance in the frequency $\\omega$ domain. For info on implementing gpflow kernels, see https://gpflow.readthedocs.io/en/master/notebooks/tailor/kernel_design.html\n", "\n", "\n", "Hint: Compute the component-wise kernels into a kernel matrix of size (N,N,Q), and sum them together to obtain (N,N) kernel matrix." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "\n", "class UnivariateSMkernel(gpflow.kernels.Kernel):\n", " \"\"\"\n", " The spectral mixture (SM) kernel function [Wilson et al 2013]\n", " \n", " Inherits the 'Kernel' class from gpflow.kernels\n", " o We need to override the constructor to include the parameters (a,mu,sigma)\n", " o We need to override the 'K' function to compute the exp*cos kernel\n", " o We need to override the 'K_diag' function\n", " \n", " Supports only 1D inputs. Wilson et al 2014 contains multivariate version of the kernel \n", " \"\"\"\n", " \n", " def __init__(self, amps=None, freqs=None, ells=None, Q=1):\n", " \"\"\"\n", " - 'amps' are the amplitudes of the mixture components \n", " - 'freqs' are the frequencies means of the mixture components \n", " - 'ells' are the frequency stdevs of the mixture components \n", " - 'Q' is the number of components\n", " \"\"\"\n", " \n", " # init as 1D kernel \n", " super().__init__(active_dims=[0]) #input_dim=1\n", "\n", " # mixture amplitudes\n", " if amps is None:\n", " amps = np.random.uniform(0,1,Q)\n", " self.amplitudes = gpflow.Parameter(amps, transform=positive())\n", " \n", " # frequency means\n", " if freqs is None:\n", " freqs = np.random.uniform(0,1,Q)\n", " self.frequencies = gpflow.Parameter(freqs, transform=positive())\n", " \n", " # frequency variances\n", " if ells is None:\n", " ells = np.random.uniform(0,1,Q)\n", " self.lengthscales = gpflow.Parameter(ells, transform=positive())\n", " \n", " def K_diag(self, X):\n", " return tf.fill(tf.stack([tf.shape(X)[0]]), tf.reduce_sum(self.amplitudes))\n", " \n", " def K(self, X1, X2=None):\n", " \"\"\" Returns a kernel matrix of size (N1,N2) where (i,j)'th element is k(xi,xj)\n", " \n", " Input: X1 of size (N1,D)\n", " X2 of size (N2,D)\n", " \"\"\"\n", " if X2 is None:\n", " X2 = X1\n", " \n", " ##############################\n", " # INSERT CODE BELOW #\n", " ##############################\n", " \n", " ## >>>> SOLUTION CODE BELOW, REMOVE BEFORE GIVING TO STUDENTS <<<<\n", " \n", " kernel = None\n", " \n", " return kernel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Task 3b__: Fit a trigonometric function using the code below with both Gaussian and SM kernel.\n", "\n", "* Why are the fits different outside the training data?\n", "\n", "* Which kernel fits the data better, and why?\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
name class transform prior trainable shape dtype value
GPR.kernel.variance ParameterSoftplus True () float646.82905
GPR.kernel.lengthscalesParameterSoftplus True () float640.117581
GPR.likelihood.varianceParameterSoftplus + Shift True () float640.00574529
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "ename": "ValueError", "evalue": "in user code:\n\n /anaconda3/lib/python3.7/site-packages/gpflow/optimizers/scipy.py:104 _tf_eval *\n loss, grads = _compute_loss_and_gradients(closure, variables)\n /anaconda3/lib/python3.7/site-packages/gpflow/optimizers/scipy.py:165 _compute_loss_and_gradients *\n loss = loss_closure()\n /anaconda3/lib/python3.7/site-packages/gpflow/models/training_mixins.py:63 training_loss *\n return self._training_loss()\n /anaconda3/lib/python3.7/site-packages/gpflow/models/model.py:57 _training_loss *\n return -(self.maximum_log_likelihood_objective(*args, **kwargs) + self.log_prior_density())\n /anaconda3/lib/python3.7/site-packages/gpflow/models/gpr.py:65 maximum_log_likelihood_objective *\n return self.log_marginal_likelihood()\n /anaconda3/lib/python3.7/site-packages/gpflow/models/gpr.py:78 log_marginal_likelihood *\n k_diag = tf.linalg.diag_part(K)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/util/dispatch.py:201 wrapper **\n return target(*args, **kwargs)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/array_ops.py:2614 matrix_diag_part\n input=input, k=k, padding_value=padding_value, align=align, name=name)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/gen_array_ops.py:5150 matrix_diag_part_v3\n align=align, name=name)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/framework/op_def_library.py:540 _apply_op_helper\n (input_name, err))\n\n ValueError: Tried to convert 'input' to a tensor and failed. Error: None values not supported.\n", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 33\u001b[0m \u001b[0mms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgpflow\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGPR\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkernel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mksm\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 34\u001b[0m \u001b[0mms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlikelihood\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariance\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massign\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m \u001b[0mopt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtraining_loss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariables\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrainable_variables\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 36\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[0;31m# print and plot\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/gpflow/optimizers/scipy.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(self, closure, variables, method, step_callback, compile, **scipy_kwargs)\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 88\u001b[0m return scipy.optimize.minimize(\n\u001b[0;32m---> 89\u001b[0;31m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minitial_params\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjac\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mscipy_kwargs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 90\u001b[0m )\n\u001b[1;32m 91\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m 608\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'l-bfgs-b'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 609\u001b[0m return _minimize_lbfgsb(fun, x0, args, jac, bounds,\n\u001b[0;32m--> 610\u001b[0;31m callback=callback, **options)\n\u001b[0m\u001b[1;32m 611\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'tnc'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 612\u001b[0m return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/lbfgsb.py\u001b[0m in \u001b[0;36m_minimize_lbfgsb\u001b[0;34m(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)\u001b[0m\n\u001b[1;32m 343\u001b[0m \u001b[0;31m# until the completion of the current minimization iteration.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 344\u001b[0m \u001b[0;31m# Overwrite f and g:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 346\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mtask_str\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstartswith\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mb'NEW_X'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 347\u001b[0m \u001b[0;31m# new iteration\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/lbfgsb.py\u001b[0m in \u001b[0;36mfunc_and_grad\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m 293\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 294\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 295\u001b[0;31m \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 296\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjac\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 297\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m 325\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 326\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 327\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 328\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 329\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, x, *args)\u001b[0m\n\u001b[1;32m 63\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 65\u001b[0;31m \u001b[0mfg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 66\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjac\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 67\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mfg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/gpflow/optimizers/scipy.py\u001b[0m in \u001b[0;36m_eval\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m 109\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 110\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_eval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m->\u001b[0m \u001b[0mTuple\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 111\u001b[0;31m \u001b[0mloss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_tf_eval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconvert_to_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 112\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mloss\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnumpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 113\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m 826\u001b[0m \u001b[0mtracing_count\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexperimental_get_tracing_count\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 827\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_name\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mtm\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 828\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 829\u001b[0m \u001b[0mcompiler\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"xla\"\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_experimental_compile\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m\"nonXla\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 830\u001b[0m \u001b[0mnew_tracing_count\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexperimental_get_tracing_count\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m_call\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m 869\u001b[0m \u001b[0;31m# This is the first call of __call__, so we have to initialize.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 870\u001b[0m \u001b[0minitializers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 871\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_initialize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0madd_initializers_to\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minitializers\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 872\u001b[0m \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 873\u001b[0m \u001b[0;31m# At this point we know that the initialization is complete (or less\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36m_initialize\u001b[0;34m(self, args, kwds, add_initializers_to)\u001b[0m\n\u001b[1;32m 724\u001b[0m self._concrete_stateful_fn = (\n\u001b[1;32m 725\u001b[0m self._stateful_fn._get_concrete_function_internal_garbage_collected( # pylint: disable=protected-access\n\u001b[0;32m--> 726\u001b[0;31m *args, **kwds))\n\u001b[0m\u001b[1;32m 727\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 728\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0minvalid_creator_scope\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0munused_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0munused_kwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m_get_concrete_function_internal_garbage_collected\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2967\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2968\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2969\u001b[0;31m \u001b[0mgraph_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_define_function\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2970\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mgraph_function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2971\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m_maybe_define_function\u001b[0;34m(self, args, kwargs)\u001b[0m\n\u001b[1;32m 3359\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3360\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_function_cache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmissed\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcall_context_key\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3361\u001b[0;31m \u001b[0mgraph_function\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_create_graph_function\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3362\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_function_cache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprimary\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcache_key\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgraph_function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3363\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/function.py\u001b[0m in \u001b[0;36m_create_graph_function\u001b[0;34m(self, args, kwargs, override_flat_arg_shapes)\u001b[0m\n\u001b[1;32m 3204\u001b[0m \u001b[0marg_names\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0marg_names\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3205\u001b[0m \u001b[0moverride_flat_arg_shapes\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0moverride_flat_arg_shapes\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3206\u001b[0;31m capture_by_value=self._capture_by_value),\n\u001b[0m\u001b[1;32m 3207\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_function_attributes\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3208\u001b[0m \u001b[0mfunction_spec\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunction_spec\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/framework/func_graph.py\u001b[0m in \u001b[0;36mfunc_graph_from_py_func\u001b[0;34m(name, python_func, args, kwargs, signature, func_graph, autograph, autograph_options, add_control_dependencies, arg_names, op_return_value, collections, capture_by_value, override_flat_arg_shapes)\u001b[0m\n\u001b[1;32m 988\u001b[0m \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moriginal_func\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf_decorator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpython_func\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 989\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 990\u001b[0;31m \u001b[0mfunc_outputs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpython_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mfunc_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mfunc_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 991\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 992\u001b[0m \u001b[0;31m# invariant: `func_outputs` contains only Tensors, CompositeTensors,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/eager/def_function.py\u001b[0m in \u001b[0;36mwrapped_fn\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m 632\u001b[0m \u001b[0mxla_context\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mExit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 633\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 634\u001b[0;31m \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mweak_wrapped_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__wrapped__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 635\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 636\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/tensorflow/python/framework/func_graph.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 975\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m# pylint:disable=broad-except\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 976\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"ag_error_metadata\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 977\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mag_error_metadata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_exception\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 978\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 979\u001b[0m \u001b[0;32mraise\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mValueError\u001b[0m: in user code:\n\n /anaconda3/lib/python3.7/site-packages/gpflow/optimizers/scipy.py:104 _tf_eval *\n loss, grads = _compute_loss_and_gradients(closure, variables)\n /anaconda3/lib/python3.7/site-packages/gpflow/optimizers/scipy.py:165 _compute_loss_and_gradients *\n loss = loss_closure()\n /anaconda3/lib/python3.7/site-packages/gpflow/models/training_mixins.py:63 training_loss *\n return self._training_loss()\n /anaconda3/lib/python3.7/site-packages/gpflow/models/model.py:57 _training_loss *\n return -(self.maximum_log_likelihood_objective(*args, **kwargs) + self.log_prior_density())\n /anaconda3/lib/python3.7/site-packages/gpflow/models/gpr.py:65 maximum_log_likelihood_objective *\n return self.log_marginal_likelihood()\n /anaconda3/lib/python3.7/site-packages/gpflow/models/gpr.py:78 log_marginal_likelihood *\n k_diag = tf.linalg.diag_part(K)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/util/dispatch.py:201 wrapper **\n return target(*args, **kwargs)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/array_ops.py:2614 matrix_diag_part\n input=input, k=k, padding_value=padding_value, align=align, name=name)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/ops/gen_array_ops.py:5150 matrix_diag_part_v3\n align=align, name=name)\n /anaconda3/lib/python3.7/site-packages/tensorflow/python/framework/op_def_library.py:540 _apply_op_helper\n (input_name, err))\n\n ValueError: Tried to convert 'input' to a tensor and failed. Error: None values not supported.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAFlCAYAAAAZA3XlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5zcVb3/8ffZ3nfTEzYhhZBCSCgJJQEEBZQi5VIUr6BcRPTq9eLPCsr1Kt4i12u9Xi+CYgMVBAtNQUBqCoYQ0gMhpJdNNtlNdjc7szPz/f3xyZfZTGY3W2bmO+X1fDzmMVtmZ052Mjvv7/l+zuc4z/MEAAAAIK4o6AEAAAAA2YaQDAAAACQgJAMAAAAJCMkAAABAAkIyAAAAkICQDAAAACQoCeJBhw8f7k2YMCGIhwYAAEABeeWVV3Z7njeivz8XSEieMGGCFi9eHMRDAwAAoIA45zYO5OcotwAAAAASEJIBAACABIRkAAAAIAEhGQAAAEhASAYAAAASEJIBAACABIRkAAAAIAEhGQAAAEhASAYAAAASEJIBAACABIRkAAAAIAEhGQAAAEhQEsSDdnVJmzcH8cjoD+f6dpvut0v2cffrZJeiovi1fykujl8DAIIVjUqxWPzav3he/DrxIh36sf95Tx93/7wnfbkNcLiiAU0KBxKSOzqkVauCeGSkU1/+eDkXv133j/3PE+/LOamszC7l5VJlZfzif628XCotTd2/AwAKRSQidXZKoZAUDksHDth7tP+1UMgmthL/Vif7e+//DU92257+9vdlMgYYjAMHJKl0QHk3kJDsnDR8eBCPjFzjeTZzEY3af/S2NvujHonEv+8H6fp6acgQqa5OqqmRKiqCHTsAZJNIxP6G7tsnNTdLLS0WgrsrKbFLcbFdqqvtmjCLXNXcPPCfDSQkA33lXPyPdm8iEam9XdqzxwK159ls85gx0ogRFqAp3QBQaA4csJCwbZu0d6+VRhQV2SRCVZVNKgBIjpCMvOAH6erq+Ne6uqQtW6S33rI3hcZGu9TXMysCIH9FoxaM33rLJg6cs7+Nw4bxtw/oD0Iy8lZpqdTQYB9Ho9KOHdKmTfZmceyx0siRzC4DyB/RqM0Yr1tnNcXV1fZ3DsDAEJJREIqL44G5s1NautQW/E2daiUZA1v3CgDB8zxp1y5p5Ur7+9bQQBkFkAqEZBScigq7hMPSsmXS+vXSjBnS0KFBjwwA+qez07pFbd9OOAZSjfkzFKyysvipyIUL7Y2mqyvYMQFAXzU3Sy++aAvyRo+mow+Qaswko+BVVVknjC1bpKYm6cQT46UZAJBtPE/auNHKKxoaCMdAujCTDMhWfA8bZrXLCxbYjpDs7AQg23ietHatBeQRIwjIQDoxkwx0U1VlZRjLl1vT/alTWdQHIDt4nrR6tbRhgzRqFO3cgHTj7R9IUFJitcobNkivvRbf3Q8AguLPIG/YYH+fCMhA+hGSgSScszeipiZrF0dQBhCkjRulN98kIAOZREgGejF8uO1Y9eqrBGUAwWhqitcgE5CBzCEkA0cwbJgF5ddesx2tACBT2tvtbNbQoewQCmQaIRnog+HDbTZn1Sq6XgDIjGjUAnJZmV0AZBYhGeijESOsNdwbbwQ9EgCFYP16af9+qbY26JEAhYmQDPSRcxaU33hD2ro16NEAyGctLfa3ZtiwoEcCFC5CMtAPRUVWerFsmb2JAUCqRaPWq722lj7tQJB4+QH9VFJib16vvCJ1dgY9GgD5ZuNG28yoqirokQCFjZAMDEBlpV0vWybFYsGOBUD+OHBAev1162YBIFiEZGCAGhqk5mZp3bqgRwIgX7zxhp2tKikJeiQACMnAIAwfbiF59+6gRwIg17W22qLghoagRwJAIiQDg1JUZG9oS5dSnwwMVjRq5QadnYXXj9zzpDVrrA6ZXfWA7MAJHWCQysvtjX3FCmn2bN7ggP7o6pJ27bLFavv2xcNxaal09NF2KS8PdoyZsGePlW+NGhX0SAD4CMlACjQ0SDt3Sps2SePHBz0aIPtFo9KWLbZILRqVamqsfMkXidhmGhs2SCefnN/9gj1PWrvWfgcAsgflFkCKDBsmrVxpO2QB6FlLi/TSS7bNe12dbdLjd4zxlZRYaK6ulhYtknbsCGasmbBnj/1OqquDHgmA7gjJQIqUlNib3LJlNjMG4FCxmC10nT/fypJGjjxyF4fycmuHtmSJhcl89MYbBGQgGxGSgRSqqbGZ5DffDHokQHbp7JQWL7ZAOGJE/zbKKC2V6uulV1+VQqH0jTEILS0W/im1ALIPIRlIsWHDbLZs796gRwJkh337pAUL7ABy5MiBbbVcUWG1u6tWpX58QVq//vBSEwDZgZAMpFhRkdVZLl1qK/eBQrZnjwXk4uLB9/8dMkTavt26QOSD9nZb8FtbG/RIACSTspDsnCt2zr3qnHs0VfcJ5KrKSgvIa9cGPRIgOE1Ntuiupkbau3e9brhhnk4/vUw33DBPGzasH9B91tXZbHI+bAe/ZYuVktA2EshOqZxJvlnS6hTeH5DThg613q9NTUGPBMi8HTusBnnIECuVuO22f9CyZVMViTyiZcse1VVXTdJZZ0mf+IT0xz9ay7e+qKyU2tpy/3UVDlt7u7q6oEcCoCcpCcnOubGSLpb041TcH5APnLOAsGxZbuzGF41md3mI59kipw0bpNWrre57z57C25ktF+zcad0ohg61mdJnnpFWrPiFpJ9Keo+koZJiOnBAevll6etfl/7+723Hub6oq7MFgLn83Pshv7g42HEA6FmqNhP5rqQvSOqxsso5d5OkmyRp5MijU/SwQHbzd+NbudI2RMi206rRqNV4btwY7+9cViaNGyc1NvavA0E67d1rwbi11Wq+y8ps5vGNN6zOdeZMugNki507pVdesYAci0lf/ar06KOSNF7SCknfkfS4Zs6cqG9/e75eeEG65x5bwPaRj1hgfte7en+MigoLmc3Nh25Akis8z/69/J8FstugZ5Kdc++V1OR53iu93c7zvLs8z5vjed6c+voRg31YIGf4u/Ft3hz0SA7V0iK9+KK0fLl9PmKEXaqrpbfekp5/3lrZBVn7GY3a7OKCBRaKR460UFRXZyFs5Eg7bT1/fv720M0lTU3xgNzebqUUjz5qofajH92tmTM/rtLSezVr1kR9/ev3asgQ6dJLpQcesOtQSLr1Vvt/eSRVVRY0c1FLi/1+KiqCHgmA3qSi3OIMSZc65zZI+o2kdznn7k3B/QJ5Y9gwacUKa4UVNM+zmeP58+1U78iRh75Zl5TYeIcOtS2DFy8OpjdtKCQ9/PB6XXHFPF12WZk+9al52rLl8FRUU2PB/uWX2e0wSDt3xmuQd+6UbrhBeu01adQo6ac/lT72seH66U9f1IIFId1zz3yNHTvp7Z8tK5P+5V+kD33IDoy++EU7S9CbmhqbSW5vT/M/LA02b7azTACy26BDsud5t3qeN9bzvAmSrpH0jOd51w56ZEAeKSmxN/Wg28LFYjYzu3KlBeHe+rP6AXr/fmnhwsyGkY4O64pw223XavXqBYpEurRs2SLdcsvntWKFLQrrXo9aUWEzi4sX28wyMmvr1vgMcnOz9LGPSZs2SVOmWEA+9tgj34dz0qc+JV18sR0gffGLR/4/V1wsbduWmn9DpoRCNmbavgHZjz7JQIZUV9sCvlWrgllwFIlYacVbb8W3A96wYb2uueZGnXrqTbr88v/QE09sO2zW2O9t628GkW7t7RaQo1Hp9dcXSxoj6duStmvNmod0/fXSe98rXXml9KtfxbsiVFfbx6+/nv4xwvi1tUuX2kFXS4v08Y/bTPIJJ0h3323/1/rKOSu3OOYYC9nf/37vt6+vt7MiudQObtcuux7IhioAMiulL1PP8571PO+9qbxPIJ8MHWqzbhs3ZvZxQyGb6duxw05/x2LS/fdLH/hAjdat+7Fisbu0ZcuX9OUvH6V3v1v61rcO3bChpsZOiS9cmN6g3NZmAXnHjvW6+eYzFIl8UtLrkv6fpJEqKdmtqVMtHG3aJH3723Zaf8cO+/khQ+x3my+bTWSzcNjKKdassVr2tWs36PLL39LWrVJl5Wp94Qtvqbo6fvtIxBZe7tplz09PJTwVFdK//7vNEj/0kD1GT0pK7MxMS0tq/23p4nnWnYVZZCA3cCwLZJBztvBs1Spp9+7MPKY/M7t/vz32nj3SP/2T9M1vSl1dIyWtk/RLST+XtFzt7dKvfy1dcYX0u9/FZ72rq62OMl1B2Q/IxcXSv//7x7R8+ZdknRBqJP1ekyd/RL/97T7dd5/0xBM2/lGj7Hd54422MYNzFqDzZbOJbOR5tkDvxRel115br1tvnad580bq+uv3KRSaKGm5Dhw4U3fc8UFJ9jw0N1s9/qhR0qxZVn7R1WUzztHo4Y8xebL04Q/bx9/4RvLb+MrLc6fkYv9+u7BgD8gNzgvgvO+UKXO8X/96ccYfF8gWoZC9Wc6bl95ZpeZm61dbVubvemb1ouvX26x2be2t2rjxDkn2d2DWrLn6whfm68474x0GzjtP+td/jdcvt7fbLOLpp6euhdX+/RaQS0stLL/3vaskHSepWdINKi39sxYsOHzqcd8+6eabrYzk6KOln//cfp9NTRbGGhtTM75sFYvZc7pjhx38+LOzxcUWHisq7Lqy0i6lpXYpK7Ov9+eUfzRqj/HGGzZzW18vfeIT87Rs2XJJf5F0uqS1kt4hqUmlpWV69tmQWlqkSZPsUlZ26Njfeis+E53YL7izU7rqKvu3ffWrVmLT07haWqxtXEmqmpqmydq1dgZk6NCgRwIUjuZm6T3vqVjueZ2z+vuzzCQDAfCDy9/+ZovUUs3zLIAsXGgzwC0t63X99efo/PMXa/16ady4sO67T/re9z6qWbNOV2lpmWbNmqvbb79X06ZJ3/2unfKurpaeesqCtV/CUF2d2tKL1la7r7Iyae3aTbriiq2ygLxS0hxJD2v69NlJf7auTvrBD2yB2KZN0le+Yv/2hgarTe5tBjLXtbZanfjLL0uvvrpen/jEPF14YZk+85l52rvXuoD4O9OtX2+b2ixebD/z3HPSX/5i10uWWFDdssVuu2ePBe+9e600YtMmu83TT9vPRyI2I1xRIa1cuUrSo7KA/JakcyXZLhlTpszWvn3SKadI06YdGpAlC+jHHCMdf7w9TuJ8TUWFtZCTpP/7v57LM4qL7XlubU3RLzZNolHrasEOe0DuYCYZCND+/VYicMopqdu4IxSycoPt2628orhYuuGGeVq27EOSPi5pvaZPv1m//OUjR7yvDRtspnbrVpuV/d73pAkT7Hvt7fZYp5wSX9zXX83NdqDgt/O6+uod6uoaLWmRpAvkXKtmzjxdt99+7yEtwxJt3Spdd53NLP/Lv0iXXWaB78QTpTFjBja2bLZtm9XqVlfbxZ7fBW9/f9asubrnnvlHvJ+uLruEwxZ+k70dOBfvHtJ95vnAAemCC5aovf1kSVslnaWqql3q6gprypTZ+uIX79Vll03qUyhcvdrCeOLGILGYdO21dsDzhS9I73tf8p9vbbXgPmPGkR8rKM3NdkDTn4WMAAaPmWQgR9XWWhBYtMhm/QaruVl66SWrdx41Kn4Ke8WKo2QBuVPSlVq37sk+3d+ECbYb2nHHWRD9yEesfZxk4ayqymYm/YVzfeV5Fopeftlm1nbulD76UR0MyC9JerekFpWUlOqee+Zr+PBJ2r3bgu+uXXZpaYmHusZG6fOft4+/+1379/tbF+dbbfK2bdKrr9ope39h3KpVh046rF7d695ObystteewocECqr+hTPfL8OF2EJMYkG++WWpvP1klJbtUUvIezZo1Wr/61Wt6+umQvvvd+br88r4FZMlqkCsq7H67KyqyenPJSml6ap9YU2MHhdn8XG/dSi0ykGsIyUDA6uosDMyfP/DFfOGwzR4vWmSlHN1rHtvbpeLiHx787IuSlr5dvuCfjveD565dh4f1YcOkH/1IOvNMm7H7x3+02V/J3vSHDLHOGatXx9ux9SYUslnQFSvi3T4+9jF77OrqVyVdIMl2XZk2bbaamuxU9cyZNoZ3vEM69VQ7CGhqiperXHCB1Xjv32+n5ysq7N++d+/AfqfZqKXFfnfDhx9afztt2lxJ50v6mqS7VVt7v370IyunSHU5z+7d9n9gyRIL0fffP0ILF654+2Cmo8Oen/7U2peW2vObbLOdc86xeuadO6XHH0/+88XF9n8vGzbrSaary0I821ADuYWQDGSBmhq7LFpkYbOvG2JEoxYyn3/eakoTd8+TpLvusi4WlZUrVVJyl2bNmquvfvVe7dxpt50928LnGWdYb9uqKguf3cNyZaX03/9tQbSjQ/rnf5YefHCHbrhhns46q0y33jpPL7+8Xi+8YGEgWVgOhax84/nn4zPd69fH651PPVX6yU8aNGvWTJWWlmnGjLn6zGfufTscH3WUBa+qKgvXxx9vY/ZbizknffazFpgeecS21K6qssfMB11dNoNcWxsPyJGIdN990tatT0t6UtJXJN2oPXsu19132+/j/POlW26Rnn22bwcxvXn9des6sWKFlbHceac0fnx8fH4N8kDqbv1txhODblFRvNPF/ff33GO8uDjegzjb7N1rs9z0RgZyCzXJQBbxuxUUFdkp6FGjku+K19FhgWD9ejtFPXSozcYl2rbNNt2IRKRf/tIWUMViFoKnT5cmTrRwmai11WamW1rsvv1QFotZWH7gAUmKyko4fizJ6mD/7//mq7XVAsvQoRb8Pc/uxw+yDQ12f4sWWZ1pe7vNAP/Xf8UD/t69ttDr5JOPPPvW2Wn35XkWIO+4Q/rtb6Wzz7ax7tpls8/de/bmolWr7EBo2DD7fPt223hjxQr7fNIkO5hobLTfxdq1LXriic3q6Jj59n2MHGmt/S6//PD63954nh14fPOb9v9t1iz73fpnLPz/UyefPLga8P37ravKiBGH/r8Mh20nvr17pR//2GrNE4XDFtTf8Y6BP366vPKKHXQykwxk3mBqkrO8YQ5QWIqKLARFItYuavVqC8l1dRYsw2GbaQuH7ba1tb3P2t15pwWHCy+0gCxZaDz2WAtVPamvl047zVbjr1plj1NZaY/5+c/b9+++u1jS3ZJmS7pZq1e/ovJyC2KxmIWpffss7JSVWSjzg8/DD1v3jGhUOvdc6etfj3c/aG62+z/ppMM7IiRTUWGzly+9ZLPVN95o9//cczabPGSIBcrJk/vwBGSp1labEfcXfa1fL33yk/ZcjhplM8VnnnlosLzhhovU0bFA0jhJ71N5+afU1DRed95pQfPcc6Wrr7azB8kOlHwrVlgHkcUH5zUuvFC67TYr65EsQO/ebf+/BrtIsrZWGj3aDqq6/78uK7Ng/9Of2gFaspBcVmY/d+BA79utZ1ooZM9Tfw5KAGQHQjKQhUpK4m+qXV02C+V5FlKrqy1EHsmmTdKf/mT39fGP29fa2mz2ry+BsajITqXX11v9aShks8DOWYnEn/7079qy5bOy2eSTNGnSNw/52cRuHVu2rNdtt31EK1e+T573j5KsI8WnPhU/Dd3cbOM78cT+9bytqrKf8bsHXH65nZr/2c+sx+7GjXZQkIunuz3PDlSqq+13v3271QQ3N9vM7X//d/IDpfhivs2SvqVo9H/0wx+G9MADVvLyxBN2aWy0mfwZMyzklpXZjO3atXa71avtXurqpM99zkJy91C9Z4+VwvR20NUfkyZZfX7iv+nKK23x3jPPWChPFjqdswOzbArJe/bYdW8HIgCyUw6+ZQCFpbTUAlJNjYXBvobHX/zCAtbFF8dPwXd0WJlFf8JiQ4OFqKqqQ7d7/sEPPqDJkz8haaOk0/TWWw/of/83Hgq6C4elT33qXq1Y8bODATmsxsb/0s03Dz4g+0aMsIC1d6+F7+Ji6cknrQwgHM7dBXy7d9vYa2rs+fv0p+13dcop0ve/3/OZhOOOm5Pw+WydeqqF6ocflv7hH+ysxdatVp7y1a/awc8//IP0mc/YYs3Vq+3+P/xh6Q9/kC666NCw19Jis78zZqQuBNbX2+x/e/uhXx892kopIhEbSzIVFXYQkU22bElde0cAmUVNMpCHmpqkSy+1coYHH7QZ4ZYWm2WdOfPIP59MJGI72+3YcWjNaEuL9K1v2ay1ZOH0pJNstrq01ELYK6903+xhqaTrVVq6+u1d9HbvtsA20IDsC4elF16wmcR/+zcb03XXSTfcEC/hyCWeZ7Oqnmf/pq99zWqDJ0yw0oNkHSQ6O+33sGPHev3nf16rNWte0fTps5P2mo5GrZxiyRJblNfUZM9zXZ39n5kzx3ZWTNa6rKPDznLMm5f61ma7dll5R2JP4Zdftg1GRoyw30Pi/xV/973zzsuOswahkPTXvx5aagQgs6hJBnCI3/7Wws6558a7D4TD8Y8HoqTE6lcrKqwm1t9KuKHBaoqvusrKG+bPt4Cz+LDj4NckfVfSLyTFNH36XEmpC8iSlQrMmGFdIN7/fgvJf/yjdNNN8RnlvtQ5Z4vmZisfGDnSdrx75BGrBf6v/zo8IO/bZ/W49fV2+4aGSfrWt+YrErGv+TXE3RUX23N6wgn9G1c4bDO9c+emp/evvxA1Ejn0/8Qpp9gBwoYN9v8scZGev/tee3t6t3vvK0otgNxGSAbyTCgk/f739vEHPmDX7e12CnuwW+IWFdkCrcpKq5NtaIiHrxNOkL7zHZsxXrLEZpCjUen3v/83bdlyr6S1kiTnnGbOnKuvfe1eNTVZoDvhhMEHZN+oUfbvrKmxwLxypZVdzJ1roTOXduBbt85KbdrarLOEZDXc3et/Y7H4wrA5cw4Nh7GY9RdevTpejz7YwBaJWPg75ZS+1cYPRHGx/RvfeOPQ2mPnpPe+1xYSPv548k4WRUV2wJANIZlSCyC3ZcEJKQCp9Je/2CnnKVPiM4Tt7dIxx6Tm/p2z2bxTTrHgtX//od+vr5fe+U7bTvjDH5Z27LhdfkCWpJKSUv34x/NVVjZJjY2pmUFOHN+0aTauq66yrz38sIXmTZtS9zjptm+f1SJXV1t98O7dVirTfWtmz7OAPHHi4QFZssA4Zox01llWl97U1POudX0Rjdo4Tjgh/dsrjxqVfAe9Cy6w5/j555PvUllRYQcGQQuF7KCMkAzkLkIykGd+9zu7ft/7LExEo/G+xak0YoRt5lFaqrd3xUsmcQHZ1Km2i96UKRb6/K2zU2noUJvlnjvXZr1fe83C5J49h299nK02brTSkI0bre1ZUZG1eutea9vUZAF52rTea3BLS21W/aSTbKZ/IFugRyL2O5wxQxo7tv8/31/V1Xb2I3HHwNGjbQOccFh66qnDf66qyoJ80FtUU2oB5D5CMpBHNm6Uli2zYPjud9vX9u+3WcR0hNHqalvYNWWKhQJ/Z7Hubr/9Xs2aNVelpWWaPn2ubrnlXp1+ui3sS1eAcM7GFIvZIi5JevRRC5LdO3Rkq1DIylXq6myL7WjUFmJOnRq/zZ49Nts6dWrff49jxtiBTVGRBcm+rtsOh+32xx9vZxEyZcKEw7tcSNaGTkq+TXVRkT3vAzkQSKVt27KrFR2A/iMkA3nkscfs+rzz4qd5w+H01uEWF1spxzveYf1y9+yxGUe/dVlFxSTdccd8Pf54SI8+Ol/vf/+kt3eNS6dhwyzEX3CBff7YY3YqfvPm9D/2YDU1WfBdu9ZmS8vLbfGh78ABK1GZObP/XRxqamyG/eij7XGShdDuWlrsQOuUUwa38HMg/P8niQde555rv5MlS5K3fHPu8DKgTAqH7TWQ67s8AoWOhXtAnohG4yH5ve+Nf62kxEoP0q2qSjruONvNb98+C1/hsAXT6mqbFU3HbHZPnLOxtLVJ48ZZOF62zGaws21Xtu48z7qH1NZK3/iGfe3qq+M1wLGYlUzMmzfwTh0lJdYve8wYW4DZ1GT3VVFh34tG4y3eRo2yco4gamtLS+M78HWvt66psYOyv/zFOpjccMOhP1dRYSG1sTGz4/W1tNjzSKkFkNuYSQbyxOLFtmCpsTHeD7i93YJQJnvGlpbaDODRR1sgHTvWakszGZB9I0bYeC66yD5/5BELLrt3Z34sfdXaagF10ybp2WdtxvTaa+Pf37vXwv+QIYN/LL9ue+5cOwsg2WN7nh1YnHGG7eoX5OKzsWOT15FffLFd//nPh3+vsrJ/5SSpRlcLID8QkoE88eijdn3xxfFQHArZTGChKimxVmJnnmnh+NlnLTht2RL0yHq2bZvN6v7iF/b55ZfH26CFQvZvmjgxdY/nnIXl6dMtFL/znXY9dergWwamQkNDvP9xd6efbuNbv94u3RUX20LDxEV/mdDVZbPYhGQg9xGSgTzQ1iY984x97M+wSRYI09XLNleMGWMz23PmWIBZsMBOh2djl4tIxAJ8KGSlBEVF0gc/GP9+a6t1lygtDW6MmVZSYrPciQvxSkqks8+2j59+OvnPBrF4r7XVXnfZsOMfgMHhZQzkgWeesWB18snxOswDB2wWLpd2mEuHqioru/AD1ZNP2uzp3r3BjiuZvXttxvQPf7DAfPbZ8TKItjYrsUh3f+JsNGaM/f9O5HcuSRaSy8ribdgyadu25LsbAsg9hGQgD/j9Yt/znvjX/HpkxDc/KS622u1wODu7XGzaZDOkDz1kn19zjV17nj2f06YV5mKw+np77hK7XJx6qi3iW7fOtqruzq9LzqRoVNqxg64WQL4gJAM5bt8+6eWX7fTuO98Z/3oslprFXflgyBDbYOT00y3ILFhgs7bJZieDEgpZl4kFC2wG9Nhj7cyAZM/xmDGZ6VKSjUpKrMtFYvlEaWn8DIFfbuQrK7MDi8HsMNhfra3xzXsA5D5CMpDjnnvOTs3Pnh3fVc+viaypCXZs2aK42Hr8nnGGff7kk3YdxOn4nvibnPzmN3Z9zTU2a+x5UmendQopZEcdZb+HROeea9c91SUfqQ90Ku3cSXkTkE8IyUCO80st/LAgWT3y0KHMaHU3erSdni8rk1591WYlt20LelRxmzdbycDq1VZe4JfOtLZanXn3PsGFqL7eDvwS27qddpqVN6xde3jXElCx0e4AACAASURBVOdsFj4TYjH7/8SBKZA/CMlADtu/X1q0yMLDu94V/3pHR2G3fkumttaC8ty5FrTmz7dWXeFw0COzg5o9e6QnnrDPL7vMNsSQbHyTJgU3tmxRWmoLMBNnhsvLbWMRKX7A6KuoyFxd8v799lyVsEUXkDcIyUAO80stTj45Xmoh2axWobd+S2b8eNupTrIWa55n7eCC1twcb/smSZdeatf79lk3i0KfRfY1NibvfezX4j///KFfr6y0320mNhVpaiIgA/mGkAzkML8O02+FJcW3w2WF/eGGD7eSi4oKacUKK2XIhpKLTZvsjEBHh3TCCdaNQ7IZ5mOOCXRoWcVfuJgYek8/3cpoli+P13ZL8U1IktUyp5K/QQ0HM0B+ISQDOaqtTVq48PCuFqGQ7UTGrNbhKipsu+W5c+3zhQttsVUkEtyYOjosrPvbK/uzyO3tdnagUDtaJFNebr+TxI1gqqqsxZ/nSS+8cOj3/PZ56dTWZkG8kDZ5AQoBIRnIUQsXWnurWbNsRznfgQOFueFEX40bJ02btlOS9L3vLdL/+3/ztHTp+iP8VPo0N0tbt0pLl1p5wPnn29fb2+lokczYsclDr1+XnFhyUVKS/sV7zc3ssAfkI17WQI7yw4AfDnxdXcw+9mbIEOmZZ66R1C7pNK1evVU33nhtYOPZtEl69ln7+PzzbVY0FLLr7nXmMEOGJK8x9l8HixYdWl5RUXFoCUY6bNlCVwsgHxGSgRwUiUgvvWQfJ4ZkiTfs3pSWSm+++ZKkRw9+5SqtWvWKotHMj6WjwzY1SSy12LfPZpGZnTxcVZX9/07cCGbECOm44+zrixbFv15RYb/jdC3e6+iwzhZ+NxIA+YM/wUAOWr7c6liPPto6NvgiEavb5A27dyedNEfSAwc/e58mT54dSJeL5mbplVfsevx4W7QXiVg4pmSmZ42Nh+++J8V333vuufjXioqs20uyrhipQKkFkL94aQM5yC+1OOss62Th6+w8tD4Zyd13372aOnWP/JKLj370N9q+PfPj2LQp3qHk0kvtuWxtlSZOZBFYb4YPV9KZfz8kv/DC4d9PXOyXKlu20EkGyFeEZCAH9VSPTEjum2OPnaSHHvqrzjnH0s1bbx2tbduSB690aW+XNm60TU2Ki6WLL7aSgGjUZkrRs9paO2OS2JXkmGPsd7d3r7X48xUXp2fxXmenHdRUVqb+vgEEj5AM5JiNG+1SV2en57vzPOqR+2rMmPjGIk8/beE0kyUXzc22YC8atXEMH261raNHW90teuZc8pIL55J3uUjX4r09e1J/nwCyByEZyDF+H9gzzji8F7LnEbD6qr7eNqGorJRWrbLZx0yWXGzcGN9h77LL7LqzM76RCHo3YkTyLcWT1SVXVNgBUKoX723eTKkFkM8IyUCO6V6P3F1XlwXksrLMjykXFRfb6fnTTrPPFyyw3fcysbFIW5st2NuwwVqanXmmBeTqatr39VV9fXxHve5OPNHKMTZssAMRyRbWpXrnvVDIDqwotQDyFyEZyCGtrdJrr1k48EsFfAcO0Fe3v8aMsRl5SXrmGQtSe/em/3F37ZKeeso+vvBCOyOwf7+F9u4LMdGz4mIrTUncWKSkxA46pENLLjwvtR0u/FILni8gfxGSgRwyf74FudmzD689DoUIyf1VX29bVPslFy0tdgo9nTxPev31eDnApZfac0rbt/4bMyb57HCyuuTiYjsQSZUtWyhtAvIdIRnIIT2VWvioj+yfoiJp0iSrTZZsq++mpuS1rqmyf7/NWre1SdOn26Yhfs9r2r71T0ODzeQm1hrPnWuh+LXX4osxKypSt9AuFLKFgIRkIL8RkoEc0dVlM8lSzyGZN+3+Gz06XnLhl0Ds3p2+x9u+Pf44l1xi15GINHZs+h4zX5WW2gK+xJKLmho72xKLxXem9HfeSwVKLYDCQEgGcsSSJRYGJk06PFBFIrZgj0V7/eeXXFRVWclFa6v01lvpeaxo1BbsvfqqBbz3vMdmlIcNo3XfQDU2Jq819ksu/G4wxcV2hiAVi/c2beKAFCgEhGQgR/hv9n6Lq+46O61LAvqvqMgWzM2da58//bRtPJHK+lVfS4v05z9becA551hA7+iwAx8MjP//PrHkwj/bsmCBnYWRbOZ3sDvvHThgpRaUNgH5j5AM5ADP670emZ32Bmf0aOld77KPH3/cgnM6eiZ334b6kkssvFVUsOByMMrL7feXGH4bG63eu73dZu8lC8mJG5D0V3Oz/f8AkP8G/VJ3zo1zzv3VObfaObfSOXdzKgYGIO7NN62H79Ch0owZh3+fnfYGp65OOvVU6y6xdauVW2zcmNqeyZ2dtnnItm32OKedZqUdkyYRugZr3LjD65Kl+AGlf4CZisV7GzdaH2YA+S8Vf5ojkj7red50SadL+qRz7rgU3C+Ag/w3+TPOsNrKZNjUYOCcs53u/NnkP/3JAvKuXal7jJ07pcces48vuSTelWH06NQ9RqEaOtR+l4klF35p0gsv2PfKywe3eK+tzUpxKioGfh8AcsegQ7Lneds9z1ty8OP9klZLahzs/QKI660eORazmUjeuAdn1Kh4SH76aVsEuX59arYyjsWkl1+2TgvFxdIVV1jYGjfOghsGxy9ZSSy5OO44K0Pavl1at84WS3Z2xmuU+2vnzp4PUgHkn5Se5HPOTZB0kqRFqbxfoJA1N0srVlho87dQ7i4ctnIB2lENTk2NNHWqBav2dlvw1doa77M7GHv2SH/4g3W3eMc7LJCHQhaSkRrJSi6KiuK77/mbt0gD23kvFrNSi7q6gY8RQG5JWUh2ztVIekjSpz3P25fk+zc55xY75xa3tqbwHCaQ51580WYz58xJXlLR2WmbKmBwnLMNPfzZ5EcesTZfb745+Ptes0Z64gn7+P3vtzA3ZAiBK5WGDUtecuHXJftnY5xLXr98JC0tdmDDhi9A4UhJSHbOlcoC8n2e5/0u2W08z7vL87w5nufNqa8fkYqHBQqCX4/s931NFIkQklNl5EgLVeXlVh6xZ4/VJe877LC/71pbrRa5udkW6c2ebSFt8uTUjRv2nI0ceXgAPu00+97KlbZJTFnZwOqSN2+mpAkoNKnobuEk/UTSas/zvj34IQHwdXZKiw4WL/W0y14sxqK9VKmqstnkc8+1zx980ALWunUDv8833ogv2Lv66njbN1r2pd64cYeXUlRWSqecYh+/+OLAOlx0dlpXErpaAIUlFTPJZ0i6TtK7nHNLD14uSsH9AgXvb3+zN+jp062OtSfs/pU648dLF1xgH99//369+91D9MEPztOrr67v1/2sX79es2adrFNOOVXLl0tFRW068cS31Npqs8i0fUu9oUNtYV00eujX/bMwzz1nM8nt7f1r79fUZGUa1P0DhSUV3S1e9DzPeZ43y/O8Ew9eHk/F4IBC5y82StbVQrI3+vJy6iRTafhwacoUqapqmWKxWkWj12nNmgX60IeuVSzW9/u59tprtXz5q5JukSTFYv+r//iPD6qoiLZv6VJSYlu2J+6W6J+FefllO+j0vL4v3ovFrC69vj61YwWQ/ZjLALJULBZfbNRTPXIoRD1yqpWW2m5tnZ1fO/iVz0kq0Zo1r2jTpr7fz+LFiyVNl3SFpE5J39GaNa9o4kQOatKpsdE6vnQ3YoR1LQmFLCj3Z+e95mYL1mVlqR8rgOxGSAay1MqV9gY9Zox07LHJb0NITo/GRmnKlB2SVkk6WtIHNX36bK1a1bdFfB0d0uTJcyR9/eBXfiJppyZPnq2xY9M1akjWMaS62l4b3XXvclFRYa+tvnjzTbs/AIWHkAxkKb+rxdln91wLGY2yHXU6NDRIX/nKL9XY+KAkqbT0P/Uv/3KfamulV16xmcWeRCLSq69K73//byVdKalD0n9o8uST9O1v38siyzRzzrqIJB7M+Gdjnn/eSpR27z7yfe3da4v8CMlAYSIkA1nKr0fuqdTCR+hKPeeks86apB/96Cs65hipq2uMFiyYqKoqq2d9+eXkvXbDYWnJEquJ/c1vhh28rx9o5szxuvXWB3XWWZMy/C8pTCNH2nPYfQHflCm2+LW52TqOdHb2frAjWVcTAjJQuAjJQBbavNm2RK6pkU4+uffbEpLTY+RIWwj2yU/a53ffbT2T6+utXvzFF6UNGywQt7dbi7AXX7QZzBdekNasqZC0Q573H1q+fIG+//1rCVwZUlZmrfxaW+Nfc+7Q2eQjbSqyd68935ypAQoXIRnIQv4s8hlnWFBLxu9s0dP3MThlZdZ39/jjbWvjtjbpG9+wmeTaWivJWLtWeuklC8XLltkBS0eH9P3v+/dysyRLaqtXvxLUP6UgjRtnr5HuO/D5IfmFF6xVXE+bisRi0qpVBGSg0BGSgSzUvR65J2xHnX7jx1vQuuUWO+3+3HPSb35j3yspsXZxI0bYJRRar5tuOkeXXbZEbW1STc1fJT3w9n3Nnj07mH9EgaqpkY466tDa5Nmzraf466/bQc/mzUra1m/bNvs5Zv6BwkZIBrJMS4u0dKmFsHnzer5dOEzv1nSrrragVVkpffnL9rXvfEd66qnDb3vbbTdoxYovSjpZ0nqNHv1VzZw5VyUlZTrttLm69957Mzl0SJo4Md4XWbKzA3Pn2scvvWTf616SIVl4XrFCGjIks2MFkH0IyUCWeeklm92aPbv3073RKNvkZoIftM4/X/roR+25ueUWK6loabHPly6VVqz4jqQLJe2SdJk2blyob35zvlatCmnhwvmaNIlFe5lWV2cHOd2D8DvfadfPPGOhedu2+PdCIXsuKyvpZQ1AopoRyDJ97WohWb9XpFddne3i1tQk3XST/c5/8APpF7+wS3m535P3JEkbJF0saZWmTp2rkhJpwoQABw8de6wF4VjMtgI/80wLwEuX2tmYDRvsOa6qst7kXV2UMQEwzCQDWeTAAWn+fPu4t3pkHyE5MyZPttrkWEz68Iet08W8ebb4KxSy1mJXXNGiGTM+otLSdZo1a64+/el7dfzxzEgGrbpaOuaY+CK9mhrp9NOtBOP5562ufMUKa+vneQRkAHHMJANZxK+TnDlTGj2659tFInaqmACWGVVV1md37VprDXfiiVZuEY1aG7HaWsm5BklPS7INKIYPt9sieMccY7PJnZ12YHniiU164YWRuuOOp/T441/R7bffq7FjKYcBcChmkoEs8rRlLJ17bu+3YzvqzBs/3k7L798f/1pxsX2t+46IHR226HLGjJ53SkRmlZRIJ5xgtcnRqPT00x+SFJHnnaNly9boK1+5NughAshChGQgS3R22mYU0pFDcjhs4QyZU1wsnXSS1az2tFNbZ6eF5JNPtpl+ZI+hQ6WpU2076tdff0bSX2UnUy+lhzWApAjJQJaYP99qkmfMkMaM6f22kQidLYJQVSWdcooF4cTWYa2tVnpx6qm05stWEydaCcyUKXMkPXTwq1dq+nR6WAM4HCEZyBJ9LbXwsWgvGA0NthNiba11vNi9265ra61zAv11s1dRkdX7f+lL92ratM2SonLuPfrCF34V9NAAZCEW7gFZoLPTtsqVCMm5wJ9Rbm+38ouyMvsasl95uXTBBZM0ZMhj+upXpb/9rVhr107QtGlBjwxAtmEmGcgCCxfaKfzp06XGxt5vG41afSw1r8GrrraZZQJybhk2zMouzjrLPn/iiWDHAyA7EZKBLPDnP9v1eecd+bZdXX7LsfSOCchnxx4rnXaadb5YvNjKZgCgO0IyELD9+21TA+ekCy448u1DIRaGAYNVW2sdL+bOtU1innoq6BEByDaEZCBgTz9tLd3mzLGd246kq4v2b8BgFRVZFxl/Z0tKLgAkIiQDAXv8cbu+6KK+3d7zpMrK9I0HKBSjRllP68pKaflyacuWoEcEIJsQkoEAbd8uLVliK+7f9a6+/xydLYDBq6+31xKzyQCSISQDAfrTn+z6nHOsU8KReJ7VLhOSgcErLZVqauILZh97zF5jACARkoHAeF7/Sy26uqzdWBGvXCAlhg6VTjhBGj5c2rRJWro06BEByBa81QIBWb5c2rDB3qRPO61vPxMOs2gPSKUhQ2yb90susc//+MdgxwMgexCSgYA8+KBdX3KJ9WrtC0IykFrV1XZW59JL7fOnnpLa2oIdE4DsQEgGAtDSYm/GzklXXNH3n4tG+1a7DKBv/N0Sx42zThedndKTTwY7JgDZgZAMBOCRR2xWeO7cI29DnYhFe0DqlJbagWc4LF12mX3t4YeDHROA7EBIBjIsFpMeesg+vvLK/v88IRlIrSFDbAb53HMtMK9YIa1bF/SoAASNkAxk2Msv26YFo0ZJZ57Z95+LxaTiYqmsLH1jAwrRkCE2k1xRIV14oX3tgQeCHROA4BGSgQz77W/t+u/+zkJvX4XDUm1tesYEFLKqqnh/5Pe/364fe0xqbQ1uTACCR0gGMuitt6Tnn7c6yMsv79/PhkK2QxiA1KqoiIfkiRNtrUAoJP3+98GOC0CwCMlABv3sZ/ZmfMkltnlBf3R10f4NSIeKCjur4wflD3zArn/7W+uhDKAwEZKBDNm2Tfrzn+3N+EMf6v/Pex6L9oB0KCqy7alDIfv89NOl8eOlnTulv/412LEBCA4hGciQX/7S+hy/+93S2LEDuw9CMpAeDQ1W9y9ZaPZnk3/1q+DGBCBYhGQgA3bvjm93e/31A7+f8vKUDAdAgvr6eEiWpIsvtvKm5culxYuDGxeA4BCSgQz4+c/tDfid75SOOab/Px+JWEDu6/bVAPqne4cLSaqslD74Qfv4zjsP/R6AwkBIBtJs61ZbAOScdOONA7uPcJjOFkA6JStluuYae90tXSotWpT5MQEIFiEZSLMf/tBmgi+8UJo6dWD3EQ7bwiIA6VFebgey3WeMq6ul666zj3/0I2aTgUJDSAbS6JVXpCeesDfgf/zHgd8P7d+A9PI7XHSvS5ak973PduRbvlyaPz+YsQEIBiEZSJNIRLrjDvv4+uulMWMGd38s2gPSq67u8JBcVRVv2fj979M3GSgkhGQgTe65R1q/3tq9DaQvciLavwHp1dAQ75Xc3dVXS42N0ptvSvffn/lxAQgGIRlIgzVrpJ/8xD6+7bbBzQL7dZDMJAPpldjhwldRIX3+8/bxXXdJu3ZldlwAgkFIBlKsrU360pds45D3vU+aM2dw99fVZQuIini1AmnV29maM8+Uzj5bam+XvvOdzI0JQHB42wVSKBaTvv51adMm6dhjpX/+58HfJ50tgMw4UknT5z5nZ3SefFJasCAzYwIQHEIykEJ33ik9/bTN/P7nf6amjpgeyUBmFBfbJiJdXcm/P2aM9NGP2sdf+5rU0pK5sQHIPEIykCL33WeL9YqLpW98Q5owITX3G40ykwxkSl1d8sV7vuuuk044wbaa/9rX7OwRgPyUkpDsnLvAObfWObfOOXdLKu4TyBWeZ9tO+3WKt94qzZ2b2sdg0R6QGfX1h7eB66642EqqamulF16wA2MA+WnQIdk5VyzpfyVdKOk4SR9wzh032PsFckEoZGUV//M/9vktt0iXX576x6H9G5AZNTV29qY3Rx0l/du/2Q59d95pGwYByD8lKbiPUyWt8zxvvSQ5534j6TJJq3r6gQ0bpCuvPPRrzh1+u8SvJbvNkX6mr/czkPtmzP27n6IiqbT00EtZmVRSEv+4tNTaMFVXH3qpqrKZm6FDsycwrlkj3X679PrrNu7bb5fOPz+1j+F59nssK0vt/QJIrq9/X844Q/rUp2yDkX/9V/v7NG9eeseWbp4ndXZKra3Svn3x685OmxDo7Dz0EgrZJRKxA4u+XLo/ln/d/eOBfD/xe0B3Rzro7U0qQnKjpM3dPt8i6bTEGznnbpJ0k302Wxs3puCRUZAqKy0sDx1q28UOHSqNGGGLavzLqFEWXNNh2zbppz+V/vhHq0dsbLTZ5OPScP7E72zRlwMXAIPXn4Pw666znsm//rX02c/aWoSzz07f2AYjFpN27JA2b5aamqSdO+OXpiZbhLhvX++lJkChSUVITvb2fdgxned5d0m6S5ImTJjjfetb3b+X5A683j/v6WupuJ+B3O9AH6sv950vY45GbdV44iUctpmIcNguHR3Wi9S/9j/et0/as0c6cEDautUuPXFOGj780ODsX446Sho9un9vhlu2SAsXSs8/b62fPM9mxj/wAenjH7fZ7nQIh6Vhw9Jz3wAO55/dikat/rg3zkmf+Yzd9oEHbMORT3/a/i4EdWAbi1kLyjVrbMfPjRvtsnlz7wsSfWVlVpddX2+LGOvqbGKiosIu5eXxjysq4r+v4mK7+B8XFR36uf81Kf676X7d/fd1pO8nfi3xGuhu7954V5r+SkVI3iJpXLfPx0ra1tsPlJWlbuU/CovnWWjes8cue/dKzc02E7J9u82UbN9un+/aZZdly5Lf19ChFpYbGizk1tTYJRazcNrWFg/jzc3xnystlc47T/rIR9L//zgctjcpAJlTV2evvcrKI9/WOQvHDQ22G9+3vy0tWWILeNN9gBuJWPnimjXxy+uv26RCMsOGSUcfbX/3Ro2SRo6061Gj7KxcfX32lLMBqVJbO/CfTUVI/pukY51zEyVtlXSNpL9Pwf0Ch3EuHmaPPrrn20UiFpC3bYuHZ/9j/3M/aPdFba106qnWteLss+0NJRNiMavHBpA5dXV2cNyXkCzZ36WbbpImTbIFfc8+Ky1ebAfSV13V9/vpjefZ364VK6SVK+2yenXy2eFRo6Rp06TJk+1A/uijpfHjaSUJ9NegQ7LneRHn3D9JekJSsaR7PM9bOeiRAYNQUhIvr0gmGrU+pzt2SPv326yxf/EXElZWWmnG2LFW83ykU6/pwswOkFn19TZD21/nnSfNmGG1yS+9JH3ve9LPfiZddJF01lnWX7kv7RyjUQvpb71lM8N+KN679/DbNjZK06dbKJ42TZo6NXMH8UC+S8VMsjzPe1zS46m4LyATiovjpxmzHT2SgcwazIHpmDHSd78rzZ8v3X23zfz++td2KSuzIOuXO/iPE4nYwjn/wH3jxuQL6BoaLITPmCEdf7wtFm5oGPhYAfQuJSEZQOrFYhbmaf8GZNZgz944Zy3izjjD6oT//Gfp5ZdtVnjZsp7XSXQ3apQ0caKVcPjBuLGRxWlAJhGSgSwVDg9uwQGAgSkvt04Mfp/ywfDLICSbLV63Lt52LRy2+y8qshnhYcOstGvCBOqHgWxASAayVDhsq88BZFZRkS2Y7epK7ZmchgZpzpzU3R+A9Br0ttQA0oP2b0BwamvZWAModIRkIEt5XmpaRwHov/r6vm2+ASB/EZKBLEb7NyAY1dW2eBZA4SIkA1mM9m9AMDhABUBIBrJQJGILhkpLgx4JUJgIyQAIyUAW6uqiBRQQpNJSO1CNRIIeCYCgEJKBLBQK0dkCCFpdHR0ugEJGSAayUCTCRiJA0AjJQGEjJANZiPZvQPDq6qz0CUBhIiQDWYrOFkCwWLwHFDZCMpCleIMGgsVrEChshGQgy0QiNotcXBz0SIDCVl4uOWflTwAKDyEZyDKhkG2JCyBYzlkrRhbvAYWJkAxkma4uOlsA2YIOF0DhIiQDWYaQDGSP+no7uwOg8BCSgSzjeSwYArJFdTU1yUChIiQDWcY52r8B2YIDVqBwEZKBLOJ5hGQgm/BaBAoXIRnIIl1dUlWVVMQrE8gKJSUWlCORoEcCINN4KwaySDhsq+kBZA8W7wGFiZAMZBF6JAPZhzZwQGEiJANZJBazcgsA2aOujnILoBARkoEsw2p6ILtUVNAGDihEhGQgyxCSgexSUWFdZwAUFkIykCViMam4WCorC3okALorK7MuF9Fo0CMBkEmEZCBLhMNsRw1kKxbvAYWHkAxkCTpbANmrro42cEChISQDWaKri5lkIFvV19PhAig0hGQgS3ieVFkZ9CgAJEOHC6DwEJKBLEJnCyA78doECg8hGcgivBED2clvA8dsMlA4CMlAFujqslKL4uKgRwIgGedszQCL94DCQUgGskAoZKvnAWSv+nrawAGFhJAMZAE6WwDZj5AMFBZCMpAFIhFCMpDtqqqoSQYKCSEZyBIs2gOyG69RoLAQkoEswRswkN3ocAEUFkIyELBYzLpalJUFPRIAvSkqkmpqqEsGCgUhGQhYOGz1yM4FPRIAR1JXR0gGCgUhGQhYKGSr5gFkv4YGeiUDhYKQDASsq4seyUCuoMMFUDgIyUDAPM922wOQ/XitAoWDkAxkATpbALmBDhdA4SAkAwHyPHvDJSQDuYEOF0DhGFRIds590zm3xjm3zDn3e+dcQ6oGBhSCri6putreeAHkBjpcAIVhsG/Nf5F0vOd5syS9LunWwQ8JKBx0tgByDx0ugMIwqJDsed6TnudFDn66UNLYwQ8JKByEZCD3VFdTkwwUglSe5L1B0p9SeH9A3vM8aykFIHewhgAoDCVHuoFz7ilJo5N868ue5/3x4G2+LCki6b5e7ucmSTdJ0siRRw9osEA+oqUUkFsqKmwdgb/wFkB+OmJI9jzvvN6+75z7sKT3SjrX83o+AeV53l2S7pKkKVPmcKIKOIhZKSC3+B0uQiFev0A+G2x3iwskfVHSpZ7ndaRmSEBh6OqyWeTi4qBHAqC/GhrocAHku8HWJP9AUq2kvzjnljrn7kzBmICCEArZGy2A3FNfT0gG8t0Ryy1643ne5FQNBCg0dLYAcldVFR0ugHzHFgZAQGIxq2sEkHtYcAvkP0IyECAW/QC5qbzc1hPEYkGPBEC6EJKBADEbBeQm56xcip33gPxFSAYCEA7brl10tgByF9tTA/mNkAwEIByW6uqCHgWAwaDDBZDfCMlAADo7pSFDgh4FgMFgS3kgvxGSgQB4npVbAMhdlZVsSw3kM0IyEBA6WwC5raTEgnJXV9AjAZAOhGQgwzzPZp/ol5YmrwAADjxJREFUbAHkPhbvAfmLkAxkWDhsm4gU8eoDct7QobbGAED+4W0ayDAW7QH5o7qa7amBfEVIBjIsHLbWUQByHx0ugPxFSAYCQGcLID/421NHo0GPBECqEZKBALBoD8gPzln5FIv3gPxDSAYyKBq1tlFlZUGPBECqNDSweA/IR4RkIINCIXtDZQMCIH/U11NuAeQjQjKQQXS2APIPi/eA/ERIBjIoGpXq6oIeBYBUqqy0vue0ggPyCyEZyCDPY9YJyDdFRXbwy+I9IL8QkoEM8Tx7M6WzBZB/2HkPyD+EZCBDQiGptpbtqIF81NAgdXUFPQoAqcTbNZAhoRCL9oB8RRkVkH8IyUCGhMOEZCBfVVWxeA/IN4RkIIOYbQLyE4v3gPxDSAYyiJAM5K9hw6QDB4IeBYBUISQDGRAOS9XVtiU1gPzU0CBFIkGPAkCqEJKBDDhwwFpEAchfnCkC8gshGcgAFu0B+c9fvBeLBT0SAKlASAYywPOs3AJA/nKOTUWAfEJIBjLAOUIyUAiGD2fxHpAvCMlAmoVCUk0Ni/aAQlBXJ0WjQY8CQCoQkoE0O3DAWkMByH+cMQLyByEZSLOuLjpbAIWivFyqrLTXPYDcRkgGMoDZJaBwjBhBXTKQDwjJQBp5nrWEon8qUDiGDWN7aiAfEJKBNOrslOrrLSgDKAw1NXaADCC38dYNpNGBA3bqFUDhqKqSiovZVATIdYRkII2iUamhIehRAMikoiL6JQP5gJAMpBE77QGFicV7QO4jJANp0tVlp13Ly4MeCYBMq62lLhnIdYRkIE06OuyUK4DCU1Nj29ETlIHcRUgG0iQcJiQDhaqkRBoyxDrcAMhNhGQgTTzPZpMAFKaRI+2MEoDcREgG0iASkcrK2EQEKGQNDdbhBkBuIiQDaXDggJVaOBf0SAAEpabG2sFRlwzkJkIykAadnWwiAhQ6vy6ZVnBAbiIkA2ngeVJdXdCjABC00aOpSwZyFSEZSLFoVCotZRMRAFaXzPbUQG5KSUh2zn3OOec552h4hYLX0WGr2qlHBlBTY2UXLOADcs+gQ7Jzbpyk8yVtGvxwgNzX2SmNGhX0KABkg6IiWsEBuSoVM8nfkfQFSazfBWT1yLW1QY8CQLYYPZpNRYBcNKiQ7Jy7VNJWz/Ne68Ntb3LOLXbOLW5t3TWYhwWyVleXVFFBf2QAcfX1tIEDclHJkW7gnHtK0ugk3/qypC9JendfHsjzvLsk3SVJU6bM4c8F8lJ7u9TYGPQoAGSTigo7u9TZaR8DyA1HDMme552X7OvOuZmSJkp6zdkKpbGSljjnTvU8b0dKRwnkiHDY6g8BoLvGRumNNwjJQC45Ykjuied5yyW9HQeccxskzfE8b3cKxgXkHM+zRTr0RwaQaOhQOlwAuYY+yUCKdHTYVtQlAz70BJCvamulsjIpEgl6JAD6KmUh2fO8Ccwio5AdOCCNGRP0KABko6Ii6aijpLa2oEcCoK+YSQZSJBaThgwJehQAstXo0bZuAUBuICQDKdDZabXIlZVBjwRAtqqrY/c9IJcQkoEUaG+Xxo0LehQAsllxsXW52L8/6JEA6AtCMpAC0aitXgeA3owZQ8kFkCsIycAgdXbayvWamqBHAiDb1dfT5QLIFYRkYJDa2qTx44MeBYBcUFRkpVn79gU9EgBHQkgGBikWs/7IANAXo0czkwzkAkIyMAgdHdb2ja4WAPqqrs4unZ1BjwRAbwjJwCC0tUkTJgQ9CgC5ZuJEulwA2Y6QDAxQNGotnSi1ANBfw4dLztEzGchmhGRggPbtswU4JSVBjwRArikrk44+WmptDXokAHpCSAYGqKvLNgYAgIEYN44FfEA2IyQDA9DWJg0bZotvAGAgamqkESOoTQayFSEZGICODmnSpKBHASDXHXOMdOBA0KMAkAwhGeinzk6pupptqAEM3pAhdmlrC3okABIRkoF+2rdPmjrVds4CgMGaMkVqbw96FAAS8TYP9IM/izxiRNAjAZAvhg6VRo5kq2og2xCSgX5obZWmTWMWGUBqTZlitcmxWNAjAeDjrR7oo/37bQMAZpEBpFpdne3euXdv0CMB4CMkA30Qjdosz7RptksWAKTa5Mm2i2coFPRIAEiEZKBP9uyxNzD6IgNIl7IyaeZMm032vKBHA4CQDBxBW5tUW0tfZADpN3KkNHGi1Nwc9EgAEJKBXnR1WUeLE06w06AAkG5TpthZq9bWoEcCFDZCMtCDSMRmc0480baPBYBMKCmxvzsS/ZOBIBGSgSSiUWn3bun446VRo4IeDYBCU1kpnXqqLeIjKAPBICQDCbq6LCDPmCGNHx/0aAAUqpoa6fTT7aC9pSXo0QCFh5AMdNPWZm9GJ59sPUsBIEi1tRaUa2qkpiYLzAAyg5AMyGaPd+60WsAzzpBGjw56RABgKiulOXNsQd/u3XYgT4s4IP1KgnpgmqXnhr5unNH9dsl+xrn417t/HCTPkzo6rN6vrEyaNUs66ii2nAaQfYqLpWOOsQP4N9+Utm61v1V1dfb3K1t4XjzAJ153v01vn/d230B/DebsSyAhuSSwaI5ER/qj09v3u/8x9MVih3/Pv47F7JLs55yzr/nXPfEDdnGxvUEkuySGcM+zx41EbMY4HI4/xrBhNjszYgQt3gBkv+pqO6CfPFnasUPavNlmlp2zv3/l5fYe6/+N9P8e+n93/Yv/99i/RKP2df96oPzH6/7Y3f8e+5MQySZLuk+kDPSxgUT19ZIUjQ3kZwOJq1VV0rx5QTwyskn3P9aJ18n+gPvX4bBduroODb6RiF26h3E/UJeW2pvHkCE281JdbbV+paVB/xYAoP+qqmyDo0mTpAMH7IxYW5u0f799HArZ38vuobeoyP4e+iG6rMz+BpaW2tf8z7uH7GQTEokhONnkBJBdIpGB/BRzughM9z+wAICBqay0y/DhQY8EyC/EEwAAACABIRkAAABIQEgGAAAAEhCSAQAAgASEZAAAACABIRkAAABIQEgGAAAAEhCSAQAAgASEZAAAACABIRkAAABIQEgGAAAAEhCSAQAAgASEZAAAACCB8zwv8w/q3H5JazP+wEiV4ZJ2Bz0IDAjPXW7j+ctdPHe5jecvt031PK+2vz9Uko6R9MFaz/PmBPTYGCTn3GKev9zEc5fbeP5yF89dbuP5y23OucUD+TnKLQAAAIAEhGQAAAAgQVAh+a6AHhepwfOXu3juchvPX+7iucttPH+5bUDPXyAL9wAAAIBsRrkFAAAAkCAjIdk5d7VzbqVzLuac63F1qHPuAufcWufcOufcLZkYG47MOTfUOfcX59wbB6+H9HC7qHNu6cHLw5keJ+KO9FpyzpU75+4/+P1FzrkJmR8lkunDc3e9c25Xt9fajUGME8k55+5xzjU551b08H3nnPv+wed3mXPu5EyPEcn14bk7xznX2u2195VMjxHJOefGOef+6pxbfTBv3pzkNv1+7WVqJnmFpCskPd/TDZxzxZL+V9KFko6T9AHn3HGZGR6O4BZJT3ued6ykpw9+nswBz/NOPHi5NHPDQ3d9fC19RNJez/MmS/qOpDsyO0ok04+/g/d3e639OKODxJH8TNIFvXz/QknHHrzcJOn/MjAm9M3P1PtzJ0kvdHvt3Z6BMaFvIpI+63nedEmnS/pkkr+d/X7tZSQke5632vO8I20ecqqkdZ7nrfc8LyzpN5IuS//o0AeXSfr5wY9/LunyAMeCI+vLa6n7c/qgpHOdcy6DY0Ry/B3McZ7nPS9pTy83uUzSLzyzUFKDc25MZkaH3vThuUOW8jxvu+d5Sw5+vF/SakmNCTfr92svm2qSGyVt7vb5Fh3+D0QwRnmet12y/4iSRvZwuwrn3GLn3ELnHEE6OH15Lb19G8/zIpJaJQ3LyOjQm77+Hbzy4OnCB51z4zIzNKQI73W5ba5z7jXn3J+cczOCHgwOd7B88CRJixK+1e/XXsp23HPOPSVpdJJvfdnzvD/25S6SfI3WGxnS2/PXj7s52vO8bc65SZKecc4t9zzvzdSMEP3Ql9cSr7fs1Jfn5RFJv/Y8L+Sc+7jsjMC70j4ypAqvvdy1RNJ4z/PanHMXSfqD7NQ9soRzrkbSQ5I+7XnevsRvJ/mRXl97KQvJnuedN8i72CKp+4zIWEnbBnmf6KPenj/n3E7n3BjP87YfPDXR1MN9bDt4vd4596zsSI6QnHl9eS35t9ninCuRVC9OM2aDIz53nuc1d/v0blFPnmt4r8tR3UOX53mPO+d+6Jwb7nne7iDHBeOcK5UF5Ps8z/tdkpv0+7WXTeUWf5N0rHNuonOuTNI1kuiQkB0elvThgx9/WNJhZwacc0Occ+UHPx4u6QxJqzI2QnTXl9dS9+f0KknPeDRNzwZHfO4SauguldXeIXc8LOlDB1fany6p1S9nQ3Zzzo321244506VZajm3n8KmXDwefmJpNWe5327h5v1+7WXspnk3jjn/k7S/0gaIekx59xSz/Pe45w7StKPPc+7yPO8iHPunyQ9IalY0j2e563MxPhwRN+Q9IBz7iOSNkm6WpKctfP7uOd5N0qaLulHzrmY7A/HNzzPIyQHoKfXknPudkmLPc97WPbH5Jfu/7d3hygVAEEch38TvIHZi1jMBpPBImi0eACL1/AG3sBzCFaj2ahFeAa1rEUs74Xvi5sGhoU/wyw789zXBPlsexXz44+9u56Zk75ec79WF1srmF9m5r46qvZn5qW6rfaqNpvNXfVQHVfP1Vt1uZ1KWf2hd6fV1cx8VO/VmeHCzjiszqunmXn8PrupDur/d8+PewAAsNildQsAANgJQjIAACyEZAAAWAjJAACwEJIBAGAhJAMAwEJIBgCAhZAMAACLT/rOQGFoNQSnAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plotting function \n", "def plotgp(m, X, y, minx, maxx):\n", " xx = np.linspace(minx, maxx, 500)[:,None]\n", " mean, var = m.predict_y(xx)\n", " plt.figure(figsize=(12, 6))\n", " plt.plot(X, y, 'k.', mew=2)\n", " plt.plot(xx, mean, 'b', lw=2)\n", " plt.fill_between(xx[:,0], mean[:,0] - 2*np.sqrt(var[:,0]), mean[:,0] + 2*np.sqrt(var[:,0]), color='blue', alpha=0.2)\n", " plt.xlim(minx, maxx)\n", "\n", "\n", "opt = gpflow.optimizers.Scipy()\n", "\n", "# data\n", "np.random.seed(1)\n", "N = 20\n", "X = np.random.rand(N,1)\n", "Y = np.sin(12*X) + 0.66*np.cos(25*X) + np.random.randn(N,1)*0.1 + 3\n", "\n", "# fit GP fit Gaussian kernel\n", "kgauss = gpflow.kernels.RBF( 1, lengthscales=0.3 ) # init lengthscale to 0.3\n", "mg = gpflow.models.GPR( (X,Y), kernel = kgauss )\n", "mg.likelihood.variance.assign(0.01)\n", "opt.minimize(mg.training_loss, variables=mg.trainable_variables)\n", "\n", "# print and plot\n", "gpflow.utilities.print_summary(mg)\n", "plotgp(mg, X, Y, -1, 2)\n", "\n", "\n", "# fit GP with SM kernel\n", "ksm = UnivariateSMkernel( Q = 3 )\n", "ms = gpflow.models.GPR( (X,Y), kernel = ksm )\n", "ms.likelihood.variance.assign(0.01)\n", "opt.minimize(ms.training_loss, variables=ms.trainable_variables)\n", "\n", "# print and plot\n", "gpflow.utilities.print_summary(ms)\n", "plotgp(ms, X, Y, -1, 2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Task 3c__ Apply the SM kernel to the CO2 dataset. Below is a template code that (usually) fits the datasets correctly due to initialising the correct 1-year frequency (1/12) for the model. It might take a few retries to find a good fit.\n", "\n", "* Modify the code to have a random initialisation of the frequencies. Can you still learn the 1-year frequency if the initialisation is random? Why or whynot?\n", "\n", "* Modify the visualisation to include both training and test points. Is the fit reasonable?\n", "\n", "* Visualise the optimised SM kernel matrix. What can you see?\n", "\n", "* Visualise the spectral density $S(\\omega) = \\sum_{q=1}^Q a_q \\mathcal{N}(\\omega | \\mu_q, \\sigma_q^2)$ using the parameter you just optimised. \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sklearn\n", "from sklearn.datasets import fetch_openml\n", "\n", "def load_mauna_loa_atmospheric_co2():\n", " ml_data = fetch_openml(data_id=41187)\n", " months = []\n", " ppmv_sums = []\n", " counts = []\n", "\n", " y = ml_data.data[:, 0]\n", " m = ml_data.data[:, 1]\n", " month_float = y + (m - 1) / 12\n", " ppmvs = ml_data.target\n", "\n", " for month, ppmv in zip(month_float, ppmvs):\n", " if not months or month != months[-1]:\n", " months.append(month)\n", " ppmv_sums.append(ppmv)\n", " counts.append(1)\n", " else:\n", " # aggregate monthly sum to produce average\n", " ppmv_sums[-1] += ppmv\n", " counts[-1] += 1\n", "\n", " months = np.asarray(months).reshape(-1, 1)\n", " avg_ppmvs = np.asarray(ppmv_sums) / counts\n", " avg_ppmvs = np.asarray(avg_ppmvs).reshape(-1, 1)\n", " return months, avg_ppmvs\n", "\n", "\n", "X,Y = load_mauna_loa_atmospheric_co2()\n", "\n", "Xtr = X[0:100]\n", "Ytr = Y[0:100]\n", "\n", "Xts = X[100:]\n", "Yts = Y[100:]\n", "\n", "\n", "# set SM kernel model\n", "ksm = UnivariateSMkernel( Q = 3, freqs = [0.083, 1.0, 0.5], ells = [0.01, 0.01, 0.01], amps = [20,20,20] )\n", "#ksm = UnivariateSMkernel( Q = 10 )\n", "m = gpflow.models.GPR( (Xtr,Ytr), kernel = ksm )\n", "m.likelihood.variance.assign(0.1)\n", "set_trainable(m.likelihood.variance, False)\n", "\n", "# optimize\n", "opt = gpflow.optimizers.Scipy()\n", "opt.minimize( m.training_loss, variables=m.trainable_variables )\n", "\n", "# print\n", "gpflow.utilities.print_summary(m)\n", "\n", "# plot\n", "plotgp(m, Xtr, Ytr, 1950., 2020.)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }