---
title: 'Diabolo Network and Deep Features'
author: 'Loic Merckel'
date: '10 Mars 2018'
output:
  html_document:
    number_sections: false
    toc: true
    highlight: tango
    theme: cosmo
    smart: true
---

<style type="text/css">
h1.title { font-weight: bold; } h1 { font-weight: normal; } .author { font-weight: normal; font-size: 1.5em; }
</style>


```{r include=FALSE}
# License -----------------

# Copyright 2017 Loic Merckel
# 
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# 
# http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
```

```{r include=FALSE}
pkgs <- c("corrplot", "caret", "data.table", "plyr",
          "xgboost", "parallel", "Metrics", "maps", 
          "ggmap", "lubridate", "fasttime", "gridExtra", 
          "geosphere", "vtreat", "parallel", "doParallel") 
for (pkg in pkgs) {
  if (! (pkg %in% rownames(installed.packages()))) { install.packages(pkg) }
  require(pkg, character.only = TRUE)
}
rm(pkgs)
rm(pkg)
```

```{r training_data_chunk, include=FALSE}
# data ---------------------------
rm(list=ls(all=TRUE))

kIsOnKaggle <- TRUE
kSaveResults <- FALSE # the kernel would timeout

X <- fread(file.path(ifelse(kIsOnKaggle, "..", "."), "input", "train.csv", fsep = .Platform$file.sep), 
           header = TRUE, data.table = TRUE, na.strings=c("NA","?", ""))

meanIsAttFullDataset <- mean(X$is_attributed)
fullDatasetSize <- nrow(X)

kMaxRows <- ifelse(kIsOnKaggle, 1000000, 200000)
if(nrow(X) > kMaxRows) {
  set.seed(1)
  X <- X[sample(1:nrow(X), kMaxRows), ]
}
rm(kMaxRows)
```

Observing how vastly imbalanced are the classes in this binary classification problem, anomaly detection techniques might prove fruitful. Throughout the kernel we explore how an autoencoder could help in distinguishing the "anomalies" (i.e., the legit clicks) from the rest.

