{ "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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "np.random.seed(8235)\n", "\n", "n_data = 10\n", "measurement_uncertainty = 0.2\n", "\n", "x_t = 10 * np.sort(np.random.rand(n_data))\n", "y_t = np.sin(x_t) + measurement_uncertainty * np.random.randn(len(x_t))\n", "\n", "x_p = np.linspace(min(x_t)-1, max(x_t)+1, 200)\n", "\n", "fig, axs = plot.subplots()\n", "fig.set_size_inches(12,7)\n", "axs.set_xlabel('x', size=15)\n", "axs.set_ylabel('y=f(x)', size=15)\n", "axs.scatter(x_t, y_t);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our model is that our outputs are jointly Gaussian distributed, i.e.,\n", "\n", "$$ P(y_p|y_t) \\propto P(y_p,y_t) \\propto \\exp\\left[-\\frac{1}{2} \\left(\\begin{array}{c} y_t \\ y_p \\end{array}\\right)^T\\Sigma^{-1}\\left(\\begin{array}{c} y_t \\ y_p \\end{array}\\right)\\right] $$\n", "\n", "Where ($y_t, y_p$) is the concatenated vector of our known and unknown outputs and $\\Sigma$ is the covariance matrix between all of the outputs.\n", "\n", "The problem is, we don't know $y_p$, so we can't compute this directly. We therefore make the _assumption_ that we can approximate the covariance between two outputs ($y_t, y_p$) as a function of the inputs ($x_t, x_p$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The standard choice for this function is the \"squared exponential\" kernel, but remember that this is just a choice!\n", "\n", "$$ \\Sigma_{ij} = \\langle y_i, y_j\\rangle \\approx \\kappa(x_i, x_j) = \\exp\\left(-\\frac{(x_i - x_j)^2 }{2l^2} \\right) + \\alpha\\delta_{ij} $$\n", "\n", "Where here $l$ is our free parameter, $\\alpha$ is our (homoskedastic) measurement uncertainty and $\\delta_{ij}$ is the Kroenecker delta function. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def kernel(xi, xj, l):\n", " \"\"\"\n", " Remember, any positive definite, symmetric function of the inputs will do here!\n", " \"\"\"\n", " \n", " squared_residual = np.square(xi-xj)\n", " covariance = np.exp(-1*squared_residual/(2*l*l))\n", " return covariance\n", "\n", "def kernel_matrix(xs, l):\n", " \"\"\"\n", " Since the matrices are symmetric, we can compute them more efficiently\n", " \"\"\"\n", " \n", " n = len(xs)\n", " \n", " M = np.zeros((n, n))\n", " for i in range(n):\n", " for j in range(i):\n", " M[i,j] = kernel(xs[i], xs[j], l)\n", " M = M + M.T\n", " for i in range(n):\n", " M[i,i] = kernel(xs[i], xs[i], l) + measurement_uncertainty**2\n", " \n", " return M" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first step is to figure out what the free parameters are (train the Gaussian process model).\n", "\n", "The $\\Sigma$ in the Likelihood function above contains covariances between the training outputs ($y_t$) and the predicted outputs ($y_p$). It is convenient to consider the different blocks of this matrix;\n", "\n", "$$ \\Sigma = \\left(\\begin{array}{cc} T & C^T \\\\ C & P \\end{array}\\right) $$\n", "\n", "* $T$ - the covariance matrix for the training outputs, made by taking pairwise covariances $\\langle y_t,y_t \\rangle$\n", "* $P$ - the covariance matrix for the outputs we want to predict, made by taking pairwise covariances $\\langle y_p,y_p \\rangle$\n", "* $C$ - the covariance matrix between the known, training outputs and the unknown outputs to predict $\\langle y_p,y_t \\rangle$\n", "\n", "Since we're still assuming that our outputs are jointly Gaussian distributed, we use a Gaussian likelihood for this also, in order to find the most likely values for our free parameter $l$.\n", "\n", "$$ P(\\alpha, l | x_t, f_t) = P(T | x_t, y_t) \\propto \\frac{1}{2|T|} \\exp\\left(-\\frac{1}{2}y_t^T T^{-1} y_t \\right) $$\n", "\n", "Where the elements of $T$ can be computed using the Kernel function, which is a function of $l$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# logs are more stable, and are monotonic so work for maximisation\n", "# since l must be positive, we pass its logarithm for optimisation later\n", "def neg_log_likelihood(log_l):\n", " \n", " l = np.exp(log_l)\n", " \n", " T = kernel_matrix(x_t, l)\n", " \n", " # multivariate Gaussian log likelihood\n", " ll = -0.5*np.log(np.linalg.det(T)) - (0.5*np.dot(y_t,np.linalg.solve(T, y_t))) - np.log(2.*np.pi)\n", " \n", " print(f\"Log Likelihood: {ll}\")\n", " \n", " return -ll" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For simplicity, I will optimise the likelihood using a scipy routine, but it's quite common to use a full Bayesian treatment here, since Gaussian likelihoods are very fast to compute." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import minimize" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimising l:\n", "\n", "Log Likelihood: -1.2038508653292523\n", "Log Likelihood: -1.2038508247602695\n", "Log Likelihood: 1.4986278795149\n", "Log Likelihood: 1.498627896855321\n", "Log Likelihood: -1.0046019674290614\n", "Log Likelihood: -1.0046021266283782\n", "Log Likelihood: 1.5969944183769633\n", "Log Likelihood: 1.5969944224813575\n", "Log Likelihood: 1.6017623612261378\n", "Log Likelihood: 1.6017623604491251\n", "Log Likelihood: 1.6019256745168962\n", "Log Likelihood: 1.6019256745553445\n", "Log Likelihood: 1.6019260806170736\n", "Log Likelihood: 1.6019260806174316\n", "Log Likelihood: 1.6019260806518476\n", "Log Likelihood: 1.6019260806518472\n", "\n", "Best l = 1.6259474735691932\n" ] } ], "source": [ "print('Optimising l:')\n", "print()\n", "l0 = np.log(0.5)\n", "result = minimize(neg_log_likelihood, l0)\n", "l_max = result.x[0]\n", "print()\n", "print(f\"Best l = {np.exp(l_max)}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have our optimised value for $l$, we can compute all of the other values of $\\Sigma$. From this, we can plug straight in to the Gaussian process update equations (the derivation of these is long, but you can see the full derivation on some loser's blog [here](http://www.concernsofadataguy.com/blog/4)), which define a multivariate Gaussian distribution from which our function is drawn\n", "\n", "$$ P(y_p|y_t) \\sim \\mathcal{N}(CT^{-1}y_t, P-CT^{-1}C^T) $$\n", "\n", "So that our predictions are given by the mean vector\n", "\n", "$$ \\mu \\approx CT^{-1}y_t $$\n", "\n", "with uncertainty summarised by the covariance matrix\n", "\n", "$$ \\Sigma \\approx P-CT^{-1}C^T$$\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# compute the parts of the matrix\n", "\n", "C = np.zeros((len(x_p), len(x_t)))\n", "for i in range(len(x_p)):\n", " for j in range(len(x_t)):\n", " C[i, j] = kernel(x_p[i], x_t[j], l_max)\n", " \n", "P = kernel_matrix(x_p, l_max)\n", "T = kernel_matrix(x_t, l_max)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# get our distribution over functions\n", "\n", "mu_p = np.dot(C, np.linalg.solve(T, y_t))\n", "cov_p = P - np.dot(C, np.dot(np.linalg.inv(T), C.T))\n", "sig_p = np.sqrt(np.diag(cov_p))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have the mean and covariance of our function, so we can show our best guess, with uncertainty, of what our approximation might look like" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAGuCAYAAACnXgEDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAACG1klEQVR4nO3deZxN9RsH8M/37rMw9j1LksieJHuRiGzZJe0o7ZJUoqL1VySRVBRZUtKCStZI9qwlW6Xsy5jt7t/fH2dmssxy587Z7r2f9+s1L8zce87DjHOf+5zn+3yFlBJERERERKQti9EBEBERERHFAibeREREREQ6YOJNRERERKQDJt5ERERERDpg4k1EREREpAOb0QHopVSpUrJq1apGh0FEREREUWzz5s0npZSlc/pazCTeVatWxaZNm4wOg4iIiIiimBDiz9y+xlYTIiIiIiIdMPEmIiIiItIBE28iIiIiIh3ETI83ERERkZp8Ph8OHz4Mt9ttdChkAJfLhUqVKsFut4f8HCbeRERERGE4fPgwihQpgqpVq0IIYXQ4pCMpJU6dOoXDhw+jWrVqIT+PrSZEREREYXC73ShZsiST7hgkhEDJkiULfLeDiTcRERFRmJh0x65wvvdMvImIiIiIdMDEm4iIiChCWa1WNGjQAHXq1EGvXr2Qnp4e9rHuvPNOLFiwAABw7733Yvfu3bk+duXKlVi3bl2Bz1G1alWcPHkyx8/XrVsXDRo0QIMGDcI6dm4ujnXq1Kn4+OOPVTt+QXBxJREREVGEiouLw7Zt2wAAAwYMwNSpU/H4449nfz0QCMBqtRb4uNOnT8/z6ytXrkRiYiKaNWtW4GPnZsWKFShVqpRqx8tycaxDhgxR/RyhYsWbiIiIKAq0bNkS+/btw8qVK3HDDTegf//+qFu3LgKBAJ588klce+21qFevHt577z0AymSOYcOGoXbt2ujUqROOHz+efaw2bdpg06ZNAIClS5eiUaNGqF+/Ptq2bYtDhw5h6tSpeOutt9CgQQOsWbMGJ06cwG233YZrr70W1157LdauXQsAOHXqFNq3b4+GDRti8ODBkFKG/Pc5P4aTJ0+iatWqAIAZM2agR48e6NChA2rUqIERI0ZkPyeUWMeMGYM33ngDALBt2zY0bdoU9erVQ/fu3XHmzJnscz/11FNo0qQJrrzySqxZsybM78qFWPEmIiIiKqRHH300u/KslgYNGmDChAkhPdbv92PJkiXo0KEDAGDDhg3YuXMnqlWrhmnTpiEpKQkbN26Ex+NB8+bN0b59e2zduhW///47duzYgWPHjqF27dq4++67LzjuiRMncN9992H16tWoVq0aTp8+jRIlSmDIkCFITEzE8OHDAQD9+/fHY489hhYtWuCvv/7CzTffjD179mDs2LFo0aIFRo8ejW+//RbTpk3L9e9www03wGq1wul04pdffsnz77tt2zZs3boVTqcTNWvWxEMPPQSXyxVSrD/++GP2ce644w5MmjQJrVu3xujRozF27Njsf3O/348NGzZg8eLFGDt2LJYtWxbS9yIvTLyJiIiIIlRGRgYaNGgAQKl433PPPVi3bh2aNGmSPV/6+++/x/bt27P7t5OTk/HHH39g9erV6NevH6xWKypUqIAbb7zxkuOvX78erVq1yj5WiRIlcoxj2bJlF/SEnzt3DikpKVi9ejW++OILAECnTp1QvHjxXP8uBWk1adu2LZKSkgAAtWvXxp9//okzZ86EFGuW5ORknD17Fq1btwYADBo0CL169cr+eo8ePQAA11xzDQ4dOhRSXPlh4k1ERERUSKFWptV2fo/3+RISErJ/L6XEpEmTcPPNN1/wmMWLF+c7Ek9KGdLYvGAwiJ9//hlxcXGXfC3ckYs2mw3BYBAALpmX7XQ6s39vtVrh9/tDjjVUWefIOr4a2ONNREREFMVuvvlmTJkyBT6fDwCwd+9epKWloVWrVpg7dy4CgQCOHDmCFStWXPLc66+/HqtWrcLBgwcBAKdPnwYAFClSBCkpKdmPa9++Pd55553sP2e9GWjVqhVmz54NAFiyZEl2D3Uoqlatis2bNwNAdrU+L6HGmiUpKQnFixfP7t/+5JNPsqvfWmHiTUREROYVDAC+DMDvATKrn1Qw9957L2rXro1GjRqhTp06GDx4MPx+P7p3744aNWqgbt26GDp0aI5JZ+nSpTFt2jT06NED9evXR58+fQAAt956KxYuXJi9YPHtt9/Gpk2bUK9ePdSuXRtTp04FADz//PNYvXo1GjVqhO+//x6VK1cOOe7hw4djypQpaNasWY4jCMON9XwzZ87Ek08+iXr16mHbtm0YPXp0yPGFQxRkdWkka9y4scxaGUtEREQmFgwC3hTA7wWCF93it1gBR4LyYbA9e/agVq1aRodBBsrpZ0AIsVlK2Tinx7PHm4iIiMzD7wEyzgIyl+p2MAC4zwE+N+BKAqxMZShysNWEiIiIjCelklCnn8496T5fwAukn1TaUIgiBN8mEhERkfHcZ5UqdkFICbiTAWEBbM78H09kMFa8iYiIyFielIIn3VmkBDLOAAGfujERaYCJNxERERnHmw54Ugt3jKzkOxhQJyYijTDxJiIiImP4PUqriBqCASX5JjIxJt5ERESkv2BQvaQ7S8BX+Op5BDl16hQaNGiABg0aoFy5cqhYsWL2n71eb57P3bRpEx5++OF8z9GsWTNVYl25ciWSkpLQsGFD1KxZE61atcI333wT0vPWrVunSgxmwMWVREREpD9vijatId5UwB6nzPuOciVLlszeIXLMmDFITEzE8OHDs7/u9/ths+Wc6jVu3BiNG+c4avoCaia9LVu2zE62t23bhm7duiEuLg5t27bN9TkrV65EYmKiam8AjMaKNxEREenL71V6u7WQNenEhL7c+g+av7Ic1UZ+i+avLMeXW/9R/Rx33nknHn/8cdxwww146qmnsGHDBjRr1gwNGzZEs2bN8PvvvwNQEtrOnTsDUJL2u+++G23atMHll1+Ot99+O/t4iYmJ2Y9v06YNevbsiauuugoDBgxA1iaMixcvxlVXXYUWLVrg4Ycfzj5uXho0aIDRo0dnbzP/9ddf47rrrkPDhg3Rrl07HDt2DIcOHcLUqVPx1ltvZe86mdPjIgkr3kRERKQfPRJjv0eZ722P0/Y8BfDl1n/w9Bc7kOFTqvz/nM3A01/sAAB0a1hR1XPt3bsXy5Ytg9Vqxblz57B69WrYbDYsW7YMo0aNwueff37Jc3777TesWLECKSkpqFmzJoYOHQq73X7BY7Zu3Ypdu3ahQoUKaN68OdauXYvGjRtj8ODBWL16NapVq4Z+/fqFHGejRo3w+uuvAwBatGiB9evXQwiB6dOn47XXXsP//vc/DBky5IJK/pkzZ3J8XKRg4k1ERET68aZeug28FtznAKsTsJjj5v7r3/2enXRnyfAF8Pp3v6ueePfq1QtWq9Jqk5ycjEGDBuGPP/6AEAI+X85jFzt16gSn0wmn04kyZcrg2LFjqFSp0gWPadKkSfbnGjRogEOHDiExMRGXX345qlWrBgDo168fpk2bFlKcWRVzADh8+DD69OmDI0eOwOv1Zh/vYqE+zqzM8dNIRERE0S8YALxp+pxLBgGfTucKwb9nc95hM7fPF0ZCQkL275977jnccMMN2LlzJ77++mu43TnPS3c6/9uAyGq1wu+/9M1RTo85P3kuqK1bt6JWrVoAgIceegjDhg3Djh078N577+UaZ6iPMysm3kRERKQPT4rSaqIXb5oyPcUEKhTLue0lt8+rJTk5GRUrKhX1GTNmqH78q666CgcOHMChQ4cAAPPmzQvpedu3b8eLL76IBx988JI4Z86cmf24IkWKICUlJfvPuT0uUjDxJiIiIu0FA4Bf5+qklEpriwk8eXNNxNkvnLQSZ7fiyZtranreESNG4Omnn0bz5s0RCKg/RSYuLg7vvvsuOnTogBYtWqBs2bJISkrK8bFr1qzJHif44IMP4u23386eaDJmzBj06tULLVu2RKlSpbKfc+utt2LhwoXZiytze1ykEIW5RRBJGjduLDdt2mR0GERERLEp46yy4FFvQgAJZTTp9d6zZ092q0Qovtz6D17/7nf8ezYDFYrF4cmba6re322E1NRUJCYmQkqJBx98EDVq1MBjjz1mdFi6yOlnQAixWUqZ46xGLq4kIiIibQX8xiTdwH9Vb1dRY85/nm4NK0ZFon2x999/HzNnzoTX60XDhg0xePBgo0MyLSbeREREpC1vSv6P0ZIvHXAkxMSmOkZ47LHHYqbCXVjs8SYiIiLtBPyAz+DJE1LqN02FKA9MvImIiEg7Zhnp50s3zYQTil1MvImIiEgbwaBxvd0Xk1JJvokMxMSbiIiItOFL13dud368aeaKh2IOE28iIiLShtkqzNJEFXiKSUy8iYiISH0+t7JpjtlwkSUZiOMEiYiISH1mq3ZnCWZOWbG71D/2uSPqHq9o+ZAedvfdd+Obb75BmTJlsHPnzgKd4ujRo3j00UexceNGOJ1OVK1aFRMmTMCVV14ZTsR4++23MWXKFDRq1AgHDx7EunXrLnnMmDFjkJiYiOHDh4d1jkjGijcRERGpK+AH/B6jo8hdlFW977zzTixdurTAz5NSonv37mjTpg3279+P3bt3Y/z48Th27FjYsbz77rtYvHgxZs+enWPSHetMl3gLIT4UQhwXQuT4lk0I0UYIkSyE2Jb5MVrvGImIiCgPZhkhmJuA19xvDAqoVatWKFGixCWf//XXX9GqVSvUrl0bFosFQgg8//zz2V9fsWIF7HY7hgwZkv25Bg0aoGXLlgCAN998E3Xq1EGdOnUwYcIEAMChQ4dQq1Yt3Hfffbj66qvRvn17ZGQoffNDhgzBgQMH0KVLF7z11ltITEzMPu64ceNQs2ZNtGvXDr///nv252fNmoUmTZqgQYMGGDx4MAKBQL7n+fjjj1GvXj3Ur18fAwcOzPM4ZmO6xBvADAAd8nnMGillg8yPF3SIiYiIiEIhZWQsYNS06i2VDxkEZACA/vPD3W43+vTpgzfeeAO7d+/GM888g+HDh2PMmDHZj9m5cyeuueaaHJ+/efNmfPTRR/jll1+wfv16vP/++9i6dSsA4I8//sCDDz6IXbt2oVixYvj8888BAFOnTkWFChWwYsWKC3ay3Lx5M+bOnYutW7fiiy++wMaNGwEAe/bswbx587B27Vps27YNVqsVs2fPzn5eTufZtWsXxo0bh+XLl+PXX3/FxIkT8z2OmZiux1tKuVoIUdXoOIiIiCgMfndkjOzzewC/F7A5VDxoVrIdvOTTQAAQFuUDQsVz5mzZsmVo1KgRmjRpAgCoV68eli5dCiFCO/dPP/2E7t27IyEhAQDQo0cPrFmzBl26dEG1atXQoEEDAMA111yDQ4cO5XmsNWvWoHv37oiPjwcAdOnSBQDw448/YvPmzbj22msBABkZGShTpkz283I6z5kzZ9CzZ0+UKlUKAFCiRAl8+umneR7HTEyXeIfoeiHErwD+BTBcSrkrpwcJIe4HcD8AVK5cWcfwiIiIYlQkVLuzeFMB26UtGmEJBpWFm3nJSsotVmjddLBz507UrVs3+89btmxBo0aNLnjM1VdfjQULFuT4fJnHmyen05n9e6vVmt0CkpecEn4pJQYNGoSXX3455PNIKS85Vn7HMRMztprkZwuAKlLK+gAmAfgytwdKKadJKRtLKRuXLl1ar/iIiIhiUzAQWb3Tfg8Q8BX+OEE/IPNJui94fCCzBUU7JUuWxPbt2wEAe/fuxRdffIG+ffte8Jgbb7wRHo8H77//fvbnNm7ciFWrVqFVq1b48ssvkZ6ejrS0NCxcuDC797ugWrVqhYULFyIjIwMpKSn4+uuvAQBt27bFggULcPz4cQDA6dOn8eeff+Z5rLZt22L+/Pk4depU9nPCOY5RIq7iLaU8d97vFwsh3hVClJJSnjQyLiIiophn1hGCefGkAPGFqHoH/Mju4S5SrmDPFdbM6nfh9OvXDytXrsTJkydRqVIljB07Fv369cNXX32FOnXqoFSpUpgzZw5Klix54emFwMKFC/Hoo4/ilVdegcvlyh4nWKNGDdx5553ZrSr33nsvGjZsmG9bSU4aNWqEPn36oEGDBqhSpUp2Al+7dm289NJLaN++PYLBIOx2OyZPnowqVarkeqyrr74azzzzDFq3bg2r1YqGDRtixowZBT6OUURetxKMktnj/Y2Usk4OXysH4JiUUgohmgBYAKUCnudfpHHjxnLTpk2axEtEREQAUo+bc9Oc/CSUAqz2gj1HSuzZvRu1rqpRuHOrlHyTMfbs2YNatWpd8DkhxGYpZeOcHm+6ircQYg6ANgBKCSEOA3gegB0ApJRTAfQEMFQI4QeQAaBvfkk3ERERaczvicykGwAyzirJd4gLDxEMAhmnkblqsnBkEJAic9ElRTvTJd5Syn75fP0dAO/oFA4RERGFIhLbTLIE/YDnHOBKCuGxmUm3Gr3hAACpvGGxiNATf4pYfHtFREREhSNlZC2qzIk3Pf+/Qw5Jtzo33WXk3i2IYeF875l4ExERUeH4MiJjdnd+3MlKcp2TgA9IP3VB0u1yWHHq1BmVku8gk+8IIqXEqVOn4HK5CvQ807WaEBERUYTxu42OQB3BgFLRtscBtjjAYlEScW+KUhG/SKWSRXD41BmcOKniYDVhYctJhHC5XKhUqVKBnsPEm4iIiMIXabO78xPwZVa1zwE2p/L7i3eizGS3WVGtbDF1z29zFm68IZkaW02IiIgofJG0U2VB+T25Jt2antMXJXcQ6BJMvImIiCh80dJmYiaec9HRM0+XYOJNRERE4Qn4VRyrR9mCAcCbanQUpAEm3kRERBQefxS3mRjNm577hBWKWEy8iYiIKDzsRdaODAK+NKOjIJUx8SYiIqKC83uVHR9JO940Vr2jDBNvIiIiKji2mWhPSvZ6Rxkm3kRERFRwbDPRhy+dO1pGESbeREREVDA+t/7zrWMVq95RhYk3ERFRpJDSHPOd2WaiL18Ge72jBLeMJyIiMiMpgYBXSbou3kHRYgVsrswPh/5xRdMW8ZFASmXCibOI0ZFQITHxJiIiMhtvGuBJzb2dIxhQHuNNA6x2wFUMsOr0ku7LMEfVPdZ40wFHIiCE0ZFQIbDVhIiIyCz8HiDtJOA+F3oPdcAHpJ9UEnU9cIt4Y8ig8qaHIhoTbyIiIqNJCbiTgfTT4W3BLiXgSVGer2U1Ohhkm4mRvNxQJ9Ix8SYiIjJSMKAkzN70wh/L79E2+eaiSmMF/RzjGOGYeBMRERnF71VaSwJe9Y4Z8AIZZ7RJvtnqYDyfCm/QyDBMvImIiIzgywAyTmszD9vvUT/59nvDa4Mhdfk9QMBvdBQUJk41odhz/oguGVQmAlgdygdXixORHrxpygJKLfk9St94XDF1judjf7Fp+NIAa5LRUVAYmHhTbPGkKi9451eYshYKCYvyAmVzGhIaEcUIT4p+E0h8GUpxwZFQuONwUaW5+NyAsyiLRRGIrSYUG6RUbrt6UnK/rSuDyqIkT4q+sRFR7HCf0y/pzuJJUdpECsOXxtndZsLRghGLiTdFv2BmQh3qSnBPqvYjuYgo9riTjRkHl1V4CAbCPwaTPPPhIsuIxMSbopuUyuKlgk4M8HsAj8b9l0QUO9zJ6owLDJcMAhlnwyso+DIKl7STNgI+LnaNQEy8Kbq5k8O/MHnTWeUhosLLOGts0p0l4FWuiQVlhtgpZ9xQJ+Iw8abopUbi7E7m2CYiCl/GWXO9gfdlFCyR9nvUnTFO6vK7lXZKihicakLRKeBTp1UkqzcyoRRXjxNRwZgt6c7iOZc5RtWe9+OCwfAq5BFs6c6jmLJqP44mu1EuyYWhraujQ51yRoeVOymV3UQLO7WGdMOKN0UfKcPvZcxJ0M9JJ0RUMBlnzJl0A/8VFPK7m+c+G1O93Ut3HsX4JXtwJNkNCeBIshvjl+zB0p1HjQ4tb2b9OaMcMfGm6ONJUZJlNfnS2XJCRPnLSmpDnaJklGAASD+V+5hBb1rMze2esmo/3L4L2zbcviCmrNpvUEQhCvj4+hRBmHhTdAn4tRmxJCWnnBBR3iIl6c4ig8rUp4srpoHYvMt3NDnn71tunzcVjhaMGOzxpujiOafd/G2/R/ngzpZEdLGspDvSqsRZrXnedCURl4GY3cOgXJILR3JIsssluQyIpoB8GYCrqNFRUAiYeFP08GVo/6LnPgckltb2HEQUWYKZleNInqkcQ5NLpJTYvHU7Dv/7L86lpCIlNRXFkpIwoEFtvLPuyAXtJi67BUNbVzcw2hDJoHKnxR4BbxJiHBNvig5S6nNrNOhXeh+5gpyIgMxe6dPqrysh1e3dtx+z5n2O2fMX4sChPy/5usViwdX1G0JWqAtXzZaoUKqo+aeanM+XzsQ7AjDxpujgTdVv9b03DbDHm2u8oJRKxSrgVapuFivgSFR+JSJtBHxK0i05R9nMjhw9huHPvoBPP1sIIQTatm6B0U89hnpX10LRIkVQpEgi/vr7Hyxa/B2+/HYpdn47A7X3r8PY999Bg0hJugHl+h8MAhYu3zMzIWOkl6tx48Zy06ZNRodBWggGgbTj+vYluoqap+qdW8VNCMAeBziK8EJMpDZfhjLjOkZeQyOR3+/H5PdnYPT41+F2ezD8oSF44N5BqFihfJ7P+/7Hlbjzgcdw8tRpjHvuKTw+bDCs1ggpYpjptSmGCSE2Sykb5/Q1VrzDJaVS7Qj6lIRHyv+qHkIAwqpUG62O/DcpoMLxpuj/4meWdpO8Km5SKgum/F5uAESkJvc5btVtcn/9fRi97xyCXzZtQYd2N2DSay/hiurVQnpu+7ZtsOPnH3H/wyMwYvRL2LjlV8z58N3ISL593EzH7Jh4F0TAp2zP6vcWbCGKxapMwrDFATaHdvHFomDAmM0DggElqXXE63/uLH6PMkUhvzcdQb+yEUZccV3CIopawaDyfynSJpfEmB+Wr0K/ex6A1+fD3A+noHePLhAFLDyULFECCz55H69PfBdPPT8OZUqXxKTXxxX4OLrLmultZXpnVvzO5CcYUBYs+DLC7yHOStK86UoC7khkAq4WLccH5sebZlziHfCHlnRn8bkBSyrgTNQ2LqJo5XNntpawn9uspJR4+X+T8OxLr+LqWjXx+Sfv48orwp9IIoTAiEcfxImTp/DGpKmoUK4cRg1/WMWINeLPAKxFjI6CcsHEOydSKpVtb7r6I5bOnwXtSuLit8II+I3dqCLoV96Q2eP0P3c4vaWeFKXtiXPIiUKXtXmWlxuUmJnX68U9w57ArHmfo3+v7pg28XUkJKhTGHn1hWdx5NhxPPPiKyhXtjTuHthPleNqxucGnEy8zYqJ9/mydj30ZWhf1fB7gLSTykIIIxK3aGCGnSS9afp//zyp4b8hzDgLJJZhvzdRKHxu5Tqj18QkCsuZM2fRY+C9WLlmHV569imMGv6wqi0hFosFH05+E8dPnMSQx0aiccP6qFentmrHV13Qr7SccH2ZKXHUAaBcXNNPA2knlERKr1uJMqgkQgVpGSCF32uOPsuAT984Aj5ldGK4ZJCLwojykzUpKOMMk26TO/Tn32h+c1esXb8Rs95/B888+YgmfdgOhwOffjAZxZKK4q4HHoPPZ/LNkriFvGnFbuIdDCi33lOPG7/Nr88NpJ/iBb4gCpN8qk2vRDZra+fCvknzpvGNHlFOggHl/1jaCXO8sac8bdryK5q264wjR4/j+4VzMKB3D03PV6pkSbz7v5ex5dcdeH3iFE3PVWhGtmFSnmIr8ZZSaSNJP60k3B4dN13JT8CntJ74Y2fb3rCZpdqdxe/RZ6tob5o6u+Ox6k10oYBfWTeRdiKz1ZBvTM3um6U/oHWnHnC5nFj7/SK0adlMl/P27NYZvbrdirGvvolde37X5ZxhkUFzvU5StthJvIMBIPWYUs0w6w+jDAIZp/lONT96bA1fUFonskGVk2VWvYkuajNM5/+JCDFl+kx07XcXal1ZA+uXfYPaV12p6/nfeWMcihZJxN0PPg6/X4ViiFaMGLVL+YqdxFsGI+OiKqXS+sIV9Dnze9SfNKMGv1vbuyfeFHXXHrDqTbHK71EKMCnHjG8zpAIJBoMY8dyLeOCJp3FL+xux8tvPUa5sGd3jKFO6FCa99hI2bN6KaR/N0v38IfO7IyPviTGxk3hHGncyE6OceEzU230+KbX7fgX82lQuWPWmWBDwKz/r6aeVZDv9tD6Tq0hVbrcb/e4eitffnoKh9wzCwtkfIjHRuB0a+9zWFS2bXYcXX5+AtDSTFsqk5BtLE2LibWbuc+ZsqzCKz23OancWrXpDvSnaHFcGeSuSokcwqKz/8KYr186sRDvthPJnv4fJdoQ6dvwEburWF/MXfo3XXngWk/83HjabsdOQhRB4+fmncfTYcUx67wNDY8mTn9d4s+Ecb7PzpCpJl6uo0ZEYz+xvQmRQGeHkULEK4/dq2/PvyzB223uivASDAOR/rYIyeN5HQGnvksH/fqWos3rtevS9eyjOnE3G3A+noM9tXY0OKVvzpk3QuUM7vDrhXQy+ayCKFy9mdEiX8nuV/zvcu8E0TJd4CyE+BNAZwHEpZZ0cvi4ATARwC4B0AHdKKbfoG6XOsmaLxxUzOhLjqDXRQ2veNMAer95FTus3GwGvcivearpLAelJXpTcQmbeZTnv16zHXfD77APk/rlLPn/x53I5B9ugYlowGMTrE9/FMy++isurVsHSz2ebctOacc+NRIMWN+G1ie/i5TGjjA7nUlnTTewuoyOhTGZ8tZ0B4B0AH+fy9Y4AamR+XAdgSuav0c2XAUACccWNjkR/Upq3t/tiwYB6VWS9Wmt86YCVd1SiVjCgvGkN+jMrwwHl/1R2lZhJLpnL/gOHMPjREfhx1U/o1e1WTJ/0BooWNecW6PXq1Eb/Xt0xcep0PDzkHpQvV9bokC7ldzPxNhHT9XhLKVcDOJ3HQ7oC+Fgq1gMoJoQor090BssafRWMsVuq3tTIuo2s1iJLvTYJ8nN8ZdTIWkiYcVbZFyDlqLJnQfpppc/Zm6ZcR/weJRGPlGlPFBP8fj9emzAZda6/ERu3/oqpb72KeTOmmjbpzjL26eHw+fwY98ZEo0PJmd/D/+cmYrrEOwQVAfx93p8PZ37uEkKI+4UQm4QQm06cOqVLcJrze5RdLgMR0HahhmAg8qa7BFWYQuLL0GdTHkD5N+bK98gUDCqLCS9eSJj188MXW4oAwWAQC79egoYt2+Op58ehY7sbsPuXlRh890BNtn9XW/XLq+LO/r3xwSdzceKkCXMNbqajn6xrch4iMfHO6X9hjq8uUsppUsrGUsrGpUuW1DgsHQX9SvIdC/+RPOciM3kozJsFI1prfCYdh0WXklK5sKedUjYFcydzYgdFpEAggIVfL0GjVu3R4/Z74PF48fkn0/HF7A9QsUJk3ch+fNhguN1uTPlgptGh5Ix3NrWT0zU5D5GYeB8GcNl5f64E4F+DYjGODCpVroyz0dt64nNH7i6eAV/4sfvS9V9I6vdE789RtMja1jz1uPKrmUdrEuXht71/4Okx41GlThP0uP0epKe78fF7b2P3hpXo0eUWo8MLS62aNXBL+7aY/P4MuN0mfN1iu4n6wrwmm3FxZX6+AjBMCDEXyqLKZCnlEYNjMo4vQ/kP5SoK2OOMjkY9waBS7Y5k3rSCL2gJBo1ZSCqlMu9VzVGIpI5gQJluw5nrFKHS0tKxet16/LBiNX5YsRo7d/8Gq9WKDu1uwISXx6Jb5w6Gz+VWwxPDBqNtl96YNe9z3DtogNHhXIjTTdQTDCr7a4S5w7iQJnsHJISYA6ANgFIAjgF4HoAdAKSUUzPHCb4DoAOUcYJ3SSk35Xfcxg3ry02rlmoVtjlYrErybYuL/PFwGWejI9FwFS1YMmvk39vmBOJLGHNuupSUyptPrTZmItJAWlo6tu3YiS2/7sDmbTuwedt27Pn9DwQCATidTrRoei063dwO/Xp2M2S7dy1JKdGoZXu4PR7s+mUlLBaTNRXY42J7LHFhSakMPQhh12eRVGGzlLJxTl8zXXYmpeyXz9clgAd1CieyBANKtdSTCljtSiJusSu/F1ZAWACzXQhy4vdER9INKJVKm0v5XuTH6L93VrtJJPyMRDufW0m6gwGjI6EYs3TnUUxZtR9Hk90ol+TC0NbV0aFOuRwfK6XE3n378dPPG/DT+g3YsHkbftu7D8HMtrWyZUrjmgb10K1TB7S8vglaNrsOcXFRdGf2IkIIDH94KG6/bxiW/LAcnW5uZ3RIF4qFdWFaCfgy20kKP/TAdBVvrcRExTsUQihJuMV6XmLuME+FPBhQFo5GU8IRSiVZSmUihdF/77hi0dWyFGmCQcCTHLlrGyiiLd15FOOX7IHb9996D5fdglEda2Un34FAAKvXrsf8hV9j4TdLcOz4CQBAqZIl0PTaRrimQb3sj/LlykbEVBI1+Xw+XF6/KWpcfjmWf/OZ0eFcKr4kYHMYHUVk8aYpRbQC5MsRVfEmjUkJSP+li/eERamM25yhV2jVFsxcMGp08qk2v0dJpPLqrTNLddOXwcTbKAE/kHEmMnZopag0ZdX+C5JuAHD7gpiyaj+ureDAhHenY9qMWTh+4iTi4+PQqX07tL+xFVo2uw5XXlE95pLsnNjtdjx0/9146vlx2LFrD+peXcvokC7kdzPxDpWUyjVZ5TsFTLxJkbXwwu8BcC4zCXcpSZgeSXjWD3i0Jh2ec/+1/1zMmx72Ig3VBbzK94IvoPryuZXbmBwJSAY6mnzpnRZfWjI2rJiDKuNWIi0tHV1uaY+BfXrilvY3Ij5ehR16o9DdA/viuXGv4/2Zs/H2ay8ZHc6F/G4A3Kk4X8GgkpNoMD2KiTflLOBTPjwpSiuK3aUs2tSq/1ejH3DTCAaUVhJHAuBIVBLbYOC/GcxmISVXvuvNm6ZsekNksHJJLhzJTL6llDi1ay3+WvYJAt4M9O/ZDaOeeBhX16ppcJTmV6pkSdzW5RZ8MvdzvDJmlLneoAQDymu71W50JOYVDGTmJNpsYsdVVJS/gFdJDFKPKa0gak5ZCPiVofNmSj61krUxTtoJJeFOO2HOvzc3WtAPk24ykaGtq8Nlt8CXlox9Cyfg4LfvIaHMZZg2ewFmT5/MpLsA7r/zdpxNTsZnX35jdCiX4jU+d1ktrxruHM3EmwrG71FG3mUl4d708Dde8aQC6Seju9Kdk2BA+Xcz68JmM74ZiEaeVCbdZCod6pRD14oe7P7waSQf2IHaHe/Ap7Nn4b5O1xsdWsRp3eJ6XHnF5Zg2Y5bRoVyK1/icSQlknNa85ZWJN4UnqyXBnawk4amZVVyfW6li58bvVdpX0k4WeJUw6UQGle8TaSdrlTyRiXy1+Ds8+9gDqFqhNHau+wG75r6CTvUrGh1WRBJC4P47b8e6XzZh5+7fjA7nQgGfORbzm42G7SXnY+JN6gj6lSpuxhmlhSLlqJKMp51UPlJPACnHlFGBnlRdfripEHgrUjs+NyvdZDrTZ85G9wH3oN7VtfDTd4tQ+6orjQ4p4g3q3wsOh8OcVe9o2StDLRlndbsTwMSbtCGlkoxnLdIM+jmxIZLwVqQ2Aj7AfdboKIgu8Nbkabjv4Sdxc9s2WP71ZyhdqqTRIUWF8xdZpqebZHJVFl7j/+NN0/WNCBNvIrpU0M9bkWrLGk/F9ioykXmfL8Ljo8agZ9fOWDTnIyQkmGgCh1FsTtUONfiugTibnIwFi75V7ZiqCHjDX58VTQJ+3dv+mHhHuaU7j6Lr5LW4bvyP6Dp5LZbuPGp0SBQp2G6inuw59XwzQ+axeu163DHkEbRsdh0+mfY27PYYHzEnLMouw/ElAGeiKods1bwpqleripmfmnAXy0CMV72zrss6F0OYeEexrO1/jyS7IQEcSXZj/JI9TL4pNLwVqR5PSuxN7yFT2/P7H+ja7y5cXrUyvpz9IVyuGJ/db7FlbqeeWe12FgEcha/+CyFwR7+eWLFmLf76+3Chj6eqWC+ueM4ZsmkfE+8oltf2v0T5ytrFkgrH71V6CIlM4uzZZHTqNRBOpwOLP5uFEiWKGx2SsYRFSbqtF+0p6EpSdm8upIF9ekJKiVnzvij0sVQVy9OrfG7Ddoxm4h3Fctr+N6/PE11ASk6fKSwpuZiSTEVKiSGPPYW/Dv+DhbM/QLWqlY0OyXj2PHZldiUpOw0XQrWqldGqeVN8PPczSDMVM2QwNu9sSqlUuw3CxDuKlUvK+dZhbp8nukSs34osLHcy+7rJVD6aNRfzvvgKLz4zAtc3aWx0OMYTAnDk0c8tBGAr/GvmoH698Psf+7Fh89ZCH0tVsZh4e1MNvS4z8Y5iWdv/ns9lt2Bo6+oGRUQRh33J4fO5OSuXTOW3vX/goRHP4sZWLTDi0QeMDsccbK7cq91ZVGg36dm1M+LiXPh4jskWWcZa4h3wG976x8Q7inWoUw6jOtZC+SQXBIDySS6M6lgLHeqUMzo0ihQBH0dOhcPgW5lEF/N4POh39wOIj4vDJ9PehtVqNTokc8ir2p3F5gQshfv3Klq0CLp37og5CxbB4zFRshv0573bdLTxnDN87ZIt/4dQJOtQpxwTbSqcgAewFL7iE1M8KWwxIVMZ98ZEbNuxC1/Pm4kK5fmaAEBJqC9eUJnrY12FrpTe0bcXPv1sIb797kf06HJLoY6lKr8bsKozPtHUfG5TVPhZ8SaivJngQhVRAn7AZ7Jd6iim7f5tL155azJu73MbOne4yehwzCOUancWFdpN2t3QEuXLlcXMOfMLfSxVxco1XueNcnLDxJuI8sY+74Ixwa1MoizBYBCDHx2BIomJeHP8GKPDMQ+rA7A5CvB4uzLruzCntFrRv2c3LPlhBU6fPlOoY6kqFnax9KYbMrM7J0y8iShvwQDHCobKJLcyibJ8+Mkc/PTzBrz+4rMoXaqk0eGYhz2MSSXhPOci/Xp2h8/nw8JvlhT6WKqK5l0spVQmmZgEE28iyh+r3qExya1MIgA4dvwEnnzuJbRq3hR33d7X6HDMxRZG60g4z7lIowZ1UaP65Ziz4MtCH0tV0Vww8KWbas0NE28iyl80X5TVYqJbmUQAMGL0S0jPyMB7E16FKOQmMFHF5sx/hGBOrDal5aQQhBDoe1sXrFizDkePHS/UsVQVrdd4KQ0fH3gxJt5ElD+2muTNZLcyiTZt+RUfz/kMjz94P666sobR4ZhLYTbEsRagLzwX/Xp2RzAYxPyFXxX6WKqRwei8zpus2g0w8SaiUMgg4Ge7Sa5MeHGn2CWlxBPPjkXpUiXx9OMPGR2OuRR2J0qbs9Ah1KpZA/Xr1sacBYsKfSxVRdtOxSasdgNMvIkoVOzzzpmUgIfVbjKPL79ZitVr1+OFUU+iaNEiRodjLlZHeG0m5z9fBf1u64b1Gzfj4KG/VDmeKqKtuOLLMGVBhIk3EYWGiXfOvGnKHQEiE/B6vRgx+iXUvupK3Duov9HhmE9h53ELUeg+bwDo06MrAGDeFyaqekfbWEETVrsB7lwZVQKBANau34jPvvwam7dtx5mzyTibfA7pGRmoVuUyXH1VTVx9VU20btEU1zdpDEth3vVT7Al4leouF2n9x6S3Mil2vTt9JvYdOIjFC2bBZuNL/AUK22aSxeoodD901SqX4fom12DO54sw0kztQNGyU7HfY9rF7vxfGQXOnk3GuDcmYtb8L3D02HG4XC40vbYR6tS6CsWLJcHpdGDfgUP4af0GfPrZQgBA5csqom+PrhjQuwfq1alt8N+AIoKUyotNQTadiHasdpOJnD2bjBdefQs33dAKHdrdYHQ45mNzqlM4sDoAFP4Nd7+e3fDwiOew+7e9qH3VlYWPSw1+tyq7dBrOxAURJt4RTEqJzxZ+jUdGjsbxEyfRrVMH9O5+Kzrd3A6JiQk5Pufs2WR8890yzFnwJd6cPA2vTXwX7dq0xIhHHkC7G1rpOnJq6c6jmLJqP44mu1EuyYWhraujQ51yup2fwhDwMPE+H7eGJxN5c/J7OHP2LF574VmOD8yJCnO4leMUfoElAPTu3gWPjnwecxZ8iRefHaHKMQstGvq8A35Tj0cUMka2Nm7csL7ctGqp0WGo5tjxE7jrgcew5IfluKZBPUyb+BoaNahXoGOcOn0aH3w8BxOmTMeRo8fQsF4dPDfiMXTr3EHzi/bSnUcxfskeuH3/VQtddgtGdazF5NvMrA4ggbvfAVAW7mScNToKIgDAyVOnUK1eU3RsdyPmz3zP6HDMRwggsax6rXJpJ1UZv3dT1z44+Off+GPrWvO8WYovGdkFFneysq+CgURShc1SysY5fY1NvhHowME/0bx9V6xa+zMmvPIC1v/4TYGTbgAoWaIERjz6IA5uX4/pk95Aalo6etx+D65pdTO+WfoDtHxTNmXV/guSbgBw+4KYsmq/ZuckFQR9SssJmfpWJsWeV9+ajPT0DIwdNdzoUMzJ6lB3fYpK00363tYV+w8ewqYtv6pyPFVE8vbxwaBSFDExJt4RZvvO3Wh+c1ecOZuM5V9/hkeG3pu9gGbpzqPoOnktrhv/I7pOXoulO4+GdEyn04l77uiP3RtWYsaUCUg+l4Jb+wxC07ad8d2ylZok4EeTc54Xmtvno0W43yPTkJLTTQDlNmY0bjZBEenI0WN45/0ZGNC7B2rV5GY5OVK7b1mlxLvHrbfAbrdjrpmmm0Ryu4kv3fTFIfZ4R5C16zegU+87UCQxAT9+Nf+CxRgXt24cSXZj/JI9ABBy64bNZsOg/r3Rv1d3zPx0Pl58fQI63NYfzZtei7FPD8eNrVuodiusXJILR3JIsssl5b/iPDU1DX8d/gf/HjmKlNQ0pKSmIjU1DRaLBQ6HHQ6HA4kJ8ShTuhTKlCqFsmVKo0iRRFXiLgw1vkem4Peo1uMYsVjtJhMZ98ZE+P1+PP/U40aHYl5Wla9ZKl0Dixcvho433YB5X3yF1198zhzTxrLGCpohloIyebUbYOIdMX7/Yx869x6EsqVL4Ycv56LyZZUu+HperRsFTersdjvuHTQAA/v2xIefzMW4/72Ndl37oGG9OnjswfvRp0cXOByFe7c/tHX1HHu8h7aunv3nlJRUbPl1B37duQs7dv2G7bv2YN+BQzh95kyBz1e2TGlcecXluLL65ahZo7ry+yuqo3q1KoX+u4RKze+RoWK90hvwmXrhDsWWP/86jGkzZuPu2/ui+uVVjQ7HnGxO9ZPIrHneKlwP+93WDV8t/h5r1v2C1i2uVyE4FUTiWEG/17QjBM/HxZUR4PTpM2ja7lacTU7GhuWLUbXKZZc85rrxPyKn76QA8MuotoU6v9vtxsdzFmDClPex5/c/UL5cWdzRtyf69OiCBvXqhF0FP3+qSZlEG26uJOFIPowNm7fil81bsfu3vQhmDvMvVbIE6tauhZo1qqPKZRVRpXIlVCxfHklFiyAxIQEJCfEAlM0jvF4fzqWk4MTJ0zh24gT+PXIUf+w/iL37D+D3P/bjxMlT2TFYrVZcecXlyozzWldm/loTNapXg91e+E0Szqfl90h3RcrF7jzvjLMRUVWh2DD0sZH4cNZc7Nu6FpdVqmh0OObkKgo4cp70VSjuc6rc/UpLS0eZK+piYJ+emDrhVRUCU4E9DogrZnQUBWOia3NeiyuZeJucz+dDhx4D8NP6DVj+9Xw0b9okx8d1nbw2x9aN8kkuLHqwuSqxBINBfL98FSa99yG+X74Kfr8fV15xObrecjOaXnsNml7bCBXK51+59Xq9OHDoT2zauh0bt2zDxi3bsHX7LrjdSvwlihdHk2saoOm1jdDkmoZoWK8OypYprVqby5kzZ/HHgYP4/Y/9+G3vPuz67Xfs2rMX+w8eyu5nt9vtuLxqZVSvVhXVq1VBtSqVUbF8OZQvVxbly5ZB8eJJSCpaNN8NKoLBIDweD9xuD/pMWY0jJ88i4HUj4M1A0OtGwJOBRKsfdzSpCK/XB6/PCyklHHYHnE4HXC4nypYujfLlyqBCuXK4rFIF1d8QhCWuOGBXYSOKSBMMAmnHTd9DSLHhn3+P4PL61+OuAX3Mk7CZUWIZwGJV/7g+N5BR8DuwOel391D8sGI1juzdZo5rvMWq/LtFCimB1GOmuTYz8UbkJt5DHxuJqR9+jI/fexsD+/bM9XF6j+c7eeoUvvhqCeZ98RXW/PwLfD7ldlv5cmVRsXw5pb+6dEkIIeB2e5DhduP0mbM4+OdfOPzPkewENz4+Do3q18W1jRqgSaMGuLZRA1xerYohY5XS09MzE/G92LXnd/yx/yD2HzqE/Qf/RGpqzlWN+Pg4xLlcsNvtsNmsEELA4/HCnZlse73qLlJxuVxoUPdqNG5YD9c3uQad2rdDUlJRVc8REkeCUkWKNZ5UwJNidBREAIBHR47GO9M+wh9b1qJa1cpGh2NOWo5ADQaA1OOqHOqrxd+ha7+78O1nn+CW9ia5A5pQGrBGSEeyN025A2ESTLwRmYn3nAUL0f+eBzHysWF4ecyofB9v1IY0brcb23bswi+btmLbjp04dvwkjp84ieMnT0JKiThXHFwuJ4olFUXVypehWpXLcHnVKmhUvy5q1axh+m2NpZQ4feYMjhw9jiNHj+Hfo8dwNvkczqWkIDn5HDLcbvj9Afj8PgSDEi6nEy6X86JfXYiPi8OB5ABW7k9Gst+KsiWScEfLmuh8TTXEuVxwOOyw2+0QQsDn88Hr9SHD7cax40q7zD9HjmLXnt+xadt2bN62HampaXA4HOjQrg16d++C27rcApdLpyq01Q4klNLnXGaSeiIieggp+h07fgLV6l2H3t27YMaUCeoe3GIFhCWzeiiVBDNSOYsATg0X16ceV+Xfx+PxoFyNBujcoR0+mTZJhcBUoPW/nZpUmquuFibeiLzE++/D/6Bes3aoVfMKrF6y0PTJKekrEAhg45ZtmL/wa8xf+DX++fcIKlYoj6cfH4Z77+gPp1OHqSOJZSNz1Xu4/B4g/bTRURABAEY+Pw6vTXwXv21ajSuvqJ7/E0IlBBBf6sJKZ8YZpa0iEmldtU0/rdpi63uHPYF5C7/C8X3bERdngoWNNicQX8LoKPIX8CmJt4lwA50IEwwGcdcDj8Hn9+HjqW8bnnRH/OzpKGS1WtH02mvw5vgx+GvXRvzw5VxcXrUyhg1/Blc0bIYZs+dpugESgNib580RgmQSp0+fweTpM9CnRxd1k24AsMdfmqg6k7TpkdaaxaZ9q4RK87wBZTOd1NQ0fPvdj6ods1ACXtP0TOfJZ+wulQXFxNuEJr33IX5c9RPeGj8WV1SvZmgsWb3jR5LdkPhv9jSTb/OwWCxod0MrrFr8BX74ci4qVaiAux54DD1uvwcnT53K/wDhiqXEOxjgCEEyjbff+wCpqWkY9cTD6h5YWJT2gotZLICrmLrn0oMeC8Ct6i2EvKFVc5QtUxpzPzfJZjqRsGGalBF3N4aJt8ns/m0vnnp+HDp3aId7B/U3Ohxu7R5BhBBod0MrrP1+Ed54aTQWf78cda9vi6XLVmhzQrNfkNUUYRUVil7nzqVg4pQP0K1zB9S9upa6B3cWyX1MqM0ROf2+WWx6JN7qVbytVit6d78V33y3DOfOmWQRt9/kSa3fA8hg/o8zESbeJhIMBnH/I08iMSEe0yf9z5DJHheL1a3dI5nFYsETDw3BhuXfomSJ4uh42wC8M+1D9U8U8Cnj9WKBSWbDEk1+fwbOJifj2eGPqntgqwNwxOf9GGcRVSu8mrJY9YlVCKWlRSV9b+sKj8eDL781yZo0s28f74+8azMTbxOZMXse1q7fiNdffA5ly5Q2OhwAuW/hHsrW7mSs+nWvxsYVi9G108146MlnMe71ier3fcdC1dvnjuypDhQ10tLS8ebk99DxphtxTcN66h48pxaTnDgipOqt0pbuIVExwb++SWNUqVzJPO0mQb95r39SRmQLIBNvkzh1+jRGjH4JLa5vgkH9exsdTrahravDZb/wx+Tird3JvOLi4vDZzGm4vc9tePalV/Hksy+om3wHIu+iV2BsMyGTmDZjFk6eOo1nn3xE3QNbrEorSSjsLlUrvJqx6TgVRMV2EyEE+vboih9WrNZ2jU5BmDW59bsjY/HnRZh4m8RTo8ch+VwKprz5CiwmGtHWoU45jOpYC+WTXBBQdsLUalMe0obdbsfMqRPx4H134n/vvIcRz72o3sFNNDdVE8FgbFT1yfTcbjdef3sKbmjZHM2uu1bdgxe0Omz2Xm9hCf2NhBpUbmnpe1tX+P1+fL5osarHDZtZCywR2gIYAW9bo9/a9RvwwSdz8OTDQ1Gn9lVGh3OJDnXKMdGOcBaLBZNeH4dgUOKNSVNRq2YN3D2wX+EPnNXnbaI3i6rypUdkRYWiz4efzMWRo8cwS4vNVQq6CNEep+zgatYWBD3bTAAl8RZCtWtF/bpX46orr8Ccz7/E4LsHqnLMQjFjn3cwaN5KfD6i9NUycgQCATw4fBQqX1YRz498wuhwKIoJIfD2ay/iphtaYchjI7Hqp5/VObBZqyFqiNCKCkUXr9eLVydOxvVNrsENrZqre3BhCS9RdSSoG4ea7AZsPmNRr+othEC/nt2weu16/PPvEdWOGzYZNF/yHYGLKrOYLvEWQnQQQvwuhNgnhBiZw9fbCCGShRDbMj9GGxGnWj6aNRe/7tiNN14cjYSEfFaUExWSzWbD/Bnv4fKqlXHbwHux/8Chwh80WlsxAj5uD0+mMGve5/jr73/w7JOPqj/tKtzqsD1eSdrNRlhU7bkOmQbtJlJKzF/4tarHDZvZCiwRNrv7fKbaMl4IYQWwF8BNAA4D2Aign5Ry93mPaQNguJSyc0GObcYt48+dS8GV17RAjerVsHrJQlOMD6TYsG//QVzXtjMqlC+LTSuXFG6LeYsNSDTHFB5VuZMBLxdWkrH8fj9qXdsaRYskYtOqpeq/TsQVD3+jGU+q0nJiJnaX8nfSmy8DyDir6iGvaXUzrFYrNqwwQa+31QEklDQ6CkUwAKQeNzqKPEXSlvFNAOyTUh6QUnoBzAXQ1eCYNPPym5Nw7PgJvDluDJNu0tUV1ath1vuTsHP3b3h+/BuFO5iZx02FKwJ3Q6PoNH/hV9h34KA21W4hCtcPbY/PfcMdo+g5zeR8GlTZ+/Xsho1btmHf/oOqH7vAzLR9fIS3AJot8a4I4O/z/nw483MXu14I8asQYokQ4urcDiaEuF8IsUkIsemEWcbyZDr059946933MbBvT1x7TQOjw6EY1PGmG3HfoAF4/e0pWPfLxsIdLNraTfzuiNsNjaJPMBjEuDfextW1aqJrp5vVP4HNWbjE2WLRZ3fIUIXbr64Gi1X11pve3W8FAMz74itVjxs2syxmNPtumvkwW+Kd0xXg4rdYWwBUkVLWBzAJwJe5HUxKOU1K2VhK2bh0SZPcIsn01PMvwWIRGD/6kjZ2It38b9zzqFypIgYNeRRpaYVoqzDLBVktEV5Roejw5TdLsfu3vXhm+MPajJlVI2k20yLLwr6RKCyV+7wrX1YJLa5vgjmff6nqccNmhj7vYCDix9iaLfE+DOCy8/5cCcC/5z9ASnlOSpma+fvFAOxCiFL6hVh46zduxvyFX2PEIw+gUsUKRodDMaxIkUR89O5b2HfgIJ4eOz78A0VTxTuCx1RR9JBS4qXXJ6BG9cvRu3sX9U8ghDqJt9VuzGLGnBgxzeR8GmxR3++2bti153fs2LVH9WMXmBkmm0RBUcRsifdGADWEENWEEA4AfQFccI9FCFFOZDa6CSGaQPk7mKuPJA9SSjz53IsoW6Y0hj80VJ+TmnHlOZlGm5bN8MjQezHpvQ+xdv2G8A4SDERPn3cEj6mi6LHkh+XYun0nnn58GKxWq/onsDrUqw47TDCRy2I1rs0kOwb1E++e3TrDarWaYwt5M6znifA2E8BkibeU0g9gGIDvAOwBMF9KuUsIMUQIMSTzYT0B7BRC/ArgbQB9pZlGs+Tj6yXf46efN2DMyCeQmKjhLTqrA3AWARJKA0XKAolllZXejgQm4nSJcc+NRMUK5fHY02MQDIbZ2xwtVeIoqKhQZJNS4sXXJqDyZRVxe5/btDmJmtVZe5yS+BrJDL3mGlS8y5QuhbatW2Du54tgilTHyMQ3CtpMAJMl3oDSPiKlvFJKWV1KOS7zc1OllFMzf/+OlPJqKWV9KWVTKeU6YyMOnd/vx8gx43HlFZfjnjtU2DUwN454ZeyPMxGwZm5OarEoY5ZcRYHEMpkJuMlWo5NhEhLiMX70SGzcsg2ffrYwvIOYof+vsAL+qLiwU2RbsXot1m/cjJGPDoPdrn4yBwCwqlwdNrrNw+jzA5ossASU6SYHDv2JXzZtUf3YBWZkgSVKiiJh/4RkztymApgxex72/P4HXn5+lHYXU3sc4ErK+zFCKAl4Qmnjb82Radze5zY0blgfI8eMD2+hpRn6/wrLx7ndZLyXXp+I8uXK4q7b+2hzAiHUr87aDSzmWGyaVJvDokEc3Tt3hMvlwqx5n6t+7AIzcqxgFLSZACEm3kKI4kKIoUKIz4UQfwshPAC8mTtIbhRCTBBCtNA41oiWnp6O51/+H65vcg2639pRm5PYXUBcsdAfb7EC8SWURJ3V75hnsVjw5vgx+OffI3hj0pSCH0AGlYpxJIuSCztFrrXrN2DFmrUY8cgDcLk0ap9Qs787i5GjBcPdAEgLFpvqh0xKKoqut7TH3M8Xwes1uMAhpTF3BaOkzQTIJ/EWQlQVQnwEZbLIc1DG/U0H8DiAwQBeBPALgKYAVmRu9T4wa/Ej/WfCu9Px75GjeO2FZ7XZLMfqCH+3Lke8Uv02y8p0MkzLZtehV7db8drEd/HPv0cKfoBIbjfxe4xfOEQxb9wbb6N0qZK4b9AA7U6iVXXYkajNcfNj1KY5OdHodXRg3544dfoMli5bocnxC8SIAkWUtJkA+Ve8dwDwA2gnpawgpewhpRwrpZwspZwupXxDSjlMStkUQDkAbwIYmflBmU6cPIVXJryDLre0R4vrr9PmJK6ihXu+xfpfXzjFtFfHPoNAIBjejpaRvMAyii7sFJk2b92OJT8sx+MP3o+EBA0nhajd3519XJv+7Ys2539rmcxAozc17W9sjTKlS+HjOQs0OX6BGDE+NpJfWy6SX+JdU0p5n5RybX4HklKeklK+J6W8GsBH6oQXHV56fQLS0tLxyphntDmBzanef3ZnESC+pPEr1Mkw1apWxn2D+mPmnM/w51+HC/ZkM20rXBBSRtWFnSLTuDcmolhSEh64907tTqJFf/f59N5Qx26CUYbn02iBpd1uR7+e3fD10h9w5sxZ1Y9fIAGfst+BXoKBqNorIs+fDinlv3l9PY/nHQ0vnOizb/9BvDt9Ju69oz9q1ayhzUmchax2X8zmAOJLmWOVOBlixCMPQAiBV956p2BPlDIyL5B+D7eIJ0Nt/XUHFn6zBI8+cC+KFi2i3Ym06O8+n82pSZ9zjixWc/V3Z9Hojc3APj3h9Xrx2ZffaHL8AtGzrTDK1t6E/LZMCPGCECLH/01CiJJCiM/UCyt6PPPiK3A47Bjz9BPanMAep81tNotFWagZX4LV7xh0WaWKuGtAH3w4a27Be70jsXLMaSZksNHjX0fxYsXw6ND7tD2RHtM/9Kp6m7U4pNEbj0YN6qJWzRr4ZJ4J2k30TIZ9MZp4A3gEwAYhRJ3zPymE6A5gN4CGagYWDTZs3or5C7/GE8MGo3y5suqfQAilNURLNqey8JJzv2POyMeGIRAI4LWJ7xbsiZFW8Q4GIy9miirrN27GN0uX4cmHhyIpSeU7mBfTYxG9PU6fjdrM1maSRaM3N0II3NG3J376eQMOHPxTk3OETK/xsVF4fS7I/4z6AM4C2CSEGCmEKCWEmA3gcwBfZH6dMkkpMWL0SyhdqiSefPgBbU6i125h58/9NmuFgVRXrWpl3NGvJ6bNmI2jx46H/kS9+/8Ky++OzL50ihqjx72O0qVK4qHBd2t7IiH0SbyF0L7qbXOa926shv/GA3r3gBACs+YbPNNbBvUZ7+ePvkXvISfeUspDUsobAYwA8DyAfwC0AHCzlHKolDJNoxgj0uLvf8Sqn37G8yMfR5EiGk0Kseu8iMViVdpPsirgZr3okWpGPf4wvF4v/jdpasGeGEljBTnNhAy0eu16/LBiNZ569EEkJmp8TbfY9btz6UjQtuqt9yLOgtBogSWgtAHe2Ko5Zsyej6DRBQ492gojsXUxHwX6yRBCJAKoB8AJ4DQAFwDOn7uI3+/HiNEvoUb1y3H/nbdrcxKrw7gRSlbbf1vPx5dQLoD5VR+E+O+DIsYV1auhX89umPLhxzh7Njn0J0bKxTLKVstTZJFS4rlxr6Fc2TIYes8d2p/QpuNeDUJoN57WYsDYwoLSsJf+7oF9cfDPv7ByzTrNzhESra/zwWDkvJYUQEEWV94AYCeADgA6A6gGYC6ABUKIWUKIMHdviT4zP52P3b/txcvPP63t1vBmYHMqSXh8CSURL1IOSCyb+VFG+bVIuQs/ipZXfo0voWzeo9cKeArLE8MGIy0tHR98Mif0J0VKMstqNxno+x9XYfXa9Rj1xEOIj9ehX1nvTdLs8drcGdV6bZMaNHxd6965I4olJRXsmqyFgFfbtsIom2aSRcgQexuFEAEoifYwKeWZ8z7fBsCHABxSykoaxKiKxg3ry02rlmp+nvT0dNRo1AKVK1XEuh++0maXSiGUhDaaqse+DMB9LjZGugmL8sbJ746YnRJbdeyOvw//i33b1sFqDfGFNKG0uTa2yEnayajZhpgiSyAQwDWtb8a5c6nYs3EVnE4dKrhFyun/uuFNB9wFuFuWH6sdSCil3vG04ssAMs5qdvhhw0dh+sdzcOT3rShevJhm58lXXDHtCoHppyO24i2SKmyWUjbO6WsFaTXpI6UccH7SDQBSypVQ2k9MMFjSeJpvDQ8oVeZoSroB5T9uQmnz3z4srKwpMVmLVZ2JEfG9fGTIvTj019/4avH3oT/J7NWKgI9JNxlm1rzP8euO3Rg/eqQ+SbfFZsy1Ru27mpFQ7QaUfnoN3TOwHzweD+Ys+FLT8+RLq+t8FE4zyRJyxTvS6VHxPnHyFKo3uB43tGyGRXNmaHei+BLRnaB6UgFPitFRqM+VpLwIXSwYUCqvJq72+/1+VG9wPapVqYyV34a4mt7mVH5Wzcp9DvByTTjpLyMjA1de0wLlypTBL8u/hcWix+i9OKU6aQS1qr9mv6ZcLOWophOTGrVsDyEENq/+TrNz5EtYgCIajEtW+06JzsKueAshWhb4ZEIkCSHqFvR50eC5l15DenqGdlvDA0q/XDQn3YBSBc4pQY1kjvjc/04Wq5KUm5jNZsOw++7Cqp9+xrbtO0N7ktm3jzd7RZ6i1ttTP8Dhf47g9Ref0yfpBvTv7z6fPU6d1y1HhM1y0LjqfffAvtjy647Qr8lakBotgIzi63N+/+PnCyHWCiHuzm/xpBCiuRBiEoA/AVyvWoQRYtv2nZg2YxYevO9O7baGBwCbCbfH1YIryZxbAYdDWABHPrdH7S7T/33vvaM/4uPj8PZ7H4T2BCnNe/H0eyKmv56iy8lTpzD+zUno3KEd2rRspt+JjUy8AcBVrHALLW1OfaeyqEHjNS4DevWA0+k0fpGl2ol3FLeZAPkn3pcD+ArK3O7jQohdQoh5QojJQog3hRAfCiFWCCHOAlgOoBKAdlLKadqGbS5SSjwycjRKFC+GMSM12ho+i1l36tKCq1h0VPddRYFQqlrOJFPPRi9evBju6NsLn372JU6cPBXak8y6MIZbxJNBXnj1LaSmpml7Z/RiwmL8QmeLJfw7exar8noQaTSueBcvXgw9bu2I2fMXIiPDwAlNal/no3xTszyzASllhpTyVQBVAXQEsAhAMSgb53QCUAtKhftJAJdJKbtLKTdpGbAZLfjyG6xeux7jnhup7epii834i6eehCh8lcRoVkfoK74tFsCp8XbRhfTQ4Lvh8Xjw0ay5oT3BjIm3lOaMi6Le9p27Mfn9GRh810BcXaumfifWcKZ0gdicBW8jFAKIKx5a8cJsdPh3H3zXQJw5exbzvvhK83PlKugHAn71jmfWO6UqyXNxpRCiMoAjUsqIX/qv1eLKjIwMXHVtKxQvloTNq74LfdRaOBwJSvU01vg9ylihSBTOSL2MM4DPvBeelh264djxk/h985rQJvfElzTXLWKNx3wR5URKida39MDu3/Zi7+afUKKEjltfOItot5lNQUmpXM9DbSXQclydHjReYCmlRJ2mNyAhPh4bVizW7Dz5UutnLBgEUo8V/jgGK8w4wYMAGgKAEGK5EOIqtYOLdK9NfBd//f0PJr7yorZJNxA7/d0Xszkj88LriA/vDkV+/eAGu2/QAPyx/wBW/fRzaE8wW/WCm+aQAT79bCHWrPsFLz8/St+kGzBPxRtQKtjxJUK7pjsSIvPafz6NN4gTQuCBewdh45Zt2Lh5m6bnypNaPdn+6L8+55d4ZwDIui/UBkAMlltzt3fffrz85jvoe1tXtG6h8XpSYTFX1VBvLnP3P+fInhDe86w2Uy+07Nm1E5KSiuL9mbNDe4KZ2jqidAtiMrdz51Iw/NkXcG2jBrjnjn76B2D0wsqLCaFUsnObyW2xKe0l0XCHV4edmQf26YnExAS8+8EMzc+VK79HnV0sY+D6nN9PxFYAE4UQP2T++SEhxJFcHiullE+pF5q5SSkx5NGRcLmceOvlsdqfMBoWGRZGVr93eoiL+oxmdxWuH99RxLTtJvHx8bi9dw9M/3gOJp0+k3/1Lqv/zwzrE7iokgzwwqtv4tjxE/hqzgz9xgdmsdrNu0mXM1F5bQt4letEMJDZBx5m0cKMrHbN77IVLVoEA/v0xEez5+GNl0ajZAmDZp37Mwr3vYuRwkh+V4D7APwFoCsACaAtgF55fMSMj+d8hhVr1uLVMc+gXNky2p8w1hNvQKn4R8ptx8LOm7XaTP13vW/QAHg8HnwyL8TNdMzSbsI2E9LZ5q3bMWHKdNwzsB+uvaaB/gGYqc0kJ1Z75vqlJKUFJZqSbkC3uw1D77kDbrcbM2bP1+V8OSpssShGCiMh71wphAgCaCql3KBtSNpQc3HlyVOncFXjVqhZozrWLP1S+wqGEEBiWfNWLfQUDAJpx809akit3dWCASDthGn/rk1uuAXpGRnY8fPy/BdZWh1AQkl9AsuN3xs5d0woKni9XjRu0xEnT53GrvUrtJ16lZtIX5wY6aRUFljqoFXH7vj3yDHs3fKT/ndWsiSWCb8tNPWEcucjChRmceX5qkFpPYl5w599EcnnUvDehNf0+eG2Oph0Z7FYzL97mVoVG4vV1C+Y9w0agF17fsf6jZvzf3DAq07/X2HESDWFzGPcGxOxY9cevDfhVWOSbsB8/d2xRgjd7jo8cM8g7D94CEuXrdDlfDkK9zrr90RN0p2fkLNGKeWf0TBWsLC+WfoDZn46HyMeeQB1aus05IVtJhdyJOiyYCUsVoe63y9HomnfdPW9rSsSEuIxfeanoT3ByNXqZt5Fk6LStu07Mf5/kzCwb0/c2rG9MUEIS+QtSo9GOr1e9ehyCypWKI8333lPl/PlKNx2kxgqjBSoXCuEsAghDgghrj7/91oFZzYnTp7CvQ8NR706tTH6qcf0O3GsjhHMjRDmXe2udn+iiaveRYokonf3WzH/y6+Rnh7CRdPIxaK+DNO27FD08fl8uPOBR1GqZAlM0GPxfW7M3t8dK3T6PjgcDjwy5B78uOonbP11hy7nvETQDwQKWKONkUWVWQraJyGg7GLpvOj3UU9KicGPjMCZs8n45L234XTq9Ne22FixyInNab47ARaNxgCauOp9R99eSE1Nw5ffhrB+IuBV+taNwEWVpKPnXnoNv+7YjSlvvqz/zO7zsc3EHDTeOv589995O4oUScQbk6bqds5LFPR660uPqcJIBO7BaoyP53yGhd8swUvPjkC9OrX1O7HZkkszcRY1V0Ja0K2QQ2Xiqner5k1RpXIlzPz0s9CeYEQCHPCrt7kDUT4Wf/8jXp0wGYPvGohunTsaGwwr3uag4/chKako7hs0APO++Ap//X1Yt/NeoKBtfTFWGGHiHYJDf/6Nh0Y8i5bNrsPjwwbre3Im3rkz08g9YQHsGiXegGmr3haLBQP79MSylWvwz7+5jfg/jxEX2BjqHSRj/X34H9wx+GHUr1sbb708xuhwWPE2CyF0XZf0yJB7AAATp36g2zkvEAyE3loYQ4sqszDxzkdGRgZuG3gvhBCYOWWi9tvCn08IXjjz4yiiJL1Gc8RrmxibuOp9R7+eCAaDmD3/i/wfHE7/X2FIGXPVFDKGz+dDv3segMfrxfwZ7yEuzuD/r2beOCcW6biBWOXLKqFPjy6YNmMWzp5N1u28F/CkhNY+4k3TPhaTMUHGYl5SSjw4fBS2/LoDs6ZNQrWqlfUNwMILZ74sFuM3XBAi/O3hC8KkVe8a1S/H9U2uwcw5nyGkfQH0TIR96YA0eIwhxYSnx4zH2vUbMW3ia7jyiupGh2PeyU+xSsc+bwB4YthgpKamYdqMWbqeN1vQn39S7cuIqUWVWZh452HaR7Pw0ax5eG7Eo8aMg7Kx2h0So8cL2uOUNwBas1i1bWcphEH9emP3b3uxZdsOLN15FF0nr8V1439E18lrsXTnRZtH6DnWz8s2E9LetI9m4X/vvIcH77sT/Xp2NzocBe+WmovO/faNGtRD29Yt8ObkaaFNndKCNzX3/RuCQcB9Tt94TIKJdy5+2bQFD414Fh3a3YDnRz5hTBBW9neHRAjAWcS48+u5oY8j0RytNRfp3f1WOJ1OvDh5JsYv2YMjyW5IAEeS3Ri/ZM+FyXcwoOwiqbUY7B0k/X23bCUeeOJpdLzpRkx45QWjw/kPE29z0bniDQDPj3wCx46fwNQPP9H93ACUVhNPLsm1+2zM3o003yu4Cew/cAhd+92FShXKY/b77+jb151Fx92uooLdZcxCVHucvuMezdBak4PixYuhS8f2WPLtN0h3X5hUu31BTFm1/8In+HTo64vB3kHS1/adu9HrzvtRp3ZNzPtoKmw2k7R3CIuuPcUUAotF9zuzLZtdh7atW+DVCZORlmZQ1duXcWmhJUZbTLIUNPEOAhgL4N+Lfh81jh0/gZt79IfP58fiBZ8YN4OV28QXnNOATXWMqLQ7Ekw5231g39vgTU/BuYOXbtxwNPmi9hK/R9uZ3gF/TF/YSXsHDv6JTr0HokhiIr6Z9zGKFNHxzld+WLQxJwPeDI15+gkcP3ESUz6Yqfu5s2WcAdJPK7+6z8Vsi0mWkBNvIURnAEJKOVZKeVQqxkopj+b75Ahx7lwKOt42AEeOHcPiBZ/gqitrGBcMbxMWnNWm3SztnOhd7c4ihL7tLSG6uW0b2OMScWr3uku+Vi7poo2FpNR2zJ8eFXWKWfsPHEKbzrchPd2NxZ99gkoVKxgd0oWYeJuTAe0mLa6/Du3atMRrE981ruotM3em9LmVO5Ex2mKSpSAV70UA/hFCvCqEuEqrgIySnp6Obv3vxo7dv+HzT6bjusaNjA2I87vDo9d4QcP7yuNNN7XA4XCgfYeOOPvHFgQ8/00ucdktGNo6hykPXo12KwsGOUKQNPPH/gNo3ek2pKdnYPnX81G/7tVGh3QpAxI8CoFBb4jGjhqOEydP4d3pMww5P12oIBlKdQDTAPQGsEsI8bMQ4j4hhAH399V15sxZ3NStL1at/RkzpkxAh3Y3GBuQsLBiES6LRZ+E2Khq9/lcScaePwdPDxmAoN8L+c92CADlk1wY1bEWOtQpd+mDpUYJsjfE+bFEBbT7t71o06knPB4Pln/9mTmTboB3TM3KoDdEza67Fu1vbI3XJr6LlJRUQ2Kg/4SceEspD0kpn5dSVgNwE4B9AN4CcEQI8YkQwuBsNTxHjh5D6063YeOWXzHvo6kY0LuH0SFxjGBhOeK1feExS6uHzWG6TXWaXXctqla+DGVO78Qvo9pi0YPNc066s6i9ADIYYLWbNLHkh+W4/qZbEQgEsOKbBahXp7bRIeXMYtNnvCkVnMViWMHmpWefwslTp/HKW+8Ycn76T1j/O6WUy6WUAwFcCWAzgAEAlgkhDgohHhNCmOseeC727T+Ilh2648ChP/Ht/I/Rs1tno0NScIxg4bmStFucao83vtqdxVnUVOMFhRDo36s7flixGkePHc//CUGVF0F6zrHaTaqSUuLNd95D59534PIqlbFh+beoU9vE3ZacZmJuBrUIXntNAwzo3QP/e+c9/PnXYUNiIEVYr9hCiNZCiBkAfgdQB8BkAO0BfAZl0snHagWolYVfL8E1bTrgzNlkLFs0Dzfd2NrokP7D24SFZ7VpU5W2WI3t7b6YxQK4zNXtNaB3DwSDQcz74qvQnqBW1TvgUxbvEKkkOfkc7hj8MJ54Ziy6deqAn75bhMqXVTI6rLzx9cPcDGwjffn5pyGEwMgx4wyLgQo21aSKEGK0EGI/gOUALgNwP4DyUsqHpJQ/SilHABgEoKs24Raez+fD8GfGosft96DmFdWxedVSNL32GqPD+g/nr6rHmaj+i5CWlfRw2eNMtRi39lVXomG9Opg9/4vQnpC12r2wPCmFPwZRphWr16Je87b49LOFGPP0E/js42lISDDnzrEXYOJtbgZ+fy6rVBHDHxqCuZ8vws8bNhkWR6wrSMX7AID7AHwK4AopZVsp5Rwp5cX3iXcB2KBWgGravnM3Wnbonr2175qlC1G1ymVGh3Uh9nerK66YeomyyRLcC7iKmaf9BUrVe+OWbdi7b3/+DwYK3yLi93BuN6kiOfkcHnnqOdx4ay84HU6s/X4Rnh/5BCyR0DfNjdfMz+CJM089+iDKlyuLx0eNhWRbniEKciW5FUAVKeVzUsqDuT1ISrlXSmm6hZaH/z2CRq1uxv6DhzD3wyl4543xcDpNmESxWqEui1Wd6R8WqzEb9ITKYgHiipumGt/3tq4QQoRe9Q4Gwq9YB4OAOzm85xJl8vl8mPz+R7iiYTO8PfUDDLv/Lmz76Xtz3RHND8cImp+BCywBIDExAeOeewrrN24O/fpMqhKx8o5HCCHvvaM/Xh37jHG7UYYioRQrFlrwpBauFSGuuLItvdn5MoCMs0ZHAQBoe2tv/Pn3YfyxdS1EqG8Iwvn5Tz/NajeFzev1Yv7Cr/Hia29h774DaNOyGf730mg0alDP6NAKzpFgujUflAODr1mBQADNbuqCA4f+xJ6Nq1CqZEnDYolWIqnCZill45y+FgH3ztRRs0Z1vD/pDXMn3ZzfrR1nYvi7WrqSIiPpBpR2GD1378zD7X16YP/BQ9iweWvoTyroVsKeFCbdFJZTp0/j5f9NQrV6TTHw/odgs9nw9byZWP71Z5GZdAO8YxopDH6dt1qtmD7pDZxNPofHnh5jaCyxKGYS78SEBKNDyB+Tbm2Fk0C7ipomkQ2ZSd4o9Lj1FjidzoLdzgx4gYwzofV7+z3KnQyiEKWlpWPugi/Rpe8glL+yIUa98DJq16yBxQtmYcfPy9G5w02h350xI76GRAYTvEGqe3UtPP34MMya9zmWLlthdDgxJWZaTRo3rC83rVpqdBh5cxZRKrOkHSkBb6oywi6/n31XUeXWbaRyn1N/g5oC6nXH/Vi19mf889sW2O0FSApszswFo7nUBnwZyt9PBlWJk6KTlBK79vyOH1asxg8rVmPV2p+Rnp6BCuXLoU+PLrhrQB/UvbqW0WGqw2IFEssYHQWFIhgEUo8ZHQU8Hg8atmyPtPR07Px5BYoUYf6hlrxaTUw3t04I0QHARABWANOllK9c9HWR+fVbAKQDuFNKuUX3QLVg1okZ0UQI5Q2OLU5ZkBfwXvoYm1OZAR7pE2ZcmZvrGDhmb0Dv7liw6BssW7kGHW+6MfQn+j1A+ikgvsSFC5GCAeX7xvYSuojb7cYf+w/it737sHX7Tmzcsg2btm7H2WRl4e2VV1yOuwb0Qc+undGy2XWwWs0zBUgVrHZHjqwFlsGAoWE4nU5Mn/QGWtzcDaNeeBmTXud8bz2YKvEWQlihbMZzE4DDADYKIb6SUu4+72EdAdTI/LgOwJTMXyMbx0Dpy2oDEkoqCVwwAMiAUgG3x0fXHHVnorJTmuecIRf5jjfdiOLFimH2/C8KlngDyq6WqceVNw8Wq/JrwMudKWOIx+NBSmoqUlLSkJKaimPHT+DIseP498hRHDl2HEeOHseRY8dw+J8j+PPvw9nj0Ww2G+rWvgq9u9+K6xo3RLs2Lc2/8U1hcaJJZLHaDU+8AaDZddfi4SH3YOKU6Wh/Y2vc2rG90SFFPbNlGE0A7JNSHgAAIcRcKJvxnJ94dwXwsVSusOuFEMWEEOWllEf0D1dFJuj5ikmxcJfB7lL+np4U3VtPnE4nenXrjNmffYG0tPTwNiCRQSDAlhItSSmRkpKKtPR0pKWlK7+e93u32wN/IACfzwe/X/nV5/fD7/crn7voa1l/9vn82Y/1er3Kn/3K571eX/bvz39MhtudnWj7fL5cYy5SJBHly5ZB+XJl0ey6xhjUvxeuqnEFataojlo1a8DlMn6dg674GhJZLHYA5thp95Uxo7B67XoMGvIotv30ffS/STWY2RLvigD+Pu/Ph3FpNTunx1QEcEniLYS4H8rumqh8WUVVA1UdL5qkJSGU1hN7nJJ8+926VY4H9O6BaTNmYdHipejfq4cu5ySFlBLHjp/AgUN/4uCff+PQn3/j2IkTOHHylPJx6jROnDyFk6dO55nkhspqtcJut8Nut8FmtcFut8Fut8OR+Tm73Q67LfNzDjvsdjtcTmf2c+w2O+LiXCiSmICiRYqgSGIiihRJUH5NTESZ0qWyk+2I2EVSL7xjGnlM9P1yuVyYP+M9NGp9M/rePRSrFn9RsDU5VCBmS7xzWk5+cXYQymOUT0o5DcA0QFlcWbjQNMbEm/RgtSu7eQaDgC9dScADhU+48tLi+ia4rFIFzJ6/kIm3hgKBAHbt+R3rftmEbTt2Yeee37Bzz+9ITr5wRGNSUlGULlkSpUuVQNXKlXBtw/ooU7oUShQvhsSEBCQkxCMhPvMjIQ4J8fFwuVyw22yw2ZRk2mazXfJnm80W2RNBIpnFZprNsyhEJmsNuqJ6NUyb8Br63fMAnn3xVbz6wrNGhxSRUlPTsH3X7jwfY7bE+zCA8/dwrwTg3zAeE1lYrSC9WSxK/7czUekz9Ls123bdYrGgf8/ueGPSVBw/cRJlSpdS/Ryxau++/fhm6TJ89+NK/LxxM1JSlPGKxZKSUPfqqzCgV3fUqlkDl1etgmpVKqNq5UqIi4szOGpSHV8/Io/ForxhCvqNjiRb357dsGLNOrw28V1c06AeevfoYnRIpnbs+Als3b4T27bvxNbtO7F1+y7sO3AQ+U0LNNU4QSGEDcBeAG0B/ANgI4D+Uspd5z2mE4BhUKaaXAfgbSllk/yObepxgjanMr2ByGjBYGYS7lY1Cd+5+zfUvf5GTHr9JQy7/27VjhuLftv7B2Z++hkWLPoW+w4cBADUqlkDbVo0Q7PrGqNZk8aoVrUyq8+xJK6Y0kZGkSXjDOAzR593FrfbjXZd+2Djll+xbNE8tGwW+bMr1HLs+AmsXLMOK9asw/LVa/HH/gPZX6tWpTIa1L0aDevXQYO6V6NL3ztzHSdoqsQbAIQQtwCYAGWc4IdSynFCiCEAIKWcmjlO8B0AHaCME7xLSrkpv+OaOvF2Jioj7ojMJOBXZp77MlQ5XP3m7RAfF4efl32tyvFiSUZGBj6Z+zk+nDUXv2zaAqvVinZtWuLWjjehU/t2qFrlsvwPQtErscyFYzcpMnhSDR33mptTp0+jefuuOH7iFNb9sAhXXVnD6JAMkZ6ejh9WrMaylWuwYs067NrzOwBlYXfr5k3Ruvn1aNywPhrUvRrFiiVd8Ny85nibLvHWiqkT7/gSsTFdgyJTMKCMIyxkZea1CZPx1PPjsG/rOlS/vKo6sUW5M2fO4t3pM/H2ex/g+ImTqHt1LQzq1wsDevdAubLcLIXAjXMimd+r7FdgQgcP/YWm7TpnF0ti5Xpz5sxZfPPdMiz8egmW/rgCGRluxMfHoUXTJrixVXPc0Ko5GtWvC5st705tJt4weeJdpBwXxpD5edOU6kyY14y/D/+DKnWaYMzTT2D0U4+rHFx0SU9PxxuTpuL1t6cgNTUNHW+6EU89+iBaNW/KFhK6kN0FxBU3OgoKh5RAylGjo8jVpi2/onWnHqhcqSKWLZqHihXKGx2SJo4dP4EvvlqMhd8swYo16+D3+1GxQnl079wB3Tt3RIvrm8DhKNgADCbeMHHibbUDCVxsRhEi4AMyzoa9IOiGzj3x75Gj+G3TGiaQOQgGg5iz4EuMHDMOh/85gtu6dMLopx5DvTq1jQ6NzMpZRGlXpMiUesJUCywvtmbdL+jUeyBKlSiBZYvm4fJqVYwOSRVutxtfL/kBM+d8hqXLViAQCKBG9cvR49aO6NHlFjRuWB8WiyXs4+eVeId/VFIHV6NTJLHagfiSymr8MAzo1R179x3A5q3bVQ4s8u0/cAhtOt2G2+8bhrKlS2P1koVY8Mn7TLopbxxFG9lMngO0bHYdln/1GZLPpaBFh27Y/dteo0MKm5QSv2zaggcefxrlr2yI3ncOxrYduzDikQew4+fl+H3zGrwy9hk0uaZhoZLu/DDxNpqVvd0UYSwWZV1CGIu5enbtDIfDgdmffaFBYJFJSokp02eiXvO22L5rDz5453/YsGIxpwlQ/jiKNvJFwPevcaP6WLX4c0gp0bx9V3y95HujQyqQf/49glfenITaTVqjadvOmPHpPNzS/kZ8v3AO/ty5AeOffxp1al+l211YJt5GY7WCIpHFqvSVFvBCVaxYEjq1b4s5C76E32/e26t6OX7iJDr06I8HnngaLZo2wc6fl+Pugf00rbZQFLHYuT4o0kVIDlCn9lVY9/1XqFblMnTpeyeeHjPe1Nfw9PR0fPrZF7i5ez9cVrsxnh77MkqVLIHpk97A0b2/Yvb0ybjpxtawWvWfBsSru5EsNqV6SBSJrPawFnXd3uc2HDt+Aj+sWK1BUJFj89btaNymA1av+wXv/u9lLP3iU1SqWMHosCiSREC1lPJhjZw3T9WqVsa6H77CfYMG4JW33sFN3fri0J9/Gx1WNiklfvr5F9z30HCUr9kQA+4dht/37cdzIx7Fvq3rsGbpl7jnjv4oWtTY8c1m27kytvCiSZHO5gQcCcrEkxB17tAOpUqWwIefzEXHm27UMDjz+njOZ7j/kREoW6YUflr6Ja5pWM/okCgSRUi1lPJhsQMBr9FRhMTlcmHa26+jxfVN8MATT6NWk9Z4+vFhePLhoYbsiiulxPaduzFnwZeY+8Ui/PnXYSQkxKNn1864s39vtGre1HR3EM0VTazhRZOigbNIgRZbOhwO3N7nNixa/B1OnjLnDFutSCnx9JjxGDTkEVzf5BpsWrmUSTeFj68h0SECi3B39OuF3zauRpeON+H58W/g6utuwJwFC+Hz+XQ5/779B/Hia2/h6uvaoEGLm/DGpKmoXfNKfPze2zi691fMmDIBbVo2M13SDTDxNhYvmhQNhABcSfk/7jx3394XPp8Pn362UKOgzMfv9+O+h4bjlbfewf133o4fvpyL0qVKGh0WRSq2KkaPCEy8AaBSxQqYN+M9/PjVfMTHx6H/PQ/i8vpN8epb7+D06TOqnisQCOCXTVvw/PjXcU2rm1GjUXOMHvc6SpcqiSlvvoKjf2zD4gWzMLBvTyQmJqh6brVxjrdRhAUoUtboKIjU404GvOkhP7xx6w7wB/zY9tMyDYMyB7fbjX73PIAvv1mKZ598FC888yTnmFPh2OOAuGJGR0FqCAaA1ONGR1EowWAQi7//EW9Nfh/LV/8Eu92O1s2botPN7VCkWn0s2u/FsXMelEtyYWjr6uhQp1yex0tPT8e2Hbuwdv1GrNuwCWvW/YJTp8/AYrHg+ibXoFunDuh7W1fTrovhBjowYeLN3cYo2kgJpJ1QXkRC8O77M/Dg8FHYsvo7NKxfV+PgjJORkYFb+wzCj6t+wsRXX8DDQ+41OiSKBq4kwBFvdBSkltTjIV87zW77zt34ZO4CfPv9j9jz+x8AAKsjDq5SFRFXsgLikkqiXd3L0OjyMrDb7DiXkoKzyedw+sxZ7D94CL/v24+//v4n+3hXXF4NzZteiw5t26D9ja1RooT5cycm3jBh4u0qqixKI4om3nSl8h2CM2fOonzNhrhvUH9Men2cxoEZw+PxoGu/u/D98lX46N23MKh/b6NDomiRUBqwcj5C1Eg/Dfg9RkehunYvfoF92zch48TfyDj1L9yn/oUv7RyAC3NPq9WKYklFUa1KZdSsUR01a1RH3dpXodl116JM6cjb3TuvxJv/a43C/m6KRo54ZcJJCFsgFy9eDN07d8Ds+Qvx+ovPweVy6RCgfrxeL3oNuh/f/bgS0ye9waSb1CMsTLqjjdURlYl3ijUJZRq2veBzUkpIvxdLHrwOXq8XSUWLIiEhPmba77gywwjcbYyimTMx5IfefXtfnDl7Fl9+Y6K7USrw+/0YcO8wfL3kB0x+YzzuuaO/0SFRNOHrR/SJ0mJcuaRLCypCCFQslYTSpUqiYoXySExMiJmkG2DibYwo/Q9GBEBZ9BXieMEbW7dA1cqX4b0ZszQOSj9SSjz4xCgsWPQN3hw/Bg/cd6fRIVG04WtI9InSN1NDW1eHy35hqumyWzC0dXWDIjIeE28j8KJJ0S7EqrfVasWQuwdi5Zp12Ln7N42D0se4NyZi2oxZGPnYMDz24P1Gh0PRiK8h0SdK74R3qFMOozrWQvkkFwSA8kkujOpYK9+pJtGMiyuNEF8SsPHCSVEu9URIvd4nT51CpVqNcfftffHumy/rEJh2Zsyeh7seeAy397kNH7/3dkzdPiWdCAEklo2YbcapAAo4kpXMK6/Flax46y1K39USXSLEqnepkiXR97au+GTeApw7l6JxUNr5/seVuO/hJ9GuTUt88M7/mHSTNqwOJt3RysLcIBYw8dabxc6LJsUGexxgsYb00AfvvROpqWn4eO5nGgeljb379qP3XUNQq2YNfP7JdDgcvKNFGmHhJnqxhSgmMPHWG1tMKJbY40J62LXXNMC1jRpg8vszEGntb8nJ59Cl752w22z4as4MFC1axOiQKJoxOYteVpsyKpKiGr/DeuNFk2KJPSHkOzwP3ncnftu7D8tX/aRxUOoJBALod88D2H/wTyz4+H1UrXKZ0SFRNBOCryHRjnc0oh4Tb73xokmxxGIBbM6QHtqnRxeULFEc77z/kcZBqWfU2Jex5IflmPTaS2jd4nqjw6Fox1bF6MccIeox8daTlRdNikH2hJAe5nK5MPiugVj07Xf4/Y99GgdVeLPnf4HXJr6LIXffgSH33GF0OBQLWA2NfiEWKihyMfHWEy+aFItsjnx/9pfuPIquk9diibcmLDY7hj33uk7BhWfTll9x70PD0ap5U0x89QWjw6FYwWpo9GOBLuox8daTle9kKUbZ43P90tKdRzF+yR4cSXbDFl8UperdgGVLv8XMH7boGGDojhw9hm4D7kbZMqWw4OP3OcGE9MPEOzbw+xzVmHjrif+ZKFbZ43JdrT9l1X64fcHsP5drcguEEBjz+tt6RRcyj8eD2wbeizNnz2LRpx+hdKmSRodEscJqV9ZMUPTj3fGoxv/FerHYeNGk2CUEYHfl+KWjye4L/uwoWgKl6rbEn5tW4MjRY3pEFxIpJYY89hR+3rAZM6dMRP26VxsdEsUSJmOxg0W6qMZMUC+c302xzpbzTO9ySZcm5OWu6wwZ9ON/k6ZqHVXIJk6Zjhmz52P0U4+hZ7fORodDsYatirGDiXdUY+KtF/5Holhncyh3fi4ytHV1uOwXXoqKlSmPG9t3xNSPPsGJk6f0ijBXPyxfhSeeGYvunTvi+ZFPGB0OxSK+hsQOIXiHI4ox8dYLqxVEObabdKhTDqM61kL5JBcEgPJJLozqWAuTX3oKbrcHz49/Q/84z7Nv/0H0uWsoal91JT5+721Y2DJGerM62KoYa/hGK2pdWn4i9bG/m0hhjwc8qZd8ukOdcuhQp9xFny2HB++7E+9M+wiD77rdkJ7qM2fOonOfOyCEwKJPP0JiYmgzyYlUxepn7LE6AKQZHQVpgNmgHtjfTaSwWAu0QcSYkU+geLEkPPzUc5BSahjYpbxeL24beB8OHPoLC2d/gMurVdH1/ETZuKlK7GHFO2ox8dYD/wMR/ceW83STnBQvXgzjnhuJ1WvX47OFX2sY1IWklHjg8aexYs1aTJ/0Blo1b6rbuYkuIARfQ2KRxZLjmhiKfEy89cD+bqL/2OMKtDPbvYP6o0HdqzH8uReQnp6uYWD/eW3CZHzwyRw8++SjuKNfL13OSZQjq4M7GcYq3i2PSky8tcb+bqILCVGgW+dWqxVvv/YS/j78L8a+8qaGgSk+/GQORo4Zjz49umDsqOGan48oT+zvjl0s2kUlZoRa4ztWokvlMtM7Ny2bXYf7Bg3AaxPfxbffLdMoKGDe54tw70PD0f7G1pg5dSInmJDxmHzFLrYYRSW+qmiN/3GILmVz5rqFfG4mvvoC6tetjYH3P4w//zqsekhfL/ket9//EFpc3wQLZ38Ap5MJDxlMCBZvYpnFwhwiCjHx1hqrFUSXKmC7CQDExcVhwcz3EQgG0GvQ/fB4PKqF8/WS79Fr0GA0qHs1vpn3MeLj41U7NlHYmHQRW42iDhNvLbG/myh3BZhukuWK6tUw490J2LhlGx55arQqIwYnv/8RuvW/G3VrX4WlX8xG0aJFCn1MIlUw8SaOkow6zAq1xFuERLkLo90EALrf2hEjHnkA7330Ce4d9gR8Pl9Ypw8EAnh81BgMG/4MOndoh5Xffo6SJUqEdSwiTTDpIk61iTocEqkltpkQ5S6r3cSXUeCnvjL2GTidDrz42gT8e/QYPps5rUC7Sh489BeGPPYUvl++Cg8NvhtvvTwWVqu1wHEQacZiZZsB/TfH3a9eax0ZixVvLbFaQZQ3e8Gmm2QRQuCFZ0Zg2sTX8cOK1Wh9Sw9s274z3+d5vV68/L9JqH1dG6zbsAnv/u9lvP3aS0y6yXzYZkJZ+LMQVVjx1gpvDxHlL6vdRAbDevp9dw5AxQrl0O+eB9CwZXvcdEMrPPnwULS7oRVE5v8/KSU2b92Or5Z8hzkLFmHfgYO4rUsnTHhlLCpVrKDm34ZIPWGsgaAoZXMCnhSjoyCVMPHWCvu7iUITZrtJllvat8WfOzZg6ocfY+LUD9C+ez+4XC4kFS2CpKJFkJKahiNHj8FisaB502sx4ZWx6HRzOxX/AkQqC2PqD0Uxq71QBQoyF6HGVIBI0Lhhfblp1VL9Thhfksk3USj8HiD9tCqH8ng8mPv5Iuzc/RvOpaQi+dw5WK1W3Ny2DW5pfyNKlSypynmINGVzAvFc6EvnyTgD+NxGRxEeV5LyZjLgBQI+5SPKiaQKm6WUjXP6GiveWuCmB0ShK2S7yfmcTicG9e+tQlBEBmJPL13M6ozMxNtVFHBk7ouQtaYnkt9EqICLK7XAiyZRwfCNKtF/2N9NF4vEnwlHgvJxMWdSWKNko0Xs/s21xN48ooKxhTfdhCjqWKyAlTej6SKRtn283aVUu3NisQDO2N2ojIm3Fji/m6hgbE5OASICIiu5In1FSlFPWABXsbwf44iP2Z91Jt5qY7WCqOCyNokginWR2FJA+oiUnw1nYmiFlKxFlzHGNIm3EKKEEOIHIcQfmb8Wz+Vxh4QQO4QQ24QQm/SOM19MHojCE+ZmOkRRg2MEKS9WG2AxeWHPYgXs8aE91moL/bFRxDSJN4CRAH6UUtYA8GPmn3Nzg5SyQW6jWgzFiyZReNiiRbGOG69Rfsy+EN0RYrU7+/E5LL6McmZKvLsCmJn5+5kAuhkXSiEweSAKj8XCN64U2yKllYCMY+afEYvtv9GBIT/HGnPXfTMl3mWllEcAIPPXMrk8TgL4XgixWQhxf14HFELcL4TYJITYdOLUKZXDzYHVoSQPRBSeGLsAE13AzEkVmUPWvgdm5EwM73kxVvXWtVlICLEMQLkcvvRMAQ7TXEr5rxCiDIAfhBC/SSlX5/RAKeU0ANMAZefKAgdcUEwaiArH5gJwzugoiPTHwg2FyuYw3wY0Flv463RsTuX5Qb+6MZmUrom3lLJdbl8TQhwTQpSXUh4RQpQHcDyXY/yb+etxIcRCAE0A5Jh4647VCqLCsViVBCTgNToSIn3Z+fpBIbK5zJd4F7Zq7YgH3LFRdDHT2+uvAAzK/P0gAIsufoAQIkEIUSTr9wDaA9ipW4R54RhBInWYffEQkRZYuKFQ2VzmWoRrsRZ+KpUtzlx/Jw2ZKfF+BcBNQog/ANyU+WcIISoIIRZnPqYsgJ+EEL8C2ADgWynlUkOivRjbTIjUwQSEYo3VriQvRKEw29hJe3zhk2aLJWZGypqmRCulPAWgbQ6f/xfALZm/PwCgvs6hhYbJApE6spKQYMDoSIj0wdcPKih7vDnaTYRQbxa3PR7wpqtzLBMzU8U7cnHXPSJ1mamaQ6Q1Jt5UUGaZbmKPV29RsNWufEQ5E3zXogA3PSBSFxMRihUWG9cHUXiMbs0QQv1RgEb/nXTAxFsNTBKI1GV1mKOaQ6Q1TjOhcBmdpNpc6q9NiIFFlnxlUwNvixOpSwhON6HYYIv+Ch9pxOjWDEeYG+bkJQZ2MGbiXVhWB1ejE2mBd5Io2lkdbDOhwjHqOml3afezG+VvRpl4FxZvExJpwxrdVQ8ivn5QoRnVbqJFtTuLXYMWFhNh4l1YrMoRaSMGbjlSDBMi6it7pAOLVf/rpM2pfYtLFOdWTLwLg20mRNpi4k3RyupQbwwbxTa1J4vkx1lE+3OoNRvchPi/vjB4m5BIW1Fc9aAYF8WJBelMjwq03uey2qJ2fxQm3oXBpIBIWxarMueYKJoItlGRyvSqemvZ230xo8claoSJd7jYZkKkDyYoFG1szqifVUw6s8dpX6SwOfUd82qPzpneTLzDxTYTIn3wzhJFG7aZkBYcGv5cCQE4i2p3/NzOGYWFFybe4eJqdCJ92LiLJUURq52bQ5E27PHaXSvt8cbMnI/CN6l8NQuHzcnV6ER6isKqB8WoKEwkyCSE0KbXW1j0mWSSE5sz6tp6mT2GI0ob/olMi+0mFA2Eha8fpC1Hgvq93s5EY3uto+z/DBPvghIWJgFEeuNiNIoGdhd/jklbQgCuJPWOZ7XrPyf8YlF2l4iJd0HxwkmkPyGidqYrxRC7wQkMxQabQ50qsRCAq1jhj1NYRuzOqSEm3gUVZe+8iCJGFF14KQbZnMYsTqPY5Cxa+IWWrmLm+ZmNotyLiXdBWO367Q5FRBdiixdFsihKHCgCWCyAqxDj/xzx5hqbbHdFzSJLJt4FEWUN/kQRxWLlG1+KTBabuZIYig32OGVhZEFZHfrP7A5FlORgTLxDJQRndxMZje0mFImMXpxGsctZpGCjAG1OIK64OdeyRcldIybeoeLsbiLjsd2EIo3FGjWVOopQzsTQ2k4cCUB8CfPmOlGyyNIkXfMRgKvRiYxntSsX32DA6EiIQuNIMGf1kGJL1nxvXzrg9wBS/vc1qx1wJEZGO5QjQYk/gjHxDoXVwS1+iczC5gS86UZHQZQ/izVqbo9TFLA5lQ8pAb9bmXpidUTWG0ObU3kDEfQbHUnYTHo/wWTYn0dkHmw3oUjBajeZkRBK+1OkbkwW4a1bTLzzw9XoROYSaRUaik2sdhNpwx4f0a8BTLzzw2o3kbkIERULbCjKsdpNpA2LJaKr3ky88yIi+5tLFLXYbkJmZrGxaEOkJUcY88lNgol3XlixIDIna4T2JlJsKMyOgUSUP4s1YtuAmXjnRljYn0dkVpbM1fhEZpM1OYKItOUowMZAJsLEOzeOBPMOkSciJjdkPkIAriSjoyCKDVZbRL4OMLPMicXK/jwis7Nx/QWZjCNBef0gIn1EYK83E++cOIuwf5TI7NhuQmZisUVkEkAU0WyOiHsdYOJ9MauDk0yIIkWELq6hKCMEEFeMBRsiIzgjq9ebiffFIuwbSBTTOFaQzMCRAFjtRkdBFJtsjojq9WbifT67S/kGElFksFgj7jYjRRmrgwUbIqM5I2eEJxPvLBYr4ORqdKKIE0GVDooyWS0mRGQsqy1i2oSZeGdxJXF8IFEkipCLLUUhVzFOMSEyiwgZjMFMEwAc8ayaEUUqi5X9taQ/ZxEu7iUykwgZBc3E22KLqN4gIsoBF1mSnhzxgJOjA4lMx5Go5HUmFtuJN0dAEUUHtpuQXmxO7k5JZFYRsHtsbCfeccV5i5ooGlisbBcj7dmcyusGEZmXzaHclTKp2E28XUl8oSaKJmw3IS1lJd28Q0pkfs6ipm05ic3E25lo6ndDRBQGm4tJEWnD7gLiS/DniyhSmLjlJPYSb0cCNzsgikYWC+9ikfqciWwvIYpENnNubhVbiberqPJBRNHJxkWWpBJhUarcJnzhJqIQORNNt/jenA0wWrDYImK+IxEVgs2pJEwyaHQkFMmsDmXiFTfHIYp8riQgGAACXqMjARBLFW/25hFFPyG4qQmFT1iUF+mEkky6iaKFEEq7mEn+T8dO4k1EsYHtJlRQQigL7hNKc+E9UTSyWID4kqaYdGJ8BEREarI5lItr0G90JGR2QgD2eKUN0STVMCLSiMWqJN8Zp4GAz7gwDDvzRYQQvYQQu4QQQSFE4zwe10EI8bsQYp8QYqSeMRJRhGDVkvJitSuLJhPKKAvumXQTxYasyreBE7BMk3gD2AmgB4DVuT1ACGEFMBlARwC1AfQTQtTWJzwiihj2eK7roAtlLbBPKKV8OBOVF2Eiii1C/DexyIDXCdO0mkgp9wCAyPsfoQmAfVLKA5mPnQugK4DdmgdIRJFDCGWElDfd6EjIKBabUtm2OpTqFqvaRHQ+Z6Ky8Zr7rK6tJ6ZJvENUEcDf5/35MIDrcnuwEOJ+APcDQOXKlbWNjIjMxZ7AxDvaCYuSUAuLkmhnfVjtvONBRPmz2pQ7YN50wJumy9ogXRNvIcQyAOVy+NIzUspFoRwih8/J3B4spZwGYBoANG7cONfHEVEUstqUaqdJZrfGrAsSYHHR587/s/jv8+f/WVgyf2+58MNiZXJNROpwxCsfPjfgS1deN6Q2aaOuibeUsl0hD3EYwGXn/bkSgH8LeUwiilaOeCCDifcFxMVJ7Hl/viTZPS9BzkqGs46R9fnzj5vT74mIIoXdpXxICfjdgN+jVMGDftUS8UhrNdkIoIYQohqAfwD0BdDf2JCIyLTscYAnRdm1LJYIAVjsSlXYYvvvV2HlgkIiovxkrRM6f7v5gB+QAWVn5KyPS5+Yb+HBNIm3EKI7gEkASgP4VgixTUp5sxCiAoDpUspbpJR+IcQwAN8BsAL4UEq5y8Cwicjs7PFK8h3NrPb/FhJaHVxISESkNqsNaqTNpkm8pZQLASzM4fP/ArjlvD8vBrBYx9CIKJI5EpRFMzlWJyKUsCiTOmxOwOpkFZuIKEKYJvEmItKEEEryHelV76xk2x5n6OYPREQUPibeRBT9HAnKSvVI7PW22pV2GXscFy0SEUU4Jt5EFP2yqt7uc0ZHEjqbE3AkAjaH0ZEQEZFKmHgTUWywx2dukGDyqrfdpSTcVrvRkRARkcqYeBNRbDB71dvmBJxFmHATEUUxJt5EFDscCcrOZGbazdJiA1xFuWCSiCgGcAYVEcUWV5I5FikKoSTciaWZdBMRxQhWvIkotlhtSg+1keMF7S7AmcT520REMYaJNxHFHmci4HcDAZ++52VbCRFRTGO5hYhik54tJ1kLOxNKMekmIophTLyJKDZZ7YCrmA7ncQDxpZRKtxl6y4mIyDBsNSGi2GV3ATIJcCerf2xhUcYDOuLVPzYREUUkJt5EFNsc8QCkuvO9HfGAowgXTxIR0QWYeBMRORKUXz0pgJThH8fmBJxFlckpREREF+GrAxERoCTfVqfSdlLQDXa4zTsREYWAiTcRURarDUgoCXjTAV9G3gm41a4snLTHs8JNREQh4asFEdHFHPHKRzAIBDxA0A9AKFNJhEWpjLN/m4iICoiJNxFRbiwWwBJndBRERBQlWLIhIiIiItIBE28iIiIiIh0w8SYiIiIi0gETbyIiIiIiHTDxJiIiIiLSARNvIiIiIiIdMPEmIiIiItIBE28iIiIiIh0w8SYiIiIi0gETbyIiIiIiHTDxJiIiIiLSARNvIiIiIiIdMPEmIiIiItIBE28iIiIiIh0w8SYiIiIi0gETbyIiIiIiHQgppdEx6EIIcQLAn0bHcZ5SAE4aHQSpjt/X6MTva3Ti9zU68fsanSLp+1pFSlk6py/ETOJtNkKITVLKxkbHQeri9zU68fsanfh9jU78vkanaPm+stWEiIiIiEgHTLyJiIiIiHTAxNs404wOgDTB72t04vc1OvH7Gp34fY1OUfF9ZY83EREREZEOWPEmIiIiItIBE28iIiIiIh0w8daZEKKDEOJ3IcQ+IcRIo+OhwhNCXCaEWCGE2COE2CWEeMTomEg9QgirEGKrEOIbo2Mh9QghigkhFgghfsv8v3u90TFR4QghHsu8Bu8UQswRQriMjonCI4T4UAhxXAix87zPlRBC/CCE+CPz1+JGxhguJt46EkJYAUwG0BFAbQD9hBC1jY2KVOAH8ISUshaApgAe5Pc1qjwCYI/RQZDqJgJYKqW8CkB98Hsc0YQQFQE8DKCxlLIOACuAvsZGRYUwA0CHiz43EsCPUsoaAH7M/HPEYeKtryYA9kkpD0gpvQDmAuhqcExUSFLKI1LKLZm/T4HyAl7R2KhIDUKISgA6AZhudCykHiFEUQCtAHwAAFJKr5TyrKFBkRpsAOKEEDYA8QD+NTgeCpOUcjWA0xd9uiuAmZm/nwmgm54xqYWJt74qAvj7vD8fBhO0qCKEqAqgIYBfDA6F1DEBwAgAQYPjIHVdDuAEgI8y24imCyESjA6Kwiel/AfAGwD+AnAEQLKU8ntjoyKVlZVSHgGUgheAMgbHExYm3voSOXyO8xyjhBAiEcDnAB6VUp4zOh4qHCFEZwDHpZSbjY6FVGcD0AjAFCllQwBpiNDb1qTI7PftCqAagAoAEoQQtxsbFdGlmHjr6zCAy877cyXwVlhUEELYoSTds6WUXxgdD6miOYAuQohDUNrCbhRCzDI2JFLJYQCHpZRZd6YWQEnEKXK1A3BQSnlCSukD8AWAZgbHROo6JoQoDwCZvx43OJ6wMPHW10YANYQQ1YQQDigLP74yOCYqJCGEgNIrukdK+abR8ZA6pJRPSykrSSmrQvm/ulxKyQpaFJBSHgXwtxCiZuan2gLYbWBIVHh/AWgqhIjPvCa3BRfMRpuvAAzK/P0gAIsMjCVsNqMDiCVSSr8QYhiA76CsuP5QSrnL4LCo8JoDGAhghxBiW+bnRkkpFxsXEhHl4yEAszOLIAcA3GVwPFQIUspfhBALAGyBMmlqK6Jki/FYJISYA6ANgFJCiMMAngfwCoD5Qoh7oLzR6mVchOHjlvFERERERDpgqwkRERERkQ6YeBMRERER6YCJNxERERGRDph4ExERERHpgIk3EREREZEOmHgTEREREemAiTcRERERkQ6YeBMRERER6YCJNxERAQCEEMWEEIeFEB9f9PmvhBB7hRDxRsVGRBQNmHgTEREAQEp5FsA9AAYKIboBgBDiLgCdANwppUw3LjoiosjHLeOJiOgCQoj3AHQD0AHACgDvSSmfMjQoIqIowMSbiIguIIRIBLAdQAUA+wBcI6X0GBsVEVHkY6sJERFdQEqZCuAbAE4AHzDpJiJSByveRER0ASFEYwA/A9gBoAqAq6WUR42Niogo8jHxJiKibEIIF4AtAA4A6A3gVwB7pJRdDA2MiCgKsNWEiIjO9xKAcgDuy5xiMghAJyHEnYZGRUQUBVjxJiIiAIAQojmA1QAGSik/Pe/zrwO4D0AdKeVho+IjIop0TLyJiIiIiHTAVhMiIiIiIh0w8SYiIiIi0gETbyIiIiIiHTDxJiIiIiLSARNvIiIiIiIdMPEmIiIiItIBE28iIiIiIh0w8SYiIiIi0sH/AcN7cZ5jLY4FAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, axs = plot.subplots()\n", "\n", "\n", "axs.scatter(x_t, y_t, label='Training Data')\n", "axs.fill_between(x_p, mu_p - sig_p, mu_p + sig_p, alpha=0.1, label=r'$1\\sigma$ Confidence')\n", "axs.plot(x_p, mu_p, label='Predicted Function', color='k')\n", "axs.set_xlim(min(x_p),max(x_p))\n", "axs.legend()\n", "\n", "axs.set_xlabel('x', size=15)\n", "axs.set_ylabel('y=f(x)', size=15)\n", "\n", "fig.set_size_inches(12,7)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.7.6 64-bit ('base': conda)", "language": "python", "name": "python37664bitbasecondab95d1cd6133a44a5a90e16cd715dd17b" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }