k-means_from_scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K-Means Clustering from Sratch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![K-Means Clustering Proccess](https://camo.githubusercontent.com/d191bbb2d798b8add5ca5d1f5b314c567eea5e76/687474703a2f2f6b6f6e756b6f69692e636f6d2f626c6f672f77702d636f6e74656e742f75706c6f6164732f323031372f30312f52756e79616e4b6d65616e732e676966)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-Means with set k=2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will go through the inner workings of k-means clustering and replicate it using only numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "k = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Initiate Cetnroids\n",
    "\n",
    "Randomly place cendtroids, on top of some existing datapoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAG5CAYAAACJG521AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XucXWV97/HPV8AQNUfBRCvhqtIiWkU79RJbLycKWBTSU08Fb9hjS7VeejWl2qMttudF8fRUrVZFS72DihpjQQlFRVuKMlgEgaiIVkJUgniJGoKB3/ljrSmbYS57wt4ze9Z83q/Xfs3ez3rW2s+ancl3r7We9TypKiRJUnfdbaEbIEmShsuwlySp4wx7SZI6zrCXJKnjDHtJkjrOsJckqeMMe0nqQ5K/SPLehW6HtDsMe+kuSPLNJDuSbE/ygyQXJXlRkr7+tpI8KcmWAbXlnUn+apY6leQnSX6c5HtJLkjyrEl1PpPk5naffpTk0iQnJ1k2h7a8Msk32vfZkuQDfazzmSTfn8v7SOqPYS/ddc+oqhXAQcCpwJ8C/7iwTZrRI6rqXsAvAO8E3pTkNZPqvLTdpwcAfwwcD5ybJLNtPMmJwPOAp7TvMwZcMMs6BwO/ChRw7Fx2RtLsDHtpQKrqh1W1EXgWcGKShwEkWZbk/yb5VpLvJnlrkuVJ7gl8AtivPQL+cZL9ktytPZL+env0/cEk+068T5Jfac8g/CDJdUlekOQk4DnA+nY7H++jvTdW1XuAFwN/luS+U9T5SVV9hiaAHwcc07NPr0+ytX28vueI/JeB86rq6+02vlNVp8/SnOcDF9N8+Tixd0F7xuLNSc5pzzZ8PsmDepavSXJJkh+2P9f0LPtMkr9qf18/TvLxJPdN8r72rMUl7ReNifpvaH+nE2c0fnWqxrZtedmkssuTrJtlP6UFYdhLA1ZVXwC20BypAvwN8PPAEcCDgdXAq6vqJ8DTgK1Vda/2sRV4ObAOeCKwH/B94M0ASQ6k+YLw98CqdpuXtWH6PuC0djvPmEOTPwbsCTx6hn36FjDes0+vAh7bvv8j2nX/vF12MfD8JK9IMpZkjz7a8Py2/e8Djkpy/0nLTwD+EtgHuAb4a4D2S9A5wBuB+wL/Dzhn0heX42nONKwGHgT8O/BPwL7A1UDvWY1L2n3aF3g/8KEke0/R3ncBz514keQR7fbP7WNfpXln2EvDsRXYtz3t/TvAH1bVTVW1Hfg/NAE0nd8FXlVVW6pqJ/AXwDOT7Elz9P4vVXVmVf2sqr5XVZfdlYZW1c+AG2kCbtZ9ap8/Bzilqm6oqm00Qfy8dnvvBV4GHAVcCNyQ5OTpNprkV2gugXywqi4Fvg48e1K1j1TVF6pqF80XgiPa8mOAr1XVe6pqV1WdCWwGer/s/FNVfb2qfkjzRenrVfUv7bY+BDyy53fx3vZ3uquq/hZYRnO5Y7KPAYcmObR9/TzgA1V1y3T7KS0kw14ajtXATTRH3/cALm1Pu/8A+GRbPp2DgI/21L8auBW4P3AATRgOTJK92vbcNEvV1T119gP+s2fZf7ZlAFTV+6rqKcB9gBcBpyQ5aprtnghsqqob29fvZ9KpfOA7Pc9/CtxrmnZMtGV1z+vv9jzfMcXriW2R5I+TXN1eEvgBcG9g5eQGt1/CPgg8t+2MeQLwnmn2T1pwhr00YEl+mSZs/pXmiHkH8NCquk/7uHfbcQ2aDmmTXQc8raf+fapq76q6vl32oCnWmW5b/TgO2AV8YYZ9OgD4JeBzbdFWmi8lEw5sy+7YoObsw4eAy4GHTbHd5cBvAk9M8p0k3wH+EHhEe2p8NpPbMdGW6/tYd3JbfpWmc+VvAvtU1X2AHwLTdUp8F80ZjrXAT6vq3+f6ntJ8MeylAUny35I8HTgLeG9VXVFVtwFvB/4uyf3aeqt7jnK/C9w3yb17NvVW4K+THNTWX5XkuHbZ+4CnJPnNJHu2nc2O6NnWA+fQ3n2TPIemP8DfVNX3pqhzjyRPpDlt/QVuvyZ9JvDnbdtWAq8G3tuu84IkxyRZ0XY2fBrwUODzUzRjHc1Zi8NpTs0fATyE5kvF8/vYjXOBn0/y7Pb38ax2W//c56+h1wqaLz3bgD2TvBr4b9NVbsP9NuBv8aheI86wl+66jyfZTnPU/SqaTmK/1bP8T2k6lV2c5EfAv9BeB66qzTTBeW172n4/4A3ARmBTu92Lgce09b8F/BrN7XA3AZfRdJCD5na/w9vtbJihvV9K8uO2Tb9N05/g1ZPqvKl97+8Crwc+DBzdfnkB+CuaDnuXA1cAX2zLAH4EvBL4FvAD4DTgxVX1r1O05USaa+rfanvtf6eqvgO8CXhO209hWu0XlKe3v4/vAeuBp/dcEpiL82iu6X+V5lLAzTSf6UzeDfwi7RcdaVSlanfP/EnS0pbk+cBJVfUrC90WaSYe2UvSbkhyD+D3gNnGEJAWnGEvSXPU9rnYRnOZ4/0L3BxpVp7GlySp4zyylySp42bs6brYrFy5sg4++OCFboYkSfPi0ksvvbGqZhqkC+hY2B988MGMj48vdDMkSZoXSSaPIDklT+NLktRxhr0kSR1n2EuS1HGGvSRJHWfYS5LUcYa9JEkdZ9hLktRxhr0kSR1n2EuS1HGGvSRJHWfYS5LUcYa9JEkdZ9hLktRxQwv7JAck+XSSq5NcmeT3p6iTJG9Mck2Sy5M8qmfZiUm+1j5OHFY7JUnqumFOcbsL+OOq+mKSFcClSc6vqqt66jwNOLR9PAZ4C/CYJPsCrwHGgGrX3VhV3x9ieyVJI2LDZjjtIti6HfZbAevXwLrDFrpVi9fQjuyr6ttV9cX2+XbgamD1pGrHAe+uxsXAfZI8ADgKOL+qbmoD/nzg6GG1VZI0OjZshpMvgOu3N0d7129vXm/YvNAtW7zm5Zp9koOBRwKfn7RoNXBdz+stbdl05VNt+6Qk40nGt23bNqgmS5IWyGkXwY5ddyzbsasp1+4ZetgnuRfwYeAPqupHkxdPsUrNUH7nwqrTq2qsqsZWrVp11xorSVpwW7fPrVyzG2rYJ9mLJujfV1UfmaLKFuCAntf7A1tnKJckddx+K+ZWrtkNszd+gH8Erq6q/zdNtY3A89te+Y8FflhV3wbOA45Msk+SfYAj2zJJUsetXwPLJ3UfX75nU67dM8ze+I8HngdckeSytuyVwIEAVfVW4Fzg14BrgJ8Cv9UuuynJa4FL2vVOqaqbhthWSQNmb2rtrol/JxP/fu69DBL4g/OaMv8tzV2qprwUviiNjY3V+Pj4QjdDWvImelP3drJaviecutb/pDU3/luaWZJLq2pstnqOoCdp4OxNrUHx39JgGPaSBs7e1BoU/y0NhmEvaeDsTa1B8d/SYBj2kgbO3tQalLn8W9qwGdacAQe/ofnpiHu3G2ZvfElL1OTe1PbG1+7q99/S5I58E0Ps9m5jKbM3viRp0VtzRhPwk61eARf9r/lvz3yxN74kacmwI9/MDHtJ0qJnR76ZGfaSpEWv3458S7UTnx30JEmLXj8d+ZZyJz7DXpLUCesOmzm0ZxqNr+th72l8SdKSsJQ78Rn2kqQlYSl34jPsJUlLwlIe2dFr9pKkJWEpj+xo2EuSlozZOvF1lafxJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjttzoRsgSbqzDZvhtItg63bYbwWsXwPrDlvoVmmxGlrYJzkDeDpwQ1U9bIrlrwCe09OOhwCrquqmJN8EtgO3AruqamxY7ZSkUbNhM5x8AezY1by+fnvzGgx87Z5hnsZ/J3D0dAur6nVVdURVHQH8GXBhVd3UU+XJ7XKDXtKSctpFtwf9hB27mnJpdwwt7Kvqs8BNs1ZsnACcOay2SNJisnX73Mql2Sx4B70k96A5A/DhnuICNiW5NMlJs6x/UpLxJOPbtm0bZlMlaV7st2Ju5dJsFjzsgWcA/zbpFP7jq+pRwNOAlyR5wnQrV9XpVTVWVWOrVq0adlslaejWr4Hlk3pULd+zKZd2xyiE/fFMOoVfVVvbnzcAHwUevQDtkqQFse4wOHUtrF4Bofl56lo752n3Leitd0nuDTwReG5P2T2Bu1XV9vb5kcApC9RESVoQ6w4z3DU4w7z17kzgScDKJFuA1wB7AVTVW9tqvw5sqqqf9Kx6f+CjSSba9/6q+uSw2ilJUtcNLeyr6oQ+6ryT5ha93rJrgUcMp1WSJPWnSwMbOYKeJEmTdG1go1HooCdJ0kjp2sBGhr0kSZN0bWAjw16SpEm6NrCRYS9J0iRdG9jIDnqSJE0y0QnP3viSJHVYlwY28jS+JEkdZ9hLktRxhr0kSR1n2EuS1HGGvSRJHWfYS5LUcYa9JEkdZ9hLktRxhr0kSR1n2EuS1HGGvSRJHWfYS5LUcYa9JEkdZ9hLktRxTnErSR21YXN35mPXXWPYS1IHbdgMJ18AO3Y1r6/f3rwGA38p8jS+JHXQaRfdHvQTduxqyrX0GPaS1EFbt8+tXN1m2EvqnA2bYc0ZcPAbmp8bNi90i+bffivmVq5uM+wldcrEterrt0Nx+7XqpRb469fA8km9spbv2ZRr6THsJXWK16ob6w6DU9fC6hUQmp+nrrVz3lJlb3xJneK16tutO8xwny+jfpujYS+pU/Zb0Zy6n6pcGqSJgL9+e3P2pNry3tscYTS+BBj2kjpl/Zo73l8OXqvW4E0ex6AmLd+xC/7iQrh512iMdeA1e0md4rVqzYep+oZM9v2bR6f/iEf2kjrHa9UatrvSB2Qh+o94ZC9J0hzN1gdk+Z5wn2W7t+4wGPaSJM3RVOMYpP05cenoL580OmMdDO00fpIzgKcDN1TVw6ZY/iTgY8A32qKPVNUp7bKjgTcAewDvqKpTh9VOSZLmauIyUT897UehN36qJvchHNCGkycAPwbePUPY/0lVPX1S+R7AV4GnAluAS4ATquqq2d5zbGysxsfHB9B6SZJGX5JLq2pstnpDO41fVZ8FbtqNVR8NXFNV11bVLcBZwHEDbZwkSUvIQl+zf1ySLyX5RJKHtmWrget66mxpy6aU5KQk40nGt23bNsy2SlLfnIxHo2Qhw/6LwEFV9Qjg74ENbXmmqDvttYaqOr2qxqpqbNWqVUNopiTNjZPxaNQsWNhX1Y+q6sft83OBvZKspDmSP6Cn6v7A1gVooiTtFifj0ahZsLBP8nNJ0j5/dNuW79F0yDs0ySFJ7g4cD2xcqHZK0lw5GY9GzTBvvTsTeBKwMskW4DXAXgBV9VbgmcCLk+wCdgDHV3NrwK4kLwXOo7n17oyqunJY7ZSkQXMyHo2aoYV9VZ0wy/I3AW+aZtm5wLnDaJckDZuT8WjUODa+JA3YXAZckeaDYS9JQ+BkPBolC32fvSRJGjKP7CVJQ7Nhs5czRoFhL0kaionBhSY6Kk4MLgQG/nzzNL4kaSgcXGh0GPaSpKFwcKHRYdhLkoZiukGEHFxo/hn2ktQRozbT3vo1zWBCvRxcaGHYQU+SOmDQneEG0YvewYVGh2EvSR0wU2e4uYbrIL84OLjQaPA0viR1wCA7w9mLvnsMe0nqgEF2hrMXffcY9pLUAYPsDGcv+u4x7CWpA9YdBqeuhdUrIDQ/T127e9fL7UXfPXbQk6SOGFRnOHvRd49hL0m6E3vRd4un8SVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjnNQHUnqiEHMQa9uMuwlqQMGOQe9usfT+JLUAc5Br5kY9pLUAc5Br5kY9pLUAc5Br5kY9pLUAc5Br5nYQU+SOsA56DUTw16SOsI56DUdT+NLktRxhr0kSR03tLBPckaSG5J8eZrlz0lyefu4KMkjepZ9M8kVSS5LMj6sNkqStBQM88j+ncDRMyz/BvDEqno48Frg9EnLn1xVR1TV2JDaJ0nSkjC0DnpV9dkkB8+wvHdcp4uB/YfVFkmSlrJRuWb/QuATPa8L2JTk0iQnzbRikpOSjCcZ37Zt21AbKWkwNmyGNWfAwW9ofm7YvNAtkrptwW+9S/JkmrD/lZ7ix1fV1iT3A85PsrmqPjvV+lV1Ou0lgLGxsRp6gyXdJU7YIs2/BT2yT/Jw4B3AcVX1vYnyqtra/rwB+Cjw6IVpoaRBc8IWaf4tWNgnORD4CPC8qvpqT/k9k6yYeA4cCUzZo1/S4uOELdL8mzXsk/x+P2VT1DkT+HfgF5JsSfLCJC9K8qK2yquB+wL/MOkWu/sD/5rkS8AXgHOq6pN97o+kEeeELdL8S9XMl7mTfLGqHjWp7D+q6pFDbdluGBsbq/Fxb8uXRtnka/bQTNhy6lqv2UtzleTSfm5Rn7aDXpITgGcDhyTZ2LNoBfC9qdeSpJk5YYs0/2bqjX8R8G1gJfC3PeXbgcuH2ShJ3eaELdL8mjbsq+o/gf8EHpfkIODQqvqXJMuB5TShL0mSRlw/HfR+BzgbeFtbtD+wYZiNkiRJg9PPrXcvAR4P/Aigqr4G3G+YjZIkSYPTT9jvrKpbJl4k2ZNmOFtJkrQI9BP2FyZ5JbA8yVOBDwEfH26zJEnSoPQT9icD24ArgN8Fzq2qVw21VZIkaWD6mQjnkVX1duDtEwVJnlFVHt1LkrQI9HNk//Ykvzjxoh1s58+H1yRJkjRI/RzZPxM4O8lzaKahfT7N5DSSJGkRmDXsq+raJMfT3Ft/HXBkVe0YesskSdJAzDQ2/hXc8Ra7fYE9gM8noaoePuzGSZKku26mI/unz1srJEnS0Mw2Nj5JHgtcWVXb29crgMNpxs2XJC1iGzY7A+FS0E9v/LcAP+55/ZO2TJK0iG3YDCdfANdvb67ZXr+9eb1h80K3TIPWT9inqv7r2n1V3UZ/vfglSSPstItgx647lu3Y1ZSrW/oJ+2uTvDzJXu3j94Frh90wSdJwbZ1movLpyrV49RP2LwLWANcDW4DHACcNs1GSpOHbb8XcyrV4zRr2VXVDVR1fVferqvtX1bOr6ob5aJwkaXjWr4Hlky7KLt+zKVe3zHrtPcnewAuBhwJ7T5RX1f8aYrskSUM20eve3vjd109Hu/cAm4GjgFOA5wBXD7NRkqT5se4ww30p6Oea/YOr6n8DP6mqdwHHAL84yzqSJGlE9BP2P2t//iDJw4B7AwcPrUWSJGmg+jmNf3qSfYD/DWwE7tU+lyRJi0A/s969o316IfDA4TZHknafQ79KU5v1NH6S+yb5+yRfTHJpktcnue98NE6S+uXQr9L0+rlmfxZwA/AbwDOBG4EPDLNRkjRXDv0qTa+fa/b7VtVre17/VZJ1w2qQJO0Oh36VptfPkf2nkxyf5G7t4zeBc4bdMEmaC4d+labXT9j/LvB+YGf7OAv4oyTbk/xomI2TpH459Ks0vX564/u9WNLIc+hXaXrOSy+pMxz6VZpaP6fxJUnSIjbUsE9yRpIbknx5muVJ8sYk1yS5PMmjepadmORr7ePEYbZTkqQumzHs2973UwZ1n94JHD3D8qcBh7aPk4C3tO+7L/Aa4DHAo4HXtEP2SpKkOZox7KvqNuBLSQ7cnY1X1WeBm2aochzw7mpcDNwnyQNoptM9v6puqqrvA+cz85cGSZI0jX466D0AuDLJF4CfTBRW1bEDeP/VwHU9r7e0ZdOV30mSk2jOCnDggbv1nUSSpE7rJ+z/cojvnynKaobyOxdWnQ6cDjA2NjZlHUmSlrJZO+hV1YXAN4G92ueXAF8c0PtvAQ7oeb0/sHWGckmSNEf9zHr3O8DZwNvaotXAhgG9/0bg+W2v/McCP6yqbwPnAUcm2aftmHdkWyZJkuaon9P4L6HpEf95gKr6WpL79bPxJGcCTwJWJtlC08N+r3Y7bwXOBX4NuAb4KfBb7bKbkryW5iwCwClVNVNHP0mSNI1+wn5nVd2SNJfRk+zJNNfPJ6uqE2ZZXjRfJqZadgZwRj/vI0mSptfPoDoXJnklsDzJU4EPAR8fbrMkSdKg9BP2JwPbgCtoZsA7F/jzYTZKkiQNTj+z3t2W5F001+wL+Ep7+l2SJC0Cs4Z9kmOAtwJfp7n//ZAkv1tVnxh24yRJ0l3XTwe9vwWeXFXXACR5EHAOYNhLkrQI9HPN/oaJoG9dC9wwpPZIkqQB6+fI/sok5wIfpLlm/z+BS5L8D4Cq+sgQ2ydJku6ifsJ+b+C7wBPb19uAfYFn0IS/YS9J0gjrpzf+b81HQyRJ0nD0c81ekiQtYoa9JEkdZ9hLktRx/Uxxe/8k/5jkE+3rw5O8cPhNkyRJg9DPkf07aeaS3699/VXgD4bVIEmSNFj9hP3KqvogcBtAVe0Cbh1qqyRJ0sD0E/Y/SXJf2jnskzwW+OFQWyVJkgamn0F1/hjYCDwoyb8Bq2hG0ZMkSYtAP4PqXJrkicAv0Mx695Wq+tnQWyZJkgain974Xwd+u6qurKovV9XPkvzzPLRNkiQNQD/X7H8GPDnJPyW5e1u2eohtkiRJA9RP2P+0qp4FXA18LslBtJ31JEnS6Oung14Aquq0JJfS3HO/71BbJUmSBqafsH/1xJOquiDJUcCJw2uSJEkapGnDPslhVbUZuD7JoyYttoOeJEmLxExH9n8EnAT87RTLCvjvQ2mRJEkaqGnDvqpOan8+ef6aI0mSBm3a3vhJfjnJz/W8fn6SjyV5YxI76EmStEjMdOvd24BbAJI8ATgVeDfNuPinD79pkiRpEGa6Zr9HVd3UPn8WcHpVfRj4cJLLht80SZI0CDMd2e+RZOLLwFrgUz3L+rllT5IkjYCZQvtM4MIkNwI7gM8BJHkwTnErSdKiMVNv/L9OcgHwAGBTVU0MkXs34GXz0ThJknTXzXg6vqounqLsq8NrjiRJGrR+JsKRJEmLmGEvSVLHDTXskxyd5CtJrkly8hTL/y7JZe3jq0l+0LPs1p5lG4fZTkmSumxot9Al2QN4M/BUYAtwSZKNVXXVRJ2q+sOe+i8DHtmziR1VdcSw2idJ0lIxzCP7RwPXVNW1VXULcBZw3Az1T6C53U+SJA3QMMN+NXBdz+stbdmdJDkIOIQ7Dtyzd5LxJBcnWTfdmyQ5qa03vm3btkG0W5KkThlm2GeKspqiDOB44OyqurWn7MCqGgOeDbw+yYOmWrGqTq+qsaoaW7Vq1V1rsSRJHTTMsN8CHNDzen9g6zR1j2fSKfyq2tr+vBb4DHe8ni9Jkvo0zLC/BDg0ySFJ7k4T6HfqVZ/kF4B9gH/vKdsnybL2+Urg8cBVk9eVJEmzG1pv/KraleSlwHnAHsAZVXVlklOA8aqaCP4TgLN6huMFeAjwtiS30XwhObW3F78kSepf7pixi9vY2FiNj48vdDMkSZoXSS5t+7fNyBH0JEnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6jjDXpKkjjPsJUnqOMNekqSOM+wlSeo4w16SpI4z7CVJ6rihhn2So5N8Jck1SU6eYvkLkmxLcln7+O2eZScm+Vr7OHGY7ZQkqcv2HNaGk+wBvBl4KrAFuCTJxqq6alLVD1TVSyetuy/wGmAMKODSdt3vD6u9i9mGzXDaRbB1O+y3AtavgXWH7X49SVK3DPPI/tHANVV1bVXdApwFHNfnukcB51fVTW3Anw8cPaR23mUbNsOaM+DgNzQ/N2ye3/c++QK4fnvzrej67c3ryW3ot54kqXuGGfarget6Xm9pyyb7jSSXJzk7yQFzXJckJyUZTzK+bdu2QbR7ThY6RE+7CHbsumPZjl1N+e7UkyR1zzDDPlOU1aTXHwcOrqqHA/8CvGsO6zaFVadX1VhVja1atWq3G7u7FjpEt27vr7zfepKk7hlm2G8BDuh5vT+wtbdCVX2vqna2L98O/FK/646KhQ7R/Vb0V95vPUlS9wwz7C8BDk1ySJK7A8cDG3srJHlAz8tjgavb5+cBRybZJ8k+wJFt2chZ6BBdvwaWT+pmuXzPpnx36kmSumdoYV9Vu4CX0oT01cAHq+rKJKckObat9vIkVyb5EvBy4AXtujcBr6X5wnAJcEpbNnIWOkTXHQanroXVK5prH6tXNK8n97Lvt54kqXtSNeWl8EVpbGysxsfH5/19vaVNkrQQklxaVWOz1RvaffZLybrDDHdJ0uhyuFxJkjrOI3vNyEsUkrT4Gfaa1sSAQRPjCEwMGAQGviQtJp7G17QWesAgSdJgGPaa1kIPGCRJGgzDXtNa6AGDJEmDYdhrWgs9YJAkaTDsoLdILESv+Int2xtfkhY3w34RWMhe8Q4YJEmLn6fxZ7FhM6w5Aw5+Q/Nzvuap72WveEnSXWHYz2DiiPr67VDcfkQ934H/9Usv4Ia3HMXNX73gDuX2ipck9cOwn8EoHFFv2rSJbacfw81f2cQNpx/Djs2b/mvZYu8VPwpnTSRpKTDsZ7DQ95lv2rSJY489ltq1k70fuDfs2sm2dxzLjs2bFn2v+FE5ayJJS4FhP4OFvM98Iuh37tzJymeu5CFnPYSVz1wJu3Zy4zuO5fi9Ny3qjnOjcNZEkpYKw34GC3Wf+eSgP+BPDyB7hAP+9ABWPnMltWsnp77kWDZt2jT7xkbUQp81kaSlxLCfwbrD4NS1sHoFhObnqWuHeyvalEF/twCQu90e+Dt37uTYY0c78Ge6Ju/ofJI0f7zPfhbzeZ/5BRdcMG3QT5gIfIAbz76RY489lnPOOYe1a9fOTyP7NNvYAOvX3HE5ODqfJA2LR/Yj5LTTTmPnzqYz3gHr7xz0E3K3cMD6A9j7gXuzc+dOXve6181zS2c32zX5hThrIklLlUf2I2T9+vVceOGF3HztzVx32nVTHtkD1G3Fdaddx83X3syyZct4xStesQCtnVk/1+QdnU+S5odH9iNk7dq1bNy4kWU98ahSAAARV0lEQVTLlnHj2Tdy3d9cR91Wd6hTtxXX/c113Hj2jSxbtoyNGzeO3Cl88Jq8JI0Sw37EHHnkkdMG/lRBf+SRRy5wi6fmjHmSNDoM+xEy0Xv9pKuP5IG/t5G97t4T+LcunqAHr8lL0ihJVc1ea5EYGxur8fHxhW7Gbpncex3gtq9t4rtvP5af3dJ02pu4Rj+IoL/gggs47bTTWL9+/UheBpAkzS7JpVU1Nls9j+xHxFS91+926JE8+CXNKf1BBv2mTZs45phj7vBTktRdhv2ImK73+k8POpJzzjmHo446inPOOWcgQT9xL//hsCgG55Ek3TWG/YiYqff62rVr+eQnP3mXT7f3Bv2LgcuBF2PgS1LXGfYjYti91ycH/ZuAPdqfBr4kdZthPyKG2Xt9qqCf+ODvhoEvSV1nb/yOu+CCCzjmmGOmDPpetwEvBd4CLFu2bCTH25ck3ZG98QXcPt7+4cDfM/0Hfrd2+USnvVEcb1+StHsM+45bv349y5Yt4yrgZTRH8FO5rV1+FYzsePuSpN1j2C8iM80PP53e8fbfQnOqfnLgTz6FP6rj7UuSdo9hv0hMjLB3/XYobp8fvp/A7x1vf3LgTxX0ozwMryRp7gz7RWK2+eFnM1Xg34pBL0lLwVDDPsnRSb6S5JokJ0+x/I+SXJXk8iQXJDmoZ9mtSS5rHxuH2c7FoJ/54WczOfAfjkEvSUvB0MI+yR7Am4Gn0XTyPiHJ4ZOq/QcwVlUPB84GTutZtqOqjmgfxw6rnYvFoOaH7w38ic54Br0kddswj+wfDVxTVddW1S3AWcBxvRWq6tNV9dP25cXA/kNsz6I2yBH2jjxysOPtS5JG256zV9ltq4Hrel5vAR4zQ/0XAp/oeb13knFgF3BqVW2YaqUkJwEnARx44IF3qcGjbGIkvdMuak7d77eiCfrdHWFv7dq19riXpCVimGGfKcqmHK4vyXOBMeCJPcUHVtXWJA8EPpXkiqr6+p02WHU6cDo0I+jd9WaPrnWHDWb4XEnS0jLM0/hbgAN6Xu8PbJ1cKclTgFcBx1bVzonyqtra/rwW+AzwyCG2VZKkzhpm2F8CHJrkkCR3B44H7tCrPskjgbfRBP0NPeX7JFnWPl8JPJ5mcDdJkjRHQzuNX1W7krwUOI9mNtUzqurKJKcA41W1EXgdcC/gQ0kAvtX2vH8I8LYkt9F8ITm1qgz7Ebdh8+D6FEiSBsdZ7wZsqQbexAh/vQP/LN9zcNP0SpLuzFnvFsBdGdJ2sburI/xJkobHsB+gpRx4gxjhT5I0HIb9AC3lwBvUCH+SpMEz7AdoKQfeIEf4kyQNlmE/QEs58NYd1nTGW72iGU1p9Qo750nSqBjmCHpLzqCHtF1sHOFPkkaTYT9gkwNvw2ZYc8bSDH9J0mgw7Ido8r3nE7figYEvSZo/hv2ATDWYzky34hn2kqT5YtgPwHRH8JODfsJSuBVPkjQ67I0/ANMdwe8x1SS/LI1b8SRJo8OwH4DpjtRvraV7K54kaXQY9gMw3ZH6xL3m83Xv+UTP/4Pf0PxcCmPyS5Jm5zX7AVi/ZuoZ3yZus5uPznj2/JckTccj+wEYhdHjlvIkPJKkmXlkPyALPXrcUp6ER5I0M4/sO2IpT8IjSZqZYd8RS3kSHknSzDyNP42pRsQb5Y5uozQJz2L73UlS1xn2U1isPdsXut8ALN7fnSR1mafxp2DP9t3n706SRo9hPwV7tu8+f3eSNHoM+ynYs333+buTpNFj2E/Bnu27z9+dJI0eO+hNYZR6ti82/u4kafSkqha6DQMzNjZW4+PjC90MSZLmRZJLq2pstnqexpckqeMMe0mSOs5r9pqWI+FJUjcY9pqSI+FJUnd4Gl9TciQ8SeoOw15TciQ8SeoOw15TciQ8SeoOw15TciQ8SeoOO+hpSo6EJ0ndMdSwT3I08AZgD+AdVXXqpOXLgHcDvwR8D3hWVX2zXfZnwAuBW4GXV9V5w2yr7mzdYYa7JHXB0E7jJ9kDeDPwNOBw4IQkh0+q9kLg+1X1YODvgL9p1z0cOB54KHA08A/t9iRJ0hwN85r9o4FrquraqroFOAs4blKd44B3tc/PBtYmSVt+VlXtrKpvANe025MkSXM0zLBfDVzX83pLWzZlnaraBfwQuG+f6wKQ5KQk40nGt23bNqCmS5LUHcMM+0xRNnmKvenq9LNuU1h1elWNVdXYqlWr5thESZK6b5hhvwU4oOf1/sDW6eok2RO4N3BTn+tKkqQ+DDPsLwEOTXJIkrvTdLjbOKnORuDE9vkzgU9VVbXlxydZluQQ4FDgC0NsqyRJnTW0W++qaleSlwLn0dx6d0ZVXZnkFGC8qjYC/wi8J8k1NEf0x7frXpnkg8BVwC7gJVV167DaKklSl6U5kO6GsbGxGh8fX+hmSJI0L5JcWlVjs9VzuFxJkjrOsJckqeMMe0mSOs6wlySp4wx7SZI6zrCXJKnjOnXrXZJtwH8OYFMrgRsHsJ1R4L6MJvdlNLkvo6kr+zKM/TioqmYdK75TYT8oScb7uW9xMXBfRpP7Mprcl9HUlX1ZyP3wNL4kSR1n2EuS1HGG/dROX+gGDJD7Mprcl9HkvoymruzLgu2H1+wlSeo4j+wlSeo4w16SpI5bcmGf5OgkX0lyTZKTp1i+LMkH2uWfT3Jwz7I/a8u/kuSo+Wz3ZH3sxx8luSrJ5UkuSHJQz7Jbk1zWPjbOb8vvrI99eUGSbT1t/u2eZScm+Vr7OHF+W35nfezL3/Xsx1eT/KBn2ah9LmckuSHJl6dZniRvbPf18iSP6lk2ap/LbPvynHYfLk9yUZJH9Cz7ZpIr2s9lwefQ7mNfnpTkhz3/ll7ds2zGf5/zqY/9eEXPPny5/fvYt102ap/JAUk+neTqJFcm+f0p6izs30tVLZkHsAfwdeCBwN2BLwGHT6rze8Bb2+fHAx9onx/e1l8GHNJuZ48R3o8nA/don794Yj/a1z9e6M9ijvvyAuBNU6y7L3Bt+3Of9vk+o7wvk+q/DDhjFD+Xtj1PAB4FfHma5b8GfAII8Fjg86P4ufS5L2sm2gg8bWJf2tffBFYu9Ocxh315EvDPU5TP6d/nQu/HpLrPAD41wp/JA4BHtc9XAF+d4v+xBf17WWpH9o8Grqmqa6vqFuAs4LhJdY4D3tU+PxtYmyRt+VlVtbOqvgFc025vIcy6H1X16ar6afvyYmD/eW5jv/r5TKZzFHB+Vd1UVd8HzgeOHlI7+zHXfTkBOHNeWrYbquqzwE0zVDkOeHc1Lgbuk+QBjN7nMuu+VNVFbVthtP9e+vlcpnNX/tYGbo77Mep/K9+uqi+2z7cDVwOrJ1Vb0L+XpRb2q4Hrel5v4c4fyH/VqapdwA+B+/a57nyZa1teSPONcsLeScaTXJxk3TAaOAf97stvtKe+zk5ywBzXnS99t6e9rHII8Kme4lH6XPox3f6O2ucyV5P/XgrYlOTSJCctUJvm6nFJvpTkE0ke2pYtys8lyT1owu/DPcUj+5mkufT7SODzkxYt6N/LnoPe4IjLFGWT7z2crk4/686XvtuS5LnAGPDEnuIDq2prkgcCn0pyRVV9fQjt7Ec/+/Jx4Myq2pnkRTRnXv57n+vOp7m053jg7Kq6tadslD6XfiyGv5U5SfJkmrD/lZ7ix7efy/2A85Nsbo9KR9UXacZL/3GSXwM2AIeyeD+XZwD/VlW9ZwFG8jNJci+aLyV/UFU/mrx4ilXm7e9lqR3ZbwEO6Hm9P7B1ujpJ9gTuTXOqqZ9150tfbUnyFOBVwLFVtXOivKq2tj+vBT5D8y10ocy6L1X1vZ72vx34pX7XnWdzac/xTDotOWKfSz+m299R+1z6kuThwDuA46rqexPlPZ/LDcBHWbjLd32pqh9V1Y/b5+cCeyVZySL9XJj5b2VkPpMke9EE/fuq6iNTVFnYv5eF7tgwnw+aMxnX0pw+neig8tBJdV7CHTvofbB9/lDu2EHvWhaug14/+/FIms44h04q3wdY1j5fCXyNhe2k08++PKDn+a8DF7fP9wW+0e7TPu3zfUd5X9p6v0DTwSij+rn0tOtgpu8Idgx37HD0hVH8XPrclwNp+uGsmVR+T2BFz/OLgKNHfF9+buLfFk0Ifqv9jPr69zkq+9EunzjYuucofybt7/fdwOtnqLOgfy9L6jR+Ve1K8lLgPJqeqWdU1ZVJTgHGq2oj8I/Ae5JcQ/OP7Ph23SuTfBC4CtgFvKTueAp21PbjdcC9gA81/Qv5VlUdCzwEeFuS22jO7JxaVVctxH5A3/vy8iTH0vzeb6LpnU9V3ZTktcAl7eZOqTue6ptXfe4LNJ2Nzqr2L701Up8LQJIzaXp2r0yyBXgNsBdAVb0VOJemh/E1wE+B32qXjdTnAn3ty6tp+ub8Q/v3squa2cnuD3y0LdsTeH9VfXLed6BHH/vyTODFSXYBO4Dj239rU/77XIBdAPraD2i+3G+qqp/0rDpynwnweOB5wBVJLmvLXknzJXIk/l4cLleSpI5batfsJUlacgx7SZI6zrCXJKnjDHtJkjrOsJckqeMMe2kE5PYZ776c5EPtEKFz3cZfJPmTOdRfl+Twub5Pl7Szp61c6HZIw2bYS6NhR1UdUVUPA24BXjQP77mOZjbHeZFkj/l6L0l3ZNhLo+dzwIMBkmxoJ/u4snfCj3Ze8i+2k51cMHkDSX6nnQRleZIHJflku53PJTksyRrgWOB17RmFByV5eZKr2gmHzppimy9I8rF2W19J8pqeZc9N8oV2W2+bCPYkP05ySpLPA4+btL07vV+Se6aZ5/ySJP+R5Li2fI8k/zfNHOaXJ3lZW762rXdFu96ytvybSf6y/R1dkeSwtvy+STa167yNqccll7pnIYcY9OHDR/OgncueZkSwjwEvbl/v2/5cDnyZZpS3VTSzZB0yqc5fAH8CvBTYyO3D715AO2wy8BjaecGBdwLP7GnD1p517jNFG18AfLttw0R7xmhG//s4sFdb7x+A57fPC/jNafb5Tu8H/B/guRNlNPOC3xN4Mc2443tO7DOwd/t7+Pm27N00E5BAMxzxy9rnvwe8o33+RuDV7fNj2vaNzLzoPnwM67GkhsuVRtjynmE2P0czbDM0QwX/evv8AJrZy1YBn62qb0Az3GbPdp5HM7HGuqr6WTsL1xpuHzYZmvkdpnI58L4kG2hmSpvK+dVOEpPkIzSzw+2imZzokvY9lgM3tPVv5Y5Tk872fkcCx/b0PdibZsjRp9DMWbFrYp+TPAL4RlV9ta37Lpq5LV7fvp6YjORS4H+0z58w8byqzkkyMYe91GmGvTQadlTVEb0FSZ5EE3KPq6qfJvkMTfiF6afA/DJwBM3MWd+guVT3g8nbnsYxNGF4LPC/kzx0Ilx7TH7fiSk631VVfzbFNm+u6eeQuNP7tdv6jar6Sm/FNN8i+pmOutfETIm3csf/6xwjXEuO1+yl0XVv4Ptt0B9GM1MWwL8DT0xyCECSfXvW+Q/gd4GNSfarZk7tbyT5n23dtEfEANuBFW353YADqurTwHqaU+j3mqJNT02yb5LlNB38/o3mMsEz08wtTrv8oJl2bIb3Ow94WRvuJJmY5ncT8KI0005P7PNm4OAkD27rPA+4cKb3BT4LPKfdxtNoZhmTOs+wl0bXJ4E9k1wOvBa4GKCqtgEnAR9J8iXgA70rVdW/0ly7P6e9rew5wAvbulcCx7VVzwJekeQ/aC4PvDfJFTRfGP6uqn4wRZv+FXgPcBnw4aoar2Z2vj8HNrVtPR94wCz7tsc07/dampnPLk/y5fY1NPPMf6st/xLw7Kq6mWbmsA+127kNeCsz+0vgCUm+SHPJ4Fuz1Jc6wVnvJPUlyQuAsap66UK3RdLceGQvSVLHeWQvSVLHeWQvSVLHGfaSJHWcYS9JUscZ9pIkdZxhL0lSx/1/hju/S16uqsgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a716af5f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eucledean distance between points\n",
    "def eucl_distance(a, b):\n",
    "    return np.linalg.norm(a-b)\n",
    "\n",
    "# import data\n",
    "dataset = np.loadtxt('durudataset.txt')\n",
    "num_instances, num_features = dataset.shape # dataset has a shape of (40, 2)\n",
    "\n",
    "# randomly select k centroids' starting ponts (from data)\n",
    "centroids = dataset[np.random.randint(num_instances, size=k)]\n",
    "\n",
    "# how much did cetroid move till last iteration\n",
    "centroids_old = np.zeros(centroids.shape) # define any \"centroid_old\" matrix\n",
    "centroids_diff = eucl_distance(centroids, centroids_old)\n",
    "\n",
    "# which cluster does each datapoint belong to\n",
    "belongs_to = np.zeros((num_instances, 1))\n",
    "\n",
    "# plot data & centroids\n",
    "colors = ['red', 'limegreen']\n",
    "plt.figure(figsize=(8,7))\n",
    "plt.scatter(dataset[:,0], dataset[:,1], label='Data', color='dodgerblue')\n",
    "for i, c in enumerate(centroids):\n",
    "    plt.scatter(centroids[i,0], centroids[i,1], color=colors[i], marker=\"D\", s=80, edgecolor='black', linewidth='2')\n",
    "plt.title('Detect DDoS Anomaly')\n",
    "plt.xlabel('Packets per second')\n",
    "plt.ylabel('Size per packet')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign Clusters\n",
    "\n",
    "List the index of the datapoints that belong to each cluster. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAGfCAYAAAByXooxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+MXOdd7/HPx4nrxS0CJ97r1Il3N703F6W96o8whEKuaIqpY9fUBlEJ525DilqtKFko9Kql4KiR8kOCVuKnS9sFrNJqSQqBInMTNw5pUa8oKR5HIWmcm8Zx7fXKTbKxS9Ni7MTx9/5xziTj8czO2fXszHlm3i9pNXOe88zsc3I8+exzzjPP44gQAAAov2W9bgAAACiG0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAk4sJeN6CZ1atXx9jYWK+bAQBAV+zbt++5iBhuV6+UoT02NqZqtdrrZgAA0BW2Dxepx+VxAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiWgb2rbX2f6K7cdtP2b7g03q2PYf2z5g+xHbV9Xtu9H2k/nPjZ0+AAAABkWRnvZpSf87Iq6U9FZJN9l+fUOdTZKuyH8mJH1KkmxfJOkWST8u6WpJt9he1aG2AwDKbHpaGhuTli3LHqene92i5LUN7Yj4dkQ8lD//nqTHJV3aUG2rpM9F5kFJP2z7tZKuk3R/RByPiO9Iul/Sxo4eAQCgfKanpYkJ6fBhKSJ7nJgguM/Tgu5p2x6T9BZJX2/YdamkI3Xbs3lZq3IAQD/bvl06ceLsshMnsnIsWuHQtv0aSX8r6Tci4vnG3U1eEvOUN3v/CdtV29W5ubmizQIAlNHMzMLKUUih0La9XFlgT0fE3zWpMitpXd32ZZKOzlN+joiYiohKRFSGh4eLNAsAUFYjIwsrRyFFRo9b0l9Iejwifr9FtV2SfikfRf5WSd+NiG9Luk/SBtur8gFoG/IyAEA/u+MOaeXKs8tWrszKsWhFetrXSLpB0k/bfjj/eaftX7H9K3mdeyUdlHRA0p9J+lVJiojjkm6TtDf/uTUvA1Byu4/t1uZHN6vyUEWbH92s3cd297pJSMn4uDQ1JY2OSrZ08cXSD/yAdMMNjCQ/D45oeou5pyqVSlSr1V43AxhYu4/t1u0zt+tknHy5bMhDunnkZm26eFMPW4Yk1UaS1w9MW7kyC/Xx8d61q0Rs74uISrt6zIgG4Bw7ju44K7Al6WSc1I6jO3rUIiSNkeQdQ2gDOMczLz6zoHJgXowk7xhCG8A51ixfs6ByYF6MJO8YQhvAOSbXTmrIQ2eVDXlIk2sne9QiJK3oSHKmPW3rwl43AED51Aab7Ti6Q8+8+IzWLF+jybWTDELD4tQGm23fnl0SHxnJArt+EFrjYLXatKf1rwejxwEAJTA2lgV1o9FR6dChbrem6xg9DgBIB4PVCiG0AQC9x2C1QghtAEDvMe1pIYQ2AKD3Gqc9HR1tPmPagI8wZ/Q4AKAcxsfnHynOCHN62gCARDAdKqENAEgEI8wJbQBAIhhhTmgDABLBCHNCGwCQiKIjzPsYo8cBAOloN8K8z9HTBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYALIXpaWlsTFq2LHucnu51i9AHLmxXwfZOST8r6dmI+B9N9n9Y0njd+10paTgijts+JOl7kl6SdDoiKp1qOACU1vS0NDEhnTiRbR8+nG1L0vh469cBbTgi5q9g/5Sk70v6XLPQbqj7Lkm/GRE/nW8fklSJiOcW0qhKpRLVanUhLwGA8hgby4K60eiodOhQt1uDBNjeV6Rj2/byeER8VdLxgr/3ekl3FqwLAP1pZmZh5UBBHbunbXulpI2S/rauOCTtsb3P9kSb10/Yrtquzs3NdapZANB9IyMLKwcK6uRAtHdJ+ueIqO+VXxMRV0naJOmm/FJ7UxExFRGViKgMDw93sFkA0GV33CGtXHl22cqVWTlwHjoZ2tvUcGk8Io7mj89K+qKkqzv4+wCgnMbHpamp7B62nT1OTTEIDeet7ejxImz/kKS3SXpPXdmrJS2LiO/lzzdIurUTvw8ASm98nJBGxxX5ytedkq6VtNr2rKRbJC2XpIj4dF7t5yXtiYj/qHvpGklftF37PX8VEV/qXNMBABgsbUM7Iq4vUOezkj7bUHZQ0psW2zAAAM7b9LS0fXs2cn9kJBtXkPAVkI5cHgcAoHT6cJIbpjEFAPSn7dtfCeyaEyey8kQR2gCA/tSHk9wQ2gCA/tSHk9wQ2gCA/tSHk9wQ2gCA/tSHk9wwehwA0L/6bJIbetoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDQJlNT0tjY9KyZdnj9HSvW4QeYj1tACir6WlpYkI6cSLbPnw425b6ao1oFEdPGwDKavv2VwK75sSJrBwDidAGgLKamVlYOfoeoQ2gnLiXK42MLKwcfY/QBlA+tXu5hw9LEa/cyx204L7jDmnlyrPLVq7MyjGQCG0A5cO93Mz4uDQ1JY2OSnb2ODXFILQBRmgDKB/u5b5ifFw6dEg6cyZ7JLCXRiK3YwhtAOXDvVx0Qy2obemGG5K4HUNoAygf7uViqdWPm5CysK5XfzumRL1wQhtA+XAvF0ut2biJRjMzpRsU6Wj866IEKpVKVKvVXjcDANCvli07t3fdaHQ0e6z1xhv3HTrUsebY3hcRlXb16GkDAAZPu/ERtdsxJRsUSWgDAAZPs3ETdvZYfzumZIMi24a27Z22n7X9jRb7r7X9XdsP5z8fq9u30fYTtg/Y/mgnGw4AwKI1Gzfx+c9nl8zrv1pXskGRRVb5+qykHZI+N0+d/xsRP1tfYPsCSZ+U9A5Js5L22t4VEfsX2VYAADpnfLz94Mba/u3bs0viIyNZYPdoUGTb0I6Ir9oeW8R7Xy3pQEQclCTbd0naKonQBgCko0i4d0mn7mn/hO1/s73b9hvyskslHamrM5uXNWV7wnbVdnVubq5DzQIAoH90IrQfkjQaEW+S9CeS/j4vd5O6LcfXR8RURFQiojI8PNyBZgHAeSrRpBqA1IHQjojnI+L7+fN7JS23vVpZz3pdXdXLJB09398HAF1Rskk1AKkDoW37EjsbJ2/76vw9j0naK+kK25fbfpWkbZJ2ne/vA4CuYKUxlFDbgWi275R0raTVtmcl3SJpuSRFxKclvVvSB2yflvSfkrZFNs3aaduTku6TdIGknRHx2JIcBQB0Wskm1QCkYqPHr2+zf4eyr4Q123evpHsX1zQA6KGRkebTV7LSGHqIGdEAoJmSTaoBSIQ2ADTHSmMooSIzogHAYCrRpBqARE8bANAO31cvDXraAIDWat9Xr339rfZ9dYmrED1ATxsA0BrfVy8VQhsA0BrfVy8VQhsA0Fqr76XzffWeILQBoEzKNuiL76uXCqENAGXR6UVKOvEHAN9XLxVn04SXS6VSiWq12utmAEB3jY01nzp1dFQ6dGhh79U46lvKesgEbinZ3hcRlXb16GkDQFl0ctAXo777EqENAGXRyUFfjPruS4Q2AJRFJwd9Meq7LxHaAFAWnRz0xajvvsQ0pgBQJp1apKT2Htu3Z5fER0aywGYQWtIIbQDoV6xS1ne4PA4AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAGXSiTWw0beYEQ0AyqJxDezDh7NtiZnNIImeNgCUB2tgow1CGwDKgjWw0QahDQBlwRrYaIPQBoCyYA1stEFoA0BZjI9LU1PS6KhkZ49TUwxCw8sYPQ4AZcIa2JgHPW0AABLRNrRt77T9rO1vtNg/bvuR/Odrtt9Ut++Q7UdtP2y72smGAwAwaIr0tD8raeM8+78l6W0R8UZJt0maatj/9oh4c0RUFtdEAAAgFbinHRFftT02z/6v1W0+KOmy828WAABo1Ol72u+TtLtuOyTtsb3P9kSHfxcAAAOlY6Ft++3KQvu36oqviYirJG2SdJPtn5rn9RO2q7arc3NznWoWgKXAohZAT3QktG2/UdKfS9oaEcdq5RFxNH98VtIXJV3d6j0iYioiKhFRGR4e7kSzACyF2qIWhw9LEa8sakFwA0vuvEPb9oikv5N0Q0R8s6781bZ/sPZc0gZJTUegA0gIi1oAPdN2IJrtOyVdK2m17VlJt0haLkkR8WlJH5N0saQ/tS1Jp/OR4mskfTEvu1DSX0XEl5bgGAB0E4taAD1TZPT49W32v1/S+5uUH5T0pnNfASBpIyPZJfFm5QCWFDOiAVgYFrUAeobQBrAwLGoB9AwLhgBYOBa1AHqCnjYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoA0DKpqelsTFp2bLscXq61y3CEmKVLwBI1fS0NDEhnTiRbR8+nG1LrMLWp+hpA0Cqtm9/JbBrTpzIytGXCG0ASNXMzMLKkTxCGwBSNTKysHIkj9AGgFTdcYe0cuXZZStXZuXoS4Q2AKRqfFyampJGRyU7e5yaYhBaH2P0OACkbHyckB4g9LQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbQLkwlzbQEqPHAZQHc2kD86KnDaA8mEsbmBehDaA8mEsbmBehDaA8mEsbmBehDaA8mEsbmBehDaA8mEsbmBejxwGUC3NpAy0V6mnb3mn7WdvfaLHftv/Y9gHbj9i+qm7fjbafzH9u7FTDAQAYNEUvj39W0sZ59m+SdEX+MyHpU5Jk+yJJt0j6cUlXS7rF9qrFNhYAgEFWKLQj4quSjs9TZaukz0XmQUk/bPu1kq6TdH9EHI+I70i6X/OHPwAAaKFTA9EulXSkbns2L2tVDgAAFqhToe0mZTFP+blvYE/Yrtquzs3NdahZAAD0j06F9qykdXXbl0k6Ok/5OSJiKiIqEVEZHh7uULMAAOgfnQrtXZJ+KR9F/lZJ342Ib0u6T9IG26vyAWgb8jIAALBAhb6nbftOSddKWm17VtmI8OWSFBGflnSvpHdKOiDphKRfzvcdt32bpL35W90aEfMNaAMAAC0UCu2IuL7N/pB0U4t9OyXtXHjTAABAPaYxBQAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIRKHQtr3R9hO2D9j+aJP9f2D74fznm7b/vW7fS3X7dnWy8QAADJIL21WwfYGkT0p6h6RZSXtt74qI/bU6EfGbdfV/TdJb6t7iPyPizZ1rMgAAg6lIT/tqSQci4mBEvCDpLklb56l/vaQ7O9E4AADwiiKhfamkI3Xbs3nZOWyPSrpc0pfriodsV20/aPvnWv0S2xN5verc3FyBZgEAMFiKhLablEWLutsk3R0RL9WVjURERdL/kvSHtv9rsxdGxFREVCKiMjw8XKBZAAAMliKhPStpXd32ZZKOtqi7TQ2XxiPiaP54UNI/6ez73QAAoKAiob1X0hW2L7f9KmXBfM4ocNs/ImmVpH+pK1tle0X+fLWkayTtb3wtAABor+3o8Yg4bXtS0n2SLpC0MyIes32rpGpE1AL8ekl3RUT9pfMrJX3G9hllfyD8bv2ocwAAUJzPzthyqFQqUa1We90MAAC6wva+fPzXvJgRDQCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYVC2/ZG20/YPmD7o032v9f2nO2H85/31+270faT+c+NnWw8AACD5MJ2FWxfIOmTkt4haVbSXtu7ImJ/Q9UvRMRkw2svknSLpIqkkLQvf+13OtL6PrL72G7tOLpDz7z4jNYsX6PJtZPadPGmRdcDAPSfIj3tqyUdiIiDEfGCpLskbS34/tdJuj8ijudBfb+kjYtr6tLafWy3Nj+6WZWHKtr86GbtPra7q7/79pnb9fSLTysUevrFp3X7zO3ntKFoPQBAfyoS2pdKOlK3PZuXNfoF24/Yvtv2ugW+VrYnbFdtV+fm5go0q3N6HYY7ju7QyTh5VtnJOKkdR3csqh4AoD8VCW03KYuG7X+QNBYRb5T0j5L+cgGvzQojpiKiEhGV4eHhAs3qnF6H4TMvPlOovGg9AEB/KhLas5LW1W1fJulofYWIOBYRp/LNP5P0o0VfWwa9DsM1y9cUKi9aDwDQn4qE9l5JV9i+3ParJG2TtKu+gu3X1m1ukfR4/vw+SRtsr7K9StKGvKxUeh2Gk2snNeShs8qGPKTJtZOLqgcA6E9tR49HxGnbk8rC9gJJOyPiMdu3SqpGxC5Jv257i6TTko5Lem/+2uO2b1MW/JJ0a0QcX4LjOC+Tayd1+8ztZ10i72YY1kZ/txsVXrQeAKA/OaLpLeaeqlQqUa1Wu/o7+SoVAKBXbO+LiEq7em172oNi08WbCGkAQKkxjSkAAIkgtAfB9LQ0NiYtW5Y9Tk/3ukUAgEXg8ni/m56WJiakEyey7cOHs21JGh/vXbsAAAtGT7vfbd/+SmDXnDiRlQMAkkJo97uZmYWVAwBKi9DudyMjCysHAJQWod3v7rhDWrny7LKVK7NyAEBSCO0u6snyn+Pj0tSUNDoq2dnj1BSD0AAgQYwe75La8p+1qVJry39KWvpJXcbHCWkA6AMD09PuSS+3Tq+X/wQApG8gQrvWy336xacVipd7ud0M7if/+Uk9edOTev5fnz+rnLWwAQBFDURo97qXu2fPHj31waf0vQe/p6c++JSe/5dXgjv5tbCZbQ0AumYgQrtVb7Ybvdw9e/Zoy5YtOvPCGQ29bkjxQuipD2XBnfxa2LXZ1g4fliJemW2N4AaAJTEQod2qN7vUvdxaYJ86dUqr371aV951pVa/e7XihdDBDx3UO596Z9orizHbGgB01UCE9uTaSQ156Kyype7lNgb2ut9aJ19grfutdVr97tU688IZ3fae27Rnz54la8OSY7Y1AOiqgQjtTRdv0s0jN+uS5ZfIsi5ZfoluHrl5yXq5TQN7mSVJXvZKcJ86dUpbtmwpd3DPd8+a2dYAoKscEb1uwzkqlUpUq9VeN2NRHnjgAW3evLlpYNeLM6Ejv3dEz939nFasWKF77rlH69ev70GL59G4QpiUzaZWm5yl3X4AQCG290VEpV29gehpd9PHP/5xnTp1SkOvG9K6jzQPbCnvcX9knYZeN6RTp07pE5/4RJdbWkC7e9bMtgYAXUVod9hHPvIRrVixQicPntSRjx9RnGl+JSPOhI58/IhOHjypFStW6MMf/nCXW1pAkXvW4+PSoUPSmTPZI4ENAEuG0O6w9evXa9euXVqxYoWeu/s5Hfm9c4O78dL4rl27yndpXOKeNQCUDKG9BDZs2NAyuJsF9oYNG3rc4hZYIQwASoXQ7rDaHOe/s/p39MY/eqOWr1j+SnC/lFBgS9yzBoCSIbQ7qHGO85d+7CVd8ftXvBzcj297vKOB/cADD+i6667TAw880KEjaIJ71gBQGoR2BzWb43zorUN68x+9+eXBaZ0K7D179mjz5s1nPQIA+huh3UGt5jI/82NndM899+i6667TPffc05HArk3e8nopjUlaAADnjdDuoPnmOF+/fr2+9KUvnfco8frA/oCkRyR9QAQ3AAwCQruDlnqO88bA3iHpgvyR4AaA/kdod9BSznHeLLBrJ2+ZCG4AGATMPZ6A+vnMGwO73hlJk5I+JZV3PnMAwDmYe7yP1OYzf72kP1Hrk7Ys318bnFbK+cwBAItGaCegNp/5fkm/pqxH3cyZfP9+qbzzmQMAFo3Q7rb51qduoX4+808puwTeGNyNl8ZLO585AGDRCO1uqq0/ffiwFJE9TkwUCu76+cwbg7tZYJd6elQAwKIUCm3bG20/YfuA7Y822f8h2/ttP2L7Adujdftesv1w/rOrk41PTrv1qdtoFtwvicAGgEHRdvS47QskfVPSOyTNStor6fqI2F9X5+2Svh4RJ2x/QNK1EfGL+b7vR8RrFtKovh09vmxZ1sNuZGdzexfUOCNa7R42gQ0Aaerk6PGrJR2IiIMR8YKkuyRtra8QEV+JiFoX8kFJly20wQOhQ+tT1/e4CWwAGBxFQvtSSUfqtmfzslbeJ2l33faQ7artB23/3CLa2D86uD71hg0bOjqfOQCg/C4sUMdNyppeU7f9HkkVSW+rKx6JiKO2Xyfpy7YfjYinmrx2QtKEJI0ssOeZjNqyltu3SzMzWQ/7jjsWvdzl+vXrGSEOAAOkSGjPSlpXt32ZpKONlWz/jKTtkt4WEadq5RFxNH88aPufJL1F0jmhHRFTkqak7J528UNIzPg4a1IDABalyOXxvZKusH257VdJ2ibprFHgtt8i6TOStkTEs3Xlq2yvyJ+vlnSNsnFTAABggdr2tCPitO1JSfcpW1RqZ0Q8ZvtWSdWI2CXpE5JeI+lvbEvSTERskXSlpM/YPqPsD4TfrR91DgAAimPBkCZ2H9utHUd36JkXn9Ga5Ws0uXayIyt1JWF6umP33AEAxRT9yleRe9oDZfex3bp95nadjJOSpKdffFq3z9wuSf0f3LUZ22oTwNRmbJMIbgAoAaYxbbDj6I6XA7vmZJzUjqM7etSiLjrPGdsAAEuL0G7wzIvPLKi8r8zMLKwcANBVhHaDNcvXLKi8r3RoxjYAwNIgtBtMrp3UkIfOKhvykCbXTvaoRV3UwRnbAACdR2g32HTxJt08crMuWX6JLOuS5Zfo5pGb+38QmpQNNpuakkZHs0VMRkezbQahAUAp8JUvAAB6rJOrfA283cd2a/Ojm1V5qKLNj27W7mO7278IAIAO43vabQz097YBAKVCT7tOsx71QH9vGwBQKvS0c6161I2BXTMQ39sGAJQKPe1cqx71shb/iQbie9sAgFIhtHOtes5ndGZwv7cNACgVQjvXqudc+552V7+3PT0tjY1Jy5Zlj9PTS/e7AADJ4J52bnLt5Dn3sGs96k0Xb+reSHFW2gIAtEBPO1eamdBYaQsA0AIzopXNsmVSs3NiS2fOdL89AIAlx4xoqWKlLQBAC4R22bDSFgCghb4O7STnDC/LSluMYAeA0unbe9qNM5xJ2WjwgVlm83w0jmCXst4+y3QCwJIY+HvazBl+HhjBDgCl1Leh3WqGM+YML2BmZmHlAICu6NvQbjXDGXOGF8AIdgAopb4N7cm1k8wZvliMYAeAUurb0C7NDGcpKssIdgDAWfp29DgAAKkY+NHjAAD0G0IbAIBEENr9jpnNAKBvsJ52P2NtbgDoK/S0+xkzmwFAXyG0+xkzmwFAXyG0+xkzmwFAXyG0+xkzmwFAXykU2rY32n7C9gHbH22yf4XtL+T7v257rG7fb+flT9i+rnNNR1vMbAYAfaXt6HHbF0j6pKR3SJqVtNf2rojYX1ftfZK+ExH/zfY2Sb8n6Rdtv17SNklvkLRW0j/a/u8R8VKnDwQtjI8T0gDQJ4r0tK+WdCAiDkbEC5LukrS1oc5WSX+ZP79b0nrbzsvviohTEfEtSQfy9wMAAAtUJLQvlXSkbns2L2taJyJOS/qupIsLvhYAABRQJLTdpKxxlZFWdYq8NnsDe8J21XZ1bm6uQLMAABgsRUJ7VtK6uu3LJB1tVcf2hZJ+SNLxgq+VJEXEVERUIqIyPDxcrPUAAAyQIqG9V9IVti+3/SplA8t2NdTZJenG/Pm7JX05sjU/d0nalo8uv1zSFZL+tTNNBwBgsLQdPR4Rp21PSrpP0gWSdkbEY7ZvlVSNiF2S/kLS520fUNbD3pa/9jHbfy1pv6TTkm5i5DgAAIvjrENcLpVKJarVaq+bAQBAV9jeFxGVdvWYEQ0AgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCJK+ZUv23OSDp/n26yW9FwHmlMGHEs5cSzl0y/HIXEsZbVUxzIaEW2nAy1laHeC7WqR77ylgGMpJ46lfPrlOCSOpax6fSxcHgcAIBGENgAAiejn0J7qdQM6iGMpJ46lfPrlOCSOpax6eix9e08bAIB+0889bQAA+kqSoW17o+0nbB+w/dEm+1fY/kK+/+u2x+r2/XZe/oTt67rZ7kYFjuNDtvfbfsT2A7ZH6/a9ZPvh/KdxffOuK3As77U9V9fm99ftu9H2k/nPjY2v7bYCx/IHdcfxTdv/XrevbOdlp+1nbX+jxX7b/uP8WB+xfVXdvtKclwLHMZ63/xHbX7P9prp9h2w/mp+Tni8fWOBYrrX93bp/Rx+r2zfvv81uK3AsH647jm/kn4+L8n1lOy/rbH/F9uO2H7P9wSZ1ev95iYikfpSt6f2UpNdJepWkf5P0+oY6vyrp0/nzbZK+kD9/fV5/haTL8/e5oMTH8XZJK/PnH6gdR779/V6fiwUey3sl7Wjy2oskHcwfV+XPV5X5WBrq/5qyNeZLd17y9vyUpKskfaPF/ndK2i3Jkt4q6eslPS/tjuMna+2TtKl2HPn2IUmre30uFnAs10r6P03KF/RvswzH0lD3XZK+XOLq4IeiAAAEWUlEQVTz8lpJV+XPf1DSN5v8f6znn5cUe9pXSzoQEQcj4gVJd0na2lBnq6S/zJ/fLWm9befld0XEqYj4lqQD+fv1QtvjiIivRMSJfPNBSZd1uY1FFTknrVwn6f6IOB4R35F0v6SNS9TOIhZ6LNdLurMrLVuEiPiqpOPzVNkq6XOReVDSD9t+rUp2XtodR0R8LW+nVO7PSpFz0sr5fM6WxAKPpeyflW9HxEP58+9JelzSpQ3Vev55STG0L5V0pG57Vuf+h325TkSclvRdSRcXfG23LLQt71P2F17NkO2q7Qdt/9xSNHABih7LL+SXlO62vW6Br+2Wwu3Jb1dcLunLdcVlOi9FtDresp2XhWj8rISkPbb32Z7oUZsW6ids/5vt3bbfkJcle05sr1QWYn9bV1za8+LslupbJH29YVfPPy8XLsWbLjE3KWscAt+qTpHXdkvhtth+j6SKpLfVFY9ExFHbr5P0ZduPRsRTS9DOIoocyz9IujMiTtn+FWVXQn664Gu7aSHt2Sbp7oh4qa6sTOeliBQ+K4XZfruy0P6fdcXX5Ofkv0i63/b/y3uIZfWQsiktv2/7nZL+XtIVSvSc5N4l6Z8jor5XXsrzYvs1yv64+I2IeL5xd5OXdPXzkmJPe1bSurrtyyQdbVXH9oWSfkjZJZwir+2WQm2x/TOStkvaEhGnauURcTR/PCjpn5T9VdgrbY8lIo7Vtf/PJP1o0dd22ULas00Nl/tKdl6KaHW8ZTsvbdl+o6Q/l7Q1Io7VyuvOybOSvqje3RIrJCKej4jv58/vlbTc9moleE7qzPdZKc15sb1cWWBPR8TfNanS+89LL2/8L+ZH2dWBg8ouS9YGY7yhoc5NOnsg2l/nz9+gsweiHVTvBqIVOY63KBt4ckVD+SpJK/LnqyU9qR4OSCl4LK+te/7zkh7Mn18k6Vv5Ma3Kn19U5mPJ6/2IsoE0Lut5qWvXmFoPetqsswfW/GsZz0uB4xhRNkblJxvKXy3pB+uef03SxpKfk0tq/66UBdlMfn4K/dss07Hk+2udpleX+bzk/40/J+kP56nT889LcpfHI+K07UlJ9ykbTbkzIh6zfaukakTskvQXkj5v+4Cyfyzb8tc+ZvuvJe2XdFrSTXH2pc2yHccnJL1G0t9k4+g0ExFbJF0p6TO2zyi7WvK7EbG/F8chFT6WX7e9Rdl/9+PKRpMrIo7bvk3S3vztbo2zL6F1VcFjkbJBNXdF/onNleq8SJLtO5WNRl5te1bSLZKWS1JEfFrSvcpGxB6QdELSL+f7SnVeChzHx5SNW/nT/LNyOrJFHdZI+mJedqGkv4qIL3X9AOoUOJZ3S/qA7dOS/lPStvzfWdN/mz04hJcVOBYp+yN9T0T8R91LS3deJF0j6QZJj9p+OC/7HWV/EJbm88KMaAAAJCLFe9oAAAwkQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEvH/ASEBCt5fHl0zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a71c1b5fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for index_data, data in enumerate(dataset):\n",
    "    \n",
    "    # an array of all centroid distances per data point (to find min)\n",
    "    dist_vec = np.zeros((k,1))\n",
    "    \n",
    "    for index_centroids, centroid in enumerate(centroids):\n",
    "        \n",
    "        # populate distance vectors of all centroids for that data point\n",
    "        dist_vec[index_centroids] = eucl_distance(centroid, data)\n",
    "    \n",
    "    # find which centroid is closest to that data point (which cluster that datapoint belongs to)\n",
    "    belongs_to[index_data] = np.argmin(dist_vec) # argmin returns indices, not values\n",
    "        \n",
    "# plot clusters\n",
    "colors = ['red', 'limegreen']\n",
    "plt.figure(figsize=(8,7))\n",
    "\n",
    "for i, b in enumerate(belongs_to):\n",
    "    plt.scatter(dataset[i,0], dataset[i,1], color=colors[int(b)])\n",
    "for i, c in enumerate(centroids):\n",
    "    plt.scatter(centroids[i,0], centroids[i,1], color=colors[i], marker=\"D\", s=80, edgecolor='black', linewidth='2')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Update centroids\n",
    "\n",
    "Find the mean of each cluster and move its corresponding centroid to that point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Cluster 1 contains datapoints [3, 4, 5, 8, 9, 12, 13, 14, 16, 17, 19, 29].\n",
      "Its updated centroid position is at [0.20385818 0.46199778].\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAGfCAYAAAByXooxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X+QHOdd5/HPR7KyYhMKZGtP/qXdde58lJO7hJgpE/AVcRDRjyyRoEgV8i3G4ZLaIngJkCubcHLFlH9UBbvq4ECBZAFXSGqxE0xCLWcvllECuSI4aOVzrNg+J4qiX6XYXkshP9hIsazv/dE99uxoZqdXmp3pZ/b9qtqamaefmX3arfFnu/vbTzsiBAAAym9FtwcAAACKIbQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAibig2wNoZO3atTE8PNztYQAA0BF79+59ISIGWvUrZWgPDw9rZmam28MAAKAjbB8q0o/D4wAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABLRMrRtr7f9OdtP237S9m806GPbf2h7v+0nbF9ds+xG21/Nf25s9woAALBcFNnTPi3pv0fEVZLeLOkm26+r67NF0pX5z5ikP5Ek2xdKuk3Sj0u6RtJttte0aewAgDKbnJSGh6UVK7LHyclujyh5LUM7Ir4REY/lz78j6WlJl9V12ybp45F5VNIP275E0iZJj0TEiYj4pqRHJG1u6xoAAMpnclIaG5MOHZIissexMYL7PC3qnLbtYUlvkvTFukWXSTpS8/po3tasHQDQy3bskObm5rfNzWXtOGeFQ9v2ayT9taTfjIhv1y9u8JZYoL3R54/ZnrE9Mzs7W3RYAIAyOnx4ce0opFBo216lLLAnI+LTDboclbS+5vXlko4t0H6WiJiIiEpEVAYGBooMCwBQVoODi2tHIUWqxy3pzyU9HRH/s0m3KUm/nFeRv1nStyLiG5IelrTR9pq8AG1j3gYA6GV33SX1989v6+/P2nHOiuxpXyvpBkk/bfvx/Ofttn/V9q/mfR6SdEDSfkl/KunXJCkiTki6Q9Ke/Of2vA1AyU0fn9bIvhFVHqtoZN+Ipo9Pd3tISMnoqDQxIQ0NSbZ00UXSD/yAdMMNVJKfB0c0PMXcVZVKJWZmZro9DGDZmj4+rTsP36mTcfLlttVerVsHb9WWi7Z0cWRIUrWSvLYwrb8/C/XR0e6Nq0Rs742ISqt+zIgG4Cw7j+2cF9iSdDJOauexnV0aEZJGJXnbENoAzvLci88tqh1YEJXkbUNoAzjLulXrFtUOLIhK8rYhtAGcZfzSca326nltq71a45eOd2lESFrRSnKmPW3pgm4PAED5VIvNdh7bqedefE7rVq3T+KXjFKHh3FSLzXbsyA6JDw5mgV1bhFZfrFad9rT2/aB6HABQAsPDWVDXGxqSDh7s9Gg6jupxAEA6KFYrhNAGAHQfxWqFENoAgO5j2tNCCG0AQPfVT3s6NNR4xrRlXmFO9TgAoBxGRxeuFKfCnD1tAEAimA6V0AYAJIIKc0IbAJAIKswJbQBAIqgwJ7QBAIkoWmHew6geBwCko1WFeY9jTxsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaALAUJiel4WFpxYrscXKy2yNCD7igVQfb90r6WUnPR8R/arD8ZkmjNZ93laSBiDhh+6Ck70h6SdLpiKi0a+AAUFqTk9LYmDQ3l70+dCh7LUmjo83fB7TgiFi4g/1Tkr4r6eONQruu7zsk/VZE/HT++qCkSkS8sJhBVSqVmJmZWcxbAKA8hoezoK43NCQdPNjp0SABtvcW2bFteXg8Ij4v6UTB33u9pPsK9gWA3nT48OLagYLadk7bdr+kzZL+uqY5JO2yvdf2WIv3j9mesT0zOzvbrmEBQOcNDi6uHSionYVo75D0TxFRu1d+bURcLWmLpJvyQ+0NRcRERFQiojIwMNDGYQFAh911l9TfP7+tvz9rB85DO0N7u+oOjUfEsfzxeUmfkXRNG38fAJTT6Kg0MZGdw7azx4kJitBw3lpWjxdh+4ckvUXSL9W0vVrSioj4Tv58o6Tb2/H7AKD0RkcJabRdkUu+7pN0naS1to9Kuk3SKkmKiI/k3X5e0q6I+Leat66T9Bnb1d/zlxHxd+0bOgAAy0vL0I6I6wv0+Zikj9W1HZD0xnMdGAAA521yUtqxI6vcHxzM6goSPgLSlsPjAACUTg9OcsM0pgCA3rRjxyuBXTU3l7UnitAGAPSmHpzkhtAGAPSmHpzkhtAGAPSmHpzkhtAGAPSmHpzkhupxAEDv6rFJbtjTBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwDKbHJSGh6WVqzIHicnuz0idBH30waAspqclMbGpLm57PWhQ9lrqafuEY3i2NMGgLLaseOVwK6am8vasSwR2gBQVocPL64dPY/QBlBOnMuVBgcX146eR2gDKJ/qudxDh6SIV87lLrfgvusuqb9/flt/f9aOZYnQBlA+nMvNjI5KExPS0JBkZ48TExShLWOENoDy4VzuK0ZHpYMHpTNnskcCe2kkcjqG0AZQPpzLRSdUg9qWbrghidMxhDaA8uFcLpZabd2ElIV1rdrTMSXaCye0AZQP53Kx1BrVTdQ7fLh0RZGO+r8uSqBSqcTMzEy3hwEA6FUrVpy9d11vaCh7rO6N1y87eLBtw7G9NyIqrfqxpw0AWH5a1UdUT8eUrCiS0AYALD+N6ibs7LH2dEzJiiJbhrbte20/b/vLTZZfZ/tbth/Pfz5Ys2yz7Wds77f9gXYOHACAc9aobuITn8gOmddeWleyosgid/n6mKSdkj6+QJ//ExE/W9tge6WkD0t6m6SjkvbYnoqIp85xrAAAtM/oaOvixuryHTuyQ+KDg1lgd6kosmVoR8TnbQ+fw2dfI2l/RByQJNv3S9omidAGAKSjSLh3SLvOaf+E7S/Znrb9+rztMklHavoczdsasj1me8b2zOzsbJuGBQBA72hHaD8maSgi3ijpjyT9Td7uBn2b1tdHxEREVCKiMjAw0IZhAcB5KtGkGoDUhtCOiG9HxHfz5w9JWmV7rbI96/U1XS+XdOx8fx8AdETJJtUApDaEtu2L7axO3vY1+Wcel7RH0pW2r7D9KknbJU2d7+8DgI7gTmMooZaFaLbvk3SdpLW2j0q6TdIqSYqIj0h6p6T32j4t6XuStkc2zdpp2+OSHpa0UtK9EfHkkqwFALRbySbVAKRi1ePXt1i+U9klYY2WPSTpoXMbGgB00eBg4+krudMYuogZ0QCgkZJNqgFIhDYANMadxlBCRWZEA4DlqUSTagASe9oAgFa4Xr002NMGADRXvV69evlb9Xp1iaMQXcCeNgCgOa5XLxVCGwDQHNerlwqhDQBortl16Vyv3hWENgCUSdmKvrhevVQIbQAoi3bfpKQdfwBwvXqpOJsmvFwqlUrMzMx0exgA0FnDw42nTh0akg4eXNxn1Vd9S9keMoFbSrb3RkSlVT/2tAGgLNpZ9EXVd08itAGgLNpZ9EXVd08itAGgLNpZ9EXVd08itAGgLNpZ9EXVd09iGlMAKJN23aSk+hk7dmSHxAcHs8CmCC1phDYA9CruUtZzODwOAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBl0o57YKNnMSMaAJRF/T2wDx3KXkvMbAZJ7GkDQHlwD2y0QGgDQFlwD2y0QGgDQFlwD2y0QGgDQFlwD2y0QGgDQFmMjkoTE9LQkGRnjxMTFKHhZVSPA0CZcA9sLIA9bQAAEtEytG3fa/t5219usnzU9hP5zxdsv7Fm2UHb+2w/bnumnQMHAGC5KbKn/TFJmxdY/nVJb4mIN0i6Q9JE3fK3RsSPRkTl3IYIAACkAue0I+LztocXWP6FmpePSrr8/IcFAADqtfuc9rslTde8Dkm7bO+1Pdbm3wUAwLLSttC2/VZlof3bNc3XRsTVkrZIusn2Ty3w/jHbM7ZnZmdn2zUsAEuBm1oAXdGW0Lb9Bkl/JmlbRByvtkfEsfzxeUmfkXRNs8+IiImIqEREZWBgoB3DArAUqje1OHRIinjlphYEN7Dkzju0bQ9K+rSkGyLiKzXtr7b9g9XnkjZKaliBDiAh3NQC6JqWhWi275N0naS1to9Kuk3SKkmKiI9I+qCkiyT9sW1JOp1Xiq+T9Jm87QJJfxkRf7cE6wCgk7ipBdA1RarHr2+x/D2S3tOg/YCkN579DgBJGxzMDok3agewpJgRDcDicFMLoGsIbQCLw00tgK7hhiEAFo+bWgBdwZ42AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoASm337t3atGmTdu/e3e2hAF1HaAMorV27dmlkZGTeI7CcEdoASmnXrl3aunWrTp06pddJOnXqlLZu3UpwY1kjtAGUTm1gv1fSE5LeK4IbILQBlEp9YO+UtDJ/JLix3BHaAEqjUWBX/ye1QgQ3QGgDKIXdu3c3DeyqRsG97KvKJyel4WFpxYrscXKy2yPCEiK0AZTC3Xff/XLR2R+p+f+cVuTLq8Vp99xzT6eGWD6Tk9LYmHTokBSRPY6NEdw9jNAGUAq33HKL+vr69JSkX5d0pkm/M/nypyT19fXp5ptv7tQQy2fHDmlubn7b3FzWjp5EaAMohQ0bNmhqakp9fX36E0njOju4z+Ttf6IssKemprRhw4ZOD7U8Dh9eXDuSR2gDKI2NGzc2De5Ggb1x48ZuDbUcBgcX147kEdoASqVRcL8kAruhu+6S+vvnt/X3Z+3oSYQ2gNKpD+43iMBuaHRUmpiQhoYkO3ucmMja0ZMcEd0ew1kqlUrMzMx0exgAuqz2um0CG73M9t6IqLTqx542gNLauHGjHnzwQW3atEkPPvgggY1l74JuDwAAFrJhw4blXSEO1GBPGwCARBDaAAAkgtAGUC7MpQ00xTltAOVRnUu7OjVndS5ticuYALGnDaBMmEsbWBChDaA8mEsbWBChDaA8mEsbWBChDaA8mEsbWBChDaA8mEsbWBDV4wDKZXSUkAaaKLSnbfte28/b/nKT5bb9h7b3237C9tU1y260/dX858Z2DRwAgOWm6OHxj0navMDyLZKuzH/GlN1FT7YvlHSbpB+XdI2k22yvOdfBAgCwnBUK7Yj4vKQTC3TZJunjkXlU0g/bvkTSJkmPRMSJiPimpEe0cPgDAIAm2lWIdpmkIzWvj+ZtzdoBAMAitSu03aAtFmg/+wPsMdsztmdmZ2fbNCwAAHpHu0L7qKT1Na8vl3RsgfazRMRERFQiojIwMNCmYQEA0DvaFdpTkn45ryJ/s6RvRcQ3JD0saaPtNXkB2sa8DQAALFKh67Rt3yfpOklrbR9VVhG+SpIi4iOSHpL0dkn7Jc1J+pV82Qnbd0jak3/U7RGxUEEbAABoolBoR8T1LZaHpJuaLLtX0r2LHxoAAKjFNKYAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkoFNq2N9t+xvZ+2x9osPz3bT+e/3zF9r/WLHupZtlUOwcPAMByckGrDrZXSvqwpLdJOippj+2piHiq2icifqum/69LelPNR3wvIn60fUMGAGB5KrKnfY2k/RFxICK+L+l+SdsW6H+9pPvaMTgAAPCKIqF9maQjNa+P5m1nsT0k6QpJn61pXm17xvajtn+u2S+xPZb3m5mdnS0wLAAAlpcioe0GbdGk73ZJD0TESzVtgxFRkfRfJf2B7X/f6I0RMRERlYioDAwMFBgWAADLS5HQPippfc3ryyUda9J3u+oOjUfEsfzxgKR/0Pzz3QAAoKAiob1H0pW2r7D9KmXBfFYVuO0fkbRG0j/XtK2x3Zc/XyvpWklP1b8XAAC01rJ6PCJO2x6X9LCklZLujYgnbd8uaSYiqgF+vaT7I6L20PlVkj5q+4yyPxA+VFt1DgAAivP8jC2HSqUSMzMz3R4GAAAdYXtvXv+1IGZEAwAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIRKHQtr3Z9jO299v+QIPl77I9a/vx/Oc9NctutP3V/OfGdg4eAIDl5IJWHWyvlPRhSW+TdFTSHttTEfFUXddPRsR43XsvlHSbpIqkkLQ3f+832zL6HjJ9fFo7j+3Ucy8+p3Wr1mn80nFtuWjLOfcDAPSeInva10jaHxEHIuL7ku6XtK3g52+S9EhEnMiD+hFJm89tqEtr+vi0RvaNqPJYRSP7RjR9fLqjv/vOw3fq2RefVSj07IvP6s7Dd541hqL9AAC9qUhoXybpSM3ro3lbvV+w/YTtB2yvX+R7ZXvM9oztmdnZ2QLDap9uheHu3bu1adMm/e6nf1cn4+S8ZSfjpHYe2zmvbeexnYX6AQB6U5HQdoO2qHv9t5KGI+INkv5e0l8s4r1ZY8RERFQiojIwMFBgWO3TjTDctWuXRkZGtGvXLu0Z36Nv//O3z+rz3IvPLfi6VTsAoLcUCe2jktbXvL5c0rHaDhFxPCJO5S//VNKPFX1vGXQ6DHft2qWtW7fq1KlTWv3a1Yrvh772/q+dFdzrVq1b8HWrdgBAbykS2nskXWn7CtuvkrRd0lRtB9uX1LzcKunp/PnDkjbaXmN7jaSNeVupdDIMawN77TvX6qr7r9Lad649K7hXe7XGL51X16fxS8e12qvntTXqBwDoTS1DOyJOSxpXFrZPS/pURDxp+3bbW/Nu77P9pO0vSXqfpHfl7z0h6Q5lwb9H0u15W6l0KgzrA3v9b6+XV1rrf3v9vOBeuWelbh289ayq8C0XbdGtg7fq4lUXy7IuXnVxw34AgN7kiIanmLuqUqnEzMxMR3/nUl9K1TCwV7xyyj/OhI783hG98MAL6uvr09TUlDZu3Ni23w8AKC/beyOi0rIfob30du/erZGRkaaBXVUf3A8++KA2bNjQhREDADqpaGgzjWkH3H333S8Xna2/pXFgS5JXWOtvWa/Vr12tU6dO6Z577unwSAEAZUZod8Att9yivr4+nTxwUkfuPqI40/joRpwJHbn7iE4eOKm+vj7dfPPN7RnA5KQ0PCytWJE9Tk6253MBAB1FaHfAhg0bNDU1pb6+Pr3wwAs68ntnB3ejc9ptOTQ+OSmNjUmHDkkR2ePYGMENAAkitDtk48aNTYN7SYvQduyQ5ubmt83NZe0AgKQQ2h3UMLhfWuKq8cOHF9cOACgtQrvD6oP76e1PL+1lXoODi2sHAJQWod0FtcFdLTpbsuuy77pL6u+f39bfn7UDAJLS8n7aaJ95E7hcsk633Xeb/vGj/6ibb7556a7HHh3NHnfsyA6JDw5mgV1tBwAkg8lVOqR6+8/au4mt9mqmIQUAMLlKvenj0xrZN6LKYxWN7BtZ8ntl1+Ne2ACA87UsDo/X7+U+++KzuvPwnZLUsb1c7oUNADhfy2JPuwx7uT17L2xmWwOAjlkWoV2GvdyevBc2s60BQEcti9Auw15uT94Lm9nWAKCjlsU57fFLxxtWbndqL7f+Xt13DN2RdlhXMdsaAHTUstjT7uZebrUI7tkXn1UoXi6C63T1+jlb6Jw1s60BQEdxnfYSG9k3omdffPas9otXXawH//ODXRjRIlTPWdceAu/vlyYmsslZWi0HABTCddolUYYiuHPW6pz16GgW0ENDkp09EtgAsGSWxTntblq3al3DPe0kLvUqcs56dJSQBoAOYU97iSV9qRfnrAGgVAjtJZb0pV7cIQwASoXD421Wf3nX+KXj2nLRljRCuh53CAOAUiG026gMc5y3HeesAaA0ODzeRmWY4xwA0LsI7TZK+vIuAEDpEdptVIY5zgEAvYvQbqOOXt7FLTEBYNmhEK2NqsVmjarH26p++tDqLTElisYAoIcx93iKhoezoK43NCQdPNjp0QAAzhNzj/cybokJAMsSoZ0iphcFgGWJ0O60dhSQMb0oACxLhHYnVQvIDh2SIl4pIFtscHNLTABYlgoVotneLOl/SVop6c8i4kN1y98v6T2STkualfTfIuJQvuwlSfvyrocjYmur39ezhWgUkAEAGihaiNbyki/bKyV9WNLbJB2VtMf2VEQ8VdPt/0qqRMSc7fdKulvSL+bLvhcRP7roNehFFJABAM5DkcPj10jaHxEHIuL7ku6XtK22Q0R8LiLyi4b1qKTL2zvMHkEBGQDgPBQJ7cskHal5fTRva+bdkqZrXq+2PWP7Uds/dw5j7B0UkAEAzkORGdHcoK3hiXDbvySpIuktNc2DEXHM9mslfdb2voj4WoP3jkkak6TBXt3z5P7UAIDzUCS0j0paX/P6cknH6jvZ/hlJOyS9JSJOVdsj4lj+eMD2P0h6k6SzQjsiJiRNSFkhWvFVSAz3pwYAnKMih8f3SLrS9hW2XyVpu6Sp2g623yTpo5K2RsTzNe1rbPflz9dKulZSbQEbAAAoqOWedkSctj0u6WFll3zdGxFP2r5d0kxETEm6R9JrJP2VbemVS7uukvRR22eU/YHwobqqcwAAUBA3DGlg+vj00t+pq6wmJznnDgAd1rbrtJeb6ePTuvPwnToZJyVJz774rO48fKck9X5wc8tPACg1pjGts/PYzpcDu+pknNTOYzu7NKIO2rHjlcCumpvL2gEAXUdo13nuxecW1d5TmLENAEqN0K6zbtW6RbX3FGZsA4BSI7TrjF86rtVePa9ttVdr/NLxLo2og5ixDQBKjdCus+WiLbp18FZdvOpiWdbFqy7WrYO39n4RmsQtPwGg5LjkCwCALit6yRd72gVMH5/WyL4RVR6raGTfiKaPT7d+EwAAbcZ12i0s6+u2AQClwp52jUZ71Mv6um0AQKmwp51rtkddH9hVy+K6bQBAqbCnnWu2R72iyX+iZXHdNgCgVAjtXLM95zM6s3yv2wYAlAqhnWu251y9Truj121PTkrDw9KKFdnj5OTS/S4AQDI4p50bv3T8rHPY1T3qLRdt6VylOHfaAgA0wZ52rjQzoXGnLQBAE8yIVjYrVkiNtoktnTnT+fEAAJYcM6KlijttAQCaILTLhjttAQCa6OnQTnLO8LLcaYsKdgAonZ49p10/w5mUVYMvm9tsno/6CnYp29vnNp0AsCSW/Tlt5gw/D1SwA0Ap9WxoN5vhjDnDCzh8eHHtAICO6NnQbjbDGXOGF0AFOwCUUs+G9vil48wZfq6oYAeAUurZ0C7NDGcpKksFOwBgnp6tHgcAIBXLvnocAIBeQ2gDAJAIQrvXMbMZAPQM7qfdy7g3NwD0FPa0exkzmwFATyG0exkzmwFATyG0exkzmwFATyG0exkzmwFATykU2rY3237G9n7bH2iwvM/2J/PlX7Q9XLPsd/L2Z2xvat/Q0RIzmwFAT2lZPW57paQPS3qbpKOS9tieioinarq9W9I3I+I/2N4u6fck/aLt10naLun1ki6V9Pe2/2NEvNTuFUETo6OENAD0iCJ72tdI2h8RByLi+5Lul7Strs82SX+RP39A0gbbztvvj4hTEfF1SfvzzwMAAItUJLQvk3Sk5vXRvK1hn4g4Lelbki4q+F4AAFBAkdB2g7b6u4w061PkvdkH2GO2Z2zPzM7OFhgWAADLS5HQPippfc3ryyUda9bH9gWSfkjSiYLvlSRFxEREVCKiMjAwUGz0AAAsI0VCe4+kK21fYftVygrLpur6TEm6MX/+Tkmfjeyen1OStufV5VdIulLSv7Rn6AAALC8tq8cj4rTtcUkPS1op6d6IeNL27ZJmImJK0p9L+oTt/cr2sLfn733S9qckPSXptKSbqBwHAODcONshLpdKpRIzMzPdHgYAAB1he29EVFr1Y0Y0AAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJKOUlX7ZnJR06z49ZK+mFNgynDFiXcmJdyqdX1kNiXcpqqdZlKCJaTgdaytBuB9szRa55SwHrUk6sS/n0ynpIrEtZdXtdODwOAEAiCG0AABLRy6E90e0BtBHrUk6sS/n0ynpIrEtZdXVdevacNgAAvaaX97QBAOgpSYa27c22n7G93/YHGizvs/3JfPkXbQ/XLPudvP0Z25s6Oe56Bdbj/bafsv2E7d22h2qWvWT78fyn/v7mHVdgXd5le7ZmzO+pWXaj7a/mPzfWv7fTCqzL79esx1ds/2vNsrJtl3ttP2/7y02W2/Yf5uv6hO2ra5aVZrsUWI/RfPxP2P6C7TfWLDtoe1++Tbp++8AC63Kd7W/V/Dv6YM2yBf9tdlqBdbm5Zj2+nH8/LsyXlW27rLf9OdtP237S9m806NP970tEJPWj7J7eX5P0WkmvkvQlSa+r6/Nrkj6SP98u6ZP589fl/fskXZF/zsoSr8dbJfXnz99bXY/89Xe7vS0WuS7vkrSzwXsvlHQgf1yTP19T5nWp6//ryu4xX7rtko/npyRdLenLTZa/XdK0JEt6s6QvlnS7tFrgWaU7AAAEf0lEQVSPn6yOT9KW6nrkrw9KWtvtbbGIdblO0v9u0L6of5tlWJe6vu+Q9NkSb5dLJF2dP/9BSV9p8P+xrn9fUtzTvkbS/og4EBHfl3S/pG11fbZJ+ov8+QOSNth23n5/RJyKiK9L2p9/Xje0XI+I+FxEzOUvH5V0eYfHWFSRbdLMJkmPRMSJiPimpEckbV6icRax2HW5XtJ9HRnZOYiIz0s6sUCXbZI+HplHJf2w7UtUsu3Saj0i4gv5OKVyf1eKbJNmzud7tiQWuS5l/658IyIey59/R9LTki6r69b170uKoX2ZpCM1r4/q7P+wL/eJiNOSviXpooLv7ZTFjuXdyv7Cq1pte8b2o7Z/bikGuAhF1+UX8kNKD9hev8j3dkrh8eSnK66Q9Nma5jJtlyKarW/Ztsti1H9XQtIu23ttj3VpTIv1E7a/ZHva9uvztmS3ie1+ZSH21zXNpd0uzk6pvknSF+sWdf37csFSfOgSc4O2+hL4Zn2KvLdTCo/F9i9Jqkh6S03zYEQcs/1aSZ+1vS8ivrYE4yyiyLr8raT7IuKU7V9VdiTkpwu+t5MWM57tkh6IiJdq2sq0XYpI4btSmO23Kgvt/1LTfG2+Tf6dpEds/798D7GsHlM2peV3bb9d0t9IulKJbpPcOyT9U0TU7pWXcrvYfo2yPy5+MyK+Xb+4wVs6+n1JcU/7qKT1Na8vl3SsWR/bF0j6IWWHcIq8t1MKjcX2z0jaIWlrRJyqtkfEsfzxgKR/UPZXYbe0XJeIOF4z/j+V9GNF39thixnPdtUd7ivZdimi2fqWbbu0ZPsNkv5M0raIOF5tr9kmz0v6jLp3SqyQiPh2RHw3f/6QpFW21yrBbVJjoe9KabaL7VXKAnsyIj7doEv3vy/dPPF/Lj/Kjg4cUHZYslqM8fq6PjdpfiHap/Lnr9f8QrQD6l4hWpH1eJOywpMr69rXSOrLn6+V9FV1sSCl4LpcUvP85yU9mj+/UNLX83Vakz+/sMzrkvf7EWWFNC7rdqkZ17CaFz2NaH5hzb+UcbsUWI9BZTUqP1nX/mpJP1jz/AuSNpd8m1xc/XelLMgO59un0L/NMq1Lvry60/TqMm+X/L/xxyX9wQJ9uv59Se7weESctj0u6WFl1ZT3RsSTtm+XNBMRU5L+XNInbO9X9o9le/7eJ21/StJTkk5LuinmH9os23rcI+k1kv4qq6PT4YjYKukqSR+1fUbZ0ZIPRcRT3VgPqfC6vM/2VmX/3U8oqyZXRJywfYekPfnH3R7zD6F1VMF1kbKimvsj/8bmSrVdJMn2fcqqkdfaPirpNkmrJCkiPiLpIWUVsfslzUn6lXxZqbZLgfX4oLK6lT/OvyunI7upwzpJn8nbLpD0lxHxdx1fgRoF1uWdkt5r+7Sk70nanv87a/hvswur8LIC6yJlf6Tvioh/q3lr6baLpGsl3SBpn+3H87b/oewPwtJ8X5gRDQCARKR4ThsAgGWJ0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARPx/cmVuzVEJo5IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a71a6df080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# relocate centroids to the mean of their respective clusters\n",
    "# find the indices of the datapoits that belong to each cluster\n",
    "for index in range(len(centroids)):\n",
    "    instances = [i for i in range(len(belongs_to)) if belongs_to[i] == index]\n",
    "    centroids[index] = np.mean(dataset[instances], axis=0)\n",
    "    \n",
    "print('\\nCluster {} contains datapoints {}.\\nIts updated centroid position is at {}.'\n",
    "      .format(index, instances, centroids[index]))\n",
    "\n",
    "# plot updated centroids\n",
    "colors = ['red', 'limegreen']\n",
    "plt.figure(figsize=(8,7))\n",
    "\n",
    "for i, b in enumerate(belongs_to):\n",
    "    plt.scatter(dataset[i,0], dataset[i,1], color=colors[int(b)])\n",
    "for i, c in enumerate(centroids):\n",
    "    plt.scatter(centroids[i,0], centroids[i,1], color=colors[i], marker=\"D\", s=80, edgecolor='black', linewidth='2')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Repeat\n",
    "\n",
    "Repeat steps 2 and 3 until centroid stabilization is achieved. Namely, when all ceontroids stay in place after an update. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAMbCAYAAABg4M3WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3X2QXNV95//PZ4QYVn6IZTQBAxoNzo9fCrx+wOnFSeyKcZSMpGgtkoqrVuyE4MQuVbxWnrwF8a4oU4VRrQNVSTbBMZ4k+jl2TSAJCanJTwhEZGdJ4uCl5SVgINiyoqeSMQNy/BBZMkLf/ePegZ6e7uk7Uk/3Pbffr6qpnj739My56hl95tz7vec6IgQAAMpvqN8DAAAAxRDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgESc0+8BtLJq1aoYGxvr9zAAAOiJvXv3PhcRI536lTK0x8bGVK/X+z0MAAB6wvbBIv04PA4AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkomNo215t+3O2n7L9hO1fbdHHtn/X9j7bj9l+a8O2621/Jf+4vts7AAAoqakpaWxMGhrKHqem+j2i5BW5NecpSf81Ir5o+1WS9tp+MCKebOizQdJl+cfbJH1C0ttsv1bSzZJqkiJ/7XREfKOrewEAKJepKWnLFun48ez5wYPZc0mamOjfuBLXcaYdEV+LiC/mn39b0lOSLm7qdo2kT0fmYUmvsf06SeskPRgRx/KgflDS+q7uAQCgfLZtezmwZx0/nrXjjC3qnLbtMUlXSvpC06aLJR1ueH4kb2vX3uprb7Fdt12fmZlZzLAAAGVz6NDi2lFI4dC2/UpJfyHp1yLiW82bW7wkFmif3xgxGRG1iKiNjIwUHRYAoIxGRxfXjkIKhbbt5coCeyoi/rJFlyOSVjc8v0TS0QXaAQBVtn27tGLF3DY7O7dNUdoZK1I9bkl/JOmpiPitNt2mJf18XkX+w5K+GRFfk/SApHHbK22vlDSetwEoOyp/cTYmJqTJSWnNmuy5LUV+oHW2KI2fqUVzRMuj1S93sN8h6e8kPS7pdN783yWNSlJE3JkH+x3KisyOS/qFiKjnr//FvL8kbY+I/6/ToGq1WtTr9cXvDYDuaK78lbJZ0+Qklb9YvLGxLKibrVkjHTjQ69GUku29EVHr2K9TaPcDoQ30Gf/JopuGhl6eZTeypdOn57cPoKKhzYpoAOaj8hfdRFFa1xDaAObjP1l0U6uitBUrsvZG1FF0RGgDmK/of7JAEY1FaXb22FwfMVtHcfBgdiidYrWWOKcNoLWpqWz1qkOHshn29u0UoWHpDHgdBYVoAIB0DHixGoVoAIB0UEdRCKENAOg/itUKIbQBAP1HsVohnNMGAKShwsVqnNMGAFQLi/4Q2gCARFCsRmgDABLBoj+ENgAgEUWK1SrunH4PAACAwiYmBiqkmzHTBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaALAUpqaksTFpaCh7nJrq94hQAed06mB7h6T/KOnZiPj3LbbfIGn2juTnSLpc0khEHLN9QNK3Jb0o6VRE1Lo1cAAorakpacsW6fjx7PnBg9lzSZqYaP86oANHxMId7B+T9B1Jn24V2k193y3p1yPix/PnByTVIuK5xQyqVqtFvV5fzEsAoDzGxrKgbrZmjXTgQK9HgwTY3ltkYtvx8HhEPCTpWMHve62kuwr2BYBqOnRoce1AQV07p217haT1kv6ioTkk7ba91/aWDq/fYrtuuz4zM9OtYQFA742OLq4dKKibhWjvlvQPEdE4K397RLxV0gZJH8wPtbcUEZMRUYuI2sjISBeHBQA9tn27tGLF3LYVK7J24Cx0M7Q3q+nQeEQczR+flXSvpKu6+P0AoJwmJqTJyewctp09Tk5ShIaz1rF6vAjb3yfpnZJ+rqHtFZKGIuLb+efjkm7pxvcDgNKbmCCk0XVFLvm6S9LVklbZPiLpZknLJSki7sy7/Yyk3RHxbw0vvUDSvbZnv8+fRMT93Rs6AACDpWNoR8S1Bfp8StKnmtr2S3rzmQ4MAICzNjUlbduWVe6PjmZ1BQkfAenK4XEAAEqngovcsIwpAKCatm17ObBnHT+etSeK0AYAVFMFF7khtAEA1VTBRW4IbQBANVVwkRtCGwBQTRVc5IbqcQBAdVVskRtm2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AaDMpqaksTFpaCh7nJrq94jQR9wwBADKampK2rJFOn48e37wYPZcqtRNMFAcM20AKKtt214O7FnHj2ftGEiENgCU1aFDi2tH5RHaAMqJc7nS6Oji2lF5hDaA8pk9l3vwoBTx8rncQQvu7dulFSvmtq1YkbVjIBHaAMqHc7mZiQlpclJas0ays8fJSYrQBhihDaB8OJf7sokJ6cAB6fTp7JHAXhqJnI4htAGUD+dy0QuzQW1L113X/nRMiQKd67QBlM/27XOvT5Y4l4vuar4GPmLu9sbTMSW6Vp6ZNoDy4Vwullqruolmhw6Vrr7C0fzXRQnUarWo1+v9HgYAoKqGhubPrputWZMFd6t+dlZn0CW290ZErVM/ZtoAgMHTqT5i9nRMyeorOoa27R22n7X9pTbbr7b9TduP5h8fadi23vbTtvfZ/nA3Bw4AwBlrdQ28nT02no4p2bXyRWban5K0vkOfv4uIt+Qft0iS7WWSPi5pg6QrJF1r+4qzGSwAAF3Rqm7iM5/JDoU3XlpXsvqKjtXjEfGQ7bEz+NpXSdoXEfslyfbdkq6R9OQZfC0AALprYqJY+Bbt1wPdOqf9I7b/yfYu22/I2y6WdLihz5G8rSXbW2zXbddnZma6NCwAOAsluj4XkLoT2l+UtCYi3izp9yT9Vd7uFn3blupFxGRE1CKiNjIy0oVhAcBZYP1zlNBZh3ZEfCsivpN/fp+k5bZXKZtZr27oeomko2f7/QCgJ0p2fS4gdSG0bV9oZyV3tq/Kv+bzkh6RdJntS22fK2mzpOmz/X4A0BOsf44S6liIZvsuSVdLWmX7iKSbJS2XpIi4U9J7JH3A9ilJ35W0ObIVW07Z3irpAUnLJO2IiCeWZC8AoNtGR7ND4q3agT4pUj1+bYftd0i6o822+yTdd2ZDA4A+Yv1zlBArogFAKyW7PheQuMsXALRXoutzAYmZNgCgE65XLw1m2gCA9prvO93n+0kPOmbaAID2uF69VAhtAEB7XK9eKoQ2AJRJ2c4fl+x+0oOO0AaAsijjeuclu5/0oCO0AaAsun3+uBuzdq5XLxVnK46WS61Wi3q93u9hAEBvDQ1lM+xmtnT69OK+VnPVt5TNkAncUrK9NyJqnfox0waAsujm+WOqviuJ0AaAsujm+WOqviuJ0AaAsujm+WOqviuJ0AaAMpmYkA4cyM5hHzhw5uefqfquJEIbAKqIqu9KYu1xAKgq7lJWOcy0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGgDIp2/20USpcpw0AZdF8Z67Z+2lLXG8NScy0AaA8uDMXOiC0AaAsuDMXOiC0AaAsuDMXOiC0AaAsuDMXOiC0AaAsuDMXOqB6HADKhDtzYQHMtAEASETH0La9w/aztr/UZvuE7cfyj8/bfnPDtgO2H7f9qO16NwcOAMCgKTLT/pSk9Qts/xdJ74yIN0n6qKTJpu3vioi3RETtzIYIAACkAue0I+Ih22MLbP98w9OHJV1y9sMCAADNun1O+32SdjU8D0m7be+1vWWhF9reYrtuuz4zM9PlYQHoKtbHBvqia9Xjtt+lLLTf0dD89og4avv7JT1o+58j4qFWr4+ISeWH1mu1WnRrXAC6jPWxgb7pykzb9psk/aGkayLi+dn2iDiaPz4r6V5JV3Xj+wHoI9bHBvrmrEPb9qikv5R0XUR8uaH9FbZfNfu5pHFJLSvQASSE9bGBvul4eNz2XZKulrTK9hFJN0taLkkRcaekj0g6X9Lv25akU3ml+AWS7s3bzpH0JxFx/xLsA4BeGh3NDom3agewpIpUj1/bYfv7Jb2/Rft+SW+e/woASdu+fe45bYn1sYEeYUU0AIvD+thA37D2OIDFY31soC+YaQMAkAhCGwCARBDaAAAkgtAGACARhDaAUtuzZ4/WrVunPXv29HsoQN8R2gBKa/fu3dq4ceOcR2CQEdoASmn37t3atGmTTp48qSsknTx5Ups2bSK4MdAIbQCl0xjYH5D0mKQPiOAGCG0ApdIc2HdIWpY/EtwYdIQ2gNJoFdiz/0kNieAGCG0ApbBnz562gT2rVXAPfFX51JQ0NiYNDWWPU1P9HhGWEKENoBRuu+22l4rOfk/t/3MayrfPFqfdfvvtvRpi+UxNZXdcO3hQisget2whuCuM0AZQCjfeeKOGh4f1pKRflnS6Tb/T+fYnJQ0PD+uGG27o1RDLZ9u2ubdIlbLn27b1ZzxYcoQ2gFJYu3atpqenNTw8rE9I2qr5wX06b/+EssCenp7W2rVrez3U8jh0aHHtSB6hDaA0xsfH2wZ3q8AeHx/v11DLYXR0ce1IHqENoFRaBfeLIrBb2r5dWrFibtuKFVk7KonQBlA6zcH9JhHYLU1MSJOT0po1kp09Tk5m7agkR0S/xzBPrVaLer3e72EA6LPG67YJbFSZ7b0RUevUj5k2gNIaHx/Xzp07tW7dOu3cuZPAxsA7p98DAICFrF27drArxIEGzLQBlAsrfAFtMdMGUB6zK3zNLhgyu8KXRHEVIGbaAMqEFb6ABRHaAMqDFb6ABRHaAMqDFb6ABRHaAMqDFb6ABRHaAMqDFb6ABVE9DqBcJiYIaaANZtoAACSiUGjb3mH7WdtfarPdtn/X9j7bj9l+a8O2621/Jf+4vlsDBwBg0BSdaX9K0voFtm+QdFn+sUXZDXlk+7WSbpb0NklXSbrZ9sozHSwAAIOsUGhHxEOSji3Q5RpJn47Mw5JeY/t1ktZJejAijkXENyQ9qIXDHwAAtNGtc9oXSzrc8PxI3taufR7bW2zXbddnZma6NCwAAKqjW6HtFm2xQPv8xojJiKhFRG1kZKRLwwIAoDq6FdpHJK1ueH6JpKMLtAMAgEXqVmhPS/r5vIr8hyV9MyK+JukBSeO2V+YFaON5GwAAWKRCi6vYvkvS1ZJW2T6irCJ8uSRFxJ2S7pP0U5L2STou6Rfybcdsf1TSI/mXuiUiFipoAwAAbRQK7Yi4tsP2kPTBNtt2SNqx+KEBAIBGrIgGAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhULb9nrbT9veZ/vDLbb/tu1H848v2/7Xhm0vNmyb7ubgAQAYJOd06mB7maSPS/pJSUckPWJ7OiKenO0TEb/e0P+XJV3Z8CW+GxFv6d6QAQAYTEVm2ldJ2hcR+yPie5LulnTNAv2vlXRXNwYHAABeViS0L5Z0uOH5kbxtHttrJF0q6bMNzefZrtt+2PZPt/smtrfk/eozMzMFhgUAwGApEtpu0RZt+m6WdE9EvNjQNhoRNUn/WdLv2P6BVi+MiMmIqEVEbWRkpMCwAAAYLEVC+4ik1Q3PL5F0tE3fzWo6NB4RR/PH/ZL+VnPPdwMAgIKKhPYjki6zfantc5UF87wqcNs/KGmlpH9saFtpezj/fJWkt0t6svm1AACgs47V4xFxyvZWSQ9IWiZpR0Q8YfsWSfWImA3wayXdHRGNh84vl/RJ26eV/YHwscaqcwAAUJznZmw51Gq1qNfr/R4GAAA9YXtvXv+1IFZEAwAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEFApt2+ttP217n+0Pt9j+Xtszth/NP97fsO1621/JP67v5uAH0a7nd2nj4xtV+2JNGx/fqF3P7+r3kAAAPXJOpw62l0n6uKSflHRE0iO2pyPiyaaufxoRW5te+1pJN0uqSQpJe/PXfqMro6+QXc/v0h1H79DXX/i6Llh+gbZetFUbzt8wr8+th27ViTghSXrmhWd066FbJWleXwBA9RSZaV8laV9E7I+I70m6W9I1Bb/+OkkPRsSxPKgflLT+zIa6tPoxg92zZ4/WrVun/3Hv/9Cth27VMy88o1C8FMbNY7jj6B0vBfasE3FCdxy9Y8nHCgDovyKhfbGkww3Pj+RtzX7W9mO277G9epGvle0ttuu26zMzMwWG1T2zM9hOodlNu3fv1saNG7V7927dtPkmPfv5Z+dsbxXGX3/h6y2/Vrt2AEC1FAltt2iLpud/LWksIt4k6W8k/fEiXps1RkxGRC0iaiMjIwWG1T29nsHu3r1bmzZt0smTJ3Xe68/T6e+d1lc/9FV96x+/NadfcxhfsPyCll+vXTsAoFqKhPYRSasbnl8i6Whjh4h4PiJO5k//QNIPFX1tGfRyBtsY2Kves0qX3325Vr1nleJ7MS+4m8N460VbdZ7Pm9N2ns/T1ovmlBIAACqqSGg/Iuky25faPlfSZknTjR1sv67h6SZJT+WfPyBp3PZK2ysljedtpdKrGWxzYK/+jdXyMmv1b6yeF9ytwnjD+Rt00+hNunD5hbKsC5dfqJtGb6IIDQAGRMfq8Yg4ZXursrBdJmlHRDxh+xZJ9YiYlvQrtjdJOiXpmKT35q89ZvujyoJfkm6JiGNLsB9nZetFW+dUZUvdn8G2DOyh7OyBh7LglqTn7nlO+z+0X++deq82XDk/jDecv4GQBoAB5YiWp5j7qlarRb1e7+n3LHLJ1Znas2ePNm7c2DKwG8Xp0OHfPKzn7nlOw8PD2rlzp9auXduVMQAAysv23oiodexHaC+9devWaffu3Trv9efp8rsvl5e1qs/LxIuhpzY/pRP7T2jdunW6//77ezhSAEA/FA1tljHtgRtvvFHDw8M6sf+EDt92WHG69R9KcTp0+LbDOrH/hIaHh3XDDTd05fuzihoAVAOh3QNr167V9PS0hoeH9dw9z+nwb84P7uZD49PT0105NN6Pa9ABAEuD0O6R8fHxtsHdKrDHx8e78n1ZRQ0AqoPQ7qGWwf3i0gW2xCpqAFAlhHaPNQf3U5ufWrLAllhFDQCqhNDug8bgni06W4rAllhFDQCqpOPiKuieOdeCv+4C3XzXzfpfn/xfuuGGG5bseuzZa82X6hp0AEDvDMx12ku5eErR799q1TWWIQUAcJ12gzJc9kQVNwDgbA1EaJchMKtaxc3CLQDQOwMR2mUIzCpWcZfhCAYADJKBCO0yBGYVq7jLcAQDAAbJQIR2vwNztgjuRJzQUP5PXoV7YZfhCAYADJKBuOSrn5c9NVeNn9bpl/5gSCGwF6q6v2D5BXrmhWfmvSblQ/4AUGYDc8lXv2x8fGPLYLtw+YXa+cadfRhRcZ0uU+MyNgDoDi75KomUDyF3Ome94fwNumn0Jl24/EJZrsQhfwAos4E4PN5PKR9CLvIHx4bzNxDSANAjzLSXWL+L4M5GGaruAQAvI7S7rHmxEUnJHkJO+Q8OAKgiDo93UXNh1uxiIzeN3lT6orNWuNkIAJQL1eNdlHKlOACgf6ge74OUK8UBAOVHaHcRhVsAgKVEaHdRLwu3uLsWAAweCtG6qFeFW+0K3hrHAACoHgrREkTBGwBUC4VoFUbBGwAMJkK7x7pxLpqCNwAYTIR2D82ei37mhWcUipfORS82uFmpDAAGE6HdQ53umlUUd9cCgMFUqHrc9npJ/1PSMkl/GBEfa9r+IUnvl3RK0oykX4yIg/m2FyU9nnc9FBGbujT25HTzXDR31wKAwdNxpm17maSPS9og6QpJ19q+oqnb/5FUi4g3SbpH0m0N274bEW/JPwY2sCXORQMAzk6Rw+NXSdoXEfsj4nuS7pZ0TWOHiPhcRBzPnz4s6ZLuDrMaOBcNADgbRUL7YkmHG54fydvaeZ+kxsqq82zXbT9s+6fbvcj2lrxffWZmpsCw0sO5aADA2ShyTtst2lquyGL75yTVJL2zoXk0Io7afr2kz9p+PCK+Ou8LRkxKmpSyxVUKjCtJnIsGAJypIjPtI5JWNzy/RNLR5k62f0LSNkmbIuLkbHtEHM0f90v6W0lXnsV4AQAYWEVC+xFJl9m+1Pa5kjZLmm7sYPtKSZ9UFtjPNrSvtD2cf75K0tslPdmtwaP7uBEJAJRXx8PjEXHK9lZJDyi75GtHRDxh+xZJ9YiYlnS7pFdK+nPb0suXdl0u6ZO2Tyv7A+FjEVH60N71/K4lv+lHGXEjEgAoN24Y0qQ5uKSswnsQCsa4EQkA9Ac3DDlD3Vq1LEXciAQAyo3QbjLIwcXiLwBQboR2k0EOLhZ/AYByI7SbDHJwsfgLAJQbhWgtNFePv+PV79Dff+vvB66aHADQG0UL0Qrd5WvQNK5axmVQAICy4PB4g1YLiwxyNTkAoFyYaefazaibA3vWIFSTAwDKhZl2rt2MeqjNP9EgVJMDAMqF0M61mzmf1umBrSYHAJQLoZ1rN3Oeveypl5dBcdMOAEArnNPObb1oa8s1x2cv7+pVpTjV6gCAdphp58qysAjV6gCAdphpN+jljLqdQV77HACwMGbaJTPIa58DABZGaJdMWdY+pxgOAMqn0ofHm9cQT2HN8Nnx9XPcFMMBQDlV9oYhzcEjZTNW7lrV2cbHN+qZF56Z137h8gu18407+zAiAKi2ojcMqezhcaqwzxzFcABQTpUNbYLnzFEMBwDlVNnQJnjOXFmK4QAAc1U2tAmeM1eWhWYAAHNVtnq8DFXYKSvDQjMAgLkqG9oSwQMAqJbKHh5HhkVSAKA6Kj3THnQskgIA1cJMu8K4Vh0AqoXQrjCuVQeAaiG0K4xr1QGgWgjtCuNadQColkKhbXu97adt77P94Rbbh23/ab79C7bHGrb9t7z9advrujd0dMIiKQBQLR2rx20vk/RxST8p6YikR2xPR8STDd3eJ+kbEfH/2N4s6Tcl/SfbV0jaLOkNki6S9De2/9+IeLHbO4LWuFYdAKqjyEz7Kkn7ImJ/RHxP0t2Srmnqc42kP84/v0fSWtvO2++OiJMR8S+S9uVfDwAALFKR0L5Y0uGG50fytpZ9IuKUpG9KOr/gawEAQAFFQtst2qJgnyKvzb6AvcV23XZ9ZmamwLAAABgsRUL7iKTVDc8vkXS0XR/b50j6PknHCr5WkhQRkxFRi4jayMhIsdEDADBAioT2I5Ius32p7XOVFZZNN/WZlnR9/vl7JH02IiJv35xXl18q6TJJ/7s7QwcAYLB0rB6PiFO2t0p6QNIySTsi4gnbt0iqR8S0pD+S9Bnb+5TNsDfnr33C9p9JelLSKUkfpHIcAIAz42xCXC61Wi3q9Xq/hwEAQE/Y3hsRtU79WBENAIBEENoAACSC0AYAIBGENgAAiSC0AQBIRCmrx23PSDp4ll9mlaTnujCcMmBfyol9KZ+q7IfEvpTVUu3LmojouLJYKUO7G2zXi5TPp4B9KSf2pXyqsh8S+1JW/d4XDo8DAJAIQhsAgERUObQn+z2ALmJfyol9KZ+q7IfEvpRVX/elsue0AQComirPtAEAqBRCGwCARCQZ2rbX237a9j7bH26xfdj2n+bbv2B7rGHbf8vbn7a9rpfjblZgPz5k+0nbj9neY3tNw7YXbT+afzTf37znCuzLe23PNIz5/Q3brrf9lfzj+ubX9lqBffnthv34su1/bdhWtvdlh+1nbX+pzXbb/t18Xx+z/daGbaV5Xwrsx0Q+/sdsf972mxu2HbD9eP6e9P32gQX25Wrb32z4OfpIw7YFfzZ7rcC+3NCwH1/Kfz9em28r2/uy2vbnbD9l+wnbv9qiT/9/XyIiqQ9l9/T+qqTXSzpX0j9JuqKpz3+RdGf++WZJf5p/fkXef1jSpfnXWVbi/XiXpBX55x+Y3Y/8+Xf6/V4scl/eK+mOFq99raT9+ePK/POVZd6Xpv6/rOwe86V7X/Lx/Jikt0r6UpvtPyVplyRL+mFJXyjp+9JpP350dnySNszuR/78gKRV/X4vFrEvV0v6/1u0L+pnswz70tT33ZI+W+L35XWS3pp//ipJX27x/1jff19SnGlfJWlfROyPiO9JulvSNU19rpH0x/nn90haa9t5+90RcTIi/kXSvvzr9UPH/YiIz0XE8fzpw5Iu6fEYiyrynrSzTtKDEXEsIr4h6UFJ65donEUsdl+ulXRXT0Z2BiLiIUnHFuhyjaRPR+ZhSa+x/TqV7H3ptB8R8fl8nFK5f1eKvCftnM3v2ZJY5L6U/XflaxHxxfzzb0t6StLFTd36/vuSYmhfLOlww/Mjmv8P+1KfiDgl6ZuSzi/42l5Z7Fjep+wvvFnn2a7bftj2Ty/FABeh6L78bH5I6R7bqxf52l4pPJ78dMWlkj7b0Fym96WIdvtbtvdlMZp/V0LSbtt7bW/p05gW60ds/5PtXbbfkLcl+57YXqEsxP6iobm074uzU6pXSvpC06a+/76csxRfdIm5RVvzdWvt+hR5ba8UHovtn5NUk/TOhubRiDhq+/WSPmv78Yj46hKMs4gi+/LXku6KiJO2f0nZkZAfL/jaXlrMeDZLuiciXmxoK9P7UkQKvyuF2X6XstB+R0Pz2/P61KhgAAAgAElEQVT35PslPWj7n/MZYll9Udk61N+x/VOS/krSZUr0Pcm9W9I/RETjrLyU74vtVyr74+LXIuJbzZtbvKSnvy8pzrSPSFrd8PwSSUfb9bF9jqTvU3YIp8hre6XQWGz/hKRtkjZFxMnZ9og4mj/ul/S3yv4q7JeO+xIRzzeM/w8k/VDR1/bYYsazWU2H+0r2vhTRbn/L9r50ZPtNkv5Q0jUR8fxse8N78qyke9W/U2KFRMS3IuI7+ef3SVpue5USfE8aLPS7Upr3xfZyZYE9FRF/2aJL/39f+nni/0w+lB0d2K/ssORsMcYbmvp8UHML0f4s//wNmluItl/9K0Qrsh9XKis8uaypfaWk4fzzVZK+oj4WpBTcl9c1fP4zkh7OP3+tpH/J92ll/vlry7wveb8fVFZI47K+Lw3jGlP7oqeNmltY87/L+L4U2I9RZTUqP9rU/gpJr2r4/POS1pf8Pblw9udKWZAdyt+fQj+bZdqXfPvspOkVZX5f8n/jT0v6nQX69P33JbnD4xFxyvZWSQ8oq6bcERFP2L5FUj0ipiX9kaTP2N6n7Idlc/7aJ2z/maQnJZ2S9MGYe2izbPtxu6RXSvrzrI5OhyJik6TLJX3S9mllR0s+FhFP9mM/pML78iu2Nyn7dz+mrJpcEXHM9kclPZJ/uVti7iG0niq4L1JWVHN35L+xuVK9L5Jk+y5l1cirbB+RdLOk5ZIUEXdKuk9ZRew+Sccl/UK+rVTvS4H9+IiyupXfz39XTkV2J6YLJN2bt50j6U8i4v6e70CDAvvyHkkfsH1K0nclbc5/zlr+bPZhF15SYF+k7I/03RHxbw0vLd37Iuntkq6T9LjtR/O2/67sD8LS/L6wjCkAAIlI8Zw2AAADidAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASMQ5/R5AK6tWrYqxsbF+DwMAgJ7Yu3fvcxEx0qlfKUN7bGxM9Xq938MAAKAnbB8s0o/D4wAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiOoa27dW2P2f7KdtP2P7VFn1s+3dt77P9mO23Nmy73vZX8o/ru70DAICSmpqSxsakoaHscWqq3yNKXpFbc56S9F8j4ou2XyVpr+0HI+LJhj4bJF2Wf7xN0ickvc32ayXdLKkmKfLXTkfEN7q6FwCAcpmakrZskY4fz54fPJg9l6SJif6NK3EdZ9oR8bWI+GL++bclPSXp4qZu10j6dGQelvQa26+TtE7SgxFxLA/qByWt7+oeAADKZ9u2lwN71vHjWTvO2KLOadsek3SlpC80bbpY0uGG50fytnbtrb72Ftt12/WZmZnFDAsAUDaHDi2uHYUUDm3br5T0F5J+LSK+1by5xUtigfb5jRGTEVGLiNrIyEjRYQEAymh0dHHtKKRQaNteriywpyLiL1t0OSJpdcPzSyQdXaAdAFBl27dLK1bMbVuxImvHGStSPW5JfyTpqYj4rTbdpiX9fF5F/sOSvhkRX5P0gKRx2yttr5Q0nrcBKDsqf3E2JiakyUlpzRrJls4/X/p3/0667jp+ns5CkZn22yVdJ+nHbT+af/yU7V+y/Ut5n/sk7Ze0T9IfSPovkhQRxyR9VNIj+ccteRuAMput/D14UIp4ufKX/2ixGBMT0oED0mc+I333u9Lzz/PzdJYc0fIUc1/VarWo1+v9HgYwuMbGsv9Ym61Zk/0nDCwGP08d2d4bEbVO/VgRDcB8VP6im/h56hpCG8B8VP6im4r+PFFH0RGhDWA+Kn/RTUV+nqijKITQBjBfc+XvmjXZc5afxJko8vPECmqFUIgGAOi/oaFsht3Mlk6f7v14eoxCNABAOqijKITQBgD0X9E6igEvViO0AQD9V+S8N8VqnNMGACSiwou0cE4bAFAtLNJCaAMAEkGxGqENAEgEi/4Q2gCARLDoj87p9wAAAChsYmKgQroZM20AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AS2FqShobk4aGssepqX6PCBVwTqcOtndI+o+Sno2If99i+w2SZu9Ifo6kyyWNRMQx2wckfVvSi5JORUStWwMHgNKampK2bJGOH8+eHzyYPZekiYn2rwM6cEQs3MH+MUnfkfTpVqHd1Pfdkn49In48f35AUi0inlvMoGq1WtTr9cW8BADKY2wsC+pma9ZIBw70ejRIgO29RSa2HQ+PR8RDko4V/L7XSrqrYF8AqKZDhxbXDhTUtXPatldIWi/pLxqaQ9Ju23ttb+nW9wKAUhsdXVw7UFA3C9HeLekfIqJxVv72iHirpA2SPpgfam/J9hbbddv1mZmZLg4LAHps+3ZpxYq5bStWZO3AWehmaG9W06HxiDiaPz4r6V5JV7V7cURMRkQtImojIyNdHBYA9NjEhDQ5mZ3DtrPHyUmK0HDWOlaPF2H7+yS9U9LPNbS9QtJQRHw7/3xc0i3d+H4AUHoTE4Q0uq7IJV93Sbpa0irbRyTdLGm5JEXEnXm3n5G0OyL+reGlF0i61/bs9/mTiLi/e0MHAGCwdAztiLi2QJ9PSfpUU9t+SW8+04EBAHDWpqakbduyyv3R0ayuIOEjIF05PA4AQOlUcJEbljEFAFTTtm0vB/as48ez9kQR2gCAaqrgIjeENgCgmiq4yA2hDQCopgouckNoAwCqqYKL3FA9DgCorootcsNMGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gBQQXv27NG6deu0Z8+efg8FXURoA0DF7N69Wxs3bpzziGogtAGgzKampLExaWgoe5yaWrD77t27tWnTJp08eVJXSDp58qQ2bdpEcFcEoQ0AZTU1JW3ZIh08KEVkj1u2tA3uxsD+gKTHJH1ABHeVENoAUFbbtknHj89tO348a2/SHNh3SFqWPxLc1UFoA0BZHTpUqL1VYM/+5z4kgrtKCG0A5bTIc7mVNDrasX3Pnj1tA3tWq+CmqjxNhDaA8lnkudzK2r5dWrFibtuKFVl77rbbbnup6Oz31P4/9aF8+2xx2u23374kQ8bSIrQBlM8izuVW2sSENDkprVkj2dnj5GTWnrvxxhs1PDysJyX9sqTTbb7U6Xz7k5KGh4d1ww03LPnw0X2ENoDyKXgudyBMTEgHDkinT2ePDYEtSWvXrtX09LSGh4f1CUlbNT+4T+ftn1AW2NPT01q7dm0PBp+QRE7HENoAyqfAuVy8bHx8vG1wtwrs8fHxfg21XGaD2pauu6796ZgSBTqhDaB8CpzLxVytgvtFEdhtNdZNSFlYN5o9HVOy+gpH80BLoFarRb1e7/cwAPTT1FT2n+ahQ9kMe/v2eYeGMV/zimiz57AJ7CZjYy8Hdjt29rPXqt+aNdnpii6xvTciah37EdoAUC2NwU1gtzE0NH923WzNmuyPxlb97KzOoEuKhjaHxwGgYsbHx7Vz506tW7dOO3fuJLBb6VQfMXs6pmT1FR1D2/YO28/a/lKb7Vfb/qbtR/OPjzRsW2/7adv7bH+4mwMHALS3du1a3X///VSJt9OqbsLOHhsvrStZfUWRmfanJK3v0OfvIuIt+cctkmR7maSPS9qg7Hr+a21fcTaDBQCgK1pdA/+Zz2SHwhsvrStwrXwvndOpQ0Q8ZHvsDL72VZL2RcR+SbJ9t6RrlNVFAADQXxMTxcK3aL8e6NY57R+x/U+2d9l+Q952saTDDX2O5G0t2d5iu267PjMz06VhAcBZKNH1uYDUndD+oqQ1EfFmZUvb/lXe7hZ925bqRcRkRNQiojYyMtKFYQHAWSjZ9bmA1IXQjohvRcR38s/vk7Tc9iplM+vVDV0vkXT0bL8fAPQE65+jhM46tG1faGcld7avyr/m85IekXSZ7Uttnytps6Tps/1+ANATrH+OEupYiGb7LklXS1pl+4ikmyUtl6SIuFPSeyR9wPYpSd+VtDmyFVtO2d4q6QFJyyTtiIgnlmQvAKDb2q2Exfrn6KMi1ePXdth+h7L7q7fadp+k+85saADQR9u3Z+ewGw+Rs/45+owV0QCglZJdnwtIBWbaADCwSnR9LiAx0wYAdML16qXBTBsA0N7s9eqz5/Znr1eXOArRB8y0AQDtcb16qRDaAID2uF69VAhtAEB7Jbuf9KAjtAGgTMpW9FWy+0kPOkIbAMqi2zcp6cYfAFyvXirOVhwtl1qtFvV6vd/DAIDeGhtrvXTqmjXSgQOL+1rNVd9SNkMmcEvJ9t6IqHXqx0wbAMqim0VfVH1XEqENAGXRzaIvqr4ridAGgLLoZtEXVd+VRGgDQFl0s+iLqu9KYhlTACiTbt2kZPZrbNuWHRIfHc0CmyK0pBHaAFBV3KWscjg8DgBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AJRJ2W7NiVLhOm0AKIvmO3PN3ppT4nprSGKmDQDlwZ250AGhDQBlwZ250AGhDQBlwZ250AGhDQBlwZ250AGhDQBl0c1bc6KSqB4HgDLhzlxYADNtAAAS0TG0be+w/aztL7XZPmH7sfzj87bf3LDtgO3HbT9qu97NgQMAMGiKzLQ/JWn9Atv/RdI7I+JNkj4qabJp+7si4i0RUTuzIQIAAKnAOe2IeMj22ALbP9/w9GFJl5z9sAAAQLNun9N+n6RdDc9D0m7be21vWeiFtrfYrtuuz8zMdHlYALqK9bGBvuha9bjtdykL7Xc0NL89Io7a/n5JD9r+54h4qNXrI2JS+aH1Wq0W3RoXgC5jfWygb7oy07b9Jkl/KOmaiHh+tj0ijuaPz0q6V9JV3fh+APqI9bGBvjnr0LY9KukvJV0XEV9uaH+F7VfNfi5pXFLLCnQACWF9bKBvOh4et32XpKslrbJ9RNLNkpZLUkTcKekjks6X9Pu2JelUXil+gaR787ZzJP1JRNy/BPsAoJdGR7ND4q3aASypItXj13bY/n5J72/Rvl/Sm+e/AkDStm+fe05bYn1soEdYEQ3A4rA+NtA3rD0OYPFYHxvoC2baAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AKRsakoaG5OGhrLHqal+jwhLiLXHASBVU1Nz77h28GD2XGJt+Ipipg0Aqdq2be4tUqXs+bZt/RkPlhyhDQCpOnRoce1IHqENAKkaHV1cO5JHaANAqrZvl1asmNu2YkXWjkoitAEgVRMT0uSktGaNZGePk5MUoVUY1eMAkLKJCUJ6gDDTBgAgEYQ2AACJILQBAEgEoQ2gXFiWE2iLQjQA5cGynMCCmGkDKA+W5QQWRGgDKA+W5QQWRGgDKA+W5QQWRGgDKA+W5QQWRGgDKA+W5QQWRPU4gHJhWU6grUIzbds7bD9r+0ttttv279reZ/sx229t2Ha97a/kH9d3a+AAAAyaoofHPyVp/QLbN0i6LP/YIukTkmT7tZJulvQ2SVdJutn2yjMdLAAAg6xQaEfEQ5KOLdDlGkmfjszDkl5j+3WS1kl6MCKORcQ3JD2ohcMfAAC00a1CtIslHW54fiRva9c+j+0ttuu26zMzM10aFgAA1dGt0HaLtligfX5jxGRE1CKiNjIy0qVhAQBQHd0K7SOSVjc8v0TS0QXaAQDAInUrtKcl/XxeRf7Dkr4ZEV+T9ICkcdsr8wK08bwNAAAsUqHrtG3fJelqSatsH1FWEb5ckiLiTkn3SfopSfskHZf0C/m2Y7Y/KumR/EvdEhELFbQBAIA2CoV2RFzbYXtI+mCbbTsk7Vj80AAAQCOWMQUAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkIhCoW17ve2nbe+z/eEW23/b9qP5x5dt/2vDthcbtk13c/AAAAySczp1sL1M0scl/aSkI5IesT0dEU/O9omIX2/o/8uSrmz4Et+NiLd0b8gAAAymIjPtqyTti4j9EfE9SXdLumaB/tdKuqsbgwMAAC8rEtoXSzrc8PxI3jaP7TWSLpX02Ybm82zXbT9s+6fbfRPbW/J+9ZmZmQLDAgBgsBQJbbdoizZ9N0u6JyJebGgbjYiapP8s6Xds/0CrF0bEZETUIqI2MjJSYFgAAAyWIqF9RNLqhueXSDrapu9mNR0aj4ij+eN+SX+ruee7AQBAQUVC+xFJl9m+1Pa5yoJ5XhW47R+UtFLSPza0rbQ9nH++StLbJT3Z/FoAANBZx+rxiDhle6ukByQtk7QjIp6wfYukekTMBvi1ku6OiMZD55dL+qTt08r+QPhYY9U5AAAoznMzthxqtVrU6/V+DwMAgJ6wvTev/1oQK6IBAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIIbQAAEkFoAwCQCEIbAIBEENoAACSC0AYAIBGENgAAiSC0AQBIBKENAEAiCG0AABJBaAMAkAhCGwCARBDaAAAkgtAGACARhDYAAIkgtAEASAShDQBAIghtAAASQWgDAJAIQhsAgEQQ2gAAJILQBgAgEYQ2AACJILQBAEgEoQ0AQCIKhbbt9baftr3P9odbbH+v7Rnbj+Yf72/Ydr3tr+Qf13dz8AAADJJzOnWwvUzSxyX9pKQjkh6xPR0RTzZ1/dOI2Nr02tdKullSTVJI2pu/9htdGX2F7Hp+l+44eoe+/sLXdcHyC7T1oq3acP6GM+4HAKieIjPtqyTti4j9EfE9SXdLuqbg118n6cGIOJYH9YOS1p/ZUJfWrud3aePjG1X7Yk0bH9+oXc/v6un3vvXQrXrmhWcUCj3zwjO69dCt88ZQtB8AoJqKhPbFkg43PD+StzX7WduP2b7H9upFvla2t9iu267PzMwUGFb39DsM7zh6h07EiTltJ+KE7jh6xxn1AwBUU5HQdou2aHr+15LGIuJNkv5G0h8v4rVZY8RkRNQiojYyMlJgWN3T7zD8+gtfL9RetB8AoJqKhPYRSasbnl8i6Whjh4h4PiJO5k//QNIPFX1tGfQ7DC9YfkGh9qL9AADVVCS0H5F0me1LbZ8rabOk6cYOtl/X8HSTpKfyzx+QNG57pe2VksbztlLpdxhuvWirzvN5c9rO83naetHWM+oHAKimjtXjEXHK9lZlYbtM0o6IeML2LZLqETEt6Vdsb5J0StIxSe/NX3vM9keVBb8k3RIRx5ZgP87K1ou26tZDt845RN7LMJyt/u5UFV60HwCgmhzR8hRzX9VqtajX6z39nlxKBQDoF9t7I6LWqV/Hmfag2HD+BkIaAFBqLGM6APp5DToAoHuYaVfc7DXos+frZ69Bl8SRBQBIDDPtiuv3NegAgO4htCuu39egAwC6h9CuuH5fgw4A6B5Cu+JYkAUAqoNCtB7qx7XgLMgCANVBaPdIP6u4uQYdAKphYA6P9/taZaq4AQBnayBCu9/3y5aqW8Xd7z+GAGCQDERol2GWW8Uq7jL8MQQAg2QgQrsMs9wqVnGX4Y8hABgkAxHa/Z7lzlaNn4gTGsr/yS9cfqFuGr0p6QKxMvwxBACDZCBCu5+z3MZDyJJ0Wqdf+t4pBPZC56z7/ccQAAyagQjtDedv0E2jN+nC5RfKck9nuSkfQu50zrqKh/wBoMwG5jrtfl2rnPIh5IX+4Gj892ThFgDojYEJ7X65YPkFLx0ab24vuyJ/cLBwCwD0zkAcHu+nlA8hc84aAMqF0O6y5sItSX07n362Uv6DAwCqiMPjXdRuffGbRm/Szjfu7PPoFo9z1gBQLsy0u+hsK8X37NmjdevWac+ePUsxvDOy4fwN2vnGnaq/ta6db9xJYANAHxHaXXQ2leK7d+/Wxo0b5zwCANCI0O6iMy3c2r17tzZt2qSTJ0/qvNefp5MnT2rTpk0ENwBgDkK7i86kcKsxsFe9Z5Uuv/tyrXrPqo7Bzd21AGDwENpdtNiV15oDe/VvrJaXWat/Y/WCwc3dtQBgMDki+j2GeWq1WtTr9X4PY0m1DOwhv7Q9TocO/+ZhPXfPcxoeHtb09LTGx8clSRsf39hywZYLl1+YZJU6AAw623sjotapHzPtPtizZ8+CgS1JHpo/456tKk95aVQAwJkjtHts1/O7tPkjm18qOlt94/zAnuUha/WNq18qTrv99tslsVIZAAwqQruHZs9Fv+a618jnWif2n9Dh2w4rTrc+RRGnQ4dvO6wT+09oeHhYN9xwgyRWKgOAQUVo99Ds4iuvvurV+oHf+gH5XOu5e57T4d+cH9ytzmmvXbtWUn9vNQoA6J9Cy5jaXi/pf0paJukPI+JjTds/JOn9kk5JmpH0ixFxMN/2oqTH866HImJTl8aenMZzzq/+kSy4v/qhr+q5e56TpJfObS9UhDaLu2sBwODpONO2vUzSxyVtkHSFpGttX9HU7f9IqkXEmyTdI+m2hm3fjYi35B8DG9jS/HPOs8E9dO7QyzPuFzsHNgBgMBU5PH6VpH0RsT8ivifpbknXNHaIiM9FxPH86cOSLunuMKuh1bno7//R79ctU7doeHhYz93znJ7a/BSBDQBoqcjh8YslHW54fkTS2xbo/z5Jjat8nGe7ruzQ+cci4q9avcj2FklbJGl0dLTAsNLT9q5ZV27Qf3j1f9CmTZteKjojsAEAzYqEdqvrkVqWO9v+OUk1Se9saB6NiKO2Xy/ps7Yfj4ivzvuCEZOSJqVscZUC40pSu3PR4+Pj2rlzp26//XbdcMMNLxWdAQAwq0hoH5G0uuH5JZKONney/ROStkl6Z0ScnG2PiKP5437bfyvpSknzQhvS2rVrCWsAQFtFzmk/Iuky25faPlfSZknTjR1sXynpk5I2RcSzDe0rbQ/nn6+S9HZJT3Zr8Og+bkQCAOXVcaYdEadsb5X0gLJLvnZExBO2b5FUj4hpSbdLeqWkP7ctvXxp1+WSPmn7tLI/ED4WEaUP7V3P75p/3nkALq+aXfzlRJyQpJduRCJpIPYfAMqOG4Y0aQ4uKVttbBAWL+FGJADQH9ww5AzNrlrW6ESc0B1H7+jTiHqHG5EAQLkR2k0GObi4EQkAlBuh3WSQg4sbkQBAuRHaTQY5uLgRCQCUG4VoLTRXj7/j1e/Q33/r7weumhwA0BtFC9EK3eVr0DSuWsZlUACAsuDweINWC4sMcjU5AKBcmGnn2s2omwN71iBUkwMAyoWZdq7djHqozT/RIFSTAwDKhdDOtZs5n9bpga0mBwCUC6Gdazdznr3sqZeXQXHTDgBAK5zTzm29aGvLNcdnL+/qVaU41eoAgHaYaefKsrAI1eoAgHaYaTfo5Yy6nUFe+xwAsDBm2iUzyGufAwAWRmiXTFnWPqcYDgDKp9KHx5vXEE9hzfDZ8fVz3BTDAUA5VfaGIc3BI2UzVu5a1dnGxzfqmReemdd+4fILtfONO/swIgCotqI3DKns4XGqsM8cxXAAUE6VDW2C58xRDAcA5VTZ0CZ4zlxZiuEAAHNVNrQJnjNXloVmAABzVbZ6vAxV2Ckrw0IzAIC5KhvaEsEDAKiWyh4eBwCgagjtimNlMwCojkofHh90rGwGANXCTLvCWGAGAKqF0K4wFpgBgGohtCuMBWYAoFoI7QpjgRkAqJZCoW17ve2nbe+z/297dx8jV1XGcfz7S1tqeAm2tGptCy2xGttEQ90QBaMgmlYMXY0krlFTDIaIxZeYmIAmmNR/NP6hMWoQpQkQbYtVTDUgVIsxsWnpQqBvWFi2CJuaUCkWK6Zk6+Mf9yy5HWe7Z0p35tzt75Ns9sw5507Os889PXPvnM7o5jbtMyVtTO07JC2qtd2S6vdLWnH6hm4T8SebmZlNLRPuHpc0DfgR8CFgBNgpaXNE7Kt1ux54MSLeImkA+A7wCUlLgQFgGfBm4A+S3hoRx093INaeP2DGzGzqyLnSvhQYiojhiHgF2AD0t/TpB+5M5U3AVZKU6jdExLGIOAAMpeczMzOzDuUs2vOB52qPR1Jd2z4RMQocAS7IPBYASTdIGpQ0eOjQobzRm5mZnUFyFm21qYvMPjnHVpURt0dEX0T0zZ07N2NYZmZmZ5acRXsEWFh7vAA4OF4fSdOB84HDmceamZlZhpxFeyewRNJiSWdRbSzb3NJnM7A6la8FtkZEpPqBtLt8MbAEePj0DN3MzOzMMuHu8YgYlXQT8AAwDVgXEXslrQUGI2IzcAdwt6QhqivsgXTsXkn3APuAUWCNd46bmZmdGlUXxGXp6+uLwcHBXg/DzMysKyQ9EhF9E/XzJ6KZmZk1hBdtMzOzhvCibWZm1hBetM3MzBqiyI1okg4Bf3uNTzMH+MdpGE4JHEuZHEt5pkoc4FhKNVmxXBQRE36yWJGL9ukgaTBnJ1yQNdYAAAWiSURBVF4TOJYyOZbyTJU4wLGUqtex+Pa4mZlZQ3jRNjMza4ipvGjf3usBnEaOpUyOpTxTJQ5wLKXqaSxT9j1tMzOzqWYqX2mbmZlNKV60zczMGqKRi7aklZL2SxqSdHOb9pmSNqb2HZIW1dpuSfX7Ja3o5rhbZcTxVUn7JO2S9EdJF9Xajkt6LP20flVq12XEcp2kQ7Uxf67WtlrSU+lndeux3ZYRy/dqcTwp6Z+1ttLysk7S85L2jNMuST9Ise6StLzWVkxeMuL4VBr/LknbJL2z1vaMpN0pJz3/JqKMWK6QdKR2Ht1aazvpudltGbF8rRbHnjQ/Zqe20vKyUNJDkp6QtFfSl9v06f18iYhG/VB9PejTwMXAWcDjwNKWPl8AbkvlAWBjKi9N/WcCi9PzTCs4jiuBs1P5xrE40uOjvc5Fh7FcB/ywzbGzgeH0e1Yqzyo5lpb+X6T6utri8pLG8z5gObBnnPargfsBAe8GdhSal4niuGxsfMCHx+JIj58B5vQ6Fx3EcgXwuzb1HZ2bJcTS0vcaYGvBeZkHLE/l84An2/w71vP50sQr7UuBoYgYjohXgA1Af0uffuDOVN4EXCVJqX5DRByLiAPAUHq+Xpgwjoh4KCJeTg+3Awu6PMZcOTkZzwpgS0QcjogXgS3AykkaZ45OY/kksL4rIzsFEfFnqu+4H08/cFdUtgOvlzSPwvIyURwRsS2NE8qeKzk5Gc9rmWeTosNYSp8rf4+IR1P5X8ATwPyWbj2fL01ctOcDz9Uej/D/f9hX+0TEKHAEuCDz2G7pdCzXU73CG/M6SYOStkv66GQMsAO5sXw83VLaJGlhh8d2S/Z40tsVi4GtteqS8pJjvHhLy0snWudKAA9KekTSDT0aU6feI+lxSfdLWpbqGpsTSWdTLWK/qlUXmxdVb6leAuxoaer5fJk+GU86ydSmrvX/rY3XJ+fYbskei6RPA33A+2vVF0bEQUkXA1sl7Y6IpydhnDlyYvktsD4ijkn6PNWdkA9kHttNnYxnANgUEcdrdSXlJUcT5ko2SVdSLdrvrVVfnnLyBmCLpL+mK8RSPUr1OdRHJV0N/AZYQkNzklwD/CUi6lflReZF0rlULy6+EhEvtTa3OaSr86WJV9ojwMLa4wXAwfH6SJoOnE91Cyfn2G7JGoukDwLfAFZFxLGx+og4mH4PA3+ielXYKxPGEhEv1Mb/U+Bducd2WSfjGaDldl9heckxXryl5WVCkt4B/Azoj4gXxuprOXkeuJfevSWWJSJeioijqXwfMEPSHBqYk5qTzZVi8iJpBtWC/fOI+HWbLr2fL7184/9UfqjuDgxT3ZYc24yxrKXPGk7ciHZPKi/jxI1ow/RuI1pOHJdQbTxZ0lI/C5iZynOAp+jhhpTMWObVyh8DtqfybOBAimlWKs8uOZbU721UG2lUal5q41rE+JuePsKJG2seLjEvGXFcSLVH5bKW+nOA82rlbcDKwnPyprHzimohezblJ+vcLCmW1D520XROyXlJf+O7gO+fpE/P50vjbo9HxKikm4AHqHZTrouIvZLWAoMRsRm4A7hb0hDVyTKQjt0r6R5gHzAKrIkTb22WFsd3gXOBX1b76Hg2IlYBbwd+Ium/VHdLvh0R+3oRB2TH8iVJq6j+7oepdpMTEYclfQvYmZ5ubZx4C62rMmOBalPNhkgzNikqLwCS1lPtRp4jaQT4JjADICJuA+6j2hE7BLwMfDa1FZWXjDhupdq38uM0V0aj+iamNwL3prrpwC8i4vddD6AmI5ZrgRsljQL/AQbSedb23OxBCK/KiAWqF+kPRsS/a4cWlxfgcuAzwG5Jj6W6r1O9ICxmvvhjTM3MzBqiie9pm5mZnZG8aJuZmTWEF20zM7OG8KJtZmbWEF60zczMGsKLtpmZWUN40TYzM2uI/wGB3OgTrUJRRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a71a09ecc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for index_data, data in enumerate(dataset):\n",
    "    \n",
    "    # an array of all centroid distances per data point (to find min)\n",
    "    dist_vec = np.zeros((k,1))\n",
    "    \n",
    "    for index_centroids, centroid in enumerate(centroids):\n",
    "        \n",
    "        # populate distance vectors of all centroids for that data point\n",
    "        dist_vec[index_centroids] = eucl_distance(centroid, data)\n",
    "    \n",
    "    # find which centroid is closest to that data point (which cluster that datapoint belongs to)\n",
    "    belongs_to[index_data] = np.argmin(dist_vec) # argmin returns indices, not values\n",
    "        \n",
    "# plot clusters\n",
    "colors = ['red', 'limegreen']\n",
    "fig = plt.figure(figsize=(8,14))\n",
    "ax = fig.add_subplot(211)\n",
    "\n",
    "for i, b in enumerate(belongs_to):\n",
    "    ax.scatter(dataset[i,0], dataset[i,1], color=colors[int(b)])\n",
    "for i, c in enumerate(centroids):\n",
    "    ax.scatter(centroids[i,0], centroids[i,1], color=colors[i], marker=\"D\", s=80, edgecolor='black', linewidth='2')\n",
    "for index in range(len(centroids)):\n",
    "    instances = [i for i in range(len(belongs_to)) if belongs_to[i] == index]\n",
    "    centroids[index] = np.mean(dataset[instances], axis=0)\n",
    "\n",
    "# plot updated centroids\n",
    "ax2 = fig.add_subplot(212)\n",
    "\n",
    "for i, b in enumerate(belongs_to):\n",
    "    ax2.scatter(dataset[i,0], dataset[i,1], color=colors[int(b)])\n",
    "for i, c in enumerate(centroids):\n",
    "    ax2.scatter(centroids[i,0], centroids[i,1], color=colors[i], marker=\"D\", s=80, edgecolor='black', linewidth='2')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Centroid stabilization - END\n",
    "\n",
    "At this point the centroids are positioned at the mean of their respective clusters and thus the K-Means search is terminated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Datapoint 0 belongs to cluster 1\n",
      "...\n",
      "Datapoint 5 belongs to cluster 1\n",
      "...\n",
      "Datapoint 10 belongs to cluster 1\n",
      "...\n",
      "Datapoint 15 belongs to cluster 1\n",
      "...\n",
      "Datapoint 20 belongs to cluster 0\n",
      "...\n",
      "Datapoint 25 belongs to cluster 0\n",
      "...\n",
      "Datapoint 30 belongs to cluster 0\n",
      "...\n",
      "Datapoint 35 belongs to cluster 0\n",
      "...\n"
     ]
    }
   ],
   "source": [
    "# datapoints per cluster\n",
    "for i in range(len(belongs_to)):\n",
    "    if i%5 == 0:\n",
    "        print('Datapoint {} belongs to cluster {}\\n...'.format(i, int(belongs_to[i])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Elbow Method\n",
    "\n",
    "What happens if we do not know k beforehand? The solution to this is to use a method called \"The Elbow Method\" and goes as follows. We start with a k=1 and we compute the sum of all eucledean distances of each datapoint from its corresponding centroid. We repeat for k=2 etc. Finally, we compare the sum-of-errors of all the different values of k and select the k that precedes the phase of very slow error decrease."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAGrCAYAAAARoORAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYHHd95/HPt3t6zhrd3RNZli17esCQA8MqDlknNsiYBUI49kmysAkhWTZms5A1CZuD7BHI9YQNYHaTQNZggrMBA+EILHGyOPggzhKzsjHEjkxGI8u2ZFkzkmYkzdnT3d/9o6rHo9EcPX1VH+/X8/TT1VXV1d+eeaTP/H71q1+ZuwsAAMQrEXcBAACAQAYAoCkQyAAANAECGQCAJkAgAwDQBAhkAACaAIGMtmFml5nZtJklK3z/UTN7WbT862b20dpWuObnvsTMjtXoWD9jZvfX4lj1/Fwze76ZHazk2GbmZpatpM4yP+/7zOz/1uv4wFoIZLScKDjnovAtPS5x9yfdPXD3QrWf4e6/6+7/thb1rlTvQKmHOtT8W5Let+z4q/1O/7CGn1c2d/+2pCkz+9E4Ph+di0BGq/rRKHxLj6fjLgjlMbPdkl4q6S9WbFr5O317DOWVfELSW2P8fHQgAhltw8z2RS25ruj1vWb2W2b2d2Z23sy+Yma7lu3/JjN7wsxOm9l/WnGsd5vZn6047pvN7EkzO7V8fzPrM7PbzWzSzA6Z2a+s1QVtZl+LFr8VtQL/1bJt7zSzcTM7YWY/u2x9j5m9L/rsk2b2x2bWt/6Pwv7AzM6a2WNmdsOyDVvN7LboM46b2W+XuvjNLGtm90XvO2Vmn96o5nUK+H0zu9/Mtq6y+UZJD7n7/EbHWcerzOxIVOfvm1ki+tyEmf3n6Pc6bmZ/Wqoh+h29M1reE/1O//2y737GzCw6/r2SbjCznipqBDaFQEa7+9eSflZSRlK3pP8ohecwJX1Y0pskXSJpp6RLNzjWD0l6rqQbJP1XM3tetP43JO2TdKXCsPmptQ7g7tdFiy+IWoGfjl5/l6StkvZIeoukPzKz7dG290p6jqSrJWWjff7rOnX+gKQjknZFtX3ezHZE226XlI+O80JJL5dU6pr/LUlfkbRd4c/iDzao+SJRIH5E0vdJerm7n11lt++V9J116i/H6yXtl/QiSa+V9G+i9T8TPV6q8PcRSCp1fd8n6SXR8vUKf0bXR6+vk/S3Hs0l7O7HJS0q/H0DDUEgo1X9hZlNRY+VXZ/L/Ym7/5O7z0n6jMJQk6Qfk/Rld/+auy9I+i+Siht85nvcfc7dvyXpW5JeEK3/CUm/6+6T7n5M0v+o4PssSvpNd1909zslTUt6btRi+zlJv+juZ9z9vKTflfSGdY41LumD0bE+rTD8fsTMhiS9UtI73H3G3ccl3bLsWIuSLpd0ibvPu/tmB4elJN0haYfC7ufZNfbbJun8KuuX/06nzOzn1vms90Y/jyclfVDSG6P1PynpA+5+xN2nJb1L0huiXpP7JP1w1Jq+TtJ/k3Rt9L7ro+3LnY9qBRqiK+4CgAq9zt3/poz9nlm2PKuwxSSFreKnShvcfcbMTtfiWCuWy3Xa3fOrHD8tqV/Sg8/2psokrTeS/HippRd5IqrxcoWheWLZsRLL6v0Vha3kb5jZpKT3u/vHNvEdsgr/SLnG3XPr7DcpaXCV9eX+TqULf8al76fo+YkV27okDbn7mJlNK/yj7IcVfte3mNlzFQbyyj+kBiVNlVkPUDVayOhUJyTtLb0ws36F3daVHmt5d/fetXaswClJc5K+2923RY+t7h6s8549y86FStJlkp5WGGILknYtO9YWd/9uSXL3Z9z959z9EoUDmj60yZHVhxSeHvirKOTW8m2FXfDVWP4zLn0/Rc+Xr9iWl3Qyen2fwt6R7qhb+j5JP62wm/7h0pvM7BKFpziq7VoHykYgo1N9VtKrzeyHzKxb0m+q8n8Pn5H0LjPbbmZ7JG00OvikwvObG3L3oqSPSLrFzDLS0oCkf7HO2zKS/oOZpczsxyU9T9Kd7n5C4Tni95vZluh877CZXR8d98fNrPSHxaQkl1S6hKysmt39Dkm/LulvzGx4jd3ukvQiM+vd6Hjr+OXo571X0s2SSue175D0i2Z2hZkFCrv3P72s9+E+hb+f0kC1eyX9gqT7V1wu9xJJd0enM4CGIJDRkdz9UUlvk/RJhS3cSUmVTs7xm9F7H5f0NwrDfr3/yN8t6fboPOlPlHH8X5V0WNLfm9m56DPWa4E+IGlEYev6dyT9mLuXuuN/WmHL7x8VfufPStodbft+SQ9E3bpfknSzuz++2Zrd/XaFP5O7zWzfKttPSrpb4WCs5f63XXgd8hfW+ZgvSnpQYav2LyXdFq3/mKT/pTBwH5c0rzBwS+5T2BVdCuT7FZ4S+Jou9JOS/nidzwdqzi481QSgWmb285Le4O7Xb7hzh4pGud+u8HxzU/0nZGbfK+lWd//BuGtBZyGQgSpZONHFlZK+rrBl+peS/tDdPxhrYQBaCqOsgep1S/qfkq5QOCr3U5I+FGtFAFoOLWQAAJoAg7oAAGgCDe2y3rVrl+/bt6+RHwkAQGwefPDBU+6eLmffhgbyvn37dPBgWbdABQCg5ZnZExvvFaLLGgCAJlB2IJtZ0sy+aWZfjl5fYWYPmNmomX06mu0IAABUYDMt5JsVzlVb8l5Jt7j7iMIZf95Sy8IAAOgkZQVyNL/tj0j6aPTaJB1QOO2eFM6487p6FAgAQCcot4X8QYW3ZivdL3anpKllE7YfU3jTdAAAUIENA9nMXi1p3N0fXL56lV1XnWHEzG4ys4NmdnBiYqLCMgEAaG/ltJCvlfQaMzuqcErAAwpbzNvMrHTZ1KV69n6kF3D3W919v7vvT6fLuhQLAICOs2Egu/u73P1Sd98n6Q0K7xH6k5LuUXijb0l6s8LboQEAgApUcx3yr0r6JTM7rPCc8m0b7A8AANawqZm63P1eSfdGy0ckXVP7kgAA6DzM1AUAQBMgkAEAaAIEMgAATYBABgCgCbRsIJ+bX9Qjx8/GXQYAADXRsoH8oXvG9PoP/Z3yheLGOwMA0ORaNpCzmUCLBdcTZ2bjLgUAgKq1dCBL0uHx6ZgrAQCgei0byMPpAUkEMgCgPbRsIA/2prR7a6/GCGQAQBto2UCWwm7rUQIZANAGWjqQh9OBxiamVSyueitmAABaRksHcjYTaDZX0Ilz83GXAgBAVVo6kEcYaQ0AaBMtHcilS59GT56PuRIAAKrT0oG8M+jR9v6UxiZoIQMAWltLB7IkjWQG6bIGALS8lg/k4ejSJ3dGWgMAWlfLB3I2E2hqdlGnZ3JxlwIAQMXaIpAlRloDAFpbywcylz4BANpBywfy7q29GuhOEsgAgJbW8oFsZhrOBAQyAKCltXwgS+F5ZAIZANDK2iaQnzk3r3Pzi3GXAgBARdojkNPhwC7ujQwAaFXtEciMtAYAtLi2COTLdvSrO5nQYea0BgC0qLYI5K5kQlfsGtDhkwQyAKA1tUUgS9FIa1rIAIAW1TaBPJwJ9NSZWc0vFuIuBQCATWubQB7JBCq69PipmbhLAQBg09omkEsjrUcZaQ0AaEFtE8hX7BpQwrj0CQDQmtomkHtTSV22o5/JQQAALaltAlkKu61Hx8/HXQYAAJvWVoE8nAn0+KkZ5QvFuEsBAGBTNgxkM+s1s2+Y2bfM7FEze0+0/uNm9riZPRw9rq5/uevLpgMtFlxPnpmNuxQAADalq4x9FiQdcPdpM0tJut/M/ira9svu/tn6lbc5I0ODksKBXVdGN5wAAKAVbNhC9lBppFQqenhdq6rQcHpAEpc+AQBaT1nnkM0saWYPSxqXdJe7PxBt+h0z+7aZ3WJmPWu89yYzO2hmBycmJmpU9uoGe1P6ri29jLQGALScsgLZ3QvufrWkSyVdY2bfI+ldkq6S9P2Sdkj61TXee6u773f3/el0ukZlr21kiDmtAQCtZ1OjrN19StK9kl7h7iei7uwFSX8i6Zo61Ldpw+lAh8enVSw2Za86AACrKmeUddrMtkXLfZJeJukxM9sdrTNJr5P0SD0LLVc2E2g2V9CJc/NxlwIAQNnKGWW9W9LtZpZUGOCfcfcvm9ndZpaWZJIelvTv6lhn2UpzWh8en9aebX0xVwMAQHk2DGR3/7akF66y/kBdKqrSyLJAvv459T9nDQBALbTVTF2StDPo0fb+lA4zhSYAoIW0XSBLYbc1d30CALQSAhkAgCbQpoE8qMnZRZ2eXoi7FAAAytKmgRwO7GIKTQBAq2jrQKbbGgDQKtoykC/Z2quB7iSBDABoGW0ZyGamYQZ2AQBaSFsGsiRl0wQyAKB1tG0gD2cCPXNuXufnF+MuBQCADbVtIJem0BybmIm5EgAANta2gbx06dNJptAEADS/tg3ky3b0qzuZ0OEJziMDAJpf2wZyVzKhK3YNaIyBXQCAFtC2gSyF3dbM1gUAaAVtHcjDmUBPnZnV/GIh7lIAAFhXWwdyNhOo6NLjpxhpDQBobm0dyCPMaQ0AaBFtHchX7BpQwrjrEwCg+bV1IPemktq7o5+R1gCAptfWgSwxpzUAoDW0fyAPBXr81IzyhWLcpQAAsKb2D+R0oFyhqCfPzMZdCgAAa2r/QGakNQCgBXROIDOnNQCgibV9IA/2pvRdW3p1+CSBDABoXm0fyFLYSqaFDABoZh0TyGPj03L3uEsBAGBVHRPIM7mCTpydj7sUAABW1TGBLDGFJgCgeXVUIHPpEwCgWXVEIO8c6Nb2/hSBDABoWh0RyGYWjrQePx93KQAArKojAlmKLn2ihQwAaFIdE8jD6UCTs4s6Pb0QdykAAFykYwJ5ZGhQEgO7AADNacNANrNeM/uGmX3LzB41s/dE668wswfMbNTMPm1m3fUvt3Jc+gQAaGbltJAXJB1w9xdIulrSK8zsxZLeK+kWdx+RNCnpLfUrs3qXbO1Vf3eSFjIAoCltGMgeKqVYKnq4pAOSPhutv13S6+pSYY2YmYbTgcaY0xoA0ITKOodsZkkze1jSuKS7JI1JmnL3fLTLMUl71njvTWZ20MwOTkxM1KLmio0w0hoA0KTKCmR3L7j71ZIulXSNpOetttsa773V3fe7+/50Ol15pTUwnAl04uy8zs8vxloHAAArbWqUtbtPSbpX0oslbTOzrmjTpZKerm1ptVca2DU2MRNzJQAAXKicUdZpM9sWLfdJepmkQ5LukfRj0W5vlvTFehVZKyPMaQ0AaFJdG++i3ZJuN7OkwgD/jLt/2cz+UdKnzOy3JX1T0m11rLMmLtvRr+5kQqNMoQkAaDIbBrK7f1vSC1dZf0Th+eSW0ZVMaN+ufo3RQgYANJmOmamrhDmtAQDNqAMDeVBPnpnV/GIh7lIAAFjSgYEcqOjS46cYaQ0AaB6dF8hpRloDAJpPxwXylekBmRHIAIDm0nGB3JtK6rId/QQyAKCpdFwgS2G3NYEMAGgmnRnImUCPn5pRvlCMuxQAACR1cCDnCkU9NTkXdykAAEjq4ECWpNGTTKEJAGgOHRnIw6WbTExwHhkA0Bw6MpC39KY0tKWHgV0AgKbRkYEsSSOZQW4yAQBoGh0byKWbTLh73KUAANC5gTycCTSTK+jE2fm4SwEAoHMDeSTDnNYAgObRsYG8dOkTgQwAaAIdG8g7B7q1rT9FCxkA0BQ6NpDNTNl0wEhrAEBT6NhAlqSRoYDJQQAATaGjA3k4HejMTE6npxfiLgUA0OE6OpCzjLQGADQJAlnMaQ0AiF9HB/IlW/vU353U6EkCGQAQr44O5ETCNJwONEYLGQAQs44OZOnZOa0BAIgTgZwJdOLsvM7PL8ZdCgCggxHI0cCusYmZmCsBAHQyAplLnwAATaDjA/nyHf1KJY1ABgDEquMDuSuZ0BW7BghkAECsOj6QpdJI6/NxlwEA6GAEsqRsOtCTZ2Y1v1iIuxQAQIcikCVlhwZVdOnoaUZaAwDiQSArbCFLYgpNAEBsCGRJV6YHZMalTwCA+GwYyGa218zuMbNDZvaomd0crX+3mR03s4ejx6vqX2599KaS2ru9n7s+AQBi01XGPnlJ73T3h8xsUNKDZnZXtO0Wd39f/cprnJFMoDFayACAmGzYQnb3E+7+ULR8XtIhSXvqXVijZTOBjkzMKF8oxl0KAKADbeocspntk/RCSQ9Eq95uZt82s4+Z2fY13nOTmR00s4MTExNVFVtPw5lAuUJRT03OxV0KAKADlR3IZhZI+pykd7j7OUkfljQs6WpJJyS9f7X3ufut7r7f3fen0+kalFwfzGkNAIhTWYFsZimFYfwJd/+8JLn7SXcvuHtR0kckXVO/MuuvFMijzNgFAIhBOaOsTdJtkg65+weWrd+9bLfXS3qk9uU1zpbelIa29NBCBgDEopxR1tdKepOkfzCzh6N1vy7pjWZ2tSSXdFTSW+tSYQNlGWkNAIjJhoHs7vdLslU23Vn7cuI1khnUnx98Su6usGMAAIDGYKauZYYzgWZyBZ04Ox93KQCADkMgL1Oa05rzyACARiOQl+HSJwBAXAjkZXYF3drWn2JOawBAwxHIy5iZsulAh7kNIwCgwQjkFbKZgBYyAKDhCOQVsplAZ2ZyOj29EHcpAIAOQiCvwMAuAEAcCOQVlgKZbmsAQAMRyCtcsrVPfakkLWQAQEMRyCskEhYO7CKQAQANRCCvgkAGADQagbyKbCbQibPzml7Ix10KAKBDEMirGI7mtOZWjACARiGQVzEyFAbyKIEMAGgQAnkVl+/oVyppnEcGADQMgbyKrmRC+3YOEMgAgIYhkNcwMhRojMlBAAANQiCvIZsO9MTpGc0vFuIuBQDQAQjkNQxnAhVdOnp6Ju5SAAAdgEBeAzeZAAA0EoG8huF0IDMCGQDQGATyGnpTSe3d3s+1yACAhiCQ15HNBMzWBQBoCAJ5HSOZQEdOzShfKMZdCgCgzRHI6xjOBMrli3pqci7uUgAAbY5AXgcjrQEAjUIgr4NABgA0CoG8ji29KQ1t6SGQAQB1RyBvIJsJdHj8fNxlAADaHIG8gWw60NjEjNw97lIAAG2MQN5ANhNoeiGvZ87Nx10KAKCNEcgbyGYGJUmjJzmPDACoHwJ5A4y0BgA0AoG8gV1Bt7b2pXR4gkAGANQPgbwBM9NIJqCFDACoqw0D2cz2mtk9ZnbIzB41s5uj9TvM7C4zG42et9e/3HhkCWQAQJ2V00LOS3qnuz9P0oslvc3Mni/p1yR91d1HJH01et2WsplAZ2ZyOjOTi7sUAECb2jCQ3f2Euz8ULZ+XdEjSHkmvlXR7tNvtkl5XryLjNszALgBAnW3qHLKZ7ZP0QkkPSBpy9xNSGNqSMmu85yYzO2hmBycmJqqrNiYjBDIAoM7KDmQzCyR9TtI73P1cue9z91vdfb+770+n05XUGLtLtvapL5XUKFNoAgDqpKxANrOUwjD+hLt/Plp90sx2R9t3SxqvT4nxSyRMw5kBWsgAgLopZ5S1SbpN0iF3/8CyTV+S9OZo+c2Svlj78ppHNh1ojEAGANRJOS3kayW9SdIBM3s4erxK0u9JutHMRiXdGL1uWyNDg3r67LymF/JxlwIAaENdG+3g7vdLsjU231DbcprXcDoc2DU2Pq0X7N0WczUAgHbDTF1lYk5rAEA9Echlunxnv1JJY05rAEBdEMhlSiUT2rdzgNswAgDqgkDehGwm0BgtZABAHRDIm5DNBHri9IwW8oW4SwEAtBkCeROymUBFlx4/NRN3KQCANkMgbwIjrQEA9UIgb8JwOpAZgQwAqD0CeRN6U0nt3d5PIAMAao5A3qRsJiCQAQA1RyBvUjYT6MipGRWKHncpAIA2QiBvUjYdKJcv6qkzs3GXAgBoIwTyJmWHGGkNAKg9AnmTSpc+jRLIAIAaIpA3aUtvSpnBHlrIAICaIpArkM0E3PUJAFBTBHIFRjKBxsan5c5IawBAbRDIFchmAk0v5PXMufm4SwEAtAkCuQLDzGkNAKgxArkCI5lBSQQyAKB2COQK7Aq6tbUvxaVPAICaIZArYGbMaQ0AqCkCuULZdDjSGgCAWiCQKzQyFOj0TE5nZnJxlwIAaAMEcoUYaQ0AqCUCuULZNIEMAKgdArlCe7b1qS+VJJABADVBIFcokTANZwY0On4+7lIAAG2AQK4CI60BALVCIFchmwn09Nl5zSzk4y4FANDiCOQqZKMpNMe4FSMAoEoEchWy0aVPoycJZABAdQjkKly+s19dCdNhWsgAgCoRyFVIJRPat2uAS58AAFUjkKs0wk0mAAA1QCBXKZsJ9MTpGS3kC3GXAgBoYQRylbKZQEWXjp6ajbsUAEAL2zCQzexjZjZuZo8sW/duMztuZg9Hj1fVt8zmleUmEwCAGiinhfxxSa9YZf0t7n519LiztmW1juF0IDMxhSYAoCobBrK7f03SmQbU0pJ6U0ldur2PFjIAoCrVnEN+u5l9O+rS3r7WTmZ2k5kdNLODExMTVXxc88qmGWkNAKhOpYH8YUnDkq6WdELS+9fa0d1vdff97r4/nU5X+HHNbWRoUEdOzahQ9LhLAQC0qIoC2d1PunvB3YuSPiLpmtqW1Vqy6UC5fFFPnWGkNQCgMhUFspntXvby9ZIeWWvfTjDMSGsAQJW6NtrBzO6Q9BJJu8zsmKTfkPQSM7takks6Kumtdayx6S1d+jQxrZdpKOZqAACtaMNAdvc3rrL6tjrU0rK29qWUGezhrk8AgIoxU1eNZDMBd30CAFSMQK6RbCbQ2Pi03BlpDQDYPAK5RkYygaYX8jp5biHuUgAALYhArpHSSGum0AQAVIJArhFuMgEAqAaBXCPpoEdbersIZABARQjkGjEzjQwNapRABgBUgECuoWw6HGkNAMBmEcg1lM0EOj2T0+RMLu5SAAAthkCuoezQs1NoAgCwGQRyDWXT0aVPTKEJANgkArmG9mzrU18qyUhrAMCmEcg1lEiYrkwP0GUNANg0ArnGRjKBDp9kti4AwOYQyDWWzQR6+uy8ZhbycZcCAGghBHKNlabQHKPbGgCwCQRyjWUzg5KY0xoAsDkEco1dvrNfXQljCk0AwKYQyDWWSia0b9cALWQAwKYQyHXAnNYAgM0ikOtgZCjQE2dmtZAvxF0KAKBFEMh1kM0EKhRdR0/Nxl0KAKBFEMh1MBzNac15ZABAuQjkOhhOBzIjkAEA5SOQ66CvO6lLt/dpdJwpNAEA5SGQ6ySbDmghAwDKRiDXSTYT6MipGRWKHncpAIAWQCDXyUhmULl8UccmGWkNANgYgVwnw9FNJkZP0m0NANgYgVwnpbs+HeauTwCAMhDIdbK1L6X0YA8DuwAAZSGQ62gkE3DXJwBAWQjkOspmwptMuDPSGgCwPgK5jrKZQNMLeZ08txB3KQCAJkcg11GWOa0BAGUikOsoOxRd+sQUmgCADWwYyGb2MTMbN7NHlq3bYWZ3mdlo9Ly9vmW2pnTQoy29XbSQAQAbKqeF/HFJr1ix7tckfdXdRyR9NXqNFcxM2QxzWgMANrZhILv71ySdWbH6tZJuj5Zvl/S6GtfVNkYygxpjchAAwAYqPYc85O4nJCl6zqy1o5ndZGYHzezgxMREhR/XurKZQKemc5qcycVdCgCgidV9UJe73+ru+919fzqdrvfHNR2m0AQAlKPSQD5pZrslKXoer11J7WUpkDmPDABYR6WB/CVJb46W3yzpi7Upp/3s2dan3lSCuz4BANZVzmVPd0j6uqTnmtkxM3uLpN+TdKOZjUq6MXqNVSQSpuF0QJc1AGBdXRvt4O5vXGPTDTWupW1lM4EOHp2MuwwAQBNjpq4GGMkEOj41p5mFfNylAACaFIHcAKWBXVyPDABYC4HcAIy0BgBshEBugMt3DqgrYQQyAGBNBHIDpJIJ7ds1oFECGQCwBgK5QbLpQGMEMgBgDQRyg2QzgZ44M6tcvhh3KQCAJkQgN0g2E6hQdB09PRN3KQCAJkQgN0hppDVTaAIAVkMgN8hwOpAZlz4BAFZHIDdIX3dSe7b1Mac1AGBVBHIDjWQCWsgAgFURyA2UzQQam5hWoehxlwIAaDIEcgNlM4Fy+aKOTc7GXQoAoMkQyA3EnNYAgLUQyA2UTQ9KElNoAgAuQiA30Nb+lNKDPbSQAQAXIZAbLJtmpDUA4GIEcoONDIU3mXBnpDUA4FkEcoNlM4HOL+R18txC3KUAAJoIgdxg2TQjrQEAFyOQG+zZS5/Ox1wJAKCZEMgNlh7s0ZbeLi59AgBcgEBuMDNTljmtAQArEMgxKM1pDQBACYEcg2wm0KnpnCZncnGXAgBoEgRyDEYy4RSa3BsZAFBCIMeAm0wAAFYikGOwZ1ufelMJAhkAsIRAjkEiYRpmTmsAwDIEcky49AkAsByBHJNsOtDxqTnNLOTjLgUA0AQI5JiUBnYdmZiJuRIAQDMgkGMyMhQG8ihzWgMARCDH5vKdA+pKGOeRAQCSCOTYpJIJXb6zn0AGAEiSuqp5s5kdlXReUkFS3t3316KoTjGSGdQ/0WUNAFBtWsgvdferCePNy2YCPXF6Vrl8Me5SAAAxo8s6RtlMoELRdfQ0I60BoNNVG8gu6Stm9qCZ3bTaDmZ2k5kdNLODExMTVX5ce2FOawBASbWBfK27v0jSKyW9zcyuW7mDu9/q7vvdfX86na7y49rLcDqQmTR6kkAGgE5XVSC7+9PR87ikL0i6phZFdYq+7qT2bOvjNowAgMoD2cwGzGywtCzp5ZIeqVVhnYI5rQEAUnWXPQ1J+oKZlY7zSXf/65pU1UGy6UBfHzutQtGVTFjc5QAAYlJxILv7EUkvqGEtHWlkKNBCvqhjk7O6fOdA3OUAAGLCZU8xY6Q1AEAikGOXTQ9KIpABoNMRyDHb2p9SerBH935nQqMnz8vd4y4JABCDquayRm3c+PwhffKBJ3XjLV/Tpdv7dOCqjF56VUY/eOVO9aaScZcHAGgAa2SLbP/+/X7w4MGGfV4reXpqTvd8Z1z3PDau+w+f0vxiUX2ppK7N7tRLr8rowFUZ7d7aF3eZAIBNMLMHy73XA4HchOYXC/r6kdO657FxffXQuI5PzUmSnrd7iw549QF4AAAMYklEQVRcldaBq4Z09d5tXCYFAE2OQG4j7q7R8Wnd/di47n5sXA8+MalC0bVjoFvXPyetl16V0fUjaW3tT8VdKgBgBQK5jZ2dXdR9oxO657Fx3fudcU3OLiqZMP2zy7frwFUZ3XBVRtlMoGjCFgBAjAjkDlEouh5+ajJqPU/o0IlzkrQ0MOzAVRm9mIFhABAbArlDMTAMAJoLgYwLBobd/di4jk0yMAwAGo1AxgU2Ghh24KqMrntOWlv7GBgGALVEIGNdDAwDgMYgkFE2BoYBQP0QyKjYegPDDlw1pJdelWZgGACUiUBGTTAwDACqQyCj5hgYBgCbRyCj7tYbGHblrgFtH+jWjv5ubetPacdAt7b1d2vHQLe296e0pTelBK1qAB1gM4HM7RdRka39Kb3mBZfoNS+4JBoYNqW7Hzupvx09pa8+Nq6p2ZwWC6v/sZcwafvKsO7v1raBlHb0d2t7f3cY6AOppW1b+lJ0jQNoa7SQURfurumFvKZmF3VmJqczszlNzeZ0ZmYxes4tbZucjR4zi8oViqsez0za1pdaCuvt/WFr+9nWd2pZKzzctq2/mxAHECtayIidmWmwN6XB3pT27ugv6z3urtlc4dmwns1pshTYMzlNRuumZnM6PjWnR46f1ZnZnHL5tUN8S28ptFNRF3oY3s+G+oXBvq0/pVQyUcsfBQCUhUBG0zAzDfR0aaCnS3t3lPced9fcYkGTs4tL4X1m5tkAD1vf4bZnzs3r0IlzOjOb0/zi6iEuSVt6u7Q9Cujt/aml7vXS8/L1W/vCcB/oTjKRCoCqEMhoaWam/u4u9Xd3ac+28q+PnssVLugqLy2v7Eo/PZ3T4fFpTc0uanohv+bxUklbCuptfcsCfCB8vX15kA90a1tf+Lq7i9Y4gBCBjI7U151UX3efLtlEiOfyRZ2dC8+Bl1rfU7NhgE/OltaHr584PauHn5rS1Oza58UlaaA7udRVvl5LfPlrRqkD7YlABsrU3ZVQerBH6cGest9TOi8+NRd2m0+VgnxuUVNRt/pSkM8t6vjUnCZnczo7t6i1xlsmTNoatbAvCPBSS3zg4i71rX0p9aeSBDnQxAhkoI6WnxffTJd6oeg6N7cYBnnUEp+cCV+XArwU5ifPzes7z5zX5GxOs7nCOrVIA91dCnq6NNCTVNDTpaC3K1zXG64PolqXLw/2rlyX1EB3F+EO1BiBDDShZMLCkeAD3bpCA2W/byFf0NmoC30pyGcXdXZuUTMLeU0v5JeepxcKmp5f1Knzueh1uC1fLO9SyIHuZBjoPV0ajMJ7+fLFIZ9U0JPSQE/ygpAn3IEQgQy0kZ6upDJbksps6a3o/e6uhXxxKZzPzy8P8LxmFgqaXliMwjzalssvLZ+Zmb0g9NeaHGal/u7kUniXWu3Pts7DIA+WWvUpDfaG27YsLYfPXLKGVkYgA1hiZupNJdWbSmpXUP658rUs5EvBXbigFX4+el499MPn41Nzml5YDN87n193cFxJbyqxFM6DvSltiYJ7sOfC4A56u6JtF64f7O1STxe3GkU8CGQAddPTlVRPkNTOoPpj5aKW+/R8XufmF3V+Pq/zK58XwuVz8/ml9SfOzi9tX+8ce0l3V2JFWF8c6Ku1zpcvc/9wVIJABtASursS2tEVTo9aqXwhDPXzF4T6imBfyF+0fuL89NK69a5HX6o1mbgopAdXtMi3RMHdm0qqL5VUbyqx9Hr58tK2LkbJtzsCGUDH6EomosvFKg/1QtGjUF890M+tFvLzeR09Nbu0PJ3Lr3lZ23q6uxLq7bo4rHtWhnpXUn3dSfVEQd7XnbzgfctDv2+VPwJ6Uwl1JxPMPtdgBDIAbEIyYdral6rq3t/Foms6l9f8YkHzuaLm8wXNLxY0lytoPl8M1y89ippbtrzWtqm5Rc2fLSwdq7RtrbneN2KmZ8O66+LWe2lbT/S6O5lQT1dCqWRC3V3RI7niedlyab+eZetTy7b3RMud1CtAIANAgyUSpi294axr9VYshiPn55YF+VwU2AuLYYDP5aKgz4d/FCxEfxSEfyA8G+4Ly/4oODu3uBT8C/nwPbl8UblCsaLW/1qSCVs10FeuK4X5hQFv6k4ml/YL/2AoHS+57P0WPSe1tS+l7710a+2+wCYQyADQxhIJi6aKbcxAM3dXoejKFaKAzhe1kC9qsVC8YF0pvC94jvZbWGVdab+FVdbl8kXNzuaVK7hy+cKK9/rSfuV4/u4tuvPmH67zT2l1BDIAoGbMTF1JU1cyoSpO1dec+7N/JCyFdL6oXKEQ/cEQrovzhi9VBbKZvULSf5eUlPRRd/+9mlQFAEANmVl4GV4TX2de8Z8CZpaU9EeSXinp+ZLeaGbPr1VhAAB0kmra5tdIOuzuR9w9J+lTkl5bm7IAAOgs1QTyHklPLXt9LFp3ATO7ycwOmtnBiYmJKj4OAID2VU0gr3Zx2EWD3d39Vnff7+770+l0FR8HAED7qiaQj0nau+z1pZKerq4cAAA6UzWB/P8kjZjZFWbWLekNkr5Um7IAAOgsFV/25O55M3u7pP+j8LKnj7n7ozWrDACADlLVdcjufqekO2tUCwAAHSu+KUkAAMASAhkAgCZAIAMA0AQIZAAAmgCBDABAEyCQAQBoAgQyAABNwNwvmn66fh9mNiHpiRoecpekUzU8Xpz4Ls2nXb6HxHdpVu3yXdrle0i1/y6Xu3tZN3JoaCDXmpkddPf9cddRC3yX5tMu30PiuzSrdvku7fI9pHi/C13WAAA0AQIZAIAm0OqBfGvcBdQQ36X5tMv3kPguzapdvku7fA8pxu/S0ueQAQBoF63eQgYAoC0QyAAANIGWDGQz+5iZjZvZI3HXUi0z22tm95jZITN71MxujrumSphZr5l9w8y+FX2P98RdU7XMLGlm3zSzL8ddSzXM7KiZ/YOZPWxmB+Oup1Jmts3MPmtmj0X/Xn4w7poqYWbPjX4Xpcc5M3tH3HVVysx+Mfo3/4iZ3WFmvXHXVAkzuzn6Do/G9ftoyXPIZnadpGlJf+ru3xN3PdUws92Sdrv7Q2Y2KOlBSa9z93+MubRNMTOTNODu02aWknS/pJvd/e9jLq1iZvZLkvZL2uLur467nkqZ2VFJ+929pSduMLPbJf2tu3/UzLol9bv7VNx1VcPMkpKOS/oBd6/lpEkNYWZ7FP5bf767z5nZZyTd6e4fj7eyzTGz75H0KUnXSMpJ+mtJP+/uo42soyVbyO7+NUln4q6jFtz9hLs/FC2fl3RI0p54q9o8D01HL1PRo/X+2ouY2aWSfkTSR+OuBZKZbZF0naTbJMndc60expEbJI21Yhgv0yWpz8y6JPVLejrmeirxPEl/7+6z7p6XdJ+k1ze6iJYM5HZlZvskvVDSA/FWUpmoi/dhSeOS7nL3lvwekQ9K+hVJxbgLqQGX9BUze9DMboq7mApdKWlC0p9EpxE+amYDcRdVA2+QdEfcRVTK3Y9Lep+kJyWdkHTW3b8Sb1UVeUTSdWa208z6Jb1K0t5GF0EgNwkzCyR9TtI73P1c3PVUwt0L7n61pEslXRN1A7UcM3u1pHF3fzDuWmrkWnd/kaRXSnpbdMqn1XRJepGkD7v7CyXNSPq1eEuqTtTt/hpJfx53LZUys+2SXivpCkmXSBows5+Kt6rNc/dDkt4r6S6F3dXfkpRvdB0EchOIzrl+TtIn3P3zcddTragr8V5Jr4i5lEpdK+k10bnXT0k6YGZ/Fm9JlXP3p6PncUlfUHierNUck3RsWa/LZxUGdCt7paSH3P1k3IVU4WWSHnf3CXdflPR5Sf885poq4u63ufuL3P06hadEG3r+WCKQYxcNhrpN0iF3/0Dc9VTKzNJmti1a7lP4D/WxeKuqjLu/y90vdfd9CrsU73b3lvurX5LMbCAaLKioi/flCrvnWoq7PyPpKTN7brTqBkktNfBxFW9UC3dXR56U9GIz64/+L7tB4TiYlmNmmej5Mkn/UjH8broa/YG1YGZ3SHqJpF1mdkzSb7j7bfFWVbFrJb1J0j9E518l6dfd/c4Ya6rEbkm3R6NGE5I+4+4tfblQmxiS9IXw/0p1Sfqku/91vCVV7BckfSLq6j0i6Wdjrqdi0XnKGyW9Ne5aquHuD5jZZyU9pLCL95tq3Wk0P2dmOyUtSnqbu082uoCWvOwJAIB2Q5c1AABNgEAGAKAJEMgAADQBAhkAgCZAIAMA0AQIZAAAmgCBDABAE/j/usCLQvTT4GkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a71dc378d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "errors = []\n",
    "\n",
    "for i in range(1, 10):\n",
    "    # keep starting state the same for all 10 examples\n",
    "    clf = KMeans(n_clusters=i, init=\"k-means++\", n_init=10, random_state=0)\n",
    "    clf.fit(dataset)\n",
    "\n",
    "    # Sum of squared distances of samples to their closest cluster center.\n",
    "    errors.append(clf.inertia_)\n",
    "\n",
    "plt.figure(figsize=(8,7))\n",
    "plt.plot(range(1, 10), errors)\n",
    "plt.title(\"Finding the best k (Elbow)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two things become obvious:\n",
    "\n",
    "1. The best k is at value 2, as the error decreases very slowly after that\n",
    "2. Why this method is called \"elbow\" method"
   ]
  }
 ],
 "metadata": {
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}