R/inference.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

# ---------------------------------------------------------------------------
# Theme A: inference. Two generics with one contract each.
#
#   cpt_confint()  "where could this changepoint be?" -- four provenances
#                  (native / nsp / bootstrap / posterior) behind one tibble,
#                  with a `source` column so the difference is visible in the
#                  data and not only in the docs.
#   cpt_test()     "is this changepoint real?" -- native tests where an
#                  engine has one, an explicitly unadjusted fallback where it
#                  does not, and a `selection_adjusted` column that never
#                  lets the two be confused.
# ---------------------------------------------------------------------------

#' Confidence intervals for changepoint locations
#'
#' Answers "where could this changepoint be?" for any result, and says which
#' of five routes it used. \code{show_ci = TRUE} in \code{autoplot()} works
#' only for the handful of engines that ship intervals of their own; this
#' generic covers the rest, and (because the four routes mean genuinely
#' different things) reports the provenance in a \code{source} column
#' rather than presenting them as interchangeable.
#'
#' @param object A \code{ggcpt} object.
#' @param level Confidence/credible level. Defaults to \code{0.95}. Honoured
#'   by the routes that compute an interval (\code{"bootstrap"},
#'   \code{"posterior"} and \code{"nsp"}) and \strong{ignored by}
#'   \code{"native"}, which reports the interval the engine already
#'   computed at whatever level it was asked for. Note that
#'   \code{method = "auto"} resolves to \code{"native"} whenever the engine
#'   supplied one, so an explicit \code{level} can go unused there too: it
#'   is reported in the \code{level} column either way, and supplying a
#'   level the answer does not carry now warns rather than passing
#'   silently. To choose the level yourself, name a computing route.
#' @param method Which route to use:
#'   \describe{
#'     \item{\code{"auto"}}{(default) native if the engine supplied
#'       intervals, else posterior if it supplied one, else conformal where
#'       the change is one between exchangeable segments (a change in mean,
#'       variance, both, or distribution), else bootstrap.}
#'     \item{\code{"native"}}{the engine's own \code{ci_lower}/\code{ci_upper}
#'       (SMUCE/HSMUCE simultaneous confidence sets, \pkg{strucchange}
#'       break-date intervals, \pkg{segmented} breakpoint intervals,
#'       \pkg{mcp} posterior intervals, \pkg{bfast} break confidence
#'       intervals, and \code{taylor}'s bootstrap confidence limits). The
#'       engines that supply them are the ones
#'       \code{\link{cpt_methods}()} marks in its \code{ci} column;
#'       \code{nsp} is marked there too but is reported under its own
#'       provenance below, because its regions are not intervals around an
#'       estimate.}
#'     \item{\code{"nsp"}}{Narrowest Significance Pursuit regions computed on
#'       the same series and matched to the changepoints. These are
#'       \emph{not} intervals around an estimate (see
#'       \code{\link{nsp_wrapper}()}), but they are the strongest guarantee
#'       available, so the mapping is reported as
#'       \code{source = "nsp_region"} and a changepoint in no region gets
#'       \code{NA}.}
#'     \item{\code{"bootstrap"}}{within-segment residual resampling (the
#'       \code{\link{cpt_stability}()} scheme), re-running the detector and
#'       taking quantiles of the re-detected location. Model-agnostic and
#'       available for every engine; it measures the sampling variability of
#'       the \emph{procedure}, conditional on the fitted segmentation, and it
#'       is not exact.}
#'     \item{\code{"conformal"}}{a distribution-free confidence set
#'       (CONCH, Hore and Ramdas 2026). Each changepoint is examined in its
#'       window between the neighbouring changepoints; every candidate
#'       split of the window is tested by permuting the observations within
#'       each side of it, and the candidates that survive at
#'       \code{1 - level} form a set that contains the true location with
#'       probability at least \code{level}, whatever the distribution of
#'       the data, provided the observations within each segment are
#'       exchangeable (independent and identically distributed is the
#'       usual case) and the window holds one change. It needs no re-run of
#'       the detector, so it serves results built with
#'       \code{\link{as_ggcpt}()} and registered methods too. The set need
#'       not be contiguous: \code{ci_lower} and \code{ci_upper} are its
#'       hull and \code{set_size} counts its members, and the p-value of
#'       every candidate is attached as the \code{"profile"} attribute.
#'       The score matches what changed: a mean shift, a variance change,
#'       both, or a rank statistic for a change in distribution; summed over
#'       coordinates for a multivariate result. \code{B} sets the number of
#'       permutations per candidate.}
#'     \item{\code{"posterior"}}{a credible interval from the engine's
#'       posterior changepoint-probability profile (\code{bcp},
#'       \code{beast}): the narrowest \emph{contiguous} set of positions
#'       around the estimate holding \code{level} of the posterior
#'       changepoint mass in that changepoint's window, where the window is
#'       bounded by the neighbouring changepoints so two of them cannot
#'       claim the same mass twice.
#'
#'       Expect these to be wide, and read the width as a statement about
#'       the profile rather than about the location. Both supplying engines
#'       put roughly two-thirds of the window's mass at the estimate itself
#'       and spread the remaining third as a thin floor across every other
#'       position, so reaching a high \code{level} means swallowing that
#'       floor. Measured on a 200-point series with one clean change:
#'       \code{level = 0.5} gives a width of 0 (the mode alone holds more
#'       than half), \code{0.8} gives 72 (\pkg{bcp}) and 91
#'       (\pkg{beast}), and \code{0.95} gives 166 and 187, which is 83\%
#'       and 94\% of the series. The requested level \emph{is} delivered in
#'       each case; what a wide interval says is that the posterior did not
#'       localise the change, not that the location is uncertain by that
#'       much. \code{cpt_confint()} warns when an interval covers more than
#'       half its window, for exactly that reason.}
#'   }
#'
#' @section How well these cover:
#' Measured over 120 replicates on a 200-point series with one changepoint
#' at 100 and a jump of three standard deviations, at a nominal level of
#' 0.95: \code{"bootstrap"} on \code{pelt} covered 0.992 of the time with
#' a mean width of 2.2; \pkg{strucchange}'s native intervals covered 1.000
#' at width 4.4; \pkg{stepR}'s (\code{smuce}) covered 0.992 at width 4.6;
#' and \code{"posterior"} on \pkg{bcp} covered 1.000 at width 157. Every
#' route is \emph{conservative} (none under-covers), and the width is
#' what separates them. The two native routes and the bootstrap are the
#' ones to quote; see the note on \code{"posterior"} above for why its
#' interval is so much wider.
#'
#' The conformal route was measured with \code{\link{cpt_confint_coverage}()}
#' on the same design (300 replicates per cell, \code{pelt}, level 0.95).
#' For a jump of three standard deviations it covered 0.960 at a mean width
#' of 0.5 (the bootstrap: 0.973 at 2.0); for a jump of one standard
#' deviation 0.970 at 18.6, where the bootstrap fell short of the level at
#' 0.926 (width 25.8); under \eqn{t_3} noise with a jump of two, 0.963 at
#' 3.5 (bootstrap 0.943 at 4.4). Its guarantee assumes exchangeable
#' observations within each segment, and under AR(1) noise with
#' \eqn{\rho = 0.3} it covered 0.920 (bootstrap 0.887): dependent noise
#' shortens every interval here, and \code{\link{cpt_assumptions}()}
#' checks for it.
#' @param B Bootstrap replicates for \code{method = "bootstrap"}, and
#'   permutations per candidate split for \code{method = "conformal"}.
#'   Defaults to \code{200}.
#' @param seed Optional seed (bootstrap and NSP are both random). The seed
#'   is scoped to this call: \code{.Random.seed} is saved and restored, so a
#'   seeded call inside a simulation loop does not pin the loop's own
#'   stream.
#' @param ... Passed to \code{\link{cpt_detect}()} on the bootstrap
#'   replicates, or to \code{\link{nsp_wrapper}()}.
#'
#' @return A tibble with one row per changepoint: \code{cp},
#'   \code{ci_lower}, \code{ci_upper}, \code{level}, \code{source}, plus
#'   \code{cp_index}/\code{ci_lower_index}/\code{ci_upper_index} on the
#'   original scale when the result carries a time index, and
#'   \code{n_replicates} for \code{method = "bootstrap"} (how many
#'   replicates actually contributed a draw for that changepoint), and
#'   \code{set_size} for \code{method = "conformal"} (how many locations
#'   the confidence set holds; an empty set has \code{NA} bounds and
#'   warns, because no single change between exchangeable segments fits
#'   that window).
#'   \code{method = "nsp"} returns \code{NA} bounds for a changepoint that
#'   falls in no region, which is a finding rather than a failure.
#' @references
#' \insertRef{hore2026conch}{ggchangepoint}
#' @seealso \code{\link{cpt_confint_coverage}()}, which measures how often
#'   each route covers; \code{\link{cpt_test}()},
#'   \code{\link{cpt_stability}()}, \code{\link{nsp_wrapper}()}.
#' @export
#' @examples
#' set.seed(2026)
#' x <- c(rnorm(80), rnorm(80, 4))
#' fit <- cpt_detect(x, method = "pelt")
#' cpt_confint(fit, method = "conformal", B = 50, seed = 1)
#' cpt_confint(fit, method = "bootstrap", B = 25, seed = 1)
cpt_confint <- function(object, level = 0.95,
                        method = c("auto", "native", "nsp", "bootstrap",
                                   "posterior", "conformal"),
                        B = 200, seed = NULL, ...) {
  if (!is_ggcpt(object)) {
    cpt_abort("`object` must be a ggcpt object.", class = "bad_argument")
  }
  method <- cpt_match_arg(method)
  # Whether the caller asked for a level, as opposed to taking the default.
  # Used at the end: only an explicit request is worth warning about.
  level_supplied <- !missing(level)
  requested <- if (level_supplied) level else NA_real_
  validate_scalar(level, "level", min = 0, max = 1,
                  min_open = TRUE, max_open = TRUE)
  validate_scalar(B, "B", min = 1)

  cp <- object$changepoints$cp
  if (length(cp) == 0) {
    return(empty_confint(object))
  }

  native <- native_bounds(object)
  has_native <- !is.null(native)
  has_post <- !is.null(posterior_prob_profile(object))

  if (method == "auto") {
    method <- if (has_native) {
      "native"
    } else if (has_post) {
      "posterior"
    } else if (conformal_possible(object)) {
      "conformal"
    } else {
      "bootstrap"
    }
    # The bootstrap route re-runs the detector by name, which a result built
    # by as_ggcpt() from an outside source cannot support. Say that in terms
    # of what the caller asked for rather than letting the bootstrap branch
    # complain about a method the caller never named.
    if (method == "bootstrap" && !bootstrap_possible(object)) {
      cpt_abort("`method = \"auto\"` has nothing to fall back on for this ",
                "result: the engine supplied no interval and no posterior, the ",
                "change is not one between exchangeable segments that the ",
                "conformal set needs, and the bootstrap route needs to re-run ",
                "the detector, but `",
                scalar_chr(object$method), "` is not a method cpt_detect() ",
                "knows. Supply the interval yourself via as_ggcpt(ci = ), or ",
                "register the detector with cpt_register_method() so it can be ",
                "re-run.", class = "capability_absent")
    }
  }

  # Outside the `auto` block on purpose: `auto` would have picked "native"
  # for the engine this matters to (strucchange supplies its own intervals),
  # so the reachable call is the EXPLICIT `method = "bootstrap"` the review
  # named, which skipped every check in that block.
  if (method == "bootstrap" && !rerun_matches_result(object)) {
    cpt_abort("This result reports `change_in = \"",
              scalar_chr(object$change_in), "\"`, which `",
              scalar_chr(object$method), "` cannot be asked for through ",
               "cpt_detect(), so the bootstrap ", "would re-run a different ",
               "model than the one that produced it and ", "report the spread ",
               "of the wrong search as this result's interval. ",
              "A formula fit keeps neither the formula nor `data` on the object. ",
              "Use the engine's own intervals with `method = \"native\"`, or ",
              "supply the interval yourself via as_ggcpt(ci = ).",
              class = "capability_absent")
  }

  if (method %in% c("bootstrap", "nsp")) {
    # Both routes re-derive the interval from `$data$value`; see rerun_dots().
    if (n_coordinates(object) > 1L) {
      rerun_dots(object, list(), paste0("cpt_confint(method = \"", method,
                                        "\")"))
    }
  }

  out <- switch(method,
    native = {
      if (!has_native) {
        cpt_abort("This result carries no engine confidence intervals. Engines ",
                  "that supply them: ",
                  paste(subset(cpt_methods(), ci %in% TRUE)$method,
                        collapse = ", "),
                  ". Use method = \"bootstrap\" for a model-agnostic interval.",
                  class = "capability_absent")
      }
      tibble::tibble(cp = cp,
                     ci_lower = as.integer(native$lower),
                     ci_upper = as.integer(native$upper),
                     level = native$level, source = native$source)
    },
    posterior = confint_posterior(object, level),
    conformal = confint_conformal(object, level, B = B, seed = seed),
    bootstrap = confint_bootstrap(object, level, B = B, seed = seed, ...),
    nsp = confint_nsp(object, level, seed = seed, ...)
  )

  # Never report an interval that does not contain the point it belongs to,
  # and never one that runs off the series.
  n <- nrow(object$data)
  out$ci_lower <- pmax(1L, pmin(as.integer(out$ci_lower), out$cp))
  out$ci_upper <- pmin(n - 1L, pmax(as.integer(out$ci_upper), out$cp))

  if (!is.null(object$index)) {
    idx <- object$index
    out$cp_index <- idx[out$cp]
    out$ci_lower_index <- idx[out$ci_lower]
    out$ci_upper_index <- idx[out$ci_upper]
  }

  # A route that reports the engine's own level cannot honour the caller's,
  # and `method = "auto"` is where that bites: an `nsp` result carries
  # regions, so "auto" resolves to "native" and reports them at the level
  # nsp_wrapper()'s `alpha` produced -- 0.9 by default. Asking for
  # `level = 0.95` there returned 0.9 with nothing said about it, and the
  # only signal was the `level` column the caller would have to inspect.
  #
  # `?cpt_confint` documented that "native" ignores `level`, which is true
  # and is not the problem: the caller who is surprised asked for "auto".
  # So the check is on the answer rather than on the route -- if the level
  # reported differs from the one requested, say so, whichever branch got
  # there. cpt_monitor() already warns on this pattern for `arl0` and
  # friends; this is the same courtesy.
  if (level_supplied) {
    got <- unique(out$level[is.finite(out$level)])
    # Most engines' own intervals carry no recorded level at all (smuce,
    # strucchange, segmented, bfast and taylor report bounds, not the level
    # they were computed at), and the check below needs a level to compare,
    # so an explicit `level = 0.8` passed through without a word there --
    # the silence the paragraph above exists to prevent.
    if (length(got) == 0 && nrow(out) > 0 && identical(method, "native")) {
      cpt_warn("`level = ", format(requested), "` was not applied: `method = ",
               "\"native\"` reports the interval the engine already ",
               "computed, and this result does not record the level it was ",
               "computed at (see the engine's own `alpha` or `conf_level`). ",
               "Use `method = \"bootstrap\"` or `\"nsp\"` for an interval ",
               "computed at the level you ask for.",
               class = "level_not_applied")
    }
    if (length(got) > 0 && !any(abs(got - requested) < 1e-9)) {
      cpt_warn("`level = ", format(requested), "` was not applied: `method = ",
               "\"", method, "\"` reports the interval the engine already ",
               "computed, at level ", paste(format(got), collapse = "/"),
               ". Use `method = \"bootstrap\"` or `\"nsp\"` for an interval ",
               "computed at the level you ask for.",
               class = "level_not_applied")
    }
  }
  out
}

