Refactor the baseline Validation Pipeline to align with GSoC Architecture

modules/benchmark/NAMESPACE

@@ -35,6 +35,7 @@ export(metric_run)
 export(metric_scatter_plot)
 export(metric_timeseries_plot)
 export(read_settings_BRR)
+export(run_benchmark)
 importFrom(ggplot2,geom_path)
 importFrom(ggplot2,geom_point)
 importFrom(ggplot2,ggplot)

modules/benchmark/R/data_intake.R

@@ -0,0 +1,40 @@
+#' Load and standardize arbitrary tabular data using a YAML mapping configuration
+#'
+#' @param data.path character, file path to the tabular data (e.g. .csv)
+#' @param mapping.path character, file path to the YAML mapping configuration
+#' @return A standardized data frame with column names mapped to PEcAn standard vocabulary
+#' @export
+#' @importFrom yaml read_yaml
+#' @importFrom dplyr rename
+load_and_map_data <- function(data.path, mapping.path) {
+  # Load the raw data (currently assuming CSV, but could be extended to NetCDF)
+  dat <- utils::read.csv(data.path, as.is = TRUE, check.names = FALSE)
+
+  # Load the YAML mapping
+  # The YAML should look like:
+  # variables:
+  #   airT: TA_F
+  #   NEE: NEE_PI
+  mapping <- yaml::read_yaml(mapping.path)
+
+  if (is.null(mapping$variables)) {
+    stop("YAML mapping must contain a 'variables' section.")
+  }
+
+  # Create a named vector for dplyr::rename (new_name = old_name)
+  rename_vector <- unlist(mapping$variables)
+
+  # Only rename columns that exist in the raw data
+  valid_renames <- rename_vector[rename_vector %in% colnames(dat)]
+
+  # Apply renaming
+  if (length(valid_renames) > 0) {
+    # dplyr::rename syntax expects: rename(df, new_name = old_name)
+    # Using tidy evaluation with !!!
+    dat <- dplyr::rename(dat, !!!valid_renames)
+  } else {
+    warning("No matching columns found in the dataset to map.")
+  }
+
+  return(dat)
+}

modules/benchmark/R/run_benchmark.R

