hands-on/040_mnist_mlp_courageous.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MNIST: multilayer perceptrons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multilayer perceptrons (MLP) have been around for more than 25 years. They are suited to both regression and classification tasks. Here we will use a classic MLP to recognize handwritten digits (MNIST data set), as well as an MLP with a number of more recent twists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import keras\n",
    "from keras.datasets import mnist\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import FunctionTransformer, OneHotEncoder\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Obtaining the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Keras' datasets module we have a handle to the MNIST dataset we want to use in this notebook.  Download the training and test set for this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(x_train, y_train), (x_test, y_test) = mnist.load_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can start doing machine learning on the data, some preparations are required.  In this case, the data set is clean, which simplifies this step considerably.  Although it would be better to create a pipeline using the scikit-learn framework, we'll do the preparation by hand in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transforming the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than use the 28 $\\times$ 28 images as input directly, we reshape each image to a 784 array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_reshaper = FunctionTransformer(lambda x: x.reshape(x.shape[0], -1),\n",
    "                                     validate=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most training algorithms work better if the input is normalized between 0 and 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_normalizer = FunctionTransformer(lambda x: x.astype(np.float32)/255.0,\n",
    "                                       validate=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both preprocessing steps can be combined into a single pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_pipeline = Pipeline([\n",
    "    ('reshaper', input_reshaper),\n",
    "    ('normalizer', input_normalizer),\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now fit and transform the training input, and transform the testing input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = input_pipeline.fit_transform(x_train)\n",
    "x_test = input_pipeline.transform(x_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output can be transformed to categorical data, i.e., one category for each digit, rather than a `uint8`. This is a one-hot encoding, so the output is now an array consisting of a single 1.0 value, and nine 0.0 values.  Before calling the `OneHotEncoder`, the output array has to be reshaped.  Note that the type of the output is now `float64`, which we can reduce to `float32` to reduce memory requirements and speed up computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_reshaper = FunctionTransformer(lambda x: x.reshape(-1, 1),\n",
    "                                      validate=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_encoder = OneHotEncoder(categories='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_type_changer = FunctionTransformer(lambda x: x.astype(np.float32),\n",
    "                                          validate=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_pipeline = Pipeline([\n",
    "    ('reshaper', output_reshaper),\n",
    "    ('binarizer', output_encoder),\n",
    "    ('type_changer', output_type_changer),\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the training output can be fit and transformed, and the testing output transformed accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = output_pipeline.fit_transform(y_train)\n",
    "y_test = output_pipeline.transform(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y_train[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the encoding of the output is now a sparse array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validation set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to make this reproducible, we have to seed the random number generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the training, we will require a validation set, so we split the training data into two sets, one for actual training, the other for validation.  Note, we don't touch the test data set at all during the training process.  The default is using 75 % of the data for training, 25 % for validation. This function will also shuffle the data set prior to splitting (hence seeding the random number generator)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train, x_val, y_train, y_val = train_test_split(x_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train.shape, x_val.shape, y_train.shape, y_val.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Verification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we've done quite some transformations, let's verify whether we didn't mess up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frame = plt.gca()\n",
    "frame.axes.get_xaxis().set_visible(False)\n",
    "frame.axes.get_yaxis().set_visible(False)\n",
    "plt.imshow(x_train[0].reshape(28, 28), cmap='gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y_train[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Everything seems fine, input/output are as expected, and we can start doing some machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classic multilayer perceptron (MLP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start off training a classic multilayer neural network to familiarie ourselves with the keras framework."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import Sequential\n",
    "from keras.layers.core import Dense, Dropout\n",
    "from keras.optimizers import SGD\n",
    "from sklearn.metrics import confusion_matrix\n",
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create a fully connected neural network with 784 input units (28 $\\times$ 28 pixels), two hidden layers with 512 units each, and an output layter with 10 units (ten categories, one per digit). For the two hidden layers, we will use a ReLU activation function, and a SoftMax for the output layer.\n",
    "\n",
    "To obtain repeatable results, we need to seed the random number generator, taking into account that TensorFlow uses its own."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1234)\n",
    "tf.set_random_seed(4958)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can can compile the model, specifying the loss function (categorical cross-entropy), the optimizer (SGD, Stochastic Gradient Descent), and the metrics (accuracy) we want to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.compile(loss='categorical_crossentropy', optimizer=SGD(), metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model_history = model.fit(x_train, y_train, batch_size=128, epochs=100,\n",
    "                          verbose=1, validation_data=(x_val, y_val))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the history of the training process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_history(network_history):\n",
    "    plt.figure()\n",
    "    plt.xlabel('Epochs')\n",
    "    plt.ylabel('Loss')\n",
    "    plt.plot(network_history.history['loss'])\n",
    "    plt.plot(network_history.history['val_loss'])\n",
    "    plt.legend(['Training', 'Validation'])\n",
    "\n",
    "    plt.figure()\n",
    "    plt.xlabel('Epochs')\n",
    "    plt.ylabel('Accuracy')\n",
    "    plt.plot(network_history.history['acc'])\n",
    "    plt.plot(network_history.history['val_acc'])\n",
    "    plt.legend(['Training', 'Validation'], loc='lower right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_history(model_history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the performance of the model on the training, validation and test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.evaluate(x_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.evaluate(x_val, y_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "model.evaluate(x_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can save the model as an HDF5 file so that it can later be reloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.save('mnist_mlp.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dropout layers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to reduce overfitting, drop out can be applied, i.e., randomly setting input values to 0 during training. We add a drop-out layer between the hidden layers, and between the last hidden layer and the output layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dropout_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dropout_model.compile(loss='categorical_crossentropy', optimizer=SGD(),\n",
    "                      metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dropout_model_history = dropout_model.fit(x_train, y_train, batch_size=128,\n",
    "                                          epochs=100, verbose=1,\n",
    "                                          validation_data=(x_val, y_val))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_history(dropout_model_history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, let's compare the performance of the model on the training, validation and test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dropout_model.evaluate(x_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dropout_model.evaluate(x_val, y_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dropout_model.evaluate(x_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dropout_model.save('mnist_mlp_dropout.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Understanding the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should try to gain some insight into the model.  There are several ways to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_model = model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace `model` by `dropout_model` in case you want to analyze the latter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import pandas as pd\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confusion matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model produces some classification errors, it would be interesting to see the type of errors.  Computing a confustion matrix is useful for that prupose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_predict_classes = analysis_model.predict_classes(x_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test_classes = np.argmax(y_test, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_test_classes, y_predict_classes)\n",
    "cm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the confusion matrix, it is clear that, e.g., 2 and 7 get confused, as well as 4 and 9.  Given the data, that should come as no surprise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm, classes,\n",
    "                          normalize=False,\n",
    "                          cmap=plt.cm.Blues):\n",
    "\n",
    "    log1p_cm = np.log1p(cm)\n",
    "    if normalize:\n",
    "        cm = cm.astype(np.float)/cm.sum(axis=1)[:, np.newaxis]\n",
    "    \n",
    "    figure, axes = plt.subplots(figsize=(6, 6))\n",
    "    axes.imshow(log1p_cm, interpolation='nearest', cmap=cmap)\n",
    "    axes.set_xticks(classes)\n",
    "    axes.set_yticks(classes)\n",
    "\n",
    "    fmt = '{0:.4f}' if normalize else '{0:d}'\n",
    "    thresh = 0.5*cm.max()\n",
    "    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n",
    "        axes.text(j, i, fmt.format(cm[i, j]),\n",
    "                  horizontalalignment=\"center\",\n",
    "                  color=\"white\" if cm[i, j] > thresh else \"black\",\n",
    "                  fontsize=8)\n",
    "\n",
    "    figure.tight_layout()\n",
    "    axes.set_ylabel('True label')\n",
    "    axes.set_xlabel('Predicted label')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_confusion_matrix(cm, range(10), normalize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sensitivity to initial conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How sensitive is the training process to the initial values of the model?  Let's train and evaluate the model several times, and observe the outcomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(1234)\n",
    "tf.set_random_seed(459845)\n",
    "names = ('train', 'val', 'test')\n",
    "inputs = (x_train, x_val, x_test)\n",
    "outputs = (y_train, y_val, y_test)\n",
    "losses = dict()\n",
    "accuracies = dict()\n",
    "for name in names:\n",
    "    losses[name] = []\n",
    "    accuracies[name] = []\n",
    "for i in range(10):\n",
    "    print(f'training model {i+ 1:d}')\n",
    "    mod = Sequential()\n",
    "    mod.add(Dense(512, activation='relu', input_shape=(784,)))\n",
    "    mod.add(Dropout(0.2))\n",
    "    mod.add(Dense(512, activation='relu'))\n",
    "    mod.add(Dropout(0.2))\n",
    "    mod.add(Dense(10, activation='softmax'))\n",
    "    mod.compile(loss='categorical_crossentropy', optimizer=SGD(),\n",
    "                metrics=['accuracy'])\n",
    "    _ = mod.fit(x_train, y_train, batch_size=128, epochs=100,\n",
    "                verbose=0, validation_data=(x_val, y_val))\n",
    "    for name, input, output in zip(names, inputs, outputs):\n",
    "        loss, accuracy = mod.evaluate(input, output)\n",
    "        losses[name].append(loss)\n",
    "        accuracies[name].append(accuracy)\n",
    "    mod.save(f'mnist_mpl_{i + 1:02d}.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Losses and accuracies are quite reproducible across runs with different initializations of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sns.stripplot(data=pd.DataFrame(losses), jitter=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.stripplot(data=pd.DataFrame(accuracies), jitter=0.1);"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}