R/CTS_3_sources.R

In OceanTracer: Marine Source Unmixing Toolkit

#' CTS Engine for 3 Sources
#'
#' Solves determined systems for 2 tracers and 3 sources to assess 
#' solution dispersion and discriminant capacity.
#'
#' @param source Data frame of sediment sources containing 3 source rows, mean values, standard deviations, and sample size ('n').
#' @param mixture Data frame of mixture samples.
#' @param iter Number of Monte Carlo iterations. Default is 1000.
#' @param seed Integer seed for reproducibility. Default is 123456.
#' @param verbose Logical. Whether to show a progress bar in the console. Default is TRUE.
#'
#' @return A \code{data.frame} ordered by \code{max_dw} (lowest dispersion first) containing:
#'   \item{id}{Pair of tracer names used for the system.}
#'   \item{w1, w2, w3}{Calculated proportion contributions for Sources 1, 2, and 3.}
#'   \item{dw1, dw2, dw3}{Dispersion uncertainties for Sources 1, 2, and 3 (half 84\% - 16\% quantile range).}
#'   \item{cons}{Proportion of Monte Carlo iterations yielding physically valid proportions between 0 and 1.}
#'   \item{max_dw}{Maximum dispersion uncertainty across all three sources.}
#'
#' @importFrom stats rt quantile
#' @importFrom utils combn txtProgressBar setTxtProgressBar
#'
#' @examples
#' # Prepare small dummy datasets for 3 sources and 2 tracers
#' src <- data.frame(
#'   ID = c("Src1", "Src2", "Src3"),
#'   Tracer1 = c(10.5, 20.2, 5.0),
#'   Tracer2 = c(2.1, 8.5, 12.3),
#'   Tracer1_SD = c(1.1, 1.8, 0.5),
#'   Tracer2_SD = c(0.3, 0.9, 1.1),
#'   n = c(10, 10, 10)
#' )
#' mix <- data.frame(ID = "Mix1", Tracer1 = 12.0, Tracer2 = 6.0)
#'
#' # Run engine
#' result <- CTS_3_sources(source = src, mixture = mix, iter = 100, verbose = FALSE)
#' head(result)
#'
#' @export
CTS_3_sources <- function(source, mixture, iter = 1000, seed = 123456, verbose = TRUE) {
  
  set.seed(seed)
  # Convert to matrix excluding ID
  source_mat <- data.matrix(source[-1])
  mixture_mat <- data.matrix(mixture[-1])
  
  ntracers <- (ncol(source_mat) - 1) / 2
  tracer_names <- colnames(source_mat)[1:ntracers]
  n_col_idx <- ncol(source_mat)
  
  df_results <- data.frame()
  
  # Combinations of 2 tracers for 3 sources
  comb_matrix <- utils::combn(ntracers, 2)
  n_comb <- ncol(comb_matrix)
  
  if (verbose) {
    pb <- utils::txtProgressBar(min = 0, max = n_comb, style = 3, width = 50, char = "=")
  }
  
  for (i in 1:n_comb) {
    t1 <- comb_matrix[1, i]
    t2 <- comb_matrix[2, i]
    
    w1 <- numeric(iter); w2 <- numeric(iter); w3 <- numeric(iter)
    sc <- 0
    
    for (l in 1:iter) {
      if (l == 1) {
        s11 <- source_mat[1, t1]; s12 <- source_mat[1, t2]
        s21 <- source_mat[2, t1]; s22 <- source_mat[2, t2]
        s31 <- source_mat[3, t1]; s32 <- source_mat[3, t2]
        m1  <- mixture_mat[t1];   m2  <- mixture_mat[t2]
      } else {
        # Variability using t-distribution and 'n' column
        deg_f <- source_mat[, n_col_idx]
        s11 <- source_mat[1, t1] + source_mat[1, t1 + ntracers] * stats::rt(1, deg_f[1]) / sqrt(deg_f[1])
        s12 <- source_mat[1, t2] + source_mat[1, t2 + ntracers] * stats::rt(1, deg_f[1]) / sqrt(deg_f[1])
        s21 <- source_mat[2, t1] + source_mat[2, t1 + ntracers] * stats::rt(1, deg_f[2]) / sqrt(deg_f[2])
        s22 <- source_mat[2, t2] + source_mat[2, t2 + ntracers] * stats::rt(1, deg_f[2]) / sqrt(deg_f[2])
        s31 <- source_mat[3, t1] + source_mat[3, t1 + ntracers] * stats::rt(1, deg_f[3]) / sqrt(deg_f[3])
        s32 <- source_mat[3, t2] + source_mat[3, t2 + ntracers] * stats::rt(1, deg_f[3]) / sqrt(deg_f[3])
        m1  <- mixture_mat[t1];   m2  <- mixture_mat[t2]
      }
      
      # Optimized Determinant (Maxima script)
      det <- -s22 * s31 + s12 * s31 + s21 * s32 - s11 * s32 + s11 * s22 - s12 * s21
      
      if (det != 0) {
        x1 <- (s21 * s32 - s22 * s31 + m1 * (s22 - s32) + m2 * (s31 - s21)) / det
        x2 <- (s12 * s31 - s11 * s32 + m1 * (s32 - s12) + m2 * (s11 - s31)) / det
        x3 <- (s11 * s22 - s12 * s21 + m1 * (s12 - s22) + m2 * (s21 - s11)) / det
        
        if (min(x1, x2, x3) >= 0 && max(x1, x2, x3) <= 1) sc <- sc + 1
        
        w1[l] <- x1; w2[l] <- x2; w3[l] <- x3
      }
    }
    
    # Uncertainty (16-84% quantiles)
    dw1 <- (stats::quantile(w1, 0.84, na.rm = TRUE) - stats::quantile(w1, 0.16, na.rm = TRUE)) / 2
    dw2 <- (stats::quantile(w2, 0.84, na.rm = TRUE) - stats::quantile(w2, 0.16, na.rm = TRUE)) / 2
    dw3 <- (stats::quantile(w3, 0.84, na.rm = TRUE) - stats::quantile(w3, 0.16, na.rm = TRUE)) / 2
    
    res <- data.frame(
      id = paste(tracer_names[t1], tracer_names[t2], sep = " "),
      w1 = w1[1], w2 = w2[1], w3 = w3[1],
      dw1 = dw1, dw2 = dw2, dw3 = dw3,
      cons = sc / iter,
      max_dw = max(dw1, dw2, dw3)
    )
    
    df_results <- rbind(df_results, res)
    if (verbose) utils::setTxtProgressBar(pb, i)
  }
  
  if (verbose) close(pb)
  
  df_results <- df_results[order(df_results$max_dw), ]
  rownames(df_results) <- NULL
  
  return(df_results)
}

# --- Aliases ---
#' @rdname CTS_3_sources
#' @export
pairs <- CTS_3_sources

#' @rdname CTS_3_sources
#' @export
CTS_seeds_pairs <- CTS_3_sources