R/triangles_dataset.R

In OceanTracer: Marine Source Unmixing Toolkit

#' Conservative Index evaluation using randomly sampled companion tracers
#'
#' For each tracer of interest, solves a series of minimal determined linear
#' systems to assess how consistently that tracer constrains the source
#' apportionment. In each iteration the tracer under evaluation is kept fixed
#' while the remaining tracers needed to complete the system are drawn at
#' random from the other measured tracers in the dataset.
#'
#' This function is called internally by \code{CI_Method()} and is not usually
#' invoked directly by the user.
#'
#' @param sources Data frame of source statistics: ID, one mean column per
#'   tracer, one SD column per tracer, and sample size \code{n}.
#' @param mixtures Data frame of mixture samples with the same tracer columns.
#' @param tracer Integer (0-based). Index of the tracer under evaluation.
#' @param iter Integer. Number of Monte Carlo iterations per mixture sample. Default 100L.
#' @param seed Integer. Random seed for reproducibility. Default \code{69512L}.
#' @param pb An optional \code{txtProgressBar} object created externally. Default \code{NULL}.
#' @param pb_offset Integer. Starting offset for the external progress bar. Default \code{0L}.
#' @param tracer_current Integer. Position of this tracer in the outer loop. Default \code{NULL}.
#' @param tracer_total Integer. Total number of tracers in the outer loop. Default \code{NULL}.
#' @param verbose Logical. Whether to print tracer progress messages to the console. Default is \code{TRUE}.
#'
#' @return A data frame with \code{nrow(mixtures) * iter} rows containing mixture \code{ID},
#'   \code{GOF}, and estimated source proportions (\code{w.*}).
#'
#' @seealso \code{\link{CI_Method}}
#' @importFrom stats rt
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
triangles_dataset <- function(sources, mixtures, tracer = 0L, iter = 100L, seed = 69512L,
                              pb = NULL, pb_offset = 0L,
                              tracer_current = NULL, tracer_total = NULL, verbose = TRUE) {
  
  set.seed(seed)
  
  # ── 1. Input Validation ───────────────────────────────────────────────────
  col_names <- colnames(sources)
  vars      <- ncol(mixtures) - 1L    # Number of tracers
  nsource   <- nrow(sources)
  npoints   <- nrow(mixtures)
  
  if (ncol(sources) != (2L + vars * 2L)) {
    stop("Wrong number of columns in source data frame.")
  }
  if (!tolower(col_names[1L]) %in% c("id")) {
    stop("First column of source data frame must be 'ID'.")
  }
  if (tracer < 0L || tracer >= vars) {
    stop("Requested tracer index is out of range.")
  }
  if (!nsource %in% c(3L, 4L)) {
    stop("triangles_dataset currently supports 3 or 4 sediment sources.")
  }
  
  # ── 2. Read Statistics ───────────────────────────────────────────────────
  src_mean <- as.matrix(sources[, 2L:(vars + 1L), drop = FALSE])
  src_sd   <- as.matrix(sources[, (vars + 2L):(vars * 2L + 1L), drop = FALSE])
  src_n    <- as.integer(sources[, ncol(sources)])
  
  mix_vals   <- as.matrix(mixtures[, 2L:(vars + 1L), drop = FALSE])
  source_ids <- as.character(sources[, 1L])
  
  # ── 3. Output Storage Setup ───────────────────────────────────────────────
  nrows_out <- npoints * iter
  out_id    <- integer(nrows_out)
  out_gof   <- numeric(nrows_out)
  out_w     <- matrix(0.0, nrow = nrows_out, ncol = nsource)
  
  # Helper: Perturb mean using Student's t distribution
  correct_sem <- function(avg, dev, n) {
    if (n > 1L) {
      return(avg + stats::rt(1L, df = n - 1L) * dev / sqrt(n))
    }
    return(avg)
  }
  
  all_tracers <- seq_len(vars) - 1L
  
  # ── 4. Progress Tracking ──────────────────────────────────────────────────
  local_pb <- is.null(pb)
  if (local_pb) {
    pb <- utils::txtProgressBar(min = 0, max = nrows_out, style = 3, width = 50, char = "=")
  }
  
  if (!is.null(tracer_current) && !is.null(tracer_total) && verbose) {
    message(sprintf("Processing Tracer %d of %d", tracer_current, tracer_total))
  }
  
  # Initialize missing row counter
  row_idx <- 1L
  
  # ── 5. Monte Carlo Simulation Loop ────────────────────────────────────────
  for (i in seq_len(npoints)) {
    
    mix_i <- mix_vals[i, ]
    t1    <- tracer + 1L    # Convert 0-based index to 1-based R matrix index
    
    for (ii in seq_len(iter)) {
      
      # Perturb source tracer stats
      sc <- matrix(0.0, nrow = nsource, ncol = vars)
      for (s in seq_len(nsource)) {
        for (v in seq_len(vars)) {
          sc[s, v] <- correct_sem(src_mean[s, v], src_sd[s, v], src_n[s])
        }
      }
      
      # Update progress bar efficiently
      current_step <- pb_offset + (i - 1L) * iter + ii
      if (current_step %% 20L == 0L || current_step == nrows_out) {
        if (local_pb || !is.null(pb)) {
          utils::setTxtProgressBar(pb, current_step)
        }
      }
      
      # ── 3 Sources Determinant Solver ──────────────────────────────────────
      if (nsource == 3L) {
        s11 <- sc[1L, t1];  s21 <- sc[2L, t1];  s31 <- sc[3L, t1]
        m1  <- mix_i[t1]
        
        other_tracers <- all_tracers[all_tracers != tracer]
        k             <- sample(other_tracers, 1L) + 1L
        s12 <- sc[1L, k];   s22 <- sc[2L, k];   s32 <- sc[3L, k]
        m2  <- mix_i[k]
        
        det <- -s22*s31 + s12*s31 + s21*s32 - s11*s32 + s11*s22 - s12*s21
        if (abs(det) > 1e-12) {
          w1 <- (s21*s32 - s22*s31 + m1*(s22 - s32) + m2*(s31 - s21)) / det
          w2 <- (s12*s31 - s11*s32 + m1*(s32 - s12) + m2*(s11 - s31)) / det
          w3 <- (s11*s22 - s12*s21 + m1*(s12 - s22) + m2*(s21 - s11)) / det
          w  <- c(w1, w2, w3)
        } else {
          w <- rep(0.0, 3L)
        }
      }
      
      # ── 4 Sources Determinant Solver ──────────────────────────────────────
      else {
        s11 <- sc[1L, t1];  s21 <- sc[2L, t1]
        s31 <- sc[3L, t1];  s41 <- sc[4L, t1];  m1 <- mix_i[t1]
        
        other_tracers <- all_tracers[all_tracers != tracer]
        k             <- sample(other_tracers, 1L) + 1L
        s12 <- sc[1L, k];   s22 <- sc[2L, k];   s32 <- sc[3L, k];  s42 <- sc[4L, k]
        m2  <- mix_i[k]
        
        remaining <- all_tracers[all_tracers != tracer & all_tracers != (k - 1L)]
        l         <- sample(remaining, 1L) + 1L
        s13 <- sc[1L, l];   s23 <- sc[2L, l];   s33 <- sc[3L, l];  s43 <- sc[4L, l]
        m3  <- mix_i[l]
        
        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 (abs(det) > 1e-12) {
          w1 <- (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
          
          w2 <- (-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
          
          w3 <- (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
          
          w4 <- (-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
          
          w <- c(w1, w2, w3, w4)
        } else {
          w <- rep(0.0, 4L)
        }
      }
      
      # ── Store Output Row ──────────────────────────────────────────────────
      out_id[row_idx]  <- i
      out_gof[row_idx] <- 1.0
      out_w[row_idx, ] <- ifelse(is.finite(w), w, 0.0)
      row_idx          <- row_idx + 1L
    }
  }
  
  # ── 6. Assemble Final Data Frame ──────────────────────────────────────────
  result <- data.frame(
    ID  = out_id,
    GOF = out_gof,
    out_w,
    stringsAsFactors = FALSE
  )
  
  col_names_w <- make.names(paste0("w.", source_ids))
  colnames(result)[3L:(2L + nsource)] <- col_names_w
  
  if (local_pb) close(pb)
  return(result)
}