{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math, numpy as np, matplotlib.pyplot as plt, operator, torch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Create data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "n_clusters=6\n",
    "n_samples =250"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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,
    "deletable": true,
    "editable": 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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Below we can see each centroid marked w/ X, and the coloring associated to each respective cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": 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": 19,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXt8lNW19797MpMbEBImAUIgJgoBBURqRMSqIIgWPBB6\n89L2VfRIbaXeevR4a5PYWlt5D7Wt2hZro+fUC/atIhV61FCCFkUMIgrKVSKShEsCIZDb3Pb7x55n\nZjLkSm5DZn0/n/k88+z9XFZC+D3rWXvttZXWGkEQBKH/Y+trAwRBEITeQQRfEAQhShDBFwRBiBJE\n8AVBEKIEEXxBEIQoQQRfEAQhShDBFwRBiBJE8AVBEKIEEXxBEIQowd7XBoSSmpqqs7Ky+toMQRCE\n04pNmzZVaa3T2jsuogQ/KyuL0tLSvjZDEAThtEIp9UVHjpOQjiAIQpQggi8IghAliOALgiBECSL4\ngiAIUUKXBV8pNUoptVYp9alSaptS6g5/+xCl1FtKqV3+bUrXzRUEQRBOle7w8D3Aj7XW5wBTgduU\nUucA9wFrtNZjgDX+fUEQBKGP6LLga60rtdYf+r8fBz4DMoD5wHP+w54D8rp6L0EQBOHU6dYYvlIq\nC5gMvA8M01pX+rsOAMO6816RTnWD5qnNLqobZAlJQRAig24TfKXUQOBvwJ1a69rQPm0Wzm1R+ZRS\ni5RSpUqp0sOHD3eXOZ2iJ8R5+XY3P9/gYvl2d7ddUxAEoSt0y0xbpZQDI/bPa61f8TcfVEqla60r\nlVLpwKGWztVaLwOWAeTm5vaJO2yJM8APJ8d2yzWvGedothUEQehruiz4SikFPAN8prVeGtK1ErgB\n+KV/+1pX79VT9IQ4OxNUtz08BEEQuoPu8PAvBr4HfKKU+sjf9gBG6F9WSt0MfAF8uxvu1SOIOAuC\nEA10WfC11v8CVCvdM7t6fUEQBKF7kJm2YYQO4LY1mCtZOIIgnG5EVHnkSCB0ABdodTDXOq7eo0m0\nK64Z58CZ0NqLjiAIQt8jgu+nukGzfLub2VnmVxI6gHvNOEez/jfLPIHj6t26zQwf6zx5IAjCaUZV\nFRQVwcKFkJra19Z0CyL4fsI9dmg+mPvUZhc/3+DivQova/Z5ASPw1Q2aRIdq9lAIFfeeSPkUBKEX\nKCqCe+813++5p29t6SZE8P1cM85BvVvzQaWXd8p9QFDQQz3/KekxgIvZWfaTBP7/bmxi6SY39W7N\nf0yJC1wXYHaWnac2u8TTF4TThYULm29DOU29/34t+J0JpzgTFIkOxTvlPi7JsFHv0ew+6qPw3aaA\nR3/NOAd3/rPRv9/EpDQbSze5ea/Cy+OXx9PgMdd6e78HFOSNdrBilxsUvPiZm99vcVPv0fzHBXE9\n/JMLgtBlUlNb9+xPU++/Xwt+Z8Mpljde3eBjaambDw54eWe/j5mZMVwzzsHy7W7W7PPijIc1+7xM\nGmpjZmYMa/Z5Wb7dTYL/t1l6UFN6MHg+wCUZJiHq3f0e/q+GhRNjxdMXhNOVefOgpMRsTyP6dVrm\nNeMcPDQ1tsMzaK2YfYI/hu/2aO4+38Hjl8fjTFDMzrIzOllR3egXcA13nB/LzMwYZmfZyRvjYGq6\njYwB5npnDVZckmHjB5McPHJJPDMzY9hwQLN0k7vVGjuS7ikIPUhVFSxZYrad6Q9vX7kSVq8229OI\nfu3hn+oM2oUTY9ly2MeafV4GxPoC7W+WedhdoxmdrJiQasI51nE5KS5W7fGw74Q5dmq6YkiCjWe3\nublslJ2UeMWkoTZyUhQJ/jTOlpBBXkHoQdoLxbTWH96+cCHU1cHhw1BQAIsXB4+L4Lh+vxb8jhIe\n63cmKB6/PD4Qr1++3c0PJ5s3hfcqvLz16SHmjR7OQ1NjmZ1l56IRHtZ96aHsQDW2gU4AHDZF3mgH\nWw75mJ1lZ/l2N0tL3Tw0NbZNIZeia4LQg7Q0EBs6AGsJeV2dabeEO/y81FQYMAAKC83+AP9rfYTH\n9ft1SKejtFTK2BL9u3MdVDf4+Nm7TRR94iLl7V/SuGQal8SVcc04RyAnf3j9Hmp+eRGJa3/BpDRw\n+zQvfuZizT4vhe82MTvL3mp4KTSMY72VSHxfELqZ1jJrnnjCCPV115l9S8hvuMGcY503b545tqAA\nduwwD4V77oHbboM334RLLoHHHjPXby901EeIh0/rXrUzQZFoVywtNQ+C46//ghOrHwUg76rLmVLw\nDz7wnsnj/7uNPb+ai+9YJXv++iiNHvBd8QDgCwzqXjTCExj4Dc8akjCOIPQCTzxhhLyuzoh2OMXF\nRrATE2HWLBOjLyoyfffeC6tWwbp1Zn/9+uDxb75pvgO8+KI55/BhI/it3auPEMHn5Fh/dYOm6BMX\nKJg+ys7U4YptL/2CSr/YAxyorOD1e68i7bu/Zc9fbsd3rDLQV/7qowxsgorrHuCJWQlcNMIbEPuW\nhF3COILQjXQ2R/666+DZZ+GLL4xg799vPPfZs5uHfl5/Pfh98uRg/+HDRvDHjjXXKi6Gyy7r9h+r\nOxDBb4Hl290s3WS8+jVfeNm8t4rDbxQ1OyaPPEqOlXDwyW+RTDLTyWMFKwL9DeufpWz69/nNpqH8\nz9wEoO03CfHsBaGbCB1gXbgwKP6LF5twTfhEqpUrjdiDEXuL0Dj8PfeYkM7ixUbsb77ZPByeeMJ8\nT0sz3vyTT5rjp0yB6dPN99CxgL5Gax0xn/PPP1/3BVX1Pv3kh026qt4X2F/yfqN+eH2jzv7DcZ3+\n1HGdlv+hjk1J14DOI0+vZa0uokhnkaWLKNJrWavzyNOAtg1O15f85iP97dfq9MPrGwLXbel+Ld07\ndF8QhE5y+LDWjz0W3ILZtnbsbbdpnZxsjjvjDK1BV51xhtb33GP6Q6iqqjJfrOuC1nPmmOMOHzbn\nzJql9fbt7d+7JXtPEaBUd0Bj+/WgbUdy2ncf9bFgRX1g0NbK2Mkb42DnUR+NPrDb4FvTclj2t2Ls\nyemUUEIZZWSRRRFFZJFFGWWUUIJtcDoZ/7GKpd+ewGWj7Px+i+eknPvQQeLwAWNZC1cQuog1QzY1\n1Xjz1kAqnDyYWlRkvPKaGjjzTLj6agoSEjj3iy/YuWQJXHSR8eoXL2bntGmce845FBQUmAHa0aMh\nN9fE+p94wtwvLc2EdO6+27wRhN67Naw3kqKito/rDjryVOitT3d7+E9+2KTTnzqun/yw6aQ+y5P+\n9so6nf7UcT3hz8f1riPewDnffb1epz91XE/7ywmd/tRxveT9Rr1kY6NO+cFfNaCzyNJrWRv4ZJGl\nAT3y9r/q9KeO6wWv1umH1zfqh9c36CUbG5t57B3x8C1bWjpGEIQwOuol5+cbrzs/P3jepZeatuxs\nnQ8a/2dETIze4ffid4AeEdKXn5JizrG2l11mvPr8fOPhd8Sz76ztbUAHPfx+HcNvazDU8qR/MMlO\n5QnN7hrNm2WeQBG1ijpNdhKcl6bYW6tp8EDVvp0ce+F2kkkmn/xm18snn7u4i4rnbsd55yo2MIYN\nlVaWjocth3zkT4tjxW43hJRWCH/7sOL5VnVOC6tSpzXrVxCEENqaUFVVZTztzZvNwCqYLBsrtv63\nv8F111FQXExhyGkVXi8zlGKZ1iyy26nweAJ9hUePgt1OwdGjpmHdOuPVr14N+fknD/i2RVs1e7qZ\nfi34bQ2GWsKOgqUz4vnNJlMB0yqi9pI/pOLVZqZt1b6dPHv7VfiOVTKdvEAYp5BC8skniyymM50V\nx1ZQ/fhcnHeuYsZXxpI/LQ5oChRcswqxJTqMbaGZO6Fpm+EPK6ssszUJTBCEENqqbFlUZMI4AC4X\njBtnwi5FRQGhrXa5eDrstDzyKNElXE0NyZ6B5DG9WWLG0x4PPxo+HGdsLPzbv8F3v2s6rrsu+GAJ\nJQIqbPZrwW+PDw6YUsjvlnvZUOkjJ8WNM0ExJT2Gqek23D7zBldWWc3TD11JwxGTemn9o5dQQg01\n3MVdRuz97b5jldT8Zi4//+RjRqckBkT/jvNjyUnxsLXKe9JCKy2lbYYK++OXxwceBoIghNGSlxw6\nYerwYePhT55sxH/OnOAEqeuuw/n226xNTmZGTQ0VGLG/gzuYz/xmTh2Y//8jYmJY6/XiPHDA3GvH\nDnjnHePhT59ubKmqMrF9MOMAEVBhMyoEv7WFSay692gTVtl80MOGA5pLMmxsqPQxOlmxuwacTid1\nU2+EkDz8FazANjidlOuf5tgLt7Pi2Ipm98yefSNjRqZS3aADJZZzUty89YU3ED4anRLb7C2krRCU\npG4KQicJFdjHHjPbqiozsGq9CVh580BOTQ1rzzyTGZ9/TgklzGd+IDEDCCRmjBg0iLXHj5Nz5pnw\n+efmOi5XsHLmvHnBSVehpRfaegvpJZTWkVOVMTc3V5eWlnb7da14eGgdm9DJVXmjTYmE6gbN77e4\nmZpuw2GDd8ot0Te/o9CZtrbB6TjvXMXUc8eS7dpD0Y+u5HiVeQMYOOd+fvSfPyHBrvjfvV721mpm\nZsbg8mne2e8je7BiwWi7lEgWhJ4kNIQCwe/V1cbjdrng7beDx2dmQkYGq3bv5urDh5uJPcBCFlJG\nGa9//evMTU+Hv/8d9u0z5+3bZ2bnWrn5hYVBL37zZtPWUpinm1BKbdJa57Z3XFR4+C15zs4EFViV\nCmB0ilndavMhE965cbydy0bZmZIewx1rGtlbq5lz60MMm+DguaI/MSV/Fdtso7l0ZAyJjrHE37aK\nusfnknjxjcz5/kNsP6LZUGni9dmDFZPSbEzPtBNrc5GToli6yR2I4wuC0AOEhnmWLAl6+yUlwVII\no0fDpZcawXa72bluHYugzcSMRa+8wtrMTHL27Wt+v+JiU39n/Hizv3kzXHyxaV+82HxfvLhPJ2F1\ni+Arpf4MXA0c0lpP8LcNAZYDWUAZ8G2t9dHuuF9n6Ug4xAr7WOGdPTU+fnFpPE9tdrG31oR5Jjhj\nSPjGA/zthtt4es9AfuC0sXCiuW69+xyqL/wAZ2oqaAIzdVPi4LKRtoDAP355PE9tbuKSDFsgji8I\nQgc4lUHP0Dg+mHPnzTPeveXh2+2wYgU733mHGW+/TYXW5DG99cQMVjCjupq1QA6Yh8Xo0TBihInh\nu1zG2y8uNiJvfS8uNqGdPqyk2V2K8yzwBPDfIW33AWu01r9USt3n3//Pbrpft1LdoPnhWw28U248\n+wGx2j/QGnwrqPfoQBG17MED2HvMR6xNBVIrEx2KvKnprNjtpsGj+cEkO1urfLxT7iPBrgKLpCzf\n7ub3W0x6lxXHFwShA5zKoGdr51x8MdTXQ0UFbN9O9eLFzHjnHSr8Dl97iRkVdXXMiInhY68Xp9MJ\nu3fDd74DAwca0Z81y6RnWnXyi4tNWx/G76GbBF9r/bZSKiuseT4w3f/9OaCECBX80AHcIQmKX1wa\nH5ile804R2Axc2uB8yFxmpEjbSyc4OB7qxrISVH8fouHdfs9gSUNZ2bG8Mgl8azY5eaDg2apQ6ti\nppUOKhk3gtAJTmXQM/ycqioTdlm92uzfcw9s24Zz6VJu+d3vKLRq4eDPxhkwgL9kjWdRWRkr6pon\nZtzi9eIcNw62bzdZP1Z5ZZfLCPzs2eZNJLSGTx/X1Om2QVu/4L8eEtKp0Von+78r4Ki13xo9NWjb\nHqEDuAsnmIHUlgZ6Sw94WfiPBqob4aGpsYHc+EtG2rhgeAzv7jdZPtlJsLfWiP6koTaWlroZnaz4\n81UJjE7p19UsBCGysWL5s2YFY+oQCBUVzJ9P4bvvAjACArH6nffcw4znn6eiogKA/GnTKLjiCvOW\nAKakMpjB2vz8Xhf4iBq01VprpVSLTxal1CJgEUBmZmZvmNMiiQ7VLG2zpYHe32xyUd0IzniYnWX3\nx+BNfv1vNrnYcMBk4+RPiwukYk5Ks3FJho13yn2s2OVuNlAsCEIP0Vq8P9zjLyoKpk/W1RkR//RT\nnq6pMTF6f/ZNzr33svbf/50Zl13GLRMmUGDVvS8sNBO5tm83Qp+fH7x+pFTIDKEn3c2DSql0AP/2\nUEsHaa2Xaa1ztda5aWlpPWhO67S24lX4ylP50+ICi5ib+LuN/5mbwMZK4+nPzIzh8cvjA+vX3n2+\ng4UTY7kgPcZcQDIwBaF3aK0gmVVQrajI5Obfe6/x0q08/cJCCs4+m4+BnEGDTNvkyfDYY+QsWsTH\n3/ymEXvrOnPmGLEfPTro7RcW9k4htFOgJz38lcANwC/929d68F5dItybb2miFsDoFBuv5iWeNOM1\n9HwrHLS01M3MTCP0CyfEktjGwuWCIHQzrcX7Q2P41iIl//qXibdbMfiSEpwAx48bQU9MDEygcq5b\nZ/LurWqczz0XvN6SJcbD70iFzD6iu9IyX8QM0KYqpfYD+Rihf1kpdTPwBfDt7rhXTxCettnWkoMt\npXiGt1mLnYfWvpF8e0HoRVorSFZUZMR5zhzIzjZFz957z3wGDDCfdetMqeSRI2HpUnA6jfe+caNZ\n2CRUzC3Rf+KJoIcfoeEcoH+XRz5VuqMcsZQ0FoRepiNlhq1jtm8PljG+9FJT1vjwYdM+evTJi5t0\n5D4dXfCkB0AWQDl1Worf98U1BEHoBB1ZSMTy/FeuDObGz5gRnAG7cqXJqQcTl7cWNwlfNKWl+4Qv\ntmIRvuhKHyJTPQVB6B90Jk/fOsbK0LFmwC5caNrAxPRXrjT7995rtgMGNJ+1G0pbYaQ+rpJpIYIv\nCEL/oCMLiezYYRYqWbo0WMI4tJKlNVGqqMjE7kOPsYQfWr9PS+mgEVAlM0BH4j699YmUGL4gCP2U\nOXOCsfnWsGLx4fH7lmL34W19FMdHYviCIAhhLF1qMnSWLm39GCu/fvXq5nH68OqbljcfGs9vLY4f\nIUhIRxCE6GHsWFi1qu1jrFTL0Fr6oYTG5MPDNb24Pu2pIIIvCIJgEb4sYUv59PPmmZr68+ZFvMCH\nIyEdQRAEC6s+TlvlEVauNOGelSt717ZuQDx8QRAEi9C0zNbi8JGUddNJomJNW0EQhP5MR8sjS0hH\nEAQhShDBFwRBiBJE8AVBEKIEEXxBiCBqtYdXXAep1Z4OtQtCZxDBF4QIothdzbOuCord1R1qF4TO\nIGmZgtDL1GoPxe5qZjmcJKnm/wVnOZwAXGgfzCuug4FjrHZrKwingnj4gtDLtOWtJyk7X48dxvue\nY4Fj2npACEJnkL8eQeginRXkWQ4njdrHMZ+bF5oquTo2rVVPf5bDGXhAAHw9dlj3/wBC1CCCLwhd\nxBLkrd4T3Bl/RodEf7evnlJvLQDxynaSkFuePiDhHKHbkJCOIHSRWQ4nuTFJlHprWwzThGbY1GoP\njzd+Qam3lvFqAJNsg7jQPrjN61viL+EcoavIX5AgdJEkZefO+DMCYZ1wQkMyAKXeWnJjkhhtS+Ql\n9wHe9xzj67HxvWmyEKWI4AtCNxAaggnHxOy9NGoflzlSAm1gwjkSqhF6CwnpCEIPk6TsxKsYXnIf\n4Jmm8oDAF7urudA+OJCJI5OrhJ6mxz18pdRVwG+AGOBPWutf9vQ9BaGvqdUeXncdBjRXxw5llsPJ\nVu+JQJy/UXt5yX2QN9xVVGpX4Dwr9GNl57SU+dNaVpCkbwrt0aN/FUqpGOBJ4ApgP/CBUmql1vrT\nnryvIPQ1xe5qXnIfACBexfD12GFcEzuMA41NjI8ZwCbPcQAqtYvcmKRmYZ0L7YNZ0rCXLb4THPN5\niFM2rAdHkrK3mqYp6ZtCe/S0GzAF2K21/hxAKfUSMB8QwRf6NbMcTg75XGz21jI+ZgAAy10H2a+b\nWO46yM1xGXziPQ4abo7LIEnZA6Gcde6jbPGdAGCDtybwBmA9OELHBGq1J+DNS/qm0B49HcPPAL4M\n2d/vbwuglFqklCpVSpUePny4h80RhN4hSdk5pF1UahfL/XH5lKMNTLIN5Oa4DN50VbPNV8c2XUfx\nwb1AqIeuudYxnEm2gVRqF2NJJF3FBh4coWMCoWmgkr4ptEefD9pqrZdprXO11rlpaWl9bY4gdBs3\nx2WQG5PEzXEZ3PLT+7g393ISPzcDt596jQcfv6ucu3Iv5+78h2jUXhbYh9KEBjTfjUsnNyYJu01R\nqV38pelAYFB3lsPJjbEjxJsXOkVPuwLlwKiQ/ZH+NkE4LTGDsYcAFSiJEDpYelx7eKapnGtih7HJ\nc5zRtgSeKPwFL//8vwB4dPY3yV31RybnjEPt+pI35i6itvIgv374EXLch/nOT+6l1GNm4K7zHKVS\nu5gT4yQhJoahxPKsq4Jj2sPgkIJqoUXWBKEtevov5ANgjFIqGyP01wLX9/A9BaHHMIOxB4FgSQQr\nFHPM5+GfnmqO4eVAYxP7dRM7fvEHdj66LHD+8cpDlM69FefvHubNHz1EY2UwjLnz0WWsws5NP72P\n/b7GQBw/yebg+tgRLGkwoZ+dHhMKatRe4lVMp8s6CNFLj4Z0tNYeYDHwBvAZ8LLWeltP3lMQepJZ\nDicL7EOZZBsYKIlwoX0wuTFJ7PTVcQwvg4nhjvhMxhz1sO/ZV5udn0ceMZVN/L9vfp/4Sjd55DXr\n//jZv/L+oS+5NX4U1zqGc61jGFfHplHsrmaL7wS5MUnk2E0s/zNvfeDerZV1EIRQetwd0FqvBlb3\n9H0Eobf4Uhvv2yqJ8L7nGKXeWhbYh5KgYrg5LoMMWzz5Iycz4s1XKJj9DU5UHiKPPO7gDuYzn0IK\nySefLLIAWMEK4tPTuGLV01QNSeBNdzUL40x+Q632cEx7AgO+g5Sdz731bPEdZ537CHfGn8HrrkOB\nrB0ITup633NMwj1CAPkrEIROUOyuptRbSzqxHPPPjh0fM4B0YnHh4+a4DNa5j2Llzf94wleZsXYt\n02fMoKSyhPnMJ4ssiigCoIwySighMX0oF676IwPGnEEdmp2eOmpjPYFJWq+6DwGwzn2Eq2OHAspv\nkQpk7TzrqiBemZd2K8xjVeSUvHwBRPAFoVNYpRD26yZedR9isLKz1XuCSlys8lRR5m1gm64LHH99\n3AhG5+Qw54mfsfwb36eQwoDYAxRSSA013PHUcwwaO56htlje8hyhChcvNFaw2lvNAvtQxtsGsM1X\nRxPaH945Tm5MElfHmsy2lnLwL7QPZoJnoGTyCAFE8AWhE7zvOcZ+3cQk20DOjjFiOj5mADsa6jiO\n1zje2hx7XHv5Sf0uDu38nNcWP0QyyeST3+x6+eRzF3fxxx/+Bz9/8xUGjs5iIDYOajcl3iMAbPUe\nx2Z59BpmxQbF3QrVhBdvs75LFU4hlD7PwxeE0wkr//3W+FGB8Mk2rxH7SbZB5NgSWeAYygL7UD70\n1LJ+xzb+e873aKw8zHSmk0UWZZSxkIWUUUYWWUxnOo2Vhymc/Q3+57NSTuADoN7/5NilG9ih6wGI\nU0omWAmnjNJa97UNAXJzc3VpaWlfmyEI7fKK6yDPuiq41jEc0DT54+7bdB3XOoYRr2JYVvkp71x0\nLfWVhwLn5ZFHCSXUUEMyyUxnOitYEegflD6URz5Yw8DUFJp8PvZ461E2OFMlkGRztLgcoiAopTZp\nrXPbO048fEE4Bax0yCbt4yX3QT731gdi95s9xxkfM4AL0kYx66bvNDtvBStoSndw6V9/S2O6o5nY\nA8y46Ts4U518I3YY18enE2+LYbuv3uTix6WL2AtdQgRfEDqBVbN+nfsIpd5a4pTixtgRjLQFY+U7\ndD3LXQfZ4juOvu97XPLA4kBffHoaU1ctY/BVX+WiVcsYmD400Pfth36Mvu97gRo51uAsQK1PauQL\nXUfcBUHoBNas2msdwwO1bJKUnaLGYMWQYcRyTewwRnni+dxbz3cffpQfY2dt0QvcvOp/KBuTRjbx\nxI45G1b9kQ1zF3H5wu+wtPDnrHMfAVQgs8aql1+uG/voJxb6EyL4gtAJwtMfrRo6ccpk0QzCxkFc\nPNtUwRGfm0pcTPbWcf1P78GzaB4NqWmg3STZHIBm4JgzuPS95dicybzvOcb1cSOa1eb5acJZPNNU\nzs1xGa2ZJAgdRkI6gtAJQjNkLG//8cYvuMwxhNyYJI7jY6SKY5uvjkpcDCPWzID1eYl1JjNRDSI3\nJonvxg0nQ8UzBDvpzlTmxKQ2e4g866qg2F1Nhi2enyacxSBll+UPhS4jgi8Ip8gshzNQx+Z9zzHu\njD+DG2NH8GDCmYy3mXo3qTYHL7kPUK6bAIizmf9ya91HWe2t4ggernA4uT4+nWJ3NTu8J9jsqWWB\nY2izCVOhDwFBOFUkpCMIp0iSsnNn/BnN1pG1Jjzdn3Amr7sOUau92FEsiE3D4VaBmbiDiAm5kgoI\n+kgVx37dhEPZAhk5tdpDo/ZxrWOYzJoVuoQIviB0gfAZrgDlvkaeaSpnlIpntacKAIfbZgqsOYZy\nxOOmUrsYbxvAxJiBgfIIAONjBgSWQLSw1se9MXZEi2mZsni50FHkr0MQupFa7eGRhs/Zr5tw23yB\niVmXOYYE6trEoXjJfZCJMQO5Pm5E4FzrwfHThIHNrtneWrWyeLnQUUTwBaEbsQqrjVRx3Bo/igx/\nfn7oYOvVsUOJVzEdDs8k+Ve3as2Lt65zoX1wYPUre7WdzUUweSEkpprj6qs4qU2ILkTwBaEbCfXG\nreUPX3cd4jN//XowXnhnPXHLi29pZSsrrGSVewAYVjSM4ntN/8X3mO3mIk5qE6ILEXxB6EbCY/qh\nSyLmxiSd0qCrNWg7yTYwsLJVSw+M0IeNfaFpm7ww2D+5hTYhupDiaYLQg7S06HlnCS3UFq9sMjgr\nnERHi6fJX40g9CBJyt5sYLY9Wsq4meVw4q6DIUVOLr7OLvF34ZSRiVeCEEG0NMEqSdkZ+adhvHu7\nGYgFMwC7fonZdoTOHi/0T8TDF4QIorUUzLHzoKzEbAE2PgHrCsFdB9MLgsfVV5k+DVy4OJiNIwO2\nAojgC0LEUF8F65+BPRNg2oWQFBK62bESdq+GrOmQek9gFUXCR+A2F5kHAUDsgKC4y4CtAF0UfKXU\nt4AC4Gxgita6NKTvfuBmwAvcrrV+oyv3EoT+zuYiWFlZTfltFcS8B3dfEczECRfsiddBxQdmG8rk\nhcbr15jLLxAUAAAgAElEQVS3gfVLROSFIF318LcCXwf+GNqolDoHuBYYD4wAipVSOVprbxfvJwj9\nlskLofYZJ1vXwfUXNg/pJKY2D8WEe/yhk6qsEM/6JSaM46ozD4fdq027dR2ZiBV9dEnwtdafASh/\nLfAQ5gMvaa2bgL1Kqd3AFOC9rtxPEPo7STY7t104rE0Brq8yXvy0e4yYW8JdfC/UHYYDmyF9Mky+\n2Ry3bz3sLYbRc5p7+xLXjz56KoafAWwI2d/vbxMEoRXaE+D6Knj/Cdi/Hj4vhjMugy/WgaceLvaf\nt+dNI+57iyExDRwDgmK/4LnmnrzE9aOPdgVfKVUMDG+h60Gt9WtdNUAptQhYBJCZmdnVywnCaUl9\nlfHWL8tvWYDrq+CFeVDuf0fOngVH95jv+98PhmbGzoPViyF1rPHuJ1wXjOmHXss6XsI70UW7gq+1\nnnUK1y0HRoXsj/S3tXT9ZcAyMDNtT+FegnDaYglt/WF4dwlcmm/a1y8x4r1jJWReAq8thOrtpi/B\nCcMnG889ORtQwVh97AD45ovBtwXHAPMpvtfE8Rc81/xNYvJCs++uC2b3SHin/9JTIZ2VwAtKqaWY\nQdsxwMYeupcgnLZY4pt5mdl31wdz7D9/04RunOOM2CedAQ1UM+Hq4IBuUy3sWwdpl1azf72Tz4tN\nvv6VS02/9baw500zaPv+EyY/3+qz7n9pPsx6TMI7/Z0uzbRVSi1QSu0HLgJWKaXeANBabwNeBj4F\n/he4TTJ0BOFkJi80YmulPTgSg+GXYZONCM8vMjH4nbkF/OrLc/nHkzvx1IM9ARqqoSl9Jw++fy7P\nFBeQmGaE/ZMXzTXqq42op082+1+uN9uL7zGhm8kLzT0uXBxsE/ovUjxNEPqQ+ip49QYj0tbAKhiR\nHjvPCLcCVlcVsORJE3NJUiN4cOpaGt7LoYqd/LdtBrU+Uxb5MvK58bICbA4T8jlzlnlLuOgek72z\nt9g8YGYU9M3PK/QMUjxNECKccLG/cqkJuSjMYOvqxUag11LAOgoD59XqCn69ZwbfGreM/961iFpv\nRaBvHYWk7IPz9haQPQua6ky7Ihj3b6gKTsgSjz66kOJpgtBHbC5q7tnvWAlvF5r4/bLzjTjXU81m\n9XSz8/LIo/FQPb/bfjU2bz155DXrf7f2aYbNqmZwZjCrx55owkUABz8xcfv3nwieI8XVogPx8AWh\ni5xqSuPkhSazRoXs1x2GD54wWTPxKZA50cn/eXstf4mdQY2rgjzyuIM7mM98Cikkn3yyyAJgBSsY\nxAh+nLOWg8VOjgwy1x00MjhQGzvA3GPf28H7wslzACRNs38iHr4gdBFLLDtSuji0LzHVCPC6QnNu\nYioMSANPgzl27ALIngGp5HD/+WsZPnQEJZRQRhlZZFFEEVlkUUYZJZQwJHEETy5ey3lfzQHAbVZU\nxHXCbK3yDF+91wzUTrguaIs1eGtl6YT/TEL/QARfELpIZ8QyvM86d+w8KCkw+fhnXGr6ju+Ds64y\naZkN7+UwP3YZNdRQGBLPByikkBpquDFnGV9bnMOBzTDigmB/U01zWyzh37EyaIvVZnnz4T+T0D+Q\nLB1B6GbaCoe01mcVOgMz27bcX+zMKp9wNHYnf3bNIIZ6fs2vA2EcgDLKuIu7wJHIralria/MIekM\nqP3CPDyyZsAUf0gn9N4Stuk/dDRLRzx8Qehmwr3l8D5rwlPVjmBIZew8k0J50T1GnK9casonVO2E\nKozYH6eC6UwPhHEWsjAQ3pnOdGrcFTxZOYMqdjJwqPHQv/03Uz0zMdVM6Cq+12zbs1Pon8igrSD0\nMlZYx5pJ+/mbMOpi8/3M2UFxtrJ0nsOIPZiBWYASSqihhru4i+lMD7Qfp4LnmEHRhR8zeaGzmQff\n2qIpQvQgHr4g9DJWfHyYf/br58WmhPFl+cFFS1z1pi/nIicXJtzS7PwVrMBLIv8+7HVUXGJA7C1m\nZt/CoARnwKPfXGTeIhRm0pWVsSNEHyL4gtALWNk5VTuCcfOv+mvYnDnLePOOAWZmbfG9QXE+azb8\n+Z8FzM/MD1xrECO4NWUtj66by4aP15KaNCLQd1VKPuftLeC9JcaTtwZeraUPYwdICCeakZCOIPQC\nVhinrKT5ylMzCpoPnlrxdUdisMqlYwB8a0IBNfvgQ57m+8lr+dHqHFLHQsOGHL6fvJbfe2dw6dBb\nOG9vASmj4dzvNB+otRY/b630sgzeRgci+ILQC1hCO3aeWZawtXTHKYuNwFv9rjqTqjlwOCzILmBe\n8o84vtnJ9ldgz/+aB4SjOofFKR9z4RwnBz+BjAvNdRJTm2f/tFb2WFa+ih5E8AWhFwhdkzY1TFTD\nBTdUdGMHmHILFtmznBwHPnsVju42bfZE8Bx18sGTJjz03hKITTTZOe2tamUtl3hpKwuvCP0LEXxB\n6GNaEmUrzDJ2nvHw979vPPfsy+HgFiP2yWdC0kgYNpGA2A+bbAaBO5qJY8X2Zz0m4ZxoQARfEPqY\nUO/fItTrT0wztW+yZ8Abd5kHQIITaj43n5EXmpz9YZPhKzeb8gzhs36h5XCNrGsbXYjgC0IEEhrz\n/+RFk7KpMStfOcdB9kwofdIIvT0xuHD5gLTmwt7W24M1SCtx++hBBF8QIhBLiNcvMTF8K70yNmRA\nN/mM5kKuONlTb+/tQcQ+uhDBF4Q+pL2UyLHzTCrn2Hkni7eVXz95YedWsJIwTvQiE68EoRcJL53c\nXhniHStN3v6OlW2f25kFTKSGTvQiHr4g9CLh4ZT2vO3Q/rbOlTCN0BGkPLIg9CJdmdV6KmWXheig\no+WRRfAFQRBOc6QeviAIgtCMLgm+UmqJUmq7UupjpdSrSqnkkL77lVK7lVI7lFJXdt1UQRBaomoH\nPD/XbAWhLbrq4b8FTNBanwvsBO4HUEqdA1wLjAeuAp5SSsV08V6CILTAG3ebTJ437u5rS4RIp0tZ\nOlrrN0N2NwDf9H+fD7yktW4C9iqldgNTgPe6cj9BEE7myqXNt4LQGt2ZlnkTsNz/PQPzALDY728T\nBKGbSR0L31nV11YIpwPtCr5SqhgY3kLXg1rr1/zHPAh4gOc7a4BSahGwCCAzM7OzpwuCIAgdpF3B\n11rPaqtfKXUjcDUwUwdzPMuBUSGHjfS3tXT9ZcAyMGmZ7ZssCIIgnApdzdK5CrgXmKe1rg/pWglc\nq5SKU0plA2OAjV25lyAIgtA1uhrDfwKIA95SSgFs0FrfqrXeppR6GfgUE+q5TWvt7eK9BEEQhC7Q\n1Syd0W30PQI80pXrC4IgCN2HzLQVBEGIEkTwBUEQogQR/M7SUAubXzVbQRCE0wgR/M6yfQ1seM5s\nBUEQTiNkAZTOMm5m860gCMJpggh+Z0lIgskL+toKQRCETiMhHUEQhChBBF8QBCFKEMEXBEGIEkTw\nBUEQogQRfEEQhChBBL8rhE7CkglZgiBEOJKWeao01MI/fwP7NgXbNjxntpK2KQgRhUs38KV7O6Mc\n44hVCX1tTp8R3R5+uFfempdutR8tD/ZvX2PEPvN8Mwlr3EyYekPbE7LkLUAQ+oQv3dv5zLWBL93b\n+9qUPiW6PXyrTAIYr3zrKihdDp5GuOA60x7qyVdsDXr0oTNuE5KC17BE3Wq3Hg7jZp58P0EQeoVR\njnHNttFKdAt+eJkEd5PZHq+GlfmQmg32uKAnP+0mGDEBsqYERdwSe4twUf9kFWxaDl9+BJcsan4/\nQRB6hViVwFmxk/vajD4nugU/vEyCPc5sK7dC7QEo3wKT8oKhmoQkSFlgPPiWPPWGWvN2cP415qGw\n+VXw+B8i5VugbKN49oIg9BnRLfjhTJwLjnhIPwfW/NqIvj0uGI6xRH/cTCPs7kYj8lboxgr95F4D\n7/7ZfM+YZB4a1nUEQRD6iOgW/ND4ekJSc49/zk+MaI+5NBimqdhqwjplG0FjQjW71kH2VHD4Qz8j\nJ0HlduPRJ2eY7ajzxLMXBKHPie4snfDa9qFZNGUbjYCXbTQPhIxJZn/1z805Chicbt4CtqwwIj8p\nzzwILLGfcbsJB1nhHcnOEQShD+n/Hn64Fx+KFZo5fig4SLtlhenLmmI8+iFnwFv/Ba46qk+4cFIJ\nSekw6itmIPZYJdUuO87yLVBXBTXlRuxryqHyU+PZtxbzFwSh15Bc/GgQ/JZSIUMfAvZ42LbctJ84\nFBxwtWLwtQegppyC17fz9PovWHvnxeRQCRtfgIM72HnwBDMeX88tMydQMAuTzXP+t2HTy+Y6IIum\nCEIEYOXiA1GbsdN/Bd8S9XDRheYPgXEzoeEY7H0fjlWaNkvs/eJdcMe/U7h6BwAzfvsea5/6CTnJ\nKez8cD0zfruBimONFL5SCnEDKXjypmA4yOuB9HEwYa549oLQx0gufhdj+EqpnymlPlZKfaSUelMp\nNSKk736l1G6l1A6l1JVdN7WTWKJupUKGhnOyphgxtx4GjjjIvtB494qg2F9+BwV/eJHCl/8VOLXi\naD0zbnuEVRu3M+Px9VQcrQv0Fb5YQsF/3GYeIEnpJpZfutxk77QWv+/I7FuZoSsIbeLSDexxbcal\nG1rtj/ZwDnR90HaJ1vpcrfV5wOvATwGUUucA1wLjgauAp5RSMV28V9uEi2JbpQ5CB2St2bVbVsCe\nf8HQsUbsJ8ylesWjPP3H3zc7NY886o/Yufqex6k/Fk8eec36n351DdXrX4Zaf6x/2Fhzr62rWi62\n9smq9hdFl4XTBaFN2iud0F5/ew+M/kKXQjpa61CXcwAmRwVgPvCS1roJ2KuU2g1MAd7ryv3aJDxW\n39bas6Ex9U9Wme/KZgZa333GbGsP4KwpZ+09s5jx8CtUHGskjzzu4A7mM59CCsknnyyyAFjBCkYM\njmftfV/DOdAFcYOM6A8aaq6vw2wE8z33mvZr8MgYgCC0SXvhmmH2LA579tPka8ClG07y8qMlvt/l\nGL5S6hHg/wDHgBn+5gxgQ8hh+/1tLZ2/CFgEkJmZeeqGdEYUQx8GYy6Fbf+AxlqIT4JpN8OW12BQ\nGiQkkzMkk7V/Op8Ztz9GSWUJ85lPFlkUUQRAGWWUUGLE/s6LyZk602T9WJOuUrPNx9NkPrnXNLex\npeyhtuwVBOEkWiudYIVyPNpDlW8/Vb79xNnMsVbfMHsWHu1mjCO3xQdGfwoHtRvSUUoVK6W2tvCZ\nD6C1flBrPQp4HljcWQO01su01rla69y0tLTO/wQWlii2J56hNNSaAVpL7L/2IJR/bGLv24vNrNtt\n/yDnrGyW/foxaqihkMJmlyikkBpqWHbDVHKmfQ1Q4Go0Ofnp40yoqGqv2W5ZYbKCwLxZeBpP/ecV\nBKFdgp67ZowjlzGO8wOibvVtbPgHu9ymKGJLgt6fKm226+FrrWd18FrPA6uBfKAcGBXSN9Lf1je0\nlotvlTgOzZuv2mv64gaZ7JqmE+x883kW/XI1ySSTT36zS+eTz13cxaLnNrD2rGxyHEdMR/0RmPOQ\nEfiGY+YhMnJSsEzDJn8qqD1evHdB6AZCPfaDnjJGOcadFOoJFe1RjnFUeys45N1nzveZOH64J9+f\nsnu6mqUzJmR3PmD9NlcC1yql4pRS2cAYYGNX7tVhGmrhgxdh44vBAdzwQc+j5fD3fDi634iwNSN2\n3ExT0TI5A5qOw8evsfPjTcx4dBUVNQ1MZzpZZFFGGQtZSBllZJHFdKZTcayRGfkvsfPgCXOP2koz\nKJw1xTxEJuXBVxcFU0XPv+bk8I4gCKeM5Yl/0PC/fObawIcNb7HXtTUg4Fb/R43/DMTxz4u/nFSb\niTbX6ZqAJx86iGuFi073cA50PYb/S6XUWMAHfAHcCqC13qaUehn4FPAAt2mtvV28V8fYvsZk3YAJ\nyUxecHJ8/90/w/4twXOGjzP1762SClc9AGUbqd7/uUm9PGZCLysws3BLKKGGGu7iLqYzPdBeUdPA\njN++x8f3X4ZzeIYprvavZca7j7HD7reD9fanXNcrvw5B6I+05s1bHvsAlUyVr5wqXzle3MTgwKvd\nDLGlc8i7j72urWTHTuBL93Zy4qZgc21idOz5pHkrGeUY12wQ19rvDzH8rmbpfKONvkeAR7py/VPC\nKpegCQp8+KDntJvA54HGE1C1B/ZvhQm1Jy2A4mx4llsuPiMw6QqC2Th/uf5CFr2whRXHVjS7/S0L\nrsA5EJOds2k5jP8anKgy4aEtr5mDNIIgdAFLkENDMqMc4xhsG8pAlQJAjTeeI/oANZ6DHNEHAIgh\nFgCvdrPX9Qm73JtItWVQ5SvH6R0RGMxt8jWQahvJMHtWuxk8p9Ogbv+baZuQFFytqjVSMuDfCk3o\np2oPHNhm3gwsIQ4R5IKrx0FsIoUrNgMwIjmRtXdcRM6wgawtGMuMgpcDbwD5+fkU3HodrP2tEXpH\nvCmLXFNuHiblW0yO/8S5PfCDC0L/JdyjH2bPAky6pdMzgmH2LEob3uCIr5J4BtBIHclqGKm2DNxH\nFQw21/HiAiBG2TlafQyS4IQ+xhA1PHDNL93b+dxjIgBbG98hKSa12WBvOKdTSmf/E/zOMGGuWeWq\naq+Jq8cPMjNtNSbOf2gXAAXfuwKGZPL0399m7a/vJKf2Qxg2lpzhZ7P20UHMuP8ZbrlmHgUFBbDq\nZ0bgS35n0jDTRgcrZo6Y0LE0TEGIcsK95pY8ektcB8amsMe1mSM+UxqlETP7vUYf5IVHXmfNcx/w\n3BtPkXHWUI7rIyTbhuIui2fB5fO5/IZcrntgDo2coNy9m7FxFzDKMQ6PdnPUeyAQFjo7dmqr3vsw\nexbV3orAAyOSiW7BbzwOX3xgBNoqwRCoc/+2GXgFOLiDgluv4UePPYMz0WHeBhqOwZYV5EzK4+NX\np+Cc6p9xO+0m86BorDWDvxPnBgU+pRuzcdqqAioIpznhXrPlXVsefai3bYVg4hlIIydIZigD7UNY\n+vDvePHR1QDceOUP+ceaVQw78wyq99Rz9cwZHK48wguPvg7AdQ/MAXTgQZMdO5FsJrKj8QMO+/Yz\n0DakxQwegIOeMg559+H0jGBgbEov/HZOneish2+VNVj3lBH7pPRgvF/5j6mtNBOnho83+xqcTufJ\n4wEHdxmxD4h6hsnnT0qHjHPbt+FU6+NIuQWhHzPKMY6zY6cGhN3KlBloSwl49lYWjRWCaeQEqbaR\nTBkwhyd+toz/efTVwPUOVhzmisuv5L9XPs2cmVdzuPJIoO+FR1/nzV99RHbsxMCDZq9rK1+6t3Nc\nH6FeH2Nr4zuB9vZsjWSiU/AtsbRi9dkXBgV7wlyTLnn+NXDFj+Gq/zQhmfC4u8O//q0V/w+l8lPz\nwNj2j9YFeau/hs7WVaf2M7RVK0gQTnPaS4UMTbEcZs8KpFaCprq6mr8/W9zs+Dzy8FZqHv7W7/FW\nclINrOV/fpUPK//FMHsWZ8dOBTSfuTbg0z4AHMpMmHTpBjY2rOKE72iHbY0kojOkY4lk1pTgilYW\nLQ36hq9pC+bBAM2zgSC4kPmkPPNQaE2QWxgg7hRSbkHo57SV/RKagun0jOArCVfwUeM/OeTdR1rS\nKNb8cw0zL5/FoYrD7dbAcqYn8/DrP6QqaTfl7sGMjZuCSzdgVw6afA3UeA6SGjOCDNtZHPZ8SZWv\nnLqGWi5OzDstRD6U6PTwLbFMyehYOQbrjcAqc2zFzyfMNbV43vovePfZYHvpckgYbB4crV174tyW\n3xwEQQBOLmng0g3saNrIjqYPADgv/nLGOM7Hoz0AnBM3jaExmQyzZzFh7Lk898ZTDEkfTAklgUmS\nRRQFJk+WUMKQ9ME8suoORo0xld2Peg8GJltZIZpU20iGO84EICduCokkUadr2Ov6pLd/JV0mOj38\nzjJuplnOcN8meOU/TQjIWgqxYqtJtyzfYkQ+dJJX6CIs1puE9QAQD10Q2qSlsghWzRu7sjPKMY5j\nvsOBrJ2j3kqqfOV4Gt04Y0Yw/KwhLP7t9Tz8rd9TSGGg4CEEa2D96ncPMGbMGOqpxU4cVb797Gj6\ngETbIDzaHUjPtLlsHPLu4+zYqWQ4ctjlLiU44Hf6IILfERKSTCG08i0mNl+1Nxg/z5piVrZKzQ4K\nuiXk1lq2FVvNwwJE5AWhDcLDOKMc49jr+gQvHtBwRsx46vSxwISoQ959DI3JBDRVPlOu64ivkiO+\nSqp31/HE7S+0WQPrv25fxqtvvYw704UbM5/msOdL6qkNFFsDRYZjdLPsIOuBc7ohgt9RQnP2L1lk\nwkFgBH5eYcvnhI4VWDn4giC0Sng6ZqhXDzA0JpMq3/7A5KvDnv0MVClkOMbg1R5qvAcZFOOkfPch\n7pv7MEcqj5FHXiCMExrDn850VlSsYN6sr/PIqjsYPeYsBtgGkxM3hcMe89aQHTsxEKcPTbmM9AlW\nraG0jpx5/rm5ubq0tLSvzRAEoY8I9/BduiHg4cdgJ82eyW7XJs6Jm8ZBT1ng4XB27FQ82m1KJdSO\n5hsX3EhFRUXgunnkBWpgJZPcrAYWgDM9mf/3QRHpqSOapWcOjcnkvPjLASK6fIJSapPWOre946Jz\n0DYcWTNWECKC8BTHWJXA2LgpjI6djF05OOz5kkPefYGCac1r3JuYepJzAN+46epm113BCmzpmp/+\n9QfY0mkm9gCzb5zGicGV7HJvYm3dSwy0DWFoTCaHvPv40r2939TEF8EHmcQkCBFO6EIm1iQn8zC4\ngLFxU4hVCaTZR5GoBlPlLmf2f57HDx+8KXD+kPTBPLrqLi64aiK/WHUHaenOQN9198/h5geu5UzH\nJBzE4aaRT5ve5bz4ywP3Op0mV7WFxPBB1owVhAjHqm8DqsWwiks38FHDGuoJvqXf/pNbSbOP4g9P\n/56HX7+NEWPSSFRJZIyBJavv5+45P2PujZfz7QeuIME2iNGxkxluP5OPGteQooYFSihb9zpd4/ah\niOCDpEgKQgQSHs+3KwefuTZgV/Zma9Ja9epDxT7VNpLs2Ak8kD+BsxemMMAZSwwOpiTM4aCnjJSJ\n6YzYOIJzh0/j06Z3OeKr9JdN0CSqQez37QAfHPUeICVmWLPB29MZEXxBECKSlgqoebQbj/YExD50\nkRKPduPVHmKUI+CZ73FtZoAzNnDNcveugHhfMfI69rg2U6drAqmdVkaQKdWgAguf25VDPHxBEISe\nInziVaiXf8x3iHPipgUeAABj46Y0O9+lG/BoN+m2szjoK8OLyeIJFe+W16tVZMdOAAh4/Z2N3Ufq\noigi+IIgRCRWxk4o4TV0QLHLXYoXN3H+iVqWwFo5/ENjMvHhJdWWQUrM8IB4t5QCaleOZtcYG3dB\ns/t3VMgjdVEUEXxBEE4brIXHrdWvtja+A0Cttyow0zbcew+toW8J+x7XZjzawy53KdXeisA12xPp\njgp5y28OfY8IviAIpxWW57/HtZkqXzmptgwGqGSw0WzVqdA3hNBZsntdW9nlLuVM+6RmufYdEemO\nCnlLbyeRgAi+IAinJZboWp46mNWn2l91ylQXiFF2zou7vFmIpj2RjlQh7ygi+IIgnJZY4uvSDRgR\nVyctfdhSvD07dmIgVt8WkTrw2hW6ZaatUurHSimtlEoNabtfKbVbKbVDKXVld9xHEAQhHKv8wti4\nC5oJc2vlEELLN7RVMqG/lFMIpcsevlJqFDAb2BfSdg5wLTAeGAEUK6VytNbert5PEAShI7QUb2+p\n/HL4MaHne7QHj3YHFkU53ekOD//XwL00X6xvPvCS1rpJa70X2A1MaelkQRCEnqCltWbDvfa21qM1\nef92drk39Rsvv0sevlJqPlCutd6iVLPVXzKADSH7+/1tgiAIfUZn0yUjNb3yVGlX8JVSxcDwFroe\nBB7AhHNOGaXUImARQGZmZlcuJQiC0CadzbKxwj79ZfC2XcHXWs9qqV0pNRHIBizvfiTwoVJqClAO\njAo5fKS/raXrLwOWgVkApTPGC4Ig9DSROmv2VDjlkI7W+hNgqLWvlCoDcrXWVUqplcALSqmlmEHb\nMcDGLtoqCILQ6/SnsE6P5OFrrbcppV4GPgU8wG2SoSMIwunI6T7ZKpRuE3ytdVbY/iPAI911fUEQ\nBKFryBKHgiAIUYIIviAIQpQggi8IQsRjlTQ2dXOEU0UEXxCEiKc/1rXpC6RapiAIEU9/So3sS0Tw\nBUGIePpTamRfIiEdQRCEKEEEXxAEIUoQwRcEQYgSRPAFQRCiBBF8QRCEKEEEXxAEIUoQwRcEQYgS\nlNaRs+aIUuow8EUXLpEKVHWTOd1NJNsGkW1fJNsGYl9XiGTbILLtC7XtDK11WnsnRJTgdxWlVKnW\nOrev7WiJSLYNItu+SLYNxL6uEMm2QWTbdyq2SUhHEAQhShDBFwRBiBL6m+Av62sD2iCSbYPIti+S\nbQOxrytEsm0Q2fZ12rZ+FcMXBEEQWqe/efiCIAhCK/QbwVdK/VgppZVSqSFt9yuldiuldiilruwj\nu36mlPpYKfWRUupNpdSISLFPKbVEKbXdb9+rSqnkSLHNb8O3lFLblFI+pVRuWF8k2HeV//67lVL3\n9YUNYfb8WSl1SCm1NaRtiFLqLaXULv82pQ/tG6WUWquU+tT/73pHpNiolIpXSm1USm3x21YYKbaF\n2BijlNqslHr9lG3TWp/2H2AU8AYmhz/V33YOsAWIA7KBPUBMH9iWFPL9duAPkWIfMBuw+7//CvhV\npNjmt+NsYCxQAuSGtPe5fUCM/75nArF+e87p7d9RmE2XAl8Btoa0PQbc5/9+n/Vv3Ef2pQNf8X8f\nBOz0/1v2uY2AAgb6vzuA94GpkWBbiI13Ay8Ar5/qv21/8fB/DdwLhA5IzAde0lo3aa33AruBKb1t\nmNa6NmR3AEEb+9w+rfWbWmuPf3cDMDJSbPPb95nWekcLXZFg3xRgt9b6c621C3jJb1efobV+GzgS\n1jwfeM7//Tkgr1eNCkFrXam1/tD//TjwGZBBBNioDSf8uw7/R0eCbQBKqZHAXOBPIc2dtu20F3yl\n1EnK724AAAKlSURBVHygXGu9JawrA/gyZH+/v63XUUo9opT6EvgO8FN/c8TY5+cm4B/+75FmWziR\nYF8k2NARhmmtK/3fDwDD+tIYC6VUFjAZ40lHhI3+kMlHwCHgLa11xNgGPI5xan0hbZ227bRY4lAp\nVQwMb6HrQeABTGiiz2jLPq31a1rrB4EHlVL3A4uB/EixzX/Mg4AHeL637LLoiH1C96C11kqpPk/L\nU0oNBP4G3Km1rlVKBfr60kattRc4zz+W9apSakJYf5/YppS6Gjiktd6klJre0jEdte20EHyt9ayW\n2pVSEzEx3C3+P5qRwIdKqSlAOSa2bzHS39Zr9rXA88BqjOD3in3t2aaUuhG4Gpip/cHA3rKtI/a1\nQq/ZF+E2dISDSql0rXWlUiod4732GUopB0bsn9dav+JvjigbtdY1Sqm1wFURYtvFwDyl1BwgHkhS\nSv3lVGw7rUM6WutPtNZDtdZZWusszGv1V7TWB4CVwLVKqTilVDYwBtjY2zYqpcaE7M4Htvu/97l9\nSqmrMK+J87TW9SFdfW5bO0SCfR8AY5RS2UqpWOBav12RxkrgBv/3G4A+e2tSxit7BvhMa700pKvP\nbVRKpVlZakqpBOAKzP/VPrdNa32/1nqkX+OuBf6ptf7uKdnWVyPOPfEByvBn6fj3H8RkUuwAvtZH\nNv0N2Ap8DPwdyIgU+zCDnV8CH/k/f4gU2/w2LMA8xJuAg8AbEWbfHEymyR5MCKrXbQiz50WgEnD7\nf283A05gDbALKAaG9KF9X8UMhH4c8jc3JxJsBM4FNvtt2wr81N/e57aF2TmdYJZOp22TmbaCIAhR\nwmkd0hEEQRA6jgi+IAhClCCCLwiCECWI4AuCIEQJIviCIAhRggi+IAhClCCCLwiCECWI4AuCIEQJ\n/x/YecAxU8guKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f15a3697d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids, data, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Mean shift"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Most people that have come across clustering algorithms have learnt about **k-means**. Mean shift clustering is a newer and less well-known approach, but it has some important advantages:\n",
    "* It doesn't require selecting the number of clusters in advance, but instead just requires a **bandwidth** to be specified, which can be easily chosen automatically\n",
    "* It can handle clusters of any shape, whereas k-means (without using special extensions) requires that clusters be roughly ball shaped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The algorithm is as follows:\n",
    "* For each data point x in the sample X, find the distance between that point x and every other point in X\n",
    "* Create weights for each point in X by using the **Gaussian kernel** of that point's distance to x\n",
    "    * This weighting approach penalizes points further away from x\n",
    "    * The rate at which the weights fall to zero is determined by the **bandwidth**, which is the standard deviation of the Gaussian\n",
    "![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": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "So here's the definition of the gaussian kernel, which you may remember from high school..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def gaussian(d, bw):\n",
    "    return np.exp(-0.5*((d/bw))**2) / (bw*math.sqrt(2*math.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    " This person at the science march certainly remembered!\n",
    "\n",
    "<img src=\"http://i.imgur.com/nijQLHw.jpg\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In our implementation, we choose the bandwidth to be 2.5. \n",
    "\n",
    "One easy way to choose bandwidth is to find which bandwidth covers one third of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def meanshift(data):\n",
    "    X = np.copy(data)\n",
    "    for it in range(5):\n",
    "        for i, x in enumerate(X):\n",
    "            dist = np.sqrt(((x-X)**2).sum(1))\n",
    "            weight = gaussian(dist, 2.5)\n",
    "            X[i] = (np.expand_dims(weight,1)*X).sum(0) / weight.sum()\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 916 ms, sys: 0 ns, total: 916 ms\n",
      "Wall time: 914 ms\n"
     ]
    }
   ],
   "source": [
    "%time X=meanshift(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can see that mean shift clustering has almost reproduced our original clustering. The one exception are the very close clusters, but if we really wanted to differentiate them we could lower the bandwidth.\n",
    "\n",
    "What is impressive is that this algorithm nearly reproduced the original clusters without telling it how many clusters there should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFNxJREFUeJzt3W9sZNd53/HvU8lWy6gFI3rheiQ5VFsxqGwQSUQILVAU\nJOzESnbhZVoEUFCgjhKQiGHXW70JJAjwzCAQkNRAFyhSWyARbw1UiGAkDa1y49iSQcJv4jpUIy8k\n2xxvbMqSKNurEdjUIKBUztMXnBkN94/I3SE5MzzfDzDYM+fcvfeRQPx49cyZq8hMJEnH39/rdwGS\npKNh4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKcXO/C+j2rne9K8fHx/tdhiQN\nlWefffa1zDyx13EDFfjj4+Osra31uwxJGioR8eJ+jrOlI0mFKC7wm83mdc1L0nFRVODXajUmJydp\nNBq75huNBpOTk9Rqtf4UJklHoJjAr9Vq1Ot1Njc3mZmZ6YR+o9FgZmaGzc1N6vW6oS/p2Coi8Nth\n39YO/fPnz3fCvs3Ql3RcHfvAbzabLC4u7pqbZZbtzW1OnTrF9uY2s8zuWl9cXLSnL+nYOfaBPzY2\nxsrKCpVKBdgJ+zOc4SxnGWecs5zlDGc6oV+pVFhZWWFsbKyfZUvSgTv2gQ8wMTHRCf1VVtlgg3HG\nOcc5xhlngw1WWe2E/cTERL9LlqQDV0Tgw07oLywssMUWdeq71urU2WKLhYUFw17SkTnqbeLFBH6j\n0WB+fp5RRqlS3bVWpcooo8zPz1+xZVOSDkM/tokXEfjdWy+nme60cR7kwU57Z5rpK7ZsStJh6Ns2\n8cwcmNe9996bB+21117LSqWSQOc1y2yOMppAjjKas8zuWq9UKvnaa68deC2SVK1Wd+VNO3OWl5ev\nyCogq9XqnucE1nIfGXvs7/DHxsaYm5vbNbfEEiOVEZaXlxmpjLDE0q71ubk5d+lIOnB93ya+n98K\nR/U6jDv8tu7fqpVKJdfX1zMzc319fddv1f38NpWkG9WdObPM5goreY5zOc54nuNcrrDS6Tp0Z9Xb\nYZ93+LFz7GCYmprKw3w8cq1WY3Fx8Yqtl+2+2dzcnN+ylXTo2pmzvbnd+U5Q2wYbPMRDjFRG9r1N\nPCKezcypPY87roFf+cyPO+PNj97aGTebzau2a641L0mH4fz585w6darznaC29maS5eVlTp48ua9z\n7Tfwj30P/3LXCnXDXtJR6dc28eICX5L6qZ/bxI9tS0eSBk2z2WRycnLXE3pnmWWVVbbYYpRRppne\ntXOwUqlw4cKFt+1C2NKRpAHT723iBr4kHaFarUa1+lbfvv3QxpMnT+56si9AtVo90J2DNx/YmSRJ\n+9IO8cu3ibef7HtY28Tt4UvSYYp4a3xZ3h7UNnF7+JI04I56m7iBL0mFsIcvSYdpgNrmB3KHHxGf\njYgfRcTzXXO3RcTTEfGd1p8/fRDXkiTdmINq6fw34P7L5h4GvpKZdwNfab2XJPXJgQR+Zn4VeP2y\n6dPA51rjz8FlD3mWJB2pw/zQ9t2Z+Wpr/APg3Yd4LUnSHo5kl07rAf1X/eQiIuYjYi0i1i5dunQU\n5UhSkQ4z8H8YEe8BaP35o6sdlJkLmTmVmVMnTpw4xHIkqWyHGfhPAR9pjT8CfOEQryVJ2sNBbcv8\nI+AvgJ+NiJcj4reA3wN+MSK+A3yw9V6S1CcH8sWrzPz1ayx94CDOL0nqnY9WkKRCGPiSVAgDX5IK\nYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAG\nviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBL\nUiEMfGlANJvN65qXrpeBLw2AWq3G5OQkjUZj13yj0WBycpJardafwnSsGPhSn9VqNer1Opubm8zM\nzHRCv9FoMDMzw+bmJvV63dBXzw498CPi/ohYj4iLEfHwYV9PGibtsG9rh/758+c7Yd9m6KtXhxr4\nEXET8F+BXwbuAX49Iu45zGtKw6LZbLK4uLhrbpZZtje3OXXqFNub28wyu2t9cXHRnr5u2GHf4d8H\nXMzM72bm3wJPAqcP+ZrSUBgbG2NlZYVKpQLshP0ZznCWs4wzzlnOcoYzndCvVCqsrKwwNjbWz7I1\nxA478G8HXup6/3JrThIwMTHRCf1VVtlgg3HGOcc5xhlngw1WWe2E/cTERL9L1hDr+4e2ETEfEWsR\nsXbp0qV+lyMduYmJCRYWFthiizr1XWt16myxxcLCgmGvnh124L8C3Nn1/o7WXEdmLmTmVGZOnThx\n4pDLkQZPo9Fgfn6eUUapUt21VqXKKKPMz89fsWVTul6HHfh/CdwdEXdFxDuBB4CnDvma0tDo3no5\nzXSnjfMgD3baO9NMX7FlU7oRNx/myTPzzYj4OPAl4Cbgs5n5wmFeUxoWzWZz19bLJZYAWGWVLbZ4\niIeYZroz3w79Cxcu+MGtbsih9/Az888ycyIz/2lmPnbY15OGxdjYGHNzc7vmllhipDLC8vIyI5WR\nTti3zc3NGfa6YX3/0FYqWa1Wo1p9q2/f3o1z8uTJXVs2AarVql+8Uk8OtaUjaW/tEF9cXNy19bK9\nZXNmZoa5uTnDXj2LzOx3DR1TU1O5trbW7zKkvmg2m1dt11xrXmqLiGczc2qv47zDl47Yh3/8V53x\nU7f+fGd8rVA37HVQ7OFLUiEMfEkqhC0d6Yh1t3Gko+QdviQVwsCXpEIY+JJUCANfkgph4EtSIQx8\nSSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJek\nQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRA9BX5E/FpEvBARfxcRU5etPRIRFyNiPSI+1FuZ\nkq5Xs9m8rnkdf73e4T8P/Bvgq92TEXEP8ADwPuB+4NMRcVOP15K0T7VajcnJSRqNxq75RqPB5OQk\ntVqtP4Wpr3oK/Mz8VmauX2XpNPBkZr6Rmd8DLgL39XItSftTq9Wo1+tsbm4yMzPTCf1Go8HMzAyb\nm5vU63VDv0CH1cO/HXip6/3LrTlJh6gd9m3t0D9//nwn7NsM/fLsGfgR8UxEPH+V1+mDKCAi5iNi\nLSLWLl26dBCnlIrUbDZZXFzcNTfLLNub25w6dYrtzW1mmd21vri4aE+/IDfvdUBmfvAGzvsKcGfX\n+ztac1c7/wKwADA1NZU3cC1JwNjYGCsrK507+VlmOcMZTnOaOnWqVBlnHIAllqhUKqysrDA2Ntbf\nwnVkDqul8xTwQETcEhF3AXcDXz+ka0lqmZiYYGVlhUqlwiqrbLDBOOOc4xzjjLPBBqusdsJ+YmKi\n3yXrCPW6LfNXI+Jl4F8C5yPiSwCZ+QLweeCbwJ8DH8vMn/RarKS9TUxMsLCwwBZb1KnvWqtTZ4st\nFhYWDPsCRebgdFGmpqZybW2t32VIQ629G2d7c5uznO20cQA22OAhHmKkMuId/jESEc9m5tRex/lN\nW+kY6d56Oc10p43zIA922jvTTF+xZVNl8A5fOiaazSaTk5O7tl7OMssqq2yxxSijTDPNEkud9Uql\nwoULF/zgdsh5hy8VZmxsjLm5uV1zSywxUhlheXmZkcrIrrAHmJubM+wLYuBLQ+oTs69z+gfP8YnZ\n1ztztVqNarXaed/ejXPy5MnO7p22arXqF68Ks+c+fEmD6cXHv0/emrz4+PeB2zrz7RBfXFzc9cFs\ne8vmzMwMc3Nzhn2B7OFLQ+oTs6/z4uPf52d++738l6XbrlhvNptXbddca17Da789fO/wpSG1E/K3\ncVlbvuNaoW7Yl8seviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RC\nGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSB\nL0mFMPAlqRA9BX5EfCoivh0RFyLiTyNitGvtkYi4GBHrEfGh3kuVJPWi1zv8p4H3Z+Yk0AAeAYiI\ne4AHgPcB9wOfjoiberyWJKkHPQV+Zn45M99svf0acEdrfBp4MjPfyMzvAReB+3q5liSpNwfZw/9N\n4Iut8e3AS11rL7fmJEl9cvNeB0TEM8A/vsrSo5n5hdYxjwJvAk9cbwERMQ/MA7z3ve+93r8uSdqn\nPQM/Mz/4dusR8RvAKeADmZmt6VeAO7sOu6M1d7XzLwALAFNTU3m1YyRJvet1l879wO8AH87M7a6l\np4AHIuKWiLgLuBv4ei/XkiT1Zs87/D38AXAL8HREAHwtM387M1+IiM8D32Sn1fOxzPxJj9eSJPWg\np8DPzH/2NmuPAY/1cn5J0sHxm7aSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4\nklQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSpEMYHf\nbDava16SjpsiAr9WqzE5OUmj0dg132g0mJycpFar9acwSTpCxz7wa7Ua9Xqdzc1NZmZmOqHfaDSY\nmZlhc3OTer1u6Es69o514LfDvq0d+ufPn++EfZuhL+m4O7aB32w2WVxc3DU3yyzbm9ucOnWK7c1t\nZpndtb64uGhPX9KxdWwDf2xsjJWVFSqVCrAT9mc4w1nOMs44ZznLGc50Qr9SqbCyssLY2Fg/y5ak\nQ3NsAx9gYmKiE/qrrLLBBuOMc45zjDPOBhusstoJ+4mJiX6XLEmH5lgHPuyE/sLCAltsUae+a61O\nnS22WFhYMOwlHXvHPvAbjQbz8/OMMkqV6q61KlVGGWV+fv6KLZuSdNwc68Dv3no5zXSnjfMgD3ba\nO9NMX7FlU5KOo8jMftfQMTU1lWtrawdyrmazyeTk5K6tl7PMssoqW2wxyijTTLPEUme9Uqlw4cIF\nP7iVNFQi4tnMnNrruJ7u8CPidyPiQkQ8FxFfjohK19ojEXExItYj4kO9XOdGjI2NMTc3t2tuiSVG\nKiMsLy8zUhnZFfYAc3Nzhr2kY6vXls6nMnMyM38OWAY+CRAR9wAPAO8D7gc+HRE39Xit61ar1ahW\n3+rbt3fjnDx5cteWTYBqteoXryQdazf38pcz82+63v4U0O4PnQaezMw3gO9FxEXgPuAvernejWiH\n+OLi4q6tl+0tmzMzM8zNzRn2ko69nnv4EfEY8O+B/wPMZOaliPgD4GuZ+d9bx/wh8MXM/OO3O1dP\nPfzPdH1r9qNLVyw3m82rtmuuNS9Jw+LAevgR8UxEPH+V12mAzHw0M+8EngA+fgOFzkfEWkSsXbp0\n6Xr/+r5dK9QNe0k3atgeu75n4GfmBzPz/Vd5feGyQ58A/m1r/ApwZ9faHa25q51/ITOnMnPqxIkT\nN/LPIElHbigfu56ZN/wC7u4a/wfgj1vj9wHfAG4B7gK+C9y01/nuvffelKRBV61Wk53PLLNSqeT6\n+npmZq6vr2elUumsVavVI6kHWMt9ZHavu3R+r9XeuQD8EnCm9UvkBeDzwDeBPwc+lpk/6fFaktR3\nw/zY9WP7xStJOmiD+oXOI/nilSSVZNgfu27gS9J1GObHrhv4knSdhvWx6wa+JF2nYX3suoEvSddh\nmB+77i4dSdond+lIUiGG/bHrBr4kXYdhfux6T49HlqQSDetj1+3hS9INGpTHru+3h+8dviTtYfnH\nn+mMT9360c542B67bg9fkgph4EtSIWzpSNIeuts4w8w7fEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4\nklQIA1+SCjFQz9KJiEvAi326/LuA1/p07V4Ma90wvLVb99Eb1tqPqu6fycwTex00UIHfTxGxtp+H\nDw2aYa0bhrd26z56w1r7oNVtS0eSCmHgS1IhDPy3LPS7gBs0rHXD8NZu3UdvWGsfqLrt4UtSIbzD\nl6RCFB/4EfG7EXEhIp6LiC9HRKVr7ZGIuBgR6xHxoX7WebmI+FREfLtV+59GxGjX2iDX/WsR8UJE\n/F1ETF22NrB1t0XE/a36LkbEw/2u51oi4rMR8aOIeL5r7raIeDoivtP686f7WePVRMSdEbESEd9s\n/Zycac0PdO0R8fcj4usR8Y1W3fXW/GDVnZlFv4B/1DX+BPB4a3wP8A3gFuAu4K+Bm/pdb1etvwTc\n3Br/PvD7Q1L3Pwd+FlgFprrmB7ruVo03ter6J8A7W/Xe0++6rlHrvwZ+AXi+a+4/AQ+3xg+3f2YG\n6QW8B/iF1vgfAo3Wz8ZA1w4EcGtr/A7gfwH/YtDqLv4OPzP/puvtTwHtDzVOA09m5huZ+T3gInDf\nUdd3LZn55cx8s/X2a8AdrfGg1/2tzFy/ytJA191yH3AxM7+bmX8LPMlO3QMnM78KvH7Z9Gngc63x\n54DZIy1qHzLz1cz8363x/wW+BdzOgNeeO37cevuO1isZsLqLD3yAiHgsIl4C/h3wydb07cBLXYe9\n3JobRL8JfLE1Hqa6uw1D3cNQ49t5d2a+2hr/AHh3P4vZS0SMAz/Pzt3ywNceETdFxHPAj4CnM3Pg\n6i4i8CPimYh4/iqv0wCZ+Whm3gk8AXy8v9W+Za+6W8c8CrzJTu0DYT91q79yp8cwsFv0IuJW4E+A\n/3jZf4UPbO2Z+ZPM/Dl2/mv7voh4/2Xrfa+7iP+nbWZ+cJ+HPgH8GVAFXgHu7Fq7ozV3ZPaqOyJ+\nAzgFfKD1wwRDUPc19L3ufRiGGt/ODyPiPZn5akS8h5070YETEe9gJ+yfyMz/0ZoeitoBMnMrIlaA\n+xmwuou4w387EXF319vTwLdb46eAByLiloi4C7gb+PpR13ctEXE/8DvAhzNzu2tpoOt+G8NQ918C\nd0fEXRHxTuABduoeFk8BH2mNPwJ8oY+1XFVEBPCHwLcy8z93LQ107RFxor1TLiL+AfCL7GTJYNXd\n70+3+/1i507ieeAC8D+B27vWHmVnV8Y68Mv9rvWyui+y009+rvV6fEjq/lV2et9vAD8EvjQMdXfV\n+Cvs7Bz5a+DRftfzNnX+EfAq8P9a/75/CxgDvgJ8B3gGuK3fdV6l7n/FTtvjQtfP9q8Meu3AJPBX\nrbqfBz7Zmh+ouv2mrSQVoviWjiSVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQ/x9r\nadY2FGvLNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f15a35d3f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We should be able to accelerate this algorithm with a GPU."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Broadcasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1, 2, 3]), (3,))"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v=np.array([1,2,3]); v, v.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[1, 2, 3],\n",
       "        [2, 4, 6],\n",
       "        [3, 6, 9]]), (3, 3))"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m=np.array([v,v*2,v*3]); m, m.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2,  4,  6],\n",
       "       [ 3,  6,  9],\n",
       "       [ 4,  8, 12]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m+v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[1],\n",
       "        [2],\n",
       "        [3]]), (3, 1))"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1=np.expand_dims(v,-1); v1, v1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2,  3,  4],\n",
       "       [ 4,  6,  8],\n",
       "       [ 6,  9, 12]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m+v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Pytorch does not support broadcasting, therefore I have replaced the operators with broadcasting versions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def unit_prefix(x, n=1):\n",
    "    for i in range(n): x = x.unsqueeze(0)\n",
    "    return x\n",
    "\n",
    "def align(x, y, start_dim=2):\n",
    "    xd, yd = x.dim(), y.dim()\n",
    "    if xd > yd: y = unit_prefix(y, xd - yd)\n",
    "    elif yd > xd: x = unit_prefix(x, yd - xd)\n",
    "\n",
    "    xs, ys = list(x.size()), list(y.size())\n",
    "    nd = len(ys)\n",
    "    for i in range(start_dim, nd):\n",
    "        td = nd-i-1\n",
    "        if   ys[td]==1: ys[td] = xs[td]\n",
    "        elif xs[td]==1: xs[td] = ys[td]\n",
    "    return x.expand(*xs), y.expand(*ys)\n",
    "\n",
    "def aligned_op(x,y,f): return f(*align(x,y,0))\n",
    "\n",
    "def add(x, y): return aligned_op(x, y, operator.add)\n",
    "def sub(x, y): return aligned_op(x, y, operator.sub)\n",
    "def mul(x, y): return aligned_op(x, y, operator.mul)\n",
    "def div(x, y): return aligned_op(x, y, operator.truediv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## GPU-accelerated mean shift in pytorch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "One advantage of pytorch is that it's very similar to numpy. For instance, the definition of `gaussian` is identical, except for the namespace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def gaussian(d, bw):\n",
    "    return torch.exp(-0.5*((d/bw))**2) / (bw*math.sqrt(2*math.pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "And the implementation of meanshift is nearly identical too!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def meanshift(data):\n",
    "    X = torch.FloatTensor(np.copy(data))\n",
    "    for it in range(5):\n",
    "        for i, x in enumerate(X):\n",
    "            dist = torch.sqrt((sub(x, X)**2).sum(1))\n",
    "            weight = gaussian(dist, 3)\n",
    "            num = mul(weight, X).sum(0)\n",
    "            X[i] = num / weight.sum()\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "This implementation actually takes longer. Oh dear! What do you think is causing this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.02 s, sys: 0 ns, total: 2.02 s\n",
      "Wall time: 1.9 s\n"
     ]
    }
   ],
   "source": [
    "%time X = meanshift(data).numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "All the computation is happening in the <tt>for</tt> loop, which isn't accelerated by pytorch. Each iteration launches a new cuda kernel, which takes time and slows the algorithm down as a whole. Furthermore, each iteration doesn't have enough processing to do to fill up all of the threads of the GPU. But at least the results are correct..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFNxJREFUeJzt3W9sZNd53/HvU8lWy6gFI3rheiQ5VFsxqGwQSUQILVAU\nJOzESnbhZVoEUFCgjhKQiGHXW70JJAjwzCAQkNRAFyhSWyARbw1UiGAkDa1y49iSQcJv4jpUIy8k\n2xxvbMqSKNurEdjUIKBUztMXnBkN94/I3SE5MzzfDzDYM+fcvfeRQPx49cyZq8hMJEnH39/rdwGS\npKNh4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKcXO/C+j2rne9K8fHx/tdhiQN\nlWefffa1zDyx13EDFfjj4+Osra31uwxJGioR8eJ+jrOlI0mFKC7wm83mdc1L0nFRVODXajUmJydp\nNBq75huNBpOTk9Rqtf4UJklHoJjAr9Vq1Ot1Njc3mZmZ6YR+o9FgZmaGzc1N6vW6oS/p2Coi8Nth\n39YO/fPnz3fCvs3Ql3RcHfvAbzabLC4u7pqbZZbtzW1OnTrF9uY2s8zuWl9cXLSnL+nYOfaBPzY2\nxsrKCpVKBdgJ+zOc4SxnGWecs5zlDGc6oV+pVFhZWWFsbKyfZUvSgTv2gQ8wMTHRCf1VVtlgg3HG\nOcc5xhlngw1WWe2E/cTERL9LlqQDV0Tgw07oLywssMUWdeq71urU2WKLhYUFw17SkTnqbeLFBH6j\n0WB+fp5RRqlS3bVWpcooo8zPz1+xZVOSDkM/tokXEfjdWy+nme60cR7kwU57Z5rpK7ZsStJh6Ns2\n8cwcmNe9996bB+21117LSqWSQOc1y2yOMppAjjKas8zuWq9UKvnaa68deC2SVK1Wd+VNO3OWl5ev\nyCogq9XqnucE1nIfGXvs7/DHxsaYm5vbNbfEEiOVEZaXlxmpjLDE0q71ubk5d+lIOnB93ya+n98K\nR/U6jDv8tu7fqpVKJdfX1zMzc319fddv1f38NpWkG9WdObPM5goreY5zOc54nuNcrrDS6Tp0Z9Xb\nYZ93+LFz7GCYmprKw3w8cq1WY3Fx8Yqtl+2+2dzcnN+ylXTo2pmzvbnd+U5Q2wYbPMRDjFRG9r1N\nPCKezcypPY87roFf+cyPO+PNj97aGTebzau2a641L0mH4fz585w6darznaC29maS5eVlTp48ua9z\n7Tfwj30P/3LXCnXDXtJR6dc28eICX5L6qZ/bxI9tS0eSBk2z2WRycnLXE3pnmWWVVbbYYpRRppne\ntXOwUqlw4cKFt+1C2NKRpAHT723iBr4kHaFarUa1+lbfvv3QxpMnT+56si9AtVo90J2DNx/YmSRJ\n+9IO8cu3ibef7HtY28Tt4UvSYYp4a3xZ3h7UNnF7+JI04I56m7iBL0mFsIcvSYdpgNrmB3KHHxGf\njYgfRcTzXXO3RcTTEfGd1p8/fRDXkiTdmINq6fw34P7L5h4GvpKZdwNfab2XJPXJgQR+Zn4VeP2y\n6dPA51rjz8FlD3mWJB2pw/zQ9t2Z+Wpr/APg3Yd4LUnSHo5kl07rAf1X/eQiIuYjYi0i1i5dunQU\n5UhSkQ4z8H8YEe8BaP35o6sdlJkLmTmVmVMnTpw4xHIkqWyHGfhPAR9pjT8CfOEQryVJ2sNBbcv8\nI+AvgJ+NiJcj4reA3wN+MSK+A3yw9V6S1CcH8sWrzPz1ayx94CDOL0nqnY9WkKRCGPiSVAgDX5IK\nYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAG\nviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBL\nUiEMfGlANJvN65qXrpeBLw2AWq3G5OQkjUZj13yj0WBycpJardafwnSsGPhSn9VqNer1Opubm8zM\nzHRCv9FoMDMzw+bmJvV63dBXzw498CPi/ohYj4iLEfHwYV9PGibtsG9rh/758+c7Yd9m6KtXhxr4\nEXET8F+BXwbuAX49Iu45zGtKw6LZbLK4uLhrbpZZtje3OXXqFNub28wyu2t9cXHRnr5u2GHf4d8H\nXMzM72bm3wJPAqcP+ZrSUBgbG2NlZYVKpQLshP0ZznCWs4wzzlnOcoYzndCvVCqsrKwwNjbWz7I1\nxA478G8HXup6/3JrThIwMTHRCf1VVtlgg3HGOcc5xhlngw1WWe2E/cTERL9L1hDr+4e2ETEfEWsR\nsXbp0qV+lyMduYmJCRYWFthiizr1XWt16myxxcLCgmGvnh124L8C3Nn1/o7WXEdmLmTmVGZOnThx\n4pDLkQZPo9Fgfn6eUUapUt21VqXKKKPMz89fsWVTul6HHfh/CdwdEXdFxDuBB4CnDvma0tDo3no5\nzXSnjfMgD3baO9NMX7FlU7oRNx/myTPzzYj4OPAl4Cbgs5n5wmFeUxoWzWZz19bLJZYAWGWVLbZ4\niIeYZroz3w79Cxcu+MGtbsih9/Az888ycyIz/2lmPnbY15OGxdjYGHNzc7vmllhipDLC8vIyI5WR\nTti3zc3NGfa6YX3/0FYqWa1Wo1p9q2/f3o1z8uTJXVs2AarVql+8Uk8OtaUjaW/tEF9cXNy19bK9\nZXNmZoa5uTnDXj2LzOx3DR1TU1O5trbW7zKkvmg2m1dt11xrXmqLiGczc2qv47zDl47Yh3/8V53x\nU7f+fGd8rVA37HVQ7OFLUiEMfEkqhC0d6Yh1t3Gko+QdviQVwsCXpEIY+JJUCANfkgph4EtSIQx8\nSSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJek\nQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRA9BX5E/FpEvBARfxcRU5etPRIRFyNiPSI+1FuZ\nkq5Xs9m8rnkdf73e4T8P/Bvgq92TEXEP8ADwPuB+4NMRcVOP15K0T7VajcnJSRqNxq75RqPB5OQk\ntVqtP4Wpr3oK/Mz8VmauX2XpNPBkZr6Rmd8DLgL39XItSftTq9Wo1+tsbm4yMzPTCf1Go8HMzAyb\nm5vU63VDv0CH1cO/HXip6/3LrTlJh6gd9m3t0D9//nwn7NsM/fLsGfgR8UxEPH+V1+mDKCAi5iNi\nLSLWLl26dBCnlIrUbDZZXFzcNTfLLNub25w6dYrtzW1mmd21vri4aE+/IDfvdUBmfvAGzvsKcGfX\n+ztac1c7/wKwADA1NZU3cC1JwNjYGCsrK507+VlmOcMZTnOaOnWqVBlnHIAllqhUKqysrDA2Ntbf\nwnVkDqul8xTwQETcEhF3AXcDXz+ka0lqmZiYYGVlhUqlwiqrbLDBOOOc4xzjjLPBBqusdsJ+YmKi\n3yXrCPW6LfNXI+Jl4F8C5yPiSwCZ+QLweeCbwJ8DH8vMn/RarKS9TUxMsLCwwBZb1KnvWqtTZ4st\nFhYWDPsCRebgdFGmpqZybW2t32VIQ629G2d7c5uznO20cQA22OAhHmKkMuId/jESEc9m5tRex/lN\nW+kY6d56Oc10p43zIA922jvTTF+xZVNl8A5fOiaazSaTk5O7tl7OMssqq2yxxSijTDPNEkud9Uql\nwoULF/zgdsh5hy8VZmxsjLm5uV1zSywxUhlheXmZkcrIrrAHmJubM+wLYuBLQ+oTs69z+gfP8YnZ\n1ztztVqNarXaed/ejXPy5MnO7p22arXqF68Ks+c+fEmD6cXHv0/emrz4+PeB2zrz7RBfXFzc9cFs\ne8vmzMwMc3Nzhn2B7OFLQ+oTs6/z4uPf52d++738l6XbrlhvNptXbddca17Da789fO/wpSG1E/K3\ncVlbvuNaoW7Yl8seviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RC\nGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSB\nL0mFMPAlqRA9BX5EfCoivh0RFyLiTyNitGvtkYi4GBHrEfGh3kuVJPWi1zv8p4H3Z+Yk0AAeAYiI\ne4AHgPcB9wOfjoiberyWJKkHPQV+Zn45M99svf0acEdrfBp4MjPfyMzvAReB+3q5liSpNwfZw/9N\n4Iut8e3AS11rL7fmJEl9cvNeB0TEM8A/vsrSo5n5hdYxjwJvAk9cbwERMQ/MA7z3ve+93r8uSdqn\nPQM/Mz/4dusR8RvAKeADmZmt6VeAO7sOu6M1d7XzLwALAFNTU3m1YyRJvet1l879wO8AH87M7a6l\np4AHIuKWiLgLuBv4ei/XkiT1Zs87/D38AXAL8HREAHwtM387M1+IiM8D32Sn1fOxzPxJj9eSJPWg\np8DPzH/2NmuPAY/1cn5J0sHxm7aSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4\nklQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSpEMYHf\nbDava16SjpsiAr9WqzE5OUmj0dg132g0mJycpFar9acwSTpCxz7wa7Ua9Xqdzc1NZmZmOqHfaDSY\nmZlhc3OTer1u6Es69o514LfDvq0d+ufPn++EfZuhL+m4O7aB32w2WVxc3DU3yyzbm9ucOnWK7c1t\nZpndtb64uGhPX9KxdWwDf2xsjJWVFSqVCrAT9mc4w1nOMs44ZznLGc50Qr9SqbCyssLY2Fg/y5ak\nQ3NsAx9gYmKiE/qrrLLBBuOMc45zjDPOBhusstoJ+4mJiX6XLEmH5lgHPuyE/sLCAltsUae+a61O\nnS22WFhYMOwlHXvHPvAbjQbz8/OMMkqV6q61KlVGGWV+fv6KLZuSdNwc68Dv3no5zXSnjfMgD3ba\nO9NMX7FlU5KOo8jMftfQMTU1lWtrawdyrmazyeTk5K6tl7PMssoqW2wxyijTTLPEUme9Uqlw4cIF\nP7iVNFQi4tnMnNrruJ7u8CPidyPiQkQ8FxFfjohK19ojEXExItYj4kO9XOdGjI2NMTc3t2tuiSVG\nKiMsLy8zUhnZFfYAc3Nzhr2kY6vXls6nMnMyM38OWAY+CRAR9wAPAO8D7gc+HRE39Xit61ar1ahW\n3+rbt3fjnDx5cteWTYBqteoXryQdazf38pcz82+63v4U0O4PnQaezMw3gO9FxEXgPuAvernejWiH\n+OLi4q6tl+0tmzMzM8zNzRn2ko69nnv4EfEY8O+B/wPMZOaliPgD4GuZ+d9bx/wh8MXM/OO3O1dP\nPfzPdH1r9qNLVyw3m82rtmuuNS9Jw+LAevgR8UxEPH+V12mAzHw0M+8EngA+fgOFzkfEWkSsXbp0\n6Xr/+r5dK9QNe0k3atgeu75n4GfmBzPz/Vd5feGyQ58A/m1r/ApwZ9faHa25q51/ITOnMnPqxIkT\nN/LPIElHbigfu56ZN/wC7u4a/wfgj1vj9wHfAG4B7gK+C9y01/nuvffelKRBV61Wk53PLLNSqeT6\n+npmZq6vr2elUumsVavVI6kHWMt9ZHavu3R+r9XeuQD8EnCm9UvkBeDzwDeBPwc+lpk/6fFaktR3\nw/zY9WP7xStJOmiD+oXOI/nilSSVZNgfu27gS9J1GObHrhv4knSdhvWx6wa+JF2nYX3suoEvSddh\nmB+77i4dSdond+lIUiGG/bHrBr4kXYdhfux6T49HlqQSDetj1+3hS9INGpTHru+3h+8dviTtYfnH\nn+mMT9360c542B67bg9fkgph4EtSIWzpSNIeuts4w8w7fEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4\nklQIA1+SCjFQz9KJiEvAi326/LuA1/p07V4Ma90wvLVb99Eb1tqPqu6fycwTex00UIHfTxGxtp+H\nDw2aYa0bhrd26z56w1r7oNVtS0eSCmHgS1IhDPy3LPS7gBs0rHXD8NZu3UdvWGsfqLrt4UtSIbzD\nl6RCFB/4EfG7EXEhIp6LiC9HRKVr7ZGIuBgR6xHxoX7WebmI+FREfLtV+59GxGjX2iDX/WsR8UJE\n/F1ETF22NrB1t0XE/a36LkbEw/2u51oi4rMR8aOIeL5r7raIeDoivtP686f7WePVRMSdEbESEd9s\n/Zycac0PdO0R8fcj4usR8Y1W3fXW/GDVnZlFv4B/1DX+BPB4a3wP8A3gFuAu4K+Bm/pdb1etvwTc\n3Br/PvD7Q1L3Pwd+FlgFprrmB7ruVo03ter6J8A7W/Xe0++6rlHrvwZ+AXi+a+4/AQ+3xg+3f2YG\n6QW8B/iF1vgfAo3Wz8ZA1w4EcGtr/A7gfwH/YtDqLv4OPzP/puvtTwHtDzVOA09m5huZ+T3gInDf\nUdd3LZn55cx8s/X2a8AdrfGg1/2tzFy/ytJA191yH3AxM7+bmX8LPMlO3QMnM78KvH7Z9Gngc63x\n54DZIy1qHzLz1cz8363x/wW+BdzOgNeeO37cevuO1isZsLqLD3yAiHgsIl4C/h3wydb07cBLXYe9\n3JobRL8JfLE1Hqa6uw1D3cNQ49t5d2a+2hr/AHh3P4vZS0SMAz/Pzt3ywNceETdFxHPAj4CnM3Pg\n6i4i8CPimYh4/iqv0wCZ+Whm3gk8AXy8v9W+Za+6W8c8CrzJTu0DYT91q79yp8cwsFv0IuJW4E+A\n/3jZf4UPbO2Z+ZPM/Dl2/mv7voh4/2Xrfa+7iP+nbWZ+cJ+HPgH8GVAFXgHu7Fq7ozV3ZPaqOyJ+\nAzgFfKD1wwRDUPc19L3ufRiGGt/ODyPiPZn5akS8h5070YETEe9gJ+yfyMz/0ZoeitoBMnMrIlaA\n+xmwuou4w387EXF319vTwLdb46eAByLiloi4C7gb+PpR13ctEXE/8DvAhzNzu2tpoOt+G8NQ918C\nd0fEXRHxTuABduoeFk8BH2mNPwJ8oY+1XFVEBPCHwLcy8z93LQ107RFxor1TLiL+AfCL7GTJYNXd\n70+3+/1i507ieeAC8D+B27vWHmVnV8Y68Mv9rvWyui+y009+rvV6fEjq/lV2et9vAD8EvjQMdXfV\n+Cvs7Bz5a+DRftfzNnX+EfAq8P9a/75/CxgDvgJ8B3gGuK3fdV6l7n/FTtvjQtfP9q8Meu3AJPBX\nrbqfBz7Zmh+ouv2mrSQVoviWjiSVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQ/x9r\nadY2FGvLNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f15797283c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## GPU batched algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "To truly accelerate the algorithm, we need to be performing updates on a batch of points per iteration, instead of just one as we were doing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def dist_b(a,b):\n",
    "    return torch.sqrt((sub(a.unsqueeze(0),b.unsqueeze(1))**2).sum(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 0.2897  0.4660  0.6124\n",
       " 0.7822  0.2356  0.2745\n",
       "[torch.FloatTensor of size 2x3]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=torch.rand(2,2)\n",
    "b=torch.rand(3,2)\n",
    "dist_b(b, a).squeeze(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# def gaussian(d, bw):\n",
    "#     return torch.exp(-0.5*((d/bw))**2) / (bw*math.sqrt(2*math.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def sum_sqz(a,axis): return a.sum(axis).squeeze(axis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def meanshift(data, bs=500):\n",
    "    n = len(data)\n",
    "    X = torch.FloatTensor(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(dist_b(X, X[s]), 2)\n",
    "            num = sum_sqz(mul(weight, X), 1)\n",
    "            X[s] = div(num, sum_sqz(weight, 1))\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "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": 50,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 40 ms, sys: 4 ms, total: 44 ms\n",
      "Wall time: 44 ms\n"
     ]
    }
   ],
   "source": [
    "%time X = meanshift(data).cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "That's more like it! We've gone from 914ms to 44ms, which is a speedup of over 2000%. Oh, and it even gives the right answer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFOFJREFUeJzt3W9sXfd93/H3d3biTvUCNoyQ5drK6GViMScg2powNmAY\nSMRt3EqI2K0FXKxY6hYkFiSb5ieFDQO596Iw1i7A9KRLChKNFqBGjSBbFU9q/tgBiTxpltJrIthJ\ndK0kLGzTTeRbcFlA1J3T7x7w3mteiRIpXZL3z+/9Ag507u93dM7XBvHh8ff+znFkJpKk0ff3+l2A\nJOlwGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQtze7wK2e8c73pETExP9LkOS\nhspzzz33WmYe3e24gQr8iYkJVldX+12GJA2ViPjLvRxnS0eSClFc4DebzZsal6RRUVTg12o1pqam\naDQaXeONRoOpqSlqtVp/CpOkQ1BM4NdqNer1Ouvr68zOznZCv9FoMDs7y/r6OvV63dCXNLKKCPx2\n2Le1Q//ChQudsG8z9CWNqpEP/GazydLSUtfYHHNsrm9y8uRJNtc3mWOua35pacmevqSRM/KBPz4+\nzvLyMpVKBdgK+9Oc5gxnmGCCM5zhNKc7oV+pVFheXmZ8fLyfZUvSvhv5wAeYnJzshP4KK6yxxgQT\nnOUsE0ywxhorrHTCfnJyst8lS9K+KyLwYSv0FxcX2WCDOvWuuTp1NthgcXHRsJd0aA57mXgxgd9o\nNFhYWGCMMapUu+aqVBljjIWFhWuWbErSQejHMvEiAn/70ssZZjptnId5uNPemWHmmiWbknQQ+rZM\nPDMHZrvvvvtyv7322mtZqVQS6GxzzOUYYwnkGGM5x1zXfKVSyddee23fa5GkarXalTftzDl//vw1\nWQVktVrd9ZzAau4hY0f+Dn98fJz5+fmusXOc40jlCOfPn+dI5QjnONc1Pz8/7yodSfuu78vE9/Jb\n4bC2g7jDb9v+W7VSqeSlS5cyM/PSpUtdv1X38ttUkm7V9syZYy6XWc6znM0JJvIsZ3OZ5U7XYXtW\n3Qh7vMOPrWMHw/T0dB7k65FrtRpLS0vXLL1s983m5+d9ylbSgWtnzub6ZueZoLY11niERzhSObLn\nZeIR8VxmTu963CgGfuWTP+rsr3/4zq65ZrO5Y7vmeuOSdBAuXLjAyZMnO88EtbUXk5w/f54TJ07s\n6Vx7DfyR7+Ff7XqhbthLOiz9WiZeXOBLUj/1c5n4SLZ0JGkQNZtNpqamut7QO8ccK6ywwQZjjDHD\nTNfKwUqlwsWLF2/YhbClI0kDpt/LxA18STpEtVqNavXNvn37pY0nTpzoerMvQLVa3deVg7fv25kk\nSdeKeHO/1UJvh/jVy8Tbb/Y9qGXi9vAl6SDtEPht+7VM3B6+JA24w14mbktHkg7SAHVRvMOXpELs\nS+BHxKci4gcR8fy2sbdHxDMR8WLrz5/aj2tJkm7Nft3h/zfgwavGHgW+nJnHgS+3PkuS+mRfAj8z\nvwL89VXDp4BPt/Y/DVe95FmSdKgOsof/zsx8tbX/V8A7D/BakqRdHMqXtq0X9O/4VXVELETEakSs\nXrly5TDKkaQiHWTgfz8i3gXQ+vMHOx2UmYuZOZ2Z00ePHj3AciSpbAcZ+E8DH2rtfwj43AFeS5K0\ni/1alvnHwJ8BPx0RL0fEbwG/C/x8RLwIPND6LEnqk3150jYzf+06U+/fj/NLknrnk7aSVAgDX5IK\nYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAG\nviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBL\nUiEMfEkqhIEvDYhms3lT49LNMvClAVCr1ZiamqLRaHSNNxoNpqamqNVq/SlMI8XAl/qsVqtRr9dZ\nX19ndna2E/qNRoPZ2VnW19ep1+uGvnp24IEfEQ9GxKWIuBwRjx709aRh0g77tnboX7hwoRP2bYa+\nenWggR8RtwH/FfhF4F7g1yLi3oO8pjQsms0mS0tLXWNzzLG5vsnJkyfZXN9kjrmu+aWlJXv6umUH\nfYd/P3A5M7+bmX8LPAWcOuBrSkNhfHyc5eVlKpUKsBX2pznNGc4wwQRnOMNpTndCv1KpsLy8zPj4\neD/L1hA76MC/C3hp2+eXW2OSgMnJyU7or7DCGmtMMMFZzjLBBGusscJKJ+wnJyf7XbKGWN+/tI2I\nhYhYjYjVK1eu9Lsc6dBNTk6yuLjIBhvUqXfN1amzwQaLi4uGvXp20IH/CnBs2+e7W2MdmbmYmdOZ\nOX306NEDLkcaPI1Gg4WFBcYYo0q1a65KlTHGWFhYuGbJpnSzDjrw/xw4HhH3RMRbgYeApw/4mtLQ\n2L70coaZThvnYR7utHdmmLlmyaZ0K24/yJNn5hsR8VHgi8BtwKcy84WDvKY0LJrNZtfSy3OcA2CF\nFTbY4BEeYYaZzng79C9evOgXt7olB97Dz8w/zczJzHxPZj5x0NeThsX4+Djz8/NdY+c4x5HKEc6f\nP8+RypFO2LfNz88b9rplff/SVipZrVajWn2zb99ejXPixImuJZsA1WrVB6/UkwNt6UjaXTvEl5aW\nupZetpdszs7OMj8/b9irZ5GZ/a6hY3p6OldXV/tdhtQXzWZzx3bN9caltoh4LjOndzvOO3zpEH3w\nR3/R2X/6zp/tmrteqBv22i/28CWpEAa+JBXClo50iK5u40iHyTt8SSqEgS9JhTDwJakQBr4kFcLA\nl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJ\nKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIXoK/Ij41Yh4ISL+LiKmr5p7LCIuR8Sl\niPhAb2VKulnNZvOmxjX6er3Dfx74V8BXtg9GxL3AQ8B7gQeBT0TEbT1eS9Ie1Wo1pqamaDQaXeON\nRoOpqSlqtVp/ClNf9RT4mfmtzLy0w9Qp4KnMfD0zvwdcBu7v5VqS9qZWq1Gv11lfX2d2drYT+o1G\ng9nZWdbX16nX64Z+gQ6qh38X8NK2zy+3xiQdoHbYt7VD/8KFC52wbzP0y7Nr4EfEsxHx/A7bqf0o\nICIWImI1IlavXLmyH6eUitRsNllaWuoam2OOzfVNTp48yeb6JnPMdc0vLS3Z0y/I7bsdkJkP3MJ5\nXwGObft8d2tsp/MvAosA09PTeQvXkgSMj4+zvLzcuZOfY47TnOYUp6hTp0qVCSYAOMc5KpUKy8vL\njI+P97dwHZqDauk8DTwUEXdExD3AceBrB3QtSS2Tk5MsLy9TqVRYYYU11phggrOcZYIJ1lhjhZVO\n2E9OTva7ZB2iXpdl/nJEvAz8c+BCRHwRIDNfAD4DfBP4AvCRzPxxr8VK2t3k5CSLi4tssEGdetdc\nnTobbLC4uGjYFygyB6eLMj09naurq/0uQxpq7dU4m+ubnOFMp40DsMYaj/AIRypHvMMfIRHxXGZO\n73acT9pKI2T70ssZZjptnId5uNPemWHmmiWbKoN3+NKIaDabTE1NdS29nGOOFVbYYIMxxphhhnOc\n68xXKhUuXrzoF7dDzjt8qTDj4+PMz893jZ3jHEcqRzh//jxHKke6wh5gfn7esC+IgS+NkFqtRrVa\n7Xxur8Y5ceJEZ/VOW7Va9cGrwhj40pB67Pjf8Ouf/Q6PHf+brvF26F+99HL7kk3Dvkz28KUh9euf\n/Q4/fPCHvO0Lb+OPfuU918w3m80d2zXXG9fwsocvjbhjj93F277wNo49tvNrqq4X6oZ9uXZ9tYKk\nwfSfXvwJ4D3wK/2uRMPCO3xJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqE\ngS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4\nklQIA1+SCtFT4EfExyPi2xFxMSL+JCLGts09FhGXI+JSRHyg91IlSb3o9Q7/GeB9mTkFNIDHACLi\nXuAh4L3Ag8AnIuK2Hq8lSepBT4GfmV/KzDdaH78K3N3aPwU8lZmvZ+b3gMvA/b1cS5LUm/3s4f8m\n8PnW/l3AS9vmXm6NSZL65PbdDoiIZ4F/uMPU45n5udYxjwNvAE/ebAERsQAsALz73e++2b8uSdqj\nXQM/Mx+40XxE/AZwEnh/ZmZr+BXg2LbD7m6N7XT+RWARYHp6Onc6RpLUu15X6TwI/Dbwwczc3Db1\nNPBQRNwREfcAx4Gv9XItSVJvdr3D38XvA3cAz0QEwFcz899l5gsR8Rngm2y1ej6SmT/u8VqSpB70\nFPiZ+U9uMPcE8EQv55ck7R+ftJWkQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCX\npEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiGKCfxm\ns3lT45I0aooI/FqtxtTUFI1Go2u80WgwNTVFrVbrT2GSdIhGPvBrtRr1ep319XVmZ2c7od9oNJid\nnWV9fZ16vW7oSxp5Ix347bBva4f+hQsXOmHfZuhLGnUjG/jNZpOlpaWusTnm2Fzf5OTJk2yubzLH\nXNf80tKSPX1JI2tkA398fJzl5WUqlQqwFfanOc0ZzjDBBGc4w2lOd0K/UqmwvLzM+Ph4P8uWpAMz\nsoEPMDk52Qn9FVZYY40JJjjLWSaYYI01VljphP3k5GS/S5akAzPSgQ9bob+4uMgGG9Spd83VqbPB\nBouLi4a9pJE38oHfaDRYWFhgjDGqVLvmqlQZY4yFhYVrlmxK0qgZ6cDfvvRyhplOG+dhHu60d2aY\nuWbJpiSNosjMftfQMT09naurq/tyrmazydTUVNfSyznmWGGFDTYYY4wZZjjHuc58pVLh4sWLfnEr\naahExHOZOb3bcT3d4UfE70TExYj4ekR8KSIq2+Yei4jLEXEpIj7Qy3Vuxfj4OPPz811j5zjHkcoR\nzp8/z5HKka6wB5ifnzfsJY2sXls6H8/Mqcz8GeA88DGAiLgXeAh4L/Ag8ImIuK3Ha920Wq1Gtfpm\n3769GufEiRNdSzYBqtWqD15JGmm39/KXM/OH2z7+JNDuD50CnsrM14HvRcRl4H7gz3q53q1oh/jS\n0lLX0sv2ks3Z2Vnm5+cNe0kjr+cefkQ8Afxb4P8As5l5JSJ+H/hqZv5R65g/BD6fmZ+90bl66uF/\ncttTsx8+d810s9ncsV1zvXFJGhb71sOPiGcj4vkdtlMAmfl4Zh4DngQ+eguFLkTEakSsXrly5Wb/\n+p5dL9QNe0m3atheu75r4GfmA5n5vh22z1116JPAv27tvwIc2zZ3d2tsp/MvZuZ0Zk4fPXr0Vv4Z\nJOnQDeVr1zPzljfg+Lb9fw98trX/XuAbwB3APcB3gdt2O999992XkjToqtVqsvWdZVYqlbx06VJm\nZl66dCkrlUpnrlqtHko9wGruIbN7XaXzu632zkXgF4DTrV8iLwCfAb4JfAH4SGb+uMdrSVLfDfNr\n10f2wStJ2m+D+kDnoTx4JUklGfbXrhv4knQThvm16wa+JN2kYX3tuoEvSTdpWF+7buBL0k0Y5teu\nu0pHkvbIVTqSVIhhf+26gS9JN2GYX7ve0+uRJalEw/radXv4knSLBuW163vt4XuHL0m7OP+jT3b2\nT9754c7+sL123R6+JBXCwJekQtjSkaRdbG/jDDPv8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1Ih\nDHxJKsRAvUsnIq4Af9mny78DeK1P1+7FsNYNw1u7dR++Ya39sOr+R5l5dLeDBirw+ykiVvfy8qFB\nM6x1w/DWbt2Hb1hrH7S6belIUiEMfEkqhIH/psV+F3CLhrVuGN7arfvwDWvtA1W3PXxJKoR3+JJU\niKIDPyJ+JyIuRsTXI+JLEVHZNvdYRFyOiEsR8YF+1rmTiPh4RHy7Vf+fRMTYtrmBrT0ifjUiXoiI\nv4uI6avmBrZugIh4sFXb5Yh4tN/13EhEfCoifhARz28be3tEPBMRL7b+/Kl+1riTiDgWEcsR8c3W\nz8np1vhA1x4RPxERX4uIb7TqrrfGB6vuzCx2A962bf8/AH/Q2r8X+AZwB3AP8B3gtn7Xe1XtvwDc\n3tr/PeD3hqF24J8CPw2sANPbxge97ttaNf1j4K2tWu/td103qPdfAj8HPL9t7D8Dj7b2H23/zAzS\nBrwL+LnW/j8AGq2fjYGuHQjgztb+W4D/BfyzQau76Dv8zPzhto8/CbS/0DgFPJWZr2fm94DLwP2H\nXd+NZOaXMvON1sevAne39ge69sz8VmZe2mFqoOtmq5bLmfndzPxb4Cm2ah5ImfkV4K+vGj4FfLq1\n/2lg7lCL2oPMfDUz/3dr//8C3wLuYsBrzy0/an18S2tLBqzuogMfICKeiIiXgH8DfKw1fBfw0rbD\nXm6NDarfBD7f2h+22tsGve5Br28v3pmZr7b2/wp4Zz+L2U1ETAA/y9bd8sDXHhG3RcTXgR8Az2Tm\nwNU98oEfEc9GxPM7bKcAMvPxzDwGPAl8tL/Vdtut9tYxjwNvsFX/QNhL3eqv3OoxDOwSvYi4E/jv\nwH+86r/EB7b2zPxxZv4MW/+1fX9EvO+q+b7XPfL/T9vMfGCPhz4J/ClQBV4Bjm2bu7s1dqh2qz0i\nfgM4Cby/9cMEA1D7Tfw7367vde9i0Ovbi+9HxLsy89WIeBdbd6IDJyLewlbYP5mZ/6M1PBS1A2Tm\nRkQsAw8yYHWP/B3+jUTE8W0fTwHfbu0/DTwUEXdExD3AceBrh13fjUTEg8BvAx/MzM1tUwNf+3UM\net1/DhyPiHsi4q3AQ2zVPEyeBj7U2v8Q8Lk+1rKjiAjgD4FvZeZ/2TY10LVHxNH2SrmI+PvAz7OV\nJ4NVd7+/3e7nxtZdxPPAReB/Andtm3ucrVUZl4Bf7HetO9R+ma2e8tdb2x8MQ+3AL7PV/34d+D7w\nxWGou1XfL7G1auQ7wOP9rmeXWv8YeBX4f61/378FjANfBl4EngXe3u86d6j7X7DV9ri47Wf7lwa9\ndmAK+ItW3c8DH2uND1TdPmkrSYUouqUjSSUx8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJ\nKsT/B/NN0ugCbls8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f15796c1e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_data(centroids+2, X, n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## course.fast.ai"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "If you found this interesting, you might enjoy the 30+ hours of deep learning lessons at [course.fast.ai](course.fast.ai). There's also a very active forum of deep learning practitioners and learners at [forums.fast.ai](forums.fast.ai). Hope to see you there! :)"
   ]
  }
 ],
 "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
}