# Internal: the interval an engine supplied itself, in whatever shape it
# supplied it. Most engines fill `ci_lower`/`ci_upper`; NSP instead reports
# an interval that provably contains a change, under `region_start` and
# `region_end`. Reading only the first pair sent every NSP result down the
# bootstrap path -- 200 re-runs of the detector, to produce a weaker
# statement than the one already sitting on the object.
#' @noRd
native_bounds <- function(object) {
  cp_tbl <- object$changepoints
  if (all(c("ci_lower", "ci_upper") %in% names(cp_tbl))) {
    return(list(lower = cp_tbl$ci_lower, upper = cp_tbl$ci_upper,
                level = NA_real_, source = "native"))
  }
  if (all(c("region_start", "region_end") %in% names(cp_tbl))) {
    lvl <- object$region_level
    return(list(lower = cp_tbl$region_start, upper = cp_tbl$region_end,
                level = if (is.null(lvl)) NA_real_ else 1 - as.numeric(lvl),
                source = "nsp_region"))
  }
  NULL
}

# Internal: can this result's detector be re-run? The bootstrap route needs
# to call cpt_detect() by name, which rules out a result assembled by
# as_ggcpt() from somewhere the package cannot reach.
#' @noRd
bootstrap_possible <- function(object) {
  method <- scalar_chr(object$method)
  if (is.na(method)) return(FALSE)
  method %in% builtin_registry()$method || !is.null(registry_get(method))
}

