homework4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "def mahalanobis_distance(x, y):\n",
    "    #X = np.vstack([x,y])\n",
    "    #V = np.cov(X.T)\n",
    "    #return np.sqrt(np.sum(np.dot(x - y, np.linalg.inv(V)) * (x - y), axis = 1))\n",
    "    return cdist(x, y, metric = \"mahalanobis\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuzzifier = 1.01\n",
    "\n",
    "def calculate_eta(assignation, distance):\n",
    "    u = assignation ** fuzzifier\n",
    "    d_squared = distance ** 2\n",
    "    return np.sum(np.dot(u, d_squared.T)) / np.sum(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_u(x, v, eta):\n",
    "    d = mahalanobis_distance(x.T, v).T\n",
    "    u = 1.0 / (1.0 + (d / eta) ** (2.0 / fuzzifier - 1.0))\n",
    "    return u, d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cluster(x, groups, centers, precision, max_iterations):\n",
    "    features, points = x.shape\n",
    "    eta = 1.0\n",
    "    v = np.zeros((max_iterations, groups, features))\n",
    "    v[0] = np.array(centers)\n",
    "    u = np.random.random((max_iterations, groups, points))\n",
    "    \n",
    "    i = 0\n",
    "    while i < max_iterations - 1:\n",
    "        u[i], distance = calculate_u(x, v[i], eta)\n",
    "        v[i + 1] = calculate_new_centers(x, u[i])\n",
    "        eta = calculate_eta(u[i], distance)\n",
    "        if np.linalg.norm(v[i + 1] - v[i]) < precision:\n",
    "            break\n",
    "        i += 1\n",
    "    plot(x, v[i], u[i - 1], groups)\n",
    "\n",
    "\n",
    "def calculate_new_centers(x, u):\n",
    "    u_temp = u ** fuzzifier\n",
    "    return u_temp.dot(x.T) / np.atleast_2d(u_temp.sum(axis=1)).T\n",
    "    \n",
    "    \n",
    "def plot(data, centers, assignation, groups):\n",
    "    fig, ax = pyplot.subplots()\n",
    "    group_num = np.argmax(assignation, axis = 0)\n",
    "    for group in range(groups):\n",
    "        ax.scatter(data[0][group_num == group], data[1][group_num == group])\n",
    "    for center in centers:\n",
    "        ax.scatter(center[0], center[1], c = 'black')      \n",
    "    pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGI1JREFUeJzt3W+QVNd95vHvs4MQiNgCSSODGbRAPFFFchGbdCQs7245ggxIOIYXcgqXazXrpYoqR4kn/2xBXBXNykmVtE4FQ62tXdZohVIuIUIUQzy2CULy+oUkpAFkBMKEMXLECDDjGiAuLGSBf3nRZ0gzt5mhu2emp7ufT1VX9/3dc7vPQaN55t577m1FBGZmZoX+Q7U7YGZm44/DwczMMhwOZmaW4XAwM7MMh4OZmWU4HMzMLMPhYGZmGQ4HMzPLcDiYmVnGhGp3oFw33XRTzJ49u9rdMDOrKXv27PlpRDQP165mw2H27Nl0d3dXuxtmZjVF0r9cTTsfVjIzs4xhw0HS45JOSTowqP6Hkg5LOijpfxbU10jqSesWF9SXpFqPpNUF9TmSdks6IulpSRNHanBmZlaeq9lzeAJYUliQ9NvAMmBeRNwO/HWq3wasAG5P23xNUpOkJuCrwD3AbcCnUluAR4G1EdEKnAZWVjooMzOrzLDhEBHfB/oHlT8LPBIR76Q2p1J9GbA5It6JiDeAHuCO9OiJiKMR8QtgM7BMkoC7ga1p+03A8grHZGZmFSr3nMOvAf85HQ76/5J+K9VnAscK2vWm2pXqNwJnIuLCoLqZmVVRubOVJgDTgAXAbwFbJM0FVKRtUDyEYoj2RUlaBawCuOWWW0rsspmZXa1y9xx6gWci72Xgl8BNqT6roF0LcHyI+k+BqZImDKoXFREbIiIXEbnm5mGn6ZrVrv1bYO0HoXNq/nn/lmr3yBpMueHwTfLnCpD0a8BE8r/otwMrJF0raQ7QCrwMvAK0pplJE8mftN4e+e8ofR64L71vO7Ct3MGY1YX9W+AfPwdnjwGRf/7HzzkgbExdzVTWp4AXgVsl9UpaCTwOzE3TWzcD7Wkv4iCwBXgd+C7wQERcTOcU/gDYARwCtqS2AA8CfyKph/w5iI0jO0SzGrPrYXj37ctr776dr5uNEeX/eK89uVwufIW01aXOqRQ/9SboPDPWvbE6I2lPROSGa+crpM3Gm+tbSqubjQKHg9l4s/Av4JrJl9eumZyvm40Rh4PZeDPv9+B318P1swDln393fb5uNkZq9q6sZnVt3u85DAyArqNdrNu7jpPnTjJ9ynQ65newdO7SUf9ch4OZ2TjVdbSLzhc6OX/xPAAnzp2g84VOgFEPCB9WMjMbp9btXXcpGAacv3iedXvXjfpnOxzMzMapk+dOllQfSQ4HM7NxavqU6SXVR5LDwcysHGNw/6uO+R1Mapp0WW1S0yQ65neM+GcN5hPSZmalGrj/1cBtTgbufwUjOsts4KRzNWYr+fYZZmalWvvBdGPEQa6fBX98IFsfR3z7DDOz0XK2t7R6DXI4mJmVqgHuf+VwMDMrVQPc/8rhYGZWqga4/5VnK5mZlaPO73/lPQczM8twOJiZWcbVfIf045JOpe+LHrzuzySFpJvSsiStl9Qjab+k+QVt2yUdSY/2gvpvSnotbbNekkZqcGZmVp6r2XN4AlgyuChpFvA7wJsF5XuA1vRYBTyW2t4APATcCdwBPCRpWtrmsdR2YLvMZ5mZ2dgaNhwi4vtAf5FVa4EvcPk3oS8Dnoy8l4CpkmYAi4GdEdEfEaeBncCStO69EfFi5C/VfhJYXtmQzMysUmWdc5D0CeCtiPjBoFUzgcJryntTbah6b5G6mZlVUclTWSVdB3wRaCu2ukgtyqhf6bNXkT8ExS233DJsX83MrDzl7Dn8KjAH+IGkHwMtwF5J08n/5T+roG0LcHyYekuRelERsSEichGRa25uLqPrZmZ2NUoOh4h4LSJujojZETGb/C/4+RFxEtgO3J9mLS0AzkbECWAH0CZpWjoR3QbsSOt+JmlBmqV0P7BthMZmZmZlupqprE8BLwK3SuqVtHKI5t8GjgI9wP8Ffh8gIvqBLwGvpMfDqQbwWeDraZsfAd8pbyhmZjZS/H0OZmYNxN/nYGZmZXM4mJlZhsPBzMwyHA5mZpbhcDAzswyHg5mZZTgczMwsw+FgZmYZDgczM8twOJiZWYbDwczMMhwOZmaW4XAwM7MMh4OZmWU4HMzMLMPhYGZmGQ4HMzPLcDiYmVmGw8HMzDKGDQdJj0s6JelAQe3Lkn4oab+kf5A0tWDdGkk9kg5LWlxQX5JqPZJWF9TnSNot6YikpyVNHMkBmplZ6a5mz+EJYMmg2k7ggxExD/hnYA2ApNuAFcDtaZuvSWqS1AR8FbgHuA34VGoL8CiwNiJagdPAyopGZGZmFRs2HCLi+0D/oNo/RcSFtPgS0JJeLwM2R8Q7EfEG0APckR49EXE0In4BbAaWSRJwN7A1bb8JWF7hmMzMrEIjcc7hvwPfSa9nAscK1vWm2pXqNwJnCoJmoG5mZlVUUThI+iJwAfjGQKlIsyijfqXPWyWpW1J3X19fqd01M7OrVHY4SGoHPg58OiIGfqH3ArMKmrUAx4eo/xSYKmnCoHpREbEhInIRkWtubi6362ZmNoyywkHSEuBB4BMR8fOCVduBFZKulTQHaAVeBl4BWtPMpInkT1pvT6HyPHBf2r4d2FbeUMzMbKRczVTWp4AXgVsl9UpaCfwv4D3ATkmvSvrfABFxENgCvA58F3ggIi6mcwp/AOwADgFbUlvIh8yfSOohfw5i44iO0MzMSqZ/PyJUW3K5XHR3d1e7G2ZmNUXSnojIDdfOV0ibmVmGw8HMzDIcDmZmluFwMDOzDIeDmZllOBzMzCzD4WBmZhkOBzMzy3A4mJlZhsPBzMwyHA5mZpbhcDAzswyHg5mZZTgczMwsw+FgZmYZDgczM8twOJiZWYbDwczMMq7mO6Qfl3RK0oGC2g2Sdko6kp6npbokrZfUI2m/pPkF27Sn9kcktRfUf1PSa2mb9ZI00oM0M7PSXM2ewxPAkkG11cCuiGgFdqVlgHuA1vRYBTwG+TABHgLuBO4AHhoIlNRmVcF2gz/LzMzG2LDhEBHfB/oHlZcBm9LrTcDygvqTkfcSMFXSDGAxsDMi+iPiNLATWJLWvTciXoyIAJ4seC8zM6uScs85vC8iTgCk55tTfSZwrKBdb6oNVe8tUjczsyoa6RPSxc4XRBn14m8urZLULam7r6+vzC6amdlwyg2Hn6RDQqTnU6neC8wqaNcCHB+m3lKkXlREbIiIXETkmpuby+y6mZkNp9xw2A4MzDhqB7YV1O9Ps5YWAGfTYacdQJukaelEdBuwI637maQFaZbS/QXvZWZmVTJhuAaSngI+BtwkqZf8rKNHgC2SVgJvAp9Mzb8N3Av0AD8HPgMQEf2SvgS8kto9HBEDJ7k/S35G1GTgO+lhZmZVpPwkodqTy+Wiu7u72t0wM6spkvZERG64dr5C2szMMhwOZmaW4XAwM7MMh4OZmWU4HMzMLMPhYGZmGQ4HMzPLcDiYmVmGw8HMzDIcDmZmluFwMDOzDIeDmZllOBzMzCzD4WBmZhkOBzMzy3A4mJlZhsPBzMwyHA5mZpbhcDAzs4yKwkHSH0s6KOmApKckTZI0R9JuSUckPS1pYmp7bVruSetnF7zPmlQ/LGlxZUMyM7NKlR0OkmYCnwNyEfFBoAlYATwKrI2IVuA0sDJtshI4HREfANamdki6LW13O7AE+JqkpnL7ZWZmlav0sNIEYLKkCcB1wAngbmBrWr8JWJ5eL0vLpPULJSnVN0fEOxHxBtAD3FFhv8zMrAJlh0NEvAX8NfAm+VA4C+wBzkTEhdSsF5iZXs8EjqVtL6T2NxbWi2xzGUmrJHVL6u7r6yu362ZmNoxKDitNI/9X/xzg/cAU4J4iTWNgkyusu1I9W4zYEBG5iMg1NzeX3mkzM7sqlRxWWgS8ERF9EfEu8AxwFzA1HWYCaAGOp9e9wCyAtP56oL+wXmQbMzOrgkrC4U1ggaTr0rmDhcDrwPPAfalNO7Atvd6elknrn4uISPUVaTbTHKAVeLmCfpmZWYUmDN+kuIjYLWkrsBe4AOwDNgBdwGZJf5lqG9MmG4G/ldRDfo9hRXqfg5K2kA+WC8ADEXGx3H6ZmVnllP/jvfbkcrno7u6udjfMzGqKpD0RkRuuna+QNjOzDIeDmZllOBzMzCzD4WBmZhkOBzMzy3A4mJlZhsPBzMwyHA5mZpbhcDAzswyHg5mZZTgczMwsw+FgZmYZDgcbE4sWLULSpceiRYuq3SUzG4LDwUbdokWL2LVr12W1Xbt2OSDMxjGHg426wcEwXN3Mqs/hYGZmGQ4HMzPLcDjYqFu4cGFJdTOrvorCQdJUSVsl/VDSIUkfkXSDpJ2SjqTnaamtJK2X1CNpv6T5Be/TntofkdRe6aBsfHn22WczQbBw4UKeffbZMfn8rqNdtG1tY96mebRtbaPraNeYfK5ZLZtQ4fbrgO9GxH2SJgLXAX8O7IqIRyStBlYDDwL3AK3pcSfwGHCnpBuAh4AcEMAeSdsj4nSFfbNxZKyCYLCuo110vtDJ+YvnAThx7gSdL3QCsHTu0qr0yawWlL3nIOm9wH8BNgJExC8i4gywDNiUmm0ClqfXy4AnI+8lYKqkGcBiYGdE9KdA2AksKbdfZoXW7V13KRgGnL94nnV711WpR2a1oZLDSnOBPuD/Sdon6euSpgDvi4gTAOn55tR+JnCsYPveVLtS3axiJ8+dLKluZnmVhMMEYD7wWER8GDhH/hDSlahILYaoZ99AWiWpW1J3X19fqf21BjR9yvSS6maWV0k49AK9EbE7LW8lHxY/SYeLSM+nCtrPKti+BTg+RD0jIjZERC4ics3NzRV03RpFx/wOJjVNuqw2qWkSHfM7qtQjs9pQdjhExEngmKRbU2kh8DqwHRiYcdQObEuvtwP3p1lLC4Cz6bDTDqBN0rQ0s6kt1cwqtnTuUjrv6mTGlBkIMWPKDDrv6vTJaLNhVDpb6Q+Bb6SZSkeBz5APnC2SVgJvAp9Mbb8N3Av0AD9PbYmIfklfAl5J7R6OiP4K+2V2ydK5Sx0GZiVSRNHD++NeLpeL7u7uanfDzKymSNoTEbnh2vkKaTMzy3A4mJlZhsPBzMwyHA5mZpbhcDAzswyHg5mZZTgczMwsw+FgZmYZDgczM8twOJiZWYbDwczMMiq98V7N+ua+t/jyjsMcP/M27586mc8vvpXlH/Z3DJmZQYOGwzf3vcWaZ17j7XcvAvDWmbdZ88xrAA4IMzMa9LDSl3ccvhQMA95+9yJf3nG4Sj0yMxtfGjIcjp95u6S6mVmjachweP/UySXVzcwaTUOGw+cX38rka5ouq02+ponPL771CluYmTWWhjwhPXDS2bOVzMyKa8hwgHxAOAzMzIqr+LCSpCZJ+yR9Ky3PkbRb0hFJT0uamOrXpuWetH52wXusSfXDkhZX2iczM6vMSJxz6AAOFSw/CqyNiFbgNLAy1VcCpyPiA8Da1A5JtwErgNuBJcDXJF1+QsDMzMZUReEgqQVYCnw9LQu4G9iammwClqfXy9Iyaf3C1H4ZsDki3omIN4Ae4I5K+mVmZpWpdM/hK8AXgF+m5RuBMxFxIS33AgMH9mcCxwDS+rOp/aV6kW0uI2mVpG5J3X19fRV23czMrqTsE9KSPg6ciog9kj42UC7SNIZZN9Q2lxcjNgAbAHK5XNE244Xv3WRmtayS2UofBT4h6V5gEvBe8nsSUyVNSHsHLcDx1L4XmAX0SpoAXA/0F9QHFG5Tk3zvJjOrdWUfVoqINRHREhGzyZ9Qfi4iPg08D9yXmrUD29Lr7WmZtP65iIhUX5FmM80BWoGXy+3XeOB7N5lZrRuN6xweBDZL+ktgH7Ax1TcCfyuph/wewwqAiDgoaQvwOnABeCAiLmbftnb43k1mVutGJBwi4nvA99LroxSZbRQR54FPXmH7vwL+aiT6Mh68f+pk3ioSBL53k5nVioa8t9Jo872bzKzWNeztM0aT791kZrXO4TBKfO8mM6tlPqxkZmYZDgczM8to2MNKvoLZzOzKGjIcfAWzmdnQGvKwkq9gNjMbWkOGg69gNjMbWkOGw5WuVPYVzGZmeQ0ZDr6C2cxsaA15QtpXMJuZDa0hwwF8BbOZ2VAa8rCSmZkNzeFgZmYZDgczM8twOJiZWYbDwczMMsoOB0mzJD0v6ZCkg5I6Uv0GSTslHUnP01JdktZL6pG0X9L8gvdqT+2PSGqvfFj1r+toF21b25i3aR5tW9voOtpV7S6ZWR2pZCrrBeBPI2KvpPcAeyTtBP4bsCsiHpG0GlgNPAjcA7Smx53AY8Cdkm4AHgJyQKT32R4RpyvoW1H1cifWrqNddL7QyfmL5wE4ce4EnS90ArB07tIq9szM6kXZew4RcSIi9qbXPwMOATOBZcCm1GwTsDy9XgY8GXkvAVMlzQAWAzsjoj8Fwk5gSbn9upKBO7G+deZtgn+/E+s397010h816tbtXXcpGAacv3iedXvXValHZlZvRuScg6TZwIeB3cD7IuIE5AMEuDk1mwkcK9isN9WuVB9R9XQn1pPnTpZUNzMrVcXhIOlXgL8H/igi/nWopkVqMUS92GetktQtqbuvr6+kftbTnVinT5leUt3MrFQVhYOka8gHwzci4plU/kk6XER6PpXqvcCsgs1bgOND1DMiYkNE5CIi19zcXFJf6+lOrB3zO5jUNOmy2qSmSXTM76hSj8ys3lQyW0nARuBQRPxNwartwMCMo3ZgW0H9/jRraQFwNh122gG0SZqWZja1pdqIqqc7sS6du5TOuzqZMWUGQsyYMoPOuzp9MtrMRkwls5U+CvxX4DVJr6banwOPAFskrQTeBD6Z1n0buBfoAX4OfAYgIvolfQl4JbV7OCL6K+hXUbV4J9auo12s27uOk+dOMn3KdDrmd1wKgKVzlzoMzGzUKKLo4f1xL5fLRXd3d7W7MWoGT1eF/KEj7yGYWSUk7YmI3HDtfIX0OOXpqmZWTQ6HccrTVc2smhwO45Snq5pZNTkcxilPVzWzamrYrwkd7wZOOl9ptpKZ2WhyOIxjnq5qZtXiw0pmZpbhcDAzswyHg5mZZTgczMwsw+FgZmYZDgczM8twOJiZWYbDwczMMmr2lt2S+oB/KXGzm4CfjkJ3xqNGGit4vPWskcYKoz/e/xgRw36VZs2GQzkkdV/NfczrQSONFTzeetZIY4XxM14fVjIzswyHg5mZZTRaOGyodgfGUCONFTzeetZIY4VxMt6GOudgZmZXp9H2HMzM7Co0RDhIWiLpsKQeSaur3Z+RIOlxSackHSio3SBpp6Qj6XlaqkvS+jT+/ZLmV6/npZM0S9Lzkg5JOiipI9XrdbyTJL0s6QdpvP8j1edI2p3G+7Skial+bVruSetnV7P/5ZLUJGmfpG+l5bodr6QfS3pN0quSulNtXP081304SGoCvgrcA9wGfErSbdXt1Yh4AlgyqLYa2BURrcCutAz5sbemxyrgsTHq40i5APxpRPw6sAB4IP03rNfxvgPcHRG/AXwIWCJpAfAosDaN9zSwMrVfCZyOiA8Aa1O7WtQBHCpYrvfx/nZEfKhg2ur4+nmOiLp+AB8BdhQsrwHWVLtfIzS22cCBguXDwIz0egZwOL3+P8CnirWrxQewDfidRhgvcB2wF7iT/IVRE1L90s81sAP4SHo9IbVTtfte4jhbyP9CvBv4FqA6H++PgZsG1cbVz3Pd7zkAM4FjBcu9qVaP3hcRJwDS882pXjf/BukQwoeB3dTxeNMhlleBU8BO4EfAmYi4kJoUjunSeNP6s8CNY9vjin0F+ALwy7R8I/U93gD+SdIeSatSbVz9PDfCd0irSK3RpmjVxb+BpF8B/h74o4j4V6nYsPJNi9RqarwRcRH4kKSpwD8Av16sWXqu6fFK+jhwKiL2SPrYQLlI07oYb/LRiDgu6WZgp6QfDtG2KuNthD2HXmBWwXILcLxKfRltP5E0AyA9n0r1mv83kHQN+WD4RkQ8k8p1O94BEXEG+B75cy1TJQ38QVc4pkvjTeuvB/rHtqcV+SjwCUk/BjaTP7T0Fep3vETE8fR8inz438E4+3luhHB4BWhNMx8mAiuA7VXu02jZDrSn1+3kj80P1O9Psx4WAGcHdl9rgfK7CBuBQxHxNwWr6nW8zWmPAUmTgUXkT9Q+D9yXmg0e78C/w33Ac5EOTteCiFgTES0RMZv8/5/PRcSnqdPxSpoi6T0Dr4E24ADj7ee52idmxujkz73AP5M/bvvFavdnhMb0FHACeJf8XxYryR933QUcSc83pLYiP2PrR8BrQK7a/S9xrP+J/G70fuDV9Li3jsc7D9iXxnsA+ItUnwu8DPQAfwdcm+qT0nJPWj+32mOoYOwfA75Vz+NN4/pBehwc+J003n6efYW0mZllNMJhJTMzK5HDwczMMhwOZmaW4XAwM7MMh4OZmWU4HMzMLMPhYGZmGQ4HMzPL+DfIR9lmVMbd/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0f6955b9b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Performance test on the Aircraft dataset\n",
    "\n",
    "groups = 3\n",
    "data = (np.array([(4.0,1940),(9,2960),(9,4630),(78,1528),(90,2040),(50,3700),(467,14815),(509,15200),(290,15700),(215,6045)])).T\n",
    "initial_centers = (data.T)[np.random.choice((data.T).shape[0], groups), :]\n",
    "cluster(data, groups, initial_centers, 1e-6, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHeNJREFUeJzt3X+QHPV55/H3o5VAAjtwh9YnGSkILmf+iKPEYo2FuUpcSIE4kkXqbH64ymfkc4Czz4d8duE6UipHR6kqcXF2LIVybAlsg3EAoWAskDli4cslNkjFShDkWMaFZWRJlqJFBAhYYLT73B8zs7sazcy3Z7un+9s9n1eVSjvf6Z5+uhmeHXV/5tvm7oiISLVMK7oAERHJnpq7iEgFqbmLiFSQmruISAWpuYuIVJCau4hIBam5i4hUkJq7iEgFqbmLiFTQ9KI2PHv2bF+wYEFRmxcRKaWdO3c+7+6DoeUSN3czGwCGgYPuvrzpuZXALcDB+tCt7n5bp9dbsGABw8PDSTcvIiKAme1Lslw3n9xXAXuAX2vz/L3u/okuXk9ERHok0Tl3M5sHLAM6fhoXEZE4JL2g+kXgM8BYh2Xeb2ZPm9lmM5vfagEzu87Mhs1seGRkpNtaRUQkoWBzN7PlwBF339lhsQeBBe6+ENgG3NFqIXff4O5D7j40OBi8HiAiIlOU5JP7xcAKM3sOuAe4xMzumryAux9199frDzcCF2RapYiIdCXY3N39Jnef5+4LgKuB77n7hyYvY2ZzJz1cQe3Cq4iIFGTKOXczuxkYdvctwA1mtgI4DrwArMymPBERmQor6jZ7Q0NDrpy7NGzdu5V1u9Zx+NXDzDl9DqsWrWLZecuKLkskOma2092HQssV9g1VkYate7ey5rE1vDb6GgCHXj3EmsfWAKjBi0yR5paRwq3btW68sTe8Nvoa63atK6gikfJTc5fCHX71cFfjIhKm5i6Fm3P6nK7GRSRMzV0Kt2rRKmYOzDxhbObATFYtWlVQRSLlpwuqUrjGRVOlZUSyo+YuUVh23jI1c5EM6bSMiEgFqbmLiFSQmruISAWpuYuIVJCau4hIBam5i4hUkJq7iEgFqbmLiFSQmruISAWpuUsmtu7dyqWbL2XhHQu5dPOlbN27teiSRPqaph+Q1HSzDZH46JO7pKabbYjER81dUtPNNkTio+YuqelmGyLxUXOX1HSzDZH46IKqpKabbYjER81dMqGbbYjERadl+oAy6CL9R5/cK04ZdJH+pE/uFacMukh/UnOvOGXQRfqTmnvFKYMu0p/U3CtOGXSR/qQLqhWnDLpIf0rc3M1sABgGDrr78qbnTgXuBC4AjgJXuftzGdYpKSiDLtJ/ujktswrY0+a5jwL/4u6/AfwF8Lm0hYk0U15fJLlEzd3M5gHLgNvaLHI5cEf9583AEjOz9OWJ1DTy+odePYTj43l9NXiR1pJ+cv8i8BlgrM3zZwP7Adz9OPAScFbq6kTqlNcX6U6wuZvZcuCIu+/stFiLMW/xWteZ2bCZDY+MjHRRpvQ75fVFupPkk/vFwAozew64B7jEzO5qWuYAMB/AzKYDZwAvNL+Qu29w9yF3HxocHExVuPQX5fVFuhNs7u5+k7vPc/cFwNXA99z9Q02LbQGuqf/8gfoyJ31yF5kq5fVFujPlnLuZ3QwMu/sW4HbgG2b2LLVP7FdnVJ8IoLy+SLesqA/YQ0NDPjw8XMi2RUTKysx2uvtQaDl9Q1WC1m5fy30/uY8xH2OaTeOKt13B6sWriy5LRDpQc5eO1m5fy73P3Dv+eMzHxh+rwYvESxOHSUf3/eS+rsZFJA5q7tLRmLf+3lq7cRGJg5q7dDTNWr9F2o2LSBz0f6h0dMXbruhqXETioAuq0lHjoqnSMiLlopy7iEiJKOfeJ6595Fq2H94+/njxnMVsvGxjgRVNzda9W/XtU4nT05vg0ZvhpQNwxjxY8llYeGV+60+RzrmXWHNjB9h+eDvXPnJtQRVNjeZql2g9vQkevAFe2g947e8Hb6iN57F+CmruJdbc2EPjsdJc7RKtR2+GN46dOPbGsdp4HuunoOYuhdNc7RKtlw50N571+imouUvhNFe7ROuMed2NZ71+CmruJbZ4zuKuxmOludolWks+CzNmnTg2Y1ZtPI/1U1BzL7GNl208qZGXMS2z7LxlrHn3GuaePhfDmHv6XNa8e43SMlK8hVfC+9bDGfMBq/39vvXJ0y5p109BOXcRkRJJmnPXJ/eS27p3K5duvpSFdyzk0s2XTik+GHqNLLYhIvnSl5hKrJEPb8QIG/lwIPEpjdBrZLENEcmfPrmXWBb58NBrKIMuUk5q7iWWRT489BrKoIuUk5p7iWWRDw+9hjLoIuWk5l5iWeTDQ6+hDLpIOemCaok1LmimmU0x9BpZbENE8qecu4hIiSjnnoEY8t0x1CDSE09vgr94O6w5s/Z3DtPg9hOdlmkjhnx3DDWI9ERjnvPGdLiNec4hl6/m9wN9cm8jhnx3DDWI9ESB85z3CzX3NmLId8dQg0hPFDjPeb9Qc28jhnx3DDWI9ESB85z3CzX3NmLId8dQg0hPFDjPeb/QBdU2Ysh3x1CDSE80Lpo+enPtVMwZ82qNXRdTM6Ocu4hIiSTNuQc/uZvZTODvgVPry2929z9tWmYlcAtwsD50q7vf1m3RcrK129dy30/uY8zHmGbTuOJtV7B68erEz0MtUtnrT/95bENEkktyWuZ14BJ3f8XMZgDfN7OH3X1703L3uvsnsi+xf63dvpZ7n7l3/PGYj40/Xr14dfB5yCcrrzy+SHyCF1S95pX6wxn1P8Wcy+kz9/3kvo7joechn6y88vgi8UmUljGzATN7CjgCfNfdd7RY7P1m9rSZbTaz+W1e5zozGzaz4ZGRkRRl94cxH+s4Hnoe8snKK48vEp9Ezd3dR939d4B5wIVm9vamRR4EFrj7QmAbcEeb19ng7kPuPjQ4OJim7r4wzVr/52mMh56HfLLyyuOLxKernLu7vwj8HfAHTeNH3f31+sONwAWZVNfnrnjbFR3HQ89DPll55fFF4pMkLTMIvOHuL5rZLGAp8LmmZea6+6H6wxXAnswr7UONi6Lt0jCh5yGfrLzy+CLxCebczWwhtdMsA9Q+6W9y95vN7GZg2N23mNmfUWvqx4EXgI+5+487va5y7iIi3Uuac9eXmDrIIrudJIeeZv0kNabdj7T7EI2nN6X/RmToNbLYhkgHmX2JqV9lkd1OkkNPs36SGtPuR9p9iEYW84eHXkNzlEtENHFYG1lkt5Pk0NOsn6TGtPuRdh+ikcX84aHX0BzlEhE19zayyG4nyaGnWT9JjWn3I+0+RCOL+cNDr6E5yiUiau5tZJHdTpJDT7N+khrT7kfafYhGFvOHh15Dc5RLREr2f2h+sshuJ8mhp1k/SY1p9yPtPkQji/nDQ6+hOcolIrqg2kYW2e0kOfQ06yepMe1+pN2HaGQxf3joNTRHuUREUUgRkRJJGoXUaRkRkQrSaZnIhb6ApJtkROahT8HOr4OPgg3ABSth+Rf6tw4pjJp7xEJfQNJNMiLz0Kdg+PaJxz468TjPxhpLHVIonZaJWOgLSLpJRmR2fr278V6JpQ4plJp7xEJfQNJNMiLjo92NV70OKZSae8RCX0DSTTIiYwPdjVe9DimUmnvEQl9A0k0yInPByu7GeyWWOqRQuqAasdAXkHSTjMg0LlYWnVKJpQ4plL7EJCJSIn0/n3sW+e/Qa+RxEwvl2LtQlhtlhDLoeexHBjcdeeDJg9zyyDP84sVjvPXMWdx42fn80TvOTlxC2vWls0o29yzy36HXyOMmFsqxd6EsN8oIZdDz2I8MbjrywJMHuen+3Rx7o5bAOfjiMW66fzdAogaddn0Jq+QF1Szy36HXyOMmFsqxd6EsN8oIZdDz2I8MbjpyyyPPjDfmhmNvjHLLI88kKiHt+hJWyeaeRf479Bp53MRCOfYulOVGGaEMeh77kcFNR37x4rGWi7QbT7pc0vUlrJLNPYv8d+g18riJhXLsXSjLjTJCGfQ89iODm4689cxZLRdpN550uaTrS1glm3sW+e/Qa+RxEwvl2LtQlhtlhDLoeexHBjcdufGy85k148RfVLNmDHDjZecnKiHt+hJWyQuqWeS/Q6+Rx00slGPvQllulBHKoOexHxncdKRx0XOqaZe060uYcu4iIiXS9zn3LOSRlb/2kWvZfnj7+OPFcxaz8bKNme2DRCiHHPumr36ei/d9ibk8zyFm84NzPs6V/+XTidf/6deu55x9mxjwMUZtGvvOuZJ//5GvZFqj9FYlz7lnoZExP/TqIRwfz5hv3bs1s9dobuwA2w9v59pHrs10XyQijQz5S/sBn8iQP70ps01s+urnWb7vzznbnmeawdn2PMv3/Tmbvvr5ROv/9GvXc95z9zCdMcxgOmOc99w9/PRr12dWo/SemnsbeWTlmxt7Q7txqYAccuwX7/sSp9mvThg7zX7Fxfu+lGj9c/ZtwuzEMbPauJSHmnsbeWTlpQ/lkGOfy/Ntxo8mWn+gzXc12o1LnNTc28gjKy9Tt3TpUsxs/M/SpUuLLimZHHLsh5jdZvysROuPtvmuRrtxiZP+a7WRR1Z+8ZzFLddrNy41S5cu5dFHHz1h7NFHHy1Hg88hx/6Dcz7OL/2UE8Z+6afwg3M+nmj9fedcSXOIzr02LuWhKGQHSsvEyZpPCE9S1Pu5K7mmZY5yiLOUlqmQpFHIYHM3s5nA3wOnUotObnb3P21a5lTgTuAC4Chwlbs/1+l1y9DcJU6lb+4iKWSZc38duMTdXzGzGcD3zexhd58c6fgo8C/u/htmdjXwOeCqKVWeUOgTcSzzoIfmfC/FfmTxSTOvOcz/179NtY0ntnyF+btu4S0+whEbZP+iG3nnipwjgHesgJ/9v4nH5/4eXLNl/GGSedCD+5HDvx5CdcYyn3tZ6uxWV6dlzOw04PvAx9x9x6TxR4A17v64mU0HDgOD3uHF03xyb57nHGrnste8ew3LzlsWfD4vzXO+N1x1/lWsXry6HPvRPLc31M4Rv2998mbQPId5w9BHW89hHthGq3PuAEvONbZ9+M1T3sYTW77C23euZtakGOExP4UfXrA2vwbf3Ngb6g2+eR50qM3J8mf/6bfGG05wP7L4bxoQqjPJfuShLHVOlvSTe6ILqmY2YGZPAUeA705u7HVnA/sB3P048BIkvDQ/BaH8eCzzoIfmfC/FfmSRy854DvNt27axZMmSE8aWnDvt5Mbe5Tbm77rlhIYIMMt+xfxdt7SuvxdaNfZJ40nmQQ/uRw5Z+1CdscznXpY6pyJRc3f3UXf/HWAecKGZvb1pkVYnQU/61G5m15nZsJkNj4yMdF9tXSg/Hku+PDTneyn2I4tcdg/mMN+2bRvuPv5n24fflHobb/HW78m3eOvceBGSzIMe3I8csvahOmOZz70sdU5FV1FId38R+DvgD5qeOgDMB6ifljkDeKHF+hvcfcjdhwYHB6dUMITz47Hky0NzvpdiP7LIZecxh3kG2zhird+TR6x1brwISeZBD+5HDln7UJ2xzOdeljqnItjczWzQzM6s/zwLWAr8uGmxLcA19Z8/AHyv0/n2tEL58VjmQQ/N+V6K/cgil53HHOYZbGP/ohs51pQPP+ansH/RjcnrSOvc3+s4nmQe9OB+5JC1D9UZy3zuZalzKpKkZeYCd5jZALVfBpvc/SEzuxkYdvctwO3AN8zsWWqf2K/uWcWE5zmPZR700JzvpdiPLOYXz2MO8wy28c4V1/ME1FMmz3PEZrP/gpzTMtds6ZiWSTIPenA/cpgzPlRnLPO5l6XOqdCXmERESqTv53OPIh8uE0K56hxy10m2kTrznME2Qsqau56KftrXrFWyuTfnwxvzqANq8EVozlU35jCHWuMLPZ9HDZyceT744jFuun83QMvMc/PzWWwjJO36ZdJP+9oLlZw4LIp8uEwI5apzyF0n2UbqzHMG2wgpc+66W/20r71QyeYeRT5cJoRy1TnkrpNsI3XmOYNthJQ5d92tftrXXqhkc48iHy4TQrnqHHLXSbaROvOcwTZCypy77lY/7WsvVLK5R5EPlwmhXHUOuesk20idec5gGyFlzl13q5/2tRcqeUE1iny4TAjlqnPIXSfZRurMcwbbCClz7rpb/bSvvaCcu4hIifR9zl0iE8h/r35gN3fv2M+oOwNmfPBd81n7R7/V1WtkIVhHoIZEuezA3PZ5ZLurkh/XsWpPzV16L5D/Xv3Abu7a/vPxxUfdxx+PN9YcsvDBOgI1JMplN89t76MTj5d/IZdsd1Xy4zpWnVXygqpEJpD/vnvH/parnTCeQxY+WEeghkS57MDc9nlku6uSH9ex6kzNXXovkP8ebXPd54TxHLLwwToCNSTKZQfmts8j212V/LiOVWdq7tJ7gfz3QJsbXp8wnkMWPlhHoIZEuezAvPN5ZLurkh/XsepMzV16L5D//uC75rdc7YTxHLLwwToCNSTKZQfmnc8j212V/LiOVWe6oCq9F8h/Ny6adkyp5JCFD9YRqCFRLjsw73we2e6q5Md1rDpTzl1EpESUc5cJecyVHpBFVvjx9St559FvM8AYo0zjibMu56Ibvp5rDVlso6y5aSkXNfeqy2Ou9IAsssKPr1/J4qPfonFtczpjLD76LR5fT6IGH0smusy5aSkXXVCtujzmSg/IIiv8zqPfpjnMYlYbz6uGLLZR5ty0lIuae9XlMVd6QBZZ4QHGuhrvRQ1ZbKPMuWkpFzX3qstjrvSALLLCo23equ3Ge1FDFtsoc25aykXNverymCs9IIus8BNnXU5zsMu9Np5XDVlso8y5aSkXXVCtujzmSg/IIit80Q1f5/H1TDktE0smusy5aSkX5dxFREpEOXfJVdrsdiz5cGXQpVlZ3xNq7pJa2ux2LPlwZdClWZnfE7qgKqmlzW7Hkg9XBl2alfk9oeYuqaXNbseSD1cGXZqV+T2h5i6ppc1ux5IPVwZdmpX5PaHmLqmlzW7Hkg9XBl2alfk9oQuqklra7HYs+XBl0KVZmd8TyrmLiJRIZjl3M5sP3AnMAcaADe6+rmmZ9wDfBn5WH7rf3fObdrDMUs61HksGN1RHsM4I5pzPwuoHdne+o1ROYnlfSHGSnJY5Dnza3XeZ2ZuBnWb2XXf/UdNy/+Duy7MvscJSzrUeSwY3VEewzgjmnM/C6gd2c9f2n48/HnUff5xng4/lfSHFCl5QdfdD7r6r/vO/AnsAvUOykHKu9VgyuKE6gnVGMOd8Fu7esb+r8V6J5X0hxeoqLWNmC4B3ADtaPH2Rmf2jmT1sZr/ZZv3rzGzYzIZHRka6LrZyUs61HksGN1RHsM4I5pzPwmib61ftxnsllveFFCtxczezNwF/A3zS3V9uenoXcI67/zbwl8ADrV7D3Te4+5C7Dw0ODk615upIOdd6LBncUB3BOiOYcz4LA823igqM90os7wspVqLmbmYzqDX2b7r7/c3Pu/vL7v5K/efvADPMbHamlVZRyrnWY8nghuoI1hnBnPNZ+OC75nc13iuxvC+kWEnSMgbcDuxx9y+0WWYO8M/u7mZ2IbVfGkczrbSKUs61HksGN1RHsM4I5pzPQuOiadFpmVjeF1KsYM7dzP4j8A/Abhi/YeWfAL8O4O5fNrNPAB+jlqw5BnzK3R/r9LrKuYuIdC+znLu7fx/oeNLQ3W8Fbk1enjRUJY8cS75bRGo0/UCBqpJHjiXfLSITNHFYgaqSR44l3y0iE9TcC1SVPHIs+W4RmaDmXqCq5JFjyXeLyAQ19wJVJY8cS75bRCbogmqBqpJHjiXfLSITNJ+7iEiJJM2567SMiEgF9fVpma17t7Ju1zoOv3qYOafPYdWiVSw7b1nRZZ2gLF9yKkudvabjILHo2+a+de9W1jy2htdGXwPg0KuHWPPYGoBoGnxZvuRUljp7TcdBYtK3p2XW7Vo33tgbXht9jXW71rVZI39l+ZJTWersNR0HiUnfNvfDrx7uarwIZfmSU1nq7DUdB4lJ3zb3OafP6Wq8CGX5klNZ6uw1HQeJSd8291WLVjFzYOYJYzMHZrJq0aqCKjpZWb7kVJY6e03HQWLStxdUGxdNY07LlOVLTmWps9d0HCQm+hKTiEiJZHazDpGqyOKGIsqxS1mouUtfyOKGIsqxS5n07QVV6S9Z3FBEOXYpEzV36QtZ3FBEOXYpEzV36QtZ3FBEOXYpEzV36QtZ3FBEOXYpE11Qlb6QxQ1FlGOXMlHOXUSkRHSzDhGRPqbmLiJSQWruIiIVpOYuIlJBau4iIhWk5i4iUkFq7iIiFRRs7mY238z+r5ntMbN/MrOTblVkNevN7Fkze9rMFvWmXBERSSLJN1SPA592911m9mZgp5l9191/NGmZ9wL/of7nXcBf1f+WlDR/uIhMRfCTu7sfcvdd9Z//FdgDNHeXy4E7vWY7cKaZzc282j7TmD/84IvHcCbmD3/gyYNFlyYikevqnLuZLQDeAexoeupsYPLE2Ac4+ReAdEnzh4vIVCVu7mb2JuBvgE+6+8vNT7dY5aRJa8zsOjMbNrPhkZGR7irtQ5o/XESmKlFzN7MZ1Br7N939/haLHAAmz506D/hF80LuvsHdh9x9aHBwcCr19hXNHy4iU5UkLWPA7cAed/9Cm8W2AB+up2YWAy+5+6EM6+xLmj9cRKYqSVrmYuA/A7vN7Kn62J8Avw7g7l8GvgP8IfAs8EvgI9mX2n80f7iITJXmcxcRKRHN5y4i0sfU3EVEKkjNXUSkgtTcRUQqSM1dRKSC1NxFRCqosCikmY0A+wrZ+MlmA88XXURAGWqEctRZhhpBdWapDDVCsjrPcffgV/wLa+4xMbPhJLnRIpWhRihHnWWoEVRnlspQI2Rbp07LiIhUkJq7iEgFqbnXbCi6gATKUCOUo84y1AiqM0tlqBEyrFPn3EVEKkif3EVEKqivmruZDZjZk2b2UIvnVprZiJk9Vf/zxwXV+JyZ7a7XcNK0mfU589eb2bNm9rSZLYq0zveY2UuTjudnC6jxTDPbbGY/NrM9ZnZR0/OxHMtQnYUeSzM7f9K2nzKzl83sk03LFH4sE9ZZ+PuyXsf/MLN/MrMfmtndZjaz6flTzeze+vHcUb/FaXfcvW/+AJ8C/hp4qMVzK4FbI6jxOWB2h+f/EHiY2q0NFwM7Iq3zPa2Oc8413gH8cf3nU4AzIz2WoToLP5aTahkADlPLWkd3LBPUWfixpHZ/6Z8Bs+qPNwErm5b5OPDl+s9XA/d2u52++eRuZvOAZcBtRdeS0uXAnV6zHTjTzOYWXVRszOzXgN+ldhcx3P1X7v5i02KFH8uEdcZkCfBTd2/+AmLhx7JJuzpjMR2YZWbTgdM4+bakl1P7pQ+wGVhSvyteYn3T3IEvAp8Bxjos8/76Pyk3m9n8Dsv1kgN/a2Y7zey6Fs+fDeyf9PhAfSxvoToBLjKzfzSzh83sN/MsDjgPGAG+Vj8Vd5uZnd60TAzHMkmdUOyxnOxq4O4W4zEcy8na1QkFH0t3Pwj8b+DnwCFqtyX926bFxo+nux8HXgLO6mY7fdHczWw5cMTdd3ZY7EFggbsvBLYx8Vszbxe7+yLgvcB/M7PfbXq+1W/vIiJPoTp3Ufsn8W8Dfwk8kHN904FFwF+5+zuAV4H/2bRMDMcySZ1FH0sAzOwUYAVwX6unW4wVEsUL1Fn4sTSzf0Ptk/m5wFuB083sQ82LtVi1q+PZF82d2n1gV5jZc8A9wCVmdtfkBdz9qLu/Xn+4Ebgg3xLH6/hF/e8jwLeAC5sWOQBM/lfFPE7+J13Phep095fd/ZX6z98BZpjZ7BxLPAAccPcd9cebqTXR5mWKPpbBOiM4lg3vBXa5+z+3eC6GY9nQts5IjuVS4GfuPuLubwD3A+9uWmb8eNZP3ZwBvNDNRvqiubv7Te4+z90XUPvn2vfc/YTflE3nB1cAe3IssVHD6Wb25sbPwKXAD5sW2wJ8uJ5OWEztn3SHYqvTzOY0zhGa2YXU3mtH86rR3Q8D+83s/PrQEuBHTYsVfiyT1Fn0sZzkg7Q/1VH4sZykbZ2RHMufA4vN7LR6LUs4ud9sAa6p//wBaj2rq0/u01OXWWJmdjMw7O5bgBvMbAVwnNpvyJUFlPTvgG/V33vTgb929/9jZv8VwN2/DHyHWjLhWeCXwEcirfMDwMfM7DhwDLi62zdnBv478M36P9P3Ah+J8FgmqbPwY2lmpwG/D1w/aSy6Y5mgzsKPpbvvMLPN1E4RHQeeBDY09aPbgW+Y2bPU+tHV3W5H31AVEamgvjgtIyLSb9TcRUQqSM1dRKSC1NxFRCpIzV1EpILU3EVEKkjNXUSkgtTcRUQq6P8Dv5jENNA8iZ0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0f692925c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Performance test on the famous iris dataset\n",
    "\n",
    "from sklearn.datasets import load_iris\n",
    "\n",
    "dataset = load_iris()\n",
    "groups = 3\n",
    "data = (np.array(dataset.data)).T\n",
    "initial_centers = (data.T)[np.random.choice((data.T).shape[0], groups, replace=False), :]\n",
    "cluster(data, groups, initial_centers, 1e-6, 10)"
   ]
  }
 ],
 "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
}