R/diagnostics.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

# ---------------------------------------------------------------------------
# Theme E: seeing the statistic.
#
# Every detector computes something -- a CUSUM path, a MOSUM statistic at
# each bandwidth, a set of random intervals with contrast values -- and 0.4.0
# kept only the argmax. These three accessors (plus their plots) give the
# discarded object back, for the engines that expose it, and say plainly
# which engines do when one does not.
# ---------------------------------------------------------------------------

#' The detector's statistic as a function of location
#'
#' Returns the criterion the engine evaluated at each position, which is what
#' explains \emph{why} a changepoint landed where it did. Available for the
#' engines that expose it
#' (\code{subset(cpt_methods(), statistic)$method} lists them), and an error
#' naming those engines for the ones that do not.
#'
#' @param object A \code{ggcpt} object.
#' @return \code{cpt_statistic()} returns a tibble with \code{index},
#'   \code{statistic}, \code{threshold} (the engine's rejection threshold, or
#'   \code{NA}) and \code{label} (what the statistic is).
#'   \code{ggcpt_statistic()} returns a two-panel ggplot: the series above,
#'   the statistic below, sharing the x axis.
#' @seealso \code{\link{cpt_solution_path}()},
#'   \code{\link{cpt_scale_space}()}.
#' @export
#' @examplesIf requireNamespace("mosum", quietly = TRUE)
#' set.seed(2026)
#' fit <- cpt_detect(c(rnorm(200), rnorm(200, 3)), method = "mosum")
#' head(cpt_statistic(fit))
#' ggcpt_statistic(fit)
cpt_statistic <- function(object) {
  if (!is_ggcpt(object)) {
    cpt_abort("`object` must be a ggcpt object.", class = "bad_argument")
  }
  out <- extract_statistic(object)
  if (is.null(out)) {
    capable <- subset(cpt_methods(), statistic %in% TRUE)$method
    cpt_abort(dropped_fit_reason(object, capable, "the statistic") %||%
                paste0("Engine `", object$method, "` does not expose a ",
                       "per-location statistic. These do: ",
                       paste(capable, collapse = ", "), "."),
              class = "capability_absent")
  }
  out
}

# Internal: why an accessor found nothing on a result whose method does
# offer it: the result carries no engine object, which is where the thing
# is read from. Its refusal used to say the engine "does not expose" it and
# then list that same engine among the ones that do. NULL when that is not
# the reason.
#' @noRd
dropped_fit_reason <- function(object, capable, what) {
  method <- scalar_chr(object$method)
  if (!is.null(object$fit) || !method %in% capable) return(NULL)
  paste0("This `", method, "` result carries no engine object (it was made ",
         "with `keep_fit = FALSE`, or read back with cpt_import()), and ",
         what, " is read from that object. Re-run the detection with ",
         "`keep_fit = TRUE`.")
}

