R/wrap-robust.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

#' DeCAFS wrapper: changes amid drift and autocorrelated noise
#'
#' Wraps \code{DeCAFS::DeCAFS()} (Romano, Rigaill, Runge and Fearnhead,
#' 2022), which detects abrupt mean changes when the underlying signal also
#' drifts (random-walk fluctuations) and the noise is AR(1)-autocorrelated,
#' the two regimes in which plain change-in-mean methods over-detect. Model
#' parameters are estimated automatically unless supplied.
#'
#' @param x A numeric vector.
#' @param penalty Penalty \eqn{\beta} for adding a changepoint. Defaults to
#'   \code{2 * log(length(x))}. \code{\link{cpt_detect}} resolves its own
#'   \code{"MBIC"} default to a stronger numeric value (on a five-changepoint
#'   series that is 3 changepoints through the dispatcher against 5 here), so
#'   pass \code{penalty} explicitly when the two must agree.
#' @param model_param Optional list of model parameters
#'   (\code{sdEta}, \code{sdNu}, \code{phi}) as accepted by
#'   \code{DeCAFS::DeCAFS()}; when \code{NULL} they are estimated from the
#'   data.
#' @param ... Additional arguments passed to \code{DeCAFS::DeCAFS()}.
#' @return A \code{ggcpt} object. The \code{data} tibble carries the
#'   estimated signal in its \code{fitted} column.
#' @references
#' \insertRef{romano2022decafs}{ggchangepoint}
#' @export
#' @examplesIf requireNamespace("DeCAFS", quietly = TRUE)
#' set.seed(2026)
#' res <- decafs_wrapper(c(rnorm(100), rnorm(100, 5)))
#' res$changepoints
#' @family changepoint engines
decafs_wrapper <- function(x, penalty = NULL, model_param = NULL, ...) {
  need_pkg("DeCAFS")
  reject_renamed_args(list(...), "decafs")
  reject_managed_args(list(...), "decafs", c(
    warningMessage = paste("the wrapper turns the engine's own advisory off",
                           "and raises the conditions R users expect",
                           "instead")))
  validate_data(x)
  data_vec <- as_uni_vector(x, "decafs")

  # "Manual" is reserved for a number the caller actually supplied
  # (penalty_descriptor() keys on is.numeric(penalty)), so reporting this
  # wrapper's own default as "Manual = 11.8" hid a documented consequence:
  # `?cpt_detect`'s @param penalty explains that the dispatcher's default
  # MBIC resolves to a STRONGER value than 2*log(n) -- 19.9 against 11.8 at
  # n = 360 -- and both paths labelled the result the same way, so the
  # object gave no way to tell which default produced it.
  pen_type <- if (is.null(penalty)) "2log(n) [wrapper default]" else "Manual"
  penalty <- resolve_penalty_model(penalty, data_vec)
  if (is.null(penalty)) {
    penalty <- 2 * log(length(data_vec))
  }
  # A vector failed with Rcpp's "Expecting a single value: [extent=2]", and
  # a negative value segmented every observation.
  validate_scalar(penalty, "penalty", min = 0)

  args <- list(data = data_vec, beta = penalty, warningMessage = FALSE, ...)
  if (!is.null(model_param)) args$modelParam <- model_param
  fit <- do.call(DeCAFS::DeCAFS, args)

  ggcpt_build(
    data_vec, as.integer(fit$changepoints),
    method = "decafs",
    change_in = "mean",
    penalty = list(type = pen_type, value = penalty),
    fit = fit,
    call = match.call(),
    fitted = as.numeric(fit$signal)
  )
}

