{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Conditional Gaussian KDE\n", "Here we go through a simple example of a conditional Gaussian KDE.\n", "Let's imagine we would like to sample from the distribution $P(x|y)$, given that samples from a full distribution $P(x, y)$ are available. Example of such scenario could be an MCMC chain, fitted for both $x$ and $y$.\n", "\n", "In general, such procedure is not possible as one would need to run MCMC from beginning, by sampling $P(x|y)$ directly, for one particular $y$. However, if we make an analytical form of the total pdf, then in principle we could obtain such conditional distribution for any $y$.\n", "\n", "`ConditionalGaussianKernelDensity` class implements these functionalities:\n", "\n", "- fitting the Gaussian KDE for the total distribution $P(x, y)$ form its samples,\n", "- slicing through it to either obtain samples of the $P(x | y)$, or calculate the exact values of the conditional probability." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "plt.rc('text', usetex=True)\n", "plt.rc('font', family='serif')\n", "\n", "from conditional_kde import ConditionalGaussianKernelDensity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example, let's create a simple 2D sample of points and fit KDE with it." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAFDCAYAAACUSq0bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcXklEQVR4nO3dzW+cWXbf8d8pii8SJapEqdXTrRmMw55BxjaSILKE8SqAMdIEWQUw2ONFNgHisP0HJN2Yv2CgTgJklzSzySpBGkR2QRakExhGxkCa0gSJ23HiNN3tGY+slyYpiRLFtzpZ6GF3qe65FIsq3qqH9f0AhFj3uXzqqSrq6Oo599xr7i4AwPFq9PsCAGAYEGwBoACCLQAUQLAFgAIItgBQwKl+X0A/jNm4T2iy35cBHJ5Z3M5sooHyRGsP3f2N6NhQBtsJTer79oN+XwaGnWX+Y+mttOup0bjr7s7hzx2cF7215Atf5I5xGwEACiDYAkABBFsAKIBgCwAFDGWCDAjlElaRXLKpm3N0cR3ZRNhx6cV7gZcwsgWAAgi2AFAAwRYACiDYAkABJMhQP11UXnVVpTUyEnfd2zv0eUcmz4TtredbhztvRq6CzEbjv8LR83WT9Grkzru9fehz4GWMbAGgAIItABRAsAWAAmoTbM3sqpl9Zma3q69bQZ81M1s0s/f7cY0AkFOnBNm0u78jvQi8ktaDPu+6+1LRqwKAQ6hNsO0IojPuvhB0a5rZjLuvdB4wszlJc5I0oThjjAGTyZ7nZg1YI/h1zs0w2NlN2hoT43HfIANv43FfmzoXtjfsSdr3bLyAvT/ZSPuengj75jSa5w/d1zeepo2Z960RtXcxq0KSvJUueN7NzIyuSqUHqJS4NrcR9pnZXCbQStK0pFUz+6jzgLvPu/s1d782qvgvCgAcl9oFW0k3cweqgLouad3MZstdEgAcrFbB1syaBxybq+7lAsDAqVWwVXWboL3BzBarbz+uHs9K0gG3GgCgOPMh3J1zyqadDR9foRfJhi7OkStHjTRyyaIgedOKkj+SRt5MN0D1x2kSS5LsG5fTxlb8OlrNzK7Nwd+zxrO49LU1GeQUduPn8/FMIis4d+vMWNh35C+/DC4i81kH7f5sM+66+Ty+tiARuff0Wfx8XfzORb9DXSXeuny+yJIv3Hb3a9Gxuo1sAaCWCLYAUADBFgAKINgCQAEEWwAooDbluijsNRfizsqV4DYsOG1mpkyuXDeYpTBy8ULY1yfTvjYeZ+u3vtlM2ja+Ffcdexy/FxMP08W8139tKuy7cyZ9L3YyFebnv8hk2z39ga3z8Xs/vZWeIzvL4Un6OiyzYPpI5nPS+qO0b9wznE0QlVpL8ULq0e+V1J9F0BnZAkABBFsAKIBgCwAFEGwBoAASZDi8Y1wbNJsMi2RKML2ZJpzscbo2rCRtfi8t1z21cTrsu30+/Wvy9K048bL2vTjVs30pTajZVu41B6W9F4LdciU9uxKXLu820yTS1KfxtT34fjNpO/1l/FmfvheURI/Ga/iO/VWm/DlYa9camXFfVMab+fzDHYVZzxYAhgvBFgAKINgCQAEEWwAogATZsOtBVVhjLK6mel2Ni9PxgbF47dudi2klU+sbcfLm+YU00bNzJf7rsHM2bXt2JX5/rnzvXtg+PpImrP76+fth34uj6Rq8f/LkrbDvvcvBxUna3Enfo41fXgr77gUFYLlqs81LaedLP4sTYa1z8V5/I2tB0itTvWfRurqZzz9KhnqXy9m+7vq52sn3Z2QLAAUQbAGgAIItABRAsAWAAgi2AFAAsxGGyevumJuZuZBbGzTK1kZrjkpSYyqYNZDpq1Nxu4+k17c7mVspNbUbV+tq49vpe/TtX70b9v3W2fWw/bcu/GnS9g+n4tkIt1a/k7S9f+U/h33/3epvhu1nR9LS1T8YTc8rSXcfNtPGX8ZlwHvBWrvPrsTr2U48iEuMR4LdeHM7G0e7EodluepyJ93c730Xf0d894CpB4FajWzNbM3MFs3s/czxWTO7YWZzpa8NAA5Sq2Ar6V13v+nuH3YeMLNZSXL3perxjdIXBwA5dQu2TTObyRy7Lmml+n5F0tUylwQAr1a3YDstadXMPgqONTseX2x/YGZzZrZsZss7iu/5AMBxqVWCzN3nJcnM1s1s1t0X2g6v60UwPuhn5yVpyqa7WDz1BHndtT17sDaoZZJbu/cfJm2nLsflpbkE2aloM8K9+Jq3p9JsmMdL1OrCn6RjkrtX4s0acx5NpUmkf/v4ctj3O+Npye+nW1fCvj++/Adh+wd/+feStrtfpuvIStLeo7RUduJJ/GaMRMvLjmXeOMu0n5tMu2bKdT3YHNI34jWKXzsBnOvfi41OVaORbTUyPejWwCf6enQ7I2nx2C8KAA6pNsFW0sfSS4mwherxYtvjmSox1txPlAHAIKjNbQR3X5d0p/paaGu/2fb9/iwFAi2AgVKnkS0A1BbBFgAKqM1tBAyAHmRlfStTwnk+ze77TlwOaZmyTNtM+49kMuJjj9NztE7Fpb0bb6fn2P08XrT780vxgtn/8vO0xubMdFq2Kknb2+lfy1+5/GXY919v/Z2w/f4XwcQciyfhXPif6etuZCpRz/08PbA3Hr/Ho3/1OD5J9DuQmTWiYNfl3GL1rZ10gfbs72xO9Lvcox16GdkCQAEEWwAogGALAAUQbAGgABJkddGjksF+8yDhIUl6FtSBjsQJKxuPk1B6lK6J2njjYtBRmnyaJmm2r8blwacfpG3nfhFfwt5ovPPrxjfTJNLzzfj1nV1JP+sv3ooX2z33eZycmo7yPJlfoTP3086tTGQYfZyuXTyxHif6FKwvLEm+FpTg7gbJLUkeJL1y69Za8PvS1Rq33QrLg/PdGdkCQAEEWwAogGALAAUQbAGgAIItABTAbIS6GIRZB11eQ5gJzp0j2olXcSZ570FcujpyIV0c2++mC3FLksbS52v+LM7sW7Ab7N6bzbDv7rm4lPTcL9JzbzUzOw3vpK/70qfxezG2+jxs37qUzl4YfRTvghw9X2P9adhXQQl1K/N56HS8Q69NpLNJ9h5mni/QyOy6HJbrdrOLbrdO6uLhAFBnBFsAKIBgCwAFEGwBoAASZCX0otT2uM5xnIm3bta53U0TL555zbkEyV5UBpop1xxRutut7gV1uYpLjEd2c+fNCNbKncicI+rrZ+JyXfO4PnQiSHopapNkT4NS6a04meZBsjAqk5Uk34yTd3ub6Tm6WaM2W/LdjT4knBnZAkABBFsAKIBgCwAFEGwBoIDaBFsza5rZVTObNbNbmT5rZrZoZu+Xvj4AOEidZiP8SJLcfd7MrpvZnLvPd/R5192X+nBtB+tFyWDmHBaUuUaZ/QOv47h0M/sh6JvLcre240x5pHE6zuK3nqeLh+d27Q1/fiMuL20EpaiS5NE1d7E4eq7s2DPnUPBacll8a8RlypGwJLYHuvlMfXcAStePoDbBtiOwzkhaDLo1zWzG3VcKXRYAHEptbiPsM7MZSauZEey0pFUz+yj4uTkzWzaz5R0F+9YDwDGqXbCVNOvu70UH3H3e3dclrZvZbHDsmrtfG1VmDysAOCa1CrZmNuvuH1bfX+04NtfZBgCDojb3bM3shqRbZvbjqumDqn3R3W9K+ljSzP6I1t0X+nOhXSSFovZc0izTbkHpajZB1sV5e5JM6+b1BX17kQhpBaWhOblEUZS8yZaXBok3Kf5MouRm7hy5ZKFah09Y5UqX/XU3oO1FAviE7B59kNoE2+oe7TtB+83qz3VJd6qv/gRaAMio1W0EAKgrgi0AFECwBYACanPPttZyya1c0iOQS26EyZRM4qWrDRiPywAnPLqpYuqmb05XiczMSrm534tufreO7TPp5rwD/HvRK4xsAaAAgi0AFECwBYACCLYAUADBFgAKYDZCr3WVVQ12Uc1kqHMzDKL1U3Mlo6iv7mYu1HfN15OMkS0AFECwBYACCLYAUADBFgAKIEHWR9FatI2zk3HnLjbly/0LGm46mNGLclT0UTfrKqMIRrYAUADBFgAKINgCQAEEWwAogGALAAUwG+GoMguCN6LdblsenyNY+NmfPYufLrObq04Fz5dZUDo8b272w2owG2EIdkA9MfhMBg4jWwAogGALAAXU6jaCmc1KWpc04+7z3R4HgH6pzci2CqRy96Xq8Y1ujgNAP9VpZHtd0n+ovl+RdFXSUhfHj66L0scoGRatOZs1Fq9ba+fOxv2303VOozJgSfKd3bRt83l8Xso9gZ6qU7Btdjy+2M1xM5uTNCdJEzrTy+sCgFeqzW0EvbgXO33U4+4+7+7X3P3aqLoYaQJAD9Qp2H6ir0evM5IWuzwOAH1Tm2Dr7guSZqrEV7MtEbZ40HEAGAR1umcrd/+w+napre3mQccBYBDUKtj2TZSFz5SudlUqe34qbTwXl89uvRX0lTS+cj9tzJTg2uMnSVuuPDj8+cwOv93u/AoMo9rcRgCAOiPYAkABBFsAKIBgCwAFkCA7jNw6roGRIOllU5lSW0t3zN1+81zYdft85qP6zuWkqbEdl9WOBgmy7Hq2a4+SJhJhwNExsgWAAgi2AFAAwRYACiDYAkABJMjadbGhYSO3AWMjTXppNK68ap1Pl3rcmxgJ+46vxcmpp29nriPQXLuQNv4yqECTNBIkzlqZtW9JnAGvxsgWAAog2AJAAQRbACiAYAsABRBsAaAAZiO0y+0eG8xSaAU71UrSqWDmQTTrQJI2304z/juT8b9/uxPxR2XpZr4ae5LZ+XcsPUcjs2tv6/7DtG9ml+C9p5k1fNmNF/hK1yNbM/tdM/ut47gYADipjnIb4Zqk3zezPTP7xMx+Yma/bWbxVgIAgO6Drbv/nrs3JF2X9LGk35C0IGnNzP7MzP6Vmf2tHl8nANTakRNk7n7H3f+Zu/+wLfj+bP/PatQbrxcIAEOmZwkyd78j6Udm9k/d/ZqZvS/pv+hF8B08UWluLxI6jfS8th0n06Lk1saV+N+/vYnM022FrWHf0Yunk7aJe2vxtZ1KfzVy5boAXu0oCbK/Xd0q+F0z+5Wgy5r01bbi75nZT17zGgGg9o5yG+E9SY8k/Z6kzzqSZL8t6eZ+x2q0u9qLCzWzppldNbNZM7uV6bNmZovVqBoABsZRbiPcdvd/I0lmNqMXwfdm9ednkv5xdeyf6EXi7LPeXKp+JEnuPm9m181szt3nO/q86+5LPXo+AOiZowTbpf1A6u4rkj7I9PuhpN+R1JPbCB2BdUbSYtCtaWYz1XUBwMA4ytSvP3f3fy7pnVf0+6G7X3f3/3jkqwtUo+nVzAh2WtKqmX0U/NycmS2b2fKOwqwSABybI89GcPff7+WFSJKZzepFwGy30hFYZ939vcw1zVfnWTezWXdf6Dg2L0lTNh3MA1BXu+jmSld1Js34b34zU++xl17G5htxVx+NL3lkM12sfOxJsIC5pI2301Li0UcX4/O20pkZvvE0c3GHL3OmhBfDaqDWRmgPjpEqgH5YfX+1SsDtH5uTtNzeBgCDojarfpnZDUm3zOy2md1WNQI2s/17tx9Xj2elVwduAChpoEa2B6luJST3id39ZvXnuqQ71ReBFsBAqc3IFgDqrDYj257rIlHTOJ0mvWw8kyDbTUtzdyfjHXPXvpv+Wzf563ENyPVv/Dxs/99rbyZtdyfTNkk6+0WaOPORLpKCo/Gvi7fi5J3vZda5BYYQI1sAKIBgCwAFEGwBoACCLQAUQLAFgAKGdzZCFzyzk27oXDoToLETz3x4fjmdHfDdZryY9x/+xUzYfnp8J2mb+m58ju37nZXQ0vM34lkVZx+mO//6gy/DvlmU5gJfYWQLAAUQbAGgAIItABRAsAWAAkiQtcusZ2uNNJFlp+Ptbn0sfUt9JF5ftnUhTW794sn5sO/ff+d/he2rO2kia+nTXw37nn+Wtlkuh7WVLrCeW8O39ZzF2IFXYWQLAAUQbAGgAIItABRAsAWAAgi2AFDA8M5G6GIn3VZQrtvYfB6fdjvtO7Kd2Rl3PF1c+x/N/DTs+2gvXcBckv7GmXRR8SWLZyPspZvrqnUqnimhvWB33cxi4L6bzqoA8DJGtgBQAMEWAAog2AJAAQRbACigVgkyM1uTtCxp0d0/DI7PSlqXNOPu8weeLFprNVeuO5Lujts4m5bJSlKUCnvyzfht3nuctv3x0yth33/x1n8L2//9k7T/5PnN+NpOjSVtE19uh32VKc3tSvR+ssYthlTdRrbvuvvNAwKt3H2penyj9MUBQE7dgm3TzOItC6Trklaq71ckXS1zSQDwanULttOSVs3so+BYs+PxxfYHZjZnZstmtrwjVqkCUNZA3bOtbgV0bpS1sn9rYP8+rJmtm9msuy+09VsPfvYr1c/OS9KUTcdVBgBwTAYq2HYEz5eY2ZykZXe/k+nyib4e3c5IWuzt1b0stwmkB8m0Vu5dnkiTRT+9+9fCrn80/UnY/p8e/s2kbXcv/g/L9pX0+Z59I06ETX26HraHctV4JMOAr9TpNsLH0kuJsIXq8WLb45kqMdbcHw0DwCAYqJHtQdx9XdKd6muhrf1m2/f7sxQItAAGSp1GtgBQWwRbACiAYAsABdTmnm0Rmey57x4+q759+UzSNpKZ1nvuf6Tls7/5Dz4P+/7XJ78Wtr8z+SBpu73x7bDv6Xvpv61jj+O1aK0Z7PK7zbq1wFExsgWAAgi2AFAAwRYACiDYAkABJMjaZcpOG9HarpnNDyf+772kbXsqXqN241vp8/3x2lth3+uX/iJs/+8PgmTY8/h1eNA8+jiznu3TZ3F7eOIuynIp7cWQYmQLAAUQbAGgAIItABRAsAWAAgi2AFAAsxHaZTLiredpve3IZLpIuCTpVPqW7kxa3HUjbfv5F5fCvu9MfRm2T40/T9rW/zz+WCfvdrFBxXhaSrx3/+Hhfz6HWQcYUoxsAaAAgi0AFECwBYACCLYAUAAJssOIkjpjo3Hf0fQtvfhH98OuW3/3zaStsRkn3v7wp78etltwaZNpzkySNL6elhifehhk6ST52qP0uYKdg6Xu1vsFhhUjWwAogGALAAUQbAGggNoEWzO7amafmdnt6utW0GfNzBbN7P1+XCMA5NQpQTbt7u9ILwKvpPWgz7vuvlT0qgDgEGoTbDuC6Iy7LwTdmmY24+4rnQfMbE7SnCRNKN0B90DBgtd7QbZeksJ8/duXw76nV9Py2Y2tTGnvs7h99HHaduH/xLvgjq8G0xRamZkEnl6bZxZMZ0Fw4NVqcxthn5nNZQKtJE1LWjWzjzoPuPu8u19z92ujCnZeAIBjNFAjWzOb1YuA2W6lY1R7U9J89PPuPl+dZ93MZg8IygBQ1EAF21cFRzNrHnBsTtKyu9/p9XUBwOuq222EaUmr7Q1mtlh9+3H1eFZ6deAGgJIGamT7KlXi672OtpvVn+uS7lRfvQ20QaJn5Ny5uG8rTSzZg7Ww6/k/Td9+b0yFfZ83M9cW/HO5dSEuqz3zZ0Fp7sbT+LwWJ+RCJMKAV6rbyBYAaolgCwAFEGwBoACCLQAUUKsEWd8EFVKtzXjR2MbpiaDx8P+mNXbiTRnPZjZrPPU0reoa2clsXHn2dNJma+thX9/dDRoziTAqyIBXYmQLAAUQbAGgAIItABRAsAWAAgi2AFAAsxEOIdpV1hpxOaudnUzafGsr7ruZrjs79f/i3W5bo3EJbmM7nTXQePQs7BvNivCdYNaBpNbz4JqZdQAcGSNbACiAYAsABRBsAaAAgi0AFECC7BB8N01keS5ZtJ5uBNl442J83l/cTfteOB/2zf2r6FNBQm5sNO4cPJ+dSUt4JeXXuQVwJIxsAaAAgi0AFECwBYACCLYAUADBFgAKYDZCAXu/vBe2N6bOJm2emQVg4+Nx+5OgNPfZZtjXg7LjVmbx8KhE2ffShcoBHA4jWwAoYGCDrZnNmtli0HbDzOYO+JnscQDol4ENtu6+0P7YzGar9qXq8Y1ujgNAPw1ssA1cl7RSfb8i6Wo3x81szsyWzWx5R/GShwBwXOqUIGt2PO6sgT3wuLvPS5qXpCmbjreq7UZmDddoHdgo2SRJ2k7LgLNPl2m3aBfcsbGwb2stLSXuKunFurXAkfUt2Fb/7Z/uaF7Zvw0QWA/6d3McAPqmb8G2857sIXyir0evM5IWuzwOAH0zsPdsqwTXtbbE14Kkmaq92ZYIWzzoOAAMAnN//duXdTNl0/59+8HxnDxYejF3z7YxERcqhDLLJobnztyz3bv3IGnjni3QO0u+cNvdr0XHBnZkCwAnSZ1mI9RWbifevY10J91Tb1wK+7YePY5PfuZMcOJ4tBpdh+cGttEolt11gSNjZAsABRBsAaAAgi0AFECwBYACSJD1WDQVq7W9Hfc9FUznasT//jXOT4Xtuw9X08ZuEla5pFfUTiIMODJGtgBQAMEWAAog2AJAAQRbACiABFmP+e7h16iNKrr2vlzr6vlGJtMKsr2nwSaQUpjgyq3bwOaOQG8xsgWAAgi2AFAAwRYACiDYAkABBFsAKIDZCH0UlvF2uWbs3sbhZz+Ep+1i9gSAo2NkCwAFEGwBoACCLQAUQLAFgAIGNtia2ayZLbY9bprZ1ar9VuZn1sxs0czeL3elvWUjI+FX/gca6VdGY2ws+QJQxsAGW3df6Gj6kaRr++1mNhf82LvuftPdPzz2CwSALtRm6pe7z7c9nJG0GHRrmtmMu68UuiwAOJSBHdnmmNmMpFV3XwoOT0taNbOPgp+bM7NlM1ve0daxXycAtOvbyNbMZvUiOLZbyQTRdrPu/l50YH/0a2brZjbbfiuiOjYvSVM27Ue/cgDoXt+CbXBP9pWqAPph9f1Vd7/TdmxO0nJ7GwAMioG9jWBmNyRdq0bA+49vmdltM7utalTcNmPh4+rxrHS0YD4IfG8v/ApnHVjjRRlv51dGa2c3+QJQhrkP3/+op2zav28/6PdlxA6YuhV63W3L2Z4c6JklX7jt7teiYwM7sgWAk4RgCwAFEGwBoIChvGdrZg8kfdHv63hNlyQ97PdF9MEwvu5hfM1SPV/3t939jejAUAbbk8DMlnM34k+yYXzdw/iapZP3urmNAAAFEGwBoACCbX3Nv7rLiTSMr3sYX7N0wl4392wBoABGtgBQAMEWAAog2NZQtTXQjcxuFSfSSdjy6LA6t4RqazvRn3nmdZ+Yz51gWzNtq5otVY9v9PeKihmaLY86V6wbls88s1LfifncCbb1c13S/rY/K5Ku9vFaSmpWu3QMo2H9zKUT9LkTbOun2fH4Yj8uog+yWx4NgWbH42H5zKUT9LkTbOtnXel2Qieeu8+7+7qk9f3/Vg+RdQ3hZy6drM+dYFs/n+jrkU5ul+ETpdqsc5j+69xp6D5z6eR97gTbmqmSCDNVkqR5iA0yT4ITseXRYXVuCTUsn3nn69YJ+9ypIAOAAhjZAkABBFsAKIBgCwAFEGwBoACCLQAUQLAFgAIItgBQAMEWAAog2AJAAQRbACjgVL8vABg01fqpN1Qt/uLuH7btkPAb7v5ev64N9cXaCEAbM2tK+rG7f1A9vi1pWdIHerHi1m1J77j7SvYkQIDbCMDL5iT9pO3xqqTp/TVVJX1AoMVRMLIF2pjZTHswNTPXi32war28H/qPYAtkVOurLrq79ftaUH/cRgDybkp6aaHuk7L5IMoj2AJtzOz9toezku60HbuqF0kyoGsEW6BSbb/yYzNrVrcQ7nR0+Z2TuiUNjh/3bIHK/rQvSZ9JWnH3pWoL7dvSi51e+3h5qDmCLQAUwG0EACiAYAsABRBsAaAAgi0AFECwBYACCLYAUADBFgAKINgCQAEEWwAo4P8DZ9zXxjXyQdEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mean1, cov1 = [0, -5], [[1, 0], [0, 5]]\n", "mean2, cov2 = [10, 5], [[5, 0], [0, 1]]\n", "data_xy = np.concatenate(\n", " (\n", " np.random.multivariate_normal(mean1, cov1, 10000), \n", " np.random.multivariate_normal(mean2, cov2, 10000)\n", " ), \n", " axis = 0\n", ")\n", "data_y = data_xy[:, 1, np.newaxis]\n", "plt.figure(figsize = (5, 5))\n", "plt.hist2d(data_xy[:, 0], data_xy[:, 1], bins = 50)\n", "plt.xlabel(\"$x$\", fontsize = 20)\n", "plt.ylabel(\"$y$\", fontsize = 20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Two most important parameters of the KDE are `whitening_algorithm` and `bandwidth`.\n", "Whitening is the process of normalizing the data in some way, here we implement three different ones:\n", "- `None` - no normalization\n", "- `\"rescale\"` - rescaling every dimension to be of unit-variance\n", "- `\"ZCA\"` - whitening along principal component axes\n", "\n", "Bandwidth on the other hand controls the size of Gaussians around each point. Higher bandwidth will give wider and smoother distributions. We implement several ways of choosing its value:\n", "- `\"scott\"` - uses approximate relation depending on the number of samples and dimensions - Scott's parameter\n", "- `\"optimized\"` - searches for the best bandwidth by minimizing KL divergence between real distribution samples and fitted KDE\n", "\n", "Besides that, we allow to specify its value explicitly." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fitting 5 folds for each of 10 candidates, totalling 50 fits\n" ] } ], "source": [ "kde = ConditionalGaussianKernelDensity(\n", " whitening_algorithm = \"rescale\",\n", " bandwidth = \"optimized\",\n", " steps = 10,\n", " cv_fold = 5,\n", " n_jobs = -1,\n", " verbose = 1,\n", ")\n", "kde = kde.fit(\n", " data_xy,\n", " features = [\"x\", \"y\"],\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can either pull samples for some conditional values of y, or calculate analytical distributions for it." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA14AAAFCCAYAAADhUww7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABX9ElEQVR4nO3deXxU9b3/8dcJSSAhkBC2hCUTQkWxuLBoq1Z/WqO2rq2iqK21tRWsG9aFXG5v1VotwlVr1aqAdrUqQqu9rVYlWupWFRhUVFQkZNgCCUkGEyD79/fHOZN9z8ycWd7PxyMPwmTmnM9kkm/mfb6bZYxBREREREREQifB7QJERERERERinYKXiIiIiIhIiCl4iYiIiIiIhJiCl4iIiIiISIgpeImIiIiIiISYgpeIiIiIiEiIKXhFOMuy5g7gsYsty8oIYjlRZyDfPxHpmtqm/lO7JOIetV1qg9yU6HYBsc754b4QyAeKgMJWX84EMoACY4y3k8cuBhYN4PR5zjn8AzhGtHvGsqzFxpgCtwsRiSRqm1yldkmkn9R2BYXaIJcoeIWYMWaZZVnrgPXYDcGq1l+3LGs2sN6yrJmtGwnLsvKBtcYYf1gLjiKWZeUBK4GlwDrsBvE0YGnge2mM8VuWtcWyrHxjTGHXRxOJL2qbes+yrJXAos7eyHVyX7VLIiGktqt7aoMim4Yahke+82+HH+5WDca8dl+a174xkU5lAIuxG+B5tGpYAowxy+j4/RURtU1dsiwrz7Kspc4V8nzsq9y9lYHaJZFQUtvVvQzUBkUk9XiFxzFAUW+vsliWNQNYG9KKYseFvbkKDay1LGtGL+8rEi/UNnXBGFOE86bEuYLeF2qXREJLbVf31AZFKPV4hUc+nVyVgTYTHFe2unkOEC9XZcKlEF3ZEWlPbZO71C6J9I/aruBQGxRm6vEKMWesbQawuou7LAaWtRtjm9/VhEfnePnOMTHGLGnVyMw0xnT7C9Tq8TNp1/VsWdZSYL3T/Rw1nBWG8gC/c5W6A2OM17KsWWEtTCSCqW0KLbVLIqGhtqt31AZFJvV4hV6HcciWZWVYlpVvWdZq7Anbvbra4PwSzTPGLDPGLAHmOL/Uz2BPoJzrNADdme00AH46XuWYi71CUDSZA8zCrjvDsqyV3Sz12tXtIvFIbVPoqF0SCR21XT1TGxSh1OMVeqdh/zLOtSwrcJsfe2zyae3v7PxiVHRxrLm0XQa1Ash0VqfxY6/u0+UvuLOiT6CrPR97xZvA12YA9LS6jbO6V0Z392lnXnc1DYQxpsiyrEWtxnh7nUZ3OfZSs+0VWZaVF6p6RKKM2qYQtAVql0RCTm1XNzWpDYpsCl6hlw8UOldSeqO7/SFWtZtImo/zS+T8wnR7jsAvv3P1ZgZtx0fnAz1OrjTGdPZL65pOJtYWAksty8qI9SVjRQZIbVOIqF0SCSm1XT0f09/uJrVBEUJDDUOoF+OQO5PR1RdaX41wrrK0Xja1L+YB3nZXN06ji4mqoeQs2by+Dx+LWz22s53XA1e1ehoaIBK31DaFjtolkdBR29UztUGRTT1eodXlPhPd8NO7LucOv9B96CrubDWgNl3kXXG6xPui225652sz+3jMQOO71LKswi6O39U5uxpuIBJP1Db10Db1h9olkZBT29VN26U2KPIpeIXWaXSzmkxnnLG5nW7UaVnWglZd67NptTSqM5Y4k95N4syj1X4WgXHI9KIhi5ThPM73qbNxzhdhX3Xyd/KwPHWxiwBqm0JC7ZJIyKnt6v5YaoMinIYahlaX+0z0lWVv4LkwsHIPHccNz+lpAmcr7X8hF9KHjQgjSEXr1YYCqxMBV3Zxf38YahKJBmqbgsB5zuutthssq10SCR21XT1TGxTB1OMVAs48pMA45BnO/xf14RdwXSfd24XAMuyrFkXGmAsty1oaGMvb1f4UXbgSu7HJdGrMI4rmUAQYY1ZZljW71Zueydi7tXe4OuVcfVoR1gJFIozapt5z3qwsbFXHYsuyCoHV7d6M5WFfFQfULomEgtqu3lMbFNksY4zbNUg7zpWXGX1Ysaer46ykF/MYLMsy2JsE9rj6TrSyLGsB9ipIMfscRUJNbVNwqV0SCQ+1XZ1TGxR+GmoYgZyrqceE4tiWZa12ftEC/1+APe431n/pjomD5ygSUmqbgk7tkkgYqO3qktqgMFPwilxL280bCJZMnMmjThfzHODUEJwnYjjfxx5XFhKRXlHbFARql0TCTm1XK2qD3KHgFaECV2eceQbBtAh7fPQC7ImjMyN14nowON+/Y/owQVZEuqG2aeDULomEn9quFmqD3KPFNSKYMabAmeS5rJ+HKKLd3gz93Bgwml3UxwmyItIDtU0DpnZJxAVqu5qpDXJJvxfXGDVqlMnNzQ1uNSLiqvXr1+81xox2u46BUNskEnvUNolIpOpL+9TvHq/c3FzWrVvX34eLSASyLMvndg0DpbZJJPaobRKRSNWX9klzvEREREREREJMwUtERERERCTEFLxERERERERCTKsaioiIiIhEgfr6enbs2EFNTY3bpcSdIUOGMGHCBJKSkvp9DAUvEREREZEosGPHDoYNG0Zubi6WZbldTtwwxlBeXs6OHTuYNGlSv4+joYYiIiIiIlGgpqaGkSNHKnSFmWVZjBw5csA9jQpeIiIiIiJRQqHLHcH4vit4iYiIiIhEo9xcsKzgfWiT75BS8BIRERERiUY+HxgTvA+fe3uVFxYWctppp7l2/nBQ8BIREREREVfl5+eTkZHhdhkhpVUNRURERESkR16vl4qKCvx+PxkZGeTn51NYWEhRURF5eXnN/1+8eDEFBQV4vV7y8vLIyMhg5cqVLF68mHXr1lFQUMDixYvbPK41v9/PsmXLmDFjBkVFRcyaNavDeaORerz6KTdXw2BFREREJH6sWLECgNmzZ5OXl0dRURFFRUXMnTuXxYsXA3bPVVFREfn5+cyePZulS5eSn5/PzJkzWbduHfn5+WRmZpKfn8/cuXOZN29eh/MsWrSI/Px88vPzWb9+fYfzRisFr34KDIFV+BIRERGReLBw4UKWLl3K5MmT8fv95OXlMXfuXPx+f5v7zZgxo8PnmZmZzbe1HlIYCHCtBXrWvF4v8+bN63DeaKXgNQDFxa7OQRQRERERCZvCwkJWrlzJ+vXrKSwsxOv1smTJkj4fp3V4Cgw3bC2wyMaMGTPIy8vrcN5opTleIiIiIiLRyOOxl4EP5vG6sXbtWsDupZo9ezZer5eMjIzm8LRq1Sry8vLwer0UFRU1h7OioiJWr17dPD8r0Ju1bt06li5dCti9XIH7LliwoE2ga3/eaGUZY/r1wFmzZpl169YFuZzoYVn2qpuBf0VigWVZ640xs9yuYyDivW0SiUVqm0RsmzZtYurUqW6XMWAXXnghK1eudLuMPuvs+9+X9klDDUVEREREJCxa94LFGw01FBERERGRsMjPz2fLli1ul+EK9XiJiIiIiIiEmIKXiIiIiIhIiCl49YI2SxYRERERkYFQ8OoFn0/7dYmIiIhIZMnNtVfYDtaHOhpCS8FLRERERCQK+Xz2tkbB+uipo8Hr9TJz5kwKCgqaN0G+8MILm/fcCnx9yZIlrFq1imXLlrFs2bIOj1+yZAmFhYWsWrWKmTNnhurb00FhYWHz5sxu0KqGPdizZw/wL+BsIM3lakRERERE3DFjxgzy8vKYM2dO88bJCxcuZMaMGW2+np+f33xbIHzNnTu306/n5eXh9/vJyMgIef35+fnNGza7QcGrB3feeSfwEDCcG274AXfccQfDhw93uywREREREdd4vV6A5gDVlYsuuoiZM2cyd+7cDl9btWoVs2fP7nDcioqK5jCWn59PYWEhRUVFzaGtsLCQxYsXU1BQgNfrJS8vj4yMDFauXMnixYtZt24dBQUFLF68uM3jWvP7/SxbtowZM2ZQVFTErFmzOpw32BS8evDvf/8bmA5M5YEHHmDYsGH84he/cLssERERERFXrFu3jvXr1zN58uQeg1egZ6y1wsJC1q1b1+n9V6xYwWmnncbs2bMpKipq/pg7dy6nnXYa+fn55OfnM2/ePPLz88nLy2PevHmsXr2aoqIi1q1bR35+PpmZmc3hafLkyR32Dlu0aBFz5sxhxowZzJs3jy1btrQ5byhojlc3Kioq2LhxI3AB8GeOOOKILn9IRERERETiQWZmJkuXLmXFihU9hhS/309eXl6b2/Lz85k7d25zMAr0ngEsXLiQpUuXMnny5ObHzp07t3lOWUDrwBf4PDMzs/m21kMX8/LyOtQZ6Fnzer3Mmzevw3lDQcGrG2+88QYAY8eeBMD06dPZsGGDmyWJiIiIiLgqEKSWL1/OhRde2O19n3nmGQoKCro8jt/vp6Kiovm2wsJCVq5cyfr16yksLMTr9TYv3tEXrcNTYLhha4FFNgLzztqfNxQ01LAbr732GjCY4uJjSEmxg9cf/vAHdu/eDWS5XZ6IiIiIxDGPx14GPpjH647X66WoqIgVK1aQl5fXHJzmzZvXPJ+qqKiIwsJCKioqmnuZAvO7vF4vXq+3+et+v59FixaxcuXK5nOsXbsWsEPZ7Nmz8Xq9zcMV8/LyWLVqFXl5ec21BMJZUVERq1evbp6fFejNWrduXfOCGoHzFxUVsWDBgjaBrv15Q8EyxvTrgbNmzTKxPuzu2GOPZe3aFIz5N5YFa9b8m5NPPpkXXniBM8/8JsbYP+z9/BaKRBzLstYbY2a5XcdAxEPbJBJv1DaJ2DZt2sTUqVPdLiMqXHjhhW0CXTB09v3vS/ukoYZdqKqqcsabntR822WXHQ3ApZdquKGIiIiISCRq3QsWSTTUsAv/+c9/aGxsBE4EAl2v6c7kPAUvEREREZFIlJ+f32EVw0ig4NWF1157jUGDBtHYeBxYFsXO7bOBIpzBtJYFGMjNheLizg4jIiIiIiKioYZdee2115ylKYfZk7icj+l33glsAfa1TO7y+VysVEREREREIp2CVycaGxt5/fV3+eSTr3X42tFHH+189kFYaxIRERERkeil4NWJHTt2ALVUVR3e4WvTp093PtM8LxGJTLm5YFl/ISPjTrdLEREREYeCVydaJuPldfhadnY2MAYFLxGJVD4fnHXW79i37363SxERkRji9XqZOXMmBQUFzRsUX3jhhc37YQW+vmTJElatWsWyZctYtmxZh8cvWbKEwsJCVq1axcyZM4NaY2FhYfPmyF1ZtWpV8+dFRUU9bgIdLApenWhZerJj8LIsC5iOgpeIRLJt27YB5dTV1bldioiIxIgZM2aQl5fHnDlzmjc1XrhwIQsWLGjz9fz8fGbPnt28cXIgfLX+euA+y5cvbw5xwZCfn09GRkaXX/f7/axevbr5/3l5eUHf76srWtWwE3bwSgQmdnGPo4BXneXmB4WtLhGR3rKDF+zZs4eJE7tqy0REJFrdcMMNvPfee0E95tFHH83999/fq/va+93iLEbXtYsuuoiZM2c2h7DWVq1axezZszvcXlhYSFFRUXNIKywsZPHixRQUFOD1esnPz28+b/v7tj/+okWLeOWVVygsLGT16tXMmzePdevWUVhYSH5+Pl6vl4KCguYwtmTJEvLz86moqOhwvIFSj1cn7OCVi8czCE/zQvKt5QL17NmzJ6x1iYj0zhfs27cPgN27d7tci4iIxJp169axdOlSCgsLe7xvoGestcLCQpYtW0ZFRUWH+xcVFVFUVMTcuXNZvHgxYPdiFRUVNfeSrVixosv7tjZ79mzy8vKa61i6dGmbXjewg2Ogh2zVqlXk5eUxY8aMNr1iwaIer07YPxx59tZc1iTAtLtHDhC4ojwurLWJiPRse/Nnxx5bAtibwGu7QRGR2NHbnqlQyMzMZOnSpcycObNNuOmM3+/v8PVAj1UgkHm93uYerLy8PObOndth+GFnPWtd3be1efPmsWzZsh575oDmHjGg0yA3UOrx6kQgeHWtdfASEYk0LW3T0qW7MUbbDYqISPAEgtTy5ct7XJjimWeeoaCgoMvj+P3+Nj1fXq+3ebGOnvTmvvn5+axYsYJZs2Z1+vjWJk+e3FxLMOedBSh4tbNv3z7Ky8tR8BKR6NXSNpWUlLhYh4iIxBKv10tRURErVqxo7sny+/3MmzcPv9/f/PXCwsLm4YRA8/wur9eL1+tt/vqqVas49dRT2/SIFRUVNQ9PzMvLY9WqVc2PCxzb6/Xi9/t7vG/A4sWL2yy4kZeXR2FhIXl5eW3uv2DBAlavXo3X62XdunVB//5ZxrQfRtc7s2bNMqEoyG3vvfees1fXKoy5ACwL2n2PLAsgneuuu5wHH3wAQ8f7iEQjy7LWG2M6XhKKIrHaNvWFZf2UQYMWM3z4cObMmcMjjzzSWVMmEjXUNonYNm3axNSpU90uI2osWbKEBQsWtBnKOBCdff/70j6px6ud7paSbytHPV4iEqG2MWHCBMaPH68eLxERiVszZsygsLAwKKErGLS4Rjs//GHPwcvjgdJSBS8RiVTbyMnJYciQIVrVUERE4lawl4MfKPV4teP3FwGZQHqX9ykuhssvV/ASkUi1nZycHLKystTjJSISY/o7TUgGJhjfd/V4ddDTioa2nJwcZxGO/SGvSESkt+yN3XeQk5NDY2Mju3fvdv5YWG6XJiIiAzRkyBDKy8sZOXIklqV2PVyMMZSXlzNkyJABHUfBq4MiYHqP98rJyXE+297t/UREwsne2L2eiRMnUltbS11dHZWVldg9+SIiEs0mTJjAjh07KCsrc7uUuDNkyBAmTJgwoGMoeLViXykuBi7o8b4twUvDDUUkcgSGQOfk5FBVVQXgzPNS8BIRiXZJSUlMmjTJ7TKknzTHq5WdO3cC9fR2qKFNwUtEIkfr4JWdnQ1oLy8REZFIoB6vVo491l7RcMyYPFJSur/vuHHjSEhIoKlJwUtEIkfr4JWcnAwoeImIiEQCBa9W9uyxg9fbb+fRUy9uUlIS48aNY8cOBS8RiRx28BpOenp688RrLSkvIiLiPg01bKOYhIQEJk6c2Kt7ezweNNRQRCLJ9u3bAXso9LBhw0hJSVGPl4iISARQ8GqjhDFjxpCY2LuOQHuel4KXiEQOu8fLvnhkWRbZ2dnq8RIREYkACl5t7GLcuHG9vrcdvLbTFLqCRET6xA5eOc3/1ybKIiIikUHBq42S5lXAesMOXnXsCV1BIiK9duDAAfbu3Uvr4KUeLxERkcig4NVGST96vDTYUEQiQ0vAamnHsrOz1eMlIiISARS8HA0NDcCefvR4wdk8FKKqRER6L7BhMgxvvi0rKwu/3w/UuFGSiIiIOBS8HKWlpYDpU/AKrH64V29oRCQCtASvYc23tbRpGm4oIiLiJgUvx65duwD6NNQwIyODlJQUhrGJ3NwQFSYi0kstwSut+basrCznMwUvERERNyl4OQJzIPrS42VZFuPHj+csHsfnC1VlIiK9U11d7XzWWY+X5nmJiIi4ScHL0Z8er8D9d4aiIBGRPupsqGFmZqbzmT/c5YiIiEgrCl4Ou8fLYuzYsX163Lhx49gVmpJERPqks+CVlhYYdljd4f4iIiISPgpeDjt4jSYxMbFPjxs/frwTvEwIqhIR6b3O5ngpeImIiEQGBS+HPdSwb8MMwe7xOgjAviBXJCLSN1VVVSQnJwPJzbclJyc7F5QUvERERNyk4OWwe7x6v7BGQMucMM30EhF3VVdXM2zYsDa3WZbl9HopeImIiLhJwcsxkB4v5whBrUdEpK+qqqo6BC8IDDes6vgAERERCZv4Dl65uWBZNFoWewI9XpbV9sPj6fYQ48ePdz5T8BIRd3UfvNTjJSIi4qa+rSQRa3w+MIay3btpys4GssH0bZGMlj1yFLxExF1VVVWtFtNooeAlIiLivvgOXo7AHl79GWqYmppKBuDXHC8Rcdnrr1dTV5feoaNewUtERMR98T3U0GEvrAH9WVwDAnFNPV4i4q66uiouuGAYxcVtb1fwEhERcZ+CFwPr8QKwZ3ntIjfXnjYmIuIOzfESERGJVBpqSOser7H9enygx8vnC1JBIiL9ouAlIiISqdTjhR28Ro0aRetNR/vCDl4l5OQ0BbEqEZHeM8YAnS+uYYcxBS8RERE3KXhhDzVs2Y+r7+yhhg28+25ZsEoSEemT2tpaoLHbHi/Tx1VbRUREJHgUvLB7vFqWhe+75i2Ud2mBDRFxR1WVvUFy18HLcPDgwTBXJSIiIgEKXgy8x0vBS0Tc1nPwgupqDTcUERFxS9wHL2MMe/bsISsrq9/HCASvnTu1l5eIuEPBS0REJLLFffDy+/00NDQwZsyYfh8jC7AsSz1eIuKaQKjqbHENBS8RERH3xX3wKi0tBWD06NH9PkYSMGbMGAUvEXGNerxEREQiW9wHr7IyeyXCgQQvgHHjxil4iYhrFLxEREQim4JXEIOX5niJiFt6E7wC9xEREZHwU/AKUvAaP368gpeIuCYQqrqb4zV7djW5ueGsSkRERAIUvIIUvLKzs51j1QehKhGRvgkMI+yux+uhh6rx+cJaloiIiDgUvMrKGDZsGIMHDx7QcVr2Adsz8KJERPrI7vFKJjk5ucPXNMdLRETEfQpeZWWMGTOG3FzwePp/nOzsbACysnZpKI+IhJ0dvDr2dgGkpKRgWZaCl4iIiIviPniVlpYyevRofD4oLu7/cQI9Xo8+WqKhPCISdt0FL8uySEtLU/ASERFxUdwHr7KysgHP74KWHi8tKS8ibrBDVceFNQIUvERERNyl4BWk4DVmzBgSEhIoKSkJQlUiIn3TXY8XKHiJiIi4La6DlyF4wSsxMZExY8aox0tEXKHgJSIiEtniOnh9AdTX1wcleIE9z0s9XiLiBgUvERGRyBbXwavM+TdYwSs7O1s9XiLiip6C17BhwxS8REREXBTXwavU+feWW0YPaCn5APV4iYhbtLiGiIhIZEt0uwA3BXq8SkvHYMzAj5ednU1paSnQQJx/a0UkjIwxGmooIiIS4eK6x6us+bPgzfEyxgB7gnI8EZHeqK2tpaGhAQUvERGRyKXgBQQreAX28gINNxSR8LF7u6B3wSsI3fsiIiLSZ3EfvIYOHQqkBOV448aNcz7TAhsiEj4tPVndz/Gye8Vqw1KTiIiItBX3wStYKxqCerxExB297fGyabihiIiIGxS8ghi8xo4di2VZqMdLRMJJwUtERCTyxXXwKiW4wSsxMZExY8agHi8RCScFLxERkcgX18Er2D1eEJjnpR4vEQmfljleCl4iIiKRKm6DlzGGMnB6qILHnuelHi8RCZ+WHq/uF9ewKXiJiIi4IW6DV3V1NbUEv8fLDl7q8RKR8NFQQxERkciX6HYBbikrs3fxCkrw8njAsgCwF5RPoMGy2n5zPR4oLh74uURE2lHwEhERiXxx2+MVCF4LFozG4xngwYqLwRgwhuyHHwaaKN25s/k2jAGfb6Ali4h0qqqqiuTkZCC5y/soeImIiLgrboNXaWkpAGVlo4PaERXYRLmkRPO8RCQ8qqurGTas694uUPASERFxW9wGr0CPF4Rijhfs2qV5XiISHl988UWPwSs1NdX5TMFLRETEDXEbvPbu3et8Forl5BW8RCR8/H4/I0aM6PY+gwYNcsKXgpeIiIgb4jp4WQwmJ2doUI87duxYwKKgYDe5uUE9tIhIp/x+P+np6T3ezx5uqOAlIiLihrgOXoZR+HxWUI+blJQEjGLfvhKtpyEiYeH3+8nIyOjxfvZwRAUvERERN8R18IJRITp6FtpEWUTCpbfBSz1eIiIi7lHwColsYHeIji0i0lbfgldVj/cTERGR4FPwCols1OMlIuHQ2NhIVVWVerxEREQinIJXSGQBu8nJMVpgQ0RC6osvvgBQj5eIiEiES3S7ADc0NDTg9/sJbY9XPV5vBaNGjQzROUREcNqy3gWvoUOHAvtDWo+IiIh0Li57vCorKzHGENrgBSUlGm4oIqHVl+Bl93gpeImIiLghLoNXy+bJCl4iEt2++U0/AN/6VjoeT/f31RwvERER98R58ArVMMAsQMFLREJvzx4/ABs2ZFBc3P197aGGtTQ0NIS6LBEREWknzoNXaHq8Jk60e7x279aS8iISan6gL0MNYf9+DTcUEREJNwWvENi2LY20tDT1eIlIGPiBvgWv6moNNxQREQm3OA9eoVtxMDs7W8FLRMLAj2VZDB8+vMd72kMN1eMlIiLihrgNXqmpqUBqyM6RlZWl4CUiYeBn+PDhJCT03Jyrx0tERMQ9cRm8ysvLGTUqVCsa2rKzszXHS0TCwN+rYYag4CUiIuKmuAxee/fuDUvwCvR45eZCLltDej4RiVe9D14aaigiIuKeRLcLcEM4gldWVhZVVVXAfny+oUBuSM8nIvFqn3q8REREokDc9niNHBm6hTXA7vGyabihiISSn/T09F7dU8FLRETEPXEbvMIx1NCmBTZEJJQ01FBERCQaxF3wqq+vZ9++fQpeIhIjtLiGiIhINIi74FVeXg4QljleNgUvEQmNpqYm4IteB68hQ4YACSxcWE1ubggLExERkQ7iLngFNk8OdfAaOXIkiYmJaI6XiITKF198AZheBy/Lshg2bCg33LAfny+kpYmIiEg7cRe8wtXjlZCQ4PR6qcdLRELD7/cD9Dp4gT3cUEMNRUREwi/ugle4erwABS8RCan+Bi8triEiIhJ+Cl4hZC+woaGGIhIa/QleQ4cOVY+XiIiIC+I2eIV6Hy8IBC/1eIlIaGiooYiISPSIy+CVlpbG4MGDQ34ue6hhGdAQ8nOJSPwJBK/ebqAMGmooIiLilrgMXgcPjgrLUsp28DLY4UtEJLj27dsHaKihiIhINIjL4NXYOCosSym3bKKseV4iEnyBHq/hw4f3+jEaaigiIuKORLcLCKfcXCgp2QuEfmENaL2JsoKXiASfHbyGOXsG9o6GGoqIiLgjrnq8fD6oqytHwUtEYoEdvDL69BgNNRQREXFHXAUvWxkQ+hUNAcaOHet8puAlIsHXn+CVlpZGbW0tWvRHREQkvOIseB0EqoGxPd0xKFJSUoB0FLxEJBT6G7xsGm4oIiISTnEWvAKrC44J4zmz0F5eIhIK/R1qaNNwQxERkXCKs+BV6vwb7uClHi8RCT47ePV+Dy9o3eOl4CUiIhJOCl4hp+AlIqGhoYYiIiLRQ8ErxIYNy8KyFLxEJLiampqcDZQz+vQ4DTUUERFxR9wGL48HPBSH/Iz//d9ZGFOFri6LSDBVVVVhjKH/PV4KXiIiIuEUh8ErFRhKcTEUMynkZ2zZy2tPyM8lIvHD5/M5n2X36XEaaigiIuKOOAxe4ZzfpU2URWTgcnPBsuyP3Fz7Nq/X63x1Rp+OpaGGIiIi7kh0u4DwUvASkejj84Ex9ueWZf/r9XpJS0ujunpKn46lHi8RERF3xGHwGh/WMyp4ichA2PO4rLY3WhbLmEktR+NhG1i9Hzad1vxZdUuK83iguHjgxYqIiEiXNNQwxEaPHk1CQgIKXiLSV7fffjtTpkwBdrW5vbGhgVo2MX/+TIpNrt0d1suPwU1NTptU3XJ783wxERGRzgWGvQeGvEvfxU3wsq8ahz94DRo0iNGjR6PgJSJ98fHHH3PXXXfx+eefA9/i4MGDzV/79NNPgQPMmNG3+V0AlmU5ww011FBERHoWCFyga3UDFTfBy97vpp5wBy8IDDdU8BKR3jHGcP3115OWlsby5cuBtcydO9e5gNSysMbMmTP7dXw7eGlxDRER6VlgnnH7EenqAeu7uJnjVVoa2MNrdNjPnZWVxfvvl4T9vCISnZ599lleeeUVRox4kCuv/BHp6bt54omfcfLJJwM/ZP369UAKhx56aL+Ob69sqOAlIiL9FwhkltXzfcUWNz1eLcFLPV4iErkaGhq48cYbOeKII6isvApjoLLyp3zlK1/hjjvuAOqcHq+jSUzs37UzDTUUEREJPwWvMLCD1x6amprCfm4RiS4ffPABPp+PgoICAoMSLMvitttuY9u2bcDv2bBhA33dv6s19XiJiIiEn4JXGNjBq57Kysqwn1tEosubb74JwEknndTm9m984xvMmjWLBBZQVVXFyJH9D16a4yUiIhJ+cRi8RoX93IG9vI44QsMNRaR7b7zxBjk5OUycOLHN7YFeryb2AfDKK/1bWAM01FBERMQNcRa8RgDJYT93dnY2ACUlCl4i0jVjDG+88QZf+9rXOv36WWedxQxg8ODBHH744f0+j4YaiohI7xnuueceTj/9dI477jjgJNauXet2UVEpzoLXGDwe8HjCe+5AjxdoZUMR6VpxcTG7du3ihBNO6PTrlmXxBLBixQqSkpL6fR4NNRQRkda6Xxr+Pm655RZ2797NsGHDgK2ceOKJwO/DWGFsiLPl5Md02IMgHFqCl3q8RKRrgfldXfV4AUwFpp533oDOo6GGIiLSWldLw69cuRK4mQsvvJCnn36ahIQELGsvJ5wwh1df/QFXX/0ucB8wpM3jjDGUlpYyYsQIkpPDP9osUsVdj5cbhg8fTgqgHi8R6c4bb7xBeno6X/7yl0N6HnuoYS0NDQ0hPY+IiESv999/n8suuwz4Gn/84x9JSAjEhlG89NJLwM088sgjwPF8/vnnADQ2NvKXv/yFY489lqysLAYPHszYsWO5//773XkSEUbBKwwsy8Ke5bXLlfOLSHR44403OO644xg0aFBIz2P3eMH+/er1EhGRzhUUFDgX6p5jyJCWHi2PB5KSEvF4/pf/+7//A4qZMmUKQ4YMITU1ldmzZ+P3+1m0aBG33347U6dO5cYbb2TNmjVuPZWIERdDDRsaGigvL8et4AUwDihSj5eIdKGiooKPPvqISy65JOTnCgSv6upq0tPTQ34+ERGJLmvWrOGll17innvu4eabR7b5WttpO+cA7/E///MYtbW1GGM49thj+fa3v918EfGmm25i5syZfOc73+H9999n1KjwrzAeKeIieNmhC9wMXs66hq6dX0Qi23/+8x+g+/ldwWJfwbSDl4iIxK/cXHt+V+uF54wxLFy4kPHjx3P11Vdz8809HSWHO+64o5vjp5Gd/TTl5V/liiuu4G9/+xtW+8lkcSIuhhq6uXlygIYaikh33n77bQYNGsQxxxzTvLpUqFZg1VBDERGBlkU1Ar1YHg8kJPydt99+m9tuu42UlJSgHL+kZDqLFi3i73//uzM/LD4peIXJOACqdYVZRDr12WefMWnSJFJTUzv8IQy21kMNRUQkDgWu8IH9r/Ox1WcxnfOAQ/jB3Lmd3qfDR+Drzlr0gUO3voDo8cBNN13LoEET+OUvfxnGJxpZFLzCJNv5t6REww1FpKMtW7aQl5cXlnNpqKGISJwLXOED+1/n443XXmMDALeQGLi93X06fAS+7vO1OXTrC4jFxWBMMo2Nt/D666/z+uuvh++5RpC4CF579uxxPnM/eJ14ooKXiHRUVFTE5MmTw3IuewNMqKqqCsv5REQkOjzwwAOMGDEC+E5Ijj9x4o+A0Zx++l0hOX6ki83g1bqP07LY9ZOfMBiAEW27RkM1gaITgeC1Z4/meYlIW5WVlVRWVoYteA0fPhxQ8BIRkRbbt2/n2Wef5corrwRSQ3KObdtSWbToRmpqXsKy1gVGJ8aN2Axerfs4jeE3qZfSmJgHWG27RkM1gaIT45o/U4+XiLS1ZcsWAG6+eXJYrgkFgtcXX3wR2hOJiEjUeOSRRzDGcPXVV/fvAB5Pr+aEXb1wIRnANzjGHp3Y3fyxVnPHYkFsBq92DhzYxaBB48LZwdXBCGDw4MEMH14SSz8/IhIEgeD1wQeTw3JNKLC4hoKXiIgAHDx4kGXLlnHeeefh6eMb5kDeyqW4V3PChhvDrffdx4sAvND9/LFWc8diQVwEL9jFeeeNC2cHVwcWkJWVxbnnlsTSz4+IBEEgeE2aNCks50tISACGsW/fvrCcT0REIttTTz1FeXk51113HdASpnozCsNeOKNv+eiaa65hypQpwE+oq6vrd93RJuaDlzEG2Mn48ePdLoVx48axa5fmeIlIW0VFRcDY5p6o8BiuHi8REcEYw4MPPsi0adM4+eSTgZYw1ZdRGIGw1psOs+TkZH71q18Bn/HQQw/1s/LoE/PBy548vp9x48b1eN9Qy87O1nLyItKB3eMVnoU1Wih4iYgIvPnmm7z33ntcd911WIE5Wv0QCGu9DWpnnnkm8E1+/vOfx83745gPXoEeJgUvEYlU7gSvdAUvERHhwQcfJCMjg+98JzRLyHfvfmpra7n22mtdOHf4xXzw2rlzJ0BEDDXMzs7G7/cDB90uRUQiRG1tLTt27ADCs3lyC/V4iYjITv7yl7/wwx/+kKFDh7pw/in8/Oc/569//St/+ctfmneEitWF6GI+eEVSj1dLDer1EhFbcXGxMxdVQw1FRCTcHqWpqYlrrrnGtQpuuukmZsyYwbXXXovPVxlrCxm2EfPB6/rrIyd4ZWcHtlFW8BIRW2BFQzeCl1Y1FBGJX/v37wce4dxzzw3bqrrteTyQlJRIScnjlJWVAbe7Uke4xHzw8vt3kp6e7lL3aVsKXiLSnpvBSz1eIiLx67HHHgPKKSgocK2GwIIcJSVH873vfQ9Yzt69e12rJ9RiPnjBrojo7YLWvW5aUl5EbD/96RZgKDk5Y8J85uFUVVXR1NQU5vOKiIjb6oB7770XOInjjjvO7XIAuOWWW4CDMb28fFwEr0hYWANg5MiRJCYmoh4vEQmoqtrCEUfk4fP1fwnf/hkOQHV1dZjPKyIibnsS2L59O/BfbpfSbOrUqcC5TvDa73Y5IREHwWtnxPR4JSQkkJWVhYKXiLQoYvLkcA8zhMzMdAAOP1zDDUVE4oW9amATVzKZo446CviG2yW1U0B5eTnwW7cLCYmYDl72EJqSiAleEJjnpeAlIoE2yp3g9cgjdo/Xzp0KXiIi8cLng7/97R80sMWZ2xXu0RY9OZ4TTjgBuJeGhga3iwm6mA5edmKuj5ihhmDP80pK2hWz+xOISO+VlpYCNeS60CAMHz7c+UwrG4qIxJNf//rXwEQuvPBCt0vp1I033gj4ePnll90uJehiOngFNk+OpB6vcePGkZa2M2b3JxCR3vM5DYHH4wn7uVuCl3q8RETixwe8+uqrwLXOugOR5+yzzyYhYTRnnfVbezNltrpdUtDEdPCKpM2TA3JycqisrCRWJw2KSO8peImISHg9QEpKChM4Hcuy99GKBB4PzfUkJydzww2XkZT0f5SWluEj1+3ygiYuglckDTWcOHGi89l2V+sQEfdt27YNUPASEZHQshfVKAOe4PLLL2c70zHG3kcrEgT28wrU84Mf/ID6+nr+/Oc/u1lW0MV08AoMNbRXEowMOTk5zmfbXK1DRNxn93ilk56eHvZzK3iJiMQPnw/uvHMZUMv111/vdjk9mjZtGsceeyy//e1vAeN2OUET08HL7vEaQ1JSktulNGvp8VLwEol3dvDK6fF+oTBs2DDnMy2uISIS+wy//e1vyc/Pd/bLinxXXHEFGzduBNa7XUrQxEHwipxhhmAPe7QsCw01FBE7eLkzwH7QoEGkpaWhHq/48umnn3LxxRdz8OBBt0sRkbDaQFFRERdffLHbhfTaxRdfTEpKCvCo26UETUwHL3uoYeQsrAGQlJTk7OWlHi+ReGfP8XJvZrM93FDBK578+te/ZsWKFXi9XrdLEZGwWsmgQYP41re+5XYhvZaens5ll10GPEFZWZnb5QSFgpcL7Hle6vESiWdffPEFfr8fBS8Jl/r6epYuXQnAeed96nI1IhIuxhhgJaeeeiojR450u5w+mT9/PlDL0qVL3S4lKGI2eNXU1Dibk7ozf6I7dvBSj5dIPPM1b+bnXhul4BVfCgsLaWraC0B5+WcuVyMi4fLee+8BWyJ2w+TuHH744QzhJH72s9/g8dS5Xc6AxWzw2r7d7lEaOTJCNihoxV5gY7tzBUJE4lFL8FKPl4THU089BWQwefJkQMFLJF6sXLkSiK5hhq09x2vAbrZte8btUgYsZoNX4E3NX/4SqT1eNezdu9ftUkTEJYE9vNwPXlrVMB4cOHCAZ599FpjNtGnTAA01FIkHxhgneH2dUaNGuV1Ov5wOzkqMv4r6TouYDV5ubkzaqcCW3JbFxPnzAdg+ZkzzbViWvbudiMQFn89HcnIyMNa1Guz9w9TjFQ+ef/55qqurgUuYMmUK8DmNjY1ulyUiIZSbCwkJG/n888/JzIy+YYYBFnDLLbcAXp55Jrp7vWI2eNk9XgmMHx8hy8kHtuQ2hpx16wDY9uyzzbdhjL27nYjEBZ/P5ww7dq8Z1lDD+PHUU085K+r+Pw499FCgrtVwVxGJRT4fLFnyEgAffni2y9UMzPe+9z3gaG655RZycg5EbX9FjAevcRG1eXJAYBPlK67QyoYi8crn87neI28HryqaXK1CQi03F5599nVKSr6JxzPI6fGCzz7TPC+RWPfyyy8zbdo058JL9Bo0aBDwANu3b2f79iVR218Rs8HL7f1xujN69GgGDx5MZaVWNhSJV9u2bYuQ4GXY72oVEmo+Xzmwl3vuOZziYpweLwUvkdh3kNdff53TTjvN7UKC5ETmzJkDLG41Tzq6xGzwsnu8Im9hDQDLspxer+j8oRGRgamrq6OkpCRCgpcGG8Y+eyGNww47DLAv/kE6n36qBTZEYtsb1NbWxkzw8nhgxYolWBb813/9l9vl9EtMBq8mAsvJR2aPF2gTZZF4tn27vZ3E7bfn4Gb2CgQvrWsY6+yAFejpsiwLOFQ9XiIxbzXJycmcdNJJbhcSFPZyCTksXHijsz3GOrdL6rOYDF67gfr6eiK1xwvs4DVo0LaonBgoIgMTWNTglVc8FBe7V4d6vOLFJyQlJZHb5g/OFPV4icS8lzn++OMZOnSo24UE1YIFC5yl8W+JuuXlYzJ4tcy1i9wer4kTJ2JMCT5fvduliEiYRcp2F/Zy8gpe0Sw3tze7kXzKIYccQmJiYqvbDmX79u0cOHAgpPWJiDv27NkDvM/pp5/udilBl56ezm233Qas4YUXXnC7nD6JyeDVMnMqsnu8mpqagF1ulyIiYVZcXAxYTJgwwdU61OMV/Xy+3uxG8mnzMMMW9sqGmzdvDlVpIuKiV155BSBm5ne1N2/ePOAQFixY4Lyfjg4xGbyipcfLpnleIvHmk08+ASYxePBgV+tQ8Ip99rD7z5sX1mihJeVFYtnLL78MZDJ9+nS3SwkJe7uon/Pxxx9HVa9XTAavbUBCwgg8nmFul9Ile3ENaB0TRSTGBMaBtfv4eMUKYGqnX+vxI4jDExW8Yt/WrVuBhk56vA4BFLxEYlFTUxMvvfQScJqz/1Wsms2ECRO477773C6k12IyePmApqYcVyet96RlkvMWN8sQkVAKjANr9dHY0MBngwcDh3f4Wq8+gtiwDRtmX5zSqoaxK7CARsfgNZTx48drqKFIDHrvvffYvXs38E23SwmxJK6//nr+9a9/8d5777ldTK/EZPBypq27XEX3UlJSnOGGn7tdioiE0datW6mtrQWmul2Ks9hCqnq8otp7FBUVdflVe1hrZ8ELxo0b57w5E5FYcuqp/wRgwoRvuFxJ6F155ZUMHTqU+++/3+1SeiUmg5c9eC9yF9YI+NKXvgToaqNIPPn444+dz9wPXrbhCl5RqKamhptvvhmYwcUXX9zl/ewerzGMGDGize0eD6xdO5ZXX90T2kJFJOz8/heYOXMm27ePdbuUkMvIyOAHP/gBTz75JCUlJW6X06OYC1779u1zhs1Edo8XwCGHHEJCwmbt5SUSRzZt2uR8FinBK13BK8pUVFQwfPgs7r33XpKSvszatWvpar6wHbzaL6xhj1j90Y+yqK9Xj5dILKmoqADe5swzz3S7lLCZP38+DQ0NPPTQQ26X0qOYC16B/XGiJXg1NZXj81W6XYqIhMnHH3/MuHHjgHS3S3GoxyvavPjii9TXf8SKFSvYtOk559a/dnpfe6hhx2GGAGPHjgXKomopZhHpnr2aYRPf/Gasz+9q8aUvfYnzzz+fhx9+mC++iOy/aDEXvHzNm5lE/lDDQw45xPlM87xE4sWmTZuYOjVSersAhmtxjShjL4hhce655zJ58mSOOuooOgteFRUV7N27l+6DVyPl5eWhLFdEwiQ3Fy655J8kJGRy7LHHul1OyHk8LRvIFxQU4Pf7WbZsmdtldSvmgldx84pfkd/jZc/xAs3zEokPxhg2bdrE4Ycf7nYprWSgPvfoYi8Bn8OQIUMAOP/884E3OyyUEVjRsLOhhhAIXrBnj+Z5icQCn6+J0aP/yZw5Z8TWMvKBhNXuo9hnYbDw+eCYY4/l68CvbrmF2t5szeLSPJ+YC16ff/45FkPJyYn8CYWTJ0/GsiwUvETiw44dO6iuro6wHq8xlLldgvSOsy/cZ08+CUxpfgNxwW23AYbs7Gft9xNWMVgW/z7+eOeBR3b6xmPsnDmAgpdIROtiP8jWH7lWMZYFY1lJWVkZZz71VFj3hAy54uJut1nxeMDC8OGYl9kFPPHYYz1vzdI8Qi68Yi54bd68GcMh+HyW26X0aMiQIVpSXiSOzJplL6xx9dWHR9DfvDGUAw0NDW4XIj3x+TBNTWxOTwemNL+BOLypiSlTppCf/1f7/QS5mKYm/nT44ZxwwgnAxE7feIx1FnpR8BKJYJ3sB9n+w0cuxsDx315BZmYm366uDuuekG4L5LLS0nxmzJjB3XffTV1dndtldSomgxcc0uP9IoWWlBeJH6Wl9lLye/ZMjaC/eWMAnLlAEunKysrYt28frf/OWZbFBRdcwL/+9S/Kyuz+S6/Xy8cff8xll13W5bE01FAkdmzZsoXnnnuOq666iqFDh7pdjkssfvGLX/D5559z7733ul1Mp2IqeDU0NLB161aiKXjZC2woeInEh01kZmYyevRotwtpxQ5epaWlLtchvWFfXASY0ub2QMD66U9/CsCf/vQnIJmrrrqoy97VjIwMIFmbKIvEgAceeIDExESuueYat0tx1Zlnnsn555/PHXfc0e3m8m6JqeBVXFzsDJf5Uo/3jRR28Kpw9l0QkdhmL6xhz+2MFApe0cReWAPaB6+pU6dyww03sHz5cuA1nnrqKeAcjBnRZe+q/XM4Rj1eIlHPz+OPP84ll1zibFcSnwJrcLzzzgMkJSVxzTXXYIxxu6w2Yip4ff65PVdqLGkuV9J7gSXlA7WLSGyyG/+PI2xhDVDwii6fffYZiYmJdLZy7+23387EiRNJSDiP0tJSRo/uephhiywFL5Go9wj79+/nJz/5iduFuCow12vnzvHceeedvPjii6xatcrtstqIqeB16aX2EIz3ucjlSnovsKR8y/AREYlFr7/+OlDOiSee6HYp7Sh4RZPNmzczefJkILHD19LS0njwwQdpavIzcuRIduzozQaqYxW8RKJMYKHD3FxYvXo1cCvnnHMORx99tLuFRZBrrrmG6dOnc9NNN7F//363y2kWU8GrsnIzw4YNc95GRIe8vDzAUo+XSIx77LHHgOFccMEFbpfSTgaJKHhFi88++4wpU6Z0+fXzzjuP+fPn87Of/Yzk5OReHFHBSyTaBBY69Pm8zj5+hzvzOgXsIYeJiYPYsOFBtm/fzt133+12Sc1iKnjBZr70pS8RSbMnemJvgJmjHi+RGOb3+1m5ciXwHVJTU90up50ERqOV7aJBE3aPV2CIelfuv/9+5s+f38ujjqW0tJSmpqYB1yci4WMvHHEmmZmZwD9JT093u6SI0bLt1wnAd/jf//3fiFloI+aCV09/kCLTlxS8RGLYk08+SU1NDfAjt0vp1BjU4xUNdgI1NTXd9nj13VgaGhqorKwM4jFFJLRKOeOMM4B6XnzxRSB+F9Toyfjxi6mtTWTatJvcLgWIoeBVX18PFEdl8Bo27DDWrv2EyFp3RUSC5bHHHmP69OnADLdL6cDjgfc5jX/8o5TcXLerke4E1jMM7t85ey8vLSkvEh2qq6uBs9i5cydZWf/g8MOndrllhMCOHeP5+c8XcPDgc3z66adulxM7wau4uBhojMrgtWTJNIz5gm1uFyIiQef1etmwYQM/+lFk9nYVF8N3WU1ubik+n9vVSHead/AKco8XaKipSLS48sorgQ0888wzlJQchzF0uWWE2ObNmwckcthhy1y/wBgzwSswVC8ag9cRRxwBwIcu1yEiwffggw+SkpLCpZde6nYpXdJQw+jwGZCamhrUfXqyaQDg1FP3uP6GRES69/LLL/P0008Dt3L22We7XU7UGDt2LBdddD6Zmb/H5zvoai0xF7wCy7NHk2nTpgGw0eU6RCS4JvAmv//9Exw8+CNGjMiI2OEgY4ADBw4AkbPkrnT0GfbfuISE4P3p3sjpAPzqV3vU4ykSwWpra7n22mud97kL3C4n6lx11VVUVFQA7u7rFWPBazijR492u5Q+s1eiyVHwEokxO3mGxEQoLr4pooeDtGzBoV6vSLaF4F9cHAEkJSVpqKFIhLvnnnvYvHkzDz30EDDE7XKizsknn+wM037U1TpiLHgdgmVF02LyrU1T8BKJIXv37gWWc+mll+KJ1K4uh4JX5DPGUAxMmjQpqMdNAMaMGaPgJRLBtm3bxl133cUFF1zgrGYofWVZFldddRXwFh988IFrdcRc8IpeR/AJgdUZRSTaPfjgg8ABCgoK3C6lRwpekW/37t3UAPfem4dlEdRhq2PHahNlkUh266230tjYyNtv3xv03/94cvnllwPJ/Pa3v3WthpgIXvv27WPr1q3ANLdLGYAjqAc8ns96vKeIRLampiYefvhh4FwOP/xwt8vpkYJX5Ats/vnCC5OCPmxVwUskcn0A/PGPf+S6665j505PRA9bj3T2ZtPn8uSTT+JWN0dMBK/333/f+Wy6q3UMjL2yYUmJBhyKRLuioiJnqOG5bpfSKy0zYxW8IpV9cTH4Qw3BDl7ax0skMv0X9loA//3f/+12KTHiMsrKynjJpbPHRPDasGGD81nkbU7ae4cCg9DahiLRr6VNio6LQalAWloaCl5BlJsLlhW0j62XXeYcNjfopY4dO9bZTsAE/dgi0n//+te/+CewcOFCp7dGBion5xvAKOZwpivnj5ngNWhQFh5PltulDMBg7PCl4CUS7bxeL0lJScCX3S6l18aMGYOCVxD5fGBM0D6Kvv99YBxDhgR/NbOxY8c684srg35sEemfxsZGbrrpJiYA1113ndvlxAyfL5nrrruEA7yC3+8P+/ljJng1Nk6PgTGvR6DgJRL9vF6vsz/fYLdL6TUFr8iVmwu///1WBg8O/jBDgPHjxzuf7QjJ8UWk78aMeZwNGzZQwwOkpKS4XU5Mueyyy4BaVq5cGfZzR33wqqmp4aOPPiJahvR07wigmKqqKrcLEZF+MsawYcMGpk+PrjZp7NixxG3wCvKwwGAvO+bzwcSJRVx0UV7Qjtlay7yxrSE5voj0TWVlJRUV/82JJ57IXq51u5yYM2vWLOAw/vCHP4T93FEfvD788EMaGxuJneCFEyRFJBrt3LmTsrIyZsyIrjmncd3jFeRhgcFfdqyOHTt2hGRhDVDwEok0t912G1DJAw88gAdfKK7nxDV7z98rePPNN3nvvffCeu6oD17RNom9e0cCrZ+TiESbwO9vtPV42cGrjKamJrdLkQ62YYwJWfAaOXIkQ4cORcFLxH3vvvuusx3JPI4++miKmRSa6zlx70ekpqbyta/9GsuyBz6EQ0wEr+HDhwOh+YMUXh5gDO+8847bhYhIP3m9XizL4qijjnK7lD6xg1cjFRUVbpciHdh7eIUqeFmW5Ry7OCTHF5He2b9/Pyec8F0aG8cxYcJdbpcT0zzs48CB77N//5OUlOzG5wvPeWMieNXWHo3HE91PxeMBDz7gq7z99ttulyMi/eT1ejnssMOcHoToYQcvnGXFJbLYPVF5eaGZ4wWBUKceLxE33XzzzTQ0fM6rr/6B7dtHuF1OTCtmEp9+Oh+o49FHHw3beaMvrbSaBN1oWXzw9tvU1k6n2BeaSc3hUlxs/xDAV/n000911VkkSkXjwhrQErz27NnjciXS0VaSkpIYN25cyM4QCF7GaC8vETc8//zzTgC4kVNOOcXtcuLClClTOPvss52hnTVhOWf0Ba9Wk6A/+/hjDgAwPYSTmsPtq4A9xldEosvevXvZvn171C2sAZCTkwPA1q3q9Yg8RXg8HgYNGhSyM9jBq5ry8vKQnUNEOrdz506+//3vc+SRRwIaYhhON9xwA2VlZcCfwnK+6Ateraxfv975LPquLndtFpZlaZ6XSBQKLKxx883T7c73KJozYw9jS+WDDz5wu5SoZozhgw8+CPIiJVtDOswQWuaPKXiLhFdjYyPf/e53OXDgACtWrCCa9n+MBV//+tf5yle+AtxJbW1tyM8X1cHr1VdfJSFhBDk5X3a7lCAaRmLiNO6+W/O8RKJNIHiVlx9td75H0aI/CQkJwBEKXgPQ2NjI/PnzOeqoo7j//vuDeOStIVtYIyDXWdJLwUskvEaOvIs1a9aQkvIbDjvsMLfLiTuWZXHHHXcA23j88cdDfr6oDV7GGF5++WWamvLx+UI3/MIN9fVfpabmHS3rLBJlNm7cCEwgMzPT7VL66Ug++OADzfPph4MHD3LRRRfx4IMPApksWHAfdXV1Az7uF198AZSHPHgFjl8c1UP1RaLLO++8w759P+e73/0uQ4deHq3LFES90047Dfgad911FwcPHgzpuaI2eG3atImdO3cCp7tdSlB5PDBy5FeBSjZv3ux2OSLSB3bwOsLtMgbgSMrLyykpKXG7kKgzf/58nn32WeBXvPDCEzQ27uTpp58e8HEDPVChHmpob8uSqR4vkTDxeGr46ld/wKBB4/nNb36Dz2dF/zIFUcreUPlOdu3aFfIVDqM2eK1evdr57DRX6wi24mJ47TV7gQ3N8xKJHvX19WzatIloDl5jx9qbuE+fruGGfVFZWcmf/vQn5s6dC9zAN77xDZKSpnH55fdgWWZAG3MGhq9OnTo1KLV24PG0rAjMJLY++mjL//v7Ea6dSEWi2LZtPwc28fzzy50LH+Ku/8epp57KokWLnJEGoRG1wevll19mypQp2JsOxxZ7jO9w7eclEkU2b97sDC2L3uC1aZNde2np+y5XEl3++Mc/UlNTw9KlVzk5xmL58puBjfzzny8NaGPO1157Dcjk8MMPD1a5bRUXt6wIzCS2TpnSdpXg/nyEaydSkSi1du1aYAk//OEPOeOMM9wuR7CvQb3yyt2UlZWxePHikJ0nKoNXbW0ta9as4fTTY2uYYYA9yf1Y3njjDbdLEZFe+n//byMA2dnRG7xGjBjBxIkTAfV49ZYxhqVLlwJfwZijm4cJXXLJJYwfP5777rtvQMd//fXXga85fxdCbRLFxcWaXywSQlVVVXznO98BxnHvvfe6XY447GtQs4BLue+++9ixY0dIzhOVweutt97iwIEDMRu8bKezceNGJkwo1qgNkSiwd+9GBg0axNat0b0qlb2PjIJXb73xxhvOENN5bW5PTk7m+9//Pq+++ipQ0a9jl5SU8PnnnwMnDbjO3sjMzKWuro6cHM3xE+mX3Nweh+JeO3w4WzZvBv5MekZGz0N3tdpGWI0ffxc1NYaJE/8nJO+/ozJ4rV69msTERK655uTY+nlsM9b+2wDs3Pk3e9SGxs6LRLiNTJkyhcGDo3sPFjt4fRKW/Uyiye7du1m1ahUVFW1D1KOPPkp6ejowp8NjzjnnHBobG4GX+nVOu7cLwhW8nnjCXtlw587isJxPJOb4fN0Ow33iT3/ij8DPbrsNOKl3Q3e12kZY7diRy4IF87GsP+Lzre/5AX0UlcHrnnv+yaBBx7F9+7DY+nlsM9b+SyQlTWPw4Gftr2nsvEiE28gRR0TvMMMAO3g18Mknn7hdSkRYv349p5xyCuPGjePCCy/khBNOcFbUhRdffJFVq1bxve99D0jt8NhjjjmG0aNHA//o17lfe+01hg4dCkzv/xPog5Yl67WyoUiwvf3221x11VWceOKJ/M///I/b5Ug3Fi5cSFZWFnBZ0JeXj7rg9R5QX/8eqakXxFZvVzseD6Smfov6+teBMrfLEZFuVFVVAVtjKHihjZSx528df/yVrFnzMcbcyqhRT7Fz505yc0/Esm7gm9/8JoceeigFBQWdPj4hIYEzzzwTeJGGhoY+n//111/n+OOPBxIH9kR6Kbd5RIWCl0gwrV27ljPOOIPs7GyKip4mKSkxpt/DRoXWo8zafWSMGMHvS0qATSxITe15OGgfRF3wWgbAELZs+V5s9Xa1U1wMr776bWeSc/+ulopIeHz00UcAMRG87NViByt4AWvWrKGubgPLl9+FMbezd+/FvPrqqzQ07AN+DczlnXfeYfz48V0e46yzzgIq+rxKbWVlJRs3buTEE08c0HPoiyFDhjiLq2wM2zlFYt3atWs5/fTTOXBgJJ9//iqJieM0gjAStB5l1snH6cYA83kI+OcLL3Q/HLQPoip4VVdX8wQAFzFixAiXqwm96dOn4/F4gGfdLkVEuvHhhx8CURy8Wl35S0xKAr7Mhnvuifv9nO655x5gDN/97ncB+9t0zDGzGDfuXVavXo3Hs5TU1JRu57/bi0Al8o9/9O0C2ptvvokxhltvPSmsV8btpa1f7Nccv7Vr1/Lyyy+zEdi/f3/QaxOJNs899xwnn3wyGRkZNDT8C2MmKnBFlbuZNm0al156Ke+++25QjhhVwWvFihVUATDX5UrCw7IsvvWtbwEvU11d7XY5ItKFjRs3AkNbDdWKMu2v/HEK/05MpKy0NG73c9q0aRMvvPACcA1DhgwBWr5NO3dOJj8/v823ras3U+np6QwefCKLFz/fpxxq79+VxIEDx4b1jdq3v/1toMpZjbH3PvjgA4477jjOOOMMjgQmT57Mnj17QlKjSKQzxnDPPfdw/vnnM23aNP7zn/8Qi/vOxjqPZwgffvh3qqszOfXUU/n3v/894GNGVfC6+uqlJDEFON7tUsLG/iNYy8qVK90uRUS6YAevL4dpr6Vw+D4NDQ38+c9/drsQ19x3331O4PrxgI/1y1+eDXyIr5dBtLq6mqeffho4jpSUlAGfvy++/vWvA2k8+2zvR1o0NTXx4x//GGNGAKuB37NnTzlf+tKtoSpTJGI1NDRw9dVXc8sttzB79mzWrFnjLNQg0ca+uJbL2LGvU109kVNO+Qa/+93vMH0cXtha1LxLWL9+PXV1a6nnGqBvE9mimT2+fwZ33HEHdXV1bpcjIu00Njby/vvvA1E6zLATHs80YBYLFvze7VJcsWfPHh577E/U1FyOxzN6wMc755xznM96F2R/9rOfsX37duCXAz53X9lh80z+9re/OUvh9+wPf/gDb731Fk1NSzAmH8P3mT//GqqrH3MuSojEh2rgvPPO49FHH6WgoICnn3467BdPJPh27BhHaem/MeYrXHHFFQwd+i0sa0+/RtNHRfAyxnDzzTcDmcD33C4nrOwr6L+kuLiY5cuXu12OiLTz5ptvOns7neZ2KUFTXAwPPfR96uvf57333nO7nLD79a9/DdTx6ac3BmWY3yGHHMKpp54KPEx9fX239127di0PPPAAP/7xj4ETBn7yfvk2paWlvVoQpKKiggULFnDCCScAlzfffuuttwLp3HjjjQO6OiwSLXbt2sVIvswLL7xEZuaj3H333TE0CkJGjx5NTs6rwL0cPPgSI0ZMxef7TZ9XrI34n4jcXBgzZgVr1qwBfokHfxwuwXk6J510Etdf/wty+MjtYkSklbPPXgUMYeLEM90uJaguueQSIJnf/e53bpcSVlVVVTz88MPA+c4Kj8Exf/58YCd//etfu7xPXV0dP/rRj8jKymLRokVBO3ffnUlycnKPww2NMVx77bVUVlY637OWtxSZmZnA7RQWFvZ5YRGRaLNx40a++tWvUkcxzz//fwwbNq/DWkPx99419vh8CRhzIx9/vIHp06cD1zJr1qw+HSPig5fPV01l5U3ADHJyfkQxk+JwRRiLRYsW0dS0h+38ze1iRMTR1NREVdVfOe+8M9i2bZjb5QSV/cb5PP785z/H1TDnZcuWsW/fPqDzvbn6y15WfjIPPPBAp1+vqanh/PPP54MPPmDXrofJyEh37Y2axzOcurpT+fWvn3W2NOncH/7wB5566ikaG3/OUUcd2VKvs0pmDt8EDufcc68hx/pwYCtkikQgj8dgWcs46qjjnZ6P1znzzDM7Xak8/t67xq6pU6dSWFjIqFEref99f58eG/HBC26nsXEX//nPb/D5BrldjCs8Hrj00uOB84Cf884777hdkoiAs7zsTmbPnu12KSFyBeXl5Tz22GMhO0NjYyPvvPMOL730Es899xy7d+8O2bl6UldXx69+9StOOeUU4JigHtsecnQtb731FuvWrWu+PTcXLKualJSzeP75F8jMfBRjznP1jVpxMTzxxHdoaCjinHPOobKyssN9Pv30U6655hrgFBoa/qttvc67Tp85hP/853EsawfbeXRgK2SKRJiPPvqIbdu+CczDmGOd92bT3S5LwsSyLMrKZlNf/3nfHmiM6dfHzJkzTagtW7bMAGbu3LktN0LIzxuJ7L88ZQbyzNixY822bdvcLkliELDO9LNNiJSPcLRNATfddJOBJFNZWdn5HaK8vYImk5+fb4YNG2a2b9/e9wN4PN2+nd4P5jwwtPoYA+btgb1F7/uHx2OMMebhhx82gHnxxRdD8tKB36SlpZkhQ84wUGrAmOxsrznyyCNNQkKC+eMf/xj8k/ZTU1OTgYdNUlKSmTRpknnsscfMW2+9Zbxer/nFL35hBg3KNTDSjB+/o8dj3XDDDc7r+1rgW91napvEbXV1debdd9819957r5k5c6bzM51qHnroIZOT02g3JWx1u0xxQV/ap4hpQDwe06ZBfvrpp41lWQa+YWpra1s/u6CeN1oEvj/wkRk+fLg56qijTGVlZYfvm8hA6M1N7zU1NZnc3FwD3+z6TlHeXtltzhYDqSYl5RzT1NQUtGOXlpaar3zlK8ayLJORsdjAmwZeMTDJwBAzatTKoJ2rN1566SWTlJRkBg8+2UBTSNpV+/v5vwYSzNChQ81FF11kEhMTTVZWlnn++eeDf8IBAmPeeustM3HixDbh2P74ilmzZk2vjlNdXW1yc3PNpEmTDHzUz1rUNkl4NTU1mXfffdfcdddd5tRTTzUpKSnNP/9JSUcbuN9MmFDa9kFR3uZL/0Rl8ApceCwvLzc//elPTWJiovna175mYH/HO8YxD1vNmDH/NImJic6bvrfj/VsiQaQ3N723bt0654/w413fKUZ+Oe+9914DmCeffDIox/vwww9NXl6eGTJkiHn22WfbfJtKS0vNcccdZwBz9913BzXsdeXtt982Q4cONUceeaSBypCfb9OmTWbOnDnGsizzve99z5SXl4f8nP0ReF3q6+vN5s2bzd///nczcuQTBkr6HEzfeustM3r0aANDzCOPPNLn11Vtk4RDU1OTycp6x8CNZtCgnOagdeSRR5rrr7/erFixwmzbtq3rpj1G2nzpm6gJXvX19Wbr1q1m9erVBh4ycJUZ7vyQXwRmIu/b3badDAmJW9B8FdLj8RhINHCzWb9+fVjeoEhs05ub3hk37jMDeQbSzIQJe7u+Y4z8EW5oaDDJyV8xkGhGjLhnQG3N3/72N5OWlmaysrJMVtbbnTbrBw8eNHCxAcwVV1zRdtRDENTW1pqdO3eaf/zjH+a73/2uSU1NNYmJk/sVKAaisbExfCfrh85Gig7k+1NSUmLgdAOYo48+2ixfvtzs37+/28c0NTWZ0tJStU0SXK1+uA+AWQ1mAZi8QI8WmBS+buAPzcOC2/wetH9vqveoca0v7ZNl37/vZs2aZVpPEO5ObW0t27dv57PPPmPTpk3cdtsG6uq8GLO53fr3w0hNPYMDB36Gx3MkPp/m1HaQm0uubw0+cgE/qXyfA/wdaAImMoRJXJ7mJWfhQnJycpg4cSLjx49n/Pjx2sRPemRZ1npjTN/WRo0wfWmb+urgwYO89dZb5OdfzKhR8Pzzz3Psscd2/QDLiplGzO/388Mf/tBZDv0sRo+ex6ZNxzNy5MhuH2eMobKykhdeeIE//elPrF69mpkzZ/Lcc88xYcL4Lr89Ho9h27bbgTuAiaSnX86bb17MpEmTSE1N7fac9fX17Nixg61bt7Jp0yY+/PBDNm/ezK5du9i1a5ezaqFtxIgRnH/++Tz++K0Yk9OXb4n0g2U18fjjv+eqq+6nvn4jkMqQISdx5535ZGVlkZSURGVlJe+//z7vv/8+b7+9iaamSgC1TXGioaGB/fv3U1NTQ21tLUlJSSQnJ5OSkkJKSgpWL1e5bGpq4osvvqCyspKKigoqKyspLy+ntLSUPddfz5aLL+aTTz5h06ZN1NbWkpiYyCmnnMLq1ZdQWfltMjIyQvtEJWb05b1Tv4NXamqqOfTQQwk8PpDkmpqaqK+vZ8uWOqCW5GQ/Bw8ebPfo8cBMhjOBRHKo4Fh28nWOS9iGNXEixcU07wat5Td7tnfvXv7+97/z4osv8swzWxjFevZ2cr/U1FSGDx9OWloaSUlJJCYmNm/uZ1lWrxsziV0bNmyI+jc3qamp5rDDDhvQMVq3a01NTTQ0NFBdXc2OHTucr03i009f7HmfpxgKXmB/Px566CEWLFhATU0NYG8qmZaWRmpqKgkJCViWRWNjI/X19dTU1FBaWtp839zcXC6//HIKCgqcN1A9f3teeOEFHnzwQV588WXsC0yQlpbWph2zLAtjDDU1Nezfv5+qqira/m1LBw4DJgDZpKePZfHi0dx6ay6lpacAyXg8+nsTDrm54PNBTo7hiSfeYMWKFSxd+goNDZ+0u+dw4CjS0r7ML35xKD/5yU/UNsWA9m1r4D1jTU0NBw8epKqqqrm96Nwg0tPTmkNYcnJy8/uYxsZGGhoaOHjwINXV1ezfv7/L7RASgITESTQ0HAZ8GTgFOBEYprZA+iwswcuyrCrg0349OHqNgk4zTayLx+cdj88Z4FBjTFRvSBWnbRPE589sPD5niM/nrbYpesXjz2s8PmeI3+fd6/YpcQAn+TTarz71lWVZ6+LtOUN8Pu94fM5gP2+3awiCuGubID5/ZuPxOUN8Pm+1TdErXn9e4+05Q3w/797eNwo2UBYREREREYluCl4iIiIiIiIhNpDgtSxoVUSPeHzOEJ/POx6fM8TG846F59Af8fi84/E5Q3w+71h4zrHwHPojHp93PD5n0PPuUb8X1xAREREREZHe0VBDERERERGREFPwEhERERERCbEBBy/Lsioty1ptWdaCYBQUqSzLmm1ZVr5lWXPdriVc4uW1hebXd3Unt8X0a97F846Z1z2WnktX4uHntDPx8NoGxGP7pLYp+sX6z2hX4uG1DYjHtgkG1j4Fo8frQmPMacaYJUE4VkSyLGs2gDGm0Pl/vrsVhU3Mv7YBxphVrf8fL695++ftiKXXPZaeSwfx8nPahZh+bVuLx/ZJbVN0i4ef0W7E9GvbWjy2TTCw9ikYwSvDsqy8IBwnkh0DFDmfFwEzXKwlnOLhte1KvL7mEFuveyw9l87o5zQ+xevrHkuveSw9l87E688oxP5r2x297j0IRvDKBCosy1oahGNFqox2/x/pRhEuiIfXtisZ7f4fL685xNbrHkvPpTMZ7f6vn9P4kNHu//HyusfSax5Lz6UzGe3+Hy8/oxD7r213Mtr9X697O4k9HcXpNsxsd3NRoBvRGLPMuZ/fsqzZXXS/RTs/Hb8HMS9OXtuu+InD1xyi63VX+6Sf0xh+bbvjJw5f92h6zdU2xefPKMTFa9sdP3rdu33dewxe3f3AOBPn1hljvP2uNDqspSXF5wGru75rbIij17YrcfeaQ/S97mqf9HPqdi0uibvXPdpec7VN8fczCnHz2nZHr3sPBjrU8BnnhIHJdDGZ6p3nledMEswIXLGKcXHx2gY4r+2sds835l/z9s+b2HrdY+m5dCpefk47EfOvbWvx2D6pbYpu8fAz2oWYf21bi8e2CQbWPlnGmNBXKCIiIiIiEse0gbKIiIiIiEiIKXiJiIiIiIiEmIKXiIiIiIhIiCl4iYiIiIiIhJiCl4iIiIiISIgpeImIiIiIiISYgpeIiIiIiEiI/X8iTJLVU8IlOAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "mini, maxi = -5, 15\n", "X = np.linspace(mini, maxi, 100) \n", "fig, axes = plt.subplots(1, 3, sharey = True, figsize = (15, 5))\n", "epsilon = 2.0\n", "for y, ax in zip([-5, 1.5, 5], axes):\n", " # plotting histogram of actual samples\n", " ax.hist(\n", " data_xy[\n", " np.logical_and(\n", " data_xy[:, 1] > y - epsilon / 2, \n", " data_xy[:, 1] < y + epsilon / 2,\n", " ), \n", " 0\n", " ], \n", " # bins = 100, \n", " density = True, \n", " label = \"samples\", \n", " histtype = \"step\",\n", " color = \"red\",\n", " )\n", " # pulling conditional samples\n", " data_x_cond_y = kde.sample(\n", " conditionals = {\"y\": y},\n", " n_samples = 10000,\n", " keep_dims = False,\n", " )\n", " # plotting histogram\n", " ax.hist(\n", " data_x_cond_y, \n", " bins = 100, \n", " density = True, \n", " label = \"KDE samples\", \n", " histtype = \"step\",\n", " color = \"blue\",\n", " )\n", " \n", " # calculating analytical probabilities along the axis\n", " x = np.stack([X, np.ones(len(X)) * y], axis = -1)\n", " probs = np.exp(kde.score_samples(x, conditional_features = [\"y\"]))\n", " # plotting the function\n", " ax.plot(X, probs, label = \"KDE analytic\", color = \"black\")\n", "\n", " ax.set_title(f\"$P(x | y = {y})$\", fontsize = 20)\n", " ax.set_xlim((mini, maxi))\n", " ax.set_yticks([])\n", " ax.set_xticks([-5, 0, 5, 10, 15])\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ZCA whitening\n", "Instead of a simple rescale of dimensions, potentially more powerful whitening can be obtained with ZCA algorithm. It scales the data along the principal axis. The method comes with a bit bigger computational overhead. Here we show results for the same example as before." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA14AAAFCCAYAAADhUww7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABcTklEQVR4nO3deXxU5fn//9dJQkgIgSwsAQITAuJSV0CxWleiVsH60YLW2tpFBe1i/X2t4K61i4KttrXVD2ity6cqS10gbiWguNSiIXUDFSRkQpBAtgkhKyHn98eZmez7zJyZOe/n45EHyWznmsxwZ65z3fd1G6ZpIiIiIiIiIsETY3cAIiIiIiIi0U6Jl4iIiIiISJAp8RIREREREQkyJV4iIiIiIiJBpsRLREREREQkyJR4iYiIiIiIBJkSrzBnGMaCQdx3iWEYKQEMJ+IM5vcnIt3T2DRwGpdE7KOxS2OQneLsDiDaed/c84EcoBDIa3N1GpACLDZNs6CL+y4B7h3E4bO9x/AM4jEi3UrDMJaYprnY7kBEwonGJltpXBIZII1dAaExyCZKvILMNM3lhmHkA5uxBoLVba83DGMesNkwjBltBwnDMHKAD0zT9IQ04AhiGEY2sApYBuRjDYjnAMt8v0vTND2GYewwDCPHNM287h9NxFk0NvWdYRirgHu7+iDXxW01LokEkcaunmkMCm+aahgaOd5/O7252wwYCztctbDjYCJdSgGWYA3AC2kzsPiYprmczr9fEdHY1C3DMLINw1jmPUOeg3WWu69S0LgkEkwau3qWgsagsKSKV2icCBT29SyLYRjTgQ+CGlH0mN+Xs9DAB4ZhTO/jbUWcQmNTN0zTLMT7ocR7Br0/NC6JBJfGrp5pDApTqniFRg5dnJWBdgscV7W5+DLAKWdlQiUPndkR6Uhjk700LokMjMauwNAYFGKqeAWZd65tCrCum5ssAZZ3mGOb092CR+/j5XgfE9M0l7YZZGaYptnjf6A2959Bh9KzYRjLgM3e8nPE8HYYygY83rPUnZimWWAYxsyQBiYSxjQ2BZfGJZHg0NjVNxqDwpMqXsHXaR6yYRgphmHkGIaxDmvBdp/ONnj/Ey00TXO5aZpLgcu8/6lXYi2gXOAdAHoyzzsAeOh8lmMBVoegSHIZMBMr7hTDMFb10Oq1u8tFnEhjU/BoXBIJHo1dvdMYFKZU8Qq+c7D+My4wDMN3mQdrbvI5HW/s/Y9R2c1jLaB9G9RKIM3bncaD1d2n2//g3o4+vlJ7DlbHG9910wF6627j7e6V0tNtOljYU0yDYZpmoWEY97aZ413gHXQfxWo121GhYRjZwYpHJMJobArCWKBxSSToNHb1EJPGoPCmxCv4coA875mUvuhpf4jVHRaS5uD9T+T9D9PjMXz/+b1nb6bTfn50DtDr4krTNLv6T2ubLhbW5gHLDMNIifaWsSKDpLEpSDQuiQSVxq7eH9PT4SKNQWFCUw2DqA/zkLuS0t0Vbc9GeM+ytG2b2h8LgYIOZzfOoZuFqsHkbdm8uR9fS9rct6ud131ntXqbGiDiWBqbgkfjkkjwaOzqncag8KaKV3B1u89EDzz0reTc6T90P0rFXXUDalci7463JN4fPZbpvdfN6Odj+gbfZYZh5HXz+N0ds7vpBiJOorGpl7FpIDQuiQSdxq4exi6NQeFPiVdwnUMP3WS64p2b2+VGnYZhLGpTWp9Hm9ao3rnEafRtEWc2bfaz8M1Dpg8DWbhM5/H+nrqa53wp1lknTxd3y1aJXQTQ2BQUGpdEgk5jV8+PpTEozGmqYXB1u89EfxnWBp63+Dr30Hne8GW9LeBso+N/yFvox0aEYaSybbchX3ci4Jpubu8JQUwikUBjUwB4n/Nmo/0GyxqXRIJHY1fvNAaFMVW8gsC7Dsk3D3m69+d7+/EfML+L8nYesBzrrEWhaZrzDcNY5pvL293+FN24BmuwSfPGmE0EraHwMU1ztWEY89p86JmCtVt7p7NT3rNPK0IaoEiY0djUd94PK7e0iWOJYRh5wLoOH8aysc6KAxqXRIJBY1ffaQwKb4ZpmnbHIB14z7xM70fHnu4eZxV9WMdgGIaJtUlgr913IpVhGIuwuiBF7XMUCTaNTYGlcUkkNDR2dU1jUOhpqmEY8p5NPTEYj20YxjrvfzTfz4uw5v1G+3+6Ex3wHEWCSmNTwGlcEgkBjV3d0hgUYkq8wteyDusGAiUN7+JRb4n5MmB2EI4TNry/x147C4lIn2hsCgCNSyIhp7GrDY1B9lDiFaZ8Z2e86wwC6V6s+dGLsBaOzgjXheuB4P39ndiPBbIi0gONTYOncUkk9DR2tdIYZB811whjpmku9i7yXD7Ahyikw94MA9wYMJJd2s8FsiLSC41Ng6ZxScQGGrv8NAbZZMDNNUaNGmVmZWUFNhoRsdXmzZvLTdMcbXccg6GxSST6aGwSkXDVn/FpwBWvrKws8vPzB3p3EQlDhmG47Y5hsDQ2iUQfjU0iEq76Mz5pjZeIiIiIiEiQKfESEREREREJMiVeIiIiIiIiQaauhiIiIiIiEeDgwYOUlJTQ0NBgdyiOk5CQQGZmJkOGDBnwYyjxEhERERGJACUlJSQnJ5OVlYVhGHaH4ximaVJRUUFJSQmTJ08e8ONoqqGIiIiISARoaGggPT1dSVeIGYZBenr6oCuNSrxERERERCKEki57BOL3rsRLRERERCQSZWWBYQTuS5t8B5USLxERERGRSOR2g2kG7stt317leXl5nHPOObYdPxSUeImIiIiIiK1ycnJISUmxO4ygUldDERERERHpVUFBAZWVlXg8HlJSUsjJySEvL4/CwkKys7P9Py9ZsoTFixdTUFBAdnY2KSkprFq1iiVLlpCfn8/ixYtZsmRJu/u15fF4WL58OdOnT6ewsJCZM2d2Om4kUsWrn7KyNP1VRERERJxnxYoVAMybN4/s7GwKCwspLCxkwYIFLFmyBLAqV4WFheTk5DBv3jyWLVtGTk4OM2bMID8/n5ycHNLS0sjJyWHBggUsXLiw03HuvfdecnJyyMnJYfPmzZ2OG6mUePWT223r9FcREREREVvccsstLFu2jClTpuDxeMjOzmbBggV4PJ52t5s+fXqn79PS0vyXtZ1S6Evg2vJV1goKCli4cGGn40YqJV4iIiIiItKrvLw8Vq1axebNm8nLy6OgoIClS5f2+3HaJk++6YZt+ZpsTJ8+nezs7E7HjVRa4yUiIiIiEolcLqsNfCAfrwcffPABYFWp5s2bR0FBASkpKf7kafXq1WRnZ1NQUEBhYaE/OSssLGTdunX+9Vm+alZ+fj7Lli0DrCqX77aLFi1ql9B1PG6kMkzTHNAdZ86caebn5wc4nDCUldVubqGB9fsyafMmd7mgqCi0cYkEgWEYm03TnGl3HIPhmLFJxEE0NolYPvvsM4488ki7wxi0+fPns2rVKrvD6Leufv/9GZ801bA3HfdH8AmTPQ9ERERERCJF2yqY02iqoYiIiIiIhEROTg47duywOwxbKPHqI7WQFxERERGRgVLi1UeaTSgiIiIiIgOlNV4iIiIiIiJBpsRLRERERCQCZWVZ3eQD9aWlNcGlxEtEREREJAJ1bL492K/eltYUFBQwY8YMFi9e7N8Eef78+f49t3zXL126lNWrV7N8+XKWL1/e6f5Lly4lLy+P1atXM2PGjGD9ejrJy8vzb85sB63xEhERERGRXk2fPp3s7Gwuu+wy/8bJt9xyC9OnT293fU5Ojv8yX/K1YMGCLq/Pzs7G4/GQkpIS9PhzcnL8GzbbQYnXAPlKsdo3WUREREScpqCgAMCfQHXn0ksvZcaMGSxYsKDTdatXr2bevHmdHreystKfjOXk5JCXl0dhYaE/acvLy2PJkiUsXryYgoICsrOzSUlJYdWqVSxZsoT8/HwWL17MkiVL2t2vLY/Hw/Lly5k+fTqFhYXMnDmz03EDTYlXF7pOqpqB+4Dzcbmskqg6HYqIiIiI0+Tn57N582amTJnSa+Llq4y1lZeXR35+fpe3X7FiBeeccw7z5s2jsLDQ/7VgwQLOOecccnJyyMnJYeHCheTk5JCdnc3ChQtZt24dhYWF5Ofnk5OTQ1pamj95mjJlSqe9w+69914uu+wypk+fzsKFC9mxY0e74waD1nh1we3uKql6G9gEPMnOnaYqXSIiIiLiSGlpaSxbtowVK1b0mqR4PB6ys7PbXZaTk8OCBQv8iZGvegZwyy23sGzZMqZMmeK/74IFC/xrynzaJny+79PS0vyXtZ26mJ2d3SlOX2WtoKCAhQsXdjpuMCjx6gPTNIEXsAqEO/nkk09sjkhERERExB6+ROrRRx9l/vz5Pd525cqVLF68uNvH8Xg8VFZW+i/Ly8tj1apVbN68mby8PAoKCvzNO/qjbfLkm27Ylq/Jhm/dWcfjBoOmGvbBRx99BOwErgOe4cUXX+TYY4+1OSoRERERcTKXy2oDH8jH60lBQQGFhYWsWLGC7Oxsf+K0cOFC/3qqwsJC8vLyqKys9FeZfOu7CgoKKCgo8F/v8Xi49957WbVqlf8YH3zwAWAlZfPmzaOgoMA/XTE7O5vVq1eTnZ3tj8WXnBUWFrJu3Tr/+ixfNSs/P9/fUMN3/MLCQhYtWtQuoet43GAwrGpO/82cOdPsbm5mpPO9gU3T+uGuO+/knnsKgceB1cyd+wyPPPIIEydm+m/DAH+PIuHEMIzNpmnOtDuOwYjmsUnEqTQ2iVg+++wzjjzySLvDiAjz589vl9AFQle///6MT5pq2As3vnmnFwJDgAsYMmQIL730kr2BiYiIiIhIJ22rYOFEiVcvXgHi4+OB872XjOT000/nzTffBFTlEhEREREJJzk5OezYsaPTui67KfHqwaFDh3gXOOmkk4Bk/+VTpkyhoaEB2G9XaCIiIiIiEkGUePXg448/pho444wz2l0+evRo73dlIY9JREREREQijxKvHrz11lsMo/OO3KNGjfJ+Vx7ymEREREREJPIo8erWQd577z2+DkybFt+uvaYqXiIiIiIi0h9KvLpVQG1tLacBbjcUFbVeM2LECIYMGYIqXiIiIiLiFAUFBcyYMYPFixf7NyieP3++fz8s3/VLly5l9erVLF++nOXLl3e6/9KlS8nLy2P16tXMmDEjoDHm5eX5N0fuzurVq/3fFxYW9roJdKAo8erWW4wYMYLjurjGMAzvdENVvERERETEGaZPn052djaXXXaZf1PjW265hUWLFrW7Picnh3nz5vk3TvYlX22v993m0Ucf9SdxgZCTk0NKSkq313s8HtatW+f/OTs7O+D7fXUnLiRHiTjNwCZOOeUs4p55pstbWNMNVfESERERkdB79NFHA75PVXZ2Ntdcc02fbmvtc9u5F0JHl156KTNmzPAnYW2tXr2aefPmdbo8Ly+PwsJCf5KWl5fHkiVLWLx4MQUFBeTk5PiP2/G2HR//3nvvZf369eTl5bFu3ToWLlxIfn4+eXl55OTkUFBQwOLFi/3J2NKlS8nJyaGysrLT4w2WKl5d2gk0ctxxx5HFznbru3xU8RIRERERJ8rPz2fZsmXk5eX1eltfZaytvLw8li9fTmVlZafbFxYWUlhYyIIFC1iyZAlgVbEKCwv9VbIVK1Z0e9u25s2b59/LKyUlhWXLlrWruoGVOPoqZKtXryY7O5vp06e3q4oFiipeXdoGwOGHH46b0ZhFnW9hJV6VtLS0KHsVERERkZDqa2UqGNLS0li2bBkzZsxol9x0xePxdLreV7HyJWQFBQX+ClZ2djYLFizoNP2wq8pad7dta+HChSxfvrzXyhzgr4gBXSZyg6WcoY2sLOsLvgBS27SN78yaatjSZaYuIiIiIhKtfInUo48+2mtjipUrV7J48eJuH8fj8bT7PF1QUOBv1tGbvtw2JyeHFStWMHPmzC7v39aUKVP8sQRy3ZmPEq823G7ry0q8pmEYRre39SVl5eVa5yUiIiIi0a+goIDCwkJWrFjhr2R5PB4WLlyIx+PxX5+Xl+efTgj413cVFBRQUFDgv3716tXMnj27XUWssLDQPz0xOzub1atX++/ne+yCggI8Hk+vt/VZsmRJu4Yb2dnZ5OXlkZ2d3e72ixYtYt26dRQUFJCfnx/w359hmuaA7jhz5kwzGAHZycqzaoDvAldimvMxDPD9inx5mGlCUVERkyf/nLfeWsRpp5/eeiORCGYYxmbTNDufEoog0Tg2iTidxiYRy2effcaRRx5pdxgRY+nSpSxatKjdVMbB6Or335/xSRWvTrZ5/z28x1v5NlFWxUtEREREJPxMnz6dvLy8gCRdgaDmGh2MGrWN2lqD+vrDerxdUlISkEhZmTobioiIiIiEm0C3gx8sVbw6+OlPv+DHP54EJPbh1qNU8RIRERGRkBnoMiEZnED83pV4tWOybds2pk2bhstlrelyUdTD7UdTXl5OFju93RBFRERERIIjISGBiooKJV8hZpomFRUVJCQkDOpxNNWwnT3U1NRwxBFHUFTkvciYDHT35h7Fvn07cJMF7pAEKCIiIiIOlZmZSUlJiZa62CAhIYHMzMxBPYYSr3a+AGDatGl9vP1oqqurgYPAkGAFJSIiIiLCkCFDmDx5st1hyABpqmE7xcTFxTFp0qQ+3XrMmNGsXw+gdV4iIiIiItI9JV7tlJGenk5MTN9+Lf/61yjq60GJl4iIiIiI9ESJVzvl/v25+mLUqFH++4mIiIiIiHTH2YlXVpbVutD3RRmj77uv/WUuV7d3V+IlIiIiIiJ94ezEy+0G0wTTpOXQIaCc0U895b8M06S1vaGVg7XNw4YOHcqwYYmk8WXIQxcRERERkcjh7MSrjaqqKqClx6mGRUXt8jAAvve9VBZzazBDExERERGRCKfEy8u3H0Lr9MG+SU1NpSoYAYmIiIiISNRQ4uV15pllJCbSr+YaACkpKXiCE5KIiIiIiEQJJV5ee/eWMXt2/xMvVbxERERERKQ3Srz8yklKSmLYsGH9uldqaiq1ADQFIygREREREYkCSrz89vV7fRdYUw0t1QGNRkREREREoocSL7+yfk8zBKviZdGEQxERERER6ZoSL7/yASVerRUvJV4iIiIiItI1JV5AY2MjsH9AUw1bK16eQIYkIiIiIiJRRIkXUF5eDvS/oyHAyJEjvd+p4iUiIiIiIl1T4kXr5skDSbyGDBlCMqDES0REREREuhNndwB2y8qC+vqBJ14AKYCmGoqIiIiISHccn3i53QBlgEF6evqAHsNa5aWKl4iIiIiIdE1TDQEr8UolLm5geWgKoIqXiIiIiIh0R4kXAOXAwKYZgipeIiIiIiLSMyVegFXxGnjilQJAAw0NDYEJR0REREREoooSL0ysxKv/e3j5+HbyqqpS1UtERERERDpT4kUd0ASkDfgRUrz/ejyewYcjIiIiIiJRR4mXvylGak836pEqXiIiIiIi0hMlXv6mGINPvFTxEhERERGRrjh+Hy/wkJgIaWkDT7xGAGCo4iUiIiIiIl1SxYsqZs+Gzz4beOIVC8AIVbxERERERKRLSryoIi4ujuHDhw/ycVJV8RIRERERkS4p8aKKlJQUDMMY5OOkqOIlIiIiIiJdUuJFFampA59m6JOamsr//q8qXiIiIiIi0pkSL2/Fa7D+8IdUamurME1z8CGJiIiIiEhUUeKFJyCJV1paGnCQ2traQT+WiIiIiIhEF0cnXi0AeAI21RC0ibKIiIiIiHTm6MSrBoCWgCZelZWVg34sERERERGJLo5OvHy1qUAkXtZUQyVeIiIiIiLSmRIvApN4aaqhiIiIiIh0R4kXgUm8EhMTgaGqeImIiIiISCdKvAhM4mVtwJymipeIhJWsLDAM618RERGxT5zdAdjJSpESSEhICNAjpirxEpGw4naDaVrJl4iIiNhHFS8GX+1qlaaphiIiIiIi0omjEy8PACkBfERVvEREREREpDNHJ17BqHjV1dXR2NgYwMcUEREREZFI5+jEywMEJPFyubwLKFIhN5fKhATr57ZfWtkuIiIiIuJYjk28Dh48SA0QkMSrqMhavU4azJ1L1ZYt1s9tv9zuwR9HREREREQikmMTr+rqau93gZxqmMr69XD22WqwISIiIiIirRzbTr61CUZgE6/6eqivV4MNERERERFp5diKV2vilRLARx0BxAKqeImIiIiISCslXgGteHkbbKCKl4iIiIiItHJs4uXxeLzfpQT4kVNRxUtERERERNpy7Bqvn/ykmgbmEvhfQRqwN8CPKSLSP8XFxbz99ttAHY880gxchjU+iYiIiB0cm3hVVnqACUF45FTgsyA8rohI3z399NP85z//AYbxyit1QAZwsc1RiYiIOJdjpxpCNVYzjEBLA/bT3NwchMcWEekbt9vNN77xDWAF48aNA7baHZKIiIijOTzxSgnC41rNOlrXkImIhFZjYyOlpaVMmjQJgCOPPBL4DNM07Q1MRETEwRyceHmAkUF4XGsNRWvXRBGR0CopKSEvz+S7352Ey+VLvKopLS21OzQRERHHcmTi1dLSAuwn0BUvlwvGj7cqXpWV6mwoIvZwu93U18OuXS6KiuCoo44CYOtWTTcUERGxiyMTr5qaGsAk0Gu8iorgo4+sxEsVLxGxS3FxMRDnXdsFEydOBJL47DM1/hEREbGLIxOv/fv3e79LCfhjp6amAgYVFRUBf2wRkb6wEq9MYmNjATAMAzhCFS8REREbOTLxam18Efg1XtYHnVQlXiJiG7fbDUzqcOlR7Nq1iwMHDtgRkoiIiOM5MvGqrq72fheM5hoAaUq8RMQWDQ0N7Nu3D3B1uOZIAD7//POQxyQiIiIOTbxaK14pQTpCuhIvEbGFNc0QOle8DiMmJkbrvERERGziyMTLqngZQHKQjqDES0TsYU0zhM4VrwSmTJmidV4iIiI2cXDiNYLgPf10Dhw4QFNTU5AeX0Ska8XFxcTHxwNjO103depUioqKtJGyiIiIDeLsDsAOrYlXcIwZk8769ZCdXUFc3DiKioJ2KBGRdoqLi8nMzKTjiSWXC37604nAAd59t5ri4hQ7whMREXEsR1a8rDVeKUF7/NdfT6e+HnbvrsA/60dEJASKi4txuTpOM7T2Gfzvfycydy7s2rUr9IGJiIg4nCMTL6viFayOhpCenu79Tuu8RCR0amtrKS8vZ9Kkjo01LFYlDECJl4iISKg5OPFKwUURXZwYHrS0tDTvd5WBf3ARkW7s3bsXgPHjx3d5fXp6OgkJCUBJCKMSERERcOAar+bmZu8GoiMpYjIUBX6R+bBhw4AEVPESkVDydVMdNWpUl9cbhsHEiRNRxUtERCT0HFfx2r9/P+vXw5gxwZtqaBgGw4enk5KixEtEQqe8vBzoPvEC33RDJV4iIiKh5rjE65hjPAC88ELwEi+A669P55prlHiJSOiUl5cTExNDSkpKt7exKl4V1NXVhSwuERERcWDi9dVX+5k9mx4/mARCero2URaR0KqoqCAtLY2YmO6H9okTJ5KYCElJJWRlhS42ERERp3Nc4gUeAEaODG7FKz09ncrKSkAblYpIaJSXl/c4zRCsxGv2bMjL26XtLkRERELIgYlXNRD8xCstLY3m5mZgf1CPIyLiU15ezmOPjcIw6LZj69ixY4mLi6OkRJ0NRUREQslxXQ3BQ2xsLElJSUE9Svu9vIKb5ImImKZJRUUFHs9MzB4K7XFxcYwbN06bKIuIiISYIyteI0eOxDCMoB5FmyiLSCjV1dXR0NAA9DzVEKzphqp4iYiIhJZjE69gU+IlIqHkayXfl8QrMzOTPXv2AM1BjUlERERaKfEKktTUVG9VTYmXiARfa+KV3uPtwKp4tbS0AF8FNSYRERFp5cDEyxP0VvIAsbGxpKamosRLREKhdfuKvlW8LLuDFo+IiIi058DEKzQVL7A6GyrxEpFQKC8v91bZU3u9besYWBPUmERERKSVoxIva+F5Q0gqXuBb56XES0SCr7y83Ftl771ZbWtX19qgxiQiIiKtHJV4VVeHZg8vHyvxKu/1diIig1VeXt6mqU/PEhMTiYmJAQ4ENygRERHxc1Ti5fF4AEJW8crIyABqqanRdB4RCa7y8nJGjep9fReAYRjeqpcSLxERkVBR4hVE48aNA6C0tDQkxxMR56qoqOhz4gUwfPhwlHiJiIiEjqMSL99Uw9BWvJR4iUhw1dXVUVdX1+ephuBLvLTGS0REJFQclXj5Kl6hWuPlS7ysjUpFRILD10q+PxUvTTUUEREJLQcmXsOIj48PyfESEhKAFFW8RCSofIlX/yteSrxERERCxYGJV0qIj5qhxEtEgqq83OqeqjVeIiIi4ctRiZe1xis00wxbjVPiJSJB5Uu8rE3b+8a3xss0zSBFJSIiIm05KvGqqqoCUkN81LGUl5dzMMRHFRHnqKioYMOGEQwdGo/L1bf7WGu8mmlqagpqbCIiImJxVOJlTTUMbcVrzJhx5OWZ7AvpUUXESaqqqqirS8U0oaiob/exKl5w4ICmG4qIiISCYxKvQ4cOeTcyTgnpcd94I4P6etBkQxEJFmsadUq/7qPES0REJLQck3j59vAKdeLl38srpEcVEScZyPrVpKQkEhMhK6uWrKyghCUiIiJtxNkdQKjYlXilpqYC8WgnLxEJiKwscLvbXeQBYDkYRp8fZjgwG7iDXGa51WBDREQk2BxT8fJtnhzqxMswDCBDFS8RCQy3G0zT/9XU2Ej93LnAyHaX9/aVtHs3zJ1L7YYNdj8jERERR3Bg4hXqdvIA45R4iUhQDPSkktZ4iYiIhJYDE68UG45uVby0X46IBFrrNOr+r/ECJV4iIiKh4pjEq7q6mri4OGCYDUfPoJG2yZ+ISGAM9KRSbGwsCQkJ1NbWBjokERER6YJjEi+Px0NKSgrQ98XngePtbFiqCYciEliDaRw0fPhwVbxERKRPsrKsHk7qhDtwDky87DAOgD171NtQRAJrMOtXlXiJiEhf+Xo7dWisK/2gxCskRgOwb98+m44vItGqurqaoUOHAgn9vu/w4cM11VBERHq1detWoMjuMCKeYxKv6upqRo60o6MhQDwpKPESkcAbzNiWlJSkipeIiPRo4sRSvva1O0hKupO6ujq7w4lojki8TNPE4/F4NzO2xxigrKzMtuOLSHSqqqoa8NimqYYiItIT0zQpKVnGJZfA2Wd7eO655+wOKaLF2R1AKNTV1fH66828+65dFS8r8SpUxUtEAqy6uprRo0cP6L5KvEREpCfvvfcekM+VV15NcXExa9asAc4BJpKV1brey+WCoiLbwowYjqh4eTwe6uuhrCzFluO7XHAzf+Ppp8twuUx1gxGRgBnMVMPhw4fT0NAANAc2KBERiXgNDQ0sX74cmMzcuXO58sorSUhIAJZjmqa/2YYabvSdYxIvS4otxy8qgjVcRX39QYqLq/XmFJGAME2T6urqATcO8m2iDGqwISIi7blcb/PkkxWMG7eQ2NhYRo4cyRVXXAF8yMcff2x3eBFJiVeIjPF/p+mGIhIYBw4c4NChQ4OqeFmUeImISHv79m3mBz9IZ/fuo/yXnXfeecAIcnNzu7zPunXr+Oijj0IUYeRxWOJl7xovixIvEQkM3+bJg2muYdE6LxERaXXo0CHgQ6ZPn45hGP7L4+PjgfPYtGkTHT/Tvvfee/z5z3/mnnvuYcKE7dpsuQuOSLwqKyuxnmo4JF5luFx6I4rI4PlOKg2+4qXES0REWm3btg2oZcaMGV1ce77331fbXFbGn/70J6ZMmUJKSgpfffUbKioqtbymA0ckXlVVVUAqdj5dayXFMGAfRUVahCgig5eT4yE3F44/fiQuV//vrzVeIiLSlfz8fCCG448/votrR/P1r38deJ2mpiaampqA+2lpaWH16sU8/vgdDBtWx29+8xsmTWpW5asNRyReVsUrze4w8O7mZXcQIhIl9u6tZu5cqKpKGVAbX1W8RESkKwUFBcARbU7QtXK54NZb55KYWMOtt97KFVdcQWLiZ6xa9VPi4sZhmlmsW/f/sX37dv7whxfV9bANR+zjVVFRAYy1OwxgNFrjJSKB48EwDJKTkwd0byVeIiIO1XYTrg6qgS8BWAFt1nf5FAEmcFM9lP4hlzOAO4BjXZ/4N/M65ZRTmDVrFitWrODMM88ERgX6GUQkVbxCSomXiARSNcnJycTGxg7o3kOGDCEuLg5NNRQRcZi2m3B1+CrYsAHmzgWmd3sbwzS5v6WFp1ta+JlpcmwXZa1rrrmGlpYWHn/8ccDaAmXixAoM47+MGbOW3bt32/DE7RX1iVdzczP79+8nPBKvMUAtdXV1dgciIlFh4Ht4ARiG4a16qeIlIiKWzZs3e5s2TenxdoZhtOt42NHYsWOZN28eb7/9NnAvV111FSUlP2Tu3DspK1vOn/70p8AGHgGiPvGyGmtA+CResG+fql4iEgieQSVe4GuwocRLRESgtraWTZs2cdJJJwHdJ1V99e1vf5vMzExgC9OmTQOu4be//S2jRl3J/fd/xvjxWwd9jEgSnYlXVpY1J9UwqBwzBnJzaZd4ea8bUBuwQbESr7IyNdgQkUDwDLiVvI9V8dJUQxERgTfeeIOGhgYuuOCCgDxefHw8Dz/8MPA0N998M/Atjj32WEpKvsV3vzuCPXtWB+Q4kSI6E68281Yr3n3XO0+1TeLlm6M6kDZggzIaUOIlIoEyuKmGgKYaiog4nK9e4XKZvPLKK0ybNo2pU6f2/4FcrtbiRpsvIyYGMFobdRgGQxMSuPCZZ4APKOriPu2+oqgXfXQmXm1YjTUA0m2Nw5IKxGmqoYgM2sGDB4HaQVe8rKmGqniJiDiVr15x8OCnPPLILv7v/wZY7Soq6rYZh8sFBta/vsvm7N8PJPD8H/7Q7f2irRe9IxIvq+PXCFwuG2YXtmMAo5V4icigeTwegACt8VLiJSLidNdf/zKXX57Mvn2nBfyxfTlZ28lm1lYo57Fx40b/37Ro54jEKzU1FTAoKrJhdmEnSrxEZPB848jo0aMH9Ti+xMs0zQBEJSIikaiiooL33nuPnJwcID6ER55NS0sLmzZtCuEx7eOIxCstLRw6GvqM0RovERk0X+I1ZsyYQT2OlXgdoqmpKQBRiYhIJHriiSeIiYlhzpw5IT5yFmPHjuW9994L8XHtEWd3AMH2xz9WEhs7zu4wAGua44ED6d4W94eAgW16KiLiO4ETmIqX1UJ46NChg45LREQizce8+eabfOc732Hs2LH+HhmhWZ5jcMopp7B27Vpqa2v9f5OiVdRXvGpqKvF4wqPiVVQETz+d7p3S47E5GhGJZFbFa+Sgk6W2iZeIiDiL1ajpYTIyMpg/fz7Q9XqsYPr6179Oc3Mz+fn5oTmgjaK64mW9mWoIj82TLenpvu6KFYRHp0URiTRZWeB2l5GcPLhphqDES0TEyV588UVgN9deezfx8aFc22VxueCoo44AUnnnnfeoqjoj5DGEUlRXvKwpfRC+iZeISP+53XDttfu45ZbBTTMEJV4iIk7V3NzsTbxOZMaMGbbEYFXXDP7615PxePJpamry7ysWRdt3+UV14lVR4Utuwi/xGj26IirfUCISCib79u0b9PouUOIlIuJEWexkyJB8nnlmPxkZA9y3K4C+/vWvA43897//9e8rFkXbd/lFdeIVXpsnW0aOHElcXBz3318RlW8oEQmF/TQ1NQ26oyEo8RIRcSI3Wfz613lceWUaJSUn2B0OxxxzDDCcjRs32h1KUDkk8QqfipdhGKSmprapxomI9FdgWsmDEi8REWeqJj8/n7POOovYWPu7bMfFxQFnedvKV9sdTtA4IPGKA5LtDqWd9PR0JV4i0med57tbreQDkXhZi6ljlXiJiDjKmxw6dIjZs2fbHUgb59Pc3AystzuQoHFA4pUKGHaH0o4SLxHpj87z3QOzhxdYVXhIUuIlIuIQ1rZGeRx++OFMnDjR7nD8XK6J5OYexfDhr3ljjD4OSLzCZ5qhjxIvERmcfSQkJDB8+PAAPZ4SLxERJ8jKgpiYnSSyNcyqXVaHww0bvsmZZ+7h448/tjucoIjqxMtKbsKnsYZPWloaDQ0NQJ3doYhIRLI6GlrVqkBQ4iUi4gRuNzz++EbOZQPf+MY37A6nk1NPPZXhw4fz2muv2R1KUER14uWbauhyWRu0hQvt5SUig7MvIOu7WinxEhFxBpO3336b6UBycnj1QABr3fHs2bO9TTbK7Q4n4KI28WpsbPR+kEinqMgqX4YLJV4iMjhlSrxERGQAPqesrIzT7A6jBxdeeKF3jdcLdocScFGbeIVjK3mf1sSrssfbiYh01ABAjRIvEREZgLcZMmQIs+wOowdjx47lzDPPBF6jujq6WstHbeLV2rwi/NZ4qeIlIgNV5v030IlXXZ3WnIqIRLOWlhbgHWbOnMkwu4Ppxfz584GDrFmzxu5QAipqE6/Wilf4JV5Dhw71blqqxEtE+mef999AtJJvlURDQ4N3/xQREYlGW7ZsAao4/fTT7Q6lV5mZmcAp5ObmEk3zMaI28WqteIXfVEPwVb2UeIlI/wQr8QKor68P4GOKiEg4eeutt4AEZs6caXcofTJ+/KWsXFlHFg/ZHUrARHXiNXToUAjTYqoSLxHpO5O8vDxgOa8AEEtaWiBPKlnj5IEDBwL4mCIiEi4aGhq8idfXSUhIsDucPtm9O5s77jiBcj6ImhkZUZt43XRTJa+/ng4Eap+bwFLiJSJ9t5Y//elPQB7Wn545xMQEcvi2NmJWgw0Rkej09ttve9fynm93KP3yrW99C6jk3XfftTuUgIjaxMvjqeDAgfCcZgi+xKuKQ4cO2R2KiISxbdu2AX9n1qxZwAoeAeCaAB/FmmqoBhsiItHp1VdfxeVyAUfYHUq/zJgxA5jASy+95G0xH9miNvGyWrWHX2MNn/T0dBITTbKyPHaHIiJhqqamhiVLlgBp3HDDDQSrgj+eSnJz4dxzNdVQRCTabN++ne3bt3P++ecTrjPBumMYBvAttm/fzhdffGF3OIMWlYmXlQ9XEK6NNcBKvGbPhpISTTcUka49//zz3kZBixk+fHjQjvNfzmDuXNi7VxUvEZFo89prrzF06FDv3liR6GySkpJ46aWX7A5k0KIy8bLO2R4k3CteFiVeItK1L774gilTpgDTgnqc1hZEWuMlIhJNamtr2bhxI2eccYZ3K6NIlMB5553Hu+++y65du+wOZlCiMvFqTWXCKPFyucAw/F+jp06F3FygrPXyrCy7oxSRMGGaJjt37iQ7Ozvox2pNvDTVUEQkmrz//vs0NjZyzjnn2B3KgLko4sc//jbr1w/jb3/7G1lZkfuxOSoTr0r/d2GUeBUVgWn6v5JbWkiYN4+0tFKyXN7L3W67oxSRMFFeXs6BAweYPHly0I8VAwwbNgzQVEMRkWjy/vvvk5qayuGHH253KANWxGRMcwR1dZezefNm3O7NEfuxOSoTr9aKV/iu8TIMg4yMDH7xi9KIfOOISHAVFhayfj3MmTMZlyv4x7MSL001FBGJFs3Nzdxzz2aefvpEJk+OrKYaXZvD+PHjgccidl+vqEy8Wite4Zt4AWRkZLB37167wxCRMLRz507q6w3q6rIoKrIuc7nAwAxKImY171DiJSISLT799FPq6urZtGlWlJzkj+Oqq64CSnj11VftDmZAojLxsipeybhcQ0JypnigMjIyKC0txdeHUUTEZ+fOncA4EhMT/ZcVFYGJ4U/EAsladK3ES0QkWmzatAmI57jjjrM7lIA58cQTgWN59tlnicS/WVGceKVTVERQPqAEytixY2lqagI8dociImGmsLAQCH5jDR9NNRQRiR6mafL+++8DxzN06FB/j7dwLkj0xuWCmBiD8eN/TE1NDfBPu0Pqt6hMvKyphmHUWKMbGRkZ3u9KbY1DRMJLXV2dtxoe/MYaPppqKCISPdxuN/v27QNmAa093sK5INEb33PYvXuKd0+ylygvL7c5qv6JysTLqniF9/ouaJt4aZ2XiLQq8v9lDF3ipYqXiEj0sKYZApxoaxzB8v3vfx9o4f/+7//87eUjocV81CVehw4d8k7cC/+K15gxY7zfqeIlIq2s9V0QyqmGvjVepqk1pyIikW7Tpk1MmzYNSLU7lKCwPkPPYcOGDbjdlf4dm8K9iUjUJV5VVVXeVhXhX/GKj48nPT0dJV4i0lZhYSEjRowglOOYNdXQpKGhIWTHlP555plnePzxx+0OQ0TCXGVlJdu3b2fWrFl2hxJUEyacx9q1JqNHv2N3KH0WZ3cAgVZZ6WsmH/4VL7AabGiqoYi0tXPnTu/GyaHbd8Waagi1tbXtOilK+NiwYQN79+7lsMMO47TTTrM7HBEJUx988AEAJ510ks2RBFdJyUR+8Yts4uI2At+yO5w+ibqKV0WFb/vk8K94gbXOKympNOznpIrIALSdeN7L1xbD4C7DYKRxDQ8+uIPf/jYbF0WdbxukllTWVEMr8ZLw09TUxL59+zAMg7/+9a8Rt6BcREJn06ZNjBkzBlcktzDsozPOOINt27axZ88eu0Ppk6hLvFr/GEVGxSsjI4Ozz67A7T5odygiEmhuN/6J5718nTry19zDfIyRGdx33ym43TkUmVmdbxukllRJSUkkJkJWVq1OBIWhr776CtM0ueKKKzh06BAPPPCA1uOJSCeNjY189NFHvPjiLGJijIhuH98Xp59+OgBvvfWWzZH0TdQlXnv37iUegBR7A+mjjIwM7x/PfXaHIiI2OXjwINXVH/Hww2fh8fyaxYsXM2nSpJDGkJSUxOzZsGnTgbBfnOxEu3fvBmDmzJksWLCATz75hJdfftnmqEQk3Hz44Yc0NTVRWTkr4tvH98WoUaM4+uijefPNNyPiZFTUJV6lpaVYTdpDtzZiMKw1XqAGGyLOtXXrVqCBGTNm2BZDWpo1Pbt1nayEE1/iNWHCBHJycjjhhBN48sknKSsrszkyEQknmzZt8k4d/5rdoYTMGWecQUlJSZuOwOEr6hKvvXv3ktH7zcKG9vIScbasLDj++M0kJsZx7LHH2hZHamoqhmG0WScr4aSkpIRRo0aRkJCAYRj89Kc/xTRNHn744X6d5S0uLub++++nqakpiNGKiB1M0+T999/3nsSLuv553Tr11FOJjY2NiOmGUZV4mabZpuIVGVJTU4mPj0cVLxFncrvhuuvyufHGo0lISLAtjtjYWFJTU9W0IUzt3r2bCRMm+H8eO3Ys3//+98nPz+ftt9/u02OYpskjjzzCW2+9RWFhYbBCFRGbfP7551RXV0d9N8OOkpOTOe6443jnnXeA8J5uGFWJ1/79+2loaGBs7zcNG4ZhqKW8iKOVsWvXLlunGfqkp6er4hWGTNP0J16+RplZWXDhhReSnZ3Nc88916eq1/jx73PffZ+SmwtnnhkZHcBEpO82btxIfHy84xIvgG984xvs3bsX2AHQbqwMJ1GVeJWWWlWjSKp4gW+dl/4IijjTZoCwSLxGjRqlilc/HDhwgKuvvppPPvkkqMeprq5mzZpafvKTTMBqbul2Q0xMDOeffz67du3qcW2D9QGkmQMH/s6CBeO58EKDvXs1y0Ikmhw6dIh33nmHk046Kfr2YnS5et2S5eRzzyUmN5d0XsIwAHcRJobVLKqr+9iUkUVV4mVlupGXeI0fPx7YE+bFUREJjnzGjBlDZmam3YGQnp6u5hr9sGXLFvbu3cunn34a1OOUlJRQXw+bN0/o1KHs1FNPJS4ujjfeeKPb+7vd8PLL/+LMM3dz1VVXkZ6ejk72iUSIPu4H+VFcHNX/+Ad/vHmateVjV/tAhmBPyKAoKup1S5Zk0+S4O+7g4qs/p6XFbN2OBbq+j03te6Mq8fJVvCJpqiHg/cDVgCb4iDhLS0sL8BHTp0/HMOzvxJqenu7dQLnB7lAiwueffw4Q9I07S0pKAMi86KLWD00AhkHyiBHMfPFFNl59NS3dfcAC1syZw1G5uZw4axbjnngCKA3PeTgi0l4f94Pc+MADJF12GaXMtdrId7UPZAj2hLTTaaedRmlpaVivYY26xCs1NZWhdgfST74F0yU2xyEioeXxeIAGJk+ebHcogDXV0KLTQD3xnYD+4Q8/Z/361pN+wWK1ko9ndElJ64cm8H9/1rvvUjV3Lh//979df8Ciia8uvJDjn3kGwzTJ+NOfgD22nvUVkcBpamrivffe45RTTgGG2B2ObU4++WRiYmK8TTbCU9QlXhkZGWSxM6IqqL4pRkq8RJzF18iiNeGxlzUFDUDrvHridkNz8yG+/e3tANx/f2lQC0dW4jW+290pZ86cSVJSUg/TDXdjmiYTJ04EYNy4cYCHhgZVNkWiwQcffEB9fT2nn3663aHYqm13w3DdTDmqEq+9e/eSkZGBm6yIqqCmpqYybNgwfsBv7A5FRELI18iiNeGxV2sCqMSrN263m8bGRn78Yxdz51bhdgcvibGmGna/BjA+Pp5TTz2Vf//7390kU7sAOiRewZ8iKSKh8cYbb5CSkmLrXpDhIicnh9LSUjZt2mR3KF2KmsSrubmZZ58t48Ybx1oLCiOIYRgsXDiBSurtDkVEQsiXeIVfxUtTDXvjW9915plnei8JznTD5uZmb+OoCT3e7tRTT6WhocEfV3u7MAzD28gJMjIySEyE7Ow9ZNF9N0QRCX+VlZV88MEHzJ49m5iYqPlYP2CnnnoqGRkZrFy5MiyrXlHzCpWVlVFfb5KXl0ER4bFeoj+s6Ya77Q5DRELI6iAYx4gRI+wOBbAqJ8nJySjx6t3nn39Oamoqxx13nPeS4OzFuGfPHm8Tlp67Xh522GEAfPnll11cW8z48eMZMsRa+zFu3Dhmz4Z//rMUN1mBDVhEQiovL4+WlhbOPfdcu0MJC7GxsXz7299m+/btfPTRR3aH00nUJF7+joZjI62nocVKvMqpr1fVS8QprIpXelh0NPSxql5KvHrz+eefc/jhh5OR4dvAJDjT9nbt2uX9rueKV3JyMmPGjGHHjh1dPYp/miFAUlISycnJmmooEuFM0+T111/n2GOP9Ve0BWbPnk1aWhqrVq2yO5ROoi7xav0jGFl8DTasRdQi4gS+xCucWImX1nj1rJo9e/ZwxBFHkJycTFJSEsFKvHyt5GFiu8t9+4m2beoxderUTolXc3Mz8FW7xAusqlewuzGKSHB9+OGH7Nu3j29+85t2hxJWhgwZwsUXX8zHH38MfGF3OO1ETeL1/e+XkpgYFzaL1PvL39mwRL0NRZzC6moYHuu7fKz1Zqp49cz6Q37EEUcAvhN+wUlidu3axejRo4GEdpf79hNt2w3+gQem8Oije5g0qdZ/mVXVOtQp8crIyFDFSyTCvfbaa4wYMYKTTz7Zv81FJHX1DqZvfvOb3qnzz9kdSjtRk3hVVOzlyivHhtWUnf6wukwZqniJOIRpmt7EK7xOFlknrzzeSol07QtiY2P966qs8Tt4iVfHpKktX+XLMCAubgpz58KuXTvb3R/osuK1b98+QK+zSCTas2cPmzZtYvbs2QwZMsS/z3IkdfUOpoSEBC655BIgn23bttkdjl/UJF5QGrHTDAHvoucMVbxEHOLAgQM0NTURnhUvX+MP6VoJ48aNIz4+HvBVvPZ5m2AEjmmalJSU9Jh4+Spfpgmffz7Ve2lrgw1f4uWbVeEzbtw4b8evfQGNWUSCr7Gxkd/97nckJibyrW99y+5wwtacOXOAZJ599lm7Q/GLisTL+mO3OwoWFmYq8RJxCF8r+XCreKWlpQGtmztLV8oYM2aM/6dx48aRmNhMbGx5QDdSLisro7GxsVPS1J2RI0d6E+fWdV5W4jWGhIT2UxWD3RRERILDNE3++te/4na7+eUvfxk225GEo8TEROAS8vPDp+oVFYmXlaw0+Kd9RKp0EvjLX3bjcgX2rKmIhJ/WxCu8/mj6/oi3xied7fOuu7JkZGQwezZ8+OGedmuuBstXrZo0aVKf7zNlyhQ6J16dK2a+TZSDNUVSRALHt34rKwtefvllbr/9Ddau/S4zZ87wTzXW2q6uTZw4h9zcZE466Rm7QwGiJPHy7VsyderUXm4Z3p5hMXPnHqS4uMzuUEQkyForSuFV8fI1KFLFq2vW9NDqdhUvX/Uo0F0Cu1uf1RMr8SqhoaGBlpYW74nJzolbamqqd6qkKl4i4c63fsvt/oDly5dTX38iLS2X+acZa21X94qLE1m16hKqqzezdetWu8OJpsQrgQkTet7nJNy1TibZ1cOtRCQaVFRUeJsBpdodSjtWa/Shqnh1o6zMOjHWtuI1atQo4uLigpJ4vfHGSEaMSO7z2ewpU6aQmGiSmLgTl2ufN1HsnLgZhuGdnv9VQGMWkeCwPusuITs7G1gUsc3k7HDhhRcCaTzxxBPeta32iYrE66abtjNixBRiYiL76bQmXlrnJRLtysvLSU1NBWLtDqUd64/5KFW8utFV4hUTE8PYsWMD3p59165d1NZO7NfZ7ClTpjB7NjzyyEZKSv7mvbTripm1dkx/b0TCXw333HMPMII777yTjttLSM+GDh0KfIfPPvuM/Px8W2OJ7EwFOHToEPv3F/KHP0T2NEOAEcCIESNQxUsk+lVUVITxouh0f4Ih7XWVeEHg98UyTbPb9Vk9SUtLIyUlhZdffhn4kHnz5gGHd3lba5bIXg4ePDjYcEUkqF6hqqoKuNXfAEn66xzGjRvHU089FfAOtP0R8YmX9YepKeIba/hYc/mVeIlEu/Ly8jBOvLLYuXMnhw4dsjuQsGMlXoZ/LZzPuHHjvFMNAzONpbq6mgMHDtDfxMswDK655hquvvpq4Al+8IMfAF1PSbIqXi0BnyIpIoHTBMBaZsyYgcs1VY00BiyO733vexQVFbFhwwbbooj4xGv79u1A5DfW8LESrxLb56CKSHBVVFR0+vAePg6nqakJdyBb9EUJa9PhNOLi4tpdnpGRQV1dHXAgIMfxNdbob+IFcPrpp3PRRRcBST3ebsKECSRSx6RJJQFtgy8igWOlCNVccskl/n371Eij/1wuOOOM03jrrSP529/+RpVNcUR84mUtNhwWBXt4WawzkDXs37/f7lBEJEjq6+upq6sL48TLmkEQLvuehBOr4jW60+Wt+2IFpnrUmnj1bQ+vgcjMzGQ2G1i5siSgbfBFJDBM0+RFAKZwzDHH2BtMhLOSVoP9+6+nsbGR/7UpjohPvKyK19So6e7i26+luLjY5khEJFh8jSvCdarhpEkZ5OYm893vKvHqqLvEa+zYsd7vvImXb+OdfnxlGUX+H78/50MSc3OBdIK1UU9iYiJpwO7duwP6uCISGJs2bcL633lJ1HzOtZvLlckLL3yXX3M27777bsiPH9GJV3NzMzt37gSiY5ohtO7X0nq2U0Sija9Ve7hWvNxug7vumkZV1Xa7Qwkrpml6X7sxna7rVPHybbzTxVeWy8TA+rft5W6yME2oqKikkh24fvlLwBjwRj0uV+/5WiZ49/oSkXBimibnnruSXH7MpEmn2h1O1Cgqgubmi6nnGJYvX47LddC/OXUoRHTi5Xa7ef31ZsaOjY7GGmB9EBs2LIE5czTnXiRaVVRUsH49HHPMqLBdJD1t2jTATUNDg92hhA2Px+PtANi54pWQkEBKSgp9mWrYuhkq1NXV8eGHH/Lxxx8DzXg8Hm699VZgv7dBxsD1ZT1Ia+KldcUi4SQ/P5+amu28zuO43eG17Uiki42NBX5AZWUlxcVv+8fjUIjr/Sbha9u2bdTXQ2Fh9FS8DMNg4cKJJCXt4je/sTsaEQmGsrIy6uuhsTGN+Hi7o+malXiZfPnllxx99NF2hxNcWVl9+qvb2mD/DquU1IFV85rUp2mBVpOO+/jOd75s00xpGDfckOjtZngPhx/edRv4QMoEamtrgWogJejHE5HemabJs88+C4zlbLuDiVKTSCE3dxLJyS9hmmfRXffXQIvoitfWrVuB1DZz66NDZmamphqKRDGrW+BY4sM16wL/Fh2OaLDhLUF1NwXQ93X62HfIZS7jx4/u8vqM3/8eKO3DtMAWHnjgAYYNK2Ht2svZvPkebr/9duAbJCcnezdIPSokT32C/ztNNxQJF/n5+d4eBpdFdoUkjLnJ5vXXL+KMMwr55JNPQnbciE687rprKykpR0XdgsOJEyd6F9/X2R2KiASBtTY12+4wejRy5EggwxmJl1fbKYBd2bu3jLlzYdu2zmu8wLfOq4zm5uZejrSaLVu2sHbttZjm5ezZcwKzZs0Cfs5DDz3Escce26f1WYHQ2jNRiZeIXXy9eLKyrGrX7NnPkJubwaRJZ9kdWlQ788wzGTFiBC+99FLIjhmxiVd5eTm1tft48snQnBUMJV+DDf0hFIlsWez0N6TzrdlsaGjgq6++AibbGVofTXNU4tW7fSQmJjJs2LAur7USL9Pb+bBr1lnsZzjttNM466zuP1SFar+e0eCtvKqzoYhd2p70+ec//0lNzZds2PBd3G7Vu4IpPj6eCy64gA8++AD4KiTHjLzEy3taYOtoa3Hz1y66qHNb3nBdrd5HrYmXphuKRDJfl7q2VRS32+1d0xPeFS/LNMrKyqiqsmurSXtkZlZiGEuZMKFj0lnG6NGju51l4etsWFrafYONJ554AhjJT37yk7CYrWFgbaSsE30i4eAjnnrqKeA0zjzzTLuDcYQLLriAmJgY4F8hOV7kJV7e0wJbH3kESCCrubnzXPsI39I7IyODuLg49IdQJHKtXLkSeKfT5dY0Q4iMitcRgLXewCl2797NuefexNy5b/PVV091uNZKvLrTW+JVWVnpXUtwLsOHDw9QxIM3YcIEhg/fHdKWyiLSnrVVxf1kZmYC14fFiRknSE1NZfr06cCbtLS0BP14kZd4eVmNNY7wtoSMLrGxsYwfPx5VvEQiU01NDU8//TSwhMcee4xDhw75r9u5cydJSUl01ZI8/EwjKyuLNWvWtOm8F80KWbRoEY2Njd5pgB912My+jDFjul7fBZCWlgYM6Tbx+ve//+39PZ4WyKAHLTMzk7POKqWp6WDIWiqLSFuHuP/++4FG73YSCXYH5Chnn302UOHd1iO4IjLxqq2tpaioCPia3aEEzaRJk4Aiu8MQkQForWody0svvcSvfvUrfPsk3XxzIStWTMblioSzmQYXXXQRRUVFIfmDZCdrv7L7iI+PZ+nSpd49tIawdu1aAKqrq4H9PVa8rDPUY7tNvK644h1ycyfhck0KePyDkZWVhWma7Nixw+5QRBxqlbeg8FNvxUtC6aSTTgKSeOONN4J+rIhMvD7//HPvWcPoa6zhY7Vy3sv+/fvtDkVE+qn1A+xirrrqKv773/8CBZimyf79RSxbNjliZkSffvrpjBw5khdffNHuUILq8ccfB0q58cYbGT9+PCNGjCA19Qx+9rM3MIxKUlJ+S2LiEGbMmNHLI2V0mXhVVFSwf/9Wnn32tLB77Y855hgMw+Cjjz6yOxQRx/nss8+AZzjrrLNwuc6MhlYFEcdqMHQq//73v70n4YInIhOvrVu3sn59LJMmTbM7lKCxNi91yB46IlHmyy+/9FZGRjBnzhxva/ZX2bNnD9DA5MmRsL7LEh8fz9y5c8nPz4+K/QV9bZvbdpssAF599VXgf9ptFp2ffyFz5zby/e/fwIUXfk5e3i/Jzu6tKYqVeHWcmvnuu+8CJqed1nmaYahax3dnxIgRZGdn8+GHH9oTgIhD1dfX8/vf/x4Yw7XXXhuybqbSlbNoaGjgP//5T1CPErGJV339FNzu6J0DO3XqVMBQ4iUSgXbs2MGUKVMAGDJkCOeccw7wvrdlLRGVeAGcf/75DBkyhDVr1tgdyqD52jbn529m8eKbcbt/yu/wdZP9XrvbZmdn87WvfY2qqioWLFjAKaec0ocjZFBXV8eBAwfaXfr2228Dk70dBNsLhw9bxx9/PJ9//jkQ3LO9ItLq6aef9m4/cWO321RIqHyNMWPGsH79+qAeJeISrwbw/nGI3vVdAAkJCYBLiZdIhKmvr+err77yJ14A5513HgDPPvssEONdwxn+fJWYlJSRbNw4m7y8PCZOrIz47nd79uxhyZIl3i5imZwF3HbbbUB8p9v+4he/4JZbbmHu3Ll9fPQM/zHaHs/6u/WNwYYeNMcff7x34+ctdoci4gjbtm0jNzeXCy64ADjS7nAEg5ycHD788EN27w7evoYRl3h9At4/DtPtDiUEDmPbtm0O6SYmEh127txJXp7J5ZdPweVtkGO1GZ9ObW0tkOmdTx7+fJUY04TKym/T0tJCSckL7fYlizwHWbp0KbGxsfzud78DbuGn0GUlCmDcuHF9rHT5jAd8nXct//jHP7yvec6Aow4ab3Z91AknMCQ3F/iw896Y/f2K5KxcJASam5t56KGHSEtL48orr7Q7HPEKxeyOiEu8vslv2bAhnmhurNHqcGpqanrcjFNEwsuOHTuor4eKiikUtdur63zvv5GwcXJnLlcGa9acwbBhr3o7/EWqJ/jyyy+54YYbGDNmjJV3YAZwjVUmxx9/PP/4xz8oKytjwoRCbrppI88/fxEuV1ogDhBY3uw63jQ56rbbgA87743ZzVeWy8TA+rfddZGblYuExKpVqygqKuL1168jKWmYmmmEAZcLUlNTePPNM8nLy6OmpiYox4m4xKuSYurqjqGrKSHRZvz4aaxfrwYbIpHE6miYQmpqaodrZnLiiScSztPNelJUBMXF85k9uyli13pZ0//WMGfOHGbNmgV48w6MAK6xMvjZz34GwF/+8he++upJLr88mQMHvh32C+aPP/54EhOLMIyqPhWtfOvllGeJ9N0HH3zAs88+y1lnnUVp6Szb13eKxTfDo6rqIpqamnjttdeCcpyISrz27dsH7Gb06BMccXaguHgS9fVD+eKLL+wORUT6yEq8pnj3dGorljvvvBOYZUNUgTFx4kROPfVUcnNzgcireq1btw4wmDdvXtCO4XJBRsZYNm78AQUFBUAB8+fP926aHd6OO+44Zs+GN9/8WMmUyEB0bJva4WuXYXD/SSeRvXYtP7vxRus+vU3ddcIH3rDi4oQTTiA3N9e7tCmwIirxsv6IwebN06Pz7IBvJbv3KzYuDpjKtoULNXdeJAI0NTVRXFwMTOl0nd0twwPlsssu4+DBg8CdQZuKEQzNzc3k5eUBMxk1alTQjuM7a1pWNodjjjkGyGDOnDlBO14gTZkyheTkZDZv3tzn+1hrkLUOWQRoLQN38VXt8fCbBQuIv+IKbtu3j3jf+v3epvVG5Qfe8HbRRRdRWVnJ2rVrA/7YEZV4WZuQjoreXb3brmT3fTGNwosvpvngQc3pEAlzbreblpYWukq8wqFleCBkZWVx2223kZhYzIgRdzJpUq3dIfVJfn4+VVVVwHkhOqLBPffcA/wpYpqpxMTEMHPmTO+2Bz2f6S0qKgIe59prrwW+F+Hr/kSC68CBA9x5552Ul5dz6623evd5lHA1ffp0Tj75ZJ544gm2bAlsp9eISbwOHTrERx99BJzQxRSeaDaNgwcPsnPnTrsDEZFefPnll97vIrOBRl/NmDGDjRtv5X/+p4hdu+6IiMrX66+/TlpaGjAzJMdzuWDIkDhcrsjam+eUU07x7kH2abe3qamp4aabbgLWMmbMGOAAK1euDFWIIhGlvr6eu+++m+LiYm699VaOOsoJzeEiQIdZZr4vF0XExBjcfvsNrFkzjqOPvo+Jxoc9Twfth4hJvLZt28aaNbWMZZzdoYTY14iJieHdd9+1OxAR6cUVV2wjNzeZSZPG2h1K0J144onceuutwE5uu+22sK54lJeXs2TJZp56KgeXKzYkx4zUCucJJ5zgrdC91+1t8vLyaGhoAB7g17/+NZDDK6+84l2HLSI+dXV13H333Wzfvp3//GcRM2fO0IqRcNHVLDPTpMjM8n6bhNt9K/PmNVDCk3iqqrqfDtoPEZN4bdiwgfr6eHZwqd2hhFgqs2bNYt26dTQ1NdkdjIj0oKrqC+6++3Dc7giryndz5q+3rxNPOgm4i90PPcQtKSnUhul+Tq+++ir19VBaem7EJUKhNnToUGbOnAn8p8s9JFtaWnj55Zc5+uijwbtdQmbm5bz4osERRzwT2mBFwlhtbS133nknDz30OWvW3MTQoV/3f0ZX34zIMGnSJBYtWgQUc+ONN+IOwHKfiEi8amtrefPNN4EzCP++UIHlcsFtt13ACy/sV9VLJIzV1dUBJRx++OF2h9J/3Zz569MXx3PXRx+xa+5cXl21qsf9nBoaGnjjjTe8zTlCo6amxtv+/lTGjo3+SmQgnHzyyUBll1uZbN68mb179zJ37lz/Zbt2jeJvf5tLVdWGgHwwEYl0NTU13HHHHd59HW/GNL/hP+nTdrjViaDwd+KJJzJ+/H08/ngzRx11Ex9//PGgHi8iEq833njDO63h/F5vG22KiqCl5Tjq68fz6quv2h2OiHRj+/btgMm0adPsDiWkXC447rhjefvt41m7dm2P7XdXrlzJAw88wB13hG5d2AsvvEBjYyNweUiOFw1OOukkIJZ///vfna7Lzc0lPT3dvw+az/z584Fh/P3vf1ePQ3G0SZP2M2LE7fzv/+70Tsf+ut0hySDt3n0YZWV/oK5uDHfffTf5+fkDfqywT7xM0+SVV17hsMMOAw6zOxxbWM1Ezuezzz6jyO5gRKRLvv32nJZ4+c7eVldb7XffeeedLm9XW1vLyy+/TFZWFtu2beOXv/yld0Pj4Kmurua669aydu1puFyTgnqsaGLtOXYs7733Xrvphrt376agoICXXz7f2zik9T7JycmMHn05v/rVZmJ4P1xmmYqElAfYtetWLrmkhPr6OzjxxBPtDkkCZNSoUUyc+DteeGESp5/+W956660up2P3JuwTry1btvDEE7v4xz+cV+1qbzbx8fG87P0pKytslk+ICL7EawLDhw+3OxSbzCAzM5OXXnqpyz9Gr732GmvX1vGXv1zPpk2/paamhiVLlgzoD1dHNTU1/PKXv2Tp0qWUlpYC1km7f/7zn9TVNVJc/B1N6emnMWNO59FH9xAT8zpZWdbarocffpj4+HjKy8/rcprUV1/N4dprM7mGk3C7QzedVCQcTJy4l1SuYdiwPdx5553AdLtDkgArLh7BgQO/ZciQqZxxxv2kpt7Z7+nVYZ94nXbayxhGEvv2nW53KLZyuZJ5880cXsfaz8zt1pZeIuHCNE1v4uWsald7BhdddBFffvklEyZstXprYG2D4XI1MW/eS8THH49pHsaQIUfy7LPX8Kc/7WDcuLcHddRDhw6xdOlSduzYwaZNm7juuuu45557+OEPf8gLL7wAnMHEiRMD8PycpbR0NnfccQIXX7wct3sHzz77LB9//DHXXXcdkNLlfeLi4rj66qux6phrQhesiM2KioooKVnEd3iUzZt/y3HHHWd3SBIkSUlJVFTcy4svXs3Bg9vJyvp5v+4f1onXli1b8Hje4eGHz8flGuro7i9FRTB8+I94k+9z//33A2rbKxIuysrKvO3UI7CxRoC4XHD++WeRm5tMXd1ztLSYuMkCoLj4DebOrWLjxnmAb+3qmVx/fTZ79z45qGYbjz/+OH/4w4e8+OJPef755bz44lncddcunnrqaHJzf8KkST8NwLNzHsMwuPHGGxk5ciRJSXfzve89R27uOZxzTk6Pf4tnzJjBiSNHAs9RPYBOmYHYJ0ckVJqbmxkz5iUmT76JpCSDJcARRxxhd1gSZHFxcVx00UWUlj7K449f3K/7hm3idfDgQf7yl78AY7jssssoKlL3F7c7gRru49ChQ8B9gNrLi4QD3/ouJyde1lqvoaxZczmnnfahtxMtVFZWAs9w2GGHceyxx/pvbxgGP/rRj0hM3Ed8/MsDmjr91ltvsWbNGurrL8I0czDNdEzzekzzUUzzJkzzfNzuhEA8PUcaOXIkixcvJienhp/9LIvGxmv71Intqi1bgCaefeSRPnXGbKiv57FHH+XTTz4JyD45IsHW2NjIW2+9xQ033EBZ2WPcddfX2LnzftquJPXt0uHkokG0S05O5kc/+lG/7hMXpFgGbdWqVZSUlAB3k5CgP5w+LppYufL/kZj4G+rrl9LcfDNxcWH7Moo4whdffOHddDbL7lBsN2fOHDZu3Mijjz4KHM1vf3sfUMfPf363t1FQq+OPP56bbprOtm0reO65k4GMPh+nvr6exx57zNt4qX9/+KTvjjjiCB588EFGjRrlfY/3bsKECcB5vPbaa3zrW99i/Pjx3d62oaGBX/3qV3z66afk5uZy3XXXcd555wUoepHA2rNnD6tWreLtt9+moaGBjIwM4HbuuuukTuOb04sF0rWwrHh98cUXrFq1ijPOOAOYYXc4YaWIyZjmLFavvg7YxIMPPkhLS4vdYYk42o03buP556fgcukkSExMDNdffz319fXAz9n2wAPAjUzOzu5yCtnV99yD8dxzwGKK+zH17Llhw6h6+mmuffBBIHZw09h0SrpHkydPJjk5uV/3ycy8nJdeGsKECU91u5d2Y2Mjv/71r9myZQs/+9nPePPN4/nmN//CqFF/D0jTFZFAKS8v56GHHuLaa69l48aNnHbaaWzefC+PPbYcl2tWp6RLpDthl3h9+OGH3H777aSnp/PnP1+jv4fduOCCC4Afctddb5Gefr93fYmIhFpRURH793/BE098TWc4vSZNmsSll15KYmItuazE5Tq52ylkE02Te4uKSEoCF5czYfxnvU5N211Swpr/+R9y/vhHpvk+oA90A2jtYhoUu3al8swzlzB37rts3fpZp2ZQpmny+9//nk8++YT33/9/fPOb5zFy5B08/PAFVFQ87930WsReDQ0N/OMf/2DhwoVs2LCBd96Zw/PPP8ovfnE98fFHY5qGhg/pl7A5PVtTU8NbU6bwWFkZmcCL/JdYqigiBdqeSFAm5udyfZvaWpPa2v8jJeW/jB79fXbvPpchQ4bYHZqII7S0tPDQQw8ByVxyySV2hxNWvvOd7zBr1iwmT57ca28El8vF9u33cccdd/C3v93M889fycUXX9zlWeS6ujqOOuqv7N8fz4sv/oAbbtCfhXB18cUX869//Yt7770XWErbqaRjxqygvPw/wDW4XGd6l3HFYprXcuutHv7nf/5GRsYE9uyZaU/w4limaVJYWMj69et58803qamp4Z13Tsfj+QEu1xgqKuyOUCKZLRWvmpoaMjLeYfToZ3nooYe4+eab+d73vsf/lpVx+OLF3FtTQwnHU2Rm6cxkD4qKoKxsHl988RC33TaVsrL/5Uc/+hFPPvlk0DcmFRF45ZVX2LZtG3BNv6diRTvDMMjOzu7zFJxx48bx4IMPArP4+9//zu23386mTZvadTzcsmULP//5z9m//1P+9a+rMc0U/VkIYwkJCdxzzz00NzcDd+LxeDBNk//85z+Ul/+DBx44m5aWC9u9foZhsGfP/8f110+mtHQp+fn5mnYog5KVRefprr4LDQPTMMg03scwNpJm3MVVMTHcMHUqry5cyHHPPsv9ubl4PDdhMpYit6Yty+AYAx3QZs6caebn5wPWWd/GxkYMw+Dww2OAGIqKYgFrDveBAwfYuXMnn376Kb/4xSc0N39Jfb1JYqJBTMxIFi4cw5NPnkBFxSxgKi6XgdutRkZdysrqcgMvE5jAy9SzgUY+op4EkpOn8dRTpzF16lQuuGAihjECt1vzkKV7hmFsNk0zok8xtx2bgqW+vp6tW7eyZMkSjjjiCH79619hml383zIMDWT95HKZFBe/xrBhT3P22TVs2DCMurpUkpLqOPtsD2PHjuWxx/4fpnmk3aFKH33xxRccccRtXHppLM3NzTQ1NZGbexiNjfd127CjvLycyZNv5cCBPYwYcTj79/9BY5MDtLS00NDQwMGDB70dnK0EPjExsU8ncUzTZP/+/ezdu5fS0lJKSkq44ooSbr+9hj//uZn9+5uBOEZygAU3JXLgwAEqKytZsaKWuXMhNzeJ3/3uOP74x+ns23cKYJ1Qc7l0gke615/PTgNOvEaNGmWed955NDQ00NTU2tY8N9f6d+5ca5G1v/HD+vWsrz+NIUwmnq8Rz1G4yWEqO3GThYsiilxnQlGR/6yE3uT9Z41L5fzzn2/xwx9upKamkMREqK8HMBg2bCgwhG9/O464OOsrNjaW2NhYDMPwf4E1gPneHx3fJ20HwLb3kcj25z//OeI/3IwdO9b8zne+E/DHNU2TQ4cO0dzcTGlpKS0tLQwfPpwHH3yQceMyus6vlHgNmMvVTHHxR6SmvsfNN9exeHESK1aM5sILL2TYsET9WiOMYXzKAw+s4ze/GUFlZToTJ55FcfHIHu/T3NxMXl4e8+atpLb27xqbooDvc4VpmjQ3N/vH1MbGRpqammhsbARg/Xrrc0tiIsyebX3OSExMJCEhgaFDhxITY03Yamlpobm5mebmZhoaGmhoaCAvz/TfNyfHYO3aMdx4Ywrx8fHExsZy6NAh7r33ILfeGs/w4cNJSUnhJz/JZseOw5gyxYVpxmroln4JSeJlGEYN8EWvN4wuo4Byu4OwgROftxOfM8DhpmlG9Jw5h45N4Mz3rBOfMzjzeWtsilxOfL868TmDc593n8enwTTX+CLSzz71l2EY+U57zuDM5+3E5wzW87Y7hgBw3NgEznzPOvE5gzOft8amyOXU96vTnjM4+3n39bZh105eREREREQk2ijxEhERERERCbLBJF7LAxZF5HDicwZnPm8nPmeIjucdDc9hIJz4vJ34nMGZzzsannM0PIeBcOLzduJzBj3vXg24uYaIiIiIiIj0jaYaioiIiIiIBJkSLxERERERkSAbdOJlGEaVYRjrDMNYFIiAwpVhGPMMw8gxDGOB3bGEilNeW/C/vuu6uCyqX/NunnfUvO7R9Fy644T3aVec8Nr6OHF80tgU+aL9PdodJ7y2Pk4cm2Bw41MgKl7zTdM8xzTNpQF4rLBkGMY8ANM087w/59gbUchE/WvrY5rm6rY/O+U17/i8vaLpdY+m59KJU96n3Yjq17YtJ45PGpsimxPeoz2I6te2LSeOTTC48SkQiVeKYRjZAXiccHYiUOj9vhCYbmMsoeSE17Y7Tn3NIbpe92h6Ll3R+9SZnPq6R9NrHk3PpStOfY9C9L+2PdHr3otAJF5pQKVhGMsC8FjhKqXDz+l2BGEDJ7y23Unp8LNTXnOIrtc9mp5LV1I6/Kz3qTOkdPjZKa97NL3m0fRcupLS4WenvEch+l/bnqR0+FmvewdxvT2Kt2yY1uHiQl8Z0TTN5d7beQzDmNdN+S3Seej8O4h6Dnltu+PBga85RNbrrvFJ79Mofm174sGBr3skveYam5z5HgVHvLY98aDXvcfXvdfEq6c3jHfhXL5pmgUDjjQyfEBrFp8NrOv+ptHBQa9tdxz3mkPkve4an/Q+tTsWmzjudY+011xjk/Peo+CY17Ynet17Mdiphiu9B/QtpovKrN77vLK9iwRTfGesopwjXlsf72s7s8PzjfrXvOPzJrpe92h6Ll1yyvu0C1H/2rblxPFJY1Nkc8J7tBtR/9q25cSxCQY3PhmmaQY/QhEREREREQfTBsoiIiIiIiJBpsRLREREREQkyJR4iYiIiIiIBJkSLxERERERkSBT4iUiIiIiIhJkSrxERERERESCTImXiIiIiIhIkP3/rSp8aHjex9EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "kde = ConditionalGaussianKernelDensity(\n", " whitening_algorithm = \"ZCA\", \n", " bandwidth = \"optimized\",\n", ")\n", "kde = kde.fit(\n", " data_xy, \n", " features = [\"x\", \"y\"], \n", ")\n", "\n", "mini, maxi = -5, 15\n", "X = np.linspace(mini, maxi, 100)\n", "fig, axes = plt.subplots(1, 3, sharey = True, figsize = (15, 5))\n", "epsilon = 2.0\n", "for y, ax in zip([-5, 1.5, 5], axes):\n", " # plotting histogram of actual samples\n", " ax.hist(\n", " data_xy[\n", " np.logical_and(\n", " data_xy[:, 1] > y - epsilon / 2, \n", " data_xy[:, 1] < y + epsilon / 2,\n", " ), \n", " 0\n", " ], \n", " # bins = 100, \n", " density = True, \n", " label = \"samples\", \n", " histtype = \"step\",\n", " color = \"red\",\n", " )\n", " # pulling conditional samples\n", " data_x_cond_y = kde.sample(\n", " conditionals = {\"y\": y},\n", " n_samples = 10000,\n", " keep_dims = False,\n", " )\n", " # plotting histogram\n", " ax.hist(\n", " data_x_cond_y, \n", " bins = 100, \n", " density = True, \n", " label = \"KDE samples\", \n", " histtype = \"step\",\n", " color = \"blue\",\n", " )\n", " \n", " # calculating analytical probabilities along the axis\n", " x = np.stack([X, np.ones(len(X)) * y], axis = -1)\n", " probs = np.exp(kde.score_samples(x, conditional_features = [\"y\"]))\n", " # plotting the function\n", " ax.plot(X, probs, label = \"KDE analytic\", alpha = 0.7, color = \"black\")\n", "\n", " ax.set_title(f\"$P(x | y = {y})$\", fontsize = 20)\n", " ax.set_xlim((mini, maxi))\n", " ax.set_yticks([])\n", " ax.set_xticks([-5, 0, 5, 10, 15])\n", "plt.legend();" ] } ], "metadata": { "interpreter": { "hash": "58ae7b3716633041a01a40c36e0d7b7460b562af78f6efaf769f0d9c360a3f13" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.13" } }, "nbformat": 4, "nbformat_minor": 4 }