Extend gap handling alignment capabilities

[ThomasKraft]

The gaps.R script and gapless_Datetimes() function create a sequence that is anchored to the minimum and maximum of the specified Datetime column. Due to minor time inaccuracies, I'm encountering frequent issues where there is a phase shift in time data sequence that occurs part way through files, e.g., 12:00:02 to 12:00:04 to 12:00:06 to 12:00:07 to 12:00:09 to 12:00:11. This results in a plot like the one below where gaps and irregular data are consistently found even during periods where data collect is actually regular:

I'd like to propose 2 alternative ideas for modifying the behavior of gap handler to handle this type of case. The first would be a procedure that does some basic phase shift detection based on some tolerance paramter to determine how far a timestamp can be from the expected sequence. Then if a phase shift is detected, align the gapless sequence based on actual timestamps instead of forcing alignment to the min/max.

Alternatively, one could simply round timestamps to the nearest epoch interval, so in the example above you would find where the phase shift starts and then round up to 12:00:02, 12:00:04, 12:00:06, 12:00:08, 12:00:10, etc.

[JZauner]

Indeed, this is an issue I have no good answer to. Many functions in LightLogR require a regular sequence in order to work properly. While it could certainly be changed to be based on a by-timepoint-difference basis, I find this unsatisfactory as well.

The phase-shifting detection you propose sounds interesting, but would it work properly if there are many such shifts? You could test it on the included sample.data.irregular dataset, that has a basic 15-second interval, but every 10th or so measurement is 16 seconds (and every 100th or so 17 seconds...) and we could see how the proposed modification would handle that. Would that be an option?

As an aside - indeed, currently if I encounter such irregular data, I either switch to a coarser interval with aggregate_Datetime(), or I round the timestamps with cut_Datetime(Datetime.colname = Datetime). In the latter case it is important to check for duplicate timestamps as a result of rounding.

[ThomasKraft]

I have now experimented with a variety of algorithms, and unfortunately I can come up with cases that make each of them fail. At the end of the day the aggregation procedure that you suggested. And I really appreciate the customization possible in aggregate_Datetime() function, with sensible defaults!

