{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Clustering" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clustering techniques are unsupervised learning algorithms that try to group unlabelled data into \"clusters\", using the (typically spatial) structure of the data itself. It has many [applications](https://en.wikipedia.org/wiki/Cluster_analysis#Applications).\n", "\n", "The easiest way to demonstrate how clustering works is to simply generate some data and show them in action. We'll start off by importing the libraries we'll be using today." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import math, matplotlib.pyplot as plt, operator, torch\n", "from functools import partial" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "torch.manual_seed(42)\n", "torch.set_printoptions(precision=3, linewidth=140, sci_mode=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_clusters=6\n", "n_samples =250" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To generate our data, we're going to pick 6 random points, which we'll call centroids, and for each point we're going to generate 250 random points about it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "centroids = torch.rand(n_clusters, 2)*70-35" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from torch.distributions.multivariate_normal import MultivariateNormal\n", "from torch import tensor" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def sample(m): return MultivariateNormal(m, torch.diag(tensor([5.,5.]))).sample((n_samples,))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "torch.Size([1500, 2])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slices = [sample(c) for c in centroids]\n", "data = torch.cat(slices)\n", "data.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below we can see each centroid marked w/ X, and the coloring associated to each respective cluster." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_data(centroids, data, n_samples, ax=None):\n", " if ax is None: _,ax = plt.subplots()\n", " for i, centroid in enumerate(centroids):\n", " samples = data[i*n_samples:(i+1)*n_samples]\n", " ax.scatter(samples[:,0], samples[:,1], s=1)\n", " ax.plot(*centroid, markersize=10, marker=\"x\", color='k', mew=5)\n", " ax.plot(*centroid, markersize=5, marker=\"x\", color='m', mew=2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABKCUlEQVR4nO2deXhV1dm37wUhM2QihCkxzFNkkEFAlIA4UFFSWqx+0iJqoQ4t1dY69X2Rt9VSrVpr1eKEtlqsCEbrLEqYEcIcmSExhEDIHJKTmfX9sc/e2edwMkASkpw893VxnXP2uDbKbz37t571LKW1RhAEQfBOOrR0AwRBEITmQ0ReEATBixGRFwRB8GJE5AVBELwYEXlBEAQvxqelG2Cna9euOjY2tqWbIQiC0KbYvn17jtY60tO+ViXysbGxJCcnt3QzBEEQ2hRKqe9r2yd2jSAIghcjIi8IguDFNJnIK6U6KqV2KqU+dv4OV0p9pZQ67PwMa6p7CYIgCA2jKSP5hcB+2++Hga+11gOAr52/BUEQhItIk4i8Uqo3cAPwmm3zTOAt5/e3gISmuJcgCILQcJoqkv8r8DvgrG1blNb6JIDzs5unE5VS85VSyUqp5Ozs7CZqjiAIggBNIPJKqRnAaa319gs5X2v9itZ6jNZ6TGSkxzRPQRAE4QJpikj+CuAmpVQa8C4wVSn1NpCllOoB4Pw83QT3EkpyYePzxqcgCEI9NFrktdaPaK17a61jgVuAb7TWc4CPgLnOw+YCHzb2XgKw62346n+Nz7qQzkAQBJp3xusS4D2l1J1AOjC7Ge/Vfhg5x/WzNszOAOCKhefuL8k1jhk5B4IimraNgiC0GppU5LXWSUCS83sucHVTXl/AEGRPou1OfZ1BfZ1AQ5HOQhBaNa2qdo3QhNTXGTT0jaA+mqqzEAShWRCRb41cjOi4oW8E9dFUnYUgCM2C1K5pjdgHVxs6gGo/rjGDrud7rtlZiFUjCK0SieRbI/bo2G6HmL/tEb4Z9Vc4YO2Smms0xELx9MYg9osgeBUi8q0Ru5ViF/ytrxhCXuGAKY8Y201RnvwQXPN/rraJJwvFLuyeBN1+PxlUFYQ2j4h8a8fFO9dunxgCXOEwttnF2DynJNfoHNAwbkGNsFeUAAom/tI4vyTXONd+v43PS1QvCG0cEfm2xLgF4BtUE22bkTYa1v7Z2GcXd3cbp9IBnQJh8sPG901/g77xcCwJfANr3hbMDkEGVQWhzSMDr20J90FOy25R51o11j4Nl1xhbPsu0egM0IbYA0QNrznXtIPW/hk2Pgcr74DiHHDkNc2griAIFx2J5NsS7h65PdJ298zt+ypL4fuNUHjc2HZ8K0x/uiZ6N8+tdBifl1wBp1KMCP9YEuQcgMNf1lzbbuGIby+0UvJKKliRfJzZY6IJD/Jt6ea0GCLybQkzOk/bAAn/8JzrbhfdKxZC9mFyt39IBBji3bETHFsDhz4ld/AcY/vG543jOwUY14idBHGz4bPfGpH+ZT8ztnka1JVsHKGVsiL5OH/67AAACyb3a+HWtBwi8m2JkXMMgT/8JST+okbo7dg7gmuf5PG5k3l1wynW/HooA2/+l2G9fPkohyp7MWVwP37+g9E83ifZOHfcAudFFASGw89sNeUC58DWpca+cfNr7iu+vXCRON/IfPaYaBwVVTgqqskrqWi30bx48m2NyMEQ3tcQek+TpUbOgQHXwuEveXzBTBZ/dpLMM5opL6ZxaO92+Oy3HNryGVN++FMycwpZ/M9veHxPLxj4A+P8zJ2GL7/yTljzZM11d71tePVrl9RUwBSrRriImJH5iuTj5+zLK6lg6dqjHM0uZunao5aoB/r68PzXhz2e016QSL4tsXWpkREDhpDXNlnqyt/x+Hu7WPxRzZK7mXklTPnBTF65XjH/4zIyz9SkYS7+YD90/gWPz7/J6DzC+xqWzrE1NRk7I+fUpF26WzV2+0gQmonZY6JdPu2YHcCWY7msOWisMLdgcj+Xc9qrRy+RfEtSX6aKuT/7sPFZUWpsD+8L1z5pfC/OMdIgB/7AEt3c12fz6pqjLpdKIAFHvj8zlpfiOBNCgtuSu69+kkxu9A+MTJvBM4yN5nU3Pm/8nvKoMQnLbtX0jTc6hq2vNPZvQxAumLGx4fSLDGLuhFgemT7YEvfwIF8WTO5HeJDveb0J2N8I2joi8i1JfQuAmPu/fNT49A0wIvi8Y3DoU2P/5r8ZGTCHPjVEt+sAIlQ+a+4bQM+QToAh8AtZyHM8RyyxPMdzLGShJfQ9u3RkzfuvE3H8U4ieCKf2wNifQ/R4SFlh3DvxF+d2RkEREH2584dGEJoaU3zNKNwUaft2gBe+OczR7BLe2pzG7DHR1jF2Zo+JdukA7JjX/uPH+1w+vcHmEbumJalv0NLcPvAHRnaL6Zv3vMzY58iDo2uge5yxb+tS6BMPQd0Y2HM0a64dwZRZPyUpP4mZzCSWWJaxDIA00kgiiZ5dfFjzMz8GnviPEZEHdgVHDhQch7yjhth3HVAzBuCeQeM+QUsQ3HC3Sc7HNrHbML+fMRTAEnF75szvZwyloiqFAVGdeWtTKs9/fcTaZ2JG9Z7aNm1oFADThkYxvm+W9empQ2hriMi3JA0t9xsYbhy38Xlj8LNvvLG90mH45gAp7zsnOmEI9fcbGdg3nleuV8xYXsBiFlsCD7CYxRRQwNsP/4yBV11mRPDZB6Hgewi9BPpNMUQ+9zDkHK4ZA4BzB1zPtwia0K4wBdlRUU2gb0ccFVUeRdgTs8dE23z2fTxz80jCg3zP8ef7RQZz1cBI/vTZAS7vE86k/hGWcJvklVTw1qY0QDN3Yh8XC8feln6Tg10+m5KWGBcQkW+tlOQaFsnhL40BT98gI1o3UyiPJdWI/bE1ED3OKFK2+11DqMP7cehMZ+Z/XEYooSxikcvlF7GI+7mf+c98wJoptzFw/VPGeQAdfGDcLyA0xrjnoU9rPP/aCpvVhuTRt3tMIXZUVPGnzw6w8OoBtdomdgxRTiU6PJBLwgNZczDbEkhPQunaIcDqfVn0mxxsCaujoprnvz4MQKCvzzkDsxeDlsjdF5FvDZjRrimoppAe/tIQ8uNbayL2hH8YtkxFKSig+6XQKQh6jYXPHwS/LgAcCr6cKb9ZRuYZTQLxxBJLGmksZjGLWEQsscQTT2JuIlOm38SaOT4MvOwqOHPSiOC/fNS4l0nK+zU1cMzOxrSP6kLy6AUnN43shfE/rfYYybpHuSuSj1sRP8CUQZFc3SfQo1AeST/J16kOfj9jKMN7Z1JaUUVucQXPfXUIgOe/PszCq/uz8OoB1v3hXAunuSPti92pQBOIvFLKH1gH+Dmv977WepFSKhz4DxALpAE3a63zG3s/r8SeimiWDzBF0SwwZrdLfIOMz7V/rqk78+I4w0sHcgP6MOWxVWQWlAOQSCIASSRRQAH3c78h8M7tmQXlTPnXWfbcdT8R/S+reYMwB4S/+l+4ZCL0nVLTER3+0hgniLRF556smaZagUpos9hFOdC3I3/67IAVSdsxvXRHRRX3XzPImsyUX1LJgVNFpK9+i0mPfsC8J99g4dUDmDY0iqVrj9LPr4iZ06+j07Bp/PmJP3D/NQNZuvaodc+FV/e33hxM4TYHbqcNjWL1viyXjqU5I21P4wLNTVNE8uXAVK11sVKqE7BBKfUZMAv4Wmu9RCn1MPAw8FAT3K9tUpc37T7Aah4z8AdGaYExP4egcOMYq378wzUCv+ttQ+D9QgCIKE3l58N9WXyq5haJJNKzR3fe/uXPmP+nf5F4JtGlCT8fqYjY/XcY8aERwZttBdfO58tHa9I33aNzsWYED5hR67ShUXy0K5OFV/e3BNo+GLvlmJG9VVp5FjAEce7EPvzmvV18+fbfKdy4HICnf3Ub0x96GYC/rEiiYMXvcRRkw8bl7PggmqND/4d1h7L52YRLCAv05aaRPVm9L8tqT15JBb95bxdrDmbXmVfvLTRa5LXWGih2/uzk/KOBmUC8c/tbQBLtWeTrEkB7tGuPjL981PDei04Yg5+ZOw2Bda8fb05U2rvCSK/0CeDxCaVQ7svitUYaWc+ePVnzySoG7nmSNXMDmfKWw5oQtejGPjx+WS50iYYv/wfStxj3L86BSb827lldAdWVNRG8JxEXa6ZdUp/FYUavS9ceddomA/jfxBQ2Hs1l3aFsXvh/l7Ei+Tjfpjpf9DU899UhSiuq2XAkm00r/mEJPEB1cR6f/flujv34N+R98BxlhTnWvleef4pvU/MoGJLA/lNFvPazsda9cosriAj2xVFRxZqD2UwZFMkvpw4AsAZpzUFdM+PGHuW3VZrEk1dKdQS2A/2BF7XW3yqlorTWJwG01ieVUt1qOXc+MB8gJiamKZrTOrkQATQj5it/B+ufqhFY30Bn3rytfnzmTkPguw6A7qMg5T0evz4KRl3Hq+99ypo1axh44n04lsTAERNY80gFU57Zyc/n/4LHrwk1rJ+i47DrXzX3z9gKwV2N78eSjLeHgde5Ztm417ORCL7dUZ/F4Z6m6KioYuNRI2rfeDTX2vfFdyfJLCjjw10nyDpjWI3VpUUU7/7C5XoJJJBUnMS+Nx8jlFCuJ8GyHgGObfiI7v2nkVcSzIPv7+ZodgkA+04WsuFIrjXwO21oFH/8eB9rDmYzvm+WlU1j2kbrDmVb7WzLBc6aROS11tXASKVUKPCBUiruPM59BXgFYMyYMd47o+ZCBDByANy2wvhuDbg6IO7HxjZ7eYHDXxq+fcI/YONfje2XzeXxa//AL3+fS0REBJxQxvaOnRjYK4g9d2girgkDtCHgcT+Gd35kZNn4hcClN597P7vVZNazMXEvXSy0C+qyOI5mFzP/n8mW0C6Y3I+8kgpKK6rZnVHIiN4hVuS8I73Q5Vy/jorygC5E3fonspY/QnVxnjWxbyYzXZIIwLAkO0d0o/OsP1DRKZh+kUH8zw1DeXX9MfpFBuPfqSOjLwlj7sRY8h0V3PnmNtJyHVzeJwxHRZWtiJnx72RYrxCuGhjZ5q2bJs2u0VoXKKWSgOuBLKVUD2cU3wM43ZT3apdk7jTE3DfQtcMws12ufdIQWLNkcKcAKMkl4oAzcwddsxIUEDHieqxVpa75P2ensrLGdz/0qfONwVw1yi1qt9ezQYsf3w6pz6r548f7OJpdQnhQJ8bGhls+fICvD9+m5tGpYwfyHRWcyHfQK9SfEH8f/Dp1xNenI9+m5gHQKbwXNz26lM+fuoekgton9nUMDueKXz5P56gYqs5q+kcG8+gHe8ksLKPq7Fm+Tc1nyqBIF4EH6NSxI89/fcQaDJ47MZZA345t3qYxaYrsmkig0inwAcA04M/AR8BcYInz88ParyLUiz1ad7d83LNd7LNQ3TN3+k4x1nXVyiiTEDfbdcaq/e0h0GYx2aN2s5MJijDq2UBNyQP7erGC12O3ajzlr/9+xlDS84xI/tmvDrHhSA4n8ktZd9gY7NxwJMeyTAAGDorkmZtH8tamVCqrq9mRXshlMSFcOWAAk599gV/fcUutE/u6X7+IsSPj+OdmY77HjvQC65jKas2UQZHWpCpT4C+LCaVvZBBaaxdfvi3bM+40RSTfA3jL6ct3AN7TWn+slNoMvKeUuhNIB2Y3wb3aLw1dBQpcrSFzW/REyD3qnDh1uauvb4q43f5xn83qqQqlnaAI41pm5C/RfLvAbtV4nD0aGcyKX0xkRfJxcovL2XAkh3WHs0nLdRAbEcjMkb2YPDCSkvIqlFL8fsZQKz9+Un9jPCivpJK/rEii5IP/rXNiX/bnL1B8zeXERgSSluuge4g/pwrLANBaMyCqM8N7h3DTyF4M6JbO7oxCQFudgjl5yttoiuyaPcAoD9tzgasbe33BSV2efkP2bXzemOQEWJk5cO4MVvBsu9ijdjv21FDJrml32KPe2rx585i8kgoigv0YGxvOC98c5vczhtIvMpila4+yNS2fR6YPpl9kMAyFLcdymTshluN5Dg4fPkz+it9TWpBNAgm1T+wrTuQfv/sZYTc/waTRl6K15lRhGdFhAew8XsjO44VMGRTJ3Il9LLsI4PI+4YzvG97mvffaaDczXqvy8ylctYqQWbPwCQtr6eY0nvOtCTNyTk3q5bgFnqN9T8v71Yd7aqhE8O0Gdz++rgJg5jHm/mduHmlUeBwKjopqK3f+ua8OseVYDt+m5lNZfZZjJ05Zg65Q/8S+0sIcKpY/wtmh/2bbqSqmDIpkQFRnXll3jLDATlZpBHvV1PF9I7j/moENfs62RrsR+cJVqzj99F8AiLjzzhZuTSOx17UBV2F1F3/77ymPeL6e+5vA+Qi1RO/tlobMDnUvTuY+s9ScjPTI9MGs3pdl1ZYBw0eP6dGNMyOuc8mTTySRLhHd+Pitj7nr5/NJPJnocs/Lp89myW1XWDnuAIezzrDmYDaT+kfgqKhm8sBuJKflM6xXCHMnxjb6OVsz7UbkQ2bNcvlsc9jr23z5aO2DsO6rNdU3CzX7cE02TeSA82+X5Ma3WxoyO9S9OBkYQjltaBTrDmUTHRbI/Cv74qio4qaRvXBUVJNfUs66wzl8m5pHzxB/QifdBmAJfcfgcK769QtMiL+GtUlrmBw/hVMnMwEIueJW/C+/hX6RwYSN8bUicPPNwaiAeZg9GQVsPJrLmNjweic+tfVZsO1G5H3CwtpmBG+Ke0WJkd1iZsmYOfGeBmHNY+yDqbWVCTY7DKjJqhGEBlBXFord4jD9+EBfH2sZvj9+vI+NR3PZeDTXynoJ9PVh7sRYfvPeLtJyHfSLDLLy66+dcx99rurHW8tep8/cP7O3OIjZ/9jE/9wwlAm/fJ4Nz/+KIfEJqNGzefrHIwDPbxFgVKA068WbnY97eYOGPmdboN2IfJvFvVaNe30bd4IiXGvPuEfa7mvCRg4yyhWYs2sFoQlwtzjsQrl07VHWHMwmNiKQa4d15ydjo62FOsyaMlMGRfL7GUNZ9OF3bDiSg69PR+bf/xBHIyeRnFVFF38fjmaX8Iu3kymrCuCqh5ZxzWX9rUj8aHYx6w/nMP/KPoB2aYu9brzZ+XjTIiHuiMg3A006yOspdTKyHnukLgvFPatm0ws1E6EEoYmoy+Kw132PCPKlX2Qw/SYHW+I/xZkrHx7ky99uHcWvlu9kw5Ec5wSmKgD6dwtmR3oBZVWa8KBO9I2OsqL2+68ZyB8/3seGIzl06qh45uaRgDHAWzOr1cCe+eOtiMg3A00xyOvSUTSl511fVo07srKTcAHUZ3EM7x3KcGdJAxN7x2AXYqMGIqTlOpjUP4LRl4Rz08iefLQrk+S0PDYezeVEfql5NIBzqcB91pKBezIKnZZQR4/tauuDq3UhIn+eNCRKb4pB3mbPBmqoeEv5YKGJMSY7HeaR6YM9RtXux248msu42DB8fTqyeOYwI5ceuP+agS7Fz+zZNP0ig1k2bxyAyxtCbXZMWx9crQsReeoWbvd97uLr6dymGORt9myghoq3pEgKTcz5CGpNdo6xdJ/7rFR7x1DbbNXa3hDstPXB1TrRWreaP6NHj9YtQc5rr+l9gwbrnNde01prXZmXp3Nee836rG2f/dzv58+3trUJinO03vBX41MQWjm5xeX6H0lHdG5xeUs3pVUCJOtadFUiec6Nmu3RurkteOpUcl9/nZBZs1yi9JBZsyjZupWStesoXLXKZV+rnmUr+e1CG6K2SLutz0a9GHRo6QZcTKry88l9/XWq8uteajZ46lQCJ06gKs+YSh1x550Uf/MNp5/+C4WrVrlcr3DVKqIefpiu997LWYfD5dr5b7/N6af/Qv7bbzfPAwlCO8ccMDVKFQieaFeRfG2DmfbtwVOncnzBL6hMT8exaTM+4eFE3HnnOdF+VX4+mQ8/TMnadQB0CAzg9NN/oUNgoO3ayvps1VG9ILQCLiQq9+YB06aiXYi8KbDBU6cC5w5m2i2Z4/feS2V6OgABY8ZY+3zCwgiZNcsS6sJVqyhZu46gyVe5XM/+PWzObXQIDLCOr2/AVhDaMxeSxujVA6ZNRLsQ+fy33yHnxRcp3rCRXs8+c46omtkwua+/TuWxVJS/P7qsjKDLL8cnLMwS5JzsbM6++Rb5K1fS88kn6XrvPYAiNy+PiPBwAKoLCqwOpfibbwieOtX6fdbh4Kyj1LpebSmS0gEI7RGJypuHdiHyZ0uNVWAcmzefMzhqJ2TWLIo3bMSxeTOBEybgP2I4R37wA4ImTOCPzz/Pqqoq/jliBL2PpZLz8ssEjRvH1ieeZN7DD/HTqyZzR2qqNQjr/pm/ciVB48dT8O/lgCZszhzrnu54VcVMQWggEpU3D1478GofZO0QEAhA4IQJ54iq/TifsDACLzPWPwm8bBRZS5ZQeSyVJ5Yu5aXcXE4VFvKzHTs4OXAgnXr15uCBg8w7dZIsh4O/fP4ZrwQF0vXuuwm/8w7OOhyE33kHXe++mw7hYVQeS6Vkw0YAHDt2AoaAe4rUQ2bNotuDv227FTMFQWg1eG0kb4+G7d64u6iax511GNOiK06dotMlMQRdeSX+I0bwyM0382JWlnX8qeJi/t/qr3h8xw4WZZ3idFWVte+vO3agb5/H/aNGUbotGdWpE2X79nM2L59OMTH0fOrP5Lz8spVuaff47e1qsxUzBUFodXilyFfl53PW4aDrvfdaAmqKprvfbUbLZx0Ocl58ybpG9vN/o+BsNStyc12unUACSaVJ3H0ig1BCSbCtSgOw/LsUFsycSZBvJ/z69yfv9TcACLpyEoEjR9JzyRIXgT/99F8o2bqVnkuWiP8uCEKT02i7RikVrZRao5Tar5T6Tim10Lk9XCn1lVLqsPPzoilY4apV5Lz4Eh0CA2qN3AtXrTon6yb0tv+H75DBdOjSmeozRfh+u5U3+/enm4/RFyaQwEIW8hzPEUssz/EcC1lIAgkARAUE8O9rryU8OIieS5YQcddddLokBoDiDRvOyc8PmTWLoMlXUbJ2HSd+8xuyX/h7vTn8giAI50NTRPJVwG+01juUUp2B7Uqpr4Dbga+11kuUUg8DDwMPNcH96qWuui/BU6dSsnUrwVOnumTdODZvJmjyVVTsN1K4ylO+o0N4GLF5+bw9cSK3Hz5M0skkZjKTWGJZxjIA0kgjiSS6+fqyrHsPeh49Rt7RY5Tt30+vZ56h55//zPF77qHq+3SjY8nNI++NN6jKzSPi53cREHcpurISx6bNODZtpkNggFg1giA0GY0Wea31SeCk8/sZpdR+oBcwE4h3HvYWkMRFEvm6PO3ib76hZO06iseNs7JufPvEEjzpCoKnTiWvV2+K160jcNw4zpaX41i/nqGTJvHyPfcw85ZbWMxiS+ABFrOYAgp4ObI3sb6+dOjShbNFRTg2bSb31dcoP3qEs3n5BE2+iuCpU0n/+XwASvfscb5xvEjXe+8lcNQoQMlgqyAITUqTevJKqVhgFPAtEOXsANBan1RKdWvKe9VGfTnm9ig//+13APAJC7cmKFWeyKAqI4Oq6GgcmzcDsGf5uyx44QVCCWURi1yut4hF3M/9LMo6xZvhQ5j8ylKyn/8bjs2bObN6NZXp6QRNvoqohx8ma8kSqjIyjBOVcpmc5amtki8vCEJjaTKRV0oFAyuBX2uti5RS9Z1injcfmA8QExNzwfc3BdEcQD3rcNAhMJDgqVMp+vgTzpY66BAQSNic26wo38y6CZ46lewX/o5j5w4cmwzbxq9ffxybN5PeuTO379nD6YpyEkgglljSSGMxi1nEImKJJZ54EqsSuSM1lQ9WrmT4wl+RmXmCyu/TCRg7loC4OIo+/oSStesInDABMHL2i7/5ptYMG5B8eUEQGo/SzlVXGnURpToBHwNfaK2fdW47CMQ7o/geQJLWelBd1xkzZoxOTk6+oDbkvv46p5/+C13vvZcOgQGcdZSS8+KLdOrbh8pjqdZxXe+9hw6BgZaoutegCZp8FT2XLAHg6NKlTH78cbLKy63zE0ggiSQKKCCUUEPgbdk13Xx8+GL2bPycufA+vXpRdeIE4XfcgU9EOCGzZlFdUMCp//sD/kMGg1Lkvf4GXe+9l8hf3ucSvQMuHYBE9oIgeEIptV1rPcbTvkZH8soI2V8H9psC7+QjYC6wxPn5YWPvVRd2G8YUxNKUvVb07D9kMB0CAjlbVuqM9EuJ/OV9Vg2aDqGhBF89Fd/u3a1rBqSm8qPgYF6yiXwiiUT5+fHPH9/Gb//zHxKrEl3a8eOQEKKGDKHCzx9dWUmps9PqEOBvReOFq1bh2LzZmFk7cYLzTG3ts0fvtRVSk8heEISG0BR2zRXAT4G9Sqldzm2PYoj7e0qpO4F0YHYT3KtW3AdbfcLCXHLSzcg3+4UXnEcYohoyaxb5K1dSeSyV0p07KTqWSofAQM46HJSsXcfvbrqJUF9fnly6FIAoPz9W3H473Y8c5c3oGG4/nm5NiHogPp67O/mCNuyYrvfeS9Dll2MvY2Desyo3j9K9e/Dp0YPACRPoMmPGOfn97jT7alGCIHgdTZFds4GamrruXN3Y6zcGT1k2YXPmWHaNSef4eErDI/AbNBCf6dNdBmUDR13GgpS9FEREsLK4mGU9e9E1aS1VQN9u3XgTuP14Oj+bdg0LFFQeS6X80CGCJl9Flxk34NenD2CMGRgdjCJszm34RIRTui2Z0m1GpF/8zTcA5Lz4Et0e/G29yxAKgiA0BK+c8VoXpvCXp6Zy4je/MSyVbckETb6Kgnf+Tdd776Vw1Sq6zLjBKoXQpaCA31VUck9UN7pGRHC2rIyK1FQif/UrwtdvYF1JMV0qKin497/x6d2bgBHDyXv9DYrHjcPPZtHYZ9SGzbmNs45Sa0DY9OrNHH53xKoRBOFC8HqRry0CzlqyBMcmI0XSTHEsirvUyrCBGjE1C5c5XnwJ7phHxbFUHJs3UzLqMiJ/eR+RwPfz7gDANzqaiLvuwic8/Jw688UbNzqvbVhFHQIDjKje2S5zfMDeOdjPt38KgiA0BK8X+fy337ZSKiN/+Utre9TDD3OqshL/IUOIuOsuAEpT9loplHYxrcrPtypHlu0/YOXPm2IN0P1//4esJUuIevhhj+3wCQuj1zPPOJcCVNZsW3vdmrqEXIqWCYJwIXi9yNuX4DOpys+n+Jtv6PXMMwDO/PpSa6UnM4XSXLjbzIYxI/6CFe9Ttn8/XWbMsK7p16ePNdBrpm/CudZKaUoKJWvX0fXee6y6NWaNexFyQRCaGq8XeXuZYRO7vw048+vvsWq4+4SFWXn3cG56pk9EuDWZyc9DiqP9WnbsSwaGzZlDGLjkxAuCIDQ1Xi/ynqLjkFmzrKX4usy4wdpm9+zdhd1+jdpsFfdz3PG0335dyaARBKGp8XqR94RPWBgdAgM5/fRfaq36WJd1Utu++uyW+vZLBo0gCE1NuxR5aJ3ZKq2xTYIgtG2apHZNU9GY2jVCw8gvyyfxSCIJ/RMI8xdLSBC8gbpq13jtQt61YV+4u7VwMduUeCSRZ7c/S+KRxGa/lyAILU+7s2suhu99vgOoF9OLT+if4PIpCIJ30+5E/mL43ucr2hfTiw/zD2Ne3Lxmv48gCK0D8eSbgeZKhRQ/XRBaBkdRId8lrWZY/DQCu4S0dHPOQTz5i4yZKtnUue7ipwtCy/Bd0mrWvbOM75JWt3RTzpt2Z9e0ZcRPF4SWYVj8NJfPtoRE8m0I008Xq0YQLg6OokK2fbQSgLE3/ahOq8Y81lFUeLGa1yBE5Nsx+WX5LEtZRn6Z59TN+vYLQkvRlIJa17XOx6ZprZaO2DXtGNPjBzxm3NS330QGhIWLjSmoYETYzXWt87FpWqulIyLfjomPjmfbqW3ER8d73N/QMYCGdgaC0FQ0VlDt2TJ1XSuwS0iDO5HzOfZiInZNOyW/LJ8nv32S9SfW81nqZx6tGU9jAJ6OS+ifwAOjH3DpDMTqEZoTU1AvNJ3Rbq3Uda3W6rOfD00SySul3gBmAKe11nHObeHAf4BYIA24WWst/+KbgfOxS8xjS6tK2XJyi7W9odG4p+M8TbCS6F642JxPLntD3gQcRYV8/tJzpO405u60xii9ITRVJP8mcL3btoeBr7XWA4Cvnb+FRlBbdOwpf762Y99IeYNntz9Lflk+d4+4m1sG3cLO0zu5rNtlVjTufq79t6eo3RMNOU6ifaGx2CPthg58NrQz+C5pNak7k+kzakyT+Owt9VbQJJG81nqdUirWbfNMIN75/S0gCXioKe7XXqktOvbknS8/sJyXd7/MpsxNPHXVU9b5KTkpAKQVpfHY+Me4Z/U9bDm5haySLJ6f+rwV5b+8+2VKq0oJ8AkgryyPN797k7yyPML9w637LEtZVuvbQ0PKJ0i0LzSU2oTZPmhqCnHf0Zez7aOV1m/7efbo/Pi+vVx/z/0ux9i/D4ufRmVZGVrRJDTlYPH50JwDr1Fa65MAWuuTSqlung5SSs0H5gPExMQ0Y3PaPnYxd7do5sXNsyLjhP4JlFaVArDl5BbmfDqHiT0n8u7Bd7ll0C0AVJ6tJLUwlQfHPkjGmQxSi1JZ+M1CUotSuXvE3dw94m62ndpGclYyY6KM2dIH8g64WDyNFWiZ3CU0lNoE0i7spjjbjwVcOoFP/vY06Xt3Edq9J6k7k/n8pefo3ncAm1cuP+f4sTf9iE7+/qx7ZxlZRw9z/T33Wx3FhZQ46Dv6co7v20vf0Zdf+F/EBdDi2TVa61eAV8CoXdPCzWlV1OW1m1HwtlPbeGLSE4T5h1nbSqtKOZB3wDo2/Uw6ZEJVcRVh/mEE+ASw/sR6nt72NE9MeoJJYZOICoxiy6ktXNnrSqb3mc7T254mOcvwIivPVjK+x3jmDp1Lpw6drGycujJzGoIUSxPqwxRUUxjdbRNz0HTbRys9RvQp33xJzKUjrE4gfe8uAPqPGU/uiXRSdyZTXVnFiGtvIG33TibePIfxP76VyrIyHEWFDIufxvF9e0ndmcx3SasZe9OPzjsidxQVsvOLj8k8uI/0vbuJHnopEb16N9VfUb00p8hnKaV6OKP4HsDpZryXV+JuZ9h/J/RPYNupbaw/sZ7EI4nMi5tnpUSag6rje4wntkssaUVpdPyqI5v+uYmxq8cyfex0HFUOiiqKuOPtO/j4wY+5dMal3P2buy2BX39iPTGdY4gMjGR71nYAOnXoxPoT6xnbfSyA9b1PSJ+W+QsSvB53QTV9bXsU7SgqpLKsjPE/vvWciD7541UAHNy4lpHXz6CivAylYeT1Myg9c4aCrFOkp+ziTF4O+ZkZnMnLYfCEK9m8cjmd/P0Ze9OPuP6e+62OZttHK2vtcOp6hi3vG28KTeXvnw/NKfIfAXOBJc7PD5vxXl6Ju51h/wzzD+OJSU9YkT5A0vEk1p9YT1zXOJdBz9sfuJ2PXv4IgGuuvoaFry4EP9i6dyupS1KpKqhi59s7+SToE/iFId4xwTGkn0nnhr43MK77OACm95nO2O5jXewVsVqE5sQ9C8ZTFP1d0mo2r1zOVbfN49j2b10i+rQ9O0nfu4v9G9Yy6IrJXDH7NuvaSe8vJz8zg5hLRzJx9hy++Mdfyc/MoKD4DFfdNs+6p/m2sHHFO2x5fzkV5WWMum5GrZaNu50zLH6aS+dysatYNkmpYaXUcoxB1q5AFrAISATeA2KAdGC21jqvrut4S6nh86EpZoua1xgQOoCntj3FH674AzFdYkg8ksj7f3+fT5d+6nK8T6gPk++fzNrn1lJVUOWyb8qdU8i+Mpvbh91uDbJeSLtkFqzQHNjtm2Pbvz1nsNT83nf05RzYtJaqsnKOJG+h4FQmfUaNsaLyxM3beGXpK9w5fjjTb/0Znfz8qKgoJyPzFPcv+Qvzbr+dhPGjqSwvp5OfHyOvn8Guzz9m88rlxFw6gsjYfmz/7yp6D4kjeuilLuJtWkdX3Tbvog2w1lVquKmya26tZdfVTXF9b6YpMkzMa4T5hZFfns99X99Hv7B+fHv0W06uOulybAIJJBUk8fWirwkllHjiSSTR2v9t4rcsmrOIAJ+AcwTaFO746HiSjifVKeBmdk9pVSn3jLzngp5LENypy4N3329aJCOuvQHVoQNDJk3l3UUP8e6X3/DVvsMAvLJ+O1H9BnL2dCbZZ4p5deMO8orO8PRzz7Fr6ACuixsIQCd/fwZdMZn9G9aSvnc3lWXlAGTsTyFjf4pl7djb01rKG7T4wGt7x92KcY+AGyKsdn8+oGMABRUFbM/azpRBU5j39Tzip8TjyHWQQAILWchMZrKYxSxiEbHEApBIIn5hfiz/73JyO+fy7PZnCfAJsDqe/LJ8HtvwGOtPrLfuBZL6KLQMdiE1I+zKsjIm3nybtd20SCrLy8nPzGDTin/xn9VrLYEHyC8u5vHX/8mChBks27KRvKIz1j7zuJ/fcjN9R1/O2n+9RkFWJgAdfAzpNCN5Tx1Na0FEvoVxzzBxH1x1F1Z7No39GqY/f1m3y/j7rr8zOHwwd8TdAcD8l+bz0i9eIik/iZnMJJZYlmFEQWmkkUQSPqE+xDwUw3/P/JdHhz5KaVUppVWlbMjYwFPbnmJk5EhrMHbB8AXEdY2jtKqU/LJ8qy32DurWwbdabwOCcCF4SlW0bzOF1MxjT9+3F73iHUZdZ1gnpv++6b13AOjcuw/bvv+3yz0SSCCpNImnl68klFASSHB5s915MoeRP/wJa//1Gqk7k4m5dAQ9Bw1l8MTJll1kplW6Dwi3FkTkWxn2yD7xSCLrT6znyl5X8uDYBwFjUHTuZ3N5furzADy97WkWDF/AjtM7rCj/1Wtfta730q6XWF26mh5ze/D9X79nMYstgQdYzGIKKGDwXYPx6e7DlpNbSDqeBMDLu18m1C+UgvIC8suNWanpZ9LZcXoHAT4B50T77mmdEuULjaG2QVb7ttwTGWQe3EePAYM5sT+FE/tTrJx28/hBV0wG4MTB/Sy4aiyvb9lDTn5+vW+2YZ2D+fjDRE7v32MIfNxIbvjVg5aI29MgzXaZE6xak9CLyLcy7JG9p2yauZ/NtSYuRQREkJyVzPdF35N+Jp1tp7axYPgClu5ZanUKnxz9hPJT5Zz+52lCCWURi1zut4hF3M/9HHntCDOemsGkUZPIK8tj26ltAAyPGM7e3L08Ou5RjhUdc2mXGe2bZQlKq0oZ32M860+s57ENj53zxiEI54Mnb9t929p/vUb63t2Edu8BQGhUzSSnLpFR7P7yExyFhQSGhJCesosBAwby0CV9+cNb/yaptPY325DAAOZPGs2u5ctQzjeFyD59Aaw0Snskb8+n3/nFx/j6+Z9zTEshIt+K8TRZaHL0ZKrTq0ktSiU8IByA0d1GA0aUvzt7N0UVRVRWV1JaVcrhI4dJ+3MalfmVXMd1xBJLGmkukUs88SQWJPLNY98w/NXhrCpYZd0v7Uwa+eX5nHSctAZQTVsGjGg/wCfA+n73iLutfHozf18Qmgp3v3vyT+8C4PIf3kLmwe/oOWgYn7/0DKk7k+kcGQVAxoF9dPTxoUu3KApOZaKAH4++lDc2JNf6ZnvHZWOI7BxMcW62y/3NiH3bf1dR6qxBY1awNDN3KsvKrKi+NRQ3E5FvQyQeSeTN797k7hF3A0bkbOawpx9NB6CoogiA1KJUKs9UWnnwgOU1JpFEAQXcz/0u2TVF2UX89ed/peeinvSO6k3v4N78+rJfW1aQvR3Pbn+Wu0fczQOjHyA+Op7PUj/j9mG3A/Dg2AfPyacXhPOlITNLAzp3JnropfgHBwPw/e7tFJwyMsoCO3emS9eunMnNoeh0FgAxcSMp6eDDB589U+eb7fvb9/L4yFEMGTqMI9s2U15STNaRQ0T1H0hoVE8KsjIJ69mbYfHTyD2Rwdp/vcbkn95lTdjq5O9Pz0HDAKzJUxdaDqGxSD35NoRZ2dEc1HzzuzcBI7KOCowiKjCKyIBIALIcWfSI7EHklEiXaySSCN3gkl9fQmVEpcsgE0C/6/oxqf8kTpWc4qreVzGi2wiXmvL5ZfmUVpVy94i7uXXwrcyLm0fS8SRe3v0y36R/w8u7X2bV4VXUh1SgFOx4qtA4LH6aS4kBT5gdwdp/vca6d5ahlSHkAFnHjtCxYyeKTmfhHxTMiGtvYOANP+R3z75A/pli4om33mznMY800qw326LScv7w1nL270uhvKSY0O49yNifwvb/ruKSEaMIjepJTNwIAL5Z9g9SdyaT+NT/4SgqtN42Mg9+R+rOZI5t/9alrRd7eUCJ5FshtU0k8uTXl1aV8u7Bd61jegT1sL6P7TGWPy37Ez8K+hE7394JQEBEAL965VdEXRJF75t7c/OMmynPM3J++83uR+W0SmJDYkFBXlmeS/YMGFH8y7tf5oHRD1jb7Smc4FrIrDa7RipQCnY8Re2BXULw9TMKhJ06dtjjgKbpzZtRc+zw0SgNYT17kX8yg7AevUhP2UVZSTEVSjF9xo1k5xlzMut7s80/U8zjr/+LZ+5bQFBAAAWnThJz6QgCO3ehICuTgi8zyc88QVjPXqTv3UXBqUyrvo2jqJCSokJi4kaeUwbBm8oaCBdIQwTQXnmytKqUz1M/55TjFJN7T7bE99bBtxLmH8bXr37NTwJ+wsZVG4n/YzyrClZxe6/bye2cS8zvYkh/Kp0uV3Zh0h2T2J61nU0nNpFenM6Wk1sI8AlwmRjlvmSg2SE9OPZB4rrGAUb5AzOnvzakAqVgpzYB9FQgzI598pPpf5s14NP37qbnoKF0HzCIU4cP4qs1N0yexGvvrrDOTySRLgF+3Dvlav69/lsSzyS6XP/yPjGEdAkmfe9uYuJG0mvgEAZdMZkTh4xiY+kpu+g5eAgTfnQrWrmWX9j+X+ON9sCmtfj6+bukfV5MRORbIfUJYH5ZPssPLLfqvd8Rdwd3xN3hUsdm+YHlLD+w3Co4dvKKk/zk+p+w7YyRNXMg7wBPXfUUpdeVkjM6h65du7oUJxsTNQaNZuuprVaBMtOasRcms3dI9pmt9RUtkwqUgp3aJhDZBzSHxU/DUVTIrs8/RitcctXtlSejh15Kz0HDqK6qpLSoCEdBAQC5x79n3g9nUl1VxbL3PwAgLDiY+VeOZvT48QTrKpau20ZBiQOAa4YO4I4fz2LqvLs5tv1bKsrL2Pz+ck4dO8zE2T8FFN0u6Wvl5dux16JXmnrHFpoTEflWiLsAuts3pmViYuaqm7n15qIfADtP77QqUt438j7+uuOvaDT3jbzPypBZcWIFV2KUGI7rGmdNdDI9/yt7XVlnsTT7pyA0NWaKopm5YtZ+zzp62CV7xRTQiF692fbRSiPS3rsbgIAuIYRE9WDzyuXcNftWevQfyKuvvMp/Ez+gQ342fUdfTsdOhhy+8e1eZkydzLAOFfToP9DqSErPnOHAhiSjPHFVFel7dxE7fFStg6id/P2tzsf+/WIjIt8GcLdvzEVB8svySStMs6wTs17M7cNud8nA2XJyC6O6jWLH6R0kZyXzwOgH2HF6h5Uhc2WvKy0/ff2J9Tww+gHr3uN7jK81312KkAkXC9OzH//jWy1rZPDEyeeUFICa0sOjb5yF0nDyyCFOHEihMMvIutEKnvjTEh747YNERERY53TvO4Ab+w3g3pcmEhEe7lID3sTM3Ol2SV9ih4+qVbjdxxgkhVKoEzNKjo+Ot1Z+umfkPSxLWca7B98l6XiSiz0S4BPgktNuX7LPvM5nqZ9ZGTKAVR/HjOKn95nusUiZvcOBxq8OJQgNwe7Z2yNnT4tv2EsPj73pR2x67x1OHEihqrKKXkPiqCwvx1FUSERERM2CHgf2k56yi6tum0d0rPFvydfPn/S9u+kzaoxV1XLMjFn4+Pu5WDT21Ejz/udbc745EZFvA5j2zbKUZedE9FDTCXiqF+Nu/ZjXcc+QsR9jWj2eqkfaOxyzoxCrRmhuzqfol/sg7qArJnNg83pOHTFWSzuxP4W8jHTL669tQQ+7z2/WrjE7DnutmtqWG2wtRcpE5NsQ7qLuLuANHcxsiI++8/TOc9In7faMp1RKQWgNuHcIx7Z/S35mhjWJyVzf1Yy+K8rLqCorx8ffDzi3CJqZuWPvBDwtIO6pg2gNiMi3QhqSJ9/Y69Z2nVsH30pKTorHsgTuFTLtn4LQWrFH5Me2f2tZL5VlZQBcMfs2qz69r58/4BqNe7KK3LfZO5XWEsFbaK1bzZ/Ro0drQes39r6h496M02/sfaNFrptXmqff2PuGzivNa9D2872O0D4oKSzQWz98X5cUFrR0U85py9YP39d/ufkGvfJPi3RJYYHL/tra3Zqexx0gWdeiq1LWoBVili+wLyTS0BIAdR3rft3ajjffGNxtmNq214YZ+ZupmkL7oqWm8TekLcPip9Fn1BjLtjGj8cAuIS7f67qGHU9lGVoLYte0QupaSKQ+u6auYz3ZPe414JvSXxdLp31zsafx2ydKuU9Qcm+L+ySrhlDX8zSkmFpL0SQLedd5A6WuB54HOgKvaa2X1HZse1zIuyGcTz76+eau25f1e2D0A5IKKbRZ7Ou+XsxFtKHlKkya1LWQd7OKvFKqI3AIuAbIALYBt2qt93k6XkS+ZZBJTYI3UFck7+3UJfLNbdeMA45orY85G/IuMBPwKPJCyyB1ZARvILBLiLWQt1BDcw+89gKO235nOLdZKKXmK6WSlVLJ2dmuq7AIgiAIjaO5RV552ObiD2mtX9Faj9Faj4mMjPRwuCAIgnChNLfIZwDRtt+9gcxmvqcgCILgpLlFfhswQCnVRynlC9wCfNTM9xQEQRCcNOvAq9a6Sil1H/AFRgrlG1rr75rznoIgCEINzT4ZSmv9KfBpc99HEARBOBcpayAIguDFiMgLgiB4MSLygiAIXoyIvCAIghcjIi8IguDFiMgLgiB4MSLygiAIXoyIvCAIghcjIi8IguDFiMgLgiB4MSLygiAIXoyIvCAIghcjIi8IguDFiMgLgiB4MSLygiAIXoyIvCAIghcjIi8IguDFiMgLgiB4MSLygiAIXkyjRF4pNVsp9Z1S6qxSaozbvkeUUkeUUgeVUtc1rpmCIAjChdDYhbxTgFnAUvtGpdRQ4BZgGNATWK2UGqi1rm7k/QRBEITzoFGRvNZ6v9b6oIddM4F3tdblWutU4AgwrjH3EgRBEM6f5vLkewHHbb8znNvOQSk1XymVrJRKzs7ObqbmCIIgtE/qtWuUUquB7h52Paa1/rC20zxs054O1Fq/ArwCMGbMGI/HCIIgCBdGvSKvtZ52AdfNAKJtv3sDmRdwHUEQBKERNJdd8xFwi1LKTynVBxgAbG2mewmCIAi10NgUyh8qpTKACcAnSqkvALTW3wHvAfuAz4F7JbNGEATh4tOoFEqt9QfAB7XsewJ4ojHXFwRBEBqHzHgVBEHwYkTkBUEQvBgReUEQBC9GRF4QBMGLEZEXBEHwYkTkBUEQvBgReUEQBC9GRF4QBMGLEZEXBEHwYkTkBUEQvBgReUEQBC9GRF4QBMGLEZEXBEHwYkTkBUEQvBgReUEQBC9GRF4QBMGLEZEXBEHwYkTkBUEQbJQWV7Djy+8pLa5o6aY0CSLygiAINvZvOsnmVUfZv+lkSzelSWjUGq9KqaeBG4EK4CgwT2td4Nz3CHAnUA38Smv9ReOaKgiC0PwMmdjD5bOt09hI/isgTms9HDgEPAKglBoK3AIMA64HXlJKdWzkvQRBEJqdgGBfLrv2EgKCfes9ti1YO40Sea31l1rrKufPLUBv5/eZwLta63KtdSpwBBjXmHsJgiC0NtqCtdMou8aNO4D/OL/3whB9kwznNkEQBK+hLVg79Yq8Umo10N3Drse01h86j3kMqALeMU/zcLyu5frzgfkAMTExDWiyIAhC68C0dloz9do1WutpWus4D39MgZ8LzABu01qbQp4BRNsu0xvIrOX6r2itx2itx0RGRjbuaQRBEGiYV17XMRfitTf2ns1FY7NrrgceAiZrrR22XR8B/1ZKPQv0BAYAWxtzL0EQhIZieuVV5dX4+HVkyMQe5wykmscA50TjDTm/tnt6ut75HNPUNNaT/zvgB3yllALYorX+hdb6O6XUe8A+DBvnXq11dSPvJQiCABgR8f5NJ13E177N9Mgry6vZvOooJw7lM+32oS7HVpZXM/aGWI9+uvv5ULsom/ftM6Kry7meaAkPv1Eir7XuX8e+J4AnGnN9QRAEd0qLK1j95j7SU/KAGvF1j5Ivu/YSSosrOHW0kPSUPHZ88T0BnX3pM6IrG1YcJj0ljwmz+nmM0MtKKjlxKJ9L43sTExduCbgn7PcdMrHHOZ2PnZbw8Jsyu0YQBKHZ2b/pJOkpecTEhbtExJ6i5IBgX7r3CyHjQD7HdmVTlF3GiUP5Hs+3Y3YCRdmlFGSVkjowh7Brgzwea7+vKfgnDuUzafYAUnfnNMjqaU5E5AVBaFPYRdUunvYo2W7dDJ/Sm9PfF1nCPmn2AFIH1i2+Y34QS1F2KVfMHkBeZkmd9or9vkMm9iB9Xy7pKXmsrTrIiQMF51hFFxupXSMIQpvCLqq1ZarYJykFBPsy7fahjL0hlm6XdME/qJMVddeW5XLySCEFWaWWwO9Zk8HW/x6rNysmINiXnv1CAYiM7kxMXDjpKXktOllKInlBENokppBXllfTyZYBU1pcQVV5NSOviaayvJrS4goCgn2pqKhm91fHyTiYR9dewaSszaT0TAVX/GgA4Hng1uwMkj9JA7AybfauyUADw6f0PidCv3RKb+s4s532NwFPg8bNiYi8IAgtxoUInns2S5UtA2bIxB7WoKwZRZsdwNHk0wCcOlJEcV45AIe3ZXHsQAbTbh1F8qdppKfkWWmT3Qb5EhDsy5CJPSg9U0HO8WL6jOjK3jUZbHOKfie/jucMpJrnmM9l31/boHFzIiIvCEKL0dC8cXtn4H5OaXGFFTnbB2XH/CCW6qqzlJ6pYM+aDIrzy63rBYX4UpxXznurX2XT/k948OjzhPr2ICYuHA0kvpHEP756iJnX/oS7/t9COvl2JONAPqm7c6ioMLLBg8L8cJypID+rhENbs1AYUTzAl699R8aBfKrKqxl3Y1+X561v0LepEZEXBKHFaGjeuHuaov3TfeATsNIkTxwo4MSBAsbeEMuwyT05/l0ePQaGEBDky6fb/8ln2/8JwNP/WciTC18jrEc0B/Yf4O+f/pb8Mzm89u7fOHG4gEWLFjFhVj+GTOzBji++B6Akv5zdXx0n/2SJFZn7+BnFdjMO5AM1tVw85dJfrIFYEXlBEFqMhuaN28W7IXnoW/97jPSUPHoOCKFbbBcqKqopzCqlKKeMLl0DePZvT1kCD1DoyOV/XpzPzRN/zfJ1z1LoyLX2fbb9n3T/VxcW/PR+ADr5GkLec0AIPQeGMXBcFIcuyaKqoprK8moGjouiqrwaDQwcF8WOL7+nsrya5E/SrPGDi4lk1wiC0OoxxTt1d46VNVNbHZjS4grS9xkiHd4ziIDOvsaA64F8ukT6k5efw6b9n7ick0ACZ4sq+cfnj6Ec1SSQ4LJ/1SfLSVq5mz1rMhgwLoqYuHDi5wxm+JTepO7OYeC4KPJOlpD8SRqpu3MYd2NfLr+xr9XeqopqYuLCqaqorrf9TY1E8oIgtBmGTOxBZXk1VeXV7FmTYWW92N8G9qzJICv1DAAFWaWMu7EvjjMVpO3KoTC7FOjIkwtf4/GlvyA7L4sEEljIQmYyk8UsZhGLiCUWgEQSCQmM4Fc3PkOwfwgKOLQ1i/SUPA5dYvjw2z5J4/i+PDIO5BMaFUCP/iHs+PJ7lyydqvJq0lPy6HZJF8v2uVh1bETkBUFoMwQE+9LJryObVx1l7A2xlmCCEcHvWZNB5iHDDw+JDGDEtGhWv7mP8B5BFGaXEhzhR3FuOd3Do/nkv59zw03Xk5SbxExmEkssy1gGQBppJJFESGAEj8x5gb59+9EtJoSKimpOHS0EIPNQPt1iuwDQNTqYDj6K9JQ8vv3oGCcOFHB8Xx7X3jXsnMFh02a6WHVsxK4RBKFNMWRiDybM6selU3q7TIoyI/vMw4X0GhxKn5FdWf+fQ6Sn5JFzvJgJs/rRe1AYYEx20rldePieJyiggMUsdrnHYhZTQAF33fgwwSoK/yBf/IM7sfur42QdKwIg83ChZcP0GhTG2SrNyGuiiYzuDBiDr3vXZACelxS0b2tO60ZEXhCENoW7YJq2x6mjhYy4JpqxN8TSs18ou746TlF2GX5BPgSH+1FZXk3h6VIA/II6cujwQf7vmYcIJZRFLHK5xyIWEUoor/13CaU+2Ub+fEU1vQeHMXhiFJ0j/AHDDkpPyWPjisNkHMgnbU8OQyf1pNfgUABOHi1sUL365lxGUEReEIQ2h10kh0w08tszDuQT2NmXcTf2ZcC4KHoPDqNH/xDKS6o4sOkUyZ+kEdWnC6FRAaSfSGPh/82l0JFLPPHEEksaacxjHmmkEUss8cRT6Mjlz8t/RfiwKrKPnyHjQD6OokrO5JYRExfOuJv60HtwGFH9uhASGWAUM9udw3V3xVlt8iTc7qJuvp00h3WjahZzannGjBmjk5OTW7oZgiC0cnZ8+T2bVx1lwqx+ludtLx9gTkYacU00OcfPcOJAAd37dyF6UDhBPeGqq8eTV5htXS+BBJJIooACQgklnngSSbT2h4dE8ruEpQwd08eluqR98HTMDbHnlFeoLd2zqUsbKKW2a63HeNonA6+CILQ5PA1aVpZXs3dNBpUV1dZkJAX06BdKWFQgaXtz2fZJGn6BHbm8/3SXPHkzi+aZ37zGE397hMTCRJf73T73Dq654TL6jOhqCbx5z8ETu5N5uICYYeH06BtqnVPXHICLWVdeRF4QhDaHu0jai4j16B8CQK/BoVRWVLPrq+N0jvC36tWUO6q5dfoCuvftwrIVfwega3g3Fi9Yyq33XsO4ySO4MWE6BcU5AEwf/TNunjafy669xHqDMEn+JI3QqACKsstI/jSNG+8beRGe/vwQkRcEoc1j1nE/caCAbn260GtQGAo44Uyn9A/2AfzxDeiIr78PU346mGG7e3HqWBHfHvmMxQuWovO6cGhrFvnf+/DLG/7Ci58/yMxrf8Lo8B+ibPcBY+btoa1ZjL0hluhh4SR/msak2QM8tu1iV510Rzx5QRC8Ak9FzEZeE03eyRLCewSx66vjLhORzPo2+5JTCfYPMQZq+4WwzRmdj7+lF+HhEWxYcZhJswfgH9SJPWsyUBg1aZI/SbPGBOrCjP5j4sKbbfEQ8eQFQfB6PBUqM6Pn/KwS8k6WuJQKriqvZtrtQ6muOmsMzPYL4dIpvTl5tJCMA/nkHqmi0C+H9JQ8Ugca1o1pCblPxKqLIRN7WEsO7t90sm2t8aqU+gMwEzgLnAZu11pnOvc9AtwJVAO/0lp/0ci2CoLQTmisxeHu2afurhFr07vIPFrApVN6c91dcS73MteE1Zw7wFtZXk1VRbW1ryFtM1emcl885GLR2Dz5p7XWw7XWI4GPgf8FUEoNBW4BhgHXAy8ppS5u6TVBENosTT05yJ6HPnxKb2LiwjlxoMBaHtA+81RhpEOaqz6Z+wKCfbn8xr4EdPYl+ZO082qbpxmvF4tGRfJa6yLbzyBqyifPBN7VWpcDqUqpI8A4YHNj7icIQvugqeu6uEf2tUXW+zedZJvTa68tt91eE742Wnqw1U6jPXml1BPAz4BCYIpzcy9gi+2wDOc2T+fPB+YDxMTENLY5giB4Ac2dR17b9evqXOxevn21J0+CfrEqTDaEeu0apdRqpVSKhz8zAbTWj2mto4F3gPvM0zxcymMaj9b6Fa31GK31mMjIyAt9DkEQ2iFNXdirNlultLiCk87qk+5C5slaas4yBedLvZG81npaA6/1b+ATYBFG5B5t29cbyDzv1gmCINRBfRFzU9km+zedJONAPjFx4Qx3ruNq4in6v5gzWuujsdk1A7TWh50/bwIOOL9/BPxbKfUs0BMYAGxtzL0EQRDcqc+7byrbxD0l005rEnRPNNaTX6KUGoSRQvk98AsArfV3Sqn3gH1AFXCv1rq6kfcSBEFwoT6BbaoB3NYu5HUhM14FQRDaOHXNeJV68oIgCF6MiLwgCIIXIyIvCILgxYjIC4IgeDEi8oIgCF6MiLwgCIIXIyIvCILgxbSqPHmlVDbGpKqLRVcg5yLerznxpmcB73oeeZbWiTc9yyVaa4/Fv1qVyF9slFLJtU0gaGt407OAdz2PPEvrxJuepS7ErhEEQfBiROQFQRC8mPYu8q+0dAOaEG96FvCu55FnaZ1407PUSrv25AVBELyd9h7JC4IgeDUi8oIgCF5MuxR5pdTTSqkDSqk9SqkPlFKhtn2PKKWOKKUOKqWua8FmNgil1Gyl1HdKqbNKqTFu+9rUswAopa53tveIUurhlm7P+aKUekMpdVoplWLbFq6U+kopddj5GdaSbWwoSqlopdQapdR+5/9jC53b29zzKKX8lVJblVK7nc+y2Lm9zT3L+dIuRR74CojTWg8HDgGPACilhgK3AMOA64GXlFIdW6yVDSMFmAWss29si8/ibN+LwHRgKHCr8znaEm9i/H3beRj4Wms9APja+bstUAX8Rms9BBgP3Ov879EWn6ccmKq1HgGMBK5XSo2nbT7LedEuRV5r/aXWusr5cwvGQuMAM4F3tdblWutU4AgwriXa2FC01vu11gc97Gpzz4LRviNa62Na6wrgXYznaDNordcBeW6bZwJvOb+/BSRczDZdKFrrk1rrHc7vZ4D9QC/a4PNog2Lnz07OP5o2+CznS7sUeTfuAD5zfu8FHLfty3Bua4u0xWdpi21uCFFa65NgCCfQrYXbc94opWKBUcC3tNHnUUp1VErtAk4DX2mt2+yznA+NXci71aKUWg1097DrMa31h85jHsN4JX3HPM3D8S2eY9qQZ/F0modtLf4s9dAW2+z1KKWCgZXAr7XWRUp5+s/U+tFaVwMjnWNwHyil4lq4SRcFrxV5rfW0uvYrpeYCM4Crdc1kgQwg2nZYbyCzeVrYcOp7llpolc9SD22xzQ0hSynVQ2t9UinVAyOSbBMopTphCPw7WutVzs1t9nkAtNYFSqkkjLGTNv0sDaFd2jVKqeuBh4CbtNYO266PgFuUUn5KqT7AAGBrS7SxCWiLz7INGKCU6qOU8sUYOP6ohdvUFHwEzHV+nwvU9vbVqlBGyP46sF9r/axtV5t7HqVUpJlFp5QKAKYBB2iDz3LeaK3b3R+MQcjjwC7nn3/Y9j0GHAUOAtNbuq0NeJYfYkTA5UAW8EVbfRZnm3+AkfF0FMOOavE2nWf7lwMngUrnf5c7gQiMzI3Dzs/wlm5nA59lEoZdtsf2b+UHbfF5gOHATuezpAD/69ze5p7lfP9IWQNBEAQvpl3aNYIgCO0FEXlBEAQvRkReEATBixGRFwRB8GJE5AVBELwYEXlBEAQvRkReEATBi/n/6HhTc/bAFm4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_data(centroids, data, n_samples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean shift" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Most people that have come across clustering algorithms have learnt about **k-means**. Mean shift clustering is a newer and less well-known approach, but it has some important advantages:\n", "* It doesn't require selecting the number of clusters in advance, but instead just requires a **bandwidth** to be specified, which can be easily chosen automatically\n", "* It can handle clusters of any shape, whereas k-means (without using special extensions) requires that clusters be roughly ball shaped." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The algorithm is as follows:\n", "* For each data point x in the sample X, find the distance between that point x and every other point in X\n", "* Create weights for each point in X by using the **Gaussian kernel** of that point's distance to x\n", " * This weighting approach penalizes points further away from x\n", " * The rate at which the weights fall to zero is determined by the **bandwidth**, which is the standard deviation of the Gaussian\n", "* Update x as the weighted average of all other points in X, weighted based on the previous step\n", "\n", "This will iteratively push points that are close together even closer until they are next to each other." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([ 9.222, 11.604])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "midp = data.mean(0)\n", "midp" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABKsElEQVR4nO2dd3xc1Zn3v8cqVrOaLcuyLFnuxhYuWDYugI0xhBbwOiFvCOzalIVQEjbJsiGQTT5v3hQSAlk2ARKWYjYQWIoDLIRginvBlhuWbdmyLFnVssqozkijkc77x517dTUaFVuSJY2e7+ejz8zccu65lN957u8857lKa40gCIIQmIwY6A4IgiAI/YeIvCAIQgAjIi8IghDAiMgLgiAEMCLygiAIAUzwQHfAzpgxY3RaWtpAd0MQBGFIsW/fvgqtdYK/fYNK5NPS0sjMzBzobgiCIAwplFKnO9sndo0gCEIAIyIvCIIQwPSZyCulgpRSB5RSH3h/xyulPlFK5Xg/4/rqWoIgCELP6MtI/iHgmO33I8BnWutpwGfe34IgCMIFpE9EXik1AbgBeMG2+WbgFe/3V4DVfXEtQRAEoef0VST/H8C/Aa22bYla61IA7+dYfycqpe5RSmUqpTLLy8v7qDuCIAgC9IHIK6VuBM5qrfedz/la6+e11hla64yEBL9pnoIgCMJ50heR/DLgJqVUPvAGsFIp9SpQppRKAvB+nu2DawkNlbDjaeNTEAShG3ot8lrrH2mtJ2it04BvAp9rrW8H3gfWeg9bC7zX22sJwMFX4ZOfGJ9dIYOBIAj074rXx4E3lVJ3AQXALf14reHDvNvbf3aGORgALHuo4/6GSuOYebdD5Oi+7aMgCIOGPhV5rfVmYLP3eyVwVV+2L2AIsj/R9qW7waC7QaCnyGAhCIOaQVW7RuhDuhsMevpE0B19NVgIgtAviMgPRi5EdNzTJ4Lu6KvBQhCEfkFq1wxG7JOrPZ1AtR/Xm0nXcz3XHCzEqhGEQYlE8oMRe3Rst0PM3/YI34z63U7Y8nhbGz2xUPw9MYj9IggBhYj8YMRupdgFf8/zhpC7nXDlj4ztpigv/yFc/bP2tok/C8Uu7P4E3X49mVQVhCGPiPxgp513rn0+MQTY7TS22cXYPKeh0hgc0LDo3jZhdzcACpZ+xzi/odI41369HU9LVC8IQxwR+aHEonshNLIt2jYjbTRs+bWxzy7uvjZOsxNCImD5I8b3nf8Jk1fAqc0QGtH2tGAOCDKpKghDHpl4HUr4TnJadovqaNVY+zRMXGZsO/KuMRigDbEHSJzTdq5pB235Nez4HbxzJ9RXgLOqbyZ1BUG44EgkP5Tw9cjtkbavZ27f1+yC0zugptDYVrgHrnuiLXo3z212Gp8Tl8GZLCPCP7UZKrIhZ2Nb23YLR3x7YZBS1eDmrcxCbslIIT4ydKC7M2CIyA8lzOg8fzus/qP/XHe76C57CMpzIPtDY9/EZRAUAqc2wYm/tYn0jqeN40PCjePSLoP0W+CjfzUi/Uv+ydjmb1JXsnGEQcpbmYX86qNsAO5dPmWAezNwiMgPJebdbgh8zkZ499ttQm/HPhBc80v4n1uhKhfGTINv/NmwXjY+CilLDXG3e/aL7vU2oiAiHv7JVlMu4nbY8ydj36J72q4rvr1wgTjXyPyWjBScbg9OdwtVDe5hG82LyA81EmZC5UlD6O2pkHYLxxwIACpyIDweJi03BP6jfzUsmBa38bn0uzDtGph+vXF8yQHj3MIvIGWRIfyRo41rbPm1cUxohFg1wgWnq8jcHABWzUrk06Nl1kAQERrMrz7KJiI0aNhG8yLyQ4k9fzIyYsAQ5s4WS13+b0aKpMsBqcugYAfsfQGqCwxhB2P/wn+Go+9B9WkYP9/IzsnZCPGTDUvn1Ka2jJ15t7elXfpaNXb7SBD6iVsyUtp92jEHgN2nKtl03HjD3L3Lp7Q7Z7h69CLyA0l3kbC5f/r1hofudhnb4ycbVgwY2S+TVxjHmKIbMQacFW3HLvkuhIYbPntTgyH6RXuh9DC0NHovZhPvhoq29Mrp17d59lc+2r5/826H3M+NgWHP820LtAThArMwLZ4pCZGsXZLG4smjLXGPjwy1Ivg/bcnt8ZOA7xPBUEZEfiDpbtLSHinnbDRWtU67xvh+4m/GMbu8kf2Jvxmie+DPhkUTOxFGBEHVKYga0xblT7jEEHloE/iwWJiyytifshQ2/9yI8iNGQ9Zbhk3jL1qPHA0pl3qfDmwLtAShj7BH33a7xvxtivDvP88ht7yBV3bl8+Q35vmN2M/lScD3iWAoIyI/kHQ3aWlun369kd1i+ubjLzH2OasgdxOMSzf27fkTTFoBkWNh/AK46CbY9hvvPm8OfEwKJM2H1hYj2lcjoLbIOC5nY9tTQHWhMWG78J+NSVtzDsB3MPJdoCUIPvjaJOdim9jF98c3zgLoIPj3Lp/Cj2+chduTxbTEUbyyM4+nPztp7TOxR/W+fVs1KxGAVbMSWTy5zPr0NyAMNUTkB5KelvuNiDeO2/G0EVVPXmFsb3YavjlA1tttE6MRY4y8+LIvjSg7YQYU7zf21RS25cuDsfo1NMKI4MuPG/587ESYcqUh8pU5xpOBOQcAHW2mcy2CJgwrTEF2uluICA3C6fb4FWF/3JKRYouqj/LkN+YRHxnaISqfkhDFFdMT+NVH2Vw6KZ7Lpo62hNukqsHNKzvzAc3apZOIjwz1O5k7ZXlUu8++ZCDmBUTkBysNlUaaZM5GY8IzNNKIyE3r5tTmNrE/tcnIhFn+Qzj0hiHU8VNglDcKyd1iCL7JiBBobTaierTR7sZHjfMARgTDom9DbGrbfIDp+XdW2KwzJI9+2GMKsdPt4VcfZfPQVdP40XUzu42SDVHOIyU+gonxEWw6Xt7OuvFnx9htlk+PljFleZQlrE53C09/lgNARGhwh4nZC8FA5O6LyA8GfCdYTSHN2WgIeeGetoh99R8NW8btAgWMuxhCIiF5Ifz9YRgZbRw38/q26L3iuPE5Mhqaag2BB2hpafPbT+8wFkvVlRoR/MZHjWuZZL3dlk9vDjamfdQVkkcveLlpXjLGf7TabyTrG+W+lVloRfwAV85I8GvV2M/98Y2zmDOhBJfbQ2W9m999cgKApz/L4aGrpvLQVdOs60NHC6e/I+0LPahAH4i8UioM2AqM9Lb3ttb6p0qpeOB/gDQgH/iG1trR2+sFJL4TrNAmiuZiJbtdEhppfG75dVvdmWcWtWXUjJkG89cCyhDvliZj+5z/Y6xqLd5nPB2UHjS2V+UZWThX/V8YPbntCeLgq8b+T34CE5fC5CvbBqKcjcY8QYItOvdnzfTVG6iEIYtdlCNCg7x568EdIlnTS3e6PXzv6hnWYiZHQzPZZ2qZljgKh9ON0+3hoaumsWpWIn/aksuqWYn8/IOjVgT/vaunt8ukeeiqqdaTgyncVQ1u61x7Fk1/R9r+5gX6m76I5JuAlVrreqVUCLBdKfURsAb4TGv9uFLqEeAR4Id9cL2hSVfetO8Eq3nM9OuNxUsZ/wyR8cYxVv34R9oE/uCrhsCPjDGOqciBP682gqa5/wgRcUYK5aJ7jQnY0zsgOMw4Njgc6kqM75t/bqxyXf3Htr5C+8Fn46Nt6Zu+0blYM4IfzKh11axE3j9YwkNXTbUE2j4Zu/uUUfTO1dwKGIK4dukkfvDmQfbkO9iT7+BQoYMv8hwsm2L8P/T0ZzmWRRMfGcLCtHhyy+vZeqKcf1oykbiIUG6aN55Pj5ZZ/alqcPODNw/6zaIZiEi7v+m1yGutNVDv/Rni/dPAzcAK7/ZXgM0MZ5HvSgDt0a49Mt74qOG91xYbwl1ywBBY3/rx5kKlw28ZKZPB4UbGDEDxbngw0/jeUAmFu43vnkZjgvabbxjXKdoD0Smw8d+hwHtMfQVc9i/GNVvc0NLcFsH7E3GxZoYl3VkcZvT6py25XttkGj95N4sduZVsPVHO7791CW9lFvJFnvdBX8PvPjmBy93C9pPlHC2tY+yoUM7WuWluMVJ1d+RWUlDl5JLUGFLiIkgbHUF+pZPff55j7T92ppYX/mmhda3Kejejo0Jxuj1sOl7OlTMS+M7KaQDWJK05qetv9exQpU88eaVUELAPmAo8o7X+QimVqLUuBdBalyqlxnZy7j3APQCpqal90Z3ByfkIoBkxX/5vbSmOaZcZ2TCf/KR9/fiSA4bAj5kG4+ZD1psQFgM3PdfW3p4/GYPGhEVGobLkBYY9M+VKQ+RrC+Hgn9uOL9pj5NiDcd7yR2D6V9pn2fjWs5EIftjRncXhm6bodHvYkWtE7TtyK619Hx8ppaS6kfcOFlNW19SujVpXMw9dNRWXu4X9BdWEh4yg0OGi0OFif0ENl6TGkBw7mu+snMa7B4oIDxlBVUMzD799iNzyBgCOltaw/WSlNfFrt3kWTy6zsmlM22jriXKrn0M5V75PRF5r3QLMU0rFAn9VSqWfw7nPA88DZGRkBO6KmvMRwIRpcNtbxndrwtUJ6V83ttnLC+RsNHz71X+EHf9hbL9kLaQutDWojI+gEAgKNVa1mnXllz9itPva14wsm5ExcPE3Ol7PbjXZ69lAx9LFwrCgK4sjt7yee/470xLae5dPoarBjcvdwqGiGuZOiLEi5/0FNe3OHRmkaPJG7o0ebZ2TFB1GaW0jwSMUnlZj//6CGtJGR/D3I2f4790FAExJiOTfb5jFf207xZSEKMJCglgwMY61S9NwON3ctX4v+ZVOLp0Uh9PtsRUxM/4/mZ0cwxXTE4a8ddOn2TVa62ql1GbgWqBMKZXkjeKTgLN9ea1hiVk8zCwQZmJmu1zzS0NgzZLBIeHtM3fQbW+CAmNQMN8qdfXPvIPKO22++4m/eZ8YzLdG+UTt9no2aPHjhyHdWTU//+AoueUNll9u+vDhocF8kVdFSNAIHE43xQ4nybFhxIQFMzIkiNDgIL7Iq7LauXJGAkdL69pt87RqYiOCSRwVxvGyevIrnRwscLAoLQ5Pq2ZqQhSP/vUwJTWNeFpb+SLPwZUzEtoJPEBIUBBPf3bSmgxeuzSNiNCgIW/TmPRFdk0C0OwV+HBgFfBr4H1gLfC49/O9zlsRusUerftaPr7ZLvZVqL6ZO5OvNN7rqlVbPRv7ilX700OEzWLyV4UycnRbPRvzTVH298UKAU9XpQYAfnzjLAqqjEj+qU9OsP1kBcUOF1tzjMnO7Scr2mXGTJ+RwJPfmMcrO/NobmmxrJg5E2JZPj0BgMljImhsbmXj0TKqnc00el92E6RgZlI0/73LWO+xv6Da6mdzi+bKGQnWoipT4C9JjWVyQiRa63a+/FC2Z3zpi0g+CXjF68uPAN7UWn+glNoFvKmUugsoAG7pg2sNX3r6Fihobw2Z21KWQmWud+HUpe19fXu5YvC/mtVfFUo7kaONtszIX6L5YYHdqvG7ejQhire+vZS3MguprG9i+8kKtuaUk1/pJG10BDfPS2b59AQamjwopfjxjbOs/PjLphrzQVUNzTz9WQ4RoUG8evelAPzyb8eodjWTEhdOocNFeMgIXM2thAUHWZOw42LCOFNj1GfSWjMtcRRzJsRw07xkpo0t4FBRDaCtQcFcPBVo9EV2zZfAfD/bK4Gretu+4KUrT78n+3Y8bSxyAqzMHOi4ghX82y72qN2OPTVUsmuGHfaotzNv3jymqsHN6KiRLEyL5/ef5/DjG2cxJSGKP23JZU++gx9dN5MpCVEwC3afqmTtkjQKq5zkVzq5bOqYdu0eKTb8+zFRoUwdG8V3Vk5jy4mz7D5V5fXZ49Fac6amkZS4cA4U1nCgsIYrZySwdukkyy4CuHRSPIsnxw95770zhs2KV4/DQc2GDcSsWUNwXNxAd6f3nGtNmHm3t6Vemi8C8Y32/b3erzt8U0Mlgh82+PrxXRUAM48x95uVIpkFTneLlTv/u09OsPtUBV/kOWhuaeV0lWGrON2edu394JoZhH6ew7TEUTy/9RQAcybEWsIdEqTYfrKKK2ckWMfERYRYpRHsVVMXTx7N966e3uP7HGoMG5Gv2bCBs0/8FoDRd901wL3pJfa6NtBeWH3F3/67s3rvvk8C5yLUEr0PW3qyOtS3OJnvylJzMdKPrpvJp0fLrNoyYPjoybFhFFc3sr+g2ltcDCvX/uU7FlHV4CanrI5Nx8uprG8kJjyYr8wex73Lp1g57oB1zGVTR+N0t7B8+lgy8x3MTo5h7dK0Xt/nYGbYiHzMmjXtPocc9iyZjY92Pgnr+7am7lahlue0ZdMkTDv3fklu/LClJ6tDfYuTgSGUq2YlsvVEOSlxEdxz+WScbg83zUvG6W7B0dDE1pwKvsirYnxMmK21tug7M7/KSnl88hvzrBWsAPtOO5iSEEVcRqgVgZtPDkYFzBy+LKpmR24lGWnx3S58GuqrYIeNyAfHxQ3NCN4Ud3dDWzExe068v0lY8xj7ZGpnZYLNAQPasmoEoQd0lYVitzhMPz4iNNh6Dd/PPzjKjtxKduRWWlkvEaHBrF2axg/ePEh+pZMpCZFWfv2itDjWLp2Ew+nmgy9L2ZFbyS1/3Mm/3zCLV3bl852V00iJi2BrTjlPfH0u4P8pAowKlGa9eHPw6eolIUM922bYiPyQxbdWjW99G18iR7evPeMbafu+EzZhhlGuwFxdKwh9gK/F4fsavk3Hy0kbHcE1s8fxfxamWC/qMCPyK2ck8OMbZ/HT946w/WQFocFBOJxufvJuFrnlDUSHBZNb3sC3X82k0aNxe1q5YnoCG+5fRnxkKLnl9WzLqeCeyycBul1f7HXjzcEnkF4S4ouIfD/Qp5O8/lInE7qxR7qyUHyzanb+vm0hlCD0EV1ZHPa676MjQ5mSEMWU5VGW+F/pzZWPjwzlP2+dz3dfP8D2kxXtFjBNHRvF/oJqGj2a+MgQkuPCraj9e1dP5+cfHGX7yQpCghRPfmMeYEzwtq1qNbBn/gQqIvL9QF9M8rYbKPrS8+4uq8YXebOTcB50Z3HMmRDLHG9JAxP7wGAXYqMGIt5UytEsmBjPTfPG8/7BEjLzq9iRW0mxw2UeDeB9VeBR65WBXxbVeC2hIL/9GuqTq10hIn+O9CRK74tJ3n7PBuqpeEv5YKGPMRY75fCj62b6jap9j92RW8mitDhCg4P4vzfPNnLpMerG24uf2bNppiRE8fIdiwDaPSF0ZscM9cnVrhCRp2vh9t3nK77+zu2LSd5+zwbqqXhLiqTQx5yLoLZl5xiv7vNdlWofGDpbrdrZE4KdoT652iVa60Hzt2DBAj0QVLzwgj46Y6aueOEFrbXWzVVVuuKFF6zPzvbZzz19zz3WtiFBfYXW2//D+BSEQU5lfZP+4+aTurK+aaC7MigBMnUnuiqRPB2jZnu0bm6LWrmSyhdfJGbNmnZResyaNTTs2UPDlq3UbNjQbt+gXmUr+e3CEKKzSHuor0a9EIwY6A5cSDwOB5UvvojH0fWrZqNWriRi6RI8VcYS6dF33UX9559z9onfUrNhQ7v2ajZsIPGRRxjzwAO0Op3t2na8+ipnn/gtjldf7Z8bEoRhjjlhapQqEPwxrCL5ziYz7dujVq6k8N5v01xQgHPnLoLj4xl9110don2Pw0HJI4/QsGUrACMiwjn7xG8ZERFha1tZn4M6qheEQcD5ROWBPGHaVwwLkTcFNmrlSqDjZKbdkil84AGaC4w3y4RnZFj7guPiiFmzxhLqmg0baNiylcjlV7Rrz/497vbbGBERbh3f3YStIAxnzieNMaAnTPuIYSHyjldfo+KZZ6jfvoPkp57sIKpmNkzliy/SfCoPFRaGbmwk8tJLCY6LswS51emk4plncbzzDuN/+UvGPHA/bdG6QUt1tTWg1H/+OVErV1q/W51OWp0uq73OUiRlABCGIxKV9w/DQuRbXd5ypbt2dZgctROzZg3123fg3LWLiCVLCJs7h5PXX0/kkiVUv/YXxjzwACGTJ9F8Ko+K554jctEizj7xW1xZhwlPv5iKZ56xJmF9Px3vvEPk4sVU/+V1QBN3++3WNX0JqIqZgtBDJCrvHwJW5O3R8Ihw42XVEUuWdBBV36g54pL5hshfMp+yxx+n+VQetVUOYm/7Fg17viBi/nzc8aMJSZ6Ap7KK8IwMGrZsJSQ5mcjlVzDmvvsYOXUqri+/JP6uOxl19dW4Dh+m+VQeDZ4WAJz7DxB3++1dDjb2T0EQhPMlYEXeHg3bvXFf+8M8rtVpLIt2nzlDyMRUIi+/nLC5cyn+7kO0VldT++HfaK2uxrU3k/CFC6n+y18AY+AAaMo5iWvvXrS7GQDX3kxUSAiNR4/RWuUgJDWV8b/5NRXPPWelW9o9fnu/hmzFTEEQBh0BKfIeh4NWp5MxDzxgCagpmr6Ruxktm367SfnT/4kKDUE3Gu+IbK2uNsT/ssuo/fQzANSoUcTfsQ4VGkJIcjKuvXtx7tpF/J13okJDGDl1KlUvvgRA5OWXETFvHuMff7ydwJ994rc07NnD+McfF/9dEIQ+p9d58kqpFKXUJqXUMaXUEaXUQ97t8UqpT5RSOd7PC6ZgNRs2UPHMs4yICO80cq/ZsKFD1k3sbd8i9KKZjIgeRUtdLQ1btqLCjJcWRCxZQtobbxA6fjytZWUA6Lo6HK+9RsOWrYwIDyd8YQZh8+aBgvGPP87ou+8mZGIqAPXbt3fIz49Zs4bI5VfQsGUrxT/4AeW//0O3OfyCIAjnQl9E8h7gB1rr/UqpUcA+pdQnwDrgM63140qpR4BHgB/2wfW6pStPO2rlShr27CFq5cp2WTfOXbuIXH4F7mNGCldT1hFGxMcZVsvkSVZWTsyaNbiLS6jfupWo5VcQf/vtVE+ZQu1HH+EpKQWg8eBBGo8dI/nJJxn/619TeP/9eE4XGANLZRVVL72Ep7KK0f98N+HpF6Obm3Hu3IVz5y5GRISLVSMIQp/Ra5HXWpcCpd7vdUqpY0AycDOwwnvYK8BmLpDId+Vp13/+OQ1btlK/aJGVdRM6KY2oy5YRtXIlVckTqN+6lYhFi2htasK5bRvh8+e3a6Mp5wSeoiIajx6l9oMPcB44YAn8iOhoWmtrce7cReV/vUBT7klaqxxELr+CqJUrKfjnewBwffml94njGcY88AAR8+cDSiZbBUHoU/rUk1dKpQHzgS+ARO8AgNa6VCk1ti+v1Rnd5Zjbo3zHq68BEBwXby1Qai4uwlNUhCclBeeuXQDUvrOBEWFhNOzaReTiJbj2ZgLQeOAgjQcOokaNAkBFR5Py/J8of/o/ce7aRd2nn9JcUEDk8itIfOQRyh5/HE9RkdERpdotzvLXV8mXFwSht/SZyCulooB3gH/RWtcqpbo7xTzvHuAegNTU1PO+vu+CpVankxEREUStXEntBx/S6nIyIjyCuNtvs6J8M+smauVKyn//B5wH9uPcadg2I6dMxblrF0GJiajgYGo//YzWsjLqXI2EL8zAXVRMS6kRveu6OkbERDNq1Soatm0j4aHvUlJSTPPpAsIXLiQ8PZ3aDz6kYctWKxvHuWsX9Z9/3mmGDUi+vCAIvUdprbs/qrtGlAoBPgA+1lo/5d12HFjhjeKTgM1a6xldtZORkaEzMzPPqw+VL77I2Sd+y5gHHmBERDitThcVzzxjLV4yGfPA/YyIiLBE1bcGTeTyKxj/+OOAUWCs5n8/oLmggOCkJDylpYwYNYrWujri77qT6jfforWurkNfIpYuwbnTeAoITk7GU1xM/J13Ejw6npg1a2iprubMz/4fYRfNBKWoevElxjzwAAnfebBd9A60GwAkshcEwR9KqX1a6wx/+3odySsjZH8ROGYKvJf3gbXA497P93p7ra6w2zCmILqyDlvRc9hFMxkRHkFro8sb6btI+M6DVg2aEbGxRF21ktBx46w2XVlZNBcUEDJ5EuN/+UtKHn2U5lN5hEyeRKvLZQn8iLFjGTkx1bJxQidNAg26uRmXd9AaER5mReM1Gzbg3LXLWHS1dIn3atraZ4/eOyukJpG9IAg9oS/smmXAPwKHlVIHvdsexRD3N5VSdwEFwC19cK1O8Z1sDY6La5eTbka+5b//vfcIQ1Rj1qzB8c47NJ/Kw3XgALWn8hgREUGr02kMEEuXkPDd71Lx3HOMeeBBKp5+mvD586nfus26Vuj48ST97GdU/fnPNGzfAdqwY8Y88ACRl16KvYyBeU1PZRWuw18SnJRExJIlRN94Y4f8fl9kJawgCOdKX2TXbMe3SlcbV/W2/d7gL8sm7vbbLbvGZNSKFbjiRzNyxnSCr7uu3aRsxPxLrFWq7sJCmgsKrCqVI2Jjaa2upvHgQcNz372b5oICmk6cIHL5FUTfeAMjJ00CjDkDY4BRxN1+G8Gj43HtzbSi//rPPweg4plnGfvwv3b7GkJBEISeEJArXrvCFP6mvDyKf/ADw1LZm0nk8iusImQ1GzYQfeMNVimE6OpqzribCZ2Uxogrr6S1sRF3Xh4J3/0uDdu2A9rIzDmVR/CECYTPnUPViy9Rv2gRI20WjX1Fbdztt9HqdFkTwqZXb+bw+yJWjSAI50PAi3xnEXDZ449bk6NmimNt+sVWhg20ialZuKzimWeJv/MO3KfycO7aRcP8S0j4zoMAnL7jTgBCU1IYfffdBMfHd6gzX79jh7dtwyoaERFuRPXefpnzA/bBwX6+/VMQBKEnBLzIO1591UqpTPjOd6ztiY88wpnmZsIuuojRd98NgCvrsJVCaRdTj8OBc/8BABqPZVv586ZYA4z7yb9T9vjjJD7yiN9+BMfFkfzkk95XASprta29bk1XQi5FywRBOB8CXuTtr+Az8Tgc1H/+OclPPgngza93WW96MlMozRd3m9kwZsRf/dbbNB47RvSNN1ptjpw0yZroNdM3oaO14srKomHLVsY8cL9Vt8ascS9CLghCXxPwIm8vM2xi97cBb379/Yx9+F8tW8fMu4eO6ZnBo+OtxUwj/aQ42tuyY39lYNzttxMH7XLiBUEQ+pqAF3l/0XHMmjXWq/iib7zB2mb37H2F3d5GZ7aK7zm++Ntvb1cyaARB6GsCXuT9ERwXx4iICM4+8dtOqz52ZZ10tq87u6W7/ZJBIwhCXzMsRR4GZ7bKYOyTIAhDmz6pXdNX9KZ2jdAzHI0O3j35LqunriYuTCwhQQgEuqpd0+s3Qw01PA4HlS++OKjewHQh+/TuyXd5at9TvHvy3X6/liAIA8+ws2suhO99rhOoF9KLXz11dbtPQRACm2EXycesWeM3vbEvsb9HdrD0ySQuLI470u8Qq6aPWb9+/TltF4QLhXjy/UB/pUKKnz44CQ4OpqWlhbVr17YT9XXr1vHKK68QFBSEx+MZuA4KvcZZW8ORzZ8ye8UqIqJjBro7HRBP/gJjpkr2da67+OmDD1PgAV555RXWrVsHtAk8QEtLC8HBw84ZDSiObP6Ura+9zJHNnw50V84Z+S9vCCF++uDCLvAmr7zyCps2baLAW47axBR6ieiHJrNXrGr3OZSQSH4IIX764GH9+vUdBH41q4klloKCAmKJZTWr2+1vaWkRj36I4aytYe/77wCw8KavdWnVmMc6a2suVPd6hIj8MMbR6ODlrJdxNPpP3exu/3Bm3bp1rF271vq9mtU8xEP8jt+RRhq/43c8xEPthH7t2rWWnSP0jr4U1K7aOhebZrBaOmLXDGNMjx/gjvQ7znm/yXCdEDaj8ldeeYXNbOZmbiaNNF7mZQDyyWczmwE6TMoKvcMUVDAi7P5q61xsmkFr6WitB83fggULtHDhOFV9St/3yX36VPUpv/urXFX6pcMv6SpXVZftvHT4JZ2+Pl2/dPil/ujmoCc1NVUDOo00vYlN1l8aaRrQqampA93FgKOhplrvee9t3VBT3evze9vWYADI1J3oqtg1wxRHo4NffvFLthVv46O8j/xaM/7mAPwdt3rqar6/4PvtJoSHi9Wzbt06y4P/KT9tt++n/NTy6MWm6VsiomO69ci7wm6tdNXWYPXZz4U+EXml1EtKqbNKqSzbtnil1CdKqRzv5/B5jr/AnIugmse+nv06u0t3W9t7mp7p7zh/g8FwSPe0p0muYAVppJFPPndwB/nkk0YaK1gBtE+vFPqHcxHk2StWccVtd3RprThra/j7s78blD77OdFZiH8uf8AVwCVAlm3bb4BHvN8fAX7dXTti13RNZ/aJP7uks2N/u/e3On19uv75rp/rZw48o3++6+f67o/v1gfLDlrH+55r/91TC6cnx/W0rcHIyy+/rDHe/2j9rWa1jiVWAzqWWL2a1R2Oefnllwe66wGF3WrZ897b+rffuEHvee/tHp/TFWZ77/zqp31i5fSnLUQXdk2fTLxqrbcqpdJ8Nt8M3jAGXgE2Az/si+sNVzqbCPWXP/969us8d+g5dpbs5DdX/MY6P6vCeNjKr83nscWPcf+n97O7dDdlDWU8vfJp3j35Li6Pi+cOPYfL4yI8OJyqxirWH1lPVWMV8WHx1nVeznq508lWM7o/n/sZCqxbt4677767XRrlu7wLQGpqKgUFBdZvk6CgIInmz5POVpzaJ03NqHzygkvZ+/471m/7eWZ0nncgk8Kjh7n2/u+1O8b+ffaKVTQ3NqIVfUJfThafC/2ZXZOotS4F0FqXKqXG+jtIKXUPcA8Y/3MInWMXc9+MljvS77CsmNVTV+PyuADYXbqb2/92O0vHL+WN42/wzRnfBKC5tZm8mjweXvgwRXVF5NXm8dDnD5FXm8d9c+/jvrn3sffMXjLLMslINFZLZ1dlt7N4eivQQ31xl8fj6bAgysyisVs5gJQ26CWdCaRd2E1xth8LtBsEPvzPJyg4fJDYcePJO5DJ35/9HeMmT2PXO693OH7hTV8jJCyMra+9TFluDtfe/z1roDifEgeTF1xK4dHDTF5w6fn/gzgPBjyFUmv9PPA8GLVrBrg7g4quUhPNKHjvmb384rJfEBcWZ21zeVxkV2VbxxbUFUCJ8T0uLI7w4HC2FW/jib1P8IvLfsG1k67lQNkBdp/ZzeXJl3PdpOt4Yu8TZJYZdYSaW5tZnLSYtbPWEjIihBUpKwDYe2av9f186Em0P9ixC709TdKeXikCf/6YgmoKo6+Hbk6a7n3/Hb8RfdbnG0m9eK41CBQcPgjA1IzFVBYXkHcgk5ZmD3OvuYH8QwdY+o3bWfz1W2lubMRZW8PsFasoPHqYvAOZHNn8KQtv+to5R+TO2hoOfPwBJcePUnD4ECmzLmZ08oS++kfULf0p8mVKqSRvFJ8EnO3HawUkvnaG/ffqqavZe2Yv24q38e7Jd7kj/Q5WpKxg75m9uDwudpfuZnHSYtKi08ivzefBeQ+y/+x+Vk9dTXVTNU6Pk1p3Ld/b/D32le1jcdJi7pt7nyXw24q3kToqlYSIBPaV7QMgZEQI24q3sXDcQgDr+6SYSQPzD2iQ4PF4rOjdzvr161mxYoVYNL3AV1DNyVV7FO2sraG5sZHFX7+1Q0Sf+YFRCfb4ji3Mu/ZG3E2NKA3zrr0RV10d1WVnKMg6SF1VBY6SIuqqKpi55HJ2vfM6IWFhLLzpa1x7//esgWbv++90OuB0dQ+73zaeFCbNz7jgefT9KfLvA2uBx72f7/XjtQISXzvD/hkXFscvLvuFFekDbC7czLbibaSPSW+X0vjuyXdJjU4lNTqV17PbHksPlR8CYFzEOHaX7mb+2Pl8lPeRIfBRqRTUFXDD5BtYNG4RANdNuo6F4xa2s1eGqtXS13Qm5CLwvcN3gZG/KPrI5k/Z9c7rXHHbHZza90W7iD7/ywMUHD7Ise1bmLFsOctuuc1q+8DHH+AoKSL14nksveV2Pv7jf+AoKUIr2mXemE8LO956jd1vv467qZH5X7mxU8vG186ZvWJVu8HlQlex7BORV0q9jjHJOkYpVQT8FEPc31RK3QUUALf0xbUCja4sGV87o7Pfphc/LXYak6InsWz8MlKjU9tNou49s5cxYWP4a+5fASyfHUApY2bp/dz3uTz5cgBWTlxpTbLa+2WP2mUVrNDfmAJr4m9y1d9K08kLLuXAxx8wduJkasvPUn2mhC1/fsGKymevWIXymsNjJ07m9KF9TF6wCHXJIpQ2zj/49w9obmoiZORI5l17o3V8yfGjNDc1se9/N3Bq/15SZl3cTrx9B6KI6Jh2g8uFpq+ya27tZNdVfdF+INMXGSZmG3Ej43A0OXjwsweZEjeFfWX7uG/ufSxIXMC24m2MChkFQOzIWH6y5CdsyNnA4YrDeFo9lDaUUlxfTH5tPvfNvQ+gg0Cbwr0iZQWbCzd3KeBmdo/L4+L+efef130Jgi9defC++02LZO41N6BGjOCiy1byxk9/iKOkCGdNDaW5J0iaNpPS3BMUHzOyzibNzyDvQCZHt2+mpqzUajckLIwZy5ZzbPsWCg4formxCYCiY1kUHcuyrB17fwZLeYMBn3gd7vhaMb4RcE+E1e7PhweFU+2uZl/ZPi5PvpxbZ97KgbMHAIgIiSA6NJrLJ1xO7MhYfpDxA549+CzPHXqOyOBIokdGW979U/ueIjw43Bp4HI0OHtv+GNuKt1nXgqGX+igEBnYhPfj3D9j1zus0Nzay9Bu3WdtNi6S5qQlHSRE73/oz1WdKiRs/geLjRynNaUtOSL5oNklTZqAVuF2NFGdb6zpJTZ/H5AWXsuXPL1BdZmQwjPC+H2DCRemkzLrY70AzWBCRH2B8LRjfyVVfYbVn09jbMP35S8Zewh8O/oGZ8TO5M/1OANKi08iuyqbMWcaEqAm8cfwNTlaf5KkVT1ltNHgaaPA08IeDf+DRSx/F5XHh8rjYXrSd3+z9DfMS5lmTsffOuZf0Mem4PC4cjQ6rL/YB6taZtxIeHC6evXDe+EtVtG8zhdTMYy84ehj91mvM/8qN7SySnW++BkDyzHTUiCC+8u1/Yfvr69tdK2nKDCvbJjZxPNAm4DOWLWfLn18g70AmqRfPZfyMWcxcupxT+75ol3/vOyE8WBCRH2TYI/t3T77LtuJtXJ58OQ8vfBgwMlrWfrSWp1c+DcATe5/g3jn3WpkzcWFx/Nc1/2W19+zBZ3nj+BsApI4yJlMBMssyeffku1w36Tp2luzkVPUp6prr2F26m82FmwF47tBzxI6MpbqpGkeTUTKhoK6A/Wf3Ex4c3iHa903rlChf6A2dTbLat1UWF1Fy/ChJ02ZSfCyL4mNZVk67efyMZcsBKD5+DEdJEflf7iNx6nTO5p/C7XIyMiISgLwDmQSHjqS6rIRx02YyIiiYGcuWk71ziyHw6fO44bsPWyJuT4M0+2UusBpMQi8iP8iwR/b+smnWfrTWWrg0Onw0mWWZnK49TUFdAXvP7OXeOffypy//ZA0KH+Z+CMDs+NkU1hcCEB0azVUpV1HVWMXPdv3MyrJZMHYBFydcTFVjFXvP7AVgzug5HK48zKOLHuVU7al2/TKjfbNmjsvjYnHSYrYVb+Ox7Y91eOIQhHPBn7ftu23Ln1+g4PAhYsclARCb2LbIKTohkUMbP8RZU0NETAwFWcYiqPxD+zmTc5yoMQm4XU6anA2U5p4gdlwS1WcMH76quBC3s4G3f/HvePMSSJg0GcBKo7RH8vZ8+gMff0DoyLAOxwwUIvKDGH+LhZanLKeloIW82jziw+MBQ5zBiPIPlR+i1l1Lc0szLo+LgvoCIoMjqXHXUOuuJTwonFp3LRWNFVamjcnFCRcTHhzOc4ees7bl1+XjaHJQ6iy1JlBNWwaMaD88ONz6ft/c+6x8ejN/XxD6Cl+/e/k/3g3Apf/wTUqOH2H8jNn8/dknyTuQyaiERACKso8SFBxM9NhEqs+UwBnj3Inp88jd9wWNdbUUH8tiwVfXkLt3N9VnSnA7GwCoryxvd30zYt/7vxtweQuhmRk0ZuZOc2OjFdXnHci0jhkoROSHEO+efJf1R9Zb2S8uj8vKYS/INWyYWnctAHm1edZ5DZ4GxgWNI3ZkLI8teoxTtadweVwkRyVzrOoYZQ1lnHGeIbsym5mjZwIwLnIcE6Im8C+X/ItlBdn78dS+p7hv7n18f8H3WZGygo/yPmLd7HUAPLzw4Q759IJwrvRkZWn4qFGkzLqYsKgoAE4f2mdF4xGjRhE9Zgx1lRXUni0DjEnUuPHJOEqKAWisqyVqTAJTLllEaOhIVv/bT/j8pT9SkHWQCRelE5OYxMm9u2hqqKfs5AkSp04nNnE81WUlxI2fwOwVq6gsLmLLn19g+T/ebS3YCgkLY/yM2QDW4qnzLYfQW0TkhxC+fr0ZOTsaHSRGGFFLq26l3FVOmbOM2fGzqWmqobGlkdyaXACOVB0hPiye5w49x+XJl1tWTUxoDLvP7AYFi8ctZveZ3Xxr5reYO3Yuc8fOtfrgaHTg8ri4b+593DrzVuLC4ng562WeO/Sc5fm7PC7iw+K7vBfJoxfs+BNAM0PGLDHgTxjtXnjegUwWf/1WUtPnUZB1kLJTJ0lNn0ft2TLCIqOYsWw5S2/5Fkc2f8qhjR9aFk99RTnHd26lsb4OZ10ttRVlzFi6nHpHJdVlJTQ11BM7LslKl5x7zQ1wCFLTjf8vPn/5jxQcPoijtIRb/98T7dI48w5kWmUMArFAmXCedCaA/vx6l8dlTawCJEUmWd8XJi3kV1f8il9+8Uvya/I54zzD4YrDLBq3iPvm3sey8cssP7/GXUNyVDK7S3cbRcwUVDVWtcueAazB5fsLvm9tt6dwQvtCZp3ZNUO5AqXQ9/gTwIjoGEJHGgXCzpzK8TuhaXrzZtScNmcBSmNE66VFxCUlU5B1kMaGesKjowFwNzUy95obOJuXi8ftpr6qksb6OmITx5O9YwtNDfXUV1biaXYDRpaN1lB9ppTUi+cSMSqa6rISqjeW4CgpJm58MgWHD1J9psSqb+OsraGhtsZKv7T3NZDKGgjnSU8E0L7a1eVx8fe8v3PGeYblE5Zb4mtG2r+54je8lPUSWRVZuJqNFbDrZq9j/9n9FNQVWBH4uMhxFNcXs7N4JwX1Bewu3U14cLiVChkXFmfVxzELk5kD0sMLHyZ9TDpglD8wc/o7Y6hXoBT6ls4E0F+BMDu+UTMYWTKT5mdQcPgQ42fMYty0GZzJOY6rttaqI2NOskaNSQAganQCkfHxVh58dOI4qooMC3REcBAFhw+Rmj6P5OkXMWPZcopPGMXGCrIOMn7mRSz52q1o1b78wr7/NermZO/cQujIsHZpnxcSEflBSHcC6Gh08Hr261a99zvT7+TO9Dvb1bF5Pft1Xs9+vV3BscuTL7cqS2ZXZfObK35jZciEB4e3OzYjMQONZs+ZPVaBsjvS77Dq45iFyewDkn1la3dFywKhAqXQd3S2gMg+oTl7xSqctTUc/PsHaEW7XHV7uYOUWRczfsZsWjzNuGprcVZXA1BZeJqIUdEs/vqtuGprqT7zIaPiRlNfUU68N+KPHTeeqRmLaXY3UVVUQPLMdFbecR+n9n2Bu6mRXW+/zplTOSy95R8BxdiJk628fDv2WvRKMyA2jYmI/CDEVwB97RvTMjExc9VNr96sVwNw4OwBqyLlvXPuxeVxodE8OO9BK0Nm/ZH1Vonh9DHp1kKn9UfWA3B58uVdFkuzfwpCX2OmKJqZK2bt97LcnHbZK6aAjk6ewN733zEi7cPGnFN4dAwxiUnseud1lnztVpbe8i1iEsZaaY7mita8A5lExMTgbmoEIGnqdGsgcdXVkb19s1Ge2OOh4PBB0ubM73QSNSQszBp87N8vNCLyQwBf+8Z8KYij0UF+Tb5lnZj1YtbNXtcuA8esMLn/7H4yyzL5/oLvW6UL7pt7H5cnX2756duKt/H9Bd+3rr04aXGn+e4yeSpcKEzPfvHXb7WskZlLl3coKQBtpYcXfHUNSkPpyRMUZ2dZtWi0av/kMDp5As7aGsZNnkbilGnt2jNrwJuYmTtjJ04mbc78ToXbd45BUiiFLjGj5BUpK6w3P90/735eznqZN46/webCze3skfDg8HY57fZX9pntfJT3kZUhA1j1ccwo/rpJ17Xz4k3sAw70/u1QgtAT7J69PXL29/INe+nhhTd9jZ1vvkZxdhaeZg/JF6XT3NRkZetYL/TIPkZB1kGuuO0Oq/3QkWEUHD7EpPkZTF5wKdk7t5Bx4xqCw0a2s2jsmUHm9c+15nx/IiI/BDDtm5ezXu4Q0UPbIOCvXoyv9WO245shYz/GtHr8VY+0DzjmQCFWjdDfnEvRL99J3BnLlpO9axtnThoFyYqPZVFVVGB5/Z290MPu85tWjjlw2GvVdPa6wcFSpExEfgjhK+rd1ZvvaTv+OHD2QIf0Sbs94y+VUhAGA74Dwql9X+AoKbIWMZnvdzWjb3dTI57GJoLDRgIdi6CZmTv2QcDfC8T9DRCDARH5QUhP8uR7225n7dw681ayKrL8liXwrZBp/xSEwYo9IjcnWbN3bqG50ZhcXXbLbVZ9+tCRYUD7aNyfVeS7zT6oDJYI3kJrPWj+FixYoAWtXzr8kk5fn65fOvzSgLRb5arSLx1+SVe5qnq0/VzbEYYHDTXVes97b+uGmuqB7kqHvux5723922/coN/51U91Q011u/2d9Xsw3Y8vQKbuRFdHDPQgI3Rk9dTV7d7Rar7ez6z22BVdHevbbmfHm08Mnb2OsKf2jBn5m6mawvDCtDSObP50oLvSoS+zV6yy3gJ1ZPOnVjQeER3T7ntXbdgxPXqnt2jZYELsmkFIVy8S6c6u6epYf3aPbw34vvTXxdIZ3lzoZfz2hVK+C5R8++K7yKondHU/A1WXpicoI9LvxwsodS3wNBAEvKC1fryzYzMyMnRmZma/9mcoci756Oeau25/rd/3F3xfUiGFIYv9va9mFsyFYqAqTJoopfZprTP87utPkVdKBQEngKuBImAvcKvW+qi/40XkBwZZ1CQEAl1F8oFOVyLf33bNIuCk1vqUtyNvADcDfkVeGBikjowQCEREx1gv8hba6O+J12Sg0Pa7yLvNQil1j1IqUymVWV7e/i0sgiAIQu/ob5FXfra184e01s9rrTO01hkJCQn93B1BEIThRX+LfBGQYvs9ASjp52sKgiAIXvpb5PcC05RSk5RSocA3gff7+ZqCIAiCl36deNVae5RSDwIfY6RQvqS1PtKf1xQEQRDa6PfFUFrrvwF/6+/rCIIgCB2RsgaCIAgBjIi8IAhCACMiLwiCEMCIyAuCIAQwIvKCIAgBjIi8IAhCACMiLwiCEMCIyAuCIAQwIvKCIAgBjIi8IAhCACMiLwiCEMCIyAuCIAQwIvKCIAgBjIi8IAhCACMiLwiCEMCIyAuCIAQwIvKCIAgBjIi8IAhCACMiLwiCEMD0SuSVUrcopY4opVqVUhk++36klDqplDqulPpK77opCIIgnA+9fZF3FrAG+JN9o1JqFvBNYDYwHvhUKTVda93Sy+sJgiAI50CvInmt9TGt9XE/u24G3tBaN2mt84CTwKLeXEsQBEE4d/rLk08GCm2/i7zbOqCUukcplamUyiwvL++n7giCIAxPurVrlFKfAuP87HpMa/1eZ6f52ab9Hai1fh54HiAjI8PvMYIgCML50a3Ia61XnUe7RUCK7fcEoOQ82hEEQRB6QX/ZNe8D31RKjVRKTQKmAXv66VqCIAhCJ/Q2hfIflFJFwBLgQ6XUxwBa6yPAm8BR4O/AA5JZIwiCcOHpVQql1vqvwF872fcL4Be9aV8QBEHoHbLiVRAEIYARkRcEQQhgROQFQRACGBF5QRCEAEZEXhAEIYARkRcEQQhgROQFQRACGBF5QRCEAEZEXhAEIYARkRcEQQhgROQFQRACGBF5QRCEAEZEXhAEIYARkRcEQQhgROQFQRACGBF5QRCEAEZEXhAEIYARkRcEQbDhqnezf+NpXPXuge5KnyAiLwiCYOPYzlJ2bcjl2M7Sge5Kn9Crd7wqpZ4Avgq4gVzgDq11tXffj4C7gBbgu1rrj3vXVUEQhP7noqVJ7T6HOr2N5D8B0rXWc4ATwI8AlFKzgG8Cs4FrgWeVUkG9vJYgCEK/Ex4VyiXXTCQ8KrTbY4eCtdMrkddab9Rae7w/dwMTvN9vBt7QWjdprfOAk8Ci3lxLEARhsDEUrJ1e2TU+3An8j/d7MobomxR5twmCIAQMQ8Ha6VbklVKfAuP87HpMa/2e95jHAA/wmnman+N1J+3fA9wDkJqa2oMuC4IgDA5Ma2cw061do7VepbVO9/NnCvxa4EbgNq21KeRFQIqtmQlASSftP6+1ztBaZyQkJPTubgRBEOiZV97VMefjtff2mv1Fb7NrrgV+CCzXWjttu94H/qKUegoYD0wD9vTmWoIgCD3F9Mo9TS0EjwzioqVJHSZSzWOADtF4T87v7Jr+2juXY/qa3nryfwBGAp8opQB2a62/rbU+opR6EziKYeM8oLVu6eW1BEEQACMiPraztJ342reZHnlzUwu7NuRSfMLBqnWz2h3b3NTCwhvS/PrpvudD56JsXnfS3DHtzvXHQHj4vRJ5rfXULvb9AvhFb9oXBEHwxVXv5tP1RynIqgLaxNc3Sr7kmom46t2cya2hIKuK/R+fJnxUKJPmjmH7WzkUZFWxZM0UvxF6Y0MzxSccXLxiAqnp8ZaA+8N+3YuWJnUYfOwMhIffl9k1giAI/c6xnaUUZFWRmh7fLiL2FyWHR4UybkoMRdkOTh0sp7a8keITDr/n2zEHgdpyF9VlLvKmVxB3TaTfY+3XNQW/+ISDy26ZRt6hih5ZPf2JiLwgCEMKu6jaxdMeJdutmzlXTuDs6VpL2C+7ZRp507sW34zr06gtd7HslmlUlTR0aa/Yr3vR0iQKjlZSkFXFFs9xirOrO1hFFxqpXSMIwpDCLqqdZarYFymFR4Wyat0sFt6QxtiJ0YRFhlhRd2dZLqUna6guc1kC/+WmIvb876lus2LCo0IZPyUWgISUUaSmx1OQVTWgi6UkkhcEYUhiCnlzUwshtgwYV70bT1ML865OobmpBVe9m/CoUNzuFg59UkjR8SrGJEeRtaUEV52bZV+bBvifuDUHg8wP8wGsTJvDm4rQwJwrJ3SI0C++coJ1nNlP+5OAv0nj/kREXhCEAeN8BM83m8Vjy4C5aGmSNSlrRtHmAJCbeRaAMydrqa9qAiBnbxkVhfUsumkSmX/LpyCrqkPa5EVLk3DVuakorGfS3DEc3lTEXq/oh4wM6jCRap5j3pd9f2eTxv2JiLwgCANGT/PG7YOB7zmuerclyvZJ2Yzr02jxtOKqc/PlpiLqHU1We5ExodRXNdFQ7aah2k1FUR2N9R5S0+PRYD0hKIyl+iGhQRRlO8g7VIHbbWSDR8aNxFnnxlHWwIk9ZSiMKB5g4wtHKMp24GlqYdFXJ7e73+4mffsaEXlBEAaMnuaN+6Yp2j99Jz4BK02yOLua4uxqFt6Qxuzl4yk8UkXS9BjCI0MZNzWWFncLBUerqC1vZMLMOOKSIql3NBKdEEZjvZusLcZC/Ywb0liyZgoXLU1i/8enAWhwNHHok0IcpQ1WZB480ii2W5TtANpqufjLpb9QE7Ei8oIgDBg9zRu3i3dP8tD3/O8pCrKqGD8thrFp0bjdLdSUuaitaCR6TDjHd5ax8IY0Fn1tGq56N4c3FVF4vIpDnxRabVUWNzBhZhwxieF43G1rOUNCDSEfPy2G8dPjmL4okRMTy/C4W2huamH6okQ8TS1oYPqiRPZvPE1zUwuZH+Zb8wcXEhF5QRAGPaZ47994utuFR656NwVHKwGIHx9J+KhQ65zohDCa3UZ1dLe7hf0bTzNp7hhKc2s4c7IWgISJowgNDyIuMYKsLSVUn3VaHj4Ywl12upbLbplGWGQIx3aWMn1RopVbHzIyyLJozP7OuzqF1PR4PO4WaxK3u4VTfYWIvCAIQ4aLlibR3NSCp6mFLzcVWYJpfxr4clMRZXl1AFSXuVj01ck469zkH6ygptxlWTOhoUHWwiXTXgFITBtFbWUjnuZWgHYCr4ATe8ooyKrixETDh9/7YT6FR6soynYQmxhO0tQY9m883S5Lx9PUQkFWFWMnRlu2z4WqYyMiLwjCkCE8KpSQkYY4L7T55GBE8F9uKqLkhCHYMQnhzF2VwqfrjxKfFElNuYuo0SOpr2yixdPKtEWJBI8MImlqDDVlLmrKXcQkhFNZ0kBpTg3RCWEARMWPJCImhLGpMbjdLZzJrQGg5ISDsWnRAIxJiWJEsKIgq4ov3j9FcXY1hUeruObu2R0mh82o/ULVsRGRFwRhSOG74tUs32v63gDJM2NJSBnFtv85QW15I60ezZI1U3CcaSB75xlKT9aQs6cMDex5P4+achcjI4OpKXcRGjGCUaPDGJsWjUJRU+4ifnwkYVEhVvsAJTk1xI+PJDU9nuQZcVQU1jPvaqPCenF2NUXZDg5vKmLRVyf7nXvobIVuX1s3IvKCIAwpfAXTtD0mzIxj7tUphHonRs1c9pGRwUTFj6S5qYWasy7vtiBctuwZM8IHKD/dAEBdZSMAsYnhFGRVEZ8UyYSZcUTFh1J8vIa6ykaqy1wUZTusGjf1jkauv38O5YV1FGdXU5pbYy3GsuMr6v1p3YjIC4Iw5PBdnWoWHUuZFc8l10zEUdZAaW4NLZ5WSk/WkL3zDADzrk7BVeemuswQZZPImDaRT5w8ivpqNw1VTSROHkWSt0yBKdyp6fHUVTZaufgAkfGh6FZjDiDvUAVfuTvdWvR0bGdpp/XqwRD1/rRuROQFQRhy+IrkqnWzLNF31bvZ+voJirIdzL06hRHBiuLsasZNjSYkNIiVay/ii/dPEZsYzpiUKABKT1YDMG5qNGOSoyg7ZUT4IaHBHPykkNT0eEvg7QXOju0stSZtM25Ia1dewd4nX7rK9e9rROQFQRhy+It8m5taOLypiGZ3iyW8CkiaEktcYgT5hyvZ+2E+X24qpMnZQnF2NUvWTAGwsnFSZsRT6p1YDYsMZu4qY5DIuD6N5OlxTJo7xiofbF5z5tJxlORUkzo7nqTJsVZ/uhLuC1lXXkReEIQhhz9f3pwUTZoaAxiTr83uFg5+Usio0WFWKmSTs4XohDAmz0uwFi8569xUFtYzbVEiKbPjrTIHBz8toDi7mrETo7n0q5Pb5ekDZH6YT2xiOLXljWT+LZ+vPjjvgv0z6Cki8oIgDHnMOu7F2dWMnRRN8ow4FFDsTacMiwoGwggNDyI0LJgr/3EmeYcq2LUhl9LcahSKomwHJ/aUcfZ0rVXHJj4pkuLsapTtOmCsvD2xx1g1mzI7nsy/5XPZLdP89u1CV530RerJC4Iw5AmPCuUrd6ezZM0UFnxlIiEjg9j7YT5j06KNFMfpcdRVNjJ90Tiu+/bF5B2qYNLcMZbXXpTtYMJMY2AoyKoiNjGcy26ZxqzLxpOaHs+0RYlWHr6nqYUTe8rI/DDfyLOfHMtXH5xHXKL/N0eZ8wefrj/abT36/kAieUEQAgJ/hcrM6NlR1kBVaUO7UsGephZWrZvFxy9kGROzU2K4+MoJlObWUJTtIGdPGcEjgyjIqiJvegWAZQn5LsTqCnv2j79Mm/6mVyKvlPp/wM1AK3AWWKe1LvHu+xFwF9ACfFdr/XEv+yoIwjChtxaHr2efd6jCEmuzMmRJbjUXXzmBr9yd3u5a5jthNR0neJubWvC4W6x9Pelbd5k2/U1v7ZontNZztNbzgA+AnwAopWYB3wRmA9cCzyqlLmzpNUEQhiz21/f1BRctTbIi7zlXTrBsGvP1gJdcM9FaPasw0iHNtz6Z+8KjQrn0q5MJHxVK5of559Q3ezsXml5F8lrrWtvPSNrKJ98MvKG1bgLylFIngUXArt5cTxCE4UFfLw7yjew7i6yP7Sxl74f5LFkzpdNVqvaa8J0x0JOtdnrtySulfgH8E1ADXOndnAzsth1W5N3m7/x7gHsAUlNTe9sdQRACgP7OI++s/a4GF7uXb3/bkz9Bv1AVJntCt3aNUupTpVSWn7+bAbTWj2mtU4DXgAfN0/w0pf1sQ2v9vNY6Q2udkZCQcL73IQjCMMQsTtZXWSud2Squere1SMpXyPxZS3Z7aKDpNpLXWq/qYVt/AT4EfooRuafY9k0ASs65d4IgCF3QXcTcV7aJWb4gNT2eOd73uJr4i/4v5IrW7uhtds00rXWO9+dNQLb3+/vAX5RSTwHjgWnAnt5cSxAEwZfuvPu+sk18UzLtDCZB90dvPfnHlVIzMFIoTwPfBtBaH1FKvQkcBTzAA1rrls6bEQRBOHe6E9i+msAd7ELeFUprv1b5gJCRkaEzMzMHuhuCIAhDCqXUPq11hr99UtZAEAQhgBGRFwRBCGBE5AVBEAIYEXlBEIQARkReEAQhgBGRFwRBCGBE5AVBEAKYQZUnr5Qqx1hUdaEYA1RcwOv1J4F0LxBY9yP3MjgJpHuZqLX2W/xrUIn8hUYpldnZAoKhRiDdCwTW/ci9DE4C6V66QuwaQRCEAEZEXhAEIYAZ7iL//EB3oA8JpHuBwLofuZfBSSDdS6cMa09eEAQh0BnukbwgCEJAIyIvCIIQwAxLkVdKPaGUylZKfamU+qtSKta270dKqZNKqeNKqa8MYDd7hFLqFqXUEaVUq1Iqw2ffkLoXAKXUtd7+nlRKPTLQ/TlXlFIvKaXOKqWybNvilVKfKKVyvJ9xA9nHnqKUSlFKbVJKHfP+N/aQd/uQux+lVJhSao9S6pD3Xv6vd/uQu5dzZViKPPAJkK61ngOcAH4EoJSaBXwTmA1cCzyrlAoasF72jCxgDbDVvnEo3ou3f88A1wGzgFu99zGUWI/xz9vOI8BnWutpwGfe30MBD/ADrfVFwGLgAe+/j6F4P03ASq31XGAecK1SajFD817OiWEp8lrrjVprj/fnbowXjQPcDLyhtW7SWucBJ4FFA9HHnqK1Pqa1Pu5n15C7F4z+ndRan9Jau4E3MO5jyKC13gpU+Wy+GXjF+/0VYPWF7NP5orUu1Vrv936vA44ByQzB+9EG9d6fId4/zRC8l3NlWIq8D3cCH3m/JwOFtn1F3m1DkaF4L0Oxzz0hUWtdCoZwAmMHuD/njFIqDZgPfMEQvR+lVJBS6iBwFvhEaz1k7+Vc6O2LvActSqlPgXF+dj2mtX7Pe8xjGI+kr5mn+Tl+wHNMe3Iv/k7zs23A76UbhmKfAx6lVBTwDvAvWutapfz9axr8aK1bgHneObi/KqXSB7hLF4SAFXmt9aqu9iul1gI3AlfptsUCRUCK7bAJQEn/9LDndHcvnTAo76UbhmKfe0KZUipJa12qlErCiCSHBEqpEAyBf01rvcG7ecjeD4DWuloptRlj7mRI30tPGJZ2jVLqWuCHwE1aa6dt1/vAN5VSI5VSk4BpwJ6B6GMfMBTvZS8wTSk1SSkVijFx/P4A96kveB9Y6/2+Fujs6WtQoYyQ/UXgmNb6KduuIXc/SqkEM4tOKRUOrAKyGYL3cs5orYfdH8YkZCFw0Pv3R9u+x4Bc4Dhw3UD3tQf38g8YEXATUAZ8PFTvxdvn6zEynnIx7KgB79M59v91oBRo9v57uQsYjZG5keP9jB/ofvbwXi7DsMu+tP2/cv1QvB9gDnDAey9ZwE+824fcvZzrn5Q1EARBCGCGpV0jCIIwXBCRFwRBCGBE5AVBEAIYEXlBEIQARkReEAQhgBGRFwRBCGBE5AVBEAKY/w/lZLkuwd4NWQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_data([midp]*6, data, n_samples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So here's the definition of the gaussian kernel, which you may remember from high school...\n", " This person at the science march certainly remembered!\n", "\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def gaussian(d, bw): return torch.exp(-0.5*((d/bw))**2) / (bw*math.sqrt(2*math.pi))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_func(f):\n", " x = torch.linspace(0,10,100)\n", " plt.plot(x, f(x))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlbUlEQVR4nO3deXxU5b3H8c8vkz2QsCRsSVgEBCOL4IgI2latFRRLxValtS5dKHXfatFrr7e9t1attZV7qb1o1WqtFBVbXFqt+64EEGQ3BoEAkrATloQkv/vHjN4YggyQ5Ewm3/fLeSXnnOec+Z2X8J3Dmec8j7k7IiKSuJKCLkBERJqXgl5EJMEp6EVEEpyCXkQkwSnoRUQSXHLQBTQmNzfXe/fuHXQZIiKtxty5cze6e15j2+Iy6Hv37k1xcXHQZYiItBpmtmp/23TrRkQkwSnoRUQSnIJeRCTBKehFRBKcgl5EJMHFFPRmNsbMlptZiZlNaWT7QDN728yqzOz6Bts6mNnjZrbMzJaa2QlNVbyIiBzYAbtXmlkImAacBpQBc8xstrsvqddsM3Al8I1GDnE38E93/6aZpQKZh121iIjELJZ+9COAEncvBTCzGcB44LOgd/dyoNzMzqy/o5llA18CLo62qwaqm6TyRkx98UNCSUa7tGTapSXTPSedfl3akdc+DTNrrrcVEYlrsQR9PrCm3nIZcHyMxz8CqAAeMLOhwFzgKnff2bChmU0CJgH07NkzxsN/3h9e/Yhd1bX7rG+flsyg/BxG9e3MqH65DC3IITmkrydEpG2IJegbuxSOdbaSZGA4cIW7v2tmdwNTgJ/tc0D36cB0gHA4fEizoSz++elU1dSxs6qGHXtqWLt1NyXllZSUVzJ31RZ+868V/OZfK+iclcpZQ3tw9rB8hhTk6GpfRBJaLEFfBhTWWy4A1sV4/DKgzN3fjS4/TiTom4WZkZ4SIj0lROd2afTOzWJ0v9zPtm/eWc3bH23i2Q/W85f3VvPgWx8zsFt7fnjSEXz9mB6k6CpfRBJQLMk2B+hvZn2iX6aeD8yO5eDu/gmwxswGRFedSr17+y2tU1YqZw7pzrTvDGfOv32VX00YjDtc99gCvnzHyzz45kqqa+qCKk9EpFlYLHPGmtkZwO+AEHC/u//SzCYDuPsfzKwbUAxkA3VAJVDk7tvN7BjgPiAVKAUucfctX/R+4XDYW2pQM3fn5eXl/OGVUt77eDN9crO4cexATivqqls6ItJqmNlcdw83ui0eJwdvyaD/lLvzyooKfvnMUkrKKxnVtzO3nj2Y3rlZLVqHiMih+KKg103pKDPj5AFd+OdVJ/Gf44/mg7XbGHP3a9z3eim1dfH3YSgiEisFfQPJoSS+e0Jv/nXNlxndN5f/emYp5/7v25Rt2RV0aSIih0RBvx/dctK576Iwvz1vKCs+2cGZU9/ghSUbgi5LROSgKei/gJlx9rACnr7yRAo7ZfCDh4r55TNLqKlVzxwRaT0U9DHo1TmLJ348igtP6MW9r6/kkgfnsG333qDLEhGJiYI+RmnJIX4xfhB3nDOEtz/axITfv8mqTfuM5CAiEncU9Afp3OMKefj7x7NpZzXfmPYm81Z/4SMBIiKBU9AfghP6dubJS0eTnZHCBfe9y+sfVgRdkojIfinoD1Gf3Cwem3wCPTtl8r0H5/DsB+uDLklEpFEK+sPQpX06f510AkMLOnD5X+Yxa15Z0CWJiOxDQX+YcjJTePj7x3NC385c99gCnpyvsBeR+KKgbwIZqSHuu/A4RvbpzHUzF/D399cGXZKIyGcU9E0kIzXEHy8OM6JPJ6756/u6Zy8icUNB34QyU5O5/+LjGN6zI1fNmK/eOCISFxT0TSwzNZk/XnwcffPa8aOH5zJf/exFJGAxBb2ZjTGz5WZWYmb7TAVoZgPN7G0zqzKz6xvZHjKz+Wb2dFMUHe9yMlJ46HsjyG2XxiUPzqGkfEfQJYlIG3bAoDezEDANGAsUARPNrKhBs83AlcCd+znMVcDSw6iz1emSnc6fv388yUlJXHT/HMp37Am6JBFpo2K5oh8BlLh7qbtXAzOA8fUbuHu5u88B9hnpy8wKgDOJTCfYpvTsnMn9F4fZvLOaH/6pmN3VtUGXJCJtUCxBnw+sqbdcFl0Xq98BNxCZS7bNGVLQgakTh7Fw7TaumjFfs1WJSIuLJegbmyE7prQys3FAubvPjaHtJDMrNrPiiorE6q1yWlFXbhlXxPNLNnDbP9rUHSwRiQOxBH0ZUFhvuQBYF+PxRwNfN7OPidzyOcXM/txYQ3ef7u5hdw/n5eXFePjW4+LRfbgoOp79E3P19KyItJxYgn4O0N/M+phZKnA+MDuWg7v7je5e4O69o/u95O4XHHK1rdzN44oY1bczN876QMMbi0iLOWDQu3sNcDnwHJGeMzPdfbGZTTazyQBm1s3MyoBrgZvNrMzMspuz8NYoJZTEtG8Pp1tOOj96eC6fbFNPHBFpfuYef18OhsNhLy4uDrqMZrP8kx1M+P2b9O/anr/+aCRpyaGgSxKRVs7M5rp7uLFtejI2AAO6tefObw3l/TVb+a+n9eWsiDQvBX1Axg7uzqQvHcHD76zSOPYi0qwU9AG64fQBjDyiEzc9+QFL1m0PuhwRSVAK+gAlh5L474nDyclI4bK/zGPHnn0eLBYROWwK+oDltU/jvycOZ9Wmndz05CLi8ctxEWndFPRxYESfTlz3tQE8tWAdj7635sA7iIgcBAV9nPjxl/tyUv9c/uOpxbpfLyJNSkEfJ5KSjN+edwwdMlK4/NF57KquCbokEUkQCvo4ktsujd+ddwwrN+7kF08tCbocEUkQCvo4M6pfLpO/3JcZc9ZognERaRIK+jh07WlHMrQghylPLGTd1t1BlyMirZyCPg6lhJKYOnEYtXXO1X99X5OViMhhUdDHqV6ds/j5+EG8t3Iz018rDbocEWnFFPRx7Jzh+Ywd1I27/rWcRWu3BV2OiLRSCvo4ZmbcevZgOmamcs1f32fPXk0uLiIHT0Ef5zpmpfLrbw3lw/JKbv/nsqDLEZFWKKagN7MxZrbczErMbEoj2wea2dtmVmVm19dbX2hmL5vZUjNbbGZXNWXxbcWXj8zjohN68cCbH/PWRxuDLkdEWpkDBr2ZhYBpwFigCJhoZkUNmm0GrgTubLC+BrjO3Y8CRgKXNbKvxGDK2KPok5vFTx5bqFEuReSgxHJFPwIocfdSd68GZgDj6zdw93J3nwPsbbB+vbvPi/6+g8ics/lNUnkbk5Ea4s5vDWX9tt388hnNSiUisYsl6POB+kMqlnEIYW1mvYFhwLv72T7JzIrNrLiiouJgD98mHNurIz+KPjX78rLyoMsRkVYilqC3RtYd1BM8ZtYOeAK42t0bHZrR3ae7e9jdw3l5eQdz+Dbl6q/2Z2C39vz0iYVs3VUddDki0grEEvRlQGG95QJgXaxvYGYpREL+EXefdXDlSUNpyZFbOJt3VmvgMxGJSSxBPwfob2Z9zCwVOB+YHcvBzcyAPwJL3f2uQy9T6huUn8OlJ/dj1vy1vLBkQ9DliEicO2DQu3sNcDnwHJEvU2e6+2Izm2xmkwHMrJuZlQHXAjebWZmZZQOjge8Cp5jZ+9HXGc12Nm3I5Sf3Y2C39tz45Ae6hSMiX8jicY7ScDjsxcXFQZcR9xat3cY3pr3J14f24K7zjgm6HBEJkJnNdfdwY9v0ZGwrVv8WzkvLdAtHRBqnoG/lLj+5HwO6tuemWYvYrgepRKQRCvpWLjU5iTu+OYTyHXv41bN6kEpE9qWgTwBDCzvww5OO4NH31vBmicbCEZHPU9AniGtOO5I+uVn89ImF7KyqCbocEYkjCvoEkZ4S4vZzhlC2ZTe/eX5F0OWISBxR0CeQEX06ccHInjzw1krmrd4SdDkiEicU9Anmp2MG0i07nSlPLKS6pi7ockQkDijoE0z79BR+efYgVmyo5PevlARdjojEAQV9AjplYFe+PrQH014uYcWGHUGXIyIBU9AnqFvOKqJdWjJTnlhIXV38DXMhIi1HQZ+gOrdL42fjipi3eit/fndV0OWISIAU9Ans7GH5nNQ/l9v/sYx1W3cHXY6IBERBn8DMjFvPHkydw8/+toh4HKlURJqfgj7BFXbK5LqvHcmLy8p5euH6oMsRkQAo6NuAS0b3YUhBDj9/arEmKRFpg2IKejMbY2bLzazEzKY0sn2gmb1tZlVmdv3B7CvNL5Rk3DZhCFt27eVWjXAp0uYcMOjNLARMA8YCRcBEMytq0GwzcCVw5yHsKy2gqEc2PzzpCGYWl/HWRxrhUqQtieWKfgRQ4u6l7l4NzADG12/g7uXuPgdoOPPFAfeVlnP1V/vTq3MmN836gD17a4MuR0RaSCxBnw+sqbdcFl0Xi5j3NbNJZlZsZsUVFRUxHl4ORnpKiFvPHszHm3Yx9cUPgy5HRFpILEFvjayLtZ9ezPu6+3R3D7t7OC8vL8bDy8Ea3S+Xc4YXMP21Upau3x50OSLSAmIJ+jKgsN5yAbAuxuMfzr7STG4+8yiyM1K4cdYH1Gp4BJGEF0vQzwH6m1kfM0sFzgdmx3j8w9lXmknHrFT+fVwR76/ZysNvfxx0OSLSzA4Y9O5eA1wOPAcsBWa6+2Izm2xmkwHMrJuZlQHXAjebWZmZZe9v3+Y6GYnd+GN68KUj8/j1c8s1PIJIgrN4fCw+HA57cXFx0GUkvDWbd/G1377G6H6duffCMGaNfaUiIq2Bmc1193Bj2/RkbBtW2CmTa087kheWlvOPRZ8EXY6INBMFfRt3yejeDMrP5pbZi9m2u+FjECKSCBT0bVxyKInbJgxhU2UVt/1jWdDliEgzUNALg/Jz+P6JfXj0vdW8W7op6HJEpIkp6AWAa047koKOGdz4pIZHEEk0CnoBIDM1mV+ePZjSip38/uWSoMsRkSakoJfPfPnIPM4els89r37E8k92BF2OiDQRBb18zs1nHkW7tGSmzFqo4RFEEoSCXj6nc7s0/v2sIuav1vAIIolCQS/7+MYx+Z8Nj7BWwyOItHoKetmHmfHLbwyizuHmJz8gHofJEJHYKeilUYWdMrn+9AG8vLyC2Qs0srRIa6agl/26eFRvhhZ24OdPLWHzzuqgyxGRQ6Sgl/0KJRm3nzOY7bv38p9PLwm6HBE5RAp6+UIDu2Vz6Vf68uT8tby8vDzockTkECjo5YAuO6Uf/bq0499mfUBlVU3Q5YjIQYop6M1sjJktN7MSM5vSyHYzs6nR7QvNbHi9bdeY2WIzW2Rmj5pZelOegDS/tOQQt58zhPXb93DHPzXCpUhrc8CgN7MQMA0YCxQBE82sqEGzsUD/6GsScE9033zgSiDs7oOAEJF5Y6WVObZXRy4e1ZuH3l7Feys3B12OiByEWK7oRwAl7l7q7tXADGB8gzbjgYc84h2gg5l1j25LBjLMLBnIBNRXr5W6/msDKOiYwZQnFmqES5FWJJagzwfW1Fsui647YBt3XwvcCawG1gPb3P35xt7EzCaZWbGZFVdUVMRav7SgrLRkbpswhNKNO/ntCyuCLkdEYhRL0Dc2Y3TDRyUbbWNmHYlc7fcBegBZZnZBY2/i7tPdPezu4by8vBjKkiCc2D+X848r5N7XSlmwZmvQ5YhIDGIJ+jKgsN5yAfveftlfm68CK929wt33ArOAUYdersSDm848ii7t0/nJ4wuoqtEtHJF4F0vQzwH6m1kfM0sl8mXq7AZtZgMXRnvfjCRyi2Y9kVs2I80s08wMOBVY2oT1SwCy01O4dcIgVmyoZNpLmqREJN4dMOjdvQa4HHiOSEjPdPfFZjbZzCZHmz0LlAIlwL3ApdF93wUeB+YBH0Tfb3pTn4S0vFMGdmXCsHx+/8pHLFq7LehyROQLWDyOTBgOh724uDjoMuQAtu3ay2m/fZVOWanMvvxEUpP1/J1IUMxsrruHG9umv5lyyHIyU7j17MEs+2QH//PSh0GXIyL7oaCXw/LVosgtnGm6hSMStxT0cthuOetoOmelct1M9cIRiUcKejlsOZkp3HbOYJZv2MHdL+gWjki8UdBLkzhlYFfODRfwh1c/Yt7qLUGXIyL1KOilyfxsXBHdczK4fuYCdlfrFo5IvFDQS5Npn57Cr78ZGQvndg1nLBI3FPTSpEb1y+XiUb158K2PebNkY9DliAgKemkGPx0zkCPysrj+sQVs27036HJE2jwFvTS5jNQQvz33GMp3VHHL3xcFXY5Im6egl2YxtLADV5zSj7+9v46nFmiuGZEgKeil2Vx2cj+GFnbg5r8tYv223UGXI9JmKeil2aSEkvjdecdQXVPHdTMXUFcXfwPoibQFCnppVn1ys/iPrxfx1kebuO+N0qDLEWmTFPTS7M4NFzLm6G78+rnlGvhMJAAxBb2ZjTGz5WZWYmZTGtluZjY1un2hmQ2vt62DmT1uZsvMbKmZndCUJyDxz8z41YTBdMpK5aoZ8/XUrEgLO2DQm1kImAaMBYqAiWZW1KDZWKB/9DUJuKfetruBf7r7QGAomkqwTeqYlcpd5x5D6cad/OLpxUGXI9KmxHJFPwIocfdSd68GZgDjG7QZDzzkEe8AHcysu5llA18C/gjg7tXuvrXpypfWZHS/XH785b48+t4anlm4PuhyRNqMWII+H1hTb7ksui6WNkcAFcADZjbfzO4zs6zDqFdauWtOO5JhPTswZdZC1mzeFXQ5Im1CLEFvjaxr2E9uf22SgeHAPe4+DNgJ7HOPH8DMJplZsZkVV1RUxFCWtEYpoSSmnj8MHK6cMZ+9tXVBlySS8GIJ+jKgsN5yAdDwUcf9tSkDytz93ej6x4kE/z7cfbq7h909nJeXF0vt0koVdsrkV+cMZv7qrdz5/PKgyxFJeLEE/Rygv5n1MbNU4HxgdoM2s4ELo71vRgLb3H29u38CrDGzAdF2pwJLmqp4ab3GDenBd47vyf++WspLyzYEXY5IQjtg0Lt7DXA58ByRHjMz3X2xmU02s8nRZs8CpUAJcC9wab1DXAE8YmYLgWOAW5uufGnNfjauiKO6Z3PtzAWs26ohEkSai7nH32Pp4XDYi4uLgy5DWsDKjTsZN/V1BnRrz19/dAIpIT3DJ3IozGyuu4cb26a/VRKoPrlZ3HbOEOat3srt/9CsVCLNQUEvgTtraA8uHtWb+95YybMfqH+9SFNT0EtcuOmMoxjWswM/eWwBH1VUBl2OSEJR0EtcSE1O4vffGU56SojJD89lZ1VN0CWJJAwFvcSN7jkZTJ04jI8qKrnhiYXEY0cBkdZIQS9xZXS/XH46ZiDPLFzPH17V+PUiTUFBL3Fn0peO4KyhPbjjuWW8ukLDYYgcLgW9xB0z4/ZzBjOga3uu+Ms8Vm3aGXRJIq2agl7iUmZqMvdeGCYpyfj+n4rZsWdv0CWJtFoKeolbhZ0yuec7x/Lxxp1c+eh8ajW5uMghUdBLXDuhb2d+MX4QLy+v4LZ/aHIykUORHHQBIgfy7eN7smLDDu59fSX9urTjvON6Bl2SSKuioJdW4eYzj6J0407+7clF5HfI5MT+uUGXJNJq6NaNtArJoSSmfXsY/bq048d/nsvyT3YEXZJIq6Ggl1ajfXoK9198HBmpIb734BzKt+8JuiSRVkFBL61Kjw4Z3H/xcWzZVc0lD86hUmPiiBxQTEFvZmPMbLmZlZjZPpN7R6cQnBrdvtDMhjfYHjKz+Wb2dFMVLm3XoPwcpn1nOMs+2cHkh+dSXaMJxkW+yAGD3sxCwDRgLFAETDSzogbNxgL9o69JwD0Ntl9FZBpCkSZx8oAu3DZhMG+UbOSGxxdQpz72IvsVyxX9CKDE3UvdvRqYAYxv0GY88JBHvAN0MLPuAGZWAJwJ3NeEdYvwrXAhPzl9AH97fx23PrtUo12K7Ecs3SvzgTX1lsuA42Nokw+sB34H3AC0/6I3MbNJRP41QM+e6ictsbn0K32p2FHFfW+sJCcjhStO7R90SSJxJ5YremtkXcNLp0bbmNk4oNzd5x7oTdx9uruH3T2cl5cXQ1kikQHQ/n1cEROG5/Obf63gwTdXBl2SSNyJ5Yq+DCist1wArIuxzTeBr5vZGUA6kG1mf3b3Cw69ZJHPS0oy7jhnCJV7aviPp5bQLj2Fbx5bEHRZInEjliv6OUB/M+tjZqnA+cDsBm1mAxdGe9+MBLa5+3p3v9HdC9y9d3S/lxTy0hySQ0lMnTiME/vlcsPjC/j7+2uDLkkkbhww6N29BrgceI5Iz5mZ7r7YzCab2eRos2eBUqAEuBe4tJnqFdmv9JQQ914Y5rjenbh25gKe/WB90CWJxAWLx54K4XDYi4uLgy5DWqmdVTVcdP97vL9mK9O+M5zTj+4WdEkizc7M5rp7uLFtejJWEk5WWjIPXHIcgwtyuOyRebqylzZPQS8JqX16Cg99bwRDCztwxaPzdc9e2jQFvSSsT8M+3KsjV//1fR4rXnPgnUQSkIJeElpWWjIPXjKC0X1z+cnjC7n/DfWzl7ZHQS8JLyM1xH0XhTn96K784ukl3PWvFRouQdoUBb20CekpIaZ9ezjfOraAqS9+yC2zF2uycWkzNJWgtBnJoSTu+OYQOmalMv21UjZs38Pd5w8jPSUUdGkizUpX9NKmmBk3nXEUPxtXxPNLNvDte99h887qoMsSaVYKemmTvn9iH37/7eEsWredc+55i9KKyqBLEmk2Cnpps8YO7s5ffnA823bv5RvT3uTNko1BlyTSLBT00qaFe3fi75eNpltOOhfe/x5/fmdV0CWJNDkFvbR5hZ0yeeLHo/hS/1xu/tsibpy1kKqa2qDLEmkyCnoRIk/R3nfRcVx2cl8efW8N5/7vO6zftjvoskSahIJeJCqUZPzk9IH84YJj+ai8knFT3+D1DyuCLkvksCnoRRoYM6gbf7tsNJ3bpXLh/e9x53PLqamtC7oskUOmoBdpRL8u7fj7ZSfyrWML+J+XS/j2ve+ybqtu5UjrFFPQm9kYM1tuZiVmNqWR7WZmU6PbF5rZ8Oj6QjN72cyWmtliM7uqqU9ApLlkpIa445tDuevcoSxet40xv3uN2QsaTpcsEv8OGPRmFgKmAWOBImCimRU1aDYW6B99TQLuia6vAa5z96OAkcBljewrEtcmDC/g2atOol+Xdlz56HyunjGfbbv2Bl2WSMxiuaIfAZS4e6m7VwMzgPEN2owHHvKId4AOZtY9OkH4PAB330Fkztn8JqxfpEX06pzFzB+dwLWnHclTC9dz2m9f5YUlG4IuSyQmsQR9PlB/xoYy9g3rA7Yxs97AMODdxt7EzCaZWbGZFVdUqKeDxJ/kUBJXntqfv106mk5ZqfzgoWKunjFfY+VI3Isl6K2RdQ3Hd/3CNmbWDngCuNrdtzf2Ju4+3d3D7h7Oy8uLoSyRYAwuyGH25Sdy1an9eXrhek79zSvMLF6jMe4lbsUS9GVAYb3lAqDhN1L7bWNmKURC/hF3n3XopYrEj9TkJK457UieufIkjshrxw2PL+S86e+wYsOOoEsT2UcsQT8H6G9mfcwsFTgfmN2gzWzgwmjvm5HANndfb2YG/BFY6u53NWnlInFgQLf2PPajE7htwmCWf7KDsXe/zi1/X8TWXbqdI/HjgEHv7jXA5cBzRL5Mnenui81ssplNjjZ7FigFSoB7gUuj60cD3wVOMbP3o68zmvokRIKUlGScP6InL1//FSaOKOThd1bxlTtf4YE3V1JdowetJHgWj/cVw+GwFxcXB12GyCFZun47//n0Et76aBM9O2Xyk9MHcObg7iQlNfZVlkjTMLO57h5ubJuejBVpYkd1z+aRHxzPA5ccR2ZqiCsenc9Z//MGLyzZoC9sJRAKepFmYGacPKALz1x5EnedO5TKqhp+8FAx46e9qcCXFqdbNyItYG9tHU/OX8vUFz+kbMtuBnRtz4+/0pdxQ7qTHNL1lhy+L7p1o6AXaUF7a+t4asE67nnlIz4sryS/QwYXjerFecf1JCcjJejypBVT0IvEmbo658Vl5fzxjVLeKd1MZmqICcPzuWBkLwZ2yw66PGmFFPQicWzxum3c/8bHPLVwHdU1dRzbqyMTR/TkjMHdyExNDro8aSUU9CKtwJad1Twxr4xH3l3Nyo07yUoNceaQ7kwYXsCI3p3UPVO+kIJepBVxd+Z8vIXH567hmYXr2VldS/ecdMYN6c5ZQ3swOD+HyEPnIv9PQS/SSu2squGFpRt4asE6Xl1Rwd5aJ79DBqcf3Y3Tj+7Ksb06qteOAAp6kYSwdVc1zy/ZwPOLP+G1DzdSXVNHTkYKXxmQxykDu3Biv1w6t0sLukwJiIJeJMFUVtXw6vIKXlpWzivLy9kUHRO/qHs2J/XPZWTfzoR7daR9urpsthUKepEEVlfnfLB2G2+UbOT1DyuYu2oLe2udUJIxqEc2x/bqRLh3R47t1ZGu2elBlyvNREEv0obsrq5l3uotvFu6iXdWbmbBmq1URUfR7JadzpCCHIYWduDoHtkU9cimS3uFfyL4oqBXJ12RBJORGmJ0v1xG98sFoLqmjqXrtzN31RYWlG1lYdk2nq83321e+zQGdmvPgK7tObJbe/p3accRee30pG4CUdCLJLjU5CSGFnZgaGGHz9Zt272Xpeu3s2Tddhav286KDTt4+J1Vn135Q+QDoE9uFr07Z9I7N4uenTIp7JhJYadMOmamqItnK6KgF2mDcjJSGHlEZ0Ye0fmzdbV1zurNuygpr+SjikpKyitZtWknLy+voKK47HP7Z6SE6NEhnR4dMuiek0637HS65qTTtX06ee3T6JKdRuesNFKT1fUzHsQU9GY2BrgbCAH3ufttDbZbdPsZwC7gYnefF8u+IhIfQklGn9ws+uRmcRpdP7etsqqGNZt3UbZlN2s272Lt1t2si76WfbKDjZVVNPZ1X3Z6Mp3bpdEpK5WOmal0zEyhY1YqORkpn72yM1LITk+mfXoK7dOTaZeWTGZqSP9iaEIHDHozCwHTgNOITAI+x8xmu/uSes3GAv2jr+OBe4DjY9xXROJcu7RkjuqezVHdGx9wraa2jorKKjZsr2LjjioqKquo2FHF5p3VbKysYlNlNWVbdrFo7V627Kr+3C2ixphBVmoyWWkhslKTyUgNkZkaIiM1mYyUJDJSQqRHX2kpSaQnR36mJYdITU4iLZREanL0FUoiOWSkhpJISU4iOclIia5LToosf/p7KMlITjKSoj9DSUaSffqTVvvhE8sV/QigxN1LAcxsBjAeqB/W44GHPNKF5x0z62Bm3YHeMewrIq1cciiJ7jkZdM/JiKn9nr21bNu9l22797J991527Klh+569VFbVULmnhsqqGnZW1bKzqobK6hr2VNeyqzqyz4ZtteypiSxX7a1lT01di83Nm2SQZJEPgiSDkEU+CMyIrvv/D4QkAyO6LfoBYcZny0aknQEYGNA5K42Zk09o8rpjCfp8YE295TIiV+0HapMf474AmNkkYBJAz549YyhLRFqrT6/Gm6pff12dU11bR1VNHVU1tVTX1LG31qmOfgjsratjb00d1bV11NQ6e2vrqKnzyCu6rtYjy7W1ddQ61NbVUVsHde7URtt69Pdad9wj7/vp7+711ke3Of+/7DjR/6j7dB8i+zmAQ3ZG83xtGstRG/u3SsO7cftrE8u+kZXu04HpEOlHH0NdIiJA5Go6PSny4QHqFtpQLEFfBhTWWy4A1sXYJjWGfUVEpBnF0vdpDtDfzPqYWSpwPjC7QZvZwIUWMRLY5u7rY9xXRESa0QGv6N29xswuB54j0kXyfndfbGaTo9v/ADxLpGtlCZHulZd80b7NciYiItIojXUjIpIAvmisGz22JiKS4BT0IiIJTkEvIpLgFPQiIgkuLr+MNbMKYNUh7p4LbGzCcloDnXPia2vnCzrng9XL3fMa2xCXQX84zKx4f988Jyqdc+Jra+cLOuempFs3IiIJTkEvIpLgEjHopwddQAB0zomvrZ0v6JybTMLdoxcRkc9LxCt6ERGpR0EvIpLgEibozWyMmS03sxIzmxJ0Pc3NzArN7GUzW2pmi83sqqBrailmFjKz+Wb2dNC1tITo1JyPm9my6P/vpp9rLs6Y2TXRP9eLzOxRM2uaqajiiJndb2blZrao3rpOZvYvM/sw+rNjU7xXQgR9vUnIxwJFwEQzKwq2qmZXA1zn7kcBI4HL2sA5f+oqYGnQRbSgu4F/uvtAYCgJfu5mlg9cCYTdfRCRIc7PD7aqZvEgMKbBuinAi+7eH3gxunzYEiLoqTeBubtXA59OQp6w3H29u8+L/r6DyF/+/GCran5mVgCcCdwXdC0twcyygS8BfwRw92p33xpoUS0jGcgws2QgkwScmc7dXwM2N1g9HvhT9Pc/Ad9oivdKlKDf3+TkbYKZ9QaGAe8GXEpL+B1wA1AXcB0t5QigAnggervqPjPLCrqo5uTua4E7gdXAeiIz1j0fbFUtpmt0dj6iP7s0xUETJehjnoQ80ZhZO+AJ4Gp33x50Pc3JzMYB5e4+N+haWlAyMBy4x92HATtpon/Ox6vofenxQB+gB5BlZhcEW1XrlihBH8sE5gnHzFKIhPwj7j4r6HpawGjg62b2MZHbc6eY2Z+DLanZlQFl7v7pv9YeJxL8ieyrwEp3r3D3vcAsYFTANbWUDWbWHSD6s7wpDpooQd/mJiE3MyNy33apu98VdD0twd1vdPcCd+9N5P/xS+6e0Fd67v4JsMbMBkRXnQosCbCklrAaGGlmmdE/56eS4F9A1zMbuCj6+0XA35vioAecHLw1aKOTkI8Gvgt8YGbvR9fd5O7PBleSNJMrgEeiFzGlwCUB19Os3P1dM3scmEekd9l8EnA4BDN7FPgKkGtmZcAtwG3ATDP7PpEPvG81yXtpCAQRkcSWKLduRERkPxT0IiIJTkEvIpLgFPQiIglOQS8ikuAU9CIiCU5BLyKS4P4PkayhaGV7dL4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_func(partial(gaussian, bw=2.5))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "functools.partial" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "partial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our implementation, we choose the bandwidth to be 2.5. \n", "\n", "One easy way to choose bandwidth is to find which bandwidth covers one third of the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def tri(d, i): return (-d+i).clamp_min(0)/i" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgJElEQVR4nO3deXhU9b3H8fc3Cwk7CAFZZREEZBEYIJC421ZqCbuCggoKIgm1e217q7V6r72t9VpN2ERQFEUFhGBdam2tTUKAsO8QWSNbQEAW2X/3D2KfiAEGmMmZ5fN6njzJnHMy8xkSPgxnZr4/c84hIiLhL8brACIiEhgqdBGRCKFCFxGJECp0EZEIoUIXEYkQcV7dcO3atV2TJk28unkRkbC0aNGiPc65pLL2eVboTZo0oaCgwKubFxEJS2a25Vz7dMpFRCRCqNBFRCKECl1EJEKo0EVEIoQKXUQkQlyw0M1sspntNrOV59hvZva8mRWa2XIz6xT4mCIiciH+PEJ/Gbj9PPt7Ai1KPkYC4y4/loiIXKwLFrpz7lPgi/Mc0huY6s7IB2qYWb1ABTzbF4eP88TcVRw5fjJYNyEiEpYCcQ69AbCt1OWikm3fYmYjzazAzAqKi4sv6cZyC/fwct5mBoybR9G+I5d0HSIikSgQhW5lbCtz1Qzn3ETnnM8550tKKvOdqxfUq0N9Jt/XhW37jpCWmcv8jXsv6XpERCJNIAq9CGhU6nJDYHsArvecbm5Vh9npKdSoFM89k+bzav4WtPKSiES7QBR6NnBvyatdkoEDzrkdAbje82qeVIXZ6Snc0DKJ385eya/fWcnxk6eDfbMiIiHrgsO5zOwN4CagtpkVAY8D8QDOufHAe8D3gULgCDAsWGHPVi0xnhfv9fHnv61j7CefsWHXQcYN6UxS1YTyiiAiEjLMq1MVPp/PBXLa4txl2/n5jGXUrFSBiUN9tGtYPWDXLSISKsxskXPOV9a+iHmnaK8O9ZkxqgcxZgwYn8ecpZ97HUlEpFxFTKEDtG1QneyMFDo0qsEj05fy9HtrOHVaT5aKSHSIqEIHqFUlgWkPdmNIcmMmfLqR4S8v5MCRE17HEhEJuogrdID42Bie6tOOp/u1I++zPfQZm0vh7oNexxIRCaqILPSvDe7amDdGJHPw6An6ZOXx99W7vI4kIhI0EV3oAL4mV5CdkUrT2pUZ8WoBWf8s1JuQRCQiRXyhA9SvUZG3R3UnrUN9/vThOjJeX6LhXiIScaKi0AES42N57q7r+FXPVry/cgf9x81j2xca7iUikSNqCh3AzHjoxuZMvr8LRfuOkJaZw7zPNNxLRCJDVBX61266pg7ZGanUqpLAkJfmM3XeZp1XF5GwF5WFDtC0dmXeGd2Dm69J4rE5q3h05gqOnTzldSwRkUsWtYUOUDUxnolDfYy55WreLNjG4In57D541OtYIiKXJKoLHSAmxvjpd68h6+5OrNlxkLQXclm2bb/XsURELlrUF/rX7mhfj5kP9yA2xhg4YR6zFhd5HUlE5KKo0EtpU78a2RkpdGpcg5+8tYz//utqTp7SohkiEh5U6GepVSWBVx/oxn3dr+LFf29imIZ7iUiYUKGXIT42hid6t+UP/dqRv3EvaVk5bNil4V4iEtpU6OcxqGtjpo9M5vCxU/TJyuVvq3Z6HUlE5JxU6BfQ+aormDsmheZ1qjDy1UU8//EGTmvRDBEJQSp0P9SrXpG3HupOv44NePaj9aS/vpjDxzTcS0RCiwrdT4nxsfz5zg781x2t+XDVTvqPy2PrXg33EpHQoUK/CGbGg9c345XhXdlx4ChpWTnkFu7xOpaICKBCvyTXt0hiTnoKSVUSuHfyAqbkbtJwLxHxnAr9EjWpXZl30lO4pVUdnpi7ml/MWK7hXiLiKRX6ZaiSEMeEIZ354a0teHtREXdNyGfXlxruJSLeUKFfppgY4yffacn4IZ1Yv+sgvV7IYcnWfV7HEpEopEIPkNvb1mPW6B4kxMdw14R83i7Y5nUkEYkyKvQAanVlNbLTU+nStCY/n7GcJ+au0nAvESk3KvQAq1m5Aq8M68rwlKZMyd3MvZMXsO/wca9jiUgUUKEHQVxsDI/1asOfBrSnYPM+0rJyWLvzS69jiUiEU6EH0UBfI6Y/lMyxE6fpNzaPD1bu8DqSiEQwvwrdzG43s3VmVmhmj5axv7qZzTWzZWa2ysyGBT5qeOrUuCZzx6TSsm5VRr22mGc/Wq/hXiISFBcsdDOLBbKAnkAbYLCZtTnrsHRgtXOuA3AT8GczqxDgrGGrbrVEpo9Mpn+nhjz/8QYeem0RhzTcS0QCzJ9H6F2BQufcRufccWA60PusYxxQ1cwMqAJ8AaixSkmMj+WZge157Adt+Mfa3fQbm8vmPYe9jiUiEcSfQm8AlH5RdVHJttIygdbAdmAF8Ihz7luv1zOzkWZWYGYFxcXFlxg5fJkZw1ObMnV4V3YfPEZaZg6fro++PwcRCQ5/Ct3K2Hb2SeDvAUuB+sB1QKaZVfvWNzk30Tnnc875kpKSLjJq5Ei5ujbZ6anUq16R+6csYNK/N2q4l4hcNn8KvQhoVOpyQ848Ei9tGDDLnVEIbAJaBSZiZGpcqxKzRvfgu22u5Km/ruGnby3j6AkN9xKRS+dPoS8EWphZ05InOgcB2WcdsxW4FcDM6gLXABsDGTQSVU6IY+w9nfjxbS2ZteRz7powj50HNNxLRC7NBQvdOXcSyAA+BNYAbznnVpnZKDMbVXLYk0APM1sBfAz80jmnlR/8EBNjPHJbCyYO7Uzh7kP0ysxh0RYN9xKRi2denbv1+XyuoKDAk9sOVet3HWTE1AJ27D/KU33acmeXRhf+JhGJKma2yDnnK2uf3ikaQlrWrcqc9BS6NbuCX8xczuNzVnJCw71ExE8q9BBTo1IFptzfhRHXN+WVeVsY+tJ8vtBwLxHxgwo9BMXFxvCbO9rw7J0dWLx1P2mZOazeruFeInJ+KvQQ1q9TQ95+qDsnTp2m/7g83luh4V4icm4q9BDXoVEN5mak0rpeVUZPW8yf/7ZOw71EpEwq9DBQp1oib4xM5i5fI174RyEjXy3g4NETXscSkRCjQg8TCXGx/KF/O55Iu5Z/rium79g8Nmm4l4iUokIPI2bGfT2a8OoDXdl76Bi9M3P4l4Z7iUgJFXoY6tG8NtkZqdSvUZFhUxYw8dPPNNxLRFTo4arRFWeGe93e9kr+5721/PjNpRruJRLlVOhhrFKFOLLu7sTPvtuSOcu2M3D8PLbv/8rrWCLiERV6mDMzMm5pwYtDfWzac5i0zBwKNn/hdSwR8YAKPULc1qYus9N7UCUhjsEv5vPGgq1eRxKRcqZCjyBX16nKnPRUujevza9mreC3szXcSySaqNAjTPVK8Uy5vwsjb2jGq/lbGDJpPnsPHfM6loiUAxV6BIqNMX79/dY8d9d1LN22n7TMXFZtP+B1LBEJMhV6BOvTsQFvj+rOaefoPy6Pd5efvRSsiEQSFXqEa9+wBtkZqbStX52M15fwxw/WckrDvUQikgo9CiRVTeD1EckM7tqIsZ98xoipBXyp4V4iEUeFHiUqxMXwP33b8WTva/l0fTF9snLZWHzI61giEkAq9ChiZgzt3oTXHuzG/iMn6J2Vyz/X7fY6logEiAo9CiU3q0V2RgqNalZi+MsLGf8vDfcSiQQq9CjVsGYlZjzcnTva1eMP76/lkelL+eq4hnuJhLM4rwOIdypViOOFwR1pXa8az/xtHRv3HGLCUB8NalT0OpqIXAI9Qo9yZkb6zVfz0n0+tuw5Qu/MHBZs0nAvkXCkQhcAbmlVl3fSU6iWGM/dL+Yzbf4WryOJyEVSoct/XF2nCu+kp5Daoja/eWclv3lnBcdPariXSLhQocs3VK8Yz0v3dWHUjc2ZNn8rQybNZ4+Ge4mEBRW6fEtsjPFoz1b8ZdB1LP98P2kv5LDycw33Egl1KnQ5p97XNWDGqB4ADBifR/YyDfcSCWV+FbqZ3W5m68ys0MwePccxN5nZUjNbZWb/CmxM8UrbBtXJHpNKuwbV+eEbS/hfDfcSCVkXLHQziwWygJ5AG2CwmbU565gawFggzTl3LTAw8FHFK7WrJDDtwWTu6daYcZ98xgOvLOTAVxruJRJq/HmE3hUodM5tdM4dB6YDvc865m5glnNuK4BzTgNCIkyFuBj+u287nurTlpwNe+iblUvhbg33Egkl/hR6A2BbqctFJdtKawnUNLNPzGyRmd1b1hWZ2UgzKzCzguLi4ktLLJ4aknwVr49I5sBXJ+iblcvHa3Z5HUlESvhT6FbGtrNPosYBnYE7gO8BvzWzlt/6JucmOud8zjlfUlLSRYeV0NC16RVkj0mlca1KPDi1gKx/Fmq4l0gI8KfQi4BGpS43BM5+uUMR8IFz7rBzbg/wKdAhMBElFDWoUZEZo3rQq319/vThOjLeWMKR4ye9jiUS1fwp9IVACzNramYVgEFA9lnHzAGuN7M4M6sEdAPWBDaqhJqKFWL5y6DreLRnK95bsYP+4+ZRtO+I17FEotYFC905dxLIAD7kTEm/5ZxbZWajzGxUyTFrgA+A5cACYJJzbmXwYkuoMDNG3dicyfd3oWjfEdIyc8nfuNfrWCJRybw69+nz+VxBQYEnty3BsbH4ECOmFrBl7xEe79WGIclXYVbWUzAicqnMbJFzzlfWPr1TVAKmWdKZ4V43tkzit3NW8atZKzh2UotmiJQXFboEVLXEeCbe6yP95uZMX7iNu1+cz+6DR72OJRIVVOgScLExxs+/14rMuzuyevuXpL2Qy7Jt+72OJRLxVOgSND9oX58ZD3cnNsYYOGEe7ywp8jqSSERToUtQXVu/OtkZKXRsVIMfv7mM//7rak6e0qIZIsGgQpegq1Ulgdce7MZ93a/ixX9vYtjLCzlwRMO9RAJNhS7lIj42hid6t+UP/dqRv3EvvbNy2LDroNexRCKKCl3K1aCujZk+MplDx07RJyuXj1ZruJdIoKjQpdx1vuoK5o5JoXmdKoyYWsDzH2/gtBbNELlsKnTxRL3qFXnroe707diAZz9aT/rrizl8TMO9RC6HCl08kxgfy7N3duA332/Nh6t20n9cHtu+0HAvkUulQhdPmRkjbmjGy8O6sn3/V6Rl5pBXuMfrWCJhSYUuIeGGlklkZ6RSu0oCQycv4OXcTVo0Q+QiqdAlZDSpXZlZo3twS6s6/G7uan45c7mGe4lcBBW6hJSqifFMGNKZH95yNW8VFDFoYj67v9RwLxF/qNAl5MTEGD/57jWMu6cT63YepFdmDks13EvkglToErJ6tqvHzId7EB8bw50T5jFjkYZ7iZyPCl1CWut61cjOSKVz45r87O1lPPmuhnuJnIsKXULeFZUrMPWBrgxLacJLOZu4f8pC9h857nUskZCjQpewEB8bw+O9ruWPA9qzYNMXpGXmsm6nhnuJlKZCl7Byp68R0x9K5qsTp+g7NpcPVu70OpJIyFChS9jp1LgmczNSaVG3KqNeW8Rzf1+v4V4iqNAlTF1ZPZE3RybTv1NDnvv7Bh6etohDGu4lUU6FLmErMT6WZwa257c/aMNHq3fRf2weW/dquJdELxW6hDUz44HUpkwd3o2dXx4lLSuHXA33kiilQpeIkNqiNtkZKdSpmsC9kxfwUo6Ge0n0UaFLxLiqVmVmjU7httZ1ePLd1fzs7eUcPaHhXhI9VOgSUaokxDHuns786LYWzFxcxF0T89ml4V4SJVToEnFiYowf3daS8UM6s2HXQXq9kMPirfu8jiUSdCp0iVi3t72SWaN7kBgfy6AJ+bxVsM3rSCJBpUKXiNbqympkZ6TQpWlNfjFjOU/MXaXhXhKx/Cp0M7vdzNaZWaGZPXqe47qY2SkzGxC4iCKXp0alCrwyrCvDU5oyJXcz905ewL7DGu4lkeeChW5msUAW0BNoAww2szbnOO5/gQ8DHVLkcsXFxvBYrzY8M7ADBVv2kZaVw9qdX3odSySg/HmE3hUodM5tdM4dB6YDvcs4bgwwE9gdwHwiATWgc0PeHJnMsROn6Tc2j/dX7PA6kkjA+FPoDYDSzyYVlWz7DzNrAPQFxp/visxspJkVmFlBcXHxxWYVCYiOjWsyd0wqLetW5eFpi3n2Iw33ksjgT6FbGdvO/u1/Dvilc+687+Jwzk10zvmcc76kpCQ/I4oEXt1qiUwfmczAzg15/uMNPPSahntJ+POn0IuARqUuNwS2n3WMD5huZpuBAcBYM+sTiIAiwZIYH8sfB7Tnd73a8I+1u+mblcvmPYe9jiVyyfwp9IVACzNramYVgEFAdukDnHNNnXNNnHNNgBnAaOfc7ECHFQk0M+P+lKa8OrwrxYeOkZaZw6frdTpQwtMFC905dxLI4MyrV9YAbznnVpnZKDMbFeyAIuWhx9W1yU5PpX6Nitw/ZQGT/r1Rw70k7JhXv7Q+n88VFBR4ctsi53L42El+9vYy3l+5k74dG/B0v3Ykxsd6HUvkP8xskXPOV9Y+vVNUpJTKCXFk3d2Jn3ynJe8s+Zw7J8xjx4GvvI4l4hcVushZYmKMH97agolDO/PZ7kP0eiGXRVu+8DqWyAWp0EXO4bvXXsns9BQqJ8QyaGI+0xds9TqSyHmp0EXOo0XdqsxJTyG5WS0enbWCx+as5ISGe0mIUqGLXECNShWYcn8XRt7QjKnztjD0pfnsPXTM61gi36JCF/FDXGwMv/5+a/7vrg4s3rqftMxcVm0/4HUskW9QoYtchL4dGzJjVHdOnXYMGDePvy7XcC8JHSp0kYvUvmENssek0KZ+NdJfX8wzH67TcC8JCSp0kUtQp2oir4/oxqAujcj8ZyEjphbw5dETXseSKKdCF7lECXGxPN2vHU/2vpZ/rS+mb1YuG4sPeR1LopgKXeQymBlDuzfh1Qe6se/ICXpn5fLJOq3xIt5QoYsEQPfmtZiTnkLDmpUY/vJCJvzrMw33knKnQhcJkEZXVGLmw93p2bYeT7+/lh+9uZSjJ8675otIQKnQRQKoUoU4Mu/uyM+/dw3Zy7YzYHwe2/druJeUDxW6SICZGek3X82ke31s3nOEtMwcFm7WcC8JPhW6SJDc2rous9N7UDUxnrtfzOf1+RruJcGlQhcJoqvrVGV2ego9mtfm1++s4L9mr+D4SQ33kuBQoYsEWfWK8Uy+vwsP3diM1/K3MuSl+ezRcC8JAhW6SDmIjTF+1bM1fxl0Hcu27ad3Zi4rP9dwLwksFbpIOep9XQNmjOrBaecYMD6Pucu2ex1JIogKXaSctWtYneyMVNo1qM6YN5bwvx+s5ZSGe0kAqNBFPJBUNYFpDyYzuGtjxn3yGQ+8spADX2m4l1weFbqIRyrExfB0v3Y81actORv20Dcrl8LdGu4ll06FLuKxIclXMe3Bbhz46gR9s3L5x9pdXkeSMKVCFwkB3ZrVIntMKo1rVeKBVwoY+0mhhnvJRVOhi4SIBjUqMmNUD37Qvj5//GAdY95YwlfHNdxL/KdCFwkhFSvE8vyg63i0Zyv+umIH/cflUbTviNexJEyo0EVCjJkx6sbmTL6vC9v2HSEtM5f5G/d6HUvCgApdJETd3KoOs9NTqFEpnnsmzefV/C1eR5IQp0IXCWHNk6owOz2FG1om8dvZK/nVLA33knPzq9DN7HYzW2dmhWb2aBn77zGz5SUfeWbWIfBRRaJTtcR4XrzXx+ibmvPGgq3c/WI+xQc13Eu+7YKFbmaxQBbQE2gDDDazNmcdtgm40TnXHngSmBjooCLRLDbG+MXtrXhhcEdWbj9AWmYOK4o03Eu+yZ9H6F2BQufcRufccWA60Lv0Ac65POfcvpKL+UDDwMYUEYBeHeozY1QPYswYMD6POUs/9zqShBB/Cr0BsK3U5aKSbefyAPB+WTvMbKSZFZhZQXFxsf8pReQ/2jaoTnZGCh0a1eCR6Ut5+r01Gu4lgH+FbmVsK/O3x8xu5kyh/7Ks/c65ic45n3POl5SU5H9KEfmGWlUSmPZgN4YmX8WETzcy/OWFHDii4V7Rzp9CLwIalbrcEPjWEGczaw9MAno75/SiWZEgi4+N4ck+bfmfvu3I+2wPfcbmUrj7oNexxEP+FPpCoIWZNTWzCsAgILv0AWbWGJgFDHXOrQ98TBE5l7u7Neb1EckcPHqCPll5/H21hntFqwsWunPuJJABfAisAd5yzq0ys1FmNqrksMeAWsBYM1tqZgVBSywi39KlyRVkZ6TStHZlRrxaQOY/Nmi4VxQyr37oPp/PFRSo90UC6eiJU/xy5nLmLN3OHe3q8aeB7alUIc7rWBJAZrbIOecra5/eKSoSQRLjY3nuruv49fdb8f7KHfQbm8e2LzTcK1qo0EUijJkx8obmTL6/C5/v/4q0zBzyPtvjdSwpByp0kQh10zV1mJOewhWVKzD0pQW8krdZ59UjnApdJII1KxnudVPLJB7PXsWjM1dw7KQWzYhUKnSRCFe1ZLjXmFuu5s2CbQyemM/uL496HUuCQIUuEgViYoyffvcasu7uxJodB+mVmcOybfu9jiUBpkIXiSJ3tK/HzId7EB8bw8AJ85i5qMjrSBJAKnSRKNOmfjWyM1Lp3LgmP317GU+9u5qTp7RoRiRQoYtEoSsqV2DqA125v0cTJuVsYtjLC9l/5LjXseQyqdBFolR8bAy/S7uWP/ZvT/7GvfTOymX9Lg33CmcqdJEod2eXRkwf2Z0jx0/RNyuXD1ft9DqSXCIVuojQ+aqazM1I5eo6VXjo1UX85e8bOK1FM8KOCl1EALiyeiJvPtSdfh0b8H9/X8/oaYs5fOyk17HkIqjQReQ/EuNj+fOdHfivO1rzt9U76T8uj617NdwrXKjQReQbzIwHr2/GK8O7suPAUdKycsgt1HCvcKBCF5EyXd8iiTnpKSRVSeDeyQuYnLNJw71CnApdRM6pSe3KvJOewq2t6vD7d1fz8xnLOXpCw71ClQpdRM6rSkIc44d05pFbWzBjURGDJuazS8O9QpIKXUQuKCbG+PF3WjJ+SCfW7zpIrxdyWLJ1n9ex5CwqdBHx2+1t6zFrdA8S4mO4a0I+MzTcK6So0EXkorS6shrZ6an4mtTkZ28v44m5qzTcK0So0EXkotWsXIGpw7syLKUJU3I3c9+UBew7rOFeXlOhi8gliYuN4fFe1/KnAe1ZuGkfaVk5rN35pdexopoKXUQuy0BfI6Y/lMyxE6fpNzaPD1bu8DpS1FKhi8hl69S4JnPHpNKyblVGvbaY//tovYZ7eUCFLiIBUbdaItNHJtO/U0P+8vEGRr22iEMa7lWuVOgiEjCJ8bE8M7A9j/2gDR+v3U2/sbls2XvY61hRQ4UuIgFlZgxPbcrU4V3ZffAYaZm5/HtDsdexooIKXUSCIuXq2mSnp3JltUTum7yASf/eqOFeQaZCF5GgaVyrErNG9+A7bery1F/X8NO3l2m4VxCp0EUkqConxDHuns78+LaWzFr8OXdNzGfnAQ33Cga/Ct3MbjezdWZWaGaPlrHfzOz5kv3LzaxT4KOKSLiKiTEeua0FE4Z2pnDXQXpl5rBoi4Z7BdoFC93MYoEsoCfQBhhsZm3OOqwn0KLkYyQwLsA5RSQCfO/aK5k1OoWK8bEMnpjP5JxNbCw+pFkwARLnxzFdgULn3EYAM5sO9AZWlzqmNzDVnXnGI9/MaphZPeec3jImIt9wzZVVyc5IYcwbS/j9u6v5/btQITaGhjUrEhtjXscrF3d1acSD1zcL+PX6U+gNgG2lLhcB3fw4pgHwjUI3s5GceQRP48aNLzariESIGpUq8Mqwrqz4/AAbdh+icPchtu07EjWvgqldJSEo1+tPoZf1T+bZf+r+HINzbiIwEcDn80XHT05EyhQTY3RoVIMOjWp4HSVi+POkaBHQqNTlhsD2SzhGRESCyJ9CXwi0MLOmZlYBGARkn3VMNnBvyatdkoEDOn8uIlK+LnjKxTl30swygA+BWGCyc26VmY0q2T8eeA/4PlAIHAGGBS+yiIiUxZ9z6Djn3uNMaZfeNr7U1w5ID2w0ERG5GHqnqIhIhFChi4hECBW6iEiEUKGLiEQI8+qdWWZWDGy5xG+vDewJYJxwoPscHXSfo8Pl3OernHNJZe3wrNAvh5kVOOd8XucoT7rP0UH3OToE6z7rlIuISIRQoYuIRIhwLfSJXgfwgO5zdNB9jg5Buc9heQ5dRES+LVwfoYuIyFlU6CIiESLsCv1CC1ZHGjNrZGb/NLM1ZrbKzB7xOlN5MLNYM1tiZu96naW8lCzdOMPM1pb8vLt7nSmYzOzHJb/TK83sDTNL9DpTMJjZZDPbbWYrS227wsw+MrMNJZ9rBuK2wqrQ/VywOtKcBH7qnGsNJAPpUXCfAR4B1ngdopz9BfjAOdcK6EAE338zawD8EPA559pyZjT3IG9TBc3LwO1nbXsU+Ng51wL4uOTyZQurQqfUgtXOuePA1wtWRyzn3A7n3OKSrw9y5i95A29TBZeZNQTuACZ5naW8mFk14AbgJQDn3HHn3H5PQwVfHFDRzOKASkToKmfOuU+BL87a3Bt4peTrV4A+gbitcCv0cy1GHRXMrAnQEZjvcZRgew74BXDa4xzlqRlQDEwpOdU0ycwqex0qWJxznwPPAFs5s5j8Aefc37xNVa7qfr2qW8nnOoG40nArdL8Wo45EZlYFmAn8yDn3pdd5gsXMfgDsds4t8jpLOYsDOgHjnHMdgcME6L/hoajknHFvoClQH6hsZkO8TRX+wq3Qo3IxajOL50yZT3POzfI6T5ClAGlmtpkzp9RuMbPXvI1ULoqAIufc1//7msGZgo9UtwGbnHPFzrkTwCygh8eZytMuM6sHUPJ5dyCuNNwK3Z8FqyOKmRlnzquucc4963WeYHPO/co519A514QzP99/OOci/pGbc24nsM3MrinZdCuw2sNIwbYVSDazSiW/47cSwU8ClyEbuK/k6/uAOYG4Ur/WFA0V51qw2uNYwZYCDAVWmNnSkm2/LlnnVSLLGGBayYOVjUTwYuvOuflmNgNYzJlXci0hQkcAmNkbwE1AbTMrAh4H/gC8ZWYPcOYft4EBuS299V9EJDKE2ykXERE5BxW6iEiEUKGLiEQIFbqISIRQoYuIRAgVuohIhFChi4hEiP8HIWEHSqUBxrgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_func(partial(tri, i=8))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = data.clone()\n", "x = data[0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([26.204, 26.349])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(torch.Size([2]), torch.Size([1500, 2]), torch.Size([1, 2]))" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.shape,X.shape,x[None].shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[ 0.000, 0.000],\n", " [ 0.513, -3.865],\n", " [-4.227, -2.345],\n", " [ 0.557, -3.685],\n", " [-5.033, -3.745],\n", " [-4.073, -0.638],\n", " [-3.415, -5.601],\n", " [-1.920, -5.686]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(x[None]-X)[:8]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[ 0.000, 0.000],\n", " [ 0.513, -3.865],\n", " [-4.227, -2.345],\n", " [ 0.557, -3.685],\n", " [-5.033, -3.745],\n", " [-4.073, -0.638],\n", " [-3.415, -5.601],\n", " [-1.920, -5.686]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(x-X)[:8]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([0.000, 3.899, 4.834, 3.726, 6.273, 4.122, 6.560, 6.002])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# rewrite using torch.einsum\n", "dist = ((x-X)**2).sum(1).sqrt()\n", "dist[:8]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([ 0.160, 0.047, 0.025, ..., 0.000, 0.000, 0.000])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight = gaussian(dist, 2.5)\n", "weight" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(torch.Size([1500]), torch.Size([1500, 2]))" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight.shape,X.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "torch.Size([1500, 1])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight[:,None].shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[ 4.182, 4.205],\n", " [ 1.215, 1.429],\n", " [ 0.749, 0.706],\n", " ...,\n", " [ 0.000, 0.000],\n", " [ 0.000, 0.000],\n", " [ 0.000, 0.000]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight[:,None]*X" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def one_update(X):\n", " for i, x in enumerate(X):\n", " dist = torch.sqrt(((x-X)**2).sum(1))\n", "# weight = gaussian(dist, 2.5)\n", " weight = tri(dist, 8)\n", " X[i] = (weight[:,None]*X).sum(0)/weight.sum()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def meanshift(data):\n", " X = data.clone()\n", " for it in range(5): one_update(X)\n", " return X" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 453 ms, sys: 0 ns, total: 453 ms\n", "Wall time: 452 ms\n" ] } ], "source": [ "%time X=meanshift(data)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUzElEQVR4nO3dUWyk13mf8eeNpDhZpC1tehNrJG0pIGJq2WEciBBctBekrcBKtaiYAELli2KrGCQCuMI68IXl6GJmCggwUDRboHHQkrA3ujCiCki8EnaTprJKwi2QRKECZy1lvaNFTNkLCtJqEjYNCCiQ/faCM9TMktw1OTMc8vD5AQPOd87MnHOw5H/PnDnfN5GZSJLK9GPD7oAkaXAMeUkqmCEvSQUz5CWpYIa8JBXs1mF3oNMHP/jBHBsbG3Y3JOlQefnll9/OzOPb1R2okB8bG2N5eXnY3ZCkQyUiXt+pzuUaSSqYIV+AZrO5q3JJR4chf8jVajUmJiZoNBpd5Y1Gg4mJCWq12nA6JulAMOQPsVqtRr1eZ3V1lenp6c2gbzQaTE9Ps7q6Sr1eN+ilI8yQP6TaAd/WDvoLFy5sBnybQS8dXYb8IdRsNllYWOgqm2GG9dV1Tp48yfrqOjPMdNUvLCy4Ri8dQYb8ITQ6Osri4iKVSgXYCPjTnOYMZxhjjDOc4TSnN4O+UqmwuLjI6OjoEHstaRgM+UNqfHx8M+iXWGKFFcYY4yxnGWOMFVZYYmkz4MfHx4fdZUlDYMgfYuPj48zPz7PGGnXqXXV16qyxxvz8vAEvHWGG/CHWaDSYm5tjhBGqVLvqqlQZYYS5ubkt2yslHR2G/CHVuU1yiqnNJZrHeGxz6WaKqS3bKyUdLXGQvv5vcnIyvXbNzTWbTSYmJrq2Sc4wwxJLrLHGCCNMMcU5zm3WVyoVLl686IevUoEi4uXMnNyuzpn8ITQ6Osrs7GxX2TnOcaxyjPPnz3Oscqwr4AFmZ2cNeGkXSrlciCF/SNVqNarV99bh27toHnrooa7tlQDVatWToaRdKOpyIZl5YG733XdfqkP1H7932+kh1WpWKpW8fPlyV/nly5ezUqlktVodcCelslSr1QQS6Prbav9NtesO0t8WsJw75Kpr8gdZ7Z903P+/Oz6s2WxuuxSzU7mk7V1/uRDYeJc8Pz/P3Nxc1+dgcHDeJbsmX7idgtyAl350pV4u5EB9M5Suc4PZu6T+al8upL01uX25kId5mDp1qlQZYwzY2OhwWC4X4kxeklpKvFxIzyEfET8RES9FxF9GxKsRUW+VfyAiXoiI11o/3997dyVpsEq7XEg/ZvLvAJ/IzF8APgY8GBEfB54AXszMe4AXW8eSdKCVdrmQnkO+tYPn71uHt7VuCTwMPN0qfxqu+8RCkg6YEi8X0pctlBFxC/Ay8LPAlzPzCxGxlpkjHY/528zcsmQTEXPAHMCJEyfue/3113vujyTt1mG+XMjAt1Bm5g8y82PAncD9EfHRXTx3PjMnM3Py+PHj/eiOJO1aqZcL6esWysxci4gl4EHgzYi4PTPfiIjbgbf62ZYk7cXYExc276986aGuuvaJTe0Tojp30XRur4SDcyLUzfRjd83xiBhp3f9J4AHgO8DzwKnWw04Bz/XaliQNWvu6UNdvk+zcXnlYAh76M5O/HXi6tS7/Y8CzmXk+Iv4EeDYiPgN8D3ikD21J0sDVajUef/zxLUsx4+PjB2INfje8do0kHXJeu0aSjihDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKtiRCvlms7mrckk67I5MyNdqNSYmJrZ8J2Oj0WBiYuLQXBtaknbjSIR8rVajXq9v+fLdzi/trdfrBr2k4hQf8u2Ab2sH/YULF7q+ygsw6CUVp+iQbzabLCwsdJXNMMP66jonT55kfXWdGWa66hcWFlyjl1SMokN+dHR08zsZYSPgT3OaM5xhjDHOcIbTnN4M+vZ3Oh6mr/aSpBspOuSh+8t3l1hihRXGGOMsZxljjBVWWGJpy5f2SlIJig952Aj6+fl51lijTr2rrk6dNdaYn5834CUV50iEfKPRYG5ujhFGqFLtqqtSZYQR5ubmtmyvlKTDrviQ79wmOcXU5hLNYzy2uXQzxdSW7ZWSVILIzGH3YdPk5GQuLy/37fWazSYTExNd2yRnmGGJJdZYY4QRppjiHOc26yuVChcvXvTDV0mHRkS8nJmT29UVPZMfHR1ldna2q+wc5zhWOcb58+c5VjnWFfAAs7OzBrykYvQc8hFxV0QsRsSliHg1Ik63yj8QES9ExGutn+/vvbu7V6vVqFbfW4dv76J56KGHurZXAlSrVU+GklSUW/vwGu8Cn8/Mv4iIfwS8HBEvAP8OeDEzvxQRTwBPAF/oQ3vbuvTPPrx5/8PfudRV1w7uhYWFrm2S7e2V09PTzM7OGvCSitP3NfmIeA747dZtKjPfiIjbgaXM/LkbPbeXNfkbhXxbs9ncdilmp3JJOgz2bU0+IsaAXwT+DPiZzHwDoPXzp3d4zlxELEfE8rVr1/rZnS12CnIDXlKp+rFcA0BE/BTw+8DnMvPvIuJHel5mzgPzsDGT32v7O83eJeko68tMPiJuYyPgv5aZf9AqfrO1TEPr51v9aEuS9KPrx+6aAL4CXMrM3+qoeh441bp/Cniu17YkSbvTj+WafwH8W+DbEfGtVtlvAl8Cno2IzwDfAx7pQ1uSpF3oOeQz8/8AOy3Af7LX15ck7V3RZ7xK0lFnyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8gVqNpu7KpdUrr6EfER8NSLeiohXOso+EBEvRMRrrZ/v70dburFarcbExASNRqOrvNFoMDExQa1WG07HJA1Fv2byvws8eF3ZE8CLmXkP8GLrWANUq9Wo1+usrq4yPT29GfSNRoPp6WlWV1ep1+sGvTRE+/1Ouy8hn5nfBP7muuKHgadb958GZvrRlrbXDvi2dtBfuHBhM+DbDHppOIbyTjsz+3IDxoBXOo7Xrqv/2x2eNwcsA8snTpxI7d7bb7+dlUolgc3bDDM5wkgCOcJIzjDTVV+pVPLtt98edtelI6NarXb9/V2+fDkzMy9fvtz191utVnf92sBy7pTNO1Xs9rbXkO+83XfffbsenDZ0/qLMMJOLLOZZzuYYY3mWs7nI4mbQd/6CSRq8zoDvDPrz589vmaDtJehvFPKD3F3zZkTcDtD6+dYA2zryxsfHWVxcpFKpsMQSK6wwxhhnOcsYY6ywwhJLVCoVFhcXGR8fH3aXpSOh2WyysLDQVTbDDOur65w8eZL11XVmrlvNXlhY6Nsa/SBD/nngVOv+KeC5AbYlNoJ+fn6eNdaoU++qq1NnjTXm5+cNeGkfjY6Obk7AYCPgT3OaM5xhjDHOcIbTnN4M+vZEbHR0tC/t92sL5e8BfwL8XERcjYjPAF8CfikiXgN+qXWsAWo0GszNzTHCCFWqXXVVqowwwtzc3JYPfSQN1jDfafdrd82nM/P2zLwtM+/MzK9kZjMzP5mZ97R+Xr/7Rn3UuU1yiqnNX5zHeGzzF2qKqS3bKyXtj2G9046NNfuDYXJyMpeXl4fdjUOn2WwyMTHRtU1yhhmWWGKNNUYYYYopznFus75SqXDx4sW+vSWUdGPtidj66vrmUk3bCiv8Br/BscqxPc3kI+LlzJzcrs7LGhRgdHSU2dnZrrJznONY5Rjnz5/nWOVYV8ADzM7OGvDSPhnqO+2dtt0M4+YWyt4Mch+upL3Zj/NYGNIWSvXRzz/985u3ndRqNarV6pYPbzo/9KlWq57tKu2jYb/Tdk3+kOgM92+f+vYNH9tsNrf9BdmpXNLgdV56pHMi1rmUA+xpInajNflbe+u2DqKdgtyAlwbnP/2bk5v3P//fz2+pbwf3wsLCtu+0p6enmZ2d7fs7bWfyktQHNwv5tkG803Z3jSQdEPv9TtvlGknqgxvN3ofJmbwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFG3jIR8SDEXE5Iq5ExBODbk+S9J6BhnxE3AJ8Gfhl4F7g0xFx7yDblCS9Z9Az+fuBK5n515n5D8AzwMMDblOS1DLokL8D+H7H8dVW2aaImIuI5YhYvnbt2oC7I0lHy6BDPrYpy66DzPnMnMzMyePHjw+4O5J0tAw65K8Cd3Uc3wmsDrhNSVLLoEP+z4F7IuLuiPhx4FHg+QG3KUlquXWQL56Z70bEvwf+GLgF+GpmvjrINiVJ7xloyANk5h8CfzjodiRJW3nGqyQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWA9hXxEPBIRr0bEDyNi8rq6L0bElYi4HBGf6q2bkqS9uLXH578C/Crw3zoLI+Je4FHgI0AF+EZEjGfmD3psT5K0Cz3N5DPzUmZe3qbqYeCZzHwnM78LXAHu76UtSdLuDWpN/g7g+x3HV1tlkqR9dNPlmoj4BvChbaqezMzndnraNmW5w+vPAXMAJ06cuFl3JEm7cNOQz8wH9vC6V4G7Oo7vBFZ3eP15YB5gcnJy2/8IJEl7M6jlmueBRyPifRFxN3AP8NKA2pIk7aDXLZS/EhFXgX8OXIiIPwbIzFeBZ4G/Av4H8Fl31kjS/utpC2Vmfh34+g51TwFP9fL6kqTeeMarJBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBWsp5CPiP8YEd+JiIsR8fWIGOmo+2JEXImIyxHxqZ57KknatV5n8i8AH83MCaABfBEgIu4FHgU+AjwI/E5E3NJjW5KkXeop5DPzf2bmu63DPwXubN1/GHgmM9/JzO8CV4D7e2lLkrR7/VyT/zXgj1r37wC+31F3tVW2RUTMRcRyRCxfu3atj92RJN16swdExDeAD21T9WRmPtd6zJPAu8DX2k/b5vG53etn5jwwDzA5ObntYyRJe3PTkM/MB25UHxGngJPAJzOzHdJXgbs6HnYnsLrXTkqS9qbX3TUPAl8A/nVmrndUPQ88GhHvi4i7gXuAl3ppS5K0ezedyd/EbwPvA16ICIA/zcxfz8xXI+JZ4K/YWMb5bGb+oMe2JEm71FPIZ+bP3qDuKeCpXl5fkkrWbDYZHR39kcv3wjNeJWkIarUaExMTNBqNrvJGo8HExAS1Wq0v7RjykrTParUa9Xqd1dVVpqenN4O+0WgwPT3N6uoq9Xq9L0FvyEvSPmoHfFs76C9cuLAZ8G39CHpDXpL2SbPZZGFhoatshhnWV9c5efIk66vrzDDTVb+wsECz2dxzm4a8JO2T0dFRFhcXqVQqwEbAn+Y0ZzjDGGOc4QynOb0Z9JVKhcXFxZ4+hDXkJWkfjY+Pbwb9EkussMIYY5zlLGOMscIKSyxtBvz4+HhP7RnykrTPxsfHmZ+fZ4016tS76urUWWON+fn5ngMeDHlJ2neNRoO5uTlGGKFKtauuSpURRpibm9uyvXIvDHlJ2ked2ySnmNpconmMxzaXbqaY2rK9cq/ivWuKDd/k5GQuLy8PuxuSNBDNZpOJiYmubZIzzLDEEmusMcIIU0xxjnOb9ZVKhYsXL97ww9eIeDkzJ7ercyYvSftkdHSU2dnZrrJznONY5Rjnz5/nWOVYV8ADzM7O9rS7ptcLlEmSrvPlX/9fm/c/+18/0VXXPrmpfUJU5y6axcXFrhOiqtVqzydDGfKStM/awb2wsNC1TbIz6GdnZ/tyWQNDXpKGoFar8fjjj29ZihkfH7/pGvxuGPKS1GfXL9HsZKcg71fAgx+8SlLRDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUsAN1gbKIuAa8Pux+7IMPAm8PuxP7zDEfDY55OP5pZh7fruJAhfxRERHLO10xrlSO+WhwzAePyzWSVDBDXpIKZsgPx/ywOzAEjvlocMwHjGvyklQwZ/KSVDBDXpIKZsjvo4h4JCJejYgfRsTkdXVfjIgrEXE5Ij41rD4OQkQ82BrXlYh4Ytj9GYSI+GpEvBURr3SUfSAiXoiI11o/3z/MPvZbRNwVEYsRcan1e326VV7suCPiJyLipYj4y9aY663yAztmQ35/vQL8KvDNzsKIuBd4FPgI8CDwOxFxy/53r/9a4/gy8MvAvcCnW+Mtze+y8W/X6Qngxcy8B3ixdVySd4HPZ+aHgY8Dn23925Y87neAT2TmLwAfAx6MiI9zgMdsyO+jzLyUmZe3qXoYeCYz38nM7wJXgPv3t3cDcz9wJTP/OjP/AXiGjfEWJTO/CfzNdcUPA0+37j8NzOxnnwYtM9/IzL9o3f9/wCXgDgoed274+9bhba1bcoDHbMgfDHcA3+84vtoqK0HJY7uZn8nMN2AjEIGfHnJ/BiYixoBfBP6MwscdEbdExLeAt4AXMvNAj9nveO2ziPgG8KFtqp7MzOd2eto2ZaXsbS15bAIi4qeA3wc+l5l/F7HdP3k5MvMHwMciYgT4ekR8dMhduiFDvs8y84E9PO0qcFfH8Z3Aan96NHQlj+1m3oyI2zPzjYi4nY2ZX1Ei4jY2Av5rmfkHreLixw2QmWsRscTGZzEHdswu1xwMzwOPRsT7IuJu4B7gpSH3qV/+HLgnIu6OiB9n4wPm54fcp/3yPHCqdf8UsNM7uUMpNqbsXwEuZeZvdVQVO+6ION6awRMRPwk8AHyHAzxmz3jdRxHxK8B/AY4Da8C3MvNTrbongV9jY8fC5zLzj4bVz36LiH8F/GfgFuCrmfnUcHvUfxHxe8AUG5edfROoAueAZ4ETwPeARzLz+g9nD62I+JfA/wa+DfywVfybbKzLFznuiJhg44PVW9iYJD+bmf8hIkY5oGM25CWpYC7XSFLBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUsP8Pfzg4GNqH4TYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_data(centroids+2, X, n_samples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Animation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def do_one(d):\n", " if d: one_update(X)\n", " ax.clear()\n", " plot_data(centroids+2, X, n_samples, ax=ax)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# create your own animation\n", "X = data.clone()\n", "fig,ax = plt.subplots()\n", "ani = FuncAnimation(fig, do_one, frames=5, interval=500, repeat=False)\n", "plt.close()\n", "HTML(ani.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## GPU batched algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To truly accelerate the algorithm, we need to be performing updates on a batch of points per iteration, instead of just one as we were doing." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(torch.Size([5, 2]), torch.Size([1500, 2]))" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bs=5\n", "X = data.clone()\n", "x = X[:bs]\n", "x.shape,X.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def dist_b(a,b): return (((a[None]-b[:,None])**2).sum(2)).sqrt()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[ 0.000, 3.899, 4.834, ..., 17.628, 22.610, 21.617],\n", " [ 3.899, 0.000, 4.978, ..., 21.499, 26.508, 25.500],\n", " [ 4.834, 4.978, 0.000, ..., 19.373, 24.757, 23.396],\n", " [ 3.726, 0.185, 4.969, ..., 21.335, 26.336, 25.333],\n", " [ 6.273, 5.547, 1.615, ..., 20.775, 26.201, 24.785]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dist_b(X, x)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "torch.Size([5, 1500])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dist_b(X, x).shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(torch.Size([1, 1500, 2]), torch.Size([5, 1, 2]), torch.Size([5, 1500, 2]))" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X[None,:].shape, x[:,None].shape, (X[None,:]-x[:,None]).shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[ 0.199, 0.030, 0.011, ..., 0.000, 0.000, 0.000],\n", " [ 0.030, 0.199, 0.009, ..., 0.000, 0.000, 0.000],\n", " [ 0.011, 0.009, 0.199, ..., 0.000, 0.000, 0.000],\n", " [ 0.035, 0.199, 0.009, ..., 0.000, 0.000, 0.000],\n", " [ 0.001, 0.004, 0.144, ..., 0.000, 0.000, 0.000]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight = gaussian(dist_b(X, x), 2)\n", "weight" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(torch.Size([5, 1500]), torch.Size([1500, 2]))" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight.shape,X.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(torch.Size([5, 1500, 1]), torch.Size([1, 1500, 2]))" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight[...,None].shape, X[None].shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "torch.Size([5, 2])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "num = (weight[...,None]*X[None]).sum(1)\n", "num.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[367.870, 386.231],\n", " [518.332, 588.680],\n", " [329.665, 330.782],\n", " [527.617, 598.217],\n", " [231.302, 234.155]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "num" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[367.870, 386.231],\n", " [518.332, 588.680],\n", " [329.665, 330.782],\n", " [527.617, 598.218],\n", " [231.302, 234.155]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "torch.einsum('ij,jk->ik', weight, X)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[367.870, 386.231],\n", " [518.332, 588.680],\n", " [329.665, 330.782],\n", " [527.617, 598.218],\n", " [231.302, 234.155]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "weight@X" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "torch.Size([5, 1])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "div = weight.sum(1, keepdim=True)\n", "div.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "tensor([[26.376, 27.692],\n", " [26.101, 29.643],\n", " [28.892, 28.990],\n", " [26.071, 29.559],\n", " [29.323, 29.685]])" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "num/div" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def meanshift(data, bs=500):\n", " n = len(data)\n", " X = data.clone()\n", " for it in range(5):\n", " for i in range(0, n, bs):\n", " s = slice(i, min(i+bs,n))\n", " weight = gaussian(dist_b(X, X[s]), 2.5)\n", "# weight = tri(dist_b(X, X[s]), 8)\n", " div = weight.sum(1, keepdim=True)\n", " X[s] = weight@X/div\n", " return X" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although each iteration still has to launch a new cuda kernel, there are now fewer iterations, and the acceleration from updating a batch of points more than makes up for it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = data.cuda()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = meanshift(data).cpu()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2 ms ± 226 µs per loop (mean ± std. dev. of 7 runs, 5 loops each)\n" ] } ], "source": [ "%timeit -n 5 _=meanshift(data, 1250).cpu()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUzElEQVR4nO3dUWyk13mf8eeNpDhZpC1tehNrJG0pIGJq2WEciBBctBekrcBKtaiYAELli2KrGCQCuMI68IXl6GJmCggwUDRboHHQkrA3ujCiCki8EnaTprJKwi2QRKECZy1lvaNFTNkLCtJqEjYNCCiQ/faCM9TMktw1OTMc8vD5AQPOd87MnHOw5H/PnDnfN5GZSJLK9GPD7oAkaXAMeUkqmCEvSQUz5CWpYIa8JBXs1mF3oNMHP/jBHBsbG3Y3JOlQefnll9/OzOPb1R2okB8bG2N5eXnY3ZCkQyUiXt+pzuUaSSqYIV+AZrO5q3JJR4chf8jVajUmJiZoNBpd5Y1Gg4mJCWq12nA6JulAMOQPsVqtRr1eZ3V1lenp6c2gbzQaTE9Ps7q6Sr1eN+ilI8yQP6TaAd/WDvoLFy5sBnybQS8dXYb8IdRsNllYWOgqm2GG9dV1Tp48yfrqOjPMdNUvLCy4Ri8dQYb8ITQ6Osri4iKVSgXYCPjTnOYMZxhjjDOc4TSnN4O+UqmwuLjI6OjoEHstaRgM+UNqfHx8M+iXWGKFFcYY4yxnGWOMFVZYYmkz4MfHx4fdZUlDYMgfYuPj48zPz7PGGnXqXXV16qyxxvz8vAEvHWGG/CHWaDSYm5tjhBGqVLvqqlQZYYS5ubkt2yslHR2G/CHVuU1yiqnNJZrHeGxz6WaKqS3bKyUdLXGQvv5vcnIyvXbNzTWbTSYmJrq2Sc4wwxJLrLHGCCNMMcU5zm3WVyoVLl686IevUoEi4uXMnNyuzpn8ITQ6Osrs7GxX2TnOcaxyjPPnz3Oscqwr4AFmZ2cNeGkXSrlciCF/SNVqNarV99bh27toHnrooa7tlQDVatWToaRdKOpyIZl5YG733Xdf6jrVf/zebbvqajUrlUpevny5q/zy5ctZqVSyWq3uQyelclSr1QQS6Prbav9NtesO0t8WsJw75Kpr8gdd7Z903P+/2z6k2WxuuxSzU7mk7V1/uRDYeJc8Pz/P3Nxc1+dgcHDeJbsmX7idgtyAl350pV4u5EB9M5S2scPsXVJ/tS8X0t6a3L5cyMM8TJ06VaqMMQZsbHQ4LJcLcSYvSS0lXi6k55CPiJ+IiJci4i8j4tWIqLfKPxARL0TEa62f7++9u5I0WKVdLqQfM/l3gE9k5i8AHwMejIiPA08AL2bmPcCLrWNJOtBKu1xIzyHf2sHz963D21q3BB4Gnm6VPw3XfWIhSQdMiZcL6csWyoi4BXgZ+Fngy5n5hYhYy8yRjsf8bWZuWbKJiDlgDuDEiRP3vf766z33R5J26zBfLmTgWygz8weZ+THgTuD+iPjoLp47n5mTmTl5/PjxfnRHknat1MuF9HULZWauRcQS8CDwZkTcnplvRMTtwFv9bEuS9mLsiQub91e+9FBXXfvEpvYJUZ27aDq3V8LBORHqZvqxu+Z4RIy07v8k8ADwHeB54FTrYaeA53ptS5IGrX1dqOu3SXZurzwsAQ/9mcnfDjzdWpf/MeDZzDwfEX8CPBsRnwG+BzzSh7YkaeBqtRqPP/74lqWY8fHxA7EGvxteu0aSDjmvXSNJR5QhL0kFM+QlqWCGvCQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFexIhXyz2dxVuSQddkcm5Gu1GhMTE1u+k7HRaDAxMXForg0tSbtxJEK+VqtRr9e3fPlu55f21ut1g15ScYoP+XbAt7WD/sKFC11f5QUY9JKKU3TIN5tNFhYWuspmmGF9dZ2TJ0+yvrrODDNd9QsLC67RSypG0SE/Ojq6+Z2MsBHwpznNGc4wxhhnOMNpTm8Gffs7HQ/TV3tJ0o0UHfLQ/eW7SyyxwgpjjHGWs4wxxgorLLG05Ut7JakExYc8bAT9/Pw8a6xRp95VV6fOGmvMz88b8JKKcyRCvtFoMDc3xwgjVKl21VWpMsIIc3NzW7ZXStJhV3zId26TnGJqc4nmMR7bXLqZYmrL9kpJKkFk5rD7sGlycjKXl5f79nrNZpOJiYmubZIzzLDEEmusMcIIU0xxjnOb9ZVKhYsXL/rhq6RDIyJezszJ7eqKnsmPjo4yOzvbVXaOcxyrHOP8+fMcqxzrCniA2dlZA15SMXoO+Yi4KyIWI+JSRLwaEadb5R+IiBci4rXWz/f33t3dq9VqVKvvrcO3d9E89NBDXdsrAarVqidDSSrKrX14jXeBz2fmX0TEPwJejogXgH8HvJiZX4qIJ4AngC/0ob1tXfpnH968/+HvXOqqawf3wsJC1zbJ9vbK6elpZmdnDXhJxen7mnxEPAf8dus2lZlvRMTtwFJm/tyNntvLmvyNQr6t2WxuuxSzU7kkHQb7tiYfEWPALwJ/BvxMZr4B0Pr50zs8Zy4iliNi+dq1a/3szhY7BbkBL6lU/ViuASAifgr4feBzmfl3EfEjPS8z54F52JjJ77X9nWbvknSU9WUmHxG3sRHwX8vMP2gVv9lapqH1861+tCVJ+tH1Y3dNAF8BLmXmb3VUPQ+cat0/BTzXa1uSpN3px3LNvwD+LfDtiPhWq+w3gS8Bz0bEZ4DvAY/0oS1J0i70HPKZ+X+AnRbgP9nr60uS9q7oM14l6agz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEM+QI1m81dlUsqV19CPiK+GhFvRcQrHWUfiIgXIuK11s/396Mt3VitVmNiYoJGo9FV3mg0mJiYoFarDadjkoaiXzP53wUevK7sCeDFzLwHeLF1rAGq1WrU63VWV1eZnp7eDPpGo8H09DSrq6vU63WDXhqi/X6n3ZeQz8xvAn9zXfHDwNOt+08DM/1oS9trB3xbO+gvXLiwGfBtBr00HEN5p52ZfbkBY8ArHcdr19X/7Q7PmwOWgeUTJ06kdu/tt9/OSqWSwOZthpkcYSSBHGEkZ5jpqq9UKvn2228Pu+vSkVGtVrv+/i5fvpyZmZcvX+76+61Wq7t+bWA5d8rmnSp2e9tryHfe7rvvvl0PThs6f1FmmMlFFvMsZ3OMsTzL2VxkcTPoO3/BJA1eZ8B3Bv358+e3TND2EvQ3CvlB7q55MyJuB2j9fGuAbR154+PjLC4uUqlUWGKJFVYYY4yznGWMMVZYYYklKpUKi4uLjI+PD7vL0pHQbDZZWFjoKpthhvXVdU6ePMn66joz161mLyws9G2NfpAh/zxwqnX/FPDcANsSG0E/Pz/PGmvUqXfV1amzxhrz8/MGvLSPRkdHNydgsBHwpznNGc4wxhhnOMNpTm8GfXsiNjo62pf2+7WF8veAPwF+LiKuRsRngC8BvxQRrwG/1DrWADUaDebm5hhhhCrVrroqVUYYYW5ubsuHPpIGa5jvtPu1u+bTmXl7Zt6WmXdm5lcys5mZn8zMe1o/r999oz7q3CY5xdTmL85jPLb5CzXF1JbtlZL2x7DeacfGmv3BMDk5mcvLy8PuxqHTbDaZmJjo2iY5wwxLLLHGGiOMMMUU5zi3WV+pVLh48WLf3hJKurH2RGx9dX1zqaZthRV+g9/gWOXYnmbyEfFyZk5uV+dlDQowOjrK7OxsV9k5znGscozz589zrHKsK+ABZmdnDXhpnwz1nfZO226GcXMLZW8GuQ9X0t7sx3ksDGkLpfro55/++c3bTmq1GtVqdcuHN50f+lSrVc92lfbRsN9puyZ/SHSG+7dPffuGj202m9v+guxULmnwOi890jkR61zKAfY0EbvRmvytvXVbB9FOQW7AS4Pzn/7Nyc37n//v57fUt4N7YWFh23fa09PTzM7O9v2dtjN5SeqDm4V82yDeabu7RpIOiP1+p+1yjST1wY1m78PkTF6SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCDTzkI+LBiLgcEVci4olBtydJes9AQz4ibgG+DPwycC/w6Yi4d5BtSpLeM+iZ/P3Alcz868z8B+AZ4OEBtylJahl0yN8BfL/j+GqrbFNEzEXEckQsX7t2bcDdkaSjZdAhH9uUZddB5nxmTmbm5PHjxwfcHUk6WgYd8leBuzqO7wRWB9ymJKll0CH/58A9EXF3RPw48Cjw/IDblCS13DrIF8/MdyPi3wN/DNwCfDUzXx1km5Kk9ww05AEy8w+BPxx0O5KkrTzjVZIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVLCeQj4iHomIVyPihxExeV3dFyPiSkRcjohP9dZNSdJe3Nrj818BfhX4b52FEXEv8CjwEaACfCMixjPzBz22J0nahZ5m8pl5KTMvb1P1MPBMZr6Tmd8FrgD399KWJGn3BrUmfwfw/Y7jq60ySdI+uulyTUR8A/jQNlVPZuZzOz1tm7Lc4fXngDmAEydO3Kw7kqRduGnIZ+YDe3jdq8BdHcd3Aqs7vP48MA8wOTm57X8EkqS9GdRyzfPAoxHxvoi4G7gHeGlAbUmSdtDrFspfiYirwD8HLkTEHwNk5qvAs8BfAf8D+Kw7ayRp//W0hTIzvw58fYe6p4Cnenl9SVJvPONVkgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpIIZ8pJUMENekgrWU8hHxH+MiO9ExMWI+HpEjHTUfTEirkTE5Yj4VM89lSTtWq8z+ReAj2bmBNAAvggQEfcCjwIfAR4EficibumxLUnSLvUU8pn5PzPz3dbhnwJ3tu4/DDyTme9k5neBK8D9vbQlSdq9fq7J/xrwR637dwDf76i72irbIiLmImI5IpavXbvWx+5Ikm692QMi4hvAh7apejIzn2s95kngXeBr7adt8/jc7vUzcx6YB5icnNz2MZKkvblpyGfmAzeqj4hTwEngk5nZDumrwF0dD7sTWN1rJyVJe9Pr7poHgS8A/zoz1zuqngcejYj3RcTdwD3AS720JUnavZvO5G/it4H3AS9EBMCfZuavZ+arEfEs8FdsLON8NjN/0GNbkqRd6inkM/Nnb1D3FPBUL68vSSVrNpuMjo7+yOV74RmvkjQEtVqNiYkJGo1GV3mj0WBiYoJardaXdgx5SdpntVqNer3O6uoq09PTm0HfaDSYnp5mdXWVer3el6A35CVpH7UDvq0d9BcuXNgM+LZ+BL0hL0n7pNlssrCw0FU2wwzrq+ucPHmS9dV1Zpjpql9YWKDZbO65TUNekvbJ6Ogoi4uLVCoVYCPgT3OaM5xhjDHOcIbTnN4M+kqlwuLiYk8fwhrykrSPxsfHN4N+iSVWWGGMMc5yljHGWGGFJZY2A358fLyn9gx5Sdpn4+PjzM/Ps8YadepddXXqrLHG/Px8zwEPhrwk7btGo8Hc3BwjjFCl2lVXpcoII8zNzW3ZXrkXhrwk7aPObZJTTG0u0TzGY5tLN1NMbdleuVfx3jXFhm9ycjKXl5eH3Q1JGohms8nExETXNskZZlhiiTXWGGGEKaY4x7nN+kqlwsWLF2/44WtEvJyZk9vVOZOXpH0yOjrK7OxsV9k5znGscozz589zrHKsK+ABZmdne9pd0+sFyiRJ1/nyr/+vzfuf/a+f6Kprn9zUPiGqcxfN4uJi1wlR1Wq155OhDHlJ2mft4F5YWOjaJtkZ9LOzs325rIEhL0lDUKvVePzxx7csxYyPj990DX43DHlJ6rPrl2h2slOQ9yvgwQ9eJalohrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkq2IG6QFlEXANeH3Y/9sEHgbeH3Yl95piPBsc8HP80M49vV3GgQv6oiIjlna4YVyrHfDQ45oPH5RpJKpghL0kFM+SHY37YHRgCx3w0OOYDxjV5SSqYM3lJKpghL0kFM+T3UUQ8EhGvRsQPI2LyurovRsSViLgcEZ8aVh8HISIebI3rSkQ8Mez+DEJEfDUi3oqIVzrKPhARL0TEa62f7x9mH/stIu6KiMWIuNT6vT7dKi923BHxExHxUkT8ZWvM9Vb5gR2zIb+/XgF+FfhmZ2FE3As8CnwEeBD4nYi4Zf+713+tcXwZ+GXgXuDTrfGW5nfZ+Lfr9ATwYmbeA7zYOi7Ju8DnM/PDwMeBz7b+bUse9zvAJzLzF4CPAQ9GxMc5wGM25PdRZl7KzMvbVD0MPJOZ72Tmd4ErwP3727uBuR+4kpl/nZn/ADzDxniLkpnfBP7muuKHgadb958GZvazT4OWmW9k5l+07v8/4BJwBwWPOzf8fevwttYtOcBjNuQPhjuA73ccX22VlaDksd3Mz2TmG7ARiMBPD7k/AxMRY8AvAn9G4eOOiFsi4lvAW8ALmXmgx+x3vPZZRHwD+NA2VU9m5nM7PW2bslL2tpY8NgER8VPA7wOfy8y/i9jun7wcmfkD4GMRMQJ8PSI+OuQu3ZAh32eZ+cAennYVuKvj+E5gtT89GrqSx3Yzb0bE7Zn5RkTczsbMrygRcRsbAf+1zPyDVnHx4wbIzLWIWGLjs5gDO2aXaw6G54FHI+J9EXE3cA/w0pD71C9/DtwTEXdHxI+z8QHz80Pu0355HjjVun8K2Omd3KEUG1P2rwCXMvO3OqqKHXdEHG/N4ImInwQeAL7DAR6zZ7zuo4j4FeC/AMeBNeBbmfmpVt2TwK+xsWPhc5n5R8PqZ79FxL8C/jNwC/DVzHxquD3qv4j4PWCKjcvOvglUgXPAs8AJ4HvAI5l5/Yezh1ZE/EvgfwPfBn7YKv5NNtblixx3REyw8cHqLWxMkp/NzP8QEaMc0DEb8pJUMJdrJKlghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkq2P8HgsA4GJpUh/wAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_data(centroids+2, X, n_samples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Homework:** implement k-means clustering, dbscan, locality sensitive hashing, or some other clustering, fast nearest neighbors, or similar algorithm of your choice, on the GPU. Check if your version is faster than a pure python or CPU version.\n", "\n", "Bonus: Implement it in APL too!\n", "\n", "Super bonus: Invent a new meanshift algorithm which picks only the closest points, to avoid quadratic time.\n", "\n", "Super super bonus: Publish a paper that describes it :D" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 2 }