R/CTS_5_sources.R

In OceanTracer: Marine Source Unmixing Toolkit

#' CTS Engine for 5 Sources
#'
#' Solves determined linear systems for 4 tracers and 5 sources using Monte Carlo 
#' simulation to assess solution dispersion and discriminant capacity.
#'
#' @param source Data frame of sediment sources containing 5 source rows, tracer mean values, standard deviations, and sample size ('n').
#' @param mixture Data frame of mixture samples containing tracer measurements.
#' @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}{Quartet of tracer names used for the system.}
#'   \item{w1, w2, w3, w4, w5}{Calculated proportion contributions for Sources 1 through 5.}
#'   \item{dw1, dw2, dw3, dw4, dw5}{Dispersion uncertainties for Sources 1 through 5 (half 84\% - 16\% quantile range).}
#'   \item{cons}{Proportion of Monte Carlo iterations yielding physically valid proportions.}
#'   \item{max_dw}{Maximum dispersion uncertainty across all five sources.}
#'
#' @importFrom stats rt quantile
#' @importFrom utils combn txtProgressBar setTxtProgressBar
#'
#' @examples
#' # Prepare small dummy datasets for 5 sources and 4 tracers
#' src <- data.frame(
#'   ID = c("Src1", "Src2", "Src3", "Src4", "Src5"),
#'   Tracer1 = c(10.5, 20.2, 5.0, 15.1, 8.4),
#'   Tracer2 = c(2.1, 8.5, 12.3, 4.2, 6.7),
#'   Tracer3 = c(0.5, 1.2, 3.4, 0.9, 2.1),
#'   Tracer4 = c(100, 250, 400, 180, 310),
#'   Tracer1_SD = c(1.1, 1.8, 0.5, 0.8, 0.6),
#'   Tracer2_SD = c(0.3, 0.9, 1.1, 0.4, 0.5),
#'   Tracer3_SD = c(0.1, 0.2, 0.3, 0.1, 0.2),
#'   Tracer4_SD = c(10, 20, 35, 15, 22),
#'   n = c(10, 10, 10, 10, 10)
#' )
#' mix <- data.frame(ID = "Mix1", Tracer1 = 12.0, Tracer2 = 6.0, Tracer3 = 1.5, Tracer4 = 220)
#'
#' # Run engine
#' result <- CTS_5_sources(source = src, mixture = mix, iter = 100, verbose = FALSE)
#' head(result)
#'
#' @export
CTS_5_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, 4)
  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]; t4 <- comb_matrix[4, i]
    
    w1_v <- numeric(iter); w2_v <- numeric(iter); w3_v <- numeric(iter); w4_v <- numeric(iter); w5_v <- numeric(iter)
    sc <- 0
    
    deg_f <- source_mat[, n_col_idx]
    
    for (l in 1:iter) {
      if (l == 1) {
        s11 <- source_mat[1, t1]; s12 <- source_mat[1, t2]; s13 <- source_mat[1, t3]; s14 <- source_mat[1, t4]
        s21 <- source_mat[2, t1]; s22 <- source_mat[2, t2]; s23 <- source_mat[2, t3]; s24 <- source_mat[2, t4]
        s31 <- source_mat[3, t1]; s32 <- source_mat[3, t2]; s33 <- source_mat[3, t3]; s34 <- source_mat[3, t4]
        s41 <- source_mat[4, t1]; s42 <- source_mat[4, t2]; s43 <- source_mat[4, t3]; s44 <- source_mat[4, t4]
        s51 <- source_mat[5, t1]; s52 <- source_mat[5, t2]; s53 <- source_mat[5, t3]; s54 <- source_mat[5, t4]
      } else {
        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])
        s14 <- source_mat[1, t4] + source_mat[1, t4 + 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])
        s24 <- source_mat[2, t4] + source_mat[2, t4 + 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])
        s34 <- source_mat[3, t4] + source_mat[3, t4 + 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])
        s44 <- source_mat[4, t4] + source_mat[4, t4 + ntracers] * stats::rt(1, deg_f[4]) / sqrt(deg_f[4])
        s51 <- source_mat[5, t1] + source_mat[5, t1 + ntracers] * stats::rt(1, deg_f[5]) / sqrt(deg_f[5])
        s52 <- source_mat[5, t2] + source_mat[5, t2 + ntracers] * stats::rt(1, deg_f[5]) / sqrt(deg_f[5])
        s53 <- source_mat[5, t3] + source_mat[5, t3 + ntracers] * stats::rt(1, deg_f[5]) / sqrt(deg_f[5])
        s54 <- source_mat[5, t4] + source_mat[5, t4 + ntracers] * stats::rt(1, deg_f[5]) / sqrt(deg_f[5])
      }
      m1 <- mixture_mat[t1]; m2 <- mixture_mat[t2]; m3 <- mixture_mat[t3]; m4 <- mixture_mat[t4]
      
      A <- matrix(c(
        1,   1,   1,   1,   1,
        s11, s21, s31, s41, s51,
        s12, s22, s32, s42, s52,
        s13, s23, s33, s43, s53,
        s14, s24, s34, s44, s54
      ), nrow = 5, ncol = 5, byrow = TRUE)
      
      B <- c(1, m1, m2, m3, m4)
      sol <- try(solve(A, B), silent = TRUE)
      
      if (!inherits(sol, "try-error")) {
        tol <- 0.05  # 5% tolerance
        if (min(sol) >= -tol && max(sol) <= (1 + tol)) sc <- sc + 1
        w1_v[l] <- sol[1]; w2_v[l] <- sol[2]; w3_v[l] <- sol[3]; w4_v[l] <- sol[4]; w5_v[l] <- sol[5]
      }
    }
    
    res <- data.frame(
      id = paste(tracer_names[t1], tracer_names[t2], tracer_names[t3], tracer_names[t4], sep = " "),
      w1 = w1_v[1], w2 = w2_v[1], w3 = w3_v[1], w4 = w4_v[1], w5 = w5_v[1],
      dw1 = (as.numeric(stats::quantile(w1_v, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w1_v, 0.16, na.rm = TRUE))) / 2,
      dw2 = (as.numeric(stats::quantile(w2_v, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w2_v, 0.16, na.rm = TRUE))) / 2,
      dw3 = (as.numeric(stats::quantile(w3_v, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w3_v, 0.16, na.rm = TRUE))) / 2,
      dw4 = (as.numeric(stats::quantile(w4_v, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w4_v, 0.16, na.rm = TRUE))) / 2,
      dw5 = (as.numeric(stats::quantile(w5_v, 0.84, na.rm = TRUE)) - as.numeric(stats::quantile(w5_v, 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", "dw5")], 1, max)
  df_results <- df_results[order(df_results$max_dw), ]
  rownames(df_results) <- NULL
  
  return(df_results)
}

# --- Aliases ---
#' @rdname CTS_5_sources
#' @export
quartets <- CTS_5_sources

#' @rdname CTS_5_sources
#' @export
CTS_seeds_quartets <- CTS_5_sources