R/triangles_virtual.R

In OceanTracer: Marine Source Unmixing Toolkit

#' Conservative Index evaluation using virtually generated 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 replaced by
#' uniformly random values in \eqn{[0, 1]} (virtual approach).
#'
#' @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 \code{100L}.
#' @param seed Integer. Random seed for reproducibility. Default \code{69512L}.
#' @param pb External \code{txtProgressBar} object. Default \code{NULL}.
#' @param pb_offset Integer. Starting offset for external progress bar. Default \code{0L}.
#' @param tracer_current Integer. Position of this tracer in outer loop. Default \code{NULL}.
#' @param tracer_total Integer. Total number of tracers in 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 ID,
#'   GOF (fixed at 1), and apportioned weights per source.
#'
#' @importFrom stats rt runif
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
triangles_virtual <- 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)
  
  # ── Input validation ───────────────────────────────────────────────────────
  col_names <- colnames(sources)
  vars      <- ncol(mixtures) - 1L      # number of tracers
  nsource   <- nrow(sources)
  npoints   <- nrow(mixtures)
  
  if (ncol(sources) != 2L + (ncol(mixtures) - 1L) * 2L) {
    stop("Wrong number of columns in source data frame")
  }
  if (!tolower(col_names[1L]) %in% c("id")) {
    stop("First data frame column must be ID")
  }
  if (tracer < 0L || tracer >= vars) {
    stop("Requested tracer is out of range")
  }
  if (!nsource %in% c(2L, 3L, 4L, 5L)) {
    stop("Number of sources not implemented")
  }
  
  # ── Matrices extraction ────────────────────────────────────────────────────
  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])
  
  # ── Output allocation ──────────────────────────────────────────────────────
  nrows_out <- npoints * iter
  out_id    <- integer(nrows_out)
  out_gof   <- rep(1.0, nrows_out)
  out_w     <- matrix(0.0, nrow = nrows_out, ncol = nsource)
  
  # Helper: SEM correction vectorization
  correct_sem_vec <- function(means, sds, ns, n_src, n_var) {
    # Degrees of freedom matrix
    df_mat <- matrix(rep(ns - 1L, n_var), nrow = n_src, ncol = n_var)
    n_mat  <- matrix(rep(sqrt(ns), n_var), nrow = n_src, ncol = n_var)
    
    t_draws <- matrix(stats::rt(n_src * n_var, df = df_mat), nrow = n_src, ncol = n_var)
    return(means + t_draws * sds / n_mat)
  }
  
  # ── Progress bar setup ─────────────────────────────────────────────────────
  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 out of %d", tracer_current, tracer_total))
  }
  
  # ── Main processing loop ───────────────────────────────────────────────────
  row_idx <- 1L
  t1 <- tracer + 1L  # Convert 0-based index to 1-based R index
  
  for (i in seq_len(npoints)) {
    mix_i <- mix_vals[i, ]
    
    for (ii in seq_len(iter)) {
      
      # Perturb source values
      sc <- correct_sem_vec(src_mean, src_sd, src_n, nsource, vars)
      
      # Prepare minimal system equations
      # Equation 1: Real tracer values
      # Equations 2..nsource: Virtual companion tracers
      A <- matrix(0.0, nrow = nsource, ncol = nsource)
      b <- numeric(nsource)
      
      # Populate Tracer 1 (Real)
      A[1L, ] <- sc[, t1]
      b[1L]   <- mix_i[t1]
      
      # Populate Virtual Tracers (if K > 2)
      if (nsource > 2L) {
        for (v_idx in 2L:nsource) {
          rnd <- stats::runif(nsource)
          mn_v <- min(rnd)
          mx_v <- max(rnd)
          
          A[v_idx, ] <- rnd
          b[v_idx]   <- mn_v + stats::runif(1L) * (mx_v - mn_v)
        }
      }
      
      # Solve linear system A %*% w = b subject to sum(w) = 1 via augmentation
      # Equality constraint: sum(w) = 1 replaced or solved directly
      w <- tryCatch({
        # For nsource == 2, simple exact scalar formula
        if (nsource == 2L) {
          det_val <- sc[1L, t1] - sc[2L, t1]
          if (abs(det_val) > 1e-12) {
            w1 <- (mix_i[t1] - sc[2L, t1]) / det_val
            c(w1, 1.0 - w1)
          } else {
            rep(0.0, 2L)
          }
        } else {
          # Solve system with qr.solve for stability
          as.vector(qr.solve(A, b))
        }
      }, error = function(e) {
        rep(0.0, nsource)
      })
      
      # Store iteration output
      out_id[row_idx]  <- i
      out_w[row_idx, ] <- w
      
      if (!local_pb) {
        utils::setTxtProgressBar(pb, pb_offset + row_idx)
      } else {
        utils::setTxtProgressBar(pb, row_idx)
      }
      
      row_idx <- row_idx + 1L
    }
  }
  
  # ── Construct final result ─────────────────────────────────────────────────
  result <- data.frame(
    ID  = out_id,
    GOF = out_gof,
    out_w,
    stringsAsFactors = FALSE
  )
  colnames(result)[3L:(2L + nsource)] <- make.names(paste0("w.", source_ids))
  
  if (local_pb) close(pb)
  return(result)
}