R/cptDist.R

In changepointGA: Changepoint Detection via Modified Genetic Algorithm

#' Comparing multiple changepoint configurations by pairwise distance
#'
#' This function is to calculate the distance metric (Shi et al., 2022) between
#' two changepoint configurations, \eqn{C_{1}=\{\tau_{1}, \ldots, \tau_{m}\}}
#' and \eqn{C_{2}=\{\eta_{1}, \ldots, \eta_{k}\}}, where \eqn{\tau}'s are the
#' changepoint locations.
#'
#' The pairwise distance was proposed by Shi et al. (2022),
#' \deqn{d(C_{1}, C_{2})=|m-k|+ \min(A(C_{1}, C_{2})),}
#' where \eqn{m} is the number of changepoints in configuration \eqn{C_{1}} and
#' \eqn{k} is the number of changepoints in configuration \eqn{C_{2}}.
#' The term \eqn{\min(A(C_{1}, C_{2}))} reflects the cost of matching
#' changepoint locations between \eqn{C_{1}} and \eqn{C_{2}} and can be calculated
#' using the linear assignment method. Details can be found in Shi et al. (2022).
#' Note: if one configuration doesn't contain any changepoints (valued \code{NULL}),
#' the distance is defined as \eqn{|m-k|}. The function can also be used
#' to examine changepoint detection performance in simulation studies. Given the
#' true changepoint configuration, \eqn{C_{true}}, used in generating the time
#' series, the calculated distance between the estimated multiple changepoint
#' configuration, \eqn{C_{est}=\{\eta_{1}, \ldots, \eta_{k}\}}, and \eqn{C_{true}}
#' can be used to evaluate the performance of the detection algorithm.
#' @param tau1 A vector contains the changepoint locations for \eqn{C_{1}}
#' for comparison. A value \code{NULL} is required if there is no changepoint detected.
#' @param tau2 A vector contains the changepoint locations for  \eqn{C_{2}}
#' for comparison. A value \code{NULL} is required if there is no changepoint detected.
#' @param N The simulated time series sample size. Two changepoint configurations
#' should have the same \eqn{n} values.
#' @return
#' \item{dist}{The calculated distance.}
#' @references{
#'   Shi, X., Gallagher, C., Lund, R., & Killick, R. (2022).
#'   A comparison of single and multiple changepoint techniques for time series data.
#'   \emph{Computational Statistics & Data Analysis}, 170, 107433.
#' }
#' @useDynLib changepointGA
#' @export
#' @examples
#' N <- 100
#'
#' # both tau1 and tau2 has detected changepoints
#' tau2 <- c(25, 50, 75)
#' tau1 <- c(20, 35, 70, 80, 90)
#' cptDist(tau1 = tau1, tau2 = tau2, N = N)
#'
#' # either tau1 or tau2 has zero detected changepoints
#' cptDist(tau1 = tau1, tau2 = NULL, N = N)
#' cptDist(tau1 = NULL, tau2 = tau2, N = N)
#'
#' # both tau1 and tau2 has zero detected changepoints
#' cptDist(tau1 = NULL, tau2 = NULL, N = N)
cptDist <- function(tau1, tau2, N) {
  m <- length(tau1)
  k <- length(tau2)

  if (is.null(tau1) & is.null(tau2)) {
    dist.res <- NA
    warning("Both configurations are NULL.")
  } else {
    if (is.null(tau1) | is.null(tau2)) {
      ACC <- 0
    } else {
      if (any(tau1 < 0 | tau1 > N)) {
        stop("First changepoint configuration invalid.")
      }
      if (any(tau2 < 0 | tau2 > N)) {
        stop("Second changepoint configuration invalid.")
      }

      costs <- matrix(0, nrow = m, ncol = k)
      for (i in 1:m) {
        for (j in 1:k) {
          costs[i, j] <- abs(tau1[i] - tau2[j]) / N
        }
      }
      if (m > k) {
        costs <- t(costs)
      }
      y <- clue::solve_LSAP(costs, maximum = FALSE)
      ACC <- sum(costs[cbind(seq_along(y), y)])
    }
    dist.res <- abs(m - k) + ACC
  }

  return(dist.res)
}