fastai/tutorials/meanshift.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering with pytorch"
   ]
  },
  {
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math, numpy as np, matplotlib.pyplot as plt, operator, torch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centroids = np.random.uniform(-35, 35, (n_clusters, 2))\n",
    "slices = [np.random.multivariate_normal(centroids[i], np.diag([5., 5.]), n_samples)\n",
    "           for i in range(n_clusters)]\n",
    "data = np.concatenate(slices).astype(np.float32)"
   ]
  },
  {
   "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": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_data(centroids, data, n_samples):\n",
    "    colour = plt.cm.rainbow(np.linspace(0,1,len(centroids)))\n",
    "    for i, centroid in enumerate(centroids):\n",
    "        samples = data[i*n_samples:(i+1)*n_samples]\n",
    "        plt.scatter(samples[:,0], samples[:,1], c=colour[i], 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": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXt8lNWd/98nmckNCJcJlySAQSGgoogERFfF1IgKLsR2\nWy/drqIrtYo3tlittSTbunax6/ZCbReraFurdreV8hO8YYNaERFKUdCAUS6SBEgCIZDb3M7vjzPP\nzDOTmVzIbZJ8369XXs/Mcz2J+Dnf53O+53uU1hpBEASh/5PQ2w0QBEEQegYRfEEQhAGCCL4gCMIA\nQQRfEARhgCCCLwiCMEAQwRcEQRggiOALgiAMEETwBUEQBggi+IIgCAMER283wE5GRobOycnp7WYI\ngiD0KbZt21attR7Z1nlxJfg5OTls3bq1t5shCILQp1BK7W/PeWLpCIIgDBBE8AVBEAYIIviCIAgD\nBBF8QRCEAYIIviAIwgBBBF8QBGGAIIIvCIIwQBDBbwd12suf3Iep097eboogCMIpI4LfDjZ4anjG\nXcEGT01vN0UQBOGUiauZtvFKgdMVtj1V6rSXDZ4aCpwu0pX86QVB6FlEddpBunLw5aTRnbpHnfby\nk6b9bPXVAXT6foIgCB2lSywdpdTTSqkjSqmdtn0jlFJvKKU+DWyHd8Wz4pW2fP4Nnhq2+urIS0zv\n9JuCIAjCqdBVHv4zwFUR+x4A3tRaTwLeDHzvt7Tl8xc4XdyclMW9Kaf1mJ0jg82CINjpEsHXWr8N\nHI3YvRB4NvD5WaCwK54Vr1iCHit6t2yhtsS+K0VaBpsFQbDTnVk6o7XWlYHPh4B+bVrbff5IwY4m\n4rGEvStFuq1OSBCEgUWPeAtaa62U0tGOKaUWA4sBxo8f3xPN6TKiZd1Ygg2hgdnIfa0N4HZVRhB0\nzWCzIAj9h+4U/MNKqUytdaVSKhM4Eu0krfUqYBVAXl5e1E4hXokm7gVOF03aR5P2U6e9pCtHCxFv\nbQBXRFoQhO6iOwV/LXAT8KPA9s/d+KxeIVo0nq4cpKhEnnFXkKISguLdpP287D7CNUmjwq6TfHxB\nEHqKLlEbpdTzwGVAhlLqILAcI/R/UErdCuwHvtYVz4onYkXj0SL6FzyHgsdvTM4KWjt/ch8W4RcE\noUfoEpXRWt8Q49DlXXH/eCaajx/ZEVzgGMprnmoqtRtQwf3RLCFBEITuQsLKTtIe0X7fe5xK7Was\nSmaOc3iwk7jAMRQIvQlI6QVBELoTUZVOYrdvyv1NPNVczjXODF72VHNrcjbZCSkUOF3s9J1kq6+O\nRxo/Z2biUF7yHgleZ4m8vfOw7xfxFwShKxAl6SR2+8ZKtfzUV89xfNAM96acxgZPDbcmZ1Pe2MRB\n3cxQf30wPz5S5K1tZ+weeVMQhAga66D0TZhyOaSm93Zreg1Rgw5gF1KghajempwNzQQj/GucGTzQ\nsIeDupkm7WOEchofX4cydC5wDOUD73E+8BznAsfQsPRO+zZaG2KJuYwNCEIEpW/C5sDE/+nX9m5b\nehER/A5gF1Igqv3y/dQzADjfMZR/b/yMg7qZTJXEJ756dvnryVRJ7NL1/LH5EMkqgY98J9nlrwfg\nqeby4PUWJyIEvj1i3pWTtwShXzDl8vDtAEUEPwqRUXSsQVbrcywRvjU5G0+TH6/W7PCfZKxKZkrC\nICp9R9nja2CXrg+em6mSuS5pdDBN07qn5f1b926PmMvkLUGIIDV9QEf2FiL4UYgU8GiCbtXAOaG9\nNGk/1zpG0aR9lPubeMtzDNBckzSK6Y50nnFXMFYlc1A349JOAHIdgziHwTSjQWuSVSKbvMd5yXOE\nKr+bg/4mrnWMYm6Si6newRKtC4LQaUTwoxAZRUeLqiMj8LzEdLZ66yjzNwYj8hSVGLzmAsdQ3vce\nD27t4wBN+HjBc4hpCYMBeNt7lBP4caoEshNS+HJSSovnQnjnI4O0giC0hahDFCItkWgWiV3Ip3oH\nh20netIAzQWOoWzw1HB+HWRnhIS7wGm8+PLqKt5I93B2wiCud45mjnMETzWXs9VXx1iVbAaBI2it\n8wHpBIQBgmTdnBKyiHkHsJc0tjoBE4GPZkhAVIcoBzcmZ3Jjchbve4/zYNH3OX/aeazcuYk67Q1W\nyly56z0eyruc3f/xq+CgbXZCCvemnMbNSVn8KC2X7ISU4DPL/U38yX0YoEVd/WhlkKUWvtAvaayD\n7S/BznUm66b0zfZf01jX/e2Lc/pN6NcTEW20MseRk6aatA9TPkGz9T+eYM+jqwBYNrcQXl9DVu4Z\nbCz9iK3zb+dE5RFOBI43L/9+cMDW/jYRa/DWjlWR0/77S6aO0C+x0itnXAezb2pf1o2kZAbpN4Lf\nkdzzU+0cohVFi5w01aT9vOA5xO7/+FVQ7AGaKqv4/pVf4Se/eoIdt9/BicrDwWN7Hl3FywmD8Tzw\nzy3ab9337MRB4CaYKdTW7y+ZOkKfpC2rxp5eGe14tOslJTNIv7F0OrK606nYHfbUzA2eGuq0N+yZ\nlsBekzSSK+qcVD6zJuz6QgrRFU3ctODL+CuaKIxY8fHt1b8n42ijEXZC9hEYAd/lq2err473vcc7\n/fsLQtxiWTU710U/bqVXRop9Yx1seR42/FdLqyfWNQOQfhPhdySiPRW7I9JaadJ+UlRC1Fz9I8NS\nyFv3P2yffzu1lYcppJB7uIeFLKSYYpaznBxyAFjDGlIyRzL75f+hekQqP27axxiVzOmJabzkOUKT\n9gPQrH1c7xxDgdPVrgqdQudoqIbtq2H6IkjL6O3WDCB0xBZCUXvOLNi3JXp0X/ombHvRfB4/Q6L5\nGPQbwe8Iljh2pB69WcnKT7P2MdE5BtDBt4SHUk/nfe/xYIeww3+CwZNO42uvPMtvrv4GGys3spCF\n5JDDalYDsI99bGQjqZmjmL3uf7hi8rmU6gYO6mYOaw9erclLTKdZ+3jJWwXAzUlZpCsHf3IfltIJ\n3cz21bDhfvP5H5ad2j2k0zgFzpkPzhQj7ttfMsK9cx1sfREO/h0O7gBPkznHEnWrM/A0meGzqfMl\nmo/BgBR8i474/mYlqwRe8Bzi5qQsLnAM5RVPNQd1M081l3NvymmA8dhPNnop1Q1MmjiJy35WxKtf\nvZNiioNiD1BMMbXUMutnP2VKbi5XJmeQ5D5Gs8/PCJzkOgbxkucInoQhXOsYRbJSwei+Sfu53jla\n7JtuZPqi8O2p0BWdxoDDsl+2v2SsGU8THCo1x1wTYOx54G0KDcJC6POsWMtyCBb9xsM/FSzf+wLH\n0GC6ZST2VMwLHEPJS0wPTp46ji+YL39Ce9npOwlAskoE4K09O9l4dxHDGMZylofddznLGcYwPr77\nEUr37OFF92HW+6qpwsOQRAdznS7GqmR2+E+QrBJICdzTWj0rRSUGraRYbRdOnbQMI9L2yLyhGt59\nzGzbw/RFULCic53GgCEydXLK5SYLRwHlO4xNM/3LpjOYeClkT4PG4yayn3Gd+fzB822nXkZL0RxA\naZsDOsK3rJ3WLJLIgmlbfXVM9KQCKjhZ6n3vcbZ7T7DDfwKa4faUsSz78C+8NP9faaqs4ioKySGH\nfewL8/Av4zLWVK5h+/zbuf/NV2gcPwiUqcHzvvc4B3UzeYnpWPYRtJ4pJPZO1xHNjomM2NuybKxO\nQ2gHkamTVqTfWAeOlHDfft8W0wmU7wBHMlSVwYFtoXvNtEX6kVk70VI0B1Da5oAV/Drt5WX3EUAx\nxzkciD6IG22A10q9vDkpi7c8R3nBc5h5iS6ciencmpxNyrF6Ns6/jaZK472vwWTsbGQjtdRyH/cZ\nsQ/sr608zNcuv5qL3nuemSPH8ZbnKHOcI8Kea5VpiByclXz77iGaHRNp8/S2ZVPTqHmx1MN1U5y4\nUlXbF/TWMxvr4KN14G02An1OhMfeWGesm7zrQr78sXLY9DTMiLIU9pTLTURfs9fc88A2SB8DdYda\nDvb+5aehzmD6tdFTNAdQ2uaAFXxjjZi0xxSVEDM6jhTYLyeNptzfxCe+kxz3h2yUZJXI1IRkhigH\n6S4Xi2+7jeLi4uDxNaxhSOYofvWr/+W737qLNRXhaZuZNy/E5XKxw2/eFFJUYgthjzZ3QLJzug57\nxB7Nw7cidsvambyg5Tk9yYulHn642Q3AHdOT4veZ9gwaC2tgdt8WI/bbXjQWjtUR/HWVGaCtrYC6\nSiPw1XshY4KxdlKHmuOjp5jr7Bk8Fh+tM2KfPS203/7mYA0Kt7eSZj8o5zBgBb/A6eK438vn/oaY\nk5li8b73ODv8J9nhP8n1zjFc7xzNR96T7PLW06R93JicRVFREc3az4/+/QcADM8aw3kv/5KySWcy\n7eVfsm3+7dQFJl9d+OCdZHz3Vk7gY1rCEM5MTAuL2K1yDLFm2gpdQ2TEHitqt87btxGufbb3MnCu\nm+IM2/bUM1M8dVznfR0aC9onfFMuN6JuRfgKY6FU7DSCPHaa8eHtmTmuCUbQk9NgzNnw+Xtw4rCx\ncVKHRp+ANTxCtK0XkMwp0dM4O2rj9APrZ8AKfrpyMDTBwQ7vyUD1Ske7Z99aKZqmBPJINnhqbLXt\nVcAuqoL7/5lcz2Eqn/kzf37zdb6Y4OJlTxWDJ53Glet/zYb5tzH15n9i2IO3ADAtYQjLUnNaPH+D\npyZYkdPK1LHsqGuSRkpxtA7Qmu/e3sycyQtg2yooWw/vr4T8oo49p6twpapglB1ptbT1vTPPvMX5\nDmz+DThU+4QvNT08g6axzlgvTXWQnmmEfex5JkK3BPXMK2D/B1D1Wei69EyYcEFI5Kdc3nrEPfFS\nOFJmtpGcio3TD6yfAa0Up7qGbHqgQJr9PlYNHasDeMFzCIDJ372dCd+8nqrMUaQAJ/GTguLiyedQ\nt+n3JLmGMZgE5jhGBK+1l1COrItj5eG3x44SWtKa724fZG2ohi0rjS5dsCRctHevhWNl5nM06Wyo\nhpduMh1CtOd0tDOwzj/tRs3aI9FFO9Jqaet7p+iM8Fm2iAJ2vWL2RU6UsoS8thySh0DzCfMWcPFi\n0ylY93njv0zE72kKdSh222XfFvMGkTU1PPo/VWumHyyiMqAF3+5/d2bw03QAWcEZsGcnDmJawmCy\nE1JIRpE8xsyQtZYrPKibyUhwMnPkWHb4T3ISPzt8JxjiSeQFz+Fgxk+T9nNjcmbUgVqrg5HB2o7R\n3ih++2p4KzAEkzTInG/39931RrNmLYl+bdl6mDgv+nM6OthrnV/v9PC75OiiHWnvtLXtFJ0RPnvx\nsxnXhWwe676W2OfMClk+42fAl+4JTcBqPA5HPoXKXeY6FeX+nibz3T4QbAm9NWYAfV7AO8qAFnw7\npzr4aS9z8LL7CC94DjMtYTA7/CeZ7khvUfnSSrW8JmkU1wB/bD7MJt8xDupmPvKd5HrnaJq136R4\nhqUchLf1xuSsU/xNBzbtTZWcvgg89ea/gCX2dt8+0sZpbcD33cfCo/mOTuqyzjvtq05yjkQXbcve\nqWnUPLHdzXVTnGGdgt3+6VUivXdrgpXCpF96m4yoV+yEi24x0fmUy6lp8OCy/newiX2Nx4nLbtlY\n97fuYx8ItjqDvA5U2uxnDOiJV12BZQW97D7CJ74GAE5PTON652iatD9YA/9P7sNc4BjKzUlZXJc0\nmp807eeE9jI3yUVCIETZ5a8nRSXylWQzEGyNBwAywaoDdHSCVCycg+CcG4yQT15gIvay9eZ7JFaH\nEHkscv+pePtWJzU224h2ax68Zd28WOoJ7is75ucb6xopO+Zv3wO7k8hCZplnwbBsE31vfhY8zeb7\ngW3w6dsAFP3wUc4991z2pEw2Yu0z/w/sOXySc7//MkUPLA1l3YC5/9T5oewda1KVNZkrmqffHvrB\nBC2J8DuJvSzyDv8J8hLTmet0BVeuSlGmT7Xq7NybchqPNe5jh/8EjU0+jmsvldpNpkpmjmNY0KdP\nUYk8464IevSxxhhkZauWdDY/3u7B79sY8uKvfTYk1vZzrQ7BUw8NVfB/N8DeDebaSx82HYWVwtnd\nufvRrJviTc28ecAHNPPb+ald/9DOsO0PxqtPGWasGzDfh2WDt5mi4mKK1+8GIL+ggJJ/v5FcfZA9\nh0+S//P3qTjeRPGqP4C7nqJZjpCfH1mioWKnsYXs+yCUotkeT1+ydAR7ITareqY9q+YCx1De8hxl\nWsIQtvrqeKBhD1NUmrlYw0HdzFiVzD0p49nlqw/eN9p6uk3aT5P2BVfcAplpG43O1sGxe/BXPg7Z\nM6G+ygzizooYwLULuHNQyPcfMTHUUZSth5zLIGNZ19ToaY1o1s09M5LYX9fEPTOSumeyVmfy0y+6\nxUTsfq+J6kdONMJ/YBtFG8qDYg9QUXWM/AeeYtXdhSz+1dtUHD0ZPFb8zDo4MpmivIj7T7k8NBZQ\n+mb0yVftFXLJ0hEsIgeArbRNaybutY6R7PObapgulRTw6jXnqCHB7ByzYlao7HLkuropKiEQ9Sd2\nyWBzf6WzJQ3sopyWES7k5R+E595PXwS1+02a5tU/hwkFJrofnA1Tv24soZzLQvfsbNvaK9j287ZU\n+iir1Wyp9LGl0tf1k7U6E/kOzzZ58ltfNJk458w37d+8hie/f0/YqYUUsvH4Rq75wQsMYxiFFAZn\nqwM8+cER7sq8EBeEZupedIuJ7K2BYCvP3z447G0yA8htCblk6fRfOmOV2MX5eucYbk7Kokn7g8XW\nbk8Zy/ve47zgqQiWPA5ZQ76YEXs0cZeZtl1PpChbWTlfvBvy8O3Hd70AjTXwyl3wr++F7KDc+ZAx\n2UT2XUV70yut896r8LH8IpMFY7d5unSyVnsi39beAjzNZuuaEDzmml1IyQ+2kP/AU1Qcb2pzTYms\noSmU/Ox+XGPNPjY9baL6ukNw1XfNvk/fDp/xa5/8ZR/c7QczamMhgh+DzlolFziGstN3kjnO4cHF\nyO0LphQ4zZ/eEu9waygxasQu4t51xBo8tXvyu9eGjucXhV9jfXbXG7EHmHilOTea199VWEI9N8cR\nzMaJFulfN8XJexU+3jzg48Isb1jn0OXZOq1FvpZ4Wlkz0PJcZ3L4FqD0TXL1QUpWLCa/6MVW15TI\nGp5Gyd0XkpvWGLp+xtdMNk9teUj8o2Xn5MwKZgLZn93XvfpYiODHoLNWyfve42z11THVO5gvJ6W0\nEOtYi7CIqPcMlvfurg/l2adlhKdfRk6cstfSsaL40wtg2i2wfyO4G2FjEUy1TSrt6hm3lkf/xHZ3\nq5G+K1UFIvtm5ua0/r95txZhi7boeGQEPXW+ScmMtFw8TeQqWPXfM7jm+ptirinxu29/k9zTCGX7\npKZD5cdmJu/4GaECbBMvNRaSRayZuv3Aq4+FCH4MOiu87e0wZNC1d7BsmoPvwueBjJorHzeZNpcu\nh4lXmfPGXxKeR28Xe9cUc61rCtR+bn7A+PxWhk/2zJD/b5+81VnxjzWRyhLvuTmOYHbOhVlehqc4\ng/vXfOoBBYummhTPrpyF26LziFbzxp4lYxfdKJH1nnVPsfjn78dcU+I+7mPxT/5AyR3nk1t3CIaM\nNPfyNsG0QvPW8MXfWs64tVfStDJ4rPb1A68+FiL43US0DiPauIAMuvYs9og7aZAR7OGBjBqf2wy4\nFqyAA++Yfcc+h5rSgHVDDam4gmK/cLU5b8i5Nez4sYsRkyEpzUzW8rnN9VkzQ4ugdEVKpl1Qo4mz\n3bt/84CPy8cnct0UZ4v9AGkO87bQlbNwW3Qe0cTTEmRPU2j2rDXRytrfWMeeA5Xk/+TdoIcfc02J\nqjXk/2IbJT//DrmNx+HV/4RDu0zRtUO7jPDbc/KtzuXAtlDOv5XB088Rwe8E7R3Ytc6z6uhDKJoX\nC6dnsYuu5bHXV8F7j8GY6XDG3IBHXwPbnzZi75oCT79RxJ83PckvlpQwcV4uZevh7R/AxLv2cPnX\n8ln0jdtYeFoRnnoT0V8UEHRvQ6gWT3tSMtsqemYXVEvI7VbM3BwH71X4uGdGEhdm+ZiVmcgdbzRy\nxrAEluY5KZzoZNpIE+FbAt+eWbg1jZrVH7nD3gyi0Wo1TbuVozEDqNMKg2mYZE0N7q85foL8u39M\nxXFTIqGtNSUqqmvJv/sxPlx2Aa7Bgd+lPmLmnX3Q1nrziFZWuR8jM207gWXHbPDUBPdFmxEbsm00\nNydlSTTfi9iXHbQ8+YvvN/vOv9Wc01AD65cYsU91wf+VFvG7TcWcoIL7ns1n7E17cE2Bzev3cHVh\nPtV1FTz2i2Ieur+IvSVw4TI4uMW8LWx6LDTLNnLZRGtGcPXu0MzgyJmykd+vm+Lke7OTwqJ2+6za\n1/d5efOASb+8Y3oSP93m5p1yP8/s8pLmUEwcnsC3ZyXz7ZnJwQ7lie1uahqjl/GweLHUw+PbPDy+\n1RP2vEisapqDtv7GiLsdy7L5aB0cDqxT60w2dorl7wf6EdewIdy2eHHY5WtYQ1pmMi+//DJpWWlh\nKZkAt922GNcZ55gvg0eacsrDAp795mfN2rjTCs1bRNMJsz9lSPjM336ORPidIJodE82Tj6x2GQ2Z\nMds1dHTZQfv5VvT/2etGrAHW1xTxFqGFbGpOVLDwxnzm+1axPnkxtc2h5S/fohjehluTijjwltk3\noaDlzNx3V0DldvNG8d5joQHifRthwSonzI5d9MwejUceq2nUNHg1S2c4g/uWX5SM29fE1Axj7ZQd\n81O8qZl7ZiSxpdJHTaOfX+7w0uDVfHtmcvA+kW8O101x0uDRYW8GMYk16Gmvc3Nwh8m7t/oZy06x\nlTQuKroB9rxF8fMbAUzq5cP/RO78+ZSUlJB/8Wwqqo6Z33Px1yh65FFYG/D5h4yCEePNm8MZF4fe\nIhIdZnukLHwlrAFCtyuLUuoq4KdAIvBrrfWPuvuZPUU0O6Y9ufLRxF0Gb7uGLSuNpeKph8uK2j4/\nlsWzdwOMuKiGHduehObQ+YUUstG3kd9zDcOaW07+2Tn4SWrr7wJcTCiAf3q+5czcTY+FvhesCJVd\nKFsPLFbc/GwSjYQXQbMXRbNEONKKebHUROCXj08M7ps4PIEXF6QZS2anm5f2eNhbB38/0khNE1yS\nHXjJtwX41ptDg1eT5lDBZ357VihtstXMnli16iPXqbWqVh4qNStZOZJDSxbWVsC871H0+K/gxJd5\n8p3PKbn3H8idZqbS5ubmUvLE98j/14e47epZFP3ocePPT1torJzZ/wJDM0NtgNBEq1ETTe7/yIkD\nxsqx6FbBV0olAr8ArgAOAh8opdZqrT/uzuf2JtYkqtai9VhvAU3aFyy4JlH+qaEjtm1hifzkBaH8\n+789BePnQPYsF68+XMI1X8un+kRFm5N/MtKz+K9/KWHvShdDJ5hJV++vNLNtrZx+a2buZ69B/g9g\n3Gzz/GufDWX/bF8NOwrarmdvOoFmdlb7eeSSlLDc+2vXNPD0ValMHG4E3eoMAFwpUNMEaQ648Uwn\nM8f4afTCj7c0s+ic0CBug0e3GC+Ym+Pg9X3mjcC6X1T/v7VcdrvwWwuQl+8wxzLPNtu6Slj/Q0gf\nRdG8Sdx1dR6urLFmYRSAxjpyx2fy4fOP4XK5oOxtM/g7fobJva/8GMZMjl4rx5HSspLmAKG7VWUW\nUKa1/hxAKfUCsBDot4IPbUfrsd4CIgumCR3ngiWhvPq2iGbnfP66ydwBOPAWFIzM5d2tJeTn57Ox\nIvbkn4z0LDa8XsKxV3LZCxzfCx/8wtxn14tmPACMnTTsNLOAyoF3QoIfOWFrdKIZfLVy6O32jWXL\njBuieGaXGSt66J0m5oxzsPyiZPbXNVJWq7nl1UaeuqSRES4XZcd8pCfBhKHwg4tT+fr/HuQELh59\nv5nThibwzkFTSTPNad4arLeKNKdibo6De//SxJsHfMEOZemM0FhCVGLZOpb4WoOlF90Cw8eZSVKj\nJpkIv3KX+Y9YV2l+sqfh8nvN/rK3YeYNwXVyXcOyYV95KM8/2kQqKxPI22Su7cd59m3R3YKfDXxh\n+34QuMB+glJqMbAYYPz48d3cnJ6hrVTLWJk5kqLZeVrz6CHc37fbOZMXGA89eZj5npkHOfnGGho7\nIpdVq1ZxzTXXxJz888df/o6Tb+dyxlXw0XNwtMy8JSQ6jT1kXwzFepZl5US2vaZRU/yXUA79xOFJ\nYfXub3nVCPqEdGOlDE+GM4Yl8MPNbl4o9VB0UTIP/7WJ7b97hLx7nuXeVa/xQmUOADuq4H//+gkH\niubimH0zo259KCj24wfDrMxEfrylOZiNY03yslI8l1+UzIVZ3rYnacXKZbcif6ukAZg1ait3wWl5\nRpTBRPLW+rcQmqXbWAcfPG+26Zkmmh8/w9TgibW2beRrXz/Os2+LXvcNtNargFUAeXl57X0Tjwti\nDbSeSqqlDNp2D3ZRh/BFTOxpklaFzGETzL6kNPPzVjF8fnAP3/nfxa1O/vnmtxZzfV0Jk6fkcrTM\npHJe/AC8918ma+fi+0Ne/u614SWXIweXXyz1hOXQW9Q0au79SxNltZqJwxSP56fw2JZm3in3U3rU\nz4R0KKvVFG1qJukvP+Lk+kcBWHXHlVz16Cu87zudCxI/55ffupqmo5Ww/lFOjElg9j9+l82Vfg6c\nhKUl5v4QPU/flaqYOLwTk7Ps6ZBWJN50wnQAObPMsSNlcPqFZsKUBiZdanz+gzvgeEVoaUQwWTgX\n3dLSmmmsM9lAKnC9M2VARvSRdLeylAPjbN/HBvb1C7pyoFUGbbuHyNx3KyPGKoBmvQ1Yx2v3w9Zf\nQNYFRmuq2cNPfptPrbui9ck/dWt4Ljmfr5eWkDkol5pSeO0+Y+UkJoWL+vRFLdthZ8EoJ/ua4e7z\nwqNoe0fwky+l4EpVzMxM5J1yP5sr/dx8toPEch/bf/dIUOwBDlVW8MdvX43jqz/lT/93DyeqQ5lF\nG59+hCsSYOnih/mg0sc75abjmDPOQYNHU9OoW10s/ZSx0iEhNAkqa6r5bkX+1taZYtaz3fS0KZMw\nZooZdK3ea7z/fVtarllrzaIF04HYZ9IOYLpb8D8AJimlJmCE/nrgxm5+Zo/RlRaM2Dkdpz11aiIt\nnisfN9spT/IaAAAgAElEQVRYdkpDtfHYpy+CmqM1vPDf+dTWGYFsa/LPseYKfuvI55v1HzJuiouF\nq83kLOuZ9me1VmBt/+8Vg+5PYr8HxtraHhlpW7nzs8coNh/SjEhRPHVJIzPveTbsfoUUsrF6I8d+\n+dWoZYXf//NqHnngHsAI4jvlfk4/oU3uvrNlJlCnyjDYB3PtJYqt7xZWNo0O7P9onRFwnxeu+Ddz\nzkfrTGll+3V2sc+eZiL8ATSTti26VfC11l6l1BLgNUxa5tNa613d+cyepCtnycqM245zKqUKLDvF\nWpDETkO1Seu0fMVxuS7uvO82iotDefhrWMPw5Cx+MP93rPjT4haTfxZ/8zbO3eviysdNls7X17V8\nxvbVcNqNmh0FHqYMcpJGeKQca0ZuZBrm6o/cPL7Nw7emOSjISQh0BBls+WsJMy/Op66q7cyipOGZ\nJN/xMj8rHcybBzwszXMyZ5yDWZmJgLtF4bVOl2GwD5jGyuSxPs+0VaGz/kTlO4zQW9k9kZk21tuC\ntfC5tU/sHKAHPHyt9XpgfXc/Rxh4nMrqUa1ds311qNCZlelzxaAiGu6Ex35hDgxLzuLrzSXMOiOX\nXywp4b5n86k5Yd4Alnx9OQtPK2LyXfDR80aj7CtkNVTDH28wWUB12sPzQ0KRcphVkqFadGBRrZTA\nJtWpwsos5ObmcvPPXmXl7Vey8XjszCLHsEyG3r2Oi8+bTO5wxbRRThZNNR3KHW808k65n2kjPWH5\n99HKMHTI5rEPmLaWLWP34KfONz869DtzYJuZuBWou0Nqutk2Hjf7L7rFnNdP69qfKjI6KPRZOrt6\nlJ2GalMg7cJl4EwzYv/XFWYmbPaEIuYAO5KfZGluCf6PcvE2wDdW5nLuP5cw5+J8zvPexuRtRWx4\nLnymrnNQqI3vrwylfM5pdHLGFaFIOVaevSWk0WroFE50kuYwaZOWQDd4NN+elcy4M3IZeuPPOPbL\nr8bMLMr6xpPo0ZNwJsAvd3j53myTDfTjwEAwQGOoQkgYsdrWIZuntWyZQNolYDz4i24xXn7OLFMT\nJ++6UD0eZ4q5T+mbsCPwxrVvi9n207r2p4oIvjCgiLSB7AuZvF1sZr5aAn1ou9nW7oUvTyzirqvu\nYudKM8ZSHVhq9eTbuXzT+yHpaS7GnA9TrwN3gxH80wvMWMHr98Ph7TBisrlmQgEU3KlIyzDi2FAN\nEzY4+c5F4VZJpMhb20iBfWK72ybQZkbuBPfnNL14d6uZRYd+czdnfGcdy649l6RtbmZlJvLEdjdH\nm0LJcqkxFCJW27qMQD18DpeGD+LaV6iacnl49o11jSL8raEtO6cfr3AViQi+MKCItHSsDmDO8lBR\nNYt5K+H315hJUsfKQL3qYtotcOIAXL3SCHV9FbhcLhprYNfvQwufW2WSP3revCWA+W4v3GZhWUn+\nHwLnh/ZHDtJa1k+DR7M0z8ncHAc/3tJMoxe+Nc1BqkPR6IWiP+/i+M/m4z5WyZWtZRYdX8PRn17D\nD1yv8IHvdPbXmZRMq9zCJdkJLDqnZcRur9kzN8fRPYunpKbDrBvCJ2plTQ3fWvvtYj3L5vs31rXv\nWf14hatIlNbxk/qel5ent27d2tvNEAYQbWX6NFQbK+bD35hIf/hEs25tWoapcGm9LQyfCOlZsP9t\nI+oQ6kjqq+Hz16Dwt+Ca2PJ5DdVQ9Ds3v0t2873ZSa3aItZKV9+bbc6xomwrVXPlXyv596/m4T9e\nGbymkMJgZtEwhoVlFgEkDM1kxo/e46B2ccnYBKa6Emj0wme1pmSDVZ6htTa01e4ux1pExSqKNvum\nlmJtnRPtmJ1+EOErpbZprfPaOk8ifGFA09rM3LSM0Hq2ngYTqR8rM8sYHtsLlz5sBF1jSjpAy1RL\na1LX0TL47FV4+4PoSyd+55tOckrbtkUi7ZMGj+aDw6bcwYulHpZcnMm7X72FN379SPCaNawhZUQm\nv3/yt9x++zdZUxWeWZT2Dzfzj+eOwZWqaPBoHt/mYXgyHGs2JRus4mtWJB/NwulSO6c9RJvAFe0c\n24IqMcV8AM28lXr4gmDDsni2rw7ff/H9xnsH+PQVI9p/+rr5bi1wElnvvr7KZOWMv8RE/Rpznb3M\nAphOpnSl4uZxsRcWsbCsHVeqClawfOTiFC4fn8jcHAeuVMXrT/6Qy255KHhNwtBMznxwHTPz53Hz\nz14lYWhm8Niy736f//xhETecaQS7cJKTy8cncixQIfSMYQk8sd3N6p3uYO39yDZY4wpt1dTvNI11\nJmq3xHv6tWaN2lj17K1CadtebFmbf4AiEb4g2IiVtpmWYUodb18Nx/ebwmi1nxvv3Z6JYx8Etrz7\nhCSTj99QHb5gukVnlz5c86mZgZs73E2qQ4GCIfMfZPAhP97Nz3D2Q+upSJ1I8aZmpp2Ri+vedZz8\n+XyW3H4bKx4pDpZssJY+vGdGEp/XNjJnnIMRKYofbna3WiytK9fEbZVIr709VswALpQWDRF8QbDR\nWqqndax6tylvfLTMZOLYOwf7IPCFy0x2jjXTNta9T2U+ARCscb+pwgj1zmp/MFtn9hjFkGu+y9du\nv5O0oSN46ws/98xIYsLQBHYcmcIbD7xH5VmjgksXvnnAxyXZZvLWvX9pYm8d8IWPqyYksnSGk0Xn\ntL6soX1rta3LB3Mjxbs9g60DyK5pDyL4gtAK0QZ1d681Ym9l5NiP28U7VrmHSNozn8ASUKsevWWj\nWDXprUqWz3/iZmeNn2Uzk3nlcy9ry4ZSftB0CFsqfUwYmsC0kQm4p4zkzQM+Vu9088Ehc3x/nZ8n\ntru5Z0ZSMGPHnp8f2RZLzKNNxuqWqD9SvCV67zAi+ILQCtHsFnuZ49eWhg/CduVkMDuWgFr16CGw\n7KBXgyYYgac6FO8c9DNztJc9x/yU15vrJ6SHvPbHt3lYOsPJnLGmQNo75X5cKXDgBPxyh/HoXypM\n44ntbnZW+4L5+dEWUo8l5l2amx/LupHovcOI4AtCK0SzWyxRf/ex6IOwp0Jb6aGWcM7NcYTVo7fW\nobWwZsY2es16tvWeJjw+P7MyEznWFMrhL5zo5PV9XgonOUlzKmZlJvLYlubg2rdW5P5OuZ+kbe5g\nJxNZLjmWdRMt6o9Ke3z4AZQn392I4Av9mvZU1GyNWBG7VYphznJTLwdMB3Cqz4k1cGsXVEtAo9Wj\nt86zSr+lOsx6tn8KROo/3Ozmjf0+ymrNhKniTc1hbwovlnp44orUMNG2Fi5v9MK0kQnMzXGErbML\noZx8gDumNEHpmxzN+RIv7Ettn3/fHjHvqHXT0bz6fpCH314kLVPo18RKs4xFQ7UR7obqtu/7diBD\nx756VnufE8n0RS1n+kLIyjFiHhsrbRIU35udFDZD9ropJtXSWjil0UvYAivWtat3usPu6UpVpDkV\nv9zhAQXFm5pbtOW6KbbsnYB4l256o11tBozIWmUSYmFZN+0VY6sTaW8qZkfP78NIhC/0azqaAfP+\nSiPk7noz4aq99z3VTBuLWG8S7fXCGz2hdfzumJ7c4vi0kQnUezSbK/1ccZoOirQrVbW68rt9gle0\nVbjCrJuAaE/J+RLfy2plvVs73eHDd/SNYAAN/orgC/2a9gyi2m0fy4BozYiIZhN112Bte73wVIcK\n29qxBmovGWte6FMjFjVZdE5ScLFy++Cs/fnWguat2jQB8R4B3DG8g79oV9LRTmQADf6K4AsDHrt/\nPmuJsWkiI3W7yMfy2zs7XtAZLNGOFlXbB3xf3+dtIeyWqIf58RGdTLsHYYW4RgRfGPBE5s5Hi9Tt\nIh/LvunsjNnOYC9xEC1jxn4sVlpla/ZRt0ykEnocEXxhwNMeO2b6IvDUG28fWp7fUG2OX7q88yma\np4q15KG1CIqd9tSvby2K77HyCUK3IoIvCK1gt2mcg0wEnzSopeBbNe0nzuuddgKtDkBEq63fEbpl\nkROhxxHBF4RWaI+VY+3bt9FMxNq+uuctHYBFU5NIcxj7ZushH0tLmng8P4W8MYmd9uDFw+8fiOAL\nQiu0x98Hc8wqkjZ5Qc+0LRK7KF+7xtTCWVrSxNs3DOqdBglxh0y8EoRWiKxx3xq715oIf/fa7m9X\na9Q0ai7OTmRCuuLx/JTebYwQV0iELwhdRGcnX3UGexbNi6UentllqlzmjUns+cYIcYsIviCcAj05\n+ao9RGbhNHjMQuM1jVrSKIUgYukIwinQ2do5p0pNo+aJ7e4Wywnaa9pYNXAe3+ppXz0bYcAgEb4g\nnAK9Zd/EyoePzKKRNEohGiL4gnAK9JZ9014hlzRKIRoi+ILQhxAhFzqDePiCIAgDBBF8QRCEAYII\nviAIwgBBBF8QBGGAIIIvCILQA7h1I5+5t+PWjb3WBhF8QRCEHuALTymfuDfzhae019ogaZmCIAg9\nwDjnlLBtbyCCLwiC0AMkqVTOSJreq23olKWjlPqqUmqXUsqvlMqLOPagUqpMKbVbKXVl55opCIIg\ndJbOevg7gS8Db9t3KqXOAq4HzgauAp5QSkmdVkEQ+gXxMAB7KnRK8LXWn2itd0c5tBB4QWvdrLXe\nC5QBszrzLEEQhHghHgZgT4Xu8vCzgc227wcD+1qglFoMLAYYP358NzVHEASh64iHAdhToc0IXym1\nQSm1M8rPwq5ogNZ6ldY6T2udN3LkyK64pSAIQrdiDcAmqdR2XxMPNlCbEb7WuuAU7lsOjLN9HxvY\nJwiCEHe4dSNfeEoZ55zSIRHvCJYNBPRatk53TbxaC1yvlEpWSk0AJgFbuulZgiAIHSIy2u4JT36c\ncwpnJs3uu3n4SqlrgZ8DI4F1Sqm/a62v1FrvUkr9AfgY8AJ3aq19nW+uIAhC54mMtnvCk4+HPHyl\ntW77rB4iLy9Pb926tbebIQhCP6c7LJyesIVioZTaprXOa+s8qaUjCEK/ob0Do9EGXU/6j7GlcR0n\n/cdO6b59IVVTSisIgtBv6OjAqD0q/7h5E0d8B6AZZqXODztm3dervTiUI2j92CP6vpCqKYIvCEK/\noaOia+8gzkq+CJoD24hj1v282sMn7s3U+CoYmjCSTz3bAIJvC73t0beFWDqCIPQb7FZNWzaMWzfS\nrBvJSMhmtCOHwQnDmZU6n8EJw3HrRrzaw+mOaXi1FzCiPtIxHicpHPEdwIeXUYnjGe3IaXf7ejsX\nXyJ8QRD6JW3ZO194SvncswOAcs+nOJQzaM/sdX/Ep55tZCRkU+0vx6EcnJE0nTL3Njw0MUgNIxEH\nR3wHcHmzGJw0vEva1N2I4AuC0OeJliETy96xzh2emElGwljSE1z48PKpexte7WVC0lSO+Q4DUO+v\n43THtOA97LZPkkrBoZyMduTwmXt7u7JzetvnF0tHEIQ+T7QMmVjlD/a6P+IT92b2NG+h2n+Q5IRU\nEgOx7zFfJXvdH1HtP0gCDho5wWHf/uC1dtvHotxTFjM7J9LCiWxTT1s8EuELgtDnsUfObt3IXvdO\nfNpDonIyIWlqmOj7MJ78oIRhjHSMC1573F/FEd8BhidmMipxPEd8B0ggkXpdy98a3+D81CvC7mN1\nMmmkc7oz9BYQLbsHYttKPWnxiOALgtDnsWfIfObezqee0AROy3+3hBhtXZMSJrLnpXwpKNQAf2t8\ng2q/KQFW7S/n3YY1zEy9Khjdj3NO4YCnlHpdy0n/seBA8d+b/mLSO2nbwulpi0cEXxCEfsU45xS8\n2huM8AcnjKCk/nkyErPZ791FRsJYJjlnMCHpnLDr7J2GsVgUAOmM5CTHqNe17Gx6J/hWkKRSmZl6\nFR83bwpL5TziO8CoxPHBc+z3tI8z9MbMXBF8QRD6FUkqlcnJM4PfS+qfp17Xor3+oFUz0jEWIDjY\n6tZNQeEenDCcLzylVPsPMipxPH7tx+/34iQZL14+cW+mynuQ81MLGJwwPPhmkORMCYvYI0U80r7p\njYwdEXxBEPo101Ly2dFUwrSUfAYnDOULTymjHTlB68WrvVR4y6jXtfibfIx0jEPXJjNqyHjOSr4I\nt26msekkwxNG8/GRbaS7BlPtP8gXntKowm3ZR/bMHSuvf5Izr4WN05MZO5KlIwhCvyBWxsuIxDHk\nD7qBEYljghbLYe8+jvgOkJGQzTFfJfW6lkSSOKmP83DR95h93sVsL/2Acs+nlLm3Ua9r2bNnN3dd\n+Agv/cfGFsIdWfZ4r3snn7g3s9e9EzDR/aeebRz3HwmecyqLqHQWEXxBEPoFkamZsToAt26k2d9o\ncvATM6j2lzNIDcOHm6cfeYHnH11PTeUxHp7/c/Z+eoAjvgPUlNWz5OqHOVp5nKcf/QO/+OGvg/eL\nLtw6bDvOOSVoJ/VmcTWxdARB6BdEWiRWB1Djq+C8lC8FBfkLTymfe80MW7/2kpEwlrHOKXx3+f08\n/+j64P2qKo+y6Mq7WPbzb/LIkp9wtPJ48NjPf7gKn/ay9PtLGO3I4bB3X5hvPyHpnODMXYuhCSMZ\nmjCqV4urSYQvCEK/IDLSjhVVj3NOYZJzBhkJYzmqD1HtP8jWyhLWPVMSdr9CCmmqaOLfvlKMv1JR\nSGHY8T88vYb3KzfwcfOmNid9WZaOQzl6vFa+HYnwBUHolySp1Ba59db+bOckjvkOc1ri2SQlpDIy\nexz/8yp866oHqa48SiGF3MM9LGQhxRSznOXkkAPAGtYwKmsk/2/DSzCqlolJM3D5sqJOvGqrzENP\nIxG+IAh9lrZKE0ROyLLO+7h5E9X+gxzxH+CYrxKAcyafx9Ov/YyMzOFsZCP72EcOOaxmNTnksI99\nbGQjIzKH8txrqxh6ehJHfAeo8n4R9gxrwLY9ZR56GhF8QRD6LO1ZZeqk/xh/bfhTIGvmI8AUP0tT\n6TTqE1T7y9nRVMKnnm2oCXXc9bN/oZZaiikOu08xxdRSy0Mr7+X8M2cHOwprYDaUnql7fbHyWIjg\nC4LQZ4mWEhnJx82baNB1QKiOzuCE4YxxTAAglSFMcJ6Lk2SqPq3jp3c/wzCGsZzlYfdZznKGMYxH\nlvyEv+xaS7W/nIyEbMC8aYx25DAqcTzZzknBOjq9Vfc+FiL4giD0WdpjlVjRPBCsimn/PNaZy+ee\nv7Pv0wN8e/6POFp5nMu4LGjjLGJR0N65jMs4Wnmcu68uImVfJsMTM/nUs40vPKXB3P7D3n1R3zyi\n2U9SLVMQBKGLcOtGDnv3MSt1XjB10hpUHekYz3F/FdnOSRypOsx354dSL9ewBoCNbKSWWu7jPi7j\nsuD+yopKbrryDt7avoFBg4YxOGEEtb7DwQlZbt1ElfcLmnUjJ/3HOOzdh1d7wpZEhJ6vlikRviAI\n/RZLUA+4P6HGV8Ex3xHebVjDJ+7NlLm3BSPyC7ILuP6Wr4Rdu4Y1JGRqnnzpp6RlpQbF3qJw0dV8\nlvY+9bqWXc3vhqVdlnvKqPaX87lnRzBtE1QL+6k9llRXIhG+IAj9EvuatbW+wxzVhzjmOxJcovCs\n5ItwebOC6ZNLv7+EY/5DwclXo7JG8uxrT/Cls+dz6VlXkZ+fT0VFBQC3fvd6vv7gNRzVhwAYocZw\nmvNMxjmncNJ/jIMBK8dBMqc5p+JKzIpaUM2ypCJr73QXEuELghCXdNbfttasrfaXM8wxmlGJ4zkv\n5UuMShwfrGtvH1wd7cjh3ofv4P6H/43RWSN59c11TMo9A4Dc3FxKSkoYlTWSGx6cx9KH7wKlgs9q\noj4o1h83b6KRkwB4aWa/Z2eb4wztyTbqCiTCFwQhLumsv23VxQfNhKRzgoI72nFa1GcAHPEd4Krv\nzGDaraPxjDrCJ+4Dwefn5uby9x1/oyG9Cq/2cNRfGczSsVfPPCv5IvxNPlLUYJr0yWCt/Lbaat92\nFyL4giDEJZ0RQWtgNts5kcPefe1+Ro2vgiO+A0wcdVbQ8rEvUp6ZMRYYi1s3htXKsc/mHZwwnNlp\nCzrUXvsEse5EBF8QhLikNRFsa7UoK3Kv8n5Btb+cKu8XTE25pEWRs8gVqYYmjCRVDeGEr4YDnk+Y\nmDSdve6dfOrZymHvflyJWYx0jKfMvY2zki8K3qenFjDpLCL4giD0Odqye6xou1k3Uu0vp9pfzsfN\nm8LWmrV3GCf9x9jSuJ4GXUcqg2nkJEc9h0hWqVgzaY/6Kznqr2SfZycemqEZZqXOD+t8rLb15LKF\nHUEEXxCEPkdbdo9bN1Hjq2Bi0gwScQKabOekYFZOZIdhn42brAbRqE8yIiHTdn9Fjbeco7oSD82k\nqXQGq+FBsbePA7Q17tAba9laiOALgtDnaMvzDkbzbhOFWwxOGg607DCsgdb0hAzGJ50ZZv0Yv97B\nualzOOD+hDp/DemJLj737OCkPhYclI3Mr49Fb6xlayGCLwhCv+Os5IugmcCatC0j6sgOI3Kg1eoY\ngKCH79VekhNSqfYeZHjiaDISsjniO8BQz0gmJ88Knt+WiPdmqWTJwxcEod8xOGE4s1LnMzhheJs5\n7la+/0n/sah5/z7tAaDGV85oRw5nJs1mQtI5DE/MBOCY71CH5gr0ZqlkifAFQejXtBVR25dCtAZ1\no0XpR/2VHPbuCx6bkDSV4/4jwRW1+kKmjgi+IAh9mrYGQaP5/fZrrI5gtCMnOKhrP4fAhNqMhLEt\nVs6KtqKWWzcG6u4rJiRNjatsHRF8QRD6JJYgR6tC2RaRA6fWdUnOlKCAW+dMcs4IFjizi3esjsZa\nvxbAoRxxFfl3SvCVUo8B/wi4gc+ARVrr2sCxB4FbAR9wt9b6tU62VRAEIUhIkPM6XHEyls1jt3fs\n2TetTe6C8I7GlHTwACruVr3qbIT/BvCg1tqrlPpP4EHgO0qps4DrgbOBLGCDUipXa+3r5POE7qa6\nGlavhkWLICOjt1sjCDGxi3ZHbZNYaZ3jnFOCXr7LmxV2TmREH6vTSFKpYVk78USnsnS01q9rrb2B\nr5uBsYHPC4EXtNbNWuu9QBkQn38BIZzVq+H++81WEOKYrsx2sTJ1AM5L+VLUN4bIbJ94WZi8I3Sl\nh38L8GLgczamA7A4GNgnxDuLFoVvBWEAEMvTt9Ob+fNdRZuCr5TaAIyJcughrfWfA+c8BHiB5zra\nAKXUYmAxwPjx4zt6udDVZGTAsmW93QpB6FH6g5i3hzYFX2td0NpxpdTNwDXA5VprHdhdDoyznTY2\nsC/a/VcBqwDy8vJ0tHMEQRC6k/aUJ+7NkghdRac8fKXUVcD9wAKtdYPt0FrgeqVUslJqAjAJ2NKZ\nZwntoLoaHnvMbHvyWkEYAPT0+rPdQWdLK6wEhgBvKKX+rpT6FYDWehfwB+Bj4FXgTsnQ6QE6M+Aq\ng7WC0Cp9cZA2kk4N2mqtJ7Zy7BHgkc7cX+ggnRlwlcFaQej3SPG0/oQ14NqR/HnLyoGOXysI/YjO\nLpreFxDBH+hEWjni5QsDlLaqavYHRPD7A7FEuj3ivWABzJtnttA+L7+1+0qHIfRR+sOgbFuI4PcH\nYom0tX/lytgivHYtrF9vttXVUF8Py5eHvPxoAt5apyCDv0IfpT8MyraFVMvsD8QacF20CKqq4Lnn\noKzM7IucVGW/dvVqKC6GFStCXr4l4PZrWxvglcFfQYhbVGiuVO+Tl5ent27d2tvNiA/aKmLW3iJn\n8+ebCH7iRPj612HJktjnR7unFFMThLhHKbVNa53X1nli6cQrbVkj7bVOHn4YpkyBK6800Xtr52dk\nhCL93bvbl70jnr0g9BnE0olX2rJGYh3fvRuWLoXHH4fJk+Gdd6C0FK67znjz9fVGnGMJuNWRrFpl\nbKD6eigqit3OaJaPIAhxiQh+vNJaEbPWbJYlS2DDBnC74Y03Wnr0998PgwaF39t+v0WLYONGYwO1\nB/HsBaHPIILfF2ktqp4+3Qj+9EBxJ3vHYYnyJZcYb//hh80bQH29sXus+z37rMnsAbjhBmPZxPLw\npbqmIPQZRPD7IgsWmCjcyp23c+utsN0s5NDCurHE2RrI3bPH2DZz5oSnYmZkhGycxx4Ty0YQ+gky\naBvvRBsUtefOR/LUUybCf+yx2AO0jz9uJltdeaX5/tZbxuaJFsEvWmTSNBcsMPe0BnNlkFYQ+hwS\n4cc7Hc2D3xKoQn366eHHq6vDbZrLLjMinpEBDQ3RB3Pt3r7VjtdfNx1KW4O5giDEHSL48U40+yaa\nb25l55xxhonYv/KV8Fz6m24KDcS++KLJ3AEj2pZtEzmYa+9srHaMCSx+tnFj69k+giDEHWLp9DQd\nzVtvzb6x32vpUnPegQPGgvnyl41Xv3u3Ee71641XP2GCEfuCAhO5W+UUli0zs3KLikJts+rsXHJJ\n+P3BdCqRlpHk5AtCXCMRfk/T0bz11uwb614bN5qMGzD+PBiRrqoyn5991mzt2TjTpxuL5913jUUz\nb17oDWDQIPM8S+TdbnNOQYG55vnno7dJcvIFIa4Rwe9pOpq33lraoz1n/rLLYN06s3/+fCP2Lhec\nfbbZt2xZy8jbEv9580xHMXNm6L4rV5r7FhSEUj3/4R/MZK5Y3r31O1kDvFKOQRDiChH8nqYzeeuR\nE64yMkz0bu2zzjn7bBOVT59uhDctjRrAtWSJEevdu80ErTvvpCYtDdf995t7LVnS0qaZPh3S0kza\n5pIl7WubRPqCEJeI4PclLCGtrw/ZLlYHYvnn9fVmu3y5uWbZMop+8xue3LuXktdfJ3fuXPNW8NZb\n7KmvJ3//fm5LS6OoqChcqJcsMc+wbCCrgmasWb72a2X2rSDEJSL4fQlLQOvrw8V19Wpj4Tz2mBH/\nFSuCQl00cSLFe/cCkP/ee5S89x65l17KHiD/ww+pqK+nOGDtFFkR/IIF4W8NVucCsaN3u8jL7FtB\niEtE8OOZaBaOFc1bImwJcEFB6Lr6emhooGjOHIrfeiu4uwLIT0xk1Y03snj7dipOnAgeC4q+PU0T\nzPOiCfuCBSEv3yq5LCIvCHGNCH480x4v3C7Azz8fzLqpAZ5MTg47tZBCNvo2cs3ttzOMYRRSyBrW\nBGJxRzMAAAwCSURBVI8/+eST3HXjjbiqqkwKZ1VV7Fz7558PDfpG5u8LghCXiODHM7G88MiOwIr0\nwWTTTJyIKyuLkrffJj8piQq3m0IKuYd7WMhCiilmOcvJIQeANawha/BgSkpKcP3616E6+G+9Bbt2\nmQyetWvD3yiWLw+NE4hXLwh9AhH8eCbSJrEsHmvWrWWrvP46vPce3HlnKJ8+O5vcggJKRo8m/7nn\n2MhGFrKQHHJYjekc9rGPjWwkKz3dDOjm5oYKr02YYGbtrl8Pn38empkb6dULgtBnkJm2fQkrurai\n7aVLja3y3nvm+EcfmVz6OXNMdL5hA7kTJ7Jq0SJqqaWY4rDbFVNMLbWs+spXyD3jDLNz5UrTabzy\nirFt5s0zYj9njulYampaXwFLEIS4RQS/L2FVrrRPjLr0UhPZFxTAOeeYDsBap/jSS9lz8CCL/+//\nGMYwlrM87HbLWc4whrF49Wr2LFhg3iAmTzYTuFwu08E8/LARfa2NXbR0adf8LlKGQRB6HBH8voAl\njtAyus7PN+L/xhuh/RdcACtWsGfSJPKfeoqKEye4jMvIIYd97GMRi9jHPnLI4TIuM9k7773Hnq98\nJSTAK1eat4lvfMN0LEqZTuXhh42NZK+5cyq0d01eQRC6Dq113PzMmDFDC1FYvlxrMFuLqiqtV6zQ\nurTUbKuqzOd587QuLdXV1dU6a8gQDQR/CinUwximAT2MYbqQwrDjWaCrrWcsW2aeCVpPnGi2K1aY\nH2v/ihWn/jtZ7a+q6sQfRhAErbUGtup2aKxE+H2Jd98NRdXWgO7ataFI2VZZ0+Vycdvtt4ddvoY1\npFHLy0OGkDbEF5aSCXBbQYEpvwCmnIJFVlZoRawFC0ykv2xZqGbOqUT6VvtlLEAQegzJ0ukLLFkC\nH3xgxHz16tgzXCP2Fa1YAUBxwA7KSkujJDmZ3GPHKAHylaIi4PcvX7aMorQ0Y+XccIO5z7JlJmtn\nwwYYPNgM2C5dar7PnRvqbKxzBUGIa0Tw+wL2ImmtVaK0p3EGVrgqSkuDO+/kyd/+1qRennEGrFhB\n7tNPU1JTQz5w20UXUTRyZEi8rc5lxQrTARQWhkonr19vBnEjO5hYNXYEQYgbxNLpK0SzcCD24OfK\nlSZjp7iYoowMPrzvPiP2GRkwciTU1JDrcvHhrbdS9PTTZlbthRearJ+HHw5lA61da9IyrRLKK1aY\nzidS1GUQVhDiHonw+xqRFk7k1oq0GxrM90CNHVdxsYncn302rI6+q7LS5NtbWUAA77wTff3cyIlg\nUiFTEPoUIvh9GfvM28ha9MuXm2jcqrFTUGDsmJUrTe0ba2Ws9evNZK3ly00nkZYWW7QjbRupkCkI\nfQoR/L6GfVnDmTONbWOtegWhxcZvuCG0OlVxsRHjuXONdVNcbCpq2hdPieW7W9bQ66+bFa+sgmlW\nhk200g/i4wtCXCKC39ewL2vodpvI/IYbjPjX15to3oraBw0KWTtpaUac7csTdiQq37DBrH41b16o\nlk8kstKVIMQ1nRq0VUr9QCn1oVLq70qp15VSWbZjDyqlypRSu5VSV3a+qQIQytiZN8+I8AcfmDII\ngwaFou9AOmZQfOfNC6VaLllijlv59m2VOFiyJFQZMy0tmOcfFXvph2hIOQVB6F3aMzsr1g+Qbvt8\nN/CrwOezgB1AMjAB+AxIbOt+MtO2A1RVaV1QEJqBGzlz1fpuzdKNNiu2qsrMzI12PNpM2M7OjrVm\n6XZmhq4gCC2gnTNtO2XpaK3rbF8HYaboAywEXtBaNwN7lVJlwCzgvc48TwhgeeXTp5soH1raM/bV\nscDYPdZn+wBvtLx665xIe6azA7OSySMIvUqnPXyl1CPAvwDHgfzA7mxgs+20g4F9QlcQmYlzySUw\nf77JvJk8OfzcjAxj99x/v9lC9FTKyEHW7hBnyeQRhF6lTcFXSm0AxkQ59JDW+s9a64eAh5RSDwJL\nIKIGb9v3XwwsBhg/fnxHLh24RAr1/PmhLJ1161o/376vNQEWcRaEfofSVu30zt5IqfHAeq311ID4\no7V+NHDsNaBIa92qpZOXl6e3bt3aJe0ZUGzeHLJoZs8OP9aRVMlY50q6pSDENUqpbVrrvLbO62yW\nziTb14VAYB081gLXK6WSlVITgEnAls48S2iFd94x5Q/eeaflsY6UPIh1rpRNEIR+QWc9/B8ppSYD\nfmA/cDuA1nqXUuoPwMeAF7hTa+3r5LOEWFhWTbTCah3x4hctMoO71gCv/R7R9guC0KfoVISvtf6K\n1nqq1vpcrfU/aq3Lbcce0VqfobWerLV+pfNNFWISq7Baa0TLibcGeIuLQ/ew7BwI3y8IQp9DZtr2\nJ6JF87Fmv8baH3mPyIwgSakUhD6LCH5/IlpmTSxLJ9b+yHu0lropCEKfosuydLoCydLpJJJNIwgD\nkh7J0hHiDMmmEQShFcTS6U90xexYeUsQhH6LRPj9Cct/74xQy1uCIPRbJMIXwpECZ4LQbxHBF8KR\nGjqC0G8RS0cQBGGAIIIvCIIwQBDBFwRBGCCI4AuCIAwQRPAFQRAGCCL4giAIAwQRfEEQhAFCXBVP\nU0pVAfVAdVvnxhkZ9L02g7S7J+mLbQZpd0/SmTafprUe2dZJcSX4AEqpre2p+hZP9MU2g7S7J+mL\nbQZpd0/SE20WS0cQBGGAIIIvCIIwQIhHwV/V2w04Bfpim0Ha3ZP0xTaDtLsn6fY2x52HLwiCIHQP\n8RjhC4IgCN1AXAm+UurflFJaKZVh2/egUqpMKbVbKXVlb7YvEqXUD5RSHyql/q6Uel0plWU7Fs/t\nfkwpVRpo+0tKqWG2Y3HZbqXUV5VSu5RSfqVUXsSxuGyzhVLqqkDbypRSD/R2e2KhlHpaKXVEKbXT\ntm+EUuoNpdSnge3w3mxjJEqpcUqpEqXUx4F/H/cE9sd7u1OUUluUUjsC7S4O7O/edmut4+IHGAe8\nBuwHMgL7zgJ2AMnABOAzILG322prc7rt893Ar/pIu+cCjsDn/wT+M97bDZwJTAY2Anm2/XHb5kD7\nEgNtOh1ICrT1rN5uV4y2XgqcD+y07VsBPBD4/ID1byVefoBM4PzA5yHAnsC/iXhvtwIGBz47gfeB\n2d3d7niK8P8buB+wDyosBF7QWjdrrfcCZcCs3mhcNLTWdbavgwi1Pd7b/brW2hv4uhkYG/gct+3W\nWn+itd4d5VDctjnALKBMa/251toNvIBpc9yhtX4bOBqxeyHwbODzs0BhjzaqDbTWlVrrvwU+nwA+\nAbKJ/3ZrrfXJwFdn4EfTze2OC8FXSi0EyrXWOyIOZQNf2L4fDOyLG5RSjyilvgC+Dnw/sDvu223j\nFuCVwOe+1G6LeG9zvLevLUZrrSsDnw8Bo3uzMa2hlMoBpmOi5bhvt1IqUSn1d+AI8IbWutvb3WNL\nHCqlNgBjohx6CPguxmaIO1prt9b6z1rrh4CHlFIPAkuA5T3awBi01e7AOQ8BXuC5nmxbLNrTZqH3\n0FprpVRcpvUppQYDfwTu1VrXKaWCx+K13VprH3BeYAztJaXU1IjjXd7uHhN8rXVBtP1KqXMw3uuO\nwH+kscDflFKzgHKMt28xNrCvx4jV7ig8B6zHCH7ct1spdTNwDXC5DhiG9HK7O/C3ttPrf+s2iPf2\ntcVhpVSm1rpSKZWJiUbjCqWUEyP2z2mt/xTYHfftttBa1yqlSoCr6OZ297qlo7X+SGs9Smudo7XO\nwbzynq+1PgSsBa5XSiUrpSYAk4AtvdjcMJRSk2xfFwKlgc/x3u6rMOMlC7TWDbZDcd3uGMR7mz8A\nJimlJiilkoDrMW3uK6wFbgp8vgmIqzctZaLEp4BPtNaP2w7Fe7tHWtlxSqlU4AqMfnRvu3t7tDrK\n6PU+Alk6ge8PYbIcdgNX93b7Itr6R2An8CHw/4DsPtLuMoyv/PfAz6/ivd3AtZhgoBk4DLwW7222\ntW8eJnvkM4w91ettitHO54FKwBP4W98KuIA3gU+BDcCI3m5nRJsvxgx2fmj79zyvD7T7XGB7oN07\nge8H9ndru2WmrSD8/3bsmAAAAIZBmH/Xc9EdJCI4gIj3pQPAhuADRAg+QITgA0QIPkCE4ANECD5A\nhOADRBw3AKVnPeNvWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f229dc9b5c0>"
      ]
     },
     "metadata": {},
     "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.\n",
    "\n",
    "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",
    "![Gaussian](http://images.books24x7.com/bookimages/id_5642/fig11-10.jpg)\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": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy import exp, sqrt, array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gaussian(d, bw): return exp(-0.5*((d/bw))**2) / (bw * math.sqrt(2*math.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " This person at the science march certainly remembered!\n",
    "\n",
    "<img src=\"images/normal.jpg\" width=400>\n",
    "\n",
    "Since all of our distances are positive, we'll only be using the right-hand side of the gaussian. Here's what that looks like for a couple of different choices of bandwidth (bw)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW5+PHPM5M9kBCyEJIQAiTsOxFEEEWRVa/rbd3r\ntS0Xt9uKdentru2vtl4V11qkWm211qooLoC4ggpC2GSHAAESlkDYCdmf3x9ngEiBDDCTk8w879fr\nvM6cbeY5UZ7vOd/zPd+vqCrGGGPCh8ftAIwxxjQuS/zGGBNmLPEbY0yYscRvjDFhxhK/McaEGUv8\nxhgTZizxG2NMmLHEb4wxYcYSvzHGhJkItwM4kZSUFM3JyXE7DGOMaTYWLly4S1VT/dm3SSb+nJwc\nCgoK3A7DGGOaDRHZ5O++VtVjjDFhxhK/McaEGUv8xhgTZppkHb8xJjxVV1dTXFxMRUWF26E0WTEx\nMWRlZREZGXnG3+FX4heR0cATgBeYoqoPn2S/c4C5wLWq+sbpHGuMMcXFxbRs2ZKcnBxExO1wmhxV\npaysjOLiYjp06HDG39NgVY+IeIFngDFAd+A6Eel+kv3+AHx4uscaYwxARUUFycnJlvRPQkRITk4+\n6zsif+r4BwKFqrpBVauA14DLT7DfXcCbQOkZHGuMMQCW9BsQiL+PP1U9mcCWesvFwKDjAskErgSG\nA+eczrGB9OTH60huEUXnNi3JS2tBq7ioYP2UMcY0W4F6uDsJuF9V6860NBKR8cB4gOzs7NM+vqa2\njilzNrC/ouboupQW0XRu04K8tBacl5vCJd3a4PHY1YQx5sSKioq49NJLWb58eUC/t6ysjGuuuYYF\nCxZwyy238PTTTwf0+0+XP4m/BGhXbznLt66+fOA1X9JPAcaKSI2fxwKgqpOByQD5+fmnPQJ8hNfD\nkl+OZNv+CtbuOEDhjoOsKz3A2h0HeXNRCS/N3UTnNi2486I8xvVqi9cKAGNMI4mJieGhhx5i+fLl\nAS9UzoQ/dfwLgDwR6SAiUcC1wLT6O6hqB1XNUdUc4A3gdlV9259jA8njETJbxTK8Sxo/HNaRP17T\nh7fvGMLSX43kiWv7ogr/84/FXPL450xdXExNbV2wQjHGNFM1NTXccMMNdOvWjWuuuYbPP/+cq666\nCoB33nmH2NhYqqqqqKiooGPHjn59Z3x8PEOHDiUmJiaYofutwSt+Va0RkTuBmThNMl9Q1RUiMsG3\n/bnTPTYwofvP6xEu75vJZb0zmLFiO09+vI67/7mUSR+t48cj8riib6Y9UDKmifnNuytYuXV/QL+z\ne0YCv7qsxyn3WbNmDX/5y18YMmQIt956KwsWLGDJkiUAzJkzh549e7JgwQJqamoYNMh5ZPnII4/w\nyiuv/Nt3DRs2jCeffDKg5xAIftXxq+oHwAfHrTthwlfVWxo61i0ejzC2V1tG90jno1U7ePITpwD4\nqrCMh67oSUyk1+0QjTEua9euHUOGDAHgxhtv5Mknn6RTp06sWrWK+fPnM3HiRGbPnk1tbS3nn38+\nAPfeey/33nuvm2GflrB8c9fjEUb2SGdEtzZM+mgtT35SyNrSgzx3Y3/aJsa6HZ4xBhq8Mg+W4+/+\nRYRhw4Yxffp0IiMjGTFiBLfccgu1tbU88sgjQIhe8Ycqj0eYOLILPTITmfjPJVz21Bc8e8MABnZo\n7XZoxhiXbN68mblz5zJ48GBeffVVhg4dyoABA7j55pu5+eabSU1NpaysjB07dtCzZ0+g+V3xWydt\nwKge6bx9xxBaxkRy/fPz+NvcIlRPu2GRMSYEdOnShWeeeYZu3bqxZ88ebrvtNgYNGsSOHTsYNmwY\nAL1796ZXr16n9WwwJyeHiRMn8te//pWsrCxWrlwZrFNokDTFBJefn69uDMSy73A1d/9zCZ+sLuU7\n+Vn87speRHqtbDSmsaxatYpu3bq5HUaTd6K/k4gsVNV8f463rFZPYmwkU27O587hubxeUMx9b3xD\nXV3TKxiNMeZshHUd/4l4PMJPRnUhOsLDo7PWkhQXxS8u7WbNPY0xIcMS/0nceVEuu8ureOHLjSS3\niOKO4bluh2SMMQFhif8kRIRfjOvO3vJqHpm5hlZxkdwwqL3bYRljzFmzxH8KHo/wx2t6s+9wNT9/\nezlJcVGM7dXW7bCMMeas2MPdBkR6PTxzfX8GZCfxo9cW88W6XW6HZIwxZ8USvx9io7z85ZZz6JTa\ngvF/K2B5yT63QzLGBEFRUdHRl7ICadasWQwYMIBevXoxYMAAPvnkkxPu9+tf/5rMzEz69u1L3759\n+eCD4PR2Y4nfT4mxkbx860ASYyO5/ZVF7K+odjskY0wzkZKSwrvvvsuyZct46aWXuOmmm0667913\n382SJUtYsmQJY8eODUo8lvhPQ1pCDE9d14+SvYf56ZvL7O1eY0JQMLpl7tevHxkZGQD06NGDw4cP\nU1lZGbRzaIg93D1N+TmtuWdkZ/44Yw3nfp3MTedaSx9jgmL6A7B9WWC/M70XjHn4lLsEu1vmN998\nk/79+xMdHX3C33/qqad4+eWXyc/P59FHHyUpKelMzvSULPGfgQnDOvH1ht089N5K+me3okdGotsh\nGWMCJJjdMq9YsYL777+fDz/88ITbb7vtNn7xi184zcl/8QvuueceXnjhhcCdnI8l/jPg8QiPfacP\nY5+cw52vLubdu4bSItr+lMYEVANX5sESrG6Zi4uLufLKK3n55Zfp1KnTCX+7TZs2Rz//8Ic/5NJL\nLw3UaX2LX3X8IjJaRNaISKGIPHCC7ZeLyDciskRECkRkaL1tRSKy7Mi2QAbvpuQW0Tx5bT82lR3i\nf9+y+n5jQsWRbpmBo90yn3/++UyaNInBgwcf7ZZ5zZo13+qW+cgD2frTkaS/d+9exo0bx8MPP3z0\nbuJEtm3bdvTz1KlTg9LCCPxI/CLiBZ4BxgDdgetEpPtxu30M9FHVvsCtwJTjtg9X1b7+9hzXXAzq\nmMzdIzozbelWXluwxe1wjDEBEIxumZ9++mkKCwt58MEHjzbVLC0tBeAHP/gBR3ojvu++++jVqxe9\ne/fm008/5fHHHw/KOTbYLbOIDAZ+raqjfMs/BVDV359i/xdUtZtvuQjIV1W/33xyq1vmM1Fbp3zv\nhfksKNrN23cMoVvbBLdDMqbZsm6Z/dMY3TJnAvUvZ4t9647/0StFZDXwPs5V/xEKfCQiC0Vk/Ml+\nRETG+6qJCnbu3OlP7E2C1yM8/t2+tIyJZOLrS6murXM7JGOMOaWAteNX1amq2hW4Anio3qahviqg\nMcAdIjLsJMdPVtV8Vc1PTU0NVFiNIrVlNL+9oiertu3n+Tkb3A7HGGNOyZ/EXwK0q7ec5Vt3Qqo6\nG+goIim+5RLfvBSYCgw842ibsNE90xnTM51JH61jw86DbodjTLNlDSVOLRB/H38S/wIgT0Q6iEgU\ncC0wrf4OIpIrvqccItIfiAbKRCReRFr61scDI4HlZx11E/Wb/+hBTISHB95aZiN3GXMGYmJiKCsr\ns+R/EqpKWVkZMTExZ/U9DTY+V9UaEbkTmAl4cR7crhCRCb7tzwFXAzeLSDVwGPiuqqqItAGm+sqE\nCOBVVZ1xVhE3YWkJMfxsXDfuf3MZry3YwvWDst0OyZhmJSsri+LiYprTc77GFhMTQ1ZW1ll9hw22\nHmCqyvXPf83ykn3MmngB6YlnVzIbY4w/bLB1F4kIv7+qF1W1dfzineV2y2qMaXIs8QdBTko8Ey/p\nzKyVO5i+fLvb4RhjzLdY4g+S7w/tQM/MBH75zgr2lle5HY4xxhxliT9IIrwe/nB1b/aUV/H/Pljl\ndjjGGHOUJf4g6pGRyA+GduD1gmKWbNnrdjjGGANY4g+6uy7OI7VlNL+etsLa9htjmgRL/EHWIjqC\n+0Z1YcmWvbyz9KQvPBtjTKOxxN8Iru6fRe+sRB6evppDlTVuh2OMCXOW+BuBxyP86rIe7NhfybOf\nFbodjjEmzFnibyQD2idxZb9Mnp+zkc1l5W6HY4wJY5b4G9H9o7viFbHmncYYV1nib0TpiTHcMbwT\nM1Zs56tCvwckM8aYgLLE38h+cH5HspJi+c27K6mx0bqMMS6wxN/IYiK9/HxcN9bsOMA/5m92Oxxj\nTBiyxO+CUT3SGdwxmUdnrWXf4Wq3wzHGhBlL/C4QEX5+aTf2llfz3Ofr3Q7HGBNm/Er8IjJaRNaI\nSKGIPHCC7ZeLyDciskRECkRkqL/HhqseGYlc0TeDF77YyLZ9h90OxxgTRhpM/CLiBZ4BxgDdgetE\npPtxu30M9FHVvsCtwJTTODZs3TOyC6owadY6t0MxxoQRf674BwKFqrpBVauA14DL6++gqgf12FBT\n8YD6e2w4a9c6jhvPbc+/Fm5h3Y4DbodjjAkT/iT+TGBLveVi37pvEZErRWQ18D7OVb/fx/qOH++r\nJioIp4GW77wol7ioCP44c43boRhjwkTAHu6q6lRV7QpcATx0BsdPVtV8Vc1PTU0NVFhNXuv4KCZc\n0JFZK3dQULTb7XCMMWHAn8RfArSrt5zlW3dCqjob6CgiKad7bLi6dWgH0lpG8/vpq21wdmNM0PmT\n+BcAeSLSQUSigGuBafV3EJFcERHf5/5ANFDmz7EG4qIi+PGIzizctIdZK3e4HY4xJsQ1mPhVtQa4\nE5gJrAJeV9UVIjJBRCb4drsaWC4iS3Ba8XxXHSc8Nhgn0tx9Jz+Ljinx/HHmGuvKwRgTVNIUqxby\n8/O1oKDA7TAa3Yzl25jw90X84epefPecbLfDMcY0IyKyUFXz/dnX3txtQkb1SKdvu1Y8Nmsth6tq\n3Q7HGBOiLPE3ISLCA2O6smN/JX+bV+R2OMaYEGWJv4k5t2My5+el8KfP1nOgwjpwM8YEniX+Jugn\nI7uwp7yaF78scjsUY0wIssTfBPVp14qR3dvw/OwN7C2vcjscY0yIscTfRE0c2ZmDVTVMnr3B7VCM\nMSHGEn8T1TU9gct6Z/Dil0XsPFDpdjjGmBBiib8J+/GIPKpq6/jTZzZYizEmcCzxN2EdU1twdf9M\n/v71JhusxRgTMJb4m7j/uTgPVeWpTwrdDsUYEyIs8TdxWUlxXDcwm9cXbGFzWbnb4RhjQoAl/mbg\nzuG5RHiFSR+vdTsUY0wIsMTfDKQlxPC9wTm8vbiEwlIbotEYc3Ys8TcT/31BJ2IjvTz+kQ3Mbow5\nO5b4m4nW8VHcOrQD73+zjdXb97sdjjGmGfMr8YvIaBFZIyKFIvLACbbfICLfiMgyEflKRPrU21bk\nW79ERMKvk/0A+sHQjrSMjmDSLLvqN8acuQYTv4h4cUbVGgN0B64Tke7H7bYRuEBVe+EMtD75uO3D\nVbWvv4MEmBNLjIvk++d3YMaK7Swv2ed2OMaYZsqfK/6BQKGqblDVKuA14PL6O6jqV6q6x7c4D2dQ\ndRMEtw7tQEJMBJM+shY+xpgz40/izwS21Fsu9q07me8D0+stK/CRiCwUkfGnH6KpLyEmkvHDOvLR\nqlKWbtnrdjjGmGYooA93RWQ4TuK/v97qoaraF6eq6A4RGXaSY8eLSIGIFOzcuTOQYYWcW4Z0ICku\nksftqt8Ycwb8SfwlQLt6y1m+dd8iIr2BKcDlqlp2ZL2qlvjmpcBUnKqjf6Oqk1U1X1XzU1NT/T+D\nMNQiOoL/vqATn63ZycJNexo+wBhj6vEn8S8A8kSkg4hEAdcC0+rvICLZwFvATaq6tt76eBFpeeQz\nMBJYHqjgw9nNg9uTHB/F47Psqt8Yc3oaTPyqWgPcCcwEVgGvq+oKEZkgIhN8u/0SSAaePa7ZZhvg\nCxFZCswH3lfVGQE/izAUFxXBbRd24ovCXXy9oazhA4wxxkdU1e0Y/k1+fr4WFFiT/4Ycrqpl2COf\n0ik1ntfGD3Y7HGOMi0Rkob9N5u3N3WYsNsrL7Rd2Yt6G3XxVuMvtcIwxzYQl/mbuuoHZpCfE8Ois\ntTTFuzdjTNNjib+Zi4n0csdFuSzctIfZ6+yq3xjTMEv8IeA7+VlktorlsQ/X2FW/MaZBlvhDQHSE\nl7suymVp8T4+WV3qdjjGmCbOEn+IuHpAFtmt43jM6vqNMQ2wxB8iIr0efnRxHiu27mfmiu1uh2OM\nacIs8YeQy/tm0DElnsdnraOuzq76jTEnZok/hER4PfxoRB5rdhzg/WXb3A7HGNNEWeIPMZf1zqBz\nmxZM+mgttXbVb4w5AUv8IcbjEe4e0Zn1Ow8xbem/daJqjDGW+EPRqB7pdGubwBMfraOmts7tcIwx\nTYwl/hDk8QgTL+lMUVk5by2yq35jzLdZ4g9RI7ql0ScrkSc+XkdlTa3b4RhjmhBL/CFKRJg4sgsl\new/zzwVbGj7AGBM2LPGHsGF5KQzs0JqnPinkcJVd9RtjHH4lfhEZLSJrRKRQRB44wfYbROQbEVkm\nIl+JSB9/jzXBIyLcO6oLOw9U8tLcIrfDMcY0EQ0mfhHxAs8AY4DuwHUi0v243TYCF6hqL+AhYPJp\nHGuC6Jyc1gzvksqfPlvP/opqt8MxxjQB/lzxDwQKVXWDqlYBrwGX199BVb9S1T2+xXlAlr/HmuC7\nZ2QX9h2uZsrsDW6HYoxpAvxJ/JlA/aeDxb51J/N9YPoZHmuCoGdmIuN6tWXKFxvZdbDS7XCMMS4L\n6MNdERmOk/jvP4Njx4tIgYgU7Ny5M5BhGeDuSzpTUV3Lnz5b73YoxhiX+ZP4S4B29ZazfOu+RUR6\nA1OAy1W17HSOBVDVyaqar6r5qamp/sRuTkNuWguu7p/F3+ZtYtu+w26HY4xxkT+JfwGQJyIdRCQK\nuBaYVn8HEckG3gJuUtW1p3OsaTw/GpGHqvLkx4Vuh2KMcVGDiV9Va4A7gZnAKuB1VV0hIhNEZIJv\nt18CycCzIrJERApOdWwQzsP4ISspjhsGtef1gi0U7TrkdjjGGJdIUxymLz8/XwsKCtwOIySVHqhg\n2B8/ZVSPdJ64tp/b4RhjAkREFqpqvj/72pu7YSatZQz/NaQD05ZuZcXWfW6HY4xxgSX+MDThgk4k\nxkby8PTVbodijHGBJf4wlBgbyZ3Dc5mzbhdz1lnTWWPCjSX+MHXT4PZktorl4emrbWB2Y8KMJf4w\nFR3h5d5RXVixdT/Tlm51OxxjTCOyxB/G/qNPBj0yEnhk5hobrMWYMGKJP4x5PMJPx3SjZO9h/jZ3\nk9vhGGMaiSX+MDc0L4Xz81J46pNC9pVbt83GhANL/IYHxnRlf0U1z35uXTkYEw4s8Rt6ZCRyZd9M\nXvyyiJK91oGbMaHOEr8BYOLIzgA89uHaBvY0xjR3lvgN4HTgdst5Oby1uJiVW/e7HY4xJogs8Zuj\n7rgwl1axkTz43gqaYud9xpjAsMRvjkqMi2TiyC7M27CbmSu2ux2OMSZIItwOIKAm9YaaSvB4Qbwg\ncuyzNxIiYnxT9LF5ZBxEt4DolhDlmx+ZYpMgtjXEtXY+eyPdPsOgu+6cdrwybxO/fX8VF3ZJIybS\n63ZIxpgAC63E33kU1FRAXR1oHWgt1NU689pqp1CoqYCqQ1Be5ixXl0PlAag6CHU1p/7+6ASnAGiR\nBvFpzrxFm2PzhLaQkAXxqeBpnjdTEV4Pv7y0O9dP+Zq/fLGRO4bnuh2SMSbA/Er8IjIaeALwAlNU\n9eHjtncFXgT6Az9T1f+rt60IOADUAjX+DhRwRsY+cubHqjqFQuVBqNzvTIf3QPluZ3708244WAp7\nimDL104BwnH14Z5ISMiAxCxIyIRW7SApxzd1cLZ5mu6V9Hm5KYzq0YZnPi3k6v5ZpCfGuB2SMSaA\nGhyBS0S8wFrgEqAYZxzd61R1Zb190oD2wBXAnhMk/nxV3eVvUM1qBK7aaji0Cw5uh/3bYH8J7NsC\n+0p8n31zrdcXjicSWmU7BUFKHiTnOlNKnlNQiLh2OkdsLitnxGOfc2nvtjz23b5uh2OMacDpjMDl\nzxX/QKBQVTf4vvw14HLgaOJX1VKgVETGnUG8zZs30lfF0xYyTjKUYW017Ct27hLqT7s3wOZ5UF1v\n/NvIOKcQSO0Kad2OTYnZjVp9lJ0cxw/O78Czn63nxsHt6Z+d1Gi/bYwJLn8Sfyawpd5yMTDoNH5D\ngY9EpBb4s6pOPo1jQ4M3Elp3cKbjqcKBbbBrHZStg7L1sGstbPoKlr1+bL/IeEjtAm16QHpvSO/l\nfI5JCFrYtw/P5Y2Fxfzm3ZVMve08PB7370SMMWevMR7uDlXVEl910CwRWa2qs4/fSUTGA+MBsrOz\nGyGsJkLEqfNPyICOF3x7W8U+2LkGSldC6WooXQGr34fFfzu2T1KOUwi07ePccbTtB/HJAQmtRXQE\n94/uyj3/WsrUxSVcPSArIN9rjHGXP4m/BGhXbznLt84vqlrim5eKyFScqqN/S/y+O4HJ4NTx+/v9\nIS0mEdoNdKYjjtwhbF8G27+B7cud+ap3j+3TKtspBDL6QUZ/yOzvNE89A1f2y+TleZv4w4zVjOqZ\nTovo0GoIZkw48udf8QIgT0Q64CT8a4Hr/flyEYkHPKp6wPd5JPDgmQZr+PYdQudRx9ZX7INtS2Hr\n4mPTyneOHARp3SErH7LOcaaUzn49M/B4hF9d1p2rnv2Kpz5ex0/HdgvOeRljGk2DiV9Va0TkTmAm\nTnPOF1R1hYhM8G1/TkTSgQIgAagTkR8D3YEUYKo4rVQigFdVdUZwTiXMxSRCh2HOdET5bti6CIoL\noHiBUxAsesnZFp0I7c6B7HOh3bmQOQCi4k741f2zk/hOfhZ/+WIjV/bPpGt68J4rGGOCr8HmnG5o\nVs05m5O6Oti93ikItnztTKW+xlmeCOc5Qbtzof15zhTX+uihew5VcfFjn5OTHMcbE+xBrzFNzek0\n57TEH+7Kdzt3A5vnOVPJQqitBMRpNdR+COQMgfZDeGN1BT/511L+35W9uH5QGD2AN6YZCHQ7fhPK\n4lo7zwqOPC+oqXSSf9GXsOkLpwXR/D8DcHVad+KT85g+fTEjO44nJTXNxcCNMWfKrvjNqdVUwbYl\nUPQFFM2hbtNXeGoqqMODJ6OP77nCBU7VUGSs29EaE7asqscET00l/3z7bbYunsmtGZtJLFsCddXg\njXYeFHe8EDoNh/Q+zbajOmOaI0v8JqgqqmsZ88QcVJUZtw8gZut82PAprP/UeckMnO6sO14AuSOg\n08VOlxbGmKCxOn4TVDGRXn57RU9umPI1z365jYkjR0DeCGfjgR2w8XOnEFj/CayY6qxP6w65FzuF\nQPvznLEQjDGusCt+c8bu/ucS3vtmK9N/dD65aSd4M1gVdqyA9R9D4ceweS7UVjkd0XUYBnmXQO4l\nkNS+8YM3JsRYVY9pFLsOVnLxo5/TuU0LXhs/GG9DbfurDjkPidfNgsJZTg+lACldnEIgbyRkD4aI\nqKDHbkyoscRvGs2bC4u5519L+dnYbvxwWEf/D1SFskKnEFj3IWz60rkbiGrpPBzuPMopCFpYk1Fj\n/GF1/KbRXNU/kxkrtvPIh2u4sEsqeW387AxOxBl4JiUPBt/u3A1s+BzWzYS1H8Kqac5+Gf197xmM\ndt4sbgKD1BjT3NkVvzlrOw9UMmrSbLKSYnnztvOI9J5lM05Vp/fRdTNh7UyniwkUWvo6pusyxnlG\nYO8NGHOUVfWYRvfBsm3c/soiJl7Smf+5OC+wX35wp1MdtHY6FH7ijFgWGee8M9BljHM3YFVCJsxZ\nVY9pdGN7teXyvhk8+fE6LuqaRs/MxMB9eYtU6HeDM9VUQtEcWDMD1s6ANR8A4nQ53WWsM6V2sSoh\nY07BrvhNwOwtr2Lk47NJioti2l1DiI7wBvcHVWHHclgz3RmZbNsSZ31SB+g6zikEss8FT5DjMKYJ\nsKoe45pPV5fyX39dwG0XduL+0V0b98f3lfjuAqY7L5HVVkFcslMV1HUcdBx+0jEHjGnuLPEbVz3w\n5je8XrCFf00YzID2rRs+IBgqD0DhR86dwNoPoXIfRMRCp4ucQqDz6ICNTWxMU3A6id+v5hciMlpE\n1ohIoYg8cILtXUVkrohUishPTudYE3p+Nq4bbRNjuef1pRyoqHYniOiW0ONKuHoK3Lcebnob+t/k\nVAe9czv8Xy789VKY9yfYs8mdGI1xSYNX/CLiBdYClwDFOGPwXqeqK+vtkwa0B64A9qjq//l77InY\nFX/zN3/jbq57fh5jeqbz1HX9kKbysFXVSf6r33emIyOQpfeCrpc6U5se9nDYNDuBvuIfCBSq6gZV\nrQJeAy6vv4OqlqrqAuD4y7sGjzWhaWCH1twzsjPvfbONv3+92e1wjhGBjH5w0c/h9rlw1yIY+VuI\njIfPHobnhsCTfWHmz2DTV1BX63bExgScP805M4Et9ZaLgUF+fr/fx4rIeGA8QHa2DesXCiYM68SC\njbt56N2V9M1qRa+sADbxDJTkTnDeXc50sNTXQug9mD8Z5j4NcSnOuwLdLnMGnImMcTtiY85akxkp\nQ1Unq2q+quanpqa6HY4JAI9HePQ7fUluEcUdry5iv1v1/f5qkQYDvgc3/Avu2wDXvOiMKbDyHXj1\nO/BIJ3j9e7DsDajY53a0xpwxf674S4B29ZazfOv8cTbHmhDQOj6Kp6/vx3f/PI/7/vUNf7qxf9Op\n7z+V6JbQ8ypnqqmEjXOcO4E1H8DKt8ET6XQb0e1S532BluluR2yM3/y54l8A5IlIBxGJAq4Fpvn5\n/WdzrAkRA9q35v7RXZmxYjt//arI7XBOX0S0M9DMZZNg4mr4/iw49zbYsxHeuxse7QpTLoEvJkHZ\nerejNaZBfrXjF5GxwCTAC7ygqr8TkQkAqvqciKQDBUACUAccBLqr6v4THdvQ71mrntCjqvzw5YV8\nvraUf004j77tWrkd0tlThZ2rYdV7zt3AkTeHU7s67wp0vdR5kNwc7nBMs2cvcJkmaV95NeOemoMq\nvHvXUFrHh9iAK3u3OFVBq951WgRpLSRkOlVBXcdBzlDwRrodpQlRlvhNk7V0y17+889z6ZOVyN9/\nMCj4/fm4pXy306X06vecYSdrDkNMIuSNgq5jnUHoo/0cu8AYP1jiN03au0u3ctc/FnNVv0we/U6f\n5vGw92yhhbbJAAAQB0lEQVRUlcOGz5xCYO0MKC8Db5TTPPRIZ3It27gdpWnmrFtm06Rd1ieDol2H\neHTWWjqkxHNXoPvvb2qi4pyr/K5jnRfCtnzte3P4PXjvx84D4swBzvYu46xbaRN0dsVvXKGq3PP6\nUt5aXMJT1/Xjsj4ZbofU+FShdJVTCKz5ALYuctbX71a63SDw2vWZaZhV9ZhmobKmlhunfM3S4n38\n44fnMqB9ktshuWv/VufN4frdSscmOYPOdxkDnS6GmAS3ozRNlCV+02zsPlTFlc9+ycGKGt6+Ywjt\nWlt/+YCvW+mPnWcCa2fC4d3OS2M5Q48NN5nU3u0oTRNiid80K4WlB7nq2S9JT4zhjdvOIyHGmjx+\nS20NFM8/djdQts5Zn9bdGXy+8xhn6EkbaSysWeI3zc5Xhbu4+YX59GnXipdvHUh8tNVrn9SuQt+d\nwAzYPBfqaiC2tVMl1HmUM9hMbAi8IGdOiyV+0yxNX7aNO/+xmHNyknjxloHERtkVbIMO74X1HzvV\nQes+hMN7wBMB7c6FziOd9waslVBYsMRvmq13lpTw438uYWhuCs/fnE9MpCV/v9XVQvGCY4XAjuXO\n+lbZzt1A3kjIOd/GHQ5RlvhNs/Z6wRbue+MbLu6axp9uHEBURJPpPbx52VcM62Y5hcCGz6C6HLzR\nkDPEKQRyL3HGI7C7gZBgid80e3+ft4mfv72c0T3Sefr6fkR4LfmfleoK2PwVrPsICmfBrrXO+qQO\nTvcRuRc7dwPRLdyN05wxS/wmJLzwxUYefG8l/9Eng8e/2xevx65MA2ZPkXM3UPiRM9ZA9SGnuWj2\nuU4hkDsC2vS0u4FmxBK/CRnPfb6eh6ev5rI+Gfzff/YO3U7d3FRTCZvnOQ+JCz8+9mwgPg06DYeO\nw525DTbTpFniNyHlz5+v5/fTV3Nep2Seu2mAtfMPtgPbnQJgw6ew/lMo3+WsT+vuKwQugvaDISre\n3TjNtwQ88YvIaOAJnMFUpqjqw8dtF9/2sUA5cIuqLvJtKwIOALVAjT+BWeI3x3trUTH3vfENuWkt\neOnWgbRJsEHPG0VdHexY5hQAGz6FTXOhttKpFmo30OlhtOOFkNnfxhpwWUATv4h4gbXAJUAxznCK\n16nqynr7jAXuwkn8g4AnVHWQb1sRkK+qu/w9AUv85kRmr93JbX9fSKu4KF669Rxy06w/+0ZXVe68\nNLbxc9jwOWxbCihEtXRaC+Wc74xF3KYneOyBfGMKdOIfDPxaVUf5ln8KoKq/r7fPn4HPVPUfvuU1\nwIWqus0Svwmk5SX7uOXF+dTUKX/5Xj4D2rd2O6TwVr4bNs72FQSfwe4NzvrYJGg/xCkEcs6HtG72\noDjIAt0ffyawpd5yMc5VfUP7ZALbAAU+EpFa4M+qOtmfwIw5kZ6Zibx12xC+9+J8rn/+a564ti+j\ne7Z1O6zwFdcaelzhTAD7SqBojtNSaONsZ8wBgLhkpyDIGerM07rbHYGLGqNDlKGqWiIiacAsEVmt\nqrOP30lExgPjAbKzsxshLNNcZSfH8caEwXz/pQIm/H0R/31BR+4d2cXa+jcFiZnQ51pnAqfZ6MY5\nsOlLKPoSVk1z1scmQfZ50P4850Fxeh8bd6AR+fOXLgHa1VvO8q3zax9VPTIvFZGpwEDg3xK/705g\nMjhVPX7Gb8JUcotoXht/Lg++t5I/f76BxZv28tT1/eyhb1OTlONM/W9ylvdudgqATV848zXvO+sj\n450eRtufB9mDnc/Waiho/Knjj8B5uHsxTjJfAFyvqivq7TMOuJNjD3efVNWBIhIPeFT1gO/zLOBB\nVZ1xqt+0On5zOt5eXMJP31pGfLSXJ67tx5DcFLdDMv7av815WLx5rtNiaMdyQEG80La3MwLZkSkx\n0+1om7RgNOccC0zCac75gqr+TkQmAKjqc77mnE8Do3Gac/6XqhaISEdgqu9rIoBXVfV3Df2eJX5z\nutbtOMBtryxi/c6DTBzRmTuG5+KxN32bn8N7Yct82DLPmRcXQM1hZ1tiO6cJadY5zpTeCyKi3Y23\nCbEXuExYOlRZw8+mLuPtJVs5Py+Fh6/uTWarWLfDMmejthq2L3MGqN/ytVMY7PfVNHujoW0fX0Ew\nADLznZ5Iw7T1kCV+E7ZUlX/M38JD763EI3Df6K7cdG57u/oPJftKoKTAuRsoLoCti4/dFcQlQ+YA\nZ8ro77xYFh8eVX+W+E3Y27K7nP+duow563YxoH0Sf7i6l73wFapqq51nAyWLYOsiZ166CqclOZCY\nDRl9oG1fyOgLbftBfLKrIQeDJX5jcK7+py4u4cH3VlJeWcsdw3O57cJO1r9/OKg84LxVXLIQti6B\nbUuOvVwGTmHQtrdTVZTe2/ncsm2zriayxG9MPbsOVvLguyuZtnQrndu04H/HduOCzqlIM/5Hbs7A\n4b1OYbBtiVM9tO0b2L3+2Pa4FOeB8ZGpTQ9I6dxs+iCyxG/MCXy8age/fncFW3Yf5tyOrbl/dFf6\nZSe5HZZxU+UB2LHCKQS2L3XmO1dDbZWz3RMJqV0hvadTEKR1c946boJ3B5b4jTmJqpo6Xv16E099\nUkjZoSpG90jnJ6O6kJtmI08Zn9pqKCuE7cudZwc7ljuFw4Ftx/aJSXQKgCMFQWoXp4CIT3WtQLDE\nb0wDDlbWMGXOBp6fvYGKmjr+c0AWt1+YS3ayDURuTuJQGexc5Tw4Ll15bF6x79g+sUmQ0uVYQZDS\nGVJynXcQPMEdRMgSvzF+2nWwkqc/KeSVrzdRW6dc0r0Ntw7pwMAOre0ZgGmYqnMnsHM17Fzrm69x\n5od3H9vPGw3JuU4hkJwHKXnQupMz2H1cYHqYtcRvzGnavq+Cl+cW8er8zewtr6ZnZgLfH9qBcb0y\nrBWQOTOHdjmD2u9aB2XrnPmudU7HdVp7bL+YVk4B0LqTc4cw7CdnVF1kid+YM3S4qpa3Fhfzwhcb\nWb/zEGkto/nuOe24vG+mPQcwgVFT5ST/3euhbH29+QYQD/z4mzP6Wkv8xpylujpl9rqdvPhlEXPW\n7aROoVdmIlf0y+SyPm1Ja2m9gJogqK054+6pLfEbE0A79lfw7tKtvL2khOUl+/EIDMlN4dLebbmw\nS5p1BW2aBEv8xgRJYekB3lniFAJbdjv9w3Rvm8CFXVIZ3jWNfu1a2YAwxhWW+I0JMlVl1bYDfLa2\nlM/W7GThpj3U1ikJMREMzUshv31rBrRPontGApFWEJhGYInfmEa273A1Xxbu4rM1pXxZWEbJXudu\nICbSQ++sVgxon8SAbKcgaJsYY01FTcBZ4jfGZdv2HWbRpr0s3LSHhZt2s2LrfmrqnH9rCTERdE1P\noGvblnRNT6BLeks6pcbTKi7K5ahNcxaMEbhGA0/gjMA1RVUfPm67+LaPxRmB6xZVXeTPsSdiid+E\nmsNVtSzfuo/V2/azevsBVm8/wJrtBzhYWXN0n4SYCNonx5OdHEdOchztW8eTmRRLm4QY0hNjaBFt\ng5GbkzudxN/g/0ki4gWeAS4BioEFIjJNVVfW220MkOebBgF/Agb5eawxIS82yss5Oa05J+fYW5qq\nSvGew6zZfoCiskNsKiunqOwQy0v2MWP5dmrrvn1R1jI6gjaJMaQnxJDWMprW8VEkxUfRut6UFBdJ\nQkwkLWMiiYn0WJWSOSF/LiEGAoWqugFARF4DLgfqJ+/LgZfVuX2YJyKtRKQtkOPHscaEJRGhXes4\n2rX+9/6Bqmvr2Lr3MCV7D7NjfwXb91X65hVs31/Bxl2H2FNeRXlV7Qm+2RHpFVrGRNIyJoKWMRHE\nRUUQH+UlLiqC2Cgvcb7PMZEeYiK9REd4iI7wEhPpzKMiPER6hSivh8gID5FeZznS6yHCI0R4PHi9\n4vsseD2CxyN4xfns9X220c+aHn8Sfyawpd5yMc5VfUP7ZPp5rDHmOJFeD+2T42mfHH/K/Sqqa9l9\nqIrdh6rYU+7MD1TUsL+i2pkfduYHKqo5VFXLroNVlFeVU15VS3lVLYeraqmqrQv6+XgEPCLO5Dn2\nWXB6JxARPL65s06c9Ue2U3/5WEHiHOv7jBxd5ywf2efEBY+cdOHfFhsUqDur1nFRvD5hcEC+61Sa\nTKWhiIwHxgNkZ2e7HI0xzUNMpJeMVrFknMWg8rV1SlVNHRXVtVTW1FFZU0tFdR1VNXVU1dZR/a3J\n2bdOlZpapbZOqa6ro7bOWa5TZ12tKrW1zryuTlHf79SpU8Xl7AeKor51CtSps+zUch3ZVm8/nGXf\n1qOjKx6pFDvyzPLY8onPuf7q459znnZzlwC2j2kZ0zgp2Z9fKQHa1VvO8q3zZ59IP44FQFUnA5PB\nebjrR1zGmADweoTYKC+xUcHtNtg0Hf68WbIAyBORDiISBVwLTDtun2nAzeI4F9inqtv8PNYYY0wj\navCKX1VrROROYCZOk8wXVHWFiEzwbX8O+ACnKWchTnPO/zrVsUE5E2OMMX6xF7iMMSYEnE47futE\nxBhjwowlfmOMCTOW+I0xJsxY4jfGmDBjid8YY8JMk2zVIyI7gU1neHgKsCuA4TQHds6hL9zOF+yc\nT1d7VU31Z8cmmfjPhogU+NukKVTYOYe+cDtfsHMOJqvqMcaYMGOJ3xhjwkwoJv7JbgfgAjvn0Bdu\n5wt2zkETcnX8xhhjTi0Ur/iNMcacQsgkfhEZLSJrRKRQRB5wO57GICIviEipiCx3O5bGICLtRORT\nEVkpIitE5EduxxRsIhIjIvNFZKnvnH/jdkyNRUS8IrJYRN5zO5bGICJFIrJMRJaISFB7qQyJqh7f\noO5rqTeoO3BdqA/qLiLDgIM44x33dDueYPON49xWVReJSEtgIXBFKP93FmdMv3hVPSgikcAXwI9U\ndZ7LoQWdiEwE8oEEVb3U7XiCTUSKgHxVDfq7C6FyxX90QHhVrQKODOoe0lR1NrDb7Tgai6puU9VF\nvs8HgFU44zqHLHUc9C1G+qbmf7XWABHJAsYBU9yOJRSFSuI/2WDvJkSJSA7QD/ja3UiCz1flsQQo\nBWapasifMzAJuA8I/kjwTYcCH4nIQt8Y5EETKonfhBERaQG8CfxYVfe7HU+wqWqtqvbFGbN6oIiE\ndLWeiFwKlKrqQrdjaWRDff+dxwB3+KpygyJUEr8/A8KbEOCr534TeEVV33I7nsakqnuBT4HRbscS\nZEOA//DVeb8GXCQif3c3pOBT1RLfvBSYilOFHRShkvhtUPcw4HvQ+Rdglao+5nY8jUFEUkWkle9z\nLE4DhtXuRhVcqvpTVc1S1Rycf8ufqOqNLocVVCIS72uwgIjEAyOBoLXWC4nEr6o1wJFB3VcBr4fD\noO4i8g9gLtBFRIpF5PtuxxRkQ4CbcK4Al/imsW4HFWRtgU9F5BucC5xZqhoWzRvDTBvgCxFZCswH\n3lfVGcH6sZBozmmMMcZ/IXHFb4wxxn+W+I0xJsxY4jfGmDBjid8YY8KMJX5jjAkzlviNMSbMWOI3\nxpgwY4nfGGPCzP8H4qv4xNFXc34AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f229bb55240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.linspace(0,5)\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, gaussian(x, 1), label='bw=1');\n",
    "ax.plot(x, gaussian(x, 2.5), label='bw=2.5')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our implementation, we choose the bandwidth to be 2.5. (One easy way to choose bandwidth is to find which bandwidth covers one third of the data, which you can try implementing as an exercise.)\n",
    "\n",
    "We'll also need to be able to calculate the distance between points - here's the function we'll use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def distance(x, X): return sqrt(((x-X)**2).sum(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try it out. (More on how this function works shortly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.41421356,  0.        ,  3.60555128])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = distance(array([2,3]), array([[1,2],[2,3],[-1,1]])); d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can feed the distances into our gaussian function to see what weights we would get in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.13598247,  0.15957691,  0.05640321])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gaussian(d, 2.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can put these steps together to define a single iteration of the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def meanshift_iter(X):\n",
    "    # Loop through every point\n",
    "    for i, x in enumerate(X):\n",
    "        # Find distance from point x to every other point in X\n",
    "        dist = distance(x, X)\n",
    "        # Use gaussian to turn into array of weights\n",
    "        weight = gaussian(dist, 2.5)\n",
    "        # Weighted sum (see next section for details)\n",
    "        X[i] = (weight[:,None]*X).sum(0) / weight.sum()\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X=meanshift_iter(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results show that, as we hoped, all the points have moved closer to their \"true\" cluster centers (even although the algorithm doesn't know where the centers actually are)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJySBRJbgECFhMVY2FVE0KmiroGi5YCHa\nel1aq1RJva0KtZVbizWJW6u2iha1F1vRet16f1eQgiuaqFehGkUsKqsiQiKQsATNnnx/f5yZYRIS\nskwm27yfj8c8ZuZ8z5zzPS1+zjef8znfY845RESk+4vp6A6IiEj7UMAXEYkSCvgiIlFCAV9EJEoo\n4IuIRAkFfBGRKNEmAd/MHjWznWa2NmTZ4Wb2qplt9L/3b4t9iYhI67TVCP8xYEq9Zb8GXnPOjQBe\n838XEZEOYm1145WZpQHLnHNj/N/XAxOdc4VmlgLkOedGtcnORESkxWIjuO2BzrlC/+evgIFN/WDA\ngAEuLS0tgl0SEel+3n///SLnXHJT60Uy4Ac555yZNfinhJllApkAw4YNIz8/vz26JCLSbZjZF81Z\nL5JVOjv8qRz87zsbWsk5t9A5l+6cS09ObvIEJSIirRTJgL8UuML/+Qrg+QjuS0REmtBWZZlPAyuB\nUWa2zcyuAn4PnGtmG4HJ/u8iItJB2iSH75y7tJGmc9pi+yIiEj7daSsiEiUU8EVEooQCvohIlFDA\nb0SJq+a5yh2UuOqO7oqISJtQwG/EiqpiHqssYEVVcUd3RUSkTbTLnbZd0eQ4X513EZGuTiP8RvS1\nWC6MH0hfa/icuL22nFvLNrO9trydeyYi0joa4bfSXyu2k19TQlW5Y1xsHybH+Ro9OTSkxFWzoqq4\nxb8TEWktRZpWuqrnYKiAodaTxyoLALgwvskJQYMC1wha+jsRkdZSwG+lwTG9uCXhaEpcNf1i4lqc\n69c1AhFpbwr4YQrk+tvrdyIiraWLti0QWptf4qr5U9kX/PDrj3irardq9kWk09MIvxkCF1jLXS1/\n+2od7w0YzE5XwS68AD+/YivfFO+hfNBoLuuZ0sG9FRFpmEb4TShx1cwv/4LHKgv439v+wNsTLuGf\nGz4JBnuAyo1f8MaEi3ni1rs00heRTksBvwGhqZsVVcXk15Sw/3eLeO6Oe/m6cCcrp2Xy9UbviWLf\nbPyC16dlUl64i5fuvJ+bsm9hWeVOBX4R6XQU8BsQOq3C5DgfCXc9Sd6dfwq2lxfuYuW0THa89Bbv\nTMuktPDA0xs3/G4hz972Bx6rLGBZ5S7NySMinYYCPgdPlDY5zseV8alMjvNRtXsf//fo03XWzyCD\nXoVVvHvRbHoVVpFBRp32Nx99isrivYDTnDwi0ml0m4Afzkg6NCgH0jinxfZjWeUuXu5dwTOvvUT/\n1EGAF+xnM5v7uI800riP+5jN7GDQ75WSzLnL/8IFA4cDcFpsv+DJQ0SkI0W8SsfMtgD7gRqg2jmX\nHon9tPbO1RJXzb7aak6I6c1psf2C21lb8zX5NSUAjBqSyInLHmbltEzyCvOYwQzSSGMRiwDYwhby\nyKNXSjITli+kfMRg1tTsZ5uroJf1UL29iHQK7TXCn+ScOzFSwR7qpmFaYkVVMYurd7Km9mv+Wb0v\nuJ2reg7mhJg+AOy1KnqPOJKxD9zMXvaSQ06dbeSQw172MvaBm+k94kiSieWEHn24JG4gp8X2Uw5f\nRDqFbpPSaWp2y4YERvfHxRzGBbFHBCcyuzB+IINjenFjQhpXxqfyq15pJGwq4KPrbyeJJLLIqrOd\nLLJIIomPrr/dX71jLK8uopf14J/V+5TDF5FOoT0CvgNWmNn7ZpbZDvtrtsDo/uPab+gXE3vQySIQ\n/G1zAW/5Sy8nMpE00tjCFmYyky1sIY00JjIxWL3z+cZNDLGenBbbr9V/eYiItDVzzkV2B2aDnXPb\nzewI4FXgOufcmyHtmUAmwLBhw07+4osvItqfUCWummWVuwDH+fFHNPjXQXFxMWPHjqWgoCC4LIMM\n8shjL3tJIomJTGQJS4LtiSlH8J2Vz5CZcqzy9yIScWb2fnNS5hEf4TvntvvfdwKLgVPrtS90zqU7\n59KTk5Mj3Z06+losl/VM4bKeqQ0G+xJXzRt9qvnx1VfVWb6EJVSmxHPq/9xPeUpcnWAPMO7Kfycz\n5ViN6kU6q7ISWL3Ye48iEQ34ZnaYmfUJfAbOA9ZGcp9tKVCxc8q8n/PrW34bXH546iBueOkZBk75\nDhOWL6RXyoET1cTfXMezt9/T4usJItKO1r0Gqx733qNIpCPSQGCxmQX29ZRz7qUI77PNhM5Zf2HO\nrfS0GB555BH+8fqrLB4aB7Vf03vEkUxYvpAPpl3DnMxryM7O7thOi0jTRp9T9z1KRDyH3xLp6eku\nPz+/o7txSMXFxbzRp5rHKgs4IaY3X7kKdrgqRuyp5o9DT+no7olIU8pKvJH96HMgoW9H96ZNNDeH\nr5xDC/l8Pia7aspdLeD4UWwKz1bu4KrBgzu6ayLSHIF0DsC4Czq2L+1MAb8V+losvSyGxyoL6GU9\nuCXh6I7ukog0V5Smc0ABv1VK/CP8S+IGqhJHpDMpK4G1y6Gqwvse1xPGTKubuknoG3Uj+wAF/BY4\n8OSrGp6p2sGV8QfKOQNtgbt1RaQdBfLy1eWQ/yzFX1fi6x3vtW1fC6lj4PhpFJdW4fP5umUevzm6\nzdQK7eHABG120N2zmgZZpIOUlcDr93t5+aoKsl/6jLF35LKhyD/KL/wY3n+WDQ/+lLHHH0f2vJsO\nrP+v5VFVj6+haAuElmnWH8WHtolIO1r3Gmx9HwafQPb9fyVn6b8AmPTHPHLnnMHIgb3ZsONrJs1/\nmYJ95eTc+XuYOorsn/0QjKi6gKsRfgscaoK21kzeJiJtIO1UGHYy2c+vIefv/xdcXLCvnEnz32b5\nv75i0vy3KdhXHmzLeWE92W/ugeFnwrCTvW00pJvdkauAL9JMpUXw9j3eu3QiG9+k+JOVPLKk7l2z\nGWRQuq8X5z/8T0r39TroyXSP/OUvFOcv9/462PJuw9vuZnfkajgq0kyrF8GKuVD1DcQdBuNmQuKA\nju6VYODrHU/udacFR/KBJ9PNYAY55JBFFmmkAd5cWKlJvci97lR8/frA+CsaL9HsZiWcGuGLNNO4\nmTD5bv9833O9E4AcWnGZ46HVlRSXRfCO/jHTIOU4Rg7sTe6cM0jt14s88oJTly9iUXBK8zzySE1K\nIHf2GYwcfjQcP+1A7r6h1E2ghLObVPIo4Is0U+IAL+gbcGaW9zmUUj4He3ZdFbevquTZdVWR20lC\nX6/sEhiZ6mPhZScc8sl0Cy8dy8iBvSExCcr3w/LbYPVz3Sp10xgFfJEWWL0I3siB+MO8E0BokA+k\nfBZfoaAfcPHoOG4eH8/Fo+Nav5H6F04bupA6wrv4uqH8MDKfWnPIJ9NlPrWGDTu+hkHHwDuPejn8\nos8PndrpJpTDF2mBwKg+8B4I8gCjpsPqR2HTC97yM27smD52Jr4E42fj4sPbSP25bwLfq8u9/JoB\nDja89waTHljlz+FPCaZxQnP4E5nIkn1LmHT/O+Se+ytG9qmAfilw6mUwaFS4h9vpKeCLtEDigAOB\nvLTIu4AbSO+sXgTF62D41IPTPRKG+hdOR58DVeVQuA62rwGg2CUy6f6VFOwtBQg+lCjwZLpf8Is6\nT6Yr2FvGpOn/zke//rZ3R27hJ1ER8JXSEamnoVx8aRHkZcOSucVsfAkWHAMvzvHSO9WlUGbFjJvp\nBf/BpzS9vfbSLhdNw9RkH+tfOE3oC3G9vGDfNwUAn5Uy6/ShdX62hCUk9itn2X+cRmK/8oOeTDfr\n6qvxnXU5pF/c7VM5AQr4IvUE0jShVTjvLoDsnGx+fM9YFnxvA8Xr4JO/e22vP7mB4aljmTU1my/f\n9k4Cob9taHttpamTSbtcNG2NkDz8s+uqeGhVMetf/9/m3+A0+hwv5z71ZjghA5KHk33xd8j60XnB\nVVL79SJ3zhlMOz7Fq95JSgi2ZV19ofewolMvhVMu7TZVOE1RSkeknvp5eoBFedm84a/6WFQ9iSss\nlwFVIynps4FHCiaxnwKeei+H7cBVk7Pr/Lah7bWV0GsIDV0zCFwsDeuiaVsJnbAsJC9/8egMTix4\nm9O3/jes63Gg/VATm4XOeHn6ld7JY9XjZN9wDQw7iUce+S9yb/kBI3vsBJxXsjn7dCY9sJJZF5xL\n9kk13j6iYDqFUAr4IvWE5ukBsrOz+dsbB0r89lPA424S32Mh//g6k/0UBNveIId+G+D7ZDe6vbbU\n1Mmk/kXT4jLHs+uqgieAwGdfgkWmg6FCL76G5ON9CcbpZ59XN9i3dH6bkO1l/wCuS1qDb1gK1CZ7\nk6f1GchI4KObzsI3+gQYODpq0jihIv6IQzObAtwP9AD+4pz7fWPrdoVHHEr3FSitHDUd1i/1gmiZ\nFXP8mLEUfnUgqGeQEbwYmERSnYuBAH1I5W83fkRysq9d78YNDeYNBfDiMsec18t5bWsNN4/3TgK3\nr6rk5vHx4VfSNEdzpyQOd+riwOyZ/gnVGDTaK9v89FX4fBWUfOWlg0JPJl18uuRO8YhDM+sBPAic\nC2wD3jOzpc65TyK5X5HWCKRHtuR5pZUAZ9zo44+X5/LTe7y0TVO37Pe1VO68MJdEfMFpGCZmt0//\nA/l6oMEA/uy6Kl7bWsM5w3rUSfG0W7qnuQ8eCfcBJQl94fSfeIF9+xoYeiL0HwwJ/bxlw04+eHQf\nJY89jHRK51Rgk3PuMwAzewaYASjgS6cTSIuMmg5pE2HYd+CJc8E3aiTzr8pl7jOTyPsmjxnMCN6y\nDwRv2e9DKj92uex+fiSlP/a29dGT4DsG3siGGYtg6Pi6+2zor4rW/kUQCNznpcXy0OrKg0b6ofl8\nX4J16sqdsG15F/ZurxvcA+9ppx48mu9mc+Y0JtJVOoOBL0O+b/MvCzKzTDPLN7P8Xbt2Rbg7InWF\nVrkEcu0DRnnvb94Gn62A9x6EUaNGcuUxCw95y/73WMgARlJbDesXQ9JRUFDoyH68km1fOhZffvD+\nA39VPHqJ4w+vVPLGotYH4UC+fskmb6S/aG1lg+2+BGPjtiLmvF7O7asqmfN6eTD4Fxd3kwf4BKp4\nzp7tfV+92Hsfd4F3Mqg/jUI3mzOnMR1elumcW+icS3fOpScnJ3d0dyTKHKpk8rv3wrcmw4QbofeZ\nG3hic+Yhb9n/B5kUsQGA8j1QWwvrJ1Sx8nuVrJ9QxdHfPXgf42Z6N2qtHFDFPy+sZMPpdcsnW1VH\n7+q915OdnU36SSfw8nvr6N8TXttaw7PrqtiwYQNjx471yhUbsGlPLZcvL2PTntrW960tHWqu+tAA\nXn+K48DJoJuP5hsS6ZTOdiD0bogh/mUincKhqlwGjILLX4UNGzYwadIkdu7xcviN3rLPEh5nEleQ\ny5DDR1LyBYwqiiP1FDiyKI6xtx28j8QBcMHjELcoloRBNUw79sB/ktu2O376XDmrY2tYWVDD/LN7\nNVpNE3rBdubx8STGWTCFs2lPLTnvVJB1ek/++/5bycnx/kKJmT8N5ixn9KiRDCzdzCnfnUzJroJg\ne/3An/NOBa9traGytpyzhsRSWu24N987QbXLRd/6mpt3r5+uCZwMAieMLnqhtlWccxF74Z1QPgOO\nAuKBNcBxja1/8sknO5HOpKioyKWmpjq88bIDXAYZLokkB7gkklwGGXXa+5Dq5lLk7vI5t+znzj16\npnPZOPf45Eb2UVrrfrSs1KU8tN89+EFFcPmN91W4lIf2u5P/9PVBbYHfPfhBRfA95aH97kfLSl1R\naa1zzrmNu2vcvz//jRv32H6X8tB+l5JxU51+Ai6mX4rr/x//4+L6pxzUlpWVVWcfG3fXuB8tK3W3\nvl3mUh7a7+75Z7l78IMKt3F3TXCddlW6z7kPnvPeW+OD55x7aIb33sUB+a4ZMTmiKR3nXDVwLfAy\n8Cnwd+fcx5Hcp0hb8vl8zJo1q86yJSyhhkSuH72MGhIPumX/JGbh8/koK4bi9XiTewGDxjW8j0D1\nzLjqHkw/4kDFzPUXxfGjingem5rQ4IyTgaqcOa+Xc15aLOcM6xFMz4A3In9rey1flULt18XsyH2s\nzu8zyKDvvjL2PHwRh+0pO+iJUH9e+AgL/q8weD1geP8YnpiWwM/G9eTm8fHMPN4r53xlS3XH3M0b\nbt49ClM7Eb/xyjn3AvBCpPcjEilzr81m69uwaIWX6hh0RCqzeufSY91IriCXx5kUvPnq2h9mMX5P\nNqdeB+/+CZKPg5X3eHn6b89tePsXj45jyxsQc3McX5QbQ/w3aQ0ZbNw9x0uVjOHglMnFo+NYWVDD\na1trmJBazfyzewXTOsVljpH9Y9hbUUt1rYPkAayZs5wd86dRu6+wyfLSmH4pHD/vBQ7r5wOqgtcD\n6tf6F5c5SqscN6THNVre2dT9AR0m3PLPLkh32oo0YfUiOHJFNjMnw8ufPMLvL8zlswUjGXw69N02\nkiu25vJ070lc+8tZnHtYNivmemWdP1zuVf8clnzocktfgjH33+J5+VWvPLO5fAlWJ8iH3lX7h/cq\neHhNFd8ZEsNb2xznDIuBgSPwzVlO8fxp5O1rvLw0pl+Kt17voymr9oJ5xvA4HlpdSWmV4973qyit\nciTGGa5sH2VrXmPwyZPxJRw4wNAg39T9AdJ+FPBFmnDgwm42xbuvI/fnPgCOPtfL1pTkjCRr7EeM\nwRcM2IHfNHdahfVLvZu90ibCgBZMw9DofPP+EfkYXwynDOpBWZVjaJ9YnokZQfVlD7Dn4YvIIScY\n7OFAeWn/yx4hduAIviqFh9dU853BMeDg3veruCHde6BJabXj9lWVPD04j7N4hm8sHrgwuK3QIN+p\n5vOJcgr4Ik0IDdrvLvDx+QqvXPO0a71l29+DTS/4ePMdqCqF8+5u+T7qVwsFbshq6Y1YxWWORWsr\ngyPzmWPieXZdFffmV3HDyXFMiNvMM09d32h56S/4Bfueup74OctJPnIExeXw1vZaTknpEbyO4Esw\nNu2pZc3OWoaefC7ffNmDZ92ZzChzwZRN/Zu8NLLvHDq8Dl+kKwlUnA85wwvEgbLKw4d7y3esbt12\nAyeVQHBv7ZTKgeD+8Jrq4Pfz0mK5eXw8u7Zu5Jkb/o3afYVMZGKwvHQmM4MP/J7IRGr3FbJr/jR2\nfLERgGF9YOaY+OBNWwCvbKnmta01vFSYyOOx53Pz+73qXLQNvclLOg+N8EVa4LRrvefZhtbtJw6A\nS5fByzd4N2u1hdZOqXzx6DhKq12wuDKYVhm2nzHTvkvtvkKg6SdC1e4rZO/90zj8NytJTTn4hsiG\n0jRK2XR+EZ8tsyU0W6ZI26lfHZOdnR28qSogpl8K/S57gH1PXR88GQT0nnoTw77/G/ZUwH+cEIsv\nIabzVdoI0PzZMpXSEWmhxh6B2FGPMWxM/bRKdnY2WVkH8vZx/VO4+sGXuDRjKkfNXc7hR6QG2348\n+7dk/Oy37Knwvq8tqu2cT86SFlFKR6SFQqdRvuBxL6XT1JOnOoPiMscRM37DrK9rWfTXv5B0/XJK\nk0ZAFZT2H8G3b3uRV2+aQo/xV1J77k3MP7sXi/5VCQYZw+N4ZUu10jZdnAK+SAuNm3lgzvzVi7wA\nH8nHGLaVQKnkOZNuwjfwKkYOGcDsk+N58bMqKmtjuOPbY9hzzhoeWNebrNN74kswfnVqz+Dvh/dX\npU1Xp4Av0kKBypxA2WRgWWcd2QeEzpc/sv9A1hbV8uJnXkXPzePjGd4/Bvon88TRHdxRiRjl8EVa\nIXGAF+xXL+pceftDCeT0h/ePYW1xLW9tr2X1ztoG5+mR7kkBX6SVWlsr3xmM8Xn/6Y87Ikb18lFE\nKR2RVuoKefvG/Gxcz2CZpUQPBXyRVuoKefvGaLqD6KSUjohIlFDAFxGJEgr4IiJRImIB38yyzWy7\nmX3of02N1L5ERKRpkb5oe59z7g8R3oeIiDSDUjoiIlEi0gH/OjP7yMweNbP+Ed6XiIgcQlgB38xW\nmNnaBl4zgIeBbwEnAoXAHxvZRqaZ5ZtZ/q5du8LpjoiIHEK7PADFzNKAZc65MYdaTw9AERFpuQ5/\nAIqZpYR8vQBYG6l9iYhI0yJZpXO3mZ2I93TNLcBPI7gvERFpQsQCvnPu8khtW0REWk5lmSIiUUIB\nX0QkSijgi4hECQV8EZEooYAvIhIlFPBFRKKEAr6ISJRQwBcRiRIK+CIiUUIBX0QkSijgi4hECQV8\nEZEooYAvIhIlFPBFRKKEAr6ISJRQwBcRiRIK+CIiUUIBX0QkSoQV8M3sIjP72MxqzSy9XttNZrbJ\nzNab2XfD66aIiIQr3GfargUuBP4rdKGZHQtcAhwHpAIrzGykc64mzP2JiEgrhTXCd8596pxb30DT\nDOAZ51yFc+5zYBNwajj7EhGR8EQqhz8Y+DLk+zb/soOYWaaZ5ZtZ/q5duyLUHRERaTKlY2YrgEEN\nNM1zzj0fbgeccwuBhQDp6eku3O2JiEjDmgz4zrnJrdjudmBoyPch/mUiItJBIpXSWQpcYmY9zewo\nYATwboT2JSIizRBuWeYFZrYNmAAsN7OXAZxzHwN/Bz4BXgJ+rgodEZGOFVZZpnNuMbC4kbY7gDvC\n2b6IiLQd3WkrIhIlFPBFRKKEAr6ISJRQwBcRiRIK+CIiHazSlbG5cjWVriyi+1HAFxHpYF9WrePT\nylV8WbUuovsJd7ZMEREJ09C40XXeI0UBX0Skg8VbAkfHj4v4fpTSERGJEgr4IiJRQgFfRCRKKOCL\niERQe5VcNocCvohIG2kouLdXyWVzKOCLiLTAoUbsgeD+YfnrwfahcaM5Jn78QSWXHTHyV8AXEWmB\nQFD/vPJfdQJ2pSuj2lVzuKWws2Yrn1euBeqWXIau3xEjf9Xhi4i0QGCkXu2q+bRyFdWumliLpcKV\n8eFXbzNswNEA7KraCsBR8WPYv7uUvX228mnlKgCOjh/XbjdbhdIIX0SkBQIj9qPix3BM/HjA8Wnl\nKu7K+QPXTbiDDRs2ArCXnWysyuepNX/i2ONH89c7nuGY+PEMjE1jc+VqgAZH/pEU7iMOLzKzj82s\n1szSQ5anmVmZmX3of/05/K6KiHQu1a6KGqpZ/Ls3+NvvFrO7cB83Tvsd2zfuAGD7xh38cuqd7Cws\n4ne33s39t/2Z7VWb6qRy2jO1E25KZy1wIfBfDbRtds6dGOb2RUQ6pS+r1rGx6n2eunMZT//uheDy\n3YX7+M20+Vz7wGUsuP4pdhfuC7b96faF1LpafnHLz4Ij/YGxaUD7pHbCGuE75z51zq1vq86IiHR2\ngeqagbFpDCgZzquPrarTnkEGtYXGrRc9TG2hkUFGnfZnH32O2H192VG9hU8rV7GjegtHx48j3hIi\n3vdI5vCP8qdz3jCz70RwPyIi7SaQgtlRvYWTUr7Ngy/ezuEp/QAv2M9mNvdxH2mkcR/3MZvZwaDv\nS0ni9uXXsaP3ukbLNSOpyYBvZivMbG0DrxmH+FkhMMyf0rkBeMrM+jay/Uwzyzez/F27drXuKERE\nIii0Zj40UH9euZbEo+Gu5b/i8JR+5JHHFraQRhqLWEQaaWxhC3nkcXhKP+5YPpu0EUdybM/Tgxd/\n22NkH9BkwHfOTXbOjWng9fwhflPhnCv2f34f2AyMbGTdhc65dOdcenJycmuPQ0QkYkIvrIbW1RdX\nbwNg0Agf1z5wGXvZSw45dX6bQw572csvH7iawSMGkkBv4q0XX9fu4d2y5Xxdu6fdjiMiKR0zSzaz\nHv7P3wJGAJ9FYl8iIpHWUPrly6p17HZfAV41zoLrnyKJJLLIqvPbLLJIIol7r/8r2zfuoIQiPih7\nlbXlb7GzZiufVLzTbscRblnmBWa2DZgALDezl/1NZwIfmdmHwP8DrnHO7Q6vqyIincfQuNEcboPY\nvnEHv5k2n92F+5jIxGAaZyYzg+mdiUykuHAv86bdz/aNOyiq3U61q2ZAzBCO7Xl6u/XZnHPttrOm\npKenu/z8/I7uhohIHZsrV/Np5SqOiR8fTOdUujJytz3HZaf9vE7pZQYZ5JHHXvaSRBITmcgSlgTb\nB6Qczv0rf01fX2+O6DGME3udHXYe38zed86lN7We7rQVEWlCQymdzyv/RVX/Es6/8uw66y5hCTEp\ncMv//AcxKa5OsAe4ZtY1nDzoLAbEDGFnzdZ2nUtHAV9EpFUMgP/M+hXX3ZwZXOpLSeLO5bM5Zcrx\n3Ll8TrBkE+DmW37DlfN+wFHxx3NSwuR2L8vU5GkiIk0IVOkU1xQEUzBHxY8h1mIZGjeaP9w6jtLa\nEhYvepE7ls/myBFDSbR+DBk1iodfHMjsqbdy4cypTJp7bJ0J1NrjweWhlMMXEWlCpSvjw/LX2Vmz\ntU4ePyCQ44/b05d+vj4U1W6vs15+YS5f9fFSN4dZEmckZrRp/X1zc/ga4YuINCHeEjix19l8WbWu\nwRTM0LjRVLtqapKrwKB/j0F11hs7aDyJFT0pqS1mTK9vt+vNVqGUwxcRaYZD3RkbbwnEWiyfVa/h\ns6o1xFocla68zo1VPWMSOClhMr1j+rd314M0whcRaQPeKL8KMIbGjQ6mgKgAX4/UOrn7jqKALyLS\nBuItgVE9Tw1+P7bn6VCBf96cXkD7Pt2qIQr4IiIR0DumP6cmTAt+78iRfYBy+CIiUUIBX0QkSijg\ni4hECQV8EZEooYAvIhIlFPBFRKKEAr6ISJRQwBcRiRIK+CIiUSLcZ9reY2brzOwjM1tsZkkhbTeZ\n2SYzW29m3w2/qyIiEo5wR/ivAmOcc2OBDcBNAGZ2LHAJcBwwBXjIzHqEuS8REQlDWAHfOfeKc67a\n/3UVMMT/eQbwjHOuwjn3ObAJOLWhbYiISPtoyxz+T4AX/Z8HA1+GtG3zLxMRkQ7S5GyZZrYCGNRA\n0zzn3PP+deYB1cCTLe2AmWUCmQDDhg1r6c9FRKSZmgz4zrnJh2o3syuB84Fz3IEH5G4HhoasNsS/\nrKHtLwQWgvdM26a7LCIirRFulc4UYC4w3TlXGtK0FLjEzHqa2VHACODdcPYlIiLhCfcBKAuAnsCr\nZgawyjnNvAdHAAAKgElEQVR3jXPuYzP7O/AJXqrn5865mjD3JSIiYQgr4Dvnhh+i7Q7gjnC2LyIi\nbUd32oqIRAkFfBGRKKGALyISJRTwRUSihAK+iEiUUMAXEYkSCvgiIlFCAV9EJEoo4IuIRAkFfBGR\nKKGALyISJRTwRUSihAK+iEiUUMAXEYkSCvjdWVER3HOP9y4iUS/cB6BIZ7VqFUyfDrt2ed9vvLFj\n+yMiHU4j/O6goZH8zJlesE9O9j6LSNRTwO8OFi2CuXO999Blw4fD974HCxY0nNZRykckqoSV0jGz\ne4DvAZXAZmCmc26vmaUBnwLr/auucs5dE86+pAFFRV5gnz7d+z59uhfAZ86E8eMhM9M7EQAcdtjB\naZ3AiQKU8hGJAuHm8F8FbnLOVZvZXcBNwH/62zY7504Mc/tyKAsWQE4OLFsGp53mBfg334RvvoHs\n7ANpnXff9d6LimDAgINPFEr5iESFcB9i/krI11XAD8LrjjRbURHk5Xmf33zTezXk44/hjTe8V2Ii\nXHstXHEFvPCC166RvUjUaMsqnZ8Az4Z8P8rMPgT2ATc7595q6EdmlglkAgwbNqwNu9NNFRV5I/vc\n3ANB/swz4fjj4YMPID4epkzxRvhvvw0rVni5/E2bvHUXLPCC/eTJGtmLRJkmA76ZrQAGNdA0zzn3\nvH+deUA18KS/rRAY5pwrNrOTgSVmdpxzrqT+RpxzC4GFAOnp6a51hxFFFi3y0jgBkyfD0097gXzl\nSm/ZL38J77xzoP2227zXlCnw2996y884w0vviEjUaDLgO+cmH6rdzK4EzgfOcc45/28qgAr/5/fN\nbDMwEsgPt8PdXiC/PnNmwwF55kwvR//ll94I/+qr4dJLKS4qwhdYJ//A/8zFo0fje+stb1T/3nte\nLn/y5AP7gkPvT0S6jXCrdKYAc4GznHOlIcuTgd3OuRoz+xYwAvgsrJ5Gi6YqZwYM8NI106Z5aZrL\nLye7qopHgFy8syqVlQBsACY9/DCzLrmE7ORkL9iPHu2N7nNyvModUKWOSJQIN4e/AOgJvGpmcKD8\n8kzgVjOrAmqBa5xzu8PcV3QI5NUby6+vX+9deB02DHw+souLCSR4JnEg6G/wfy+oqSHnSS/Tlj16\ntHdCeeklyMqquw/l80W6PfNnYTqF9PR0l5+vrM8hnXuudyEWyD7lFHLee69OcyreBZFMoKDeT7NO\nP53sc8/1Rvd3360RvUg3YWbvO+fSm1pPd9p2NePGAVAMPPLBB3WaMsiglCTOB0pJIoOMOu2PvPMO\nxcXFXrDXiF4k6mjytK5m7lxITMT3xBPkfvaZl7bBC/azmc0MZpBDDllkkUYaAEtYQipeuseXkKCR\nvUiU0gi/qwlctP3+9xmJF8RTgTzy2MIW0khjEYtII40tbCGPPC/Yx8d7F3RXr/auAzQ1h47m2RHp\ndhTwu6KiIu+u2bPOYiRezn4ve8khp85qOeSwl70sBEaecAJMnerl/6+91vtLYcGCxvfR0IRsItKl\nKeB3RYGbr0pL2ZCQQCaQRBJZZNVZLYsskkgi04wNv/oVPP64l78fNcpbITe38RH8zJnK9Yt0Mwr4\nXdH06TB6NBvee49JZWUUABOZGEzjzGRmML0zkYkUOMekWbPYsHmz9/vERO/9zTcbH8EPGODl+nUz\nlki3oYu2XdHSpRSvW8ekuDgKqqoA78IseLn8vezlF/yCiUwMLi8oKWHSeefxUUkJvsmTvWCemKgR\nvEgU0Qi/K5o+Hd/Uqcw65ZQ6i5ewhET2sqxPHxLZGwz2AbMuvxxfII+fnOxd/A2M4HWRVqTbU8Dv\nipYuhRdeIPuMM8gaPjy4OPWww8idMIFpr7xC7oQJpIb8JAvIPvLIA3n8+iN7XaQV6faU0umKAsH6\nm2/I3rQJJk/mkXffJbekhJEXXADjxzNy6VJy776bSU88wazvf5/s9esP/UDzpqZ0EJEuT1MrdGUh\nM2sW796N7/nnD5r1sri4GN+jjx6YIE1TKoh0O82dWkEBv7sKnWYZDtTcX3utKm9EupnmBnyldLqr\n+tMsZ2d3aHdEpOMp4HdXysmLSD2q0umqmiqj1I1TIlKPAn5XpTJKEWkhpXS6KqVsRKSFwhrhm9lt\nZvaRmX1oZq+YWWpI201mtsnM1pvZd8PvqtShlI2ItFC4KZ17nHNjnXMnAsuAWwDM7FjgEuA4YArw\nkJn1CHNfIiIShrACvnOuJOTrYUCgqH8G8IxzrsI59zmwCTg1nH2JiEh4ws7hm9kdwI+BfcAk/+LB\nwKqQ1bb5l4mISAdpcoRvZivMbG0DrxkAzrl5zrmhwJPAtS3tgJllmlm+meXvCsz1IiIiba7JEb5z\nbnIzt/Uk8ALexIzbgaEhbUP8yxra/kK8p/SRnp7eeeZ5EBHpZsKt0hkR8nUGsM7/eSlwiZn1NLOj\ngBHAu+HsS0REwhNuDv/3ZjYKqAW+AK4BcM59bGZ/Bz4BqoGfO+dqwtyXiIiEoVPNlmlmu/BOHB1p\nANAdHvvUXY4DdCydlY6l8zjSOZfc1EqdKuB3BmaW35xpRju77nIcoGPprHQsXY/m0hERiRIK+CIi\nUUIB/2ALO7oDbaS7HAfoWDorHUsXoxy+iEiU0AhfRCRKKODTvaZ5NrN7zGyd/3gWm1lSSFtXO5aL\nzOxjM6s1s/R6bV3qWADMbIq/v5vM7Ncd3Z+WMLNHzWynma0NWXa4mb1qZhv97/07so/NZWZDzSzX\nzD7x//ua7V/eJY+nRZxzUf8C+oZ8vh74s//zscAaoCdwFLAZ6NHR/W3iWM4DYv2f7wLu6sLHcgww\nCsgD0kOWd8Vj6eHv57eAeH//j+3ofrWg/2cCJwFrQ5bdDfza//nXgX9rnf0FpAAn+T/3ATb4/011\nyeNpyUsjfLrXNM/OuVecc9X+r6vw5jGCrnksnzrn1jfQ1OWOBa9/m5xznznnKoFn8I6jS3DOvQns\nrrd4BvC4//PjQEa7dqqVnHOFzrkP/J/3A5/izebbJY+nJRTw/czsDjP7Evgh/ge54P0j+DJkta42\nzfNPgBf9n7v6sYTqisfSFfvclIHOuUL/56+AgR3ZmdYwszRgHPBPusHxNCVqnmlrZiuAQQ00zXPO\nPe+cmwfMM7Ob8KZ5zmrXDrZAU8fiX2ce3jxGT7Zn31qqOccinZ9zzplZlyr5M7PewP8Cc5xzJWYW\nbOuKx9McURPwXYSneW5PTR2LmV0JnA+c4/wJSbrosTSiUx5LE7pin5uyw8xSnHOFZpYC7OzoDjWX\nmcXhBfsnnXPP+Rd32eNpLqV06F7TPJvZFGAuMN05VxrS1OWO5RC64rG8B4wws6PMLB7vmc9LO7hP\n4VoKXOH/fAXQJf4iM28o/1fgU+fcvSFNXfJ4WqSjrxp3hhfemX4t8BHwD2BwSNs8vOqK9cC/dXRf\nm3Esm/ByxR/6X3/uwsdyAV6uuwLYAbzcVY/F3+epeBUhm/FSVh3epxb0/WmgEKjy/39yFeADXgM2\nAiuAwzu6n808lm/jFWZ8FPLfydSuejwteelOWxGRKKGUjohIlFDAFxGJEgr4IiJRQgFfRCRKKOCL\niEQJBXwRkSihgC8iEiUU8EVEosT/Bz+DVm45lMV8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f229bb731d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By repeating this a few times, we can make the clusters more accurate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def meanshift(data):\n",
    "    X = np.copy(data)\n",
    "    # Loop through a few epochs\n",
    "    # A full implementation would automatically stop when clusters are stable\n",
    "    for it in range(5): X = meanshift_iter(X)\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.17 s, sys: 0 ns, total: 1.17 s\n",
      "Wall time: 1.17 s\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. (In the chart below we are offsetting the centroids a bit to the right, otherwise we couldn't be able to see the points since they're now on top of each other)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFOJJREFUeJzt3XFsXed53/HvE7s2oKYDW0aIfW1rdDBzqBMQaU14W5EG\nJGI0TiTETIAMXrqhcQZyLrJNyz9BDAHjZQsDazNMG9bOKblFKNB0RoDUsislbayB7FAUXkptDmHH\n0Y2cyLBNr5FuIKQBBzuOn/3BS+ZeiRIpXV5e8rzfD3DAc9/38J73heSfjp/73nMiM5EkVd/b+j0A\nSdLOMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9Jhbix3wNo9453vCOHhob6PQxJ\n2lNOnz59ITP3b3bcrgr8oaEhFhcX+z0MSdpTIuKlrRxnSUeSClFU4DebzWtql6QqKSbw6/U6IyMj\nNBqNjvZGo8HIyAj1er0/A5OkHVJE4NfrdWZmZlheXmZ8fHw99BuNBuPj4ywvLzMzM2PoS6q0ygf+\nWtivWQv9kydProf9GkNfUpVVOvCbzSZzc3MdbRNMsLK8wqFDh1hZXmGCiY7+ubk5a/qSKqnSgT84\nOMj8/Dy1Wg1YDfvDHOYoRxliiKMc5TCH10O/VqsxPz/P4OBgP4ctST1R6cAHGB4eXg/9BRY4xzmG\nGOIYxxhiiHOcY4GF9bAfHh7u95AlqScqH/iwGvqzs7Nc5CIzzHT0zTDDRS4yOztr2EuqtCICv9Fo\nMDU1xQADTDPd0TfNNAMMMDU1ddmSTUmqksoHfvvSyzHG1ss4D/HQenlnjLHLlmxKUtVEZvZ7DOtG\nR0dzO++l02w2GRkZ6Vh6OcEECyxwkYsMMMAYYxzn+Hp/rVZjaWnJD24l7RkRcTozRzc7rtJX+IOD\ng0xOTna0Hec4+2r7OHHiBPtq+zrCHmByctKwl1RJlQ58WP3i1fT0T+v2a6txDh482LFkE2B6etov\nXkmqrF11e+ReWQvxubm5jqWXa0s2x8fHmZycNOwlVVqla/iXajabG5ZrrtQuSXuBNfwNXCnUDXtJ\nJSgq8CWpZJWp4X/kR/9nff+pt/9SH0ciSbtTzwM/Is4Bfwv8BHhzK3UmSdL226kr/PHMvLBD55Ik\nbaAyJR3LOJK2qtQVezvxoW0CpyLidERM7cD5JOmKSn6+9U4E/vsy873Ah4BPR8T72zsjYioiFiNi\n8fz58zswHEmlKv351j0P/Mx8tfXz+8ATwL2X9M9m5mhmju7fv7/Xw5FUKJ9v3ePAj4ifjYifW9sH\nfg14rpfnlKRL+XzrVb2+wn8n8JcR8U3gG8DJzPyzHp9Tkjr4fOtVRd1LR1LZ1mr1K8sr62G/5hzn\n+AyfYV9t3557vrX30pGkS5T+fGsDX1IxSn++tYEvqQg+39oavqQCVP351tbwJanF51uvMvAlFcHn\nW1fo5mlSL5R6k62qKv351l7hS1dQ8k22qqxer7O0tHTZ0svh4WGWlpYq/edq4EsbKP0mW1VX6vOt\nDXzpEt5kS1Vl4EttvMmWqszAl9p4ky1VmYEvXWJtxUatVmOBhfVvYR7j2Pq3MxdYWA/7qt53ZS/6\n0WMTvPXYBD96bGLzgwtk4EsbKOUmW1cqRe3VEtU+VkNtX78HsksZ+NIGSrjJVhWXna4Ab7V+agOZ\nuWu2e+65J6V+O3PmTNZqtQRygomcZz6PcSyHGMpjHMt55nOCiQSyVqvlmTNn+j3kazY9PZ3AZXNo\nnzuQ09PT/R2otgRYzC1kbN9Dvn0z8NVvFy5c6Ai8tdAfYCCBHGBgPezbA/PChQv9HvqWtYd9+xxO\nnDhx2dwN/b1hq4FvSUdqU/WbbLnstHBb+Vdhpzav8LVbVLnkUULJqjRs8Qrf++FLV1Cv1y+7yRb8\n9PYKe/kmW1V9tmuptno//J4HfkTcD/wn4Abgv2bmv7vSsQa+dpsq3y3z5MmTHDp0aP07BmvWngB1\n4sQJDh482McRaqt2xQNQIuIG4PeBDwF3A/8kIu7u5Tml7VTVm2yVsOxUl+v1h7b3Amcz87uZ+Qbw\nOPBAj88p6Sp8tmu5ev0AlNuAl9tevwL8gx6fU7pu//gfvsFffeINfuWPb+LLz9zU7+Fsu2az2XHH\nz7UVR2vPdv0Mn+l4tuta6O+VZ7vq6vq+LDMipiJiMSIWz58/3+/hqHB/9Yk3eOtnVn9WUdWXnerq\neh34rwJ3tL2+vdW2LjNnM3M0M0f379/f4+FIV/crf3wTb/vx6s+q8tmu5erpKp2IuBFoAB9gNej/\nGvhEZj6/0fGu0pF2TpWXnZZmNy3L/DDwH1ldlvnFzHz0Ssca+NLOqvKy05JsNfB7/aEtmflV4Ku9\nPo+ka1fVZafaWN8/tJUk7QwDX6qg9z32HU79/jTve+w7/R6KdpGel3Qk7bz6W3/EfW/75urTQC55\nYpfKZeBLFVR/2z+Ft1Z/3tfvwWjXMPClCvrL37wLmDHs1cEaviQVwsCXpEIY+JJUCANfkgph4EtS\nIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpED0L/Iio\nR8SrEfFsa/twr84lSdpcr594dTQz/32PzyFJ2gJLOpJUiF4H/r+KiKWI+GJE/PxGB0TEVEQsRsTi\n+fPnezwcSSpXZOb1/3LEKeCWDbqOAM8AF4AEfhu4NTM/dbX3Gx0dzcXFxesejySVKCJOZ+boZsd1\nVcPPzPu2OJg54EQ355IkdaeXq3RubXv5UeC5Xp1LkrS5Xq7S+d2IeC+rJZ1zwL/o4bkkSZvoWeBn\n5j/r1XtLkq6dyzIlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLA\nl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klSIrgI/Ij4eEc9HxFsR\nMXpJ3yMRcTYizkTEB7sbpiSpWzd2+fvPAR8D/qC9MSLuBh4E3g3UgFMRMZyZP+nyfJKk69TVFX5m\nvpCZZzboegB4PDNfz8zvAWeBe7s5lySpO72q4d8GvNz2+pVWmySpTzYt6UTEKeCWDbqOZOaT3Q4g\nIqaAKYADBw50+3aSpCvYNPAz877reN9XgTvaXt/eatvo/WeBWYDR0dG8jnNJkragVyWdp4AHI+Lm\niLgTuAv4Ro/OJUnagm6XZX40Il4B/hFwMiL+HCAznwe+DHwL+DPg067QkaT+6mpZZmY+ATxxhb5H\ngUe7eX9J0vbxm7aSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JPVJs9m8pvZuGfiS1Af1ep2R\nkREajUZHe6PRYGRkhHq9vu3nNPAlaYfV63VmZmZYXl5mfHx8PfQbjQbj4+MsLy8zMzOz7aFv4EvS\nDloL+zVroX/y5Mn1sF+z3aFv4EvSDmk2m8zNzXW0TTDByvIKhw4dYmV5hQkmOvrn5ua2raZv4EvS\nDhkcHGR+fp5arQashv1hDnOUowwxxFGOcpjD66Ffq9WYn59ncHBwW85v4EvSDhoeHl4P/QUWOMc5\nhhjiGMcYYohznGOBhfWwHx4e3rZzG/iStMOGh4eZnZ3lIheZYaajb4YZLnKR2dnZbQ17MPAlacc1\nGg2mpqYYYIBppjv6pplmgAGmpqYuW7LZLQNfknZQ+9LLMcbWyzgP8dB6eWeMscuWbG6HyNw9j5Ed\nHR3NxcXFfg9Dknqi2WwyMjLSsfRyggkWWOAiFxlggDHGOM7x9f5arcbS0tJVP7iNiNOZObrZ+b3C\nl6QdMjg4yOTkZEfbcY6zr7aPEydOsK+2ryPsASYnJ12lI0l7Ub1eZ3r6p3X7tdU4Bw8e7FiyCTA9\nPb2tX7zq6pm2kqRrtxbic3NzHUsv15Zsjo+PMzk5ue23Vuiqhh8RHwfqwC8C92bmYqt9CHgBONM6\n9JnMfHiz97OGL6kkzWZzw3LNldqvZKs1/G6v8J8DPgb8wQZ9L2bme7t8f0mqrCuF+nbV7C/VVeBn\n5gsAEbE9o5Ek9UwvP7S9MyKejYi/iIhf7eF5JElbsOkVfkScAm7ZoOtIZj55hV97DTiQmc2IuAc4\nHhHvzswfbvD+U8AUwIEDB7Y+cknapU786LH1/UNv/80+jqTTpoGfmfdd65tm5uvA66390xHxIjAM\nXPaJbGbOArOw+qHttZ5LkrQ1PSnpRMT+iLihtf8u4C7gu704lyRpa7r60DYiPgr8Z2A/cDIins3M\nDwLvB34rIn4MvAU8nJk/6Hq0krQH7KYyTrtuV+k8ATyxQftXgK90896SpO3lrRUkqRAGviQVwsCX\npEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkq\nhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCtFV4EfE5yPi2xGxFBFPRMRAW98jEXE2Is5ExAe7\nH6okqRvdXuE/DbwnM0eABvAIQETcDTwIvBu4H/gvEXFDl+eSJHWhq8DPzK9n5putl88At7f2HwAe\nz8zXM/N7wFng3m7OJUnqznbW8D8FfK21fxvwclvfK622y0TEVEQsRsTi+fPnt3E4kqR2N252QESc\nAm7ZoOtIZj7ZOuYI8CbwpWsdQGbOArMAo6Ojea2/L0namk0DPzPvu1p/RHwSOAR8IDPXAvtV4I62\nw25vtUmS+qTbVTr3A58FPpKZK21dTwEPRsTNEXEncBfwjW7OJUnqzqZX+Jv4PeBm4OmIAHgmMx/O\nzOcj4svAt1gt9Xw6M3/S5bkkSV3oKvAz8+9dpe9R4NFu3l+StH38pq0kFcLAl6RCGPiSVAgDX5IK\nYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBv8s0m81rapek\nrTLwd5F6vc7IyAiNRqOjvdFoMDIyQr1e78/AJFWCgb9L1Ot1ZmZmWF5eZnx8fD30G40G4+PjLC8v\nMzMzY+hLum4G/i6wFvZr1kL/5MmT62G/xtCXdL0M/D5rNpvMzc11tE0wwcryCocOHWJleYUJJjr6\n5+bmrOlLumYGfp8NDg4yPz9PrVYDVsP+MIc5ylGGGOIoRznM4fXQr9VqzM/PMzg42M9hS9qDugr8\niPh8RHw7IpYi4omIGGi1D0XE/4uIZ1vbF7ZnuNU0PDy8HvoLLHCOcwwxxDGOMcQQ5zjHAgvrYT88\nPNzvIUvag7q9wn8aeE9mjgAN4JG2vhcz872t7eEuz1N5w8PDzM7OcpGLzDDT0TfDDBe5yOzsrGEv\n6bp1FfiZ+fXMfLP18hng9u6HVKZGo8HU1BQDDDDNdEffNNMMMMDU1NRlSzYlaau2s4b/KeBrba/v\nbJVz/iIifnUbz1M57UsvxxhbL+M8xEPr5Z0xxi5bsilJ1yIy8+oHRJwCbtmg60hmPtk65ggwCnws\nMzMibgbenpnNiLgHOA68OzN/uMH7TwFTAAcOHLjnpZde6mpCe02z2WRkZKRj6eUEEyywwEUuMsAA\nY4xxnOPr/bVajaWlJT+4lQRARJzOzNHNjtv0Cj8z78vM92ywrYX9J4FDwK9n61+PzHw9M5ut/dPA\ni8CGxefMnM3M0cwc3b9//5YnWBWDg4NMTk52tB3nOPtq+zhx4gT7avs6wh5gcnLSsJd0zbpdpXM/\n8FngI5m50ta+PyJuaO2/C7gL+G4356qyer3O9PRP6/Zrq3EOHjzYsWQTYHp62i9eSbouN3b5+78H\n3Aw8HREAz7RW5Lwf+K2I+DHwFvBwZv6gy3NV2lqIz83NdSy9XFuyOT4+zuTkpGEv6bptWsPfSaOj\no7m4uNjvYfRVs9ncsFxzpXZJ2rYavnbWlULdsJfUrW5LOroWq2WvVbvo/6wklcErfEkqhIEvSYWw\npLOTLONI6iOv8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKsSuunlaRJwHXgLeAVzo83B6\nrepzrPr8wDlWQVXm93czc9MHiuyqwF8TEYtbufPbXlb1OVZ9fuAcq6Dq87uUJR1JKoSBL0mF2K2B\nP9vvAeyAqs+x6vMD51gFVZ9fh11Zw5ckbb/deoUvSdpmuyrwI+K3I2IpIp6NiK9HRK2t75GIOBsR\nZyLig/0c5/WKiM9HxLdbc3wiIgba+vb8/AAi4uMR8XxEvBURo5f0VWWO97fmcDYiPtfv8WyHiPhi\nRHw/Ip5ra/uFiHg6Ir7T+vnz/RxjtyLijoiYj4hvtf6OHm61V2qeV5WZu2YD/k7b/r8GvtDavxv4\nJnAzcCfwInBDv8d7HfP7NeDG1v7vAL9Tpfm15vKLwN8HFoDRtvZKzBG4oTX2dwE3teZ0d7/HtQ3z\nej/wy8BzbW2/C3yutf+5tb+ve3UDbgV+ubX/c0Cj9feyUvO82rarrvAz84dtL38WWPuA4QHg8cx8\nPTO/B5wF7t3p8XUrM7+emW+2Xj4D3N7ar8T8ADLzhcw8s0FXVeZ4L3A2M7+bmW8Aj7M6tz0tM/8n\n8INLmh8A/rC1/4fAxI4Oaptl5muZ+b9b+38LvADcRsXmeTW7KvABIuLRiHgZ+HXg37aabwNebjvs\nlVbbXvYp4Gut/SrO71JVmWNV5rEV78zM11r7/xd4Zz8Hs50iYgj4JeB/UeF5XmrHH3EYEaeAWzbo\nOpKZT2bmEeBIRDwC/EtgekcH2KXN5tc65gjwJvClnRzbdtnKHFUtmZkRUYklfRHxduArwL/JzB9G\nxHpflea5kR0P/My8b4uHfgn4KquB/ypwR1vf7a22XWez+UXEJ4FDwAeyVTRkD80PrunPsN2emuNV\nVGUeW/E3EXFrZr4WEbcC3+/3gLoVET/Dath/KTP/pNVcuXleya4q6UTEXW0vHwC+3dp/CngwIm6O\niDuBu4Bv7PT4uhUR9wOfBT6SmSttXZWY3yaqMse/Bu6KiDsj4ibgQVbnVkVPAb/R2v8NYE//31us\nXsr/N+CFzPwPbV2VmudV9ftT4/aN1X95nwOWgD8FbmvrO8Lq6ogzwIf6PdbrnN9ZVuu/z7a2L1Rp\nfq15fJTVuvbrwN8Af17BOX6Y1RUeL7Jaxur7mLZhTv8deA34cevP758Dg8D/AL4DnAJ+od/j7HKO\n72N1IchS23+DH67aPK+2+U1bSSrErirpSJJ6x8CXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDw\nJakQ/x+A3jxEXbKixgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f229b8918d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Broadcasting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "How did our distance function `sqrt(((x-X)**2).sum(1))` work over a matrix without us writing any loops? The trick is that we used *broadcasting*. The term broadcasting was first used by Numpy, although is now used in other libraries such as [Tensorflow](https://www.tensorflow.org/performance/xla/broadcasting) and Matlab; the rules can vary by library.\n",
    "\n",
    "From the [Numpy Documentation](https://docs.scipy.org/doc/numpy-1.10.0/user/basics.broadcasting.html):\n",
    "\n",
    "    The term broadcasting describes how numpy treats arrays with \n",
    "    different shapes during arithmetic operations. Subject to certain \n",
    "    constraints, the smaller array is “broadcast” across the larger \n",
    "    array so that they have compatible shapes.Broadcasting provides a \n",
    "    means of vectorizing array operations so that looping occurs in C\n",
    "    instead of Python. It does this without making needless copies of \n",
    "    data and usually leads to efficient algorithm implementations.\n",
    "    \n",
    "In addition to the efficiency of broadcasting, it allows developers to write less code, which typically leads to fewer errors.\n",
    "\n",
    "Operators (+,-,\\*,/,>,<,==) are usually element-wise. Here's some examples of element-wise operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([12, 14,  3]), array([False,  True,  True], dtype=bool))"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([10, 6, -4])\n",
    "b = np.array([2, 8, 7])\n",
    "\n",
    "a + b, a < b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now this next example clearly can't be element-wise, since the second parameter is a scalar, not a 1d array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True, False], dtype=bool)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a > 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "So how did this work? The trick was that numpy automatically *broadcast* the scalar `0` so that had the same `shape` as a. We can manually see how numpy broadcasts by using `broadcast_to()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(0, a.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(0, a.shape).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here's another example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11,  7, -3])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "It works with higher-dimensional arrays too, for instance 2d (matrices):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6],\n",
       "       [7, 8, 9]])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.array([[1, 2, 3], [4,5,6], [7,8,9]]); m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2,  4,  6],\n",
       "       [ 8, 10, 12],\n",
       "       [14, 16, 18]])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 2, 2],\n",
       "       [2, 2, 2],\n",
       "       [2, 2, 2]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(2, m.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We can use the same trick to broadcast a vector to a matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 20, 30])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = np.array([10,20,30]); c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 22, 33],\n",
       "       [14, 25, 36],\n",
       "       [17, 28, 39]])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m + c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's see what numpy has done with `c` in this case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 20, 30],\n",
       "       [10, 20, 30],\n",
       "       [10, 20, 30]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(c, m.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Interesting - we see that it has duplicated `c` across rows. What if `c` was a column vector, i.e. a 3x1 array?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10],\n",
       "       [20],\n",
       "       [30]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Indexing an axis with None adds a unit axis in that location\n",
    "cc = c[:,None]; cc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 12, 13],\n",
       "       [24, 25, 26],\n",
       "       [37, 38, 39]])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m + cc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's see what numpy has done with `c` in this case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 10, 10],\n",
       "       [20, 20, 20],\n",
       "       [30, 30, 30]])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(cc, m.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Note that numpy isn't actually replicating the memory of the axes being broadcast - it's simply looping over the same locations multiple times. This is very efficient both for compute and memory.\n",
    "\n",
    "The behaviour of numpy's broadcasting seems quite intuitive, but you'll want to remember the explicit broadcasting rules to use this technique effectively:\n",
    "\n",
    "When operating on two arrays, Numpy/PyTorch compares their shapes element-wise. It starts with the **trailing dimensions**, and works its way forward. Two dimensions are **compatible** when\n",
    "\n",
    "- They are equal, or\n",
    "- One of them is 1.\n",
    "\n",
    "When axes have the same dimension, no broadcasting is required. Any axes of dimension 1 are replicated to match the other array.\n",
    "\n",
    "Arrays do not need to have the same number of dimensions. For example, if you have a 256 x 256 x 3 array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:\n",
    "\n",
    "    Image  (3d array): 256 x 256 x 3\n",
    "    Scale  (1d array):             3\n",
    "    Result (3d array): 256 x 256 x 3\n",
    "    \n",
    "Numpy will insert additional unit axes as required to make the array few fewer dimensions math. So in this case the Scale array would be first reshaped automatically to 1x1x3, and then broadcast to 256 x 256 x 3. The [numpy documentation](https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html#general-broadcasting-rules) includes several examples of what dimensions can and can not be broadcast together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We can now see how our `distance()` function works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 1],\n",
       "       [0, 0],\n",
       "       [9, 4]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=array([2,3])\n",
    "b=array([[1,2],[2,3],[-1,1]])\n",
    "c=(a-b)**2; c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  2],\n",
       "       [ 2,  3],\n",
       "       [-1,  1]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.13598247,  0.15957691,  0.05640321])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w=gaussian(sqrt(c.sum(1)), 2.5); w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "...and we can also now pull apart our weighted average:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3,), (3, 2))"
      ]
     },
     "execution_count": 183,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.shape, b.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.13598247],\n",
       "       [ 0.15957691],\n",
       "       [ 0.05640321]])"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w[:,None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.13598247,  0.27196495],\n",
       "       [ 0.31915382,  0.47873074],\n",
       "       [-0.05640321,  0.05640321]])"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w[:,None]*b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.13288486,  2.29313827])"
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(w[:,None]*b).sum(0) / w.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GPU-accelerated mean shift in pytorch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at using [PyTorch](http://pytorch.org/), a Python framework for dynamic neural networks with GPU acceleration, which was released by Facebook's AI team.\n",
    "\n",
    "PyTorch has two overlapping, yet distinct, purposes.  As described in the [PyTorch documentation](http://pytorch.org/tutorials/beginner/blitz/tensor_tutorial.html):\n",
    "\n",
    "<img src=\"images/what_is_pytorch.png\" alt=\"pytorch\" style=\"width: 80%\"/>\n",
    "\n",
    "The neural network functionality of PyTorch is built on top of the Numpy-like functionality for fast matrix computations on a GPU. Although the neural network purpose receives much more attention, both are very useful. Today we'll use PyTorch to accelerate our meanshift algorithm by running it on the GPU.\n",
    "\n",
    "If you want to learn more PyTorch, you can try this [introductory tutorial](http://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html) or this [tutorial to learn by examples](http://pytorch.org/tutorials/beginner/pytorch_with_examples.html).\n",
    "\n",
    "One advantage of pytorch is that it's very similar to numpy. For instance, in fact, our definitions of `gaussian` and `distance` and `meanshift_iter` are identical. So we'll simply import PyTorch's alternate implementations of the two numpy functions we use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from torch import exp, sqrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then we'll use the exact same code as before, but first convert our numpy array to a GPU PyTorch tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def meanshift(data):\n",
    "    X = torch.from_numpy(np.copy(data)).cuda()\n",
    "    for it in range(5): X = meanshift_iter(X)\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try it out..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.59 s, sys: 8 ms, total: 1.6 s\n",
      "Wall time: 1.63 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSJJREFUeJzt3VGMY9d93/Hvv7LjZKqHsaiFIloyaKBiAMVgbGDgumgf\nhomNyNlFNAmQQnko1K3BQQ0XGAcBArkGSvLBgIEAnZekbknUGz0YMQwkHgu7aV1ZmIFRoLE9apyt\nZFuUkoxheWRpQ2OQCgMolXP6MCR3uLuzO7Mkh0Oe7wcYDHnOJe85XOmHO3+ee2+klJAk5eMfTXsA\nkqTTZfBLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMvOOaQ/gsPvvvz+VSqVpD0OS\nZsrzzz//tymlc8fd/kwFf6lUYnt7e9rDkKSZEhE/OMn2lnokKTMGv6Q76na7J2rX2WbwS7qtRqNB\npVKh0+kMtXc6HSqVCo1GYzoD010bOfgj4mcj4lsR8ZcR8WJENHvt90XEsxHxcu/3u0cfrqTT1Gg0\naDab7O7uUq1WB+Hf6XSoVqvs7u7SbDYN/xkzjiP+t4BfTin9EvAB4LGI+DDwFPBcSukR4Lnec0kz\noh/6ff3wv3LlyiD0+wz/2TJy8KcDb/aevrP3k4DHgad77U8DK6PuS9Lp6Ha7tNvtobYVVtjf3efC\nhQvs7+6zcsP/0u1225r/jBhLjT8i7omI7wBvAM+mlL4JPJBSeq23yY+BB8axL0mTVygU2NzcpFgs\nAgehv8Ya66xTosQ666yxNgj/YrHI5uYmhUJhmsPWMY0l+FNKP00pfQB4CPhQRLz/hv7EwV8BN4mI\n1YjYjojta9eujWM4ksagXC4Pwn+LLXbYoUSJS1yiRIkddthiaxD65XJ52kPWMY11VU9KaQ/YBB4D\nXo+IBwF6v9844jWtlNJSSmnp3Lljn3gm6RSUy2VarRZ77NGkOdTXpMkee7RaLUN/xoxjVc+5iFjs\nPf454KPA94FngCd7mz0JfHXUfUk6XZ1Oh9XVVRZZpE59qK9OnUUWWV1dvWmpp862cRzxPwhsRsRV\n4Nsc1PgvA58DPhoRLwMf6T2XNCMOL9lcZnlQ3rnIxUHZZ5nlm5Z66uyLg/L72bC0tJS8Vo80fd1u\nl0qlMrRkc4UVtthijz0WWWSZZTbYGPQXi0WuXr3qF7xTEBHPp5SWjru9Z+5KukmhUKBWqw21bbDB\nQnGBy5cvs1BcGAp9gFqtZujPCINf0i01Gg3q9et1/f7qnfPnzw8t9QSo1+uewDVDztRlmSWdLf0w\nb7fbQ0s2+0s9q9UqtVrN0J8x1vil3H3+0Bm4n9i45SbdbveWZZyj2nW6rPFLGrujwt3Qn00GvyRl\nxhq/lLsjyjuaXx7xS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqM\nwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmayCv5ut3uidkmaR9kEf6PRoFKp0Ol0hto7\nnQ6VSoVGozGdgUnSKcsi+BuNBs1mk93dXarV6iD8O50O1WqV3d1dms2m4S8pC3Mf/P3Q7+uH/5Ur\nVwah32f4S8rBXAd/t9ul3W4Pta2wwv7uPhcuXGB/d58VVob62+22NX9Jc22ug79QKLC5uUmxWAQO\nQn+NNdZZp0SJddZZY20Q/sVikc3NTQqFwjSHLUkTNdfBD1Aulwfhv8UWO+xQosQlLlGixA47bLE1\nCP1yuTztIUvSRM198MNB+LdaLfbYo0lzqK9Jkz32aLVahr6kLIwc/BHxcERsRsR3I+LFiFjrtd8X\nEc9GxMu93+8efbh3p9PpsLq6yiKL1KkP9dWps8giq6urNy31lKR5NI4j/reB300pPQp8GPhkRDwK\nPAU8l1J6BHiu9/zUHV6yuczyoLxzkYuDss8yyzct9ZSkeRUppfG+YcRXgT/o/SynlF6LiAeBrZTS\nL9zutUtLS2l7e3tsY+l2u1QqlaElmyussMUWe+yxyCLLLLPBxqC/WCxy9epVv+CVNDMi4vmU0tJx\ntx9rjT8iSsAHgW8CD6SUXut1/Rh44IjXrEbEdkRsX7t2bZzDoVAoUKvVhto22GChuMDly5dZKC4M\nhT5ArVYz9CXNtbEFf0TcC/wJ8KmU0t8d7ksHf1bc8k+LlFIrpbSUUlo6d+7cuIYz0Gg0qNev1/X7\nq3fOnz8/tNQToF6vewKXpLn3jnG8SUS8k4PQ/2JK6U97za9HxIOHSj1vjGNfd6Mf5u12e2jJZn+p\nZ7VapVarGfqSsjByjT8iAnga+ElK6VOH2n8f6KaUPhcRTwH3pZR+73bvNe4a/4263e4tyzhHtUvS\nLDhpjX8cR/z/HPhXwP+JiO/02v498DngyxHxceAHwL8cw76OVPz8m4PHu5+495bbHBXuhr6knIwc\n/Cml/wnEEd2/Mur7S5LGK4szdyVJ143ly92z4KjyjiRpmEf8kpQZg1+SMmPwS1JmDH5JyozBL0mZ\nMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmD\nX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZGUvw\nR8QXIuKNiHjhUNt9EfFsRLzc+/3ucexLd9btdk/ULikv4zri/yPgsRvangKeSyk9AjzXe64JazQa\nVCoVOp3OUHun06FSqdBoNKYzMElnxliCP6X0DeAnNzQ/Djzde/w0sDKOfelojUaDZrPJ7u4u1Wp1\nEP6dTodqtcru7i7NZtPwlzI3yRr/Ayml13qPfww8MMF9Za8f+n398L9y5cog9PsMfylvp/Llbkop\nAelWfRGxGhHbEbF97dq10xjO3Ol2u7Tb7aG2FVbY393nwoUL7O/us3LDH1ztdtuav5SpSQb/6xHx\nIEDv9xu32iil1EopLaWUls6dOzfB4cyvQqHA5uYmxWIROAj9NdZYZ50SJdZZZ421QfgXi0U2Nzcp\nFArTHLakKZlk8D8DPNl7/CTw1QnuK3vlcnkQ/ltsscMOJUpc4hIlSuywwxZbg9Avl8vTHrKkKRnX\ncs4/Bv4X8AsR8WpEfBz4HPDRiHgZ+EjvuSaoXC7TarXYY48mzaG+Jk322KPVahn6UubioPx+Niwt\nLaXt7e1pD2Nm9Vfv7O/uD8o8fTvs8Dv8DgvFBY/4pTkTEc+nlJaOu71n7s6Jw0s2l1kelHcucnFQ\n9llm+aalnpLy4xH/HOh2u1QqlaElmyussMUWe+yxyCLLLLPBxqC/WCxy9epVv+CV5oBH/BkqFArU\narWhtg02WCgucPnyZRaKC0OhD1Cr1Qx9KVMG/5xoNBrU6/XB8/7qnfPnzw8t9QSo1+uewCVl7B3T\nHoDGpx/m7XZ76Avc/lLParVKrVYz9KXMWeOfQ91u95ZlnKPaJc22k9b4PeKfMb/+5l8MHj9z7wdv\nuc1R4W7oSwJr/JKUHYNfkjJjqWfGHFXekaTj8ohfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbg\nl6TMGPzSHOh2uydqV94MfmnGNRoNKpXKTXdV63Q6VCoVr8aqmxj80gxrNBo0m82bbql5+FaczWbT\n8NcQg1+aUf3Q7+uH/5UrVwah32f46zCDX5pB3W6Xdrs91LbCCvu7+1y4cIH93X1WWBnqb7fb1vwF\nGPzSTCoUCkO31FxhhTXWWGedEiXWWWeNtUH492/F6T0ZBAa/NLP6t9QsFotsscUOO5QocYlLlCix\nww5bbA1Cv38rTsngl2ZYuVym1Wqxxx5NmkN9TZrssUer1TL0NcTgl2ZYp9NhdXWVRRapUx/qq1Nn\nkUVWV1dvWuqpvBn80ow6vGRzmeVBeeciFwdln2WWb1rqKUVKadpjGFhaWkrb29vTHoZ05nW7XSqV\nytCSzRVW2GKLPfZYZJFlltlgY9BfLBa5evWqX/DOoYh4PqW0dNztPeKXZlChUKBWqw21bbDBQnGB\ny5cvs1BcGAp9gFqtZugLMPilmdVoNKjXr9f1+6t3zp8/P7TUE6Ber3sClwa82bo0w/ph3m63h5Zs\n9pd6VqtVarWaoa8h1vilOdDtdm9ZxjmqXfPlzNX4I+KxiHgpIl6JiKcmvT9pXjTj+s+dHBXuhr5u\nZaLBHxH3AH8IfAx4FPjtiHh0kvuUJN3epI/4PwS8klL665TS3wNfAh6f8D4lSbcx6S933wP88NDz\nV4F/OuF9SnOhfna+ftOcmfpyzohYjYjtiNi+du3atIcjSXNv0sH/I+DhQ88f6rUNpJRaKaWllNLS\nuXPnJjwcSdKkg//bwCMR8b6I+BngCeCZCe9TknQbE63xp5Tejoh/B3wNuAf4QkrpxUnuU5J0exM/\nczel9GfAn016P5Kk45n6l7uSpNNl8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmDX5Iy\nY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZMfglKTMG\nvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlJmRgj8ifisi\nXoyIf4iIpRv6Ph0Rr0TESxHxq6MNU5I0Lu8Y8fUvAL8J/JfDjRHxKPAE8ItAEfh6RJRTSj8dcX+S\npBGNdMSfUvpeSumlW3Q9DnwppfRWSulvgFeAD42yL0nSeEyqxv8e4IeHnr/aa7tJRKxGxHZEbF+7\ndm1Cw5Ek9d2x1BMRXwd+/hZdn0kpfXXUAaSUWkALYGlpKY36fpKk27tj8KeUPnIX7/sj4OFDzx/q\ntUmSpmxSpZ5ngCci4l0R8T7gEeBbE9qXJOkERl3O+RsR8Srwz4ArEfE1gJTSi8CXge8C/x34pCt6\nJOlsGGk5Z0rpK8BXjuj7LPDZUd5fkjR+nrkrSVPS7XZP1D4uBr8kTUGj0aBSqdDpdIbaO50OlUqF\nRqMxsX0b/JJ0yhqNBs1mk93dXarV6iD8O50O1WqV3d1dms3mxMLf4JekU9QP/b5++F+5cmUQ+n2T\nCn+DX5JOSbfbpd1uD7WtsML+7j4XLlxgf3efFVaG+tvt9thr/ga/JJ2SQqHA5uYmxWIROAj9NdZY\nZ50SJdZZZ421QfgXi0U2NzcpFApjHYfBL0mnqFwuD8J/iy122KFEiUtcokSJHXbYYmsQ+uVyeexj\nMPgl6ZSVy2VarRZ77NGkOdTXpMkee7RarYmEPhj8knTqOp0Oq6urLLJInfpQX506iyyyurp601LP\ncTH4JekUHV6yuczyoLxzkYuDss8yyzct9RynSOnsXAl5aWkpbW9vT3sYkjQR3W6XSqUytGRzhRW2\n2GKPPRZZZJllNtgY9BeLRa5evXrbL3gj4vmU0tKRG9zAI35JOiWFQoFarTbUtsEGC8UFLl++zEJx\nYSj0AWq1mqt6JGmWNRoN6vXrdf3+6p3z588PLfUEqNfrEzmBa9SbrUuSTqgf5u12e2jJZn+pZ7Va\npVarTeySDdb4JWlKut3uLcs4R7Uf5aQ1fo/4T2hc/1CS5tflNz8/eHzh3k8cud1RmTHpLLHGfwLT\nvIyqJI2LwX9M076MqiSNi6WeYzjqMqqtVovV1dWbLqPaf42kPN2uvHMWeMR/B2flMqqSNC4G/x2c\nlcuoStK4GPzHcBYuoypJ42LwH9O0L6MqSeNi8B/TtC+jKknjYvAfw1m4jKokjYuXbLiDSV1GVZLG\nxcsyj9lZuYyqJI2LwX8MZ+EyqpI0Lp65e0zTvoyqJI2LNf6+iOuPb/OZeHVOSWeNNf4Jm9ZlVCVp\nXEYK/oj4/Yj4fkRcjYivRMTiob5PR8QrEfFSRPzq6EOVJI3DqEf8zwLvTylVgA7waYCIeBR4AvhF\n4DHgP0XEPSPua7JSuv4jSXNspOBPKf2PlNLbvad/DjzUe/w48KWU0lsppb8BXgE+NMq+JEnjMc4a\n/78B/lvv8XuAHx7qe7XXJkmasjsu54yIrwM/f4uuz6SUvtrb5jPA28AXTzqAiFgFVgHe+973nvTl\nkqQTumPwp5Q+crv+iPjXwAXgV9L1taE/Ah4+tNlDvbZbvX8LaMHBcs47D1mSNIpRV/U8Bvwe8Osp\npf1DXc8AT0TEuyLifcAjwLdG2ZckaTxGPXP3D4B3Ac/GwQlQf55S+rcppRcj4svAdzkoAX0ypfTT\nEfclSRqDkYI/pfRPbtP3WeCzo7y/JGn8PHNXkjJzpq7VExHXgB+c4CX3A387oeHMitw/A+ef9/zB\nz+B+4B+nlM4d9wVnKvhPKiK2T3JhonmU+2fg/POeP/gZ3M38LfVIUmYMfknKzKwHf2vaAzgDcv8M\nnL9y/wxOPP+ZrvFLkk5u1o/4JUknNJPBn/sNYCLityLixYj4h4hYuqFv7ucPB5cL6c3xlYh4atrj\nOQ0R8YWIeCMiXjjUdl9EPBsRL/d+v3uaY5ykiHg4IjYj4ru9//7Xeu1ZfAYR8bMR8a2I+Mve/Ju9\n9hPPfyaDn3m6AczdeQH4TeAbhxtzmX9vTn8IfAx4FPjt3tzn3R9x8O962FPAcymlR4Dnes/n1dvA\n76aUHgU+DHyy9++ey2fwFvDLKaVfAj4APBYRH+Yu5j+TwZ/7DWBSSt9LKb10i64s5s/BnF5JKf11\nSunvgS9xMPe5llL6BvCTG5ofB57uPX4aWDnVQZ2ilNJrKaX/3Xv8f4HvcXCfjyw+g3Tgzd7Td/Z+\nEncx/5kM/ht4A5jrcpl/LvM8jgdSSq/1Hv8YeGCagzktEVECPgh8k4w+g4i4JyK+A7wBPJtSuqv5\nj3p1zomZ9A1gzrrjzF86LKWUImLul+lFxL3AnwCfSin9Xe/KwMD8fwa9qxx/oPe95lci4v039B9r\n/mc2+Cd9A5iz7k7zP8LczP8OcpnncbweEQ+mlF6LiAc5OBKcWxHxTg5C/4sppT/tNWf1GQCklPYi\nYpOD73xOPP+ZLPV4A5gj5TL/bwOPRMT7IuJnOPhC+5kpj2langGe7D1+Epjbvwbj4ND+vwLfSyn9\nx0NdWXwGEXGuv4IxIn4O+Cjwfe5m/imlmfvh4EvLHwLf6f3850N9nwH+CngJ+Ni0xzqh+f8GB3Xt\nt4DXga/lNP/ePH+NgxVdf8VB+WvqYzqFOf8x8Brw/3r//h8HChys5HgZ+Dpw37THOcH5/wsOvsy8\neuj//V/L5TMAKsBf9Ob/AvAfeu0nnr9n7kpSZmay1CNJunsGvyRlxuCXpMwY/JKUGYNfkjJj8EtS\nZgx+ScqMwS9Jmfn/nGYAT84PDosAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5004c588d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%time X = meanshift(data).cpu().numpy()\n",
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works, but this implementation actually takes longer. Oh dear! What do you think is causing this?\n",
    "\n",
    "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. To use the GPU effectively, we need to process a *batch* of data at a time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GPU batched algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To process a batch of data, we need batched versions of our functions. Here's a version of `distance()` that works on batches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def distance_b(a,b): return sqrt(((a[None,:] - b[:,None]) ** 2).sum(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 0.4714  0.6740  0.9486\n",
       " 0.5486  0.7749  0.5401\n",
       "[torch.FloatTensor of size 2x3]"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=torch.rand(2,2)\n",
    "b=torch.rand(3,2)\n",
    "distance_b(b, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the two parameters to `distance_b()` have a unit axis added in two different places (`a[None,:]` and `b[:,None]`). This is a handy trick which effectively generalizes the concept of an 'outer product' to any function. In this case, we use it to get the distance from every point in `a` (our batch) to every point in `b` (the whole dataset).\n",
    "\n",
    "Now that we have a suitable distance function, we can make some minor updates to our meanshift function to handle batches of data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def meanshift(data, bs=500):\n",
    "    n = len(data)\n",
    "    X = torch.from_numpy(np.copy(data)).cuda()\n",
    "    for it in range(5):\n",
    "        for i in range(0,n,bs):\n",
    "            s = slice(i,min(n,i+bs))\n",
    "            weight = gaussian(distance_b(X, X[s]), 2.5)\n",
    "            num = (weight[:,:,None] * X).sum(1)\n",
    "            X[s] = num / weight.sum(1)[:,None]\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  9.22296143,   2.36360407],\n",
       "       [  9.3774662 ,  -0.18868986],\n",
       "       [  9.21319103,   0.2641083 ],\n",
       "       ..., \n",
       "       [-16.89390182, -18.92529297],\n",
       "       [-17.19715118, -19.02577019],\n",
       "       [-16.48303604, -19.09796906]], dtype=float32)"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  9.2230   2.3636\n",
       "  9.3775  -0.1887\n",
       "  9.2132   0.2641\n",
       "        ⋮         \n",
       "-16.8939 -18.9253\n",
       "-17.1972 -19.0258\n",
       "-16.4830 -19.0980\n",
       "[torch.cuda.FloatTensor of size 1500x2 (GPU 0)]"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "torch.from_numpy(np.copy(data)).cuda()"
   ]
  },
  {
   "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": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 24 ms, sys: 4 ms, total: 28 ms\n",
      "Wall time: 26.4 ms\n"
     ]
    }
   ],
   "source": [
    "%time X = meanshift(data).cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's more like it! We've gone from 2000ms to 26ms, which is a speedup of over 7000%. Oh, and it even gives the right answer!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSBJREFUeJzt3WFs4/d93/H3d06aVvMDxfTBNWMfGGBmATdgE0DIMmwP\nxDZBnd6hVgt0cB8M3i2gsCADlKJA4SzASD4IEKDA9KRdNhLL1Q+CBgHayMZdt8wxJAQD1iTymmp2\nEtNuqyCOHPvKQOgMAe6c/vZAJE+8O91JJ1IU+Xu/AEHk7/cn/78fz/7gry9///8/UkpIkvLxjyY9\nAEnS6TL4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmYMfknKjMEvSZl516QHcND999+fSqXSpIch\nSVPlhRde+NuU0rmjbn+mgr9UKrG5uTnpYUjSVImIHxxne0s9kpQZg1/SHXW73WO162wz+CXdVqPR\noFKp0Ol0hto7nQ6VSoVGozGZgemunTj4I+JnI+JbEfGXEfFSRDR77fdFxHMR8Urv93tPPlxJp6nR\naNBsNtnZ2aFarQ7Cv9PpUK1W2dnZodlsGv5TZhRH/G8Dv5xS+iXgg8BjEfER4Cng+ZTSI8DzveeS\npkQ/9Pv64X/16tVB6PcZ/tPlxMGf9r3Ve/ru3k8CHgee7rU/DSyddF+STke326Xdbg+1LbHE3s4e\nFy9eZG9nj6Ub/pdut9vW/KfESGr8EXFPRHwHeBN4LqX0TeCBlNLrvU1+DDwwin1JGr9CocD6+jrF\nYhHYD/0VVlhllRIlVlllhZVB+BeLRdbX1ykUCpMcto5oJMGfUvppSumDwEPAhyPiAzf0J/b/CrhJ\nRCxHxGZEbF67dm0Uw5E0AuVyeRD+G2ywzTYlSlzmMiVKbLPNBhuD0C+Xy5Meso5opKt6Ukq7wDrw\nGPBGRDwI0Pv95iGvaaWUFlJKC+fOHfnEM0mnoFwu02q12GWXJs2hviZNdtml1WoZ+lNmFKt6zkXE\nfO/xzwEfA74PPAs82dvsSeCZk+5L0unqdDosLy8zzzx16kN9derMM8/y8vJNSz11to3iiP9BYD0i\ntoBvs1/jvwJ8HvhYRLwCfLT3XNKUOLhkc5HFQXnnEpcGZZ9FFm9a6qmzL/bL72fDwsJC8lo90uR1\nu10qlcrQks0llthgg112mWeeRRZZY23QXywW2dra8gveCYiIF1JKC0fd3jN3Jd2kUChQq9WG2tZY\nY644x5UrV5grzg2FPkCtVjP0p4TBL+mWGo0G9fr1un5/9c6FCxeGlnoC1Ot1T+CaImfqssySzpZ+\nmLfb7aElm/2lntVqlVqtZuhPGWv8Us6+cODs20+uHbpZt9u9ZRnnsHadLmv8kkbusHA39KeTwS9J\nmbHGL+XsNuUdzS6P+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZMfglKTMGvyRl\nxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzWQV/t9s9VrskzaJsgr/RaFCpVOh0OkPt\nnU6HSqVCo9GYzMAk6ZRlEfyNRoNms8nOzg7VanUQ/p1Oh2q1ys7ODs1m0/CXlIWZD/5+6Pf1w//q\n1auD0O8z/CXlYKaDv9vt0m63h9qWWGJvZ4+LFy+yt7PHEktD/e1225q/pJk208FfKBRYX1+nWCwC\n+6G/wgqrrFKixCqrrLAyCP9iscj6+jqFQmGSw5aksZrp4Acol8uD8N9gg222KVHiMpcpUWKbbTbY\nGIR+uVye9JAlaaxmPvhhP/xbrRa77NKkOdTXpMkuu7RaLUNfUhZOHPwR8XBErEfEdyPipYhY6bXf\nFxHPRcQrvd/vPflw706n02F5eZl55qlTH+qrU2eeeZaXl29a6ilJs2gUR/zvAL+bUnoU+AjwqYh4\nFHgKeD6l9AjwfO/5qTu4ZHORxUF55xKXBmWfRRZvWuopSbMqUkqjfcOIZ4A/6P0sppRej4gHgY2U\n0i/c7rULCwtpc3NzZGPpdrtUKpWhJZtLLLHBBrvsMs88iyyyxtqgv1gssrW15Re8kqZGRLyQUlo4\n6vYjrfFHRAn4EPBN4IGU0uu9rh8DDxzymuWI2IyIzWvXro1yOBQKBWq12lDbGmvMFee4cuUKc8W5\nodAHqNVqhr6kmTay4I+Ie4E/AT6dUvq7g31p/8+KW/5pkVJqpZQWUkoL586dG9VwBhqNBvX69bp+\nf/XOhQsXhpZ6AtTrdU/gkjTz3jWKN4mId7Mf+l9KKf1pr/mNiHjwQKnnzVHs6270w7zdbg8t2ewv\n9axWq9RqNUNfUhZOXOOPiACeBn6SUvr0gfbfB7oppc9HxFPAfSml37vde426xn+jbrd7yzLOYe2S\nNA2OW+MfxRH/Pwf+FfB/IuI7vbZ/D3we+EpEfAL4AfAvR7CvQxW/8Nbg8c4n773lNoeFu6EvKScn\nDv6U0v8E4pDuXznp+0uSRiuLM3clSdeN5Mvds+Cw8o4kaZhH/JKUGYNfkjJj8EtSZgx+ScqMwS9J\nmTH4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZ\ng1+SMmPwS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmRlJ\n8EfEFyPizYh48UDbfRHxXES80vv93lHsS3fW7XaP1S4pL6M64v8j4LEb2p4Cnk8pPQI833uuMWs0\nGlQqFTqdzlB7p9OhUqnQaDQmMzBJZ8ZIgj+l9A3gJzc0Pw483Xv8NLA0in3pcI1Gg2azyc7ODtVq\ndRD+nU6HarXKzs4OzWbT8JcyN84a/wMppdd7j38MPDDGfWWvH/p9/fC/evXqIPT7DH8pb6fy5W5K\nKQHpVn0RsRwRmxGxee3atdMYzszpdru02+2htiWW2NvZ4+LFi+zt7LF0wx9c7Xbbmr+UqXEG/xsR\n8SBA7/ebt9oopdRKKS2klBbOnTs3xuHMrkKhwPr6OsViEdgP/RVWWGWVEiVWWWWFlUH4F4tF1tfX\nKRQKkxy2pAkZZ/A/CzzZe/wk8MwY95W9crk8CP8NNthmmxIlLnOZEiW22WaDjUHol8vlSQ9Z0oSM\najnnHwP/C/iFiHgtIj4BfB74WES8Any091xjVC6XabVa7LJLk+ZQX5Mmu+zSarUMfSlzsV9+PxsW\nFhbS5ubmpIcxtfqrd/Z29gZlnr5ttvkdfoe54pxH/NKMiYgXUkoLR93eM3dnxMElm4ssDso7l7g0\nKPsssnjTUk9J+fGIfwZ0u10qlcrQks0llthgg112mWeeRRZZY23QXywW2dra8gteaQZ4xJ+hQqFA\nrVYbaltjjbniHFeuXGGuODcU+gC1Ws3QlzJl8M+IRqNBvV4fPO+v3rlw4cLQUk+Aer3uCVxSxt41\n6QFodPph3m63h77A7S/1rFar1Go1Q1/KnDX+GdTtdm9ZxjmsXdJ0O26N3yP+KfPrb/3F4PGz937o\nltscFu6GviSwxi9J2TH4JSkzlnqmzGHlHUk6Ko/4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmYM\nfknKjMEvzYBut3usduXN4JemXKPRoFKp3HRXtU6nQ6VS8WqsuonBL02xRqNBs9m86ZaaB2/F2Ww2\nDX8NMfilKdUP/b5++F+9enUQ+n2Gvw4y+KUp1O12abfbQ21LLLG3s8fFixfZ29ljiaWh/na7bc1f\ngMEvTaVCoTB0S80lllhhhVVWKVFilVVWWBmEf/9WnN6TQWDwS1Orf0vNYrHIBhtss02JEpe5TIkS\n22yzwcYg9Pu34pQMfmmKlctlWq0Wu+zSpDnU16TJLru0Wi1DX0MMfmmKdTodlpeXmWeeOvWhvjp1\n5plneXn5pqWeypvBL02pg0s2F1kclHcucWlQ9llk8aalnlKklCY9hoGFhYW0ubk56WFIZ16326VS\nqQwt2VxiiQ022GWXeeZZZJE11gb9xWKRra0tv+CdQRHxQkpp4ajbe8QvTaFCoUCtVhtqW2ONueIc\nV65cYa44NxT6ALVazdAXYPBLU6vRaFCvX6/r91fvXLhwYWipJ0C9XvcELg14s3VpivXDvN1uDy3Z\n7C/1rFar1Go1Q19DrPFLM6Db7d6yjHNYu2bLmavxR8RjEfFyRLwaEU+Ne3/SrGjG9Z87OSzcDX3d\nyliDPyLuAf4Q+DjwKPDbEfHoOPcpSbq9cR/xfxh4NaX01ymlvwe+DDw+5n1Kkm5j3F/uvg/44YHn\nrwH/9OAGEbEMLAOcP39+zMORpkf97Hz9phkz8eWcKaVWSmkhpbRw7ty5SQ9HkmbeuIP/R8DDB54/\n1GuTJE3IuIP/28AjEfH+iPgZ4Ang2THvU5J0G2Ot8aeU3omIfwd8DbgH+GJK6aVx7lOSdHtjP3M3\npfRnwJ+Nez+SpKOZ+Je7kqTTZfBLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPw\nS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8k\nZcbgl6TMGPySlBmDX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzJwo+CPityLipYj4h4hY\nuKHvMxHxakS8HBG/erJhSpJG5V0nfP2LwG8C/+VgY0Q8CjwB/CJQBL4eEeWU0k9PuD9J0gmd6Ig/\npfS9lNLLt+h6HPhySuntlNLfAK8CHz7JviRJozGuGv/7gB8eeP5ar02SNGF3LPVExNeBn79F12dT\nSs+cdAARsQwsA5w/f/6kbydJuoM7Bn9K6aN38b4/Ah4+8PyhXtut3r8FtAAWFhbSXexLknQM4yr1\nPAs8ERHviYj3A48A3xrTviRJx3DS5Zy/ERGvAf8MuBoRXwNIKb0EfAX4LvDfgU+5okeSzoYTLedM\nKX0V+OohfZ8DPneS95ckjZ5n7kpSZgx+SZqQbrd7rPZRMfglaQIajQaVSoVOpzPU3ul0qFQqNBqN\nse3b4JekU9ZoNGg2m+zs7FCtVgfh3+l0qFar7Ozs0Gw2xxb+Br8knaJ+6Pf1w//q1auD0O8bV/gb\n/JJ0SrrdLu12e6htiSX2dva4ePEiezt7LLE01N9ut0de8zf4JemUFAoF1tfXKRaLwH7or7DCKquU\nKLHKKiusDMK/WCyyvr5OoVAY6TgMfkk6ReVyeRD+G2ywzTYlSlzmMiVKbLPNBhuD0C+XyyMfg8Ev\nSaesXC7TarXYZZcmzaG+Jk122aXVao0l9MHgl6RT1+l0WF5eZp556tSH+urUmWee5eXlm5Z6jorB\nL0mn6OCSzUUWB+WdS1walH0WWbxpqecoRUpn50rICwsLaXNzc9LDkKSx6Ha7VCqVoSWbSyyxwQa7\n7DLPPIssssbaoL9YLLK1tXXbL3gj4oWU0sKhG9zAI35JOiWFQoFarTbUtsYac8U5rly5wlxxbij0\nAWq1mqt6JGmaNRoN6vXrdf3+6p0LFy4MLfUEqNfrYzmB60SXZZYkHV8/zNvt9tCSzf5Sz2q1Sq1W\nG9slG6zxS9KEdLvdW5ZxDms/zHFr/B7xH9Oo/qEkza4rb31h8PjivZ88dLvDMmPcWWKN/xgmeRlV\nSRoVg/+IJn0ZVUkaFUs9R3DYZVRbrRbLy8s3XUa1/xpJebpdeecs8Ij/Ds7KZVQlaVQM/js4K5dR\nlaRRMfiP4CxcRlWSRsXgP6JJX0ZVkkbF4D+iSV9GVZJGxeA/grNwGVVJGhUv2XAH47qMqiSNipdl\nHrGzchlVSRoVg/8IzsJlVCVpVDxz94gmfRlVSRoVa/x9Edcf3+Yz8eqcks4aa/xjNqnLqErSqJwo\n+CPi9yPi+xGxFRFfjYj5A32fiYhXI+LliPjVkw9VkjQKJz3ifw74QEqpAnSAzwBExKPAE8AvAo8B\n/yki7jnhvsYrpes/kjTDThT8KaX/kVJ6p/f0z4GHeo8fB76cUno7pfQ3wKvAh0+yL0nSaIyyxv9v\ngP/We/w+4IcH+l7rtUmSJuyOyzkj4uvAz9+i67MppWd623wWeAf40nEHEBHLwDLA+fPnj/tySdIx\n3TH4U0ofvV1/RPxr4CLwK+n62tAfAQ8f2OyhXtut3r8FtGB/OeedhyxJOomTrup5DPg94NdTSnsH\nup4FnoiI90TE+4FHgG+dZF+SpNE46Zm7fwC8B3gu9k+A+vOU0r9NKb0UEV8Bvst+CehTKaWfnnBf\nkqQROFHwp5T+yW36Pgd87iTvL0kaPc/claTMnKlr9UTENeAHx3jJ/cDfjmk40yL3z8D55z1/8DO4\nH/jHKaVzR33BmQr+44qIzeNcmGgW5f4ZOP+85w9+Bnczf0s9kpQZg1+SMjPtwd+a9ADOgNw/A+ev\n3D+DY89/qmv8kqTjm/YjfknSMU1l8Od+A5iI+K2IeCki/iEiFm7om/n5w/7lQnpzfDUinpr0eE5D\nRHwxIt6MiBcPtN0XEc9FxCu93++d5BjHKSIejoj1iPhu77//lV57Fp9BRPxsRHwrIv6yN/9mr/3Y\n85/K4GeWbgBzd14EfhP4xsHGXObfm9MfAh8HHgV+uzf3WfdH7P+7HvQU8HxK6RHg+d7zWfUO8Lsp\npUeBjwCf6v275/IZvA38ckrpl4APAo9FxEe4i/lPZfDnfgOYlNL3Ukov36Iri/mzP6dXU0p/nVL6\ne+DL7M99pqWUvgH85Ibmx4Gne4+fBpZOdVCnKKX0ekrpf/ce/1/ge+zf5yOLzyDte6v39N29n8Rd\nzH8qg/8G3gDmulzmn8s8j+KBlNLrvcc/Bh6Y5GBOS0SUgA8B3ySjzyAi7omI7wBvAs+llO5q/ie9\nOufYjPsGMGfdUeYvHZRSShEx88v0IuJe4E+AT6eU/q53ZWBg9j+D3lWOP9j7XvOrEfGBG/qPNP8z\nG/zjvgHMWXen+R9iZuZ/B7nM8yjeiIgHU0qvR8SD7B8JzqyIeDf7of+llNKf9pqz+gwAUkq7EbHO\n/nc+x57/VJZ6vAHMoXKZ/7eBRyLi/RHxM+x/of3shMc0Kc8CT/YePwnM7F+DsX9o/1+B76WU/uOB\nriw+g4g411/BGBE/B3wM+D53M/+U0tT9sP+l5Q+B7/R+/vOBvs8CfwW8DHx80mMd0/x/g/269tvA\nG8DXcpp/b56/xv6Krr9iv/w18TGdwpz/GHgd+H+9f/9PAAX2V3K8AnwduG/S4xzj/P8F+19mbh34\nf//XcvkMgArwF735vwj8h177sefvmbuSlJmpLPVIku6ewS9JmTH4JSkzBr8kZcbgl6TMGPySlBmD\nX5IyY/BLUmb+P7NI+QJg5YZcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5005108278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## course.fast.ai"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "If you found this interesting, you might enjoy the 30+ hours of deep learning lessons at [course.fast.ai](http://course.fast.ai). There's also a very active forum of deep learning practitioners and learners at [forums.fast.ai](http://forums.fast.ai). Hope to see you there! :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}