03-transfer-xray.Rmd

---
title: "Deep Learning in R: Part 3 (Fine-tuning)"
output: html_document
---

## Install dependencies

One this chunk once manually if needed. It won't be executed when knitting the file or selecting "run all chunks".

```{r install_deps, eval = FALSE}
# Install magick if not already installed.
install.packages("magick")
```

ImageMagick is an open source image processing library in C++ that is supported in many different languages; we are using an R package that interfaces with the C++ code. See details in the [magick vignette](https://cran.r-project.org/web/packages/magick/vignettes/intro.html).

## Load packages

Typically when conducting statistical analysis we want to ensure reproducible results. In other words, if I run the analysis file twice I should get exactly the same estimated accuracy each time. This is important when publishing research and when doing experiments.

Unfortunately deep learning in Keras and Tensorflow is not easily reproducible. RStudio recommends that we run the `use_session_with_seed()` function immediately after loading the keras package. We may need to also disable both GPU support and use of multiple CPU cores to guarantee reproducible results. Realistically that would be too slow in many cases, so we may have to settle with not necessarily having 100% reproducible results. See more details in the [RStudio Keras FAQ](https://keras.rstudio.com/articles/faq.html#how-can-i-obtain-reproducible-results-using-keras-during-development).

```{r load_packages}

library(keras)

# Set the seed immediately after loading keras library.
# If the library was already loaded, we may want to run .rs.restartR() to restart
# our R session.
# This is broken right now.
#use_session_with_seed(1, disable_gpu = FALSE, disable_parallel_cpu = FALSE)

library(dplyr)
library(ggplot2)
# Gives us the tf object for certain manipulations.
library(tensorflow)

```

## Download imaging data

In parts 3 and 4 we don't store the image data in the GitHub repository directly, because
Git is intended to track changes text files and doesn't work as well for binary files like images ((see one opinionated article on this)[https://robinwinslow.uk/2013/06/11/dont-ever-commit-binary-files-to-git/]). Instead we will download a zip file and unzip it programmatically. This keeps the repository size smaller and quicker to clone or download.

```{r download_data}
local_file = "data-raw/medical_images.zip"

# Download the zip file if we don't already have it in our working directory.
# It is 13.5 MB in size.
if (!file.exists(local_file)) {
  download.file("https://github.com/paras42/Hello_World_Deep_Learning/raw/9921a12c905c00a88898121d5dc538e3b524e520/Open_I_abd_vs_CXRs.zip", local_file)
}

# Unzip the local file into the data-raw directory, if we haven't already.
if (!file.exists("data-raw/Open_I_abd_vs_CXRs")) {
  unzip(local_file, exdir = "data-raw")
}  
```

`data-raw` is a common folder name to store "raw" data, i.e. prior to conducting any cleaning or processing. This could include .csv files or other exports generated by another system (SAS files, survey results, etc.). Files in this directory should be treated as "read-only" to ensure that any analysis can always be re-run. Then if we want to save processed files we would keep them in another directory, like "data" or "data-processed".

## Setup paths

We will use the `dirs` list to organize our different directories into a few variables:

  * `dirs$base` - path to our main directory with the images unzipped.
  * `dirs$train` - path to the training set folder.
  * `dirs$val` - path to the validation set folder.
  
It's a good practice to organize multiple settings into a single list object. For one, it keeps the environment more organized, and it also makes it easier to save those settings into an RData file for posterity.

```{r setup_paths}
# Organize our directories into a list.
dirs = list(base = "data-raw/Open_I_abd_vs_CXRs/")

# Don't include a "/" at the end of these, because list.files() will add later.
dirs$train = paste0(dirs$base, "TRAIN")
dirs$val = paste0(dirs$base, "VAL")

# Print out our directory configuration.
dirs
```

## Examine paths

How many images are in each directory?

```{r examine_paths}

# Images are organized into two subdirectories (1 per type of image), so
# we set recursive = TRUE to go into each of the subdirectories.
length((train_files = list.files(dirs$train, recursive = TRUE, full.names = TRUE)))

# Review the first two file elements. Both are abdominal xrays (note the subdirectory).
train_files[1:2]

# Challenge:
# How many files are in our validation directory?
# What are the 2nd and 3rd elements in our validation file list?

# Hint:
# length((val_files = list.files(_________, ________, ________)))
```

## Plot images

```{r plot_images}

# The magick package is one way to display images.
library(magick)

# Display an image in R
print(image_read(train_files[5]))
# Open a system review to see the image.
image_browse(image_read(train_files[5]))

# Alternatively we can use cowplot to print images inside of a ggplot.
library(cowplot)

# Plot it with ggplot + cowplot + magick packages.
ggdraw() + draw_image(train_files[1])

# Plot the second image also, this time adding a title and removing extra items.
ggdraw() + draw_image(train_files[2]) + ggtitle("2nd image") + theme_minimal() +
  theme(axis.text = element_blank(), panel.grid = element_blank())

###
# Challenge: plot the third image and put the filename in the title.
# Bonus hint: basename() will remove all directories from a file path.
###

```

## Pre-trained network + Fine tuning

In order to reach high levels of accuracy, neural networks generally need to be very deep, which means many hidden layers. This allows the network to build a complex understanding of the input data because deeper layers build upon earlier layers to "engineer" new features that it has learned are related to predicting accurately.

However deep neural networks consist of millions of links (weights) between neurons, and if we have only 100 observations the network doesn't have sufficient information to adjust the weights accurately. As a result deep neural networks built from scratch require a large number of observations: preferably millions, but possibly down to 50k (depending on the task).

What has emerged as a way around this limitation is "pre-training + fine-tuning". First, we train a complex (deep) neural network on a large dataset, which allows the layers to calibrate their weights well. Our hope is that the features it engineers on that dataset will be applicable to our new, smaller dataset. This is known as the "base model", and Keras supplies 10 such pre-trained models for our usage.

Next we apply that pre-trained neural network to our smaller dataset. We remove the last layer or two which are more specific to the original dataset, and add on a new layer or two that are randomly initialized. **Fine-tuning** is then re-running the network on our smaller dataset and updating primarily those new layers to our small dataset; we may also update the "base" neural network weights slightly.

We are using Inception v3 to jump-start our deep learning, but [many others are available in Keras](https://keras.rstudio.com/articles/applications.html).

## Setup data and core model

```{r setup_model}
# Dimensions of our images as expected by the neural architecture.
img_width = img_height = 299L
batch_size = 5L

train_datagen = keras::image_data_generator(rescale = 1/255)

val_datagen = keras::image_data_generator(rescale = 1/255)

train_gen =
  train_datagen$flow_from_directory(dirs$train,
                                    # Images will be resized to this target size.
                                    target_size = c(img_width, img_height),
                                    batch_size = batch_size,
                                    class_mode = "binary")

val_gen =
  val_datagen$flow_from_directory(dirs$val,
                                  target_size = c(img_width, img_height),
                                  batch_size = batch_size,
                                  class_mode = "binary")

# This will download the inception weights the first time it is run (~84 MB)
base_model = keras::application_inception_v3(include_top = FALSE,
                                             pooling = "avg",
                                             input_shape = c(img_width, img_height, 3L))
# Outputs an 8x8x2048 tensor.
base_model$output_shape
?application_inception_v3

summary(base_model)
```

## Add custom layer for our task

```{r}

# first: train only the top layers (which were randomly initialized)
# i.e. freeze all convolutional InceptionV3 layers
# This is not working - possible bug in RStudio's Keras package.
freeze_weights(base_model)

# Add custom layer to inception.
model_top = base_model$output %>%
  #layer_global_average_pooling_2d() %>%
  layer_dense(units = 128, activation = "relu") %>%
  layer_dropout(0.5) %>%
  layer_dense(units = 1, activation = "sigmoid")

# this is the model we will train
model = keras_model(inputs = base_model$input, outputs = model_top)

length(model$layers)

# Manually freeze the original inception layers, just train the last 3 layers.
freeze_weights(model, 1, length(model$layers) - 3)
summary(model)

# Compile the model (should be done *after* setting layers to non-trainable)
model %>%
  compile(optimizer =
            optimizer_adam(#lr = 0.00001,
                           #lr = 0.0005,
                           lr = 0.0001,
                           # epsilon is a configuratoin setting for the adam optimizer.
                           epsilon = 1e-08),
          # Or can use string version: loss = "binary_crossentropy".
          loss = loss_binary_crossentropy,
          metrics = "accuracy")

(num_train_samples = length(train_files))
num_validation_samples = 10L
```

## Fit model

```{r fit_model}

# Train the model on the new data for a few epochs
history = model %>%
  fit_generator(train_gen,
                steps_per_epoch = as.integer(num_train_samples / batch_size),
                epochs = 5,
                validation_data = val_gen,
                validation_steps = as.integer(num_validation_samples / batch_size))

# Review fitting history.
plot(history)
```

## Train full model

Now we unfreeze the entire model and train for just a few epochs to customize the entire model to our specific dataset.

```{r}
# Unfreeze_weights() seems to require that we explicitly specify the layers.
unfreeze_weights(model, 1, length(model$layers))

model %>%
  compile(optimizer =
            # Notice that our learning rate is lower here.
            optimizer_adam(lr = 0.00001,
                           epsilon = 1e-08),
          loss = loss_binary_crossentropy,
          metrics = "accuracy")

# Train the full set of layers but only for a few epochs. 
history = model %>%
  fit_generator(train_gen,
                steps_per_epoch = as.integer(num_train_samples / batch_size),
                epochs = 4,
                validation_data = val_gen,
                validation_steps = as.integer(num_validation_samples / batch_size))
```

Our validation loss is now lower than the training loss - any ideas why that would be?

Challenges:
1. Go back to the first model and try modifying some settings: number of hidden units, additional hidden layers, and/or dropout rate. How do they each affect your results?
2. Review the help page for optimizer_adam and try changing a setting or two. Click to the other optimizers and try a different optimizer.

## Data augmentation

As optional homework, here is sample code for a more complex image data generator. This code uses **data augmentation**, which applies random small perturbations to the original images to 1) approximate having a larger sample size, and 2) encourage the network to be less sensitive to the original image characteristics. As a result, this augmentation step should allow us to achieve better performance compared to not doing it.

```{r data_augmentation}
train_datagen =
  keras::image_data_generator(rescale = 1 / 255,
                              shear_range = 0.2,
                              zoom_range = 0.2,
                              rotation_range = 20,
                              width_shift_range = 0.2,
                              height_shift_range = 0.2,
                              horizontal_flip = TRUE)

```

## Acknowledgements

Material drawn from (https://github.com/paras42/Hello_World_Deep_Learning) and (https://keras.rstudio.com).