{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A simple demonstration of coregionalisation\n", "\n", "This notebook shows how to construct a multi-output GP model using GPflow. We will consider a regression problem for functions $f: \\mathbb{R}^D \\rightarrow \\mathbb{R}^P$. We assume that the dataset is of the form $(X_1, f_1), \\dots, (X_P, f_P)$, that is, we do not necessarily observe all the outputs for a particular input location (for cases where there are fully observed outputs for each input, see [Multi-output Gaussian processes in GPflow](./multioutput.ipynb) for a more efficient implementation).\n", "\n", "For this problem, we model $f$ as a *coregionalised* Gaussian process, which assumes a kernel of the form:\n", "\n", "$$\\textrm{cov}(f_i(x), f_j(y)) = k_1(x, y) \\times B[i, j].$$\n", "\n", "The covariance of the $i$th function at $x$ and the $j$th function at $y$ is a kernel applied at $x$ and $y$, times the $i, j$th entry of a positive definite matrix $B$. This is known as the **intrinsic model of coregionalisation (ICM)** (Bonilla and Williams, 2008).\n", "\n", "To make sure that B is positive-definite, we parameterise it as:\n", "\n", "$$B = W W^\\top + \\textrm{diag}(\\kappa).$$\n", "\n", "To build such a model in GPflow, we need to perform the following two steps:\n", "\n", " * Create the kernel function defined previously, using the `Coregion` kernel class.\n", " * Augment the training data X with an extra column that contains an integer index to indicate which output an observation is associated with. This is essential to make the data work with the `Coregion` kernel." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import gpflow\n", "import numpy as np\n", "import matplotlib\n", "%matplotlib inline\n", "\n", "matplotlib.rcParams['figure.figsize'] = (12, 6)\n", "plt = matplotlib.pyplot\n", "from gpflow.test_util import notebook_niter\n", "np.random.seed(123)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data preparation\n", "We start by generating some training data to fit the model with. For this example, we choose the following two correlated functions for our outputs:\n", "$$\n", "\\begin{align}\n", "y_1 &= \\sin(6x) + \\epsilon_1, \\qquad \\epsilon_1 \\sim \\mathcal{N}(0, 0.009) \\\\\n", "y_2 &= \\sin(6x + 0.7) + \\epsilon_2, \\qquad \\epsilon_2 \\sim \\mathcal{N}(0, 0.01) \\\\\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAD8CAYAAABevCxMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X181PWd9/vXN0Mm5AYtJFE5VIRUscXSrZijoRZsL+2i0Utk22o5rGCvVmSlx6qn6+l19py6uo+9Hl22K61tXEXbKhyWit2HNwcR1nusTdRAu1JRbuSm6oWQhKhMJmSSyff88ctvMjOZSSbJ3P1m3s/HYx7MzY+Znz9DPvP9fj/fz8dYaxERERHvKMn1CYiIiMjoKHiLiIh4jIK3iIiIxyh4i4iIeIyCt4iIiMcoeIuIiHiMgreIiIjHKHiLiIh4jIK3iIiIx0zI9QkkU1NTY2fMmJHr0xAREcmaHTt2tFtra0c6Lm+D94wZM2htbc31aYiIiGSNMeZwKsdp2lxERMRjFLxFREQ8RsFbRETEYxS8RUREPEbBO9NefxACbYOPA23OcyIiImOUt9nmBeH1B2HLD+CNh2D5Zue5R66Ctnec+xfemLtzExERz1LwzqTZ1ziBu+0duK/BeS7YDrWfdV4TEREZA02bZ1JVrTPirqhxgnaw3bm/fLPzWrppil5EpCho5F0oNEUvIlI0NPLOpECbE0DdEbc7An/kqtgRcjrMvsaZjnen6O9rcO5ril5EpOAoeGfS7icGA+jNLc7NDbC7n0jvZ2V7il5ERHJG0+aZ5E5Vz75mMIAu3+wEbk1ji4jIGGnknWkX3hg78q2qzUzgzuYUvYiI5JSCd6HI5hS9iIjklKbNC4Wm6EVEioaCdyGJD9KZmqIXEZGc0rS5iIiIxyh4i4iIeIyCt9eoBKqISNFT8PYStwSqu/3L3R625QfJA7iCvYhIwVHw9pJUSqBGB2s32N9/cerBXkRE8p6Ct5eMVAI1fmR+1nzw+SFwFH72F6p3LiJSIBS8C0n8yPyRKyEcAuOD3q7s1TvXVL2ISEYpeOeLVALeSCVQE43My6fAxFMH36O3G7o6hv+c8ZzzphtGvy4vIiKjkpbgbYz5lTHmmDHmT0leN8aYe40x+40xbxpj5qbjcwtGqoloYymBevJj6D7uBPTSSmcEvnY+HH1nfIE12TnvfhyqTh92XX5d8yHaAz2Rt2oP9LCu+dDoPj/J+6zasCMt7y0iks/SVWHtYeAXwLokr18BnDNwuwj414E/BZyg9sZDgwEPnFFz/Nr0SCVQ40fmfd0Q6nLWvZc/7Ry/dr4zlf7L/wITyhN/zjjO+aPKOvq//jBTHvsr572Bbv9kygem6tc1H+JHT77F+ubDbFzh/L0la1vYdywAwLJ5M1L6+ETvc8VPt9MWCPH6oU6e+f78Mb+3iEi+M9ba9LyRMTOAzdbazyd47QHgJWvtxoHHe4CvWGuPJHu/+vp629rampZz84RAmxMEBwIeFTXOyHo0a9PuaLj2s05Qh4FM86Nw6Z0w/3ZnxP3QV6E3OPbPSXLO3f7JfPmT/0FdTSWPhm+npNt5vt2ewgtffYJrv3IB7YGeSECtrvQD0NEV4pzTqti4ooGaqrJhP3Jd8yEa50wFBgNzpd/HxFIfHV0h/L4SQuF+qiv9nOwN0xUKR94bYMuuIwriIpK3jDE7rLX1Ix2XrTXvacB7UY/fH3iucIy0Zp2NJK4Lb4TGnwwmpFXVwkUrndfefHTw8/sGp5XpG2ENfBTnXebz8cUpYf7xk/9OSXc7xzmFdnsKNeYTLnv9u3QcfZ+aqjI2rmhgckUpHV0hOrpCVFf6Uw7cP3ryLZasbQGgaelcfMbQFQpH3mfzLV+mutJPR1eIrlAYnzE0LXVWaZasbeFHT76laXQR8by8akxijFkBrACYPn16js8mzusPxk5XB9oGp6vdEe8bDw2OeB+5yplSdg33+uxrYqe7YTARbbSZ4fGNSM6/3gncbe9A04XOGrgNOxnoE8qcafW182HFK1BZnfp595xw3jfqnEuC7awt/3t8Je3s7Z/GzaFbuHrif/LXFS1MCR7g57++lyXf+weOd4X4pLsv9f+mAY1zprK++TD7jgVYuGY7/dYSjpo5OtkbpjMYivk7YWu57oFmSoyJjPDdkbuIiFdp2jwViaaj3SDW+JPB4Nv2TmzwTXR8otd3PzH8+4+3M1j8lLzxwcrfO/fdNfDSSiiNWgP/wnUwqxF+u9w5j/IpEO5xgr37+vN3DZ7zH9Y7j4EXqecfT36T+/z3MqvkA16e/je8eayff/noEiZXlPJJdx9ha/EZwynlE+gM9qY8bd4e6GHhmu10dA0G6ckVpYT6+ukKhTGAHXju4+5e+qN+vKsr/Wy7bQGg6XMRyU/5Nm3+FLBsIOu8Afh4uMCdd0aqbDZS8ZSRXgdnTTr68ReuS0/gTmTiqc4o+/TPOiNuNwvdPS83MP92OXzjESdwdx93Are/0jnP+bfHTtGffz2d5WcBMMfuZdPEf2RWyQfs7Z/G7Xu/wMnz/xvVlX46g72RwP3MrfN59vZLOOe0KvYdC7Bl1+h/JHzG8OhN83h81cX4jMECxjiv9cd9L+23luNdoaTT56s27GDv0RORx3uPnmDVhh2jPicRkUxL11axjUAzcK4x5n1jzHeMMSuNMQMLrmwBDgD7gQeBm9PxuVmTSvAdK3dU/+ajzmM3Y3xgFDtuifaGdx8f3OJVWe2MuKPNahz8svLwFdDdOfiaL2pkfOGNg//9VbU8cf5DkTXuKXzCcU5hSej/poNTiXdK+QSmVPoja+B3LzpvxJGwm+zmrm9X+n2ErWXVhp1MqfTzzK3z8fsM1kJnsBffQBQfiOV0Bnu54qevsO9YYMj0+aoNO3h614dcde/v2Hv0BHuPnuCqe3/H07s+VAAXkbyTljVva+2SEV63wKp0fFZeig+QELtmDclf/8YjsaN697V0lTCN3hsePyX/h/WRdetu/2TKfD5Kgu3w2+VsPOtuFnXeSEX3cQAshp7SU5noBv4EX1y+/aWZBN/wQa/zuN/CqeWlfKvhM/zHW0cjQRecDPMla1siU+WpTGFv2XUkEnjjt5lt2XWExjlTmTSxNDKlHraWmTUVNM6Zyv/b8mc+7nZG/RV+H4vnTouZov/+ZbN4dvcxQuF+Fq7ZPvDfDH5fCd+/bNZYr76ISEaowloqRqpsNlLxlOFeP/xK5kb1kDgDfflm57mySdD2Dh9V1vHlT/4H1/nuoa/6XGh7h8PN/05/72BW+sdUcdWJ/85HlXWJi8IMXKOK3s6YLPMH7d9TYz5h37EAtZOcUfa22xZEpsrvfDJhXZ+Els2bwd2LzosE/OhRe+OcqTGjcvdLwicn+2h68V0CJwcT5Hp6+1m9dU/MtPms0yex+ZYvR9bMLc6IffMtX2bW6ZNGfdlFRDIpr7LN89Zwo1c34xySF09xDVdcJZPi182ramOeC591BVM2vMsbxwJcXnEH/5Xn+Svfy1SZk3QxkRATmMwJHqr4Bf3f2ABtrw99z4FrdKjkTL4e/L+YXOHngf47+Qzvs2Xnb7nsc9/kubePsWrDTjauaKBp6VyWPtjC07s+5KLmQyknj8Uf547a1zUfSjoq95UYwv02Mn0etha/r4SGuuqUL6GISD5JW7Z5uuVVtjkMv1VsPNxRfbJM9Ew2EIkSncV9ve8/+IfSh3mXT3Ptyb8DYNPEf+QzvD9sEt1rj/6Ym/8wnSmnTWPjigZMVxv/9ut7+ZePLuGOy8/l8Z0fjLk4SyrcAi7ue7UHerjzyT/x9K4P8RkT2Vbm3r/j8nOpKpvAsnkzImvcoXB/JMi70+YafYtItqSaba7gnWsjbUPLRLb5gOhg1x7o4Wv3vExn0Fmwvqn8eZ5lHge6KwA4u6KbJ756jKr5w+caJgqg7ras+G1e7tatdATu4dz30n5+8fw+gr39AJxaPoHrG2aw7a0P2XcswN2LzuO1Ax08vevDSLAGIsH8yjln0LT0goyeo4gIpB68NW2ea6lOuadZdG3wpqVzuWl9ayRwV5SW8ED3pYCzX7rEGPZ3weI3zmPj+T3DBttk09rDytSsBs6Xh8d3fkCwtz8y4v64u49/feldwtZGss6XzZsBG3bw/ctmRUbZm2/5Mj97bq8Ct4jkHSWs5YPoLVcwZE06ExrnTI0kjS1uepWD7U6t87qaSv7bl+six924oC4mwWwse7Fh6DYvt4Tpxl/8PxltIRqdof7MrfOZXFEKOOvelX5fzJT9RXXVTBmY0geYUunnIq2Li0ge0si7SLmZ2tHT2JMrStm0ch41VWVUlPkAuPkrZwOwcUXDuKqSJdvm9fCxL3J9dR2fytBWOfd83T3dJW4FF2BiqS9yP13dzkREskHBuwANt+48nOjA5gZtV6p7sZOJDqLueblfCD415z+GdlRLY7Keu94ePfKH2L3m8XXT3dfPOa2KQE8f7YGeUV9PEZFM0bR5gYnuvNUe6OG+l/Zz7f3NkXKg7YGeyJ+JprHdv5cJy+bNiFkvH+8XgtGIHvlvu23BkKUAdybCvQ7udTmrupzVW/dErsveoye44qfb1Z1MRHJKI+8CEz2CXLD6RYKhMABnTang6CcnI1PBrx3oGLZaWdZGlSNVp0vj6BsSj/yT/bee7A3z3Ntt+H0lzvX8pxc42ddPv4XaKj8NddWsa059j7qISLooeBeYmqoyFs+dxuqteyKBG+C940GaXnwXgHNOq+KuRZ/norrqUQWzjEi1AE4aDJcJn2xa3e8rIRR2tpi5W81KDPxi6VxWbdipNXERyQnt8y5AztTuKzG9rl2TK0p59vZLMr63elQyuFUsVe5yQ6KZiPLSEroHArdrckXpqFqZioikIt9agkqWtAd6WLVhJ2E7WA40WnRSWt7IwVa5eMnqpn/vq2cT6hv6Jagz2MvkitJhA7ebW+By8w1ERMZLwXs4rz8YW3c80JaWvceZ5CZm1dVUckp57KpIeWlJxpPSvCw+oQ7g0db3Iv3HTy0vjXkt1Bc7Go8WnzjoTssr0U1E0kHBOxm3bGmGiodkyrJ5M7jj8nOxOJXEYPB/8tRTy5lZUzGuYisZk+iL0qYbcvrlacuuI7Sd6KF2Uhn/tuIiplQOBm+/z9AVCif9IhTo6Ytc64VrtvO1e15O2EdcRGQslLCWzOxr4I2HMtdnO4OqyiZwsD3IOadVsXjuNC773OmR5KroZhx5w/2i9MZDg0lr918MgaPw59/Dyled59xENsjKtHp0hvqWXUc42B5kZk0F36w/k2vrz0yanb+u+RCrt+6hrqaSyRWD/cWBIX3ERUTGovAT1saTDBVoG1o85OaWrHX6Go+xFmrJiWSd1Xx+CIdy2m0tWqrX1J0i33csEOkPDk6W+tZbFzDr9En5/f9DRHJGCWvg2anvdMhlQZRRq6p1AnJFjROg3T3fK14Z+lyOAjekfk1rqspoWjoXnzFEfzXut3DT+lb2Hj2h9W8RGZfCDt6zr3FGau7U930Ng3uKR5r6ji8e4gaRR66CV+7xXCJbQfFAIuFzbx+N2arn5vgfbA+yuOlV9h0LUFvl1/q3iIxJYQfvZCO6+NFbomCw5W8HA/3NLc7N/SLw/F1FOZrPmGRflNbOH/rc/Rfn/WyK24YUnP3glX5fzAi8KxTGZwxtgVD+JQ6KiCcoYS1RspS7/jp7MTT+c2yf7T+shzcf9WQiW95KVGXNTVirOt1JWPvDevjjBujY7zzX9g787C+gtDzvrn98B7XjXSGuvPcVesODITxsLXU1lRp5i8iYFHbwTqVu9nBZ5dGBG5z782+H86/PaBesouImFMLgn7ufcAL2lr+FqV9wAvfzd0H12TD/B84xv1sDvV3OLc+uf6I66lNPncifj3fHHGfJz2RREcl/hT1tHj2ii5/63v2Ec0yqU+s5UtBVuqITCt3A7U6B734CZlzsBO0/bnACd8d+J2i/8hOw4eHfO8eik9u27DrCn49344uqbuczhoPtQU2bi8iYFPbI290OFr1VbPnm8dXNzlIXLBis0rW++fCQettQAM0wUtlL775ePgUwsUG7fAqYkoxd/3RpnDOV1Vv3EOjpi2l6UlU2QdPmIjImhT3yhpHrZg+XVR6dxOZKZTSfBuuaD9FQVx3pOf21e17mktUvFlaVrpFmPaJf7z4O0dPMxgc3PJOx659Odz75JwI9ffh9JWxc0cDGFQ34fSUEevq488k/5fr0RMSDCnvknYrRtqTMxGg+TnSHq6alc7nugWY6g70AVPp9xdfFykbXEDdOklpvEH673Ln2ab7+6XbXos/z+qFO2k44xVsAQuF+aieVcdeiz+f47ETEiwq/wloq8qAlZbToCl2TK0r5KNgbGXPmZUvPsUpWWS1R1jk4o20bhopaKD/VWQNv/EneBu1o7YEeFq7ZHimVWl3pZ9ttC6ipKvNWNTwRyShVWBuNPGhJGc1tR+n2jLY4RT7cxwXTFWykJYgtf+sEbp/fyTJfttm5H2yD6lmxgdsDhVsSUfcxERkLTZvnseiWk5+qKOXRm+ZFGowUxMhspCWI2dc4jUkCR2HHw84tHHL2eV997+DfGW6vfvTn5IgbkDu6QjEJa0vWttC0dC61Vf5I97GTvWG6QuFI9TWNwkUkEU2b58BI06TR0+aVfh/+CSV0Bnsja+AtBzqK55d5Ks1hRpp+z3EGenQOQ/yugSvnnMHTuz7EZ0xMOVWAOy4/l8d3fhA5rmnpBTk4exHJplSnzTXyzrJUtn/FV+iKPqaoAneq3Kz0PC2ck6hoy8YVDWzZdYTGOVN5+8gJDrR3Dfl7D7z8Lh9391HqMzy960Muaj6k//ciAih4Z13jnKmsbz4cmSYFZwo1evvXcL/si+qXdxb31Gda/P83tyNZe6AnaaW1j7v7MBApqxro6cvwWYqIVyhhLcvcZLRKv4+OrlBkHbRp6dyYalueaumZKfEJbfNWOZXW3IQ2NylttHv188iWXUc42B6Mqb4WzQ3rM2squLb+zOydmIjkNY28c2BT63t0hQYrhfVby8r1OyJTp0UXpJOJTmjb/cRgffNL73Sec9e5D706ur36eaRxzlR+/sJ+2k70MLmilE+6+4asfQP8+OtfAJxlF/18iIiCd5a1B3p4rPW9yGMDdAZ76Qz2qstUItEB3C2V2tzk3KIbyMy4OKOFczJly64jtJ3o4ZzTqlg8dxqrt+5JeNzK9Tv4VEUpB9uDvHagQ8lrIkVOwTvL3GnSmTUVfDQQtF3fqP90YRRfyYSRktLig3SO9+qnKj6/YefhTp57+9iQ49wveH5fiZLXRCQ9a97GmMuNMXuMMfuNMT9M8PoNxpg2Y8wfB27fTcfnetGyeTO4e9F5PHB9PSVR65yVfp/WNItUdH7DglnOzMHMmgoevalhyD/QULi/cGrbi8iYjTt4G2N8QBNwBTAbWGKMmZ3g0EettV8cuD003s/1kvi2ng111Sx96LVIslp1pZ+uULhwKqdlgoeT0kbD/XL32Mov8ZnaKk6tKI15fXJFafHVtheRIdIx8r4Q2G+tPWCtDQG/ARal4X0LQqLyl0sfbKHtRA+1k8rYdtsCtt22INI9TP2dk8hSN7d84E6HL1nbQmewl+g89E+6+zjeFSqsvu4iMmrpWPOeBrwX9fh94KIEx33dGLMA2AvcZq19L/4AY8wKYAXA9OnT03BquZdsX3dtlZ8N372ouPdxj0YWurnlE7dQj99XQig8WCY3bC1L1rZwSvkEDrYHAe1OEClG2drn/f8BM6y1XwCeBR5JdJC1dq21tt5aW19b650CHMNx93VXV/pj9nU/c+sCZp0+KeY4/RIeQTobyOR5I5Nl82Zw5ZwzCIX7qaupZGZNBeDsTujoCnGwPai1b5Eilo7g/QEQnWn16YHnIqy1HdZadzH3IUD7XCR33EYm918MR98ZXE/f8gPYdMPQY3MU5JuWXsDdi85j08p5PLbyS1RX+iNFW6or/Vr7Fili6QjebwDnGGNmGmP8wLeAp6IPMMZEDw+uBt5Ow+d6QnxHKXcEruS0HJp9jdOZLHAU7v8SNF3orJ37/LD78cHg7AZ5NykuOshnKYDHV9oTEYE0BG9rbR/wPWAbTlDeZK19yxhztzHm6oHDbjHGvGWM+U/gFuCG8X6uV0Q3GVFyWp6oqoXrnwLjAxuG7uOAcdqN1n7WCe7g/Okmxd3X4NzcpDn3mCzQF0ARiaeWoFkwUgtQyYFAmzPi7j4++Fz5FFj1+tB2oyO1JM2w4VqK3r3oPP0ciRQQtQTNI8k6SkmOuNPf7ojbXUk++TF0deRdtzJ1mROReOoqJsXH3TPu8wPWGXG7U+jrrx5MUMujwjDqMici0RS8pfhceCPMXjy4xr3qdVj5+8EkNrfoSxEVhhERb9Ga9zhpPdvDXn8wtuhLoG1o0ZdUjhERSZNU17wVvMdBiUQiIpJOSljLgmSlT1X5SkREMklr3uOQrPSpKl+JiEgmKXinIL6lpzo6iYhILmnafATuuvb65sND1rUDPX08vvODyIgbiFS+0uhbREQyRcF7BMOtawOR0qfxgV0Z5yIikikK3iNw17UXrtlOR1cIiO3oVFU2QZWvREQkqxS8x0mlT0VEJNuUsDYCdXSSQqLkS5HCoOA9ArX0lBivPxhb1zzQlrXe3uPlJl+6XzzdL6Y/evItBXARj9G0+QjU0UkiXn8QtvwA3ngIlm92nnvkKqfWOeR9ydTGOVP5+Qv7hyRf1k4qU1EhEY/RyDsF6ugkgFPj3G1Mcl+Dc3Mbl8y+JtdnN6Itu47QdqIHnzGRokI+Y2g70aNZJBGPUfAWSVVVrTPidluDuq1Cl2/Oux7giTTOmcrMmgrCUf0MwtYys6ZCI28Rj1HwFikiBpPScyKS3xS8RVIVaHPWuN0RtzsCf+Sq2CS2PLVl1xEOtHfhM4PB2mcMB9q7NG0u4jEK3iKp2v3E4Br3zS3OzV0D3/1Ers9uRI1zplJb5SdsbWTbY9haqvw+GuqqI8dp+5hI/ivabPN1zYdiMsjbAz3KIJfhudnks68ZXONevtkJ3HmeaQ4DCWuBUEw53yt+up22QIilD73GM9+fDwyW+IWhRYhEJD8YG5W8kk/q6+tta2trRt7b3e+aqCb53YvO0y8sKVjxX1r3Hj3B0odeo+1ET0xzHfffhprriGSXMWaHtbZ+xOOKMXi7xSn2HQvoF5YUvfZAz5Da/dtuW6B/ByI5kGrwLso1b7fZiFvq1C19qsAtIiJeUJTBWyStPFwyVbX7RbypKIO3fmFJ2rglU93tYu52si0/8EQAV+1+EW8qymzz6F9Y8QlryjiXUZl9jVPr3C2ZCs7eb4+UTI2v3b+u+RBNS+fScqCDZfNmaBeGSJ4qyuCtZiOSFq8/6ATo5ZudwB1sd573UMlUGPz3EL8LIzqxM/o4Ecm9opw2BzUbkXGKni7v6gDbP/hab/fwfy9P18cb50yNTJkvXLOdhWu2R2aoVPtcJL8UbfAWGZfoDmP3fwm6jzvPGx/0diUumZrn6+PahSHiHQreImPhdhjzV4INO8+VT4GVv09eMtWDLUVP9oY5PrD/G1Q6VSRfKHiLjMeE8sH7pgQqq52g3viToSVT87ylaPwujEq/j65QmKvu/R17j56IvP6jJ99SABfJMQVvkbEYrsMYeKLWebz4bWOPr7oYv6+EULifxU2vcsnqF2PWwDUKF8mdosw2Fxm36A5jyzc7zz1y1eB0eaLgHR/wYTDg58HoO34XRk1VGZtv+TKLm16lK+QsDfiMoWnpXEANTERyScFbZCzG0mFsLAE/y+KD8JRKPxNLfZHgHbaW6x5opsSYSD8AZaKLZF9agrcx5nLgZ4APeMha++O418uAdcAFQAdwnbX2UDo+WyRnEq1pDxeAPdZSNH4NvN9aOoO9dAZ7AaeByeK504b8HdVLEMm8ca95G2N8QBNwBTAbWGKMmR132HeATmvt2cAa4J/G+7kinnThjbHT47ufiM00d/d958F+8Pg18EdvmofPmMjrJ3vDrN66J1JWWAltItmTjpH3hcB+a+0BAGPMb4BFwO6oYxYBfz9w/7fAL4wxxuZrP1KRbHD3fb/x0NBpdEj+fJZG6dFr4ACrNuwkbC2Vfh8TS310dIXw+0oiRV0ATaWLZEk6gvc04L2ox+8DFyU7xlrbZ4z5GKgG2tPw+SLelKwuevXZzv08qJceXTo1WT+A8tKSSC/wyRWlLJ47TUVdRDIsrxLWjDErgBUA06dPz/HZiGSYu+87vi76t7c69/OoXnqifgCL505j9dY9hPoGJ9A+6e5j9dY9VJVN0Lq3SAalY5/3B8CZUY8/PfBcwmOMMROAU3ES12JYa9daa+uttfW1tbkvWiGSdX3dTq30PBTfD+Cyz52O31dC2FoMYHCy0f2+EhrqqnN2niLFIB3B+w3gHGPMTGOMH/gW8FTcMU8BywfufwN4QevdUvTi9337KyHUBWvnwy8vc543PudYdz94fL30HGo50EEo3I/PGCxgcfaBh8L9tBzIzy8gIoVi3MHbWtsHfA/YBrwNbLLWvmWMudsYc/XAYb8Eqo0x+4HbgR+O93NFPC963/fNLfCdF8Dnh3AIOg85gduGnTXw6rMT10vPoWXzZnDH5edySvng6tsp5RO44/JzAWfbmEvV2ETSKy1r3tbaLcCWuOd+FHX/JPDNdHyWSMGI3/ddVQsrXoGH/ovTmcyGY9fA82w/eHugh8d3fkBnsJfqSj/gZJv/+ncHaQuEWN98eEhyG6gam0g65FXCmkjRiQ/GldVQWu4E72gjFYDJgeh94PFBunZSmbaQiWSQgrdIvsjz2ufxEmWgb1zRwJZdR2icM5WFa7ZHtpCpL7hIeqmrmEi+iF8Dv7kleW/wPBGfgV5TVaZpcZEs0MhbJF94rPZ5MvE10cGZNl+ytkWjb5E00chbJJ/E1z7Pw7XukcTXRN922wLOOa2KfccCbNl1JNenJ1IQNPIWkbRKVo0t+jV1HxMZH428RfJZHnQXG4votfB1zYdYvXUPj+/8QN3HRNJEI2+RfDVS1zGPTKc3zpnK+ubD2jomkkYaeYvkq9nXDGab39e3OLCqAAAVwklEQVTg3Nxs9Cx3Fxsrd2S9cUUDlX6njaibyNa0dK7WwEXGSCNvkXyVrOtYHu75TmRd8yF+9ORbrG8+zF+edzpdoXDkte7eMDetb+VgexBQ1TWR0dLIWyRfJFrf/sP63J3PODXOmRrJMv/Xl96NeS0YCnOwPcjMmgpNnYuMgYK3SD5w17fdzmGBNvj15fD8XYMV1ypq8rK7WDI1VWWR6fL+JD0Ev1l/pvZ9i4yBgrdIPki0vt2x33mt+mzPVFxLZGKpL3LfRD1f6fdxbf2Z2T8hkQKg4C2SD9z1bXd07Y62L73T6Srmdh1bvhkaf+KJTPPoSmuTK0ojfb99xjC5opSuUJgla1tiWoeKSGoUvEXy2fnXe7biWnSltRsX1BG2Fr+vhLC13LigTlXXRMZB2eYi+cBjHcVSEV9prapsAg111bQc6GDZvBlcW3+mqqyJjJGxNkkmSY7V19fb1tbWXJ+GSHa4CWu1nx1akMUj0+QiMn7GmB3W2vqRjtPIWyQfFEhHsbFY13wopg666p6LjExr3iL5ogA6io2WW8jFTVxT3XPxgnXNh2ISLdsDPVn/edXIW0RyRnXPxQuiZ4fcL5w/f2E/z3x/PgBL1raw71gAyF61QI28RbzKox3HormFXKor/TF1zzeuaFDxFskL8bNDDXXV+H0ltJ3o4ZLVL7JwzfbIropsfuHUyFvEiwqk45hIvks0OxQK9+Mzhq5QmK5QOCdfODXyFvGiAug4BrGFXKor/ZER+BU/e4W9R0/EHKc1cMmFRLNDkytKOaU8t2NfBW8RL0pWkc1je8KjC7lsu20B225bQO2kMtpO9LD0QSWxSW7FJ6a5PunuozPYG/OFM9vVAjVtLiI5E1/IBWDDdy9i6YMttAVCSmKTnHHXun/96kEMho6uED5jCFsbqRa4cUUDUyr9kYS1bG5x1MhbxIviK7J5rONYtGXzZsSsFc46fRLP3LogZpqy0u+jaencmL3gGoVLJrktbQ+2BznQ3hUJ3DNrKqiu9BMK99NyoCMyrX73ovOyWptAwVvEi3Y/MbjG7eGOY6nqCoVZuX6HptEla6LXugHC1lJd6eexlV9i220LYoJ1TVVZ1osKadpcxIsKuCJbfBJbMNRHd28/B9q7+No9LwPQGezVNLrkTC6CdTyNvEW86sIbnWDtTpNX1TrB3GN7veNFJ7F9Z/5Munv7KfU5ncA7g710BnsBWDx3mvaCS8Yk2wmRL21sNfIW8aoC3esdncQG8PjODyLVq1w+Y7jsc6dn+9SkiER/idy4ogEgJ4lpyairmIhXuUlrbe/EthF1O5N5aMvYcPYePcEVP32F8MDvKgNYiPxS1ehbMiUXTXNS7SqmaXMRr0q21/sL18Ue58GyqdGee/toJHADlPt9nDWlIjICUua5ZEr8Toh8WOt2adpcpJD0dsPzd8GbjxbEVHp7oIfHd34AwOSKUkJ9/XSFwvzPj7tZ9dXP0Dhnak6aQojkmoK3iFfF7/UG577PP1g21X3OY2VTXfHrjse7Qlx17+8Ihft5+NVD/Ob191TARYqSgreIV0Xv9Y4fZfsrnaANniyb6oqvwFZTVcbmW77M4qZXI00hEhVwyYeEIpFMGteatzFmijHmWWPMvoE/Jyc5LmyM+ePA7anxfKaIDLjwRmj8yWBgdtfAL70TJpSn9h4eaCsav+44pdLPxFJf5HFXKMxN61tVwEWKynhH3j8EnrfW/tgY88OBx/9nguO6rbVfHOdniUi8RGvYbz46dCr9kauGjr49uNUsfu9tv7V0Bns52B7ka/e8TIkxmkaXojDebPNFwCMD9x8BvLeoJlJIRlM21YNtReO7kD17+yVUVznlKzuDvZF2jSrgIoVuvCPv0621Rwbufwgkq5ow0RjTCvQBP7bWFl7xZZF8MJqyqe40+30Nnlkfj18DX9d8iI5AKLL3G5x2jau37qGqbILWvSUne7WzYcTgbYx5DjgjwUt/F/3AWmuNMckqvpxlrf3AGFMHvGCM2WWtfTfBZ60AVgBMnz59xJMXkQQSBek8nAIfq+hfug111fh9JYTC/ZiB59x2jQ111Tk5P8kfblvP9c2Hh1RJA29vLRxx2txae5m19vMJbk8CR40xUwEG/jyW5D0+GPjzAPAScH6S49Zaa+uttfW1tfn5zV+kYBRAW9GWAx2Ewv34jMHijL59xhAK9/Oz5/bG1KBWMZfi47b13HcswMI121m4Zntk2cXrORHjXfN+Clg+cH858GT8AcaYycaYsoH7NcDFwO5xfq6IjFcBtBVdNm8Gd1x+LqeUD04inlI+gdMm+Xl614eRJhKvHezgK//8krLQi0x0W0+3N3x1pb8gyuqOd837x8AmY8x3gMPAtQDGmHpgpbX2u8DngAeMMf04XxZ+bK1V8BbJtQJoK+pWYOsM9kb6Lnd0hSKv7zsW4NJ/eZmPu51OZFVlEzw/4hIBNSYREQ9z1zSjOz99/b7fc/h4MOHxj97UwEUztRZeLNythfuOBWK+3OVzU5tUG5OowpqIeFZ89jnAv9/8JX796kGaXozNiT21vJTP1FZl+xQlh/K9red4aOQtIgXntYMdXPdAy5Dnzzmtiqalc2k50OHpX9ySOq9tFUt15K3gLSIFJb7/d7zaKj9tgRBXzjmDpqUXZPnsRIanft4iUpR+9txewtbiM4a1yy5gckVp5LUSA22BEH5fCU/v+lCZ5+JZCt4iUlCall7AlXPO4Jlb5/OXs8/g0ZvmRQq49FswQCjcXxB7faV4KXiLFDMPdBUbi6alFzDr9EmA04XsU1GjbwtMrijN22xjkVQo21ykWHmwq9houVuFOoO9Q+qfH+8KKXiLZ2nkLVKsPNhVbLTcrUJ+X0lkxO0zhrC1LH2wJaZ8qnjTuuZDRVkGV8FbpFi5XcXcmuZujfM87io2WsvmzeDKOWdE1rifvf0Snrl1fiTjfMuuIyO/ieQtt0iPWwbXnWkphjK4mjYXkYLWtPQCLmoe3OtbU1XGM7cuyOu9vpKaxjlTWd98ONJ4BAYrqBV6MqL2eYsUK7erWNs7UD4FTIkz+q4+G764FObfnuszFBlRe6CHhWu2R2raV1f62XbbAs/mM2ift4gMz+0qBlAxBZY/7QTujv3w/F0FkXUuUqgUvEWK1YU3wqV3DgbsR66E4HHntQJKWpPC5a5xu60+3daf7hp4IVPwFilm82+Hb28dTFrrPl5wSWtSuKIbj2y7bQHbblvAOadVRRqPFDIlrImIiCcl6iq3cUVDUSQjauQtUszcpDV3m5g7An/kqtjKayJ5atm8GTHJaTVVZQUfuEHBW6S4uUlrtZ+Fm1ucm1u4ZfcTuT67nCvWAiCS/zRtLlLM3BKos68ZXONevtkJ3AVQHnU83AIg65sPs3FFA5ta3+Ox1vc42B4EnKnaYpielfyk4C1S7OKDdFVt0QduiC0AcsnqF+kKhQGYWVNBQ101S9a2sO9YAEABXLJO0+YiIgnUVJWxcUUD1ZX+SOAG+CjYGwncdTWVkUpemlKXbFLwFhFJkQE6g72Ral52oE9ZMdXUTrfovIJ1zYfYe/RE5BrqC1FymjYXEUkgvgBIv7V0Bnsjr/uM4WB7kIVrtnOyN0xXKBypqd0e6NF6eAqi8woWz53G6q178PtKCIX7CfT08fjOD7Q0kYSCt4hIAtEFQJqWzmXl+h2R4F3p99EVCuMzJjIK9xnDX553Ose7QqzasJN9xwIEevqoKpugwJNEdF7Bg9sP4DOGULgfnzE8uP0AncHeomgyMhYK3iIiCUQXANmy6wgH2ruoq6nkG/Wf5tr6M7n2/mYOtHdFjg9bS9OL73L/SwcIW0t1pT8mO10BfCg3ryC6sYjBuZadwV6qK/1sXNHg2SYjmaTgLSKShBtw4yt5tQd6Iuvd7ijcFR7o1NjRFSqa9pSjsa75UExFtONdIU72hof/SzKEEtZERFLgVvJa13yITQMj6nNOq+JX3/5fmVxRmvDvVJSWaOQYxV3jdhuH7D16gqvu/R1doTAVfh8+Y7A4SxCTK0qLpsnIWCh4i4ikyA0+j+/8gDsuP5empXP54b+/GZPIFi3Y28+m1veye5J5rHHO1EjjkIVrtrO46VVC4X78vhK+ffEMwtbi95UQtpYbF9QVTZORsTB2YIon39TX19vW1tZcn4aISISbgb7vWIDqSj9AZK02mZk1FTy28ksARZmBHj9NvvfoCRY3vRpZaqj0+3h81cXMOn0S65oP0VBXTcuBDpbNm1GUWfvGmB3W2vqRjtOat4hIihIlWE2uKKWkxNARGBrEJ5TAwfYgm1rfK+htT/EB2g26QEyJWYCb1rfG5AhMLPUxZeCLkHtdZp0+CSieJiNjoWlzEZFxCPX10xEIUVdTSUVp7K/Uvn5nZPnLVw5Gtp25yWuF0vQkfh07umBNoKcvZpr8a/e8HMm+n1xRSnWlX+vaY6SRt4hIiuILtwAD90sJ91uCvf0xBV18xtAVcgq4uNuetuw6QqCnj9Vb97C++TBNS+fy3NtHPTcyd0fb8TXgJ5b6Iln219afybX1Z8bMVADU1VSyaeU8gMgyRLFNj4+XRt4iIimKLtyy7bYFbLttAeecVkVHVy+Hjwcjzz97+yXMrKmIbBtzbWp9jx89+RaPtb5HdaXfGZH+dDurt+5h37EAVWUTOPeMSXlfHjR6tA3QtHRu5IuK+8UmWZZ9pd/H/ddfQE1VWWQZ4u5F5ylwj5JG3iIiKYrf7w1ERtPRz7cHejAYgJgR+m9b32dmTUVk6hjAje8GCPT08b+tfY2wtXlTHjTRenb8dHi/tUO+qLjHJpqpWLVhZyS4a117bDTyFhEZBXe/t8sNPtHPuxXZ4kfoB9q7+Gb9mUPWxgHc0Be2NlIeNH6dPNuSrWev3rqHxXOnRdas3a1y8evYm1rfSzhToe1f46eRt4hImg03Qg/09BHs7R/276dSHjRZhnc6R7HR69kL12wHiKxnB3vC9EeNtkuApRedxQ0Xz4isY1eVTeDuReclvA4abY+PgreISAbEB6eaqjIa50zlm/f/ftzvHd2Ny92C5QbMRJ89Vom2xlVX+ll43hn84sX9gDPa7ukNE+zt5xcv7qeizDdsgNY0eXqMa9rcGPNNY8xbxph+Y0zSTeXGmMuNMXuMMfuNMT8cz2eKiHjVll1HONgeZPqUcoxJftxI5UEb50yltsofGREvXLOdfccC1Fb5szLFbold3y4r9cU8VoDOvPGuef8J+Ctge7IDjDE+oAm4ApgNLDHGzB7n54qIeM6yeTO4e9F5XP0X04jP73Lro/uMGbE86JZdR2gLhCItSTu6nPttgVBa15LjE87c9ez/eOsoq776GSZXlNIZ7KUz2MvkilLuuPxcbv7K2Wn7fEluXNPm1tq3AcxwXyHhQmC/tfbAwLG/ARYBu8fz2SIiXtQ4Z2pki9XkilJCff10hcJ8qqKUC86azB2XfzZSHvTa+jMTTj83zpnKw68eGtKStK6mMq0j7+itcfHT84vnTqMk6nd/iTFcW39m2j5bhpeNNe9pQHRl/veBixIdaIxZAawAmD59eubPTEQky4YLiN++eCazTp+UUnnQ+KnrZM+NR7LEO7fca/wWsCVrW9RFLUtGDN7GmOeAMxK89HfW2ifTeTLW2rXAWnAak6TzvUVE8sFwmeiprhO7a+fuFDs40+0H24Npz+ROlHhXVTYh6RcQZZJnx4jB21p72Tg/4wMgei7l0wPPiYgUpUQBcTQBr3HOVH7+wn7aTvTEjHxrJ5VlJWEtHV9AZHyyUaTlDeAcY8xMY4wf+BbwVBY+V0SkIG3ZdYS2Ez1Dip+0nejJWvGTZMVqJDvGteZtjFkM/ByoBZ42xvzRWrvQGPO/AA9ZaxuttX3GmO8B2wAf8Ctr7VvjPnMRkSKlka8Ym6AebT6or6+3ra2tuT4NERGRrDHG7LDWJq2b4lJtcxEREY9R8BYREfEYBW8RERGPUfAWERHxGAVvERERj8nbbHNjTBtwOI1vWQO0p/H9ipWu4/jpGqaHrmN66DqOXzqv4VnW2tqRDsrb4J1uxpjWVNLvZXi6juOna5geuo7poes4frm4hpo2FxER8RgFbxEREY8ppuC9NtcnUCB0HcdP1zA9dB3TQ9dx/LJ+DYtmzVtERKRQFNPIW0REpCAUXPA2xlxujNljjNlvjPlhgtfLjDGPDrz+mjFmRvbPMv+lcB1vN8bsNsa8aYx53hhzVi7OM5+NdA2jjvu6McYaY5Txm0Aq19EYc+3Az+Nbxph/y/Y5ekEK/6anG2NeNMb8YeDfdWMuzjOfGWN+ZYw5Zoz5U5LXjTHm3oFr/KYxZm7GTsZaWzA3nJaj7wJ1gB/4T2B23DE3A/cP3P8W8Giuzzvfbilex68CFQP3/0bXcfTXcOC4ScB2oAWoz/V559stxZ/Fc4A/AJMHHp+W6/POt1uK13Et8DcD92cDh3J93vl2AxYAc4E/JXm9EXgGMEAD8FqmzqXQRt4XAvuttQestSHgN8CiuGMWAY8M3P8tcKkxxmTxHL1gxOtorX3RWhsceNgCfDrL55jvUvlZBPgH4J+Ak9k8OQ9J5TreCDRZazsBrLXHsnyOXpDKdbTAKQP3TwX+ZxbPzxOstduB48McsghYZx0twKeMMVMzcS6FFrynAe9FPX5/4LmEx1hr+4CPgeqsnJ13pHIdo30H59umDBrxGg5MqZ1prX06myfmMan8LM4CZhljXjXGtBhjLs/a2XlHKtfx74G/Nsa8D2wB/vfsnFpBGe3vzjGbkIk3leJhjPlroB64JNfn4iXGmBLgHuCGHJ9KIZiAM3X+FZwZoO3GmDnW2o9yelbeswR42Fr7L8aYecB6Y8znrbX9uT4xGarQRt4fAGdGPf70wHMJjzHGTMCZHurIytl5RyrXEWPMZcDfAVdba3uydG5eMdI1nAR8HnjJGHMIZ33sKSWtDZHKz+L7wFPW2l5r7UFgL04wl0GpXMfvAJsArLXNwEScmt2SupR+d6ZDoQXvN4BzjDEzjTF+nIS0p+KOeQpYPnD/G8ALdiDTQCJGvI7GmPOBB3ACt9YYhxr2GlprP7bW1lhrZ1hrZ+DkDVxtrW3NzenmrVT+TT+BM+rGGFODM41+IJsn6QGpXMc/A5cCGGM+hxO827J6lt73FLBsIOu8AfjYWnskEx9UUNPm1to+Y8z3gG042ZW/sta+ZYy5G2i11j4F/BJnOmg/TuLBt3J3xvkpxev4z0AV8NhAvt+frbVX5+yk80yK11BGkOJ13Ab8pTFmNxAG/tZaq9m0KClex/8DeNAYcxtO8toNGtjEMsZsxPmiWDOQG3AnUApgrb0fJ1egEdgPBIFvZ+xc9P9GRETEWwpt2lxERKTgKXiLiIh4jIK3iIiIxyh4i4iIeIyCt4iIiMcoeIuIiHiMgreIiIjHKHiLiIh4zP8PMmcZLkfgv5IAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# make a dataset with two outputs, correlated, heavy-tail noise. One has more noise than the other.\n", "X1 = np.random.rand(100, 1) # Observed locations for first output\n", "X2 = np.random.rand(50, 1) * 0.5 # Observed locations for second output\n", "\n", "Y1 = np.sin(6*X1) + np.random.randn(*X1.shape) * 0.03\n", "Y2 = np.sin(6*X2+ 0.7) + np.random.randn(*X2.shape) * 0.1\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.plot(X1, Y1, 'x', mew=2)\n", "plt.plot(X2, Y2, 'x', mew=2);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data formatting for the coregionalised model\n", "We add an extra column to our training dataset that contains an index that specifies which output is observed." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Augment the input with ones or zeros to indicate the required output dimension\n", "X_augmented = np.vstack((np.hstack((X1, np.zeros_like(X1))), np.hstack((X2, np.ones_like(X2)))))\n", "\n", "# Augment the Y data\n", "Y_augmented = np.vstack((np.hstack((Y1, np.zeros_like(Y1))), np.hstack((Y2, np.ones_like(Y2)))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Building the coregionalisation kernel\n", "We build a coregionalisation kernel with the Matern 3/2 kernel as the base kernel. This acts on the leading ([0]) data dimension of the augmented X values. The `Coregion` kernel indexes the outputs, and acts on the last ([1]) data dimension (indices) of the augmented X values. To specify these dimensions, we use the built-in `active_dims` argument in the kernel constructor. To construct the full multi-output kernel, we take the product of the two kernels (for a more in-depth tutorial on kernel combination, see [Manipulating kernels](./kernels.ipynb))." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "output_dim = 2 # Number of outputs\n", "rank = 1 # Rank of W\n", "k1 = gpflow.kernels.Matern32(1, active_dims=[0]) # Base kernel\n", "coreg = gpflow.kernels.Coregion(1, output_dim=output_dim, rank=rank, active_dims=[1]) # Coregion kernel\n", "coreg.W = np.random.rand(output_dim, rank) # Initialise the W matrix\n", "kern = k1 * coreg " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** By default, the `W` matrix is initialised with zeros; however, this is a saddle point in the objective, so the value of `W` is not optimised to fit the data. Hence, re-initialising the matrix to random entries should give a more accurate result." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Constructing the model\n", "The final element in building the model is to specify the likelihood for each output dimension. To do this, use the `SwitchedLikelihood` object in GPflow." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# This likelihood switches between Gaussian noise with different variances for each f_i:\n", "lik = gpflow.likelihoods.SwitchedLikelihood([gpflow.likelihoods.Gaussian(), gpflow.likelihoods.Gaussian()])\n", "\n", "# now build the GP model as normal\n", "m = gpflow.models.VGP(X_augmented, Y_augmented, kern=kern, likelihood=lik, num_latent=1)\n", "# Here we specify num_latent=1 to avoid getting two outputs when predicting as Y_augmented is 2-dimensional\n", "\n", "# fit the covariance function parameters\n", "gpflow.train.ScipyOptimizer().minimize(m, maxiter=notebook_niter(1000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's it: the model is trained. Let's plot the model fit to see what's happened." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAD8CAYAAABevCxMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4HOW1/z+zVVvUV7LkKrl3cAEXjCkhsTG9GYgTQwI4uZCEhCQkNwUSbpIL/AghCXaIIRQ7YFouzZgOBoN7l7tkS7blImnVd1db5/398WqLZLmA2sp+P8+jx9LM7MyoeL9zznvO92hCCBQKhUKhUPQcDN19AwqFQqFQKL4cSrwVCoVCoehhKPFWKBQKhaKHocRboVAoFIoehhJvhUKhUCh6GEq8FQqFQqHoYSjxVigUCoWih6HEW6FQKBSKHoYSb4VCoVAoehim7r6BE+FyuURBQUF334ZCoVAoFF3Chg0b3EKInJMdl9TiXVBQwPr167v7NhQKhUKh6BI0Tdt/KseptLlCoVAoFD0MJd4KhUKhUPQwlHgrFAqFQtHDSOo1b4VCoVAojkcoFKK8vBy/39/dt/KlSUlJoW/fvpjN5q/0eiXeXcHaJ2Hk1eBsLiD0VMGO1+HcO7r3vhQKhaIHU15eTmpqKgUFBWia1t23c8oIIaiurqa8vJzCwsKvdA6VNu9s1j4Jy34Gz10uRdtTJT9f9jO5T6FQKBRfCb/fT3Z2do8SbgBN08jOzm5XxkCJd2cz8mrIGQ5Vu2DBZPlRtUtuG3l1x15r7ZPy4SCKp0o9ICgUitOanibcUdp730q8OxtnDtyyFOwu8Lnlh90ltzlP2od/6qgIX6FQKM4YlHifLnRlhK9QKBQKhBBMmzaNd955J7btlVdeYebMmXz3u98lNzeX0aNHd8q1lXh3NtEIOBpxRyPwaITcUXRVhK9QKBQKQKa+n3jiCe655x78fj8ej4df/epXzJ8/n1tvvZV33323066tqs07mx2vxyPgW5bKbc9dLrepinOFQqHo0YwePZorrriChx56CK/Xy9y5cxk0aBCDBg2irKys066rxLuziYpzYqvYLUs7XrhbR/gQj/BV9K1QKE5zCn75dqect+zBy056zP3338/48eOxWCxdNo9DiXdX0FqknTlfXbiP1zMOx4/wNy2G8+9pebyK+BUKhaJDcDgc3HjjjTidTqxWa5dcU4l3TyJaUb7uqbhAP3sZuHfDJb+Hr90PgUbwVoIjF0ZcIcV79T9g3LcALS7ooARcoVCcNpxKhNyZGAwGDIauKyNT4p3s6DoEGsBfD66hkN5Xiu/fzgIhIOQDezZ4KqSIl3wEG5+DsTdC0AuaUYr5oyPBaIGQF9J6Q1pfqNwF1lRISQers+PuWTnKKRQKRaeixDtZWPskjLwKDCbwVUPFDtj5JmQPBj0EGEDoMOpaGXkHvfJ1ZjuMmQ1oYM8Ba7p8/bp/yf0iIvfpIflhsoNrGJR+CofWy31BD9SUwsgr5cNBar4U3pQMOJGRQFsiveznsOO1ltmBNqL9RavKmDUmH5dTppjcngDLio4wd0rBKf/IOuIcCoVC0RncfPPNLF++HLfbTd++ffn973/Pbbfd1mHn14QQ7T+Jpj0NXA5UCiGOaWrTpJXMX4FZgA+4VQix8WTnnThxouiqxf9uw1cDqx6HFX8GhwtG3yC3F70sRXj4FeAaAp6j0FgBjUehaqcU8vZgtEDvcTJq378K/LUwbBb0GtMs+AJMNsgqgKzBUtQdrriYR1P4bayxRxy5GL2VYHehCzA0uePHOXNYtKqM+97YzpBcJ0vmTQbg5oWrKa708MBVo05JfE90jntnDuPOCwcDStAVitOZnTt3MmLEiO6+ja9MW/evadoGIcTEk722oyLvZ4HHgUXH2X8pMKT5YxLwj+Z/Tx9OJVW89kkYcaUUx+oSKPsCytfKSNjuAq8bNjwjhTkSlCnvXW+d/NoGkxRXS6oU18odoIdbHaQBCQ9qkSAcXNPykLLP5X0ZTDDoYrBlQdGrkFkIFgeYHZCWL6P0CXNldB01hQHwudmj9+Eh8d88afsdBp8bA+AWafw/5x/5OWm4gFlj8nnmi1KKKz3M+MtnAFR7gwzJdTJrTP4Jv9VotD1rTD6LV+2nuNLDBQ9/QorZSLU3CMCr68uZPbEfz3xRyovrDlLtkdtnjclXQq5QKE4LOkS8hRCfaZpWcIJDrgIWCRnmr9Y0LUPTtHwhxJGOuH6XcaJK79aFZK1TxSsehY9+D58+KKNrzQBbX5LtXHUHwNj8qwgnGNWLCBiMsvjM2QvCARl127Jg4m1ynXr1P2RU3mcC5J8Nq+ZL4bY4IRKCSECey5YF594u0+k+97Hfm8Es19YDDfLryh3EBN/sgFHXyO/38Ab5+sMbYdKd8NH9sfPpZge/0X/O3hojdSlBsppPrQEf7qxk48LVsShZQ0bwUcHNdlhYMm9yLAXeFtFoe/Gq/SyZN5n5c8Zz6WMr8AYjeIMRAMxGjX1uL9Mf/gRf8zaHxciwvNRYZO4JhHFaTUrEFQpFj6VD0uYAzeK99Dhp86XAg0KIz5u//gj4hRDimJy4pmnzgHkA/fv3n7B///4Oub92c4I0MV+7Xwpx1a6WPdauoTDzYajeAzX7oOgVmQo3WqTAnij1rRlh0vcha6AU8ChlK6RIW1Obr1MDhzdJ8T6wCorfk6nwMbNloL3pOVnUBnJ9XABhn4zQE3/3JpvcFvJJIddDx783azpMuRNCTbD2n/Lf6O3puegYGWg4QlgYMGnye9xHX27w/xocOehCUOsLtcgF2MwGrh3XF384Qr0/RFMwQjCkE4joBMM6QkBE6ByobiIY0TEaNDQE4YQfYavcQguMmkZECApddjSkwLdO0d/1/AbuvmQoQ3vJn+2eikb++uEe5s+ZcPyfhUKh6DZU2jyJEEIsBBaCXPPu5tuJM/LqNtPE5AyHcd+WHwsmx6NaaxoMukSKrdEs16uNFrkvEoyf15YFTTUyEhe6FNiwX0bdGxfBBfeCKUWKaqgJMgvA31x9DjIqHniBvN6EW6XYD7oIbNlgMMCwGbDjTdj2ajyqjgq32dYs5k3yA2TEPuUuwuufw+Q9im6wYtADLX8WgXr47BH5UBH2o5vsCCEwRpooMFQCxIR7r56HhsZAQzmzjGtY7P1G7DSJv9ymkM7zaw+c8q8joh/7p3GiP5aIEGhAnS9ErS9ETqq1RYr+ruc38HbRUT7YUcnSH00D4PK/fU4wosPzG5SAKxSKpKKrxPsQ0C/h677N23oOUe/wRIFO9A6vPwR6JH58JAT+Otj/BVRsay4Ca4XZIVPqqxfI/mytuUdQ6DLyDjbCvk+g/1RIzZMtXs5ekJImRdaaKh8MEhk+69jr9BotK9ejCCGr2L/5MtQdQF98LQZkCKvrOrWNPm5u/BUvaL/BpTe0OJVHWKk196Jf+EAsOveE4K3IFM43bKW/Qf5sTJqOW6QxO3g/gBTuyDdoi0y7GW8wQjCsMyDbTr9MGyajAZNBw2jQMGpa7IHDH9JZUVJFKBKX6pH5aQQjEUrdXiInSGYIiEX8VY0BXl5/MJY+v/uSoXywo5JgRI+twwvAYjRw9yVDj39ShUKh6Aa6SrzfBH6gadqLyEK1+h633n08/HWwfyW8fY+s2DbZQYRlJLvlheaDNJlC91TK4y3NPdVBD6x/BsbcABsXSzEM+aQoX/E3WTw25S5ZLPZVOZ5tKvDGhlIu2/0rTOj4SMEnLLjCDdRtfA09+EMM1rgS6kKjDgdZmofDwQiXBf/AD82vU6BVMNxwkDmmjwGICA2jFhdWDR03mbyszaCt2DjdZuKacX0wGTVKKr1MHJCJ0ahhMmiYDQaMRg2z0YDRoBEIhfn7x3sJRQSpKSaCYZ1AWKfOF+SO8wfy5/f34AnKQj2zEUJtPC+RcBePfVBMMKJztN7PvTOHs/RH05jxl89i+zVg6Y+mxdLoCoVCkSx0iHhrmrYEuBBwaZpWDtwPmAGEEE8Ay5BtYiXIVrHvdMR1u5TWIigi8vOnZ0DuSPm5yRbvp47Sa5QsUKsogm3/kZHz+LkQ8Mh2MG8lhAKy+KypVr7GaIWC82H0te2/7xMMRjGseBiTcTfh7GEcufwlbl64mn+b/4ehhnJesP6JLM0DxFPgtbqTOuFkqOEQN6du5ovUmym1NnFI7GNE4xryvdtjwh0RGi6tgSWWP/Ld8K84GM5iQLqRuqBGfVO8En7e+QO5y7m8uRBwZPxn3Yapy6JVZRyp97fZHrZ+fw2eYJgcp4Xn75jM+9uP8sj7ewCwmAwEw8eG5MHmMH3B8r28sr6cLIe5xeOFAN7cfJifzRjWvt+BQqE4LRFCcP755/PrX/+aSy+9FJAjQZ944gl0XaeiogJN05g3bx533313h167wwrWOoOk6vOOFqxl9JdiHPLB1pehqRqcebLiO4prKPSZKNeuC6fLbUEP7P1YtoXlDIf8s6Rg73yr7WK3hL7oL0tr8xLPigUsi0xi9oXN67aeKjybXuGadaOYXP1/rLRM40gkFV8wQjb1/M78LFcY17BH78P/RabxYWQ8Cyx/Y6jhEB/1+g4OY4QGa29y9SpS7WZsadmk127Hsed1KvR0hMFMHvGK9g8i43nXdQvfG1DBzRtHUB2UBXgOi5E1X9+H86P/brsQcNYjbQr48YxZWu9bsLyEjftr+XBnZaxgLRGDBn0zbRyp97dIw7dm+hAXj954tjKDUSiSjGQoWNu2bRs33HADmzZtIhwOM27cOJ599llsNhvjx4+nsbGRCRMm8PrrrzNy5MgWr21PwZoS75MR9MlK8fL1sOMNKTJpvaFyuxRjT4U8zmCGvudIsU7Nk9v0sOzd1kOyMG3AVGm4Yk1Iw56oir0N8ToZrc1LXl5/kFfXl8eqq6O9zjef258law7wwNIdhPVjRW2O4X3ejUyiinRsZgPZWgMXRVbyafrVLJw7gUE5TszBBqjeK3vVvVVwZCsvi0u4uMCCq3YLFL8PjYflj8Js5xn9Mh71ziTFbALNQHXQyDmuMEusf8RUvbvDHl4ScXsCXPrYZ1R5gmTazdT5QrHoOt1mYsE3x1Ne5+PRD4qpaAgc9zyF2XZe+a+p1HiDzHlyNVWe4CkbyigUis4hGcQb4N5778XhcOD1eklNTeW3v/1ti/1XXXUVP/jBD/j617/eYrsS744m6IW6g7LQrHKXLCCzpkrzkgMrpZlJtHLbmgYF06QwR9eyQ37wVcmis95ny9autN7HtxrtQC9wtycQSyXbLcZYr/OALDuXn5XPm5sPc7C26bip5CvG5uFKTeGVdeV4gmGyHBbe+sF5WM3G47ugCSGd345skf3fehhsmbIgz71Htq/V7AOgHiemwRfRlDuOm9cMoNhr4+HpVmZv+17LQsA7V3fIGNPow8xAlwOBoNTtaxGFRx9yUlNM3Pr0OjYfrKUp1HbVm6G5SF9ALD2/el+1EnCFoptoIX6/S++ci/yu/qSHeL3eFiNBEyeLlZWVMX36dLZt20ZaWlqL151WrWLdgh6REXLDYTi6VZqmRFupUvNkL3XJh9KRLLqenZoPAy+E3hPiBiv+hnhB2pAZkDf61IrNOnBkqMtp5ZrxfXj43d0x4QbYX+Nj/id7Y18Hw7JXOqILrCYNXUAoItha3sDi24Zz10WDuf+Nbfz+qtGxdPGSeZPbThdrmnReS2v+mbj3SO/0+gPgzIUpPwT3Hiq3vkdu0z4oeQvHgY95vfBrvBGZymxR2tKcpgOJ3qsnEObhd3czJNfJ/Dnj+XBnBa9tPERxpSf2Pf39m+OY8ZfPaArJVj67xQBosZ9jYoJidJ80bn9uPQdqfC2uk4jyXlcozgyONxLU4/Fw3XXX8dhjjx0j3O3lzBPvcCDeJ+1zg7sE6vbHDVOsqTJKRpPR4vp/yagymmzNHiwFfdR18Ui6qVa2ejl7wdjZkD0kLujdwCUjevHn9/Ycs8YbZXSfNEwGjc0H6+mVZuXxb46nd3oKtz6zjuJKD8v3VDF3SsExvc0up/XkwmNOgfyxsj2tthRKV8ifry2L3IvuksK+512o24+j5E2+af4QDhlkTYAppXkwixuemSl/jlf+vd0Zieg9O62mmJgO7ZXK7In9YmIazVhUe4NkO2Q/frU3SEG2HbcngCfQsnT9k90yS2AxatR4gtT6gmTaLbH90b7xqBtcYro98Z4UCkUHcQoRcmfSeiRoKBTiuuuuY86cOVx7bQcUH7fizBHvkF/6hnurmg1RhBRfi1OKbtTFLBKGQxug9DOoPxh/fe/xMjVe9CpUfw6OXpA7XBqnZAyAUVdDRoE0RulG3J4Adzy3/rjCnWLSmDkqj+lDclhTWsM14/ucPLL+KhgMkD1IftQfkmY1VbvBngVTfwTVxbD7HSnsUTQjjL4edr4h2+SqS+CJ8+D7X8j97ZxF3vr7SnwYWVZ0hOJKT5uV7CeYq0YwInjso2Ke+nwfUwe5+N70QRQdquPtoqNYjAaKKz1Mf+hj/GEdXch0++SB2SxaVaYEXKE4TRFCcNtttzFixAjuueeeTrnG6S/e0fVkk1UKd0oGHNksW7ES8TfI9ez9X8goGqTbGZqcge3eIz+CHrkm68iRw0AKp8t/k4Dth+v50ZJN7G9O5bbGbNTwhwVvbD7MTef2545+A1vsP6XI+quQ3gfOuklmMMo+l7UEtsy4iG99SWYvQl7Y+Kws/gMp5p4KmH+ufOCKFrKNvLrDbzH6fSemuefPGc+1C1biCYQxahrOFBP1TfE2QJMBxvbNYG+Vl/qmEO/vqOCDHRUMzHGQ47TEomxf8xq6QYPH54znruc3UlzpaXFdhUJx+vDFF1+wePFixowZw9lnnw3An/70J2bNasNE6ytyeot3tJJ73VNw0xJZOb7p+Xhb14BpULNXCvbhzcRS46m9od+58vM+E+DTh6Rog1wHL5gG4+bI9DC0q8CsI9hb5eEPS3ewfHcVAmkuYjYaYn3MUc/vfpl2BKLFOm+XkpoHY66XDzz7V8KRrZCSDhf+Wgr6pkVy+SJaV9B/isyCNNXIr21ZsgJ9x+sdVuCXSOufx+p91XgCYXJSrTz+zXH88j9bY+JtMWoEI4Iab5CfzxjC+rI6ig7VU+r2srdKzlpv7bWuC/j+4g3SntVpOekENYVC0XP43e9+F/t82rRpdHYx+Okt3ol+5P+6JO4P7siVafRPH2rZn220wtlzpH/46gVy3zEDRLS48cnxJoh1EYfrmvjD0h28t71CendrML5fBkaDxrqyWvpm2pg9sS8zR+fHor17Zw7r/olazly5zFA4HT56QGY9LPa4h3uU/Z+3fF3IB0t/DLuWwpon5O9q3Lc77eefGI0vKzpCqdtHocvODRP7MXtiv1haXRcaf7hmDHuONvDJrirWltVQdKi+RcFglFpfCIMGVZ6gKl5TKBRfmdO/VcxT1dKP3ND8vBKdd21Nk+vZFdtle1eidak9WzqpNdXJNK9mlOcxWuRwkU7oS26L1lXL+91efv7qFjYfrCPYbC4ytk86Y/qmk5tq5bzBLraU13HV2X2Su9I5mhnJKpSFhA2HOfFsMI7db3bIdHsn/vyjnEr1uK4Lth+u55tPrqYxcBx/VmTr3o3n9uPOCwd32v0qFKc7ydLn/VVRrWInQw8f+7lrmEzL5o2RxWqDv3Zserz3eCj5oG0DFbOj7QElHciiVWWxFqfFq/bzt5vH8eA7O1lR7I61LQ3rlcqEARlkO61MGZjNuP6Z2CxGJhZktThXp61nt4fEzAjIhyMRkenxcFPzqNHWYt5K2EPeeDq9E4UbTlzwFsVg0Nh0sI7GQASX08KdFwzif9/ddYyD28FaHw+/u7v7syAKRQ9HCIF2PA+NJKa9gfPpLd5RP3J/XcvtdheM+1ZLpzOghTAYrbJNadfSluurtyyFTYth1XwpHJ1E1Fyk0GWnMNtOcaWHS/+6IrbfZNCYOboXA7IdjO+fybmFWaSmmI9/wmSk9aQ2EZG/m3nLoXIHfHC/zIxUFLVMp7cm3ARv/xQu+Ll8GINurUNonW5vy3pVF7KAbVt5PW5PILkzJApFkpKSkkJ1dTXZ2dk9SsCFEFRXV5OSkvKVz3F6i3fiUI6Lfyvbv3a8IdeyEyvO/fWw8u/SWS0lvbnXuBoWX912RLf1pWOndD13eYdGf7PG5LN41X6KKz0YW/1Nmo0aN5/bj3MLspk62EWWw9L2SXoqphQYOkP22b/7Szhhs5YmI/Sdb0DZZ3DRb2VL36u3dksdQpREAf/7xyVUNQbaLGB7eUM5r28+xKLvnMvgvFQu/esKqhoDrNlXzfw5E5SYKxQnoG/fvpSXl1NVVdXdt/KlSUlJoW/fr96pdHqLd/RNO9oqdniTHLHZWrj3LZcCnD0EvvOO3B5Nj7eO3k4wpetUI71TWTsNhXX6ZNgorvTQOnCzmAzcdM4ARvbuWMeeLud440qjD0KFF8XrC9rCYGq5JNJUC+//Wmp9qEn6yHdCW9mXYVnREaoaAxS67NT5QtT64q1mualWKhsDBCOCm55ag9kAUWfWt4uOMmp5ScwF7qV1B3j7R9O76btQKJITs9lMYWFhd99Gt3B6izfExdQf9SJPlcIdCcq+Y4cLLv+LtD9tnR5vS4wTHwhOdmwbRNPhUectiJuBANwwoR+PvL+bRavK2ky3pqWYaPCHufvFTSyZNzn2ANAjOdmDEDQLdxtFbGa7TKVnFkqhjnYNhJvi+wd/A9zFMpti6p7sxNwpBXgCYV5dX06tL0S2w4IuBLW+EKkpJmaO6sWi1QeAuHBH+cfyEhr9ETRg++FG7np+wzGudwqF4szk9Bfv1gghZ2hHwjB0puzjNpq+nL94O7zIE9PhFzz8CRaTgVpfiME5DqoaA4z/n/djgzEsCb3aBk2mWbOdFlxOa/f1anckp/IgFGiUSxrRXm/NKM1deo2CpXdLC9YRV8miw+2vERP5cEBG5ns/lFPPhl0qR7V2w7qY02pin9vbpnvbdRP6kpZipMF/bGV6Y/M2gfxbmDu1QDmzKRQK4ExoFYvib4DPH5XinTNMCrc96+Sv6wTcngAXPPwJ3uY+YKfVSIbdQnmtjBrTUkwUuuxsKW+gb6aNGyf2Y8bovOTq1e4Komn1ql2yotxfL4vaNKOMpptqIL0fDL4Edr4p6xRapNk1GPINuQbeVAs5I2DoN7rl997WUsnL6w/G0uInapDTgBe/N5nfvLaN4koPl43JUxG4QnGa0qWtYpqmzQT+ChiBp4QQD7bafyvw/4BDzZseF0I81RHXPmWMZnDmS3e0XiO7PAJr/eZtSqhC8wQieAJSuB0WI9eM682Fw3IpqfR2nvd4T6B1Wt1bDU9Mbe69r5Hr5Hd8AluWSOG2ZcPE78jXrnpcTiorfg+ObpGe6XVlsGoBDLkE+kzs0uExbbWZOa0miis9sQxLus3cwn41igC+88xafEEdi9HA20VHmbRKReAKxZlMuyNvTdOMwB7g60A5sA64WQixI+GYW4GJQogffJlzd9s87w4mus4dHUd5x6J17K9uOuY4q1Hjt5ePZNbY3qdfBflXJXHWuadK+pxHU+iJc7/XPgl9z5FjWxuPyD78fZ/IYSjRfvwBU+WoVn+9dHkbcUW3+9JHp48NyXXGRrmejOjfkZolrlCcfpxq5N0RI7DOBUqEEPuEEEHgReCqDjjvacOsMfkMyXVSXOnhir9/3kK4DQkJAIPBwMwx+Uq4Ezn3jrhwP3d5POK2u+KV6Z4qeVzvs+Gc2+UauIjILMv0n8vUuWaUfuorHm3uMPgM/m8ebH1FWuV6quS+tU+2vP7aJ+W+KJ6qY49pB/PnTOCBq0axZN5k7rxwMKN6S+8Bo6bx7++eg7nV/1CzUePRG8Zy1/Mbue+N7SxaVdZh96JQKHoOHZE37AMkzM6kHJjUxnHXaZo2HRml/0QIcbCNY9A0bR4wD6B///4dcHvdj8tp5aFrxzB74WoC4fiwELvViDcQIS3FRFgX+IIRbl64uudXkXcGp9qiZzBCn3FyFGnx+3C0CAovgPxxsPVFOYJ0/b/i563cIYef7F7WcjzpuXe0HGzTiT72idHz2z+azl3Pb2Du1AJ+89q2YyrQQxHBVQtWogsZgavhJgrFmUlXDZ9+CygQQowFPgCeO96BQoiFQoiJQoiJOTmda3fZFQRCER5Yup3ZC1cT1uNLFEYNvIEIg3KcfPyzC/ns3oti0fmyoiPdeMdJyrl3wKxH4kY4UXe2WY+0LaIpaXKC2dnflGvfQoepP4RR18iitihNtbDmH3HhzhoY7w0febV8WKjaJV3gFkyOP0B0Yv/4/DkT2H20MbYeDrKIMUr0z2jigEyspu6dH69QKLqHjljzngL8Tggxo/nr/wYQQvzvcY43AjVCiPSTnbunr3l/sP0ov359G5WNgdg2q0nDZDDgDUbIcVp4/o7JDO0lU6XKTauTCPpg78dQvk4Om4mEZJFbdXHL40w2mPxf0jo3e5Dc1nqwTeI6eyeTuB4+f854rn9iJQ1N4RbHjMhP5X+uHMWEgqweZQ+pUCjapivXvNcBQzRNK9Q0zQLcBLzZ6mYSc3tXAjs74LpJwaJVZbg9cXF2ewI8/vEe5v5rDXcs3kBlYwCHxQhAlsPCE3Mm8PHPLmRIrpMqT5DV+6pjr03K4SGnAxY7jLgcJtwiHdnCTXKcqLHV0oQeBoMZNi6CnUubB6N0H9H18PlzxnPX8xtpaAqTaTeTYoqL9M4jjXznufX8+f3dbVaqKxSK05MO6fPWNG0W8BiyVexpIcQfNU17AFgvhHhT07T/RYp2GKgB/ksIsetk5032yDuxinzJvMnous6lf/2caq/sMzYaNM4blM3QPCfVnhA/vHgwA3PkyFEVZXcToSbY/jq8+9/grz12v8EsLXSNZtAMUPQq1OztsvGvbdH67wxg9j9Xsa/Ki9NqxNM8enTigEx+fdkIzu6XoaJwhaKHcqqR95lj0tIJuD2BmFNWWooJTyAcW48ckG1n+hAXOalWvjEyj2F5qeoNNVmIFqKBnNM+8TZYuxACzRa6BpNcGw94YM87kNEfbvtAVqxHC9aOt9beSbRl8vLWlsOxxlw+AAAgAElEQVQUZNv5+8clbDpQhwBcTgu3TyvkW1MKcFrPPANFhaKno8S7gzneMJGpg7KY9dfPCSb4kH9jZC8G5tg5q28m04fm4FBvosnH2ifBVysL1Cq2ySr1I1tk9XlDs5dQ7kjIGCDFO2+M7AuPhLpt1OjxKHV7+cfyEt7ddpQGfxijQWP6EBe/unQ4Q/J6+PAaheIMQ4l3B9JW2vKmf66ipMqLxWQgGI7386SYNOZNH8T1E/rRP9veTXes+FLUHZT2qt5qSOstHdmKXpEpdmsqnPXNuKXqqGukvW6SUd8U4v82HuSldeXsOtoIQL9MGz+4eDBXnd2HFLOxm+9QoVCcCkq8O5DE9Hi2w0IootPgb1n1azUaMBigKaQzONfJi6pXu2cRDsL+z6F0BVjTAAGb/i3XuwEKpsOgr8l18n6TYPDX5JjZJCKiC9aV1fDMF6V8tsdNUyiC1WTg8rH5/GLmcHLTUrr7FhUKxUlQ4t3BuD0Bvv7opy3mMZsMGmFdkGEzcedFg5k+NIcfvrCJ4koPD1w1ShWj9UTqy2VBW1MtOHtB6afSwEXokJovU+muoZDRD0ZfJ9fBkyyNfqS+iSVrDvD65kMcqJEV88PynEwuzCY31cqN5/bH5bTGhqOcEUNuFIoeghLvDsQXCPPI+7t5+ouy2DaTQePGc/pyuM7P7dMKmTLIhcGgqSry04FwEMo+g9LPZbq8qRY2LQZvs02qJRXOnSf7x/e8DXUHuryA7WT4QxGW767k36v3s3pfTQuDoIJsOwvnTuR7i9dT6vYBqIdNhSJJUOLdAYQjOq9vPsT/e3c3Fc1GKwYNzEYDgbBOfnoKL82bTP9sR7fdo6ITqSmVM8JDfrBlyHXwwxvlPs0o0+YhHzjz4Ja3IGdo995vG5S6vSxaVcZbWw7j9gTbPGagy8HL35+ilnkUiiSgK01aTjvCEZ0Ve6q4ZsFKfvbK1phwp1pNfHtyf3540WAKsu0cqfezfE/VSc6m6LFkFcKk70HuMPBUwJjZMPZGuU9EpHCbUuCsG2Hbf+DIVjkvPokodDn4+Yxh3H/FKMb0ObbyPNNuVsKtUHxVAo0tBxd1IWdsD1NbrV9vbj7MmD5pPPZhMWtKZaoxarQihGBYXiqzxuQztm8GN03qr9LjZwIWB4y6FrKHwK6lct3bbJfCDdI3vXSFnGS27T/g3gNDZ4LV2b33nYDdYuKKs3qTl25lzpNrWrQ1NjSF2e/2HNMCqf6uFYqTUH9I2iy7hsLIK7v88mekeEdbvxav2s+SeZOp8Qa59em1HK73YzEaCEZk69fI/DSmD81BQ3pIXzAsl3SbGVBWpmcUmgb5Y2WU/dwVUrjNdmmnGgnKVHptKZz9LXDvhtoy2VIW9UdPAtyeAPe+urWFcANEhOC6J1bz6OyxTB+aG+uqANTft0JxPI5shR1vgB6SxazdwBkp3rPG5PPcyjKKKz1Mf+hjwrqIvakFIzp9MmzMGNULm8VIqtXMjNF5FLrUuvYZT+mn4DkK6f1h2KVyctm6p6GpWha1rXocBl0MvcfB2z+FyXfCwAvB1P3z2ZcVHYkVpw10OfjVZcP50ZJN+ILyjeeel7diMxtoCulq1KhCcTwiYTnkaP8XMgsX9HTbrZxR4n2orokyt5ddRxo4f4iLgzU+fAkDk9NSTMwcnU+/TBtNoQiTB2ZzbmEWVpMyuFAQryYfeTVEArDt/2Tr2P4VYHbIiHzvR1Lk9TCgyVGjo64FR3a33no0ivYEwsye2A+X08p7P57Or18rYkVxNQLpUWDQ4DeXjVBpdIWiNUEvbHtNej+k95OujEq8Ox9PIMyLaw+gadDYFOazYneLFKLFqPHdaQX4Q4Ish4Wvj+pFbqoytVC0Iirga5+UqXFHDlRuh6Ya6YmuR5qFG1mN7quFtf+Ugt9rpHzdyKvjQ008VV3WJ95agPtlOXj0xnFc/MjymOmQLuA7z67jFzOHcc24vsx5ao1KoysUnirY+qIsUMvo3913A5xB4i2E4Gh9E1vK69lTEX9aMhs1zAYNX0hn0aoDPHbjWUwbnIPBoIaIKI5DdLDJuqfkdDF7NrwwOy7aBpP8fM87zZH3NbD1JfDXw4Zn4q+D+KAT6PI+8ahzYIM/TJbdjC8YwR/W0QX87zu7eeT9PYQiQqXRFWc2Nftgy0ty0mBq8vw/OCPEe2+Vh9+9uZ0VxW5AjuqM6AKXw8Kt5xVQ1Rjg3e1HqWgIUFbtY/pQJdyKEzDyainAVbtgweTmgpWEQjCTFUbcIAtaKndI17az58iCN2de/HUQHzE68uou/zaWFR2huNIT8+wXQnDV419wuN4PQCgi0IBbpgwgy25hwfISAO68cDCgUuqKM4Dy9bLLxO6SnSdJxBkh3iaDxsqSasxGjamDXEwb7GLroXp6Z1gJhnWun9CPuy4azLvbj6o3IsXJcebIyHnBZCm+UWxZUsj9dVD8AUz6L9jxmnxyX/MEDJ0Bo6+HDU/HX2d3dels8ESif+uJLZPPfvdcrp7/Bb6gnBEugN++sZ2nvyhlX3PBG8Dsif1UZbri9EXXYd/HUPoZpPUBY/cXnbamQ8Rb07SZwF8BI/CUEOLBVvutwCJgAlAN3CiEKOuIa58KA7IdPHT9GHYebmBQbirhiE7fjBQGuBzMHJVHhl3+YtQbkOIroxnh1ndkYdozl0J1sWwhO+cO2Lccit+HPe+Ca1iLIL27Sfybd3sC3PX8RnzBCJl2M3W+EAJ5u4nC/eRn+/jXilKqvUFyUq1MHpjd4hwqGlf0aMJB2PkWHC2SnSWG5CxYbrfDmqZpRmA+cCkwErhZ07SRrQ67DagVQgwG/gI81N7rfllmjMrDZjHS4A9xtMHPRcNzuWFCv5hwKxSnjKdKrlX73PF0mojAq7fI/d95By59GC76tTym/2TphW6yyT7wsA/MNvnhc8MzM+MuTZ4qePnWlq5Nniq5zt7JJKbRP7jnAt77yXTaKv2o9YWo9gZxWIxUNUrBd3sCsTX0+97YzqJVZZ1+vwpFhxPwwObn5XJXRvIKN3RM5H0uUCKE2AegadqLwFXAjoRjrgJ+1/z5q8DjmqZpoouN1cO6AAHfmjyA3hm2rry04nRix+ty3Tpn+LGFZ9HK8Unfk9sz+knnNasTCqdD8Xtyux6B3NFwZCNUl8DGZ2H8rfDEedKK9cBK+P4XLc8NnVrU1lYaPd1mjk3SMxogkuBHoQPZdgvFlR5m/OUzgFg0rgrcFD0OXw1sfgECDZDet7vv5qR0hHj3AQ4mfF0OTDreMUKIsKZp9UA24KaLsJmNXDQ8l7F9MrBZkvdpStEDSOz3jq5V37K07Zav1Dw453Yo+UgKtsEE9QdkSu7IxrjV6uePwep/gK9arq95KuLFcE01XVbUFhXwaBRd6wuR7bCgCxETcQ2ZSm8KRmhqXhuv9sqhJ0ZNo6pRpc4VPYzGo7Dp+fjo3x5A0g0m0TRtnqZp6zVNW19V1XGG7yajgUmF2Uq4FR3DuXe0LDJz5hw/KjZZYfgsGHMD9BoFI6+BEVcCmhRug0maPfiq5fFpvWXxm88thRvkQJQuLGpLTKG/95Pp3DF9YGzfnEn9mTAgo83XRYQg22E+Zh1cpdEVSUttGax/BjQDOFzdfTenTEdE3oeAfglf923e1tYx5ZqmmYB0ZOHaMQghFgILQY4E7YD7UyiSg7zR0tSl6GXIGSbXyzctlv7IMTT5ZkLCYrNmhKGzuvRWW6fQo+1hIFvFGvwh/vL+bp5duf+Y+rtqb4jvL97Ay9+fAqCq0hXJS8UOKHoV7JlgSZ5hQqdCR4j3OmCIpmmFSJG+Cfhmq2PeBG4BVgHXAx939Xq3QpEUpPaCibfJ3tGmdWBNjUfXsYQ0Cf9q8WK4Lm4pay20iQKelmKmwOU4buH8Pre3xTq4MnpRJB3lG2Hnm+DsBeae56bZ7rS5ECIM/AB4D9gJvCyE2K5p2gOapkXnpP0LyNY0rQS4B/hle6+rUPRYLHbZ720wSuG2Zcte0rak0JYJ2YPjxXBJxGVje5PjlN0aNrMxMVeA1WSg2huk2hsk22FhybzJama4InkoWwk734C0/B4p3NBBfd5CiGXAslbb7kv43A/c0BHXUihOCwwGuPRBWYUe8IIto3ke+G6532iJi7s9C752f5fbp56MZUVHqPLIqPqFOyax5UAddy7ZRDCsEwjHy9L9oQg13qAadqLofnQd9n0ihwel9QNjz/UpS7qCNYXijOLi38C0H8uitvSE0pFIEAacJ3vDq0vi1oxd1PN9KsydUsADV41iybzJ5KSmcMmoPJbcMQmrqeXbijcY4fK/fc6eikbVC67oPvSINEoq/Uyar/Rg4YYzxB5VoUhqUnvBObfJFLmIyKi7+D0o+VDuN9nktobDsPiabhtk0hato+fthxsIhHV6pVoJ6yLWQhaM6Fz598+xWYzU+kIMdDnwBMLdcMeKM5JIGHa9BYe3SPMVrefHrUq8FYpkYMuLMGwWpKTB4c2ykG3bfwAB4SZ495fw/m9kS1k3DTI5FRKr1J9bWcrfP94b2+cP6/ib0+n1TSEefnc3TqtJpc8VnUs4KOtFKndK06TTQLhBpc0Viu4nOmL039dC/6mQPw72fgyIuD1j2C+F2+KAm5Z0yyCTU2XulAJcTiu3TC1kUE7bk5hUBbqiSwj5oegVma1KP32EG5R4KxTdz8irZTRdtQv+MRWW/TReqDbqelr0fOs6bH0RvG3aJCQVLqeVf3xrAkbtWIN0DXj0hrGqAl3ReYSaYOtLUFPaLNyn16hnJd4KRXcTHTFqd0lXtejAk2sWwt4PkS1kzW884SZYsxBW/g0ajnTnXZ8SH+6sINKGpYMAvvPceoorGrv+phSnP0GvtDutL4f0Pt19N52CEm+FIlk5WiQtU+0uOZXMkSu3+2vh8CZY/7ScFZ6kuD0BXtsozRbbeqNxe4Jc9rfPWfhpCbouYq9RVeiKduFvgI2LwVslrYY7kEUlNtz+eAR/tMnAUztO8IJORIm3QtHdtB4xGo3At74k+7tveVuaSYz7Vnza0aENMrrYuAiObO3e+z8OUX/0HKcFHRjoclDosgOQapW1ssGIzp/e2c13n1vH3spG1UamaB9NdfL/hL9ODgVqJ4livajExn2bU7n0wyzcfo19dREu/3wAf9hg7pa/V1VtrlB0NycaMWpNhV7DIS1PzhkeMxt2vS3NXNYuhLO+KavSA43gGgI73kiKFjJoWXm+rOhIrDgt+vmCT0p4Ye0B/CGd5bur+GxPFbpAFbIpvhq+GjkrIByQlqftJCrWi/fZWDK9lsk5ASwGJ1V+I+ctyyasQwQDNqPeLX+vWjJbjE+cOFGsX7++u29Doeh81j7ZcsSop+rYEaMBD2x5ARqOyqK1hub5P2NvAmsa7F4GDeUw65GkEfCTsc/tYdZfV+APxR3Zbp1aQP8sO1ee3Vu5silOjcYKKdwA9uwTH3sKLCqxMTknwF1rMihuMJFp0Qnq4A0baFGDguC6fl4euXM2WgcVxGmatkEIMfFkx6nIW6FIBlqLbVsjRq1OOPtbcirZ8MtlejDcJIXcaIVIQM4iHjqz6+67naSlmLFbTPhDwdi2Z1eWAfDcyjL+c+dUQE0mU5yAuoOw6d9gSpE2w+0kGnEPSbMxf1IdN36aRW0wusKcKNzgNIS4rE9Thwn3l0GteSsUPQmLXUba+WfB+FvAYJbbIwHpwjbqGhmB++u79z5PgahVak3z8JJ0W8tYYn+Njwse/oQLHv4kNlt81ph8VdSmiFOzDzY+J/0POkC4AWb19TMkLUxxg4nZy7OoDSYKs/zcaQyTagzj0S38vigTtyfQIdf+MijxVih6GuYUGHMDuAa33B4JykK24g9g3b+g8Wj33N8pEi1oG5Lr5L2fTOejn17IoBw5U9lukeY03mAEbzCCBsyfMx6AG55YqYraFHB0u6wqT8mQtSEdhCtFsGR6LU5jhLqQgahga81T/0yazi8H7uOZc8sZ5PCz32tiWVHXt22qtLlC0RPx18NnD4MeAjQZdUcCcHCN3K8ZpUnF2NmykC0JSSxoi65tv/S9ySwrOsLgHAdzn15HuLmFTABXPf45KWbpjQ4ob/QzFSHg4FqZYUrNk+nyDqTSJ/jV+hQ8kXhsa9Ei/GJAMX8rH0R92EzYkcfEPn5esh9gWeOgblnKUeKtUPREdrwOnspm0Q42u0dpxGaC1+yFfpOloYvBBDP/lJQOU63f9FxOK7PG5HPzwtWEdYHVZIiNF20K6TQ1F7YNyLIze6KcwqaK2c4gdB1KPoL9KyCtLxjNHXbqiA5P79T4654sPBEpjUZ0rJqOT5j4d9VAXrqwnrVuK3MH+wFwWSPM7a2f6LSdhhJvhaInEi1mG3A+PHeZ7AuPoUF1Max/Skbj3ioINsJlj4LJ0i23+2VITKe/cMcklhUd5v43d7Y4xh+KAPF1c1XMdgYQDsrJYEeL5EjPqO9/B1Dk1vnlhlS2N8a9+HPMAX5RuJdBOQ5+vrUPJY2WZuFu6rDrtod2ibemaVnAS0ABUAbMFkLUtnFcBChq/vKAEOLK9lxXoVAgBdxT1WpjNPrWZNsYyNYZi1NWp4+5Pl7Ycyrtad1AYjod4NmV+485pqIxwNQHP8KoaTSF9BbFbCoKPw0JNELRq812p/07LIvkC+n8eYuJRftdhISBVFOYOb0rqPHDjN5NTOmfid1i5MXMWpaVpySNcEP7I+9fAh8JIR7UNO2XzV//oo3jmoQQZ7fzWgqFIpHWzmzQHIEnpM9BTlYy26TV6tqFch28+AM5yWzdU8caw0DSCPiC5SWUun0AZNrNhCI6noCMuoNhAYgWxWwqCj8N8VTCliWyhiPqMNgBfH4owq82ZXDAbwPg/Kx6rs05RKo5wsihaeSlp2FofkZwpYikEm5ov3hfBVzY/PlzwHLaFm+FQtHRtOXM9sxMqC6JH6MZIeSFVY/L1rKMAbDhWeg7KT7JbMFkeazPnXSzwp3NNqoDXQ5e/v4UAK6e/wXltfE3UgFc9rcVOK0man0hCl12Vcx2ulC9V9oEm1K+tGvaohIbs/r6caU0++b7NZaVp3DNAC9/2mDmpfIcdDRclhDf7XuYQms9ealWhvfOxGbquJR8Z9Fe8e4lhIjWyB8FjvfTTdE0bT0QBh4UQrzezusqFIpodBxNfSem0FMy5Hp3qAlMNmnmsv5fMPp66HsOHFwth5188sf4erndJR8CkmhWeOuKdLcngNkow6HUFBMefxgBhCKCWl8Ip9VIZUOAh9/dDcCdFw5mT0UjD7+7k+lDc1U03lPQdShfC7vfBYdL9nF/CVpbmwLc/FkmxQ0m/rIthdqwLHSbkVPLFdmHSTHpjMxPJz/dRvKVdbbNScVb07QPgbYc3n+d+IUQQmiadjyv1QFCiEOapg0EPtY0rUgIsfc415sHzAPo37//yW5PoTizSUxv73hdRt3RSLy6BF6aI9Pl0Si76GU5K3zopXJ/yN99936KJArusqIjlLp9DMl1smTeZLYerOO7z8UtlKMpdYAX1hxgwoBMvvXkGkK64MOdVcecT5GEhPyw+205cCetz1eqKJ/V18/ifTaKG0zM+CAbXRfUhoyAoDZsppclyB39D9PH3ECvVCsjeqf1iGg7kXZ5m2uathu4UAhxRNO0fGC5EGLYSV7zLLBUCPHqyc6vvM0Vii9J6yK08vXwyZ9g4EVQsU2Kt9Ahbww0VoK3Qq6HGy2ydzwq/EkUfbdm0aqyWDFb4hr3yRiQZec/d06N9ZQrkhCvG7Y2P2Cm9m5XYZrbrzHjg2yqAy29yGbk1HKl6xBWo2BkXhr5Ge2Itn3VkD0YRnXcUtOpepu312HtTeCW5s9vAd5o40YyNU2zNn/uAs4DumkCqkJxmnPuHS2Ft+9EuOIx2SqWOwLOuUP6oB8tksLtyIXz74WxN4MzT0bnO5J7VWvulAJcTmuspazQZSfTfuLoLC3FFPNJX7C8RLmzJSOVu+TDZ8gnI+52VpSHdIE3FP9aQ3B3wQGudh0kP9XMeYNd9G6PcHcz7V3zfhB4WdO024D9wGwATdMmAt8XQtwOjAD+qWmajnxYeFAIocRboegqMvrDhFvl1KXUfJj6Q/kmGaiXUbiIQO5wMN8ko/P8syASBmNy20DMnVKAJxDm1fXl1PpCZDss6ELEHNgSafCHKalo5FevbWOf29viHIpuJhyEkg+lO6AzF8z2dp9yVYXGrV9kE9BlfGrWdELCwCtHclk4WWNUjrXHinaUdkXeQohqIcTXhBBDhBCXCCFqmrevbxZuhBArhRBjhBBnNf/7r464cYVC8SVI6w0TviOF2myHaT+WY0R9bvjiL1BbGu8H3/2OLG6rP9Tdd31SnFYT+9zemD/6E9+ecNxjb3pyTUy4C112NTM8GfBUyr+1QxvkQ2Y7hTuiwxPbDHz7cyncBgTf71fOHwbvpK8twOFAChvrM3q8cINyWFMozhycuTIC37gIDqyBQINc6w56YdV8Kdz+Ohh9nRz0sPZJ6D8ZCs7r0MEPHUnravT739gW2zcgy87t5xXw27daJvosRo0n505Ua9/dia7Dkc2w621Zc/EV+rdbt4JtcWv8cE0qB5qk1/kAm5/b++wn3xpkRH4aFw9r4J0kM1ppD+0qWOtsVMGaQtEJ+Gpg9QLY8IwsEDKYQG/ui7Y4Yfq9kJIGegQ8FYAGAy+APhPkRLMk567nN1DocnLreQUAXPjIJ3j8kRbHuJwWfnzJUK6f0JcUc8+qMu7xNNVJ0Xbvkcs4pi//EBWfuR3m39NqeXWfxp93ZaOjYdF0bu9/hNG2anJTrYzMT4tNqetwurFgTYm3QnEm0lQLq/8Bnz8KkVZrxOn95bp4tEUnEpTpTbMNBn8deo1K+vVwkL7ns59Y1WKNuzVj+qTz35cOZ9LAbIyG0yGZmsQIIQsld70ti9G+pOlKIm6/FuvbNmk6YSFXgK0GnV8U7qWv1c+IvFT6ZNrp1F9rD642VygUPRFbJiBaCrexOQKqPwDL/xf8Dc3bLTKtabbJSvRVf4fDW44V/SRjWdGRFmvc7/9kOoWu+JqqBhQdqufWZ9dx5783UFLZ2E13egbgq4EtL8G2V6WBUDuEGyDDojOrVz0gYsKdYojwwMCdjMqIcN5gF/2yOlm4uxkVeSsUZyKeqmOtVKNoBlmFnpIB59x+7Hpk0CsjjpQ0GHgx9BrZoaMZO5JFq8rwBMLMntgv5tD28vqDRHSBPxjhjS2HY1ar6TYz14zrzZ0XDiY3LfmXB3oEkbAsRiv5UE4Bc+S2uwVsQ6XgNxud7PS0dF1LNYZZNPkQZ+XZuk60Vdq8bZR4KxSdxNon5WCS7MHyDagpYRjg4K9LUa8tBYMZxt4o+8VbE/TJ11qcMPhrSS3ibRGO6GzYX8Oi1QdYWeKOtZj1yUjh1qmFCATXju+Ly2ll0aoyJg/MZvW+auZOKVDTy06F+kOwayk0HoXUPJnBaQd1AcGfNll4tTwDHQ0DAh0NhyGM0QANYRND0sIsmV4bK2LrdLpRvJN/4UqhUHQ8reeBRzGYoXA6DJkB216Rvbeb/w11+6VzW+IMZYtdfgR9sOMNKPkIBl0MeaN6hIibjAYmDXQxpFcaH++s4PXNh1m/v4ZDdX7+uEzOD//X56V8a1J//vxBMRajgWBExxMI89rGQ2p62fHwN0DpZ9LdLyVNtoC1g7AOi/do/H13BjUhCxqCYQ4fu7128q1+/nlOJXnpNuaskGvgyTa6s7NQkbdCcaYSHSlatat5pKiQkYTDBVPvlhH1gVWw7T+yPzxroJxMlpLe9vlCPlm9bk2FwZdA7sgeUdgW5UC1j7e2HObDnRUUHaonrMv3RqMGupDTy4yaRppNTi+LTjqLpuPP+Eg8EoLDm2SKHCEd+7RTK6s63gSw3ilN/GlrKvt8cmxngc3P3D6HyTd72NTUi28Pg36pWovXdKlwq8hboVB0OW2OFL0Uqoth/0oYOgMGTJUGL+ufgZp9sOLPUsCzBx17PrNdRllBL2x/DfZ9IkU8Z3jLiD1J6Z9t53sXDGT60ByWFR3h410V7K7wEEmIbyIJDm6ieWb6guUlvLL+IKVuX2x9fVnREWaNyT8zBF3X5TJL8XuyMM3Z65Tav6KCvaw8hfs2p/JMiY0bCpqYXeDn6o8yKG8yA9JfIMMU5obelYy3u3FYjYzMz+JSZ8uF7WScud2ZqMhboTiTaT3IxFMlbVT1sEyh2zLk9kCjnANes1dGU8Mvl8NOTlR8FPTIyMSeLUXcNQwMPaPBJRjWWb3PzbzFG/CH9GP2a8hI3GEx4g3Ge8j7Z9nQ0Nhf4yPHaaHKE+SBq0advgJedwCKP4D6g2DLOmUzn8Q+7fmT6vjeqgxKPTKWjK5lA1g0nSvyqrkgvRK7UWdQrpN+mXZMyVJGrgrW2kaJt0LRTTRWSBMXk02uW4I0bdn9Nuz9WH6dNxbOulm2kJ2IQKOcEuXIkcVw2YOTXsTdnkBsYllqiolGf/ikr4kKeiLR0aWnnZtb41HY9ylU7ZQ2u7bMk74kMTXu9mvMXp7JPo+JLItOMCLwRBKzM4KLs+uYlX2UNFOYQpeD/ll2rKYk+7tRfd4KhSKpSO0F4+fKdexAc/+zwQgjroSJt4EpBY5ulWn0hiMnPpc1FTIGyLahzS/Amn9A5U75dZISnVg2JNfJf10olwiMzVkG43GCvrbCoGvG94mtiZ8Wk8waj8LWV2DNP6G2TBr6HEe4F5XYcPu12Of3bU7l0g+zYtuiyxE1QUMr4QaHIcKs7COMyjUzbYiLIbnO5BPubkateSsUirZJ6w3jviW90DVNFrCBnAV+/s9gw9PQcFgONjnrJug9/sTnS0mTH4FGObPZlgGFF8pRpab2tRF1NK09051WE8PyUvnB8xup8gRP+TwTBmS2iOKj547OJI9G5Elf8NZwGMq+gMod8vRaQQEAACAASURBVMEtvW+sGK2tYrP7N6Xy9qEUFu+zsWR6LZNzAlgMTqr8RqYty0ZAbOJXInZNtn01Rkw8cWQILw6tw25O3uxwd6LS5gqF4sTU7pcCbssES4IxRiQoRfhQ8//RwgtkZH6qxWlBryxwMqfAgPMgf2zSDkCJ4vYEuODhT2Lr3A6LFCBv8Nh1cZCp9HSbibqmMENynVwzvg8b99fy4c5KhuQ6mT9nPG9sPsTL6w5S5Qly78xh3Hnh4Ni1ulXQhZDRddkKqCmVou1wgWY4ptis0BmOFZtFbUtzUiJU+Y1kW3UiuqAuZETmJ1qmLhLtTfvZQ/xzSh13r5PneODsxuQuQlPV5gqFImnJHADj5sDGf8sIPDq20WiBs+fI/dtfg9JPob5cjh61Ok9+XotDfoQDch197yfQ+2xpCJOa17nfUztIMceL1IQAXxsFbVEEUNcUxmLUGJmfysPv7gbAbNQorvTwjb981uL4l9YeBGD2xH7HROtdRiQM1cUs+mgDs9LLcDlTIL0f7oCBZXul89x9m1NZvM/G/El1FDptlHpMPLwtlfm7HHjDBganhrlnZAP3rMugOpAYYbcU7lyzn4fHHmJ9YzZvH3ZS6jGzvjqFJdNrz5h+7a+KirwVCsWp4S6Gzc+Do9ex08Vq9slq9ECDjNDPuR3S+ny58+sR8FVBOChT9gOmQvaQpEmpJ6a/sx3ynqq9MoVuNmqEIoI0mwmPP4x+nLfVtoraotv7ZNgor5NiFa1iP17BW6ek3QMeqNgOZStYtNfOfbsHxBzLgFhEfUlegL0eI6UeE9lWnbAO9aG4QFs0HacpQk2otVGPjLqdhjCaARrDJgY4QvznorpYEVuPE+yeWrCmadoNmqZt1zRN1zTtuBfTNG2mpmm7NU0r0TTtl+25pkKh6CZcQ2DMbDkmNBxouS9rIJz/U9nn3VQLX/wVDm/+cuc3GKWxR7RXvOg/siCu+ENZKNXNJBaxvfeT6bz3k+kMyZUZhlBEMCTXycc/vZAl8ybHitusJi1WaGXQ2hZukNsrGuOi5Q1GyHZY2hTuu57fwH1vbOfmhatxewLsqWjk0sc+4743tn/5ojgh5Hr2zqXyd1b8HpgdzBpiZ0hamOIGEzM+yGbGB9mxVPiHR63oApymCNUBQwvhBggKAzUhM4bm79ZmiGAhAmiYNZ1F51fz8cxahqSF2e81s6xcPgieaX3a7aVdkbemaSMAHfgn8DMhxDFhsqZpRmAP8HWgHFgH3CyE2HGy86vIW6FIQo4UQdGrkN77WL/qSAiKXobydfLrId+AoTNP2WnrGCJB8LllKteZC/2mgGvwqaXlO4HjRbwQL25btKqM+97YTk6qlXfuPh9dF9y4cDWlbi/j+qez6UD9KV3LajJw10WDGNc/k0E5ThxWE//ZUM4DS3fErFoz7WYamsJEhCDHaeGdH08/YVtaOKIT1gWhgA9RVYx2YBV4jhIxmAlbs9E1IxEh0HUZCc9d3Yf6sKxhsBoiDHd42e5xxtaoWyL4/+3deXTc5Xno8e8zM5rRLtmSJUu2ZHnfwTaKAWO2mLA4KYQAKdm4JDlxktv0tElaGg5/3Jz03pO2NLlt05xLaJI2zU0aCA0XJ7hhSULMZoIN2NjYxrvlfdW+zPbeP56RLQvJGnlGM/MbPZ9zBs+Mfszv5edBz+993+d9Xr84Ys5HiT9CV6yAKYW9/OPSY0wsCfHHL9Zwss9/bh7bkz3twby+zltEXmD44H018HXn3C2J1w8COOe+OdLnWvA2JkcdegO2P6VD44PrmDsH+16Ad9YCTteDL/lEUlW3LqqvA3padd69erZmt0+YlvrnjoGhgvzjG1vOVWIbLZ9AaShAcSjAma4w4eiF8+wCfP76GdRXFOH3CT6fjgZEYnHC0cQjEqWw9zgTO3ZR2b2X3piPM1TQRildMT/t0QCt0QJaI/rnyXCQ9ljyNeprCnr502kHefTwNFp6Qtxa38P/XNb5npKnng7Wg+V5wtoUoGXA60PAlRk4rzFmrExdBi4CO9bpsiHfgF8lIlp9rXQyvPEjXQ/+8imdBy+eeOnnDJXpw8U1Me7kTr1xqFkIkxdBRUPOzI8PnnvuX27WH7hF9B5nKKVBP3Mnl7H1SDt90TiFAR+90TjtvVHahykW44BHfr/3gvck8Q9f4oiY608Wm5J4JMeHQxI96spAhPklHbzadv7vsdAXpyoU43BPIbHSOp5c1T5kkM7LYfFYRKsRZsGIwVtEngeGSv18yDn3VLobJCJrgDUAjY2p7UZjjBlDDVdCNAK7n9ViHYOXiNXMh5Vfhte/Dx1H4KVvQ/NndH48FeLTkqvFVfqL8+QOOLYZxK+boUxepDcUI1V+y7D7rm6isy/KD17cx+mu8AWlVatKgty7XGui7zvVzYpZ1Xz22hk8s+0YTVXF9EZ0N7P2nggnO/rYsO/Mez7fJxAK+IkkhsYdgAM9gwbuAolT4HMUiKPYH6M0EKfEH6MkEKMyEGNCMMrEggi7OwtZd7KKumAvX23Sm4JvHZjB0b5CdnZrxb3KYBy/wOk+HyF/nAcWdnJ/IjjnXZAeyvFt8PbPdQ+AxXdn/PQjBm/n3E0pnuMw0DDg9dTEe8Od71HgUdBh8xTPbYwZS03XQKQHDrykiWaD57ZLa+GaL2sP/NROePW7sOhumHZ1es7vC+hcOGi2+undcOxt7dpWNuqQfWWj9vgvVoc9Q0pDAU53hc9lkT++sYUnNh5i76kuassL+fkXVlyQNb56cR3xuKO9N0J7T5SDZzr5k59qIqAAoQIfvRENojEHH5jcxd2lWymKtiI+oa+gknigGBEffnH0lwR3QDwxtx1LPHQXNc0IX1jkqCoIc3N9D9MrSygNBVg5q52/3gxPHy4cMgu9dLwUU+k+Ddt+ocEb4PAmHUbJ8PcrE3PeATRhbRUatF8HPu6c2zbS59qctzEeEI9rzfNDmxIBfIhfYvGYzpHvS6xrblwBiz5y4XB7Orm4zpH3l3YNlmivvGoWlNddWGwmw1JZ5jUwGe4nn1pEefQUn/jZPva0+7ip6hTfmLWLDkrojIfoi8WJxuJEYhqcHbqvqYgggN8nFBb4KSzwESrwUxL0UxwMUBTU9wqGqT8/3Paded/bjkVhz290y9N4RHMtpl8Pl90Li+9K22kykrAmIncC3wEmAa3AW865W0SkHvi+c2514rjVwD8AfuCHzrn/lcznW/A2xiPiMS3UcuIdnXseTstrOtQYj8KE6VrQpX/jk7EU7YPeVl1DDlBaDVVzYEKT9twz0YZU9bRC5wn+/ZW9rK7cT3X0BCCcihaz7kwt982ODPuvOs53DrM//uBBp/fAlseg64S+rr8CFtyu32MvZ5uPFQvexnhINFEutfWAFlkZTutB2PgD6G2DwgoN4BOaMtZMnEtsuNKuvSlB67ZPaILKJt0etXhiVnvnWjr2tG7RenYfnD0IkU40Ay2giXsFxTkxFZDXwt2wfS20bNDXJTWw+B5d7QDeXyo2Vix4G+MxkR6tg97Ten4ueih9Hbrl6Jm9mmi24MPQtDJ7wSjap/uPR7o5VwctUKS7q5VP0e1MgyUaNIOlmgyXalujfRqkw51a3azzuBZMaT8C0d7zxxUU6TlzcElcXju6GbY+od9V8eue9LM+AP4BUz15vlTMGDNeFBTpDmMb/1UrrQ23z3OoDK76E+3V7Ps9bPtPOLtX5w+zEaQCocR5q86/F4voxilth/W5yIVrvAKFicBarM99Ac24F7/+6eJ6rHM6rRDt1qAc6dWbhFj/MHfiZsHn1950YQX4qzN8Acw5fZ0atI8mKgROnKGVBXOs3r4Fb2NMehVWaFGWjT/QXstwO4X5/LDwTh2u3vwfcORN7XVe8enc+EXpLwB/hf73DNYfkONR6OuC3vZEUHfngzac752L73xQ9wc16I9Vsp65dEfe0sAd7tS/p/l/pDveXWqFwDFk3x5jTPqVTtIAvulfNUhdbM11/VIoq9NjO4/revBFd0PD8sy1d7REdPjUb79C80K4S4P2kTf1ddVsuOyPdQvUHJV7txPGmPxQ2aC/ADtPaI3yiymbDCu/AlOa9djNP9UdzAZvgGJMup3YDr//Ow3c/qDeOF71xZwO3GA9b2PMWJo0F+athh2/GroK20CBkPbWq2frxieHXtfM9WX3Xzx73ZhLEe2D7b/UAkOgSxeXfCLng3Y/C97GmLHVsBx62uDAy8MXceknomVXKxth04+g85gOo8+/I7vZ6Ca/tLXAGz/Wddvih7m3wcz35+Tc9nAseBtjxt6sVbqu+sQ7Wnd8JGV1cO1XYOsvdI3ttv/UGuaXfyxr24GaPODisPcF2PE0uJhunrP0k8l9J3OMd24zjDHe5fNp5m7FVOg4lty/4w/qsrNl92vC24ltsP7vdDcxY0arpxVee0SXJ7oYNF2rN4geDNxgwdsYkymBoFanCpVB16nk/736JXDdA7retq8dXvs/ujHESElwxvQ7vhXWPwyn3tWCN+/7HCy6S28QPcqCtzEmc0KlsOTj+ry3Lfl/r2gCXP0lmLta5yX3rYcXv637ehsznFhEp15e/z5EumDSPL0RrF2Y7ZalzIK3MSaziifqPGO4Sx/JEh/Mvhmu+XOtMd15DF7634ldnmJj117jTZ3H4eV/gP3r9bsz/3ZYvsYbm9AkwYK3MSbzyut0Prv71OjXclc2wnV/AdNW6tzljl/BK/+kv6yNcU53r3vxW9B+GIqr9IbPY9nkI8mf/xJjjLdUzYQFd0LHUd3dazT8QVh8Nyz/vJYvbT2gc5p7fqcZxWZ8ivTCmz/WcruxMNQvg2v/Um/48owFb2NM9tRfrkPh7S2XNvRdMx+u/ytdSx6Pwvan4JXvJJ/RbvJH60F48WE48kZipcLHYOmnoKAw2y0bE7bO2xiTXdNW6Nz3wVegomH0Q5sFxXD5x2Hy5bDlMd3/+sWHYfYtOlRqG4DkNxfXvId3f63Py6fAsvugtDbbLRtTKfW8ReQeEdkmInERGXb/URHZLyJvi8hbImIbdBtjzhPRvZLrLtftNy9V7UK44WvQcJX24neu03nP1oPpa6vJLd2ndaRl5zoN3NOvg2u+nPeBG1LveW8FPgJ8L4ljb3TOjWJxpzFm3PD5YN6HINytPefyKZf2OQXFmgg3ZZn2wjuOakZ600qY+8G8HUIdd5yDQ3/Q9f7RPghV6BLESXOz3bKMSann7Zzb7pyzckfGmNT5C2DRR7RkZapz1tVzdC58xvu1Z7//RXjhm3B08/m9to039ZyFPzyqSWnRPp0uuf6BcRW4IXMJaw54VkQ2iciaDJ3TGOM1BUXacy4sh66TqX2WPwgLbodrvwqV06CvTfcMf/1fRlfhzeQG5+DAK/D7v4GT2xPflY/DFfdDsCTbrcu4EYfNReR5YPIQP3rIOfdUkudZ6Zw7LCI1wHMissM5t36Y860B1gA0NuZfer8xZgShUt2acdO/QfcZLeqSivIpcM2f6S/+Hb/SzVFO7dRe+aybdCtSk9s6j+s2sad36evaRVpqt7Aiu+3KInFpGEISkReAv3DOjZiMJiJfBzqdc38/0rHNzc1u40bLbzNmXOo6BRt/qD3odP2S7m2HHb/UvcIBCithwR1Qt8S2G81F0T7Y9azuBOZi2sNedBfULc2Nv6/u01A1CxZ+OG0fKSKbnHPDJoD3G/NhcxEpEZGy/ufAzWiimzHGDK+kWtfpRrqhrzM9n1lYrr36FX8G5VOhtxXe+JFWaDu7Lz3nMKlzTvMTXvgm7PmNBu6Gq+CGB7XwSi4E7ixLdanYnSJyCLgaeFpEnkm8Xy8i6xKH1QIvichm4A/A0865X6dyXmPMOFFep8G2t1Uz0dNl4nTdDnLxPbrL1Nl98PI/6lC9zYdn19l98Op3ND+ht1W37LzmzzUXImh7ufdLy7D5WLFhc2MMAKd2w1v/F0pq07/cK9Krvbu9L0A8AuKHxith1s1QVJnec5nhdR6HHU/DsS36OlgCc27TIj65WpM8i8PmVnrIGJP7qmfB4o/Clse1N57OfZgLCmHeB2HaNVrs49DrmtzW8ho0rtCktnGcGDXmOo5phbQjb2ihFX8QZtygCYW2Ln9YFryNMd5Qu0B7ONue1Axyf0F6P7+oUgt9zHy/lto8+pauDz+4ARqu1IBSUp3ec45nbS2w6zk49jbgtHfdeDXMudVulpJgwdsY4x31SyDWBzvWacKZfwx+hZVN1rXD7Uc0iB/bAgdeggMv6/lnrtJ5WDN68Rgc36o3Rad363u+gN4czXy/bt9pkmLB2xjjLQ1X6haiu56Bikbw+cfmPOX10PwZaD8Ke38LhzfBkTf1MXGGDrNPvnxsbiDyTW87HHoN9r+sSWigw+PTVsCMG62nfQnsW2eM8Z6mayAW0aA6lgEczme8z12tSW0HN8CZvfoI/kJvJhquHBebYYxKLAzHtmoOwckdaKFNoGSS1pqfulyrpJlLYsHbGONNM67XAH7gJahsHPuM5KIJsPBOmHub9sIPvKxD63t+q4+KBpjarOuQQ2Vj25ZcFQvDiR061XB8K0R79X3xQ80C7WlPmpu72eMeYsHbGONNIjBrFbio9oYzEcABAoU6ZN64AloPwMFXtaBIW4s+3nkKqmZD3WUw+bL8D+Q9rdqzPvEOnNiuy+36VTRCw/ugfqmt0U4zC97GGO/y+WD2LRCP6/BsZUPmenUiMKFJH4vuguPb4PBGDWCndurj7Sd0fnzyIqiep8lwXq8OFu6Es/vh1C44uRM6B+0AV9moNy2TL4PSmqw0cTyw4G2M8TafT4eyicOhTZkN4P38Qe1d1i+FcJcOGR/drMHtzB598BSEynXYuGqWBrnS2tweQo6FdWqg/TCcPaDVzwbv9uYP6Tr8SfOgdrEVtskQC97GGO/z+TShLB7X9dkVU7MXFIMl55PYIj06nHxyhwbyvnYdIejfGMUf0rZWTNVeakmNBvRQWWZ76LGwloXtPJ54nNCg3Xmcc4lm/XwFeuNRNROq58KEabrcy2SUXXFjTH7w+WH+h4A4HN2S3QDer6AIplyhD+eg46gOp5/dD60HoefsgJ75AP6QJsgVVeqfhRU6Zxws0UdBsfb2/UEtVnNB8HR6rlj4/CPapyMC4U7d5CXcqXPVPWe0DeFhNn4Rn95MlCduMCbO0AI5Y5ndb5JiwdsYkz98fph/O+CDo29qBni2A3g/EV07Xl5//r2+Dg3iHUe1t9t1Qnu7kR6dSx48nzxmbfNBcbX2/ktr9c+yOp2jT2cpWpM2FryNMfnlXA+c3Avgg4XKoHahPvo5B9GeRM/4rD562xI950TvOdKty+RiEc3ujiUyvPuH2kUG9MyDEAhBQQmEShI9+FLdy7xoQqJnX56718gMyYK3MSb/9AdwAY5keQ58tER0WLyg+MJeujEDeOTbbIwxo+Tzw7wPaQZ4a4vuWGVMnrDgbYzJX/0BvGG5zi1bADd5IqXgLSIPi8gOEdkiIk+KyJAL/ETkVhHZKSK7ReRrqZzTGGNGxefTbSYbr9YAHo9lu0XGpCzVnvdzwCLn3GXAu8CDgw8QET/wXeA2YAHwMRFZkOJ5jTEmeT4fzLkFpq2ENgvgxvtSCt7OuWedc9HEyw3AUJvcLgd2O+f2OufCwM+AO1I5rzHGjJoIzL5Jt6BsO6jbihrjUemc8/4M8F9DvD8FaBnw+lDiPWOMySwRmHkjzL4V2g9pARNjPGjEpWIi8jwweYgfPeSceypxzENAFPhJqg0SkTXAGoDGxsZUP84YY96raYWufd6+VouRBELZbpExozJi8HbO3XSxn4vI/cCHgFXOOTfEIYeBhgGvpybeG+58jwKPAjQ3Nw/1ecYYk7qpV2jQ3voElEzSddXGeESq2ea3Ag8Atzvnuoc57HVgtohMF5EgcC+wNpXzGmNMWkxeBEs+Ad1ntFSpMR6R6pz3PwNlwHMi8paIPAIgIvUisg4gkdD2JeAZYDvwuHNuW4rnNcaY9KieDc2f1pKj3Wey3RpjkpJSeVTn3Kxh3j8CrB7weh2wLpVzGWPMmKmYCs2fhbd+ohuElNZku0XGXJRVWDPGGIDSSdoDL6zQvayNyWEWvI0xpl9hBSy7T/ettnKqJodZ8DbGmIGCxbD4Hpj6Pmg9cH67TWNyiAVvY4wZzB+AubfB3A9C+2FNZjMmh1jwNsaYoYhA45W6lKynFXrOZrtFxpxjwdsYYy5m0hxY/jnwF0D7URiyFpUxmWXB2xhjRlJao0vJJk63bUVNTrDgbYwxyQgWw2UfhRnXQ1sLhLuy3SIzjlnwNsaYZPn8uivZ0k9BuFMLuhiTBRa8jTFmtKpnwZVfgLJaG0Y3WWHB2xhjLkVRJSz5JEy/DtoO2cYm40m4S//Oo71Q0TDy8WMgpdrmxhgzrvkDOow+cTpse1LXhJfVgVi/KO+4OHSfhkgvFE/QOgCT5kGoNCvNseBtjDGpmtAEyz8Pu56FI29Aaa3tD54vYmHoOqnBu2aBVt6raABfdm/QLHgbY0w6BIth4Yeheg5sXwu9bVA6WYu9GO8Jd+kWsYEgTFsJ9ZdD0YRst+ocC97GGJNOtQt0i9F3fw3Ht+kaceuFe0dvmz6KKmH+7VA7HwKhbLfqPSx4G2NMuhWW6+YmkxfD9l9Cb2uiF25z4Tmr+wyEO6CsHuauhqqZujQwR6UUvEXkYeCPgDCwB/i0c651iOP2Ax1ADIg655pTOa8xxuQ8EaiZr73wPb+Fw29qb66wItstM/2cSyShdUPlNFj0Ec1f8MBUR6o97+eAB51zURH5W+BB4K+GOfZG59ypFM9njDHeEiqDBXdA3RLY8StdF15am5NDseOGc9B9CiI9UD0bmq6Fyuws+bpUKQVv59yzA15uAO5OrTnGGJOnJkzTjPQjb8Cu54G4BnGfzV5mjItD12kN2jVzdY1+eX22W3VJ0vmt+Qzw2DA/c8CzIuKA7znnHk3jeY0xxhv8AWhYrsPpBzZAy6savEtrbT58LPX3tMM9UDMvEbTrst2qlIwYvEXkeWDyED96yDn3VOKYh4Ao8JNhPmalc+6wiNQAz4nIDufc+mHOtwZYA9DY2JjEf4IxxnhMqAzmfACmXgH71sPRt8Af0sx0C+LpM3BOu3o2zLjBsz3twcSluDetiNwPfB5Y5ZzrTuL4rwOdzrm/H+nY5uZmt3HjxpTaZ4wxOa/jOBx8FY5u0d55SU1OZzrnPOeg5wz0dWrW+MwbNXHQA0RkUzJJ3almm98KPABcP1zgFpESwOec60g8vxn4RirnNcaYvFJWqwVemq6Fltfg8CZAoKTaEttGq/u0FliZOF23cK3MzxHcVOe8/xkIoUPhABucc18QkXrg+8651UAt8GTi5wHgp865X6d4XmOMyT8lVTBvNTRdA0c3w8EN0NmrS8xCZdluXe7q72mHu3TJ1+J7NEEwj6U8bD6WbNjcGDOuRcNwcgfsf0kTrnwF2hu3DHXl4lp3PNoHVbOgaaX2tD2wTns4GRk2N8YYM4YCQai7TCu1tR+BI29qjzwe0ypuoXJPB6pLFoto0I5HoXYRTFvh+ezx0bLgbYwxuU4EKqboY9YqOLVL58VbW0CAwglZ25oyo/o6tO64LwANV0L9Up1qGIcseBtjjJcUFGlvvO4y6GmF07vh0OvnA3mwTHvl+bLkLNqnSWjxmE4ZLLhDd24rKMx2y7LKgrcxxnhVUSVMbdZH9xk4exCOvw1n92sSlz+gtdS9tqtZtA96zkI8AoEiHRavWZBYBz8OpwmGYMHbGGPyQfFEfUxZAuFu6DgKZ/Zpwltbi9a59PkhWArBktxaRx6P6ZB4uFNfB0v0hqR6NlQ05FZbc4QFb2OMyTfBYi1OUjUTZt8Eve3QeQI6jsCZvdB2GIhrQAftmRcUa4LcWA+3x2MQ7dEbjGifns/n0yVe1St1V6+Sauthj8CCtzHG5LvCcn1Uz9K63vGYDkv3tGrWdvth6DyuvfVzy4dFl2L5fLpEzV8A4tdesPj1/cHiMXAx/TMe02HvWCTxcaKf7fNpBbmJs6ByqtZ1L66y3vUoWfA2xpjxxufX3m1JtQb0fvE4RLq0rGi4S2uCh7s0yIc7tccc7YNYWDf5cGiSnIg+8QcTPfhCTawLlUNRIhO+oFgLzYTKhw78ZlQseBtjjFE+XyLAWjW3XGe3P8YYY4zHWPA2xhhjPMaCtzHGGOMxFryNMcYYj7HgbYwxxniMBW9jjDHGYyx4G2OMMR5jwdsYY4zxGHHnSuHlHhE5CRxI40dWA6fS+HnjkV3D1Nk1TJ1dw/Sw65i6dF/Dac65SSMdlNPBO91EZKNzrjnb7fAyu4aps2uYOruG6WHXMXXZuoY2bG6MMcZ4jAVvY4wxxmPGW/B+NNsNyAN2DVNn1zB1dg3Tw65j6rJyDcfVnLcxxhiTD8Zbz9sYY4zxvLwM3iJyq4jsFJHdIvK1IX4eEpHHEj9/TUSaMt/K3JbENfyKiLwjIltE5DciMi0b7cxlI13DAcfdJSJORCzrd5BkrqGIfDTxXdwmIj/NdBtzXRL/LzeKyO9E5M3E/8+rs9HOXCYiPxSREyKydZifi4j8U+IabxGRZWPeKOdcXj0AP7AHmAEEgc3AgkHH/HfgkcTze4HHst3uXHokeQ1vBIoTz79o13D01zBxXBmwHtgANGe73bn0SPJ7OBt4E5iQeF2T7Xbn0iPJa/go8MXE8wXA/my3O9cewHXAMmDrMD9fDfwXIMBVwGtj3aZ87HkvB3Y75/Y658LAz4A7Bh1zB/CjxPMngFUiIhlsY64b8Ro6537nnOtOvNwATM1wG3NdMt9DgL8G/hbozWTjPCKZa/g54LvOubMAzrkTGW5jrkvmGjqgPPG8AjiSwfZ5gnNuPXDmIofcAfy7UxuAShGpG8s25WPwngK0DHh9KPHekMc456JAG1CVkdZ5QzLXBQaxmgAAAixJREFUcKDPoned5rwRr2FiaK3BOfd0JhvmIcl8D+cAc0TkZRHZICK3Zqx13pDMNfw68EkROQSsA/40M03LK6P9nZmywFh+uMl/IvJJoBm4Pttt8RIR8QHfBu7PclO8LoAOnd+Ajv6sF5HFzrnWrLbKWz4G/Jtz7lsicjXwYxFZ5JyLZ7thZnj52PM+DDQMeD018d6Qx4hIAB0qOp2R1nlDMtcQEbkJeAi43TnXl6G2ecVI17AMWAS8ICL70XmytZa0doFkvoeHgLXOuYhzbh/wLhrMjUrmGn4WeBzAOfcqUIjW6zbJS+p3ZjrlY/B+HZgtItNFJIgmpK0ddMxa4L8lnt8N/NYlsg4MkMQ1FJGlwPfQwG3zjO910WvonGtzzlU755qcc01o3sDtzrmN2WluTkrm/+X/h/a6EZFqdBh9byYbmeOSuYYHgVUAIjIfDd4nM9pK71sL3JfIOr8KaHPOHR3LE+bdsLlzLioiXwKeQTMtf+ic2yYi3wA2OufWAj9Ah4Z2o0kI92avxbknyWv4MFAK/DyR63fQOXd71hqdY5K8huYikryGzwA3i8g7QAz4S+ecjaIlJHkNvwr8i4h8GU1eu986MxcSkf9AbxKrE7kB/wMoAHDOPYLmCqwGdgPdwKfHvE32d2SMMcZ4Sz4OmxtjjDF5zYK3McYY4zEWvI0xxhiPseBtjDHGeIwFb2OMMcZjLHgbY4wxHmPB2xhjjPEYC97GGGOMx/x/PmNNKBsGDacAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def plot_gp(x, mu, var, color, label):\n", " plt.plot(x, mu, color=color, lw=2, label=label)\n", " plt.fill_between(x[:, 0], \n", " (mu - 2*np.sqrt(var))[:, 0], \n", " (mu + 2*np.sqrt(var))[:, 0], \n", " color=color, alpha=0.4)\n", "\n", "def plot(m):\n", " plt.figure(figsize=(8, 4))\n", " xtest = np.linspace(0, 1, 100)[:,None]\n", " line, = plt.plot(X1, Y1, 'x', mew=2)\n", " mu, var = m.predict_f(np.hstack((xtest, np.zeros_like(xtest))))\n", " plot_gp(xtest, mu, var, line.get_color(), 'Y1')\n", "\n", " line, = plt.plot(X2, Y2, 'x', mew=2)\n", " mu, var = m.predict_f(np.hstack((xtest, np.ones_like(xtest))))\n", " plot_gp(xtest, mu, var, line.get_color(), 'Y2')\n", " \n", " plt.legend()\n", "\n", "plot(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the plots, we can see:\n", "\n", " - The first function (blue) has low posterior variance everywhere because there are so many observations, and the noise variance is small. \n", " - The second function (orange) has higher posterior variance near the data, because the data are more noisy, and very high posterior variance where there are no observations (x > 0.5). \n", " - The model has done a reasonable job of estimating the noise variance and lengthscales.\n", " - The model recognises the correlation between the two functions and is able to suggest (with uncertainty) that because x > 0.5 the orange curve should follow the blue curve (which we know to be the truth from the data-generating procedure).\n", " \n", "The covariance matrix between outputs is as follows:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "B = [[3.10552457 2.9025773 ]\n", " [2.9025773 4.83607491]]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAD+1JREFUeJzt3X+s3XV9x/Hna4W2UycUukGDyo9IVAwK2hQRozgQkD8KiWSWbKMsEKaTLdG4iCFBg1uG7g8WM502iKLZgMmm1g3G+CFxCRatG1DBAaUsg4riKGIYWGx974/z7fL1em977z0fzrnn5vlIbs73fD7fz7nvbwqvfM/3nO99p6qQpFZ+bdwFSFpcDBVJTRkqkpoyVCQ1ZahIaspQkdTUUKGS5KAktyR5qHtcMcN+u5Pc3f1s7I0fmeSuJFuTXJ9k6TD1SBq/Yc9ULgFuq6qjgdu659N5rqqO637W9sY/DlxZVa8EngIuGLIeSWOWYb78luQB4OSqejzJKuCOqnrVNPs9U1UvmTIW4MfAoVW1K8mJwEer6vR5FyRp7PYbcv0hVfV4t/1D4JAZ9lueZDOwC7iiqr4KHAz8pKp2dfs8Bhw20y9KchFwEcCS7P/GFy89aMjSNUq18/lxl6A5+Bn/y/O1M/NZu89QSXIrcOg0U5f2n1RVJZnptOfwqtqe5Cjg9iRbgKfnUmhVbQA2AByw/NA68RXnzWW5xmz31kfGXYLm4K66bd5r9xkqVXXqTHNJfpRkVe/tzxMzvMb27nFbkjuA44F/AA5Msl93tvIyYPs8jkHSAjLshdqNwPpuez3wtak7JFmRZFm3vRI4Cbi/BhdzvgGcs7f1kibLsKFyBfCOJA8Bp3bPSbI6yVXdPq8BNie5h0GIXFFV93dzHwI+kGQrg2ssnxuyHkljNtSF2qp6EjhlmvHNwIXd9p3AsTOs3wasGaYGSQuL36iV1JShIqkpQ0VSU4aKpKYMFUlNGSqSmjJUJDVlqEhqylCR1JShIqkpQ0VSU4aKpKYMFUlNGSqSmjJUJDVlqEhqylCR1JShIqmpF7ztaZLjknwryX1J7k3y7t7cF5I80muJetww9Ugav1G0PX0WOK+qXgucAfxVkgN783/aa4l695D1SBqzYUPlLOCabvsa4OypO1TVg1X1ULf9Awa9gX5zyN8raYEaNlRm2/YUgCRrgKXAw73hP+/eFl25pz+QpMk1qrandB0MvwSsr6pfdMMfZhBGSxm0NP0QcPkM6/+/l/Ly/V66r7IljclI2p4meSnwz8ClVbWp99p7znJ2Jvk88MG91PFLvZT3Vbek8RhF29OlwFeAL1bVDVPmVnWPYXA95ntD1iNpzEbR9vR3gLcC50/z0fHfJtkCbAFWAn82ZD2SxiyDPumT5YDlh9aJrzhv3GVoDnZvfWTcJWgO7qrb+GntyHzW+o1aSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNdUkVJKckeSBJFuT/Err0yTLklzfzd+V5Ije3Ie78QeSnN6iHknjM3SoJFkCfAp4J3AMcG6SY6bsdgHwVFW9ErgS+Hi39hhgHbCnz/Knu9eTNKFanKmsAbZW1baqeh64jkGP5b5+z+UbgFO6Xj9nAddV1c6qegTY2r2epAnVIlQOAx7tPX+sG5t2n6raBTwNHDzLtcCg7WmSzUk2P7/7uQZlS3ohTMyF2qraUFWrq2r10iW/Pu5yJM2gRahsB17ee/6ybmzafZLsBxwAPDnLtZImSItQ+Q5wdJIju77J6xj0WO7r91w+B7i9Bq0RNwLruk+HjgSOBr7doCZJY7LfsC9QVbuSXAzcDCwBrq6q+5JcDmyuqo3A54AvJdkK7GAQPHT7/T1wP7ALeF9V7R62JknjYy9ljYS9lCeLvZQlLRiGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmRtX29ANJ7k9yb5Lbkhzem9ud5O7uZ+ofzJY0YYb+w9e9tqfvYNAM7DtJNlbV/b3d/gNYXVXPJnkv8Ang3d3cc1V13LB1SFoYRtL2tKq+UVXPdk83MejvI2kRGlXb074LgJt6z5d37Uw3JTl7pkW2PZUmw9Bvf+Yiye8Bq4G39YYPr6rtSY4Cbk+ypaoenrq2qjYAG2DQomMkBUuas1G1PSXJqcClwNqq2rlnvKq2d4/bgDuA4xvUJGlMRtL2NMnxwGcZBMoTvfEVSZZ12yuBkxh0K5Q0oUbV9vQvgZcAX04C8N9VtRZ4DfDZJL9gEHBXTPnUSNKEaXJNpapuBG6cMnZZb/vUGdbdCRzbogZJC4PfqJXUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqalRtT09P8mPe+1NL+zNrU/yUPezvkU9ksZnVG1PAa6vqounrD0I+AiDXkAFfLdb+9SwdUkaj5G0Pd2L04FbqmpHFyS3AGc0qEnSmLT4a/rTtT09YZr93pXkrcCDwPur6tEZ1k7bMjXJRcBFAMt5Ebu3PtKgdI3KzT+4e9wlaA7WnP7svneawagu1H4dOKKqXsfgbOSaub5AVW2oqtVVtXp/ljUvUFIbI2l7WlVP9lqdXgW8cbZrJU2WUbU9XdV7uhb4frd9M3Ba1/50BXBaNyZpQo2q7emfJFkL7AJ2AOd3a3ck+RiDYAK4vKp2DFuTpPFJVY27hjl7aQ6qE3LKuMvQHHihdrKsOf1RNt/zs8xnrd+oldSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpqVG1Pb2y1/L0wSQ/6c3t7s1tnLpW0mQZSdvTqnp/b/8/Bo7vvcRzVXXcsHVIWhjG0fb0XODaBr9X0gLUIlTm0rr0cOBI4Pbe8PIkm5NsSnL2TL8kyUXdfpt/zs6ZdpM0Zi16Kc/FOuCGqtrdGzu8qrYnOQq4PcmWqnp46sKq2gBsgEGLjtGUK2muRtL2tGcdU976VNX27nEbcAe/fL1F0oQZSdtTgCSvBlYA3+qNrUiyrNteCZwE3D91raTJMaq2pzAIm+vql1sivgb4bJJfMAi4K/qfGkmaPE2uqVTVjcCNU8Yum/L8o9OsuxM4tkUNkhYGv1ErqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJTrdqeXp3kiSTfm2E+ST7ZtUW9N8kbenPrkzzU/axvUY+k8Wl1pvIF4Iy9zL8TOLr7uQj4G4AkBwEfAU5g0OnwI0lWNKpJ0hg0CZWq+iawYy+7nAV8sQY2AQcmWQWcDtxSVTuq6ingFvYeTpIWuFF1KJypNepcWqZexOAsh+W86IWpUtLQJuZCbVVtqKrVVbV6f5aNuxxJMxhVqMzUGnUuLVMlTYBRhcpG4LzuU6A3AU9X1eMMuhqe1rU/XQGc1o1JmlBNrqkkuRY4GViZ5DEGn+jsD1BVn2HQvfBMYCvwLPAH3dyOJB9j0I8Z4PKq2tsFX0kLXKu2p+fuY76A980wdzVwdYs6JI3fxFyolTQZDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTY2q7envdu1OtyS5M8nre3P/1Y3fnWRzi3okjc+o2p4+Arytqo4FPgZsmDL/9qo6rqpWN6pH0pi0+sPX30xyxF7m7+w93cSgv4+kRWgc11QuAG7qPS/gX5N8t2ttKmmCjaqXMgBJ3s4gVN7SG35LVW1P8lvALUn+s2v4PnWtvZSlCTCyM5UkrwOuAs6qqif3jFfV9u7xCeArwJrp1ttLWZoMIwmVJK8A/hH4/ap6sDf+4iS/sWebQdvTaT9BkjQZRtX29DLgYODTSQB2dZ/0HAJ8pRvbD/i7qvqXFjVJGo9RtT29ELhwmvFtwOt/dYWkSeU3aiU1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdTUqHopn5zk6a5f8t1JLuvNnZHkgSRbk1zSoh5J4zOqXsoA/9b1Sz6uqi4HSLIE+BTwTuAY4NwkxzSqSdIYNAmVrqPgjnksXQNsraptVfU8cB1wVouaJI3HKNuenpjkHuAHwAer6j7gMODR3j6PASdMt7jf9hTYeWvdsBibjq0E/mfcRbwQlqxatMe2WI/rVfNdOKpQ+Xfg8Kp6JsmZwFeBo+fyAlW1AdgAkGRz14xsUVmsxwWL99gW83HNd+1IPv2pqp9W1TPd9o3A/klWAtuBl/d2fVk3JmlCjaqX8qHpepsmWdP93ieB7wBHJzkyyVJgHbBxFDVJemGMqpfyOcB7k+wCngPWVVUBu5JcDNwMLAGu7q617MuGFnUvQIv1uGDxHpvHNUUG/29LUht+o1ZSU4aKpKYmIlSSHJTkliQPdY8rZthvd+9WgAV7wXdftyYkWZbk+m7+riRHjL7KuZvFcZ2f5Me9f6MLx1HnXM3iNpQk+WR33PcmecOoa5yPYW6v2auqWvA/wCeAS7rtS4CPz7DfM+OudRbHsgR4GDgKWArcAxwzZZ8/Aj7Tba8Drh933Y2O63zgr8dd6zyO7a3AG4DvzTB/JnATEOBNwF3jrrnRcZ0M/NNcX3cizlQYfHX/mm77GuDsMdYyrNncmtA/3huAU/Z8JL+ALdpbLmrft6GcBXyxBjYBByZZNZrq5m8WxzUvkxIqh1TV4932D4FDZthveZLNSTYlWajBM92tCYfNtE9V7QKeBg4eSXXzN5vjAnhX9xbhhiQvn2Z+Es322CfRiUnuSXJTktfOZsEo7/3ZqyS3AodOM3Vp/0lVVZKZPgc/vKq2JzkKuD3Jlqp6uHWtmrevA9dW1c4kf8jgbOy3x1yTZjav22sWTKhU1akzzSX5UZJVVfV4d1r5xAyvsb173JbkDuB4Bu/zF5LZ3JqwZ5/HkuwHHMDgG8gL2T6Pq6r6x3AVg2tli8GivN2kqn7a274xyaeTrKyqvd5AOSlvfzYC67vt9cDXpu6QZEWSZd32SuAk4P6RVTh7s7k1oX+85wC3V3flbAHb53FNuc6wFvj+COt7IW0Ezus+BXoT8HTv7frE2svtNXs37ivQs7xKfTBwG/AQcCtwUDe+Griq234zsIXBpw5bgAvGXfdejudM4EEGZ1GXdmOXA2u77eXAl4GtwLeBo8Zdc6Pj+gvgvu7f6BvAq8dd8yyP61rgceDnDK6XXAC8B3hPNx8Gf2zs4e6/vdXjrrnRcV3c+/faBLx5Nq/r1/QlNTUpb38kTQhDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrq/wAi4riG5hRTWAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "B = coreg.W.value @ coreg.W.value.T + np.diag(coreg.kappa.value)\n", "print('B =', B)\n", "plt.imshow(B);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "Bonilla, Edwin V., Kian M. Chai, and Christopher Williams. \"Multi-task Gaussian process prediction.\" _Advances in neural information processing systems_. 2008." ] } ], "metadata": { "anaconda-cloud": {}, "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.8" } }, "nbformat": 4, "nbformat_minor": 2 }