R/CTS_4_sources.R

In OceanTracer: Marine Source Unmixing Toolkit

#' CTS Engine for 4 Sources
#'
#' Solves determined systems for 3 tracers and 4 sources to assess 
#' solution dispersion and discriminant capacity.
#'
#' @param source Data frame of sediment sources containing 4 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}{Triplet of tracer names used for the system.}
#'   \item{w1, w2, w3, w4}{Calculated proportion contributions for Sources 1 through 4.}
#'   \item{dw1, dw2, dw3, dw4}{Dispersion uncertainties for Sources 1 through 4 (half 84\% - 16\% quantile range).}
#'   \item{cons}{Proportion of Monte Carlo iterations yielding physically valid proportions.}
#'   \item{max_dw}{Maximum dispersion uncertainty across all four sources.}
#'
#' @importFrom stats rt quantile
#' @importFrom utils combn txtProgressBar setTxtProgressBar
#'
#' @examples
#' # Prepare small dummy datasets for 4 sources and 3 tracers
#' src <- data.frame(
#'   ID = c("Src1", "Src2", "Src3", "Src4"),
#'   Tracer1 = c(10.5, 20.2, 5.0, 15.1),
#'   Tracer2 = c(2.1, 8.5, 12.3, 4.2),
#'   Tracer3 = c(0.5, 1.2, 3.4, 0.9),
#'   Tracer1_SD = c(1.1, 1.8, 0.5, 0.8),
#'   Tracer2_SD = c(0.3, 0.9, 1.1, 0.4),
#'   Tracer3_SD = c(0.1, 0.2, 0.3, 0.1),
#'   n = c(10, 10, 10, 10)
#' )
#' mix <- data.frame(ID = "Mix1", Tracer1 = 12.0, Tracer2 = 6.0, Tracer3 = 1.5)
#'
#' # Run engine
#' result <- CTS_4_sources(source = src, mixture = mix, iter = 100, verbose = FALSE)
#' head(result)
#'
#' @export
CTS_4_sources <- function(source, mixture, iter = 1000, seed = 123456, verbose = TRUE) {
  set.seed(seed)
  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)
  
  comb_matrix <- utils::combn(ntracers, 3)
  n_comb <- ncol(comb_matrix)
  
  df_results <- data.frame()
  
  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]; t3 <- comb_matrix[3, i]
    w1 <- numeric(iter); w2 <- numeric(iter); w3 <- numeric(iter); w4 <- numeric(iter)
    sc <- 0
    
    for (l in 1:iter) {
      deg_f <- source_mat[, n_col_idx]
      if (l == 1) {
        s11 <- source_mat[1, t1]; s12 <- source_mat[1, t2]; s13 <- source_mat[1, t3]
        s21 <- source_mat[2, t1]; s22 <- source_mat[2, t2]; s23 <- source_mat[2, t3]
        s31 <- source_mat[3, t1]; s32 <- source_mat[3, t2]; s33 <- source_mat[3, t3]
        s41 <- source_mat[4, t1]; s42 <- source_mat[4, t2]; s43 <- source_mat[4, t3]
      } else {
        # Monte Carlo sampling using t-distribution
        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])
        s13 <- source_mat[1, t3] + source_mat[1, t3 + 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])
        s23 <- source_mat[2, t3] + source_mat[2, t3 + 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])
        s33 <- source_mat[3, t3] + source_mat[3, t3 + ntracers] * stats::rt(1, deg_f[3]) / sqrt(deg_f[3])
        s41 <- source_mat[4, t1] + source_mat[4, t1 + ntracers] * stats::rt(1, deg_f[4]) / sqrt(deg_f[4])
        s42 <- source_mat[4, t2] + source_mat[4, t2 + ntracers] * stats::rt(1, deg_f[4]) / sqrt(deg_f[4])
        s43 <- source_mat[4, t3] + source_mat[4, t3 + ntracers] * stats::rt(1, deg_f[4]) / sqrt(deg_f[4])
      }
      m1 <- mixture_mat[t1]; m2 <- mixture_mat[t2]; m3 <- mixture_mat[t3]
      
      det <- s21 * (s32 * s43 - s33 * s42) - s11 * (s32 * s43 - s33 * s42) - s31 * (s22 * s43 - s23 * s42) + s11 * (s22 * s43 - s23 * s42) + s31 * (s12 * s43 - s13 * s42) - s21 * (s12 * s43 - s13 * s42) + (s22 * s33 - s23 * s32) * s41 - (s12 * s33 - s13 * s32) * s41 + (s12 * s23 - s13 * s22) * s41 - s11 * (s22 * s33 - s23 * s32) + s21 * (s12 * s33 - s13 * s32) - (s12 * s23 - s13 * s22) * s31
      
      if (det != 0) {
        x1 <- (s21 * (s32 * s43 - s33 * s42) - s31 * (s22 * s43 - s23 * s42) + (s22 * s33 - s23 * s32) * s41 + m1 * (-s32 * s43 + s22 * s43 + s33 * s42 - s23 * s42 - s22 * s33 + s23 * s32) + m2 * (s31 * s43 - s21 * s43 - s33 * s41 + s23 * s41 + s21 * s33 - s23 * s31) + m3 * (-s31 * s42 + s21 * s42 + s32 * s41 - s22 * s41 - s21 * s32 + s22 * s31)) / det
        x2 <- (-s11 * (s32 * s43 - s33 * s42) + s31 * (s12 * s43 - s13 * s42) - (s12 * s33 - s13 * s32) * s41 + m1 * (s32 * s43 - s12 * s43 - s33 * s42 + s13 * s42 + s12 * s33 - s13 * s32) + m2 * (-s31 * s43 + s11 * s43 + s33 * s41 - s13 * s41 - s11 * s33 + s13 * s31) + m3 * (s31 * s42 - s11 * s42 - s32 * s41 + s12 * s41 + s11 * s32 - s12 * s31)) / det
        x3 <- (s11 * (s22 * s43 - s23 * s42) - s21 * (s12 * s43 - s13 * s42) + (s12 * s23 - s13 * s22) * s41 + m1 * (-s22 * s43 + s12 * s43 + s23 * s42 - s13 * s42 - s12 * s23 + s13 * s22) + m2 * (s21 * s43 - s11 * s43 - s23 * s41 + s13 * s41 + s11 * s23 - s13 * s21) + m3 * (-s21 * s42 + s11 * s42 + s22 * s41 - s12 * s41 - s11 * s22 + s12 * s21)) / det
        x4 <- (-s11 * (s22 * s33 - s23 * s32) + s21 * (s12 * s33 - s13 * s32) - (s12 * s23 - s13 * s22) * s31 + m1 * (s22 * s33 - s12 * s33 - s23 * s32 + s13 * s32 + s12 * s23 - s13 * s22) + m2 * (-s21 * s33 + s11 * s33 + s23 * s31 - s13 * s31 - s11 * s23 + s13 * s21) + m3 * (s21 * s32 - s11 * s32 - s22 * s31 + s12 * s31 + s11 * s22 - s12 * s21)) / det
        
        tol <- 0.05 # 5% tolerance
        if (min(x1, x2, x3, x4) >= -tol && max(x1, x2, x3, x4) <= (1 + tol)) sc <- sc + 1
        w1[l] <- x1; w2[l] <- x2; w3[l] <- x3; w4[l] <- x4
      }
    }
    
    res <- data.frame(
      id = paste(tracer_names[t1], tracer_names[t2], tracer_names[t3], sep = " "),
      w1 = w1[1], w2 = w2[1], w3 = w3[1], w4 = w4[1],
      dw1 = (as.numeric(stats::quantile(w1, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w1, 0.16, na.rm = TRUE))) / 2,
      dw2 = (as.numeric(stats::quantile(w2, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w2, 0.16, na.rm = TRUE))) / 2,
      dw3 = (as.numeric(stats::quantile(w3, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w3, 0.16, na.rm = TRUE))) / 2,
      dw4 = (as.numeric(stats::quantile(w4, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w4, 0.16, na.rm = TRUE))) / 2,
      cons = sc / iter
    )
    df_results <- rbind(df_results, res)
    if (verbose) utils::setTxtProgressBar(pb, i)
  }
  
  if (verbose) close(pb)
  
  df_results$max_dw <- apply(df_results[, c("dw1", "dw2", "dw3", "dw4")], 1, max)
  df_results <- df_results[order(df_results$max_dw), ]
  rownames(df_results) <- NULL
  
  return(df_results)
}

# --- Aliases ---
#' @rdname CTS_4_sources
#' @export
triplets <- CTS_4_sources

#' @rdname CTS_4_sources
#' @export
CTS_seeds_triplets <- CTS_4_sources