Given our severe lack of computational resources at the moment (basically, what Kaggle offers with its [kernels](https://www.kaggle.com/wiki/Scripts)), we only consider a tiny fraction of the entire dataset (`r round(nrow(X) / fullDatasetSize * 100, 4)`% of the full dataset), especially the validation set for model selection is most likely too small, leading to overly optimistic model performance.


```{r include=FALSE}
X[, ip := as.factor(ip)]
X[, app := as.factor(app)]
X[, device := as.factor(device)]
X[, os := as.factor(os)]
X[, channel := as.factor(channel)]
X[, is_attributed := as.factor(is_attributed)]

X[, attributed_time := NULL]

responseFactorLevels <- levels(X$is_attributed)
```

```{r include=FALSE}
set.seed(1)
trainIndexes <- createDataPartition(X$is_attributed, p = 0.8, list = FALSE)
# DO NOT CHANGE ORDER
Xv <- X[-trainIndexes, ]
X <- X[trainIndexes, ]
rm(trainIndexes)
```

# Features Tinkering

```{r include=FALSE}
toIgnore <- c("id", "is_attributed")
```

## Date

```{r data_train_chunk, echo=TRUE, warning=FALSE, results='hide', message=FALSE}
X[, click_time := fastPOSIXct (click_time)]
X[, wday := lubridate::wday(click_time)]
X[, hour := lubridate::hour (click_time)]
```

```{r data_valid_chunk, include=FALSE}
Xv[, click_time := fastPOSIXct (click_time)]
Xv[, wday := lubridate::wday(click_time)]
Xv[, hour := lubridate::hour (click_time)]

X[, click_time := NULL]
Xv[, click_time := NULL]
```

## Convert Categorical Variables to Numerical

We rely on the [vtreat](https://cran.r-project.org/web/packages/vtreat/index.html) package to encode the high cardinality categorical variables via the so-called impact encoding (also referred to as likelihood encoding).

```{r feature_processing_train_chunk, echo=TRUE, warning=FALSE, results='hide', message=FALSE}
set.seed(1)
if(!kIsOnKaggle){
  cl <- makeCluster(parallel::detectCores())
  registerDoParallel(cl)
  getDoParWorkers()
} else { cl <- NULL }

set.seed(1)
cfe <- mkCrossFrameCExperiment(
  dframe = X,
  varlist = setdiff(names(X), toIgnore),
  outcomename = "is_attributed",
  outcometarget = 1,
  ncross = 8,
  parallelCluster = cl)
if(!is.null(cl)) { stopCluster(cl) }
rm(cl)

X <- as.data.table(cfe$crossFrame)
X[, click_id := as.numeric(NA)]
```

```{r include=FALSE}
toIgnore <- c("click_id")
featuresTreatment <- cfe$treatments
rm(cfe)
```

```{r feature_processing_valid_chunk, echo=TRUE, warning=FALSE, results='hide', message=FALSE}
# validation set
Xv <- as.data.table(prepare(featuresTreatment, Xv))
Xv[, click_id := as.numeric(NA)]
```


# Deep Features

## Convert to H2O Frames

```{r h2o_init_chunk, include=FALSE}
if (! ("h2o" %in% rownames(installed.packages()))) { install.packages("h2o") }
require("h2o")

tryCatch(
  if (h2o.clusterIsUp()) {
    h2o.shutdown(prompt=FALSE)
    Sys.sleep(5)
  }, error = function(e) {
    
  })

h2o.init(nthreads = parallel:::detectCores(), 
         max_mem_size = "15g", min_mem_size = "1g")
h2o.removeAll()
h2o.clusterStatus()
```

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
train <- as.h2o(X)
predictors <- setdiff (names(X), c(toIgnore))
```

```{r include=FALSE}
valid <- as.h2o(Xv)
rm(X, Xv)
```


## Autoencoder (Diabolo Network)

```{r autoencoder_chunk, echo=TRUE, warning=FALSE, results='hide', message=FALSE}
hyperParamsAutoencoder = list( 
  hidden = list(c(40, 25, 40), c(40, 20, 40)), #c(30, 18, 30)  
  activation = c("Tanh"))

gridAutoencoder <- h2o.grid(
  x = predictors,
  autoencoder = TRUE,
  training_frame = train,
  hyper_params = hyperParamsAutoencoder,
  search_criteria = list(strategy = "Cartesian"),
  algorithm = "deeplearning",
  grid_id = "grid_autoencoder", 
  reproducible = TRUE, 
  seed = 1,
  variable_importances = TRUE,
  categorical_encoding = "AUTO",
  score_interval = 10,
  epochs = 800,
  adaptive_rate = TRUE,
  standardize = TRUE,
  ignore_const_cols = FALSE)
```

```{r include=FALSE}
rm (hyperParamsAutoencoder)
```

The following table summarizes the grid results (it is sorted increasingly by 'mse'):
```{r echo=FALSE}
sortedGridAutoencoder <- h2o.getGrid(
  "grid_autoencoder", 
  sort_by = "mse", decreasing = FALSE)
tmpDf <- as.data.frame(sortedGridAutoencoder@summary_table)
knitr::kable(head(tmpDf[, -grep("model_ids", colnames(tmpDf))]), row.names = TRUE)
rm(tmpDf)
```

```{r include=FALSE}
bestAutoencoder <- h2o.getModel(sortedGridAutoencoder@model_ids[[1]])

# to free some memory
i <- 2
while(TRUE) {
  if(i > length(sortedGridAutoencoder@model_ids)) { break }
  h2o.rm(h2o.getModel(sortedGridAutoencoder@model_ids[[i]]))
  i <- i+1
}
rm(i)

bestAutoencoderErr <- as.data.frame(
  h2o.anomaly(
    bestAutoencoder, 
    train, 
    per_feature = FALSE))
```

Considering the "best" autoencoder (i.e., the one with the lowest 'mse', which is the one with the hidden layers [`r bestAutoencoder@parameters$hidden`]), the two following figures illustrate the fact that it performs rather well; only a limited portion of the input signal could not be reconstructed. 

```{r  echo=FALSE, warning=FALSE, results='hide', message=FALSE, fig.width=10, fig.height=3.5}
plotReconstructionError <- function (error) {
  cut <- 0.5 * sd(error)
  sortedErr <- sort(error)
  sortedErrFrame <- data.frame (index = seq(0, length(sortedErr)-1), error = sortedErr)
  ylim <- c(min(sortedErrFrame$error), max(sortedErrFrame$error))
  xlim <- c(min(sortedErrFrame$index), max(sortedErrFrame$index))
  # could do as in https://stackoverflow.com/questions/11838278/plot-with-conditional-colors-based-on-values-in-r
  plot (x = sortedErrFrame$index[which(sortedErrFrame$error <= cut)], 
        y = sortedErrFrame$error[which(sortedErrFrame$error <= cut)], 
        type = "o", col="forestgreen", lwd=1, ylim = ylim, xlim = xlim, ylab="mse",
        main = "Reconstruction Error",
        xlab = "Sorted Index")
  par(new=TRUE)
  plot (x = sortedErrFrame$index[which(sortedErrFrame$error > cut)], 
        y = sortedErrFrame$error[which(sortedErrFrame$error > cut)],  
        type = "o", col="firebrick3", xaxt='n', yaxt='n', ann=FALSE, ylim = ylim, xlim = xlim, xlab="")
  return (0)
}

layout(matrix(c(1,2,2), 1, 3, byrow = TRUE))
plotReconstructionError (bestAutoencoderErr$Reconstruction.MSE)

# https://stackoverflow.com/questions/21858394/partially-color-histogram-in-r
h <- hist(x = bestAutoencoderErr$Reconstruction.MSE, breaks = 100, plot = FALSE)
cuts <- cut(h$breaks, c(-Inf, 0.5 * sd(bestAutoencoderErr$Reconstruction.MSE), Inf))
plot(h, col = c("forestgreen","firebrick3")[cuts], main = "Reconstruction Error", xlab = "mse", lty="blank")
rm(h, cuts)
```

```{r include=FALSE}
rm(bestAutoencoderErr)
```


## Deep Features Visualization

```{r  include=FALSE}
plotDeepFeatures <- function(data, maxPlot = 16, ncol = 4) {
  count <- 1
  plotList <- list()
  n <- (ncol(data) - 1)
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      plotList[[paste0("p", count)]] <- ggplot(
        data, 
        aes_string(
          x = paste0("DF.L", layer, ".C", i), 
          y = paste0("DF.L", layer, ".C", j), 
          color = "is_attributed")) +
        geom_point(alpha = 0.5, aes(colour = is_attributed)) +
        theme(legend.position = 
                ifelse(count == min((n-1)*n, maxPlot), "right", "none")) +     
        labs(color="is_at")
      
      count <- count + 1
      if (count > maxPlot) {
        break
      }
    }
    if (count > maxPlot) {
      break
    }
  }
  grid.arrange(grobs = as.list(plotList), ncol = ncol)
}
```


### Second Layer

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
layer <- 2
```

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
deepFeature2 <- h2o.deepfeatures(bestAutoencoder, train, layer = layer)
```

```{r include=FALSE}
data <- as.data.frame(deepFeature2)
data$is_attributed <- as.factor(as.vector(train$is_attributed))

summary(data)
```

```{r  second_layer_plot_chunk, echo=FALSE, warning=FALSE, results='hide', message=FALSE, fig.width=10, fig.height=10}
plotDeepFeatures(data, 16)
```

```{r include=FALSE}
rm (deepFeature2, data)
```


### Third Layer

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
layer <- 3
```

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
deepFeature3 <- h2o.deepfeatures(bestAutoencoder, train, layer = layer)
```

```{r include=FALSE}
data <- as.data.frame(deepFeature3)
data$is_attributed <- as.factor(as.vector(train$is_attributed))

summary(data)
```

```{r third_layer_plot_chunk, echo=FALSE, warning=FALSE, results='hide', message=FALSE, fig.width=10, fig.height=10}
plotDeepFeatures(data, 16)
```

```{r include=FALSE}
rm (deepFeature3, data)
```


# Predictions Using Deep Features and GBM

We use the second layer of the autocencoder (`bestAutoencoder`) and the gradient boosting machine algorithm offered by H<small>2</small>O (`h2o.gbm`). Those are arbitrary choices.

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
layer <- 2
```


## Get Deep Features

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
deepFeatureTrain <- h2o.deepfeatures(bestAutoencoder, train, layer = layer)
deepFeatureTrain[["is_attributed"]] <- as.h2o(as.factor(as.vector(train$is_attributed)))
```

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
deepFeatureValid <- h2o.deepfeatures(bestAutoencoder, valid, layer = layer)
deepFeatureValid[["is_attributed"]] <- as.h2o(as.factor(as.vector(valid$is_attributed)))
```


## Grid Search

```{r gbm_chunk, echo=TRUE, warning=FALSE, results='hide', message=FALSE}
deepfeatures <- setdiff(names(deepFeatureTrain), c("is_attributed"))

hyperParamsGbm <- list(
  max_depth = c(3, 4),
  histogram_type = c("UniformAdaptive", "QuantilesGlobal", "RoundRobin"),
  min_split_improvement = c(0, 1e-6))

gridGbm <- h2o.grid(
  hyper_params = hyperParamsGbm,
  search_criteria = list(strategy = "Cartesian"),
  algorithm = "gbm",
  grid_id = "grid_gbm", 
  x = deepfeatures, 
  y = "is_attributed", 
  training_frame = deepFeatureTrain, 
  validation_frame = deepFeatureValid,
  nfolds = 0,
  ntrees = 400,                                        
  learn_rate = 0.05,                                                         
  learn_rate_annealing = 0.99,                                               
  max_runtime_secs = 1200,                              
  stopping_rounds = 5, 
  stopping_tolerance = 1e-5, 
  stopping_metric = "logloss", 
  score_tree_interval = 10,                                                
  seed = 1)

sortedGridGbm <- h2o.getGrid("grid_gbm", sort_by = "logloss", decreasing = FALSE)
```

```{r echo=FALSE}
tmpDf <- as.data.frame(sortedGridGbm@summary_table)
knitr::kable(head(tmpDf[, -grep("model_ids", colnames(tmpDf))]), row.names = TRUE)
rm(tmpDf)
```


## Best Model & Performance on the Validation Set

```{r echo=TRUE, warning=FALSE, results='hide', message=FALSE}
bestGbmModel <- h2o.getModel(sortedGridGbm@model_ids[[1]])
perf <- h2o.performance(bestGbmModel, valid = TRUE)
```

```{r echo=FALSE, warning=FALSE, results='show', message=FALSE}
perf
```

```{r include=FALSE}
h2o.rm(
  c(h2o.getId(train), h2o.getId(valid), 
    h2o.getId(deepFeatureTrain), h2o.getId(deepFeatureValid)))
rm(train, valid, deepFeatureTrain, deepFeatureValid)
```


<!--
## Test Set & Submission
-->

```{r test_data_chunk, include=FALSE}
if(kSaveResults){
  Xt <- fread(file.path(ifelse(kIsOnKaggle, "..", "."), 
                        "input", 
                        ifelse(kIsOnKaggle, "test.csv", "test.csv"), 
                        fsep = .Platform$file.sep), 
              header = TRUE, data.table = TRUE, na.strings=c("NA","?", ""))

  Xt[, is_attributed := NA]
  
  Xt[, ip := as.factor(ip)]
  Xt[, app := as.factor(app)]
  Xt[, device := as.factor(device)]
  Xt[, os := as.factor(os)]
  Xt[, channel := as.factor(channel)]
  
  Xt[, click_time := fastPOSIXct (click_time)]
  Xt[, wday := lubridate::wday(click_time)]
  Xt[, hour := lubridate::hour (click_time)]
  
  Xt[, click_time := NULL]
}
```

```{r predict_test_and_save_chunk, include=FALSE}
if(kSaveResults){
  step <- 500000
  predTest <- c()
  for (k in 1:((nrow(Xt) + step) %/% step)) {
    testk <- Xt[((k-1) * step + 1):(min(nrow(Xt), k * step)), ]
    testk <- as.data.table(
      cbind(click_id = testk$click_id, prepare(featuresTreatment, testk)))
    
    testkH2o <- as.h2o(testk)
    
    rm(testk)
    
    deepFeatureTest <- h2o.deepfeatures(bestAutoencoder, testkH2o, layer = layer)
    
    pred <- h2o.predict (bestGbmModel, newdata = deepFeatureTest)
    predTest <- c(predTest, as.vector(pred$predict))
    h2o.rm(c(h2o.getId(testkH2o), h2o.getId(pred), h2o.getId(deepFeatureTest)))
    
    # the case where nrow(Xt) is a multiple of step
    if (nrow(Xt) <= k * step) { break }
  }
  rm(k, step, pred)
  #rm(featuresTreatment)
  #rm(Xt)
  
  fwrite(data.table(click_id = Xt$click_id, is_attributed = predTest), 
         file = file.path(".", paste0("output-gbm-", Sys.Date(), ".csv"), 
                          fsep = .Platform$file.sep), 
         row.names = FALSE,  quote = FALSE)
}
```

```{r echo=FALSE, warning=FALSE, results='show', message=FALSE}
if(kSaveResults){
  # those two quantities should be roughly equal... Otherwise it is rather suspecious...
  mean(predTest == "1")
  meanIsAttFullDataset
}
```