Correct ALDEx rank tails and harden numerical boundaries after fuzz review

.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: |

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:

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

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

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

proofs/agda/CompositionalDA/DiscreteTails.agda

@@ -0,0 +1,27 @@
+-- SPDX-License-Identifier: MPL-2.0
+-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+{-# 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

proofs/agda/reject/LostAtom.agda

@@ -0,0 +1,9 @@
+-- SPDX-License-Identifier: MPL-2.0
+-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+-- EXPECT: 0 != 1
+{-# OPTIONS --safe --without-K #-}
+module reject.LostAtom where
+open import Data.Nat using (_∸_)
+open import Relation.Binary.PropositionalEquality using (_≡_; refl)
+false-inclusive-lower : 6 ∸ 6 ≡ 1
+false-inclusive-lower = refl

src/Core/bh.jl

@@ -14,14 +14,18 @@ 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)
     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

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,22 @@ 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.
+    # Widen narrow floats before division, while leaving large integers and
+    # arbitrary-precision values in their original scaling domain.
+    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) ||
+        throw(CompositionError("closed proportions are not representable as positive Float64 values"))
+    return out
 end
 
 """
@@ -81,6 +108,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)

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,56 @@ 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
+        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 * 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