R/CI_Method.R

In OceanTracer: Marine Source Unmixing Toolkit

#' Conservativeness Index (CI) for individual tracer evaluation
#'
#' Quantifies how consistently each tracer constrains the mixture apportionment
#' when paired with randomly chosen complementary tracers.
#'
#' @param data Data frame in standard OceanTracer format (sources + mixture).
#' @param points Integer. Number of Monte Carlo draws per tracer. Default 2000.
#' @param Means Logical or NULL (default). Auto-detected if NULL.
#' @param triangles Character. Strategy to select accompanying tracers ("virtual" or "dataset").
#' @param seed Integer. Random seed for reproducibility. Default 123456L.
#'
#' @return A list containing tracer output data frames, the original dataset, 
#'   and the final summary CI index table.
#' @export
CI_Method <- function(data, points = 2000, Means = NULL, triangles = "virtual", seed = 123456L) {
  set.seed(seed)
  
  # ── Helper: Detect averaged format (Mean_/SD_ or mean_/sd_) ────────────────
  is_averaged_internal <- function(d) {
    if (sum(grepl("^n$", tolower(colnames(d)))) != 1L) return(FALSE)
    if ((ncol(d) - 3L) %% 2L != 0L) return(FALSE)
    ntracer <- (ncol(d) - 3L) / 2L
    nms     <- colnames(d)[-c(1L:2L)][1L:(2L * ntracer)]
    v1      <- nms[grepl("^[Mm]ean_", nms)]
    v2      <- nms[grepl("^[Ss][Dd]_", nms)]
    
    length(v1) == ntracer && length(v2) == ntracer &&
      all(gsub("^[Mm]ean_", "", v1) == gsub("^[Ss][Dd]_", "", v2))
  }
  
  if (is.null(Means)) {
    Means <- is_averaged_internal(data)
  }
  
  source_n <- nrow(inputSource(data))
  if (!source_n %in% c(2L, 3L, 4L, 5L)) {
    stop("ERROR: source_n = ", source_n, " is not implemented. Supported dimensions: 2, 3, 4, 5.")
  }
  
  # ── Prepare sources and mixtures input ────────────────────────────────────
  if (isTRUE(Means)) {
    sources1   <- data[-nrow(data), 3L:ncol(data), drop = FALSE]
    id_sources <- seq_len(nrow(data) - 1L)
    id         <- paste0("S", id_sources)
    sources    <- cbind(id = id, sources1)
    mixtures   <- inputSample(data[, 1L:((ncol(data) / 2L) + 1L), drop = FALSE])
    n_tracers  <- as.integer(((ncol(data) + 1L) / 2L) - 2L)
  } else {
    sources    <- inputSource(data)
    mixtures   <- inputSample(data)
    n_tracers  <- as.integer(ncol(data) - 2L)
  }
  
  # ── Run solver algorithm ───────────────────────────────────────────────────
  total_ops <- n_tracers * nrow(mixtures) * points
  pb <- utils::txtProgressBar(min = 0, max = total_ops, style = 3, width = 50, char = "=")
  
  results <- vector("list", n_tracers)
  for (i0 in seq_len(n_tracers) - 1L) {
    pb_offset <- i0 * nrow(mixtures) * points
    
    if (triangles == "virtual") {
      results[[i0 + 1L]] <- triangles_virtual(
        sources, mixtures, tracer = i0, iter = points, seed = seed,
        pb = pb, pb_offset = pb_offset,
        tracer_current = i0 + 1L, tracer_total = n_tracers
      )
    } else if (triangles == "dataset") {
      results[[i0 + 1L]] <- triangles_dataset(
        sources, mixtures, tracer = i0, iter = points, seed = seed,
        pb = pb, pb_offset = pb_offset,
        tracer_current = i0 + 1L, tracer_total = n_tracers
      )
    } else {
      warning("Strategy '", triangles, "' is not recognized.")
    }
  }
  close(pb)
  
  # Standardize weight column numeric formatting
  w_cols <- 3L:(2L + source_n)
  for (k in seq_along(results)) {
    results[[k]][, w_cols] <- lapply(results[[k]][, w_cols, drop = FALSE], function(x) as.numeric(as.character(x)))
  }
  
  # ── Compute Conservativeness Index (CI) ───────────────────────────────────
  # Filter solutions within acceptable range [-0.1, 1.1]
  prct_in <- numeric(n_tracers)
  CI_list <- vector("list", n_tracers)
  
  for (i in seq_along(results)) {
    res_df <- results[[i]]
    weights <- res_df[, w_cols, drop = FALSE]
    
    # Mask valid bounds per weight column
    valid_mask <- apply(weights, 1L, function(row) all(row > -0.1 & row < 1.1))
    
    # Calculate percentage inside strictly valid simplex limits [0, 1]
    in_simplex_mask <- apply(weights, 1L, function(row) all(row > 0.0 & row < 1.0))
    prct_in[i]      <- (sum(in_simplex_mask) / points) * 100.0
    
    # Build output matrix structure for tracer output
    if (!any(valid_mask)) {
      weights[] <- -1.0
    }
    
    # Combine metadata columns (ID, GOF) and weights
    CI_list[[i]] <- cbind(res_df[, 1L:2L, drop = FALSE], weights)
  }
  
  # ── Prepare output object ──────────────────────────────────────────────────
  if (isTRUE(Means)) {
    names_CI_R <- names(data)[3L:ceiling(ncol(data) / 2L)]
    names_CI_R <- gsub("Mean_", "", names_CI_R)
  } else {
    names_CI_R <- names(data)[3L:ncol(data)]
  }
  
  CI_summary <- data.frame(
    tracer   = names_CI_R,
    CI_score = prct_in,
    stringsAsFactors = FALSE
  )
  
  # Reference benchmarks derived from orthogonal standard simplexes
  ref_max <- switch(as.character(source_n),
                    "2" = 100.0,
                    "3" = 45.0,
                    "4" = 13.0,
                    "5" = 2.5,
                    100.0
  )
  
  CI_summary$CI_score_norm <- round((CI_summary$CI_score / ref_max) * 100.0, 2L)
  
  # Assign Qualitative Categories & Visual Status Indicators
  CI_summary$Category <- cut(
    CI_summary$CI_score_norm,
    breaks = c(-Inf, 20, 40, 60, 80, Inf),
    labels = c("Highly Non-Conservative", "Poor", "Moderate", "Conservative", "Highly Conservative"),
    right  = FALSE
  )
  
  CI_summary$Status <- cut(
    CI_summary$CI_score_norm,
    breaks = c(-Inf, 20, 40, 60, 80, Inf),
    labels = c("\U0001F534", "\U0001F7E0", "\U0001F7E1", "\U0001F7E2", "\U0001F31F"),
    right  = FALSE
  )
  
  # Attach DB and summary table into final result list
  CI_out <- CI_list
  CI_out$DB <- data
  CI_out$CI <- CI_summary
  
  return(CI_out)
}