# Internal: refuse, before any work, a result whose detector cannot be
# re-run. cpt_influence() and cpt_sensitivity() used to find out one
# perturbation at a time: every re-run failed, and the caller got an
# engine error or a grid of failed rows about "Unknown method = custom".
#' @noRd
refuse_unrerunnable <- function(object, what) {
  if (!bootstrap_possible(object)) {
    cpt_abort("`", what, "` re-runs the detector, and `",
              scalar_chr(object$method), "` is not a method cpt_detect() ",
              "knows: a result built with as_ggcpt() records someone else's ",
              "changepoints. Register the detector with cpt_register_method() ",
              "to make it re-runnable.", class = "capability_absent",
              data = list(method = scalar_chr(object$method)))
  }
  invisible(TRUE)
}

# Internal: can the detector be re-run to reproduce THIS result? Being in the
# registry is not enough. `strucchange_wrapper(y ~ x1 + x2, data = d)` returns
# `change_in = "regression"`, and neither the formula nor `data` is
# recoverable from the object -- `$data$value` is the response alone. So the
# bootstrap resampled the response and re-ran an INTERCEPT-ONLY breakpoint
# search, on a series whose changepoints came from a multi-regressor model,
# and reported the spread of the wrong search as the interval of the right
# one. bootstrap_possible() said TRUE, because the method is in the registry.
#
# The discriminator is the MODEL CONTEXT, not the registry: `"regression"` is
# the one label in cpt_change_in_levels() that describes a fit to covariates
# the result does not carry. (`"covariance"` and `"network"` results keep
# their full input in `$data_wide`, so those re-run correctly.) Keyed on the
# label rather than on `supports`, because `supports` answers a different
# question -- what you may ASK cpt_detect() for -- and widening it to make
# this label "legal" would claim a request route that does not exist.
#' @noRd
rerun_matches_result <- function(object) {
  !identical(scalar_chr(object$change_in), "regression")
}

# ---------------------------------------------------------------------------
# Re-running a finished result.
#
# Six tools take a finished `ggcpt` and call cpt_detect() again on (a
# perturbation of) its series: the bootstrap route here, cpt_influence(),
# cpt_leverage(), cpt_sensitivity() and cpt_select(). Each used to rebuild
# the request from `object$method` and `object$change_in` alone, and three
# things went wrong in all of them at once:
#
#   * the fit's PENALTY was dropped, so a `penalty = 40` pelt fit was
#     bootstrapped, perturbed and re-segmented under the default MBIC;
#   * `change_in` was passed back verbatim, but a result records what it
#     DETECTED, which is not always something cpt_detect() can be ASKED for:
#     kcp says "running mean", envcpt names its selected model ("trend"),
#     wbsts says "var" for its one request type, sn says "acf". Every re-run
#     then failed in match.arg(), and the failures were swallowed -- a
#     zero-width interval blamed on "the detector found no changepoints",
#     and an influence table in which every observation destroyed the
#     segmentation;
#   * a multivariate result was re-run on `$data$value`, its first
#     coordinate (or the cross-sectional mean), so the re-fit answered a
#     different question from the one that produced the changepoints.
#
# These helpers are the one place those three are decided.
# ---------------------------------------------------------------------------

# Internal: the `change_in` to request when re-running `object`, or NULL
# when the label it records cannot be turned back into a request.
#' @noRd
rerun_change_in <- function(object) {
  method <- scalar_chr(object$method)
  label <- scalar_chr(object$change_in)
  support <- method_change_in_support()[[method]]
  registered <- !is.null(registry_get(method))
  requestable <- function(v) {
    if (length(v) != 1L || is.na(v) || !v %in% cpt_change_in_levels()) {
      return(FALSE)
    }
    # cpt_detect() accepts "mean" for every built-in method and routes it to
    # the native change type; a registration accepts only what it declared.
    if (registered) return(v %in% support)
    identical(v, "mean") || is.null(support) || v %in% support
  }
  if (requestable(label)) return(label)
  # kcp labels a result by its running statistic, and the two statistics
  # cpt_detect() can ask for map back one to one.
  stripped <- sub("^running ", "", label)
  if (!identical(stripped, label) && requestable(stripped)) return(stripped)
  # An engine with a single request type records its own label for what
  # that request found (wbsts' "var", envcpt's selected model), so that one
  # request is what produced the result.
  if (length(support) == 1L && requestable(support)) return(support)
  NULL
}

# Internal: the `penalty` to request when re-running `object`, or NULL to
# leave the engine's default. Only a number the fit was given, or a name
# cpt_detect() itself understands, is forwarded: an engine-specific
# descriptor ("alpha", "threshold", "prob_threshold", ...) describes a tuning
# argument `penalty` does not reach.
#' @noRd
rerun_penalty <- function(object) {
  pen <- object$penalty
  if (!is.list(pen)) return(NULL)
  type <- scalar_chr(pen[["type", exact = TRUE]])
  value <- suppressWarnings(as.numeric(pen[["value", exact = TRUE]]))
  if (identical(type, "Manual") && length(value) == 1L && is.finite(value)) {
    return(value)
  }
  if (!is.na(type) && type %in% c("MBIC", "BIC", "SIC", "AIC",
                                  "Hannan-Quinn", "None", "sSIC", "MDL")) {
    return(type)
  }
  NULL
}

# Internal: how a fit handled missing values and what constrained it, as the
# cpt_detect() arguments that reproduce that. The re-running tools replayed
# the method, change type, penalty and family and none of these, so on a fit
# made with `na_action = "omit"` every bootstrap replicate failed on the
# gaps (and cpt_select() stopped in base R), and a `fixed` changepoint was
# re-estimated like any other. Positions are the fit's own; a tool that
# re-runs on a reversed, shortened or halved series maps them with
# map_rerun_positions().
#' @noRd
rerun_constraints <- function(object) {
  out <- list()
  if (!is.null(object$diagnostics$na_omitted)) {
    out$na_action <- "omit"
  } else if (anyNA(object$data$value)) {
    out$na_action <- "engine"
  }
  cons <- object$constraints
  if (length(cons$fixed)) out$fixed <- as.integer(cons$fixed)
  if (!is.null(cons$within) && length(cons$within)) {
    out$within <- unname(as.matrix(cons$within))
  }
  if (!is.null(cons$min_segment)) out$min_segment <- cons$min_segment
  if (!is.null(cons$min_effect)) out$min_effect <- cons$min_effect
  out
}

