Merge pull request #3535 from blesson07asd/blesson-GSoC

joint_ensemble sampling

Author: mdietze

base/workflow/NEWS.md

@@ -1,3 +1,14 @@
+### Joint Ensemble Sampling Implementation
+- **Behavior Change**: Ensemble runs now use shared input samples across all sites instead of independent sampling per site
+- **Affected Components**: 
+  - `get.parameter.samples.R` 
+  - `runModule.run.write.configs.R`
+  - `run.write.configs.R` 
+  - `ensemble.R`
+- **New Default**: The sampling method default has changed from "uniform" to "random"
+- **Impact**: This ensures consistency across sites in ensemble runs but may produce different results compared to previous versions
+
+
 # PEcAn.workflow 1.9.0
 
 * PEcAn.workflow is now distributed under the BSD 3-clause license instead of the NCSA Open Source license.

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

@@ -5,9 +5,10 @@
 #' \code{write.config.*} function for your specific ecosystem model
 #' (e.g. write.config.ED2, write.config.SIPNET).
 #'
+#'
 #' @param settings a PEcAn settings list
 #' @param write should the runs be written to the database?
-#' @param ens.sample.method how to sample the ensemble members('halton' sequence or 'uniform' random)
+#' @param input_design input indices for samples 
 #' @param posterior.files Filenames for posteriors for drawing samples for ensemble and sensitivity
 #'    analysis (e.g. post.distns.Rdata, or prior.distns.Rdata)
 #' @param overwrite logical: Replace output files that already exist?
@@ -22,9 +23,13 @@
 #' @export
 #'
 #' @author David LeBauer, Shawn Serbin, Ryan Kelly, Mike Dietze
