R/wrap-highdim.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

#' inspect wrapper: high-dimensional changepoints via sparse projection
#'
#' Wraps \code{InspectChangepoint::inspect()} (Wang and Samworth, 2018). For
#' a \eqn{p}-variate series whose mean changes in an unknown sparse subset of
#' coordinates, the algorithm computes the CUSUM transformation, finds the
#' optimal sparse projection direction via a convex relaxation, and locates
#' changepoints on the projected univariate series, recursing via wild binary
#' segmentation.
#'
#' @param x A numeric matrix or data frame with one row per time point and
#'   one column per coordinate.
#' @param lambda Regularisation parameter of the sparse projection; when
#'   \code{NULL} the engine default \eqn{\sqrt{\log(p \log n)/2}} is used.
#' @param threshold Detection threshold; when \code{NULL} it is computed by
#'   Monte Carlo (via the engine).
#' @param ... Additional arguments passed to
#'   \code{InspectChangepoint::inspect()}.
#' @return A \code{ggcpt} object. The changepoints tibble carries a
#'   \code{strength} column (the maximum projected CUSUM statistic). The
#'   first coordinate is used for \code{cp_value} and the univariate plot
#'   line; the full matrix is kept for the faceted multivariate
#'   \code{autoplot()}. Coordinates that are constant carry no changepoint
#'   information and would make the engine's variance rescaling undefined, so
#'   they are dropped (with a warning) before detection and an all-constant
#'   matrix returns an empty result; the dropped coordinates are still kept
#'   for plotting, and reported locations always refer to the original rows.
#' @references
#' \insertRef{wang2018inspect}{ggchangepoint}
#' @export
#' @examplesIf requireNamespace("InspectChangepoint", quietly = TRUE)
#' set.seed(2026)
#' X <- cbind(c(rnorm(80), rnorm(80, 3)), c(rnorm(80), rnorm(80, -2)),
#'            rnorm(160))
#' res <- inspect_wrapper(X)
#' res$changepoints
#' @family changepoint engines
inspect_wrapper <- function(x, lambda = NULL, threshold = NULL, ...) {
  need_pkg("InspectChangepoint")

  validate_data(x)
  X <- as_mv_matrix(x)
  data_vec <- as.numeric(X[, 1])

  # A flat coordinate makes the engine's variance rescaling divide by zero,
  # which fails with "missing value where TRUE/FALSE needed" even when the
  # other coordinates carry a real change.
  X_fit <- drop_constant_cols(X, "inspect")
  if (is.null(X_fit)) {
    return(ggcpt_build(data_vec, integer(0), method = "inspect",
                       change_in = "mean",
                       penalty = list(type = "threshold", value = NA_real_),
                       call = match.call(), data_wide = mv_data_wide(X)))
  }

  # The engine divides each coordinate by mad(diff(x)) / sqrt(2), which is
  # zero whenever more than half of a coordinate's successive differences
  # are equal: a 0/1 alternation, a noiseless step, a count series with
  # long runs. That failed with base R's "missing value where TRUE/FALSE
  # needed". A flat coordinate is dropped above because it carries nothing;
  # one of these moves, and may carry the change, so it is refused by name
  # rather than dropped.
  no_scale <- vapply(seq_len(ncol(X_fit)), function(j) {
    stats::mad(diff(X_fit[, j])) == 0
  }, logical(1))
  if (any(no_scale)) {
    nm <- colnames(X_fit) %||% paste0("V", seq_len(ncol(X_fit)))
    cpt_abort("`inspect` rescales each coordinate by the median absolute ",
              "deviation of its successive differences, and that is zero for ",
              paste(nm[no_scale], collapse = ", "), ": more than half of the ",
              "differences are identical. Drop or jitter ",
              if (sum(no_scale) > 1) "those coordinates" else "that coordinate",
              ", or use a method that does not rescale this way, such as ",
              "`ecp` or `kcp`.", class = "input_error")
  }
  # Forwarded unchecked, `NA` failed with "missing value where TRUE/FALSE
  # needed", a vector with "the condition has length > 1", a string ran
  # (as a lambda, with the default; as a threshold, finding nothing), and a
  # negative threshold put a changepoint at nearly every observation.
  if (!is.null(lambda)) validate_scalar(lambda, "lambda", min = 0)
  if (!is.null(threshold)) validate_scalar(threshold, "threshold", min = 0)
  args <- list(x = t(X_fit), ...)
  if (!is.null(lambda)) args$lambda <- lambda
  if (!is.null(threshold)) args$threshold <- threshold

  # The engine prints Monte Carlo progress; keep the console clean. That
  # took stdout only, which is not all of it: `inspect()` and
  # `sparse.svd()` both call requireNamespace("RSpectra") without
  # `quietly = TRUE`, and RSpectra is merely *suggested* by
  # InspectChangepoint -- so on a machine holding the engine but not
  # RSpectra (the Windows CI runner, for one) every call writes "Loading
  # required namespace" and "Failed with error: there is no package called
  # 'RSpectra'". The engine handles the absence itself, falling back to
  # base::svd, so it is a loading diagnostic rather than a problem -- but
  # it is not ours to print, and thirteen repetitions of it are what
  # truncated a CI test log down to nothing else.
  #
  # suppressMessages() is not enough and was tried first: only "Loading
  # required namespace" is a condition, while requireNamespace() writes the
  # "Failed with error" line straight to stderr. Capturing the message
  # stream is what silences it. Warning *conditions* still reach the
  # caller and errors still propagate -- both verified -- because only the
  # stream is redirected, not the condition system.
  ignore <- utils::capture.output(
    inner <- utils::capture.output(
      fit <- do.call(InspectChangepoint::inspect, args),
      type = "message"))

  cp_mat <- fit$changepoints
  cp_indices <- if (is.null(cp_mat)) integer(0) else as.integer(cp_mat[, "location"])

  ggcpt_build(
    data_vec, cp_indices,
    method = "inspect",
    change_in = "mean",
    penalty = list(type = "threshold",
                   value = if (!is.null(threshold)) threshold else NA_real_),
    fit = fit,
    call = match.call(),
    extra_cp_cols = if (length(cp_indices) > 0) {
      list(strength = as.numeric(cp_mat[, "max.proj.cusum"]))
    },
    data_wide = mv_data_wide(X)
  )
}