# Internal: per-engine statistic extraction. Returns NULL when the engine
# keeps nothing usable, so the caller can produce the "these engines do"
# message in one place.
#' @noRd
extract_statistic <- function(object) {
  n_full <- nrow(object$data)
  # A fit made with `na_action = "omit"` ran its engine on the observed
  # values only, so the engine's statistic is on that shorter series; it
  # is put back at the observed positions (NA at the gaps), where it used
  # to be padded as though a window had trimmed it, and read shifted.
  na_map <- object$na_map
  n <- if (!is.null(na_map)) length(na_map) else n_full
  fit <- object$fit
  method <- object$method

  mk <- function(stat, label, threshold = NA_real_) {
    stat <- as.numeric(stat)
    if (length(stat) != n) {
      # Some engines report a statistic on a shorter grid, because a moving
      # window trims BOTH ends. Left-aligning the pad therefore shifted
      # every value left by the bandwidth -- exactly the mis-alignment this
      # branch exists to prevent. Centre it instead: the two ends get the
      # same number of NAs (the extra one, when the shortfall is odd, goes
      # at the end), which is where a symmetric window has no value to
      # report. A statistic that is genuinely a prefix of the series is the
      # rarer case, and off-by-half-a-bandwidth in the readable direction.
      full <- rep(NA_real_, n)
      m <- min(n, length(stat))
      off <- (n - m) %/% 2L
      full[off + seq_len(m)] <- stat[seq_len(m)]
      stat <- full
    }
    if (!is.null(na_map)) {
      stat <- replace(rep(NA_real_, n_full), na_map, stat)
    }
    tibble::tibble(index = seq_len(n_full), statistic = stat,
                   threshold = as.numeric(threshold)[1], label = label)
  }

  # Sparse statistics (one value per candidate location) come back as a
  # length-n vector of zeros with the candidate values filled in, so the
  # panel is comparable with the dense ones.
  sparse <- function(loc, val, label, threshold = NA_real_) {
    stat <- rep(0, n)
    keep <- loc >= 1 & loc <= n & is.finite(val)
    if (any(keep)) {
      # A location can be the argmax of several intervals; keep the largest.
      ag <- stats::aggregate(list(v = abs(val[keep])),
                             by = list(i = loc[keep]), FUN = max)
      stat[ag$i] <- ag$v
    }
    mk(stat, label, threshold)
  }

  # `[[` with exact = TRUE throughout: `$` on a list partially matches, so
  # a renamed upstream field could silently resolve to a longer sibling
  # (`$stat` -> `statistics`) and be plotted as if it were the right one.
  fld <- function(nm) if (is.list(fit)) fit[[nm, exact = TRUE]] else NULL
  if (method == "mosum" && !is.null(fld("stat"))) {
    return(mk(fld("stat"), "MOSUM statistic",
              fld("threshold.value") %||% NA_real_))
  }
  if (method == "wbs" && !is.null(fld("res"))) {
    res <- fld("res")
    return(sparse(res[, "cpt"], res[, "CUSUM"],
                  "WBS CUSUM (max over intervals)"))
  }
  if (method == "not" && !is.null(fld("contrasts"))) {
    ct <- fld("contrasts")
    return(sparse(ct[, "arg.max"], ct[, "max.contrast"],
                  "NOT contrast (max over intervals)"))
  }
  if (method %in% c("pettitt", "buishand", "snht") &&
      inherits(fit, "htest") && !is.null(fld("data"))) {
    return(mk(fld("data"), paste0(toupper(substr(method, 1, 1)),
                               substring(method, 2), " statistic")))
  }
  if (method == "amoc") {
    return(mk(amoc_lr_profile(object$data$value),
              "AMOC standardised CUSUM profile"))
  }
  # np.mojo computes a per-location statistic and keeps it as `test.stat`
  # (and its cutoff as `threshold.val`); the registry did not claim it, so
  # nothing read it.
  if (method == "npmojo" && is.numeric(fld("test.stat"))) {
    thr <- suppressWarnings(as.numeric(fld("threshold.val") %||%
                                         NA_real_))[1]
    return(mk(fld("test.stat"), "NP-MOJO statistic", thr))
  }
  if (method == "nsp" && !is.null(object$regions)) {
    stat <- rep(0, n)
    reg <- object$regions
    if (nrow(reg) > 0 && !is.null(reg$value)) {
      for (i in seq_len(nrow(reg))) {
        stat[reg$start[i]:reg$end[i]] <- pmax(
          stat[reg$start[i]:reg$end[i]], abs(reg$value[i])
        )
      }
    }
    return(mk(stat, "NSP interval deviation",
              object$threshold_used %||% NA_real_))
  }
  prob <- posterior_prob_profile(object)
  if (!is.null(prob)) {
    return(mk(prob, "Posterior changepoint probability"))
  }
  diag_stat <- object$diagnostics[["statistic", exact = TRUE]]
  if (!is.null(diag_stat)) {
    return(mk(diag_stat[["statistic", exact = TRUE]] %||% diag_stat,
              diag_stat[["label", exact = TRUE]] %||% "Detector statistic",
              diag_stat[["threshold", exact = TRUE]] %||% NA_real_))
  }
  NULL
}

# Internal: the AMOC single-changepoint scan statistic for a change in mean
# with common variance. AMOC *is* this statistic maximised, so drawing it is
# drawing the method rather than inventing a diagnostic.
#
# It is a standardised CUSUM divided by the WHOLE-SERIES sd, not a
# log-likelihood ratio, and it used to be labelled as one. The argmax is
# unaffected -- sd(y) is a constant in k -- so the peak location was always
# honest; the y-axis values were not, because under a real change sd(y) is
# inflated by the change itself and the profile is systematically
# compressed. The label now says which quantity is drawn.
#' @noRd
amoc_lr_profile <- function(y) {
  n <- length(y)
  cs <- cumsum(y)
  tot <- cs[n]
  k <- seq_len(n - 1L)
  # Standardised CUSUM: |sqrt(n/(k(n-k))) * (S_k - k/n * S_n)|
  stat <- abs(sqrt(n / (k * (n - k))) * (cs[k] - k * tot / n))
  s <- stats::sd(y)
  if (is.finite(s) && s > 0) stat <- stat / s
  c(stat, NA_real_)
}

