gp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "test"
    ]
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.linalg import inv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Header 1\n",
    "## Header 2\n",
    "### Header 3\n",
    "#### Header 4\n",
    "\n",
    "Some body text\n",
    "\n",
    "![Image of Yaktocat](https://octodex.github.com/images/yaktocat.png)\n",
    "\n",
    "* List one\n",
    "* List two\n",
    "* List three"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol><li>One</li><li>Two</li></ol>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<ol><li>One</li><li>Two</li></ol>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "code",
     "python"
    ]
   },
   "outputs": [],
   "source": [
    "def compute_kernel(X, Y, kernel_fn):\n",
    "    vectorized_binary_fn = np.vectorize(kernel_fn, signature='(),()->()')\n",
    "    x_indices, y_indices = np.meshgrid(np.arange(X.shape[0]), np.arange(Y.shape[0]))\n",
    "    return vectorized_binary_fn(X[x_indices], Y[y_indices])\n",
    "\n",
    "def sum(kernels):\n",
    "    def kernel(x1, x2):\n",
    "        s = 0\n",
    "        for k in kernels:\n",
    "            s += k(x1, x2)\n",
    "        return s\n",
    "    return kernel\n",
    "\n",
    "def gaussian(b):\n",
    "    def kernel(x1, x2):\n",
    "        return np.exp((-0.5)*(x1-x2)**2/b**2)\n",
    "    return kernel\n",
    "\n",
    "def polynomial(n):\n",
    "    def kernel(x1, x2):\n",
    "        return (1 + 0.01*x1*x2)**n\n",
    "    return kernel\n",
    "\n",
    "def periodic(p):\n",
    "    def kernel(x1, x2):\n",
    "        return np.cos((x1-x2)/p)**2\n",
    "    return kernel\n",
    "\n",
    "def linear(a):\n",
    "    def kernel(x1, x2):\n",
    "        return a*(x1-x2)\n",
    "    return kernel\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(100)\n",
    "# print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x10cc83630>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGtFJREFUeJztnV+M5VZ9xz+/md0lJAglCwWW3bQJUkpBVDQoogGqChFQgSLCQ6gCCK1QqrzQEigSJPQBKlEJJMSfhyrSqilKK9RA06iJUgRCITz0ZWEDqEBCIA0o2SSQUCVQpc3O3nt/fbA91/bY4zO+9rHv+PuRVnvtc3x8xnPt7znHv+9vzN0RQkyLjaE7IISIj258ISaIbnwhJohufCEmiG58ISaIbnwhJohufCEmyEo3vpm92czuN7MHzOz6rjolhOgXaxvAY2abwE+ANwGnge8A73L3e7vrnhCiDw6scOyrgQfc/UEAM7sFuBKovfGff3jTL7rwYFDjTvUDqWpvVrdclt/OHnDlOgussN/T7cI+t0KZbx+7UWgjqVtsZ7t9t511szql9pfbFLYLffBy+xX9z/blji/8PJUX07KGiv+XyyvKrO4X4AF1qspr6u7Wxs4yL+6v/gIldUrHVF+DUnt1x+QPKu/bsR1SN6ta0ancx2d4mi0/U/yFV7DKjX8UeDi3fRr4w3IlM7sWuBbgt48e4NtfvzCo8bkvAFiUfvIFi1wdL+w7mx4zJ9u/ZCute9YplG15cvOeTW/is76c/Zz1TQCe8YNp3c20brp/cahQXjhmcbBQdibbXizrnlkcKP1/sLSdtjVfHrM1T8q20rIz6fbZ+WZavrnsy3yjUDZPt2ezZHuRbi/my+/JYpb+/Nm+bDu9YDZb1s0+W1p3Y57un5fLtw9hY5aWbbdH8djZsm523Mbcq+uWygvtp/vKx5b3J2VpnUVWd1GsO8t9kxaLQh2y9udpZ+ZZ+XzHMVnZcjup47PcDz0vlWXtZsK1fZ5l+/k6J/0uQljlxq96qux4nrr7CeAEwGWvPCd4XrFp6Reu9ADYyC9LZN+e9GY9mPXIs1/U8nSHzAr7sgfAoe02sjYDOpfV3dgi7dxOsm4uStsjYvl1y3cuu96lHyB9APiB+l9hNvpYtrbzom5fjuzkB0r7d23XK+suCr+07HtipbrFvuR/Zdvn3O5T8Wdf5HqV9dsPZK0t0rNu5lov3gj5h1kVxd7PKsuym9s2N3e2n68TeIet8nU8DeTl+xjw6ArtCSEisYrifwe4xMwuBh4Brgbe3UmvctQpP+TUv1H5IXsUBis/NKt/WfmTDhZZK+WHcofrlB/q1b9e+aGsuHXKXzy+3G618ufrhCp/ZTs1yp982ij0u0n5lz1pVv78caHKX2g/v9FA6xvf3Wdm9hfA14FN4B/c/Udt2xNCxGMVxcfdvwp8taO+CCEisdKNH5PykB8qFvzqhvyF4wKH/Mui8CE/1C/47Rw1jnzY3zDkh8YFv51Dfihf1Lohf6Gstt2db292Lvg1Dfkr+pLt3jHkX9YKHfIvW+9nyA+lBb9Z05c1YYRfPSFE36yN4mdsKz/Uv+orKT/Uv+qrU/6krIdXfTvFY5SP32blX+4LVf58a43KD42v+srKX1W3SfmLxzcoP9S+6qtT/uTM1a/6ulB+KC74mUnxhRA1RFV8x5n7oqjaK9AY5FOYr+8e5FNWflCQD9Qrf/Kp+yCfghIHBvkURxRNQT5F5S/WbVB+aAzyKSt/0truQT6rKH+hLDfvb2KEXzUhRN9En+Mv8G21jab8EBzee6gwR1J4b8b+C+/NB4J1H96bX6voI7y3VvnDpvhj/IoJIfpmkFX97bDb3Dv5LtS/m/DeSMaepINF1kr5ITS8t53yQ2h4byxjT+U5OzT2LHuyirFHq/pCiBoir+on3vlMfQte+w7n/StF+cUy9kB4lN9IH8/BUX6tjD0QGuUXz9hT0ZdsdwfGnmXrK6z4a44vhKhDN74QEyT64t7cfXuYnX/dVl7w63TIn2t3NMaefF0Ze7ZpFd4bydiTP74PY09y5hWDfBSyK4SoI3rI7oLF0jyTU8odC36xg3xk7FmZLo09Sc2w8N5Yxp5i3e6NPUlrq4b3SvGFEDVEf5131hc5Jd0ZUjsa5c/1T8aevbFfjT1JWVh4bxtjD4SH99Yqv17nCSHqiL+qnzPpFFfQi4o7uPLn+iJjTzv2n7Fnec7x5u3XHF8IUcMAIbvZJwrhsTvm/Q3KD2tu7AHl7aer8N44xp58nT6MPcmnsPDeWmOP5vhCiDp04wsxQeIO9d3Zcq98DbZjwa9pyJ87Zi0dfcsi5e1ntfDeWI6+Yq/7cPQta7V39GlxTwhRwwABPNmn6tdgocoP623sScoU3ttFeG8sY09V3dHl7dfinhCijgGy7OZVr/41WJPywwiCfFYw9oCCfGC48N42xp6kLCy8d7C8/bLlCiHqiKr4C4wt39hWuLOFB3eNMtYpP4wnvLeVsWe5T8ofP7y3jbGn2G73xp7COVcw9oQwwq+AEKJv4q/qs7H9MMy/0y7P+xuVH0YT3tvG2APK219FvLz94zP2VJ5zj8Ye16q+EKKORsU3swuBfwReRPK4OeHuXzCzw8CXgYuAnwN/5u5P7taWY5z1jUpv4c55/+7KD+OJ8mtn7IHgKL82xh4Ij/Ib6eO/37z9ezf2FI8vt9uFsaeiL9nuYGNPd6v6M+DD7v4y4HLg/Wb2cuB64C53vwS4K90WQqwBjTe+uz/m7t9NP/8PcB9wFLgSuDmtdjPwjr46KYTolj0t7pnZRcClwEnghe7+GCQPBzN7QdPxScju0lxQZSoOHvLDeMJ72xh7Csf1YOzJ15WxZ5tVjD2Fstp22xt78se3NvZ0vbhnZs8B/hX4oLv/Zg/HXWtmp8zs1FP/PQ89TAjRI0GKb2YHSW76L7n7benuX5rZkVTtjwCPVx3r7ieAEwAv/f1z/Bk/WHOS7IDkv2blz1UeOLy3jbGn2F8Ze0IY3NgD48/b35Xim5kBNwH3uftnc0V3AMfTz8eB28NOKYQYmhDFfx3wXuAHZvb9dN/HgE8BXzGza4CHgHc2NeQYW/k5fhWByp+UTSW8V8aePMrbv4uxJ9Ck03jju/t/UP+1uiLoLEKIUTFAyO5mMXaijkblXxZ2Gd6rvP1ZGcWytVB+mEp4b52xJ+TWquunEGKfE9mWu8Ezi0NLpVlB+aGn8F7l7S+yVsoP/Rp7YPR5+2XSEULUET0RxzN+MPeIyilNk/qvYuyB4Cg/pe+eSJRfK2MPhEb5tTH2FHvd3tgTwgh/dUKIvtGNL8QEifyXdEomncIYJ3DBr5WxJ1d5LMaeXLvK2989fRh78q31YeypqrtnY48W94QQdUQO4DGeWRyse9eR0MGrPuXtT1CQT7fGnqRmWHjvUHn7lXNPCFHLwK/zKiulZV0q/7JwOsae5T4p/3oYe4rttjT26C/pCCHqiLyqb5zJz/F3j3JI98cN790vxh5Q3v4qpmDsCWGEvxohRN/En+Mvcqm3Qt4Jd6D80HPe/lEaeyA4vLeNsQfCw3tHKC9tjD2wl/DevRt7iseX2w0z9mhVXwhRi258ISZI9ACeM4uaUzYND1dx9OXqTsbRVziuB0dfvu5UHH0w/rz9GuoLIeqIu7jnuyh+RqjyQy/hvfvF2FPsr4w9IeyXvP0hjPDyCyH6ZoA5fs1f0ikT8kpIxp6Owntl7MmzDuG9dcYevc4TQtQyvjl+mcGUf1m4lsaeXF9k7GnHWob3SvGFEHVE/0s6ZxYNfzuvjn1m7IExhfe2MPaA8vbTtfJDF3n7QxjhZRZC9E30Of4z88BV/TraGHsgPMovlrEnd9zgUX7K2x/M6PP2a44vhKhDN74QEyR6AM/WvMNT9hHeG8nYk3RvHOG9ytu/d8aat18BPEKIWqLn3Ntq+zpvNwYO721j7Em6OZLwXuXtb80Y8/aHMMJLKYTom2DFN7NN4BTwiLu/zcwuBm4BDgPfBd7r7lu7tbHAONPlHL/MOhl7YDzhvcrbvzKjMfb0MMe/Drgvt/1p4HPufgnwJHDNHtoSQgxIkPya2THgT4G/Bf7KzAx4A/DutMrNwCeAG3drxx3OznuY45eJHN7bytgDweG9YzT2gPL2VzEKY08AoZfu88BHcud4HvCUu2enOg0crTrQzK41s1Nmdmr26/8L75kQojcaFd/M3gY87u73mNnrs90VVSsfW+5+AjgBcN7vHvGtGIqf0Sa8N5KxB5S3v5I24b2jVn7oJ7y32tgT+h4/ZKj/OuDtZvZW4BzguSQjgPPN7ECq+seAR8NOKYQYmsZnpbvf4O7H3P0i4Grgm+7+HuBu4Kq02nHg9t56KYTolFXerX0UuMXMPgl8D7ip6YBkcW+gcVloeG8sRx8Eh/eO0tFXOK4HR1++7lQcfbB63v4Oh/rL5t2/BXwr/fwg8Oq9HC+EGAfRQ3ajvM7bjdEYe3KV19DYU+yvjD0hRMnbL5OOEKKOqIoPMB9qjl9mncJ7x2jsyfVPxp690X94bzMjvCxCiL6JPsefzQae45cZTPmXhWtp7Mn1RcaedvQR3qtEHEKIWiIrPizmG+zBSxCPyMYeUN5+UN5+6D5vfwgjvAxCiL6Jvqq/mBvZ82bUyg+9GntAefvz29Hz9o9U8lbO2685vhCiDt34QkyQuEN9NxazDarGW6Me9vdh7MnVVd7+5Xa08N79auyxsEW/Ef7IQoi+iaz4wNyqn1rrsODXpbEHlLc/+QAMEN67j409IYzwRxVC9E18xc/N8SufWuus/IUy5e1X3v7uCDb26HWeEKKO6AE8LMjNRXaxI0r5lbe/z7z9+0b5oTCCluILIeqI/h7fZrYMMZzVPLVy22ul/Pl9lMqUt79w7sHz9u9XY4/e4wsh6tCNL8QEif46z2bLsVjBVVRa8FvLIT8ob3+OUeftb+Poy+8bEdULfrszwh9DCNE30V/n2bz6UbxzwW935YeRq3+X4b3K2599AJS3v44ZYHqdJ4SoI6rim8PGvPhXQMqEKz+sxbx/HcJ7x2jsyfVvNMaepKPFspFJp+l1nhCijgHm+PmH5irKv9wn5Vfe/iUK7w1hzborhOiCgd7jZ0/3JXXqX6f8yac1fNc/wvDeSebtb6P8SQeLjEz5NccXQtQSX/HnUBWS1jTv33fGHlDe/hzR8/a3MfbA6PP26z2+EKKWoBvfzM43s1vN7Mdmdp+ZvcbMDpvZN8zsp+n/F/TdWSFEN4QO9b8AfM3drzKzQ8C5wMeAu9z9U2Z2PXA98NHdGjGHjVl+lLTzNUzwkB+mE96rvP3Fcw9l7MnXHamxx4LmhQFdMrPnAn8M3ATg7lvu/hRwJXBzWu1m4B2teiqEiE6I4r8EeAL4opm9ErgHuA54obs/BuDuj5nZC0JOaItE9aH8sCw+mbsJ8tkn4b1tjD2gvP3JB2A6xp4uF/cOAK8CbnT3S4GnSYb1gR2xa83slJmdmv3v06GHCSF6JETxTwOn3f1kun0ryY3/SzM7kqr9EeDxqoPd/QRwAuDZL7rQbbY860ZOfnfO+7tU/vwZ1lj5C2XK2y9jT9WpOprju/svgIfN7KXpriuAe4E7gOPpvuPA7XvvphBiCEJX9f8S+FK6ov8g8D6Sh8ZXzOwa4CHgnSENJbbcnWffOe/vTvmTT2sY5DNmYw8Eh/fK2JOVUSzrQflDQ3aDbnx3/z5wWUXRFXvokxBiJERPxGGz6led5Xl/k/IndSYS3jtCYw+se97+FsYeGH3e/s7e4wsh9h+68YWYIAO58xKqXraFDvnzx3cZ3jvqIT8ob3+O6I6+ZdFo8/ZvyI8vhKgjes69jblXKvSOh1+j8kM/4b37xNgD4Qt+E8zb38bYk5SNO7w3MGJXii/EFBkg594yrHA15QcZe1iP8N59YuyB8Qf5aI4vhKhloL+kk5A3FJRVukn5C2Uy9ow7vHffGHuW+8aq/ArgEULUMvBf0llSN++vU37oK7x3nxh78vsolSlvfytjD4w/b7/m+EKIWqKv6uff41fPpgOVH/qJ8tsvxh4Ij/Jrofyw3nn72xl7IDjKr42xB8Kj/Gouieb4QohadOMLMUHiv86bARXD+Z0jmd2H/Pm6XYb37htjD4SH904xb38bY0/huB6MPfm6LY09neXcE0LsPwZ4nedsVKwUlRf8mpS/qu5w4b0TNPbAWuftb2PsKfZ3nMYe/ZlsIUQtkV/nefo6L2GjYsLYh/IXymTsWQ9jD6x5eO8wxh7N8YUQtQyUZTdT7yXleX+T8ufrLNso1o2Vt1/GnhTl7WdoY48CeIQQtQyQbNNhW83z6k26L0z588d3Gd4rY0+pzkDGHlDe/koajD1d/rVcIcQ+Y6Bkmxn5yL3ivL9J+Yt1u4zyk7Fnx77CfqXvzp97bOm7taovhKhFN74QE2SALLtevfhWWvBrGvInZZTqdhHeG8fYAyMf9ncZ3qu8/dkHoN/w3o3K93w7keILMUHiBvAAtvBtqatcfAtU/qRsnY09yxbXUvkLZcrbP5a8/YFv86T4QkyR6CYdmy3gQPq8yUndTj1sUn7oI7xXeftLrEN4774x9iz3tVV+hewKIWoJUnwz+xDw5yTPnR8A7wOOALcAh4HvAu91963aRiAXsps+xQ7knjuleX+T8if7egjvjWTsST6tYZDPCMN794uxB1bP299ZIg4zOwp8ALjM3V8BbAJXA58GPufulwBPAtcEnVEIMTihc/wDwLPN7CxwLvAY8Abg3Wn5zcAngBt3bcVhY7aoXuEuzfublD9f1mV4bzxjz7LFtVR+CA/vbaH80HPe/lEaeyA4vLfG2NPZqr67PwJ8BniI5Ib/NXAP8JS7Z9/N08DRquPN7FozO2Vmp87Ong7slhCiT0KG+hcAVwIXAy8GzgPeUlG18rnu7ifc/TJ3v+zggfNW6asQoiNChvpvBH7m7k8AmNltwGuB883sQKr6x4BHm5tyWCxyQ+eKcWPgkD8p6z68N5qjD5S3v4pW4b37xNFXOK6do6/L13kPAZeb2blmZsAVwL3A3cBVaZ3jwO1BZxRCDE6j4rv7STO7leSV3Qz4HnAC+HfgFjP7ZLrvpqa2kpx7C7wyWKYkH43KD+MJ7927sQf2Et47QWMP9Jq3f4zGnmJ/2xl7Qhf3glb13f3jwMdLux8EXh14HiHEiIiegYe5Y6kEeGWwTJjy52rQZXivjD0rIGNPJTHz9isDjxCilvgmnfkcZxNgW/mBHfP+JuWHfsJ7ZezpABl7KomVtz8EKb4QEyR66i3mi9yDe3O7qDzvb1R+6CW8N5axp6rdDBl7GhjI2ANjCu+tXtdQIg4hRC2RV/Udm80rnYXleX+j8kMvUX6xjD3F42XsWSnKL5axJ3fc4FF+NcYeC/xTOlJ8ISaIbnwhJkj8xb3FApstNzPKC37NQ/5kC+g0vDeesWd5zj6MPTDyYX+X4b2RjD1J98YR3ltn7NHinhCiluiLe8xz2jcrlCT7trd3V37oK7w3lrFn2ReF91ItQT2G97Yx9iTdHEl4767Gnmak+EJMkEHm+FWU5/1Nyg/9hPfGMvZU9UXhvURX/qRsDcN7a4w9muMLIWoZYI4/37VGqPJDP+G9sYw9xX4X+9KF8ief1jDIJ1J4bytjDwSH9w5l7LFAzZfiCzFBottyfTYLeiY1KX+yr/vw3mjGHlDe/irahPdGMvbAmuTtD0CKL8QE0Y0vxAQZxI/v6SBzlSF/sq+P8N44jr7COZW3v5rQ8N5Yjj4IDu8dytGn13lCiFoGfZ3nOa1pelKtYuyB8PDeWMaeynMqb381ozH25CqP2NgTghRfiAkSfY7v83mlRoXO+9sYeyA8vFfGnhHP+9cpvHdAY08IUnwhJkj8v6Tjjqfz/FWUH3oK741k7MnXqT2fjD3VDKb8y8KxGnsUsiuEqCWq4nuq9raZKLLnVvjLz6kulT/ZFxjeG8nYUzy+e2NPVbsZMvY00EF475B5+0OQ4gsxQQZ5j7+tzJtLRa6b98eP8otj7CnW7d7YUzxexp4ulB/WK2//bkjxhZgguvGFmCCDBvAUwm9LC35dDvnzx43F2FPVjvL2d0Cfxp5c3XXK21+FFF+ICTKaAJ7ygl+T8leVlYmVt7+NsScpU3hvb/Rh7IG1z9u/PI8QYnKYe7iVb+WTmT0BPA38KtpJV+P5rE9fYb36u059hfXp7++4+281VYp64wOY2Sl3vyzqSVuyTn2F9ervOvUV1q+/TWioL8QE0Y0vxAQZ4sY/McA527JOfYX16u869RXWr7+7En2OL4QYHg31hZgg0W58M3uzmd1vZg+Y2fWxzhuKmV1oZneb2X1m9iMzuy7df9jMvmFmP03/v2DovmaY2aaZfc/M7ky3Lzazk2lfv2xmh4buY4aZnW9mt5rZj9Nr/JqxXlsz+1D6Hfihmf2zmZ0z5mvbhig3vpltAn8HvAV4OfAuM3t5jHPvgRnwYXd/GXA58P60j9cDd7n7JcBd6fZYuA64L7f9aeBzaV+fBK4ZpFfVfAH4mrv/HvBKkn6P7tqa2VHgA8Bl7v4KYBO4mnFf273j7r3/A14DfD23fQNwQ4xzr9Dn24E3AfcDR9J9R4D7h+5b2pdjJDfLG4A7SYJJfwUcqLrmA/f1ucDPSNeUcvtHd22Bo8DDwGGSgOg7gT8Z67Vt+y/WUD+7mBmn032jxMwuAi4FTgIvdPfHANL/XzBczwp8HvgIy6jz5wFPuXsWOj+ma/wS4Angi+nU5O/N7DxGeG3d/RHgM8BDwGPAr4F7GO+1bUWsG786oe4IMbPnAP8KfNDdfzN0f6ows7cBj7v7PfndFVXHco0PAK8CbnT3S0nCtgcf1leRrjNcCVwMvBg4j2SKWmYs17YVsW7808CFue1jwKORzh2MmR0kuem/5O63pbt/aWZH0vIjwOND9S/H64C3m9nPgVtIhvufB843s8xxOaZrfBo47e4n0+1bSR4EY7y2bwR+5u5PuPtZ4DbgtYz32rYi1o3/HeCSdGX0EMliyR2Rzh2EmRlwE3Cfu382V3QHcDz9fJxk7j8o7n6Dux9z94tIruU33f09wN3AVWm1UfQVwN1/ATxsZi9Nd10B3MsIry3JEP9yMzs3/U5kfR3ltW1NxEWTtwI/Af4L+OuhFzcq+vdHJMO3/wS+n/57K8nc+S7gp+n/h4fua6nfrwfuTD+/BPg28ADwL8Czhu5frp9/AJxKr++/AReM9doCfwP8GPgh8E/As8Z8bdv8U+SeEBNEkXtCTBDd+EJMEN34QkwQ3fhCTBDd+EJMEN34QkwQ3fhCTBDd+EJMkP8HFaW3BVZhaXoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f_mean = np.zeros([100])\n",
    "def kernel(x1, x2):\n",
    "    return gaussian(50)(x1, x2)\n",
    "K = compute_kernel(x, x, kernel) #+ np.identity(100)*0.1**2\n",
    "plt.imshow(K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x116b69978>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4VOXB/vHvk52EJYSELRAChD2yhl3rrqhVEbGVzV3crfW1orbV2qqtYt+qVVAUi1rAhcWNiltVVDQQthCWsAYIhJAEspCFkOT5/ZHojxeBRJjJmTlzf66LC2bmZM59OOHm5JnnnGOstYiIiLsEOR1AREQ8T+UuIuJCKncRERdSuYuIuJDKXUTEhVTuIiIupHIXEXEhlbuIiAup3EVEXCjEqRXHxsbaxMREp1YvIuKXVqxYkW+tjatvOcfKPTExkbS0NKdWLyLil4wxOxqynIZlRERcSOUuIuJCKncRERdSuYuIuJDKXUTEhVTuIiIupHIXEXEhx+a5i4j4G2stucWH2J5fSkHpIfaXVlJYdpjqGssPNyxtHhFCTFQYLaPCSIiJJCEmktDgxj+OVrmLiBxHbnEFy7bvJy1rP6uzi9i67yAHD1X9rPcICTJ0ahVJ3w7RDE6MYXBiS7rGNSUoyHgpdd16vfruIiJ+xFpLxu5iPlm/l0/W5ZKZWwJAZFgwfTu04MqB8SS1bkrn2KbENQsnJiqM6MhQQoIMxhhqaiwlFVUcKKukoPQQWfllbMs/yKbcg3y9OZ+Fq3YDcP3IRB65tI9Xt0XlLiIBb29RBfNXZvNO2i6yCsoIMjCkcwy/v7gXQ7vE0Ktd8wYNrQQFGVpEhtIiMpTE2CgGdYr58TVrLTsKylietZ8ucU29uTmAyl1EAtiKHQeYsWQrn67PpcbC0M4x3H5WEuf1bkNMVJhH12WMITE2isTYKI++7/Go3EUkoFhr+XJTHtO+2MLyrAO0aBLKLWd25dcpHRuteBuDyl1EAsbyrP08tXgjy7MOEB/dhEcu7c2vUjoSFe6+KnTfFomIHGVnQRl//nA9n23IJa5ZOH8ZncyvUzoSFuLeU31U7iLiWhWHq5n+5Vamf7WVkCDD7y7swfUjE4kMc3/1uX8LRSQgpW4r4P756ewoKOPSfu35/cW9aNsiwulYjabecjfGdAReB9oCNcAMa+2zRy1jgGeBi4Ey4Dpr7UrPxxURObGyyiqeWpzJrKVZJMREMuemoYxIinU6VqNryJF7FfA/1tqVxphmwApjzKfW2vVHLHMR0K3u11Bget3vIiKNZm12EXfNXUlWQRnXjUjk/lE9AmII5ljq3WprbQ6QU/fnEmPMBiAeOLLcLwdet9Za4HtjTLQxpl3d14qIeJW1lpnfbOfJxRuJbRrO3JuHMbxrK6djOepn/ZdmjEkEBgCpR70UD+w64nF23XMqdxHxqqKyw9z79mo+37iP83u3YerYvkRHevYEJH/U4HI3xjQF5gP3WGuLj375GF9ij37CGDMZmAyQkJDwM2KKiPxU5t4SJr+Rxp7Cch69rA/XDO9E7UeA0qBJnsaYUGqLfba1dsExFskGOh7xuAOw5+iFrLUzrLUp1tqUuLi4k8krIgLAovQcrpj2LeWV1bw5eTjXjkhUsR+h3nKvmwkzE9hgrf3f4yz2PnCNqTUMKNJ4u4h4g7WWZz/bzB1zVtKrXXM+vOt0BnVq6XQsn9OQYZmRwCRgrTFmdd1zDwEJANbaF4H/UDsNcgu1UyGv93xUEQl0lVU1PLAgnQUrdzNmYDx/HXMa4SHBTsfySQ2ZLfMNxx5TP3IZC9zhqVAiIkcrKj/MLW+k8f22/dx7fnfuOidJwzAnEJgTQEXEr+wrruCaV5exNe8gz/y6P6MHxDsdyeep3EXEp+0oKGXSzGXkHzzEq9cN5oxumozRECp3EfFZmXtLmPBKKlU1Ncy5eRj9O0Y7HclvqNxFxCet21PExFdSCQsJYu7Nw+nWppnTkfyKyl1EfM7a7CImzkwlKiyYOTcPc9UdkhqLyl1EfMqaXYVMnJlK84hQ3pw8jI4xkU5H8ksqdxHxGev2FHHNq8uIjgzlzcnDiY9u4nQkv+Xee0yJiF/ZlFvCpJnLaodibhqmYj9FKncRcdz2/FLGv5xKSJBh9s0aivEElbuIOGpvUQUTX0mlxlrm3DyUzvrw1CNU7iLimMKySibNTKWo/DCvXT+EpNaa7ugpKncRcURZZRXXz1rOjoIyZlwziNM6tHA6kquo3EWk0VVV13DH7JWs2VXIc+MGMKJr4N3A2ts0FVJEGpW1lj++t44vMvN4/IpkRiW3dTqSK+nIXUQa1bQvtzJ32U5uP6srE4Z2cjqOa6ncRaTRvLtqN1M/zmR0//b87sIeTsdxNZW7iDSK5Vn7uX9eOsO6xPDU2H660YaXqdxFxOt2FpRxyxsriG/ZhBcnDiIsRNXjbfobFhGvKq44zA2vLae6xjLz2hSiI8OcjhQQVO4i4jXVNZa75qwiK7+U6RMH0iWuqdORAoamQoqI1zy1eCNfbcrjiStO01z2RqYjdxHxindX7ealJduYNKwT44cmOB0n4KjcRcTj0rMLmTI/nSGdY3j40t5OxwlIKncR8aiCg4e49Y0VxDYNZ/qEgYQGq2acoDF3EfGYquoa7pq7ioLSSubfNoJWTcOdjhSwVO4i4jFTP8lk6dYCpo7tS3K8rvLoJP28JCIe8dHaHF76ahsThiZwVUpHp+MEPJW7iJyybXkH+d28dPp3jNYHqD5C5S4ip6S8sprbZ68kNNgwbcJAwkOCnY4kaMxdRE7Rw+9lkJlbwr+uG0z76CZOx5E6OnIXkZP2dtou3lmRzZ1nJ3FWj9ZOx5EjqNxF5KRk7i3h4fcyGNG1Ffec193pOHIUlbuI/GxllVXcOWclTcNDeObq/gQH6drsvqbecjfGvGqM2WeMyTjO6y2MMR8YY9YYY9YZY673fEwR8SV/en8dW/IO8syvB9C6WYTTceQYGnLkPgsYdYLX7wDWW2v7AWcBfzfG6ILNIi61cFU2b6dlc8dZSZzeTVd69FX1lru1dgmw/0SLAM1M7T2zmtYtW+WZeCLiS7LyS/nDwgwGJ7bknvO6OR1HTsATUyGfB94H9gDNgF9ba2s88L4i4kMqq2q4+81VhAQH8ezVAwjRBcF8mif2zoXAaqA90B943hjT/FgLGmMmG2PSjDFpeXl5Hli1iDSWv3+SSXp2EU9eeZrms/sBT5T79cACW2sLsB3oeawFrbUzrLUp1tqUuLg4D6xaRBrD15vzeGnJNsYPTWBUcjun40gDeKLcdwLnAhhj2gA9gG0eeF8R8QEFBw9x79tr6Na6KX+8RNeN8Rf1jrkbY+ZSOwsm1hiTDTwChAJYa18E/gLMMsasBQwwxVqb77XEItJorLVMmb+WorLDvH7DEJqE6box/qLecrfWjqvn9T3ABR5LJCI+Y+6yXXy2IZc/XNKLXu2O+VGa+Ch93C0ix7Q17yB//nAdZ3SL5YaRnZ2OIz+Tyl1EfqKyqoZ73lxNk9Bgnr6qH0G6vIDf0SV/ReQnnvt8M2t3F/HixIG0aa7LC/gjHbmLyP+xYsd+pn25hasGddC0Rz+mcheRH5UequK3b62hfXQT3S7Pz2lYRkR+9Nii9ew6UMZbk4fTLCLU6ThyCnTkLiIA/HdjLnOX7WLyGV0Y0jnG6ThyilTuIsKB0kqmzF9Lz7bNuPcC3VXJDTQsIyL84b0MCssqmXX9YMJDdBaqG+jIXSTAvb9mD4vSc7jnvO70ad/C6TjiISp3kQCWW1zBH9/NYEBCNLf8oovTccSDVO4iAcpay4ML1nKoqpqnr+qnm2+4jPamSIB6Z0U2/924j/sv7EnXuKZOxxEPU7mLBKA9heX85YP1DOkcw3UjEp2OI16gchcJMLXXaE+n2lqeHquLgrmVyl0kwMxZtpOvN+fz4MW9SGgV6XQc8RKVu0gAyT5QxhOLNjAyqRUThyY4HUe8SOUuEiCstTwwfy0AfxvTF2M0HONmKneRADFn2U6+2VI7HNMxRsMxbqdyFwkARw7HTNBwTEBQuYu43A8nK4GGYwKJyl3E5d5avouvN+fzgIZjAorKXcTF9hSW89iiDQzv0ooJQzQcE0hU7iIu9cNwTHWN5ckr++pkpQCjchdxqXkrsvlqUx5TRvXQyUoBSOUu4kK5xRX85cP1DEmM4ZrhiU7HEQeo3EVcxlrL7xdmcKiqhifHajgmUKncRVzmg/QcPtuQy30X9KBzbJTTccQhKncRFyk4eIg/vb+Ofh2jueH0zk7HEQep3EVc5JH313GwooqpY/sSrOGYgKZyF3GJj9ft5cP0HO4+N4nubZo5HUccpnIXcYGissP84d0Merdrzi1ndnU6jviAEKcDiMipe2zRevaXVvKv6wYTqhtdCw04cjfGvGqM2WeMyTjBMmcZY1YbY9YZY77ybEQROZGvNuXxzopsbj2zC8nxLZyOIz6iIf/FzwJGHe9FY0w0MA24zFrbB7jKM9FEpD4HD1Xx0IK1dI2L4q5zujkdR3xIveVurV0C7D/BIuOBBdbanXXL7/NQNhGpx1OLN7KnqJynxvYjIjTY6TjiQzwxONcdaGmM+dIYs8IYc40H3lNE6rFs+35e/24H141IZFCnlk7HER/jiQ9UQ4BBwLlAE+A7Y8z31tpNRy9ojJkMTAZISNDlR0VOVsXhaqbMT6djTBN+d2EPp+OID/LEkXs2sNhaW2qtzQeWAP2OtaC1doa1NsVamxIXF+eBVYsEpn98uont+aU8OaYvkWGa9CY/5Ylyfw84wxgTYoyJBIYCGzzwviJyDGt2FfLy19sYN6QjI5JinY4jPqre//KNMXOBs4BYY0w28AgQCmCtfdFau8EYsxhIB2qAV6y1x502KSInr7Kqhinz02ndLIIHL+7ldBzxYfWWu7V2XAOWmQpM9UgiETmuaV9uYePeEmZem0LziFCn44gP06lsIn4ic28JL3yxhcv7t+fcXm2cjiM+TuUu4geqqmu4f94amkeE8silfZyOI35AH7OL+IFXv93Omuwi/jluADFRYU7HET+gI3cRH7ct7yB//2QTF/Ruwy/7tnM6jvgJlbuID6upsUyZn054SBCPjU7GGN2AQxpG5S7iw/6duoPlWQf44y9707p5hNNxxI+o3EV81K79Zfzto42c0S2WsYM6OB1H/IzKXcQHWWt5cMFaDPC3K/tqOEZ+NpW7iA96a/kuvtmSz4MX9yI+uonTccQPqdxFfExOUTmPL9rAsC4xjB+iq6fKyVG5i/gQay0PLVhLVY3lySv7EhSk4Rg5OSp3ER8yf+VuvsjM474Le9CpVZTTccSPqdxFfERucQV//mAdgxNbcv2IRKfjiJ9TuYv4AGstv1+4lkNVNTw1tp+GY+SUqdxFfMB7q/fw2YZ9/O7CHnSO1XCMnDqVu4jD9pVU8KcP1jEwIZrrR3Z2Oo64hMpdxEG1wzEZlFdWM/WqfgRrOEY8ROUu4qB3V+/m0/W53HdBD7rGNXU6jriIyl3EIbnFFfzp/fUM6tSSG07XcIx4lspdxAE/nKxUcbiaqWP7ajhGPE7lLuKA+St38/nG2tkxXTQcI16gchdpZHsKy3n0/XUMSYzR7BjxGpW7SCOytvbOSlU1lqlXaThGvEflLtKI5izbydeb83nokl66dox4lcpdpJHsLCjjiUUbOD0plolDdSlf8S6Vu0gjqK6x3PfOGoKM4cmxurOSeF+I0wFEAsHMb7axLGs/f7+qn+6sJI1CR+4iXpa5t4SnP97EhX3aMGZgvNNxJECo3EW8qLKqhnvfXk2ziBCeuOI0DcdIo9GwjIgXPfv5JtbtKealSYNo1TTc6TgSQHTkLuIly7P2M/3LrfwqpQMX9mnrdBwJMCp3ES8oqTjMb99aTYeWkTx8aR+n40gA0rCMiBc8+sF69hSW886tw2karn9m0vjqPXI3xrxqjNlnjMmoZ7nBxphqY8xYz8UT8T8frc1h3opsbj8riUGdYpyOIwGqIcMys4BRJ1rAGBMMPAl87IFMIn4rp6icBxaspW+HFtx9bjen40gAq7fcrbVLgP31LHYXMB/Y54lQIv6opsZy71trqKyq4dmrBxAWoo+0xDmn/N1njIkHrgBebMCyk40xacaYtLy8vFNdtYhPefnrbXy3rYA/XdabzrG6KJg4yxOHFs8AU6y11fUtaK2dYa1NsdamxMXFndTKtuUdZOIrqewtqjiprxfxhozdRTz9SSYXJbflVykdnY4j4pFyTwHeNMZkAWOBacaY0R5432PaXVjOqp0HuOz5b1izq9BbqxFpsIOHqrhr7ipaRYXrLFTxGadc7tbaztbaRGttIjAPuN1a++4pJzuOM7rFMf/2EYSFBPGrl77jvdW7vbUqkQZ55L117Cgo5Zmr+9MyKszpOCJAw6ZCzgW+A3oYY7KNMTcaY241xtzq/XjH1rNtc967YyT9OkTzmzdX849PN2GtdSqOBLB3V+1m/sps7jynG8O6tHI6jsiP6j27wlo7rqFvZq297pTS/Aytmobz75uG8uCCtTz7+WayCkp58sq+RIQGN1YECXA7Ckr5w7sZDE5syd3nJDkdR+T/8OtT58JCgnj6qr50iYti6seZ7D5QzoxrUojRj8biZYeqqrljzkqCgwzPXD2AkGBNexTf4vffkcYY7jg7iRfGD2Tt7iLGTPuW7fmlTscSl/vrfzaSsbuYp3XzDfFRfl/uP7ikbzvm3DyM4ooqxkz7lrSs+s67Ejk5H63NYdbSLG48vTPn927jdByRY3JNuQMM6tSShbePIDoyjPGvpLIoPcfpSOIyOwvKuH9+Ov06RjNlVE+n44gcl6vKHaBTqygW3DaCvvEtuGPOSl75epvTkcQlKg5Xc/ucFRjg+XG6vID4Nld+d7aMCuPfNw3louS2PLZoA49+sI7qGk2VlFPz5w/Xk7G7mL//qj8dYyKdjiNyQq4sd4CI0GCeHz+Q60cm8q9vs7hr7koqDtd7hQSRY1qwMps5qTu59cyuGmcXv+DXUyHrExxkeOTSPrRv0YTH/7OB/IPLeHlSCi0iQ52OJn4kc28JDy1cy9DOMdx3QXen44g0iGuP3I908y+68Ny4AazaeYCrXlrKnsJypyOJnygqP8wtb6TRLCKUf47XfHbxHwHznXpZv/a8dv0QcgorGDNtKZtyS5yOJD6u9vrsq8k+UM60CQNp3SzC6UgiDRYw5Q4wIimWt28dTo21jJ2+lGXbNRdeju/Zzzfz+cZ9PHxpbwYn6nZ54l8CqtwBerVrzoLbRxDbLJyJM1NZnKG58PJTn63P5dnPN3PlwA5MGtbJ6TgiP1vAlTtAh5aRzL91BMntm3Pb7JW88f0OpyOJD9mcW8I9b60mOb45j1+RrOuzi18KyHKH2rnws28axrk9W/PHdzN4+uNMXTZYKCyr5KbX04gIDWbGpBRdZVT8VsCWO0CTsGBenDiIqwd35PkvtnD/vHQOV9c4HUscUlVdwx1zVpJTWMFLkwbRXhcEEz/m6nnuDRESHMRfx5xG6+YRPPf5ZvIPHuKFCQOJDAv4v5qA89iiDXy7pYCnr+rHoE4tnY4jckoC+sj9B8YY7j2/O09ccRpfbcpj3IzvKTh4yOlY0oheW5rFrKVZ3HxGZ8YO6uB0HJFTpnI/wvihCbw0KYXM3BKunL6UHQW6Lnwg+GLjPh79YB3n927DAxf1cjqOiEeo3I9yfu82zL5pGEXlhxkzbSlrdhU6HUm8aP2eYu6cs5Je7Zrz7NX9CQ7SzBhxB5X7MQzq1JJ5t40gMjyYq2d8z3835jodSbxgT2E5N8xaTrOIUGZeO1ifs4irqNyPo2tcU+bfNoKuraO4+fUVzEnd6XQk8aCi8sNc969lHDxUxavXDaZtC11aQNxF5X4CrZtF8Nbk4ZzRLZaHFq5l6scbNRfeBQ5VVTP59TS255fy0qRB9G7f3OlIIh6ncq9HVHgIr1yTwrghCbzwxVbueWs1h6p0XXh/VV1jufftNaRu38/Usf0YmRTrdCQRr9AgYwOEBAfxxBXJdGjZhKkfZ/54kkvLqDCno8nPYK3l4fcyWJSew4MX9WT0gHinI4l4jY7cG8gYwx1nJ/HcuAGszi5kzPSlZOVrqqQ/efqTTGbX3U3pljO7Oh1HxKtU7j/TZf3aM+emoRSWVXLFtG9J3VbgdCRpgBlLtvLCF1sZNySBKaN6OB1HxOtU7ichJTGGhbePpGVUGBNnpjJvRbbTkeQEXluaxRP/2cglfdvx2Ghd5VECg8r9JCXGRrHwtpEM6RzDfe+s4cnFG6mp0UwaXzM7dQePvF979ukzv9ZJShI4VO6noEVkKLOuH8L4oQlM/3Irk99YwcFDVU7HkjpvL9/F7xdmcE7P1jw/fgChuv+pBBB9t5+i0OAgHh+dzKOX9eGLzH2MmfYtOwvKnI4V8Oak7mTKgnTO6BbLtAkDCQ/RddklsKjcPcAYw7UjEnn9hiHkFh/ishe+4evNeU7HClizvt3OQwvXcmb3OF6+RjfckMCkcvegkUmxvH/nSNo0i+DaV5fx4ldbdUZrI5uxZCt/+mA95/duw0uTBqnYJWCp3D2sU6soFtw+gotOa8ffPtrInXNWaRy+EVhreXLxxh9nxWgoRgJdveVujHnVGLPPGJNxnNcnGGPS634tNcb083xM/xIVHsLz4wbw4EU9+Sgjh8ue/4ZNuSVOx3KtquoapsxPZ/qXWxk/NIHnrtaHpyIN+RcwCxh1gte3A2daa/sCfwFmeCCX3zPGcMuZXZl90zCKy6u4/PlvWbBS8+E9rbyymttmr+TttGzuPieJx0cna7qjCA0od2vtEmD/CV5faq09UPfwe0D3KDvC8K6t+M/dp9O3QwvufXsN989bQ1mlhmk8YV9xBVfP+I7PNuTy6GV9uPeCHjpBSaSOp392vRH46HgvGmMmG2PSjDFpeXmBM5ukdfMIZt80lDvPTuKdFdlc+s9v2JBT7HQsv7ZxbzGjX/iWTbkHmTEphWtHJDodScSneKzcjTFnU1vuU463jLV2hrU2xVqbEhcX56lV+4WQ4CDuu7AHs28cSnFFFZe/8C3/+na7zmo9CR+v28vY6d9RbS3v3Dqc83u3cTqSiM/xSLkbY/oCrwCXW2t1Ja0TGJEUy+LfnMHpSbE8+sF6rnl1GTlF5U7H8gvVNZanP87kljdW0CUuinfvGElyfAunY4n4pFMud2NMArAAmGSt3XTqkdyvVdNwZl6bwhNXnMaKHQe48B9LWLgqW3PiT+BAaSU3vrac57/Ywq9SOvD2LcNp16KJ07FEfJapr1CMMXOBs4BYIBd4BAgFsNa+aIx5BbgS2FH3JVXW2pT6VpySkmLT0tJOPrlLZOWXcu/bq1m5s5Bzerbm8SuSVVpHSd1WwD1vrSb/4CEeubQPE4Ym6INTCVjGmBUN6dh6y91bVO7/X3WNZdbSLKZ+vJHQoCCmXNST8UMSCArwKX1V1TX8879b+Od/N5MQE8k/xw3ktA4ahpHApnL3QzsKSnlwwVqWbi2gX8doHh+dHLBjyptzS7hvXjprdhUyZmA8f748mabhuiukiMrdT1lreW/1Hh5btJ79pZVMHNaJe87rTkyA3K+1qrqGl7/ezj8+3URUeDB/vjyZS/u1dzqWiM9oaLnrUMjHGGMYPSCes3u05ulPMvn39zt4d9VufnNedyYN60RYiHtPq0/L2s8f31vHhpxiLkpuy19GJxPbNNzpWCJ+SUfuPi5zbwmPLVrP15vzSYiJ5Lfnd+OyfvGuOsV+X0kFTy3OZN6KbNq1iODhX/ZmVHJbfWgqcgwalnERay1fbspj6uJM1ucU071NU35zbndGJbf165IvqTjMy0u28co32zlcXcNNZ3ThrnOSiAzTD5Qix6Nyd6GaGstHGXv5308z2ZpXSpfYKG45swujB8T71eVtDx6qYvb3O3hpyTb2l1ZySd923HdBDzrHRjkdTcTnqdxdrLrGsjhjL9O+3MK6PcXENQtn3OCOjB/aibYtIpyOd1z5Bw/x2tIsXluaRXFFFacnxfK7C3vQr2O009FE/IbKPQBYa/l6cz6zlmbxReY+gozh3J6tGTOwA2f3jPOJo3lrLanb9zM7dSeLM3I4XG25sE8bbjsrif4qdZGfTbNlAoAxhl90j+MX3ePYWVDG7NQdzF+5m0/W59KiSSgXJbflgj5tGNE1tlFvN2etZUNOCR+m72HR2hx2FJTRPCKEicM6MWFoJ5JaN220LCKBSkfuLlNVXcM3W/JZuGo3n63PpbSymiahwYxMimV411YM7RxDr3bNPf5BbGFZJd9vK2DJ5ny+2ZzPzv1lBAcZRnRtxeX94/ll33a6n6mIB+jIPUCFBAdxVo/WnNWjNYeqqvl+234+W5/LV5vy+GxDLgBNw0Po1a4Zvdo1p2fb5nSMaUKHlpG0axFxwgK21lJcUUX2gTJ2FpSxvaCU9XuKWbu7iB0FZT++97Aurbj1zK6MSm4bMCdfifgaHbkHkD2F5SzP2k9a1gE25BSzcW/JT27eHREaRLOIUJqFh2AMWAvV1lJSUUVR+WGqj7r+fHx0E/p2aMFpHVqQ0imGAQnRun+piBfpyF1+on10Ey7vH8/l/eOB2qmVe4rKyT5Qzu4D5ewpLKe44jAlFVWU1JW+AYKMoVlECC0jw4iODKV9dBMSYiLp1CqSZhGhDm6RiByPyj2ABQUZOrSMpEPLSKejiIiH6ednEREXUrmLiLiQyl1ExIVU7iIiLqRyFxFxIZW7iIgLqdxFRFxI5S4i4kKOXX7AGJMH7DjJL48F8j0Yx18E4nYH4jZDYG53IG4z/Pzt7mStjatvIcfK/VQYY9Iacm0FtwnE7Q7EbYbA3O5A3Gbw3nZrWEZExIVU7iIiLuSv5T7D6QAOCcTtDsRthsDc7kDcZvDSdvvlmLuIiJyYvx65i4jICfhduRtjRhljMo0xW4wxDzidxxuMMR2NMV8YYzYYY9YZY35T93yMMeZTY8zmut9bOp3VG4wxwcaYVcaYD+sedzbGpNZt91vGGFfdu88YE22MmWeM2Vi3z4cHwr42xvy27vs7wxgz1xgT4cZ9bYx51Ri+WbkTAAADMElEQVSzzxiTccRzx9y/ptZzdf2WbowZeLLr9atyN8YEAy8AFwG9gXHGmN7OpvKKKuB/rLW9gGHAHXXb+QDwubW2G/B53WM3+g2w4YjHTwL/qNvuA8CNjqTynmeBxdbankA/arfd1fvaGBMP3A2kWGuTgWDgaty5r2cBo4567nj79yKgW92vycD0k12pX5U7MATYYq3dZq2tBN4ELnc4k8dZa3OstSvr/lxC7T/2eGq39bW6xV4DRjuT0HuMMR2AS4BX6h4b4BxgXt0irtpuY0xz4BfATABrbaW1tpAA2NfU3gmuiTEmBIgEcnDhvrbWLgH2H/X08fbv5cDrttb3QLQxpt3JrNffyj0e2HXE4+y651zLGJMIDABSgTbW2hyo/Q8AaO1cMq95BrgfqKl73AootNb+cCdvt+3zLkAe8K+6oahXjDFRuHxfW2t3A08DO6kt9SJgBe7e10c63v71WMf5W7mbYzzn2uk+xpimwHzgHmttsdN5vM0Y80tgn7V2xZFPH2NRN+3zEGAgMN1aOwAoxWVDMMdSN8Z8OdAZaA9EUTskcTQ37euG8Nj3u7+VezbQ8YjHHYA9DmXxKmNMKLXFPttau6Du6dwffkSr+32fU/m8ZCRwmTEmi9oht3OoPZKPrvvRHdy3z7OBbGttat3jedSWvdv39XnAdmttnrX2MLAAGIG79/WRjrd/PdZx/lbuy4FudZ+oh1H7Acz7DmfyuLpx5pnABmvt/x7x0vvAtXV/vhZ4r7GzeZO19kFrbQdrbSK1+/a/1toJwBfA2LrFXLXd1tq9wC5jTI+6p84F1uPyfU3tcMwwY0xk3ff7D9vt2n19lOPt3/eBa+pmzQwDin4YvvnZrLV+9Qu4GNgEbAV+73QeL23j6dT+KJYOrK77dTG148+fA5vrfo9xOqsX/w7OAj6s+3MXYBmwBXgHCHc6n4e3tT+QVre/3wVaBsK+Bh4FNgIZwBtAuBv3NTCX2s8VDlN7ZH7j8fYvtcMyL9T121pqZxOd1Hp1hqqIiAv527CMiIg0gMpdRMSFVO4iIi6kchcRcSGVu4iIC6ncRURcSOUuIuJCKncRERf6fzBMZfJEbWPLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113bafcc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f1 = np.random.multivariate_normal(f_mean, K)\n",
    "plt.plot(f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 8) (8, 8)\n",
      "k_x_x.shape (100, 100)\n",
      "k_x_D.shape (100, 8)\n",
      "kk.shape (100, 100)\n",
      "cov.shape (100, 100)\n",
      "[ 1.00018483  0.99668407  0.98270263  0.9388855   0.83215529  0.63348757\n",
      "  0.36045938  0.10705399  0.001999    0.10704719  0.36035731  0.632714\n",
      "  0.82790428  0.9206692   0.92066913  0.82790405  0.63271349  0.36035654\n",
      "  0.10704655  0.001999    0.10704655  0.36035653  0.6327134   0.82790352\n",
      "  0.9206669   0.92066166  0.82788424  0.63267313  0.36029834  0.10700324\n",
      "  0.001999    0.1067278   0.3572308   0.61696102  0.77227849  0.77226615\n",
      "  0.61692286  0.35717309  0.10668456  0.001999    0.10668453  0.3571727\n",
      "  0.61692086  0.77225896  0.77225896  0.61692085  0.3571727   0.10668452\n",
      "  0.001999    0.10668452  0.3571727   0.61692085  0.77225896  0.77225896\n",
      "  0.61692085  0.3571727   0.10668452  0.001999    0.10668452  0.3571727\n",
      "  0.61692085  0.77225896  0.77225896  0.61692086  0.3571727   0.10668453\n",
      "  0.001999    0.10668456  0.35717309  0.61692288  0.77226623  0.77227873\n",
      "  0.61696153  0.35723157  0.10672845  0.001999    0.10701068  0.36040117\n",
      "  0.6334472   0.83213548  0.93887804  0.98270042  0.99668355  1.00018473\n",
      "  1.0008767   1.00098507  1.00099855  1.00099989  1.00099999  1.001       1.001\n",
      "  1.001       1.001       1.001       1.001       1.001       1.001       1.001\n",
      "  1.001       1.001     ]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0, 100)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7IAAADGCAYAAADv27y+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XecVNX5+PHPnZ2Z7b2yBXaBpfcuqIAIodjRaOwxlsSW3r+JX01MLDF+k5+J0RhjiZJYYi8IKqCIwNJZyu6yC2zvOzu7O33O748zW0Dq9oXn/Xpd7sww5c7dOzP3Oc85zzGUUgghhBBCCCGEEAOFqa83QAghhBBCCCGEOB0SyAohhBBCCCGEGFAkkBVCCCGEEEIIMaBIICuEEEIIIYQQYkCRQFYIIYQQQgghxIAigawQQgghhBBCiAHlpIGsYRgZhmF8ahjGXsMwcg3D+G7g9jjDMFYZhpEfWMce5/E3Be6TbxjGTd39BoQQQgghhBBCnF2Mk80jaxjGIGCQUmqrYRiRwBbgMuBmoE4p9ZBhGD8DYpVSPz3qsXFADjANUIHHTlVK1Xf7OxFCCCGEEEIIcVY4aUZWKVWulNoauGwH9gJpwKXA84G7PY8Obo/2NWCVUqouELyuAhZ3x4YLIYQQQgghhDg7ndYYWcMwMoHJwEYgWSlVDjrYBZKO8ZA0oLjD9ZLAbUIIIYQQQgghRKeYT/WOhmFEAK8D31NKNRqGcUoPO8Ztx+zLbBjG7cDtAOHh4VNHjRp1qpsmhBBCCCGEEGKA2LJlS41SKrErz3FKgaxhGBZ0EPuSUuq/gZsrDcMYpJQqD4yjrTrGQ0uAeR2upwNrjvUaSqmngacBpk2bpnJyck7pDQghhBBCCCGEGDgMwzjU1ec4larFBvAPYK9S6o8d/uttoLUK8U3AW8d4+EpgkWEYsYGqxosCtwkhhBBCCCGEEJ1yKmNk5wA3ABcYhrE9sCwFHgIWGoaRDywMXMcwjGmGYTwDoJSqA34DbA4sDwRuE0IIIYQQQgghOuWk0+/0BelaLIQQQgghhBBnJsMwtiilpnXlOU6rarEQQgghhBBCCNHXJJAVQgghhBBCCDGgSCArhBBCCCGEEGJAkUBWCCGEEEIIIcSAIoGsEEIIIYQQQogBRQJZIYQQQgghhBADigSyQgghhBBCCCEGFAlkhRBCCCGEEEIMKBLICiGEEEIIIYQYUCSQFUIIIYQQQggxoEggK4QQQgghhBBiQJFAVgghhBBCCCHEgGI+2R0Mw3gWuAioUkqNC9z2H2Bk4C4xQINSatIxHnsQsAM+wKuUmtZN2y2EEEIIIYQQ4ix10kAWeA54Anih9Qal1NWtlw3DeAywneDx85VSNZ3dQCGEEEIIIYQQoqOTBrJKqXWGYWQe6/8MwzCArwMXdO9mCSGEEEIIIYQQx9bVMbLnAZVKqfzj/L8CPjIMY4thGLef6IkMw7jdMIwcwzByqquru7hZQgghhBBCCCHOVF0NZL8BrDjB/89RSk0BlgB3GYZx/vHuqJR6Wik1TSk1LTExsYubJYQQQgghhBDiTNXpQNYwDDNwBfCf491HKVUWWFcBbwAzOvt6QgghhBBCCCEEdC0jeyGwTylVcqz/NAwj3DCMyNbLwCJgdxdeTwghhBBCCCGEOHkgaxjGCmADMNIwjBLDML4V+K9rOKpbsWEYqYZhvB+4mgx8bhjGDmAT8J5S6sPu23QhhBBCCCGEEGejU6la/I3j3H7zMW4rA5YGLhcCE7u4fUIIIYQQQgghxBG6WuxJCCGEEEIIIYToVRLICiGEEEIIIYQYUCSQFUIIIYQQQggxoEggK4QQQgghhBBiQJFAVgghhBBCCCHEgCKBrBBCCCGEEEKIAUUCWSGEEEIIIYQQA4oEskIIIYQQQgghBhQJZIUQQgghhBBCDCgSyAohhBBCCCGEGFAkkBVCCCGEEEIIMaCcNJA1DONZwzCqDMPY3eG2/zUMo9QwjO2BZelxHrvYMIz9hmEUGIbxs+7ccCGEEEIIIYQQZ6dTycg+Byw+xu2PK6UmBZb3j/5PwzCCgL8AS4AxwDcMwxjTlY0VQgghhBBCCCFOGsgqpdYBdZ147hlAgVKqUCnlBv4NXNqJ5xFCCCGEEEIIIdp0ZYzs3YZh7Ax0PY49xv+nAcUdrpcEbjsmwzBuNwwjxzCMnOrq6i5slhBCCCGEEEKIM1lnA9kngWHAJKAceOwY9zGOcZs63hMqpZ5WSk1TSk1LTEzs5GYJIYQQQgghhDjTdSqQVUpVKqV8Sik/8Hd0N+KjlQAZHa6nA2WdeT0hhBBCCCGEEKJVpwJZwzAGdbh6ObD7GHfbDGQbhpFlGIYVuAZ4uzOvJ4QQQgghhBBCtDKf7A6GYawA5gEJhmGUAPcB8wzDmITuKnwQuCNw31TgGaXUUqWU1zCMu4GVQBDwrFIqt0fehRBCCCGEEEKIs4ah1HGHrfaZadOmqZycnL7eDCGEEEIIIYQQ3cwwjC1KqWldeY6uVC0WQgghhBBCCCF6nQSyQgghhBBCCCEGFAlkhRBCCCGEEEIMKBLICiGEEEIIIYQYUCSQFUIIIYQQQggxoEggK4QQQgghhBBiQJFAVgghhBBCCCHEgCKBrBBCCCGEEEKIAUUCWSGEEEIIIYQQA4oEskIIIYQQQgghBhQJZIUQQgghhBBCDCgnDWQNw3jWMIwqwzB2d7jtUcMw9hmGsdMwjDcMw4g5zmMPGoaxyzCM7YZh5HTnhgshhBBCCCGEODudSkb2OWDxUbetAsYppSYAecDPT/D4+UqpSUqpaZ3bRCGEEEIIIYQQot1JA1ml1Dqg7qjbPlJKeQNXvwTSe2DbhBBCCCGEEEKIr+iOMbK3AB8c5/8U8JFhGFsMw7i9G15LCCGEEEIIIcRZztyVBxuG8UvAC7x0nLvMUUqVGYaRBKwyDGNfIMN7rOe6HbgdYPDgwV3ZLCHOfH4fOG3gbNBrrwuUAuU/cqH1NqWXIxx9vZUBhgkM9BoDjNbbTGANh5BoCImB4CgI6tLXiBBCCCGEEKet02eghmHcBFwELFDqK2fIACilygLrKsMw3gBmAMcMZJVSTwNPA0ybNu14Z9hCnPl8XqjJg/IderEVtwetjsDa1djXW9nOGhkIbANLaAwkjoSUCTBoIsRmgUkKpAshhBBCiO7TqUDWMIzFwE+BuUqpluPcJxwwKaXsgcuLgAc6vaXi7NNSBxW7oHI3VO8HdxN4nOBpAW9g7XGCxwFeh85Sxg6BuGEQPzywDNXXQ49ZWLvveV1Qtac9aC3fAZW5+v0BWMIgNlNnP6PSIXlcezY0NKY9eDSHtGdM25YOWdTWLOvRjKOuq8A/rVncoy/7feBuDgTWHZeG9sv1ByF/Ffg9+jmDoyBlvA5qW5f47IGTyfW6wF4BPjdEDoLgiL7eIiGEEEKIs95JzyQNw1gBzAMSDMMoAe5DVykORncXBvhSKfVtwzBSgWeUUkuBZOCNwP+bgZeVUh+e0la57FC2DUJj9RIcpU/KxZnJ74O6QqjYCRW7deBasRvsZe33CUvQAZslDCwhYAmF4BS9bl1AB1Elm2H36xzRdTY8UQe0SaNhzCWQeX7fBVKuJtj7Dmx/CQ5vAH+gblpwlA7ypt/aIeAbDqagvtnOrvC6oGqv/pu2Bug5/9QNDqD/jqOWwaRrIWtu377Hpioo3qiDVXt5h3WlXjvqjrx/SLRuVIhOg6i0wDpwPWEkRCb3zfsQQgghTsTnBUe9/l1z1HdIEAQSAp7A0pos8LrA59EN036fPl9pWzpcP2JYk2of6oQ68nYI3Hacyx0du7OnGEiWP6N76PUg4zi9gvvUtNQglXN7h6yHEaSzT6FxOrANi4OYIRA3tH2JGQxma99ttDg9jnrI+wj2vwcFH+svUwCTGRJG6Axe8jhIGQfJ4yEi8fSe3+PUQW3dAagtCCyBYNnVCOFJMO4KGH8VpE3t+YYSvx8OrYcdKyD3TfA06y63Yy6B1Mk6aI3JPLO74Pq8+u9QvkMH8Llv6ExuZCpMvEYHtQnZvbMtjWW6MWHPW3DoC9p+QI0giEyBiGSdfY1MCayTdda7sRRspYF1iV6ODnQHTYIRi2HE1/TlM/lvKoQQou/5/brxv64Qag/odXM1tNTq3m2OOn3ZaTv15wwK1r97QRZ9bmYy6wSAqeMSpH83TUFH9v5q6xVmHFlro7Ub2AkvdyRJrAFt8e90jHYchmFs6er0rP0zkJ0wWuW88lig1ajjEmhBaq6B+kPgtrc/yDBBdEYgsM3SmayMWTpAGChdGM90thLY9z7se1cHdX4vRKTAyMWQPkMHrYmjwBzcc9vgcUL+R7DrVchbCT6XDijHX6WXxBHd+3p1hbDj3zqAbTisx5OOuxwmXQcZM8/ungZeF+z/ALa/DAWrQfkgfboOaMde0f3dwRsOw563dfBasknfljgaxlwK2Qt1Y1hY/Olnh90tOjBuLIHSrfq4KtmkW6DDk2DEIh3YDp0v3ZKFEEJ0ntcN5dt1j6e6A+1Ba11Re48ngCCrbpANi9NJoLA4/ft2xOUYfU5iCdG9pMwh7T3czCEDszeYGFDO3EB22jSVk5Nz4jsppQPausIjl/oivXbU6/tZI2HwLMg8FzLPk8C2t1Xt05mvfe/qL1/QGddRy2DURZA6pe8yVk6b3rZdr0LROh14DJoIE66ByddDSFTnntfnhT1vQs6zOmDHgKHzdPA6ahlYw7rxTZwh7BWw8xUd1Fbv1S3Boy+Cc+7SGfNOP2+lbkTY8xaUbdW3pYzXwevoS7u/4aJVc60OzvM+1D0OXDZ9YpF5Loy9HMZ/XZ88CCGEEMfjdevfroOfwcHP4fDG9oA1yKob4uOGQvywDuthetiL9AYS/dzZHciejL1CBxEHP4eD66Fmv77dGtEe2GbN1d06z+asWE/w+6FgFXzx//SXL+hM26hlMHJZzwUPXWGvgN3/hV2v6PHZwVEw7RaY+W2IGnRqz+FxwvZ/wfo/Q8Mh/WMy+TodGEen9ez2nymU0g0e21+Gnf/RjQ1D58G5P4Cs80/9s1qZCxv+ohspfG7dYDLmUt2V+wTdXHqEzwOHv9RB7f4PdCt6eCLMuAOmf0u3jgshhBAnClyTx+lz1yFzdKN7dLpkTcWAJoHs6Wiq0l8KrcFt9T59e+JomHozTLxaj78Vnedx6uBjw190w0FkKsy8Q49/jEzp6607dWXbdDC65009BmTC1TD73uMH4E4bbP4HfPkkNFdB2jQ47wcwYom0iHaFy66z2hv+Ak2VJ9+vSsGBT2DDE3ptDtUNCbPu1K3U/YFS+gTli/+nu7ibQ3X2/5w7ez/AFkII0T9U7IKtL+ieSc4GfVtr4NoavEqjpzjDSCDbFU3VsP992Po8lG7R4wHGXKaD2sGzJEt7OpprdCC3+e+6uEDKeDjnHt2FciAX4Kor0kHUtn/pFtGRy1iXdC0/3xxGWYODsdEu/jh4PSMOv6ILSA27QGcOM8+V46c7eZy6wvP6P+lMd+JocjJu5ge5wyi2uRkSbeaPY/KZUvoyVOXqcUEzbtcZ9f78w1+1F754Qjf+KB+Mvlg3mKR36TtdCCHEQOC0wa7XYNuLugE9yKp/B8Zcps8j+vPvlxDdQALZ7lK+Uwe0OwMBScLIQJb2GvkiOZGaAp392rFCl2rPXgTn3H16XUAHguYa2PR33F88idVjY7N/BAX+NK4I+hwLXsrTFpF20S8gdVJfb+mZzeeF3P/SuOphouwFHPYnsto/lYuCviTJaMAWNYLo+d+D8Vf2bMGw7mavgI1PQc4/9InN4Nlw/g9h+IV9vWW9x92sp0Aq2667YgPt0zIcNT1DUDCkTYGMGRAc2Ucb3M+4W/S0Y4e/bM/mHJcBCcNhyLm6SviZ9F3dWX6/7qV1aL2usXEyYXH6c5o2Vca6t2qp09Xoy7bp4RxtOhxfrcda62d48Cw9ndnZRCm9n7a+qCv3ex2QNBam3AgTvi7nnOKsIoFsd3M36y+WLc/pk4KgYD1Fy/xf6IqmQnPa4NPfwaa/6663E6+GWXdB0qi+3rIua3Z5Kbc5sDu9NLm8NDm92APrv63eyWLPx9wW9D7JRh2v+87nKd9F1AZn8N0F2YRagwi3mgmzBhFmNRMTZiE9NpToUAuGnCyeMqfHR0m9g3Kbg2aXD4fHi8Ptp8Xtxenx8dTaAmZ6NnOX+S0mmwpY5xvP333L2Gmdwq8vHkt0qIXoMIteh1qICbMQbB4A44hcdp393/BXsB3WlY6/9rv+0y26O7nseuzXwc908FC2rX0+5eNq/Qx1mCpp0EQYMlsvg885e04CnY068D+0XteAKNum53nEOHlw7/fq+SFBz889ZHag6+JsfUJ9NgyH8Hmhcpeeeuvgejj8RXuBSEv4yccduuzoRhWrHvKQOUfvv/QZZ09l8sZyvd8OBZaqPfp2I0hP1wJHzQPa4bLPo68bJt2Da8icwDL7zP0Me5z63HLz3/U0dNZI3eg65QZdw0HOEcRZSALZnlSxW3/pbPuXrmY7+2449/tndwZAKT2VzKpf6y7E026BeT+DiKS+3rLT4vL6OFzbQmFNMwdrmjlY20xhtV5XNrpO+FjD3Ehw4kqCrNU4ay/A1zSSk81zFhlsJi02lPTYMNJjQ0mPDSUjTl8eHBdGZIilG99d/+f2+jlc10JxfQsl9Q5KAuvSegcl9Q5qmk78NwDA5MIatxZrxF5cDbPwNEwHjn0CbhiQERvGiOQIspMj9TopkuFJEYRY+mGA63XDxr/B2od1ZuOcu+C8Hw3sE2SfFwo/1dXBD36u5xJWPt0QljolEAicCxnTdUG+jnMKHn2C57JD8Sad1Tj0BZTk6Gm0QNc8GDIbhs6FkUvbT6gHOqWgcA3kr9LBa8VO/bvUuv9ag9GMGSfPcCmls46H1rcHcrbD+v9ConWDwJA5uhE3Or3H31qvqT+oC/od+kJnrlun74vN0u+3NRiNGXLyoKKl7siGhI7H86BJ+nlGLNbrMyVA8Th0lf+itXoftmaurRF6Krkhs/V+TJty8h4xHodOFhxcr/dhyWbdqwsgaYx+nqFz9T4c6J9hv0+fN635PdiK9b6aerMuPmgN7+utE6JPSSDbG2yl8PH9ehxbRDIs+DVMvPbsaLXuqGIXvP9jffKYNg2W/UFXfO7nGp0edpXY2F7cwM6SBnLLGiltcBzRUBwfbiUzIZyswJIeG0pUqIXIYDMRIWYigs1YzT6WPPcgrohVYPhQnkhM1ga8TcOJarmCVfdcg8Pto9nlpcXto8Xto67ZFQjU2oO14roWmt2+I7YxJSqE4UkRX1niw60DOpPr9Pg4UN1EQZVe8iubKKhu4mBNM15/+x/AEmSQFnNkoJ8eG8ag6BAiQsyEWc2EWoIItQRhMSvOf/IRHBHvYTI34XfHY7LW4nOmEmG/kle/eS02h+eIparRRUFVE3mVdoo6vLbJgMFxYWQnRzJ6UBTTM2OZOiSWMGs/mZ7LXgGr74cdL0PkILjwft31bCAdE+7mQJb5L3p88xEZrDk68OrqyZzHqTOSrYFZ8SYdpESl6arjU28auN0XvW7Y/bouDlaVq2s5pE9vz0KnT++ek+GGw3BoAxz6XO/D2gIdlI1bDrPv0Vmzgap0iy7et/dtHfwnjmoPuobMhqjUrr+Gyx4IbAPZydItgWrpk/X+G33pwJ32r7kWNj8Dm56GlhpdFHPw7PZjMGVC19+b16Xn4D4UmGWieBN4mvVneNZ3YMpNnZ8Or68opeuwfPyA7raeOhku/F9dhV8IAUgg27tKcuDDn+mWw5QJsPj3ugX8TOdo0C2Jm57WP2AX/i9Mur5fBvJOj4+95Y3sKG5gR4mNHSUNFFY3t/1/VkI449KiGZoQztDEcDLjw8lMCCc69PgtvkopVh5ayeM5j1PWXIa/aRwtFYtRnhgssRsJTlyNKcjJ5dmXc/eku0kMSzzhNiqlaGjx6KC2voWDtc0UBAK8gqomWjoEuTFhFoYnRjAsMSIQaIeRmRDOkLhwQq39I5Po9yvKG50crGluy3AX1TRTUNVEcX1LW4OByYDM+HCGJUWQnaTf05D4MNJjw0iKDMZkOnFwppTi89LP+eOWP1LQUIByZNJSsRS/MwNz1E5Ckt7HsNhYkrmEH0z7ASnhx66S7fH5OVjTTF6lDmzzq+zkVTZRWN2EX4HZZDAuLZqZWXHMyIpj2pA4osP6OCNQkqMbkcq26tb8JQ/3/0akpmr9nbH577rLZvoMfUKfvRAsoT372j6vrlr9xZ9112VrpA5mZ31n4GQYnTbI+afOzNvLdZZq9j0w9oreGZNZf1CP297yvA4ohl2gX3/o/IHRkOL366rgX/xZN3AER8O0b+oicL0xFZq7JVDB/wndKBAzWA+/mXz9wOlZUXsAvvwrbHtJj+McsVjXwBgyp+d//32ewGc4MIVfcJT++838dvc0PPS0g+th9f9CySaIHw4X/EpnYAfCZ0eIXiSBbG9TSreOr7oPGktg9CWw8AGIy+rrLet+bd2If6WLHU27BS74n341fqXZ5SXnUD0bC2vZWFTHzpIGPD59PCdGBjMxPYZJGdFMzIhhQlrMaQckubW5PLLpEbZWbWVk7Eh+Mv0nlFWk8ejK/ZQ1OEiNCeWuBamU8DYr9q3AYrLwrXHf4saxNxJqPv2TdaUU5TYn+VXtWcwDVU0U1jRR0+Q+4r4pUSFkJoSRlaAD8tSYUOLDrcRHBBMXbiU2zII5qGsnG0opmlxe6prd1DS5qW1yUdPk5lBdoEt2jQ7EXV5/22NCLCYy48PbssqtXXgzE8I6PU51f91+Hst5jA3lGxgcOZjvT/0+9tpR/OGjvLa/w3cXDqHK9AHP5T6HyTBxy7hbuHnszYSYT+2kv8nlZeuhejYW1bKpqI4dxTbcPj+GAaNSopiZFcesoXHMHp5AVF90Bff7dWZ29f/qz+OUG2DBfRCe0PvbciK1B/TJZ2sBuJHLYM69uqhLXyjbpitD576hTyLHXqGHiQya2DfbczK2Ej2N15bndVY5a66uZD18Qd+cBDvq2wPqpkpIHq8D2nFX9M8unx6nngv8iyf0FHBR6Xpqqyk39s2wIL8f8j7QGeHiLyEkRs8dPeMOiEzu/e05FcWbdAPA3nf133jC1TqA7asaGKVb9XfKnjf1+NvxV+ljMHlM32zPiVTs1j348j/SvWjm/Uw3/A/UbLwQPUwC2b7icegfys//qAtnnHMXzPv5wKqUeiINh+GNb+uW7H7UjbjR6SHnYB0bC+v4sqiO3aU2fH6F2WQwPj2aGZlxTB4cw8SMGFKiQjrdLbeqpYo/b/0zbx94m9iQWO6dfC+XDb+MoBMUADnceJjHtzzO6sOrSQ5L5rtTvsuyocswGd3Tct3o9HCopoWi2tYgsrntcn2L5yv3NwyICbUQFwhu48OtBJuP3JajP/lK6depbQ1am924OwSprSxBBoPjwshKiGjLErd2y06ODDlpdvVUVbVU8cS2J3iz4E2igqP49oRvc/XIq7Gc4AS6tKmUx3IeY9WhVaRFpPHDaT/kwsEXnvax4PT42F7cwKaiOjYV1bHlUD0Ojw+zyWDKkFjmjUxk7ohExgyK6t3u304brH1EBxYh0bD0UR2c9XVLf0kOrP+/9pPfidfoKbiON/dyb2s4DF/+TVendzfpAPHc7+lMY39QuQc+fxxy/6s/iGMv1yfr/aUSutcFu17VAUX1vvYun9Nu6R/j/JyNugfAxqf0XN4p42H2d2HsZf0n4C7epPff3nfaA8Q539MVpPuaUrD/A/0ZLt7YPwPu+oO6EN62F3WxsuEX6v2XdV5fb5nufr3qV7D9Zd0F+twf6Oy/Nayvt0yIfq3XAlnDMJ4FLgKqlFLjArfFAf8BMoGDwNeVUvXHeOxNwP8Erv5WKfX8yV6v3weyrRrL9fiHHS/riauX/2PgV+7NfQPe/q4eS/S1B2HyDX3Wjdjl9bHlUD2f5dfweX4NuWU2/EoHUpMyYpiZFc/MoXFMGRxLeHDXWzy9fi8v7X2Jv2z/C16/l+vHXM/t428nwnrqXcG2VG7h0c2Pkluby/iE8fxq1q8YHT+6y9t2IrYWD5V2pw5Am11tGdS6Iy678fi+GpQeHf5EtQa/4cEkRFiJj9CX4yKsJIQHEx9hJSkyuMvZ3hPx+D38a8+/eHLHk3j9Xq4ddS23TbiN6OBTH+e4qXwTD21+iPz6fGamzOTX5/yawVGdrzzu8fnZdriBtXlVrNlfTW5ZIwBJkcHMHZHIvJFJnDs8ofe6IVfthTfv1N2NR18My/7YN0XXbCXw0f/o742QaJh+a/86+T2aowG2/FMHtU0VkP01PUykrypDO+p1BfjNz4AlTI8FnPXt/lsl3++HgtXt3baj0vXvRF91m1RKT5u36lc6Yzz8Qt0AkDW37xt3jqdjl12/B2bdCXN/0neFJKv26qELBz/Tx905d8Ok6/pvF+iWOj1V2candNHJURfpYzA2s/e3pfX4W/lz3cg460447wd6GJYQ4qR6M5A9H2gCXugQyD4C1CmlHjIM42dArFLqp0c9Lg7IAaahE0BbgKnHCng7GjCBbKv9H8Jbd+mW/kW/1Sdz/fVH9HjcLXoM8Nbn9dx4y//R612mlVIcqG7ms/xqPsuv4cvCWlrcgSzY4FhmDYtnViBw7e5qs/vr9nPfF/eRW5vLvPR5/GT6T8iIyujUc/mVn/cK3+OxnMeod9Vz/ejruWvSXYRZpHX2ZHZU7+CBDQ+QV5/H/Iz5/Hj6j8mI7Nzfwev38lrea/x565/x+D18d8p3uXb0td2SJa9qdLI2r5q1efpYtTk8BJkMpgyO4YJRySwYnUR2UkTPZmt9Xj0G79Pf6azY0kd1cZ7e+O7xunR26bPHdKPXud/XJ8D99eT3aF6XPhFurQw9+159Atpb2UW/H7b/S3cVd9TrzOb8X/aroRsndegLHQBV7taB49JHIXFk771++U744CeBAoRTYcmjkD61916/q5qqdDfUbf+CiBR97jD+yt47d3DaYM3DundHcCQs+BVMuXngdIP1OGHjk7D2UV0xes73dC+Lnh6D36quEN79ga7Gnj4dLv4TJI/tnde2jAmfAAAgAElEQVQW4gzRq12LDcPIBN7tEMjuB+YppcoNwxgErFFKjTzqMd8I3OeOwPWnAvdbcaLXGnCBLOgfpTfvhIJVupX/0r9AxIkL//QbFbvhtVugJk//EMz/Za91x7I7PXyWX8O6QEBQ2uAAdGGm87ITOD87kVnD4onohozrsbh8Lp7a8RT/3P1PooKj+PnMn/O1IV/rlgCk0d3In7b8iVfyXiElPIVfzvwl8zLmdX2jz0CN7kb+vPXPvLL/FZLCkvjFzF9wweDu6fZZ2VzJA18+wLqSdUxOmswDsx8gMzqzW54bwOvzs6OkgTX7q/l0fxW7S3W2Nj02lAWjkrhgdDIzs+J6bqqf6v36u6c0R2cnlv2xZzOieSt1o1ddYSAb8juIHdJzr9eTGsth9X26ME9vZRdLt+gAsHSLLt619NH+O2b3ZHxeyHkWPv2trlA989sw96c9W2G2pU433uT8I1CA8H6dQeyHBQhPSfFmeP9HUL5dF1Ja8gikjOu51/P79fHeOo3e1Jt1MaLw+J57zZ5kK9UZ+d2v64zy134Po5b13GfY59ENiGseApMFLrxPN0SdbO7hfkophdvnx+X14/T4cHn8uLw+nIG1y+PH61d4/X68PoXPr/D6O671//uVfi6lwK+OvK5ovR54TdQRM0e0xiFHhyP9b+CjOF3XTM8gKer4tUr6OpBtUErFdPj/eqVU7FGP+REQopT6beD6rwCHUuoPx3j+24HbAYJThk+dcPeTxIRZiA5tXaxEh1qICdNLanT7XJzd0a20Wyilx+l89Cv9Q37Zk7pKZ3+llO7StvKXEBoDlz8Fw+b3+MtW2Z2s3lPFR3sq+KKgFrfPT2SImTnDEjhvhA5eM+J6Pnu5tXIr931xHwcbD3LJsEv48bQfExMSc/IHnqbtVdu5f8P9FDQUsGDwAn4242fHrap7tlFKsfLgSh7e/DB1zjquG30dd026i3BL92bGlFK8W/guv9/0e9w+N/dMvofrR19/wnHPnVVhc/Lp/io+3lvF5wXVOD1+wqxBzBmewIJRScwflUTyCb7YO8Xv0ydXnzyox2UteUQXRenOk7m6Qvjw55D3IcRn6+rJwxd03/P3paOzi0se6f5hIs01OgO39UXdDXzhA3qc5EDrvXMsTdUdsotJsPA33T9VlN+vx0d+fL/OYk+/Deb//Mzoxun3wdYX9FAlpw1m3KbrboR28+9R+Q59nBdv7Ff1L7pF0Wc6Q1+1B4Yt0N9PCdnd+xolW+Cde/X3xKiL9PdEb1TBPgWtMyJU2p1U2JxUNbpocLhpdHhpdHpodHhodHoDa0/b7Q6P7ysBpBDd5d17zmVc2vGHhQ2EQPbHQPBRgWyLUuqxE71WWvY4df3vX26bB7LBoT+ENoeHJpf3K/ePC7eSHhtKRus8lHFhDI4LY1xqFPERfVCAqXIPvH6rnvdvxh2w8P7e6+5yqppr4e279Txn2Yt00N2DFVALq5v4aE8lH+VWsK24AaX0HJ6LxiSzcEwyU4fE9ui4y46a3E3839b/4z/7/0NaRBq/nvVrZqfN7tHX9Pg9vJD7An/b8TdMhol7p9zLNSOv6ZFAaqAothfz4MYHWV+6nrHxY/n1Ob9mTHzPVqKsbqnmgS8fYE3xGiYkTuA3c37D0OihPfZ6To+PDYW1fLK3io/3VlJmcwIwelBUW8GoqUNisXTXsV+dB2/dqacJG7lUZ2ejBnXtOd3N8Nkf9bjIIKvOuM38Npit3bPN/YXPq8fPfvKb7s0u9tTz9kclW3R2sWwrDD5Hn+gPmtDNzztbZ7F7MmvZV1rq4NMHdZY7NE6fO3THvPUtdfDJb/XzhsXrRpSJ3xi4Wezj8Xl14/ynD+qinOfcCef/uOvjj112vf82PgWRKfr4G31x92zzKXJ7/Ryua9HFHmubKWtwUtkYWOxOKhtdxyzOaDJ07YuoEAtRoWa9DlyODLEQZg0ixBJEsNlEsCWIkA7rEEsQVrMJS5AJs8kgyGRgDjICl9tvCzIZmAwDwwCTYWAywMDAMOk6HK3/Z2Ac0bbVernj7Uc3ffVqMUXR7UzGif+GfR3I9knXYo/PT32Lm9J6R9tcnCX1DorrWtpuc3coapMeG8rE9BgmpOtpWMalRfdYN9UjN9SpW46//Cskjoblz/SfH96iz+C/t+kMwcIHdPXJHviyyKu089b2UlbmVlJQ1QTAuLQoFo1JYdHYZEYmR/b6l9S6knU8sOEBqlqquH7M9dw96e5eHbvaMXgbEz+GX836FeMS+slx0UtcPhcv7nmRv+34G2aTmXsm39OrQb1SiveL3uf3m36Pw+Pg7sl3c+OYG3v89ZVS7Kuws2Z/NWv2V7HlUD1evyIy2Myc4QnMHZnIvJGJDIruYqOX36e/dz75LZjMupDMzO+cfuDp9+kpuD59EBpLdfbwwvu7Hhj3d801OjO29QUIT9RZv8k3dm7sYMHHuodOVS5kna/HcQ70goAnc/T438k36OEqnenubiuBj38DO//dN+NI+0rHzGnqZFj0IGTOOf3n8bp1F+y1D+vKzjNu11PCdHemt79pqoLV9+vjMCIF5v9Cz+F7ut/xfr+eSuzj+/VzTv8WLPi1LmzXA5RSlNQ7KKxppqi6iYO1LW3zs5fUt+DvcLoeagkiJTqE5KhgkqNCSIkKISlKX0+JCiEpMoS4CCvh1iAJBkW/1deB7KNAbYdiT3FKqZ8c9Zg4dIGnKYGbtqKLPdWd6LW6MkbW71dUN7korG5mV2kDO4pt7ChpoKTeEdgmGJ4YwcSMGCZlxDB3RA93Yy1YrcevORp0N57JN/Tdj7Dfr6cM+vRBiBuqCzp18/QOlY1O3t5exhvbStlT3kiQyWBmVpzOvI5NIS2mbzLThxsP81jOY3xS/AnDY4Zz/+z7mZDYDZmCTlBKsfLQSh7e9DA1jhoWDlnInRPvZHhsP5iGoQd5/B7eLHiTp3Y8RWVLJQuHLOSn039KcnjfVLitcdTw2y9/y8eHP2Zc/Dh+OeuXvdqoYHd6WF9Qq4tG7a9qy9aOTI7k3OwEZmTFMT0zjrjwTmY+aw/oYQN5H0D8cFj80KkNdVAK8lfpQKQqV59If+13MKRney30O6Vb9P47vAESRsKi3+jeK6fy/V2xW4/bO/AJxAzRDYZ9Vdm3rzjqdTGhzX+HoGCY8109h++pFNRyNurpiL78qz4eZ30Hzv9R31X27QtK6bGsHz+gG5JGLtMZ2lPpLqsU7HlLf4bri2DoPP0ZPtuKERVvhpW/gJJNOqGw8P5T/wwf2qBrAZRv192wFz8EGdO7bdNsLR72VTSyv9LOvgo7+8obyatsOqLXYbg1iKxEPVd865zxWYnhZMWHExNmkQBVDHi9WbV4BTAPSAAqgfuAN4FXgMHAYeAqpVSdYRjTgG8rpW4NPPYW4BeBp3pQKfXPk71eTxR7qm1ysbNEB7U7ihvYWWKjttkN6BPH+aOSWDA6ickZMd3fxbWpWmdACz/V3YSWPdb784s56vXcsHkf6nknL/lzt50U2J0ePtxdwZvbS/niQC1KwcSMGC6flMpFE1NJ6Ivu3QE2l42ndj7Fin0rsJqs3DbhNm4ac9MJ5yLtLXa3nRf3vMgLe16gxdPCsqHL+M7E73Rpmpj+yOf38X7R+/x1+18paSphYuJE7p18LzMGzejrTWsbo/vQpoeoc9Zxefbl3Dv5XuJDe7fwiVKK/Kom1u6vZk1eFTkH63EFuomNSI5gRlYcM7LimZEZR0r0aY6vzV+lT8hqC2DEYn1Ce7zpZkq3wKr79FQcsVk6+zD28h4PwLw+f9tQEn9bgZCvFgdRSveGjA61EBNqJcRi6tmTOaVg33u6ME7dAZ1VXfTb4xdnaizXhY+2vaSzNnN/oqvY9+Ac40opGp1eGlrcNLR48Pq/2r3waOHBZmJCrcSEWXquCFmr2gM6o7XnLZ0du+B/YNK1x86O+Tyw5TldSKelRvcCuOB/enw6IqfHR32Lm/pmDw7PV4cvHc0SZCI2zEp0mIXIYHPPHoMehw7oP3tcz5867RadVT3eUKDiTboBpmQTJI3R45WHL+jRz7Dfr7C7vNha9PhLf4fzSqNDZ9HWTTAZBlGhZuLCrYRaejhjqBTsfVsH9XWFkHmeblhKm3Ls+9cf0gXgct+AyFQd/I67stPdsF1eHweqmtlf2ci+Cjv7K+zsK7dT0ehsu090qIWRKZGMSolkZEokwxMjyEoMJzEiWIJVcUbr1Yxsb+qNqsVKKQprmvl0XxWf7KtiU1EdXr8iJszC3BGJXDAqibkjEokJ66ZxYH4frH1Ed/FJGgNff6H3JkIv2w6v3AiNZboq54zbu/yj5vX5WZtXzRvbSlm1pxKX18+Q+DAunZTGZZNSGZrYt9NwePweXt3/Kk/ueBKby8YV2Vdw9+S7SQjtuXHAndXgbODZ3GdZsXcFHr+Hy4Zfxh0T7mBQxMDuwulXfj4+/DF/2fYXDtgOMDpuNHdPvpvz0s7rdz/OTe4mnt75NC/ufZGQoBDunHQn14y6Boupbxo8XF4fu0ttbCyqY1NRHTkH69ta6ofEhzEjM44JGTGMSIpgRHIksSfL2nrdgakqHtHTzZxzN+9Ef4OHPi6hrMHB9KgGHk94i7SylRCWoE+Up9zUpXGwbq+fykYnpQ0Oym0Oym1O6prc1Ld4dNDQ4qa+WV+3OTydeg2r2URMaxHAUB1YxATmQk6OCiE1JoSU6FAGRYeQEBFMkKmTx53PAzn/hDW/142CE69hZfKtPLCukbIGB8Oi4YnMzxhV+Ly+78w74Lwfdno6Hb9fUdPsosLmpKzBSbnNQYXNSU2TWwesDr0PGwL7zufv/O94iEUHZTFhVmJCLcSGW4gJ0/NFp0aHMigmhEHRoaTGhBBm7cKwnMMb4aNf6vHbSWNZP+x7/GRbAmUNDlKjQ3hsYimzDvxJN7hknqcz4F0oRNTk8lLe4KDM5qTC5qCswUmV3UlDS/u+a73sOsbYwlMVZDLajkG9H/X+i4+w6v0XHUJqTCgp0SHEh1s7/93XVA1rH9LHoSUMzvsBb4deysOrD331MxyRHGgwuK7T1XS9Pj+VdlfbvquwOSm3OY84/mwtuoZJQ4ubzh6CVrOJ2A77LjbMSmy4lbgwKynR+jM8KDqU1OhQokK70GjQ+hle+xC01MK4K/ko5Xbu/7yZsgYHw6PhicFrGFn4PBgmPYPD7HtOeUqu1m7B+yvs7K+0s7e8kf0VdopqmvEGdo4lyGBYYgSjB0UxMhC0jk6JIjlKAlZxdpJAths1Oj18llfDJ/uqWLO/itpmNyYDpmXGcemkVJaNH9Q9QW3Bx7oQlM8Dl/4/ne3oKUrpcV7v/1i33l71fJe7xhyqbeaVnGJezSmhyu4iNszCxRNTuWxyGpMzYvr8y1gpxWeln/GHnD9QZCtiZspMfjz9x4yM68X5DTupuqWaZ3Y9w6t5rwJw1YiruG3Cbf0y+D6R1r/BE9ueYG/dXoZGD+WuSXdx4ZALu2UO155UZCvi4c0Ps750PUOjh/LTGT9ldmrfd6n1+vzsq7AHAttaNhXVUd/SHvwlRgYzIlkHtXqJIDs5kqiQowJxe4XOTOxYQaWK5Q/eqxhnFHFt0Ce4MVM86hZGXfHLE/bWUEpnX6oaXVTbXVTZnVTbXW0BV1kgeKhpcn2lGmaoJYi4cH3CqtdW4gIn/nHhVqJCzQSZTBjwleIgrbf5/Pr7uqHFQ4MjcDIduNwanNQ1u4+olQBgNhl6LFl0CIOiW8eUBZMUeeT6hBk2pw0++yO+DX/F41P8w7eEUpXI982vkWjYKEldTPqVD51wDu4Wt953lY1OquwuvTQ6qWh0Ut7gpLxRB60e35E7zxpkIiEiEHAeETRZjghEreYTf8YU0OzydgjmdGNCa0a3PnC9LtBjqaOoEDOpMTo4GxQTSlKkHp+XFKn3X3JUMPEnajBQCva8SfN7vyK8pZh1vvG85FvAN80rmWXaiz1iKJEX/073HDjO38Dr81PT5KbKriuzVgbWVfb2oL/c5sTuPDKzahgQFwiSYsP0TAixYRZiA8djTOB6WLD5KwVnjuby+o/YX63BXH1ze2BX0+T66t/QbNL7LloHZynR7fsuKSqY5MD6hFny6jydMdz/PmUqnse9yxlpFHNj0Ed4MHN41K2MvuIXx53TWSmFzeEJHHd6v7Vermh0tO3DarvrK8Fp6+c3Nry94Sg2sO9aA/ioEHPb3/+IKVY6PI/Pr3tg1Ld4Ao1Z7cdgx/XRDTRhVj0utLWBoPUYTIoMJilwHCZGBp+4eJ7TBuv/hG/9E/h8Pp73LaJApfFD86skGQ0Up19ExlXHrkbscPsorm/hcG0Lh+v0Uty6rm/B6TmyPktrhnVkShSjUiLJSgjvvsJ+QpwBJJDtIT6/YmdJA5/sq+KD3RUUVDVhCTKYNzKJyyalsWB0Ute6Y9lK4NVv6q4/M7+tu/50dwVQd4uu8rj9JRg6Xxeb6mRVYqfHx0d7Kvn3psN8caAWkwHzRibx9WkZLBid1G++mPPr83l086NsKN/AkKgh/HDqD5mXMa/Pg+vTVd5UzlM7n+LNgjexmCxcnn05Fw29iPEJ4/v1e2l0N7L60Gpez3+dndU7SY9I585Jd7I0a+mAqs6slGJdyToe3vwwxfZiLsi4gB9N/xEZkRl9vWltlFKU25zkVdoDSxP5gbXD42u7X2SwOXCyaW2bziw2zErRtk/5ifoHE0xFeJWJf/vm8yfvclwhCdwxdxgurx93YHF5fdQ1u6mytweuHU/YWoVZg9qCnNToUH05JqRDVq+LGb3T3D91zW7KbYFMUqOT8gZHW1apolGfrB/rfYRYTDqwiAwmLtxKmDWI0EBlz9DA8u5nG7nDt4Irgj4HYLN/BL/zXMeB4NHceE4mDo9PL269NDgC+6/Rhf0YlfetQSaSo4MZ1HqC3rbupmxeJ7i8PiptLsoC2eAym0MH2rb2YPFYwa7JgPiI4LYAIzzYrPebVe+7EEsQ/1qfx2WeD7jH/AYxRjM1KorHvVfyvmUhN8wZjrN13wX2Y5PTGzj2XNQ2f7WBBCAhwtoWHKYGgpyO+zI5KuSkQX53OjqrXhHYZ0dniI8OdgEiQ8xtDQTRoRa936ztx1+oNYita9/hu/7nmWAqwqcM/uObx+PeK3GGJHL9rCE43D69Hz0+Wtw+aptcbQ0nx6pu2zFIbN2HKR0+u4OiQ4kK6eFu1B34/Eo3kLUdd+1Bdus+rLIf+1iIC7e2BbXx4VZCrUEEm4MItpgINgcRYjHxxprN3OZdwZVB6zAZim3+4TzguYF862gum5xKk9OL3enF7vLq469Jf/91FG4NYnB8OIPjQhkcF0ZWQkRbprVXiooKMcBJINsLlFLkljXy1vZS3tpeRpXdRWSwmcXjUrhschqzhsZ3rrua161bVb/8qy4kcNVzENNNJ8q1B3RX4srdepqHuT/tVPeifRWN/GdzMW9sK6WhxUN6bChXT8vgymnpXa+s2o1ya3J5Lvc5Pjr0EeGWcL4z8TtcM/KafjEOtisONx7mbzv+xsqDK3H73WREZrA0aylLhy7t0SljTofL52JdyTreK3yPdSXr8Pg9DI4czE1jb+Ly7Mv7rGtud3D73Lyw5wWe3vk0Pr+PG8bcwA1jbuj18bOnw+9XlDY42oLbKrsTW4eska1Dxi3IUs2wmA9pcmVQ1ng+cORJvjXIhCXIINgSRGyYpUPmUp8gtgZ7SVHBJEaG9OpJbnfomFluzSpXNrZm9/Rt9c0enN72oMrp8bUHHoaHjOjVBJttHKhfhvLpLHaQyWgL2EKtJkItQUSFWEiOCtH7rUPmrWOwMpD2XSu31091k84otwZJ1Y16KpAqu5PqJhctbh/ODkFpx8aDyNA9pEV+SUnjeTQ524sYHRn4mgizmo841o7IBEcFkxBxkixcP6WUor7F055dDuzHtmPR7sIemOvT4fbj9PhocXvbMqUmcwPDYt6nxZNIqe0CQP/OW4NMhFhMRzQetAZ3SW37LeSIbOZADLyOzs63fm5bM8zVdid1LW6cHj8ujw9noHGunY/06E8JtVZyoH4pfq+eQTIu3EpkiJmIYL1EhpiJDbMyJD6MjMD0jkPiw4mVgktCdIkEsr3M51d8WVjLG9tK+XB3BU0uL8lRwVw8QXetHZsadfpfanvegjfv0lM7XPY3GLm4axuZ+ya8fY8OXK/4+6lVKe2gocXNOzvKeG1rKTuKG7AGmVg0Nplrpg9m9rB4TJ0dY9bN/MrP56Wf81zuc2yu2EyEJYKrRlzFLeNuISbkzJpawO62s/rQat4reo9N5ZtQKEbHjWbZ0GUszlzc61V/fX4fOZU5vFf4HqsPrcbusRMfEs+SrCUsG7qMsfFjz6gf98rmSh7f+jjvFb6H1WTlkuGXcNOYm8iMzuzrTTttSim2Vm3l1jcewxOyG8PQ3/9+VwKu2vkkMIu1P16ANaiHiygNYA1OOwuefghn+KeYzHYAlN+Cp2E6se6FbPjJctl3J+Dz+Zn1p7/QFPoR5rCDbbd7m4cT6VjEl9+/A9OZNr9pN1JKUdRQzPIVD+IJ3Yhh0j0w/J4Y3LXnkaDOY8PPlvTxVvZffr/C7naw4KlHcYStxmRtAECpILy2ScR4FvLlj67r460U4uwggWwfcnp8fLy3ije2lbI2rwqPTzEsMZzLJqVxyaRUhsSfWoEA4MgM6pjLYPHvISr19Dao4TB88FPY/74ukPH1F0650qPH52ft/mpe31rCx3urcPv8jEqJ5Mqp6VwxJb3z03/0ALfPzXuF7/F87vMcsB0gOSyZG8bcwPLs5URY+7bAVG+obqnmg6IPeL/ofXJrczEwmJ4ynTlpc8iOySY7NpvksORuPZF2ep0U2grJr89nT+0eVh9aTZWjijBzGBcOuZBlWcuYMWgGZtPAa9E/HYW2Ql7c8yJvF7yNx+9hXsY8vjnum0xKnNTvAxeP38PKgyt5cc+L7KndQ1hQFE3V02mpnUFQ6GGsCZ8QFFJOrDWFe6fewaXDLh3wPRq6W4OzgZf3vcxLe1+i0d2Iv2U4zup5KG8U1vi1mKO3YTaZuGTYxdwy7pYB2dDRk3x+H6sPr+aZXc+wr24fyhODq/Y8vI0TMEdvIzj+MwyznQkJE7h1/K3MzZjb78fU97b8+nye2fUMHx78EIMgPA1Taak+D5O1BmvCJ5jDDhEeFM2tE2/m6pFXE2k9i6YrOgUtnhZezXuV53Ofp9pRjXIMwVEzH78rCWvc51hiNmOYPMxNn8s3x32TKUlT+v13uxADmQSy/URDi5v3d+npZzYV6SlyJw+O4bJJaSybMOjUpp/xumD9n+GzP4DJAhf8EqbfpjO1J+Lz6O7Jax7S1+f9DGbdCSc5CW3tMv361hLe3l5GbbOb+HArl0xKZfmU9M5ll3uQzWXj1bxXeXnvy1Q7qhkRO4Kbx97M4qzFA7r7alcU2Yr4oOgDPij6gIONB9tuj7RGtgW12THZDI8dTnZsNlHWqBM+n1/5KbGXkF+fT15DHvn1+eTX53PYfhi/0t2xgoOCOSf1HJYNXca89HmEmE9zKpgzQI2jhn/v+zf/3v9vbC4bExIncPPYm7kg44J+Nxa49XOzYt8KqlqqyIrO4vrR13PxsItZuauOR1fup6zBwaCYEC45p4Ed9tfYXbublPAUbhl3C1dkX0FwUN9Nn9UfVLVU8ULuC7yS9woOr4P5GfO5dfytFJbEt+2/1JhQbpsfQxkf8t/8/+L2uVmUuYhbx9/KqLhRff0W+pTH5+Hdwnd5dvezHGw8SGZUJreMuwVf42QeX1XYtv++tzALFbGZZ3c/S2lTKcNjhnPr+Fv5WubXzvhGspPZWb2Tv+/6O2uK1xBqDuXqkVdzw5gb+GK/54hj8KpzPexzvsH60vVEWCL4xqhvcN3o6/r1cIje0OhuZMXeFfxr779ocDUwM2Umt0+4ndLyVP7wUV7b/rtzQQqN1rWs2LuCelc9ExIncMvYW5iXMa/ffbcLcSaQQLYfKm1w8Pb2Mt7aXsq+CjtBJoNzhyewbMIgzstOOPnY0roiXWW4YBWkjIeL/g/Sj/M3PvwlvPt9qNoDI5fCkodPmIX1+xV7yhtZl1/N29vL2FdhxxpkYsHoJJZPSWfuyMR+Nc7Ir/zkVOTw9oG3+ejQRzi8Ds4ZdA43j72Zc1LP6VeBdl+zuWw68GzIp6C+oG1t99jb7mMcVYtTcfzPvoFBemQ6I2JHtAXE2bHZDI4cLD/oAS2eFt468BYv5L5ASVMJGZEZXDPyGuYPnt+nhaG8fi+bKzaz8uBK3i96H4fXwaxBs7hxzI3MSZtzwiyXUoovyr7gqZ1Psa1qG4mhidw89maWj1hOuOU0epmcAYpsRby450XeLHgTn/KxJGsJ3xr3LbJjs0/4uBpHDS/tfYl/7/s3TZ4mzk07l1vG3cK05Gln1XeW3W3nrYK3eC73OSpbKhkdN5pbx9/KgsELTvgd4vV7+fDgh/xj1z8oaCggPSKdb477JhcNvYgwSy/Pv96H/MrPhrIN/HP3P9lYsZHo4GiuG3Ud146+lujg6BM+dk/tHp7Z9QyrD60mOCiYK7Kv4Pox1/ergnW9ocZRw8t7X2bFvhU0eZqYmz6X2ybcxsTE48wDHeDwOnir4C2ez32ekqYSMqMyuW70dSzJWnLSfS+EOHUSyPZz+yvsvLm9lLe3l1Ha4ABgWGI4c4YnMGd4ArOGxhMdeoxsYusE3h/8VE+XMfVmuPA+CNWFCGip04Witr4AUemw9BEYteyY21Bhc/JZfjWfF9TweX4NtYEqkxMzYrhyShoXT0ztvrlyu0lhQyHvFL7Du4XvUtFcQbglnEVDFnHt6GvP+uzG6VBKUdlSSV59HgUNBbR4Wr5yn6NPrAeFDyI7JpthMcPOqpPGrvD5fXx8+GOez32enTU7ARgaPZS56eD655sAABr+SURBVHM5P/18JiVN6vGMksfnYUP5BlYfWs0nxZ9gc9kINYeyaMgibhhzw2lPP6WUYnPFZp7a+RSbKjYRZg5j6dClXDniSsbGj+2hd9H3XD4Xqw6t4vW818mpzNFVw4dfzs3jbj7tIKDR3ch/9v2HF/e8SL2rnsyoTJZnL+eS4ZcQF9K5+WX7O6UUO6p38Hr+66w8uBKH18HU5KncNv42ZqfOPq1A3q/8rClewzO7nmFXzS7CLeEszdLH4Jj4MT34LvpWdUs1bxa8yev5r1PaVEpiaCI3jb2Jq0ZcddrfyYW2Qp7d9SzvFb6HV3mZNWgWy0csZ0HGgjN26IDP72ND+QZey3uNNcVr8Cs/C4cs5LYJt532+YPX72X14dU8t/s5cmtzsZgszE2fy0XDLuL8tPPP2H0oRG+RQHaA8PsV+yrsrC+oYf2BGjYW1uHw+DAZMCE9hnMDge2kjBhCrR1aql12+PT3sPFJXJYYHlHX09Ds5n+sLxNtNGM6506Y+7O2+eKUUjQ6vGwtruezvBo+y68mv6oJgISIYM7LTuC87ATOHZ5AUlT/6hJa56zjg6IPeOfAO+TW5hJkBDE7dTYXD7uYeRnzCDX3nyrJQhxPcWMx60rXsbZ4LZsrN+P1e4m0RnJu6rmcn3E+56ae223FyJxeJ+vL1rP60GrWFq/F7rETbglnXsY8Fg5eyOy02d3yudlZvZNX817lw6IPcfqcjI4bzZUjrmRp1tIzZlz6gYYDvJb3Gu8UvoPNZSMjMoMrsq/gsuGXdXkeZ4fXwapDq3gt7zW2VW3DbDKzYPAClmcvZ+agmWfEOFCby8a7he/yWt5rFDQUEGYOY0nWEq4ccSXjEsZ16bmVUmyv3s5rea+x8uBKXD7XGXcM+vw+1pet57W811hXsg6f8jEjZQZXjriSBYMXYA3qWmNzZXMlbxS8wX/z/0t5czlxIXFcMuwSlmcvP2PGch/9HmODY7ls+GUsH7GcIVFDuvTcSin21e3jncJ3eK/wPeqcdcQEx7A4czEXD7u430+NJ0R/JYHsAOX2+tl2uJ71BTV8XlDDjhJb28TfkSHmIyZIT4oMxlK1m4VFjzDZlA/AFn82v/Z9i4RhU4kIMVPVYbqD1qkNgs0mZmTFBYLXREalRPa7L9oGZwNrS9ay6tAq1peux6u8jI4bzcXDLmZJ1pIun0AK0ZeaPc1sKNvA2pK1rCtZR52zDpNhYtj/b+9eg+Mq8zuPf/+639GtZcm2fLdkGIIsj3HsZcJ4HEKAYWICTsFUdjM1lS0qyaR2spvs1mxeTGpTm61K1W5uS3YoKjPJJJUN7BhmxgVMZoaYSkIBBiMJsMEWtvFF2FK3LrZurZa6+78vTquRZRnLbiF1y78P1dXnnD7q87g5enR+/VxO9UaayptoLGuksXzGo6yRFeUr0hetU8kpBqIDRMYjRKIRIuMRwtEwkfEIfeN9dIY7icajVBVV8YXmL3DvunvZ2bQz44veqxmZHOGFUy/wve7v0T3UTWlBaRBWNgdhJdvql2uZK2Des+YeHml5hB2NOz6VgDk7MK+uWM0jLY8sSGBebNMzYO/v3s9Pz/yUWCLG7XW380jLI9y//v5PpSv68OQwL5x6gf3d+9Pn4H3r7mNfy76cDBO9Y718/4Pv89yJ5+gd66W2pJaHNj3Ew5sfzjh8zWV2a2XCE2xfsZ19Lfu4Z+09OTcePpFM8MpHrwRfAHz0LyQ9yc6mnexr2cee5j2fSotpPBnn1fOv8vzJ5zl47iCxRIx1Vet4cMOD7Fmzh3W3rFu283a4O0lPkvAE8WSchCdIJBMkPIETvDY9n8b0suPpn4OPhzQ5zvTopvS2WXnkk4Y/LYalPv5ysKZyzSfOpbJsg2xja6P/1t/8FpVFlVQVVQXPxVVUFVYFz0VVNJQ1LJuJZkYmpnj91CDdfSOX3Y9v+t5osXgSY4rPVf0AgFeGH8IpxID19eXp+xFO31evtbGSO9fVUlKYfWMZe0Z6ePncyxw8e5COcAdJT9JY3sj96+/nSxu+dM3xZyK5KOlJjvYf5Z97/pnjg8fpHe+ld6yXi7GLV+xbW1KLYQxODF7xhzTP8qgvqae+rJ7b6m7jF9b8Anc23bmoF07uzpH+I+z/YD8/+vBHRONRWmta+dyqz9He0E5bqC0rb4EVT8Y5PnScrnAXHX0dvHbhNUYmR5aky28sEeOlMy/x7AfP8mbvm+RbPtsbt7OtYRtbG7bSFmrLujHJ7k7PaA9d4S46w5280fsGZ4bPUFFYwRc3fJF9LfsWbejH9Dn47AfPpseBb7hlA3c23sm2hm20N7TTVNG0KGW5HiOTI7wdeZvOcCcdfR10hDtwd3at3MW+ln3sXr170bqrRsYj/PDkD3m2+1l6Rnu4pfgWdjXtYmvDVtob2mmpacm6SbbcnXMj5+gMd9IZ7uSVj16hb7yPupI6fnnzL/Pwpodprlq8ccDTt8Y7cPIAh/uCxpfCvEI2VQcTLLbWtNJS20JrTSs1JTWfalncnVgixsjkSPCYGmFscoyx+BhjU8FjfGo8eI4Hz9F4lIn4BJOJSWKJGJPJyfRyLBFjKjFFLBFLB9a4xz/Vf4MsP888+MwnDgVZ0iBrZq3AMzM2bQC+6e5/NmOf3cAPgQ9Tm55z9z+81ntXb6r2tj9qY3hymHjy6r84NcU1NJYHrRhN5U3pVo2miiaaypsW/DYkSyEyHuGu//M/KKx+g7zCYQCSU1VMXdxB/OIOTv33Ly9xCT+Zu3N86DgHzx7k4NmDHB86DsDmms3sad7DnjV7uLX21pz//yRyI6LxKH1jfelgO/0AaChrIFQWoqG0gfqyehpKG6gtqc2qybZGJ0d58cMXOXDyAEf7j6YvdNbfsp6toeCCuK2hjfVV6xf9d3xsaiwdGjrDnbwTeYdoPJirYGX5SrY3buehTQ8t+SRMpy+d5rkTz/Ha+dfoHuom6UnyLI/WmtZ0qGhvaKexvHFRyzWVnOL44HE6+jroigThtT/aD0BFYQVtDW3ct+4+7l1775KOpx+bGuPFD1/kJ6d/wjuRdxiPB3MBrChbQXtD+5IFM3fno9GP6Ax3BuE/0smJoRM4Tr7l01LTws+t/jke3vwwqypWLVq5Zkt6kkMXDnHg5AHe7H2TvvE+AMoKyrgjdEf6y5U7Qncs+pcrk4lJ3ht4L/3lSVeki8GJ4M4QlYWVfHbFZ9m7aS+fb/78kreCnh89z1t9bwWz/g91c3zoePr3BSBUGqKlpoXG8kaK84spKSihJL+E4oLiYD21XJJfQsIT6ZAZjUfnfIxOjqYD63R4nUpOXbOcBVZAWWEZ5YXllBWUUVJQQlF+EUX5RRTnB2UpzCukOL84vb0gr4ACKyA/L598y6cgr4B8yyfP8ijIKyDP8si3fMwMw8izPAzD7PLl6Qkn088ztn38dHldPHt9senaNDM7Gnd84gRpWdMia2b5wEfAz7r7mRnbdwO/5+4PXs/7TXctdncmEhOMTI4wHBtmeHKYkckRLsYu0jfel77ouzB2gb6xvstmaIWgIt5YvZGN1RvZVL2JDbdsYFP1JhrLG7P65JzusvX0sad56cxLxD1OfHQzk0O7ACiqeZ2Cim7wPO5ddw+PbXlsyS/GZhqaGOJQ7yFeP/86r55/lQtjFzCM9oZ29qzZw57mPYv6ramIfPqi8ShH+4/SFemiK9xFV6SLS7FLAFQXV/Mz9T/DyoqVhEpD1JfWpx91pXXUldRdV0uUuzM8Ofxxt+tU1+tINEJ4PMy5kXOXhcKWmhbaG9rTF+WLHQrna3RylHci79AZuTJ8N5Y3srF6I6HSUPozDJVdvjzfrqGzP7/+aH/6M+yP9tM33sexwWPpY6+qWBUEwlAQDDdVb8qqL1SmxZNxuoe6Pw6P4c7LgtmW2i2sKFtBfVn9nJ9hVdH8bzsXjUfpH+9Pn3/90f7LPs8TQycIR8MAlBeW0xZqS4fqO+rvyNrJ9C6MXqAj3JH+DLuHunE8/XvUXNkcfG6zPr9QWYjq4up5d8l3d0amRuiP9tM/3p8+BweiA0SiEXpGenhv4D0mk8EElc2VzekeH+0N7Wys3pj148sHogN0D3Vf9hiIDjCRmCAWj6X/bfNRkl9CaUEpJQXBc0VRRdBbsbAqvTyzF2NFYQUVRRWUFQShtbywnLLCMoryirLmWlEkm4LsvcAfuPtds7bvJoMge73GpsbSwfb86HlOXjzJyYsnOXHxBAMTA+n9ygvL2XhLEHC31G7h1rpbaa1pXfI/LONT4zx/6nmeOf4M3UPdVBZWsnfTXuqSu/lfLwwRnUqk9y0tG+Lznz3BkZGXuBS7xIZbNvBo66N8aeOXFv0m6LFEjI6+Dl6/8DqvnX8tuNk9TmVhJXc23snu5t3cvfrum/5ediI3k6QnOT18Ogi14S6ODBwhPB5Oh9vZqourqSmpIY/g4tSn/5vxN8pxphJT9Ef757wIrCisoL60nqbyJtoa2mgPtXNH6I6cnRBoZnfoznAn50bO0T/ez8DEAAlPXLF/ZVElpfmfPMGX41yKXZrz8ysrKCNUFgSULbVb0uF1RfmKBfs3LbYLoxfSrfLdQ93poDkd0mcqzi8Owuw1WoHG4+OMTo1esb3ACqgrrSNUGqK5qjndzTlbg/98jEyOBF+uhDt5t/9desd6iYxHrmg4gODfX1tSm275niswGUbCEwxODBJLxK54vSiviFBZiBVlK7i9/vZ0q3qujSGfj0Qyke7GG0vEmIhPEEvEyLO8dGidbjXN9tAuciOyKch+B+hw9ydmbd8NPAv0AOcJQu3Ra73fpzHZ08WJi5y89HGwPXnxJB8MfcBQbCgoK8baqrVsqd0ShNvaW9lSt+VTHTPl7nw4/CGdfZ10hDs4ePYgo1OjtNa08tiWx3hg/QPpcP2Dzo8uu/H5f/7FVh5qX8VEfIIfn/4xTx97miMDRygtKOXu1XcHlX9oKy21LQva5cbdiUQjdA91c2zwGG9ceIOOcAexRIyCvALaQm3satrFzpU7+UzdZ7JujI2ILK3JxCQD0YGgJSbaT/9Ef7pVZro+njazS9p0uCjI+zgszGwNCpWGlvzLyMWSSCYYig2lWwFntqjOp3thVXHVnC27N8vn5+6MTY3N2ZI6MnllQJutOL84HfhnfobX0yKZ66Lx6BWfXXg8zODEYHpiH+CKL6IgGOtfW1Kb7pWR7qVRVk9lYfZNTCkin46sCLJmVkQQUj/j7n2zXqsCku4+amYPAH/u7nPO5mNmjwOPA6xZs+azZ86cmWu3BTV9n81jg8d4f/B9jg0c49jgMc6PnU/vEyoNsbZqLWur1rKmag1rK4Pn5srm655saioxxXuD76WDa2e4Mz3ZS21JLbtW7uLR1kfZGtp6QxX5kf4jfK/7e7x6/tX0OLuS/JLLvtVsC7XN+4besUSMUxdPcXzoeNAtZjDoGjPzYnNT9SZ2Nu1k18pdbF+x/aa5EBIRERERkRuTLUF2L/A1d793HvueBra7e/8n7bfUt9+5FLvEscEg1HYPdXN2+CxnR86mJxmY1ljeyNrKtTSUNWBmuPvlU4vz8beR4fEw7/a/m+5Ks7ZqbXrMVntDO2ur1i7ot5C9Y710Rbp4O/w2XeEujg0eu2wilrqSOpKeJO5xEsnEFctTySl6x3rT3ddK8kvYVL2JltoWWmo+fsw3FIuIiIiIiED2BNmngR+7+1/P8Voj0OfubmY7gP3AWr/GQZc6yF7N8OQw54bPcWb4DGdGzgQBd/gskWgk3QVu2uxZ2aqLq9nasJVtDdtoa2hb9PEeMydieTv8NqNTo+RbPvl5qZnnLDXzXF5+sJyXx8rylemp49dUrsnZMT4iIiIiIpI9ljzImlkZcA7Y4O6XUtt+A8DdnzSz3wZ+E4gDUeA/ufur13rfbA2yIiIiIiIikpmFCLIZzcTj7uNA3axtT85YfgJ4YvbPiYiIiIiIiNyom2N6PREREREREVk2FGRFREREREQkpyjIioiIiIiISE5RkBUREREREZGcoiArIiIiIiIiOUVBVkRERERERHKKgqyIiIiIiIjkFAVZERERERERySkKsiIiIiIiIpJTFGRFREREREQkpyjIioiIiIiISE5RkBUREREREZGcknGQNbPTZvaumXWZ2eE5Xjcz+wszO2Fm75jZtkyPKSIiIiIiIjevggV6ny+4e/9VXrsf2Jx6/CzwrdSziIiIiIiIyHVbjK7Fe4G/9cDrQLWZNS3CcUVERERERGQZWogg68BPzOwtM3t8jtdXAedmrPektl3GzB43s8NmdjgSiSxAsURERERERGQ5Wogge5e7byPoQvw1M7t71us2x8/4FRvcn3L37e6+PRQKLUCxREREREREZDnKOMi6+/nUcxj4PrBj1i49QPOM9dXA+UyPKyIiIiIiIjenjIKsmZWbWeX0MnAvcGTWbgeAX0vNXrwTuOTuFzI5roiIiIiIiNy8Mp21eAXwfTObfq//6+7/aGa/AeDuTwIvAg8AJ4Bx4KsZHlNERERERERuYhkFWXc/BbTNsf3JGcsOfC2T44iIiIiIiIhMW4zb74iIiIiIiIgsGAVZERERERERySkKsiIiIiIiIpJTFGRFREREREQkpyjIioiIiIiISE5RkBUREREREZGcoiArIiIiIiIiOUVBVkRERERERHKKgqyIiIiIiIjkFAVZERERERERySkKsiIiIiIiIpJTbjjImlmzmb1sZu+b2VEz+/oc++w2s0tm1pV6fDOz4oqIiIiIiMjNriCDn40Dv+vuHWZWCbxlZj919/dm7fev7v5gBscRERERERERSbvhFll3v+DuHanlEeB9YNVCFUxERERERERkLgsyRtbM1gHtwKE5Xt5lZm+b2Y/M7DOf8B6Pm9lhMzsciUQWolgiIiIiIiKyDGUcZM2sAngW+B13H571cgew1t3bgP8N/OBq7+PuT7n7dnffHgqFMi2WiIiIiIiILFMZBVkzKyQIsX/v7s/Nft3dh919NLX8IlBoZvWZHFNERERERERubpnMWmzAt4H33f1PrrJPY2o/zGxH6ngDN3pMERERERERkUxmLb4L+HfAu2bWldr2+8AaAHd/EtgH/KaZxYEo8Ji7ewbHFBERERERkZvcDQdZd38FsGvs8wTwxI0eQ0RERERERGS2BZm1WERERERERGSxKMiKiIiIiIhITlGQFRERERERkZyiICsiIiIiIiI5RUFWREREREREcoqCrIiIiIiIiOQUBVkRERERERHJKQqyIiIiIiIiklMUZEVERERERCSnKMiKiIiIiIhITlGQFRERERERkZySUZA1s/vM7LiZnTCzb8zxerGZPZN6/ZCZrcvkeCIiIiIiIiI3HGTNLB/4S+B+4Dbgy2Z226zdfh0YcvdNwJ8Cf3yjxxMRERERERGBzFpkdwAn3P2Uu08CTwN7Z+2zF/huank/8PNmZhkcU0RERERERG5ymQTZVcC5Ges9qW1z7uPuceASUJfBMUVEREREROQmV5DBz87Vsuo3sE+wo9njwOOp1ZiZHcmgbCLZoB7oX+pCiCwAncuyHOg8luVC57IsB62ZvkEmQbYHaJ6xvho4f5V9esysALgFGJzrzdz9KeApADM77O7bMyibyJLTeSzLhc5lWQ50HstyoXNZlgMzO5zpe2TStfhNYLOZrTezIuAx4MCsfQ4AX0kt7wMOuvucLbIiIiIiIiIi83HDLbLuHjez3wZ+DOQD33H3o2b2h8Bhdz8AfBv4OzM7QdAS+9hCFFpERERERERuXpl0LcbdXwRenLXtmzOWJ4BfuYG3fiqTcolkCZ3HslzoXJblQOexLBc6l2U5yPg8NvX0FRERERERkVySyRhZERERERERkUWXVUHWzO4zs+NmdsLMvrHU5RGZLzNrNrOXzex9MztqZl9Pba81s5+a2Qep55qlLqvItZhZvpl1mtnzqfX1ZnYodR4/k5rgTySrmVm1me03s2OpunmX6mTJNWb2H1PXFUfM7B/MrER1suQCM/uOmYVn3lL1anWwBf4ilQHfMbNt8zlG1gRZM8sH/hK4H7gN+LKZ3ba0pRKZtzjwu+5+K7AT+Frq/P0G8E/uvhn4p9S6SLb7OvD+jPU/Bv40dR4PAb++JKUSuT5/Dvyju28B2gjOadXJkjPMbBXwH4Dt7n47weSqj6E6WXLD3wD3zdp2tTr4fmBz6vE48K35HCBrgiywAzjh7qfcfRJ4Gti7xGUSmRd3v+DuHanlEYILplUE5/B3U7t9F3hoaUooMj9mthr4IvBXqXUD9gD7U7voPJasZ2ZVwN0Ed0/A3Sfd/SKqkyX3FAClZlYAlAEXUJ0sOcDd/4XgrjUzXa0O3gv8rQdeB6rNrOlax8imILsKODdjvSe1TSSnmNk6oB04BKxw9wsQhF2gYelKJjIvfwb8FyCZWq8DLrp7PLWuullywQYgAvx1qpv8X5lZOaqTJYe4+0fA/wTOEgTYS8BbqE6W3HW1OviGcmA2BVmbY5umVJacYmYVwLPA77j78FKXR+R6mNmDQNjd35q5eY5dVTdLtisAtgHfcvd2YAx1I5Yckxo/uBdYD6wEygm6YM6mOlly3Q1da2RTkO0BmmesrwbOL1FZRK6bmRUShNi/d/fnUpv7prtGpJ7DS1U+kXm4C/glMztNMLxjD0ELbXWqWxuobpbc0AP0uPuh1Pp+gmCrOllyyT3Ah+4ecfcp4Dng36A6WXLX1ergG8qB2RRk3wQ2p2ZiKyIYzH5gicskMi+pcYTfBt539z+Z8dIB4Cup5a8AP1zssonMl7v/V3df7e7rCOrgg+7+q8DLwL7UbjqPJeu5ey9wzsxaU5t+HngP1cmSW84CO82sLHWdMX0eq06WXHW1OvgA8Gup2Yt3ApemuyB/EnPPnt4IZvYAwbf/+cB33P2PlrhIIvNiZp8D/hV4l4/HFv4+wTjZ/wesIfiD9CvuPnvgu0jWMbPdwO+5+4NmtoGghbYW6AT+rbvHlrJ8ItdiZlsJJi0rAk4BXyX4Al91suQMM/tvwKMEd0foBP49wdhB1cmS1czsH4DdQD3QB/wB8APmqINTX9Q8QTDL8TjwVXc/fM1jZFOQFREREREREbmWbOpaLCIiIiIiInJNCrIiIiIiIiKSUxRkRUREREREJKcoyIqIiIiIiEhOUZAVERERERGRnKIgKyIiIiIiIjlFQVZERERERERyioKsiIiIiIiI5JT/D6BHPBNqso91AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c143080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class GaussianProcess(object):\n",
    "    def compute_kernel(X, Y, kernel_fn):\n",
    "        x_indices, y_indices = np.meshgrid(np.arange(X.shape[0]), np.arange(Y.shape[0]))\n",
    "        return kernel_fn(X[x_indices], Y[y_indices])\n",
    "\n",
    "    def __init__(self, x, y, kernel_fn, mean_fn, sigma=0):\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.kernel_fn = np.vectorize(kernel_fn, signature='(),()->()')\n",
    "        self.mean_fn = mean_fn\n",
    "        self.sigma = sigma\n",
    "        self.K = compute_kernel(self.x, self.x, self.kernel_fn)\n",
    "        self.inv_K = inv(self.sigma*np.identity(self.x.shape[0]) + self.K)\n",
    "        \n",
    "    def mean(self, x):\n",
    "        k_x_D = compute_kernel(self.x, x, self.kernel_fn)\n",
    "        print(k_x_D.shape, self.inv_K.shape)\n",
    "        return self.mean_fn(x) + np.dot(np.dot(k_x_D, self.inv_K), (self.y - self.mean_fn(self.x)))\n",
    "\n",
    "    def covariance(self, x):\n",
    "        k_x_x = compute_kernel(x, x, self.kernel_fn)\n",
    "        k_x_D = compute_kernel(self.x, x, self.kernel_fn)\n",
    "        print('k_x_x.shape',k_x_x.shape)\n",
    "        print('k_x_D.shape',k_x_D.shape)\n",
    "        kk = np.dot(np.dot(k_x_D, self.inv_K), k_x_D.T)\n",
    "        print('kk.shape', kk.shape)\n",
    "        cov = k_x_x - kk + self.sigma*np.identity(x.shape[0])\n",
    "        print('cov.shape', cov.shape)\n",
    "        return cov\n",
    "        \n",
    "        \n",
    "inst_p = np.round(10+1*np.random.normal(size=8)).astype(np.int)\n",
    "t = np.cumsum(inst_p)\n",
    "gp = GaussianProcess(t, inst_p, sigma=0.001, kernel_fn=gaussian(3), mean_fn=lambda x: 10+np.zeros_like(x))\n",
    "mm = gp.mean(x)\n",
    "cov = gp.covariance(x)\n",
    "\n",
    "plt.figure(figsize=(16,3))\n",
    "plt.scatter(t, inst_p)\n",
    "plt.plot(x, mm)\n",
    "# plt.plot(x, cov*10)\n",
    "stds = np.sqrt(cov.diagonal())\n",
    "print(cov.diagonal())\n",
    "plt.plot(x, mm + 2*stds)\n",
    "plt.plot(x, mm - 2*stds)\n",
    "# plt.plot(x, np.random.multivariate_normal(mm, cov))\n",
    "# plt.patch(x, mm + 2*stds, alpha=0.2, line_width=0, color='grey', legend='Uncertainty')\n",
    "\n",
    "plt.ylim(0,20)\n",
    "plt.xlim(0,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Attachments",
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}