neural-networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural networks with keras"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "---\n",
    "\n",
    "The last decade or so has seen an explosion in interest in neural networks, spurred on by huge successes that various neural network architectures have enjoyed across a number of classification tasks. This notebook shows how to use the Keras neural network library from within R to fit various kinds of neural networks. Keras is a high level library for specifying and training neural networks. It provides a relatively user-friendly way to interact with 3 powerful Python neural network libraries (`Theano`, `TensorFlow`, and `CNTK`) that provide efficient (low-level) tools for working with neural networks. \n",
    "\n",
    "The goal of this notebook is to get you up and running with Keras on a couple of \"normal\" prediction problems. In later notebooks we'll look at some extensions like tricks to boost performance and specific architectures - convolutional neural networks - that have are particularly useful for certain kinds of prediction tasks.\n",
    "\n",
    "This notebook is meant to be read after going through three spreadsheet-based explanations of three core concepts in neural networks.\n",
    "\n",
    "1. **Layers**: Neural networks are alternating sequences of linear algebra (matrix multiplication) and non-linear transformations. Each pair of matrix multiplication and non-linear transformation of the result constitutes a layer of the network.\n",
    "2. **Backpropagation**: a way to train the weight parameters of a neural network, involving passing inputs \"forward\" through the network to get predictions and hence errors, and then passing these errors \"backwards\" through the network to adjust the parameters. This happens multiple times until some stopping criterion is met.\n",
    "3. **Stochastic gradient descent** (or other optimization algorithms): This specifies precisely how errors should be used to update the weights during the backwards pass of the backpropagation algorithm. \n",
    "\n",
    "These are contained as worksheets in the accompanying Excel file *\"neural-networks-on-a-spreadsheet.xlsx\"*.\n",
    "\n",
    "#### Sources and acknowledgements\n",
    "\n",
    "There are a lot of excellent sources of neural network information on the web. Some sources that I have used for this notebook are:\n",
    "\n",
    "* Keras interface to RStudio: https://keras.rstudio.com/ (the MNIST example borrows heavily from here).\n",
    "* Keras documentation: https://keras.io/\n",
    "* The [fast.ai](http://www.fast.ai/) neural network course is really excellent and I have pretty shamelessly used their spreadsheet examples of the [matrix multiplication step](https://www.youtube.com/watch?v=e3aM6XTekJc&feature=youtu.be&t=3810) and [SGD and related optimizers](https://www.youtube.com/watch?v=V2h3IOBDvrA&feature=youtu.be&t=688)\n",
    "* Backpropagation is nicely explained [here](https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting your own installation of Keras working in R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We're starting by working on a web server - an Amazon Web Services EC2 \"instance\" - that already has everything you need pre-installed. This is just to save setup time during the workshop. Later we'll see how you can set up your own AWS instance, so you can keep working this way after the workshop is finished. But you'll probably also want to get Keras installed locally on your own laptop. This section contains a few pointers on how to do that.\n",
    "\n",
    "> A possible source of confusion is that the \"main\" (original) Keras is a Python library. Until quite recently, if you wanted to use Keras, you had to use Python. What we'll be using is an R package *also* called **keras**. This essentially allows you to access the \"main\" (Python) Keras from within R, without using Python.\n",
    "\n",
    "The **keras** R package creates a virtual (Python) environment - a named, isolated, working copy of Python with its own files, directories, and paths so that you can work with Python or Python libraries like Theano and Tensorflow from within R. The **keras** package will create the virtual environment for you, but it does mean that you need to have a working version of Python already installed.\n",
    "\n",
    "Supported versions include:\n",
    "* Anaconda Python (Recommended): https://www.continuum.io/downloads\n",
    "* Python Software Foundation: https://www.python.org/downloads/release/python-353/\n",
    "\n",
    "Once you have a Python installation, you can do the rest from R/RStudio. See the instructions [here](https://keras.rstudio.com/). You need to (within R/RStudio):\n",
    "\n",
    "1. install the keras R package from CRAN as follows:\n",
    "```\n",
    "install.packages(\"keras\")\n",
    "```\n",
    "2. Install both the core Keras library as well as the TensorFlow backend:\n",
    "```\n",
    "library(keras)\n",
    "install_keras()\n",
    "```\n",
    "\n",
    "This will provide you with default CPU-based installations of Keras and TensorFlow. If you want a more customized installation, e.g. if you want to take advantage of NVIDIA GPUs, see the documentation for `install_keras()`. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Revisiting predicting the presence of *Aloe dichotoma*\n",
    "\n",
    "In this section we revisit an example we looked in the previous sessions. In that notebook we built a classifier using various tree-based methods to predict either presence/absence or abundance of *Aloe dichotoma*. Here we'll do the same using a few varieties of feed-forward neural network. \n",
    "\n",
    "First we load the packages and data we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sys.setenv(KERAS_BACKEND = \"tensorflow\")\n",
    "library(keras)\n",
    "library(tidyverse)\n",
    "\n",
    "load(\"data/aloe.RData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first need to split our data along two \"dimensions\": \n",
    "\n",
    "* We need to divide our data into a training data set and a test data set. This is already done but usually you would need to do this here.\n",
    "* We need to separate the predictor variables/features (the **X** matrix) from the outcome variable (the **y** vector)\n",
    "\n",
    "Here we first split the training data frame into features `x_train` and outcome `y_train`, converting both to a matrix. The `x_train` matrix has 20 columns (one for each feature), while `y_train` is a single column vector containing the response. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_train <- as.matrix(aloe_pa %>% filter(train_id != 3) %>% select(-train_id, -present)) \n",
    "y_train <- as.matrix(aloe_pa %>% filter(train_id != 3) %>% select(present), ncol = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_test <- as.matrix(aloe_pa %>% filter(train_id == 3) %>% select(-train_id, -present))\n",
    "y_test <- as.matrix(aloe_pa %>% filter(train_id == 3) %>% select(present), ncol = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# check dimensions\n",
    "dim(x_train)\n",
    "dim(y_train)\n",
    "dim(x_test)\n",
    "dim(y_test)\n",
    "y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scaling the predictors\n",
    "\n",
    "In neural networks variables with very different scales can cause problems for the gradient calculations used in stochastic gradient descent, causing it either to slow down or fail to improve. Its always a good idea to standardise each predictor -- the easiest way to do this is by mean-variance standardiation i.e. making the mean of each column equal to zero and its variance equal to one. \n",
    "\n",
    "Note that we standardise both the training and the test data using the means and standard deviations calculated *on the training data*. That is, we do not independently standardise the test set. This is because we build the model based on the standardised training data, so the test set should be directly comparable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# means of predictors in training set \n",
    "means <- apply(x_train, 2, mean)\n",
    "# standard deviations of predictors in training set \n",
    "sds <- apply(x_train, 2, sd)\n",
    "# standardisation of training data by mean and std dev\n",
    "x_train <- t((t(x_train) - means)/ sds)\n",
    "# standardisation of test data by mean and std dev (note OF TRAINING DATA)\n",
    "x_test <- t((t(x_test) - means)/ sds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic feed-forward neural network\n",
    "\n",
    "We'll start by creating a simple 3-layer, fully connected network. This network has a single hidden layer, in addition to the input and output layers. We'll then add some more complexity by adding layers. In each case we use a very similar model building and testing strategy:\n",
    "\n",
    "* Create the model\n",
    "* Define and compile the model\n",
    "* Train/fit the model\n",
    "* Evaluate the model\n",
    "\n",
    "We go through each of these steps on more detail below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create the model\n",
    "\n",
    "A Keras \"model\" is a way to organize layers. We use the \"sequential\" model, in which each layer of the model follows sequentially from the previous one, so that we build up a linear stack of layers using the pipe (%>%) operator. The other kind of Keras model is a \"functional\" model, which is useful for more complex cases and which we won't cover here.\n",
    "\n",
    "We begin by creating an empty sequential model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model <- keras_model_sequential()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the model\n",
    "\n",
    "We define a model by sequentially adding layers. Each layer consists of \n",
    "\n",
    "* a way of connecting up the neurons from the previous layer with neurons in this layer,\n",
    "* an activation function, which takes as input the pre-activation (coming from the input neurons or the previous layer of neurons) and calculates an output. This can be interpreted as the rate of how much a neuron is firing.  The simplest form of such a function is that of a binary function, where the neuron is either firing or not.\n",
    "\n",
    "In Keras, you can specify both of these parts together, or you can separate them out. Thus below\n",
    "\n",
    "```\n",
    "layer_dense(units = 128, input_shape = c(20)) %>% layer_activation('relu')\n",
    "``` \n",
    "\n",
    "is entirely equivalent to \n",
    "\n",
    "```\n",
    "layer_dense(units = 128, input_shape = c(20), activation = 'relu')\n",
    "``` \n",
    "\n",
    "Below we define a model with a single hidden layer i.e. a three-layer feed-forward neural network. Note that we don't need to reassign the piped result back to `model` (i.e. `model <- model %>% ...`). This happens automatically. Take care with this - if you want to start with a fresh model you need to do this with a new call to `keras_model_sequential()`!\n",
    "\n",
    "The first layer we add specifies both the shape of the input data and the number of neurons we want in the first hidden layer (128). The shape of the input data is specified in the `input_shape` argument i.e. a length 20 numeric vector representing the tf-idf values of the tweet. A dense layer (`layer_dense`) specifies a fully connected layer, in which each of these 20 input neurons will be connected up with each of the 128 neurons in the hidden layer. This gives 2560 parameters to be estimated - with bias terms for each hidden neuron, this makes 2688 parameters. \n",
    "\n",
    "We use a rectified linear unit or \"relu\" activation function to non-linearly transform the pre-activation values (i.e. the weighted sum of inputs to each neuron in this layer). The role of the activation function is to add non-linearity to the neural network and hence provide the flexibility that give neural networks their predictive power (otherwise we'd just have a series of matrix multiplications that we could reduce to a single one). You can read more about activations functions [here](https://www.quora.com/What-is-the-role-of-the-activation-function-in-a-neural-network-How-does-this-function-in-a-human-neural-network-system/answer/Sebastian-Raschka-1). The plot below shows a few popular activation functions:\n",
    "\n",
    "<img src=\"img/activations.png\" />\n",
    "\n",
    "Together, the call to `layer_dense()` and `layer_activation()` define the hidden layer.\n",
    "\n",
    "The only other layer outputs a length 1 numeric vector (probabilities for the \"present\" class) using a sigmoid activation function. Each of the 128 neurons in the hidden layer are connected to the output neuron, giving another 129 parameters to be estimated (including bias). We end by using the `summary()` function to print the details of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% \n",
    "  layer_dense(units = 128,                  # number of neurons in the hidden layer\n",
    "              input_shape = c(20)) %>%     # dimension of input array\n",
    "  layer_activation('relu') %>%             # use a rectified linear unit as an activation function in the hidden layer\n",
    "  layer_dense(units = 1) %>%               # adds an output layer to the network\n",
    "  layer_activation('sigmoid')              # use sigmoid activation function in the output layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "summary(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compile the model \n",
    "\n",
    "Keras is built on top of Theano and TensorFlow, both packages that take the symbolic *computational graph* of the model and compile it in an efficient implementation that can then be used to train and evaluate the model. \n",
    "\n",
    "When compiling a model, Keras asks you to specify your **loss function** and your **optimizer**.  Here we use logarithmic loss as the loss function, which for two-class problems is called *binary crossentopy* by Keras (*categorical crossentropy* for multi-class classification problems). Cross-entropy is a measure of how different two probability distributions are. In this case we have two distributions (for each observation), a predicted distribution and a target distribution. \n",
    "\n",
    "* Our predictions are a probability distribution across the two classes (e.g. \"we're 80% confident the site contains an aloe (present) and 20% sure it doesn't (absent)\").\n",
    "* The target is a probability distribution with 100% for the correct category, and 0 for the other category. \n",
    "\n",
    "For example, say we have a site we know to be an absence, and that our model allocates a probability of 0.8 to that class. Then:\n",
    "* observed \"probability\" for class 1 = 100%, for class 2 = 0%\n",
    "* predicted probability for class 1 = 80%, predicted probability for class 2 = 20%\n",
    "* log loss $= - (1\\times\\log(0.8) + 0\\times\\log(0.2))$\n",
    "\n",
    "Total log loss would be calculated by summing over all observations in the training set.\n",
    "\n",
    "The optimizer helps determine how quickly the model learns, how resistent it is to getting \"stuck\" or \"blowing up\".  The chosen optimization algorithm will help to find optimized values for the weights which minimizes the error (or maximizes the accuracy). Popular choices for the optimizer include stochastic gradient descent, Adam and RMSprop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% compile(\n",
    "  optimizer = 'adam',\n",
    "  loss = 'binary_crossentropy',\n",
    "  metrics = c('accuracy')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have done the define and compile steps together, using the pipe, if we wanted to:\n",
    "\n",
    "```\n",
    "model %>% \n",
    "  layer_dense(units = 128, activation = 'relu', input_shape = c(20)) %>% \n",
    "  layer_dense(units = 1, activation = 'sigmoid') %>% \n",
    "  compile(\n",
    "    optimizer = 'SGD',\n",
    "    loss = 'binary_crossentropy',\n",
    "    metrics = c('accuracy')\n",
    "  )\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the model\n",
    "\n",
    "Next we train the model on our training data. We use **batch training** with a batch size of 32.  This approach sends groups of 32 observations through the network at a time.  Adjustments to the weights and biases (deltas) from each observation are accumulated and updated once after each batch. Watch the loss, which is the objective function which we are minimising, and the estimated accuracy of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% fit(x_train, y_train, epochs = 20, batch_size = 32) %>% plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note if we want to run for (say) 50 *further* epochs, we just make another call to `fit()`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% fit(x_train, y_train, epochs = 50, batch_size = 32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluate the model \n",
    "\n",
    "Once the model is trained, we can evaluate its performance on the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% evaluate(x_test, y_test, batch_size=32, verbose = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate predictions on new data (e.g. test data)\n",
    "\n",
    "This is done in a straightforward way by passing new data to `predict_classes()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% predict_classes(x_test) %>% head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding more layers = Deep learning\n",
    "\n",
    "Adding additional layers is straightforward in Keras. In this secton we build a model with two (rather than one) hidden layers. This model has the same number of hidden nodes as the previous model, half in each of the hidden layers. We use the same steps as before, so just give the headings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model2 <- keras_model_sequential()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define and compile the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model2 %>% \n",
    "  layer_dense(units = 64, activation = 'relu', input_shape = c(20)) %>% \n",
    "  layer_dense(units = 64, activation = 'relu') %>% \n",
    "  layer_dense(units = 1, activation = 'sigmoid') %>% \n",
    "  compile(\n",
    "    loss = 'binary_crossentropy',\n",
    "    optimizer = 'rmsprop',\n",
    "    metrics = c('accuracy')\n",
    "  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "summary(model2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model2 %>% fit(x_train, y_train, epochs = 25, batch_size = 32) %>% plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluate the model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model2 %>% evaluate(x_test, y_test, batch_size=32, verbose = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding dropout\n",
    "\n",
    "Dropout is a technique for avoiding overfitting in neural networks. The idea is that randomly selected neurons (together with their connections) are ignored during training. In each layer, we specify the proportion of neurons to exclude - in the example below we use a dropout rate of 20%. You can read more about dropout (although that is pretty much all there is to it!) in the original paper [here](http://www.jmlr.org/papers/volume15/srivastava14a/srivastava14a.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model3 <- keras_model_sequential()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define and compile the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model3 %>% \n",
    "  layer_dense(units = 64, activation = 'relu', input_shape = c(20)) %>%\n",
    "  layer_dropout(rate = 0.2) %>% \n",
    "  layer_dense(units = 64, activation = 'relu') %>% \n",
    "  layer_dropout(rate = 0.2) %>% \n",
    "  layer_dense(units = 1, activation = 'sigmoid') %>% \n",
    "  compile(\n",
    "    loss = 'binary_crossentropy',\n",
    "    optimizer = 'rmsprop',\n",
    "    metrics = c('accuracy')\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model3 %>% fit(x_train, y_train, epochs = 25, batch_size = 32) %>% plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluate the model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model3 %>% evaluate(x_test, y_test, batch_size=32, verbose = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another tutorial: the MNIST dataset\n",
    "\n",
    "No introduction to neural networks is complete without the MNIST dataset! :) The MNIST dataset is a handwritten digit recognition problem, in which handwritten digits (represented by 28x28 grayscale images) must be classified into one of the 10 digits (0,1,...,9). I included it here to give a quick introduction to multi-class classification. The dataset is included with Keras and can be accessed using the `dataset_mnist()` function. \n",
    "\n",
    "We first load the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mnist <- dataset_mnist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data preparation\n",
    "\n",
    "As before, we split the data into test and training data, and identify features and outcome variables. The `x_train` and `x_test` data are each 3-d arrays (images,width,height) of grayscale values. The `y_train` and `y_test` data consist of integers from 0-9.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_train <- mnist$train$x\n",
    "y_train <- mnist$train$y\n",
    "x_test <- mnist$test$x\n",
    "y_test <- mnist$test$y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll plot a few of the images. Before we can do that, we need to transform the data into \"long format\", so that we can use **ggplot** for the plotting (thanks to Rob Schlegel for the code below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tranform wide to long, in preparation for ggplot\n",
    "x_train_long <- as.data.frame(x_train) %>%\n",
    "  mutate(picture = 1:nrow(.)) %>%\n",
    "  slice(1:5) %>% # Choose number of pictures to visualise\n",
    "  gather(-picture, key = \"index\", value = \"shade\") %>%\n",
    "  arrange(picture) %>%\n",
    "  group_by(picture) %>%\n",
    "  mutate(index = 1:784,\n",
    "         y = rep(rev(seq(1:28)), 28),\n",
    "         x = rep(1:28, each = 28))\n",
    "\n",
    "# Visualise\n",
    "ggplot(data = x_train_long, aes(x = x, y = y, fill = shade)) +\n",
    "  geom_raster() +\n",
    "  scale_fill_gradient(low = \"white\", high = \"black\") +\n",
    "  coord_equal(expand = 0) +\n",
    "  facet_wrap(~picture, nrow = 1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our model we are going to treat each pixel independently, by \"unravelling\" the 28x28 matrix into a vector of length 784. We also convert the grayscale values from integers ranging between 0 to 255 into floating point values ranging between 0 and 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# reshape\n",
    "dim(x_train) <- c(nrow(x_train), 784)\n",
    "dim(x_test) <- c(nrow(x_test), 784)\n",
    "# rescale\n",
    "x_train <- x_train / 255\n",
    "x_test <- x_test / 255"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output variable is an integer taking on integers 0 through 9, so this is a multi-class classification problem. It is good practice to use **one-hot encoding** for the class values. This transforms the vector of class integers into a binary matrix:\n",
    "\n",
    "```\n",
    "0 -> [1, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
    "1 -> [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]\n",
    "2 -> [0, 0, 1, 0, 0, 0, 0, 0, 0, 0]\n",
    "etc.\n",
    "```\n",
    "One-hot encoding is implemented with the Keras `to_categorical()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_train <- to_categorical(y_train, 10)\n",
    "y_test <- to_categorical(y_test, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now carry out the same model buildings steps as before.\n",
    "\n",
    "#### Create the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model <- keras_model_sequential() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define the model\n",
    "\n",
    "Here we fit a sequential model with two hidden layers and dropout. The first hidden layer consists of 128 neurons, the second 64. Both hidden layers use relu activations. \n",
    "\n",
    "The final output layer now has 10 neurons, one for each of the classes/digits. We use a softmax activation in the output layer. The softmax activation has the property that it forces the outputs of a neural network to sum to 1, so that they can represent valid probabilities for each class. The transformation is given by:\n",
    "\n",
    "$$g(\\textbf{a}) = \\text{softmax}(\\textbf{a}) = \\left[ \\frac{\\exp{(a_1)}}{\\Sigma_c \\exp{(a_c)}} \\ldots \\frac{\\exp{(a_C)}}{\\Sigma_c \\exp{(a_c)}} \\right]^T $$\n",
    "\n",
    "The transformed values - the activations - sum to 1 across classes and are stricly positive. The predicted class is whichever class receives the highest predicted probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% \n",
    "  layer_dense(units = 128, activation = 'relu', input_shape = c(784)) %>% \n",
    "  layer_dropout(rate = 0.4) %>% \n",
    "  layer_dense(units = 64, activation = 'relu') %>%\n",
    "  layer_dropout(rate = 0.3) %>%\n",
    "  layer_dense(units = 10, activation = 'softmax')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compile the model\n",
    "\n",
    "Here our predictions are probability distributions across the ten different digits (e.g. \"we're 80% confident this image is a 3, 10% sure it's an 8, 5% it's a 2, etc.\"), and the target is a probability distribution with 100% for the correct category, and 0 for everything else. We use categorical cross-entropy, an extension of the binary cross-entropy we looked at in the previous example to multi-class classification problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model %>% compile(\n",
    "  loss = 'categorical_crossentropy',\n",
    "  optimizer = 'rmsprop',\n",
    "  metrics = c('accuracy')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Train the model\n",
    "\n",
    "Use the fit() function to train the model for 5 epochs using batches of 128 images. Here we specify a `validation_split`, which further splits the \"training\" data into true \"training\" data and validation data that can be used to do finetuning, in the sense of choosing hyperparameters of the neural network, like the dropout rate. This should *not* be done on the test data, which should always be kept as a totally independent and unseen dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "history <- model %>% fit(\n",
    "  x_train, y_train, \n",
    "  epochs = 5, batch_size = 32, \n",
    "  validation_split = 0.2\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The history object returned by `fit()` includes loss and accuracy metrics which we can plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Evaluate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model %>% evaluate(x_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at some of the images the model incorrectly predicted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get IDs of incorrectly predicted test cases\n",
    "preds <- model %>% predict_classes(x_test)              # get predicted classes for test data\n",
    "wrong_preds <- tibble(id = 1:nrow(y_test),              # make tibble with obs, pred, and id\n",
    "                      obs = mnist$test$y, \n",
    "                      preds = preds) %>%\n",
    "    filter(obs != preds)                                # identify incorrect classifications\n",
    "\n",
    "# create a subset of the test data containing only wrong predictions\n",
    "x_test_wrong <- x_test[wrong_preds$id, ]\n",
    "\n",
    "# tranform wide to long, in preparation for ggplot\n",
    "x_test_wrong_long <- as.data.frame(x_test_wrong) %>%\n",
    "  mutate(picture = 1:nrow(.)) %>%\n",
    "  slice(1:5) %>% # Choose number of pictures to visualise\n",
    "  gather(-picture, key = \"index\", value = \"shade\") %>%\n",
    "  arrange(picture) %>%\n",
    "  group_by(picture) %>%\n",
    "  mutate(index = 1:784,\n",
    "         y = rep(rev(seq(1:28)), 28),\n",
    "         x = rep(1:28, each = 28))\n",
    "\n",
    "# visualise\n",
    "ggplot(data = x_test_wrong_long, aes(x = x, y = y, fill = shade)) +\n",
    "  geom_raster() +\n",
    "  scale_fill_gradient(low = \"white\", high = \"black\") +\n",
    "  coord_equal(expand = 0) +\n",
    "  facet_wrap(~picture, nrow = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try this!\n",
    "\n",
    "Once you're done with predicting the presence of *Aloe dichotoma*, try and predict abundance (the total number of trees, or the log total number)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}