R/ecp.R

In ggchangepoint: Combines Changepoint Analysis with 'ggplot2'

#' ecp wrapper
#'
#' The ecp package provides a non-parametric way to detect changepoints. Unlike
#' the changepoint package, it does not assume raw data to have any formal
#' distribution. This wrapper function wraps two functions from the ecp package,
#' i.e., \code{e.divisive()} and \code{e.agglo()}. Users can use either function
#' by switching the \code{algorithm} argument. Before using the wrapper function,
#' seed should be set for the sake of reproducibility.
#'
#' @param data A numeric vector (for univariate) or matrix/data.frame (for
#'   multivariate).
#' @param algorithm Either \code{divisive} or \code{agglo}. \code{divisive} is
#'   the default.
#' @param min_size Minimum number of observations between change points. By
#'   default is 2. This argument is only applied when \code{algorithm =
#'   "divisive"}.
#' @param seed Optional. A seed for reproducibility of the stochastic
#'   permutation test. 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 ... Extra arguments to pass on either from \code{e.divisive()} or
#'   \code{e.agglo()}.
#'
#' @return A tibble includes which point(s) is/are the changepoint along with
#'   raw changepoint value corresponding to that changepoint. Changepoint
#'   locations follow the \code{ecp} package convention: the first index of the
#'   right segment. When no changepoint is found, an empty tibble is returned
#'   (0 rows). The upstream fit is not retained (\code{$fit} is
#'   \code{NULL} on a \code{ggcpt} from \code{cpt_detect(method = "ecp")}),
#'   because \code{ecp::e.agglo()}'s cluster-progression matrix is quadratic
#'   in the series length; call the \pkg{ecp} functions directly if you need
#'   their full output.
# `@import tibble`, `ecp` and `Rdpack` removed: no unqualified call reaches
# any of them. `RdMacros: Rdpack` plus `Imports: Rdpack` in DESCRIPTION is
# what makes `\insertRef` resolve. The NAMESPACE import was not folklore,
# though: without one, R CMD check reports "Namespace in Imports field not
# imported from: 'Rdpack'". Rdpack's own remedy, `importFrom(Rdpack,
# reprompt)`, now sits in the package block in R/ggchangepoint.R.
#' @references
#' \insertRef{james2014ecp}{ggchangepoint}
#' @export
#'
#' @examples
#' set.seed(2022)
#' ecp_wrapper(c(rnorm(100,0,1),rnorm(100,0,10)))
#' ecp_wrapper(c(rnorm(100,0,1),rnorm(100,10,1)))
#'
#' @family changepoint engines
ecp_wrapper <- function(data,
                        algorithm = "divisive",
                        min_size = 2,
                        seed = NULL,
                        ...){
  reject_renamed_args(list(...), "ecp")
  # `min_size` reaches ecp as `min.size`, so without a check here a bad value
  # was reported as "min.size must be an integer greater than 1" -- naming
  # the engine's argument rather than the one the caller passed -- and
  # `min_size = NA` reached an `if` and answered "missing value where
  # TRUE/FALSE needed". Only meaningful for algorithm = "divisive", which is
  # where ecp itself applies it, but the check is cheap and the message is
  # about the argument either way.
  validate_scalar(min_size, "min_size", min = 2)

  algorithm <- cpt_match_arg(algorithm, c("divisive", "agglo"))

  if (!is.numeric(data) && !is.data.frame(data) && !is.matrix(data)) {
    cpt_abort("`data` must be a numeric vector, matrix, or data.frame.",
              class = "bad_type")
  }
  # ecp absorbs non-finite values rather than refusing them, and what it
  # returns is wrong rather than merely missing: on a 180-point series with
  # one changepoint at 90, twenty NAs lose the changepoint entirely, and an
  # all-NA second half reports two changepoints at 12 and 14 that the data
  # does not contain. `cpt_wrapper()` and `cpt_detect()` have both always
  # refused this input; only the ecp route was open. as_mv_matrix() handles
  # the rectangular case so a non-numeric *column* is named, as everywhere
  # else in the package.
  finite_check <- if (is.matrix(data) || is.data.frame(data)) {
    as_mv_matrix(data, arg = "data")
  } else {
    data
  }
  if (anyNA(finite_check) || any(!is.finite(finite_check))) {
    stop_nonfinite(finite_check, "data")
  }

  local_seed(seed)

  estimates <- if (algorithm == "divisive") {
    ecp::e.divisive(as.matrix(data), min.size = min_size, ...)$estimates
  } else {
    ecp::e.agglo(as.matrix(data), ...)$estimates
  }

  # Keep genuine changepoints (2..n in ecp's right convention) rather than
  # stripping the first/last entries positionally: with no changepoints,
  # estimates is c(1, n+1) and `estimates[2:(length - 1)]` would evaluate to
  # the reversed boundaries; and e.agglo's wrap-around case omits the
  # boundary markers entirely, so positional stripping would drop genuine
  # changepoints.
  n_obs <- if (is.matrix(data) || is.data.frame(data)) nrow(data) else length(data)
  cp <- sort(unique(estimates[estimates > 1 & estimates <= n_obs]))

  if (length(cp) == 0) {
    return(tibble::tibble(cp = integer(), cp_value = numeric()))
  }

  # Check matrix/data.frame BEFORE the numeric test: is.numeric() is TRUE for a
  # numeric matrix, so testing it first would column-major flatten the matrix
  # and return the wrong cp_value.
  if (is.matrix(data) || is.data.frame(data)) {
    # For multivariate input, use the first column's value at the changepoint
    data_mat <- as.matrix(data)
    cp_value <- data_mat[cp, 1]
  } else {
    cp_value <- data[cp]
  }

  tibble::tibble(cp = cp, cp_value = cp_value)
}



