From a8c0acd231e618219b18d4e20d88e30e8f7b63ba Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:14:52 +0000 Subject: [PATCH 1/9] Audit numerical boundaries and correct inclusive ALDEx rank tails Add seeded core/statistical fuzz, exhaustive rank distributions, independent R checks and a finite-tail Agda identity with a rejection control. Keep pristine-reference defect controls explicit. --- .github/workflows/ci.yml | 19 +++- docs/ALDEx2-SPEC.md | 41 +++++-- ext/CompositionalDARCallExt.jl | 11 +- proofs/PROOF-STATUS.md | 11 ++ proofs/agda/CompositionalDA/All.agda | 1 + .../agda/CompositionalDA/DiscreteTails.agda | 27 +++++ proofs/agda/reject/LostAtom.agda | 9 ++ src/Core/bh.jl | 4 +- src/Core/composition.jl | 33 +++++- src/Core/zeros.jl | 37 ++++-- src/Methods/ALDEx2.jl | 73 +++++++++--- src/Methods/ANCOMBC.jl | 4 + test/fuzz_ancombc.jl | 53 +++++++++ test/fuzz_core.jl | 107 ++++++++++++++++++ test/fuzz_statistics.jl | 90 +++++++++++++++ test/oracle/adapter_boundaries.jl | 23 ++++ test/oracle/aldex2_tail_control.jl | 18 +++ test/oracle/ancombc_oracle.jl | 1 + test/oracle/runoracle.jl | 16 +++ test/runtests.jl | 3 + 20 files changed, 531 insertions(+), 50 deletions(-) create mode 100644 proofs/agda/CompositionalDA/DiscreteTails.agda create mode 100644 proofs/agda/reject/LostAtom.agda create mode 100644 test/fuzz_ancombc.jl create mode 100644 test/fuzz_core.jl create mode 100644 test/fuzz_statistics.jl create mode 100644 test/oracle/adapter_boundaries.jl create mode 100644 test/oracle/aldex2_tail_control.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 09dc0ed..61e18d7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -3,7 +3,7 @@ name: CI on: push: - branches: [main, 'fix/issue-23-reference-pins'] + branches: [main, 'fix/issue-23-reference-pins', 'audit/issue-23-fuzz-correctness'] tags: ['v*'] pull_request: permissions: @@ -95,6 +95,7 @@ jobs: if (!requireNamespace(p, quietly = TRUE)) BiocManager::install(p, ask = FALSE, update = FALSE) stopifnot(packageVersion(p) == Sys.getenv("REFERENCE_VERSION")) if (p == "ANCOMBC") stopifnot(packageVersion("CVXR") == "1.0.15") + if (p == "ALDEx2" && !requireNamespace("testthat", quietly=TRUE)) install.packages("testthat") cat(p, as.character(packageVersion(p)), "\n")' - uses: julia-actions/setup-julia@4c0cb0fce8556fdb04a90347310e5db8b1f98fb9 # v2 with: @@ -114,11 +115,23 @@ jobs: for l in ~/.julia/scratchspaces/*/*/build.log; do echo "== $l"; cat "$l"; done julia --project=test/oracle -e 'using RCall; println(rcopy(R"R.version.string"))' julia --project=test/oracle -e 'using CompositionalDA, RCall; @assert Base.get_extension(CompositionalDA, :CompositionalDARCallExt) !== nothing; println("extension loaded")' - - name: Native vs reference agreement + - name: Pristine ALDEx2 discrete-tail defect control if: matrix.reference == 'ALDEx2' run: | export R_HOME="$(R RHOME)" - julia --project=test/oracle test/oracle/${{ matrix.script }} + julia --project=test/oracle test/oracle/aldex2_tail_control.jl + - name: ALDEx2 conformance with the pinned inclusive-tail repair + if: always() && matrix.reference == 'ALDEx2' + run: | + mkdir -p /tmp/aldex-tail-repair + for source in stats.fast.R clr_ttest.r; do + curl --fail --location --retry 3 "https://raw.githubusercontent.com/hyperpolymath/ALDEx_bioc/c1217796f551373cfd942d5e272cd56e05ce77a9/R/$source" -o "/tmp/aldex-tail-repair/$source" + done + for source in test-directional-tails.R test-aldex.ttest.R; do + curl --fail --location --retry 3 "https://raw.githubusercontent.com/hyperpolymath/ALDEx_bioc/c1217796f551373cfd942d5e272cd56e05ce77a9/tests/testthat/$source" -o "/tmp/aldex-tail-repair/$source" + done + export R_HOME="$(R RHOME)" + ALDEX2_REPAIR_DIR=/tmp/aldex-tail-repair julia --project=test/oracle test/oracle/runoracle.jl - name: Pristine ANCOMBC known-defect control (three expected failures) if: matrix.reference == 'ANCOMBC' run: | diff --git a/docs/ALDEx2-SPEC.md b/docs/ALDEx2-SPEC.md index dbe4332..93d8ce3 100644 --- a/docs/ALDEx2-SPEC.md +++ b/docs/ALDEx2-SPEC.md @@ -36,6 +36,18 @@ agent that was barred from the reference source, R and web search wrote `src/Methods/ALDEx2.jl` from this document (spec committed alone as `87a6316`, before the code). +### Correctness rereview, 2026-11-07 + +The historical clean-room implementation above describes the original work. +This rereview inspected the upstream R source to diagnose and repair the +directional-tail defect there. The Julia correction is derived from the +finite-distribution CDF identity and extends the existing native recurrence; +no R source was copied or translated into the MPL implementation. It intentionally +differs from pristine ALDEx2 1.42.0 at this defect. CI first demonstrates the +pristine counterexample, then explicitly loads the pinned proposed upstream +repair for conformance. A green repaired-reference job is not pristine-release +equivalence or a proof of posterior calibration/FDR control. + ## Notation ALDEx2's CLR instances and descriptive columns (`diff_btw`, `diff_win`, @@ -67,10 +79,12 @@ number. - The property test P6 calls the function without this flag and gets the default. -Test group 1 only matters for the exact Wilcoxon p-value, because that -distribution is discrete. The rule is part of the reference behaviour: on the -soil dataset, whose first sample is in `L2`, using `L1` instead moves `wi_pvalue` -by up to 0.0097. +Both directional tails now include the observed statistic. Swapping test +groups swaps the two tails and leaves the final minimum unchanged on identical +instances. The earlier 1.42.0 compatibility rule instead used `1 - upper` for +the discrete lower tail; that drops the observed atom and can depend on which +condition occurs first. The exhaustive two-versus-two allocation test detects +this defect: an extreme allocation must report 1/3 in either direction. ## Per-instance one-sided p-values @@ -85,18 +99,23 @@ For each instance `k` and feature `j`, let `a` be the values of test group 1 and distribution of the Mann–Whitney statistic for sizes `(n1, n2)`. - Otherwise use the normal approximation **without continuity correction**: `pR = 1 − Φ((U − n1·n2/2) / sqrt(n1·n2·(n1+n2+1)/12))`. - - CLR instances are continuous, so ties have probability zero. If any do - occur, count a tie as ½ in `U`; the result in that case is not - oracle-checked. + - Floating-point ties can occur. Count each cross-group tie as ½ in `U` + and use the normal approximation with variance + `n1*n2/12 * (N+1 - sum(t^3-t)/(N*(N-1)))`, where the `t` are pooled tie + block sizes. A zero variance is refused explicitly. + - Compute the inclusive lower tail `P(U′ ≤ U)` separately. For the + symmetric exact distribution it is the upper tail at `n1*n2-U`. + Normal and Welch tails are evaluated directly in each direction to avoid + subtractive cancellation at small probabilities. ## Consistent-sign expected p-values -For a per-instance upper-tail p-value `p` (either `pW` or `pR`), define the two -doubled directional values: +For separately evaluated per-instance inclusive tails `upper` and `lower`, +define the two doubled directional values: ``` -g = min(1, 2p) # "group 1 higher" -l = min(1, 2(1 − p)) # "group 1 lower" (note: 1 − p, not the lower tail) +g = min(1, 2 * upper) # "group 1 higher" +l = min(1, 2 * lower) # "group 1 lower", including the observed atom ``` Then, per test: diff --git a/ext/CompositionalDARCallExt.jl b/ext/CompositionalDARCallExt.jl index 9e0c63b..a2bedef 100644 --- a/ext/CompositionalDARCallExt.jl +++ b/ext/CompositionalDARCallExt.jl @@ -11,7 +11,7 @@ module CompositionalDARCallExt using CompositionalDA using CompositionalDA: Provenance, DAResult, validate_counts, CompositionError, - validate_two_groups, validate_mc + validate_two_groups, validate_mc, float_count_matrix using DataFrames using RCall @@ -31,8 +31,8 @@ function CompositionalDA.aldex2_reference(counts::AbstractMatrix{<:Integer}, fea levels = design.levels keep = vec(any(>(0), counts; dims = 1)) count(keep) >= 2 || throw(CompositionError("ALDEx2 needs at least two nonzero features")) + reads = permutedims(float_count_matrix(counts[:, keep])) # R integers are only 32-bit df = runner() do - reads = permutedims(Float64.(counts[:, keep])) # R integers are only 32-bit @rput reads conds = [g == levels[1] ? "0" : "1" for g in design.groups]; lv = ["0", "1"] @rput conds lv @@ -73,9 +73,14 @@ function CompositionalDA.ancombc_reference(counts::AbstractMatrix{<:Integer}, fe 0 <= prevalence < 1 || throw(ArgumentError("prevalence must be in [0,1)")) isfinite(pseudocount) && pseudocount > 0 || throw(ArgumentError("pseudocount must be finite and positive")) isfinite(tolerance) && 0 < tolerance < 1 || throw(ArgumentError("tolerance must be in (0,1)")) + pseudocount = Float64(pseudocount) + tolerance = Float64(tolerance) + isfinite(pseudocount) && pseudocount > 0 || throw(ArgumentError("pseudocount must be representable as a positive finite Float64")) + 0 < tolerance < 1 || throw(ArgumentError("tolerance must be representable as a Float64 strictly between 0 and 1")) max_iterations > 0 || throw(ArgumentError("max_iterations must be positive")) + negative_lower_bound && !structural_zeros && throw(ArgumentError("negative_lower_bound requires structural_zeros")) + reads = permutedims(float_count_matrix(counts)) table, version = runner() do - reads = permutedims(Float64.(counts)) fnames = design.features conds = [g == design.levels[1] ? "reference" : "contrast" for g in design.groups] @rput reads fnames conds prevalence pseudocount structural_zeros negative_lower_bound conserve tolerance max_iterations diff --git a/proofs/PROOF-STATUS.md b/proofs/PROOF-STATUS.md index f90dc51..d4c8a4e 100644 --- a/proofs/PROOF-STATUS.md +++ b/proofs/PROOF-STATUS.md @@ -41,6 +41,17 @@ The exact commands, exit codes and counts of the most recent run are in ## What is proved +### Finite directional tails (correctness rereview) + +`CompositionalDA.DiscreteTails` proves on natural-number outcome counts that +subtracting an inclusive upper-tail count from the total leaves the strict +lower count, and adding the observed atom restores the inclusive lower count. +The six-outcome, one-atom counterexample is checked by normalization; the +`LostAtom` rejection control must fail for claiming the strict count is one. +Exhaustive Julia tests connect this identity to the implemented rank CDF for +all group sizes two through four. This is not a proof of arbitrary-size +floating-point execution, posterior calibration, or FDR control. + ### Additional method identities (v0.3 development) `CompositionalDA.PosteriorP` proves over exact rationals that the sum of diff --git a/proofs/agda/CompositionalDA/All.agda b/proofs/agda/CompositionalDA/All.agda index 0515b43..c2830dd 100644 --- a/proofs/agda/CompositionalDA/All.agda +++ b/proofs/agda/CompositionalDA/All.agda @@ -18,3 +18,4 @@ import CompositionalDA.BenjaminiHochberg import CompositionalDA.BenjaminiHochberg.StepUp import CompositionalDA.PosteriorP import CompositionalDA.BiasModel +import CompositionalDA.DiscreteTails diff --git a/proofs/agda/CompositionalDA/DiscreteTails.agda b/proofs/agda/CompositionalDA/DiscreteTails.agda new file mode 100644 index 0000000..76eca22 --- /dev/null +++ b/proofs/agda/CompositionalDA/DiscreteTails.agda @@ -0,0 +1,27 @@ +-- SPDX-License-Identifier: MPL-2.0 +-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +{-# OPTIONS --safe --without-K #-} +module CompositionalDA.DiscreteTails where + +open import Data.Nat using (ℕ; zero; suc; _+_; _∸_) +open import Relation.Binary.PropositionalEquality using (_≡_; refl; cong) + +-- A finite outcome space is partitioned into below, observed, and above. +-- Subtracting the inclusive upper count from the total leaves only below. +cancel-prefix : ∀ upper below → (upper + below) ∸ upper ≡ below +cancel-prefix zero below = refl +cancel-prefix (suc upper) below = cancel-prefix upper below + +inclusive-lower : ∀ below observed above → + (((observed + above) + below) ∸ (observed + above)) + observed ≡ + below + observed +inclusive-lower below observed above = + cong (_+ observed) (cancel-prefix (observed + above) below) + +-- The extreme allocation in a two-versus-two rank test has one observed +-- outcome among six. Complementing its upper tail loses that one outcome. +lost-atom : (6 ∸ 6) + 1 ≡ 1 +lost-atom = refl + +strict-lower : 6 ∸ 6 ≡ 0 +strict-lower = refl diff --git a/proofs/agda/reject/LostAtom.agda b/proofs/agda/reject/LostAtom.agda new file mode 100644 index 0000000..2c5bc5c --- /dev/null +++ b/proofs/agda/reject/LostAtom.agda @@ -0,0 +1,9 @@ +-- SPDX-License-Identifier: MPL-2.0 +-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +-- EXPECT: 0 != 1 +{-# OPTIONS --safe --without-K #-} +module LostAtom where +open import Data.Nat using (_∸_) +open import Relation.Binary.PropositionalEquality using (_≡_; refl) +false-inclusive-lower : 6 ∸ 6 ≡ 1 +false-inclusive-lower = refl diff --git a/src/Core/bh.jl b/src/Core/bh.jl index 8d1cf0c..9fe0529 100644 --- a/src/Core/bh.jl +++ b/src/Core/bh.jl @@ -14,8 +14,10 @@ function bh_adjust(p::AbstractVector{<:Real}) idx = findall(!isnan, p) m = length(idx) m == 0 && return out + # Validate before conversion: an invalid high-precision value just above + # one can otherwise round to 1.0 and be silently accepted. + all(v -> 0 <= v <= 1, p[idx]) || throw(ArgumentError("p-values must lie in [0, 1]")) q = Float64.(p[idx]) - all(v -> 0 <= v <= 1, q) || throw(ArgumentError("p-values must lie in [0, 1]")) ord = sortperm(q; rev = true) running = 1.0 adj = similar(q) diff --git a/src/Core/composition.jl b/src/Core/composition.jl index e601fb5..987b347 100644 --- a/src/Core/composition.jl +++ b/src/Core/composition.jl @@ -24,6 +24,13 @@ end validate_counts(::AbstractMatrix) = throw(CompositionError("counts must be an integer matrix (raw counts, not normalised values)")) +function float_count_matrix(counts) + validate_counts(counts) + out = Float64.(counts) + all(isfinite, out) || throw(CompositionError("counts are not representable as finite Float64 values")) + return out +end + # Both native methods and their reference adapters share this boundary. In # particular, converting missing labels to the string "missing" is forbidden. function validate_two_groups(counts, features, groups) @@ -33,6 +40,10 @@ function validate_two_groups(counts, features, groups) length(groups) == n || throw(CompositionError("group labels do not match samples")) any(ismissing, features) && throw(CompositionError("feature names cannot be missing")) any(ismissing, groups) && throw(CompositionError("group labels cannot be missing")) + any(isnothing, features) && throw(CompositionError("feature names cannot be nothing")) + any(isnothing, groups) && throw(CompositionError("group labels cannot be nothing")) + any(g -> g isa Real && !isfinite(g), groups) && + throw(CompositionError("numeric group labels must be finite")) f = string.(features) g = string.(groups) any(isempty ∘ strip, f) && throw(CompositionError("feature names cannot be empty")) @@ -46,8 +57,14 @@ function validate_two_groups(counts, features, groups) levels = string.(original_levels) allunique(levels) || throw(CompositionError("distinct group labels must have distinct display names")) length(levels) == 2 || throw(CompositionError("exactly two groups are required")) - i1, i2 = (findall(==(level), groups) for level in original_levels) + # unique uses isequal, including its treatment of signed zero and mixed + # numeric types. Membership and the R labels must use that same partition. + i1, i2 = (findall(isequal(level), groups) for level in original_levels) min(length(i1), length(i2)) >= 2 || throw(CompositionError("each group needs at least two samples")) + length(i1) + length(i2) == n && isempty(intersect(i1, i2)) || + throw(CompositionError("group labels do not define a disjoint sample partition")) + g[i1] .= levels[1] + g[i2] .= levels[2] return (; features = f, groups = g, levels, i1, i2) end @@ -60,12 +77,19 @@ end """ closure(x) -> Matrix{Float64} -Divide each row of a strictly positive samples × features matrix by its row sum, +Divide each row of a nonempty strictly positive samples × features matrix by its row sum, so that each row sums to 1. """ function closure(x::AbstractMatrix{<:Real}) + all(>(0), size(x)) || throw(CompositionError("closure needs a nonempty matrix")) all(v -> isfinite(v) && v > 0, x) || throw(CompositionError("closure needs strictly positive finite parts")) - return x ./ sum(x; dims = 2) + # Scale before summing: both integer totals and finite floating-point + # totals can overflow even though every closed proportion is representable. + scaled = Float64.(x ./ maximum(x; dims = 2)) + out = scaled ./ sum(scaled; dims = 2) + all(v -> isfinite(v) && v > 0, out) || + throw(CompositionError("closed proportions are not representable as positive Float64 values")) + return out end """ @@ -81,6 +105,7 @@ Proved properties (`proofs/PROOF-STATUS.md`, P1): - relabelling features permutes the columns in the same way. """ function clr(x::AbstractMatrix{<:Real}) + all(>(0), size(x)) || throw(CompositionError("clr needs a nonempty matrix")) all(v -> isfinite(v) && v > 0, x) || throw(CompositionError("clr needs strictly positive finite parts; apply a ZeroPolicy first")) l = log.(x) return l .- mean(l; dims = 2) @@ -92,4 +117,4 @@ end Validate raw counts, apply the zero policy, then take the centred log-ratio. """ clr(policy::ZeroPolicy, counts::AbstractMatrix{<:Integer}) = - clr(replace_zeros(policy, validate_counts(counts))) + clr(replace_zeros(policy, counts)) diff --git a/src/Core/zeros.jl b/src/Core/zeros.jl index f182b0c..182174f 100644 --- a/src/Core/zeros.jl +++ b/src/Core/zeros.jl @@ -31,7 +31,9 @@ struct Pseudocount <: ZeroPolicy α::Float64 function Pseudocount(α::Real = 0.5) (isfinite(α) && α > 0) || throw(ArgumentError("pseudocount must be finite and > 0, got $α")) - new(Float64(α)) + stored = Float64(α) + isfinite(stored) && stored > 0 || throw(ArgumentError("pseudocount must be representable as a positive finite Float64")) + new(stored) end end @@ -47,39 +49,50 @@ struct MultiplicativeReplacement <: ZeroPolicy δ::Float64 function MultiplicativeReplacement(δ::Real = 0.65) (isfinite(δ) && 0 < δ < 1) || throw(ArgumentError("δ must be in (0, 1), got $δ")) - new(Float64(δ)) + stored = Float64(δ) + 0 < stored < 1 || throw(ArgumentError("δ must be representable as a Float64 strictly between 0 and 1")) + new(stored) end end """ replace_zeros(policy, counts) -> Matrix{Float64} -Apply a zero policy to a validated samples × features count matrix and return a +Validate a samples × features count matrix, apply a zero policy and return a strictly positive matrix. For `MultiplicativeReplacement` the result is closed: each row sums to 1. For the other policies, rows keep their count scale. """ function replace_zeros end -replace_zeros(::RefuseZeros, counts::AbstractMatrix{<:Integer}) = - any(iszero, counts) ? - throw(CompositionError("zeros present and the zero policy is RefuseZeros")) : - Float64.(counts) +function replace_zeros(::RefuseZeros, counts::AbstractMatrix{<:Integer}) + out = float_count_matrix(counts) + any(iszero, out) && throw(CompositionError("zeros present and the zero policy is RefuseZeros")) + return out +end -replace_zeros(p::Pseudocount, counts::AbstractMatrix{<:Integer}) = Float64.(counts) .+ p.α +function replace_zeros(p::Pseudocount, counts::AbstractMatrix{<:Integer}) + out = float_count_matrix(counts) .+ p.α + all(v -> isfinite(v) && v > 0, out) || throw(CompositionError("pseudocount addition is not representable as positive finite Float64 values")) + return out +end function replace_zeros(p::MultiplicativeReplacement, counts::AbstractMatrix{<:Integer}) + validate_counts(counts) out = Matrix{Float64}(undef, size(counts)) for i in axes(counts, 1) row = @view counts[i, :] - total = sum(row) - total > 0 || throw(CompositionError("sample $i has zero total count")) - r = p.δ / total + scale = maximum(row) + scaled = Float64.(row ./ scale) + scaled_total = sum(scaled) nz = count(iszero, row) + r = nz == 0 ? 0.0 : Float64(p.δ / scale / scaled_total) + nz == 0 || r > 0 || throw(CompositionError("sample $i: zero replacement underflows Float64")) shrink = 1 - nz * r shrink > 0 || throw(CompositionError("sample $i: replacement δ too large for $nz zeros")) for j in axes(counts, 2) - out[i, j] = iszero(row[j]) ? r : shrink * row[j] / total + out[i, j] = iszero(row[j]) ? r : shrink * scaled[j] / scaled_total end end + all(v -> isfinite(v) && v > 0, out) || throw(CompositionError("replacement proportions are not representable as positive finite Float64 values")) return out end diff --git a/src/Methods/ALDEx2.jl b/src/Methods/ALDEx2.jl index 1d9c0be..ed2aa2e 100644 --- a/src/Methods/ALDEx2.jl +++ b/src/Methods/ALDEx2.jl @@ -22,14 +22,16 @@ function dirichlet_clr_instances(rng::AbstractRNG, counts::AbstractMatrix{<:Inte validate_counts(counts) mc_samples >= 2 || throw(ArgumentError("mc_samples must be at least 2")) isfinite(prior) && prior > 0 || throw(ArgumentError("prior must be finite and positive")) + stored_prior = Float64(prior) + isfinite(stored_prior) && stored_prior > 0 || throw(ArgumentError("prior must be representable as a positive finite Float64")) + shapes = float_count_matrix(counts) .+ stored_prior + all(isfinite, shapes) || throw(CompositionError("Dirichlet shape is not representable in Float64")) n, d = size(counts) out = Array{Float64,3}(undef, n, d, mc_samples) g = Vector{Float64}(undef, d) for i in 1:n, k in 1:mc_samples for j in 1:d - shape = Float64(counts[i, j]) + Float64(prior) - isfinite(shape) || throw(CompositionError("Dirichlet shape is not representable in Float64")) - g[j] = rand(rng, Gamma(shape, 1.0)) + g[j] = rand(rng, Gamma(shapes[i, j], 1.0)) end # Flooring a zero draw would silently change the posterior. Refuse an # unrepresentable draw instead. ALDEx2 reports CLR values in log2 units. @@ -41,12 +43,23 @@ function dirichlet_clr_instances(rng::AbstractRNG, counts::AbstractMatrix{<:Inte end """ - welch_upper_p(a, b) -> Float64 + welch_tails(a, b) -> (upper, lower) -One-sided Welch t-test p-value `P(T ≥ t)` for `mean(a) > mean(b)`, with -Welch–Satterthwaite degrees of freedom. +Both one-sided Welch t-test p-values, evaluated separately with +Welch–Satterthwaite degrees of freedom to preserve small probabilities. """ -welch_upper_p(a, b) = pvalue(UnequalVarianceTTest(a, b); tail = :right) +function welch_tails(a, b) + variance = var(a) / length(a) + var(b) / length(b) + isfinite(variance) && variance > 0 || throw(CompositionError("Welch test has zero or non-finite sampling variance")) + try + test = UnequalVarianceTTest(a, b) + return pvalue(test; tail=:right), pvalue(test; tail=:left) + catch error + error isa DomainError && throw(CompositionError("Welch test is not estimable at the available numerical precision")) + rethrow() + end +end +welch_upper_p(a, b) = first(welch_tails(a, b)) """ mann_whitney_u(a, b) -> Float64 @@ -85,18 +98,40 @@ function mann_whitney_upper_tail(n1::Integer, n2::Integer) end """ - wilcoxon_upper_p(a, b, tail) -> Float64 + wilcoxon_tails(a, b, tail) -> (upper, lower) One-sided Wilcoxon rank-sum p-value for `a` above `b`. With `tail` from [`mann_whitney_upper_tail`](@ref) it is exact, `P(U' ≥ U)`; with `tail === nothing` it is the normal approximation without continuity correction. """ -function wilcoxon_upper_p(a, b, tail) +function wilcoxon_tails(a, b, tail) u = mann_whitney_u(a, b) - tail === nothing || return tail[Int(ceil(u)) + 1] n1 = length(a); n2 = length(b) - return ccdf(Normal(), (u - n1 * n2 / 2) / sqrt(n1 * n2 * (n1 + n2 + 1) / 12)) + pooled = sort(vcat(a, b)) + all(isfinite, pooled) || throw(CompositionError("Wilcoxon observations must be finite")) + correction = 0.0 + first = 1 + while first <= length(pooled) + last = first + while last < length(pooled) && pooled[last+1] == pooled[first] + last += 1 + end + size = Float64(last - first + 1) + correction += size^3 - size + first = last + 1 + end + if tail !== nothing && correction == 0 + # Both tails include the observed atom. 1-upper is the STRICT lower + # tail and gives direction-dependent, anti-conservative rank p-values. + return tail[Int(u) + 1], tail[n1*n2 - Int(u) + 1] + end + n = n1 + n2 + variance = n1*n2 / 12 * (n+1 - correction / (n*(n-1))) + isfinite(variance) && variance > 0 || throw(CompositionError("Wilcoxon test has zero or non-finite rank variance")) + z = (u - n1*n2/2) / sqrt(variance) + return ccdf(Normal(), z), ccdf(Normal(), -z) end +wilcoxon_upper_p(a, b, tail) = first(wilcoxon_tails(a, b, tail)) """ effect_pairs(rng, n1, n2) -> NamedTuple @@ -152,14 +187,19 @@ function aldex2_from_instances(inst::AbstractArray{<:Real,3}, i1::AbstractVector tail = length(t1) < 50 && length(t2) < 50 ? mann_whitney_upper_tail(length(t1), length(t2)) : nothing acc = Dict(k => zeros(d) for k in (:wg, :wl, :wgq, :wlq, :rg, :rl, :rgq, :rlq)) pw = Vector{Float64}(undef, d); pr = Vector{Float64}(undef, d) + lw = similar(pw); lr = similar(pr) for k in 1:mc for j in 1:d a = inst[t1, j, k]; b = inst[t2, j, k] - pw[j] = welch_upper_p(a, b) - pr[j] = wilcoxon_upper_p(a, b, tail) + pw[j], lw[j] = welch_tails(a, b) + pr[j], lr[j] = wilcoxon_tails(a, b, tail) + all(p -> isfinite(p) && 0 <= p <= 1, (pw[j], lw[j])) || + throw(CompositionError("Welch test is not estimable for feature column $j at Monte Carlo instance $k")) + all(p -> isfinite(p) && 0 <= p <= 1, (pr[j], lr[j])) || + throw(CompositionError("Wilcoxon test is not estimable for feature column $j at Monte Carlo instance $k")) end - for (p, g, l, gq, lq) in ((pw, :wg, :wl, :wgq, :wlq), (pr, :rg, :rl, :rgq, :rlq)) - up = min.(1.0, 2 .* p); dn = min.(1.0, 2 .* (1 .- p)) + for (greater, less, g, l, gq, lq) in ((pw, lw, :wg, :wl, :wgq, :wlq), (pr, lr, :rg, :rl, :rgq, :rlq)) + up = min.(1.0, 2 .* greater); dn = min.(1.0, 2 .* less) acc[g] .+= up; acc[l] .+= dn acc[gq] .+= bh_adjust(up); acc[lq] .+= bh_adjust(dn) end @@ -207,7 +247,8 @@ function aldex2(counts::AbstractMatrix{<:Integer}, features::AbstractVector, insertcols!(table, 1, :feature => design.features[keep]) prov = Provenance(:aldex2, :native; seed, mc_samples, test = "welch+wilcoxon", contrast = "$(levels[2]) vs $(levels[1])", prior = 0.5, - log_base = 2, excluded_all_zero = design.features[.!keep]) + log_base = 2, directional_tails = "inclusive; tie-adjusted normal for tied ranks", + excluded_all_zero = design.features[.!keep]) return DAResult(table, prov) end diff --git a/src/Methods/ANCOMBC.jl b/src/Methods/ANCOMBC.jl index 6138790..5c7dabf 100644 --- a/src/Methods/ANCOMBC.jl +++ b/src/Methods/ANCOMBC.jl @@ -150,6 +150,10 @@ function ancombc(counts::AbstractMatrix{<:Integer}, features::AbstractVector, isfinite(prevalence) && 0 <= prevalence < 1 || throw(ArgumentError("prevalence must be in [0,1)")) isfinite(pseudocount) && pseudocount > 0 || throw(ArgumentError("pseudocount must be finite and positive")) isfinite(tolerance) && 0 < tolerance < 1 || throw(ArgumentError("tolerance must be in (0,1)")) + pseudocount = Float64(pseudocount) + tolerance = Float64(tolerance) + isfinite(pseudocount) && pseudocount > 0 || throw(ArgumentError("pseudocount must be representable as a positive finite Float64")) + 0 < tolerance < 1 || throw(ArgumentError("tolerance must be representable as a Float64 strictly between 0 and 1")) max_iterations > 0 || throw(ArgumentError("max_iterations must be positive")) negative_lower_bound && !structural_zeros && throw(ArgumentError("negative_lower_bound requires structural_zeros")) present = counts .> 0 diff --git a/test/fuzz_ancombc.jl b/test/fuzz_ancombc.jl new file mode 100644 index 0000000..7fcca88 --- /dev/null +++ b/test/fuzz_ancombc.jl @@ -0,0 +1,53 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +@testset "Seeded ANCOM-BC metamorphic and likelihood checks" begin + for seed in (23, 20261107, 0xcafe) + rng = StableRNG(seed) + n, d = 20, 40 + x = round.(Int, permutedims(range(1000,5000; length=d)) .* exp.(0.25 .* randn(rng,n,d))) + x[11:20,1:3] .*= 8 + x[11:20,4:5] .= max.(1,x[11:20,4:5] .÷ 4) + f = string.(1:d); g = repeat([1,2]; inner=10) + result = ancombc(x,f,g) + @test all(isfinite,result.table.effect) + @test all(0 .<= result.table.pvalue .<= result.table.qvalue .<= 1) + @test all(result.table.se .> 0) + @test result.table == ancombc(x,f,g).table + perm = randperm(rng,n) + reordered = ancombc(x[perm,:],f,g[perm]) + reversed = ancombc(x,f,3 .- g) + # Scaling counts AND the pseudocount adds a constant to every log + # count. The estimator removes it exactly in real arithmetic. + scaled = ancombc(7 .* x,f,g; pseudocount=7) + for col in (:effect,:se,:pvalue,:qvalue) + @test reordered.table[!,col] ≈ result.table[!,col] atol=1e-7 + @test scaled.table[!,col] ≈ result.table[!,col] atol=1e-7 + @test reversed.table[!,col] ≈ (col == :effect ? -1 : 1) .* result.table[!,col] atol=1e-7 + end + y = log.(x .+ 1) + moments = CompositionalDA.ancombc_moments(y,1:10,11:20) + fit = CompositionalDA.ancombc_bias(moments.difference,moments.variance) + responsibilities, likelihood = CompositionalDA.mixture_responsibilities( + moments.difference,moments.variance,fit.centres,fit.extra,fit.mixing) + @test vec(sum(responsibilities;dims=2)) ≈ ones(d) + @test likelihood ≈ fit.likelihood + @test issorted(fit.centres) + @test sum(fit.mixing) ≈ 1 + @test all(fit.extra .>= 0) && fit.extra[2] == 0 + # Homoscedastic alternative component has an analytic optimum. + for _ in 1:64 + values = randn(rng,8); variance = rand(rng)+0.01 + weights = rand(rng,8); centre = randn(rng) + expected = max(sum(weights .* (values .- centre).^2)/sum(weights)-variance,0) + actual = CompositionalDA.variance_component(values,fill(variance,8),weights,centre) + @test actual ≈ expected atol=2e-7 rtol=2e-7 + end + sparse = hcat(x,vcat(ones(Int,1),zeros(Int,19))) + bounded = ancombc(sparse,vcat(f,"rare"),g; prevalence=0,structural_zeros=true,negative_lower_bound=true) + @test bounded.table.status[end] == "structural_zero_both" + @test isnan(bounded.table.pvalue[end]) && isnan(bounded.table.qvalue[end]) + @test_throws ArgumentError ancombc(x,f,g; negative_lower_bound=true) + @test_throws CompositionError ancombc(fill(big(10)^400,4,4),string.(1:4),[1,1,2,2]) + end +end diff --git a/test/fuzz_core.jl b/test/fuzz_core.jl new file mode 100644 index 0000000..69b31c2 --- /dev/null +++ b/test/fuzz_core.jl @@ -0,0 +1,107 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +@testset "Adversarial core boundaries" begin + @test closure(fill(typemax(Int), 1, 2)) ≈ fill(0.5, 1, 2) + @test closure(fill(floatmax(Float64), 1, 2)) ≈ fill(0.5, 1, 2) + @test closure(fill(big(10)^1000, 1, 2)) ≈ fill(0.5, 1, 2) + for empty in (zeros(0, 2), zeros(1, 0)) + @test_throws CompositionError closure(empty) + @test_throws CompositionError clr(empty) + end + @test_throws CompositionError closure([BigFloat("1e-1000") BigFloat(1)]) + @test replace_zeros(MultiplicativeReplacement(), fill(typemax(UInt64), 1, 2)) ≈ fill(0.5, 1, 2) + @test replace_zeros(MultiplicativeReplacement(), fill(big(10)^1000, 1, 2)) ≈ fill(0.5, 1, 2) + @test_throws CompositionError replace_zeros(MultiplicativeReplacement(), [big(0) big(10)^1000]) + for policy in (RefuseZeros(), Pseudocount(), MultiplicativeReplacement()) + @test_throws CompositionError replace_zeros(policy, [-1 2]) + @test_throws CompositionError replace_zeros(policy, zeros(Int, 1, 2)) + end + @test_throws CompositionError replace_zeros(Pseudocount(1e308), fill(big(10)^308, 1, 2)) + @test_throws CompositionError replace_zeros(RefuseZeros(), fill(big(10)^1000, 1, 2)) + @test_throws ArgumentError Pseudocount(BigFloat("1e-400")) + @test_throws ArgumentError Pseudocount(BigFloat("1e400")) + @test_throws ArgumentError MultiplicativeReplacement(BigFloat(1) - big(2.0)^(-200)) + @test_throws ArgumentError bh_adjust([BigFloat(1) + big(2.0)^(-200)]) + @test_throws ArgumentError bh_adjust([-BigFloat("1e-400")]) + + x = [2 4; 3 5; 4 2; 5 3] + signed = CompositionalDA.validate_two_groups(x, ["a", "b"], [-0.0, -0.0, 0.0, 0.0]) + @test signed.i1 == [1, 2] && signed.i2 == [3, 4] + mixed = CompositionalDA.validate_two_groups(x, ["a", "b"], Any[2, 2.0, 10, 10.0]) + @test mixed.groups == ["2", "2", "10", "10"] + for missing_label in (missing, nothing, NaN, Inf, -Inf) + @test_throws CompositionError CompositionalDA.validate_two_groups(x, ["a", "b"], + Any[missing_label, missing_label, 1, 1]) + end + @test_throws CompositionError CompositionalDA.validate_two_groups(x, ["a", nothing], [1, 1, 2, 2]) +end + +@testset "Seeded count and composition fuzz" begin + for seed in (23, 20261107, 0xcafe), case in 1:128 + @testset "seed=$seed case=$case" begin + rng = StableRNG(UInt64(seed) + UInt64(case) * 1009) + T = (Int8, Int16, Int64, UInt64, BigInt)[mod1(case, 5)] + n, d = rand(rng, 1:8), rand(rng, 2:12) + x = T.(rand(rng, 0:100, n, d)) + if T != BigInt && case % 3 == 0 + x[rand(rng, 1:n), rand(rng, 1:d)] = typemax(T) + end + for i in 1:n + any(>(0), x[i, :]) || (x[i, 1] = one(T)) + end + @test validate_counts(x) === x + positive = max.(x, one(T)) + exact = BigFloat.(positive) + reference = Float64.(exact ./ sum(exact; dims=2)) + @test closure(positive) ≈ reference rtol=2e-15 + scale = exp2(rand(rng, -300:300)) + @test closure(Float64.(positive) .* scale) ≈ reference rtol=2e-15 + @test clr(Float64.(positive) .* scale) ≈ clr(Float64.(positive)) atol=2e-12 + + policy = MultiplicativeReplacement(rand(rng, (0.1, 0.5, 0.65, 0.9))) + totals = sum(BigFloat.(x); dims=2) + replacements = BigFloat(policy.δ) ./ totals + shrinks = 1 .- sum(iszero.(x); dims=2) .* replacements + if any(<=(0), shrinks) + @test_throws CompositionError replace_zeros(policy, x) + else + expected = [Float64(iszero(x[i,j]) ? replacements[i] : + shrinks[i] * BigFloat(x[i,j]) / totals[i]) for i in 1:n, j in 1:d] + result = replace_zeros(policy, x) + @test result ≈ expected rtol=3e-15 + @test all(>(0), result) + @test all(isapprox.(sum(result; dims=2), 1; atol=3e-15)) + end + @test replace_zeros(Pseudocount(0.5), x) == Float64.(x) .+ 0.5 + if any(iszero, x) + @test_throws CompositionError replace_zeros(RefuseZeros(), x) + else + @test replace_zeros(RefuseZeros(), x) == Float64.(x) + end + end + end +end + +@testset "Seeded BH fuzz against a high-precision threshold definition" begin + for seed in (23, 20261107, 0xcafe), case in 1:256 + rng = StableRNG(UInt64(seed) + UInt64(case) * 1009) + p = round.(rand(rng, rand(rng, 0:32)); digits=rand(rng, 1:5)) + if !isempty(p) + p[rand(rng, eachindex(p))] = rand(rng, (0.0, 1.0, NaN)) + end + q = bh_adjust(p) + finite = findall(!isnan, p) + bigp = BigFloat.(p[finite]) + for (k, i) in enumerate(finite) + # Independent definition by all candidate rejection thresholds; + # no sorted cumulative-minimum loop is reused from the implementation. + expected = minimum(vcat(BigFloat[1], + [length(bigp) * t / count(<=(t), bigp) for t in bigp if t >= bigp[k]])) + @test q[i] ≈ Float64(expected) atol=2e-15 + end + @test all(isnan, q[findall(isnan, p)]) + permutation = randperm(rng, length(p)) + @test isequal(bh_adjust(p[permutation]), q[permutation]) + end +end diff --git a/test/fuzz_statistics.jl b/test/fuzz_statistics.jl new file mode 100644 index 0000000..b354c5b --- /dev/null +++ b/test/fuzz_statistics.jl @@ -0,0 +1,90 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +@testset "Exhaustive small-sample rank distributions" begin + for n1 in 2:4, n2 in 2:4 + n = n1+n2 + allocations = [[i for i in 1:n if !iszero(mask & (1 << (i-1)))] + for mask in 0:((1 << n)-1) if count_ones(mask) == n1] + statistics = [sum(a .> b' ) for a in allocations for b in [setdiff(1:n, a)]] + tail = CompositionalDA.mann_whitney_upper_tail(n1, n2) + reported = Float64[] + exact_p = Rational{Int}[] + for (a, u) in zip(allocations, statistics) + b = setdiff(1:n, a) + greater, less = CompositionalDA.wilcoxon_tails(a, b, tail) + expected_greater = count(>=(u), statistics) // length(statistics) + expected_less = count(<=(u), statistics) // length(statistics) + @test greater ≈ Float64(expected_greater) atol=2e-15 + @test less ≈ Float64(expected_less) atol=2e-15 + reversed = CompositionalDA.wilcoxon_tails(b, a, + CompositionalDA.mann_whitney_upper_tail(n2, n1)) + @test greater ≈ reversed[2] atol=2e-15 + @test less ≈ reversed[1] atol=2e-15 + push!(reported, min(1, 2greater, 2less)) + push!(exact_p, min(1, 2expected_greater, 2expected_less)) + end + # Every equiprobable allocation is enumerated: this is an exact finite + # null check, not a Monte Carlo estimate of type-I error. + for alpha in unique(exact_p) + @test count(p -> p <= Float64(alpha)+1e-14, reported) / length(reported) <= Float64(alpha)+1e-14 + end + end + tails = CompositionalDA.wilcoxon_tails([1., 2], [3., 4], + CompositionalDA.mann_whitney_upper_tail(2, 2)) + @test min(1, 2tails[1], 2tails[2]) ≈ 1/3 +end + +@testset "Directional tails, ties and numerical refusals" begin + greater, less = CompositionalDA.welch_tails([1e8, 1e8+1, 1e8+2, 1e8+3], [1., 2, 3, 4]) + @test 0 < greater < 1e-20 + opposite = CompositionalDA.welch_tails([1., 2, 3, 4], [1e8, 1e8+1, 1e8+2, 1e8+3]) + @test opposite[2] == greater + @test less == opposite[1] + @test_throws CompositionError CompositionalDA.welch_tails(ones(3), ones(3)) + @test_throws CompositionError CompositionalDA.wilcoxon_tails(ones(3), ones(3), nothing) + tied = CompositionalDA.wilcoxon_tails([1., 1, 2], [2., 3, 3], + CompositionalDA.mann_whitney_upper_tail(3, 3)) + # Midrank U = 1/2 and tie-adjusted variance = 24/5, without continuity correction. + @test tied[2] ≈ CompositionalDA.ccdf(CompositionalDA.Normal(), 4 / sqrt(24/5)) + @test sum(tied) ≈ 1 + x = [2 4; 3 5; 4 2; 5 3] + @test_throws ArgumentError dirichlet_clr_instances(StableRNG(1), x; prior=BigFloat("1e-400")) + @test_throws ArgumentError dirichlet_clr_instances(StableRNG(1), x; prior=BigFloat("1e400")) + @test_throws CompositionError aldex2(fill(big(10)^300, 4, 2), ["a", "b"], [1, 1, 2, 2]; mc_samples=2) + @test_throws ArgumentError ancombc(x, ["a", "b"], [1, 1, 2, 2]; pseudocount=BigFloat("1e-400")) + @test_throws ArgumentError ancombc(x, ["a", "b"], [1, 1, 2, 2]; tolerance=BigFloat("1e-400")) +end + +@testset "Seeded ALDEx2 input and instance fuzz" begin + for seed in (23, 20261107, 0xcafe), case in 1:16 + @testset "seed=$seed case=$case" begin + rng = StableRNG(UInt64(seed)+UInt64(case)*1009) + n1, n2, d = rand(rng, 2:6), rand(rng, 2:6), rand(rng, 2:10) + x = rand(rng, 1:100, n1+n2, d) + features = ["f$j" for j in 1:d] + groups = vcat(fill(2, n1), fill(10, n2)) + first = aldex2(x, features, groups; seed=case, mc_samples=8) + @test isequal(first.table, aldex2(x, features, groups; seed=case, mc_samples=8).table) + @test first.provenance.parameters[:contrast] == "10 vs 2" + for (p, q) in ((:pvalue, :qvalue), (:wi_pvalue, :wi_qvalue)) + @test all(isfinite, first.table[!, p]) + @test all(0 .<= first.table[!, p] .<= first.table[!, q] .<= 1) + end + inst = dirichlet_clr_instances(StableRNG(case), x; mc_samples=8) + rows1, rows2 = collect(1:n1), collect(n1+1:n1+n2) + direct = CompositionalDA.aldex2_from_instances(inst, rows1, rows2; rng=StableRNG(case)) + reverse = CompositionalDA.aldex2_from_instances(inst, rows1, rows2; + rng=StableRNG(case), first_is_i1=false) + for column in (:pvalue, :qvalue, :wi_pvalue, :wi_qvalue) + @test direct[!, column] ≈ reverse[!, column] atol=2e-14 + end + permutation = randperm(rng, d) + permuted = CompositionalDA.aldex2_from_instances(inst[:, permutation, :], rows1, rows2; + rng=StableRNG(case)) + for column in propertynames(direct) + @test isequal(permuted[!, column], direct[permutation, column]) + end + end + end +end diff --git a/test/oracle/adapter_boundaries.jl b/test/oracle/adapter_boundaries.jl new file mode 100644 index 0000000..8146632 --- /dev/null +++ b/test/oracle/adapter_boundaries.jl @@ -0,0 +1,23 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +@testset "Reference adapter validation precedes the R runner" begin + x = [10 20 30 40; 11 22 29 41; 30 20 10 42; 31 21 11 43] + features = string.(1:4); groups = [1,1,2,2] + calls = Ref(0) + runner = f -> (calls[] += 1; error("invalid input entered R runner")) + for method in (aldex2_reference,ancombc_reference) + for bad in (fill(-1,4,4),zeros(Int,4,4),fill(big(10)^400,4,4)) + @test_throws CompositionError method(bad,features,groups;runner) + end + for labels in ([1,1,2,missing],[NaN,NaN,2,2],[1,1,1,2],[1,1,1,1]) + @test_throws CompositionError method(x,features,labels;runner) + end + @test_throws CompositionError method(x,["a","a","b","c"],groups;runner) + end + @test_throws ArgumentError aldex2_reference(x,features,groups;seed=big(2)^40,runner) + @test_throws ArgumentError aldex2_reference(x,features,groups;mc_samples=1,runner) + @test_throws ArgumentError ancombc_reference(x,features,groups;negative_lower_bound=true,runner) + @test_throws ArgumentError ancombc_reference(x,features,groups;pseudocount=big"1e-1000",runner) + @test_throws ArgumentError ancombc_reference(x,features,groups;tolerance=big"1e-1000",runner) + @test calls[] == 0 +end diff --git a/test/oracle/aldex2_tail_control.jl b/test/oracle/aldex2_tail_control.jl new file mode 100644 index 0000000..bfc7421 --- /dev/null +++ b/test/oracle/aldex2_tail_control.jl @@ -0,0 +1,18 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +using RCall, Test + +@testset "Pristine ALDEx2 discrete-tail defect control" begin + @test rcopy(String,R"as.character(packageVersion('ALDEx2'))") == "1.42.0" + R""" + rank_data <- rbind(1:4,4:1) + upper <- ALDEx2:::wilcox.fast(rank_data,c(1,1,2,2),FALSE) + legacy <- pmin(1,2*upper,2*(1-upper)) + inclusive <- vapply(seq_len(2),function(i) wilcox.test(rank_data[i,1:2], + rank_data[i,3:4],alternative="two.sided")$p.value,0) + """ + @test rcopy(Vector{Float64},R"inclusive") ≈ [1/3,1/3] + @test rcopy(Float64,R"legacy[1]") < 1e-12 + @test rcopy(Float64,R"legacy[2]") ≈ 1/3 + println("KNOWN-DEFECT CONTROL: pristine ALDEx2 1.42.0 drops the observed rank atom; legacy values = ",rcopy(R"legacy")) +end diff --git a/test/oracle/ancombc_oracle.jl b/test/oracle/ancombc_oracle.jl index f5686c7..0c97b7f 100644 --- a/test/oracle/ancombc_oracle.jl +++ b/test/oracle/ancombc_oracle.jl @@ -3,6 +3,7 @@ using CompositionalDA, DataFrames, RCall, Statistics, Test include("datasets.jl") +include("adapter_boundaries.jl") # Diagnostic comparison to an explicitly pinned upstream repair. The ordinary # job leaves this unset and still tests the unmodified Bioconductor release. diff --git a/test/oracle/runoracle.jl b/test/oracle/runoracle.jl index f9c6c51..d58c4d4 100644 --- a/test/oracle/runoracle.jl +++ b/test/oracle/runoracle.jl @@ -11,7 +11,23 @@ using CompositionalDA, DataFrames, RCall, Statistics, Test +if haskey(ENV,"ALDEX2_REPAIR_DIR") + repair_dir = ENV["ALDEX2_REPAIR_DIR"] + @rput repair_dir + R""" + repair_environment <- new.env(parent=asNamespace("ALDEx2")) + sys.source(file.path(repair_dir,"stats.fast.R"),envir=repair_environment) + sys.source(file.path(repair_dir,"clr_ttest.r"),envir=repair_environment) + for (symbol in c("t.fast","wilcox.fast","aldex.ttest")) + assignInNamespace(symbol,repair_environment[[symbol]],ns="ALDEx2") + for (test_file in c("test-directional-tails.R","test-aldex.ttest.R")) + testthat::test_file(file.path(repair_dir,test_file),stop_on_failure=TRUE) + """ + println("REFERENCE OVERRIDE: inclusive directional tails from ",repair_dir) +end + include("datasets.jl") +include("adapter_boundaries.jl") const MC = 512 const SLACK = 2.0 diff --git a/test/runtests.jl b/test/runtests.jl index 901425f..20669bb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,10 +15,13 @@ include("mutants.jl") Aqua.test_all(CompositionalDA) end include("unit_core.jl") + include("fuzz_core.jl") + include("fuzz_statistics.jl") include("property_clr.jl") include("property_bh.jl") include("property_aldex2.jl") include("boundary_aldex2.jl") include("method_ancombc.jl") + include("fuzz_ancombc.jl") include("posterior_aggregation.jl") end From 6c906da6b5b256db3fc3dd197077fa8da2b62e70 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:21:49 +0000 Subject: [PATCH 2/9] Check normal-approximation boundaries and the exact lost-atom rejection Name the Agda negative-control module by its actual path so the gate verifies the intended false equality; compare native tails at 49/50 samples and tied ranks against independent R stats. --- proofs/agda/reject/LostAtom.agda | 2 +- test/oracle/runoracle.jl | 26 ++++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/proofs/agda/reject/LostAtom.agda b/proofs/agda/reject/LostAtom.agda index 2c5bc5c..cafcebf 100644 --- a/proofs/agda/reject/LostAtom.agda +++ b/proofs/agda/reject/LostAtom.agda @@ -2,7 +2,7 @@ -- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) -- EXPECT: 0 != 1 {-# OPTIONS --safe --without-K #-} -module LostAtom where +module reject.LostAtom where open import Data.Nat using (_∸_) open import Relation.Binary.PropositionalEquality using (_≡_; refl) false-inclusive-lower : 6 ∸ 6 ≡ 1 diff --git a/test/oracle/runoracle.jl b/test/oracle/runoracle.jl index d58c4d4..e4c33bc 100644 --- a/test/oracle/runoracle.jl +++ b/test/oracle/runoracle.jl @@ -29,6 +29,32 @@ end include("datasets.jl") include("adapter_boundaries.jl") +@testset "Native directional branches against independent R stats" begin + rng = StableRNG(23) + for (n1,n2,tied) in ((2,3,false),(49,49,false),(49,50,false), + (50,49,false),(50,50,false),(3,4,true),(50,50,true)) + a = tied ? Float64.(rand(rng,1:4,n1)) : randn(rng,n1) + b = tied ? Float64.(rand(rng,1:4,n2)) : randn(rng,n2) .+ 0.5 + @rput a b + R""" + rank_tails <- vapply(c("greater","less"),function(direction) + suppressWarnings(wilcox.test(a,b,correct=FALSE,alternative=direction)$p.value),0) + welch_tails <- vapply(c("greater","less"),function(direction) + t.test(a,b,alternative=direction)$p.value,0) + """ + cache = n1 < 50 && n2 < 50 ? CompositionalDA.mann_whitney_upper_tail(n1,n2) : nothing + @test collect(CompositionalDA.wilcoxon_tails(a,b,cache)) ≈ rcopy(Vector{Float64},R"rank_tails") atol=1e-12 + @test collect(CompositionalDA.welch_tails(a,b)) ≈ rcopy(Vector{Float64},R"welch_tails") atol=1e-12 + end + x,f,g,_ = mock_dataset() + canonical = aldex2_reference(x,f,g;mc_samples=8,seed=23) + numeric_groups = [g[i] == first(g) ? (isodd(i) ? 2 : 2.0) : (isodd(i) ? 10 : 10.0) for i in eachindex(g)] + # The fixture starts in its lower sorted level; equivalent mixed numeric + # labels must not split the R group encoding. + @test first(g) == first(sort(unique(g))) + @test isequal(aldex2_reference(x,f,numeric_groups;mc_samples=8,seed=23).table,canonical.table) +end + const MC = 512 const SLACK = 2.0 const FLOOR = (effect = 0.05, pvalue = 0.01, qvalue = 0.01, From 031b75cfc5335d16127388b9ee34bb80962ca552 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:50:09 +0000 Subject: [PATCH 3/9] Update src/Core/composition.jl Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- src/Core/composition.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Core/composition.jl b/src/Core/composition.jl index 987b347..147f3bd 100644 --- a/src/Core/composition.jl +++ b/src/Core/composition.jl @@ -117,4 +117,4 @@ end Validate raw counts, apply the zero policy, then take the centred log-ratio. """ clr(policy::ZeroPolicy, counts::AbstractMatrix{<:Integer}) = - clr(replace_zeros(policy, counts)) + clr(replace_zeros(policy, validate_counts(counts))) From 38e59234fd3ec29a1ec8905be44da0b9603c6d73 Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Sat, 7 Nov 2026 19:56:21 +0000 Subject: [PATCH 4/9] Preserve positive proportions in narrow-float closure and near-limit zero replacement --- src/Core/composition.jl | 5 ++++- src/Core/zeros.jl | 6 ++++++ test/fuzz_core.jl | 21 +++++++++++++++++++++ 3 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/Core/composition.jl b/src/Core/composition.jl index 147f3bd..e3f7569 100644 --- a/src/Core/composition.jl +++ b/src/Core/composition.jl @@ -85,7 +85,10 @@ function closure(x::AbstractMatrix{<:Real}) all(v -> isfinite(v) && v > 0, x) || throw(CompositionError("closure needs strictly positive finite parts")) # Scale before summing: both integer totals and finite floating-point # totals can overflow even though every closed proportion is representable. - scaled = Float64.(x ./ maximum(x; dims = 2)) + # Widen narrow floats before division, while leaving large integers and + # arbitrary-precision values in their original scaling domain. + wide = eltype(x) <: Union{Float16,Float32} ? Float64.(x) : x + scaled = Float64.(wide ./ maximum(wide; dims = 2)) out = scaled ./ sum(scaled; dims = 2) all(v -> isfinite(v) && v > 0, out) || throw(CompositionError("closed proportions are not representable as positive Float64 values")) diff --git a/src/Core/zeros.jl b/src/Core/zeros.jl index 182174f..adde511 100644 --- a/src/Core/zeros.jl +++ b/src/Core/zeros.jl @@ -88,6 +88,12 @@ function replace_zeros(p::MultiplicativeReplacement, counts::AbstractMatrix{<:In r = nz == 0 ? 0.0 : Float64(p.δ / scale / scaled_total) nz == 0 || r > 0 || throw(CompositionError("sample $i: zero replacement underflows Float64")) shrink = 1 - nz * r + if abs(shrink) <= 8eps(Float64) + # Near cancellation, use the exact count total and stored δ so + # rounding r cannot erase a positive residual. BigInt avoids overflow. + total = sum(BigInt, row) + shrink = Float64((total - nz * Rational{BigInt}(p.δ)) / total) + end shrink > 0 || throw(CompositionError("sample $i: replacement δ too large for $nz zeros")) for j in axes(counts, 2) out[i, j] = iszero(row[j]) ? r : shrink * scaled[j] / scaled_total diff --git a/test/fuzz_core.jl b/test/fuzz_core.jl index 69b31c2..18b65d2 100644 --- a/test/fuzz_core.jl +++ b/test/fuzz_core.jl @@ -1,10 +1,31 @@ # SPDX-License-Identifier: MPL-2.0 # SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# A policy that could hide invalid counts if the CLR boundary skipped validation. +struct AbsoluteCountsPolicy <: ZeroPolicy end +CompositionalDA.replace_zeros(::AbsoluteCountsPolicy, counts::AbstractMatrix{<:Integer}) = + Float64.(abs.(counts)) .+ 1 + @testset "Adversarial core boundaries" begin @test closure(fill(typemax(Int), 1, 2)) ≈ fill(0.5, 1, 2) @test closure(fill(floatmax(Float64), 1, 2)) ≈ fill(0.5, 1, 2) @test closure(fill(big(10)^1000, 1, 2)) ≈ fill(0.5, 1, 2) + for T in (Float16, Float32) + narrow = [nextfloat(zero(T)) floatmax(T)] + expected = Float64.(BigFloat.(narrow) ./ sum(BigFloat.(narrow))) + @test all(>(0), closure(narrow)) + @test closure(narrow) ≈ expected rtol=2e-15 + end + @test_throws CompositionError clr(AbsoluteCountsPolicy(), [-1 2]) + @test clr(AbsoluteCountsPolicy(), [1 2]) ≈ clr([2.0 3.0]) + near_one = replace_zeros(MultiplicativeReplacement(prevfloat(1.0)), [3 2 0 0 0 0 0]) + @test all(>(0), near_one) + @test sum(near_one) ≈ 1.0 + @test near_one[1, 1] / near_one[1, 2] ≈ 3 / 2 + @test near_one[1, 1] ≈ 3 / 5 * (1 - prevfloat(1.0)) + many_zeros = reshape(vcat([1025, 1025], zeros(Int, 2050)), 1, :) + @test all(>(0), replace_zeros(MultiplicativeReplacement(prevfloat(1.0)), many_zeros)) + @test_throws CompositionError replace_zeros(MultiplicativeReplacement(prevfloat(1.0)), [1 0 0]) for empty in (zeros(0, 2), zeros(1, 0)) @test_throws CompositionError closure(empty) @test_throws CompositionError clr(empty) From 911989080145b9ab43ce06a8c1c8dfe6712ca2e3 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 20:08:31 +0000 Subject: [PATCH 5/9] Preserve BH pointwise bounds and widen abstractly typed float inputs A predeclared calibration pilot exposed multiply/divide cancellation that could round a BH result below its input. Form the scale first and test continuous random probabilities. Independently audit the narrow-float repair on Matrix{Real}, preserving arbitrary-precision inputs. Complete Windows Julia suite: 20,341 assertions passed. --- src/Core/bh.jl | 4 +++- src/Core/composition.jl | 2 +- test/fuzz_core.jl | 12 ++++++++++++ 3 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/Core/bh.jl b/src/Core/bh.jl index 9fe0529..9e51554 100644 --- a/src/Core/bh.jl +++ b/src/Core/bh.jl @@ -23,7 +23,9 @@ function bh_adjust(p::AbstractVector{<:Real}) adj = similar(q) for (r, k) in enumerate(ord) rank = m - r + 1 - running = min(running, q[k] * m / rank) + # Form the scale first: (p*m)/m can round below p even at rank m. + # m/rank >= 1, so this multiplication preserves q >= p in Float64. + running = min(running, (m / rank) * q[k]) adj[k] = running end out[idx] = adj diff --git a/src/Core/composition.jl b/src/Core/composition.jl index e3f7569..6762886 100644 --- a/src/Core/composition.jl +++ b/src/Core/composition.jl @@ -87,7 +87,7 @@ function closure(x::AbstractMatrix{<:Real}) # totals can overflow even though every closed proportion is representable. # Widen narrow floats before division, while leaving large integers and # arbitrary-precision values in their original scaling domain. - wide = eltype(x) <: Union{Float16,Float32} ? Float64.(x) : x + wide = map(v -> v isa Union{Float16,Float32} ? Float64(v) : v, x) scaled = Float64.(wide ./ maximum(wide; dims = 2)) out = scaled ./ sum(scaled; dims = 2) all(v -> isfinite(v) && v > 0, out) || diff --git a/test/fuzz_core.jl b/test/fuzz_core.jl index 18b65d2..68eee7c 100644 --- a/test/fuzz_core.jl +++ b/test/fuzz_core.jl @@ -15,6 +15,7 @@ CompositionalDA.replace_zeros(::AbsoluteCountsPolicy, counts::AbstractMatrix{<:I expected = Float64.(BigFloat.(narrow) ./ sum(BigFloat.(narrow))) @test all(>(0), closure(narrow)) @test closure(narrow) ≈ expected rtol=2e-15 + @test closure(Matrix{Real}(narrow)) ≈ expected rtol=2e-15 end @test_throws CompositionError clr(AbsoluteCountsPolicy(), [-1 2]) @test clr(AbsoluteCountsPolicy(), [1 2]) ≈ clr([2.0 3.0]) @@ -58,6 +59,17 @@ CompositionalDA.replace_zeros(::AbsoluteCountsPolicy, counts::AbstractMatrix{<:I @test_throws CompositionError CompositionalDA.validate_two_groups(x, ["a", nothing], [1, 1, 2, 2]) end +@testset "BH preserves pointwise bounds without multiply-divide cancellation" begin + @test all(bh_adjust([0.1,0.3,0.7]) .>= [0.1,0.3,0.7]) + for seed in (23,20261107,0xcafe) + rng=StableRNG(seed) + for _ in 1:256 + probabilities=rand(rng,rand(rng,1:100)) + @test all(probabilities .<= bh_adjust(probabilities) .<= 1) + end + end +end + @testset "Seeded count and composition fuzz" begin for seed in (23, 20261107, 0xcafe), case in 1:128 @testset "seed=$seed case=$case" begin From cbb8ac4d85b4259c8122ed83c2cce8ec6c34362a Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 20:30:58 +0000 Subject: [PATCH 6/9] Record launch blockers and add reproducible calibration and benchmark experiments Keep registration and scientific approval pending. Document the predeclared adverse small-unbalanced ANCOM-BC null result, the distinction between algebraic proof and calibration, and actual fleet/standards gaps. Add a pinned repaired-reference investigation with convergence receipts, separate benchmark environment and machine-specific measurement harness. --- .github/workflows/ci.yml | 13 +++ EXPLAINME.adoc | 62 ++++++++++++++ README.md | 9 +++ SECURITY.md | 25 ++++++ benchmark/Project.toml | 11 +++ benchmark/README.md | 27 +++++++ benchmark/benchmarks.jl | 41 ++++++++++ benchmark/calibration_pilot.jl | 116 +++++++++++++++++++++++++++ benchmark/calibration_reference.R | 44 ++++++++++ benchmark/export_calibration_case.jl | 26 ++++++ docs/CALIBRATION-PILOT-PROTOCOL.adoc | 64 +++++++++++++++ docs/LAUNCH-READINESS.adoc | 92 +++++++++++++++++++++ 12 files changed, 530 insertions(+) create mode 100644 EXPLAINME.adoc create mode 100644 SECURITY.md create mode 100644 benchmark/Project.toml create mode 100644 benchmark/README.md create mode 100644 benchmark/benchmarks.jl create mode 100644 benchmark/calibration_pilot.jl create mode 100644 benchmark/calibration_reference.R create mode 100644 benchmark/export_calibration_case.jl create mode 100644 docs/CALIBRATION-PILOT-PROTOCOL.adoc create mode 100644 docs/LAUNCH-READINESS.adoc diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 61e18d7..7f014ca 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -143,3 +143,16 @@ jobs: curl --fail --location --retry 3 https://raw.githubusercontent.com/hyperpolymath/ANCOMBC/6d094731ed85e9b0458053ad3510e9bdc238e152/R/ancombc_bias_correct.R -o /tmp/ancombc-variance-repair.R export R_HOME="$(R RHOME)" ANCOMBC_ESTIMATOR_SOURCE=/tmp/ancombc-variance-repair.R julia --project=test/oracle test/oracle/ancombc_oracle.jl + - name: Investigate the small-unbalanced calibration failure + if: matrix.reference == 'ANCOMBC' + run: | + julia --project=test/oracle benchmark/export_calibration_case.jl /tmp/calibration/small-unbalanced-inputs + Rscript benchmark/calibration_reference.R /tmp/calibration/small-unbalanced-inputs /tmp/ancombc-variance-repair.R + - name: Preserve calibration investigation evidence + if: always() && matrix.reference == 'ANCOMBC' + uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # v4.6.2 + with: + name: small-unbalanced-calibration + path: /tmp/calibration + if-no-files-found: error + retention-days: 30 diff --git a/EXPLAINME.adoc b/EXPLAINME.adoc new file mode 100644 index 0000000..b9fa97e --- /dev/null +++ b/EXPLAINME.adoc @@ -0,0 +1,62 @@ +// SPDX-License-Identifier: CC-BY-SA-4.0 +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) += What CompositionalDA can establish + +This package analyses raw integer counts in samples-by-features orientation. +Its current methods are unpaired two-group ALDEx2 and original ANCOM-BC. +It does not implement ANCOM-BC2, covariate adjustment or repeated measures. +Zero handling, the contrast direction, seeds and numerical policies belong in +the analysis record. Swapping groups changes the effect direction. + +== Evidence and its limits + +[cols="2,3,3",options="header"] +|=== +|Claim |Owning evidence |What remains outside it +|Exact algebraic identities |`proofs/PROOF-STATUS.md`, safe Agda modules and negative controls +|Floating-point refinement, sorter correctness, FDR and biological truth +|Numerical method agreement |`test/oracle/`, pinned R/Bioconductor releases and explicitly repaired reference functions +|Equivalence to pristine defective references, every input or later upstream releases +|Regression resistance |`test/`, seeded fuzz/property tests and planted mutants +|Exhaustive path/state coverage or the absence of unknown defects +|Real processed microbiome inputs |`docs/DATASETS.md`, archived source hashes and gut/soil fixtures +|Known ground truth, prospective clinical validity or general calibration +|Exploratory operating characteristics |`docs/CALIBRATION-PILOT-PROTOCOL.adoc`, `benchmark/calibration_pilot.jl` +|A calibration certificate; the small-unbalanced ANCOM-BC null pilot is adverse +|Performance |`benchmark/benchmarks.jl`, machine-specific raw measurements +|A speed advantage over R or an agreed latency/service-level guarantee +|=== + +The 2026-11-07 calibration pilot fixed six regimes and 32 independently seeded +datasets per regime before measurement. Of 384 method fits, 365 returned and +19 refused. Rates conditional on returned fits must be accompanied by refusal +counts. Under the small-unbalanced complete null, original ANCOM-BC made false +discoveries in 19/29 successful fits at 0.05 (Wilson 95% interval 0.473–0.801). +The unconditional observed count is 19/32; refusals are not successful negative +findings. This requires investigation against the repaired R reference and a +reviewed support restriction or estimator correction before a calibration claim. +Do not fix it by tuning the pilot seeds, removing adverse regimes or quietly +substituting another method. + +ALDEx2's pilot is much less adverse, but 32 datasets cannot certify 5% FDR. +All-null FDR is the probability of any rejection. With signal present, that +probability is family-wise error; mean false discovery proportion estimates +FDR instead. The reported Wilson intervals are for event probabilities, not +mean FDP, taxa treated as independent replications, or arbitrary populations. + +At only two ALDEx2 Monte-Carlo instances, a separate 4,096-seed probe produced +58 infinite and three undefined effect estimates. Its p/q bounds still held. +The WebUI refuses non-finite effects; this is an estimator usability limitation +to address explicitly, not permission to coerce undefined effects to zero. + +== Trust and release decisions + +Julia, GMP/MPFR, BLAS, dependencies, the operating system, reference software, +data provenance and proof checkers remain in the trusted environment. +Agda algebra does not independently verify those implementations. No second +Rocq/Lean witness or full executable refinement is claimed. + +General registration is explicitly on hold pending the human maintainer's review +of the changed implementation. A signature authenticates a commit, not its +mathematical or scientific correctness. Version 0.3.0 in Project.toml is not a +stable-release grade. No CRG/TRG/RSR certificate is asserted. diff --git a/README.md b/README.md index bc115d1..8522287 100644 --- a/README.md +++ b/README.md @@ -17,6 +17,15 @@ backend. > gate now checks all statistics and shared Monte-Carlo instances. ANCOM-BC's > precise model and limitations are described in `docs/ANCOMBC-SPEC.md`. +> **Launch review remains pending.** The predeclared calibration pilot found +> severe false-discovery inflation for original ANCOM-BC with four versus +> sixteen samples under a complete null: 19 of 29 successful fits made at least +> one discovery at q <= 0.05; three additional fits refused to return a result. +> This is a release blocker for a general calibration claim. See +> [`docs/LAUNCH-READINESS.adoc`](docs/LAUNCH-READINESS.adoc), the reproducible +> protocol and scripts in `benchmark/`, and [`EXPLAINME.adoc`](EXPLAINME.adoc). +> Passing numerical reference comparisons does not establish FDR control. + The v0.3 implementation, tests, proofs and integration work were developed with OpenAI Codex assistance. General registration requires a human maintainer to review and understand these contributions; automated checks and code review diff --git a/SECURITY.md b/SECURITY.md new file mode 100644 index 0000000..4933997 --- /dev/null +++ b/SECURITY.md @@ -0,0 +1,25 @@ + + + +# Security and numerical integrity + +This is a development package; no released version currently carries a +security support commitment. Report exploitable defects through GitHub's +private vulnerability-reporting facility if enabled. Otherwise contact the +maintainer at the address in Project.toml before publishing exploit details. +Ordinary numerical defects can be reported as issues with a minimal synthetic +reproducer, package/tool versions and expected versus observed results. +Do not attach confidential or identifiable biological data. + +The native API is a local scientific library, not a sandbox for untrusted jobs. +Callers must bound matrix dimensions, count bit widths, Monte-Carlo work, +concurrency and memory at their service boundary. Input validation and explicit +nonconvergence errors do not create a hard CPU or memory limit. R reference +execution loads installed R packages and must run in a separately controlled +environment when inputs or dependencies are untrusted. + +Aqua checks Julia package quality. It is not a vulnerability-database scanner. +The launch ledger must separately record dependency advisories, secret scanning, +workflow permissions/action pins, release provenance and unresolved findings. +No successful comprehensive security audit is claimed here. Statistics that +fail calibration are safety/integrity concerns even without an attacker. diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000..edd428a --- /dev/null +++ b/benchmark/Project.toml @@ -0,0 +1,11 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CompositionalDA = "bb1787b7-c25e-80cd-8559-27e3f8d35816" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" + +[compat] +BenchmarkTools = "1" +Distributions = "0.25" +julia = "1.10" diff --git a/benchmark/README.md b/benchmark/README.md new file mode 100644 index 0000000..535f89e --- /dev/null +++ b/benchmark/README.md @@ -0,0 +1,27 @@ + + + +# Exploratory launch measurements + +From the repository root, prepare the separate environment: + +```sh +julia --project=benchmark -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()' +julia --project=benchmark benchmark/calibration_pilot.jl /path/to/pilot-output +julia --project=benchmark benchmark/benchmarks.jl /path/to/benchmark-output +``` + +Read `docs/CALIBRATION-PILOT-PROTOCOL.adoc` before interpreting the output. +Seeds, scenarios and comparisons were fixed before results. Refusals are +separate from successful fits. Neither script certifies calibration, speed, +stability or package readiness. Keep generated measurements outside source. +Benchmark timings exclude compilation/warmup and retain trial counts; ten +trials provide an exploratory baseline, not precise tail-latency estimates. + +The optional CI investigation exports the small-unbalanced complete-null +inputs and native results, then runs pinned ANCOMBC 2.12.0 with the explicitly +identified variance-estimator repair. It records R convergence rather than +treating an iteration-cap result as a valid fit. A successful investigation +job means the experiment completed, not that the adverse calibration passed. +The raw artifact expires after 30 days; archive it with its SHA-256 before +using it in a release decision. diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 0000000..8d2f243 --- /dev/null +++ b/benchmark/benchmarks.jl @@ -0,0 +1,41 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +using BenchmarkTools, CompositionalDA, Statistics, TOML, Dates +include("calibration_pilot.jl") +destination=only(ARGS) +mkpath(destination) +rows=NamedTuple[] +for (n,d) in ((20,40),(20,400),(100,400)) + regime=(name="benchmark",n1=n÷2,n2=n÷2,depth1=10000,depth2=10000,sd=.35,sparse=false,signals=4) + x,f,g,_=simulate(regime,20261107+n+d;d) + for method in (:closure,:clr,:bh,:aldex2,:ancombc) + p=collect(range(.001,.999;length=d)) + call = method==:closure ? (() -> closure(x)) : + method==:clr ? (() -> clr(Pseudocount(1),x)) : + method==:bh ? (() -> bh_adjust(p)) : + method==:aldex2 ? (() -> aldex2(x,f,g;mc_samples=128,seed=123)) : + (() -> ancombc(x,f,g;structural_zeros=false)) + try + call() # compilation/warmup is excluded and labelled as such + trial=@benchmark $call() samples=10 evals=1 seconds=15 + push!(rows,(method=String(method),samples=n,features=d,status="measured", + trials=length(trial.times),median_ns=median(trial.times),p95_ns=quantile(trial.times,.95), + bytes=trial.memory,allocations=trial.allocs)) + catch error + error isa CompositionError || rethrow() + push!(rows,(method=String(method),samples=n,features=d,status="refused", + trials=0,median_ns=NaN,p95_ns=NaN,bytes=0,allocations=0)) + println("REFUSED ",method," ",n,"x",d,": ",sprint(showerror,error)) + end + end +end +open(joinpath(destination,"benchmarks.csv"),"w") do io + println(io,join(keys(first(rows)),',')) + foreach(row->println(io,join(values(row),',')),rows) +end +metadata=Dict("julia"=>string(VERSION),"threads"=>Threads.nthreads(),"cpu"=>Sys.CPU_NAME, + "machine"=>Sys.MACHINE,"word_size"=>Sys.WORD_SIZE,"date_utc"=>string(now(UTC)), + "source_head"=>strip(read(`git -C $(@__DIR__) rev-parse HEAD`,String)), + "benchmarktools"=>string(pkgversion(BenchmarkTools)),"warmup_excluded"=>true, + "claim"=>"exploratory warm-process baseline; no reference speedup or acceptance budget") +open(io->TOML.print(io,metadata),joinpath(destination,"benchmarks-metadata.toml"),"w") diff --git a/benchmark/calibration_pilot.jl b/benchmark/calibration_pilot.jl new file mode 100644 index 0000000..81415cc --- /dev/null +++ b/benchmark/calibration_pilot.jl @@ -0,0 +1,116 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# Run with the package project; results are written only to the supplied directory. +using CompositionalDA, Random, Statistics, Dates, Printf, SHA, TOML +using CompositionalDA: StableRNG +using Distributions: Multinomial + +const REGIMES = [ + (name="balanced_null",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=0), + (name="depth_imbalance_null",n1=10,n2=10,depth1=1000,depth2=10000,sd=0.35,sparse=false,signals=0), + (name="small_unbalanced_null",n1=4,n2=16,depth1=1000,depth2=1000,sd=0.35,sparse=false,signals=0), + (name="sparse_null",n1=6,n2=6,depth1=200,depth2=200,sd=0.7,sparse=true,signals=0), + (name="sparse_signal",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=4), + (name="dense_balanced_signal",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=20), +] + +function simulate(regime,seed;d=40) + rng=StableRNG(seed) + baseline=regime.sparse ? exp.(range(-4,0;length=d)) : collect(range(1000,5000;length=d)) + logeffect=zeros(d) + half=regime.signals÷2 + logeffect[1:half].=log(4) + logeffect[half+1:2half].=-log(4) + x=Matrix{Int}(undef,regime.n1+regime.n2,d) + for i in axes(x,1) + condition=i>regime.n1 + weights=baseline .* exp.(regime.sd .* randn(rng,d) .+ condition .* logeffect) + x[i,:]=rand(rng,Multinomial(condition ? regime.depth2 : regime.depth1,weights/sum(weights))) + end + return x,string.(1:d),vcat(fill(0,regime.n1),fill(1,regime.n2)),logeffect +end + +function wilson(events,total) + total==0 && return (NaN,NaN) + z=1.959963984540054; p=events/total; scale=1+z^2/total + centre=(p+z^2/(2total))/scale + radius=z*sqrt(p*(1-p)/total+z^2/(4total^2))/scale + return (max(0,centre-radius),min(1,centre+radius)) +end + +function run_pilot(destination) + mkpath(destination) + protocol=joinpath(@__DIR__,"..","docs","CALIBRATION-PILOT-PROTOCOL.adoc") + metadata=Dict("protocol_sha256"=>bytes2hex(sha256(read(protocol))), + "julia"=>string(VERSION),"package_version"=>string(pkgversion(CompositionalDA)), + "source_head"=>strip(read(`git -C $(@__DIR__) rev-parse HEAD`,String)), + "started_utc"=>string(now(UTC)),"replicates_per_regime"=>32, + "mc_samples"=>128,"seed_formula"=>"20261107 + 10000*regime_index + replicate", + "status"=>"exploratory pilot; not a calibration certificate") + records=NamedTuple[] + failures=NamedTuple[] + for (regime_index,regime) in enumerate(REGIMES), replicate in 1:32 + seed=20261107+10000regime_index+replicate + x,features,groups,truth=simulate(regime,seed) + for method in (:aldex2,:ancombc) + start=time_ns() + fitted=try + method==:aldex2 ? aldex2(x,features,groups;mc_samples=128,seed=seed+1_000_000) : + ancombc(x,features,groups;prevalence=0.1,structural_zeros=false,conserve=true) + catch error + if error isa CompositionError + push!(failures,(regime=regime.name,replicate,method=string(method),seed, + reason=sprint(showerror,error),seconds=(time_ns()-start)/1e9)) + continue + end + rethrow() + end + elapsed=(time_ns()-start)/1e9 + table=fitted.table + @assert all(0 .<= table.pvalue .<= table.qvalue .<= 1) + ids=parse.(Int,table.feature) + for alpha in (0.05,0.10) + discovered=table.qvalue .<= alpha + false_calls=count(discovered .& iszero.(truth[ids])) + true_calls=count(discovered .& .!iszero.(truth[ids])) + push!(records,(regime=regime.name,replicate,method=string(method),seed,alpha, + tested=length(ids),excluded=40-length(ids),false_calls,true_calls, + fdp=false_calls/max(false_calls+true_calls,1), + power=regime.signals==0 ? NaN : true_calls/regime.signals,seconds=elapsed)) + end + end + replicate==32 && println("Completed ",regime.name) + flush(stdout) + end + open(joinpath(destination,"calibration-pilot.csv"),"w") do io + println(io,join(string.(keys(first(records))),',')) + for record in records + println(io,join(values(record),',')) + end + end + metadata["finished_utc"]=string(now(UTC)) + metadata["successful_fits"]=length(records)÷2 + metadata["refused_fits"]=length(failures) + metadata["refusals"]=[Dict(string(k)=>v for (k,v) in pairs(f)) for f in failures] + open(io->TOML.print(io,metadata),joinpath(destination,"calibration-pilot-metadata.toml"),"w") + open(joinpath(destination,"calibration-pilot-summary.md"),"w") do io + println(io,"# Exploratory calibration pilot\n\n32 independent datasets per regime. Metrics condition on successful fits; refusals are separate. No calibration certification.\n") + println(io,"| Regime | Method | q threshold | Fits/32 | Mean FDP | Mean power | Any false discovery: fraction [95% Wilson interval] |") + println(io,"| --- | --- | --- | --- | --- | --- | --- |") + for regime in REGIMES,method in ("aldex2","ancombc"),alpha in (0.05,0.10) + rows=filter(r->r.regime==regime.name && r.method==method && r.alpha==alpha,records) + isempty(rows) && (println(io,"| $(regime.name) | $method | $alpha | 0/32 | — | — | — |");continue) + events=count(r->r.false_calls>0,rows) + lo,hi=wilson(events,length(rows)) + power=regime.signals==0 ? "N/A" : @sprintf("%.3f",mean(r.power for r in rows)) + @printf(io,"| %s | %s | %.2f | %d/32 | %.3f | %s | %.3f [%.3f, %.3f] |\n", + regime.name,method,alpha,length(rows),mean(r.fdp for r in rows),power,events/length(rows),lo,hi) + end + println(io,"\nFor complete-null regimes, any false discovery equals the false-discovery-rate event. For signal regimes that event is FWER, not mean FDP. The interval covers the event probability among successful fits, not arbitrary populations or the mean FDP. See the fixed protocol and refusal records.") + end + println("Successful fits: ",metadata["successful_fits"],"; refusals: ",length(failures)) +end + +if abspath(PROGRAM_FILE)==(@__FILE__) + run_pilot(isempty(ARGS) ? error("Supply an output directory") : only(ARGS)) +end diff --git a/benchmark/calibration_reference.R b/benchmark/calibration_reference.R new file mode 100644 index 0000000..ae2ba82 --- /dev/null +++ b/benchmark/calibration_reference.R @@ -0,0 +1,44 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +args <- commandArgs(trailingOnly=TRUE) +input <- args[[1]] +loadNamespace("ANCOMBC") +stopifnot(as.character(packageVersion("ANCOMBC")) == "2.12.0") +replacement <- new.env(parent=asNamespace("ANCOMBC")) +sys.source(args[[2]], envir=replacement) +assignInNamespace(".bias_em",replacement$.bias_em, ns="ANCOMBC") +results <- list() +for (replicate in seq_len(32)) { + reads <- as.matrix(read.delim(file.path(input,paste0("replicate-",replicate,".tsv")),row.names=1,check.names=FALSE)) + meta <- data.frame(condition=factor(c(rep("reference",4),rep("contrast",16)),levels=c("reference","contrast")), + unused=seq_len(20),row.names=colnames(reads)) + diagnostics <- list() + trace(".bias_em",where=asNamespace("ANCOMBC"),print=FALSE,exit=quote({ + .GlobalEnv$diagnostics[[length(.GlobalEnv$diagnostics)+1L]] <- c(iterations=iterNum,residual=epsilon) + })) + fit <- tryCatch(ANCOMBC::ancombc(data=reads,aggregate_data=reads,meta_data=meta, + tax_level=NULL,formula="condition",p_adj_method="BH",pseudo=1,prv_cut=.1,lib_cut=0, + group="condition",struc_zero=FALSE,neg_lb=FALSE,tol=1e-7,max_iter=20000, + conserve=TRUE,global=FALSE,n_cl=1,verbose=FALSE),error=function(e)e) + untrace(".bias_em",where=asNamespace("ANCOMBC")) + if(inherits(fit,"error")) { + results[[replicate]] <- data.frame(replicate,status=conditionMessage(fit),false_calls=NA,iterations=NA,residual=NA,native_false_calls=NA,max_q_gap=NA) + } else { + diag <- do.call(rbind,diagnostics) + converged <- nrow(diag)>0 && all(diag[,"iterations"]<20000) && all(diag[,"residual"]<=1e-7) + q <- fit$res$q_val$conditioncontrast + native_file <- file.path(input,paste0("native-",replicate,".csv")) + native_calls <- max_q_gap <- NA_real_ + if (file.exists(native_file)) { + native <- read.csv(native_file,colClasses=c(feature="character")) + ids <- if("taxon" %in% names(fit$res$q_val)) as.character(fit$res$q_val$taxon) else rownames(fit$res$q_val) + stopifnot(setequal(ids,native$feature),all(is.finite(q))) + native_calls <- sum(native$qvalue<=.05) + max_q_gap <- max(abs(native$qvalue-q[match(native$feature,ids)])) + } + results[[replicate]] <- data.frame(replicate,status=if(converged) "converged" else "unconverged", + false_calls=sum(q<=.05,na.rm=TRUE),iterations=max(diag[,"iterations"]),residual=max(diag[,"residual"]),native_false_calls=native_calls,max_q_gap) + } + write.csv(do.call(rbind,results),file.path(dirname(input),"small-unbalanced-repaired-R.csv"),row.names=FALSE) + cat("Finished replicate",replicate,"\n"); flush.console() +} diff --git a/benchmark/export_calibration_case.jl b/benchmark/export_calibration_case.jl new file mode 100644 index 0000000..14b87c3 --- /dev/null +++ b/benchmark/export_calibration_case.jl @@ -0,0 +1,26 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +include("calibration_pilot.jl") +destination=only(ARGS) +mkpath(destination) +for replicate in 1:32 + x,features,groups,truth=simulate(REGIMES[3],20261107+30000+replicate) + open(joinpath(destination,"replicate-$replicate.tsv"),"w") do io + println(io,join(vcat("feature",["s$i" for i in axes(x,1)]),'\t')) + for j in axes(x,2) + println(io,join(vcat(features[j],string.(x[:,j])),'\t')) + end + end + try + fitted=ancombc(x,features,groups;prevalence=.1,structural_zeros=false,conserve=true) + open(joinpath(destination,"native-$replicate.csv"),"w") do io + println(io,"feature,pvalue,qvalue,se,effect") + for row in eachrow(fitted.table) + println(io,join((row.feature,row.pvalue,row.qvalue,row.se,row.effect),',')) + end + end + catch error + error isa CompositionError || rethrow() + write(joinpath(destination,"native-refused-$replicate.txt"),sprint(showerror,error)) + end +end diff --git a/docs/CALIBRATION-PILOT-PROTOCOL.adoc b/docs/CALIBRATION-PILOT-PROTOCOL.adoc new file mode 100644 index 0000000..2c54380 --- /dev/null +++ b/docs/CALIBRATION-PILOT-PROTOCOL.adoc @@ -0,0 +1,64 @@ +// SPDX-License-Identifier: CC-BY-SA-4.0 +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) += Calibration pilot protocol + +This protocol is fixed before examining pilot results. It is exploratory +evidence for planning a larger independent study, not a calibration certificate +or a release acceptance test. No tolerance is tuned to observed results. + +== Estimands and regimes + +Generate independent samples with 40 absolute-abundance components, lognormal +biological variation and multinomial sampling at a specified library depth. +Use baseline abundances spanning 1,000–5,000, with log-scale SD 0.35 unless +specified. Conditions are labelled 0 and 1. Signal scenarios apply equal +numbers of fourfold increases and fourfold decreases, so the sum of log fold +changes is zero. In that restricted construction, population CLR differences +and absolute log-abundance differences have the same zero/nonzero support. +This alignment is intentional; it does not make these estimands equivalent +on arbitrary data. + +Run 32 independent datasets per regime, seeds determined before execution: + +* Balanced null: 10/10 samples, depth 10,000, no group effect. +* Depth imbalance null: 10/10, depths 1,000/10,000, no abundance effect. +* Small unbalanced null: 4/16, depth 1,000, no abundance effect. +* Sparse null: 6/6, depth 200, baseline abundances geometrically spanning + exp(-4)–1, biological log SD 0.7, no abundance effect. +* Sparse signal: 10/10, depth 10,000, two features up fourfold and two down + fourfold; 36 null features. +* Dense balanced signal: 10/10, depth 10,000, ten features up and ten down + fourfold; 20 null features. Treat this as sensitivity to the bias estimator's + predominantly-null identification premise. + +ALDEx2 uses 128 Monte Carlo instances and a separately recorded seed. ANCOM-BC +uses its documented conservative setting, prevalence 0.1, structural-zero +decisions disabled, tolerance 1e-7 and maximum 1,000 iterations. Refusals are +recorded separately and are never interpreted as successful zero discoveries. +Analyse reported q-values at 0.05 and 0.10. Record true/false discoveries, +false-discovery proportion, power, tested/excluded features and runtime. + +== Uncertainty and interpretation + +Report the fit-success fraction and metrics conditional on successful fits. +For each null regime, report the probability of at least one false discovery +(equal to FDR under the complete null) with a 95% Wilson interval over +independent datasets, not over correlated features. For signal scenarios +report mean false-discovery proportion and power, with transparent sample +counts; avoid a narrow-confidence claim from this pilot. + +Thirty-two replicates cannot certify a 5% false-positive bound: even zero +events leaves a substantial upper confidence limit. The study does not cover +all distributions, batch effects, arbitrary zero mechanisms, repeated +measures, ecological identification, external deployment or user acceptance. +The six scenarios are not six unrelated external CRG validation targets. +Reference-program agreement is evaluated separately on the existing real +gut/soil fixtures and mock data; it is not substituted for calibration here. + +== Provenance and reuse + +The estate's proven-tests-and-benches framework was inspected first. It +provides evidence/coverage abstractions, not a ready-made Julia compositional +calibration experiment for this estimand. This domain-specific pilot therefore +records its own assumptions, seeds, results and limitations. Its synthetic +data contain no human records. Preserve the protocol hash alongside results. diff --git a/docs/LAUNCH-READINESS.adoc b/docs/LAUNCH-READINESS.adoc new file mode 100644 index 0000000..dafb6cc --- /dev/null +++ b/docs/LAUNCH-READINESS.adoc @@ -0,0 +1,92 @@ +// SPDX-License-Identifier: CC-BY-SA-4.0 +// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) += Launch readiness: evidence before labels + +Assessment date: 2026-11-07. Candidate numerical source: +`911989080145b9ab43ce06a8c1c8dfe6712ca2e3`. +This is the author's assessment, not independent scientific review. +Registration remains on hold. No alpha/beta/stable or RSR certification is granted. + +== Source authority + +Estate standards were inspected at `fe57d352f0b1a50e5cb07b4aac2581caf355388e`, +gitbot-fleet at `46d26a4b084c40847d9a2f2ef98305502b7320bc`, and +proven-tests-and-benches at `e82c54ea6ef0e004f759f3533ebe8373e517f9e3`. +Refresh these before release. Joshua's CONTRIBUTING governs the WebUI contribution. +Current estate metadata canon uses one chora DEED, not retired A2ML/SCM files. +Conflicting old examples must not be copied as if authoritative. + +CRG assesses component readiness; TRG imposes additional language-toolchain +requirements. This scientific library is not a compiler/runtime merely because +it is written in Julia. SemVer, General eligibility and readiness grades are +separate. Six simulation regimes do not meet CRG's diverse external-target bar. +The taxonomy is 18 test categories by 14 aspects; the companion evidence +framework still documents a 17-by-14 schema. Preserve that migration gap. + +== Gates and work allocation + +[cols="2,3,3",options="header"] +|Gate or owner |Required evidence |Current disposition +|Finishbot (finishingbot) |README/EXPLAINME, install example, licences, About/topics, links, changelog and release description all agree with source +|About/topics corrected; explanation/security/launch docs proposed. Wiki, changelog, native estate metadata and final release packet remain pending. +|Seambot |Counts/orientation/groups, zero policy, Julia/R, package/WebUI, JSON/CSV, registry/install seams each have owners and positive/negative contract tests +|Substantial adapter and export tests exist. Current fleet parser list has no Julia parser; do not report automated Julia seam coverage without an implemented bridge. +|Echidnabot |Claim-to-theorem map, checker versions, safe/axiom audit, negative controls, counterexamples and executable refinement limits +|Safe Agda and negative controls pass. Independent proof witnesses and floating-point refinement remain pending; backend inventory is not a successful proof run. +|Rhodibot / Hypatia |Canonical metadata checks, repository policy and machine-readable findings with source hashes +|Estate conformance not established; retired metadata examples are an unresolved shared-guidance problem. +|Cipherbot |Actual dependency/secret/security scans, pinned CI provenance and triaged findings +|No complete scanner receipt yet; Aqua is quality evidence only. +|Glambot / Accessibilitybot |Representative branding, correct alt text, legibility at small sizes, text alternatives and no invented assurance marks +|Social concept approved by user; production export dimensions/size and avatar/hero review remain pending. +|Sustainabot |Measured runtime, allocations, scaling, work budgets and reproducible hardware/runtime context +|Benchmark harness proposed; no environmental-impact claim or performance advantage established. +|RobotRepoAutomaton |Execute authorised repairs with regression evidence and preserve provenance +|May implement fixes; cannot approve its own scientific or human-review gate. +|Panicbot |Escalate severe integrity/security failures without suppressing evidence +|Adverse small-unbalanced ANCOM-BC calibration must block a general calibration claim. +|Human maintainer and independent scientist |Review actual diff by hand, assess applicability, data rights, uncertainty and accepted claims +|Pending; automated or neural reviews do not count as those decisions. +|=== + +These are proposed objectives for the real fleet roles, not claims that a bot +ran or supports every language/backend. The fleet's finishing documentation +still mentions retired six-file SCM metadata; it needs reconciliation with +current canon before that objective can be marked complete. + +== Stage boundaries + +* Research/pre-alpha: explicit scope/refusals, reproducible baseline, no unsupported + calibration/clinical/stable claim. Adverse findings stay visible. +* Alpha candidate: applicable CRG minimum and canonical repository requirements, + tested installation and user workflow, signed provenance, known-bug treatment, + and human approval of the exact bounded claim. Test count alone is insufficient. +* Beta candidate: independent reproduction and scientific review; adequately + powered, prespecified calibration; relevant external datasets/targets; mutation, + fuzz and seam campaigns; measured performance/support budgets and regression policy. +* Stable candidate: satisfy the applicable estate V1 gate, including all claimed + proof/test/benchmark/security obligations and artifact/claim parity. A green CI + icon or version number cannot substitute for these requirements. + +General submission additionally needs correct name/UUID/compat/licence, clean +installation without an unregistered dependency, declared AI assistance and the +human review General requires. Registration review, the waiting period, tagging +and release publication remain distinct events. Other registries need their own +adapters for namespace, licence, provenance, binaries, signing and review rules. + +== Risk and uncertainty + +Prioritise severity, exposure, detectability and confidence; retain the raw +evidence rather than treating an ordinal risk score as a probability. Calibration +failure outranks cosmetic completeness. Define a Six Sigma defect and opportunity +before collecting DPMO; dependent test assertions are not independent customer +opportunities. No Six Sigma capability is inferred from these tests. + +Budget independent datasets, not just taxon-level observations. Record model +misspecification, selection/refusal bias, Monte-Carlo uncertainty, software drift +and human applicability review separately. Each release claim needs an owner, +acceptance criterion fixed before measurement, receipt, expiry/retest trigger, +remaining uncertainty and a rollback/withdrawal decision. + +The comprehensive independent-review prompt and 18-by-14 ledger are maintained +in the review packet. They must remain capable of rejecting this assessment. From 0a400bbb056443017017f4757b3ef1f9d445fca8 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 20:49:37 +0000 Subject: [PATCH 7/9] Record repaired-reference calibration failure and preserved conformance counterexample The completed 32-dataset investigation found 19 false-discovery events in converged repaired-R fits. Of 29 paired native successes, 28 agree closely and replicate 27 has a q-value gap about 0.5041. Preserve both findings as unresolved launch blockers; green experiment CI establishes execution, not calibration or universal conformance. --- EXPLAINME.adoc | 10 ++++++++-- docs/LAUNCH-READINESS.adoc | 2 +- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/EXPLAINME.adoc b/EXPLAINME.adoc index b9fa97e..df4a606 100644 --- a/EXPLAINME.adoc +++ b/EXPLAINME.adoc @@ -33,8 +33,14 @@ datasets per regime before measurement. Of 384 method fits, 365 returned and counts. Under the small-unbalanced complete null, original ANCOM-BC made false discoveries in 19/29 successful fits at 0.05 (Wilson 95% interval 0.473–0.801). The unconditional observed count is 19/32; refusals are not successful negative -findings. This requires investigation against the repaired R reference and a -reviewed support restriction or estimator correction before a calibration claim. +findings. The follow-up repaired-R experiment converged for all 32 datasets and +also made false discoveries in 19. Of 29 paired successful native fits, 28 had +maximum q-value differences below 8.5e-6, but replicate 27 differed by 0.5041 +and changed the discovery decision. Numerical conformance is therefore not +established across this expanded regime. A reviewed support restriction or +estimator correction, and fresh prespecified validation, are needed before a +calibration claim. The raw receipt is CI run 42411431729, artifact +`small-unbalanced-calibration`, at source `cbb8ac4d85b4259c8122ed83c2cce8ec6c34362a`. Do not fix it by tuning the pilot seeds, removing adverse regimes or quietly substituting another method. diff --git a/docs/LAUNCH-READINESS.adoc b/docs/LAUNCH-READINESS.adoc index dafb6cc..0a3a9e9 100644 --- a/docs/LAUNCH-READINESS.adoc +++ b/docs/LAUNCH-READINESS.adoc @@ -44,7 +44,7 @@ framework still documents a 17-by-14 schema. Preserve that migration gap. |RobotRepoAutomaton |Execute authorised repairs with regression evidence and preserve provenance |May implement fixes; cannot approve its own scientific or human-review gate. |Panicbot |Escalate severe integrity/security failures without suppressing evidence -|Adverse small-unbalanced ANCOM-BC calibration must block a general calibration claim. +|Adverse small-unbalanced ANCOM-BC calibration occurs in native and repaired R results. Expanded comparison also finds one large q-value disagreement. Both remain launch blockers for the corresponding claims. |Human maintainer and independent scientist |Review actual diff by hand, assess applicability, data rights, uncertainty and accepted claims |Pending; automated or neural reviews do not count as those decisions. |=== From 41c4d66b77898fe9842d19ea42f485e1db4fff27 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 20:59:41 +0000 Subject: [PATCH 8/9] Update benchmark/README.md Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> --- benchmark/README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/benchmark/README.md b/benchmark/README.md index 535f89e..383b37b 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -23,5 +23,6 @@ inputs and native results, then runs pinned ANCOMBC 2.12.0 with the explicitly identified variance-estimator repair. It records R convergence rather than treating an iteration-cap result as a valid fit. A successful investigation job means the experiment completed, not that the adverse calibration passed. -The raw artifact expires after 30 days; archive it with its SHA-256 before -using it in a release decision. +The workflow requests 30 days of artifact retention, but repository policy +may shorten it. Check the artifact expiry and archive it with its SHA-256 +before using it in a release decision. From ba228b2823f30656b5231300b7a7d7fb6a93374a Mon Sep 17 00:00:00 2001 From: "coderabbitai[bot]" <136622811+coderabbitai[bot]@users.noreply.github.com> Date: Sat, 7 Nov 2026 21:02:40 +0000 Subject: [PATCH 9/9] Capture resolved environments and source provenance for benchmark runs, reject dirty timing runs, and clarify calibration guidance --- .gitignore | 1 + EXPLAINME.adoc | 9 +++--- benchmark/README.md | 23 ++++++++++++-- benchmark/benchmarks.jl | 4 ++- benchmark/calibration_pilot.jl | 5 ++- benchmark/export_calibration_case.jl | 1 + benchmark/run_environment.jl | 26 ++++++++++++++++ benchmark/test_run_environment.jl | 46 ++++++++++++++++++++++++++++ 8 files changed, 107 insertions(+), 8 deletions(-) create mode 100644 benchmark/run_environment.jl create mode 100644 benchmark/test_run_environment.jl diff --git a/.gitignore b/.gitignore index 16b9253..f07ef2d 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ *.agdai _build/ /test/oracle/Manifest.toml +/benchmark/Manifest.toml diff --git a/EXPLAINME.adoc b/EXPLAINME.adoc index df4a606..1c06f1a 100644 --- a/EXPLAINME.adoc +++ b/EXPLAINME.adoc @@ -45,10 +45,11 @@ Do not fix it by tuning the pilot seeds, removing adverse regimes or quietly substituting another method. ALDEx2's pilot is much less adverse, but 32 datasets cannot certify 5% FDR. -All-null FDR is the probability of any rejection. With signal present, that -probability is family-wise error; mean false discovery proportion estimates -FDR instead. The reported Wilson intervals are for event probabilities, not -mean FDP, taxa treated as independent replications, or arbitrary populations. +All-null FDR is the probability of any rejection. With signal present, +family-wise error is the probability of at least one false rejection; mean +false discovery proportion estimates FDR. The reported Wilson intervals are +for event probabilities, not mean FDP, taxa treated as independent +replications, or arbitrary populations. At only two ALDEx2 Monte-Carlo instances, a separate 4,096-seed probe produced 58 infinite and three undefined effect estimates. Its p/q bounds still held. diff --git a/benchmark/README.md b/benchmark/README.md index 535f89e..2a7a4e9 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -14,7 +14,19 @@ julia --project=benchmark benchmark/benchmarks.jl /path/to/benchmark-output Read `docs/CALIBRATION-PILOT-PROTOCOL.adoc` before interpreting the output. Seeds, scenarios and comparisons were fixed before results. Refusals are separate from successful fits. Neither script certifies calibration, speed, + + + + stability or package readiness. Keep generated measurements outside source. +Each output directory includes the active environment's `Project.toml` and +resolved `Manifest.toml`, including transitive dependencies such as StatsFuns. +Retain both with the results; instantiate the active project before running. +Timing benchmarks refuse a known dirty source tree (including untracked files). +When Git provenance is unavailable, scripts warn and record `source_head` and +`source_tree` as `unknown`; those runs need a separate source archive to reproduce. + + Benchmark timings exclude compilation/warmup and retain trial counts; ten trials provide an exploratory baseline, not precise tail-latency estimates. @@ -23,5 +35,12 @@ inputs and native results, then runs pinned ANCOMBC 2.12.0 with the explicitly identified variance-estimator repair. It records R convergence rather than treating an iteration-cap result as a valid fit. A successful investigation job means the experiment completed, not that the adverse calibration passed. -The raw artifact expires after 30 days; archive it with its SHA-256 before -using it in a release decision. + + + + +The workflow requests 30 days of artifact retention; repository policy may +shorten that period. Reviewers should check the actual expiry before archiving +the raw artifact with its SHA-256 for a release decision. + + diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 8d2f243..50ba0a2 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -4,6 +4,8 @@ using BenchmarkTools, CompositionalDA, Statistics, TOML, Dates include("calibration_pilot.jl") destination=only(ARGS) mkpath(destination) +provenance=source_provenance(pkgdir(CompositionalDA);require_clean=true) +save_run_environment(destination) rows=NamedTuple[] for (n,d) in ((20,40),(20,400),(100,400)) regime=(name="benchmark",n1=n÷2,n2=n÷2,depth1=10000,depth2=10000,sd=.35,sparse=false,signals=4) @@ -35,7 +37,7 @@ open(joinpath(destination,"benchmarks.csv"),"w") do io end metadata=Dict("julia"=>string(VERSION),"threads"=>Threads.nthreads(),"cpu"=>Sys.CPU_NAME, "machine"=>Sys.MACHINE,"word_size"=>Sys.WORD_SIZE,"date_utc"=>string(now(UTC)), - "source_head"=>strip(read(`git -C $(@__DIR__) rev-parse HEAD`,String)), "benchmarktools"=>string(pkgversion(BenchmarkTools)),"warmup_excluded"=>true, "claim"=>"exploratory warm-process baseline; no reference speedup or acceptance budget") +merge!(metadata,provenance) open(io->TOML.print(io,metadata),joinpath(destination,"benchmarks-metadata.toml"),"w") diff --git a/benchmark/calibration_pilot.jl b/benchmark/calibration_pilot.jl index 81415cc..f21677c 100644 --- a/benchmark/calibration_pilot.jl +++ b/benchmark/calibration_pilot.jl @@ -4,6 +4,7 @@ using CompositionalDA, Random, Statistics, Dates, Printf, SHA, TOML using CompositionalDA: StableRNG using Distributions: Multinomial +include("run_environment.jl") const REGIMES = [ (name="balanced_null",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=0), @@ -40,13 +41,15 @@ end function run_pilot(destination) mkpath(destination) + provenance=source_provenance(pkgdir(CompositionalDA)) + save_run_environment(destination) protocol=joinpath(@__DIR__,"..","docs","CALIBRATION-PILOT-PROTOCOL.adoc") metadata=Dict("protocol_sha256"=>bytes2hex(sha256(read(protocol))), "julia"=>string(VERSION),"package_version"=>string(pkgversion(CompositionalDA)), - "source_head"=>strip(read(`git -C $(@__DIR__) rev-parse HEAD`,String)), "started_utc"=>string(now(UTC)),"replicates_per_regime"=>32, "mc_samples"=>128,"seed_formula"=>"20261107 + 10000*regime_index + replicate", "status"=>"exploratory pilot; not a calibration certificate") + merge!(metadata,provenance) records=NamedTuple[] failures=NamedTuple[] for (regime_index,regime) in enumerate(REGIMES), replicate in 1:32 diff --git a/benchmark/export_calibration_case.jl b/benchmark/export_calibration_case.jl index 14b87c3..3bee2a2 100644 --- a/benchmark/export_calibration_case.jl +++ b/benchmark/export_calibration_case.jl @@ -3,6 +3,7 @@ include("calibration_pilot.jl") destination=only(ARGS) mkpath(destination) +save_run_environment(destination) for replicate in 1:32 x,features,groups,truth=simulate(REGIMES[3],20261107+30000+replicate) open(joinpath(destination,"replicate-$replicate.tsv"),"w") do io diff --git a/benchmark/run_environment.jl b/benchmark/run_environment.jl new file mode 100644 index 0000000..15ef3cf --- /dev/null +++ b/benchmark/run_environment.jl @@ -0,0 +1,26 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +using Pkg + +function source_provenance(source; require_clean=false) + head, state = try + head = strip(read(`git -C $source rev-parse HEAD`, String)) + status = read(`git -C $source status --porcelain --untracked-files=all`, String) + (head, isempty(status) ? "clean" : "dirty") + catch error + @warn "Source revision could not be captured; Git provenance is unavailable" exception=(error, catch_backtrace()) + ("unknown", "unknown") + end + require_clean && state == "dirty" && + error("Refusing to benchmark a dirty source tree; commit source changes and keep outputs outside the repository") + return Dict("source_head" => head, "source_tree" => state) +end + +function save_run_environment(destination) + # EnvCache selects the manifest used by Pkg, including version-specific names. + environment = Pkg.Types.EnvCache() + isfile(environment.manifest_file) || + error("A resolved manifest is required; instantiate the active project before running") + cp(environment.project_file, joinpath(destination, "Project.toml"); force=true) + cp(environment.manifest_file, joinpath(destination, "Manifest.toml"); force=true) +end diff --git a/benchmark/test_run_environment.jl b/benchmark/test_run_environment.jl new file mode 100644 index 0000000..4b6aae3 --- /dev/null +++ b/benchmark/test_run_environment.jl @@ -0,0 +1,46 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +using Test, TOML +include("run_environment.jl") + +@testset "Run provenance" begin + mktempdir() do source + run(`git -C $source init -q`) + write(joinpath(source,"source.jl"),"original\n") + run(`git -C $source add source.jl`) + run(`git -C $source -c user.name=Test -c user.email=test@example.invalid commit -qm baseline`) + head=strip(read(`git -C $source rev-parse HEAD`,String)) + @test source_provenance(source;require_clean=true)==Dict("source_head"=>head,"source_tree"=>"clean") + write(joinpath(source,"source.jl"),"modified\n") + @test source_provenance(source)["source_head"]==head + @test_throws ErrorException source_provenance(source;require_clean=true) + run(`git -C $source checkout -- source.jl`) + write(joinpath(source,"untracked.jl"),"new source\n") + @test_throws ErrorException source_provenance(source;require_clean=true) + withenv("PATH"=>"") do + provenance=@test_logs (:warn,r"Source revision could not be captured") source_provenance(source;require_clean=true) + @test provenance==Dict("source_head"=>"unknown","source_tree"=>"unknown") + end + end + mktempdir() do source + provenance=@test_logs (:warn,r"Source revision could not be captured") source_provenance(source) + @test provenance["source_head"]=="unknown" + end +end + +@testset "Resolved environment snapshot" begin + mktempdir() do destination + save_run_environment(destination) + environment=Pkg.Types.EnvCache() + @test read(joinpath(destination,"Project.toml"))==read(environment.project_file) + @test read(joinpath(destination,"Manifest.toml"))==read(environment.manifest_file) + manifest=TOML.parsefile(joinpath(destination,"Manifest.toml")) + @test haskey(only(manifest["deps"]["StatsFuns"]),"version") + end + mktempdir() do project + Pkg.activate(project) do + write(joinpath(project,"Project.toml"),"[deps]\n") + @test_throws ErrorException save_run_environment(project) + end + end +end