Exercise-09.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MNIST handwritten digits classification with MLPs\n",
    "\n",
    "In this notebook, we'll train a multi-layer perceptron model to classify MNIST digits using [TensorFlow](https://www.tensorflow.org/) (version $\\ge$ 2.0 required) with the [Keras API](https://www.tensorflow.org/guide/keras/overview).\n",
    "\n",
    "First, the needed imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from pml_utils import show_failures\n",
    "\n",
    "import tensorflow as tf\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, Activation, Dropout, Flatten\n",
    "from tensorflow.keras.utils import plot_model, to_categorical\n",
    "\n",
    "from IPython.display import SVG, display\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "\n",
    "print('Using Tensorflow version: {}, and Keras version: {}.'.format(tf.__version__, tf.keras.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check if we have GPU available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if tf.test.is_gpu_available():\n",
    "    from tensorflow.python.client import device_lib\n",
    "    for d in device_lib.list_local_devices():\n",
    "        if d.device_type == 'GPU':\n",
    "            print('GPU', d.physical_device_desc)\n",
    "else:\n",
    "    print('No GPU, using CPU instead.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MNIST data set\n",
    "\n",
    "Next we'll load the MNIST handwritten digits data set using TensorFlow's own tools.  First time we may have to download the data, which can take a while.\n",
    "\n",
    "#### Altenative: Fashion-MNIST\n",
    "\n",
    "Alternatively, MNIST can be replaced with Fashion-MNIST, which can be used as drop-in replacement for MNIST.   Fashion-MNIST contains images of 10 fashion categories:\n",
    "\n",
    "Label|Description|Label|Description\n",
    "--- | --- |--- | ---\n",
    "0|T-shirt/top|5|Sandal\n",
    "1|Trouser|6|Shirt\n",
    "2|Pullover|7|Sneaker\n",
    "3|Dress|8|Bag\n",
    "4|Coat|9|Ankle boot\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras.datasets import mnist, fashion_mnist\n",
    "\n",
    "## MNIST:\n",
    "(X_train, y_train), (X_test, y_test) = mnist.load_data()\n",
    "## Fashion-MNIST:\n",
    "#(X_train, y_train), (X_test, y_test) = fashion_mnist.load_data()\n",
    "\n",
    "nb_classes = 10\n",
    "\n",
    "X_train = X_train.astype('float32')\n",
    "X_test = X_test.astype('float32')\n",
    "X_train /= 255.0\n",
    "X_test /= 255.0\n",
    "\n",
    "# one-hot encoding:\n",
    "Y_train = to_categorical(y_train, nb_classes)\n",
    "Y_test = to_categorical(y_test, nb_classes)\n",
    "\n",
    "print()\n",
    "print('MNIST data loaded: train:',len(X_train),'test:',len(X_test))\n",
    "print('X_train:', X_train.shape)\n",
    "print('y_train:', y_train.shape)\n",
    "print('Y_train:', Y_train.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training data (`X_train`) is a 3rd-order tensor of size (60000, 28, 28), i.e. it consists of 60000 images of size 28x28 pixels. `y_train` is a 60000-dimensional vector containing the correct classes (\"0\", \"1\", ..., \"9\") for each training sample, and `Y_train` is a [one-hot](https://en.wikipedia.org/wiki/One-hot) encoding of `y_train`.\n",
    "\n",
    "Let's take a closer look. Here are the first 10 training digits (or fashion items for Fashion-MNIST):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pltsize=1\n",
    "plt.figure(figsize=(10*pltsize, pltsize))\n",
    "\n",
    "for i in range(10):\n",
    "    plt.subplot(1,10,i+1)\n",
    "    plt.axis('off')\n",
    "    plt.imshow(X_train[i,:,:], cmap=\"gray\")\n",
    "    plt.title('Class: '+str(y_train[i]))\n",
    "    print('Training sample',i,': class:',y_train[i], ', one-hot encoded:', Y_train[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-layer perceptron (MLP) network\n",
    "\n",
    "### Activation functions\n",
    "\n",
    "Let's start by plotting some common activation functions for neural networks. `'relu'` stands for rectified linear unit, $y=\\max(0,x)$, a very simple non-linearity we will be using in our MLP network below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(-4,4,.01)\n",
    "plt.figure()\n",
    "plt.plot(x, np.maximum(x,0), label='relu')\n",
    "plt.plot(x, 1/(1+np.exp(-x)), label='sigmoid')\n",
    "plt.plot(x, np.tanh(x), label='tanh')\n",
    "plt.axis([-4, 4, -1.1, 1.5])\n",
    "plt.title('Activation functions')\n",
    "plt.legend(loc='best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialization\n",
    "\n",
    "Let's now create an MLP model that has multiple layers, non-linear activation functions, and optionally dropout layers for regularization.\n",
    "\n",
    "We first initialize the model with `Sequential()`.  Then we add a `Dense` layer that has 28*28=784 input nodes (one for each pixel in the input image) and 20 output nodes. The `Dense` layer connects each input to each output with some weight parameter. \n",
    "\n",
    "Next, the output of the dense layer is passed through a ReLU non-linear activation function.\n",
    "\n",
    "Commented out is an alternative, more complex, model that you can also try out.  It uses more layers and dropout.  `Dropout()` randomly sets a fraction of inputs to zero during training, which is one approach to regularization and can sometimes help to prevent overfitting.\n",
    "\n",
    "The output of the last layer needs to be a softmaxed 10-dimensional vector to match the groundtruth (`Y_train`).  This means that it will output 10 values between 0 and 1 which sum to 1, hence, together they can be interpreted as a probability distribution over our 10 classes.\n",
    "\n",
    "Finally, we select *categorical crossentropy* as the loss function, select [*Adam*](https://keras.io/optimizers/#adam) as the optimizer, add *accuracy* to the list of metrics to be evaluated, and `compile()` the model.  Adam is simply a an advanced version of stochastic gradient descent, note there are [several different options](https://keras.io/optimizers/) for the optimizer in Keras that we could use instead of *adam*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model initialization:\n",
    "model = Sequential()\n",
    "\n",
    "# A simple model:\n",
    "model.add(Dense(units=20, input_dim=28*28))\n",
    "model.add(Activation('relu'))\n",
    "\n",
    "# A bit more complex model:\n",
    "#model.add(Dense(units=50, input_dim=28*28))\n",
    "#model.add(Activation('relu'))\n",
    "#model.add(Dropout(0.2))\n",
    "\n",
    "#model.add(Dense(units=50))\n",
    "#model.add(Activation('relu'))\n",
    "#model.add(Dropout(0.2))\n",
    "\n",
    "# The last layer needs to be like this:\n",
    "model.add(Dense(units=10, activation='softmax'))\n",
    "\n",
    "model.compile(loss='categorical_crossentropy', \n",
    "              optimizer='adam', \n",
    "              metrics=['accuracy'])\n",
    "print(model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The summary shows that there are 15,910 parameters in total in our model.\n",
    "\n",
    "For example for the first dense layer we have 785x20 = 15,700 parameters as the weight matrix is of size 785x20 (not 784, as there's an additional bias term).\n",
    "\n",
    "We can also draw a fancier graph of our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_model(model, show_shapes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning\n",
    "\n",
    "Next, we'll train our model.  Notice how the interface is similar to scikit-learn: we still call the `fit()` method on our model object.\n",
    "\n",
    "An *epoch* means one pass through the whole training data, we'll begin by running training for 10 epochs.\n",
    "\n",
    "The `reshape()` function flattens our 28x28 images into vectors of length 784.  (This means we are not using any information about the spatial neighborhood relations of pixels.  This setup is known as the *permutation invariant MNIST*.)  \n",
    "\n",
    "You can run code below multiple times and it will continue the training process from where it left off.  If you want to start from scratch, re-initialize the model using the code a few cells ago. \n",
    "\n",
    "We use a batch size of 32, so the actual input will be 32x784 for each batch of 32 images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "epochs = 10\n",
    "\n",
    "history = model.fit(X_train.reshape((-1,28*28)), \n",
    "                    Y_train, \n",
    "                    epochs=epochs, \n",
    "                    batch_size=32,\n",
    "                    verbose=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now see how the training progressed. \n",
    "\n",
    "* *Loss* is a function of the difference of the network output and the target values.  We are minimizing the loss function during training so it should decrease over time.\n",
    "* *Accuracy* is the classification accuracy for the training data.  It gives some indication of the real accuracy of the model but cannot be fully trusted, as it may have overfitted and just memorizes the training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,3))\n",
    "plt.plot(history.epoch,history.history['loss'])\n",
    "plt.title('loss')\n",
    "\n",
    "plt.figure(figsize=(5,3))\n",
    "plt.plot(history.epoch,history.history['accuracy'])\n",
    "plt.title('accuracy');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference\n",
    "\n",
    "For a better measure of the quality of the model, let's see the model accuracy for the test data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "scores = model.evaluate(X_test.reshape((-1,28*28)), Y_test, verbose=2)\n",
    "print(\"%s: %.2f%%\" % (model.metrics_names[1], scores[1]*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now take a closer look at the results using the `show_failures()` helper function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the first 10 test digits the MLP classified to a wrong class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = model.predict(X_test.reshape((-1,28*28)))\n",
    "\n",
    "show_failures(predictions, y_test, X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `show_failures()` to inspect failures in more detail. For example, here are failures in which the true class was \"6\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_failures(predictions, y_test, X_test, trueclass=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute the confusion matrix to see which digits get mixed the most, and look at classification accuracies separately for each class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "print('Confusion matrix (rows: true classes; columns: predicted classes):'); print()\n",
    "cm=confusion_matrix(y_test, np.argmax(predictions, axis=1), labels=list(range(10)))\n",
    "print(cm); print()\n",
    "\n",
    "print('Classification accuracy for each class:'); print()\n",
    "for i,j in enumerate(cm.diagonal()/cm.sum(axis=1)): print(\"%d: %.4f\" % (i,j))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model tuning\n",
    "\n",
    "Modify the MLP model.  Try to improve the classification accuracy, or experiment with the effects of different parameters.  If you are interested in the state-of-the-art performance on permutation invariant MNIST, see e.g. this [recent paper](https://arxiv.org/abs/1507.02672) by Aalto University / The Curious AI Company researchers.\n",
    "\n",
    "You can also consult the Keras documentation at https://keras.io/.  For example, the Dense, Activation, and Dropout layers are described at https://keras.io/layers/core/."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}