---
title: "Neural Networks - Reuters - Densily connected NN"
subtitle: "TMA4268 Statistical Learning V2018"
author: "Thiago G. Martins, Department of Mathematical Sciences, NTNU"
date: "Spring 2018"
output: #3rd letter intentation hierarchy
#  beamer_presentation:
###    incremental: true # or >* for one at a time
#  slidy_presentation:
#    font_adjustment: +1  
   prettydoc::html_pretty:
    theme: cayman
    highlight: github
    toc: true
    toc_depth: 4
   # prettydoc::html_pretty:
   #  theme: architect
   #  highlight: github
#  pdf_document:
  #  toc: true
  #  toc_depth: 2

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
library(keras)
```

## The Reuters dataset

The objective here is to classify short news stories into one of 46 topics available.

### Preparing the data

Here, we use the multi-assignment operator (`%<-%`) from the `zeallot` package to unpack the list into a set of distinct variables.

```{r, message=FALSE, warning=FALSE}
reuters <- dataset_reuters(num_words = 10000)
c(c(train_data, train_labels), c(test_data, test_labels)) %<-% reuters
```

```{r, message=FALSE, warning=FALSE}
length(train_data)
length(test_data)
```

As with the IMDB reviews, each example is a list of integers (word indices):

```{r, message=FALSE, warning=FALSE}
train_data[[1]]
train_labels[[1]]
```

You can vectorize the data with the exact same code as in the IMDB example

```{r, message=FALSE, warning=FALSE}
vectorize_sequences <- function(sequences, dimension = 10000) {
  results <- matrix(0, nrow = length(sequences), ncol = dimension)
  for (i in 1:length(sequences))
    results[i, sequences[[i]]] <- 1
  results
}

x_train <- vectorize_sequences(train_data)          
x_test <- vectorize_sequences(test_data)            
```

Vectorize the labels:

```{r, message=FALSE, warning=FALSE}
one_hot_train_labels <- to_categorical(train_labels)
one_hot_test_labels <- to_categorical(test_labels)
```

### Building the model

* The dimensionality of the output space (46 classes) is much larger.

Information bottleneck

* Each layer can only access information present in the output of the previous layer. 
* Each layer can potentially become an information bottleneck. 
* A 16-dimensional intermediate layer may be too limited to learn to separate 46 different classes: 
* Such small layers may act as information bottlenecks, permanently dropping relevant information.

For this reason we will use larger layers. Let's go with 64 units.

```{r, message=FALSE, warning=FALSE}
model <- keras_model_sequential() %>%
  layer_dense(units = 64, activation = "relu", input_shape = c(10000)) %>%
  layer_dense(units = 64, activation = "relu") %>%
  layer_dense(units = 46, activation = "softmax")
```
  
### Compiling the model

The best loss function to use in this case is `categorical_crossentropy`.

```{r, message=FALSE, warning=FALSE}
model %>% compile(
  optimizer = "rmsprop",
  loss = "categorical_crossentropy",
  metrics = c("accuracy")
)
```

### Validating your approach

Let’s set apart 1000 samples in the training data to use as a validation set.

```{r, message=FALSE, warning=FALSE}
val_indices <- 1:1000

x_val <- x_train[val_indices,]
partial_x_train <- x_train[-val_indices,]

y_val <- one_hot_train_labels[val_indices,]
partial_y_train = one_hot_train_labels[-val_indices,]
```

Now, let’s train the network for 20 epochs.

```{r, message=FALSE, warning=FALSE}
history <- model %>% fit(
  partial_x_train,
  partial_y_train,
  epochs = 20,
  batch_size = 512,
  validation_data = list(x_val, y_val)
)
```

```{r, message=FALSE, warning=FALSE}
plot(history)
```

The network begins to overfit after nine epochs. Let’s train a new network from scratch for nine epochs and then evaluate it on the test set.

```{r, message=FALSE, warning=FALSE}
model <- keras_model_sequential() %>%
  layer_dense(units = 64, activation = "relu", input_shape = c(10000)) %>%
  layer_dense(units = 64, activation = "relu") %>%
  layer_dense(units = 46, activation = "softmax")

model %>% compile(
  optimizer = "rmsprop",
  loss = "categorical_crossentropy",
  metrics = c("accuracy")
)

history <- model %>% fit(
  partial_x_train,
  partial_y_train,
  epochs = 9,
  batch_size = 512,
  validation_data = list(x_val, y_val)
)
```

```{r, message=FALSE, warning=FALSE}
results <- model %>% evaluate(x_test, one_hot_test_labels)

results
```

This approach reaches an accuracy of ~ 79%. With a balanced binary classification problem, the accuracy reached by a purely random classifier would be 50%. But in this case it’s closer to 18%, so the results seem pretty good, at least when compared to a random baseline:

```{r, message=FALSE, warning=FALSE}
test_labels_copy <- test_labels
test_labels_copy <- sample(test_labels_copy)
length(which(test_labels == test_labels_copy)) / length(test_labels)
```

### Predictions on new data

```{r, message=FALSE, warning=FALSE}
predictions <- model %>% predict(x_test)
```

Each entry in predictions is a vector of length 46:

```{r, message=FALSE, warning=FALSE}
dim(predictions)
```

The coefficients in this vector sum to 1:

```{r, message=FALSE, warning=FALSE}
sum(predictions[1,])
```

The largest entry is the predicted class—the class with the highest probability:

```{r, message=FALSE, warning=FALSE}
which.max(predictions[1,])
```

### Recommended exercise 2

Use a vector of integers as labels instead of the one-hot encoding used above. Remember that this choice will impact the loss function used to train the model.


