Expose launch calibration blockers and add reproducible assurance experiments

.github/workflows/ci.yml

@@ -143,3 +143,16 @@ jobs:
           curl --fail --location --retry 3 https://raw.githubusercontent.com/hyperpolymath/ANCOMBC/6d094731ed85e9b0458053ad3510e9bdc238e152/R/ancombc_bias_correct.R -o /tmp/ancombc-variance-repair.R
           export R_HOME="$(R RHOME)"
           ANCOMBC_ESTIMATOR_SOURCE=/tmp/ancombc-variance-repair.R julia --project=test/oracle test/oracle/ancombc_oracle.jl
+      - name: Investigate the small-unbalanced calibration failure
+        if: matrix.reference == 'ANCOMBC'
+        run: |
+          julia --project=test/oracle benchmark/export_calibration_case.jl /tmp/calibration/small-unbalanced-inputs
+          Rscript benchmark/calibration_reference.R /tmp/calibration/small-unbalanced-inputs /tmp/ancombc-variance-repair.R
+      - name: Preserve calibration investigation evidence
+        if: always() && matrix.reference == 'ANCOMBC'
+        uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # v4.6.2
+        with:
+          name: small-unbalanced-calibration
+          path: /tmp/calibration
+          if-no-files-found: error
+          retention-days: 30

.gitignore

@@ -2,3 +2,4 @@
 *.agdai
 _build/
 /test/oracle/Manifest.toml
+/benchmark/Manifest.toml

EXPLAINME.adoc

@@ -0,0 +1,69 @@
+// SPDX-License-Identifier: CC-BY-SA-4.0
+// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath)
+= What CompositionalDA can establish
+
+This package analyses raw integer counts in samples-by-features orientation.
+Its current methods are unpaired two-group ALDEx2 and original ANCOM-BC.
+It does not implement ANCOM-BC2, covariate adjustment or repeated measures.
+Zero handling, the contrast direction, seeds and numerical policies belong in
+the analysis record. Swapping groups changes the effect direction.
+
+== Evidence and its limits
+
+[cols="2,3,3",options="header"]
+|===
+|Claim |Owning evidence |What remains outside it
+|Exact algebraic identities |`proofs/PROOF-STATUS.md`, safe Agda modules and negative controls
+|Floating-point refinement, sorter correctness, FDR and biological truth
+|Numerical method agreement |`test/oracle/`, pinned R/Bioconductor releases and explicitly repaired reference functions
+|Equivalence to pristine defective references, every input or later upstream releases
+|Regression resistance |`test/`, seeded fuzz/property tests and planted mutants
+|Exhaustive path/state coverage or the absence of unknown defects
+|Real processed microbiome inputs |`docs/DATASETS.md`, archived source hashes and gut/soil fixtures
+|Known ground truth, prospective clinical validity or general calibration
+|Exploratory operating characteristics |`docs/CALIBRATION-PILOT-PROTOCOL.adoc`, `benchmark/calibration_pilot.jl`
+|A calibration certificate; the small-unbalanced ANCOM-BC null pilot is adverse
+|Performance |`benchmark/benchmarks.jl`, machine-specific raw measurements
+|A speed advantage over R or an agreed latency/service-level guarantee
+|===
+
+The 2026-11-07 calibration pilot fixed six regimes and 32 independently seeded
+datasets per regime before measurement. Of 384 method fits, 365 returned and
+19 refused. Rates conditional on returned fits must be accompanied by refusal
+counts. Under the small-unbalanced complete null, original ANCOM-BC made false
+discoveries in 19/29 successful fits at 0.05 (Wilson 95% interval 0.473–0.801).
+The unconditional observed count is 19/32; refusals are not successful negative
+findings. The follow-up repaired-R experiment converged for all 32 datasets and
+also made false discoveries in 19. Of 29 paired successful native fits, 28 had
+maximum q-value differences below 8.5e-6, but replicate 27 differed by 0.5041
+and changed the discovery decision. Numerical conformance is therefore not
+established across this expanded regime. A reviewed support restriction or
+estimator correction, and fresh prespecified validation, are needed before a
+calibration claim. The raw receipt is CI run 42411431729, artifact
+`small-unbalanced-calibration`, at source `cbb8ac4d85b4259c8122ed83c2cce8ec6c34362a`.
+Do not fix it by tuning the pilot seeds, removing adverse regimes or quietly
+substituting another method.
+
+ALDEx2's pilot is much less adverse, but 32 datasets cannot certify 5% FDR.
+All-null FDR is the probability of any rejection. With signal present,
+family-wise error is the probability of at least one false rejection; mean
+false discovery proportion estimates FDR. The reported Wilson intervals are
+for event probabilities, not mean FDP, taxa treated as independent
+replications, or arbitrary populations.
+
+At only two ALDEx2 Monte-Carlo instances, a separate 4,096-seed probe produced
+58 infinite and three undefined effect estimates. Its p/q bounds still held.
+The WebUI refuses non-finite effects; this is an estimator usability limitation
+to address explicitly, not permission to coerce undefined effects to zero.
+
+== Trust and release decisions
+
+Julia, GMP/MPFR, BLAS, dependencies, the operating system, reference software,
+data provenance and proof checkers remain in the trusted environment.
+Agda algebra does not independently verify those implementations. No second
+Rocq/Lean witness or full executable refinement is claimed.
+
+General registration is explicitly on hold pending the human maintainer's review
+of the changed implementation. A signature authenticates a commit, not its
+mathematical or scientific correctness. Version 0.3.0 in Project.toml is not a
+stable-release grade. No CRG/TRG/RSR certificate is asserted.

