{
 "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.\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": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math, matplotlib.pyplot as plt, operator, torch\n",
    "torch.manual_seed(1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "centroids = torch.randint(-35, 35, (n_clusters, 2)).float()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.distributions.multivariate_normal import MultivariateNormal\n",
    "from torch import tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample(m): return MultivariateNormal(m, torch.diag(tensor([5.,5.]))).sample((n_samples,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "slices = [sample(c) for c in centroids]\n",
    "data = torch.cat(slices)"
   ]
  },
  {
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_data(centroids, data, n_samples):\n",
    "    for i, centroid in enumerate(centroids):\n",
    "        samples = data[i*n_samples:(i+1)*n_samples]\n",
    "        plt.scatter(samples[:,0], samples[:,1], s=1)\n",
    "        plt.plot(centroid[0], centroid[1], markersize=10, marker=\"x\", color='k', mew=5)\n",
    "        plt.plot(centroid[0], centroid[1], markersize=5, marker=\"x\", color='m', mew=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABJ00lEQVR4nO2deXxU1fn/3wdCQhbITliSkABhV0ACAi6EVVHUFLUtLlW0xbpSq1bEfov0V1usVWvdcaG2WlwRFWtFkLBvYV+FQAIkhBCSSUIykw3O74879zIzmYRAEkKG5/165TUzdz3X5XOe+znPeY7SWiMIgiD4Jq2auwGCIAhC0yEiLwiC4MOIyAuCIPgwIvKCIAg+jIi8IAiCD+PX3A1wJSoqSickJDR3MwRBEFoUGzduPK61jva274IS+YSEBNLT05u7GYIgCC0KpdTB2vaJXSMIguDDiMgLgiD4MCLygiAIPoyIvCAIgg/TaCKvlGqtlNqslFro/B2hlPpeKbXP+RneWPcSBEEQ6kdjRvLTgN0uv6cDS7TWScAS529BEAThPNIoIq+UigWuB95x2XwT8L7z+/tAamPcSxAEQag/jRXJ/x34HXDKZVuM1joXwPnZwduJSqmpSql0pVR6fn5+IzVHuJCwlxSz4avPsZcUN3dTBOGio8Eir5SaCBzTWm88l/O11nO01sla6+ToaK8TtoQWzs60xSz/cC470xY3d1ME4aKjMWa8XgHcqJS6DmgLtFdKfQDkKaU6aa1zlVKdgGONcC+hBdIvZazb59liLylmZ9pi+qWMJah9aGM2TRB8ngZH8lrrp7TWsVrrBODnwA9a6zuAr4C7nIfdBXzZ0HsJLZOg9qEMufHmsxJoV4unrjcBsYIEoW6asnbNbOATpdS9wCHg1ia8l9BCqS1KN4Ud6n4TcD1uyI03n4cWC0LLolFFXmudBqQ5vxcAYxrz+oLvUZtIuwq7+SYANTuFhlpBguDryIxXoVnp3Ksf4Z1j6dyrn2W9FORks/m7hVSVl+M4cYINX31Ozo97mD/7GdZ/9bmbdWMK/c60xXVaNmLrCBcrF1SpYeHCxoyiuw2+nAMb19U6EFqQk82yf7/DyDt/SWSX2DqvlbV1M7Yj2Xz35t/pNeIq1n42j+1Lv8d2JBuAowf2kbk5nfDOsdiOZHOyqppht0ymqrwce0kxQe1D62XZiK0jXKyIyAv1xhTKw7u2k7nZqPvvTTCX/fsda/+k6c94vdbm7xay9rN59EsZR2D7UGxHslEaEgclk7k5nfhLBtKlZx96XTGSuL6X0LlXP9Z98REj7/wlBzauY/mHc2nTti1DbryZboMv5/Cu7XQbfLnbPVw7G7F1hIsVEXmh3pgC2W3w5cT1vcSrYNpLionsEs/J6mpG3vnLWgdWlTY+TxzPx1FSTOKgZAZeOxHHiROcrK6iQ9duDLx2IkHtQ4nsEmtd9/u3XyUqrivDb55s3f/HVcvI3JxORGw8wU77Jqh9aI3OxrVDkrRM4WJBRF6ot+C5DoDWZsPsTFtM+sL5XH37FCK7xLLB6aEDlnfeL2Usva4YSc7eXbSLjCb+kgGMvPOXlvVyaPtWdqdvICj0tN9eVV5O+sL5lFVUkrN7B8NumUxQ+1DsJcXk/GiUTMrPPMDGHVsA4w1j5J2/BLA+XZ+1qrycNZ/Ps44VBF9FRF5oVL/a0xZx/XS9D8Ch7Vut7wv++kdSf/cH+qWM5c1/f8hXy1fRM2U8VeULWfP5POL7D4SIDrzw/n+4PDGe4TefbvuhHVtIHJRsWTnmPSO7xNawi8w2DLtlcg1vXxB8ERF5od5+dX0ifm/pjt0GX259ut6nqrycrG2byd23h6KjR/jhvTdZtGsv7332BQD3zfg9f5w6hfhLBrBx9UreXbuVEkcF3+/aR//la8j5cTfhnbsw/ObJDLx24lk/qyn4prcvCL6IiLzgJsx1cbYRv3l81rZNHNq+laryckb89HZr/4if3o5WkLtvD6EdYpj33WK+WLXO2l/iqOAPc+by69QbeHvlRmwnSq19L732Gjv6JnFN/55cffsUgtqHullDQ268GXtJMVv+txCtYNA1hr/v+qyuYwwbvvpc/HnBJxGRF+rN2WaomMfZi4s5tH0rh3dtJ+2D96iuqMCWm83oKfcz6JqJ+Ae0JTc7hyXvfuh2fiqppDnS+Ou8zwgjjFRSWcACa3/6oVwmXTMOe3Ex9pJit/bZS4r55uXnOeT06JWGNm3behXyPauXsfYz8ecF30QmQwn15kw1aOwlxaz+5ENWffqh26Sj/qPHkzgomezdO9j49Xy2LvqGQ9u3suzf71jX9Fea+1OG0T6wLWAI/DSm8RIvkUACL/ES05hGqnNZgvaBATw8IYWynIOkL5zPhi8/d2uf6dUDhHXsTFVFRY36N+abhtJw9e1TJL1S8Ekkkr9IOdcUwromRG3530IrY8U/wBDr5R/OxV5czMnqagaMvx6Aoxl7cZSWMHD8DWz46nO6Db6c/RvXEt0uhIfGX82ri5aT5kjjJm4igQTmYlgwWWSRRhrtAwN4YPQVtFca1ao1+tRJcvfvdWtjVXk5yRMnkbt/Lzm7d5C7fy/DbpnsJuSepRMEwRcRkb9IOZeMGlcLJGvrZitSNs/Xyjgu/pKBbmJqevKt/fzo2C2JvAP7ANj43/kc2r6VrG2bKD9xAj//ACKAWwb3572V6cxiliXwALOYRRFFPDL+WqICWtO2XTvKT5wAICquK6s+/RCljXas/XweV98+hfi+l5Czewc5u3cQ3/cSt46tvmMRgtCSEZG/SDmXGaCuFkh0YjcSBgxyO9/0100RtZcUU1lRTninWEAx8s5fsmf1MgBCO8RQVV4BQIeu3enSsy9Z2zazbVM6n23cThhhzGSm2/1nMpNHeZRP1m7m+d8+Qvnh/ZSfOEFYx8608Q+wfPW27doxYPz1bm3TCjR4zdmXKF7wZUTkL1LOJYrtlzKWyopylMaajVrXNXemLbaE15wcNeiaieTtN+rRFB/LI77/QNoEBGA/UcK2Tem8kbaWEkcFqUwggQSyyGIWs5jJTBJIIIUUFuQt4JE//YWpVw6mU1QkRUeP0CYggGG3TGbzt19TfuIEB7dtZuy99wNYGT32kmKrE5JaNsLFgtJaN3cbLJKTk3V6enpzN0NoAK5ePxg1ajw7BXtJsbW9srKCjV/PR7ULZdZ/PqfEUWFdK5VU0kijiCLCCDME3iW7pn1gAI+Nv5peAy/jht8+RVD7UBa/+wZbF31Dv5RxRHaJtdrhGrW7plb2HjGyzmJrgtASUEpt1Fone9snkbxQJ2c7QOsZIV9x6+01jglqH2ptX/WpkTZ56fAruXrbLhZu2GIdt4AFhLcL4Z4ByXy+aScL7AvcrnN5YjzBAf6crK7mf6+/xMg7f8mIW28jNLoDlRXlLP9wLpUV5dabg9mmnWmLWeP07M1iZ+Y+QfA1ROSFOqmPrVGQk80Pc9+kQ9du9B89Hjidq36mDmLQNcZMVaVh8vgxVJSV8f0uY2A2KjyMJ3/2E1ThMaLbBfPu2m0ct9kAuCXlKoZ1aE9YTGda+fnVKESW8+Mefly9AkdJCZmb00kclFxryQXP74LgS4jIC3VSnwHaZf9+h0Pbt3Bo+xaCQkOt2ab/e/0lMjenc3jXdre6Mq6CH9Q+FP+Atiz/cC5d+vTjmv496dg9if+tWss//u8pctakATBw6DBeGT+R+2b8nmtHDGNoZAiJg5K59oFHcZw4wfdvv0qlo5yCnGwiu8Sy7ouPsB3Jpl1klJUDb97Xc+xAInjBlxGRF+qkPgO0I+/8JSerq+nQtZubB24u9mFG2eanZ2aLeY5ZM/7uO3/JS20DCWzjxzcniji0Ywtx/S5h0DUTiYyK4rJRY906jKD2oSgge88OFr/9Kj97ZrZVefLyn/ycIz/ubKJ/OoJw4SMiLzSYyC6x3Pr7P7ltc60Lc2DjOrca9Oakqaytm7l+2hNWR7Lq0w/J3JxOTPckBl0zkZ1pixl9z6/Zs3qZVX9+3B13AxAaFEzxp5/hP2kSfuHhmOkD2qVNk6Y/U6OejSBcbDRY5JVSbYHlQIDzep9prWcqpSKAj4EEIAv4qdba1tD7CS0Db7XnzU9z0tShHVusiH5n2mKqnXnzSruPBZh2jmu1yOL58zn2/N+M6957L+N+9ZC1CpQrZzMfoNpmo3j+fEKdHYcg+AINTqFUSikgWGtdqpRqA6wEpgGTgEKt9Wyl1HQgXGv9ZF3XkhRK36TaZiPvo3lkR7Tnkgk3AEYJhKqKCtoEBBgzVD+bx7BbJlszVl1TGwG3tMydaYvpPWgIlYu+x3/8OPZs3lDv7J+6hLzg3Xc59vzfiHrwQVoFBYrYCy2GulIoG1ygTBuYNWDbOP80cBPwvnP7++CsLCVcNNhLitnw1efkfTSPLf/5Nyvnf2QVCGvTti1tAgJY8/k8q0DYoGsm0qZtW9Z+No8DG9dZxcY8C48t/3AuezZvIPLee9mzeYNb4bFqm42Cd9+l2ub9pdF8AyieP7/GvtBJk+jwxOOAdjvmTNdsTgrLKnlr2X4KyyqbuynCBUqjePJKqdbARqAH8JrWep1SKkZrnQugtc5VSnWo5dypwFSA+Pj4xmiOcIFgCvKVk37OwNvuJCqivdts02G3TGb4zZPRihqDsJ72ihmB9x4/jlMOB51z8qi22Ujq1pOitu1J6tbT6FRemE3YFwsBQ7TNqB0MgQ8ZPdra54lfeDiR995Ltc1Gq6Agt/NcraELiU/TD/OXb/cAcN/I7s3cGuFCpFFEXmt9EhiolAoDvlBK9T+Lc+cAc8CwaxqjPcKFQbfBl3N413Z6XJlCZJdYuji3e1udyT+grRWtmymY5kIe/idPcWT6dMqWLSeioJDItKWUHMikrX8AZevXE7tuMyWvvEre9eNJ37eT5J9MpJ9T4E1xBjj2/N8oW7+ezrNnWzaMN/vGFHsTU+y9dQzNza3JcW6f50JhWSWfph+2rmF+jwj2b5Q2Cs1Lo2bXaK2LlFJpwLVAnlKqkzOK7wQca8x7CRc+BzauI3NzOnF9L6l14e/aIvdNCz5j3TdfUFFURJ8qKFu2nOCRV+PYvp2qA5m0iY8ndNIkQkaPJg+ImT6dqMgI61p+7UNriHPpqlWULVuO7YMPiH74YaB+Ubqn6F9IRAT7n3MEb4q7vfIkLy/ZZ22XNwPfosGevFIq2hnBo5QKBMYCe4CvgLuch90FfNnQewkti34pY70uxmFG7zvTFrv57a7ed/nuXQCUfPcd1YWFRD34AJ1nzybwEuMlsd3YMfiFh9M6LIzgoUNpHRZGUPtQBl01Gsenn1Fts+EXHm5ZNgD+iYkAVBcVWW0xffiQ0aMvWN+9qTht9WiemtCbW5PjuDU5zvou+AaNEcl3At53+vKtgE+01guVUmuAT5RS9wKHgFsb4V5CC6K2iVSde/UjvHMsnXv1c9vuGlVf/sA09O+n0yF9K4UbtxP14IMUz59P2K234hcZWatfbv4+ZbcDULZ+A44NGyhduQpdXQVAxd69FLz7rmXRRN57r5VZA+5evmd2jS+lWY7tG8PaAwXcOLAL3aNDrO0SwfsWDRZ5rfU2YJCX7QXAmIZeX2heznUFqbowSw6s++IjJk1/xtruaq/4hYcz6u+vYfvgA0BhZrwA1uBowbvvWgOpIaNHk/f88zi2buPUT1IpW7ceh0s6rn3NGhx9+tA+Ph7/+HhnR+Ag+uGHrPPL1q8nZPToOi2cC2kQ1tVLPxv//LRNU83SH/MZ1i2P7iNDznxiE7RFaHpkxqtQJ01Rd92csOQ5ccnT+/YLD7e8czPjJWT0aPJfeRX7pk3Y16yhuqAQx7ZtFH74H6qPHOHV4/l8VlLCP2PjSPA3xMa/d2/2H8vjzoVfc0toKNOqjIi+MjeXQ/fdR8z06ZT+8ANly5ZTOnRonQOtrp3B+cZTSN9fncXLS/ZRUFZJZLC/m8DWJbqmTTNtTFKt1kxd53vbJ1k+Fy4i8kKdNMWMUbPkQH1wnUjVd/w48v70LGXLlgMQPPJqynfvtiL2V4/n83pBAQB3Hz7EP7v3IAHYl5nJXbt2cqy62ti/fRsPRUVT+sMPnCoq4nBmFu3GjiHqwQfdLBzXZwoZPdo43u6wOoMAL5F8U9o5NYXUSEbbmVPMqv0FLttrF93CskrsldVMG5PEjQM789WWHN5fncVdIxKs88b2jeFPC3ex9Md8t/MLyyp5f3Umaw8UsC7Thr2ymrtGJFrnwJmzfCTiP/+IyAt1cjYrSJ2LleEqigC2Dz4ENOF33IFfeDjF8+ez5T//Zk/nSMq3b6fDsuUEDR9O0GWDCL/jDk4WFXHkqad4cd06S+ABjlVXc/f+DGbFxTNzfwbHqqutfa8XFNA6NJT7i4poFRFO1aFDFL43lw5PPF5DmM1nKlu/nrJly4l68AE6PPF4remUTWnn3Joch72yGnvlSQrLKrlrRCJB/n4MSYiART9SUFbJ/vxSFu/KqyG6rhbNy0syeGpCbxbvyuPlJRnW9TcetLEy4zgr9uWzMqOA7tHB1nUA55uD6/FFQJaVmVOfCF4i/vOPiLzQaJxLPrmrKJ6yOzj+2msAtAoKIvLeewmdNImBlRVERbSn74irqUzsaUXVACeLi8nbv59PsrPdrptKKmnVadyfeYAwwkj1WFXq87IyHvzNVNpXVlCxL4PAS/oTOmkS1TabW0djPkvI6NGWlVNXhN6UOfURwf4E+fvxl2/3EOTfmvtGdue+kd15a9l+Vu0vYNX+AvblnbAi8FuT43h/dSbmmMbLSzJqWDT2ympr/8qM4wCUOKroGhHE/vwyZn65k76d2hHo74ej0ugoL08Mx69VK1ZmHGdw17CzysZpjLx+4eyQ5f+EZsUQVWNwtdpWSNF/5hE4JJnYf/yj1glLZiZM8MirqczKourgIQ4GtmVKZiZHi4tJJZVpTKuxPuzLvMwCFtChTRs+m3IPvfv3tzqVDk88XiPLxnhjuIzwO253E/bmzLDxZneYNgoobhzYmcW78rg1Oc4tak6IDGJ83xh+ndIDcJ/wtD+/lBnzt3G40MGR4nLrXhHBbSgsq7J+X54YgV8rxR9T+xMe5C+2ywWELP8nXLD4hYfTKiiIY8//jaDhwwEIHno5fuHhlpiaEX51YSEVGRlE3X8/QcOHU7ZsOQG9e9MqLIwr/voc31dVMfaGG0hzpHETN5FAAnMxBo2zyCKNNDr4+fHP2Diili3DXlVJ2G23UZmZScjo0VTbbJyyO4i45x7K1q/DvmYN9jVraBUUSOS991KRmUne7NkEdO9B4XvvAec3w6YuPzvI38/a3n1kiNN7P8nUqxL5ZnsuWQV2tmYXExHsz1vL9vOn+euxV54kyL81K/YdZ12mMT+gXwRc0b8bWw7ZWJ9l4/LEcAbEhrEr94QV6X+1JYdHx/VyvilkAZq7RhhzEET4LzxE5IVmx80S+eGHGjnwpg9eumoV9tVrcGzdRptYYwZtxR4jUq3ct49OwDMRkdyfk80sZlkCDzCLWRRRxBsxsUbWjb8/9tVrACO9smThQloFBXH8tdfo8MTjnCo5AUCrsDCrPXmzZ1O2bDm6ssrNl/cW2TdFtP/+6kxeXpKBvbKaR8f1srabEbur9QKKl5fsY9qYJFopo7az1pr9+aW89sJfKF71NT9Ev8L20mCmXpXI9pwijuccZPEbM+j/4P0MvPI21mfZqDqpQSn6dmpP1clTrMssdN7DuK/px2/LLubS2LCz8ueF84OIvNDsuM5MdRVFz7z5kNGjybr9dk4V2qhwmbUaNHw4oZMmsWfTJp4pLCCMMGYy0+0eM5nJozzKrJJi3o+PI67Qhl9sLLrK8Jmri4qo3LSZiHvuIXTSJAIHD+bIjBnETJ9utSvq/vupPHyY6GmPEDRwoHGezWbV1YHTkX3TDMAqj0+D0/766fIECZFBTL26G6A5bHPQPTqY6RP6MPb2h8j63igO++1z93PZ/S9xbf8RrEjfxs55T3GytJC//uVPjJ5cCPHXselQEZsOFQEwbUwPRvfuwJCECO54Zx2RwW2ICw8kPMiPpT/mk9QhhGljkrBXVlNYVinR/AWCiLxwQeAqiq6C75rKWLLwG9pfO4GKvXsJ6NWTVgEBtAoMIvyO29mbkcHYG28kz+EglWtIIKGGJ59CCgtOLOAX27Yzt0MHegYH49iwgcDkZMpWrqTq4CGq8o4S+atfEjRwID3++183j/6U3UHVgUzKVqy0RL54/nyrro7rYGtTDMDeOLAz27KLuHFgZ7ftZv0as9zwl1tyyCqwsy/vBC/8dKBl5Uz+9WOWwAOcLC1k0xuP8lDpUdb/+zlOlhZa+36Y9zqhV9joOu4urunXkU6hbblxYBcW78rjxe/3WtaNQSAAu3JLGNw1gpeXZBDk71cjdVOsnOZBRF64IHAVRW9RcPH8+dYgadSDD1iTpAAKCgoYM348ec5SBmYWTRppFFHEozxqCLxze15ZKXcfLufLsFBCgYp9ezlVXGKkUx7IJOexxwgaZAy4urbLGCAG+6ZNbrVxzP2utkxTFDVbvCuvzhmqEcH+PDquJzcO7MyfFu7i9xP7EuGcJPX6/zax4uuP3I5PJZW00jTWvPEkYYSRQqpbBpJj+yIKB9/A4t2KT389gsW78vjLt3uYepXhv5c4KtmWU0J0uwDiI4JYmVFA307tGdUr2i31EiR1sjlpcIEyQWgMTFE0hdMzFz100iRrYNa+aTPVNtvp0gYlJdwxbJjb9RawgICQav51++0EBFW6iRfALeHhhNodAJwqLqFN13gSPvyQ4JFXY1+9huOvvUbx/Plu7Qq/4w5j/5o1HJk+3RJ6c39TM7ZvDKN6RTMkIaLOhUK6R4cwd8pQqx7Np+mHeWdDAeE/fZbAsGgAKwPpJV4igQRe4iWmMY1U59o+rUMiGPPYq0RFRlJYVsXML3dYA7mB/n78Y/Ig/nnP5YzqFc2mQ0UopZg2pgeB/oZ1s3hXnlubpPBZ8yEiLzQLda225OrRm2JePH8+Hf/wf5bIFs+fb0X8hx98kHsys3jihhusa3QMDeWj664jOX0jc2M6EhMUZO17IDKShyIiaRMfT1un7dJu7FgCEhPpPHs2EffcQ9Dw4QQOHuzWRr/wcDrPnk3wyKspW7b8vK8cZUbyr/ywj798u4dP0w/XOMbbSlFj+8aQEBlEm4gujH3sVdpFdCCNNLLIsjKQTHsrjTTaRXRg/BOvs8vejrjwQK7sEUnfTqG8vGQfn23K5uUl+3gzbT+PfbKFh0cnMapXNCszjhPk78ddIxLcfHkT01ISq+b8I3aN0CycaWDSdaZpYP/+HH/tdQA6z55t+fUni4qwff45VQcyCR55NX+ePRv/Z57hnXfe4fOpUwn/6mv8unShh1LMBabkHeXmwEAeiorGL7YLcW+9SeuwMJc8fSMy94uMwL5mDcf929QYUAUI7N+fgB49OGW3Wx3Q+ShaZkbBY/vGMKxbnteo2LOezZCECB79eAsHC+10jQhiWyEMvv1J0l55rNYMpF8++iYl0XGQaWNbTglX9ogEoGtEEAcL7UQEt2FrdhHrMgupOnmKwV3DuTQ2zGrPtuwilv6YX8OXNxF//vwiIi80C54Dk2YOesz06QQkJhI6aZJVSiCgRw+CR15NyOjRbl538fz5lsCbqz09eumlXN+5Cz3i4rCPGG6lSfa5YgRLu3Un8MABHBs20KZTZ0oWLiT8jjsAZfn90Q8/5HWWq4kxNvC6Fc27LhPY1CtHuS4QYnry+/NLLf/dsGeMyY1bD9tYl2kjIdIQZoCRvaIpPXaYN574S50ZSHOfm8EjL3zAFd0jqT6lWZlRwMqMAqZe1Y3vdh7lYKGdXjEhtGmtqDqprTIJAI/M28zKjONc2SOyVmtG/Pnzi4i80Cx4DkyaOeh5QPxbb1nWiDEZyl6jKNjpiUtTaBV42ooJnTSJHnYHoGnbp48l8qDg448JfvABWgUFUrZsOY4NGzidV4716do21yJk5j2jHnyA9hMnupU5aK6yw6cLie1i7pShVj2bgtJK1mXauDopmpOnjnHY5qAgO4t3ZkyhsuQ415JaewZS6QL+M/Nevvp2EQsPwpEiB8O7RRLo34pr+ndkzvIDtPX3o03rVqzMyOfKHlHYK6t5f3WmlXUzuGt4rVG6lDY4v4jICxcEMdOnW8v4mdS2sDaczrY5HVEHuqwCpTn+2uuGr56cjCM9ncBLLqFdykhCJ00iHMj57WPY16zBrFHjeX3PwmlmZ3P8tdfp8MTjBCQmeq1Ceb4wI/i7hicAcNfwBG5/ey39uoTysyFxfLzhMJcnRtC2TStG9e7A3B+2M/fpR6gsMUT4jBlIR3MZP24s7W57idaB7fE/ZGP/xjKmjenBUxN6W7XoR/WK5tLYUKsuzrQxSbjOgHXF1aY53xG8o7SS3atz6TOiE4Eh/mfc7ktI7RqhReJZAtg19TLqwQdx7NjurBr5oNUBADWEu7YZqWZ+fIcnHgewrmteq7lXhZoyd70lsnOnDLV+A4QHtcFmd685sy6zkK5Z37D84zfcrtM6JIKIax+mdPFrOIqOu+17+PGnyOtxA1prpk/ow4asQstHP5fFv99atp+/fLuHpyb0Pu8iv2nRQdbM38/wSd25bHzXM24/W5q7s5DaNUKLx7NMgGuUb+LqjYdTU8RN4T7ltHNAcbKoqEZn4Zn/7nr95hZ3k99P7Avs4uHRSby1bD8Pj06isvoUhwrtHLY5rOO6RgTRNSKIdZmFXHPnwwyIC+OVv/0FgDbtIkm86zkqgjvSJvxZ/L/4A8XHjdTH3z31e44l3cC6H/OZNibJTeCh5gLirjXnaxP85rRp+ozo5PZ5pu1ny+7VuayZvx+gQZ1FUyAiL7QIastgsX3wAcdfe51TdjvRDz/sts/TJzcF27RdAGuFKXOQ1zzP02dv7qX+PDFz4V2j41duu4w30/azNdtGr5j2HDheysqM0zX2A/1b8Y/n/wzAu++8w8hH/8EuezsSIoO4acxIrnhgMTdddw2Dx99Mx1F38PGSDEb1igZ0rQuQnM0KUZ4dQ2Nypkg6MMTfq/i6bve8xtlE543VWTQFIvJCi6D2DBbv9Vy84Rr92zdvxr56DW379Cbkyiu8ZtJc6Liu8mSWFp6z4gDTxiQR5N+au65IsAZmR/WKtnzyab+bQX78KG67qi9/X7KP34xJYv7mHG4c2JfH3/ySV1bnMQblNnnJLI3gijdBb65o/VwiaVcRB1j8z10c2lFoXcO8ZnXFSfwCWnsVe/MaiQOiGvFpGpcGi7xSKg74F9AROAXM0Vq/rJSKAD4GEoAs4Kda66adLSL4LLVlsITfcbub5+4Nb1ZPlxdeqGHneGbSNFfNeKhfLrlRBdJIXzTLF4CxEIgpvi/8dGCN68z8cidrcirZ89VObPYqZn6901k3fhcv/HQgIaE172vWxnlr2X5rn2vevuv25kiLPJdI2rVjADi0o5D4/hE1rlVVcbLWDsS8Rs5eG4d2FHJ4VyEdu4dy6ajYC2YgtzEi+WrgMa31JqVUO2CjUup74G5gidZ6tlJqOjAdeLIR7icIFvVJX/Rm9ZzpvPM1wak26pNL7hk1uxYqc60v73l+307tWJlxHJu9yloYJCK4DQ+PTqp3m8w3B8/FSZor7702O8YbntG3a8fgGq2b13SUVtLGGcm7nt9nRCdrW6ceoZTkO8jeYyN7j402Aa0vGG++wSKvtc4Fcp3fTyildgNdgJuAFOdh7wNpiMgLzcC5TFY6XxOcaqM+tkdtUfOZoulfp/Qg0L81Gw8WsTLjON2jg9mfX8ayvcd45Ydij+UDszBTIl3bZM6stVeetBYBv1Dy3uvy0h2llTVsGU/bZtOig17F3sTTGrpsfFc2LTpIUZ6D2N7hdOweSuKAKNZ9fQAFXNLMUX2jevJKqQRgELAOiHF2AGitc5VSHWo5ZyowFSA+Pr4xmyMIwLlVhGzOCU5Qt1A3tCyAUa2yl3PZwCwclSeZeGkrQFn+vRmhuy4K8sJPB7q06fQEsuayaGqjLn9+9+pcy5ZJHBDFpkUHqao4Sfo3WdYxpv0y9u6+AGxbmu0m1n1GdMJ+opLDuwpJHBBFeEwwiQOiyNlr48pbkwiPCWbTooPWNf3OENU3dfplo4m8UioE+Bz4jda6RKkzD4QBaK3nAHPAyJNvrPYIgq/SWPaIsTB4a15ess8aZA3yb+3mudsrT5KeVcjSH/P5NP2wdT9zZm1zRu+1iaOnp+7NXukzopPVGQwcF0ds73DsJyrpd2VnDu8q5NCOQrYvzcYvoLUl1rn7ixn/y34EhvhTcLiU7D02vnl1Gzc/OZjMrcc5tKOQzJ7HCR8fTJ8RnaiqOInizOMETZ1+2Sgir5RqgyHwH2qt5zs35ymlOjmj+E7Asca4lyBc7DRmBovrtTwjcrM+vefEJ3Nfc0XvpmhXV5xkg1OAXcWxLnvFFNzysirsJyqJ7R1OdeUpy0s/WXkSgNLyYrTz+EO7CsjZU0T2HhsrvtxOVHQUoTGBZO+xUZzvsPz9nL02OvUIteyey2/odsZn8Ox4moLGyK5RwLvAbq31iy67vgLuAmY7P79s6L0EQWg8ga2v7dNcgl5b+0zRTr4+gSHXJ2A/Ucn6rw/U6n17i97NbBiAkoLTk8cO7izko/+9xbqMb5kw7QcCQ/wZObkXS/+9h0O5mfz81z9jWK/rePKxGQSHB4DW2I6WcnhXIdl7jOTBQzsKa6Rder51ePP1m4rGiOSvAO4Etiultji3zcAQ90+UUvcCh4BbG+FegiA0EhdCVkxd1NY+T9He4MX79hxMBSN6d5yopFOPUNpHtmXAuDjLdunUIxSAj757i283/guACROv4duF33FwVTVb0nfwj68fo9hewLcb/0X799tyVcLPANiz2pglHBYTSPJ1CXTpGV4j7dIU9aqKk7QJaO01s6epaIzsmpXUPhNlTEOvLwhC03ChV4OsrX2udkxt3rdrxB7TtT0bvslyi95zM4oJiwnkiluTaOWnuPLWJH73mxl8uvht6xpHjhxh7Lix3DzkET5Z/XeK7adnD3/83RwcIyoZf8mdRMYGUV2hKcpzkJtRXCPt0lFaSVXFSQaOi+Po/mIr4j9fKZZSoEwQBJ/CUVrJ9qXZHNlfRM6eIpKvT0AB9tJKbEfKiOgczKFdhZTklxMaHUhxvoOQTprfvjiZwpJ86zqppFpVOo01cFPclpGMCI3md6lv0SWuI8X5DuL7RzD27r41LCOzCFqX3mHk7Ckitne4NYDbWNRVoEyW/xMEoUXgKK1k06KDOEq9r21rYlo4nbuHMeT6BMBI+Ny57AhH9hVTUlDO2Cl9ie8fQWxfYzZzaa7iwQnPExZi2ChnWgM3OiKGF2bMJSoykuJ8B2ExgVx5a5Il3K5t7TOiE8MndSc6rh0AUXEh7F6de8bnaCxE5AVBaBGYFszu1bl1HmeK6iWjYq0USAX0G9mZtsF+HNpRSG5GMWPv7ksb/9aERAYAEN85kYev/xsRodF1roEb3i6K+699jrIDgZSXVdM2xI+iPAeZW4/jKDUGgRe9s9Nqa3lZFTl7bXQbFM2AcXHs35TPmvn72bY0+zz8UxORFwThAqO2iN0Ub2+Dla7nmJ69OXFpyPUJaKA4z2GJcuKAKHavzmXL94cpLaggLCaQCns1MWFxPDH1jxRRxCxmud3DXAP3Z1c8Sod2xjhB22A/xtzd15pctW1pNhu+ySJ7j43Y3sYA7PJ5ezm0o5D0/2Zhyy3jREE5AEf3F5+XaF5EXhCEC4raInZX8T6bczSQ/k0WoTGBhMUEUl5azcJXt3KiwMHAcXEkX5/A6Lv6EBIZQF7RYZ574/9qXQM3jDA+WvkSx+2HCQ7z57oHL6XwSBmHdhSy8tN9VDvz7GN7h9Opeyjp32QRFRdCfP8Irrw1iStvTaJL7zA6JYWSvcfGtqXZbFp0EFteWb2sqHNBSg0LgnBBcS6Tg+o6x0z9CwrxZ9ITg5n//EaK8hzsWHaE4ZO6kzggiv++vo39+zOsNMnUM6yB+9o3T/Lrcc+xf3MM/v6tie0dzqEdhUR0CrYEvW1wG68lilN/c5nbhC7PvP3GzrqR7BpBEHwa16qTe9fnUWor5/DuQkLCAxhzV19WfrqPXemZzP5sKkVlp5dAPFN2TWhQJK/M+JSSw5ouvcOIjmtH/uET5OwpokvvMDp3DyNpaAyZW4+TOCCKzK3Ha11jtrb99aWu7BoReUEQfBbXSVGe9eMBoruGYC+pAjSffP+2NRHKJCI0miem/pEX332G/MI8t30/u2YqL73+HOu/yrQ8+Ow9NtqG+FFeWg0YE6SK8oz0ykM7Chu8lmxtSAqlIAgXJa5efZ8RnRgwLo5OPUKJ6dYegPyDpZTZKvDza82f/vxHfvfYU9a50ZExPHjt8/SLH8avxz1HdGSMtW/C4F9wVcLPyM0opmN3Y7ZsZFwIsb3DKS+tplOPUGJ7h1vlhyM6BZN8fUKzLA8onrwgCD6LOSO2usIYEA1q509uRrFT6NtRXXmKSsdJxkzpQ6duYVw2/s8Ehvjz9ttvs+DTbzi6CatUwS1PLmbM2DEMSbyWe376MB27h1JdcZKeQ2Os2a3blmZbM1qj4kKIigsh//AJtnx/mOGTunu1YlpMqWFBEITmojahDAzxp01Aa9bM34+fs2bM7lVHyM0oto5Jvj6BTt3CrN/PPPMMDz/8MBmrizm0I4sOXdtz+Q3dcJRW8uE//ktY+3AuGRVrvSW41szpOTSGjPQ8cjOKyc0oJr5/BDl7itzq13u2sUWUGhYEQWhO6hJKz4JmpoVSVVlN3oETVFeerNFJBAW04+j+Q8DpfPbdq3PZt6KY4ZOMWbFVFScZ4mLBOEorWfnpPoryHO4Drz2P1xgTcG3jBV9qWBAEobmpSyhd12o1C4X5+bemuvIkeQdO0Ma/dQ0B3r06l+w9NkKjA6189p5DY8jZa7MmUqV/k0VYTCBJQ2Os8sHmqlOuNWzCxwfX2cazWZ/2XBCRFwShxVMfoTSF2cx0GXJ9Qo0ZtOZ3cxGQ9pFtKV7mQIHb6k99RnRi96ojFOU5WPnpPm54aKCbiHvz1ptazGtDRF4QhIsCVwE3hdpVjF0F2BT05Fo6gsAQf6574FJWfrqPK29NAuoW8aYeXK0LEXlBEHwaV4E1Rdi0UGrD01rxPB8gPCaYGx4aeMZ7elsJ6nwiefKCIPg09a1e6YprnRzP8+tT8tjznLqKqzU1EskLguDTNDR7xVtUf6aovM+ITlRXnKSq4qRbZczmQEReEASf5lwEtq6aMvXpNAJD/PFz5ue3ccmjbw5E5AVBEDxwXSPWszpkfTuNps5/ry+NIvJKqfeAicAxrXV/57YI4GMgAcgCfqq1tjXG/QRBEJoSb5k4Z0tzWjSuNNbA6z+Baz22TQeWaK2TgCXO3xcPZQWw6mXjUxCEFoUp0OExwbUuVNJSaBSR11ovBwo9Nt8EvO/8/j44V8C9WNjyAXz/B+NTEAShmWjKFMoYrXUugPOzg7eDlFJTlVLpSqn0/Pz8JmzOOXKuEfnAO2DcH43Ppr6XIAhCLTR7nrzWeo7WOllrnRwdHd3czanJuUbkwZGGwG/5oP6ibd5r/Vsi9oIgNApNmV2Tp5TqpLXOVUp1Ao414b2aDjMSP5uI3MQUbYArpp3eXlZg7Bt4h9EZeN6r0u79PEEQhLOkKUX+K+AuYLbz88smvFfTERx57kJbWwdRm/ib9L8F/IPOrWMRBEFwobFSKOcBKUCUUiobmIkh7p8ope4FDgG3Nsa9WhS1dRBnEv+slZD6pnuULwiCcA40ishrrSfXsmtMY1y/RePNmqlL/LNWwr5Fxjli1QiC0ECafeDV56lr4NYzmyY40ojgzzYrRxAEoRakrEFTU9fArTdv/lzHAGobzBUE4aJGRL4xqEtg6xLtujqAsgJYPwfQMPS+Mwu32WFU2k8P2orYC8JFj4h8Y+AZkXsKNBi57ygYOtVdfCvtxj7zOLOzWD8Hls02tvkHn865r028rfTLsvqnX0r0Lwg+j4h8Y+AZkW/5wF2gAZY95zxYw6gZ3o+rtBu/K+3GcQDdUqDndfD5PXAgzRBx83xPkTY7GLNTOBNnSuUUBKHFIyLvDW8R7tlYMj2vg/0/QMdLT4vt3u/g4CqodJw+buAdpwV94B3OaB84vBYm/A1Qxr4dnxkCD85tThrq6TdkopcgCC0C38muacy6L94yYs6mvMHe/xqiHBxl/DbFG0Dp0+0MjoRRTxmReXAkdB8HQVHGuYtmANp4A6gqg/groOsIY6KUybnUx3HF7BDEqhEEn8V3IvnGtB68Rbj1jXrLCozofOSTLt6606rplgLZm+DQKuOYoVONDqDSAf6BcHgd2I8bQr9vEXS+zBDxSrtxDhgdSHQDM3EEQbho8B2Rb0zrwZt4unreS/9iRNdtgg2hhtNWjumzJ413nuj01iO6QUQSpL9Ngf0UkWjnsc9Ztyjoci2REd2MqD20y+lB2kMbYOs8CIkxrCBPZABVEIRa8B2Rb+qo1hTSyjI3YcY/yPg03yI8Z60OvQ+ObDZ+h8XzTFo5b29txdLbxtAzqYfh3R9IY2/BSUa99Dm/GtSaZ1IOwEjnGiurXjaOKTpo/K1/C/Yvge5jIOUp47llAFUQhFrwHZFvakwhHTkduo2CA0uNT09LJzgSxv/Z+B03wjhv/J8h4UqeWZDBrGWVAIy6fhJLly6l583vsfef0xj10kccKTnJrGUnAXhmpD59zxGPGNfreCns+tIQ+8IDEBRh+PkygCoIQi2IyNcXVyEdOrWmPeIaQe/9rxG52zLh+D6oLOOZ9/7HrA9WWIccOXKEUaNGMWfOHKbO/JojJSetfbOWVcI/l/DM3950v/f6ORAYYYg8YGXaiDcvCEIt+E52TVPjmoniKqqrXob8faczZvL3wf6lEDvUEPiwBAq2/8DbX650u1wqqdiP2Jk4cSL2glakeqyO+PYXSynY8Onpe65/y/D6czcbWTYjHgF7AfzrJuOe3pCVpgThokcieW/UZyCzrAAW/NqI2J2+OpV2yFph5MO3DTWOK8oikiyWPtiNUa8d4MgJTSqpTGMaN3ETs5jFTGaSQAIAC1hA53aKpfdGEjnEtTqzS358dD/Y+63RiQB8+zh0H12zveLVC8JFj4i8N+ojjls+MAQ+aTxE93ZOVtJY2TTlxUZGTbcx0CaQnkfSWXrXUUb9q4K0kjRu4iYSSGAucwHIIos00ugcHsjSu4LpGVppTIIy69AMnQpVdji6DfJ3GQIfGget2xhZO97q0ItXLwgXPWLXeKM+k4x6XmfkvYd1hSqH8b37OOiSDJ0GQtswY3A0PM6YAHVwNT0vu4o5t3SkiCJmMcvtcrOYRRFFzPlFH0Pgk8Ybbwvf/wFW/t0p3NroTMITjf39fmLcIzjS+L1vkfv6sDLZSRAuekTkvVEfcTRntW542/g7kAYr/gqr/2Hk0JcXQVSS0VEc3WGckl/F1Pn5hBHGTGa6XW4mMwkjjKnv72Rvt3uNiDx/l7Fz95eGaDuvQ0kO3P4pDPqFIe79bzldhx51bguPC4Lgk4jInysD7zBmtSb/CrpeYQyEjv+zIbQ3vmGI78/mGR3FhOfZGzKcUS9u40hRBSmkkEACWWQxhSlkkUUCCaSQwpGiCkY99h57cwqgQx/jXkUHDdEe9pAxG3bYQ8b2HZ8Z0fuOz053TEOnyqIjgiBYiCdfF+YAbM/rYPO/DD98wt8gOslZd2aGYY2kvw0JVxnRfc/rDNHtPMjIYy8roGDDp4x6YTNH8vIBY3AVII00iijiUR4lhRRr+5HjxYxKGcm23w8iEoy0ybgRxpuC/Tj873eQcwsUHzHaaT9+us2STikIggsi8nVl0rgurL1vkbFt0QzDKjHPNevUVJUZKY4b/2n45GCVGY5c/xy/GtaJWQuOWZc2s2g+mBjI1IXFLDixwO3WvxrXl8hjq4zI3X7cEPjxfz6de7/sOWNgFyBvt9HZSFkDQRA8aHKRV0pdC7wMtAbe0VrPbup7nhV1ZdKYpYCryiC0KxTuOz2btazgdI33kU8adWzAEPiuVxhZL2admayVPDOkK9gyrRmvndsplv5+FD0HjmBpSg6jfvcBR4qrAJg50p9nbhsBeSHQPg5s+40MnqAImPIdrPo75KRDeHcIi4fIHs5VoVxqzQuCINDEIq+Uag28BowDsoENSqmvtNa7mvK+Z0VdaYbBkUYK47LZhs898W+n9235wL3Gu1morKrMGCA9sNSwb8B4Cxj5JM88+1dYkMHbc//F0hfvpefNTwPQc+81LL3Tn1H/hl9dl8wz13YwMnLM6yeNNwZ0g6OcA8JRcHC18TdyulG90myHIAiCC00dyQ8FMrTWBwCUUh8BNwEXjsifycOurRNwXfDDrBY56inDNjmw1MiscYnk6X8rRCfxDC/zsDpFZFzU6YJnx/fRs1dvtr1zNZE73oUDOyHucuMNodJhCH7ny9zbUpZvdCZVZcb9ksaf7mgEQRCcNHV2TRfgsMvvbOe2lkOd6ZSaGtHzwDsMwT2+z4jkzTo2ZlTf8zoiB1xr5NZ//wfj/HF/hJ/NIzIyGkY8bETnQ6ca1ktIFKx+BY5scm9TcLQh7m2CjfNdJ0EJgiA4aepI3pt/oN0OUGoqMBUgPj6+iZvTiLjWgvcPcl9+L/VNo5hYZRl0GWpE9XEjjP1m2mNUb2c5Ye1eh37cH93fLFxLF69/6/T6ra5RvYi7IAi10NQinw3EufyOBY64HqC1ngPMAUhOTnbrAC5oBt5hiDiqppVjevnf/8EQ+OP7jOyY2z/F6uP8A08f47rwtrdrpb5pCPzh9Ub0DkZHIKmSgiCcgaYW+Q1AklIqEcgBfg7c1sT3PD+YefImnqmYplib+e1mVs7Q+9xFHU6fU5toB0ca53h6/YIgCGegST15rXU18BDwHbAb+ERrvbMp79lseC70bYp2/BAjgo9Oct/uWrI4OPJ0WWDXssWueHr9nkhZYUEQvNDkefJa6/8CXlTJx3DLfDmHNVe9TbwyI3vzelf9zvjtLZKXssKCIHhBZrw2Fp4LiZyt4JqdRM/rIOFKdzvHFHCz0mTClRDtZeKW66cgCAIi8k3D2Qiua9Rvdgi1Cbi3DsBEatYIguAFEfmm4GwE19Nm8Wb1uF7PswMQBEGoAyk13Jy4Fjgzo3PPAVxBEIQGIJF8c+I6AUqW7BMEoQkQkW9OvAm6eOuCIDQiIvLNiQi6IAhNjHjygiAIPoyIvCAIgg8jIi8IguDDiMgLgiD4MCLygiAIPoyIvCAIgg8jIi8IguDDiMgLgiD4MCLygiAIPoyIvCAIgg8jIi8IguDDiMgLgiD4MCLygiAIPoyIvCAIgg/TIJFXSt2qlNqplDqllEr22PeUUipDKfWjUuqahjVTEARBOBcaWk9+BzAJeMt1o1KqL/BzoB/QGVislOqptT7ZwPsJgiAIZ0GDInmt9W6t9Y9edt0EfKS1rtBaZwIZwNCG3EsQBEE4e5rKk+8CHHb5ne3cVgOl1FSlVLpSKj0/P7+JmiMIgnBxcka7Rim1GOjoZdfTWusvazvNyzbt7UCt9RxgDkBycrLXYwRBEIRz44wir7Ueew7XzQbiXH7HAkfO4TqCIAhCA2gqu+Yr4OdKqQClVCKQBKxvonsJgiAItdDQFMqfKKWygeHAN0qp7wC01juBT4BdwP+AByWzRhAE4fzToBRKrfUXwBe17HsWeLYh1xcEQRAahsx4FQRB8GFE5AVBEHwYEXlBEAQfRkReEATBhxGRFwRB8GFE5AVBEHwYEXlBEAQfRkReEATBhxGRFwRB8GFE5AVBEHwYEXlBEAQfRkReEATBhxGRFwRB8GFE5M8TtnIbc3fMxVZua+6mCIJwESEif55YkLGAFze+yIKMBc3dFEEQLiIaVE9eqD+pPVLdPgVBEM4HEsk3Et7sGFu5jde3vM7rW14HDIFfkLGgQZaN2D6CIJwNEsk3EqYdAzCl/xRr2xtb3wAg0C8QwDrGFPzUHqmEtw13u5at3FZjn7nNUe2wrlnXNQRBEEBEvtHwZsek9kjFUe3wun3ennm8sfUNHNUOJvee7CbWtXUYL258kfsH3M9vB//WEnjP4wRBEFwRka8H3iJrb/um9J9i2SnmsQ8MfMDteG9i7CnW9ekwwtuGi88vCMIZaZDIK6WeB24AKoH9wBStdZFz31PAvcBJ4BGt9XcNa2rT4ynm3iwST5F2jcgfGPiAJdiOageBfoGWAM/bMw9HtYPi8mI2HtvIU0OfsiJyE/N7eNvwGvcJbxtOoF8gL258kUC/QKb0n2IJvVg2giDURkMj+e+Bp7TW1Uqp54CngCeVUn2BnwP9gM7AYqVUT631yQber0nxjKi9WSSeHC096vZpHuOodljXAqxOwmTGyhm8MvoVnl75NPddeh+Oagfv7XgPMPz7yb0nA0bnADC592SvkbtYNoIg1EWDRF5rvcjl51rgFuf3m4CPtNYVQKZSKgMYCqxpyP2aGtMScVQ7sJXb3ES1Nptmfd56AJYcXkLHLR2Z3HuyZdsE+gVydfjVRERGUFheyI7jO+gQ2IHVuavJL8jn/1b9H5klmew4vgNbhXu2zI7jO+gf1d9t4NZV3M37p8SlWG0UBEHwpDE9+XuAj53fu2CIvkm2c1sNlFJTgakA8fHxjdics8ebJVJbdGzaNL3CetG2dVtKKkt4Y+sblhjP2zOPb978hqe/epq0tDQA0vPSSY5J5omuT3DbQ7cx/BfD4SrILMlkWKdh9I7oTXF5MUuzl7IiZwX9o/pzd7+72VOwh5S4FLeoHZAIXhCEM3JGkVdKLQY6etn1tNb6S+cxTwPVwIfmaV6O196ur7WeA8wBSE5O9nrM+aS2wUxPv94cBP2x6EcAS6QLywt5+IeHWTRnEflf5gOQkpJC8qxk8IdVW1fxyfOfYC+wM/eluTzi/wgxo2PoHdGbSUmTmPbDNIoqiohvZ3R4gX6BrD26lm8zvwXg/gH3e/XxBUEQvHFGkddaj61rv1LqLmAiMEZrbYp0NhDnclgscORcG9kY1JUh44q3QU+oOcBq5r13Ce7CuIRxXN7xcmasnIGtwkbeF3mWwAPk5uayaPoiOv6iI0f/dZSKwgpr3z+e+wfRe6JZ+5O17CncQ2ZJJmEBYYQHhPPG1jes8QBz8Pe3g39rtV8ieEEQzkSDZrwqpa4FngRu1FrbXXZ9BfxcKRWglEoEkoD1DblXQ2ns2jGTe0/mqi5XkVOWQ0TbCP664a+Gr14GtmXu/noqqQQWBnLw7wcJLAwklVS3/bZlNqpLq8kqziI5JpmiiiK2Ht/K4A6DjfN7pDK59+RaB38FQRBqo6FlDV4F2gHfK6W2KKXeBNBa7wQ+AXYB/wMebO7MmtQeqTVEsq4SAZ77JiRO4KouVzEhcQJgRPzPXvmsdc3/d8X/I75dPNf1u44bnr+BdlHtjPuSyjSm8RIvkUACL/ES05hmCX1gZCBDZg7BL8SPo/ajaK0Z1mkYAPmOfN7Y+gYLMha4pUtmFmd6LaEg5Q4EQfCkodk1PerY9yzwbEOu35h4s2HqSj90zXd3VDv4/uD35JTm0D+qvzXBybymrdzGpmObuL7b9byx9Q2GxQ+j4+Mdccx2kFaUxk3cRAIJzGUuAFlkkUYa7aPb85u3f8P8ovnWfZVSzLh8Bs9veJ4VOSu4qstVbrn2b2x9g9VHVrM2d61bLr7rs0juvCAIJhf1jNe6Blkd1Q7uH3A/AP/c+c8a57p6/J759IXlhaztuJaRj45kycwlzGKWJfAAs5hFEUUsnLuQ2CGxrFu+jrC2Yfi18iM9L535++bTPaw7jmoH3cO6W+eZg70J7RMY0XmEWy6+67N4jh8IgnDxclGL/JkGWe8fcD+Te0/GUe2w8tYnJE7g9S2vszlvM2uPGlmirgJbVFHErDWzSKpKYsmrSwgjjJnMdLv+TGbyKI/yq1/9ipQ/p5CjcugW1o3uYd3Zmr+VxQcXk12aDRhplxFtI5jSf4o12JtRlEF423AmJE6wIvnankUQhIubi1rk60N423AeS37M+v1C+gtWZD+s0zBS4lKsWakAs9bMYvXW1WTOzqS6qJpUUkkggSyymMUsZjKTBBJIIYUFuQtY/NRixv5lLE8MecJKk8wuzbZSMl0nQU3uPZkdx3ewImcF6XnpVi6/J5N7T64xeUoQhIsTEXkvTEicwI7jO5iQOKFG6uX249ut405UnuChJQ9x6MQhwMhpryipsAQeYAELAEgjjSKKeJRHDYF3bs8/ms+Xj3/JfcPuA+DufndbZQ08/XRzsNfsVGoTcc+ovr7po4Ig+B4i8l5IO5zGipwVDOk4BHCfWXpJ1CVszNsIwM6CnQDEt4tndPxoCssL6R/fn93X72brh1ut6y1gAX5hfnS9uyv5/85nQcECt/sFXxnM/235P2wVNn47+Ld12i7eKlueCfHoBeHiRUTeC7WV+gW4p/89BPoF4qh2UF5dTlZxFjOGzSDtcJrVGVx191UcKT1iTYhqH92emMdiuKzfZZzqf4rl/7ecvNw8AAbcPoAuk7pw6MQhrupyFSlxKW6lil2pazERidIFQfCGiLwXPO0O1+/hbcOtRT5S4lJIO5xGWECYW733CYkT6P9Mf77p+A2bv9rML175Bd+UfkNReRE5fjmM/vNofpjxA3Fj43jnb++w6sgq67w/r/sza3PXsix7GUM7DnWzbbyVMa4tDdRWbrNsHdcBWkEQLi5E5M8SW7mNp1c+zYqcFWw4usH6fPbKZ92skAcGPsADbz5AwbMFfHz4Y9gKpzjF4A6DeXTwo7Rt35YNJzbw1ra3WJGzgt8O/i1ph9NYm2tk7GzM28jGvI1sPraZv179V7dFQmpLnXTFc+lBybwRhIsTEfmzZEHGAmuS0hNDngBgRc4Knl75NE8MeYK0w2lWhJ/aI5USvxI2520mNiSW7NJseoT1YNOxTWw4scG6Rv+o/jiqHUxInGBNvtqev52NxzayNnctCzIWWIuEuJYxrit1sralBwVBuLgQkXdSX2/bs8b8s1c+y++W/Y4VOSuoOlXF2ty1fLHvCzJLMtlwdAP2ajsb8zbSL7IfbVq14b5L7yO+fbzbNVzLG5tvA652i6dI1ycn/lwGaAVB8D1E5J1487vrM5AZ3jacQTGDWHt0LQntE8gryyOzJJPE9omsyFlBckwyYFgmOwt2sunYJgZ0GOAm0t4sF0+Rdl0kxHxLkIFWQRDOhIi8E29+d30W3YbTk48c1Q4ySzItG8bVunEVZ09cI/Pa3ijM+5rjALW1TxAEwRUReSfe/G5v1BZ1e/PKE0MTAWp81obroC7gNdpPiUthSMch1m9JoRQEoS4aWmrY5zAFuzbBPNP++lBbWWDXQd3afPiwgLAa5zRmnXxBEHwLieQbkbpKF9fnuLoWDq/t3NpSKAVBEEBEvlGpr+DWdpyr7VPbrFfPc6X6pCAIdSEi3wi4pjt6KyzmyZmEua43AhF1QRDOBhH5RsAsAAaNM7tULBhBEBoLEflGZFinYY0izBKtC4LQWIjINwKui3RIGqMgCBcSDUqhVEr9P6XUNqXUFqXUIqVUZ5d9TymlMpRSPyqlrml4Uy9cGiOtUhAEoSloaJ7881rrS7XWA4GFwB8AlFJ9gZ8D/YBrgdeVUq0beC9BEAThLGmQyGutS1x+BgPa+f0m4COtdYXWOhPIAIY25F6CIAjC2dNgT14p9SzwC6AYGOXc3AVY63JYtnObt/OnAlMB4uPjG9ocQRAEwYUzRvJKqcVKqR1e/m4C0Fo/rbWOAz4EHjJP83Ip7WUbWus5WutkrXVydHT0uT6HIAiC4IUzRvJa67H1vNZ/gG+AmRiRe5zLvljgyFm3ThAEQWgQDc2uSXL5eSOwx/n9K+DnSqkApVQikASsb8i9BEEQhLOnoZ78bKVUL+AUcBD4NYDWeqdS6hNgF1ANPKi1PtnAewmCIAhnidLaq1XeLCil8jE6i3MhCjjeiM1pTnzpWUCe50LHl57Hl54F6v88XbXWXgc1LyiRbwhKqXStdXJzt6Mx8KVnAXmeCx1feh5fehZonOeRRUMEQRB8GBF5QRAEH8aXRH5OczegEfGlZwF5ngsdX3oeX3oWaITn8RlPXhAEQaiJL0XygiAIggci8oIgCD6MT4i8UupxpZRWSkW5bGtx9ex9rT6/Uup5pdQe5zN9oZQKc9nXop5HKXWrUmqnUuqUUirZY1+LehYTpdS1zjZnKKWmN3d7zhal1HtKqWNKqR0u2yKUUt8rpfY5P1vMIg9KqTil1FKl1G7nf2vTnNsb9kxa6xb9h1Ej5zuMSVRRzm19ga1AAJAI7AdaN3db6/Es7V2+PwK82cKfZzzg5/z+HPBcS30eoA/QC0gDkl22t7hncba7tbOt3QB/5zP0be52neUzXA1cBuxw2fZXYLrz+3Tzv7mW8Ad0Ai5zfm8H7HX+99WgZ/KFSP4l4He4V7lskfXstY/V59daL9JaVzt/rsUoVAct8Hm01ru11j962dXinsXJUCBDa31Aa10JfITxLC0GrfVyoNBj803A+87v7wOp57NNDUFrnau13uT8fgLYjVGivUHP1KJFXil1I5Cjtd7qsasLcNjld6317C80lFLPKqUOA7fjXGmLFvw8LtwDfOv87gvPY9JSn6WltvtMxGitc8EQTaBDM7fnnFBKJQCDgHU08Jku+IW8lVKLgY5edj0NzMCwBGqc5mXbBZErWtfzaK2/1Fo/DTytlHoKoz7/TFrw8ziPeRqjUN2H5mlejm/256nPs3g7zcu2Zn+WetBS2+3zKKVCgM+B32itS5Ty9q+q/lzwIq9rqWevlLoEwwPd6vyHEAtsUkoN5QKuZ1/b83ihRdTnP9PzKKXuAiYCY7TTVOQCfZ6z+HfjygX5LPWgpbb7TOQppTpprXOVUp2AY83doLNBKdUGQ+A/1FrPd25u0DO1WLtGa71da91Ba52gtU7A+I/2Mq31UVpoPXtfq8+vlLoWeBK4UWttd9nVIp+nFlrqs2wAkpRSiUopf+DnGM/S0vkKuMv5/S6gtjewCw5lRKvvAru11i+67GrYMzX3iHIjjkxn4cyucf5+GiN74EdgQnO3r57P8DmwA9gGfA10aeHPk4Hh+25x/r3ZUp8H+AlGIFEB5AHftdRncWn3dRgZHPsxLKlmb9NZtn8ekAtUOf/d3AtEAkuAfc7PiOZu51k8z5UYltk2l/9nrmvoM0lZA0EQBB+mxdo1giAIwpkRkRcEQfBhROQFQRB8GBF5QRAEH0ZEXhAEwYcRkRcEQfBhROQFQRB8mP8Pc+azxJ/I8GAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids, data, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "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": "markdown",
   "metadata": {},
   "source": [
    "So here's the definition of the gaussian kernel, which you may remember from high school..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def gaussian(d, bw): return torch.exp(-0.5*((d/bw))**2) / (bw*math.sqrt(2*math.pi))\n",
    "\n",
    "x = torch.linspace(0,10,100)\n",
    "plt.plot(x, gaussian(x,2.5));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " This person at the science march certainly remembered!\n",
    "\n",
    "<img src=\"http://i.imgur.com/nijQLHw.jpg\" width=400>"
   ]
  },
  {
   "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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([0.0000, 1.4130, 3.2164, 2.8909, 4.5990, 3.1394, 3.9166, 5.3368])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = data.clone()\n",
    "x = data[0]\n",
    "dist = torch.sqrt(((x-X)**2).sum(1))\n",
    "dist[:8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([1.5958e-01, 1.3602e-01, 6.9749e-02,  ..., 3.5634e-09, 1.7959e-10,\n",
       "        2.7274e-16])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight = gaussian(dist, 2.5)\n",
    "weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([1500]), torch.Size([1500, 2]))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight.shape,X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[-2.1094e-01,  3.9275e+00],\n",
       "        [-2.9346e-01,  3.1928e+00],\n",
       "        [-1.0870e-01,  1.4929e+00],\n",
       "        ...,\n",
       "        [-5.0652e-08,  1.1388e-07],\n",
       "        [-2.6280e-09,  6.0299e-09],\n",
       "        [-4.7550e-15,  1.0221e-14]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(weight[:,None]*X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def meanshift(data):\n",
    "    X = data.clone()\n",
    "    for it in range(5):\n",
    "        for i, x in enumerate(X):\n",
    "            dist = torch.sqrt(((x-X)**2).sum(1))\n",
    "            weight = gaussian(dist, 2.5)\n",
    "            X[i] = (weight[:,None]*X).sum(0)/weight.sum()\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 567 ms, sys: 0 ns, total: 567 ms\n",
      "Wall time: 567 ms\n"
     ]
    }
   ],
   "source": [
    "%time X=meanshift(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that mean shift clustering has almost reproduced our original clustering. The one exception are the very close clusters, but if we really wanted to differentiate them we could lower the bandwidth.\n",
    "\n",
    "What is impressive is that this algorithm nearly reproduced the original clusters without telling it how many clusters there should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWrklEQVR4nO3db4xjV3nH8e9DQkOnlDhxhhAnUR3EGvIHF8Q0ouqLeiARS3fEmlSUINSuFmQrUoimiBdNFKm2VUWKhMpULSCwRZYtDYSoJJPtDH+SbD2KkApk0kZLNpt1VjA0K0dkYpiEdiTohqcvxnbsnZnMztoer8/9fSTL955j+54zkX65e/z4XnN3REQkTK8b9gBERGRwFPIiIgFTyIuIBEwhLyISMIW8iEjAzh/2ADpdcsklnkwmhz0MEZGR8sQTT7zo7uMb9Z1TIZ9MJllcXBz2MERERoqZ/WyzPi3XiIgETCEvO6bRaGyrXUR6p5CXHVEsFkmn09Rqta72Wq1GOp2mWCwOZ2AigVPIy8AVi0VKpRL1ep3Jycl20NdqNSYnJ6nX65RKJQW9yAAo5GWgWgHf0gr6+fn5dsC3KOhF+k8hLwPTaDSoVCpdbVmyrNZXmZqaYrW+SpZsV3+lUtEavUgfKeRlYOLxONVqlUQiAawF/DTTzDBDkiQzzDDNdDvoE4kE1WqVeDw+xFGLhEUhLwOVSqXaQb/AAksskSTJAQ6QJMkSSyyw0A74VCo17CGLBEUhLwOXSqUol8ussEKJUldfiRIrrFAulxXwIgOgkJeBq9Vq5PN5YsQoUOjqK1AgRox8Pr+uvFJEeqeQl4HqLJPMkGkv0exnf3vpJkNmXXmliPSHnUu3/5uYmHBduyYcjUaDdDrdVSaZJcsCC6ywQowYGTLMMtvuTyQSHDlyRF++imyDmT3h7hMb9elMXgYmHo+Ty+W62maZZSwxxtzcHGOJsa6AB8jlcgp4kT5SyMtAFYtFCoVX1+FbVTR79uzpKq8EKBQK+jGUSJ+dU5caljC1grtSqXSVSbbKKycnJ8nlcgp4kQHQmrzsmEajseFSzGbtInJmtCYv54TNglwBLzI4CnkRkYAp5EVEAqaQl7459ctf0vjKVzj1y18Oeygi0qTqGumblx54gK8+8DV4+EEAPvPNuSGPSER0Ji99c+FNN8Eb39je//uPTg1xNCICCnnpo/Mvumhd2z/u/+gQRiIiLQp56avxq97Wtf9/q/87pJGICCjkpc/+6u5/AKyrbfXll4YyFhFRyMsAfOab/8Yb3hRr7x9deHR4gxGJuJ6ra8zsDcBjwAXNz/tXdy+Y2cXAN4EksAT8hburti4ibq38C6svv8TRhUe5NnPDsIcjEln9OJP/NfA+d/9D4F3AbjN7L3A7cNjddwGHm/sSIWNvupA/+tCfM/amC4c9FJHI6jnkfc3/NHdf33w4sBc42Gw/CGR7PZZIS6PR2Fa7SFT1ZU3ezM4zsyeBF4BH3P2HwKXu/jxA8/nNm7w3b2aLZra4vLzcj+FI4IrFIul0et2tAmu1Gul0WpcsFunQl5B391fc/V3AFcD1ZnbdNt5bdvcJd58YHx/vx3AkYMVikVKptO6esJ33ki2VSgp6kaa+Vte4+wqwAOwGfm5mlwE0n1/o57FkNK2+/BKPH/rWWZVVtgK+pRX08/Pz7YBvUdCLrOk55M1s3Mxize3fBW4AngEOAfuaL9sHPNTrsWT0HV14lMfuPcB3vzizraBvNBpUKpWutixZVuurTE1NsVpfJXva1z6VSkVr9BJ5/TiTvwyomtkR4HHW1uTngLuBG83sWeDG5r5E3LWZG7jq3RP89L8W+fIt+zj2/cfO6H3xeLzrnrBZskwzzQwzJEkywwzTTLeDvnUvWd2QRKJOt/+THbf68kt8KfdxHHjdK6+w+6klrn7m2Bm9t7X2vlpfbQd8yxJLfJpPM5YY67qXrEjodPs/OaeMvelC0j95nte98grv/NnaVzXH3nH1Gb03lUpRLpdZYYUSpa6+EiVWWKFcLivg+0BlqmFQyMtQvPUd17L7qSUu/9Xqtt5Xq9XI5/PEiFGg0NVXoECMGPl8fl15pWyPylTDoZCXoXjr1/552+/pLJPMkCFJkiWW2M9+llgiSZIMmXXllbI9KlMNjLufM4/3vOc9LtHy9Nvf0X68lhdffNETiYSz9mtqBzxL1mPEHPAYMc+S7epPJBL+4osv7tBMwlAoFLr+hq2/49zc3Lq/P+CFQmHYQxZ3BxZ9k1zVmbwM1dXPHGs/Xks8HieXy3W1zTLLWGKMubk5xhJjzDLb1Z/L5VRdsw0qUw3UZuk/jIfO5GUrnWeaiUTCjx8/7u7ux48f7zrT1Bnm2en8O2bJepWqH+CAJ0n6AQ54lWr7X0ydf38ZLl7jTF4llDJyisUilUplXZlka804l8tpvbgHKlMdPa9VQqmQl5HUaDQ2XIrZrF22Z35+nqmpKZIkOcCBdnvrS+65uTn27NkzxBFKJ9XJS3A2C3IFfO9UphoWhbyItKlMNTxarhERYG2pK51Od13NM0uWBRZYYYUYMTJkuqqYEokER44c0b+ghkzLNSKyJZWphqnnG3mLyOhJ3j7f3l66+9UvUFtVSa3r9reu5plKpahWq13X7S8UCqpiGgEKeRHp0gru08tUO4NeZaqjQyEvIusUi0Vuu+22dUsxqVRKa/AjRiEvEkGdSzSbUZlqGPTFq4hIwBTyIiIBU8iLiARMIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIjIkGx2f9x+3je355A3syvNrGpmx8zsqJlNN9svNrNHzOzZ5vNFvQ9XRCQMxWKRdDq97pr8tVqNdDrdt2sD9eNM/hTwGXe/GngvcKuZXQPcDhx2913A4ea+iEjkFYtFSqXSupuvdN60pVQq9SXoew55d3/e3f+zuf0r4BhwObAXONh82UEg2+uxRERGXSvgW1pBPz8/33UpZ6AvQd/XNXkzSwLvBn4IXOruz8Pa/wiAN2/ynryZLZrZ4vLycj+HIyJyTmk0GlQqla62LFlW66tMTU2xWl8le9r5cKVS6WmNvm8hb2ZvBL4F/LW7v3ym73P3srtPuPvE+Ph4v4YjInLOicfjVKtVEokEsBbw00wzwwxJkswwwzTT7aBv3bSllyt/9iXkzez1rAX8ve7+QLP552Z2WbP/MuCFfhxLRGSUtW6+kkgkWGChfYP0Axxo3zh9gYWuu3L1oh/VNQZ8BTjm7p/r6DoE7Gtu7wMe6vVYIiIhSKVSlMtlVlihRKmrr0SJFVYol8s9Bzz050z+T4C/BN5nZk82H38G3A3caGbPAjc290VEIq9Wq5HP54kRo0Chq69AgRgx8vn8uvLKs9GP6prvu7u5e9rd39V8fNvdG+7+fnff1Xz+Rc+jFREZcZ1lkhky7SWa/exvL91kyKwrrzxb5u59GnrvJiYmfHFxcdjDEBEZiEajQTqd7iqTzJJlgQVWWCFGjAwZZplt9ycSiS3vq2tmT7j7xEZ9uqyBiMgOicfj5HK5rrZZZhlLjDE3N8dYYqwr4AFyudzwq2tEROTMFItFCoVX1+FbVTR79uzpKq8EKBQKPf8Y6vye3i0iIut84ZZ/b2/f+qX3retvBXelUukqk2yVV05OTpLL5fpyWQOFvIjIEBSLRW677bZ1SzGpVGrLNfjt0HKNiMiQbBbk/Qp40Jm8iEjfbbREMyw6kxcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYH0JeTO7x8xeMLOnOtouNrNHzOzZ5vNF/TiWiIicuX6dyX8V2H1a2+3AYXffBRxu7ouIyA7qS8i7+2PAL05r3gscbG4fBLL9OJaIiJy5Qa7JX+ruzwM0n9+80YvMLG9mi2a2uLy8PMDhiIhEz9C/eHX3srtPuPvE+Pj4sIcjIhKUQYb8z83sMoDm8wsDPJaIiGxgkCF/CNjX3N4HPDTAY4mIyAb6VUL5DeA/gLeb2Ukz+yRwN3CjmT0L3NjcFxGRHXR+Pz7E3T+2Sdf7+/H5IiJydob+xauIiAyOQl5EJGAKeRGRgCnkRUQCppBvajQa22oXERkFCnmgWCySTqep1Wpd7bVajXQ6TbFYHM7ARER6FPmQLxaLlEol6vU6k5OT7aCv1WpMTk5Sr9cplUoKehEZSZEO+VbAt7SCfn5+vh3wLQp6ERlFkQ35RqNBpVLpasuSZbW+ytTUFKv1VbKnXR25UqlojV5ERkpkQz4ej1OtVkkkEsBawE8zzQwzJEkywwzTTLeDPpFIUK1WicfjQxy1iMj2RDbkAVKpVDvoF1hgiSWSJDnAAZIkWWKJBRbaAZ9KpYY9ZBGRbYl0yMNa0JfLZVZYoUSpq69EiRVWKJfLCngRGUmRD/larUY+nydGjAKFrr4CBWLEyOfz68orRURGQaRDvrNMMkOmvUSzn/3tpZsMmXXllSIio8LcfdhjaJuYmPDFxcUdOVaj0SCdTneVSWbJssACK6wQI0aGDLPMtvsTiQRHjhzRl68ick4xsyfcfWKjvsieycfjcXK5XFfbLLOMJcaYm5tjLDHWFfAAuVxOAS8iI6UvNw05ZxUv7Nh+aX1388dNrR9EdVbRVKvVrh9EFQoF/RhKREZO2CF/BlrBXalUusokO4M+l8sp4EVkJIW9Jr/FmXynRqOx4VLMZu0iIueK11qTD/tMfotg77RZkCvgRWSURfaLVxGRKFDIi4gETCEvIhIwhbyISMAU8iIiARt4yJvZbjM7bmYnzOz2QR9PREReNdCQN7PzgC8AHwSuAT5mZtcM8pgiIvKqQZ/JXw+ccPefuPtvgPuAvQM+poiINA065C8HnuvYP9lsazOzvJktmtni8vLygIcjIhItgw5526Ct6zoK7l529wl3nxgfHx/wcEREomXQIX8SuLJj/wqgvslrRUSkzwYd8o8Du8zsKjP7HeBm4NCAjykiIk0DvUCZu58ys08B3wPOA+5x96ODPKaIiLxq4FehdPdvA98e9HFERGQ9/eJVRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgPYW8mX3EzI6a2W/NbOK0vjvM7ISZHTezD/Q2TBERORvn9/j+p4CbgC93NprZNcDNwLVAAnjUzFLu/kqPxxMRkW3o6Uze3Y+5+/ENuvYC97n7r939p8AJ4PpejiUiIts3qDX5y4HnOvZPNtvWMbO8mS2a2eLy8vKAhiMiEk1bLteY2aPAWzboutPdH9rsbRu0+UYvdPcyUAaYmJjY8DUiInJ2tgx5d7/hLD73JHBlx/4VQP0sPkdERHowqOWaQ8DNZnaBmV0F7AJ+NKBjiYjIJnotofywmZ0E/hiYN7PvAbj7UeB+4Gngu8CtqqwREdl5PZVQuvuDwIOb9N0F3NXL54uISG/0i1cRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYBFLuQbjca22kVERlmkQr5YLJJOp6nVal3ttVqNdDpNsVgczsBERAYkMiFfLBYplUrU63UmJyfbQV+r1ZicnKRer1MqlRT0IhKUSIR8K+BbWkE/Pz/fDvgWBb2IhKSnkDezz5rZM2Z2xMweNLNYR98dZnbCzI6b2Qd6HulZajQaVCqVrrYsWVbrq0xNTbFaXyVLtqu/UqlojV5EgtDrmfwjwHXungZqwB0AZnYNcDNwLbAb+KKZndfjsc5KPB6nWq2SSCSAtYCfZpoZZkiSZIYZppluB30ikaBarRKPx4cxXBGRvuop5N39YXc/1dz9AXBFc3svcJ+7/9rdfwqcAK7v5Vi9SKVS7aBfYIEllkiS5AAHSJJkiSUWWGgHfCqVGtZQRUT6qp9r8p8AvtPcvhx4rqPvZLNtHTPLm9mimS0uLy/3cTjdUqkU5XKZFVYoUerqK1FihRXK5bICXkSCsmXIm9mjZvbUBo+9Ha+5EzgF3Ntq2uCjfKPPd/eyu0+4+8T4+PjZzOGM1Go18vk8MWIUKHT1FSgQI0Y+n19XXikiMsq2DHl3v8Hdr9vg8RCAme0DpoCPu3sryE8CV3Z8zBVAnSHpLJPMkGkv0exnf3vpJkNmXXmliMios1dz+SzebLYb+Bzwp+6+3NF+LfB11tbhE8BhYJe7v/JanzcxMeGLi4tnPZ6NNBoN0ul0V5lkliwLLLDCCjFiZMgwy2y7P5FIcOTIEX35KiIjwcyecPeJjfp6XZP/PPD7wCNm9qSZfQnA3Y8C9wNPA98Fbt0q4AclHo+Ty+W62maZZSwxxtzcHGOJsa6AB8jlcgp4EQnC+b282d3f9hp9dwF39fL5/dL6cVPrB1GdVTTVarXrB1GFQkE/hhKRYPQU8ueSdx58Z3v7x/t+vK6/FdyVSqWrTLIz6HO5nAJeRILS05p8v/WyJr9VyLc0Go0Nl2I2axcROdcNck1+5GwW5Ap4EQlRMMs1r3X2LiISVZE7kxcRiRKFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwBTyIiIBO6d+8Wpmy8DPdvCQlwAv7uDxziVRnjtEe/5RnjuEOf8/cPcNb8hxToX8TjOzxc1+Chy6KM8doj3/KM8dojd/LdeIiARMIS8iErCoh3x52AMYoijPHaI9/yjPHSI2/0ivyYuIhC7qZ/IiIkFTyIuIBCxyIW9mf2dmR5o3Hn/YzBIdfXeY2QkzO25mHxjmOAfFzD5rZs80/wYPmlmsoy/o+ZvZR8zsqJn91swmTusLeu4tZra7OccTZnb7sMczaGZ2j5m9YGZPdbRdbGaPmNmzzeeLhjnGQYtcyAOfdfe0u78LmAP+FsDMrgFuBq4FdgNfNLPzhjbKwXkEuM7d00ANuAMiM/+ngJuAxzobIzJ3mnP6AvBB4BrgY825h+yrrP037XQ7cNjddwGHm/vBilzIu/vLHbu/B7S+ed4L3Ofuv3b3nwIngOt3enyD5u4Pu/up5u4PgCua28HP392PufvxDbqCn3vT9cAJd/+Ju/8GuI+1uQfL3R8DfnFa817gYHP7IJDdyTHttMiFPICZ3WVmzwEfp3kmD1wOPNfxspPNtpB9AvhOczuK82+JytyjMs+tXOruzwM0n9885PEMVDD3eO1kZo8Cb9mg6053f8jd7wTuNLM7gE8BBcA2eP1I1pduNf/ma+4ETgH3tt62wetHbv5nMveN3rZB28jN/QxEZZ7SIciQd/cbzvClXwfmWQv5k8CVHX1XAPU+D21HbDV/M9sHTAHv91d/KBHE/Lfx375TEHM/A1GZ51Z+bmaXufvzZnYZ8MKwBzRIkVuuMbNdHbsfAp5pbh8CbjazC8zsKmAX8KOdHt+gmdlu4G+AD7n7akdXJOa/iajM/XFgl5ldZWa/w9qXzYeGPKZhOATsa27vAzb7F14QgjyT38LdZvZ24LesXdb4FgB3P2pm9wNPs7aMcau7vzK8YQ7M54ELgEfMDOAH7n5LFOZvZh8G/gkYB+bN7El3/0AU5g7g7qfM7FPA94DzgHvc/eiQhzVQZvYNIANcYmYnWftX+93A/Wb2SeC/gY8Mb4SDp8saiIgELHLLNSIiUaKQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRg/w/7CCvH/YdeUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the computation is happening in the <tt>for</tt> loop, which isn't accelerated by pytorch. Each iteration launches a new cuda kernel, which takes time and slows the algorithm down as a whole. Furthermore, each iteration doesn't have enough processing to do to fill up all of the threads of the GPU. But at least the results are correct...\n",
    "\n",
    "We should be able to accelerate this algorithm with a GPU."
   ]
  },
  {
   "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": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dist_b(a,b): return torch.sqrt(((a[None]-b[:,None])**2).sum(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[0.6161, 0.7434, 0.7351, 0.9002, 0.5875, 0.5845, 0.2929, 0.1938],\n",
       "        [1.1132, 0.2402, 0.9845, 0.4507, 1.0699, 0.6556, 0.6886, 0.7938],\n",
       "        [0.0261, 0.9847, 0.3508, 1.0109, 0.0595, 0.5418, 0.4450, 0.4208],\n",
       "        [0.4530, 0.8696, 0.1635, 0.7858, 0.4366, 0.4354, 0.6171, 0.7125],\n",
       "        [0.6937, 0.4877, 0.4028, 0.3711, 0.6562, 0.2454, 0.5457, 0.6985]])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X=torch.rand(8,2)\n",
    "x=torch.rand(5,2)\n",
    "dist_b(X, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[1.9947e-01, 1.5541e-01, 5.4735e-02,  ..., 2.2132e-13, 2.0775e-15,\n",
       "         1.6831e-24],\n",
       "        [1.5541e-01, 1.9947e-01, 1.1173e-01,  ..., 3.1427e-13, 2.0317e-15,\n",
       "         9.7350e-25],\n",
       "        [5.4735e-02, 1.1173e-01, 1.9947e-01,  ..., 3.5939e-16, 9.4595e-19,\n",
       "         3.9723e-29],\n",
       "        [7.0175e-02, 3.1377e-02, 1.8560e-02,  ..., 5.2249e-18, 3.3886e-20,\n",
       "         3.0925e-30],\n",
       "        [1.4180e-02, 3.1085e-03, 1.1476e-03,  ..., 8.5815e-20, 8.2003e-22,\n",
       "         9.4978e-32]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bs=5\n",
    "X = data.clone()\n",
    "x = X[:bs]\n",
    "weight = gaussian(dist_b(X, x), 2)\n",
    "weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([5, 1500]), torch.Size([1500, 2]))"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weight.shape,X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([5, 2])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num = (weight[...,None]*X[None]).sum(1)\n",
    "num.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([5, 1])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "div = weight.sum(1, keepdim=True)\n",
    "div.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[-0.5274, 24.3688],\n",
       "        [-1.0042, 23.6233],\n",
       "        [-0.7959, 22.5723],\n",
       "        [ 1.0092, 24.4519],\n",
       "        [ 1.8587, 25.1916]])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num/div"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fastcore.all import chunked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "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)\n",
    "            num = (weight[...,None]*X[None]).sum(1)\n",
    "            div = weight.sum(1, keepdim=True)\n",
    "            X[s] = num/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": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = meanshift(data).cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.25 ms ± 290 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit -n 1 X = meanshift(data).cpu()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVrElEQVR4nO3db4hj13nH8e8TO0k7balseexYXlMZugq1HdUBxcT0RTWxjTfZIasETB1CWTZBIuAs05AXtTFUEsVgCMm8aBJSiXqzL5K4Lm3Gy4zzZ72dIRTyx7PFnXhjr7wkk3qRsScK07QMONh5+mIkWdrReHdG0mh07u8Dg+49R9I9x4afrx+de6+5OyIiEqZ3jHoAIiIyPAp5EZGAKeRFRAKmkBcRCZhCXkQkYFePegCdrrvuOk8mk6MehojIWDl79uyv3H2yV9++CvlkMsny8vKohyEiMlbM7Jfb9alcIyISMIW87JlGo7GjdhHpn0Je9kSpVCKdTlOr1braa7Ua6XSaUqk0moGJBE4hL0NXKpUol8vU63WmpqbaQV+r1ZiamqJer1MulxX0IkOgkJehagV8SyvoFxYW2gHfoqAXGTyFvAxNo9GgWq12teXIsVHfYHp6mo36BjlyXf3ValU1epEBUsjL0MTjcRYXF0kkEsBmwM8wwyyzJEkyyywzzLSDPpFIsLi4SDweH+GoRcKikJehSqVS7aBfYolVVkmS5AQnSJJklVWWWGoHfCqVGvWQRYKikJehS6VSVCoV1lmnTLmrr0yZddapVCoKeJEhUMjL0NVqNQqFAjFiFCl29RUpEiNGoVDYsrxSRPqnkJeh6lwmmSXbLtEc41i7dJMlu2V5pYgMhu2nx/9lMhnXvWvC0Wg0SKfTXcskc+RYYol11okRI0uWOeba/YlEgpWVFf34KrIDZnbW3TO9+nQmL0MTj8fJ5/NdbXPMMZGYYH5+nonERFfAA+TzeQW8yAAp5GWoSqUSxeJbdfjWKprDhw93La8EKBaLuhhKZMD21a2GZfx88a+m29uf/+f5nu9pBXe1Wu1aJtlaXjk1NUU+n1fAiwyBQl72RKlU4vjx41tKMalUSjV4kSFSuUb2zHZBroAXGR6dyUtftivRiMj+oDN5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYH2HvJn9npn9xMz+y8zOmVm52X6tmZ02s5ear9f0P1wREdmJQZzJvw58yN3/HLgDOGRmHwQeAs64+0HgTHNfRET2UN8h75v+r7n7zuafA0eAk832k0Cu32OJiMjODKQmb2ZXmdlzwGvAaXf/MXCDu78C0Hy9fpvPFsxs2cyW19bWBjEciYBGo7GjdpGoGkjIu/ub7n4HcAC408xu38FnK+6ecffM5OTkIIYjgSuVSqTT6S3Pg63VaqTTad2XXqTDQFfXuPs6sAQcAl41sxsBmq+vDfJYEk2lUolyubzlwd+dDwwvl8sKepGmQayumTSzWHP794F7gBeBU8DR5tuOAk/1eyyJtlbAt7SCfmFhoR3wLQp6kU2DOJO/EVg0sxXgWTZr8vPAY8C9ZvYScG9zX2RXGo0G1Wq1qy1Hjo36BtPT02zUN8hd8tt+tVpVjV4ir++Hhrj7CvD+Hu0N4O5+v18ENp8e1XoebL1eJ0eOGWY4whHKlClSJEkSgDnm2g8M11OnJOp0xauMjdaDvxOJBEssscoqSZKc4ARJkqyyyhJL7YBvPTBcdkcrmMKgkJexkkqlqFQqrLNOmXJXX5ky66xTqVQU8H3SCqZwKORlrNRqNQqFAjFiFCl29RUpEiNGoVDYEk5y5bSCKSwKeRkbnSGTJdsu0RzjWLt0kyW7JZzkymkFU3jM3Uc9hrZMJuPLy8ujHobsQ41Gg3Q63RUyOXIsscQ668SIkSXLHHPt/kQiwcrKin58vUL6Zzy+zOysu2d69elMXsZCPB4nn893tc0xx0Rigvn5eSYSE13hA5DP5xU+O9BawZRIJADaK5hmmSVJkllmmWGmvVRVK5jGg0Je9pX5uz7AF+//CPN3fWBLX6lUolh8qw7fCpnDhw93hRNAsVhUKWEXtIIpPAp52VfOH5iEd7xj87WHVtBfGjKd4aSA749WMIVFNXnZV+bv+gDnD0zy3otrTP/w2W3f12g0epYJtmuXK9f6gXujvtEu1bSsssrn+BwTiQmdye8jqsnL2Jj+4bN8/l+eftuAB7YNcgV8f7SCKTw6kxcRQKtrxpnO5EXksrSCKUx936BMRMZP8qGF9vbqY4fb260frFsXRHX+wN15gzjQCqZxoZAXkS6t4K5Wqz1XME1NTZHP5xXwY0IhLyJblEoljh8/vqUUk0qlVIMfMwp5kQjqLNFsRyuYwqAfXkVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVERqTRaOyofTf6Dnkzu9nMFs3sBTM7Z2YzzfZrzey0mb3UfL2m/+GKiIShVCqRTqe33JO/VquRTqcHdm+gQZzJvwF83t3/DPgg8KCZ3Qo8BJxx94PAmea+iEjklUolyuXyloevdD60pVwuDyTo+w55d3/F3f+zuf2/wAvATcAR4GTzbSeh+Yh3EZEIawV8SyvoFxYWum7lDAwk6AdakzezJPB+4MfADe7+Cmz+hwC4fpvPFMxs2cyW19bWBjkcEZF9pdFoUK1Wu9py5NiobzA9Pc1GfYPcJefD1Wq1rxr9wELezP4Q+Ffgb9z9N1f6OXevuHvG3TOTk5ODGo6IyL4Tj8dZXFwkkUgAmwE/w0z7gemzzDLDTDvoWw9t6efOnwMJeTN7J5sB/w13/7dm86tmdmOz/0bgtUEcS0RknLUevpJIJFhiqf2A9BOcaD84fYmlrqdy9WMQq2sM+CfgBXf/UkfXKeBoc/so8FS/xxIRCUEqlaJSqbDOOmXKXX1lyqyzTqVS6TvgYTBn8n8B/DXwITN7rvn3EeAx4F4zewm4t7kvIhJ5tVqNQqFAjBhFil19RYrEiFEoFLYsr9yNQayu+Q93N3dPu/sdzb+n3b3h7ne7+8Hm66/7Hq2IyJjrXCaZJdsu0RzjWLt0kyW7ZXnlbpm7D2jo/ctkMr68vDzqYYiIDEWj0SCdTnctk8yRY4kl1lknRowsWeaYa/cnEonLPlfXzM66e6ZXn25rICKyR+LxOPl8vqttjjkmEhPMz88zkZjoCniAfD4/+tU1IiJyZUqlEsXiW3X41iqaw4cPdy2vBCgWi31fDHV1X58WEZEuX/nMv7e3H/zah3q+pxXc1Wq1a5lka3nl1NQU+Xx+ILc1UMiLiIxAqVTi+PHjW0oxqVTqsjX4nVC5RkRkRLYL8kEFPOhMXkRkoLYr0YyKzuRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCdhAQt7MHjez18zs+Y62a83stJm91Hy9ZhDHEhGRKzeoM/mvA4cuaXsIOOPuB4EzzX0REdlDAwl5d/8B8OtLmo8AJ5vbJ4HcII4lIiJXbpg1+Rvc/RWA5uv1vd5kZgUzWzaz5bW1tSEOR0Qkekb+w6u7V9w94+6ZycnJUQ9HRCQowwz5V83sRoDm62tDPJaIiPQwzJA/BRxtbh8FnhrisUREpIdBLaH8FvBD4L1mdtHMPg08BtxrZi8B9zb3RURkD109iC9x909s03X3IL5fRER2Z+Q/vIqIyPAo5EVEAqaQFxEJmEJeRCRgCvmmRqOxo3YRkXGgkAdKpRLpdJpardbVXqvVSKfTlEql0QxMRKRPkQ/5UqlEuVymXq8zNTXVDvparcbU1BT1ep1yuaygF5GxFOmQbwV8SyvoFxYW2gHfoqAXkXEU2ZBvNBpUq9Wuthw5NuobTE9Ps1HfIHfJ3ZGr1apq9CIyViIb8vF4nMXFRRKJBLAZ8DPMMMssSZLMMssMM+2gTyQSLC4uEo/HRzhqEZGdiWzIA6RSqXbQL7HEKqskSXKCEyRJssoqSyy1Az6VSo16yCIiOxLpkIfNoK9UKqyzTplyV1+ZMuusU6lUFPAiMpYiH/K1Wo1CoUCMGEWKXX1FisSIUSgUtiyvFBEZB5EO+c5lklmy7RLNMY61SzdZsluWV4qIjAtz91GPoS2Tyfjy8vKeHKvRaJBOp7uWSebIscQS66wTI0aWLHPMtfsTiQQrKyv68VVE9hUzO+vumV59kT2Tj8fj5PP5rrY55phITDA/P89EYqIr4AHy+bwCXkTGykAeGrKvlf64Y/t/uruaFze1LojqXEWzuLjYdUFUsVjUxVAiMnbCD/nLaAV3tVrtWibZGfT5fF4BLyJjKfya/NucyXdqNBo9SzHbtYuI7BdvV5MP/0z+bYK903ZBroAXkXEW2R9eRUSiQCEvIhIwhbyISMAU8iIiAVPIi4gEbOghb2aHzOy8mV0ws4eGfTwREXnLUEPezK4CvgJ8GLgV+ISZ3TrMY4qIyFuGfSZ/J3DB3X/u7r8FngCODPmYIiLSNOyQvwl4uWP/YrOtzcwKZrZsZstra2tDHo6ISLQMO+StR1vXfRTcveLuGXfPTE5ODnk4IiLRMuyQvwjc3LF/AKhv814RERmwYYf8s8BBM7vFzN4FPACcGvIxRUSkaag3KHP3N8zss8D3gKuAx9393DCPKSIibxn6XSjd/Wng6WEfR0REttIVryIiAVPIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhIwhbyISMAU8iIiAesr5M3sfjM7Z2a/M7PMJX0Pm9kFMztvZvf1N0wREdmNq/v8/PPAx4F/7Gw0s1uBB4DbgATwjJml3P3NPo8nIiI70NeZvLu/4O7ne3QdAZ5w99fd/RfABeDOfo4lIiI7N6ya/E3Ayx37F5ttW5hZwcyWzWx5bW1tSMMREYmmy5ZrzOwZ4D09uh5x96e2+1iPNu/1RnevABWATCbT8z0iIrI7lw15d79nF997Ebi5Y/8AUN/F94iISB+GVa45BTxgZu82s1uAg8BPhnQsERHZRr9LKD9mZheBu4AFM/segLufA54EfgZ8F3hQK2tERPZeX0so3f3bwLe36XsUeLSf7xcRkf7oilcRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYBFLuQbjcaO2kVExlmkQr5UKpFOp6nVal3ttVqNdDpNqVQazcBERIYkMiFfKpUol8vU63WmpqbaQV+r1ZiamqJer1MulxX0IhKUSIR8K+BbWkG/sLDQDvgWBb2IhKSvkDezL5jZi2a2YmbfNrNYR9/DZnbBzM6b2X19j3SXGo0G1Wq1qy1Hjo36BtPT02zUN8iR6+qvVquq0YtIEPo9kz8N3O7uaaAGPAxgZrcCDwC3AYeAr5rZVX0ea1fi8TiLi4skEglgM+BnmGGWWZIkmWWWGWbaQZ9IJFhcXCQej49iuCIiA9VXyLv79939jebuj4ADze0jwBPu/rq7/wK4ANzZz7H6kUql2kG/xBKrrJIkyQlOkCTJKqsssdQO+FQqNaqhiogM1CBr8p8CvtPcvgl4uaPvYrNtCzMrmNmymS2vra0NcDjdUqkUlUqFddYpU+7qK1NmnXUqlYoCXkSCctmQN7NnzOz5Hn9HOt7zCPAG8I1WU4+v8l7f7+4Vd8+4e2ZycnI3c7gitVqNQqFAjBhFil19RYrEiFEoFLYsrxQRGWeXDXl3v8fdb+/x9xSAmR0FpoFPunsryC8CN3d8zQGgzoh0LpPMkm2XaI5xrF26yZLdsrxSRGTc2Vu5vIsPmx0CvgT8pbuvdbTfBnyTzTp8AjgDHHT3N9/u+zKZjC8vL+96PL00Gg3S6XTXMskcOZZYYp11YsTIkmWOuXZ/IpFgZWVFP76KyFgws7PununV129N/svAHwGnzew5M/sagLufA54EfgZ8F3jwcgE/LPF4nHw+39U2xxwTiQnm5+eZSEx0BTxAPp9XwItIEK7u58Pu/qdv0/co8Gg/3z8orYubWhdEda6iWVxc7Logqlgs6mIoEQlGXyG/n7zv5Pva2z89+tMt/a3grlarXcskO4M+n88r4EUkKH3V5Aetn5r85UK+pdFo9CzFbNcuIrLfDbMmP3a2C3IFvIiEKJhyzdudvYuIRFXkzuRFRKJEIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwPbVFa9mtgb8cg8PeR3wqz083n4S5blDtOcf5blDmPP/E3fv+UCOfRXye83Mlre7FDh0UZ47RHv+UZ47RG/+KteIiARMIS8iErCoh3xl1AMYoSjPHaI9/yjPHSI2/0jX5EVEQhf1M3kRkaAp5EVEAha5kDezvzezleaDx79vZomOvofN7IKZnTez+0Y5zmExsy+Y2YvNfwbfNrNYR1/Q8zez+83snJn9zswyl/QFPfcWMzvUnOMFM3to1OMZNjN73MxeM7PnO9quNbPTZvZS8/WaUY5x2CIX8sAX3D3t7ncA88DfAZjZrcADwG3AIeCrZnbVyEY5PKeB2909DdSAhyEy838e+Djwg87GiMyd5py+AnwYuBX4RHPuIfs6m/9OOz0EnHH3g8CZ5n6wIhfy7v6bjt0/AFq/PB8BnnD31939F8AF4M69Ht+wufv33f2N5u6PgAPN7eDn7+4vuPv5Hl3Bz73pTuCCu//c3X8LPMHm3IPl7j8Afn1J8xHgZHP7JJDbyzHttciFPICZPWpmLwOfpHkmD9wEvNzxtovNtpB9CvhOczuK82+JytyjMs/LucHdXwFovl4/4vEMVTDPeO1kZs8A7+nR9Yi7P+XujwCPmNnDwGeBImA93j+W60svN//mex4B3gC+0fpYj/eP3fyvZO69PtajbezmfgWiMk/pEGTIu/s9V/jWbwILbIb8ReDmjr4DQH3AQ9sTl5u/mR0FpoG7/a0LJYKY/w7+3XcKYu5XICrzvJxXzexGd3/FzG4EXhv1gIYpcuUaMzvYsftR4MXm9ingATN7t5ndAhwEfrLX4xs2MzsE/C3wUXff6OiKxPy3EZW5PwscNLNbzOxdbP7YfGrEYxqFU8DR5vZRYLv/wwtCkGfyl/GYmb0X+B2btzX+DIC7nzOzJ4GfsVnGeNDd3xzdMIfmy8C7gdNmBvAjd/9MFOZvZh8D/gGYBBbM7Dl3vy8Kcwdw9zfM7LPA94CrgMfd/dyIhzVUZvYtIAtcZ2YX2fy/9seAJ83s08B/A/ePboTDp9saiIgELHLlGhGRKFHIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhKw/weTxCkIdcRX9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}