Update Sobol GSA to use sensobol package and validate ensemble design handling

CHANGELOG.md

@@ -9,6 +9,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 ## Unreleased
 
 ### Added
+- Added `sensobol` based global sensitivity design and postprocessing support in `PEcAn.uncertainty`, including saved Sobol design and index outputs.
 - Added PEcAn.PEPRMT model, including a demo run with example data
 - Add `format_try_for_ma()` and `try_trait_mapping()` to `PEcAn.data.remote` to convert trait data from the external TRY database into the tabular format required by the PEcAn meta-analysis module (#3717).
 - Add function `qsub_sda()` for submitting SDA batch jobs by splitting a large number of sites into multiple small groups of sites (#3634).

base/workflow/R/run.write.configs.R

@@ -36,9 +36,9 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
   # number of (pft, trait, quantile) combinations
   if (!is.null(input_design) && "ensemble" %in% names(settings)) {
     if (nrow(input_design) != ensemble.size) {
-      stop(
+      PEcAn.logger::logger.error(
         "input_design has ", nrow(input_design), " rows, but ensemble.size is ",
-        ensemble.size, ".The design matrix must have exactly one row per run."
+        ensemble.size, ". The design matrix must have exactly one row per run."
       )
     }
   }
@@ -137,6 +137,26 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
         !is.null(input_design) && 
         "param" %in% colnames(input_design)) {
       trait_sample_indices <- input_design[["param"]]
+      if (is.null(trait.samples) || length(trait.samples) == 0) {
+        PEcAn.logger::logger.error(
+          "samples.Rdata does not contain trait.samples required for input_design$param"
+        )
+      }
+      trait_sample_indices <- as.integer(trait_sample_indices)
+      if (any(is.na(trait_sample_indices)) || any(trait_sample_indices < 1L)) {
+        PEcAn.logger::logger.error("input_design$param must contain positive integer indices")
+      }
+      parameter_bank_size <- PEcAn.uncertainty::trait_sample_bank_size(trait.samples)
+      if (parameter_bank_size == 0L) {
+        PEcAn.logger::logger.error(
+          "samples.Rdata does not contain usable trait.samples required for input_design$param"
+        )
+      }
+      if (any(trait_sample_indices > parameter_bank_size)) {
+        PEcAn.logger::logger.error(
+          "input_design$param includes indices beyond the available parameter sample bank"
+        )
+      }
       ensemble.samples <- list()
       for (pft in names(trait.samples)) {
         pft_traits <- trait.samples[[pft]]
@@ -148,6 +168,10 @@ run.write.configs <- function(settings, ensemble.size, input_design, write = TRU
         )
         names(ensemble.samples[[pft]]) <- names(pft_traits)
       }
+    } else if ("ensemble" %in% names(settings) && !is.null(input_design)) {
+      PEcAn.logger::logger.error(
+        "input_design must include a `param` column selecting rows from trait.samples"
+      )
     } else {
       # use pre-generated samples
       ensemble.samples <- existing_data$ensemble.samples

base/workflow/tests/testthat/test.run.write.configs.sobol.R

@@ -0,0 +1,83 @@
+test_that("runModule.run.write.configs uses input_design row count", {
+  settings <- PEcAn.settings::Settings(
+    ensemble = list(
+      size = 3,
+      samplingspace = list(parameters = list(method = "uniform"))
+    ),
+    database = list(bety = list(write = FALSE)),
+    pfts = list(list(posterior.files = "post.distns.Rdata"))
+  )
+  input_design <- data.frame(param = seq_len(5))
+  captured <- new.env(parent = emptyenv())
+
+  mockery::stub(
+    runModule.run.write.configs,
+    "PEcAn.workflow::run.write.configs",
+    function(settings, ensemble.size, input_design, write, posterior.files, overwrite) {
+      captured$ensemble.size <- ensemble.size
+      captured$input_design <- input_design
+      list(
+        ensemble = list(ensemble.id = 123),
+        pfts = settings$pfts
+      )
+    }
+  )
+
+  result <- runModule.run.write.configs(
+    settings = settings,
+    input_design = input_design
+  )
+
+  expect_equal(captured$ensemble.size, 5)
+  expect_identical(captured$input_design, input_design)
+  expect_equal(result$ensemble$ensemble.id, 123)
+})
+
+test_that("run.write.configs validates param against the shortest trait sample bank", {
+  withr::with_tempdir({
+    trait.samples <- list(
+      pftA = list(SLA = seq_len(4)),
+      pftB = list(SLA = seq_len(3))
+    )
+    sa.samples <- list()
+    runs.samples <- list()
+    env.samples <- list()
+    ensemble.samples <- list()
+    save(
+      trait.samples,
+      sa.samples,
+      runs.samples,
+      env.samples,
+      ensemble.samples,
+      file = file.path(getwd(), "samples.Rdata")
+    )
+
+    settings <- list(
+      outdir = getwd(),
+      database = list(),
+      model = list(type = "FAKE"),
+      ensemble = list(size = 2),
+      pfts = list(
+        list(name = "pftA", posteriorid = NULL),
+        list(name = "pftB", posteriorid = NULL)
+      )
+    )
+    run_write_configs <- PEcAn.workflow::run.write.configs
+    mockery::stub(
+      run_write_configs,
+      "PEcAn.logger::logger.error",
+      function(...) stop(paste(...), call. = FALSE)
+    )
+
+    expect_error(
+      run_write_configs(
+        settings = settings,
+        ensemble.size = 2,
+        input_design = data.frame(param = c(1L, 4L)),
+        write = FALSE,
+        overwrite = FALSE
+      ),
+      "input_design\\$param includes indices beyond the available parameter sample bank"
+    )
+  })
+})

modules/uncertainty/DESCRIPTION

@@ -43,7 +43,8 @@ Imports:
     purrr,
     randtoolbox,
     rlang,
-    sensitivity
+    sensobol,
+    tibble
 Suggests:
     mockery,
     testthat (>= 1.0.2),

modules/uncertainty/NAMESPACE

@@ -33,6 +33,7 @@ export(sd.var)
 export(sensitivity.analysis)
 export(sensitivity.filename)
 export(spline.truncate)
+export(trait_sample_bank_size)
 export(write.ensemble.configs)
 export(write.sa.configs)
 importFrom(dplyr,"%>%")

modules/uncertainty/NEWS.md

@@ -1,5 +1,7 @@
 # PEcAn.uncertainty 1.9.0.9000
 
+* Switched Sobol global sensitivity workflows to `sensobol`, including explicit design metadata and saved Sobol index outputs.
+* Treat all inputs (met, initial conditions, events, etc..) as independent Sobol factors; removed parent-child filter that previously confounded events with meteorology.
 * run.ensemble.analysis() now respects `settings$modeloutdir` rather than assuming an `out/` folder inside `settings$outdir` (@Akash-paluvai, #3722).
 * Added `generate_OAT_SA_design()` for creating input design matrices for
 sensitivity analysis. This function ensures non-parameter inputs

modules/uncertainty/R/compute_sobol_indices.R

@@ -1,45 +1,145 @@
 #' Compute Sobol indices from a finished PEcAn run
 #'
-#' Loads model outputs from a Sobol ensemble, calculates summary
-#' statistics for a chosen variable, feeds them to \code{sensitivity::tell()},
-#' and returns the updated Sobol object.
+#' Loads standardized ensemble output for a Sobol run, computes first-order and
+#' total-order Sobol indices with \code{sensobol}, and saves both the Sobol design
+#' metadata and the computed indices using PEcAn-style filenames.
 #'
-#' @param outdir     PEcAn run output directory that contains runs.txt
-#' @param sobol_obj  object produced by PEcAn.uncertainty::generate_joint_ensemble_design()
-#' @param var        Variable name to summarise (default "GPP").
-#' @param stat_fun   Summary statistic applied to var default mean .
+#' First-order indices quantify the share of output variance attributable to a
+#' factor alone, while total-order indices summarize the full contribution of
+#' that factor including its interactions with other factors. PEcAn uses
+#' \code{sensobol} for these variance-based estimators; see Saltelli et al. (2008)
+#' for methodological background and Puy et al. (2022) for package details.
 #'
-#' @return           sobol_obj 
-#' .
+#' This function handles one output variable at a time. To compute indices for
+#' multiple variables, call it in a loop (see examples).
+#'
+#' @param outdir PEcAn run output directory containing \code{ensemble.output.*.Rdata}
+#'   files.
+#' @param sobol_obj object produced by
+#'   \code{PEcAn.uncertainty::generate_joint_ensemble_design(..., sobol = TRUE)}.
+#' @param var Variable name to summarize (default \code{"GPP"}).
+#'
+#' @examples
+#' \dontrun{
+#'   # single variable
+#'   result <- compute_sobol_indices(outdir, sobol_obj, var = "GPP")
+#'
+#'   # multiple variables
+#'   vars <- c("GPP", "NPP", "TotSoilCarb")
+#'   all_results <- purrr::map_dfr(vars, function(v) {
+#'     compute_sobol_indices(outdir, sobol_obj, var = v) |>
+#'       dplyr::mutate(variable = v)
+#'   })
+#' }
+#'
+#' @return A tibble of Sobol first-order and total-order indices with attached
+#'   factor metadata.
+#' @references Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni,
+#'   J., Gatelli, D., et al. (2008). Global Sensitivity Analysis: The Primer.
+#'   John Wiley & Sons.
+#'
+#'   Puy, A., Lo Piano, S., Saltelli, A., and Levin, S. A. (2022).
+#'   sensobol: An R Package to Compute Variance-Based Sensitivity Indices.
+#'   Journal of Statistical Software, 102(5), 1-37.
+#'   \doi{10.18637/jss.v102.i05}
 #' @export
 compute_sobol_indices <- function(outdir,
                                   sobol_obj,
-                                  var = "GPP",
-                                  stat_fun = mean) {
-
-
-
-  runs_file <- file.path(outdir, "runs.txt")
-  if (!file.exists(runs_file)) {
-    PEcAn.logger::logger.error("runs.txt not found in ", outdir)
-   }
-  run_ids <- readLines(runs_file) 
-
-
-
-  # Load outputs and compute response vector y
-  y <- vapply(run_ids, function(rid) {
-    fpath <- file.path(outdir, rid)
-    out   <- PEcAn.utils::read.output(runid = rid, outdir = fpath)
-    if (!is.list(out) || !var %in% names(out)) {
-      PEcAn.logger::logger.error("Variable '", var, "' missing in output for run ", rid)
-    }
-    stat_fun(out[[var]], na.rm = TRUE)
-  }, numeric(1))
-
-  # Compute Sobol indices
-  sobol_obj <-sensitivity::tell(sobol_obj, y)
-
-  # Return the updated object
-  return(invisible(sobol_obj))
+                                  var = "GPP") {
+  if (is.null(sobol_obj$backend) || sobol_obj$backend != "sensobol") {
+    PEcAn.logger::logger.error(
+      "compute_sobol_indices expects a sensobol design object returned by ",
+      "generate_joint_ensemble_design(..., sobol = TRUE)"
+    )
+  }
+
+  output_files <- list.files(
+    outdir,
+    pattern = "^ensemble\\.output\\..*\\.Rdata$",
+    full.names = TRUE
+  )
+  if (length(output_files) == 0) {
+    PEcAn.logger::logger.error("No ensemble.output.*.Rdata files found in ", outdir)
+  }
+
+  output_var <- vapply(
+    strsplit(basename(output_files), "\\."),
+    function(x) if (length(x) >= 4) x[[4]] else NA_character_,
+    character(1)
+  )
+  matched_files <- output_files[output_var == var]
+
+  if (length(matched_files) == 0) {
+    PEcAn.logger::logger.error(
+      "No standardized ensemble output found for variable '", var,
+      "' in ", outdir
+    )
+  }
+  if (length(matched_files) > 1) {
+    PEcAn.logger::logger.error(
+      "Multiple standardized ensemble outputs found for variable '", var,
+      "' in ", outdir, ". Please keep only one matching file."
+    )
+  }
+
+  output_file <- matched_files[[1]]
+  output_env <- new.env(parent = emptyenv())
+  load(output_file, envir = output_env)
+  if (is.null(output_env$ensemble.output)) {
+    PEcAn.logger::logger.error(
+      "Object `ensemble.output` missing from standardized output file ",
+      output_file
+    )
+  }
+
+  y <- as.numeric(unlist(output_env$ensemble.output, use.names = FALSE))
+  expected_length <- sobol_obj$N * (length(sobol_obj$params) + 2L)
+  if (length(y) != expected_length) {
+    PEcAn.logger::logger.error(
+      "Standardized ensemble output has ", length(y),
+      " values but expected ", expected_length,
+      " for Sobol design size N = ", sobol_obj$N
+    )
+  }
+
+  sobol_indices_result <- sensobol::sobol_indices(
+    matrices = sobol_obj$matrices,
+    Y = y,
+    N = sobol_obj$N,
+    params = sobol_obj$params,
+    first = sobol_obj$first,
+    total = sobol_obj$total,
+    order = "first",
+    boot = FALSE
+  )
+
+  sobol_results <- tibble::as_tibble(sobol_indices_result$results)
+  if (!is.null(sobol_obj$factor_metadata)) {
+    factor_metadata <- tibble::as_tibble(sobol_obj$factor_metadata) |>
+      dplyr::rename(parameters = .data$factor)
+    sobol_results <- dplyr::left_join(
+      sobol_results,
+      factor_metadata,
+      by = "parameters"
+    )
+  }
+
+  sobol_design <- sobol_obj
+  save(
+    sobol_design,
+    file = file.path(
+      outdir,
+      sub("^ensemble\\.output\\.", "sobol.design.", basename(output_file))
+    )
+  )
+  save(
+    sobol_indices_result,
+    sobol_results,
+    file = file.path(
+      outdir,
+      sub("^ensemble\\.output\\.", "sobol.indices.", basename(output_file))
+    )
+  )
+
+  return(sobol_results)
 }

modules/uncertainty/R/ensemble.R

@@ -598,21 +598,20 @@ input.ens.gen <- function(settings, ensemble_size, input, method = "sampling", p
   if (input == "parameters") return(NULL)
 
   input_path <- settings$run$inputs[[tolower(input)]]$path
-  if (is.null(input_path)) {
-    PEcAn.logger::logger.error(
-      "No paths found for input", sQuote(input), "in settings$run$inputs"
-    )
+  if (is.null(input_path) || length(input_path) == 0) {
+    PEcAn.logger::logger.error("Input ", sQuote(input), " has no paths specified")
   }
 
   if (!is.null(parent_ids)) {
     samples$ids <- parent_ids$ids
-    out.of.sample.size <- length(samples$ids[samples$ids > length(input_path)])
-    # sample for those that our outside the param size -
-    # for example, parent id may send id number 200 but we have only 100 sample for param
-    samples$ids[samples$ids %in% out.of.sample.size] <- sample(
-      seq_along(input_path),
-      out.of.sample.size,
-      replace = TRUE)
+    overflow_positions <- which(samples$ids > length(input_path))
+    if (length(overflow_positions) > 0) {
+      samples$ids[overflow_positions] <- sample(
+        seq_along(input_path),
+        length(overflow_positions),
+        replace = TRUE
+      )
+    }
   } else if (tolower(method) == "sampling") {
     samples$ids <- sample(
       seq_along(input_path),
@@ -623,6 +622,14 @@ input.ens.gen <- function(settings, ensemble_size, input, method = "sampling", p
       seq_along(input_path),
       length.out = ensemble_size)
   }
+  samples$ids <- as.integer(samples$ids)
+  if (any(is.na(samples$ids)) ||
+      any(samples$ids < 1L) ||
+      any(samples$ids > length(input_path))) {
+    PEcAn.logger::logger.error(
+      "Invalid sampled ids generated for input ", sQuote(input)
+    )
+  }
   #using the sample ids
   samples$samples <- input_path[samples$ids]