@@ -0,0 +1,129 @@
+##' Run a simple benchmark pipeline
+##'
+##' Takes two validated dataframes, aligns by time,
+##' computes metrics, and returns a results table with a plot.
+##'
+##' @param model_df data.frame with columns: time (POSIXct), value (numeric)
+##' @param obs_df   data.frame with columns: time (POSIXct), value (numeric)
+##' @param metrics character vector of metrics to compute. Options: "RMSE", "MAE"
+##' @param tolerance_secs nearest-neighbor time tolerance in seconds (default 1 hour)
+##' @param method alignment method: "nearest" or "interpolate"
+##'
+##' @return list with: metrics (data.frame), aligned (data.frame), plot (ggplot)
+##' @export
+##' @author Anshul Jain
+run_benchmark <- function(model_df, obs_df,
+                          metrics = c("RMSE", "MAE"),
+                          tolerance_secs = 3600,
+                          method = "nearest") {
+
+  # Stage 1: Validate schema
+  bm_validate(model_df, obs_df)
+
+  # Stage 2: Align by time
+  aligned <- align_by_time(model_df, obs_df, tolerance_secs = tolerance_secs)
+
+  # Stage 3: Compute metrics via registry
+  results <- compute_metrics(aligned, metrics)
+
+  # Stage 4: Plot
+  plot <- plot_time_series(aligned)
+
+  list(metrics = results, aligned = aligned, plot = plot)
+}
+
+##' Validate benchmark input dataframes
+##'
+##' @param model_df data.frame with columns: time (POSIXct), value (numeric)
+##' @param obs_df   data.frame with columns: time (POSIXct), value (numeric)
+##' @return invisible(TRUE)
+bm_validate <- function(model_df, obs_df) {
+  for (df in list(model_df, obs_df)) {
+    if (!inherits(df$time, "POSIXct"))
+      stop("Column 'time' must be POSIXct, got: ", class(df$time))
+    if (!is.numeric(df$value))
+      stop("Column 'value' must be numeric, got: ", class(df$value))
+  }
+  invisible(TRUE)
+}
+
+##' Align model and observation data frames by nearest time
+##'
+##' @param model_df data.frame with columns: time (POSIXct), value
+##' @param obs_df   data.frame with columns: time (POSIXct), value
+##' @param tolerance_secs max allowed time difference in seconds
+##'
+##' @return data.frame with columns: time, model, obs
+align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) {
+  # Sort both dataframes by time to ensure findInterval works correctly
+  model_df <- model_df[order(model_df$time), ]
+  obs_df <- obs_df[order(obs_df$time), ]
+
+  # For each model time, find the interval in obs_time it falls into
+  idx <- findInterval(model_df$time, obs_df$time, all.inside = TRUE)
+
+  # findInterval returns index i where obs[i] <= model_time < obs[i+1]
+  # We check both i and i+1 to see which one is the absolute nearest
+  idx_next <- pmin(idx + 1, nrow(obs_df))
+
+  diff_current <- abs(as.numeric(difftime(model_df$time, obs_df$time[idx], units = "secs")))
+  diff_next <- abs(as.numeric(difftime(model_df$time, obs_df$time[idx_next], units = "secs")))
+
+  # Select the index of the closest observation
+  nearest_idx <- ifelse(diff_current <= diff_next, idx, idx_next)
+  time_diffs <- pmin(diff_current, diff_next)
+
+  # Filter by our time tolerance
+  valid <- time_diffs <= tolerance_secs
+
+  # Construct the aligned base data.frame
+  aligned <- data.frame(
+    time = model_df$time[valid],
+    model = model_df$value[valid],
+    obs = obs_df$value[nearest_idx][valid]
+  )
+
+  return(aligned)
+}
+
+##' Compute benchmark metrics
+##'
+##' @param aligned data.frame with columns: time, model, obs
+##' @param metrics character vector of metric names
+##' @return data.frame with columns: metric, value
+compute_metrics <- function(aligned, metrics = c("RMSE", "MAE", "R2")) {
+  # Future-proofing: Functions in the registry now accept the full aligned dataframe
+  # This aligns with the decoupled metric architecture introduced in PR #3888
+  METRIC_REGISTRY <- list(
+    RMSE = function(dat) sqrt(mean((dat$model - dat$obs)^2, na.rm = TRUE)),
+    MAE  = function(dat) mean(abs(dat$model - dat$obs), na.rm = TRUE),
+    R2   = function(dat) {
+      if (exists("metric_R2", where = asNamespace("PEcAn.benchmark"), mode = "function")) {
+        return(PEcAn.benchmark::metric_R2(dat))
+      }
+      # Fallback if PR #3888 is not yet merged
+      numer <- sum((dat$obs - mean(dat$obs, na.rm=T)) * (dat$model - mean(dat$model, na.rm=T)), na.rm=T)
+      denom <- sqrt(sum((dat$obs - mean(dat$obs, na.rm=T))^2, na.rm=T)) * sqrt(sum((dat$model - mean(dat$model, na.rm=T))^2, na.rm=T))
+      (numer / denom)^2
+    }
+  )
+
+  results <- lapply(toupper(metrics), function(m) {
+    if (!m %in% names(METRIC_REGISTRY)) stop("Unknown metric: ", m)
+    METRIC_REGISTRY[[m]](aligned)
+  })
+
+  data.frame(metric = toupper(metrics), value = unlist(results, use.names = FALSE))
+}
+
+##' Plot model vs observations time series
+##'
+##' @param aligned data.frame with columns: time, model, obs
+##' @return ggplot object
+plot_time_series <- function(aligned) {
+  ggplot2::ggplot(aligned, ggplot2::aes(x = .data$time)) +
+    ggplot2::geom_line(ggplot2::aes(y = .data$model, color = "Model")) +
+    ggplot2::geom_line(ggplot2::aes(y = .data$obs, color = "Obs")) +
+    ggplot2::labs(color = "", y = "value", title = "Model vs Observations") +
+    ggplot2::theme_bw()
+}