# Internal: move the positional constraints of a re-run request (`fixed`,
# `within`) onto a transformed series of `n_new` observations. `map` takes
# a changepoint position on the fit's series to one on the new series;
# `NA` drops it. A fixed changepoint that no longer fits is dropped, and a
# window is clamped to the new series (dropped if nothing of it is left).
#' @noRd
map_rerun_positions <- function(dots, map, n_new) {
  if (length(dots$fixed)) {
    f <- map(dots$fixed)
    f <- sort(unique(f[!is.na(f) & f >= 1 & f <= n_new - 1]))
    dots$fixed <- if (length(f)) as.integer(f) else NULL
  }
  if (!is.null(dots$within)) {
    w <- as.matrix(dots$within)
    rows <- lapply(seq_len(nrow(w)), function(i) {
      b <- sort(map(w[i, ]))
      if (length(b) != 2L || anyNA(b)) return(NULL)
      b <- c(max(b[1], 1), min(b[2], n_new - 1))
      if (b[1] > b[2]) NULL else b
    })
    rows <- Filter(Negate(is.null), rows)
    dots$within <- if (length(rows)) do.call(rbind, rows) else NULL
  }
  dots
}

# Internal: the request a re-run makes, with the caller's `...` winning over
# anything derived from the result (the precedence cpt_detect() gives `...`
# over its own derived arguments). `what` names the calling function in the
# errors.
#' @noRd
rerun_dots <- function(object, dots, what) {
  if (n_coordinates(object) > 1L) {
    cpt_abort("`", what, "` re-runs the detector on the result's series, but ",
              "this `", scalar_chr(object$method), "` result has ",
              n_coordinates(object), " coordinates and the re-run would see ",
               "only ", "`$data$value` (the first coordinate, or the ",
               "cross-sectional ", "mean): a different problem from the one ",
               "that produced these ", "changepoints. Re-run the multivariate ",
               "detector yourself, or use ", "an interval the engine supplied ",
               "(`cpt_confint(method = ", "\"native\")`).",
              class = "unsupported")
  }
  # A formula fit keeps neither the formula nor `data`, so a re-run would
  # search the response alone: an intercept-only breakpoint search standing
  # in for the regression that found these changepoints. cpt_influence(),
  # cpt_leverage(), cpt_sensitivity() and cpt_select() did exactly that
  # (a minute of re-fits, reported as the regression's influence), while
  # the tools that checked first refused.
  if (!rerun_matches_result(object)) {
    cpt_abort("`", what, "` re-runs the detector, and this result was ",
              "fitted from a formula: the covariates are not stored on it, ",
              "so a re-run would search the response alone, a different ",
              "model from the one that found these changepoints. Re-run the ",
              "original cpt_detect(y ~ ..., data = ) call on the perturbed ",
              "data yourself.", class = "capability_absent",
              data = list(method = scalar_chr(object$method)))
  }
  if (is.null(dots[["change_in", exact = TRUE]])) {
    ci <- rerun_change_in(object)
    if (is.null(ci)) {
      cpt_abort("`", what, "` re-runs `", scalar_chr(object$method), "`, but ",
                "this result records `change_in = \"",
                scalar_chr(object$change_in), "\"`, which cpt_detect() cannot ",
                "be asked for, so the re-run would fit a different model. Pass ",
                "`change_in` through `...`, together with the engine argument ",
                "that produced this result (for example `family`, ",
                "`running_stat` or `parameter`).", class = "unsupported")
    }
    dots$change_in <- ci
  }
  if (is.null(dots[["penalty", exact = TRUE]])) {
    dots$penalty <- rerun_penalty(object)
  }
  # ...and the family: a Poisson fit re-run without it would be re-run as a
  # Gaussian one, a different model answering a different question.
  if (is.null(dots[["family", exact = TRUE]]) && !is.null(object$family) &&
      is.null(registry_get(scalar_chr(object$method)))) {
    dots$family <- object$family
  }
  # ...and its missing-value handling and constraints.
  cons <- rerun_constraints(object)
  for (nm in setdiff(names(cons), names(dots))) dots[[nm]] <- cons[[nm]]
  dots
}

# Internal: refuse a detection request up front, for the tools that call
# cpt_detect() inside tryCatch() and would otherwise turn a typo into a
# result. cpt_select(x, method = "wbz") built a ladder in which every rung
# errored, reported K = 0 labelled "wbz", and blamed the method for
# "tuning itself"; cpt_label_error_curve() returned an all-NA curve, after
# which cpt_learn_penalty() said the labels were satisfied at every
# penalty. Resolves the name the way cpt_detect() does (a registered name
# first, then partial matching against the built-ins) and returns it.
#' @noRd
check_detect_request <- function(method, change_in = "mean") {
  if (!is.character(method) || length(method) != 1L || is.na(method) ||
      !nzchar(method)) {
    cpt_abort("`method` must be a single method name; see cpt_methods().",
              class = "bad_argument")
  }
  if (method %in% planned_methods()$method) {
    cpt_abort("`", method, "` is planned but not wired in this release: see ",
               "the ", "\"planned\" rows of `cpt_methods()`.",
              class = "planned_method")
  }
  entry <- registry_get(method)
  if (is.null(entry)) {
    hit <- pmatch(method, builtin_registry()$method)
    if (is.na(hit)) {
      cpt_abort("`", method, "` is not a method cpt_detect() knows. See ",
                "cpt_methods() for the available ones, or register a detector ",
                "with cpt_register_method().", class = "unknown_method")
    }
    method <- builtin_registry()$method[hit]
  }
  if (!is.character(change_in) || length(change_in) != 1L ||
      !change_in %in% cpt_change_in_levels()) {
    cpt_abort("`change_in` must be one of ",
              paste0("\"", cpt_change_in_levels(), "\"", collapse = ", "), ".",
              class = "bad_argument")
  }
  if (!is.null(entry)) {
    if (!change_in %in% entry$change_in) {
      cpt_abort("`change_in = \"", change_in, "\"` is not supported by the ",
                "registered method `", method, "`. Supported: ",
                paste(entry$change_in, collapse = ", "), ".",
                class = "unsupported")
    }
  } else {
    validate_method_change_in(method, change_in)
  }
  method
}

#' @noRd
empty_confint <- function(object) {
  out <- tibble::tibble(cp = integer(), ci_lower = integer(),
                        ci_upper = integer(), level = numeric(),
                        source = character())
  if (!is.null(object$index)) {
    out$cp_index <- object$index[0]
    out$ci_lower_index <- object$index[0]
    out$ci_upper_index <- object$index[0]
  }
  out
}

