04. Dimension Reduction and Images.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sandbox calcs from lecture!\n",
    "# import numpy as np\n",
    "# k = np.array([[37.1,-6.55],\n",
    "#               [-6.55,43.99]])\n",
    "# eigenvalues, eigenvectors = np.linalg.eig(k)\n",
    "# print (eigenvectors, eigenvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Dimensionality Reduction with a Simple Dataset\n",
    "Here we are going to start with the iris dataset. This is a very simple dataset that is four dimensions. Details about this data are rampant on the web. \n",
    "- Here is a good explanation from UCI: https://archive.ics.uci.edu/ml/datasets/Iris \n",
    "- And even more details from Wikipedia: http://en.wikipedia.org/wiki/Iris_flower_data_set \n",
    "\n",
    "Let's get started. We will use scikit for this analysis, which has a built-in loading mechanism for the iris dataset.\n",
    "\n",
    "## Dependencies\n",
    "- `jupyter labextension install @jupyter-widgets/jupyterlab-manager`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets as ds\n",
    "\n",
    "iris = ds.load_iris() # wow that was simple to load up!\n",
    "\n",
    "print (iris['DESCR']) # a description of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's convert it to a pandas data frame and visualize it\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df = pd.DataFrame(iris.data,columns=iris.feature_names)\n",
    "# add in the class targets and names\n",
    "df['target'] = iris.target.astype(np.int)\n",
    "print (df.info())\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set()\n",
    "sns.pairplot(df, hue=\"target\", size=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using PCA with sklearn\n",
    "Now let's use PCA to find the two \"best\" dimensions of this data\n",
    "These are linear transforms to help project the features into something more understandable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "target_names = iris.target_names\n",
    "\n",
    "pca = PCA(n_components=2)\n",
    "X_pca = pca.fit(X).transform(X) # fit data and then transform it\n",
    "\n",
    "# print the components\n",
    "\n",
    "print ('pca:', pca.components_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "cmap = sns.set(style=\"darkgrid\") \n",
    "\n",
    "# this function definition just formats the weights into readable strings\n",
    "# you can skip it without loss of generality to the Data Science content\n",
    "def get_feature_names_from_weights(weights, names):\n",
    "    tmp_array = []\n",
    "    for comp in weights:\n",
    "        tmp_string = ''\n",
    "        for fidx,f in enumerate(names):\n",
    "            if fidx>0 and comp[fidx]>=0:\n",
    "                tmp_string+='+'\n",
    "            tmp_string += '%.2f*%s ' % (comp[fidx],f[:-5])\n",
    "        tmp_array.append(tmp_string)\n",
    "    return tmp_array\n",
    "  \n",
    "plt.style.use('ggplot')\n",
    "# now let's get to the Data Analytics!\n",
    "pca_weight_strings = get_feature_names_from_weights(pca.components_, iris.feature_names) \n",
    "\n",
    "# create some pandas dataframes from the transformed outputs\n",
    "df_pca = pd.DataFrame(X_pca,columns=[pca_weight_strings])\n",
    "\n",
    "from matplotlib.pyplot import scatter\n",
    "\n",
    "# scatter plot the output, with the names created from the weights\n",
    "ax = scatter(X_pca[:,0], X_pca[:,1], c=y, s=(y+2)*10, cmap=cmap)\n",
    "plt.xlabel(pca_weight_strings[0]) \n",
    "plt.ylabel(pca_weight_strings[1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Biplots \n",
    "Biplots are visualizations that allow you to see how much of each principle component is represented by each feature. They are used many times in survey data results to display the responses of multiple survey questions. An example biplots used in this fashion is available here: https://datasciencebe.files.wordpress.com/2015/03/pca-skillscores.png\n",
    "\n",
    "![biplot](https://datasciencebe.files.wordpress.com/2015/03/pca-skillscores.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Manipulated example from https://github.com/teddyroland/python-biplot/blob/master/biplot.py\n",
    "\n",
    "def biplot(pca, dat, title=''):\n",
    "    \n",
    "    import plotly\n",
    "    from plotly.graph_objs import Bar, Line\n",
    "    from plotly.graph_objs import Scatter, Layout\n",
    "    from plotly.graph_objs.scatter import Marker\n",
    "    from plotly.graph_objs.layout import XAxis, YAxis\n",
    "    plotly.offline.init_notebook_mode() # run at the start of every notebook\n",
    "    \n",
    "    # 0,1 denote PC1 and PC2; change values for other PCs\n",
    "    xvector = pca.components_[0] \n",
    "    yvector = pca.components_[1]\n",
    "\n",
    "    tmp = pca.transform(dat.values)\n",
    "    xs = tmp[:,0] \n",
    "    ys = tmp[:,1]\n",
    "\n",
    "    annotations = [Scatter(x=xs, y=ys, mode ='markers',\n",
    "                           marker=Marker(color=iris.target),\n",
    "                           name='PCA Trans. Data')]\n",
    "    for i in range(len(xvector)):\n",
    "        txt = list(dat.columns.values)[i]\n",
    "        annotations.append(\n",
    "                Scatter( \n",
    "                    x=[0, xvector[i]*max(xs)],\n",
    "                    y=[0, yvector[i]*max(ys)],\n",
    "                    mode='lines+text',\n",
    "                    text=['', txt],\n",
    "                    name=txt,\n",
    "                ))\n",
    "    \n",
    "    plotly.offline.iplot({\n",
    "        \"data\": annotations,\n",
    "        \"layout\": Layout(xaxis=XAxis(title='Principal Component One'), \n",
    "                         yaxis=YAxis(title='Principal Component Two'),\n",
    "                        title=title)\n",
    "    })\n",
    "\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "X = iris.data\n",
    "pca = PCA(n_components=4)\n",
    "pca.fit(X) \n",
    "biplot(pca,pd.DataFrame(iris.data,columns=iris.feature_names),'Iris Biplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explained Variance\n",
    "Recall that the explained variance is a factor of the eigen values. We can access these via the following equation:\n",
    "\n",
    "$$ r_q=\\frac{\\sum_{j=1}^q \\lambda_j}{\\sum_{j=1}^p \\lambda_j} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# manipulated from Sebastian Raschka Example (your book!)\n",
    "# also from hi blog here: http://sebastianraschka.com/Articles/2015_pca_in_3_steps.html\n",
    "\n",
    "\n",
    "def plot_explained_variance(pca):\n",
    "    import plotly\n",
    "    from plotly.graph_objs import Bar, Line\n",
    "    from plotly.graph_objs import Scatter, Layout\n",
    "    from plotly.graph_objs.scatter import Marker\n",
    "    from plotly.graph_objs.layout import XAxis, YAxis\n",
    "    plotly.offline.init_notebook_mode() # run at the start of every notebook\n",
    "    \n",
    "    explained_var = pca.explained_variance_ratio_\n",
    "    cum_var_exp = np.cumsum(explained_var)\n",
    "    \n",
    "    plotly.offline.iplot({\n",
    "        \"data\": [Bar(y=explained_var, name='individual explained variance'),\n",
    "                 Scatter(y=cum_var_exp, name='cumulative explained variance')\n",
    "            ],\n",
    "        \"layout\": Layout(xaxis=XAxis(title='Principal components'), yaxis=YAxis(title='Explained variance ratio'))\n",
    "    })\n",
    "        \n",
    "\n",
    "pca = PCA(n_components=4)\n",
    "X_pca = pca.fit(X)\n",
    "plot_explained_variance(pca)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PCA With Image Data\n",
    "This is code manipulated from Olivier Grisel's eigen face classification demonstration. You can find the original notebook here: http://nbviewer.ipython.org/github/ogrisel/notebooks/blob/master/Labeled%20Faces%20in%20the%20Wild%20recognition.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fetch the images for the dataset\n",
    "# this will take a long time the first run because it needs to download\n",
    "# after the first time, the dataset will be save to your disk (in sklearn package somewhere) \n",
    "# if this does not run, you may need additional libraries installed on your system (install at your own risk!!)\n",
    "from sklearn.datasets import fetch_lfw_people\n",
    "\n",
    "lfw_people = fetch_lfw_people(min_faces_per_person=20, resize=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get some of the specifics of the dataset\n",
    "X = lfw_people.data\n",
    "y = lfw_people.target\n",
    "names = lfw_people.target_names\n",
    "\n",
    "n_samples, n_features = X.shape\n",
    "_, h, w = lfw_people.images.shape\n",
    "n_classes = len(names)\n",
    "\n",
    "print(np.sum(~np.isfinite(X)))\n",
    "print(\"n_samples: {}\".format(n_samples))\n",
    "print(\"n_features: {}\".format(n_features))\n",
    "print(\"n_classes: {}\".format(n_classes))\n",
    "print(\"Original Image Sizes {} by {}\".format(h,w))\n",
    "print (125*94) # the size of the images are the size of the feature vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So basically each feature vector is a giant image with the rows of the image just stacked one after the other into a giant vector. The image sizes are 125 pixels by 94 pixels. This gives us 125x94=11750 pixels per image.\n",
    "\n",
    "So we are using each pixel location in the image as a separate feature. Wow. That's a lot of features. The process is as follows:\n",
    "\n",
    "<img src=\"https://i.ytimg.com/vi/_lY74pXWlS8/maxresdefault.jpg\" alt=\"Drawing\" width=\"500\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a helper plotting function\n",
    "def plot_gallery(images, titles, h, w, n_row=3, n_col=6):\n",
    "    \"\"\"Helper function to plot a gallery of portraits\"\"\"\n",
    "    plt.figure(figsize=(1.7 * n_col, 2.3 * n_row))\n",
    "    plt.subplots_adjust(bottom=0, left=.01, right=.99, top=.90, hspace=.35)\n",
    "    for i in range(n_row * n_col):\n",
    "        plt.subplot(n_row, n_col, i + 1)\n",
    "        plt.imshow(images[i].reshape((h, w)), cmap=plt.cm.gray)\n",
    "        plt.title(titles[i], size=12)\n",
    "        plt.xticks(())\n",
    "        plt.yticks(())\n",
    "\n",
    "plot_gallery(X, names[y], h, w) # defaults to showing a 3 by 6 subset of the faces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reconstruct Using Full PCA\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets do some PCA of the features and go from 1850 features to 20 features\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "n_components = 300\n",
    "print (\"Extracting the top %d eigenfaces from %d faces\" % (\n",
    "    n_components, X.shape[0]))\n",
    "\n",
    "pca = PCA(n_components=n_components)\n",
    "%time pca.fit(X.copy())\n",
    "eigenfaces = pca.components_.reshape((n_components, h, w))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_explained_variance(pca)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenface_titles = [\"eigenface %d\" % i for i in range(eigenfaces.shape[0])]\n",
    "plot_gallery(eigenfaces, eigenface_titles, h, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstruct_image(trans_obj,org_features):\n",
    "    low_rep = trans_obj.transform(org_features)\n",
    "    rec_image = trans_obj.inverse_transform(low_rep)\n",
    "    return low_rep, rec_image\n",
    "    \n",
    "idx_to_reconstruct = 4    \n",
    "X_idx = X[idx_to_reconstruct]\n",
    "low_dimensional_representation, reconstructed_image = reconstruct_image(pca,X_idx.reshape(1, -1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(1,2,1)\n",
    "plt.imshow(X_idx.reshape((h, w)), cmap=plt.cm.gray)\n",
    "plt.title('Original')\n",
    "plt.grid()\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(reconstructed_image.reshape((h, w)), cmap=plt.cm.gray)\n",
    "plt.title('Reconstructed from Full PCA')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now Let's Use Randomized PCA\n",
    "Of course, the methods that are used to calculate eigenvectors (the components) do not scale well with the number of features (in this case, the image dimensions). Recall that PCA works on the covariance matrix and the covariance matrix size is not affected by the number of examples, only the number of features. Randomized PCA helps to mitigate this by formulating the eigenvectors of a lower rank matrix using randomized projections. \n",
    "\n",
    "More information here: \n",
    "- http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html\n",
    "- ftp://ftp.math.ucla.edu/pub/camreport/cam08-60.pdf\n",
    "- Halko, et al., 2009 (arXiv:909) A randomized algorithm for the decomposition of matrices Per-Gunnar Martinsson, Vladimir Rokhlin and Mark Tygert"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets do some PCA of the features and go from 1850 features to 300 features\n",
    "\n",
    "n_components = 300\n",
    "print (\"Extracting the top %d eigenfaces from %d faces\" % (\n",
    "    n_components, X.shape[0]))\n",
    "\n",
    "rpca = PCA(n_components=n_components,svd_solver='randomized')\n",
    "%time rpca.fit(X.copy())\n",
    "eigenfaces = rpca.components_.reshape((n_components, h, w))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenface_titles = [\"eigenface %d\" % i for i in range(eigenfaces.shape[0])]\n",
    "plot_gallery(eigenfaces, eigenface_titles, h, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernel PCA \n",
    "Warning: the block of code below takes a long time to run. Kernel methods typically take some time, so I would advise subsampling the data before using the kernel. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.decomposition import KernelPCA\n",
    "\n",
    "n_components = 300\n",
    "print (\"Extracting the top %d eigenfaces from %d faces, not calculating inverse transform\" % (n_components, X.shape[0]))\n",
    "\n",
    "kpca = KernelPCA(n_components=n_components, kernel='rbf', \n",
    "                fit_inverse_transform=False, gamma=12, # very sensitive to the gamma parameter,\n",
    "                remove_zero_eig=True)  \n",
    "kpca.fit(X.copy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "#  THIS  TAKES A LONG TIME TO RUN\n",
    "from sklearn.decomposition import KernelPCA\n",
    "\n",
    "n_components = 300\n",
    "print (\"Extracting the top %d eigenfaces from %d faces, ALSO getting inverse transform\" % (n_components, X.shape[0]))\n",
    "\n",
    "kpca = KernelPCA(n_components=n_components, kernel='rbf', \n",
    "                fit_inverse_transform=True, gamma=12, # very sensitive to the gamma parameter,\n",
    "                remove_zero_eig=True)  \n",
    "kpca.fit(X.copy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the above operation takes a long time to save the inverse transform parameters\n",
    "# so let's save out the results to load in later!\n",
    "import pickle\n",
    "\n",
    "pickle.dump(kpca, open( 'large_data/kpca.p', 'wb' ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load the Kpca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "kpca = pickle.load(open( 'large_data/kpca.p', 'rb' ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# widgets example\n",
    "from ipywidgets import widgets  # make this interactive!\n",
    "\n",
    "widgets.interact(lambda x: print(x),x=(0,5,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "# warnings.simplefilter('ignore', DeprecationWarning)\n",
    "# warnings.simplefilter(\"always\",DeprecationWarning)\n",
    "\n",
    "\n",
    "\n",
    "def plt_reconstruct(idx_to_reconstruct):\n",
    "    idx_to_reconstruct = np.round(idx_to_reconstruct)\n",
    "    \n",
    "    reconstructed_image = pca.inverse_transform(pca.transform(X[idx_to_reconstruct].reshape(1, -1)))\n",
    "    reconstructed_image_rpca = rpca.inverse_transform(rpca.transform(X[idx_to_reconstruct].reshape(1, -1)))\n",
    "    reconstructed_image_kpca = kpca.inverse_transform(kpca.transform(X[idx_to_reconstruct].reshape(1, -1)))\n",
    "    \n",
    "    \n",
    "    plt.figure(figsize=(15,7))\n",
    "    \n",
    "    plt.subplot(1,4,1)\n",
    "    plt.imshow(X[idx_to_reconstruct].reshape((h, w)), cmap=plt.cm.gray)\n",
    "    plt.title(names[y[idx_to_reconstruct]])\n",
    "    plt.grid()\n",
    "    \n",
    "    plt.subplot(1,4,2)\n",
    "    plt.imshow(reconstructed_image.reshape((h, w)), cmap=plt.cm.gray)\n",
    "    plt.title('Full PCA')\n",
    "    plt.grid()\n",
    "    \n",
    "    plt.subplot(1,4,3)\n",
    "    plt.imshow(reconstructed_image_rpca.reshape((h, w)), cmap=plt.cm.gray)\n",
    "    plt.title('Randomized PCA')\n",
    "    plt.grid()\n",
    "    \n",
    "    plt.subplot(1,4,4)\n",
    "    plt.imshow(reconstructed_image_kpca.reshape((h, w)), cmap=plt.cm.gray)\n",
    "    plt.title('Kernel PCA')\n",
    "    plt.grid()\n",
    "    \n",
    "    \n",
    "widgets.interact(plt_reconstruct,idx_to_reconstruct=(0,n_samples-1,1),__manual=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "# Image Processing Crash Course\n",
    "Now let's look at other Methods for extracting features from images.\n",
    "\n",
    "Let's start by calculating gradients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.io import imshow\n",
    "\n",
    "idx_to_reconstruct = int(np.random.rand(1)*len(X))\n",
    "img  = X[idx_to_reconstruct].reshape((h,w))\n",
    "imshow(img)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.filters import sobel_h, sobel_v\n",
    "\n",
    "gradient_mag = np.sqrt(sobel_v(img)**2 + sobel_h(img)**2 ) \n",
    "imshow(gradient_mag)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Feature: DAISY Bag of Features Model\n",
    "As discussed in lecture, DAISY Features are a means of looking at histograms of edges weights using gradients at various orientations (if you are new to this, that is not a satisfying explanation, please look back at the notes/slides for the course.). \n",
    "\n",
    "http://scikit-image.org/docs/dev/_images/sphx_glr_plot_daisy_001.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.feature import daisy\n",
    "\n",
    "# lets first visualize what the daisy descripto looks like\n",
    "features, img_desc = daisy(img,step=40, radius=10, rings=3, histograms=5, orientations=8, visualize=True)\n",
    "imshow(img_desc)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's understand how to use it\n",
    "features = daisy(img,step=10, radius=10, rings=2, histograms=4, orientations=8, visualize=False)\n",
    "print(features.shape)\n",
    "print(features.shape[0]*features.shape[1]*features.shape[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a function to tak in the row of the matric and return a new feature\n",
    "def apply_daisy(row,shape):\n",
    "    feat = daisy(row.reshape(shape),step=10, radius=10, rings=2, histograms=6, orientations=8, visualize=False)\n",
    "    return feat.reshape((-1))\n",
    "\n",
    "%time test_feature = apply_daisy(X[3],(h,w))\n",
    "test_feature.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "0.027 * len(X) # approximate how long it may run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply to entire data, row by row,\n",
    "# takes about a minute to run\n",
    "%time daisy_features = np.apply_along_axis(apply_daisy, 1, X, (h,w))\n",
    "print(daisy_features.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "# find the pairwise distance between all the different image features\n",
    "%time dist_matrix = pairwise_distances(daisy_features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "# find closest image to current image\n",
    "idx1 = 5\n",
    "distances = copy.deepcopy(dist_matrix[idx1,:])\n",
    "distances[idx1] = np.infty # dont pick the same image!\n",
    "idx2 = np.argmin(distances)\n",
    "\n",
    "plt.figure(figsize=(7,10))\n",
    "plt.subplot(1,2,1)\n",
    "imshow(X[idx1].reshape((h,w)))\n",
    "plt.title(\"Original Image\")\n",
    "plt.grid()\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "imshow(X[idx2].reshape((h,w)))\n",
    "plt.title(\"Closest Image\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import fixed\n",
    "# put it together inside a nice widget\n",
    "def closest_image(dmat,idx1):\n",
    "    distances = copy.deepcopy(dmat[idx1,:]) # get all image diatances\n",
    "    distances[idx1] = np.infty # dont pick the same image!\n",
    "    idx2 = np.argmin(distances)\n",
    "    \n",
    "    distances[idx2] = np.infty\n",
    "    idx3 = np.argmin(distances)\n",
    "    \n",
    "    plt.figure(figsize=(10,16))\n",
    "    plt.subplot(1,3,1)\n",
    "    imshow(X[idx1].reshape((h,w)))\n",
    "    plt.title(\"Original Image \"+names[y[idx1]])\n",
    "    plt.grid()\n",
    "\n",
    "    plt.subplot(1,3,2)\n",
    "    imshow(X[idx2].reshape((h,w)))\n",
    "    plt.title(\"Closest Image  \"+names[y[idx2]])\n",
    "    plt.grid()\n",
    "    \n",
    "    plt.subplot(1,3,3)\n",
    "    imshow(X[idx3].reshape((h,w)))\n",
    "    plt.title(\"Next Closest Image \"+names[y[idx3]])\n",
    "    plt.grid()\n",
    "    \n",
    "widgets.interact(closest_image,idx1=(0,n_samples-1,1),dmat=fixed(dist_matrix),__manual=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gabor Kernels for Feature Extraction\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from skimage.filters import gabor_kernel\n",
    "from scipy import ndimage as ndi\n",
    "from scipy import stats\n",
    "\n",
    "# prepare filter bank kernels\n",
    "kernels = []\n",
    "for theta in range(4):\n",
    "    theta = theta / 4. * np.pi\n",
    "    for sigma in (1, 3):\n",
    "        for frequency in (0.05, 0.25):\n",
    "            kernel = np.real(gabor_kernel(frequency, theta=theta,\n",
    "                                          sigma_x=sigma, sigma_y=sigma))\n",
    "            kernels.append(kernel)\n",
    "\n",
    "            \n",
    "# compute the filter bank and take statistics of image\n",
    "def compute_gabor(row, kernels, shape):\n",
    "    feats = np.zeros((len(kernels), 4), dtype=np.double)\n",
    "    for k, kernel in enumerate(kernels):\n",
    "        filtered = ndi.convolve(row.reshape(shape), kernel, mode='wrap')\n",
    "        _,_,feats[k,0],feats[k,1],feats[k,2],feats[k,3] = stats.describe(filtered.reshape(-1))\n",
    "        # mean, var, skew, kurt\n",
    "        \n",
    "    return feats.reshape(-1)\n",
    "\n",
    "idx_to_reconstruct = int(np.random.rand(1)*len(X))\n",
    "\n",
    "gabr_feature = compute_gabor(X[idx_to_reconstruct], kernels, (h,w))\n",
    "gabr_feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# takes ~3 minutes to run entire dataset\n",
    "%time gabor_stats = np.apply_along_axis(compute_gabor, 1, X, kernels, (h,w))\n",
    "print(gabor_stats.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "# find the pairwise distance between all the different image features\n",
    "%time dist_matrix_gabor = pairwise_distances(gabor_stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "widgets.interact(closest_image,idx1=(0,n_samples-1,1),dmat=fixed(dist_matrix_gabor),__manual=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Is gabor filtering a good feature extraction method for this dataset? **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}