{
 "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": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import operator\n",
    "import torch\n",
    "\n",
    "from fastai.core import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_clusters=6\n",
    "n_samples =250"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To generate our data, we're going to pick 6 random points, which we'll call centroids, and for each point we're going to generate 250 random points about it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "centroids = 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": {},
   "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+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXt4VdWZ/z/r3HIBQkISQhLAQLkEhVJqQNSiolERKaDTGS9tR6mVqXfrb3Sqtk2odpzqlLZqL4PTUmfaCrajlAqtigWlImqQoqBBIvcQIAmXkOu5rd8f6+xz9jk5uUDunvfzPOc5++y9z15vQviud73rXe9SWmsEQRCETz+OvjZAEARB6B1E8AVBEBIEEXxBEIQEQQRfEAQhQRDBFwRBSBBE8AVBEBIEEXxBEIQEQQRfEAQhQRDBFwRBSBBcfW2AnaysLF1QUNDXZgiCIAwotmzZUqO1zu7ovn4l+AUFBZSVlfW1GYIgCAMKpdS+ztwnIR1BEIQEQQRfEAQhQRDBFwRBSBBE8AVBEBIEEXxBEIQEQQRfEAQhQRDBFwRBSBBE8AWhN6mpgSeeMO+C0MuI4AtCb7J8OTzwgHkXhF6mX620FYRPPYsWRb93RE2N6RwWLYKsrJ6zS0gIxMMXhN4kKwvuv7/z4i0jAqEbEQ9fELqT7vbIT3dEIAjtIB6+IHQn3e2RxxsRyMSvcIaIhy8I3cmZeuSnMzKwOhUwnYEgdBIRfEHoTiyP/HRpT8RjOwMJ8whnSJcFXymVDLwBJIWe9wetdYlSagywAsgEtgBf1Vp7u9qeIPQ7ziRuX1MDTz9tju+8s30Rt3cGixZJ1o5wxnSHh98CXKq1rldKuYG/KaX+DNwH/EhrvUIp9QvgFuDn3dCeIPQvziTEsnw5LFlijt98E557ru3v2juD0xkJCEIMXZ601Yb60Ed36KWBS4E/hM4/CyzsaluC0OfEmzBdtAgef7xzIZaaGigthepquOgic27duviTvFZbEJm4ba8tSeEUOqBbYvhKKScmbDMO+CnwCXBCa+0P3XIQyO+OtgShT4n1sC2vev78jr3rmhq46SZYu9Z8LimBKVPg5Zdh1qy226quhh07YOlSmDixcyMBQYiH1rrbXkA6sB74AlBhOz8K2N7GdxYDZUDZ6NGjtSD0a6qrtX78cfOutTkGrefONe+PPx7/vurqyD3FxVqXlESfmzu3dRtvvWXOX3RR63vitSEkLECZ7oRGd2uWjtb6hFJqPXA+kK6Ucmnj5Y8EKtv4zjJgGUBRUZHuTnsEoduxwiqWN2950/Pnw/Tp0NAQ8frtI4Hly41nP3cuPPuseU5NDZxzDni9xnu3vtfQYOL7c+ea79x/PwwebO6xI+mZwmnSHVk62YAvJPYpwOXADzCe/pcwmTo3AX/saluC0C+IFVpLbAcNMucHDYruCJ54woRuxo2DMWOgttZk6Lz5ponfl5TA6tURob//fiP23/kOXHJJ22GiRYvMd6xOxupEZOJWaIPu8PBzgWdDcXwH8LzW+iWl1IfACqXUo8BW4Jfd0JYg9D3xhDb2PBjhfuIJ0wmMGwcVFeb18svmHUwHUFMTEfqSkkhHcMkl8T13u6jbOxlrJCFev9AGXRZ8rfX7wLQ453cDM7r6fEHod2RltRZai3ffNWEY67zl6e/bBxUV1Kank2mJPcCePfDBB9QCmamp5nvr1kFhIYweDZMmGRGfOTPyndi8/PbeBcGGrLQVhNPBnpUDRlitc/v2GbG/+OLWgvuVr1C6YgXP1NayvqiICW43HDoE+/bx8ahRzE5K4taaGkrvvBM2bDDPuesuk6Hz1a/C4sXxV9rGruw905W+QmLQmZnd3nqde+653Tx3LQjdjD0rJzZTJzMzkoVTXm7uuf9+rUGXZGZqzPoUnTdkiN4JWt9xh945dqzOS00NXyu54w6TwVNSovWKFVpnZ2t9yy2ts3ssrHbKy/vglyH0F+iLLB1B+NSzaFHEA7fi6A0NUFxsQjHZ2fDII6Zcwrp1UF9P6ZgxLNmzJ/yIQ6dOMXvIEJa98AKLq6o4ZHv8kp/+FIDS4mJzoroaRo6MZOysWxeZFF6+HF55xZwDWLOmd34HwoBFBF8QToesLJNWaYm9VSKhpAQ8HiPKGzfCtGmwbh21p07xzP79UY9YyEI2nNrAvFNVpJPOQi5hFavC158B7lq3jkyv13Qkc+aYC+ecA6mp0SUW7r/ftBubsikIcRDBF4TTob3KlbW15nj+fMjMhL/9jcy33mI9MBs4hBH7e7iHBSxgCUsooYQCCgBYxSryhg9n/Re/SObrr8Mbb5jn+Xzw+uumUyktjbRnvUv6pdBJZAMUQTidDUVi69XYNyhZvdp4+I8/brJstm8HYEJSEuvz8shLTWUDG9jLXgooYDnLKaCAvexlAxvIA9aPHs2E1183aZuZmaYNHbMeUXLthTNEPHxBsER8w4bIKti2sC+oKi2FxkbzOTUVbrjBHP/gBxFvPyMDhgxhwv79LMvNZV5jFUtYwnIiBc6WsIQTnOA3ublMKCszJ7OzTfy+uNiEh2bPNvMC7dkrHYHQEZ2Z2e2tl2TpCH2CvaaNVQunI6zMHPvL+u6KFVo7neZcqA7OzpQUnZecrNNJ18tZrtezPvxaznKdTrrOGzxY7/zyl7UeN07rP//ZZOQUF7fOCoq116qpU1Jyej+D8KkBydIRhE4SOxHbFvZNS+bMMd73xInGuwfjkVsx9kDAvE+ZwscffMDs48dDMfw54TCOPYZ/CZewqn4Vs//0J9bX1TFh82aziGvdOrMAa+1a8+ynn25t79NPR1bqdrZMs5CYdKZX6K2XePhCv8bu1ceOCOzXSkq0vuMOrceO1TVFRTovlGNvvRayUKeTrgGdTrpeyMKo63lDhuiaO+80z7r4Yq3T083x6NGmnfLy6CqZlmdfXCyVMxMUxMMXhG7GyrkHE6+3CptZ16qrYetWc+2++2D3bjJ37+bWmDz8VawiD/jNyJEsPnaMVY2ropq59bzzyLQmbCsr4cQJc5yfH4ndWzX177/fxPatkg5PPx3J05c4vhCDCL4gdJasrEjIBiIboDzxhBHY7GwTglm92uTFe70wcSKlO3fCnj2ENjQ02Tjnn8+E1atZP38+s996K7z4qmTcOEqfe858sIqojRkDX/oS3HKLefb8+dGdjT3E09AgxdOENhHBF4Su0F4hswsvNCth33qL0uRkaG7mGaVYrzUTvvAFWL6cCUuXsv6++5i9bRu3ZmVRes015rtZWeb769aZAmtgxD5e/Rw7N9wQXZ5ZEOx0Ju7TWy+J4QsDjrZq2Vhx9ZhXjRVrt7JvQnMANfb77btmWffF3N8Kaw5BMnQSEiSGLwg9TE2NidWvXdt27fqRI6G+PhyHzywuNp67taNVyBPPtHLsIdo7nxaqPP7II3DFFdHX7Hn3UhZZ6AQi+IJwpti3LbSE1hLhG26ITKSCSeG88EJz/rnnTJmEO++MTKxa8wP21E8w8wMWzz0XvauVtUMWRO+8ZUcWYwk2RPAF4UyJV8/GHtN/9tmIeFse/E03RcovxBNgqxgbGAG3dspat85cs+9qVVxsOo72vHrZAUuwIYIvCGdKvMnT2E7AyuqpqYmIvX1EEIs99ROM2FvlFawtFe0lmq+4onOlICTUIyCCLwg9iz38Yol9e/V6srLMaMDaVcte+96+pWJnVgZbzxPPXggh1TIF4XSw8u537oxfYTO2mqY9pPL4422Lvb1ip/Wd1asjlTgXLYoum2Cv0ikInUQ8fEE4HeyVKu2rXS3sIZSaGuPZx07Qtvfc2GdYtOWpy6SscBp0WfCVUqOA/wFyMLVAlmmtf6KUGgasBAqAvcA/aa2Pd7U9QehT7OWR7atdLezC/MQTZgK2rQnaeM/taGFVLDIpK5wGyuTsd+EBSuUCuVrr95RSQ4AtwELgZuCY1vo/lFLfAjK01v/W3rOKiop0mVUPXBAGOp3xvrvqoYuHLwBKqS1a66KO7utyDF9rXaW1fi90fAr4CMgHFgDPhm57FtMJCELi0Jk4e2zMvyfaEIQQ3RrDV0oVANOAt4EcrXVV6NJhTMhHEAQ7kjYp9CLdlqWjlBoM/B9wr9a6zn4tVOshbuxIKbVYKVWmlCqrrq7uLnMEoe84nT1yxUMXepFuEXyllBsj9r/VWr8QOn0kFN+34vxH431Xa71Ma12ktS7Kzs7uDnMEoW/paphGEHqI7sjSUcAvgY+01kttl1YDNwH/EXr/Y1fbEoQBgYRphH5Kd2TpfAHYCHwABEOnH8LE8Z8HRgP7MGmZx9p7lmTpCIIgnD6dzdLpsoevtf4boNq4fFlXny8IgiB0D1JaQRAEIUEQwRcEQUgQRPAFQRASBBF8QRCEBEEEXxAEIUEQwRcEQUgQRPAFQRASBBF8QRCEBEEEXxAEIUEQwRcEQUgQRPAFQRASBBF8QRCEBEEEXxAEIUEQwRcEQUgQRPAFQRASBBF8QRCEBEEEXxAEIUEQwRcEQUgQRPAFQRASBBF8QRCEBEEEXxAEIUHoFsFXSv1KKXVUKbXddm6YUupVpdSu0HtGd7QlCIIgnBnd5eH/GpgTc+5bwGta6/HAa6HPgiAIQh/RLYKvtX4DOBZzegHwbOj4WWBhd7QlCIIgnBk9GcPP0VpXhY4PAzk92JYgCILQAb0yaau11oCOd00ptVgpVaaUKquuru4NcwRBEBKSnhT8I0qpXIDQ+9F4N2mtl2mti7TWRdnZ2T1ojiAIQmLTk4K/GrgpdHwT8McebEsQBEHogO5Ky3wOeAuYqJQ6qJS6BfgP4HKl1C6gOPRZEARB6CNc3fEQrfUNbVy6rDueLwiCIHQdWWkrCIKQIIjgC4IgJAgi+IIgCAmCCL4gCEKCIIIvCIKQIIjgC4IgJAgi+IIgCAmCCL4gCEKCIIIvCIKQIIjgC4IgJAgi+IIgCAmCCL4gCEKCIIIvCIKQIIjgC4IgJAgi+IIgCAmCCL4gCEKCIIIvCIKQIIjgC4IgJAgi+IIgCAmCCL4gCEKCIIIvCIKQIPS44Cul5iildiqlKpRS3+rp9gRBEIT49KjgK6WcwE+Bq4CzgRuUUmf3ZJuCIAhCfHraw58BVGitd2utvcAKYEEPtykIgiDEoacFPx84YPt8MHROEARB6GX6fNJWKbVYKVWmlCqrrq7ua3MEQRA+tfS04FcCo2yfR4bOhdFaL9NaF2mti7Kzs3vYHEEQhMSlpwX/XWC8UmqMUsoDXA+s7uE2BUEQhDi4evLhWmu/UupO4GXACfxKa72jJ9sUBEEQ4tOjgg+gtV4LrO3pdgRBEIT26fNJW0EQBKF3EMEXBEFIEETwBUEQEgQRfEEQhARBBF8QBCFBEMEXBEFIEETwBUEQEgQRfEEQhARBBF8QBCFBEMEXBEFIEETwBUEQEgQRfEEQhARBBF8QBCFBEMEXBEFIEETwBUEQEgQRfEEQhARBBF8QBCFBEMEXBEFIEETwBUFIaOq0nxe8R6jT/k9le3ZE8AVBSGjW+Wr5tfcQ63y1nf5OV0T7TNrrLnp8E3NBEIT+TLE7M+q9M1iiDXCtJ6fH2+suuuThK6X+USm1QykVVEoVxVx7UClVoZTaqZS6smtmCoIg9AxpysW1nhzSVOf932J3Jjd78toU7fZGAGfSXnfR1Ra3A9cC/2U/qZQ6G7geOAfIA9YppSZorQNdbE8QBKHPsUS7LboyAuhJuuTha60/0lrvjHNpAbBCa92itd4DVAAzutKWIAjC6dLZWLv9vjrt53cth/hdS1WnY/Sx7XQ0AugrempMkQ9stn0+GDonCILQa7TnaddpP+t8tRS7M3nJW80K32GadRDQrPAdASBZObjWkxN1b2wopk77+XHzPsoCdeF2OhoB9BUdCr5Sah0wIs6lh7XWf+yqAUqpxcBigNGjR3f1cYIgCGHamyC1OoNmHeCjQAMAH/hPEVAagHMcg8Lfs+7dHqjn3uSzokR/na+WskAdRc40it2Z1Gk/L3mrAc08z/DwPfE6i96mw9a11sVn8NxKYJTt88jQuXjPXwYsAygqKtJn0JYgCEJc4nnalrd+nmtoSOwb2RasZ6RKYodugJAKTXCkhoW62J3JVn8dZYE6XvIe5cakvFbPAQUQHi0AVASbGOdICY8Y+trr76nuZjXwO6XUUsyk7XjgnR5qSxAEodPYwzzJysm24CmG4uTrSfmUBxo4pQMcDDaDUlHhoLHOVLYF62kJ9Qixz/m19xDJyoHVY+QqD2WBOsY5UrnenUOzDoZj/H3l8XepNaXUNcBTQDawRin1d631lVrrHUqp54EPAT9wh2ToCILQHyh2Z9KsgzTrABe7h/Fnbw1H8PJc82HO96TTojXbgvWMdaTyD6eSOWf4IL7X9AnD8QDwsb+B/6rdwbwRnwk/D6BZB6gOetkbaGKuM4skh4MkFPM82eHOwYSOzPOh9z3+rmbpvKi1Hqm1TtJa52itr7Rd+77W+jNa64la6z933VRBEITO01aGTppy0UKQFb4j/KR5H0GCABzDFxJlI8arHv0hDxUV85PtmygL1LFPN1HkTOPtjz/k3qLZfOU7D0R56RXBJtb4a9ihG9garONF31FAkaZcFLszKXKmsS14im3BeqY6Boc9/t4stSClFWx4dROfeLfi1U19bYogCF2kvRIGuwONAJQHG6nGz0iVxHTnUACqtZed//4LVn//Rxw6dIhnr/oq9bv2McExiHn7W9h69Tdorqpmw78/xaLv/Fu4rbJAHYWOVHLxMF6lArDVV8fvWkzY597ks7jePYLr3TlMcg5ihe8w63y1vVpqQUor2DjgK+cjr8km/YxnWh9bIwhCVzjPNZSt/jpOBo0XbXnilcFm/Ggmkso4VyoeFEnKycXuDI62ePntI4/z8WPLws9prDrKO1f/Czf/1y8o/sbtnKg6Er72wveXUuoawq3f/RbbA/WMUsm8GDxKXdB46ztpZKevkYpgE/cmn8WNSbmAGX0kK2dU9lBv5OyLh29jlLuQSZ6ZjHIX9rUpgiCcJlZopDLYzAveI7zuO862YD0v+o9Gec+/aD7IjmADJ5SPIcpJknKywneYt/0nubY+mQO/fjHquQtZiLvKy03zr0UfamYhC6OuP/XMf/HS4QrKAnUkKUWu8tAQChOd4xjEVMdgygJ1UTbYyyv0ZqkFEXwbHpXCZzzT8KiUvjZFEITTxAqN/LKlMpQ9o8MhFLv3PNaRDMAR7WOF7wiv+48x15lFsw7wfpqDmWuWkZKbDRixv4d7+BE/ooACfsSPuId7wqKfmjucqS/9nJb0IeQqD7VBH2naCPc5jkHcmTyaSc7BrWyAvimTLCGdHsSrmzjgK2eUu1A6EUHoYSxBPc81lMn+wVETqpa4Frsz+YekESQpBy1oNvmOU6W9HAvU0hLQTHUM5ryJZ8OaZbx39TfYULWBBSyggAKWsxyAvexlAxsYnDucf/3L82wZm8YO3UCV9lIVOAbAVMdg7k8ZwzpfLSt8h7nZk9fKg++Lejvi4fcg1pzAAV95X5siCJ96rNBIviOZYncm63y1UXnv1sRomnJxY1Iei5LyucCVDkALmiE42Basx4Vi8Piz+NyTD3OCEyxhSVQ7S1jCCU5w989+yDenXMjNnjzuSx7NVMdgJoYmayc5B4ezc9qqqXOeayhFzjTOcw3t4d9MBBH8HsQ+JyAZQILQe8RmvsQKb2Wwme81fcIF7vSwSH/ekcZUxxDcQU39rn28ffcjpJNOCSVRzy6hhHTSeeaOB/j2Bxs4zzWUXEcy01xpfD05nyJnGhe7M9qtvwPwuu8YZYE6Xvcd6+HfRgQJ6XQz9jCOHcvbrw0c4nPJl0qIRxB6kPZq6NRpP99v2s1B3UJTc4A6zCjgBAG2BU9Rv2sfb129mOaqauawkAIK2MtelrCEEkoooIBLuIRVh1bxzJwbObJ2OQsmTWOF7whTHYPZFqzH1xwMpV62V1JBxbz3PCL43Yw9tRMIH49yF1IbOMTRwH4O+Mqj0j4l1i8I3YsV3rFi9806EBZfgIO6haE4CaCp0l5GqiS+kTySlw5/wgOhPHuAVawCYAMbOMEJvsk3jdiHzjdXVbN67iJGvr2G64fnU6cDbAvWsy1YzyTn4FbhHLvXP8+TTbJy9GoJZQnptIM9DNNWSCb2vD2MYz/2qBQ+l3xp3LRPifULQvdjlS02E6OK690jaNaBcOz8JAEmOQdzjmMQQ0O+7w25E5lzy1einrOKVTTnupnx+58QyE0Ki73FuJuvZVeGi48CDSSFJHWqYzDzPNmt0i2tUNOPm/cB9PrOVwkv+O3F1u1C3JYox563p3bGpnm2lfaZ4cxlkEonw5krsX5BiENtk+ZnW73UNnW+oK69bPHF7gwqgo2s8B3hbf9J7k0+i5s9eVzhyeSk9rNDN/B0836+1fgxTf/2ZT7/4G3h56TmDue7L/8f93zxS6z566sMyY2EZ2Y9dCdjHvoXALYF60lSKrSSdnBcm6wSC/a8/N5Mz0z4kE57q2tzXAXUBg6R4yrAo5Lxax9+7cerm8KibXnrXVmsVeHdQoM+QYV3C5nOPFntKwgxrCz38ehmLwC3T/N06jv2OL4l/lYNGzDe9e9aqjioW8jBTU3QyxF8jFRJlH7vMW4iSMWvX+C8Nf/FxIkTTIho4nBKXvk9j175j1y06EYWfvv/mcYUoUJpw8NefEWwsVXt/DTl4pakfGghnJ3Tm+mZCS/47Qn2Ef9ejgb2k+nPY5S7kJPBao4G9uNSLka5C8Nx99MV5tiY/dlJF0ALnJ10AR6V3KY9gpCoXFfojnrvDPZa+Oe5hprSB45kVvgOh3eyskoZj3AkhWvi35M8mpXeI4x56F+47BuLGJKVGSXOrxcM5l/efokP05286D/K9e4RtOggHwUbuFj7Oc81NNzBWPF6e7bO2/6TlAXqmOwfzLWe5HYnmLubhBd8K8wSD3tncMBXztHAftwk0Rg8RUXLVnb7t9ESbCLJkUKOq4Aj/r1hEW9vIjZ2VDHYkcGMlKvD10+nA5EJXyERyExR3D7NEw7tXFHgYlWFDzQsmuIhMyU608XsOmWqVc7zZIdFdpwjJWoi9WL3MD70NZG8Zxj/VJDG/FSz3aEZDQxhUs6gSNmFkDg36wAt2ZqJWpOknIDmRf9RAJ5u3s9J7eegbgnvgBXrwccKfG9uh5iQgt+WSNYHj7O9eSNpzizGhWLtlviOchey31dOgz7BPv8OUlUaAHXBGmr8lezzfkgjdfi1j4lJM9pNw4wdVXRFtKXgm5BIWKGdFeU+Kk4Y7zzVrVqFecwK18i+tJaHf7F7GPmh0goAb/tP8j51HKjxcKs3h7RpLiyvf5IzlXkeU2KhWQeoDDbztv8kAC/6jnKzJy+cCdSiNbuDjfi15qBuYaRKCodz+lLgY0lIwY8nkl7dxLtNf6FBn6AmWIkTdzh0Y31nvKeIHS0bSVVpnNTVDHeOZpznXJqa19OgTwAQCOX05rgK2B8aFVhpmHZht0YNOa4CPmzZxNHA/ih7Okt3zCEIwkDhukI3bx0K8Nr+ALNGOpie44wb5rE8cVBRMXwrjAJmFNCsAyxQOZxKT+e6CeY58zzDw5Usrfj7Ct+R0FaIp7jePaJVuuVQh4v7k8ZwSvv5ZUsl13lywmGc/kRCCn48kTwQ8t5TSWOEewwB7WOXrwy/9uNSLj7ybmaQSsdHC27lYZJ7Zli0q2oOMiIzj0bqcOLCq5v4sGVT+HxLKK3T3tH4tZ9dvjKq/QepCR5kuHP0GYl2eyEpQfi0kZmi+PGlyaws93FFgYtX9kYyW2JXtt6YlBe+Fi9Obo0CbvbkccvU1PD51h648fj9BLnencM8T3bcVEswIZvvpnyGF7xHwueAXq+Z0xYJKfjxRHKUuxC/9gOaMZ4p7PF+ELqiyXEVUO0/gFuloINBJiTNYJhzBF7dxI++9zN+v/yPvPjqSvRoTbX/IAH8bC1/l4evfpIrbj6fGx+aR33wuJmcDbVV0bIVgGQGMdw5OjRha8I5sat1JUYvCBGseP7PtnqjMnfay3bxNTs5XJ6Br9AJof9GVpinvVo2lcFmPvDXk6Pc7Ag2MN01lDTliupc4nUm7Z3rSz41gt/VyUuPSgl78gbFWNdUAvjZ3ryRmmBl+N7jgSqGOUfwr9+9l59+/78BWFj8JR5Zcxf543P45MMKvjvvp9RWHee5x9aSqtJY+OBFBJuDZDjNH6Iz5CE0U09NoJJgc4AMZy5jPJPZ493eanQBEqMXBDuxmTuf9w/jreoAn88ZBjGZm/HSOmOzZeLxy5ZKdugGgPAkLLT26mM7mNhRQl979hafGsHvjslLy6O2wi3DHLkcC1YBkEoa2a5ReFQyo9yFlJaW8tSjkV1xqquO8fDVT3LHk9fz9N2/41jVyfC1X/77CgAWPngRNcGDuJSbMZ4puJSbHFdBuEOpCVbiUpFJI9ASoxeENrA8fTALs76z3s9r+9M4PNrPjy91kZmiqG3S4fAPRKd1xnrhdq/d1+xkZbmPayfk4XMEGetM5R9sq2J7M5WyO/nUrLQ93d2q4q1otUI9+e5xDHeOxh/apsyBk0bqaNKnGOOZwqljjTzzzLKo5y1kIYEqzff+8ecEq1SrXXH+tHwddbX1ZDnyyXEVRIVsGvUpALIcIxnlLmSMZwqTPDMZ45nSqU1ZZHWukOisLPfx2v4Amcnw2v4A9/61OSz2j272smqXr9V37DtN1TZpHtl7NFxh0/reGx87eSR1PIuS8lstoLK89t7exKQrfGoE/3R3q4otiWAXTWvBVSPGSw8SYJBK52hgP+81vUrN4N0885elDMs1sb+OdsUZljuU5S8/xbjhZzM5eRaVvgo+8m5mj3c7H7ZsolHXMUilMzn5C2F7OvOzWDbv8X4gtXiETy2xZRXilVmYkeskMxlqmwmLvuXZXzbaSZMfHt3s5ZpVjVQcD7Z6xspyHyvXD2JitcmTt75njQzaojc3IO8OuhTSUUo9AXwR8AKfAIu0NvmJSqkHgVuAAHC31vrlLtrarcSGSuwhoRxXARXerfhoIYUhjHRPJN89Lir0MnRsNv++5l4euvrH7e6KMyx3KP/Sj6TaAAAgAElEQVS55iGSxwbDq3atkM3xQBUTkmYQbAkwyJHO+00bOKYP49d+JiZN7/BnsGwe7y6SvXiFfo3lbV9X6G61SKojYuPv9s9XFLh4eGMz++qCYbGvbYZx6cosztplPP8JGYpx6YqKE5olm1o4P8/Jo5u9NPo0qW4VEvYUrhuRRppS/Gavl9f2Bzg/z8+4jLZLOQy00E5XPfxXgcla688CHwMPAiilzgauB84B5gA/U0o5u9hWt2IfEXh1E37tY7y7iFHuQip9FfhoASCJVAL4qPTtIs2RFf6+VzeRPz6H//fUre3uivPNJ7/GqPEjaNQnGaTSwyGbLEc+NcFKPm55lzRnFvv8OzimD4e+3bkCUVYYa4xnsuzFK/RrLJFeWd46tGLRVoE0u7dd26Rp9GvuO9fNFQUuvvaXJjZWBtl/yoj8Ty5NDgv7K3v94VLzKW7Fr+akcNloJ/ec66HRr7ltqotNlf5wyOf2aZEVu9cVuvn2TE/cHH+7nb25AXl30CXB11q/onU4eLUZGBk6XgCs0Fq3aK33ABXAjK601ZMc8JWzy7cF0BzwlRPQ5o8yVaVxgiPs9m0z1xWMdU8ly5HPtJTLafwEfnjXM+3uivOju3/FJ7t2k0oaGY4c3mt6Fa9uJsM5AoCa4EHqAjXh72U5RpLvHt+pmHyibrreWANvPmHehYFBewJqEa9TqG0yHvlr+wO8stfP8u1elpb5QMEre/1UnNCMSVNcP9FJVori6a0tVJzQXDbaLMhaNNnDt2d6WDjOzSt7/fz40mTeqQqwtMzHx8c1mw+HOhcVLeTWhHC80UhnOq/+SnfG8L8G/Dl0nA8csF07GDrXCqXUYqVUmVKqrLq6uhvN6RgrBp7jKmCSZyag+Mi7GadyM8kzkxkpcxnvPpdhjlwAnLg4O+kCZqbO51DFUe646tscqzrJJVwS3hVnEYvYy97wrjjHqk7y3at/ytFPTnAwsJOaYCUftmxijGcK493nMtY9lUEqnXSVwzBHLpOTv8AR/16JybfD1uWw7gHzLgwM2hNQIMpzt3cK1mSsJeDW4LfJp2n0ae4rcrP62lQq6zWbq4JsrtLMGungx5cmk5miwu2uqjAi/bOt3rB3PyHDwc3nOJk10sHCce5Oj0Li2TlQ6HAcopRaB4yIc+lhrfUfQ/c8DPiB356uAVrrZcAygKKios4Xu+4GrBh4tf8gaQ4TgxvrnoqVDulRKUxMmoFXN4UWYim8uona2lpmz55NTZXZi7KjXXGOVtVwx1UP89RbDzM8c0R4kdXEpBnsbHmH3YFtoV8G4QJsIKmYbTFtUfS7MPBZWe5jaZmPb8+MdAq1TRFRXzTZnF843s27RwJsPRJg82Edvn9yloONlUGC9bVMzxkR1bHUNmnWf3QEGMaaT3zsr4fLRjt5bb+fy0Y72XgwyCt7/Z2qyGnZeV+R+4znJPqSDj18rXWx1npynJcl9jcD84Ava60twa4ERtkeMzJ0rl8xyl3IcOdoaoIH2e3fxm7/NuqDx9nl2xLlXZuQiWKXr4yKlq0cGLSNy2+eGfWsVawiNS+FH/5fCY5c3WpXnCtvvpC0zMEMdgxlsCMDMCOM4wFT4GmYI5fx7nPDHU0ihmo6S2oWXHi/eRcGEE11sPVF8x5DvJDPynIfS7f4SHUpjjdrvrqmiec+8rLxYJDNhzVj0ghn0dwwyUPGG4/R+Pj5zEreG35GbZPmq8++z1/unsmpl/6d/fUwJg1GDVHcNtXNhAzFbVNdNPqNdLU3CrHbiWZAhnW6mqUzB3gAuFhr3Wi7tBr4nVJqKZAHjAfe6UpbPYG17WCFdyvH/IdxoDjLPRkwmTrRmD+IumANNcFKbnnoBpIZzPLHngcgMzeDZ15eSrDgOD9a+12+Ofd74cVXNzw4lxsfmkcKQ5icPCv8xAO+cmqCB8lyjCTDmRPOu4+HlEEWBjzlr8HmZ2nwa551zYvyju2LqCDau7+u0M2tLzexuSrIyRYHs0Y62HgwyJ46+MpLjQwf5KB85aPsev4xAK696lJe+PNf2RIYw86dH/Pyt64ieLKK+rWPMdQDe654iD07/CEvPxD29lNdratuxmIv05zqVgMurNPVqeWngSTgVaUUwGat9Te01juUUs8DH2JCPXdorQNdbOuM6EgoPSqFJJXCCW08bZdvO0cD+xnqy8al3OHvRa+M/Rs1wYNc+9AlNFPPK7/exPfX3M2wsalkOD9DoNDH/77yC2664g4uv3kmNz40jzSyyHEXhDc4AfvKXh+7fFtwKXebq4SlDLJwxjTVGbEtvAxS0vrOjsLLAFjpm9Xm7lVW+majT7N0iwmdLP/Ay4E6s0uVS2l+VpzK8u1e/lBuwjM7VjxK/drHws84dOgQVxRfivNLP6F55d0ET1aFr1WueozCJMWNd3+HGyZ5OD/Pz4xcJ+BlRq6T/3ynhSa/JsWtwmGkeMR2UAMFFYnC9D1FRUW6rKysW5/5iXcrH3k3M8kzs02hNDH67YAm3z2eI/69tOgmdvu2MdY9lSSVwih3IV7dzIctm0hhCPsCOwDIcuQz9NRoDg/+iKnJs6n272eXbwtZjnx2V+9kRFYejbqOLMdIaoIH49oRr1OKPdefPfzGGjOBOm2RhFn6JVtfhM3PwsybjOjaxb+XOgN7Hj7QKv4dK/SzRjqYnOlge02QjZVG7DOTYflVKRSNMBnej2xq4advHubYY+fjOx4R9YUsDM+npZMeNZ8G4Biay3d/X0Z2VhbXFUYma2flO8JtAXx7pqfNDqm/xe6VUlu01kUd3Tcwkke7QGcmQM0EamSh02BPBjtbTASqLlATLpxmlTK24u2gGOOZzAFXOQ3eExwPVGEl/qY5szh3xAi8wWZOBY8xyDGUDGdO3E1P4hHr0ffnMshW1gyY2LrQzwh51hTMgL/+BPZvMZ+nXRMOs4Q/Q490ArGLp2KF1Lp+X5E7HLIB2FgZZOYIxZFG2FOneacqwJihDlaW+7hqrIvtNdmkP/UyK745h7rqQ+FV7wtYwBKWUEIJBRQAZp7NMTSXMQ+soSUpg0c3e3nrUICSC5IAqDhuMsxHpMI1402e/8+2eqM6KatDsn6OgcanXvDPVCjtIRwrc6Yl2ERN8CDpjuFMTDLLCmIXbRk0Xt3MYf8eGjETVMf8VUzyzAx753ZBB1qFazqbqdMfPP+uZs3ICKGHSUkzYr71RSP2o8+NdAKx7xC/E4gltlPooJPoKAPm+oImPnfoVQrHXQ46mY0Hg0zOdDB9hJNNlQH21AWZle/gigIX9/61mdf2B8gfBJUNAGNIuv0lPE9ezYbjba96d2fkknffGhozxvPX/QFm5TvCq2lvn+bhuj8ZwT/caEI2r+z1hzspINwhdbSeoD/zqRf80yFWPC3x9biTOeArZ7RnEkmOlCgRNqWMtzDeXRQWXJdyhxZymSqbI1xjcIbmA6w2rElh+7Psx53tqAZybN8Sel8DvB5aqCwjhB7ELu6WKFudQex9/mbwNcPxStj7TmshtzqFQ9vh0ns67CQyOcXtvAZcBrTuEIbt/SsX7P8N5DlZNGVheEJ0ZbmPzVXG258+wskre/28tj8QXk1rMSRvPP7rn+T4z/+RJSwJiz1EVr1/+aHl/HXQeAD2n4Kg1pw7HA6eCvDIphY+M9SBL6CZNjx6F63Y4/4UyjldEkbwO+MJtyWe7Yuqjnk3wt2im6gL1DA5eVZUGubfm/8adzvDrpZ07qmc/c54350N6VjPmjgfdq6OCP1FJVD8eOsRgnj+3UCs5z3tmkh6ZMEMqHjD/OlOuTriqW9fA1XlULkNqivMqMDfDK7kyHMKLzNiv39L5PkQPVKw09GowfZ9+4TodYVuapuCbK8JcsloFxsO+LltqpGtBp+fqgZI88CDY/dz24N3t7nq/Zt8kxceu53Bd60hc9R4TnrhYL3mYD1sORrAlPxqHbe3H1uji5ILkhiXMTDrTiaM4HfGE25LPNsTVSv0E+udW7tbxdpwNLD/jLczjEdPx/Y7I+YdhXRqdsLL90H2OfDWE7B3A1SsNUJ/cYnRm3iiHtu2dABngF1orQlbfzOUrYwINpipJ1dy5BpAej5MvhoCfjjwdziyE5pOgjvJ/KNd8DXImxwt8ierzDzBBV+DjNDi+qY6M1oouq7tDsHqRGLCQma1rIONlX48W7xR3v1tU11srwmSfqqCu6+7iqZjVVzJwvCqd3sM/xIuYdWxVfDU1Tzxh3X88cRZNAXA7TA1ePbWwbThjnZDNVaJB2jhf6/uX4kTnSVhBL+zk7fxxLM9UbWuWWUa2htB2G3ob5k2bdGZ+Ly1EKotXr7PCLy3HsbNhYu+AwWXmGe+/TS8EQrnzC6NFvVpi8DbYEYC1nmZHD5N7J63Jf7nXmcydgpmQMYoqPwAPlgDLfWQPxWmLoSaPcbD3/RLOGFbM3l0F1SZDDU++RvMeciIs5UJlJwGzXWmkxj1uUi7W1aaNtubBG5jFGCJsFlk1RLOnb99WhJP/62K7/3znHDqZUer3puOVfH1a4vJfOgtLhifzeaqICluxeaqIMVnuVplDdlDOGZytyU8yTsQSRjB72lPuKMRRFshpdM931vYhber4nrlUvDVw8mDcGK3EXvLYz/4prnH1wjrS83n3evMuQvvB8+giMiDGRVMWyTefoc01RkRVxgv3fKgwQj93ncgeQikDIWaT8x55TAi73TBrMWw6VdmBDDiHNB+cLhg2Ggj+ElDTEewcZkR9tyzzYjgRCWk5ULQHz2ysL+3we7c2byd7uO83NmMtZ23h3jsopuZohiUnknqhTdH5eFb2TgZNz5D08q7jWdvI/mCm3EMzoRQSro1OVzbFOQ/32lh0RRPq6wigHEZjjY9+/6arhnLwAxE9UM62nErdsOVMz3fW5xpgbLGGnj1Afify00oByBrIjg8RuzTxxqv3RLs3euM1w/G09+9DsYWR+6ZON9c9zaaeL9nkBF4y74Xb5KqmXGxvOqyleYYIjH8ve8YMbbCJzkTzXUdNB76/i3mnkvvobbwWsgZb8I5VTugpQFcSdSOvQJGTjWe/OZn4d3fGbHPn0rjmIugagfekZ/vUOTtlGxJ4v+dmEvJlqQ2SyVbk7av7DUZNYsmeygpLaHwnx4M3+MYmkvmvWuYcP4c/vnJlxmSlRu+Nnjug2TNf4jrJzqZluPkviI3t09LItWl+Pk2P0u3+MLC/e2ZnnBqZqwdsQyUCpoJ4+H3NG2NINrLyrF/Pp15g94gXign1qu2fwYTnjnwJuwJeegv3wdfXmOOc6eZ8zqzljeWZKKAz8yBzEIT4nnvhVogk7HFMPJCI/6H3oX86SYclD89MrHbWGM6hLHF5trW5ZFRiHj+IQovM3FzRUR0rQncglCl8oIZZoJ2xCTzssI1Q3Mh92xK772VZ/60kfVP/isTANyDoOJvfHykjtkP/Au3zp5E6eWjzAggcwwc3AYjClnB5RzDy0XayYzmU5GRArSd5km09x7Pw4bo9M7aJs3y7V42HfRz8pKHGFwP/s2/Zshda/jC5yZyQb6TpWVjSL5jDQ0/vpqc2TfDFQ/xtSkuUlyKpVt83DbVFd4Zq9GvQUcycW6f5uE/321haZmPRr/mX6e3HcrpTOG1/oAIfjcTG4rpKNRzJvMGvUG8uLw9hj5tkfGuK9ZGrlux+JEXANrE7Hf9Bd55yoj6c1tLWb35GW5gPZoJvPEI1JbDyn/7mKUfz2ZO8a3c/1wpYMS+Yq2Z6B1bbDz8GXcau958wrR1UQmMvSK6U0rYOH9sNk5Kmsm8sbx7aB0j3/pi9ARt/meN4J+sovS+b7DkuQ0AzL7rB6y/41wmZMLHR+qZ/eM3OXSymSUvvgctDZR+ORcmXW4mc30tXM+rbM93MqPyD7BpT+vc/3g2A+P2vsb/XnoZpKSQkRwR0IrjQZZsaglnx1gdwM+2hmrjh5j7jW8z9s47eHbvUNwOWDjOTZMP1nwyHsfDb3HbhTlkppiJ2eXbTWeyvTbIxoNmtGAJujW6sJdj7mhPooFSakEEv5uJFfi+9tS7E7vXv3W5EeRxc03Y5b1fwtAxcHIP1HwEnjSo2we/32ImXf/wUSmr9pgeYUXabL5+1XrGzZnAR9s/5rEtsznWcIjlh5cw9FswbV8pQ8+CMcWmvd3rInH96h2m87C8/VgvPmFLJ8fmxaektRZ4e359U13k88evm3BMSjrkT6X0xS1hsQc4dPgos/9zPcv++TwW/3YLh042h68tWWvidqXjLgJ3MpStJBWYce51MCo0MZw9DvwtZk7BSv+022xhs9UuoNZCK3t2jFWX/rapLpr8mk9Oau6fnsSGA1nM8gXYWBnk4Y1mhLO/Hi4/ezi3TzOCvrLcx8JxblJdZmtDe2lk67o1ulg0xTMgi6S1hQh+NxMr8H3tqXcndq/fEtSJ8+GPi+DgW+azexA0Hzev5AzzvjmjlL/siWwBWVN3iLnXzubeK5fx5P7FnOJQ+NqPf7mEi4HZlAImFDS2GHKmwZGtRviDXhP2eefpiNdvt9HqkBIqrBMvL77ppMm6sUI4KWnGU92y0ghwylAzoTvuIhN2yRhF7ZvP88zqv0U9eiEL2XByA/Oeep100lnIwqjaNM+8uY+7jhwiMzPTzAc4XDD+okhapjvZtGkdW2Gdwsto8GtW+maxcLybYdbPEYM91BNbc+fbM02n8OsdXjwOk7Z537luPE4V6iRgVn5kQ5SfbfXGmYyN9szt4ZmB4rl3FhH8bubTJPB2YmPjlvi/+URE7FMy4ZrfwItfgaZaGP5ZKH+9ljL/M1HPWshCNhzawAPL58UVkPd4huKRdzHjhkw0RuRHTDOi7/eB3xsJH7kHtR168jWY6wkh/ClpkRWvVirkttDvdNOvIl6/lUBSu8fE3A/+HXIKTd78R6+SOX4a65+cyex7f8ShmhMd1qbJG5rM+nsvJPNEORy07VhX8QZMv8EcWyMJTav5hJX6IpZu8VFYvYELLr0ibtqmPTvGEux4JQ6uKHBxfl7EW7/9VbPf7fRcZ9Retfb3eHzaRN6OCL7QingTn23FxqctgpP7oOJluOopOPoB3PASvPGICb1M2JjJtVPWM/cfZnOssePiVsNS87gtaz2B/Zkc3gqBUIj2oxfheIWZ5K0tN+Ge0RfGD91Y57wNCRbPt5dJsDz8PW9HvP5p1xiP3pVs0ijra0Kib1tRC0wA1v/bHGY/uooNJ9uuTZOXnsr6e85nQs5gqK+GQdmm7YbqSMzbitOP+jxseR6aT0WFm64r0kzOb2TG/j/ABz6YcQM01dGwfR0r9UUsmJIZleYY631bxHrrtU2a6blOpo9wsmhyRLy7Q8wHSgpmPCQtU2hFvJTMaYuis2SsTcRTs2Du03D3LiP26x4wYl+xFvZvNPdXPDWBGxvXM4Q8NrAhvOfvcpaHV0VuYANDyOPGxvW4908gOcNk9ex/HTLGRYt9xjiT6tnWPJo1+jjvzvglGxKClDQTsqmrip4wtTqFqg9N3H5IjgnDDBluFlydMxfSRjBhUBPLbpzKCU6whCVRj7Zq0yz7xc+Z8LnzooW+odq0N+Vqc86K069/0nQom35lzhdeBjNvYtDk4lA9eiKjj/LXGFT2P1RuWRdOc7QmUqHjXakgshVhqlu1uretlM/OMlBSMOMhHn4v09cLqjpDvIlPe2zcXuzMHi+37h89CwJeOLEPVn7JiPbonAl88cgyfse8Notb3chvyDIJgDQfj7RdeI3RlInzYe2dpiN4t8JcO/SuWdi1c3Xr0E1HK4A/9cQueHr3OdNLjr/ITNzmTzULrU4dMXn2M28y99Ud5mNGsvgPb7Rbm2bxXfey/sXfMcHxFhRMB5cHXEnm+VZoyWo792zj4Z/7TyY7qPCyyGjEGnXY7G3wa/L1RSywlSZua9MUaO11txe66ehZHTFQUjDjIYLfy/R1dcvO5KnHCqX1HW9DJB3S8pxjQz1WuuaedZF8fIAjTR/zJxa3KyB/YjE3sT4s+jmfMyt0J10Lo0JbCI++0Dx39MXgdEfSQq33hBb4WOwhHnsKphW+yZ8Kg7NNOGZIjhHck1V83DSI2T9YwaHq4yxsrzZN9Spmf/FLrL97JhOmXwxXfyfSlj07yLLh6u+0vmbZGVNHZ9D0a/ma7Ue5osDFW4cC4T1sY4kV8fZCN10V7IEc4xfB72X6Ok3zdPPUG2si+fYXl7ROh7SPBuz3jimGzIkmzHOk6WN++O5sTnGofQFhFc8yOyz6J/cZT/+NR8wCrsYa46BeXGKyc6yfZ+L8SG2ehKetuvThlMwW89nXbLx7zyDzOXkwrPshtcfrmP29F8Kplx3Vpjl0vIHZP9nE+89cROa7zxlv3T6ysC/2qnjDtB+viFoH1TStFbbn5/lbZdVAxx2Cne4U7IEWzxfB72X6OovndPLU7QKeWQiTbzCxczv20cCbT0Ry86951ojxGz+tZZlndjj1siMBOcUhnmU2t/E+HM8kNdtM/kLoeUtMp2PvcBIuBbM92hLOlDSTNfPuc8bTH5JjznsbTDmFalNPJ3NwNrdeeg5LXtwS/uoqVpGXncFvHryfxU/8hlVV0bVpbr1gNJnV70H1e+aEvYyy5dHbK3PGK6LWQb2djrzy2A6ht4S4q+Gh3kYEP8E4nbi2tbjKmizduRqy2vnuxPmm9PGVS007E+fDqF9lMrX8Vl63TfytYhVpzjxuSf8Nz9cujkrJBJg3/VYmpWaChn1vmMnfUTPN3EBmoXm32xg7Ykno8godFSqz5ilPHTElEVwumH4j7H7LZPTUVVF6+UhQiiUvmP2l87KGsn7dOiZ89vOsv+x6ZhdfzqHqYwCUfG0BpV88xxRLGzHJPN++AMxesG34uOjUTIgekbRTdqEjrzy2Q+gtIR5o8XxnaWlpX9sQZtmyZaWLFy/uazOEEKlZcHy3EfDMiUZA3amR6zU7zQggr8jc+87TsPW/ITXThFi2Locdz8Gk9EvwNcM+XgdgCHl8d9Z6/nXV+Yxp+iKvv/8CLfoUAF+5oISvzynl/Wdh2HgYfxVUvgP558FfH4bKt+BUFXz2y8aGrImQmh1t2zs/NZ1AaraJ+ScU7iTInWTeLZrqYPtaUz4hZ4KZXM2dDOd+CYIBGPlZ852AF7LHgyuJS7Ibwelh1+E6E6OfMAFyJ5F5aDNfTDvACx+e4L4b51K6+J9MNtAnb5oJ22nXQt0R480f2w9jZppqmilpkD/FvOy2/f0FeOe3xqb8KWf8Y6e6FdNznaS6jTc/Nt1BZrKZvLXO9QSx7fYVS5YsqSotLV3W0X1d8vCVUo8AC4AgcBS4WWt9SCmlgJ8Ac4HG0Pn3utKW0PvsXG08fHs54zefiHjOVp17MDH22LIjo2fBsHGmPPpsSlEKtuhn+HrKehremMDO1ZBUNYGvBtbzv67ZLLzwVsa8Xoq63ISFKtZC3UEzunj5PtPxQOQd4o9YEra8Qjya6lpvXG4tiLJPoIIJ9cy8KVwaufSey7jrC38mc8SoyGrdwsuYALz/1fFklr9oVtCOnGquVW4LV9kMt2nl/7dFJ2vVnC4DeWK1J1Fan/lvWimVprWuCx3fDZyttf6GUmoucBdG8M8DfqK1Pq+j5xUVFemysrIztkfoXmJDI28+YTzni0vMCtbRs8yE6pVLjacde/9vr44urgamWuaXf5PJ0Q9CE721Rsyn3F1L7VuZ+BpNOrYGPKmQPxNe/iYsWG7COu1V7Ey48E1b2MMkVkx/9LmR1bbx7oPW3ym6Do5WRIqf2b9vzQXk2CZ1RkwyHn5Hm5q31W57m6MI7aKU2qK1Lurovi55+JbYhxhEpJ9eAPyPNr3JZqVUulIqV2td1ZX2hN4l1nuOXcFa/Hik/HG8+69cCs0n4FQlFFwG6aNgxp2ZgMneaaw1owgzwZvJG0sinj2Y9M+jHxgP34rjWztkeRvMDlkJWx2zPWK3NbTeYwXVnrbZZPuvHDsP8OoPjehvXxMZHVj/008cghYTjmPU51rvn9uRffa0TaHH6fKkrVLq+8A/AyeB2aHT+cAB220HQ+daCb5SajGwGGD06NFdNUfoQSxBb6wx2Xzt1coH4/Unp8PBTdBwGBb+0py3Rgpbf2XE3JroBTPRG/Sa+jmK1uEZK1LqbzTPmTg/+rpAa5HvjKC2J8IjCk24xh4MGP152PFnE69LG2Fi9daELZgQ0KZfRe9tG88+oVfpUPCVUuuAEXEuPay1/qPW+mHgYaXUg8CdELOqpgO01suAZWBCOqfzXaFviK1ICdG18WO9fPt7zU745BU46yKTgZNZ2HqeYOSFMOpCkwYa24nMuNOEk3yJVifndOisyNuxZ9NYK2GtEcGUq02VS7tAb3neiH16PiwMbS/4wZpIjr193sBakNUV+4RuoUPB11oXd/JZvwXWYgS/EhhluzYydE74lGAPpUAk/z7W086aGB32efk+s1J2TLEJCU2cHymLYD3XCu1ApISDJer2UYZ7kHj23UK8mD9Er4SNFegLvhZ5t/Lt7RuV268L/YauZumM11rvCn1cAFgbsK4G7lRKrcBM2p6U+P2nCyvnfuJ8k4YJnZs4tXv81iKurJh5gr0bTAeSN73t4mcJXyenO4kX84/19GMnYTPyzSTu9jWR+jz272fkt/bshT6nqzH8/1BKTcSkZe4DvhE6vxaToVOBScsUP+xTRmzKZkfia4/x2z3+WFKzIqt0JfOml4gX84+teRPbKXywBo6Um9LKEL2xidBv6WqWzj+0cV4Dd3Tl2UL/5nRz3U8nm0a8914mXsgmdmLV/l7+WmQHq5FTzQYqMgE7IOhSHn53I3n4n14kX74f0F5u/Ok+54M1JmVq8tWSP98P6JU8fEHoLOK19wM6qEjZaVLSzM5UwoBDBF8QEgXJf094RPAFIVGQ/PeER/a0FQRBSBBE8AVBEBIEEXxBEIQEQQRfEOvfNoAAAAUbSURBVAQhQRDBFwRBSBBE8AVBEBIEEXxBEIQEoV+VVlBKVWOKsLVHFlDTC+Z0JwPRZhiYdg9Em2Fg2j0QbYaBaXdHNp+ltc7u6CH9SvA7g1KqrDM1I/oTA9FmGJh2D0SbYWDaPRBthoFpd3fZLCEdQRCEBEEEXxAEIUEYiIK/rK8NOAMGos0wMO0eiDbDwLR7INoMA9PubrF5wMXwBUEQhDNjIHr4giAIwhkwIARfKfWIUup9pdTflVKvKKXyQueVUupJpVRF6Prn+9pWO0qpJ5RS5SHbXlRKpduuPRiye6dS6sq+tNOOUuoflVI7lFJBpVRRzLV+abOFUmpOyLYKpdS3+tqetlBK/UopdVQptd12bphS6lWl1K7Qe0Zf2hiLUmqUUmq9UurD0N/HPaHz/dZupVSyUuodpdS2kM1LQufHKKXeDv2drFRKefra1ngopZxKqa1KqZdCn7tut9a637+ANNvx3cAvQsdzgT9jNlubCbzd17bG2H0F4Aod/wD4Qej4bGAbkASMAT4BnH1tb8i2ScBEYANQZDvfb20O2ecM2TQW8IRsPbuv7WrD1ouAzwPbbeceB74VOv6W9bfSX15ALvD50PEQ4OPQ30S/tTukC4NDx27g7ZBOPA9cHzr/C+C2vra1DfvvA34HvBT63GW7B4SHr7Wus30cBFgTDwuA/9GGzUC6Uiq31w1sA631K1prf+jjZmBk6HgBsEJr3aK13gNUADP6wsZYtNYfaa13xrnUb20OMQOo0Frv1lp7gRUYm/sdWus3gGMxpxcAof0HeRZY2KtGdYDWukpr/V7o+BTwEZBPP7Y7pAv1oY/u0EsDlwJ/CJ3vVzZbKKVGAlcD/x36rOgGuweE4AMopb6vlDoAfBn4buh0PnDAdtvB0Ln+yNcwoxEYWHZb9Heb+7t9HZGjta4KHR8GcvrSmPZQShUA0zAec7+2OxQW+TtwFHgVMwo8YXPE+uvfyY+BB4Bg6HMm3WB3vxF8pdQ6pdT2OK8FAFrrh7XWo4DfAnf2rbUROrI7dM/DgB9je5/TGZuFvkObMXu/TJ9TSg0G/g+4N2bk3S/t1loHtNafw4yuZwCFfWxShyil5gFHtdZbuvvZ/WZPW611cSdv/S2wFigBKoFRtmsjQ+d6jY7sVkrdDMwDLgv9h4A+tvs0ftd2+vx33QH93b6OOKKUytVaV4XCkkf72qBYlFJujNj/Vmv9Quh0v7cbQGt9Qim1HjgfE/p1hbzl/vh3ciEwXyk1F0gG0oCf0A129xsPvz2UUuNtHxcA5aHj1cA/h7J1ZgInbcPLPkcpNQczLJuvtW60XVoNXK+USlJKjQHGA+/0hY2nQX+3+V1gfCiTwQNcj7F5oLAauCl0fBPwxz60pRWhGPIvgY+01kttl/qt3UqpbCszTimVAlyOmXtYD3wpdFu/shlAa/2g1nqk1roA83f8V631l+kOu/t6JrqTs9X/B2wH3gf+BOTryCz8TzFxuQ+wZZX0hxdmYvMA8PfQ6xe2aw+H7N4JXNXXttrsugYTH2wBjgAv93ebbfbNxWSPfAI83Nf2tGPnc0AV4Av9rm/BxGhfA3YB64BhfW1njM1fwIRr3rf9Pc/tz3YDnwW2hmzeDnw3dH4sxlmpAH4PJPW1re38DJcQydLpst2y0lYQBCFBGBAhHUEQBKHriOALgiAkCCL4giAICYIIviAIQoIggi8IgpAgiOALgiAkCCL4giAICYIIviAIQoLw/wGGA0sgGddO8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": {},
   "outputs": [],
   "source": [
    "from numpy import exp, sqrt, array, abs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "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": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FdX5+PHPk5sNQoCQBAgJWQhhXyWCCIIrsqho1Yq71ZavVmtbXKrfb3+1xa/furRKVaxQtWqtpS614oogiKAiBASRJRBCgLCGsEPWe5/fH3OBSANcwk0mufd5v17zyp2ZM3OfofWZueecOUdUFWOMMeEhwu0AjDHGNBxL+sYYE0Ys6RtjTBixpG+MMWHEkr4xxoQRS/rGGBNGLOkbY0wYsaRvjDFhxJK+McaEkUi3AzhWUlKSZmZmuh2GMcY0KYsXL96pqsknK9fokn5mZiZ5eXluh2GMMU2KiGwIpJxV7xhjTBixpG+MMWHEkr4xxoSRRlenb4wJX1VVVRQXF1NeXu52KI1WbGwsaWlpREVF1en4gJK+iIwE/gR4gBdU9dHjlLsSeAs4U1Xz/NseBG4DvMDdqjqjTpEaY0JecXEx8fHxZGZmIiJuh9PoqCqlpaUUFxeTlZVVp3OctHpHRDzAZGAU0AO4VkR61FIuHvg58HWNbT2AcUBPYCTwnP98xhjzH8rLy0lMTLSEfxwiQmJi4mn9EgqkTn8gUKCqhapaCUwDxtZS7mHgMaBmNGOBaapaoarrgQL/+YwxplaW8E/sdP99Akn6qcCmGuvF/m01gzgD6KiqH5zqscHi9Sm//3AVb+Rt4puNu9lfXlUfX2OMMU3aaTfkikgE8CRwy2mcYzwwHiA9Pb1O59i+r5yXvyyiotp3ZFuHVrF0bhdPTtsWXNCtLYOz7WejMeb4ioqKuOSSS/juu++Cet7S0lKuuuoqFi1axC233MKzzz4b1POfikCS/magY431NP+2w+KBXsBn/oTaHpguIpcFcCwAqjoVmAqQm5tbp5naO7RuxsqJI9m06xBrdxxg7Y79rN3u/P3716W8OH89AzIS+Nn5nRneJdmSvzGmwcTGxvLwww/z3XffBf2GcqoCqd5ZBOSISJaIROM0zE4/vFNV96pqkqpmqmomsAC4zN97ZzowTkRiRCQLyAEWBv0q/DwRQmZSHBf1aMdPz+3MU9f04/2fncPS34zg4ct7sW1vObf8dRGXT/6CWSu3o1qn+4sxJoRVV1dz/fXX0717d6666irmzp3LD37wAwDeffddmjVrRmVlJeXl5XTq1Cmgc8bFxTF06FBiY2PrM/SAnPRJX1WrReQuYAZOl82XVHWFiEwE8lR1+gmOXSEibwArgWrgTlX1Bin2gMVGebjxrAyuye3Iv5YUM/mzAn78ah49O7Tk3ou7cl7Xtg0dkjHmJH733gpWbtkX1HP26NCShy7tecIy+fn5vPjiiwwZMoRbb72VRYsWsXTpUgDmzZtHr169WLRoEdXV1QwaNAiAJ554gr///e//ca5hw4bx9NNPB/UaTldAdfqq+iHw4THbfnOcsuces/4I8Egd4wuq6MgIxg1M58oBaby7dAuT5xTwo78u4q7zOvPLi7rgibAqH2PCXceOHRkyZAgAN9xwA08//TTZ2dmsWrWKhQsXMmHCBD7//HO8Xi/nnHMOAPfddx/33Xefm2EHLCzfyI3yRHDVgDQu7ZvCQ++u4Nk5BazYspdJ4/rTqlnd3nIzxgTXyZ7I68ux7X0iwrBhw/joo4+Iioriwgsv5JZbbsHr9fLEE08AIfikH6piIj38/ge96ZXait9OX8HYZ+cz9aZcurSLdzs0Y4xLNm7cyFdffcXgwYN5/fXXGTp0KAMGDOCmm27ipptuIjk5mdLSUrZv306vXr2ApvWkH/YDrokIN5yVwT/Gn8WBCi9XTP6Cj7/b6nZYxhiXdO3alcmTJ9O9e3d2797NHXfcwaBBg9i+fTvDhg0DoE+fPvTu3fuUegFmZmYyYcIEXn75ZdLS0li5cmV9XcIJSWPrwZKbm6tuTaKybW85//XaYpZt2sPd5zv1/Na105iGs2rVKrp37+52GI1ebf9OIrJYVXNPdmzYP+nX1L5VLP8cfxZXD0jj6dkFPDlzjdshGWNMUIV1nX5tYqM8PH5VHyJEeGZ2AQnNo7l1aN1GszPGmMbGkn4tRIRHrujF3rIqJr6/koS4KK7on+Z2WMYYc9qseuc4Ij0RTBrXj8GdErn3zW+ZvXq72yEZY8xps6R/ArFRHqbeNIAeKS2547UlLCra5XZIxhhzWizpn0R8bBQv/+hMUls349aXF7Fqa3BfCzfGmIZkST8AiS1i+NuPBxEXHclNLy1k854yt0MyxtSDoqKiIy9cBdPMmTMZMGAAvXv3ZsCAAcyePbvWcr/97W9JTU2lX79+9OvXjw8//LDWcqfDkn6AUls342+3DaSs0stdry+hyus7+UHGGAMkJSXx3nvvsXz5cl555RVuvPHG45b95S9/ydKlS1m6dCmjR48OeiyW9E9BTrt4Hr2yN99s3MMfZuS7HY4xph7Ux9DK/fv3p0OHDgD07NmTsrIyKioq6u0aTsS6bJ6iS/p04Mt1pUz5vJBBndpwfrd2bodkTGj66AHYtjy452zfG0Y9esIi9T208ttvv80ZZ5xBTExMrd//7LPP8uqrr5Kbm8sf//hHEhIS6nKlx2VP+nXwm0t60K19PPe8sYyte61+35hQcuzQyvPnz691aOV58+Z9b2jlw1UyNZdjE/6KFSv41a9+xZQpU2r97jvuuIN169axdOlSUlJSuOeee4J+ffakXwexUR4mX38Glz4zn7v/8Q3/+MlZRHrs/mlMUJ3kiby+1NfQysXFxVxxxRW8+uqrZGdn1/rd7dodrTn4yU9+wiWXXBKsyzoioEwlIiNFJF9ECkTkgVr23y4iy0VkqYjMF5Ee/u2ZIlLm375URJ4P9gW4JTu5Bf93RW8WFe22MXqMCSGHh1YGjgytfM455zBp0iQGDx58ZGjl/Pz87w2tfKIn/T179jBmzBgeffTRI78iarN169ERft9555166Ul00qQvIh5gMjAK6AFcezip1/C6qvZW1X7A48CTNfatU9V+/uX2YAXeGFzeP5Vrcjvy3GfrmLumxO1wjDFBUB9DKz/77LMUFBQwceLEI90xd+zYAcCPf/xjDo8sfP/999O7d2/69OnDnDlzeOqpp4J+fScdWllEBgO/VdWL/esPAqjq749T/lrgJlUdJSKZwPuqGvDtys2hleuirNLL5ZO/YOeBCj78+Tm0a+n+xMfGNFU2tHJg6nto5VRgU431Yv+2Y7/wThFZh/Okf3eNXVki8o2IzBWRc2r7AhEZLyJ5IpJXUtK0npibRXuYfH1/DlV6efBfy2ls8xMYY0xNQWt9VNXJqpoN/Ar4tX/zViBdVfsDE4DXRaRlLcdOVdVcVc1NTk4OVkgNpnPbeO4Z0YXZq3fw3rc265YxpvEKJOlvBjrWWE/zbzueacDlAKpaoaql/s+LgXVAl7qF2rj9aEgWfdNa8bvpK9h9sNLtcIxpsuzX8omd7r9PIEl/EZAjIlkiEg2MA6bXLCAiOTVWxwBr/duT/Q3BiEgnIAcoPK2IGylPhPDolX3YW1bFw++7M/elMU1dbGwspaWllviPQ1UpLS0lNrbubYcn7aevqtUichcwA/AAL6nqChGZCOSp6nTgLhG5EKgCdgM3+w8fBkwUkSrAB9yuqiE7PnH3lJbccW42z8wuYGz/VIZ3aXpVVca4KS0tjeLiYppa215Dio2NJS2t7pM62cToQVZe5WX00/OoqPLxyS+HERdj778ZY+qfTYzuktgoD49d2YfNe8r44yf20pYxpnGxpF8Pzsxsw41nZfDXL9fzzcbdbodjjDFHWNKvJ/eP7Er7lrE88PZyKqtt7H1jTONgSb+exMdG8b+X9yJ/+36en7vO7XCMMQawpF+vLujejjG9U3juswKbYtEY0yhY0q9n/z2mO6rw+w9XuR2KMcZY0q9vqa2bcfvwbN7/disL14fsKwrGmCbCkn4DuH14NimtYvndeyvw+hrXexHGmPBiSb8BNIv28ODo7qzYso838jad/ABjjKknlvQbyKV9UjgzM4E/zMhnb1mV2+EYY8KUJf0GIiI8dGlPdh2q5JlP17odjjEmTFnSb0C9UltxTW5HXv6yiIIdB9wOxxgThizpN7B7L+5KsygP//uBDb9sjGl4lvQbWFKLGH5+YQ6f5ZcwZ/UOt8MxxoQZS/ouuGlwJp2S4nj4/ZVUeW1cHmNMw7Gk74LoyAgeHN2dwp0HrQunMaZBBZT0RWSkiOSLSIGIPFDL/ttFZLmILBWR+SLSo8a+B/3H5YvIxcEMvim7sHtbcjMSmDRrLYcqq90OxxgTJk6a9P1z3E4GRgE9gGtrJnW/11W1t6r2Ax4HnvQf2wNnTt2ewEjgucNz5oY7EeHB0d0o2V/Bi/PWux2OMSZMBPKkPxAoUNVCVa0EpgFjaxZQ1X01VuOAw2MNjAWmqWqFqq4HCvznM8CAjDaM6NGOKZ8XUnqgwu1wjDFhIJCknwrUrHgu9m/7HhG5U0TW4Tzp332Kx44XkTwRyQu3CZHvH9mVQ5XVPDunwO1QjDFhIGgNuao6WVWzgV8Bvz7FY6eqaq6q5iYnJwcrpCahc9t4fpjbkdcWbGDTrkNuh2OMCXGBJP3NQMca62n+bcczDbi8jseGpV9c2AVPhPCHT/LdDsUYE+ICSfqLgBwRyRKRaJyG2ek1C4hITo3VMcDhwWWmA+NEJEZEsoAcYOHphx1a2reK5dYhWby7dAvfbd7rdjjGmBB20qSvqtXAXcAMYBXwhqquEJGJInKZv9hdIrJCRJYCE4Cb/ceuAN4AVgIfA3eqqrcerqPJ+6/h2bRuHsVjH692OxRjTAgT1cY1qUdubq7m5eW5HYYrXphXyP9+sIrXbhvE0Jwkt8MxxjQhIrJYVXNPVs7eyG1Ebjgrg9TWzXj041X4bIYtY0w9sKTfiMRGeZhwURe+27yPj1dsczscY0wIsqTfyFzeP5XObVvw5Mw1Np+uMSboLOk3Mp4IYcJFXSjYcYB3l1rvVmNMcFnSb4RG9mxPzw4tmTRrrQ29bIwJKkv6jVBEhHDPiC5s3HWIN/OK3Q7HGBNCLOk3Uud1bcsZ6a15ZvZayqvs1QZjTHBY0m+kRIR7R3Rl695yXv96o9vhGGNChCX9RuzszkmcnZ3Ic58V2EQrxpigsKTfyN0zois7D1Ty8pdFbodijAkBlvQbuQEZCZzfrS1T5hayt6zK7XCMMU2cJf0m4J4RXdhbVsWL821aRWPM6bGk3wT07NCKMb1TeHFeIbsOVrodjjGmCbOk30T88qIcyqq8TPl8nduhGGOaMEv6TUTntvGM7ZfKq19uoGS/TaJujKmbgJK+iIwUkXwRKRCRB2rZP0FEVorItyLyqYhk1NjnFZGl/mX6sceawN19QQ6VXh9T5trTvjGmbk6a9EXEA0wGRgE9gGtFpMcxxb4BclW1D/AW8HiNfWWq2s+/XIaps6ykOK7on8rfFmxg+75yt8MxxjRBgTzpDwQKVLVQVStxJj4fW7OAqs5R1UP+1QU4E6CbenD3+TlU+5Q/f2ZP+8aYUxdI0k8FNtVYL/ZvO57bgI9qrMeKSJ6ILBCRy+sQo6khPbE5Vw9I4/WvN7J1b5nb4RhjmpigNuSKyA1ALvBEjc0Z/nkbrwMmiUh2LceN998Y8kpKSoIZUki66/zOKMrkOQVuh2KMaWICSfqbgY411tP8275HRC4E/ge4TFWPdC9R1c3+v4XAZ0D/Y49V1amqmququcnJyad0AeEoLaE515zZkX8u2kTx7kMnP8AYY/wCSfqLgBwRyRKRaGAc8L1eOCLSH5iCk/B31NieICIx/s9JwBBgZbCCD2d3ntcZQXh2tj3tG2MCd9Kkr6rVwF3ADGAV8IaqrhCRiSJyuDfOE0AL4M1jumZ2B/JEZBkwB3hUVS3pB0FKq2ZcNyidNxcXs6H0oNvhGGOaCFFtXJNv5+bmal5entthNAnb95Uz7PE5XNq3A3+4uq/b4RhjXCQii/3tpydkb+Q2Ye1axnLDWRn8a0kxhSUH3A7HGNMEWNJv4m4fnk1MpIc/fbrW7VCMMU2AJf0mLjk+hpvOzmD6si2s3b7f7XCMMY2cJf0Q8F/Dsmke5WHSLHvaN8acmCX9ENAmLpofDcnig+VbWbV1n9vhGGMaMUv6IeIn53QiPjaSp2aucTsUY0wjZkk/RLRqHsWPh3bik5XbWV681+1wjDGNlCX9EPKjoZm0ahbFU7Psad8YUztL+iGkZWwU44d1YvbqHSzZuNvtcIwxjZAl/RBzy9mZJMZFW92+MaZWlvRDTFxMJLcPz2be2p0sXL/L7XCMMY2MJf0QdMNZGSTHx/DkzHy3QzHGNDKW9ENQs2gPPz03mwWFu/iyYKfb4RhjGhFL+iHq2oHppLSK5Q+f5NPYRlI1xrjHkn6Iio3ycNf5nVmycQ9z8nec/ABjTFiwpB/CfpjbkfQ2zXlixhp8PnvaN8YEmPRFZKSI5ItIgYg8UMv+CSKyUkS+FZFPRSSjxr6bRWStf7k5mMGbE4vyRDDhoi6s2rqPD7/b6nY4xphG4KRJX0Q8wGRgFNADuFZEehxT7BsgV1X7AG8Bj/uPbQM8BAwCBgIPiUhC8MI3J3Np3w50bRfPk5+sodrrczscY4zLAnnSHwgUqGqhqlYC04CxNQuo6hxVPeRfXQCk+T9fDMxU1V2quhuYCYwMTugmEJ4I4Z4RXSjceZB/LdnsdjjGGJcFkvRTgU011ov9247nNuCjOh5r6sFFPdrRt2NrJs1aQ0W11+1wjDEuCmpDrojcAOQCT5ziceNFJE9E8kpKSoIZkgFEhPtGdGXL3nJe/3qj2+EYY1wUSNLfDHSssZ7m3/Y9InIh8D/AZapacSrHqupUVc1V1dzk5ORAYzenYEjnRAZ3SmTynAIOVVa7HY4xxiWBJP1FQI6IZIlINDAOmF6zgIj0B6bgJPyancJnACNEJMHfgDvCv800MBHh3ou7svNAJX/9osjtcIwxLjlp0lfVauAunGS9CnhDVVeIyEQRucxf7AmgBfCmiCwVken+Y3cBD+PcOBYBE/3bjAsGZCRwYfe2TJm7jr2HqtwOxxjjAmlsr+jn5uZqXl6e22GErJVb9jH66XnceV42913cze1wjDFBIiKLVTX3ZOXsjdww06NDSy7t24GX5hexY1+52+EYYxqYJf0wdO+ILlT7fDw1a63boRhjGpgl/TCUkRjH9YMyeCNvEwU79rsdjjGmAVnSD1M/O78zzaI8PPaxTbRiTDixpB+mElvEcPvwTsxcuZ28IutQZUy4sKQfxm4dmkXb+Bj+78NVNtGKMWHCkn4Yax4dyYSLurBk4x5mrNjmdjjGmAZgST/MXTUgjc5tW/D4x/lU2dDLxoQ8S/phLtITwQMju1G48yDTFm06+QHGmCbNkr7hgu5tGZjZhj/NWsOBChuMzZhQZknfICI8OLobOw9U8pfPC90OxxhTjyzpGwD6pycwund7/jKv0IZnMCaEWdI3R9x/cTeqvD4en2EvbBkTqizpmyMyk+K4dWgWby0uZtmmPW6HY4ypB5FuBxA0h3bB5EEgERDhAfFARISzLh7wRENULETGQmTM0b9RzSEmHqJbOH9jWkBMS2dpluAszdtAbGvwhM4/1/HcdV5n3l68md+9t4K37zgbEXE7JGNMEIVOFouIhG5jQL3g8zl/1Qc+r/PZWwXV5VBdARX74WCJ87nyEFTud7bpSfqpx7aCZm2gRTtokez/2w5atHX+xqdAqzRonghNNFnGx0Zx/8iu3P/Wt0xftoWx/Wwee2NCSUCTqIjISOBPgAd4QVUfPWb/MGAS0AcYp6pv1djnBZb7Vzeq6mWcgGuTqKhCVZmT/Cv2Q8VeKNsNh3ZD2S7/511wqBQO7oADO+DAdmf7sSJjoWUqtEqFlmnQOh0SMo8uLdo5v0IaKZ9PGTv5C0r2VzD73uE0jw6dZwNjQlWgk6ic9L9mEfEAk4GLgGJgkYhMV9WVNYptBG4B7q3lFGWq2i+gqN0kAtHNnSW+XeDHVVc4vxr2b4d9m51lb/HRv4Wfwf6tQI2ba2QstM6ANlmQ2NlZknIgMcf51eDyr4SICOGhS3tw1fNf8fxn65gwoqur8RhjgieQR7iBQIGqFgKIyDRgLHAk6atqkX9f+L3HHxnjVOm0SgMG1F6mugL2bILdRbB7vf9vEewqhHVzwFtxtGxMS+cGkNwd2naDtt2dzy07NOjNIDezDZf17cCUzwu5OrcjHds0b7DvNsbUn0CSfipQ8/38YmDQKXxHrIjkAdXAo6r671M4NjRExkBSZ2c5ls8HezdB6VooXQc718LOfFj7CSx97Wi5mFbODaBdT2jfG9r3cdaj6y8ZPzCqG5+s3MajH61m8vVn1Nv3GGMaTkNU1mao6mYR6QTMFpHlqrquZgERGQ+MB0hPT2+AkBqRiAhIyHCWzhd+f9/BUihZBTtqLMvfhLwXnf0S4VQNte8NKX2hQ39I6QexLYMSWofWzbhjeGeemrWGGwtLOatTYlDOa4xxTyBJfzPQscZ6mn9bQFR1s/9voYh8BvQH1h1TZiowFZyG3EDPHfLiEiFuKGQOPbpNFfZsgG3LYdt3zt9NC+G7t4+WScxxbgAd+kPqGc4NIapZnUIYP6wTb+Rt4nfvreT9nw3FE9E0eyUZYxyBJP1FQI6IZOEk+3HAdYGcXEQSgEOqWiEiScAQ4PG6Bmtw6vUP9wLqfunR7QdLYcs3R5ei+bD8DWdfRKTzayDtTP+SCwlZAbURNIv28ODobtz1+je8/vUGbhycWR9XZYxpIIF22RyN0yXTA7ykqo+IyEQgT1Wni8iZwDtAAlAObFPVniJyNjAF8OG8/TtJVV880Xe51mUzFO3fBpuXQPEiZ9m8BKoOOvuaJ0HHQZB+lrOk9IPI6FpPo6rc8OLXfFu8l0/vGU7b+NgGvAhjTCAC7bIZUNJvSJb065HP67QLbM6DjV/DpgVODyJwupF2OMO5AWQMgfRBzhvKfoUlBxg5aR4X92rPM9f2d+kCjDHHY0nfBObADti4wFk2LYCty8BX7QxdkdIXModAxlBIP4tJX+xg0qy1vHLrQIZ3SXY7cmNMDZb0Td1UHIDihVD0BWz4AjYvBm8lSAS+9n35584s8iJ688jdPyY2Lji9hIwxp8+SvgmOqjKnPaDoCyiah2/TQiJ8VXglEk/HMyFrGGQNdxqIj9MmYIypf5b0Tf2oPMjzr72OrP+cm1M2EFuy3BmoLrqF0xaQfR50Og+Su7o+nIQx4SRoY+8Y8z3RcVz9w5u44MkMPo2IZ9p93YjYMN8ZY2jdHFg7wykXn+Ik/84XQPb5zvDUxhjXWdI3pyyxRQz/Pao797/9LW+tOMAPz7z06DsDezY6yb9wDqz5CJa9DojzoljnC52bQGpuWMxNYExjZNU7pk5UlWumLGDNjv18OmE4iS1i/rOQzwtblkLBLFj3qdM2oD5nXoJO50HOCOdGcCqjmhpjamV1+qberd2+n9FPz2N07xT+NC6Avvtlu6FwLhTMhLWz4MA2Z3tKX+cGkDMCUgc4M58ZY06JJX3TIP40ay1PzVrDn68/g1G9UwI/UNUZN6hgJqyd6YwfpF5n1rHOF0KXiyH7AmjWuv6CNyaEWNI3DaLK6+MHz33J5j1lzPjFMJLja6nmCUTZblg3G9Z84gwrXbbLeUEsfTB0GQFdRkJSF+sRZMxxWNI3DWbt9v2MeWY+53ZJZsqNA05/MnWf13kpbM3Hzk1gu3+2zYQs6DrKWdIHgyfq9IM3JkRY0jcN6i+fF/LIh6t48od9+cEZacE9+d5i5waQ/xGs/9x5QzimFeRcCF1GQc5FVg1kwp4lfdOgvD5l3NSvWL1tP5/8chgpreo2fv9JVRxwuoPmf+y8E3CwxBk6OuNs6DoGuo50hp02JsxY0jcNbkPpQUZOmkduZgKv3jrw9Kt5TuZwNVD+h86vgJLVzva2PaHbaOg62nk/wNoBTBiwpG9c8bcFG/h///6OR67oxfWDMhr2y0vXOck//0PY+JXzTkDLVKcNoNsYZ7RQGx/IhChL+sYVqspNLy1k8YbdfPTzc8hIjHMnkIOlTvXP6g+g4FOoLvO3A1zk3AByLvrefAHGNHWBJv2IAE82UkTyRaRARB6oZf8wEVkiItUictUx+24WkbX+5ebAL8E0RSLCY1f2wSPChDeWUe31uRNIXCL0uw7G/R3uL4Rx/4AelzrtAW/9CB7vBH+/Gha/DPu3uxOjMS446ZO+iHiANcBFQDHOnLnXqurKGmUygZbAvcB0VX3Lv70NkAfkAgosBgao6u7jfZ896YeGd5du5ufTlnL78GweGNXN7XCO8nlh09fOL4DV78PuIkCg40DnF0C3SyAx2+0ojTllwRxlcyBQoKqF/hNPA8YCR5K+qhb59x37WHcxMFNVd/n3zwRGAv8I4HtNEza2Xypfr9/F83PXMTArgfO7NZLxdSI8Tk+fjLNhxP/CjpXODWDVezDzN86S3B26X+LcBFL6WUOwCSmBVO+kAptqrBf7twUioGNFZLyI5IlIXklJSYCnNo3dby7pQY+Ulkx4Yxmb95S5Hc5/EoF2PWH4/XD7PPjFchj5GMQlwbw/wtRz4ale8OH9/vcDqt2O2JjTFlCdfn1T1amqmququcnJNvdqqIiN8vDc9WdQ7VXuen0JldUu1e8HqnU6nHU73PI+3FsAY5+DlD6w5BV45VL4Q2d45w7nl0HlIbejNaZOAkn6m4GONdbT/NsCcTrHmhCQmRTHY1f24ZuNe3j849VuhxO4uETofz1c+w+nIfiHf4OciyH/A5h2HTyRDdOuh2XT4NAut6M1JmCB1OkvAnJEJAsnYY8Drgvw/DOA/xORBP/6CODBU47SNGlj+qSwcH0GL8xfz8CsNozo2d7tkE5NdBz0uMxZvFVQNN9pBF79ofNXPJA51GkE7jYGWgVa+2lMwwuon76IjAYmAR7gJVV9REQmAnmqOl1EzgQhvGNSAAARg0lEQVTeARKAcmCbqvb0H3sr8N/+Uz2iqn890XdZ753QVFHt5ernv6Jo50E+uPscOrZp7nZIp8/ng63fwKr3nSqfnfnO9g79/T2BLrW5gk2DsZezTKOzadchRj89j6ykON74r8HERoXYZCk71zq9gFZ/AJv9/x9uk320K2jamRDRKJrRTAiypG8apZkrtzP+b3mM6Z3C0+P6ExERok/B+7Y6w0Gs/sDp+eOrgri2R4eEyBoOUbFuR2lCiCV902hNmbuO33+0mrvP78yEEV3dDqf+le91Zgdb/b4zTWTlfoiKcyaJ73aJM0lMs4STn8eYEwjmy1nGBNX4YZ0oLDnI07MLyEyKC/74+41NbCvofZWzVFdA0TznF0D+R7BqutMQnHG28wug6ygbGtrUK3vSN66orPZxs39gttd+PIiBWW3cDqnh+Xyw5RunG+jqD6FklbP9yNDQoyClv7UDmIBY9Y5p9PYequKKP3/B7oOVvPPTIWQmuTQiZ2Oxq9BJ/vkfwcYvnaGh41OcSeK7jIJOwyGqnianMU2eJX3TJBTtPMjlz31Bm+bRvPPTIbRqbvPeAs4LX2tmwJqPnKGhKw9AZDPIPs/5BZBzMcQ3kvGMTKNgSd80GQvX7+L6FxaQm9GGV24dSHSkVWd8T3WF80JY/kfOXMF7/cNZpQ6ALiOdXwLt+9j7AGHOkr5pUv61pJgJbyxjZM/2PHtdfyI9lvhrpQrbv/NPFP+xM10kCvEd/NVAIyFrGESHwMtv5pRY0jdNzovz1/Pw+yu5rG8HnrqmH55Q7cMfTAd2ON1B13wM62b7q4FiIfMcyBnhdAe13kBhwZK+aZL+/Nk6Hvt4NVcNSOPxK/uE7stb9aG6AjZ84b8JzIBd65ztSV2d6SFzRkD6YJsnOERZ0jdN1qRZa5g0ay3XDUrnkct7IVZXXTel65zkv3YGFH3hvBUc3cJ5GzjnImdpFeLvSIQReznLNFk/vyCHimoff/5sHTGREfzmkh6W+OsiMRsG/9RZKg44w0EUzHR+CeR/4JRJ7g45F0L2Bc6vABsaIuRZ0jeNjohw/8Vdqaz28eL89URHRvDAyG6W+E9HTAvnha9uo53G4JL8ozeAr6fAl884XUKzznFuAJ0vdG4a9m8ecizpm0ZJRPj1mO5UVHuZMrcQn095cFR3q+MPBhFo281Zzv4ZVB50uoQWfAoFs2DtJ065VumQfS5kn+9UCTUPw7emQ5AlfdNoiQgTL+uFR4S/zFvPjv0VPHFVX+vHH2zRcf7unhc767uLnOS/bg6s+DcseRUQZ56A7POg03nQcSBExrgZtamjQCdRGQn8CWcSlRdU9dFj9scArwIDgFLgGlUtEpFMYBXgn12CBap6+4m+yxpyzbFUlT/PXcfjH+czpHMiz98wgPhYe3O3QXirYcsSpzvoujlQvAjU61QFZQyGTuc6vwLa97ExglwWtN47IuIB1gAXAcU40ydeq6ora5T5KdBHVW8XkXHAFap6jT/pv6+qvQIN3JK+OZ63Fxfzq7e/JaddPK/86EzatrRGxwZXvtfpCbR+LhR+BiX+eY+btXHaA7KGQeYwSMqx9oAGFszeOwOBAlUt9J94GjAWWFmjzFjgt/7PbwHPirW6mSC7ckAaSfEx3PHaYq547kteuXUgndu2cDus8BLb6miDMDiTxaz/3LkBrJ8LK991trdo78wbfPhGkJBlN4FGIpCknwpsqrFeDAw6XhlVrRaRvUCif1+WiHwD7AN+rarzTi9kE86Gd0nmn+MH86OXF3LV81/yl5tyOTPTGhhd0zIF+l7jLKrOSKFF82D9POdm8N1bTrn4DpA5xLkRZAy1nkEuqu+G3K1AuqqWisgA4N8i0lNV99UsJCLjgfEA6enp9RySaep6p7XiX3cM4ea/LuTaqQt4YFQ3bhuaZV063SbiJPPEbBhwi3MT2LnGSf4bvoDCubD8Tadsi3aQMcSZPCZ9MLTtYW0CDSSQOv3BwG9V9WL/+oMAqvr7GmVm+Mt8JSKRwDYgWY85uYh8BtyrqsettLc6fROovWVV3PfmMj5ZuZ2RPdvz+NV9aGkNvI2XKpQWON1DN3zhtA3s3+Lsi20FHc9yGofTz4YO/ax30CkKZkNuJE5D7gXAZpyG3OtUdUWNMncCvWs05P5AVX8oIsnALlX1ikgnYJ6/3K7jfZ8lfXMqVJUX5q3n0Y9Xk5bQjOeuP4OeHVq5HZYJhCrs2QAbvnImjdnwFZSudfZ5YpwuoumDoKN/iUtyN95GLqhj74jIaGASTpfNl1T1ERGZCOSp6nQRiQX+BvQHdgHjVLVQRK4EJgJVgA94SFXfO9F3WdI3dZFXtIu7Xv+GXYcqmXhZT645s6NV9zRFB3fCxq9g4wLY9DVsWeqMGQTQJttJ/mm5kHamUyXksVeNDrMB10zYKT1QwS/+uZR5a3dyeb8OPHRpTxLibETJJq2q3JlHeNPX/mUhHNrp7IuKg9QznJtAaq4zqUzLFHfjdZElfROWvD7l2dkFPDN7La2bR/HQpT25pE+KPfWHClXnjeHiPNic57wstvXbo78G4js4N4LUM5ybQIf+TntBGLCkb8Layi37eOBf3/Jt8V4u7N6Why/vRUorm1Q8JFWVw7ZvYfMSZyaxzYuPziUATrVQh37ODSClH6T0hdiW7sVbTyzpm7BX7fXx8pdF/OGTfCIjIvjVqG5cPzDdBm0LB2W7nWqhzYuddoEtS2Ff8dH9bbKd5J/SxxlCIqVvk28otqRvjN/G0kP89zvLmV+wkwEZCfx6THf6pye4HZZpaAd3+m8A38DWpU610N6NR/fHpzg3gPa9oX0vaNcL2nSCCI97MZ8CS/rG1KCqvLW4mMc+Xs3OA5Vc3LMd913czYZxCHeHdjkTzW/91qki2vqt80KZep39kc2gbXdo19O5CbTr4fQaaoS/CizpG1OLgxXVvDBvPVM/X0dZlZerB3TkFxflWH2/OaqqHHbmw7bvYPsK2L7c+VxW4/WiuGTnZtC2ByR38y9dXZ1zwJK+MSdQeqCCZ+cU8NqCDUSIcMvZmdx2ThZt423kTlMLVTiwHXas8i8rj36uOni0XFxbJ/knd3VuBEk5kJgDLTvU+1hDlvSNCcCmXYd4auYa3lm6maiICC7t24Fbh2baW70mMD6f0y5Qssb5dVCy2pmKsmQNVOw9Wi4qDpI6OzeApBxI7Oy0FyRmB61LqSV9Y05BYckBXv6yiDfziimr8nJWpzbcNrQT53dri8d6+5hTpQr7tzntA6VrYWfB0c97NgE18m5cstObKDHb6VY68Cd1+kpL+sbUwd5DVUxbtJFXvixiy95yMhKbM+7MdMb260CH1lbvb4Kgqgx2rXfeJShd5wxCt6vQ+ZyUA7e8X6fTWtI35jRUeX3MWLGNl78oIm/DbgAGZbXh8v6pjO6VQqvmNpqnqQfe6jqPJ2RJ35gg2VB6kHeXbuHfSzdTWHKQaE8E53VLZnTvFIblJNv4PqZRsKRvTJCpKss37+Xf32zhvW+3ULK/ggiBfh1bc17XtpzbtS09O7S0N36NKyzpG1OPvD5lWfEePssvYW7+DpYVOz01klrEcE5OEgMyEhiQkUCXdvHWEGwahCV9YxrQzgMVfL6mhDn5JXy1rpSdByoAiI+JpF96awZkJNA/PYEeKS1JjrcZoUzwWdI3xiWqyqZdZSzeuIu8ot0s3rCb/O37OfyfWmJcNN1S4unWviVd28fTrX08WUlxxNtUj+Y0BJr0A2omFpGRwJ9wZs56QVUfPWZ/DPAqMAAoBa5R1SL/vgeB2wAvcLeqzjiF6zCmyRER0hObk57YnCv6pwGwv7yK5cV7WbVtP/nb9rF6237+/vUGyqt8R45rExdNRmJzMto0Jz0xjow2zUlNaEb7lrG0bxVLbFTTGPjLNG4nTfoi4gEmAxcBxcAiEZmuqitrFLsN2K2qnf1z5D4GXCMiPYBxQE+gAzBLRLqoHh7NyJjwEB8bxdmdkzi789GBurw+ZUPpQdZs309R6SE2lB5kQ+khFhXt5t1lWzj2R3irZlG0bxlLu1axtI2PoU1ctLM0jybB/zmheRTxsVHEx0baTcLUKpAn/YFAgaoWAojINGAsUDPpjwV+6//8FvCsOFMVjQWmqWoFsF5ECvzn+yo44RvTdHkihE7JLeiU/J8jfVZUeyneXcbWPeVs21fO9n3lbNt79HPB9v2UHqykotpXy5kd0ZERtIyNPHITaB7toXl0JM2iPcTV+Bwb6SE2KoKYyAhiog5/9hDtiSAqMoIojzifjyyCJ0KI8kTgiRAiI8T/N4KICOe6IsTZ5hGx3kyNTCBJPxXYVGO9GBh0vDKqWi0ie4FE//YFxxybWudojQkTMZEespNbkF3LDaGmQ5XV7DpYeWTZc6iK/eVV7CuvZl95FfvKqtlfXsX+8mrKKr1s31dOWaWXQ5VeDlVWc6jSS7Wv/tv1nBuBU/UVIRAhzo1BcMYhO7z98F8QZ/vh/d9bP3oTcY71f0aObHPWD5ep/aYjx1mpyy0qWNNxdk9pyTPX9g/KuY6nUUwlLyLjgfEA6enpLkdjTNPRPDqS5tGRpCU0r/M5qr0+Kr0+yqt8VFR7qajyUV7tpapaqfT6qKqxVFYrVV4fXp9S7VO8Pp//r1LlVXw+xavO+uHPPp/iU/Cp81dV8ani9YGiqH+bUrMMcGRfjXJwpNpL0SND2By+bR3umHJ0vfZrrrm5ZmeWOt3+gnjP7JhQ/0N9BJL0NwMda6yn+bfVVqZYRCKBVjgNuoEci6pOBaaC03sn0OCNMacv0hNBpCeC5vZicViICKDMIiBHRLJEJBqnYXb6MWWmAzf7P18FzFbn9jkdGCciMSKSBeQAC4MTujHGmFN10id9fx39XcAMnC6bL6nqChGZCOSp6nTgReBv/obaXTg3Bvzl3sBp9K0G7rSeO8YY4x57OcsYY0JAoC9nBVK9Y4wxJkRY0jfGmDBiSd8YY8KIJX1jjAkjlvSNMSaMNLreOyJSAmw4jVMkATuDFE5TEW7XHG7XC3bN4eJ0rjlDVZNPVqjRJf3TJSJ5gXRbCiXhds3hdr1g1xwuGuKarXrHGGPCiCV9Y4wJI6GY9Ke6HYALwu2aw+16wa45XNT7NYdcnb4xxpjjC8UnfWOMMccRMklfREaKSL6IFIjIA27HU99E5CUR2SEi37kdS0MRkY4iMkdEVorIChH5udsx1TcRiRWRhSKyzH/Nv3M7poYgIh4R+UZE3nc7loYiIkUislxElopIvY06GRLVO/7J29dQY/J24NpjJm8PKSIyDDgAvKqqvdyOpyGISAqQoqpLRCQeWAxcHuL/OwsQp6oHRCQKmA/8XFUXnOTQJk1EJgC5QEtVvcTteBqCiBQBuapar+8mhMqT/pHJ21W1Ejg8eXvIUtXPceYuCBuqulVVl/g/7wdWEeJzLqvjgH81yr80/Se1ExCRNGAM8ILbsYSiUEn6tU3eHtLJINyJSCbQH/ja3Ujqn7+qYymwA5ipqqF+zZOA+wGf24E0MAU+EZHF/nnD60WoJH0TRkSkBfA28AtV3ed2PPVNVb2q2g9njumBIhKy1XkicgmwQ1UXux2LC4aq6hnAKOBOfxVu0IVK0g9oAnbT9Pnrtd8G/q6q/3I7noakqnuAOcBIt2OpR0OAy/z129OA80XkNXdDahiqutn/dwfwDk61ddCFStIPZPJ208T5GzVfBFap6pNux9MQRCRZRFr7PzfD6ayw2t2o6o+qPqiqaaqaifPf8WxVvcHlsOqdiMT5OycgInHACKBeeuaFRNJX1Wrg8OTtq4A3VHWFu1HVLxH5B/AV0FVEikXkNrdjagBDgBtxnv6W+pfRbgdVz1KAOSLyLc7DzUxVDZtujGGkHTBfRJYBC4EPVPXj+viikOiyaYwxJjAh8aRvjDEmMJb0jTEmjFjSN8aYMGJJ3xhjwoglfWOMCSOW9I0xJoxY0jfGmDBiSd8YY8LI/wd/jwuymt+eLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": {},
   "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.41421, 0.     , 3.60555])"
      ]
     },
     "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.13598, 0.15958, 0.0564 ])"
      ]
     },
     "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": {},
   "outputs": [],
   "source": [
    "def meanshift_inner(x, X, bandwidth):\n",
    "    # Find distance from point x to every other point in X\n",
    "    dist = distance(x, X)\n",
    "\n",
    "    # Use gaussian to turn into array of weights    \n",
    "    weight = gaussian(dist, bandwidth)\n",
    "    # Weighted sum (see next section for details)\n",
    "    return (weight[:,None]*X).sum(0) / weight.sum()\n",
    "    \n",
    "def meanshift_iter(X, bandwidth=2.5):\n",
    "    return np.array([meanshift_inner(x, X, bandwidth) for x in X])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "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+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9//HXJxsJSwgMIYawBMsSRUWuKaJWC4JLwQraxaX3Fqi31Na69mpdm6R1abU/FbX2Fq1oWxdsr6AVWisKxaqIoYiisinIkgBhIAmQPfn+/jgzk0kIJDBZmffz8chjMuecnPM9lb7nO5/zPd9jzjlEROTYF9PRDRARkfahwBcRiRIKfBGRKKHAFxGJEgp8EZEoocAXEYkSCnwRkSihwBcRiRIKfBGRKBHX0Q0I169fP5eZmdnRzRAR6VJWrly52zmX2tx2nSrwMzMzyc/P7+hmiIh0KWb2RUu2U0lHRCRKKPBFRKKEAl9EJEoo8EVEokTEgW9miWa2wsxWm9nHZpYXWD7UzN4zs41mNs/MEiJvroiIHK3W6OFXAuc650YDpwIXmtk44FfAQ865YcBe4KpWOJaIiByliAPfefYH3sYHfhxwLvCXwPJngGmRHktERI5eq9TwzSzWzD4AdgGvA58Bxc65msAm24CM1jiWSJe3ezc88ID3KtKOWiXwnXO1zrlTgYHAWCCrpX9rZrPMLN/M8ouKilqjOSKd29y5cMst3ms4fRBIG2vVO22dc8VmtgQ4A0gxs7hAL38gsP0QfzMHmAOQnZ2tJ6rLsW/mzIavQcEPggMHoEcPb32/fu3fPjlmtcYonVQzSwn8ngScB3wKLAG+GdhsOvBypMcS6ZIa99z79YObb/Zew9fNnAn33+9t09Q3AJEItUZJJx1YYmYfAu8DrzvnXgV+CtxkZhsBH/D7VjiWSNfTVAln3TqYMsUL+MbrLrwQJk+Giy9WmUdaVcQlHefch8CYJpZ/jlfPF4lOu3d7QX7xxd778BLOj38MixfDp59CTo63zfTpsGiRF/aLFsH48d62t9zivd58c7s2X449nWq2TJFjQjDoDxyAvDxYuhSeeQb69cO/fj2+J5+ETz7xtt20CQD/rFn4li0Dnw+uvdYL+/APiMb1fpGjoKkVRFpbsIRTVgbDhnm99fvvJ/e88zjltNNY/8ADUFDgbTt0KOv/8Q9OWbaMXAC/Hx591Av4YJknWO8XiZB6+CKtbeZMr3f/9tuwcSMAuU8+Sd7evQBMwBvRMCIujvWlpUzYtIkCIA9g6FBy77oLrrjCK/kUFUFqqkbsSKtQD1+ktQWDefFiOPNMcvv0CYU9QAFe6C+sqWGC309B2J/mbdrkBf7ixd6CVas0YkdajQJfpLWtWwfPPguA//PPeSIs7AGmMY0yUrgIKCOFaY1mHXni3XfxX3mlVw76yU+8kTyq4UsrUElHpLXddJNXyvH58O3YwRK8Hn0BXthfz/VMZSp55JFDDplkArCABQwAlhw4gO/1171yzqOPwsKFHXcuckxRD18kEk2Nk7/rLsjKgqlTARiBV7MfACxlKZvZTCaZzGUumWSymc0sZakX9sCIjAwv7IcNgwcfbP9zkmOWAl8kEsEROdOn14f+W2/B2rVeHb53b8AL/TlAMcXkeZdnQ/LIo5hi5gS240tf8lZkZnqvuvFKWokCXyQSM2fW3ygVvLB69tnQvTts2QIlJQCsB2YBKaSQQ06DXeSQQwopzApsx3/+p7fPxYu9femirbQSBb5IJPr1826qCr+w+otfeGPwARITWU99DX8840NlnJnMDJV3xjM+NHpn/bXXwqhRXkmnqMi7GauoSL18iZg513kmqMzOznb5+fkd3QyRyKxb502dMGQI/r//nVO2b28w9HIa01jKUoopJoUUxjOeBSwIrR8AfAj4br4Z5s8PjeXn/vs1vYI0ycxWOueym9tOPXyR1jZyJDz/PLz9Nr7t2/l+fHyD1QtYQHcr4VWgO8UNwh7g+4Bv0iSvlPOd73gLJ03S0EyJmAJfpC3MnetduB02jNzvfpechITQqgHAEueYMmkSS+bNY0DYuhwg9/jj4b//26vfjxvn9eyff97bQBdwJQIahy/SFoK98aIieOABcq+6Cl56iSdKS1lSW+uNxjnrLEZ8+9ssefhhJrz7Lt8Hbz6dzz+Ha67x5tW58UZvRk3wwl4zZ0oEFPgirW33bnjssYbLNmwgd+9err3qKnxffAFjxnh1/t27GZGQ4NXszzgD4uNh2TJvDP+yZTBxore/fv0O/aQskRZS4Iu0trlzvWmRwZvr/v77vZ7+smX4Bg6EJ59sOIXyP/+Jb/Jkb7RP8O+Ds2XecgsMGVI/Y6Z69hIBBb5IawvOlgleLz74KMPgrJdQH+bBD4Tw2TCDoa4evbQyDcsUaUvBnnzj6Y0PtVzkKGhYpkhn0NTzbKHhg8xF2olKOiJtSWUZ6UQU+CJtSRdapRNRSUdEJEoo8EVEooQCX0QkSijwRUSiRMSBb2aDzGyJmX1iZh+b2fWB5X3N7HUz2xB47RN5c0VE5Gi1Rg+/BviJc+5EYBxwjZmdCNwKvOGcGw68EXgvIiIdJOLAd84VOuf+Hfh9H/ApkAFMBQKTg/AMMC3SY4mIyNFr1Rq+mWUCY4D3gDTnXGFg1Q4g7RB/M8vM8s0sv6ioqDWbIyIiYVot8M2sJ/B/wA3OudLwdc6bsKfJSXucc3Occ9nOuezU1NTWao6IiDTSKoFvZvF4Yf+sc+6lwOKdZpYeWJ8O7GqNY4mIyNFpjVE6Bvwe+NQ592DYqleA6YHfpwMvR3osERE5eq0xl85ZwH8BH5nZB4FltwO/BF40s6uAL4Bvt8KxRETkKEUc+M65fwF2iNUTI92/iIi0Dt1pKyISJRT4IiJRQoEvIhIlFPgiIlFCgS8iEiUU+CIiUUKBLyISJRT4IiJRQoEvIhIlFPgiIlFCgS8iEiUU+CIiUUKBLyISJRT4IiJRQoEvIhIlFPgiIlFCgS8iEiUU+CIiUUKBLyISJRT4IiJRQoEvIhIlFPgiIlFCgS8iEiUU+CIiUaJVAt/MnjKzXWa2JmxZXzN73cw2BF77tMaxRETk6LRWD/9p4MJGy24F3nDODQfeCLwXEZEO0iqB75xbBuxptHgq8Ezg92eAaa1xLBEROTptWcNPc84VBn7fAaS14bFERKQZ7XLR1jnnANfUOjObZWb5ZpZfVFTUHs0REYlKbRn4O80sHSDwuqupjZxzc5xz2c657NTU1DZsjohIdGvLwH8FmB74fTrwchseS0REmtFawzKfB94FRprZNjO7CvglcJ6ZbQAmBd6LiEgHiWuNnTjnrjjEqomtsX8REYmc7rQVEYkSCnwRkSihwBcRiRIKfBGRKKHAFxGJEgp8EZEoocAXEYkSCnwRkSihwBcRiRIKfBGRKKHAFxGJEgp8EZEoocAXEYkSCnwRkSihwBcRiRIKfBGRKKHAFxGJEgp8EZEoocAXEYkSCnwRkSihwBcRiRIKfBGRKKHAFxGJEgp8EZEoocAXEYkSbR74Znahma0zs41mdmtbH09ERJrWpoFvZrHAb4CvAScCV5jZiW15TBERaVpb9/DHAhudc58756qAF4CpbXxMERFpQlsHfgawNez9tsAyERFpZx1+0dbMZplZvpnlFxUVdXRzRESOWW0d+NuBQWHvBwaWhTjn5jjnsp1z2ampqW3cHBGR6NXWgf8+MNzMhppZAnA58EobH1NERJoQ15Y7d87VmNmPgdeAWOAp59zHbXlMERFpWpsGPoBzbhGwqK2PIyIih9fhF21FRKR9KPBFRKKEAl9EJEoo8EVEooQCX0QkSijwRUSihAJfRCRKKPBFRKKEAl9EJEoo8EVEooQCX0QkSijwRUSihAJfRCRKKPBFRKKEAl9EJEoo8EVEooQCX0QkSijwRUSihAJfRCRKKPBFRKKEAl9EJEoo8EVEooQCX0QkSijwRUSihAJfRCRKRBT4ZvYtM/vYzOrMLLvRutvMbKOZrTOzCyJrpoiIRCouwr9fA1wK/C58oZmdCFwOjAIGAIvNbIRzrjbC44mIyFGKqIfvnPvUObeuiVVTgRecc5XOuU3ARmBsJMcSEZHItFUNPwPYGvZ+W2DZQcxslpnlm1l+UVFRGzVHRESaLemY2WLguCZW3eGceznSBjjn5gBzALKzs12k+xMRkaY1G/jOuUlHsd/twKCw9wMDy0REpIO0VUnnFeByM+tmZkOB4cCKNjqWiEiHKXU1vFS1k1JXc0TrWrK+tUU6LPMSM9sGnAEsNLPXAJxzHwMvAp8Afweu0QgdETkWLa7283RVAYur/fj9/ibXzd+xATg44MP/tj1EOkpnvnNuoHOum3MuzTl3Qdi6e5xzX3LOjXTO/S3ypoqIdLzGoX16XG+yY5N57e6HGHryKN5a+1Fo20nxPr66+QA3nDaBb9z5E16t2tUg4CfF+5iRMIBJ8b52aXuk4/BFRKJKsFde4eqopI5/Ve/hX/c+xvr75gBw4bmTuPO1vzBkxJcYtMlP3vnfZF/hTl6650G21VUw7PYfkGrxvFS1k0nxPi5NSGu3tivwRUSOQLA3XuFqmV+9i3X3/m8o7AHKCnfx8wu+wSmP3Mma6+6hrHBXaN2K+x6nmBoev/2HHMDxVtVeEi2GEXE9+EZCGsnWtpGsuXRERFqo1NWwuNrPpHgfFyX0p/+eCrY8Pb/BNtOYRmJhNSu+dT0JhVVMY1qD9Vuens8+/14APqOcj90B5lfvapc6vgJfRKSFguWchyu+oLCuAn/fRM5YOIfE9FTAC/vruZ6HeIhMMnmIh7ie60Ohn5ieyhkL5xDjSwEgk0TSiGdKXL92qeMr8EVEWqDU1VBSV0O6JZBfW8qvyzdTC/QcPiQU+ktZymY2k0kmc5lLJplsZjNLWRoK++ThQziDZC6PT2NcfAo7qSY1JqHNyzmgwBcRaVapq+Hhii+YX7OLQlfFQOtGz7BLoD2HD+GUR+6kmGLyyGvwt3nkUUwxpzxyJ72HD6EOGBLfnSu7DeC0uF4MtG6Miu3RLuehwBcRacbiaj/5taWkk8ComB5sc5UcH5uEBdbv3/AFH153NymkkENOg7/NIYcUUvjwurup2bAFgEpXy0tVO/lTZSHbXCXzqna2y3ko8EVEmjEp3kd2bDKFVBGHcUl8f9bU7cfhhf27U2ZRUVjEeMaHyjgzmRkq74xnPBWFRSye8n0yP/fTzWJ5uqqAGhyjYnowyLq1y922GpYpItKMZIvjqm4Z7CivZHXdfgAKXRV1/uJQ2AMsYAEAS1lKMcXcyI2MZ3xoeUVhEfO+9l1uWr2KT7v3ZHXdfkbH9GR+TRG9Y+LbfEy+evgiIs0odTX8vnI721wlA60bA2MSAejh68PgGZc02HYBC6hJT+CcPz9CRXp8KOyDhs/4Bv9OhuNjuwPea3vdbavAFxFpRngNf5urxDnHQOtGOY6Rt19N1m2zQtv2Tk/jywt/x5lfu4BJC58MDdkEGHHbLM6+/VqeriqgG8aMhAF8IyGNS9vhpitQSUdEpFmT4n1UuFo+qt1PYV0V210V21wlo6wHGEy+6w6es268P/cFrvzbH/ENz2R13X5OGpnFgEV/5E+Tv8vAGdPIvv0ahsQkcnJMLy5KSG2XkA+nwBcRaUayxZFosXxcd4DRMT05PjaJE2K789X4vvy+cjuLav3MvOsWhv3gcrb1TeIrsT2Jtxjya0uZceKpXLv63zzevYRtrpJFtX5mxA5o97CHKAn8KlfO1uq1DIrPIsGSOro5ItKFBKdTOD2uNwAVro4XqncwsiiNuv4l5LtSRsf04vO6cnb3TWJ0TE8AruqWwUk1PZkU7yM5NY5fujRerSoCXLvNjtlYVAT+1uq1fFq1HIAvJYzp4NaISFcSnE4B4NKENEpdDasKHPOW9CDjtHhmnDCAClfLC9U7GWjdOD62Oy9U7yDRYhqMukm2OK7slt5RpwFEyUXbQfFZnJAwjkHxWS3avsqV81nVKqpceRu3TEQ6kr/c8fiqKvzlBz9OOzjv/elxvZmRMIDT43rzUuAGqVuPS+fWMUn814gkRsX24J81xWRZd7a5ytDF2I7qxR9OVPTwEyypxT37KlfOBxVvsqt2C/7aAk5NPFdlIJFj1Ly11dy9vAqAH41JCC0PTqWQX1vKmtr93JA4pGFPPykttP1dB7ZQ6CqBhFDQd0R9viU6Z6s6QLDOX+Nq2FW7hR6Wwq7aLXxQ8SanJp4LoOsAIseYy7LiG7wGBYdhDrRu5NeWhqZEBg7quV+fOJjZFVu4PnEwI2N7hr4ZdMbg71ytaQPhF2yrXAWfVL7DsITT2FtbyKD4LPx+P2XJRdS4GjZU5zM8/jROSBhHXEkyO3uuZVftFrZWr6WyrpzPa1ZTVruP7rG9FPwix4C9FY53C2o5PzMOX5KFlgdD/fS43rxXU8Lpcb1Dod84xEfG9uTxHieG3jeu+Xcmx3wNP3jBdmv1Wj6pfIddtVtYVbGYT6uWc/WdMzjllNEs/vhlwHFCwjiGJpxM7eYenHnqOfz1l+8wPD6bGldNcZ331Jqium2h/YlI1xSs3d/xVgVvbKkl753KBuuTLY5LE9LIiEnk0oQ03qsp4emqAl6tKmrwPNumtPdzao/EMR34Va6cGlfD8XGj2Vezl9LaPSRbKnWulufufZW5977I7sI93DFlNis/XUFlXTkfr1vDORO+QkFBAXf//F5+c/eTbKheSUpMf/rHDubUxHObvACsC70iXUewdl9d5zg7I4acM7s1WN/4Ym4wxMHxdFUBv9i8q8kLvVD/YdHZyjlwjAf+1uq1bKjOZ7/by7a6dVSwn33Oz9x7X+T5+xaFtvMXFnPt5DuZ9+qfOO/cC9hZUBRa9+jdc/jzva+THOvdHp1g3fhSwpiDyjnh3yREpHO7LCueszNiWF7oODXD8WGPoga99uAHwry11UB9iF+U0J+RRWnMW9IjtK4r6XwfQUepqZurgr3wtLhMEiqT2Fa7jhJ/Ka89/XaDv53GNJYWLuXn3/otKaQwjWkNJjx69ek3+doP/kpPXxJV5ZXEWRzJMf0Y1m3MQcdq6dBPEek4viTjy+mxvLW9jh0pxbxb5ZVsgzX38Iu5/nLHvLXVXJYVjy8pjluPS2fomOqDLvR2BcdMD7+5HnYVXqkl2deTexfeQN9076655p5B2Te9N/cuvIGePi/Y97k97K7bzuc1qxscK8GSGBSfxdbqtSrriHQBM09K4KbT4jmuJIVvW8Oauy/J+NGYBHxJ1qC33zD87bDj+DujiALfzB4ws7Vm9qGZzTezlLB1t5nZRjNbZ2YXRN7Uw2vq5qrgh8CairfYVbsFIxaAjOFpodA/3DMo+6b35v6FtzJk+MDQPmupJpEeHB83OnSsYP1+U9VHKuuIdBG+JKN7vPHoe47S9X2orojl1ysq+fX7lQ0C/LKseO4cl8BlWfEHlXoav+/szLmj/2Qys/OBN51zNWb2KwDn3E/N7ETgeWAsMABYDIxwztUebn/Z2dkuPz//qNvTWLDMExxSCQbUn+/7f/+In3/rt6GwDwo+qeZnf/4hX77w5IP2292S+Ur3S0PlnM+qVvFp1XKGx2cTZ3EasinSBWzcW8eNb5ZRVG48NimRFYW1oZuwxqUbY/rHkhRvzDwpITRks6kefvj7jmJmK51z2c1tF1EP3zn3D+dCVzqWA8Gu8FTgBedcpXNuE7ARL/zbVbDMUkUlMcQSHvbbN+zkseueO+wzKB+77jm2bzj4WZNlrrRBLz747WJowklNXtAVkc4n751KVu6CLfscD7xfSVm1Y1y6F4nLCx2/XV3Dg/nVDXrvwVIPwOOr6u/Q7ciwPxKtWcP/HvC3wO8ZwNawddsCy9rd1uq1bKtZSx31Xy52bPCGYu4pLDnsMyj3FJZw+5SHKdhQP2qnJz5OSBhHWlxmaBhmcOoGBb1I5xVeb9+4t44D1Y5RfWFwLyOjBzy4spox/WO4KTuey0fGMrgXXD4ylrJqd1CNvquVcoKaHaVjZouB45pYdYdz7uXANncANcCzR9oAM5sFzAIYPHjwkf75QRqP1hkUn0VhzSaK67yeerm/mp9OeYA9hSVA88+g3FNYwh1TZjP73Vvp6/NxStI57K0tZHv1BjZUrwQ0A6dIZxYsu5RVOx5cWc22fbW8uK6Wshro0w32VjoGlAZ76Eb3OGP7fseWfZCws46N62rpHm8N5tq5LCueshoX+jA4Znr4zrlJzrmTmvgJhv0M4CLgO67+gsB2YFDYbgYGljW1/znOuWznXHZqampTmxyRxqN1EiyJsUlfY0jcKLpbMsNTT+KCGWc1+JsFLCAm3fGzP/+QmHR30DMovznzYo7rl8HYpMnsrS0MTLVsRzQDp4h0jLlrqrh7eRV7KrybrP7ywS7KarwrensDN9hm9TW+klJMeU0ddy+v4kspxsTBsTw4ITF0wTacL8n7YHhwZXWX6uVHNA7fzC4EbgG+6pwrC1v1CvCcmT2Id9F2OLAikmO11KD4LGpcDTWuOlRuSbAkusf0oqymlIy4EVx359UAoZuv+qX35e6F1zJs+DAyFqZz25T/F/oG8J3bLmbKbWMpcyWh+XeCx2lcwtGDVkQ6F3+54/1Cr5z75pZaPn7+Hsrefpp+NywkNm04GT3gshPi2frZBv7vh+fzpQtmwPjbAeOMATEM7R1D9nGxTe77UBOvdWaR3nj1GNANeN3MAJY75652zn1sZi8Cn+CVeq5pboROa0mwJOIsjk+rlhNn8aFyS3hQD+UkcnKOo2/MAObPXchrb/6d/YM2s7tuO2NHnsWDCxP5nyn38bUZ4/n27efR3ZLJiBsRCvJDlXD0oBWRjtV41My8tdW8tb2Oob2ND5+9h/2L7gNg98NT8N2wkEGnjuTsbpv56jXnU1dSyIYX76Pnfniz1x1s2eeNRwkv5YQLv4DbVUQ0LLO1tdawzCPpafv9fnw+H+sq3w/MlpnNyG5fxu/3YynVrK5YwujECfSM6R3aJzQ9VbJ6+CJto7nhj43r9Ddlx9M9zsjqG0PuO5VU//1elv/h3gZ/k9g3nftm/y+3XX81FXsKG6zrOfk2Jl99Jyf5Yg4amtkZtXRY5jEztUK4I3ngic/n3V03NOEk4iwuNPomrU8mn1Su5IArZm9tYVjt3tNUT/5IjisiLXeoB5U0Xn9Tdjw3ZcfzfmEtb22vY3Av2Fzop+TVuQ22n8Y0lu5Zyo3/NZUUUriw0XQqNcufZuNFP+CtbV4+dI+z0I1XHT3mPhLHZOAfjWBYB2+i8tcWsKt2C/1jBze4MHuo30Wk7TRXLz8/M453C2qZNiyeBRu8Ms5pacY6vyOmp48v5yzi/bzJVOwpDE2nMpWp5JFHDjlkkgl4AzjiUtLpde1CCvHRMw5OSo1hbHosN7zpTaUMhy7zdHbHZEknEsGyTFpcJjtrNqs8I9JJhZd5gj38YSlGdloML6yrJTkBSr0vBQzuBZ9v3ID/4Skkl5SH5s4K2sxmbuRGSnsn4bthIXFpwxsc6+yMGN7aXsfZGTE8fl5Sp+vht8udtseiYE+/Z0wf3Uwl0pHKS2HVfPbsLWlygrLwm58uy4pnWIqxsdjxXmEdUB/2Kd1gyz6ITxtO7ysfoZhi8shrsK888iimmN5XPtIg7EenwtkDY8jo5UXlSf1iO13YHwkFvoh0TmvfgOXPsPad1w++q7W8lOk1r3L3aRWhmvpTFyYxcXAssycmcnaGF21Dk43iSvAlQvXODZQ8d91hp1Mpee46anZuAKB3AozpH8db2+rYvs/7EEkKFMG72iyZQQr8FtDTrETa3kEhmjURxk0n68zzDr75ae0b9Mj/A9+LfyvU4x7WJ4Y/Tkki+7hYHj8viTvHJTB7YiLDUoydX2yg5JEp1JUUHnY6lbqSQvwPTyFp7wZKquCj3bUMSzF+MDqBO8clMPNkr3Z/zE6tIBpfL9IeDhqJk5QMYy6hL/CjPo02zhwLBWu81yYEx8g/vqqK9dt2s//RKVTt9YZeNjedSl1JIbsemsLEh95jR1kftuyDuWuq+eOU+vJuV7zpCtTDb5Gm5tqXyJTthrcf8F4ligXq9JSXNph3/lD85Y5fv1/JP996B7ashM0rGnwzaPwt4bKseP5n/HEMPX9mg/0sYAGlvZPo88M/U9o76aDpVH446/ucPyqNLftgWIod9Mzb8AekdCXq4beAxte3vlVzYfEt3u9n3dyxbZFGyku9+nnmWNi8wiutJCW3zbECdXoA35hLmh3uOG9tNQ/mV9OXrzBnsHFm1sTQN4OyasfqorrQ0Mng6B0c7D3nNk4EPnnBu/kqpnc6Q25eSEXf4Rx320K2/noK+3d73wBycnLIzc1t8KFxqPnwuxoFvnSIMTMbvkonEgzhtYuheDvUVEBc4uE/APZuh3eegjO/B4m9vH0c6oMi+IGSNdH7gfrXZgRnqcT1ZeTJ34Ak47IsL5jLahxvbKll4uDYBkM1b8oOfHOYcTe/Oj6O3z3xJFc//hqx/Yfx9Mc1XHH2idSM/Ad3zzifCy+7itzcXKDpqROauwGss9M4fBFpaO92+Pu9Xtgnp0O3HlC0EXqk4t+5Hd+AwfCVWbD6Zeg3FE44D/+ff4av1g+DT4P+wyB/Hoya7O2j31AYc6kX/uWl8OZsrxwzbjqMueTwbQn/cGjmW0ZzT6MKvj8vtZThA/s1WA/w++U7uGrccYftuTf+m87S29c4fGkXjWvxfr+/ye0OtVw6mfJSr6devB1SMqC00At7IHfeW5xyzxLWb9gIi+6G7ath9QLW//ZqTrnlOXIXbYARE2DbGgBK1q8IbcObs+vDe8tK74Ohca8+rJ4fEvq28UZoUUuHRIbX2f3ljhverODu5VW8XNjroKdV+ZKMWyakNxvch3q4eVehko60WNlur/Y+ZiZ07+cte+8xWJYHVQfgn+TyxBNPsGTJEgb2HcF7j0FNGewoX89P/zKBWT/4fujrsnRSHy30AjljNIy9Et79A/g/J3f+KvIWrQNgwsNvs+SGsxiR1pP1O/cz4eG3KSipIO/VTyDuWnLPHwLAZ9W9+A92Q49Ub59vzvahfGLNAAALd0lEQVRKPtB0jz0Q7iu2HmBUend6nDSp/kMhc6z3YZA1kXlrExuUVRpPnBZcHm7e2upQuQdHq5RluuJIHQW+tFhTF1qD/aFnlubyzD+9uxdPP3kCP+q/hIRtI9jNep5hAvsoIC8vjy1vw/2P5bLulYYfHNLBykthzUJvqGPQ1n/Djo/JfXVtKOwBCkoqmPDw28y5cjSznltNQUlFaF3egg+gqpw7Lz+HuB4jYOcmGHYW7N3qhT7AudfXh314ySZzLJvXfkj+9jLGbv8LxJlX8hlziRf2gYu7l2VNC7zWl1Ua1OqbCODG4dw93iIOak2PHCHV8Duvst2w4jHvMfCn/9gL6rLd8K/74Tcv5vLKFw1vVe/FAC7tPoeXymaxj4IG6y5IyeGM4lyGToLBZ8FJV6APgI4WFqh06wWV+2DkRPwfLeOUnFcbhPo0poXGsKeQ0mAMO8CA3ol8eMcEfKdPhV79vd75p6/Dpve8ElF47T543HHTvffLn2FFxjfre/hNfTA0+mbQ1UfOtAbV8KVVrZoL/8zzevSr5sLudTB/OrzxgJ8lXzzRYNtpTCOWMp4pu4hYypjGtAbr3yl+gqokP5sWe/t87Sbvm8OqhjPYSnvKmgjZl3mlnMp93rJtH+JLqGbJDWcxoHciQGimyeDkYw/xENdzfei/8YDeiSy54Sx8PROgpMAL9o3LvDp+aaG3/5qK+jp94G7a0Iid7MsYe1xsw7CH0E1YTV247apj4juCAl9aZMxMmHQ/VJd54fy3H8PGRdC7l4/pLKEXA4DmA6EXA5jOEhLKfXRLgTNuhnPuguMnwYEi3YjVYZKS4ctXwHk/8UbmAFAL8T0YkdYzFPpLWRqaimAuc0NTFCxlaSjsR6T1hIQe8OUrvd0EiwgZoyE9yxvBE7wIGx7kScne8M+V8xpcpJXWo8CXZoVfrI3r7i3rOxJ8WVC9D/oxIhT6hwuEYNj3YwQAlcWwYxV89nf4fDG8+4B6+R0qWDaZeIM3QufAHqg+AMCItJ7MuXL0YWeanHPlaEak94H+I72r+Fv/7ZVsBv+HF/b9hsKwc+p79E2Nygnv8UurU+BLyKGmOwherF01F06+AoZOgqKPwL8W+p8KFuuF/teZc9hA+DpzQmEftGmx1wH8ag6ck6MbsTpUcAhk4Scw7T5IHxVatX7nfmY9t/qwM03Oem416wv3QmwcDBwN1ZXe/la+WD88c/OK+h79moXe+jUL63d2mNKNRE6BLyHhwR5uzEwvkEu+gOcu8kL6i2UwbDJkTQVXC7tZz1+ZddhA+Cuz2M36BusGf9X7EInvUX8xWDpIsHedOdYL4UApJnzo5eFmmiwoqWDC7HdZ/8F7sG01+DfBaZd5QzFHT/N6+eGTnblGr9LmFPgSEqzTN+5ld+/n/X/y/d/A3o2Q5KuvvW97G2JOqh96ebhA2EcBzzCB3awn1rsGSMZYb4SOLtp2AknJXui/85RXZ9/xMf7k4Ux4dEVolM4CFjCb2dzIjaGnRM1mdmiUTkFxGRMeXYE/ebgX+vGJ0CcDknp7vfzNK+qPd/IU7wPm5CkdcbZRSePwJaR7v0NPZBYc/5Dkg3I/JHSHl2fC1rV+fhc7ITT0srmpZ4Oh/8OKD+mOj4TumlenUwneCZsxGo7LwnfyFL7/SRx59/4ytMkCFjAgJYk/XXE6s/68jgW7G840+f1Lz8M3+Yb6eXeg6TlzguUbaTcahy8tErxwO/Jir0de8oXX408eAktqc/nrtoPH4X+dOfyVg8fhf31gDqdty6XvMLjqXZVxOpVDjHfPveO2UOgP6NODJY/ewojSf7M+9Rwm/CCHgj37Acj59pnkjk/1LvpeeLvXu5c2p3H40qqCvf9+I73X3YEbL2Pj4bRtuVwyrL5un2zeaJxT+kzhp6cuISV+QGjdzdfk8NTiXIZNhiteVdh3Ooe4aJr7zbHkTB7phf114xjRqxZGT2PE4HSWPHIzA3onknPxyV7YJyZ7c/G881QHnYQcinr4clR2r/NumDrnLtjylleOuf+xXH7z0BNcXrqEkVkjuHyB9wGx+r31TDrfm0vnnvtzO7rp0lLhvX2AtW/gTx6Ob+38+gnQtqyEUZPxr1yELxGvZz/hOm9kzpnfUw+/nbS0hx9R4JvZL4CpQB2wC5jhnCswMwNmA5OBssDyfze3PwV+17d1vZ8tL/sOmibB7/fj8/k6rmFy5MKnPQivtTd+QMrWD7wLsonJ3nBOhXy7a6/AT3bOlQZ+vw440Tl3tZlNBq7FC/zTgdnOudOb258CX6QTaelc9OEPP1HYd4iWBn5Eo3SCYR/Qg/oRtVOBPzjv02S5maWYWbpzrjCS44lIO2rpKJo+GTDlrrZvj0Qs4mGZZnYP8F2gBJgQWJwBbA3bbFtg2UGBb2azgFkAgwcPjrQ5IiJyCM2O0jGzxWa2pomfqQDOuTucc4OAZ4EfH2kDnHNznHPZzrns1NTUIz8DERFpkWZ7+M65SS3c17PAIiAH2A4MCls3MLBMREQ6SETj8M1seNjbqcDawO+vAN81zzigRPV7EZGOFWkN/5dmNhJvWOYXwNWB5YvwRuhsxBuWqZvmRUQ6WKSjdL5xiOUOuCaSfYuISOvS1AoiIlGiU02tYGZFeKWhzqYfcKw8fE/n0jkdK+dyrJwHdK1zGeKca3aYY6cK/M7KzPJbchdbV6Bz6ZyOlXM5Vs4Djq1zCVJJR0QkSijwRUSihAK/ZeZ0dANakc6lczpWzuVYOQ84ts4FUA1fRCRqqIcvIhIlFPiHYWYPmNlaM/vQzOabWUrYutvMbKOZrTOzCzqynS1hZt8ys4/NrM7Mshut62rncmGgrRvN7NaObs+RMLOnzGyXma0JW9bXzF43sw2B1z4d2caWMrNBZrbEzD4J/Nu6PrC8y52PmSWa2QozWx04l7zA8qFm9l7g39o8M0vo6LZGQoF/eK8DJznnTgHWA7cBmNmJwOXAKOBC4HEzi+2wVrbMGuBSYFn4wq52LoG2/Qb4GnAicEXgHLqKp/H+dw53K/CGc2448EbgfVdQA/zEOXciMA64JvDfoiueTyVwrnNuNHAqcGFgHrBfAQ8554YBe4GrOrCNEVPgH4Zz7h/OuZrA2+V4s36CN1HcC865SufcJrw5g8Z2RBtbyjn3qXNuXROrutq5jAU2Ouc+d85VAS/gnUOX4JxbBuxptHgq8Ezg92eAae3aqKPknCsMPrrUObcP+BTvuRdd7nycZ3/gbXzgxwHnAn8JLO8S53I4CvyW+x7wt8Dvh3rAS1fU1c6lq7W3JdLCZpPdAaR1ZGOOhpllAmOA9+ii52NmsWb2Ad7zuV8HPgOKwzp9Xf7fWsRPvOrqzGwxcFwTq+5wzr0c2OYOvK+vz7Zn245US85FOjfnnDOzLjV0zsx6Av8H3OCcKzWz0LqudD7OuVrg1MC1uvlAVgc3qdVFfeA394AXM5sBXARMdPVjWDvlA16O4GE14TrluRxGV2tvS+wMPvPZzNLxephdgpnF44X9s865lwKLu+z5ADjnis1sCXAGkGJmcYFefpf/t6aSzmGY2YXALcDFzrmysFWvAJebWTczGwoMB1Z0RBtbQVc7l/eB4YHREwl4F5xf6eA2ReoVYHrg9+lAl/g2Zl5X/vfAp865B8NWdbnzMbPU4Cg8M0sCzsO7JrEE+GZgsy5xLoflnNPPIX7wLmBuBT4I/Pxv2Lo78Gp864CvdXRbW3Aul+DVICuBncBrXfhcJuONmvoMr1zV4W06grY/DxQC1YH/HlcBPrzRLBuAxUDfjm5nC8/lK3gXNj8M+//I5K54PsApwKrAuawBfhZYfjxeB2gj8GegW0e3NZIf3WkrIhIlVNIREYkSCnwRkSihwBcRiRIKfBGRKKHAFxGJEgp8EZEoocAXEYkSCnwRkSjx/wHkPo84yNCNPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": {},
   "outputs": [],
   "source": [
    "def meanshift(X, it=0, max_it=5, bandwidth=2.5, eps=0.000001):\n",
    "    # perform meanshift once\n",
    "    new_X = meanshift_iter(X, bandwidth=bandwidth)\n",
    "    # if we're above the max number of allowed iters\n",
    "    # or if our approximations have converged\n",
    "    if it >= max_it or abs(X-new_X).sum()/abs(X.sum()) < eps:\n",
    "        return new_X\n",
    "    else:\n",
    "        return meanshift(new_X, it+1, max_it, bandwidth, eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.09 s, sys: 0 ns, total: 1.09 s\n",
      "Wall time: 1.09 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": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE5VJREFUeJzt3W9sZFd5x/Hvk2xIWdHKYLYhk0R12satlsoNwopA7Ys1f0poVlkXQRVeVCmNbIqoFAISDUWq7ReVQEhNkVoCdiHNC9SUUli2G/4lkVeoUgs4NGwTwjpLWJTEgTim2xJZCk3y9IWvnfGud73rmfF47vl+pJHvnHN97znyzE/XZ86cG5mJJKn+Luh2AyRJ28PAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBViV6sHiIhfAL4BXFwd7/OZORERVwJ3Af3A/cAfZebPz3asV77ylTkwMNBqkySpKPfff//Tmblns/1aDnzgWeANmflMRFwE/FtEfAV4P3BbZt4VEZ8EbgJuP9uBBgYGmJuba0OTJKkcEfGjc9mv5SGdXPFM9fSi6pHAG4DPV+V3AqOtnkuStHVtGcOPiAsj4gHgKeAe4AfAycx8rtrlceCydpxLqpOlpaXzKpda0ZbAz8znM/Nq4HLgGuA3z/V3I2I8IuYiYm5xcbEdzZF6wuTkJENDQ8zPz68rn5+fZ2hoiMnJye40TLXV1lk6mXkSmAVeD/RFxOpnBJcDT5zhd6Yzczgzh/fs2fQzB6kWJicnmZqaYmFhgZGRkbXQn5+fZ2RkhIWFBaampgx9tVXLgR8ReyKir9p+KfBm4GFWgv/t1W43Al9q9VxSHayG/arV0L/77rvXwn6Voa92ascV/qXAbEQcBb4N3JOZh4E/B94fEcdZmZr56TacS+ppS0tLzMzMrCsbZZTlhWX279/P8sIyo6fMb5iZmXFMX23R8rTMzDwKvGaD8kdZGc+XVOnv72d2dnbtSn6UUW7mZg5wgCmmmGCCAQYAOMhBGo0Gs7Oz9Pf3d7fhqgW/aStts8HBQWZnZ2k0GhzhCCc4wQAD3MEdDDDACU5whCNrYT84ONjtJqsmDHypCwYHB5menuYkJ5lial3dFFOc5CTT09OGvdrKwJe6YH5+nvHxcfroY4KJdXUTTNBHH+Pj46dN2ZRaYeBL26x56uU+9q0N47yLd60N7+xj32lTNqVWRWZ2uw1rhoeH07V0VGdLS0sMDQ2tm3o5yihHOMJJTtJHH/vYx0EOrtU3Gg2OHj3qB7c6o4i4PzOHN9vPK3xpG/X39zM2Nrau7CAH2d3YzeHDh9nd2L0u7AHGxsYMe7WFgS9ts8nJSSYmXhy3X52Nc911163N3lk1MTHhF6/UNu1YHlnSRiJe3D5l6HQ1xGdmZtZNvVydsjkyMsLY2Jhhr7ZyDF/qlLME/qqlpaUNh2vOVC5txDF8qQecKdQNe3WCQzpSp+yg/54l8Apfkoph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klSIlgM/Iq6IiNmI+F5EPBQRN1flr4iIeyLikerny1tvriRpq9pxhf8c8IHM3Au8DnhvROwFbgXuy8yrgPuq55KkLmk58DPzycz8TrX9M+Bh4DLgAHBntdudwGir55IkbV1bx/AjYgB4DfBN4JLMfLKq+jFwSTvPJUk6P20L/Ih4GfAvwPsy83+b6zIzgQ3v6BwR4xExFxFzi4uL7WqOJOkUbQn8iLiIlbD/bGZ+oSr+SURcWtVfCjy10e9m5nRmDmfm8J49e9rRHEnSBtoxSyeATwMPZ+ZfN1UdAm6stm8EvtTquSRJW7erDcf4HeCPgP+KiAeqsr8APgJ8LiJuAn4E/GEbziVJ2qKWAz8z/w2IM1S/sdXjS5Law2/aSlIhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFaEvgR8RnIuKpiHiwqewVEXFPRDxS/Xx5O84lSdqadl3h/wNw7SlltwL3ZeZVwH3Vc0lSl7Ql8DPzG8BPTyk+ANxZbd8JjLbjXJKkrenkGP4lmflktf1j4JIOnkuStIlt+dA2MxPIjeoiYjwi5iJibnFxcTuaI0lF6mTg/yQiLgWofj610U6ZOZ2Zw5k5vGfPng42R5LK1snAPwTcWG3fCHypg+eSJG2iXdMy/xH4d+A3IuLxiLgJ+Ajw5oh4BHhT9VyS1CW72nGQzHznGare2I7jS5Ja5zdtJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCdDzwI+LaiDgWEccj4tZOn0+StLGOBn5EXAj8HfBWYC/wzojY28lzSpI21ukr/GuA45n5aGb+HLgLONDhc0qSNtDpwL8MeKzp+eNV2ZqIGI+IuYiYW1xc7HBzJKlcXf/QNjOnM3M4M4f37NnT7eZIUm11OvCfAK5oen55VSZJ2madDvxvA1dFxJUR8RLgBuBQh88pSdrArk4ePDOfi4g/A74GXAh8JjMf6uQ5JUkb62jgA2Tml4Evd/o8kqSz6/qHtpKk7WHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhWgr8iHhHRDwUES9ExPApdR+KiOMRcSwi3tJaMyVJrdrV4u8/CLwN+FRzYUTsBW4AXg00gHsjYjAzn2/xfJKkLWrpCj8zH87MYxtUHQDuysxnM/OHwHHgmlbOJUl1s7S0dF7lrerUGP5lwGNNzx+vyiRJwOTkJENDQ8zPz68rn5+fZ2hoiMnJybafc9PAj4h7I+LBDR4H2tGAiBiPiLmImFtcXGzHISVpR5ucnGRqaoqFhQVGRkbWQn9+fp6RkREWFhaYmppqe+hvOoafmW/awnGfAK5oen55VbbR8aeBaYDh4eHcwrkkqWeshv2q1dCfnp5mfHychYWFtbrV/doV/J0a0jkE3BARF0fElcBVwLc6dC5J6glLS0vMzMysKxtllOWFZfbv38/ywjKjjK6rn5mZaduYfqvTMv8gIh4HXg/cHRFfA8jMh4DPAd8Dvgq81xk6kkrX39/P7OwsjUYDWAn7m7mZ27iNAQa4jdu4mZvXQr/RaDA7O0t/f39bzh+ZO2cUZXh4OOfm5rrdDEnqqNWx+uWF5bWwX3WCE9zCLexu7GZ2dpbBwcFNjxcR92fm8Gb7+U1bSdpmg4ODTE9Pc5KTTDG1rm6KKU5ykunp6XMK+/Nh4EvSNpufn2d8fJw++phgYl3dBBP00cf4+PhpUzZbZeBL0jZqnnq5j30MMMAJTvAu3sUJTjDAAPvYd9qUzXZwDF+StsnS0hJDQ0Prpl6OMsoRjnCSk/TRxz72cZCDa/WNRoOjR4+e9YNbx/AlaYfp7+9nbGxsXdlBDrK7sZvDhw+zu7F7XdgDjI2NtW2WjoEvSdtocnKSiYkXx+1Xp15ed91166ZsAkxMTLT127atrpYpSTpPqyE+MzOzburl4OAgs7OzjIyMMDY21valFRzDl6Q2uv7H/8kLu+GCZTj0qtecdd+lpaUNh2vOVH4mjuFLUhe8sBsuuGDl52bOFOrtGrM/lYEvSW10wTK88MLKz53GMXxJaqO1YZxf6m47NuIVviQVwsCXpEIY+JJUiGICf7tvFixJO00Rgd+NmwVL0k5T+8Dv1s2CJWmnqXXgn+lmwXffffda2K8y9CXVXW0Dv9s3C5aknaa2gd/tmwVL0k5T28CHF1eeazQaHOHI2t1k7uCOtbvMHOHIWti3+/6RkrST1DrwoXs3C5aknab2gd+tmwVL0k5T68Dv5s2CJWmnqe0NUDp1s2BJ2mmKvwFKt28WLEk7TW0DH7p7s2BJva9ua3C1FPgR8bGI+H5EHI2IL0ZEX1PdhyLieEQci4i3tN7Uszv8zO1rj2aroX/q1MvmKZuGvaRT1XINrszc8gP4PWBXtf1R4KPV9l7gu8DFwJXAD4ALNzvea1/72tyqf/3ZJ9YeG3n66afPq1xSuSYmJhJIIBuNRh47diwzM48dO5aNRmOtbmJiorsNrQBzeQ6Z3dIVfmZ+PTOfq57+B3B5tX0AuCszn83MHwLHgWtaOVertvtmwZJ6U53X4GrnGP6fAF+pti8DHmuqe7wq65j9L3vP2kOStqLua3BtGvgRcW9EPLjB40DTPh8GngM+e74NiIjxiJiLiLnFxcXz/XVJapu6r8G1a7MdMvNNZ6uPiD8G9gNvrMaSAJ4Armja7fKqbKPjTwPTsDIPf/MmS1LnrE7oGBkZ4cjCEQ5wYG0NLqCn1+BqdZbOtcAHgeszc7mp6hBwQ0RcHBFXAlcB32rlXJK0Xeq6BlerY/h/C/wicE9EPBARnwTIzIeAzwHfA74KvDczn2/xXJK0Leq6Blers3R+PTOvyMyrq8efNtX9VWb+Wmb+RmZ+5WzHkaSdos5rcNV2LR1JOl+9ugZX8WvpSNLZXD9zP4/+/TjXz9y/Vlb3NbgMfElF+pv4FL/6f0/xN/GpdeV1XoPLwJdUpPflu3n0ol/mffnu0+rqugaXY/iS2mZpaWnD4Y0zle90vdIfx/DVFXVbTlbnro6rS9ZtDS4DX21Txze8zs3qgmOnTlVsnuLYawuN1dK5LKm5XY9WlkdWd/XacrJqn+a/ffNr4PDhw+v+9r4GOodzXB656yHf/DDwe5Nv+HI9/fTTp/2NRxnNPvoSyD76cpTR014b3oeivc418B3SUUvqvpyszq7uq0vWjYGvlviGV/NUxSMcWVt+4A7uWFuWoFdXl6wbA18t8w2vuq4uWTcGvtrCN3zZ6rq6ZN0Y+GoL3/DlqvPqknVj4KtlvuHLtbS0tO7G3gc5yMf5OLdwCyc4wS3cwsf5+NqCY6uvAT+07w4DXy3xDV+2/t0XMbb/d9eV1Wl1ybox8NWSui8nq/WeuX2UF24f5Znbq6m237+PyaufZWL8D9f2qdPqkrVzLpP1t+vhF696l9+0LcPznziQ+YkDKz8zM5f/J/M7X8hc/p+cmJhY97dftfoa8G/fOZzjF69cLVNtMzk5yczMzGlTL1fH+MfGxry663HP3D7KbmAZeNl7Dp5W3yurS9bNua6WaeDrnE3Fi9sTZ3jZ+IaXtp/LI6sr6racrFQnBr4kFWJXtxug3nGmYRxJvcErfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCrGjllaIiEXgRx049CuBpztw3G6zX73FfvWWXurXr2Tmns122lGB3ykRMXcu60z0GvvVW+xXb6ljvxzSkaRCGPiSVIhSAn+62w3oEPvVW+xXb6ldv4oYw5cklXOFL0nFq3XgR8THIuL7EXE0Ir4YEX1NdR+KiOMRcSwi3tLNdp6viHhHRDwUES9ExPApdT3bL4CIuLZq+/GIuLXb7dmqiPhMRDwVEQ82lb0iIu6JiEeqny/vZhu3IiKuiIjZiPhe9Rq8uSrv6b5FxC9ExLci4rtVv6aq8isj4pvV6/GfIuIl3W5rK2od+MA9wG9l5hAwD3wIICL2AjcArwauBT4RERd2rZXn70HgbcA3mgt7vV9VW/8OeCuwF3hn1ade9A+s/A2a3Qrcl5lXAfdVz3vNc8AHMnMv8DrgvdXfqNf79izwhsz8beBq4NqIeB3wUeC2zPx14L+Bm7rYxpbVOvAz8+uZ+Vz19D+Ay6vtA8BdmflsZv4QOA5c0402bkVmPpyZxzao6ul+sdLW45n5aGb+HLiLlT71nMz8BvDTU4oPAHdW23cCo9vaqDbIzCcz8zvV9s+Ah4HL6PG+5YpnqqcXVY8E3gB8virvuX6dqtaBf4o/Ab5SbV8GPNZU93hV1ut6vV+93v7NXJKZT1bbPwYu6WZjWhURA8BrgG9Sg75FxIUR8QDwFCujAz8ATjZdNPb867Hnb3EYEfcCr9qg6sOZ+aVqnw+z8q/oZ7ezba04l36pd2VmRkTPTpGLiJcB/wK8LzP/NyLW6nq1b5n5PHB19VnfF4Hf7HKT2q7nAz8z33S2+oj4Y2A/8MZ8cQ7qE8AVTbtdXpXtGJv16wx2fL820evt38xPIuLSzHwyIi5l5Uqy50TERayE/Wcz8wtVcS36BpCZJyNiFng90BcRu6qr/J5/PdZ6SCcirgU+CFyfmctNVYeAGyLi4oi4ErgK+FY32thmvd6vbwNXVTMjXsLKB9CHutymdjoE3Fht3wj03H9qsXIp/2ng4cz866aqnu5bROxZncUXES8F3szK5xOzwNur3XquX6fJzNo+WPnQ8jHggerxyaa6D7MyRncMeGu323qe/foDVsYTnwV+AnytDv2q2v/7rMyo+gErw1ddb9MW+/GPwJPA/1V/q5uAflZmsDwC3Au8otvt3EK/fpeVDzOPNr2vfr/X+wYMAf9Z9etB4C+r8l9l5aLpOPDPwMXdbmsrD79pK0mFqPWQjiTpRQa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mF+H9q5qfTLaUnxQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": 18,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([12, 14,  3]), array([False,  True,  True]))"
      ]
     },
     "execution_count": 18,
     "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": 19,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True, False])"
      ]
     },
     "execution_count": 19,
     "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": 20,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.broadcast_to(0, a.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 22,
     "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": 23,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11,  7, -3])"
      ]
     },
     "execution_count": 23,
     "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": 24,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [4, 5, 6],\n",
       "       [7, 8, 9]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = np.array([[1, 2, 3], [4,5,6], [7,8,9]]); m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2,  4,  6],\n",
       "       [ 8, 10, 12],\n",
       "       [14, 16, 18]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 2, 2],\n",
       "       [2, 2, 2],\n",
       "       [2, 2, 2]])"
      ]
     },
     "execution_count": 27,
     "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": 28,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10, 20, 30])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = np.array([10,20,30]); c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 22, 33],\n",
       "       [14, 25, 36],\n",
       "       [17, 28, 39]])"
      ]
     },
     "execution_count": 29,
     "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": 30,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 20, 30],\n",
       "       [10, 20, 30],\n",
       "       [10, 20, 30]])"
      ]
     },
     "execution_count": 30,
     "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": 31,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10],\n",
       "       [20],\n",
       "       [30]])"
      ]
     },
     "execution_count": 31,
     "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": 32,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[11, 12, 13],\n",
       "       [24, 25, 26],\n",
       "       [37, 38, 39]])"
      ]
     },
     "execution_count": 32,
     "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": 33,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[10, 10, 10],\n",
       "       [20, 20, 20],\n",
       "       [30, 30, 30]])"
      ]
     },
     "execution_count": 33,
     "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": 34,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 1],\n",
       "       [0, 0],\n",
       "       [9, 4]])"
      ]
     },
     "execution_count": 34,
     "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": 35,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  2],\n",
       "       [ 2,  3],\n",
       "       [-1,  1]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.13598, 0.15958, 0.0564 ])"
      ]
     },
     "execution_count": 36,
     "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": 37,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3,), (3, 2))"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.shape, b.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.13598],\n",
       "       [0.15958],\n",
       "       [0.0564 ]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w[:,None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.13598,  0.27196],\n",
       "       [ 0.31915,  0.47873],\n",
       "       [-0.0564 ,  0.0564 ]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w[:,None]*b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.13288, 2.29314])"
      ]
     },
     "execution_count": 40,
     "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": 41,
   "metadata": {},
   "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": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def meanshift_iter_torch(X, bandwidth=2.5):\n",
    "    out = torch.stack([meanshift_inner(x, X, bandwidth) for x in X], 0)\n",
    "    return to_gpu(out.cuda())\n",
    "\n",
    "def meanshift_torch(X_torch, it=0, max_it=5, bandwidth=2.5, eps=0.000001):\n",
    "    new_X = meanshift_iter_torch(X_torch, bandwidth=bandwidth)\n",
    "    if it >= max_it or abs(X_torch-new_X).sum()/abs(X_torch.sum()) < eps:\n",
    "        return new_X\n",
    "    else:\n",
    "        return meanshift_torch(new_X, it+1, max_it, bandwidth, eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try it out..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.57 s, sys: 96 ms, total: 1.66 s\n",
      "Wall time: 1.66 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE3RJREFUeJzt3W9sZFd5x/Hvk2wIrGhlcLZpJonqtI1BS+WCsCJQ+2LNnxKaVdZUbRVeVCGN7BZRaQmV2lCk2n5RqagSAaklYBfSvEBNKYVluymFJPIqqlQITkm3CWGdJSxK4kC8blclWik0ydMXvnbGu971rmfG47nn+5FGvnPO9b3nyDM/XZ85c25kJpKk+ruo2w2QJG0NA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUiB2tHiAiXg08CFxaHe9LmTkREdcA9wD9wMPA72XmT891rMsuuywHBgZabZIkFeXhhx8+kZm7Ntqv5cAHXgDekZnPR8QlwL9FxNeAjwB3ZOY9EfEZ4FbgznMdaGBggLm5uTY0SZLKERE/PJ/9Wh7SyWXPV08vqR4JvAP4UlV+NzDa6rkkSZvXljH8iLg4Ih4BngPuA74PnMzMF6tdngaubMe5pDpZWlq6oHKpFW0J/Mx8KTPfDFwFXAe88Xx/NyLGI2IuIuYWFxfb0RypJ0xOTjI0NMT8/Pya8vn5eYaGhpicnOxOw1RbbZ2lk5kngVng7UBfRKx8RnAV8MxZfmc6M4czc3jXrg0/c5BqYXJykqmpKRYWFhgZGVkN/fn5eUZGRlhYWGBqasrQV1u1HPgRsSsi+qrt1wDvBh5nOfh/u9rtZuCrrZ5LqoOVsF+xEvr33nvvativMPTVTu24wr8CmI2II8C3gfsy8xDwp8BHIuIYy1MzP9eGc0k9bWlpiZmZmTVlo4xyauEUe/fu5dTCKUZPm98wMzPjmL7aouVpmZl5BHjLOuVPsjyeL6nS39/P7Ozs6pX8KKPsZz/72McUU0wwwQADABzgAI1Gg9nZWfr7+7vbcNWC37SVttjg4CCzs7M0Gg0Oc5jjHGeAAe7iLgYY4DjHOczh1bAfHBzsdpNVEwa+1AWDg4NMT09zkpNMMbWmboopTnKS6elpw15tZeBLXTA/P8/4+Dh99DHBxJq6CSboo4/x8fEzpmxKrTDwpS3WPPVyD3tWh3Fu4ZbV4Z097DljyqbUqsjMbrdh1fDwcLqWjupsaWmJoaGhNVMvRxnlMIc5yUn66GMPezjAgdX6RqPBkSNH/OBWZxURD2fm8Eb7eYUvbaH+/n7GxsbWlB3gADsbOzl06BA7GzvXhD3A2NiYYa+2MPClLTY5OcnExCvj9iuzcW644YbV2TsrJiYm/OKV2qYdyyNLWk/EK9unDZ2uhPjMzMyaqZcrUzZHRkYYGxsz7NVWjuFLnXKOwF+xtLS07nDN2cql9TiGL/WAs4W6Ya9OcEhH6pRt9N+zBF7hS1IxDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IK0XLgR8TVETEbEd+NiMciYn9V/vqIuC8inqh+vq715kqSNqsdV/gvAn+cmbuBtwEfiojdwO3AA5l5LfBA9VyS1CUtB35mPpuZ/1Ft/wR4HLgS2AfcXe12NzDa6rkkSZvX1jH8iBgA3gJ8C7g8M5+tqn4EXN7Oc0mSLkzbAj8iXgv8E/DhzPzf5rrMTGDdOzpHxHhEzEXE3OLiYruaI0k6TVsCPyIuYTnsv5CZX66KfxwRV1T1VwDPrfe7mTmdmcOZObxr1652NEeStI52zNIJ4HPA45n5iaaqg8DN1fbNwFdbPZckafN2tOEYvwb8HvBfEfFIVfZnwF8CX4yIW4EfAr/bhnNJkjap5cDPzH8D4izV72z1+JKk9vCbtpJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIdoS+BHx+Yh4LiIebSp7fUTcFxFPVD9f145zSZI2p11X+H8HXH9a2e3AA5l5LfBA9VyS1CVtCfzMfBD479OK9wF3V9t3A6PtOJckaXM6OYZ/eWY+W23/CLi8g+eSJG1gSz60zcwEcr26iBiPiLmImFtcXNyK5khSkToZ+D+OiCsAqp/PrbdTZk5n5nBmDu/atauDzZGksnUy8A8CN1fbNwNf7eC5JEkbaNe0zL8H/h14Q0Q8HRG3An8JvDsingDeVT2XJHXJjnYcJDPff5aqd7bj+JKk1vlNW0kqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEB0P/Ii4PiKORsSxiLi90+eTJK2vo4EfERcDfwO8F9gNvD8idnfynJKk9XX6Cv864FhmPpmZPwXuAfZ1+JySpHV0OvCvBJ5qev50VbYqIsYjYi4i5hYXFzvcHEkqV9c/tM3M6cwczszhXbt2dbs5klRbnQ78Z4Crm55fVZVJkrZYpwP/28C1EXFNRLwKuAk42OFzSpLWsaOTB8/MFyPij4CvAxcDn8/Mxzp5TknS+joa+ACZ+S/Av3T6PJKkc+v6h7aSpK1h4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIVoK/Ij4nYh4LCJejojh0+o+GhHHIuJoRLyntWZKklq1o8XffxT4LeCzzYURsRu4CXgT0ADuj4jBzHypxfNJkjappSv8zHw8M4+uU7UPuCczX8jMHwDHgOtaOZck1c3S0tIFlbeqU2P4VwJPNT1/uiqTJAGTk5MMDQ0xPz+/pnx+fp6hoSEmJyfbfs4NAz8i7o+IR9d57GtHAyJiPCLmImJucXGxHYeUpG1tcnKSqakpFhYWGBkZWQ39+fl5RkZGWFhYYGpqqu2hv+EYfma+axPHfQa4uun5VVXZesefBqYBhoeHcxPnkqSesRL2K1ZCf3p6mvHxcRYWFlbrVvZrV/B3akjnIHBTRFwaEdcA1wIPdehcktQTlpaWmJmZWVM2yiinFk6xd+9eTi2cYpTRNfUzMzNtG9NvdVrm+yLiaeDtwL0R8XWAzHwM+CLwXeBfgQ85Q0dS6fr7+5mdnaXRaADLYb+f/dzBHQwwwB3cwX72r4Z+o9FgdnaW/v7+tpw/MrfPKMrw8HDOzc11uxmS1FErY/WnFk6thv2K4xznNm5jZ2Mns7OzDA4Obni8iHg4M4c32s9v2krSFhscHGR6epqTnGSKqTV1U0xxkpNMT0+fV9hfCANfkrbY/Pw84+Pj9NHHBBNr6iaYoI8+xsfHz5iy2SoDX5K2UPPUyz3sYYABjnOcW7iF4xxngAH2sOeMKZvt4Bi+JG2RpaUlhoaG1ky9HGWUwxzmJCfpo4897OEAB1brG40GR44cOecHt47hS9I209/fz9jY2JqyAxxgZ2Mnhw4dYmdj55qwBxgbG2vbLB0DX5K20OTkJBMTr4zbr0y9vOGGG9ZM2QSYmJho67dtW10tU5J0gVZCfGZmZs3Uy8HBQWZnZxkZGWFsbKztSys4hi9JbXTjj77DyzvholNw8Offcs59l5aW1h2uOVv52TiGL0ld8PJOuOii5Z8bOVuot2vM/nQGviS10UWn4OWXl39uN47hS1IbrQ7j/Gx327Eer/AlqRAGviQVwsCXpEIUE/hbfbNgSdpuigj8btwsWJK2m9oHfrduFixJ202tA/9sNwu+9957V8N+haEvqe5qG/jdvlmwJG03tQ38bt8sWJK2m9oGPryy8lyj0eAwh1fvJnMXd63eZeYwh1fDvt33j5Sk7aTWgQ/du1mwJG03tQ/8bt0sWJK2m1oHfjdvFixJ201tb4DSqZsFS9J2U/wNULp9s2BJ2m5qG/jQ3ZsFS+p9dVuDq6XAj4i/iojvRcSRiPhKRPQ11X00Io5FxNGIeE/rTT23Q8/fufpothL6p0+9bJ6yadhLOl0t1+DKzE0/gN8AdlTbHwc+Xm3vBv4TuBS4Bvg+cPFGx3vrW9+am/XPP/n06mM9J06cuKBySeWamJhIIIFsNBp59OjRzMw8evRoNhqN1bqJiYnuNrQCzOV5ZHZLV/iZ+Y3MfLF6+k3gqmp7H3BPZr6QmT8AjgHXtXKuVm31zYIl9aY6r8HVzjH83we+Vm1fCTzVVPd0VdYxe1/7wdWHJG1G3dfg2jDwI+L+iHh0nce+pn0+BrwIfOFCGxAR4xExFxFzi4uLF/rrktQ2dV+Da8dGO2Tmu85VHxEfAPYC76zGkgCeAa5u2u2qqmy9408D07A8D3/jJktS56xM6BgZGeHwwmH2sW91DS6gp9fganWWzvXAnwA3ZuappqqDwE0RcWlEXANcCzzUyrkkaavUdQ2uVsfw/xr4GeC+iHgkIj4DkJmPAV8Evgv8K/ChzHypxXNJ0pao6xpcrc7S+eXMvDoz31w9/rCp7i8y85cy8w2Z+bVzHUeStos6r8FV27V0JOlC9eoaXMWvpSNJ53LjzMM8+bfj3Djz8GpZ3dfgMvAlFemT8Vl+8f+e45Px2TXldV6Dy8CXVKQP5x/w5CU/x4fzD86oq+saXI7hS2qbpaWldYc3zla+3fVKfxzDV1fUbTlZnb86ri5ZtzW4DHy1TR3f8Do/KwuOnT5VsXmKY68tNFZL57Ok5lY9WlkeWd3Va8vJqn2a//bNr4FDhw6t+dv7GugcznN55K6HfPPDwO9NvuHLdeLEiTP+xqOMZh99CWQffTnK6BmvDe9D0V7nG/gO6agldV9OVudW99Ul68bAV0t8w6t5quJhDq8uP3AXd60uS9Crq0vWjYGvlvmGV11Xl6wbA19t4Ru+bHVdXbJuDHy1hW/4ctV5dcm6MfDVMt/w5VpaWlpzY+8DHOBTfIrbuI3jHOc2buNTfGp1wbGV14Af2neHga+W+IYvW91Xl6wbA18t8Q1flufvHOXlO0d5/s5XptrWeXXJ2jmfyfpb9fCLV73Lb9qW4aVP78v89L7ln6eZmJhY87dfsfIa8G/fOZznF69cLVNtMzk5yczMzBlTL1fG+MfGxry663HP3znKTuAU8NoPHjijvldWl6yb810t08DXeZuKV7YnzvKy8Q0vbT2XR1ZX1G05WalODHxJKsSObjdAveNswziSeoNX+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFWJbLa0QEYvADztw6MuAEx04brfZr95iv3pLL/XrFzJz10Y7bavA75SImDufdSZ6jf3qLfart9SxXw7pSFIhDHxJKkQpgT/d7QZ0iP3qLfart9SuX0WM4UuSyrnCl6Ti1TrwI+KvIuJ7EXEkIr4SEX1NdR+NiGMRcTQi3tPNdl6oiPidiHgsIl6OiOHT6nq2XwARcX3V9mMRcXu327NZEfH5iHguIh5tKnt9RNwXEU9UP1/XzTZuRkRcHRGzEfHd6jW4vyrv6b5FxKsj4qGI+M+qX1NV+TUR8a3q9fgPEfGqbre1FbUOfOA+4FcycwiYBz4KEBG7gZuANwHXA5+OiIu71soL9yjwW8CDzYW93q+qrX8DvBfYDby/6lMv+juW/wbNbgceyMxrgQeq573mReCPM3M38DbgQ9XfqNf79gLwjsz8VeDNwPUR8Tbg48AdmfnLwP8At3axjS2rdeBn5jcy88Xq6TeBq6rtfcA9mflCZv4AOAZc1402bkZmPp6ZR9ep6ul+sdzWY5n5ZGb+FLiH5T71nMx8EPjv04r3AXdX23cDo1vaqDbIzGcz8z+q7Z8AjwNX0uN9y2XPV08vqR4JvAP4UlXec/06Xa0D/zS/D3yt2r4SeKqp7umqrNf1er96vf0buTwzn622fwRc3s3GtCoiBoC3AN+iBn2LiIsj4hHgOZZHB74PnGy6aOz512PP3+IwIu4Hfn6dqo9l5lerfT7G8r+iX9jKtrXifPql3pWZGRE9O0UuIl4L/BPw4cz834hYrevVvmXmS8Cbq8/6vgK8sctNarueD/zMfNe56iPiA8Be4J35yhzUZ4Crm3a7qirbNjbq11ls+35toNfbv5EfR8QVmflsRFzB8pVkz4mIS1gO+y9k5per4lr0DSAzT0bELPB2oC8idlRX+T3/eqz1kE5EXA/8CXBjZp5qqjoI3BQRl0bENcC1wEPdaGOb9Xq/vg1cW82MeBXLH0Af7HKb2ukgcHO1fTPQc/+pxfKl/OeAxzPzE01VPd23iNi1MosvIl4DvJvlzydmgd+uduu5fp0hM2v7YPlDy6eAR6rHZ5rqPsbyGN1R4L3dbusF9ut9LI8nvgD8GPh6HfpVtf83WZ5R9X2Wh6+63qZN9uPvgWeB/6v+VrcC/SzPYHkCuB94fbfbuYl+/TrLH2YeaXpf/Wav9w0YAr5T9etR4M+r8l9k+aLpGPCPwKXdbmsrD79pK0mFqPWQjiTpFQa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mF+H/Y7qdg1/eFMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X_torch = to_gpu(torch.from_numpy(X))\n",
    "%time X = meanshift_torch(X_torch).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": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_b(a,b): return sqrt(((a[None,:] - b[:,None]) ** 2).sum(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " 0.6297  0.8461  0.9324\n",
       " 0.1305  0.2712  0.1951\n",
       "[torch.FloatTensor of size 2x3]"
      ]
     },
     "execution_count": 45,
     "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": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "def meanshift_gpu(X, it=0, max_it=5, bandwidth=2.5, eps=0.000001):\n",
    "    weights = gaussian(distance_b(X, X), bandwidth)\n",
    "    num = (weights[:,:,None] * X).sum(1)\n",
    "    X_new = num / weights.sum(1)[:,None]\n",
    "    \n",
    "    if it >= max_it or abs(X_new - X).sum()/abs(x.sum()) < eps:\n",
    "        return X_new\n",
    "    else:\n",
    "        return meanshift_gpu(X_new, it+1, max_it, bandwidth, eps)"
   ]
  },
  {
   "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": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 60 ms, sys: 60 ms, total: 120 ms\n",
      "Wall time: 124 ms\n"
     ]
    }
   ],
   "source": [
    "X_torch = to_gpu(torch.from_numpy(data))\n",
    "%time X = meanshift_gpu(X_torch).cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's more like it! We've gone from 1660ms to 124ms, which is a speedup of 13.5! Oh, and it even gives the right answer!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE5VJREFUeJzt3W9sZFd5x/Hvk2xIWdHKYLYhk0R12satlsoNwopA7Ys1f0poVlkXQRVeVCmNbIqoFAISDUWq7ReVQEhNkVoCdiHNC9SUUli2G/4lkVeoUgs4NGwTwjpLWJTEgTim2xJZCk3y9IWvnfGud73rmfF47vl+pJHvnHN97znyzE/XZ86cG5mJJKn+Luh2AyRJ28PAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBViV6sHiIhfAL4BXFwd7/OZORERVwJ3Af3A/cAfZebPz3asV77ylTkwMNBqkySpKPfff//Tmblns/1aDnzgWeANmflMRFwE/FtEfAV4P3BbZt4VEZ8EbgJuP9uBBgYGmJuba0OTJKkcEfGjc9mv5SGdXPFM9fSi6pHAG4DPV+V3AqOtnkuStHVtGcOPiAsj4gHgKeAe4AfAycx8rtrlceCydpxLqpOlpaXzKpda0ZbAz8znM/Nq4HLgGuA3z/V3I2I8IuYiYm5xcbEdzZF6wuTkJENDQ8zPz68rn5+fZ2hoiMnJye40TLXV1lk6mXkSmAVeD/RFxOpnBJcDT5zhd6Yzczgzh/fs2fQzB6kWJicnmZqaYmFhgZGRkbXQn5+fZ2RkhIWFBaampgx9tVXLgR8ReyKir9p+KfBm4GFWgv/t1W43Al9q9VxSHayG/arV0L/77rvXwn6Voa92ascV/qXAbEQcBb4N3JOZh4E/B94fEcdZmZr56TacS+ppS0tLzMzMrCsbZZTlhWX279/P8sIyo6fMb5iZmXFMX23R8rTMzDwKvGaD8kdZGc+XVOnv72d2dnbtSn6UUW7mZg5wgCmmmGCCAQYAOMhBGo0Gs7Oz9Pf3d7fhqgW/aStts8HBQWZnZ2k0GhzhCCc4wQAD3MEdDDDACU5whCNrYT84ONjtJqsmDHypCwYHB5menuYkJ5lial3dFFOc5CTT09OGvdrKwJe6YH5+nvHxcfroY4KJdXUTTNBHH+Pj46dN2ZRaYeBL26x56uU+9q0N47yLd60N7+xj32lTNqVWRWZ2uw1rhoeH07V0VGdLS0sMDQ2tm3o5yihHOMJJTtJHH/vYx0EOrtU3Gg2OHj3qB7c6o4i4PzOHN9vPK3xpG/X39zM2Nrau7CAH2d3YzeHDh9nd2L0u7AHGxsYMe7WFgS9ts8nJSSYmXhy3X52Nc911163N3lk1MTHhF6/UNu1YHlnSRiJe3D5l6HQ1xGdmZtZNvVydsjkyMsLY2Jhhr7ZyDF/qlLME/qqlpaUNh2vOVC5txDF8qQecKdQNe3WCQzpSp+yg/54l8Apfkoph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klSIlgM/Iq6IiNmI+F5EPBQRN1flr4iIeyLikerny1tvriRpq9pxhf8c8IHM3Au8DnhvROwFbgXuy8yrgPuq55KkLmk58DPzycz8TrX9M+Bh4DLgAHBntdudwGir55IkbV1bx/AjYgB4DfBN4JLMfLKq+jFwSTvPJUk6P20L/Ih4GfAvwPsy83+b6zIzgQ3v6BwR4xExFxFzi4uL7WqOJOkUbQn8iLiIlbD/bGZ+oSr+SURcWtVfCjy10e9m5nRmDmfm8J49e9rRHEnSBtoxSyeATwMPZ+ZfN1UdAm6stm8EvtTquSRJW7erDcf4HeCPgP+KiAeqsr8APgJ8LiJuAn4E/GEbziVJ2qKWAz8z/w2IM1S/sdXjS5Law2/aSlIhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFaEvgR8RnIuKpiHiwqewVEXFPRDxS/Xx5O84lSdqadl3h/wNw7SlltwL3ZeZVwH3Vc0lSl7Ql8DPzG8BPTyk+ANxZbd8JjLbjXJKkrenkGP4lmflktf1j4JIOnkuStIlt+dA2MxPIjeoiYjwi5iJibnFxcTuaI0lF6mTg/yQiLgWofj610U6ZOZ2Zw5k5vGfPng42R5LK1snAPwTcWG3fCHypg+eSJG2iXdMy/xH4d+A3IuLxiLgJ+Ajw5oh4BHhT9VyS1CW72nGQzHznGare2I7jS5Ja5zdtJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCdDzwI+LaiDgWEccj4tZOn0+StLGOBn5EXAj8HfBWYC/wzojY28lzSpI21ukr/GuA45n5aGb+HLgLONDhc0qSNtDpwL8MeKzp+eNV2ZqIGI+IuYiYW1xc7HBzJKlcXf/QNjOnM3M4M4f37NnT7eZIUm11OvCfAK5oen55VSZJ2madDvxvA1dFxJUR8RLgBuBQh88pSdrArk4ePDOfi4g/A74GXAh8JjMf6uQ5JUkb62jgA2Tml4Evd/o8kqSz6/qHtpKk7WHgS1IhDHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhWgr8iHhHRDwUES9ExPApdR+KiOMRcSwi3tJaMyVJrdrV4u8/CLwN+FRzYUTsBW4AXg00gHsjYjAzn2/xfJKkLWrpCj8zH87MYxtUHQDuysxnM/OHwHHgmlbOJUl1s7S0dF7lrerUGP5lwGNNzx+vyiRJwOTkJENDQ8zPz68rn5+fZ2hoiMnJybafc9PAj4h7I+LBDR4H2tGAiBiPiLmImFtcXGzHISVpR5ucnGRqaoqFhQVGRkbWQn9+fp6RkREWFhaYmppqe+hvOoafmW/awnGfAK5oen55VbbR8aeBaYDh4eHcwrkkqWeshv2q1dCfnp5mfHychYWFtbrV/doV/J0a0jkE3BARF0fElcBVwLc6dC5J6glLS0vMzMysKxtllOWFZfbv38/ywjKjjK6rn5mZaduYfqvTMv8gIh4HXg/cHRFfA8jMh4DPAd8Dvgq81xk6kkrX39/P7OwsjUYDWAn7m7mZ27iNAQa4jdu4mZvXQr/RaDA7O0t/f39bzh+ZO2cUZXh4OOfm5rrdDEnqqNWx+uWF5bWwX3WCE9zCLexu7GZ2dpbBwcFNjxcR92fm8Gb7+U1bSdpmg4ODTE9Pc5KTTDG1rm6KKU5ykunp6XMK+/Nh4EvSNpufn2d8fJw++phgYl3dBBP00cf4+PhpUzZbZeBL0jZqnnq5j30MMMAJTvAu3sUJTjDAAPvYd9qUzXZwDF+StsnS0hJDQ0Prpl6OMsoRjnCSk/TRxz72cZCDa/WNRoOjR4+e9YNbx/AlaYfp7+9nbGxsXdlBDrK7sZvDhw+zu7F7XdgDjI2NtW2WjoEvSdtocnKSiYkXx+1Xp15ed91166ZsAkxMTLT127atrpYpSTpPqyE+MzOzburl4OAgs7OzjIyMMDY21valFRzDl6Q2uv7H/8kLu+GCZTj0qtecdd+lpaUNh2vOVH4mjuFLUhe8sBsuuGDl52bOFOrtGrM/lYEvSW10wTK88MLKz53GMXxJaqO1YZxf6m47NuIVviQVwsCXpEIY+JJUiGICf7tvFixJO00Rgd+NmwVL0k5T+8Dv1s2CJWmnqXXgn+lmwXffffda2K8y9CXVXW0Dv9s3C5aknaa2gd/tmwVL0k5T28CHF1eeazQaHOHI2t1k7uCOtbvMHOHIWti3+/6RkrST1DrwoXs3C5aknab2gd+tmwVL0k5T68Dv5s2CJWmnqe0NUDp1s2BJ2mmKvwFKt28WLEk7TW0DH7p7s2BJva9ua3C1FPgR8bGI+H5EHI2IL0ZEX1PdhyLieEQci4i3tN7Uszv8zO1rj2aroX/q1MvmKZuGvaRT1XINrszc8gP4PWBXtf1R4KPV9l7gu8DFwJXAD4ALNzvea1/72tyqf/3ZJ9YeG3n66afPq1xSuSYmJhJIIBuNRh47diwzM48dO5aNRmOtbmJiorsNrQBzeQ6Z3dIVfmZ+PTOfq57+B3B5tX0AuCszn83MHwLHgWtaOVertvtmwZJ6U53X4GrnGP6fAF+pti8DHmuqe7wq65j9L3vP2kOStqLua3BtGvgRcW9EPLjB40DTPh8GngM+e74NiIjxiJiLiLnFxcXz/XVJapu6r8G1a7MdMvNNZ6uPiD8G9gNvrMaSAJ4Armja7fKqbKPjTwPTsDIPf/MmS1LnrE7oGBkZ4cjCEQ5wYG0NLqCn1+BqdZbOtcAHgeszc7mp6hBwQ0RcHBFXAlcB32rlXJK0Xeq6BlerY/h/C/wicE9EPBARnwTIzIeAzwHfA74KvDczn2/xXJK0Leq6Blers3R+PTOvyMyrq8efNtX9VWb+Wmb+RmZ+5WzHkaSdos5rcNV2LR1JOl+9ugZX8WvpSNLZXD9zP4/+/TjXz9y/Vlb3NbgMfElF+pv4FL/6f0/xN/GpdeV1XoPLwJdUpPflu3n0ol/mffnu0+rqugaXY/iS2mZpaWnD4Y0zle90vdIfx/DVFXVbTlbnro6rS9ZtDS4DX21Txze8zs3qgmOnTlVsnuLYawuN1dK5LKm5XY9WlkdWd/XacrJqn+a/ffNr4PDhw+v+9r4GOodzXB656yHf/DDwe5Nv+HI9/fTTp/2NRxnNPvoSyD76cpTR014b3oeivc418B3SUUvqvpyszq7uq0vWjYGvlviGV/NUxSMcWVt+4A7uWFuWoFdXl6wbA18t8w2vuq4uWTcGvtrCN3zZ6rq6ZN0Y+GoL3/DlqvPqknVj4KtlvuHLtbS0tO7G3gc5yMf5OLdwCyc4wS3cwsf5+NqCY6uvAT+07w4DXy3xDV+2/t0XMbb/d9eV1Wl1ybox8NWSui8nq/WeuX2UF24f5Znbq6m237+PyaufZWL8D9f2qdPqkrVzLpP1t+vhF696l9+0LcPznziQ+YkDKz8zM5f/J/M7X8hc/p+cmJhY97dftfoa8G/fOZzjF69cLVNtMzk5yczMzGlTL1fH+MfGxry663HP3D7KbmAZeNl7Dp5W3yurS9bNua6WaeDrnE3Fi9sTZ3jZ+IaXtp/LI6sr6racrFQnBr4kFWJXtxug3nGmYRxJvcErfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCrGjllaIiEXgRx049CuBpztw3G6zX73FfvWWXurXr2Tmns122lGB3ykRMXcu60z0GvvVW+xXb6ljvxzSkaRCGPiSVIhSAn+62w3oEPvVW+xXb6ldv4oYw5cklXOFL0nFq3XgR8THIuL7EXE0Ir4YEX1NdR+KiOMRcSwi3tLNdp6viHhHRDwUES9ExPApdT3bL4CIuLZq+/GIuLXb7dmqiPhMRDwVEQ82lb0iIu6JiEeqny/vZhu3IiKuiIjZiPhe9Rq8uSrv6b5FxC9ExLci4rtVv6aq8isj4pvV6/GfIuIl3W5rK2od+MA9wG9l5hAwD3wIICL2AjcArwauBT4RERd2rZXn70HgbcA3mgt7vV9VW/8OeCuwF3hn1ade9A+s/A2a3Qrcl5lXAfdVz3vNc8AHMnMv8DrgvdXfqNf79izwhsz8beBq4NqIeB3wUeC2zPx14L+Bm7rYxpbVOvAz8+uZ+Vz19D+Ay6vtA8BdmflsZv4QOA5c0402bkVmPpyZxzao6ul+sdLW45n5aGb+HLiLlT71nMz8BvDTU4oPAHdW23cCo9vaqDbIzCcz8zvV9s+Ah4HL6PG+5YpnqqcXVY8E3gB8virvuX6dqtaBf4o/Ab5SbV8GPNZU93hV1ut6vV+93v7NXJKZT1bbPwYu6WZjWhURA8BrgG9Sg75FxIUR8QDwFCujAz8ATjZdNPb867Hnb3EYEfcCr9qg6sOZ+aVqnw+z8q/oZ7ezba04l36pd2VmRkTPTpGLiJcB/wK8LzP/NyLW6nq1b5n5PHB19VnfF4Hf7HKT2q7nAz8z33S2+oj4Y2A/8MZ8cQ7qE8AVTbtdXpXtGJv16wx2fL820evt38xPIuLSzHwyIi5l5Uqy50TERayE/Wcz8wtVcS36BpCZJyNiFng90BcRu6qr/J5/PdZ6SCcirgU+CFyfmctNVYeAGyLi4oi4ErgK+FY32thmvd6vbwNXVTMjXsLKB9CHutymdjoE3Fht3wj03H9qsXIp/2ng4cz866aqnu5bROxZncUXES8F3szK5xOzwNur3XquX6fJzNo+WPnQ8jHggerxyaa6D7MyRncMeGu323qe/foDVsYTnwV+AnytDv2q2v/7rMyo+gErw1ddb9MW+/GPwJPA/1V/q5uAflZmsDwC3Au8otvt3EK/fpeVDzOPNr2vfr/X+wYMAf9Z9etB4C+r8l9l5aLpOPDPwMXdbmsrD79pK0mFqPWQjiTpRQa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mF+H9q5qfTLaUnxQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": {},
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}