# Internal: credible interval from a posterior changepoint-probability
# profile. The mass is allocated within the window bounded by the
# neighbouring changepoints, so two nearby changepoints cannot claim the
# same probability twice, and the interval is the narrowest set of
# consecutive positions around `cp` holding `level` of that window's mass.
#' @noRd
confint_posterior <- function(object, level) {
  prob <- posterior_prob_profile(object)
  if (is.null(prob)) {
    cpt_abort("This result carries no posterior changepoint-probability ",
               "profile; ", "bcp and beast supply one. Use method = ",
               "\"bootstrap\" instead.", class = "capability_absent")
  }
  cp <- object$changepoints$cp
  n <- length(prob)
  bounds <- c(0L, cp, n)

  lo <- integer(length(cp))
  hi <- integer(length(cp))
  for (i in seq_along(cp)) {
    left <- if (i == 1L) 1L else as.integer(floor((cp[i - 1] + cp[i]) / 2)) + 1L
    right <- if (i == length(cp)) {
      n
    } else {
      as.integer(floor((cp[i] + cp[i + 1]) / 2))
    }
    win <- seq.int(max(1L, left), max(max(1L, left), right))
    p <- prob[win]
    total <- sum(p, na.rm = TRUE)
    if (!is.finite(total) || total <= 0) {
      lo[i] <- cp[i]; hi[i] <- cp[i]; next
    }
    # Grow outwards from the changepoint until `level` of the window's mass
    # is covered: a contiguous interval, which is what a reader expects a
    # location interval to be.
    centre <- match(cp[i], win)
    if (is.na(centre)) centre <- which.max(p)
    l <- centre; r <- centre
    covered <- p[centre]
    while (covered / total < level && (l > 1L || r < length(win))) {
      take_left <- if (l > 1L && r < length(win)) {
        p[l - 1L] >= p[r + 1L]
      } else {
        l > 1L
      }
      if (take_left) {
        l <- l - 1L; covered <- covered + p[l]
      } else {
        r <- r + 1L; covered <- covered + p[r]
      }
    }
    lo[i] <- win[l]; hi[i] <- win[r]
    # A wide interval here is not a wide uncertainty about the location, it
    # is a diffuse profile: measured, bcp and beast both put about
    # two-thirds of a window's posterior changepoint mass at the estimate
    # and the rest as a thin floor over every other position, so `level`
    # acts less like a confidence level than like a switch between the mode
    # (width 0 at 0.5) and most of the window (width 166 of 199 at 0.95).
    # The construction is right and the requested level is delivered -- what
    # is missing is any way for the reader to tell this apart from a bug,
    # standing next to a bootstrap interval two observations wide.
    span <- hi[i] - lo[i] + 1L
    if (span > length(win) / 2 && length(win) > 4L) {
      at_mode <- p[centre] / total
      cpt_warn("The posterior interval for the changepoint at ", cp[i],
               " covers ", span, " of the ", length(win), " positions in its ",
                "window. That is the profile, not the ", "location: ",
               format(100 * at_mode, digits = 2), "% of the window's ",
                "posterior changepoint mass sits at the ", "estimate itself ",
                "and the rest is spread thinly, so reaching ", format(level),
               " requires most of the window. Compare ",
               "`method = \"bootstrap\"`, and see the `\"posterior\"` note ",
               "in ?cpt_confint.", class = "wide_interval")
    }
  }
  tibble::tibble(cp = cp, ci_lower = lo, ci_upper = hi, level = level,
                 source = "posterior")
}

# Internal: the signal the resampling tools perturb around. Residuals are
# resampled within segments and added back to it, so it has to carry the
# change the fit found. The segment means do for a change in mean, variance
# or distribution, where anything constant within a segment cancels out of
# a within-segment resample. They do not for a change in slope or in
# seasonality: there the residuals against a segment mean still hold the
# trend or the seasonal cycle, and resampling them shuffled it away. On a
# clean slope change (cpop, a rise to 10 and back, noise sd 0.5) every
# bootstrap replicate came back flat, the interval collapsed onto the
# point with "the detector found no changepoints in 20 of 20 draws", and
# cpt_stability() scored the change 0. So those two use the engine's
# fitted signal, or a least-squares line per segment when a slope engine
# supplies none (not, nsp).
#' @noRd
resampling_signal <- function(object) {
  seg <- object$segments
  step <- rep(seg$param_estimate, times = seg$n)
  ci <- scalar_chr(object$change_in)
  if (!ci %in% c("slope", "seasonality")) return(step)
  v <- object$data$value
  eng <- object$data[["fitted"]]
  if (!is.null(eng) && length(eng) == length(v)) return(as.numeric(eng))
  if (!identical(ci, "slope")) return(step)
  segment_lines(v, seg$start, seg$end, step)
}

# Internal: a least-squares line through each segment `start[i]..end[i]` of
# `v`, on position, ignoring gaps; `fallback` where a segment has fewer
# than two observations to fit.
#' @noRd
segment_lines <- function(v, start, end, fallback) {
  out <- fallback
  for (i in seq_along(start)) {
    idx <- seq.int(start[i], end[i])
    ok <- idx[!is.na(v[idx])]
    if (length(ok) < 2L) next
    co <- stats::coef(stats::lm.fit(cbind(1, ok), v[ok]))
    if (all(is.finite(co))) out[idx] <- co[1] + co[2] * idx
  }
  out
}

# Internal: the fitted signal the residual diagnostics are taken against,
# or NULL for the segment means (residual_dependence()'s own default): the
# engine's signal where the result has one, and for a slope fit without
# one the line per segment above. Against a segment's level, a slope
# fit's residuals carry the trend, which the Ljung-Box test reads as
# dependent noise.
#' @noRd
residual_signal <- function(fit) {
  eng <- fit$data[["fitted"]]
  if (!is.null(eng) && length(eng) == nrow(fit$data)) return(as.numeric(eng))
  if (identical(scalar_chr(fit$change_in), "slope")) {
    return(resampling_signal(fit))
  }
  NULL
}

# Internal: bootstrap interval. Resamples residuals within the fitted
# segments (preserving the regime structure), re-runs the detector, and for
# each original changepoint records the nearest re-detection. Quantiles of
# that distribution are the interval.
#' @noRd
confint_bootstrap <- function(object, level, B = 200, seed = NULL, ...) {
  method <- object$method
  # Inherit the change type from the result, not just the method.
  #
  # Every entry point that takes raw data forwards `change_in` to
  # cpt_detect() explicitly; every one that takes a finished `ggcpt` read
  # `object$method` and stopped there, leaving `change_in` at its own
  # default of "mean" -- so a meanvar or var fit was silently re-detected as
  # a change in the MEAN. `object$change_in` was on the object the whole
  # time and never read.
  #
  # The symptom pointed away from the cause. On a pure variance change the
  # mean detector finds nothing, every replicate is discarded, and the
  # caller gets a zero-width interval plus a warning blaming the detector
  # for finding no changepoints.
  #
  # `...` still wins, so an explicit `change_in` overrides the object --
  # same precedence cpt_detect() gives `dots` over `derived_args_for()`.
  if (isTRUE(object$registered) && is.null(registry_get(method))) {
    cpt_abort("This result came from a registered method (`", method,
              "`) that is no longer registered, so it cannot be re-run for a ",
              "bootstrap interval.", class = "capability_absent")
  }
  if (!bootstrap_possible(object)) {
    cpt_abort("A bootstrap interval re-runs the detector, but `", method,
              "` is not a method cpt_detect() knows. Use method = \"native\" or ",
              "\"posterior\" if the engine supplied one, or register the ",
              "detector with cpt_register_method() so it can be re-run.",
              class = "capability_absent")
  }
  # ...and the penalty, and a change type cpt_detect() can be asked for;
  # see rerun_dots().
  refuse_resampling_family(object, "cpt_confint(method = \"bootstrap\")")
  dots <- rerun_dots(object, list(...), "cpt_confint(method = \"bootstrap\")")

  data_vec <- object$data$value
  n <- length(data_vec)
  cp <- object$changepoints$cp
  seg <- object$segments
  fitted_step <- resampling_signal(object)
  resid <- data_vec - fitted_step
  seg_id <- rep(seq_len(nrow(seg)), times = seg$n)

  local_seed(seed)

  draws <- matrix(NA_real_, nrow = B, ncol = length(cp))
  # A replicate the detector FAILED on is not one on which it found nothing,
  # and counting the two together blamed the detector for an error it
  # raised: the one diagnostic this bootstrap gave was the wrong one.
  n_failed <- 0L
  first_error <- NULL
  once_per_kind(for (b in seq_len(B)) {
    resampled <- resid
    for (s in seq_len(nrow(seg))) {
      # Observed positions only: a gap stays where it is, rather than
      # being resampled into the replicate as a missing residual.
      idx <- which(seg_id == s & !is.na(resid))
      # `sample.int()` on the index, not `sample()` on the values: R's
      # classic pitfall is that `sample(x, n)` means `sample.int(x, n)` when
      # `x` is a single number >= 1, so a one-observation segment resamples
      # `1:round(resid)` instead of the residual itself. It is currently
      # safe only by accident -- a length-1 segment's residual against its
      # own mean is exactly 0, and `0 >= 1` is FALSE -- which couples this
      # bootstrap to `param_estimate` staying the exact segment mean.
      resampled[idx] <- resid[idx][sample.int(length(idx), length(idx),
                                              replace = TRUE)]
    }
    rep_cp <- tryCatch(
      do.call(cpt_detect,
              c(list(fitted_step + resampled, method = method),
                dots))$changepoints$cp,
      error = function(e) {
        n_failed <<- n_failed + 1L
        if (is.null(first_error)) first_error <<- conditionMessage(e)
        NULL
      }
    )
    if (length(rep_cp) == 0) next
    # Nearest re-detection to each original changepoint. A replicate that
    # detects nothing contributes nothing rather than being counted as a
    # detection at the original location, which would shrink the interval.
    draws[b, ] <- vapply(cp, function(k) {
      rep_cp[which.min(abs(rep_cp - k))]
    }, numeric(1))
  })

  a <- (1 - level) / 2
  lo <- integer(length(cp))
  hi <- integer(length(cp))
  n_eff <- integer(length(cp))
  for (j in seq_along(cp)) {
    d <- draws[, j]
    d <- d[is.finite(d)]
    n_eff[j] <- length(d)
    if (length(d) == 0) {
      lo[j] <- cp[j]; hi[j] <- cp[j]; next
    }
    q <- stats::quantile(d, probs = c(a, 1 - a), names = FALSE, type = 1)
    lo[j] <- as.integer(floor(q[1]))
    hi[j] <- as.integer(ceiling(q[2]))
  }
  if (n_failed == B) {
    cpt_abort("The detector could not be re-run on any of the ", B,
              " bootstrap ", "replicates, so there is no interval to report. ",
               "The first error ", "was: ", first_error, class = "engine_error")
  }
  if (n_failed > 0L) {
    cpt_warn("The detector failed on ", n_failed, " of ", B, " bootstrap ",
             "replicate(s), which are excluded; the first error was: ",
             first_error, class = "replicates_failed")
  }
  n_empty <- sum(B - n_failed - n_eff)
  if (n_empty > 0) {
    cpt_warn("The detector found no changepoints in ", n_empty, " of ",
             (B - n_failed) * length(cp), " (replicate, changepoint) ",
             "draws; those draws are excluded, so the interval is based on ",
             "fewer than `B` replicates. `n_replicates` records how many were ",
             "used.", class = "replicates_failed")
  }

  tibble::tibble(cp = cp, ci_lower = lo, ci_upper = hi, level = level,
                 source = "bootstrap", n_replicates = n_eff)
}

