R/CTS_2_sources.R

In OceanTracer: Marine Source Unmixing Toolkit

#' CTS Engine for 2 Sources
#'
#' Solves determined systems for 1 tracer and 2 sources to assess 
#' solution dispersion and discriminant capacity.
#'
#' @param source Data frame of sediment sources containing source means, standard deviations, and sample sizes.
#' @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}{Tracer name identifier.}
#'   \item{w1, w2}{Calculated proportion contributions for Source 1 and Source 2.}
#'   \item{dw1, dw2}{Dispersion uncertainties for Source 1 and Source 2 (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 both sources.}
#'
#' @importFrom stats rt quantile
#' @importFrom utils txtProgressBar setTxtProgressBar
#'
#' @examples
#' # Prepare small dummy datasets
#' src <- data.frame(
#'   ID = c("Src1", "Src2"),
#'   Tracer1 = c(10.5, 20.2),
#'   Tracer1_SD = c(1.1, 1.8),
#'   n = c(10, 10)
#' )
#' mix <- data.frame(ID = "Mix1", Tracer1 = 15.0)
#'
#' # Run engine
#' result <- CTS_2_sources(source = src, mixture = mix, iter = 100, verbose = FALSE)
#' head(result)
#'
#' @export
CTS_2_sources <- function(source, mixture, iter = 1000, seed = 123456, verbose = TRUE) {
  
  set.seed(seed)
  # Convert to matrix excluding the first column (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) # Column for 'n'
  
  df_results <- data.frame()
  
  if (verbose) {
    pb <- utils::txtProgressBar(min = 0, max = ntracers, style = 3, width = 50, char = "=")
  }
  
  for (i in 1:ntracers) {
    w1 <- numeric(iter)
    w2 <- numeric(iter)
    sc <- 0 # Consistency counter
    
    for (l in 1:iter) {
      if (l == 1) {
        # Mean values for the first iteration
        s11 <- source_mat[1, i]
        s21 <- source_mat[2, i]
        m1  <- mixture_mat[i]
      } else {
        # Variability using t-distribution and the 'n' column
        s11 <- source_mat[1, i] + source_mat[1, i + ntracers] * stats::rt(1, source_mat[1, n_col_idx]) / sqrt(source_mat[1, n_col_idx])
        s21 <- source_mat[2, i] + source_mat[2, i + ntracers] * stats::rt(1, source_mat[2, n_col_idx]) / sqrt(source_mat[2, n_col_idx])
        m1  <- mixture_mat[i]
      }
      
      # System determinant (s21 - s11)
      det <- s21 - s11
      
      if (det != 0) {
        x1 <- (s21 - m1) / det
        x2 <- (m1 - s11) / det
        
        # Check physical feasibility (0-1)
        if (min(x1, x2) >= 0 && max(x1, x2) <= 1) sc <- sc + 1
        
        w1[l] <- x1
        w2[l] <- x2
      }
    }
    
    # Calculate dispersion (Dw) using 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
    
    res <- data.frame(
      id = tracer_names[i],
      w1 = w1[1],
      w2 = w2[1],
      dw1 = dw1,
      dw2 = dw2,
      cons = sc / iter,
      max_dw = max(dw1, dw2)
    )
    
    df_results <- rbind(df_results, res)
    if (verbose) utils::setTxtProgressBar(pb, i)
  }
  
  if (verbose) close(pb)
  
  # Sort by lowest dispersion (most discriminant)
  df_results <- df_results[order(df_results$max_dw), ]
  rownames(df_results) <- NULL
  
  return(df_results)
}

# --- Compatibility Aliases ---

#' @rdname CTS_2_sources
#' @export
singles <- CTS_2_sources

#' @rdname CTS_2_sources
#' @export
CTS_seeds_singles <- CTS_2_sources