#' @rdname cpt_statistic
#' @export
ggcpt_statistic <- function(object) {
  st <- cpt_statistic(object)
  idx_vals <- plot_index(object)
  x_lab <- plot_index_label(object)
  lab <- st$label[1]

  top <- tibble::tibble(index = idx_vals, y = object$data$value,
                        panel = "Series")
  bottom <- tibble::tibble(index = idx_vals, y = st$statistic,
                           panel = lab)
  both <- rbind(top, bottom)
  both$panel <- factor(both$panel, levels = c("Series", lab))

  p <- ggplot2::ggplot(both, ggplot2::aes(index, y)) +
    ggplot2::geom_line(data = both[both$panel == "Series", ],
                       colour = "grey40") +
    ggplot2::geom_line(data = both[both$panel == lab, ],
                       colour = "#0072B2", na.rm = TRUE) +
    ggplot2::facet_grid(panel ~ ., scales = "free_y", switch = "y") +
    ggplot2::labs(x = x_lab, y = NULL,
                  title = paste0("Detector statistic (", object$method, ")"))

  if (nrow(object$changepoints) > 0) {
    p <- p + ggplot2::geom_vline(
      xintercept = idx_vals[object$changepoints$cp],
      colour = "blue", linetype = "dashed", linewidth = 0.4
    )
  }
  thr <- st$threshold[1]
  if (is.finite(thr)) {
    p <- p + ggplot2::geom_hline(
      data = tibble::tibble(yint = thr, panel = factor(lab,
                                                       levels = levels(both$panel))),
      ggplot2::aes(yintercept = yint),
      linetype = "dotted", colour = "grey30"
    )
  }
  p
}

#' The solution path of a search-based detector
#'
#' The order in which candidate changepoints entered the model, with the
#' contrast (or split criterion) at each step. Binary segmentation splits
#' recursively, WBS/WBS2/NOT/TGUH rank random intervals; in every case the
#' final answer is a prefix of a path, and seeing the path shows how
#' decisively each changepoint beat the next.
#'
#' @param object A \code{ggcpt} object from an engine with a solution path
#'   (\code{subset(cpt_methods(), path)$method}).
#' @return \code{cpt_solution_path()} returns a tibble with \code{step},
#'   \code{cp}, \code{contrast} and, for interval-based searches,
#'   \code{start}/\code{end} of the interval that proposed it, plus a
#'   \code{selected} flag marking the changepoints in the final model.
#'   \code{ggcpt_solution_path()} draws it.
#'
#'   \code{contrast} is the engine's own ordering criterion, and it is a
#'   \emph{different quantity} per engine: a penalty value for
#'   \code{binseg}/\code{segneigh} (from
#'   \code{changepoint::pen.value.full()}), \eqn{|CUSUM|} for \code{wbs},
#'   \eqn{|}\code{max.contrast}\eqn{|} for \code{not}, and
#'   \pkg{breakfast}'s candidate criterion for \code{wbs2}/\code{tguh}.
#'   The values order the candidates within one result; they are not
#'   comparable across engines, and the plot legend names the quantity
#'   rather than calling all of them "Contrast".
#'
#'   For \code{wbs2} and \code{tguh} the path is \strong{recomputed} with
#'   \pkg{breakfast}, because their fit objects do not keep the candidate
#'   list. \code{wbs2}'s search is randomised, so its path is a second
#'   search of the same series rather than a record of the first: it can
#'   differ between calls, and \code{selected} can be \code{FALSE}
#'   throughout if the recomputed candidates miss the fit's own
#'   changepoints. Every other engine's path is read off the fit.
#' @seealso \code{\link{cpt_statistic}()}, \code{\link{cpt_crops}()} for the
#'   penalty path of an optimal-partitioning method.
#' @export
#' @examples
#' set.seed(2026)
#' fit <- cpt_detect(c(rnorm(150), rnorm(150, 3)), method = "binseg")
#' cpt_solution_path(fit)
cpt_solution_path <- function(object) {
  if (!is_ggcpt(object)) {
    cpt_abort("`object` must be a ggcpt object.", class = "bad_argument")
  }
  out <- extract_solution_path(object)
  if (is.null(out)) {
    capable <- subset(cpt_methods(), path %in% TRUE)$method
    cpt_abort(dropped_fit_reason(object, capable, "the solution path") %||%
                paste0("Engine `", object$method, "` does not expose a ",
                       "solution path. These do: ",
                       paste(capable, collapse = ", "), ". For the penalty ",
                       "path of an optimal-partitioning method see ",
                       "cpt_crops()."),
              class = "capability_absent")
  }
  out
}