#' ocd wrapper: online high-dimensional changepoint detection
#'
#' Wraps the \code{ocd} package (Chen, Wang and Samworth, 2022): online
#' multiscale detection of a mean change in a high-dimensional stream, with
#' worst-case detection-delay guarantees and per-observation cost independent
#' of history. The detector assumes standardised data with known pre-change
#' mean; this wrapper estimates the baseline mean and standard deviation
#' from an initial training window, then monitors the remainder of the
#' series, resetting after each declaration so multiple changes can be
#' found.
#'
#' @param x A numeric matrix or data frame with one row per time point and at
#'   least two columns. The \code{ocd} detector is inherently
#'   high-dimensional and cannot be constructed for a single coordinate, so
#'   univariate input is rejected; use a univariate engine
#'   (see \code{\link{cpt_methods}()}) for one series.
#' @param train Number of initial observations used to estimate the baseline
#'   mean/sd (not monitored). Defaults to
#'   \code{max(20, floor(0.2 * n))}, capped at \code{n/2}.
#' @param thresh Threshold specification passed to
#'   \code{ocd::ChangepointDetector()}; \code{"MC"} (default) calibrates by
#'   Monte Carlo, which is what makes this the slowest wrapper (see the
#'   timing note below). Supplying the three thresholds directly, as a named
#'   numeric vector \code{c(diag =, off_d =, off_s =)}, skips calibration
#'   altogether.
#' @param patience Target average run length to false alarm. Defaults to
#'   \code{5000}.
#' @param beta Assumed lower bound on the squared Euclidean norm of the mean
#'   change. Defaults to \code{1}.
#' @param mc_reps Monte Carlo repetitions for threshold calibration.
#'   Defaults to \code{100}. The cost is linear in this and grows with the
#'   number of coordinates; see the timing note below.
#' @param ... Additional arguments passed to
#'   \code{ocd::ChangepointDetector()}.
#' @return A \code{ggcpt} object. Because the detector is online, it
#'   declares a change without estimating where it began. \code{cp} is the
#'   declaration time \emph{minus one}: the latest changepoint consistent
#'   with the alarm in this package's convention (the last observation
#'   before the change), and exact when the change is declared on the first
#'   observation after it; a later declaration makes \code{cp} late by the
#'   detection delay. The \code{declared_at} column holds the declaration
#'   time itself. (Before 0.6.0 \code{cp} was the declaration time, one
#'   position late even at zero delay, which the convention test found.)
#'   Compare \code{\link{cpm_wrapper}()}, whose engine estimates the
#'   location separately from its \code{detection_time}.
#' @section How long this takes:
#' Nearly all of the run time is \code{ocd}'s Monte Carlo threshold
#' calibration, which happens before a single observation is read. It is
#' linear in \code{mc_reps} and grows with the number of coordinates. Timed
#' on one Linux x86-64 machine at \code{mc_reps = 5}, construction took
#' about 10 s at \eqn{p = 3}, 22 s at \eqn{p = 10} and 113 s at
#' \eqn{p = 50}; raising \code{mc_reps} scales it linearly, so at
#' \eqn{p = 3} it was 38 s at \code{mc_reps = 20} and 189 s at the default
#' \code{mc_reps = 100}. The practical reading is that the default costs
#' \emph{minutes} rather than seconds even for a handful of coordinates,
#' and better than half an hour at \eqn{p = 50}. Another machine will give
#' different absolute numbers; the linearity in \code{mc_reps} is the part
#' to plan around. Monitoring the observations afterwards is cheap by
#' comparison: 0.37 s for a thousand of them at \eqn{p = 3}. Lower
#' \code{mc_reps} while exploring (the example below uses 2, which
#' measures 3.8 s), or pass \code{thresh} directly to skip
#' calibration entirely, which brings the same fit down to a tenth of a
#' second.
#'
#' @references
#' \insertRef{chen2022ocd}{ggchangepoint}
#' @export
#' @examplesIf requireNamespace("ocd", quietly = TRUE)
#' \donttest{
#' set.seed(2026)
#' X <- rbind(matrix(rnorm(60 * 3), 60), matrix(rnorm(40 * 3, 3), 40))
#' # `mc_reps = 2`, not the default 100 and not the 5 this example used to
#' # pass: the calibration is linear in `mc_reps` and is nearly all of the
#' # cost, so 5 measured 9.7 s here against CRAN's 5 s budget and 2
#' # measures 3.8 s for the same answer. Neither is a calibration you
#' # would trust (see the timing section above).
#' res <- ocd_wrapper(X, mc_reps = 2)
#' res$changepoints
#' }
#' @family changepoint engines
ocd_wrapper <- function(x, train = NULL, thresh = "MC", patience = 5000,
                        beta = 1, mc_reps = 100, ...) {
  need_pkg("ocd")
  reject_renamed_args(list(...), "ocd")
  # 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(patience, "patience", min = 1)
  validate_scalar(mc_reps, "mc_reps", min = 1)
  # `beta = 0` or `NA` failed with "missing value where TRUE/FALSE needed",
  # and a vector ran. ocd's beta is a lower bound on the size of the mean
  # change, so it has to be positive.
  validate_scalar(beta, "beta", min = 0, min_open = TRUE)

  validate_data(x)
  X <- if (is.matrix(x) || is.data.frame(x)) as_mv_matrix(x)
       else matrix(as.numeric(x), ncol = 1)
  n <- nrow(X)
  p <- ncol(X)
  # ocd::ChangepointDetector() fails to construct at dim = 1 under this
  # wrapper's default thresh = "MC", where the Monte Carlo threshold search
  # ends in "'dims' cannot be of length 0": the method projects across
  # coordinates, so it needs at least two. (With an explicit numeric
  # `thresh` the constructor does succeed at dim = 1, which is why the
  # shape is checked here rather than left to the engine -- the failure
  # depends on an argument the caller may never have set.) Say so rather
  # than surfacing the engine's internal error.
  if (p < 2) {
    cpt_abort("Method `ocd` is high-dimensional and needs at least two ",
              "coordinates, but `x` has ", p, ". See cpt_methods() for ",
               "univariate methods.", class = "wrong_dimension")
  }
  data_vec <- as.numeric(X[, 1])

  if (is.null(train)) {
    train <- max(20L, floor(0.2 * n))
  }
  # `NA` and strings reached `min()` below and failed with "missing value
  # where TRUE/FALSE needed", and `c(2, 3)` trained on two observations.
  validate_scalar(train, "train", min = 2)
  train <- min(as.integer(train), floor(n / 2))
  if (train < 2) {
    cpt_abort("`train` must be at least 2 observations.",
              class = "short_series")
  }

  estimate_baseline <- function(rows) {
    m <- colMeans(X[rows, , drop = FALSE])
    s <- apply(X[rows, , drop = FALSE], 2, stats::sd)
    s[s == 0 | is.na(s)] <- 1
    list(mean = m, sd = s)
  }
  base <- estimate_baseline(seq_len(train))

  detector <- ocd::ChangepointDetector(dim = p, method = "ocd",
                                       thresh = thresh, patience = patience,
                                       beta = beta, MC_reps = mc_reps, ...)

  # The engine prints "Changepoint declared at time = ..." on detection.
  declared <- integer(0)
  i <- train + 1
  while (i <= n) {
    z <- (X[i, ] - base$mean) / base$sd
    utils::capture.output(detector <- ocd::getData(detector, z))
    if (!identical(ocd::status(detector), "monitoring")) {
      declared <- c(declared, i)
      # Stop monitoring rather than re-train on a window the baseline cannot
      # be estimated from. `i + 2 > n` only guaranteed TWO remaining
      # observations, and a per-coordinate sd from two points is noise for a
      # high-dimensional stream -- `s[s == 0 | is.na(s)] <- 1` catches only
      # the fully degenerate case, so everything after such a re-training
      # was standardised by a scale estimate that was wrong rather than
      # missing, and spurious later declarations are the expected result.
      # The re-training window now has to clear the same floor `train` was
      # validated against.
      if (n - i < train) break
      detector <- ocd::reset(detector)
      # Re-estimate the baseline from a window after the declaration:
      # keeping the pre-change baseline would re-declare immediately on the
      # shifted regime.
      new_train <- seq(i + 1, min(i + train, n))
      base <- estimate_baseline(new_train)
      i <- new_train[length(new_train)] + 1
    } else {
      i <- i + 1
    }
  }

  # A declaration at observation i says the change is at or before i - 1
  # in the "left" convention; reporting i itself was one late even when the
  # alarm fired on the first post-change observation.
  ggcpt_build(
    data_vec, declared - 1L,
    method = "ocd",
    change_in = "mean",
    penalty = list(type = "patience", value = patience),
    fit = detector,
    call = match.call(),
    extra_cp_cols = if (length(declared) > 0) {
      # When the alarm fired, one after `cp`, because ocd declares a change
      # without estimating where it began -- see @return. Kept as a column
      # so the online engines have one name for "when the alarm fired"
      # whether or not the engine also estimates a location.
      list(declared_at = as.integer(declared))
    },
    data_wide = mv_data_wide(X)
  )
}

