diff --git a/README.md b/README.md index f3f6593..6ef7e33 100644 --- a/README.md +++ b/README.md @@ -34,9 +34,10 @@ q = bh_adjust([0.01, 0.04, 0.2]) # equals R's p.adjust(, "BH") `proofs/agda` holds machine-checked proofs (Agda 2.7.0.1, `--safe --without-K`, no postulates) of the algebraic properties the methods rely on. Each CLR theorem is mirrored by a Julia property test. Each property test is in turn shown to reject -a planted mutant (`test/mutants.jl`). The Benjamini–Hochberg proofs cover the -scaling step only, and its Julia property test checks a different property; the -table says which is which. Statistical properties such as calibration, +a planted mutant (`test/mutants.jl`). The Benjamini–Hochberg step-up (P5: q ≥ p, +q ≤ 1, order-preserving, at the original indices) is proved over exact +rationals, assuming the input was sorted correctly; the sort itself is not +verified, and nothing about floating point is proved. Statistical properties such as calibration, consistency and FDR control are **not** proved. They are to be tested against the R reference implementations; that oracle is not built yet. The full table is in [`proofs/PROOF-STATUS.md`](proofs/PROOF-STATUS.md). diff --git a/proofs/PROOF-STATUS.md b/proofs/PROOF-STATUS.md index 9c4766e..90ae738 100644 --- a/proofs/PROOF-STATUS.md +++ b/proofs/PROOF-STATUS.md @@ -1,62 +1,265 @@ - - + # Proof status -What is **proved** (machine-checked in Agda), what is only **tested** (Julia -property and unit tests), and what is still **open**. A row says only what its -column names; nothing else is implied. - -Gate: Agda 2.7.0.1, `--safe --without-K`, no postulates, agda-stdlib pinned by SHA -(`proofs/bootstrap.sh`). `proofs/tests/axiom-audit.sh` checks every module is -`--safe`, has no postulate, FFI, unsound flag, trusted primitive or hole, and is -reachable from `CompositionalDA/All.agda`. `proofs/tests/gate-selftest.sh` shows -the gate rejects the three controls in `proofs/agda/reject/` for their stated -reasons, and that each of the seven audit checks fires on a planted defect. - -## Scope of the proofs - -The CLR theorems are stated over an abstract commutative ring with a -logarithm-like homomorphism and a reciprocal `ι` of the dimension `n` -(`n · ι = 1`). They are proved about **exact algebra**. The Julia code computes -in `Float64`, so its tests check the same properties **up to a tolerance**. The -proofs say the definitions are right. They do not say the floating-point code -equals the definitions. - -## Core layer - -| Property | Proved (Agda) | Tested (Julia) | Mutant the test rejects | -|---|---|---|---| -| P1a — CLR components sum to zero | `CLR.clr-sum-zero` | `property_clr.jl` "P1a" | `mutant_clr_uncentred`, `mutant_clr_first_ref` | -| P1b — CLR is scale-invariant | `CLR.clr-scale-invariant` (via `clr-perturb`, `clr-constant`) | "P1b" | `mutant_clr_uncentred` | -| P1c — CLR is permutation-equivariant | `CLR.clr-permute` | "P1c" | `mutant_clr_position_weighted` | -| CLR turns perturbation into addition | `CLR.clr-perturb` | not tested separately | — | -| CLR is **not** invariant under an arbitrary perturbation (n = 2) | `CLR.Instance.clr-not-shift-invariant` | — | — | -| `LogHom` and `Reciprocal` are jointly inhabited for every n ≥ 1, so the CLR theorems are not vacuous (an exact-rational model, not real logarithms) | `CLR.Instance.additive-log`, `CLR.Instance.reciprocal` | — | — | -| BH scaling step is exactly `M·n / ((j+1)(d+1))` | `BenjaminiHochberg.bhScale-numerator`, `bhScale-denominator` | — | — | -| BH q-value is monotone in the family size `M` | `BenjaminiHochberg.monotone-in-family-size` | — | — | -| P5 — BH output ≥ raw, ≤ 1, order-preserving | **not proved** | `property_bh.jl` "P5" | `mutant_bh_no_cummin` | -| BH equals R `p.adjust(, "BH")` | not proved | `unit_core.jl` known answers | — | - -P5 and the Agda BH theorems are **different properties**. The Agda file proves -facts about one scaling step; P5 is about the whole adjusted vector, including -the running minimum. Neither substitutes for the other. - -## Open obligations - -- **BH non-negativity** of the scaled value. Mechanical, not conceptual: - `+ M * + n` does not reduce to `+ (M ℕ.* n)`, so the stdlib non-negativity - lemmas do not apply without a transport. It is not asserted anywhere. -- **CLR beyond algebra**: nothing about the real numbers, `exp`, continuity, - floating point, or the geometry of the simplex is proved. -- **The running minimum** (step-up) of BH is not modelled in Agda. -- **Toolchain**: the proofs need agda-stdlib on the 3.0 development line, pinned - by SHA in `bootstrap.sh`. They do not compile against v2.4 or v2.1. - -## Methods - -Not proved, and not claimed. Calibration, consistency and FDR control are -statistical properties, and the plan is to test them against the Bioconductor -reference implementations (ALDEx2, ANCOMBC). That oracle does **not exist yet**: -there is no RCall extension and no `test/oracle_*.jl`. `src/Methods/ALDEx2.jl` is -work in progress and is not yet included in the module. +Machine-checked proofs for CompositionalDA.jl, in **Agda 2.7.0.1** with +**agda-stdlib pinned at `2ffa8b7d4e8e818717ad643d184f055a4d1b0447`**, under +`--safe --without-K`. Three words are kept apart throughout this file: + +- **proved** — type-checked by Agda with zero postulates, zero holes and + `--safe` on every module (the axiom audit fails the build otherwise). A + theorem may take *named assumptions* as hypotheses or module parameters; the + table says which. +- **tested** — exercised by a Julia property test in `test/`, which is also + shown to reject a planted mutant in `test/mutants.jl`. +- **not done** — neither of the above. Listed, not implied. + +> **The stdlib pin is a development SHA, not a release.** The branch's own +> library file declares `name: standard-library-3.0` and no `v3.0` tag exists. +> `proofs/bootstrap.sh` normalises the name and refuses any other revision. + +## Gates + +| Gate | Command | Result | +| --- | --- | --- | +| Type-check every proof, then audit | `proofs/bootstrap.sh` | PASS — exit 0, no warnings | +| Axiom audit on its own | `proofs/tests/axiom-audit.sh` | PASS — 7/7 modules reachable from `CompositionalDA.All`, no postulates, no FFI, no unsound flags, no holes, all `--safe` | +| Gate self-test (does the gate reject anything?) | `proofs/tests/gate-selftest.sh` | PASS — see "Last verified" below for the count | + +The CI job `proofs` in `.github/workflows/ci.yml` runs exactly those first and +third commands with `AGDA_BIN` and `Agda_datadir` set to the pinned release +tarball. `bootstrap.sh` vendors the stdlib at the pinned SHA into +`proofs/.vendor` (git-ignored) and writes `proofs/.agda-libraries`; the +self-test reads that file, so the two scripts cannot disagree about which +library they checked against. + +The exact commands, exit codes and counts of the most recent run are in +"Last verified" at the end of this file. + +## What is proved + +Every module is reachable from `proofs/agda/CompositionalDA/All.agda`. + +### P1 — the centred log-ratio transform + +Module `CompositionalDA.CLR`, parametrised by a commutative ring `R`, a +commutative monoid `M` (the compositions; `×ₘ` is Aitchison's perturbation), +a `LogHom L : M → R` and a `Reciprocal ρ` for the part count `n`. With +`clr x i = log (x i) − ι · Σⱼ log (x j)`: + +| Theorem | Statement | Named assumptions | Julia test mirrored | +| --- | --- | --- | --- | +| `clr-sum-zero` | `Σᵢ clr x i ≈ 0#` | `R`, `L`, `ρ` | `test/property_clr.jl` "P1a clr sums to zero" (`prop_clr_sums_to_zero`; mutants `mutant_clr_uncentred`, `mutant_clr_first_ref`) | +| `clr-scale-invariant` | `clr (scale k x) i ≈ clr x i` for every `k : Pos` | `R`, `L`, `ρ` | "P1b clr scale invariant" (`prop_clr_scale_invariant`; mutant `mutant_clr_uncentred`) | +| `clr-permute` | `clr (rearrange (π ⟨$⟩ʳ_) x) i ≈ clr x (π ⟨$⟩ʳ i)` for every permutation `π` | `R`, `L`, `ρ` | "P1c clr permutation equivariant" (`prop_clr_permutation_equivariant`; mutant `mutant_clr_position_weighted`) | +| `clr-perturb` | `clr (perturb x y) i ≈ clr x i + clr y i` | `R`, `L`, `ρ` | — (lemma for P1b; not separately tested) | +| `clr-constant` | `clr (replicate n k) i ≈ 0#` | `R`, `L`, `ρ` | — (lemma for P1b) | +| `centre-sum-zero` | `Σ (centre v) ≈ 0#` | `R`, `ρ` | — (lemma for P1a) | + +Module `CompositionalDA.CLR.Instance` shows the assumptions are satisfiable, +so the theorems above are not vacuous, and supplies the negative control: + +| Theorem | Statement | Named assumptions | +| --- | --- | --- | +| `additive-log` | `LogHom +-*-commutativeRing +-0-commutativeMonoid` — the monoid `(ℕ, +)`, the ring `ℚᵘ`, `log a = mkℚᵘ (+ a) 0` (`a/1`) | none | +| `reciprocal` | `Reciprocal +-*-commutativeRing (suc m)` for every `m`, `ι = mkℚᵘ (+ 1) m` (`1/(m+1)`) | none | +| `clr-not-shift-invariant` | `¬ (CLR₂.clr (CLR₂.perturb x₀ y₀) 0F ≃ CLR₂.clr x₀ 0F)` at `n = 2` | none — this is why `clr-scale-invariant` needs the perturbing composition constant | + +### P5 — Benjamini–Hochberg step-up + +Module `CompositionalDA.BenjaminiHochberg.StepUp` models the loop in +`src/Core/bh.jl` on exact rationals `ℚᵘ`: + +``` +scale m k p = p · m/(k+1) -- q[k] * m / rank, rank = k+1 +stepUp m run (p ∷ ps) = (run ⊓ scale m k p) ∷ stepUp m (run ⊓ scale m k p) ps +adjust {m} = stepUp m 1ℚᵘ -- running = 1.0, m = length +``` + +Theorems on the **sorted** sequence (what the loop sees): + +| Theorem | Statement | Named assumptions | +| --- | --- | --- | +| `adjust-≤1` | `lookup (adjust ps) i ≤ 1` | none | +| `adjust-descending` | `Linked _≥_ (adjust ps)` — the output never rises along the sorted order | none | +| `adjust-≥` | `lookup ps i ≤ lookup (adjust ps) i` — **q ≥ p** | input descending; every p in `[0, 1]` | +| `adjust-monotone` | `lookup ps i ≤ lookup ps j → lookup (adjust ps) i ≤ lookup (adjust ps) j` — **p_i ≤ p_j ⇒ q_i ≤ q_j**, at any two positions, ties included | input descending; every p `≥ 0` | +| `stepUp-tie` | positions holding the same p-value as the head receive the head's adjusted value | descending, `≥ 0` | +| `scale-≥` | `p ≤ scale m k p` when `k+1 ≤ m` and `0 ≤ p` | as stated | +| `scale-antitone` | `scale m (k+1) p ≤ scale m k p` when `0 ≤ p` | as stated | +| `stepUp-≤run`, `stepUp-descending`, `stepUp-≥`, `stepUp-tie-from`, `stepUp-monotone`, `desc-lookup`, `1≤ratio`, `ratio-antitone` | the generalised forms the four `adjust-*` theorems are instances of | as stated in the module | + +Theorems at the **original** indices (what `test/property_bh.jl` reads), in +the sub-module `Unsorted p π` where `π : Permutation m m` stands for +`ord = sortperm(q; rev = true)` and `bh k = lookup (adjust sorted) (π ⟨$⟩ˡ k)` +is the `adj[k] = running` write-back: + +| Theorem | Statement | Named assumptions | Julia test mirrored | +| --- | --- | --- | --- | +| `bh-≤1` | `bh k ≤ 1` | none | `test/property_bh.jl` "P5 bh monotone" (`prop_bh_monotone`: `q .<= 1`) | +| `bh-≥` | `lookup p k ≤ bh k` | `Linked _≥_ sorted`; `∀ k → 0 ≤ p k`; `∀ k → p k ≤ 1` | same (`q .>= p .- 1e-15`); also the known-answer vector in `test/unit_core.jl` (`bh_adjust([0.01, 0.04, 0.03, 0.005, 0.5, 0.04, NaN, 0.2])`) | +| `bh-monotone` | `lookup p i ≤ lookup p j → bh i ≤ bh j` | `Linked _≥_ sorted`; `∀ k → 0 ≤ p k` | same (`q[i] <= q[j] + 1e-15` whenever `p[i] <= p[j]`); mutant `mutant_bh_no_cummin` | +| `sorted-lookup` | `lookup sorted (π ⟨$⟩ˡ k) ≡ lookup p k` | none | — | + +Module `CompositionalDA.BenjaminiHochberg` keeps the earlier per-rank +scaling as a single `mkℚᵘ`: + +| Theorem | Statement | Named assumptions | Julia test mirrored | +| --- | --- | --- | --- | +| `bhScale-numerator`, `bhScale-denominator` | the numerator is `M·n`, the denominator `(j+1)(d+1)` (by `refl`) | none | — | +| `monotone-in-family-size` | `M ≤ M′ → bhScale M j n d ≤ bhScale M′ j n d` | `NonNegative n` | — (no Julia test varies the family size) | + +### Support + +`CompositionalDA.Prelude` (exact-rational lemmas: `common-denominator-+`, +`whole-is-one`, `zero-over`, `numerator-monotone-≤`, `≡⇒≃`, `+-cong-≃ʳ`) and +`CompositionalDA.LogHom` (the two assumption records) carry no claims of +their own. + +## Named assumptions + +Everything a theorem above takes as given, and how each is discharged — or +not. "Discharged by Julia" means a line of Julia code establishes it at run +time and that line is tested, not proved. + +| Assumption | Where it appears | Discharged by | +| --- | --- | --- | +| `LogHom R M` — `log (a ×ₘ b) ≈ log a + log b` | module parameter of `CLR` | Proved satisfiable by `CLR.Instance.additive-log`. That the real logarithm on positive reals is one is a fact about ℝ, **not formalised**. | +| `Reciprocal R n` — `n · ι ≈ 1` | module parameter of `CLR` | Proved satisfiable by `CLR.Instance.reciprocal` for every `n ≥ 1`. `n = 0` has no reciprocal and no composition. | +| `Linked _≥_ sorted` — π sorts `p` descending | hypothesis of `Unsorted.bh-≥`, `Unsorted.bh-monotone` | Julia: `ord = sortperm(q; rev = true)` in `bh_adjust`. **Not proved**; `sortperm` is not re-implemented in Agda (`residue/benjamini-hochberg.residue` R-BH-3). | +| `0 ≤ p ≤ 1` for every p-value | hypotheses of `adjust-≥`, `adjust-monotone` and the `Unsorted` theorems | Julia: `bh_adjust` raises `ArgumentError` for any finite value outside `[0, 1]` before the loop runs; tested by `@test_throws ArgumentError bh_adjust([0.5, 1.5])` in `test/unit_core.jl`. **Not proved** (R-BH-4). The hypotheses are load-bearing: without `0 ≤ p` the tie argument fails; without `p ≤ 1`, `running = 1` undercuts the largest p-value. | +| `NonNegative n` | `monotone-in-family-size` | A p-value numerator is a count. Stated, not derived. | +| Rank bound `k + 1 ≤ m` | internal to `stepUp-≥` | **Proved** at the call site (`adjust-≥` passes `≤-refl`; the recursion passes `<⇒≤`). Not an assumption of any exported theorem. | + +## Tested, not proved + +Each item here is covered by a Julia test and by nothing in `proofs/`. + +- **`sortperm` and the write-back.** That `ord` is a permutation sorting the + non-NaN p-values descending, and that `adj[k] = running` writes sorted + position `r` to original index `ord[r]` — i.e. that the Julia `ord` is the + Agda `π` and `adj[k]` is `bh k`. Tested by `prop_bh_monotone` at the + original indices. +- **NaN exclusion.** `ℚᵘ` has no NaN. The Agda model is the loop over the + NaN-free vector with `m` its length; that NaN entries are excluded, that `m` + counts only the rest, and that NaN-in gives NaN-out are Julia and tested: + the generator `pvectors` in `test/properties.jl` plants a NaN in roughly one + vector in five, `prop_bh_monotone` asserts + `all(isnan, q[findall(isnan, p)])`, and `test/unit_core.jl` asserts + `all(isnan, bh_adjust([NaN, NaN]))` and a known-answer vector containing a + NaN. +- **The `[0, 1]` refusal.** The `ArgumentError` in `bh_adjust` is what makes + the Agda hypotheses true at run time. The guard is tested + (`@test_throws ArgumentError bh_adjust([0.5, 1.5])` in `test/unit_core.jl`); + the proofs assume its postcondition. +- **Float64 ↔ `ℚᵘ`.** `q[k] * m / rank` in Float64 is a rounding of + `scale m k p`; `min` on floats is `_⊓_` on rationals only up to that + rounding. IEEE-754 is not modelled. The `1e-15` slack in `prop_bh_monotone` + (`q .>= p .- 1e-15`, `q[i] <= q[j] + 1e-15`) is the test's allowance for + this, not a proved bound. Likewise the real-valued `clr` in `src/` against + the parametric theorems: `test/property_clr.jl` checks numerically with a + tolerance. +- **Ties in floating point.** `adjust-monotone` covers ties exactly on + rationals; whether two Float64 p-values that are equal stay equal through + `* m / rank` is a floating-point question and is tested only. The generator + does plant ties: `pvectors` rounds to 1–3 decimal digits (up to 40 entries) + and, in the vectors where it does not plant a NaN, copies `p[1]` over a + random position. +- **`monotone-in-family-size`** is proved but has **no Julia mirror**: no + property test varies `m`. +- **R reference agreement** (ALDEx2/ANCOMBC through the optional RCall + extension) is a test oracle only. + +## Not done + +- **FDR control.** That the q-values control the false discovery rate at the + nominal level is a theorem about a probability measure over a sampling + model. Nothing in `proofs/` bears on it and nothing in `test/` tests it + against a ground truth; the README says the same. R-BH-7. +- **`bhScale` non-negativity in the expanded-denominator form** (R-BH-1): + open for `bhScale`; the product form `scale` has it via `scale-≥`. + `scale-is-bhScale`, which would connect the two modules, is not written. +- **A verified sort** replacing the `Linked _≥_ sorted` assumption (R-BH-3). +- **Real numbers.** No `LogHom` instance over ℝ, no `exp`, no continuity, no + simplex geometry (R-CLR-1, R-CLR-5). +- **Floating-point models** of either method (R-CLR-4, R-BH-6). + +The per-obligation records, with status and the shape of a closing proof, +are in `proofs/residue/clr.residue` and +`proofs/residue/benjamini-hochberg.residue`. + +## The gate's own evidence + +`proofs/tests/gate-selftest.sh` mutates a throwaway copy of the tree and +requires the gate to reject each mutation; a mutation whose `sed` pattern no +longer matches is a failure, not a pass. Each audit control must also fail with +the finding its row names, so one audit check cannot stand in for another; there +is one control per check in `axiom-audit.sh` (seven). The controls: + +| Control | What is planted | Must be rejected by | +| --- | --- | --- | +| clr: sum-to-zero claimed as sum-to-one | `≈ 0#` → `≈ 1#` in `clr-sum-zero` | type-checker | +| clr: centring removed | `clr x = map log x` | type-checker | +| bh: family-size monotonicity reversed | `≤ℚ` operands swapped | type-checker | +| stepup: running minimum dropped | `stepUp` emits `scale m k p` without `⊓` (mirrors `mutant_bh_no_cummin`) | type-checker | +| stepup: q ≥ p reversed | conclusion of `adjust-≥` flipped | type-checker | +| stepup: monotonicity reversed | conclusion of `adjust-monotone` flipped | type-checker | +| stepup: scale written as rank/m | `mkℚᵘ (+ m) k` → `mkℚᵘ (+ suc k) m` | type-checker | +| audit: module dropped from the gate entry | `StepUp` import deleted from `All.agda` | axiom audit | +| audit: postulate injected | `postulate cheat : ∀ {A : Set} → A` | axiom audit | +| audit: `--safe` removed | from `CLR.agda` | axiom audit | +| audit: FFI pragma injected | `{-# FOREIGN GHC … #-}` appended to `Prelude.agda` | axiom audit | +| audit: unsound flag added | `--type-in-type` beside `--safe` in `Prelude.agda` | axiom audit | +| audit: trustMe mentioned | a `trustMe` reference appended to `Prelude.agda` | axiom audit | +| audit: hole left in a module | `hole = {! !}` appended to `Prelude.agda` | axiom audit | +| reject/FakeCLR, reject/MissingReciprocal, reject/Postulate | the standing negative controls | type-checker, with exit 42 **and** the error named on each file's `-- EXPECT:` line | +| pristine tree | nothing | must be **accepted** | + +## Reproducing + +```sh +proofs/bootstrap.sh # vendor the stdlib at the pin, type-check, audit +proofs/tests/gate-selftest.sh # prove the gate rejects planted defects +proofs/tests/axiom-audit.sh # the audit alone +``` + +`bootstrap.sh` needs exactly `Agda version 2.7.0.1` as `AGDA_BIN` or `agda` +on `PATH`, with `Agda_datadir` pointing at the release tarball's `data/` +directory; it prints the install steps and fails if the version differs. It +does not install Agda (no package manager, no Python) and never writes into +the toolchain's data directory. **An absent prover is a failure, never a +skip** — in CI and locally alike. + +## Last verified + +2026-10-02, one run of the unmodified CI interface, from an **empty** +`proofs/.vendor` (so the stdlib was fetched at the pin and type-checked from +scratch), with the release binary the CI job downloads: + +``` +AGDA_BIN=…/Agda-v2.7.0.1/bin/agda Agda_datadir=…/Agda-v2.7.0.1/data +proofs/bootstrap.sh → exit 0 (7 modules reachable from CompositionalDA.All; axiom-audit: clean) +proofs/tests/gate-selftest.sh → exit 0 (gate-selftest: 18/18 controls behaved correctly) +proofs/tests/axiom-audit.sh → exit 0 (auditing 7 module(s); axiom-audit: clean) +``` + +Wall-clock 14:18:26Z → 14:22:32Z, almost all of it the stdlib. The self-test +was re-run after rebasing onto PR #1, which added the four remaining audit +controls; a planted wrong reason on one of them was reported as a FAIL +(17/18), so the reason check is live. The 18 +controls are the 14 mutations, the 3 `reject/` files and the pristine tree +listed above; each mutation was confirmed to have changed its file before the +gate was run, and each `reject/` file was refused with exit 42 and the error +its `-- EXPECT:` line names. + +The Julia side, same day: + +``` +julia --project -e 'using Pkg; Pkg.test()' → exit 0 (Test Summary: CompositionalDA | Pass 43 | Total 43) +``` + +That run includes `test/property_clr.jl` (P1a/P1b/P1c) and +`test/property_bh.jl` (P5), and the mutant checks in `test/mutants.jl`. diff --git a/proofs/agda/CompositionalDA/All.agda b/proofs/agda/CompositionalDA/All.agda index bb6c1b5..2411273 100644 --- a/proofs/agda/CompositionalDA/All.agda +++ b/proofs/agda/CompositionalDA/All.agda @@ -15,3 +15,4 @@ import CompositionalDA.LogHom import CompositionalDA.CLR import CompositionalDA.CLR.Instance import CompositionalDA.BenjaminiHochberg +import CompositionalDA.BenjaminiHochberg.StepUp diff --git a/proofs/agda/CompositionalDA/BenjaminiHochberg/StepUp.agda b/proofs/agda/CompositionalDA/BenjaminiHochberg/StepUp.agda new file mode 100644 index 0000000..9238c47 --- /dev/null +++ b/proofs/agda/CompositionalDA/BenjaminiHochberg/StepUp.agda @@ -0,0 +1,280 @@ +-- SPDX-License-Identifier: MPL-2.0 +-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +-- +-- Benjamini–Hochberg step-up, as the running minimum that `src/Core/bh.jl` +-- computes (P5). +-- +-- The Julia loop, after sorting the p-values into DESCENDING order, is +-- +-- running = 1.0 +-- for (r, k) in enumerate(ord) # r = 1 is the largest p-value +-- rank = m - r + 1 # m, m-1, …, 1 +-- running = min(running, q[k] * m / rank) +-- adj[k] = running +-- end +-- +-- `stepUp` below is that loop, on exact rationals. The rank of the head is +-- the number of p-values still to process, so it counts down `m, m-1, …, 1` +-- exactly as `m - r + 1` does, and `running` starts at `1ℚᵘ`. `adjust` fixes +-- the family size to the length of the vector, which is what `m = length(idx)` +-- does after the NaNs are dropped. +-- +-- What is proved about `adjust`, each as a statement about `lookup` at a +-- position of the sorted sequence, and then re-indexed through the sorting +-- permutation so that it reads at the ORIGINAL index the way the Julia +-- property test `test/property_bh.jl` reads it: +-- +-- * `adjust-≤1` — every adjusted value is ≤ 1. (no hypothesis) +-- * `adjust-descending`— the output is descending too. (no hypothesis) +-- * `adjust-≥` — q_i ≥ p_i. Needs: input descending, 0 ≤ p ≤ 1. +-- * `adjust-monotone` — p_i ≤ p_j ⇒ q_i ≤ q_j, at ANY two positions, +-- ties included. Needs: input descending, 0 ≤ p. +-- +-- The hypotheses are not formalities. `0 ≤ p ≤ 1` is the `ArgumentError` +-- guard in `bh_adjust`; without `0 ≤ p` the scaling `p · m/rank` runs the +-- wrong way, and without `p ≤ 1` the initial `running = 1` already undercuts +-- the largest p-value. "Input descending" is what `sortperm(q; rev = true)` +-- establishes; it is assumed here, not re-proved, and that assumption is +-- recorded in `proofs/PROOF-STATUS.md`. +-- +-- Ties are the part a test on distinct values would never exercise: equal +-- p-values occupy consecutive ranks and so receive DIFFERENT scaled values, +-- and the proof that they nevertheless receive the same adjusted value +-- (`stepUp-tie`) is where `0 ≤ p` is actually used. +-- +-- Everything is over `ℚᵘ`, the unnormalised rationals, so no gcd step and no +-- rounding sits between the arithmetic and the proof. `_⊓_` on `ℚᵘ` does +-- not reduce on variables, so every step below goes through the lattice +-- lemmas of `Data.Rational.Unnormalised.Properties`, never through `refl`. + +{-# OPTIONS --without-K --safe #-} + +module CompositionalDA.BenjaminiHochberg.StepUp where + +open import Data.Nat.Base as ℕ using (ℕ; zero; suc; z≤n; s≤s) +open import Data.Nat.Properties as ℕₚ using (n≤1+n; <⇒≤) +open import Data.Integer.Base as ℤ using (ℤ; +_; +≤+) +open import Data.Integer.Properties as ℤₚ using () +open import Data.Fin.Base using (Fin; zero; suc; toℕ) +open import Data.Fin.Permutation using (Permutation; _⟨$⟩ʳ_; _⟨$⟩ˡ_; inverseʳ) +open import Data.Vec.Base using (Vec; []; _∷_; lookup; tabulate) +open import Data.Vec.Properties using (lookup∘tabulate) +open import Data.Vec.Relation.Unary.All as All using (All; []; _∷_) +open import Data.Vec.Relation.Unary.All.Properties using (tabulate⁺) +open import Data.Vec.Relation.Unary.Linked using (Linked; []; [-]; _∷_) +open import Data.Sum.Base using (inj₁; inj₂) +open import Relation.Binary.PropositionalEquality + using (_≡_; refl; sym; trans; cong; subst; subst₂) +open import Data.Rational.Unnormalised.Base + using (ℚᵘ; mkℚᵘ; 0ℚᵘ; 1ℚᵘ; _≃_; _≤_; _≥_; _*_; _⊓_; *≤*; nonNegative) +open import Data.Rational.Unnormalised.Properties + using ( ≃-refl; ≃-sym; ≃-trans + ; ≤-refl; ≤-reflexive; ≤-trans; ≤-antisym + ; *-identityʳ; *-congʳ; *-monoʳ-≤-nonNeg + ; p⊓q≤p; p⊓q≤q; ⊓-glb; p≤q⇒p⊓q≃p ) + +------------------------------------------------------------------------ +-- The model + +-- `p · m/(k+1)`: the Julia `q[k] * m / rank` with `rank = k + 1`. `ℚᵘ` +-- stores a denominator-minus-one, so the fraction `m/(k+1)` is `mkℚᵘ (+ m) k`. +scale : ℕ → ℕ → ℚᵘ → ℚᵘ +scale m k p = p * mkℚᵘ (+ m) k + +-- The loop. `run` is Julia's `running`. The head of a vector of length +-- `suc k` has rank `suc k`, i.e. `scale m k`. +stepUp : (m : ℕ) → ℚᵘ → ∀ {k} → Vec ℚᵘ k → Vec ℚᵘ k +stepUp m run [] = [] +stepUp m run {suc k} (p ∷ ps) = + (run ⊓ scale m k p) ∷ stepUp m (run ⊓ scale m k p) ps + +-- `bh_adjust` on the sorted, NaN-free p-values: family size = length, +-- `running` starts at 1. +adjust : ∀ {m} → Vec ℚᵘ m → Vec ℚᵘ m +adjust {m} = stepUp m 1ℚᵘ + +------------------------------------------------------------------------ +-- The scaling factor + +-- `m/(k+1) ≥ 1` whenever the rank `k+1` does not exceed the family size. +1≤ratio : ∀ m k → suc k ℕ.≤ m → 1ℚᵘ ≤ mkℚᵘ (+ m) k +1≤ratio m k sk≤m = *≤* + (ℤₚ.≤-trans (ℤₚ.≤-reflexive (ℤₚ.*-identityˡ (+ suc k))) + (ℤₚ.≤-trans (+≤+ sk≤m) + (ℤₚ.≤-reflexive (sym (ℤₚ.*-identityʳ (+ m)))))) + +-- A smaller rank gives a larger factor: `m/(k+2) ≤ m/(k+1)`. +ratio-antitone : ∀ m k → mkℚᵘ (+ m) (suc k) ≤ mkℚᵘ (+ m) k +ratio-antitone m k = *≤* (ℤₚ.*-monoˡ-≤-nonNeg (+ m) (+≤+ (n≤1+n (suc k)))) + +-- Scaling a non-negative p-value by a factor ≥ 1 cannot shrink it. +scale-≥ : ∀ m k p → suc k ℕ.≤ m → 0ℚᵘ ≤ p → p ≤ scale m k p +scale-≥ m k p sk≤m 0≤p = + ≤-trans (≤-reflexive (≃-sym (*-identityʳ p))) + (*-monoʳ-≤-nonNeg p {{nonNegative 0≤p}} (1≤ratio m k sk≤m)) + +-- Moving a non-negative p-value one rank down (towards rank 1) scales it up. +scale-antitone : ∀ m k p → 0ℚᵘ ≤ p → scale m (suc k) p ≤ scale m k p +scale-antitone m k p 0≤p = + *-monoʳ-≤-nonNeg p {{nonNegative 0≤p}} (ratio-antitone m k) + +------------------------------------------------------------------------ +-- The running minimum, position by position + +-- The output never rises: each value is the previous one `⊓` something. +stepUp-descending : ∀ m run {k} (ps : Vec ℚᵘ k) → Linked _≥_ (stepUp m run ps) +stepUp-descending m run [] = [] +stepUp-descending m run (p ∷ []) = [-] +stepUp-descending m run {suc (suc k)} (p ∷ p′ ∷ ps) = + p⊓q≤p (run ⊓ scale m (suc k) p) (scale m k p′) + ∷ stepUp-descending m (run ⊓ scale m (suc k) p) (p′ ∷ ps) + +-- Nothing exceeds the initial `running`. +stepUp-≤run : ∀ m run {k} (ps : Vec ℚᵘ k) (i : Fin k) → + lookup (stepUp m run ps) i ≤ run +stepUp-≤run m run {suc k} (p ∷ ps) zero = p⊓q≤p run (scale m k p) +stepUp-≤run m run {suc k} (p ∷ ps) (suc i) = + ≤-trans (stepUp-≤run m (run ⊓ scale m k p) ps i) (p⊓q≤p run (scale m k p)) + +-- In a descending vector a later position holds a smaller-or-equal value. +desc-lookup : ∀ {k} {xs : Vec ℚᵘ k} → Linked _≥_ xs → + ∀ (i j : Fin k) → toℕ i ℕ.≤ toℕ j → lookup xs j ≤ lookup xs i +desc-lookup lnk zero zero _ = ≤-refl +desc-lookup {xs = x ∷ x′ ∷ xs} (x′≤x ∷ lnk) zero (suc j) _ = + ≤-trans (desc-lookup lnk zero j z≤n) x′≤x +desc-lookup {xs = x ∷ xs} (_ ∷ lnk) (suc i) (suc j) (s≤s i≤j) = + desc-lookup lnk i j i≤j + +-- q_i ≥ p_i, given: the rank of the head does not exceed the family size, +-- the head is below the incoming `running`, the input is descending and +-- non-negative. The last three are what the recursion has to re-establish +-- for the tail: the new `running` is ≥ the head, which is ≥ the next value. +stepUp-≥ : ∀ m run {k} (p : ℚᵘ) (ps : Vec ℚᵘ k) → suc k ℕ.≤ m → p ≤ run → + Linked _≥_ (p ∷ ps) → All (0ℚᵘ ≤_) (p ∷ ps) → + ∀ (i : Fin (suc k)) → lookup (p ∷ ps) i ≤ lookup (stepUp m run (p ∷ ps)) i +stepUp-≥ m run {k} p ps sk≤m p≤run _ (0≤p ∷ _) zero = + ⊓-glb p≤run (scale-≥ m k p sk≤m 0≤p) +stepUp-≥ m run {suc k} p (p′ ∷ ps) sk≤m p≤run (p′≤p ∷ lnk) (0≤p ∷ all) (suc i) = + stepUp-≥ m (run ⊓ scale m (suc k) p) p′ ps (<⇒≤ sk≤m) + (≤-trans p′≤p (⊓-glb p≤run (scale-≥ m (suc k) p sk≤m 0≤p))) + lnk all i + +-- Ties. If the value at position `j` is ≥ the head, then (the input being +-- descending) every value from the head to `j` equals the head, and the +-- running minimum does not move across that block: each later scaled value +-- is ≥ the head's scaled value (`scale-antitone`), which is already ≥ the +-- minimum. So the adjusted value at `j` IS the adjusted value at the head. +stepUp-tie : ∀ m run {k} (p : ℚᵘ) (ps : Vec ℚᵘ k) → + Linked _≥_ (p ∷ ps) → All (0ℚᵘ ≤_) (p ∷ ps) → + ∀ (j : Fin (suc k)) → p ≤ lookup (p ∷ ps) j → + lookup (stepUp m run (p ∷ ps)) j ≃ run ⊓ scale m k p +stepUp-tie m run p ps _ _ zero _ = ≃-refl +stepUp-tie m run {suc k} p (p′ ∷ ps) (p′≤p ∷ lnk) (0≤p ∷ 0≤p′ ∷ all) (suc j) p≤pⱼ = + ≃-trans (stepUp-tie m q p′ ps lnk (0≤p′ ∷ all) j (≤-trans p′≤p p≤pⱼ)) + (p≤q⇒p⊓q≃p q≤s′) + where + q : ℚᵘ + q = run ⊓ scale m (suc k) p + pⱼ≤p′ : lookup (p′ ∷ ps) j ≤ p′ + pⱼ≤p′ = desc-lookup lnk zero j z≤n + p′≃p : p′ ≃ p + p′≃p = ≤-antisym p′≤p (≤-trans p≤pⱼ pⱼ≤p′) + q≤s′ : q ≤ scale m k p′ + q≤s′ = ≤-trans (p⊓q≤q run (scale m (suc k) p)) + (≤-trans (scale-antitone m k p 0≤p) + (≤-reflexive (*-congʳ (≃-sym p′≃p)))) + +-- The tie lemma, started at any position `i` rather than at the head. +stepUp-tie-from : ∀ m run {k} (ps : Vec ℚᵘ k) → Linked _≥_ ps → All (0ℚᵘ ≤_) ps → + ∀ (i j : Fin k) → toℕ i ℕ.≤ toℕ j → lookup ps i ≤ lookup ps j → + lookup (stepUp m run ps) i ≃ lookup (stepUp m run ps) j +stepUp-tie-from m run (p ∷ ps) lnk all zero j _ p≤pⱼ = + ≃-sym (stepUp-tie m run p ps lnk all j p≤pⱼ) +stepUp-tie-from m run {suc k} (p ∷ ps) (_ ∷ lnk) (_ ∷ all) (suc i) (suc j) (s≤s i≤j) pᵢ≤pⱼ = + stepUp-tie-from m (run ⊓ scale m k p) ps lnk all i j i≤j pᵢ≤pⱼ + +-- p_i ≤ p_j ⇒ q_i ≤ q_j at any two positions. Either `j` is at or before +-- `i`, and the output is descending; or `i` is at or before `j`, in which +-- case the input (also descending) is constant between them and the tie +-- lemma makes the outputs equal. +stepUp-monotone : ∀ m run {k} (ps : Vec ℚᵘ k) → Linked _≥_ ps → All (0ℚᵘ ≤_) ps → + ∀ (i j : Fin k) → lookup ps i ≤ lookup ps j → + lookup (stepUp m run ps) i ≤ lookup (stepUp m run ps) j +stepUp-monotone m run ps lnk all i j pᵢ≤pⱼ with ℕₚ.≤-total (toℕ i) (toℕ j) +... | inj₁ i≤j = ≤-reflexive (stepUp-tie-from m run ps lnk all i j i≤j pᵢ≤pⱼ) +... | inj₂ j≤i = desc-lookup (stepUp-descending m run ps) j i j≤i + +------------------------------------------------------------------------ +-- `bh_adjust` on the sorted sequence + +-- Every adjusted value is ≤ 1. +adjust-≤1 : ∀ {m} (ps : Vec ℚᵘ m) (i : Fin m) → lookup (adjust ps) i ≤ 1ℚᵘ +adjust-≤1 {m} = stepUp-≤run m 1ℚᵘ + +-- The adjusted sequence is descending, as the raw one was. +adjust-descending : ∀ {m} (ps : Vec ℚᵘ m) → Linked _≥_ (adjust ps) +adjust-descending {m} = stepUp-descending m 1ℚᵘ + +-- q_i ≥ p_i for a descending sequence of p-values in [0, 1]. +adjust-≥ : ∀ {m} (ps : Vec ℚᵘ m) → Linked _≥_ ps → + All (0ℚᵘ ≤_) ps → All (_≤ 1ℚᵘ) ps → + ∀ (i : Fin m) → lookup ps i ≤ lookup (adjust ps) i +adjust-≥ [] _ _ _ () +adjust-≥ (p ∷ ps) lnk all (p≤1 ∷ _) i = stepUp-≥ _ 1ℚᵘ p ps ℕₚ.≤-refl p≤1 lnk all i + +-- p_i ≤ p_j ⇒ q_i ≤ q_j for a descending sequence of non-negative p-values. +adjust-monotone : ∀ {m} (ps : Vec ℚᵘ m) → Linked _≥_ ps → All (0ℚᵘ ≤_) ps → + ∀ (i j : Fin m) → lookup ps i ≤ lookup ps j → + lookup (adjust ps) i ≤ lookup (adjust ps) j +adjust-monotone {m} = stepUp-monotone m 1ℚᵘ + +------------------------------------------------------------------------ +-- `bh_adjust` at the ORIGINAL indices +-- +-- `ord = sortperm(q; rev = true)` is a permutation π: sorted position `r` +-- holds the raw p-value at original index `π ⟨$⟩ʳ r`. `adj[k] = running` +-- writes sorted position `r`'s value back to `ord[r]`, so the adjusted value +-- at original index `k` is the sorted-position value at `π ⟨$⟩ˡ k`. That π +-- really sorts `p` descending is the named assumption `sorted-desc`; the +-- Agda does not re-implement `sortperm`. + +module Unsorted {m} (p : Vec ℚᵘ m) (π : Permutation m m) where + + sorted : Vec ℚᵘ m + sorted = tabulate (λ r → lookup p (π ⟨$⟩ʳ r)) + + -- The adjusted p-value at original index `k`. + bh : Fin m → ℚᵘ + bh k = lookup (adjust sorted) (π ⟨$⟩ˡ k) + + -- Reading the sorted sequence at `k`'s sorted position gives `p` at `k`. + sorted-lookup : ∀ k → lookup sorted (π ⟨$⟩ˡ k) ≡ lookup p k + sorted-lookup k = + trans (lookup∘tabulate (λ r → lookup p (π ⟨$⟩ʳ r)) (π ⟨$⟩ˡ k)) + (cong (lookup p) (inverseʳ π)) + + -- q_k ≤ 1 at every original index. + bh-≤1 : ∀ k → bh k ≤ 1ℚᵘ + bh-≤1 k = adjust-≤1 sorted (π ⟨$⟩ˡ k) + + -- P5, second half: q_k ≥ p_k, given π sorts `p` descending and every + -- p-value lies in [0, 1]. + bh-≥ : Linked _≥_ sorted → + (∀ k → 0ℚᵘ ≤ lookup p k) → (∀ k → lookup p k ≤ 1ℚᵘ) → + ∀ k → lookup p k ≤ bh k + bh-≥ sorted-desc 0≤p p≤1 k = + subst (_≤ bh k) (sorted-lookup k) + (adjust-≥ sorted sorted-desc + (tabulate⁺ (λ r → 0≤p (π ⟨$⟩ʳ r))) + (tabulate⁺ (λ r → p≤1 (π ⟨$⟩ʳ r))) + (π ⟨$⟩ˡ k)) + + -- P5, first half: p_i ≤ p_j ⇒ q_i ≤ q_j at any two original indices, given + -- π sorts `p` descending and every p-value is non-negative. + bh-monotone : Linked _≥_ sorted → (∀ k → 0ℚᵘ ≤ lookup p k) → + ∀ i j → lookup p i ≤ lookup p j → bh i ≤ bh j + bh-monotone sorted-desc 0≤p i j pᵢ≤pⱼ = + adjust-monotone sorted sorted-desc + (tabulate⁺ (λ r → 0≤p (π ⟨$⟩ʳ r))) + (π ⟨$⟩ˡ i) (π ⟨$⟩ˡ j) + (subst₂ _≤_ (sym (sorted-lookup i)) (sym (sorted-lookup j)) pᵢ≤pⱼ) diff --git a/proofs/residue/benjamini-hochberg.residue b/proofs/residue/benjamini-hochberg.residue new file mode 100644 index 0000000..b3a498d --- /dev/null +++ b/proofs/residue/benjamini-hochberg.residue @@ -0,0 +1,125 @@ +# SPDX-License-Identifier: CC-BY-SA-4.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Residue — Benjamini–Hochberg (P5) +# +# See clr.residue for the record format. Two modules bear on this file: +# `CompositionalDA.BenjaminiHochberg` the per-rank scaling `bhScale` +# `CompositionalDA.BenjaminiHochberg.StepUp` the running-minimum step-up +# that `src/Core/bh.jl` computes +# +# The archived MetaManifold version of this file recorded "ℚᵘ has no ⊓ lemmas" +# and "q_i ≥ p_i is FALSE in general and must not be attempted". Both were +# true of what was being modelled there and are false of what is modelled +# here: the stdlib at the pinned SHA has the full `_⊓_` lattice for `ℚᵘ`, and +# `q_i ≥ p_i` is a theorem of the STEP-UP (running minimum from `running = 1`, +# rank ≤ family size, p ∈ [0,1]), which is what the Julia code computes. It +# remains false of the bare per-rank scaling, where `M/(j+1)` can be < 1. + +R-BH-1 + statement: | + ∀ M j n d → 0ℚᵘ ≤ℚ bhScale M j (+ n) d + (the scaled value of a non-negative p-value is non-negative, in the + `bhScale` form with the denominator expanded) + status: OPEN for `bhScale`; CLOSED for the product form by + `StepUp.scale-≥` (`0 ≤ p → p ≤ scale m k p`, which with `0 ≤ p` gives + `0 ≤ scale m k p` by `≤-trans`), under the rank bound `suc k ≤ m`. + why: | + `+ M * + n` does not reduce to `+ (M ℕ.* n)`; `ℤ._*_` unfolds to a + sign-times-absolute-value form, and the non-negativity lemmas in + `Data.Integer.Properties` are stated about the reduced form. The StepUp + module avoids the obstacle by keeping the scaling as a `ℚᵘ` product + (`p * mkℚᵘ (+ m) k`) and using `*-monoʳ-≤-nonNeg` with a `NonNegative` + instance built from the `0 ≤ p` hypothesis. + impact: LOW. A negative q-value would need a negative p-value, which + `bh_adjust` refuses with an `ArgumentError`. + closes-with: | + `scale-is-bhScale : scale m k (mkℚᵘ n d) ≃ bhScale m k n d` + (a `*-comm`/`+-comm` rearrangement of the denominator), after which + R-BH-1 follows from `scale-≥` without the rank bound being needed for + non-negativity alone. + +R-BH-2 + statement: | + The BH step-up running minimum: adjusted values are ≤ 1, descending in + rank, ≥ the raw p-value, and monotone (p_i ≤ p_j ⇒ q_i ≤ q_j) including + across ties. + status: CLOSED (by `StepUp.adjust-≤1`, `StepUp.adjust-descending`, + `StepUp.adjust-≥`, `StepUp.adjust-monotone`, and at original indices by + `StepUp.Unsorted.bh-≤1`, `StepUp.Unsorted.bh-≥`, + `StepUp.Unsorted.bh-monotone`) + why: n/a + impact: n/a + closes-with: n/a + +R-BH-3 + statement: | + `ord = sortperm(q; rev = true)` returns a permutation π such that the + sorted sequence is descending (`Linked _≥_ sorted` in + `StepUp.Unsorted`). + status: NOT ATTEMPTED (named assumption) + why: | + `sortperm` is not re-implemented in Agda. Every `Unsorted` theorem takes + `Linked _≥_ sorted` as an explicit hypothesis and so says: IF π sorts + descending THEN the property holds at the original indices. That + Julia's `sortperm` with `rev = true` produces such a π is tested, not + proved. + impact: | + LOW. `sortperm` is Base Julia and the property test + `test/property_bh.jl` evaluates the end-to-end result at the original + indices. + closes-with: | + A verified sort (e.g. insertion sort with a `Linked _≥_` output proof and + a permutation witness) substituted for π. Not planned. + +R-BH-4 + statement: | + Every p-value handed to the step-up lies in [0, 1]. + status: DISCHARGED BY A JULIA REFUSAL (not proved) + why: | + `bh_adjust` raises `ArgumentError` for any finite value outside [0, 1] + before the loop runs. The Agda hypotheses `All (0ℚᵘ ≤_)` and + `All (_≤ 1ℚᵘ)` are exactly what that guard establishes. The guard + itself is tested, not proved. + impact: | + The hypotheses are load-bearing: without `0 ≤ p` the scaling runs the + wrong way and `adjust-monotone` fails; without `p ≤ 1` the initial + `running = 1` undercuts the largest p-value and `adjust-≥` fails. + closes-with: n/a (a Julia-side guard cannot be proved in Agda). + +R-BH-5 + statement: | + NaN entries are excluded before the step-up and remain NaN afterwards; + `m` is the count of non-NaN entries. + status: OUT OF SCOPE + why: | + `ℚᵘ` has no NaN. The Agda model is the loop over the NaN-free vector, + with `m` its length. The exclusion and the write-back are Julia code, + tested by `prop_bh_monotone` (which asserts NaN-in ⇒ NaN-out). + impact: LOW. + closes-with: n/a + +R-BH-6 + statement: | + The floating-point `bh_adjust` agrees with the exact `adjust` to within + the `1e-15` tolerance the property test uses. + status: OUT OF SCOPE + why: | + IEEE-754 semantics are not modelled. `q[k] * m / rank` in Float64 is a + rounding of `scale m k p`; the tolerance in `test/properties.jl` is the + test's allowance for that, not a proved bound. + impact: MEDIUM in principle; bounded by the property test and the planted + `mutant_bh_no_cummin`. + closes-with: A verified floating-point model. Not planned. + +R-BH-7 + statement: | + The q-values control the false discovery rate at the nominal level. + status: OUT OF SCOPE + why: | + FDR control is a theorem about a probability measure over a sampling + model. This development is about exact rational arithmetic. Nothing in + `proofs/` bears on it, and `proofs/PROOF-STATUS.md` says so in its own + words. + impact: n/a (the claim is not made). + closes-with: n/a diff --git a/proofs/residue/clr.residue b/proofs/residue/clr.residue new file mode 100644 index 0000000..ae18a3b --- /dev/null +++ b/proofs/residue/clr.residue @@ -0,0 +1,96 @@ +# SPDX-License-Identifier: CC-BY-SA-4.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Residue — the centred log-ratio transform (P1) +# +# What `proofs/agda/CompositionalDA/CLR.agda` does NOT prove, each with a +# status and what would close it. A residue entry is a named gap, not a +# claim; `proofs/PROOF-STATUS.md` is the table of what IS proved. +# +# Format (one record per obligation): +# R-- +# statement: what would be proved +# status: OPEN | NOT ATTEMPTED | OUT OF SCOPE | CLOSED (by ) +# why: the obstacle, or the reason it is out of scope +# impact: what trusting the unproved claim would cost +# closes-with: the shape of a closing proof + +R-CLR-1 + statement: | + The CLR theorems hold for the real-valued transform + `clr(x)_i = log x_i − (1/n) Σ_j log x_j` on positive reals. + status: OUT OF SCOPE (by construction) + why: | + The Agda development is parametric: `R` is any commutative ring, `log` is + any monoid-to-ring homomorphism (`LogHom`), and `1/n` is any reciprocal + witness (`Reciprocal`). The real numbers, `exp`/`log` as functions, and + the `n · (1/n) = 1` identity over ℝ are not constructed here. The theorems + therefore say: ANY structure satisfying the two records has these three + properties. That the reals (with the natural logarithm) are such a + structure is a fact about ℝ, not about CLR, and is not formalised. + impact: | + LOW for the algebra, which is exactly what the parametric theorems cover. + The step from "the reals satisfy LogHom" to the Julia `clr` is covered by + `test/property_clr.jl`, which is tested, not proved. + closes-with: | + A `LogHom` instance over a constructive real-number ring with a proved + logarithm homomorphism. Not planned. + +R-CLR-2 + statement: | + The named assumptions are satisfiable: there is at least one non-trivial + `LogHom` and one `Reciprocal`. + status: CLOSED (by `CompositionalDA.CLR.Instance.additive-log` and + `CompositionalDA.CLR.Instance.reciprocal`) + why: | + A theorem whose hypotheses nothing satisfies is vacuous. The instance + module takes the monoid `(ℕ, +)` (`Data.Nat.Properties.+-0-commutativeMonoid`) + and the ring `ℚᵘ` (`+-*-commutativeRing`), with `log a = mkℚᵘ (+ a) 0`, + i.e. `a/1` — additive by definition, so the "logarithm" law is a + denominator calculation — and a reciprocal `1/(m+1)` for every `suc m`. + It type-checks `CLR` at `n = 2` against them (`module CLR₂`). + impact: n/a + closes-with: n/a + +R-CLR-3 + statement: | + CLR is NOT invariant under an arbitrary perturbation (only under a + constant one). + status: CLOSED (by `CompositionalDA.CLR.Instance.clr-not-shift-invariant`, + and exercised as the negative control `proofs/agda/reject/FakeCLR.agda`) + why: | + `clr-scale-invariant` needs the perturbing composition to be constant. A + proof that the hypothesis can be dropped must not exist; the instance + module proves its negation at `n = 2`, and the gate self-test checks that + Agda rejects the attempt. + impact: n/a + closes-with: n/a + +R-CLR-4 + statement: | + Floating-point `clr` in `src/` agrees with the exact transform to within + a stated tolerance. + status: OUT OF SCOPE + why: | + IEEE-754 semantics are not modelled. `test/property_clr.jl` checks the + three properties numerically with a tolerance, which is evidence about the + implementation and not a proof about it. + impact: | + MEDIUM in principle (a rounding-order bug would be invisible to the Agda + layer); bounded in practice by the property tests and their planted + mutants in `test/mutants.jl`. + closes-with: | + A verified floating-point model. Not planned. + +R-CLR-5 + statement: | + Anything about the Aitchison geometry of the simplex (distances, the ILR + basis, the isometry). + status: NOT ATTEMPTED + why: | + Not required by P1. The three P1 properties are the ones the downstream + methods rely on. + impact: LOW. + closes-with: | + An inner-product space over `R` and a proof that `clr` is an isometry + onto the zero-sum hyperplane. diff --git a/proofs/tests/gate-selftest.sh b/proofs/tests/gate-selftest.sh index 0661078..12a73d4 100755 --- a/proofs/tests/gate-selftest.sh +++ b/proofs/tests/gate-selftest.sh @@ -2,130 +2,272 @@ # SPDX-License-Identifier: MPL-2.0 # SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) # -# Show that the proof gate rejects what it claims to reject. +# Prove that the proof gate can fail. # -# A gate that has never been seen to fail is not evidence of anything. This -# script plants known defects and requires the gate to catch every one: +# A gate that has never been observed to reject anything is not evidence. This +# script takes a throwaway copy of the proof tree, breaks it in one specific way +# at a time, and requires the gate to reject it. If any mutation is *accepted*, +# or if a mutation failed to apply at all, this script exits non-zero. # -# A. Each negative control in proofs/agda/reject/ must FAIL to type-check, -# and Agda's output must match the control's own `-- EXPECT: ` -# line. A control that fails for some other reason (a missing import, a -# typo) does not count: it would keep "passing" after the defect it -# guards against had been reintroduced. -# B. Each of the seven checks in axiom-audit.sh must flag a defect planted -# in a scratch copy of the proof tree. The unmodified copy must pass -# first, so a planted failure cannot be an artefact of the copy. +# Three kinds of control: +# 1. planted defects in a theorem — the type-checker must reject them; +# 2. planted defects in the tree's hygiene (a module dropped from the entry +# point, a postulate, `--safe` removed) — the axiom audit must reject them; +# 3. the standing negative controls in `proofs/agda/reject/` — Agda must +# reject each one with the exact error its `-- EXPECT:` line names, so that +# a rejection for the wrong reason (a missing import, say) does not pass +# for the right one. +# And finally the pristine tree must still be accepted, so that a gate which +# rejects everything cannot pass this test either. # -# Prerequisite: proofs/bootstrap.sh has been run (it writes the libraries file -# and vendors the standard library). Every failure exits non-zero. +# It runs in CI. A gate whose self-test is skipped is a gate nobody can trust, +# so there is no flag to turn this off. +# +# Usage: proofs/tests/gate-selftest.sh +# +# Environment (same interface as proofs/bootstrap.sh): +# AGDA_BIN path to the agda binary (default: `agda` on PATH) +# PROOFS_VENDOR where bootstrap.sh vendored the stdlib (default: proofs/.vendor) -set -euo pipefail +set -uo pipefail PROOFS_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" -AGDA_DIR="$PROOFS_DIR/agda" -LIBRARIES_FILE="$PROOFS_DIR/.agda-libraries" -AGDA="${AGDA_BIN:-agda}" - -failures=0 - -# Record one finding on stderr and keep going, so a run reports every gap. -fail() { printf 'gate-selftest: FAIL: %s\n' "$*" >&2; failures=$((failures + 1)); } - -# Print a progress line on stdout. -note() { printf 'gate-selftest: %s\n' "$*"; } - -[[ -f "$LIBRARIES_FILE" ]] \ - || { printf 'gate-selftest: FATAL: %s missing; run proofs/bootstrap.sh first\n' "$LIBRARIES_FILE" >&2; exit 1; } -command -v "$AGDA" >/dev/null 2>&1 || [[ -x "$AGDA" ]] \ - || { printf 'gate-selftest: FATAL: agda not found (set AGDA_BIN)\n' >&2; exit 1; } - -SCRATCH="$(mktemp -d)" -trap 'rm -rf "$SCRATCH"' EXIT - -############################################################################## -# A. Negative controls must fail, for the stated reason -############################################################################## - -# Type-check one reject/ control and require a non-zero exit whose output -# matches the regex on the control's `-- EXPECT:` line. -check_reject() { - local file="$1" rel expect out rc - rel="${file#"$AGDA_DIR"/}" - expect="$(sed -n 's/^-- EXPECT: //p' "$file" | head -1)" - if [[ -z "$expect" ]]; then - fail "$rel has no '-- EXPECT:' line, so its failure reason cannot be checked" - return - fi - set +e - out="$(cd "$AGDA_DIR" && "$AGDA" --library-file="$LIBRARIES_FILE" --safe --without-K "$rel" 2>&1)" - rc=$? - set -e - if [[ $rc -eq 0 ]]; then - fail "$rel type-checked; the gate accepts a defect it must reject" - elif ! grep -Eq -- "$expect" <<<"$out"; then - fail "$rel failed, but not for the expected reason /$expect/: $(tail -3 <<<"$out" | tr '\n' ' ')" +NAMESPACE="CompositionalDA" +LIB_BASENAME="compositionalda.agda-lib" +VENDOR="${PROOFS_VENDOR:-$PROOFS_DIR/.vendor}" + +# Print a fatal error on stderr and exit 1. Every setup failure goes through +# here; a self-test that cannot run must not look like a self-test that passed. +die() { printf 'gate-selftest: FATAL: %s\n' "$*" >&2; exit 1; } + +# Set AGDA from AGDA_BIN or PATH, the same way bootstrap.sh does. There is no +# vendored-binary fallback: the estate installs Agda from the release tarball, +# and an absent prover is a failure here, not a skip. +resolve_agda() { + if [[ -n "${AGDA_BIN:-}" ]]; then + AGDA="$AGDA_BIN" + elif command -v agda >/dev/null 2>&1; then + AGDA="$(command -v agda)" else - note "rejected as expected: $rel" + die "agda not found (run proofs/bootstrap.sh, or set AGDA_BIN)" fi + [[ -x "$AGDA" ]] || die "$AGDA is not executable" } -mapfile -t REJECTS < <(find "$AGDA_DIR/reject" -name '*.agda' -type f | sort) -[[ ${#REJECTS[@]} -gt 0 ]] || fail "no negative controls under $AGDA_DIR/reject" -for f in "${REJECTS[@]}"; do check_reject "$f"; done - -############################################################################## -# B. Every axiom-audit check must fire on a planted defect -############################################################################## - -# Make a fresh scratch copy of the audited tree (agda/ and tests/) in $1. -fresh_copy() { - rm -rf "$1" - mkdir -p "$1" - cp -r "$AGDA_DIR" "$1/agda" - cp -r "$PROOFS_DIR/tests" "$1/tests" - rm -rf "$1/agda/_build" +# Set STDLIB_LIBS to the library files the gate depends on, read from the +# `proofs/.agda-libraries` that bootstrap.sh wrote (every line except the +# project's own library file). Falls back to the conventional vendored +# location if that file is absent, and fails if neither exists. +resolve_stdlib() { + STDLIB_LIBS=() + local libs_file="$PROOFS_DIR/.agda-libraries" line + if [[ -f "$libs_file" ]]; then + while IFS= read -r line; do + [[ -z "$line" ]] && continue + [[ "$(basename "$line")" == "$LIB_BASENAME" ]] && continue + [[ -f "$line" ]] || die "$libs_file names $line, which does not exist" + STDLIB_LIBS+=("$line") + done < "$libs_file" + elif [[ -f "$VENDOR/agda-stdlib/standard-library.agda-lib" ]]; then + STDLIB_LIBS+=("$VENDOR/agda-stdlib/standard-library.agda-lib") + fi + [[ ${#STDLIB_LIBS[@]} -gt 0 ]] \ + || die "standard library not found; run proofs/bootstrap.sh first" } -# Run the audit inside the scratch copy $1; prints its combined output and -# returns its exit status. -run_audit() { "$1/tests/axiom-audit.sh" 2>&1; } +resolve_agda +resolve_stdlib -CLEAN="$SCRATCH/clean" -fresh_copy "$CLEAN" -if ! out="$(run_audit "$CLEAN")"; then - fail "the unmodified copy fails the audit, so planted results would mean nothing: $out" -else - note "positive control: the unmodified copy passes the audit" -fi +WORK="$(mktemp -d)" +trap 'rm -rf "$WORK"' EXIT + +# Replace the working copy with a fresh copy of the pristine proof tree and +# the test scripts. Interface files (`_build/`, `*.agdai`) are removed from +# the copy so that nothing checked against the pristine sources can be reused +# for a mutated one. +reset_tree() { + rm -rf "$WORK/agda" "$WORK/tests" + cp -r "$PROOFS_DIR/agda" "$WORK/agda" + cp -r "$PROOFS_DIR/tests" "$WORK/tests" + rm -rf "$WORK/agda/_build" + find "$WORK/agda" -name '*.agdai' -delete + [[ -f "$WORK/agda/$LIB_BASENAME" ]] || die "missing $PROOFS_DIR/agda/$LIB_BASENAME" +} + +reset_tree + +# Point Agda at the *copy*: without this the library file would resolve +# `CompositionalDA.All` back to the pristine tree and every mutation below would +# be checked against unmutated sources — a self-test that tests nothing. +{ + printf '%s\n' "$WORK/agda/$LIB_BASENAME" + printf '%s\n' "${STDLIB_LIBS[@]}" +} > "$WORK/libraries" -TARGET="agda/CompositionalDA/Prelude.agda" - -# Plant one defect with the shell snippet $2 (run inside the copy), then -# require the audit to exit non-zero with a finding matching the regex $3. -plant() { - local name="$1" mutate="$2" expect="$3" dir="$SCRATCH/$1" out - fresh_copy "$dir" - (cd "$dir" && eval "$mutate") - if out="$(run_audit "$dir")"; then - fail "audit check '$name' missed its planted defect" - elif ! grep -Eq -- "$expect" <<<"$out"; then - fail "audit check '$name' failed, but without the expected finding /$expect/: $out" +# Type-check the copy's entry module with the gate's exact flags. Output goes +# to $WORK/out.txt; the exit status is the verdict. +run_gate() { + ( cd "$WORK/agda" && "$AGDA" --library-file="$WORK/libraries" --safe --without-K "$NAMESPACE/All.agda" ) \ + >"$WORK/out.txt" 2>&1 +} + +# Run the copied axiom audit, which resolves its own tree from its location +# and so audits the copy. +run_audit() { "$WORK/tests/axiom-audit.sh" >"$WORK/out.txt" 2>&1; } + +total=0 +failed=0 + +# expect_reject [reason-regex] +# +# Restore (relative to proofs/agda) in the copy, apply the mutation +# (a shell snippet receiving the file as $1), and require to fail. +# Counts a FAIL if the mutation did not change the file (a stale pattern +# would otherwise pass silently as "rejected") or if the runner accepted the +# mutated tree. With a [reason-regex], the runner must also fail WITH that +# finding: a rejection for some other reason is a FAIL. Restores the file +# afterwards. +expect_reject() { + local name="$1" runner="$2" target="$3" mutator="$4" reason="${5:-}" + local file="$WORK/agda/$target" + total=$((total + 1)) + + cp "$PROOFS_DIR/agda/$target" "$file" + local before after + before="$(sha256sum "$file" | awk '{print $1}')" + bash -c "$mutator" mutator "$file" + after="$(sha256sum "$file" | awk '{print $1}')" + if [[ "$before" == "$after" ]]; then + printf 'gate-selftest: FAIL %-48s mutation did not apply (stale pattern?)\n' "$name" + failed=$((failed + 1)) + cp "$PROOFS_DIR/agda/$target" "$file" + return + fi + + if $runner; then + printf 'gate-selftest: FAIL %-48s gate ACCEPTED a broken proof\n' "$name" + sed 's/^/ | /' "$WORK/out.txt" | head -6 + failed=$((failed + 1)) + elif [[ -n "$reason" ]] && ! grep -qE -- "$reason" "$WORK/out.txt"; then + printf 'gate-selftest: FAIL %-48s rejected, but not for the expected reason\n' "$name" + printf ' | expected /%s/\n' "$reason" + sed 's/^/ | /' "$WORK/out.txt" | head -6 + failed=$((failed + 1)) else - note "audit caught planted defect: $name" + printf 'gate-selftest: ok %-48s rejected\n' "$name" fi + cp "$PROOFS_DIR/agda/$target" "$file" } -plant safe "sed -i 's/ --safe//' $TARGET" 'does not declare --safe' -plant postulate "printf '\npostulate\n bogus : Set\n' >> $TARGET" 'postulate block' -plant ffi "printf '\n{-# FOREIGN GHC import Data.List #-}\n' >> $TARGET" 'FOREIGN/COMPILE/BUILTIN' -plant unsound "sed -i 's/--safe/--safe --type-in-type/' $TARGET" 'unsound flag' -plant trust "printf '\n-- uses trustMe\n' >> $TARGET" 'trustMe/primTrust' -plant hole "printf '\nhole = {! !}\n' >> $TARGET" 'contains a hole' -plant reach "printf '{-# OPTIONS --without-K --safe #-}\nmodule CompositionalDA.Orphan where\n' > agda/CompositionalDA/Orphan.agda" \ - 'not reachable' - -if [[ $failures -gt 0 ]]; then - printf 'gate-selftest: %d gap(s) in the gate\n' "$failures" >&2 - exit 1 +# expect_reject_for_reason +# +# Type-check one standing negative control from proofs/agda/reject/ and +# require BOTH that Agda exits 42 (a type error, as opposed to 1 for a +# missing file or library) AND that its output matches the extended regex on +# the control's `-- EXPECT:` line. A control with no EXPECT line is a FAIL: +# an unexplained rejection is not a control. +expect_reject_for_reason() { + local target="$1" name file pattern rc + name="reject/$(basename "${target%.agda}")" + file="$WORK/agda/$target" + total=$((total + 1)) + + pattern="$(sed -nE 's/^-- EXPECT: (.*)$/\1/p' "$file" | head -1)" + if [[ -z "$pattern" ]]; then + printf 'gate-selftest: FAIL %-48s has no -- EXPECT: line\n' "$name" + failed=$((failed + 1)) + return + fi + + ( cd "$WORK/agda" && "$AGDA" --library-file="$WORK/libraries" --safe --without-K "$target" ) \ + >"$WORK/out.txt" 2>&1 + rc=$? + if [[ $rc -eq 0 ]]; then + printf 'gate-selftest: FAIL %-48s gate ACCEPTED a negative control\n' "$name" + failed=$((failed + 1)) + elif [[ $rc -ne 42 ]]; then + printf 'gate-selftest: FAIL %-48s exited %d, not 42 (not a type error)\n' "$name" "$rc" + sed 's/^/ | /' "$WORK/out.txt" | head -6 + failed=$((failed + 1)) + elif ! grep -qE -- "$pattern" "$WORK/out.txt"; then + printf 'gate-selftest: FAIL %-48s rejected, but not for the expected reason\n' "$name" + printf ' | expected /%s/\n' "$pattern" + sed 's/^/ | /' "$WORK/out.txt" | grep -v '^ *| *Checking' | head -6 + failed=$((failed + 1)) + else + printf 'gate-selftest: ok %-48s rejected for the expected reason\n' "$name" + fi +} + +# --- the type-checker must reject wrong mathematics ------------------------- +# +# Each mutator is a shell snippet receiving the target file as $1. If a snippet +# stops matching (because the proof was rewritten), `expect_reject` reports the +# mutation as unapplied rather than silently passing. + +expect_reject "clr: sum-to-zero claimed as sum-to-one" run_gate "$NAMESPACE/CLR.agda" \ + 'sed -i "s|clr-sum-zero : ∀ (x : Vector Pos n) → sum (clr x) ≈ 0#|clr-sum-zero : ∀ (x : Vector Pos n) → sum (clr x) ≈ 1#|" "$1"' + +expect_reject "clr: centring removed from the transform" run_gate "$NAMESPACE/CLR.agda" \ + 'sed -i "s|^clr x = centre (map log x)|clr x = map log x|" "$1"' + +expect_reject "bh: family-size monotonicity reversed" run_gate "$NAMESPACE/BenjaminiHochberg.agda" \ + 'sed -i "s| bhScale M j n d ≤ℚ bhScale M′ j n d| bhScale M′ j n d ≤ℚ bhScale M j n d|" "$1"' + +# Mirrors test/mutants.jl `mutant_bh_no_cummin`: scale each rank, never take +# the running minimum. +expect_reject "stepup: running minimum dropped" run_gate "$NAMESPACE/BenjaminiHochberg/StepUp.agda" \ + 'sed -i "s| (run ⊓ scale m k p) ∷ stepUp m (run ⊓ scale m k p) ps| scale m k p ∷ stepUp m run ps|" "$1"' + +expect_reject "stepup: q ≥ p reversed" run_gate "$NAMESPACE/BenjaminiHochberg/StepUp.agda" \ + 'sed -i "s| ∀ (i : Fin m) → lookup ps i ≤ lookup (adjust ps) i| ∀ (i : Fin m) → lookup (adjust ps) i ≤ lookup ps i|" "$1"' + +expect_reject "stepup: monotonicity reversed" run_gate "$NAMESPACE/BenjaminiHochberg/StepUp.agda" \ + 'sed -i "s| lookup (adjust ps) i ≤ lookup (adjust ps) j| lookup (adjust ps) j ≤ lookup (adjust ps) i|" "$1"' + +expect_reject "stepup: scale written as rank/m" run_gate "$NAMESPACE/BenjaminiHochberg/StepUp.agda" \ + 'sed -i "s|^scale m k p = p \* mkℚᵘ (+ m) k|scale m k p = p * mkℚᵘ (+ suc k) m|" "$1"' + +# --- the axiom audit must reject an unchecked or unsound module ------------- + +expect_reject "audit: module dropped from the gate entry" run_audit "$NAMESPACE/All.agda" \ + 'sed -i "/import CompositionalDA.BenjaminiHochberg.StepUp/d" "$1"' 'not reachable' + +expect_reject "audit: postulate injected" run_audit "$NAMESPACE/BenjaminiHochberg/StepUp.agda" \ + 'printf "postulate cheat : ∀ {A : Set} → A\n" >> "$1"' 'postulate block' + +expect_reject "audit: --safe removed" run_audit "$NAMESPACE/CLR.agda" \ + 'sed -i "s|--without-K --safe|--without-K|" "$1"' 'does not declare --safe' + +expect_reject "audit: FFI pragma injected" run_audit "$NAMESPACE/Prelude.agda" \ + 'printf "\n{-# FOREIGN GHC import Data.List #-}\n" >> "$1"' 'FOREIGN/COMPILE/BUILTIN' + +expect_reject "audit: unsound flag added" run_audit "$NAMESPACE/Prelude.agda" \ + 'sed -i "s|--safe|--safe --type-in-type|" "$1"' 'unsound flag' + +expect_reject "audit: trustMe mentioned" run_audit "$NAMESPACE/Prelude.agda" \ + 'printf "\n-- uses trustMe\n" >> "$1"' 'trustMe/primTrust' + +expect_reject "audit: hole left in a module" run_audit "$NAMESPACE/Prelude.agda" \ + 'printf "\nhole = {! !}\n" >> "$1"' 'contains a hole' + +# --- the standing negative controls must fail for their stated reason ------ + +while IFS= read -r control; do + expect_reject_for_reason "reject/$(basename "$control")" +done < <(find "$PROOFS_DIR/agda/reject" -maxdepth 1 -name '*.agda' -type f | sort) + +# --- the gate must still accept the pristine tree -------------------------- + +total=$((total + 1)) +reset_tree +if run_gate && run_audit; then + printf 'gate-selftest: ok %-48s accepted\n' "pristine tree" +else + printf 'gate-selftest: FAIL %-48s gate REJECTED the real proofs\n' "pristine tree" + sed 's/^/ | /' "$WORK/out.txt" | head -12 + failed=$((failed + 1)) fi -note "OK: ${#REJECTS[@]} negative control(s) rejected, 7 audit checks fire" + +printf 'gate-selftest: %d/%d controls behaved correctly\n' "$((total - failed))" "$total" +[[ $failed -eq 0 ]] || exit 1