#' Self-normalisation wrapper (SNSeg)
#'
#' Wraps \code{SNSeg::SNSeg_Uni()} (Zhao, Jiang and Shao, 2022):
#' self-normalised segmentation with nested local windows. Self-normalisation
#' avoids estimating the long-run variance, is robust to temporal dependence,
#' and detects changes in general parameters (mean, variance, quantiles,
#' autocorrelation, or bivariate correlation) within one framework.
#'
#' @param x A numeric vector (or a two-column matrix for
#'   \code{parameter = "bivcor"}).
#' @param parameter Which parameter to test for changes: \code{"mean"},
#'   \code{"variance"}, \code{"acf"}, or \code{"bivcor"} (bivariate
#'   correlation). Defaults to \code{"mean"}.
#' @param confidence Confidence level of the self-normalised test, one of
#'   0.9, 0.95, 0.99, 0.995 or 0.999. Defaults to \code{0.9}.
#' @param grid_size Grid size controlling the local-window sweep; when
#'   \code{NULL} the engine's default is used.
#' @param ... Additional arguments passed to \code{SNSeg::SNSeg_Uni()}.
#' @return A \code{ggcpt} object. About 20 observations are needed for the
#'   nested local windows at the default \code{grid_size}; a constant series
#'   returns an empty result rather than an engine error.
#' @references
#' \insertRef{zhao2022snseg}{ggchangepoint}
#' @export
#' @examplesIf requireNamespace("SNSeg", quietly = TRUE)
#' \donttest{
#' set.seed(2026)
#' res <- sn_wrapper(c(rnorm(150), rnorm(150, 3)))
#' res$changepoints
#' }
#' @family changepoint engines
sn_wrapper <- function(x, parameter = c("mean", "variance", "acf", "bivcor"),
                       confidence = 0.9, grid_size = NULL, ...) {
  need_pkg("SNSeg")
  reject_renamed_args(list(...), "sn")
  # Forwarded to the engine, which reported a bad value from deep inside
  # itself -- "missing value where TRUE/FALSE needed", "negative length
  # vectors are not allowed", "NAs in foreign function call" and the like,
  # none of which names the argument. Measured across all 64 wrapper
  # argument slots; these are the ones that needed it.
  validate_scalar(confidence, "confidence", min = 0, max = 1, min_open = TRUE, max_open = TRUE)
  reject_managed_args(list(...), "sn", c(
    plot_SN = paste("SNSeg would draw to the active device; this wrapper",
                    "returns a ggcpt and leaves plotting to",
                    "`autoplot()`")))
  parameter <- cpt_match_arg(parameter)
  # `NA`, a vector or a string failed inside the engine without naming it.
  if (!is.null(grid_size)) {
    validate_scalar(grid_size, "grid_size", min = 0, min_open = TRUE)
  }

  validate_data(x)
  is_mv <- is.matrix(x) || is.data.frame(x)
  if (parameter == "bivcor") {
    if (!is_mv || ncol(as.matrix(x)) != 2) {
      cpt_abort("`parameter = \"bivcor\"` requires a two-column matrix.",
                class = "wrong_dimension")
    }
    X <- as_mv_matrix(x)
    input <- X
    data_vec <- as.numeric(X[, 1])
  } else {
    if (is_mv) {
      cpt_abort("`parameter = \"", parameter, "\"` requires a numeric vector; ",
                "use `parameter = \"bivcor\"` for bivariate input.",
                class = "wrong_dimension")
    }
    input <- as.numeric(x)
    data_vec <- input
  }

  change_lab <- switch(parameter,
    mean = "mean", variance = "var", acf = "acf", bivcor = "correlation"
  )

  # A flat series leaves the self-normalised statistic undefined (the engine
  # fails with "missing value where TRUE/FALSE needed"); it also plainly has
  # no changepoint.
  if (any(constant_cols(as.matrix(input)))) {
    return(ggcpt_build(data_vec, integer(0), method = "sn",
                       change_in = change_lab,
                       penalty = list(type = "confidence", value = confidence),
                       call = match.call(),
                       data_wide = if (parameter == "bivcor") {
                         mv_data_wide(as.matrix(input))
                       }))
  }

  # Below roughly 20 observations there is no room for the nested local
  # windows and the engine fails with "only 0's may be mixed with negative
  # subscripts"; say what is actually wrong.
  fit <- tryCatch(
    SNSeg::SNSeg_Uni(ts = input, paras_to_test = parameter,
                     confidence = confidence, grid_size = grid_size,
                     plot_SN = FALSE, ...),
    error = function(e) {
      if (grepl("only 0's may be mixed with negative subscripts",
                conditionMessage(e), fixed = TRUE)) {
        cpt_abort("`sn` needs a longer series: ", length(data_vec),
                  " observations leave no room for the self-normalisation ",
                  "windows (about 20 are needed at the default `grid_size`).",
                  class = "short_series")
      }
      # A constant RUN, not a constant series: the guard above catches a
      # column that never moves, but a series with a long enough flat
      # stretch anywhere in it puts a zero variance inside one
      # self-normalisation window, and the engine answers with base R's
      # "missing value where TRUE/FALSE needed". Measured on noise with a
      # leading run of zeros, the smallest breaking run tracks the window
      # size: 6 at n = 60, 10 at n = 100, 14 at n = 150, 20 at n = 200 --
      # so this is diagnosed from the engine's failure rather than
      # predicted, which avoids reverse-engineering `grid_size` and cannot
      # refuse a series the engine would have handled.
      if (grepl("missing value where TRUE/FALSE needed",
                conditionMessage(e), fixed = TRUE)) {
        runs <- rle(as.numeric(data_vec))$lengths
        if (max(runs) >= 2L) {
          cpt_abort("`sn` could not self-normalise this series. Its longest ",
                     "run ", "of identical values is ", max(runs), " of ",
                    length(data_vec), " observations, which leaves a window ",
                     "with ", "zero variance; roughly a tenth of the series ",
                     "is enough to ", "do it. Use `parameter = \"variance\"` ",
                     "on a series that does ", "vary, jitter the ties, or ",
                     "pick a method that tolerates flat ",
                    "stretches (`pelt`, `binseg`, `pettitt`).",
                    class = "input_error")
        }
      }
      stop(e)
    }
  )

  ggcpt_build(
    data_vec, as.integer(fit$est_cp),
    method = "sn",
    change_in = change_lab,
    penalty = list(type = "confidence", value = confidence),
    fit = fit,
    call = match.call(),
    data_wide = if (parameter == "bivcor") mv_data_wide(input)
  )
}