{"cells":[{"metadata":{},"cell_type":"markdown","source":"**........... Please cite and acknowledge if you use this code in scientific publications..............\n**\n\nThis kernel has the feature generation and various cross-validation strategies with gradient boosters (LightGBM, Catboost, XGBoost) implemented in R. The features here gave the best cross-validation score and were part of our 10th place solution. The model training used in our solution was implemented by Giba (see his kernel), but I'm including my functions here and their usage in case they are useful to anybody."},{"metadata":{"trusted":false},"cell_type":"code","source":"library(devtools)\ninstall_github(\"andrewuhl/RollingWindow\")\nlibrary(RollingWindow)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"install.packages(\"tsfeatures\")","execution_count":null,"outputs":[]},{"metadata":{"_execution_state":"idle","_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","trusted":false},"cell_type":"code","source":"library(caret)\nlibrary(catboost)\nlibrary(tictoc)\nlibrary(tsfeatures)\nlibrary(rlist)\nlibrary(furrr)\nlibrary(wmtsa)\nlibrary(RcppRoll)\nlibrary(tuneR)\nlibrary(seewave)\nlibrary(Matrix)\nlibrary(tidyverse) \nlibrary(data.table)\nlibrary(lightgbm)\n\n#' count the number of peaks\n#' @param x the input numeric vector of accoustic data\n#' @param window the window size to use\n#' @return a scalar numeric for the number of detected peaks\n#' @export\nnumber_peaks <- function(x, window){\n    rollmax <- as.numeric(RollingMax(x, window))\n    npeak <- sum(x > lead(rollmax, n = window) & x > lag(rollmax), na.rm = T)\n    data.frame(npeak) %>% set_names(paste0(\"num_peaks_\", window))\n}\n\n#' count the number of peaks with several windows\n#' @param x the input numeric vector of accoustic data\n#' @param windows the window sizes to use\n#' @return a data frame with the number of peaks at different windows threshold\n#' @export\npeakFeatures <- function(x, windows = c(10, 50, 100, 500)){\n    map_dfc(windows, ~number_peaks(x, window=.x))\n}\n\n#' Collect statistical features from accoustic numeric data\n#' @param x the input numeric vector of accoustic data\n#' @param quant the quantiles to compute\n#' @return a data frame of mean, sd, min, max, IQR, RMS and selected quantiles\n#' @export\ncollectStats <- function(x, quant = c(0.05, 0.1,0.2,0.5,0.8, 0.95,0.9)){\n    data.frame(t(quantile(x,quant))) %>%\n      set_names(paste0(\"Q\",quant))   %>%\n      mutate(mean = mean(x), sd = sd(x), \n             min  = min(x), max = max(x), \n             IQR  = IQR(x), RMS = sqrt(mean(x^2)))\n}\n\n#' Compute number of crossing points at the mean and selected quantiles devided by segment length \n#' @param x the input numeric vector of accoustic data\n#' @param q the quantiles at which to compute crossing points\n#' @return a data frame with computed fraction of crossing points\n#' @export\ncrosspoints <- function(x, q = c(0.01,0.05, 0.1,0.15, 0.2,0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8,0.85, 0.9, 0.95, 0.99)){\n    require(dplyr)\n    require(purrr)\n    midlines <- c(mean(x),as.numeric(quantile(x,q)))\n    out <- map_dfc(midlines, function(midline,x){\n        lagx <- lag(x, default = 0)\n        sum(x > midline & lagx <= midline  | x <= midline & lagx > midline)/length(x)},\n        x = x)\n    colnames(out) <- paste0(\"crossPoints_\",c(\"M\",q))\n    out\n}\n\n#' Computes spectrum features: the mean and variance at frequency bands from start to end, 1 to start, and end to 75000\n#' @param x the input numeric accoustic data\n#' @param start first frequency\n#' @param end last frequency\n#' @param by the step size between start and end\n#' @return a data frame with spectrum features\n#' @export\nspecFeatures <- function(x, start = 1000, end = 20000, by = 1000){\n    s <- spectrum(x, plot = F)\n    spec <-  s$spec\n    bandStart <- seq(start, end, by)\n    bandEnd <- bandStart + by - 1\n    bandStart <- c(1, bandStart, end + by + 1)\n    bandEnd <- c(start - 1, bandEnd, 75000)\n    features <- map2_dfc(bandStart, bandEnd,\n                  ~ data.frame(meanBand = log(mean(spec[.x:.y])),\n                               varBand = log(var(spec[.x:.y]))))\n    return(features)                           \n}\n\n\n\n#' The feature extraction engine. It first denoise the data, then extract the features\n#' @param x the input numeric accoustic data\n#' @param time time to failure input to be used with training data, keep NULL with test data\n#' @param train whether the input is training or testing data\n#' @param wavelet the wavelet to be used for both denoising and decomposition\n#' @return a data frame with extracted features\n#' @export\nextractFeatures <- function(x, time = NULL, train = T, wavelet = \"l14\" ){\n    # first denoise the data. There are thresholding options here. I'm using the default\n    # the choice of wavelet length is important\n    x <- wavShrink(x, wavelet = wavelet)\n    \n    # Give segments with time reset a zero ttf to be removed later\n    if(train){\n      if( sum(diff(time)>0)){\n        time <- 0 \n      }\n    }\n    \n    #xsplit <- splitSegment(x)\n\n    # decomposition of the signal.\n    dwt1 <- wavMODWT(x, n.levels = 11, wavelet = \"l14\")\n    dwt2 <- wavMODWT(x, n.levels = 11, wavelet = \"haar\")\n    dwt3 <- wavMODWT(x, n.levels = 11, wavelet = \"s6\")\n    # collect stats on the signal, derivative of it, absolute value of it, and its decomposed parts\n    # I'm skipping dwt$data$d1 here because it seems to be mostly noise\n    \n    \n    # statistical features were omitted in our solution\n    #features1 <- map_dfc(list(x, diff(x),abs(x),\n    #                         dwt1$data$d2, dwt1$data$d3,dwt1$data$d4, dwt1$data$d6, dwt1$data$d8, dwt1$data$d10,\n    #                         dwt2$data$d2, dwt2$data$d3,dwt2$data$d4, dwt2$data$d6, dwt2$data$d8, dwt2$data$d10,\n    #                         dwt3$data$d2, dwt3$data$d3,dwt3$data$d4, dwt3$data$d6, dwt3$data$d8, dwt3$data$d10), collectStats)\n    \n    # collect some features inspired by time-series analysis\n    # auto-correlation and partial autocorrelation in addition to entropy and stability\n    features2 <- tsfeatures(as.ts(x), \n                    features = c(\"entropy\", \"acf_features\", \"stability\", \"pacf_features\"))\n\n    # compute crossing points on the signal and its decomposed parts\n    features3 <- map_dfc(list(x, dwt1$data$d2, dwt1$data$d3,dwt1$data$d4, dwt1$data$d6, dwt1$data$d8, dwt1$data$d10,\n                                 dwt2$data$d2, dwt2$data$d3,dwt2$data$d4, dwt1$data$d6, dwt2$data$d8, dwt2$data$d10,\n                                 dwt3$data$d2, dwt3$data$d3,dwt3$data$d4, dwt1$data$d6, dwt3$data$d8, dwt3$data$d10), crosspoints)\n    \n    # compute spectrum features on the signal\n    features4 <- map_dfc(list(x),specFeatures)\n    \n    \n    features5 <- map_dfc(list(x),peakFeatures)\n    # join everything together\n    features <- data.frame(features2, features3, features4, features5)\n    \n    # annotate ttf if the input is training data\n    if(train){\n        features <- features %>% mutate(ttf = time[length(time)])\n    }\n    return(features)\n}\n\n#' Computes mean absolute error\n#' @param x real value\n#' @param pred predicted value\n#' @return a scalar numeric\n#' @export\nmae = function(x,pred){\n    return(mean(abs(x-pred)))\n}\n\n#' Utility function for staged prediction for catboost\n#' @param model catboost model\n#' @param test test data\n#' @return a matrix prediction for each iteration\n#' @export\nstagedPredict <- function(model, test){\n    ntree <- model$tree_count\n    walker <- catboost.staged_predict(model, \n         test, verbose = FALSE, \n         prediction_type = \"RawFormulaVal\", \n         ntree_start = 0, ntree_end = 0, \n         eval_period = 1, thread_count = -1)\n    p <- c()\n    for(i in 1:ntree){\n        p <- rbind(p, walker$nextElem())\n    }\n    return(p)\n}\n\n#' Extracting features with augmentation. Each two consecutive segments overlap, but not with the next two (useful for shuffling)\n#' @param signal the accoustic signal\n#' @param time time to failure\n#' @return data frame with two rows with the extracted features from two partially overlapping segments with the same ID\n#' @export\nextractFeaturesAugment <- function(signal, time){\n    map2_dfr(list(head(signal, 150000), tail(signal, 150000)),\n             list(head(time, 150000), tail(time, 150000)), ~ extractFeatures(.x, .y)) %>%\n    mutate(ID = c(1,2))\n}\n\n#' train data reading engine. Reads the data and run ExtractFeatures on segments of 150k\n#' @param skip the number of rows to skip from the input file\n#' @param nmax the number of rows to read. Should be a multiple of segment size\n#' @param segmentSize the size of the segment\n#' @param augment whether to use stride in extracting features, segmentSize should be more than 150k.\n#' @return a data frame with features \n#' @export\nprocessData <- function(skip, nmax, segmentSize = 150000, augment = F){\n  if(augment){\n     fread(\"../input/LANL-Earthquake-Prediction/train.csv\", skip = skip, \n          nrows = nmax) %>% set_names(c(\"signal\",\"time\")) %>%\n          mutate(idx = rep(c(1:(nmax/segmentSize)),each = segmentSize)) %>%\n          group_by(idx) %>% do(extractFeaturesAugment(.$signal, .$time)) %>%\n          ungroup() %>% select(-idx) %>% as.data.frame()\n  } else {\n     fread(\"../input/LANL-Earthquake-Prediction/train.csv\", skip = skip, \n          nrows = nmax) %>% set_names(c(\"signal\",\"time\")) %>%\n          mutate(idx = rep(c(1:(nmax/segmentSize)),each = segmentSize)) %>%\n          group_by(idx) %>% do(extractFeatures(.$signal, .$time)) %>%\n          ungroup() %>% select(-idx) %>% as.data.frame()\n\n  }\n}\n\n#' test data reading engine. Reads test data segments and run ExtractFeatures\n#' @param ID the ID of the test segment\n#' @return data frame with extracted features\n#' @export\nbuildTestSet <- function(ID){\n    test <- fread(paste0(\"../input/LANL-Earthquake-Prediction/test/\",ID,\".csv\")) %>%\n            select(signal = acoustic_data)\n    return(extractFeatures(test$signal, train = F))\n}\n\n#' A generic function to build data pools for LGB, catboost and XGB models\n#' @param data data frame of input features\n#' @param label the label to predict\n#' model the model type can be catboost, lgbm or xgboost\n#' @return a suitable data pool to the model type\n#' @export\nmakePool <- function(data, label, model = c(\"catboost\",\"lgbm\",\"xgboost\")){\n    model <- match.arg(model)\n    if(model == \"catboost\"){\n      catboost.load_pool(data = data, label = label)\n    } else if(model == \"lgbm\"){\n       lgb.Dataset(Matrix(as.matrix(data),sparse=T), label = label)\n    } else if(model == \"xgboost\"){\n       xgb.DMatrix(data = as.matrix(data), label = label)\n    }\n}\n\n#' Internal function to fit a model on one fold in cross-validation setting.\n#' The function scales the training folds and apply the scaling on the validation fold\n    \n#' @param folds a list of data frames with each is a fold\n#' @param params a list of training paramters. Should match the model type\n#' @param index the index of the validation fold\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#'                Currently only works with LGBM. Set to FALSE with non-shuffled CV\n#' @return a list with the folloing elements: 'model' the fitted model, 'pred' validation fold predictions and\n#'         'trainData' a data frame for the scaled training folds, to be used to scale in later predictions\n#' @export    \nevalFold <- function(folds, params, index, modelType = c(\"catboost\",\"lgbm\",\"xgboost\"), useBest = T){\n    modelType <- match.arg(modelType)\n    train <- bind_rows(list.remove(folds, index))\n    x_train <- scale(select(train, -ttf))\n    y_train <- train$ttf\n    valid <- folds[[index]]\n    \n    # scale the data\n    x_valid <- scale(select(valid, -ttf),\n        center = attr(x_train,\"scaled:center\"),\n        scale = attr(x_train, \"scaled:scale\"))\n    y_valid <- valid$ttf\n   \n    if(modelType == \"catboost\"){\n      validFold <- makePool(data = x_valid, label = y_valid, model = modelType)\n      trainFold <- makePool(data = x_train, label = y_train, model = modelType)\n      model <- catboost.train(trainFold,validFold, params = params)\n      pred <- stagedPredict(model, validFold)\n      pred <- as.numeric(pred[nrow(pred),])\n      mae <- mae(pred,folds[[index]]$ttf)\n    \n    } else if(modelType == \"lgbm\"){\n      trainFold <- lgb.Dataset(Matrix(as.matrix(x_train),sparse=T), label = y_train)\n      validFold <- lgb.Dataset.create.valid(trainFold, data = Matrix(as.matrix(x_valid),sparse=T), label = y_valid)\n      valids <- list(eval = validFold)\n      model <- lgb.train(params, trainFold, params$num_iterations, valids, verbose = F)\n      \n      if(useBest){\n          pred <- predict(model, Matrix(as.matrix(x_valid),sparse=T))\n      } else {\n         pred <- predict(model, Matrix(as.matrix(x_valid),sparse=T),num_iteration = model$current_iter()) \n      }\n      mae <- mae(pred,folds[[index]]$ttf)\n        \n    } else if(modelType == \"xgboost\"){\n      trainFold <- makePool(data = x_train, label = y_train, model = modelType)\n      validFold <- makePool(data = x_valid, label = y_valid, model = modelType)\n      watchlist <- list(train = trainFold, test = validFold)\n      nrounds <- params$nrounds\n      capture.output(model <- xgb.train(params, trainFold, watchlist = watchlist, nrounds = nrounds, \n                         early_stopping_round = params$early_stopping_round), file = 'NUL')\n      pred <- predict(model, validFold)\n      mae <- mae(pred, folds[[index]]$ttf)\n    }    \n    cat(paste(\"MAE for fold\",index,\"is\",mae,\"\\n\"))\n    return(list(model = model, pred = pred, trainData = x_train))\n}\n\n#' The cross validation engine function. Loops over folds and calls evalFold. Prints MAE for each fold,\n#'    the whole training data, and long and short quakes\n#' @param folds a list of data frames each is a fold \n#' @param params a list of training paramters. Should match the model type\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param quake a vector of quake IDs    \n#' @param objective convenience paramter to pass the objective for training\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#' @return a list of 'models' is a list of models, 'pred' out of fold predictions, \n#'    'trainData' a list of data frames of input training data scaled according to which fold is left out\n#' @export    \ncrossValidate <- function(folds, params, modelType = c(\"catboost\",\"lgbm\", \"xgboost\"), quake, objective = \"\", useBest = T){\n    modelType <- match.arg(modelType)\n    if(objective != \"\"){\n        if(modelType == \"lgbm\"){\n            params$objective <- objective\n        } else if(modelType == \"catboost\"){\n            params$loss_function <- objective\n        } else if(modelType == \"xgboost\"){\n            params$objective <- objective\n        }\n    }\n    foldResult <- map(c(1:length(folds)),\n                      ~evalFold(folds, index = .x, params = params, modelType = modelType, useBest = useBest))\n    models <- map(foldResult, \"model\")\n    pred <- unlist(map(foldResult, \"pred\"))  \n    trainData <- map(foldResult, \"trainData\")\n    target <- bind_rows(folds)$ttf\n    MAE <- mae(pred, target)\n    \n    # this is how I am classifying the quakes. Hard-coded. \n    quakeType <- ifelse(quake %in% c(1,2,3,5,8,11,12,15,17), \"Type1\",\"Type2\")\n    maeType1 <- mae(pred[quakeType == \"Type1\"], target[quakeType == \"Type1\"])\n    maeType2 <- mae(pred[quakeType == \"Type2\"], target[quakeType == \"Type2\"])\n    cat(paste(\"final MAE is\",MAE,\"\\n\"))\n    cat(paste(\"final Long quake MAE is\",maeType1,\"\\n\"))\n    cat(paste(\"final Short quake MAE is\",maeType2,\"\\n\"))\n    return(list(models = models, pred = pred, trainData = trainData))\n}\n\n#' A generic prediction function that scales the data according to training data then apply suitable predict function\n#' @param model the fitted model\n#' @param trainData scaled training data to be used for applyin the scaling\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#'    Set to FALSE for non-shuffled CV\n#' @return a vector of predictions\n#' @export\n testPredict <- function(model, test, trainData, modelType = c(\"catboost\",\"lgbm\",\"xgboost\"), useBest = T){\n     modelType <- match.arg(modelType)\n     test <- scale(test, center = attr(trainData,\"scaled:center\"),\n        scale = attr(trainData, \"scaled:scale\"))\n     if(modelType == \"catboost\"){\n         testPool <- catboost.load_pool(data = test)\n         return(catboost.predict(model,testPool))\n     } else if(modelType == \"lgbm\"){\n         if(useBest){\n            return(predict(model, Matrix(as.matrix(test),sparse=T))) \n         } else {\n            return(predict(model, Matrix(as.matrix(test),sparse=T), num_iteration = model$current_iter())) \n         }\n         \n     } else if(modelType == \"xgboost\"){\n         return(predict(model, as.matrix(test)))\n     }\n }\n\n#' Creates shuffled stratified KFold indices\n#' @param data data frame of training data \n#' @param k number of folds\n#' @param seed the random seed\n#' @return a vector of fold indices \n#' @export         \ncreateKFolds <- function(data, k, seed){\n    byFold <- round(nrow(data)/k)\n    splits <- rep(c(1:byFold),each = k, length.out=nrow(data))\n    set.seed(seed)\n    folds <- as.numeric(unlist(map(split(splits,splits), ~ sample(1:(length(.x))))))\n    return(ifelse(folds>k, folds-k, folds))\n}\n\n#' The engine to run cross-validation strategy\n#' @param trainSet data frame containing training data. Must be sorted by ttf in each quake.\n#'     Each reset is assumed a new quake\n#' @param validationStrategy the validation strategy to be applied. \n#'     Possible values \"shuffledKFold\",\"quakeWise\",\"KFold\", \"repeatedKFold\"\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#'    Set to FALSE for non-shuffled CV\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param objective convenience paramter to pass the objective for training\n#' @param k number of folds\n#' @param seed the random seed     \n#' @param repeats the number of repeats if strategy is repeatedKFold\n#' @return a list of 'models' is a list of models, 'pred' out of fold predictions, \n#'    'trainData' a list of data frames of input training data scaled according to which fold is left out\n#'    The out-of-fold prediction in 'pred' is always sorted as the input data.frame regarding of the strategy\n#' @export         \nperformCV <- function(trainSet, validationStrategy = c(\"shuffledKFold\",\"quakeWise\",\"KFold\", \"repeatedKFold\"),useBest = T,\n                      modelType = c(\"lgbm\",\"catboost\",\"xgboost\"), params, objective, seed = 123, k = 5, repeats = 5){\n  modelType <- match.arg(modelType)\n  validationStrategy <- match.arg(validationStrategy)\n  \n  # count how many quakes in input training data\n  nQuake <- sum((diff(trainSet$ttf)>0)) + 1\n  \n  # get the quake IDs\n  quake <- unlist(map2(c(1:nQuake),diff(c(0,which(diff(trainSet$ttf)>0),nrow(trainSet))), ~ rep(.x, each = .y)))\n  \n  # just a running ID to be used for sorting\n  id <- c(1:nrow(trainSet))\n  \n  if(validationStrategy == \"shuffledKFold\"){\n     # split the input by quake\n     tr <- split(trainSet, quake)  \n     \n     # create shuffled KFold indices by quake. This achieves quake stratification\n     kIndex <- as.numeric(unlist(map(tr, ~ createKFolds(.x, k = k, seed = seed))))                      \n     \n     # split the data by folds\n     train <- split(trainSet, kIndex)\n     quakeID <- unlist(split(quake, kIndex))\n     shuffledID <- unlist(split(id, kIndex))\n     \n     # run cross validation\n     cv <- crossValidate(train, params = params, modelType = modelType,\n              quake = quakeID, objective = objective)\n     \n     # sort oof predictions to match the input\n     cv$pred <- cv$pred[order(shuffledID)]\n     \n  } else if(validationStrategy == \"quakeWise\"){\n     # just split by quake and run cross validation\n     train <- split(trainSet, quake)\n     cv <- crossValidate(train, params = params, modelType = modelType,\n              quake = quake, objective = objective, useBest = useBest)\n              \n  } else if(validationStrategy == \"KFold\") {\n     # split the data by normal KFold without shuffling\n     kIndex <- rep(c(1:k),each = round(nrow(trainSet)/k), length.out = nrow(trainSet))\n     train <- split(trainSet, kIndex)\n     cv <- crossValidate(train, params = params, modelType = modelType,\n              quake = quake, objective = objective, useBest = useBest)\n  \n  } else if(validationStrategy == \"repeatedKFold\"){\n     # KFold repeated N times\n     # each repeat has the fold indices shifted\n     kIndex <- rep(c(1:k),each = round(nrow(trainSet)/k), length.out = nrow(trainSet)) \n     shifter <- function(x, n = 1) {\n         if (n == 0) x else c(tail(x, -n), head(x, n))\n     }\n     shiftLag <- round((nrow(trainSet)/k)/repeats)\n     kIndices <- map(c(0:(repeats-1)), ~shifter(kIndex, .x * shiftLag))\n     trainSetSplits <- map(kIndices, ~split(trainSet, .x))\n     idSplits <- map(kIndices, ~ unlist(split(id, .x)))\n     quakeSplits <- map(kIndices, ~ unlist(split(quake, .x)))\n     cv <- map2(trainSetSplits,quakeSplits, ~ crossValidate(.x, params = params, modelType = modelType,\n              quake = .y, objective = objective, useBest = useBest))\n     cv <- map2(cv, idSplits,  function(x,y){x$pred <- x$pred[order(y)]; x})        \n  }\n  return(cv)\n}\n\n#' A generic prediction function that can be applied on performCV output\n#' @param cvlist the output from performCV\n#' @param test the test set for prediction\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#'    Set to FALSE for non-shuffled CV\n#' @param repeated whether performCV was run using repeatedKFold strategy or not    \n#' @return a vector of predictions\n#' @export\nCVpredict <- function(cvlist, test, repeated = F, modelType = c(\"lgbm\",\"catboost\",\"xgboost\"), useBest = T){\n    modelType <- match.arg(modelType)\n    if(repeated){\n        rowMeans(map_dfc(cvlist, ~ rowMeans(map2_dfc(.x$models, .x$trainData,\n                 ~ testPredict(model =.x, trainData = .y, modelType = modelType, test =  test, useBest = useBest)))))\n    } else {\n        rowMeans(map2_dfc(cvlist$models, cvlist$trainData,\n                 ~ testPredict(model =.x, trainData = .y, modelType = modelType, test =  test, useBest = useBest)))\n    }\n}\n\n#' A function to plot distribution of a feature between train and test stratified by quake\n#' @param train training set data frame\n#' @param test test set data frame\n#' @param feature name of the feature\n#' @param binsPerRange a histogram binning parameter\n#' @return a ggplot\n#' @export\nplotDist <- function(train, test, feature, binsPerRange = 100){\n  require(ggthemes)\n  nQuake <- sum((diff(train$ttf)>0)) + 1\n  train$quake <-  factor(unlist(map2(c(1:nQuake),diff(c(0,which(diff(train$ttf)>0),nrow(train))), ~ rep(.x, each = .y))))\n  binwidth <- diff(range(train[,feature]))/binsPerRange\n\n  ggplot() + geom_freqpoly(data = train,aes_string(x = feature, y = \"..density..\",color = \"quake\"),binwidth = binwidth)  + \n      geom_freqpoly(data = test, aes_string(x = feature,y = \"..density..\"), binwidth = binwidth,color = \"black\", size = 1) + \n      geom_freqpoly(data = train, aes_string(x = feature,y = \"..density..\"), binwidth = binwidth, color = \"red\", size = 1) +\n      scale_colour_tableau(\"Tableau 20\") + theme_minimal()\n}    \n    \n#' A function to run CV strategy on training data after removing certain quakes for validation\n#' @param trainSet training set data frame\n#' @param testSet test set data frame, to predict on the fly and conserve memory\n#' @param validQuakes the quakes to use for validation\n#' @param validationStrategy the validation strategy to be applied. \n#'     Possible values \"shuffledKFold\",\"quakeWise\",\"KFold\", \"repeatedKFold\"\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#'    Set to FALSE for non-shuffled CV\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param objective convenience paramter to pass the objective for training\n#' @param k number of folds\n#' @param seed the random seed     \n#' @param repeats the number of repeats if strategy is repeatedKFold\n#' @return a list: pred is out of fold predictions and testPred the prediction on test set\n#' @export   \nInternalCV <- function(trainSet, testSet, validationStrategy = c(\"shuffledKFold\",\"quakeWise\",\"KFold\", \"repeatedKFold\"),useBest = T,\n                      modelType = c(\"lgbm\",\"catboost\",\"xgboost\"), params, \n                       objective = \"\", seed = 123, k = 5, repeats = 1, validQuakes = c(12:17)){\n     \n    modelType <- match.arg(modelType)\n    validationStrategy <- match.arg(validationStrategy)\n  \n  # count how many quakes in input training data\n    nQuake <- sum((diff(trainSet$ttf)>0)) + 1\n  \n  # get the quake IDs\n   quake <- unlist(map2(c(1:nQuake),diff(c(0,which(diff(trainSet$ttf)>0),nrow(trainSet))), ~ rep(.x, each = .y)))\n  \n   train <- trainSet[!quake %in% validQuakes,]\n   valid <- trainSet[quake %in% validQuakes,]\n   validPredict <- select(valid, -ttf)\n   \n   cv <- performCV(train, validationStrategy = validationStrategy, useBest = useBest, seed=seed,modelType = modelType,\n                      k = k, objective = objective, params = params, repeats = repeats)\n   if(repeats > 1){\n       repeated <- TRUE\n   } else { \n       repeated <- FALSE\n   }\n   p <- CVpredict(cv, validPredict, repeated = repeated, modelType = modelType, useBest = useBest) \n   testPred <- CVpredict(cv, testSet, repeated = repeated, modelType = modelType, useBest = useBest)\n   cat(paste0(\"MAE for Nested CV in quake \",paste(validQuakes , collapse = \" \"), \" is \",mae(p, valid$ttf),\"\\n\"))\n   return(list(pred = p, testPred = testPred)) \n}\n\n#' A function to run nested CV. The inner CV loop can be from any strategy. The outer loop is by quake\n#' @param trainSet training set data frame\n#' @param testSet test set data frame, to predict on the fly and conserve memory\n#' @param validationStrategy the validation strategy to be applied. \n#'     Possible values \"shuffledKFold\",\"quakeWise\",\"KFold\", \"repeatedKFold\"\n#' @param useBest whether to use the best iteration (default) or the last iteration (if FALSE). \n#'    Set to FALSE for non-shuffled CV\n#' @param modelType the model type can be catboost, lgbm or xgboost\n#' @param objective convenience paramter to pass the objective for training\n#' @param k number of folds\n#' @param seed the random seed     \n#' @param repeats the number of repeats if strategy is repeatedKFold\n#' @return a list: pred is out of fold predictions and testPred the prediction on test set\n#' @export       \nNestedCV <- function(trainSet, testSet, validationStrategy = c(\"shuffledKFold\",\"quakeWise\",\"KFold\", \"repeatedKFold\"),useBest = T,\n                      modelType = c(\"lgbm\",\"catboost\",\"xgboost\"), params, \n                       objective = \"\", seed = 123, k = 5, repeats = 1){\n    nQuake <- sum((diff(trainSet$ttf)>0)) + 1\n    models <- map(c(1:nQuake), ~ InternalCV(trainSet, testSet, validationStrategy = validationStrategy, \n                      useBest = useBest, seed=seed, modelType = modelType,\n                      k = k, objective = objective, params = params, repeats = repeats, validQuakes = .x))\n    pred <- unlist(map(models, \"pred\"))\n    testPred <- rowMeans(map_dfc(models, \"testPred\"))\n    cat(paste(\"Final Nested CV is \", mae(pred, trainSet$ttf),\"\\n\"))\n    return(list(pred = pred, testPred = testPred))\n}    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"if(F){\ncat(\"Extracting Features ...\")    \ntic()\n# the batches of data to read, for conserving memory and running in parallel\nskips <- c(1,seq(37500000, 600000000, by = 37500000))\nnmax <- c(rep(37500000,length(skips)-1), 29100000)\n\nplan(multiprocess)\n# loop over skips and nmax and run processData\n# bind outputs by row\n# remove ttf == 0 segments where resets happen\ntrain <- future_map2_dfr(skips, nmax, processData) %>% \n          as.data.frame() %>%\n          filter(ttf > 0)\n\n# remove features with near zero variance\nnzv <- nearZeroVar(train)\nif(length(nzv)>0){\n   train <- train[,-nzv] \n}\nwrite.csv(train, \"train.csv\", row.names = F, quote = F)\n\n# extract features from test segments\nsampleSub <- read.csv(\"../input/LANL-Earthquake-Prediction/sample_submission.csv\", stringsAsFactors=F)\ntestFeatures <- future_map_dfr(sampleSub$seg_id, buildTestSet)\ntestFeatures <- testFeatures[,colnames(testFeatures) %in% colnames(train)]\nwrite.csv(testFeatures, file = \"test.csv\", row.names = F, quote = F)\ntoc()\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"if(F){\nlgbmParams <- list(metric=\"mae\",\n                   verbosity = -1,\n                   max_depth = -1,\n                   seed = 123,\n                   feature_fraction = 0.2,\n                   num_leaves = 41,\n                   min_data_in_bin = 25,\n                   max_bin = 350,\n                   min_gain_to_split = 0.2,\n                   num_iterations = 500000,\n                   early_stopping_round = 100,\n                   #lambda_l2 = 0.01,\n                   #lambda_l1 = 0.01,\n                   learning_rate=0.01)\n\nmodel <- performCV(train, validationStrategy = \"shuffledKFold\", useBest = T,seed=1234,\n                      k = 10, objective = \"gamma\", params = lgbmParams)\n\nsampleSub <- read.csv(\"../input/LANL-Earthquake-Prediction/sample_submission.csv\", stringsAsFactors=F)\n# submission is average of 10 bags\nsub_lgb <- sampleSub %>% \n       mutate(time_to_failure = CVpredict(model, testFeatures))\nwrite.csv(sub_lgb, file = \"sub_lgb.csv\", row.names = F, quote = F)\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"sessionInfo()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"codemirror_mode":"r","file_extension":".r","mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.6.0"}},"nbformat":4,"nbformat_minor":1}