diff --git a/modules/uncertainty/NEWS.md b/modules/uncertainty/NEWS.md index ac0c0ee6484..869f9780b03 100644 --- a/modules/uncertainty/NEWS.md +++ b/modules/uncertainty/NEWS.md @@ -1,5 +1,10 @@ # PEcAn.uncertainty 1.9.0.9000 +## Fixed + +* `flux.uncertainty()`: removed a duplicate assignment of `E2` (the same line appeared twice in sequence with no effect except wasted computation). +* `flux.uncertainty()`: fixed a crash that occurred when all positive-flux magnitude bins contained `NA` error estimates. In that case the positive-slope model `mp` was never assigned, but the `else` branch for the negative-slope fallback unconditionally referenced `mp$coefficients[1]`, producing an "object 'mp' not found" error. The fallback now uses `else if (exists("mp", inherits = FALSE))` so it only runs when `mp` was actually fitted. + * Multiple bugfixes in `input.ens.gen()` handling of parent ids (#3783): - No longer skips inputs that have a parent but no sampling method. - Argument `parent_ids` now accepts integer vectors even if not wrapped in a list. diff --git a/modules/uncertainty/R/flux_uncertainty.R b/modules/uncertainty/R/flux_uncertainty.R index ef18e1ad31d..65bf39e6164 100644 --- a/modules/uncertainty/R/flux_uncertainty.R +++ b/modules/uncertainty/R/flux_uncertainty.R @@ -105,23 +105,22 @@ flux.uncertainty <- function(measurement, QC = 0, flags = TRUE, bin.num = 10, ## would be better to fit a two line model with a common intercept, but this ## is quicker to implement for the time being E2 <- errBin - errBin[zero] - E2 <- errBin - errBin[zero] intercept <- errBin[zero] - - return.list <- list(mag = magBin, - err = errBin, - bias = biasBin, + + return.list <- list(mag = magBin, + err = errBin, + bias = biasBin, n = nBin, intercept = intercept) - - if(!all(is.na(E2[pos]))){ + + if (!all(is.na(E2[pos]))) { mp <- stats::lm(E2[pos] ~ magBin[pos] - 1) return.list$slopeP <- mp$coefficients[1] - } - if(!all(is.na(E2[neg]))){ + } + if (!all(is.na(E2[neg]))) { mn <- stats::lm(E2[neg] ~ magBin[neg] - 1) return.list$slopeN <- mn$coefficients[1] - }else{ + } else if (exists("mp", inherits = FALSE)) { return.list$slopeN <- mp$coefficients[1] } diff --git a/modules/uncertainty/tests/testthat/test.flux_uncertainty.R b/modules/uncertainty/tests/testthat/test.flux_uncertainty.R new file mode 100644 index 00000000000..d8f96965434 --- /dev/null +++ b/modules/uncertainty/tests/testthat/test.flux_uncertainty.R @@ -0,0 +1,78 @@ +## Tests for flux.uncertainty() +## +## The bug fixed in this PR: +## When all flux measurements are negative (all bins fall below zero), +## the positive-flux linear model `mp` is never assigned. The original +## code unconditionally referenced `mp` in the else-branch for slopeN, +## causing an "object 'mp' not found" error at runtime. +## +## Fix: the else-branch is now `else if (exists("mp", inherits = FALSE))` +## so slopeN falls back to the positive model only when it was actually fit. + +## Helper: create paired flux measurements where even-indexed values +## (which become the "magnitude" in flux.uncertainty) span a requested range. +make_flux_series <- function(range_lo, range_hi, n_pairs = 80, sd_noise = 0.2, seed = 1) { + set.seed(seed) + base <- seq(range_lo, range_hi, length.out = n_pairs) + # interleave: odd positions = one reading, even positions = paired reading + meas <- as.vector(rbind( + base + stats::rnorm(n_pairs, sd = sd_noise), + base + stats::rnorm(n_pairs, sd = sd_noise) + )) + QC <- rep(0L, length(meas)) + list(measurement = meas, QC = QC) +} + +test_that("flux.uncertainty returns cleanly when all fluxes are negative (mp never assigned)", { + d <- make_flux_series(range_lo = -20, range_hi = -2) + + # Before the fix this would throw "object 'mp' not found" + result <- suppressMessages( + capture.output( + val <- flux.uncertainty(d$measurement, QC = d$QC, bin.num = 5, minBin = 3), + type = "output" + ) + ) + + expect_type(val, "list") + expect_true(all(c("mag", "err", "bias", "n", "intercept") %in% names(val)), + label = "required list elements present") + # No positive bins -> slopeP should not be set + expect_false("slopeP" %in% names(val), + label = "slopeP absent when no positive-flux bins exist") +}) + +test_that("flux.uncertainty fits both slopes for mixed positive/negative flux", { + d <- make_flux_series(range_lo = -15, range_hi = 15) + + capture.output( + val <- flux.uncertainty(d$measurement, QC = d$QC, bin.num = 8, minBin = 3), + type = "output" + ) + + expect_type(val, "list") + expect_true("slopeP" %in% names(val), label = "slopeP present for mixed data") + expect_true("slopeN" %in% names(val), label = "slopeN present for mixed data") + expect_true(is.numeric(val$slopeP) && length(val$slopeP) == 1) + expect_true(is.numeric(val$slopeN) && length(val$slopeN) == 1) +}) + +test_that("flux.uncertainty uses positive-slope fallback for slopeN when neg bins all NA", { + # Data spans a zero crossing but stays almost entirely positive. + # With bin.num=8 and minBin=5 the tiny negative portion falls below the + # minimum bin size, so errBin[neg] is NA for every negative bin. + # The else-if branch therefore assigns slopeN = mp$coefficients[1], + # making slopeN equal to slopeP. + d <- make_flux_series(range_lo = -0.5, range_hi = 20, seed = 3) + + capture.output( + val <- flux.uncertainty(d$measurement, QC = d$QC, bin.num = 8, minBin = 5), + type = "output" + ) + + expect_type(val, "list") + expect_true("slopeP" %in% names(val), label = "slopeP present when positive bins have data") + expect_true("slopeN" %in% names(val), label = "slopeN set via mp fallback when neg bins empty") + expect_equal(val$slopeN, val$slopeP, + label = "slopeN equals slopeP when neg bins unavailable and mp is the fallback") +})