LMB_Python_DNN_Part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deep Neural Networks Using the Keras Python Module\n",
    "\n",
    "This course introduces the basics of working with deep neural networks (DNNs) using the [Keras](https://keras.io/) Python module. It is intentional that this course does not go into the mathematics that underpin the function of the DNNs. Rather, the aim of the course is to provide practical experience of working with DNNs using high-level instructions, focusing on a general encoding, training and evaluation pipeline that can be applied in many subject areas.\n",
    "\n",
    "This course uses two examples, the first relates to the identification of handwritten numerals (a standardized dataset) and a second, more complex example, relating to the classification of biological sequences (proteins in this case). Various simple DNN architectures will be covered, but this course will not cover the 2D convolutional neural networks commonly employed in image analysis; this will be covered by a separate, dedicated course.\n",
    "\n",
    "To get a feel for how deep neural networks actually work we highly recommend Google Tensorflow's [Neural Network Playground](https://playground.tensorflow.org).  \n",
    "\n",
    "\n",
    "## The Keras module\n",
    "\n",
    "Here we make extensive use of the Keras library because it is a fairly simple, high-level means to implement DNNs, that sits on top of an efficient underlying machine-learning framework (here Google's TensorFlow). This provides fast computation, running on both a regular computer processor (CPU) or graphics cards (GPU), while at the same time being fairly easy to program for the most commonly used types of DNN.\n",
    "\n",
    "Installing the Keras Python module so that it works on CPUs is generally straightforward, for example it is easily installed on Linux systems using the `pip install keras` command. However, installing Keras for GPUs can be a little more complex as it involves having the correct hardware drivers and GPU toolkit (low-level API) installed for your graphics card (usually made by NVIDIA). For installing the GPU version, the general procedure is to install an up-to-date [NVIDA graphics driver](https://www.nvidia.com/drivers), and the [CUDA toolkit](https://developer.nvidia.com/cuda-downloads), being careful to make sure that your driver and CUDA versions are compatible. Then you can install the GPU version of [Tensorflow](https://www.tensorflow.org/); this is easily achieved on Linux using `pip install tensorflow-gpu` (see [here](https://www.tensorflow.org/install/pip) for further help). A final installation of Keras with `pip install keras` should then detect the GPU-enabled libraries.\n",
    "\n",
    "For the simple DNNs we use in this course the Keras the CPU-only version is generally sufficient and it will make good use of parallel processors. However, the GPU implementation is somewhat faster (caveat to having a modern, compatible graphics card) and makes training much easier in many real-world situations, and especially so with larger DNNs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The workflow\n",
    "\n",
    "This tutorial will spend some time to present various steps that you would implement in a typical DNN workflow. These may be roughly broken down as:\n",
    "\n",
    "* Preparing main training and test data\n",
    "* Constructing neural networks\n",
    "* Training neural networks to minimise loss\n",
    "* Avoiding over-training\n",
    "* Optimising and network architecture and parameters\n",
    "* Objectively assessing predictive performance\n",
    "\n",
    "\n",
    "<BR CLEAR=\"left\">\n",
    "<img src=\"images_ML/03_08.png\" style=\"width: 1000px;\" align=\"left\"/>\n",
    "<BR CLEAR=\"left\">\n",
    "\n",
    "### Glossary of terms\n",
    "\n",
    "**Training data**\n",
    "\n",
    "Training data is the information that is used to update the neural network (via its internal weights that determine the strength of inter-neuron connections) so that it performs some predictive function. For the examples used here, we will be performing supervised machine learning; we have data for which a correct answer is known. Hence, the training data will consist of inputs and the corresponding known, hopefully correct, answer for the outputs. By learning to connect inputs to known outputs the hope is that the network is then sufficiently general to make reliable predictions on unseen data.\n",
    "\n",
    "**Test data**\n",
    "\n",
    "Test data refers to data that is used to assess the predictive performance of the network for the purposes of optimising the network architecture and parameters. Test data, just like the main training data, consists of input data paired with the known, correct output data. However, training is separate from the directly used training data; the test data is not used to update DNN weights. However, the test data does influence the formulation of the DNNs predictive model, and so in the strictest sense it cannot be used to test the performance of the network in a completely objective way. \n",
    "\n",
    "**Validation data**\n",
    "\n",
    "Validation data also comprises of inputs paired with known outputs, but has not been used in any way to train the network or optimise it's architecture; it should be completely unseen by the network. Analysis of entirely separate validation data at the end is the only completely fair way to test the predictive performance of the network. During this course we will be lazy and not use properly separate validation data, but this is the wayyou should work for scientific papers.\n",
    "\n",
    "**Labels**\n",
    "\n",
    "Labels refer to categorical output data; used by a network for classification purposes. This relates to both the true, known category values that come paired with the training/test data and also the predicted outputs for previously unseen input.\n",
    "\n",
    "**Array shape**\n",
    "\n",
    "The shape of an array (coming from the NumPy module) refers to the sizes of the array's independent axes. For linear, 1D arrays this simply means their length. For higher-dimensional arrays, i.e. with more axes, this means the number of rows, columns and layers etc. For example an image array of 512 by 256 pixels with three colour channels (red, green, blue) would have a shape of `(512, 256, 3)`.\n",
    "\n",
    "**Loss**\n",
    "\n",
    "Loss represents how close the neural network output from the ideal output from the known training data, i.e. the error. This is scored using a loss function that is selected to be appropriate to the task at hand, e.g. for categorical output whether only one or multiple outputs can be true/selected.\n",
    "\n",
    "**Accuracy**\n",
    "\n",
    "Accuracy represents the proportion of predictions that are, to a certain precision, correct.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started - Handwriting recognition\n",
    "\n",
    "Our first practical examples involve recognising handwritten numerals. A selection of these is shown here:\n",
    "\n",
    "![MNIST examples](https://upload.wikimedia.org/wikipedia/commons/2/27/MnistExamples.png)\n",
    "\n",
    "As you might guess, the task we will require of the DNNs is to take a black-and-white image of a handwritten number and classify it as one of the ten discrete number classes, 0-9. This is an example set that is automatically provided with Keras. These data simplify the task of character recognition by having the training data pre-prepared; the numbers are scaled and centred on to a fixed size pixel grid and have normalised image intensities (going from white to black). In a real-world application you would probably want to use a 2D convolutional network for handwriting recognition as this kind of network is position and (somewhat) scale independent. However, using pre-aligned images means we can provide a tutorial that works with the simplest of neural network architectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing modules\n",
    "\n",
    "The first commands in this notebook perform some set-up for smoother running of the course. These control how the graphical content for graphs and charts will be displayed. Also, certain warnings, which are harmless but would otherwise be distracting are surppressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we perform some imports to enable access to the NumPy module (which deals with array data) and the MatPlotLib module which will handle the drawing of graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  # Load the NumPy module, assign it the name \"np\" for convienence\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (16,8) # Set plot size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last imports relate to the `Keras` library for machine learning. The first line below imports various sub-modules that will be used to construct the DNNs. The second line is the`mnist` example dataset of handwritten digits. And the last import is a handy function to help prepare the input data into a standard form compatible with the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras import models, layers, regularizers\n",
    "from keras.datasets import mnist\n",
    "from keras.utils import to_categorical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we actually load the example data and extract the pre-separated training and test data. Then for each of these two sections we separate the input (image) data from the output (numerical categories) data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, test = mnist.load_data()\n",
    "train_images, train_labels = train\n",
    "test_images, test_labels = test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can inspect the contents of the dataset using the `.imshow()` function of MatPlotLib. These data are greyscale images if 28 by 28 pixels, which corresponds to a total of 784 separate (though often correlated) input features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.imshow(train_images[1], cmap='Greys')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can use these data they must be adapted to work with a DNN that expects a 1D input vector and the input and output values to be floating point numbers in the range 0.0 to 1.0 (or sometimes -1.0 to 1.0), or at least approximately so."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, to prepare the input data we change the shapes of the data arrays to convert what is effectively a list of image matrices (each with rows and columns) into a list of flat, 1D vectors. The DNNs we use will sill be able to detect correlations between the pixels even if the data is rearranged because we well be using a fully connected network (as long as the same arrangement is preserved for all input data)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(train_images.shape)\n",
    "print(test_images.shape)\n",
    "n_train, w, h = train_images.shape\n",
    "n_test, w, h = test_images.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, the last two axes are combined into one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size = w * h\n",
    "train_images = train_images.reshape(n_train, size)\n",
    "test_images  = test_images.reshape(n_test, size)\n",
    "print(train_images.shape)\n",
    "print(test_images.shape)\n",
    "print(train_images.min(), train_images.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input image data is then scaled so that the pixel intensity values, which are initially in the range 0 to 255, are put into the range 0.0 to 1.0, i.e. to work at a scale that the DNNs are tuned to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(train_images.min(), train_images.max())\n",
    "\n",
    "train_images = train_images.astype(float)/255.0\n",
    "test_images  = test_images.astype(float)/255.0\n",
    "\n",
    "print(train_images.min(), train_images.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output data must also be prepared for use with the network. As we can see, these are a flat 1D array containing the numbers 0-9; i.e. labels representing number categories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(train_labels.shape)\n",
    "print(train_labels[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we want the output data to be in the range 0.0 to 1.0 (not 0 to 9). You might be tempted to simply scale the number labels. However, this will not work well, because we are dealing with discrete categories; the labels are separate and the \"0\" numeral is not closest to \"1\" in the handwriting sense.\n",
    "\n",
    "Instead what we do is convert the number categories into a binary matrix. Each output will be a length 10 vector full of mostly `0.0` with a `1.0` in one column indicating which (number) category is selected. Our output data is easily converted into this form using the `to_categorical()` function from Keras."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_labels = to_categorical(train_labels)\n",
    "test_labels = to_categorical(test_labels)\n",
    "print(train_labels.shape)\n",
    "print(train_labels[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding fully-connected neural network layers\n",
    "\n",
    "At last we come to a point where we actually start to construct a neural network. We will begin by constructing a somewhat old-fashioned three-layer feed-forward neural network. Firstly we construct the neural network framework using the \"Sequential\" model provided by Keras: this is an easy means for successively adding the networks layers, each of which is connected to the previous layer. The other \"functional\" way of constructing the NN allows more complex interconnections, but we will not need that here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn1 = models.Sequential()\n",
    "print(nn1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the framework for our first neural network `nn1` is initialised, we can then add its first layer. The first layer is a `Dense` layer, which means it is a fully connected one; all inputs connect to all nodes:\n",
    "This is the layer that will accept the flattened image data.\n",
    "\n",
    "We specify that the layer should have a specified number of neuron nodes (here `32`) and that it should expect input of length `size` (which came from image width times height). The `activation` parameter is set to `relu` ([rectified linear unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks)). Setting the activation means we are specifying how the neural nodes respond to the sum of their input weights. The `relu` option generates a linear response if positive and none when negative. This type of activation function is very commonly used with DNNs and enables better training of deep networks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layer = layers.Dense(32, input_shape=[size], activation='relu')\n",
    "nn1.add(layer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then add a second, output layer. Here we do not need to specify the data shape, as Keras can deduce this from the number of nodes in the previous layer. The layer has `10` nodes: one for each category of output. The `softmax` activation here is commonly used for binary categorical output (where one number should be 1.0 and the rest 0.0); it is a scaled exponential that acts as a sensitive trigger for positive values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn1.add(layers.Dense(10, activation='softmax'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object representing the neural network has a handy `.summary()` function that described the construction of its various layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(nn1.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should be noted that the size of the input vector is `784` fixed (by the size of the image) and that the size of the output, `10` is fixed by the number of categories. Also, the number of nodes diminishes between these two values (`256` then `32`): this funnels the decision-making nodes into fewer states. Setting an optimal number of nodes can be achieved by trying different values and seeing which work best. When tweaking the parameters in this way it can be quickly established that, beyond a certain number, adding more nodes does not improve training and slows things down. Having fewer nodes will improve speed, but too few will degrade predictive performance.\n",
    "\n",
    "Lastly we finalize the network's construction using the `.compile()` function. Here we are specifying three options. The first, `optimiser` specifies which method (from a standard set) will be used to improve the network weights during training. The `loss` option states how we score comparisons between the network's generated output and the actual known training output. This is generally set according to the type of output. Here `'categorical_crossentropy'` is the standard choice for categorical selection where only one category can be correct. The last `metrics` option simply specifies which parameters we would like to follow during the training procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn1.compile(optimizer='rmsprop', loss='categorical_crossentropy',\n",
    "            metrics=['accuracy'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loss functions\n",
    "\n",
    "There are many different loss functions that are availabe in Keras. Athough you are fairly free to try different ones, each is created for a particular kind of output and often works best with a specific activation function being present in the last layer.\n",
    "\n",
    "Some of the available loss functions include:\n",
    "\n",
    "* `binary_crossentropy` : this is used when is a single true/false output *or* if there are multiple true/false outputs where multiple outputs can be true at once. This is usually paired with a *sigmoidal* activation function in the output layer.\n",
    "\n",
    "* `categorical_crossentropy` : this is used for selecting one amongst multiple categorical classes, i.e. where only one output can be true. There should be an output node for each class. This is usually paired with a *softmax* output layer activation function.\n",
    "\n",
    "* `mean_squared_error` : this is often the first choice for regression problems. This can be used with a *linear* activation function if the range of outputs is unbounded or a *sigmoidal* activation function if bounded."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training \n",
    "\n",
    "Next we actually strain the netral network `nn1` by using the `.fit()` function on the training input and (correct) target training output. Here `epochs` relates to the number of repeat training iterations, over the whole training dataset and `batch_size` states how many input examples to sample each time the network is updated. The validation data is naturally the test input and test output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "history = nn1.fit(train_images, train_labels, epochs=10, batch_size=1024,\n",
    "                 validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the function call above captures `history` which records the accuracy values etc. during the training procedure. (We will inspect this later on)\n",
    "\n",
    "Although the training produces textual output to indicate how accurate the predictions are. The network can be evaluated at any time, using any paired input and output data, via the `.evaluate()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_loss, test_acc = nn1.evaluate(test_images, test_labels)\n",
    "\n",
    "msg = 'Acc: {:.2f}%'.format(test_acc*100.0)\n",
    "print(msg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully you will see an accuracy of around 94-95%, which is pretty good considering that out network only has two fairly narrow layers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the trained network we can also, obviously, use it to make predictions. Here we extract and display a random image from the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import randint\n",
    "idx = randint(0, len(test_images))\n",
    "img = test_images[idx]\n",
    "plt.imshow(img.reshape(w,h), cmap='Greys')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can run the tet data through the neural network with the `.predict()` function and see what the prediction for this image is. Note that the network output is categorical binary: a length ten vector where only one value will be near one, the rest zero. Hence, to extract the digit class we find the index of the maximum (i.e. selected) value with NumPy's `.argmax()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = nn1.predict(test_images)\n",
    "print('Binary out:', np.round(out[idx], 3))\n",
    "print('Digit out:', np.argmax(out[idx]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we return to the history data that we recorded when the neural network was trained. \n",
    "\n",
    "So that we can easily inspect the training history, a helper function `plot_training_histroy()` is provided below. This will take any number of history objects and plot graphs for how the loss (the mismatch between predictions and training output) and the accuraccy change during changing, for both the main training data and the test data (which was nut directly trained with)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "\n",
    "def plot_training_history(*histories):\n",
    "    cmap = cm.get_cmap('tab10')\n",
    "    fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "    plot_options = {'linewidth':2, 'alpha':0.8} # A dictionary of inputs for all charts\n",
    "    \n",
    "    for i, history in enumerate(histories):\n",
    "        hd = history.history\n",
    "        n = np.arange(len(hd['loss'])) + 1\n",
    "        plot_options['color'] = cmap(float(i % 10)/10)\n",
    "        \n",
    "        ax1.plot(n, hd['loss'], label='Train %d' % i,\n",
    "                 linestyle='--', **plot_options)\n",
    "        ax1.plot(n, hd['val_loss'], label='Test %d' % i,\n",
    "                 **plot_options)\n",
    "        ax1.set_title('Loss')\n",
    "        \n",
    "        ax2.plot(n, hd['accuracy'], label='Train %d' % i,\n",
    "                 linestyle='--', **plot_options)\n",
    "        ax2.plot(n, hd['val_accuracy'], label='Test %d' % i,\n",
    "                 **plot_options)\n",
    "        ax2.set_title('Accuracy')\n",
    "        \n",
    "    ax1.legend()    \n",
    "    ax2.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The history plotting function is then run on the `history` object generated earlier. You can see that the loss is reduced as the accuraccy increases, for both segments of data. Though the accuracy levels of at about 97-98% the training could be run for longer (e.g. more epochs) to get a slight improvement. However, while training for longer increases the accuracy for the training data it plateaus earlier for the test data. This is an indication that the network is becoming **over trained**, i.e. too specialised to the traiing data, and so lo longer general enough to give the same accuracy on unseen data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding layers\n",
    "\n",
    "Expanding upon this first network, we can see what happens if we introduce more layers. Construction is similar to before, albeit now with three layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn1 = models.Sequential()\n",
    "nn1.add(layers.Dense(256, activation='relu', input_shape=[size])) # New, wider input layer\n",
    "nn1.add(layers.Dense(32, activation='relu'))                      # Middle layer\n",
    "nn1.add(layers.Dense(10, activation='softmax'))\n",
    "nn1.compile(optimizer='rmsprop', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history1 = nn1.fit(train_images, train_labels, epochs=10, batch_size=1024,\n",
    "                 validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the three-layer network we get even better predictive performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history, history1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Layer widths\n",
    "\n",
    "Next we look at what happens when the number of nodes in the network layers are adjusted. Here we create a similar network and try using double the number of nodes in the middle layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn2 = models.Sequential()\n",
    "nn2.add(layers.Dense(256, activation='relu', input_shape=[size]))\n",
    "nn2.add(layers.Dense(64, activation='relu')) # Initially 32\n",
    "nn2.add(layers.Dense(10, activation='softmax'))\n",
    "nn2.compile(optimizer='rmsprop', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history2 = nn2.fit(train_images, train_labels, epochs=10, batch_size=1024,\n",
    "                 validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison of the history objects shows that although the loss minimisation is generally slightly better for the training data there is no real improvement in the accuracy of the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history1, history2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running with far fewer nodes trains faster, as we might expect, but shows a decrease in the training accuracy. Indeed this is not much better than the intial two-layer network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn2 = models.Sequential()\n",
    "nn2.add(layers.Dense(32, activation='relu', input_shape=[size])) # Was 256\n",
    "nn2.add(layers.Dense(32, activation='relu'))\n",
    "nn2.add(layers.Dense(10, activation='softmax'))\n",
    "nn2.compile(optimizer='rmsprop', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history2 = nn2.fit(train_images, train_labels, epochs=10, batch_size=1024,\n",
    "                 validation_data=(test_images, test_labels))\n",
    "\n",
    "plot_training_history(history1, history2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Over-training\n",
    "\n",
    "We will now run the first three-layer network again for more iterations (`epochs=20`) to better show the overtraining in the network; as the loss and accuracy for the training data diverges from the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn1 = models.Sequential()\n",
    "nn1.add(layers.Dense(256, activation='relu', input_shape=[size]))\n",
    "nn1.add(layers.Dense(32, activation='relu')) # Initially 32\n",
    "nn1.add(layers.Dense(10, activation='softmax'))\n",
    "nn1.compile(optimizer='rmsprop', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history1 = nn1.fit(train_images, train_labels, epochs=20, batch_size=1024,\n",
    "                 validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In later iterations (epoch > 10) there is a clear divergence in both the loss and accuracy of the test and compare to the training data. Also, we can see that the results for the test data are more noisy than the training data. Hence, in essence out network has become too specialised toward the training data and loses predictive performance on unseen data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Going  deep\n",
    "\n",
    "It may be tempting to improve the network by opting for a deeper architecture, with more layers. However, the deeper the network the more the capacity for overtraining and the slower and more difficult the training. Here for example, having three fairly wide internal layers does not improve predictions and indeed makes the overtraining slightly worse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn2 = models.Sequential()\n",
    "nn2.add(layers.Dense(256, activation='relu', input_shape=[size]))\n",
    "nn2.add(layers.Dense(128, activation='relu'))\n",
    "nn2.add(layers.Dense(128, activation='relu'))\n",
    "nn2.add(layers.Dense(128, activation='relu'))\n",
    "nn2.add(layers.Dense(10, activation='softmax'))\n",
    "nn2.compile(optimizer='rmsprop', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history2 = nn2.fit(train_images, train_labels, epochs=20, batch_size=1024,\n",
    "                   validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history1, history2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Activations\n",
    "\n",
    "As is hopefully becoming obvious, there are a vast number of tweaks that can be made to the DNN architecture. However, it shouldn't take too long to see what is to complex or too simple as long as we take care to analyse the training history. Similarly, we can test the other high-level parameters. For example here we test the 'sigmoid' activation function: historically this was often the favoured option for three-layer networks before training deep networks became practical.\n",
    "\n",
    "As we can see this activation function trains the DNN somewhat more slowly than `relu`, and in the end (if we train for further iterations) has a worse accuracy. Interestingly, however, the overtraining is less severe; this is probably because the sigmoidal function can attenuate large weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn3 = models.Sequential()\n",
    "nn3.add(layers.Dense(256, activation='sigmoid', input_shape=[size]))\n",
    "nn3.add(layers.Dense(32, activation='sigmoid'))\n",
    "nn3.add(layers.Dense(10, activation='softmax'))\n",
    "nn3.compile(optimizer='rmsprop', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history3 = nn3.fit(train_images, train_labels, epochs=20, batch_size=1048,\n",
    "                   validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history2, history3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimisers\n",
    "\n",
    "We can also play with different training optimisers: the routines which adjust the network weights (in a back-propagation manner) according to internal gradients. When choosing an optimiser we not only have to think about the prediction accuracy, but also the convergence efficiency and the stability; sometimes the DNN weights can bounce around rather than converge smoothly. Here we test the `adam` optimiser and show that it works fairly well with fast convergence and good accuracy. Indeed, this is often a good intial one to try for many fully-connected DNNs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn3 = models.Sequential()\n",
    "nn3.add(layers.Dense(256, activation='relu', input_shape=[size]))\n",
    "nn3.add(layers.Dense(32, activation='relu'))\n",
    "nn3.add(layers.Dense(10, activation='softmax'))\n",
    "\n",
    "opt = 'adam' #'sgd' is bad\n",
    "nn3.compile(optimizer=opt, loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history3 = nn3.fit(train_images, train_labels, epochs=20, batch_size=1048,\n",
    "                   validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this new optimiser we can see that training has become a bit more efficient and a little smoother. There is a perhaps a small improvement in accuracy, but overtraining is still evident. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history2, history3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regularisation\n",
    "\n",
    "One way of reducing over-training is to add a regulariser to the layers. This encourages the network to adopt a more general set of internal weights by penalising weights from becoming too large; large weights are an indication of specialisation to specific inputs. Here we add an 'L2' regulariser which penalises based on the square of the weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg = regularizers.l2(0.001)\n",
    "\n",
    "nn4 = models.Sequential()\n",
    "nn4.add(layers.Dense(256, activation='relu', input_shape=[size], kernel_regularizer=reg))\n",
    "nn4.add(layers.Dense(32, activation='relu', kernel_regularizer=reg))\n",
    "nn4.add(layers.Dense(10, activation='softmax'))\n",
    "nn4.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history4 = nn4.fit(train_images, train_labels, epochs=20, batch_size=1024,\n",
    "                   validation_data=(test_images, test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see this does help with overtraining and doesn't affect test accuracy much (if we can train for longer)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_training_history(history3, history4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dropout\n",
    "\n",
    "A generally better way to avoid overtraining is by using a technique called *dropout*. What this does is to temporarily remove a random selection of weights (and hence input nodes) from the networks inter-layer connections. This prevents the weights becoming too specialised toward the training data; as random connections in the network are omitted the remaining connections are forced toward sparser, more general/average solutions with less co-dependency between particular nodes. \n",
    "\n",
    "Within Keras this is simply achieved by adding extra `Dropout` layers, and for each specifying the proportion of input connections to randomly omit (values from 0.2 to 0.4 is generally a good starting point). As we can see, this technique slows training but removes the over-training and maintains a good test accuracy. Indeed, the network may be a little under-trained (so we could train for longer)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn4 = models.Sequential()\n",
    "nn4.add(layers.Dense(256, activation='relu', input_shape=[size]))\n",
    "nn4.add(layers.Dropout(0.4))\n",
    "nn4.add(layers.Dense(32, activation='relu'))\n",
    "nn4.add(layers.Dropout(0.4))\n",
    "nn4.add(layers.Dense(10, activation='softmax'))\n",
    "nn4.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "history4 = nn4.fit(train_images, train_labels, epochs=40, batch_size=1024,\n",
    "                   validation_data=(test_images, test_labels))\n",
    "\n",
    "plot_training_history(history3, history4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outlier analysis\n",
    "\n",
    "After making a variety of tweaks that improve the DNN accuracy and reduce overtraining we are consistently reaching a test set accuracy of around 98%. We may be fairly satisfied with the result. However, we may learn something if we investigate the images for which we still cannot make good predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss, acc = nn4.evaluate(test_images, test_labels)\n",
    "print('Loss: {:.3f} Accuracy:{:.3f}'.format(loss, acc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do this we need to find the images where the DNN output prediction has a mismatch to the known labels. We can do this by converting the binary, categorical arrays to the index of the highest-scoring/largest value and test whether these are the same or different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_pred = nn4.predict(test_images)\n",
    "out_nums = out_pred.argmax(axis=1)\n",
    "known_nums = test_labels.argmax(axis=1)\n",
    "print(out_nums[:10])\n",
    "print(known_nums[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing whether the outputs differ involves the `!=` test (not equals), so we can pull out the indices where they are indeed different using NumPy's `.nonzero()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "known_nums = test_labels.argmax(axis=1)\n",
    "different = out_nums != known_nums\n",
    "bad = different.nonzero()[0]  # Indices of True values: where arrays are different\n",
    "print(bad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then inspect some of the images where the known and predicted digit differ. To the human eye, some of these errors might be deemed understandable, but we might to improve others by adjusting the DNN (perhaps moving to a convolutional architecture)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10 # Limit to first five errors\n",
    "indices = bad[:n] \n",
    "fig, axarr = plt.subplots(1, n) # n columns of images\n",
    "\n",
    "for a, idx in enumerate(indices):\n",
    "    img = test_images[idx,:].reshape(w,h)\n",
    "    p = out_nums[idx]\n",
    "    k = known_nums[idx]\n",
    "    axarr[a].imshow(img, cmap='Greys')\n",
    "    axarr[a].set_title('Pred: {} True: {}'.format(p, k))\n",
    "    \n",
    "plt.show()"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}