# Internal: NSP regions matched to the changepoints. A changepoint inside a
# region takes that region; one outside every region gets NA rather than
# being handed the nearest, because "no region covers this" is a finding.
#' @noRd
confint_nsp <- function(object, level, seed = NULL, ...) {
  need_pkg("nsp")
  cp <- object$changepoints$cp
  alpha <- 1 - level
  fit <- nsp_wrapper(object$data$value, alpha = alpha, seed = seed, ...)
  reg <- fit$regions

  lo <- rep(NA_integer_, length(cp))
  hi <- rep(NA_integer_, length(cp))
  if (!is.null(reg) && nrow(reg) > 0) {
    for (i in seq_along(cp)) {
      inside <- which(reg$start <= cp[i] & reg$end >= cp[i])
      if (length(inside) > 0) {
        # The narrowest covering region is the sharpest statement available.
        j <- inside[which.min(reg$end[inside] - reg$start[inside])]
        lo[i] <- reg$start[j]
        hi[i] <- reg$end[j]
      }
    }
  }
  if (anyNA(lo)) {
    cpt_inform(sum(is.na(lo)), " of ", length(cp), " changepoint(s) fall in ",
                "no ", "NSP region at global level ", format(alpha),
               "; those rows are NA. A changepoint in no region is not ",
               "supported by NSP at this level.")
  }
  out <- tibble::tibble(cp = cp, ci_lower = lo, ci_upper = hi, level = level,
                        source = "nsp_region")
  # The clipping in cpt_confint() would turn an NA row into `cp`; keep the NA.
  out$ci_lower[is.na(lo)] <- NA_integer_
  out
}

#' Test detected changepoints
#'
#' Attaches a test to each detected changepoint (\code{type = "jump"}) or to
#' each fitted segment (\code{type = "segment"}), using the engine's own
#' test where it has one and an explicitly \emph{unadjusted} two-sample test
#' where it does not.
#'
#' @param object A \code{ggcpt} object.
#' @param type \code{"jump"} (one row per changepoint: is the change at this
#'   location real?) or \code{"segment"} (one row per segment: does this
#'   segment differ from the one before it?).
#' @param correction Multiple-testing correction applied across the rows, one
#'   of the \code{\link[stats]{p.adjust}} methods (\code{"none"},
#'   \code{"bonferroni"}, \code{"holm"}, \code{"BH"}, ...). Defaults to
#'   \code{"none"}; a \code{p_adjusted} column is added when it is not.
#' @param relevance Optional smallest change worth acting on, in noise
#'   standard deviations. The null hypothesis becomes "the change is no
#'   larger than \code{relevance}" rather than "there is no change", so a
#'   statistically detectable but practically negligible shift is not
#'   significant. Applies to the Welch route (the changepoints of a series
#'   fit); the \code{p_value} column is then the relevance test's, and a
#'   \code{relevance} column records the threshold in the data's units.
#'
#' @section Selection bias: read this before quoting a p-value:
#' Testing a changepoint at a location that was \emph{chosen because the data
#' looked like it changed there} is circular, and the resulting p-values are
#' anti-conservative, often severely. The \code{selection_adjusted} column
#' records, per row, whether the test accounts for that:
#' \itemize{
#'   \item \code{TRUE} for \pkg{segmented}'s Davies test, which is built
#'     for a nuisance parameter present only under the alternative. It is
#'     one \emph{global} test of "is there a breakpoint", not a test per
#'     breakpoint, so on a multi-break fit every row carries the same
#'     statistic and p-value, and the method string says so.
#'   \item \code{FALSE} for the generic Welch two-sample fallback, which
#'     compares the segments either side of the changepoint as if the
#'     location had been fixed in advance. Useful as a descriptive effect
#'     size with a scale attached; not a valid significance test for the
#'     existence of the change. The exception is a changepoint that
#'     \emph{was} fixed in advance, through \code{cpt_detect(fixed = )}:
#'     its row is \code{TRUE}.
#'   \item The fallback compares means whatever the fit detected, so on a
#'     fit of a change in variance, slope or distribution it tests a
#'     different change from the one found, and warns (class
#'     \code{ggchangepoint_assumption}).
#'   \item \code{FALSE} for \pkg{strucchange}'s route as well, which is
#'     the Chow F evaluated \emph{at} each estimated break date. The Chow
#'     statistic's reference distribution assumes the date was fixed in
#'     advance, so quoting it at a date the Bai-Perron dynamic program chose
#'     is exactly the circularity this column exists to flag. Reporting
#'     it is conventional in that literature, which does not make it
#'     adjusted. The selection-adjusted objects there are the sup-type
#'     statistics and the Bai-Perron critical values.
#' }
#' Two further limits worth knowing. \code{type = "segment"} is
#' \strong{always} the unadjusted Welch test: the native routes above apply
#' only to \code{type = "jump"}, so a \pkg{segmented} fit tested
#' per-segment does not use Davies' test. And the split above is by
#' \emph{engine and type}, not by engine alone.
#' For a guarantee that survives selection, use \code{\link{nsp_wrapper}()}
#' (regions with exact global coverage) or \code{\link{cpt_confint}()} with
#' \code{method = "nsp"}. The canonical post-detection tests of Jewell,
#' Fearnhead and Witten (2022) are implemented in \code{ChangepointInference},
#' which is not on CRAN; \code{\link{cpt_register_method}()} is the supported
#' way to bring it in.
#'
#' @return A tibble with columns \code{cp} (or \code{seg_id}),
#'   \code{estimate}, \code{statistic}, \code{p_value}, \code{method} and
#'   \code{selection_adjusted}, plus \code{cp_index} on the original scale
#'   when the result carries a time index (as \code{\link{cpt_confint}()}
#'   does), and \code{p_adjusted} when \code{correction} is not
#'   \code{"none"}. \code{estimate} is the change in the segment mean
#'   (the breakpoint's slope change for \pkg{segmented}), and \code{NA}
#'   for a regression-mode \pkg{strucchange} fit, whose change is a vector
#'   of coefficients; its Chow test is run on the full regression model.
#'   A multivariate result is refused: both routes test one series.
#' @seealso \code{\link{cpt_confint}()}, \code{\link{nsp_wrapper}()}.
#' @export
#' @examples
#' set.seed(2026)
#' fit <- cpt_detect(c(rnorm(60), rnorm(60, 4)), method = "pelt")
#' cpt_test(fit)
cpt_test <- function(object, type = c("jump", "segment"),
                     correction = "none", relevance = NULL) {
  if (!is_ggcpt(object)) {
    cpt_abort("`object` must be a ggcpt object.", class = "bad_argument")
  }
  type <- cpt_match_arg(type)
  correction <- cpt_match_arg(correction, stats::p.adjust.methods)
  if (!is.null(relevance)) validate_scalar(relevance, "relevance", min = 0)
  # Every route below tests `$data$value`, the result's one univariate
  # series; for a multivariate result that is its first coordinate (or the
  # cross-sectional mean), which need not carry the change at all.
  if (n_coordinates(object) > 1L) {
    cpt_abort("`cpt_test()` compares the segments either side of each ",
              "changepoint on a single series, but this `",
              scalar_chr(object$method), "` result has ", n_coordinates(object),
              " coordinates and the test would see only `$data$value` (the ",
              "first coordinate, or the cross-sectional mean). Test a ",
              "coordinate you care about with cpt_test(as_ggcpt(",
              "fit$changepoints$cp, <that column>)).",
              class = "wrong_dimension")
  }

  out <- if (type == "segment") {
    test_segments(object)
  } else if (!is.null(relevance)) {
    naive_jump_test(object)
  } else {
    native <- native_jump_test(object)
    native %||% regression_jump_test(object) %||% naive_jump_test(object)
  }
  if (!is.null(relevance) && nrow(out) > 0) {
    out <- apply_relevance(out, object, relevance)
  }

  # On the original scale too, when there is one. `cpt_confint()` has
  # reported `cp_index` since the time index was added; this table is the
  # other half of the same question about the same changepoints and reported
  # positions only, so a user with dates got a date from one function and an
  # integer from the other.
  if (!is.null(object$index) && "cp" %in% names(out)) {
    # No `nrow(out) > 0` guard: `index[integer(0)]` is a zero-length value
    # of the index's own type, which is what keeps the column set the same
    # whether or not anything was found -- `empty_confint()` does the same,
    # and it is what makes rbind() over several fits line up.
    out$cp_index <- object$index[out$cp]
  }
  if (nrow(out) > 0 && correction != "none") {
    out$p_adjusted <- stats::p.adjust(out$p_value, method = correction)
    attr(out, "correction") <- correction
  }
  # The fallback compares segment MEANS, whatever the fit detected: on a
  # clean variance change it reported p = 0.31, and on a slope change (a
  # rise and a fall with equal means) p = 0.90, each beside a changepoint
  # the detector had no doubt about.
  welch <- nrow(out) > 0 && any(grepl("^Welch", out$method))
  ci_fit <- scalar_chr(object$change_in)
  if (welch && !is.na(ci_fit) && !ci_fit %in% c("mean", "regression")) {
    cpt_warn("`cpt_test()` compares the segment means either side of each ",
             "changepoint (a Welch t-test), but this fit detected a change in ",
             "`", ci_fit, "`. The test measures a shift in the mean only, so ",
             "a large p-value does not argue against the change the detector ",
             "found.", class = "assumption", data = list(change_in = ci_fit))
  }
  if (nrow(out) > 0 && any(!out$selection_adjusted)) {
    cpt_warn("`selection_adjusted` is FALSE for ", sum(!out$selection_adjusted),
             " of ", nrow(out), " row(s): the ", "changepoint locations were ",
              "chosen from these data, so those ", "p-values are ",
              "anti-conservative. See the selection-bias section ",
             "of ?cpt_test.", class = "selection")
  }
  out
}

