---
title: "Neural Networks - IMDB - 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 IMDB dataset

The objective here is to classify a movie review as either positive or negative.

### Preparing the data

The data has already been preprocessed: the reviews (sequences of words) have been turned into sequences of integers, where each integer stands for a specific word in a dictionary.

* The argument `num_words = 10000` keep only the top 10,000 most frequently occurring words in the training data.

```{r, message=FALSE, warning=FALSE, cache=TRUE}
imdb <- dataset_imdb(num_words = 10000)

train_data <- imdb$train$x
train_labels <- imdb$train$y
test_data <- imdb$test$x
test_labels <- imdb$test$y
```

* Each review is a list of word indices.
* The labels are lists of 0s and 1s, where 0 stands for negative and 1 stands for positive.
* The first review in the list:

```{r, message=FALSE, warning=FALSE}
str(train_data[[1]])

train_labels[[1]]
```

#### Turning sequence of integers back to english

Below is the code to turn the reviews from sequence of integers back to english.

```{r, message=FALSE, warning=FALSE}
word_index <- dataset_imdb_word_index()  
reverse_word_index <- names(word_index)                                    
names(reverse_word_index) <- word_index
decoded_review <- sapply(train_data[[1]], function(index) {                
  word <- if (index >= 3) reverse_word_index[[as.character(index - 3)]]
  if (!is.null(word)) word else "?"
})
```

#### Turning sequence of integers to tensor format

* The `vectorize_sequences` below will produce a tensor of rank 2 of the form `(samples, features)`
* Each sample is represented by a feature vector of the size of the dictionary being used with values equal to 1 if a particular word is present and 0 if the particular word is absent.

```{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)
```

```{r, message=FALSE, warning=FALSE}
y_train <- as.numeric(train_labels)
y_test <- as.numeric(test_labels)
```

### Model definition

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

### Model compilation

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

### Validating your approach

Create a validation set by setting apart 10,000 samples from the original training data.

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

x_val <- x_train[val_indices,]
partial_x_train <- x_train[-val_indices,]
y_val <- y_train[val_indices]
partial_y_train <- y_train[-val_indices]
```

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

Note that the call to `fit()` returns a history object. Let's take a look at it:

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

The history object includes parameters used to fit the model (`history$params`) as well as data for each of the metrics being monitored (`history$metrics`).

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

* You can customize all of this behavior via various arguments to the `plot()` method. 
* We can create custom visualization by using `as.data.frame()` method on the history to obtain a data frame with factors for each metric as well as training versus validation:

```{r, message=FALSE, warning=FALSE}
history_df <- as.data.frame(history)
head(history_df)
```

This fairly naive approach achieves an accuracy of 88%. With state-of-the-art approaches, you should be able to get close to 95%.

### Predicting on new data

```{r, message=FALSE, warning=FALSE}
model %>% predict(x_test[1:10,])
```

### Fighting overfitting

#### Reducing the network's size

Let's try a smaller network:

```{r, message=FALSE, warning=FALSE}
history <- keras_model_sequential() %>%
  layer_dense(units = 4, activation = "relu", input_shape = c(10000)) %>%
  layer_dense(units = 4, activation = "relu") %>%
  layer_dense(units = 1, activation = "sigmoid") %>%
  compile(optimizer = "rmsprop",
    loss = "binary_crossentropy",
    metrics = c("accuracy")) %>%
  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)
```

And a bigger network:

```{r, message=FALSE, warning=FALSE}
history <- keras_model_sequential() %>%
  layer_dense(units = 512, activation = "relu", input_shape = c(10000)) %>%
  layer_dense(units = 512, activation = "relu") %>%
  layer_dense(units = 1, activation = "sigmoid") %>%
  compile(optimizer = "rmsprop",
    loss = "binary_crossentropy",
    metrics = c("accuracy")) %>%
  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)
```

#### Adding weight regularization

Adding L2 weight regularization to the model:

```{r, message=FALSE, warning=FALSE}
history <- keras_model_sequential() %>%
  layer_dense(units = 16, kernel_regularizer = regularizer_l2(0.001),
              activation = "relu", input_shape = c(10000)) %>%
  layer_dense(units = 16, kernel_regularizer = regularizer_l2(0.001),
              activation = "relu") %>%
  layer_dense(units = 1, activation = "sigmoid") %>%
  compile(optimizer = "rmsprop",
    loss = "binary_crossentropy",
    metrics = c("accuracy")) %>%
  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)
```

#### Adding dropout

Let's add two dropout layers in the IMDB network to see how well they do at reducing overfitting.

```{r, message=FALSE, warning=FALSE}
history <- keras_model_sequential() %>%
  layer_dense(units = 16, activation = "relu", input_shape = c(10000)) %>%
  layer_dropout(rate = 0.5) %>%
  layer_dense(units = 16, activation = "relu") %>%
  layer_dropout(rate = 0.5) %>%
  layer_dense(units = 1, activation = "sigmoid") %>%
  compile(optimizer = "rmsprop",
    loss = "binary_crossentropy",
    metrics = c("accuracy")) %>%
  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)
```

