mirror of
https://github.com/zhigang1992/fastai.git
synced 2026-04-24 05:14:58 +08:00
1377 lines
108 KiB
Plaintext
1377 lines
108 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Clustering with pytorch"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Clustering techniques are unsupervised learning algorithms that try to group unlabelled data into \"clusters\", using the (typically spatial) structure of the data itself.\n",
|
|
"\n",
|
|
"The easiest way to demonstrate how clustering works is to simply generate some data and show them in action. We'll start off by importing the libraries we'll be using today."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"%matplotlib inline\n",
|
|
"import math, numpy as np, matplotlib.pyplot as plt, operator, torch"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Create data"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"n_clusters=6\n",
|
|
"n_samples =250"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"To generate our data, we're going to pick 6 random points, which we'll call centroids, and for each point we're going to generate 250 random points about it."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"centroids = np.random.uniform(-35, 35, (n_clusters, 2))\n",
|
|
"slices = [np.random.multivariate_normal(centroids[i], np.diag([5., 5.]), n_samples)\n",
|
|
" for i in range(n_clusters)]\n",
|
|
"data = np.concatenate(slices).astype(np.float32)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 32,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Save data so we can show consistent results\n",
|
|
"np.save('tmp/meanshift_data.npy', data)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"data = np.load('tmp/meanshift_data.npy')"
|
|
]
|
|
},
|
|
{
|
|
"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": 9,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def plot_data(centroids, data, n_samples):\n",
|
|
" colour = plt.cm.rainbow(np.linspace(0,1,len(centroids)))\n",
|
|
" for i, centroid in enumerate(centroids):\n",
|
|
" samples = data[i*n_samples:(i+1)*n_samples]\n",
|
|
" plt.scatter(samples[:,0], samples[:,1], c=colour[i], s=1)\n",
|
|
" plt.plot(centroid[0], centroid[1], markersize=10, marker=\"x\", color='k', mew=5)\n",
|
|
" plt.plot(centroid[0], centroid[1], markersize=5, marker=\"x\", color='m', mew=2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvXt8VNW5//9ek5nJhTshhCQQAoaLgqWUgKhVQdFS4Ej0\nnNZLT39KrZxTL/XS6mm1bYitvcip9qK2B4+i/Z1WaU+FcoSqBQneQA0iihogSgRCCCQQAuQyt/X9\nY82e2TOZ3EhCJpnn/XrlNdmX2XttNJ/17M961rOU1hpBEASh/+Po7QYIgiAIZwYRfEEQhARBBF8Q\nBCFBEMEXBEFIEETwBUEQEgQRfEEQhARBBF8QBCFBEMEXBEFIEETwBUEQEgRnbzfAzogRI3ReXl5v\nN0MQBKFPsW3bthqtdUZ758WV4Ofl5VFaWtrbzRAEQehTKKU+68h5YukIgiAkCCL4giAICYIIviAI\nQoIggi8IgpAgiOALgiAkCCL4giAICYIIviAIQoIggi8IgpAgiOALgtA6jfWwfbX5FPo8IviCILRO\n2UbY+oz5FPo8cVVaQRCEOGPyZZGfQp9GBF8QhNZJHQzTr+rtVgjdhFg6giAICYIIviAIQoIggi8I\ngpAgiOALgiAkCCL4giAICYIIviD0FWQSlNBFRPAFoa8QL5OgznTHIx1dtyF5+ILQVzhTk6Aa602n\nMvkyk4cfjdXxwJnJ0T/T9+vHiOALQl+huyZBtSfoH6yDbavA2wSzrmv5nTM9+1Zm+3YbYukIQn8n\n2hJpzxpSUZ/R37E6nlidRU9wpu/Xj5EIXxD6O9GWSGsRsxXF518MzpTI452Jstt7gxB6DRF8Qeir\ndFRY7WJtfSdvVuR3G+vhlV/Dvm3hc63jYGweFfvyLdpkv4547nFFlwVfKZUCvAokB6/3v1rrIqXU\ncGAVkAdUAF/VWh/r6v0EIeGxRNvbZLx26JiwNp2AN58yYrz/PajcEfbpyzaa/bkzwmJvvRVA+D6H\ny+HSO1rvYKzr5EwDX5Npq0T5cUN3RPjNwKVa65NKKRfwulLq78DVwEat9c+VUt8Dvgf8RzfcTxAS\nG0uMC66B2TeEo3B79F7xdjh637kOSlfBZ6VQ9aER46HZRvAP7jTfs78FpA6O3G46YTqIgM+IednG\n1jsY63u+JnNPZ4pE+XFElwVfa62Bk8FNV/BHA4uBOcH9zwAliOALQtfJm2WEOv9iGJYT3m91BAd3\nRloq3mbzu99nPkeMg8N7zO9VHxoL5tI7jFjvXGf+es9dGBbqso2mcyi4BsYWRFpD0XaSNcDaWN9y\nHEDodbrFw1dKJQHbgHzgMa31W0qpTK11VfCUQ0BmK99dCiwFyM3N7Y7mCEL/puJtI+jZU2GYLXq2\nxDVvljlmbbuSzWfW2TB+tom+qz40+wZlmmu98msYmW+icoAj5XDBN8y98maZ72jbPdrz6VtLIZUB\n3V6lWwRfa+0HPq+UGgqsVkpNjTqulVK6le+uAFYAFBQUxDxHEAQbbWXM+Jpgz6smQrcEdepCcKZQ\nO3gC6Uc/MBF/5iSo3gVpQ6llMOn7tsHR/ZA5GdBGzP0+E9mDida3PgOuFLO9bxsMzYGsc0zKZ0cF\nXCZR9SrdmqWjta5TSm0C5gPVSqksrXWVUioLONyd9xKEhCBWRNxa9Fy2sWWEXv4qeJtZ9uQanljz\nCptun8XEzIFG8IHd729j7uPvcvNFZ7FsHnDysIn6wVg/I8bB3q1me1phuJOxbKNtf+5cRo5MoupV\nuiNLJwPwBsU+Fbgc+AWwFrgB+Hnw829dvZcgJBwdiYgb64337m2GKQtg/7tGhI8fhONVLHuhjOL1\nuwCY+5stbPr2+UzM8LO7+iRzf/sWB4+epPj5UtCaZZePhhPVkDIYcj4HbzwBx4POrDMl3Olceodp\nW9Y5ZjtvVseeR5ZM7FW6Y6ZtFrBJKfU+8A7wD631Cxihv1wptQeYF9wWBKE97DNjJ1/WMhPHPmvW\nynsvXQU71sDxSqg/ZAQ7SuwBDh5rYO5jpax7bRtzf/UGB4+eDB0rXr2NZSW1kDwQmurhzSfDYg8m\n2j9WCet+bDJ3pl8FVR+ZzqXi7TPxLyN0kS4Lvtb6fa31dK3157TWU7XWDwT312qtL9NaT9Baz9Na\nH+16cwWhn2MJeKwyBscqYc33I8siWHnvg7PMtt8Ho6ZAUz21abk8sbUy4vKFFNJQo1j0mxIajqdQ\nSGHE8SdeKqW29qixdXKmwcCgvTN4FEy/OpzH/+ZTZn90h2Q9g1S3jEuklo4gxBPRE6DsvPkU1FWa\nwVK7Fz77BrjsTrO/6kNwGqc2fcggNt1+HtkZwwAj9ndwB4/wCHnk8QiPcAd3hEQ/OzODTQ/eQPpA\nNwwcAR+uh7EzjNiP+YK53wXfMB3BoJHwzrNmX3Sdm3gp4yy0QARfEOIJS8Av+IYRTHuUfME3TEcw\n/76WA7j73jWdweAsmHm9GWA9WcPEzIFsKvoq2SOGUkIJFVSQRx4rWUkeeVRQQQklZA8fxKb/+ysT\nc9LNdQM+k3efMshYRB+uN+0ZlgNjPg8f/t3YSNGi3lhvZu8WXCMDs3GICL4gxCPlr7aMkoflwMIf\nRk62svAFJ1fVVxlf/dj+0ODrRMchVnxlEnXUUUxxxNeKKaaOOlZccw4Td/8pnJ9fvcuUUZhwsek8\nRk8LD8xOvgxmXBNb1Ms2mjIM9gFeIW6Q4mmCEE9YdsiMa1p649E01rcsajYwExqPw4yvmu0ZX2X3\n6+tY+twPGcpQiiiKuEQRRdzFXSz90w42jRzIxBkXQvo4qN1rrKWMfNN5HNhhBmZTgnV27Hn+dqxZ\nwB3N2hHOKBLhC0I8YVk6VmkDu6jGqmu/bVU4935oDpysNtk6Gx6GC77B7nrF3Dt+ycFjDcxhTsjG\nWcKSkL0zhzkcPN7E3N++ze5R8yB1CHxxqWmHInZRNfubh71d1ixgydqJSyTCF4R4wl6LJnoGa6y6\n9t4mI8oa4+EnucHvgRPV1L74KHO/+wcOHjwIwBrWAFBCCXXUcRd3MYc5of0Hj55g7sLFvH/PbNLn\n0bImTnRRNQt7u2RiVVwjEb4gxCOxIml7CqQ1A/fchTDzOuO1D80xYj8oEzInkT52Ejf/6zURl13D\nGtKG+3jhO/NJy1Ahsbe4+bJzSZ/z9bBgR682FWv1KXu7ZHWquEYEXxDiEWtg1KopD5FiGt0hVLxt\nIvzcGfDPy2HcbNixhmUXD6NowaTQZbOHprLp9vNYOPcCNv3k62QPSQkdK7rycyy7eFh4wLWj+fQi\n8n0GsXQEIR5JHWwKlW19pmVNeSv1cVqh+YxVz94aPJ3xVZYV5cPIP/LE2lfZ9OzvmVj7Bsz4KhOn\nnmJTTQVzHy7h5ss/z7JHVpgMH+ta9nLL1qInUu2yTyOCLwjxSrQfbomtL7jSVe4MM0DqCnYI9k7B\nXkJ51nUsm3Udt3+yk/S3f2feBKo+goM7mThM8f5P/4X0rz9iBHzUpMj7W0XSLJFvbRlEEf8+gVg6\nghCv2K2Sxnr4xy9NxO1tNjnwg0aaWa+xUiCjLaHG+rDYD84y+2d8FYbmkO6vNUJ+rLKlhZORH863\n/2BdePnC1jJ2WkPKLcQFEuF3I7WNmlVlXq6Z7CI9tSMrPgtCB7FWnQKzoIkzBT4MpmNWvB25EAqY\n4mafvG4E3hn06a2yDGNnhpcfLPxZZNRuL3VspX3OvsHsqy6LvEdrGTmxbB+pgx8XJLTgxxLoroj2\nqjIvP9nqAeCW6e5ub6+QwNhTMKcuNPuiV6GyaKyHF3/asu4OmLeBPa9GzpLNyDerXeVfHLlSlv2z\nbKOZfDU0x3Q81rq2rdXljxZ3SdeMCxLa0rEEelWZt819HeWayS5+MNvNNZNd3dlMIVGx2yCpg2HW\ndUbsLQtl5nUmLTO65k7ZxrDYz/12+PzpV5m3gW2rTNkE61yrFMKwnLBAb18d/g6E6+PMv6/9GcCx\nKmhKJk9ckNARviXMdoGOta+jpKcqieyF7iNWpPzBOiPQ3ibTAbQXTcearBU9EAsm8rcmekV/x27t\nDMtpaR9FI4ucxC0JLfixBPpMiLZ4/UKHiGWDqKjPWOfYBTf6eOrg8GpV9olS21eHRT66Ho7YMf2G\nhLZ0upPaRs3j2z3UNuqY2xblxwJctabhtG2j022P0AeJZYNMXWgibcvHb88qiXW8vdmyrdXDaToh\nmTZ9nISO8GPR0eg7+rzoAdvobev8zft9lNdp8ocqrpns6rFoXwaQ+yk9ZZe09VZgn4DVmQXLhbhD\nBD+KWEIZS5Tt510z2UWDV3N3gavFGID1aZ3/rWku3EmKoguSSU9VPL7d0yPC3JWxCCGBiZVSaff5\n7Vk8Qp9DBD+KWEIZqxOwn7eqzMvD27xclpsEhDuIK/KcoY7Cfr49krf2X5Hn5PHtnm6L9GUAuR/S\nmbIG1rl5s4w109HZsLEGge3Rf3sDtkJck3CC356FEksoY3UC9vOumexiy0E/G/f5WbnTw47DATbu\n87P5gI/XDgTYfMDH4/NSI94YVn7godEHqU5Ycm5LC0gQWtCZyUvRNoy9Hk5byABtvybhBL8zwmqP\n1KP3r9zpAW3EGmDaSAfTMhygYeM+P5flJpGRagZMXzsQYFWZN3Q/643AIs2lxIIR2qczYmy3YSCc\nhmnVuG/tTUFSKvs1CSf4lt/e4NPUNuo2Z9hanYMVvW856OdXl6YYwS41gp3mMuc+XOqNmHTV6IO1\n5eac3EGRA7RX5Dlp8OpQhG/d07KHJF1TiElnxNh+rj0NE6TMQQLTZcFXSo0B/gBkYiZ6r9Ba/1op\nNRxYBeQBFcBXtdbHunq/rpKeqkhzqVCUn+YMC21rPr2J8JvZuM8fEuTaRs3OGj9X5DkZlqLYf6iG\nDRVDqG0MkOpS/G6HEfvAyVqGZaQD8Pj2Zn63w8eGCi9PzE+LEPXaRs3NLzWwtUpT2xjglunJIv5C\n9xDdUYhtk7B0Rx6+D/iO1vocYDZwq1LqHOB7wEat9QRgY3A7LrBKIKDhJ1s93PlKE1fkOUMRupXD\nDkb884c5+NWlKdw9w0WDT7P3eIAXK3y8VhlgTbmX3/6imF//awGv79jN73b4QMPdBS7ST+zhyIPn\n8/rKn7KqzMvO2gAAWw/pFjn4q8q8bK0yFtDOmkCoA1q509NqPr3k2gvtEqtKpZQ5SFi6HOFrrauA\nquDvJ5RSHwM5wGJgTvC0Z4AS4D+6er/uwBpwrW3U7DhiBlg9/iZmjjJZNpbYNnh1yF9PT1WgjHWz\neo+PvceNyP7t8Qd4/emfAuD79UIu/PF6Cid8jsCRch5YvojA8SpOrv8Zfxnl4Kyv3EdKEnw5T9k6\nlma2Hw4webjixilOPqkL8OBFKQxLMVF9g1e3miba4NMha0kGeoWYiH0j2OhWD18plQdMB94CMoOd\nAcAhjOUTV6SnKn51aQp3vtLExn1+XqsMRAygNvgixbbRa0T+vFGKumbNvv/9Ka+v/1noet66Kl79\nwQJ+5n2cv/z0Fk7UVIWOlTz1IKWHAgxadB/v1RDKwf/dDh8AW6vgB7Pd/PTi8JJzVqdkbxOEO6S7\nZ0ixNqEdxL4RbHRbaQWl1EDgr8CdWuuIuddaa43x92N9b6lSqlQpVXrkyJHuak6HsUT/7gIXd89w\nhXLnr8hz0uiFi3IcoSyd1OAAbfagJJ6+uBH/1qcjrlVIIYOPN/LU3VeRVNNIIYURxxveeJrAyVou\nGe2gttEMHN84JYnZWQ6+Nc0ZU7ittxH74HKDT3P3DBdLznVHHBOEFoh9I9joFsFXSrkwYv9HrfXz\nwd3VSqms4PEs4HCs72qtV2itC7TWBRkZGd3RnE6Tnqr47sxklpzrpvjNZn6y1UPxm838boeX1yoD\n3FvSyNVrTnG0MSi0U93MPCuDe558GceQLMCI/R3cwSM8Qh55PMIj3MEdIdFPGZ7FhQ+sY9ZZGaQ6\nFSt3eni41MvoQUk8dEkKu4+ZsYG2PPnaRs2drzTxcKmXNJcSoRcEoVN0R5aOAp4EPtZaP2w7tBa4\nAfh58PNvXb1XT7Nyp4eN+/yMG6y4Y4abicN8vLjXy9ZDRoC3VvlCs2kBvv2lKfD0y/zXLV+ipKqE\nxSwmjzxWshKACioooYRRWdnMWvZ33vGPh+oApdUBvjUtPEh8y4ZGXjsQ4NPjTaGxgVie/KoybyjH\nX2wcQRA6S3dE+BcCXwcuVUq9F/xZgBH6y5VSe4B5we0eo7MZK7HOtzz6vfWat6v8pKcq9tZD7kBz\nPHeQmVR15ytNoRz+4sKpPPnECuqoo5jiiHsUU0wddUxe+lve8Y9naJSGW6LeELzvYLd5g2hNzK3s\nol9dmtIipVOydQRBaI/uyNJ5nXB17mjO2EhRZ0sTxDo/1Wkew/Lt15R7uXuGi8IJLl6u8HFFnpPi\nN00+/soPPKS5FDOS9rJ06VKGMpQiiiLuUUQRd3EXr/7qNnK+u4669AmhYztrA6Fsm23VRqh3HIER\nqYFW29xafRwpyyAIQkfoN/Xw21peMFYEHOv8Jee6+cFsN49fnsrLFb6QV54/zBGRj/+D2W4afbDs\nbx8yZ+5cDh48yBzmkEceFVSwhCVUUEEeecxhDoHjVVT+50IcR/YAZuatVW7hmskuLsox/xnGDTZv\nELdsaOQ/327ucMQuSysKgtARlEmgiQ8KCgp0aWlpt1/XKkH8g9nuDkfA9lILQItZr8X/OMgDXykg\ncDycellIISWUUEcdQxnKHOawhjWh444hWWTcv4VvnJfJ6EGO0PXCxdQ02w8H2FplovzLcpP41aUp\nMe8vCIJgoZTaprUuaO+8fhPht0VnI+DWFjexz4697YtZZF92Y8T31rCG5Kw0Jn7nL9QPSY0Qe4D8\nL92IY2A6Hxzx0eALd7RWuYff7fDhCv4XGTdYhUo5rNxpOqxbNjSKTy8IwmmTEMXT2qoNH6tompUy\n2eDTfHdmckRNnf98uzlU9OyZhx/gWh8cWWsmX2VnZ/PNx17iiao80u9cR9NjC0OTr+6570eMLLyP\nh0u9bDsM2w57SXOqFnV7ZmUl4d7m4Y4Zbt6u8nPNZBcrPzD+/GsHAqz8wMN3ZyX33D+WIAj9loQQ\n/LaIVUbBmiL2TpWf2kbNsSbNloN+ahsDoZmxAJflBnDOv49zBygOlzzNNx97iUumT2bDqSYumXI2\ncxe+wuL5l5I0+0by/uV+rshz8s4hP2cNUQxPdYREvvxYgOI3mym6IJmXK3xs3Ofn/Gx/qDNYcq6b\ndw6ZmcCtDo8LgiC0Q8ILfqwyCoUTXKz9xBRHW1XmDZVH9gQcXJTj4LXKALOzHEwcppg4zAnTfsi2\nhf/OE1XD2NDQzN56GH08wOhB4xn8vS1cMjmDBq/m2Y+9vHYgwMxRLr47MxylW5k/0Bzy7KMXW3n8\n8tSIMQXo+Pq7giAIIIIfUUjNXiq5vE6HJjjNykris/om7pmZzLghDlO4zKt5eJuX/KGK8joNDGP2\nKIVPw97jMDXdRPANvlG8U+Xn4W1eZo8yomzl+4MR7YnDHHj8OrTObWtr6UbbUpKOKQhCZ0howW9N\nUKPXn11V5qe8zkzGGjfEjKrOyXWy9hMf5XWanAFQeQr2n9RUnjTZNdedbZYtbPRq80YwSoXsmFSn\nCle89Gp+t8MsnpI/zBHRLqtTgdiCLqtkCYLQGRIiS6c1YmXfQMuCZfYsH+s7d7zSRHmd5qLRDsYM\nMudVnjSTtn51aQovV/j4yVYPL1b4AXA5FVurzFtDxBq2ylTJtBYxt8Tefqw1QY9upyAIQlskdITf\n0Qg5Ovq3FicHmJmZROEEF/+6roF9J2DqiKTQcoWW939ZblJoQNa+nKF1PatUsmXPRB+zEM9eEISu\nkNCC31a6JsQW2PRUxczMJF47EOCiHAdLzjUR9sLxzmAGjw59r+iCZM7PDot8/rDwvaLvHS3ysXx8\n8ewFQegKCS347dGawC451x25EhbhOjypTtXq99qK0DtSJ0c8e0EQuoIIfhu0JrCxxNneCUR/3+J0\nIvTWIn9BEITOIoJ/mkRH69FifDpZNa1ZSCLygiB0BwmZpdPR+vGtZfG0dqy967aXVdPW/QRBELpK\nQkb4HbVW2orIYx2zrrvloL/FIiUdQTx6QRB6koQojxxNT6U3WmvObtzn71Qp5kSkXvvY4K1lniud\nwSoh4w5B6DakPHIbRFsr3bVEYHqqCi2QIlF622zw1vK05yAbvLW93RRBSBgktKJ1i+d03gRkkLVj\nzHOlR3wKgtDziODTuncuE516jsHKydXuzN5uhiAkFAlp6UTTWvZMX10rtl77eN5TTb32tX9yN1yj\nO+4nCELPI4LfBn21OFl3+OOduYb48YLQNxBLpx/Slj/e0eyYznjs0edKBo4gxCcS4fdDLH88lti2\nFo1H2zKxrtGadRN9rkT8ghCfdEv4pZR6ClgEHNZaTw3uGw6sAvKACuCrWutj3XE/4fQ5zzmEnf6T\nnOccAoSj8Sbt5zlvNU3aT4pK4jznEN7yHY+I0i0hB9occJUMHEGIT7orwn8amB+173vARq31BGBj\ncFvoYdobQH3Ld5xSfz1v+Y4DYRH/2H+Ka12jAMXTnoM82VzZIkqf50rnRnd2u0Le1huGIAi9R7f8\nRWqtX1VK5UXtXgzMCf7+DFAC/Ed33E9onfai8HmudJp0gCbtp177mOdKZ6f/JKX+eqY7BzPPlU6K\ncnCecwhTfQPFlxeEfkRP/uVmaq2rgr8fAiTp+gxgt2zqtY8XPEdo1gGSlWKReySDlZMU5Qh1Cikq\niUWuERwKNJPnSGGDt/a07RzpFAQhvjkjf5Vaa62Uilm3QCm1FFgKkJubeyaa06+xLBtvk1mCcUfg\nZMTxZjS7fae4ypmBZd9k4aYKD79squAEgVDEb/fzm7Sfa12jQhF/ZaCJJ5sruSk5hxxHCtB6pyAd\ngSDEBz3511etlMrSWlcppbKAw7FO0lqvAFaAKZ7Wg+1JCOwWDcAUNQAUTHSk0Yxmtdf8Z0jVSdyZ\nPJYU5eBIwMM6Xw0nCDBaJXONO5Op/oE06UDIxz+gm7nRnR0S7CebK809muFHqWeF7m3/tOjoYK8g\nCD1LTwr+WuAG4OfBz7/14L36HV2JivMdaYxRKSQrFRL5c5MGAaY/zcLNTck5ocHVippqBg0axcf+\nk+wInGSbr54UlcTn6gOUDxxMqb+egqTBEUJ+U3IONMMi1wh+2FDOeEcK/5xs3gCi2y1ZO4IQH3RL\nlo5S6llgCzBJKXVAKXUTRugvV0rtAeYFt4UOcrq57Bu8tTznPcQQh5Prk7NJxswS3u6t5wP/Sa5y\nZrB8wKSQDbNs2TIunPYFJlcc4+ykgVzrygQUj364hUs+X0Ddz1dSkDQ41EGA6Yw2e4+S70hjtecw\nOwInWO07wgZvbUS7rYwhQLJ2BCEO6K4snetaOXRZd1w/EbFHxW1F+9HHoqPpRe6RlAcajf0SACeO\n0PkPF/+E4uJiAC6/9DKmvfA7bptyPhmfHqZ04b9zoqqaRx54kIneI0xd9gBXu8Ne/XNeI+QLnSMA\nxXhHCuc5h7DZezTk9b/gOcJz3kM06QDXJ2ediX82QRDaQGbaxin2XPa2Zsf+qumziGP2KpRWdH1n\nyliucmYwzTGI8Y4UnvYcZPEP7g6JPcDRg4d4Z+G/sXHd31l46TxOVFWHju3+2Qpef/C3ofz+ea50\npjkGAjBIJfHjtHyWpIzmLd9xnvNWUx5oCH5Thz6lwJog9D4i+H2Aea50rnWNokn7qQw0hYRzg7eW\nUn89o1VyaOashb2TGKycLEkZzY/T8vnn5FGcU+dn68rnIs4vpBBnlYfHr15CUpWHQgojjj/930+y\nouojXvAcAWB8UhpT1ACadYDKQBN/aq7ieMDHFDWAUn89f20+BCiudWWyyD1Syi0IQhwgpmqcEm3V\nWLnzIXsG0xFs99WzI3CSzd6jXJ+cHfqe1QGc5xzC857qkMWzwVvL7dlT4e/P8usvX0dT1REKKeQO\n7mAxiymmmCKKyCMPgDWsYXj2KL774iq2pA+iWfvZ4K0NZft86DvFft0calMmppT07kADHwZOMc0x\nkEXIwK0gxAMi+HFKdCqjJZTWDNjznEPY4K1lfFJaMNdehSweK4ceFL9v2h+Ri29dc89Z6Zy/bgVb\nFi6lpKqExSwmjzxWshKACioooYS0rJFc9sITDMwfC/6jbPXXc7crl2mOgewInGSaYyA3JefgbdLs\nCJygGi/THAPxBddK3hE4yQZvLVe7MyUlUxB6GRH8OCU6IrZ78/NczpCwX+vKDNW3sSyegqTBgOI5\n7yGAFimV81zpvOo5infCWOb9ppgXvnILxRSHxB6gmGLqqGPWb35N04Qc3vTXkYyiSjfz66Z93J86\nnrd8x0Ozcv81eRQ0a8Y70khWDp7zHmKKGoBTqRZ2kyAIvYN4+HFKeyWOS/31oYFTS8ybtJ+rXCMZ\no5Jp1gGuco3kWlcmd6aMZbByRnQa05yDyPrkMCXfXsZQhlJEUcQ9iihiKEPZ9e2f4t+zn1MEaEaT\njOKAbubloBe/2XuMpz0H+XXTPnYETjLE4WSRO4Mb3dmc6xzEjsDJUKE2QRB6FxH8Pog1iAvwnLc6\nlP/+nLea/YEmVvuOsNp3mCHK5OLb8+ef91Tz1+Zq/v+PS/nTl/8/TlYdZg5zyCOPCipYwhIqqCCP\nPOYwh+NV1WxeeDMn93xGhnJxkWMoAG8GhX6z9yhTHAM4oJspSBrMlKQB/KrpM85zDgkJv/j2ghAf\niKXTB7EGcXcETrawa85zDiHfmwqoiP12f3/8MQ9bFi6lqcpk3KxhDQAllFBHHXdxF3OYE9p/quow\nWxYuJWfL81RlmMXcq/HiBKrwMJJkbnRnc55zCA82fsoB3RwquSC+vSDEDxLh91Gs2vTRdk2OI4Xr\nk7O5PjkrIrK3xL4gaTDTMkaTe+NVEddbwxr8WW5m/eXXeLLcIbG3GHfj1RxPHxDMvBlEJi58wCCS\n+PeU0VztzuQt33EO6GZGq2RTeiGI5ODHHw018MZy8ykkDiL4fZTOLDLygudwyPO/KTmHZOXgez/6\nIbPvuzXIQVh0AAAgAElEQVR0zqCsTP688UUWLlzIeev+i0FZ4ch8yvf/jfz7/o1pjkFc68rkntQ8\nLnAOA2CeMz1UpuE85xAKkgZzf+p4Tmoft5z6iF3+k5KD3w10RqA7cu72lbDhXvMpJA5i6fRDWpZi\nMPV0zk4aGJwNe4gb3dnc/qP7OKq9HHh6DQXrfs/a3GTuT85h6pSB7H35L/zgin9m6o1fIeu+b5KF\nm7OTBrDInQFAslJc6xoV2oZwaeapvoG85KmhCg8PN37G8gGTgMRe5Lyhxojr9CWQNqLz37cEGuDC\nezp2rvcUuAbEvuf0JZGfQmKQGH9t/ZC2RDM6h3+RO4MU5WCeK50T2sc7vuO84z3OjSnZ/GzZA2Tf\neR9/GHCKA7qZt3zHudqdSf2UdEZs28zskWNY5almjCOF57yHSFEOmnSA57zVXOsaxQntY3ljBeMd\nKVzhNqoyz5UeKrn8BefgiOwge/t2+k+GLKn+TluCbe8MrHPtIt1QA55TcElRpEC31olY53hOtX7P\ntBHtdxxC/6P//6X1U6JF3axudRhQXOIK2i1ROfyVgSYeaPiEKjwArPJUm1r2I+EcWwcSWtwkcxw5\njhR+lDqQeu1jiHJynnMIv2/aH2yF5snmSnYETrAjcIIhDldI2K9LziLD4Y6ZoWOv2W9NyurvtBVR\nW51BRQnkzITNwRJHliBvXwmvFsO8hyKF3fqe5xS4bZG8JeYNNeH9ggAi+H2W6IlZ9gqWKcoRU0Sf\nbK4MiX2WckcMrNqjcGuA1764iXX8eU91KDtokXskJ7QPb5NmvCMlQtyjo3qIXCXrzpSxoQ6mv9Oe\nnTN9iRH78vWQPdMIu12kW+ssrG1vVCRvv59E8YIdEfw+il1Q67WPJh3gKmcGycElCa36OXa7xJRA\nCDA+KY0rXOls9h4DNIvcIwFCNXjGqGS8jkHclJzTZvllKzvox2n5oXbEuq/19rHZV0eVLWUz1mIp\n/ZFoO6dmF7x0N1z8Q9j3mhHmq56B1x+CA2/Alx+N7Bhi2S+WqE+6Et59EsbPM7/Hup8gWPTfv7IE\nwlr05EZ3digKf9pzMLQmrSWoOY4Ufpw2ATClk63SC6AoDzRQ6q8PWS03urPJcaSErrXBW8v9qeMZ\n1IYwt7aUof3tw56ymShLH0ZH6Otvg70boKYM6j6FT182Ir/7/6C2zHQGX1vX+vUaamD1DeaNwHoz\nANi1FkbcIwOyQuuI4PcDou2dcKmFQKuLips3gpEkKwXoUI7+Tck55HtTadKBUO17a03bJ5srmZo0\nsFWRPs85hJ3+ky1q58xzpYeKuS1yZyTc0odpI4z4WjZL1nQj+AMzg4K/wYh8bRmkT4YvPQz7t8Lf\nlsDilTBmtnkr+PttkDkd3GlG5PMXmHOzZ5o8rOlL2raPupopJPR9RPD7CG1l5UT75dZ2vfaFsnPs\nRL8RmPPCbwIpKomnPQdDYwH3p44Pee+DosTajj0t01ody2rP9cnZLc6P5fP3ZWIJqrXPc8oMvAJc\neC+kZRgL5oNnjVifNd8c+9LDkJYOKy+ChiNG9G/72HQIn24wPxcXhX3+tBEwd1kr94ny88XqEUTw\n+winY3+0JqhtVeKMddxk6pwVOt7a/RMlYm+NaEG1Wy+XFBmh9p6KPG5R/qI5L2+O2W44Ajq9li89\nks6mZTB0HIw+HxwuyJlfy8TZkf/G1r0vKYoc9LW3SaweQQS/j9CdYtpeZH26kXd/i9g7S7Sgbl8Z\ntl5m3RYpvq4BRtTfXG62rQ6h4YhZGPLTS5bxj11PMOmVTexePhGAoeOh/NPd3H3JXG6542YefGhZ\nxL29p8x37W8Y9jZJ7r2gtNbtn3WGKCgo0KWlpb3djIQlEWfA9iTRFo+1feoIbFkOuRfDvldh3DxY\n8KixbcrXwyaWsRnjy4wYlM2DV26icctEyj/dzTPM5QTmTe+fRhfx1IZlpKVH2jnj5kHuhcYmevXH\nYZtI/Pv+i1Jqm9a6oN3zRPAFCysjx/L2hZ6hZJmZXHX+PeBKMx6+xoj1e+OWsWZvccT5Q5zZLPSt\n4P9YGhJ7i6vyi/jWVct4czkU3Ap1e8NZO0PGwfG9pgM46wrzdpG/wKSAiuj3Lzoq+FI8TQhhVeBM\nVA/+TGD59hcXwRduMmmVm4uhfj9kzqtly4knIs4vpBDla+BPLCKJhhaLy2+pf4KyN0xRut3/ByOn\nwNhLzDHtN59Z001kn7/AdAZvPdqyuJpUz0wM5L1dCJHoHnxPEp1Fk7/ARPb7Npvju/4GjbXpfPfi\nTTzwqrFt2ltcPmtUNiWbN/HxL9I5+ibU7zNjAjNvhSSXyegZng+Trzb3tlI4D7xhjn3yMoy5EM67\nTTJ4EoUej/CVUvOVUruUUuVKqe/19P2EziG16g09HeFagrr/DWOxlK83Ns759xhRbqw1ncA1v5jI\nD87fxBBXNiWUhFYfW8nK0KpkJZSExH64nsiuv5l7OFPN5+Gd5toAR8uNj7/hXjMxC4zYD8s3cwFe\nLQ57+9ElHYT+R49G+EqpJOAx4HLgAPCOUmqt1vqjnryv0HESZbZre5xuhNvaZKboCpieU0bo924w\nIg/gbYAv3mt+tq+E3ItM3n1j2UQWsoI/sajVxeX/57//h9HDJ/LURaazcLjgrC9DUy34vFC52QwK\nj5sLU68zkb33VLhDy/8SpH0tMqtHIvv+T09H+LOAcq31p1prD/AcsLiH7yl0AvHtDacb4ba2kIi1\nf/UNxjO3JkNdcA8c2m6Ef8tyc15DrbFX/vdaM9v2xJDdrHctbXNx+RuvX8pfi3dTW2bEPuCFXc+b\ndM/KN825lti/dDf4GsxYQe0ucyxtBMxZZiZtgfj3iUJPe/g5wH7b9gHgPPsJSqmlwFKA3NzcHm6O\nEI349obTjXCtDmLSlUY0J11prJNJV4br3OTMDA+YghH73Etg9CwTdVu1dQBq2M0zx8MevmXj2D38\nOcxhTf0abv+vudw6ehPuAxNJGQZTr4fzbjf305jc/79eZywc7ynToVnts3ds4t8nDr0+aKu1XgGs\nAJOW2cvNEYROYXUUbywP17S3hP2qZ0x0D6Yypt8DjcfN9uhZprzChnvNIOvBUjhWV8sfHHM5ETAW\nW3uLyx/3HuS3B+byLd6HY+kMHQsjJpnI3SJzuhH87FlhMR8RJerTlxjLybJ8JGWz/9LTgl8JjLFt\njw7uE4R+hT3Sz5sT9sXdA4yoV74TjuIBnGnmnOOfwc7noLkO0kjn/ME381JdOA9/DWsY6s7m7nP+\nhyfeW9picfkvcDPZ49IZfT68u8KMA4yZHT7+xXthQEbbVpW9na4BEuX3Z3pa8N8BJiilxmGE/lrg\n+h6+pyCcceyWkD2CtncEOTPNQK0zzaRCpo2ATzeaQdeUYWY27Pnly0jOhbX7jOgPIpvvTNmEf/tE\nbmATf1BzqdfmDaBwXBH/cvYyytdDc725zuqvw7f3mHt2pjqmvQSDVNXsv/So4GutfUqp24CXgCTg\nKa31hz15T0GIJ+wdgWW12AV18cpwGeT0fDPIy/plHAd2DnyC335rEyOYyJvbYQQTuSVjE4/XzOWS\n7Jv5zZ+WUf6i6Uhq98DOP5nMHAtrsNhzKjw425F2WvYUSLTf3+hxD19rvR5Y39P3EYS+QvQg6W0f\nm99rdkHAAznnw9wty7j9xtsZnpaOt8Fk93y8Go6VT+TfeZ8v35TOvtfCa90OCq5WmZYe7lC8DWaf\n6mT7pKpm/0VKKwjCaXK6k7XsKaD2a1g1790DTPXMU2XpvFps0jddaXD2VaY+Thrpofx56zruNHNt\nd1q4Q3GlmeOzbutc+6xoX+yc/kevZ+kIQl+lo+mM9vVnrZTIWPbJlx42Ef7wSWZG7t4NZrLWqOmw\nL7h9SZEZWI2eLDXrtvB+C/HghWhE8AXhNGioMd74JUXtWx9Wx2BP2bSEetKVZtLVqSPGjhl/RbgD\nGD8PRl9oft+7IVxXP5aIR88jEO9diIUIviCcBttXhv3zjmbA2FM2LT541oj53g1hW+b8e8zvVslk\na3GUWF68ZNQInUE8fEE4DTpTisGKvkdMaumNWyI+fp4R+M3FJm9+1m1m0HXcPFMewT3AHGuthEP0\nfkGIhUT4gnAadFexMct7txY0tyyi7SvNYC1ElkKIzpOPLu0gkb7QFiL4gtCL2Esz2C2iWGvUtpYn\nby/tYO0ThFiI4AtCHBCd+25Vs2zvXHsGkP2YIMRCBF8Q4oDTtYik0qXQGUTwBaGPYRd5mRUrdAbJ\n0hGEOKGjM3ftGULtzYqVxckFOyL4ghAndDTFsjOlDyRtU7Ajlo4gxAk9Yc+I5SPYEcEXhDihJxYS\nl8XJBTti6QiCICQIIviCIAgJggi+IAhCgiCCLwiCkCCI4AuCICQIIviCIAgJggi+IAhCgiCCLwiC\nkCCI4AuCICQIXRJ8pdRXlFIfKqUCSqmCqGPfV0qVK6V2KaW+1LVmCoIgCF2lqxH+TuBq4FX7TqXU\nOcC1wBRgPvC4Uiqpi/cSBEHoNjy6kU882/Hoxt5uyhmjS4Kvtf5Ya70rxqHFwHNa62at9V6gHJjV\nlXsJgiB0J/u9ZXzs2cp+b1lvN+WM0VPF03KArbbtA8F9giAIccEY1+SIz0SgXcFXSm0ARsU4dL/W\n+m9dbYBSaimwFCA3N7erlxMEQegQbpXKWe7pvd2MM0q7gq+1nnca160Exti2Rwf3xbr+CmAFQEFB\ngT6NewmCIHQJj25kv7eMMa7JuFVqr12jp+mptMy1wLVKqWSl1DhgAvB2D91LEAShU0QP2Fp+/ntN\nr7Q6iNveIG9fGBPoalrmVUqpA8D5wDql1EsAWusPgT8DHwEvArdqrf1dbawgCEJ7WMJ8MnCsVYGO\nFucxrsmMTMrlsH8f+71leHQju5rfZlfzOy06Bbug2zuBMa7JnO2eHddjAl0atNVarwZWt3LsQeDB\nrlxfEAQhFm3ZJ5Yw1/oPcti/D6CFVz/GNRmf9tEcaOSjpjdJUi7OSb6AdF82Y1yT2e8tY493W+j8\nSckzYw7yWvey7hHvYwKyxKEgCH2OaKG1k+nMo9Z/kHz3DNL92RECbe8onMrJHm9p6JhTOTnLPR2P\nbsSnfQxXoziqDwFmaNEa5LWi+jGuya1m+sSrny+CLwhCn6OtlMpqXwWH/ftI92e36AzsHYWJ8r14\ndBMnAkdpDjSGhHqPt5QJrgLSyQEUHt0YEu6ORPVtdUi9iQi+IAh9jrZSKmN1BpaQZzrzQsfcKpVx\n7nN5r+kVjgaqOBqoItmRGvH9vZ4P2OMt5Zi/ii+kXo5bRR5vLZKP1xx/KZ4mCEJc09kSCFZn4Fap\nocHXdxs38LFnK5XePfi0j72eD0Jifdi/j6GOTNLUYIYlZYVEfb+3DD8+AGoClTGzb1rLzLG3oSvP\n0t1IhB9P1NTAypWwZAmMGNHbrRGEuKAr9oh98HW4IwtQEb59jmsCtf6D+AJeGnQ95Z5tzEpdyF7P\nTvZ4SxnuyGK8cxpJyhmK1u3tscYLrDcH6NiA8uk8S3cgEX48sXIl3Huv+RQEATCiOjIpN0JUo2kt\nch7jmkwyaQCcChwnx5UfFH6o9R2k0ruH8sMfMTRpJCOTcjkn+YLgNzX1tSc5GqgiSZm4eK9nJycD\nx/BpHxNcMxjjmhwaL6j0lofu31Y+fm+nbkqEH08sWRL5KQhCeBDWl81A97AWxz26kfeaXomZgulW\nqThUEmhopoFqXwXpSTnGs9dV/OnH6/jLyr/xzEvj+OI5l1Hp3QMojn/q5Y7Lfs4/LZnHN+7LpE5X\nA1Drr+RooIqxzikRYwI+7Y0YDLZ/2untcg4S4XeFmhpYvtx8dgcjRsA990TaOd19D0GIc6Kj9fai\nYsuHH5mUG/OcqckX4SKFLHUWPu0jx5XPBFcBL/78XR578L85fPAIX7vimzz//kr2eLdR8tE6Fl52\nJTVVR1n50z/z+INPAjDCkQPapGge8e/nY89W3ml8kUxnHuPc54ba2Jp/35ln7ilE8LvCmbBgxOYR\nEoxoS6Q9AbU6hM+nXArAruZ32NX8dkg8TwaO4qUJv8PLHm8p1b4Knv3ZOh578L9D1zhadZy7FjxA\n2cuV/GjR49RUHQ0de/Zn63nx5+/yhdTLOSv5C7hIYaJ7FgPUUE7pOj5qfrPFQHFnxftMlWUQS6cr\n9KQFYw3gXnllx+5hP3/tWhn4FfosHZnMBEQMjFpZNT7tDQ3SHg8c4fMpl4bOz3Tmke7Lxnl8ML9/\n4vGIaxdSSElVCff884MMZSiFFLKGNaHjq1eup/iun/NZ2k68NHHQt5uZqfP5qPlNm+9PqF0fe7bi\n0z6cwcHe9qL9M5XGKYLfFSwLxqK1LJvTyb6xInuIvEd755eUwPr1Hf+eIMQZrfnc9gwXICLbxTo2\nwVXABFcBx/xVHPbvY69nZ4Toul0pvDfwFYpfuIUfLXqcwwePUEghd3AHi1lMMcUUUUQeeQCsYQ2Z\n2Rls2rSJ9PR0kgMXQDOck3wBAx3DmJW6sEU7rbINx/xV1AQqOyT8Z8rbF8HvTloT6c6KN3T+7cE6\n78orYc4cGfgV+h2xouBMZx6feLZHTKjy6CaO+asY75wG6FBdnc+nXBry+6dPnskzLz3O1674JiVV\nJSxmMXnksRJjnVZQQQklZGZn8Pw/VlGZs40hvmROBo5yTvIFocHdce6pMUX8mP8QNYFKRjhGh9oA\nvT/rVmkdPyXoCwoKdGlpafsnxittRfiPPmp+v+46Y7lY1otYMIJw2nzi2c7Hnq2c7Z4dEtOtDWup\nCVSSpgYzK3UBHzW/yWH/PgaooZyTfAGfeXcy1jWVvZ4dbH9pN7dcfU+E2AMsYQkVVPD71b8k//KR\nNOoTOHHjw0OqGkSjPgEQcd/oNgFMcM1gnPvcHq+ro5TaprUuaO88ifDPFG+8ARs2hD8t66WrFoxM\n1hISmFhR/wA1lBoqadD1VHr3MMQxkpP+utAA6yldx6lAPbt37+KHt/6GoQyliKKI6xZRxF3cxX23\nPMBv1j3AkAkuHCQB0KhPMFRl4lAOGvwn2NX8DuPcUwGTq+/RjQxXWQxyDAcU0PuRvYUIfnfSlqWz\nYYP5fdIkcLvhhz801kssC8Z6I2hogLQ0uO221sX8dOwiQegnxPK+3Y5U8MMIx2j8+PjUu4OhypRO\nGOf+HHs971O95yj3L/wNtVXHKKSQPPKooCLCw5/DHNZUreHuRT9m5Yu/xTHOG7qHA2Vy+akCv6m0\nCURW38TFZ94POeavYlhSVqv2z5lEBN/idCNl+/fsvrs9a+bUKSPGaWlGxNevh5kzYdkyc360UK9c\nCcXFse9nWULWvU6dgqKi8LbdPrrttvD15A1A6Ad0pOzwOPfU0CDpXs8HANTpakY6cqn1VfLpkTJu\n//KDHK06DhDKximhhDrquIu7jNgH9x8+WMONX7qVX235Hunp6XhpBqUY65jC8UCNmdSlG8l1nY1P\n+/BrL0nKSY5rAjTDYf8+agKVofLLvYkIvsXpRsrR37O+u2yZEe1162DzZrN/wAAj+AB//KMR70mT\nIjuHJ5+E11+H88+HL3whLNJWB/DOO2EraOZMs/+hh8Ln2TuLAQPMp7wBCP0E++Imn0+5NCT6JwPH\nQimSAx3DQsI6zn0uoEJZOwFHgMHpAylcMp+nfroqdN01rCE9axg/+s23eOzbz7Kmak3EfZd8cwn5\nI88h3z2Dcs82Dvv34UxyhWbg1nmrSVapTEqeGfG9c5IvINDkZ3DSiLionCmCbxErOu9IVBwrm6am\nxnj1EJqZx9tvh4V/8mQoK4O776b2D38g3eo0fvc72LsXgFog/cgRuOoquOmm8HXffRfGjQu/JRQV\nmSi/psa0dckSsw2mA3n2WXPOlVeaGbsS6Qt9DI9uZK9nJ6BDxc6spQgtYbcGZmkmIlXSHRRhj26k\n3LOdOl81Y5OmcPP9OTTqEzz7MzN+NjxrCA+u+zZ5E3LJXjcyZPcAXPf9Bdxw39UkO1Jxq2SGODIY\n4hhJjiufId6ReAKNnNLHY9b6qfZVUBOoJMM5ptftHBDBD2PPqV++vONRcXQuPhhLZcMGuPhiOO88\nmDvXRPabNxtbZ80auPtulg0bxhM5OWz6wx+YOG9eyOffnZzMXJ+Pm8vLWbZ8ublmRobpRLZsMdvn\nnw+bNoFS5rrvvAMPP2zsHsvzX748/Aawdm34mZYsEZtH6DNYC5IAOJWLz6dcyl7PzuDiJWZhknOS\nw/nxsXCrVE4GjnFUH6I50MQpXce3f/BvDHdk8/zKdTz10m/InZDNUW8VORMyeXDdt/nRwsf4pyXz\nKPz+xdQHaqnxHWCft4xTuo6RSbmMU1MZ557Ke02vUBM4QLWvokWtn3iriy+CH4vumkHrdhvRXbDA\nDNJ++CHMnw9r17IsK4viJ02Njrlf/zqbvv51JgK73W7mNjdzELBc/GXW28Gtt0JFBZSXQ3U1fPqp\nOSEjI5zpY8/4saL9U6eMfWQ9kwz0Cn2ITGce1d4KUIpMZx5ulYpTOfnYsxWncnGWe3qLSVCxvP58\n9wxONB4jSbsY75pGvns6s39yJcV31ZKens6u5rep4zAAORMy+f3WH+MarnGSTIpjICPIoSZQiYuU\n4KSuDzgeOBJK+YwV4fd2sbRoRPBjEStq7wzWYKnl169fDx6PieA9HpZt2IB9SPagx8PcP/+ZFcBS\nj4eDtmPFAF4vy6x2bdliBPuzz+Cxx2D8eCP8CxaYCH/mzEiLZ8AAI+4DBoSfSapyCn2Ial+FWVtW\nE4qi21t1Klbd+WP+Kho5QSMncPuTQ+emp6fj0Y0c8x8CYKjKxKlcTMyZyY6mTZzSdRzwlTHWOYVG\nfYpTuo4RjhyO+aupCRwI1dSJFeHHGyL4PYEltMXFxj+fPh3++lcAaidP5omSEvD5QqcXUkjJiRIW\nQbCOx5yIOh5PfPABt99zD+lWRwJw++0wdmzLiVvRAi/iLvRxrLVnQYWE3h452yc6WftidQjDkrJI\nZSCNnGSwIz10fau8spkZm8MAx1BOBeoAyEwaS5XPRyMnORWoC9k5Qxwj2eMtDdXQr/ZVxI1t0xYi\n+N1FTY3xyrdvNx6+XWivu85E4enppN92G5suuoi511/PQb+/3Toe2YMHs6m+nvSMjLAvb9XMeeaZ\n8ECtfSavZeNYUX6stE+xdIQ+ghl4ndXq8WifPDritzqEkUm5NHKSkUm55CeHbRar3MIANZTBSSP4\n1LsDgIamjTToesY7p5HsSDXWkk3Y7fVx4j2ytxDB7w5qauCGG8L+eWEhvPaaEdNdu+CTT8z+2lpY\nu5aJwCa/n7nJyZQ0t17HI3vwYDatWsXEF16AFSsgNxdefhkuucTcyxp4td97wIDYNo4difqFfkS0\nTx5t50RXy7S2P/FsZ4xrMmNck0OZP0mcxVjnFI749jNcZdGg60my5c/bhT2evPmO0iXBV0otB/4J\n8ACfAEu01nXBY98HbgL8wLe11i91sa09S1dKFKxcaQT34ovh4EGTcvnoo0ZwV682qZa5ueHJWbW1\nTHzqKVaUlbGIZoopjqjjUUwxddTxP/Uw8YMPYONGM1D79a+D1wvz5pm3CSuyX7/edAL2GbvRUb6d\nro5RCEIcEx3x2zsEtyslWEbZF8r8Ocs9PVRYzZqs1UA9oxzjONsZufCK9fZgj/bjId2yo3Q1wv8H\n8H2ttU8p9Qvg+8B/KKXOAa4FpgDZwAal1ESttb+L9+s5OmJz2CdIPfus2XfbbeFsGGuQNs2socm9\n95qceYATJ4zdMmIEPPoou8vKWOp2M9ST1modj6XuBjaNHcvEyy6DI0fMG0JGBnznO/DBB5FtmzMn\nPHMX2o/ypQaP0E+xC3y0vRMuozyjjVW0TP2bpGAGkB37xK9YSyrGO10SfK31y7bNrcC/BH9fDDyn\ntW4G9iqlyoFZwJau3K9H6YjN8eijZiD25ZfDtXEsQbUGaSGcMQNGWJ97zoj13XfDM8+w++WXmYvJ\nzilkQet1PDxrmPu1r7HJ52NiUVG48Novf2k+rQlWRUXhzCCLK680Pr+1gEo0bXVw0hkIcU5HSixA\n6/ZO9Pfsi5aAZoKrIFQQzX5PawHzHNeECHvodNrWG3Snh/8NwJqrnIPpACwOBPfFL52xOaZPhwsv\nNBG9ZZssWWKi8NWrw376zJkmdfLWW+Gll2DcOGofeoi5W7aEUi/bq+Nx0OdjblIS7x84QLrXG77/\nFVeY+y1fbgQfImfSrl0bno07YEBL8W6rg5NBXSHOiZV2GYu27J1Y51krZp3tnh1artASb2sCmDmW\n0uW29Qpa6zZ/gA3Azhg/i23n3A+sJlxf/1HgX23HnwT+pZXrLwVKgdLc3Fwd1xw5ovVDD5lPrbUu\nKtIazKfW5hhoPXmy+Zw3zxy75x6zHTy3KD9fAxE/2UlJ+gXQ2Uq1OFZkfRe0zs/XuqzMtGHevMj7\ngNYLFphjR46YfdY5Dz10+s8pCHFGc6BBlze/q5sDDZ061tnvlje/q//vxOOhffZPa3/099q7f08A\nlOp2tFwbJWj/pDYvADdirJo0277vY7x9a/sl4Pz2rjVjxowe/UfpduyCbwlsUZER5AULwkJrnTdv\nXkiMiy64ICz2TqfeFRT0XaCzk5PDYj96tNbjx4cF3xJ165r5+ebznnsi76l1uAOyOgFBSACixVhr\nHVOs32p4oV1Rbk282+oYeoOOCn5Xs3TmA/cCl2itG2yH1gJ/Uko9jBm0nQC83ZV7xSXXXWdq2Myf\nH06NLCoydop91ut117WwVZYVFsKbb/JEcjKbmpuZmJ8P5eVMzM9n029/y9zCQm5ubmbZgQPGVtm+\n3VTW3LUrPHO3qMjYSsuXh/P/7Zk69k/x4oUEIVb9GrvNYk/DtBdgi0VrFlD0/nirmdMqHekVWvsB\nyoH9wHvBn9/bjt2PSdXcBXy5I9fr9Qg/2spozdqItkusyNoeeT/0UDjCjmWnBN8CalatMte59dbw\nm6LqMxUAAAXoSURBVILWumbLFq0vvljrSy4x51ps2aJ1RkbYOtqyJTKyFztGEFoQHZH3hu3Sk3Cm\nLJ3u/Ol1wY8W6NYE29pviXxZWVhk7YJr/W4/Hn0NS6wtW8Z+nv3+1rWsTsblirRr7OeIZy8ICYUI\n/unQ2QjfFpG3SSwRtncGlohbXr3N6w/d33pzuPXW8KDw5Mnh6L+tzuV02icIQp9BBL+76Gj029Z5\n7XUkR45ERvqxxDd6gDj6XtEZQz3xjIIgxCUdFXyppdMeHc1Jb+u86Bz/6HNHjDCF0OwLl0fnx992\nW+TAb3fmx0upBUFICETw26OjhcZam9kaa9ZqrGuOGGFKI1jnW9+1L0jelijbOwRBEIQYiOC3R6zo\nN1rEa2pM2YT1601apHV+dBVN+9KC9mvar2eP/iFyQfK2BF+idEEQ2kEE/3SItmSsapkLFkRG2NH7\nW7N9rBo9p06Fa+LYq17atwVBEE4TEfzTIdbEJuszul6NJdixzo9FdKRur4DZEdoqfCZF0QQhoRHB\nPx2iRbk1OyXWkoOxzovlv5+uOLc1eCxF0QQhoRHB72lOJ6qH0xfn6PvZOw5Z6UoQEhqrumVcUFBQ\noEtLS3u7GfFBd9kv1hq4Dz0kUb0g9FOUUtu01gXtnScRfrzSXVk3EtULghBEBL+/I+magiAEcfR2\nAwRBEIQzgwi+IAhCgiCCLwiCkCCI4AuCICQIIviCIAgJggi+IAhCgiCCLwiCkCCI4AuCICQIcVVa\nQSl1BPist9vRBUYANb3diG6mvz1Tf3se6H/PJM/TecZqrTPaOymuBL+vo5Qq7Ug9i75Ef3um/vY8\n0P+eSZ6n5xBLRxAEIUEQwRcEQUgQRPC7lxW93YAeoL89U397Huh/zyTP00OIhy8IgpAgSIQvCIKQ\nIIjgdwNKqeVKqTKl1PtKqdVKqaG2Y99XSpUrpXYppb7Um+3sKEqpryilPlRKBZRSBVHH+tzzWCil\n5gfbXa6U+l5vt6ezKKWeUkodVkrttO0brpT6h1JqT/BzWG+2sTMopcYopTYppT4K/v92R3B/X36m\nFKXU20qpHcFnKg7uj4tnEsHvHv4BTNVafw7YDXwfQCl1DnAtMAWYDzyulErqtVZ2nJ3A1cCr9p19\n+HkItvMx4MvAOcB1wefpSzyN+Xe38z1go9Z6ArAxuN1X8AHf0VqfA8wGbg3+N+nLz9QMXKq1ngZ8\nHpivlJpNnDyTCH43oLV+WWvtC25uBUYHf18MPKe1btZa7wXKgVm90cbOoLX+WGu9K8ahPvk8QWYB\n5VrrT7XWHuA5zPP0GbTWrwJHo3YvBp4J/v4MUHhGG9UFtNZVWut3g7+fAD4Gcujbz6S11ieDm67g\njyZOnkkEv/v5BvD34O85wH7bsQPBfX2Vvvw8fbntbZGpta4K/n4IyOzNxpwuSqk8YDrwFn38mZRS\nSUqp94DDwD+01nHzTLKmbQdRSm0ARsU4dL/W+m/Bc+7HvKb+8Uy27XToyPMIfQuttVZK9bm0O6XU\nQOCvwJ1a63qlVOhYX3wmrbUf+HxwLG+1Umpq1PFeeyYR/A6itZ7X1nGl1I3AIuAyHc51rQTG2E4b\nHdzX67T3PK0Qt8/TAfpy29uiWimVpbWuUkplYaLKPoNSyoUR+z9qrZ8P7u7Tz2Shta5TSm3CjLvE\nxTOJpdMNKKXmA/cCV2qtG2yH1gLXKqWSlVLjgAnA273Rxm6iLz/PO8AEpdQ4pZQbM/i8tpfb1B2s\nBW4I/n4D0GfezpQJ5Z8EPtZaP2w71JefKcPK0lNKpQKXA2XEyzNpreWniz+Ywcv9wHvBn9/bjt0P\nfALsAr7c223t4PNchfG4m4Fq4KW+/Dy2ti/AZFF9grGuer1NnWz/s0AV4A3+97kJSMdkfewBNgDD\ne7udnXieL2IGNN+3/e0s6OPP9Dlge/CZdgI/Cu6Pi2eSmbaCIAgJglg6giAICYIIviAIQoIggi8I\ngpAgiOALgiAkCCL4giAICYIIviAIQoIggi8IgpAgiOALgiAkCP8PUyM1UStP2x8AAAAASUVORK5C\nYII=\n",
|
|
"text/plain": [
|
|
"<matplotlib.figure.Figure at 0x7f500506f438>"
|
|
]
|
|
},
|
|
"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",
|
|
"\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": 34,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from numpy import exp, sqrt, array"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 35,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def gaussian(d, bw): return exp(-0.5*((d/bw))**2) / (bw * math.sqrt(2*math.pi))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
" This person at the science march certainly remembered!\n",
|
|
"\n",
|
|
"<img src=\"images/normal.jpg\" width=400>\n",
|
|
"\n",
|
|
"Since all of our distances are positive, we'll only be using the right-hand side of the gaussian. Here's what that looks like for a couple of different choices of bandwidth (bw)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 36,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW5+PHPM5M9kBCyEJIQAiTsOxFEEEWRVa/rbd3r\ntS0Xt9uKdentru2vtl4V11qkWm211qooLoC4ggpC2GSHAAESlkDYCdmf3x9ngEiBDDCTk8w879fr\nvM6cbeY5UZ7vOd/zPd+vqCrGGGPCh8ftAIwxxjQuS/zGGBNmLPEbY0yYscRvjDFhxhK/McaEGUv8\nxhgTZizxG2NMmLHEb4wxYcYSvzHGhJkItwM4kZSUFM3JyXE7DGOMaTYWLly4S1VT/dm3SSb+nJwc\nCgoK3A7DGGOaDRHZ5O++VtVjjDFhxhK/McaEGUv8xhgTZppkHb8xJjxVV1dTXFxMRUWF26E0WTEx\nMWRlZREZGXnG3+FX4heR0cATgBeYoqoPn2S/c4C5wLWq+sbpHGuMMcXFxbRs2ZKcnBxExO1wmhxV\npaysjOLiYjp06HDG39NgVY+IeIFngDFAd+A6Eel+kv3+AHx4uscaYwxARUUFycnJlvRPQkRITk4+\n6zsif+r4BwKFqrpBVauA14DLT7DfXcCbQOkZHGuMMQCW9BsQiL+PP1U9mcCWesvFwKDjAskErgSG\nA+eczrGB9OTH60huEUXnNi3JS2tBq7ioYP2UMcY0W4F6uDsJuF9V6860NBKR8cB4gOzs7NM+vqa2\njilzNrC/ouboupQW0XRu04K8tBacl5vCJd3a4PHY1YQx5sSKioq49NJLWb58eUC/t6ysjGuuuYYF\nCxZwyy238PTTTwf0+0+XP4m/BGhXbznLt66+fOA1X9JPAcaKSI2fxwKgqpOByQD5+fmnPQJ8hNfD\nkl+OZNv+CtbuOEDhjoOsKz3A2h0HeXNRCS/N3UTnNi2486I8xvVqi9cKAGNMI4mJieGhhx5i+fLl\nAS9UzoQ/dfwLgDwR6SAiUcC1wLT6O6hqB1XNUdUc4A3gdlV9259jA8njETJbxTK8Sxo/HNaRP17T\nh7fvGMLSX43kiWv7ogr/84/FXPL450xdXExNbV2wQjHGNFM1NTXccMMNdOvWjWuuuYbPP/+cq666\nCoB33nmH2NhYqqqqqKiooGPHjn59Z3x8PEOHDiUmJiaYofutwSt+Va0RkTuBmThNMl9Q1RUiMsG3\n/bnTPTYwofvP6xEu75vJZb0zmLFiO09+vI67/7mUSR+t48cj8riib6Y9UDKmifnNuytYuXV/QL+z\ne0YCv7qsxyn3WbNmDX/5y18YMmQIt956KwsWLGDJkiUAzJkzh549e7JgwQJqamoYNMh5ZPnII4/w\nyiuv/Nt3DRs2jCeffDKg5xAIftXxq+oHwAfHrTthwlfVWxo61i0ejzC2V1tG90jno1U7ePITpwD4\nqrCMh67oSUyk1+0QjTEua9euHUOGDAHgxhtv5Mknn6RTp06sWrWK+fPnM3HiRGbPnk1tbS3nn38+\nAPfeey/33nuvm2GflrB8c9fjEUb2SGdEtzZM+mgtT35SyNrSgzx3Y3/aJsa6HZ4xBhq8Mg+W4+/+\nRYRhw4Yxffp0IiMjGTFiBLfccgu1tbU88sgjQIhe8Ycqj0eYOLILPTITmfjPJVz21Bc8e8MABnZo\n7XZoxhiXbN68mblz5zJ48GBeffVVhg4dyoABA7j55pu5+eabSU1NpaysjB07dtCzZ0+g+V3xWydt\nwKge6bx9xxBaxkRy/fPz+NvcIlRPu2GRMSYEdOnShWeeeYZu3bqxZ88ebrvtNgYNGsSOHTsYNmwY\nAL1796ZXr16n9WwwJyeHiRMn8te//pWsrCxWrlwZrFNokDTFBJefn69uDMSy73A1d/9zCZ+sLuU7\n+Vn87speRHqtbDSmsaxatYpu3bq5HUaTd6K/k4gsVNV8f463rFZPYmwkU27O587hubxeUMx9b3xD\nXV3TKxiNMeZshHUd/4l4PMJPRnUhOsLDo7PWkhQXxS8u7WbNPY0xIcMS/0nceVEuu8ureOHLjSS3\niOKO4bluh2SMMQFhif8kRIRfjOvO3vJqHpm5hlZxkdwwqL3bYRljzFmzxH8KHo/wx2t6s+9wNT9/\nezlJcVGM7dXW7bCMMeas2MPdBkR6PTxzfX8GZCfxo9cW88W6XW6HZIwxZ8USvx9io7z85ZZz6JTa\ngvF/K2B5yT63QzLGBEFRUdHRl7ICadasWQwYMIBevXoxYMAAPvnkkxPu9+tf/5rMzEz69u1L3759\n+eCD4PR2Y4nfT4mxkbx860ASYyO5/ZVF7K+odjskY0wzkZKSwrvvvsuyZct46aWXuOmmm0667913\n382SJUtYsmQJY8eODUo8lvhPQ1pCDE9d14+SvYf56ZvL7O1eY0JQMLpl7tevHxkZGQD06NGDw4cP\nU1lZGbRzaIg93D1N+TmtuWdkZ/44Yw3nfp3MTedaSx9jgmL6A7B9WWC/M70XjHn4lLsEu1vmN998\nk/79+xMdHX3C33/qqad4+eWXyc/P59FHHyUpKelMzvSULPGfgQnDOvH1ht089N5K+me3okdGotsh\nGWMCJJjdMq9YsYL777+fDz/88ITbb7vtNn7xi184zcl/8QvuueceXnjhhcCdnI8l/jPg8QiPfacP\nY5+cw52vLubdu4bSItr+lMYEVANX5sESrG6Zi4uLufLKK3n55Zfp1KnTCX+7TZs2Rz//8Ic/5NJL\nLw3UaX2LX3X8IjJaRNaISKGIPHCC7ZeLyDciskRECkRkaL1tRSKy7Mi2QAbvpuQW0Tx5bT82lR3i\nf9+y+n5jQsWRbpmBo90yn3/++UyaNInBgwcf7ZZ5zZo13+qW+cgD2frTkaS/d+9exo0bx8MPP3z0\nbuJEtm3bdvTz1KlTg9LCCPxI/CLiBZ4BxgDdgetEpPtxu30M9FHVvsCtwJTjtg9X1b7+9hzXXAzq\nmMzdIzozbelWXluwxe1wjDEBEIxumZ9++mkKCwt58MEHjzbVLC0tBeAHP/gBR3ojvu++++jVqxe9\ne/fm008/5fHHHw/KOTbYLbOIDAZ+raqjfMs/BVDV359i/xdUtZtvuQjIV1W/33xyq1vmM1Fbp3zv\nhfksKNrN23cMoVvbBLdDMqbZsm6Z/dMY3TJnAvUvZ4t9647/0StFZDXwPs5V/xEKfCQiC0Vk/Ml+\nRETG+6qJCnbu3OlP7E2C1yM8/t2+tIyJZOLrS6murXM7JGOMOaWAteNX1amq2hW4Anio3qahviqg\nMcAdIjLsJMdPVtV8Vc1PTU0NVFiNIrVlNL+9oiertu3n+Tkb3A7HGGNOyZ/EXwK0q7ec5Vt3Qqo6\nG+goIim+5RLfvBSYCgw842ibsNE90xnTM51JH61jw86DbodjTLNlDSVOLRB/H38S/wIgT0Q6iEgU\ncC0wrf4OIpIrvqccItIfiAbKRCReRFr61scDI4HlZx11E/Wb/+hBTISHB95aZiN3GXMGYmJiKCsr\ns+R/EqpKWVkZMTExZ/U9DTY+V9UaEbkTmAl4cR7crhCRCb7tzwFXAzeLSDVwGPiuqqqItAGm+sqE\nCOBVVZ1xVhE3YWkJMfxsXDfuf3MZry3YwvWDst0OyZhmJSsri+LiYprTc77GFhMTQ1ZW1ll9hw22\nHmCqyvXPf83ykn3MmngB6YlnVzIbY4w/bLB1F4kIv7+qF1W1dfzineV2y2qMaXIs8QdBTko8Ey/p\nzKyVO5i+fLvb4RhjzLdY4g+S7w/tQM/MBH75zgr2lle5HY4xxhxliT9IIrwe/nB1b/aUV/H/Pljl\ndjjGGHOUJf4g6pGRyA+GduD1gmKWbNnrdjjGGANY4g+6uy7OI7VlNL+etsLa9htjmgRL/EHWIjqC\n+0Z1YcmWvbyz9KQvPBtjTKOxxN8Iru6fRe+sRB6evppDlTVuh2OMCXOW+BuBxyP86rIe7NhfybOf\nFbodjjEmzFnibyQD2idxZb9Mnp+zkc1l5W6HY4wJY5b4G9H9o7viFbHmncYYV1nib0TpiTHcMbwT\nM1Zs56tCvwckM8aYgLLE38h+cH5HspJi+c27K6mx0bqMMS6wxN/IYiK9/HxcN9bsOMA/5m92Oxxj\nTBiyxO+CUT3SGdwxmUdnrWXf4Wq3wzHGhBlL/C4QEX5+aTf2llfz3Ofr3Q7HGBNm/Er8IjJaRNaI\nSKGIPHCC7ZeLyDciskRECkRkqL/HhqseGYlc0TeDF77YyLZ9h90OxxgTRhpM/CLiBZ4BxgDdgetE\npPtxu30M9FHVvsCtwJTTODZs3TOyC6owadY6t0MxxoQRf674BwKFqrpBVauA14DL6++gqgf12FBT\n8YD6e2w4a9c6jhvPbc+/Fm5h3Y4DbodjjAkT/iT+TGBLveVi37pvEZErRWQ18D7OVb/fx/qOH++r\nJioIp4GW77wol7ioCP44c43boRhjwkTAHu6q6lRV7QpcATx0BsdPVtV8Vc1PTU0NVFhNXuv4KCZc\n0JFZK3dQULTb7XCMMWHAn8RfArSrt5zlW3dCqjob6CgiKad7bLi6dWgH0lpG8/vpq21wdmNM0PmT\n+BcAeSLSQUSigGuBafV3EJFcERHf5/5ANFDmz7EG4qIi+PGIzizctIdZK3e4HY4xJsQ1mPhVtQa4\nE5gJrAJeV9UVIjJBRCb4drsaWC4iS3Ba8XxXHSc8Nhgn0tx9Jz+Ljinx/HHmGuvKwRgTVNIUqxby\n8/O1oKDA7TAa3Yzl25jw90X84epefPecbLfDMcY0IyKyUFXz/dnX3txtQkb1SKdvu1Y8Nmsth6tq\n3Q7HGBOiLPE3ISLCA2O6smN/JX+bV+R2OMaYEGWJv4k5t2My5+el8KfP1nOgwjpwM8YEniX+Jugn\nI7uwp7yaF78scjsUY0wIssTfBPVp14qR3dvw/OwN7C2vcjscY0yIscTfRE0c2ZmDVTVMnr3B7VCM\nMSHGEn8T1TU9gct6Z/Dil0XsPFDpdjjGmBBiib8J+/GIPKpq6/jTZzZYizEmcCzxN2EdU1twdf9M\n/v71JhusxRgTMJb4m7j/uTgPVeWpTwrdDsUYEyIs8TdxWUlxXDcwm9cXbGFzWbnb4RhjQoAl/mbg\nzuG5RHiFSR+vdTsUY0wIsMTfDKQlxPC9wTm8vbiEwlIbotEYc3Ys8TcT/31BJ2IjvTz+kQ3Mbow5\nO5b4m4nW8VHcOrQD73+zjdXb97sdjjGmGfMr8YvIaBFZIyKFIvLACbbfICLfiMgyEflKRPrU21bk\nW79ERMKvk/0A+sHQjrSMjmDSLLvqN8acuQYTv4h4cUbVGgN0B64Tke7H7bYRuEBVe+EMtD75uO3D\nVbWvv4MEmBNLjIvk++d3YMaK7Swv2ed2OMaYZsqfK/6BQKGqblDVKuA14PL6O6jqV6q6x7c4D2dQ\ndRMEtw7tQEJMBJM+shY+xpgz40/izwS21Fsu9q07me8D0+stK/CRiCwUkfGnH6KpLyEmkvHDOvLR\nqlKWbtnrdjjGmGYooA93RWQ4TuK/v97qoaraF6eq6A4RGXaSY8eLSIGIFOzcuTOQYYWcW4Z0ICku\nksftqt8Ycwb8SfwlQLt6y1m+dd8iIr2BKcDlqlp2ZL2qlvjmpcBUnKqjf6Oqk1U1X1XzU1NT/T+D\nMNQiOoL/vqATn63ZycJNexo+wBhj6vEn8S8A8kSkg4hEAdcC0+rvICLZwFvATaq6tt76eBFpeeQz\nMBJYHqjgw9nNg9uTHB/F47Psqt8Yc3oaTPyqWgPcCcwEVgGvq+oKEZkgIhN8u/0SSAaePa7ZZhvg\nCxFZCswH3lfVGQE/izAUFxXBbRd24ovCXXy9oazhA4wxxkdU1e0Y/k1+fr4WFFiT/4Ycrqpl2COf\n0ik1ntfGD3Y7HGOMi0Rkob9N5u3N3WYsNsrL7Rd2Yt6G3XxVuMvtcIwxzYQl/mbuuoHZpCfE8Ois\ntTTFuzdjTNNjib+Zi4n0csdFuSzctIfZ6+yq3xjTMEv8IeA7+VlktorlsQ/X2FW/MaZBlvhDQHSE\nl7suymVp8T4+WV3qdjjGmCbOEn+IuHpAFtmt43jM6vqNMQ2wxB8iIr0efnRxHiu27mfmiu1uh2OM\nacIs8YeQy/tm0DElnsdnraOuzq76jTEnZok/hER4PfxoRB5rdhzg/WXb3A7HGNNEWeIPMZf1zqBz\nmxZM+mgttXbVb4w5AUv8IcbjEe4e0Zn1Ow8xbem/daJqjDGW+EPRqB7pdGubwBMfraOmts7tcIwx\nTYwl/hDk8QgTL+lMUVk5by2yq35jzLdZ4g9RI7ql0ScrkSc+XkdlTa3b4RhjmhBL/CFKRJg4sgsl\new/zzwVbGj7AGBM2LPGHsGF5KQzs0JqnPinkcJVd9RtjHH4lfhEZLSJrRKRQRB44wfYbROQbEVkm\nIl+JSB9/jzXBIyLcO6oLOw9U8tLcIrfDMcY0EQ0mfhHxAs8AY4DuwHUi0v243TYCF6hqL+AhYPJp\nHGuC6Jyc1gzvksqfPlvP/opqt8MxxjQB/lzxDwQKVXWDqlYBrwGX199BVb9S1T2+xXlAlr/HmuC7\nZ2QX9h2uZsrsDW6HYoxpAvxJ/JlA/aeDxb51J/N9YPoZHmuCoGdmIuN6tWXKFxvZdbDS7XCMMS4L\n6MNdERmOk/jvP4Njx4tIgYgU7Ny5M5BhGeDuSzpTUV3Lnz5b73YoxhiX+ZP4S4B29ZazfOu+RUR6\nA1OAy1W17HSOBVDVyaqar6r5qamp/sRuTkNuWguu7p/F3+ZtYtu+w26HY4xxkT+JfwGQJyIdRCQK\nuBaYVn8HEckG3gJuUtW1p3OsaTw/GpGHqvLkx4Vuh2KMcVGDiV9Va4A7gZnAKuB1VV0hIhNEZIJv\nt18CycCzIrJERApOdWwQzsP4ISspjhsGtef1gi0U7TrkdjjGGJdIUxymLz8/XwsKCtwOIySVHqhg\n2B8/ZVSPdJ64tp/b4RhjAkREFqpqvj/72pu7YSatZQz/NaQD05ZuZcXWfW6HY4xxgSX+MDThgk4k\nxkby8PTVbodijHGBJf4wlBgbyZ3Dc5mzbhdz1lnTWWPCjSX+MHXT4PZktorl4emrbWB2Y8KMJf4w\nFR3h5d5RXVixdT/Tlm51OxxjTCOyxB/G/qNPBj0yEnhk5hobrMWYMGKJP4x5PMJPx3SjZO9h/jZ3\nk9vhGGMaiSX+MDc0L4Xz81J46pNC9pVbt83GhANL/IYHxnRlf0U1z35uXTkYEw4s8Rt6ZCRyZd9M\nXvyyiJK91oGbMaHOEr8BYOLIzgA89uHaBvY0xjR3lvgN4HTgdst5Oby1uJiVW/e7HY4xJogs8Zuj\n7rgwl1axkTz43gqaYud9xpjAsMRvjkqMi2TiyC7M27CbmSu2ux2OMSZIItwOIKAm9YaaSvB4Qbwg\ncuyzNxIiYnxT9LF5ZBxEt4DolhDlmx+ZYpMgtjXEtXY+eyPdPsOgu+6cdrwybxO/fX8VF3ZJIybS\n63ZIxpgAC63E33kU1FRAXR1oHWgt1NU689pqp1CoqYCqQ1Be5ixXl0PlAag6CHU1p/7+6ASnAGiR\nBvFpzrxFm2PzhLaQkAXxqeBpnjdTEV4Pv7y0O9dP+Zq/fLGRO4bnuh2SMSbA/Er8IjIaeALwAlNU\n9eHjtncFXgT6Az9T1f+rt60IOADUAjX+DhRwRsY+cubHqjqFQuVBqNzvTIf3QPluZ3708244WAp7\nimDL104BwnH14Z5ISMiAxCxIyIRW7SApxzd1cLZ5mu6V9Hm5KYzq0YZnPi3k6v5ZpCfGuB2SMSaA\nGhyBS0S8wFrgEqAYZxzd61R1Zb190oD2wBXAnhMk/nxV3eVvUM1qBK7aaji0Cw5uh/3bYH8J7NsC\n+0p8n31zrdcXjicSWmU7BUFKHiTnOlNKnlNQiLh2OkdsLitnxGOfc2nvtjz23b5uh2OMacDpjMDl\nzxX/QKBQVTf4vvw14HLgaOJX1VKgVETGnUG8zZs30lfF0xYyTjKUYW017Ct27hLqT7s3wOZ5UF1v\n/NvIOKcQSO0Kad2OTYnZjVp9lJ0cxw/O78Czn63nxsHt6Z+d1Gi/bYwJLn8Sfyawpd5yMTDoNH5D\ngY9EpBb4s6pOPo1jQ4M3Elp3cKbjqcKBbbBrHZStg7L1sGstbPoKlr1+bL/IeEjtAm16QHpvSO/l\nfI5JCFrYtw/P5Y2Fxfzm3ZVMve08PB7370SMMWevMR7uDlXVEl910CwRWa2qs4/fSUTGA+MBsrOz\nGyGsJkLEqfNPyICOF3x7W8U+2LkGSldC6WooXQGr34fFfzu2T1KOUwi07ePccbTtB/HJAQmtRXQE\n94/uyj3/WsrUxSVcPSArIN9rjHGXP4m/BGhXbznLt84vqlrim5eKyFScqqN/S/y+O4HJ4NTx+/v9\nIS0mEdoNdKYjjtwhbF8G27+B7cud+ap3j+3TKtspBDL6QUZ/yOzvNE89A1f2y+TleZv4w4zVjOqZ\nTovo0GoIZkw48udf8QIgT0Q64CT8a4Hr/flyEYkHPKp6wPd5JPDgmQZr+PYdQudRx9ZX7INtS2Hr\n4mPTyneOHARp3SErH7LOcaaUzn49M/B4hF9d1p2rnv2Kpz5ex0/HdgvOeRljGk2DiV9Va0TkTmAm\nTnPOF1R1hYhM8G1/TkTSgQIgAagTkR8D3YEUYKo4rVQigFdVdUZwTiXMxSRCh2HOdET5bti6CIoL\noHiBUxAsesnZFp0I7c6B7HOh3bmQOQCi4k741f2zk/hOfhZ/+WIjV/bPpGt68J4rGGOCr8HmnG5o\nVs05m5O6Oti93ikItnztTKW+xlmeCOc5Qbtzof15zhTX+uihew5VcfFjn5OTHMcbE+xBrzFNzek0\n57TEH+7Kdzt3A5vnOVPJQqitBMRpNdR+COQMgfZDeGN1BT/511L+35W9uH5QGD2AN6YZCHQ7fhPK\n4lo7zwqOPC+oqXSSf9GXsOkLpwXR/D8DcHVad+KT85g+fTEjO44nJTXNxcCNMWfKrvjNqdVUwbYl\nUPQFFM2hbtNXeGoqqMODJ6OP77nCBU7VUGSs29EaE7asqscET00l/3z7bbYunsmtGZtJLFsCddXg\njXYeFHe8EDoNh/Q+zbajOmOaI0v8JqgqqmsZ88QcVJUZtw8gZut82PAprP/UeckMnO6sO14AuSOg\n08VOlxbGmKCxOn4TVDGRXn57RU9umPI1z365jYkjR0DeCGfjgR2w8XOnEFj/CayY6qxP6w65FzuF\nQPvznLEQjDGusCt+c8bu/ucS3vtmK9N/dD65aSd4M1gVdqyA9R9D4ceweS7UVjkd0XUYBnmXQO4l\nkNS+8YM3JsRYVY9pFLsOVnLxo5/TuU0LXhs/GG9DbfurDjkPidfNgsJZTg+lACldnEIgbyRkD4aI\nqKDHbkyoscRvGs2bC4u5519L+dnYbvxwWEf/D1SFskKnEFj3IWz60rkbiGrpPBzuPMopCFpYk1Fj\n/GF1/KbRXNU/kxkrtvPIh2u4sEsqeW387AxOxBl4JiUPBt/u3A1s+BzWzYS1H8Kqac5+Gf197xmM\ndt4sbgKD1BjT3NkVvzlrOw9UMmrSbLKSYnnztvOI9J5lM05Vp/fRdTNh7UyniwkUWvo6pusyxnlG\nYO8NGHOUVfWYRvfBsm3c/soiJl7Smf+5OC+wX35wp1MdtHY6FH7ijFgWGee8M9BljHM3YFVCJsxZ\nVY9pdGN7teXyvhk8+fE6LuqaRs/MxMB9eYtU6HeDM9VUQtEcWDMD1s6ANR8A4nQ53WWsM6V2sSoh\nY07BrvhNwOwtr2Lk47NJioti2l1DiI7wBvcHVWHHclgz3RmZbNsSZ31SB+g6zikEss8FT5DjMKYJ\nsKoe45pPV5fyX39dwG0XduL+0V0b98f3lfjuAqY7L5HVVkFcslMV1HUcdBx+0jEHjGnuLPEbVz3w\n5je8XrCFf00YzID2rRs+IBgqD0DhR86dwNoPoXIfRMRCp4ucQqDz6ICNTWxMU3A6id+v5hciMlpE\n1ohIoYg8cILtXUVkrohUishPTudYE3p+Nq4bbRNjuef1pRyoqHYniOiW0ONKuHoK3Lcebnob+t/k\nVAe9czv8Xy789VKY9yfYs8mdGI1xSYNX/CLiBdYClwDFOGPwXqeqK+vtkwa0B64A9qjq//l77InY\nFX/zN3/jbq57fh5jeqbz1HX9kKbysFXVSf6r33emIyOQpfeCrpc6U5se9nDYNDuBvuIfCBSq6gZV\nrQJeAy6vv4OqlqrqAuD4y7sGjzWhaWCH1twzsjPvfbONv3+92e1wjhGBjH5w0c/h9rlw1yIY+VuI\njIfPHobnhsCTfWHmz2DTV1BX63bExgScP805M4Et9ZaLgUF+fr/fx4rIeGA8QHa2DesXCiYM68SC\njbt56N2V9M1qRa+sADbxDJTkTnDeXc50sNTXQug9mD8Z5j4NcSnOuwLdLnMGnImMcTtiY85akxkp\nQ1Unq2q+quanpqa6HY4JAI9HePQ7fUluEcUdry5iv1v1/f5qkQYDvgc3/Avu2wDXvOiMKbDyHXj1\nO/BIJ3j9e7DsDajY53a0xpwxf674S4B29ZazfOv8cTbHmhDQOj6Kp6/vx3f/PI/7/vUNf7qxf9Op\n7z+V6JbQ8ypnqqmEjXOcO4E1H8DKt8ET6XQb0e1S532BluluR2yM3/y54l8A5IlIBxGJAq4Fpvn5\n/WdzrAkRA9q35v7RXZmxYjt//arI7XBOX0S0M9DMZZNg4mr4/iw49zbYsxHeuxse7QpTLoEvJkHZ\nerejNaZBfrXjF5GxwCTAC7ygqr8TkQkAqvqciKQDBUACUAccBLqr6v4THdvQ71mrntCjqvzw5YV8\nvraUf004j77tWrkd0tlThZ2rYdV7zt3AkTeHU7s67wp0vdR5kNwc7nBMs2cvcJkmaV95NeOemoMq\nvHvXUFrHh9iAK3u3OFVBq951WgRpLSRkOlVBXcdBzlDwRrodpQlRlvhNk7V0y17+889z6ZOVyN9/\nMCj4/fm4pXy306X06vecYSdrDkNMIuSNgq5jnUHoo/0cu8AYP1jiN03au0u3ctc/FnNVv0we/U6f\n5vGw92yhhbbJAAAQB0lEQVRUlcOGz5xCYO0MKC8Db5TTPPRIZ3It27gdpWnmrFtm06Rd1ieDol2H\neHTWWjqkxHNXoPvvb2qi4pyr/K5jnRfCtnzte3P4PXjvx84D4swBzvYu46xbaRN0dsVvXKGq3PP6\nUt5aXMJT1/Xjsj4ZbofU+FShdJVTCKz5ALYuctbX71a63SDw2vWZaZhV9ZhmobKmlhunfM3S4n38\n44fnMqB9ktshuWv/VufN4frdSscmOYPOdxkDnS6GmAS3ozRNlCV+02zsPlTFlc9+ycGKGt6+Ywjt\nWlt/+YCvW+mPnWcCa2fC4d3OS2M5Q48NN5nU3u0oTRNiid80K4WlB7nq2S9JT4zhjdvOIyHGmjx+\nS20NFM8/djdQts5Zn9bdGXy+8xhn6EkbaSysWeI3zc5Xhbu4+YX59GnXipdvHUh8tNVrn9SuQt+d\nwAzYPBfqaiC2tVMl1HmUM9hMbAi8IGdOiyV+0yxNX7aNO/+xmHNyknjxloHERtkVbIMO74X1HzvV\nQes+hMN7wBMB7c6FziOd9waslVBYsMRvmq13lpTw438uYWhuCs/fnE9MpCV/v9XVQvGCY4XAjuXO\n+lbZzt1A3kjIOd/GHQ5RlvhNs/Z6wRbue+MbLu6axp9uHEBURJPpPbx52VcM62Y5hcCGz6C6HLzR\nkDPEKQRyL3HGI7C7gZBgid80e3+ft4mfv72c0T3Sefr6fkR4LfmfleoK2PwVrPsICmfBrrXO+qQO\nTvcRuRc7dwPRLdyN05wxS/wmJLzwxUYefG8l/9Eng8e/2xevx65MA2ZPkXM3UPiRM9ZA9SGnuWj2\nuU4hkDsC2vS0u4FmxBK/CRnPfb6eh6ev5rI+Gfzff/YO3U7d3FRTCZvnOQ+JCz8+9mwgPg06DYeO\nw525DTbTpFniNyHlz5+v5/fTV3Nep2Seu2mAtfMPtgPbnQJgw6ew/lMo3+WsT+vuKwQugvaDISre\n3TjNtwQ88YvIaOAJnMFUpqjqw8dtF9/2sUA5cIuqLvJtKwIOALVAjT+BWeI3x3trUTH3vfENuWkt\neOnWgbRJsEHPG0VdHexY5hQAGz6FTXOhttKpFmo30OlhtOOFkNnfxhpwWUATv4h4gbXAJUAxznCK\n16nqynr7jAXuwkn8g4AnVHWQb1sRkK+qu/w9AUv85kRmr93JbX9fSKu4KF669Rxy06w/+0ZXVe68\nNLbxc9jwOWxbCihEtXRaC+Wc74xF3KYneOyBfGMKdOIfDPxaVUf5ln8KoKq/r7fPn4HPVPUfvuU1\nwIWqus0Svwmk5SX7uOXF+dTUKX/5Xj4D2rd2O6TwVr4bNs72FQSfwe4NzvrYJGg/xCkEcs6HtG72\noDjIAt0ffyawpd5yMc5VfUP7ZALbAAU+EpFa4M+qOtmfwIw5kZ6Zibx12xC+9+J8rn/+a564ti+j\ne7Z1O6zwFdcaelzhTAD7SqBojtNSaONsZ8wBgLhkpyDIGerM07rbHYGLGqNDlKGqWiIiacAsEVmt\nqrOP30lExgPjAbKzsxshLNNcZSfH8caEwXz/pQIm/H0R/31BR+4d2cXa+jcFiZnQ51pnAqfZ6MY5\nsOlLKPoSVk1z1scmQfZ50P4850Fxeh8bd6AR+fOXLgHa1VvO8q3zax9VPTIvFZGpwEDg3xK/705g\nMjhVPX7Gb8JUcotoXht/Lg++t5I/f76BxZv28tT1/eyhb1OTlONM/W9ylvdudgqATV848zXvO+sj\n450eRtufB9mDnc/Waiho/Knjj8B5uHsxTjJfAFyvqivq7TMOuJNjD3efVNWBIhIPeFT1gO/zLOBB\nVZ1xqt+0On5zOt5eXMJP31pGfLSXJ67tx5DcFLdDMv7av815WLx5rtNiaMdyQEG80La3MwLZkSkx\n0+1om7RgNOccC0zCac75gqr+TkQmAKjqc77mnE8Do3Gac/6XqhaISEdgqu9rIoBXVfV3Df2eJX5z\nutbtOMBtryxi/c6DTBzRmTuG5+KxN32bn8N7Yct82DLPmRcXQM1hZ1tiO6cJadY5zpTeCyKi3Y23\nCbEXuExYOlRZw8+mLuPtJVs5Py+Fh6/uTWarWLfDMmejthq2L3MGqN/ytVMY7PfVNHujoW0fX0Ew\nADLznZ5Iw7T1kCV+E7ZUlX/M38JD763EI3Df6K7cdG57u/oPJftKoKTAuRsoLoCti4/dFcQlQ+YA\nZ8ro77xYFh8eVX+W+E3Y27K7nP+duow563YxoH0Sf7i6l73wFapqq51nAyWLYOsiZ166CqclOZCY\nDRl9oG1fyOgLbftBfLKrIQeDJX5jcK7+py4u4cH3VlJeWcsdw3O57cJO1r9/OKg84LxVXLIQti6B\nbUuOvVwGTmHQtrdTVZTe2/ncsm2zriayxG9MPbsOVvLguyuZtnQrndu04H/HduOCzqlIM/5Hbs7A\n4b1OYbBtiVM9tO0b2L3+2Pa4FOeB8ZGpTQ9I6dxs+iCyxG/MCXy8age/fncFW3Yf5tyOrbl/dFf6\nZSe5HZZxU+UB2LHCKQS2L3XmO1dDbZWz3RMJqV0hvadTEKR1c946boJ3B5b4jTmJqpo6Xv16E099\nUkjZoSpG90jnJ6O6kJtmI08Zn9pqKCuE7cudZwc7ljuFw4Ftx/aJSXQKgCMFQWoXp4CIT3WtQLDE\nb0wDDlbWMGXOBp6fvYGKmjr+c0AWt1+YS3ayDURuTuJQGexc5Tw4Ll15bF6x79g+sUmQ0uVYQZDS\nGVJynXcQPMEdRMgSvzF+2nWwkqc/KeSVrzdRW6dc0r0Ntw7pwMAOre0ZgGmYqnMnsHM17Fzrm69x\n5od3H9vPGw3JuU4hkJwHKXnQupMz2H1cYHqYtcRvzGnavq+Cl+cW8er8zewtr6ZnZgLfH9qBcb0y\nrBWQOTOHdjmD2u9aB2XrnPmudU7HdVp7bL+YVk4B0LqTc4cw7CdnVF1kid+YM3S4qpa3Fhfzwhcb\nWb/zEGkto/nuOe24vG+mPQcwgVFT5ST/3euhbH29+QYQD/z4mzP6Wkv8xpylujpl9rqdvPhlEXPW\n7aROoVdmIlf0y+SyPm1Ja2m9gJogqK054+6pLfEbE0A79lfw7tKtvL2khOUl+/EIDMlN4dLebbmw\nS5p1BW2aBEv8xgRJYekB3lniFAJbdjv9w3Rvm8CFXVIZ3jWNfu1a2YAwxhWW+I0JMlVl1bYDfLa2\nlM/W7GThpj3U1ikJMREMzUshv31rBrRPontGApFWEJhGYInfmEa273A1Xxbu4rM1pXxZWEbJXudu\nICbSQ++sVgxon8SAbKcgaJsYY01FTcBZ4jfGZdv2HWbRpr0s3LSHhZt2s2LrfmrqnH9rCTERdE1P\noGvblnRNT6BLeks6pcbTKi7K5ahNcxaMEbhGA0/gjMA1RVUfPm67+LaPxRmB6xZVXeTPsSdiid+E\nmsNVtSzfuo/V2/azevsBVm8/wJrtBzhYWXN0n4SYCNonx5OdHEdOchztW8eTmRRLm4QY0hNjaBFt\ng5GbkzudxN/g/0ki4gWeAS4BioEFIjJNVVfW220MkOebBgF/Agb5eawxIS82yss5Oa05J+fYW5qq\nSvGew6zZfoCiskNsKiunqOwQy0v2MWP5dmrrvn1R1jI6gjaJMaQnxJDWMprW8VEkxUfRut6UFBdJ\nQkwkLWMiiYn0WJWSOSF/LiEGAoWqugFARF4DLgfqJ+/LgZfVuX2YJyKtRKQtkOPHscaEJRGhXes4\n2rX+9/6Bqmvr2Lr3MCV7D7NjfwXb91X65hVs31/Bxl2H2FNeRXlV7Qm+2RHpFVrGRNIyJoKWMRHE\nRUUQH+UlLiqC2Cgvcb7PMZEeYiK9REd4iI7wEhPpzKMiPER6hSivh8gID5FeZznS6yHCI0R4PHi9\n4vsseD2CxyN4xfns9X220c+aHn8Sfyawpd5yMc5VfUP7ZPp5rDHmOJFeD+2T42mfHH/K/Sqqa9l9\nqIrdh6rYU+7MD1TUsL+i2pkfduYHKqo5VFXLroNVlFeVU15VS3lVLYeraqmqrQv6+XgEPCLO5Dn2\nWXB6JxARPL65s06c9Ue2U3/5WEHiHOv7jBxd5ywf2efEBY+cdOHfFhsUqDur1nFRvD5hcEC+61Sa\nTKWhiIwHxgNkZ2e7HI0xzUNMpJeMVrFknMWg8rV1SlVNHRXVtVTW1FFZU0tFdR1VNXVU1dZR/a3J\n2bdOlZpapbZOqa6ro7bOWa5TZ12tKrW1zryuTlHf79SpU8Xl7AeKor51CtSps+zUch3ZVm8/nGXf\n1qOjKx6pFDvyzPLY8onPuf7q459znnZzlwC2j2kZ0zgp2Z9fKQHa1VvO8q3zZ59IP44FQFUnA5PB\nebjrR1zGmADweoTYKC+xUcHtNtg0Hf68WbIAyBORDiISBVwLTDtun2nAzeI4F9inqtv8PNYYY0wj\navCKX1VrROROYCZOk8wXVHWFiEzwbX8O+ACnKWchTnPO/zrVsUE5E2OMMX6xF7iMMSYEnE47futE\nxBhjwowlfmOMCTOW+I0xJsxY4jfGmDBjid8YY8JMk2zVIyI7gU1neHgKsCuA4TQHds6hL9zOF+yc\nT1d7VU31Z8cmmfjPhogU+NukKVTYOYe+cDtfsHMOJqvqMcaYMGOJ3xhjwkwoJv7JbgfgAjvn0Bdu\n5wt2zkETcnX8xhhjTi0Ur/iNMcacQsgkfhEZLSJrRKRQRB5wO57GICIviEipiCx3O5bGICLtRORT\nEVkpIitE5EduxxRsIhIjIvNFZKnvnH/jdkyNRUS8IrJYRN5zO5bGICJFIrJMRJaISFB7qQyJqh7f\noO5rqTeoO3BdqA/qLiLDgIM44x33dDueYPON49xWVReJSEtgIXBFKP93FmdMv3hVPSgikcAXwI9U\ndZ7LoQWdiEwE8oEEVb3U7XiCTUSKgHxVDfq7C6FyxX90QHhVrQKODOoe0lR1NrDb7Tgai6puU9VF\nvs8HgFU44zqHLHUc9C1G+qbmf7XWABHJAsYBU9yOJRSFSuI/2WDvJkSJSA7QD/ja3UiCz1flsQQo\nBWapasifMzAJuA8I/kjwTYcCH4nIQt8Y5EETKonfhBERaQG8CfxYVfe7HU+wqWqtqvbFGbN6oIiE\ndLWeiFwKlKrqQrdjaWRDff+dxwB3+KpygyJUEr8/A8KbEOCr534TeEVV33I7nsakqnuBT4HRbscS\nZEOA//DVeb8GXCQif3c3pOBT1RLfvBSYilOFHRShkvhtUPcw4HvQ+Rdglao+5nY8jUFEUkWkle9z\nLE4DhtXuRhVcqvpTVc1S1Rycf8ufqOqNLocVVCIS72uwgIjEAyOBoLXWC4nEr6o1wJFB3VcBr4fD\noO4i8g9gLtBFRIpF5PtuxxRkQ4CbcK4Al/imsW4HFWRtgU9F5BucC5xZqhoWzRvDTBvgCxFZCswH\n3lfVGcH6sZBozmmMMcZ/IXHFb4wxxn+W+I0xJsxY4jfGmDBjid8YY8KMJX5jjAkzlviNMSbMWOI3\nxpgwY4nfGGPCzP8H4qv4xNFXc34AAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<matplotlib.figure.Figure at 0x7f5004f15da0>"
|
|
]
|
|
},
|
|
"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": 37,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def distance(x, X): return sqrt(((x-X)**2).sum(1))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Let's try it out. (More on how this function works shortly)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 38,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([ 1.41421356, 0. , 3.60555128])"
|
|
]
|
|
},
|
|
"execution_count": 38,
|
|
"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": 39,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([ 0.13598247, 0.15957691, 0.05640321])"
|
|
]
|
|
},
|
|
"execution_count": 39,
|
|
"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": 40,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def meanshift_iter(X):\n",
|
|
" # Loop through every point\n",
|
|
" for i, x in enumerate(X):\n",
|
|
" # Find distance from point x to every other point in X\n",
|
|
" dist = distance(x, X)\n",
|
|
" # Use gaussian to turn into array of weights\n",
|
|
" weight = gaussian(dist, 2.5)\n",
|
|
" # Weighted sum (see next section for details)\n",
|
|
" X[i] = (weight[:,None]*X).sum(0) / weight.sum()\n",
|
|
" return X"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 41,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"X=meanshift_iter(data)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The results show that, as we hoped, all the points have moved closer to their \"true\" cluster centers (even although the algorithm doesn't know where the centers actually are)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 42,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4lOWd//H3l5yAKkQmHBIOBksiRY1So+thWxNFRXEl\nuu162Gt/mFrodj2gbfFQ3WbSWutPt+B66HZBi/y69VQP0RbWKjTRHkQbRRElhigRJBFkMAQkkAP3\n749nZpiBcJAcZ57P67pyzcxzPzPPc1/oZ+75zv3cY845REQk+Q3o6xMQEZHeocAXEfEJBb6IiE8o\n8EVEfEKBLyLiEwp8ERGfUOCLiPiEAl9ExCcU+CIiPpHa1ycQKysry+Xm5vb1aYiIJJQ33nhjs3Nu\n+MH261eBn5ubS3V1dV+fhohIQjGzjw5lP5V0RER8QoEvIuITCnwREZ9Q4IuI+IQCX0TEJxT4IiI+\nocAXEfEJBb6IiE8o8EXk0LU0w4pnvdtkOI7PKPBF5NDVLIPli7zbZDiOz3R5aQUzGwi8AmSEX+8p\n51yZmQ0DngBygXrgn5xzn3X1eCLShyaeE3+b6Mfxme4Y4e8CznbOnQicBEw1s9OAW4Blzrk8YFn4\nsYgkskFDYPIl3m0yHMdnuhz4zrM9/DAt/OeA6cCi8PZFQElXjyUiIoevW2r4ZpZiZm8Bm4CXnHOv\nASOdc43hXT4BRu7nubPMrNrMqj/99NPuOB0REelEtwS+c67DOXcSMAY41cyO36vd4Y36O3vufOdc\noXOucPjwgy7nLCLd4YvMgondV7NnElq3rofvnGsys0pgKrDRzLKdc41mlo03+heR/iAyCwYIjfs6\ngUBgn11CoZC3PbLv+re8hg1ve7eTL+mts5Vu0uURvpkNN7PM8P1BwLlADfA8MCO82wzgua4eS0S6\nycRz4LQZBJ96nYKCAmpra+Oaa1e+ScFX8gjO/EfIngRDRnlBv+FtGHey93yN9hNOd4zws4FFZpaC\n9wbypHPu92b2KvCkmV0NfAT8UzccS0S6w6AhBJ96nfI77wKguLiYyspK8vPzqV35JsXFRTRs2Ub5\nQ89A8yaCReFPANnHwYgJ3v3IyL9hFZw9WzNqEkCXA985txKY3Mn2EKBJtCL9UPC2W6NhD9DQ0EBx\nURHzZ1/KrLt+RUNTS7St/Mk/w/ZjCZZMht0dUP0EpA70RvkNq2DdG174q8TT7+lKWxGfCYVCLFgw\nP25bCSXsaGzholseZEdTBiV7zaJe8JePCDVthY01e0o6g4Z4I/vTZugCqQShwBfxk5ZmAuteofLZ\n35Az7EjAC/vZzGYe88gll3nMYzazo6GfM3QglTecSeCIdDhyZHz5RhdIJRQFvoifvLMYli8iP/Qq\nldedSs7QgVRRRT315JLLQhaSSy711FNFVTTs80cNhVHHwbQfKdwTmAJfxE8sfJs1nvwLr2b+dy+g\niSbKKY/brZxymmhi/pUnkj/yCHAdkJIKda9oVk4CU+CLJLvY6ZPHT4PCyyA1g9qBxzJr/ktkkkkZ\nZXFPKaOMTDKZ9ejb1HYM90b3G972vrDVCpYJS4EvkuxWeWUcVi32Hm+qo3bJwxSfcw4NW7ZTRFG0\njFNKabS8U0QRDVt3Ulz+W2obQt5zh2RD7ql91xfpkm690lZE+iEXc1uzjNB7r1J833IaPvscgAoq\nAKiiiiaauJEbKaIour1h81aKb1vEyruvJNDcCPWvw1GagpmINMIXSXYnTPOmTp4wDSaeQ2DS6cw8\nfUzcLhVUMHjoTn7/b2cwOLM1GvYRM//xfALHfw1Ovswb4esK24SkwBdJdrFTJ8Nz54N3/Iyyqy6K\n7pKTlUnlzecz7fjhVP7fb5OTvWdx27JvX0rwO5fB2xWQNtAb4evXqBKSSjoiftHS7IV0uAYfnPtL\nSL2eBc9VUrl0KfkjvgSV95Ff/G9Ujsin+KofMPOCUwl+tQPad+17gZUutko4CnwRv4hd+2bdGwAE\nFzzNdXeFV8Vc8Sw0bYDG98g//19Y+eyRBDJ2h0f2GfFLJ2gZhYSkwBfxi8iIPHuSdxse6UeXRo79\nHdlBQwgU/x/vU8GgoRrNJwnV8EX8IlLLb3zPG+HXvx4/R3/QEC/Ya5bt+UJWSyckFY3wRfwmdiQf\n80MoTL5k38eSVBT4In4TGbXDnvCPTLWMXFSlEk5SUuCL+Fkk/Fc8q5G9DyjwRSS+zCNJS4EvIvFl\nHklamqUjIuITCnwREZ9Q4IuI+IQCX0TEJxT4IiI+ocAXEfEJBb6IiE8o8EVEfMI3gR9qcfxiRSuh\nFnfwnUVEkpBvAv+JmjbuWN7KEzVtfX0qIiJ9osuBb2ZjzazSzN4zs3fNbHZ4+zAze8nM1oRvj+r6\n6R6+yyamcftp6Vw2Ma3T9r0/AegTgYgkm+4Y4bcD33fOTQJOA64xs0nALcAy51wesCz8uM8EBhn/\nNjmdwCCLC/NQKATs+wngP6oauGN5Kz//285DDn69SYhIf9blwHfONTrn3gzf3wasBkYD04Hweqss\nAkq6eqyuigTywndauWN5K8VX385xJxRw81OrCLXs5nuFaVw2MY3XVr7P3H8+hW2/v5Pn6joOuRSk\nspGI9GfdWsM3s1xgMvAaMNI51xhu+gQY2Z3HOlSxo+5IILe0Q9rSO3nnsTvZ2NjAf3z7PO5/cTUt\nbY7Lf7WSrxUV0/pZI9uX/Ix1T9/JOeNS9lsKinWwspGISF/qtuWRzewI4GngBudcs5lF25xzzsw6\nrXOY2SxgFsC4ceO663SiIiEfatnNqs27ueq4VB6//8ese+Zn0X12b22k5cFp/M/2+6l76Dp2b22M\ntm1f8jOOOj6NJ3Juiwb5EzVtXDYxjcAgiztWpGwkItIfdUvgm1kaXtj/xjn3THjzRjPLds41mlk2\nsKmz5zrn5gPzAQoLC7u9+B0J6ZfXt/OnDbtZ2/ApNf+7MG6fEkqo2lxF7c+/QSaZFFFCBRXR9vkP\nLSAz61tANgB3LG8FULiLSELpjlk6BjwMrHbOzY1peh6YEb4/A3iuq8c6HJFR90+/NpBzxqUwZdJI\nAjcsJjXTC+8SSpjNbOYxj1xymcc8ZjObkvBXDgOGZjPkusWcMSGL83JT2dHu+N7JaSrbiEjC6Y4a\n/pnAvwBnm9lb4b8LgbuAc81sDTAl/LjHHGyGzISjBnDv2QMZNtAoPvlYjpq9mAFDs6miinrqySWX\nhSwkl1zqqaeKKgYMzfbeHEbm0e6MF+vbmVvdxuA0r5SjGTkikki6Y5bOn51z5pwrcM6dFP5b4pwL\nOefOcc7lOeemOOe2dMcJ709nM2Tipl+2OG74407mvtHGKSNT+LuCPIZeeR9NNFFOedxrlVNOE00M\nvfI+UkfmAfDOpt3UfdYRHd1rRo6IJJqk+U3bSIklttQSCeWIZes6+NroAWCQ0/IhWx+9nkwyKaMs\n7rXKKONGbmTro9eTGh7h73Lw+Psd3H5aCoFB1unxRET6s6RZWiH2wqqIyDTJ2Nr7Kdkp3L34PR65\nbiq7tzZSRFG0jFNKabS8U0QRu7c2Erp3GkO2rgFg5GDY0e6o+2z3fmfqiIj0V0kzwu9M5E3gFyta\nmVvdxvdOTiO0eTNb75tG62fe1MvIbJwqqmiiiRu5kSKKott3b23k4/+YxhE3v8rYkVnMrW7jb40d\n/GnDbkAzdUQkcSTNCP9AIiN9DBbVDyX99Kvi2iuooHnoII767m9pHjoobkomwBmXfIsfTcnmlFEp\nAByfNUAXWIlIwknqEX5EZKQfanHggMIyNp2Uxj13/hiA1Mxsjprt1epTb1hM6N5p0Yuvzv32bTx2\n30+ia/AEBg1QKUdEEpIvRvgRgUFG6QnpDE41br49SFlZGTk5OXzrgT+QOjKPE4cP4LrzvsIPHnqR\nUdk5lJWV8eKCO6LhHvmy9omatujMH03NFJFE4YsRfmQdndjplADBYJDrrruOzwYcxY6/7qLsjAwm\nHDUAOJ6bilcSCAT2eY0dbY65b+yZiqmrbkUkUfgi8GND/rzcVF5t6OC8XK/rgUCAJ1a0smxdB6fn\ntDPhqPTo9s5e43uF+y6Qplq+iCQCXwR+7Jz5J2ra9gn3Q5lTH7tPbP1eI3sRSRTmXP+pPxcWFrrq\n6uoePUZseUdfvIpIMjCzN5xzhQfbL+m/tN37i9XOLtASEfGDpA98rXkjIuJJ+hq+1rwREfEkfeDr\nV6hERDxJX9IRERGPAl9ExCcU+CIiPqHAFxHxCQW+iIhPKPBFRHxCgS8i4hMKfBERn1Dgi4j4hAJf\nRMQnFPgiIj6hwBcR8QkFvoiITyjwRUR8QoEvIuIT3RL4ZvYrM9tkZqtitg0zs5fMbE349qjuOJaI\niBye7hrhPwJM3WvbLcAy51wesCz8WERE+ki3BL5z7hVgy16bpwOLwvcXASXdcSzpPc2unWdaN9Ls\n2vv6VESkG/RkDX+kc64xfP8TYGQPHkt6wNK2EI+0NrC0LdTXpyIi3aBXftPWOefMzHXWZmazgFkA\n48aN643TkUM0JS0Qdysiia0nR/gbzSwbIHy7qbOdnHPznXOFzrnC4cOH9+DpyBc1xFK5NH0kQyzp\nf+texBd6MvCfB2aE788AnuvBY/mW6uwicqi6a1rmY8CrwLFm9rGZXQ3cBZxrZmuAKeHH0s32V2cP\nheIfR94Y6jdv7M3TE5F+pLtm6VzhnMt2zqU558Y45x52zoWcc+c45/Kcc1Occ3vP4pFuMCUtwFXp\nOXF19mAwSEFBAbW1tdFtS9tCPPDuq5x84kkEg8E+OFMR6Wu60jbJBINBysvLaWhooKi4OBr6uWtD\nvH3Rd9nS8Anl5eVc9u8/YMPunSoHifiIvo1LcJGSzqqO7bTctYi7fvyTaFtjQwNfLy6i6P4gL1z7\nI7Y27innPHnHz9m0u5Ujby0F4NJ0zZoVSXYa4SeI/X05OyUtQGHKEP66aR3zH1oQ11ZCCbsaWnji\nH7+DNe6iZK9r395a+CTTmtP5u9ShGumL+IACP0FERvL37vwoWorZsHsnS9tCXJ0xmlnZk6hY9iKD\nsr2prSWUMJvZzGMeueQyj3nMZnY09Idmj2Ty4l9SPcTxYqsusBLxA5V0EsSUtACrOrZT3dEMu6C6\no5kV7c28vXs7O10HV2bk8MwxcNri+bw6bRZVjVVMZzq55LKQhQDUU08VVQzMHs6Zi/+blLxxbKSN\nd9q3cXnaKHa6Dppdu+bdiyQpjfATxBBL5YaBR3NVeg5XZ4zmqvQcjkkZDEB1ezMzP1/Fx+07+XL+\nBAruu50mmiinPO41yimniSYK7rudlLw9VzVvpZ3VHZ/zeNtGft/6aa/2S0R6jwI/gUSufB09YCCX\npo/kvLQAYyyDOtfCRtfG0t1bOObDLbx3/U/JJJMyyuKeX0YZmWSy8vo72L7mIwAygGGk8fbubeG9\nOl0BQ0SSgAI/gb3WvpWP3S6OZTAjLY2TP2zm3qmXs71xE0UUkUsu9dRTSin11JNLLkUUsbPxU16d\nNovtaz6iFXifHWRgTBkwjIvSR/R1t0Skh5hz/WdEV1hY6Kqrq/v6NBLGht07eXjXBq7OGM3Azz6n\noKCAhoaGaHsJJVRRRRNNZJJJEUVUUBFtH5g9nLNefYLBgUzagWzSOSvtKMC4KH24avkiCcLM3nDO\nFR5sP43wE9hr7Vup7mjmtfatBAIBZs6cGddeQQU7s9M49bf/ya7stLiwBygqvZKjs0bQDgwlha+m\nDuHxto083vYJ97Ss1YVZIklGQ7gEtvfyxcFgkGbXzrwf/xSAI7JHcMGSh9k5YTRfCs/e2dnofSn7\nT7d/n5ab/5lG10phyhBuGHg0APUdLbzrPuft3dt5eNcGb1YQujBLJBko8BNMs2tnaVuIKWmB6Je4\nseaW30EGA5j/0AJe+uMyJuTn8/SujdQe+yXO/cNT3D/1CmbNnMn3ym7n962b2Lt8c0LqEbzb9jkn\nDjiCqzNGc3z7EVoPXyRJqIafYJ5p3cgjrQ1clZ5zwFF3KBQiEPCC+tFdjTze9gmXp43i/O3p0e2d\n2fsNRUT6v0Ot4ev/6ARzqL9CFR/qLnp7oLAHOv3UICLJQYGfYA4nkC9KH8FAS1FpRsTnFPg+oFG7\niICmZYqI+IYCX0TEJxT4IiI+ocAXEfEJBb6IiE8o8EVEfEKBL5IgQqF9f4Jyx2ZYEgyxY3MfnJAk\nHAW+SD+1YzP85R7vNhgMUlBQQG1tbdw+z91dy+XlBVx7RTBuf5HO6MIrkX7q9Qfg5XJYsDjIope9\nn6ssOquYqpcryc/P5+3Xarnuv4rZRgMLl5bT8Q045uUg9VVwySIYnNW35y/9j0b4Iv2UAyrZE/YA\njZ80UFxczOLFiznn3GJC2/f84M3/e7mctycEqVsCKxb2/vlK/6fAF+lHYssyE64M8e6RC+LaSyhh\nR8MOLrroIjq27aCEkrj2VzYu4KQ5ISaX9uZZS6JQ4Iv0IysWwtKbvNux+QHu+9dKjiQH8MJ+NrOZ\nxzxyyWUe85jN7GjoZw3J4YptlTSvCPDaA6rly75UwxfpRyIj88mlXmCPGpzPNSMreXBjMVVUMZ3p\n5JLLQryaTT31VFHFkeRwwzGVHLE1nw+XwodLoW0HnHd3H3ZG+p0eH+Gb2VQze9/M6szslp4+nkh/\ndaBZNJE2gDPneLdPXQGvlMOQHfn8A/NpoolyyuOeV045TTTxD8yn/a18mtZCykCvbeOKHuyMJKQe\nDXwzSwEeBC4AJgFXmNmknjymSH8VW645WNuKhbB2qXe/YVstv2MWmWRSRlnc88ooI5NMfscstg3x\npmx27ITBw6HoJz3ZG0lEPT3CPxWoc8596JxrBR4HpvfwMUX6pcmlMOVuOv1Cde+2yaUwfgpsppZF\neFMviygil1zqqaeUUuqpJ5dciihiGw080lFM+im1HH0W7PgU1v2pd/sn/V9P1/BHA+tjHn8M/F3s\nDmY2C5gFMG7cuB4+HZG+MzhrT7nmUPYtfjDE9wuL2bbNm3pZQQUAVVTRRBM3ciNFFEW3hz5vYG5N\nMW9Wr2TdcwHN1JF99PksHefcfOdcoXOucPjw4X19OiJ9orNyz9j8ALP+dWbcfhVU0MFgruT3dDA4\nGvYRF351JmPzA5w5Rxdeyb56OvA3AGNjHo8JbxORsB2bofVzOKts33LPhcODnBVTtz+SHGZQyddP\nmcYPTqhkaEpOtG3GWWX811PBXjprSUQ9Hfh/A/LMbLyZpQOXA8/38DFFEsqKhd5snLQv7Tsqn1wK\nP5wT5B/GlJGZlsNd36jklCn5HFME7p18/ut7leTk5FBWVsYjVUGN6uWAerSG75xrN7NrgT8AKcCv\nnHPv9uQxRRJN7Be1exuc5c2l/9LwIF+56TqG7ghQtxRGTYYJF8K5V+ez8uaVBAKB3j1pSUg9fuGV\nc24JsKSnjyOSqA7ly1zvzSDAsRdDbhG0fQ51S7z7Z85R2Muh0ZW2Igkg9k0ha45X90/7UuefCkT2\nR4EvkoC+yBRPkYg+n5YpIiK9Q4EvIuITCnwREZ9Q4IuI+IQCX0TEJxT4IiI+ocAXEfEJBb6IiE8o\n8EVEfEKBLyLiEwp8ERGfUOCLiPiEAl9ExCcU+CIiPqHAFxHxCQW+iIhPKPBFRHxCgS8i4hMKfBER\nn1Dgi4j4hAJfRMQnFPgiIj6hwBcR8QkFvoiITyjwRUR8QoEvIuITXQp8M/ummb1rZrvNrHCvtlvN\nrM7M3jez87t2miIi0lWpXXz+KuBS4L9jN5rZJOBy4DggB1hqZvnOuY4uHk9ERA5Tl0b4zrnVzrn3\nO2maDjzunNvlnFsL1AGnduVYIiLSNT1Vwx8NrI95/HF4m4iI9JGDlnTMbCkwqpOm25xzz3X1BMxs\nFjALYNy4cV19ORER2Y+DBr5zbsphvO4GYGzM4zHhbZ29/nxgPkBhYaE7jGOJiMgh6KmSzvPA5WaW\nYWbjgTzg9R46loiIHIKuTsu8xMw+Bk4HFpvZHwCcc+8CTwLvAS8A12iGjohI3+rStEzn3LPAs/tp\n+ynw0668voiIdB9daSsi4hMKfBERn1Dgi4j4hAJfRMQnFPgiIl3U6lr4oHUFra6lr0/lgBT4IiJd\ntL6thtWty1nfVtPXp3JAXV0tU0TE98amTQRgZGouH7SuYGzaRNJtUB+f1b40whcROQyxZZxtW3bw\n5fTJbGyvjxvph0KhPj7LeAp8EZHDECnj/OBHN1BQUMCq91eyY/c2BtsQjhgwjKWrKigoOIFgMBj3\nvL6s9yvwRUS+oFbXQrtr54W73uT+O+bT0NDA2cXF/HX1H9nhmnlh1W+57LyraGhopLy8nB+W3Rx9\nbl/W+xX4IiIH0NmIfH1bDeXlQR786UPRbZ82buG2af9J9Qur+f6Fd7KlcWu07Wc/vpt/L7sN8Or9\nX0k/LVr3700KfBGRA4gdkUfCP3XrEP646G9x+5VQQkcjlH/zfjoaHSWUxLX/csEv+ejTD1jb+g7t\nrr03uxClwBcR2Y9I6SYv7WTGpk1kbes7rG5dzs4hIV6p/DNZ2cMAL+xnM5t5zCOXXOYxj9nMjob+\nsOyh3LH4OjYeUcOatjdY01YdV9Lprbq+Al9EZD/Wt9Wwpq0aMNa2riLU0RBuMXLzxnLP4psZlj2U\nKqqop55cclnIQnLJpZ56qqhiWPZQHnnhQYomTWNSxhnkpZ1MXlphXEmnt+r6CnwRkf2I1NvBsaat\nmi27GxmRMo7RaRN4s+UlsvKO5Nr7rqSJJsopj3tuOeU00cS1913JoGM6aHUtrGtdDcD49OPj5un3\nVl1fF16JiOxHug3iy+mTw6UWAxzj009gfVsNm3dvYMOajTxw/aNkkkkZZXHPLaOMG7mRB65/lJzF\nI9iZ93lc+7EZp+5znJ6mEb6IyEGk2yCOzTiF8eknsLb1HXa5FprrOvjhtHvZ0riVIoqiZZxSSqPl\nnSKK2NK4lR9Ou5cNazZipIRf0fqkH+Zc//nd8MLCQlddXd3XpyEiso9W18JbO//Ipo51NIe2M/v0\nu9jcuCXaXkIJVVTRRBOZZFJEERVURNuHZQ/l6dcWMnZkbrcvvWBmbzjnCg+2n0b4IiKd2HvmzPq2\nGjZ1rGMwQzhp1JnMnPntuP0rqGBAtuNHv/0uA7JdXNgDlJRO5YwxU/ly+uQ+W2dHgd+dNm+Ge+7x\nbkUkoe09c2Zs2kRGpIxjB81k2CCCwSDX3T4ruv+w7KHctfgHnDL1BO5cfAPDsodG26649UK+/+/X\n9/mCagr87rRwIdx0k3crIglt75kz6TaIkwaeHd22vq2G824+ietun8WInOHcufgGRueNAGB03kh+\nvviHZGUP46pbv8mVP7yor8r2cVTD/yI2b/bCvLQUsrIAbzW8QCCwT3vIbM92EUl4ra6F9W010fp7\n7OM3G//M5iF1AKSRwdjUiUzImMy2LTvYfMQHrGl7g7y0k+Nm5nQn1fB7wl4j+GAwSEFBAbW1tV57\nVhbMmUPtli0UFBTss0qeiCSuvUs8kamU6TaIr2b/PceknkjWgDGcOfgSJg08g3QbRCAQYHz6CeSl\neVnc17+IpXn4X0RpqXd78cUEi4oof/llAIqLi6msrCR/2DBqb7mF4l//mobWVsrLy+Gllwieey5c\ne230U4GIJJ5IaWd/F0dlDBjEVzOm7FOnT7dBpFoqq1uXk2ppvTLffn9U0jkMwTPPpPyvf43bljNq\nFPOLi5n12GM07LV/GRCcOBEqKiAQ2KcsJCKJ7YPWFaxuXc5X0k/rNND3Lgd1N5V0uiJ2ts1eM29C\ntbUseO21uN1LKGHHJzu56LHH2EHmPqvkLTAjVFMDF1wAd9/tlYWuuEKzeUSSxMGWRogt//QllXQ6\nE6nVR0Tuz5lD4LnnqOzooDglhYaOjugqedOZTjnllFFGLrmANy83JyWFyo4OAgBr18LTT3uvtXSp\nd5w5c3qzZyLSA3praYSuUuB3JlKrj9zG3r/4YvKrqqhcv57id96hiiqmMz26Sh4QXSUvB6gcM4b8\no4+Gjz7y/j78EMaPh298I/71RUR6mEo6B5OV5QXzwoXw/vvwne/AkiXk19czHw64St58IP+TT+CV\nVyA3F04/3dth7VoYPHhPDV8XbIlIL+hS4JvZPWZWY2YrzexZM8uMabvVzOrM7H0zO7/rp9qLIiWd\nGTP2zK2/6SZvpk14Zk7ttm3Mgv2ukpdJJrOA2owMb+PLL8PGjQc+ni7YEpEe1NUR/kvA8c65AqAW\nuBXAzCYBlwPHAVOBX5hZyn5fpb8pLYULL4QlS/bMqCkrg8+95U1rgWKgAQ64Sl4DUNzcTO2YMd7r\nfvghTJjg1e2vvTb+eHffrRKPiPSoLtXwnXMvxjxcDnwjfH868Lhzbhew1szqgFOBV7tyvF6TlQVz\n53r3TzjBG+kfdxy8+ioh9oQ9EF0gKbJK3o3cGLdKXgNQvGkTK6+9lkBNjfdl7bvv7ns8fXkrIj2s\nO7+0/RbwRPj+aLw3gIiPw9v2YWazgFkA48aN68bTOQybN8MDD+x5vGSJNyqvqYHWVigrI/CXvzBz\n6dK4qn0FFeQA/wPMsq1UuPhV8ma2thJ44QX49a8hPX3PJ4dIyHeyZIOISHc7aEnHzJaa2apO/qbH\n7HMb0A785ouegHNuvnOu0DlXOHz48C/69O61cCGUl3t/4JVZ5s3zyjDHHuuVYR57jODdd1N2zTXR\np+WMGEHlKacwbfx4Kp0jJ+Yly8aNIwhQV+cF+ty5+5ZvVMMXkV5w0BG+c27KgdrN7CrgIuAct+ey\n3Q3A2JjdxoS39W+lpdE6fXQphHvu8cK6rg6OPjoa1MHwSHzBkCFU/ulP5A8bBg88QP7mzVT+7ncU\nr1vHzAkTCNbVwde/Dg0N3ieF55/ft3zT2TRQEZFu1qWlFcxsKjAXOMs592nM9uOAR/Hq9jnAMiDP\nOddxoNfrl0srxJZ5rr12z2h8yhQ480xCV15JID/fe2O46abol72hKVMIPPCAF/CRIF+4EC6+2NsW\nuVUZR0S66FCXVuhqDf8BIAN4ycwAljvn/tU5966ZPQm8h1fqueZgYd9vZWVB7KqXpaVQVeXV4c87\nzwv7yHY1cMbOAAAD8ElEQVTwgryoiEAkyGNH83Pm7HljiLxGZLuISA/T4mmHoytfskaeqxG+iHST\nQx3hK/BFRBKcVssUEZE4CnwREZ9Q4IuI+IQCX0TEJxT4IiI+ocAXEfEJBb6IiE8o8EVEfKJfXXhl\nZp8CH/XBobMAP/2+oPqb3PzWX/Bfn/fu79HOuYMuN9yvAr+vmFn1oVyllizU3+Tmt/6C//p8uP1V\nSUdExCcU+CIiPqHA98zv6xPoZepvcvNbf8F/fT6s/qqGLyLiExrhi4j4hG8D38zuMbMaM1tpZs+a\nWWZM261mVmdm75vZ+X15nt3FzL5pZu+a2W4zK9yrLen6G2FmU8P9qjOzW/r6fLqbmf3KzDaZ2aqY\nbcPM7CUzWxO+Paovz7E7mdlYM6s0s/fC/z3PDm9Pyj6b2UAze93M3g73tzy8/bD669vAB14CjnfO\nFQC1wK0AZjYJuBw4DpgK/MLMUvrsLLvPKuBS4JXYjUncX8L9eBC4AJgEXBHubzJ5BO/fLdYtwDLn\nXB7e70kn0xtdO/B959wk4DTgmvC/abL2eRdwtnPuROAkYKqZncZh9te3ge+ce9E51x5+uBwYE74/\nHXjcObfLObcWqMP7MfaE5pxb7Zx7v5OmpOxv2KlAnXPuQ+dcK/A4Xn+ThnPuFWDLXpunA4vC9xcB\nJb16Uj3IOdfonHszfH8bsBoYTZL22Xm2hx+mhf8ch9lf3wb+Xr4F/G/4/mhgfUzbx+FtySqZ+5vM\nfTuQkc65xvD9T4CRfXkyPcXMcoHJwGskcZ/NLMXM3gI2AS855w67v6k9dI79gpktBUZ10nSbc+65\n8D634X1M/E1vnltPOJT+ir8455yZJd1UPDM7AngauME512xm0bZk67NzrgM4Kfw947Nmdvxe7Yfc\n36QOfOfclAO1m9lVwEXAOW7P/NQNwNiY3caEt/V7B+vvfiRsfw9BMvftQDaaWbZzrtHMsvFGhknD\nzNLwwv43zrlnwpuTus8AzrkmM6vE+87msPrr25KOmU0FbgIuds7tiGl6HrjczDLMbDyQB7zeF+fY\nS5K5v38D8sxsvJml4305/Xwfn1NveB6YEb4/A0iaT3fmDeUfBlY75+bGNCVln81seGQGoZkNAs4F\najjc/jrnfPmH9+XkeuCt8N8vY9puAz4A3gcu6Otz7ab+XoJXw94FbAT+kMz9jenbhXizsD7AK231\n+Tl1c/8eAxqBtvC/79VAAG/mxhpgKTCsr8+zG/v793hfWq6M+X/3wmTtM1AArAj3dxXwo/D2w+qv\nrrQVEfEJ35Z0RET8RoEvIuITCnwREZ9Q4IuI+IQCX0TEJxT4IiI+ocAXEfEJBb6IiE/8f3aYeuD0\ne86BAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<matplotlib.figure.Figure at 0x7f5004dd4dd8>"
|
|
]
|
|
},
|
|
"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": 43,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def meanshift(data):\n",
|
|
" X = np.copy(data)\n",
|
|
" # Loop through a few epochs\n",
|
|
" # A full implementation would automatically stop when clusters are stable\n",
|
|
" for it in range(5): X = meanshift_iter(X)\n",
|
|
" return X"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 44,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 1.28 s, sys: 0 ns, total: 1.28 s\n",
|
|
"Wall time: 1.28 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": 45,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSJJREFUeJzt3VGMY9d93/Hvv7LjZKqHsaiFIloyaKBiAMVgbGDgumgf\nhomNyNlFNAmQQnko1K3BQQ0XGAcBArkGSvLBgIEAnZekbknUGz0YMQwkHgu7aV1ZmIFRoLE9apyt\nZFuUkoxheWRpQ2OQCgMolXP6MCR3uLuzO7Mkh0Oe7wcYDHnOJe85XOmHO3+ee2+klJAk5eMfTXsA\nkqTTZfBLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMvOOaQ/gsPvvvz+VSqVpD0OS\nZsrzzz//tymlc8fd/kwFf6lUYnt7e9rDkKSZEhE/OMn2lnokKTMGv6Q76na7J2rX2WbwS7qtRqNB\npVKh0+kMtXc6HSqVCo1GYzoD010bOfgj4mcj4lsR8ZcR8WJENHvt90XEsxHxcu/3u0cfrqTT1Gg0\naDab7O7uUq1WB+Hf6XSoVqvs7u7SbDYN/xkzjiP+t4BfTin9EvAB4LGI+DDwFPBcSukR4Lnec0kz\noh/6ff3wv3LlyiD0+wz/2TJy8KcDb/aevrP3k4DHgad77U8DK6PuS9Lp6Ha7tNvtobYVVtjf3efC\nhQvs7+6zcsP/0u1225r/jBhLjT8i7omI7wBvAM+mlL4JPJBSeq23yY+BB8axL0mTVygU2NzcpFgs\nAgehv8Ya66xTosQ666yxNgj/YrHI5uYmhUJhmsPWMY0l+FNKP00pfQB4CPhQRLz/hv7EwV8BN4mI\n1YjYjojta9eujWM4ksagXC4Pwn+LLXbYoUSJS1yiRIkddthiaxD65XJ52kPWMY11VU9KaQ/YBB4D\nXo+IBwF6v9844jWtlNJSSmnp3Lljn3gm6RSUy2VarRZ77NGkOdTXpMkee7RaLUN/xoxjVc+5iFjs\nPf454KPA94FngCd7mz0JfHXUfUk6XZ1Oh9XVVRZZpE59qK9OnUUWWV1dvWmpp862cRzxPwhsRsRV\n4Nsc1PgvA58DPhoRLwMf6T2XNCMOL9lcZnlQ3rnIxUHZZ5nlm5Z66uyLg/L72bC0tJS8Vo80fd1u\nl0qlMrRkc4UVtthijz0WWWSZZTbYGPQXi0WuXr3qF7xTEBHPp5SWjru9Z+5KukmhUKBWqw21bbDB\nQnGBy5cvs1BcGAp9gFqtZujPCINf0i01Gg3q9et1/f7qnfPnzw8t9QSo1+uewDVDztRlmSWdLf0w\nb7fbQ0s2+0s9q9UqtVrN0J8x1vil3H3+0Bm4n9i45SbdbveWZZyj2nW6rPFLGrujwt3Qn00GvyRl\nxhq/lLsjyjuaXx7xS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqM\nwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmayCv5ut3uidkmaR9kEf6PRoFKp0Ol0hto7\nnQ6VSoVGozGdgUnSKcsi+BuNBs1mk93dXarV6iD8O50O1WqV3d1dms2m4S8pC3Mf/P3Q7+uH/5Ur\nVwah32f4S8rBXAd/t9ul3W4Pta2wwv7uPhcuXGB/d58VVob62+22NX9Jc22ug79QKLC5uUmxWAQO\nQn+NNdZZp0SJddZZY20Q/sVikc3NTQqFwjSHLUkTNdfBD1Aulwfhv8UWO+xQosQlLlGixA47bLE1\nCP1yuTztIUvSRM198MNB+LdaLfbYo0lzqK9Jkz32aLVahr6kLIwc/BHxcERsRsR3I+LFiFjrtd8X\nEc9GxMu93+8efbh3p9PpsLq6yiKL1KkP9dWps8giq6urNy31lKR5NI4j/reB300pPQp8GPhkRDwK\nPAU8l1J6BHiu9/zUHV6yuczyoLxzkYuDss8yyzct9ZSkeRUppfG+YcRXgT/o/SynlF6LiAeBrZTS\nL9zutUtLS2l7e3tsY+l2u1QqlaElmyussMUWe+yxyCLLLLPBxqC/WCxy9epVv+CVNDMi4vmU0tJx\ntx9rjT8iSsAHgW8CD6SUXut1/Rh44IjXrEbEdkRsX7t2bZzDoVAoUKvVhto22GChuMDly5dZKC4M\nhT5ArVYz9CXNtbEFf0TcC/wJ8KmU0t8d7ksHf1bc8k+LlFIrpbSUUlo6d+7cuIYz0Gg0qNev1/X7\nq3fOnz8/tNQToF6vewKXpLn3jnG8SUS8k4PQ/2JK6U97za9HxIOHSj1vjGNfd6Mf5u12e2jJZn+p\nZ7VapVarGfqSsjByjT8iAnga+ElK6VOH2n8f6KaUPhcRTwH3pZR+73bvNe4a/4263e4tyzhHtUvS\nLDhpjX8cR/z/HPhXwP+JiO/02v498DngyxHxceAHwL8cw76OVPz8m4PHu5+495bbHBXuhr6knIwc\n/Cml/wnEEd2/Mur7S5LGK4szdyVJ143ly92z4KjyjiRpmEf8kpQZg1+SMmPwS1JmDH5JyozBL0mZ\nMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmD\nX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZGUvw\nR8QXIuKNiHjhUNt9EfFsRLzc+/3ucexLd9btdk/ULikv4zri/yPgsRvangKeSyk9AjzXe64JazQa\nVCoVOp3OUHun06FSqdBoNKYzMElnxliCP6X0DeAnNzQ/Djzde/w0sDKOfelojUaDZrPJ7u4u1Wp1\nEP6dTodqtcru7i7NZtPwlzI3yRr/Ayml13qPfww8MMF9Za8f+n398L9y5cog9PsMfylvp/Llbkop\nAelWfRGxGhHbEbF97dq10xjO3Ol2u7Tb7aG2FVbY393nwoUL7O/us3LDH1ztdtuav5SpSQb/6xHx\nIEDv9xu32iil1EopLaWUls6dOzfB4cyvQqHA5uYmxWIROAj9NdZYZ50SJdZZZ421QfgXi0U2Nzcp\nFArTHLakKZlk8D8DPNl7/CTw1QnuK3vlcnkQ/ltsscMOJUpc4hIlSuywwxZbg9Avl8vTHrKkKRnX\ncs4/Bv4X8AsR8WpEfBz4HPDRiHgZ+EjvuSaoXC7TarXYY48mzaG+Jk322KPVahn6UubioPx+Niwt\nLaXt7e1pD2Nm9Vfv7O/uD8o8fTvs8Dv8DgvFBY/4pTkTEc+nlJaOu71n7s6Jw0s2l1kelHcucnFQ\n9llm+aalnpLy4xH/HOh2u1QqlaElmyussMUWe+yxyCLLLLPBxqC/WCxy9epVv+CV5oBH/BkqFArU\narWhtg02WCgucPnyZRaKC0OhD1Cr1Qx9KVMG/5xoNBrU6/XB8/7qnfPnzw8t9QSo1+uewCVl7B3T\nHoDGpx/m7XZ76Avc/lLParVKrVYz9KXMWeOfQ91u95ZlnKPaJc22k9b4PeKfMb/+5l8MHj9z7wdv\nuc1R4W7oSwJr/JKUHYNfkjJjqWfGHFXekaTj8ohfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbg\nl6TMGPzSHOh2uydqV94MfmnGNRoNKpXKTXdV63Q6VCoVr8aqmxj80gxrNBo0m82bbql5+FaczWbT\n8NcQg1+aUf3Q7+uH/5UrVwah32f46zCDX5pB3W6Xdrs91LbCCvu7+1y4cIH93X1WWBnqb7fb1vwF\nGPzSTCoUCkO31FxhhTXWWGedEiXWWWeNtUH492/F6T0ZBAa/NLP6t9QsFotsscUOO5QocYlLlCix\nww5bbA1Cv38rTsngl2ZYuVym1Wqxxx5NmkN9TZrssUer1TL0NcTgl2ZYp9NhdXWVRRapUx/qq1Nn\nkUVWV1dvWuqpvBn80ow6vGRzmeVBeeciFwdln2WWb1rqKUVKadpjGFhaWkrb29vTHoZ05nW7XSqV\nytCSzRVW2GKLPfZYZJFlltlgY9BfLBa5evWqX/DOoYh4PqW0dNztPeKXZlChUKBWqw21bbDBQnGB\ny5cvs1BcGAp9gFqtZugLMPilmdVoNKjXr9f1+6t3zp8/P7TUE6Ber3sClwa82bo0w/ph3m63h5Zs\n9pd6VqtVarWaoa8h1vilOdDtdm9ZxjmqXfPlzNX4I+KxiHgpIl6JiKcmvT9pXjTj+s+dHBXuhr5u\nZaLBHxH3AH8IfAx4FPjtiHh0kvuUJN3epI/4PwS8klL665TS3wNfAh6f8D4lSbcx6S933wP88NDz\nV4F/OuF9SnOhfna+ftOcmfpyzohYjYjtiNi+du3atIcjSXNv0sH/I+DhQ88f6rUNpJRaKaWllNLS\nuXPnJjwcSdKkg//bwCMR8b6I+BngCeCZCe9TknQbE63xp5Tejoh/B3wNuAf4QkrpxUnuU5J0exM/\nczel9GfAn016P5Kk45n6l7uSpNNl8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmDX5Iy\nY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZMfglKTMG\nvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlJmRgj8ifisi\nXoyIf4iIpRv6Ph0Rr0TESxHxq6MNU5I0Lu8Y8fUvAL8J/JfDjRHxKPAE8ItAEfh6RJRTSj8dcX+S\npBGNdMSfUvpeSumlW3Q9DnwppfRWSulvgFeAD42yL0nSeEyqxv8e4IeHnr/aa7tJRKxGxHZEbF+7\ndm1Cw5Ek9d2x1BMRXwd+/hZdn0kpfXXUAaSUWkALYGlpKY36fpKk27tj8KeUPnIX7/sj4OFDzx/q\ntUmSpmxSpZ5ngCci4l0R8T7gEeBbE9qXJOkERl3O+RsR8Srwz4ArEfE1gJTSi8CXge8C/x34pCt6\nJOlsGGk5Z0rpK8BXjuj7LPDZUd5fkjR+nrkrSVPS7XZP1D4uBr8kTUGj0aBSqdDpdIbaO50OlUqF\nRqMxsX0b/JJ0yhqNBs1mk93dXarV6iD8O50O1WqV3d1dms3mxMLf4JekU9QP/b5++F+5cmUQ+n2T\nCn+DX5JOSbfbpd1uD7WtsML+7j4XLlxgf3efFVaG+tvt9thr/ga/JJ2SQqHA5uYmxWIROAj9NdZY\nZ50SJdZZZ421QfgXi0U2NzcpFApjHYfBL0mnqFwuD8J/iy122KFEiUtcokSJHXbYYmsQ+uVyeexj\nMPgl6ZSVy2VarRZ77NGkOdTXpMkee7RarYmEPhj8knTqOp0Oq6urLLJInfpQX506iyyyurp601LP\ncTH4JekUHV6yuczyoLxzkYuDss8yyzct9RynSOnsXAl5aWkpbW9vT3sYkjQR3W6XSqUytGRzhRW2\n2GKPPRZZZJllNtgY9BeLRa5evXrbL3gj4vmU0tKRG9zAI35JOiWFQoFarTbUtsEGC8UFLl++zEJx\nYSj0AWq1mqt6JGmWNRoN6vXrdf3+6p3z588PLfUEqNfrEzmBa9SbrUuSTqgf5u12e2jJZn+pZ7Va\npVarTeySDdb4JWlKut3uLcs4R7Uf5aQ1fo/4T2hc/1CS5tflNz8/eHzh3k8cud1RmTHpLLHGfwLT\nvIyqJI2LwX9M076MqiSNi6WeYzjqMqqtVovV1dWbLqPaf42kPN2uvHMWeMR/B2flMqqSNC4G/x2c\nlcuoStK4GPzHcBYuoypJ42LwH9O0L6MqSeNi8B/TtC+jKknjYvAfw1m4jKokjYuXbLiDSV1GVZLG\nxcsyj9lZuYyqJI2LwX8MZ+EyqpI0Lp65e0zTvoyqJI2LNf6+iOuPb/OZeHVOSWeNNf4Jm9ZlVCVp\nXEYK/oj4/Yj4fkRcjYivRMTiob5PR8QrEfFSRPzq6EOVJI3DqEf8zwLvTylVgA7waYCIeBR4AvhF\n4DHgP0XEPSPua7JSuv4jSXNspOBPKf2PlNLbvad/DjzUe/w48KWU0lsppb8BXgE+NMq+JEnjMc4a\n/78B/lvv8XuAHx7qe7XXJkmasjsu54yIrwM/f4uuz6SUvtrb5jPA28AXTzqAiFgFVgHe+973nvTl\nkqQTumPwp5Q+crv+iPjXwAXgV9L1taE/Ah4+tNlDvbZbvX8LaMHBcs47D1mSNIpRV/U8Bvwe8Osp\npf1DXc8AT0TEuyLifcAjwLdG2ZckaTxGPXP3D4B3Ac/GwQlQf55S+rcppRcj4svAdzkoAX0ypfTT\nEfclSRqDkYI/pfRPbtP3WeCzo7y/JGn8PHNXkjJzpq7VExHXgB+c4CX3A387oeHMitw/A+ef9/zB\nz+B+4B+nlM4d9wVnKvhPKiK2T3JhonmU+2fg/POeP/gZ3M38LfVIUmYMfknKzKwHf2vaAzgDcv8M\nnL9y/wxOPP+ZrvFLkk5u1o/4JUknNJPBn/sNYCLityLixYj4h4hYuqFv7ucPB5cL6c3xlYh4atrj\nOQ0R8YWIeCMiXjjUdl9EPBsRL/d+v3uaY5ykiHg4IjYj4ru9//7Xeu1ZfAYR8bMR8a2I+Mve/Ju9\n9hPPfyaDn3m6AczdeQH4TeAbhxtzmX9vTn8IfAx4FPjt3tzn3R9x8O962FPAcymlR4Dnes/n1dvA\n76aUHgU+DHyy9++ey2fwFvDLKaVfAj4APBYRH+Yu5j+TwZ/7DWBSSt9LKb10i64s5s/BnF5JKf11\nSunvgS9xMPe5llL6BvCTG5ofB57uPX4aWDnVQZ2ilNJrKaX/3Xv8f4HvcXCfjyw+g3Tgzd7Td/Z+\nEncx/5kM/ht4A5jrcpl/LvM8jgdSSq/1Hv8YeGCagzktEVECPgh8k4w+g4i4JyK+A7wBPJtSuqv5\nj3p1zomZ9A1gzrrjzF86LKWUImLul+lFxL3AnwCfSin9Xe/KwMD8fwa9qxx/oPe95lci4v039B9r\n/mc2+Cd9A5iz7k7zP8LczP8OcpnncbweEQ+mlF6LiAc5OBKcWxHxTg5C/4sppT/tNWf1GQCklPYi\nYpOD73xOPP+ZLPV4A5gj5TL/bwOPRMT7IuJnOPhC+5kpj2langGe7D1+Epjbvwbj4ND+vwLfSyn9\nx0NdWXwGEXGuv4IxIn4O+Cjwfe5m/imlmfvh4EvLHwLf6f3850N9nwH+CngJ+Ni0xzqh+f8GB3Xt\nt4DXga/lNP/ePH+NgxVdf8VB+WvqYzqFOf8x8Brw/3r//h8HChys5HgZ+Dpw37THOcH5/wsOvsy8\neuj//V/L5TMAKsBf9Ob/AvAfeu0nnr9n7kpSZmay1CNJunsGvyRlxuCXpMwY/JKUGYNfkjJj8EtS\nZgx+ScqMwS9Jmfn/nGYAT84PDosAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<matplotlib.figure.Figure at 0x7f5004dfb5c0>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"plot_data(centroids+2, X, n_samples)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Broadcasting"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"How did our distance function `sqrt(((x-X)**2).sum(1))` work over a matrix without us writing any loops? The trick is that we used *broadcasting*. The term broadcasting was first used by Numpy, although is now used in other libraries such as [Tensorflow](https://www.tensorflow.org/performance/xla/broadcasting) and Matlab; the rules can vary by library.\n",
|
|
"\n",
|
|
"From the [Numpy Documentation](https://docs.scipy.org/doc/numpy-1.10.0/user/basics.broadcasting.html):\n",
|
|
"\n",
|
|
" The term broadcasting describes how numpy treats arrays with \n",
|
|
" different shapes during arithmetic operations. Subject to certain \n",
|
|
" constraints, the smaller array is “broadcast” across the larger \n",
|
|
" array so that they have compatible shapes.Broadcasting provides a \n",
|
|
" means of vectorizing array operations so that looping occurs in C\n",
|
|
" instead of Python. It does this without making needless copies of \n",
|
|
" data and usually leads to efficient algorithm implementations.\n",
|
|
" \n",
|
|
"In addition to the efficiency of broadcasting, it allows developers to write less code, which typically leads to fewer errors.\n",
|
|
"\n",
|
|
"Operators (+,-,\\*,/,>,<,==) are usually element-wise. Here's some examples of element-wise operations:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 46,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(array([12, 14, 3]), array([False, True, True], dtype=bool))"
|
|
]
|
|
},
|
|
"execution_count": 46,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"a = np.array([10, 6, -4])\n",
|
|
"b = np.array([2, 8, 7])\n",
|
|
"\n",
|
|
"a + b, a < b"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Now this next example clearly can't be element-wise, since the second parameter is a scalar, not a 1d array."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 47,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([ True, True, False], dtype=bool)"
|
|
]
|
|
},
|
|
"execution_count": 47,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"a > 0"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"So how did this work? The trick was that numpy automatically *broadcast* the scalar `0` so that had the same `shape` as a. We can manually see how numpy broadcasts by using `broadcast_to()`."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 48,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(3,)"
|
|
]
|
|
},
|
|
"execution_count": 48,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"a.shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 49,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([0, 0, 0])"
|
|
]
|
|
},
|
|
"execution_count": 49,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"np.broadcast_to(0, a.shape)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 50,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(3,)"
|
|
]
|
|
},
|
|
"execution_count": 50,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"np.broadcast_to(0, a.shape).shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here's another example."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 51,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([11, 7, -3])"
|
|
]
|
|
},
|
|
"execution_count": 51,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"a + 1"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"It works with higher-dimensional arrays too, for instance 2d (matrices):"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 52,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[1, 2, 3],\n",
|
|
" [4, 5, 6],\n",
|
|
" [7, 8, 9]])"
|
|
]
|
|
},
|
|
"execution_count": 52,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"m = np.array([[1, 2, 3], [4,5,6], [7,8,9]]); m"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 53,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[ 2, 4, 6],\n",
|
|
" [ 8, 10, 12],\n",
|
|
" [14, 16, 18]])"
|
|
]
|
|
},
|
|
"execution_count": 53,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"m * 2"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 54,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"(3, 3)"
|
|
]
|
|
},
|
|
"execution_count": 54,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"m.shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 55,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[2, 2, 2],\n",
|
|
" [2, 2, 2],\n",
|
|
" [2, 2, 2]])"
|
|
]
|
|
},
|
|
"execution_count": 55,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"np.broadcast_to(2, m.shape)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can use the same trick to broadcast a vector to a matrix:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 56,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([10, 20, 30])"
|
|
]
|
|
},
|
|
"execution_count": 56,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"c = np.array([10,20,30]); c"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 57,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[11, 22, 33],\n",
|
|
" [14, 25, 36],\n",
|
|
" [17, 28, 39]])"
|
|
]
|
|
},
|
|
"execution_count": 57,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"m + c"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Let's see what numpy has done with `c` in this case:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 58,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[10, 20, 30],\n",
|
|
" [10, 20, 30],\n",
|
|
" [10, 20, 30]])"
|
|
]
|
|
},
|
|
"execution_count": 58,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"np.broadcast_to(c, m.shape)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Interesting - we see that it has duplicated `c` across rows. What if `c` was a column vector, i.e. a 3x1 array?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 59,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[10],\n",
|
|
" [20],\n",
|
|
" [30]])"
|
|
]
|
|
},
|
|
"execution_count": 59,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"# Indexing an axis with None adds a unit axis in that location\n",
|
|
"cc = c[:,None]; cc"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 60,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[11, 12, 13],\n",
|
|
" [24, 25, 26],\n",
|
|
" [37, 38, 39]])"
|
|
]
|
|
},
|
|
"execution_count": 60,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"m + cc"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Let's see what numpy has done with `c` in this case:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 61,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[10, 10, 10],\n",
|
|
" [20, 20, 20],\n",
|
|
" [30, 30, 30]])"
|
|
]
|
|
},
|
|
"execution_count": 61,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"np.broadcast_to(cc, m.shape)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"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": {},
|
|
"source": [
|
|
"We can now see how our `distance()` function works:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 62,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[1, 1],\n",
|
|
" [0, 0],\n",
|
|
" [9, 4]])"
|
|
]
|
|
},
|
|
"execution_count": 62,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"a=array([2,3])\n",
|
|
"b=array([[1,2],[2,3],[-1,1]])\n",
|
|
"c=(a-b)**2; c"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 64,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[ 1, 2],\n",
|
|
" [ 2, 3],\n",
|
|
" [-1, 1]])"
|
|
]
|
|
},
|
|
"execution_count": 64,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"b"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 63,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([ 0.13598247, 0.15957691, 0.05640321])"
|
|
]
|
|
},
|
|
"execution_count": 63,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"w=gaussian(sqrt(c.sum(1)), 2.5); w"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"...and we can also now pull apart our weighted average:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 183,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"((3,), (3, 2))"
|
|
]
|
|
},
|
|
"execution_count": 183,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"w.shape, b.shape"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 184,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[ 0.13598247],\n",
|
|
" [ 0.15957691],\n",
|
|
" [ 0.05640321]])"
|
|
]
|
|
},
|
|
"execution_count": 184,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"w[:,None]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 185,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[ 0.13598247, 0.27196495],\n",
|
|
" [ 0.31915382, 0.47873074],\n",
|
|
" [-0.05640321, 0.05640321]])"
|
|
]
|
|
},
|
|
"execution_count": 185,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"w[:,None]*b"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 186,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([ 1.13288486, 2.29313827])"
|
|
]
|
|
},
|
|
"execution_count": 186,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"(w[:,None]*b).sum(0) / w.sum()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## GPU-accelerated mean shift in pytorch"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Let's now look at using [PyTorch](http://pytorch.org/), a Python framework for dynamic neural networks with GPU acceleration, which was released by Facebook's AI team.\n",
|
|
"\n",
|
|
"PyTorch has two overlapping, yet distinct, purposes. As described in the [PyTorch documentation](http://pytorch.org/tutorials/beginner/blitz/tensor_tutorial.html):\n",
|
|
"\n",
|
|
"<img src=\"images/what_is_pytorch.png\" alt=\"pytorch\" style=\"width: 80%\"/>\n",
|
|
"\n",
|
|
"The neural network functionality of PyTorch is built on top of the Numpy-like functionality for fast matrix computations on a GPU. Although the neural network purpose receives much more attention, both are very useful. Today we'll use PyTorch to accelerate our meanshift algorithm by running it on the GPU.\n",
|
|
"\n",
|
|
"If you want to learn more PyTorch, you can try this [introductory tutorial](http://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html) or this [tutorial to learn by examples](http://pytorch.org/tutorials/beginner/pytorch_with_examples.html).\n",
|
|
"\n",
|
|
"One advantage of pytorch is that it's very similar to numpy. For instance, in fact, our definitions of `gaussian` and `distance` and `meanshift_iter` are identical. So we'll simply import PyTorch's alternate implementations of the two numpy functions we use:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 66,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"from torch import exp, sqrt"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"And then we'll use the exact same code as before, but first convert our numpy array to a GPU PyTorch tensor."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 72,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def meanshift(data):\n",
|
|
" X = torch.from_numpy(np.copy(data)).cuda()\n",
|
|
" for it in range(5): X = meanshift_iter(X)\n",
|
|
" return X"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Let's try it out..."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 73,
|
|
"metadata": {
|
|
"scrolled": false
|
|
},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 1.59 s, sys: 8 ms, total: 1.6 s\n",
|
|
"Wall time: 1.63 s\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSJJREFUeJzt3VGMY9d93/Hvv7LjZKqHsaiFIloyaKBiAMVgbGDgumgf\nhomNyNlFNAmQQnko1K3BQQ0XGAcBArkGSvLBgIEAnZekbknUGz0YMQwkHgu7aV1ZmIFRoLE9apyt\nZFuUkoxheWRpQ2OQCgMolXP6MCR3uLuzO7Mkh0Oe7wcYDHnOJe85XOmHO3+ee2+klJAk5eMfTXsA\nkqTTZfBLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMvOOaQ/gsPvvvz+VSqVpD0OS\nZsrzzz//tymlc8fd/kwFf6lUYnt7e9rDkKSZEhE/OMn2lnokKTMGv6Q76na7J2rX2WbwS7qtRqNB\npVKh0+kMtXc6HSqVCo1GYzoD010bOfgj4mcj4lsR8ZcR8WJENHvt90XEsxHxcu/3u0cfrqTT1Gg0\naDab7O7uUq1WB+Hf6XSoVqvs7u7SbDYN/xkzjiP+t4BfTin9EvAB4LGI+DDwFPBcSukR4Lnec0kz\noh/6ff3wv3LlyiD0+wz/2TJy8KcDb/aevrP3k4DHgad77U8DK6PuS9Lp6Ha7tNvtobYVVtjf3efC\nhQvs7+6zcsP/0u1225r/jBhLjT8i7omI7wBvAM+mlL4JPJBSeq23yY+BB8axL0mTVygU2NzcpFgs\nAgehv8Ya66xTosQ666yxNgj/YrHI5uYmhUJhmsPWMY0l+FNKP00pfQB4CPhQRLz/hv7EwV8BN4mI\n1YjYjojta9eujWM4ksagXC4Pwn+LLXbYoUSJS1yiRIkddthiaxD65XJ52kPWMY11VU9KaQ/YBB4D\nXo+IBwF6v9844jWtlNJSSmnp3Lljn3gm6RSUy2VarRZ77NGkOdTXpMkee7RaLUN/xoxjVc+5iFjs\nPf454KPA94FngCd7mz0JfHXUfUk6XZ1Oh9XVVRZZpE59qK9OnUUWWV1dvWmpp862cRzxPwhsRsRV\n4Nsc1PgvA58DPhoRLwMf6T2XNCMOL9lcZnlQ3rnIxUHZZ5nlm5Z66uyLg/L72bC0tJS8Vo80fd1u\nl0qlMrRkc4UVtthijz0WWWSZZTbYGPQXi0WuXr3qF7xTEBHPp5SWjru9Z+5KukmhUKBWqw21bbDB\nQnGBy5cvs1BcGAp9gFqtZujPCINf0i01Gg3q9et1/f7qnfPnzw8t9QSo1+uewDVDztRlmSWdLf0w\nb7fbQ0s2+0s9q9UqtVrN0J8x1vil3H3+0Bm4n9i45SbdbveWZZyj2nW6rPFLGrujwt3Qn00GvyRl\nxhq/lLsjyjuaXx7xS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqM\nwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmayCv5ut3uidkmaR9kEf6PRoFKp0Ol0hto7\nnQ6VSoVGozGdgUnSKcsi+BuNBs1mk93dXarV6iD8O50O1WqV3d1dms2m4S8pC3Mf/P3Q7+uH/5Ur\nVwah32f4S8rBXAd/t9ul3W4Pta2wwv7uPhcuXGB/d58VVob62+22NX9Jc22ug79QKLC5uUmxWAQO\nQn+NNdZZp0SJddZZY20Q/sVikc3NTQqFwjSHLUkTNdfBD1Aulwfhv8UWO+xQosQlLlGixA47bLE1\nCP1yuTztIUvSRM198MNB+LdaLfbYo0lzqK9Jkz32aLVahr6kLIwc/BHxcERsRsR3I+LFiFjrtd8X\nEc9GxMu93+8efbh3p9PpsLq6yiKL1KkP9dWps8giq6urNy31lKR5NI4j/reB300pPQp8GPhkRDwK\nPAU8l1J6BHiu9/zUHV6yuczyoLxzkYuDss8yyzct9ZSkeRUppfG+YcRXgT/o/SynlF6LiAeBrZTS\nL9zutUtLS2l7e3tsY+l2u1QqlaElmyussMUWe+yxyCLLLLPBxqC/WCxy9epVv+CVNDMi4vmU0tJx\ntx9rjT8iSsAHgW8CD6SUXut1/Rh44IjXrEbEdkRsX7t2bZzDoVAoUKvVhto22GChuMDly5dZKC4M\nhT5ArVYz9CXNtbEFf0TcC/wJ8KmU0t8d7ksHf1bc8k+LlFIrpbSUUlo6d+7cuIYz0Gg0qNev1/X7\nq3fOnz8/tNQToF6vewKXpLn3jnG8SUS8k4PQ/2JK6U97za9HxIOHSj1vjGNfd6Mf5u12e2jJZn+p\nZ7VapVarGfqSsjByjT8iAnga+ElK6VOH2n8f6KaUPhcRTwH3pZR+73bvNe4a/4263e4tyzhHtUvS\nLDhpjX8cR/z/HPhXwP+JiO/02v498DngyxHxceAHwL8cw76OVPz8m4PHu5+495bbHBXuhr6knIwc\n/Cml/wnEEd2/Mur7S5LGK4szdyVJ143ly92z4KjyjiRpmEf8kpQZg1+SMmPwS1JmDH5JyozBL0mZ\nMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmD\nX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZGUvw\nR8QXIuKNiHjhUNt9EfFsRLzc+/3ucexLd9btdk/ULikv4zri/yPgsRvangKeSyk9AjzXe64JazQa\nVCoVOp3OUHun06FSqdBoNKYzMElnxliCP6X0DeAnNzQ/Djzde/w0sDKOfelojUaDZrPJ7u4u1Wp1\nEP6dTodqtcru7i7NZtPwlzI3yRr/Ayml13qPfww8MMF9Za8f+n398L9y5cog9PsMfylvp/Llbkop\nAelWfRGxGhHbEbF97dq10xjO3Ol2u7Tb7aG2FVbY393nwoUL7O/us3LDH1ztdtuav5SpSQb/6xHx\nIEDv9xu32iil1EopLaWUls6dOzfB4cyvQqHA5uYmxWIROAj9NdZYZ50SJdZZZ421QfgXi0U2Nzcp\nFArTHLakKZlk8D8DPNl7/CTw1QnuK3vlcnkQ/ltsscMOJUpc4hIlSuywwxZbg9Avl8vTHrKkKRnX\ncs4/Bv4X8AsR8WpEfBz4HPDRiHgZ+EjvuSaoXC7TarXYY48mzaG+Jk322KPVahn6UubioPx+Niwt\nLaXt7e1pD2Nm9Vfv7O/uD8o8fTvs8Dv8DgvFBY/4pTkTEc+nlJaOu71n7s6Jw0s2l1kelHcucnFQ\n9llm+aalnpLy4xH/HOh2u1QqlaElmyussMUWe+yxyCLLLLPBxqC/WCxy9epVv+CV5oBH/BkqFArU\narWhtg02WCgucPnyZRaKC0OhD1Cr1Qx9KVMG/5xoNBrU6/XB8/7qnfPnzw8t9QSo1+uewCVl7B3T\nHoDGpx/m7XZ76Avc/lLParVKrVYz9KXMWeOfQ91u95ZlnKPaJc22k9b4PeKfMb/+5l8MHj9z7wdv\nuc1R4W7oSwJr/JKUHYNfkjJjqWfGHFXekaTj8ohfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbg\nl6TMGPzSHOh2uydqV94MfmnGNRoNKpXKTXdV63Q6VCoVr8aqmxj80gxrNBo0m82bbql5+FaczWbT\n8NcQg1+aUf3Q7+uH/5UrVwah32f46zCDX5pB3W6Xdrs91LbCCvu7+1y4cIH93X1WWBnqb7fb1vwF\nGPzSTCoUCkO31FxhhTXWWGedEiXWWWeNtUH492/F6T0ZBAa/NLP6t9QsFotsscUOO5QocYlLlCix\nww5bbA1Cv38rTsngl2ZYuVym1Wqxxx5NmkN9TZrssUer1TL0NcTgl2ZYp9NhdXWVRRapUx/qq1Nn\nkUVWV1dvWuqpvBn80ow6vGRzmeVBeeciFwdln2WWb1rqKUVKadpjGFhaWkrb29vTHoZ05nW7XSqV\nytCSzRVW2GKLPfZYZJFlltlgY9BfLBa5evWqX/DOoYh4PqW0dNztPeKXZlChUKBWqw21bbDBQnGB\ny5cvs1BcGAp9gFqtZugLMPilmdVoNKjXr9f1+6t3zp8/P7TUE6Ber3sClwa82bo0w/ph3m63h5Zs\n9pd6VqtVarWaoa8h1vilOdDtdm9ZxjmqXfPlzNX4I+KxiHgpIl6JiKcmvT9pXjTj+s+dHBXuhr5u\nZaLBHxH3AH8IfAx4FPjtiHh0kvuUJN3epI/4PwS8klL665TS3wNfAh6f8D4lSbcx6S933wP88NDz\nV4F/OuF9SnOhfna+ftOcmfpyzohYjYjtiNi+du3atIcjSXNv0sH/I+DhQ88f6rUNpJRaKaWllNLS\nuXPnJjwcSdKkg//bwCMR8b6I+BngCeCZCe9TknQbE63xp5Tejoh/B3wNuAf4QkrpxUnuU5J0exM/\nczel9GfAn016P5Kk45n6l7uSpNNl8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlBmDX5Iy\nY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZMfglKTMG\nvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8kZcbgl6TMGPySlJmRgj8ifisi\nXoyIf4iIpRv6Ph0Rr0TESxHxq6MNU5I0Lu8Y8fUvAL8J/JfDjRHxKPAE8ItAEfh6RJRTSj8dcX+S\npBGNdMSfUvpeSumlW3Q9DnwppfRWSulvgFeAD42yL0nSeEyqxv8e4IeHnr/aa7tJRKxGxHZEbF+7\ndm1Cw5Ek9d2x1BMRXwd+/hZdn0kpfXXUAaSUWkALYGlpKY36fpKk27tj8KeUPnIX7/sj4OFDzx/q\ntUmSpmxSpZ5ngCci4l0R8T7gEeBbE9qXJOkERl3O+RsR8Srwz4ArEfE1gJTSi8CXge8C/x34pCt6\nJOlsGGk5Z0rpK8BXjuj7LPDZUd5fkjR+nrkrSVPS7XZP1D4uBr8kTUGj0aBSqdDpdIbaO50OlUqF\nRqMxsX0b/JJ0yhqNBs1mk93dXarV6iD8O50O1WqV3d1dms3mxMLf4JekU9QP/b5++F+5cmUQ+n2T\nCn+DX5JOSbfbpd1uD7WtsML+7j4XLlxgf3efFVaG+tvt9thr/ga/JJ2SQqHA5uYmxWIROAj9NdZY\nZ50SJdZZZ421QfgXi0U2NzcpFApjHYfBL0mnqFwuD8J/iy122KFEiUtcokSJHXbYYmsQ+uVyeexj\nMPgl6ZSVy2VarRZ77NGkOdTXpMkee7RarYmEPhj8knTqOp0Oq6urLLJInfpQX506iyyyurp601LP\ncTH4JekUHV6yuczyoLxzkYuDss8yyzct9RynSOnsXAl5aWkpbW9vT3sYkjQR3W6XSqUytGRzhRW2\n2GKPPRZZZJllNtgY9BeLRa5evXrbL3gj4vmU0tKRG9zAI35JOiWFQoFarTbUtsEGC8UFLl++zEJx\nYSj0AWq1mqt6JGmWNRoN6vXrdf3+6p3z588PLfUEqNfrEzmBa9SbrUuSTqgf5u12e2jJZn+pZ7Va\npVarTeySDdb4JWlKut3uLcs4R7Uf5aQ1fo/4T2hc/1CS5tflNz8/eHzh3k8cud1RmTHpLLHGfwLT\nvIyqJI2LwX9M076MqiSNi6WeYzjqMqqtVovV1dWbLqPaf42kPN2uvHMWeMR/B2flMqqSNC4G/x2c\nlcuoStK4GPzHcBYuoypJ42LwH9O0L6MqSeNi8B/TtC+jKknjYvAfw1m4jKokjYuXbLiDSV1GVZLG\nxcsyj9lZuYyqJI2LwX8MZ+EyqpI0Lp65e0zTvoyqJI2LNf6+iOuPb/OZeHVOSWeNNf4Jm9ZlVCVp\nXEYK/oj4/Yj4fkRcjYivRMTiob5PR8QrEfFSRPzq6EOVJI3DqEf8zwLvTylVgA7waYCIeBR4AvhF\n4DHgP0XEPSPua7JSuv4jSXNspOBPKf2PlNLbvad/DjzUe/w48KWU0lsppb8BXgE+NMq+JEnjMc4a\n/78B/lvv8XuAHx7qe7XXJkmasjsu54yIrwM/f4uuz6SUvtrb5jPA28AXTzqAiFgFVgHe+973nvTl\nkqQTumPwp5Q+crv+iPjXwAXgV9L1taE/Ah4+tNlDvbZbvX8LaMHBcs47D1mSNIpRV/U8Bvwe8Osp\npf1DXc8AT0TEuyLifcAjwLdG2ZckaTxGPXP3D4B3Ac/GwQlQf55S+rcppRcj4svAdzkoAX0ypfTT\nEfclSRqDkYI/pfRPbtP3WeCzo7y/JGn8PHNXkjJzpq7VExHXgB+c4CX3A387oeHMitw/A+ef9/zB\nz+B+4B+nlM4d9wVnKvhPKiK2T3JhonmU+2fg/POeP/gZ3M38LfVIUmYMfknKzKwHf2vaAzgDcv8M\nnL9y/wxOPP+ZrvFLkk5u1o/4JUknNJPBn/sNYCLityLixYj4h4hYuqFv7ucPB5cL6c3xlYh4atrj\nOQ0R8YWIeCMiXjjUdl9EPBsRL/d+v3uaY5ykiHg4IjYj4ru9//7Xeu1ZfAYR8bMR8a2I+Mve/Ju9\n9hPPfyaDn3m6AczdeQH4TeAbhxtzmX9vTn8IfAx4FPjt3tzn3R9x8O962FPAcymlR4Dnes/n1dvA\n76aUHgU+DHyy9++ey2fwFvDLKaVfAj4APBYRH+Yu5j+TwZ/7DWBSSt9LKb10i64s5s/BnF5JKf11\nSunvgS9xMPe5llL6BvCTG5ofB57uPX4aWDnVQZ2ilNJrKaX/3Xv8f4HvcXCfjyw+g3Tgzd7Td/Z+\nEncx/5kM/ht4A5jrcpl/LvM8jgdSSq/1Hv8YeGCagzktEVECPgh8k4w+g4i4JyK+A7wBPJtSuqv5\nj3p1zomZ9A1gzrrjzF86LKWUImLul+lFxL3AnwCfSin9Xe/KwMD8fwa9qxx/oPe95lci4v039B9r\n/mc2+Cd9A5iz7k7zP8LczP8OcpnncbweEQ+mlF6LiAc5OBKcWxHxTg5C/4sppT/tNWf1GQCklPYi\nYpOD73xOPP+ZLPV4A5gj5TL/bwOPRMT7IuJnOPhC+5kpj2langGe7D1+Epjbvwbj4ND+vwLfSyn9\nx0NdWXwGEXGuv4IxIn4O+Cjwfe5m/imlmfvh4EvLHwLf6f3850N9nwH+CngJ+Ni0xzqh+f8GB3Xt\nt4DXga/lNP/ePH+NgxVdf8VB+WvqYzqFOf8x8Brw/3r//h8HChys5HgZ+Dpw37THOcH5/wsOvsy8\neuj//V/L5TMAKsBf9Ob/AvAfeu0nnr9n7kpSZmay1CNJunsGvyRlxuCXpMwY/JKUGYNfkjJj8EtS\nZgx+ScqMwS9Jmfn/nGYAT84PDosAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<matplotlib.figure.Figure at 0x7f5004c588d0>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"%time X = meanshift(data).cpu().numpy()\n",
|
|
"plot_data(centroids+2, X, n_samples)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"It works, but this implementation actually takes longer. Oh dear! What do you think is causing this?\n",
|
|
"\n",
|
|
"Each iteration launches a new cuda kernel, which takes time and slows the algorithm down as a whole. Furthermore, each iteration doesn't have enough processing to do to fill up all of the threads of the GPU. To use the GPU effectively, we need to process a *batch* of data at a time."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## GPU batched algorithm"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"To process a batch of data, we need batched versions of our functions. Here's a version of `distance()` that works on batches:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 74,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def distance_b(a,b): return sqrt(((a[None,:] - b[:,None]) ** 2).sum(2))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 75,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"\n",
|
|
" 0.4714 0.6740 0.9486\n",
|
|
" 0.5486 0.7749 0.5401\n",
|
|
"[torch.FloatTensor of size 2x3]"
|
|
]
|
|
},
|
|
"execution_count": 75,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"a=torch.rand(2,2)\n",
|
|
"b=torch.rand(3,2)\n",
|
|
"distance_b(b, a)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Note how the two parameters to `distance_b()` have a unit axis added in two different places (`a[None,:]` and `b[:,None]`). This is a handy trick which effectively generalizes the concept of an 'outer product' to any function. In this case, we use it to get the distance from every point in `a` (our batch) to every point in `b` (the whole dataset).\n",
|
|
"\n",
|
|
"Now that we have a suitable distance function, we can make some minor updates to our meanshift function to handle batches of data:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 76,
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"def meanshift(data, bs=500):\n",
|
|
" n = len(data)\n",
|
|
" X = torch.FloatTensor(np.copy(data)).cuda()\n",
|
|
" for it in range(5):\n",
|
|
" for i in range(0,n,bs):\n",
|
|
" s = slice(i,min(n,i+bs))\n",
|
|
" weight = gaussian(distance_b(X, X[s]), 2.5)\n",
|
|
" num = (weight[:,:,None] * X).sum(1)\n",
|
|
" X[s] = num / weight.sum(1)[:,None]\n",
|
|
" return X"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 80,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[ 9.22296143, 2.36360407],\n",
|
|
" [ 9.3774662 , -0.18868986],\n",
|
|
" [ 9.21319103, 0.2641083 ],\n",
|
|
" ..., \n",
|
|
" [-16.89390182, -18.92529297],\n",
|
|
" [-17.19715118, -19.02577019],\n",
|
|
" [-16.48303604, -19.09796906]], dtype=float32)"
|
|
]
|
|
},
|
|
"execution_count": 80,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"data"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 81,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"\n",
|
|
" 9.2230 2.3636\n",
|
|
" 9.3775 -0.1887\n",
|
|
" 9.2132 0.2641\n",
|
|
" ⋮ \n",
|
|
"-16.8939 -18.9253\n",
|
|
"-17.1972 -19.0258\n",
|
|
"-16.4830 -19.0980\n",
|
|
"[torch.cuda.FloatTensor of size 1500x2 (GPU 0)]"
|
|
]
|
|
},
|
|
"execution_count": 81,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"torch.FloatTensor(np.copy(data)).cuda()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Although each iteration still has to launch a new cuda kernel, there are now fewer iterations, and the acceleration from updating a batch of points more than makes up for it."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 78,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"CPU times: user 24 ms, sys: 4 ms, total: 28 ms\n",
|
|
"Wall time: 26.4 ms\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%time X = meanshift(data).cpu().numpy()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"That's more like it! We've gone from 2000ms to 26ms, which is a speedup of over 7000%. Oh, and it even gives the right answer!"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 79,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFSBJREFUeJzt3WFs4/d93/H3d06aVvMDxfTBNWMfGGBmATdgE0DIMmwP\nxDZBnd6hVgt0cB8M3i2gsCADlKJA4SzASD4IEKDA9KRdNhLL1Q+CBgHayMZdt8wxJAQD1iTymmp2\nEtNuqyCOHPvKQOgMAe6c/vZAJE+8O91JJ1IU+Xu/AEHk7/cn/78fz/7gry9///8/UkpIkvLxjyY9\nAEnS6TL4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmYMfknKjMEvSZl516QHcND999+fSqXSpIch\nSVPlhRde+NuU0rmjbn+mgr9UKrG5uTnpYUjSVImIHxxne0s9kpQZg1/SHXW73WO162wz+CXdVqPR\noFKp0Ol0hto7nQ6VSoVGozGZgemunTj4I+JnI+JbEfGXEfFSRDR77fdFxHMR8Urv93tPPlxJp6nR\naNBsNtnZ2aFarQ7Cv9PpUK1W2dnZodlsGv5TZhRH/G8Dv5xS+iXgg8BjEfER4Cng+ZTSI8DzveeS\npkQ/9Pv64X/16tVB6PcZ/tPlxMGf9r3Ve/ru3k8CHgee7rU/DSyddF+STke326Xdbg+1LbHE3s4e\nFy9eZG9nj6Ub/pdut9vW/KfESGr8EXFPRHwHeBN4LqX0TeCBlNLrvU1+DDwwin1JGr9CocD6+jrF\nYhHYD/0VVlhllRIlVlllhZVB+BeLRdbX1ykUCpMcto5oJMGfUvppSumDwEPAhyPiAzf0J/b/CrhJ\nRCxHxGZEbF67dm0Uw5E0AuVyeRD+G2ywzTYlSlzmMiVKbLPNBhuD0C+Xy5Meso5opKt6Ukq7wDrw\nGPBGRDwI0Pv95iGvaaWUFlJKC+fOHfnEM0mnoFwu02q12GWXJs2hviZNdtml1WoZ+lNmFKt6zkXE\nfO/xzwEfA74PPAs82dvsSeCZk+5L0unqdDosLy8zzzx16kN9derMM8/y8vJNSz11to3iiP9BYD0i\ntoBvs1/jvwJ8HvhYRLwCfLT3XNKUOLhkc5HFQXnnEpcGZZ9FFm9a6qmzL/bL72fDwsJC8lo90uR1\nu10qlcrQks0llthgg112mWeeRRZZY23QXywW2dra8gveCYiIF1JKC0fd3jN3Jd2kUChQq9WG2tZY\nY644x5UrV5grzg2FPkCtVjP0p4TBL+mWGo0G9fr1un5/9c6FCxeGlnoC1Ot1T+CaImfqssySzpZ+\nmLfb7aElm/2lntVqlVqtZuhPGWv8Us6+cODs20+uHbpZt9u9ZRnnsHadLmv8kkbusHA39KeTwS9J\nmbHGL+XsNuUdzS6P+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPwS1JmDH5JyozBL0mZMfglKTMGvyRl\nxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzWQV/t9s9VrskzaJsgr/RaFCpVOh0OkPt\nnU6HSqVCo9GYzMAk6ZRlEfyNRoNms8nOzg7VanUQ/p1Oh2q1ys7ODs1m0/CXlIWZD/5+6Pf1w//q\n1auD0O8z/CXlYKaDv9vt0m63h9qWWGJvZ4+LFy+yt7PHEktD/e1225q/pJk208FfKBRYX1+nWCwC\n+6G/wgqrrFKixCqrrLAyCP9iscj6+jqFQmGSw5aksZrp4Acol8uD8N9gg222KVHiMpcpUWKbbTbY\nGIR+uVye9JAlaaxmPvhhP/xbrRa77NKkOdTXpMkuu7RaLUNfUhZOHPwR8XBErEfEdyPipYhY6bXf\nFxHPRcQrvd/vPflw706n02F5eZl55qlTH+qrU2eeeZaXl29a6ilJs2gUR/zvAL+bUnoU+AjwqYh4\nFHgKeD6l9AjwfO/5qTu4ZHORxUF55xKXBmWfRRZvWuopSbMqUkqjfcOIZ4A/6P0sppRej4gHgY2U\n0i/c7rULCwtpc3NzZGPpdrtUKpWhJZtLLLHBBrvsMs88iyyyxtqgv1gssrW15Re8kqZGRLyQUlo4\n6vYjrfFHRAn4EPBN4IGU0uu9rh8DDxzymuWI2IyIzWvXro1yOBQKBWq12lDbGmvMFee4cuUKc8W5\nodAHqNVqhr6kmTay4I+Ie4E/AT6dUvq7g31p/8+KW/5pkVJqpZQWUkoL586dG9VwBhqNBvX69bp+\nf/XOhQsXhpZ6AtTrdU/gkjTz3jWKN4mId7Mf+l9KKf1pr/mNiHjwQKnnzVHs6270w7zdbg8t2ewv\n9axWq9RqNUNfUhZOXOOPiACeBn6SUvr0gfbfB7oppc9HxFPAfSml37vde426xn+jbrd7yzLOYe2S\nNA2OW+MfxRH/Pwf+FfB/IuI7vbZ/D3we+EpEfAL4AfAvR7CvQxW/8Nbg8c4n773lNoeFu6EvKScn\nDv6U0v8E4pDuXznp+0uSRiuLM3clSdeN5Mvds+Cw8o4kaZhH/JKUGYNfkjJj8EtSZgx+ScqMwS9J\nmTH4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZ\ng1+SMmPwS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmRlJ\n8EfEFyPizYh48UDbfRHxXES80vv93lHsS3fW7XaP1S4pL6M64v8j4LEb2p4Cnk8pPQI833uuMWs0\nGlQqFTqdzlB7p9OhUqnQaDQmMzBJZ8ZIgj+l9A3gJzc0Pw483Xv8NLA0in3pcI1Gg2azyc7ODtVq\ndRD+nU6HarXKzs4OzWbT8JcyN84a/wMppdd7j38MPDDGfWWvH/p9/fC/evXqIPT7DH8pb6fy5W5K\nKQHpVn0RsRwRmxGxee3atdMYzszpdru02+2htiWW2NvZ4+LFi+zt7LF0wx9c7Xbbmr+UqXEG/xsR\n8SBA7/ebt9oopdRKKS2klBbOnTs3xuHMrkKhwPr6OsViEdgP/RVWWGWVEiVWWWWFlUH4F4tF1tfX\nKRQKkxy2pAkZZ/A/CzzZe/wk8MwY95W9crk8CP8NNthmmxIlLnOZEiW22WaDjUHol8vlSQ9Z0oSM\najnnHwP/C/iFiHgtIj4BfB74WES8Any091xjVC6XabVa7LJLk+ZQX5Mmu+zSarUMfSlzsV9+PxsW\nFhbS5ubmpIcxtfqrd/Z29gZlnr5ttvkdfoe54pxH/NKMiYgXUkoLR93eM3dnxMElm4ssDso7l7g0\nKPsssnjTUk9J+fGIfwZ0u10qlcrQks0llthgg112mWeeRRZZY23QXywW2dra8gteaQZ4xJ+hQqFA\nrVYbaltjjbniHFeuXGGuODcU+gC1Ws3QlzJl8M+IRqNBvV4fPO+v3rlw4cLQUk+Aer3uCVxSxt41\n6QFodPph3m63h77A7S/1rFar1Go1Q1/KnDX+GdTtdm9ZxjmsXdJ0O26N3yP+KfPrb/3F4PGz937o\nltscFu6GviSwxi9J2TH4JSkzlnqmzGHlHUk6Ko/4JSkzBr8kZcbgl6TMGPySlBmDX5IyY/BLUmYM\nfknKjMEvzYBut3usduXN4JemXKPRoFKp3HRXtU6nQ6VS8WqsuonBL02xRqNBs9m86ZaaB2/F2Ww2\nDX8NMfilKdUP/b5++F+9enUQ+n2Gvw4y+KUp1O12abfbQ21LLLG3s8fFixfZ29ljiaWh/na7bc1f\ngMEvTaVCoTB0S80lllhhhVVWKVFilVVWWBmEf/9WnN6TQWDwS1Orf0vNYrHIBhtss02JEpe5TIkS\n22yzwcYg9Pu34pQMfmmKlctlWq0Wu+zSpDnU16TJLru0Wi1DX0MMfmmKdTodlpeXmWeeOvWhvjp1\n5plneXn5pqWeypvBL02pg0s2F1kclHcucWlQ9llk8aalnlKklCY9hoGFhYW0ubk56WFIZ16326VS\nqQwt2VxiiQ022GWXeeZZZJE11gb9xWKRra0tv+CdQRHxQkpp4ajbe8QvTaFCoUCtVhtqW2ONueIc\nV65cYa44NxT6ALVazdAXYPBLU6vRaFCvX6/r91fvXLhwYWipJ0C9XvcELg14s3VpivXDvN1uDy3Z\n7C/1rFar1Go1Q19DrPFLM6Db7d6yjHNYu2bLmavxR8RjEfFyRLwaEU+Ne3/SrGjG9Z87OSzcDX3d\nyliDPyLuAf4Q+DjwKPDbEfHoOPcpSbq9cR/xfxh4NaX01ymlvwe+DDw+5n1Kkm5j3F/uvg/44YHn\nrwH/9OAGEbEMLAOcP39+zMORpkf97Hz9phkz8eWcKaVWSmkhpbRw7ty5SQ9HkmbeuIP/R8DDB54/\n1GuTJE3IuIP/28AjEfH+iPgZ4Ang2THvU5J0G2Ot8aeU3omIfwd8DbgH+GJK6aVx7lOSdHtjP3M3\npfRnwJ+Nez+SpKOZ+Je7kqTTZfBLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzBj8kpQZg1+SMmPw\nS1JmDH5JyozBL0mZMfglKTMGvyRlxuCXpMwY/JKUGYNfkjJj8EtSZgx+ScqMwS9JmTH4JSkzBr8k\nZcbgl6TMGPySlBmDX5IyY/BLUmYMfknKjMEvSZkx+CUpMwa/JGXG4JekzJwo+CPityLipYj4h4hY\nuKHvMxHxakS8HBG/erJhSpJG5V0nfP2LwG8C/+VgY0Q8CjwB/CJQBL4eEeWU0k9PuD9J0gmd6Ig/\npfS9lNLLt+h6HPhySuntlNLfAK8CHz7JviRJozGuGv/7gB8eeP5ar02SNGF3LPVExNeBn79F12dT\nSs+cdAARsQwsA5w/f/6kbydJuoM7Bn9K6aN38b4/Ah4+8PyhXtut3r8FtAAWFhbSXexLknQM4yr1\nPAs8ERHviYj3A48A3xrTviRJx3DS5Zy/ERGvAf8MuBoRXwNIKb0EfAX4LvDfgU+5okeSzoYTLedM\nKX0V+OohfZ8DPneS95ckjZ5n7kpSZgx+SZqQbrd7rPZRMfglaQIajQaVSoVOpzPU3ul0qFQqNBqN\nse3b4JekU9ZoNGg2m+zs7FCtVgfh3+l0qFar7Ozs0Gw2xxb+Br8knaJ+6Pf1w//q1auD0O8bV/gb\n/JJ0SrrdLu12e6htiSX2dva4ePEiezt7LLE01N9ut0de8zf4JemUFAoF1tfXKRaLwH7or7DCKquU\nKLHKKiusDMK/WCyyvr5OoVAY6TgMfkk6ReVyeRD+G2ywzTYlSlzmMiVKbLPNBhuD0C+XyyMfg8Ev\nSaesXC7TarXYZZcmzaG+Jk122aXVao0l9MHgl6RT1+l0WF5eZp556tSH+urUmWee5eXlm5Z6jorB\nL0mn6OCSzUUWB+WdS1walH0WWbxpqecoRUpn50rICwsLaXNzc9LDkKSx6Ha7VCqVoSWbSyyxwQa7\n7DLPPIssssbaoL9YLLK1tXXbL3gj4oWU0sKhG9zAI35JOiWFQoFarTbUtsYac8U5rly5wlxxbij0\nAWq1mqt6JGmaNRoN6vXrdf3+6p0LFy4MLfUEqNfrYzmB60SXZZYkHV8/zNvt9tCSzf5Sz2q1Sq1W\nG9slG6zxS9KEdLvdW5ZxDms/zHFr/B7xH9Oo/qEkza4rb31h8PjivZ88dLvDMmPcWWKN/xgmeRlV\nSRoVg/+IJn0ZVUkaFUs9R3DYZVRbrRbLy8s3XUa1/xpJebpdeecs8Ij/Ds7KZVQlaVQM/js4K5dR\nlaRRMfiP4CxcRlWSRsXgP6JJX0ZVkkbF4D+iSV9GVZJGxeA/grNwGVVJGhUv2XAH47qMqiSNipdl\nHrGzchlVSRoVg/8IzsJlVCVpVDxz94gmfRlVSRoVa/x9Edcf3+Yz8eqcks4aa/xjNqnLqErSqJwo\n+CPi9yPi+xGxFRFfjYj5A32fiYhXI+LliPjVkw9VkjQKJz3ifw74QEqpAnSAzwBExKPAE8AvAo8B\n/yki7jnhvsYrpes/kjTDThT8KaX/kVJ6p/f0z4GHeo8fB76cUno7pfQ3wKvAh0+yL0nSaIyyxv9v\ngP/We/w+4IcH+l7rtUmSJuyOyzkj4uvAz9+i67MppWd623wWeAf40nEHEBHLwDLA+fPnj/tySdIx\n3TH4U0ofvV1/RPxr4CLwK+n62tAfAQ8f2OyhXtut3r8FtGB/OeedhyxJOomTrup5DPg94NdTSnsH\nup4FnoiI90TE+4FHgG+dZF+SpNE46Zm7fwC8B3gu9k+A+vOU0r9NKb0UEV8Bvst+CehTKaWfnnBf\nkqQROFHwp5T+yW36Pgd87iTvL0kaPc/claTMnKlr9UTENeAHx3jJ/cDfjmk40yL3z8D55z1/8DO4\nH/jHKaVzR33BmQr+44qIzeNcmGgW5f4ZOP+85w9+Bnczf0s9kpQZg1+SMjPtwd+a9ADOgNw/A+ev\n3D+DY89/qmv8kqTjm/YjfknSMU1l8Od+A5iI+K2IeCki/iEiFm7om/n5w/7lQnpzfDUinpr0eE5D\nRHwxIt6MiBcPtN0XEc9FxCu93++d5BjHKSIejoj1iPhu77//lV57Fp9BRPxsRHwrIv6yN/9mr/3Y\n85/K4GeWbgBzd14EfhP4xsHGXObfm9MfAh8HHgV+uzf3WfdH7P+7HvQU8HxK6RHg+d7zWfUO8Lsp\npUeBjwCf6v275/IZvA38ckrpl4APAo9FxEe4i/lPZfDnfgOYlNL3Ukov36Iri/mzP6dXU0p/nVL6\ne+DL7M99pqWUvgH85Ibmx4Gne4+fBpZOdVCnKKX0ekrpf/ce/1/ge+zf5yOLzyDte6v39N29n8Rd\nzH8qg/8G3gDmulzmn8s8j+KBlNLrvcc/Bh6Y5GBOS0SUgA8B3ySjzyAi7omI7wBvAs+llO5q/ie9\nOufYjPsGMGfdUeYvHZRSShEx88v0IuJe4E+AT6eU/q53ZWBg9j+D3lWOP9j7XvOrEfGBG/qPNP8z\nG/zjvgHMWXen+R9iZuZ/B7nM8yjeiIgHU0qvR8SD7B8JzqyIeDf7of+llNKf9pqz+gwAUkq7EbHO\n/nc+x57/VJZ6vAHMoXKZ/7eBRyLi/RHxM+x/of3shMc0Kc8CT/YePwnM7F+DsX9o/1+B76WU/uOB\nriw+g4g411/BGBE/B3wM+D53M/+U0tT9sP+l5Q+B7/R+/vOBvs8CfwW8DHx80mMd0/x/g/269tvA\nG8DXcpp/b56/xv6Krr9iv/w18TGdwpz/GHgd+H+9f/9PAAX2V3K8AnwduG/S4xzj/P8F+19mbh34\nf//XcvkMgArwF735vwj8h177sefvmbuSlJmpLPVIku6ewS9JmTH4JSkzBr8kZcbgl6TMGPySlBmD\nX5IyY/BLUmb+P7NI+QJg5YZcAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<matplotlib.figure.Figure at 0x7f5005108278>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"plot_data(centroids+2, X, n_samples)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## course.fast.ai"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"source": [
|
|
"If you found this interesting, you might enjoy the 30+ hours of deep learning lessons at [course.fast.ai](course.fast.ai). There's also a very active forum of deep learning practitioners and learners at [forums.fast.ai](forums.fast.ai). Hope to see you there! :)"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"anaconda-cloud": {},
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.6.2"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|