{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook, I present a practical example of what is actually happening when you train a Gaussian process model. I don't suggest you ever use this implementation to do actual machine learning. I have coded everything with readability in mind, at the cost of efficiency. There are also a bunch of cool linear algebra and numerical stability tricks that I have also skipped in the interest of simplicity.\n", "\n", "The leading libraries for working with GPs (in Python) are;\n", "* [GPy](https://sheffieldml.github.io/GPy/) - Probably the most popular\n", "* [George](https://george.readthedocs.io/en/latest/) - Written by an astrophysicist I had beers with once, so maybe I'm biased, but this is the one I used to use as it gives the user a lot of flexibility and plays well with stochastic sampling methods\n", "* [Scikit-Learn](https://scikit-learn.org/stable/modules/gaussian_process.html) - Used to be considered a poor implementation, not sure if that is still true" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, some notation. We have;\n", "\n", "* $x_t$ - the inputs for our training data\n", "* $y_t$ - the outputs for our training data\n", "* $x_p$ - the inputs for which we want to predict the output\n", "* $y_p$ - the (unknown) outputs that we want to predict\n", "\n", "We generate some training data;" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAGuCAYAAACnXgEDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdqklEQVR4nO3df7Ddd13n8dfbNOhFVoM2KkmrrTs1WBU2mmVRdOWXpsUfjY7rwrqILNJhFxQdJ0uzOyPj6oyy2XXUFS0d7KLjj45CJnSdahZQl1kVJRAklBrMFKFJqg1i/HkX0vDeP+5NuQk3P25y7+ece+/jMZPJOZ/vtyfvzpxpn/fk8/2e6u4AAAAr69MmPQAAAKwHwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMcM2kBxjl2muv7RtuuGHSYwAAsIa9613v+kh3b17s2LoJ7xtuuCEHDx6c9BgAAKxhVfWhCx2z1QQAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGmLrwrqq7q+qRqnrfBY5/V1W9d/7XH1TVU0fPCAAASzV14Z3kDUluucjxDyb5+u5+SpIfTXLXiKEAAOBqXDPpAc7X3W+vqhsucvwPFjx9R5LrVnwoAAC4StP4ifdSvCTJb13oYFXdXlUHq+rgyZMnB44FAADnmrpPvC9XVT0rc+H9tRc6p7vvyvxWlB07dvSg0QAAsv/Q8ew9cCQnTs1my6aZ7N65Lbu2b530WEzQqgzvqnpKktcnubW7/2rS8wAALLT/0PHs2Xc4s6fPJEmOn5rNnn2Hk0R8r2OrbqtJVX1hkn1JXtjdH5j0PAAA59t74Mhj0X3W7Okz2XvgyIQmYhpM3SfeVfVrSZ6Z5NqqOpbk1Uk2Jkl335nkh5N8bpKfq6okebS7d0xmWgCAT3Xi1OyS1lkfpi68u/sFlzj+vUm+d9A4AABLtmXTTI4vEtlbNs1MYBqmxarbagIAMO1279yWmY0bzlmb2bghu3dum9BETIOp+8QbAGC1O3sBpbuasJDwBgBYAbu2bxXanMNWEwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABhAeAMAwADCGwAABhDeAAAwgPAGAIABhDcAAAwgvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABhAeAMAwADCGwAABhDeAAAwwNSFd1XdXVWPVNX7LnC8qupnqupoVb23qr5y9IwAALBUUxfeSd6Q5JaLHL81yU3zv25P8vMDZgIAgKsydeHd3W9P8tGLnHJbkl/qOe9IsqmqnjRmOgAAuDJTF96XYWuShxY8Pza/9imq6vaqOlhVB0+ePDlkOAAAWMxqDO9aZK0XO7G77+ruHd29Y/PmzSs8FgAAXNhqDO9jSa5f8Py6JCcmNAsAAFyW1Rje9yb57vm7mzw9yd9098OTHgoAAC7mmkkPcL6q+rUkz0xybVUdS/LqJBuTpLvvTHJfkuclOZrkH5O8eDKTAgDA5Zu68O7uF1zieCd5+aBxAABgWazGrSYAALDqCG8AABhAeAMAwADCGwAABhDeAAAwgPAGAIABhDcAAAwgvAEAYICp+wIdAGD12H/oePYeOJITp2azZdNMdu/cll3bt056LJhKwhsAuCL7Dx3Pnn2HM3v6TJLk+KnZ7Nl3OEnENyzCVhMA4IrsPXDkseg+a/b0mew9cGRCE8F0E94AwBU5cWp2Seuw3glvAOCKbNk0s6R1WO+ENwBwRXbv3JaZjRvOWZvZuCG7d26b0EQw3VxcCQBckbMXULqrCVwe4Q0AXLFd27cKbbhMtpoAAMAAwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBggGsmPcB6s//Q8ew9cCQnTs1my6aZ7N65Lbu2b530WAAArDDhPdD+Q8ezZ9/hzJ4+kyQ5fmo2e/YdThLxDQCwxtlqMtDeA0cei+6zZk+fyd4DRyY0EQAAowjvgU6cml3SOgAAa4fwHmjLppklrQMAsHYI74F279yWmY0bzlmb2bghu3dum9BEAACM4uLKgc5eQOmuJgAA64/wHmzX9q1CGwBgHbLVBAAABpi68K6qW6rqSFUdrao7Fjn+2VX1v6rqT6rq/qp68STmBACApZiq8K6qDUlem+TWJDcneUFV3XzeaS9P8v7ufmqSZyb571X1uKGDAgDAEk1VeCd5WpKj3f1gd388yT1JbjvvnE7yT6qqkjwhyUeTPDp2TAAAWJppC++tSR5a8PzY/NpCP5vkS5OcSHI4ySu7+xOLvVhV3V5VB6vq4MmTJ1diXgAAuCzTFt61yFqf93xnkvck2ZLknyX52ar6rMVerLvv6u4d3b1j8+bNyzknAAAsybSF97Ek1y94fl3mPtle6MVJ9vWco0k+mOTJg+YDAIArMm3h/c4kN1XVjfMXTD4/yb3nnfPhJM9Jkqr6/CTbkjw4dEoAAFiiqfoCne5+tKpekeRAkg1J7u7u+6vqZfPH70zyo0neUFWHM7c15VXd/ZGJDQ0AAJdhqsI7Sbr7viT3nbd254LHJ5J84+i5AADgakzbVhMAAFiThDcAAAwgvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABjgmkkPwPTaf+h49h44khOnZrNl00x279yWXdu3TnosAIBVSXizqP2HjmfPvsOZPX0mSXL81Gz27DucJOIbAOAK2GrCovYeOPJYdJ81e/pM9h44MqGJAABWN+HNok6cml3SOgAAFye8WdSWTTNLWgcA4OKEN4vavXNbZjZuOGdtZuOG7N65bUITAQCsbi6uXCGr/Y4gZ2ddzf8OAADTRHivgLVyR5Bd27euqnlX2mr/YQoAmCxbTVaAO4KsPWd/mDp+ajadT/4wtf/Q8UmPBgCsEsJ7BbgjyNrjhykA4GoJ7xXgjiBrjx+mAICrJbxXgDuCrD1+mAIArpbwXgG7tm/Nj3/7V2TrpplUkq2bZvLj3/4VLsRbxfwwBQBcLXc1WSHuCLK2uL0iAHC1hDdcJj9MAQBXw1YTAAAYQHgDAMAAUxfeVXVLVR2pqqNVdccFznlmVb2nqu6vqv8zekYAAFiqqdrjXVUbkrw2yTckOZbknVV1b3e/f8E5m5L8XJJbuvvDVfV5ExkWAACWYNo+8X5akqPd/WB3fzzJPUluO++cf5NkX3d/OEm6+5HBMwIAwJJNW3hvTfLQgufH5tcW+pIkT6yq36uqd1XVd1/oxarq9qo6WFUHT548uQLjAgDA5ZmqrSZJapG1Pu/5NUm+Kslzkswk+cOqekd3f+BT/sHuu5LclSQ7duw4/3VWvf2HjruvNADAKjFt4X0syfULnl+X5MQi53yku/8hyT9U1duTPDXJp4T3Wrb/0PHs2Xc4s6fPJEmOn5rNnn2Hk0R8AwBMoWnbavLOJDdV1Y1V9bgkz09y73nnvDnJ11XVNVX1+CT/IskDg+ecuL0HjjwW3WfNnj6TvQeOTGgiAAAuZqo+8e7uR6vqFUkOJNmQ5O7uvr+qXjZ//M7ufqCqfjvJe5N8Isnru/t9k5t6Mk6cml3SOgAAkzVV4Z0k3X1fkvvOW7vzvOd7k+wdOde02bJpJscXiewtm2YmMA0AAJcybVtNuEy7d27LzMYN56zNbNyQ3Tu3TWgiAAAuZuo+8ebynL2A0l1NAABWB+G9iu3avlVoAwCsEraaAADAAMIbAAAGEN4AADCA8AYAgAGuOLyrasOlzwIAAJLLDO+qemJV/fuqelNVPVRVH0vy8ar6m6p6Z1X9VFV97QrPCgAAq9ZFbydYVTckeXWS5yf56yTvSPL6JB9J8rEkm5LckOTpSV5eVQ8m+bEkv9zdvVJDAwDAanOp+3gfTnJPkud29+9f7MSq+twk35HkjiTXJfnxZZkQAADWgEuF97buPnE5L9Tdf5XkdUleV1VfcNWTAQDAGnLRPd6XG92L/HN/cWXjAADA2nTZdzWpqv9SVYt+Ql5Vn1tVv7F8YwEAwNqylNsJvjLJH1fVly9crKpvS/L+JNuXczAAAFhLlhLeT01yKsnBqrqjqq6tql9J8qYk++aPAwAAi7jUxZWP6e4/T/Lsqvr+JK9J8iNJ/iLJzu5+y8qMBwAAa8OSvrmyqp6Q5ClJPj3JR5N8RpInrMBcAACwpizl4spnJXlfkluSfHOSGzN3j+83VtUvV9UTV2ZEAABY/Zbyifdbk/x+kq/o7vu6+/919yuTPCfJ12Tuy3YAAIBFLCW8/3V3f1d3//XCxe7+vcxtP/nN5RwMAADWkqVcXPnGixz7+yQvW5aJAABgDbroJ95V9XVLfcGq+uyq+oorHwkAANaeS201+fWq+v2q+neXuniyqp5RVf8jyYeSfPWyTQgAAGvApbaafHGS70/y6iSvq6oPZO7OJh9J8rEkmzJ3d5PtSWaS3Jfkud19cKUGBgCA1eii4d3ds0leU1X/NXN3L3l2kq9K8uTM3cP7o0mOJPnVJG/u7kdWdlwAAFidLhreVfWFSR7u7tOZu53gW4dMBQAAa8yl9nh/MHPbSFJVv1NVT175kQAAYO25VHjPJnn8/ONnJvmsFZ0GAADWqEtdXHkoyU9X1Vvmn39fVT18gXO7u1+1fKMBAMDacanwfmmSvUluS9KZu8DyYxc4t5MIbwAAWMSl7mryp0m+JUmq6hNJdnX3H48YDAAA1pLL/sr4zN2v+8RKDQIAAGvZZYd3d39oJQcBAIC17FJ3NTlHVX1aVT1YVV+28PFKDQcAAGvFksI7SSW5Icmnn/cYAAC4iKWGNwAAcAWENwAADCC8AQBggKkL76q6paqOVNXRqrrjIuf986o6U1XfMXI+AAC4ElMV3lW1Iclrk9ya5OYkL6iqmy9w3muSHBg7IQAAXJmpCu8kT0tytLsf7O6PJ7knc19Xf77vS/KmJI+MHA4AAK7UtIX31iQPLXh+bH7tMVW1Ncm3JbnzUi9WVbdX1cGqOnjy5MllHRQAAJZiqeH9iSQ/krmvjl/4eLnUImt93vOfSvKq7j5zqRfr7ru6e0d379i8efNyzAcAAFfksr8yvqq+Ocl93f0jC5Z/5ELnX6FjSa5f8Py6fGrY70hyT1UlybVJnldVj3b3/mWeBQAAls1SPvF+c5LjVfWaqnryCs3zziQ3VdWNVfW4JM9Pcu/CE7r7xu6+obtvSPLGJP9BdAMAMO2WEt7/NMldSb4zyf1V9YdV9dKq+qzlGqa7H03yiszdreSBJL/e3fdX1cuq6mXL9ecAAMBo1X3+FurL+Ieqnp3kxZm7yLGS7Etyd3f/7vKOt3x27NjRBw8enPQYAACsYVX1ru7esdixK7qrSXf/Tne/MMmXJHlXku9K8taq+mBV/WBVXfbecQAAWA+uKLyr6uur6g1JjiT58sx96c03JvmNzF1w+UvLNSAAAKwFS7mryRcledH8rxuS/F6S25Ps6+6PzZ/2tqr6wyS/vLxjAgDA6raULSEPZu7Wfm/I3H7uD17gvPuT/PFVzgUAAGvKUsL7W5L8dnd/4mIndfcHkjzrqqYCAIA15rLDu7vvW8lBAABgLbuiiysBAIClEd4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABhAeAMAwADCGwAABhDeAAAwgPAGAIABhDcAAAwgvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABhAeAMAwABTF95VdUtVHamqo1V1xyLHv6uq3jv/6w+q6qmTmBMAAJZiqsK7qjYkeW2SW5PcnOQFVXXzead9MMnXd/dTkvxokrvGTgkAAEs3VeGd5GlJjnb3g9398ST3JLlt4Qnd/Qfd/dfzT9+R5LrBMwIAwJJNW3hvTfLQgufH5tcu5CVJfutCB6vq9qo6WFUHT548uUwjAgDA0k1beNcia73oiVXPylx4v+pCL9bdd3X3ju7esXnz5mUaEQAAlu6aSQ9wnmNJrl/w/LokJ84/qaqekuT1SW7t7r8aNBsAAFyxafvE+51JbqqqG6vqcUmen+TehSdU1Rcm2Zfkhd39gQnMCAAASzZVn3h396NV9YokB5JsSHJ3d99fVS+bP35nkh9O8rlJfq6qkuTR7t4xqZkBAOByVPeiW6jXnB07dvTBgwcnPQYAAGtYVb3rQh8KT9tWEwAAWJOENwAADCC8AQBgAOENAAADCG8AABhgqm4nCEyf/YeOZ++BIzlxajZbNs1k985t2bV966THAoBVR3gDF7T/0PHs2Xc4s6fPJEmOn5rNnn2Hk0R8A8AS2WoCXNDeA0cei+6zZk+fyd4DRyY0EQCsXsIbuKATp2aXtA4AXJjwBi5oy6aZJa0DABcmvIEL2r1zW2Y2bjhnbWbjhuzeuW1CEwHA6uXiSuCCzl5A6a4mAHD1hDdwUbu2bxXaALAMbDUBAIABhDcAAAwgvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABhAeAMAwADCGwAABhDeAAAwgPAGAIABhDcAAAwwdeFdVbdU1ZGqOlpVdyxyvKrqZ+aPv7eqvnIScwIAwFJMVXhX1YYkr01ya5Kbk7ygqm4+77Rbk9w0/+v2JD8/dEgAALgCUxXeSZ6W5Gh3P9jdH09yT5LbzjvntiS/1HPekWRTVT1p9KAAALAU0xbeW5M8tOD5sfm1pZ6TJKmq26vqYFUdPHny5LIOCgAASzFt4V2LrPUVnDO32H1Xd+/o7h2bN2++6uEAAOBKTVt4H0ty/YLn1yU5cQXnAADAVJm28H5nkpuq6saqelyS5ye597xz7k3y3fN3N3l6kr/p7odHDwoAAEtxzaQHWKi7H62qVyQ5kGRDkru7+/6qetn88TuT3JfkeUmOJvnHJC+e1LwAAHC5piq8k6S778tcXC9cu3PB407y8tFzAQDA1Zi2rSYAALAmCW8AABhAeAMAwADCGwAABhDeAAAwgPAGAIABhDcAAAwgvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4AwDAAMIbAAAGEN4AADCA8AYAgAGENwAADCC8AQBgAOENAAADCG8AABjgmkkPAAAAy2X/oePZe+BITpyazZZNM9m9c1t2bd866bGSCG8AANaI/YeOZ8++w5k9fSZJcvzUbPbsO5wkUxHftpoAALAm7D1w5LHoPmv29JnsPXBkQhOdS3gDALAmnDg1u6T10YQ3AABrwpZNM0taH014AwCwJuzeuS0zGzecszazcUN279w2oYnO5eJKAADWhLMXULqrCQAArLBd27dOTWifz1YTAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAaYmvKvqc6rqLVX1Z/O/P3GRc66vqt+tqgeq6v6qeuUkZgUAgKWamvBOckeSt3X3TUneNv/8fI8m+aHu/tIkT0/y8qq6eeCMAABwRaYpvG9L8ovzj38xya7zT+juh7v73fOP/y7JA0mm80aNAACwwDSF9+d398PJXGAn+byLnVxVNyTZnuSPLnLO7VV1sKoOnjx5cjlnBQCAJRn6zZVV9dYkX7DIof+8xNd5QpI3JfmB7v7bC53X3XcluStJduzY0Uv5MwAAYDkNDe/ufu6FjlXVX1bVk7r74ap6UpJHLnDexsxF9690974VGhUAAJbVNG01uTfJi+YfvyjJm88/oaoqyS8keaC7f3LgbAAAcFWmKbx/Isk3VNWfJfmG+eepqi1Vdd/8Oc9I8sIkz66q98z/et5kxgUAgMs3dKvJxXT3XyV5ziLrJ5I8b/7x/01Sg0cDAICrNk2feAMAwJolvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMMDU3McbWFv2HzqevQeO5MSp2WzZNJPdO7dl1/atkx4LACZGeAPLbv+h49mz73BmT59Jkhw/NZs9+w4nifgGYN2y1QRYdnsPHHksus+aPX0mew8cmdBEADB5whtYdidOzS5pHQDWA+ENLLstm2aWtA4A64HwBpbd7p3bMrNxwzlrMxs3ZPfObROaCAAmz8WVwLI7ewGlu5oAwCcJb2BF7Nq+VWgDwAK2mgAAwADCGwAABhDeAAAwgPAGAIABhDcAAAwgvAEAYADhDQAAAwhvAAAYQHgDAMAAwhsAAAYQ3gAAMIDwBgCAAaq7Jz3DEFV1MsmHzlu+NslHJjAO08X7gMT7gE/yXiDxPmDOlbwPvqi7Ny92YN2E92Kq6mB375j0HEyW9wGJ9wGf5L1A4n3AnOV+H9hqAgAAAwhvAAAYYL2H912THoCp4H1A4n3AJ3kvkHgfMGdZ3wfreo83AACMst4/8QYAgCGENwAADLAuw7uqbqmqI1V1tKrumPQ8TEZVXV9Vv1tVD1TV/VX1yknPxORU1YaqOlRVvznpWZiMqtpUVW+sqj+d/+/CV096Jsarqh+c/3/C+6rq16rqMyY9E2NU1d1V9UhVvW/B2udU1Vuq6s/mf3/i1fwZ6y68q2pDktcmuTXJzUleUFU3T3YqJuTRJD/U3V+a5OlJXu69sK69MskDkx6CifrpJL/d3U9O8tR4P6w7VbU1yfcn2dHdX55kQ5LnT3YqBnpDklvOW7sjydu6+6Ykb5t/fsXWXXgneVqSo939YHd/PMk9SW6b8ExMQHc/3N3vnn/8d5n7n+zWyU7FJFTVdUm+KcnrJz0Lk1FVn5XkXyb5hSTp7o9396mJDsWkXJNkpqquSfL4JCcmPA+DdPfbk3z0vOXbkvzi/ONfTLLrav6M9RjeW5M8tOD5sYitda+qbkiyPckfTXgUJuOnkvzHJJ+Y8BxMzhcnOZnkf85vOXp9VX3mpIdirO4+nuS/JflwkoeT/E13/+/JTsWEfX53P5zMfWCX5POu5sXWY3jXImvuqbiOVdUTkrwpyQ90999Oeh7GqqpvTvJId79r0rMwUdck+cokP9/d25P8Q67yr5RZfeb3796W5MYkW5J8ZlX928lOxVqyHsP7WJLrFzy/Lv4aad2qqo2Zi+5f6e59k56HiXhGkm+tqj/P3NazZ1fVL092JCbgWJJj3X32b73emLkQZ315bpIPdvfJ7j6dZF+Sr5nwTEzWX1bVk5Jk/vdHrubF1mN4vzPJTVV1Y1U9LnMXTdw74ZmYgKqqzO3nfKC7f3LS8zAZ3b2nu6/r7hsy99+D3+lun3CtM939F0keqqpt80vPSfL+CY7EZHw4ydOr6vHz/494Tlxku97dm+RF849flOTNV/Ni11z1OKtMdz9aVa9IciBzVyvf3d33T3gsJuMZSV6Y5HBVvWd+7T91932TGwmYoO9L8ivzH8o8mOTFE56Hwbr7j6rqjUnenbk7Xx2Kr45fN6rq15I8M8m1VXUsyauT/ESSX6+ql2TuB7N/dVV/hq+MBwCAlbcet5oAAMBwwhsAAAYQ3gAAMIDwBgCAAYQ3AAAMILwBAGAA4Q0AAAMIbwAAGEB4A5AkqapNVXWsqn7pvPV7q+oDVfX4Sc0GsBYIbwCSJN19KslLkrywqnYlSVW9OMk3Jfme7v7HyU0HsPr5yngAzlFVr0uyK8ktSX43yeu6+1UTHQpgDRDeAJyjqp6Q5L1JtiQ5muSruvtjk50KYPWz1QSAc3T33yf5zSSfnuQXRDfA8vCJNwDnqKodSf4wyeEkX5Tky7r7LyY7FcDqJ7wBeExVfUaSdyd5MMl3JvmTJA9097dOdDCANcBWEwAW+rEkX5DkpfN3MXlRkm+qqu+Z6FQAa4BPvAFIklTVM5K8PckLu/tXF6zvTfLSJF/e3ccmNR/Aaie8AQBgAFtNAABgAOENAAADCG8AABhAeAMAwADCGwAABhDeAAAwgPAGAIABhDcAAAzw/wHs2wJBn7roHgAAAABJRU5ErkJggg==\n", "text/plain": [ "