notebooks/3._Cluster_GPS_data_with_PCA_and_tSNE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Welcome\n",
    "Welcome to clustering GPS data using visualization techniques.\n",
    "Robert de Munter, Ahmet Taspinar and Rob Romijnders collaborated to make this repo.\n",
    "\n",
    "## Goal\n",
    "\n",
    "  * We aim to visualize and understand messy data such as GPS trajectories. The data consist of 18000 trajectories in GPS coordinates. Somehow we need insight into this data\n",
    "  * As a second goal, we aim to understand dimensionality reduction. Techniques such as PCA and tSNE allow for 2D visualizations of high dimensional data.\n",
    "  \n",
    "## Target\n",
    "We assume little knowledge of ML or Python in this notebook. Coding challenges require maximum two lines of code. Theoretical questions require insight, no mathematics. If you already know PCA/tSNE, then you could better join another group. If data excites you and you want to learn visualization techniques, then this is the ideal notebook for you!\n",
    "  \n",
    "## Approach\n",
    "\n",
    "  * __First__ take time to read the different functions. Run the standard implementations and discuss with groupmates what you see/expect.\n",
    "  * __Second__ we have some assignments to challenge your understanding. \n",
    "\n",
    "## Assignments\n",
    "\n",
    "  * Discuss with your group the differences between PCA and tSNE: \n",
    "    * What do they do? \n",
    "    * When to choose one over the other?\n",
    "  * Look at the PCA plots: \n",
    "    * What insight do you gain on the data?\n",
    "  * Consider the code for PCA. Explain to your groupmates how it works: \n",
    "    * What does the abbreviation _SVD_ mean? \n",
    "    * Why are we using SVD? (bonus question: and why truncated SVD?)\n",
    "  * Look at the tSNE plots: \n",
    "    * What insight do you gain on the data?\n",
    "  * Consider the code for tSNE. Explain to your groupmates how it works: \n",
    "    * What does the _perplexity_ parameter do?\n",
    "    * Bonus question: why do we define the _gradient norm_ and _number of iterations_?\n",
    "  * Coding challenge\n",
    "    * Pick either tSNE or PCA and improve the plot: what story does the data tell you? how can you clarify the story?\n",
    "    * Tune the hyperparameters of the techniques. For PCA, play with plotting different _principal components_. For tSNE play with the _perplexity_ and _angle_.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[data here](https://www.dropbox.com/s/xm5vtriolvmnw4h/03_cluster_data.tar.gz?dl=0) or [here](https://www.dropbox.com/s/p1gdt070a6hl67u/03_cluster_data.zip?dl=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Import libraries and define some globally useful parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "walk=r,train=g,subway=b,taxi=y,bus=m,bike=c,UNK=k\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append('../cluster')\n",
    "import numpy as np\n",
    "from sklearn.manifold import TSNE\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import zscore\n",
    "from util_cluster import some_cluster\n",
    "import glob\n",
    "import pandas as pd\n",
    "from collections import OrderedDict\n",
    "\n",
    "INPUT_FOLDER = '../processed_data/'\n",
    "IM_FOLDER = '../cluster/im/'\n",
    "HOT = OrderedDict()\n",
    "HOT['walk']= np.array([1, 0, 0, 0, 0, 0])\n",
    "HOT['train']= np.array([0, 1, 0, 0, 0, 0])\n",
    "HOT['subway']= np.array([0, 0, 1, 0, 0, 0])\n",
    "HOT['taxi']= np.array([0, 0, 0, 1, 0, 0])\n",
    "HOT['bus']= np.array([0, 0, 0, 0, 1, 0])\n",
    "HOT['bike']= np.array([0, 0, 0, 0, 0, 1])\n",
    "HOT['UNK']= np.array([0, 0, 0, 0, 0, 0])\n",
    "\n",
    "\n",
    "colors = ['r', 'g', 'b', 'y', 'm', 'c', 'k']\n",
    "lbl_clr = zip(HOT.keys(), colors)\n",
    "title = ''\n",
    "for lbl, clr in lbl_clr:\n",
    "    title += lbl + '=' +clr + ','\n",
    "print(title[:-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Functions to munge the data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def munge():\n",
    "    def clean_label(label):\n",
    "        # Pandas treats NaN's as floats. So in the next line we, effectively, check for NaN\n",
    "        if not isinstance(label,float):\n",
    "          return label.lstrip(',').rstrip(',').replace(',,', ',')\n",
    "\n",
    "\n",
    "    headers_metadf = ['trajectory_id', 'start_time', 'end_time', 'v_ave', 'v_med', 'v_max', 'v_std', 'a_ave', 'a_med', 'a_max',\n",
    "                      'a_std','labels']\n",
    "\n",
    "    list_df_metadata = []\n",
    "\n",
    "    for file in glob.glob(INPUT_FOLDER + \"*_metadata.csv\"):\n",
    "        df_metadata = pd.read_csv(file, index_col=0)\n",
    "        list_df_metadata.append(df_metadata)\n",
    "\n",
    "    \n",
    "    df_metadata = pd.concat(list_df_metadata).dropna(subset=['v_ave', 'v_med', 'v_max', 'v_std', 'a_ave', 'a_med', 'a_max','a_std'])\n",
    "    \n",
    "    X = df_metadata.as_matrix(['v_ave', 'v_med', 'v_max', 'v_std', 'a_ave', 'a_med', 'a_max','a_std'])\n",
    "    y = df_metadata['labels'].values\n",
    "\n",
    "    N = X.shape[0]\n",
    "    D = 6\n",
    "\n",
    "    # Rewrite the label strings to hot-vectors\n",
    "    # Note that labels are not mutually exclusive!\n",
    "    Y = np.zeros((N,D),dtype=np.int16)\n",
    "    for iy in range(N):\n",
    "        lbl = y[iy]\n",
    "        # Pandas treats NaN's as floats. So in the next line we, effectively, check for NaN\n",
    "        if not isinstance(lbl, float):\n",
    "            for key, value in HOT.items():\n",
    "                if key in lbl:\n",
    "                    Y[iy] += value\n",
    "    return X,Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Functions to preprocess data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def remove_outliers(X,y):\n",
    "    \"\"\"Start on visualization.\n",
    "    Remove all data outside (MAX_STD) standard deviations\"\"\"\n",
    "    MAX_STD = 2.0\n",
    "    mean = np.mean(X, 0)\n",
    "    std = np.std(X, 0)\n",
    "    mask = []\n",
    "    for i, x in enumerate(X):\n",
    "        if np.isnan(np.sum(x)):\n",
    "            continue\n",
    "        elif np.any(np.abs((x - mean) / std) > MAX_STD):\n",
    "            continue\n",
    "        else:\n",
    "            mask.append(i)\n",
    "    X = X[mask]\n",
    "    y = y[mask]\n",
    "    X = zscore(X)\n",
    "    return X,y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### The working horses: PCA and tSNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def PCA(X,y):\n",
    "    \"\"\"PCA\"\"\"\n",
    "    PCA_model = TruncatedSVD(n_components=2)\n",
    "    reduced = PCA_model.fit_transform(X)\n",
    "    return reduced,y\n",
    "\n",
    "def tSNE(X,y,perp=30):\n",
    "    \"\"\"Subsample to prevent memory error\n",
    "    - perp: the perplexity for the tSNE algo\"\"\"\n",
    "    N = X.shape[0]\n",
    "    ind = np.random.choice(N,10000)\n",
    "    X = X[ind]\n",
    "    y = y[ind]\n",
    "\n",
    "    tSNE_model = TSNE(verbose=2,perplexity=perp,min_grad_norm=1E-07,n_iter=300,angle=0.6)\n",
    "    reduced_tsne = tSNE_model.fit_transform(X)\n",
    "    return reduced_tsne,y\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Our cockpit\n",
    "Here we load the data and preprocess.\n",
    "Next we run the algorithms and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----Munge the data----\n",
      "----Start PCA----\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEICAYAAABVv+9nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FPX9x/HXJ+GUQwJE5BREBA8ONSIoWqqCCK1YtYpX\nwWrRIvU+8ELqQRW11pOfeAFWOYpSkYqKqChU1CA3aFUE5L7lDJDk8/tjJrKEJGSTTTZZ3s/HI4/M\ntd/5zHd3Pvud78zOmLsjIiKJJyneAYiISMlQghcRSVBK8CIiCUoJXkQkQSnBi4gkKCV4EZEEVS4S\nvJkNMrN/hsNNzczNrEK84zoQM2tiZtvMLDnesUTLzD4xs2viHceBmNlwM3soTuteYGadS7D8uG1b\ncUXus/nM/6XuDrTswSoW73+5SPDxYGZ9zGxaccpw92XuXt3ds2IVlxSdmXU2s+WxKs/dj3P3T2JV\n3sGkrNZd2Hg8Kte0yAZm53CZ53MtM83M+oTD++QOM6tpZtPN7E0zq1QKm/ELJXigqEcD5bFlLiLF\nth240syaHmhBM0sBpgBLgUvcfXfJhravEk/wZnaVmb0TMf6dmf0rYvwnM2tnZk+Fw1vMbKaZnV7I\n8i80syVmdnwUMeV081xtZsuAj3LNPwb4P6Bj2MWyOZw+3MyGmtm7ZrYd+LWZ9TCzWWHcP5nZoDzW\nUyEc/8TMHgy/zbea2QdmVrewcUeU+wczW2pmG8zsvnD7z85n2T5mtjhc349mdnk4fZ/D4ny6vpqb\n2Zfhtr1tZrXDZUeY2a3hcMPwddeH483NbKOZJZlZiplNNLN1ZrYpHG4ULvd7M5uZK9ZbzOztaOsD\nqGtmk8NtnGpmR+RRD9WASUCD8D3dZmYNzKy9mX1uZpvNbJWZPZvTyjKzU81svZk1DsfbhtvRKhzP\nt97zYmaHm9kOM6sTMe3EsH4qRrNteb1fFtGtZmZHhcv/HG7DmMLGGVGem1m/cJ/dGn52m5vZf8PP\nxFgruEVaxczGhK/92szaRpSdZ92ZWUUzG2Vhazf8HA0wsx/Cz/vYnM9hlNvSKYx7c7if9om2jAib\ngeHA/QdYZyrwMTAfuMLdM4u6QjOrYWYfm9nTZmaFfV1ptOCnAqeHb1QDoBLQEcDMjgSqA3OBr4B2\nQG3gDeBfZlaloILN7CrgUeBsd58fTttcwN+AXEX8CjgGOCdyorsvAq4DPg+7WGpFzL4MeBioAUwj\n+Db/A1AL6AH82czOLyDsy4CrgMPCurgtYnsOGLuZHQs8D1wO1AcOBRrmUz/VgKeBc929BnAqMLuA\n2HL7A/DHcD2ZYVkQvKedw+FfAYuBMyLGP3P3bILP16vAEUATYCfwbLjcBKCZBV+mOa4ERoaxDyio\nPnLFeTnwIFA33L7Xc2+Iu28HzgVWhu9pdXdfCWQBN4ev7QicBfQLX/Nf4AVghJlVBf4J3Ofu3+Qu\n38wuO8D718TdVwOfABfn2ubR7r4nd5mF3bZ8PAh8AKQAjYBnImKdW0Ccz+cq5xzgJKADcAcwDLgC\naAwcD1xaQAw9gX+xd5/+dwFfZIR1/G9gF3Bx2Nr9C3A+weeqAbAJeC7iNYXZZ44g+HJ/BkglyDPR\n7Ad5eRi40Mxa5jO/NsF7/Tnwx3B/yIk5ms82FjQIpgDT3f0Gj+b+Mu5e4n/AT8CJQC+CD8iXQCuC\nRDchn9dsAtqGw4OAf4bDTQEnSIwLgUZFiCenjCMLWKYPMC3XtOHAyAOU/Q/gyVzrqRCOfwLcG7Fs\nP+C9KGMfCIyKGD8E2E3wJZd72WoErY0Lgaq55v1SpwXE+kjE/GPD9SQDzcP3J4ngSOdaYHm43Ajg\nlnxibwdsihgfCjwcDh8Xllk5yvoYTpAgc8arEyTtxnks2zknzgLKuwkYHzFeEZgJzAPeAyxi3pK8\n6v0A5V9CsKMS1uVqoH2025b7/Yp4z64Jh0cS7GtR7x8R5TlwWsT4TODOiPEngH/k89pBwIyI8SRg\nFXB67roLl51A0HB4OlcdLwLOihivD+yJ3O5CbMddke9pIbb5qPz2lcjPEDAEGBMOTwP6hMN9gK1h\nnKcUo/6HA68QHAHcXpQySqsPPqfFd0Y4/AnBN/KvwnHM7DYzWxQeUm4maJkW1H1xO/CcuxfnpNlP\nxX2NmZ0SHjqtM7OfCVr+BcW9OmJ4B8FOG40GkTG4+w5gQ14LetBqvSSMaZWZ/cfC7oVCitzWpQTJ\nrq67/0Bw5NIOOB2YCKwMWzOR7+khZvaCBd1JW4BPgVq299zFCOCy8JDzSmCsu++KIr794nT3bcBG\ngno6IDM72oKuo9VhjIOJeP88aFkPJ2itPuHhnlcMbwPHmlkzoAvws7t/WcDyRd22OwADvrTgipU/\nFjHeNRHDO/MYL+jzGxl7NrCc/GPvALQhaFRE1vERwPiI1u0igi+5eoXeguAL8YdCLptF8DmPVJEg\nWef2KHBOZNdThDkEjdBJZnZCYQPNQw+gKkFDKmqlneBPD4enEpHgLehvv4Pg0DXFgy6Rnwk+oPnp\nCtxrZhdGTrS9/at5/d2dq4yCdtb85uWe/gZB66Oxux9K8EYUuo+sCLGvIjjkznlNVaBOngUC7v6+\nu3chaPl8A7wYztpO0PrPcXgeL28cMdyE4EO+PhyfClwEVHL3FeF4b4IugZzD31uBlgStmJrs7cax\nMLYZBEcFpxN0Xb0WsV13F1Qf+cVpZtUJDo9X5lUdeUwbSlAvLcIY7ybi/TOzhgR9ra8CT5hZ5TzK\nwMwuP8D71yTc5gxgLEE3x5WR25yP/LZtezg5z/fQ3Ve7+5/cvQHBEdbzFl4dEib8/OIsUiIpROxJ\nBJ/bvN4XCLqT/gZMMbPI5P0TQRdjrYi/KuFnrrD7zE8ER52FsYzg6ChSM4IGzj7cfQPBEfuDeRXk\n7k8BjwCTLeIcYZSf7RcJjhzftaDLNSqlmeB/TdBNsBz4DOhGkJhmEfRnZwLrgApmNhCoeYAyF4Rl\nPGdm5+VM9L39q3n9DS6oQAtO/PQJR9cAjezAlzXVADa6e4aZtSdIVEVSyNjHAb+14ARgJYLDx8iE\n1NnMPByuZ2Y9ww/GLmAbkNMXOBs4w4Jr9Q8lOIzN7QozO9bMDgEeAMb53ks+pwL9CVrlEByV9Sfo\n1spZpgZBK2+zBSfG8jopNZKgX36Pu/9yaZm7Dy6oPnKV0d2Ck2iVCHa2Ge7+U1gHn9jeE99rgDrh\n9uaoAWwBtoVHN3+OqEsjaL2/DFxN8OWa3878+gHev2W5trkPcB77fqnlnDhteqBtc/d1wAqC9yg5\nbKE3jyjr9xae0Cbo+nLC996DSxTzi/O6vLavMMLYO0dMOsnMLrDgRPBNBJ/BGfm93t2HEDSYptje\niw/+D3jY9p5cTjWznhGvKcw+8zpwtpldbGYVzKyOmbULy+tjZksiwhhD0HBsZMF5w7OB3xLsd3n5\nO8G5rWPymhlu01PAhxb210f52YZgv/oWeCds0OUVd55KJcG7+/8Ikstn4fgWghNz08Nk8D7Bt9T/\nCL4pMyhE94m7zwF+A7xoZucWJ8ZwB6rD3g/gRwRfIqvNbH2+Lwz60R8ws60E/eNjixPHgbj7AoIT\nT6MJEs42YC3BzgNBq+m/4XAScAtBq2kjwRHTn8NyJhN8mOcS9K1OzGN1rxEkuNVAFeCGiHlTCZJj\nToKfRtCa/DRimX8QHF6uJ6jX9/JZx/EEJzCL6g2CL4+NBCcEr4iY1xiYDuDBydFRwGILDvkbEBxG\nX0bQZ/oiQZ3kuIHgZPh9YbfBVcBVVsgrvPLj7tMJku3X7h7ZMmxM8PlfUcht+xNBV+UGgnMY/42Y\ndzLwRdginADc6O6LixN3QSy40mgrwbmKHG8TdBFuIjhaucDzP5kMgLs/SHCi9cOwUfAUQfwfhPvY\nDOCUaGILv1y7ExxRbiRo3OR0q/zy+Qg9QFCP08K4hwCXe3gRRx5lbwmXyffKnnCbXiL44irskUTk\n6x3oS9DF9bYFF5/kjjtPVvwuxcRgZp2A6929oKsCypzwsH0zQRfDj2b2EvAvd38/zqEVStgiWQuc\n6O7fxbjsRgT9+qfGstxYMLOPgDfc/aWIafcC69z9hfhFVjRmdgVwnLvndSRYZpnZBwRffoviHUs0\nChu3Enw5ZGa/JbhsygiuZDiFIEGWuzfTzG4BfuPuZ8Y7ltJiZicDkwnO22yNdzySuMr8/VwkTz0J\nujYMSAd6ldPkvoRgGwr63UBCMbMRBNt7o5K7lDS14EVEEpTuRSMikqDi0kVTt25db9q0aTxWLSJS\nbs2cOXO9u6cWdvm4JPimTZuSnp4ej1WLiJRbZrbfD64Koi4aEZEEpQQvIpKglOBFRBKUEryISIJS\nghcRSVBK8CIiCUoJXkQkQSnBi4iUgo8++ogTTjiB//znP6W2TiV4EZFS8OGHHzJ79mzeey+vxyKU\nDN1NUkSkFNxzzz20bt2aHj16lNo6leBFREpBtWrVuPTS0n2ekLpoREQSlBK8iEiCUoIXEUlQSvAi\nIglKCV5EJEEpwYuIJCgleBGRBFXoBG9mVczsSzObY2YLzOyv4fRmZvaFmX1vZmPMrFLJhSsiIoUV\nTQt+F3Cmu7cF2gHdzKwD8CjwpLsfBWwCro59mCIiEq1CJ3gPbAtHK4Z/DpwJjAunjwDOj2mEIiJS\nJFH1wZtZspnNBtYCk4EfgM3unhkushxomM9r+5pZupmlr1u3rjgxi4hIIUSV4N09y93bAY2A9kCr\nKF47zN3T3D0tNTU1yjBFRCRaRbqKxt03Ax8DHYFaZpZz07JGwIoYxSYiIsUQzVU0qWZWKxyuCnQB\nFhEk+ovCxXoDb8c6SBERiV40twuuD4wws2SCL4ax7j7RzBYCo83sIWAW8HIJxCkiIlEqdIJ397nA\nCXlMX0zQHy8iImWIfskqIpKglOBFRBKUEryISIJSghcRSVBK8CIiCUoJXkQkQSnBi4gkKCV4EZEE\npQQvIpKglOBFRBKUEryISIJSghcRSVBK8CIiCUoJXkQkQSnBi4gkKCV4EZEEFc0j+xqb2cdmttDM\nFpjZjeH0QWa2wsxmh3/dSy5cEREprGge2ZcJ3OruX5tZDWCmmU0O5z3p7o/HPjwRESmqaB7ZtwpY\nFQ5vNbNFQMOSCkxERIqnSH3wZtaU4PmsX4ST+pvZXDN7xcxS8nlNXzNLN7P0devWFSlYEREpvKgT\nvJlVB94EbnL3LcBQoDnQjqCF/0Rer3P3Ye6e5u5pqampxQhZREQKI6oEb2YVCZL76+7+FoC7r3H3\nLHfPBl4E2sc+TBERiVY0V9EY8DKwyN3/HjG9fsRivwPmxy48EREpqmiuojkNuBKYZ2azw2l3A5ea\nWTvAgSXAtTGNUEREiiSaq2imAZbHrHdjF46IiMSKfskqIpKglOBFRBKUEryISIJSghcRSVBK8CIi\nCUoJXkQkQSnBi4gkKCV4EZEEpQQvIpKglOBFRBKUEryISIJSghcRSVBK8CIiCUoJXkQkQSnBi4gk\nKCV4EZEEFc0j+xqb2cdmttDMFpjZjeH02mY22cy+C/+nlFy4IiJSWNG04DOBW939WKADcL2ZHQsM\nAKa4ewtgSjguIiJxVugE7+6r3P3rcHgrsAhoCPQERoSLjQDOj3WQIiISvSL1wZtZU+AE4Augnruv\nCmetBurl85q+ZpZuZunr1q0rympFRCQKUSd4M6sOvAnc5O5bIue5uwOe1+vcfZi7p7l7WmpqapGC\nFRGRwosqwZtZRYLk/rq7vxVOXmNm9cP59YG1sQ1RRESKIpqraAx4GVjk7n+PmDUB6B0O9wbejl14\nIiJSVBWiWPY04EpgnpnNDqfdDTwCjDWzq4GlwMWxDVFERIqi0Ane3acBls/ss2ITjoiIxIp+ySoi\nkqCU4EVEEpQSvIhIglKCFxFJUErwIiIJSgleRCRBKcGLiCQoJXgRkQSlBC8ikqCU4EVEEpQSvIhI\nglKCFxFJUErwIiIJSgleRCRBKcGLiCQoJXgRkQQVzSP7XjGztWY2P2LaIDNbYWazw7/uJROmiIhE\nK5oW/HCgWx7Tn3T3duHfu7EJS0REiqvQCd7dPwU2lmAsIiISQ7Hog+9vZnPDLpyU/BYys75mlm5m\n6evWrYvBakVEpCDFTfBDgeZAO2AV8ER+C7r7MHdPc/e01NTUYq5WREQOpFgJ3t3XuHuWu2cDLwLt\nYxOWiIgUV7ESvJnVjxj9HTA/v2VFRKR0VSjsgmY2CugM1DWz5cD9QGczawc4sAS4tgRiFBGRIih0\ngnf3S/OY/HIMYxERkRjSL1lFRBKUEryISIJSghcRSVBK8CIiCUoJXkQkQSnBi4gkKCV4EZEEpQQv\nIpKglOBFRBKUEryISIJSghcRSVBK8CIiCUoJXkQkQSnBi4gkKCV4EZEEpQQvIpKglOBFRBJUoRO8\nmb1iZmvNbH7EtNpmNtnMvgv/p5RMmCIlyB1WroSMjHhHIhJT0bTghwPdck0bAExx9xbAlHBcpPzI\nyICzzoKGDeGII2DOnHhHJBIzhU7w7v4psDHX5J7AiHB4BHB+jOKSOMj27HiHUPreegs+/jgYXrsW\nHnggvvGIxFBx++DrufuqcHg1UC+/Bc2sr5mlm1n6unXrirlaiaVdmbu4YMwFVHywIs2fbs68NfPi\nHVLpqVKl4HGRcixmJ1nd3QEvYP4wd09z97TU1NRYrVZyyc7O5sMPP2Tr1q2Ffs3zXz3P+G/Gk+3Z\nLN60mL4T+5ZghGXM+edDnz5QqRK0bQuDB8c7IpGYKW6CX2Nm9QHC/2uLH5IUx9ChQ+nSpQv9+vUr\neMHJk2HUKNi2jfU71u8za932g+gIKykJXn0Vdu2C2bODfniRBFHcBD8B6B0O9wbeLmZ5Ukxt27al\nSZMmdOzYMf+F7r8funaFyy6D00/niqMvomblmr/M/nPan0shUhEpaRUKu6CZjQI6A3XNbDlwP/AI\nMNbMrgaWAheXRJBSeJ06dWLp0qUFLzRy5N7h2bM5ZnUms66dxYeLP+SLd7+g6uyqcGrJxikiJc+C\nrvPSlZaW5unp6aW+Xgl16wbvvx8MV6gADz0Ed97JDz/8wFFHHQXArl27qFSpUhyDFJHczGymu6cV\ndvlCt+AlgYwcGXTPTJkCmZkwYACkpNC8b1+GDBlCSkqKkrtIAtCtCg5Ghx0GnTvvO+2rrwC4/fbb\nueaaa0o/JhGJOSX4g1XXrlCxYjBsBj16xDceEYk5ddEcrNq3h+nT4aOP4OST4cwz4x2RiMSYEvzB\n7OSTgz8RSUjqohERSVBK8HLQW7NmDSNHjiRDtwuWBKMELwe9m266id69e/PMM8/EOxSRmFIfvBz0\nzj//fH788UfOOuuseIciElP6JauISDkR7S9Z1UVTTLt37yYzMzPeYYiI7EcJvhg2bdpEkyZNaN26\nNdnZwdOQHn/8cUZG3sxLRCROlOCLITMzk507d7J9+3bcnW+//Zbbb7+d3r17s3r16niHV2b9/PPP\n3HvvvcycOTPeoYgkNCX4YkhNTeXHH39kwYIFJCcn06JFCw4//HAAPs55zmdRuMO110LVqtCuHSxZ\nEpuAy4iRI0fy8MMPc+edd8Y7FJGEpqtoiql27dq/DCclJfHOO+8wdepULrjggqIX+s47MGxYMDxn\nTnC3x9Gjixlp2XHRRRcxa9YsLr300niHIpLQlOBjLC0tjbS0Qp/kztvOnfuO79hRvPLKmPr16/PK\nK6/EOwyRhBeTLhozW2Jm88xstpnp+sfiOv/8vTf/qls3eMSeiEiUYtmC/7W7rz/wYnJAlSsHD8Ve\nuRLq1An64kVEoqSTrAXYvHlziZa/fv16Jk2aRJ4/NktKgkaNlNxFpMhileAd+MDMZppZ37wWMLO+\nZpZuZunr1q2L0WpLzuuvv05KSgqDBg0qsXVcddVVdO/evcDr5rdv386YMWPYunVricUhIokpVgm+\nk7ufCJwLXG9mZ+RewN2HuXuau6elpqbGaLVFt3x5cLHKjz/mPX/Pnj1A8EvVHDVq1CApKYkZM2YU\neb2TJ09m6NChuDtnnnkmLVu2pF27dvku/8ADD9CrVy8GDhxY5HWKyEHK3WP6BwwCbitomZNOOsnj\naeZM95o13cG9alX3jz/Oe7lVq1Z5dnb2L+MERyr+xBNPFHndtWrVcsBnzJhRqOU/+OADb9Omjf/n\nP/8p8jpFJDEA6R5FPi72SVYzqwYkufvWcLgr8EBxyy1Jzz0HW7YEwzt3wpNP7v8MauCXHy3lGDt2\nLHPnzuWWW24p8roHDhzIwoULadu2baGW79KlC3PmzCny+kTk4BWLq2jqAePNLKe8N9z9vRiUGzsZ\nGdC/P3z5JXTvTs0afwPsl9k1awb/x48fz4oVK+jfv/9+RUyYMIHp06fz0EMPFSuUm2++uVivj+Tu\nbMrYREqVFML6FxH5RbETvLsvBgrXHI2XRx+Fl18OhufN455/tGZ6jRP5ausxNE1+g18fnwH8kUsu\nuYQ9e/bQuXNnjj/++H2KuOuuu1i4cCGnn346F154YelvQy4/bPyBrv/syuJNi2l9WGs+uPIDDq9+\n+IFfKCIHjYPjl6wrVuwzWufjf9Fw6xXspCL/y9rD1QPgjAvPYMiQISxbtoxWrVrtV8STTz7J1KlT\n6d69e2lFXaA7P7yTxZsWAzBv7TwemPoAz/d4Ps5RiUhZUi4T/N/+9jf++9//Mm7cOCpXrnzgFyQF\nFwutBaaY8dInnzACMFK5kUy8STUaN27MjTfeyKZNm6hQYf9q6dq1K127do3thhTDll1bChwXESmX\nP3S69957mThxIh+fdhr06QNLl/4yb8tXW1j+7HK2frIMxo/n69Gj6f7aa1xMLRqRxGXuXLxnC+n0\npDk/8CZr+GbH52zfXpn+/ftTp04devbsyc8//xy/DSyEmzvcTKXkSgBUq1iN/u33P28gIge3ctmC\n73XxxeweN44zZs6EmTMhPR3mz2fDexuY95t5zD9+Cz1W/5kaa1YyBpgEQHDDru7nQss7nD9d8hC7\n1lYBYOH6erz6Kr88tGPChAm89dZbXHXVVXHZvsI4t8W5zLluDvPWzCOtQRrNUprFOyQRKWPKZYI/\n46STuHX0aKYDDYCjvvuOv/7wA5/98C1ZbUfx3OYvaLZmJT8BmypWIXlPJr/jfE6u14Smx38KpEPN\nPUGfTahiRXjuuefo2rUr06ZNK97tfktJq7qtaFV3//MFIiJQDh+6vWXLFs477zzWTp3KIcBMoB+Q\nfdNN/F/PniSvWcux/f7MCZUqsWrtWiZnZ9OQhtzDPRzDMZCUBU/dyKw9qdxxx3tkZlbm1FPhgw+g\nWrW969mzZw/33HMPLVu25Oqrr47FZouIFEvCP3T7hhtuYPnU5XQmiYrhtApA1/ALI/nbb/jbxo28\nuno1z2Vn8zDQlBV8wifBwtnJMK819ev/SPXqmzm7RxaffbZvcgf4+uuveeyxx2J63bqISGkqdwn+\nlHan8BiPcTi/5TOMxcBTQCZwyJBn6fHULFrQgiSgBXA38ClwI8eFJTjL53fiTzeMY3ONCcy8+kt+\n2p3BpEmTuOaaa8i5EdrJJ5/M4MGD9WAKESm3yleC37qVtvc9RW2q0ZsarOYoagEPA19Nm0O1SQMZ\nv3E8x7CIj2lDFvAFwQ1kdvE23/M9YNiZs9le5Xz4qS+bZvybSRs2MHjwYF5++WX+/e9/A8Hj9+66\n6y4uuuiiuG2uiEhxlKsEn/XEExy17Tuq8RwZVZZSw9ZyEnAvxnguZB0tAMgmmVfozyCgA8GNcWpa\nCtsrbGJE9Wep3PYlDt0a/vhpwQKOrFqVxx9/nPvuu4/MzEwuvfRSNmzYEJ+NFBGJkXJzFc3u7Gze\n/edYvqlcneN3TaRyRja3AWuoAmTwM8OpxEh28znQng85iy5AdUuigmfTv9oXzNs2haRsOPZH2LwZ\nMOMPv/0tXWvXhlNO4ZRTTqF169bMnz+f4zofx6yUWdSsXJMHf/0gjWo2iuv2i4hEq9wk+BlPP0vT\nZd/zv8xM+uH89MucDGoBW4FssjmM19hJKx7iOtaeeipzl23m/rXf49u2ATDhYqjVoQavvnoPbdp0\n4cQTT9xnPYMGDeLKP1zJwGcH4hcFVxh9sfwLFvRboBt6iUi5Um66aDLGjmd8ixbc6ZHJPbAZyABS\ngfv5gFtIoQeTyd62g7nLF+G7d0OrVvylShXSklqRljaHif/+nB49erBs2bJ9ymrSpAkZOzPwjXsv\nH120fhGbM0r28X0iIrFWbhI8O7dwycKFBS6yBejH/xhENocDA+bOpjMOKbX5fdeu/C4jg8NnGFX/\ncAffvv02a1evZtNXX+1Txsknn8yEYRPo3aw3ned3BqBlnZbUqlKrZLZLRKSElJsumsY7N5Bfik2q\nUIFsM3YBP+zZQ/NwejJwM8b0zZtY+fTTnAmMS6nNhePGMRVYB7QcMQIibv/784yfqfHnGvTJ7APA\nOTXO4bKbLlP3jIiUO+WmBd9sxBskp1SjClX2m5edmQmNm5C5Zw/Tc82rUcHJdCejUiWOql+fVpdc\nAkBtoCVArsS96f1NeObe7pluy7rR5NAmsd0YEZFSEJMEb2bdzOxbM/vezAbEoszcqiRVpPWw4dzV\n6qm8F0g7iSpt2tC4eXMWNW5MhhkTgXMy4SGg14A+fLdyJcf16wdhkqdxYxg8GIAvv/yS2267jV3N\ndu1TbPWjkktic0RESlyxE7yZJQPPAecCxwKXmtmxxS13P5Mn43UyOXX7G/vN6msGn39Oxty5nHXZ\nZRz74ousP7Ey44HdwF+Bp14aQ2ZmJiQn46NGcddtt3FG06a8MG0aAPfffz9PPPEE7254lxZ/2kEK\n6TRiLM0m9oQ1a2K+OSIiJS0WffDtge/DR/dhZqOBnkDBZ0SjdfLJ1MlcxGGZW/ebNc6dSilV2X3q\nJfirr2KPD+HWGruYlASdK8KMXbB81c80afIsFSveRM2au5g//31gHnPmzOHaa6/l7rvvplmzZpx5\n5hXYf26lLa8Fha8AJk0K7jsvIlKOxKKLpiHsc+Xi8nDaPsysr5mlm1l6zv1eotKlCytXVKduxcWk\n5pq1EcgXUVM6AAAOTUlEQVSe+z/6TjiGGltXUtV3MXdLJbZmw3epyWQY4GewatVNLFsG8+dXITV1\nKr169eK114JEfvrpp3Pooc/Trl0qDaeM5D4e2LuCI4+MPl4RkTgrtato3H0YMAyC2wUXpYz0hd9Q\nf9lmsvOYNwjYvdOpmTyHrEseY+m3S2HGVA6peiT4d8C+J0q3bk1h1KhRv4wvXgyPPLJ3/kPcx5/S\nZtPkj2fDGWcUJVwRkbiKRQt+BdA4YrxROC3mjtgc/NjouDzm9QcOZw0rth3P6hHr2TljKmbG8uXL\nAejcOZN69fZ+NVx55b6vz8rav8yf/vEM9QYN4rzzzovRFoiIlJ5YJPivgBZm1szMKgG9gAkxKHc/\nfvHF7D6kAtcM3n/eXdTiBt4ERlGRlcHyfiG7d18FtKJhw6a0bfsn4CZ69hzDY4/t25ffogVcd93e\n8f794dBDN7JhwwYWL15cEpsjIlKiip3g3T2ToAH9PrAIGOvuC4pbbl6OnjaNASe157g8oh7KZjJZ\nR/2UB9jDbOAQ4GWysp4GFjF37mDOOONIatUawdtv96J//377lzEUFiyAhQvhmWfg+OOP55tvvuGz\nzz4ric0RESlRMbkO3t3fdfej3b25uz8cizLzsql5c/5ZqzK7GuQ9/+ha0O+Wb8KxHUBrgh6jDSxY\nYNx99z288847NG7cmA4dOuRZxrHHwjHH7B0/6qijSElJieFWiIiUjnLzS1aA3u3bs77fHWRlV95v\nXh3q8Kv20KkTnHfekXTp0pPgKvifgQzOOSf40WqnTp1YtmwZ119/PStXrqRTp0489NBDpb0pIiIl\nrlwl+CdatKB6796MfGa/qzDZyAY21oB3J8GECYuZM+dzUlOPxuwPXHTRKsaN27+8uXPnMn36dMaO\nHVvs2N577z369+/P5s2666SIlA3mXqQrFoslLS3N08OHZEcrOTkZz3YMz/NyyUOAXSSTRRb169dn\n1apV1KpVi02bNu2z3PLlyxk5ciQNGjTgtNNOo0WLFkWKJ0f79u356quvGD58OL179y5WWSIieTGz\nme6eVtjly1ULHuCOO+4gCSf3lenVgXrA7cBdnEnN5PqsX38lUJXNm1uwdu3afZZ/8MEHueeee5g3\nb16xkzvAkCFDuOOOO7jggguKXVZ59+mnnzJw4EB27NgR71BEDmrl5nbBOb788kuygHdIokZEG/5e\n4BpO5lDmUoHJ3JeVzOVZaYxjB5DNrFk7OOecveVcddVVrF27liuuuCImcXXu3JnOnTvHpKzy7qab\nbmLWrFm0bNmSyy+/PN7hiBy0yl0LPuea9A0css/0ZcC3LGEwu/geqEQWL3JNODeJ77+vts/yHTp0\nYPz48ZxwwgklH/RB5v7776dv375079493qGIHNTKXYLv1q0bAKmV972t7wsGPUjmfqAF8ApQid3h\n3AzGj7+ZIUOGlGaoB62ePXvywgsv6PJSkTgrdwn+vvvu47TTTuP7iyvtMz3LIbPC3tv6GjA66QKo\n9TxwFlOmPFXkBL9lyxamTJlCdnZep3VFRMqmcpfgGzRowDvvvMPMmseRlCv67VlOO67iaX7PIh7l\nT8kPYX36U7/5j/Tp04e33nqrSOu87rrrOPvss3nhhRdisAUiIqWj3CX4rVu3kpSURNX273Jcrrv4\nVkiG2dzMDYzlMW4ne8+j+GfO1uytvPrqq5xRyLtCLlq0iAsvvJBPP/0UgI4dO9KkSRPatGkT680R\nESkx5SrB79q1i1atWtG8eXOatdzJRefs203TqDZAe0j+FOgHSS9AE+jZqWdU6xk9ejRvvfUWw4YN\nA+Avf/kLS5cu5bTTTovNhoiIlIJydZmkmVGpUiVa7d5Ns8WzqNFs9z7zV22pAGRA1nUkJy2i92nV\neeXdbRxz3zF5F5iPG264gcqVK9OrV68YRi8iUrrKVQu+UqVKfHfrrXy4YQOHXXYeZLXaZ/6ePTWA\n84BFZGXDwPkVefONN7ntttv2WW7y5Mlcd911rF+/Ps/11KlTh7vvvpsj9SQnESnHylWCB6gwdCgW\n3l5h2WPf7DMvK2sTwa3oa1EBSO7ShQsuuIDKlfe9OdmDDz7ICy+8wJtvvlk6QYuIxEG5S/A02fvo\nvd55NsCfAc4nE3i9Vq08i3j00UcZMGCAumBEJKGVqz54AF56Cfr1g1WrmF6/PkzI/fCoV7jk2CQW\nLYRfvftunkV07NiRjh07lnysIiJxVKwWvJkNMrMVZjY7/Cv536Y3bEjGmDFc3qIFdO263+yjrRqj\nF85kDtAhDnfKFBEpK2LRgn/S3R+PQTmFNnDgQN544w3eGfXOfvM6+hHANLxCJezxUg1LRKRMKX99\n8EDv3r05/PDDyfKs/ebVrnQYw48azJ5vl4P62EXkIBaLBN/fzOaa2Stmlu/dpcysr5mlm1n6unXr\nirXC4447jlWrVlHvkHokkURb2gJQi1osqLeAjv+8gOcnvE61atWYOHFi1OVv3TqbpUv/xtq1eTwG\nSkSknDhggjezD81sfh5/PYGhQHOgHbAKeCK/ctx9mLunuXtaampqTIJ/5uVnyCabOcyhIhXZzGY+\nWvkRtTfXZtn8ZezYsYMlS5ZEVeaWLel8/XUHfvzxbhYu/D1Lluh5rSJSPh2wD97dzy5MQWb2IhB9\nc7kYJk7du7rkKsnUrlibdjvasaDbAi48/EIuf+dyTvrNSVGVuX79W7jvvRXx2rWjaNr03pjFLCJS\nWop7FU39iNHfAfOLF050br75Zrp160b79u0ZP34867etZ0HWAgD2rN5Dyozo70depUrTAsdFRMqL\n4l5FM8TM2gEOLAGuLXZEUTj66KOZNGkSACtXrqRu1bo02bH3h1AVa1eMusz69a9h+/Z5/PDDa2Rl\nHU7HjsNiFq+ISGkqVoJ39ytjFUhxNWjQgKXfL2XBBQvYNnsbtc+tTcPrG0ZdjlkScAPnnvssZlvY\nvbte7IMVESkF5e+XrAXYU2MPr5/6Ol3u78IZ3Qp37/e8tGjRgkceeYS6detSoUJCVZGIHETM4/Br\nz7S0NE9PT495uWPGjKFXr160adOGOXPmFOo1GRkZfPfdd7Ru3Trm8YiIxJKZzXT3tMIuXy5/6JSf\nHj16cOutt0b17NWrr76aNm3aMHr06BKMTESk9CVU/0P16tV5PMrbExx55JFUq1aNBg0alFBUIiLx\nkVBdNCIiieyg7qIREZG9lOBFRBKUEryISIJSghcRSVBK8CIiCUoJXkQkQSnBi4gkKCV4EZEEFZcf\nOpnZOmBpDIqqC6yPQTmlQbGWDMVaMhRryShurEe4e6EfiReXBB8rZpYeza+64kmxlgzFWjIUa8ko\n7VjVRSMikqCU4EVEElR5T/Dl6Xl6irVkKNaSoVhLRqnGWq774EVEJH/lvQUvIiL5UIIXEUlQ5TbB\nm1k3M/vWzL43swHxjieSmTU2s4/NbKGZLTCzG8Pptc1sspl9F/5PiXesOcws2cxmmdnEcLyZmX0R\n1u8YM6sU7xgBzKyWmY0zs2/MbJGZdSyr9WpmN4fv/3wzG2VmVcpKvZrZK2a21szmR0zLsx4t8HQY\n81wzO7EMxPpY+BmYa2bjzaxWxLy7wli/NbNz4h1rxLxbzczNrG44XuL1Wi4TvJklA88B5wLHApea\n2bHxjWofmcCt7n4s0AG4PoxvADDF3VsAU8LxsuJGYFHE+KPAk+5+FLAJuDouUe3vKeA9d28FtCWI\nuczVq5k1BG4A0tz9eCAZ6EXZqdfhQLdc0/Krx3OBFuFfX2BoKcWYYzj7xzoZON7d2wD/A+4CCPez\nXsBx4WueD/NFaRnO/rFiZo2BrsCyiMklX6/uXu7+gI7A+xHjdwF3xTuuAuJ9G+gCfAvUD6fVB76N\nd2xhLI0IdugzgYmAEfzarkJe9R3HOA8FfiS8OCBiepmrV6Ah8BNQm+DZxxOBc8pSvQJNgfkHqkfg\nBeDSvJaLV6y55v0OeD0c3icXAO8DHeMdKzCOoEGyBKhbWvVaLlvw7N15ciwPp5U5ZtYUOAH4Aqjn\n7qvCWauBenEKK7d/AHcA2eF4HWCzu2eG42WlfpsB64BXw+6kl8ysGmWwXt19BfA4QYttFfAzMJOy\nWa858qvHsr6//RGYFA6XuVjNrCewwt3n5JpV4rGW1wRfLphZdeBN4CZ33xI5z4Ov7Lhfo2pmvwHW\nuvvMeMdSCBWAE4Gh7n4CsJ1c3TFlqF5TgJ4EX0oNgGrkceheVpWVejwQM7uHoEv09XjHkhczOwS4\nGxgYj/WX1wS/AmgcMd4onFZmmFlFguT+uru/FU5eY2b1w/n1gbXxii/CacB5ZrYEGE3QTfMUUMvM\nKoTLlJX6XQ4sd/cvwvFxBAm/LNbr2cCP7r7O3fcAbxHUdVms1xz51WOZ3N/MrA/wG+Dy8AsJyl6s\nzQm+5OeE+1gj4GszO5xSiLW8JvivgBbhFQmVCE6qTIhzTL8wMwNeBha5+98jZk0AeofDvQn65uPK\n3e9y90bu3pSgHj9y98uBj4GLwsXKSqyrgZ/MrGU46SxgIWWwXgm6ZjqY2SHh5yEn1jJXrxHyq8cJ\nwB/Cqz46AD9HdOXEhZl1I+hWPM/dd0TMmgD0MrPKZtaM4ATml/GIEcDd57n7Ye7eNNzHlgMnhp/l\nkq/X0jz5EOMTGd0Jzp7/ANwT73hyxdaJ4PB2LjA7/OtO0Lc9BfgO+BCoHe9Yc8XdGZgYDh9JsGN8\nD/wLqBzv+MK42gHpYd3+G0gpq/UK/BX4BpgPvAZULiv1CowiODewhyDpXJ1fPRKcdH8u3NfmEVwZ\nFO9Yvyfov87Zv/4vYvl7wli/Bc6Nd6y55i9h70nWEq9X3apARCRBldcuGhEROQAleBGRBKUELyKS\noJTgRUQSlBK8iEiCUoIXEUlQSvAiIgnq/wES64rRlsWS1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f598859fe48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('----Munge the data----')\n",
    "X, y = munge()\n",
    "X,y = remove_outliers(X,y)\n",
    "N,D = y.shape\n",
    "\n",
    "#PCA\n",
    "print('----Start PCA----')\n",
    "X_red,y_red = PCA(X,y)\n",
    "plt.figure()\n",
    "y_color,y_s = some_cluster(y_red, colors)\n",
    "plt.scatter(X_red[:, 0], X_red[:, 1],c = y_color,s=y_s ,marker='o', linewidths=0)\n",
    "plt.title(title)\n",
    "plt.savefig(IM_FOLDER+'pca.png')\n",
    "plt.show(block=True)\n",
    "# If the plot doesnt show after above line, then change the MATPLOTLIB settings for your interpreter\n",
    "\n",
    "plt.figure()\n",
    "f4, ax4 = plt.subplots(1, D)\n",
    "for d in range(D):\n",
    "    ax4[d].scatter(X_red[:, 0], X_red[:, 1],c =y_red[:,d] ,marker='*', linewidths=0)\n",
    "    ax4[d].set_title('Mode%3i'%d)\n",
    "plt.setp([a.get_xticklabels() for a in ax4], visible=False)\n",
    "plt.setp([a.get_yticklabels() for a in ax4], visible=False)\n",
    "plt.savefig(IM_FOLDER+'pca_table.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "----Start tSNE----\n",
      "----------tSNE on perplexity 300----------\n",
      "[t-SNE] Computing pairwise distances...\n",
      "[t-SNE] Computing 901 nearest neighbors...\n",
      "[t-SNE] Computed conditional probabilities for sample 1000 / 10000\n",
      "[t-SNE] Computed conditional probabilities for sample 2000 / 10000\n",
      "[t-SNE] Computed conditional probabilities for sample 3000 / 10000\n",
      "[t-SNE] Computed conditional probabilities for sample 4000 / 10000\n",
      "[t-SNE] Computed conditional probabilities for sample 5000 / 10000\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-ad7af058ba13>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mper\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m300\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m400\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m500\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'-'\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'tSNE on perplexity %i'\u001b[0m\u001b[0;34m%\u001b[0m\u001b[0mper\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'-'\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m     \u001b[0mX_red\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_red\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtSNE\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mper\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m     \u001b[0my_color\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msome_cluster\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my_red\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolors\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-4-59f7050a8ad1>\u001b[0m in \u001b[0;36mtSNE\u001b[0;34m(X, y, perp)\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0mtSNE_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTSNE\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mperplexity\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mperp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmin_grad_norm\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1E-07\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn_iter\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m300\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mangle\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m     \u001b[0mreduced_tsne\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtSNE_model\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_transform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     17\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mreduced_tsne\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/sklearn/manifold/t_sne.py\u001b[0m in \u001b[0;36mfit_transform\u001b[0;34m(self, X, y)\u001b[0m\n\u001b[1;32m    882\u001b[0m             \u001b[0mEmbedding\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mtraining\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlow\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mdimensional\u001b[0m \u001b[0mspace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    883\u001b[0m         \"\"\"\n\u001b[0;32m--> 884\u001b[0;31m         \u001b[0membedding\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    885\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0membedding_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0membedding\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    886\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0membedding_\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/sklearn/manifold/t_sne.py\u001b[0m in \u001b[0;36m_fit\u001b[0;34m(self, X, skip_num_points)\u001b[0m\n\u001b[1;32m    764\u001b[0m                 \u001b[0mneighbors_nn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mneighbors_nn\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    765\u001b[0m             P = _joint_probabilities_nn(distances, neighbors_nn,\n\u001b[0;32m--> 766\u001b[0;31m                                         self.perplexity, self.verbose)\n\u001b[0m\u001b[1;32m    767\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    768\u001b[0m             \u001b[0mP\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_joint_probabilities\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdistances\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mperplexity\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/sklearn/manifold/t_sne.py\u001b[0m in \u001b[0;36m_joint_probabilities_nn\u001b[0;34m(distances, neighbors, desired_perplexity, verbose)\u001b[0m\n\u001b[1;32m     93\u001b[0m     \u001b[0mneighbors\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneighbors\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mint64\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     94\u001b[0m     conditional_P = _utils._binary_search_perplexity(\n\u001b[0;32m---> 95\u001b[0;31m         distances, neighbors, desired_perplexity, verbose)\n\u001b[0m\u001b[1;32m     96\u001b[0m     \u001b[0mm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"All probabilities should be finite\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     97\u001b[0m     \u001b[0;32massert\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconditional_P\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32msklearn/manifold/_utils.pyx\u001b[0m in \u001b[0;36msklearn.manifold._utils._binary_search_perplexity (sklearn/manifold/_utils.c:2677)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32msklearn/manifold/_utils.pyx\u001b[0m in \u001b[0;36msklearn.manifold._utils._binary_search_perplexity (sklearn/manifold/_utils.c:2415)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "#tSNE\n",
    "print('----Start tSNE----')\n",
    "for per in [300,400,500]:\n",
    "    print('-'*10+'tSNE on perplexity %i'%per + '-'*10)\n",
    "    X_red, y_red = tSNE(X, y, per)\n",
    "    plt.figure()\n",
    "    y_color, y_s = some_cluster(y_red, colors)\n",
    "    #color\n",
    "    plt.scatter(X_red[:, 0], X_red[:, 1], c=y_color, s=y_s, marker='o',\n",
    "                linewidths=0)\n",
    "    plt.title(title)\n",
    "    plt.savefig(IM_FOLDER+'tSNE%i.png' % per)\n",
    "    #black-white\n",
    "    plt.scatter(X_red[:, 0], X_red[:, 1], c='k', s=y_s, marker='o',\n",
    "                linewidths=0)\n",
    "    plt.title(title)\n",
    "    plt.savefig(IM_FOLDER+'tSNE_bw%i.png'%per)\n",
    "    #\n",
    "    plt.figure()\n",
    "    f3, ax3 = plt.subplots(1, D)\n",
    "    for d in range(D):\n",
    "        ax3[d].scatter(X_red[:, 0], X_red[:, 1],c =y_red[:,d] ,marker='*', linewidths=0)\n",
    "        ax3[d].set_title('Mode%3i'%d)\n",
    "    plt.setp([a.get_xticklabels() for a in ax3], visible=False)\n",
    "    plt.setp([a.get_yticklabels() for a in ax3], visible=False)\n",
    "    plt.savefig(IM_FOLDER+'tsne_table.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": 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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}