fix(uncertainty): resolve SA median manifest lookup (#3882)

docker/depends/pecan_package_dependencies.csv

@@ -252,6 +252,7 @@
 "ncdf4","*","models/stics","Imports",FALSE
 "ncdf4","*","modules/assim.sequential","Imports",FALSE
 "ncdf4","*","modules/data.remote","Imports",FALSE
+"ncdf4","*","modules/uncertainty","Suggests",FALSE
 "ncdf4",">= 1.15","base/utils","Imports",FALSE
 "ncdf4",">= 1.15","base/visualization","Imports",FALSE
 "ncdf4",">= 1.15","models/biocro","Imports",FALSE

modules/uncertainty/DESCRIPTION

@@ -47,6 +47,7 @@ Imports:
     tidyr
 Suggests:
     mockery,
+    ncdf4,
     testthat (>= 1.0.2),
     withr
 License: BSD_3_clause + file LICENSE

modules/uncertainty/NEWS.md

@@ -10,6 +10,13 @@
 sensitivity analysis. This function ensures non-parameter inputs
 (met, IC, soil) are held constant, which is required for valid 
 variance decomposition in one-at-a-time (OAT) sensitivity analysis (#3729)
+* Fixed `read.sa.output()` dropping the median (q50) row when using the
+  manifest-based OAT sensitivity path. `write.sa.configs()` writes a single
+  shared median run with `pft_name = "NA"` and `trait = "NA"`, but the lookup
+  previously required an exact per-trait/per-PFT match, so the cell was always
+  `NA`. The fix adds a fallback: when `quantile == "50"` and no exact match is
+  found, the function now resolves the shared median entry before emitting a
+  warning.
 * `get.parameter.samples()` now consistently accepts PFTs with no `outdir` specified. These previously failed when no database connection was available.
 * Improved row alignment and text scaling in `plot_variance_decomposition()`.
 * `plot_variance_decomposition()` gains argument `order_by`, sorting highest variance to the top by default. To reproduce old behavior, use `order_by = "none"`.

modules/uncertainty/R/sensitivity.R

@@ -38,18 +38,42 @@ read.sa.output <- function(traits, quantiles, pecandir, outdir, pft.name = "",
   for (trait in traits) {
     for (quantile in quantiles) {
       # We look for the row that matches the current pft, trait, and quantile.
+      # Use which() to handle NA values in manifest columns (read.csv converts
+      # the string "NA" to actual R NA, and comparison with NA yields NA, not FALSE).
       subset_df <- manifest[
-        manifest$type == "Sensitivity" & 
-        manifest$pft_name == pft.name &
-        manifest$trait == trait &
-        as.character(manifest$quantile) == as.character(quantile), 
+        which(
+          manifest$type == "Sensitivity" & 
+          manifest$pft_name == pft.name &
+          manifest$trait == trait &
+          as.character(manifest$quantile) == as.character(quantile)
+        ), 
       ]
 
       if (nrow(subset_df) == 1) {
          run.id <- subset_df$run_id
       } else if (nrow(subset_df) > 1) {
          PEcAn.logger::logger.warn("Multiple runs found for", trait, quantile, "- using the last one.")
          run.id <- utils::tail(subset_df$run_id, 1)
+      } else if (as.character(quantile) == "50") {
+        # The median run is written as a single shared entry with pft_name = "NA"
+        # and trait = "NA". Fall back to that shared row when no exact match exists.
+        # Note: read.csv() converts the string "NA" to actual R NA values,
+        # so we must use is.na() rather than == "NA" for the comparison.
+        median_df <- manifest[
+          manifest$type == "Sensitivity" &
+          is.na(manifest$pft_name) &
+          is.na(manifest$trait) &
+          as.character(manifest$quantile) == "50",
+        ]
+        if (nrow(median_df) >= 1) {
+          PEcAn.logger::logger.debug(
+            "Using shared median run for", trait, "quantile 50"
+          )
+          run.id <- utils::tail(median_df$run_id, 1)
+        } else {
+          PEcAn.logger::logger.warn("No run found in manifest for", trait, quantile)
+          next # Skip this quantile
+        }
       } else {
          PEcAn.logger::logger.warn("No run found in manifest for", trait, quantile)
          next # Skip this quantile

modules/uncertainty/tests/testthat/test-read_sa_output_median.R

@@ -0,0 +1,104 @@
+# Regression test for the shared-median manifest fallback in read.sa.output()
+#
+# Background: write.sa.configs() writes the median (q50) run as a single shared
+# entry with pft_name = "NA" and trait = "NA".  read.sa.output() previously did
+# a strict per-trait/per-PFT lookup, so the median row always failed to match
+# and the 50th-quantile cell was returned as NA.
+#
+# The fix adds a fallback: when quantile == "50" and no exact match is found,
+# look for the shared median entry before giving up.
+#
+# Reference: https://github.com/pecanproject/pecan/issues/3882
+
+setup_logger_for_testing <- function() {
+  PEcAn.logger::logger.setUseConsole(TRUE, FALSE)
+  PEcAn.logger::logger.setLevel("WARN")
+}
+
+test_that("read.sa.output resolves median (q50) via shared manifest fallback", {
+  skip_if_not_installed("ncdf4")
+
+  withr::with_tempdir({
+    setup_logger_for_testing()
+    on.exit(PEcAn.logger::logger.setUseConsole(TRUE, TRUE), add = TRUE)
+
+    pft   <- "temperate.coniferous"
+    trait <- "growth_resp_factor"
+    yr    <- 2004
+
+    # One shared median row (pft_name = "NA", trait = "NA") + no trait-specific row
+    # for quantile 50.  This is the exact situation that triggered the original bug.
+    median_run_id <- "SA-median--1"
+    manifest <- data.frame(
+      type     = "Sensitivity",
+      pft_name = "NA",
+      trait    = "NA",
+      quantile = "50",
+      run_id   = median_run_id,
+      stringsAsFactors = FALSE
+    )
+    write.csv(manifest, "runs_manifest.csv", row.names = FALSE)
+
+    run_outdir <- file.path(getwd(), median_run_id)
+    dir.create(run_outdir, recursive = TRUE)
+
+    nc_path <- file.path(run_outdir, paste0(yr, ".nc"))
+    nc_obj <- ncdf4::nc_create(
+      nc_path,
+      list(ncdf4::ncvar_def("NPP", "kg m-2 s-1", list(), missval = NA_real_))
+    )
+    ncdf4::ncvar_put(nc_obj, "NPP", 1.23)
+    ncdf4::nc_close(nc_obj)
+
+    out <- PEcAn.uncertainty::read.sa.output(
+      traits     = trait,
+      quantiles  = "50",
+      pecandir   = getwd(),
+      outdir     = getwd(),
+      pft.name   = pft,
+      start.year = yr,
+      end.year   = yr,
+      variable   = PEcAn.utils::convert.expr("NPP")$variable.eqn
+    )
+
+    expect_false(
+      is.na(out[["50", trait]]),
+      label = "median (q50) output should not be NA"
+    )
+  })
+})
+
+
+test_that("read.sa.output still warns when no median fallback row exists", {
+  withr::with_tempdir({
+    setup_logger_for_testing()
+    on.exit(PEcAn.logger::logger.setUseConsole(TRUE, TRUE), add = TRUE)
+
+    # q50 row with concrete pft_name and trait: no exact match and no NA/NA fallback
+    manifest <- data.frame(
+      type     = "Sensitivity",
+      pft_name = "temperate.coniferous",
+      trait    = "other_trait",
+      quantile = "50",
+      run_id   = "SA-other--1",
+      stringsAsFactors = FALSE
+    )
+    write.csv(manifest, "runs_manifest.csv", row.names = FALSE)
+
+    out <- expect_output(
+      PEcAn.uncertainty::read.sa.output(
+        traits     = "growth_resp_factor",
+        quantiles  = "50",
+        pecandir   = getwd(),
+        outdir     = getwd(),
+        pft.name   = "temperate.coniferous",
+        start.year = 2004,
+        end.year   = 2004,
+        variable   = PEcAn.utils::convert.expr("NPP")$variable.eqn
+      ),
+      regexp = "No run found in manifest"
+    )
+
+    expect_true(is.na(out[["50", "growth_resp_factor"]]))
+  })
+})