#' Geometrically-inspired multivariate changepoint wrapper (geomcp)
#'
#' Wraps \code{changepoint.geo::geomcp()} (Grundy, Killick and Mihaylov,
#' 2020): each multivariate observation is mapped to its distance from, and
#' angle to, a reference point, and univariate PELT is run on the two mapped
#' series. Distance changes capture shifts in magnitude, angle changes
#' capture shifts in orientation/correlation structure.
#'
#' @param x A numeric matrix or data frame with one row per time point.
#' @param penalty Penalty for the univariate PELT runs (a
#'   \pkg{changepoint}-style character penalty). Defaults to \code{"MBIC"}.
#' @param mapping Which mapped series' changepoints to report:
#'   \code{"both"} (union, default), \code{"distance"}, or \code{"angle"}.
#' @param ... Additional arguments passed to
#'   \code{changepoint.geo::geomcp()}.
#' @return A \code{ggcpt} object whose changepoints tibble carries a
#'   \code{mapping} column (\code{"distance"} or \code{"angle"}; a location
#'   found in both is labelled \code{"both"}).
#' @references
#' \insertRef{grundy2020geomcp}{ggchangepoint}
#' @export
#' @examplesIf requireNamespace("changepoint.geo", quietly = TRUE)
#' set.seed(2026)
#' X <- rbind(matrix(rnorm(100 * 4), 100), matrix(rnorm(100 * 4, 2), 100))
#' res <- geomcp_wrapper(X)
#' res$changepoints
#' @family changepoint engines
geomcp_wrapper <- function(x, penalty = "MBIC",
                           mapping = c("both", "distance", "angle"), ...) {
  need_pkg("changepoint.geo")
  mapping <- cpt_match_arg(mapping)

  validate_data(x)
  X <- as_mv_matrix(x)
  # changepoint.geo answers "Univariate changepoint analysis is not
  # supported" on a single column -- clearer than base R, but it names
  # neither the method nor the argument, so a reader with several methods in
  # flight cannot tell which call failed. Match the other eight
  # multivariate-only engines.
  if (ncol(X) < 2) {
    cpt_abort("Method `geomcp` is high-dimensional and needs at least two ",
              "coordinates, but `x` has ", ncol(X), ". See cpt_methods() for ",
               "univariate methods.", class = "wrong_dimension")
  }
  data_vec <- as.numeric(X[, 1])

  fit <- changepoint.geo::geomcp(X, penalty = penalty, ...)

  dist_cp <- as.integer(fit@dist.cpts)
  ang_cp <- as.integer(fit@ang.cpts)

  cp_indices <- switch(mapping,
    both = sort(unique(c(dist_cp, ang_cp))),
    distance = sort(unique(dist_cp)),
    angle = sort(unique(ang_cp))
  )

  map_lab <- vapply(cp_indices, function(cp) {
    in_d <- cp %in% dist_cp
    in_a <- cp %in% ang_cp
    if (in_d && in_a) "both" else if (in_d) "distance" else "angle"
  }, character(1))

  ggcpt_build(
    data_vec, cp_indices,
    method = "geomcp",
    change_in = "distribution",
    penalty = list(type = penalty, value = NA_real_),
    fit = fit,
    call = match.call(),
    extra_cp_cols = if (length(cp_indices) > 0) list(mapping = map_lab),
    data_wide = mv_data_wide(X)
  )
}