{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Working with the whole trainig set on my 16G laptop is difficult. I am breaking the whole set by session into a few (4 seems good enough) pieces to work with the data.\n\nThere is some feature engineering applied that I found useful so far.","metadata":{}},{"cell_type":"code","source":"suppressMessages({\n    library(arrow)\n    library(dplyr)\n    library(tsibble)\n    library(Rcpp)\n    library(slider)\n})","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-10-05T22:41:56.304678Z","iopub.execute_input":"2023-10-05T22:41:56.363070Z","iopub.status.idle":"2023-10-05T22:41:56.388372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr <- read_parquet(\"../input/child-mind-institute-detect-sleep-states/train_series.parquet\")","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:42:09.353936Z","iopub.execute_input":"2023-10-05T22:42:09.355448Z","iopub.status.idle":"2023-10-05T22:42:36.938677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr$series_id <- as.factor(dtr$series_id)\nids <- levels(dtr$series_id)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:42:36.942397Z","iopub.execute_input":"2023-10-05T22:42:36.944118Z","iopub.status.idle":"2023-10-05T22:43:41.062190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(dtr)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:43:41.064771Z","iopub.execute_input":"2023-10-05T22:43:41.066274Z","iopub.status.idle":"2023-10-05T22:44:16.024472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Helper function in C++ (for better performance) to calculate date hour. Make use of OMP as well, calculations can be done in parallel.","metadata":{}},{"cell_type":"code","source":"src <- '#include <Rcpp.h>\nusing namespace Rcpp;\n// [[Rcpp::plugins(openmp)]]\n// [[Rcpp::export]]\nNumericVector calc_hour(const std::vector< std::string > &inp){\n  std::vector<double> outp(inp.size());\n  #pragma omp parallel\n  for(std::size_t i=0; i < inp.size(); ++i) {\n    const char *p = inp[i].c_str();\n    double h = strtod(p+11, NULL);\n    outp[i] = h + strtod(p+14, NULL)/60.;\n  }\n  return Rcpp::wrap(outp);\n}'\n\nsourceCpp(code = src)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:44:16.026898Z","iopub.execute_input":"2023-10-05T22:44:16.028114Z","iopub.status.idle":"2023-10-05T22:44:23.923699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Helper function (in C++ for better performance) to evaluate sleep/awake state","metadata":{}},{"cell_type":"code","source":"src <- '#include <Rcpp.h>\nusing namespace Rcpp;\n// [[Rcpp::export]]\nIntegerVector fill_event_gaps(IntegerVector inp, const IntegerVector ts, int max_hrs) {\n  IntegerVector outp(inp.size());\n  const int ionset = 1;\n  const int iwakeup = 2;\n  const int iasleep = 3; \n  const int iawake = 4;\n  int current = NA_INTEGER;\n  int initial = NA_INTEGER;\n  int initial_index = -1;\n  int last_index = 0;\n  for(int i=0; i<inp.length(); i++){\n    if( inp[i] == NA_INTEGER ) {\n      outp[i] = current;\n    } else {\n      if( initial_index < 0 ) {\n        initial_index = i;\n        initial = inp[i]==ionset?iawake:iasleep;\n      }\n      current = inp[i]==iwakeup?iawake:iasleep;\n      if( (ts[i] - ts[last_index]) > max_hrs*3600 ) {\n        for(int j=last_index+1; j<i; j++){\n          outp[j] = NA_INTEGER;\n        }\n      }\n      last_index = i;\n      outp[i] = inp[i];\n    }\n  }\n  if( initial_index > 0 && (ts[initial_index] - ts[0]) < max_hrs*3600 ) {\n    for(int i=0; i<initial_index; ++i){\n      outp[i] = initial;\n    }\n  }\n  if( last_index != 0 && (ts[inp.size()-1] - ts[last_index]) > max_hrs*3600 ) {\n    for(int j=last_index+1; j<inp.size(); j++){\n      outp[j] = NA_INTEGER;\n    }\n  }\n  CharacterVector ch = {\"onset\",\"wakeup\",\"asleep\",\"awake\"};\n  outp.attr(\"class\") = \"factor\";\n  outp.attr(\"levels\") = ch;\n  return(outp);\n}'\n\nsourceCpp(code = src)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:44:23.927996Z","iopub.execute_input":"2023-10-05T22:44:23.929361Z","iopub.status.idle":"2023-10-05T22:44:29.133336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Calculate standard deviation over a sliding window","metadata":{}},{"cell_type":"code","source":"src <- '#include <Rcpp.h>\nusing namespace Rcpp;\n// [[Rcpp::export]]\nNumericVector win_sd(const NumericVector inp, const int half_window) {\n  NumericVector outp(inp.size());\n  std::vector<double> sqrs(inp.size());\n  for(std::size_t i=0; i< half_window; ++i) {\n    outp[i] = NA_REAL;\n    outp[inp.size()-i-1] = NA_REAL;\n  }\n  for(std::size_t i=0; i < inp.size(); ++i) {\n    sqrs[i] = inp[i]*inp[i];\n  }\n  const int N = (half_window*2 + 1);\n  const double nn = sqrt((double)(N*(N-1)));\n  double s2 = 0, s1 = 0;\n  for(std::size_t i=0; i <= 2*half_window; ++i) {\n    s2 = s2 + sqrs[i];\n    s1 = s1 + inp[i];  \n  }\n  for(std::size_t i=half_window; i < inp.size()-half_window; ++i) {\n    double v = N*s2 - s1*s1;\n    if(v <= 0) {\n      outp[i] = 0.;\n    } else {\n      outp[i] = sqrt(v)/nn;\n    }\n    s2 = s2 - sqrs[i-half_window] + sqrs[i+half_window+1];\n    s1 = s1 - inp[i-half_window] + inp[i+half_window+1];\n  }\n  return outp;\n}'\n\nsourceCpp(code = src)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:44:29.137326Z","iopub.execute_input":"2023-10-05T22:44:29.138864Z","iopub.status.idle":"2023-10-05T22:44:34.359744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Read recorded events data","metadata":{}},{"cell_type":"code","source":"te <- read.csv(\"../input/child-mind-institute-detect-sleep-states/train_events.csv\")\nte$series_id <- as.factor(te$series_id)\nte$event <- as.factor(te$event)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:44:34.362825Z","iopub.execute_input":"2023-10-05T22:44:34.364206Z","iopub.status.idle":"2023-10-05T22:44:34.432980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Drop series without any recorded events","metadata":{}},{"cell_type":"code","source":"no_steps_ser <- te %>% group_by(series_id) %>% summarise(no_steps=all(is.na(step) & timestamp==\"\")) %>% \n  filter(no_steps) %>% pull(series_id)\nte <- te %>% filter( !(series_id %in% no_steps_ser) ) %>% mutate(series_id=droplevels(series_id))","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:44:34.435594Z","iopub.execute_input":"2023-10-05T22:44:34.437037Z","iopub.status.idle":"2023-10-05T22:44:34.491432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr <- dtr %>% filter( !(series_id %in% no_steps_ser) )\ndtr$series_id <- droplevels(dtr$series_id)\ngc(verbose=FALSE)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:44:34.493862Z","iopub.execute_input":"2023-10-05T22:44:34.495147Z","iopub.status.idle":"2023-10-05T22:46:16.567248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Add and transform features. Work on one split at a time, otherwise I keep running out of memory.","metadata":{}},{"cell_type":"code","source":"add_transform_feature <- function(dtr) {\n  dtr$hour <- calc_hour(dtr$timestamp)\n  dtr$timestamp <- as.POSIXct(dtr$timestamp, format=\"%Y-%m-%dT%H:%M:%S%z\")\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(status=fill_event_gaps(event, as.integer(timestamp), 24))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(enmo_diff = difference(enmo), anglez_diff = difference(anglez))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(enmo_sd1m = win_sd(enmo, 6))    \n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(enmo_sd10m = win_sd(enmo, 6*10))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(anglez_sd1m = win_sd(anglez, 6))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(anglez_sd10m = win_sd(anglez, 6*10))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(enmo_sd2m = win_sd(enmo, 6*2))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(anglez_sd2m = win_sd(anglez, 6*2))   \n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(enmo_sqrt_sd1m = win_sd(sqrt(enmo), 6))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(anglez_sqrt_sd1m = win_sd(sqrt(abs(anglez)), 6))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(anglez_sqrt_sd10m = win_sd(sqrt(abs(anglez)), 6*10))\n  dtr <- dtr %>% group_by(series_id) %>% arrange(step) %>% \n    mutate(enmo_sqrt_sd10m = win_sd(sqrt(enmo), 6*10))\n  dtr\n}","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:47:48.617323Z","iopub.execute_input":"2023-10-05T22:47:48.618847Z","iopub.status.idle":"2023-10-05T22:47:48.633157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr <- left_join(dtr, te[,c(\"series_id\",\"step\",\"night\",\"event\")])\ngc(verbose=FALSE)","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:46:16.612327Z","iopub.execute_input":"2023-10-05T22:46:16.615536Z","iopub.status.idle":"2023-10-05T22:46:55.690605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"parts = 4\nids <- levels(dtr$series_id)\npart <- sample(1:parts,length(ids),replace=TRUE)\nrow_parts <- part[as.integer(dtr$series_id)]","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:48:25.357517Z","iopub.execute_input":"2023-10-05T22:48:25.359045Z","iopub.status.idle":"2023-10-05T22:48:26.698245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc(verbose=FALSE)\nfor(i in 1:parts) {\n  cat(\"working on piece\", i, \"\\n\")\n  dd <- add_transform_feature(dtr[row_parts==i,])\n  write_parquet(dd,paste0(\"preproc_train_series_\",i,\".parquet\"),compression=\"zstd\")\n  rm(dd)\n  gc(verbose=FALSE)\n}","metadata":{"execution":{"iopub.status.busy":"2023-10-05T22:48:33.374930Z","iopub.execute_input":"2023-10-05T22:48:33.376561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}