#' @noRd
extract_solution_path <- function(object) {
  fit <- object$fit
  method <- object$method
  # The engine's locations are on the series it saw: the observed values
  # of a fit made with gaps, mapped back here (see extract_statistic()).
  na_map <- object$na_map
  n <- if (!is.null(na_map)) length(na_map) else nrow(object$data)
  selected <- object$changepoints$cp
  # Exact extraction, for the same reason as in extract_statistic().
  pfld <- function(nm) if (is.list(fit)) fit[[nm, exact = TRUE]] else NULL

  finish <- function(cp, contrast, start = NA_integer_, end = NA_integer_) {
    keep <- !is.na(cp) & cp >= 1 & cp < n
    cp <- cp[keep]
    contrast <- contrast[keep]
    if (length(start) > 1) start <- start[keep]
    if (length(end) > 1) end <- end[keep]
    if (length(cp) == 0) return(NULL)
    if (!is.null(na_map)) {
      back <- function(v) {
        out <- rep(NA_integer_, length(v))
        ok <- !is.na(v) & v >= 1 & v <= length(na_map)
        out[ok] <- na_map[v[ok]]
        out
      }
      cp <- back(cp)
      start <- back(start)
      end <- back(end)
    }
    tibble::tibble(step = seq_along(cp), cp = as.integer(cp),
                   contrast = as.numeric(contrast),
                   start = as.integer(start), end = as.integer(end),
                   selected = cp %in% selected)
  }

  if (method %in% c("binseg", "segneigh") && inherits(fit, "cpt.range")) {
    full <- changepoint::cpts.full(fit)
    if (is.null(dim(full))) full <- matrix(full, nrow = 1)
    pens <- as.numeric(changepoint::pen.value.full(fit))
    # Row i of cpts.full holds the segmentation with i changepoints. For
    # binary segmentation the rows are nested, so the one new entry per row
    # is the split that step added. Segment Neighbourhood is NOT nested --
    # it re-solves the dynamic program at each K, so consecutive rows can
    # differ by more than one changepoint -- and keeping only `new[1]`
    # discarded the rest, leaving the segneigh path missing candidates with
    # a `step` numbering that did not match its own rows. Emitting all of
    # `new` is a no-op for binseg, where `new` has length one.
    order_cp <- integer(0)
    contrast <- numeric(0)
    prev <- integer(0)
    for (i in seq_len(nrow(full))) {
      row <- as.integer(full[i, ][!is.na(full[i, ])])
      new <- setdiff(row, prev)
      if (length(new) > 0) {
        order_cp <- c(order_cp, new)
        contrast <- c(contrast,
                      rep(if (i <= length(pens)) pens[i] else NA_real_,
                          length(new)))
      }
      prev <- row
    }
    return(finish(order_cp, contrast))
  }
  # strucchange computes the optimal segmentation for every number of
  # breaks (the Bai-Perron table) and the result kept only the chosen one.
  # Like SegNeigh's, the solutions are not nested, so each row contributes
  # the breaks it adds; the contrast is that solution's residual sum of
  # squares.
  if (method == "strucchange" && inherits(fit, "breakpoints") &&
      requireNamespace("strucchange", quietly = TRUE)) {
    sm <- tryCatch(summary(fit), error = function(e) NULL)
    bpm <- if (!is.null(sm)) sm$breakpoints
    if (is.matrix(bpm) && nrow(bpm) > 0) {
      rss <- as.numeric(sm$RSS["RSS", ])
      order_cp <- integer(0)
      contrast <- numeric(0)
      prev <- integer(0)
      for (i in seq_len(nrow(bpm))) {
        row <- as.integer(bpm[i, ][!is.na(bpm[i, ])])
        new <- setdiff(row, prev)
        if (length(new) > 0) {
          order_cp <- c(order_cp, new)
          contrast <- c(contrast, rep(if (i + 1L <= length(rss)) rss[i + 1L]
                                      else NA_real_, length(new)))
        }
        prev <- row
      }
      return(finish(order_cp, contrast))
    }
  }
  if (method == "wbs" && !is.null(pfld("res"))) {
    res <- pfld("res")
    ord <- order(-abs(res[, "CUSUM"]))
    res <- res[ord, , drop = FALSE]
    # One row per distinct candidate location, strongest first.
    dup <- duplicated(res[, "cpt"])
    res <- res[!dup, , drop = FALSE]
    return(finish(res[, "cpt"], abs(res[, "CUSUM"]),
                  res[, "s"], res[, "e"]))
  }
  if (method == "not" && !is.null(pfld("contrasts"))) {
    ct <- pfld("contrasts")
    ord <- order(-abs(ct[, "max.contrast"]))
    ct <- ct[ord, , drop = FALSE]
    dup <- duplicated(ct[, "arg.max"])
    ct <- ct[!dup, , drop = FALSE]
    return(finish(ct[, "arg.max"], abs(ct[, "max.contrast"]),
                  ct[, "start"], ct[, "end"]))
  }
  if (method %in% c("wbs2", "tguh") &&
      requireNamespace("breakfast", quietly = TRUE)) {
    # The only branch that RECOMPUTES rather than reading `object$fit`,
    # because breakfast's fit object does not keep the candidate list. wbs2
    # is randomised (see the reproducibility note in ?wbs2_wrapper), so two
    # calls on the same object can return different paths, the recomputed
    # candidates need not contain the fit's own changepoints -- in which
    # case `selected` comes back all FALSE and the plot shows a path with
    # nothing selected -- and any `seed` given to the original cpt_detect()
    # is not in scope here. Say so rather than presenting a second search
    # as the first one's path.
    if (identical(method, "wbs2")) {
      cpt_warn("`wbs2`'s solution path is recomputed here: breakfast's fit ",
               "does not keep its candidate list, and the search is ",
               "randomised, so this path is a second search of the same ",
               "series and need not contain every changepoint in the ",
               "result. Read it as which splits the method considers ",
               "strong, not as the exact path behind this fit.",
               class = "warning")
    }
    sol <- tryCatch(
      if (method == "wbs2") {
        breakfast::sol.wbs2(object$data$value)
      } else {
        breakfast::sol.tguh(object$data$value)
      },
      error = function(e) NULL
    )
    if (!is.null(sol) && !is.null(sol$cands)) {
      cands <- sol$cands
      return(finish(cands[, 3], cands[, 4], cands[, 1], cands[, 2]))
    }
  }
  diag_path <- object$diagnostics[["solution_path", exact = TRUE]]
  if (!is.null(diag_path)) {
    # A registered detector's path used to be returned untouched, bypassing
    # the contract every built-in branch goes through: one without `cp`
    # reached ggcpt_solution_path() and died at `path$cp_x <-
    # idx_vals[path$cp]`, one with an out-of-range `cp` got NA x-positions,
    # and one without `selected` broke the colour scale. extract_statistic()
    # already defaults its three fields the same way.
    diag_path <- tibble::as_tibble(diag_path)
    if (!"cp" %in% names(diag_path)) {
      cpt_abort("A registered method's `solution_path` needs a `cp` column ",
                "holding the candidate locations; this one has ",
                if (ncol(diag_path) == 0) "no columns" else
                  paste0("(", paste(names(diag_path), collapse = ", "), ")"),
                ". The columns are described in the `@return` of ",
                "?cpt_solution_path.", class = "engine_error")
    }
    fill <- function(nm, default) {
      v <- if (nm %in% names(diag_path)) diag_path[[nm]] else default
      if (length(v) == 1L) rep(v, nrow(diag_path)) else v
    }
    return(finish(diag_path$cp, fill("contrast", NA_real_),
                  fill("start", NA_integer_), fill("end", NA_integer_)))
  }
  NULL
}

