{"cells":[{"metadata":{},"cell_type":"markdown","source":"# R for LANL 2019: feature generation\n\n## Load packages"},{"metadata":{"trusted":true},"cell_type":"code","source":"library(keras)\nlibrary(data.table)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Define function for Data-aggregation / feature generation\n\n Aggregate the input data in groups data-rows with the Test-Batchsize (150000 rows):\n - Distribution-parameters based on raw input and absolute values (averaged).\n - Percentage of values higher than x * standard deviation (x= $2^{-1}\\dots2^8$ ).\n - Peak-Counts, width a.s.o. from the pracma-package - applied on the sum of 10-100 absolute values of the input data\n - Autocorrelation of the raw input."},{"metadata":{"trusted":true},"cell_type":"code","source":"EQ_Aggregate <- function(x, sdInput, pInput, mInput) {\n  x <- x - mInput\n  x10abs <- colMeans(matrix(abs(x), nrow = 10))\n  x100abs <- colMeans(matrix(abs(x), nrow = 100))\n  PeakInfo <- function(x, nups = 1, minpeakheight = sdInput/2) {\n    y <- pracma::findpeaks(x, nups = nups, minpeakheight = minpeakheight)\n    if (is.null(y)) {\n      PeakZus <- rep(0, 11)\n    } else { # Output of PeakInfo: number of peaks and some quantiles of the peak height and width\n      PeakZus <- c(nrow(y), \n                   quantile(y[, 1], probs = c(0.05, 0.25, 0.5, 0.75, 0.95), names = FALSE), \n                   quantile(y[, 4] - y[, 3], probs = c(0.05, 0.25, 0.5, 0.75, 0.95), names = FALSE) \n      )\n    }\n    return(PeakZus)\n  }\n  c(sd(x)/sdInput,\n    mean(x100abs),\n    quantile(x, probs = c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95), names = FALSE),\n    sapply(2^(-2:8), function(g, sd = sdInput, h = x) sum(abs(h) > g*sd)),\n    PeakInfo(x10abs, 3, 2*sdInput),\n    PeakInfo(x100abs, 1, sdInput/8)[1:6],\n    PeakInfo(x100abs, 1, 2*sdInput)[1:6],\n    as.vector(acf(x, lag.max = 50,plot = FALSE)$acf)[2:51]\n    )\n}\ntestBatchsize <- 150000","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load input data"},{"metadata":{"trusted":true},"cell_type":"code","source":"system.time(trainEQ <- fread(\"../input/LANL-Earthquake-Prediction/train.csv\", header = TRUE))\ncat(paste(\"number of rows:\", nrow(trainEQ)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First calculate the input for the aggregation function above(mean, peak, and standard deviation) then output the number of parameters.\nFinally chose a row to start with and determine the last row to be used for training and the number of batches"},{"metadata":{"trusted":true},"cell_type":"code","source":"mInput <- mean(trainEQ$acoustic_data)\npInput <- max(abs(trainEQ$acoustic_data))\nsdInput <- sd(trainEQ$acoustic_data)\nnParameter <- length(EQ_Aggregate(trainEQ$acoustic_data[1:testBatchsize], sdInput, pInput, mInput))\ncat(paste(\"number of aggregated parameters:\", nParameter))\nset.seed(1)\nfirstRow <- 1 + round(runif(1, 0, 100))\nlastRow <- (((nrow(trainEQ) - firstRow) %/% testBatchsize ) - 1) * testBatchsize + firstRow\nnBatches <- ((lastRow - firstRow) / testBatchsize) + 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature generation\n\nWrite a matrix as an input for neural network model - based on applying the above function in a loop on batches of the same size as the Test-Data. Check out how long it took by system-time. Tried parallel processing (e.g. future, future.apply, foreach) as well: no improvement - mainly due to RAM-limit of 16GB."},{"metadata":{"trusted":true},"cell_type":"code","source":"system.time({\n    trainMatrix <- matrix(NA, nrow = nBatches, ncol = nParameter)\n    l <- 0\n    for (j in seq(firstRow, lastRow, by = testBatchsize)) {\n        k <- j + testBatchsize - 1\n        l <- l+1\n        trainMatrix[l, ] <- EQ_Aggregate(trainEQ$acoustic_data[j:k], sdInput, pInput, mInput)\n    }\n    trainMatrix <- scale(trainMatrix) \n    colMeansTrain <- attr(trainMatrix, \"scaled:center\")\n    colStdAbwTrain <- attr(trainMatrix, \"scaled:scale\")\n    trainMatrix <- trainMatrix[, colStdAbwTrain != 0]\n})\ncat(paste(\"dimensions of the final training matrix:\", dim(trainMatrix)[1], \"rows, \", dim(trainMatrix)[2], \"columns.\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Generate Target-Values and remove earthquake-containing rows\n\nThe target values are taken from the center value of the batches of time_to_failure-row with test-batch length - equals the median as the data is sorted decreasingly. Batches in which an earthquake actually occurred are removed (collected in the vector removeRow first). The criterion for earthquake within the batch is simple: time_to_failure in the beginning is lower then in the end of the batch."},{"metadata":{"trusted":true},"cell_type":"code","source":"trainTarget <- vector(mode = \"double\", length = nBatches)\nremoveRow <- NULL\nl <- 0\nfor (j in seq(firstRow, lastRow, by = testBatchsize)) {\n    k <- j + testBatchsize - 1\n    l <- l+1\n    trainTarget[l] <- trainEQ$time_to_failure[j + testBatchsize/2]\n    if (trainEQ$time_to_failure[j] < trainEQ$time_to_failure[k]) removeRow <- c(removeRow, l)\n}\ntrainMatrix <- trainMatrix[-removeRow, ]\ntrainTarget <- trainTarget[-removeRow]\ncat(paste(\"number of removed rows:\", length(removeRow)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Build and train a neural networks with Keras\n\nSet up a function to generate neural networks:"},{"metadata":{"trusted":true},"cell_type":"code","source":"library(keras)\nbuildModel <- function(NumberHiddenLayer = 1, UnitsHiddenLayer = 32, lr = 0.001, reg = 0.001, droprate = 0.5) {\n    kerModell <- keras_model_sequential() \n    layer_gaussian_noise(kerModell, stddev = 0.1, input_shape = dim(trainMatrix)[2])\n    layer_dense(kerModell, units = dim(trainMatrix)[2], activation = \"relu\", kernel_regularizer = regularizer_l2(l=reg))\n    layer_dropout(kerModell, rate = droprate)\n    if (NumberHiddenLayer >= 1) {\n        for (i in 1:NumberHiddenLayer) {\n            layer_dense(kerModell, units = max(4, UnitsHiddenLayer), activation = \"relu\", kernel_regularizer = regularizer_l2(l=reg))\n            layer_dropout(kerModell, rate = droprate)}\n    }\n    layer_dense(kerModell, units = max(2, UnitsHiddenLayer/2), activation = \"linear\", kernel_regularizer = regularizer_l2(l=reg))\n    layer_dropout(kerModell, rate = droprate)\n    layer_dense(kerModell, units = 1, activation = \"linear\")\n    compile(kerModell, loss = \"mse\", optimizer = optimizer_rmsprop(lr = lr, decay = 1e-6), metrics = list(\"mean_absolute_error\") )\n    kerModell\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Build a list of models for the final network including 5fold cross-validation. Train these models and show the evolution of the MSE/MAE for training and validation data each:"},{"metadata":{"trusted":true},"cell_type":"code","source":"set.seed(1)\nTrainVal <- caret::createFolds(trainTarget, k = 5)\nEQ_Modell_Keras <- list()\nearly_stop <- callback_early_stopping(monitor = \"val_loss\", patience = 10)\nfor (i in 1:5) {\n    EQ_Modell_Keras[[i]] <- buildModel(NumberHiddenLayer = 3, UnitsHiddenLayer = 128, lr = 0.001, reg = 0.050, droprate = 0.35)\n    trainRows <- unlist(TrainVal[1:5 != i])\n    valRows <- TrainVal[[i]]\n    fit_history <- fit(EQ_Modell_Keras[[i]], \n                       trainMatrix[trainRows, ], \n                       trainTarget[trainRows],\n                       epochs = 100,\n                       validation_data = list(trainMatrix[valRows,], \n                                              trainTarget[valRows]),\n                       callbacks = list(early_stop),\n                       verbose = 0)\n    plot(fit_history, main = paste(\"MSE and MAE-history for keras model\", i))\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Build a LightGBM-Model\n\nDirectly build 5 models for CV - similar to Keras-NN above. List the MSE and best iteration for the 5 \"submodels\":"},{"metadata":{"trusted":true},"cell_type":"code","source":"library(lightgbm)\nEQ_Modell_LGBM <- list()\nlgbParams <- list(objective = \"regression\", metric = \"l2\")\nfor (i in 1:5) {\n    trainRows <- unlist(TrainVal[1:5 != i])\n    valRows <- TrainVal[[i]]\n    actTrain <- lgb.Dataset(trainMatrix[trainRows, ], label = trainTarget[trainRows], free_raw_data = FALSE)\n    actValid <- lgb.Dataset.create.valid(actTrain, trainMatrix[valRows, ], label = trainTarget[valRows])\n    EQ_Modell_LGBM[[i]] <- lgb.train(params = lgbParams, data = actTrain, valids = list(train = actTrain, test = actValid), obj = \"regression\",\n                                     nrounds = 1000, learning_rate = 0.1, min_data_in_leaf = 20, num_leaves = 255,\n                                     verbose = 0, eval_freq = 50, early_stopping_rounds = 10)\n    cat(paste(\"LGBM Model\", i, \"best score:\", round(EQ_Modell_LGBM[[i]]$best_score, 3), \"Best iteration:\", EQ_Modell_LGBM[[i]]$best_iter, \"\\n\"))\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Predict the submission data\n\nThe test data comes as single files for each batch - so generate a file list and submission-dataframe first, then read the files one by one, predict from the five models for each file, and add to the submission-dataframe for keras, lightGBM and a combination of both - by using the median of the single predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"FileList <- list.files(\"../input/LANL-Earthquake-Prediction/test\", \".csv\")\nsubmission_keras <- submission_lgbm <- submission_combi <- data.frame(seg_id = substr(FileList, 1, 10), time_to_failure = rep(0, length(FileList)))\nfor (i in 1:length(FileList)){\n    testInput <- fread(paste0(\"../input/LANL-Earthquake-Prediction/test/\", FileList[i]), header = TRUE)\n    testInput <- matrix(EQ_Aggregate(testInput$acoustic_data, sdInput, pInput, mInput), nrow = 1)\n    testInput <- scale(testInput, center = colMeansTrain, scale = colStdAbwTrain)\n    testInput <- matrix(testInput[, colStdAbwTrain != 0], nrow = 1)\n    testPredict_keras <- sapply(EQ_Modell_Keras, FUN = function(x) max(c(0, predict(x, testInput))))\n    submission_keras[i, 2] <- median(testPredict_keras)\n    testPredict_lgbm <- sapply(EQ_Modell_LGBM, FUN = function(x) max(c(0, predict(x, testInput))))\n    submission_lgbm[i, 2] <- median(testPredict_lgbm)\n    submission_combi[i, 2] <- median(c(testPredict_keras, testPredict_lgbm))\n}\nhead(cbind(submission_keras, submission_lgbm[, 2], submission_combi[, 2]))\nwrite.table(submission_keras, file = \"submission_keras.csv\", quote = FALSE, sep = \",\", row.names = FALSE)\nwrite.table(submission_lgbm, file = \"submission_lgbm.csv\", quote = FALSE, sep = \",\", row.names = FALSE)\nwrite.table(submission_combi, file = \"submission_combi.csv\", quote = FALSE, sep = \",\", row.names = FALSE)","execution_count":null,"outputs":[]}],"metadata":{"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"nbformat":4,"nbformat_minor":1}