modules/benchmark/README.md

@@ -1,4 +1,39 @@
+## Quickstart: run_benchmark()
 
+`run_benchmark()` is a simple entry point that loads model output and
+observations, aligns them by time, computes metrics, and returns a plot.
+
+### Input format
+
+Both input files must be CSV with two columns:
+- `time` — timestamp (e.g. `2020-01-01 00:00:00`)
+- `value` — numeric variable value
+
+### Usage
+```r
+library(PEcAn.benchmark)
+
+res <- run_benchmark(
+  model_path = "inst/testdata/sample_model.csv",
+  obs_path   = "inst/testdata/sample_obs.csv"
+)
+
+# View metrics
+print(res$metrics)
+#   metric     value
+# 1   RMSE 0.1322876
+# 2    MAE 0.1250000
+
+# View plot
+res$plot
+```
+
+### Parameters
+
+- `model_path` — path to model output CSV
+- `obs_path` — path to observations CSV  
+- `metrics` — vector of metrics to compute: `"RMSE"`, `"MAE"` (default: both)
+- `tolerance_secs` — max time difference for matching (default: 3600 seconds)
 # PEcAn.benchmark
 
 <!-- badges: start -->

modules/benchmark/inst/testdata/sample_model.csv

@@ -0,0 +1,5 @@
+time,value
+2020-01-01 00:00:00,1.0
+2020-01-01 01:00:00,2.0
+2020-01-01 02:00:00,3.0
+2020-01-01 03:00:00,4.0

modules/benchmark/inst/testdata/sample_obs.csv

@@ -0,0 +1,5 @@
+time,value
+2020-01-01 00:00:00,1.1
+2020-01-01 01:00:00,1.9
+2020-01-01 02:00:00,3.2
+2020-01-01 03:00:00,3.9

modules/benchmark/tests/testthat/test-run_benchmark.R

@@ -0,0 +1,40 @@
+library(testthat)
+
+model_df <- data.frame(
+  time  = as.POSIXct(seq(0, 3600*3, by = 3600), origin = "1970-01-01", tz = "UTC"),
+  value = c(1, 2, 3, 4)
+)
+obs_df <- data.frame(
+  time  = as.POSIXct(seq(0, 3600*3, by = 3600), origin = "1970-01-01", tz = "UTC"),
+  value = c(1.1, 1.9, 3.2, 3.9)
+)
+
+test_that("run_benchmark returns correct structure", {
+  res <- run_benchmark(model_df, obs_df, metrics = c("RMSE", "MAE"))
+  expect_true("metrics" %in% names(res))
+  expect_true("aligned" %in% names(res))
+  expect_true("plot"    %in% names(res))
+  expect_equal(nrow(res$metrics), 2)
+})
+
+test_that("bm_validate rejects bad input", {
+  bad_df <- data.frame(time = c("2023-01-01"), value = c(1.0))
+  expect_error(bm_validate(bad_df, obs_df), "POSIXct")
+})
+
+test_that("compute_metrics returns correct values", {
+  aligned <- data.frame(
+    time  = model_df$time,
+    model = c(1, 2, 3, 4),
+    obs   = c(1, 2, 3, 4)
+  )
+  res <- compute_metrics(aligned, c("RMSE", "MAE"))
+  expect_equal(res$value[res$metric == "RMSE"], 0)
+  expect_equal(res$value[res$metric == "MAE"],  0)
+})
+
+test_that("align_by_time matches exact timestamps", {
+  aligned <- align_by_time(model_df, obs_df)
+  expect_equal(nrow(aligned), 4)
+  expect_true(all(c("time", "model", "obs") %in% names(aligned)))
+})