-run.write.configs <- function(settings, write = TRUE, ens.sample.method = "uniform", 
+
+run.write.configs <- function(settings, input_design, write = TRUE,  
                               posterior.files = rep(NA, length(settings$pfts)), 
                               overwrite = TRUE) {
+  
+
+
   ## Skip database connection if settings$database is NULL or write is False
   if (!isTRUE(write) && is.null(settings$database)) {
     PEcAn.logger::logger.info("Not writing this run to database, so database connection skipped")
@@ -93,13 +98,23 @@ run.write.configs <- function(settings, write = TRUE, ens.sample.method = "unifo
   scipen <- getOption("scipen")
   options(scipen = 12)
 
-  PEcAn.uncertainty::get.parameter.samples(settings, posterior.files, ens.sample.method)
   samples.file <- file.path(settings$outdir, "samples.Rdata")
   if (file.exists(samples.file)) {
     samples <- new.env()
     load(samples.file, envir = samples) ## loads ensemble.samples, trait.samples, sa.samples, runs.samples, env.samples
     trait.samples <- samples$trait.samples
-    ensemble.samples <- samples$ensemble.samples
+    trait_sample_indices <- input_design[["param"]]
+    ensemble.samples <- list()
+    for (pft in names(trait.samples)) {
+        pft_traits <- trait.samples[[pft]]
+        ensemble.samples[[pft]] <- as.data.frame(
+          lapply(
+            names(pft_traits),
+            function(trait) pft_traits[[trait]][trait_sample_indices]
+          )
+        )
+      names(ensemble.samples[[pft]]) <- names(pft_traits)
+    }
     sa.samples <- samples$sa.samples
     runs.samples <- samples$runs.samples
     ## env.samples <- samples$env.samples
@@ -167,6 +182,7 @@ run.write.configs <- function(settings, write = TRUE, ens.sample.method = "unifo
                                                           ensemble.samples = ensemble.samples, 
                                                           settings = settings,
                                                           model = model, 
+                                                          input_design = input_design,
                                                           write.to.db = write)
 
     # Store output in settings and output variables
@@ -190,6 +206,8 @@ run.write.configs <- function(settings, write = TRUE, ens.sample.method = "unifo
        file = file.path(settings$outdir, "samples.Rdata"))
   PEcAn.logger::logger.info("parameter values for runs in ", file.path(settings$outdir, "samples.RData"))
   options(scipen = scipen)
-  
-  return(invisible(settings))
+  invisible(settings)
+  return(settings)
 }
+
+  

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

@@ -2,22 +2,34 @@
 #'
 #' @param settings a PEcAn Settings or MultiSettings object
 #' @param overwrite logical: Replace config files if they already exist?
+#' @param input_design the input indices for samples 
 #' @return A modified settings object, invisibly
 #' @importFrom dplyr %>%
 #' @export
-runModule.run.write.configs <- function(settings, overwrite = TRUE) {
 
+
+runModule.run.write.configs <- function(settings, overwrite = TRUE, input_design = NULL) {
+  
   if (PEcAn.settings::is.MultiSettings(settings)) {
     if (overwrite && file.exists(file.path(settings$rundir, "runs.txt"))) {
       PEcAn.logger::logger.warn("Existing runs.txt file will be removed.")
       unlink(file.path(settings$rundir, "runs.txt"))
     }
-    return(PEcAn.settings::papply(settings, runModule.run.write.configs, overwrite = FALSE))
+    if (is.null(input_design)) {
+      ensemble_size <- settings$ensemble$size
+      input_design <- PEcAn.uncertainty::generate_joint_ensemble_design(settings=settings[1],ensemble_size=ensemble_size)
+      }
+    return(PEcAn.settings::papply(settings, runModule.run.write.configs, overwrite = FALSE,input_design=input_design))
   } else if (PEcAn.settings::is.Settings(settings)) {
     # double check making sure we have method for parameter sampling
     if (is.null(settings$ensemble$samplingspace$parameters$method)) {
       settings$ensemble$samplingspace$parameters$method <- "uniform"
     }
+      if (is.null(input_design)) {
+      ensemble_size <- settings$ensemble$size
+      input_design <- PEcAn.uncertainty::generate_joint_ensemble_design( settings = settings, ensemble_size = ensemble_size )
+    }
+    
 
 
     #check to see if there are posterior.files tags under pft
@@ -27,9 +39,9 @@ runModule.run.write.configs <- function(settings, overwrite = TRUE) {
     return(PEcAn.workflow::run.write.configs(
       settings = settings,
       write = isTRUE(settings$database$bety$write), # treat null as FALSE
-      ens.sample.method = settings$ensemble$samplingspace$parameters$method,
       posterior.files = posterior.files,
-      overwrite = overwrite
+      overwrite = overwrite,
+      input_design = input_design
     ))
   } else {
     stop("runModule.run.write.configs only works with Settings or MultiSettings")

base/workflow/man/run.write.configs.Rd

@@ -580,6 +580,9 @@ Tags related to ensemble setup are:
 * `size` : (required) the number of runs in the ensemble.
 * `samplingspace`: (optional) Contains tags for defining how the ensembles will be generated.
 
+Shared sampling design: In multi-site workflows, PEcAn now generates one joint ensemble design and reuses the same input sample indices (parameters, meteorology, etc.) across all sites to ensure consistent draws; previously, inputs were sampled independently per site. This change does not introduce new XML tags and applies whenever multiple sites are processed together.
+The joint sampling design is created once at the start of configuration using generate_joint_ensemble_design(), called from the run configuration module, and the resulting indices are threaded through to write.ensemble.configs().
+
 Each piece in the sampling space can potentially have a method tag and a parent tag. Method refers to the sampling method and parent refers to the cases where we need to link the samples of two components. When no tag is defined for one component, one sample will be generated and used for all the ensembles. This allows for partitioning/studying different sources of uncertainties. For example, if no met tag is defined then, one met path will be used for all the ensembles and as a result the output uncertainty will come from the variability in the parameters. At the moment no sampling method is implemented for soil and vegetation.
 Available sampling methods for `parameters` can be found in the documentation of the `PEcAn.utils::get.ensemble.samples` function.
 For the cases where we need simulations with a predefined set of parameters, met and initial condition we can use the restart argument. Restart needs to be a list with name tags of `runid`, `inputs`, `new.params` (parameters), `new.state` (initial condition), `ensemble.id` (ensemble ids), `start.time`, and `stop.time`.

base/workflow/man/runModule.run.write.configs.Rd

@@ -682,11 +682,12 @@ write.config.SIPNET <- function(defaults, trait.values, settings, run.id, inputs
         soilWFrac <- ncdf4::ncvar_get(IC.nc, "SoilMoistFrac")
         if (!is.na(soilWFrac) && is.numeric(soilWFrac)) {
           param[param[, 1] == "soilWFracInit", 2] <- sum(soilWFrac) / 100
+          ## litterWFracInit fraction
+          litterWFrac <- soilWFrac
         }
       }
-      ## litterWFracInit fraction
-      litterWFrac <- soilWFrac
-
+      
+           
       ## snowInit cm water equivalent (cm = g / cm2 because 1 g water = 1 cm3 water)
       if (ic_has_ncvars[["SWE"]]) {
         snow <- ncdf4::ncvar_get(IC.nc, "SWE")

book_source/03_topical_pages/03_pecan_xml.Rmd

@@ -3,6 +3,7 @@
 export(ensemble.filename)
 export(ensemble.ts)
 export(flux.uncertainty)
+export(generate_joint_ensemble_design)
 export(get.change)
 export(get.coef.var)
 export(get.elasticity)

models/sipnet/R/write.configs.SIPNET.R

@@ -80,11 +80,12 @@ read.ensemble.output <- function(ensemble.size, pecandir, outdir, start.year, en
 ##' @export
 ##' @author David LeBauer, Istem Fer
 get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples, 
-                                 method = "uniform", param.names = NULL, ...) {
-  
-  if (is.null(method)) {
-    PEcAn.logger::logger.info("No sampling method supplied, defaulting to uniform random sampling")
-    method <- "uniform"
+                                 method = "random", param.names = NULL, ...) {
+
+  # Define supported methods
+  supported_methods <- c("random", "uniform", "halton", "sobol", "lhc")
+  if (!method %in% supported_methods) {
+    stop("Invalid sampling method")
   }
 
   ## force as numeric for compatibility with Fortran code in halton()
@@ -139,48 +140,59 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 
 
     ensemble.samples <- list()
-    
+    sampled.indices <- list()
 
     col.i <- 0
     for (pft.i in seq(pft.samples)) {
       ensemble.samples[[pft.i]] <- matrix(nrow = ensemble.size, ncol = length(pft.samples[[pft.i]]))
+      sampled.indices[[pft.i]] <- matrix(nrow = ensemble.size, ncol = length(pft.samples[[pft.i]]))
 
       # meaning we want to keep MCMC samples together
       if(length(pft.samples[[pft.i]])>0 & !is.null(param.names)){ 
         if (method == "halton") {
-          same.i <- round(randtoolbox::halton(ensemble.size) * length(pft.samples[[pft.i]][[1]]))
+          same.i <- floor(randtoolbox::halton(ensemble.size) * length(pft.samples[[pft.i]][[1]]))+1
         } else if (method == "sobol") {
-          same.i <- round(randtoolbox::sobol(ensemble.size, scrambling = 3) * length(pft.samples[[pft.i]][[1]]))
+          same.i <- floor(randtoolbox::sobol(ensemble.size, scrambling = 3) * length(pft.samples[[pft.i]][[1]]))+1
         } else if (method == "torus") {
-          same.i <- round(randtoolbox::torus(ensemble.size) * length(pft.samples[[pft.i]][[1]]))
+          same.i <- floor(randtoolbox::torus(ensemble.size) * length(pft.samples[[pft.i]][[1]]))+1
         } else if (method == "lhc") {
-          same.i <- round(c(PEcAn.emulator::lhc(t(matrix(0:1, ncol = 1, nrow = 2)), ensemble.size) * length(pft.samples[[pft.i]][[1]])))
+          same.i <- floor(c(PEcAn.emulator::lhc(t(matrix(0:1, ncol = 1, nrow = 2)), ensemble.size) * length(pft.samples[[pft.i]][[1]])))+1
         } else if (method == "uniform") {
           same.i <- sample.int(length(pft.samples[[pft.i]][[1]]), ensemble.size)
-        } else {
-          PEcAn.logger::logger.info("Method ", method, " has not been implemented yet, using uniform random sampling")
-          # uniform random
-          same.i <- sample.int(length(pft.samples[[pft.i]][[1]]), ensemble.size)
+        } else if (method == "random") {
+            PEcAn.logger::logger.info("Using random row sampling for MCMC draws")
+           same.i <- sample(nrow(pft.samples[[pft.i]][[1]]), ensemble.size, replace = TRUE)
+        }
+        else {
+          PEcAn.logger::logger.error("Sampling method %s is not recognized", method)
+        
         }
 
       }
 
       for (trait.i in seq(pft.samples[[pft.i]])) {
         col.i <- col.i + 1
-        if(names(pft.samples[[pft.i]])[trait.i] %in% param.names[[pft.i]]){ # keeping samples
-          ensemble.samples[[pft.i]][, trait.i] <- pft.samples[[pft.i]][[trait.i]][same.i]
-        }else{
+        if (names(pft.samples[[pft.i]])[trait.i] %in% param.names[[pft.i]]) { 
+             ensemble.samples[[pft.i]][, trait.i] <- pft.samples[[pft.i]][[trait.i]][same.i]
+             sampled.indices[[pft.i]][, trait.i] <- same.i
+       }else{
+          # Extract original trait values
+          trait.values <- pft.samples[[pft.i]][[trait.i]]
+          sampled.values <- stats::quantile(trait.values, random.samples[, col.i])
+
           ensemble.samples[[pft.i]][, trait.i] <- stats::quantile(pft.samples[[pft.i]][[trait.i]],
                                                                   random.samples[, col.i])
-        }
-      }  # end trait
-      ensemble.samples[[pft.i]] <- as.data.frame(ensemble.samples[[pft.i]])
-      colnames(ensemble.samples[[pft.i]]) <- names(pft.samples[[pft.i]])
-    }  #end pft
-    names(ensemble.samples) <- names(pft.samples)
-    ans <- ensemble.samples
+          sampled.indices[[pft.i]][, trait.i] <- sapply(sampled.values, function(val) {which.min(abs(trait.values - val)) })
+      }   
+    }  
+          ensemble.samples[[pft.i]] <- as.data.frame(ensemble.samples[[pft.i]])
+          colnames(ensemble.samples[[pft.i]]) <- names(pft.samples[[pft.i]])
+    
+  }  #end pft
+   names(ensemble.samples) <- names(pft.samples)
+   ans <- ensemble.samples
   }
-  return(ans)
+    return(list(ans,sampled.indices))
 } # get.ensemble.samples
 
 
@@ -190,6 +202,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' Given a pft.xml object, a list of lists as supplied by get.sa.samples, 
 ##' a name to distinguish the output files, and the directory to place the files.
 ##'
+##' @param input_design the input indices for samples 
 ##' @param defaults pft
 ##' @param ensemble.samples list of lists supplied by \link{get.ensemble.samples}
 ##' @param settings list of PEcAn settings
@@ -211,7 +224,7 @@ get.ensemble.samples <- function(ensemble.size, pft.samples, env.samples,
 ##' @importFrom rlang .data
 ##' @export
 ##' @author David LeBauer, Carl Davidson, Hamze Dokoohaki
-write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, 
+write.ensemble.configs <- function(defaults, ensemble.samples, settings, model, input_design ,
                                    clean = FALSE, write.to.db = TRUE, restart = NULL, samples = NULL, rename = FALSE) {
 
 
@@ -314,18 +327,21 @@ for (input_tag in names(settings$run$inputs)) {
     samp.ordered <- samp[c(order, names(samp)[!(names(samp) %in% order)])]
     if(is.null(samples)){
        #performing the sampling
-       samples<-list()
-       # For the tags specified in the xml I do the sampling
-      for(i in seq_along(samp.ordered)){
-         myparent<-samp.ordered[[i]]$parent # do I have a parent ?
-         #call the function responsible for generating the ensemble
-         samples[[names(samp.ordered[i])]] <- input.ens.gen(settings=settings,
-                                                         input=names(samp.ordered)[i],
-                                                         method=samp.ordered[[i]]$method,
-                                                         parent_ids=if( !is.null(myparent)) samples[[myparent]]) # if I have parent then give me their ids - this is where the ordering matters making sure the parent is done before it's asked
-       }
+      samples <- list()
+      input_tags <- names(settings$run$inputs)
+
+      for (input_tag in input_tags) {
+           if (input_tag %in% colnames(input_design)) {
+                  input_paths <- settings$run$inputs[[input_tag]]$path
+                  input_indices <- input_design[[input_tag]]
+
+                 samples[[input_tag]] <- list(
+                   samples = lapply(input_indices, function(idx) input_paths[[idx]])
+                 )
+    }
+
      }
-    
+    }
     # if there is a tag required by the model but it is not specified in the xml then I replicate n times the first element 
     required_tags%>%
       purrr::walk(function(r_tag){