Final Project Notebook Submission/notebook/VanlunenDaniel_NautiyalAniruddh_MunjuluriAnusha_ChallaUsha_FinalModel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Facial Keypoint Detection Project\n",
    "\n",
    "The objective of this project was to predict facial keypoints using Kaggle's Facial Keypoints Detection dataset. (https://www.kaggle.com/c/facial-keypoints-detection). Input consists of a list of pixels in grayscale and keypoint coordinates of 7000 images. Images are 96x96 pixels. Each predicted keypoint is specified by an (x,y) real-valued pair in the space of pixel indices. There are 15 keypoints, which represent elements of the face such as left_eye_center, right_eye_center, etc. and final goal was to predicts these keypoint positions for the test dataset.\n",
    "\n",
    "We explored various neural network structures (dense, convolutions, vgg16, inceptionV3 etc.) for our predictions and ultimately trained two models for project submission. Our final solution's predictions submitted to Kaggle produced a RMSE error of 1.71947 on the unseen test data labels. This placed us 4th on the public Kaggle leaderboard.\n",
    "\n",
    "## Group Members:\n",
    "- Vanlunen, Daniel\n",
    "- Nautiyal, Aniruddh\n",
    "- Munjuluri, Anusha\n",
    "- Challa, Usha\n",
    "\n",
    "### Notes about the directory structure:\n",
    "\n",
    "The zip file submitted contains 3 folders:\n",
    "\n",
    "<b>1. data:</b> This folder contains the competition data. *Idlookuptable*, *test*, *training*, and *samplesubmission* files in this folder are needed for this notebook to run. Training data should be unzipped before running the notebook.\n",
    "\n",
    " <b>2. notebook:</b> This folder contains this main subission notebook. It also contains a subfolder with a screenshot (img) of Kaggle score, and another subfolder (augmentdata) that defines a custom image generating class.\n",
    "\n",
    " <b>3. saved-models:</b> This folder contains the saved models weights of final models in .h5 files. It also has a subfolder that contains the training history in pickle files.\n",
    "\n",
    "### Remarks and Learnings: \n",
    "\n",
    "<b>1. Splitting data into two groups: </b>Initial data exploration showed that training images come from two different marking schemes. Hence, training two separate models for images that come with different set of keypoints (Group 1: 8 keypoints or less and Group 2: 9 keypoints or more, upto 30) reduced our training/validation loss.\n",
    "\n",
    "<b>2. Augmenting data: </b> Augmentating images (rotating, flipping etc.) was a nice way to artificially increase the variance in our training dataset and improve the model.\n",
    "\n",
    "<b>3. Convolutions vs Dense: </b>Convolutions significantly improved the results.\n",
    "\n",
    "<b>4. Choice of Optimizers: </b> Different optimizers worked differently under different circumstances -  Nadam converged faster but eventually Adam and Nadam reached the same loss. However, Adam did achieve a better minimum loss for simple models without convolutions. These optimizers gave a much better loss validation than using SGD optimizer. https://towardsdatascience.com/types-of-optimization-algorithms-used-in-neural-networks-and-ways-to-optimize-gradient-95ae5d39529f.\n",
    "\n",
    "<b>5. Normalization: </b>Normalization was helpful (both of the features and the output) and helped prevent divergence.\n",
    "\n",
    "<b>6. Batch Normalization: </b>Batch normalization worked well after the activation. https://github.com/ducha-aiki/caffenet-benchmark/blob/master/batchnorm.md\n",
    "\n",
    "<b>7. Dropout and Batch Normalization: </b> Dropout and batch normalization together seemed to cause weird differences between training and validation error due to \"variance shift\". One fix was to only apply dropout after all the batch normalization layers had been applied. http://arxiv.org/abs/1801.05134v1\n",
    "\n",
    "<b>8. Adding additional layers: </b> After bringing validation loss more in line with training loss through normalization and dropout, adding further complexity through more layers seemed to further reduce the validation loss.\n",
    "\n",
    "<b>9. Activation Function: </b> Though ReLu is the most popular activation function, Exponential Linear Units (ELU) in some cases give better results because their mean value is closer to zero and thus they allow less bias to propagate through the network. https://arxiv.org/pdf/1511.07289.pdf\n",
    "\n",
    "<b>10. Using pre-tuned models: </b> Though training pre-tuned models (\"from the zoo\") is likely a good choice in real world applications, our results training InceptionNetV3 and VGG16 did not show loss lower than our simpler architectures. This could be due to the small size of our training dataset and need more sample data to train these complex architectures.\n",
    "\n",
    "<b>11. Strides and Convolution Window Sizes: </b> Large strides and convolution window sizes generally worked worse than smaller convolutions with a stride of 1. Smaller windows were in line with the ideas presented in the InceptionV3 model paper of replacing larger windows with multiple small windows. https://arxiv.org/pdf/1512.00567.pdf\n",
    "\n",
    "<b>12. Training Time: </b> Depth and additional data significantly increased training time.\n",
    "\n",
    "<b>13. Use of GPUs: </b> GPUs gave a ~30x speedup over CPU training. They were definitely necessary to quickly train, test and evaluate various models.\n",
    "\n",
    "### Further exploration:\n",
    "Given more time, there are a number of other things we would have evaluated:\n",
    "- Use cross validation to see if our best architectures hold up when the seed value is not 42.\n",
    "- Testing more advanced, inception based architectures.\n",
    "- Utilizing the images with missing keypoints (total keypoints != 30 or 8). Our solution used only images with 30 or 8 keypoints train models.\n",
    "- Training a separate model for the eye keypoints that uses the data in both groups of images, since it appears the two labelling schemes were the same for the center of the eye keypoints.\n",
    "- Tuning the amount of data augmentation used.\n",
    "- Importing layers and hyperparameter settings from the other group's model, to see if those settings offer improvement over the current group's model.\n",
    "- Use Ensemble architecture to split dataset based on number of available keypoints, train on separate models and combine results.\n",
    "\n",
    "\n",
    "### Location of other models tried:\n",
    "https://drive.google.com/file/d/1dK069Q08DIZ6g_HIi4osF20DjV2zktC3/view?usp=sharing\n",
    "https://drive.google.com/drive/folders/1IcfVLCy_btNzYmqzvKzyVPFQQkKrYVvZ?usp=sharing\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pandas.io.parsers import read_csv\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "import os\n",
    "import pickle\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from keras.models import Sequential, load_model, Model\n",
    "from keras.layers import Dense, Conv2D, MaxPooling2D, Flatten, Dropout, BatchNormalization\n",
    "from keras.optimizers import Nadam\n",
    "from keras.callbacks import EarlyStopping, Callback, History\n",
    "\n",
    "from augmentdata.CustImageDataGenerator import CustImageDataGenerator,CustNumpyArrayIterator\n",
    "\n",
    "from tensorflow.python.client import device_lib\n",
    "\n",
    "print(device_lib.list_local_devices()) # confirm using GPU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TRAIN_DATA = '../data/training.csv'                # train dataset downloaded from Kaggle\n",
    "TEST_DATA = '../data/test.csv'                     # test dataset downloaded from Kaggle\n",
    "IMAGE_ROWS = 96                                    # number of pixel rows\n",
    "IMAGE_COLS = 96                                    # number of pixel columns\n",
    "INPUT_SHAPE = (IMAGE_ROWS, IMAGE_COLS, 1)          # shape of input to first Conv2D layer of the model\n",
    "RETRAIN = False                                    # bool to load and use existing saved models\n",
    "VERBOSE_TRAIN = True                               # bool to show/hide progress while training a model\n",
    "NUM_KEYPOINTS = 30                                 # maximum no. of facial keypoints for any image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EDA\n",
    "\n",
    "We began by examining images from train dataset, which is organised as (x,y) coodinates of each of the 15 facial features for a total of 30 keypoints. The 31st column contains 96x96 image pixel array data coded as raw grayscale values from 0 to 255. There are a total of about 7000 images in the train dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Setting up the data for EDA\n",
    "\n",
    "df_train = pd.read_csv(TRAIN_DATA)\n",
    "Y = np.array(df_train[df_train.columns.difference(['Image'])])\n",
    "X = df_train['Image']\n",
    "labels = list(df_train.columns.difference(['Image']))\n",
    "\n",
    "img_vec_len = IMAGE_ROWS*IMAGE_COLS                             # images pixel grid size\n",
    "\n",
    "imgArray = np.zeros((X.shape[0], img_vec_len), dtype=int)       # temporary array to save each image as numpy array\n",
    "\n",
    "idx=0\n",
    "for i in X.keys(): \n",
    "    imgArray[idx] = np.fromstring(X[i], dtype=int, sep=' ')\n",
    "    idx += 1\n",
    "X = np.reshape( imgArray, (X.shape[0], IMAGE_ROWS, IMAGE_COLS, 1) )\n",
    "print(\"Total images in train dataset: \", X.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Function to subplot a group of images, and label the ones with missing keypoints distinctly\n",
    "\n",
    "def plot_images(images, points, type='actual', subplotting=False, gridRows=0, gridCols=0, \n",
    "                imageIndices=1, subtitles=True, title=None, labelsList=[] ):\n",
    "    \n",
    "    plt.figure(figsize=(4*gridCols, 4*gridCols))    \n",
    "    img_nums = images.shape[0]\n",
    "    points_nums = points.shape[0]\n",
    "    \n",
    "    if ( ( img_nums != points_nums) | ( img_nums != imageIndices.shape[0] ) ):\n",
    "        raise ValueError(\"Mismatch in number of images and keypoints' rows passed to plot_images().\")\n",
    "    \n",
    "    \n",
    "    for thisImg in range(0, gridRows*gridCols ):\n",
    "        \n",
    "        if subplotting:\n",
    "            plt.subplot(gridRows, gridCols, thisImg + 1)\n",
    "            noKeypNums  = np.isnan(points[thisImg]).sum()\n",
    "            \n",
    "            if subtitles:\n",
    "                if( noKeypNums == 0 ):                            # no missing keypoints (group1)\n",
    "                    plt.title(\"#: \" + str(imageIndices[thisImg]) + \n",
    "                              \",  Points: \" + str(NUM_KEYPOINTS - noKeypNums), color='k')\n",
    "                \n",
    "                elif( ( noKeypNums > 0) & (noKeypNums < 22 ) ):   # (1,21) missing keypoints\n",
    "                    plt.title(\"#: \" + str(imageIndices[thisImg]) + \n",
    "                              \",  Points: \" + str(NUM_KEYPOINTS - noKeypNums), color='m')\n",
    "                \n",
    "                elif( ( noKeypNums == 22 ) ):                     # 22 missing keypoints (group2)\n",
    "                    plt.title(\"#: \" + str(imageIndices[thisImg]) + \n",
    "                              \",  Points: \" + str(NUM_KEYPOINTS - noKeypNums), color='b')\n",
    "                \n",
    "                else:                                             # > 22 missing keypoints\n",
    "                    plt.title(\"#: \" + str(imageIndices[thisImg]) + \n",
    "                              \",  Points: \" + str(NUM_KEYPOINTS - noKeypNums), color='r')\n",
    "        \n",
    "        plt.imshow(np.reshape(images[thisImg,:],(96,96)), cmap = 'gray')\n",
    "\n",
    "        x = 0\n",
    "        for idx in range(0, points[thisImg].shape[0]):\n",
    "            label = labelsList[idx]\n",
    "            if label[-1]=='x':\n",
    "                x = points[thisImg, idx]\n",
    "            else:\n",
    "                if label in ['left_eye_center_y',\n",
    "                             'left_eye_inner_corner_y', \n",
    "                             'left_eye_outer_corner_y', \n",
    "                             'left_eyebrow_inner_end_y', \n",
    "                             'left_eyebrow_outer_end_y',\n",
    "                             'mouth_left_corner_y'\n",
    "                            ]:\n",
    "                    if(type=='actual'):\n",
    "                        plt.plot(x, points[thisImg, idx], 'c<')\n",
    "                    else:\n",
    "                        plt.plot(x, points[thisImg, idx], 'c*')\n",
    "                        \n",
    "                elif label in ['right_eye_center_y',\n",
    "                             'right_eye_inner_corner_y', \n",
    "                             'right_eye_outer_corner_y', \n",
    "                             'right_eyebrow_inner_end_y', \n",
    "                             'right_eyebrow_outer_end_y',\n",
    "                              'mouth_right_corner_y']:\n",
    "                    if(type=='actual'):\n",
    "                        plt.plot(x, points[thisImg, idx], 'r>')\n",
    "                    else:\n",
    "                        plt.plot(x, points[thisImg, idx], 'r*')\n",
    "                \n",
    "                else:\n",
    "                    if(type=='actual'):\n",
    "                        plt.plot(x, points[thisImg, idx], 'mo')\n",
    "                    else:\n",
    "                        plt.plot(x, points[thisImg, idx], 'm*')\n",
    "                    \n",
    "        plt.axis('off')\n",
    "    \n",
    "    if (title != None):\n",
    "        plt.suptitle(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Function to plot an array of image indices\n",
    "\n",
    "idx_max = df_train.shape[0]                      # all images\n",
    "grid_cols = 4                                    # grid columns size for a subplot of images\n",
    "grid_rows = 4                                    # grid rows size for images subplot\n",
    "subImgNum = grid_cols * grid_rows\n",
    "\n",
    "def plot_img_group( thisGroup, dataset='train', denorm=False, thisLabels=labels, thisSubTitle=True, thisTitle=None ):\n",
    "\n",
    "    if(dataset == 'train'):\n",
    "        thisX = X\n",
    "        thisY = Y\n",
    "    elif(dataset == 'group1'):\n",
    "        thisX = X1\n",
    "        thisY = Y1\n",
    "    elif(dataset == 'group2'):\n",
    "        thisX = X2\n",
    "        thisY = Y2     \n",
    "    \n",
    "    thisSubsetX = np.zeros( (subImgNum, IMAGE_ROWS, IMAGE_COLS, 1), dtype=float)\n",
    "    thisSubsetY = np.zeros( (subImgNum, thisY.shape[1]), dtype=float)\n",
    "    img_indices = np.zeros( (subImgNum, 1), dtype=int)\n",
    "\n",
    "    img_sub = 0                                  # local iterator for images in subplot\n",
    "    flushed = False\n",
    "    for img in thisGroup:\n",
    "        \n",
    "        if( ( (img_sub + 1 )  % subImgNum ) != 0 ):\n",
    "            thisSubsetX[img_sub,:] = thisX[img-1,:]\n",
    "            if denorm:\n",
    "                thisSubsetY[img_sub,:] = 48*thisY[img-1,:] + 48\n",
    "            else:\n",
    "                thisSubsetY[img_sub,:] = thisY[img-1,:]\n",
    "            img_indices[img_sub] = img\n",
    "            img_sub += 1\n",
    "            flushed = False\n",
    "            \n",
    "        else:\n",
    "            thisSubsetX[img_sub,:] = thisX[img-1,:]\n",
    "            if denorm:\n",
    "                thisSubsetY[img_sub,:] = 48*thisY[img-1,:] + 48\n",
    "            else:\n",
    "                thisSubsetY[img_sub,:] = thisY[img-1,:]\n",
    "            img_indices[img_sub] = img\n",
    "            \n",
    "            # plot when all images for the subplot are accumulated\n",
    "            plot_images(images=thisSubsetX, points=thisSubsetY, subplotting=True, \n",
    "                        gridRows=grid_rows, gridCols=grid_cols, imageIndices=img_indices, \n",
    "                        subtitles=thisSubTitle, title=thisTitle, labelsList=thisLabels )\n",
    "            \n",
    "            # reset subplot indexing pointer and subplot image/keypoints buckets\n",
    "            img_sub = 0\n",
    "            flushed = True\n",
    "            thisSubsetX = np.zeros( (subImgNum, IMAGE_ROWS, IMAGE_COLS, 1), dtype=int)\n",
    "            thisSubsetY = np.zeros( (subImgNum, Y.shape[1]), dtype=float)\n",
    "    \n",
    "    if not flushed:                          # for images leftover from partial subplot grid\n",
    "\n",
    "        thisGridRows = ( (img_sub - 1) // grid_rows ) + 1\n",
    "        if( thisGridRows > 1 ):\n",
    "            thisGridCols = grid_cols\n",
    "        else:\n",
    "            thisGridCols = img_sub\n",
    "\n",
    "        plot_images(images=thisSubsetX, points=thisSubsetY, subplotting=True, \n",
    "                    gridRows=thisGridRows, gridCols=thisGridCols, imageIndices=img_indices, \n",
    "                    subtitles=thisSubTitle, title=thisTitle, labelsList=thisLabels )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Images in train dataset were broadly split into two groups. The first group had more than 8 keypoints, up to maximum of 30, and starting at index 1(or row 1 in .csv) upto index 2284. The images after index 2284 had a maximum of 8 keypoints. \n",
    "\n",
    "We identified approx. 320 images for in-depth examination after an initial review of the dataset. These images had some peculiarities, examples of these are discussed below: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Peculiarity in Group1 Images:\n",
    "\n",
    "##### Eyes closed: \n",
    "For identifying the center of eye keypoint, these images had an approximate labelling of that keypoint, without providing a high level of feature extraction for the model e.g. the absence contrast between the white sclera and the relatively darker eyelids or iris. Since the dataset had more than 300+ such images, we assumed these image types to be common in the test dataset. Hence, were retained in the training dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "imgs_eyes_closed    = np.array([152, 242, 288, 1676])\n",
    "\n",
    "print(\"\\nSample images with eyes closed\")\n",
    "plot_img_group( imgs_eyes_closed )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Peculiarity in Group2 Images:\n",
    "\n",
    "##### Highly blurred:\n",
    "Group 2 had a higher number of images that were extremely blurred and failed to provide any significant feature extractions. These images were not useful in training the model, and were dropped from our training dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "imgs_high_blur    = np.array([2322, 2574, 2916, 6605])\n",
    "\n",
    "print(\"\\nSample images that are highly blurred\")\n",
    "plot_img_group( imgs_high_blur )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Peculiarities common in both groups:\n",
    "\n",
    "##### Partially/fully covered keypoints:\n",
    "Some images had part of the faces, near eyes, partially covered with hair, hat, or other form of facial obstruction (by objects like props, sunglasses, partially covered in dark shadow). These images were marked with approximate location of obscured keypoints. Including these imagees would not provide any benefits in training the model. Therefore, most of these were eliminated from training dataset. \n",
    "\n",
    "##### A second partial/full face:\n",
    "Some images had more than one face either partially or fully appearing in the foreground/background. Although some of these were labelled correctly, we assumed these could mislead the model, especially in cases where more than one face is prominent. We decided to drop images where the primary face was not sufficiently prominent.\n",
    "\n",
    "##### Wrong/bad labels:\n",
    "These images had wrong or bad labelling and were removed from the train dataset. \n",
    "\n",
    "##### Missing keypoints:\n",
    "Images that did not fall into Group 1 or 2 (images where identified keypoints were not 30 or 8) were either out of frame, had the face viewed sideways, or were overshadowed and hence did not contain all keypoints. While these images were valid for training, we found that the model fitting function expected a coherent specification of keypoints for all input samples. Developing a solution to accurately handle all these scenarios requires additional models. Due to time constraints, we decided to limit our scope to include only images with 30 or 8 keypoints.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "imgs_obstructed     = np.array([1697, 1947, 5628, 6755])\n",
    "imgs_2nd_face       = np.array([1820, 2064, 4264, 4491])\n",
    "imgs_bad_labels     = np.array([1748, 1878, 2200, 6493])\n",
    "imgs_miss_keyps     = np.array([1709, 1875, 1882, 2322])\n",
    "\n",
    "print(\"\\nSample images with partial or fully covered/obscured keypoints\")\n",
    "plot_img_group( imgs_obstructed )\n",
    "\n",
    "print(\"\\nSample images with a second partial/full face\")\n",
    "plot_img_group( imgs_2nd_face )\n",
    "\n",
    "print(\"\\nSample images with bad/wrong labelling\")\n",
    "plot_img_group( imgs_bad_labels )\n",
    "\n",
    "print(\"\\nSample images with missing keypoints\")\n",
    "plot_img_group( imgs_miss_keyps )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It was evident from EDA that a different labeling scheme was used for images in the two groups. For example, images in group 2 had keypoints marked underneath the nose while images in group 1 had tip of the nose as a keypoint. Group2 also marks center of the bottom lip, in the center of the mouth, often on teeth.\n",
    "\n",
    "These differences in keypoint location identification reconfirmed our approach to develop 2 separate models, one for group2 images with 8 keypoints, and another for group1 with 30.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Final indices of images to be dropped from the dataset\n",
    "\n",
    "IDX_BAD_IMAGES = np.array( [1621, 1862, 1748, 1878, 1927, 2200, 2431, 2584, 2647, \n",
    "                            2671, 2765, 4198, 1627, 1628, 1637, 1957, 4477, 1820, \n",
    "                            2064, 2089, 2091, 2109, 2195, 4264, 4491, 6490, 6493, \n",
    "                            6494, 1655, 2096, 2454, 3206, 3287, 5628, 5653, 6754, \n",
    "                            6755, 2321, 2322, 2414, 2428, 2462, 2574, 2584, 2663, \n",
    "                            2691, 2694, 2830, 2910, 2916, 3126, 3176, 3291, 3299, \n",
    "                            3361, 4061, 4483, 4484, 4494, 4766, 4809, 4837, 4880, \n",
    "                            4905, 5068, 5362, 5566, 5868, 6535, 6538, 6588, 6605, \n",
    "                            6659, 6724, 6733, 6753, 6758, 6766, 6907 ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading the data\n",
    "This section loads data for two the groups (Separate models were trained on each of these groups because they appeared to have different labelling schemes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Function to clean up the train dataset, normalize it, drop bad images & labels, \n",
    "# and finally split the dataset into 2, for group1 and group2 modelling.\n",
    "\n",
    "def loaderV2(test=False, seed=None, keeplabels=None):\n",
    "    \n",
    "    if seed:\n",
    "        np.random.seed(seed)\n",
    "    fileloc = TEST_DATA if test else TRAIN_DATA\n",
    "    \n",
    "    df = read_csv(fileloc)\n",
    "    \n",
    "    df['Image'] = df['Image'].apply(lambda x: np.fromstring(x, sep=' '))\n",
    "    \n",
    "    if keeplabels:\n",
    "        df = df[list(keeplabels) + ['Image']]\n",
    "        \n",
    "    X = np.vstack(df['Image'])\n",
    "    \n",
    "    if not test:                                       # process train dataset\n",
    "        Y = df[df.columns.difference(['Image'])].values\n",
    "        Y = Y.astype(np.float32)\n",
    "        \n",
    "        # remove rows having bad images or labels\n",
    "        X = np.delete( X, (IDX_BAD_IMAGES - 1), axis=0 )\n",
    "        Y = np.delete( Y, (IDX_BAD_IMAGES - 1), axis=0 )\n",
    "        \n",
    "        # normalize - by pixel across the whole dataset subtract mean and divide by stdev\n",
    "        X = X - np.tile(np.mean(X,axis=0),(X.shape[0],1))\n",
    "        X = X / np.tile(np.std(X,axis=0),(X.shape[0],1))\n",
    "        X = X.astype(np.float32)\n",
    "    \n",
    "        Y = (Y - 48) / 48                     # this helps, but tanh on output doesnt\n",
    "        shuffle = np.random.permutation(np.arange(X.shape[0]))\n",
    "        X, Y = X[shuffle], Y[shuffle]\n",
    "    \n",
    "        X = X.reshape(-1, 96, 96, 1)\n",
    "        \n",
    "        # split X and Y into dataset for model1 (more than 8 keypoints) and model2 (less than 8 keypoints)\n",
    "        X_model1 = np.zeros( (X.shape[0], 96,96,1), dtype=np.float32)\n",
    "        X_model2 = np.zeros( (X.shape[0], 96,96,1), dtype=np.float32)\n",
    "        Y_model1 = np.zeros( Y.shape, dtype=float)\n",
    "        Y_model2 = np.zeros( Y.shape, dtype=float)\n",
    "        tempIdx1 = 0\n",
    "        tempIdx2 = 0\n",
    "        \n",
    "        for thisIdx in range(0, Y.shape[0]):\n",
    "            numKeyps  = NUM_KEYPOINTS - np.isnan(Y[thisIdx]).sum()\n",
    "            \n",
    "            if( ( numKeyps > 8 ) ):\n",
    "                X_model1[tempIdx1] = X[thisIdx,:,:]\n",
    "                Y_model1[tempIdx1] = Y[thisIdx,:]\n",
    "                tempIdx1 = tempIdx1 + 1\n",
    "            else:\n",
    "                X_model2[tempIdx2] = X[thisIdx,:,:]\n",
    "                Y_model2[tempIdx2] = Y[thisIdx,:]\n",
    "                tempIdx2 = tempIdx2 + 1\n",
    "    \n",
    "        # remove empty rows\n",
    "        drop_idx1 = []\n",
    "        drop_idx2 = []\n",
    "        \n",
    "        for idx in range(0, X.shape[0]):\n",
    "            if( (np.all(Y_model1[idx] == 0)) | (np.isnan(Y_model1[idx]).sum() != 0) ):\n",
    "                drop_idx1.append(idx)\n",
    "            if( (np.all(Y_model2[idx] == 0)) | (np.isnan(Y_model2[idx]).sum() != 22) ):\n",
    "                drop_idx2.append(idx)\n",
    "        \n",
    "        X_model1 = np.delete( X_model1, np.array(drop_idx1), axis=0 )\n",
    "        Y_model1 = np.delete( Y_model1, np.array(drop_idx1), axis=0 )\n",
    "        X_model2 = np.delete( X_model2, np.array(drop_idx2), axis=0 )\n",
    "        Y_model2 = np.delete( Y_model2, np.array(drop_idx2), axis=0 )            \n",
    "        \n",
    "        # remove empty columns, setup lists of labels\n",
    "        labels = df.columns.difference(['Image'])\n",
    "        labels1 = labels\n",
    "        labels2 = []\n",
    "        drop_idx3 = []\n",
    "        for idx in range(0, Y.shape[1]):\n",
    "            if( (np.all(Y_model2[:,idx] == 0)) | (np.isnan(Y_model2[:,idx]).sum() != 0) ):\n",
    "                drop_idx3.append(idx)\n",
    "            else:\n",
    "                labels2.append(labels[idx])\n",
    "        Y_model2 = np.delete( Y_model2, np.array(drop_idx3), axis=1 ) \n",
    "        \n",
    "        # return the original dataset and the group splits\n",
    "        return X_model1, Y_model1, labels1, X_model2, Y_model2, labels2, X, Y, labels\n",
    "    \n",
    "    else:                                          # process test dataset\n",
    "        Y = None\n",
    "        \n",
    "        # normalize - by pixel across the whole dataset subtract mean and divide by stdev\n",
    "        X = X - np.tile(np.mean(X,axis=0),(X.shape[0],1))\n",
    "        X = X / np.tile(np.std(X,axis=0),(X.shape[0],1))\n",
    "        X = X.reshape(-1, 96, 96, 1)\n",
    "        labels = df.columns.difference(['Image'])\n",
    "        \n",
    "        return X, Y, labels\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X1, Y1, labels1,   X2, Y2, labels2,   X, Y, labels = loaderV2(seed=42)\n",
    "\n",
    "print(\"Group1 data Y1 shape: \", Y1.shape, \", X1 shape: \", X1.shape)\n",
    "print(\"Group2 data Y2 shape: \", Y2.shape, \", X2 shape: \", X2.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Fitting\n",
    "\n",
    "### Model fitting function\n",
    "Below is general model fitting function to train model architectures used in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fit_model(model, data, modelname,\n",
    "              generator=None,retrain=RETRAIN,\n",
    "              epochs=10000, patience=1000,optimizer=Nadam()):\n",
    "    # check if the user wants to retrain or if the saved model doesn't exist\n",
    "    if retrain or not os.path.exists('../saved-models/' + modelname + '.h5'):\n",
    "        # data setup\n",
    "        X_train = data[0]\n",
    "        y_train = data[2]\n",
    "        if len(data) == 4:\n",
    "            valid_dat = (data[1], data[3])\n",
    "        else:\n",
    "            valid_dat = None\n",
    "\n",
    "        # default optimization routine to use Nadam and minimize mse\n",
    "        model.compile(loss='mse', optimizer=optimizer)\n",
    "        \n",
    "        # set an early stopping criteria\n",
    "        if valid_dat:\n",
    "            earlystop = EarlyStopping(monitor='val_loss',\n",
    "                                     patience=patience,\n",
    "                                     verbose=1,\n",
    "                                     mode=\"auto\")\n",
    "            \n",
    "        else:\n",
    "            earlystop = EarlyStopping(monitor='loss',\n",
    "                                     patience=patience,\n",
    "                                     verbose=1,\n",
    "                                     mode=\"auto\")\n",
    "        \n",
    "        callbacks = [earlystop]\n",
    "        \n",
    "        # fit the model (with a data generator if present)\n",
    "        if generator:\n",
    "            history = model.fit_generator(generator,\n",
    "                        epochs=epochs,\n",
    "                        steps_per_epoch=data[0].shape[0]//32,\n",
    "                        callbacks=callbacks,\n",
    "                        validation_data=valid_dat\n",
    "             )\n",
    "        else:\n",
    "            history = model.fit(X_train, y_train,\n",
    "                                epochs=epochs,\n",
    "                                batch_size=32,\n",
    "                                callbacks=callbacks,\n",
    "                                validation_data=valid_dat,\n",
    "                                verbose=VERBOSE_TRAIN\n",
    "                     )\n",
    "        # save the model weights\n",
    "        model.save('../saved-models/'+ modelname + '.h5')\n",
    "        \n",
    "        # save the model loss history\n",
    "        with open('../saved-models/histories/'+modelname+'_hist',\n",
    "                  'wb') as file_pi:\n",
    "            pickle.dump(history.history, file_pi)\n",
    "        history = history.history\n",
    "    \n",
    "    # if the user doesnt want to retrain, load the weights and model history\n",
    "    else:\n",
    "        model = load_model('../saved-models/'+modelname+'.h5')\n",
    "        history = pickle.load(open( \"../saved-models/histories/\" + modelname + '_hist',\n",
    "                                   \"rb\" ))\n",
    "        \n",
    "    return history, model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data augmenting function\n",
    "Keras offers an image data generator class that automatically transforms input data images. Data augmentation helps artifically increase the size of training data set by adding variance to training images without distorting their meaning. This keras data generator does not change the labels of the images. \n",
    "\n",
    "In our case, when the image was changed, we needed the keypoint locations to adjust accordingly. Therefore, we had to create a custom subclass of the Keras ImageDataGenerator that not only transforms the images, but also transforms the labels. \n",
    "\n",
    "Our custom generator implements below transformations to image data:\n",
    "1. randomly rotates images by up to 5 degrees\n",
    "2. flips images horizontally\n",
    "3. translates images by up to 5% of the total height/weight of the image.  \n",
    "\n",
    "These transformations were only applied if they would not move the keypoints outside of the image. Refer to CustImageDataGenerator.py file for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generator for the models that predict 30 keypoints\n",
    "datagen = CustImageDataGenerator(\n",
    "    rotation_range=5. #degrees\n",
    "     ,horizontal_flip=True\n",
    "     ,width_shift_range=.05 # percent of image width\n",
    "     ,height_shift_range=.05 # percent of image height\n",
    "    ).flow(X1,Y1,whichlabels=list(labels1), batch_size=32)\n",
    "\n",
    "# Generator for the models that predict 8 keypoints\n",
    "g2_datagen = CustImageDataGenerator(\n",
    "    rotation_range=5. #degrees\n",
    "     ,horizontal_flip=True\n",
    "     ,width_shift_range=.05 # percent of image width\n",
    "     ,height_shift_range=.05 # percent of image height\n",
    "    ).flow(X2,Y2,whichlabels=list(labels2), batch_size=32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build and train Group 1 Images' Model (more than 8 and upto 30 keypoints to predict)\n",
    "**Note**: In order to retrain the model, retrain=True should be passed. **Caution:** this will overwrite the old saved model in the saved-models directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "g1_model3_new = Sequential()\n",
    "\n",
    "g1_model3_new.add(Conv2D(32,\n",
    "                 (6, 6),\n",
    "                 activation='relu',\n",
    "                input_shape=INPUT_SHAPE))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "g1_model3_new.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "g1_model3_new.add(Conv2D(filters=64,\n",
    "                 kernel_size=(5, 5),\n",
    "                 activation='relu'))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "g1_model3_new.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "g1_model3_new.add(Conv2D(filters=256,\n",
    "                 kernel_size=(4, 4),\n",
    "                 activation='relu'))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "g1_model3_new.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "g1_model3_new.add(Conv2D(filters=64,\n",
    "                 kernel_size=(3,3),\n",
    "                 activation='relu'))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "g1_model3_new.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "g1_model3_new.add(Conv2D(filters=128,\n",
    "                 kernel_size=(2, 2),\n",
    "                 activation='relu'))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "g1_model3_new.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "\n",
    "g1_model3_new.add(Flatten())\n",
    "g1_model3_new.add(Dense(500, activation = \"relu\"))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "\n",
    "g1_model3_new.add(Dense(500, activation = \"relu\"))\n",
    "g1_model3_new.add(BatchNormalization())\n",
    "g1_model3_new.add(Dropout(.3))\n",
    "\n",
    "g1_model3_new.add(Dense(30))\n",
    "\n",
    "print(g1_model3_new.summary())\n",
    "g1_model3_new_hist, g1_model3_new = fit_model(g1_model3_new, [X1, None, Y1],\n",
    "                                'g1_CNN_aug_addedLayers',datagen,retrain=RETRAIN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build and train Group 2 Images' Model (maximum 8 keypoints to predict)\n",
    "**Note**: In order to retrain the model, retrain=True should be passed. **Caution:** this will overwrite the old saved model in the saved-models directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "g2_model4 = Sequential()\n",
    "\n",
    "g2_model4.add(Conv2D(filters=64,\n",
    "                 kernel_size=(6, 6),\n",
    "                 strides=1,\n",
    "                 activation='elu',\n",
    "                 input_shape=INPUT_SHAPE))\n",
    "g2_model4.add(BatchNormalization())\n",
    "g2_model4.add(Dropout(.1))\n",
    "\n",
    "g2_model4.add(Conv2D(filters=128,\n",
    "                 kernel_size=(5, 5),\n",
    "                 strides=1,\n",
    "                 activation='elu'))\n",
    "g2_model4.add(BatchNormalization())\n",
    "g2_model4.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "g2_model4.add(Dropout(.2))\n",
    "\n",
    "g2_model4.add(Conv2D(filters=256,\n",
    "                 kernel_size=(4, 4),\n",
    "                 activation='elu'))\n",
    "g2_model4.add(BatchNormalization())\n",
    "g2_model4.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "g2_model4.add(Dropout(.2))\n",
    "\n",
    "g2_model4.add(Conv2D(filters=512,\n",
    "                 kernel_size=(3, 3),\n",
    "                 activation='elu'))\n",
    "g2_model4.add(BatchNormalization())\n",
    "g2_model4.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "g2_model4.add(Dropout(.3))\n",
    "\n",
    "g2_model4.add(Conv2D(filters=512,\n",
    "                 kernel_size=(2, 2),\n",
    "                 activation='elu'))\n",
    "g2_model4.add(BatchNormalization())\n",
    "g2_model4.add(MaxPooling2D(pool_size=(2, 2)))\n",
    "g2_model4.add(Dropout(.4))\n",
    "\n",
    "g2_model4.add(Flatten())\n",
    "g2_model4.add(Dense(500, activation = \"elu\"))\n",
    "g2_model4.add(BatchNormalization())\n",
    "g2_model4.add(Dropout(.4))\n",
    "\n",
    "g2_model4.add(Dense(100, activation = \"elu\"))\n",
    "g2_model4.add(BatchNormalization())\n",
    "\n",
    "g2_model4.add(Dense(8))\n",
    "\n",
    "print(g2_model4.summary())\n",
    "\n",
    "g2_model4_hist, g2_model4 = fit_model(g2_model4,[X2, None, Y2],\n",
    "                                'g2_CNNv2_aug', generator=g2_datagen\n",
    "                                     ,retrain=RETRAIN, patience=100)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictions\n",
    "### Load test images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "out_images, _ , _ = loaderV2(test=True, seed=None, keeplabels=None)\n",
    "print(out_images.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create predictions\n",
    "Predictions with both models: g1 predicts all 30 keypoints, g2 predicts 8 keypoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "g1_prediction = g1_model3_new.predict(out_images)\n",
    "g2_prediction = g2_model4.predict(out_images)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Prepare predictions into submission format\n",
    "#### Adding keypoints per test image to be predicted\n",
    "This is necessary to know which model's predictions should be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# aggregate the number of keypoints needed for each distinct image in test dataset\n",
    "\n",
    "IdLookupTable = read_csv('../data/IdLookupTable.csv')\n",
    "\n",
    "rowIDs = np.array(IdLookupTable['RowId'])\n",
    "# keypoints to be predicted for each imageID\n",
    "imgs_numKeyps = np.zeros(max(IdLookupTable['ImageId']), dtype=int)\n",
    "\n",
    "thisImgKeyps = 0\n",
    "for rowIdx in range(0, rowIDs.shape[0]):\n",
    "    thisImgID = IdLookupTable.loc[rowIdx,'ImageId']\n",
    "    imgs_numKeyps[thisImgID - 1] += 1\n",
    "\n",
    "print(\"Total test images: {}\".format(imgs_numKeyps.shape) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Use appropriate model for predicting each image's keypoints\n",
    "If the test data requests 8 or fewer keypoints use the group 2 model, else use group 1 model. We tried using mean of the predictions from both models, for common keypoints, as left_eye_center and right_eye_center; but these yielded higher RMSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "label_locs1 = {}\n",
    "label_locs2 = {}\n",
    "for i, label in enumerate(labels1):\n",
    "    label_locs1[label]=i\n",
    "\n",
    "for i, label in enumerate(labels2):\n",
    "    label_locs2[label]=i\n",
    "    \n",
    "IdLookupTable['test'] = IdLookupTable['FeatureName'].replace(label_locs1)\n",
    "\n",
    "thisRowId = 0\n",
    "modelUsed = np.zeros(IdLookupTable['RowId'].shape[0], dtype=int)\n",
    "\n",
    "for imgIdx in range(0, imgs_numKeyps.shape[0]):\n",
    "    \n",
    "    if( imgs_numKeyps[imgIdx] > 8):\n",
    "        for keypIdx in range(0, imgs_numKeyps[imgIdx]):\n",
    "            map_label_idx1 = label_locs1[IdLookupTable.loc[thisRowId, 'FeatureName']]\n",
    "            IdLookupTable.loc[thisRowId, 'Location'] =  (48*g1_prediction[imgIdx, map_label_idx1]) + 48\n",
    "            modelUsed[thisRowId] = 1\n",
    "            thisRowId += 1\n",
    "    else:\n",
    "        for keypIdx in range(0, imgs_numKeyps[imgIdx]):\n",
    "            map_label_idx2 = label_locs2[IdLookupTable.loc[thisRowId, 'FeatureName']]\n",
    "            IdLookupTable.loc[thisRowId, 'Location'] =  (48*g2_prediction[imgIdx, map_label_idx2]) + 48\n",
    "            \n",
    "            modelUsed[thisRowId] = 2\n",
    "            thisRowId += 1\n",
    "\n",
    "IdLookupTable['Location'] = (IdLookupTable['Location'].\n",
    "                             where(IdLookupTable['Location']<=96, 96).\n",
    "                             where(IdLookupTable['Location']>=0, 0)\n",
    "                            )\n",
    "IdLookupTable.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Creating the submission"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Submission = IdLookupTable[['RowId','Location']]\n",
    "Submission.to_csv(path_or_buf='Final_FacialKeypoints.csv',\n",
    "                  index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation\n",
    "Our final model's predictions submitted to Kaggle produce a RMSE error of 1.71947 with the unseen test data labels. This places us 4th on the public Kaggle leaderboard. This is much better than the baseline RMSE we set at 3.138, by just using the mean value of the labels as our prediction. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](img/Score.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}