R/consistency_error.R

In OceanTracer: Marine Source Unmixing Toolkit

# Internal helper — not exported
# Computes the normalized CTS error for each tracer given pre-processed
# source and mixture matrices and a solution weight vector.
.cts_core <- function(source, mixture, solution) {
  
  # Identify numeric tracer columns (ignoring non-numeric identifier columns)
  num_cols <- vapply(source, is.numeric, logical(1L))
  src_num  <- data.matrix(source[, num_cols, drop = FALSE])
  mix_num  <- data.matrix(mixture[, num_cols, drop = FALSE])
  
  # Calculate tracer count assuming standard mean/SD pairing layout
  n_cols   <- ncol(src_num)
  n_tracer <- as.integer(n_cols / 2L)
  
  src_means <- src_num[, 1L:n_tracer, drop = FALSE]
  src_sds   <- src_num[, (n_tracer + 1L):n_cols, drop = FALSE]
  mix_means <- mix_num[, 1L:n_tracer, drop = FALSE]
  
  tracer_names <- colnames(src_means)
  
  # Calculate bounds for range normalization across source distributions
  mx  <- apply(src_means + src_sds, 2L, max)
  mn  <- apply(src_means - src_sds, 2L, min)
  rng <- mx - mn
  
  # Prevent division by zero for invariant/depleted tracers
  rng[rng == 0] <- 1.0
  
  # Vectorized min-max normalization
  src_means_norm <- sweep(sweep(src_means, 2L, mn, "-"), 2L, rng, "/")
  mix_means_norm <- sweep(sweep(mix_means, 2L, mn, "-"), 2L, rng, "/")
  
  # Predicted mixture: dot product of normalized source matrix and solution weights
  predicted <- as.vector(crossprod(src_means_norm, solution))
  
  # Compute absolute consistency error per tracer
  err <- abs(predicted - mix_means_norm[1L, ])
  
  data.frame(
    tracer            = tracer_names,
    consistency_error = as.numeric(err),
    stringsAsFactors  = FALSE
  )
}


#' Consistency error of a tracer selection for an apportionment solution
#'
#' Assesses the mathematical consistency of a tracer selection by computing
#' the normalized error between the predicted and observed tracer concentrations
#' in the mixture. A low error (typically below 0.05) for all tracers indicates
#' a consistent selection. High errors point to non-conservative tracers or
#' multiple partial solutions in the dataset.
#'
#' The error for each tracer is the absolute difference between the observed
#' mixture concentration and the virtual mixture reconstructed from the
#' apportionment solution, after normalizing by the range of the source
#' confidence intervals.
#'
#' @param data A data frame containing the characteristics of sediment sources
#'   and one mixture, in the standard fingerprinting format.
#' @param solution Either a numeric vector of apportionment weights (one per
#'   source, in the same order as they appear in \code{data}), or a single-row
#'   data frame with columns named \code{w1}, \code{w2}, ... as returned by
#'   the unmixing functions.
#'
#' @return A data frame with columns \code{tracer} and \code{consistency_error},
#'   one row per tracer.
#'
#' @references
#' Latorre, B., Lizaga, I., Gaspar, L., & Navas, A. (2021). A novel method for
#' analysing consistency and unravelling multiple solutions in sediment
#' fingerprinting. \emph{Science of The Total Environment}, \emph{789}, 147804.
#'
#' @export
consistency_error <- function(data, solution) {
  
  source   <- inputSource(data)
  mixture  <- inputSample(data)
  source_n <- nrow(source)
  
  # Strip Mean_/mean_ prefixes to standardize column headers
  colnames(source)  <- gsub("^[Mm]ean_", "", colnames(source))
  colnames(mixture) <- gsub("^[Mm]ean_", "", colnames(mixture))
  
  # Accept either a named data frame or a numeric vector
  if (is.data.frame(solution)) {
    w_cols <- grep("^w[0-9]+$", colnames(solution), value = TRUE)
    if (length(w_cols) == 0L) {
      stop("'solution' data frame must contain columns named w1, w2, ...")
    }
    sol <- as.numeric(solution[1L, w_cols])
  } else if (is.numeric(solution)) {
    sol <- solution
  } else {
    stop("'solution' must be a numeric vector or a data frame with w1, w2, ... columns.")
  }
  
  if (length(sol) != source_n) {
    stop("Length of solution (", length(sol), ") does not match number of sources (", source_n, ").")
  }
  
  if (any(!is.finite(sol))) {
    stop("Solution vector contains NA, NaN, or non-finite elements.")
  }
  
  .cts_core(source, mixture, sol)
}


# --- Backwards Compatibility Aliases ---
#' @rdname consistency_error
#' @export
CTS_error <- consistency_error

#' @rdname consistency_error
#' @export
cts_error <- consistency_error