From 25993176ca06d7b659c7dcccb3526e379a1124bd Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" Date: Sat, 24 Oct 2026 11:21:45 +0000 Subject: [PATCH 1/7] bench(null_calibration): a package test that cannot run must say why The first native-R run reported every glmgampoi_ql taxon as "failed" with no reason recorded, which is indistinguishable from a fit problem when the CSV is read later - the silent skip the bench's own header forbids. Both package tests now put the error message in a note column, print it, and the cell loop warns when a requested test produced no usable p-value at all. Refs #261 Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- bench/null_calibration/calibration.R | 45 +++++++++++++++++++++++----- 1 file changed, 38 insertions(+), 7 deletions(-) diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R index 60e85d4..a6fc9dc 100644 --- a/bench/null_calibration/calibration.R +++ b/bench/null_calibration/calibration.R @@ -330,6 +330,15 @@ fit_nb_candidates <- function(counts, offset, n_per, want) { # "PValue". A match that missed either spelling would mark every taxon "failed" # and look like a fit problem, so the column is resolved case-insensitively and # a missing column is an error in THIS file, not a data outcome. +# A package test that errors must say why. Returning a frame of "failed" rows +# with no reason is the silent skip this bench forbids: it is indistinguishable +# from a fit problem when the results are read later. +loud_failure <- function(out, test, message) { + cat(" ", test, "failed for every taxon: ", message, "\n", sep = "") + out$note <- message + out +} + p_value_column <- function(tbl, caller) { wanted <- c("pvalue", "p.value", "pval", "p_value") hit <- match(wanted, tolower(names(tbl))) @@ -342,7 +351,8 @@ p_value_column <- function(tbl, caller) { fit_glmgampoi_ql <- function(counts, offset, n_per) { g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) - out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_) + out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_, + note = NA_character_) screens <- vapply(seq_len(ncol(counts)), function(j) nzchar(screen(counts[, j], g)), logical(1)) res <- tryCatch(suppressWarnings({ @@ -354,9 +364,11 @@ fit_glmgampoi_ql <- function(counts, offset, n_per) { # alike must never hand it a label. glmGamPoi::test_de(gp, contrast = c(0, 1), sort_by = "none") }), error = function(e) e) - if (inherits(res, "error")) return(out) + if (inherits(res, "error")) return(loud_failure(out, "glmgampoi_ql", conditionMessage(res))) pcol <- p_value_column(res, "glmGamPoi::test_de") - if (length(res[[pcol]]) != ncol(counts)) return(out) + if (length(res[[pcol]]) != ncol(counts)) { + return(loud_failure(out, "glmgampoi_ql", "test_de returned the wrong number of rows")) + } p <- as.numeric(res[[pcol]]) ok <- !screens & is.finite(p) & p >= 0 & p <= 1 out$status[ok] <- "ok" @@ -366,7 +378,8 @@ fit_glmgampoi_ql <- function(counts, offset, n_per) { fit_edger_ql <- function(counts, offset, n_per) { g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) - out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_) + out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_, + note = NA_character_) screens <- vapply(seq_len(ncol(counts)), function(j) nzchar(screen(counts[, j], g)), logical(1)) res <- tryCatch(suppressWarnings({ @@ -382,9 +395,11 @@ fit_edger_ql <- function(counts, offset, n_per) { qlf <- edgeR::glmQLFTest(fit, coef = 2) edgeR::topTags(qlf, n = Inf, sort.by = "none")$table }), error = function(e) e) - if (inherits(res, "error")) return(out) + if (inherits(res, "error")) return(loud_failure(out, "edger_ql", conditionMessage(res))) pcol <- p_value_column(res, "edgeR::topTags") - if (length(res[[pcol]]) != ncol(counts)) return(out) + if (length(res[[pcol]]) != ncol(counts)) { + return(loud_failure(out, "edger_ql", "topTags returned the wrong number of rows")) + } p <- as.numeric(res[[pcol]]) ok <- !screens & is.finite(p) & p >= 0 & p <= 1 out$status[ok] <- "ok" @@ -487,10 +502,12 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { if ("glmgampoi_ql" %in% enabled) fits$glmgampoi_ql <- fit_glmgampoi_ql(counts, offset, n_per) if ("edger_ql" %in% enabled) fits$edger_ql <- fit_edger_ql(counts, offset, n_per) rows <- do.call(rbind, lapply(names(fits), function(t) { + f <- fits[[t]] + note <- if ("note" %in% names(f)) f$note else rep(NA_character_, nrow(f)) data.frame(truth = truth, shape = shape, offset = off_m, n_per = n_per, taxa = n_taxa, sim = s, taxon = seq_len(ncol(counts)), true_log2fc = sim$true_lfc, test = t, - status = fits[[t]]$status, pvalue = fits[[t]]$pvalue, + status = f$status, pvalue = f$pvalue, note = note, stringsAsFactors = FALSE) })) write.table(rows, f, sep = ",", row.names = FALSE, @@ -500,5 +517,19 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { } cat(sprintf("cell %s done: %d sims skipped (<5 testable taxa), %d offsets refused\n", cell, skipped_sims, refused_offsets)) + # A requested test that produced no usable p-value in the whole cell is a + # broken package path, not a data outcome: say so here and in REPORT.md terms. + cell_rows <- tryCatch(read.csv(f, stringsAsFactors = FALSE), error = function(e) NULL) + if (!is.null(cell_rows) && nrow(cell_rows) > 0) { + for (t in unique(cell_rows$test)) { + usable <- sum(!is.na(cell_rows$pvalue[cell_rows$test == t])) + if (usable == 0) { + reasons <- if ("note" %in% names(cell_rows)) + unique(stats::na.omit(cell_rows$note[cell_rows$test == t])) else character(0) + cat(sprintf("WARNING %s: test %s produced no p-value in this cell%s\n", + cell, t, if (length(reasons)) paste0(": ", paste(reasons, collapse = "; ")) else "")) + } + } + } } cat("calibration run complete\n") From 8784a117f516090eec968848a7bc1eebc30f3d34 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" Date: Sat, 24 Oct 2026 11:50:40 +0000 Subject: [PATCH 2/7] bench(null_calibration): fix glmGamPoi test_de and surface its failure reason Every glmgampoi_ql taxon failed in the native-R run. The cause is the call, not the data: test_de's sort_by names a column to sort by, or NULL for input order, so 'sort_by = "none"' makes the call error every time. It is dropped; rows keep input order, which is what the taxon-by-taxon comparison assumes. The reason was invisible because it only reached the per-cell CSV inside an artifact. The lane now publishes the distinct failure reasons from the raw cells as warnings, so a broken package path is readable where a run is read. Refs #261 Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .../null-calibration-comparators.yml | 19 +++++++++++++++++++ bench/null_calibration/calibration.R | 11 ++++++++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/.github/workflows/null-calibration-comparators.yml b/.github/workflows/null-calibration-comparators.yml index 7cea561..3453951 100644 --- a/.github/workflows/null-calibration-comparators.yml +++ b/.github/workflows/null-calibration-comparators.yml @@ -233,6 +233,25 @@ jobs: bash scripts/ci/run-with-log.sh "$RUNNER_TEMP/summarise.log" -- \ Rscript --no-init-file bench/null_calibration/summarise.R + # A test that produced no usable p-value leaves its reason only in the + # per-cell CSV, which is an artifact a reader may not be able to open. + # The distinct reasons are therefore published as annotations too. + - name: Report any test that could not run + if: always() + run: | + set -uo pipefail + shopt -s nullglob + for csv in bench/null_calibration/results/raw/pvalues_*.csv; do + if grep -q ",failed," "$csv"; then + cell="$(basename "$csv" .csv)" + while IFS= read -r line; do + body="$(printf '%s' "$line" | sed -e 's/%/%25/g' | awk '{printf "%s%%0A", $0}')" + echo "::warning title=${cell}: a requested test failed::${body}" + done < <(awk -F, 'NR > 1 && $10 == "failed" && $12 != "NA" && $12 != "" { + key = $9 " :: " $12; if (!seen[key]++) print key }' "$csv") + fi + done + - name: Publish the summaries as annotations if: always() run: | diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R index a6fc9dc..165c1ca 100644 --- a/bench/null_calibration/calibration.R +++ b/bench/null_calibration/calibration.R @@ -360,9 +360,14 @@ fit_glmgampoi_ql <- function(counts, offset, n_per) { size_factors = exp(offset), overdispersion_shrinkage = TRUE) # Numeric contrast weights on the design columns ((Intercept), gB): - # test_de parses a TEXT contrast as R code, so bench and production - # alike must never hand it a label. - glmGamPoi::test_de(gp, contrast = c(0, 1), sort_by = "none") + # test_de parses a TEXT contrast as R code, so bench and production alike + # must never hand it a label. + # + # No sort_by argument. test_de's sort_by names a COLUMN to sort by, or + # NULL for input order; the string "none" is neither, so passing it makes + # test_de error on every call. Rows then keep input order, which is what + # the taxon-by-taxon comparison below assumes. + glmGamPoi::test_de(gp, contrast = c(0, 1)) }), error = function(e) e) if (inherits(res, "error")) return(loud_failure(out, "glmgampoi_ql", conditionMessage(res))) pcol <- p_value_column(res, "glmGamPoi::test_de") From b0910e7e60594f36d052aca3177995536de95234 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 24 Oct 2026 12:31:50 +0000 Subject: [PATCH 3/7] docs(null_calibration): record native comparator evidence Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- bench/null_calibration/REPORT.md | 85 +++++++++---------- bench/null_calibration/UPSTREAM_DRAFT.md | 69 ++++++++++----- ..._comparator_provenance_run_40091522843.csv | 3 + ...ive_comparator_summary_run_40091522843.csv | 13 +++ 4 files changed, 105 insertions(+), 65 deletions(-) create mode 100644 bench/null_calibration/results/native_comparator_provenance_run_40091522843.csv create mode 100644 bench/null_calibration/results/native_comparator_summary_run_40091522843.csv diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md index 8d926f5..64b43d0 100644 --- a/bench/null_calibration/REPORT.md +++ b/bench/null_calibration/REPORT.md @@ -1,15 +1,18 @@ # Null-calibration bench report (issue #193) -This report records the one run that has actually been executed, superseding the -reduced 2026-10-22 run in this file's history. Anything not listed here has not -been run. The committed summary CSVs in `results/` -(`calibration_summary.csv`, `power_summary.csv`) are the source of truth for -every number below; the per-cell run records are `results/provenance_summary.csv` -(generated from the per-cell provenance files, not hand-written). The raw -p-value files under `results/raw/` are gitignored and reproducible from the -recorded seeds and arguments. - -## Run identity +This report combines the original 2026-10-23 WebR simulation study and a +supplemental native-R/Bioconductor comparator run completed on 2026-10-24. +The reduced 2026-10-22 run remains superseded. The committed +`calibration_summary.csv`, `power_summary.csv`, and `provenance_summary.csv` +cover the original WebR cells. The native comparator +run is preserved separately in +`results/native_comparator_summary_run_40091522843.csv` and +`results/native_comparator_provenance_run_40091522843.csv`, captured from the +successful check-run annotations for [Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843). +Raw p-value files remain outside Git; the run's `null-calibration-raw` artifact +is retained through 2026-11-23. + +## Original WebR baseline run identity - **Date:** 2026-10-23 (UTC, per-cell provenance records). - **Executor:** WebR 0.6.0 (R 4.6.0, wasm32) through `webr_run.mjs`, because the @@ -21,12 +24,12 @@ recorded seeds and arguments. - **Inputs pinned by digest:** `bench/null_calibration/calibration.R` md5 `2717e709ff718e57e602eeaff2410e71`; extracted `_FIT_R` block md5 `cc9efdfdd6e4efe2f9a393ba4d178ce0`. -- **Not available in this sandbox:** `glmGamPoi` and `edgeR` (compiled - Bioconductor packages). `wald_trend`, `glmgampoi_ql` and `edger_ql` were - therefore not run; their code paths in `calibration.R` are unexecuted here. - §"What was not run" gives the commands for a native R/Bioconductor run. -- **Tests run in every cell:** `wald_ml` (the production `_FIT_R` block, - `da_use_gp = FALSE`), `lrt` and `ql_f` (both from one full `glm.nb` fit, plus +- **Not available in the WebR baseline:** `glmGamPoi` and `edgeR` (compiled + Bioconductor packages). The baseline did not run `wald_trend`, + `glmgampoi_ql` or `edger_ql`; the supplemental native-R run below measures + them for original/TSS/null cells at n=4 and n=12. +- **Tests run in every baseline cell:** `wald_ml` (the production `_FIT_R` + block, `da_use_gp = FALSE`), `lrt` and `ql_f` (both from one full `glm.nb` fit, plus one intercept-only fit shared by `lrt`). - **Seeds:** `seed0 = 20261020`; per-sim seed = `(seed0 + cell_key × 1000003 + sim) mod 2³¹−1`, where `cell_key` is a hash of @@ -167,40 +170,34 @@ Readings: |t| > 2.45. The calibration at that size is bought with power, and no measured option is a free lunch at 3–4 samples per group. -## What was not run, and why +## Native R/Bioconductor comparator run (issue #261) + +The comparator lane completed successfully in [Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843), triggered by tag `null-calibration-run4-2026-10-24` at measured commit `bc435feb4614c29f4339fe578caba0c938b6e4c9`. It ran 1,000 null simulations × 40 taxa at original shape/TSS offset for n=4 and n=12 per group, using the same per-cell seed scheme as the WebR study. Its summary and provenance are preserved in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. + +The run used the final two post-PR fixes from that tag lineage: the loud failure-reason reporting commit `17ce140` and `test_de` `sort_by` fix `bc435fe`. They were cherry-picked onto this work branch as `2599317` and `8784a11`; the resulting `calibration.R` md5 matches the measured run (`12a90679f74d2e1a9ab28cec643b5e76`). -`glmGamPoi` and `edgeR` are compiled Bioconductor packages and are not -installable in this sandbox (no native R, and the WebAssembly package repository -does not carry them). Their bench code paths are therefore unexecuted here: +| n/group | test | package | p-values (failed fits) | α=0.05 (Wilson 95% CI) | verdict | +|---:|---|---|---:|---|---| +| 4 | `glmgampoi_ql` | glmGamPoi 1.22.0 | 39,904 (89) | **0.0525** (0.0503–0.0547) | anti-conservative, narrowly | +| 4 | `edger_ql` | edgeR 4.8.2 | 39,904 (89) | **0.0687** (0.0663–0.0713) | anti-conservative | +| 4 | `wald_trend` | glmGamPoi 1.22.0 | 39,904 (89) | **0.0838** (0.0811–0.0865) | anti-conservative | +| 12 | `glmgampoi_ql` | glmGamPoi 1.22.0 | 40,000 (0) | 0.0499 (0.0478–0.0520) | not flagged | +| 12 | `edger_ql` | edgeR 4.8.2 | 40,000 (0) | **0.0712** (0.0687–0.0737) | anti-conservative | +| 12 | `wald_trend` | glmGamPoi 1.22.0 | 40,000 (0) | **0.0764** (0.0739–0.0791) | anti-conservative | -```sh -Rscript -e 'renv::restore()' # native R with the lockfile -Rscript bench/null_calibration/calibration.R \ - --shapes=original --offsets=tss --n-per=4,12 --sims=1000 --taxa=40 \ - --tests=glmgampoi_ql,edger_ql,wald_trend -Rscript bench/null_calibration/summarise.R -``` +The verdict rule remains the one declared above: a Wilson 95% lower bound above 0.05 is anti-conservative; an interval containing 0.05 is “not flagged”, not proof of exact calibration. The glmGamPoi QL result is just above nominal at n=4 (its lower bound is 0.0503) and not flagged at n=12. edgeR QL is anti-conservative at both sizes. `wald_trend` is the glmGamPoi dispersion-trend Wald comparator, **not** glmGamPoi's QL test. -They skip *loudly* when absent (a stdout line plus a provenance record), so the -committed results never silently imply they ran. Issue #193's own measurement -remains the reference for glmGamPoi's QL F-test (0.046 and 0.048 at 4 and 12 -per group); this bench adds a dispersion-scaled QL F test that can be run -without Bioconductor. Both land within about a point of nominal, against 0.14 -and 0.075 for the production Wald test at the same sizes. +Provenance: native R 4.5.0 (2025-04-11), x86_64 Linux, MASS 7.3.65, glmGamPoi 1.22.0, edgeR 4.8.2; `seed0=20261020`; `calibration.R` md5 `12a90679f74d2e1a9ab28cec643b5e76`; extracted production `_FIT_R` md5 `cc9efdfdd6e4efe2f9a393ba4d178ce0`. The run tested all six methods in its two cells; the full 12-row summary is retained in the run-specific CSV. The raw p-value artifact remains on Actions with 30-day retention, through 2026-11-23. + +## What remains unmeasured + +The successful native run fills issue #261's requested glmGamPoi QL and edgeR QL null-calibration cells (and also records `wald_trend`) for n=4 and n=12 only. Those package comparators have not been measured here at other group sizes, shapes, offsets, or under the power/alternative design. The original 16-cell WebR study still uses 400 simulations for the shape/offset/power cells; that limitation is explicit above. The native null-only result does not establish power, FDR under alternatives, broad robustness, or a production-default recommendation. ## Checklist mapping (#193) -Done: the 1,000-sim direct replication at six group sizes with Wilson intervals -at α = 0.01/0.05/0.10 and KS tests; constant, sparse and unequal-library shapes; -RLE as well as TSS offsets; the LRT and an additional quasi-likelihood F test -paired on identical datasets; power and BH-FDR at \|log2FC\| = 1; the bench lives -in `bench/null_calibration/` and never runs in CI. - -Not done here: glmGamPoi QL F, edgeR QL F and the production `glmGamPoi` -dispersion-trend Wald option, all of which need native R + Bioconductor -(commands above); the full 1,000-sim count on the shape/offset/power cells -(400 each — sufficient to resolve a 4-point difference at α = 0.05, and stated -above rather than implied). +Done: the 1,000-sim direct replication at six group sizes with Wilson intervals at α = 0.01/0.05/0.10 and KS tests; constant, sparse and unequal-library shapes; RLE as well as TSS offsets; the LRT and tagwise quasi-likelihood F test paired on identical datasets; power and BH-FDR at |log2FC| = 1; and the native-R glmGamPoi QL, edgeR QL, and `wald_trend` null cells at n=4 and n=12. The bench remains opt-in and never runs in ordinary CI. + +Not done here: power/alternative calibration or other shapes/offsets for the native Bioconductor comparators; the full 1,000-simulation count on the original shape/offset/power cells (400 each — sufficient to resolve a 4-point difference at α = 0.05, as stated rather than implied); and any independent scientific review or upstream default decision. ## Consequence for the production caveat diff --git a/bench/null_calibration/UPSTREAM_DRAFT.md b/bench/null_calibration/UPSTREAM_DRAFT.md index 51ac320..21a7801 100644 --- a/bench/null_calibration/UPSTREAM_DRAFT.md +++ b/bench/null_calibration/UPSTREAM_DRAFT.md @@ -1,9 +1,10 @@ # Draft upstream report: null calibration of the `glm.nb` Wald test -*Status: draft prepared downstream for review before anything is filed -upstream. Every number is from `REPORT.md` in this directory, which defines -exactly what ran, where, and what did not run. This replaces the earlier draft -based on the reduced 2026-10-22 run.* +*Status: draft prepared downstream for maintainer review; it has not been filed +upstream. It includes the completed native comparator run #40091522843. +`REPORT.md` and the run-specific summary/provenance CSVs define exactly what ran, +where, and what remains unmeasured. This replaces the earlier draft based on +the reduced 2026-10-22 run and does not recommend an automatic default change.* ## Summary @@ -28,6 +29,21 @@ Two cheaper alternatives, measured on identical datasets: 0.055 at 4, 0.058–0.060 through 8, 0.056 at 20. It remains ~1 point above nominal, not at nominal. +## Native Bioconductor comparators + +The successful [native-R Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843) measured the same seeded null designs at original shape/TSS, n=4 and n=12, with 1,000 simulations × 40 taxa per cell (R 4.5.0, glmGamPoi 1.22.0, edgeR 4.8.2): + +| n/group | comparator | α=0.05 (Wilson 95% CI) | verdict | +|---:|---|---|---| +| 4 | glmGamPoi QL F (`glmgampoi_ql`) | **0.0525** (0.0503–0.0547) | narrowly anti-conservative | +| 4 | edgeR QL F (`edger_ql`) | **0.0687** (0.0663–0.0713) | anti-conservative | +| 4 | glmGamPoi trend Wald (`wald_trend`) | **0.0838** (0.0811–0.0865) | anti-conservative | +| 12 | glmGamPoi QL F (`glmgampoi_ql`) | 0.0499 (0.0478–0.0520) | not flagged | +| 12 | edgeR QL F (`edger_ql`) | **0.0712** (0.0687–0.0737) | anti-conservative | +| 12 | glmGamPoi trend Wald (`wald_trend`) | **0.0764** (0.0739–0.0791) | anti-conservative | + +“Not flagged” means the Wilson interval includes 0.05; it is not proof of exact calibration. Rejection rates are conditional on successful fits (89 failed fits among 39,993 attempted at n=4; none among 40,000 at n=12). The native comparator lane tested nulls only: it did not measure power, BH-FDR under alternatives, other sample sizes, or other shapes/offsets. The complete 12-row output and per-cell provenance are retained in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. + The inflation is not specific to one data shape: constant dispersion, sparse zero-heavy tables, unequal library sizes and RLE offsets all reproduce it (4/group Wald 0.112–0.142; 12/group 0.070–0.081). RLE is not a remedy; at @@ -47,9 +63,11 @@ Under a power design (25% of taxa at |log2FC| = 1, 400 sims), BH at 0.05: | 12 | QL F (`ql_f`) | 0.433 | 0.063 | The default's apparent power advantage at small n is an artefact of an -uncontrolled FDR. The QL F test is the only measured option whose FDR is near -nominal, but at 4/group it is nearly powerless after BH because its F reference -has only six residual degrees of freedom (`F(1,6)` ≈ 5.99, |t| > 2.45). +uncontrolled FDR. Among the three methods included in this power design, the +tagwise QL F test is closest to nominal FDR, but at 4/group it is nearly +powerless after BH because its F reference has only six residual degrees of +freedom (`F(1,6)` ≈ 5.99, |t| > 2.45). The native Bioconductor comparator run +was null-only and did not measure their power or FDR under alternatives. ## Options for the upstream maintainer's consideration @@ -61,17 +79,22 @@ has only six residual degrees of freedom (`F(1,6)` ≈ 5.99, |t| > 2.45). but still anti-conservative everywhere, slower, and it fails more often at tiny n. It is a mitigation, not a calibration fix; the numbers above quantify the trade-off. -3. **Quasi-likelihood F test below a group-size threshold.** Closest to nominal - at every size and the only option with near-nominal BH FDR, but it changes - critical values substantially at 3–4 samples per group and costs almost all - BH power there. glmGamPoi's own QL F-test (the calibrated path in the - downstream issue) is the version worth considering rather than a tagwise - reimplementation, and it was not measured in this sandbox (see `REPORT.md` - §"What was not run" for the native-R command). - -The downstream recommendation is **1 now, 2 and 3 as an upstream decision**. -A default change would need the glmGamPoi QL F row measured on native R first; -this bench provides the code path and the protocol for it. +3. **Consider glmGamPoi's QL F test below a group-size threshold.** The native + null run places it closest to nominal among the two Bioconductor QL tests: + 0.0525 at 4/group (Wilson 95% CI 0.0503–0.0547, narrowly anti-conservative) + and 0.0499 at 12/group (0.0478–0.0520, not flagged). edgeR QL is higher at + 0.0687 and 0.0712. The separate `wald_trend` row is a glmGamPoi + dispersion-trend Wald test, not its QL test; the production fixed-θ Wald + path must not be described as glmGamPoi QL. These comparator cells are + null-only: they do not measure power/BH behaviour, other sample sizes, or + other shapes/offsets. The tagwise `ql_f` power results therefore cannot be + transferred to glmGamPoi QL. The measured row makes option 3 discussable, + not an automatic default recommendation. + +The downstream recommendation remains **1 now, 2 and 3 as an upstream +maintainer decision**. The native run fills the requested QL calibration cells, +but it does not establish a default threshold or quantify the power trade-off +for glmGamPoi QL. No production default change is proposed by this draft. ## Reproducing @@ -83,6 +106,10 @@ Rscript bench/null_calibration/calibration.R --shapes=original --offsets=tss \ Rscript bench/null_calibration/summarise.R ``` -The committed `results/` CSVs carry every rate and interval quoted above; -`results/provenance_summary.csv` carries the executor, the package versions and -the input digests for each cell. +The committed `results/` CSVs carry the baseline rates and intervals. The +native comparator output and provenance are recorded separately in +`results/native_comparator_summary_run_40091522843.csv` and +`results/native_comparator_provenance_run_40091522843.csv`; the source is [Actions +run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843), +tag `null-calibration-run4-2026-10-24`. The raw Actions artifact is retained +through 2026-11-23. diff --git a/bench/null_calibration/results/native_comparator_provenance_run_40091522843.csv b/bench/null_calibration/results/native_comparator_provenance_run_40091522843.csv new file mode 100644 index 0000000..a25efe8 --- /dev/null +++ b/bench/null_calibration/results/native_comparator_provenance_run_40091522843.csv @@ -0,0 +1,3 @@ +"cell","date_utc","args","seed0","cell_key","executor","calibration_md5","fit_r_md5","mass_mode","glmGamPoi","edgeR","tests_run","r_version","platform" +"null_original_tss_n12_t40","2026-10-24 12:06:12 UTC","--shapes=original --offsets=tss --n-per=4,12 --sims=1000 --taxa=40 --tests=wald_ml,lrt,ql_f,glmgampoi_ql,edger_ql,wald_trend --truth=null","20261020","488430","R version 4.5.0 (2025-04-11)","12a90679f74d2e1a9ab28cec643b5e76","cc9efdfdd6e4efe2f9a393ba4d178ce0","installed MASS 7.3.65","1.22.0","4.8.2","wald_ml,lrt,ql_f,glmgampoi_ql,edger_ql,wald_trend","R version 4.5.0 (2025-04-11)","x86_64-pc-linux-gnu" +"null_original_tss_n4_t40","2026-10-24 11:52:49 UTC","--shapes=original --offsets=tss --n-per=4,12 --sims=1000 --taxa=40 --tests=wald_ml,lrt,ql_f,glmgampoi_ql,edger_ql,wald_trend --truth=null","20261020","940291","R version 4.5.0 (2025-04-11)","12a90679f74d2e1a9ab28cec643b5e76","cc9efdfdd6e4efe2f9a393ba4d178ce0","installed MASS 7.3.65","1.22.0","4.8.2","wald_ml,lrt,ql_f,glmgampoi_ql,edger_ql,wald_trend","R version 4.5.0 (2025-04-11)","x86_64-pc-linux-gnu" diff --git a/bench/null_calibration/results/native_comparator_summary_run_40091522843.csv b/bench/null_calibration/results/native_comparator_summary_run_40091522843.csv new file mode 100644 index 0000000..aea8422 --- /dev/null +++ b/bench/null_calibration/results/native_comparator_summary_run_40091522843.csv @@ -0,0 +1,13 @@ +"truth","shape","offset","n_per_group","n_taxa_setting","test","n_sims","n_attempted","n_pvalues","n_failed","n_boundary","fail_rate","rate_01","rate_01_lo","rate_01_hi","rate_05","rate_05_lo","rate_05_hi","verdict_05","rate_10","rate_10_lo","rate_10_hi","ks_D","ks_p","mean_p_null" +"null","original","tss",4,40,"edger_ql",1000,39993,39904,89,0,0.00222538944315255,0.0166900561347233,0.0154788401059137,0.01799431729848,0.0687399759422614,0.066298813340031,0.0712641632114681,"anti-conservative",0.1264033680834,0.123178859989955,0.129699799689849,0.0358925874707971,4.45754819037285e-45,0.476462218351502 +"null","original","tss",4,40,"glmgampoi_ql",1000,39993,39904,89,0,0.00222538944315255,0.0109512830793905,0.00997619219156174,0.0120205239106234,0.052450882117081,0.0503062992178485,0.054681625602491,"anti-conservative",0.101318163592622,0.0983957875809282,0.10431729243289,0.00539507158196134,0.195782827354291,0.501771815218351 +"null","original","tss",4,40,"lrt",1000,39993,38478,1515,0,0.0378816292851249,0.0320962627995218,0.030381335324203,0.0339046074747647,0.101148708352825,0.0981756413272795,0.104201406225251,"anti-conservative",0.169889287384999,0.166170023962165,0.173674457565967,0.0888304287577666,3.7725013555616e-264,0.441804862475957 +"null","original","tss",4,40,"ql_f",1000,39993,38897,1096,0,0.0274047958392719,0.00925521248425328,0.00835086608108642,0.0102564810012879,0.0545543358099596,0.0523410510614498,0.0568555960937881,"anti-conservative",0.111705272900224,0.108613096510505,0.114874137513051,0.0265496992477194,3.06254627624567e-24,0.483428398794458 +"null","original","tss",4,40,"wald_ml",1000,39993,38897,1096,536,0.0274047958392719,0.0646836516955035,0.0622820140379626,0.0671712643480973,0.139650872817955,0.136241760867281,0.14313115373471,"anti-conservative",0.203280458647196,0.199310489016049,0.207309030398846,0.110379139142975,0,0.430179469668388 +"null","original","tss",4,40,"wald_trend",1000,39993,39904,89,0,0.00222538944315255,0.0280924218123496,0.0265160499044766,0.0297596437117289,0.0837760625501203,0.0810976425888806,0.0865346124836889,"anti-conservative",0.136527666399358,0.133193840316964,0.139931466900361,0.0370445799173282,5.45559297940209e-48,0.487364640551281 +"null","original","tss",12,40,"edger_ql",1000,40000,40000,0,0,0,0.017875,0.0166220864984174,0.0192205077758517,0.071175,0.068696263702048,0.0737360940675417,"anti-conservative",0.13105,0.127778399720934,0.13439245878566,0.040734883798201,4.46657515315624e-58,0.473569758835047 +"null","original","tss",12,40,"glmgampoi_ql",1000,40000,40000,0,0,0,0.0102,0.00926128596865429,0.0112327823238777,0.04985,0.0477601075870381,0.0520263457446958,"not flagged",0.10105,0.0981345874981169,0.104042032643381,0.004913321488283,0.289074033658789,0.498653883472153 +"null","original","tss",12,40,"lrt",1000,40000,39985,15,0,0.000375,0.0143303738902088,0.0132112365828849,0.0155428212210844,0.0629736151056646,0.0606343669360676,0.065396827641314,"anti-conservative",0.121245467050144,0.1180824163193,0.124481286577902,0.0299948859648441,1.13294579102053e-31,0.480516955373746 +"null","original","tss",12,40,"ql_f",1000,40000,39987,13,0,0.000325,0.0123790231825343,0.0113411587666502,0.0135105678429054,0.0576687423412609,0.0554260807775291,0.0599963832277284,"anti-conservative",0.11311176132243,0.11004445661532,0.116253393809745,0.017810937776997,1.91836030034179e-11,0.488360413567825 +"null","original","tss",12,40,"wald_ml",1000,40000,39987,13,0,0.000325,0.0219071198139395,0.0205176447996964,0.0233884445640075,0.074949358541526,0.0724091830369924,0.0775711934698831,"anti-conservative",0.132192962712882,0.128908522375355,0.135548065008805,0.0373066051831862,9.14657842092707e-49,0.47654269902081 +"null","original","tss",12,40,"wald_trend",1000,40000,40000,0,0,0,0.0265,0.0249708723288788,0.0281200654753554,0.076425,0.0738618961794178,0.0790694533040637,"anti-conservative",0.12555,0.122338827066967,0.128833087739359,0.0274444762975241,1.35583467199326e-26,0.499242926997882 From d6b89e39335e0d767f597844f2e0300267650a66 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 24 Oct 2026 12:35:59 +0000 Subject: [PATCH 4/7] refactor(null_calibration): share comparator result handling Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- bench/null_calibration/REPORT.md | 4 +- bench/null_calibration/calibration.R | 108 +++++++++++++++------------ 2 files changed, 62 insertions(+), 50 deletions(-) diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md index 64b43d0..6a06f61 100644 --- a/bench/null_calibration/REPORT.md +++ b/bench/null_calibration/REPORT.md @@ -174,7 +174,9 @@ Readings: The comparator lane completed successfully in [Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843), triggered by tag `null-calibration-run4-2026-10-24` at measured commit `bc435feb4614c29f4339fe578caba0c938b6e4c9`. It ran 1,000 null simulations × 40 taxa at original shape/TSS offset for n=4 and n=12 per group, using the same per-cell seed scheme as the WebR study. Its summary and provenance are preserved in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. -The run used the final two post-PR fixes from that tag lineage: the loud failure-reason reporting commit `17ce140` and `test_de` `sort_by` fix `bc435fe`. They were cherry-picked onto this work branch as `2599317` and `8784a11`; the resulting `calibration.R` md5 matches the measured run (`12a90679f74d2e1a9ab28cec643b5e76`). +The run used the final two post-PR fixes from that tag lineage: the loud failure-reason reporting commit `17ce140` and `test_de` `sort_by` fix `bc435fe`. They were cherry-picked onto this work branch as `2599317` and `8784a11`. The measured file digests are recorded from the run itself: `calibration.R` md5 `12a90679f74d2e1a9ab28cec643b5e76`; extracted production `_FIT_R` md5 `cc9efdfdd6e4efe2f9a393ba4d178ce0`. + +After measurement, the package-result normalization and cell-warning logic were refactored in this PR in response to Sonar findings (shared handling instead of duplicated blocks and nested control flow). The estimator calls, simulation seeds, and formulas are unchanged; the report remains explicitly tied to the measured tag/code digest above. | n/group | test | package | p-values (failed fits) | α=0.05 (Wilson 95% CI) | verdict | |---:|---|---|---:|---|---| diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R index 165c1ca..91344ae 100644 --- a/bench/null_calibration/calibration.R +++ b/bench/null_calibration/calibration.R @@ -325,14 +325,15 @@ fit_nb_candidates <- function(counts, offset, n_per, want) { out } -# Both package tests come back as a data.frame whose p-value column is named by -# the package: glmGamPoi's test_de returns "pval", edgeR's topTags table returns -# "PValue". A match that missed either spelling would mark every taxon "failed" -# and look like a fit problem, so the column is resolved case-insensitively and -# a missing column is an error in THIS file, not a data outcome. -# A package test that errors must say why. Returning a frame of "failed" rows -# with no reason is the silent skip this bench forbids: it is indistinguishable -# from a fit problem when the results are read later. +# Stable test identifiers are constants because they are used by the +# availability, failure-reporting, and result-assembly paths. +TEST_GLMGAMPOI_QL <- "glmgampoi_ql" +TEST_EDGER_QL <- "edger_ql" + +# Both package tests return a data.frame whose p-value column is named by the +# package: glmGamPoi's test_de returns "pval", edgeR's topTags table returns +# "PValue". Resolve it case-insensitively; a missing column is an error in this +# file, not a data outcome. loud_failure <- function(out, test, message) { cat(" ", test, "failed for every taxon: ", message, "\n", sep = "") out$note <- message @@ -349,10 +350,25 @@ p_value_column <- function(tbl, caller) { names(tbl)[hit[1]] } +package_fit_result <- function(res, test, caller, screens, n_taxa, row_error) { + out <- data.frame(status = rep("failed", n_taxa), pvalue = NA_real_, + note = NA_character_) + if (inherits(res, "error")) { + return(loud_failure(out, test, conditionMessage(res))) + } + pcol <- p_value_column(res, caller) + if (length(res[[pcol]]) != n_taxa) { + return(loud_failure(out, test, row_error)) + } + p <- as.numeric(res[[pcol]]) + ok <- !screens & is.finite(p) & p >= 0 & p <= 1 + out$status[ok] <- "ok" + out$pvalue[ok] <- p[ok] + out +} + fit_glmgampoi_ql <- function(counts, offset, n_per) { g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) - out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_, - note = NA_character_) screens <- vapply(seq_len(ncol(counts)), function(j) nzchar(screen(counts[, j], g)), logical(1)) res <- tryCatch(suppressWarnings({ @@ -369,22 +385,12 @@ fit_glmgampoi_ql <- function(counts, offset, n_per) { # the taxon-by-taxon comparison below assumes. glmGamPoi::test_de(gp, contrast = c(0, 1)) }), error = function(e) e) - if (inherits(res, "error")) return(loud_failure(out, "glmgampoi_ql", conditionMessage(res))) - pcol <- p_value_column(res, "glmGamPoi::test_de") - if (length(res[[pcol]]) != ncol(counts)) { - return(loud_failure(out, "glmgampoi_ql", "test_de returned the wrong number of rows")) - } - p <- as.numeric(res[[pcol]]) - ok <- !screens & is.finite(p) & p >= 0 & p <= 1 - out$status[ok] <- "ok" - out$pvalue[ok] <- p[ok] - out + package_fit_result(res, TEST_GLMGAMPOI_QL, "glmGamPoi::test_de", screens, + ncol(counts), "test_de returned the wrong number of rows") } fit_edger_ql <- function(counts, offset, n_per) { g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B")) - out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_, - note = NA_character_) screens <- vapply(seq_len(ncol(counts)), function(j) nzchar(screen(counts[, j], g)), logical(1)) res <- tryCatch(suppressWarnings({ @@ -400,27 +406,38 @@ fit_edger_ql <- function(counts, offset, n_per) { qlf <- edgeR::glmQLFTest(fit, coef = 2) edgeR::topTags(qlf, n = Inf, sort.by = "none")$table }), error = function(e) e) - if (inherits(res, "error")) return(loud_failure(out, "edger_ql", conditionMessage(res))) - pcol <- p_value_column(res, "edgeR::topTags") - if (length(res[[pcol]]) != ncol(counts)) { - return(loud_failure(out, "edger_ql", "topTags returned the wrong number of rows")) + package_fit_result(res, TEST_EDGER_QL, "edgeR::topTags", screens, + ncol(counts), "topTags returned the wrong number of rows") +} + +report_cell_failures <- function(results_file, cell) { + cell_rows <- tryCatch(read.csv(results_file, stringsAsFactors = FALSE), + error = function(e) NULL) + if (is.null(cell_rows) || nrow(cell_rows) == 0) return(invisible(NULL)) + for (test_name in unique(cell_rows$test)) { + test_rows <- cell_rows[cell_rows$test == test_name, , drop = FALSE] + if (any(!is.na(test_rows$pvalue))) next + reasons <- if ("note" %in% names(test_rows)) { + unique(stats::na.omit(test_rows$note)) + } else { + character(0) + } + reason <- if (length(reasons)) paste0(": ", paste(reasons, collapse = "; ")) else "" + cat(sprintf("WARNING %s: test %s produced no p-value in this cell%s\n", + cell, test_name, reason)) } - p <- as.numeric(res[[pcol]]) - ok <- !screens & is.finite(p) & p >= 0 & p <= 1 - out$status[ok] <- "ok" - out$pvalue[ok] <- p[ok] - out + invisible(NULL) } ## ------------------------------------------------------------------- main enabled <- want for (t in want) { - if (t %in% c("wald_trend", "glmgampoi_ql") && !have_gp) { + if (t %in% c("wald_trend", TEST_GLMGAMPOI_QL) && !have_gp) { cat("SKIPPED test", t, ": glmGamPoi is not available in this R library\n") enabled <- setdiff(enabled, t) } - if (t == "edger_ql" && !have_edger) { + if (t == TEST_EDGER_QL && !have_edger) { cat("SKIPPED test edger_ql: edgeR is not available in this R library\n") enabled <- setdiff(enabled, t) } @@ -504,8 +521,12 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { nb <- fit_nb_candidates(counts, offset, n_per, nb_want) for (t in names(nb)) fits[[t]] <- nb[[t]] } - if ("glmgampoi_ql" %in% enabled) fits$glmgampoi_ql <- fit_glmgampoi_ql(counts, offset, n_per) - if ("edger_ql" %in% enabled) fits$edger_ql <- fit_edger_ql(counts, offset, n_per) + if (TEST_GLMGAMPOI_QL %in% enabled) { + fits[[TEST_GLMGAMPOI_QL]] <- fit_glmgampoi_ql(counts, offset, n_per) + } + if (TEST_EDGER_QL %in% enabled) { + fits[[TEST_EDGER_QL]] <- fit_edger_ql(counts, offset, n_per) + } rows <- do.call(rbind, lapply(names(fits), function(t) { f <- fits[[t]] note <- if ("note" %in% names(f)) f$note else rep(NA_character_, nrow(f)) @@ -522,19 +543,8 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { } cat(sprintf("cell %s done: %d sims skipped (<5 testable taxa), %d offsets refused\n", cell, skipped_sims, refused_offsets)) - # A requested test that produced no usable p-value in the whole cell is a - # broken package path, not a data outcome: say so here and in REPORT.md terms. - cell_rows <- tryCatch(read.csv(f, stringsAsFactors = FALSE), error = function(e) NULL) - if (!is.null(cell_rows) && nrow(cell_rows) > 0) { - for (t in unique(cell_rows$test)) { - usable <- sum(!is.na(cell_rows$pvalue[cell_rows$test == t])) - if (usable == 0) { - reasons <- if ("note" %in% names(cell_rows)) - unique(stats::na.omit(cell_rows$note[cell_rows$test == t])) else character(0) - cat(sprintf("WARNING %s: test %s produced no p-value in this cell%s\n", - cell, t, if (length(reasons)) paste0(": ", paste(reasons, collapse = "; ")) else "")) - } - } - } + # A requested test with no usable p-value is a broken package path, not a + # data outcome. Report the reason without adding control-flow nesting here. + report_cell_failures(f, cell) } cat("calibration run complete\n") From 8161ab0eca7c9ce305e9b0775e188b69eb8705a8 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 24 Oct 2026 12:39:13 +0000 Subject: [PATCH 5/7] docs(null_calibration): qualify n=4 denominator Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- bench/null_calibration/REPORT.md | 2 ++ bench/null_calibration/UPSTREAM_DRAFT.md | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md index 6a06f61..477823b 100644 --- a/bench/null_calibration/REPORT.md +++ b/bench/null_calibration/REPORT.md @@ -187,6 +187,8 @@ After measurement, the package-result normalization and cell-warning logic were | 12 | `edger_ql` | edgeR 4.8.2 | 40,000 (0) | **0.0712** (0.0687–0.0737) | anti-conservative | | 12 | `wald_trend` | glmGamPoi 1.22.0 | 40,000 (0) | **0.0764** (0.0739–0.0791) | anti-conservative | +**n=4 denominator note:** the run summary records 39,904 p-values and 89 failed fits (39,993 rows total), leaving seven of the 40,000 planned taxon-simulation positions unaccounted for by those counts. Rejection rates use the 39,904 reported p-values. Do not treat 39,993 as the number attempted or classify the seven missing positions as successes or failures. At n=12, all 40,000 p-values are reported with no failed fits. + The verdict rule remains the one declared above: a Wilson 95% lower bound above 0.05 is anti-conservative; an interval containing 0.05 is “not flagged”, not proof of exact calibration. The glmGamPoi QL result is just above nominal at n=4 (its lower bound is 0.0503) and not flagged at n=12. edgeR QL is anti-conservative at both sizes. `wald_trend` is the glmGamPoi dispersion-trend Wald comparator, **not** glmGamPoi's QL test. Provenance: native R 4.5.0 (2025-04-11), x86_64 Linux, MASS 7.3.65, glmGamPoi 1.22.0, edgeR 4.8.2; `seed0=20261020`; `calibration.R` md5 `12a90679f74d2e1a9ab28cec643b5e76`; extracted production `_FIT_R` md5 `cc9efdfdd6e4efe2f9a393ba4d178ce0`. The run tested all six methods in its two cells; the full 12-row summary is retained in the run-specific CSV. The raw p-value artifact remains on Actions with 30-day retention, through 2026-11-23. diff --git a/bench/null_calibration/UPSTREAM_DRAFT.md b/bench/null_calibration/UPSTREAM_DRAFT.md index 21a7801..15881c4 100644 --- a/bench/null_calibration/UPSTREAM_DRAFT.md +++ b/bench/null_calibration/UPSTREAM_DRAFT.md @@ -42,7 +42,7 @@ The successful [native-R Actions run 40091522843](https://github.com/hyperpolyma | 12 | edgeR QL F (`edger_ql`) | **0.0712** (0.0687–0.0737) | anti-conservative | | 12 | glmGamPoi trend Wald (`wald_trend`) | **0.0764** (0.0739–0.0791) | anti-conservative | -“Not flagged” means the Wilson interval includes 0.05; it is not proof of exact calibration. Rejection rates are conditional on successful fits (89 failed fits among 39,993 attempted at n=4; none among 40,000 at n=12). The native comparator lane tested nulls only: it did not measure power, BH-FDR under alternatives, other sample sizes, or other shapes/offsets. The complete 12-row output and per-cell provenance are retained in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. +“Not flagged” means the Wilson interval includes 0.05; it is not proof of exact calibration. At n=4, the summary reports 39,904 p-values and 89 failed fits (39,993 rows total); seven of the 40,000 planned taxon-simulation positions are not accounted for by those counts. The reported rejection rate uses the 39,904 p-values. Do not treat 39,993 as the number attempted or classify the seven unaccounted positions as successes or failures. At n=12, the summary reports 40,000 p-values and no failed fits. The native comparator lane tested nulls only: it did not measure power, BH-FDR under alternatives, other sample sizes, or other shapes/offsets. The complete 12-row output and per-cell provenance are retained in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. The inflation is not specific to one data shape: constant dispersion, sparse zero-heavy tables, unequal library sizes and RLE offsets all reproduce it From a2e994760a232c7dd054b9e7c7f0c68e47925748 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 24 Oct 2026 12:49:49 +0000 Subject: [PATCH 6/7] fix(null_calibration): harden CSV output and resume validation Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .../null-calibration-comparators.yml | 39 +++++++--- bench/null_calibration/REPORT.md | 22 +++--- bench/null_calibration/UPSTREAM_DRAFT.md | 24 +++--- bench/null_calibration/calibration.R | 74 ++++++++++++++----- 4 files changed, 107 insertions(+), 52 deletions(-) diff --git a/.github/workflows/null-calibration-comparators.yml b/.github/workflows/null-calibration-comparators.yml index 3453951..fcd35fe 100644 --- a/.github/workflows/null-calibration-comparators.yml +++ b/.github/workflows/null-calibration-comparators.yml @@ -239,17 +239,38 @@ jobs: - name: Report any test that could not run if: always() run: | - set -uo pipefail + set -euo pipefail shopt -s nullglob for csv in bench/null_calibration/results/raw/pvalues_*.csv; do - if grep -q ",failed," "$csv"; then - cell="$(basename "$csv" .csv)" - while IFS= read -r line; do - body="$(printf '%s' "$line" | sed -e 's/%/%25/g' | awk '{printf "%s%%0A", $0}')" - echo "::warning title=${cell}: a requested test failed::${body}" - done < <(awk -F, 'NR > 1 && $10 == "failed" && $12 != "NA" && $12 != "" { - key = $9 " :: " $12; if (!seen[key]++) print key }' "$csv") - fi + cell="$(basename "$csv" .csv)" + python3 - "$csv" "$cell" <<'PY' + import csv + import sys + + path, cell = sys.argv[1:] + expected = ["truth", "shape", "offset", "n_per", "taxa", "sim", + "taxon", "true_log2fc", "test", "status", "pvalue", "note"] + seen = set() + with open(path, encoding="utf-8", newline="") as stream: + reader = csv.DictReader(stream) + if reader.fieldnames not in (expected, expected[:-1]): + raise SystemExit(f"{path}: unexpected CSV header {reader.fieldnames!r}") + for record_number, row in enumerate(reader, start=2): + if None in row or any(value is None for value in row.values()): + raise SystemExit(f"{path}: malformed CSV record {record_number}") + if row.get("status") != "failed": + continue + note = row.get("note") + if not note or note == "NA": + continue + key = (row.get("test", ""), note) + if key in seen: + continue + seen.add(key) + message = f"{key[0]} :: {note}" + message = message.replace("%", "%25").replace("\r", "%0D").replace("\n", "%0A") + print(f"::warning title={cell}: a requested test failed::{message}") + PY done - name: Publish the summaries as annotations diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md index 477823b..224172a 100644 --- a/bench/null_calibration/REPORT.md +++ b/bench/null_calibration/REPORT.md @@ -172,22 +172,22 @@ Readings: ## Native R/Bioconductor comparator run (issue #261) -The comparator lane completed successfully in [Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843), triggered by tag `null-calibration-run4-2026-10-24` at measured commit `bc435feb4614c29f4339fe578caba0c938b6e4c9`. It ran 1,000 null simulations × 40 taxa at original shape/TSS offset for n=4 and n=12 per group, using the same per-cell seed scheme as the WebR study. Its summary and provenance are preserved in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. +The comparator lane completed successfully in [Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843), triggered by tag `null-calibration-run4-2026-10-24` at measured commit `bc435feb4614c29f4339fe578caba0c938b6e4c9`. It was configured for 1,000 null simulations × 40 taxa at original shape/TSS offset for n=4 and n=12 per group, using the same per-cell seed scheme as the WebR study. Its summary and provenance are preserved in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. The run used the final two post-PR fixes from that tag lineage: the loud failure-reason reporting commit `17ce140` and `test_de` `sort_by` fix `bc435fe`. They were cherry-picked onto this work branch as `2599317` and `8784a11`. The measured file digests are recorded from the run itself: `calibration.R` md5 `12a90679f74d2e1a9ab28cec643b5e76`; extracted production `_FIT_R` md5 `cc9efdfdd6e4efe2f9a393ba4d178ce0`. -After measurement, the package-result normalization and cell-warning logic were refactored in this PR in response to Sonar findings (shared handling instead of duplicated blocks and nested control flow). The estimator calls, simulation seeds, and formulas are unchanged; the report remains explicitly tied to the measured tag/code digest above. +After measurement, this PR refactored package-result normalization and cell-warning logic in response to Sonar findings, and hardened CSV serialization plus resume validation (schema and cell-key checks) so package error notes remain valid data. The comparator workflow now reads failure rows with a CSV-aware parser. These post-run changes do not alter estimator calls, simulation seeds, or formulas. At this PR revision the working branch's `calibration.R` md5 is `d689c158fcf2c606ed392af005a69d8b`, not the measured run's `12a90679f74d2e1a9ab28cec643b5e76`; the results below remain tied to the measured tag/code digest above. -| n/group | test | package | p-values (failed fits) | α=0.05 (Wilson 95% CI) | verdict | -|---:|---|---|---:|---|---| -| 4 | `glmgampoi_ql` | glmGamPoi 1.22.0 | 39,904 (89) | **0.0525** (0.0503–0.0547) | anti-conservative, narrowly | -| 4 | `edger_ql` | edgeR 4.8.2 | 39,904 (89) | **0.0687** (0.0663–0.0713) | anti-conservative | -| 4 | `wald_trend` | glmGamPoi 1.22.0 | 39,904 (89) | **0.0838** (0.0811–0.0865) | anti-conservative | -| 12 | `glmgampoi_ql` | glmGamPoi 1.22.0 | 40,000 (0) | 0.0499 (0.0478–0.0520) | not flagged | -| 12 | `edger_ql` | edgeR 4.8.2 | 40,000 (0) | **0.0712** (0.0687–0.0737) | anti-conservative | -| 12 | `wald_trend` | glmGamPoi 1.22.0 | 40,000 (0) | **0.0764** (0.0739–0.0791) | anti-conservative | +| n/group | test | package | summary rows (`n_attempted`) | p-values (failed fits) | α=0.05 (Wilson 95% CI) | verdict | +|---:|---|---|---:|---:|---|---| +| 4 | `glmgampoi_ql` | glmGamPoi 1.22.0 | 39,993 | 39,904 (89) | **0.0525** (0.0503–0.0547) | anti-conservative, narrowly | +| 4 | `edger_ql` | edgeR 4.8.2 | 39,993 | 39,904 (89) | **0.0687** (0.0663–0.0713) | anti-conservative | +| 4 | `wald_trend` | glmGamPoi 1.22.0 | 39,993 | 39,904 (89) | **0.0838** (0.0811–0.0865) | anti-conservative | +| 12 | `glmgampoi_ql` | glmGamPoi 1.22.0 | 40,000 | 40,000 (0) | 0.0499 (0.0478–0.0520) | not flagged | +| 12 | `edger_ql` | edgeR 4.8.2 | 40,000 | 40,000 (0) | **0.0712** (0.0687–0.0737) | anti-conservative | +| 12 | `wald_trend` | glmGamPoi 1.22.0 | 40,000 | 40,000 (0) | **0.0764** (0.0739–0.0791) | anti-conservative | -**n=4 denominator note:** the run summary records 39,904 p-values and 89 failed fits (39,993 rows total), leaving seven of the 40,000 planned taxon-simulation positions unaccounted for by those counts. Rejection rates use the 39,904 reported p-values. Do not treat 39,993 as the number attempted or classify the seven missing positions as successes or failures. At n=12, all 40,000 p-values are reported with no failed fits. +The `n_attempted` summary field is `nrow(d)`: it counts rows present in the per-test output, not the 40,000 planned taxon-simulation positions. The `taxa=40` setting is the raw design size, not an observed post-filter testing denominator. At n=4, 39,904 p-values plus 89 failed fits account for 39,993 rows; seven of the 40,000 positions in the configured 1,000 × 40 design are absent from the summary, and their disposition is unresolved. Rejection rates use the 39,904 reported p-values; do not assign the seven absent positions to either successes or failures. At n=12, all 40,000 rows are reported as p-values with no failed fits. The verdict rule remains the one declared above: a Wilson 95% lower bound above 0.05 is anti-conservative; an interval containing 0.05 is “not flagged”, not proof of exact calibration. The glmGamPoi QL result is just above nominal at n=4 (its lower bound is 0.0503) and not flagged at n=12. edgeR QL is anti-conservative at both sizes. `wald_trend` is the glmGamPoi dispersion-trend Wald comparator, **not** glmGamPoi's QL test. diff --git a/bench/null_calibration/UPSTREAM_DRAFT.md b/bench/null_calibration/UPSTREAM_DRAFT.md index 15881c4..ad5b851 100644 --- a/bench/null_calibration/UPSTREAM_DRAFT.md +++ b/bench/null_calibration/UPSTREAM_DRAFT.md @@ -31,18 +31,18 @@ Two cheaper alternatives, measured on identical datasets: ## Native Bioconductor comparators -The successful [native-R Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843) measured the same seeded null designs at original shape/TSS, n=4 and n=12, with 1,000 simulations × 40 taxa per cell (R 4.5.0, glmGamPoi 1.22.0, edgeR 4.8.2): - -| n/group | comparator | α=0.05 (Wilson 95% CI) | verdict | -|---:|---|---|---| -| 4 | glmGamPoi QL F (`glmgampoi_ql`) | **0.0525** (0.0503–0.0547) | narrowly anti-conservative | -| 4 | edgeR QL F (`edger_ql`) | **0.0687** (0.0663–0.0713) | anti-conservative | -| 4 | glmGamPoi trend Wald (`wald_trend`) | **0.0838** (0.0811–0.0865) | anti-conservative | -| 12 | glmGamPoi QL F (`glmgampoi_ql`) | 0.0499 (0.0478–0.0520) | not flagged | -| 12 | edgeR QL F (`edger_ql`) | **0.0712** (0.0687–0.0737) | anti-conservative | -| 12 | glmGamPoi trend Wald (`wald_trend`) | **0.0764** (0.0739–0.0791) | anti-conservative | - -“Not flagged” means the Wilson interval includes 0.05; it is not proof of exact calibration. At n=4, the summary reports 39,904 p-values and 89 failed fits (39,993 rows total); seven of the 40,000 planned taxon-simulation positions are not accounted for by those counts. The reported rejection rate uses the 39,904 p-values. Do not treat 39,993 as the number attempted or classify the seven unaccounted positions as successes or failures. At n=12, the summary reports 40,000 p-values and no failed fits. The native comparator lane tested nulls only: it did not measure power, BH-FDR under alternatives, other sample sizes, or other shapes/offsets. The complete 12-row output and per-cell provenance are retained in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. +The successful [native-R Actions run 40091522843](https://github.com/hyperpolymath/MetaManifold-WebUI/actions/runs/40091522843) used the same seeded null designs at original shape/TSS, n=4 and n=12, configured for 1,000 simulations × 40 taxa per cell (R 4.5.0, glmGamPoi 1.22.0, edgeR 4.8.2): + +| n/group | comparator | summary rows (`n_attempted`) | p-values (failed fits) | α=0.05 (Wilson 95% CI) | verdict | +|---:|---|---:|---:|---|---| +| 4 | glmGamPoi QL F (`glmgampoi_ql`) | 39,993 | 39,904 (89) | **0.0525** (0.0503–0.0547) | narrowly anti-conservative | +| 4 | edgeR QL F (`edger_ql`) | 39,993 | 39,904 (89) | **0.0687** (0.0663–0.0713) | anti-conservative | +| 4 | glmGamPoi trend Wald (`wald_trend`) | 39,993 | 39,904 (89) | **0.0838** (0.0811–0.0865) | anti-conservative | +| 12 | glmGamPoi QL F (`glmgampoi_ql`) | 40,000 | 40,000 (0) | 0.0499 (0.0478–0.0520) | not flagged | +| 12 | edgeR QL F (`edger_ql`) | 40,000 | 40,000 (0) | **0.0712** (0.0687–0.0737) | anti-conservative | +| 12 | glmGamPoi trend Wald (`wald_trend`) | 40,000 | 40,000 (0) | **0.0764** (0.0739–0.0791) | anti-conservative | + +“Not flagged” means the Wilson interval includes 0.05; it is not proof of exact calibration. The summary's `n_attempted` field is computed as `nrow(d)`, the rows present in each per-test output; the `taxa=40` setting is the raw design size, not an observed post-filter testing denominator. At n=4, the summary reports 39,904 p-values and 89 failed fits (39,993 output rows total); seven positions from the configured 40,000 taxon-simulation design are absent from the summary, and their disposition is unresolved. Rejection rates use the 39,904 reported p-values; do not assign the seven absent positions to successes or failures or treat the 40,000 design count as an observed testing denominator. At n=12, the summary reports 40,000 p-values and no failed fits. The native comparator lane tested nulls only: it did not measure power, BH-FDR under alternatives, other sample sizes, or other shapes/offsets. The complete 12-row output and per-cell provenance are retained in `results/native_comparator_summary_run_40091522843.csv` and `results/native_comparator_provenance_run_40091522843.csv`. The inflation is not specific to one data shape: constant dispersion, sparse zero-heavy tables, unequal library sizes and RLE offsets all reproduce it diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R index 91344ae..9015256 100644 --- a/bench/null_calibration/calibration.R +++ b/bench/null_calibration/calibration.R @@ -429,6 +429,50 @@ report_cell_failures <- function(results_file, cell) { invisible(NULL) } +RESULT_COLUMNS <- c("truth", "shape", "offset", "n_per", "taxa", "sim", + "taxon", "true_log2fc", "test", "status", "pvalue", "note") + +validate_resume_cell <- function(results_file, provenance_file, cell, cell_key) { + if (!file.exists(results_file)) return(integer(0)) + prev <- tryCatch(read.csv(results_file, stringsAsFactors = FALSE, + check.names = FALSE, fill = FALSE), + error = function(e) e) + if (inherits(prev, "error")) { + stop("cell ", cell, " has an unreadable CSV (", conditionMessage(prev), + "); delete ", results_file, " and ", provenance_file, + " before rerunning", call. = FALSE) + } + if (!identical(names(prev), RESULT_COLUMNS)) { + stop("cell ", cell, " has an incompatible CSV schema; expected ", + paste(RESULT_COLUMNS, collapse = ", "), " but found ", + paste(names(prev), collapse = ", "), "; delete ", results_file, + " and ", provenance_file, " before rerunning", call. = FALSE) + } + if (!file.exists(provenance_file)) { + stop("cell ", cell, " has rows but no provenance file; delete ", + results_file, " before rerunning", call. = FALSE) + } + old_key <- tryCatch( + grep("^cell key:", readLines(provenance_file, warn = FALSE), value = TRUE), + error = function(e) e) + if (inherits(old_key, "error")) { + stop("cell ", cell, " has unreadable provenance (", conditionMessage(old_key), + "); delete ", results_file, " and ", provenance_file, + " before rerunning", call. = FALSE) + } + if (length(old_key) != 1L) { + stop("cell ", cell, " provenance must contain exactly one cell key; delete ", + results_file, " and ", provenance_file, " before rerunning", call. = FALSE) + } + old_key_value <- trimws(sub("^cell key:[ ]*", "", old_key[1])) + if (!identical(old_key_value, as.character(cell_key))) { + stop("cell ", cell, " has a mismatched cell key (found ", old_key_value, + ", expected ", cell_key, "); delete ", results_file, " and ", + provenance_file, " before rerunning", call. = FALSE) + } + unique(as.integer(prev$sim[!is.na(prev$sim)])) +} + ## ------------------------------------------------------------------- main enabled <- want @@ -461,7 +505,7 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { if (truth == "power") sprintf("_lfc%s", lfc) else "") f <- file.path(out_dir, paste0("pvalues_", cell, ".csv")) prov <- file.path(out_dir, paste0("provenance_", cell, ".txt")) - writeLines(c( + provenance_lines <- c( paste("cell:", cell), paste("date:", format(Sys.time(), tz = "UTC", usetz = TRUE)), paste("args:", paste(argv, collapse = " ")), paste("seed0:", seed0), paste("cell key:", stable_key(cell)), @@ -474,27 +518,14 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { paste("edgeR:", if (have_edger) as.character(packageVersion("edgeR")) else "absent"), paste("tests run:", paste(enabled, collapse = ",")), capture.output(sessionInfo()) - ), prov) + ) cat(sprintf("cell %s: %d sims x %d taxa, tests %s\n", cell, n_sims, n_taxa, paste(enabled, collapse = "+"))) header_done <- file.exists(f) - done_sims <- integer(0) - if (header_done) { - # Resuming is only safe when the existing rows came from the same seed - # scheme. The cell key is recorded in the provenance file; a mismatch means - # a stale cell from a different convention, not a resumable one. - prev <- tryCatch(read.csv(f, stringsAsFactors = FALSE), error = function(e) NULL) - if (!is.null(prev) && "sim" %in% names(prev)) done_sims <- unique(as.integer(prev$sim)) - if (file.exists(prov)) { - old_key <- grep("^cell key:", readLines(prov, warn = FALSE), value = TRUE) - if (length(old_key) && !identical(trimws(sub("^cell key:[ ]*", "", old_key[1])), - as.character(stable_key(cell)))) { - stop("cell ", cell, " has rows from a different seed scheme (key ", - sub("^cell key:[ ]*", "", old_key[1]), ", now ", stable_key(cell), - "); delete ", f, " and ", prov, " before rerunning") - } - } - } + # Validate the existing schema and provenance before refreshing its record or + # appending rows; unknown files cannot safely be resumed. + done_sims <- validate_resume_cell(f, prov, cell, stable_key(cell)) + writeLines(provenance_lines, prov) skipped_sims <- 0L refused_offsets <- 0L if (length(done_sims)) cat(sprintf(" resuming: %d sims already present\n", length(done_sims))) @@ -536,8 +567,11 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { status = f$status, pvalue = f$pvalue, note = note, stringsAsFactors = FALSE) })) + # Package error messages can contain commas or newlines; quote all text + # fields so both read.csv and CSV-aware workflow consumers preserve them. write.table(rows, f, sep = ",", row.names = FALSE, - col.names = !header_done, append = header_done, quote = FALSE) + col.names = !header_done, append = header_done, quote = TRUE, + qmethod = "double") header_done <- TRUE if (s %% 25 == 0) { cat(sprintf(" %s sim %d/%d\n", cell, s, n_sims)); flush.console() } } From f9ce2d9b77b4dc3b4987a2a4c57ff65d970b84aa Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 24 Oct 2026 12:58:35 +0000 Subject: [PATCH 7/7] test(null_calibration): cover CSV failure handling Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- .github/workflows/ci.yml | 6 + .../null-calibration-comparators.yml | 29 +---- bench/null_calibration/REPORT.md | 2 +- bench/null_calibration/calibration.R | 38 ++++-- bench/null_calibration/report_failures.py | 77 +++++++++++++ .../tests/test_null_calibration_failures.py | 89 ++++++++++++++ test/unit/test_null_calibration_csv.R | 109 ++++++++++++++++++ 7 files changed, 310 insertions(+), 40 deletions(-) create mode 100644 bench/null_calibration/report_failures.py create mode 100644 scripts/tests/test_null_calibration_failures.py create mode 100644 test/unit/test_null_calibration_csv.R diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6938cdf..92e9901 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -293,6 +293,12 @@ jobs: - name: Check DADA2 CSV assembly run: Rscript --vanilla test/unit/test_dada2_tables.R + - name: Check null-calibration CSV handling + run: Rscript --vanilla test/unit/test_null_calibration_csv.R + + - name: Check null-calibration failure annotations + run: python3 -m unittest discover -s scripts/tests -p 'test_null_calibration_failures.py' -v + # Saved on failure too, since the packages built before a failure are the ones # the next attempt should reuse. Root writes the cache, so it is made readable first. - name: Make the renv cache readable diff --git a/.github/workflows/null-calibration-comparators.yml b/.github/workflows/null-calibration-comparators.yml index fcd35fe..4e57f3e 100644 --- a/.github/workflows/null-calibration-comparators.yml +++ b/.github/workflows/null-calibration-comparators.yml @@ -243,34 +243,7 @@ jobs: shopt -s nullglob for csv in bench/null_calibration/results/raw/pvalues_*.csv; do cell="$(basename "$csv" .csv)" - python3 - "$csv" "$cell" <<'PY' - import csv - import sys - - path, cell = sys.argv[1:] - expected = ["truth", "shape", "offset", "n_per", "taxa", "sim", - "taxon", "true_log2fc", "test", "status", "pvalue", "note"] - seen = set() - with open(path, encoding="utf-8", newline="") as stream: - reader = csv.DictReader(stream) - if reader.fieldnames not in (expected, expected[:-1]): - raise SystemExit(f"{path}: unexpected CSV header {reader.fieldnames!r}") - for record_number, row in enumerate(reader, start=2): - if None in row or any(value is None for value in row.values()): - raise SystemExit(f"{path}: malformed CSV record {record_number}") - if row.get("status") != "failed": - continue - note = row.get("note") - if not note or note == "NA": - continue - key = (row.get("test", ""), note) - if key in seen: - continue - seen.add(key) - message = f"{key[0]} :: {note}" - message = message.replace("%", "%25").replace("\r", "%0D").replace("\n", "%0A") - print(f"::warning title={cell}: a requested test failed::{message}") - PY + python3 bench/null_calibration/report_failures.py "$csv" "$cell" done - name: Publish the summaries as annotations diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md index 224172a..2627b56 100644 --- a/bench/null_calibration/REPORT.md +++ b/bench/null_calibration/REPORT.md @@ -176,7 +176,7 @@ The comparator lane completed successfully in [Actions run 40091522843](https:// The run used the final two post-PR fixes from that tag lineage: the loud failure-reason reporting commit `17ce140` and `test_de` `sort_by` fix `bc435fe`. They were cherry-picked onto this work branch as `2599317` and `8784a11`. The measured file digests are recorded from the run itself: `calibration.R` md5 `12a90679f74d2e1a9ab28cec643b5e76`; extracted production `_FIT_R` md5 `cc9efdfdd6e4efe2f9a393ba4d178ce0`. -After measurement, this PR refactored package-result normalization and cell-warning logic in response to Sonar findings, and hardened CSV serialization plus resume validation (schema and cell-key checks) so package error notes remain valid data. The comparator workflow now reads failure rows with a CSV-aware parser. These post-run changes do not alter estimator calls, simulation seeds, or formulas. At this PR revision the working branch's `calibration.R` md5 is `d689c158fcf2c606ed392af005a69d8b`, not the measured run's `12a90679f74d2e1a9ab28cec643b5e76`; the results below remain tied to the measured tag/code digest above. +After measurement, this PR refactored package-result normalization and cell-warning logic in response to Sonar findings, and hardened CSV serialization plus resume validation (schema and cell-key checks) so package error notes remain valid data. The comparator workflow now reads failure rows with a CSV-aware parser. These post-run changes do not alter estimator calls, simulation seeds, or formulas. At this PR revision the working branch's `calibration.R` md5 is `3ed22ba8b04f695d531a45925abbc13f`, not the measured run's `12a90679f74d2e1a9ab28cec643b5e76`; the results below remain tied to the measured tag/code digest above. | n/group | test | package | summary rows (`n_attempted`) | p-values (failed fits) | α=0.05 (Wilson 95% CI) | verdict | |---:|---|---|---:|---:|---|---| diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R index 9015256..767c1d6 100644 --- a/bench/null_calibration/calibration.R +++ b/bench/null_calibration/calibration.R @@ -410,10 +410,22 @@ fit_edger_ql <- function(counts, offset, n_per) { ncol(counts), "topTags returned the wrong number of rows") } +RESULT_COLUMNS <- c("truth", "shape", "offset", "n_per", "taxa", "sim", + "taxon", "true_log2fc", "test", "status", "pvalue", "note") + report_cell_failures <- function(results_file, cell) { - cell_rows <- tryCatch(read.csv(results_file, stringsAsFactors = FALSE), - error = function(e) NULL) - if (is.null(cell_rows) || nrow(cell_rows) == 0) return(invisible(NULL)) + if (!file.exists(results_file)) return(invisible(NULL)) + cell_rows <- tryCatch(read.csv(results_file, stringsAsFactors = FALSE, + check.names = FALSE, fill = FALSE), + error = function(e) e) + if (inherits(cell_rows, "error")) { + stop("cell ", cell, " has an unreadable result CSV (", + conditionMessage(cell_rows), ")", call. = FALSE) + } + if (!identical(names(cell_rows), RESULT_COLUMNS)) { + stop("cell ", cell, " has an incompatible result CSV schema", call. = FALSE) + } + if (nrow(cell_rows) == 0) return(invisible(NULL)) for (test_name in unique(cell_rows$test)) { test_rows <- cell_rows[cell_rows$test == test_name, , drop = FALSE] if (any(!is.na(test_rows$pvalue))) next @@ -429,9 +441,6 @@ report_cell_failures <- function(results_file, cell) { invisible(NULL) } -RESULT_COLUMNS <- c("truth", "shape", "offset", "n_per", "taxa", "sim", - "taxon", "true_log2fc", "test", "status", "pvalue", "note") - validate_resume_cell <- function(results_file, provenance_file, cell, cell_key) { if (!file.exists(results_file)) return(integer(0)) prev <- tryCatch(read.csv(results_file, stringsAsFactors = FALSE, @@ -473,6 +482,17 @@ validate_resume_cell <- function(results_file, provenance_file, cell, cell_key) unique(as.integer(prev$sim[!is.na(prev$sim)])) } +write_cell_rows <- function(rows, results_file, header_done) { + if (!identical(names(rows), RESULT_COLUMNS)) { + stop("cell result rows do not match the expected CSV schema", call. = FALSE) + } + # Package error messages can contain commas, quotes, or newlines. Use RFC-style + # double-quote escaping so R and the workflow's CSV reader preserve each field. + write.table(rows, results_file, sep = ",", row.names = FALSE, + col.names = !header_done, append = header_done, quote = TRUE, + qmethod = "double") +} + ## ------------------------------------------------------------------- main enabled <- want @@ -567,11 +587,7 @@ for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) { status = f$status, pvalue = f$pvalue, note = note, stringsAsFactors = FALSE) })) - # Package error messages can contain commas or newlines; quote all text - # fields so both read.csv and CSV-aware workflow consumers preserve them. - write.table(rows, f, sep = ",", row.names = FALSE, - col.names = !header_done, append = header_done, quote = TRUE, - qmethod = "double") + write_cell_rows(rows, f, header_done) header_done <- TRUE if (s %% 25 == 0) { cat(sprintf(" %s sim %d/%d\n", cell, s, n_sims)); flush.console() } } diff --git a/bench/null_calibration/report_failures.py b/bench/null_calibration/report_failures.py new file mode 100644 index 0000000..c62d4c3 --- /dev/null +++ b/bench/null_calibration/report_failures.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: AGPL-3.0-only +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +"""Publish distinct native-calibration failure notes as Actions warnings.""" + +from __future__ import annotations + +import csv +from pathlib import Path +import sys + +RESULT_FIELDS = [ + "truth", + "shape", + "offset", + "n_per", + "taxa", + "sim", + "taxon", + "true_log2fc", + "test", + "status", + "pvalue", + "note", +] +LEGACY_FIELDS = RESULT_FIELDS[:-1] + + +def _escape_annotation(text: str) -> str: + """Escape data characters reserved by the GitHub Actions command protocol.""" + return text.replace("%", "%25").replace("\r", "%0D").replace("\n", "%0A") + + +def failure_annotations(csv_path: str | Path, cell: str) -> list[str]: + """Return one warning for each distinct non-empty (test, failure-note) pair.""" + annotations: list[str] = [] + seen: set[tuple[str, str]] = set() + with Path(csv_path).open(encoding="utf-8", newline="") as stream: + reader = csv.DictReader(stream) + if reader.fieldnames not in (RESULT_FIELDS, LEGACY_FIELDS): + raise ValueError(f"unexpected CSV header {reader.fieldnames!r}") + for record_number, row in enumerate(reader, start=2): + if None in row or any(value is None for value in row.values()): + raise ValueError(f"malformed CSV record {record_number}") + if row.get("status") != "failed": + continue + note = row.get("note") + if not note or note == "NA": + continue + key = (row.get("test", ""), note) + if key in seen: + continue + seen.add(key) + message = _escape_annotation(f"{key[0]} :: {note}") + annotations.append( + f"::warning title={cell}: a requested test failed::{message}" + ) + return annotations + + +def main(argv: list[str] | None = None) -> int: + args = sys.argv[1:] if argv is None else argv + if len(args) != 2: + print(f"usage: {Path(sys.argv[0]).name} RESULTS.csv CELL", file=sys.stderr) + return 2 + csv_path, cell = args + try: + for annotation in failure_annotations(csv_path, cell): + print(annotation) + except (OSError, UnicodeError, csv.Error, ValueError) as error: + print(f"{csv_path}: {error}", file=sys.stderr) + return 2 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/tests/test_null_calibration_failures.py b/scripts/tests/test_null_calibration_failures.py new file mode 100644 index 0000000..69b9c92 --- /dev/null +++ b/scripts/tests/test_null_calibration_failures.py @@ -0,0 +1,89 @@ +# SPDX-License-Identifier: MPL-2.0 +"""CSV parsing and annotation escaping for null-calibration package failures.""" +import csv +from pathlib import Path +import sys +import tempfile +import unittest + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT / "bench/null_calibration")) +from report_failures import LEGACY_FIELDS, RESULT_FIELDS, failure_annotations + + +class FailureAnnotationTests(unittest.TestCase): + def write_rows(self, path, fields, rows): + with path.open("w", encoding="utf-8", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=fields, lineterminator="\n") + writer.writeheader() + writer.writerows( + {field: row.get(field) for field in fields} for row in rows + ) + + def row(self, *, test, status, note=None): + return { + "truth": "null", + "shape": "original", + "offset": "tss", + "n_per": 4, + "taxa": 40, + "sim": 1, + "taxon": 1, + "true_log2fc": 0, + "test": test, + "status": status, + "pvalue": "NA" if status == "failed" else 0.01, + "note": note, + } + + def test_distinct_failures_are_csv_aware_deduplicated_and_escaped(self): + note_one = 'comma, quote " and newline\ncontinued' + note_two = "second 50% and carriage\rreturn\ncontinued" + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "pvalues_smoke.csv" + self.write_rows( + path, + RESULT_FIELDS, + [ + self.row(test="duplicate_test", status="failed", note=note_one), + self.row(test="duplicate_test", status="failed", note=note_one), + self.row(test="duplicate_test", status="failed", note=note_two), + self.row(test="missing_note", status="failed", note="NA"), + self.row(test="ok_test", status="ok"), + ], + ) + cell = "null_original_tss_n4_t40_smoke" + self.assertEqual( + failure_annotations(path, cell), + [ + f"::warning title={cell}: a requested test failed::duplicate_test :: " + 'comma, quote " and newline%0Acontinued', + f"::warning title={cell}: a requested test failed::duplicate_test :: " + "second 50%25 and carriage%0Dreturn%0Acontinued", + ], + ) + + def test_legacy_header_without_notes_and_missing_note_are_ignored(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "legacy.csv" + self.write_rows( + path, + LEGACY_FIELDS, + [self.row(test="old_test", status="failed")], + ) + self.assertEqual(failure_annotations(path, "legacy_cell"), []) + + def test_extra_csv_field_is_rejected(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory) / "malformed.csv" + path.write_text( + ",".join(RESULT_FIELDS) + "\n" + + "null,original,tss,4,40,1,1,0,edger_ql,failed,NA,reason,extra\n", + encoding="utf-8", + ) + with self.assertRaisesRegex(ValueError, "malformed CSV record"): + failure_annotations(path, "bad_cell") + + +if __name__ == "__main__": + unittest.main() diff --git a/test/unit/test_null_calibration_csv.R b/test/unit/test_null_calibration_csv.R new file mode 100644 index 0000000..6f9f707 --- /dev/null +++ b/test/unit/test_null_calibration_csv.R @@ -0,0 +1,109 @@ +# Regression checks for the null-calibration CSV writer, resumer and warning readers. +# Run from the repository root: Rscript --vanilla test/unit/test_null_calibration_csv.R +expressions <- parse("bench/null_calibration/calibration.R") +wanted <- c("RESULT_COLUMNS", "report_cell_failures", "validate_resume_cell", + "write_cell_rows") +loaded <- character() +for (expr in expressions) { + if (is.call(expr) && identical(expr[[1]], as.name("<-")) && + is.symbol(expr[[2]]) && as.character(expr[[2]]) %in% wanted) { + eval(expr) + loaded <- c(loaded, as.character(expr[[2]])) + } +} +stopifnot(setequal(loaded, wanted)) + +cell <- "null_original_tss_n4_t40_smoke" +double_quote <- intToUtf8(34) +note_multiline <- paste0("comma, quote ", double_quote, " and newline", intToUtf8(10), "continued") +note_control <- paste0("second reason 50% and newline", intToUtf8(10), "continued") +rows <- data.frame( + truth = rep("null", 5), shape = rep("original", 5), offset = rep("tss", 5), + n_per = rep(4L, 5), taxa = rep(40L, 5), sim = rep(1L, 5), + taxon = seq_len(5), true_log2fc = rep(0, 5), + test = c("duplicate_test", "duplicate_test", "duplicate_test", "missing_note", "ok_test"), + status = c("failed", "failed", "failed", "failed", "ok"), + pvalue = c(NA_real_, NA_real_, NA_real_, NA_real_, 0.01), + note = c(note_multiline, note_multiline, note_control, NA_character_, NA_character_), + stringsAsFactors = FALSE +) +stopifnot(identical(names(rows), RESULT_COLUMNS)) + +csv_path <- tempfile(fileext = ".csv") +provenance_path <- tempfile(fileext = ".txt") +legacy_path <- tempfile(fileext = ".csv") +bad_csv_path <- tempfile(fileext = ".csv") +write_cell_rows(rows, csv_path, header_done = FALSE) +write_cell_rows(rows[1, , drop = FALSE], csv_path, header_done = TRUE) +round_trip <- read.csv(csv_path, stringsAsFactors = FALSE, check.names = FALSE, fill = FALSE) +stopifnot(nrow(round_trip) == nrow(rows) + 1L, + identical(names(round_trip), RESULT_COLUMNS), + identical(round_trip$note[1:3], c(note_multiline, note_multiline, note_control)), + is.na(round_trip$note[4]), round_trip$status[5] == "ok") + +# R's per-cell diagnostic reader must retain special characters, deduplicate +# repeated reasons, and still warn when a failed row has no note. +r_messages <- capture.output(report_cell_failures(csv_path, cell)) +r_log <- paste(r_messages, collapse = "\n") +stopifnot(sum(grepl(paste0("WARNING ", cell, ": test duplicate_test"), + r_messages, fixed = TRUE)) == 1L, + grepl(note_multiline, r_log, fixed = TRUE), + grepl(paste0("WARNING ", cell, ": test missing_note produced no p-value in this cell"), + r_log, fixed = TRUE), + !grepl("ok_test", r_log, fixed = TRUE)) +file.create(bad_csv_path) +read_error <- tryCatch({ + report_cell_failures(bad_csv_path, cell) + "" +}, error = function(e) conditionMessage(e)) +stopifnot(grepl("unreadable result CSV", read_error, fixed = TRUE)) + +# The workflow calls this same CSV-aware Python reader. It emits distinct +# reasons once, skips NA notes and escapes annotation control characters. +python <- Sys.which("python3") +stopifnot(nzchar(python)) +reader_script <- "bench/null_calibration/report_failures.py" +annotations <- suppressWarnings(system2( + python, + args = c(shQuote(reader_script), shQuote(csv_path), shQuote(cell)), + stdout = TRUE, stderr = TRUE +)) +reader_status <- attr(annotations, "status") +stopifnot(is.null(reader_status) || identical(reader_status, 0L)) +expected_annotations <- c( + paste0("::warning title=", cell, + ": a requested test failed::duplicate_test :: comma, quote ", double_quote, + " and newline%0Acontinued"), + paste0("::warning title=", cell, + ": a requested test failed::duplicate_test :: second reason 50%25 and newline%0Acontinued") +) +stopifnot(identical(unname(annotations), expected_annotations)) + +# Refuse legacy 11-column CSVs and mismatched provenance before any append. +legacy <- rows[, RESULT_COLUMNS[-length(RESULT_COLUMNS)], drop = FALSE] +write.table(legacy, legacy_path, sep = ",", row.names = FALSE, quote = TRUE, + qmethod = "double") +writeLines("cell key: 123", provenance_path) +legacy_error <- tryCatch({ + validate_resume_cell(legacy_path, provenance_path, cell, 123) + "" +}, error = function(e) conditionMessage(e)) +stopifnot(grepl("incompatible CSV schema", legacy_error, fixed = TRUE)) + +stopifnot(identical(validate_resume_cell(csv_path, provenance_path, cell, 123), 1L)) +writeLines("cell key: 999", provenance_path) +key_error <- tryCatch({ + validate_resume_cell(csv_path, provenance_path, cell, 123) + "" +}, error = function(e) conditionMessage(e)) +stopifnot(grepl("mismatched cell key", key_error, fixed = TRUE)) +unlink(provenance_path) +provenance_error <- tryCatch({ + validate_resume_cell(csv_path, provenance_path, cell, 123) + "" +}, error = function(e) conditionMessage(e)) +stopifnot(grepl("no provenance file", provenance_error, fixed = TRUE)) +stopifnot(identical(validate_resume_cell(tempfile(), provenance_path, cell, 123), integer(0))) +unlink(c(csv_path, provenance_path, legacy_path, bad_csv_path)) + +cat("Null-calibration CSV round-trip, warning reader and resume checks passed\n")