#' Plot for the ecp package
#'
#' The plot for changepoints detected by the ecp package is a line plot for the
#' raw data and the vertical lines representing each changepoint. The x-axis is
#' the row number of the raw data in the original data vector. The plot inherits
#' ggplot2, meaning users can add ggplot2 functions on top the changepoint plot
#' for customization.
#'
#' @inheritParams ecp_wrapper
#' @inheritParams ggcptplot
#'
#' @return A line plot with data points along with the vertical lines
#'   representing changepoints.
#' @export
#'
#' @examples
#' ggecpplot(c(rnorm(100,0,1),rnorm(100,0,10)))
#' ggecpplot(c(rnorm(100,0,1),rnorm(100,10,1)))
#'
#' @family plotting
ggecpplot <- function(data,
                      algorithm = "divisive",
                      min_size = 2,
                      ...,
                      cptline_alpha = 1,
                      cptline_color = "blue",
                      cptline_type = "solid",
                      cptline_linewidth = 0.5,
                      cptline_size = lifecycle::deprecated(),
                      index = NULL,
                      show_points = NULL,
                      show_line = TRUE){

  if (lifecycle::is_present(cptline_size)) {
    lifecycle::deprecate_soft("0.2.0", "ggecpplot(cptline_size)", "ggecpplot(cptline_linewidth)")
    cptline_linewidth <- cptline_size
  }

  n_obs <- if (is.matrix(data) || is.data.frame(data)) nrow(data) else length(data)
  if (is.null(show_points)) {
    show_points <- n_obs <= 500
  }

  result <- ecp_wrapper(data, algorithm, min_size, ...)

  # The line plot is univariate; for multivariate input draw the first
  # column (and say so) rather than crashing or silently mis-plotting.
  plot_vec <- if (is.matrix(data) || is.data.frame(data)) {
    if (ncol(as.matrix(data)) > 1) {
      cpt_inform("Multivariate input: plotting the first column. ",
                 "Use autoplot(cpt_detect(data, method = \"ecp\")) for a ",
                 "faceted multivariate plot.")
    }
    as.numeric(as.matrix(data)[, 1])
  } else {
    as.numeric(data)
  }

  ggcptplot_internal(plot_vec, result, cptline_alpha, cptline_color,
                     cptline_type, cptline_linewidth, index,
                     show_points, show_line)
}