#' @rdname cpt_solution_path
#' @param max_steps Longest prefix of the path drawn. Defaults to \code{40}:
#'   a randomised search proposes hundreds of candidates and only the head of
#'   the ranking is readable.
#' @export
ggcpt_solution_path <- function(object, max_steps = 40) {
  path <- cpt_solution_path(object)
  validate_scalar(max_steps, "max_steps", min = 1)
  path <- utils::head(path, as.integer(max_steps))
  idx_vals <- plot_index(object)
  x_lab <- plot_index_label(object)
  path$cp_x <- idx_vals[path$cp]
  has_interval <- !all(is.na(path$start))
  if (has_interval) {
    path$start_x <- idx_vals[pmax(1L, pmin(path$start, length(idx_vals)))]
    path$end_x <- idx_vals[pmax(1L, pmin(path$end, length(idx_vals)))]
  }

  p <- ggplot2::ggplot(path, ggplot2::aes(cp_x, step))
  if (has_interval) {
    p <- p + ggplot2::geom_segment(
      ggplot2::aes(x = start_x, xend = end_x, y = step, yend = step,
                   colour = selected),
      linewidth = 0.5, alpha = 0.7
    )
  }
  p +
    ggplot2::geom_point(ggplot2::aes(colour = selected, size = contrast)) +
    ggplot2::scale_y_reverse(breaks = scales_int_breaks(path$step)) +
    ggplot2::scale_colour_manual(values = c(`FALSE` = "grey60",
                                            `TRUE` = "#D55E00"),
                                 name = "In final model") +
    # The `contrast` column carries a different quantity per engine -- a
    # penalty value for binseg/segneigh, |CUSUM| for wbs, |max.contrast| for
    # not, breakfast's own candidate criterion for wbs2/tguh -- and one
    # legend title of "Contrast" claimed they were all the same thing.
    ggplot2::scale_size_continuous(range = c(0.6, 3),
                                   name = contrast_label(object$method)) +
    ggplot2::labs(x = x_lab, y = "Step in the solution path",
                  title = paste0("Solution path (", object$method, ")"),
                  subtitle = if (has_interval) {
                    "Bars show the interval that proposed each candidate"
                  } else {
                    NULL
                  })
}