One final note: In your tutorial for VEETs (https://tscnlab.github.io/ZaunerEtAl_JVis_2026/supplement.html#veet-data-spectral-data-processing) you aggregate to 5 second intervals (there's a small typo, in the text it says 10 sec but in the code it is 5), and by default aggregate_Datetime() will take the mean. Later you replace 0 distances with Infinity. Taking the mean with the zero values still in is probably not the best strategy.

I don't have a perfect alternative, but here is what I have developed: For each 8×8 grid distance cell over the specified window, I aggregate two separate things.

1. occupancy / detection proportion
2. conditional distance summary statistics when detected (proportion of detected, mean/median/min distance of those non-zero, etc.)

This gives the user a good range of faithful metrics that can be used in zero-inflated/hurdle models or to capture average/median distances in cases specifically across observations that are not zero (therefore avoiding the zero averaging problem). I append those metrics to each grid cell in a standard format to make it easy to access.

I'm attaching a modified version of your function here if you would like to review. Note that this allows the user to specify a parameter in the function call (modality_ToF) that when = TRUE, tells the function to use my procedures on columns starting with a "dist" string while leaving others to the default procedures. I can PR this from my fork if desired, but I realize that you may already have implemented something like this and I am just not seeing it! Please accept my apologies if that is the case.

#' Aggregate Datetime data
#'
#' Condenses a `dataset` by aggregating the data to a given (shorter) interval
#' `unit`. [aggregate_Datetime()] is opinionated in the sense that it sets
#' default handlers for each data type of `numeric`, `character`, `logical`,
#' `factor`, `duration`, `time`, and `datetime`. These can be overwritten by the
#' user. Columns that do not fall into one of these categories need to be
#' handled individually by the user (`...` argument) or will be removed during
#' aggregation. If no unit is specified the data will simply be aggregated to
#' the most common interval (`dominant.epoch`), which is most often not an
#' aggregation but a rounding.)
#'
#' Summary values for type `POSIXct` are calculated as the mean, which can be
#' nonsensical at times (e.g., the mean of Day1 18:00 and Day2 18:00, is Day2
#' 6:00, which can be the desired result, but if the focus is on time, rather
#' then on datetime, it is recommended that values are converted to times via
#' [hms::as_hms()] before applying the function (the mean of 18:00 and 18:00 is
#' still 18:00, not 6:00).
#'
#' @inheritParams cut_Datetime
#' @param
#' numeric.handler,character.handler,logical.handler,factor.handler,datetime.handler,duration.handler,time.handler
#' functions that handle the respective data types. The default handlers
#' calculate the `mean` or `median` for `numeric`, `POSIXct`, `duration`, and
#' `hms`, and the `mode` for `character`, `factor` and `logical` types.
#' @param unit Unit of binning. See [lubridate::round_date()] for examples. The
#' default is `"dominant.epoch"`, which means everything will be aggregated to
#' the most common interval. This is especially useful for slightly irregular
#' data, but can be computationally expensive. `"none"` will not aggregate the
#' data at all.
#' @param modality_ToF Logical. If `FALSE`, numeric columns are aggregated with
#' `numeric.handler`, as before. If `TRUE`, numeric columns whose names begin
#' with `"dist"` are treated as time-of-flight distance cells in which zero
#' means no valid return, not a true distance. Each `dist*` column is expanded
#' into detection and conditional distance summaries: `p_detected`, `any_detected`, `dist_mean`, `dist_median`,
#' `dist_min`, and `dist_max`. Numeric columns that do not begin
#' with `"dist"` are aggregated with `numeric.handler`.
#' @param ... arguments given over to [dplyr::summarize()] to handle columns
#' that do not fall into one of the categories above. Be careful with partial
#' matching of argument names. E.g., creating a new column `n`, without having
#' changed the default for `numeric.handler` will match `n` to that argument.
#' You can avoid this be either explicitly supplying the `numeric.handler`
#' argument (using the default again), or you change the variable name, e.g,
#' to `n.`, which avoids the partial matching.
#'
#' @return A `tibble` with aggregated `Datetime` data. Usually the number of
#' rows will be smaller than the input `dataset`. If `modality_ToF = FALSE` and
#' the handler arguments capture all column types, the number of columns will be
#' the same as in the input `dataset`. If `modality_ToF = TRUE`, each numeric
#' ToF column whose name begins with `"dist"` is expanded into multiple summary
#' columns, while other numeric columns retain the standard numeric aggregation.
#' @export
#'
#' @examples
#' # dominant epoch without aggregation
#' sample.data.environment %>%
#'   dominant_epoch()
#'
#' # dominant epoch with 5 minute aggregation
#' sample.data.environment %>%
#'   aggregate_Datetime(unit = "5 mins") %>%
#'   dominant_epoch()
#'
#' # ToF aggregation from 2 second data to 30 second bins
#' # tof_data %>%
#' #   aggregate_Datetime(unit = "30 secs", modality_ToF = TRUE)
#'
aggregate_Datetime <- function(dataset,
                               unit = "dominant.epoch",
                               Datetime.colname = Datetime,
                               type = c("round", "floor", "ceiling"),
                               numeric.handler = mean,
                               character.handler = \(x) names(which.max(table(x, useNA = "ifany"))),
                               logical.handler = \(x) mean(x) >= 0.5,
                               factor.handler = \(x) factor(names(which.max(table(x, useNA = "ifany")))),
                               datetime.handler = mean,
                               duration.handler = \(x) lubridate::duration(mean(x)),
                               time.handler = \(x) hms::as_hms(mean(x)),
                               modality_ToF = FALSE,
                               ...) {

  # Initial Checks ----------------------------------------------------------
  # Match input arguments
  type <- match.arg(type)

  Datetime.colname.str <- colname.defused({{ Datetime.colname }})
  Datetime.colname.defused <- colname.defused({{ Datetime.colname }}, as_string = FALSE)

  stopifnot(
    "dataset is not a dataframe" = is.data.frame(dataset),
    "Datetime.colname must be part of the dataset" =
      Datetime.colname.str %in% names(dataset),
    "Datetime.colname must be a Datetime" =
      lubridate::is.POSIXct(dataset[[Datetime.colname.str]]),
    "numeric.handler must be a function" = is.function(numeric.handler),
    "character.handler must be a function" = is.function(character.handler),
    "logical.handler must be a function" = is.function(logical.handler),
    "factor.handler must be a function" = is.function(factor.handler),
    "datetime.handler must be a function" = is.function(datetime.handler),
    "modality_ToF must be TRUE or FALSE" =
      is.logical(modality_ToF) && length(modality_ToF) == 1 && !is.na(modality_ToF)
  )

  # ToF helpers -------------------------------------------------------------
  # VEET/ToF zero values encode no valid return/out-of-range, not distance 0.
  # When modality_ToF = TRUE, only numeric columns with names beginning with
  # "dist" are treated this way. Other numeric columns retain numeric.handler.
  # Therefore, each selected distance cell is summarized as a two-part process:
  #   1. detection / occupancy in the aggregation window
  #   2. distance conditional on a valid positive return
  ToF_detected <- function(x) {
    !is.na(x) & x > 0
  }

  ToF_positive <- function(x) {
    x[ToF_detected(x)]
  }

  ToF_empty_na <- function(x, fun) {
    x <- ToF_positive(x)
    if (length(x) == 0) {
      NA_real_
    } else {
      fun(x)
    }
  }

  ToF_iqr <- function(x) {
    x <- ToF_positive(x)
    if (length(x) < 2) {
      NA_real_
    } else {
      stats::IQR(x)
    }
  }

  ToF.handlers <- list(
    p_detected = \(x) mean(ToF_detected(x)),
    any_detected = \(x) any(ToF_detected(x)),
    dist_mean = \(x) ToF_empty_na(x, mean),
    dist_median = \(x) ToF_empty_na(x, stats::median),
    dist_min = \(x) ToF_empty_na(x, min),
    dist_max = \(x) ToF_empty_na(x, max)
  )

  # Function ----------------------------------------------------------
  unit2 <- unit |> as.character()

  if (unit2 != "none") {
    dataset <-
      dataset %>%
      cut_Datetime(
        Datetime.colname = !!Datetime.colname.defused,
        unit = unit,
        type = type,
        group_by = TRUE
      ) %>%
      dplyr::summarize(
        ..., # allow for additional functions

        # If modality_ToF is FALSE, preserve original numeric behavior.
        # If TRUE, only columns starting with "dist" are expanded into ToF
        # detection and conditional-distance summaries. Other numeric columns
        # are still aggregated with numeric.handler.
        dplyr::across(
          dplyr::where(is.numeric) & !dplyr::starts_with("dist"),
          numeric.handler,
          .names = "{.col}"
        ),

        dplyr::across(
          dplyr::where(is.numeric) & dplyr::starts_with("dist"),
          if (modality_ToF) ToF.handlers else numeric.handler,
          .names = if (modality_ToF) "{.col}_{.fn}" else "{.col}"
        ),

        dplyr::across(dplyr::where(is.character), !!character.handler),
        dplyr::across(dplyr::where(is.logical), !!logical.handler),
        dplyr::across(dplyr::where(is.factor), !!factor.handler),
        dplyr::across(dplyr::where(lubridate::is.POSIXct), !!datetime.handler),
        dplyr::across(dplyr::where(lubridate::is.duration), !!duration.handler),
        dplyr::across(dplyr::where(hms::is_hms), !!time.handler),
        .groups = "drop_last"
      ) %>%
      dplyr::select(-{{ Datetime.colname }}) |>
      dplyr::rename(Datetime = Datetime.rounded)
  }

  # Return ----------------------------------------------------------
  dataset
}