22# '
33# ' Creates an input design matrix for sensitivity analysis where non-parameter
44# ' inputs (met, IC, soil, etc.) are held constant while parameters vary
5- # ' one-at-a-time across quantiles. This differs from ensemble design where
6- # ' all inputs vary together.
5+ # ' one-at-a-time across quantiles. This differs from ensemble design, where all
6+ # ' inputs vary together.
77# '
8- # ' @details
9- # ' ## Settings requirements
10- # '
11- # ' This function directly uses:
12- # ' \itemize{
13- # ' \item \code{settings$outdir} - Output directory path for samples.Rdata
14- # ' \item \code{settings$pfts} - List of PFTs (extracts \code{posterior.files})
15- # ' \item \code{settings$ensemble$samplingspace} - Input types to include in design
16- # ' }
17- # '
18- # ' When \code{sa_samples = NULL}, settings is passed to
19- # ' \code{\link{get.parameter.samples}} which additionally requires:
20- # ' \itemize{
21- # ' \item \code{settings$sensitivity.analysis} - SA quantile configuration
22- # ' \item \code{settings$database$bety} - Database connection (optional)
23- # ' \item \code{settings$host$name} - Host name for dbfile.check (optional)
24- # ' }
8+ # ' Parameter samples are drawn in memory via \code{\link{load_pft_posteriors}}
9+ # ' and \code{\link{get_parameter_samples}} (mirroring
10+ # ' \code{\link{generate_joint_ensemble_design}}), or reused when a \code{samples}
11+ # ' bundle is supplied. The design is built from the quantile-based
12+ # ' \code{sa.samples} within that bundle.
2513# '
26- # ' ## OAT design logic
27- # ' For sensitivity analysis, we must isolate the effect of each
28- # ' parameter by holding all other inputs constant. The param column contains
29- # ' sequential indices (1, 2, 3, ...) matching the SA run order in
30- # ' \code{write.sa.configs}. All other columns (met, ic, soil, etc.) are set to 1,
31- # ' meaning the first input file is always used.
14+ # ' @param settings PEcAn settings object. Uses \code{settings$pfts},
15+ # ' \code{settings$sensitivity.analysis$quantiles} (the SA quantiles, when
16+ # ' sampling here), and \code{settings$ensemble$samplingspace} (the input types
17+ # ' that form the design columns).
18+ # ' @param samples Optional pre-computed parameter samples (a list containing at
19+ # ' least \code{sa.samples}, as returned by \code{\link{get_parameter_samples}}).
20+ # ' When supplied these are used directly; when \code{NULL} (default) they are
21+ # ' sampled in memory.
3222# '
33- # ' Note on internal dependencies
23+ # ' @return A list with \code{X}, a data.frame with one row per SA run and one
24+ # ' column per input type (the \code{param} column holds sequential run indices,
25+ # ' every other column is held at 1), and \code{samples}, the parameter bundle
26+ # ' used.
3427# '
35- # ' If sa_samples is NULL we hand off to get.parameter.samples(), which does
36- # ' the work of finding and loading parameter distributions.
37- # '
38- # ' In practice it:
39- # ' - uses pft$posterior.files directly when it is defined (an Rdata file with
40- # ' post.distns or prior.distns),
41- # ' - otherwise figures out an output directory from pft$outdir or, if needed,
42- # ' via pft$posteriorid in the database,
43- # ' - then looks in that directory for post.distns.Rdata, falling back to
44- # ' prior.distns.Rdata,
45- # ' - and, for MCMC posteriors, looks up trait.mcmc*.Rdata linked to the same
46- # ' posteriorid or a trait.mcmc.Rdata file in that directory.
47- # '
48- # ' @param settings PEcAn settings object. See details for required elements.
49- # ' @param sa_samples Optional. Pre-loaded SA samples (named list with one
50- # ' element per PFT, each a matrix with quantiles as rows and traits as columns).
51- # ' If NULL (default), samples are generated via \code{get.parameter.samples}.
52- # '
53- # ' @return list with component X: a data.frame with columns for each input type
54- # ' and one row per SA run. Non-parameter columns are all 1 (constant).
55- # '
56- # ' @examples
57- # ' \dontrun{
58- # ' # Generate SA design for a multi-site run
59- # ' sa_design <- generate_OAT_SA_design(settings)
60- # '
61- # ' # View the design matrix
62- # ' print(sa_design$X)
63- # ' # param met ic soil
64- # ' # 1 1 1 1 1 # Median run
65- # ' # 2 2 1 1 1 # trait1 @ q=2.3%
66- # ' # 3 3 1 1 1 # trait1 @ q=15.9%
67- # ' # 4 4 1 1 1 # trait1 @ q=84.1%
68- # ' # ...
69- # '
70- # ' # With pre-loaded sa_samples (skips get.parameter.samples call)
71- # ' load("samples.Rdata")
72- # ' sa_design <- generate_OAT_SA_design(settings, sa_samples = sa.samples)
73- # ' }
74- # ' @export
75- # ' @author Akash B V
28+ # ' @author Akash B V, Om Kapale
7629# ' @importFrom rlang %||%
30+ # ' @export
31+ generate_OAT_SA_design <- function (settings , samples = NULL ) {
7732
78- generate_OAT_SA_design <- function (settings , sa_samples = NULL ) {
79-
80- samples_file <- file.path(settings $ outdir , " samples.Rdata"
81-
82- if (is.null(sa_samples )) {
83-
33+ # Generate parameter samples in memory (or use the ones passed in), mirroring
34+ # generate_joint_ensemble_design. A sensitivity analysis needs the quantile-
35+ # based sa.samples, so we request those and skip the ensemble draw.
36+ if (is.null(samples )) {
8437 posterior.files <- settings $ pfts %> %
8538 purrr :: map_chr(" posterior.files" .default = NA_character_ )
86-
87- # generate parameter samples - sa.samples created from quantiles
88- PEcAn.uncertainty :: get.parameter.samples(
89- settings ,
90- posterior.files = posterior.files
39+ loaded <- load_pft_posteriors(settings , posterior.files )
40+ samples <- get_parameter_samples(
41+ pft_names = loaded $ pft_names ,
42+ prior_distns_list = loaded $ prior_distns_list ,
43+ trait_mcmc_list = loaded $ trait_mcmc_list ,
44+ ensemble.size = settings $ ensemble $ size %|| % 1 ,
45+ ens.sample.method = settings $ ensemble $ samplingspace $ parameters $ method %|| % " uniform"
46+ sa_quantiles = settings $ sensitivity.analysis $ quantiles ,
47+ do_ensemble = FALSE ,
48+ independent = loaded $ independent
9149 )
50+ }
9251
93- samples_env <- new.env()
94- load(samples_file , envir = samples_env )
95- sa_samples <- samples_env $ sa.samples
96-
97- if (is.null(sa_samples )) {
98- PEcAn.logger :: logger.severe(
99- " sa.samples not found in samples.Rdata."
100- " Ensure sensitivity.analysis is configured in settings."
101- )
102- }
52+ sa_samples <- samples $ sa.samples
53+
54+ if (is.null(sa_samples ) || length(sa_samples ) == 0 ) {
55+ PEcAn.logger :: logger.severe(
56+ " sa.samples are empty."
57+ " Ensure sensitivity.analysis quantiles are configured in settings."
58+ )
10359 }
104-
105- # calculate total number of SA runs
106- # 1 median + ( traits * non-median quantiles) per PFT
60+
61+ # Total number of SA runs: 1 median run plus, for each PFT,
62+ # (n traits) * (n non-median quantiles).
10763 MEDIAN <- " 50"
108- num_sa_runs <- 1 # start with median run
109-
64+ num_sa_runs <- 1
65+
11066 for (pft_name in names(sa_samples )) {
11167 if (pft_name == " env" next
112-
68+
11369 pft_samples <- sa_samples [[pft_name ]]
11470 n_traits <- ncol(pft_samples )
11571 quantile_names <- rownames(pft_samples )
11672 n_non_median <- sum(quantile_names != MEDIAN )
117-
118- # add runs for this pft: (traits) * (non-median quantiles)
73+
11974 num_sa_runs <- num_sa_runs + (n_traits * n_non_median )
12075 }
121-
122- # get input types from samplingspace
76+
77+ # Input types come from the sampling space; parameters map to the "param" column.
12378 samp <- settings $ ensemble $ samplingspace
12479 input_types <- names(samp )
12580 input_types [input_types == " parameters" <- " param"
126-
81+
12782 if (! " param" %in% input_types ) {
12883 input_types <- c(" param" input_types )
12984 }
130-
131- # build design matrix
132- # key difference from ensemble design:
133- # - ensemble: all columns get random/quasi-random indices
134- # - SA (OAT): param column = sequential index, ALL other columns = 1
135- #
136- # the "1" means: use the FIRST (and only) input file for that type.
137- # this ensures all SA runs use the SAME met, same ic, etc.
13885
86+ # OAT design: the param column carries sequential indices matching the SA run
87+ # order, and every other input column is held constant at 1 (the first input
88+ # file), so each run isolates a single parameter.
13989 design_list <- list ()
140-
90+
14191 for (input_type in input_types ) {
14292 if (input_type == " param"
143- # sequential indices map to SA run order
144- # 1 = median run
145- # 2 = first (pft, trait, quantile) combination
146- # 3 = second (pft, trait, quantile) combination
147- # ...
14893 design_list [[input_type ]] <- seq_len(num_sa_runs )
14994 } else {
150- # all other inputs constant(always use first input file)
15195 design_list [[input_type ]] <- rep(1L , num_sa_runs )
15296 }
15397 }
154-
98+
15599 design_matrix <- data.frame (design_list )
156100
157- return (list (X = design_matrix ))
101+ return (list (X = design_matrix , samples = samples ))
158102}
0 commit comments