# Internal: what the `contrast` column of a solution path actually holds,
# which depends on the engine. Used for the plot legend so the reader is not
# told a penalty value and a CUSUM are the same quantity.
#' @noRd
contrast_label <- function(method) {
  switch(method %||% "",
    binseg = "Penalty value",
    segneigh = "Penalty value",
    wbs = "|CUSUM|",
    not = "|Contrast|",
    wbs2 = "Split criterion",
    tguh = "Split criterion",
    "Contrast")
}

# Internal: integer axis breaks without pulling in the scales package.
#' @noRd
scales_int_breaks <- function(v) {
  r <- range(v)
  by <- max(1L, ceiling(diff(r) / 10))
  seq(r[1], r[2], by = by)
}

#' Scale space: the statistic across bandwidths
#'
#' Sweeps a multiscale detector's bandwidth and returns the statistic at
#' every (location, bandwidth) pair. The resulting heatmap answers the
#' question a single-bandwidth fit cannot: \emph{at which resolutions does
#' this feature exist?} A change that is significant only at a wide bandwidth
#' is a slow shift; one that appears only at a narrow bandwidth is a spike.
#'
#' @param x A numeric vector (or, for \code{method = "npmojo"}, a matrix
#'   with rows as time points), or any \code{ggcpt} object. The sweep is
#'   computed here rather than read off the result, so a fit from any
#'   engine is accepted and the scale space is the one belonging to its
#'   series: from a \code{mosum} or \code{npmojo} result the engine is
#'   adopted as well (see \code{method}), and from anything else the
#'   default \code{"mosum"} sweep is drawn over the same series.
#' @param bandwidths Integer vector of bandwidths to sweep. Defaults to a
#'   geometric grid between \code{max(5, n/50)} and \code{n/4}.
#' @param method Which engine to sweep: \code{"mosum"} (default) or
#'   \code{"npmojo"}. Taken from \code{x} when it is a \code{ggcpt} produced
#'   by one of them.
#' @param ... Additional arguments passed to the engine at each bandwidth.
#'
#' @return \code{cpt_scale_space()} returns a tibble with \code{index},
#'   \code{bandwidth}, \code{statistic}, \code{threshold},
#'   \code{significant} (the statistic at that location and bandwidth
#'   exceeds the threshold) and \code{detected} (the engine reported a
#'   changepoint there at that bandwidth). The two differ: a location can
#'   clear the threshold without surviving the engine's own pruning. \code{ggcpt_scale_space()} draws the heatmap with
#'   the accepted changepoints overlaid.
#' @seealso \code{\link{cpt_statistic}()}.
#' @export
#' @examplesIf requireNamespace("mosum", quietly = TRUE)
#' set.seed(2026)
#' x <- c(rnorm(300), rnorm(300, 2))
#' ss <- cpt_scale_space(x, bandwidths = c(20, 40, 80))
#' head(ss)
cpt_scale_space <- function(x, bandwidths = NULL,
                            method = c("mosum", "npmojo"), ...) {
  obs_pos <- NULL
  if (is_ggcpt(x)) {
    if (missing(method) && x$method %in% c("mosum", "npmojo")) {
      method <- x$method
    }
    # A multivariate result keeps its coordinates: npmojo sweeps them all,
    # and flattening to the first one would sweep a different series from
    # the one that was detected on.
    series <- if (!is.null(x$data_wide) && n_coordinates(x) > 1) {
      as.matrix(x$data_wide[, setdiff(names(x$data_wide),
                                      c("index", "index_value")),
                            drop = FALSE])
    } else {
      x$data$value
    }
    # A fit made with `na_action = "omit"` keeps its gaps, and the sweep
    # handed them to the engine: mosum is measured to lose changepoints on
    # a missing value without saying so. The sweep runs on the observed
    # rows and reports their original positions.
    observed <- if (is.matrix(series)) stats::complete.cases(series) else
      !is.na(series)
    if (!all(observed)) {
      obs_pos <- which(observed)
      series <- if (is.matrix(series)) series[obs_pos, , drop = FALSE] else
        series[obs_pos]
    }
  } else {
    validate_data(x)
    series <- x
  }
  method <- cpt_match_arg(method)

  # Shape first, engine second: telling someone to install `mosum` for input
  # `mosum` cannot accept anyway is the wrong complaint, and it makes the
  # error depend on what happens to be installed.
  #
  # mosum is univariate; npmojo is not, and refusing its matrix here would
  # make the scale-space view unavailable for the one multiscale engine that
  # needs it most.
  if (method == "mosum") {
    # Pass `method`, not a literal: as_uni_vector() builds "Method `%s` is
    # univariate", and a literal "scale_space" named this function rather
    # than the engine the caller actually asked for -- so someone who wrote
    # `method = "mosum"` was told about a "method" they had never heard of.
    series <- as_uni_vector(series, method)
  } else if (is.data.frame(series)) {
    series <- as_mv_matrix(series, arg = "series")
  }
  n <- NROW(series)

  if (is.null(bandwidths)) {
    lo <- max(5, floor(n / 50))
    hi <- max(lo + 1, floor(n / 4))
    bandwidths <- unique(as.integer(round(
      exp(seq(log(lo), log(hi), length.out = 8))
    )))
  }
  # Checked before the coercion: text or Inf became NA with a warning and
  # was reported as "No usable bandwidth" for a series that has plenty.
  validate_grid(bandwidths, "bandwidths", min = 1)
  bandwidths <- sort(unique(as.integer(bandwidths)))
  bandwidths <- bandwidths[bandwidths >= 2 & 2 * bandwidths < n]
  if (length(bandwidths) == 0) {
    cpt_abort("No usable bandwidth: a moving window needs `2 * G < n`, and ",
               "the ", "series has ", n, " observations.",
              class = "short_series")
  }
  need_pkg(if (method == "mosum") "mosum" else "CptNonPar")

  # Why the engine refused is kept, not discarded. Every bandwidth failing
  # used to report "No bandwidth produced a usable fit", which reads as a
  # statement about the bandwidths -- and it is the wrong diagnosis for the
  # commonest cause. `...` goes straight to the engine, so a single
  # unsupported argument makes every fit fail at once:
  # `cpt_scale_space(x, bandwidths = c(20, 40))` works and
  # `cpt_scale_space(x, bandwidths = c(20, 40), index = dates)` did not,
  # because `index` is not one of `mosum::mosum()`'s arguments. Telling
  # someone their bandwidths are unusable when the engine said "unused
  # argument (index = ...)" sends them to re-tune the sweep.
  why <- NULL
  keep_why <- function(e) {
    if (is.null(why)) {
      # R's own "unused argument (index = c(...))" deparses the value, so a
      # rejected 180-point index arrives as 180 numbers. Keep the first
      # clause, which is the part that names what was wrong.
      msg <- gsub("[[:space:]]+", " ", conditionMessage(e))
      why <<- if (nchar(msg) > 120) paste0(substr(msg, 1, 117), "...") else msg
    }
    NULL
  }
  rows <- lapply(bandwidths, function(G) {
    if (method == "mosum") {
      fit <- tryCatch(mosum::mosum(series, G = G, ...), error = keep_why)
      if (is.null(fit)) return(NULL)
      # Exact [[ ]] for the same reason as the np.mojo branch below: mosum
      # returns both `threshold` (the rule) and `threshold.value` (the
      # number), and the same file argues the case twice and then read this
      # one with a bare `$`.
      fm <- function(nm) if (is.list(fit)) fit[[nm, exact = TRUE]] else NULL
      thr <- as.numeric(fm("threshold.value") %||% NA_real_)
      stat <- as.numeric(fm("stat") %||% rep(NA_real_, n))
      cp <- as.integer(fm("cpts") %||% integer(0))
    } else {
      fit <- tryCatch(CptNonPar::np.mojo(series, G = G, ...),
                      error = keep_why)
      if (is.null(fit)) return(NULL)
      # np.mojo names its per-location statistic `test.stat`, and
      # `threshold` is the *rule* ("bootstrap"), not the number -- that is
      # `threshold.val`. Reading `$threshold` here would put a character
      # string where a cutoff belongs and make every cell insignificant.
      fl <- function(nm) if (is.list(fit)) fit[[nm, exact = TRUE]] else NULL
      stat <- as.numeric(fl("test.stat") %||% rep(NA_real_, n))
      thr <- suppressWarnings(as.numeric(fl("threshold.val") %||% NA_real_))[1]
      cp <- as.integer(fl("cpts") %||% integer(0))
    }
    if (length(stat) != n) {
      full <- rep(NA_real_, n)
      full[seq_len(min(n, length(stat)))] <- stat[seq_len(min(n, length(stat)))]
      stat <- full
    }
    tibble::tibble(index = seq_len(n), bandwidth = G, statistic = stat,
                   threshold = thr,
                   significant = is.finite(stat) & is.finite(thr) &
                     stat > thr,
                   detected = seq_len(n) %in% cp)
  })
  rows <- Filter(Negate(is.null), rows)
  if (length(rows) == 0) {
    detail <- if (is.null(why)) {
      paste0(": the engine returned nothing usable at any of the ",
             length(bandwidths), " bandwidth(s) tried.")
    } else {
      paste0(". All ", length(bandwidths), " bandwidth(s) failed, and `",
             method, "` gave the same reason each time: ", why,
             " Arguments in `...` go straight to the engine, so if one of ",
             "those is unexpected, check it against ?",
             if (method == "mosum") "mosum::mosum" else "CptNonPar::np.mojo",
             " rather than re-tuning the sweep.")
    }
    cpt_abort("No bandwidth produced a usable fit for method `", method, "`",
              detail, class = "engine_error")
  }
  out <- do.call(rbind, rows)
  if (!is.null(obs_pos)) {
    # Back to original positions, with the gaps as empty cells so the
    # heatmap's grid stays regular.
    out$index <- obs_pos[out$index]
    full <- expand.grid(index = seq_len(length(observed)),
                        bandwidth = unique(out$bandwidth))
    key <- paste(full$index, full$bandwidth)
    m <- match(key, paste(out$index, out$bandwidth))
    filled <- out[m, , drop = FALSE]
    filled$index <- full$index
    filled$bandwidth <- full$bandwidth
    filled$threshold <- out$threshold[match(full$bandwidth, out$bandwidth)]
    filled$significant[is.na(m)] <- FALSE
    filled$detected[is.na(m)] <- FALSE
    out <- filled[order(filled$bandwidth, filled$index), , drop = FALSE]
  }
  attr(out, "method") <- method
  out
}

