R/crops.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

#' CROPS: the full penalty path of a penalised changepoint method
#'
#' Runs PELT once per \emph{distinct} optimal segmentation as the penalty
#' ranges over \code{[pen_min, pen_max]}, using the CROPS algorithm of
#' Haynes, Eckley and Fearnhead (2017) as implemented by the
#' \pkg{changepoint} package. Instead of committing to one penalty, the
#' analyst sees every segmentation the data admits along the path, together
#' with its cost, and picks the elbow.
#'
#' @param x For \code{cpt_crops()}, a numeric vector, or a \code{ggcpt}
#'   fit whose series is used (with its change type, when it is one CROPS
#'   sweeps); for the \code{print()} and \code{tidy()} methods, a
#'   \code{ggcpt_path} object.
#' @param change_in What to detect change in: \code{"mean"}, \code{"var"},
#'   or \code{"meanvar"}. Defaults to \code{"mean"}.
#' @param pen_min,pen_max The penalty interval to sweep. Default to
#'   \code{log(n)} and \code{10 * log(n)}.
#' @param ... Additional arguments passed to the underlying
#'   \code{changepoint::cpt.mean()}, \code{cpt.var()}, or
#'   \code{cpt.meanvar()} call.
#' @return A \code{ggcpt_path} object: a list with a \code{solutions} tibble
#'   (one row per distinct segmentation: \code{penalty}, \code{n_cpts},
#'   \code{cost}, and a \code{cpts} list-column), the \code{data}, and
#'   metadata. \code{penalty} is the \strong{lower end} of the penalty
#'   interval on which that segmentation is optimal: CROPS returns the
#'   breakpoints of the penalty axis, so \eqn{K} segmentations come with
#'   \eqn{K + 1} boundaries and each row is optimal from its own
#'   \code{penalty} up to the next row's. Methods: \code{print()}, \code{tidy()}, and
#'   \code{autoplot()} (elbow plot by default;
#'   \code{type = "path"} for penalty vs. number of changepoints;
#'   \code{type = "segmentations"} for the faceted segmentations).
#' @references
#' \insertRef{haynes2017crops}{ggchangepoint}
#'
#' \insertRef{killick2014changepoint}{ggchangepoint}
#' @export
#' @examples
#' set.seed(2026)
#' x <- c(rnorm(100), rnorm(100, 3), rnorm(100, -1))
#' path <- cpt_crops(x)
#' path
#' ggplot2::autoplot(path)
#' ggplot2::autoplot(path, type = "segmentations")
cpt_crops <- function(x, change_in = c("mean", "var", "meanvar"),
                      pen_min = NULL, pen_max = NULL, ...) {
  # `...` reaches changepoint::cpt.mean()/cpt.var()/cpt.meanvar(), and this
  # function pins three of their arguments -- the two that make the call
  # CROPS at all, and the interval it sweeps. Passing any of them gave R's
  # raw "formal argument \"penalty\" matched by multiple actual arguments",
  # which names neither this function nor what to use instead. Same
  # treatment as the fifty wrappers get.
  reject_managed_args(list(...), "crops", c(
    penalty = paste("CROPS *is* the penalty regime this function runs; a",
                    "single penalty value is what `cpt_detect()` takes"),
    method = paste("the penalty path is computed by PELT; the other search",
                   "methods do not produce one"),
    pen.value = "use `pen_min` and `pen_max` to set the interval to sweep"))
  obs_pos <- NULL
  if (is_ggcpt(x)) {
    # The series a fit carries, as every other post-detection tool reads
    # it; the fit's own change type when the caller did not choose one.
    if (n_coordinates(x) > 1L) {
      cpt_abort("`cpt_crops()` sweeps one series, and this fit has ",
                n_coordinates(x), " coordinates.", class = "wrong_dimension")
    }
    if (missing(change_in) &&
        scalar_chr(x$change_in) %in% c("mean", "var", "meanvar")) {
      change_in <- scalar_chr(x$change_in)
    }
    x <- x$data$value
    # A fit made with `na_action = "omit"` keeps its gaps: the sweep runs on
    # the observed values, and its segmentations come back in the original
    # positions. It used to stop at "`x` must be finite".
    if (anyNA(x)) {
      full_series <- x
      obs_pos <- which(!is.na(x))
      x <- x[obs_pos]
    }
  }
  change_in <- cpt_match_arg(change_in)

  validate_data(x)
  data_vec <- as_uni_vector(x, "crops")
  n <- length(data_vec)

  if (is.null(pen_min)) pen_min <- log(n)
  if (is.null(pen_max)) pen_max <- 10 * log(n)
  validate_scalar(pen_min, "pen_min", min = 0, min_open = TRUE)
  validate_scalar(pen_max, "pen_max", min = 0, min_open = TRUE)
  if (pen_min >= pen_max) {
    cpt_abort("`pen_min` must be strictly smaller than `pen_max`.",
              class = "bad_argument")
  }

  cpt_fun <- switch(change_in,
    mean = changepoint::cpt.mean,
    var = changepoint::cpt.var,
    meanvar = changepoint::cpt.meanvar
  )

  # The CROPS run prints progress; keep the console clean.
  utils::capture.output(
    fit <- cpt_fun(data_vec, method = "PELT", penalty = "CROPS",
                   pen.value = c(pen_min, pen_max), ...)
  )

  pens <- as.numeric(changepoint::pen.value.full(fit))
  cpts_mat <- changepoint::cpts.full(fit)
  if (is.null(dim(cpts_mat))) {
    cpts_mat <- matrix(cpts_mat, nrow = 1)
  }

  cpts_list <- lapply(seq_len(nrow(cpts_mat)), function(i) {
    v <- cpts_mat[i, ]
    v <- as.integer(v[!is.na(v)])
    v[v >= 1 & v < n]
  })
  # `pen.value.full()` returns the penalty-axis BREAKPOINTS, not one value
  # per segmentation: K + 1 boundaries delimiting K intervals of penalty on
  # each of which one segmentation is optimal. Measured on a 240-point
  # series swept over [2, 200]: length(pens) == 23 against
  # nrow(cpts.full()) == 22, and pens[1] is the lower end of the swept
  # range. So row i is optimal on [pens[i], pens[i + 1]) and its label is
  # pens[i] -- which is what the old length-reconciliation happened to
  # keep, by truncating the tail, without saying why. Derive it instead,
  # and keep the reconciliation only as a backstop in case a later
  # changepoint release changes the convention.
  if (length(pens) == length(cpts_list) + 1L) {
    pens <- pens[seq_along(cpts_list)]
  } else if (length(pens) > length(cpts_list)) {
    pens <- pens[seq_along(cpts_list)]
  } else if (length(pens) < length(cpts_list)) {
    pens <- c(pens, rep(NA_real_, length(cpts_list) - length(pens)))
  }

  # Match the changepoint package's cost conventions so the elbow is the
  # same objective CROPS optimised: change-in-var fixes the mean at the
  # global estimate (as cpt.var does); meanvar re-estimates both.
  global_mean <- mean(data_vec)
  seg_cost <- function(cps) {
    breaks <- c(0L, cps, n)
    cost <- 0
    for (i in seq_len(length(breaks) - 1)) {
      seg <- data_vec[(breaks[i] + 1):breaks[i + 1]]
      m <- length(seg)
      if (change_in == "mean") {
        cost <- cost + sum((seg - mean(seg))^2)
      } else {
        mu <- if (change_in == "var") global_mean else mean(seg)
        v <- sum((seg - mu)^2) / m
        if (!is.finite(v) || v <= 0) v <- .Machine$double.eps
        cost <- cost + m * log(v)
      }
    }
    cost
  }
  costs <- vapply(cpts_list, seg_cost, numeric(1))

  solutions <- tibble::tibble(
    penalty = pens,
    n_cpts = vapply(cpts_list, length, integer(1)),
    cost = costs,
    cpts = cpts_list
  )
  solutions <- solutions[order(solutions$n_cpts), , drop = FALSE]
  data_tbl <- tibble::tibble(index = seq_len(n), value = data_vec)
  if (!is.null(obs_pos)) {
    solutions$cpts <- lapply(solutions$cpts, function(cp) obs_pos[cp])
    data_tbl <- tibble::tibble(index = seq_along(full_series),
                               value = full_series)
  }

  structure(
    list(
      solutions = solutions,
      data = data_tbl,
      change_in = change_in,
      pen_range = c(pen_min, pen_max),
      fit = fit,
      call = match.call()
    ),
    class = "ggcpt_path"
  )
}