README.md

@@ -17,6 +17,15 @@ backend.
 > gate now checks all statistics and shared Monte-Carlo instances. ANCOM-BC's
 > precise model and limitations are described in `docs/ANCOMBC-SPEC.md`.
 
+> **Launch review remains pending.** The predeclared calibration pilot found
+> severe false-discovery inflation for original ANCOM-BC with four versus
+> sixteen samples under a complete null: 19 of 29 successful fits made at least
+> one discovery at q <= 0.05; three additional fits refused to return a result.
+> This is a release blocker for a general calibration claim. See
+> [`docs/LAUNCH-READINESS.adoc`](docs/LAUNCH-READINESS.adoc), the reproducible
+> protocol and scripts in `benchmark/`, and [`EXPLAINME.adoc`](EXPLAINME.adoc).
+> Passing numerical reference comparisons does not establish FDR control.
+
 The v0.3 implementation, tests, proofs and integration work were developed
 with OpenAI Codex assistance. General registration requires a human maintainer
 to review and understand these contributions; automated checks and code review

SECURITY.md

@@ -0,0 +1,25 @@
+<!-- SPDX-License-Identifier: CC-BY-SA-4.0 -->
+<!-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) -->
+
+# Security and numerical integrity
+
+This is a development package; no released version currently carries a
+security support commitment. Report exploitable defects through GitHub's
+private vulnerability-reporting facility if enabled. Otherwise contact the
+maintainer at the address in Project.toml before publishing exploit details.
+Ordinary numerical defects can be reported as issues with a minimal synthetic
+reproducer, package/tool versions and expected versus observed results.
+Do not attach confidential or identifiable biological data.
+
+The native API is a local scientific library, not a sandbox for untrusted jobs.
+Callers must bound matrix dimensions, count bit widths, Monte-Carlo work,
+concurrency and memory at their service boundary. Input validation and explicit
+nonconvergence errors do not create a hard CPU or memory limit. R reference
+execution loads installed R packages and must run in a separately controlled
+environment when inputs or dependencies are untrusted.
+
+Aqua checks Julia package quality. It is not a vulnerability-database scanner.
+The launch ledger must separately record dependency advisories, secret scanning,
+workflow permissions/action pins, release provenance and unresolved findings.
+No successful comprehensive security audit is claimed here. Statistics that
+fail calibration are safety/integrity concerns even without an attacker.

benchmark/Project.toml

@@ -0,0 +1,11 @@
+# SPDX-License-Identifier: MPL-2.0
+# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath)
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CompositionalDA = "bb1787b7-c25e-80cd-8559-27e3f8d35816"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+
+[compat]
+BenchmarkTools = "1"
+Distributions = "0.25"
+julia = "1.10"

benchmark/README.md