#' @rdname cpt_scale_space
#' @export
ggcpt_scale_space <- function(x, bandwidths = NULL,
                              method = c("mosum", "npmojo"), ...) {
  # Passing `method` straight through would make `missing(method)` FALSE
  # inside cpt_scale_space(), so a result from npmojo would be swept with
  # mosum. Resolve the default here instead.
  if (missing(method) && is_ggcpt(x) && x$method %in% c("mosum", "npmojo")) {
    method <- x$method
  }
  ss <- cpt_scale_space(x, bandwidths = bandwidths, method = method, ...)
  eng <- attr(ss, "method")
  idx_vals <- if (is_ggcpt(x)) plot_index(x) else seq_len(max(ss$index))
  x_lab <- if (is_ggcpt(x)) plot_index_label(x) else "Index"
  ss$index_x <- idx_vals[ss$index]

  det <- ss[ss$detected, , drop = FALSE]

  p <- ggplot2::ggplot(ss, ggplot2::aes(index_x, factor(bandwidth),
                                        fill = statistic)) +
    ggplot2::geom_raster(na.rm = TRUE) +
    ggplot2::scale_fill_viridis_c(option = "D", na.value = "grey95",
                                  name = "Statistic") +
    ggplot2::labs(x = x_lab, y = "Bandwidth (G)",
                  title = paste0("Scale space (", eng, ")"),
                  subtitle = "Points mark the changepoints accepted at each bandwidth")
  if (nrow(det) > 0) {
    p <- p + ggplot2::geom_point(data = det, colour = "white", size = 1.6,
                                 shape = 21, fill = "#D55E00")
  }
  p
}