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

@@ -70,9 +70,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."
       )
     }
   }
@@ -171,6 +171,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]]
@@ -182,6 +202,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"
+    )
+  })
+})

docker/depends/Dockerfile

@@ -23,6 +23,7 @@ RUN apt-get update \
        librdf0-dev \
        libnetcdf-dev \
        libudunits2-dev \
+       libgsl-dev \
        libgl1-mesa-dev \
        libglu1-mesa-dev \
     && rm -rf /var/lib/apt/lists/*

docker/depends/pecan_package_dependencies.csv

@@ -585,7 +585,7 @@
 "RPostgreSQL","*","models/biocro","Suggests",FALSE
 "RSQLite","*","base/db","Suggests",FALSE
 "rstac","*","models/peprmt","Suggests",FALSE
-"sensitivity","*","modules/uncertainty","Imports",FALSE
+"sensobol","*","modules/uncertainty","Imports",FALSE
 "sessioninfo","*","base/all","Suggests",FALSE
 "sf","*","modules/assim.sequential","Suggests",FALSE
 "sf","*","modules/data.atmosphere","Imports",FALSE
@@ -674,6 +674,7 @@
 "tibble","*","modules/data.land","Imports",FALSE
 "tibble","*","modules/data.remote","Suggests",FALSE
 "tibble","*","modules/meta.analysis","Suggests",FALSE
+"tibble","*","modules/uncertainty","Imports",FALSE
 "tictoc","*","modules/assim.sequential","Suggests",FALSE
 "tidyr","*","base/db","Imports",FALSE
 "tidyr","*","models/ed","Imports",FALSE

modules/uncertainty/DESCRIPTION

@@ -43,7 +43,8 @@ Imports:
     purrr,
     randtoolbox,
     rlang,
-    sensitivity,
+    sensobol,
+    tibble,
     tidyr
 Suggests:
     mockery,

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