Refactor the baseline Validation Pipeline to align with GSoC Architecture

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

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)

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)
@@ -21,9 +22,11 @@ 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)
@@ -34,10 +37,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)

modules/benchmark/R/data_intake.R

@@ -0,0 +1,48 @@
+#' 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 YAML mapping first to know which variables we need
+  # The YAML should look like:
+  # variables:
+  #   airT: TA_F
+  #   NEE: NEE_PI
+  mapping <- yaml::read_yaml(mapping.path)
+
+  if (is.null(mapping$variables)) {
+    PEcAn.logger::logger.severe("YAML mapping must contain a 'variables' section.")
+  }
+
+  # 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)]
+
+  # 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 {
+    PEcAn.logger::logger.warn("No matching columns found in the dataset to map.")
+  }
+
+  return(dat)
+}

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,55 +18,59 @@ 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
   time.col <- list()
   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
+}

modules/benchmark/R/metric_Coverage.R

@@ -0,0 +1,22 @@
+#' @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)) {
+    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))
+}

modules/benchmark/R/metric_PMU.R

@@ -0,0 +1,26 @@
+#' @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)) {
+    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))
+}

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