{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In this post, we'll produce an animation of the k-means algorithm. The k-means algorithm is a very useful clustering tool. It allows you to cluster your data into a given number of categories." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The algorithm, as described in [Andrew Ng's Machine Learning class over at Coursera](https://www.coursera.org/learn/machine-learning/lecture/93VPG/k-means-algorithm) works as follows:\n", "\n", "- initialize $k$ cluster centroids\n", "- repeat the following:\n", " - for each point, compute which centroid is nearest to it\n", " - for each centroid, move its location to the mean location of the points assigned to it" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A word of caution before going on: in this post, we will write pure numpy based functions, based on the numpy array object. This has advantages but also disadvantages. In particular:\n", "\n", "- the code becomes efficient and fast, due to the fact that numpy supports vector operations that are coded in C\n", "- at the expense of being readable, which is usually what Python code is\n", "\n", "To follow along, a working knowledge of numpy is therefore necessary.\n", "\n", "To implement the algorithm, we will start by defining a dataset to work with. We choose a dataset containing three clusters, with a little bit of variance around each cluster center. We construct the point cloud by stacking shifted random numbers:" ] }, { "cell_type": "code", "execution_count": 101, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# necessary imports\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import seaborn as sns; sns.set()" ] }, { "cell_type": "code", "execution_count": 102, "metadata": { "collapsed": false }, "outputs": [], "source": [ "points = np.vstack(((np.random.randn(150, 2) * 0.75 + np.array([1, 0])),\n", " (np.random.randn(50, 2) * 0.25 + np.array([-0.5, 0.5])),\n", " (np.random.randn(50, 2) * 0.5 + np.array([-0.5, -0.5]))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the dataset as well as the cluster centers:" ] }, { "cell_type": "code", "execution_count": 103, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 103, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFXCAYAAACRLCZbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FPX9x/FXgOByBPGgoiISFUerVu3PArGFCIhiUQSp\nJSKBqGC1SFFUQGi1HljBqlWqoIIGQjViBSSCSDgMtg1Y61HPqZbgjUUQCcKSEPP7Y5LNtdnM7s7u\nzO6+n48Hjwe7Ozvzne9O5jPfO626uhoRERFxTyu3EyAiIpLqFIxFRERcpmAsIiLiMgVjERERlykY\ni4iIuEzBWERExGVtovmyYRg+oAQ4CGgLPG+a5i1OJExERCRVRFUyNk3TD/Q3TfMM4EdAf8MwfuZI\nykRERFJE1NXUpmnurflvW6A1sDPafYqIiKSSqKqpAQzDaAW8DhwPzDVN872oUyUiIpJC0pyaDtMw\njIOBl4Bppmm+7MhORUREUkDUJeNapml+axjGSuAs4OVg21RXV1enpaU5dUgRERGvsxX0ou1NfThw\nwDTNXYZhtAMGAbc3m6K0NLZvL4/mkCmhS5cM5ZNNyit7lE/2Ka/sUT7Z06VLhq3toi0ZHwksrGk3\nbgUUmKa5Lsp9ioiIpJSogrFpmm8DP3YoLSIiIilJM3CJiIi4TMFYRETEZQrGIiIiLlMwFhERcZmC\nsYiIiMsUjEVERFymYCwiIuIyBWMRERGXKRiLiIi4TMFYRETEZQrGIiIiLlMwFhERcZmCsYiIiMsU\njEVERFymYCwiIuIyBWMRERGXKRiLiIi4TMFYRETEZQrGIiIiLlMwFhERcZmCsYiIiMsUjEVERFym\nYCwiIuIyBWMRERGXKRiLiIi4TMFYRETEZQrGIiIiLlMwFhERcZmCsYiIiMsUjEVERFymYCwiIuIy\nBWMRERGXKRiLiIi4TMFYRETEZQrGIiIiLlMwFhERcVkbtxMgksr8fj+FhRsByMnph8/nczlFIuIG\nBWMRl/j9fkaOXEZp6RUALFv2JM88M1wBWSQFqZpaxCWFhRtrAnE6kE5paV6glCwiqUXBWERExGUK\nxiIuycnpR1bWk0AFUEFWVj45Of3cTpaIuEBtxiIu8fl8PPPMcAoLiwDIyVF7sUiqiioYG4ZxDLAI\n+AFQDTxmmuZDTiRMJBX4fD7y8s5zOxki4rJoq6krgRtM0zwF6ANMMAzj5OiTJSIikjqiCsamaW4z\nTfPNmv/vAd4HjnIiYSIiIqnCsTZjwzB6AGcCm53ap4jYpwlERBKXI8HYMIyOwF+BSTUlZBGJI00g\nIpLY0qqrq6PagWEY6cALwIumaf6phc2jO5iIBDVv3kquvfY8rAlEACqYO7eYa64Z4mayRATS7GwU\nbW/qNGAB8J6NQAzA9u3l0RwyJXTpkqF8skl5ZSkv9wd9rzZvlE/2Ka/sUT7Z06VLhq3tou1N/VNg\nNNDfMIw3av4NjnKfIhImTSAiktiiKhmbpvk3NIuXiOs0gYhIYtMMXCJJQhOIiCQulWpFRERcpmAs\nIiLiMgVjERERlykYi4iIuEzBWERExGUKxiIiIi7T0CYRkRjTIh7SEgVjEZEY0iIeYoeqqUVEYqiw\ncGNNIE4H0iktzQuUkkVqKRiLiIi4TMFYJMn4/X7y89eQn78Gv7/pak4SX1rEQ+xQm7FIEmmufRLs\nLeMmztMiHmKHSsYiSUTtk95Uu4hHXt55CsQSlIKxiIiIyxSMRZKI2idFEpPajEWSiNonRRKTgrFI\nkqltnxSRxKFqahEREZcpGIuIiLhMwVhERMRlCsYiIiIuUzAWERFxmXpTizRDa9CKSLwoGIsEoTVo\nRSSeVE0tEoTmeBaReFLJWEQ8Qc0CkspUMpakFunavrGc41nrDTdV2ywwZcpQpkwZysiRy5Q3klJU\nMpakFU27b6zmeFZbdHANmwWoaRYo0rSekjJUMpakFW27byzWoFVbtIgEo2AskmCSsZpbSz9Kqkur\nrq6O5/Gqt28vj+fxElKXLhkon+wJlVd1VcJ5AGRl5bteJRxtmhpXc2dl2avmDueacqsjlVc6cOnv\nzx7lkz1dumSk2dlOwdiDdJHb11JeeeUGX180acrPX8OUKUOpbVuFCmbPbrlt1e41FWmwTyb6+7NH\n+WSP3WCsDlyS1Ly4tq8X01RLHalE3KE2Y5EEorZVkeSkkrFIAonVkKtaOTn9WLbsyQZt2jk5wx3b\nv4gEp2AskmAiqeau7YENodupYx3sRSQ4BWORJOf3+7n00mcpKckFWp5oxMtt2iLJSm3GIkmusHBj\nTSDWRCMiXqWSsYh4jheHpInEkoKxJIxUvEE7cc45Of1YubKAkpLRgPc7ZWn+bklFCsaSEJq7QUOG\nuwmLIaeCks/nY/XqUcyZkxidsjTWWVJR1G3GhmE8YRjGV4ZhvO1EgkSCScUFFpw851gseiEiznGi\nA9eTwGAH9iMiktATmyTjIh4SH1EHY9M0XwG+cSAtIs1K5Bt0pMI552QKArVjnWfPLmL27KKEaS+u\nbVaYMmUoU6YMZeTIZQn/W0j8qM1YEkIqTkZh95yTscNTIo51Vlu3REPjjCVhpGK7p8/nC5SGCws3\nBi1ppWJ7ukiyiXvJuEuX5O396iTlU8v8fj/z5q0EIC9vYFIG6MazZ61cWcDq1aManGtGRtPzzsjw\nNbmGdE3ZF0leTZw4pMEQsuzsxUycOCopr8tauqac48h6xoZh9ACKTNM8rYVNtZ6xDVontGVeWXc3\n1mOf7axfXJcXeYA1jrhxXuiasi+avEqlsfC6puyJ23rGhmE8DWQDhxmG8Slwq2maT0a7X4k9r904\nwkmPF9rnommrdTLvU7E93asSsa1bvCHqYGya5mVOJETiy2udfryWHjsifSAI51ztLmkYjyDgtYe3\nxryevmgl+/mlOnXgSlFe6/QTbnoSeahTOOfqlWE+Xh+24/X01RfJMLREOj+JjIKxJKTaIDV3brFr\nQSpeDwRe6EUer4e3SMdLe+3hsjmRBtVEOT+JnIJxErJzQ/NKybI2rZWVFfTpsyCs9Ph8Pq65Zojj\nQcpuQIi01OqVvPeaVCj9KahKczTpR5Kx2x7phU4/jdPau/djzJy5lPT09JimJ1TbW7ht15G01Xoh\n78Nlt+06GtF0yotH+tyU7OcnCsZJJ5wbmlOdfloKbs191jitmzePZ8SI2PaIbinYxquXdqL1uvX6\nA4TX01cr0qCaKOcnkVMwlibC6bUZKri1FPgqKyub7C/Ye07ywpCoRBXrB4hoS3+J8IATTVBNhPOT\nyKnNOMlE2x4ZbrtdqDawltvHqoGFgbTCopr3nBNuhyC157rHKz3HY80LHfLEe1QyTjLRVmfFs+SY\nnt4WuBAornnnl6SnF4f4RniClcwXLrwgZOlL1YHuClb60/haZygfvU3BOAnFszorVNViS9WO1udP\nx6xTSrAHi+XLi1oMtolcHeiFG66TaYjXZDD10zxx4hBH9+0FiTipTqpRME5hwW6a4bbbhSpJtlTK\ndKsUmsjBNpja37GyspIVK3awefPVgDs3XKdv+vGoqWmc5pUrCygouCipApX6SnifgnEKqR98hw3r\nxdixLwa9adYGyMrKCqAThYUbQ5ZwQgW3lgJfLANjKgwHaRxIrDb4KsAX1vSckZRkg30vEW/6jdNc\nUjLa82mW5KNg7AGNb2rg/LJkjW/ac+feT1nZr6l/0ywoWMr48UMCJWQ7JRwvVIs2JxXafxsHEhiD\n1QZvr6o10pJsc99zWio8UMWD8tH7HFlCMQxaQrGRYEsBrl8/hvJyZ4f4BFuKD14ELg68zsy8lZKS\nG/H5fGEu3Rf7ZQybC/qpvoxb8N91FfDzBkspNpdPdn5nu8edPbuo3kNcHhB8OcdwxfqBr/ESlNnZ\ni5Oumhqcz8dU/9uzK25LKEp0glXr5ecXM2KEc8Np/H4/paXvAUMbvN++/XPs3dsf8AGLKSubTmFh\ncaMbcQVQDuziu+92s3//fg466KBm0x6L6r1k7XzixM2xcYmnT58nuOiig0lPL3KlJiAWtRGxbuNv\nnOaJE0c5/jDsBcnWVyLZKBgnubpAdj1We+KYmk8WsXfvw8AjwInAKKAV1dXVvPvuO+zc+SadOk1n\n9+4twAEAbr8d7rgjjTPP/DH9+vXnjTe+BdKAn2MF9NhIxHbIljj1gNE0+I0Iax+RVl+G+l4i3vTr\np9nn8wUNxl5ukpHEp2pqlzWuIsvKyne0mrphdaIfuA84FTgfK4DWVWueeuptpKW9yNtvvxXGEY7H\nap+8k6yswpiUWENVpTpVVRbvG22k1cORCpVPu3btYurUfABmzcqjc+fOtvaZrMEpWF7Fs0kmUaia\n2h5VUyeI5qr1nArGVo/oVVg/9UDgFKySbHpgm2HD3mLv3kKKi5cQ7OGsc+fO+Hzt2L/fz65duxpt\n81/gITIzS1m0aKmtm5Pf76egYC2vvfYRZ53Vk9zcgSG/F+vOJ16tBo8m2NntFOj3+2t61U8DYNs2\n++ceTQk40QJ5MtbOiLcoGHtArKr1/H4/RUXfAr+oeWchP/nJHlq3XsCmTVcCVmA74ojtPProMw3S\nc+6559O+fTdOPPEMrr764sDN8ptvdnL77Q/y1FM7gOeBbwEoK/snkydPZMGCRS2m6dJLl7B580HA\nNJYtg6KiBSxZ0nz1aqx7Rbtxo23pASOaB4Rg312/fkzQbd0490jOLdGCt0i4NDd1Eiss3MimTVdR\nOzc0jGHYsK4sWTIiMP/vHXecxPz58wLfOfvsn1FcvJEdO37OkiUPcNddYxvMT33IIYdyzz23kJX1\nU+Ad4JLAd4uKlrNu3ZoW07R5czdgbCBdmzZd2eKarsk2n29L8zBHs+5tsO/m56+LxWlEJNxz88I6\nx5qzXGJNwTjFpKenBwLb2LGDuPXW6VRVVQFWIP7rX1dQWvpJyJtlXSDZzKxZlzF8+C8Cn1133XXs\n2rUrorSFu6iDU9y60XrhASMRgkw0DyZOSZVFLMQ9CsZJrKUb7WeffUpp6d8Dr++55z7atLHXclEb\nSK644nx+97vbadXKGu60Y8c2hg37U7PBdNiwXmRkvE791Zp69ZrPsGG9QpZ+YhmovXijjSZIBvtu\nXt7AoNu6ce6J8AAQjBceniR5qTe1h9S2i2Vk+BgypJcjf/D15y2GatLT2zJsWC+WL3+VTz75D3/+\n83QAfvjDU3n55X8EvtPSxA312/AqKyuYMeNpoKjm08XMnh28HdzqRXwesBYwgROYObOS9PS2zfYu\nDtWTNZl7dEbbgat+J7nJk4e5NnY22HnUXZcVQBrp6enNnqOd69FJyXxNOUn5ZE9K96ZOtM4etTfO\n+fM/oaxsMuDc0IlgU1veddcsdu+eBGQCVjDevv1/VFVV0bp16xY7TDUOjpmZ9wOf1ztqO0KvS+wD\najsrVfDaa7M466wTmmxVWvoeOTn9UrYna7Qd+154oZzSUquT3Jo17ix+EKqzlt0pVxN5WtNEuxeJ\ne5KumtoLnT3CUZveGTPa1wRi59vFGre57d49BdiINUbYujls3/4/8vMXBL4Tqkqu8f7Kygzg9cDn\nZ5zxTrPVjjk5/ejTZwG1VZSwiGXLrqeo6Ft6936syfsjRy6rKdXX8gOrKC19z9O/q9sa/0bW4gfx\nbWcNlo7613U4bcGJWEWcaPcicVfSBWMvdPYIR1163aikOBiYFnh19913sGlTaYvfahgc3wN+E3j1\nk58MZMWKa0PeMKurq4B7sebGHg1ksGnTlRxxxHZOP/16rOpu632rarK6po1xN7AYGMKyZdM8d3Nz\nqwOa2BPv3yfR7kXirqQLxolrIFBAtJ1agt1wGneY6dRpNtAXqKB378Pp0SMTgPLy3Qwf/nNuu20G\nn3/+WbP7X7FiB/Bn4C7gJ8CngDU5yMKFj4cMxNbQpquBM2g8jeaKFWfy1lsPAh9hBWor/enpbXnm\nmeEMH/4g9YdEeWnIjtdKQY1/8+zsxa50kgrVWSueHbm89vuINJZ0Hbji3dkjWg3T6ycz825uuKEP\nw4b9NOw0h+roVL/DTGXlAd56aytnnXUCubnn8vbbb5GbO5KdO3c22N9JJ53MgAGDOPFEA5/Px/79\n+1mxYj3r1pVhVUvXXTvt2rVj3rwnuOCC0Ev31U0DWQU8hVUCBlhU8//aKTpfBLbTu/d+nn32l82u\nJDV3rrOLakQq3tNbhuL3+/n44618/PFW1q59naqqAwwefCZVVa1p2zadQw45lB49MjnkkEPjlp7m\n2k3j1aYazu/j5BSriXQvCpc6cNmTsh24Eq2zR9P03sgxx3SxfZE37tXcXEenph1m/Lz55kwgjdzc\ngaxd+woTJlzdYKjTBx+8zwcfvN9iGrp27c7TTxdyyimntrhtw5mnfkFm5q2cdpqPFStupOFiE22A\nMQwdWjfFZrBZq/LyxrB9e3lKdpKprq7mww//w4YNa3n//ffYurWMrVvL+PLLL5pMa1pQ0PT7Bx/c\nmR49MunRI5Pjjz+Bvn2z+clPetO2bVtH0xmqI1oiLiphV6Ldi8RdSVcyTgZ2nziD9WouK/s1dfMQ\nN3z6b1oqzQXqStDV1dWce+50PvzwY6wOXhUhjt4a6ENmZnfWrbuPjh072j6/pvMm06AEYbULW6tI\nBVtDuf53u3TJYMCARQ3yYNy47uTmnhvVXM6R10pY5xCrUtCePXv42982sm5dMRs2rOWTTz52dP8d\nO2bQt282AwcOYuDAQRx9dDdH9++WcH4flfjsUT7ZY7dkrGDsQXYv8mBVb5mZt1JWdgfQ8Ibj9/uZ\nNGkey5ZNA9YA59G4yg6ot7/vgGL69XuCo47qxL59+/D5fHTq1IkzzzyL7dvTadeuQ8ixoeEEN2t4\n1zrmz3+TsrLpgM9WQHvuuY1ce23Dc4EXycr6KuK5nCMdVhbLKtdPP/2ERx99mMWLF7F373cht23V\nqhXduh3D0Ud345tv9tG6dRuOP74rfn8F+/fvZ9u2bWzZ8t+aMb6hDRw4iAkTJvHTn/YlLc3WPcWz\n7P4+CjL2KJ/sSdlq6lQ3btwZpKc3rBZruqbx4Tb21AH4ORdeWBV2NWIkCwH4fD7Gjx9Cbu5ACguL\nG6Q/fG3CGovs1DjmWFS5vv32v3n44Qd5/vmlgWlL6+vQoSP9+p3Dz37Wl+OPP4EePTLp1q0733//\nPSNHLuODD6zf4PDDrXHGYNVCmGYe8CWnnPJnrryyO//+91usX1/Mp59+0mD/69YVs25dMaeffibX\nXTeJIUOG2p6lzWuSuUpcEl9i/lUJELwNNTe3aQArKFhLaekRWFXPvwDWcNhhd7Jjx28D36tdMciJ\npQqt4HYZVgkcSktzKCwstnUjtHPDrF/CGT9+IIsWPRmkijuxbdv2JdOm3cSqVUVNPjvhhJ4MHjyE\ngQMHNdvGm5+/psEDhjXO2NpX3fvH8u67M6mqKuLee/Oorq7mo48+ZN26NaxdW8wrr7wcaHt+6603\nGD8+j+OOO5577/0Tfftmx+rURVKSgnECs9NBxO/3M3/+J8DkmncKgF8weXJVkxI04EiHE6v68xms\nYUgAC6ms7BD2foJpXOpeubKAhQsv4Nlnl9ar4m5Fnz5PUFl5MPn5a1qsMo71esnhqK6u5umnF/O7\n391CefnuBp/17ZvNhAm/oX//c2NSZZyWlkbPnifSs+eJXHPNdWzZ8hFz5z7MM8/8JTAMaMuW/zJi\nxEWMHj2WO+64m44dg6+TLCLhUZuxBznZFtNcu3JJyY0x69n5+OMrmTFjRINjzpy5lPHjQw97aixY\nG1+oISr15+FesWJHzXhme23AXpi2cPfub7nppkksX760wfsXX3wJ1103idNPP9P2vhp3WMrOXtyg\nmjqcjmbbt2/niSce4/HH57F797eB94877ngefzyf00473Xa6EoHaQu1RPtmjNmNp1rhxZ8Q02KSn\np9t6r77me1g3bHcOpbaKOz9/TU0gbtgGXDvPde0x6udB4+rxr77axujRDwCwePENHHFE15DHjtbO\nnTu45JKLeO+9dwLvZWYexwMP/Jmzz/5Z2PtrXGsyfvwlPP64NUHKwoUXsHx5Xe0HWA9t1uumDyJd\nunRh6tQZ5OVdxbRpN7Fy5QrAKiUPGTKIv/zlWVVbi0RJJWMPcvKJ042JB8I9ZrDezBde2KlJ6Xr2\n7KJ6Y6WtfdeW+OrvO1jpeebMpbzwwm5bPaa/+mobZ5xRQFXV7wBo3fpO3nwzN2YB+dtvdzFixFD+\n/e83A+/l5uZx55330L59+6j37/f7yc0toqSk4VC2hp377Pckf/bZQqZMmcx33+0BoH379hQWLqNP\nn6zA8dyuZYiGSnz2KJ/ssVsy1nSYSaj+lJhA3NerDXeN3GBz+L722ofNbn/hhZ0YPnwWM2c+x+rV\no5rsO9g0i1Bte57g0aMfqAnE1rZVVb8NlJKdtmdPOTk5lwQCcVpaGg8++Aj33feQI4EYrPy1AnF0\nizXUuvTSHIqLSzjyyKMA2Lt3L6NG/YJ//eufmnZSJEKqpk4yzQ0rSrQhHWeddQLbtjXsVDVs2AUN\nzm3btieZPLnpd4N1bPPiBP3V1dVcccVo/vWv1wLv3X//HC67bHSIb3nDCSf0ZOnSIoYOvYDt2//H\nnj3ljBx5Cb/+9R8oLR1Hqi13KRItlYwTUKjVZ7ywUky4paNgJdnc3HOblK6XL3+1ybk1t1BE4yX3\nwlmUYPHiG2jd+s7Atq1b38XixTdElSfB/PWvz1BSsiHw+g9/+COXXz4m6LbRrDiUk9OP7Ozgi5AE\ny5dhw3rZOtbxx/fkueeKOOywwwCrA1pR0ZNhpU1ELCoZJ5hIJtSIt3An0Qg1RMupElU48wQfcURX\n3nwzl9GjrwdCd+CKtH109+5v+f3vfxt4/atfTeCqq64Ouu2uXbs4//yHKSvLAgawbNnTYf3mPp+P\n1atHMWdO03NvnC/Dhl3A2LEv2r6+TjrpZPLzn+aii6zf6d13/8lJJ03mgw/uA9wdJiaSSBSME0xL\ngc5LY2ad1txCEeXllaG/WCOcGZiOOKIrxcWzQm4TzYPR7Nl3s337/wDo2vVIpk6d3uwxzj//KcrK\n7qx5pyCsSVRq2V2sofFkIXaqmXv37sOoUbk89ZS1GsXevUXMnHk26elttTiCiE0KxknGCyvFhPtA\nYDeoNXdudoOx0yKdRnPnzh0sWPBY4HWoyTMKCzdSVjaZup7ho4FV0SfeYb/97e2sWlXErl27+OST\nj2nffmezVe4i0pTajBOMnbbPxu2lsRCqDdNub+rafUyaNK8mqFUBaygt/QEFBWuDHjce5xY+P7AS\nWGVr8YUNG9YF5pk+/fQzufjiS5ru0e/n8cdf4C9/WQesqDmGJTNzU7Pt3dEKp229/jXQsWNHrrvu\n+sBnxcUvxSR9Iskq6pKxYRiDgT9hrak33zTN0HV7EhUvlHztlGRbqhJuuI82WMHmWWqXdZw//35y\nc/0xHw8dzXjYnJx+PPfcY2zefBC1U38WFS1oMd1r164J/P+CC4Y0mdrS7/dz6aVLavb7YM27TwI5\nZGbO5aWXJsQsX+xeX8GugZkzB3HXXb8HoKRkAxUVFY6vjSySrKIqGRuG0Rr4MzAY+CFwmWEYJzuR\nMGleOKXDaHrhNseJHtsN93E+cDdWILb2WVZ2Q0x7gTsxHtbn8zF06GFYgdhK96ZNV4ZMd1VVFRs2\n1JX6zz236QNLYeFGNm/u1mC/kMfw4Q9SUnIlnTt3Diud4Wp8fQW7hoJdA//85xd0794DgO++28Om\nTf+IaTpFkkm01dS9gI9M09xqmmYlUAhcHH2yxAmJMwGDD/i/kFs4/VDh1BCwlqb5bOzjj8vYuXMn\nAIcffjinnvoj29/Nyvqha7Ugdq6htLQ0+vcfGHj9+uuvBd1ORJqKNhgfDXxa7/VnNe+JB8RqzHE4\n7Yp299G791f06bMg6D69/FARbl7s2bMn8P+uXY+iVaumf4I5Of3o3fszrLWnrf326fNEzNqJQ2nu\nGmruvI866qjAd7/77ru4p1ckUUXbZhz2xNZdumjJNTucyKeMjKalqIwMX5N9W6VOa/KMvLyBNkpf\nGaxfP4b8/OKaDks/YOXKV21+t+k+rONeBVDv9ZjAvubNa9prefnyVYFSaXjHtUycOISVKwsoKbFm\nu8rOXszEiU2n1gz/PMaE3Ee7dq0D/+/QoV0zv3MGL798FY8//iL/+Mf9nH32SYwfnxd1qTiSa6q5\na+iYY7oEPe9DD+0U2K5Vq+8T9u89UdMdb8on50QbjD8Hjqn3+his0nGzNLF4y5yagH3IkF5kZTUc\nYjRkyPAG+27cEWfRIvtjZYcM6RXxd2uNGGGV9mqHJ9V/XfteeXnTUvADD2wKjL2N5LgABQUX1euo\ndFGDY4YrWLqD2b+/7v+7d+8J+Tvn5JxLTs65tvbbkoyMdObMWVmzX/ud1Zq7hj79dHuDzm+16du+\n/ZvAd9PS2iTE33vjjnzHHNMlIdLtNi0UYY/dB5Zog/FrQE/DMHoAXwAjgcui3Kc4xE7P2EjHykb7\n3XA0HrecmfkAZWUzoj5uOJOAOKV26kiArVu3sG/fPtq1axfTY1q9s58NrNoUzuQkwa6huolIJjfZ\n3wcfvB/47qGHHhZ0n14SrFf4+vUaHy3xF1WbsWmaB4DrgJeA94BnTNN8P/S3JJ68OS43PI3HLY8b\n1x2r01fiOfLIozjhhJ4A7Nu3j9LSv8X8mKFWbbKj/jUE1EzNWTsRSd3+Kisrefnl9YHv9e17jpOn\nERPB2sSbm+/cbbEYGSHeEfU4Y9M0XwRedCAt4oJops+M59Sb9Uuxfr+fF16I7XFjuSbvwIHn8dFH\n1hKRxcUvcfbZfRNm/V9rRrCsoJ/985+b2b37WwCOProbJ52kUY5OSYQ56SU6adXVYffBika12hha\nFu+2mJYCT6jP3VpIvva4GRk+hgzp5ehxG9/4srKcvfFt3Pgyv/jFUAB+8IMj6NFjBq++ek1MjgXW\n+eTmFgU6q2Vl5Ud8jPz8NUyZch7WBC3W/jIzH6Ck5EqmT7+ZxYsXAjB27FXce6/za0A7fb3V/dZ5\ngJU369cy7D3/AAAek0lEQVTbn+88Xqx8H0rdtKgVzJ7t7tKUajO2p0uXjLSWt1Iw9iQvXeSxDkx2\njh/q5huLvIr2xtdSmisqKjj9dIMdO3bUvDMVuCfksewEodptrB7uaaSnpwe2jbQDV7BjWNdDDrCe\nzMxNvPTSBLZu3cL55/en9n6yZMlyzjlnQETHaPnYzl6LidCBS8E4cdkNxlooQkKKVyetYBKxas5O\nmtu2bcvvfncH118/oead+4ErgRMb7Kc2QAwb1ivosoZAkG0uA56hdnrO2m27dMlw5Der69BlDWnK\nybmRtm3bMnXq5EAgHjhwENnZ/aM+VmOxuhbd6MgXrmRejU0sCsYx5FYVbrJw60GgpRtfqN/Vbppz\nci6noCCff/3rn0AlMA5YSVbW0wwbdgG//OVzbNpkjb1+5JE/snXrdQ32+cQThSxc+GWgR/PcufdT\nVvZrYCN102jWHX/ixCHk568JmuZwNQ5ec+f+mTfeeB2wHjRmzpzdZL5tCU+wa8ztOeklthSMHVa/\nqrCo6NvADdULpbrGf+BAiw8LqfhEHurG51RpvVWrVsyefT+DBmXz/fffA69wyinnsnjxMgoLN9Zc\nN1ZA3br1Rqw+knUzzc6ZU8KOHfMC25SV3VCzTdM/6crKCgYPfiqioU0teeqpAm67rW4t5uuuu57j\njjve1nfDfVhNlWsx1DXm9RK8RE5txg5q+Ee0Cvg5kbTxxKItpvEfeJ8+C6iurmLz5quB0O1vTpXw\nw91PsM41jdMY73arltru7KS5vvvum8WsWTMDr63lFE/j+efrxlFDBYcddhM7dvwRgGOPnc3HH1cD\nP8JaZMMHVJCZeStlZdOBJcCYwPEvvLATM2aMaDbNkfrrX59hwoSrA9XTP/7x/7Fs2Spb46Yjbf+N\nR22T222hXmwfDsbtfEoUajN2QcMqSm9lbePq002brsQqSQWvTt23bx8ff7yVTz7Zyt69e+nQoZI2\nbdqwYcM6unU7hszMTDp2tD8VXiQlSp/Px8KFFzB1qrUq56xZ0U8JGWvhVidOnjyFPXv28PDD1lKJ\nzz+/lFNP/Qw4HLi6ZqtFTJz4f7RvX0RlZQWPP94auKnms4XALwNLKy5fXkxlZQdgaU0HruGOr35V\nXV3NokVPNmgnPu200yksXGp7ApNImyBUOpRk5a2IkVQGYt0o60oo3q5S+4i//30Va9Y8xjvvvM22\nbV+2+I3DDz+ck08+hXPOGcjAgYM4+eQfNttWGMnN1+/313RKmgbAtm3uV/XbqSoNJ2CkpaVx6613\nUFGxn8cfnwfAO++8Sps2/+HAgXeB/vTuvZ8rr7wYn89Hfv4atm69ibpS0xgOO+xaXnrp93Tu3Dno\ncXNy+jWYhzuaa/Gbb3YyZcpknn9+aeC9k08+hSVLltO58yG291NZ2XToULD3UlGqVMdLQwrGDmr4\nR9SK3r33M3RoXQklllXA4aUN+vR5ggMH9vLaa/cBjwIf8vzz4e3z66+/5pVXSnjllRLuvPNWjjrq\naEaPHsuVV453ZCpEpztwOZHXsehIk5aWxl13zaK6upr58x8F4MCBXcBj+HybeeihJ1soXQ8Kucax\nz+dj9epRzJkTeZqrq6tZsWIZt9xyM19/vT3w/sknn8Kzzz7fYJpPm3uk/sMqLGDz5q8aDMdKVeqs\nlZrUZuywcG74zbWbxWqcY23a9u37jq+/fpNFi57k2293Bd22devWHHNMd449tgedOx9C69at+f77\nKsrLy2uqrz+moqIi6HfbtWvHqFG5TJp0I127HtnoXPNqzrXliSeaazvLyekXyOOJE4fYmqCh7viX\nYY2PLeWllyaEDGJu+M1vfkth4RNA3VKLrVq1Z+LEa/nVrybQsWPHsPMRIm/fq6ys5Pnnl/LII3N4\n551/N/hs7NiruOOOuyOaW7tu8pCNWMswbgfcGcvemNpC7VE+2aNJPxJAc8Hm5ptHxOwif+mlF7np\npkl89dW2Bu+3a9eOn/2sHwMGDKJfv3Po0SMzsERhMFVVVXz22af84x9/Y926Yl5+eX1gKsRaGRmd\nuP32mVx++RjS0tIc6cC1cOEFDcbcZmcXUFBwUYv7ajhzlNWrODPzfkpKrvRUqcNK59dYY4U3NPjM\n5/ORk3M5eXnjePXVzwH7Jfxwb5zl5bt56qkCHn30ET777NMGn3XteiSzZz/A4ME/t72/xhr+tquA\nIXilw1KyBxmnauOSPZ+cog5c0sCBAwe49dZbAtWgtXr0yOSaa64jJ+dy2rdvb3t/rVu35thje3Ds\nsT247LLRVFZWsmLFMh5++KFACaq8fDeTJ09kw4Z1PPTQXDp06BDWDTZYdV3jquuSktFNqq6bv9ms\nxwrEdcOB4jWBiV05Of2YO/c+yspWA08BtwNbgdqFAhaQn7+Anj1PZMCAQWza1IasrJ9y0EEHRXXc\n6upq3nvvXdavX8v69cVs3lzKgQMHGmzTrl07xoy5gptvvoVOnQ6O6nj1f9vS0vdYtmxIVPsTexJx\nIp1UoZKxi5qrunW6mrqqqorf/OZann22sN67XTnhhCEUF99Fhw4dHDtWdXU1Gzas5ZZbbqasbEvg\n/b59s1m8eEnUywXaH1rUsMoTIDv7XsrKZjb47syZSxk/3luBYNeuXTVLFN4AHKBnz0kcdNCrTaqJ\na7Vv357/+79eHH/88fTocRw9emTSo0cm3bt3p1279rRu3ZouXTL46qtvqaio4KuvtrF1axlbt5ZR\nVraFrVvLeOONf/Hll18E3f9hhx3GuHHXcMUV42KyLGIkTRixlMwlPieHTSVzPjlJ1dQJIlgpzsmL\n/Pvvv+emmyYFJvC3DAfmAx1jVh24b98+brttOvn5CwLvDRhwLgsXPh1VKa7xjTs7e3GDaupQN5tH\nHlnK73+/Dbih5rNFzJzZgfHjLww7DbHudNf4GAcddBB/+9tGHnvsEUpKNoS1hF6rVq1o1apVk5Ju\nS0499UeMGXMFI0eOsv0Q5ff7KShYx2uvfchZZ51Abu65tvLHS7PVJXOQUTCOP1VTJ4hYj5t8+unF\njQLxOOAxIA2r40xstGvXjtmzH+Doo7sxc+btAKxfv5YHHpjNtGm/a/H7zd2cG1ddT5w4yvYKO+3b\ndwR+DRTXvPNL0tOLQ3wjeLriUc0X7Lro2zebvn2z2bdvH//4xyusW1fMunXFDWoggvn+++9rZvkK\nrV27DgwceC6DBp3PgAHncsQRXcNKs9/vbzCN57JlC1mxYgnPPvvLFvNH44fjQ8OmvEslYw9y6onz\nm292kpX1Y3bu3AnAJZdcyhdfDKiZ8COy6sBISjD33HMn999/L2DNXbxx4yaOO+6EkMewOztT47wK\nVeXpRHWoF2dHKivbwgcfvF9T9bwlUP385ZdfUFFRQf2/8fT0dA455NBAVXa3bsewYsXXfPTROOAM\nsrIKolxesWHewCpmz65KqECb7CU+deCKL5WMhT/+8Z5AIO7W7Rjuv38OrVq1inj8YqSlwilTZlBS\nsoF//es1KioquO22GRQUPNPs9tGMLw41RjNZx29mZh5HZuZxQT+rrq6mqqqKww7rwK5d/iaTsuTn\nr+Gjj+oCaGlpDpMmPUhW1g9dry6W2FAthDe1cjsB4jy/388TT6zm6afrOmzdcccfaN++feAPMS/v\nvLBvtA2DZHpNkGx5qsVWrVoxa9b9gddr165pMgzKSaHOMZrzB6skkZX1JFapr6Kmmq+fMwmPgbS0\nNNq0aUPbtm1trKTkB5awbNk0pkwZysiRy8Jqm87J6UefPguozRtYRO/en3s6f+ywerGvIT9/TVj5\nIRIOBeMkU1t6nTatO3v2WBN6HHroYVxwgbs9hn/0ozM4/fQzAat3d0nJhma3jSbgxfrGWVu6nj27\niNmzixJ+WEjDvF5N3fKL9h+2avl8PpYsGcHMmUsZPnwWM2d2sNVe7GW1f09TpgyN6AFFxC5VUyeZ\nutLrQ4H3evQ4hdatW0e972g7fwwcOIi33noDgL///RUuumhY0O0irU52s3NVIqptOzz/fB9ffDGV\njz8ux5p8I3I+n4/x44cwfrwzaXSbW2tqS+pRME5a3wT+d/jhRzqyx2jbXI8/vq7T1q5d34TYMrKA\n58aN00tDcsLR+MEF3gfuAAoAa0GJzMwHyMm5MuQ+EvHcRbxIwTjJ1JVe64aynHJKd8f2H02psP7w\nmrS08FpIam/8lZUVQFpgQQGwv4yj0xJlNiO/38+8eRspL/cHgmbjBxfojbUu8iisoV+VjBvXvdlz\nSZRzj5aGAkm8qM04DInQkaO29Hrxxf8LvPef/5gupqjOu+++E/j/4Yd3sf29una785gx4ztmzBgR\ntP0u3p2rCgrWRdShLZ5q8+7aa89roc1zAHA/1i1hEH36bCM399xm9xtpZ754cPLvNNn6CIh3qWRs\nUyKVBHw+H1OmXMvzz1uzX5WUbKCiooK2bdu6mq5169YE/n/OOf1tf6/uxr+Gug5GVjV0fn4xI0ZY\nATeeQ5f8fj/z578JjIjJ/p3SXNV90yU1n2Lw4K689dYszjqrJ7m5I2zknR/rNwHoG6MzsKeu5qSS\nFSt2sHnz1YAzf6fJ0kdAvE3B2KZE68hxwgk96d69B598spU9e8p56qkC8vKuci0969ev5cMP/wNY\nN7ezz47NzTteN87Cwo2Ulc0gnDZWL2n64GIn+NYZNqwXd931ILt3TwWgU6dZDBs2NiZpbUnT9u+F\nQBXgc/3vVO3qYpeqqZNUWloaOTmjAq/vvvt2duzY4Upa9u/fz/TpNwdeDxs2IqwFI+qqn/ti3Wjr\nqqHz8gY6ndww1G9jXRWyjdUtoaruoxlzvXz5qzWB2Kqm3r17CsuXv+p08m1pXGUOY4B1rqSlPg2L\nknAoGNvkVHtkPNudJ0yYRPfuPQBrJaBx48awd+/emB0v2LlVVVVx/fUT2LLlvwB06nQwv/3t7WHt\nt67drpiZMzswc+ZS19vv6q4Hq401K+t/IdtY7YjFtVGbd3PnFrueZ/FViduTsni5XV28R9XUNjnR\nHhnvdud27drxhz/M5vLLfwlYY3vz8kaxaFGh48cMdm5PP30xt9xyE889tySw3S23/JYf/OAHYe/f\ni+12F16YQdeutW2s0f2Osbw2fD4f11wzxNF5hOPRy9huFW/T9u8nuOiig0lPL0qaKU8l+WmhiDiy\nu8iA0xOwz5nzJ+6889bA69NOO51586wF6p3S9Ny20rPnxXz4Yd0avHl5VzFr1v02pmW0x+/3s3Ll\nqw2G7ES6n3Da9cJZyMKuWC9AEYtJ/VvKt2jaS8PNYyfbZp3KK6+t0+w0LRRhjxaKSDGhbkYTJ17P\n/v1+Zs++G4C3336LAQN+ysSJNzB+/DUccsihDqbkO6y1km/jww/r5p++7LLR3HPPfU0CcaQ3UadK\nko33s3TpgppSVdtm05NonfliJVRtRbS/T7h57MWak2RdmERiQ23GcRSrcbB2OorceONU/vCHewPD\nm/bv388f/3gPZ555CjNmTGHLlo+iSsM55/SkW7fhwLHA9YAViNPS0rj++psCK0aFm+7mONUe13g/\nmzZdyYwZ7ePe4SbRFqBoidpLLdEuTCKpQ8E4jmI1gYCdG19aWhpXXfUrXnrp5cCCDQB7937H44/P\no0+fH9O79xlMn34zq1ev4qOPPmT//v1Bj3fgwAHKyrawfn0xd955G+eccza9ep3OZ5+tBOp6bGdm\nHsfSpS8wffqtQefG9u4Nuw2h0hOLwKnJJRqK18NJIkzkI6lB1dRx5nZ12imnnMqLL67jueeW8Mgj\nD/H+++8FPisr28L8+Y8yf/6jgBXAjz66G506HUx6ejoHDhygvLyczz//lKqqqmaP0b37sVxzzQQu\nv3xsWEOYwuFUB6LG+4FF1I4bbk6sqh/dvjZqOdH+Gu3vE48q3uaq0t2cYlVSlzpweVC4HSMi7ShS\nXV3Nhg1refLJ+Wzc+DL79u2LOM3p6en06XM2ubl5XHjhxbRp0/JzXrB0L1x4QWC8akuBwOkOXHWz\nN40PpCdZSqh2ryknO6cFC+p2A308JstortPczTePUMckG9SByx67HbgUjD2k9gaUkeFjyJBeYfc+\njebm5ff7KS39O+vXr+Xdd99m69YyPv/8M5q7Po44ois9emRy0kk/ZMCAc+nbtx8dO4Zfoqif7mHD\nejF27IthBYKWbgiR9JROxhmT7N44Y9mr226gj0Vv9WAUjKOjYGyPelMnmGhvQNFWcfp8Pvr3H0j/\n/nUzWu3fv5/PP/+U777bS2VlBenp6fh87Tj66G60b98+4mM1Pm5tuvPz1zjaSzmSHr3R5mOyBnMn\n2O0hHa/e6lqRSbxEwdgjvDhc5qCDDuK4405oeUMX1XbAgabBz8k8tRNkE2kxkeakUoDS0CPxEgVj\nCSme1bzhBgK/38+llz5LSUkuELvgZzfIhlohKVFKy7EMUHZ/33g+EHil05yI2ow9wouz9UQyC1K0\nbX3hBPOW2jedylO77ajBtps58zleeKE85u2fLfFK+56XOnA1xyt55XXKJ3vUZpxg6pdIrA5c7leZ\nhVvN60S1sJMllXhXQwYr0UEnzzU/uMnu76sSq6QaTfrhIbU3oGuuGeJ6IE4EOTn9yM4uINTEEE7M\ngGR3AopgE3ekp6fX28IPrKK09D1NMCEiDaia2oO8Uv0TbjWvG1XtGRnpzJmzEohtdWb0c2jnAEuA\nsUD8q6ujvaZSqZe4V/7+vE75ZE/MxxkbhnEp8HvgJOAnpmm+buNrCsY2xCvA2BHuTXjXrl1MnZoP\nwKxZeXTu3Dmm6UuEG4Lf72fSpHksWzaNWK3K1JJo8ile4369IhGuKS9QPtkTjzbjt4HhwKNR7EMa\niVcPYbvCabvz+/01k3ZMA2DbtuS+advl8/nIyvohy5a5nZLIeHHYnUiyibjN2DTND0zT/I+TiRHr\nxmcFYq8tntAy7y78EDt2Fxpwe1UmLYgg4m3qTS1xkYxtjuFM8uHmBBPR1rak0kQgIm4J2WZsGEYx\n0DXIR9NN0yyq2WYDcKPdNuOIUplArBLIOgDy8gZGNEf04MFPUVJirRyUnb2Y1atHJUTwai7tQM37\nuTXvFyTMOYUyb95Krr32POq3A8+dW8w11wxxM1lNOJHOaK9rkRQWfZuxaZqDnElLnWRu8G9cUlq0\nKLI209WrRzFnTm0J6iLKyyspL690PL2xUFBwUb3Sn5X2/Pw19areoaRkNHPmONPm6GYnkvLyptW9\n5eV+z13jTqVzxIh+Nd9NnOsxEuqYZI/yyZ4uXewtoOPUOGNbkT/ZOdVm6sTYWLckctrD5XY7sF12\nxmOLiLsibjM2DGM48BBwOLDSMIw3TNO8wLGUSdJI1jbHRFlowOfzNapt8WY6RVKZJv1wkFOTXiRj\n9U+sOnAlY17FgvLJPuWVPconezQ3tQsSpaTkBs017K5QS01KfCTjiAJxjoKxwxR0oqMblvO8NpFM\nKkqGta4ltrRQRIry4iQQtTesKVOGMmXKUEaOXOaZtCWyRJ5IJlmk4oQ4Eh4F4xTk1aCnG5YE48UH\nRxGnKRinIAW91BLp0CYvBEGvPjiGK1GGwYl71GYsnpGsQ6DcFsnQJq+0cSbLIhXq3CktUTBOQV4N\nerphxU64HQuTJQh6iTp3SigKxinIy0FPNyypz6sPjiJO06QfHqTB9PYpr+wJN5+cmsDGCfEe7qZr\nyh7lkz2a9ENEIual2hPVlkgqUDBOQJoYQ+JBQVAkfhSME4xXermKiIhzNM44wWiMsIi7vDD+WpKP\nSsYiIjapZkpiRSXjBBPLmXz0xC+JIti1Gsn1G+53VDMlsaKScYKJVS9XPfFLogh2rS5ceAFjx74Y\n1vXb3DUPGTE/B5HGVDJOQLW9XPPyznMsWOqJXxJFsGt16tT8sK/fSK55zTEtsaKSsYiITV4afy3J\nRSVjAfTEL4kj2LU6a1Ze2NdvpNd8LGqmRDQdpge5Nc1cIk4moin57Em2fAp2rUZy/Qb7TrLlVawo\nn+yxOx2mgrEH6SK3T3llj/LJPuWVPcone+wGY1VTi4iIuEzBWERExGUKxiIiIi5TMBYREXGZgrGI\niIjLFIxFRERcpmAsIiLiMgVjERERlykYi4iIuEwLRYhIykrEKWAlOSkYi0hK0nrG4iWqphaRlKQ1\nvMVLFIxFRERcpmAskqD8fj/5+WvIz1+D3+93OzkJR2t4i5eozVgkATXX3qkOSPb5fD6eeWY4hYVF\nAOTkKP/EPQrGIgmoYXsnNe2dReTlneduwhKMz+dTnoknqJpaRETEZQrGIgkokvZOtTGLeJeqqUUS\nULjtnWpjFvE2BWORBBVOe2d+/jq1MYt4mKqpRUREXBZxydgwjHuBC7Earf4LXGGa5rdOJUxEnJOX\nN5BFi56ktDQPoKaNebi7iRKRgGiqqdcAU03T/N4wjHuAW4BpziRLRJykMbUi3hZxMDZNs7jey83A\niOiTIyKxojG1It7lVJvxlcAqh/YlIiKSUkKWjA3DKAa6BvloummaRTXbzAAqTNN8KgbpExERSXpp\n1dXVEX/ZMIw8YDww0DRNO7MIRH4wERGRxJNmZ6NoelMPBm4Gsm0GYgC2by+P9JApo0uXDOWTTcor\ne5RP9imv7FE+2dOlS4at7aJpM54DdASKDcN4wzCMR6LYl4iISMqKpjd1TycTIiIikqo0HaYktNrF\nD8BaPEFjZ0UkESkYS8Ly+/1ceumzlJTkAlr8QEQSl+amloRVWLixJhCnA+k1ix9sdDtZIiJhUzAW\nERFxmYKxJKycnH5kZxdgrVVSUbP4QT+3kyUiEja1GUvC8vl8rF49ijlztPiBiCQ2BWNJaFr8QESS\ngaqpRUREXKZgLCIi4jIFYxEREZcpGIuIiLhMwVhERMRlCsYiIiIuUzAWERFxmYKxiIiIyxSMRURE\nXKZgLCIi4jIFYxEREZcpGIuIiLhMwVhERMRlCsYiIiIuUzAWERFxmYKxiIiIyxSMRUREXKZgLCIi\n4jIFYxEREZcpGIuIiLhMwVhERMRlCsYiIiIuUzAWERFxmYKxiIiIyxSMRUREXKZgLCIi4jIFYxER\nEZcpGIuIiLhMwVhERMRlCsYiIiIuUzAWERFxmYKxiIiIyxSMRUREXKZgLCIi4jIFYxEREZe1ifSL\nhmHcCQwFqoEdQJ5pmp86lTAREZFUEU3JeLZpmqebpnkGsBy4zaE0iYiIpJSIg7FpmuX1XnYEvo4+\nOSIiIqkn4mpqAMMwZgK5wF6gjyMpEhERSTEhg7FhGMVA1yAfTTdNs8g0zRnADMMwpgEPAFfEII0i\nIiJJLa26ujrqnRiG0R1YZZrmqdEnSUREJLVE3GZsGEbPei8vBt6IPjkiIiKpJ5o24z8YhmEAVcB/\ngWudSZKIiEhqcaSaWkRERCKnGbhERERcpmAsIiLiMgVjERERl0U16UckDMO4F7gQqMDq+HWFaZrf\nxjsdXmcYxqXA74GTgJ+Ypvm6uynyFsMwBgN/AloD803TnOVykjzJMIwngCHA/0zTPM3t9HiVYRjH\nAIuAH2DNt/+YaZoPuZsqbzIMwweUAAcBbYHnTdO8xd1UeZdhGK2B14DPTNO8qLnt3CgZrwFOMU3z\ndOA/gH7E4N4GhgMb3U6I19Rc3H8GBgM/BC4zDONkd1PlWU9i5ZOEVgncYJrmKVizCU7QNRWcaZp+\noH/NugQ/AvobhvEzl5PlZZOA97Ae8poV92BsmmaxaZrf17zcDHSLdxoSgWmaH5im+R+30+FRvYCP\nTNPcappmJVCINdZdGjFN8xXgG7fT4XWmaW4zTfPNmv/vAd4HjnI3Vd5lmubemv+2xaqd2ulicjzL\nMIxuwM+B+UBaqG3jXk3dyJXA0y6nQRLP0UD95To/A3q7lBZJMoZh9ADOxCosSBCGYbQCXgeOB+aa\npvmey0nyqgeAm4FOLW0Yk2Dc0pzWNdvMACpM03wqFmlIBHbySYLS4HiJCcMwOgJ/BSbVlJAliJra\nzTMMwzgYeMkwjHNM03zZ5WR5imEYF2L11XjDMIxzWto+JsHYNM1BoT43DCMPq+g+MBbHTxQt5ZM0\n63PgmHqvj8EqHYtEzDCMdOA5YLFpmsvdTk8iME3zW8MwVgJnAS+7nByvORsYahjGzwEf0MkwjEWm\naY4JtrEbvakHYxXbs2s6AkjLQrY1pKDXgJ411YlfACOBy1xNkSQ0wzDSgAXAe6Zp/snt9HiZYRiH\nAwdM09xlGEY7YBBwu8vJ8hzTNKcD0wEMw8gGbmouEIM7vannAB2BYsMw3jAM4xEX0uB5hmEMNwzj\nU6yenSsNw3jR7TR5hWmaB4DrgJeweik+Y5rm++6mypsMw3ga+AdwomEYnxqGoWVOg/spMBqrZ/Ab\nNf/UCz24I4H1hmG8idWuXmSa5jqX05QIQjavaW5qERERl2kGLhEREZcpGIuIiLhMwVhERMRlCsYi\nIiIuUzAWERFxmYKxiIiIyxSMRUREXKZgLCIi4rL/BwHwT/vyMIPsAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.scatter(points[:, 0], points[:, 1])\n", "ax = plt.gca()\n", "ax.add_artist(plt.Circle(np.array([1, 0]), 0.75/2, fill=False, lw=3))\n", "ax.add_artist(plt.Circle(np.array([-0.5, 0.5]), 0.25/2, fill=False, lw=3))\n", "ax.add_artist(plt.Circle(np.array([-0.5, -0.5]), 0.5/2, fill=False, lw=3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now write a function that initializes k centroids by randomly selecting them from the data points." ] }, { "cell_type": "code", "execution_count": 104, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def initialize_centroids(points, k):\n", " \"\"\"returns k centroids from the initial points\"\"\"\n", " centroids = points.copy()\n", " np.random.shuffle(centroids)\n", " return centroids[:k]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try the function:" ] }, { "cell_type": "code", "execution_count": 105, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-0.48840901, -0.14088659],\n", " [-0.53984176, 0.53868892],\n", " [ 0.98947381, 0.30255966]])" ] }, "execution_count": 105, "metadata": {}, "output_type": "execute_result" } ], "source": [ "initialize_centroids(points, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's represent the results using a plot:" ] }, { "cell_type": "code", "execution_count": 106, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 106, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFXCAYAAACRLCZbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4XNV57/GvEbKnJAKeEgVIII2eQlZza9zU8SUXu2Bw\nTHyJXMcgfHxRXJyQJn4Il5iLmzQtdVM7BELcxJxgjGRxiIwDcmwMxgKDnNPK5nCAPHnidOVyRBvK\n5Ti0gPtwxlKEzh97RhqN5rJn9p5Ze8/8Ps/D82BpZs/aa7b2u9da71pr0sjICCIiIuLOSa4LICIi\nUu8UjEVERBxTMBYREXFMwVhERMQxBWMRERHHFIxFREQcOznIm40xCaAPmAJMBn5krb0xjIKJiIjU\ni0AtY2ttErjAWjsV+GPgAmPMx0MpmYiISJ0I3E1trX0j9b+TgQbgP4IeU0REpJ4E6qYGMMacBDwN\n/CGw1Vp7NHCpRERE6siksJbDNMacBjwC3GCtfSKUg4qIiNSBwC3jNGvta8aYfcA04IlcrxkZGRmZ\nNGlSWB8pIiISdb6CXtBs6rcBv7PWvmqM+T3gYuBv8pZo0iSOHTse5CPrQnNzk+rJJ9WVP6on/1RX\n/qie/GlubvL1uqAt47OBztS48UlAl7X2sYDHFBERqSuBgrG19qfAh0Mqi4iISF3SClwiIiKOKRiL\niIg4pmAsIiLimIKxiIiIYwrGIiIijikYi4iIOKZgLCIi4piCsYiIiGMKxiIiIo4pGIuIiDimYCwi\nIuKYgrGIiIhjCsYiIiKOKRiLiIg4pmAsIiLimIKxiIiIYwrGIiIijikYi4iIOKZgLCIi4piCsYiI\niGMKxiIiIo4pGIuIiDimYCwiIuKYgrGIiIhjCsYiIiKOKRiLiIg4pmAsIiLimIKxiIiIYwrGIiIi\njikYi4iIOKZgLCIi4piCsYiIiGMKxiIiIo4pGIuIiDimYCwiIuKYgrGIiIhjJ7sugEg9SyaTdHcf\nAqCtbTaJRMJxiUTEBQVjEUeSySSXXdZDf/9nAejpuZudO5coIIvUIXVTizjS3X0oFYgbgUb6+9tH\nW8kiUl8UjEVERBxTMBZxpK1tNrNm3Q0MAoPMmtVBW9ts18USEQc0ZiziSCKRYOfOJXR37wWgrU3j\nxSL1KlAwNsacC+wA3g6MAN+31n4njIKJ1INEIkF7+zzXxRARx4J2Uw8BV1tr3w/MBL5ojHlv8GKJ\niIjUj0DB2Fr7krX22dT//xfwc+AdYRRMRESkXoQ2ZmyMeTfwJ8CRsI4pIv5pARGR+AolGBtj3gr8\nELgq1UIWkSrSAiIi8TZpZGQk0AGMMY3Ag8DD1tpvF3l5sA8TkZzuuGMfX/jCPLwFRAAG2bq1lyuv\nXOCyWCICk/y8KGg29STgLuCoj0AMwLFjx4N8ZF1obm5SPfmkuvIcP57M+bN03aie/FNd+aN68qe5\nucnX64JmU38MWAFcYIx5JvXf/IDHFJESaQERkXgL1DK21v5PtIqXiHNaQEQk3rQCl0iN0AIiIvGl\nVq2IiIhjCsYiIiKOKRiLiIg4pmAsIiLimIKxiIiIYwrGIiIijmlqk4hIhWkTDylGwVhEpIK0iYf4\noW5qEZEK6u4+lArEjUAj/f3to61kkTQFYxEREccUjEVqTDKZpKPjAB0dB0gmJ+7mJNWlTTzED40Z\ni9SQfOOT4G8bNwmfNvEQP9QyFqkhGp+MpvQmHu3t8xSIJScFYxEREccUjEVqiMYnReJJY8YiNUTj\nkyLxpGAsUmPS45MiEh/qphYREXFMwVhERMQxBWMRERHHFIxFREQcUzAWERFxTNnUInloD1oRqRYF\nY5EctAetiFSTuqlFctAazyJSTWoZi0gkaFhA6plaxlLTyt3bt5JrPGu/4YnSwwLr1y9m/frFXHZZ\nj+pG6opaxlKzgoz7VmqNZ41F5zZ+WIDUsMBeLespdUMtY6lZQcd9K7EHrcaiRSQXBWORmKnFbm5t\n/Sj1btLIyEg1P2/k2LHj1fy8WGpubkL15E+huhrrEm4HYNasDuddwkHLlN3NPWuWv27uUq4pV4lU\nUUng0t+fP6onf5qbmyb5eZ2CcQTpIvevWF1F5QafKUiZOjoOsH79YtJjqzDI5s3Fx1b9XlPlBvta\nor8/f1RP/vgNxkrgkpoWxb19o1imNCVSibihMWORGNHYqkhtUstYJEYqNeUqra1tNj09d48b025r\nWxLa8UUkNwVjkZgpp5s7nYENhcepKx3sRSQ3BWORGpdMJlm2bBd9fSuB4guNRHlMW6RWacxYpMZ1\ndx9KBWItNCISVWoZi0jkRHFKmkglKRhLbNTjDTqMc25rm82+fV309a0Aop+UpfW7pR4pGEss5LtB\nQ5PbglVQWEEpkUiwf/9ytmyJR1KW5jpLPQo8ZmyM2W6MedkY89MwCiSSSz1usBDmOVdi0wsRCU8Y\nCVx3A/NDOI6ISKwXNqnFTTykOgIHY2vtj4H/DKEsInnF+QZdrlLOuZaCQHqu8+bNe9m8eW9sxovT\nwwrr1y9m/frFXHZZT+y/C6kejRlLLNTjYhR+z7kWE57iONdZY90ShOYZS2zU47hnIpEYbQ13dx/K\n2dKqx/F0kVpT9ZZxc3PtZr+GSfVUXDKZ5I479gHQ3j63JgN09upZ+/Z1sX//8nHn2tQ08bybmhIT\nriFdU/6VU1fr1i0YN4Vszpx7WLdueU1el2m6psITyn7Gxph3A3uttR8s8lLtZ+yD9gktLir77lZ6\n7rOf/YvH6qId8OYRZ9eFrin/gtRVPc2F1zXlT9X2MzbG/ACYA5xhjPkN8DVr7d1BjyuVF7UbRynl\nicL4XJCx2jDrvh7H06MqjmPdEg2Bg7G19vIwCiLVFbWkn6iVx49yHwhKOVe/WxpWIwhE7eEtW9TL\nF1Stn1+9UwJXnYpa0k+p5YnzVKdSzjUq03yiPm0n6uXLVM40tDidn5RHwVhiKR2ktm7tdRakqvVA\nEIUs8mo9vJU7XzpqD5f5lBtU43J+Uj4F4xrk54YWlZZluqxDQ4PMnHlXSeVJJBJceeWC0IOU34BQ\nbqs1KnUfNfXQ+lNQlXy06EeN8TseGYWkn+yyzpjxfTZufIDGxsaKlqfQ2FupY9fljNVGoe5L5Xfs\nOoggSXnVKJ9LtX5+omBcc0q5oYWV9FMsuOX7XXZZjxxZy9Kllc2ILhZsq5WlHbes26g/QES9fGnl\nBtW4nJ+UT8FYJigla7NQcCsW+IaGhiYcL9fPwhSFKVFxVekHiKCtvzg84AQJqnE4PymfxoxrTNDx\nyFLH7QqNgRUfHxsBOkfLCjtSPwtPqQlBGs91JyqZ45UWhYQ8iR61jGtM0O6sarYcGxsnAwuB3tRP\nLqWxsbfAO0qTq2Xe2XlJwdaXugPdytX60/zacKgeo03BuAZVsztr2bKP0d39LZ5++jPAO5k1q2s0\nuBXrdvR+/4OKJaXkerDYvXtv0WAb5+7AKNxwwyxDtRaDySzzunULQj12FMRxUZ16o2Bcx3LdNP2O\n2504cYLHb/kGiUd72fYLy3DDV3nq7ecw5cOfZtKkPweKtzJdtULjHGxzSX+PQ0ND7NnzCkeOfA5w\nc8MN+6ZfjZ6a7DLv29dFV9eimgpUypWIPgXjOpIZfFtbp7N69cM5b5rpADk0NAicSnf3oXEtnBMn\nTvDg6uWsOdjL5IzjT3vxOQa/ezvbf36UhZ33MmXKlKKBr5KBsR6mg2QHEm8MfhhIlLQ8Zzkt2Vzv\ni+NNP7vMfX0rIl9mqT1K4IqAclcdKvUzMhOzPvnJe+nvv5zM5KqurseAsT10H3zwOBs2LJ2QyPX4\nLd+YEIjTJgNrDvby+Lf+oSLnUYp6SAjKTpKDVcBjvt9f7kIb1VqgQwl14VA9Rp+CsWPVuqll37QH\nBq4GDo57zbZtz45+dr5M6OHhYRKP5g7EaZOBxKMHGB4eDq385T6w1Gfm6hB+b7jlrgiV731h3/Sr\n8UCVXeY5c+6puUBVDw+mcaduasdydet1dPSydGl4N4NkMkl//1Fg8bifn3LK/bzxxgVAAriHgYGb\n6O7uLdg99+KLL3D+L23RzzzvF5YXX3yBc845N1jhqd3kkzASnbK74mfO3M6iRafR2LjXSSZ4JfIA\nKj3Gn13mdeuWc/x4Zee7u1BruRK1RsG4xo0Fsi/jjSeuSv1mB2+88V3ge8B7gOVkdpTkG2/97W+P\nVbP4QG0mn4T1gDEx+C0t6RjljqsXel8cb/qZZU4kEjmDcRQy1aV2KRg7luum1t6+KrQn8/GBbAXw\nTeADqf9PAFcDDwEnTbih5mrhnH32O3jivPOZevRnBT/3V+8xfOLsd4RyDtVQ7RttmA8YQYJfIpGg\ns/MSrr9+EwCbNrX7Ovd6m49dq70zEh0Kxo7lu6mFFYy9jOiH8L7qucD7gU+RDgIAS5b8lFmzhnNO\nPcq+yTc0NPC/z3wPf370Z3nHjQeB5EXzaGhoyPn7ZDJJV9ejPPXUr5g27XxWrpxb8KZW6azoqN5o\ngzwgZL8XmvK+zsuqvwGAl17yf+5BHgLi1sqsxd4ZiRYF4wioVLdeMplk797XgM+kftLJRz7yXzQ0\n3MXhw2sAL7DdfvuVE26GhW6Wv3/xpfy3x9/gf/DIhIA8CGy/8GIWXXdj3jItW3YfR45MAW6gpwf2\n7r2L++7L371a6VaYixttsQeMIA8Iud578OCqnK91ce7lnFvcgrdIqRSMa1h39yEOH/4LxlrBq2ht\nfYCVK+cWDGzFbpYrVszlRz96hVlHPsSneIRF/IyTG0/i18aQvGgei667kcmTc7ebu7sPceTIOWS2\nzg8fXlM0AMRxHLKQYg8YQYJkNZICgyj13KLQc1EPc9bFLQXjOtPY2Fg0sBW7WSYSCXbtupTu7kO8\n+eaf8uac8xkE/k/fLznppAbefPPNssrmqvXj6kYbhQeMOASZKHQR19sYuVSfgnENq+SNNjOQlNJy\naW2dzs03b+P48WOkM7unT99Ga+uigseoZKCO4o02yHdXSlKgi3OPwwNALlF4eJLaNWlkJNwt64oY\nOXbseDU/L1bSAaepKcGCBdNDuSlmrlsMIzQ2Tqa1dTq7dz8J5A5sY8G1HfBultnBNTM4Dg0NsmHD\nZxjrDh9k8+bcLZeOjgOsXz8PeBSwwHls3DhEY+Nk1q9fnPMY2cF+1qyxQN3c3EStXlNBE7gyk+Su\nuabV2dzZXOcxdl0OApNobGzMe45+rscw1fI1FSbVkz/NzU2T/LyuJlvGcUv2SN84t237NwYGrgHG\nB5wg0ktbZgazv/u7Tbz++lVAImcrtlhrKTs4trTcCiTJzNAuUiog3RIa5KmnNjFt2nkTXtXff5S2\nttmR6KZ0IWhL7MEHj9Pf7yXJHTjgZvODQr0m2ddlvh6VKPZc+BW3e5G4U3PLYVZrecmwpMu7YcMp\nqUBc2rKEfmQvXfj66+uBQwU/p9AykrmW1mxp2YifJRDb2mYzc+Zdo6+FHfT0fJm9e19jxozvT/j5\nZZf1pFr1aUngIfr7j0b6e3Ut+zvyNj8I53oKUo7M662UpTjjuKxp3O5F4lbNBeNy19p1Zay88emk\nGB8cPVdcMdX3urcjI8N4i488jLf4SBOHD6/hzDOP8aEPfRnYO/pzr2tyJLV28OvAPcACenpuiNzN\nrRobfkj5qv39xO1eJG7VXDCOr7lAF0EX2M91w8leCP/UUzcDnyjrc5LJJHv2vIK3tKZ3vJkzt7Ny\n5VxfLRdvatPngKl405vGXrtnz5/wk5/cDvwKL1B75W9snMzOnUtYsuR2YDWZN7eODv87FFVS1FpB\nUdn8oNDGEdXcSShq349ItppL4Kp2skdQ48ubpKXl77n66pm0tn6s5DIXSnTKTJgZGvodP/nJc0yb\ndh4rV15U0ud4CViL8fbMfQwYYuPG/8fatQvLeP+9eC1ggB2MLdE5iBeMjzFjxgl27bqURCKR8d6x\nJK+tW6MxfzZX2fIlsVVL5njlunULIpXA5ed3YSrl+wkrMSlu96JSKYHLn7pN4IpbssfE8l7Luec2\n+77Is7Oa8yU6TUyYSfLssxuBSeOWo/R/c0wAC4BBGhv3+j7f8dNaPkNLy9f44AcT7NlzLZmtZO/S\nXMXixQ+MliHflJ1jx44rSSYHP5sfVLscpfwu7uJ2LxK3aq5lXAv8PnHmymoeGPhLxtYhHv/0P7FV\nuhIYa0EDeVvWEz+zPfWa0p/2J66bzLhjeuPC3i5S2a2X7Pc2Nzdx4YU7xtXBFVe8q+QWf9AWWtRb\nQfXeiinl+6n3uvJL9eSP35axgnEE+b3Ic3W9tbR8jYGBvwXG33CSySRXXXUHPT03AAeAeWR32QG+\nuvL8BK5Sg5s3vesxtm17loGBm4CEr4B2//2H+MIXxp8LPMysWS+XvZZzudPKojaNJV83ddTKWS1+\nz1tBxh/Vkz91201d7664Yupot3G6W2zinsZvC/QZxboWy1lLOJFIsHbtgtS62b3jyl+6kwOv5VzO\nPOYodblmfwf79nnzjIGSv5taEaXvRySbgnGM5RpDXbly4o21q+tR+vvPxJtb/BngAGeccTOvvPJX\no+9LL0cYxjKFXnC7HK8FDv39bXR39/q6Efq5YWa2cNauncuOHXfn6OKub9kPGN48Y+8hrR4XUBGJ\nOgXjGPOTIJJMJtm27d+Aa1I/6QI+wzXXDE9oQQOhJJx4SxzuxJuGBNDJ0NBbSj5OLrlafJ2dl7Br\n1wMZXdwnMXPmdoaGTqOj40DRrti4rpUsIrVDY8YRFOZYTL5x5b6+ayvWNXnnnfvYsGHpuM/cuPEB\n1q5dUNJxco3xFZqikrkO9549r6TmM/sbA661cdTshKU5c+7J6qb2fh61RLMo0FioP6onfzRmLHld\nccXUit58GxsnrlGd62eZ8mdYjx/bLCTdxd3RcSAViMd3xabXuU5/RvZ63JldtS+//BIrVtwGwD33\nXM2ZZ55V8LOjJrvXZO3aP+fOO70FUjo7L2H37rHeD/Ae2rx/x/9BRCSOFIxrXL5x5Wp/ZqFu31wJ\nXwsXnppzbDP72N7KUouKlmloaMh34tLLL7/E1KldDA/fDsDUqTfz7LMrYxmQ0z0Gra0P0NfnTWXL\nPPdyku2KqbVeBpFq0HKYNShzSUzwxoH9rhsdhnSrzO9n5lrD96mnfpn39QsXnsqSJZvYuPF+9u9f\nPuHYuZZZhBHf6wSvWHEbw8NfHX3t8PBfjbaSsw0PD/P887/h+ed/w/DwcN4yu9TdfSgViINt1uCH\nlp0UKY9axjUmX0snbtmy06adx0svjW9dt7ZeMu7cXnrpbq65ZuJ7cyW2hb1A/4kTJ3j8lm+QeLSX\n839pAfjx+YbkRRdzwXU3MmXKlFA/Ly7qdbtLkaCUwBVBxRIjCnUDRmGN5FIX0ci3OhIw7jy7uw+V\nvTZ1KSswjXVTe1O/Ghr+blw39YkTJ3hw9XLWHOxlctZ7B4HtF17Mws57Qw3IQbp+k8kkK1fupa/P\nWwc8ezGY7HrxxpSfLOuzonD9BaXEJH9UT/4ogatGVWKML2ylto4KTdEK6yZeyjrBZ555Fs8+u5IV\nK74MTEzgevyWb+QMxACTgTUHe9n5rX9g/k1/HUrZX331VT75ye8yMDALuJCenh+U9J0nEgn271/O\nli0Tzz27XlpbL2H16ofLvr40TUykPBozjpliY3zV3Jau2nKdW3v7XN/vL2WD+jPPPIve3k309m4a\nF4iHh4dJPJo7EKdNBhKPHghlDDmZTPLJT97LwMDNeFtO7kotolJat3uhc8/83e7dTwYaQy41X0BE\nPGoZ15go7BQTRjZ1rpt4vnOr5m5EL774wugYcSHn/cLy4osvcM455wb6vO7uQwwMXMNYt+8K4KFA\nx6w0LTspUjoF45jxE+iqcTMsNIbp94EgfYz+/qP0999Aeo/k/v6309X1aM49kuvlRu9tnPEo9933\n49RPFpDeYrKl5TBtbddW5HNLeZDSFCaR8ARO4DLGzAe+DTQA26y1mwq8XAlcPgRJ4KqGMHY5Gn+M\nh4ALgV2kt3VsabmVvr41RY8ZJImknHocHh7mxxfNZtnPflrwdbs+8EE+0XuIhoaGssq1bNl9HDky\nhbElRe8G2mhp2cojjyzn9NNPL+mYpdST3125wtjpKoqUmOSP6skfvwlcgcaMjTENwD8C84H3AZcb\nY94b5JhSXCljn5lzjsOa7xnG3NTxx/gk8Pd4gdg75sDA1aFPR8pU7nzYhoYGkhddzGCB1wwCyYvm\nlRWIwaubI0fOwQvEjan/2lmy5Hb6+taUHIhLlX195bqGwp6fLFLvgiZwTQd+Za19zlo7BHQDnw5e\nLAlDfBZgSAB/WvAVYT9UBAkmF1x3I9svzB2QB4GbzVQuvO7GwGXMNmvW+5z1gkT/GhKJt6DB+J3A\nbzL+/XzqZxIBlWq9hJGxnX2MGTNeZubMu3IeM2oBYcqUKSzsvJedX76W7vd/gCc5mcNM5qt8iGXv\nvIT2fQ8yeXKhfOvC2tpmM2PG83h7T3v1MXPmdidZ8fmuoVrO2hdxIWgCV8kDzs3NTQE/sj6EUU9N\nTRNbUU1NiQnH9lqd3iYC7e1zfbS+mjh4cBUdHb2p7RLfzr59T/p878RjeJ/7FwAZ/141eqw77pg4\nb3n37odGN58o7XM969YtYN++rtGFMObMuYd16yYurVmo/Ctvu4Xh4U38+te/5oc//CfOPv0MNqwp\nPnTg59hPPPEX3Hnnw/zzP9/KRz/6R6xd2x74uOVcU/muoXPPbc76/lbVxHhxmu5T/qiewhMogcsY\nMxP4urV2furfNwJvFkjiUgKXD2ElRvhZdSpIIk61knjybQPpzb0t/3NdJ8JVU1NTI1u27ANKO1e/\nq6PFue6yr4Nzz21WYpIPSuDyx28CV9BgfDJggbnAC8CTwOXW2p/neYuCsQ9hXuTFAk6Q5QurtfRh\ndkBoabmNgYG/BNJP5fFbcrGaxpbD9DLVS314yb6GxhYiuaas40VJrgfKgwdXVXXuelwpGPtTleUw\nrbW/M8Z8CXgEb2rTXQUCsThQC/Nys+ctDw29iw0b4nfjd2X8rk2lb96QeQ15gfi7qV6J+G8GkWvp\n1o4Of+udV1s99eTUo8CLflhrHwYeDqEs4kCQtYSruQ5xdkB48MHKfm41b3xxusl6K4LNcl2MuhOH\nNeklGO3aFEHV7v4pFgwK/d5VIEl/blNTggULpof6udVc0KIan1Vo16ZSeUMT8/AWaPGO19Jym68F\nWsIQ9vWWa0w8it3UUdwNS93U/lRlzLgMCsY+ROkid73SUrGbbyXqKuiNr5SA4fez/K6K1d19KJXh\nPonGxsbR15abwJXrM7zroQ04SEvLYR555IsVX4hk/GeHey3GIYFLwTi+tIWihMLlZvFx7JoLq8yZ\nAaK1dXrObQ2BHK+5HNhJehnN9Gubm5tC+c7Gxu+9KU1tbddW7fuo1LUYh7wKbU1Z+xSMKyhOY4FR\n5OpBoNiNr9D3WmqZc31Wa+slXHrp/Rw+7M29/t73buG557407pjbt3fT2fniaEbz1q23pjLMDzG2\njObY569bt4COjgM5y1yqOASvuMt1jbnejU0qS8E4ZJldhXv3vjZ6Q41Cqy77DxyKzxWtxyfyQje+\nsFvruT6rq+ux1HXjBdTnnrsWL0dybKXZLVv6eOWVO0ZfMzBwdeo1E/+kh4YGmT//3tGpTVG4FqH0\nh9V6uRYLXWN6CKpdGjMO0cSdiD5FOWM8lRiLyf4DnznzLkZGhjly5HNA4fG3sFr4pR7Hz6Il1R63\nKjZ256fMxXz+89+mp+eGcZ9xxhnX8cortwDwB3+wmX/91xHgj/E22UgwthDKTcB9wKrRz1+48FQ2\nbFiat8wulDv+W43eJtdjoVEcH87FdT3FhcaMHRjfRRmtqs3uPj18eA1eS6p4d2oYT+TltCgTiQSd\nnZdw/fXegm6bNgVfErLSwuhOnDbtPHp6OkkHVNjBunV/yimn7GVoaJA772wArkv9rhO4NLW14hfZ\nvbuXoaG3AA+kEriWRHI3pXKHINQ6lFoVdKMIyWsumQv91/tC+uVsWpFMJlm9+mF6em6gp+cGVq9+\n2PmOQX42SChli8tcli37OO9+twW+CexmxowTrFnzKdrb59HYOJnnnruOsa0VV3HGGVeP7nHc3j6P\ntWsXsnbtgtHPb2ubzZw5XQXLXG1DQxOnDuX6WT3SJhz1KVrNt5gbP6Z1EjNmnGDx4rEWSiW7gEsr\nG8ycuT3VTe1tBBjF8bewE7jCqOtKJ9KkH0Cee+5rwEHOOOMxtm1bV/Azrrnm4oJTixKJBPv3L2fL\nligl/4zgPaymW/93ceTIy+OmY9UrJWvVJ40Zh6yUG36+cbNKzXMsJ4ErzM8udSw139hZW9vs0XKv\nW7fA1wINY59/Od782P6qzY8txfhFNbyEq5aWW0cX1Sh3TDpq43tj53kIrwV4DHAzlz1b1OoqqlRP\n/mjRjxjIF2y+8pWlNXmRh5HA1dl5ybg5t3PmdNHVtajosYoFuajwynkyhZL/ymnhR+3GOf67fQhY\nQFQSlqJWV2ELqzeu1uspLErgkpxczn0uNfkmV3dddtd1X9+KCV3X+c/xIF4gHpsOFLUNDtraZrN1\n67cYGPhU3tfUQhJT5nfb33+Unp4FrotUF+K4kE69UAKXQ9VO1Ej/Ia5fv5j16xdz2WU9zhOiiik1\nGSrfOba1zaal5Z8mvD5qSUOJRIJHHvkiLS23UusJPOnv9vbbr1TCUpWUk0gp1aGWsUPVTtRwubRl\nWLIT0ebMuYe2tkWjvy90jqtXv5+vf/1W4OrUq3cAbym5DJXuXTj99NPp61sTywSeZDJJV9djPPXU\nL5k27TxWrryoaNmVsCSiYOxcLXQ5VkK+gJd94163brnvHXZOOeWtwF8CvamfXEpjY2+Bd+QuVzW6\n+ap5XYS5qEvmMp49PZ3s2XMfu3Zd6isg6++g8uplFbM4UgJXBFUqMSKM1aGqMeZcyupM2XVV6BzD\nOP+4rI6ULd81NVYnq4B/58Mf3sX997fzlreU3mOQq27gITZvHo58/WSq9cQkJXBVlxK4ZIKg3YHV\nahUG6U6lCnTpAAARQ0lEQVQvdI7qDp3onnse40T/r/kqM1nAz+FpuH/qt3npwx/l89//Dqeddprr\nIkrI1AsRTQrGNajQk2+QP8S4jDkXOsegN6Ja6uY7ceIEb3Ruop9nmZzx8xmvvczg4z186SO/YP3/\n2u87ILe1zeaBB+5KLbUKsIMZM07Q1nZp6GWvJu2+JtWgYFxjamHqQpCAV+kbZy21rh+/5Rt81Y4P\nxGmTgX989Wds+NxVrN/Z4et4iUSC++5bSlfXAxkJXJ+Obf1Abfw9STxozDiCgozFVHJMM4wx11I+\ny09QzayrcncCqgfZ19Tw8DBPzP04bUd/VvB9d5z9blqffoaGhoZKFzEyMusqrjkC1aAxY380Ziyh\nq2arsJzuZBfd6HHtwnzuuQH+8Kgt+rrpv32BF198gXPOOXfC7+J67iJRpGBcYyo9pukq+SN94x8a\nGgQmjW4oAE1VL0tmmeLQhZlMJrnjjkMcP54cDZp79hzh0z7W/JmU55k+LuceVC3lCEi0KRiXIA4t\ngVoa00wbv8nDTmA14AWAgwdXjb6u2jfOrq7HIp/Qli9onnba77OP9zKDnxR8/6/eY/jE2e+Y8PMo\nJ/OF+Xdai39PEk0Kxj7FqSVQa1MXxm78B/AC8VgA6OjoZelSb+nEat44k8kk27Y9CyytyPHDki9o\nLl9+AV/acjeD//6TnAlc4M0STl40r8B4cRLvOwH4RKjlLtVYz8kQe/a8wpEjnwPC+Tuttb8niSat\nTe2T1nSNh1LXsi5Xd/chBgY2AF2k11RuabktNmsqJxIJvvXEndxspjKY4/eDwPYLL+bC627M+f7W\n1umceurtwDxgHqee+h1aW6dXsMT5Za5HvmHDUo4cmQIME4W/02QySUfHATo6DkR+HXhxS8FYIm9s\nQ41P4G1IP7ahQHv7XIclSwDL8ZbXfIgrrnhX5HpKCm1Gctppp/H5R3vZ+eVr2fWBD/LM5Mk8M3ky\nuz7wQXZ++VoW7fgBkyfnbjfv3v0kr79+PemH09dfX8/u3U9W67TGyX5QhlXAY07KkimOG7OIO+qm\n9ims8cg4jDuXq1LnNtb93MvQ0FuAB1IJXF73o9+1qcM0/nq4mFmzOli5Mtj4dCXqL113+/b1phK4\nxnfZTpkyhfk3/TXD1/8VL774AgCfOPsdNTCVaYixhw83CVdRHleX6FEw9imM8cg4jTuXqtLnFsVx\nu4ULmzjrrE1Mm3Y+K1cGO9dK1l8ikeDKKxcUnBPa0NCQc/pSPtVIlvP7cJJdlpkzt7No0Wk0Nu5V\nwpXEhhb9qCK/CwjEcTK9i8URkskk+/Y9OW7KTrnHKaVFWonFRSpdf5W4porVW5CWfql1HGavQlh1\nVc1FclyI433KBS36UWfi2v1dbrnDaklmH+eBB+5Ktaom5y2Puh89hXorgn4/pdZxFHtONC1KSqEE\nrioqlEwTRBQSRco5tyDlDiu7Pfs4hw+vYcOGU6pej5W6NlzR7ANPtbL7Jf4UjKso/aS8efNeNm/e\nG1qXVRRufOWcWxTKndvJFCpPJQJnpa6NuKrWw4mmHklUqJu6yqLYnRaWap5bWAlE2ceBHcCKgu+p\nVPdjVK6NMIY8gn4/1ejizdeV7nKJValfSuCKoFITI+KaKJKr3J2dl4zOVy0WCMJO4BpbvWntaHni\nUI9++L2mwkxOyxXU/Qb6auRA5Eua+8pXlioxyQclcPnjN4FLwThC0jegpqYECxZMLzn7NO4JXK2t\n01m9+uGSAkGxG0I5mdJxrMdi/N44q7MFZ+Hvt1pbYSoYB6Ng7I+yqWMm6A0oKl2cpcosd0fHgVCz\nlMvJ6A1aj7UazMPgN0O6Wtnq2pFJokQJXBER3WSmaCuUgBNmnfpJ9IlCVntQtZbVXYiS5iRK1DKW\ngqrZzVtqSyWZTLJs2S76+lYClVvRzG8LO1+Lrq1tdmxay5VMnPL7/VazxRrXHiWpPRozjogoJmGV\nswpS0LG+UoJ5sfHNsOrU7zhqrtdt3Hg/Dz54vOLjn8VEZXwvSglc+USlrqJO9eSPxoxjJrNF4iVw\nue8yK3XsLoyxvjBbKtVeASlXiw5O1WpdGfx+v2qxSr3RmHGEpG9AV165wHkgjoO2ttnMmTO2n3Cu\n8c0wVkDyO46aawyysbEx4xVJ4CH6+4/GbixZRCpL3dQRFJXun1K7eV10tTc1NbJlyz6gst2ZwdfQ\nbgPuA1YD1e+uDnpN1VOWeFT+/qJO9eRPxecZG2OWAV8H/gj4iLX2aR9vUzD2oVoBxo9Sb8Kvvvoq\n11/fAcCmTe2cfvrpFS1fHG4IyWSSq666g56eG6jmrlaZgtRTteb9RkUcrqkoUD35U40x458CS4D/\nHuAYkqVaGcJ+lTJ2l0wmU4t23ADASy/V9k3br0QiwaxZ76Onx3VJyqNdqkQqr+wxY2vtv1hrfxFm\nYcS78XmBOH7zjetxrrTfjQZcz9/Vhggi0aZsaqmKWhxzLGWFL5d72wbtbdFKVSKVV3DM2BjTC5yV\n41c3WWv3pl7zOHCt3zHjskoZI14L5DEA2tvnlrXZ/fz599LX5+0cNGfOPezfvzwWwStf2YHUz1em\nft4Vm3Mq5I479vGFL8wjcxx469ZerrxygctiTRBGOYNe1yJ1LPiYsbX24nDKMqaWB/yzW0o7dpQ3\nZrp//3K2bEm3oBZx/PgQx48PhV7eSujqWpTR+vPK3tFxIKPrHfr6VrBlSzhjji6TSI4fn9jde/x4\nMnLXeFjlXLp0duq98bkey6HEJH9UT/40N/vbkjOseca+In+tC2vMNIy5sa7Eueylcj0O7Jef+dgi\n4lbZY8bGmCXAd4C3AfuMMc9Yay8JrWRSM2p1zNHlOHApEolEVm9LNMspUs+06EeIwlr0oha7fyqV\nwFWLdVUJqif/VFf+qJ780drUDsSlpeSC1hp2Kz21CWonmz1uanFGgYRHwThkCjrB6IYVvqgtJFOP\nSpkGJ/VJG0XUqSguApG+Ya1fv5j16xdz2WU9kSlbnMV5IZlaUY8L4khpFIzrUFSDnm5YkksUHxxF\nwqZgXIcU9OpLuVObohAEo/rgWKq4TIMTdzRmLJFRq1OgXCtnalNUxjhrZZMKJXdKMQrGdSiqQU83\nrMopNbGwVoJglCi5UwpRMK5DUQ56umFJpqg+OIqETYt+RJAm0/unuvKn1HoKawGbMFR7upuuKX9U\nT/5o0Q8RKVuUek/UWyL1QME4hrQwhlSDgqBI9SgYx0xUslxFRCQ8mmccM5ojLOJWFOZfS+1Ry1hE\nxCf1TEmlqGUcM5VcyUdP/BIXua7Vcq7fUt+jnimpFLWMY6ZSWa564pe4yHWtdnZewurVD5d0/ea7\n5qGp4ucgkk0t4xhKZ7m2t88LLVjqiV/iIte1ev31HSVfv+Vc81pjWipFLWMREZ+iNP9aaotaxgLo\niV/iI9e1umlTe8nXb7nXfCV6pkS0HGYEuVpmLo6LiWhJPn9qrZ5yXavlXL+53lNrdVUpqid//C6H\nqWAcQbrI/VNd+aN68k915Y/qyR+/wVjd1CIiIo4pGIuIiDimYCwiIuKYgrGIiIhjCsYiIiKOKRiL\niIg4pmAsIiLimIKxiIiIYwrGIiIijmmjCBGpW3FcAlZqk4KxiNQl7WcsUaJuahGpS9rDW6JEwVhE\nRMQxBWORmEomk3R0HKCj4wDJZNJ1cWJHe3hLlGjMWCSG8o13KgHJv0Qiwc6dS+ju3gtAW5vqT9xR\nMBaJofHjnaTGO/fS3j7PbcFiJpFIqM4kEtRNLSIi4piCsUgMlTPeqTFmkehSN7VIDJU63qkxZpFo\nUzAWialSxjs7Oh7TGLNIhKmbWkRExLGyW8bGmG8CC/EGrX4NfNZa+1pYBROR8LS3z2XHjrvp728H\nSI0xL3FbKBEZFaSb+gBwvbX2TWPMPwA3AjeEUywRCZPm1IpEW9nB2Frbm/HPI8DS4MURkUrRnFqR\n6AprzHgN8FBIxxIREakrBVvGxphe4Kwcv7rJWrs39ZoNwKC19t4KlE9ERKTmTRoZGSn7zcaYdmAt\nMNda62cVgfI/TEREJH4m+XlRkGzq+cBXgDk+AzEAx44dL/cj60Zzc5PqySfVlT+qJ/9UV/6onvxp\nbm7y9bogY8ZbgLcCvcaYZ4wx3wtwLBERkboVJJv6/DALIiIiUq+0HKbEWnrzA/A2T9DcWRGJIwVj\nia1kMsmyZbvo61sJaPMDEYkvrU0tsdXdfSgViBuBxtTmB4dcF0tEpGQKxiIiIo4pGEtstbXNZs6c\nLry9SgZTmx/Mdl0sEZGSacxYYiuRSLB//3K2bNHmByISbwrGEmva/EBEaoG6qUVERBxTMBYREXFM\nwVhERMQxBWMRERHHFIxFREQcUzAWERFxTMFYRETEMQVjERERxxSMRUREHFMwFhERcUzBWERExDEF\nYxEREccUjEVERBxTMBYREXFMwVhERMQxBWMRERHHFIxFREQcUzAWERFxTMFYRETEMQVjERERxxSM\nRUREHFMwFhERcUzBWERExDEFYxEREccUjEVERBxTMBYREXFMwVhERMQxBWMRERHHFIxFREQcUzAW\nERFxTMFYRETEMQVjERERxxSMRUREHFMwFhERcUzBWERExLGTy32jMeZmYDEwArwCtFtrfxNWwURE\nROpFkJbxZmvth6y1U4HdwF+HVCYREZG6UnYwttYez/jnW4HfBi+OiIhI/Sm7mxrAGLMRWAm8AcwM\npUQiIiJ1pmAwNsb0Amfl+NVN1tq91toNwAZjzA3AbcBnK1BGERGRmjZpZGQk8EGMMe8CHrLWfiB4\nkUREROpL2WPGxpjzM/75aeCZ4MURERGpP0HGjL9hjDHAMPBr4AvhFElERKS+hNJNLSIiIuXTClwi\nIiKOKRiLiIg4pmAsIiLiWKBFP8phjPkmsBAYxEv8+qy19rVqlyPqjDHLgK8DfwR8xFr7tNsSRYsx\nZj7wbaAB2Gat3eS4SJFkjNkOLAD+r7X2g67LE1XGmHOBHcDb8dbb/7619jtuSxVNxpgE0AdMASYD\nP7LW3ui2VNFljGkAngKet9Yuyvc6Fy3jA8D7rbUfAn4B6EvM7afAEuCQ64JETeri/kdgPvA+4HJj\nzHvdliqy7sarJylsCLjaWvt+vNUEv6hrKjdrbRK4ILUvwR8DFxhjPu64WFF2FXAU7yEvr6oHY2tt\nr7X2zdQ/jwDnVLsMcWCt/Rdr7S9clyOipgO/stY+Z60dArrx5rpLFmvtj4H/dF2OqLPWvmStfTb1\n//8F/Bx4h9tSRZe19o3U/07G6536D4fFiSxjzDnAp4BtwKRCr616N3WWNcAPHJdB4uedQOZ2nc8D\nMxyVRWqMMebdwJ/gNRYkB2PMScDTwB8CW621Rx0XKapuA74CnFrshRUJxsXWtE69ZgMwaK29txJl\niAM/9SQ5aXK8VIQx5q3AD4GrUi1kySHVuznVGHMa8Igx5s+stU84LlakGGMW4uVqPGOM+bNir69I\nMLbWXlzo98aYdrym+9xKfH5cFKsnyevfgXMz/n0uXutYpGzGmEbgfuAea+1u1+WJA2vta8aYfcA0\n4AnHxYmajwKLjTGfAhLAqcaYHdbaVble7CKbej5es31OKhFAiis41lCHngLOT3UnvgBcBlzutEQS\na8aYScBdwFFr7bddlyfKjDFvA35nrX3VGPN7wMXA3zguVuRYa28CbgIwxswBrssXiMFNNvUW4K1A\nrzHmGWPM9xyUIfKMMUuMMb/By+zcZ4x52HWZosJa+zvgS8AjeFmKO621P3dbqmgyxvwA+GfgPcaY\n3xhjtM1pbh8DVuBlBj+T+k9Z6LmdDRw0xjyLN66+11r7mOMyxUHB4TWtTS0iIuKYVuASERFxTMFY\nRETEMQVjERERxxSMRUREHFMwFhERcUzBWERExDEFYxEREccUjEVERBz7//ai6xD5sTZiAAAAAElF\nTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.scatter(points[:, 0], points[:, 1])\n", "centroids = initialize_centroids(points, 3)\n", "plt.scatter(centroids[:, 0], centroids[:, 1], c='r', s=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's define a function that returns the closest centroid for each point. We will use numpy broadcasting to do this. If this code looks overly complicated, don't worry, we'll try to explain why it works in the next section." ] }, { "cell_type": "code", "execution_count": 107, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def closest_centroid(points, centroids):\n", " \"\"\"returns an array containing the index to the nearest centroid for each point\"\"\"\n", " distances = np.sqrt(((points - centroids[:, np.newaxis])**2).sum(axis=2))\n", " return np.argmin(distances, axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can test the code like so:" ] }, { "cell_type": "code", "execution_count": 108, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0,\n", " 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,\n", " 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,\n", " 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,\n", " 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2,\n", " 1, 0, 2, 2, 2, 2, 0, 0, 2, 2, 0, 0, 0, 0, 2, 2, 0, 0, 0, 2, 2, 0, 0,\n", " 0, 2, 2, 2, 0, 2, 1, 2, 0, 0, 0, 0, 0, 2, 0, 1, 1, 1, 1, 2, 1, 1, 1,\n", " 1, 1, 2, 1, 2, 1, 1, 1, 0, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 2, 1, 1,\n", " 2, 0, 1, 1, 2, 1, 2, 0, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1], dtype=int64)" ] }, "execution_count": 108, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c = initialize_centroids(points, 3)\n", "closest_centroid(points, c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So why does the previous function work? It makes use of what I would call some advanced numpy broadcasting tricks. \n", "\n", "To investigate this broadcasting into more detail, let's look at some intermediate steps in the computing process. First, our variable `c` denotes our centroids that we want to work with. Let's look at their coordinates:" ] }, { "cell_type": "code", "execution_count": 109, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-0.04883724, 0.6316724 ],\n", " [ 0.27231992, 0.32268518],\n", " [-0.76809563, 0.82267701]])" ] }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the interesting things with numpy is that we can extend an array by a new dimension using the `np.newaxis` command like this:" ] }, { "cell_type": "code", "execution_count": 110, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[[-0.04883724, 0.6316724 ]],\n", "\n", " [[ 0.27231992, 0.32268518]],\n", "\n", " [[-0.76809563, 0.82267701]]])" ] }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c_extended = c[: , np.newaxis, :]\n", "c_extended" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our array, which had a shape of 3x2 becomes an array of the following shape:" ] }, { "cell_type": "code", "execution_count": 111, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(3L, 1L, 2L)" ] }, "execution_count": 111, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c_extended.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We've just added a dimension in the middle of the array. This allows us to substract this array from an existing point `p`, due to the fact that numpy applies broadcasting rules to array of non-matching sizes which allow for efficient operations (this is described in detail in the broadcasting section of the following document: [http://nbviewer.ipython.org/url/www.astro.washington.edu/users/vanderplas/Astr599_2014/notebooks/11_EfficientNumpy.ipynb](http://nbviewer.ipython.org/url/www.astro.washington.edu/users/vanderplas/Astr599_2014/notebooks/11_EfficientNumpy.ipynb)).\n", "\n", "In particular, the rule used here for non-matching dimensions is:\n", "\n", "> If the two arrays differ in their number of dimensions, the shape of the array with fewer dimensions is padded with ones on its leading (left) side." ] }, { "cell_type": "code", "execution_count": 112, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 1.13251313, -0.49784654])" ] }, "execution_count": 112, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p = points[0]\n", "p" ] }, { "cell_type": "code", "execution_count": 113, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[[-1.18135037, 1.12951893]],\n", "\n", " [[-0.86019321, 0.82053172]],\n", "\n", " [[-1.90060876, 1.32052355]]])" ] }, "execution_count": 113, "metadata": {}, "output_type": "execute_result" } ], "source": [ "c_extended - p" ] }, { "cell_type": "code", "execution_count": 114, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(3L, 1L, 2L)" ] }, "execution_count": 114, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(c_extended - p).shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The broadcasting allows us to keep the first dimension as the point dimension, which allows us to generalize the trick for more than one point, as demonstrated by the next line:" ] }, { "cell_type": "code", "execution_count": 115, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[[ 1.18135037e+00, -1.12951893e+00],\n", " [ 3.19833710e-01, -5.72897144e-02],\n", " [ 1.66808263e+00, 4.60938179e-01],\n", " [ 1.45933022e+00, -1.33264298e+00]],\n", "\n", " [[ 8.60193206e-01, -8.20531719e-01],\n", " [ -1.32345648e-03, 2.51697498e-01],\n", " [ 1.34692546e+00, 7.69925391e-01],\n", " [ 1.13817306e+00, -1.02365577e+00]],\n", "\n", " [[ 1.90060876e+00, -1.32052355e+00],\n", " [ 1.03909210e+00, -2.48294333e-01],\n", " [ 2.38734102e+00, 2.69933560e-01],\n", " [ 2.17858861e+00, -1.52364760e+00]]])" ] }, "execution_count": 115, "metadata": {}, "output_type": "execute_result" } ], "source": [ "points[:4] - c_extended" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This then allows us to apply square root, square and sum operations that efficiently reduce the number of dimensions of the matrix:" ] }, { "cell_type": "code", "execution_count": 116, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1.63444233, 0.32492417, 1.73059633, 1.97625455],\n", " [ 1.18878284, 0.25170098, 1.55144878, 1.53078707],\n", " [ 2.31432412, 1.06834567, 2.40255307, 2.6585241 ]])" ] }, "execution_count": 116, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sqrt(((points[:4] - c_extended)**2).sum(axis=2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we can obtain the index of the closest centroid using the `np.argmin` function:" ] }, { "cell_type": "code", "execution_count": 117, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([1, 1, 1, 1], dtype=int64)" ] }, "execution_count": 117, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.argmin(np.sqrt(((points[:4] - c_extended)**2).sum(axis=2)), axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The last step in the algorithm is to move the centroids to the mean location associated with it:" ] }, { "cell_type": "code", "execution_count": 118, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def move_centroids(points, closest, centroids):\n", " \"\"\"returns the new centroids assigned from the points closest to them\"\"\"\n", " return np.array([points[closest==k].mean(axis=0) for k in range(centroids.shape[0])])" ] }, { "cell_type": "code", "execution_count": 119, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-0.06971844, 0.66614582],\n", " [ 0.85831935, -0.18948577],\n", " [-0.82669027, 0.23483006]])" ] }, "execution_count": 119, "metadata": {}, "output_type": "execute_result" } ], "source": [ "move_centroids(points, closest_centroid(points, c), c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can visualize these first two steps in the following way:" ] }, { "cell_type": "code", "execution_count": 120, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 120, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFXCAYAAACRLCZbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvX14FPW99/+KMWFrC9JjacWnmvuIY+1ppf1RQmqFoyg+\nAGlyUTV6jKSIv9Jaih4VVFq1pbSFWiz1eOOpiInhaNBqUiAoBGmDPQ1w+Pnwy63tqG18rFqPpwg9\np0vSNPcfs5vdbGZ3Z3dnZ2Z336/r4rrI7uzMZ2fnPZ/vfL+fh7KhoSGEEEII4R9H+G2AEEIIUerI\nGQshhBA+I2cshBBC+IycsRBCCOEzcsZCCCGEz8gZCyGEED5zZC4fNgwjBHQDY4BK4Oemad7shmFC\nCG+RnoXwj5yejE3TDANnm6Y5Gfg0cLZhGF9wxTIhhKdIz0L4R87T1KZp/k/kv5VAOfBfue5TCOEP\n0rMQ/pDTNDWAYRhHAE8Dfw+sM03zhZytEkL4gvQshD+UuVUO0zCMo4HtwE2maf7SlZ0KIXxBehbC\nW3J+Mo5imub7hmF0AlOAX9ptMzQ0NFRWVubWIYUoZnwVSjo9S8tCOMaRUHKNpv4I8FfTNA8YhvEB\n4Dzg20ktKivj3XcP5XJIV5gwYazsCJANssPeDq/JRM/SsuyQHc5tcEKuT8YTgZbIOtMRQKtpmk/m\nuE8hhD9Iz0L4RE7O2DTNXuCzLtkihPAR6VkI/1AFLiGEEMJn5IyFEEIIn5EzFkIIIXxGzlgIIYTw\nGTljIYQQwmfkjIUQQgifkTMWQgghfEbOWAghhPAZOWMhhBDCZ+SMhRBCCJ+RMxZCCCF8Rs5YCCGE\n8Bk5YyGEEMJn5IyFEEIIn5EzFkIIIXxGzlgIIYTwGTljIYQQwmfkjIUQQgifkTMWQgghfEbOWAgh\nhPAZOWMhhBDCZ+SMhRBCCJ+RMxZCCCF8Rs5YCCGE8Bk5YyGEEMJn5IyFEEIIn5EzFkIIIXxGzlgI\nIYTwGTljIYQQwmfkjIUQQgifkTMWQgghfEbOWAghhPAZOWMhhBDCZ+SMhRBCCJ+RMxZCCCF8Rs5Y\nCCGE8Jkj/TZAuEc4HKatbTcADQ3TCYVCPlskhMgW6bm0kDMuEsLhMJde2k5Pz5cBaG+/n02b6iVg\nIQoQ6bn00DR1kdDWtjsi3Aqggp6epuFRtRCisJCeSw85YyGEEMJn5IyLhIaG6dTU3A/0A/3U1DTT\n0DDdb7OEEFkgPZceWjMuEkKhEJs21dPWtgWAhgatLwlRqEjPpUdOztgwjBOBB4CPAkPAT03T/Ikb\nhonMCYVCNDXNSrnN4OAgb731BwAmTjzOC7NEgSA9BwsnehbFQ67T1APAdaZpfhKYBlxjGMYncjdL\nuM3hw4d5YuXtPHXudMqmfYayaZ/hqXOn89gtt3D48GG/zRPBQHoWwidyejI2TfNt4O3I//9sGMZv\ngOOA37hgm3CJw4cPs3X+5SzY1UVl3OuTn++l//leNuzdz5yWBxkzZoxvNgr/kZ6F8A/XArgMwzgZ\n+Ayw1619Cnf4xR3fH+WIo1QCC3Z18Ysf/cBrs0SAkZ4Lg3A4THPzDpqbdxAOh/02R+SAK87YMIwP\nAT8Dlpim+Wc39incYXBwkNBOe0ccpRII7dzB4OCgV2aJACM9FwbRwiBLl9aydGktl17aLodcwJQN\nDQ3ltAPDMCqArcDjpmn+OM3muR1MZMxrr73GnyZN4oz+/pTbPVtZyd+99BInnXSSR5aJNJT5cdAM\n9Cwt+8w993Ty1a/OwioMAtDPunVdLFo020+zxGgcaTnXaOoy4D7gBQeOGIB33z2UyyFdYcKEsSVj\nx3vv/dnRlTAU2fYDH/DnvJTSb+LUDq/JVM9BOU+lasehQ/ZPwaV6PoJqh1Mt5zpNfSZwBXC2YRjP\nRP5dkOM+hYtMnHgcL00y0m738qmGUp2E9FxA2BUGaWqa6bdZIktyjab+FariFWjKy8sJn3se/c/3\nJl037gfC586ivLzcS9NEwJCeC4tkhUEOHRrw2TKRDarAVQKcfcPNbOjttY2o7gc2nHMec2+42Q/T\nhBA5oMIgxYOccQkwZswY5rQ8yKYf/YDQzh2c8qIJWFPT5bVzmfu1f6ayMlW8tRBCiHwiZ1wijBkz\nhgtuuY3BZd8cLod51sTjOPbY8b4HOAghRKkjZ1xilJeXc8IJJ/pthhBCiDgUrCGEEEL4jJyxEEII\n4TNyxkIIIYTPaM1YFAzhcJi2tt2AVfBAzdaFKEyk5dHIGYuCIFoUv6fnywC0t9/Ppk31ErEQBYa0\nbI+mqUVB0Na2OyLeCqCCnp6m4ZG1EKJwkJbtkTMWQgghfEbOuAQpxIbkdkXxGxqm+22WEL5TaHqW\nlu3RmnGJkbhe09nZSmvr3MCv1yQrii9EKVOIepaW7ZEzLjFGrtdAd/cVtLVtKYhi8yqKL8RIClXP\n0vJoNE0thBBC+IyccYmRuF4zY8ZGrdcIUaBIz8WDpqlLjMT1msWLL1czciEKFOm5eJAzLkHi12tC\noZDEK0QBIz0XB5qmFkIIIXxGzlgIIYTwGTljIYQQwmfkjIUQQgifkTMWQgghfEbR1AFHfT+FKB6k\nZ5EMOeMAo76fQhQP0rNIhaapA4z6fgpRPEjPIhV6MhauoOk3IYoH6dl79GQcAJL1I82k76efPU2j\n029Ll9aydGktl17aXhB9VYXIB9KzyAY9GftMqnUkp30//V6LSmzjZk2/Bb+NmxBuIz2LbNGTsc+k\nW0eK1p1tapqVVIxaixIiGEjPIlvkjEuAfE95ZTL9JoTIDem5OCkbGhry8nhD7757yMvj2TJhwliC\nYsfrr78bmZJqAqCmpjnjKanYtNbofSROedXUjJzySjwX2QZu5BrwEaTfJCB2lPltQxqkZRs7UmnR\nKdKzewTBDqdaljMOgB1uRC4m20dz8w6WLq0luv4D/axeHVv/iT8X6YSeT4L2m/iNnLEzAvR7jdCR\n9By838VHGxxpWQFcASC+H6mf+1DghhC5Iz2LbNCacZGj9R8higfpuXjRk3GR4zSdwnpvOu3t949Y\nq2poqPfKVCFEGqTn4kVrxiVsRzgcprNzH4cOhYfXpXJZ7xocHOStt/4AwMSJx1FeXu74s36fiwDa\noTVjBwTo9/LdDrf1nAtBOB9BsUNrxiIlqQoLZLqmdPjwYX5xx/cJ7exi0ksmAE9NMgifex5n33Az\nY8aMcd1+IUQMN/Us/EFrxiWKW4UFDh8+zNb5l3Pp2jVc/Hwvk/v7mdzfz8XP93Lp2jVsnX85hw8f\ndt1+IUQMFQopfOSMRU784o7vs2BXF5U271UCC3Z18Ysf/cBrs4QQOeJnfexSRM7YQ7y+uFMdz42o\nzMHBQUI77R1xlEogtHMHg4ODmZovRKApNj0nHkvNIrxFa8YeYbems2vXlZ4eLz7hPxqV2dnZFQn4\nyLwYwFtv/WF4jTgVp7xo8tZbf+CEE07M/IsIEUAS9dXZ2Upr69y8BUh5oed4lKPsPTk/GRuGscEw\njHcMw+h1w6BixW5Np7n5SU+Pl7iGFAqFWLRodsqi9aJ0kJadk6iv7u4r8rpGKz0XP25MU98PXODC\nfkSBMXHicbw0yUi73cunGkyceFxWx9C6ladIywLIX3ER6Tk5OTtj0zSfAv7kgi1Fjd3F3dQ009Pj\nNTRMd1UM5eXlhM89j/4U2/QD4XNnZZRzHEXrVt4iLTsnUV8zZmzMayUsOz3X1U3Nm2OLTnuvXr2F\n1au3uFLTWnpOjdaMPSJZ5ZxDhwY8Ox7g+rr12TfczIbeXtuI6n5gwznnMfeGm7Pat9atRFBJ1Nfi\nxZfnTct2x6uru5D58x8ftYYMY109pptak55TI2fsIV4n4IdCIRoaptPWtpu2tt0MDPSPEkNzcxfz\n5mU/oh8zZgxzWh5k049+QGjnDk550QroevlUg/C5s5h7w81UVqaKtxaiMInXcz4H1nY88sivbB3b\njTfO88wG4S6eO+MJE9wbueWCn3ZYU8VW8FZT08y8RmBefPEjdHc3AjBp0o+BMLH2axa5n4uxNN55\nB4ODg7z55psAfPr44zOemk60Y/Hi2XR2ttLdfQUAM2ZsZPHiy/MenBKUazToBOU8+W1HOBzmnns6\ngfzp2YmWx46N9TQOAtJzZrhSm9owjJOBLaZpfirNpiVfzzbXHqOZ1Jq1631aVXUrfX3fiRy7mV27\nrvR0RJ+MZL+J17V1g1DLNmKHL7WppeXMyEXPbmt506Z6TjxxQlCuX+k5ZoM3takNw3gImAEcYxjG\n68Ctpmnen+t+g0A+moTnsm6SLtcw8Vh2LFw4mYoKb9at3UC1db2jmLUMuevZ7vPZ6jlfWg460nNy\ncnbGpmle5oYhQSOdWLLdx5w547K2KZXw7Y7V0nLhqBZqjY35eQr3Yj9BOU6xUqxahtz1nOzz2SIt\nB+M4QUIBXElwI/LPbh9z5jxKTY37PUbtjtXRscVx79NE3BiMuLmf+P3ZidTt44jiIlc9J/t8PnoG\nl4qWo/tM1HOpalm1qT2moqJyOH9v3bqujC6ybBLxo9NCmVblybQLTLL8ZTe7yaTKU1TXGuEH0ZSj\ndeu6MsrHDbKWwV7PbmssmZ5LVct6Mo4jfpRWVzc1pxFvOBxmYKCfqqo19PVdN2IfUVGlCy6wGzUm\nGx3nY4TuFK9GsspTFE5J1E6u+qirm8q6daO1DLEylMWg5aitbk7JJyOZnksVOeMIydZpOjpynRYK\nU1V1KwsXTk67xpM4GLBL6k8WAJFK3NmQyQ0hlZP06sbi9w1MBIdkziSXad758x+nr+9rwONUVfXQ\n0nJN2s9H9TwwMMDmze+xd+//O8KeIGoZvJ2Sd8PeYkHOOEKydZpUT112EY/Wfl6gp+daotMsfX3f\noaJiS1pHHH8DsUbhXyOTp0A3IxXduiG4eWNJJVK3b2CicEk1OHSq57q6qXR07ANIKJbzRfr6LnR0\nb4jXM7QAg0BIWo6QTM+lqmU5Y4fYOd54sT366E8pKytnz56rgFos8V0BWBfRwEDq9KHEG4g1HfY4\n8EXXv4tTnN4Q0o1k3bqxpBOp0iaEU1LrOcx3v7uWgweXAXDyyXeQWGAjUz3DlUAXMNvNr+GYTLSR\nbtDrxSChFLUsZxwh1QWYLEUpXmx79x6PJbR48W0DLgIeAD6Y9NjhcJienhewfo7ziTrwqqoe+vou\nHGVP0PByJFuKIhWZkW5wmF7Pj3Pw4GXAW8DxvPLK9cBtwLcje8hOzzBALFgrmFoG6dkv5IwjpLoA\n7aa9jj12lYO9/jbymUuoqOiy3SJ2Y7gp8koLcAk1NW20tFyT1Zq1H8SLKhqJCcHJESzFvMVSJZ0z\nSa7nw3yW25hNG7Oxyrp28gk6OZ+n+QzWky1ko+dp0x5k7tyjqajYEngtQ0zPQdVNUO3KBTnjODIZ\npU2Zcgpvvx0bfVdXv0lZ2X3s2bMAgHHjVnPw4DeAUEbBT3Al9fWrWLt2UV5GjaOn55zXbXUigEwj\nq6P7HDs2xOzZU/NW17cU8xZLmUy1M3nyx/nIL6bxLweeH9F9rJrn+CbPcc3Rp7P+/T3AmJz07Db5\n1HM2upGes0fOOA1WitIAVVXfoq/vFqLOtbGxnsZG4kbflwCxv+vq5tPR0RV5L7MLpabm9LxcWAcO\nHOD88++mr68GOIf29occt1B0KoBM0o9yrdPtFKVECRgZ3Vxd/VP27r0asKaxT/zjC9yU4IijVAJ3\nv/8C489u4qMXXl4yes5UN9JzbsgZ22CfkjCPj398FWec8Veqq08D7Eff8X+7EfzkFuFwmPPPf5C+\nvhWRV1rp6WmguflJ5s2bnnKUHA6HWbLknsjUm3sCKFZRiWARDodpbX2S9eufpa9vORBi2rT7uP32\nh3juudf47Gf/ng88tNPWEUepBKa++xJnNc5M2Y0s6HqOFdV4QXoOGCXjjOMvzsWLk0c0Jk9JgFdf\nPZZXX53P5s2wdWs0D9lKf8h23cKrYIm2tt309f0zsemzK7ACzCpSjpJj733a0XGCmCMYRJtE9jjV\ncnTb2LU9D2gFLmfPnst555119PXdRHv7a+zDTHvcU140eeutP3DCCScm3SbIeob4iPH0t/6g6iao\nduWKKy0UM8CXtmuJF+eMGa20ts61FYldq7JY4MasEa9bLcyskamTKZnE0aqX7c6StWB74YXbueuu\nzlHvrV5tjWhjnxsEHsQSfaxlW7J1YyfBFbHfpSntPnPFqU1BaLkWscOXFooZEHgtQyo9DxDLfniN\nPUyimv6Ux36mshL2PMMJJ5zoq5YhOz0Dca+HgY1YWR/JtZdJoJT0bI9nLRQLgcTpk+7uKzKcPvkL\ndqfKWqvJbj2lvf1+x+s7Tkl1gSaOJquq7mT79vRVhGKEgMuBbdTX96YMSHEaPBP/FGEFfCiFQqQm\nUy3b5wP/Bfj/iOX8Hk8np1HN/5/y2C+fanDWxOM80TJ4oedLqK9fRU3N6Umf4DPRjfScG2oUkUBd\n3VTGjv0B0QLuVk7hb4HDQKyw+8kn/yjy/06sUWZq7IqfNzc/6ZrdqZooQEwoq1dvYfXqLXR3L2D8\n+PFA6qL1I987gpqaP7oaGRoV1aJFs4enxe0aTgiRKeFwmM2b3yNet5aeXwS+SlXVmshrg/z7uA+n\nfC7uB8LnzqK8vDzvWo7ankzPg4OD/Od/vsuaNZ/jBz/ocKzn0a+3sXbtoowbT6QiXs+AtJwBJfFk\nnDiKnDFjIw0Nc2237ejYx6FDp2FVvzoSa1r2CGAb1dUD1NY+BsDPf/4hXnmlLvKpFqqrDw9HVLtB\npnl0ToIn0tXCbW19lP37X2bKlEmj3vOiAECxpiwI98hEy21tuyPBl1tJ1HNV1a1s327l8Q8MDNDe\nXss/7Q/xb2wfFcjVD2w45zzm3nBzVjZnkxNrp+eNGx/jhHf+g9DOLia9ZBIC/tckg/C551FWNmP4\ns6FQiJaWC1m2zKqFsGpV0/AxpeXgUhLOONGhLF58OYcO2ZezGxjoB34PLCV+zSV+ara5eQf79i0i\nPpewtvax4QvtwIEDLFvWDFhCGD9+vG3QQVPTlbZ2uHMhR6sA2UdH2wXAbN16iJ6em2hvtwLUosdM\n5sTdTrxXNKZIRyZatggDJnA98XpeuHAy48ePH46L2L9/Mfu5ht/zHS6ik9ry5zmy/AjMSafyHx/5\ne3439gz+q6WLxsaZPmgZ4H3e/PGtXP/Hl0YMFiY/30v/872seGIXR81fxnXX1Q03togWHnn77fRa\njtrqlp6bm5+UljOkJJwxjHwqDIVCSYWzZcv7wDewoqitdaBp0zawalWsp6bdOlRFhXXRHThwgClT\nWjh40BLCk0+uYv/++YwfP952VGpnRzZOaeQNIsy4cT+hvd1yrPbR0dbNobPTCoDJNacw8RjFVh1H\nBAcnWgZrycmqMb2EeD3X1DRz8cUXDleJi+l5DE+zkqf5JoNLW6itrebOJb9iX/cHgPn8/OewZct9\nPPzwvLxq2drnSD3XVH6JNQmOOEol8C3zWWpueZbt2/+HWbOOyviY0rP/aM04jra23ZFGD2OJpgrU\n1n6PBx6Yzfz5jw+v32ze/B7Tpt2H3RrrsmXNkQLz1nrSwYNLh5+SozcRN9doosSvCdfXrx1hQ3xz\n7sT1LisAJvPG3ckagKdbu05FNg3XhUhGR8e+iA5ieq6vX0VLy4Vp9NzKV786j6ee+j379p0EzCd6\nne/Zs4C2tt151TKM1HNd3Z2c2/9+2jzoi9hOd/dl7N//UsbHc1vPTU0zpeUMkTNOSgi4iC984VN0\ndOwbcaHu3Xs1c+cezerVW1i58jHmzBk7fOHa8fvfv5VREEO2Til6g6ipOd3pl8z5mIkkE7UTEoPM\ntMYk3MPSc03N6Y71nKwzU0/PC54MMKN6/sQnPsZsfpN2+9m8ALzJlCmnuOYIs9WztJw5JTNN7YRk\nyeR2F19FRSUNDdNHTe3cdVc9mzevYHDwm5EtV/Dcczfz3HPH8thj90WKxVemrCOba9BUqqT4ZAEw\nmR4zk3OVCcWYsiD8IVc9f+5zd/OhD/0Xf/7zu0SnuGED7e3X8vbbD40o+pOs+IgbAZC1tdUc8f2/\nOdq2urqdxsamhFK96Y+ZDz1Ly5lREkU/ojgpYh4tnbd372+Bv1Fd/QkuvvgLkYCIJsC6UK1oxWba\n2z9NrE1aP/X1q2hvXwTcF9njAuA/sHIa+7GiOi+ipsbKTUwdfJL7dwX7AK7W1p3s3/8yZ5/9Kerq\nzszoBhErF9oPlAFDQBkVFRXU1U0dda6cjIqDkJwfMDtU9CMFThsS5KbnDmAmsBMrCOwUoByoJ7Ho\nT7riI7kwODjIL2d+gYYXnk+53Y+OPpam3/0f/va3VBPao4nX88DAAM899ypTpkyisXEmQFaFPAKk\nI9/tcKrlwDpjt4MGrKLqD0ZKyCWvmBUOh7nkkkcja8cQTVtqba0dHgXHHE60ZKZVYg+OiDjjWM1X\nywFvA+qIVf+xHPO6dV3Mm+f9OkouBd0TP1td/VPKysqHz1dNTXZlQtOJxqsgkiCIN2JH0Thjv7Qc\nPXb2en6CkT3K+4E1wLVYK3yPA18cfi9atS4fPLHydi5duybpunE/UMMyxs44NaNBwUg9hxk3Lhpv\nEjuvQMa/n/Q8wobCrcDldo6aVVT97sgoNnWEYSyIK5a2tHfvNjo69g1v29y8I6FNmhUcUlPzR1at\nauLNN+9l376FkfcewCq9dxB4GPgSVqGQAQYGPB0IDZNLClHiZ/fuPZ74G1ZPTxMdHe7elJSzWLj4\nqWXITc/V1W8yNJSo5Wuwyki+DXwFS8sAZ2X1fZxy9g03s6G3lwW7umzzoP+JC3ma70A3GaUQjdTz\njrjAz5HnVXrOP4EM4MolCCjZ/qzSlfmjvr53+IJ6553fY42au7CEfRXwY6ycx59h1biezSOP/GfW\nlWkKuVJVpra7fT0I7ygELdsFatXX99LaWktZ2RDwQyw9X4EV53ElcALwr1hansXRR99FXd3UrI7v\nRA9jxoxhTsuDbLr2eh75h0/xTGUlz1RWcs/Ek6lhGT+jA1LGW+eHbO5D0rM9gXTG+eEcrOknK8Kw\nqurOURGGsd7Fy7GeZK3yedXVb47Ytq5uKuPGrRre17hxq4er3LS17ebVV88CLsJ6YoyO9t4D/o74\nNImnnpqf1UWYS/oQ5BbhmfjZ6uo3k6Z55cN2IZxqubl5BwMDA0ydes/wtol6jpXMbBneZtq0Daxd\nu4iOjn2RCl6TsfQc/+S2C4ilD77//o3D096ZkIkexowZwwW33MZZXbthzzOw5xku+HUPY2r+V2SL\n/khApvOlr5F6PmvEfU1a9pZATlO73SLL2t9D9PQ0ANs45ZR9bNv21VFBTfGt1k466fscffSbVFVN\n5I47rhqxrZW/uIRoN6eDB79BR0dX3FTOOVhPwksif28EVgD/DFyc9feIkmulqlwKuo+ODr0kYpOz\nyM1ktt9447yknynWlmmlgP9ahqlT72H27BW88cZ71NdXs2DBF4e3j5XMHCTazWnu3KMTruGZjNZz\nbdbfIZ5stFxeXj6ijWNmFclGkqjnurr5dHTkpuXMCppIz1EC6Yzdrocc25/lPBcvvnX4gh3ZbPta\nou3FXnvtOOA6ent30dt7N5s3X8Xjj1tdXaxprRCxri+xEvOxm8XHidXDtYK7YBYnn3wHr7xyPRDm\nlFN+wMDApwmHw56vl0TTDrIJcLBLWchnCoOX9bGFu3ipZbD0vGTJPZH+21Yf8n37PoBV3nYXBw/2\nMG/e52lri9cyxPTcT0VF1NZ4p5Go52gw13UERcvR/2eaoZGo53ynI0nP9gQ2mjqfRB1Q4gjamqq6\nAngSmA48AjQCEArdTjh8ExCKiyBeAIwO94+mDq1Z8xveey+ab7wR+BIrV1oBH+vXvzYcDVpVtYbt\n2y8f7roS3Ueq1CS3+oZ6HW2YzHav+8EmIwjRlxE7iiaaOp/E/16j9dwKfBg4l1y0HE37uffeVyID\nabD0PIfa2n+ht/eDgdBy4vnIN6lsD5COfLej4FOb8kn0B7JvPL4N+GvktYsS3oulJa1c+SgVFVbA\nRHRdJVFwVgrG3fT1TQPOoaamLTIi3G3bGLy7+3rb+tF2qRu5pgbE52nOnPnpjFORcsHO9iCIBoIh\n3ogdcsYOiP+97PX8Q+AfyFXLgK2e58wZx/Ll8/BTy/H7CIXKOXToL8NFTPzQMgRKR77bUdCpTX5S\nW/sszz13JK+++t9YAranoqJyeDonWaj++PHj6e6+PnKxdqWcjunrm0Zb226ammbl1A7RCYn2jhu3\nKrIGHvIkzUCVeYRXzJ79F3bu3M7hw7lpORQK2erZLgDTSy3b2WzN8M2hvf0habmAKKFo6tHYRRVX\nVxu8+uoS4DPAd4mPmLZyCUdHGVqCuwzYAeygp6dhWKSJBeWtiO1+jjlmBbEIz41YQV/ekJhacPDg\nUmA3SjMQhYydnj/3uUkcPnwW8VoOhb5NNlqG0Xquq5vqq5ZjNsf0bKVe7ZaWC4ySfjK2CyRobd0J\nbMJKQboQuI3a2jHcccdVdHR0DW8XP9q0ykJGPwPQwsDAB0cdb+QI9gCh0HWEw+cBX6Kqah0DAycR\nDocVbShEFoyODL6Q889/ECuLIQzcDnyWm24yOOqo3LQMDPcNfu+9xcBtQDXSssiWknbGYDfNUkZ8\nLjB8m+rqx4YbkdsT/xmwRqaPjdhidJTnBMLh71Nbu4be3j309S1n+fIQW7da02L5jDZMvEGMG7ea\ngwe/QexJQTcLUZjE67m5eUcksCqq5dupqrqVBQuuT6EnZ1oemYExFssZP8Exx1xHX9+PPdOytc+R\nerYqhV0iLRcYJe+ME6moqHD0mpPPxEdibtnyPnv23BR5N1r7NkR5efmo0n5Llqxi7dpFNDXN4p13\n3mbu3NsA2LjxOj72sWNz+HYW0SeI1tbH6O19hdNOO4mKimhZv3G0te1WE3FRlCxcODnldZ1K/9Es\niVgmRC2xDAyrReN775UR7cbmhZZhpJ6fe+53hMOH+exn26mokJYLiZKOprYjm1QDu89EG5hbU9Lb\nGB3NuY140itvAAAgAElEQVQZM/6LWbM+wPLlX0p473Fqat7hrrs+T3X1wwwOfguA8vIVPPtsoysi\nTgz6mDbtPoaGBiMFEFI30shHgfcgRD0GzA5FUzvACy1HmyVYr38MOy3DRVRV3Ulf39eItUb1Rssj\n7R4dlOm1liFQOvLdDqU2pSDZDxR7kh0AhjJKD0i8sEemL3Vi1bC1CorAE5xxxg5++cvvcejQwAjx\nWwEgVpGQ449fwJtv3k+88M8441q6ulbl8vWBZGkgqbvQ5NLtKR1BEE3A7JAzdoAXWg6FQnF62YGd\nlq+66jwuuGDKiNaMXmkZkuk5lr7lpZYhUDry3Q6nWi7paOp44uusLl8+j61bD7k4WpyJNZ11EEug\ns3nuubXU1VlrUZs21VNfvwrLGVrT10FEBd5FIZB/LbeSqOVHH/3T8HSxtCyyoaSdcXzHkdbWnVlf\nnHYF0+vqpsalWRxBdfVhamvXEB8c1t19RSRRP8TatYuoqXkH6yexAql+/vObKC+PpU2Ul3+XjRuv\ny/o7xhdxT0wDmTZtA9XVb2DdZDqoqro16y40QvhB9Fq3AiUvwy0tR6OiLb0cAXyJY465jqBoGUbr\nOZaKeZCqqlsZGBhQE4eAU7IBXInTNFVVa7CmnVIHa9lhl9jf0bFlVEOFtrbdbN5sv49k9VqffbaR\nK664FrCCPo4+ejzNzTsi26Qe7YfDYS6++GH27j0BgEcffZhHHrmEUCg0fLzOzi4OHQrT0DAv0it2\nHX19/0xf32zmzx85daU0DRFURuo5MbDKOamKdMTrc2DgPJYvt9+H11qOP2ZUz3V183nkkc5IsNkK\nli9nOLo7FApJywGkZJ1xouj6+q6jqupW+vq+A+Tn4kwUgNXubO7w+3bVbD72sWOH15Uybcrd2rqT\nvXvHEK0ktndvC62tO7n66jnDx1u0aPbwmorVKzaaCjK6WpAKvIugkqhnKyXJCqzKh5YvvvhMtm4N\njpajx4zXc0VFZVI9S8vBo2SdsR0LF06O69ji/OK0G2XW1V1oK7Zc2p1l2rJs//6XgZuIv0Ht37+K\nq692fMhRqPydKBTq63upqRnMWcsNDfW2zrOl5cLhdoPSssiVknXGdqJrbHQu2sSIy9GVvJ6MiG0Q\neJKeno8Oj2RzaXcWOTpWZ6mBSMUge/v6+w8Dq4FPAhcARzBlyqSke011I8pXCoQQbmB37a5duygr\nPcc72egTtbUOfS1WNDX09DSM6GEeNC1btts/JDidGhfeknNqk2EYF2B13i4H1pummSpWP1DpENk6\nmXRpAeFwmBkzfkRf33LiW7dVVa2hu3tBys4m6dqtWetGY4iW65s27T4efnhe2u3gfqZOPczPftYw\nYttEGxKPD+Q1BSKZHX4RIDs8T20qZC1D9nmzqfQce+/DwPvEl8lcufKDw9PE6VKsEm3Kh5bt7Ig/\nfl3d1LjaB/nTsp0dfhEEOzxJbTIMoxz4F6yh2unAZYZhfCKXfXpJYtF3SB2xGCVdWoC19roc+D6W\nI7a26+u7LmVUZ7JIznh7a2uPIT6Kc8+eBaP22da2OxLoEV/Ws4kvfnFCWuFFz0k0VzqXyFRROBS6\nliE/eh75XryersQqnZmcVHr2Qsvx56ShYTrLljVHipYMjvqewn9yTW2aCrxsmuYrpmkOAG3EqkYU\nHOmcYWaEsE7PiCPQ0/NC0huDk9y/dKU5U+H0s/Hnob39JqzC+UqLKHKKSsvgtp4TtRNm//6Xchq0\ne6FliJ0HS8sXAQ8iPQePXJ3x8cDrcX+/EXmtIHGaCG/Xqi2+DVvs/XOwUiz6gYOMG/cT2ttvyunG\nkO7Y0W2snOEW4vOIE7dzeh6sp4Ankh5PFAVFpWVwR8+x986iELVsdx6slK8npOWAkWsAV8YLzhMm\njE2/kQfY2TF27Ohpn1ConEcftQTc1DQzMjU0ll27ruTeezv59a9NPv/505gwYWzctJH1fnPzkwwM\nHANs49e//i1tbcuIj5788pfXcP/91w5/bvHi2XR2ttLdfQVgpUssXnx5wnRUdN9dEZuutJmuGssv\nf3kV9977OL/+9Ro+//nTuPrqpqTTWonnwu48NDT8lhkzKpIczx2CfG2UAEWlZXBLz9lpGZzo2X0t\n252PZHpOtNdtgn59BI2cArgMw5gG3G6a5gWRv28G/pYi8CNwQR/xJBaJnzZtQ9LmCZnWdk1WC7qm\n5p1RwV9eRS6Hw2E6O/dFin6MDC6JnYcwVVXfY+HCyTQ2znTVnvjvunjx7CyjUd0lCAEfETs8DeAq\nNi1D/vTsVMtRG6Rn/wiCnj1pFGEYxpGAiVWw9Q/APuAy0zR/k+QjBSHg6AU1MNA/qqNStOC6nSAT\ni7En7jdZQ4hUn8sXTiLCW1t3cu+9r/DKKzcA9tGebh1/xoxWWlvn+p5qEQTxRuzw2hkXnZYhP3oO\nmpZH2pRczxs2bOOuu37He+99y3YbN48vPY+wIf/R1KZp/hX4OrAdeAHYlEK8BUF8RGZFRaWr+w1S\nEfl062mWiMoijjgW7dna+mRejh+t7ZsLTiJnhT3FqGXIj56DpmVwtj7e0vKbiCN2PzPCbT2XopZz\nbhRhmubjpmkapmmeYprm990wKig4C+xIHnyRSLIi8qk+l+lF6eZFvH//S45eCwLuRs6WJsWsZXBX\nz4WmZSvdsianfXhFqWpZ/YzTYFcEIz6JvqNj3/B7mVb7GTs2xOzZU5N+LtN16Uy2d9J4/d57t7J8\n+X9jRVMDPDCi0IGT7wj25ybx+DNmbMxpWivTZYNkBGFaK2KH+hk7INPfy209B0HLI7dvimw/Us+W\nPmZhFSGyAsqqqu4cUYTIyfeE/OvZLS1DMPTsVMslWw7TKfH1WzMt7h5P4sXc1DQr7YWSrn5t4j7t\ntm9tfWw4JzF682lt3cn+/S9z/vkf57LLthEO29fvbWw8l82bH2bv3m0AVFcfprExfeqpk/OUWKg+\n09q+QmSDG3r2S8ttbVuG34tuB870bJXGfIiengZgG1VVe9i+/RrH9y7pOf/IGWdApsXdoyS7mCH7\nkHu7fc6ZM27UdmvWdPHeez8GYN26FQwODvDaa58AbqK9Hc46q5l/+7c6W1GGQiEeeeSSOPFf4ki8\nTs9T/I0x+9q+RGxTSziRGdno2U8t/+pXvaxb9x+R6n6Z6TnmLK00qoaG6x0/tXqt51LVcs5rxkEk\naIv/TosPJJK4jjVt2gYGBvppbt4R14gitk8YGrE9PBBxxG3Ag/T1reS11/4RiH5ukKeemsCSJfck\nPU92JQaDRPS3bmvbTUvLhaxevYXVq7fkreau8BZp2do+FFrJ5s3X09e3gmz1HHQtQ2yGYM6csaxc\n+VhJabnonowPHDjA+ec/GOnjmdlUcjzRi2JgYAAYoqKikrq6qZ6O2OKnfgYGBti8eTCSmgFVVd8C\n5o3YvqKikk2b6vn617/H5s2fIdZcfT5W1GcFsZ88jFUWr5H29tm8/XZu5wlia0lejWxzWTYQwSef\nWm5omO7pE1i2Wm5tfSxudiv69F0aes5nI4sgUlQBXLFuSStItfjvJDcx/qKwStFdSk3NQ5H2asmD\nPOwu5mTBFSeeOMFxcMHooIaDVFWto6/vuhH7DIVCfOUrP47UoY0vSvBz4KjI/98FPgrMTnme0pGu\n200mxQ6yCbRwM9AjFzvyQakHcHmhZWt6mZTXafIuZk1A/rVsbXskVk3p2PZwFzAZOBP4GfARpGd3\n7HCbkgzgcit8P3GNxIom7qKnp4mODieFPUY/qbW0XMiyZVYxo1WrUpezc0aIhQtPoqIi1nc1us8p\nU06hvb2FWBT07Vg/9bcif3+XiRN38dZbs3OyINVaUvz6kRCZ4oWWo9dqpnr2UssW5wCtWDNdYeBH\nwDcj77Xw4Q+/yAkn7Ke3Nzc9j5wul569pgjXjKMXrrXWUlV1p2fF0JOtJ4XDYebPf5z29ptob7+J\n+fMfz3j9yy4PsrHxXNs1oMbGc/nc5/4M/BD4ISeddAjLEUcLxS/nK1+5KOM86USsab/0r+WLbHK9\nRSHhn5bBXs+trU96qmWrUcQDwDjgh3z4w1/HcsSxJi5LlnyKzs4VOWkhHA6zfv2zGX0Ptyl1PRfV\nk3Gy8H2wpkCi26SLfExcI4EHgEuyXitJ9gR5440j14lSTQUlpg7YpSLFU1ZWDvxD5K9BrBF1rO3a\nUUd9kE2bLqSzsytSyzabtZkhrGm/S4BdwF4GBk7LcB/Zk+k5EYVDqlScxBrI6ffjnpb373+Jnp7Y\nEpB3Wj4fqKCs7G3gADBh+P2jjjpqeJ/Z6jnWg70V+BKwi2OO2UVd3Y2O95Erpa7nolozhlRrPLF1\nkF27rkwbdp8s6CPVxZFsbbitbfeotZD6+lU8+OCNw3aEw2EuvvjhSCNxqK5+g0cecZZKlEisWMf8\nyCstfPzjv+XVV783wq5QKJTTmoq1xjMdaAauBaCqao3jQgLxBGFtJ2B2lPSaMaSLv3BeAzkbLUc/\nl6jnOXPGjqpvXV+/ilmzPjtc9CP2ucuAXVRV9bB9+zWMHz8+43MwslhHIwCh0O2Ew8uA0KjiHtle\nv7H12v8G7iWq52yDqAKkI9/t8KRRRBZ4XrXHLihg3bou5s3Lz/SHkwAua3R+KTNmPDJ8I7FzoE6r\nXSViF8BVW/s9vvCFT42wC3K7WJ0G2TjBC9E4CUIJgngjdpS8M7YjH0E+qUgXwDVu3GoOHlyC5Rjv\njxt8j3Sg2Q5S7QO4rAFATc3po67jbK/f2D0q98DOXOxwitOAsiDo2ZNGEcVMtvmNdrl8owvLXwGM\npbv7Clpbd9LcvIOHH34KuJT4taD9+1/OyvYpUyaNeq26+jTXcwxDoRALF052ZV/5plTr3YrccpUT\n9RzV8urVW6ivXxVxxGOx1pQbWLLkHnp6XgCewHLElp77+q7LqnFCQ8N0qqp6Rr1eU3O661q27lG9\nruwvnxSrloveGdsFBTQ1zUz5mXz82KFQiJqa07FGuFEBhVm//jWWLq3luefWApuw1nYtzjjj5KyO\n1dg4k2nT7iO+wEBjY+rvnC2NjTMLIugi22INIlgk6nnGjI1pmzPkQ8tNTbMieo5pGR4eDuwKhZ4l\nXsuQXXBjKBRi+/ZrqKpaQ741Fmt+EWw9F6uWiyqAy45kQQGp1oyzLXuZjsRgklNO+QEvv/xtRqZd\nbMNy2M1UVIwMNHM6NRMKhXj44Xlx39m9HsStrU+yf/9LTJlyCo2N55Z80IXwlkxrIOdLy5Co5yew\nlpis44TDtwE3AHdEtn4A+ODwZzPJ2x0/fjzd3Qtc1VgyLUvP/lH0zhhG1kz1g3jhWUVDtkTsmso3\nvpG49W+xBN1ARUXXiH1kUm0qm++c6gYRDoe55JJH2bPnqsjxW9i8+eHhILOg5yCWar3bYsTNmubZ\nYKfnp59+iba2xMjuGUBUw5cM6zmbynGZaiwXLQddz8Wq5aIP4LLDedWeJmB0O7JMRrWpqtqMHVvB\nOec8kBAM8g3soiTzFbgSPRfpStHZHR+2sXr1oCvCVQDXKDsUwOWAXLUc3SZXPU+YMHaElqdN28DQ\n0CB791496rj5DEKbMGEsr7/+rq9ajtqhAK5hG0qvApdbpJqqyXRUm65KVfxx6urm09ER7ari7fRQ\nzM5B4El6ej5Ka+vOrKK5g0rQR/zCfdJNu7ql5xtvnJdwnHmR7b2f7pWWC5OiD+DKhPioS8A2+tjt\n4IH4aM3x48cnjXj2pjpNtNj8LGA269e/Nhzs0tAwfURQGDxAdfWbrtoQtA49onBxomVwV892kdfJ\nqmol0/Lg4CBvvPE6b7zxOoODg1nZYeGvlkF6zhQ9GUfIVwegVOsb0Ys1ul26NaN8BlY0NEzn7rtX\n8OqrK4mO+q10jNhT/MMPz6O19bG4oI8vumaDOjAJt8jntZROz04DLBO1XFZWxhMrbye0s4tJL5kA\nPDXJIHzueZx9w82MGTMmIxvXrVtBX58/WgbpORu0ZhzB6TqOkzWoRJIVAmls3EJ3d2NkP/5crNE1\npg0bHmfFim4GB+8hH2tZg4ODvPXWHwCYOPE4ysvLR7z/6KO7+epXZ+Xl2JkQhDWmiB1aM3ZALloG\n9/Qci//IvP3f4cOH2Tr/chbs6qIy4b1+YMM55zGn5UFHDnnChLG89NLrVFd/kz/9aT3x52Dlyse4\n+urcmkk4RXoeYYPWjPNBNk+odusbbW27I444s5SLTFuZpSN2M/oY8GNi3WHgmGO+S0PD4pyOefjw\nYX5xx/ddGfEL4TZu6bm5OXnHo2REdfXHx/+Nlb8Y7YgBKoEFu7rY9KMfcMEtt6X9PuFwmPPPf5A/\n/WkOVt34aOc2K7XK7fuHcA854wiZhMtnGzxg5fbtZP/+lyPrQfEjR2efd3PqJxwO841v3E1Pzyew\nmj6EgMux0jEGOPNMa7tsj5lsxD/5+V76n+9lQ2/v8Ii/qWkmDzxQfOkKwnsyTX3JNg0wquUpUyYx\ndmxmg8qYlq/kW6ywdcRRKoHQzh0MLvvmqBml+P21te3m6adfoq/vGmA3cAHxqVWwzbOpY+k5cwp+\nmjqbkV6yqYt8jhpjjSDGEK0/PW7cqqSpTHa4mRIxuun6Bqx4vqhtq9m/fz4dHfuyPuYTK2/n0rVr\nkt5o+oFN117PBbfcNjxd7veoPQjTWhE7Sm6aupC1fOaZG/jrX/9qm8pkR0zLb7GHSVTTn/KYz1RW\nwp5nOOGEE23tGanlFqzOS48RneWaNm0Dc+ceParJRS73j1TnV3oeYUPxT1O7/aSYz3D5trbdkY5M\nsYLvBw8ujSv4nt5uN/sHj266/mWgg2OOuYEzz/w77rjjqrguM2FgR+T/Zzna/+DgIKGd9lNvUeJH\n/FCc6QrCGYWu5X//9y+zcuVjzJvnbLrbzb7fo7UcreT3JY455gYWL/4cCxbMizjGzLUcT7Ry1/r1\nz0ZaLoaS/lbSc2YUdGpTMdQozazge7R/cCwlwXrNHerrTZ555lusX3/9sCOuq5vKuHFrsabUZzFu\n3E+oq5uadl9vvfUHTnnRTLvdKS+aw4FdqVCaRHFTDFquqKjIoBlLVMsT6CR9D/CXTzWYOPE4x7bU\n1/eyenUXzzzzLb72tTpCoVDWWo4SHTAtXz4v0qntEWAwq99Keh5NQTvjQqKhYTrV1W8Q70zPOqsl\no9y+iopKrM5OXZF/l0Rey86exFzHtWsXjbqJdHTs4+DBZURvkgcPLqWjY5+jYwz+7W9Z2ZZIsXZp\nEYWJu1reRSdGyknqfiB87qyk68XJtJw4KMhFyzB6wGRNgT/p+PNRpGd7CtoZpyuE4cXoy+kxQqEQ\njzxyCStXfpD6+lWsXPkYO3Y0ZjQNZ33fh4DzgPOoqWmzvQE4sSkaRbpuXRerV29xPZBj4sTj2P/R\nE9Ju52TEXwxPTSI1ToraBEXPbmv5ae7j6+M/aeuQ+4EVxmTe+OiUtFpevXoL69Z1eZwiOZBxASLp\n2Z6CXjN2s2xlNmTTvOHqq+dw9dWxvzMtcj9nzliOPXYVU6ZMorFx9LEysSkUCrFo0eyUAQ6ZRqbG\nB3Z8uHY2/f96d8oArlQjflE6ZFu20uol7A6ZaicXLYdCIVpaLmTZslUALFnRwab19xDauWN4ecec\ndCoPHTiRzeYjcEslW7aktqepaVbKgKVctNzQMH3U56uq7mThwpNs70Micwo+mjoZqSKP3Yqwyza6\nOXqRjx0bYvbsqbYO1a5ISKri79na5ORcxNtTVzd1eGrLrhtMvI3V1T/ln8b8jK/s/mXSYgZzH3iI\nysrKlHZkU5ghW4IQfRmxo+SiqVOR7Lq+8cZ5rv1euei5s3Mfhw6FbaOGM9FzRUXFcAzFjh3/h5tu\nqs/IHieNM+JtAWwjnpPZl2z7TOwoNT2XRDS1m+SjmEYmF3my99etW8PChScB5K0vazqio+50Tw6J\nUZ17917NF793DJs++/+MGPG/fKpB+NxZzL3hZior069557sUqCgu8pHW5FTPiZqw0/P27ZfT0bEv\nqZ6j6UtHHPGbnO1OJD7COZXtqRrc5HrPkZ7tKeg141Rk0lgh24CCZMdItb906yWJ7/f1Xcfy5Uex\nfv1rWGkJ7n3vTMlmrefIIyu44JbbOKtrN+x5BvY8w1ldu7ngltscOeIoyYrui9LA6XWdS3CQV3o+\n//y7GRhInVecyXfOFj/XbqXn0RTtk3Emo69Uo8BsjtHcvMPlp9gj6eu7jqqqW+nr+w6QfL3Hz1Fn\nqjWp8vLyUQUL4p82Fi/2pmauKEycXtfZajnVMdzWc1/fNCBMTU3q9dugajkZ0nNuFK0zBm+SzjM9\nRrqLPPF92IhVohIWLpxMRUV6YebyvVNN8aWzPZObR+IUWWdnK62tczVKFkkpVD2vW7eGvr7rIq9s\nBL5ERUWXI63k+p2z1XOmAwHpOXeKNoArFYmL+m4HFIzcX5iqqu+xcOFkGhtnDgdvpAvgilW5uQWn\n5TKzIf5cOAkSSxXMBc6CO8Dd0p5uEoSAj4gdCuByQL61PHqf9npOFcB14MABzj//7sgT8TnU1LR5\nErCUiZ6tafMyKioqRgSZSc+u2OBIyyXjjBOnUBLTEPIRwNXaupP161+jr++fgdFiiF4oqYJD8l3b\nNf5iza71XCxyuqysnD17rrL9rolIvGntkDNOQaoBbb4CuFLpOdEJ2kVPe1GnORs92zntlpYLmT//\nccctIaXnlDYomjqKkykUt6fAQqEQFRWVEeEmX2tKFdGYyia/W6GNjpw+HpiN03W1xCmyGTM20tAw\nN/+Gi4In3RNfPqaz3dBzoWi5p6eJZctW0dNzE9Kzd5SEM0682Lq7r/AsLSgZ4XCYe+7ZzY4dT9PT\ncy3pLvrE6eH4UatbBU2yCdrIlsQ1qcWLL8+4AIooTXIJ0soHUS0fOhRmYKA/rW1eaBmk50KjJJyx\nW2Q6nZxMDCNHz7OwatxegdVP2P64I3MVvxUp1F4BhOnp+ShLltxjW1s6EzIJ2kj8btXVb1JWdh97\n9iwY8V3THS96k8qmGpkQueCGnuvqLhyhzaqqNVgpiPZ9ylNrGXp6GliyZG2kk1tuT8lO9Wz3vVat\nauLttzNz5NJzbpTEmnFiUMeMGRszjvRLVZEmXRGPRGHbra9YLc8usg04Gb19B9aU8CDwINBoe2wn\n5LKm4rSaT77tcJMA2aE14yS4EaTllp7b2naP0nJiCmL855Nr2RpYW9HW822P7YRsr1+317kDpCPf\n7dCacRxuTKEkmxqz/p98Wsrp+lV9fS81NYMOcwnPoapqDX19BpYj9me6zu67+R2wIYqfeD1bAVyZ\nT+vmU89OUxAtolq+DngCyxH7V2Uv3WsifxRtBa5EglTxxWnLs+Tbt7F9++XU1/dGtggDncA2R5V9\nhCh0onpetGi2r3q203Jj48yk95pkWl69eov0XOKUxDR1ItlMXSSbGgOymjJLl5tot73dNNLFFz/M\n3r1jiE5tTZt2Hw8/PM/xDSrZufA6wjMI00kBs0PT1A7IZVrWLT27oeXo69Jz8dmR9zxjwzAuBm4H\nTgM+Z5rm0w4+5ruAMxVO4mftLugDBw6wbFkzAKtWNTF+/HhH+3PjQrn33k6WL59Htvl9djY47RDl\nJkEQTcDs8NQZZ6HngtZy9PN2A9zW1ifZv/8lpkw5hcbGcx3t163rRnouPju8WDPuBeqBf81hH56S\na49juzWUcDgcSU24CYC3387/hR5PRYV91GYuBC11RHhCQenZjX7liXpO3Ofbb99PY6O7dqdDei5d\nsl4zNk3zt6ZpvuimMfkmH11KvO58Eg6HaW7eQXPzDsLhcE6dXaz8yM7hfYnSpdD0XAxaBulZxCiJ\naGqv6el5AXB/bSbZ00A2nV1SPVl4WSxAiCCTLy2D9CwSGBoaSvrv1FNP7Tr11FN7bf7NjdvmF6ee\neupnU+0n7l/e+ctf/jK0bt3WoXXrtg795S9/GfXejBn3DcHhITg8NGPGfaO2yeZ48fs8+ujvDsHB\nIeh3Zf/xrFu3dQj6h2Ao8u/w0Lp1W/Oyr1TnUXiCEz1l9M9lPeedYtby0JD0XEI40mfKJ2PTNM9z\n2/nnczE9cXT4wAOj15FaW+fS2dkVCfqYy6FDAzlXimltnUtb2xZ6el6gvf1aYCxgld28667kazOZ\nBhccOjR66unQoXBW59TJvubNmx55PfdzlI4gBFoEzQ63cVvP0nKMbK4b6Tn/BMEOp1p2K884EGkY\nTtZ8QqEQixbNdjXfOBoIUlNzOslKWrpBLutJ+dyXKDp813OxaxmkZzGSrNeMDcOoB34CfAToNAzj\nGdM0L3TNsgIk32szmTb8drKv2JOFdxHgInhIzyPxYp1VehbxFFXRD6c1a/M5dZFJcn1QplD8tkF2\n2Nrh+9NpGqRlj+zIBNkRPDvyXvQjS/JeKMCJgILwAwXBjlyLJriJ03OR70pCfv8mcXaUtDMGaTlT\nCk3PXlQFC8LvUrKNIoJW3NzvxuHJcKNogtcUos0ie6Rl5xSaNgrNXi8omUYR2ZCYkJ/N5y+9tJ2l\nS2tZurSWSy9tD0wyvh8FDnKlEG0WwaCYtQyFp41Cs9cL5IyT4Ib4SvWCy/XGJ4SbSMvZIy17h5xx\nEopdfKlSIXIRYD6fIJS+IbKh2LUMybUhLRcORbdmHCSCXIYuWSpErms5+SxK72YqiBCZEGQtg72e\nAWm5gJAzToIb4gv6BRctmhAfbRj0Di9BC+oRwacUtAyj9dzcvENaLiDkjJPglvhK7YIL+hOEKD2k\n5eyQlr2l6PKMnRCE3LOg2JFog9NiC6nIJgUkCOciYHaUfJ6xEwL0ewXODr+0nGiHnwTBjpLNMxa5\n4cZTRKk9QQgRRKTlwqKknXGQk/j9RAIUhYj0PBppuXAoWWesCjBCFA/Ssyh0SjbPuBRyDwsVFRoQ\nmSI9BxNp2Tkl+2QsgomecIQoDqTlzCjZJ+NMKsBodJcd0fN2772d3Hvv1lHnz+686glHZIP0nF/i\nz9YCwpwAAAmvSURBVNmBAwdsz1/ieZWWM6Nkn4ydRhpqdJcdiecNWoA5tLc/xKZNo6sDdXa20to6\n1x9jRcEjPeePxHP23e+u4uDBJUBo+PzBaD3PmvUBv0wuSEr2yRhikYZNTbOSilGju+xIPG9wJbB7\n+Pwlvt/dfQVtbbtVs1ZkjfScHxLP2cGDS4HdxJ8/Oz1DmbScASX7ZCyCSSGUHRRCpKeiokJazoCS\nfjJ2gp7UsiPxvMEDwFnD5y/x/RkzNg6fVydPOEJkg/ScOYnnbNy41cBZxJ+/ZHqWlp2jcpgOyFcx\ngYCUasubDdHzNjAwAAxRUVE54vzFn9fFi2dz6NBAXuzIhCD8JhE7VA7TAdn8XvnQc4Cum7zYEX/O\n6uqm0tGxD0B6dmaDIy3LGZe4HUGwQXbY2iFn7IAA/V6yQ3Yks8GRljVNLYQQQviMnLEQQgjhM3LG\nQgghhM/IGQshhBA+I2cshBBC+IycsRBCCOEzcsZCCCGEz8gZCyGEED4jZyyEEEL4jBpFCN9ILJ8n\nhChcpOfckDMWvpDYIzXaz1jF5IUoPKTn3NE0tfCFZP2MhRCFh/ScO3LGQgghhM/IGZcI4XCY5uYd\nNDfvIBwO+21Oyn7GQojUSM/Fh9aMS4DE9Zz29vvZtKne1/WcUCjEpk31tLVtAWDx4ssD0f9UiKAj\nPRcnejIuARLXc3p6mgKxnhMKhWhqmkVT0ywFegjhEOm5OJEzFkIIIXxGzrgESFzPqalpHl7PCdra\nkxAiNcn0HA6HueeeTmm5QNGacQmQuJ7T0GCtL4XDYS6++BG6uxuBYKw9CSFSY6dnIHDryCIz5IxL\nhOh6TjxtbbsjjrgCILL2tGXUdkKIYJGo5+bmHXHryNJyIaJpaiGEEMJnsn4yNgzjh8AcrIWL3wFf\nNk3zfbcME/mnoWE6nZ2tdHdfARBZe6r32SrhB9JzYdPQMJ329vvp6WkCpOVCJJcn4x3AJ03TPAN4\nEbjZHZOEV4RCIZ544nJWr97C6tVbtMZU2kjPBUx0HXndui5puUDJ+snYNM2uuD/3AvNyN0d4jd1a\nsig9pOfCJxQKsWjRbN5995DfpogscGvNeAGwzaV9CSH8RXoWwmNSPhkbhtEFHGvz1i2maW6JbLMc\n6DdN88E82CeEcAnpWYjgUjY0NJT1hw3DaAKuBmaapukkyzz7gwlRWpR5fcAM9SwtC+EMR1rOJZr6\nAuBGYIZDRwwQiPWMCRPGyo4A2SA77O3wkmz0HJTzJDtkR5DtcKrlXNaM7wI+BHQZhvGMYRj/O4d9\nCSH8RXoWwkdyiaae5KYhQgj/kJ6F8BeVw/SBcDhMW9tuxo4NMXv2VOUDClHAWA0adnPoUJiGhunS\ns8gKOWOPSWwMXlOjgu5CFCqJelaDBpEtqk3tMUFtDC6EyBzpWbiFnLEQQgjhM3LGHpOsMbgQovCQ\nnoVbaM3YY+Ibg1sBXFpfEqJQieq5s7MrEsAlPYvskDP2gWhzhiAkpAshckMNGoQbaJpaCCGE8Bk5\nYyGEEMJn5IyFEEIIn5EzFkIIIXxGzlgIIYTwGTljIYQQwmfkjIUQQgifkTMWQgghfEbOWAghhPAZ\nOWMhhBDCZ+SMhRBCCJ+RMxZCCCF8Rs5YCCGE8Bk5YyGEEMJn5IyFEEIIn5EzFkIIIXxGzlgIIYTw\nGTljIYQQwmfkjIUQQgifkTMWQgghfEbOWAghhPAZOWMhhBDCZ+SMhRBCCJ+RMxZCCCF8Rs5YCCGE\n8Bk5YyGEEMJn5IyFEEIIn5EzFkIIIXxGzlgIIYTwGTljIYQQwmfkjIUQQgifkTMWQgghfEbOWAgh\nhPAZOWMhhBDCZ+SMhRBCCJ85MtsPGoaxAqgFhoD3gCbTNF93yzAhhHdIz0L4Sy5PxqtN0zzDNM3J\nQAdwm0s2CSG8R3oWwkeydsamaR6K+/NDwH/mbo4Qwg+kZyH8JetpagDDMFYCjcD/ANNcsUgI4QvS\nsxD+kdIZG4bRBRxr89YtpmluMU1zObDcMIybgDuBL+fBRiGEC0jPQgSXsqGhoZx3YhjGScA20zT/\nIXeThBB+Ij0L4T1ZrxkbhjEp7s8vAs/kbo4Qwg+kZyH8JZc14+8bhmEAg8DvgK+6Y5IQwgekZyF8\nxJVpaiGEEEJkjypwCSGEED4jZyyEEEL4jJyxEEII4TM5Ff3IBsMwfgjMAfqxAkW+bJrm+x7bcDFw\nO3Aa8DnTNJ/2+PgXAD8GyoH1pmmu8vL4ERs2ALOBP5qm+Smvjx9nx4nAA8BHseoi/9Q0zZ/4YEcI\n6AbGAJXAz03TvNlrOyK2lAP7gTdM05zrhw1OCIKWI3ZIzwHQs7Sc1B5HevbjyXgH8EnTNM8AXgT8\nOEm9QD2w2+sDR36YfwEuAE4HLjMM4xNe2wHcH7HBbwaA60zT/CRW1adr/DgfpmmGgbMjtZk/DZxt\nGMYXvLYjwhLgBawbWpAJgpZBeoZg6FlatseRnj13xqZpdpmm+bfIn3uBE3yw4bemab7o9XEjTAVe\nNk3zFdM0B4A2rLxOTzFN8yngT14f18aOt03TfDby/z8DvwGO88mW/4n8txLrKee/vLbBMIwTgIuA\n9UCZ18fPhCBoOWKH9BwAPUvLo8lEz55PUyewAHjIZxu85nggvjXdG0C1T7YECsMwTgY+g3Vj9+P4\nRwBPA38PrDNN8wUfzLgTuBEY58Oxc6EUtQzSsy3S8jCO9ZwXZ5yuBm5km+VAv2maD/plg08EferR\nFwzD+BDwM2BJZFTtOZGnvMmGYRwNbDcM4x9N0/ylV8c3DGMO1prfM4Zh/KNXx01FELTs1A6fkJ4T\nkJYtMtVzXpyxaZrnpXrfMIwmrEf3mfk4vhMbfORN4MS4v0/EGk2XLIZhVACPAhtN0+zw2x7TNN83\nDKMTmAL80sNDfx6oNQzjIiAEjDMM4wHTNK/00IYRBEHLTuzwEek5Dml5BBnp2Y9o6guwHttnRBba\n/cbrdbn9wKTINM4fgEuByzy2ITAYhlEG3Ae8YJrmj3204yPAX03TPGAYxgeA84Bve2mDaZq3ALdE\n7JkB3OCnI05HALUM0rNvSMsjyVTPfkRT34XVvLzLMIxnDMP4314bYBhGvWEYr2NF/HUahvG4V8c2\nTfOvwNeB7VgRdptM0/yNV8ePYhjGQ8CvgVMNw3jdMAy/2uWdCVyBFfH4TOSfH1GhE4FdhmE8i7XO\ntcU0zSd9sCOeoE+B+q5lkJ4hMHqWllOTUs+qTS2EEEL4jCpwCSGEED4jZyyEEEL4jJyxEEII4TNy\nxkIIIYTPyBkLIYQQPiNnLIQQQviMnLEQQgjhM3LGQgghhM/8Xx7rw2IELJpjAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.subplot(121)\n", "plt.scatter(points[:, 0], points[:, 1])\n", "centroids = initialize_centroids(points, 3)\n", "plt.scatter(centroids[:, 0], centroids[:, 1], c='r', s=100)\n", "\n", "plt.subplot(122)\n", "plt.scatter(points[:, 0], points[:, 1])\n", "closest = closest_centroid(points, centroids)\n", "centroids = move_centroids(points, closest, centroids)\n", "plt.scatter(centroids[:, 0], centroids[:, 1], c='r', s=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can animate this type of plot using a module called JSAnimation:" ] }, { "cell_type": "code", "execution_count": 121, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 121, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from JSAnimation import IPython_display\n", "from matplotlib import animation\n", "\n", "# create a simple animation\n", "fig = plt.figure()\n", "ax = plt.axes(xlim=(-4, 4), ylim=(-4, 4))\n", "centroids = initialize_centroids(points, 3)\n", "\n", "def init():\n", " return\n", "\n", "def animate(i):\n", " global centroids\n", " closest = closest_centroid(points, centroids)\n", " centroids = move_centroids(points, closest, centroids)\n", " ax.cla()\n", " ax.scatter(points[:, 0], points[:, 1], c=closest)\n", " ax.scatter(centroids[:, 0], centroids[:, 1], c='r', s=100)\n", " return \n", "\n", "animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=10, interval=200, blit=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The previous animation shows how the points in the dataset are changing their assigned centroid and how the centroid itself moves." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also try running the same type of animation using more clusters than intended in the original data:" ] }, { "cell_type": "code", "execution_count": 122, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " Once \n", " Loop \n", " Reflect \n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 122, "metadata": {}, "output_type": "execute_result" } ], "source": [ "fig = plt.figure()\n", "ax = plt.axes(xlim=(-4, 4), ylim=(-4, 4))\n", "centroids = initialize_centroids(points, 7)\n", "\n", "def init():\n", " return\n", "\n", "def animate(i):\n", " global centroids\n", " closest = closest_centroid(points, centroids)\n", " centroids = move_centroids(points, closest, centroids)\n", " ax.cla()\n", " ax.scatter(points[:, 0], points[:, 1], c=closest)\n", " ax.scatter(centroids[:, 0], centroids[:, 1], c='r', s=100)\n", " return \n", "\n", "animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=30, interval=200, blit=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So that's it for today. I hope you have enjoyed this little notebook that demosntrated how to implement the k-means algorithm and animate it in the browser." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This post was entirely written using the IPython notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at [20150717_Kmeans.ipynb](http://nbviewer.ipython.org/urls/raw.github.com/flothesof/posts/master/20150717_Kmeans.ipynb)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.8" } }, "nbformat": 4, "nbformat_minor": 0 }