diff --git a/README.md b/README.md index ac95ba2..f3f6593 100644 --- a/README.md +++ b/README.md @@ -4,14 +4,16 @@ # CompositionalDA.jl Native-Julia compositional differential abundance for microbiome count tables. -Methods: ALDEx2, ANCOM-BC, and a Songbird-style multinomial regression. +Methods: ALDEx2 and ANCOM-BC. Multinomial and Dirichlet-multinomial regression +live in a separate package, +[`CompositionalCounts.jl`](https://github.com/hyperpolymath/CompositionalCounts.jl). It needs no Python and no TensorFlow. R is needed only for the optional reference backend. > **Status: v0.1.0, core layer only.** Closure, the centred log-ratio, zero > policies, Benjamini–Hochberg and the result/provenance types are implemented. -> The methods land one per release: ALDEx2 (v0.2), ANCOM-BC (v0.3), multinomial -> (v0.4). Nothing here claims a method that is not yet in `src/`. +> The methods land one per release: ALDEx2 (v0.2), then ANCOM-BC (v0.3). +> Nothing here claims a method that is not yet in `src/`. ## Input contract @@ -30,17 +32,20 @@ q = bh_adjust([0.01, 0.04, 0.2]) # equals R's p.adjust(, "BH") ## Proved and tested are kept apart `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 theorem is +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`). Statistical properties such as calibration, -consistency and FDR control are **not** proved. They are tested against the R -reference implementations. The full table is in +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, +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). ## The R bridge is temporary by design -The optional reference backend calls Bioconductor ALDEx2/ANCOMBC through RCall. -It ships as a Julia package extension and exists for two purposes: to serve as +**Planned, not yet in the tree.** The optional reference backend will call +Bioconductor ALDEx2/ANCOMBC through RCall, shipped as a Julia package extension. +It will exist for two purposes: to serve as the test oracle, and to help researchers move off R. It has a written sunset criterion; see [`docs/R-SUNSET.md`](docs/R-SUNSET.md). diff --git a/docs/R-SUNSET.md b/docs/R-SUNSET.md index e7961b5..99e308c 100644 --- a/docs/R-SUNSET.md +++ b/docs/R-SUNSET.md @@ -3,6 +3,9 @@ # Sunset criterion for the R reference backend +> **Status:** the extension described here is planned and not yet in the tree: +> there is no `ext/` directory and no `test/oracle_*.jl` yet. + The R backend (`ext/CompositionalDARCallExt.jl`, loaded only when RCall is present) is a bridge, not a product. It never owns R state. Callers pass an `evaluator` function, so a host application keeps its own R session and lock. diff --git a/proofs/PROOF-STATUS.md b/proofs/PROOF-STATUS.md new file mode 100644 index 0000000..9c4766e --- /dev/null +++ b/proofs/PROOF-STATUS.md @@ -0,0 +1,62 @@ + + + +# 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. diff --git a/proofs/agda/CompositionalDA/BenjaminiHochberg.agda b/proofs/agda/CompositionalDA/BenjaminiHochberg.agda index 5b8159b..1b38d8e 100644 --- a/proofs/agda/CompositionalDA/BenjaminiHochberg.agda +++ b/proofs/agda/CompositionalDA/BenjaminiHochberg.agda @@ -79,7 +79,7 @@ monotone-in-family-size M M′ j n d M≤M′ {{np}} = (ℤₚ.≤-trans (ℤₚ.*-monoˡ-≤-nonNeg n {{np}} (ℤ.+≤+ M≤M′)) (ℤₚ.≤-reflexive (ℤₚ.*-comm n (+ M′))))) -- Non-negativity of the scaled value is left as an open obligation, recorded in --- `proofs/residue/benjamini-hochberg.residue`. The obstacle is mechanical, not +-- `proofs/PROOF-STATUS.md` (Open obligations). The obstacle is mechanical, not -- conceptual: `+ M * + n` does not reduce to `+ (M ℕ.* n)`, so the non-negativity -- lemmas in `Data.Integer.Properties`, which are stated about the reduced form, -- do not apply directly, and the sign/absolute-value form `ℤ._*_` actually diff --git a/proofs/agda/CompositionalDA/CLR.agda b/proofs/agda/CompositionalDA/CLR.agda index 4b2ad80..65d62bb 100644 --- a/proofs/agda/CompositionalDA/CLR.agda +++ b/proofs/agda/CompositionalDA/CLR.agda @@ -22,7 +22,7 @@ -- `clr-constant` (a constant composition has zero CLR) that (b) is built from. -- -- What is NOT here: anything about real numbers, `exp`, continuity, floating --- point, or the geometry of the simplex. See `proofs/residue/clr.residue`. +-- point, or the geometry of the simplex. See `proofs/PROOF-STATUS.md`. {-# OPTIONS --without-K --safe #-} diff --git a/proofs/bootstrap.sh b/proofs/bootstrap.sh index 8ed14df..7ecd28a 100755 --- a/proofs/bootstrap.sh +++ b/proofs/bootstrap.sh @@ -42,7 +42,7 @@ AGDA_RELEASE_URL="https://github.com/agda/agda/releases/tag/v${AGDA_VERSION}" # development line towards 3.0 — the branch self-identifies as # `standard-library-3.0` but no `v3.0` tag has been cut. They do NOT compile # against the latest release (v2.4) or against v2.1. See -# proofs/residue/toolchain.residue. +# proofs/PROOF-STATUS.md (Open obligations). STDLIB_VERSION="2ffa8b7d4e8e818717ad643d184f055a4d1b0447" REPO_ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" diff --git a/proofs/tests/axiom-audit.sh b/proofs/tests/axiom-audit.sh old mode 100644 new mode 100755 diff --git a/proofs/tests/gate-selftest.sh b/proofs/tests/gate-selftest.sh new file mode 100755 index 0000000..0661078 --- /dev/null +++ b/proofs/tests/gate-selftest.sh @@ -0,0 +1,131 @@ +#!/usr/bin/env bash +# 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. +# +# 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. 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. +# +# Prerequisite: proofs/bootstrap.sh has been run (it writes the libraries file +# and vendors the standard library). Every failure exits non-zero. + +set -euo 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' ' ')" + else + note "rejected as expected: $rel" + fi +} + +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" +} + +# 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; } + +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 + +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" + else + note "audit caught planted defect: $name" + fi +} + +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 +fi +note "OK: ${#REJECTS[@]} negative control(s) rejected, 7 audit checks fire"