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/5] 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/5] 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/5] 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/5] 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/5] 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