index.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Machine Learning with Scikit-Learn: Feature Engineering\n",
    "## Jake VanderPlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The code in these cells is runable. \n",
    "# Click on this cell, then press Shift+Enter to run it, \n",
    "# or click the Run button in the toolbar.\n",
    "\n",
    "print(\"Hello, World!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In previous notebooks, we discussed the fact that machine learning requires datasets to be represented by an n_samples x n_features matrix, but it is not always clear how real-world data can be expressed this way. In this notebook, we will discuss a few specific examples of engineering features for use in machine learning algorithms, including:\n",
    "\n",
    "- One-hot encoding for categorical data\n",
    "- Frequency-based encoding for textual data\n",
    "- Histogram-of-gradient (HOG) features for image data\n",
    "\n",
    "In each of these, I will give a brief example of applying these features to a dataset in the course of a machine learning algorithm.\n",
    "\n",
    "At the end of this notebook, you will:\n",
    "- Understand several common approaches to feature engineering\n",
    "- Gain exposure to more examples of Scikit-learn’s API applied to real-world datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll start with the standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "plt.style.use('seaborn')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Categorical Features\n",
    "\n",
    "One common type of non-numerical data is *categorical* data.\n",
    "For example, imagine you are exploring some data on housing prices, and along with numerical features like \"price\" and \"rooms\", you also have \"neighborhood\" information.\n",
    "For example, your data might look something like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = [\n",
    "    {'price': 850000, 'rooms': 4, 'neighborhood': 'Queen Anne'},\n",
    "    {'price': 700000, 'rooms': 3, 'neighborhood': 'Fremont'},\n",
    "    {'price': 650000, 'rooms': 3, 'neighborhood': 'Wallingford'},\n",
    "    {'price': 600000, 'rooms': 2, 'neighborhood': 'Fremont'}\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might be tempted to encode this data with a straightforward numerical mapping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "{'Queen Anne': 1, 'Fremont': 2, 'Wallingford': 3};"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that this is not generally a useful approach in scikit-learn: the package's models make the fundamental assumption that numerical features reflect algebraic quantities.\n",
    "Thus, such a mapping would imply, for example, that *Queen Anne < Fremont < Wallingford*, or even that *Wallingford - Queen Anne = Fremont*, which (niche demographic jokes aside) does not make much sense.\n",
    "\n",
    "In this case, one proven technique is to use *one-hot encoding*, which effectively creates extra columns indicating the presence or absence of a category with a value of 1 or 0, respectively.\n",
    "When your data comes as a list of dictionaries, scikit-learn's ``DictVectorizer`` will do this for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction import DictVectorizer\n",
    "vec = DictVectorizer(sparse=False, dtype=int)\n",
    "vec.fit_transform(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the neighborhood column has been expanded into three separate columns, representing the three neighborhood labels, and that each row has a 1 in the column associated with its neighborhood.\n",
    "With these categorical features thus encoded, you can proceed as normal with fitting a scikit-learn model.\n",
    "\n",
    "To see the meaning of each column, you can inspect the feature names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vec.get_feature_names()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.DataFrame(vec.fit_transform(data),\n",
    "             columns=vec.get_feature_names())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is one clear disadvantage of this approach: if your category has many possible values, this can *greatly* increase the size of your dataset.\n",
    "However, because the encoded data contains mostly zeros, a sparse output can be a very efficient solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vec = DictVectorizer(sparse=True, dtype=int)\n",
    "vec.fit_transform(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many (though not yet all) of the scikit-learn estimators accept such sparse inputs when fitting and evaluating models. ``sklearn.preprocessing.OneHotEncoder`` and ``sklearn.feature_extraction.FeatureHasher`` are two additional tools that scikit-learn includes to support this type of encoding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Text Features\n",
    "\n",
    "Another common need in feature engineering is to convert text to a set of representative numerical values.\n",
    "For example, most automatic mining of social media data relies on some form of encoding the text as numbers.\n",
    "One of the simplest methods of encoding data is by *word counts*: you take each snippet of text, count the occurrences of each word within it, and put the results in a table.\n",
    "\n",
    "For example, consider the following set of three phrases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sample = ['problem of evil',\n",
    "          'evil queen',\n",
    "          'horizon problem']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a vectorization of this data based on word count, we could construct a column representing the word \"problem,\" the word \"evil,\" the word \"horizon,\" and so on.\n",
    "While doing this by hand would be possible, the tedium can be avoided by using scikit-learn's ``CountVectorizer``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "\n",
    "vec = CountVectorizer()\n",
    "X = vec.fit_transform(sample)\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a sparse matrix recording the number of times each word appears; it is easier to inspect if we convert this to a ``DataFrame`` with labeled columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.DataFrame(X.toarray(), columns=vec.get_feature_names())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some issues with this approach, however: the raw word counts lead to features that put too much weight on words that appear very frequently, and this can be suboptimal in some classification algorithms.\n",
    "One approach to fix this is known as *term frequency-inverse document frequency* (*TF–IDF*), which weights the word counts by a measure of how often they appear in the documents.\n",
    "The syntax for computing these features is similar to the previous example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "vec = TfidfVectorizer()\n",
    "X = vec.fit_transform(sample)\n",
    "pd.DataFrame(X.toarray(), columns=vec.get_feature_names())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Text Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at classifying some text based on this approach.\n",
    "\n",
    "We'll use the *20 newsgroups* dataset, available from the web via scikit-learn's dataset tools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_20newsgroups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "categories = ['talk.religion.misc', 'soc.religion.christian',\n",
    "              'sci.space', 'comp.graphics']\n",
    "train = fetch_20newsgroups(subset='train', categories=categories)\n",
    "test = fetch_20newsgroups(subset='test', categories=categories)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data consist of the full text content of messages posted to early internet newsgroups. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(train.data[5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The target is a set of numbers encoding which newsgroup the messages belong to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train.target[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train.target_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the ``TfidfVectorizer`` and put it in a pipeline with a ``MultinomialNB`` classifier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import TfidfVectorizer\n",
    "from sklearn.naive_bayes import MultinomialNB\n",
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "model = make_pipeline(TfidfVectorizer(), MultinomialNB())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now fit the model to the training data, and predict the labels of the test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.fit(train.data, train.target)\n",
    "labels = model.predict(test.data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's evaluate this with a confusion matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "import seaborn as sns\n",
    "\n",
    "mat = confusion_matrix(test.target, labels)\n",
    "sns.heatmap(mat.T, square=True, annot=True, fmt='d', cbar=False, cmap='Blues',\n",
    "            xticklabels=train.target_names, yticklabels=train.target_names)\n",
    "plt.xlabel('true label')\n",
    "plt.ylabel('predicted label');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More interestingly, we can now use this classifier to predict, given some phrase, what general topic the phrase might be referring to.\n",
    "\n",
    "We'll use the following quick convenience function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def predict_category(s, train=train, model=model):\n",
    "    pred = model.predict([s])\n",
    "    return train.target_names[pred[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predict_category('sending a payload to the ISS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predict_category('discussing islam vs atheism')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "predict_category('determining the screen resolution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that this is nothing more sophisticated than a simple probability model for the (weighted) frequency of each word in the string; nevertheless, the result is striking.\n",
    "Even a very naive algorithm, when used carefully and trained on a large set of high-dimensional data, can be surprisingly effective."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Features\n",
    "\n",
    "We've seen an approach to image classification that relies on the pixels themselves as features.\n",
    "The problem is that the relationships between pixels are often more meaningful: for example, lines and gradients tell us much more about the content of the image, and are much more robust to lighting levels and other incidental differences between images.\n",
    "\n",
    "For more complicated tasks involving images, it can be useful to derive other features that are more informative.\n",
    "There are many approaches to this, but the one we'll demo here is known as *histograms of gradients* (HOG).\n",
    "HOG involves the following steps:\n",
    "\n",
    "1. Optionally prenormalize images. This leads to features that resist dependence on variations in illumination.\n",
    "2. Convolve the image with two filters that are sensitive to horizontal and vertical brightness gradients. These capture edge, contour, and texture information.\n",
    "3. Subdivide the image into cells of a predetermined size, and compute a histogram of the gradient orientations within each cell.\n",
    "4. Normalize the histograms in each cell by comparing to the block of neighboring cells. This further suppresses the effect of illumination across the image.\n",
    "5. Construct a one-dimensional feature vector from the information in each cell.\n",
    "\n",
    "A fast HOG extractor is built into the scikit-image project, and we can try it out relatively quickly and visualize the oriented gradients within each cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from skimage import data, color, feature\n",
    "import skimage.data\n",
    "\n",
    "image = color.rgb2gray(data.chelsea())\n",
    "hog_vec, hog_vis = feature.hog(image, visualise=True)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(12, 6),\n",
    "                       subplot_kw=dict(xticks=[], yticks=[]))\n",
    "ax[0].imshow(image, cmap='gray')\n",
    "ax[0].set_title('input image')\n",
    "\n",
    "ax[1].imshow(hog_vis, cmap='binary')\n",
    "ax[1].set_title('visualization of HOG features');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HOG in Action: A Simple Face Detector\n",
    "\n",
    "Using these HOG features, we can build up a simple facial detection algorithm with any scikit-learn estimator; here we will use a linear support vector machine (refer back to _Machine Learning with Scikit-Learn: Support Vector Machines_ if you need a refresher on this).\n",
    "The steps are as follows:\n",
    "\n",
    "1. Obtain a set of image thumbnails of faces to constitute \"positive\" training samples.\n",
    "2. Obtain a set of image thumbnails of non-faces to constitute \"negative\" training samples.\n",
    "3. Extract HOG features from these training samples.\n",
    "4. Train a linear SVM classifier on these samples.\n",
    "5. For an \"unknown\" image, pass a sliding window across the image, using the model to evaluate whether that window contains a face or not.\n",
    "6. If detections overlap, combine them into a single window.\n",
    "\n",
    "Let's go through these steps and try it out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Obtain a Set of Positive Training Samples\n",
    "\n",
    "Let's start by finding some positive training samples that show a variety of faces.\n",
    "We have one easy set of data to work with—-the Labeled Faces in the Wild dataset, which can be downloaded by scikit-learn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_lfw_people\n",
    "faces = fetch_lfw_people()\n",
    "positive_patches = faces.images\n",
    "positive_patches.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us a sample of 13,000 face images to use for training."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Obtain a Set of Negative Training Samples\n",
    "\n",
    "Next, we need a set of similarly sized thumbnails that *do not* have a face in them.\n",
    "One way to do this is to take any corpus of input images and extract thumbnails from them at a variety of scales.\n",
    "Here we can use some of the images shipped with scikit-image, along with scikit-learn's ``PatchExtractor``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from skimage import data, transform\n",
    "\n",
    "imgs_to_use = ['camera', 'text', 'coins', 'moon',\n",
    "               'page', 'clock', 'immunohistochemistry',\n",
    "               'chelsea', 'coffee', 'hubble_deep_field']\n",
    "images = [color.rgb2gray(getattr(data, name)())\n",
    "          for name in imgs_to_use]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.image import PatchExtractor\n",
    "\n",
    "def extract_patches(img, N, scale=1.0, patch_size=positive_patches[0].shape):\n",
    "    extracted_patch_size = tuple((scale * np.array(patch_size)).astype(int))\n",
    "    extractor = PatchExtractor(patch_size=extracted_patch_size,\n",
    "                               max_patches=N, random_state=0)\n",
    "    patches = extractor.transform(img[np.newaxis])\n",
    "    if scale != 1:\n",
    "        patches = np.array([transform.resize(patch, patch_size, mode='constant')\n",
    "                            for patch in patches])\n",
    "    return patches\n",
    "\n",
    "negative_patches = np.vstack([extract_patches(im, 1000, scale)\n",
    "                              for im in images for scale in [0.5, 1.0, 2.0]])\n",
    "negative_patches.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have 30,000 suitable image patches that do not contain faces.\n",
    "Let's take a look at a few of them to get an idea of what they look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(6, 10)\n",
    "for i, axi in enumerate(ax.flat):\n",
    "    axi.imshow(negative_patches[500 * i], cmap='gray')\n",
    "    axi.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our hope is that these would sufficiently cover the space of non-faces that our algorithm is likely to see."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Combine Sets and Extract HOG Features\n",
    "\n",
    "Now that we have these positive samples and negative samples, we can combine them and compute HOG features.\n",
    "This step takes a little while, because the HOG features involve a nontrivial computation for each image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from itertools import chain\n",
    "X_train = np.array([feature.hog(im)\n",
    "                    for im in chain(positive_patches,\n",
    "                                    negative_patches)])\n",
    "y_train = np.zeros(X_train.shape[0])\n",
    "y_train[:positive_patches.shape[0]] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are left with 43,000 training samples in 1,215 dimensions, and we now have our data in a form that we can feed into scikit-learn!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Training a Support Vector Machine\n",
    "\n",
    "Next, we use the tools we have been exploring in this Oriole to create a classifier of thumbnail patches.\n",
    "For such a high-dimensional binary classification task, a linear support vector machine is a good choice.\n",
    "We will use scikit-learn's ``LinearSVC``, because in comparison to ``SVC``, it often has better scaling for a large number of samples.\n",
    "\n",
    "First, though, let's use a simple Gaussian naive Bayes classifier to get a quick baseline:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "cross_val_score(GaussianNB(), X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that on our training data, even a simple naive Bayes algorithm gets us upward of 90% accuracy.\n",
    "Let's try the support vector machine, with a grid search over a few choices of the C parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.svm import LinearSVC\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "grid = GridSearchCV(LinearSVC(), {'C': [1.0, 2.0, 4.0, 8.0]})\n",
    "grid.fit(X_train, y_train)\n",
    "grid.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take the best estimator and retrain it on the full dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = grid.best_estimator_\n",
    "model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ``fit`` step finds our best estimator for detecting the presence of a face in a small image thumbnail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Find Faces in a New Image\n",
    "\n",
    "Now that we have this model in place, let's grab a new image and see how the model does.\n",
    "We will use one portion of the astronaut image for simplicity (see a discussion of this below in \"Caveats and Improvements,\" and run a sliding window over it and evaluate each patch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test_image = skimage.data.astronaut()\n",
    "test_image = skimage.color.rgb2gray(test_image)\n",
    "test_image = skimage.transform.rescale(test_image, 0.5, mode='constant')\n",
    "test_image = test_image[:160, 40:180]\n",
    "\n",
    "plt.imshow(test_image, cmap='gray')\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's create a window that iterates over patches of this image, and compute HOG features for each patch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sliding_window(img, patch_size=positive_patches[0].shape,\n",
    "                   istep=2, jstep=2, scale=1.0):\n",
    "    Ni, Nj = (int(scale * s) for s in patch_size)\n",
    "    for i in range(0, img.shape[0] - Ni, istep):\n",
    "        for j in range(0, img.shape[1] - Ni, jstep):\n",
    "            patch = img[i:i + Ni, j:j + Nj]\n",
    "            if scale != 1:\n",
    "                patch = transform.resize(patch, patch_size)\n",
    "            yield (i, j), patch\n",
    "            \n",
    "indices, patches = zip(*sliding_window(test_image))\n",
    "patches_hog = np.array([feature.hog(patch)\n",
    "                        for patch in patches])\n",
    "patches_hog.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can take these HOG-featured patches and use our model to evaluate whether each patch contains a face:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "labels = model.predict(patches_hog)\n",
    "labels.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that out of nearly 2,000 patches, we have found 30 detections.\n",
    "Let's use the information we have about these patches to show where they lie on our test image, drawing them as rectangles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.imshow(test_image, cmap='gray')\n",
    "ax.axis('off')\n",
    "\n",
    "Ni, Nj = positive_patches[0].shape\n",
    "indices = np.array(indices)\n",
    "\n",
    "for i, j in indices[labels == 1]:\n",
    "    ax.add_patch(plt.Rectangle((j, i), Nj, Ni, edgecolor='red',\n",
    "                               alpha=0.3, lw=2, facecolor='none'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the detected patches overlap and found the face in the image!\n",
    "Not bad for a few lines of Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caveats and Improvements\n",
    "\n",
    "If you dig a bit deeper into the preceding code and examples, you'll see that we still have a bit of work before we can claim a production-ready face detector.\n",
    "There are several issues with what we've done, and several improvements that could be made. In particular:\n",
    "\n",
    "### Our training set, especially for negative features, is not very complete\n",
    "\n",
    "The central issue is that there are many face-like textures that are not in the training set, and so our current model is very prone to false positives.\n",
    "You can see this if you try out the above algorithm on the *full* astronaut image: the current model leads to many false detections in other regions of the image.\n",
    "\n",
    "We might imagine addressing this by adding a wider variety of images to the negative training set, and this would probably yield some improvement.\n",
    "Another way to address this is to use a more directed approach, such as *hard negative mining*.\n",
    "In hard negative mining, we take a new set of images that our classifier has not seen, find all the patches representing false positives, and explicitly add them as negative instances in the training set before retraining the classifier.\n",
    "\n",
    "###  Our current pipeline searches at only one scale\n",
    "\n",
    "As currently written, our algorithm will miss faces that are not approximately 62×47 pixels.\n",
    "This can be straightforwardly addressed by using sliding windows of a variety of sizes, and resizing each patch with ``skimage.transform.resize`` before feeding it into the model.\n",
    "In fact, the ``sliding_window()`` utility used here is already built with this in mind.\n",
    "\n",
    "###  We should combine overlapped detection patches\n",
    "\n",
    "For a production-ready pipeline, we would prefer not to have 30 detections of the same face, but to somehow reduce overlapping groups of detections down to a single detection.\n",
    "This could be done via an unsupervised clustering approach (mean shift clustering is one good candidate for this), or via a procedural approach like *non-maximum suppression*, an algorithm common in machine vision.\n",
    "\n",
    "###  The pipeline should be streamlined\n",
    "\n",
    "Once we address these issues, it would also be nice to create a more streamlined pipeline for ingesting training images and predicting sliding-window outputs.\n",
    "This is where Python as a data science tool really shines: with a bit of work, we can take our prototype code and package it with a well-designed object-oriented API that gives the user the ability to easily use it.\n",
    "I will leave this as a proverbial \"exercise for the reader.\"\n",
    "\n",
    "###  More recent advances: deep learning\n",
    "\n",
    "Finally, I should add that HOG and other procedural feature extraction methods for images are no longer state-of-the-art techniques.\n",
    "Instead, many modern object detection pipelines use variants of deep neural networks: one way to think of neural networks is that they are an estimator that determines optimal feature extraction strategies from the data, rather than relying on the intuition of the user.\n",
    "An intro to these deep neural net methods is conceptually (and computationally!) beyond the scope of this notebook, although open tools like Google's <a href=\"https://www.tensorflow.org/\" target=\"_blank\">TensorFlow</a> have recently made deep learning approaches much more accessible than they once were."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "In this notebook we've taken a look at a few approaches to *feature engineering*, the process of taking raw data and turning it into the ``n_samples`` x ``n_features`` matrix that we need in order to model the data with scikit-learn. Feature engineering can range from simple (mapping categories to binary columns) to complex (specific operations on images to find gradients and other useful features). This lesson showed a few examples to illustrate how this might be done in practice."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}