{"cells":[{"metadata":{"_uuid":"7e47797c8cf1b5ff170e5094a2a65f4ae35b024e","_execution_state":"idle","trusted":false},"cell_type":"markdown","source":"# Easy earthquake prediction with random forests\nThis is a simple approach for predicting the next earthquake based on basic feature extraction. I am using the kernel to evaluate the effect of overlapping chunk windows during training for the discussion in https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/79247. \n\nWe will use random forest since they are easy to tune and work well with a lot of highly correlated features. \n\nSo let's start by loading a couple of packages, defining some functions as well as constants that we will use later."},{"metadata":{"trusted":true,"_uuid":"6a668801131d6fb854954ab6869f5036ba5658b6"},"cell_type":"code","source":"#===============================================================\n# PACKAGES\n#===============================================================\n\nlibrary(data.table)\nlibrary(roll)\nlibrary(ranger)\nlibrary(tidyverse)\n\n#===============================================================\n# FUNCTIONS\n#===============================================================\n\n# Mean absolute error\nmae <- function(y, pred) {\n  mean(abs(y - pred))\n}\n\n# Calculates the coefficient of an AR(1) process z\nar1 <- function(z) {\n  cor(z[-length(z)], z[-1])  \n}\n\n# Calculates a buntch of statistics from vector x\nunivariate_stats <- function(x, tag = NULL, p = c(0, 0.25, 0.75, 1)) {\n  x <- x[!is.na(x)]\n \n  out <- c(\n    mean = mean(x),\n    sd = sd(x),\n    setNames(quantile(x, p = p, names = FALSE), paste0(\"q\", p)),\n    ar1 = ar1(x))\n  \n  if (is.null(tag)) \n    return(out)\n  \n  names(out) <- paste(names(out), tag, sep = \"_\")\n  out\n}\n\n# Feature extraction on vector x. Basically calls \"univariate_stats\" on differently transformed x\ncreate_X <- function(x, rolling_windows = c(10, 1000)) {\n  stats_full <- univariate_stats(x = x, \"full\", p = c(0, 1, 5, 10, 25, 50, 75, 90, 95, 99, 100) / 100)\n  stats_abs <- univariate_stats(x = abs(x), \"abs\")\n\n  # Rolling versions of x\n  x_mat <- as.matrix(x, ncol = 1)\n  roll_sd_k <- lapply(rolling_windows, function(k) roll_sd(x_mat, width = k))\n  \n  # Derive stats from rolling versions\n  stats_roll_sd <- Map(univariate_stats, roll_sd_k, tag = paste(\"roll_sd\", rolling_windows, sep = \"_\"))\n  \n  c(stats_full, stats_abs, unlist(stats_roll_sd))\n}\n\n#===============================================================\n# CONSTANTS\n#===============================================================\n\n# Length of test data sets\nn_test <- 150000\n\n# By how much do we shift the time window of n_test rows within earthquake?\nstride <- 75000\n\n# Positions of earthquakes (see answers in https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/77390). Used to create contiguous folds for cross-validation and train/validation split\nearthquakes <- c(\n    5656573,\n   50085877,\n  104677355,\n  138772452,\n  187641819,\n  218652629,\n  245829584,\n  307838916,\n  338276286,\n  375377847,\n  419368879,\n  461811622,\n  495800224,\n  528777114,\n  585568143,\n  621985672) + 1","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"01a7e7bd0c0ac503b69f2400e3118693228fb776"},"cell_type":"markdown","source":"So let's create one data set per earthquake as follows:\n\n1. Load all rows before an earthquake occurs.\n\n2. For chunks of size 150'000 (like test data), we extract a couple of features and store them in a row of a matrix X. The response is stored in a vector y. \n\n3. We move by \"stride\" positions and repeat steps 2 & 3 until the earthquake happens."},{"metadata":{"trusted":true,"_uuid":"95706522f9b033260878ae2d4203ac14b5f64521"},"cell_type":"code","source":"# Figure out colnames of input as well as number of features\nraw <- fread(file.path(\"..\", \"input\", \"train.csv\"), nrows = 150000, data.table = FALSE)\nnames_input <- names(raw)\nnames_features <- names(create_X(raw[[1]]))\n\n# If not yet created, make directory to save folds\nfold_dir <- file.path(\"strides\", stride)\nif (!dir.exists(fold_dir)) {\n  dir.create(fold_dir, recursive = TRUE)  \n}\n\nfor (prep_fold in seq_along(earthquakes)) {# prep_fold <- 1\n cat(prep_fold, \"\\n\")\n  \n  # Read data between two earthquakes\n  raw <- fread(file.path(\"..\", \"input\", \"train.csv\"), \n               nrows = c(earthquakes[1], diff(earthquakes))[prep_fold], \n               skip = c(0, earthquakes)[prep_fold] + 1)\n  setnames(raw, names_input)\n  \n  # How many times do we calculate features for this data chunk?\n  n_steps <- (nrow(raw) - n_test) %/% stride\n  \n  # Init feature matrix and vector of response\n  y <- numeric(n_steps)\n  X <- matrix(NA, nrow = n_steps, ncol = length(names_features), dimnames = list(NULL, names_features))\n  \n  # Loop through chunk and build up y and X\n  pb <- txtProgressBar(0, n_steps, style = 3)\n  \n  for (i in seq_len(n_steps)) {\n    setTxtProgressBar(pb, i)\n    from <- 1 + stride * (i - 1)\n    to <- n_test + stride * (i - 1)\n    X[i, ] <- create_X(raw$acoustic_data[from:to])\n    y[i] <- raw$time_to_failure[to]\n  }\n\n  save(y, X, file = file.path(fold_dir, paste0(\"fold_\", prep_fold, \".RData\")))\n}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3009326c561692722095c126f1e4724dd88e3976"},"cell_type":"markdown","source":"Now, the data preparation is over and we can load all saved .RData files to create the full data set. Each .RData (= each earthquake) will serve as one fold in our leave-one-earthquake-out cross-validation strategy."},{"metadata":{"trusted":true,"_uuid":"0c86fc5e7ab7eba50affef35e4bebcc33919ca8d"},"cell_type":"code","source":"for (i in seq_along(earthquakes)) {\n  load(file.path(\"strides\", stride, paste0(\"fold_\", i, \".RData\")))\n  \n  fold <- rep(i, length(y))\n  \n  if (i == 1) {\n    X_mat <- X\n    y_vec <- y\n    fold_vec <- fold \n  } else {\n    X_mat <- rbind(X_mat, X)\n    y_vec <- c(y_vec, y)\n    fold_vec <- c(fold_vec, fold)\n  }\n}","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9d742609052434aaba89d37abcc7ad82b327edde"},"cell_type":"markdown","source":"Now, we will use above mentioned cross-validation strategy to evaluate the quality of our random forest model. "},{"metadata":{"trusted":true,"_uuid":"92bfad8eb7b3c240735f72f462e975549ada9c12"},"cell_type":"code","source":"form <- reformulate(colnames(X_mat), \"label\")\nfullDF <- data.frame(label = y_vec, X_mat)\n\n# cross-validation\nm_fold <- length(earthquakes)\ncv <- numeric(m_fold)\npb <- txtProgressBar(0, m_fold, style = 3)\n\nfor (j in seq_along(cv)) { # j <- 1\n  setTxtProgressBar(pb, j)\n  fit <- ranger(form, fullDF[fold_vec != j, ], seed = 3564 + 54 * j, verbose = 0)\n  cv[j] <- mae(fullDF[fold_vec == j, \"label\"], predict(fit, fullDF[fold_vec == j, ])$predictions)\n}\n\n# Resulting score\nmean(cv)\nweighted.mean(cv, w = tabulate(fold_vec, nbins = m_fold))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"daae25c27374b7d9cd5fe41e82a9dd876621350a"},"cell_type":"markdown","source":"We are almost done. Just retrain the random forest on the full data set for submission and have a quick look at the variable importance."},{"metadata":{"trusted":true,"_uuid":"afb60b6dcfe8e74d20e83a5603f64122c9ea5167"},"cell_type":"code","source":"# retrain on full data for submission\nfit_rf <- ranger(form, fullDF, importance = \"impurity\", seed = 345)\n\n# Variable importance\npar(mar = c(5, 10, 1, 1))\nbarplot(importance(fit_rf) %>% sort %>% tail(60), horiz = T, las = 1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a9f6434b6149da1a3fc6563992b36c41fa08cc24"},"cell_type":"markdown","source":"Preparing the data for submission is slow. As long as you do not change the features (but e.g. only the stride or the fitting algo), you can reload the prepared data `all_test2` from disk and skip the painful step."},{"metadata":{"trusted":true,"_uuid":"9e29689bd481d334d3f6f45d4d0a98de50bfd219"},"cell_type":"code","source":"# submission\nsubmission <- fread(file.path(\"..\", \"input\", \"sample_submission.csv\"))\n\n# Load each test data and create the feature matrix. Takes 2-3 minutes and can be skipped\n# if playing with strides (and not with features)\nload_and_prepare <- function(file) {\n  seg <- fread(file.path(\"..\", \"input\", \"test\", paste0(file, \".csv\")))\n  create_X(seg$acoustic_data)\n}\nall_test <- lapply(submission$seg_id, load_and_prepare)\nall_test2 <- do.call(rbind, all_test)\nsave(all_test2, file = \"test_prep.RData\")\n# load(\"test_prep.RData\") \ndim(all_test2) # 2624   35\n\nsubmission$time_to_failure <- predict(fit_rf, data.frame(all_test2))$prediction\n\nhead(submission)\n\n# Save\nfwrite(submission, paste0(\"submission_rf_stride_\", stride, \".csv\"))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a168ed699801a74edfbeadcd71018d788ee1cfbf"},"cell_type":"markdown","source":""}],"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}