R/Ternary_Plot.R

In OceanTracer: Marine Source Unmixing Toolkit

#' Ternary diagrams with the CI_Method function results
#'
#' Visualizes Conservativeness Index (CI) predictions using ternary diagrams. Incorporates 
#' advanced hierarchical clustering to automatically group and color-code structurally 
#' similar apportionment solutions.
#'
#' @param data List. The output generated by the CI_Method function.
#' @param tracers Numeric vector. Indices of the tracers to display. Default 1:4.
#' @param n_row Integer. Number of rows in the output grid layout. Default 1.
#' @param n_col Integer. Number of columns in the output grid layout. Default 2.
#' @param Means Logical. If TRUE, indicates the input utilizes mean and SD formats rather than raw observations. Default FALSE.
#' @param cluster Logical. If TRUE, applies hierarchical clustering to group structurally similar ternary solutions. Default FALSE.
#' @param dendro_lvls Integer. Number of branch levels to cut the clustering dendrogram. Default 5.
#' @param vrtl_sol Logical. If TRUE, overlays a theoretical or target virtual solution on the plot. Default FALSE.
#' @param vrtl Numeric vector. Apportionment values for the theoretical solution. Default c(0.05, 0.9, 0.05).
#' @param seed Integer. Random seed for reproducible cluster ordering when cluster = TRUE. Default NULL.
#' 
#' @return Ternary diagrams from the individual tracer predictions.
#' @export
Ternary_Plot <- function(data, tracers = 1:4, n_row = 1, n_col = 2, Means = FALSE, 
                         cluster = FALSE, dendro_lvls = 5, vrtl_sol = FALSE, 
                         vrtl = c(0.05, 0.9, 0.05), seed = NULL) {
  
  # Helper to filter valid ternary coordinates
  .filt3 <- function(df) {
    if (is.null(df) || ncol(df) < 3 || nrow(df) == 0) return(data.frame())
    df <- as.data.frame(df)
    valid <- complete.cases(df) & is.finite(rowSums(df)) &
      df[,1] > -0.1 & df[,1] < 1.1 &
      df[,2] > -0.1 & df[,2] < 1.1 &
      df[,3] > -0.1 & df[,3] < 1.1
    return(df[valid, , drop = FALSE])
  }
  
  meta_idx <- length(data) - 1
  sources  <- nrow(inputSource(data[[meta_idx]]))
  
  if (!sources %in% c(3, 4)) {
    stop("Ternary_Plot currently only supports datasets with 3 or 4 sources.")
  }
  
  # Extract tracer names from metadata table
  meta_df <- as.data.frame(data[[meta_idx]])
  colnames(meta_df) <- sub("^DB\\.", "", colnames(meta_df))
  
  if (Means) {
    cols_to_keep <- !(seq_len(ncol(meta_df)) %in% c(1, 2, ncol(meta_df)) | grepl("^D", colnames(meta_df)))
    t_names <- gsub("^[Mm]ean_", "", colnames(meta_df)[cols_to_keep])
  } else {
    t_names <- gsub("^[Mm]ean_", "", colnames(meta_df)[3:ncol(meta_df)])
  }
  
  # Setup clustering if requested
  color_palette <- rep("blue", length(t_names))
  clusters_levels <- NULL
  
  if (cluster) {
    if (!is.null(seed)) set.seed(seed)
    
    if (Means) {
      synth_df <- head(meta_df, -1)
      # Extract sample size if present, fallback to 30
      sample_n <- if ("n" %in% colnames(synth_df)) synth_df$n[1] else 30
      
      synthetic_data <- synth_df %>%
        rowwise() %>%
        mutate(across(starts_with("SD_"), 
                      ~ list(rnorm(sample_n, mean = get(paste0("Mean_", sub("^SD_", "", cur_column()))), sd = .)),
                      .names = "{.col}_synthetic")) %>%
        ungroup() %>%
        select(ends_with("_synthetic")) %>%
        tidyr::unnest(cols = everything()) %>%
        rename_with(~ gsub("^SD_|_synthetic$", "", .x))
      
      cor_mat <- cor(synthetic_data, use = "pairwise.complete.obs")
    } else {
      num_cols <- meta_df %>% select(-(1:2)) %>% select_if(is.numeric)
      cor_mat  <- cor(num_cols, use = "pairwise.complete.obs")
    }
    
    hc <- hclust(as.dist(1 - cor_mat), method = "ward.D2")
    clusters_levels <- cutree(hc, k = dendro_lvls)
    palette_colors  <- colorRampPalette(c("darkgreen", "orange", "purple"))(dendro_lvls)
    
    # Reorder tracer data based on hierarchical clustering
    data[1:(length(data)-2)] <- data[hc$order]
    t_names <- t_names[hc$order]
  }
  
  # Graphics Setup 
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  
  if (sources == 3) {
    par(mar = c(0, 0.1, 3, 0.1), mfcol = c(n_row, n_col))
  } else {
    par(mar = c(0, 0.1, 3, 0.1), mfcol = c(6, n_col))
  }
  
  for (i2 in tracers) {
    result <- data[[i2]]
    if (is.null(result) || ncol(result) < sources) next
    
    # Assign point color
    if (cluster && !is.null(clusters_levels) && t_names[i2] %in% names(clusters_levels)) {
      pt_col <- scales::alpha(palette_colors[clusters_levels[t_names[i2]]], 0.5)
    } else {
      pt_col <- scales::alpha("blue", 0.5)
    }
    
    if (sources == 3) {
      test <- .filt3(result[, 3:5])
      test_vrt <- data.frame(x = vrtl[1], y = vrtl[2], z = vrtl[3])
      
      TernaryPlot(xlim = c(-0.6, 0.6), tip.cex = 1, grid.lines = 4, grid.minor.lines = 1)
      if (nrow(test) > 0) TernaryPoints(test, col = pt_col, cex = 0.1)
      if (vrtl_sol) TernaryPoints(test_vrt, col = scales::alpha("red", 0.9), cex = 0.8, pch = 23, bg = "black", lwd = 2)
      title(t_names[i2], cex.main = 3.5)
      
    } else if (sources == 4) {
      w <- result[, 3:6]
      
      # 6 binary combinations for 4-source projection
      combos <- list(
        list(data = cbind(w[,1], w[,2], w[,3] + w[,4]), vrt = c(vrtl[1], vrtl[2], vrtl[3] + vrtl[4])),
        list(data = cbind(w[,2], w[,3], w[,1] + w[,4]), vrt = c(vrtl[2], vrtl[3], vrtl[1] + vrtl[4])),
        list(data = cbind(w[,3], w[,4], w[,1] + w[,2]), vrt = c(vrtl[3], vrtl[4], vrtl[1] + vrtl[2])),
        list(data = cbind(w[,4], w[,1], w[,2] + w[,3]), vrt = c(vrtl[4], vrtl[1], vrtl[2] + vrtl[3])),
        list(data = cbind(w[,1], w[,3], w[,2] + w[,4]), vrt = c(vrtl[1], vrtl[3], vrtl[2] + vrtl[4])),
        list(data = cbind(w[,2], w[,4], w[,1] + w[,3]), vrt = c(vrtl[2], vrtl[4], vrtl[1] + vrtl[3]))
      )
      
      for (idx in seq_along(combos)) {
        test     <- .filt3(combos[[idx]]$data)
        test_vrt <- as.data.frame(matrix(combos[[idx]]$vrt, nrow = 1))
        
        TernaryPlot(xlim = c(-0.6, 0.6), tip.cex = 1, grid.lines = 4, grid.minor.lines = 1)
        if (nrow(test) > 0) TernaryPoints(test, col = pt_col, cex = 0.1)
        if (vrtl_sol) TernaryPoints(test_vrt, col = scales::alpha("red", 0.9), cex = 0.8, pch = 23, bg = "black", lwd = 2)
        if (idx == 1) title(t_names[i2], cex.main = 3.5)
      }
    }
  }
}

#' @rdname Ternary_Plot
#' @export
Ternary_diagram <- Ternary_Plot

#' @rdname Ternary_Plot
#' @export
ternary_diagram <- Ternary_Plot