#' @rdname cpt_crops
#' @param object A \code{ggcpt_path} object (for \code{autoplot()}).
#' @param type Plot type for \code{autoplot()}: \code{"elbow"} (cost against
#'   number of changepoints, the classic CROPS diagnostic), \code{"path"}
#'   (number of changepoints against penalty), or \code{"segmentations"}
#'   (the data faceted by solution, with that solution's changepoints
#'   drawn).
#' @param max_facets Maximum number of solutions shown by
#'   \code{type = "segmentations"}. Defaults to \code{12}.
#' @export
autoplot.ggcpt_path <- function(object,
                                type = c("elbow", "path", "segmentations"),
                                max_facets = 12, ...) {
  type <- cpt_match_arg(type)
  sol <- object$solutions

  if (type == "elbow") {
    return(
      ggplot2::ggplot(sol, ggplot2::aes(n_cpts, cost)) +
        ggplot2::geom_line(color = "grey50") +
        ggplot2::geom_point(size = 2) +
        ggplot2::labs(x = "Number of changepoints",
                      y = "Segmentation cost",
                      title = "CROPS penalty path: cost elbow")
    )
  }

  if (type == "path") {
    sol_p <- sol[is.finite(sol$penalty), , drop = FALSE]
    return(
      ggplot2::ggplot(sol_p, ggplot2::aes(penalty, n_cpts)) +
        ggplot2::geom_step(direction = "vh") +
        ggplot2::geom_point(size = 2) +
        ggplot2::labs(x = "Penalty", y = "Number of changepoints",
                      title = "CROPS penalty path")
    )
  }

  # type == "segmentations"
  sol_f <- utils::head(sol, max_facets)
  panels <- do.call(rbind, lapply(seq_len(nrow(sol_f)), function(i) {
    d <- object$data
    d$panel <- paste0(sol_f$n_cpts[i], " changepoint(s)")
    d$.order <- i
    d
  }))
  cp_panels <- do.call(rbind, lapply(seq_len(nrow(sol_f)), function(i) {
    cps <- sol_f$cpts[[i]]
    if (length(cps) == 0) return(NULL)
    data.frame(panel = paste0(sol_f$n_cpts[i], " changepoint(s)"),
               cp = cps, .order = i)
  }))
  panels$panel <- stats::reorder(panels$panel, panels$.order)

  p <- ggplot2::ggplot(panels, ggplot2::aes(index, value)) +
    ggplot2::geom_line(color = "grey50") +
    ggplot2::facet_wrap(~panel) +
    ggplot2::labs(x = "Index", y = "Value",
                  title = "CROPS penalty path: candidate segmentations")

  if (!is.null(cp_panels)) {
    cp_panels$panel <- factor(cp_panels$panel, levels = levels(panels$panel))
    p <- p + ggplot2::geom_vline(
      data = cp_panels,
      ggplot2::aes(xintercept = cp),
      color = "blue", linewidth = 0.4
    )
  }
  p
}

#' @rdname cpt_crops
#' @export
print.ggcpt_path <- function(x, ...) {
  cat("ggcpt_path (CROPS penalty path)\n")
  # Wider than the default: "Distinct segmentations" outgrows it, and a
  # label that outgrows the pad is what left this header ragged before.
  w <- 23L
  cat_field("Change in", x$change_in, w)
  cat_field("Penalty range",
            paste0("[", format(x$pen_range[1], digits = 4), ", ",
                   format(x$pen_range[2], digits = 4), "]"), w)
  cat_field("Series length", nrow(x$data), w)
  cat_field("Distinct segmentations", nrow(x$solutions), w)
  cat("\n")
  print(x$solutions[, c("penalty", "n_cpts", "cost")],
        n = min(nrow(x$solutions), 10))
  invisible(x)
}

#' @rdname cpt_crops
#' @export
tidy.ggcpt_path <- function(x, ...) {
  x$solutions
}