misc/gaussian_process_regression_basic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "05DE16FF3E9B483788D3B192E7AF6F73"
   },
   "source": [
    "# Gaussian Process Regression in Pytorch\n",
    "\n",
    "*Thomas Viehmann*, <tv@lernapparat.de>\n",
    "\n",
    "Bayesian time!\n",
    "\n",
    "I could not write a Gaussian Process introduction better than [Rasmussen and Williams in the canonical book](http://www.gaussianprocess.org/gpml/), so you will just follow the link.\n",
    "\n",
    "This notebook touches selected aspects of Chapter 2 (Regression) and Section 5.4 (Model Selection for GP regression, we use the marginal likelihood).\n",
    "\n",
    "As usual, the best part is that Pytorch will do the gradients for us.\n",
    "\n",
    "(Note that the explanations are in no way finished.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "cell_id": "C56FA1DC11504021833EAF55BEF3E07D"
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot\n",
    "%matplotlib inline\n",
    "import IPython\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch.autograd import Variable\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "7D0580E460BC4D258B714A255A53612C"
   },
   "source": [
    "Let's have a very cheap cholesky layer.\n",
    "\n",
    "If you want to be a hero, write a autograd triangular solve function. I think that would be both more efficient and more stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "cell_id": "9BF8B14EE53A47A08E2A7082A831858F"
   },
   "outputs": [],
   "source": [
    "class Cholesky(torch.autograd.Function):\n",
    "    @staticmethod\n",
    "    def forward(ctx, a):\n",
    "        l = torch.potrf(a, False)\n",
    "        ctx.save_for_backward(l)\n",
    "        return l\n",
    "    @staticmethod\n",
    "    def backward(ctx, grad_output):\n",
    "        l, = ctx.saved_variables\n",
    "        # Gradient is l^{-H} @ ((l^{H} @ grad) * (tril(ones)-1/2*eye)) @ l^{-1}\n",
    "        # TODO: ideally, this should use some form of solve triangular instead of inverse...\n",
    "        linv =  l.inverse()\n",
    "        \n",
    "        inner = torch.tril(torch.mm(l.t(),grad_output))*torch.tril(1.0-Variable(l.data.new(l.size(1)).fill_(0.5).diag()))\n",
    "        s = torch.mm(linv.t(), torch.mm(inner, linv))\n",
    "        # could re-symmetrise \n",
    "        #s = (s+s.t())/2.0\n",
    "        \n",
    "        return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "B0752DA711874B3A9E068AB2F952EC8A"
   },
   "source": [
    "So let's do Gaussian Process Regression. We use the basic RBF kernel but with automatic relevance detection and a gaussian likelihood (i.e. noisy data).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "cell_id": "8E992D5F44DA489C89BDEF3E536700EA"
   },
   "outputs": [],
   "source": [
    "LOG2PI = numpy.log(2*numpy.pi)\n",
    "class ARDRBFKernel(nn.Module):\n",
    "    def __init__(self, n, lambda0=None, sigma0=None):\n",
    "        super(ARDRBFKernel, self).__init__()\n",
    "        self.n = n\n",
    "        if lambda0 is None:\n",
    "            lambda0 = torch.ones(n)*0.05\n",
    "        self.lam = nn.Parameter(lambda0)\n",
    "        if sigma0 is None:\n",
    "            sigma0 = torch.ones(1)*0.1\n",
    "        elif numpy.isscalar(sigma0):\n",
    "            sigma0 = torch.FloatTensor(1).fill_(sigma0)\n",
    "        self.sigma = nn.Parameter(sigma0)\n",
    "    def forward(self, x, x2=None):\n",
    "        # todo stabilize\n",
    "        x = x/self.lam.unsqueeze(0).expand_as(x)\n",
    "        if x2 is None:\n",
    "            x2 = x\n",
    "        else:\n",
    "            x2 = x2/self.lam.unsqueeze(0).expand_as(x2)\n",
    "        x_  = x.unsqueeze(1).expand(x.size(0),x2.size(0),x.size(1))\n",
    "        x2_ = x2.unsqueeze(0).expand(x.size(0),x2.size(0),x.size(1))\n",
    "        sqdists = ((x_-x2_)**2).sum(2) #.squeeze(2)\n",
    "        return self.sigma.view(1,1).expand_as(sqdists)*(-sqdists/2).exp()\n",
    "\n",
    "class GP(nn.Module):\n",
    "    def __init__(self, K, x, y, var0=None):\n",
    "        super(GP, self).__init__()\n",
    "        if var0 is None:\n",
    "            var0 = torch.ones(1)\n",
    "        elif numpy.isscalar(var0):\n",
    "            var0 = torch.FloatTensor(1).fill_(var0)\n",
    "        self.var = nn.Parameter(var0)\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.K = K\n",
    "    def predict(self, xstar=None):\n",
    "        Kxx_inv = torch.inverse(self.K(self.x)+Variable(torch.eye(self.x.size(0)))*self.var)\n",
    "        Kxxstar = self.K(self.x, xstar)\n",
    "        Kxstarxstar = self.K(xstar)\n",
    "        mu = torch.mv(Kxxstar, torch.mv(Kxx_inv,self.y))\n",
    "        cov = torch.Kxstarxstar-torch.mm(Kxxstar.t(), torch.mm(Kxx_inv), Kxxstar)\n",
    "        return mu, cov\n",
    "    def forward(self, xstar=None):\n",
    "        Kxx_noise = self.K(self.x)+Variable(torch.eye(self.x.size(0)))*self.var\n",
    "        L = Cholesky.apply(Kxx_noise)\n",
    "        Linv = L.inverse()\n",
    "        alpha = torch.mv(Kxx_noise.inverse(),self.y)\n",
    "        if xstar is not None:\n",
    "            Kxxstar = self.K(self.x, xstar)\n",
    "            Kxstarxstar = self.K(xstar)\n",
    "            mu = torch.mv(Kxxstar.t(), alpha)\n",
    "            v = torch.mm(Linv,Kxxstar)\n",
    "            cov = Kxstarxstar-torch.mm(v.t(),v)\n",
    "        #print (self.K.lam.data[0], (0.5*self.y.dot(alpha)).data[0],L.diag().log().sum().data[0],self.y.size(0)/2.0*LOG2PI)\n",
    "        neg_logp = 0.5*self.y.dot(alpha)+L.diag().log().sum()+self.y.size(0)/2.0*LOG2PI\n",
    "        if xstar is not None:\n",
    "            return mu,cov,neg_logp\n",
    "        else:\n",
    "            return neg_logp\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "582380B3DDE44415B06F652C6940AF71"
   },
   "source": [
    "Let's have a regression example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "cell_id": "7D9AA2647C2C4AD6A42AA374EFF07F28"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa246b15668>]"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEyxJREFUeJzt3X2MY1d9xvHnIVleqk1JtWuLKNnNFBFQzagl1EqDEG0c\nUJWmKJFoXAWJNyntKmGmBYFU0SKlkKoCKhUqtJ7SRUENlELWAdFtlPwRyl3xombpbNiEiZdEGwrN\niqg2G9iwKqQs/fUPe6zZycz6esb2HZ/5fiQr1+Ozvj+fzDxz59xzz3VECACQlucVXQAAYPQIdwBI\nEOEOAAki3AEgQYQ7ACSIcAeABBHuAJAgwh0AEkS4A0CCLixqx7t3746ZmZmidg8AU+no0aM/jIjS\noHaFhfvMzIwWFxeL2j0ATCXb38/TjmEZAEgQ4Q4ACSLcASBBhDsAJIhwByBJajQaarfb/eftdluN\nRqPAirAZhc2WAbB1NBoNzc/Pa2FhQVmWSZJqtZparZYkaW5ursjysAGEOwDV63UtLCyo1WppdnZW\nktTpdFSpVFSv1wuuDhvBsAwAlctlZVmmUqmkTqejTqejUqmkLMtULpeLLg8bQLgDQIIIdwBqt9uq\n1Wr9I/blI/harXbOSVZMD8IdgJrNplqtliqVipaWlrS0tKRKpaJWq6Vms1l0edgATqgC6M+Gqdfr\n/TH2LMvUbDaZKTOlHBGF7LharQYLhwHAcGwfjYjqoHYMywBAggh3AEgQ4Q5gXSxJML04oQpgTSxJ\nMN0IdwBrYkmC6cawDIA1sSTBdCPcASBBucPd9gW2v2X73jVee4Htu22fsH3E9swoiwQweSxJMN2G\nOXJ/l6Tj67x2i6QfRcTLJH1M0kc2WxiAYrEkwXTLdULV9mWSfl/SX0t6zxpNbpT0gd72PZL223YU\ndfkrgE1jSYLplne2zN9J+jNJF63z+qWSnpSkiDhr+7SkXZJ+uOkKARRmdYiXy2WCfUoMHJax/UZJ\n7Yg4er5ma3ztOUfttvfZXrS92Ol0higTADCMPGPur5V0g+3vSfq8pGtt/9OqNicl7ZEk2xdKerGk\np1e/UUQciIhqRFRLpdKmCgcArG9guEfEn0fEZRExI+lmSV+JiLesanZI0tt72zf12jDeDgAF2fA8\nd9t32L6h9/ROSbtsn1D3hOv7RlHcSqxxAQD5DbX8QEQclnS4t337iq//TNLYrkdmjQsAGM5UrC3D\nGhcAMJypWH6ANS4AYDhTEe4AgOFMRbizxgUADGcqwp01LgBgOFNxQpU1LgBgOC7qWqNqtRqLi4uF\n7BsAppXtoxFRHdRuKoZlAADDIdwBIEGEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLc\nE8LdqgAsm4q1ZTAYd6sCsBLhngjuVgVgJYZlEsHdqgCsRLgDQIII90RwtyoAKxHuieBuVQBW4oRq\nIrhbFYCVuBMTAEwR7sQEJIQL1DCsgeFu+4W2v2n7YduP2v7gGm3eYbtj+1jv8UfjKRfYfpYvUFs+\nOb588nx+fp6Ax7ryjLk/K+naiDhje4ekr9u+PyIeXNXu7oiYH32JwPbGBWrYiIFH7tF1pvd0R+9R\nzEA9sA1xgRo2IteYu+0LbB+T1Jb0QEQcWaPZH9h+xPY9tves8z77bC/aXux0OpsoGwBwPrnCPSJ+\nERGvknSZpKtsz65q8q+SZiLi1yV9WdJd67zPgYioRkS1VCptpm5g2+ACNWzEULNlIuLHkg5Lum7V\n109FxLO9p5+U9JsjqQ4AF6hhQ/LMlinZvri3/SJJb5D0nVVtLlnx9AZJx0dZJLCdzc3Naf/+/f0x\n9mazqYMHD2r//v2am5tjWiTWlGe2zCWS7rJ9gbq/DA5GxL2275C0GBGHJP2p7RsknZX0tKR3jKtg\nYDtavsp4eVpkpVJRlmX9IRvW7cdqXKEKTJGVYb583mp5WiSzZ7YHrlAFEsS0SORFuANAggh3YIow\nLRJ5Ee7AFGFaJPJiPXdgirBuP/JitgwATBFmywDANka4A0CCCHcASBDhDkwJbrWHYTBbBpgCy2vK\nLCwsKMsySWJNGZwX4Q5MAW61h2ExLANMAdaUwbAIdwBIEOEOTAHWlMGwCHdgCrCmDIbFCVVgCrCm\nDIbF2jIAMEVYWwYAtjHCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACRoYLjbfqHtb9p+2Pajtj+4\nRpsX2L7b9gnbR2zPjKNYAEA+eY7cn5V0bUT8hqRXSbrO9tWr2twi6UcR8TJJH5P0kdGWCQAYxsBw\nj64zvac7eo/Vl7XeKOmu3vY9kl5v2yOrEgAwlFxj7rYvsH1MUlvSAxFxZFWTSyU9KUkRcVbSaUm7\nRlkoACC/XOEeEb+IiFdJukzSVbZnVzVZ6yj9OYvW2N5ne9H2YqfTGb5aAEAuQ82WiYgfSzos6bpV\nL52UtEeSbF8o6cWSnl7j3x+IiGpEVEul0oYKBgAMlme2TMn2xb3tF0l6g6TvrGp2SNLbe9s3SfpK\nFLXcJAAg13rul0i6y/YF6v4yOBgR99q+Q9JiRBySdKekz9g+oe4R+81jqxgAMNDAcI+IRyRducbX\nb1+x/TNJ3IIdALYIrlAFgAQR7gCQIMIdABJEuANAggh3AEgQ4Q4ACSLcASBBhDsAJIhwB4AEEe4A\nkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7AExAo9FQu93uP2+322o0GmPbX547MQEANqHRaGh+\nfl4LCwvKskySVKvV1Gq1JElzc3Mj3yfhDgBjVq/XtbCwoFarpdnZWUlSp9NRpVJRvT6em9gxLAMA\nY1Yul5VlmUqlkjqdjjqdjkqlkrIsU7lcHss+CXcASBDhDgBj1m63VavV+kfsy0fwtVrtnJOso0S4\nA8CYNZtNtVotVSoVLS0taWlpSZVKRa1WS81mcyz75IQqAIzZ8myYer3eH2PPskzNZnMsM2UkyREx\nljcepFqtxuLiYiH7BoBpZftoRFQHtWNYBgASRLgDQIIGhrvtPbYz28dtP2r7XWu0ucb2advHeo/b\nx1MuACCPPCdUz0p6b0Q8ZPsiSUdtPxARrVXtvhYRbxx9iQCAYQ08co+IpyLiod72TyQdl3TpuAsD\nAGzcUGPutmckXSnpyBovv8b2w7bvt/3KEdQGANig3PPcbe+U9AVJ746IZ1a9/JCkyyPijO3rJX1J\n0hVrvMc+Sfskae/evRsuGgBwfrmO3G3vUDfYPxsRX1z9ekQ8ExFnetv3Sdphe/ca7Q5ERDUiqqVS\naZOlAwDWk2e2jCXdKel4RHx0nTYv6bWT7at673tqlIUCAPLLMyzzWklvlfRt28d6X/sLSXslKSI+\nIekmSbfZPivpp5JujqIufQUADA73iPi6JA9os1/S/lEVBQDYHK5QBYAEEe4AkCDCHQASRLgDQIII\ndwBIEOEOAAki3AEgQYQ7ACSIcAeABBHuOEej0VC73e4/b7fbajQaBVYEYCNyL/mL9DUaDc3Pz2th\nYUFZlkmSarWaWq3uTbfm5uaKLA/AEAh39NXrdS0sLKjVaml2dlaS1Ol0VKlUVK/XC64OwDAYlkFf\nuVxWlmUqlUrqdDrqdDoqlUrKskzlcrno8gAMgXAHgAQR7uhrt9uq1Wr9I/blI/harXbOSVYAWx/h\njr5ms6lWq6VKpaKlpSUtLS2pUqmo1Wqp2WwWXR6AIXBCFX3Ls2Hq9Xp/jD3LMjWbTWbKAFPGRd0N\nr1qtxuLiYiH7BoBpZftoRFQHtWNYBgASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABA0Md9t7\nbGe2j9t+1Pa71mhj2x+3fcL2I7ZfPZ5yAQB55LlC9ayk90bEQ7YvknTU9gMR0VrR5vckXdF7/Jak\nv+/9FwBQgIFH7hHxVEQ81Nv+iaTjki5d1exGSZ+OrgclXWz7kpFXCwDIZagxd9szkq6UdGTVS5dK\nenLF85N67i8AAMCE5A532zslfUHSuyPimdUvr/FPnrNoje19thdtL3Y6neEqBQDklivcbe9QN9g/\nGxFfXKPJSUl7Vjy/TNIPVjeKiAMRUY2IaqlU2ki9AIAc8syWsaQ7JR2PiI+u0+yQpLf1Zs1cLel0\nRDw1wjqxTTUajXNuFNJut9VoNAqsCJgOeY7cXyvprZKutX2s97je9q22b+21uU/SdyWdkPRJSe8c\nT7nYDpYDvdFoaH5+Xq973ev04Q9/uH+nqPn5eQIeGID13LGlLAd6pVLRwYMH9aY3vUmPP/64JGnn\nzp06c+aMKpUKN+3GtsV67phK9Xq9f2u/Wq2mU6dO9V87c+aMSqUSwQ7kQLhjSymXy8qyrH9z7lOn\nTql72gfAMAh3bHkRoV27dvUDv1arnXOSFcBzEe7YUpZPmnY6He3cubP/9d27dyvLsv6QTbPZLLBK\nYOsj3LGlNJtNtVotVSoVPfHEE/rQhz6kV7ziFXrsscd0+PBhZVmm/fv3a25uruhSgS2N2TLYchqN\nhur1ev+kabvdVrPZJNAB5Z8tQ7gDwBRhKiQAbGOEOwAkiHAHgAQR7gCQIMIdABJEuANAggh3AEgQ\n4Q4ACSLcASBBhDsAJIhwB4AEEe4AkCDCHQASRLgDQIIIdwBIEOEOAAki3AEgQYQ7ACSIcAeABA0M\nd9ufst22vbTO69fYPm37WO9x++jLBAAM48Icbf5R0n5Jnz5Pm69FxBtHUhEAYNMGHrlHxFclPT2B\nWgAAIzKqMffX2H7Y9v22Xzmi9wQAbFCeYZlBHpJ0eUScsX29pC9JumKthrb3SdonSXv37h3BrgEA\na9n0kXtEPBMRZ3rb90naYXv3Om0PREQ1IqqlUmmzuwYArGPT4W77Jbbd276q956nNvu+AICNGzgs\nY/tzkq6RtNv2SUl/KWmHJEXEJyTdJOk222cl/VTSzRERY6sYADDQwHCPiDcPeH2/ulMlAQBbBFeo\nYiwajYba7Xb/ebvdVqPRKLAiYHsZxWwZ4ByNRkPz8/NaWFhQlmWSpFqtplarJUmam5srsjxgWyDc\nMXL1el0LCwtqtVqanZ2VJHU6HVUqFdXr9YKrA7YHhmUwcuVyWVmWqVQqqdPpqNPpqFQqKcsylcvl\nossDtgXCHQASRLhj5Nrttmq1Wv+IffkIvlarnXOSFcD4EO4YuWazqVarpUqloqWlJS0tLalSqajV\naqnZbBZdHrAtcEIVI7c8G6Zer/fH2LMsU7PZZKYMMCEu6mLSarUai4uLhewbAKaV7aMRUR3UjmEZ\nAEgQ4Q4ACSLcASBBhDsAJIhwB4AEFTZbxnZH0vcl7Zb0w0KK2Dq2ex9s988v0QcSfSDl64PLI2Lg\nrewKC/d+AfZinmk9KdvufbDdP79EH0j0gTTaPmBYBgASRLgDQIK2QrgfKLqALWC798F2//wSfSDR\nB9II+6DwMXcAwOhthSN3AMCITSzcbV9n+zHbJ2y/b43XX2D77t7rR2zPTKq2Scjx+d9ju2X7Edv/\nZvvyIuocp0F9sKLdTbbDdnIzJ/L0ge0/7H0vPGr7nydd47jl+FnYazuz/a3ez8P1RdQ5LrY/Zbtt\ne2md1237473+ecT2qze0o4gY+0PSBZKekPRSSc+X9LCkyqo275T0id72zZLunkRtW+jz1yT9Um/7\ntpQ+f94+6LW7SNJXJT0oqVp03QV8H1wh6VuSfqX3vFx03QX0wQFJt/W2K5K+V3TdI+6D35b0aklL\n67x+vaT7JVnS1ZKObGQ/kzpyv0rSiYj4bkT8r6TPS7pxVZsbJd3V275H0utte0L1jdvAzx8RWUT8\nT+/pg5Ium3CN45bne0CS/krS30j62SSLm5A8ffDHkhoR8SNJiojUbl2Vpw9C0i/3tl8s6QcTrG/s\nIuKrkp4+T5MbJX06uh6UdLHtS4bdz6TC/VJJT654frL3tTXbRMRZSacl7ZpIdeOX5/OvdIu6v7lT\nMrAPbF8paU9E3DvJwiYoz/fByyW93PY3bD9o+7qJVTcZefrgA5LeYvukpPsk/clkStsyhs2LNU3q\nTkxrHYGvnqaTp820yv3ZbL9FUlXS74y1osk7bx/Yfp6kj0l6x6QKKkCe74ML1R2auUbdv96+Zns2\nIn485tomJU8fvFnSP0bE39p+jaTP9Prg/8Zf3pYwkiyc1JH7SUl7Vjy/TM/9U6vfxvaF6v45dr4/\nXaZJns8v22+Q9H5JN0TEsxOqbVIG9cFFkmYlHbb9PXXHGg8ldlI178/Bv0TEzyPiPyU9pm7YpyJP\nH9wi6aAkRcS/S3qhumuubBe58mKQSYX7f0i6wvav2n6+uidMD61qc0jS23vbN0n6SvTOLiRg4Ofv\nDUn8g7rBnto4qzSgDyLidETsjoiZiJhR97zDDRGR0r0Y8/wcfEndk+uyvVvdYZrvTrTK8crTB/8l\n6fWSZPvX1A33zkSrLNYhSW/rzZq5WtLpiHhq6HeZ4Bni6yU9ru6Z8vf3vnaHuj/AUvd/YFPSCUnf\nlPTSos9qT/jzf1nSf0s61nscKrrmSffBqraHldhsmZzfB5b0UUktSd+WdHPRNRfQBxVJ31B3Js0x\nSb9bdM0j/vyfk/SUpJ+re5R+i6RbJd264nug0eufb2/054ArVAEgQVyhCgAJItwBIEGEOwAkiHAH\ngAQR7gCQIMIdABJEuANAggh3AEjQ/wODCPxB2Js5WQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa24709a780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 12\n",
    "X = torch.rand(N,1)\n",
    "Y = (torch.sin(12*X) + 0.6*torch.cos(25*X) + torch.randn(N,1)*0.1 + 3).squeeze(1)\n",
    "pyplot.figure()\n",
    "pyplot.plot(X.numpy(), Y.numpy(), 'kx', mew=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "cell_id": "37D848D2C2374C19894E433C3C5F1A58",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final loss 15.55954647064209\n"
     ]
    }
   ],
   "source": [
    "K = ARDRBFKernel(1)\n",
    "gp = GP(K, Variable(X), Variable(Y), var0=0.1)\n",
    "\n",
    "def closure():\n",
    "    opt.zero_grad()\n",
    "    nlogp = gp.forward()\n",
    "    nlogp.backward()\n",
    "    return nlogp\n",
    "\n",
    "opt = torch.optim.LBFGS(gp.parameters(), lr=0.1)\n",
    "for i in range(10):\n",
    "    opt.zero_grad()\n",
    "    nlogp = gp.forward()\n",
    "    nlogp.backward()\n",
    "    opt.step(closure)\n",
    "print (\"Final loss\",nlogp.data[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "EF88856D8CAF446380AC9BFEADCE4647"
   },
   "source": [
    "And we can plot with the credible interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "cell_id": "0912147330F74C12B0E6E1E38FBD8C60"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa244e59240>]"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXlcW3X2//86gYQ9YSesYStQ6GqZ6oyjUh21blW0am3r\nVrXjNq3185txGUc/o99Z9KP2Y8c6H2uta2u1Kuo4LqNjGR1nbG1rV1qWAmELhJ1AgJDk/fsDwlAK\n5Ca5yc3yfj4ePBrgnXvPLcnJued9zusQYwwcDofD8S9kUhvA4XA4HPHhzp3D4XD8EO7cORwOxw/h\nzp3D4XD8EO7cORwOxw/hzp3D4XD8EO7cORwOxw/hzp3D4XD8EO7cORwOxw8JlurE8fHxLDMzU6rT\nczgcjk+yf//+DsZYgr11kjn3zMxM7Nu3T6rTczgcjk9CRFoh63hahsPhcPwQ7tw5HA7HD+HOncPh\ncPwQ7tw5HA7HD+HOncPhcPwQ7tw5HA7HDxHk3ImonoiOENFBIjqtfpFG2URENUR0mIjOEN9UDofD\n4QjFkTr3JYyxjml+dwmAWWNfZwL489i/HA6Hw5EAsZqYrgTwOhsdyPodEUUTUTJjTCfS8b2WwcFB\ndHZ2wmAwYHh4GIwxyOVyREZGIiYmBkqlEkQktZkcDifAEOrcGYC/ERED8CJjbMuk36cCaJzwfdPY\nz/zSuTPG0NXVhcbGRvT09Ey5pqOjA/X19QgLC0N6ejrUajVkMr7FweFwPINQ5342Y6yFiBIBfEFE\nJxhjX0/4/VShKZv8AyJaC2AtAGRkZDhsrDdgNBpRXV2N7u5uQesHBwdRVVWF5uZm5OfnQ6lUutlC\nDsc3MJvN6Ovrw8DAAEwmExhjCA4ORlhYGJRKJcLCwqQ20acR5NwZYy1j/+qJqAzAYgATnXsTgPQJ\n36cBaJniOFsAbAGA4uLi05y/N8MYg06nQ01NDaxWq8PPHxgYwA8//ICcnBykpqbyVA0nIGGMoaOj\nAzqdDt3d3RjN5E5NeHg41Go1kpOTIZfLPWilf2DXuRNRBAAZY8ww9vgiAI9PWvYRgHuJaCdGN1J7\n/SnfbrVaUVVVhdbWVpeOwxhDTU0NjEYjZs2axR08J2BgjKG9vR11dXUYHBwU9Byj0Yja2lpotVpk\nZGQgPT2dpzYdQEjkngSgbMwRBQPYwRj7jIjuBADG2P8B+ATApQBqABgB3Ooecz3PyMgIjh49it7e\nXtGO2dLSArPZjNmzZ3MHz/F7hoaGUFlZKTiVORmLxYK6ujq0tbVh9uzZiIqKEtlC/4Rmui1yJ8XF\nxczbJX9NJhMOHz6M/v5+txxfrVYjPz+fO3iO39Le3o7KykqYzWZRjkdEyMvLQ3JysijH80WIaD9j\nrNjeOsn03L0dk8mEgwcPwmg0uu0cra2tCAkJQVZWltvOweFIAWMM9fX10GoFSY87dNzKykoYjUZk\nZ2fzwGgGeAJrCkZGRnDo0CG3OnYbWq0Wer3e7efhcDyF1WpFZWWl6I59Io2NjaiqqppxQzbQ4c59\nEhaLBUeOHMHAwIDHznnixAmPno/DcRdWqxUVFRUuFx8IQafTobq6mjv4aeDOfQJWqxXHjh1DX1+f\nJOe1WCwePS+HIyY2x97RMZ1Kifi0tLSgvr7eY+fzJbhzH4MxhqqqKnR1dUlyfqPRiJMnT0pybg7H\nVWzvH086dhtarRY6nd9UXosGd+5jaLVaj9xKzkRLS4tkHy4cjivU1dVJ+v6pqqoStVzZH+DOHaNV\nK95yaydm2RiH4wl0Oh0aGhoktYExhoqKCphMJknt8CYC3rl3d3ejsrJSajPGGR4e9poPGg7HHr29\nvaiqqpLaDACj750TJ07wDdYxAtq5DwwM4NixY173YmhqaoLBYJDaDA5nRoaHh73u/dPV1YXm5map\nzfAKAta5Dw8P48iRI16bAuElXhxvxlYZ441pkNraWo/0qHg7AenczWYzjhw5gqGhIalNmZa+vj7e\n3MTxWurr6712A9NqtfL0DALQudtqyt2lFyMmtbW1vPad43V0d3dLvoFqj76+voBPzwSUc7darTh+\n/LjT6nSeZnh4GE1NTVKbweGMYzKZcPz4canNEERdXZ1X3527m4Bx7jbBofb2dqlNcYiGhgavzGty\nAg9bo5KvvB4tFktANwYGhHO35eDa2tqkNsVhLBaL198CcwKDtrY2STpQXaG9vT1gGwP93rlbLBZU\nVFT4pGO30dzcjOHhYanN4AQww8PDqKmpkdoMp3B2NKav49fO3TZsw9eijckwxtwqn8rhzIQtHeOt\nZcP2MBqNAak947fO3WAw4MCBA15bruUoOp0uoDeHONKh1+vR2dkptRkuUV9f77MfTs4i2LkTURAR\n/UBEH0/xu1uIqJ2IDo593S6umcJhjKGxsREHDhzwK2fIGOO5d47HMZlMPpuOmcjIyEjAvX8cGbO3\nHsBxAMppfv82Y+xe101yHoPBgOrqao/rsXsKnU4HjUaDkJAQqU3hBAgnT57EyMiI1GaIQlNTE1JT\nUwPm/SMocieiNACXAdjqXnOcw2Aw4NixY9i/f7/fOnZgNHrnde8cT9HV1eXThQiTsVqtAbV3JTRy\n/18AvwIQNcOaa4joXABVADYwxhonLyCitQDWAkBGRoaDpv4HxhgGBwfR2dkJvV4fUCJbLS0tyMjI\ngFwul9oUjh9jsVi8Ru1RTHQ6HTIyMhAaGiq1KW7HbuRORJcD0DPG9s+w7C8AMhlj8wB8CeC1qRYx\nxrYwxooZY8UJCQlOGQyM3iru3bsXJ0+eDCjHDoy+6QK9rZrjfhoaGvxqz8pGIFWeCUnLnA1gGRHV\nA9gJ4HwienPiAsZYJ2PMVoj9EoBFolo5iUCsWZ1Ic3Mz15zhuA2j0ejXm4+tra0YHByU2gy3Y9e5\nM8YeYoylMcYyAawA8BVjbPXENUSUPOHbZRjdeOW4iZGREb/KhXK8B1tNuz8rKgZK5ZnTde5E9DgR\nLRv7dh0RHSOiQwDWAbhFDOM409PU1OTXb0CONOj1evT09EhthttpbW31y7TTRBwphQRjrBxA+djj\nRyf8/CEAD4lpGGdmjEYjurq6EBcXJ7UpHD/BbDYHjNCWrR9m1qxZUpviNvy2QzUQ4GWRHDGpr6/3\nGcVHMdDpdH59vdy5+zDd3d0YGBiQ2gyOHzAwMBBwwYLVavXryjPu3H2cQHtDcsSHMYbq6mqpzZAE\nf648487dx2lra/Ob9nCONATKJupUmM1mv1WM5M7dx7FarWhtbZXaDI6PEkibqNPhr5Vn3Ln7Ac3N\nzX754uS4H61WO76pWFZWdsp84e7ubpSVlUllmscYGhryufGbQnCoFDLQsViA48eV+P77WFRVRaGx\nMQwGgxwWCyEy0oyEhGHk5Rkwb14vFi/uRFiYZzpph4aGeFkkx2EmbqKWlZVh06ZN+PDDD7Fx40YA\nwIYNG8Zb9UtLSyWz0xM0NTUhMTFRajNEhTt3AfT2BuPDD1Px178mQ68PhUzGkJk5gNzcfkRHjyA4\nmKGvLxhtbaH45JNkvP9+GkJDLTjvvHZcf30DsrKMbrexpaWFO3eOYBhjqKmpGb/jKykpwYcffgit\nVos1a9YAAHp6eqDRaFBSUiL4uFYroNeHoLExHH19cgwNySCXMyiVI0hKGkZGhhFBQd53l9nX14e+\nvj4oldMpmvse3LnPwNCQDO++m4adOzMwMBCM4uIu/PzntSgu7oJSOfVUF7OZcPSoEn//exK+/DIJ\nn3+uRkmJHnfddRKJie6bg9rZ2YmhoaGAULvjuE57e/spKZiYmBhs3LgRa9asGd9cjY6OxsaNGxET\nEzPjsQyGYPzjHwnYuzcWBw7EYGBgereiUFiQn2/Aj3/ciXPO6UBamvdovDQ1NaGwsFBqM0SDO/dp\nOHpUiSefLEBTUzh+8pMO3H57raAIPDiYYcGCXixY0Ivbb6/Fe++l4Z130vHdd3G47bZaXHNNM4jc\nY7NOp0NWVpZ7Ds7xGywWiyibqHV1Edi5Mx3l5QkwmYKQlDSEkhI98vIMyMgwIjZ2BCEhFpjNMvT2\nytHcHIrq6igcPBiNLVtysGVLDhYu7EZpaTPOPrsDMol3ANvb2zE8POw3wzxIqo244uJitm/fPqee\nW1VVhZaWFpEtGsVqBV59NRPbt2uQmDiEX/6yEmec4VqZWGtrCJ57Lg/ffReHs87qxAMPnEB0tPjl\niwqFAmeddRZkUr9LOF5NbW3tacJZ3d3d4zn26OhoAP9Jy0yO3ltbQ/HSS1nYvTsRoaEWXHhhGy67\nTIdZs/oFBy5tbSH48sskfPRRCvT6UOTmGnDHHXX40Y+63Bb8CCEjIwPZ2dnSGSAAItrPGCu2t457\ngQkMDAThN7+ZgzfeyMRFF7Xi5Zf3uezYAUCtHsbvf38E69ZVYf/+GNx55yLU14eLYPGpmEwmnx9k\nzHEvAwMDaGw8bY4OysvLodVqodFosG3bNmzbtg0ajQZarRbl5eUAAIuF8M47abj11h/hX/+Kxw03\nNOCtt77Dhg3VyMsT7tgBIClpGKtWNWDHjj14+OEKDAwE44EH5uGRR+agvV26yFmn0/lNUxOP3Mfo\n6pLjl7+cj/r6CNx7bzWuuqrFLRFEZWUkHn54LoaHg/DEE0excKG4zSOxsbGYN2+eqMfk+AeMMRw6\ndGjahqWysjKUlJSMR+nd3d0oLy9HaWkp9PoQPP54IY4dU+EnP+nA+vXVou4hjYwQ3nsvDa++momg\nIIZf/KIaS5dKI2tdUFAAtVotybmFwCN3B2hrC8H69QvR0hKGJ588jNJS9zh2AMjP78fmzQcQHz+M\nBx6Yh++/n3mzylG6urr8XsqU4xzt7e0zdqKWlpaekn6xPf7qK4Y77ihGXV0ENmz4BsXFT4heHCCX\nM6xY0YiXX/4eubn9ePLJ2XjyyXwMDXneRfmL3kzAO/f2dgXuu28BursV+J//OYTi4m77T3IRtXoY\nzz33AzIyjPj1r+eK7uD9tZ2a4zxmsxk1NTUOPcdW+/7EE/ciJqYJTz31Bd5//0Zs2rTJbc1NqalD\nePbZQ7jxxnp8/rka99xzBvR6z6ZpDAYD+vr6PHpOdxDQzr2vLxi/+tV89PbK8cwzhzBnjuf+oCqV\nGc88cxAZGUY88sgcVFSIV1/b2trKO1Y5p1BXV+eQvK3VCmi1NwMoBFCBnp4z8MgjK8fz8o7UvjtK\nUBDDmjX1+OMfD6OtLRR3330Gqqsj3Xa+qfCH6D1gnfvgoAwPPzwXzc1h+N3vjiI/3/ODtlUqM55+\n+hDi4014+OE5aG4OE+W4w8PD6OrqEuVYHN/HYDA45KysVuCZZ/Lx4YdnYOnSHYiOjkZvbw96enoE\n176LweLF3fjTn35AUBDDunULsX9/tNvPaUOv1/u81rtg505EQUT0AxF9PMXvQojobSKqIaI9RJQp\nppFiwxjw5JMFOH5ciUceqRB9U9MRoqNH8OSTh8EY4Ve/moeeHrkox+WpGQ7wn5moQrFagY0b8/DJ\nJ8m48cZ63HFHrRuts09W1gA2bz6A5ORBPPyw+CnM6WCM+bwgnyOR+3pMP/j6NgDdjLFcABsBPOmq\nYe5k+/YM/OMfibjjjlqce26H1OYgLW0Qv//9EbS3j1YkWCyu7+Z2dnb6fOTBcZ2WlhYYDMLuShkD\nnn8+Fx9/nIJVq7S46qofcP/9G8Yj9ujoaPT09GDDhg2ndLe6m/h4EzZuPIT09EH8+tdzsXdvrEfO\n29LS4tPpTUHOnYjSAFwGYOs0S64E8NrY43cBXEAkZSvC9Pz733HYti0LF1zQhuuvP73eVyqKivqw\nYUMVfvghBtu2Zbp8PMYY2tqkKSXjeAfDw8OorRUeee/YkYGysjRcf30DbrutDv/4h/3ad0+hUo3g\nmWcOIjNzAI88MgcHD7o/RWMT5PNVhEbu/wvgVwCmkzlMBdAIAIwxM4BeAF6nYtXaGoLf/74Aubn9\n+OUvKyXthJuKSy5pxeWXt2DHDg3++c94l4+n0+l8OvLguMbJkycFN+R88UUitm7Nxs9+1oaf/7wW\nRKOlkevWrRvPsdv0Z9atWyeJSqRtjyolZRCPPDIHNTURbj+nuzrhPYFd505ElwPQM8b2z7Rsip+d\n5lWIaC0R7SOifZ7WT7ZYCL/7XSGsVsJjjx1DSIhn5Hgd5Re/qEZBQR/++McCtLW5VgJmNBoF35Jz\n/IvOzk7o9XpBaw8ejMZTTxVgwYJu/PKXJ04JeqaqfZdS/lepNOOppw4jPNyMBx6YB53OvUJ5NkE+\nX0RI5H42gGVEVA9gJ4DziejNSWuaAKQDABEFA1ABOO1+hjG2hTFWzBgrTkhIcMlwIUwcPvDGGxoc\nPTqMc875b6Smeu8fS6FgeOSRClito5u+Vhc/g3x9U4jjOBaLRfBM1La2EDz2WBFSUwfx+OPHoFB4\n/51eYuIwnnzyMEZGZHjggXno73ev/qGvRu92nTtj7CHGWBpjLBPACgBfMcZWT1r2EYCbxx4vH1sj\n6avE1oCxYcMGfPutBa+/HoaIiHPw+ee/8/rpMqmpQ7j33hr88EMM3n03zaVj6fV6v9HK4Aijvr5e\nULRpMsnw2GNFsFgITzxxFFFRU8tYeyNZWUY88cRRtLSE4oknZsOdL/HW1lZYXY2yJMDpOnciepyI\nlo19+zKAOCKqAXA/gAfFMM4VSkpKxjd/Hn30BhDNxcBAldsbMMTikktacc457di6Ndul3KLZbEZH\nh/QVQRzP0N/fP6Uw2FRs2pSLykolHnzwONLTvUdXXSjz5/di3boa7N0bh5decp+So8lk8sn3kEPO\nnTFWzhi7fOzxo4yxj8YeDzHGrmWM5TLGFjPGpC2OxX+GD4SExMFqbYfV2u7RBgxXIQL+67+qEBU1\ngqeeKnCpPJKnZgIDxhgqKysFrf3ssyT89a8pWLlSi5/+1HeVRJcta8GyZc14++0MfPFFktvO44up\nGb/uUD1xIgrDw757iSrVCNatq0Z1dRR27XI+PdPd3e2zm0Ic4TQ3NwvaQG9uDsWmTbMwf34P1qyp\n84Bl7uXee2swb14PnnkmD3V17qmg6enpgdHo/nGZYuK7ns8Oen0vHn10HYB2qFTSNWC4yrnnduDs\nszvwyiuZaG52vjKA17z7N0NDQ6irs++ozebRqrGgIIaHHz6OoCAPGOdm5HKGxx6rQHi4Bb/9bSEG\nB93j1nwtevdb5/7ss4dgNh9HUlI2XnlF2gYMVyAC1q+vglzO8Oyz+XB2m5qLifkvjDFUV1cL2jh/\n4w0Njh9X4v77q9w609fTxMaa8OtfV6ChIRzPPZfnlnO0trb6VHGCXzr3lpZQHDjwCPLyfo8///lp\nr2jAcIWEBBPuuKMWBw7EOJ1XHBwc5DXvfkpHR4egCVxHjyrx5psaXHxxK5Ys8WyfiSdYtKgHN96o\nxeefq/HZZ+IP2zCbzYJ7B7wBv3Tuzz+fi+BgK/7f/zvPqxowXOGKK1pQUNCHF1/MhtHo3L0031j1\nP8xms6CadpNJhqeeKkBCwjB+8QthNfC+yE031WPhwm4899ws0VRWJ+JLqRm/c+7/+lcc/v3veNxy\nSz0SEvxHOEsmG+1e7eoKweuva5w6hl6v98l6Xc701NbWChKIe+01DRobw/H//X+ViIjwndSCowQF\nAQ89dALBwVb84Q8Fote/GwwGn7kD9ivnPjwsw/PP50KjGcDVV/u+2P5kCgsNWLpUh/feS0NDg+NR\nCa959y96e3sFRZLV1ZHYuTMDS5fqPDJpTGoSEoaxfn01jh1TYefODNGP7yuDPPzKue/alQadLgz3\n3VeN4GD/3Dy8445ahIRYsXlzrlObq7xqxj+wWq2CdNrNZsJTT+UjOtqEu+8+6QHLvIMLLtCjpESP\nV1/NRE2NuFOc9Ho9RkZGRD2mO/Ab597VJceOHRk455x2LFgg3fANdxMbO4JbbqnH3r1x+O47x4U3\nuc67f9DU1ISBgQG76955Jx01NVG4775qn5IXcBUi4L77qqBUjuD3vy+AySSeq7NarT6xf+U3zv21\n1zJhMskknxzjCa66qhnp6Ua8+GK2U52rvrTjzzmdwcFB1NfX212n04Xitdc0OOecdpxzTuCl41Qq\nM375y0rU1UXizTfFTc80Nzd7fWmxXzj3hoZwfPxxCpYta/FJjQxHCQ5muOOOWmi1EQ6XfJWVlaGi\nomL8e71ej82bN4ttIsdN2GrahWyMb96cC5lstIMzUDnrrC5cdFErduzIQG2teN2rvjDIwy+c+4sv\nZiMszIKbbtJKbYrH+OlPO1BU1ItXXskU3JFnU8pcu3Yt6urqoNfrsWTJEtx7773cwfsIHR0dgpzK\nnj2x+PbbeNx4o9avmpWc4e67axAVZcb//E++qNUzTU1N4h3MDfi8cz94UIV//SseK1c2IDra+zc5\nxIIIuPPOk+jsDMG776YLes5Epczi4mLMmTMHFRUVKCwsxLXXXutmizmuIrymnfCnP+UiPd2Ia691\nbZRkZGQkEhMTkZaWhvT0dCQnJyMmJgZyuTiD3D2BSmXGvffW4MQJJcrKUkU7bnd3t1frzbhX5d7N\nMAZs2ZKDxMQhXHONd3+KuoM5c/pwzjnteOutdFx+eQtiYmb+cLN16a5Zs2Y8+ktISMDu3buRmJjo\nCZM5LlBfXy9oM/ydd9LR3ByOp546BLnc8bywUqlESkoK4uLipnXijDEYjUZ0dHSgtbUVg4PenQ49\n/3w9vvgiCS+/nI2f/rQTarU4QnpNTU3Iy3OP3IGr+HTk/t13cTh+XImbbtJ67dg8d3P77bUYHg7C\n669nSm0Kx4309/cLSgO0tYXgzTdHN1F/9CPHatojIiIwb948LFy4EGq1esbonIgQEREBjUaDxYsX\nY+7cuVAqlQ6dz5MQARs2VIGIYePGWU5rNE2mtbXVa8sifda5W63Atm2ZSEkZxMUXe39ZkrvIyBjE\npZfq8PHHyWhtnXnmand3NzZs2ICenh5ER0cjNjYW7e3tWLJkCa+g8WIYY4Jq2gFgy5bRoRX33OPY\nJqpGo8GiRYsQGxsLcnByPBEhLi4OCxcuRFFREUJCXJv96y6SkoZx22112Ls3Dl9/7foAemC0LFKn\n04lyLLHxWef+zTfxqKmJws031/ttw5JQbrxRCyLgzTdnliUoLy+HVquFRqPBtm2jSpmFhYWoqKjA\nrl27PGQtx1FaW1vR19dnd11FRRS++ioJ113XiKQkYZuoQUFBmDdvHrKysiCTueYOiAgJCQlYvHgx\nUlPFy22LyVVXtSAnpx8vvJArmjRwU1OTV8p62L06Igolor1EdIiIjhHRb6dYcwsRtRPRwbGv291j\n7igWC/Dqq1nIyBjABRfwjsvExGFcfnkLPv00eUbN99LSUqxbt258GpVKpcLOnTvx/PPP45577vGg\nxRyhjIyMoLZ2tHdj4sB3YPROzDYPmDHgz3/ORUyMCStWCNtEVSgUWLhwIWJjY0W1OSgoCLNmzcK8\nefO8buM1KIhh/foq6PWhdoMhoZhMJq/s/Bby0TUM4HzG2HwACwAsJaKzplj3NmNswdjXVlGtnMQn\nn0Shvj4Ct9xS7xfDBsRg1aoGBAdb8cYbmTOuKy0tPUUp02KxcMfuxdTV1WFkZOSUge/d3d3jKbZN\nmzahrKwM33wTj6NHVbj11jqEh9uv95PL5ViwYAEiI8VtzZ9IbGwsiouLvS4XP3duHy6+uBXvvJPu\nlEbTVDQ2NnpdU5Nd585G6R/7Vj72JdlVmM3A88/HITu7H+ed53+a1M4SF2fClVe24Isvkhx6wfb0\n9PARfF6KwWAYFwabWMa6Zs0arFmzZjzFdvbZS7BlSzY0mgFceqn9/aegoCDMnz8f4eHh7r4EhISE\nYMGCBVCrxddXd4W1a08iJMSKP/1JnM1Vo9EoSFPfkwhKOhFREBEdBKAH8AVjbM8Uy64hosNE9C4R\nCSu8doLt24H6egVuvbUOLqYI7RISEoLk5GTk5eVh/vz5WLRoEYqLi7FgwQLMnj0bGRkZiI6OdngD\nyl3ccEMDFAorXnst06HneeMtZaBj60S1YStjtY2LtG2Kb9y4EV9/XYTm5nDcdddJBAXZ91RFRUVu\njdgnI5PJkJ+fj8zMTI+d0x6xsSO47bY67NsXK9rmqlar9aroXVCdO2PMAmABEUUDKCOiOYyxoxOW\n/AXAW4yxYSK6E8BrAM6ffBwiWgtgLQBkZDin9XD11UBDQxvOPtt9n5JxcXFIT0+HSqUS5LjNZjO6\nurrQ1tYm6ad3TMwISkubsXNnOlav1iIrS1iDRWtrKzIyMrzmQ4oz+oErZBO1vz8Yb7yRiUWLurB4\nsf3O1ezsbNFz7EIgImRmZkKhUAiu/HE3y5a14K9/TcYLL+TizDO7EBrq2qaowWBAV1cX4uIcF/Rz\nBw7FvoyxHgDlAJZO+nknY8y2Pf8SgEXTPH8LY6yYMVackJDghLlAVBRw/fW9cIcfCgsLw4IFCzB3\n7lyHIvLg4GAkJiZi7ty5OPPMM5Gamupy5YGzXH99A8LCLA5F73wEn3dhNpvHN1FtTC5jDQsLQ09P\nD+6997/Q19eFlSv34YMPymY8ri1okZKUlBTMnj1bUhtsBAUxrFtXDb0+FO+8I87/S319vddE70Kq\nZRLGInYQURiAnwE4MWlN8oRvlwE4LqaRniApKQnFxcWIjo526ThhYWGYNWsWzjzzTKSkpIhknXBU\nKjNKS5vx9dcJ0GqF51R9QcI0UNBqtad1ok4sY12+fDkGBwcRHCxHX181cnM3YdOmteObq1OhUCiQ\nn5/vFXdnSUlJKCwslNoMAMC8eb0499x2vPVWBjo6FC4fz2AweE3uXUh4mQxgNxEdBvA9RnPuHxPR\n40S0bGzNurEyyUMA1gG4xT3muofMzEwUFBQgSMTSm5CQEOTl5UlSLbB8eRNCQqzYsUN46ouP4PMO\njEbjlJ2oE8tYL730Umg0GpjNIwAiodf/edzxl5SUTHnc/Px8KBSuOy+xSExM9JoI/uc/PwmLhfDy\ny1miHK+2ttYronch1TKHGWMLGWPzGGNzGGOPj/38UcbYR2OPH2KMFTHG5jPGljDGTsx8VO8hKysL\nmZmZbovTEyeuAAAgAElEQVRoIiMjsXDhQuTm5nosVRMdPYLLL2/Bl18moaVl+rr3iZjNZq+JOAKZ\nmpqaaR2DrYw1JiYGDz30AoAEAP3o6+se31ydWOZqQ61We00eeCJJSUnIz8+X2gykpIxqU332WTIq\nK13faDYajV7RteqzHapikJaWBo1GnEaGmSAipKWlYdGiRR6rUrj++kYEBTGHZkjy1Iy0dHZ2CtYI\nf/ttYTliuVyOnJwcV8xyK8nJycjOzpbaDKxerUV0tAkvvODc+MrJ1NXVwWyWdvJVwDr3mJgYj7/o\nIyIisHDhQiQnJ9tf7CLx8SZccokOn32mRnu7MK0PPoJPOqxWK2pqhOnBHDkyhN27VwBoR3R09Hh5\npK3BaSI5OTle1yU6mYyMDKSlpUlqQ0SEBbfeWofDh6NFKY0cGRkRNC3LnQSkc1coFJg9e7Ykm0tB\nQUHIz89HXl6e289/ww2NsFhIcJQH8Jp3qWhubhYsm7tp00EAFUhPzxzXCLI1OJWXl4+vU6lUSEpK\nco/BIpOTkwNnK+jE4rLLWpGd3Y8XX8wRZeZqU1OTpFVoAencCwoKJN9cSklJwfz58xEc7D5JfbV6\nCBde2IaPP05Gd7ew6I2nZjyPyWQSHOVVVUWipuZBLF78Wzz33LPjOfiNGzdi3bp1KC0tHV+bm5vr\nFdUxQiAiFBQUQKVSSWZDUBDD3XefhE4XhvffF0f47MSJE5IVKgScc1er1ZI0cUxFdHQ0Fi5c6FaJ\n1JUrtTCZZHj3XWG3vQMDA7zm3cPU19fDInD+27ZtWVAqR/Cb3yw5ZfM0JibmFMeuVqsRFRUluq3u\nJCgoCEVFRQgNFVYE4A4WLerGmWd2Yvv2DBgMrgdeAwMDkqVnAsq5BwcHe93mUkREBM444wy36Xxk\nZAyipKQdH3yQKvjFyqN3z9Hf3z+uH2OPI0dU2LNnJ6666gdERo5+GExUhrQhk8mQlSVOWZ+nUSgU\nmDt3rqhlyY5yxx21GBgIdqiUeCYaGhokGaYdUM49KyvLKzeXbOJK7qqkWbVKC6MxGB9+KKypqq2t\njde8ewDGGE6ePCl4/dNPfwHgXpSXXz+lMqSNtLQ0rx2YIYSIiAgUFRVJdv6cnAFceGEb3nsvDXq9\nOP+PFRUVHp+3GjDOPTw83CNVKs6iUCgwf/58tzj4nJwBnHlmJ95/Pw3Dw/b/5Lzm3TN0dnaeVt0y\nHceOKdHQcAtiYnLR0FB/mjKkrXkpODhYcokBMYiNjZX0LvvWW+sAAK+8kinK8cxmMw4fPozhYWFD\nVMQgYJy7GJNm3I1cLsf8+fMREREh+rFXrGhEd7cCn38uTHqVp2bci9VqdShqf/11DVSqaGze/PSU\nypC2/Ht6erpX3p06Q1pammTVPmr1MEpLm/G3v6lRVyfO+3FoaAgHDx70mMS2d3s7kYiMjER8vDiy\nnu7G5uDDwsQZImBj/vweFBT04e230yFk766zs9OjUUag0dLSIrj08fjxKOzdG4frr2+cUblQLpd7\n7Xg7ZyAi5OXlSbYxvHKlFmFhFmzdKt7+xeDgII4f94z0VkA4d41G4zMlYcB/UjRi5k2JRvXeW1rC\n8PXXwuqJec27e3C0weWNNzRQKkdw3nnHTlGGnNy8lJ6e7tbSWimwVdBIcTeiUpmxcmUD/vWveBw+\nLF6Jpqf2s/zeuYeHh/tM1D6R0NBQzJs3T9Q369lndyAtzYidOzMEtVjrdDqvEEDyN7RareDW9Kqq\nSPz73/G49tpG7Nnz99MGnNual7755htJVEg9QWhoqGQqktdc04T4+GG8+GK2KLIEnsTvnXtaWppP\nRe0TiYiIwJw5c0SzPyhoVHOmqioKBw7YlzYeHBwUNDCCIxyj0Yjm5mbB619/PRORkaNDWCYPOJ/Y\nvLR+/Xq/i9onEhMTI4kGTUiIFbfcUoeKChW++ca3gkS/du5yudxn2q+nIzo6GgUFBaId76KL2hAb\nOyxYUMwb1O38CUfkYGtqIvHtt/FYvrwJERGjGyWTB5zHxMRg+fLlfpVrn4709HRJ7sKXLm2DRjOA\nbduyBO1XeQt+7dyTk5MlbYYQi6SkJNHmTyoUVixf3oR9+2JRXW2/7FKv10uubucv9PT0oKOjQ/D6\nN97QICLCjGuumTnST01N9euo3YZNokDsYgN7BAUx3HprHbTaCHz5pe8Ei37t3P0pB6nRaEQTVrri\nCh3Cw83YudN+PbTVaoVerxflvIEMY0yw6iMA1NVF4OuvE3D11U2IjJz+w1Umk0muqOhJgoODUVRU\n5PGy5nPP7cCsWQa89lomRkZ8I83rt849Li5OUo0KsbFFLWI0OUVGmnHFFS0oL08UNMyDp2Zcp62t\nDf39/YLXv/GGBuHhZixffvpUpomo1WrJRfA8TWRkJPLy8jx6TiLgttvqoNOF4dNPvbcZciJ+69y9\nuRvVWWxlYWLcgi9f3gSZjGHXLvvRu8Fg4GJiLmCxWE4beD0T9fXhKC9PQGlpM5TKmVNi/tCN6gxq\ntdrj7/HFi7swZ04v3nhDI6jTW2qEDMgOJaK9RHRobE7qb6dYE0JEbxNRDRHtIaJMdxgrFIVC4TXK\nj2ITFhYmSllYfLwJF17Yhk8/VaOnx34NMY/enaexsdGhIShvvqlBSIgV1147c9SekJDg8fyzN5Gb\nm+uxyWaALXqvRUdHiGCdJikR8vEzDOB8xth8AAsALCWisyatuQ1AN2MsF8BGAE+Ka6ZjJCUleb3U\ngCvExsaKssG6YkUjTCaZIO3qtrY2wbK0nP8wPDyMhoYGwesbGsKwe3cirrqqGSrVyIxrAzVqt2G7\nk/Vk0cSCBb0oLu7Cjh0ZMBq9u1hDyIBsxhizJQvlY1+Ta7muBPDa2ON3AVxAEhaXq9XC9FN8GY1G\nM+UwZEfIyDDi7LM78MEHqRgcnPmFarFYeMeqE9TV1TnUkbh9uwZyuRXXXdc44zqVSgWlUumqeT5P\nWFiYqKXCQlizpg69vQq89553b2QLCm+JKIiIDgLQA/iCMbZn0pJUAI0AwBgzA+gFcNq4dSJaS0T7\niGhfe3u7a5ZPQ2RkpFuEt7wNIsLs2bNd3kxbsaIRBoMcf/2r/fwlT804hsFgcEiArbk5DF9+mYRl\ny1oQEzNz1B5IFTL2SEhI8Oj/x+zZBpx9dgfefjsdfX3eW4IqyLkzxiyMsQUA0gAsJqI5k5ZMFaWf\n1qnBGNvCGCtmjBW7a16irzctOYJCoXA5/15U1Id583qwa1cazOaZb7b4xqpwHC19BIA338xAcLAV\nK1bMHLWHhob6pKSGO8nOzvboncyaNXUwGoMcmk/saRxKTDPGegCUA1g66VdNANIBgIiCAagAeH70\nCCD5kF1PEx0dDY1G49IxVqxogF4fiq++SrS71pHW+UCmo6MDvb29gtfrdKH429/UuOIKHWJjZ958\n9WVJDXchk8lQWFjosWau7OwBLFmix/vvp6GryztLUYVUyyQQUfTY4zAAPwNwYtKyjwDcPPZ4OYCv\nmASKUyqVyq9q24Wi0WhcilrOOqsLWVn92Lkz3a44kl6vx8jIzCmDQMdRrXYA2L49A0FBDCtWzLz5\nGhQUFBB7Ss7gaYGxW2+th8kkw/bt4ozjExshkXsygN1EdBjA9xjNuX9MRI8T0bKxNS8DiCOiGgD3\nA3jQPebOTKBF7TZkMhlmz57tdNUA0Wjuva4uEt99N3MJqdVq5YM87NDU1OTQQIbW1hB8/rkal12m\nQ3z8zFG7Wq0OCKkBZ4mNjXX5TlYoaWmDWLq0FX/5Swra2rxvrKGQapnDjLGFjLF5jLE5jLHHx37+\nKGPso7HHQ4yxaxljuYyxxYwx4R0bIhKozh0YrRqYNWuW088//3w9kpKGBAmKNTc3cyngaTCZTNBq\ntQ495623MkAErFxpv2SSb6TaJzMzE9HR9lVPxeCmm+oBjFY5eRt+UwyuVCp9eiiwGCQlJTn9ARcc\nzLB8eSMOH47GsWMzp3iGhoYkmebuC9TV1TnUD6DXh+CTT5JxySU6JCTMPPkqLi4uoJuWhEJEKCws\n9IgsQ1LSMC67TIdPPlFDp/OulLDfOPdAjtpt2MaSOfuivuyyVkRFjeCtt4RF75xTMRgMDpeL2vY5\nbrjBftQeCLK+YqFQKFBUVOSRjeeVK7WQyUb1gLwJv3HuvDRsFLlcjvz8fKeeGxZmwVVXNePbb+Oh\n1YbPuLarqwtGo9Gp8/gjzpQ+dnQo8PHHKbj44jao1TNH7eHh4S43rQUaKpXKIwM+EhJMWLasBZ9/\nrkZzs/fcWfmFc4+IiOC3qxOIi4tzuqLi6qubERJiEVS/29Q0s/ZJINHe3u5Q6SMA7NyZAYuFsGqV\n/Rx9amoqL390grS0NI/c1a9c2QC53IrXX/ee6N0vnHtc3GnNsAFPbm6uU3sQ0dEjuOSSVnzxRRLa\n22dO77S2tvKySIxKMzha+tjZqcBf/pKMiy5qRUrKzJU1QUFBAdWcJyZEhPz8fLcHf7GxJlx1VTO+\n/DIJDQ3eEWj6hXPnKZnTCQ4Odjo9c911jbBaya52htVqRUtLi1Pn8CcaGhowPDxzWmUyb7+dDrNZ\nhtWr7efaefmjawQHB2POnDluFxNcsaIRCoUVr7+e6dbzCMXnnbtcLkdUVJTUZnglsbGxTqVnkpOH\nUFKix1/+koL+/pmdSnNzs0PCWP7G4OCgQ6qPANDVJcdHH6XgZz9rQ2rqoN31fCPVdSIiItwuMBYd\nPYKrr27CV18lor5+5j0rT+Dzzj0uLo7nImcgJyfHqeqZG25ohNEYbFe32mQyBbRaZE1NjcM1/++8\nk46RERlWr7afa4+JiUF4uPSOwh9ITEx0e5/Addc1IjTUgldfzXTreYTg887dX4dyiIVcLndqJFlu\nbj+Ki7vw3ntpMJlmfpk0NjYGZFNTZ2cnOjs7HXpOT48cH36YivPP1yM9nUftniY7Oxsqlcptx1ep\nzJg79/f4xz+AkydH1Wm7u7tRVlbmtnNOh087dyLizl0A8fHxTlUM3HBDA7q7Ffjss5k384xGIzo6\nOpw1zyexWCyorq52+Hm7dqVheFhY1B4aGsqLBURGJpOhqKjIbQ1OZWVl2Lv3MRAtwUsvRaC7uxsb\nNmzApk2bPO7gfdq5K5VKvtEkkFmzZjn8f7VwYQ/y8/vwzjvpsNd02dDQEFDRe0NDg0P6MQDQ2xuM\nsrJULFmih0Zjv0cgJSWFpxzdgEKhwJw5c9zyf1tSUgKNRgPGKrBnTwluvvk2aLVaaDQalJSUiH6+\nmfBp586jduEoFArk5OQ49Byi0dx7c3M4/vnPmSN/g8GA7u5uV0z0GYxGo8ObqADw7rvpGBoKEhS1\ny2Qyvxzy7i0olUqn0pX2iImJwcaNG6FSRQNoh8HQjejoaGzcuNHjTWjcuQcQarXaYUGln/60HWlp\nRrz1ln05YEcFs3wRxhiqqqocvkvp6wvG+++n4txz25GVZT9qT0xMhFxuf3A5x3mSk5ORkuKeQdcT\nbwrsDcFxFz7r3OVyuUcnn/sDNu0ZR25Hg4JGKwAqK5U4cGDmD4be3l6/j97b2trQ09Pj8PPeey8N\nRmMwbrpJ2Aegu5wO51Ryc3NF3WC15dh7enqgUkWDKAH9/aM/8/R7w2ede0xMDM9HOkF4eLjDetcX\nX9yGuLhhvPmm/efV19f7be7dZDI5rB8DAP39wXjvvTScc047srMH7K6Pioriw689hNgbrOXl5eM5\n9lde2YbVqz8GUAitVovy8nJRziEUn3XuPCXjPBkZGQ61YysUo3M9Dx6MweHDM0c5/hy919TUwGw2\nO/y8XbvSMDAgPGrn5Y+eRcwN1tLSUqxbt248x37DDUNQKj9HevqTKC0tFcFa4fisc+cKec4jk8kc\n3ky6/PIWxMSYBMma1tXV+V303tHRAb1e7/Dzenvl2LUrDeedp0dubr/d9XK5nMtXS4BSqXRarmMy\npaWl4/4pLMyKVauG0Nj4Kxw65L76+qkQMkM1nYh2E9FxIjpGROunWFNCRL1EdHDs61H3mDtKRERE\nwA/mcJWYmBiHxKhCQ6247rpG7NsXi4qKmVMGBoMB7e3trproNYyMjKCqqsqp57711miFzC231Ata\nn5yc7PS4RI5rqNVqt9w1LVvWgtjYYY93rQqJ3M0A/osxNhvAWQDuIaKpptB+wxhbMPb1uKhWToJH\n7eKQk5PjUO37lVe2QKkcESRrWldX5zeaMydPnoTJNPNs06no7FTggw9S8bOftSEzU5j2Pd9IlZac\nnBzRR/SFhlqxcmUDDh6MwQ8/eGb8HyBshqqOMXZg7LEBwHEAkiYFuXa7OCgUCmRlZQleHxZmwbXX\nNmLPnjhUVs5cqTQ4OOgXipEdHR1ODwTfvj0DIyMy3HxzvaD18fHxCA31rlFtgYZMJkNhYaHomYEr\nrtAhPn4Yr7ySabekWCwcyrkTUSaAhQD2TPHrHxPRISL6lIiKRLBtJjvcefiAIiUlxSFVzdLSZkRG\njuCNNzLtrq2vr/dpvXeTyYTKykqnntvWFoKPP07BJZfokJoqrJOVR+3egW2DVUyJYIXCitWrtThy\nJBp79nhGxVaw9UQUCeA9APcxxvom/foAAA1jbD6APwH4YJpjrCWifUS0z59ysr6MrfZdKBERFixf\n3oRvv41HTU3EjGvNZjPq6+tdtFAaGGOorKx0+sPJtvF8443CKmT4GD3vIioqSvQO1ksu0SEpaQgv\nvpjikehdkHMnIjlGHft2xtj7k3/PGOtjjPWPPf4EgJyITpugwRjbwhgrZowV84oA7yEqKsqhqPHq\nq5sREWEWVPfe3NyM/n77VSLeRktLi8OKjzaam8Pw6afJuPzyFiQlCRviwcfoeR9qtVrUuymFgo1F\n75H47DPRDjstQqplCMDLAI4zxp6dZo16bB2IaPHYcZ17Z3AkISsrS3C7e1SUGVdf3YR//CPRbvQO\nwKl2fSnp7+93qlnJxmuvaSCXWwVNWQL4GD1vJjc3V9SGsqVLW5GePoQDB0Q75LQIidzPBnAjgPMn\nlDpeSkR3EtGdY2uWAzhKRIcAbAKwgvnSu5kDuVzu0KT4a69tQmTkCLZts78h29fX5zObq2azGceO\nHXP6w6imJhJffpmE0tJmxMYKq7DhY/S8F9sGq1g6P8HBDG+9VYFf/1qUw818LnsLGGP/BDDj/SJj\n7HkAz4tlFEca1Go1dDod+vomb6mcTlSUGStWNGLr1mwcO6ZEUdHMz6mtrUVcXJxXV4PY8uyDg/aH\naEzHiy9mIyrKjFWrhKtG8o5U7yY0NBSzZ8/G4cOHRTqeZ+Jen+1Q5YgPEWHWrFmC1199dRNiYkzY\nujXL7gaRxWJBZWWlV6dnmpqaXGq++v77GOzbF4vVq7WIjBQmU8DH6PkGsbGxyMjIkNoMh+DOnXMK\njmyujrZWa3HwYAwOHLBf6dHd3Y3m5mZXTXQLXV1dOHnypNPPt1qBF1/MQXLyIK68Uvg1unumJ0c8\nMjMz3TqiT2y4c+echiObq1dc0YLExCFB0Tsw2u3pbdUzAwMDOHbsmEvH+PLLJJw8GYnbb6+DQiHs\n7iQsLIwL4PkQMpkMs2fP9hl5CO7cOafhyOaqQsFw8831OHFCiW+/tT/vkzGGY8eOOaWu6A5MJhOO\nHDkCi705gjMeQ4aXX85Cfn4fSkqEi4vx8kffIzQ01C0TnNwBd+6cKVGr1YJLwC6+uA3p6Ua8/HK2\n3VmrwKg0wYkTJyTPv5vNZhw+fNjhWaiTef/9VOj1ofj5z2shtKkxKCgIarXapfNypCEpKcknlDu5\nc+dMyVSbq2VlZadotXd3d6OsrAxBQQy33VaL+voIfPqpsLmfHR0dqKurE9VmR7BYLDhy5IjLKaKu\nLjnefFODs87qxMKFwic08fJH3yYvL8/rxyBy586ZlqioqPEhzWVlZdi0adP4uDDbOLFNmzahrKwM\n557bgblze7BtWxaMRmE5yYaGBknq322Ovbe31+Vjbd2ajeFhGe66y7GmJ17+6NvI5XKvT89w586Z\nkezsbAQHB6OkpAQajQZarRZr1qzBmjVrxseJlZSUgAi4++6T6O5WYMcO4SVjVVVVaGtrc+MVnMrI\nyAgOHTrk1BzUyZw4EYXPPlPjmmuakJEhvDY+Li6Olz/6AQkJCYiPP01lxWvgzp0zI3K5HFlZWYiJ\nicHGjRsRHR2Nnp4e9PT0IDo6enycGAAUFBhw4YWteOeddLS2CpdMPX78uEcc/NDQEH744QdBTVr2\nsFqBP/0pFzExJsHiYDZ4+aP/kJubK6p6pJh4p1UcryIlJQWRkTPrt9u4/fY6EDFs3SpcygAYdfCN\njY1u22Tt6enB/v37YTQKG5phjy++SEJFhQpr19YiIkJ4pU1ERITowyA40hEaGorMzEypzZgS7tw5\ndiEiREdHY8OGDeMRuy2Ct+XgbSQmDuP66xvx978n4ehRxwSXTp48icrKSpfKEifDGINWq8XBgwdF\n05YfGAjCli3ZKCzsxYUXOnbHkZ6ezssf/Yy0tDSvHCDEnTtHEJ9//vl4jn3btm3Ytm3beA6+vLz8\nlLU33NCIhIQh/O//5sFiccyRtba2Yv/+/aKkTgwGAw4cOCB6Vc7Wrdno7lZg3boawaWPwOgQiMTE\nRFFt4UiPTCZDbm6u1GacBq/F4gjinnvugdlsRlZW1nj9+8aNG1FeXo7S0tJT1oaFWXDPPTX47/+e\ng7KyVCxf3uTQuYxGIw4cOAC1Wo3MzEyHxcYGBweh1WqdHo83E8eOKfHhhym4+upm5OcbHHpuamqq\n1+ZnOa4RGxs7fjfrLZBUjSTFxcVs3759kpyb4zxNTU2CtM4ZAx56aC4OH1bhtdf2IiHB8QHTwGhK\nKD4+HklJSYiJiZm29XtkZATd3d1oa2tzesiGPUZGCD//+SIMDATjlVe+R3i48PSRTCbDj3/8Y6+v\njeY4T19fHw4IEGqPiorCokWLnD4PEe1njBXbW8cjd45DpKSkQKfTYWBgYMZ1RMAvflGNNWt+hBde\nyMVjj1U4dT7GGNrb29He3g4iQmRkJEJDQ8cbgEZGRmA0GkXbKJ2Jt99OR11dJH73uyMOOXYASE5O\n5o7dz1EqlYiPj0dHR4fUpgDgOXeOg8hkMsGywKmpQ1i1qgHl5YnYs8d1gSzGGAwGA9rb26HT6aDT\n6dDR0eERx97QEIbXX8/Eeefp8ZOfOH5nkJ6e7garON5GVpb94TWegjt3jsNER0cLHgu3YkUDNJoB\nPP10Pvr7ffNG0WIh/OEPsxEWZsEvfuH4+L3ExESvHlLCEY+IiAiv0Z0RMkM1nYh2E9FxIjpGROun\nWENEtImIaojoMBGd4R5zOd5Cdna2IOlThYLhwQdPoKtLgc2bczxgmfhs356BEyeUuO++KsTFOb53\n4GtDHjiuodHYHxzvCYRE7mYA/8UYmw3gLAD3EFHhpDWXAJg19rUWwJ9FtZLjdYSEhAi+BS0oMGDF\nigZ89lkyvvvOt/TLKysj8frrGlxwQRuWLHF8SlNsbKzgBjCOfxAZGekVOv12nTtjTMcYOzD22ADg\nOIDJqkdXAnidjfIdgGgiEiYPyPFZUlJSEBERIWjtzTfXIzNzAM88k4/eXt9IzwwOyvCHP8xGTMwI\n1q+vduoYPGoPTLzh7+5Qzp2IMgEsBLBn0q9SATRO+L4Jp38AcPwMmUwmWBlPoWB46KHj6OmR48kn\nCwRNbZISxoCNG/PQ0BCOBx88jqgox4eLqFQqLjUQoKhUKkRFRUlqg2DnTkSRAN4DcB9jbHL74FRt\niKe9fYloLRHtI6J9rgwi5ngPKpVK8NCJvLx+3HnnSfz73/F4913XxLOm05YXi08+UeOLL9S4+eZ6\nLFrkXGOKN0RvHGkgIskF4gQ5dyKSY9Sxb2eMvT/FkiYAE2u90gCcJtTNGNvCGCtmjBV7y44yx3Vs\nssBCuPrqZvz0p+3YsiUbx487FtnYHLpNW37dunXYsWPHadryrlJTE4FNm2ahuLgLq1c7pvhoIyoq\nyivyrhzpSEhIkLS3QUi1DAF4GcBxxtiz0yz7CMBNY1UzZwHoZYzpRLST48UoFArk5AirhCECfvWr\nSsTFmfDb3xahq0vYi3/isJAFCxYgLS0NTU1NeOmll7Bq1apTtOVdobtbjkcfnQOl0oyHHz4OZ2ch\nazQaLhAW4MhkMqSkpEh3fgFrzgZwI4Dziejg2NelRHQnEd05tuYTALUAagC8BOBu95jL8VYcmbka\nFWXGb397FD09o47UZLL/Mpw4LOT+++8/RVhscHDwNG15ZzCZZPjNb+agq0uBJ544ipgY51QkIyMj\nERdnf1g4x/+xTTKTAiHVMv9kjBFjbB5jbMHY1yeMsf9jjP3f2BrGGLuHMZbDGJvLGOOiMQEGESEv\nL09wtJqf34+HHjqBY8dUePrpPLsbrJOHhfT19YkaGVutwFNP5ePYMRUeeugECgocEwWbSGZmJo/a\nOQBG9d6l+qDnHaoc0YiMjHRoE+m889qxZk0dvvhCja1bHW/bZoxBqVROqy0v/DjA88/n4u9/T8Lt\nt9fivPOc3+znUTtnMlJF79y5c0TFUYne1au1WLasGTt2aPDmm9NXl9g2TXt6ek4ZjKBSqfDss89O\nqy1vD8aALVuyUVaWhmuvbcTKlQ0OPX8yWVlZPGrnnEJsbKwkG6vcuXNEJSgoyKGp8ETA+vXVuPDC\nVrz8cjbeeWfqyL+8vHx803T79u24/fbbkZ6ejsbGRhw8eBAbN27EunXrTtOWnwnGgK1bs7BzZwaW\nLWvGXXedhCt+WalU8goZzmnIZDLB5cJi4hutghyfIjY2FomJidDr9YLWy2TAAw9UwmSS4c9/zkVv\nr3xsFut/1ticdklJCWJiYrBq1SpceumlpwwLccSxWyyEZ5/NwyefJOOyy1qwfn21S44d4FE7Z3qS\nkpLQ2Nhof6GIcOfOcQu5ubno6uqC2SysszMoiOE3v6mAUmnGjh0adHSE4P77qxASYh1fM9l5x8TE\nOLGVqmgAAAp4SURBVOTQbfT1BeN3v5uNvXvjcOON9bj11nqXHXtsbKxLlToc/yYyMhKRkZHo7+/3\n2Dl5WobjFhQKhWDddxtBQcCGDVW49dY6/O1vatxzzxlobhZ38HBlZRTWri3GgQMx2LChEmvWuO7Y\nAe/S8eZ4J0JlssWCO3eO20hMTHQ4B00E3HSTFn/4w2Ho9SG4445FePfdVFgcG3x0GsPDMrz0Uhbu\nvvsMMAY899wPWLZMnD67pKQkyXVEON6Pp4ejc+fOcRu22nchuu+TOeusLmzZsg9z5/Zi8+ZZuOuu\nRdizJ9ZhwTGLhfD550m49dYfYccODS66qBUvvbQPhYXO17FPRCaT8aidI4iQkBCPCsnxnDvHrYSG\nhiI3NxeVlZUOP1etHsYf/3gEu3cnYMuWHDz44Dzk5hqwdGkrzjuvHfHxUw/OYAxobAzH7t0J+Owz\nNVpbwzBrlgFPP33QaRGw6UhPT+dTljiCSUhIQGtrq0fORUwi7dXi4mK2bx9vZA0EGGM4cuQIurq6\nnD7GyAjhb39T48MPU1BdPZoCSU01IjPTiISEYQQHW2EyydDaGora2kh0dIQAABYt6sJVV7Xg7LM7\nRMmtT0ShUGDx4sWCRdM4HJPJhKNHj+KMM5wfVkdE+xljxfbW8Vclx+0QEfLz8/H9998Lrp6ZjFzO\ncNllOlx2mQ51deH4/vtYHDmiQnNzGA4dUsFqJQQHMyQlDWH+/B4sWNCDH/2oC0lJwyJfzX/Iycnh\njp3jEAqFwmOpGf7K5HiEkJAQzJo1C8ePH3f5WFlZRmRlGXHddU0iWOYcKpXK4xtkHP/AU3LnfEOV\n4zGSkpL8wiE6KpLG4UzEUzN1uXPneJRZs2YhJCREajNcIj09XfDsWA5nMp4KCrhz53gUuVyOwsJC\nqc1wmrCwMGg0GqnN4HDswp07x+OoVCqfrQ3Pz893qm6fw/E03LlzJCEjI8PntFhSU1M92oTC4bgC\nd+4cSSAizJ4922fy72FhYcjOzpbaDA5HMEIGZG8jIj0RHZ3m9yVE1Dthvuqj4pvJ8UcUCgWKioq8\nvurE9kHE0zEcX0JI5P4qgKV21nwzYb7q466bxQkUlEqlQ8M9pCArK0vw8G8Ox1sQMiD7awDO941z\nOHZITk5Genq61GZMSVxcnNfaxuHMhFg59x8T0SEi+pSIiqZbRERriWgfEe1rb3d+CDHH/8jOzkZ8\nfLzUZpxCWFgYCgoKvD5txOFMhRjO/QAADWNsPoA/AfhguoWMsS2MsWLGWLGnWnA5voEtr+0t6Y+g\noCDMmTNHksHGHI4YuOzcGWN9jLH+scefAJATkXeFYByfICgoCHPnzpW8+5OIUFRUJLkdHI4ruOzc\niUhNY/etRLR47Jidrh6X459s3rz5lMHZer0emzdvHv9eLpdj/vz5CA8Pl8I8AEBBQYHDE6Q4HG/D\nriokEb0FoARAPBE1AXgMgBwAGGP/B2A5gLuIyAxgEMAKJpVIPMer2bx5M+6991688MIL2L17NwBg\nyZIlqKioAADcc889AEZLJBcsWIBDhw5hYGDAozbm5eV5fNYlh+MO+LAOjsfQ6/Xjzty259Le3o7C\nwkLs3r37NMXIkZERHD16FL29vR6xLy8vDykpKR45F4fjLEKHdfAOVY7HSExMxO7du5GQkID29na0\nt7cjISFhSscO/CdF4+5IWiaToaioiDt2jl/BnTvHq5HJZCgoKEBubq5bShJDQkKwYMECjw1Q4HA8\nBXfuHI9hS8vYInZbBL9kyZJTNlknQ0RIS0vDGWecIWoFS0JCAoqLi72m/JLDERPu3DkeY9euXaio\nqEBhYSGOHj2Ko0ePorCwEBUVFdi1a5fd50dFRWHRokUuzy4NDQ1FUVERioqKeB07x2/hM1Q5HsNW\nDXPttdeO59h3796NXbt2jf/OHjKZDOnp6UhOTkZLSwtaWlowNDQk6LkRERFIS0tDUlISZDIe13D8\nG14tw/FpGGMwGAzo6uqCwWDA4OAgRkZGAIxuyIaFhUGpVCIuLg4RERFcSoDj8witluGRO8enISIo\nlUqeN+dwJsHvTTkcDscP4c6dw+Fw/BDu3DkcDscP4c6dw+Fw/BDu3DkcDscP4c6dw+Fw/BDu3Dkc\nDscP4c6dw+Fw/BDu3DkcDscPkUx+gIjaAWidfHo8gA4RzfEF+DUHBvyaAwNXrlnDGLOrUS2Zc3cF\nItonRFvBn+DXHBjwaw4MPHHNPC3D4XA4fgh37hwOh+OH+Kpz3yK1ARLArzkw4NccGLj9mn0y587h\ncDicmfHVyJ3D4XA4M+DVzp2IlhJRJRHVENGDU/w+hIjeHvv9HiLK9LyV4iLgmu8nogoiOkxEfyci\njRR2iom9a56wbjkRMSLy+coKIddMRNeN/a2PEdEOT9soNgJe2xlEtJuIfhh7fV8qhZ1iQUTbiEhP\nREen+T0R0aax/4/DRHSGqAYwxrzyC0AQgJMAsgEoABwCUDhpzd0A/m/s8QoAb0tttweueQmA8LHH\ndwXCNY+tiwLwNYDvABRLbbcH/s6zAPwAIGbs+0Sp7fbANW8BcNfY40IA9VLb7eI1nwvgDABHp/n9\npQA+BUAAzgKwR8zze3PkvhhADWOsljFmArATwJWT1lwJ4LWxx+8CuIB8e0im3WtmjO1mjBnHvv0O\nQJqHbRQbIX9nAHgCwFMAhE3D9m6EXPMdADYzxroBgDGm97CNYiPkmhkA27xEFYAWD9onOoyxrwF0\nzbDkSgCvs1G+AxBNRMlind+bnXsqgMYJ3zeN/WzKNYwxM4BeAHEesc49CLnmidyG0U9+X8buNRPR\nQgDpjLGPPWmYGxHyd84DkEdE3xLRd0S01GPWuQch1/zfAFYTUROATwD8wjOmSYaj73eH8OYB2VNF\n4JNLe4Ss8SUEXw8RrQZQDOA8t1rkfma8ZiKSAdgI4BZPGeQBhPydgzGaminB6N3ZN0Q0hzHW42bb\n3IWQa74BwKuMsWeI6McA3hi7Zqv7zZMEt/ovb47cmwCkT/g+Daffpo2vIaJgjN7KzXQb5O0IuWYQ\n0c8A/BrAMsbYsIdscxf2rjkKwBwA5URUj9Hc5Ec+vqkq9LX9IWNshDFWB6ASo87eVxFyzbcBeAcA\nGGP/BhCKUQ0Wf0XQ+91ZvNm5fw9gFhFlEZECoxumH01a8xGAm8ceLwfwFRvbqfBR7F7zWIriRYw6\ndl/PwwJ2rpkx1ssYi2eMZTLGMjG6z7CMMbZPGnNFQchr+wOMbp6D/v/27VA3YSCO4/j33oEnmNue\nAE8ywRNMYTC8wxwPgOMZJubQEzMYFAnJBCETmIag0IhOXMWCGYHuDi7fT1LTNM3/32t/aXrXEDrE\nzzTfSats1zk9b4EeQAjhkRju+6RVpjUDBs2qmS5wqOu6au3suWeU/5ht7gNr4iz7a7NvTHy4IQ7+\nO7ABFsBD7poT9PwB7IBls81y1/zfPZ8c+8mdr5Y5c5wDMAG+gBXwkrvmBD0/AXPiSpol8Jy75iv7\nfQMq4Eh8Sx8CI2D0a4ynzfVYtX1f+4eqJBXolj/LSJIuZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXI\ncJekAhnuklSgH2s1VHaHhI0nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa246ed8dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xstar = torch.linspace(0,1,1000)\n",
    "mu, vsq, nll = gp.forward(xstar=Variable(xstar.unsqueeze(1)))\n",
    "cred_size = vsq.diag()**0.5*2\n",
    "pyplot.plot(xstar.numpy(),mu.data.numpy(),'b')\n",
    "pyplot.fill_between(xstar.numpy(),mu.data.numpy()+cred_size.data.numpy(), mu.data.numpy()-cred_size.data.numpy(),facecolor='0.75')\n",
    "pyplot.plot(X.numpy(), Y.numpy(), 'kx', mew=2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "B2A048522FDA4F3C8666F773BA8F809F",
    "collapsed": true
   },
   "source": [
    "I appreciate your feedback at <tv@lernapparat.de>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cell_id": "513CA294C85944AC9417C7C16DCEB8EC"
   },
   "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}