{"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":"code","source":"suppressMessages({\n    library(arrow)\n    library(dplyr)\n    library(ggplot2)\n})","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-09-23T16:38:40.024850Z","iopub.execute_input":"2023-09-23T16:38:40.027238Z","iopub.status.idle":"2023-09-23T16:38:41.142810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Read the data\n\nWhen the data set was split into 4 groups all four splits somehow preserved all factor levels in each file.\nDo not need empty series. Removing empty factor levels from the data set.","metadata":{}},{"cell_type":"code","source":"dtr1 <- read_parquet(\"../input/sleep-states-train-set-4-pieces-split/train_series_1.parquet\")\ndtr1$series_id <- as.character(dtr1$series_id)\ndtr1$timestamp <- as.POSIXct(dtr1$timestamp, format=\"%Y-%m-%dT%H:%M:%S%z\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:38:41.145847Z","iopub.execute_input":"2023-09-23T16:38:41.189983Z","iopub.status.idle":"2023-09-23T16:39:54.066452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"te <- read.csv(\"../input/child-mind-institute-detect-sleep-states/train_events.csv\")\nte$event <- as.factor(te$event)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:39:54.080469Z","iopub.execute_input":"2023-09-23T16:39:54.082542Z","iopub.status.idle":"2023-09-23T16:39:54.141699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the `train_events.csv` data set there are series with events but no steps or timestamps. \nNeed to ignore those.","metadata":{}},{"cell_type":"code","source":"no_steps_ser <- te %>% \n  group_by(series_id) %>% summarise(no_steps=all(is.na(step) & timestamp==\"\")) %>% filter(no_steps) %>% pull(series_id)\nte <- te %>% filter( !(series_id %in% no_steps_ser) ) \ndtr1 <- dtr1 %>% filter( !(series_id %in% no_steps_ser) ) ","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:39:54.145943Z","iopub.execute_input":"2023-09-23T16:39:54.148257Z","iopub.status.idle":"2023-09-23T16:39:57.865636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Enhance the time series with event annotations","metadata":{}},{"cell_type":"code","source":"dtr1e <- left_join(dtr1, te[,c(\"series_id\",\"step\",\"night\",\"event\")])\nsummary(dtr1e)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:39:57.868316Z","iopub.execute_input":"2023-09-23T16:39:57.869777Z","iopub.status.idle":"2023-09-23T16:40:14.104136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Fill in the state gaps with asleep and awake statuses. The code to do that is writtent in C++. Loops written in R usually do not perform well.\nPeriods of status change longer then 24 hours are not getting classified. ","metadata":{}},{"cell_type":"code","source":"library(Rcpp)\n\nsrc <- '#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 ) {\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-09-23T16:40:36.760806Z","iopub.execute_input":"2023-09-23T16:40:36.762443Z","iopub.status.idle":"2023-09-23T16:40:44.958101Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr1e <- dtr1e %>% group_by(series_id) %>% arrange(step) %>% \n  mutate(status=fill_event_gaps(event, as.integer(timestamp),24))","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:40:44.962399Z","iopub.execute_input":"2023-09-23T16:40:44.964510Z","iopub.status.idle":"2023-09-23T16:40:52.941560Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Fill in the night number gaps in a similar way to status gaps","metadata":{}},{"cell_type":"code","source":"src <- '#include <Rcpp.h>\nusing namespace Rcpp;\n// [[Rcpp::export]]\nIntegerVector fill_night_gaps(IntegerVector inp){\n  IntegerVector outp(inp.size());\n  int current = 1;\n  for(int i=0; i<inp.length(); ++i){\n    if( inp[i] > current ) {\n      current = inp[i];\n    }\n    outp[i] = current;\n  }\n  return(outp);\n}'\n\nsourceCpp(code = src)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:40:59.146904Z","iopub.execute_input":"2023-09-23T16:40:59.148505Z","iopub.status.idle":"2023-09-23T16:41:04.566425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr1e <- dtr1e %>% group_by(series_id) %>% arrange(step) %>% \n  mutate(filled_night=fill_night_gaps(night))","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:41:04.569237Z","iopub.execute_input":"2023-09-23T16:41:04.570785Z","iopub.status.idle":"2023-09-23T16:41:10.937192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(dtr1e)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:41:10.939960Z","iopub.execute_input":"2023-09-23T16:41:10.941418Z","iopub.status.idle":"2023-09-23T16:41:23.518341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualisations\n\nSelect a few random series to look at","metadata":{}},{"cell_type":"code","source":"sample_ser <- sample(unique(dtr1e$series_id),4)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:41:39.697154Z","iopub.execute_input":"2023-09-23T16:41:39.702470Z","iopub.status.idle":"2023-09-23T16:41:40.408592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Look at the full timeline. Subsampling the data to make plots run faster.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser) %>% sample_n(100000) %>% \n  ggplot() + geom_point(aes(timestamp,anglez,col=status)) +\n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:41:43.523107Z","iopub.execute_input":"2023-09-23T16:41:43.524954Z","iopub.status.idle":"2023-09-23T16:42:10.660898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser) %>% \n  sample_n(100000) %>% ggplot() + geom_point(aes(timestamp,enmo,col=status)) + \n  scale_y_sqrt() +\n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:42:10.663708Z","iopub.execute_input":"2023-09-23T16:42:10.665214Z","iopub.status.idle":"2023-09-23T16:42:36.444247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Look at the first few days of each series","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser, step<100000) %>% \n  ggplot() + geom_point(aes(timestamp,anglez,col=status)) + \n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:42:36.448959Z","iopub.execute_input":"2023-09-23T16:42:36.451578Z","iopub.status.idle":"2023-09-23T16:43:00.973500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser, step<100000) %>% \n  ggplot() + geom_point(aes(timestamp,enmo,col=status)) + \n  scale_y_sqrt()  + facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:43:00.976730Z","iopub.execute_input":"2023-09-23T16:43:00.978417Z","iopub.status.idle":"2023-09-23T16:43:25.413437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Use time series tools","metadata":{}},{"cell_type":"code","source":"suppressMessages({\n   library(tsibble)\n   library(feasts)\n   library(fable)    \n})","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:45:52.679369Z","iopub.execute_input":"2023-09-23T16:45:52.681149Z","iopub.status.idle":"2023-09-23T16:45:53.450267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr_ts <- tsibble(dtr1e, key=series_id, index=step)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:48:16.055524Z","iopub.execute_input":"2023-09-23T16:48:16.057455Z","iopub.status.idle":"2023-09-23T16:48:25.425477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr_ts %>% filter(series_id==sample_ser[[1]]) %>% head(100000) %>%\n  gg_lag(anglez, geom = \"point\", lags=1:9)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:48:25.428264Z","iopub.execute_input":"2023-09-23T16:48:25.429753Z","iopub.status.idle":"2023-09-23T16:49:13.970469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr_ts %>% filter(series_id %in% sample_ser) %>% \n  ACF(anglez, lag_max=100) %>% autoplot()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:49:13.973235Z","iopub.execute_input":"2023-09-23T16:49:13.974744Z","iopub.status.idle":"2023-09-23T16:49:17.049076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr_ts %>% filter(series_id %in% sample_ser) %>% \n  ACF(enmo, lag_max=100) %>% autoplot()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:49:17.051794Z","iopub.execute_input":"2023-09-23T16:49:17.053236Z","iopub.status.idle":"2023-09-23T16:49:20.117643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There seems to be nothing to learn from overall series autocorrelation.\nTake a look at the periods of sleep and awake for a few observed nights.","metadata":{}},{"cell_type":"code","source":"dtr_ts2 <- dtr1e %>% filter(series_id %in% sample_ser, filled_night<10) %>% \n  tsibble(key=c(series_id, status, filled_night), index=timestamp)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:49:20.122385Z","iopub.execute_input":"2023-09-23T16:49:20.124967Z","iopub.status.idle":"2023-09-23T16:49:22.318336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=10)\nfill_gaps(dtr_ts2, anglez=NA) %>% \n  filter(filled_night<3,status %in% c(\"awake\",\"asleep\")) %>% \n  ACF(anglez, lag_max=100) %>% autoplot()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:49:22.321373Z","iopub.execute_input":"2023-09-23T16:49:22.323063Z","iopub.status.idle":"2023-09-23T16:49:27.491233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=10)\nfill_gaps(dtr_ts2, enmo=NA) %>% \n  filter(filled_night<3,status %in% c(\"awake\",\"asleep\")) %>% \n  ACF(enmo, lag_max=100) %>% autoplot()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:49:27.494099Z","iopub.execute_input":"2023-09-23T16:49:27.495806Z","iopub.status.idle":"2023-09-23T16:49:32.418098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Difference time series","metadata":{}},{"cell_type":"code","source":"dtr1e <- dtr1e %>% arrange(series_id, step) %>% \n  mutate(enmo_diff = difference(enmo), anglez_diff = difference(anglez))\ndtr_ts <- dtr1e %>% filter(series_id %in% sample_ser, filled_night<10) %>% \n  tsibble(key=c(series_id, status, filled_night), index=timestamp)","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:45:58.549757Z","iopub.execute_input":"2023-09-23T16:45:58.551935Z","iopub.status.idle":"2023-09-23T16:46:11.233125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"enmo_diff overall timeline","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser) %>% sample_n(100000) %>% \n  ggplot() + geom_point(aes(timestamp,enmo_diff,col=status)) + \n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:46:11.235967Z","iopub.execute_input":"2023-09-23T16:46:11.237397Z","iopub.status.idle":"2023-09-23T16:46:35.517693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"enmo_diff first few days","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser, step<100000) %>% \n  ggplot() + geom_point(aes(timestamp,enmo_diff,col=status)) + \n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:46:35.520794Z","iopub.execute_input":"2023-09-23T16:46:35.522473Z","iopub.status.idle":"2023-09-23T16:46:58.538522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"anglez_diff overall timeline","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser) %>% sample_n(100000) %>% \n  ggplot() + geom_point(aes(timestamp,anglez_diff,col=status)) + \n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:46:58.541448Z","iopub.execute_input":"2023-09-23T16:46:58.542985Z","iopub.status.idle":"2023-09-23T16:47:23.226267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"anglez_diff first few days","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\ndtr1e %>% filter(series_id %in% sample_ser, step<100000) %>% \n  ggplot() + geom_point(aes(timestamp,anglez_diff,col=status)) + \n  facet_wrap(~series_id, scale=\"free\")","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:47:23.229084Z","iopub.execute_input":"2023-09-23T16:47:23.230591Z","iopub.status.idle":"2023-09-23T16:47:46.588452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Autocorrelation plots","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\nfill_gaps(dtr_ts, anglez_diff=NA) %>% \n  filter(filled_night<3,status %in% c(\"awake\",\"asleep\")) %>% \n  ACF(anglez_diff, lag_max=20) %>% autoplot()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:47:46.592956Z","iopub.execute_input":"2023-09-23T16:47:46.595202Z","iopub.status.idle":"2023-09-23T16:47:51.439931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=8)\nfill_gaps(dtr_ts, enmo_diff=NA) %>% \n  filter(filled_night<3,status %in% c(\"awake\",\"asleep\")) %>% \n  ACF(enmo_diff, lag_max=20) %>% autoplot()","metadata":{"execution":{"iopub.status.busy":"2023-09-23T16:47:51.442815Z","iopub.execute_input":"2023-09-23T16:47:51.444366Z","iopub.status.idle":"2023-09-23T16:47:56.513312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}