# Internal: the engine's own test, when it has one. NULL otherwise.
#' @noRd
native_jump_test <- function(object) {
  fit <- object$fit
  cp <- object$changepoints$cp
  if (length(cp) == 0 || is.null(fit)) return(NULL)

  if (inherits(fit, "breakpoints") &&
      requireNamespace("strucchange", quietly = TRUE)) {
    return(strucchange_jump_test(object, fit))
  }
  if (inherits(fit, "segmented") &&
      requireNamespace("segmented", quietly = TRUE)) {
    return(segmented_jump_test(object, fit))
  }
  NULL
}

# Internal: a Chow F test at each estimated break date. `sctest()` with
# `type = "Chow"` is strucchange's own single-break test; running it at the
# dates the Bai-Perron dynamic program selected is the standard reporting
# convention in that literature -- but the convention does not make the
# p-value selection-adjusted, and this column used to claim it was. The
# Chow F's reference distribution assumes the break date was fixed IN
# ADVANCE; evaluating it at a date chosen because the data looked like it
# broke there is precisely the circularity `selection_adjusted` exists to
# flag. The selection-adjusted objects in this framework are the sup-type
# statistics (supF / Fstats) and the Bai-Perron critical values, not a
# pointwise Chow F, so the flag is FALSE and the method string says
# "unadjusted" out loud.
#' @noRd
strucchange_jump_test <- function(object, fit) {
  cp <- object$changepoints$cp
  y <- object$data$value
  d <- data.frame(.y = y)
  # A formula fit's breaks are in its regression COEFFICIENTS, and the test
  # has to be of that model. It used to be `.y ~ 1` whatever the fit, so a
  # break in a slope from 2 to -1 (Chow F = 772 on the fitted model) was
  # reported as F = 0.26, p = 0.61: an intercept-only test of the response,
  # i.e. no evidence at all for the break the method had just dated. The
  # breakpoints object keeps its own response and design matrix.
  X <- fit[["X", exact = TRUE]]
  yy <- fit[["y", exact = TRUE]]
  full <- is.matrix(X) && ncol(X) > 1L && length(yy) == nrow(X)
  rows <- lapply(cp, function(k) {
    tst <- tryCatch(
      if (full) {
        strucchange::sctest(yy ~ X - 1, type = "Chow", point = k)
      } else {
        strucchange::sctest(.y ~ 1, data = d, type = "Chow", point = k)
      },
      error = function(e) NULL
    )
    if (is.null(tst)) return(NULL)
    left <- y[seq_len(k)]
    right <- y[seq.int(k + 1L, length(y))]
    tibble::tibble(cp = k,
                   estimate = if (full) NA_real_ else mean(right) - mean(left),
                   statistic = as.numeric(tst$statistic),
                   p_value = as.numeric(tst$p.value),
                   method = if (full) {
                     "Chow F at estimated break, full regression model (unadjusted)"
                   } else {
                     "Chow F at estimated break (unadjusted)"
                   },
                   selection_adjusted = FALSE)
  })
  rows <- Filter(Negate(is.null), rows)
  if (length(rows) == 0) return(NULL)
  do.call(rbind, rows)
}

