From 4b30f2422d11bba35b16193aafde2551d03f5b01 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:58 -0400 Subject: [PATCH 01/30] add period mean, contrast, bind and transform operators --- R/observation_operator.R | 354 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 329 insertions(+), 25 deletions(-) diff --git a/R/observation_operator.R b/R/observation_operator.R index 8303290..92bfbc3 100644 --- a/R/observation_operator.R +++ b/R/observation_operator.R @@ -1,48 +1,72 @@ -# the observation operator: turn a model ensemble's output into the prediction -# matrix G the calibration compares to y. reads each ensemble member's run, -# samples the calibrated variable at each observation date carried in meta, -# converts it from the model's unit to the observation's unit, and lays the -# members out in run id order so G rows align with the parameter ensemble. +# the observation operators: everything that maps between model output, raw +# observation slots, and the fitted target. +# +# harvest_output_to_G reads a model ensemble into the prediction matrix G aligned +# to the raw slots. period_mean_contrast and contrast_target contract raw slots +# into fitted quantities (a period-mean level plus treatment contrasts; a +# per-date treatment contrast), and bind_obs stacks targets. each contraction is +# linear and records itself as a matrix in `transform`, applied identically to +# observations and to G, so the fitted quantity is the same operation on both +# sides by construction. ##' @title Harvest a model ensemble into the prediction matrix G ##' @name harvest_output_to_G ##' @author Akash BV ##' -##' @description For each ENS-- run under `out_root`, reads -##' `variable`, samples it at the observation date of each of that treatment's -##' slots (the midpoint of the slot's date window in meta), converts to the +##' @description For each ENS-- run under `out_root`, reads the +##' model outputs the treatment's slots need over that treatment's run window, +##' samples each slot at the midpoint of its date window, converts to the ##' observation unit, and assembles the J x P matrix aligned to the observation ##' slots. Assumes every expected run has finished; a missing run output fails ##' loud in read.output rather than being silently dropped. ##' ##' @param out_root the model output directory holding the ENS-* run dirs. -##' @param meta observation meta (slot, treatment_id, min_date, max_date). -##' @param variable the model output variable to read. -##' @param start_year,end_year the run year range. -##' @param from_unit the variable's model unit (udunits string). -##' @param to_unit the observation unit to convert to (udunits string). +##' @param meta observation meta (slot, treatment_id, variable, min_date, max_date). +##' @param var_map named list keyed by observation `variable`, each +##' `list(model_var, from, to)`: the model output to read, its unit, and the +##' observation unit. +##' @param run_window integer matrix (2 x n_treatments) of first and last run year, +##' columns named by treatment. Each treatment is read over its own window; a +##' joint run spans different periods per site. ##' @return matrix (members x slots) named by observation slot, member-ordered. ##' @export -harvest_output_to_G <- function(out_root, meta, variable, start_year, end_year, - from_unit, to_unit) { +harvest_output_to_G <- function(out_root, meta, var_map, run_window) { run_dirs <- list.files(out_root, pattern = "^ENS-") rows <- lapply(run_dirs, function(rid) { treat <- sub("^ENS-[0-9]+-", "", rid) member <- as.integer(sub("^ENS-0*([0-9]+)-.*", "\\1", rid)) + md <- meta[meta$treatment_id == treat, ] + if (nrow(md) == 0L) return(NULL) # a run with no matching slots + if (!treat %in% colnames(run_window)) { + PEcAn.logger::logger.severe("no run window for treatment ", treat) + } + win_years <- run_window[, treat] + model_vars <- unique(vapply(md$variable, function(v) var_map[[v]]$model_var, character(1))) o <- PEcAn.utils::read.output( runid = rid, outdir = file.path(out_root, rid), - start.year = start_year, end.year = end_year, - variables = variable, dataframe = TRUE, verbose = FALSE + start.year = win_years[1], end.year = win_years[2], + variables = model_vars, dataframe = TRUE, verbose = FALSE ) - md <- meta[meta$treatment_id == treat, ] - obs_dates <- as.Date(md$min_date) + - (as.Date(md$max_date) - as.Date(md$min_date)) / 2 dates <- as.Date(o$posix) - native <- vapply(obs_dates, - function(d) o[[variable]][which.min(abs(dates - d))], - numeric(1)) - tibble::tibble(member = member, slot = md$slot, - value = PEcAn.utils::ud_convert(native, from_unit, to_unit)) + vals <- vapply(seq_len(nrow(md)), function(i) { + vm <- var_map[[md$variable[i]]] + col <- o[[vm$model_var]] + a <- as.Date(md$min_date[i]); b <- as.Date(md$max_date[i]) + # an observation outside the run window must fail here. substituting the + # nearest available date returns a finite number from the wrong year and makes + # a dry run look clean when the slot was never simulated. + if (a < min(dates) || b > max(dates)) { + PEcAn.logger::logger.severe( + "observation slot ", md$slot[i], " spans ", as.character(a), " to ", + as.character(b), ", outside the ", treat, " run window ", + as.character(min(dates)), " to ", as.character(max(dates)), + ". Extend the run or drop the observation; do not substitute a neighbor." + ) + } + native <- col[which.min(abs(dates - (a + (b - a) / 2)))] + PEcAn.utils::ud_convert(native, vm$from, vm$to) + }, numeric(1)) + tibble::tibble(member = member, slot = md$slot, value = vals) }) long <- do.call(rbind, rows) wide <- tidyr::pivot_wider(long, names_from = "slot", values_from = "value") @@ -51,3 +75,283 @@ harvest_output_to_G <- function(out_root, meta, variable, start_year, end_year, rownames(G) <- wide$member G } + +##' shrink correlations toward zero by the least amount that restores positive +##' definiteness, leaving the diagonal exactly as estimated. a mean over K years +##' gives an empirical covariance of rank at most K - 1, so more quantities than +##' years is singular by construction and the EnKF's Cholesky of cov(G) + Sigma +##' has no reason to succeed. the marginal variances are well estimated from K +##' years; the correlations are not, so only they are damped. +##' @keywords internal +.shrink_to_pd <- function(S, label = "", tol = 1e-8) { + if (min(eigen(S, symmetric = TRUE, only.values = TRUE)$values) > tol) return(S) + D <- diag(diag(S), nrow = nrow(S)) + for (lambda in seq(0.05, 1, by = 0.05)) { + Sl <- (1 - lambda) * S + lambda * D + if (min(eigen(Sl, symmetric = TRUE, only.values = TRUE)$values) > tol) { + PEcAn.logger::logger.info( + "period mean covariance for ", label, " was singular (rank ", + qr(S)$rank, " of ", nrow(S), "); correlations shrunk by ", lambda, + " toward the diagonal, marginal variances kept" + ) + dimnames(Sl) <- dimnames(S) + return(Sl) + } + } + PEcAn.logger::logger.severe( + "period mean covariance for ", label, + " is not positive definite even with correlations fully removed" + ) +} + +##' @title Period mean level and treatment contrasts +##' @name period_mean_contrast +##' @author Akash BV +##' +##' @description Collapses a per treatment per year series into one level slot for +##' the control and one contrast slot per remaining treatment, all on the mean +##' over the period. A bijection of the per treatment means: the control level +##' carries the net rate, the contrasts carry the treatment effects, and one +##' level slot keeps a calibrated initial state out of the contrasts. +##' +##' The covariance is the empirical covariance of the annual series divided by +##' the number of years (the standard error of the mean), which carries the +##' shared-control structure of the contrasts. +##' +##' @param obs a build_obs target list(y, Sigma, meta). +##' @param variable the per treatment per year variable to collapse. +##' @param control treatment id every contrast is taken against. +##' @param years optional integer vector restricting the period; defaults to every +##' year present. +##' @param new_variable name for the resulting variable. +##' @return an obs list carrying one level slot and one contrast slot per +##' treatment, with the contraction recorded in `transform`. +##' @export +period_mean_contrast <- function(obs, variable, control, years = NULL, + new_variable = paste0(variable, "_periodmean")) { + meta <- obs$meta + sub <- meta[meta$variable == variable, , drop = FALSE] + if (nrow(sub) == 0L) { + PEcAn.logger::logger.severe("no slots for variable '", variable, "'") + } + if (!is.null(years)) { + sub <- sub[sub$obs_year %in% as.integer(years), , drop = FALSE] + } + if (!control %in% sub$treatment_id) { + PEcAn.logger::logger.severe( + "control treatment '", control, "' is not in variable '", variable, + "'. Present: ", paste(unique(sub$treatment_id), collapse = ", ") + ) + } + + yrs <- sort(unique(sub$obs_year)) + trts <- unique(sub$treatment_id) + others <- setdiff(trts, control) + + # a mean over an unbalanced panel is not the same quantity across treatments, + # and the contrasts would not be paired by year. refuse rather than average + # whatever is present. + cells <- table(sub$treatment_id, sub$obs_year) + if (any(cells != 1L)) { + bad <- which(cells != 1L, arr.ind = TRUE) + PEcAn.logger::logger.severe( + "period mean needs exactly one slot per treatment per year; ", + nrow(bad), " cell(s) violate that, e.g. treatment ", + rownames(cells)[bad[1, "row"]], " year ", colnames(cells)[bad[1, "col"]] + ) + } + MIN_YEARS <- 3L + if (length(yrs) < MIN_YEARS) { + PEcAn.logger::logger.severe( + "period mean needs at least ", MIN_YEARS, " years to estimate its own ", + "standard error; got ", length(yrs) + ) + } + + slot_of <- function(t, y) sub$slot[sub$treatment_id == t & sub$obs_year == y] + series <- vapply(trts, function(t) { + vapply(yrs, function(y) obs$y[[slot_of(t, y)]], numeric(1)) + }, numeric(length(yrs))) + dimnames(series) <- list(as.character(yrs), trts) + + # columns of the fitted quantity: control level, then each contrast paired by year + X <- cbind(series[, control, drop = FALSE], + series[, others, drop = FALSE] - series[, control]) + level_slot <- paste0(new_variable, "__", control, "__level") + contrast_slots <- paste0(new_variable, "__", others, "__vs_", control) + colnames(X) <- c(level_slot, contrast_slots) + + n_yr <- length(yrs) + y_new <- colMeans(X) + Sigma <- .shrink_to_pd(stats::cov(X) / n_yr, label = new_variable) + + first <- sub[match(c(control, others), sub$treatment_id), , drop = FALSE] + new_meta <- tibble::tibble( + slot = colnames(X), + variable = new_variable, + sitename = first$sitename, + treatment_id = c(control, paste0(others, "_vs_", control)), + obs_year = NA_integer_, + min_date = min(sub$min_date), + max_date = max(sub$max_date), + min_depth = first$min_depth, + max_depth = first$max_depth, + units = first$units, + observation_level = "period_mean", + n_rep = n_yr, + value = unname(y_new), + var_obs = unname(diag(Sigma)) + ) + + # no model variable is "the period mean minus the control", so the operation is + # recorded as a matrix and applied to G. averaging and differencing are both + # linear; one matrix does both. + tmat <- matrix(0, ncol(X), nrow(meta), dimnames = list(colnames(X), meta$slot)) + for (y in yrs) tmat[level_slot, slot_of(control, y)] <- 1 / n_yr + for (k in seq_along(others)) { + for (y in yrs) { + tmat[contrast_slots[k], slot_of(others[k], y)] <- 1 / n_yr + tmat[contrast_slots[k], slot_of(control, y)] <- + tmat[contrast_slots[k], slot_of(control, y)] - 1 / n_yr + } + } + + list(y = stats::setNames(unname(y_new), colnames(X)), + Sigma = Sigma, meta = new_meta, transform = tmat) +} + +##' @title Contract per treatment observations into a treatment contrast +##' @name contrast_target +##' @author Akash BV +##' +##' @description One slot per date, the difference between a treatment and its +##' control on the dates both were measured. Chamber flux data supports +##' treatment comparisons rather than absolute magnitudes, which is what this +##' contraction fits. +##' +##' @param obs a build_obs target list(y, Sigma, meta). +##' @param variable the per treatment variable to contract. +##' @param treatment,control treatment ids to compare. +##' @param new_variable name for the resulting variable. +##' @return an obs list carrying one slot per shared date, with the contraction +##' recorded in `transform`. +##' @export +contrast_target <- function(obs, variable, treatment, control, + new_variable = paste0(variable, "_contrast")) { + meta <- obs$meta + sub <- meta[meta$variable == variable, , drop = FALSE] + if (nrow(sub) == 0L) { + PEcAn.logger::logger.severe("no slots for variable '", variable, "'") + } + present <- unique(sub$treatment_id) + absent <- setdiff(c(treatment, control), present) + if (length(absent) > 0L) { + PEcAn.logger::logger.severe(paste0( + "contrast needs both treatments; missing: ", paste(absent, collapse = ", "), + ". Present: ", paste(present, collapse = ", ") + )) + } + + a <- sub[sub$treatment_id == treatment, , drop = FALSE] + b <- sub[sub$treatment_id == control, , drop = FALSE] + dates <- intersect(a$min_date, b$min_date) + unpaired <- setdiff(union(a$min_date, b$min_date), dates) + if (length(unpaired) > 0L) { + PEcAn.logger::logger.severe(paste0( + length(unpaired), " date(s) do not carry both treatments and cannot form a ", + "contrast: ", paste(utils::head(unpaired, 5), collapse = ", ") + )) + } + + var_obs <- diag(obs$Sigma) + + rows <- lapply(dates, function(d) { + ra <- a[a$min_date == d, , drop = FALSE] + rb <- b[b$min_date == d, , drop = FALSE] + va <- obs$y[[ra$slot]] + vb <- obs$y[[rb$slot]] + r <- ra + r$slot <- paste0(new_variable, "__", d) + r$variable <- new_variable + r$treatment_id <- paste0(treatment, "_vs_", control) + r$value <- va - vb + # the two cells are independent measurements, so their variances add + r$var_obs <- var_obs[[ra$slot]] + var_obs[[rb$slot]] + r + }) + new_meta <- dplyr::bind_rows(rows) + Sigma <- diag(new_meta$var_obs, nrow = nrow(new_meta)) + dimnames(Sigma) <- list(new_meta$slot, new_meta$slot) + + tmat <- matrix(0, nrow(new_meta), nrow(meta), + dimnames = list(new_meta$slot, meta$slot)) + for (d in dates) { + row <- paste0(new_variable, "__", d) + tmat[row, a$slot[a$min_date == d]] <- 1 + tmat[row, b$slot[b$min_date == d]] <- -1 + } + list(y = stats::setNames(new_meta$value, new_meta$slot), + Sigma = Sigma, meta = new_meta, transform = tmat) +} + +##' @title Combine observation targets into one +##' @name bind_obs +##' @author Akash BV +##' +##' @description Stacks several obs lists into a single target with a +##' block-diagonal covariance. Targets from different variables are assumed +##' independent of each other, which is what block diagonal encodes; +##' correlations within a target are carried through from the input blocks. +##' +##' @param ... obs lists, each list(y, Sigma, meta, transform). +##' @return a single obs list with the transforms combined. +##' @export +bind_obs <- function(...) { + parts <- list(...) + slots <- unlist(lapply(parts, function(p) p$meta$slot), use.names = FALSE) + if (anyDuplicated(slots) > 0L) { + PEcAn.logger::logger.severe( + "duplicate slot name(s) across targets: ", + paste(unique(slots[duplicated(slots)]), collapse = ", ") + ) + } + meta <- dplyr::bind_rows(lapply(parts, function(p) tibble::as_tibble(p$meta))) + n <- length(slots) + Sigma <- matrix(0, n, n, dimnames = list(slots, slots)) + for (p in parts) { + Sigma[p$meta$slot, p$meta$slot] <- p$Sigma[p$meta$slot, p$meta$slot] + } + # a target without a transform cannot be combined with contracted ones: the + # forward applies one map to one G. build_obs already stacks raw targets. + no_tf <- vapply(parts, function(p) is.null(p$transform), logical(1)) + if (any(no_tf)) { + PEcAn.logger::logger.severe("target(s) ", paste(which(no_tf), collapse = ", "), + " carry no transform") + } + raw_slots <- unique(unlist(lapply(parts, function(p) colnames(p$transform)))) + tmat <- matrix(0, length(slots), length(raw_slots), + dimnames = list(slots, raw_slots)) + for (p in parts) { + tmat[rownames(p$transform), colnames(p$transform)] <- p$transform + } + list(y = stats::setNames(unlist(lapply(parts, function(p) as.numeric(p$y))), slots), + Sigma = Sigma, meta = meta, transform = tmat) +} + +##' @title Apply a target transform to a prediction matrix +##' @name apply_transform +##' @author Akash BV +##' +##' @description Puts model predictions on the raw slots through the same linear +##' map the observations went through, so the fitted quantity is the same +##' operation on both sides. +##' +##' @param G prediction matrix (members x raw slots). +##' @param transform matrix mapping raw slots (columns) to fitted slots (rows). +##' @return matrix (members x fitted slots), columns named by fitted slot. +##' @export +apply_transform <- function(G, transform) { + out <- G[, colnames(transform), drop = FALSE] %*% t(transform) + colnames(out) <- rownames(transform) + out +} From 98e76103404b6f62bf2214f9aade00f8f409927e Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:58 -0400 Subject: [PATCH 02/30] build the observation target from multiple curated variables --- R/observations.R | 283 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 195 insertions(+), 88 deletions(-) diff --git a/R/observations.R b/R/observations.R index 731c8c8..9383ae7 100644 --- a/R/observations.R +++ b/R/observations.R @@ -1,114 +1,221 @@ -# build a calibration target (y, Sigma, meta) from curated cal/val observations. +# build a calibration target (y, Sigma, meta) from the curated cal/val observations. # -# the curated data holds replicate-level measurements. the Gaussian likelihood -# covariance Sigma is the variance of each (site, treatment, year) cell mean, -# computed from the replicate spread rather than read from a column, with a -# relative floor so a cell whose replicates happen to agree cannot drive its -# variance to zero and make Sigma singular. y is the per-cell replicate mean; -# meta carries the date and depth window per slot so the observation operator can -# align model output to y. - -##' @title Read the curated cal/val observations table -##' @name read_cal_val_observations -##' -##' @description Reads the single observations tsv exported from the cal/val -##' workbook, drops empty trailing sheet rows, and strips thousands-separator -##' commas from the numeric columns so coercion does not silently produce NA. -##' -##' @param cal_val_dir directory holding the observations tsv export. -##' @return tibble of observation rows, value and study_year coerced to numeric. +# the estimator aligns by slot name only, so this layer owns everything about the +# curated table: it reads replicate level records only, every record must carry +# its target's declared unit, two unit bases never average into one cell, the +# likelihood variance comes from replicate spread or the reported standard error, +# and cells are keyed by depth so two depth increments cannot collapse into one +# mean. +# +# dates come from min_date/max_date, never from study_year, because study_year can +# be a study offset at one site and a calendar year at another. + +##' read the observations table from a cal-val-data checkout. ##' @keywords internal read_cal_val_observations <- function(cal_val_dir) { - f <- list.files(cal_val_dir, pattern = "observations\\.tsv$", full.names = TRUE) - if (length(f) != 1L) { - PEcAn.logger::logger.severe( - "expected exactly one observations tsv in ", cal_val_dir, ", found ", length(f) - ) + f <- file.path(cal_val_dir, "data", "observations.csv") + if (!file.exists(f)) { + PEcAn.logger::logger.severe("curated observations not found at ", f) } - strip_commas <- function(x) as.numeric(gsub(",", "", as.character(x))) - readr::read_tsv(f, show_col_types = FALSE) |> - dplyr::filter(!is.na(variable), variable != "") |> - dplyr::mutate(value = strip_commas(value), study_year = strip_commas(study_year)) + readr::read_csv(f, show_col_types = FALSE, progress = FALSE) |> + dplyr::filter(!is.na(variable)) } -##' @title Build the calibration target from cal/val observations +##' variance of each cell mean: replicate spread where the cell has replicates, +##' otherwise the reported standard error. no relative floor -- it would override +##' a reported error by orders of magnitude. +##' @keywords internal +cell_variance <- function(cells) { + v <- dplyr::case_when( + cells$n_rep > 1L & is.finite(cells$cell_sd) & cells$cell_sd > 0 ~ cells$cell_sd^2 / cells$n_rep, + is.finite(cells$reported_se) & cells$reported_se > 0 ~ cells$reported_se^2, + TRUE ~ NA_real_ + ) + + # small-n sd^2/n lands near zero often enough that one cell can take most of a + # target's inverse-variance weight; pooled_cv instead estimates one relative + # measurement error across the target's cells and scales it by each cell's mean. + pooled <- cells$variance_model == "pooled_cv" + for (v_name in unique(cells$variable[pooled])) { + idx <- pooled & cells$variable == v_name + rel <- ifelse(is.finite(cells$cell_sd[idx]) & cells$cell_sd[idx] > 0, + cells$cell_sd[idx], cells$reported_se[idx]) / abs(cells$cell_mean[idx]) + cv <- stats::median(rel, na.rm = TRUE) + if (!is.finite(cv) || cv <= 0) { + PEcAn.logger::logger.severe( + "target ", v_name, " declares pooled_cv but no cell has usable replicate ", + "spread to pool from" + ) + } + v[idx] <- (cv * abs(cells$cell_mean[idx]))^2 / pmax(cells$n_rep[idx], 1L) + } + v +} + +##' @title Build the calibration target from the curated cal/val observations ##' @name build_obs ##' @author Akash BV ##' ##' @description Assembles the observation vector `y`, its diagonal likelihood -##' covariance `Sigma`, and the per-slot `meta` for one variable at one or more -##' sites. Each (site, treatment, year) cell becomes one slot: `y` is the -##' replicate mean, `Sigma` the variance of that mean floored relative to its -##' magnitude. Nothing here is variable- or site specific; the caller names them. +##' covariance `Sigma`, and the per-slot `meta`, stacking every target into one +##' vector. A cell is one (variable, site, treatment, year, depth) group; the slot +##' name carries all of them so the estimator, which aligns by name only, never has +##' to know what a site or a variable is. ##' -##' @param cal_val_dir directory of the cal/val tsv export. -##' @param target_var the cal/val variable to calibrate to. -##' @param sites character vector of sitenames to include. -##' @param rel_var_floor relative floor on each cell variance: the variance is at -##' least `(rel_var_floor * value)^2`, so an agreeing cell cannot make Sigma -##' singular. Set 0 to disable. -##' @return list(y, Sigma, meta): `y` named numeric (length P), `Sigma` a P x P -##' diagonal variance matrix named to match `y`, `meta` one row per slot. +##' @param cal_val_dir the cal-val-data checkout root. +##' @param targets list of target specs: `variable`, `sites`, `units` (the unit +##' every source record must carry), and optionally `source_variables` (raw names +##' feeding it, default `variable`), `years` (inclusive c(first, last) filter), +##' `variance` ("replicate", the default, or "pooled_cv"), and `cell_period` +##' ("year", the default, or "date" for an episodic sub-annual flux). +##' @return list(y, Sigma, meta). ##' @export -build_obs <- function(cal_val_dir, target_var, sites, rel_var_floor = 0.05) { - raw <- read_cal_val_observations(cal_val_dir) |> - dplyr::filter(variable == target_var, sitename %in% sites, - observation_level == "replicate") - if (nrow(raw) == 0L) { +build_obs <- function(cal_val_dir, targets) { + raw_all <- read_cal_val_observations(cal_val_dir) + + cells <- dplyr::bind_rows(lapply(targets, function(tg) { + rows <- target_rows(raw_all, tg) + rows |> + dplyr::summarize( + cell_mean = mean(value, na.rm = TRUE), + cell_sd = stats::sd(value, na.rm = TRUE), + n_rep = sum(!is.na(value)), + reported_se = mean(reported_se, na.rm = TRUE), + min_date = as.character(min(obs_date_start)), + max_date = as.character(max(obs_date_end)), + variable = dplyr::first(target_variable), + units = dplyr::first(target_units), + variance_model = dplyr::first(variance_model), + observation_level = dplyr::first(observation_level), + .by = c(cell, sitename, treatment_id, obs_year, cell_period, min_depth, max_depth) + ) + })) + + cells$var_obs <- cell_variance(cells) + bad <- !is.finite(cells$var_obs) | cells$var_obs <= 0 + if (any(bad)) { PEcAn.logger::logger.severe( - "no replicate-level ", target_var, " rows for site(s) ", - paste(sites, collapse = ", ") + sum(bad), " observation cell(s) have no usable variance (no replicate ", + "spread, no reported standard error): ", + paste(utils::head(cells$cell[bad]), collapse = ", "), + ". Give them a reported error in the curated source or drop them from the target." ) } + cells <- dplyr::arrange(cells, variable, sitename, treatment_id, obs_year) - cells <- raw |> - dplyr::summarize( - cell_mean = mean(value, na.rm = TRUE), - cell_sd = stats::sd(value, na.rm = TRUE), - n_rep = sum(!is.na(value)), - min_date = dplyr::first(min_date), - max_date = dplyr::first(max_date), - min_depth = dplyr::first(min_depth), - max_depth = dplyr::first(max_depth), - .by = c(sitename, treatment_id, study_year) - ) |> - dplyr::arrange(sitename, treatment_id, study_year) |> - dplyr::mutate( - var_mean = dplyr::if_else(n_rep > 1L, cell_sd^2 / n_rep, NA_real_), - var_obs = pmax(var_mean, (rel_var_floor * abs(cell_mean))^2, na.rm = TRUE) - ) - - if (any(!is.finite(cells$var_obs))) { - bad <- paste0(cells$treatment_id, "_y", cells$study_year)[!is.finite(cells$var_obs)] + if (any(duplicated(cells$cell))) { PEcAn.logger::logger.severe( - "non-finite observation variance for: ", paste(utils::head(bad), collapse = ", ") + "duplicate observation slot(s): ", + paste(utils::head(cells$cell[duplicated(cells$cell)]), collapse = ", ") ) } - slot <- paste(cells$sitename, cells$treatment_id, - paste0("y", cells$study_year), sep = "__") - y <- stats::setNames(cells$cell_mean, slot) + y <- stats::setNames(cells$cell_mean, cells$cell) Sigma <- diag(cells$var_obs, nrow = length(y)) - dimnames(Sigma) <- list(slot, slot) - - meta <- tibble::tibble( - slot = slot, - sitename = cells$sitename, - treatment_id = cells$treatment_id, - study_year = cells$study_year, - min_date = cells$min_date, - max_date = cells$max_date, - min_depth = cells$min_depth, - max_depth = cells$max_depth, - n_rep = cells$n_rep, - value = cells$cell_mean, - var_obs = cells$var_obs - ) + dimnames(Sigma) <- list(cells$cell, cells$cell) + + meta <- dplyr::select(cells, slot = cell, variable, sitename, treatment_id, + obs_year, min_date, max_date, min_depth, max_depth, units, + observation_level, n_rep, value = cell_mean, var_obs) PEcAn.logger::logger.info( - length(y), " observation slots for ", target_var, " across ", - dplyr::n_distinct(cells$treatment_id), " treatments, ", - dplyr::n_distinct(cells$study_year), " years" + length(y), " observation slots across ", dplyr::n_distinct(meta$variable), + " variable(s) and ", dplyr::n_distinct(meta$sitename), " site(s)" ) list(y = y, Sigma = Sigma, meta = meta) } + +##' filter, unit-check and key one target's rows (see build_obs for the spec). +##' @keywords internal +target_rows <- function(raw_all, tg) { + stopifnot(!is.null(tg$variable), !is.null(tg$sites), !is.null(tg$units)) + + srcs <- if (is.null(tg$source_variables)) tg$variable else tg$source_variables + rows <- dplyr::filter(raw_all, variable %in% srcs, + sitename %in% as.character(tg$sites), + observation_level == "replicate") + if (nrow(rows) == 0L) { + PEcAn.logger::logger.severe("no rows for target ", tg$variable, + " at site(s) ", paste(tg$sites, collapse = ", ")) + } + + # dates come from min_date/max_date; study_year is not a usable key here. + rows <- rows |> + dplyr::mutate( + obs_date_start = as.Date(min_date), + obs_date_end = as.Date(max_date) + ) + undated <- is.na(rows$obs_date_start) | is.na(rows$obs_date_end) + if (any(undated)) { + PEcAn.logger::logger.severe( + sum(undated), " row(s) of target ", tg$variable, + " have no min_date/max_date, so the operator cannot place them in a run window. ", + "Exclude them from the target or date them in the curated source." + ) + } + rows$obs_year <- as.integer(format(rows$obs_date_start, "%Y")) + if (!is.null(tg$years)) { + rows <- dplyr::filter(rows, obs_year >= tg$years[1], obs_year <= tg$years[2]) + if (nrow(rows) == 0L) { + PEcAn.logger::logger.severe("target ", tg$variable, ": year filter removed every row") + } + } + + rows$target_variable <- tg$variable + rows$target_units <- tg$units + rows$variance_model <- if (is.null(tg$variance)) "replicate" else tg$variance + + # an episodic sub-annual flux keyed by year averages its peaks away before the + # estimator sees them, so a target measured on discrete dates declares + # cell_period = "date" and keeps one cell per measurement date. + period <- if (is.null(tg$cell_period)) "year" else tg$cell_period + rows$cell_period <- switch(period, + year = paste0("y", rows$obs_year), + date = as.character(rows$obs_date_start), + PEcAn.logger::logger.severe("unknown cell_period '", period, + "' for target ", tg$variable, "; use 'year' or 'date'") + ) + rows$cell <- paste(tg$variable, rows$sitename, rows$treatment_id, rows$cell_period, + paste0("d", rows$min_depth, "_", rows$max_depth), sep = "__") + + off <- unique(rows$reported_units[rows$reported_units != tg$units]) + if (length(off) > 0L) { + PEcAn.logger::logger.severe( + "target ", tg$variable, " declares '", tg$units, "' but records carry [", + paste(off, collapse = ", "), "]; convert or correct the curated source" + ) + } + + rows$reported_se <- suppressWarnings(as.numeric(rows$stat)) + rows +} + +##' @title Split an observation target into fitted and validation parts +##' @name subset_obs +##' @author Akash BV +##' +##' @description Keeps slots in the run without letting them into the likelihood; +##' a held-out variable is still predicted and still scored, which is what a +##' validation target is. +##' +##' @param obs a build_obs target list(y, Sigma, meta). +##' @param keep logical vector over rows of `obs$meta`, TRUE to retain. +##' @return an obs list of the same shape carrying only the kept slots. +##' @export +subset_obs <- function(obs, keep) { + if (length(keep) != nrow(obs$meta)) { + PEcAn.logger::logger.severe( + "keep must be one logical per observation slot: got ", length(keep), + " for ", nrow(obs$meta), " slots" + ) + } + if (!any(keep)) { + PEcAn.logger::logger.severe("subset_obs would leave no slots") + } + slots <- obs$meta$slot[keep] + list( + y = obs$y[slots], + Sigma = obs$Sigma[slots, slots, drop = FALSE], + meta = obs$meta[keep, , drop = FALSE] + ) +} From b5bb9f866cf8dddcbcbcb030b0c68afd115c3e5d Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:58 -0400 Subject: [PATCH 03/30] support joint runs with per treatment windows and a calibrated state --- R/forward_sipnet.R | 183 ++++++++++++++++++++++++++++++--------------- 1 file changed, 121 insertions(+), 62 deletions(-) diff --git a/R/forward_sipnet.R b/R/forward_sipnet.R index 507d90c..831d514 100644 --- a/R/forward_sipnet.R +++ b/R/forward_sipnet.R @@ -2,11 +2,10 @@ # fwd(U, iteration) -> G the calibration calls. only this file knows sipnet and # the pecan run machinery; the estimator in method_eki.R knows neither. # -# one iteration is one ensemble where only the calibrated parameters change; met, -# events, and the pools we do not calibrate are built once and held fixed. a call -# writes the proposal into the sample object and, when an initial state is -# calibrated, the ic files, runs the pecan config and model steps, and harvests -# the output into G. +# one iteration is one ensemble where only the calibrated parameters change: met, +# events, and uncalibrated pools are pinned to one member, otherwise the prediction +# spread measures the input draw and cov(U, G) in the kalman gain is sampling +# noise. input uncertainty belongs in a separate forward pass. # # launching is left to pecan and the prepared host block. runModule_start_model_runs # submits through the settings host (qsub, sge_array_launcher.sh, Njobmax, qstat) @@ -26,31 +25,41 @@ ##' @param settings a prepared PEcAn multisite settings object (the forward run). ##' @param obs the build_obs target list(y, Sigma, meta); names(y) are the slots. ##' @param n_particles ensemble size J. -##' @param harvest_var the model output variable to compare to the observations. -##' @param from_unit the model unit of harvest_var (udunits string). -##' @param to_unit the observation unit to convert the harvest to. -##' @param soil_pft name of the PFT whose traits U overwrites. +##' @param var_map named list keyed by observation variable, each +##' `list(model_var, from, to)` (see harvest_output_to_G): the crosswalk from +##' each observed variable to its model output and units. +##' @param soil_pfts character vector of soil PFT names that share the calibrated +##' rates; the same proposal column is written into each (see inject_traits). ##' @param state_prefix column prefix marking calibrated initial-state entries in ##' U; when present each is written into the initial condition per particle. ##' @param state_pool the initial condition pool variable the state writes into. ##' @param base_out_dir parent directory for the per iteration model output. +##' @param fixed_traits trait names pinned through the run dir `default.param`, +##' dropped from the baseline sample so the pinned value is not overwritten by +##' the PFT posterior median. +##' @param raw_obs the untransformed target the model output is harvested against; +##' its `transform` (the linear map from raw to fitted slots) is applied to G so +##' model and observations are the same quantity. NULL fits the raw slots. ##' @return function(U, iteration) -> matrix (J, P) aligned to names(obs$y). ##' @export -make_forward_sipnet <- function(settings, obs, n_particles, - harvest_var, from_unit, to_unit, - soil_pft = "soil", +make_forward_sipnet <- function(settings, obs, n_particles, var_map, + soil_pfts, state_prefix = "soilInit.", state_pool = "soil_organic_carbon_content", - base_out_dir = settings$outdir) { + base_out_dir = settings$outdir, + fixed_traits = character(0), + raw_obs = NULL) { + transform <- raw_obs$transform + if (!is.null(raw_obs) && is.null(transform)) { + PEcAn.logger::logger.severe( + "raw_obs carries no transform; the fitted target cannot be reached from the ", + "model output without one" + ) + } + harvest_meta <- if (is.null(raw_obs)) obs$meta else raw_obs$meta obs_order <- names(obs$y) meta <- obs$meta - start_year <- as.integer(format(as.Date(settings[[1]]$run$start.date), "%Y")) - end_year <- as.integer(format(as.Date(settings[[1]]$run$end.date), "%Y")) - - inputs <- settings[[1]]$run$inputs - n_met <- length(inputs$met$path) - n_events <- length(inputs$events$path) - n_ic <- length(inputs$poolinitcond$path) + window <- run_window(settings) # the multisite site ids, iterated in settings order treatments <- vapply(settings, function(x) x$run$site$id, character(1)) @@ -62,10 +71,7 @@ make_forward_sipnet <- function(settings, obs, n_particles, function(i) settings[[i]]$run$inputs$poolinitcond$path[[1]]), treatments) - baseline <- baseline_trait_samples(settings$pfts, n_particles) - # config + model steps run from the prepared run dir: the host's launcher path - # is relative to it. - run_dir <- dirname(settings$outdir) + baseline <- baseline_trait_samples(settings$pfts, n_particles, fixed_traits) function(U, itr) { out_itr <- file.path(base_out_dir, paste0("itr", itr)) @@ -80,8 +86,7 @@ make_forward_sipnet <- function(settings, obs, n_particles, state_cols <- grep(paste0("^", state_prefix), colnames(U), value = TRUE) trait_cols <- setdiff(colnames(U), state_cols) - ensemble.samples <- inject_traits(baseline, soil_pft, U[, trait_cols, drop = FALSE]) - write_samples_rdata(ensemble.samples, file.path(s$outdir, "samples.Rdata")) + ensemble.samples <- inject_traits(baseline, soil_pfts, U[, trait_cols, drop = FALSE]) # calibrated initial state: write one ic per (site, particle) with the pool # set to that particle's proposal and point the run's poolinitcond at them, so @@ -91,41 +96,60 @@ make_forward_sipnet <- function(settings, obs, n_particles, ic_paths <- write_state_ensemble(U, state_cols, state_prefix, state_pool, treatments, file.path(out_itr, "IC_files"), ic_template_paths) - s <- repoint_poolinitcond(s, treatments, ic_paths) + s <- repoint_poolinitcond(s, treatments, ic_paths, n_particles) poolinitcond_idx <- seq_len(n_particles) } else { - poolinitcond_idx <- rep_len(seq_len(n_ic), n_particles) + poolinitcond_idx <- rep(1L, n_particles) } - input_design <- tibble::tibble( - param = seq_len(n_particles), - poolinitcond = poolinitcond_idx, - met = rep_len(seq_len(n_met), n_particles), - events = rep_len(seq_len(n_events), n_particles) + # a design's `param` column indexes into the samples it was drawn with, so the + # two arrive together; supplying both also keeps pecan from generating its own + # design, which pins met and events to the first member. + input_design <- list( + design_matrix = data.frame( + param = seq_len(n_particles), + poolinitcond = poolinitcond_idx, + met = 1L, + events = 1L + ), + samples = list( + ensemble.samples = ensemble.samples, + trait.samples = lapply(ensemble.samples, as.list), + sa.samples = NULL, + runs.samples = list(), + env.samples = list() + ) ) - old_wd <- setwd(run_dir) - on.exit(setwd(old_wd), add = TRUE) s <- PEcAn.workflow::runModule.run.write.configs(s, input_design = input_design) PEcAn.workflow::runModule_start_model_runs(s, stop.on.error = FALSE) - G <- harvest_output_to_G(s$modeloutdir, meta, harvest_var, - start_year, end_year, from_unit, to_unit) + G <- harvest_output_to_G(s$modeloutdir, harvest_meta, var_map, window) + if (!is.null(transform)) G <- apply_transform(G, transform) + missing <- setdiff(obs_order, colnames(G)) + if (length(missing) > 0L) { + PEcAn.logger::logger.severe( + length(missing), " observation slots have no forward output (e.g. ", + paste(utils::head(missing, 5), collapse = ", "), + "); every treatment's ensemble runs must finish before harvest" + ) + } G[, obs_order, drop = FALSE] } } -##' fixed baseline trait samples: every parameter at its prior median (the -##' PEcAn.priors::get.sample p = 0.5 of the post.distns row, so it matches the -##' family the pft carries), replicated over particles, one data.frame per pft. -##' the calibrated columns are overwritten by U; the rest stay fixed so the -##' prediction spread reflects only the estimated parameters. +##' baseline trait samples: every parameter at its prior median, replicated over +##' particles, one data.frame per pft. calibrated columns are overwritten by U. +##' traits named in `fixed` are dropped so the run dir default.param value stands; +##' a trait left in here overwrites it. ##' @keywords internal -baseline_trait_samples <- function(pfts, n_particles) { +baseline_trait_samples <- function(pfts, n_particles, fixed = character(0)) { out <- list() for (pft in pfts) { e <- new.env() load(pft$posterior.files, envir = e) pd <- get(ls(e)[[1]], envir = e) + keep <- setdiff(rownames(pd), fixed) + pd <- pd[keep, , drop = FALSE] med <- vapply(seq_len(nrow(pd)), function(i) { PEcAn.priors::get.sample(pd[i, c("distn", "parama", "paramb")], p = 0.5) }, numeric(1)) @@ -137,27 +161,26 @@ baseline_trait_samples <- function(pfts, n_particles) { out } -##' overwrite the pft's calibrated trait columns with the proposal U. +##' write the proposal U into each named soil pft: the calibrated rates are one +##' shared quantity, not one per pft. fails if a named pft is absent rather than +##' silently calibrating a subset. ##' @keywords internal -inject_traits <- function(baseline, soil_pft, U_traits) { +inject_traits <- function(baseline, soil_pfts, U_traits) { + missing_pfts <- setdiff(soil_pfts, names(baseline)) + if (length(missing_pfts) > 0L) { + PEcAn.logger::logger.severe( + "soil PFT(s) not present in the prepared settings: ", + paste(missing_pfts, collapse = ", "), "; the run carries ", + paste(names(baseline), collapse = ", ") + ) + } es <- baseline - for (nm in colnames(U_traits)) es[[soil_pft]][[nm]] <- U_traits[, nm] + for (pft in soil_pfts) { + for (nm in colnames(U_traits)) es[[pft]][[nm]] <- U_traits[, nm] + } es } -##' write samples.Rdata in the object run.write.configs expects. -##' @keywords internal -write_samples_rdata <- function(ensemble.samples, file) { - trait.samples <- lapply(ensemble.samples, as.list) - pft.names <- names(ensemble.samples) - trait.names <- lapply(ensemble.samples, names) - sa.samples <- NULL - runs.samples <- list() - env.samples <- list() - save(ensemble.samples, trait.samples, sa.samples, runs.samples, - pft.names, trait.names, env.samples, file = file) -} - ##' write the per-(site, particle) initial condition ensemble for a calibrated ##' state. for each site, read its fixed pools once from an existing ic, then ##' write one ic per particle with `state_pool` set to that particle's proposal @@ -182,11 +205,47 @@ write_state_ensemble <- function(U, state_cols, prefix, state_pool, treatments, } ##' point each site's poolinitcond path at the freshly written per particle ics, -##' so particle j uses its own ic at every site and the design indexes 1:J. +##' so particle j uses its own ic at every site and the design indexes 1:J. sites +##' without a calibrated state (no `ic_paths` entry) keep their template ic, +##' recycled to J paths so the shared design column stays in range; overwriting +##' them with an empty list silently drops every one of their run dirs. ##' @keywords internal -repoint_poolinitcond <- function(settings, treatments, ic_paths) { +repoint_poolinitcond <- function(settings, treatments, ic_paths, n_particles) { for (i in seq_along(treatments)) { - settings[[i]]$run$inputs$poolinitcond$path <- as.list(ic_paths[[treatments[i]]]) + t <- treatments[i] + if (!is.null(ic_paths[[t]])) { + settings[[i]]$run$inputs$poolinitcond$path <- as.list(ic_paths[[t]]) + } else { + have <- unlist(settings[[i]]$run$inputs$poolinitcond$path, use.names = FALSE) + if (length(have) == 0L) { + PEcAn.logger::logger.severe( + "block ", t, " has neither a calibrated state column nor a pinned ", + "poolinitcond path; every block needs an initial condition" + ) + } + settings[[i]]$run$inputs$poolinitcond$path <- + as.list(rep_len(have, n_particles)) + } } settings } + +##' @title Run years per treatment from a multisite settings object +##' @name run_window +##' @author Akash BV +##' +##' @description First and last run year of each treatment, for reading model +##' output over that treatment's own window; a joint run spans different +##' periods per site. +##' +##' @param settings a PEcAn multisite settings object. +##' @return integer matrix (2 x n_treatments), columns named by treatment. +##' @export +run_window <- function(settings) { + win <- vapply(settings, function(x) { + c(as.integer(format(as.Date(x$run$start.date), "%Y")), + as.integer(format(as.Date(x$run$end.date), "%Y"))) + }, integer(2)) + colnames(win) <- vapply(settings, function(x) x$run$site$id, character(1)) + win +} From c315e4bb65f55bfb68a49a225e18067d0f676c5d Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:58 -0400 Subject: [PATCH 04/30] share one prior across soil pfts and anchor the state on observations --- R/priors.R | 104 ++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 87 insertions(+), 17 deletions(-) diff --git a/R/priors.R b/R/priors.R index ab7e6b8..33e28e1 100644 --- a/R/priors.R +++ b/R/priors.R @@ -1,18 +1,35 @@ # priors for the calibrated parameters, plus the initial ensemble draw. # # priors are not hand tuned. they come from a pft's meta-analysis posterior -# (post.distns, read straight off disk), from an -# explicit specification (a biologically plausible, zero bounded distribution -# written down), or, for a calibrated initial state, anchored to an observation -# with its measurement uncertainty as the spread. each returns the same dist_list -# the transport maps (transport.R) and the sampler consume. +# (post.distns, read straight off disk), from an explicit specification, or, for +# a calibrated initial state, anchored to an observation with its measurement +# uncertainty as the spread. each returns the same dist_list the transport maps +# (transport.R) and the sampler consume. records carry the PEcAn (distn, parama, +# paramb) triple, so PEcAn.priors does the drawing here and the estimating +# upstream (fit.dist for samples, prior.fn for elicited quantiles); this file +# only constructs dist_lists from sources that already exist. ##' distn family -> support the transport map uses. positive-support families get -##' a log map, beta a logit, normal an identity, and a uniform its own bounds -##' (handled by prior_from_specs, not here). +##' a log map, beta and unif a logit onto their own bounds, normal an identity. +##' +##' there is deliberately no default branch: a family falling through to +##' c(-Inf, Inf) estimates a bounded parameter on the whole real line, and the +##' posterior can leave the physical range without anything objecting. a new +##' family must state its support rather than inherit an unbounded one. +##' @param distn family name. +##' @param parama,paramb distribution parameters; uniform takes its support from them. ##' @keywords internal -.distn_support <- function(distn) { - switch(as.character(distn), +.distn_support <- function(distn, parama = NULL, paramb = NULL) { + distn <- as.character(distn) + if (identical(distn, "unif")) { + if (length(parama) != 1L || length(paramb) != 1L || anyNA(c(parama, paramb))) { + PEcAn.logger::logger.severe( + "uniform prior needs parama and paramb to define its support" + ) + } + return(c(parama, paramb)) + } + switch(distn, weibull = , lnorm = , gamma = , @@ -21,7 +38,10 @@ geom = c(0, Inf), beta = c(0, 1), norm = c(-Inf, Inf), - c(-Inf, Inf) + PEcAn.logger::logger.severe( + "no declared support for distribution family '", distn, + "'. Add it to .distn_support rather than letting it default to the real line." + ) ) } @@ -56,14 +76,52 @@ prior_from_postdistns <- function(params, post_distns_path) { ) } stats::setNames(lapply(params, function(p) { + if (!p %in% rownames(pd)) { + PEcAn.logger::logger.severe("no posterior for '", p, "' in ", post_distns_path) + } row <- pd[p, ] list(param_name = p, len = 1L, - constraint = .distn_support(row$distn), + constraint = .distn_support(row$distn, row$parama, row$paramb), distn = as.character(row$distn), parama = row$parama, paramb = row$paramb) }), params) } +##' @title Prior for a rate shared across several PFTs +##' @name prior_from_shared_postdistns +##' @author Akash BV +##' +##' @description A rate written into several PFTs is one calibrated quantity, so +##' there must be one prior for it. This reads the requested traits from each +##' named PFT's posterior and requires them to agree before returning a single +##' record; disagreement is an error rather than a quiet choice of the first. +##' +##' @param params character vector of PEcAn trait names to calibrate. +##' @param posterior_files named character vector of post.distns paths, one per soil PFT. +##' @return a dist_list, one record per trait. +##' @export +prior_from_shared_postdistns <- function(params, posterior_files) { + stopifnot(length(posterior_files) >= 1L, !is.null(names(posterior_files))) + per_pft <- lapply(posterior_files, function(f) prior_from_postdistns(params, f)) + reference <- per_pft[[1]] + for (i in seq_along(per_pft)[-1]) { + for (p in params) { + a <- reference[[p]][c("distn", "parama", "paramb")] + b <- per_pft[[i]][[p]][c("distn", "parama", "paramb")] + if (!isTRUE(all.equal(a, b))) { + PEcAn.logger::logger.severe( + "prior for '", p, "' differs between soil PFTs '", names(posterior_files)[1], + "' (", a$distn, " ", a$parama, ", ", a$paramb, ") and '", + names(posterior_files)[i], "' (", b$distn, " ", b$parama, ", ", b$paramb, + "). One shared calibrated rate needs one prior; reconcile the PFTs or ", + "calibrate them separately." + ) + } + } + } + reference +} + ##' @title Priors from an explicit specification ##' @name prior_from_specs ##' @author Akash BV @@ -82,7 +140,7 @@ prior_from_specs <- function(specs) { distn <- as.character(s$distn) a <- as.numeric(s$parama) b <- as.numeric(s$paramb) - constraint <- if (identical(distn, "unif")) c(a, b) else .distn_support(distn) + constraint <- .distn_support(distn, a, b) list(param_name = p, len = 1L, constraint = constraint, distn = distn, parama = a, paramb = b) }), names(specs)) @@ -99,17 +157,29 @@ prior_from_specs <- function(specs) { ##' into the state. The observation is converted from its reported unit to the ##' model's unit with PEcAn.utils::ud_convert. ##' -##' @param meta observation meta (treatment_id, study_year, value, var_obs). +##' @param meta observation meta (variable, treatment_id, obs_year, value, var_obs). ##' @param prefix column prefix marking the state entries (e.g. "soilInit."). +##' @param variable the observed variable whose slots anchor the state; the joint +##' meta can span several variables and an unscoped anchor would take the wrong one. ##' @param from_unit unit the observation is reported in (udunits string). ##' @param to_unit unit the model state is in (udunits string). -##' @param anchor_year study_year to anchor on; defaults to the earliest. +##' @param anchor_year observation year to anchor on; defaults to the earliest. ##' @return a dist_list keyed , in site order. ##' @export -state_prior_from_obs <- function(meta, prefix, from_unit, to_unit, +state_prior_from_obs <- function(meta, prefix, from_unit, to_unit, variable, anchor_year = NULL) { - if (is.null(anchor_year)) anchor_year <- min(meta$study_year) - base <- meta[meta$study_year == anchor_year, ] + # scope to one variable before anchoring: a joint meta can span several + # variables and sites, and an unscoped selection would anchor the state + # on whichever variable sorts first. + base <- meta[meta$variable == variable, ] + if (nrow(base) == 0L) { + PEcAn.logger::logger.severe("no ", variable, " slots to anchor an initial state on") + } + if (is.null(anchor_year)) anchor_year <- min(base$obs_year) + base <- base[base$obs_year == anchor_year, ] + if (nrow(base) == 0L) { + PEcAn.logger::logger.severe("no ", variable, " slots in anchor year ", anchor_year) + } base <- base[order(base$treatment_id), ] center <- PEcAn.utils::ud_convert(base$value, from_unit, to_unit) spread <- PEcAn.utils::ud_convert(sqrt(base$var_obs), from_unit, to_unit) From 32e3afadf5d02d68ddb0eab87381c5e783edc506 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:58 -0400 Subject: [PATCH 05/30] add slot layout, trace, validation and treatment effect figures --- R/plots.R | 291 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 281 insertions(+), 10 deletions(-) diff --git a/R/plots.R b/R/plots.R index 363c861..1d49428 100644 --- a/R/plots.R +++ b/R/plots.R @@ -10,14 +10,14 @@ ensemble_long <- function(G, meta) { tibble::as_tibble(G, rownames = "member") |> dplyr::mutate(member = as.integer(member)) |> tidyr::pivot_longer(-member, names_to = "slot", values_to = "value") |> - dplyr::left_join(dplyr::distinct(meta, slot, treatment_id, study_year), + dplyr::left_join(dplyr::distinct(meta, slot, treatment_id, obs_year), by = "slot") } ##' observation table for plotting (mean and sd band from the cell variance). ##' @keywords internal obs_long <- function(meta) { - dplyr::transmute(meta, treatment_id, study_year, obs = value, + dplyr::transmute(meta, treatment_id, obs_year, obs = value, obs_sd = sqrt(var_obs)) } @@ -37,24 +37,24 @@ obs_long <- function(meta) { plot_ensembles_vs_truth <- function(G, meta, ylab = "value", title = NULL) { el <- ensemble_long(G, meta) band <- el |> - dplyr::summarise(q05 = stats::quantile(value, 0.05), + dplyr::summarize(q05 = stats::quantile(value, 0.05), q95 = stats::quantile(value, 0.95), m = mean(value), - .by = c(treatment_id, study_year)) + .by = c(treatment_id, obs_year)) ol <- obs_long(meta) ggplot2::ggplot() + ggplot2::geom_ribbon(data = band, - ggplot2::aes(study_year, ymin = q05, ymax = q95), + ggplot2::aes(obs_year, ymin = q05, ymax = q95), fill = "grey75", alpha = 0.55) + ggplot2::geom_line(data = el, - ggplot2::aes(study_year, value, group = member), + ggplot2::aes(obs_year, value, group = member), color = "steelblue", alpha = 0.35, linewidth = 0.3) + - ggplot2::geom_line(data = band, ggplot2::aes(study_year, m), linewidth = 0.7) + + ggplot2::geom_line(data = band, ggplot2::aes(obs_year, m), linewidth = 0.7) + ggplot2::geom_pointrange(data = ol, - ggplot2::aes(study_year, obs, ymin = obs - obs_sd, ymax = obs + obs_sd), + ggplot2::aes(obs_year, obs, ymin = obs - obs_sd, ymax = obs + obs_sd), color = "firebrick", linewidth = 0.3, size = 0.25) + ggplot2::facet_wrap(~treatment_id, ncol = 2) + - ggplot2::labs(x = "study year", y = ylab, title = title) + + ggplot2::labs(x = "year", y = ylab, title = title) + ggplot2::theme_minimal(base_size = 11) } @@ -96,9 +96,280 @@ plot_param_densities <- function(U_list, pdf_file, width = 7, height = 4) { print( ggplot2::ggplot(df, ggplot2::aes(value, color = stage)) + ggplot2::geom_density(linewidth = 1) + - ggplot2::labs(title = p, x = "value", y = "density", color = NULL) + + # no page title: the x label already names the parameter and its unit, and a + # title that repeats the axis is noise + ggplot2::labs(x = param_label_unit(p), y = "density", color = NULL) + ggplot2::theme_minimal(base_size = 12) ) } invisible(pdf_file) } + +##' axis and strip labels for calibrated parameters; unmapped names fall back to +##' the trait name with underscores as spaces. +##' @keywords internal +param_label <- function(x) { + labels <- c( + som_respiration_rate = "SOM respiration rate", + fracLitterRespired = "Litter fraction respired", + n_volatilization_rate = "N volatilization rate", + soil_respiration_Q10 = "Soil respiration Q10", + litterBreakdownRate = "Litter breakdown rate" + ) + trait <- sub("^[^.]+\\.", "", x) + out <- labels[trait] + fallback <- is.na(out) + out[fallback] <- gsub("_", " ", trait[fallback]) + unname(out) +} + +##' units as SIPNET reads them, not always as documented: baseSoilResp is read in +##' per year and divided by 365 at setup in sipnet.c, so a calibrated value is a +##' per year rate; wrong on an axis is a factor of 365. +##' @keywords internal +param_unit <- function(x) { + units <- c( + som_respiration_rate = "g C g-1 soil C yr-1", + fracLitterRespired = "unitless", + n_volatilization_rate = "unitless", + soil_respiration_Q10 = "unitless", + litterBreakdownRate = "yr-1" + ) + trait <- sub("^[^.]+\\.", "", x) + out <- units[trait] + out[is.na(out)] <- "" + unname(out) +} + +##' parameter label with its unit appended, for a facet strip or page label. +##' @keywords internal +param_label_unit <- function(x) { + u <- param_unit(x) + lab <- param_label(x) + ifelse(nzchar(u) & u != "unitless", paste0(lab, " (", u, ")"), lab) +} + +##' @title Parameter trace across tempering steps +##' @author Akash BV +##' +##' @description Ensemble mean and 5-95 % spread of each calibrated parameter at the +##' prior and after each tempering step, with its declared support drawn in where +##' finite. Shows whether a parameter settles inside its support or moves onto a +##' bound, and whether it walks there or jumps. +##' +##' @param trace named list of parameter matrices, prior first then one per step. +##' @param prior the dist_list the run used, for the support lines. +##' @return a ggplot object. +##' @export +plot_param_trace <- function(trace, prior = NULL) { + long <- dplyr::bind_rows(lapply(seq_along(trace), function(i) { + U <- trace[[i]] + dplyr::bind_rows(lapply(colnames(U), function(k) { + tibble::tibble(step = i - 1L, param = k, + m = mean(U[, k]), + lo = stats::quantile(U[, k], 0.05), + hi = stats::quantile(U[, k], 0.95)) + })) + })) + long$label <- param_label_unit(long$param) + + bounds <- NULL + if (!is.null(prior)) { + bounds <- dplyr::bind_rows(lapply(names(prior), function(k) { + lim <- prior[[k]]$constraint + tibble::tibble(param = k, label = param_label_unit(k), bound = lim[is.finite(lim)]) + })) + } + + p <- ggplot2::ggplot(long, ggplot2::aes(step, m)) + if (!is.null(bounds) && nrow(bounds) > 0) { + p <- p + ggplot2::geom_hline(data = bounds, + ggplot2::aes(yintercept = bound), color = "grey70", linetype = 2) + } + p + + ggplot2::geom_ribbon(ggplot2::aes(ymin = lo, ymax = hi), fill = "grey75", alpha = 0.55) + + ggplot2::geom_line(linewidth = 0.7) + + ggplot2::geom_point(size = 2) + + ggplot2::facet_wrap(~label, scales = "free_y", + labeller = ggplot2::label_wrap_gen(26)) + + ggplot2::labs(x = "tempering step", y = "value") + + ggplot2::theme_minimal(base_size = 11) +} + +##' @title Held-out validation: predicted against observed +##' @author Akash BV +##' +##' @description Ensemble mean prediction against the observation for slots kept out +##' of the likelihood, with the 5-95 % ensemble range and a 1:1 line. +##' +##' Panels are on free scales: variables held out together can span very different +##' ranges, and a shared scale would flatten the smaller one. +##' +##' @param G prediction matrix covering the validation slots. +##' @param meta observation meta for those slots, carrying `value` and `units`. +##' @return a ggplot object. +##' @export +plot_validation <- function(G, meta) { + G <- G[, meta$slot, drop = FALSE] + d <- tibble::tibble( + variable = meta$variable, + obs = meta$value, + m = colMeans(G), + lo = apply(G, 2, stats::quantile, 0.05), + hi = apply(G, 2, stats::quantile, 0.95) + ) + d$label <- label_variable(d$variable) + # the axis unit comes from the observations the model output was converted onto; + # slots spanning more than one unit cannot share an axis pair + unit <- unique(meta$units) + if (length(unit) != 1L) { + PEcAn.logger::logger.severe( + "validation slots span more than one unit: ", paste(unit, collapse = ", ") + ) + } + # a log axis silently drops every non-positive point, and a treatment contrast is + # signed by construction, so the scale follows the data rather than the variable. + positive <- all(c(d$obs, d$m, d$lo, d$hi) > 0, na.rm = TRUE) + scales <- if (positive) { + list(ggplot2::scale_x_log10(), ggplot2::scale_y_log10()) + } else { + list(ggplot2::geom_hline(yintercept = 0, color = "grey85", linewidth = 0.3), + ggplot2::geom_vline(xintercept = 0, color = "grey85", linewidth = 0.3)) + } + + ggplot2::ggplot(d, ggplot2::aes(obs, m)) + + scales + + ggplot2::geom_abline(slope = 1, intercept = 0, color = "grey70") + + ggplot2::geom_linerange(ggplot2::aes(ymin = lo, ymax = hi), + color = "grey50", linewidth = 0.3) + + ggplot2::geom_point(color = "firebrick", size = 2) + + ggplot2::facet_wrap(~label, scales = "free") + + ggplot2::labs(x = paste0("observed (", unit, ")"), + y = paste0("predicted (", unit, ")")) + + ggplot2::theme_minimal(base_size = 11) +} + +##' readable variable label: the variable name with underscores as spaces. +##' @keywords internal +label_variable <- function(x) { + gsub("_", " ", x) +} + +##' @title Fitted target slots against observations +##' @name plot_target_slots +##' @author Akash BV +##' +##' @description A contracted target has no year axis (obs_year is NA by +##' construction), so its slots go on a categorical axis: member points, +##' ensemble mean, observation +/- sd. +##' +##' @param G prediction matrix covering the slots in `meta`. +##' @param meta observation meta of the fitted slots. +##' @param ylab y axis label. +##' @param title optional factual page label for a multi page pdf. +##' @return a ggplot object. +##' @export +plot_target_slots <- function(G, meta, ylab = "value", title = NULL) { + lab <- sub("_vs_", " - ", meta$treatment_id) + d <- data.frame( + slot = factor(rep(lab, each = nrow(G)), levels = lab), + value = as.vector(G[, meta$slot, drop = FALSE]) + ) + m <- data.frame(slot = factor(lab, levels = lab), + value = colMeans(G[, meta$slot, drop = FALSE])) + o <- data.frame(slot = factor(lab, levels = lab), obs = meta$value, + obs_sd = sqrt(meta$var_obs)) + ggplot2::ggplot() + + ggplot2::geom_jitter(data = d, ggplot2::aes(slot, value), + width = 0.12, height = 0, color = "steelblue", alpha = 0.35, size = 0.8) + + ggplot2::geom_point(data = m, ggplot2::aes(slot, value), shape = 95, size = 8) + + ggplot2::geom_pointrange(data = o, + ggplot2::aes(slot, obs, ymin = obs - obs_sd, ymax = obs + obs_sd), + color = "firebrick", linewidth = 0.4, size = 0.3) + + ggplot2::labs(x = NULL, y = ylab, title = title) + + ggplot2::theme_minimal(base_size = 11) + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 30, hjust = 1)) +} + +##' @title Measured against modeled treatment effect, before and after calibration +##' @name plot_treatment_effect +##' @author Akash BV +##' +##' @description Per year treatment minus control from the prior and posterior +##' forward ensembles against the measured contrast, as the absolute effect and +##' as percent of control. Bands are the 5-95 % ensemble range around the mean; +##' measurement bars are one standard deviation, the two arms' variances added +##' (delta method for the relative form). +##' +##' @param G_prior,G_post prediction matrices on the raw slots (members x slots). +##' @param meta raw observation meta covering the paired slots. +##' @param variable the per treatment per year variable the effect is on. +##' @param treatment,control treatment ids to compare. +##' @return a ggplot object. +##' @export +plot_treatment_effect <- function(G_prior, G_post, meta, variable, treatment, control) { + sub <- meta[meta$variable == variable & + meta$treatment_id %in% c(treatment, control), , drop = FALSE] + yrs <- sort(unique(sub$obs_year)) + slot_of <- function(t, y) sub$slot[sub$treatment_id == t & sub$obs_year == y] + ok <- vapply(yrs, function(y) { + length(slot_of(treatment, y)) == 1L && length(slot_of(control, y)) == 1L + }, logical(1)) + if (length(yrs) == 0L || !all(ok)) { + PEcAn.logger::logger.severe( + "treatment effect needs one ", treatment, " and one ", control, + " slot per year of ", variable, "; unpaired: ", + paste(yrs[!ok], collapse = ", ") + ) + } + + summarize_effect <- function(G, stage) { + dplyr::bind_rows(lapply(yrs, function(y) { + gt <- G[, slot_of(treatment, y)] + gc <- G[, slot_of(control, y)] + eff <- list(absolute = gt - gc, relative = 100 * (gt - gc) / gc) + tibble::tibble( + year = y, stage = stage, form = names(eff), + m = vapply(eff, mean, numeric(1)), + lo = vapply(eff, stats::quantile, numeric(1), 0.05), + hi = vapply(eff, stats::quantile, numeric(1), 0.95) + ) + })) + } + bands <- dplyr::bind_rows(summarize_effect(G_prior, "before calibration"), + summarize_effect(G_post, "after calibration")) + bands$stage <- factor(bands$stage, levels = c("before calibration", "after calibration")) + + mt <- sub[match(paste(treatment, yrs), paste(sub$treatment_id, sub$obs_year)), ] + mc <- sub[match(paste(control, yrs), paste(sub$treatment_id, sub$obs_year)), ] + vt <- mt$value + vc <- mc$value + obs <- tibble::tibble( + year = rep(yrs, 2), + form = rep(c("absolute", "relative"), each = length(yrs)), + obs = c(vt - vc, 100 * (vt - vc) / vc), + obs_sd = c(sqrt(mt$var_obs + mc$var_obs), + 100 * sqrt(mt$var_obs / vc^2 + vt^2 * mc$var_obs / vc^4)) + ) + + form_labels <- c(absolute = paste0("absolute effect (", unique(sub$units), ")"), + relative = "relative effect (% of control)") + ggplot2::ggplot() + + ggplot2::geom_hline(yintercept = 0, color = "grey70", linewidth = 0.3) + + ggplot2::geom_ribbon(data = bands, + ggplot2::aes(year, ymin = lo, ymax = hi, fill = stage), alpha = 0.4) + + ggplot2::geom_line(data = bands, ggplot2::aes(year, m, color = stage), + linewidth = 0.7) + + ggplot2::geom_pointrange(data = obs, + ggplot2::aes(year, obs, ymin = obs - obs_sd, ymax = obs + obs_sd), + color = "firebrick", linewidth = 0.4, size = 0.3) + + ggplot2::scale_fill_manual(values = c("before calibration" = "grey75", + "after calibration" = "steelblue"), name = NULL) + + ggplot2::scale_color_manual(values = c("before calibration" = "grey40", + "after calibration" = "steelblue4"), name = NULL) + + ggplot2::facet_wrap(~form, scales = "free_y", + labeller = ggplot2::labeller(form = form_labels)) + + ggplot2::labs(x = "year", y = NULL) + + ggplot2::theme_minimal(base_size = 11) + + ggplot2::theme(legend.position = "bottom") +} From f7ccf68e8f771c86dd56d2badc9324a0fb12a17c Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 06/30] split the fit scores by variable --- R/scores.R | 36 +++++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 13 deletions(-) diff --git a/R/scores.R b/R/scores.R index 14ebd54..abc1d83 100644 --- a/R/scores.R +++ b/R/scores.R @@ -15,12 +15,14 @@ crps_sample <- function(ens, obs) { ##' @author Akash BV ##' ##' @description Scores a J x P prediction ensemble G against the observations in -##' meta: RMSE and bias of the ensemble mean, coverage and mean width of the 90% -##' band, and the mean CRPS. Columns of G are matched to observations by slot. +##' meta, per variable: RMSE and bias of the ensemble mean, coverage and mean +##' width of the 90% band, and the mean CRPS. Columns of G are matched to +##' observations by slot. ##' ##' @param G matrix (J x P), columns named by observation slot. -##' @param meta observation meta (slot, value = obs mean). -##' @return one row tibble: rmse, bias, coverage, mean_width, crps. +##' @param meta observation meta (slot, variable, value = obs mean). +##' @return tibble, one row per variable: slots, rmse, bias, coverage, +##' mean_width, crps. ##' @export score_iteration <- function(G, meta) { G <- G[, meta$slot, drop = FALSE] @@ -28,15 +30,23 @@ score_iteration <- function(G, meta) { q05 <- apply(G, 2, stats::quantile, 0.05) q95 <- apply(G, 2, stats::quantile, 0.95) obs <- meta$value - crps <- mean(vapply(seq_along(obs), function(j) crps_sample(G[, j], obs[j]), - numeric(1))) - tibble::tibble( - rmse = sqrt(mean((m - obs)^2)), - bias = mean(m - obs), - coverage = mean(obs >= q05 & obs <= q95), - mean_width = mean(q95 - q05), - crps = crps - ) + crps <- vapply(seq_along(obs), function(j) crps_sample(G[, j], obs[j]), numeric(1)) + + # scored per variable: a pooled rmse over a joint target would average + # quantities in different units, which is not a quantity. `variable` also names + # the unit, so each row is internally consistent. + score_group <- function(idx) { + tibble::tibble( + slots = length(idx), + rmse = sqrt(mean((m[idx] - obs[idx])^2)), + bias = mean(m[idx] - obs[idx]), + coverage = mean(obs[idx] >= q05[idx] & obs[idx] <= q95[idx]), + mean_width = mean(q95[idx] - q05[idx]), + crps = mean(crps[idx]) + ) + } + dplyr::bind_rows(lapply(split(seq_along(obs), meta$variable), score_group), + .id = "variable") } ##' @title Per iteration score table From 48adf04d1d6f796daa395b4601ee1f614f0c22f8 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 07/30] refresh the nse globals --- R/calibration-package.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/R/calibration-package.R b/R/calibration-package.R index 5d364f3..984db99 100644 --- a/R/calibration-package.R +++ b/R/calibration-package.R @@ -4,9 +4,12 @@ # column names used non-standardly inside dplyr / ggplot2 verbs; declared so R CMD # check does not flag them as undefined globals. utils::globalVariables(c( - "variable", "sitename", "observation_level", "value", "study_year", + "variable", "sitename", "observation_level", "value", "obs_year", "treatment_id", "min_date", "max_date", "min_depth", "max_depth", - "cell_mean", "cell_sd", "n_rep", "var_mean", "var_obs", + "cell", "cell_period", "reported_se", "reported_units", + "target_variable", "target_units", "variance_model", + "obs_date_start", "obs_date_end", "var_obs", "cell_mean", "n_rep", "member", "slot", "q05", "q95", "m", "obs", "obs_sd", "stage", + "step", "lo", "hi", "bound", "label", "year", "form", "prior_mean", "prior_sd", "post_mean", "post_sd" )) From 9b3a82eaefa4a6de06abbaf32763eedc770c7f14 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 08/30] require the r version the scripts use --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index ef36d6e..185fd3c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,7 +14,7 @@ Description: Method-agnostic ensemble calibration of process-based ecosystem License: BSD_3_clause + file LICENSE Encoding: UTF-8 Depends: - R (>= 4.1.0) + R (>= 4.4.0) Imports: PEcAn.logger, PEcAn.priors, From 079559dcde2573181a5133e231886338d5757601 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 09/30] read targets from the config list --- scripts/010_prepare_observations.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scripts/010_prepare_observations.R b/scripts/010_prepare_observations.R index dd0fb9d..dbe74c1 100644 --- a/scripts/010_prepare_observations.R +++ b/scripts/010_prepare_observations.R @@ -11,8 +11,7 @@ config <- config::get(file = args$config) obs <- build_obs( cal_val_dir = file.path(config$scc, config$observations$dir), - target_var = config$observations$target_var, - sites = as.character(config$observations$sites) + targets = config$observations$targets ) dir.create(config$cache_dir, showWarnings = FALSE, recursive = TRUE) From b0e177a40899b3a3de75df397027ee3de3cfac03 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 10/30] contract raw observations into the fitted target --- scripts/012_build_target.R | 48 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 scripts/012_build_target.R diff --git a/scripts/012_build_target.R b/scripts/012_build_target.R new file mode 100644 index 0000000..3c76bac --- /dev/null +++ b/scripts/012_build_target.R @@ -0,0 +1,48 @@ +#!/usr/bin/env Rscript +# contract the raw curated observations into the fitted target. config$target is +# a list of contraction specs, each dispatched on its `type` (period_mean or +# contrast), so adding a variable or a site pair is a config edit, not a code +# edit. +# +# kept separate from 010 because 010's job is to read the curated data faithfully. +# what we choose to fit is a modeling decision and belongs where it can be read. + +library(calibration) + +args <- optparse::parse_args(optparse::OptionParser(option_list = list( + optparse::make_option(c("-c", "--config"), default = "config.yml", + help = "project config yaml [default: %default]") +))) +config <- config::get(file = args$config) + +if (length(config$target) == 0L) { + PEcAn.logger::logger.severe( + "config declares no target contractions; drop the target block to fit the ", + "raw slots, or list entries with type period_mean or contrast" + ) +} + +raw <- readRDS(file.path(config$cache_dir, "obs.rds")) + +contract <- function(tg) { + switch(tg$type, + period_mean = period_mean_contrast( + raw, variable = tg$variable, control = tg$control, + years = if (is.null(tg$years)) NULL else seq(tg$years[[1]], tg$years[[2]]), + new_variable = tg$new_variable + ), + contrast = contrast_target( + raw, variable = tg$variable, treatment = tg$treatment, control = tg$control, + new_variable = tg$new_variable + ), + PEcAn.logger::logger.severe("unknown target type '", tg$type, "'") + ) +} +obs <- do.call(bind_obs, lapply(config$target, contract)) + +out <- file.path(config$cache_dir, "target.rds") +saveRDS(obs, out) + +PEcAn.logger::logger.info(length(obs$y), " target slots across ", + dplyr::n_distinct(obs$meta$variable), " variables") +PEcAn.logger::logger.info("cached target -> ", out) From bc2d73e514f39284744570e027ba80d8bcc9093d Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 11/30] build the multisite settings from a blocks table --- scripts/015_build_settings.R | 215 +++++++++++++++++++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 scripts/015_build_settings.R diff --git a/scripts/015_build_settings.R b/scripts/015_build_settings.R new file mode 100644 index 0000000..920f8ce --- /dev/null +++ b/scripts/015_build_settings.R @@ -0,0 +1,215 @@ +#!/usr/bin/env Rscript + +# build the multisite PEcAn settings for a calibration run, plus a run dir +# default.param carrying the pinned parameters and one template initial +# condition per block. +# +# read a template, expand it with createMultiSiteSettings, fix up the per block +# bits with papply, write with write.settings. blocks come from a table, so +# adding or dropping a treatment is a data edit. + +library(PEcAn.settings) + +options <- list( + optparse::make_option(c("-c", "--config"), default = "config.yml", + help = "project config yaml"), + optparse::make_option("--blocks", default = NULL, + help = "block table; defaults to the config's forward$blocks"), + optparse::make_option("--workspace", default = NULL, + help = "run workspace; defaults to the config's forward$workspace"), + optparse::make_option("--obs", default = NULL, + help = "cached obs.rds; every fitted treatment must have a block") +) |> + purrr::modify(\(x) { x@help <- paste(x@help, "[default: %default]"); x }) + +args <- optparse::OptionParser(option_list = options) |> optparse::parse_args() +config <- config::get(file = args$config) + +blocks_file <- args$blocks %||% config$forward$blocks +workspace <- args$workspace %||% config$forward$workspace +obs_file <- args$obs %||% file.path(config$cache_dir, "obs.rds") + +# papply is chatty at DEBUG and drowns the real messages +PEcAn.logger::logger.setLevel("INFO") + +# absolute paths for everything the model reads or writes: PEcAn's restart code +# changes working directory, so a relative input silently resolves somewhere +# else. the one exception is host$modellauncher$binary, which stays +# ./scripts/sge_array_launcher.sh and resolves because qsub runs with -cwd from +# the workspace. +abs_path <- function(path) { + if (substr(path, 1, 1) != "/") path <- file.path(getwd(), path) + normalizePath(path, mustWork = FALSE) +} +workspace <- abs_path(workspace) +# prepared inputs are shared across runs, so a second workspace can point +# prepared_root at an existing one rather than duplicating them. +prepared_root <- abs_path(config$forward$prepared_root %||% workspace) + +blocks <- utils::read.csv(blocks_file, stringsAsFactors = FALSE) +required <- c("block_id", "lat", "lon", "veg_pft", "soil_pft", + "run_start", "run_end", "met_src", "ic_src", "events_src") +missing_cols <- setdiff(required, names(blocks)) +if (length(missing_cols) > 0L) { + PEcAn.logger::logger.severe("blocks table is missing: ", + paste(missing_cols, collapse = ", ")) +} +if (anyDuplicated(blocks$block_id) > 0L) { + PEcAn.logger::logger.severe( + "block_id must be unique; it becomes site$id and the observation operator ", + "matches a treatment's output on it. Duplicated: ", + paste(unique(blocks$block_id[duplicated(blocks$block_id)]), collapse = ", ") + ) +} + +obs <- readRDS(obs_file) +fitted <- unique(obs$meta$treatment_id) +orphans <- setdiff(fitted, blocks$block_id) +if (length(orphans) > 0L) { + PEcAn.logger::logger.severe( + length(orphans), " fitted treatment(s) have no run block: ", + paste(orphans, collapse = ", "), ". block_id must equal treatment_id." + ) +} +idle <- setdiff(blocks$block_id, fitted) +if (length(idle) > 0L) { + PEcAn.logger::logger.warn( + length(idle), " block(s) carry no fitted slot: ", paste(idle, collapse = ", ") + ) +} + +# createMultiSiteSettings copies every non-id column of this frame into run$site, +# so anything the per block fixups need is carried here rather than looked up again. +site_info <- data.frame( + id = blocks$block_id, + lat = blocks$lat, + lon = blocks$lon, + name = blocks$block_id, + veg_pft = blocks$veg_pft, + soil_pft = blocks$soil_pft, + block_start = blocks$run_start, + block_end = blocks$run_end, + met_src = blocks$met_src, + ic_src = blocks$ic_src, + events_src = blocks$events_src, + stringsAsFactors = FALSE +) + +settings <- read.settings(config$forward$template) + +## pinned parameters ----------------------------------------------------------- +# write.config.SIPNET reads default.param and then overwrites it row by row from +# the PFT trait samples, so pinning a parameter takes both halves: the value here +# and the trait dropped from the baseline sample (make_forward_sipnet's +# fixed_traits). one half alone silently loses to the PFT posterior. +default_param <- file.path(workspace, "sipnet.default.param") +# same version mapping write.config.SIPNET applies, so the pinned file is built +# from the template that run would otherwise have read +rev_num <- numeric_version(sub("^v", "", settings$model$revision, ignore.case = TRUE), + strict = FALSE) +if (is.na(rev_num)) { + PEcAn.logger::logger.severe("cannot parse model revision '", + settings$model$revision, "' as a version") +} +rev_str <- if (rev_num >= "2.0") "v2" else "v1" +stock_template <- system.file(paste0("template.param_", rev_str), package = "PEcAn.SIPNET") +if (!nzchar(stock_template)) { + PEcAn.logger::logger.severe("no stock template.param_", rev_str, + " for revision ", settings$model$revision) +} +param <- utils::read.table(stock_template, stringsAsFactors = FALSE) +for (nm in names(config$fixed_params)) { + sipnet_name <- config$fixed_params[[nm]]$sipnet + value <- config$fixed_params[[nm]]$value + hit <- param[[1]] == sipnet_name + if (!any(hit)) { + PEcAn.logger::logger.severe("no '", sipnet_name, "' row in ", stock_template) + } + param[hit, 2] <- value + PEcAn.logger::logger.info("pinned ", sipnet_name, " = ", value, " (", nm, ")") +} +dir.create(workspace, recursive = TRUE, showWarnings = FALSE) +# provision the workspace the launcher convention expects: modellauncher$binary is +# ./scripts/sge_array_launcher.sh relative to the workspace (see abs_path note above), +# and a workspace without it submits array jobs that die before the model runs. +launcher_src <- file.path(prepared_root, "scripts", "sge_array_launcher.sh") +launcher_dir <- file.path(workspace, "scripts") +if (!file.exists(file.path(launcher_dir, "sge_array_launcher.sh"))) { + if (!file.exists(launcher_src)) { + PEcAn.logger::logger.severe("no launcher at ", launcher_src, + "; the workspace cannot submit array jobs without it") + } + dir.create(launcher_dir, recursive = TRUE, showWarnings = FALSE) + file.copy(launcher_src, launcher_dir) + Sys.chmod(file.path(launcher_dir, "sge_array_launcher.sh"), "0755") + PEcAn.logger::logger.info("provisioned ", launcher_dir, "/sge_array_launcher.sh") +} +utils::write.table(param, default_param, quote = FALSE, sep = "\t", + row.names = FALSE, col.names = FALSE) +settings$model$default.param <- default_param +settings$model$binary <- abs_path(config$forward$binary) + +for (i in seq_along(settings$pfts)) { + settings$pfts[[i]]$posterior.files <- abs_path(settings$pfts[[i]]$posterior.files) + if (!is.null(settings$pfts[[i]]$outdir)) { + settings$pfts[[i]]$outdir <- abs_path(settings$pfts[[i]]$outdir) + } +} + +## per block fixups ------------------------------------------------------------ +# the blocks table declares the exact met, template ic, and events file per +# block, relative to the staging trees; the calibration pins each input to that +# one member, so no ic/met/events spread crosses the particles. a calibrated +# initial state overrides its pool per particle in the forward; a block without +# one runs the template ic as-is. +input_file <- function(root, rel, what, id) { + f <- file.path(root, rel) + if (!file.exists(f)) { + PEcAn.logger::logger.severe("block ", id, " ", what, " not found: ", f) + } + f +} + +set_block <- function(s) { + site <- s$run$site + # getRunSettings puts the dates inside run, alongside site and inputs + s$run$start.date <- site$block_start + s$run$end.date <- site$block_end + s$run$site$met.start <- site$block_start + s$run$site$met.end <- site$block_end + s$run$site$site.pft <- list(veg = site$veg_pft, soil = site$soil_pft) + s +} + +set_inputs <- function(s) { + site <- s$run$site + s$run$inputs$met$path <- list(path1 = input_file( + file.path(config$forward$sa_root, "inputs", "met"), site$met_src, "met", site$id)) + s$run$inputs$poolinitcond$path <- list(path1 = input_file( + file.path(config$forward$sa_root, "inputs", "IC"), site$ic_src, "ic", site$id)) + s$run$inputs$poolinitcond$ensemble <- 1 + s$run$inputs$events$path <- list(path1 = input_file( + file.path(prepared_root, "inputs"), site$events_src, "events", site$id)) + s +} + +settings <- settings |> + createMultiSiteSettings(site_info) |> + papply(set_block) |> + papply(set_inputs) + +settings$ensemble$size <- config$eki$n_particles +settings$ensemble$start.year <- as.integer(format(as.Date(min(blocks$run_start)), "%Y")) +settings$ensemble$end.year <- as.integer(format(as.Date(max(blocks$run_end)), "%Y")) + +out_dir <- file.path(workspace, "output") +settings$outdir <- out_dir +settings$modeloutdir <- file.path(out_dir, "out") +settings$rundir <- file.path(out_dir, "run") +settings$host$outdir <- file.path(out_dir, "out") +settings$host$rundir <- file.path(out_dir, "run") + +out_file <- file.path(workspace, "settings.xml") +write.settings(settings, outputfile = basename(out_file), outputdir = dirname(out_file)) +PEcAn.logger::logger.info("wrote ", out_file, ": ", nrow(blocks), " blocks, ", + "default.param ", default_param) From 3022c2772590b413d836869157fcb2110aa8c92a Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 12/30] build the shared prior from the settings pfts --- scripts/020_build_priors.R | 37 ++++++++++++++++++++++++++++++++----- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/scripts/020_build_priors.R b/scripts/020_build_priors.R index ecd1c4a..47a2045 100644 --- a/scripts/020_build_priors.R +++ b/scripts/020_build_priors.R @@ -10,13 +10,31 @@ args <- optparse::parse_args(optparse::OptionParser(option_list = list( ))) config <- config::get(file = args$config) -settings <- PEcAn.settings::read.settings( - file.path(config$scc, config$forward$run_dir, "settings.xml") +settings_file <- file.path( + config$forward$workspace %||% file.path(config$scc, config$forward$run_dir), + "settings.xml" +) +settings <- PEcAn.settings::read.settings(settings_file) +# config$soil_pft is a vector: the calibrated rates are one shared quantity written +# into every soil PFT. select explicitly and require every named PFT to be present. +soil_pft_names <- as.character(config$soil_pft) +pft_names <- vapply(settings$pfts, `[[`, character(1), "name") +missing_pfts <- setdiff(soil_pft_names, pft_names) +if (length(missing_pfts) > 0L) { + PEcAn.logger::logger.severe( + "soil PFT(s) named in the config are not in the settings: ", + paste(missing_pfts, collapse = ", "), "; settings carry ", + paste(pft_names, collapse = ", ") + ) +} +posterior_files <- stats::setNames( + vapply(settings$pfts[match(soil_pft_names, pft_names)], + `[[`, character(1), "posterior.files"), + soil_pft_names ) -soil_pft <- Filter(function(p) p$name == config$soil_pft, settings$pfts)[[1]] prior <- c( - prior_from_postdistns(config$priors$post_distns_params, soil_pft$posterior.files), + prior_from_shared_postdistns(config$priors$post_distns_params, posterior_files), prior_from_specs(config$priors$specified) ) @@ -25,12 +43,21 @@ state <- config$priors$state if (!is.null(state)) { obs <- readRDS(file.path(config$cache_dir, "obs.rds")) prior <- c(prior, state_prior_from_obs( - obs$meta, prefix = state$prefix, + obs$meta, prefix = state$prefix, variable = state$variable, from_unit = state$from_unit, to_unit = state$to_unit, anchor_year = state$anchor_year )) } +# a parameter declared in two sources would sample two columns under one name +dup <- names(prior)[duplicated(names(prior))] +if (length(dup) > 0L) { + PEcAn.logger::logger.severe( + "parameter(s) declared by more than one prior source: ", + paste(unique(dup), collapse = ", ") + ) +} + out <- file.path(config$cache_dir, "prior.rds") saveRDS(prior, out) PEcAn.logger::logger.info("cached prior (", length(prior), " params) -> ", out) From 318c94bff94c166406b54ee8d4ddfae79cbde274 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 13/30] run the calibration and keep the raw predictions --- scripts/030_calibrate.R | 129 +++++++++++++++++++++++++++++++--------- 1 file changed, 100 insertions(+), 29 deletions(-) diff --git a/scripts/030_calibrate.R b/scripts/030_calibrate.R index 84eda51..4b61d28 100644 --- a/scripts/030_calibrate.R +++ b/scripts/030_calibrate.R @@ -6,52 +6,123 @@ library(calibration) args <- optparse::parse_args(optparse::OptionParser(option_list = list( optparse::make_option(c("-c", "--config"), default = "config.yml", - help = "project config yaml [default: %default]") + help = "project config yaml [default: %default]"), + optparse::make_option("--dry-run", action = "store_true", default = FALSE, + help = "one iteration, result written to --out [default: %default]"), + optparse::make_option("--particles", type = "integer", default = NULL, + help = "override ensemble size [default: config]"), + optparse::make_option("--iterations", type = "integer", default = NULL, + help = "override tempering steps [default: config]"), + optparse::make_option("--out", default = NULL, + help = "output root; defaults to the settings outdir [default: %default]") ))) config <- config::get(file = args$config) options(warn = 1) -options(error = quote({ - if (!interactive()) q(status = 1) -})) -obs <- readRDS(file.path(config$cache_dir, "obs.rds")) +# when the config declares a target transform, 012 has written the fitted quantity +# to target.rds and that is what gets calibrated. obs.rds is the raw cache 012 +# reads from, so fitting it would silently calibrate the raw per-treatment cells +# instead of the contractions the transform produced. a missing target.rds means +# 012 has not been run against this config. +target_file <- file.path(config$cache_dir, "target.rds") +if (!is.null(config$target)) { + if (!file.exists(target_file)) { + PEcAn.logger::logger.severe( + "config declares a target transform but ", target_file, " does not exist; ", + "run 012_build_target.R against this config first" + ) + } + obs <- readRDS(target_file) + # the model is harvested on the raw slots and put through the same transform + raw_obs <- readRDS(file.path(config$cache_dir, "obs.rds")) + raw_obs$transform <- obs$transform + PEcAn.logger::logger.info("fitting the transformed target: ", nrow(obs$meta), " slots (", + paste(unique(obs$meta$variable), collapse = ", "), ")") +} else { + obs <- readRDS(file.path(config$cache_dir, "obs.rds")) + raw_obs <- NULL +} prior <- readRDS(file.path(config$cache_dir, "prior.rds")) -# fit only study years at or after fit_from_study_year. the earliest year stays -# the initial condition; the state prior in 020 anchors there, so it is dropped -# from the fitted target here but not from the run. -fit_from <- if (is.null(config$fit_from_study_year)) { - min(obs$meta$study_year) -} else { - config$fit_from_study_year +# drop each site's establishment years: config$fit maps sitename -> the first +# observation year to fit, so the pre-fit years stay the initial condition (the state +# prior in 020 anchors there) but are not fitted. study-year scales differ across +# sites, so the rule is per site; a site absent from config$fit fits every year. +unknown <- setdiff(names(config$fit), obs$meta$sitename) +if (length(unknown) > 0L) { + PEcAn.logger::logger.severe("config$fit names site(s) with no slots: ", + paste(unknown, collapse = ", ")) } -keep <- obs$meta$slot[obs$meta$study_year >= fit_from] -obs$y <- obs$y[keep] -obs$Sigma <- obs$Sigma[keep, keep, drop = FALSE] -obs$meta <- obs$meta[obs$meta$study_year >= fit_from, ] -PEcAn.logger::logger.info("fitting study_year >= ", fit_from, ": ", length(obs$y), " slots") - -settings <- PEcAn.settings::read.settings( - file.path(config$scc, config$forward$run_dir, "settings.xml") +keep <- rep(TRUE, nrow(obs$meta)) +for (site in names(config$fit)) { + keep <- keep & !(obs$meta$sitename == site & obs$meta$obs_year < config$fit[[site]]) +} +obs <- subset_obs(obs, keep) +PEcAn.logger::logger.info(length(obs$y), " slots after per-site establishment drop") + +# variables held out of the likelihood but still predicted and reported. a target +# the model cannot reach cannot be calibrated through: the estimator answers by +# pushing a parameter to its bound for reasons unrelated to the process it stands +# for. held out, the variable is still simulated and scored. +obs_all <- obs +val_vars <- config$validation_variables +if (!is.null(val_vars)) { + is_val <- obs$meta$variable %in% as.character(val_vars) + if (!any(is_val)) { + PEcAn.logger::logger.severe( + "validation_variables match no observation: ", paste(val_vars, collapse = ", ") + ) + } + obs <- subset_obs(obs_all, !is_val) + PEcAn.logger::logger.info( + "fitting ", length(obs$y), " slots; holding ", sum(is_val), + " out as validation (", paste(unique(obs_all$meta$variable[is_val]), collapse = ", "), ")" + ) +} + +settings_file <- file.path( + config$forward$workspace %||% file.path(config$scc, config$forward$run_dir), + "settings.xml" ) +settings <- PEcAn.settings::read.settings(settings_file) + +n_particles <- args$particles %||% config$eki$n_particles +n_iterations <- args$iterations %||% (if (args$`dry-run`) 1L else config$eki$n_iterations) +base_out_dir <- args$out %||% settings$outdir forward <- make_forward_sipnet( - settings = settings, obs = obs, n_particles = config$eki$n_particles, - harvest_var = config$forward$harvest_var, - from_unit = config$forward$from_unit, to_unit = config$forward$to_unit, - soil_pft = config$soil_pft, - state_prefix = config$priors$state$prefix, - state_pool = config$forward$state_pool + settings = settings, obs = obs, n_particles = n_particles, + var_map = config$forward$var_map, + soil_pfts = as.character(config$soil_pft), + state_prefix = config$priors$state$prefix %||% "soilInit.", + state_pool = config$forward$state_pool, + fixed_traits = names(config$fixed_params), + base_out_dir = base_out_dir, + raw_obs = raw_obs ) control <- calibration_control( - method = "eki", n_particles = config$eki$n_particles, - n_iterations = config$eki$n_iterations, seed = config$eki$seed + method = "eki", n_particles = n_particles, + n_iterations = n_iterations, seed = config$eki$seed ) result <- calibrate(obs, prior, forward, control) +result$obs_all <- obs_all + +# validation and figure predictions cost no extra model runs: itr1 holds the +# prior forward and itr(n + 1) the posterior forward, both already on disk. +window <- run_window(settings) +val_meta <- if (is.null(raw_obs)) obs_all$meta else raw_obs$meta +itr_dir <- function(k) file.path(base_out_dir, paste0("itr", k), "out") +result$raw_meta <- val_meta +result$G_raw_prior <- harvest_output_to_G(itr_dir(1L), val_meta, + config$forward$var_map, window) +result$G_raw_post <- harvest_output_to_G(itr_dir(n_iterations + 1L), val_meta, + config$forward$var_map, window) +result$G_validation <- if (is.null(obs_all$transform)) result$G_raw_post else + apply_transform(result$G_raw_post, obs_all$transform) -out <- file.path(config$cache_dir, "result.rds") +out <- file.path(if (args$`dry-run`) base_out_dir else config$cache_dir, "result.rds") saveRDS(result, out) PEcAn.logger::logger.info("cached result -> ", out) From 7015bc93bc14c706cbb15dff60ad8fc66ccdfe48 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 14/30] write figures and scores per variable --- scripts/040_plot.R | 99 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 87 insertions(+), 12 deletions(-) diff --git a/scripts/040_plot.R b/scripts/040_plot.R index 14bc2a0..93bff5f 100644 --- a/scripts/040_plot.R +++ b/scripts/040_plot.R @@ -1,28 +1,103 @@ #!/usr/bin/env Rscript -# plot and score the calibration result +# plot and score the calibration result. +# +# figures go to the run workspace, not the repo: they are run artifacts and belong +# beside the result they were made from. nothing here writes a caption -- the +# figure carries labs(x, y) only and the story goes in the writeup. library(calibration) args <- optparse::parse_args(optparse::OptionParser(option_list = list( optparse::make_option(c("-c", "--config"), default = "config.yml", - help = "project config yaml [default: %default]") + help = "project config yaml [default: %default]"), + optparse::make_option("--result", default = NULL, + help = "result.rds; defaults to the config cache [default: %default]"), + optparse::make_option("--figdir", default = NULL, + help = "output directory [default: /figures]") ))) config <- config::get(file = args$config) -result <- readRDS(file.path(config$cache_dir, "result.rds")) -dir.create("figures", showWarnings = FALSE, recursive = TRUE) +result_file <- args$result %||% file.path(config$cache_dir, "result.rds") +figdir <- args$figdir %||% file.path( + config$forward$workspace %||% file.path(config$scc, config$forward$run_dir), + "figures" +) +result <- readRDS(result_file) +dir.create(figdir, showWarnings = FALSE, recursive = TRUE) -# per iteration prediction matrices, the tempered forwards through the posterior +# the tempered forwards, prior through posterior G_list <- c(lapply(result$eki$eki_list, function(s) s$G), list(result$eki$G)) names(G_list) <- paste("iteration", seq_along(G_list)) -save_iterations_pdf(G_list, result$obs$meta, - file.path("figures", "ensembles_vs_truth.pdf"), - ylab = config$observations$target_var) +# ensembles against observations, one page per iteration, site, and variable: +# a page shares one axis, so it holds one unit. +meta <- result$obs$meta +groups <- unique(meta[c("sitename", "variable")]) +grDevices::pdf(file.path(figdir, "ensembles_vs_truth.pdf"), width = 9, height = 7) +for (nm in names(G_list)) { + for (g in seq_len(nrow(groups))) { + sm <- meta[meta$sitename == groups$sitename[g] & + meta$variable == groups$variable[g], , drop = FALSE] + plot_fn <- if (all(is.na(sm$obs_year))) plot_target_slots else plot_ensembles_vs_truth + print(plot_fn( + G_list[[nm]][, sm$slot, drop = FALSE], sm, + ylab = unique(sm$units), + title = paste(groups$sitename[g], groups$variable[g], nm) + )) + } +} +grDevices::dev.off() plot_param_densities(list(prior = result$U0, posterior = result$U), - file.path("figures", "param_prior_posterior.pdf")) + file.path(figdir, "param_prior_posterior.pdf")) +readr::write_csv(param_shift(result$U0, result$U), + file.path(figdir, "param_shift.csv")) -scores <- score_table(G_list, result$obs$meta) -readr::write_csv(scores, file.path("figures", "scores.csv")) -PEcAn.logger::logger.info("wrote figures/ (ensembles_vs_truth, param_prior_posterior, scores)") +# parameter trace: does a parameter settle inside its support, or move onto a bound, +# and does it walk there or jump. eki_list holds the update in unconstrained space, +# so it is mapped back before plotting. +pm <- result$eki$par_map +trace <- c(list(result$U0), lapply(result$eki$eki_list, function(s) pm$inv(s$U))) +prior <- readRDS(file.path(config$cache_dir, "prior.rds")) +ggplot2::ggsave(file.path(figdir, "param_trace.png"), + plot_param_trace(trace, prior), + width = 8, height = 3.5, dpi = 300) + +# measured against modeled treatment effect, before and after calibration +fe <- config$figures$treatment_effect +if (!is.null(fe)) { + if (is.null(result$G_raw_prior)) { + PEcAn.logger::logger.info("skipping treatment effect figure: this result ", + "carries no raw harvests") + } else { + ggplot2::ggsave(file.path(figdir, "treatment_effect.png"), + plot_treatment_effect(result$G_raw_prior, result$G_raw_post, + result$raw_meta, variable = fe$variable, + treatment = fe$treatment, control = fe$control), + width = 9, height = 4, dpi = 300) + } +} + +# score_iteration already splits by variable, so this is iteration x variable: a +# target that degrades while another improves stays visible rather than averaged away +scores <- score_table(G_list, meta) +readr::write_csv(scores, file.path(figdir, "scores.csv")) + +# held-out validation, when the config keeps variables out of the likelihood +val_vars <- as.character(config$validation_variables) +if (length(val_vars) > 0L) { + if (is.null(result$G_validation)) { + PEcAn.logger::logger.info("skipping validation figure: this result carries ", + "no validation predictions") + } else { + vm <- result$obs_all$meta + vm <- vm[vm$variable %in% val_vars, , drop = FALSE] + ggplot2::ggsave(file.path(figdir, "validation_predicted_vs_observed.png"), + plot_validation(result$G_validation, vm), + width = 8, height = 4, dpi = 300) + readr::write_csv(score_iteration(result$G_validation, vm), + file.path(figdir, "scores_validation.csv")) + } +} + +PEcAn.logger::logger.info("wrote ", figdir) From b8a5e4d085937a679f7aeff6a7bb621b4f48ccba Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 15/30] test the target operators --- tests/testthat/test_observation_operator.R | 216 +++++++++++++++++++++ 1 file changed, 216 insertions(+) create mode 100644 tests/testthat/test_observation_operator.R diff --git a/tests/testthat/test_observation_operator.R b/tests/testthat/test_observation_operator.R new file mode 100644 index 0000000..044c9a5 --- /dev/null +++ b/tests/testthat/test_observation_operator.R @@ -0,0 +1,216 @@ +# Five years, three treatments, none of them a rigid offset of another: a contrast with +# no year to year spread is degenerate and would make the covariance singular for a +# reason the real data does not have. +fake_series_obs <- function() { + yrs <- 2005:2009 + trts <- c("c", "a", "b") + vals <- c(10, 12, 11, 14, 13, + 16, 17, 18, 19, 20, + 15, 16, 14, 18, 16) + meta <- tibble::tibble( + slot = paste0("SOC__", rep(trts, each = length(yrs)), "__", rep(yrs, length(trts))), + variable = "SOC", sitename = "salinas", + treatment_id = rep(trts, each = length(yrs)), + obs_year = rep(yrs, length(trts)), + min_date = paste0(rep(yrs, length(trts)), "-10-01"), + max_date = paste0(rep(yrs, length(trts)), "-10-01"), + min_depth = 0, max_depth = 30, + units = "Mg C ha-1", observation_level = "replicate", n_rep = 4L, + value = vals, var_obs = 1 + ) + Sigma <- diag(1, nrow(meta)); dimnames(Sigma) <- list(meta$slot, meta$slot) + list(y = stats::setNames(vals, meta$slot), Sigma = Sigma, meta = meta) +} + +series <- function(t) fake_series_obs()$y[paste0("SOC__", t, "__", 2005:2009)] + +test_that("period_mean_contrast returns one level plus one contrast per treatment", { + out <- period_mean_contrast(fake_series_obs(), "SOC", control = "c") + expect_equal(length(out$y), 3L) + expect_equal(out$meta$treatment_id, c("c", "a_vs_c", "b_vs_c")) + expect_equal(unname(out$y[1]), mean(series("c"))) + expect_equal(unname(out$y[2]), mean(series("a") - series("c"))) + expect_true(all(out$meta$variable == "SOC_periodmean")) + expect_equal(unique(out$meta$n_rep), 5L) +}) + +test_that("the transform reproduces the fitted values from the raw slots", { + obs <- fake_series_obs() + out <- period_mean_contrast(obs, "SOC", control = "c") + raw <- obs$y[colnames(out$transform)] + expect_equal(as.numeric(out$transform %*% raw), as.numeric(out$y)) +}) + +test_that("a constant added to every year moves the level and leaves contrasts alone", { + obs <- fake_series_obs() + out <- period_mean_contrast(obs, "SOC", control = "c") + shifted <- obs; shifted$y <- shifted$y + 7 + out2 <- period_mean_contrast(shifted, "SOC", control = "c") + expect_equal(unname(out2$y[1]), unname(out$y[1]) + 7) + expect_equal(unname(out2$y[-1]), unname(out$y[-1])) +}) + +test_that("uncertainty is the standard error of the period mean across years", { + out <- period_mean_contrast(fake_series_obs(), "SOC", control = "c") + n <- 5 + expect_equal(sqrt(out$Sigma[1, 1]), stats::sd(series("c")) / sqrt(n)) + expect_equal(sqrt(out$Sigma[2, 2]), + stats::sd(series("a") - series("c")) / sqrt(n)) + # and it is the paired spread, not the sum of the two levels' spreads + expect_lt(out$Sigma[2, 2], + stats::var(series("a")) / n + stats::var(series("c")) / n) +}) + +test_that("the covariance is positive definite so the EnKF can factor it", { + out <- period_mean_contrast(fake_series_obs(), "SOC", control = "c") + expect_true(isSymmetric(unname(out$Sigma))) + expect_gt(min(eigen(out$Sigma, symmetric = TRUE, only.values = TRUE)$values), 0) + expect_false(inherits(try(chol(out$Sigma), silent = TRUE), "try-error")) +}) + +test_that("shrinking to positive definiteness leaves the marginal variances alone", { + S <- matrix(c(4, 2, -4, + 2, 1, -2, + -4, -2, 4), 3, 3) # rank 1, singular by construction + out <- calibration:::.shrink_to_pd(S, label = "test") + expect_equal(diag(out), diag(S)) + expect_gt(min(eigen(out, symmetric = TRUE, only.values = TRUE)$values), 0) + # off-diagonals are damped toward zero, never sign flipped + expect_true(all(sign(out[upper.tri(out)]) == sign(S[upper.tri(S)]))) + expect_true(all(abs(out[upper.tri(out)]) <= abs(S[upper.tri(S)]))) +}) + +test_that("an unbalanced panel is refused rather than averaged", { + obs <- fake_series_obs() + drop <- obs$meta$treatment_id == "a" & obs$meta$obs_year == 2007 + obs$meta <- obs$meta[!drop, ] + obs$y <- obs$y[obs$meta$slot] + obs$Sigma <- obs$Sigma[obs$meta$slot, obs$meta$slot, drop = FALSE] + expect_error(period_mean_contrast(obs, "SOC", control = "c"), + "exactly one slot per treatment per year") +}) + +test_that("too short a period is refused, since the mean has no estimable error", { + expect_error( + period_mean_contrast(fake_series_obs(), "SOC", control = "c", years = 2005:2006), + "needs at least" + ) +}) + +test_that("an absent control is an error naming what is present", { + expect_error(period_mean_contrast(fake_series_obs(), "SOC", control = "nope"), + "Present: ") +}) + +test_that("contrast_target pairs treatments by date and refuses unpaired dates", { + meta <- tibble::tibble( + slot = c("n__t__d1", "n__t__d2", "n__c__d1", "n__c__d2"), + variable = "n2o", sitename = "s", + treatment_id = c("t", "t", "c", "c"), + obs_year = 2019L, + min_date = c("2019-01-01", "2019-02-01", "2019-01-01", "2019-02-01"), + max_date = c("2019-01-01", "2019-02-01", "2019-01-01", "2019-02-01"), + units = "g ha-1 day-1", value = c(3, 5, 1, 4), var_obs = c(1, 1, 1, 1) + ) + S <- diag(1, 4); dimnames(S) <- list(meta$slot, meta$slot) + obs <- list(y = stats::setNames(meta$value, meta$slot), Sigma = S, meta = meta) + + ct <- contrast_target(obs, "n2o", treatment = "t", control = "c", + new_variable = "n2o_c") + expect_equal(unname(ct$y), c(2, 1)) + # independent cells: variances add + expect_equal(unname(diag(ct$Sigma)), c(2, 2)) + # the difference form records its contraction + expect_equal(as.numeric(ct$transform %*% obs$y[colnames(ct$transform)]), + as.numeric(ct$y)) + + short <- subset_obs(obs, obs$meta$slot != "n__c__d2") + expect_error(contrast_target(short, "n2o", "t", "c"), "cannot form a") +}) + +test_that("bind_obs stacks targets block diagonally and refuses duplicate slots", { + obs <- fake_series_obs() + a <- period_mean_contrast(obs, "SOC", control = "c", new_variable = "pm1") + b <- period_mean_contrast(obs, "SOC", control = "c", new_variable = "pm2") + both <- bind_obs(a, b) + expect_equal(length(both$y), length(a$y) + length(b$y)) + # off blocks are zero, within blocks carried through + expect_equal(unname(both$Sigma[a$meta$slot, b$meta$slot]), + matrix(0, nrow(a$meta), nrow(b$meta))) + expect_equal(unname(both$Sigma[a$meta$slot, a$meta$slot]), unname(a$Sigma)) + expect_error(bind_obs(a, a), "duplicate slot") +}) + +# the operator reads real model output, so these mock read.output and assert on +# the window it is asked for and on what happens at the window edge. + +make_meta <- function(slot, treat, a, b, variable = "SOC_stock") { + tibble::tibble(slot = slot, variable = variable, treatment_id = treat, + min_date = a, max_date = b) +} + +test_that("each treatment is harvested over its own run window", { + skip_if_not_installed("PEcAn.utils") + out_root <- withr::local_tempdir() + dir.create(file.path(out_root, "ENS-00001-early")) + dir.create(file.path(out_root, "ENS-00001-late")) + + asked <- list() + fake_read <- function(runid, outdir, start.year, end.year, variables, ...) { + asked[[sub("^ENS-[0-9]+-", "", runid)]] <<- c(start.year, end.year) + days <- seq(as.Date(paste0(start.year, "-01-01")), + as.Date(paste0(end.year, "-12-31")), by = "day") + data.frame(posix = days, TotSoilCarb = rep(start.year / 1000, length(days))) + } + testthat::local_mocked_bindings(read.output = fake_read, .package = "PEcAn.utils") + + meta <- rbind( + make_meta("early_slot", "early", "2005-06-01", "2005-06-01"), + make_meta("late_slot", "late", "2020-06-01", "2020-06-01") + ) + win <- cbind(early = c(2005L, 2011L), late = c(2017L, 2023L)) + vm <- list(SOC_stock = list(model_var = "TotSoilCarb", from = "kg/m2", + to = "kg/m2")) + + G <- harvest_output_to_G(out_root, meta, vm, win) + + expect_equal(asked$early, c(2005L, 2011L)) + expect_equal(asked$late, c(2017L, 2023L)) # not the first block's window + expect_equal(unname(G[1, "early_slot"]), 2.005) + expect_equal(unname(G[1, "late_slot"]), 2.017) +}) + +test_that("an observation outside the run window fails instead of returning a neighbor", { + skip_if_not_installed("PEcAn.utils") + out_root <- withr::local_tempdir() + dir.create(file.path(out_root, "ENS-00001-site")) + + fake_read <- function(runid, outdir, start.year, end.year, variables, ...) { + days <- seq(as.Date("2017-01-01"), as.Date("2023-12-31"), by = "day") + data.frame(posix = days, TotSoilCarb = seq_along(days) / 1000) + } + testthat::local_mocked_bindings(read.output = fake_read, .package = "PEcAn.utils") + + win <- cbind(site = c(2017L, 2023L)) + vm <- list(SOC_stock = list(model_var = "TotSoilCarb", from = "kg/m2", + to = "kg/m2")) + + inside <- make_meta("in_slot", "site", "2020-06-01", "2020-06-01") + expect_silent(harvest_output_to_G(out_root, inside, vm, win)) + + # 2024 is past the run end: nearest-date substitution would answer with + # 2023-12-31 and look clean + outside <- make_meta("out_slot", "site", "2024-06-01", "2024-06-01") + expect_error(harvest_output_to_G(out_root, outside, vm, win), + "outside the") +}) + +test_that("apply_transform is the same linear map the observations went through", { + obs <- fake_series_obs() + out <- period_mean_contrast(obs, "SOC", control = "c") + G <- rbind(obs$y, obs$y * 2) + fitted <- apply_transform(G, out$transform) + expect_equal(unname(fitted[1, ]), unname(out$y)) + expect_equal(colnames(fitted), out$meta$slot) +}) + From 2ab5b8555f23cc902bfbaada7580bf552204c9f7 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 16/30] test the shared rate injection and run windows --- tests/testthat/test_forward_sipnet.R | 38 ++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 tests/testthat/test_forward_sipnet.R diff --git a/tests/testthat/test_forward_sipnet.R b/tests/testthat/test_forward_sipnet.R new file mode 100644 index 0000000..2cb657a --- /dev/null +++ b/tests/testthat/test_forward_sipnet.R @@ -0,0 +1,38 @@ +# a shared soil rate reaching only one PFT returns a plausible number rather +# than an error, so these assert the write reaches every named PFT and that an +# absent one is loud. + +test_that("inject_traits writes the shared rate into every soil PFT", { + baseline <- list( + soil = data.frame(som_respiration_rate = rep(0.01, 3), kCN = rep(80, 3)), + soil_rice = data.frame(som_respiration_rate = rep(0.01, 3), kCN = rep(80, 3)), + soil_nfixer = data.frame(som_respiration_rate = rep(0.01, 3), kCN = rep(80, 3)), + annual_crop = data.frame(SLA = rep(19, 3)) + ) + U <- matrix(c(0.02, 0.03, 0.04), ncol = 1, + dimnames = list(NULL, "som_respiration_rate")) + out <- inject_traits(baseline, c("soil", "soil_rice", "soil_nfixer"), U) + + for (pft in c("soil", "soil_rice", "soil_nfixer")) { + expect_equal(out[[pft]]$som_respiration_rate, c(0.02, 0.03, 0.04)) + expect_equal(out[[pft]]$kCN, rep(80, 3)) # uncalibrated column untouched + } + expect_equal(out$annual_crop$SLA, rep(19, 3)) # veg PFT untouched +}) + +test_that("inject_traits fails when a named soil PFT is absent", { + baseline <- list(soil = data.frame(som_respiration_rate = rep(0.01, 2))) + U <- matrix(0.02, nrow = 2, ncol = 1, + dimnames = list(NULL, "som_respiration_rate")) + expect_error(inject_traits(baseline, c("soil", "soil_rice"), U), "soil_rice") +}) + +test_that("run_window reads each treatment's own years from the settings", { + fake <- function(id, a, b) list(run = list(site = list(id = id), + start.date = a, end.date = b)) + settings <- list(fake("early", "2005-01-01", "2011-12-31"), + fake("late", "2017-01-01", "2023-12-31")) + w <- run_window(settings) + expect_equal(w[, "early"], c(2005L, 2011L)) + expect_equal(w[, "late"], c(2017L, 2023L)) +}) From 3bf13b0d28059c8ddd3a2c1a5bf1f2cc77f0f77e Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 17/30] test the shared prior and state anchor --- tests/testthat/test_priors.R | 41 ++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 tests/testthat/test_priors.R diff --git a/tests/testthat/test_priors.R b/tests/testthat/test_priors.R new file mode 100644 index 0000000..0be2ac9 --- /dev/null +++ b/tests/testthat/test_priors.R @@ -0,0 +1,41 @@ +# one calibrated rate shared across PFTs needs one prior, and a state anchored +# on a joint meta must scope to its variable; both failure modes are silent +# without these contracts. + +test_that("a shared prior is returned only when the soil PFTs agree", { + skip_if_not_installed("PEcAn.priors") + make_post <- function(path, parama) { + post.distns <- data.frame(distn = "weibull", parama = parama, paramb = 0.0112, + row.names = "som_respiration_rate") + save(post.distns, file = path) + path + } + d <- withr::local_tempdir() + agree <- c(soil = make_post(file.path(d, "a.Rdata"), 2.21), + soil_rice = make_post(file.path(d, "b.Rdata"), 2.21)) + pr <- prior_from_shared_postdistns("som_respiration_rate", agree) + expect_equal(pr$som_respiration_rate$parama, 2.21) + + differ <- c(soil = make_post(file.path(d, "c.Rdata"), 2.21), + soil_rice = make_post(file.path(d, "e.Rdata"), 9.99)) + expect_error(prior_from_shared_postdistns("som_respiration_rate", differ), + "differs between soil PFTs") +}) + +test_that("state_prior_from_obs anchors on one variable, not the whole joint meta", { + meta <- tibble::tibble( + slot = paste0("s", 1:4), + variable = c("SOC_stock", "SOC_stock", "N2O_flux", "N2O_flux"), + treatment_id = c("t1", "t2", "t1", "t2"), + obs_year = c(2005, 2005, 2005, 2005), + value = c(40, 44, 0.3, 0.5), + var_obs = c(4, 4, 0.01, 0.01) + ) + pr <- state_prior_from_obs(meta, prefix = "soilInit.", from_unit = "Mg/ha", + to_unit = "kg/m2", variable = "SOC_stock", + anchor_year = 2005) + expect_length(pr, 2) # only the two SOC treatments + expect_named(pr, c("soilInit.t1", "soilInit.t2")) + expect_error(state_prior_from_obs(meta, "soilInit.", "Mg/ha", "kg/m2", + variable = "AbvGrndWood"), "AbvGrndWood") +}) From a58feced498513d6e502bb4868eeb077a660ca86 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 18/30] move the salinas config to a numbered directory and repair it --- .../{salinas_soc => 1_salinas_soc}/config.yml | 30 +++++++++---------- 1 file changed, 15 insertions(+), 15 deletions(-) rename examples/{salinas_soc => 1_salinas_soc}/config.yml (56%) diff --git a/examples/salinas_soc/config.yml b/examples/1_salinas_soc/config.yml similarity index 56% rename from examples/salinas_soc/config.yml rename to examples/1_salinas_soc/config.yml index 6d1671e..c0d8875 100644 --- a/examples/salinas_soc/config.yml +++ b/examples/1_salinas_soc/config.yml @@ -1,17 +1,21 @@ # salinas organic cropping systems soil carbon calibration. this is the worked # example that runs the package on one real dataset; the soc and salinas # specifics live here, not in R/. run the scripts with -# --config examples/salinas_soc/config.yml. +# --config examples/1_salinas_soc/config.yml. default: scc: /projectnb/dietzelab/ccmmf - cache_dir: cache + cache_dir: /projectnb/dietzelab/ccmmf/usr/akash/salinas_socs/cache soil_pft: soil observations: - dir: usr/akash/cal_val - target_var: SOC_stock_Mg_ha - sites: - - salinas_socs + dir: usr/akash/cal-val-data + targets: + # 2003-2004 is the establishment drop, a one time disturbance the model + # should not be pushed to reproduce; the fit starts at the 2005 stock + - variable: SOC_stock_Mg_ha + sites: [salinas_socs] + units: Mg C ha-1 + years: [2005, 2011] priors: # traits with a soil pft meta-analysis posterior are read from it. @@ -21,25 +25,21 @@ default: specified: soil_respiration_Q10: {distn: unif, parama: 1.4, paramb: 3.0} turn_over_time: {distn: unif, parama: 0.13, paramb: 1.2} - # per site initial soil carbon, anchored to the year 2 observation (the first - # fitted year), converted from the reported Mg/ha to the model kg/m2. + # per system initial soil carbon, anchored to the earliest fitted observation + # (2005), converted from the reported Mg/ha to the model kg/m2. state: prefix: soilInit. + variable: SOC_stock_Mg_ha from_unit: Mg/ha to_unit: kg/m2 - anchor_year: 2 forward: run_dir: usr/akash/salinas_socs - harvest_var: TotSoilCarb - from_unit: kg/m2 - to_unit: Mg/ha + var_map: + SOC_stock_Mg_ha: {model_var: TotSoilCarb, from: kg/m2, to: Mg/ha} state_pool: soil_organic_carbon_content eki: n_particles: 50 n_iterations: 3 seed: 556688 - - # fit study_year 2 onward; year 0-1 is not a target. - fit_from_study_year: 2 From 593557541de1a917b3401b457a5bc820286519da Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 19/30] move the salinas readme and match it to the current mechanism --- .../{salinas_soc => 1_salinas_soc}/README.md | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) rename examples/{salinas_soc => 1_salinas_soc}/README.md (58%) diff --git a/examples/salinas_soc/README.md b/examples/1_salinas_soc/README.md similarity index 58% rename from examples/salinas_soc/README.md rename to examples/1_salinas_soc/README.md index 832f990..0327229 100644 --- a/examples/salinas_soc/README.md +++ b/examples/1_salinas_soc/README.md @@ -9,13 +9,13 @@ systems. everything soc and salinas specific lives here; `R/` stays generic. from the package root, with the pecan modules loaded: ``` -Rscript scripts/010_prepare_observations.R --config examples/salinas_soc/config.yml -Rscript scripts/020_build_priors.R --config examples/salinas_soc/config.yml -Rscript scripts/030_calibrate.R --config examples/salinas_soc/config.yml -Rscript scripts/040_plot.R --config examples/salinas_soc/config.yml +Rscript scripts/010_prepare_observations.R --config examples/1_salinas_soc/config.yml +Rscript scripts/020_build_priors.R --config examples/1_salinas_soc/config.yml +Rscript scripts/030_calibrate.R --config examples/1_salinas_soc/config.yml +Rscript scripts/040_plot.R --config examples/1_salinas_soc/config.yml ``` -010 caches the target, 020 the prior, 030 runs the ensemble kalman inversion and +010 caches the observations, 020 the prior, 030 runs the ensemble kalman inversion and caches the result, 040 writes the figures and scores. ## what it calibrates @@ -25,24 +25,26 @@ caches the result, 040 writes the figures and scores. - `soil_respiration_Q10` and `turn_over_time` from the soil pft bety priors, set in the config, not tuned to the data. - per system initial soil carbon (`soilInit.`), anchored to each system's - year 2 observation and free to update. + 2005 observation and free to update. launching is left to pecan and the prepared host block in the forward run's settings (qsub, sge_array_launcher.sh, Njobmax, qstat), used exactly as written. ## what we fit, and why -- **fitting window is study_year 2 to 8.** the measured soil carbon drops sharply - in the first year, between 2003 and 2004, as the cropping systems are - established. that first year drop is a one-time disturbance, not gradual - decomposition, and a first order soil model should not be pushed to reproduce a - change that large in a single step. so year 0 to 1 is left out of the fitted - target: the model is initialized at the year-2 stock and the fit runs from - study_year 2 onward. -- **initial soil carbon** (`soilInit`) is anchored to each system's year-2 +- **fitting window is 2005 to 2011.** the measured soil carbon drops sharply + between 2003 and 2004 as the cropping systems are established. that first + year drop is a one time disturbance, not gradual decomposition, and a first + order soil model should not be pushed to reproduce a change that large in a + single step. so the target's `years` filter starts the fit at 2005. +- **initial soil carbon** (`soilInit`) is anchored to each system's 2005 measured stock, converted from the reported Mg/ha to the model's kg/m2. - **depth:** the observation is a 0-30 cm stock and the harvest reads the whole model soil pool (`TotSoilCarb`); the two are aligned by anchoring `soilInit` to the measured 0-30 cm stock. meta carries the depth window. - **treatments as sites:** each treatment is one site in the multisite run, so `soilInit.` maps to that system's site id in the settings. + +run records live beside the run output, not in the repo. the estimator itself +is unit tested (`tests/testthat/`); the full sipnet run is validated on the +cluster. From 18ceacaec6c0eb3b33a85b761f20b7ec8b71ceb9 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 20/30] add the joint calibration configuration --- examples/2_joint_soc_n2o/config.yml | 103 ++++++++++++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 examples/2_joint_soc_n2o/config.yml diff --git a/examples/2_joint_soc_n2o/config.yml b/examples/2_joint_soc_n2o/config.yml new file mode 100644 index 0000000..869d096 --- /dev/null +++ b/examples/2_joint_soc_n2o/config.yml @@ -0,0 +1,103 @@ +# joint cal/val calibration: salinas SOC fitted as a period mean level plus +# treatment contrasts, modesto N2O contracted to date-matched treatment +# contrasts and held out as validation. +# +# Rscript scripts/010_prepare_observations.R -c examples/2_joint_soc_n2o/config.yml +# Rscript scripts/012_build_target.R -c examples/2_joint_soc_n2o/config.yml +# Rscript scripts/015_build_settings.R -c examples/2_joint_soc_n2o/config.yml +# Rscript scripts/020_build_priors.R -c examples/2_joint_soc_n2o/config.yml +# Rscript scripts/030_calibrate.R -c examples/2_joint_soc_n2o/config.yml +# Rscript scripts/040_plot.R -c examples/2_joint_soc_n2o/config.yml +default: + scc: /projectnb/dietzelab/ccmmf/usr/akash + cache_dir: /projectnb/dietzelab/ccmmf/usr/akash/cal_val_joint/cache_pass5 + soil_pft: [soil, soil_nfixer] + + observations: + dir: cal-val-data + targets: + - variable: SOC_stock_salinas + sites: [salinas_socs] + source_variables: [SOC_stock_Mg_ha] + units: Mg C ha-1 + years: [2005, 2011] + variance: pooled_cv + # raw per date per treatment cells; 012 contracts them to per date contrasts + - variable: N2O_chamber_daily + sites: [modesto_almond_usda] + source_variables: [N2O_flux_g_N_ha_d] + units: g N2O-N ha-1 day-1 + cell_period: date + variance: pooled_cv + + # how the raw targets become the fitted target: one contraction per entry (012) + target: + - type: period_mean + variable: SOC_stock_salinas + control: socs_sys1 + years: [2005, 2011] + new_variable: SOC_periodmean_salinas + - type: contrast + variable: N2O_chamber_daily + treatment: compost + control: no_compost + new_variable: N2O_contrast_modesto + + priors: + post_distns_params: + # maps to SIPNET litterBreakdownRate; read from the soil pft posterior + - turn_over_time + specified: + # BETY prior 1000000095 on som_respiration_rate ("SIPNET default. + # Reference temp = 0C"); the statewide weibull posterior does not cover an + # intensively tilled irrigated vegetable soil and is not this system's prior + som_respiration_rate: + distn: unif + parama: 0.003 + paramb: 0.6 + # 0 to 1 is the support of a fraction; the BETY unif(0.4, 0.6) traces to an + # uncited 2017 placeholder (prior 1000000246) + fracLitterRespired: + distn: unif + parama: 0.0 + paramb: 1.0 + # calibrated initial state, one entry per treatment, prior anchored on that + # treatment's anchor year observation at its measurement error rather than + # pinned to it: the anchor sample is a measurement, not the truth. the forward + # writes one IC per (treatment, particle). + state: + prefix: "soilInit." + variable: SOC_stock_salinas + anchor_year: 2005 + from_unit: Mg ha-1 + to_unit: kg m-2 + + # held out of the likelihood, still simulated and scored + validation_variables: + - N2O_contrast_modesto + + fixed_params: + soil_respiration_Q10: {sipnet: soilRespQ10, value: 2.0} + + forward: + workspace: /projectnb/dietzelab/ccmmf/usr/akash/cal_val_joint/pass5 + prepared_root: /projectnb/dietzelab/ccmmf/usr/akash/cal_val_joint + blocks: examples/2_joint_soc_n2o/blocks.csv + template: examples/2_joint_soc_n2o/template.xml + sa_root: /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821 + binary: /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/tools/sipnet_d2fc7a2 + var_map: + SOC_stock_salinas: + {model_var: TotSoilCarb, from: kg/m2, to: Mg/ha} + N2O_chamber_daily: + {model_var: N2O_flux, from: kg m-2 s-1, to: g ha-1 day-1} + state_pool: soil_organic_carbon_content + + # measured vs modeled treatment effect figure; the curated validation pair + figures: + treatment_effect: + variable: SOC_stock_salinas + treatment: socs_sys2 + control: socs_sys1 + + eki: {n_particles: 50, n_iterations: 3, seed: 556688} From bf22c1f333b0f625a3f1dd1a938de003ebb4d9bf Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 21/30] declare each block's met, ic and events inputs --- examples/2_joint_soc_n2o/blocks.csv | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 examples/2_joint_soc_n2o/blocks.csv diff --git a/examples/2_joint_soc_n2o/blocks.csv b/examples/2_joint_soc_n2o/blocks.csv new file mode 100644 index 0000000..fdaa798 --- /dev/null +++ b/examples/2_joint_soc_n2o/blocks.csv @@ -0,0 +1,11 @@ +block_id,lat,lon,veg_pft,soil_pft,run_start,run_end,met_src,ic_src,events_src +compost,37.6273,-121.0893,temperate.deciduous_almond,soil,2018-01-01,2019-12-31,37.5N_121W/ERA5.1.2018-01-01.2019-12-31.clim,modesto/9cb08ca2174bede7/IC_site_9cb08ca2174bede7_1.nc,modesto/compost/events_ens_01.in +no_compost,37.6273,-121.0893,temperate.deciduous_almond,soil,2018-01-01,2019-12-31,37.5N_121W/ERA5.1.2018-01-01.2019-12-31.clim,modesto/9cb08ca2174bede7/IC_site_9cb08ca2174bede7_1.nc,modesto/no_compost/events_ens_01.in +socs_sys1,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys1/events.in +socs_sys2,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys2/events.in +socs_sys3,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys3/events.in +socs_sys4,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys4/events.in +socs_sys5,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys5/events.in +socs_sys6,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys6/events.in +socs_sys7,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys7/events.in +socs_sys8,36.62,-121.53,annual_crop_row,soil,2005-01-01,2011-12-31,36.5N_121.5W/ERA5.1.2005-01-01.2011-12-31.clim,salinas/9f296becf416ce87/IC_site_9f296becf416ce87_1.nc,salinas_e25/socs_sys8/events.in From 990a3f8f3a0aabd36ad4572e5c09cc3673b88143 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 22/30] commit the run settings template --- examples/2_joint_soc_n2o/template.xml | 128 ++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 examples/2_joint_soc_n2o/template.xml diff --git a/examples/2_joint_soc_n2o/template.xml b/examples/2_joint_soc_n2o/template.xml new file mode 100644 index 0000000..d829487 --- /dev/null +++ b/examples/2_joint_soc_n2o/template.xml @@ -0,0 +1,128 @@ + + + + joint cal/val calibration template; per block run sections, input paths, dates, and output dirs are inserted by scripts/015_build_settings.R + -1 + + + + + temperate.deciduous_almond + + 20 + + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/temperate.deciduous_almond/post.distns.Rdata + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/temperate.deciduous_almond + + + soil + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/soil/post.distns.Rdata + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/soil + + + annual_crop_row + + 20 + + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/annual_crop_row/post.distns.Rdata + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/annual_crop_row + + + annual_crop_alfalfa + + 20 + + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/annual_crop_alfalfa/post.distns.Rdata + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/annual_crop_alfalfa + + + soil_nfixer + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/soil_nfixer/post.distns.Rdata + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/soil_nfixer + + + annual_crop_corn + + 20 + + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/annual_crop_corn/post.distns.Rdata + /projectnb/dietzelab/ccmmf/usr/akash/calval_sa_v3_20260821/pfts/annual_crop_corn + + + + + TotSoilCarb + NEE + HeteroResp + NPP + AbvGrndWood + Qle + SoilMoistFrac + CH4_flux + N2O_flux + N_fixation + + + uniform + + + sampling + + + sampling + + + sampling + + + + + + + 99000000003 + SIPNET + 2.1.0 + FALSE + ./sipnet + + 0 + 1 + 1 + 1 + 1 + + + + localhost + output/out + output/run + qsub -V -cwd -N @NAME@ -j y -o /dev/null -S /bin/bash + Your job-array ([0-9]+).* + qstat -j @JOBID@ || echo DONE + + ./scripts/sge_array_launcher.sh + -t 1-@NJOBS@ + 140 + + module load R/4.4.3 + + output + output/out + output/run + + + + + + + + + + + + + + + + + From 4f3b29da436354270c2016960e8d9e7557fe9563 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 23/30] describe the joint example --- examples/2_joint_soc_n2o/README.md | 47 ++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 examples/2_joint_soc_n2o/README.md diff --git a/examples/2_joint_soc_n2o/README.md b/examples/2_joint_soc_n2o/README.md new file mode 100644 index 0000000..afc7611 --- /dev/null +++ b/examples/2_joint_soc_n2o/README.md @@ -0,0 +1,47 @@ +# Joint cal/val calibration + +One inversion over the cal/val treatments the curated record supports, rather +than a per site pass. The fitted target is the salinas SOC period mean level +plus its treatment contrasts; the modesto N2O treatment contrast is held out of +the likelihood and scored as validation. Physiological parameters are held at +strong priors; the compost events carry an amendment quality +scaling documented beside the prepared inputs in the workspace. + +This directory holds configuration only, the same shape as +`examples/1_salinas_soc`. The code is the generic numbered scripts in +`../../scripts`, and the run artifacts live in the workspace. + +| file | what it is | +|---|---| +| `config.yml` | targets, priors, pinned parameters, PFTs, EKI settings, workspace path | +| `template.xml` | whole-run PEcAn settings; expanded per block by `015_build_settings.R` | +| `blocks.csv` | one row per fitted treatment: dates, PFTs, and the exact met, template IC, and events file each block runs | + +## Running it + +```sh +CFG=examples/2_joint_soc_n2o/config.yml +Rscript scripts/010_prepare_observations.R -c $CFG # raw observation cache +Rscript scripts/012_build_target.R -c $CFG # contract to the fitted target +Rscript scripts/015_build_settings.R -c $CFG # settings + default.param + template ICs +Rscript scripts/020_build_priors.R -c $CFG +Rscript scripts/030_calibrate.R -c $CFG # add --dry-run --particles 3 for a proof pass +Rscript scripts/040_plot.R -c $CFG +``` + +Run from the workspace: the settings carry the array launcher as +`./scripts/sge_array_launcher.sh` and `qsub` runs with `-cwd`, so the working +directory has to be the workspace (`scripts/` sits beside `settings.xml` in +every working run tree). + +## Workspace + +`/projectnb/dietzelab/ccmmf/usr/akash/cal_val_joint/`, on geo. + +``` +template.xml global sections; expanded per block by 015 +pass5/ run workspace: settings.xml, sipnet.default.param, + scripts/, per iteration output +inputs/ events prepared from the curated management record +cache_pass5/ obs.rds, target.rds, prior.rds, result.rds +``` From 54c9670529279d49cd303b522356e074ea104158 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 24/30] document every config key --- examples/README.md | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 examples/README.md diff --git a/examples/README.md b/examples/README.md new file mode 100644 index 0000000..2993b10 --- /dev/null +++ b/examples/README.md @@ -0,0 +1,45 @@ +# Examples + +Each directory is one complete calibration configuration for the shared +numbered scripts in `../scripts`. Start with `1_salinas_soc` (one site, one +variable, a prepared run) and move to `2_joint_soc_n2o` (joint multi-variable +target, settings build, holdout validation). The package code in `R/` carries +nothing site or variable specific; everything a run needs is declared here. + +## Config reference + +Every key the scripts read from a `config.yml`, grouped as in the files. +Optional keys are marked; everything else is required by the script that +consumes it. + +| key | consumed by | meaning | +|---|---|---| +| `scc` | 010, 020, 030, 040 | root that relative data and run paths hang off | +| `cache_dir` | all | directory for `obs.rds`, `target.rds`, `prior.rds`, `result.rds` | +| `soil_pft` | 020, 030 | soil PFT name(s); a vector shares the calibrated rates across all of them | +| `observations.dir` | 010 | cal-val-data checkout, relative to `scc` | +| `observations.targets[]` | 010 | one entry per observed variable: `variable`, `sites`, `units`; optional `source_variables`, `years`, `variance` ("replicate" or "pooled_cv"), `cell_period` ("year" or "date") | +| `target[]` (optional) | 012, 030 | contractions from raw slots to the fitted target, one per entry, dispatched on `type`: `period_mean` (`variable`, `control`, optional `years`, `new_variable`) or `contrast` (`variable`, `treatment`, `control`, `new_variable`). Omit the block to fit the raw slots | +| `priors.post_distns_params` | 020 | traits read from the soil PFT meta-analysis posterior | +| `priors.specified` | 020 | explicit priors, one entry per trait: `distn`, `parama`, `paramb` | +| `priors.state` (optional) | 020, 030 | calibrated initial state: `prefix`, `variable`, `from_unit`, `to_unit`, optional `anchor_year` (default earliest) | +| `fit` (optional) | 030 | sitename -> first observation year to fit; earlier years stay initial condition only | +| `validation_variables` (optional) | 030, 040 | fitted-target variables held out of the likelihood, still predicted and scored | +| `fixed_params` | 015, 030 | parameters pinned in `default.param`, one entry per trait: `sipnet` (model name), `value` | +| `forward.workspace` | 015, 020, 030, 040 | run workspace holding `settings.xml`; `forward.run_dir` (relative to `scc`) is the prepared-run alternative | +| `forward.blocks` | 015 | blocks table: one row per treatment with dates, PFTs, and the exact met, IC, and events file | +| `forward.template` | 015 | whole-run PEcAn settings template, expanded per block | +| `forward.prepared_root` (optional) | 015 | root of prepared inputs and the launcher; defaults to the workspace | +| `forward.sa_root` | 015 | staging tree holding `inputs/met` and `inputs/IC` | +| `forward.binary` | 015 | SIPNET binary the run pins | +| `forward.var_map` | 030 | per observed variable: `model_var`, `from`, `to` -- the crosswalk from model output to observation units | +| `forward.state_pool` | 030 | initial condition pool the calibrated state writes into | +| `figures.treatment_effect` (optional) | 040 | measured vs modeled effect figure: `variable`, `treatment`, `control` | +| `eki` | 015, 030 | `n_particles`, `n_iterations`, `seed` | + +## Script sequence + +`010` observations, `012` fitted target (only when `target` is declared), +`015` settings build (only when the run is not already prepared), `020` prior, +`030` calibration, `040` figures and scores. Each takes +`--config /config.yml`. From 99f83d82fde28ca497b15cc010e4c0904ea8b1a3 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:01:59 -0400 Subject: [PATCH 25/30] report the pass 5 run --- examples/2_joint_soc_n2o/report.qmd | 233 ++++++++++++++++++++++++++++ 1 file changed, 233 insertions(+) create mode 100644 examples/2_joint_soc_n2o/report.qmd diff --git a/examples/2_joint_soc_n2o/report.qmd b/examples/2_joint_soc_n2o/report.qmd new file mode 100644 index 0000000..7932e84 --- /dev/null +++ b/examples/2_joint_soc_n2o/report.qmd @@ -0,0 +1,233 @@ +--- +title: "Joint SOC and N2O calibration: run report" +date: today +format: + html: + self-contained: true + toc: true + theme: cosmo + number-sections: true +execute: + echo: false + warning: false + message: false +--- + +```{r} +#| label: setup +#| include: false +library(calibration) +library(dplyr) +library(ggplot2) +library(knitr) + +config <- config::get(file = file.path(here::here(), "examples/2_joint_soc_n2o/config.yml")) +cache <- config$cache_dir +result <- readRDS(file.path(cache, "result.rds")) +target <- readRDS(file.path(cache, "target.rds")) +obs <- readRDS(file.path(cache, "obs.rds")) +prior <- readRDS(file.path(cache, "prior.rds")) + +G_list <- c(lapply(result$eki$eki_list, function(s) s$G), list(result$eki$G)) +names(G_list) <- paste("iteration", seq_along(G_list)) +shift <- param_shift(result$U0, result$U) +scores <- score_table(G_list, result$obs$meta) + +theme_set(theme_minimal(base_size = 12) + + theme(panel.grid.minor = element_blank(), legend.position = "bottom")) +``` + +This report interprets one executed calibration of SIPNET soil parameters +against the curated cal/val observations, run with the configuration in this +directory. The run itself is reproduced by the script sequence in the README; +nothing here reruns the model: every figure and number is read from the run's +cached artifacts under `` `r cache` ``. + +# What was calibrated + +The fitted target is soil organic carbon at the Salinas organic cropping +systems trial, contracted from `r sum(obs$meta$variable == "SOC_stock_salinas")` +per system per year stocks (8 systems, 2005 to 2011) into +`r sum(result$obs$meta$variable == "SOC_periodmean_salinas")` slots: the control +system's period mean level and one period mean contrast per remaining system. +The Modesto almond N2O treatment contrast +(`r sum(target$meta$variable == "N2O_contrast_modesto")` date-matched slots) is +held out of the likelihood and scored as validation. + +Three soil parameters are estimated (`som_respiration_rate`, +`fracLitterRespired`, `turn_over_time`) together with one initial soil carbon +stock per system (`soilInit.*`), each anchored on that system's 2005 observation +at its measurement error rather than pinned to it. `soil_respiration_Q10` is +pinned at 2.0 in the run's `default.param`. The ensemble runs +`r result$control$n_particles` particles through +`r result$control$n_iterations` tempered ensemble Kalman updates +(seed `r result$control$seed`). + +# Does the calibrated model reproduce the observed target? + +```{r} +#| label: fitted-slots +#| fig-cap: "Fitted target slots before (left) and after (right) calibration: the control period mean level and the seven system contrasts. Blue points are the ensemble members, the black dash their mean, red the observation with one standard deviation." +#| fig-width: 10 +#| fig-height: 4.5 +sm <- result$obs$meta +library(patchwork) +p1 <- plot_target_slots(G_list[[1]][, sm$slot, drop = FALSE], sm, + ylab = unique(sm$units), title = "before calibration") +p2 <- plot_target_slots(G_list[[length(G_list)]][, sm$slot, drop = FALSE], sm, + ylab = unique(sm$units), title = "after calibration") +p1 + p2 +``` + +```{r} +#| label: scores +soc <- scores[scores$variable == "SOC_periodmean_salinas", ] +kable(soc[, c("iteration", "rmse", "bias", "coverage", "crps")], digits = 2, + caption = "Fit of the SOC target by tempering step. Iteration 1 is the prior forward run; the last row is the posterior.") +``` + +The prior ensemble misses the target with an RMSE of +`r round(soc$rmse[1], 1)` `r unique(sm$units)`; the calibrated ensemble closes +that to `r round(soc$rmse[nrow(soc)], 1)` with coverage +`r round(soc$coverage[nrow(soc)], 2)` of observations inside the 90 percent +band. + +# What did the data constrain? + +```{r} +#| label: shift-table +kable(shift[, c("param", "prior_mean", "post_mean", "post_sd", "var_reduction")], + digits = 3, + caption = "Prior to posterior shift per parameter. var_reduction is the fractional loss of prior variance; a value near zero means the data did not inform the parameter.") +``` + +```{r} +#| label: som-density +#| fig-cap: "Prior and posterior for som_respiration_rate, the one parameter the target strongly constrains." +#| fig-width: 6 +#| fig-height: 3 +d <- rbind( + data.frame(stage = "prior", value = result$U0[, "som_respiration_rate"]), + data.frame(stage = "posterior", value = result$U[, "som_respiration_rate"]) +) +ggplot(d, aes(value, color = stage)) + + geom_density(linewidth = 1) + + labs(x = "som_respiration_rate (g C g-1 soil C yr-1)", y = "density", color = NULL) +``` + +```{r} +#| label: som-numbers +#| include: false +sr <- result$U[, "som_respiration_rate"] +srq <- quantile(sr, c(0.025, 0.975)) +vr <- shift$var_reduction[shift$param == "som_respiration_rate"] +``` + +The target constrains one process rate: `som_respiration_rate` collapses from +its uniform prior to `r round(mean(sr), 3)` +[`r round(srq[1], 3)`, `r round(srq[2], 3)`], a +`r round(100 * vr)` percent variance reduction, and it does so in the first +tempering step (below). `turn_over_time` and `fracLitterRespired` barely +narrow: the period mean SOC target carries little information about how carbon is +partitioned on its way to the soil pool, only about how fast the pool turns +over. The initial stocks stay at their observation-anchored priors, which is +what the anchoring is for. + +```{r} +#| label: trace +#| fig-cap: "Ensemble mean and 5-95 percent range of each parameter across tempering steps; dashed lines are the declared support bounds. Initial stocks are in kg C m-2, in model units." +#| fig-width: 10 +#| fig-height: 6 +pm <- result$eki$par_map +trace <- c(list(result$U0), lapply(result$eki$eki_list, function(s) pm$inv(s$U))) +plot_param_trace(trace, prior) +``` + +# Does the model capture the compost effect? + +The curated treatment pair table names systems 1 and 2 as the compost +comparison: the same quadrennial cover crop with and without 7.6 Mg/ha compost +per crop. The effect is computed per ensemble member and year as treatment +minus control, from the raw per year stocks of the prior (itr1) and posterior +(itr4) forward runs. + +```{r} +#| label: effect-harvest +#| include: false +em <- obs$meta[obs$meta$variable == "SOC_stock_salinas" & + obs$meta$treatment_id %in% c("socs_sys1", "socs_sys2"), ] +vm <- config$forward$var_map["SOC_stock_salinas"] +win <- cbind(socs_sys1 = c(2005L, 2011L), socs_sys2 = c(2005L, 2011L)) +out <- file.path(config$forward$workspace, "output") +G_prior <- harvest_output_to_G(file.path(out, "itr1", "out"), em, vm, win) +G_post <- harvest_output_to_G(file.path(out, "itr4", "out"), em, vm, win) +eff <- function(G) mean(colMeans(G[, em$slot[em$treatment_id == "socs_sys2"]]) - + colMeans(G[, em$slot[em$treatment_id == "socs_sys1"]])) +meas <- with(em, mean(value[treatment_id == "socs_sys2"] - + value[treatment_id == "socs_sys1"])) +``` + +```{r} +#| label: effect-figure +#| fig-cap: "Measured against modeled compost effect (system 2 minus system 1), before and after calibration, as the absolute effect and as percent of control. Bands are the 5-95 percent ensemble range; red points the measured contrast with one standard deviation." +#| fig-width: 9 +#| fig-height: 4 +plot_treatment_effect(G_prior, G_post, em, "SOC_stock_salinas", + treatment = "socs_sys2", control = "socs_sys1") +``` + +Averaged over the period, the measured effect is +`r round(meas, 1)` Mg C/ha; the prior ensemble puts it at +`r round(eff(G_prior), 1)` and the calibrated ensemble at +`r round(eff(G_post), 1)`. The calibration was fit to the contrasts, so +matching their mean is expected; that the posterior band also tracks the +year to year shape in both the absolute and relative form is the substantive +check. + +# The held-out N2O contrast + +```{r} +#| label: validation +#| fig-cap: "Held-out validation: predicted against measured Modesto N2O treatment contrast per measurement date, with the 5-95 percent ensemble range and a 1:1 line. Axes are linear because a contrast is signed." +#| fig-width: 8 +#| fig-height: 4 +vmeta <- result$obs_all$meta +vmeta <- vmeta[vmeta$variable == "N2O_contrast_modesto", , drop = FALSE] +plot_validation(result$G_validation, vmeta) +``` + +```{r} +#| label: val-scores +kable(score_iteration(result$G_validation, vmeta)[, c("variable", "slots", "rmse", "bias", "coverage")], + digits = 2, caption = "Scores of the held-out contrast at the posterior.") +``` + +The model predicts a positive compost minus control N2O contrast on most dates +while several measured contrasts are negative. That sign structure is why the +variable is held out: a target the model cannot reach cannot be calibrated +through; the estimator would answer by pushing a nitrogen parameter to its +bound for reasons unrelated to the process it stands for. Held out, the +mismatch stays visible and scored instead of being absorbed. + +# Limitations + +- The compost carbon inputs carry an amendment quality scaling sourced from the + literature, not fitted; its uncertainty is not propagated here. +- One process rate dominates the fit. With a single site's period mean and + contrasts, `turn_over_time` and `fracLitterRespired` are close to + unidentifiable, and their posteriors should be read as priors. +- The SOC inference rests on one trial; transfer to other soils and systems is + an assumption until a second site joins the target. +- The establishment years (2003 to 2004) are excluded by design: the first year + drop is a one time disturbance a first order soil model should not be forced + to reproduce. +- The N2O response has the wrong sign structure at Modesto; nitrogen parameters + were not calibrated and the held-out scores quantify, not fix, that gap. + +# Reproducing this report + +The run is reproduced by the script sequence in this directory's README; this +document then renders from the run's caches with +`quarto render examples/2_joint_soc_n2o/report.qmd`. The only computation here +beyond reading caches is the compost effect harvest, which reads the stored +model output of the first and last iteration. From db5f43ab43bfef2c833afc72bb91d51480e48026 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:02:00 -0400 Subject: [PATCH 26/30] ignore rendered reports --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index f483eee..44facfe 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,7 @@ .Rhistory .RData .Rproj.user/ + +# rendered reports, regenerated from the qmd +examples/*/report.html +examples/*/report_files/ From 92ebaac4940fd526d57af855a66d2c5fb7669404 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:02:00 -0400 Subject: [PATCH 27/30] update the vignette for the multi target api --- vignettes/calibration_demo.qmd | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/vignettes/calibration_demo.qmd b/vignettes/calibration_demo.qmd index 8c0e3e6..f70f282 100644 --- a/vignettes/calibration_demo.qmd +++ b/vignettes/calibration_demo.qmd @@ -147,15 +147,17 @@ the general package. ## The observations -`build_obs()` reads the curated cal/val data and builds the target for one -variable at the sites you name. Each site, treatment, and year becomes one slot; -the value is the replicate mean and the variance comes from the replicate spread. +`build_obs()` reads the curated cal/val data and stacks a target from a list of +`targets`, one per variable with the sites that carry it. Each variable, site, +treatment, and year becomes one slot; the value is the replicate mean and the +variance comes from the replicate spread, or the reported standard error where a +cell has no replicates. A single-variable target is a one-entry list; a joint multi-variable +target adds more entries. ```{r sipnet-obs, eval = FALSE} obs <- build_obs( cal_val_dir = "/projectnb/dietzelab/ccmmf/usr/akash/cal_val", - target_var = "SOC_stock_Mg_ha", - sites = "salinas_socs" + targets = list(list(variable = "SOC_stock_Mg_ha", sites = "salinas_socs")) ) ``` @@ -169,13 +171,14 @@ Priors come from three sources, combined into one list: ```{r sipnet-prior, eval = FALSE} prior <- c( - prior_from_postdistns("som_respiration_rate", soil_pft$posterior.files), + prior_from_postdistns("som_respiration_rate", "pfts/soil/post.distns.Rdata"), prior_from_specs(list( soil_respiration_Q10 = list(distn = "unif", parama = 1.4, paramb = 3.0), turn_over_time = list(distn = "unif", parama = 0.13, paramb = 1.2) )), state_prior_from_obs(obs$meta, prefix = "soilInit.", - from_unit = "Mg/ha", to_unit = "kg/m2", anchor_year = 2) + variable = "SOC_stock_Mg_ha", + from_unit = "Mg/ha", to_unit = "kg/m2") ) ``` @@ -183,8 +186,8 @@ prior <- c( `make_forward_sipnet()` wraps a prepared PEcAn multisite run as the `forward` function. It runs the ensemble through SIPNET, using the settings and the SGE -launcher exactly as prepared, and harvests the variable back to the observation -units. +launcher exactly as prepared, and harvests each observed variable back to its +observation units through the `var_map` crosswalk. ```{r sipnet-forward, eval = FALSE} settings <- PEcAn.settings::read.settings( @@ -193,8 +196,10 @@ settings <- PEcAn.settings::read.settings( forward <- make_forward_sipnet( settings = settings, obs = obs, n_particles = 50, - harvest_var = "TotSoilCarb", from_unit = "kg/m2", to_unit = "Mg/ha", - soil_pft = "soil", state_prefix = "soilInit.", + var_map = list( + SOC_stock_Mg_ha = list(model_var = "TotSoilCarb", from = "kg/m2", to = "Mg/ha") + ), + soil_pfts = "soil", state_prefix = "soilInit.", state_pool = "soil_organic_carbon_content" ) ``` @@ -203,7 +208,7 @@ forward <- make_forward_sipnet( This launches the SIPNET ensemble on the cluster, so it is not run here. In practice you run it through the numbered scripts (`scripts/010` to `040`) against -a config; see `examples/salinas_soc/`. +a config; see `examples/1_salinas_soc/`. ```{r sipnet-run, eval = FALSE} result <- calibrate( @@ -241,9 +246,9 @@ prior to posterior shift in each parameter. To calibrate something else, you change the four pieces, not the estimator: -- point `build_obs()` at a different `target_var` and `sites`, +- give `build_obs()` the `targets` you want (variables and their sites), - build a prior for the parameters you want, -- give `make_forward_sipnet()` (or your own `forward`) the variable and units, +- give `make_forward_sipnet()` (or your own `forward`) a `var_map` for those variables, - keep `calibrate()` as is. The estimator only ever sees `y`, `Sigma`, and the model predictions, so it does From 6991b2e23fb537276ad2b351a53c3b74bcd9ff95 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:02:00 -0400 Subject: [PATCH 28/30] point the readme at the numbered example --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 16a1175..da95bc6 100644 --- a/README.md +++ b/README.md @@ -75,7 +75,7 @@ Rscript scripts/030_calibrate.R --config Rscript scripts/040_plot.R --config ``` -`examples/salinas_soc/` is a worked example (soil carbon at the Salinas organic +`examples/1_salinas_soc/` is a worked example (soil carbon at the Salinas organic cropping systems), and `vignettes/calibration_demo.qmd` walks through the whole thing step by step. From 6e3fcfff37925216d11d39a0f09de38b430abc78 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:02:00 -0400 Subject: [PATCH 29/30] keep run debris out of the package build --- .Rbuildignore | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.Rbuildignore b/.Rbuildignore index 5dc7b01..c0065ee 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -6,3 +6,8 @@ ^\.git$ ^examples$ ^vignettes$ +^analysis$ +^cache_ +^figures_ +^harmonize$ +^cal_val_manual$ From 173f5a851a8e793b0b6bb029170025afe738e931 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Fri, 4 Sep 2026 12:02:00 -0400 Subject: [PATCH 30/30] regenerate the documentation --- NAMESPACE | 11 ++++++++ man/apply_transform.Rd | 24 ++++++++++++++++ man/baseline_trait_samples.Rd | 20 ++++++------- man/bind_obs.Rd | 23 +++++++++++++++ man/build_obs.Rd | 29 +++++++++---------- man/cell_variance.Rd | 16 +++++++++++ man/contrast_target.Rd | 36 +++++++++++++++++++++++ man/dot-distn_support.Rd | 17 +++++++---- man/dot-shrink_to_pd.Rd | 22 +++++++++++++++ man/harvest_output_to_G.Rd | 30 ++++++++------------ man/inject_traits.Rd | 10 +++++-- man/label_variable.Rd | 12 ++++++++ man/make_forward_sipnet.Rd | 29 +++++++++++-------- man/param_label.Rd | 14 +++++++++ man/param_label_unit.Rd | 12 ++++++++ man/param_unit.Rd | 16 +++++++++++ man/period_mean_contrast.Rd | 44 +++++++++++++++++++++++++++++ man/plot_param_trace.Rd | 25 ++++++++++++++++ man/plot_target_slots.Rd | 28 ++++++++++++++++++ man/plot_treatment_effect.Rd | 30 ++++++++++++++++++++ man/plot_validation.Rd | 26 +++++++++++++++++ man/prior_from_shared_postdistns.Rd | 25 ++++++++++++++++ man/read_cal_val_observations.Rd | 12 ++------ man/repoint_poolinitcond.Rd | 12 ++++++-- man/run_window.Rd | 22 +++++++++++++++ man/score_iteration.Rd | 10 ++++--- man/state_prior_from_obs.Rd | 16 +++++++++-- man/subset_obs.Rd | 24 ++++++++++++++++ man/target_rows.Rd | 12 ++++++++ man/write_samples_rdata.Rd | 12 -------- 30 files changed, 521 insertions(+), 98 deletions(-) create mode 100644 man/apply_transform.Rd create mode 100644 man/bind_obs.Rd create mode 100644 man/cell_variance.Rd create mode 100644 man/contrast_target.Rd create mode 100644 man/dot-shrink_to_pd.Rd create mode 100644 man/label_variable.Rd create mode 100644 man/param_label.Rd create mode 100644 man/param_label_unit.Rd create mode 100644 man/param_unit.Rd create mode 100644 man/period_mean_contrast.Rd create mode 100644 man/plot_param_trace.Rd create mode 100644 man/plot_target_slots.Rd create mode 100644 man/plot_treatment_effect.Rd create mode 100644 man/plot_validation.Rd create mode 100644 man/prior_from_shared_postdistns.Rd create mode 100644 man/run_window.Rd create mode 100644 man/subset_obs.Rd create mode 100644 man/target_rows.Rd delete mode 100644 man/write_samples_rdata.Rd diff --git a/NAMESPACE b/NAMESPACE index 1dfedc7..76a63f5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,19 +1,30 @@ # Generated by roxygen2: do not edit by hand +export(apply_transform) +export(bind_obs) export(build_obs) export(calibrate) export(calibration_control) +export(contrast_target) export(get_par_map_funcs) export(harvest_output_to_G) export(make_forward_sipnet) export(param_shift) +export(period_mean_contrast) export(plot_ensembles_vs_truth) export(plot_param_densities) +export(plot_param_trace) +export(plot_target_slots) +export(plot_treatment_effect) +export(plot_validation) export(prior_from_postdistns) +export(prior_from_shared_postdistns) export(prior_from_specs) export(run_eki) +export(run_window) export(sample_initial_ensemble) export(save_iterations_pdf) export(score_iteration) export(score_table) export(state_prior_from_obs) +export(subset_obs) diff --git a/man/apply_transform.Rd b/man/apply_transform.Rd new file mode 100644 index 0000000..cf0338d --- /dev/null +++ b/man/apply_transform.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observation_operator.R +\name{apply_transform} +\alias{apply_transform} +\title{Apply a target transform to a prediction matrix} +\usage{ +apply_transform(G, transform) +} +\arguments{ +\item{G}{prediction matrix (members x raw slots).} + +\item{transform}{matrix mapping raw slots (columns) to fitted slots (rows).} +} +\value{ +matrix (members x fitted slots), columns named by fitted slot. +} +\description{ +Puts model predictions on the raw slots through the same linear + map the observations went through, so the fitted quantity is the same + operation on both sides. +} +\author{ +Akash BV +} diff --git a/man/baseline_trait_samples.Rd b/man/baseline_trait_samples.Rd index e0428d9..124efb5 100644 --- a/man/baseline_trait_samples.Rd +++ b/man/baseline_trait_samples.Rd @@ -2,19 +2,17 @@ % Please edit documentation in R/forward_sipnet.R \name{baseline_trait_samples} \alias{baseline_trait_samples} -\title{fixed baseline trait samples: every parameter at its prior median (the -PEcAn.priors::get.sample p = 0.5 of the post.distns row, so it matches the -family the pft carries), replicated over particles, one data.frame per pft. -the calibrated columns are overwritten by U; the rest stay fixed so the -prediction spread reflects only the estimated parameters.} +\title{baseline trait samples: every parameter at its prior median, replicated over +particles, one data.frame per pft. calibrated columns are overwritten by U. +traits named in `fixed` are dropped so the run dir default.param value stands; +a trait left in here overwrites it.} \usage{ -baseline_trait_samples(pfts, n_particles) +baseline_trait_samples(pfts, n_particles, fixed = character(0)) } \description{ -fixed baseline trait samples: every parameter at its prior median (the -PEcAn.priors::get.sample p = 0.5 of the post.distns row, so it matches the -family the pft carries), replicated over particles, one data.frame per pft. -the calibrated columns are overwritten by U; the rest stay fixed so the -prediction spread reflects only the estimated parameters. +baseline trait samples: every parameter at its prior median, replicated over +particles, one data.frame per pft. calibrated columns are overwritten by U. +traits named in `fixed` are dropped so the run dir default.param value stands; +a trait left in here overwrites it. } \keyword{internal} diff --git a/man/bind_obs.Rd b/man/bind_obs.Rd new file mode 100644 index 0000000..42dd73a --- /dev/null +++ b/man/bind_obs.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observation_operator.R +\name{bind_obs} +\alias{bind_obs} +\title{Combine observation targets into one} +\usage{ +bind_obs(...) +} +\arguments{ +\item{...}{obs lists, each list(y, Sigma, meta, transform).} +} +\value{ +a single obs list with the transforms combined. +} +\description{ +Stacks several obs lists into a single target with a + block-diagonal covariance. Targets from different variables are assumed + independent of each other, which is what block diagonal encodes; + correlations within a target are carried through from the input blocks. +} +\author{ +Akash BV +} diff --git a/man/build_obs.Rd b/man/build_obs.Rd index c0ee68c..2a46a58 100644 --- a/man/build_obs.Rd +++ b/man/build_obs.Rd @@ -2,31 +2,28 @@ % Please edit documentation in R/observations.R \name{build_obs} \alias{build_obs} -\title{Build the calibration target from cal/val observations} +\title{Build the calibration target from the curated cal/val observations} \usage{ -build_obs(cal_val_dir, target_var, sites, rel_var_floor = 0.05) +build_obs(cal_val_dir, targets) } \arguments{ -\item{cal_val_dir}{directory of the cal/val tsv export.} +\item{cal_val_dir}{the cal-val-data checkout root.} -\item{target_var}{the cal/val variable to calibrate to.} - -\item{sites}{character vector of sitenames to include.} - -\item{rel_var_floor}{relative floor on each cell variance: the variance is at -least `(rel_var_floor * value)^2`, so an agreeing cell cannot make Sigma -singular. Set 0 to disable.} +\item{targets}{list of target specs: `variable`, `sites`, `units` (the unit +every source record must carry), and optionally `source_variables` (raw names +feeding it, default `variable`), `years` (inclusive c(first, last) filter), +`variance` ("replicate", the default, or "pooled_cv"), and `cell_period` +("year", the default, or "date" for an episodic sub-annual flux).} } \value{ -list(y, Sigma, meta): `y` named numeric (length P), `Sigma` a P x P - diagonal variance matrix named to match `y`, `meta` one row per slot. +list(y, Sigma, meta). } \description{ Assembles the observation vector `y`, its diagonal likelihood -covariance `Sigma`, and the per-slot `meta` for one variable at one or more -sites. Each (site, treatment, year) cell becomes one slot: `y` is the -replicate mean, `Sigma` the variance of that mean floored relative to its -magnitude. Nothing here is variable- or site specific; the caller names them. + covariance `Sigma`, and the per-slot `meta`, stacking every target into one + vector. A cell is one (variable, site, treatment, year, depth) group; the slot + name carries all of them so the estimator, which aligns by name only, never has + to know what a site or a variable is. } \author{ Akash BV diff --git a/man/cell_variance.Rd b/man/cell_variance.Rd new file mode 100644 index 0000000..a2630aa --- /dev/null +++ b/man/cell_variance.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observations.R +\name{cell_variance} +\alias{cell_variance} +\title{variance of each cell mean: replicate spread where the cell has replicates, +otherwise the reported standard error. no relative floor -- it would override +a reported error by orders of magnitude.} +\usage{ +cell_variance(cells) +} +\description{ +variance of each cell mean: replicate spread where the cell has replicates, +otherwise the reported standard error. no relative floor -- it would override +a reported error by orders of magnitude. +} +\keyword{internal} diff --git a/man/contrast_target.Rd b/man/contrast_target.Rd new file mode 100644 index 0000000..5b12800 --- /dev/null +++ b/man/contrast_target.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observation_operator.R +\name{contrast_target} +\alias{contrast_target} +\title{Contract per treatment observations into a treatment contrast} +\usage{ +contrast_target( + obs, + variable, + treatment, + control, + new_variable = paste0(variable, "_contrast") +) +} +\arguments{ +\item{obs}{a build_obs target list(y, Sigma, meta).} + +\item{variable}{the per treatment variable to contract.} + +\item{treatment, control}{treatment ids to compare.} + +\item{new_variable}{name for the resulting variable.} +} +\value{ +an obs list carrying one slot per shared date, with the contraction + recorded in `transform`. +} +\description{ +One slot per date, the difference between a treatment and its + control on the dates both were measured. Chamber flux data supports + treatment comparisons rather than absolute magnitudes, which is what this + contraction fits. +} +\author{ +Akash BV +} diff --git a/man/dot-distn_support.Rd b/man/dot-distn_support.Rd index 4315797..a786b6d 100644 --- a/man/dot-distn_support.Rd +++ b/man/dot-distn_support.Rd @@ -3,14 +3,19 @@ \name{.distn_support} \alias{.distn_support} \title{distn family -> support the transport map uses. positive-support families get -a log map, beta a logit, normal an identity, and a uniform its own bounds -(handled by prior_from_specs, not here).} +a log map, beta and unif a logit onto their own bounds, normal an identity.} \usage{ -.distn_support(distn) +.distn_support(distn, parama = NULL, paramb = NULL) +} +\arguments{ +\item{distn}{family name.} + +\item{parama, paramb}{distribution parameters; uniform takes its support from them.} } \description{ -distn family -> support the transport map uses. positive-support families get -a log map, beta a logit, normal an identity, and a uniform its own bounds -(handled by prior_from_specs, not here). +there is deliberately no default branch: a family falling through to +c(-Inf, Inf) estimates a bounded parameter on the whole real line, and the +posterior can leave the physical range without anything objecting. a new +family must state its support rather than inherit an unbounded one. } \keyword{internal} diff --git a/man/dot-shrink_to_pd.Rd b/man/dot-shrink_to_pd.Rd new file mode 100644 index 0000000..dfc2d43 --- /dev/null +++ b/man/dot-shrink_to_pd.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observation_operator.R +\name{.shrink_to_pd} +\alias{.shrink_to_pd} +\title{shrink correlations toward zero by the least amount that restores positive +definiteness, leaving the diagonal exactly as estimated. a mean over K years +gives an empirical covariance of rank at most K - 1, so more quantities than +years is singular by construction and the EnKF's Cholesky of cov(G) + Sigma +has no reason to succeed. the marginal variances are well estimated from K +years; the correlations are not, so only they are damped.} +\usage{ +.shrink_to_pd(S, label = "", tol = 1e-08) +} +\description{ +shrink correlations toward zero by the least amount that restores positive +definiteness, leaving the diagonal exactly as estimated. a mean over K years +gives an empirical covariance of rank at most K - 1, so more quantities than +years is singular by construction and the EnKF's Cholesky of cov(G) + Sigma +has no reason to succeed. the marginal variances are well estimated from K +years; the correlations are not, so only they are damped. +} +\keyword{internal} diff --git a/man/harvest_output_to_G.Rd b/man/harvest_output_to_G.Rd index 4db310c..2b02e9a 100644 --- a/man/harvest_output_to_G.Rd +++ b/man/harvest_output_to_G.Rd @@ -4,36 +4,28 @@ \alias{harvest_output_to_G} \title{Harvest a model ensemble into the prediction matrix G} \usage{ -harvest_output_to_G( - out_root, - meta, - variable, - start_year, - end_year, - from_unit, - to_unit -) +harvest_output_to_G(out_root, meta, var_map, run_window) } \arguments{ \item{out_root}{the model output directory holding the ENS-* run dirs.} -\item{meta}{observation meta (slot, treatment_id, min_date, max_date).} +\item{meta}{observation meta (slot, treatment_id, variable, min_date, max_date).} -\item{variable}{the model output variable to read.} +\item{var_map}{named list keyed by observation `variable`, each +`list(model_var, from, to)`: the model output to read, its unit, and the +observation unit.} -\item{start_year, end_year}{the run year range.} - -\item{from_unit}{the variable's model unit (udunits string).} - -\item{to_unit}{the observation unit to convert to (udunits string).} +\item{run_window}{integer matrix (2 x n_treatments) of first and last run year, +columns named by treatment. Each treatment is read over its own window; a +joint run spans different periods per site.} } \value{ matrix (members x slots) named by observation slot, member-ordered. } \description{ -For each ENS-- run under `out_root`, reads -`variable`, samples it at the observation date of each of that treatment's -slots (the midpoint of the slot's date window in meta), converts to the +For each ENS-- run under `out_root`, reads the +model outputs the treatment's slots need over that treatment's run window, +samples each slot at the midpoint of its date window, converts to the observation unit, and assembles the J x P matrix aligned to the observation slots. Assumes every expected run has finished; a missing run output fails loud in read.output rather than being silently dropped. diff --git a/man/inject_traits.Rd b/man/inject_traits.Rd index d50ef59..37e2d91 100644 --- a/man/inject_traits.Rd +++ b/man/inject_traits.Rd @@ -2,11 +2,15 @@ % Please edit documentation in R/forward_sipnet.R \name{inject_traits} \alias{inject_traits} -\title{overwrite the pft's calibrated trait columns with the proposal U.} +\title{write the proposal U into each named soil pft: the calibrated rates are one +shared quantity, not one per pft. fails if a named pft is absent rather than +silently calibrating a subset.} \usage{ -inject_traits(baseline, soil_pft, U_traits) +inject_traits(baseline, soil_pfts, U_traits) } \description{ -overwrite the pft's calibrated trait columns with the proposal U. +write the proposal U into each named soil pft: the calibrated rates are one +shared quantity, not one per pft. fails if a named pft is absent rather than +silently calibrating a subset. } \keyword{internal} diff --git a/man/label_variable.Rd b/man/label_variable.Rd new file mode 100644 index 0000000..882d741 --- /dev/null +++ b/man/label_variable.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{label_variable} +\alias{label_variable} +\title{readable variable label: the variable name with underscores as spaces.} +\usage{ +label_variable(x) +} +\description{ +readable variable label: the variable name with underscores as spaces. +} +\keyword{internal} diff --git a/man/make_forward_sipnet.Rd b/man/make_forward_sipnet.Rd index c06a1af..7099b01 100644 --- a/man/make_forward_sipnet.Rd +++ b/man/make_forward_sipnet.Rd @@ -8,13 +8,13 @@ make_forward_sipnet( settings, obs, n_particles, - harvest_var, - from_unit, - to_unit, - soil_pft = "soil", + var_map, + soil_pfts, state_prefix = "soilInit.", state_pool = "soil_organic_carbon_content", - base_out_dir = settings$outdir + base_out_dir = settings$outdir, + fixed_traits = character(0), + raw_obs = NULL ) } \arguments{ @@ -24,13 +24,12 @@ make_forward_sipnet( \item{n_particles}{ensemble size J.} -\item{harvest_var}{the model output variable to compare to the observations.} +\item{var_map}{named list keyed by observation variable, each +`list(model_var, from, to)` (see harvest_output_to_G): the crosswalk from +each observed variable to its model output and units.} -\item{from_unit}{the model unit of harvest_var (udunits string).} - -\item{to_unit}{the observation unit to convert the harvest to.} - -\item{soil_pft}{name of the PFT whose traits U overwrites.} +\item{soil_pfts}{character vector of soil PFT names that share the calibrated +rates; the same proposal column is written into each (see inject_traits).} \item{state_prefix}{column prefix marking calibrated initial-state entries in U; when present each is written into the initial condition per particle.} @@ -38,6 +37,14 @@ U; when present each is written into the initial condition per particle.} \item{state_pool}{the initial condition pool variable the state writes into.} \item{base_out_dir}{parent directory for the per iteration model output.} + +\item{fixed_traits}{trait names pinned through the run dir `default.param`, +dropped from the baseline sample so the pinned value is not overwritten by +the PFT posterior median.} + +\item{raw_obs}{the untransformed target the model output is harvested against; +its `transform` (the linear map from raw to fitted slots) is applied to G so +model and observations are the same quantity. NULL fits the raw slots.} } \value{ function(U, iteration) -> matrix (J, P) aligned to names(obs$y). diff --git a/man/param_label.Rd b/man/param_label.Rd new file mode 100644 index 0000000..ed3571f --- /dev/null +++ b/man/param_label.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{param_label} +\alias{param_label} +\title{axis and strip labels for calibrated parameters; unmapped names fall back to +the trait name with underscores as spaces.} +\usage{ +param_label(x) +} +\description{ +axis and strip labels for calibrated parameters; unmapped names fall back to +the trait name with underscores as spaces. +} +\keyword{internal} diff --git a/man/param_label_unit.Rd b/man/param_label_unit.Rd new file mode 100644 index 0000000..7eba394 --- /dev/null +++ b/man/param_label_unit.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{param_label_unit} +\alias{param_label_unit} +\title{parameter label with its unit appended, for a facet strip or page label.} +\usage{ +param_label_unit(x) +} +\description{ +parameter label with its unit appended, for a facet strip or page label. +} +\keyword{internal} diff --git a/man/param_unit.Rd b/man/param_unit.Rd new file mode 100644 index 0000000..d1a1824 --- /dev/null +++ b/man/param_unit.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{param_unit} +\alias{param_unit} +\title{units as SIPNET reads them, not always as documented: baseSoilResp is read in +per year and divided by 365 at setup in sipnet.c, so a calibrated value is a +per year rate; wrong on an axis is a factor of 365.} +\usage{ +param_unit(x) +} +\description{ +units as SIPNET reads them, not always as documented: baseSoilResp is read in +per year and divided by 365 at setup in sipnet.c, so a calibrated value is a +per year rate; wrong on an axis is a factor of 365. +} +\keyword{internal} diff --git a/man/period_mean_contrast.Rd b/man/period_mean_contrast.Rd new file mode 100644 index 0000000..e7094fe --- /dev/null +++ b/man/period_mean_contrast.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observation_operator.R +\name{period_mean_contrast} +\alias{period_mean_contrast} +\title{Period mean level and treatment contrasts} +\usage{ +period_mean_contrast( + obs, + variable, + control, + years = NULL, + new_variable = paste0(variable, "_periodmean") +) +} +\arguments{ +\item{obs}{a build_obs target list(y, Sigma, meta).} + +\item{variable}{the per treatment per year variable to collapse.} + +\item{control}{treatment id every contrast is taken against.} + +\item{years}{optional integer vector restricting the period; defaults to every +year present.} + +\item{new_variable}{name for the resulting variable.} +} +\value{ +an obs list carrying one level slot and one contrast slot per + treatment, with the contraction recorded in `transform`. +} +\description{ +Collapses a per treatment per year series into one level slot for + the control and one contrast slot per remaining treatment, all on the mean + over the period. A bijection of the per treatment means: the control level + carries the net rate, the contrasts carry the treatment effects, and one + level slot keeps a calibrated initial state out of the contrasts. + + The covariance is the empirical covariance of the annual series divided by + the number of years (the standard error of the mean), which carries the + shared-control structure of the contrasts. +} +\author{ +Akash BV +} diff --git a/man/plot_param_trace.Rd b/man/plot_param_trace.Rd new file mode 100644 index 0000000..600fc82 --- /dev/null +++ b/man/plot_param_trace.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{plot_param_trace} +\alias{plot_param_trace} +\title{Parameter trace across tempering steps} +\usage{ +plot_param_trace(trace, prior = NULL) +} +\arguments{ +\item{trace}{named list of parameter matrices, prior first then one per step.} + +\item{prior}{the dist_list the run used, for the support lines.} +} +\value{ +a ggplot object. +} +\description{ +Ensemble mean and 5-95 % spread of each calibrated parameter at the + prior and after each tempering step, with its declared support drawn in where + finite. Shows whether a parameter settles inside its support or moves onto a + bound, and whether it walks there or jumps. +} +\author{ +Akash BV +} diff --git a/man/plot_target_slots.Rd b/man/plot_target_slots.Rd new file mode 100644 index 0000000..97a9dca --- /dev/null +++ b/man/plot_target_slots.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{plot_target_slots} +\alias{plot_target_slots} +\title{Fitted target slots against observations} +\usage{ +plot_target_slots(G, meta, ylab = "value", title = NULL) +} +\arguments{ +\item{G}{prediction matrix covering the slots in `meta`.} + +\item{meta}{observation meta of the fitted slots.} + +\item{ylab}{y axis label.} + +\item{title}{optional factual page label for a multi page pdf.} +} +\value{ +a ggplot object. +} +\description{ +A contracted target has no year axis (obs_year is NA by + construction), so its slots go on a categorical axis: member points, + ensemble mean, observation +/- sd. +} +\author{ +Akash BV +} diff --git a/man/plot_treatment_effect.Rd b/man/plot_treatment_effect.Rd new file mode 100644 index 0000000..c417727 --- /dev/null +++ b/man/plot_treatment_effect.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{plot_treatment_effect} +\alias{plot_treatment_effect} +\title{Measured against modeled treatment effect, before and after calibration} +\usage{ +plot_treatment_effect(G_prior, G_post, meta, variable, treatment, control) +} +\arguments{ +\item{G_prior, G_post}{prediction matrices on the raw slots (members x slots).} + +\item{meta}{raw observation meta covering the paired slots.} + +\item{variable}{the per treatment per year variable the effect is on.} + +\item{treatment, control}{treatment ids to compare.} +} +\value{ +a ggplot object. +} +\description{ +Per year treatment minus control from the prior and posterior + forward ensembles against the measured contrast, as the absolute effect and + as percent of control. Bands are the 5-95 % ensemble range around the mean; + measurement bars are one standard deviation, the two arms' variances added + (delta method for the relative form). +} +\author{ +Akash BV +} diff --git a/man/plot_validation.Rd b/man/plot_validation.Rd new file mode 100644 index 0000000..2fe4bb1 --- /dev/null +++ b/man/plot_validation.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{plot_validation} +\alias{plot_validation} +\title{Held-out validation: predicted against observed} +\usage{ +plot_validation(G, meta) +} +\arguments{ +\item{G}{prediction matrix covering the validation slots.} + +\item{meta}{observation meta for those slots, carrying `value` and `units`.} +} +\value{ +a ggplot object. +} +\description{ +Ensemble mean prediction against the observation for slots kept out + of the likelihood, with the 5-95 % ensemble range and a 1:1 line. + + Panels are on free scales: variables held out together can span very different + ranges, and a shared scale would flatten the smaller one. +} +\author{ +Akash BV +} diff --git a/man/prior_from_shared_postdistns.Rd b/man/prior_from_shared_postdistns.Rd new file mode 100644 index 0000000..653b57b --- /dev/null +++ b/man/prior_from_shared_postdistns.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/priors.R +\name{prior_from_shared_postdistns} +\alias{prior_from_shared_postdistns} +\title{Prior for a rate shared across several PFTs} +\usage{ +prior_from_shared_postdistns(params, posterior_files) +} +\arguments{ +\item{params}{character vector of PEcAn trait names to calibrate.} + +\item{posterior_files}{named character vector of post.distns paths, one per soil PFT.} +} +\value{ +a dist_list, one record per trait. +} +\description{ +A rate written into several PFTs is one calibrated quantity, so + there must be one prior for it. This reads the requested traits from each + named PFT's posterior and requires them to agree before returning a single + record; disagreement is an error rather than a quiet choice of the first. +} +\author{ +Akash BV +} diff --git a/man/read_cal_val_observations.Rd b/man/read_cal_val_observations.Rd index b91d091..2462d62 100644 --- a/man/read_cal_val_observations.Rd +++ b/man/read_cal_val_observations.Rd @@ -2,19 +2,11 @@ % Please edit documentation in R/observations.R \name{read_cal_val_observations} \alias{read_cal_val_observations} -\title{Read the curated cal/val observations table} +\title{read the observations table from a cal-val-data checkout.} \usage{ read_cal_val_observations(cal_val_dir) } -\arguments{ -\item{cal_val_dir}{directory holding the observations tsv export.} -} -\value{ -tibble of observation rows, value and study_year coerced to numeric. -} \description{ -Reads the single observations tsv exported from the cal/val -workbook, drops empty trailing sheet rows, and strips thousands-separator -commas from the numeric columns so coercion does not silently produce NA. +read the observations table from a cal-val-data checkout. } \keyword{internal} diff --git a/man/repoint_poolinitcond.Rd b/man/repoint_poolinitcond.Rd index e353ace..05ef222 100644 --- a/man/repoint_poolinitcond.Rd +++ b/man/repoint_poolinitcond.Rd @@ -3,12 +3,18 @@ \name{repoint_poolinitcond} \alias{repoint_poolinitcond} \title{point each site's poolinitcond path at the freshly written per particle ics, -so particle j uses its own ic at every site and the design indexes 1:J.} +so particle j uses its own ic at every site and the design indexes 1:J. sites +without a calibrated state (no `ic_paths` entry) keep their template ic, +recycled to J paths so the shared design column stays in range; overwriting +them with an empty list silently drops every one of their run dirs.} \usage{ -repoint_poolinitcond(settings, treatments, ic_paths) +repoint_poolinitcond(settings, treatments, ic_paths, n_particles) } \description{ point each site's poolinitcond path at the freshly written per particle ics, -so particle j uses its own ic at every site and the design indexes 1:J. +so particle j uses its own ic at every site and the design indexes 1:J. sites +without a calibrated state (no `ic_paths` entry) keep their template ic, +recycled to J paths so the shared design column stays in range; overwriting +them with an empty list silently drops every one of their run dirs. } \keyword{internal} diff --git a/man/run_window.Rd b/man/run_window.Rd new file mode 100644 index 0000000..7841150 --- /dev/null +++ b/man/run_window.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/forward_sipnet.R +\name{run_window} +\alias{run_window} +\title{Run years per treatment from a multisite settings object} +\usage{ +run_window(settings) +} +\arguments{ +\item{settings}{a PEcAn multisite settings object.} +} +\value{ +integer matrix (2 x n_treatments), columns named by treatment. +} +\description{ +First and last run year of each treatment, for reading model + output over that treatment's own window; a joint run spans different + periods per site. +} +\author{ +Akash BV +} diff --git a/man/score_iteration.Rd b/man/score_iteration.Rd index 8ec43f3..1eceadd 100644 --- a/man/score_iteration.Rd +++ b/man/score_iteration.Rd @@ -9,15 +9,17 @@ score_iteration(G, meta) \arguments{ \item{G}{matrix (J x P), columns named by observation slot.} -\item{meta}{observation meta (slot, value = obs mean).} +\item{meta}{observation meta (slot, variable, value = obs mean).} } \value{ -one row tibble: rmse, bias, coverage, mean_width, crps. +tibble, one row per variable: slots, rmse, bias, coverage, + mean_width, crps. } \description{ Scores a J x P prediction ensemble G against the observations in -meta: RMSE and bias of the ensemble mean, coverage and mean width of the 90% -band, and the mean CRPS. Columns of G are matched to observations by slot. +meta, per variable: RMSE and bias of the ensemble mean, coverage and mean +width of the 90% band, and the mean CRPS. Columns of G are matched to +observations by slot. } \author{ Akash BV diff --git a/man/state_prior_from_obs.Rd b/man/state_prior_from_obs.Rd index 15b10c4..2e28fd2 100644 --- a/man/state_prior_from_obs.Rd +++ b/man/state_prior_from_obs.Rd @@ -4,10 +4,17 @@ \alias{state_prior_from_obs} \title{Per-site prior for a calibrated initial state, anchored to an observation} \usage{ -state_prior_from_obs(meta, prefix, from_unit, to_unit, anchor_year = NULL) +state_prior_from_obs( + meta, + prefix, + from_unit, + to_unit, + variable, + anchor_year = NULL +) } \arguments{ -\item{meta}{observation meta (treatment_id, study_year, value, var_obs).} +\item{meta}{observation meta (variable, treatment_id, obs_year, value, var_obs).} \item{prefix}{column prefix marking the state entries (e.g. "soilInit.").} @@ -15,7 +22,10 @@ state_prior_from_obs(meta, prefix, from_unit, to_unit, anchor_year = NULL) \item{to_unit}{unit the model state is in (udunits string).} -\item{anchor_year}{study_year to anchor on; defaults to the earliest.} +\item{variable}{the observed variable whose slots anchor the state; the joint +meta can span several variables and an unscoped anchor would take the wrong one.} + +\item{anchor_year}{observation year to anchor on; defaults to the earliest.} } \value{ a dist_list keyed , in site order. diff --git a/man/subset_obs.Rd b/man/subset_obs.Rd new file mode 100644 index 0000000..63f0644 --- /dev/null +++ b/man/subset_obs.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observations.R +\name{subset_obs} +\alias{subset_obs} +\title{Split an observation target into fitted and validation parts} +\usage{ +subset_obs(obs, keep) +} +\arguments{ +\item{obs}{a build_obs target list(y, Sigma, meta).} + +\item{keep}{logical vector over rows of `obs$meta`, TRUE to retain.} +} +\value{ +an obs list of the same shape carrying only the kept slots. +} +\description{ +Keeps slots in the run without letting them into the likelihood; + a held-out variable is still predicted and still scored, which is what a + validation target is. +} +\author{ +Akash BV +} diff --git a/man/target_rows.Rd b/man/target_rows.Rd new file mode 100644 index 0000000..fb8765a --- /dev/null +++ b/man/target_rows.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/observations.R +\name{target_rows} +\alias{target_rows} +\title{filter, unit-check and key one target's rows (see build_obs for the spec).} +\usage{ +target_rows(raw_all, tg) +} +\description{ +filter, unit-check and key one target's rows (see build_obs for the spec). +} +\keyword{internal} diff --git a/man/write_samples_rdata.Rd b/man/write_samples_rdata.Rd deleted file mode 100644 index 77706ab..0000000 --- a/man/write_samples_rdata.Rd +++ /dev/null @@ -1,12 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/forward_sipnet.R -\name{write_samples_rdata} -\alias{write_samples_rdata} -\title{write samples.Rdata in the object run.write.configs expects.} -\usage{ -write_samples_rdata(ensemble.samples, file) -} -\description{ -write samples.Rdata in the object run.write.configs expects. -} -\keyword{internal}