SpectralClusteringIntroduction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tnotargi/anaconda3/lib/python3.6/site-packages/matplotlib/font_manager.py:280: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  'Matplotlib is building the font cache using fc-list. '\n"
     ]
    }
   ],
   "source": [
    "#Documentation : http://docs.scipy.org/doc/scipy/reference/sparse.html\n",
    "\"\"\"\n",
    "Sparse matrices can be used in arithmetic operations: they support addition, subtraction, multiplication, division, and matrix power.\n",
    "Advantages of the CSR format\n",
    "\n",
    "        efficient arithmetic operations CSR + CSR, CSR * CSR, etc.\n",
    "        efficient row slicing\n",
    "        fast matrix vector products\n",
    "Disadvantages of the CSR format\n",
    "\n",
    "        slow column slicing operations (consider CSC)\n",
    "        changes to the sparsity structure are expensive (consider LIL or DOK)\n",
    "\"\"\"\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.sparse import dia_matrix\n",
    "import scipy.sparse.linalg as spl\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SpectralGraphTheory\n",
    "\n",
    "## Introduction\n",
    "This notebook intends to show some experimental work about spectral graph theory, and in this case mainly about spectral clustering\n",
    "We will try to reproduce some of the interesting principle presented in the paper by Andrew Ng (On spectral clustering: Analysis and an algorithm.)\n",
    "\n",
    "## Graph as matrices\n",
    "In this notebook we will only be using undirected graphs, whose structure can be represented using an adjacency matrix $A \\in \\mathbb{R}^{n\\times n}$ that can be designed as follow for, $n$ vertices $s_i, i = 0,1,\\dots n-1$:\n",
    "$$\n",
    "    A_{ij} =\n",
    "    \\begin{cases}\n",
    "        e^{-\\frac{\\|s_i-s_j\\|}{2\\sigma^2}},& \\text{if } i\\neq j\\\\\n",
    "        0              & \\text{otherwise}\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "For any vertex $s_i$ we define its degree $d_j$ as the sum of the $j^{th}$ row of the matrix $A$.\n",
    "From this we can also define $D = diag(d_0, d_1,\\dots,d_{n-1})$ \n",
    "\n",
    "## Laplacian and normalized laplacian\n",
    "\n",
    "We can define the Laplacian operator , which in our case is also known as \"normalized laplacian\" L:\n",
    "$$\n",
    "    L = D^{-\\frac{1}{2}}AD^{-\\frac{1}{2}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fe07e348ef0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+MHdd13793H3cBL6km4iMt07L3UQbUAnKRGPZCsA3D\nVWI3odkGSgLXsLFWVMPAgssYaFO4AIUFEqHAoon9R8vEVlRWYExrNxYCFI6FhI4iKQGc1nDipaGf\ndWnLCpeW4NYSaVuhVwFl8vaPmfHOm70/zp25M3Nm5nyAwXtvfry5M3PnnHvPOfdcpbWGIAiCIGTM\ntF0AQRAEgReiGARBEIQpRDEIgiAIU4hiEARBEKYQxSAIgiBMIYpBEARBmEIUgyAIgjCFKAZBEARh\nClEMgiAIwhR72i5AGQ4cOKAPHz7cdjEEQRA6xblz517WWh/07RdFMSilTgP41wC+r7X+54btCsBJ\nAEcBbAP4t1rrb6TbjqTbRgAe0Fr/ru98hw8fxubmZoyiC4IgDAal1BZlv1impM8BOOLY/gEAt6bL\nMoA/BACl1AjAZ9PttwH4iFLqtkhlEgRBEEoQRTForb8C4LJjlzsBfF4nfA3AzyqlDgG4HcBzWuvn\ntdZXATyU7isIgiC0RFPO55sBfDf3+4V0nW29IAiC0BKdiUpSSi0rpTaVUpsvvfRS28URBEHoLU0p\nhhcBvDn3+03pOtv6XWitT2mtF7XWiwcPep3qgiAIQkmaUgwPA/gNlfBOAD/SWn8PwNcB3KqUukUp\nNQfgw+m+giAIQkvEClf9AoA7ABxQSr0A4HcAzAKA1vp+AGeRhKo+hyRc9WPptp8opT4B4BEk4aqn\ntdbPxiiTIAiCUI4oikFr/RHPdg3gNy3bziJRHILQOzY2gNVV4OJFYGEBWFsDlpbaLpUguOmM81kQ\nOLGxARw+DMzMJJ8bG+Z9lpeBrS1A6+Rzedm8ryBwQhSDIBDIK4IDB4CPfcwv8FdXge3t6XXb28n6\nmOWxKSZBKIsoBqE3UFvxoQL1+HHgrrt2FMGlS8Brr03vUxT4GxvJ/iYuXqRdjw3piQi1o7Xu3PKO\nd7xDC8NgfV3ryURrpZLP9XX7fvPzWieiMlnm56f3p+xj+l+lpo+xLUrZz5NfJpNq92QyMf/vaOS/\nT8KwAbCpCTK2dSFfZhHF0F/yimA81np2libIbcIyL4Qp+1D/1/U/rmN8iogCRVHFOI/QP6iKQUxJ\nAhuKJhKKySbDZp7Z2toxHZUx7VDNPvPzScSR75hTp6pHJS0s+Pcx3SfxSwhURDEIbDA5a03kBW8m\n7LQ276vUjqKx4RK0rm3jcfL/k8m0wLcdM5nECVVdW0sUkY/ifRK/hEBFFINAoonWJrV1ngnevLAz\noZRbIQDTLX0TJiGsFLCyArz8MnD9OnDhwrTANx3jO08IS0uJIppMkrKMRub98gqqzggpoYdQ7E3c\nFvExNIvJmTo7m/gAYjo7Kfb8vO3ctb/vv0LKTXWAU48p83++cxWfz9zc9PPxOcyFYQBxPguxCBXY\nZVhfTwRZ8X+LAi5/DpsTNhN2tnIr1Z5jtkxkFPV/XU57272qGiEldAtRDEI0qOGamZAJbRHbwjvH\nY/exvigjV6hpWwKxTGRU6P10KUSJXBo2ohiE0hQFkakl72qJh7aIy8blU87FzYTi6+UUKXM/XYo8\npglL6B6iGIRS2OzVRdOETeiUaRFXicv3tabLlKdOQssTc+yFmI0EUQwCmbxwtQnp8Xjahj03Zxbc\nvhaxSZBTB5GVEWx12fTLElqe0B5Gdo7i85mbkx6CQFcMEq46cIrx7Vqb97t8OQnLvH49CdM8fXon\nXDIfx2+L4V9YsMfSHz0aHpdPpRjaWRxz0DSh5XHdTxfF52h7roAMfBMMULQHt0V6DPGI3Vp3tYhd\nJo58T2I0itdj6DoxfTam+8etRyXUC8SUJFCgRhyFCAqb3Z9qFhFhNU1oVFKI+Un8EcOCqhhUsm+3\nWFxc1Jubm20XoxccPmwfOZwxHifmo7rONZkkZqo8MvNZeULu88yM2cykVGI2FPqFUuqc1nrRt18U\nH4NS6ohS6rxS6jml1AnD9v+olHoiXZ5RSl1TSu1Pt11QSj2dbhNp3zC+vDvz88DJk/Wdy5YqYmlp\nx6dRTDkhuAm5zxQfRl3zXAiMoXQrXAuAEYDvAHgLgDkATwK4zbH/rwD4q9zvCwAOhJxTTEk0QuYy\nyEccxU51UaZMQjWo93llxWxKWlnZ+R9fuo2VFTH9dQU0ZUpSSr0LwL1a619Of9+TKpz/bNn/jwH8\ntdb6v6e/LwBY1FqTjRViSvKTRQDlE6fNz7cbkSPww2d2opgabckKRyPgzBmpb5xo0pR0M4Dv5n6/\nkK4zFWoewBEA/yO3WgN4TCl1Tim1HKE8AiSbpkDDFgKcraeECNvalteuSWrvrtL0OIZfAfC/tNaX\nc+veo7V+G4APAPhNpdR7TQcqpZaVUptKqc2XXnqpibJ2Gt8LLwiA3cewf797OxVpjHSTGIrhRQBv\nzv1+U7rOxIcBfCG/Qmv9Yvr5fQBfBHC76UCt9Smt9aLWevHgwYOVC90XbE6/sgOjhGGxtgbMze1e\n/8orSV2iTgqklH2bNEa6RwzF8HUAtyqlblFKzSER/g8Xd1JK/QyAfwHgS7l1e5VSN2TfAfwSgGci\nlGkQ2EYSHz8OXLmye/+Yk8UI/WBpCbjhht3rX3staekXR2qPx8Ds7PS+8/PAsWO0CYOEblBZMWit\nfwLgEwAeAfBNAH+itX5WKXVMKXUst+uvAfhLrfWPc+tuAvA/lVJPAvg7AH+utf6LqmUaCjY/wv33\nJ/Ml5xmPxfEsmLl82bw+a+nnQ4dffhn4oz/andLjvvsSR3NRaczOSmOki+yJ8Sda67MAzhbW3V/4\n/TkAnyusex7Az8cowxCxddFNzsB9+0QpCGYWFsyRR7aW/tKSvS4VTUouE5PAF0mi12FCuuhi5xVs\nxJqjenUVuHp1et3Vq3TnswyS44Mohg5jm6jehNh5BRuxMtCGRMIVlcDx42Z/mSiHdhDF0GGyF3o8\n3lk3P292DoqdtyI9b87GSEFCjYQzBU3cf7+Mu+GEKIaWMcmbUBn06qs733/8453oEQ7zD/QCV/hX\nj5VFKFSTlClowjZITkygLUHJm8Ft6UuuJFMemtlZ++xoxWNdcylI2uSIUCetkARBemVlZz6N0Wgn\n51Ieaqp3IMnJJMQDMoMbf0wtp9de2+3AK3ap8w1YG71saTVtzsnO50sWlGF6UAPqUWxsJCGr164l\nv69dS34XL1v8XR2Aoj24LX3pMYS0nPKTrFAasL3rMays7L5hthZ6jBSupu5cyIMa4GxD1El/Qm6t\na25rIRzIDG78oVooii+XT6H0Tv6sr9svmiJ1yiiQkIdjKs8Ap0YLmTmueOvH48HdrlYQxdAByvoY\nfL6FXikFrd0XXJQ6Ic3W4o2em/PPQ+pr3lLnMe3hxBRVdGFIh1AojyiGjmCSDz6Z0VsrRehk0Sap\nQ2222pqombcztMeg1LSn1SUlKQ+wg4qjbL00HVe8nUIcRDH0nA7KjWmKF+CaBswmZPMt9Axqs9Ul\n5LPyUQzhZbS3r4wd1vxl6uUArW6tIYpB4IutiRjSwrY1KalC1SXsi103am/FVJaQHlDWqxmYpAzx\nTQjVoCoGCVcVmidkhNPW1s7+WV7nyQR48MEkpWcR03Dw171u93757aZzLqeTCV64AKyvl0smlB9O\nvLaWXMeM45XL4jgHNsuSzB3CEIr24LZIj6HjhMbphppUKJ7M9fXE0+/rOeT3L2O7W1+3+zPyy+zs\nzn8OrMfQYctZ54CYkgS2uHwGVPOSjdDQ1pCIp1BCAvbzQ3wHKCk77zPrCFTFIKakBsgPgD1wIFkG\nMhjWjC0t7C/+4nSKT+0wL9lYXbUfVzTFZKaeycS8f1VbhslkZiM/W44t3SnQ25HUriR+AxtAzgOK\n9uC2dKnH4Gs0dr4hWLapRzH3hEQjZYSEtuavoY4WeojJjOLIjlHGjjXNbZe9stKpy2ADxJTEg16n\nr6girCh29BCzkO9/XcokO1dsSUMdC5EfWBf6Xy5lZxog0zETFdXqyPwy2CCKgQmURmNnw/Jsb+14\n7Bey1BjF0JvGabQU1ceQdzzbCM03YVIAHcw7EbPTJTSsGAAcAXAewHMAThi23wHgRwCeSJffph5r\nWrqkGCiNxtGoo11i6ls7O5sIpfxFUlvAZSJ0OJlLYiUFst0HU+UpM2qbKSGXwvgy2NCYYgAwAvAd\nAG8BMAfgSQC3Ffa5A8CflTnWtHRJMYQm6exUl7hsornMSEwxa3TQ/OGk7GguakWihMZ2qKntsiZ2\n6DLYQFUMMaKSbgfwnNb6ea31VQAPAbizgWM7QTHAZDzemV0tG6+Vp1PTGZqiiyhsbwNnz9ImGo41\nITEXyo7mKt4HU+UBgEuX7P8xHpcbqNciS0uJ2PehFHD0aP3lGQwU7eFaAHwQwAO533cB+ExhnzsA\nXAbwFIAvA3gr9djctmUAmwA2FxYW6lKojdKLVABUU0mnLzIisXpAodlfs3NwMrMRkQn04gFm4xi+\nAWBBa/1zAP4AwJ+G/oHW+pTWelFrvXjw4MHoBaxC2TjrXqYC+NCHaL2ITl9kBWL1gELv3913J+co\nDhgA2A8SsA17KdKp3jZ3KNrDtQB4F4BHcr/vAXCP55gLAA6UOVZrXj6GKg3A1s3nVVuPlCDz8Zg2\nibUQRqjzKp+51dWrY/psilXVdakMi88GNOh83gPgeQC3YMeB/NbCPm8AoNLvtwO4CEBRjjUtnBRD\n1bQ2rfXsqVrJVsD19Z1Z330X30HzRSeg5mHKTHfr61rv2UNXIoxxKQemuo0FjSmG5Fw4CuBbSCKM\nVtN1xwAcS79/AsCzqeD/GoB3u471LZwUg8tPwFoeUgeY2XoErtbqUP0HbeHL+ZQtMzN0JcIc08D5\njum2VmhUMTS9cFIMrjFerKMsKZ5vW2vU1lOQt7J9fBKTsuQT+mnNroVDsaJ1QLe1AlUxSBK9ipgc\nY9nvYv40Vs4xn+d7Y8Me+njtmv1/mYc/9p777kvmqrAlBgxlYyOZm2JrK5G52VwVLTqpKbkJhxrb\nEAtRDBXJB5kASXj59rZdprKZa8Wm0TKh7tJgthj60ajbYwz6QhZ5ZArdoZDP9GqSwi23cHzv0Px8\nMqaBebAVa0QxRGBpaUfOuhrTAKOWjC9s0vX2LS+blcqZM6IUOFG2suWPYzibnOuyJpMkMvfMmd2d\nnOPHRVmQodibuC2cfAwZFN8fKx+Dibwt2eZHyOzPzOzOggGTMX5ubidv1Xi8exY7aurzFv1I6+u7\nI6CBnVyEkpHVDsT53Cwufx9L2VkU7L5Io/ybJ3QHnwKnbDdlq82UQ0v1wZWLUDKy2hHF0ACUBjbL\niud62V0LZd4AoX/km+FMmt2uoDrJyGqHqhjEx1CSYrCGzbdw5QpDW6bJoai1/7irVxmFVQlRKOZz\nMRni81OgFuvJ9jbw0Y82brR3BdVRU2i4/mfwULQHt4VDj8HWKjGNIWJny6wS5z60JpbusTuFMiAg\nX3l99abBiu4buE+xlLJ7LxsAYkqql1DZysqkVHYeBXYXUg2KwG89n1WdUOtB9swp6TcarB+hCru3\nCj4AUQyBhFQaV5qgTjS0KS3FmZleJ7+jCnxXUE7nBQ21dZNVXopiYFXR7XT+2ZVEFEMAIa1Cl0xl\nNa0uNRrF9YL3+O2hPieXDDTVg337diJBi7OZsiOkx7C+Tt+XOb3uBXoQxRBASKi2a+rdbB6U1itd\nSCEYxqnXjU/GZc+xzCyZLiXCTvBQfQyUUGa2F7mbAVb5nyKKIYCQmdQo+7be0A7ptrDQZPUSOsnc\n7Kx5AFXVhSJ4Gq87Ji9tsQC+nmV2cR2pM72YObEkohgCiNFjYNPacDWHbTW/dU1mhyK3TPvl14fM\nZ1P34vNdsdTRLl9E64ULh/07XCOiGAKwjfdaWaHty+LlzXC17jpW86tYOrJnUiUAq47FVVfYCiy2\nBZsm/7yz4BDqYG5W73CNiGIIxJTG3uWAZtrA7lXrjirUXaPOq05NUMdik6dsTRwdkKS+oJDMb5TV\niU4EB9SAKIZAOtIo8uOaOahj1C3UYzqXAbpfwiboWddB1q0hfyOC/cRZDUFVDJISI4VhduFy2PIB\nXLrUuVzDdacr2LfPvT2fkXxlZfr3+963O83Cnj3AeOw/r9bmR+GbIqNVsrQY168nn8zSq/ve00uX\n2E0rwRuK9vAtAI4AOA/gOQAnDNuXADwF4GkAXwXw87ltF9L1T4CozaTH4IFh0rMy1O04diU/HI3c\nZQtJiWJbTI+CecOcLWV9Sa2b6RqGKmNjKIURgO8AeAuAOQBPAritsM+7AdyYfv8AgL/NbbsA4EDI\nOesa+dy7rmbHtJ1JKMYeT1C8DSsr5m2mwIM8scxco1GHlcDKyo5mHY38N61GOjPwtGWaVAzvAvBI\n7vc9AO5x7H8jgBdzv9kohnzlGY87+KJqPS1dO9RMcinmOiKL8kq/jHyjlqlqD6J1bF2Yshq1gaJm\nzzET/GwGnjKgScXwQQAP5H7fBeAzjv0/Wdj/71Mz0jkAy5RztpkSgzVU2wvDZpKrcxNbKWSj1KtQ\nl5nLl4OpUVOT68Vw2eCY2sOYFqtRWCoGAL8A4JsAxrl1N6efr0/NUO+1HLsMYBPA5sLCQtSb1TGL\nix2KFGWq8XwdnFiCN+bl5wVNaFJFXxlNZW68AVNWWxdvRot1TpTBNOxMSQB+LvVF/FPHf90L4JO+\nc7aZEoM1PgmaJcZjiCsHVYhQ3bPHbr6p0zxoE9qh/hHXVNuNzxLoejFCH0wLrSxbr66zZuIINKkY\n9gB4HsAtOefzWwv7LKQRS+8urN8L4Ibc968COOI7Z5spMVhD6TEwxSZYQ2WPyWHdlCCwOc+bSMlR\nC64Xw+ZjsC0ttLJcSplpx7l2GlMMyblwFMC30h7BarruGIBj6fcHAPwg9SX8NCw1jWR6Ml2ezY71\nLXX4GGZnpytOZ+a9Lw7nLF5IBxSDy2kYMvqZKy6TfIyltmt32a7W18M86w23snwZdFsoEgsaVQxN\nL3UohuKo1U7Me296cV3DbxmOfvbZzUNa3Jypu+dQa8FNRvqQqIAWmueU4nXOVBwBUQwBdNaU5Ep/\n0RFNR7n3FCcv+2elaa3YMksr106JCGjR40spXhfqTGyoikFSYqDD6TBsBbx8GTh9ejqHw+nT7NIY\nAPZL2NpK0kbMzCRpC9bWkmwMZ84wThvhYWkpeRQxae3afflKlGo1fYaveF2pM61B0R7clsH2GKgz\nzrAruB3bvXdl8uhyCGJMk1Kr1+67kJbroC2Vfuv3rWUgpiQ6rMPaXHmPTFONdSzcwpTu3GYG6JC+\nc+J6pKzNR0Vs+UqY1MEuNyDqQhRDICzrOKV5OR53rva7BGPHMnlUwtZb8qXvbl3uFiWubVo9gR2i\nGErAzqTUw9AKiq7rsoM5hJDewt69TORuT200Q+ldUBWDOJ9zsHNCU05c96QFkVld3Z0Xv8i1a911\nMIewfz9933/8x/rK4WVjYycS4O67dz9ArZPPrS1geRk4fnxnfwZzgOSLbyrOxkZS7K2t5FKyy+jQ\n1CXxoWgPbov0GNKlM6PwdqBGOQ7BOlE2nXijpqQy3nJGc4BQBq+ye+9rBGJKCoddllXTyLv8wnTA\nmiszKHUEcB9fyiJllELj9ydWetuWHqhN+eZfnd7kSiNAVQxiSsqxtAScOjUd/n/qVMvh/1rbt12+\n3Fw5CLi65Nm2a9do/8V+DEkEZiq8fY3dn1gnaumBXrrkX2+zxrKw0vrsYHVB0R7clrI9hs45mHyt\nNWbNaleXPLThyezSolN1FHSMOSVIxJo+r6UH6ipSBjtLQY0Fg5iSpmH78F24DPIMC+/qkofIEIaX\nFp0YFpra75PJQA/Q4mmZPFCKKSm7VHaNxhqcH6IYClDvMasK4pqkgEXNnSbGnAosBhU2QKzJh2pt\niLtycblG3TN6iTqbIFPrWpwfohgKUO4xu16FrUBMQ3aqpHtoeS753dQs3GJZaGp1kLq0l60nwaQu\nZpjm5mD6+uxGegw8egwsw9ZMo0xZaa9pqkx32cpl2GbXqbmZ2QnTve2FcE2RxwhTQ6VTWWTEx1C/\nYqA0vltplYXCUnvtpqxztdHLMCVqci0RBV8MU1IjPoaQ+T5YvShhfhxmr88OkXuuohgMUBrf7CtN\nB4Kuq5iUGruM9fVy0jkSVZ3PjZlAqBl92b0oYY+X0etTK6IYCFBeThbmjXy3hnEioZCpOFuXLWUL\nGglfj2p+nkFW9aKB3jeVJzN7TC96DJFpVDEAOALgPIDnAJwwbFcAfj/d/hSAt1OPNS2xFIPPfNSK\nY6pMc5uBkTTGPAONXkaZ3sLMTNQi+GRsq8EQvlH3xYWZf0Frep1k8Po0RmOKAcAIwHcAvAXAHIAn\nAdxW2OcogC+nCuKdAP6WeqxpqbvH0GrrIbQlmyUWapnOmEaqFjhiIV2nyWgt8jPk/jCWrL5eLNPI\n79poUjG8C8Ajud/3ALinsM9/A/CR3O/zAA5RjjUtsRQDu/BUrcu1ZBn0g8s6U1u732W7OBELTB18\n1QrUB8o63nMHlu+6jRpbA1TFECNX0s0Avpv7/UK6jrIP5djaYJkbqUyCFgaJhcrmlWnsfhdzzmQn\nD52EeXs7yR0egZMngdnZ6XWzs8n61qE80MkkyYW+usomxbaN/LsOAKPRzqNkVWQuOcAp2sO1APgg\ngAdyv+8C8JnCPn8G4D25348DWKQcm9u2DGATwObCwkI0DcqOMi1ZBj0GW4vMFcSiVIP2cltz0TZ1\nn2uJGMLCaJDwND4fw9wc+zE1JmxjG8bjlp+Bz+YV6R2HmJI6TEh4D6MX0TZWzDRItlGdZruXe/eW\nMycxUMSN4FKamTTt2P1hGYlIaQxGaow0qRj2AHgewC3YcSC/tbDPv8K08/nvqMealt4rhgxq2BTb\nZmeCKzSzkfjxWImJGu3mNASl7pQJiGgZ22WFuE4ag3p/uzbADUnU0beQRBitpuuOATiWflcAPptu\nfxrAoutY3zIYxeDrWrYe00in1QiwWJPNACwiwKLhM7H5UgKwkKphl0WtCo3qtpB7LCkxRDF4u5gs\nRkHRaD0mP4ZSYBEuVAHqKOayJjYmjRJXI4TqwmPZY4hQOFEMfaHscGIG3fkirVq8qmat6+D82lPE\nGIFoW0YjVmZMX9aYfD0cjxkk1bO1mmp4t6mKQab25M7SEnDhQhJPG0JD8xJSZh7M9rnrruT3gw8m\nl9RoWPDJk+H3ME+VYzmwuprEZ9bB9evJ0vhDNRMyVee+fcDHP95yyLotbt4WSt3Eu03RHtyWQfUY\nMlyTprRko6GYh0z7tDZYOzSbaqv2hcjEcMAzztOVx+c6YeGSs4Xw1ZxiH2JK6gl5U1Lx5TY5Chvs\nzleZ46K14J4y4xYidOFbJ3RuheKS5ahnIVX92F4JFmlwQiaKiDyrkCiGPmBrbmc1ueUXkmLLZdsA\nL6sgGNz3UtiaypReVP6amYdG+2CRtb7FtK+iGPoAi+aNnarRHywa4GUilpi2kr3ke5+ZWYgipDoS\noltlSEajr1SLE0WIYugDLJo3dmLEi7NodIYqhrwk6VoL2tRD8Amq0ajtUnvx+Q4oFtnGkB6DKIZK\nsGjeuKk6wpRFA7xsODAbTyaR9XV3Y8N1vYxZX3f7xdlZZEN8DJELJ4qhD3RN8OQIlbWtjh1bWQlX\nDDZJ1JTipvRUqAPaKNfKFJ/JMrt0du0rSlRSDe+5KIY+UHSQjsc0AcBAcZQZT9ValJJLyJdd6ngO\nIfaQ0AcwHicjnk3bGPsYfA0QV1YPJhbZRhHF0HWovQXGvYr1dbusYdGCK6O92rKRUcqa73aFdtlM\nKXBnZlgrBa3dj8Pn72JkkW0MUQxdh1qbmdf6EPnUeAsuVHiWGSAXMnekaYBT9pvaqynj5LHty6QO\nubA9wvxtZ9x2ahxRDF2H2v9l3k8OkU+Ny6EyhStjp6cMVIqVyyiTiDEyyjKpQy6oCWLHYwaT8TCA\nqhgkVxJXqAlfQhLDtMD+/fR9r1xpeAZD6j1SKpnCsizb28D997una4yVy+jateS/jx4F5uer/ReT\nOuTClmYImJ4h89Il4NVXW8rT1UUo2oPbMogeQw98DFq7MzubtjVa9JBopKyZGXPSn3wXKeb/5sub\nNZlDnD3M6lAZmFtY3dQYTAIxJfUAagVh3Gd2Wbpaf3lDzS1lBCzVVBNzMqHif2esrLh9FczSZ1MI\nHUfD1jrW0Ag8UQxDhGHvwSX8W395y7bSXZNYhyx5DUjxMZgGQVH+uwjDelKGMiPvG0+WR23Y+Z59\npIKLYhgCxZ4Co+gSSgOo9ZfXFdLieknH4zgDx0xmQdu5Z2Z2zuPLiEoR8szGvriKUyZTauu6L6QA\nlN5ipNZSI4oBwH4AjwL4dvp5o2GfNwP4awD/G8CzAP5dbtu9AF4E8ES6HKWcVxSDDotiabj/TE1B\nYLsE2zi+RgpKSfZUvJ9lEvHNzJjLVNXX0Fo+8/L4Iots20Jmamtc94W0eijPvEs9BgCfAnAi/X4C\nwO8Z9jkE4O3p9xsAfAvAbXpHMXwy9LyiGHSribjKFs1UjOLg7sZbdzbp4ZI8pgsp02swUdXX0DGl\noLW7vpTd1jpUOyll5H3XfAwAzgM4pHcUwHnCMV8C8C+1KIZqUAVFC3MVh/oO2L7gpgFttpc0NGLJ\nln+oyniGVhNOladMJ4l9DkNKpXY965qy/DWlGH6Y+67yvy37HwZwEcA/0TuKYQvAUwBOm0xRpmXw\niiFECLUgLFyNXlM9b90J7SLEHhE6Mtp1zlBJyUYihuOz2rnkKzNXyQ4UrUUZth2ZaIoBwGMAnjEs\ndxYVAYAfOP5nH4BzAH49t+4mACMAMwDWAJx2HL8MYBPA5sLCQi03rTOwzjPhb/RS349iKD6rF99G\nscC28FZfdyjUVMj+xtjx1ZfW508oi6/yttAiYmVKAjAL4BEA/8HxX4cBPEM57+B7DKzzTCT4/Lf5\nwJ7xeHe7ADzvAAAQpUlEQVQEaMemGLZT1t5BNSmx6FZVx9dJyszwoxH7vH50WrChNqUYPl1wPn/K\nsI8C8HkA/9Ww7VDu+28BeIhy3sErBluFYti0ouqwubnd4/LY+h5CKdvtyR/nmommJ3SoWsehBSdJ\nU4phDODxNFz1MQD70/VvBHA2/f4eADr1I0yFpQJ4EMDT6baH84rCtQxeMdgqlC9RWwtUCZ5i7Xto\nGtae1jisr8cZv9cpGraVNqIY2loGrxi07ozxPcRPXhT4vekxxKIjz7xIyABg6qDyQTYOIkBVDCrZ\nt1ssLi7qzc3NtoshEFGKtt9kkmS+zNjYSDJk5pOOzs/vZM9cXQUuXkySgK6tScZMjpieIQCMx8DJ\nk9PP7PDhJBsqhWJdEWgopc5prRd9+0nabaF2JhP/PvPz5szWr3vdzvfx2JxSeWsLuOsu4P3vT4TL\nzEzy2WgKb8GILZv4pUvJMzx+fOeZUZWCra50ho0N/hWV0q3gtogpiQAjs4PJPG5yNvuOoWSr6LEJ\nvpP4zIhlBrd1+pm27CuC+BgGDENHZaieKpOVVfwR/IidTbzzz7NlxxlVMYgpqWtQuqGm/vv2drK+\nJZaWEpvw9evuGbSyy7OZFTKfApWtLd499r6ztlZ+Irmib6rzJiQgqcAh61tCFEOXyDx5rikiAXfl\nY2zfzF+ejYWFZNbKEFy3SqiXbOrN8Xj3NltQwmSSPLMHH9w9ZSeLAIMq7xDzqXh/CqVbwW0ZrCmJ\n2g217TceszMx5fGZHUJ9DL00RTCFYios7tPJke1UM60ra6/4GEQxRCUkla+p8tlSQzORli7fQf7d\n8u3n80EIYfiEfhVZxyhGgkbZrKn5G9LiRYti6COhEx0UKx/zocRVO0T5/Xx5mjojiFqmSpJQX3uD\nvVIo+w4xHpkpiqGPVO2GMq6wWof10osjZIvTTlBHXLM3XbQMpcqUaW8wDJybpkqvm3EDTBRDXwlt\nZuX3t6UxZfM20m3VxZw6c3O79/UpBWZ6kSVVGsjFqbFj9DIaw6YAKH46xhcnikGwjyzLzxHQ2ATL\n8ahqcmLYkGsVkzL2Ofl9JvXZWXNCvKy6MW5Uu3OAK5V4zV15wBl3h0QxCL3NY1zFB8+0IdcaIULd\nVV2KysU1BTb7OAifg6pKVFLLiGIQejtEuKwPfjzeLfCK73SV95mpLHASGvpLvS5f1WMdOe0qPGuN\n5kcUg+Buttma3B2QaJSeuiuM3Ca8q4ZdshV0DkLaDiFmHp/CUYqxInU5TVjbwPyIYhg6IbOemPr6\nbN5SM2UEvG0eoxB7ug1O/sYQgRvSYxiPw8rgMuOxbmCXyeDI+oJ2EMUwdGwVeN++3hveQ1wrptG3\nZRqDvoZkU63j0J4L1Q8Tqhiy/zZ1WjvQ7mA7crkqohiGjktSUfJKdKRrbCLEPGKbSjlWjyETqK4I\nnbKYZFeZBm0+yKaOKsHWZKR1ucKxviA3ohiGTozhwR2t/KEOVZ8wLOtjoCxlG5u2hmuoUA8pd4c7\nkWY63vovQyOKAcB+AI8C+Hb6eaNlvwsAngbwRL5g1OOLiygGApRK72pa+8J3GGO6dNulUlrKIect\no5TKCFzbeWzXUxx3kOl8yvV37PHT6bi/oAxNKYZPATiRfj8B4Pcs+10AcKDs8cVFFAMRX5fX9mLM\nzHT+hSleui2T58qK3/QUOsA8VDGUMdG4zuPL3xbas+lYh5FOxyOMytCUYjgP4FD6/RCA85b9bIqB\ndHxxEcUQiVj2iI5g05NVW8xlzUhl9e3Kil+I58du5JMGhkYwd6gtEI70GGpTDD/MfVf534X9/j41\nI50DsBx6fLp9GcAmgM2FhYXabtzgiOXB7DBU84/t8qv4NLJGK7VV7lIKprEcVRRWL81HeUKyNnbU\n31YkmmIA8BiAZwzLnUVBDuAHlv+4Of18PYAnAbxXFxSD6/jiIj2GmhmYU44qQG0dpjLmo7KC2OUT\nKKYTD1FY+dQ/IYqq8/iEfs/eBVampMIx9wL4ZNnjtSiGZuhRK4kCRZCG9hgmk3IpJ1yEKJlB9wyo\n9besH66jveemFMOnC87jTxn22Qvghtz3rwI4Qj3etIhiEOrC1foP8TFkwja0N+Fz41CjiPKtf1PP\notc6P8REVDZyr6P+tqYUwxjA42m46WMA9qfr3wjgbPr9Lan56EkAzwJY9R3vW0QxNIwv/0SPpIwr\nTY4L222oElZqwuVjoPQcetk7KEJt5VcZ6yM9Bn6LKIYGcbWqyuRfYK5EYpuUXXmbyp6nOB3Avn12\n2dWBWx4faiufksekszk9zIhiEOJQxoBuak11yIkXW5i60u7EOE+Hbm0zxOgx2CISOjixVR5RDEIc\nXHaKEPtrz7rk3Bhkz0Br84XH8DH0tL6KYhDiYDOGj0b+l4cyHLijTjyBAT4zZ5WopJ45nTOoikEl\n+3aLxcVFvbm52XYxhoFS9m3r68DyMrC9vbNufh44dSr5XtxmYjIBLlyoXExhgBw+DGxt7V4fo07V\n+d8topQ6p7Ve9O0300RhhA4zmdjXLy0lSmAySRTIZJL8XloCVlf9SmF+Hlhbi19mYRhcvBi2PoS1\ntaR+5hlQfRXFILjxvSBLS0kL6vr15HNpKVnvejmLSkQQyrB/v3n9wkL1/3Y1egaAKAbBTdkXxPZy\nTia7lYgg+NjYSMw7MzPJ5/HjwCuvmPe9ciXZvyq2Rs8AEB+DUA8bG3b/w4BeMCECprqkVOIOtiF1\nzYj4GIR2GXhXXIiIyV/la9BubyfHCaUQxSDUR7ErDkybA2J094X+UDQXZfWjrDM5hhN6oOxpuwDC\nQCiaA7a2kt+A9CIEd/1YWDCHjvrMSTGc0ANFegxCfEwtP5M5QLr7Qoarftgi444d2wmnLo63GVBo\naR2IYhDikrX8traS1lzW8jO1+ADp7gsJrjEJNn/VffclJkqtgQcfFH9WRCQqSYiLbcToaARcu7Z7\nfcdHkgoWsl7ixYuJSWdtzS2oezrSmBsSlSS0g63ld+1a+ZGkNqekwBNbr9H13AY+0pgbohiEuLgG\nthXNAXffnbQqXQLfJGQ++lFg3z6aohCl0jxl/EkS3swLSqY9botkV2VMzGkVtaZNmmybfIByjsHm\nq66RnmYm7QOQtNtCa1CELTXfPXXSZFOefEpacJnhJj6hcxmIcm6MRhQDgP0AHkUyZ/OjAG407PPP\nADyRW14B8O/TbfcCeDG37SjlvKIYegC1VUnpMdhaoy6lMpmYp23swWQsjZEX6ONxsiil9d695vu6\nsmL+D1HOjUFVDFV9DCcAPK61vhXA4+nvoqnqvNb6bVrrtwF4B4BtAF/M7fJfsu1a67MVyyN0BZsv\norh+bc09J4Tr/1wDnLa2gEuX7NsEN0Xfz6VLyaI18OMfm485a3i9ZXwLS6oqhjsBnEm/nwHwq579\n3wfgO1prefOGDjUKZWkpGcjkwha9YjoHBaWacVJzdYxTykWZb6OIKWKtzjkVhPJQuhW2BcAPc99V\n/rdl/9MAPpH7fS+ALQBPpdt2maJMi5iSekKIbdlm9hmN3Me55u8N9VmEltl37UUTytzcjjmmLVs7\n1bRD9f1U9QMJUUEsHwOAxwA8Y1juLCoCAD9w/M8cgJcB3JRbdxOAEZKeyxqA047jlwFsAthcWFio\n+/4J3Kga6RKqHEz/6xOaIUqjSrRVnVAFdej9rBI5JkQjmmJwHgycB3Ao/X4IwHnHvncC+EvH9sMA\nnqGcV3oMA6Rqy9IkgJTSet8++v+6yhAq4KpEW8Umr9CoitJ0vWV7PxKV1BhNKYZPAziRfj8B4FOO\nfR8C8LHCukO5778F4CHKeUUxDJAYLcuVld3Cb3Y2EWKU/3X1WnxKoyj4qkRblcEmfH0C3qWgbFFJ\nItzZ0pRiGCOJRvp2anLan65/I4Czuf32ArgE4GcKxz8I4OnUx/BwXlG4FlEMA6Vqy9ImjMdj2v+6\nhL+rtW1SaCsrcQRysbymbS6lytWkJdRCI4qhrUUUg1CKqn6KMgJ2NLIL/GKLe3bWL5BdZbBtsznu\nbetN5RR6gSgGQSgSIwImxCRTFPQ+ZVR1xHiZ6CvpKQwKqmKQtNvCcDBNKl910vh8eun9+5N1ly8n\n3//hH4CrV83HlU0nPTOTiO0i2SDA2O+zpL3uFZJ2WxCKxM7gaRr9++qryaQx+/bZlUKVdNKZ8jGt\nt430Ho/LDfQDZKDZQJEegyCUxTW5zMWL9tb7+np5ZXTggDmVx3gMnDxp7xEB0xPnXLliTwmSR3oM\nvUJ6DIJQN650Dq55KarMMXD5sn29q0e0tJQI+OvXk8+TJ/29CJkoZ7CIYhCEsrgSAdY1I5kv+WBR\nAdiUkEmJrKzIRDkCAFEMglAel/Cva0aymAqnqETuu4+mVITeI4pBEMriE/7U1nvMcwpCBMT5LAiC\nMBDE+SwIgiCUQhSDIAiCMIUoBkEQBGEKUQyCIAjCFKIYBEEQhClEMQiCIAhTiGIQBEEQphDFIAiC\nIEzRyQFuSqmXABjSWno5AODlyMXpG3KPaMh9oiH3iUZT92mitT7o26mTiqEsSqlNyqi/ISP3iIbc\nJxpyn2hwu09iShIEQRCmEMUgCIIgTDE0xXCq7QJ0ALlHNOQ+0ZD7RIPVfRqUj0EQBEHwM7QegyAI\nguCht4pBKfVvlFLPKqWuK6Ws3n6l1BGl1Hml1HNKqRNNlpEDSqn9SqlHlVLfTj9vtOx3QSn1tFLq\nCaXUYCbD8NUPlfD76fanlFJvb6OcbUK4R3copX6U1p0nlFK/3UY520YpdVop9X2l1DOW7WzqUm8V\nA4BnAPw6gK/YdlBKjQB8FsAHANwG4CNKqduaKR4bTgB4XGt9K4DH0982fkFr/TZOYXV1QqwfHwBw\na7osA/jDRgvZMgHv0N+kdedtWuv/1Ggh+fA5AEcc29nUpd4qBq31N7XW5z273Q7gOa3181rrqwAe\nAnBn/aVjxZ0AzqTfzwD41RbLwg1K/bgTwOd1wtcA/KxS6lDTBW0ReYeIaK2/AuCyYxc2dam3ioHI\nzQC+m/v9QrpuSNyktf5e+v3/ArjJsp8G8JhS6pxSarmZorUOpX4MvQ5Rr//dqXnky0qptzZTtM7B\npi7taeOksVBKPQbgDYZNq1rrLzVdHq647lP+h9ZaK6VsYWrv0Vq/qJR6PYBHlVL/J20BCYKPbwBY\n0FpfUUodBfCnSMwlAlM6rRi01u+v+BcvAnhz7veb0nW9wnWflFL/Tyl1SGv9vbTb+n3Lf7yYfn5f\nKfVFJCaEvisGSv0YRB1y4L1+rfUrue9nlVL3KaUOaK0lh9I0bOrS0E1JXwdwq1LqFqXUHIAPA3i4\n5TI1zcMA7k6/3w1gV09LKbVXKXVD9h3ALyFx7vcdSv14GMBvpBEl7wTwo5xpbgh475FS6g1KKZV+\nvx2J3LnUeEn5w6YudbrH4EIp9WsA/gDAQQB/rpR6Qmv9y0qpNwJ4QGt9VGv9E6XUJwA8AmAE4LTW\n+tkWi90GvwvgT5RSH0eSsfZDAJC/T0j8Dl9M3+09AP5Ya/0XLZW3MWz1Qyl1LN1+P4CzAI4CeA7A\nNoCPtVXeNiDeow8CWFFK/QTAqwA+rAc4slYp9QUAdwA4oJR6AcDvAJgF+NUlGfksCIIgTDF0U5Ig\nCIJQQBSDIAiCMIUoBkEQBGEKUQyCIAjCFKIYBEEQhClEMQiCIAhTiGIQBEEQphDFIAiCIEzx/wF/\n91YS2JxvawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe07e4178d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Let's first define some funny manifold in R2\n",
    "N = 256\n",
    "#norm = np.log(np.linspace(np.exp(0.1),np.exp(1),N))\n",
    "norm = np.linspace(0.1,1,N)\n",
    "angle = 2*np.pi*np.linspace(0,1,N)*1.5 #we want a loop of 1.5 times the circle\n",
    "angle2 = 1*np.pi+2*np.pi*np.linspace(0,1,N)*1.5 #we want a loop of 1.5 times the circle\n",
    "\n",
    "#defines 2 sets of points, that we will try to cluster later, and add noise\n",
    "NbOfClass = 2\n",
    "noise = np.random.rand(2,N)*0.1\n",
    "points=norm*[np.cos(angle),np.sin(angle)]+noise\n",
    "noise = np.random.rand(2,N)*0.1\n",
    "points2=norm*[np.cos(angle2),np.sin(angle2)]+noise\n",
    "\n",
    "#Simple example for isotropic gaussian distribution\n",
    "#sigmaNoise = 0.3\n",
    "#points = np.array([np.random.normal(-1,sigmaNoise,N),np.random.normal(-1,sigmaNoise,N)])\n",
    "#points2 = np.array([np.random.normal(1,sigmaNoise,N),np.random.normal(1,sigmaNoise,N)])\n",
    "\n",
    "#Point cloud\n",
    "plt.scatter(points[0,:],points[1,:], color='b')\n",
    "plt.scatter(points2[0,:],points2[1,:], color='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we will concatenate those points\n",
    "points = np.concatenate( (points, points2), axis=1 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Threshold: Below this value we consider that vertices are not connected\n",
    "thresh= 1e-4\n",
    "#Sigma of the gaussian distance metric\n",
    "sigma = 5e-2\n",
    "sigmaSquare = sigma**2\n",
    "#will contain the value of the edges\n",
    "data=[]\n",
    "rowIdx=[]\n",
    "colIdx=[]\n",
    "\n",
    "#Now we are going to fill the graph matrix, neglecting small terms\n",
    "for idxI in range(points.shape[1]):\n",
    "    for idxJ in range(idxI+1,points.shape[1]):\n",
    "        dist = np.exp(-np.linalg.norm(points[:,idxI]-points[:,idxJ])**2/(2*sigmaSquare))\n",
    "        if( dist > thresh ):\n",
    "            data.append(dist)\n",
    "            data.append(dist)\n",
    "            rowIdx.append(idxI)\n",
    "            rowIdx.append(idxJ)\n",
    "            colIdx.append(idxJ)\n",
    "            colIdx.append(idxI)\n",
    "\n",
    "#Practical implementation using scipy support for sparse matrices\n",
    "A = csr_matrix((np.array(data), (np.array(rowIdx), np.array(colIdx))), dtype=np.float32)\n",
    "\n",
    "#Now the diagonal matrix\n",
    "D = dia_matrix((np.array(A.sum(axis=0)),0),shape=A.shape)\n",
    "\n",
    "#We can define the laplacian: L = D^-0.5 A D^-0.5\n",
    "L = D.power(-0.5)*A*D.power(-0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering algorithm initialization\n",
    "We will now implement the spectral clustering algorithm from Andrew Ng and Michael Jordan, which reads as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Spectral Dimensionality reduction that will be used, must be equal to nb of class\n",
    "SpectralDimensionality = NbOfClass\n",
    "\n",
    "# Get the two first eigenvalues of L, a symmetric positive definite matrix\n",
    "# See http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html and\n",
    "# http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.svds.html\n",
    "[w,X] = spl.eigsh(L, k=2)\n",
    "\n",
    "#with w : array Array of k eigenvalues\n",
    "#with X : array An array representing the k eigenvectors.\n",
    "#The column X[:, i] is the eigenvector corresponding to the eigenvalue w[i].\n",
    "\n",
    "#As stated in Andrew Ng paper, we will now perform row normalization of these eigen values\n",
    "Y=np.divide(X.T,np.linalg.norm(X,axis=1)).T\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run K-Means unsupervised clustering on spectral data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "NbIterKMean = 250\n",
    "iterIdx = 0\n",
    "\n",
    "#Initialize centroids\n",
    "centroids = Y[np.random.randint(points.shape[1],size=NbOfClass),:]\n",
    "#Initialize the previous minimal distance map, with cluster 0\n",
    "distmap = np.linalg.norm(Y-centroids[0,:],axis=1)\n",
    "#Initialize \"belonging vector\" that contains for each sample its class\n",
    "belong = np.zeros([points.shape[1]],dtype=np.int)\n",
    "#Do we need to change cluster ?\n",
    "changeCluster = np.zeros([points.shape[1]],dtype=np.bool)\n",
    "changeCluster.fill(True)\n",
    "\n",
    "while( changeCluster.any() and iterIdx < NbIterKMean):\n",
    "    #Assign points to clusters\n",
    "    for clusterIndex in range(NbOfClass):\n",
    "        #Compute the current distance between every point and the current centroid\n",
    "        distance = np.linalg.norm(Y-centroids[clusterIndex,:],axis=1)\n",
    "        #Does we need to change cluster (boolean vector)\n",
    "        changeCluster = distance < distmap\n",
    "        #if so update belonging matrix, and distmap\n",
    "        belong[changeCluster] = clusterIndex\n",
    "        distmap = np.minimum(distance, distmap)\n",
    "\n",
    "    #update centroids\n",
    "    for clusterIndex in range(NbOfClass):\n",
    "        centroids[clusterIndex,:] = np.mean(Y[belong==clusterIndex,:])\n",
    "        \n",
    "    iterIdx += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate classification results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fe05065add8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+MHdd1379nH3cRkXQT8ZGWadm7lAG1gFykhk2otmG4\nCpomNINASeAaNtaKKglYcBUDbgoXpUAgEQosmsj9o3RsWmUNJrJ2YyFA4VhI5CiSGsBtDSdeBfpZ\nl7Zkc2kJbiSRtgWbBvSDt3/MjHbe7P1x7vw8M3M+wOC9N2/ee3fm3Tm/7rnnkjEGiqIoipIx13UD\nFEVRFFmoYlAURVFmUMWgKIqizKCKQVEURZlBFYOiKIoygyoGRVEUZQZVDIqiKMoMqhgURVGUGVQx\nKIqiKDPs6roBZdi/f785dOhQ181QFEXpFY8++uhLxpgDoeN6qRgOHTqEzc3NrpuhKIrSK4hoi3Nc\nLaEkIjpDRC8Q0VOO94mIPkNEzxDRE0T07tx7R4jobPre8TraoyiKopSnrjGGPwFwxPP+hwBcm24r\nAD4PAEQ0AfC59P3rAHyMiK6rqU2KoihKCWpRDMaYrwG46DnkRgBfNAnfAPALRHQQwPUAnjHGfNcY\n8wqA+9JjFUVRlI5oKyvpagDfz71+Lt3n2q8oiqJ0RG/SVYlohYg2iWjzxRdf7Lo5iqIog6UtxfA8\ngLfnXr8t3efavwNjzGljzGFjzOEDB4LZVoqiKEpJ2lIM9wP47TQ76b0AfmyM+QGAbwK4loiuIaIF\nAB9Nj1UURVE6opZ5DET0JQA3ANhPRM8B+H0A8wBgjLkbwAMAjgJ4BsAlALek771GRJ8A8CCACYAz\nxpin62iToiiKUo5aFIMx5mOB9w2A33G89wASxaEow2RjAzhxAjh/HlhcBNbWgOXlrlulKE56M/is\nKCLZ2AAOHQLm5pLHjY2d76+sAFtbgDHJ48rKzuMURRCqGBSFS1EJ3H57WOifOAFcujT7PZcuJfvr\nbo8qG6UmVDEow4IrLGOF6sYGcMsts0rg85/3C/2NjeQ4G+fP88/J1R71RJSmMMb0bnvPe95jlBGx\nvm7M0pIxRMnj+rr7uN27jUlEZbLt3r3zeO5xeabT2eN9G5H9N/Lb0lK1a7K0ZP/eySR8nZTRAmDT\nMGQsJcf2i8OHDxutrjpg8oO1+/YBL78MvPrq9vu7dwOnT+8cwD10yG6hLy0B587FH5eHiN/+paXk\n0eUtuNofw9xcogp81PE7yqAgokeNMYdDx2koSZFFMURy4cKsUgDcMXpXeGZrazZs1FR4B0iE8dqa\n/7vqENaLi+FjbNdJxyUUBqoYFFnYBmtt5AVvJuxcFjTRbCzeZf37hO10at+/d2/iIRAlj5nQd33X\n0lI9FvzaWqKEQhSvk45LKAxUMShxNG1xcq32TPDmhZ0Nop0Kw5idyiGz9F2cPAksLMzuW1gA7r47\nCT9dvpw8ZkLfJrhDvxHD8nKihDKlNJnYj8srqCYzpJRhwRmIkLbp4HNH2AZU5+eTgdm6Bjxdg6qu\ngWLf8aHvim03dxCcc3zsd3HaVvxvFhZmz9E3WK6MAjAHnzsX8mU2VQwdESu0Y1lft2f/FAVc/vuJ\n/MLO1+Yqba1Cmawo7vdmymY6TZR28Zo0kSGl9AZVDEr9uASLS9DEWMWu9M7p1P85l+DPt6HJtNEy\nhNpsI9bDcP1G8T/sSjkqnaCKQamGTRBxPIZM+MRaxWXz8jm/s77ub2vbhLycImU8DJ8SrzOEpfQK\nVQxKeVyCaHXVb33nBU+sVczxRlzCkGNNl7HSmyK2LWXaLul8FTGoYlDiCXkFec8hi2MvLNiFt88q\nruKNlBVsTcX122hLrIeR/Ubxv1lYUA9h5KhiUOIIxeJdgshlrbsE/XRazRupEvqpOxOoCjFtKTsm\nURx8np/3lxORcm2UxlDFoMTBsdhjrHWXVeyqOVT0RiaT6m0YCmW8nRhlIsmbUhpFFYMSByfGHyso\nbFYoNyyiwmqWWIs+Jvyk4xGjgasYtIiekuCrIQQkJSFeeqm537EVsNOVz8oTc51dBfmIkhndymBo\ntYgeER0horNE9AwRHbe8/++J6LF0e4qIXieifel754joyfQ9lfZd4au9s3t3UhKiqd9xlYpYXraX\nm1DCxFxnV12n/P6m1rlQZMJxK3wbgAmAZwG8A8ACgMcBXOc5/tcB/I/c63MA9sf8poaSIoidZJaF\nFbIYfxMDkTrQ2Q7c67y6ag8lra5uf0+o3Mb6uob/egDaCiUR0fsA3GmM+dX09R2pwvlPjuP/FMDf\nGGP+W/r6HIDDxhh2nEJDSUyyAnP5wmlao18pEgo7hcKMQNKvrrgiKZNeZDIB7rlH+5wA2gwlXQ3g\n+7nXz6X7bI3aDeAIgP+e220APExEjxLRSg3tUTK0mqbCwVXRNtvPqXh76ZJdKQDA669ree+e0XbZ\n7V8H8L+NMRdz+z5gjHkXgA8B+B0i+qDtg0S0QkSbRLT54osvttHW/hO64RUFcI8xzM0lwpyzKFAI\nNUh6RR2K4XkAb8+9flu6z8ZHAXwpv8MY83z6+AKALwO43vZBY8xpY8xhY8zhAwcOVG704LAN+nEG\nFRXFlXiQWfpHj/IWBZpO/cepQdIb6lAM3wRwLRFdQ0QLSIT//cWDiOjnAfwLAF/J7dtDRG/KngP4\nFQBP1dCmcWFbmevWW+3ppXUuFqMMg2zRH9tiP5cuAQ88MLso0HQKzM/PHpdlrrm+B1CDpEdUVgzG\nmNcAfALAgwC+BeDPjDFPE9ExIjqWO/Q3Afy1MeanuX1XAfhfRPQ4gL8D8JfGmL+q2qbRYRtLeOUV\n4Kc/nd03nerAs2Jnedk9Z+H8+dnU4ZdeAv74j+1Lmi4vJ0ZK7Ap5iih0gtsQcE1QKmKb3KQoGTGT\n4lzYMuGIgGPHgFOn6milUoFWJ7gpHcN10TXGq/ioY51qm/dqTBKO4qAT5ESgimEI+GYt59EYr+Ij\nG2uwhYi4xGTCFZXA7bfvHCvTNNdOUMUwBLIbejp1H6Mx3sqMwpitWoaEmwlnS5i4+26ddyMEVQwS\ncEmcWEn0s5/Nvs4GAMtYfsoMNjmWGbOjUBhcuOEoV8jJRmjWtVI/nLoZ0rZB1UqKWUbTVneGs+qa\nUhnXJZ6b27lv9OWBVle3a21NJts1l/JwyrznS4WP+oLWB3Q9hp7gkjichWrKrrrWY7qov7e+zpdh\nxb9pdPUCuYX0uEu5qoFTK6oY+kKM5VQU9HWvuiac1dWdl8tlndclkDm612fkjq7YKHfRH9fFGYmB\n0xWqGPpCFY8hpFQGJIV8i79xZU4ZBRJr2ObbNMqF0WJWjrNd/FFetPZQxdAXqowxhMYWBqIUjPGf\nalHmVDFaiWZD4rEOXf5vCsnIQYaZqgr2GLdQiUYVQ59wSYiQ5BhgrMJ1yj4BXZQ5XKPVJcPyY52x\nHsN0Gv5sdl6hv66XiqNKn+RoaqUSqhjGQi+lxzb55k+nxszP22UKR4hncI1WjrLhju/H6u1QG3ut\n88v2SQ0jNY4qBkU83IFdl4XtMibrSozJO26cIR/XOcZ4QJlXM0oZGTM+oZSCqxh0gpvSGbY5Tjay\nagpXXLG9bzoF7r3XXpfNNhE8/9mMtbWdRUDzZBPYlpeTlSnLlBEqTiQGkklwxtiPzyYIj3KNJV0/\nRAyqGJTO4Aq5ffsSIZ1fObI4ydvGyy9vP79wAbjlltlZycvLSdFPl3LIV2OoUkYomxlNBNx0k38i\n79GjyeMoZWQdRfyUeuC4FdI2DSUNA87A7u7dydhDbFjF9ZnpdOexvglsVaMYsfMgBjHGUIWej5lJ\nBxpKEka+oM7+/ck28uI6rqKwe/bMWuUXL+48BvB7HK516W37l5eT37JR1ULnhssysnNyeSjAwOsy\nuYr4aUGqduFoD2lb7zyGkNnYc1OwipHHSVuPsf4zfFa56xyasNBj50H4vKC62tg7o9w316dXJ9I9\n0KwkQQy4dEVVYcXJvimjGMp8pgmBGTMPYmHB/5uxmUq28+lliMqXq9yrE+keVQyS4JiNPU3Jc92z\n0ylPyHIyFMtkMa6vJ4I2RvA2gSvN1nY+8/P+9sVWm7ApgDLjNZ0T43aJPpHuaVUxADgC4CyAZwAc\nt7x/A4AfA3gs3X6P+1nb1jvFwDEbJ5NeusTce3Z+PhFKxVPkWMFlc/qlhEzqKgnkK6vFva69tEti\nTkb0iXRPa4oBwATAswDeAWABwOMArisccwOAvyjzWdvWO8UQm5rSI5e4bJG57BS5pSF6F/4IUNYL\n4szCLvN/iDa0fYWnenUi3dOmYngfgAdzr+8AcEfhGJdiCH7WtvVOMRizs/ZDZj6XnVIrhLJlqfOn\nyLHspVj/dVGHF+TqOr5tOu2pko05QfEn0x1tKoYPA/hC7vVNAD5bOOYGABcBPAHgqwDeyf1s7r0V\nAJsANhcXF5u8du0ygDIARaHtimP3+BRrpw4vKNY7yHtpvVOyLk2qS+hFwVUMbc1j+HsAi8aYXwTw\nRwD+PPYLjDGnjTGHjTGHDxw4UHsDK1M2z3qAU1w/8hH7/IQiPT7FylSZSZ0Re/1uvjn5fttUAfHT\nBFyzoq+8cuex+SnrSjk42sO3oUQ4CMA5APvLfNYYgaGkuksNt2jxVLUeOSnm0+nODCE16qpTZVZ1\nZoBn4aheZH7aOqvPbVJ2gBZDSbsAfBfANdgeQH5n4Zi3AKD0+fUAzgMgzmdtmzjFULUUZke+PVcn\n+ZaL4A6R9DJ80QPyQj4UWspKlBeVdEiRiCamHrvSnmJIfgtHAXwbSYbRiXTfMQDH0uefAPB0Kvi/\nAeD9vs+GNnGKoadLdXH0mc8j0CV6ZbG+Hh7fsYXke/sfhk64F5qtXVpVDG1v4hRD1aW6OoIz7u3L\nne+9tTlQbGVGymzFWeKi7Btu7q4yA1cxaBG9OvCVC7ZVURMyOMYZ93aViH79dff3aqXkbjl1Klmr\nwlUYsAwbG0np862tROpubW2vV9EJnOqEY85uqApHe0jbxHkMxthH9HyzvwRYMxxnJjZXfjIR4Qwp\nKVU8B47n2JlnGDqxXufmNgfUY2iZ5eVtzyEzp7e23KvACLBmOCmTIc+g+Pqee+JSLpVmqdLN8p8V\nt6Kc78TyNcqLbs6tt2rJewaqGOrE5t4as1M5CIq1uJaezO6b/PKYebJ7r0oevtI8tijnwkLyv2b/\n2+pqeOE0cdNt1taSEykyP5+8t7xsvx9feSVZlENEPEwwHLdC2iYylGSM370V6M4WvWxbptH8/M5w\nUqgKqCKLOkqO+KrEdtalQ6Vitb7SDqBZSS3BKV4jsOPFlIMubl2Ur1a6xzdnopNEu1BaHbfCo4Dx\nvrbgKgYNJVWhmKphC8gTJe8Li2e6ol4cXnlFRFKVUhO2chi2fVnYcWlpZ1+5dCkpudFqFw/Ft1xr\nx3K/Z8xwtIe0TYzHEEryF2FW2ama5z4iI8sYM9zkFpvnuLCQhAtdXZeTDNRZ42012/O1WXwnNgKg\noaQW8Lmy4vL7ZonxsgWfRj0EpL7gOYqViVlPg1OeqPW+Eauxh6rhmahiKENsp/ENfgkvp82ZODqZ\n2AekhyIUjTEsqR/S8X2WNTGeo8/eEdjFw/T5jyuJKoZYYs1CVzWyLGVHgMfAyTThzMEb7P3DrALo\nuz4+BZtfj0nqdYvxGLjrbPTCmxyyG+hBFUMssYLcdXxWYEZAOW3uzwvQYe3DcJlCNdoywc8VrBLl\nDneMwbavL+doZZSdXhVDPLGhH87xHZraoRTvPKMwniKXmbsMMv9mft0rBHft4iuFGLnTdrex/V7M\n5ep8PkMZhId6m0IVQyx1eQwCLI71df9N7PqM2HARR3L56uLErmiTbt/DUrTg52y+aytVSfvGIrpu\nWykE379NooohFpfwcC0uLvUONv64ce/6fZVYR/Z/xATSc9vroEYUg2/muFR5JbVdO3AVs+RM5RZy\n/zaJKoYyuILKrg4j1MwelHVXUqjPSK6SkzZ+jD2NKIbM3rAhNcLRCznq8wxt1Vb7kB1QM6oYytIb\n08hNaFy8V9Sx4oxv83z/a5jzfnRhIb4seX6zIbn7CbWDtgkZEdNpD7Rbs3AVg5bEKCKuvnA8tkoA\nRElRSWGVOcI0Xa7AGOdbc7g8Uz12dXW2muxttyUlI/LMz7sr0hax/Re+NZ+6pliJV1wl3dA9euGC\n2EWzxMHRHqENwBEAZwE8A+C45f1lAE8AeBLA1wH8s9x759L9j4GpzdRjCCOu4FlZSg4c17JNJt6m\nubpKzLrKtv9CvGUulbJhx67jdC3ClbF1KIUJgGcBvAPAAoDHAVxXOOb9AK5Mn38IwN/m3jsHYH/M\nbzaqGHoRTOXTOz3nyjayjf3UFWbas8e+f3XV29S6fn4y6a8SKP41rlyN1hrjG2OIyeEeKG0qhvcB\neDD3+g4Ad3iOvxLA87nXchQDN6NBOHnZ2isjyaeUqw5Ch0z21dXt/zyrBRKA26S9e/lKRKIN4ssC\ntk1667Qsu+8eHpjRV4Y2FcOHAXwh9/omAJ/1HP+pwvHfS8NIjwJY4fymiJIYQuFGXkQaST73pgml\nUHGB6qaiXJkcc4WT2gw1ldXV06nQcNjI43QiFQOAXwLwLQDT3L6r08c3p2GoDzo+uwJgE8Dm4uJi\n/VesdzEXOxwZKlbfhVycOqVvTRchL2eqZCjZmmdrbtv2i++2iPlLOu1zI1cGecSFkgD8YjoW8Y89\n33UngE+FflNESQyhhG7WrPCbSFxSKFbi7trlL+zTUCDcJbRj6in5TtdXmqIp+6VMZXnX1omNFTtx\ndeC0qRh2AfgugGtyg8/vLByzmGYsvb+wfw+AN+Wefx3AkdBviiiJIRTOzSoWl2SNlT75iUx5SduC\ntegaO28jsaoJfLfF+nqczu7ExvLdEGJd5+ZoTTEkv4WjAL6degQn0n3HABxLn38BwA/TsYQ30lLT\nTKbH0+3p7LOhrbExhqKV2ZNV74uTOUNVMMXiSnHhmqaB9NIuaXL8vMlT94WuXIPPPp3dOiIb1R2t\nKoa2t8YUQ3F9hR6seu8qJeS6D8TOfg5JIK7JLZz19eYmczfZZluIPkbRdWach1yanoWKq6KKIZae\nhpJ85S96pedC158zyiv8v8pYXa1fKXRx6lwF1+l4r8QL1yFcxaAlMTJ6WgrD1byLF4EzZ2ZLOJw5\nI7CMQYbrRLa2kroTJ04kdSEuXwbuuUdu3QgGp07xy2ZwmJ/v5tQ51UqIOi6fsbTkfq9HfaZ1ONpD\n2jZmj6Ho1g9mMmdsnm3PUxDrGpDuMrmGcw6d90PNSpoBGkqKxNaBBC1N5aoKMT+/M2TUu2QL18mJ\nlDT14av08XM/149L0YuaXD03IupEFUMZhPbykGUmdpapD9+19m0DGyz0jRH5lg7tWui60nJ71w9H\nhiqGKggLK4WiLL2TlZwYRM8HmLn4dKIrFbRroWv7++bnB7DmzQg0G1cx6OCzDWED0aGfbXrJgto5\ncWJnXfwir7/e6wFmLvv2ud979dX22sFhYyNZQ+LjH9/59736arLcgTFJvsCttwL79yd5AyLWAMka\n72rQxgawspI0PjuJlRUBDe8IjvaQtqnHMGup9c6w4YSOiJK8zoFbcLHlMroKJVUdLO+8VlJo8qqw\ne74poKGkCgirtLq+7r7hRE5YC5UG5dZRGNhNaaOsoG378tQxa7uzv9OlffM3z0BqpYXgKgYNJdlY\nXgZOn56dBHD6tMhJABcvdt2CAj6XPHvv9dd53yV8DkkdFJcGjaHNy1PHb3X2d164EN7viscKiNOG\nomCNwNEe0rbKHkPPBpl81po4o9rnkveiHGd7+DxBzlZxOYmodtZRUryzv9PXqPxJCooSNNUsaCjJ\ngdAO4MMXkhfXbJ9LHiNFhP8ndVBHeKbpy1R2bKEY0u/07+SEkowRaTDWPfShisEF50oL6yC+XHdx\n1LGmwkhmpdZVTK9JS9z3d7qWyg6tQNc6PS2QaUz9Qx+qGFyErrRAj6JXeeNV0leYay23RdPCra4y\n3E2Oj8bOsxApb23l3HuS8aYeQ1uKIXSlhaatFddcEF0Go8p6lx2dSFEJrK42b2TWVYK7C49hbs6+\nX5wX67KqRN9A2+gYQ1uKwXWlMwuiC7MsEqG6aydlR1dbPJGYMk1NCL6qSqGLMQbfeh+CbpOEGLdM\n3A2UUKfnqorBh808FF8mcptepFxXCSm1dCJlm1gnZSe4ZV2yrawkTkVfYbdJwojrcNlQxRBDyKro\nwM206a5erFNTxxqWLZ1I2WbW2RX27vV3u67Lqhc9KlcIqYlrUwsD8BjqpFXFAOAIgLMAngFw3PI+\nAfhM+v4TAN7N/axtq10x+KyKDgamyliyIkKkdSwy0OKJlI3v1yk/QkK2y1wIWzKPbxM3vmAMv0+K\nuIGapzXFAGAC4FkA7wCwAOBxANcVjjkK4KupgngvgL/lfta2teYxCKyNZNuyskKdU4en0OLNWaW5\ndcH5ja5SP2Ouj2i5GvJi25opKIA2FcP7ADyYe30HgDsKx/xXAB/LvT4L4CDns7atdsUgLEW1jCUr\nwgsua4J3mIlUxsEhqq+53LlXXdCLNZ1jEHaf+2jKGOAqhjpqJV0N4Pu518+l+zjHcD7bPMJqI5Up\nzyKirFDZujItXetizZnsp33LAtswJqkcXgcnTyZrNueZn0/2dw3n71xaSiqhnzghqMS2i/x9DgCT\nSVI//MQJUY0WUQGcoz18G4APA/hC7vVNAD5bOOYvAHwg9/oRAIc5n829twJgE8Dm4uJiPepTKGUs\nWREeg88i881naMFi8zUtNmU18xrqbJvEuVahMYaFBXtCn1AjfBvBM0Z9t0od9zg0lNRvYpJ7RN2I\nLim3utqpVnNdyz17yoWTRCjiFvApzUyW9u76cG6sjjIRfX2xDmOkTcWwC8B3AVyD7QHkdxaO+TXM\nDj7/Hfeztm0MiiHDF+ctyl6plucbtGF+O6irLlGLTk5rcPpNmYQIEdhOLmbwpEU417hXE9yQZB19\nG0mG0Yl03zEAx9LnBOBz6ftPAjjs+2xoG5NiCJXczisF8S59h9lfddUlyrahEAqxZTI19vqI8Bhc\nJ8eNG7as3bjXWUtiqGIIupdZJxGWcWunQ+1VtjqH+GsaCXcWc9kQmyiDxFeWmHNyAj2Gqk1TxTAg\nQuMNPqtOjEuf0WG8q0r5CXFCrwR1zD90bZOJwBCm76YQWJUy5v9puuy2Lu3ZA5aXgXPnkkxaG+fP\nd7wyIWftweyYm25KXt97b3JSLaYEnzwJ7N5d7TtuvlnkCq8sTpxIsjOb4PLlZGv5L/XDvSn27gVu\nu63zdHVb1vx0aj+28fuaoz2kbWPzGDJ84aLOojScH3aZQh0syFMmNbXD6EKt1DEAL7pOVxFfJWUB\nA3I255lT31PHGFQxGGNmQ0nFmzvfSTqJ0nAGN3xxsA5nPZdREOJCcxHErq1g+6uEyFQ+tptCwIBc\nzDIRda4ppIphINg6UKYcRMRzOSviCTbDyygIEde9BDEGtO+cxadFhxAwINdV0VdVDANBgHHjJza+\nJdQMj1UOoq1kD3mDOQsLZSEL3/mKKNLIJaS5BNxUXS0ToYphIAgwbvz4xhi4ZpGA6pZl0lkzOdI3\nC9pW6mJhwa8cJ5OuW80kNOblchFb1vTqMahiqIQA4yaMSzLGmEUCTHDfojku5dyLiYUFfHMXfOcr\nGk5Ot7AkiC6WolbFMBD6KHjeIHa6cYe1ptfXkxszprm+uoBtKG6Op2I7xndOroFo0R4DJ2SZXQBB\nVhYnK6nu+1wVw0Aoer4u40ZkOKPMjKqOMpR8Qr7s1tT/wI2G2C5/6DxdXoPoMQZukSHxcdnmUcUw\nALjegmivImasoQPrrYzu6jJCFmpv3umqoz7U3JxwpWBM+A/pVd2YZlHFMAC4/bgX/Z0rpVq23mKF\n5+pq/ESx2LF12ySn7DXHsykzxCNET5djUJUmm0UVwwDger698JCFSqIY4Zk1rczEOFvFUluoqY56\nRpkiqsNjENWHXHBLxBa1rJiYa3uoYhgAg/EY1tf5ErjlDJEYRyZrVtmSGrainkWDta7y4NyJa1xl\nKB7XSO7IPYQiqhgGwCDGGIzxSzvbaGeLjY/JRsrkTd2L/uSFb53fna8AQRROR+3wb2gG8RaTm6aS\nSVQxDARuB8kfN52KWL52G5+0E3DzxngAZQQs5zsz6vIYit+bsbrqH6cQWT47ROw8GqHxMW5NtCqo\nYhgpIr0Hn/AXcPOWsdJj5zz4trwOLDOfgvO9RUT2kzKUmXnfQeYbZ74Jp15VVVQxjIh8x3NNUOrE\ne+aYQAJu3jJW+nQ6e7OXHXeYn98pKEKlKfKPro0j5MXOfbE1yrU/tlZXy9qP2wROH6zDVmpFMQDY\nB+AhAN9JH6+0HPN2AH8D4P8AeBrAJ3Pv3QngeQCPpdtRzu+qYthGbJ06W8NsZWEFlClw3bw+AV28\nnmVqLWWnWaSOcQYRQj4W1x/hq/XNqe7bofbj2j2c/7w3HgOAuwAcT58fB/CHlmMOAnh3+vxNAL4N\n4DqzrRg+Ffu7qhi24Vq7rXsMMZ6AgMJmrqQW1w1rO42yXkOROsYZeonrxH2rAwnwOH1wIqWcmfe9\nGmMAcBbAQbOtAM4yPvMVAP/KqGKoTIyV2roFGTt2IPQGt01oc92kZTKWbDWIqs5lECIT44m9eD2o\nYhjq1r7/uol1V9pSDD/KPaf8a8fxhwCcB/CPzLZi2ALwBIAztlCUbVPFECc85uY6aCB3NmqGgEFo\nFzHRiDIzo12/WUYpCJKJ8cS6SnnpKm6wJCGkt3xOUhOnUZtiAPAwgKcs241FRQDgh57v2QvgUQC/\nldt3FYAJgDkAawDOeD6/AmATwObi4mL9V6xnxN5DrRPSXNyZXcVkfGE3vg3bZNsyayWXkZPCL42f\n0JhBTzWgr/u2bQ+JCiUBmAfwIIB/5/muQwCe4vyuegzlSjm0TqguQ17ou27+3i0ybKdMxCPGKxTg\nWNWD7yTzKVniK/vxaDuC2pZi+HRh8PkuyzEE4IsA/ovlvYO5578L4D7O76picHcokYaVT4sVk/aL\ngVWhYw9Kgc+mAAAQGElEQVRlKOP4FD/jGuDu4eWw06uOXZ22h0jaUgxTAI+k6aoPA9iX7n8rgAfS\n5x8AYNJxhJm0VAD3Angyfe/+vKLwbaoY/Jl94qIuZWPHxogee+gC4WOt1Vlf37mMmThXuF7ajJS2\nohi62lQxJPQm9B6bZsOpETEQoVCG3vzveWJqu3Cnfo/UOKiCKgZFFjFpNsUaEaHF3XsnJUdGzCTG\nGO9yxMZBWbiKYQ6K0gbLy8DSUvi43buBtbXZfVdcsf18OgVOn06+b2MDuPVWYGsrERVbW8DNNwP7\n9wNzc8ChQ8kxSrecOAFcurRz/4ULwMoKcPvtyX81N5f8hxxs/aRHbGxsn7LIbsrRHtI29Rh4iDOm\nbZbjwoK7FGzIW+BMNR5UAL6nhFLoytQA6fF/2uU4ETSUNG7EDlLGaKvQ+IKGHPpBnbXEB/B/djls\nxlUMGkrqKSFX1Oa9X7qU7O+U5WXg3Dng8uXkcXl55zHZybnCCufPx/3m1pZwv33grK0loZ8yEM2+\n7nkICXB339hu3Sgc7SFtG7vHwPEGfFme4kJMeWIK03NXzBloDnyvcBVKDFUpFNxZyzatDx5D50K+\nzDZ2xcDpWK5jOOsOd0oo7BA7xhBTHlWpTkha2uqFiO6QdjjGmW9pCR1jUMVQO9xSvjFrDIiRk6Fl\nQDluUXYxfEpGc+Dj4JjHZSWeYK/ARZmqqRKyrFUxDBiuK2rrfOInEsf42ZxjXZpwbq5XgqhTuAK/\nTIxEuFJwNS90H0mdl6mKYcBUcUWldtg3iDk5zmIJmtJaHW6nibU6xKbOJfiaF7okUg0wVQwDJ9bQ\nChUxFXIvJpQNWxDtrLrJzZEXoxkFwpVyvsUFbP+lcCvFF3Ytu85C16emikF5g1CiT4vLK9cH987j\n5tB3bcpJwKaQ19f5i0lwa2JlHU6qWW38FVyy5uXX2ShWApfqDKliUN6AIxsldNoouEKFK6y6NuW6\nxnad5ufdlU5965tmysW3kLHwTAjOMiJls5K6RBWD8gaDjKbE+Or5O3Q63Sns6ryjJUoDDjGzk7nr\nToY6nuDcaV/TOWMMUlHFoBhj/JEAl8HdC3lWNZHcJbyrxACkxg84xNQr4oZ6QspG8GxL3zwgY0RH\nwbyoYlDYURSXpy/kHnVTRsC7VjLymYFcU1CSGRkrcJsodx3qgILN674OLodQxaA4O+/cHG+RLOmd\n3Ivr5LnrSpcxBUNmZFvWcRnPhWtFZFZ+TFtsYwk9sDyaciy7RBWDwqqXVFUWiiUmNMKJtVXxGLLv\nL7apuL51WYoSrOygbj7NxrfV0UZBErRs0wSfkhNVDEqlmkp5edW3zm+MiQuN1GUldxG7i/lNn6Yf\nafZWXy3/srSiGADsA/AQgO+kj1c6jjsH4EkAj+Ubxv18cVPFwIPT6X2GdU/rmyW4JsCV9RhifpeT\nrlmX0C07NsD1Mnr55/Pp61hBWdpSDHcBOJ4+Pw7gDx3HnQOwv+zni5sqBj4hi9+XfdH7m6Z48i5N\nxx1jCLlMxd+L9RrKxO64IbNiBbcyns3AlIIx/c0uKktbiuEsgIPp84MAzjqOcykG1ueLmyqG+vB5\nFYO8aXwprLECtvi9XA+lDo/BNahb1OzZ3I388qkc76C31kAcvTd+ImlLMfwo95zyrwvHfS8NIz0K\nYCX28+n7KwA2AWwuLi42duHGiEtWju2mYVv5tgvAzYLybdzaJOvr4bSy7LuqjHuElOEA4I4x9Has\nrUBtigHAwwCesmw3FgU5gB86vuPq9PHNAB4H8EFTUAy+zxc39RjaYWwDc1Epm0V8CiDGQudcYK4C\n85Wd8HkZwPb4SJ+lIJOQ0B/SfcBVDME1n40xv2yM+aeW7SsA/oGIDgJA+viC4zueTx9fAPBlANen\nb7E+r3TD8jJw+jSwtJQsvbu0lLy2LdM8CPIn7GNxkbcPSL5r715+GzgLc3MXB750CbhwgXfs7t3A\nyZPJOtzGAK+9ljy61uXuCaG10TmIXT+9STjaw7UB+DRmB4/vshyzB8Cbcs+/DuAI9/O2TT0GpXFC\nxXKKlBmsifFI8tSRipsfgxioV8CtmlJl/fS+gZbGGKYAHkGSbvowgH3p/rcCeCB9/g4k4aPHATwN\n4ETo86FNFUP7hGaBDk7GhIrl2IgdrOGWs7b9DmfqetbeocRBIqkyjyf2mL7QimLoalPF0C4+qyo2\n/tobJVJnYNlXt6lKwb78+MGePe6qsb256PXCsfQ5VUx6WtHDiioGpTZ8FlOMNdW7Qbw6BaovTbbp\n3xgpVb0BVy5CLxe2SlHFoNSGL1IRE38dkkuuyMGnc6uMMQyxv3IVQzArSVEmE/d+VzJOtj+fFbK1\nZT+Wm2SjKEU2NoCVlaRvGZM8rqwk+zlZdb5jXP1yDP2VEiXSLw4fPmw2Nze7bsZoIHK/t76e3Ij5\ndL7du5ObC9j5no2lpSQrUlFiOXTIbnDU0aea/O6uIKJHjTGHQ8epx6AEcaX1Ly35LS5b/neR3buB\ntbX626yMgyat+rW1pH/mGUt/VcWgBAndIMvLiQV1+fLsfCjfzTmKCXNK4+zbZ9/vCnHGMLoJnjlU\nMShByt4gvsnARSWiKCGKs5hvvx14+WX7sT/5SblZzkVcRs/Q0TEGpTGygUHb+MNYbjClHmx9iSgZ\ncHahfW0nOsagdM6YXXGlXmzjVSGbdvD1jBpEFYPSKEVXHKhe1EwZLq6id2UHk8eQWtoEu7pugDIe\niuGALOccUC9C8fePxUV76mgonFTHIPQYUY9BaYyi9ffJT46wfLHCxlfe2pUZd+zYdjp1cb7NWFJL\nm0AVg9IIthmprqUB1N1XAP+cBNd41alT20tI3HuvjmfVhWYlKY3gmjVqo88zSRU/GxuJxX/+fBLW\nWVtzC+shzjSWhmYlKZ3C9QJi3P06VuNS2sNXx8jGmGcaS0MVg9IIrkG/6XTW3b/55sSiDAn7jQ3g\nlltmhczHP558T0hJqELphtglMTW9WRCcEqzSNi27LZ+6llXMCK1p7/octx26jEH9+Eq1K90AXY9B\n6ZqQwI2pd+9TCr7PhX6jd4sH9QTfUteu9QxUQTdPK4oBwD4ADyFZs/khAFdajvknAB7LbS8D+Lfp\ne3cCeD733lHO76piGAYxi/xwFIPtc67fyASUyxPp82IsbZIX5tNpshG5l7MmKu/ZKdXhKoZKWUlE\ndBeAi8aYPyCi46li+A+e4yepIvjnxpgtIroTwE+MMf855nc1K2kYxGSh7N/vTnf1fS4mO6pIhVtj\nFNjqF3GwXVfNSGqHtrKSbgRwT/r8HgC/ETj+XwJ41hhT8lZVhkRMFsrJk+6V5Hyfs/0GB6L2Bqnz\ng+P79ydb1wPlnAF7znobRVxre4x5tTSRcNwK1wbgR7nnlH/tOP4MgE/kXt8JYAvAE+l7O0JRtk1D\nScMhJq7sCvtMJv7P+dbvjR2ziG1zCNeC812GU7hhHV+YLvZchri+skRQ1xgDgIcBPGXZbiwqAgA/\n9HzPAoCXAFyV23cVgAkSz2UNwBnP51cAbALYXFxcbPr6KQKJGZOwEascbN9bd5YTp01tC0eukOa0\nfTLhXQcdY2iH2hSD98PAWQAH0+cHAZz1HHsjgL/2vH8IwFOc31WPYZxUtSpd1vncHP97685y4ljd\nbaR35pUZtx11ezualdQ8bSmGTwM4nj4/DuAuz7H3AbilsO9g7vnvAriP87uqGMZJHVbl+np4ToTv\ne0Nei0txTKd2odemx+ASvCEB72uHKytJBbtM2lIMUwCPIElXfRjAvnT/WwE8kDtuD4ALAH6+8Pl7\nATyZjjHcn1cUvk0Vw3ipw6p0CWNO2CPkMXDj7pniKWt1+66D7T2fUuUoJw3rDINWFENXmyoGpQpV\nxipCXkvMOEY+/BRjdfva4HrP5SVxvCe1/ocDVzFodVVldFTNmfdVDI3J7SdKVraLxdd+oPy8DRu6\nbvKw0OqqiuKgahXP4nKly8vbef833QRccUVSLBDYuXhMnrKri/ly/uvO+9eFlMaJKgZldNRdxbNY\nXvrCBeBnP0uUg8shr1JOet8+935fVdsyE/0AnWQ2RlQxKKPEZvWXxVVe2lfCo6nwjMsbOnlypzLM\nvJoQum7y+FDFoCgVibWol5aqKYWLF937fd5QURmePBn2InShnHGiikFRKhITvqlD0Lp+L9vP9YZs\nSmR1VRfKUVQxKEplYsI3dQjaOpfALCqRU6fqC7Ep/WVX1w1QlL6TCU9XCmvdwjX0e4pSFZ3HoCiK\nMhJ0HoOiKIpSClUMiqIoygyqGBRFUZQZVDEoiqIoM6hiUBRFUWZQxaAoiqLMoIpBURRFmaGX8xiI\n6EUAZavO7wfwUo3NGSp6ncLoNeKh14lHG9dpyRhzIHRQLxVDFYhokzPBY+zodQqj14iHXicekq6T\nhpIURVGUGVQxKIqiKDOMUTGc7roBPUGvUxi9Rjz0OvEQc51GN8agKIqi+Bmjx6AoiqJ4GLxiIKJ/\nTURPE9FlInKO+BPRESI6S0TPENHxNtvYNUS0j4geIqLvpI9XOo47R0RPEtFjRDSauuehvkEJn0nf\nf4KI3t1FO7uGcZ1uIKIfp/3nMSL6vS7a2SVEdIaIXiCipxzvi+hLg1cMAJ4C8FsAvuY6gIgmAD4H\n4EMArgPwMSK6rp3mieA4gEeMMdcCeCR97eKXjDHvkpJW1zTMvvEhANem2wqAz7faSAFE3EP/M+0/\n7zLG/MdWGymDPwFwxPO+iL40eMVgjPmWMeZs4LDrATxjjPmuMeYVAPcBuLH51onhRgD3pM/vAfAb\nHbZFGpy+cSOAL5qEbwD4BSI62HZDO2bs9xALY8zXAFz0HCKiLw1eMTC5GsD3c6+fS/eNhauMMT9I\nn/8/AFc5jjMAHiaiR4lopZ2mdQ6nb4y9/wD8a/D+NETyVSJ6ZztN6xUi+tIg1nwmoocBvMXy1glj\nzFfabo9EfNco/8IYY4jIlar2AWPM80T0ZgAPEdH/TS0gReHw9wAWjTE/IaKjAP4cSchEEcYgFIMx\n5pcrfsXzAN6ee/22dN9g8F0jIvoHIjpojPlB6ra+4PiO59PHF4joy0jCB0NXDJy+Mfj+wyB4DYwx\nL+eeP0BEp4hovzFG6yhtI6IvaSgp4ZsAriWia4hoAcBHAdzfcZva5H4AN6fPbwaww8sioj1E9Kbs\nOYBfQTKwP3Q4feN+AL+dZpS8F8CPc6G5sRC8TkT0FiKi9Pn1SOTPhdZbKhsRfWkQHoMPIvpNAH8E\n4ACAvySix4wxv0pEbwXwBWPMUWPMa0T0CQAPApgAOGOMebrDZrfNHwD4MyK6DUnV2o8AQP4aIRl3\n+HJ6X+8C8KfGmL/qqL2t4eobRHQsff9uAA8AOArgGQCXANzSVXu7gnmdPgxglYheA/AzAB81I5th\nS0RfAnADgP1E9ByA3wcwD8jqSzrzWVEURZlBQ0mKoijKDKoYFEVRlBlUMSiKoigzqGJQFEVRZlDF\noCiKosygikFRFEWZQRWDoiiKMoMqBkVRFGWG/w/jWxRAyjIl9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe051a49b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Assign the resulting classification to vectors\n",
    "part1 = points[:,belong==0]\n",
    "part2 = points[:,belong==1]\n",
    "\n",
    "#Point cloud\n",
    "plt.scatter(part1[0,:],part1[1,:], color='b')\n",
    "plt.scatter(part2[0,:],part2[1,:], color='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## References\n",
    "\n",
    "Book :\n",
    "Lectures on Spectral Graph Theory by FanR.K.Chung: http://www.math.ucsd.edu/~fan/cbms.pdf\n",
    "\n",
    "Article:\n",
    "Ng, A. Y., Jordan, M. I., & Weiss, Y. (2002). On spectral clustering: Analysis and an algorithm. Advances in neural information processing systems, 2, 849-856.\n",
    "http://ai.stanford.edu/~ang/papers/nips01-spectral.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}