# Internal: segmented's own inference. The Davies test addresses exactly the
# nuisance-parameter-under-the-alternative problem that makes a naive test
# invalid, so it is the selection-adjusted route; the per-break estimate and
# statistic come from the U-terms of the fitted model.
#' @noRd
segmented_jump_test <- function(object, fit) {
  cp <- object$changepoints$cp
  # Davies' test asks whether the BASE model (the line before any
  # breakpoint) needs a breakpoint. Handed the segmented fit itself it asks
  # something else, whether one MORE breakpoint is needed, which on a
  # correctly fitted series is "no": a p-value near 1 labelled as the test
  # of the breakpoint that was found.
  base <- attr(fit, "ggcpt_base")
  seg_var <- object$regression$seg_z %||% ".t"
  if (is.null(base)) {
    base <- stats::lm(.y ~ .t,
                      data = data.frame(.y = object$data$value,
                                        .t = seq_len(nrow(object$data))))
    seg_var <- ".t"
  }
  dav <- tryCatch(
    segmented::davies.test(base,
                           seg.Z = stats::as.formula(paste("~", seg_var))),
    error = function(e) NULL
  )
  # A Davies test that could not be computed must not become a row that
  # claims a selection-adjusted p-value of NA: fall back to the honest,
  # unadjusted route instead.
  if (is.null(dav) || !is.finite(as.numeric(dav$p.value))) return(NULL)
  tt <- tryCatch(summary(fit)$Ttable, error = function(e) NULL)
  u_rows <- if (!is.null(tt)) grep("^U[0-9]+", rownames(tt)) else integer(0)

  tibble::tibble(
    cp = cp,
    estimate = if (length(u_rows) == length(cp)) {
      as.numeric(tt[u_rows, 1])
    } else {
      rep(NA_real_, length(cp))
    },
    statistic = if (!is.null(dav)) {
      rep(as.numeric(dav$statistic), length(cp))
    } else {
      rep(NA_real_, length(cp))
    },
    p_value = if (!is.null(dav)) {
      rep(as.numeric(dav$p.value), length(cp))
    } else {
      rep(NA_real_, length(cp))
    },
    # Davies' test is a single test of "is there A breakpoint", not a test
    # per breakpoint, so on a three-break fit these rows carry the same
    # statistic and p-value three times -- which reads as three
    # per-changepoint tests. It IS the right selection-adjusted object; the
    # row-wise presentation was the problem, so the method string names the
    # scope.
    method = "Davies test (segmented; global, one test per fit)",
    selection_adjusted = TRUE
  )
}

# Internal: a Chow F test at each break of a regression fit whose design
# the result carries (a formula fit through fastcpd), on the rows between
# the neighbouring breaks. Unadjusted, like every pointwise test at a
# data-chosen date. NULL for a result that is not a regression.
#' @noRd
regression_jump_test <- function(object) {
  X <- object$regression$X
  cp <- object$changepoints$cp
  if (is.null(X) || !length(cp)) return(NULL)
  y <- object$data$value
  cp_fit <- cp
  if (!is.null(object$na_map)) {
    y <- y[object$na_map]
    cp_fit <- match(cp, object$na_map)
  }
  n <- length(y)
  k <- ncol(X)
  bounds <- c(0L, cp_fit, n)
  rss <- function(r) {
    sum(stats::lm.fit(X[r, , drop = FALSE], y[r])$residuals^2)
  }
  rows <- lapply(seq_along(cp_fit), function(i) {
    left <- (bounds[i] + 1L):bounds[i + 1L]
    right <- (bounds[i + 1L] + 1L):bounds[i + 2L]
    f <- NA_real_
    pv <- NA_real_
    if (length(left) > k && length(right) > k) {
      r1 <- rss(left) + rss(right)
      r0 <- rss(c(left, right))
      dfd <- length(left) + length(right) - 2 * k
      f <- ((r0 - r1) / k) / (r1 / dfd)
      pv <- stats::pf(f, k, dfd, lower.tail = FALSE)
    }
    tibble::tibble(
      cp = cp[i], estimate = NA_real_, statistic = f, p_value = pv,
      method = "Chow F at estimated break, full regression model (unadjusted)",
      selection_adjusted = FALSE)
  })
  do.call(rbind, rows)
}

# Internal: replace each Welch row's p-value with the test of a relevance
# hypothesis, H0: |change| <= relevance * sigma, one-sided on |t|.
#' @noRd
apply_relevance <- function(out, object, relevance) {
  y <- object$data$value
  thr <- relevance * noise_sd(y)
  cp <- object$changepoints$cp
  bounds <- c(0L, cp, length(y))
  out$p_value <- vapply(seq_len(nrow(out)), function(i) {
    k <- match(out$cp[i], cp)
    if (is.na(k)) return(NA_real_)
    a <- y[(bounds[k] + 1L):bounds[k + 1L]]
    b <- y[(bounds[k + 1L] + 1L):bounds[k + 2L]]
    if (length(a) < 2L || length(b) < 2L) return(NA_real_)
    va <- stats::var(a) / length(a)
    vb <- stats::var(b) / length(b)
    se <- sqrt(va + vb)
    dfw <- se^4 / (va^2 / (length(a) - 1) + vb^2 / (length(b) - 1))
    stats::pt((abs(mean(b) - mean(a)) - thr) / se, dfw, lower.tail = FALSE)
  }, numeric(1))
  out$relevance <- thr
  out$method <- paste0(out$method, ", relevance test (H0: |change| <= ",
                       format(signif(thr, 3)), ")")
  out
}

# Internal: the generic fallback. A Welch two-sample t-test across each
# changepoint, comparing the stretches bounded by the neighbouring
# changepoints. Honest about what it is: `selection_adjusted = FALSE`.
#' @noRd
naive_jump_test <- function(object) {
  cp <- object$changepoints$cp
  y <- object$data$value
  n <- length(y)
  if (length(cp) == 0) {
    return(tibble::tibble(cp = integer(), estimate = numeric(),
                          statistic = numeric(), p_value = numeric(),
                          method = character(),
                          selection_adjusted = logical()))
  }
  bounds <- c(0L, cp, n)
  fixed <- fixed_changepoints(object)
  rows <- lapply(seq_along(cp), function(i) {
    left <- y[seq.int(bounds[i] + 1L, bounds[i + 1L])]
    right <- y[seq.int(bounds[i + 1L] + 1L, bounds[i + 2L])]
    tt <- tryCatch(stats::t.test(right, left), error = function(e) NULL)
    # A location fixed in advance was not chosen from these data, which is
    # the one thing `selection_adjusted` records.
    pre <- cp[i] %in% fixed
    tibble::tibble(
      cp = cp[i],
      # na.rm: t.test() drops a fit's gaps itself, and the estimate beside
      # its statistic came back NA for every `na_action = "omit"` fit.
      estimate = mean(right, na.rm = TRUE) - mean(left, na.rm = TRUE),
      statistic = if (is.null(tt)) NA_real_ else as.numeric(tt$statistic),
      p_value = if (is.null(tt)) NA_real_ else as.numeric(tt$p.value),
      method = if (pre) "Welch two-sample t (location fixed in advance)" else
        "Welch two-sample t (unadjusted)",
      selection_adjusted = pre
    )
  })
  do.call(rbind, rows)
}

# Internal: the changepoints a fit was given through `fixed`, as positions.
#' @noRd
fixed_changepoints <- function(object) {
  cp <- object$changepoints
  if (!"fixed" %in% names(cp)) return(integer(0))
  as.integer(cp$cp[cp$fixed %in% TRUE])
}

# Internal: one row per segment, comparing it with the segment before.
#' @noRd
test_segments <- function(object) {
  seg <- object$segments
  y <- object$data$value
  if (nrow(seg) < 2) {
    return(tibble::tibble(seg_id = integer(), estimate = numeric(),
                          statistic = numeric(), p_value = numeric(),
                          method = character(),
                          selection_adjusted = logical()))
  }
  fixed <- fixed_changepoints(object)
  rows <- lapply(seq.int(2L, nrow(seg)), function(i) {
    prev <- y[seq.int(seg$start[i - 1], seg$end[i - 1])]
    cur <- y[seq.int(seg$start[i], seg$end[i])]
    tt <- tryCatch(stats::t.test(cur, prev), error = function(e) NULL)
    pre <- seg$end[i - 1] %in% fixed
    tibble::tibble(
      seg_id = seg$seg_id[i],
      estimate = mean(cur, na.rm = TRUE) - mean(prev, na.rm = TRUE),
      statistic = if (is.null(tt)) NA_real_ else as.numeric(tt$statistic),
      p_value = if (is.null(tt)) NA_real_ else as.numeric(tt$p.value),
      method = if (pre) "Welch two-sample t (location fixed in advance)" else
        "Welch two-sample t (unadjusted)",
      selection_adjusted = pre
    )
  })
  do.call(rbind, rows)
}