Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions modules/uncertainty/NEWS.md
Original file line number Diff line number Diff line change
@@ -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.
Expand Down
19 changes: 9 additions & 10 deletions modules/uncertainty/R/flux_uncertainty.R
Original file line number Diff line number Diff line change
Expand Up @@ -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]
}

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

updated code is going to leave slopeN undefined if the new if is false. Doing so silently is going to cause issues downstream. Seems like there should be an additional else that sets slopeN to something, as well as a similar set of cases for the earlier slopeP. One option would be a sensible default (e.g., mean value from analyzing many site-years of data), another would be a NA or something like that. If we do set the parameter to NA, it would be good to check what existing PEcAn code calls this function to see what additional checks would be useful there.


Expand Down
78 changes: 78 additions & 0 deletions modules/uncertainty/tests/testthat/test.flux_uncertainty.R
Original file line number Diff line number Diff line change
@@ -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")
})
Loading