{"cells":[{"metadata":{},"cell_type":"markdown","source":"# R for LANL 2019 feature generation: basic R vs. data.table\n\nThe LANL 2019 was my first Kaggle competition - as a beginner I struggled with some basic topics like \"how to load and calculate with that amount of data?\". Here I want to show a short comparison of base-R and data.table for loading the training data and two approaches how to get reasonable calculation speed.\n\n## Define function for Data-aggregation / feature generation\n\nAggregate 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\n\n### Read with base-R \n\nLoad the data: trying to read the whole file at once leads to a kernel shutdown - so only the first 6E7 rows (~10% of the overall data) are loaded timing test - this already shows the speed-difference."},{"metadata":{"trusted":true},"cell_type":"code","source":"system.time({\n    trainEQ <- read.csv(\"../input/LANL-Earthquake-Prediction/train.csv\", header = TRUE, nrows = 6E7, colClasses = c(\"integer\", \"numeric\"))\n})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The data is loaded and takes ~2GB RAM - this is critical for 18GB RAM being only 10% of the overall data - the rason why the kernal shuts down when trying to load all the data.\n\n### Read with data.table\n\nWhen parsing csv-Data data.table is usually the fastest way - compare to the speed above:"},{"metadata":{"trusted":true},"cell_type":"code","source":"library(data.table)\nsystem.time(trainEQ <- fread(\"../input/LANL-Earthquake-Prediction/train.csv\", header = TRUE))\ncat(paste(\"number of rows:\", nrow(trainEQ)))\nsystem.time({# first: inputevaluate mean, peak, sd. Then: output the number of parameters.\n    mInput <- mean(trainEQ$acoustic_data)\n    pInput <- max(abs(trainEQ$acoustic_data))\n    sdInput <- sd(trainEQ$acoustic_data)\n    nParameter <- length(EQ_Aggregate(trainEQ$acoustic_data[1:testBatchsize], sdInput, pInput, mInput))\n})\ncat(paste(\"number of aggregated parameters:\", nParameter))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Confirmed: reading the data is way faster (same order of magnitude as it took with read.csv for 3% of the data) - actually with that amount of data the limits of the base-R read.csv in combination with the available RAM for this kernel are almost reached.\n\nContinue with the data.table-Version and prepare some limits for first/last row to use in the following:"},{"metadata":{"trusted":true},"cell_type":"code","source":"set.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 with data.table\n\nWhen working on the \"Predicting Molecular Properties\" 2019 competition I found data.table very helpful for the data-aggregation required there. \n\nDoing all the calculation in a single step results in too high RAM-consumption. I had to split the calculation in steps of 50 batches (1.5E5 rows for each batch - 500 and 1000 was still to high for running as a kernal) to limit the RAM-consumption and thus avoiding kernel-crashes.\n\nSo let's check how it works with these restrictions:"},{"metadata":{"trusted":true},"cell_type":"code","source":"system.time({\n    trainEQ[, batchNum:=as.integer((1:nrow(trainEQ)-firstRow) %/% testBatchsize)]\n})\nsystem.time({\n    trainMatrix <- matrix(NA, nrow = nBatches, ncol = nParameter)\n    for (i in 0:(nBatches %/% 50)) {\n        fRow <- firstRow + i*50*testBatchsize\n        lRow <- min(fRow + 50*testBatchsize - 1, lastRow + testBatchsize - 1)\n        tmpMatrix <- trainEQ[fRow:lRow, as.list(EQ_Aggregate(acoustic_data, sdInput, pInput, mInput)), batchNum]\n        trainMatrix[seq(from = i*50 + 1, length.out = nrow(tmpMatrix)), ] <- as.matrix(tmpMatrix[, 2:(nParameter + 1)])\n    }\n    trainEQ[,batchNum:=NULL]; rm(fRow, lRow, tmpMatrix)\n})\ncat(paste(\"dimensions of the final matrix:\", dim(trainMatrix)[1], \"rows, \", dim(trainMatrix)[2], \"columns.\"))\ntrainMatrix_DT <- trainMatrix","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Feature generation in for-loop\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.\n\nGeneral recommendation that is included in the loop: generate the empty objects with appropriate dimensions (or larger) first and fill them later.\n"},{"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})\ncat(paste(\"dimensions of the final matrix:\", dim(trainMatrix)[1], \"rows, \", dim(trainMatrix)[2], \"columns.\"))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The data.table way takes almost the same time as the loop - if you don't take into account the a column with the batch number has to be generated in the beginning (and removed later). But: we are talking about <20s for that column - that's not really critical. However, getting the calculation stable enough for the kernel was the real problem with the data.table-approach.\n\nFinally: a quick sanity check - do both approaches yield equal results? Sum up all different entries (should be 0):"},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"sum(trainMatrix_DT != trainMatrix)\nrm(trainMatrix_DT)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looks OK: all entries are equal.\n\n## 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":"system.time({\n    trainTarget <- vector(mode = \"double\", length = nBatches)\n    removeRow <- NULL\n    l <- 0\n    for (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    }\n    trainMatrix <- trainMatrix[-removeRow, ]\n    trainTarget <- trainTarget[-removeRow]\n})\ncat(paste(\"number of removed rows:\", length(removeRow)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Summary Feature generation\n\nFor reading GB-size data: use data.table. Depending on which calculation you are planning to do both a simple loop with well prepared objects for the final data and data.table can be fastest - here it was the for-loop. Keep in mind the ratio of RAM and data-size for calculation - parallel calculation only makes sense if there's enough RAM available.\n\n## Quick EDA for all parameters\n\nQucik exploratory data analysis: just check the output parameters for impact on the time to failure:\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"par(mfrow = c(2,2))\nfor (i in 1:ncol(trainMatrix)) {\n    plot(trainTarget~trainMatrix[, i], xlab = paste(\"parameter\", i), ylab = \"time to failure\")\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Summary EDA\n\nAlready with these example feautures one can see some correlations to the time to failure."}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}