@@ -0,0 +1,52 @@
+<!-- SPDX-License-Identifier: CC-BY-SA-4.0 -->
+<!-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) -->
+
+# Exploratory launch measurements
+
+From the repository root, prepare the separate environment:
+
+```sh
+julia --project=benchmark -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+julia --project=benchmark benchmark/calibration_pilot.jl /path/to/pilot-output
+julia --project=benchmark benchmark/benchmarks.jl /path/to/benchmark-output
+```
+
+Read `docs/CALIBRATION-PILOT-PROTOCOL.adoc` before interpreting the output.
+Seeds, scenarios and comparisons were fixed before results. Refusals are
+separate from successful fits. Neither script certifies calibration, speed,
+<!-- coderabbit-ai-generated:v1 -->
+
+<!-- coderabbit-ai-generated-section:v1:2daa59cb87a49803898dd374077604af362283cb0f06dd91a39c6483729ff4fc -->
+
+stability or package readiness. Keep generated measurements outside source.
+Each output directory includes the active environment's `Project.toml` and
+resolved `Manifest.toml`, including transitive dependencies such as StatsFuns.
+Retain both with the results; instantiate the active project before running.
+Timing benchmarks refuse a known dirty source tree (including untracked files).
+When Git provenance is unavailable, scripts warn and record `source_head` and
+`source_tree` as `unknown`; those runs need a separate source archive to reproduce.
+
+<!-- /coderabbit-ai-generated:v1 -->
+Benchmark timings exclude compilation/warmup and retain trial counts; ten
+trials provide an exploratory baseline, not precise tail-latency estimates.
+
+The optional CI investigation exports the small-unbalanced complete-null
+inputs and native results, then runs pinned ANCOMBC 2.12.0 with the explicitly
+identified variance-estimator repair. It records R convergence rather than
+treating an iteration-cap result as a valid fit. A successful investigation
+job means the experiment completed, not that the adverse calibration passed.
+<!-- coderabbit-ai-generated:v1 -->
+
+<!-- coderabbit-ai-generated-section:v1:293a501349b835ccb7aa6a0ef06715bbeee5f83376f52e981bd78fa617317411 -->
+
+<!-- coderabbit-ai-generated:v1 -->
+
+<!-- coderabbit-ai-generated-section:v1:782763b315c6b798efa73f3b136076ffe1bc009ef35e82ac28614d5de8ad0002 -->
+
+The workflow requests 30 days of artifact retention, but repository policy
+may shorten it. Check the artifact expiry and archive it with its SHA-256
+before using it in a release decision.
+
+<!-- /coderabbit-ai-generated:v1 -->
+
+<!-- /coderabbit-ai-generated:v1 -->

benchmark/benchmarks.jl

@@ -0,0 +1,43 @@
+# SPDX-License-Identifier: MPL-2.0
+# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath)
+using BenchmarkTools, CompositionalDA, Statistics, TOML, Dates
+include("calibration_pilot.jl")
+destination=only(ARGS)
+mkpath(destination)
+provenance=source_provenance(pkgdir(CompositionalDA);require_clean=true)
+save_run_environment(destination)
+rows=NamedTuple[]
+for (n,d) in ((20,40),(20,400),(100,400))
+    regime=(name="benchmark",n1=n÷2,n2=n÷2,depth1=10000,depth2=10000,sd=.35,sparse=false,signals=4)
+    x,f,g,_=simulate(regime,20261107+n+d;d)
+    for method in (:closure,:clr,:bh,:aldex2,:ancombc)
+        p=collect(range(.001,.999;length=d))
+        call = method==:closure ? (() -> closure(x)) :
+            method==:clr ? (() -> clr(Pseudocount(1),x)) :
+            method==:bh ? (() -> bh_adjust(p)) :
+            method==:aldex2 ? (() -> aldex2(x,f,g;mc_samples=128,seed=123)) :
+            (() -> ancombc(x,f,g;structural_zeros=false))
+        try
+            call() # compilation/warmup is excluded and labelled as such
+            trial=@benchmark $call() samples=10 evals=1 seconds=15
+            push!(rows,(method=String(method),samples=n,features=d,status="measured",
+                trials=length(trial.times),median_ns=median(trial.times),p95_ns=quantile(trial.times,.95),
+                bytes=trial.memory,allocations=trial.allocs))
+        catch error
+            error isa CompositionError || rethrow()
+            push!(rows,(method=String(method),samples=n,features=d,status="refused",
+                trials=0,median_ns=NaN,p95_ns=NaN,bytes=0,allocations=0))
+            println("REFUSED ",method," ",n,"x",d,": ",sprint(showerror,error))
+        end
+    end
+end
+open(joinpath(destination,"benchmarks.csv"),"w") do io
+    println(io,join(keys(first(rows)),','))
+    foreach(row->println(io,join(values(row),',')),rows)
+end
+metadata=Dict("julia"=>string(VERSION),"threads"=>Threads.nthreads(),"cpu"=>Sys.CPU_NAME,
+    "machine"=>Sys.MACHINE,"word_size"=>Sys.WORD_SIZE,"date_utc"=>string(now(UTC)),
+    "benchmarktools"=>string(pkgversion(BenchmarkTools)),"warmup_excluded"=>true,
+    "claim"=>"exploratory warm-process baseline; no reference speedup or acceptance budget")
+merge!(metadata,provenance)
+open(io->TOML.print(io,metadata),joinpath(destination,"benchmarks-metadata.toml"),"w")

