From 7df98d90d54a9c6a200d5a5d14d3118ef492d005 Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Mon, 23 Mar 2026 22:57:38 +0530 Subject: [PATCH 01/27] feat(benchmark): add run_benchmark MVP with alignment, metrics, and test data --- modules/benchmark/R/run_benchmark.R | 72 +++++++++++++++++++ .../benchmark/inst/testdata/sample_model.csv | 5 ++ .../benchmark/inst/testdata/sample_obs.csv | 5 ++ 3 files changed, 82 insertions(+) create mode 100644 modules/benchmark/R/run_benchmark.R create mode 100644 modules/benchmark/inst/testdata/sample_model.csv create mode 100644 modules/benchmark/inst/testdata/sample_obs.csv diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R new file mode 100644 index 00000000000..8cc99bf6765 --- /dev/null +++ b/modules/benchmark/R/run_benchmark.R @@ -0,0 +1,72 @@ +##' Run a simple benchmark pipeline +##' +##' Loads model output and observations, aligns by time, +##' computes RMSE and MAE, and returns a results table with a plot. +##' +##' @param model_path path to model output CSV file (must have 'time' and 'value' columns) +##' @param obs_path path to observations CSV file (must have 'time' and 'value' columns) +##' @param metrics character vector of metrics to compute. Options: "RMSE", "MAE" +##' @param tolerance_secs nearest-neighbor time tolerance in seconds (default 1 hour) +##' +##' @return list with: metrics (data.frame), aligned (data.frame), plot (ggplot) +##' @export +##' +##' @author Your Name +run_benchmark <- function(model_path, obs_path, + metrics = c("RMSE", "MAE"), + tolerance_secs = 3600) { + + # --- Load data --- + model_df <- read.csv(model_path, stringsAsFactors = FALSE) + obs_df <- read.csv(obs_path, stringsAsFactors = FALSE) + + # --- Ensure time column is POSIXct --- + model_df$time <- as.POSIXct(model_df$time, tz = "UTC") + obs_df$time <- as.POSIXct(obs_df$time, tz = "UTC") + + # --- Align by nearest time --- + aligned <- align_by_time(model_df, obs_df, tolerance_secs = tolerance_secs) + + # --- Compute metrics --- + results <- list() + for (m in toupper(metrics)) { + results[[m]] <- switch(m, + "RMSE" = sqrt(mean((aligned$model - aligned$obs)^2, na.rm = TRUE)), + "MAE" = mean(abs(aligned$model - aligned$obs), na.rm = TRUE), + stop("Unknown metric: ", m) + ) + } + metrics_df <- data.frame(metric = names(results), + value = unlist(results, use.names = FALSE)) + + # --- Plot --- + plot <- ggplot2::ggplot(aligned, ggplot2::aes(x = time)) + + ggplot2::geom_line(ggplot2::aes(y = model, color = "model")) + + ggplot2::geom_line(ggplot2::aes(y = obs, color = "obs")) + + ggplot2::labs(color = "", y = "value", title = "Model vs Observations") + + list(metrics = metrics_df, aligned = aligned, plot = plot) +} + + +##' 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) { + aligned <- do.call(rbind, lapply(seq_len(nrow(model_df)), function(i) { + diffs <- abs(as.numeric(difftime(obs_df$time, model_df$time[i], units = "secs"))) + nearest <- which.min(diffs) + if (diffs[nearest] <= tolerance_secs) { + data.frame(time = model_df$time[i], + model = model_df$value[i], + obs = obs_df$value[nearest]) + } else { + NULL + } + })) + aligned +} diff --git a/modules/benchmark/inst/testdata/sample_model.csv b/modules/benchmark/inst/testdata/sample_model.csv new file mode 100644 index 00000000000..cf5df29d50d --- /dev/null +++ b/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 diff --git a/modules/benchmark/inst/testdata/sample_obs.csv b/modules/benchmark/inst/testdata/sample_obs.csv new file mode 100644 index 00000000000..f24d28f7c26 --- /dev/null +++ b/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 From b2731cd5ad5968aff3bbda4077da6ded1323d1e1 Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Mon, 23 Mar 2026 23:16:03 +0530 Subject: [PATCH 02/27] feat(benchmark): add tests and update README with quickstart --- modules/benchmark/README.md | 35 +++++++++++++++++++ .../tests/testthat/test-run_benchmark.R | 21 +++++++++++ 2 files changed, 56 insertions(+) create mode 100644 modules/benchmark/tests/testthat/test-run_benchmark.R diff --git a/modules/benchmark/README.md b/modules/benchmark/README.md index a8fd53648d8..69bdb5e2f3d 100644 --- a/modules/benchmark/README.md +++ b/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 diff --git a/modules/benchmark/tests/testthat/test-run_benchmark.R b/modules/benchmark/tests/testthat/test-run_benchmark.R new file mode 100644 index 00000000000..746e6f370d7 --- /dev/null +++ b/modules/benchmark/tests/testthat/test-run_benchmark.R @@ -0,0 +1,21 @@ +library(testthat) + +test_that("run_benchmark basic works", { + model <- data.frame( + time = as.POSIXct(seq(0, 3600*3, by = 3600), origin = "1970-01-01", tz = "UTC"), + value = c(1, 2, 3, 4) + ) + obs <- 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) + ) + tmp1 <- tempfile(fileext = ".csv") + tmp2 <- tempfile(fileext = ".csv") + write.csv(model, tmp1, row.names = FALSE) + write.csv(obs, tmp2, row.names = FALSE) + + res <- run_benchmark(tmp1, tmp2, metrics = c("RMSE", "MAE")) + expect_true("metrics" %in% names(res)) + expect_true("aligned" %in% names(res)) + expect_true(nrow(res$metrics) == 2) +}) From 6769cf917110e7aefc99864fc6e599d1481a6ef4 Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 25 Mar 2026 12:07:39 +0530 Subject: [PATCH 03/27] docs(benchmark): add roxygen man page for run_benchmark --- modules/benchmark/man/run_benchmark.Rd | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 modules/benchmark/man/run_benchmark.Rd diff --git a/modules/benchmark/man/run_benchmark.Rd b/modules/benchmark/man/run_benchmark.Rd new file mode 100644 index 00000000000..cb742bc35e6 --- /dev/null +++ b/modules/benchmark/man/run_benchmark.Rd @@ -0,0 +1,19 @@ +\name{run_benchmark} +\alias{run_benchmark} +\title{Run a simple benchmark pipeline} +\usage{ +run_benchmark(model_path, obs_path, metrics = c("RMSE", "MAE"), tolerance_secs = 3600) +} +\arguments{ +\item{model_path}{path to model output CSV file} +\item{obs_path}{path to observations CSV file} +\item{metrics}{character vector of metrics to compute} +\item{tolerance_secs}{nearest-neighbor time tolerance in seconds} +} +\value{ +list with: metrics (data.frame), aligned (data.frame), plot (ggplot) +} +\description{ +Loads model output and observations, aligns by time, +computes RMSE and MAE, and returns a results table with a plot. +} From f2aa206322fea2773ed60e91f9e9e1a4b67550d5 Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 1 Apr 2026 12:11:23 +0530 Subject: [PATCH 04/27] docs(benchmark): add man pages and update NAMESPACE --- modules/benchmark/NAMESPACE | 1 + modules/benchmark/man/align_by_time.Rd | 21 +++++++++++++++++++++ modules/benchmark/man/run_benchmark.Rd | 16 ++++++++++++---- 3 files changed, 34 insertions(+), 4 deletions(-) create mode 100644 modules/benchmark/man/align_by_time.Rd diff --git a/modules/benchmark/NAMESPACE b/modules/benchmark/NAMESPACE index 503822f77bb..5d98b4cc2cf 100644 --- a/modules/benchmark/NAMESPACE +++ b/modules/benchmark/NAMESPACE @@ -41,3 +41,4 @@ importFrom(ggplot2,ggplot) importFrom(ggplot2,labs) importFrom(magrittr,"%>%") importFrom(rlang,.data) +export(run_benchmark) diff --git a/modules/benchmark/man/align_by_time.Rd b/modules/benchmark/man/align_by_time.Rd new file mode 100644 index 00000000000..33e313bb5b5 --- /dev/null +++ b/modules/benchmark/man/align_by_time.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R +\name{align_by_time} +\alias{align_by_time} +\title{Align model and observation data frames by nearest time} +\usage{ +align_by_time(model_df, obs_df, tolerance_secs = 3600) +} +\arguments{ +\item{model_df}{data.frame with columns: time (POSIXct), value} + +\item{obs_df}{data.frame with columns: time (POSIXct), value} + +\item{tolerance_secs}{max allowed time difference in seconds} +} +\value{ +data.frame with columns: time, model, obs +} +\description{ +Align model and observation data frames by nearest time +} diff --git a/modules/benchmark/man/run_benchmark.Rd b/modules/benchmark/man/run_benchmark.Rd index cb742bc35e6..d7d7470694d 100644 --- a/modules/benchmark/man/run_benchmark.Rd +++ b/modules/benchmark/man/run_benchmark.Rd @@ -1,3 +1,5 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R \name{run_benchmark} \alias{run_benchmark} \title{Run a simple benchmark pipeline} @@ -5,10 +7,13 @@ run_benchmark(model_path, obs_path, metrics = c("RMSE", "MAE"), tolerance_secs = 3600) } \arguments{ -\item{model_path}{path to model output CSV file} -\item{obs_path}{path to observations CSV file} -\item{metrics}{character vector of metrics to compute} -\item{tolerance_secs}{nearest-neighbor time tolerance in seconds} +\item{model_path}{path to model output CSV file (must have 'time' and 'value' columns)} + +\item{obs_path}{path to observations CSV file (must have 'time' and 'value' columns)} + +\item{metrics}{character vector of metrics to compute. Options: "RMSE", "MAE"} + +\item{tolerance_secs}{nearest-neighbor time tolerance in seconds (default 1 hour)} } \value{ list with: metrics (data.frame), aligned (data.frame), plot (ggplot) @@ -17,3 +22,6 @@ list with: metrics (data.frame), aligned (data.frame), plot (ggplot) Loads model output and observations, aligns by time, computes RMSE and MAE, and returns a results table with a plot. } +\author{ +Your Name +} From d53d3726b757c98245ff5ab63996f8b7067e4618 Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 1 Apr 2026 12:25:56 +0530 Subject: [PATCH 05/27] fix(benchmark): fix NAMESPACE export and run_benchmark.Rd usage format --- modules/benchmark/NAMESPACE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/benchmark/NAMESPACE b/modules/benchmark/NAMESPACE index 5d98b4cc2cf..fc8f6f45eda 100644 --- a/modules/benchmark/NAMESPACE +++ b/modules/benchmark/NAMESPACE @@ -35,10 +35,10 @@ 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) importFrom(ggplot2,labs) importFrom(magrittr,"%>%") importFrom(rlang,.data) -export(run_benchmark) From d56416b948845caaaf132288f8c2be70ef3edfef Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 1 Apr 2026 12:42:15 +0530 Subject: [PATCH 06/27] fix(benchmark): use multi-line usage format and fix author name --- modules/benchmark/R/run_benchmark.R | 2 +- modules/benchmark/man/run_benchmark.Rd | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index 8cc99bf6765..c4b7004c1e0 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -11,7 +11,7 @@ ##' @return list with: metrics (data.frame), aligned (data.frame), plot (ggplot) ##' @export ##' -##' @author Your Name +##' @author Anshul Jain run_benchmark <- function(model_path, obs_path, metrics = c("RMSE", "MAE"), tolerance_secs = 3600) { diff --git a/modules/benchmark/man/run_benchmark.Rd b/modules/benchmark/man/run_benchmark.Rd index d7d7470694d..b1998b4ea5d 100644 --- a/modules/benchmark/man/run_benchmark.Rd +++ b/modules/benchmark/man/run_benchmark.Rd @@ -4,7 +4,12 @@ \alias{run_benchmark} \title{Run a simple benchmark pipeline} \usage{ -run_benchmark(model_path, obs_path, metrics = c("RMSE", "MAE"), tolerance_secs = 3600) +run_benchmark( + model_path, + obs_path, + metrics = c("RMSE", "MAE"), + tolerance_secs = 3600 +) } \arguments{ \item{model_path}{path to model output CSV file (must have 'time' and 'value' columns)} @@ -23,5 +28,5 @@ Loads model output and observations, aligns by time, computes RMSE and MAE, and returns a results table with a plot. } \author{ -Your Name +Anshul Jain } From 4df7daf2eb4742d2251d69bf824c18c3cfc39544 Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 1 Apr 2026 12:46:40 +0530 Subject: [PATCH 07/27] refactor(benchmark): dataframe-first API, add bm_validate/compute_metrics/plot_time_series, update tests --- modules/benchmark/R/run_benchmark.R | 91 ++++++++++++------- .../tests/testthat/test-run_benchmark.R | 49 +++++++--- 2 files changed, 93 insertions(+), 47 deletions(-) diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index c4b7004c1e0..9897b59b015 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -1,53 +1,51 @@ ##' Run a simple benchmark pipeline ##' -##' Loads model output and observations, aligns by time, -##' computes RMSE and MAE, and returns a results table with a plot. +##' Takes two validated dataframes, aligns by time, +##' computes metrics, and returns a results table with a plot. ##' -##' @param model_path path to model output CSV file (must have 'time' and 'value' columns) -##' @param obs_path path to observations CSV file (must have 'time' and 'value' columns) +##' @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_path, obs_path, +run_benchmark <- function(model_df, obs_df, metrics = c("RMSE", "MAE"), - tolerance_secs = 3600) { - - # --- Load data --- - model_df <- read.csv(model_path, stringsAsFactors = FALSE) - obs_df <- read.csv(obs_path, stringsAsFactors = FALSE) + tolerance_secs = 3600, + method = "nearest") { - # --- Ensure time column is POSIXct --- - model_df$time <- as.POSIXct(model_df$time, tz = "UTC") - obs_df$time <- as.POSIXct(obs_df$time, tz = "UTC") + # Stage 1: Validate schema + bm_validate(model_df, obs_df) - # --- Align by nearest time --- + # Stage 2: Align by time aligned <- align_by_time(model_df, obs_df, tolerance_secs = tolerance_secs) - # --- Compute metrics --- - results <- list() - for (m in toupper(metrics)) { - results[[m]] <- switch(m, - "RMSE" = sqrt(mean((aligned$model - aligned$obs)^2, na.rm = TRUE)), - "MAE" = mean(abs(aligned$model - aligned$obs), na.rm = TRUE), - stop("Unknown metric: ", m) - ) - } - metrics_df <- data.frame(metric = names(results), - value = unlist(results, use.names = FALSE)) + # Stage 3: Compute metrics via registry + results <- compute_metrics(aligned, metrics) - # --- Plot --- - plot <- ggplot2::ggplot(aligned, ggplot2::aes(x = time)) + - ggplot2::geom_line(ggplot2::aes(y = model, color = "model")) + - ggplot2::geom_line(ggplot2::aes(y = obs, color = "obs")) + - ggplot2::labs(color = "", y = "value", title = "Model vs Observations") + # Stage 4: Plot + plot <- plot_time_series(aligned) - list(metrics = metrics_df, aligned = aligned, plot = plot) + 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 ##' @@ -70,3 +68,32 @@ align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { })) 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")) { + METRIC_REGISTRY <- list( + RMSE = function(x, y) sqrt(mean((x - y)^2, na.rm = TRUE)), + MAE = function(x, y) mean(abs(x - y), na.rm = TRUE) + ) + results <- lapply(toupper(metrics), function(m) { + if (!m %in% names(METRIC_REGISTRY)) stop("Unknown metric: ", m) + METRIC_REGISTRY[[m]](aligned$model, aligned$obs) + }) + 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 = time)) + + ggplot2::geom_line(ggplot2::aes(y = model, color = "Model")) + + ggplot2::geom_line(ggplot2::aes(y = obs, color = "Obs")) + + ggplot2::labs(color = "", y = "value", title = "Model vs Observations") + + ggplot2::theme_bw() +} diff --git a/modules/benchmark/tests/testthat/test-run_benchmark.R b/modules/benchmark/tests/testthat/test-run_benchmark.R index 746e6f370d7..765dfe1cb02 100644 --- a/modules/benchmark/tests/testthat/test-run_benchmark.R +++ b/modules/benchmark/tests/testthat/test-run_benchmark.R @@ -1,21 +1,40 @@ library(testthat) -test_that("run_benchmark basic works", { - model <- data.frame( - time = as.POSIXct(seq(0, 3600*3, by = 3600), origin = "1970-01-01", tz = "UTC"), - value = c(1, 2, 3, 4) - ) - obs <- 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) - ) - tmp1 <- tempfile(fileext = ".csv") - tmp2 <- tempfile(fileext = ".csv") - write.csv(model, tmp1, row.names = FALSE) - write.csv(obs, tmp2, row.names = FALSE) +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) +) - res <- run_benchmark(tmp1, tmp2, metrics = c("RMSE", "MAE")) +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(nrow(res$metrics) == 2) + 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))) }) From bac2138624163affd224c13e73b2677a4e4fd3ee Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 1 Apr 2026 13:02:18 +0530 Subject: [PATCH 08/27] docs(benchmark): add man pages for bm_validate, compute_metrics, plot_time_series --- modules/benchmark/man/bm_validate.Rd | 19 +++++++++++++++++++ modules/benchmark/man/compute_metrics.Rd | 19 +++++++++++++++++++ modules/benchmark/man/plot_time_series.Rd | 17 +++++++++++++++++ modules/benchmark/man/run_benchmark.Rd | 17 ++++++++++------- 4 files changed, 65 insertions(+), 7 deletions(-) create mode 100644 modules/benchmark/man/bm_validate.Rd create mode 100644 modules/benchmark/man/compute_metrics.Rd create mode 100644 modules/benchmark/man/plot_time_series.Rd diff --git a/modules/benchmark/man/bm_validate.Rd b/modules/benchmark/man/bm_validate.Rd new file mode 100644 index 00000000000..11a4fd2eb36 --- /dev/null +++ b/modules/benchmark/man/bm_validate.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R +\name{bm_validate} +\alias{bm_validate} +\title{Validate benchmark input dataframes} +\usage{ +bm_validate(model_df, obs_df) +} +\arguments{ +\item{model_df}{data.frame with columns: time (POSIXct), value (numeric)} + +\item{obs_df}{data.frame with columns: time (POSIXct), value (numeric)} +} +\value{ +invisible(TRUE) +} +\description{ +Validate benchmark input dataframes +} diff --git a/modules/benchmark/man/compute_metrics.Rd b/modules/benchmark/man/compute_metrics.Rd new file mode 100644 index 00000000000..a9a66fe98e9 --- /dev/null +++ b/modules/benchmark/man/compute_metrics.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R +\name{compute_metrics} +\alias{compute_metrics} +\title{Compute benchmark metrics} +\usage{ +compute_metrics(aligned, metrics = c("RMSE", "MAE")) +} +\arguments{ +\item{aligned}{data.frame with columns: time, model, obs} + +\item{metrics}{character vector of metric names} +} +\value{ +data.frame with columns: metric, value +} +\description{ +Compute benchmark metrics +} diff --git a/modules/benchmark/man/plot_time_series.Rd b/modules/benchmark/man/plot_time_series.Rd new file mode 100644 index 00000000000..0710201ed1a --- /dev/null +++ b/modules/benchmark/man/plot_time_series.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R +\name{plot_time_series} +\alias{plot_time_series} +\title{Plot model vs observations time series} +\usage{ +plot_time_series(aligned) +} +\arguments{ +\item{aligned}{data.frame with columns: time, model, obs} +} +\value{ +ggplot object +} +\description{ +Plot model vs observations time series +} diff --git a/modules/benchmark/man/run_benchmark.Rd b/modules/benchmark/man/run_benchmark.Rd index b1998b4ea5d..3d15828a2b4 100644 --- a/modules/benchmark/man/run_benchmark.Rd +++ b/modules/benchmark/man/run_benchmark.Rd @@ -5,27 +5,30 @@ \title{Run a simple benchmark pipeline} \usage{ run_benchmark( - model_path, - obs_path, + model_df, + obs_df, metrics = c("RMSE", "MAE"), - tolerance_secs = 3600 + tolerance_secs = 3600, + method = "nearest" ) } \arguments{ -\item{model_path}{path to model output CSV file (must have 'time' and 'value' columns)} +\item{model_df}{data.frame with columns: time (POSIXct), value (numeric)} -\item{obs_path}{path to observations CSV file (must have 'time' and 'value' columns)} +\item{obs_df}{data.frame with columns: time (POSIXct), value (numeric)} \item{metrics}{character vector of metrics to compute. Options: "RMSE", "MAE"} \item{tolerance_secs}{nearest-neighbor time tolerance in seconds (default 1 hour)} + +\item{method}{alignment method: "nearest" or "interpolate"} } \value{ list with: metrics (data.frame), aligned (data.frame), plot (ggplot) } \description{ -Loads model output and observations, aligns by time, -computes RMSE and MAE, and returns a results table with a plot. +Takes two validated dataframes, aligns by time, +computes metrics, and returns a results table with a plot. } \author{ Anshul Jain From 1a774479d66df44c030fc5d098240b6d67af935b Mon Sep 17 00:00:00 2001 From: anshul23102 Date: Wed, 1 Apr 2026 13:16:23 +0530 Subject: [PATCH 09/27] fix(benchmark): use .data$ in aes() to fix R CMD check NOTE --- modules/benchmark/R/run_benchmark.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index 9897b59b015..03128bd4cd7 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -91,9 +91,9 @@ compute_metrics <- function(aligned, metrics = c("RMSE", "MAE")) { ##' @param aligned data.frame with columns: time, model, obs ##' @return ggplot object plot_time_series <- function(aligned) { - ggplot2::ggplot(aligned, ggplot2::aes(x = time)) + - ggplot2::geom_line(ggplot2::aes(y = model, color = "Model")) + - ggplot2::geom_line(ggplot2::aes(y = obs, color = "Obs")) + + 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() } From 0887a471282bda4694aa87e79a68a8e0fa28bfb7 Mon Sep 17 00:00:00 2001 From: ayushman1210 Date: Sun, 31 May 2026 15:47:37 +0530 Subject: [PATCH 10/27] Phase 2: Refactor benchmark pipeline and add data intake API --- modules/benchmark/R/data_intake.R | 40 +++++++++++++++++++ modules/benchmark/R/run_benchmark.R | 59 +++++++++++++++++++++-------- 2 files changed, 83 insertions(+), 16 deletions(-) create mode 100644 modules/benchmark/R/data_intake.R diff --git a/modules/benchmark/R/data_intake.R b/modules/benchmark/R/data_intake.R new file mode 100644 index 00000000000..1a9e97e5178 --- /dev/null +++ b/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) +} diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index 03128bd4cd7..ed9aa03037e 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -55,18 +55,32 @@ bm_validate <- function(model_df, obs_df) { ##' ##' @return data.frame with columns: time, model, obs align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { - aligned <- do.call(rbind, lapply(seq_len(nrow(model_df)), function(i) { - diffs <- abs(as.numeric(difftime(obs_df$time, model_df$time[i], units = "secs"))) - nearest <- which.min(diffs) - if (diffs[nearest] <= tolerance_secs) { - data.frame(time = model_df$time[i], - model = model_df$value[i], - obs = obs_df$value[nearest]) - } else { - NULL - } - })) - aligned + # Ensure data.table is available + if (!requireNamespace("data.table", quietly = TRUE)) { + stop("Package 'data.table' is required for high-performance alignment.") + } + + # Convert to data.table and explicitly name value columns + dt_model <- data.table::data.table(time = model_df$time, model = model_df$value) + # Keep original obs time in a separate column to calculate difference after join + dt_obs <- data.table::data.table(time = obs_df$time, obs = obs_df$value, obs_time = obs_df$time) + + # Set keys for fast roll-join + data.table::setkey(dt_model, time) + data.table::setkey(dt_obs, time) + + # Perform rolling join to nearest time + aligned <- dt_obs[dt_model, roll = "nearest"] + + # Filter by tolerance + aligned[, time_diff := abs(as.numeric(difftime(obs_time, time, units = "secs")))] + aligned <- aligned[time_diff <= tolerance_secs] + + # Format output to match expected schema + aligned <- aligned[, .(time, model, obs)] + data.table::setDF(aligned) # Convert back to base data.frame + + return(aligned) } ##' Compute benchmark metrics @@ -74,15 +88,28 @@ align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { ##' @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")) { +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(x, y) sqrt(mean((x - y)^2, na.rm = TRUE)), - MAE = function(x, y) mean(abs(x - y), na.rm = TRUE) + 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$model, aligned$obs) + METRIC_REGISTRY[[m]](aligned) }) + data.frame(metric = toupper(metrics), value = unlist(results, use.names = FALSE)) } From e24195bbdda5c138f8d4878b74611746234f6197 Mon Sep 17 00:00:00 2001 From: ayushman1210 Date: Wed, 10 Jun 2026 23:06:14 +0530 Subject: [PATCH 11/27] used baseR findInterval() funct --- modules/benchmark/R/run_benchmark.R | 41 ++++++++++++++++------------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index ed9aa03037e..b027ea1939f 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -55,30 +55,33 @@ bm_validate <- function(model_df, obs_df) { ##' ##' @return data.frame with columns: time, model, obs align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { - # Ensure data.table is available - if (!requireNamespace("data.table", quietly = TRUE)) { - stop("Package 'data.table' is required for high-performance alignment.") - } + # 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) - # Convert to data.table and explicitly name value columns - dt_model <- data.table::data.table(time = model_df$time, model = model_df$value) - # Keep original obs time in a separate column to calculate difference after join - dt_obs <- data.table::data.table(time = obs_df$time, obs = obs_df$value, obs_time = obs_df$time) + # 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)) - # Set keys for fast roll-join - data.table::setkey(dt_model, time) - data.table::setkey(dt_obs, time) + 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"))) - # Perform rolling join to nearest time - aligned <- dt_obs[dt_model, roll = "nearest"] + # 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 tolerance - aligned[, time_diff := abs(as.numeric(difftime(obs_time, time, units = "secs")))] - aligned <- aligned[time_diff <= tolerance_secs] + # Filter by our time tolerance + valid <- time_diffs <= tolerance_secs - # Format output to match expected schema - aligned <- aligned[, .(time, model, obs)] - data.table::setDF(aligned) # Convert back to base data.frame + # 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) } From 67d3f33c35acf5074df2022903ac6aa1dd596181 Mon Sep 17 00:00:00 2001 From: ayushman1210 Date: Tue, 23 Jun 2026 17:43:40 +0530 Subject: [PATCH 12/27] Address maintainer review comments for validation framework PR - Refactor run_benchmark to use 'obvs' naming convention and point plots - Update bm_validate to explicitly check for missing time/value columns - Replace stop() and warning() with PEcAn.logger equivalents - Abstract METRIC_REGISTRY into an extensible package environment - Add NSE metric and MEF alias - Update README to reflect correct data.frame inputs - Manually export load_and_map_data and add man page - Add test coverage for missing columns, R2 metric, and tolerance drops - Add info log for time alignment tolerance dropping --- modules/benchmark/NAMESPACE | 5 ++ modules/benchmark/R/data_intake.R | 4 +- modules/benchmark/R/pecan_bench.R | 3 +- modules/benchmark/R/run_benchmark.R | 78 ++++++++++++------- modules/benchmark/README.md | 17 ++-- modules/benchmark/man/load_and_map_data.Rd | 19 +++++ .../tests/testthat/test-run_benchmark.R | 48 ++++++++++-- 7 files changed, 131 insertions(+), 43 deletions(-) create mode 100644 modules/benchmark/man/load_and_map_data.Rd diff --git a/modules/benchmark/NAMESPACE b/modules/benchmark/NAMESPACE index fc8f6f45eda..33b69823cf7 100644 --- a/modules/benchmark/NAMESPACE +++ b/modules/benchmark/NAMESPACE @@ -13,6 +13,7 @@ export(clean_settings_BRR) export(create_BRR) export(define_benchmark) export(format_wide2long) +export(load_and_map_data) export(load_csv) export(load_data) export(load_rds) @@ -34,11 +35,15 @@ export(metric_residual_plot) export(metric_run) export(metric_scatter_plot) export(metric_timeseries_plot) +export(pecan_metric_registry) export(read_settings_BRR) +export(register_metric) export(run_benchmark) +importFrom(dplyr,rename) importFrom(ggplot2,geom_path) importFrom(ggplot2,geom_point) importFrom(ggplot2,ggplot) importFrom(ggplot2,labs) importFrom(magrittr,"%>%") importFrom(rlang,.data) +importFrom(yaml,read_yaml) diff --git a/modules/benchmark/R/data_intake.R b/modules/benchmark/R/data_intake.R index 1a9e97e5178..2681001a3cc 100644 --- a/modules/benchmark/R/data_intake.R +++ b/modules/benchmark/R/data_intake.R @@ -18,7 +18,7 @@ load_and_map_data <- function(data.path, mapping.path) { mapping <- yaml::read_yaml(mapping.path) if (is.null(mapping$variables)) { - stop("YAML mapping must contain a 'variables' section.") + PEcAn.logger::logger.severe("YAML mapping must contain a 'variables' section.") } # Create a named vector for dplyr::rename (new_name = old_name) @@ -33,7 +33,7 @@ load_and_map_data <- function(data.path, mapping.path) { # Using tidy evaluation with !!! dat <- dplyr::rename(dat, !!!valid_renames) } else { - warning("No matching columns found in the dataset to map.") + PEcAn.logger::logger.warn("No matching columns found in the dataset to map.") } return(dat) diff --git a/modules/benchmark/R/pecan_bench.R b/modules/benchmark/R/pecan_bench.R index 598a25c4260..79306153658 100644 --- a/modules/benchmark/R/pecan_bench.R +++ b/modules/benchmark/R/pecan_bench.R @@ -35,8 +35,7 @@ pecan_bench <- function(comp_run, bench_id, imp_limit, high_limit) { logic_check <- validation_check(comp_run) if (isFALSE(logic_check)) { - print("The results were found to be invalid") - stop() + PEcAn.logger::logger.severe("The results were found to be invalid") } # The observed values against which the runs will be compared to. Should be the same for all diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index b027ea1939f..25d214fc8f4 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -39,10 +39,15 @@ run_benchmark <- function(model_df, obs_df, ##' @return invisible(TRUE) bm_validate <- function(model_df, obs_df) { for (df in list(model_df, obs_df)) { + if (!"time" %in% names(df)) + PEcAn.logger::logger.severe("Missing required column: 'time'") + if (!"value" %in% names(df)) + PEcAn.logger::logger.severe("Missing required column: 'value'") + if (!inherits(df$time, "POSIXct")) - stop("Column 'time' must be POSIXct, got: ", class(df$time)) + PEcAn.logger::logger.severe(paste0("Column 'time' must be POSIXct, got: ", class(df$time)[1])) if (!is.numeric(df$value)) - stop("Column 'value' must be numeric, got: ", class(df$value)) + PEcAn.logger::logger.severe(paste0("Column 'value' must be numeric, got: ", class(df$value)[1])) } invisible(TRUE) } @@ -53,7 +58,7 @@ bm_validate <- function(model_df, obs_df) { ##' @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 +##' @return data.frame with columns: model, obvs, time 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), ] @@ -76,41 +81,62 @@ align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { # Filter by our time tolerance valid <- time_diffs <= tolerance_secs + n_kept <- sum(valid) + n_dropped <- length(valid) - n_kept + PEcAn.logger::logger.info(sprintf("Time alignment kept %d points and dropped %d points outside of tolerance (%d secs)", n_kept, n_dropped, 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] + obvs = obs_df$value[nearest_idx][valid], + time = model_df$time[valid] ) return(aligned) } +##' Metric Registry for PEcAn.benchmark +##' @export +pecan_metric_registry <- new.env(parent = emptyenv()) + +##' Register a new benchmark metric +##' +##' @param name Character name of the metric +##' @param func Function that takes an aligned dataframe and returns a numeric value +##' @export +register_metric <- function(name, func) { + assign(toupper(name), func, envir = pecan_metric_registry) +} + +# Pre-populate default metrics +register_metric("RMSE", function(dat) sqrt(mean((dat$model - dat$obvs)^2, na.rm = TRUE))) +register_metric("MAE", function(dat) mean(abs(dat$model - dat$obvs), na.rm = TRUE)) +register_metric("R2", function(dat) { + if (exists("metric_R2", where = asNamespace("PEcAn.benchmark"), mode = "function")) { + return(PEcAn.benchmark::metric_R2(dat)) + } + numer <- sum((dat$obvs - mean(dat$obvs, na.rm=T)) * (dat$model - mean(dat$model, na.rm=T)), na.rm=T) + denom <- sqrt(sum((dat$obvs - mean(dat$obvs, na.rm=T))^2, na.rm=T)) * sqrt(sum((dat$model - mean(dat$model, na.rm=T))^2, na.rm=T)) + (numer / denom)^2 +}) +register_metric("NSE", function(dat) { + # Nash-Sutcliffe Efficiency + 1 - (sum((dat$obvs - dat$model)^2, na.rm = TRUE) / sum((dat$obvs - mean(dat$obvs, na.rm = TRUE))^2, na.rm = TRUE)) +}) +register_metric("MEF", get("NSE", envir = pecan_metric_registry)) + ##' Compute benchmark metrics ##' -##' @param aligned data.frame with columns: time, model, obs +##' @param aligned data.frame with columns: model, obvs, time ##' @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) + if (!exists(m, envir = pecan_metric_registry)) { + PEcAn.logger::logger.severe(paste0("Unknown metric: ", m)) + } + func <- get(m, envir = pecan_metric_registry) + func(aligned) }) data.frame(metric = toupper(metrics), value = unlist(results, use.names = FALSE)) @@ -118,12 +144,12 @@ compute_metrics <- function(aligned, metrics = c("RMSE", "MAE", "R2")) { ##' Plot model vs observations time series ##' -##' @param aligned data.frame with columns: time, model, obs +##' @param aligned data.frame with columns: model, obvs, time ##' @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::geom_point(ggplot2::aes(y = .data$obvs, color = "Obs")) + ggplot2::labs(color = "", y = "value", title = "Model vs Observations") + ggplot2::theme_bw() } diff --git a/modules/benchmark/README.md b/modules/benchmark/README.md index 69bdb5e2f3d..af6174bf921 100644 --- a/modules/benchmark/README.md +++ b/modules/benchmark/README.md @@ -1,21 +1,22 @@ ## 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. +`run_benchmark()` is a simple entry point that takes model output and +observations dataframes, 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`) +Both input dataframes must have the following two columns: +- `time` — timestamp (`POSIXct`, e.g. `2020-01-01 00:00:00`) - `value` — numeric variable value ### Usage ```r library(PEcAn.benchmark) +# Assuming model_df and obs_df are loaded data.frames res <- run_benchmark( - model_path = "inst/testdata/sample_model.csv", - obs_path = "inst/testdata/sample_obs.csv" + model_df = model_df, + obs_df = obs_df ) # View metrics @@ -30,8 +31,8 @@ res$plot ### Parameters -- `model_path` — path to model output CSV -- `obs_path` — path to observations CSV +- `model_df` — data.frame containing model output +- `obs_df` — data.frame containing observations - `metrics` — vector of metrics to compute: `"RMSE"`, `"MAE"` (default: both) - `tolerance_secs` — max time difference for matching (default: 3600 seconds) # PEcAn.benchmark diff --git a/modules/benchmark/man/load_and_map_data.Rd b/modules/benchmark/man/load_and_map_data.Rd new file mode 100644 index 00000000000..ac3249a505b --- /dev/null +++ b/modules/benchmark/man/load_and_map_data.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data_intake.R +\name{load_and_map_data} +\alias{load_and_map_data} +\title{Load and standardize arbitrary tabular data using a YAML mapping configuration} +\usage{ +load_and_map_data(data.path, mapping.path) +} +\arguments{ +\item{data.path}{character, file path to the tabular data (e.g. .csv)} + +\item{mapping.path}{character, file path to the YAML mapping configuration} +} +\value{ +A standardized data frame with column names mapped to PEcAn standard vocabulary +} +\description{ +Load and standardize arbitrary tabular data using a YAML mapping configuration +} diff --git a/modules/benchmark/tests/testthat/test-run_benchmark.R b/modules/benchmark/tests/testthat/test-run_benchmark.R index 765dfe1cb02..55c86b7402a 100644 --- a/modules/benchmark/tests/testthat/test-run_benchmark.R +++ b/modules/benchmark/tests/testthat/test-run_benchmark.R @@ -18,23 +18,61 @@ test_that("run_benchmark returns correct structure", { }) 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") + bad_type_df <- data.frame(time = c("2023-01-01"), value = c(1.0)) + expect_error(bm_validate(bad_type_df, obs_df), "POSIXct") + + missing_time_df <- data.frame(timestamp = as.POSIXct("2023-01-01"), value = c(1.0)) + expect_error(bm_validate(missing_time_df, obs_df), "Missing required column: 'time'") + + missing_val_df <- data.frame(time = as.POSIXct("2023-01-01"), val = c(1.0)) + expect_error(bm_validate(missing_val_df, obs_df), "Missing required column: 'value'") }) 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) + obvs = c(1, 2, 3, 4) ) - res <- compute_metrics(aligned, c("RMSE", "MAE")) + res <- compute_metrics(aligned, c("RMSE", "MAE", "NSE", "R2")) expect_equal(res$value[res$metric == "RMSE"], 0) expect_equal(res$value[res$metric == "MAE"], 0) + expect_equal(res$value[res$metric == "NSE"], 1) + expect_equal(res$value[res$metric == "R2"], 1) +}) + +test_that("register_metric extends the registry", { + aligned <- data.frame( + time = as.POSIXct("2023-01-01", tz="UTC"), + model = 1, + obvs = 1 + ) + register_metric("CUSTOM", function(dat) 999) + res <- compute_metrics(aligned, c("CUSTOM")) + expect_equal(res$value[res$metric == "CUSTOM"], 999) }) 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))) + expect_true(all(c("model", "obvs", "time") %in% names(aligned))) + expect_equal(aligned$model, c(1, 2, 3, 4)) + expect_equal(aligned$obvs, c(1.1, 1.9, 3.2, 3.9)) +}) + +test_that("align_by_time drops points outside tolerance", { + model_df_tol <- data.frame( + time = as.POSIXct(c(0, 3600), origin = "1970-01-01", tz = "UTC"), + value = c(1, 2) + ) + obs_df_tol <- data.frame( + time = as.POSIXct(c(10, 5000), origin = "1970-01-01", tz = "UTC"), + value = c(1.1, 2.1) + ) + + # Tolerance of 20 seconds should keep the first point (diff=10s) but drop the second (diff=1400s) + aligned <- align_by_time(model_df_tol, obs_df_tol, tolerance_secs = 20) + expect_equal(nrow(aligned), 1) + expect_equal(aligned$model, 1) + expect_equal(aligned$obvs, 1.1) }) From 4eaf2219e6d1545b179039aaea7d2432e1036ed2 Mon Sep 17 00:00:00 2001 From: ayushman1210 Date: Tue, 23 Jun 2026 18:31:44 +0530 Subject: [PATCH 13/27] Expand dataframe passthrough and add PMU/Coverage metrics --- modules/benchmark/NAMESPACE | 2 + modules/benchmark/R/metric_Coverage.R | 21 ++++++++++ modules/benchmark/R/metric_PMU.R | 25 ++++++++++++ modules/benchmark/R/run_benchmark.R | 35 ++++++++++++++--- .../tests/testthat/test-run_benchmark.R | 38 +++++++++++++++++++ 5 files changed, 115 insertions(+), 6 deletions(-) create mode 100644 modules/benchmark/R/metric_Coverage.R create mode 100644 modules/benchmark/R/metric_PMU.R diff --git a/modules/benchmark/NAMESPACE b/modules/benchmark/NAMESPACE index 33b69823cf7..1cbb0236d2e 100644 --- a/modules/benchmark/NAMESPACE +++ b/modules/benchmark/NAMESPACE @@ -30,6 +30,8 @@ export(metric_R2) export(metric_RAE) export(metric_RMSE) export(metric_cor) +export(metric_PMU) +export(metric_Coverage) export(metric_lmDiag_plot) export(metric_residual_plot) export(metric_run) diff --git a/modules/benchmark/R/metric_Coverage.R b/modules/benchmark/R/metric_Coverage.R new file mode 100644 index 00000000000..52302dce834 --- /dev/null +++ b/modules/benchmark/R/metric_Coverage.R @@ -0,0 +1,21 @@ +##' @name metric_Coverage +##' @title Prediction Interval Coverage +##' @export +##' @param dat dataframe with columns `model_q025` and `model_q975` +##' @param ... ignored +##' @details +##' Measures the fraction of observations that fall within the model's +##' stated 95% prediction interval. + +metric_Coverage <- function(dat, ...) { + if (!"model_q025" %in% names(dat) || !"model_q975" %in% names(dat)) { + PEcAn.logger::logger.severe("Metric Coverage requires 'model_q025' and 'model_q975' columns in the dataset.") + } + + PEcAn.logger::logger.info("Metric: Prediction Interval Coverage") + + valid <- !is.na(dat$obvs) & !is.na(dat$model_q025) & !is.na(dat$model_q975) + covered <- dat$obvs[valid] >= dat$model_q025[valid] & dat$obvs[valid] <= dat$model_q975[valid] + + return(mean(covered)) +} diff --git a/modules/benchmark/R/metric_PMU.R b/modules/benchmark/R/metric_PMU.R new file mode 100644 index 00000000000..f91a87f2e77 --- /dev/null +++ b/modules/benchmark/R/metric_PMU.R @@ -0,0 +1,25 @@ +##' @name metric_PMU +##' @title Predictive Model Uncertainty (PMU) +##' @export +##' @param dat dataframe with columns `obs_se` and `obs_n` +##' @param ... ignored +##' @details +##' Calculates the pooled measurement uncertainty. +##' Requires `obs_se` (standard error) and `obs_n` (replicate counts) in the observation data. + +metric_PMU <- function(dat, ...) { + if (!"obs_se" %in% names(dat) || !"obs_n" %in% names(dat)) { + PEcAn.logger::logger.severe("Metric PMU requires 'obs_se' and 'obs_n' columns in the dataset.") + } + + PEcAn.logger::logger.info("Metric: Predictive Model Uncertainty (PMU)") + + # Calculate pooled standard error + # Formula: sqrt(sum(SE^2 * n) / sum(n)) + valid <- !is.na(dat$obs_se) & !is.na(dat$obs_n) + + se2_n <- (dat$obs_se[valid]^2) * dat$obs_n[valid] + pooled_var <- sum(se2_n) / sum(dat$obs_n[valid]) + + return(sqrt(pooled_var)) +} diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index 25d214fc8f4..407bfb70d15 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -85,12 +85,23 @@ align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { n_dropped <- length(valid) - n_kept PEcAn.logger::logger.info(sprintf("Time alignment kept %d points and dropped %d points outside of tolerance (%d secs)", n_kept, n_dropped, tolerance_secs)) - # Construct the aligned base data.frame - aligned <- data.frame( - model = model_df$value[valid], - obvs = obs_df$value[nearest_idx][valid], - time = model_df$time[valid] - ) + # Construct the aligned base data.frame without dropping metadata + # Rename value columns to prevent collision and fit convention + names(model_df)[names(model_df) == "value"] <- "model" + names(obs_df)[names(obs_df) == "value"] <- "obvs" + + # Prevent time collision if obs_df carries it forward + if ("time" %in% names(obs_df)) { + names(obs_df)[names(obs_df) == "time"] <- "obs_time" + } + + model_sub <- model_df[valid, , drop = FALSE] + obs_sub <- obs_df[nearest_idx[valid], , drop = FALSE] + + # Drop overlapping columns from obs to cleanly cbind + obs_sub <- obs_sub[, !(names(obs_sub) %in% names(model_sub)), drop = FALSE] + + aligned <- cbind(model_sub, obs_sub) return(aligned) } @@ -124,6 +135,18 @@ register_metric("NSE", function(dat) { 1 - (sum((dat$obvs - dat$model)^2, na.rm = TRUE) / sum((dat$obvs - mean(dat$obvs, na.rm = TRUE))^2, na.rm = TRUE)) }) register_metric("MEF", get("NSE", envir = pecan_metric_registry)) +register_metric("PMU", function(dat) { + if (exists("metric_PMU", where = asNamespace("PEcAn.benchmark"), mode = "function")) { + return(PEcAn.benchmark::metric_PMU(dat)) + } + metric_PMU(dat) +}) +register_metric("COVERAGE", function(dat) { + if (exists("metric_Coverage", where = asNamespace("PEcAn.benchmark"), mode = "function")) { + return(PEcAn.benchmark::metric_Coverage(dat)) + } + metric_Coverage(dat) +}) ##' Compute benchmark metrics ##' diff --git a/modules/benchmark/tests/testthat/test-run_benchmark.R b/modules/benchmark/tests/testthat/test-run_benchmark.R index 55c86b7402a..17f95ee9c79 100644 --- a/modules/benchmark/tests/testthat/test-run_benchmark.R +++ b/modules/benchmark/tests/testthat/test-run_benchmark.R @@ -76,3 +76,41 @@ test_that("align_by_time drops points outside tolerance", { expect_equal(aligned$model, 1) expect_equal(aligned$obvs, 1.1) }) + +test_that("align_by_time passes through metadata columns", { + model_df_meta <- data.frame( + time = as.POSIXct(c(0, 3600), origin = "1970-01-01", tz = "UTC"), + value = c(1, 2), + model_q025 = c(0.5, 1.5), + model_q975 = c(1.5, 2.5) + ) + obs_df_meta <- data.frame( + time = as.POSIXct(c(0, 3600), origin = "1970-01-01", tz = "UTC"), + value = c(1.1, 1.9), + obs_se = c(0.1, 0.2), + obs_n = c(3, 3) + ) + + aligned <- align_by_time(model_df_meta, obs_df_meta) + expect_true(all(c("model_q025", "model_q975", "obs_se", "obs_n") %in% names(aligned))) +}) + +test_that("metric_Coverage calculates correct fraction", { + aligned <- data.frame( + model = c(1, 2, 3), + obvs = c(1, 4, 3), + model_q025 = c(0.5, 1.5, 2.5), + model_q975 = c(1.5, 2.5, 3.5) + ) + + expect_equal(metric_Coverage(aligned), 2/3) +}) + +test_that("metric_PMU calculates correct pooled uncertainty", { + aligned <- data.frame( + obs_se = c(0.1, 0.2), + obs_n = c(3, 5) + ) + + expect_equal(metric_PMU(aligned), sqrt(0.23/8)) +}) From f7e0c9950f0dbae43f816a34eb7aef505a2b98a4 Mon Sep 17 00:00:00 2001 From: ayushman1210 Date: Thu, 25 Jun 2026 21:34:05 +0530 Subject: [PATCH 14/27] Refactor Phase 2: Standardize Roxygen docs and integrate NetCDF loader with YAML intake --- modules/benchmark/R/data_intake.R | 16 +++-- modules/benchmark/R/load_netcdf.R | 61 ++++++++++--------- modules/benchmark/R/metric_Coverage.R | 17 +++--- modules/benchmark/R/metric_PMU.R | 17 +++--- modules/benchmark/R/run_benchmark.R | 84 +++++++++++++-------------- 5 files changed, 105 insertions(+), 90 deletions(-) diff --git a/modules/benchmark/R/data_intake.R b/modules/benchmark/R/data_intake.R index 2681001a3cc..9a061d26bc7 100644 --- a/modules/benchmark/R/data_intake.R +++ b/modules/benchmark/R/data_intake.R @@ -7,10 +7,7 @@ #' @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 + # Load the YAML mapping first to know which variables we need # The YAML should look like: # variables: # airT: TA_F @@ -23,6 +20,17 @@ load_and_map_data <- function(data.path, mapping.path) { # Create a named vector for dplyr::rename (new_name = old_name) rename_vector <- unlist(mapping$variables) + required_vars <- unname(rename_vector) + + # Load the raw data based on file extension + if (grepl("\\.nc$", data.path, ignore.case = TRUE)) { + # If NetCDF, we only load the variables requested in the YAML mapping + # Assuming standard NA string representations for now, can be expanded via YAML + dat <- load_x_netcdf(data.path, format = list(na.strings = c("-9999", "-9999.0", "NA")), site = NULL, vars = required_vars) + } else { + # Default to CSV + dat <- utils::read.csv(data.path, as.is = TRUE, check.names = FALSE) + } # Only rename columns that exist in the raw data valid_renames <- rename_vector[rename_vector %in% colnames(dat)] diff --git a/modules/benchmark/R/load_netcdf.R b/modules/benchmark/R/load_netcdf.R index 2e956ffee85..7ac9e20dc0d 100644 --- a/modules/benchmark/R/load_netcdf.R +++ b/modules/benchmark/R/load_netcdf.R @@ -1,11 +1,12 @@ -##' Load from netCDF -##' -##' @param data.path character vector or list -##' @param format list -##' @param site list -##' @param vars character -##' @author Istem Fer -##' @export +#' Load from netCDF +#' +#' @param data.path character vector or list of paths to NetCDF files +#' @param format list containing format info (e.g. na.strings) +#' @param site list site info (optional) +#' @param vars character vector of variables to load +#' @return A data.frame with the requested variables and a `time` column (POSIXct) +#' @author Istem Fer, Anshul Jain +#' @export load_x_netcdf <- function(data.path, format, site, vars = NULL) { data.path <- sapply(data.path, function(x) dir(dirname(x), basename(x), full.names = TRUE)) nc <- lapply(data.path, ncdf4::nc_open) @@ -17,7 +18,9 @@ load_x_netcdf <- function(data.path, format, site, vars = NULL) { dat <- as.matrix(as.data.frame(dat)) # we need to replace filling/missing values with NA now we don't want these values to go into unit # conversion - dat[dat %in% as.numeric(format$na.strings)] <- NA + if (!is.null(format$na.strings)) { + dat[dat %in% as.numeric(format$na.strings)] <- NA + } dat <- as.data.frame(dat) colnames(dat) <- vars # deal with time @@ -25,47 +28,49 @@ load_x_netcdf <- function(data.path, format, site, vars = NULL) { for (i in seq_along(nc)) { dims <- names(nc[[i]]$dim) time.var <- grep(pattern = "time", dims, ignore.case = TRUE) - time.col[[i]] <- ncdf4::ncvar_get(nc[[i]], dims[time.var]) - t.units <- ncdf4::ncatt_get(nc[[i]], dims[time.var])$units + if (length(time.var) == 0) { + PEcAn.logger::logger.error("No 'time' dimension found in NetCDF file.") + } + time.col[[i]] <- ncdf4::ncvar_get(nc[[i]], dims[time.var[1]]) + t.units <- ncdf4::ncatt_get(nc[[i]], dims[time.var[1]])$units # If the unit has if of the form * since YYYY-MM-DD * with "-hour" timezone offset # This is a feature of the met produced by met2CF if(stringr::str_detect(t.units, "ince\\s[0-9]{4}[.-][0-9]{2}[.-][0-9]{2}.*\\s-\\d+")){ unit2 <- stringr::str_split_fixed(t.units,"\\s-",2)[1] offset <- stringr::str_split_fixed(t.units,"\\s-",2)[2] %>% as.numeric() - date_time <- suppressWarnings(try(lubridate::ymd((unit2)))) - if(is.na(date_time)){ - date_time <- suppressWarnings(try(lubridate::ymd_hms(unit2))) + date_time <- suppressWarnings(try(lubridate::ymd((unit2)), silent = TRUE)) + if(inherits(date_time, "try-error") || is.na(date_time)){ + date_time <- suppressWarnings(try(lubridate::ymd_hms(unit2), silent = TRUE)) } - if(is.na(date_time)){ + if(inherits(date_time, "try-error") || is.na(date_time)){ PEcAn.logger::logger.error("All time formats failed to parse. No formats found.") } t.units <- paste(stringr::str_split_fixed(t.units," since",2)[1], "since", date_time - lubridate::hms(paste(offset,":00:00"))) }else if(stringr::str_detect(t.units, "ince\\s[0-9]{4}[.-][0-9]{2}[.-][0-9]{2}.*")){ unit2 <- stringr::str_split_fixed(t.units,"\\s-",2)[1] - date_time <- suppressWarnings(try(lubridate::ymd((unit2)))) - if(is.na(date_time)){ - date_time <- suppressWarnings(try(lubridate::ymd_hms(unit2))) + date_time <- suppressWarnings(try(lubridate::ymd((unit2)), silent = TRUE)) + if(inherits(date_time, "try-error") || is.na(date_time)){ + date_time <- suppressWarnings(try(lubridate::ymd_hms(unit2), silent = TRUE)) } - if(is.na(date_time)){ + if(inherits(date_time, "try-error") || is.na(date_time)){ PEcAn.logger::logger.error("All time formats failed to parse. No formats found.") } t.units <- paste(stringr::str_split_fixed(t.units," since",2)[1], "since", date_time) } # for heterogenous formats try parsing ymd_hms - date.origin <- suppressWarnings(try(lubridate::ymd_hms(t.units))) + date.origin <- suppressWarnings(try(lubridate::ymd_hms(t.units), silent = TRUE)) # parsing ymd - if (is.na(date.origin)) { - date.origin <- lubridate::ymd(t.units) + if (inherits(date.origin, "try-error") || is.na(date.origin)) { + date.origin <- suppressWarnings(try(lubridate::ymd(t.units), silent = TRUE)) } # throw error if can't parse time format - if (is.na(date.origin)) { + if (inherits(date.origin, "try-error") || is.na(date.origin)) { PEcAn.logger::logger.error("All time formats failed to parse. No formats found.") } time.stamp.match <- gsub("UTC", "", date.origin) - t.units <- gsub(paste0(" since ", time.stamp.match, ".*"), "", - t.units) + t.units <- gsub(paste0(" since ", time.stamp.match, ".*"), "", t.units) # need to change system TZ otherwise, lines below keeps writing in the current time zone Sys.setenv(TZ = 'UTC') foo <- as.POSIXct(date.origin, tz = "UTC") + PEcAn.utils::ud_convert(time.col[[i]], t.units, "seconds") @@ -75,8 +80,8 @@ load_x_netcdf <- function(data.path, format, site, vars = NULL) { # while reading Ameriflux for example however the model timesteps are more regular and the last # value can be '2006-12-31 23:30:00'.. this will result in cutting the last value in the # align_data step - dat$posix <- round(as.POSIXct(do.call("c", time.col), tz = "UTC"), "mins") - dat$posix <- as.POSIXct(dat$posix) + dat$time <- round(as.POSIXct(do.call("c", time.col), tz = "UTC"), "mins") + dat$time <- as.POSIXct(dat$time) lapply(nc, ncdf4::nc_close) return(dat) -} # load_x_netcdf +} diff --git a/modules/benchmark/R/metric_Coverage.R b/modules/benchmark/R/metric_Coverage.R index 52302dce834..d5943ad5ad1 100644 --- a/modules/benchmark/R/metric_Coverage.R +++ b/modules/benchmark/R/metric_Coverage.R @@ -1,11 +1,12 @@ -##' @name metric_Coverage -##' @title Prediction Interval Coverage -##' @export -##' @param dat dataframe with columns `model_q025` and `model_q975` -##' @param ... ignored -##' @details -##' Measures the fraction of observations that fall within the model's -##' stated 95% prediction interval. +#' @name metric_Coverage +#' @title Prediction Interval Coverage +#' @export +#' @param dat dataframe with columns `model_q025` and `model_q975` +#' @param ... ignored +#' @return A numeric value representing the fraction of observations that fall within the 95% prediction interval. +#' @details +#' Measures the fraction of observations that fall within the model's +#' stated 95% prediction interval. metric_Coverage <- function(dat, ...) { if (!"model_q025" %in% names(dat) || !"model_q975" %in% names(dat)) { diff --git a/modules/benchmark/R/metric_PMU.R b/modules/benchmark/R/metric_PMU.R index f91a87f2e77..1ab473a992b 100644 --- a/modules/benchmark/R/metric_PMU.R +++ b/modules/benchmark/R/metric_PMU.R @@ -1,11 +1,12 @@ -##' @name metric_PMU -##' @title Predictive Model Uncertainty (PMU) -##' @export -##' @param dat dataframe with columns `obs_se` and `obs_n` -##' @param ... ignored -##' @details -##' Calculates the pooled measurement uncertainty. -##' Requires `obs_se` (standard error) and `obs_n` (replicate counts) in the observation data. +#' @name metric_PMU +#' @title Predictive Model Uncertainty (PMU) +#' @export +#' @param dat dataframe with columns `obs_se` and `obs_n` +#' @param ... ignored +#' @return A numeric value representing the pooled measurement uncertainty. +#' @details +#' Calculates the pooled measurement uncertainty. +#' Requires `obs_se` (standard error) and `obs_n` (replicate counts) in the observation data. metric_PMU <- function(dat, ...) { if (!"obs_se" %in% names(dat) || !"obs_n" %in% names(dat)) { diff --git a/modules/benchmark/R/run_benchmark.R b/modules/benchmark/R/run_benchmark.R index 407bfb70d15..42b7faa415b 100644 --- a/modules/benchmark/R/run_benchmark.R +++ b/modules/benchmark/R/run_benchmark.R @@ -1,17 +1,17 @@ -##' 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 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, @@ -32,11 +32,11 @@ run_benchmark <- function(model_df, obs_df, 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) +#' 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 (!"time" %in% names(df)) @@ -52,13 +52,13 @@ bm_validate <- function(model_df, obs_df) { 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: model, obvs, time +#' 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: model, obvs, time 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), ] @@ -106,15 +106,15 @@ align_by_time <- function(model_df, obs_df, tolerance_secs = 3600) { return(aligned) } -##' Metric Registry for PEcAn.benchmark -##' @export +#' Metric Registry for PEcAn.benchmark +#' @export pecan_metric_registry <- new.env(parent = emptyenv()) -##' Register a new benchmark metric -##' -##' @param name Character name of the metric -##' @param func Function that takes an aligned dataframe and returns a numeric value -##' @export +#' Register a new benchmark metric +#' +#' @param name Character name of the metric +#' @param func Function that takes an aligned dataframe and returns a numeric value +#' @export register_metric <- function(name, func) { assign(toupper(name), func, envir = pecan_metric_registry) } @@ -148,11 +148,11 @@ register_metric("COVERAGE", function(dat) { metric_Coverage(dat) }) -##' Compute benchmark metrics -##' -##' @param aligned data.frame with columns: model, obvs, time -##' @param metrics character vector of metric names -##' @return data.frame with columns: metric, value +#' Compute benchmark metrics +#' +#' @param aligned data.frame with columns: model, obvs, time +#' @param metrics character vector of metric names +#' @return data.frame with columns: metric, value compute_metrics <- function(aligned, metrics = c("RMSE", "MAE", "R2")) { results <- lapply(toupper(metrics), function(m) { if (!exists(m, envir = pecan_metric_registry)) { @@ -165,10 +165,10 @@ compute_metrics <- function(aligned, metrics = c("RMSE", "MAE", "R2")) { data.frame(metric = toupper(metrics), value = unlist(results, use.names = FALSE)) } -##' Plot model vs observations time series -##' -##' @param aligned data.frame with columns: model, obvs, time -##' @return ggplot object +#' Plot model vs observations time series +#' +#' @param aligned data.frame with columns: model, obvs, time +#' @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")) + From 014fb3e2ec71dfdef5fca44cf6c236749d9a9864 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:43:05 -0400 Subject: [PATCH 15/27] update namespace --- modules/benchmark/NAMESPACE | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/benchmark/NAMESPACE b/modules/benchmark/NAMESPACE index 1cbb0236d2e..07067c6c4c0 100644 --- a/modules/benchmark/NAMESPACE +++ b/modules/benchmark/NAMESPACE @@ -22,16 +22,16 @@ export(load_x_netcdf) export(match_timestep) export(mean_over_larger_timestep) export(metric_AME) +export(metric_Coverage) export(metric_Frechet) export(metric_MAE) export(metric_MSE) +export(metric_PMU) export(metric_PPMC) export(metric_R2) export(metric_RAE) export(metric_RMSE) export(metric_cor) -export(metric_PMU) -export(metric_Coverage) export(metric_lmDiag_plot) export(metric_residual_plot) export(metric_run) From e1edfbe0975481ed51b59775b26b21f9b87f3c0f Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:43:26 -0400 Subject: [PATCH 16/27] update align_by_time.Rd --- modules/benchmark/man/align_by_time.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/benchmark/man/align_by_time.Rd b/modules/benchmark/man/align_by_time.Rd index 33e313bb5b5..02b9032b24c 100644 --- a/modules/benchmark/man/align_by_time.Rd +++ b/modules/benchmark/man/align_by_time.Rd @@ -14,7 +14,7 @@ align_by_time(model_df, obs_df, tolerance_secs = 3600) \item{tolerance_secs}{max allowed time difference in seconds} } \value{ -data.frame with columns: time, model, obs +data.frame with columns: model, obvs, time } \description{ Align model and observation data frames by nearest time From 09ac45a78fe5b49c4ba82a100c52443fcbd75e1b Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:44:13 -0400 Subject: [PATCH 17/27] update compute_metrics.Rd --- modules/benchmark/man/compute_metrics.Rd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/benchmark/man/compute_metrics.Rd b/modules/benchmark/man/compute_metrics.Rd index a9a66fe98e9..41a8d74a916 100644 --- a/modules/benchmark/man/compute_metrics.Rd +++ b/modules/benchmark/man/compute_metrics.Rd @@ -4,10 +4,10 @@ \alias{compute_metrics} \title{Compute benchmark metrics} \usage{ -compute_metrics(aligned, metrics = c("RMSE", "MAE")) +compute_metrics(aligned, metrics = c("RMSE", "MAE", "R2")) } \arguments{ -\item{aligned}{data.frame with columns: time, model, obs} +\item{aligned}{data.frame with columns: model, obvs, time} \item{metrics}{character vector of metric names} } From 00e9c7d92c7c3946569a2582adfce5007f1d2ce8 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:44:40 -0400 Subject: [PATCH 18/27] update load_x_netcdf.Rd --- modules/benchmark/man/load_x_netcdf.Rd | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/modules/benchmark/man/load_x_netcdf.Rd b/modules/benchmark/man/load_x_netcdf.Rd index c02d2b55229..fcb1afc6f23 100644 --- a/modules/benchmark/man/load_x_netcdf.Rd +++ b/modules/benchmark/man/load_x_netcdf.Rd @@ -7,17 +7,20 @@ load_x_netcdf(data.path, format, site, vars = NULL) } \arguments{ -\item{data.path}{character vector or list} +\item{data.path}{character vector or list of paths to NetCDF files} -\item{format}{list} +\item{format}{list containing format info (e.g. na.strings)} -\item{site}{list} +\item{site}{list site info (optional)} -\item{vars}{character} +\item{vars}{character vector of variables to load} +} +\value{ +A data.frame with the requested variables and a `time` column (POSIXct) } \description{ Load from netCDF } \author{ -Istem Fer +Istem Fer, Anshul Jain } From d97792d40dc28fbb7a7a070f13b35f3d6fbaa8cf Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:45:07 -0400 Subject: [PATCH 19/27] plot_time_series.Rd --- modules/benchmark/man/plot_time_series.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/benchmark/man/plot_time_series.Rd b/modules/benchmark/man/plot_time_series.Rd index 0710201ed1a..51aae9201f4 100644 --- a/modules/benchmark/man/plot_time_series.Rd +++ b/modules/benchmark/man/plot_time_series.Rd @@ -7,7 +7,7 @@ plot_time_series(aligned) } \arguments{ -\item{aligned}{data.frame with columns: time, model, obs} +\item{aligned}{data.frame with columns: model, obvs, time} } \value{ ggplot object From 1cbeec9ce72ec2f5047ed1456a58c12a981623a2 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:45:34 -0400 Subject: [PATCH 20/27] add metric_Coverage.Rd --- modules/benchmark/man/metric_Coverage.Rd | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 modules/benchmark/man/metric_Coverage.Rd diff --git a/modules/benchmark/man/metric_Coverage.Rd b/modules/benchmark/man/metric_Coverage.Rd new file mode 100644 index 00000000000..d4f891b9e49 --- /dev/null +++ b/modules/benchmark/man/metric_Coverage.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metric_Coverage.R +\name{metric_Coverage} +\alias{metric_Coverage} +\title{Prediction Interval Coverage} +\usage{ +metric_Coverage(dat, ...) +} +\arguments{ +\item{dat}{dataframe with columns `model_q025` and `model_q975`} + +\item{...}{ignored} +} +\value{ +A numeric value representing the fraction of observations that fall within the 95% prediction interval. +} +\description{ +Prediction Interval Coverage +} +\details{ +Measures the fraction of observations that fall within the model's +stated 95% prediction interval. +} From 98d7eac4e5baded2e6effdae8351e8bfd30aec08 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:45:52 -0400 Subject: [PATCH 21/27] metric_PMU.Rd --- modules/benchmark/man/metric_PMU.Rd | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 modules/benchmark/man/metric_PMU.Rd diff --git a/modules/benchmark/man/metric_PMU.Rd b/modules/benchmark/man/metric_PMU.Rd new file mode 100644 index 00000000000..c28330a1e9f --- /dev/null +++ b/modules/benchmark/man/metric_PMU.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metric_PMU.R +\name{metric_PMU} +\alias{metric_PMU} +\title{Predictive Model Uncertainty (PMU)} +\usage{ +metric_PMU(dat, ...) +} +\arguments{ +\item{dat}{dataframe with columns `obs_se` and `obs_n`} + +\item{...}{ignored} +} +\value{ +A numeric value representing the pooled measurement uncertainty. +} +\description{ +Predictive Model Uncertainty (PMU) +} +\details{ +Calculates the pooled measurement uncertainty. +Requires `obs_se` (standard error) and `obs_n` (replicate counts) in the observation data. +} From e627d97d232a8a11abc02d598323880194c484d1 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:46:22 -0400 Subject: [PATCH 22/27] add pecan_metric_registry.Rd --- modules/benchmark/man/pecan_metric_registry.Rd | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 modules/benchmark/man/pecan_metric_registry.Rd diff --git a/modules/benchmark/man/pecan_metric_registry.Rd b/modules/benchmark/man/pecan_metric_registry.Rd new file mode 100644 index 00000000000..30701d12129 --- /dev/null +++ b/modules/benchmark/man/pecan_metric_registry.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R +\docType{data} +\name{pecan_metric_registry} +\alias{pecan_metric_registry} +\title{Metric Registry for PEcAn.benchmark} +\format{ +An object of class \code{environment} of length 7. +} +\usage{ +pecan_metric_registry +} +\description{ +Metric Registry for PEcAn.benchmark +} +\keyword{datasets} From 21d749fc9de8966c750895bafbf7d972daeb4260 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:46:42 -0400 Subject: [PATCH 23/27] add register_metric.Rd --- modules/benchmark/man/register_metric.Rd | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 modules/benchmark/man/register_metric.Rd diff --git a/modules/benchmark/man/register_metric.Rd b/modules/benchmark/man/register_metric.Rd new file mode 100644 index 00000000000..ed9db7abdc4 --- /dev/null +++ b/modules/benchmark/man/register_metric.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/run_benchmark.R +\name{register_metric} +\alias{register_metric} +\title{Register a new benchmark metric} +\usage{ +register_metric(name, func) +} +\arguments{ +\item{name}{Character name of the metric} + +\item{func}{Function that takes an aligned dataframe and returns a numeric value} +} +\description{ +Register a new benchmark metric +} From 756ec13a6e305709f653531dd54221d15e7f68db Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 04:47:06 -0400 Subject: [PATCH 24/27] update load_x_netcdf.Rd --- modules/benchmark/man/load_x_netcdf.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/benchmark/man/load_x_netcdf.Rd b/modules/benchmark/man/load_x_netcdf.Rd index fcb1afc6f23..3fdc9fc44ae 100644 --- a/modules/benchmark/man/load_x_netcdf.Rd +++ b/modules/benchmark/man/load_x_netcdf.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2: do not edit by hand +load_x_netcdf.R% Generated by roxygen2: do not edit by hand % Please edit documentation in R/load_netcdf.R \name{load_x_netcdf} \alias{load_x_netcdf} From 0df39d795883f75ce062449ab2a8967e7a6665de Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 06:34:51 -0400 Subject: [PATCH 25/27] update docker depends --- docker/depends/pecan_package_dependencies.csv | 1 + 1 file changed, 1 insertion(+) diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv index 75142c527f8..7ee8e35d125 100644 --- a/docker/depends/pecan_package_dependencies.csv +++ b/docker/depends/pecan_package_dependencies.csv @@ -749,5 +749,6 @@ "xtable","*","base/utils","Suggests",FALSE "xts","*","models/peprmt","Suggests",FALSE "xts","*","modules/data.atmosphere","Imports",FALSE +"yaml","*","modules/benchmark","Imports",FALSE "zoo","*","modules/benchmark","Imports",FALSE "zoo","*","modules/data.atmosphere","Imports",FALSE From dc4eecdb6f14d0f56ea300445e00b3475954b482 Mon Sep 17 00:00:00 2001 From: divne7022 Date: Mon, 29 Jun 2026 06:43:02 -0400 Subject: [PATCH 26/27] add yaml to desc --- modules/benchmark/DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/modules/benchmark/DESCRIPTION b/modules/benchmark/DESCRIPTION index f200f140925..f8ada65d1d7 100644 --- a/modules/benchmark/DESCRIPTION +++ b/modules/benchmark/DESCRIPTION @@ -43,7 +43,8 @@ Imports: utils, grDevices, XML (>= 3.98-1.4), - zoo + zoo, + yaml Suggests: PEcAn.data.land, testthat (>= 2.0.0) From 656e29846193e8779289fe7a64b8c6739f4b5f3f Mon Sep 17 00:00:00 2001 From: divne7022 Date: Tue, 30 Jun 2026 08:47:26 -0400 Subject: [PATCH 27/27] remove hand edited load_x_netcdf.Rd and auto gen --- modules/benchmark/man/load_x_netcdf.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/benchmark/man/load_x_netcdf.Rd b/modules/benchmark/man/load_x_netcdf.Rd index 3fdc9fc44ae..fcb1afc6f23 100644 --- a/modules/benchmark/man/load_x_netcdf.Rd +++ b/modules/benchmark/man/load_x_netcdf.Rd @@ -1,4 +1,4 @@ -load_x_netcdf.R% Generated by roxygen2: do not edit by hand +% Generated by roxygen2: do not edit by hand % Please edit documentation in R/load_netcdf.R \name{load_x_netcdf} \alias{load_x_netcdf}