benchmark/calibration_pilot.jl

@@ -0,0 +1,119 @@
+# SPDX-License-Identifier: MPL-2.0
+# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+# Run with the package project; results are written only to the supplied directory.
+using CompositionalDA, Random, Statistics, Dates, Printf, SHA, TOML
+using CompositionalDA: StableRNG
+using Distributions: Multinomial
+include("run_environment.jl")
+
+const REGIMES = [
+    (name="balanced_null",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=0),
+    (name="depth_imbalance_null",n1=10,n2=10,depth1=1000,depth2=10000,sd=0.35,sparse=false,signals=0),
+    (name="small_unbalanced_null",n1=4,n2=16,depth1=1000,depth2=1000,sd=0.35,sparse=false,signals=0),
+    (name="sparse_null",n1=6,n2=6,depth1=200,depth2=200,sd=0.7,sparse=true,signals=0),
+    (name="sparse_signal",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=4),
+    (name="dense_balanced_signal",n1=10,n2=10,depth1=10000,depth2=10000,sd=0.35,sparse=false,signals=20),
+]
+
+function simulate(regime,seed;d=40)
+    rng=StableRNG(seed)
+    baseline=regime.sparse ? exp.(range(-4,0;length=d)) : collect(range(1000,5000;length=d))
+    logeffect=zeros(d)
+    half=regime.signals÷2
+    logeffect[1:half].=log(4)
+    logeffect[half+1:2half].=-log(4)
+    x=Matrix{Int}(undef,regime.n1+regime.n2,d)
+    for i in axes(x,1)
+        condition=i>regime.n1
+        weights=baseline .* exp.(regime.sd .* randn(rng,d) .+ condition .* logeffect)
+        x[i,:]=rand(rng,Multinomial(condition ? regime.depth2 : regime.depth1,weights/sum(weights)))
+    end
+    return x,string.(1:d),vcat(fill(0,regime.n1),fill(1,regime.n2)),logeffect
+end
+
+function wilson(events,total)
+    total==0 && return (NaN,NaN)
+    z=1.959963984540054; p=events/total; scale=1+z^2/total
+    centre=(p+z^2/(2total))/scale
+    radius=z*sqrt(p*(1-p)/total+z^2/(4total^2))/scale
+    return (max(0,centre-radius),min(1,centre+radius))
+end
+
+function run_pilot(destination)
+    mkpath(destination)
+    provenance=source_provenance(pkgdir(CompositionalDA))
+    save_run_environment(destination)
+    protocol=joinpath(@__DIR__,"..","docs","CALIBRATION-PILOT-PROTOCOL.adoc")
+    metadata=Dict("protocol_sha256"=>bytes2hex(sha256(read(protocol))),
+        "julia"=>string(VERSION),"package_version"=>string(pkgversion(CompositionalDA)),
+        "started_utc"=>string(now(UTC)),"replicates_per_regime"=>32,
+        "mc_samples"=>128,"seed_formula"=>"20261107 + 10000*regime_index + replicate",
+        "status"=>"exploratory pilot; not a calibration certificate")
+    merge!(metadata,provenance)
+    records=NamedTuple[]
+    failures=NamedTuple[]
+    for (regime_index,regime) in enumerate(REGIMES), replicate in 1:32
+        seed=20261107+10000regime_index+replicate
+        x,features,groups,truth=simulate(regime,seed)
+        for method in (:aldex2,:ancombc)
+            start=time_ns()
+            fitted=try
+                method==:aldex2 ? aldex2(x,features,groups;mc_samples=128,seed=seed+1_000_000) :
+                    ancombc(x,features,groups;prevalence=0.1,structural_zeros=false,conserve=true)
+            catch error
+                if error isa CompositionError
+                    push!(failures,(regime=regime.name,replicate,method=string(method),seed,
+                        reason=sprint(showerror,error),seconds=(time_ns()-start)/1e9))
+                    continue
+                end
+                rethrow()
+            end
+            elapsed=(time_ns()-start)/1e9
+            table=fitted.table
+            @assert all(0 .<= table.pvalue .<= table.qvalue .<= 1)
+            ids=parse.(Int,table.feature)
+            for alpha in (0.05,0.10)
+                discovered=table.qvalue .<= alpha
+                false_calls=count(discovered .& iszero.(truth[ids]))
+                true_calls=count(discovered .& .!iszero.(truth[ids]))
+                push!(records,(regime=regime.name,replicate,method=string(method),seed,alpha,
+                    tested=length(ids),excluded=40-length(ids),false_calls,true_calls,
+                    fdp=false_calls/max(false_calls+true_calls,1),
+                    power=regime.signals==0 ? NaN : true_calls/regime.signals,seconds=elapsed))
+            end
+        end
+        replicate==32 && println("Completed ",regime.name)
+        flush(stdout)
+    end
+    open(joinpath(destination,"calibration-pilot.csv"),"w") do io
+        println(io,join(string.(keys(first(records))),','))
+        for record in records
+            println(io,join(values(record),','))
+        end
+    end
+    metadata["finished_utc"]=string(now(UTC))
+    metadata["successful_fits"]=length(records)÷2
+    metadata["refused_fits"]=length(failures)
+    metadata["refusals"]=[Dict(string(k)=>v for (k,v) in pairs(f)) for f in failures]
+    open(io->TOML.print(io,metadata),joinpath(destination,"calibration-pilot-metadata.toml"),"w")
+    open(joinpath(destination,"calibration-pilot-summary.md"),"w") do io
+        println(io,"# Exploratory calibration pilot\n\n32 independent datasets per regime. Metrics condition on successful fits; refusals are separate. No calibration certification.\n")
+        println(io,"| Regime | Method | q threshold | Fits/32 | Mean FDP | Mean power | Any false discovery: fraction [95% Wilson interval] |")
+        println(io,"| --- | --- | --- | --- | --- | --- | --- |")
+        for regime in REGIMES,method in ("aldex2","ancombc"),alpha in (0.05,0.10)
+            rows=filter(r->r.regime==regime.name && r.method==method && r.alpha==alpha,records)
+            isempty(rows) && (println(io,"| $(regime.name) | $method | $alpha | 0/32 | — | — | — |");continue)
+            events=count(r->r.false_calls>0,rows)
+            lo,hi=wilson(events,length(rows))
+            power=regime.signals==0 ? "N/A" : @sprintf("%.3f",mean(r.power for r in rows))
+            @printf(io,"| %s | %s | %.2f | %d/32 | %.3f | %s | %.3f [%.3f, %.3f] |\n",
+                regime.name,method,alpha,length(rows),mean(r.fdp for r in rows),power,events/length(rows),lo,hi)
+        end
+        println(io,"\nFor complete-null regimes, any false discovery equals the false-discovery-rate event. For signal regimes that event is FWER, not mean FDP. The interval covers the event probability among successful fits, not arbitrary populations or the mean FDP. See the fixed protocol and refusal records.")
+    end
+    println("Successful fits: ",metadata["successful_fits"],"; refusals: ",length(failures))
+end
+
+if abspath(PROGRAM_FILE)==(@__FILE__)
+    run_pilot(isempty(ARGS) ? error("Supply an output directory") : only(ARGS))
+end