gam_terms/inference/structure_evidence.rs
1//! Anytime-valid structure discovery: e-process gates, universal-inference
2//! atom tests, e-BH error control, and KL-optimal steering probes.
3//!
4//! Interpretability as sequential experimental design.
5//!
6//! # The thesis
7//!
8//! The dictionary-learning stack (#974–#981) DISCOVERS structure: atom
9//! birth/death/fission/fusion (#976), geometry adjudication — circle vs
10//! clusters vs line (#907), feature binding (#975). Today those decisions
11//! are made by evidence heuristics (likelihood-ratio ladders, BIC-flavored
12//! gates). Three facts, never previously combined, say that is not merely
13//! informal but WRONG in a specific, fixable way — and that fixing it
14//! upgrades the whole capstone from observational description to
15//! error-controlled experimental science:
16//!
17//! 1. **Atom existence is a NON-REGULAR testing problem.** "Does a K+1-th
18//! dictionary atom exist?" is testing K vs K+1 mixture components — the
19//! textbook boundary/loss-of-identifiability case where the classical
20//! likelihood-ratio χ² asymptotics FAIL (the null sits on the boundary
21//! of the alternative; the nuisance parameters of the new atom vanish
22//! under the null — Davies' problem). Every SAE/dictionary paper that
23//! thresholds a likelihood improvement is running this broken test.
24//! **Universal inference** (Wasserman–Ramdas–Balakrishnan 2020) is the
25//! modern resolution: a split-likelihood-ratio e-value that is valid in
26//! finite samples with NO regularity conditions whatsoever — exactly
27//! the irregular regime atom birth lives in.
28//!
29//! 2. **Discovery happens on streams with optional stopping.** Dictionaries
30//! are trained until the features "look right" — data-dependent
31//! stopping that p-hacks any fixed-sample test by construction.
32//! **E-processes** (nonnegative supermartingales under the null,
33//! `E[E_τ] ≤ 1` at every stopping time) are immune: by Ville's
34//! inequality `P(sup_t E_t ≥ 1/α) ≤ α`, the guarantee survives stopping
35//! whenever you like, peeking included, streaming corpora (#973)
36//! included. Evidence is a RUNNING PRODUCT, resumable across shards.
37//!
38//! 3. **This laboratory is INTERVENTIONAL.** The steering primitive with
39//! output dosimetry and a validity radius is landed
40//! (`crate::inference::steering`); the per-token output-Fisher harvest
41//! (#980) gives the local information geometry of the model's output.
42//! So "which probe next?" is not a vibe — it is OPTIMAL EXPERIMENTAL
43//! DESIGN: choose the steering intervention that maximizes the expected
44//! log-growth of the e-process deciding a contested structural claim.
45//! Under the local Gaussian output-Fisher model, that growth rate IS a
46//! KL divergence with a closed form (below). The model chooses its own
47//! next experiment, optimally, inside its certified validity radius.
48//!
49//! The deliverable shape: **every discovered atom ships an e-value; every
50//! dictionary ships an e-BH FDR certificate over its claimed structure;
51//! contested claims get design-optimal probes until evidence resolves.**
52//! No other interpretability stack has finite-sample, optional-stopping-
53//! safe error control over its discovered structure. This module is the
54//! statistical substrate; the SAE structure search plugs its gates in.
55//!
56//! The instruments, bottom-up: [`EProcess`] (running evidence, Ville
57//! semantics) → [`PredictablePluginEProcess`] (streaming universal
58//! inference) → [`AtomBirthGate`] + [`run_atom_birth_gate`] (the K vs K+1
59//! gate with demote-never-reject [`GateVerdict`]s; the runner enforces
60//! the predictability contract by call order) → [`StructureLedger`] (one
61//! e-process per claim, serializable across #973 shards) →
62//! [`StructureLedger::certify`] (the e-BH [`StructureCertificate`],
63//! shipped beside the gauge report via
64//! `crate::terms::sae::identifiability::dictionary_report`) →
65//! [`plan_probe_for_contested_claim`] (the design loop: contested claims
66//! get a [`ProbePlan`] whose δ runs through
67//! `crate::inference::steering::steer_delta` and whose per-hypothesis
68//! μ₀/μ₁ come from `crate::inference::steering::predicted_response`).
69//!
70//! # The math, fixed here so implementations cannot drift
71//!
72//! **E-value / e-process.** A nonnegative statistic `E` with `E_{H0}[E] ≤ 1`.
73//! Calibration to tests: reject at level α when `E ≥ 1/α` (Markov). An
74//! e-PROCESS compounds multiplicatively: `E_t = Π_{s≤t} e_s` where each
75//! `e_s` is conditionally valid given the past (`E[e_s | F_{s−1}] ≤ 1`
76//! under H0). Ville: `P_{H0}(∃t: E_t ≥ 1/α) ≤ α` — anytime validity.
77//! Always accumulate in log space; evidence products underflow doubles.
78//!
79//! **Universal inference (the atom-birth test).** Split the data (or the
80//! token stream) into D₀ (evaluation) and D₁ (estimation). Fit the
81//! K+1-atom alternative on D₁ by ANY method (the production fitter, warm
82//! starts, GPU, anything — no conditions). Fit the K-atom null by
83//! CONSTRAINED MLE ON D₀ (this is the one side that must be honest). Then
84//!
85//! ```text
86//! E = L(θ̂₁ ; D₀) / sup_{θ ∈ H0} L(θ ; D₀)
87//! ```
88//!
89//! satisfies `E_{H0}[E] ≤ 1` in finite samples, mixtures and boundaries
90//! and all (the proof is three lines of Markov + tower; no asymptotics).
91//! Sequential version: at each batch t, the alternative plug-in is fit on
92//! data BEFORE t (predictable), the null sup is over the batch; the
93//! product is an e-process. This is `SplitLikelihoodEValue` /
94//! `PredictablePluginEProcess` below, generic over log-likelihood
95//! closures so the SAE stack passes its own (manifold likelihoods,
96//! superposition-aware residual models #974, whatever exists then).
97//!
98//! **Bayes factors are e-values (the #907 bridge).** A Bayes factor
99//! `BF = ∫ L(θ) dΠ₁(θ) / L_{H0}` with a FIXED (data-independent) prior Π₁
100//! and a SIMPLE (or sup-dominated) null has `E_{H0}[BF] ≤ 1` — the #907
101//! geometry-adjudication harness (circle vs clusters vs line, with its
102//! discrete-mixture null) is therefore ONE PRIOR-FREEZE away from anytime
103//! validity. The integration contract: route its per-batch BFs through
104//! [`EProcess::absorb`] instead of comparing a final BF to a threshold,
105//! and geometry claims inherit optional-stopping safety for free.
106//!
107//! **e-BH (the dictionary certificate).** Given e-values e_1..e_m for m
108//! structural claims (one per atom/edge/binding), sort descending and
109//! reject the top k* where `k* = max{ k : e_(k) ≥ m/(α·k) }`. This
110//! controls FDR ≤ α under ARBITRARY dependence between the e-values
111//! (Wang–Ramdas 2022) — no independence assumptions about atoms sharing
112//! tokens, which is good because they all share every token. That
113//! arbitrary-dependence robustness is WHY the certificate uses e-BH and
114//! not a p-value BH: the p-version needs PRDS, which atom statistics
115//! flagrantly violate.
116//!
117//! **Design-optimal probing.** For a contested claim with competing
118//! structural hypotheses H₀/H₁ (e.g. "feature f is one curved atom" vs
119//! "two flat atoms"), a candidate steering intervention δ (within its
120//! certified validity radius, `steering.rs`) produces predicted output
121//! distributions P₀^δ, P₁^δ. The expected per-observation log-growth of
122//! the likelihood-ratio e-process under H₁ is exactly `KL(P₁^δ ‖ P₀^δ)`,
123//! so the optimal next experiment is
124//!
125//! ```text
126//! δ* = argmax_{δ : ‖δ‖ ≤ r_valid} KL(P₁^δ ‖ P₀^δ)
127//! ≈ argmax_δ ½ (μ₁(δ) − μ₀(δ))ᵀ F (μ₁(δ) − μ₀(δ))
128//! ```
129//!
130//! under the local Gaussian model with output-Fisher metric F (#980
131//! harvest) — the SAME quadratic form the steering dosimetry already
132//! computes, repurposed: dosimetry measures nats delivered, design
133//! maximizes nats of DISCRIMINATION. Probes that maximize raw effect are
134//! not optimal; probes that maximize the *disagreement between the
135//! hypotheses' predictions*, weighted by output information, are. (Full
136//! KL-optimal design over the probe manifold is a research arc; the greedy
137//! quadratic rule below is the correct first instrument and is exact in
138//! the local model.)
139//!
140//! # What this kills
141//!
142//! - Birth/death gates that are threshold heuristics → replaced by
143//! anytime-valid tests with declared error rates (#976's "detectable,
144//! correctable misspecification" becomes a literal hypothesis test).
145//! - "We trained until the features looked interpretable" → optional
146//! stopping is SAFE; the certificate survives it.
147//! - "We found N features" → "we found N features at FDR ≤ α, certificate
148//! attached, reproducible from the e-value ledger."
149//! - Probe selection by intuition → probe selection by information.
150
151use ndarray::{Array1, Array2};
152use serde::{Deserialize, Serialize};
153
154/// Running anytime-valid evidence against one null hypothesis, in log
155/// space. Multiplicative absorption of conditionally-valid e-values;
156/// Ville's inequality converts the running product into a sequential test
157/// that survives optional stopping. Serializable so evidence is resumable
158/// across corpus shards (#973): persist, reload, keep absorbing.
159#[derive(Clone, Debug, Serialize, Deserialize)]
160pub struct EProcess {
161 /// log E_t — the running log-evidence. Starts at 0 (E_0 = 1).
162 log_e: f64,
163 /// Number of absorbed batches (the ledger length).
164 steps: usize,
165 /// Running maximum of log E_t — Ville's inequality applies to the
166 /// supremum, so a claim once proven at level α stays proven even if
167 /// later evidence retreats (evidence is not p-hackable in reverse).
168 log_e_max: f64,
169}
170
171impl EProcess {
172 pub fn new() -> Self {
173 Self {
174 log_e: 0.0,
175 steps: 0,
176 log_e_max: 0.0,
177 }
178 }
179
180 /// Absorb one conditionally-valid e-value (NOT in log space; must be
181 /// ≥ 0; `E[e | past] ≤ 1` under H0 is the caller's contract — e.g. a
182 /// universal-inference batch ratio or a fixed-prior Bayes factor).
183 pub fn absorb(&mut self, e_value: f64) -> Result<(), String> {
184 if e_value.is_nan() || e_value < 0.0 {
185 return Err(format!("e-value must be in [0, ∞], got {e_value}"));
186 }
187 self.absorb_log(e_value.ln())
188 }
189
190 /// Absorb a batch e-value supplied in log space (the only numerically
191 /// honest interface for long streams).
192 pub fn absorb_log(&mut self, log_e_value: f64) -> Result<(), String> {
193 let next_log_e = checked_log_e_sum(self.log_e, log_e_value)?;
194 self.log_e = next_log_e;
195 self.steps += 1;
196 if self.log_e > self.log_e_max {
197 self.log_e_max = self.log_e;
198 }
199 Ok(())
200 }
201
202 pub fn log_evidence(&self) -> f64 {
203 self.log_e
204 }
205
206 pub fn steps(&self) -> usize {
207 self.steps
208 }
209
210 /// Anytime-valid rejection at level α: by Ville,
211 /// `P_{H0}(sup_t E_t ≥ 1/α) ≤ α`, so crossing 1/α at ANY time —
212 /// including data-dependent stopping times — proves the claim with
213 /// type-I error ≤ α. Uses the running supremum: once crossed, always
214 /// rejected.
215 pub fn rejects_at(&self, alpha: f64) -> bool {
216 alpha > 0.0 && self.log_e_max >= -(alpha.ln())
217 }
218
219 /// The e-value to hand to [`e_benjamini_hochberg`] for the
220 /// dictionary-level FDR certificate (current evidence, not the sup —
221 /// e-BH's guarantee is stated for e-values at the chosen stopping
222 /// time).
223 pub fn current_e_value_log(&self) -> f64 {
224 self.log_e
225 }
226}
227
228impl Default for EProcess {
229 fn default() -> Self {
230 Self::new()
231 }
232}
233
234fn checked_log_e_sum(current: f64, increment: f64) -> Result<f64, String> {
235 if current.is_nan() {
236 return Err("EProcess invariant violation: current log evidence is NaN".to_string());
237 }
238 if increment.is_nan() {
239 return Err("log e-value must not be NaN".to_string());
240 }
241 if current.is_infinite()
242 && increment.is_infinite()
243 && current.is_sign_positive() != increment.is_sign_positive()
244 {
245 return Err(format!(
246 "cannot combine opposing infinite log e-values: current {current}, increment {increment}"
247 ));
248 }
249 Ok(current + increment)
250}
251
252/// One universal-inference (split-likelihood-ratio) e-value: finite-sample
253/// valid with NO regularity conditions — the correct instrument for atom
254/// birth (K vs K+1 components, boundary/Davies regime where χ² fails).
255///
256/// `log_lik_alternative_on_eval`: log-likelihood of the EVALUATION fold
257/// under the alternative fitted on the ESTIMATION fold (any fitter — the
258/// production manifold/SAE fit, warm-started, GPU; zero conditions on it).
259/// `log_lik_null_sup_on_eval`: the SUPREMUM of the evaluation-fold
260/// log-likelihood over the NULL model class (the honest side: a real
261/// constrained fit on the eval fold, e.g. the K-atom dictionary refit on
262/// D₀). Then `log E = ℓ_alt(D₀) − sup_{H0} ℓ(D₀)` and `E_{H0}[E] ≤ 1`
263/// exactly.
264#[derive(Clone, Copy, Debug, PartialEq)]
265pub enum SplitLikelihoodError {
266 UndefinedLogLikelihood {
267 alternative: f64,
268 null_supremum: f64,
269 },
270}
271
272impl std::fmt::Display for SplitLikelihoodError {
273 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
274 match self {
275 Self::UndefinedLogLikelihood {
276 alternative,
277 null_supremum,
278 } => write!(
279 f,
280 "split-likelihood evidence requires defined log-likelihoods; alternative={alternative}, null_supremum={null_supremum}"
281 ),
282 }
283 }
284}
285
286impl std::error::Error for SplitLikelihoodError {}
287
288pub fn split_likelihood_log_e_value(
289 log_lik_alternative_on_eval: f64,
290 log_lik_null_sup_on_eval: f64,
291) -> Result<f64, SplitLikelihoodError> {
292 // NaN is a failed likelihood evaluation, not neutral evidence. Recasting
293 // it as E=1 would hide a broken source and make the serialized certificate
294 // impossible to audit.
295 if log_lik_alternative_on_eval.is_nan() || log_lik_null_sup_on_eval.is_nan() {
296 return Err(SplitLikelihoodError::UndefinedLogLikelihood {
297 alternative: log_lik_alternative_on_eval,
298 null_supremum: log_lik_null_sup_on_eval,
299 });
300 }
301 // Both models assigning exact zero density is a common-null event. The
302 // likelihood ratio is immaterial there; its conservative continuous
303 // extension is E=1 (zero log evidence). Positive infinite likelihoods do
304 // not describe normalized finite-row models and remain an error.
305 if log_lik_alternative_on_eval == f64::NEG_INFINITY
306 && log_lik_null_sup_on_eval == f64::NEG_INFINITY
307 {
308 return Ok(0.0);
309 }
310 if log_lik_alternative_on_eval == f64::INFINITY && log_lik_null_sup_on_eval == f64::INFINITY {
311 return Err(SplitLikelihoodError::UndefinedLogLikelihood {
312 alternative: log_lik_alternative_on_eval,
313 null_supremum: log_lik_null_sup_on_eval,
314 });
315 }
316 Ok(log_lik_alternative_on_eval - log_lik_null_sup_on_eval)
317}
318
319/// Sequential universal inference over a stream of batches with a
320/// PREDICTABLE plug-in: at batch t the alternative parameters were fit
321/// using only data before t, the null sup is taken on batch t, and the
322/// per-batch ratios compound into an e-process. This is the streaming /
323/// optional-stopping form the corpus-scale pipeline (#973 shards) needs —
324/// evidence is resumable: serialize `EProcess`, keep absorbing on the next
325/// shard.
326#[derive(Clone, Debug, Serialize, Deserialize)]
327pub struct PredictablePluginEProcess {
328 pub process: EProcess,
329}
330
331impl PredictablePluginEProcess {
332 pub fn new() -> Self {
333 Self {
334 process: EProcess::new(),
335 }
336 }
337
338 /// Absorb one batch. The caller guarantees the alternative was fit
339 /// WITHOUT batch-t data (predictability — this is what makes the
340 /// product a supermartingale; violating it voids the guarantee, which
341 /// is why the SAE integration must hand this function the PREVIOUS
342 /// shard's fitted dictionary, never the current one).
343 pub fn try_absorb_batch(
344 &mut self,
345 log_lik_alternative_prefit: f64,
346 log_lik_null_sup_on_batch: f64,
347 ) -> Result<(), String> {
348 let log_e =
349 split_likelihood_log_e_value(log_lik_alternative_prefit, log_lik_null_sup_on_batch)
350 .map_err(|error| error.to_string())?;
351 self.process.absorb_log(log_e)
352 }
353}
354
355impl Default for PredictablePluginEProcess {
356 fn default() -> Self {
357 Self::new()
358 }
359}
360
361/// The anytime-valid verdict on one structural claim. Deliberately
362/// two-valued — there is NO "rejected" arm. Demote-never-reject (#969
363/// philosophy): an e-process that has not crossed 1/α has failed to prove
364/// the claim, not disproven it; the claim stays contested, keeps its
365/// evidence, and earns a design-optimal probe budget instead of being
366/// silently dropped (or worse, silently accepted the way a threshold gate
367/// accepts whatever clears it).
368#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
369pub enum GateVerdict {
370 /// The running supremum crossed 1/α: the claim is proven with type-I
371 /// error ≤ α, permanently (Ville applies to the sup, so later evidence
372 /// retreat cannot un-prove it).
373 Certified { log_e: f64 },
374 /// Not (yet) proven. Carries the CURRENT log-evidence — the value the
375 /// dictionary certificate's e-BH consumes, and the state a probe loop
376 /// resumes from.
377 Contested { log_e: f64 },
378}
379
380/// The atom-birth gate (#976's threshold comparison, replaced): a
381/// universal-inference e-process over corpus shards deciding "does the
382/// K+1-th atom exist?", the boundary/Davies-regime question where the χ²
383/// gate every dictionary paper runs is broken.
384///
385/// Per shard t the integration contract is exactly the work plan's:
386/// - `log_lik_alternative_prefit`: the K+1-atom dictionary fit on shards
387/// BEFORE t (the PREVIOUS shard's fit — predictability is the one rule;
388/// handing in the current shard's fit voids the guarantee), evaluated on
389/// shard t. Any fitter, warm starts, GPU — no conditions.
390/// - `log_lik_null_sup_on_shard`: the K-atom dictionary REFIT on shard t
391/// (the honest constrained sup on the evaluation data).
392///
393/// The gate never rejects: [`GateVerdict::Contested`] is the only
394/// alternative to certification, and a contested atom's next move is a
395/// probe plan ([`plan_probe_for_contested_claim`]), not deletion.
396#[derive(Clone, Debug, Serialize, Deserialize)]
397pub struct AtomBirthGate {
398 pub test: PredictablePluginEProcess,
399 /// The level the certificate is claimed at; fixed at construction so a
400 /// verdict can never be shopped across α after seeing the evidence.
401 alpha: f64,
402 /// First-passage time: the shard count at which the running supremum
403 /// first crossed the single-hypothesis Ville bar `ln(1/α)`. `None` until
404 /// the crossing happens. This is the realized *time-to-certification* —
405 /// the design-relevant quantity (`expected_resolution_budget`) — recorded
406 /// SEPARATELY from the running evidence so that absorption can continue
407 /// past the crossing without losing the crossing time. By Ville the
408 /// crossing is permanent (it is a property of the running sup), so this is
409 /// monotone: once `Some`, never reset. `#[serde(default)]` keeps gates
410 /// persisted before this field was added (#973 resumability) loadable.
411 #[serde(default)]
412 certified_at_step: Option<usize>,
413}
414
415impl AtomBirthGate {
416 pub fn new(alpha: f64) -> Result<Self, String> {
417 if !(alpha > 0.0 && alpha < 1.0) {
418 return Err(format!(
419 "AtomBirthGate: alpha must be in (0,1), got {alpha}"
420 ));
421 }
422 Ok(Self {
423 test: PredictablePluginEProcess::new(),
424 alpha,
425 certified_at_step: None,
426 })
427 }
428
429 pub fn alpha(&self) -> f64 {
430 self.alpha
431 }
432
433 /// The realized time-to-certification: the shard count at which the
434 /// running supremum first crossed `ln(1/α)`, or `None` if it never did.
435 /// This is the first-passage time the design budget
436 /// ([`expected_resolution_budget`]) predicts — distinct from
437 /// [`EProcess::steps`] (total shards absorbed), which keeps growing after
438 /// the crossing because absorption does not stop (continuing to accumulate
439 /// is what lets the dictionary-level e-BH certificate clear its higher
440 /// multiplicity bar `ln(m/(α·k))`).
441 pub fn certified_at_step(&self) -> Option<usize> {
442 self.certified_at_step
443 }
444
445 /// Absorb one shard's split-likelihood ratio (see type-level contract).
446 pub fn try_absorb_shard(
447 &mut self,
448 log_lik_alternative_prefit: f64,
449 log_lik_null_sup_on_shard: f64,
450 ) -> Result<(), String> {
451 self.test
452 .try_absorb_batch(log_lik_alternative_prefit, log_lik_null_sup_on_shard)?;
453 // Record the first-passage time once, the step the running sup first
454 // crosses the single-hypothesis bar. `rejects_at` reads the running
455 // maximum, so this latches permanently even if later evidence retreats.
456 if self.certified_at_step.is_none() && self.test.process.rejects_at(self.alpha) {
457 self.certified_at_step = Some(self.test.process.steps());
458 }
459 Ok(())
460 }
461
462 pub fn absorb_shard(
463 &mut self,
464 log_lik_alternative_prefit: f64,
465 log_lik_null_sup_on_shard: f64,
466 ) {
467 self.try_absorb_shard(log_lik_alternative_prefit, log_lik_null_sup_on_shard)
468 .expect("AtomBirthGate received invalid log evidence");
469 }
470
471 pub fn verdict(&self) -> GateVerdict {
472 if self.test.process.rejects_at(self.alpha) {
473 GateVerdict::Certified {
474 log_e: self.test.process.log_evidence(),
475 }
476 } else {
477 GateVerdict::Contested {
478 log_e: self.test.process.log_evidence(),
479 }
480 }
481 }
482}
483
484/// Run the atom-birth gate over a shard stream with the predictability
485/// contract enforced BY CONSTRUCTION: on each shard the alternative is
486/// evaluated strictly before it is refit with that shard, so the plug-in
487/// is always predictable and the product is always a supermartingale
488/// under H0. This is the orchestration the SAE structure search calls;
489/// the closures are the only fitter-specific surface.
490///
491/// - `alternative_log_lik(alt, shard)`: evaluation-fold log-likelihood of
492/// shard under the CURRENT alternative state (fit on prior shards only —
493/// guaranteed here by call order).
494/// - `null_sup_log_lik(shard)`: the honest constrained sup — the K-atom
495/// null REFIT on this shard (any fitter; this side must genuinely
496/// maximize over H0 on the shard, an under-maximized null inflates the
497/// e-value and voids validity).
498/// - `refit_alternative(alt, shard)`: fold the shard into the alternative.
499///
500/// All three operations are fallible. An undefined likelihood or non-converged
501/// refit aborts the gate; it is never converted into neutral evidence.
502///
503/// `initial_alternative` is the K+1 fit from data BEFORE the stream (or a
504/// prior-driven init; validity never depends on its quality — a bad init
505/// only costs power). Absorbs EVERY shard's evidence into the e-process — it
506/// does NOT stop at the single-hypothesis Ville crossing. An earlier
507/// early-stop ("the crossing is permanent; further shards only cost compute")
508/// conflated two distinct bars:
509/// * the per-move VERDICT ([`AtomBirthGate::verdict`] via
510/// [`EProcess::rejects_at`]) tests the running SUPREMUM against the
511/// single-hypothesis bar `ln(1/α)`, where the crossing is indeed
512/// permanent and the realized crossing step is captured separately in
513/// [`AtomBirthGate::certified_at_step`]; but
514/// * the dictionary-level CERTIFICATE ([`StructureLedger::certify`] →
515/// [`e_benjamini_hochberg`]) consumes the CURRENT log e-value
516/// ([`EProcess::current_e_value_log`]) under a multiplicity correction
517/// whose bar `ln(m/(α·k))` is strictly higher for m > 1 claims.
518/// Stopping at the single-hypothesis bar banked just enough evidence for the
519/// move's own verdict but starved the FDR certificate, so a genuinely real
520/// atom (e.g. 0.8 nats/shard × 10 shards = 8 nats) failed e-BH confirmation
521/// against a second claim (bar `ln(2/0.05) ≈ 3.69`) even though it cleared
522/// its own `ln(1/0.05) ≈ 3.00` bar — it carried ample evidence but the early
523/// stop threw the surplus away. Continuing to absorb is always sound: the
524/// e-process is a supermartingale under H0, so additional
525/// conditionally-valid increments preserve validity, and the predictability
526/// contract is untouched (the alternative is still evaluated before it is
527/// refit with each shard). Returns the gate (verdict + resumable evidence,
528/// with `certified_at_step` holding the realized time-to-certification) and
529/// the final alternative state, which has seen the whole stream.
530pub fn run_atom_birth_gate<S, A>(
531 alpha: f64,
532 initial_alternative: A,
533 shards: impl IntoIterator<Item = S>,
534 mut alternative_log_lik: impl FnMut(&A, &S) -> Result<f64, String>,
535 mut null_sup_log_lik: impl FnMut(&S) -> Result<f64, String>,
536 mut refit_alternative: impl FnMut(A, &S) -> Result<A, String>,
537) -> Result<(AtomBirthGate, A), String> {
538 let mut gate = AtomBirthGate::new(alpha)?;
539 let mut alt = initial_alternative;
540 for shard in shards {
541 let log_lik_alt = alternative_log_lik(&alt, &shard)?;
542 let log_lik_null = null_sup_log_lik(&shard)?;
543 gate.try_absorb_shard(log_lik_alt, log_lik_null)?;
544 alt = refit_alternative(alt, &shard)?;
545 }
546 Ok((gate, alt))
547}
548
549/// e-BH: FDR control over m structural claims under ARBITRARY dependence
550/// (Wang–Ramdas). Input: per-claim log-e-values. Output: indices of
551/// rejected (i.e. CONFIRMED-STRUCTURE) claims, FDR ≤ α.
552///
553/// Sort e-values descending; reject the top k* where
554/// `k* = max{ k : e_(k) ≥ m/(α·k) }`.
555///
556/// This is the "dictionary certificate": run one e-process per claimed
557/// atom (and per claimed binding edge, #975), call this at the chosen
558/// stopping time, and the dictionary ships with an FDR-controlled list of
559/// which of its claimed structures are statistically real. No
560/// independence assumptions — atoms sharing every token is fine; that is
561/// exactly the case p-value BH cannot legally handle (PRDS violation) and
562/// e-BH can.
563#[derive(Clone, Copy, Debug, PartialEq)]
564pub enum EBhError {
565 InvalidAlpha { alpha: f64 },
566 InvalidLogEvidence { claim: usize, value: f64 },
567}
568
569impl std::fmt::Display for EBhError {
570 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
571 match self {
572 Self::InvalidAlpha { alpha } => {
573 write!(
574 f,
575 "e-BH alpha must be finite and lie in (0, 1); got {alpha}"
576 )
577 }
578 Self::InvalidLogEvidence { claim, value } => {
579 write!(
580 f,
581 "e-BH log evidence for claim {claim} must be finite or -infinity; got {value}"
582 )
583 }
584 }
585 }
586}
587
588impl std::error::Error for EBhError {}
589
590pub fn e_benjamini_hochberg(log_e_values: &[f64], alpha: f64) -> Result<Vec<usize>, EBhError> {
591 if !(alpha.is_finite() && alpha > 0.0 && alpha < 1.0) {
592 return Err(EBhError::InvalidAlpha { alpha });
593 }
594 let m = log_e_values.len();
595 if m == 0 {
596 return Ok(Vec::new());
597 }
598 // Negative infinity is meaningful: it is the exact zero e-value. NaN has
599 // no statistical meaning, and positive infinity cannot be serialized in
600 // the certificate artifact. Both fail the whole certificate rather than
601 // being silently recoded or capped to a different claim.
602 if let Some((claim, &value)) = log_e_values
603 .iter()
604 .enumerate()
605 .find(|(_, value)| value.is_nan() || **value == f64::INFINITY)
606 {
607 return Err(EBhError::InvalidLogEvidence { claim, value });
608 }
609 let mut order: Vec<usize> = (0..m).collect();
610 order.sort_by(|&a, &b| log_e_values[b].total_cmp(&log_e_values[a]));
611 let m_f = m as f64;
612 let mut k_star = 0usize;
613 for (rank0, &idx) in order.iter().enumerate() {
614 let k = (rank0 + 1) as f64;
615 // e_(k) ≥ m / (α k) ⟺ log e_(k) ≥ log m − log α − log k
616 if log_e_values[idx] >= m_f.ln() - alpha.ln() - k.ln() {
617 k_star = rank0 + 1;
618 }
619 }
620 order.truncate(k_star);
621 Ok(order)
622}
623
624/// Per-claim e-BH verdict: the descending-`log_e` rank the rule used, the
625/// confirmation flag, and the anytime-valid evidence budget
626/// `max(0, (ln m − ln α − ln k) − log_e)` measured against that SAME rank.
627#[derive(Clone, Debug, PartialEq)]
628pub struct EBhClaimVerdict {
629 /// 1-based rank of the claim in descending `log_e` order (stable on ties,
630 /// so equal `log_e` keep ascending index order).
631 pub rank: usize,
632 /// Whether the e-BH rule confirms the claim at the given `alpha`.
633 pub confirmed: bool,
634 /// Nats of additional evidence the claim would need to cross its own
635 /// rank's threshold; `0.0` once confirmed.
636 pub evidence_remaining_nats: f64,
637}
638
639/// The per-claim report companion to [`e_benjamini_hochberg`]: one owner for
640/// the rank/threshold arithmetic, so a report surface cannot drift from the
641/// rule it reports on (#2470 — the Python bindings carried a second copy of
642/// the ranking and threshold formula next to the rule itself).
643pub fn e_bh_claim_verdicts(
644 log_e_values: &[f64],
645 alpha: f64,
646) -> Result<Vec<EBhClaimVerdict>, EBhError> {
647 let confirmed: std::collections::HashSet<usize> =
648 e_benjamini_hochberg(log_e_values, alpha)?.into_iter().collect();
649 let m = log_e_values.len();
650 let mut order: Vec<usize> = (0..m).collect();
651 order.sort_by(|&a, &b| log_e_values[b].total_cmp(&log_e_values[a]));
652 let mut rank_of = vec![0usize; m];
653 for (rank0, &idx) in order.iter().enumerate() {
654 rank_of[idx] = rank0 + 1;
655 }
656 let m_f = m as f64;
657 Ok((0..m)
658 .map(|i| {
659 // e_(k) ≥ m / (α k) ⟺ log e_(k) ≥ log m − log α − log k
660 let threshold = m_f.ln() - alpha.ln() - (rank_of[i] as f64).ln();
661 EBhClaimVerdict {
662 rank: rank_of[i],
663 confirmed: confirmed.contains(&i),
664 evidence_remaining_nats: (threshold - log_e_values[i]).max(0.0),
665 }
666 })
667 .collect())
668}
669
670/// What one structural claim asserts about the dictionary. One e-process
671/// runs per claim; the kinds mirror the discovery stack's claim surface:
672/// atom existence (#976 birth), binding edges (#975), geometry
673/// adjudication (#907). `Custom` keeps the ledger open to claim types
674/// that do not exist yet without an enum churn per new discovery gate.
675#[derive(Clone, Debug, PartialEq, Eq, Hash, Serialize, Deserialize)]
676pub enum ClaimKind {
677 /// "Atom `atom` is statistically real" — the K vs K+1 birth claim.
678 AtomExists { atom: usize },
679 /// "Atoms `a` and `b` are bound" — a #975 binding edge.
680 BindingEdge { a: usize, b: usize },
681 /// "Atom `atom`'s latent geometry is `kind`" (e.g. "circle",
682 /// "clusters", "line") — a #907 adjudication claim.
683 GeometryKind { atom: usize, kind: String },
684 /// Any other structural claim, labeled.
685 Custom { label: String },
686}
687
688/// One claim plus its running evidence.
689#[derive(Clone, Debug, Serialize, Deserialize)]
690pub struct StructuralClaim {
691 pub kind: ClaimKind,
692 pub evidence: EProcess,
693}
694
695/// The dictionary's claim ledger: every structural claim the discovery
696/// stack makes, each with its own e-process. Serializable — evidence
697/// resumes across corpus shards (#973) by persisting the ledger, not by
698/// refitting. Calling [`StructureLedger::certify`] at ANY data-dependent
699/// stopping time yields a valid certificate; that is the entire point.
700#[derive(Clone, Debug, Default, Serialize, Deserialize)]
701pub struct StructureLedger {
702 claims: Vec<StructuralClaim>,
703}
704
705impl StructureLedger {
706 pub fn new() -> Self {
707 Self { claims: Vec::new() }
708 }
709
710 /// Register a claim and return its ledger index. Idempotent on the
711 /// claim kind: re-registering an existing claim (a resumed shard loop
712 /// re-announcing its claim surface) returns the existing index and
713 /// PRESERVES its accumulated evidence — a fresh e-process here would
714 /// silently discard the stream's history.
715 pub fn register(&mut self, kind: ClaimKind) -> usize {
716 if let Some(idx) = self.claims.iter().position(|c| c.kind == kind) {
717 return idx;
718 }
719 self.claims.push(StructuralClaim {
720 kind,
721 evidence: EProcess::new(),
722 });
723 self.claims.len() - 1
724 }
725
726 /// Absorb one conditionally-valid log e-value for claim `idx` (a
727 /// universal-inference shard ratio, a frozen-prior log-BF — the
728 /// caller's contract is per-source, documented on the producing gate).
729 pub fn absorb_log(&mut self, idx: usize, log_e_value: f64) -> Result<(), String> {
730 let n = self.claims.len();
731 let claim = self.claims.get_mut(idx).ok_or_else(|| {
732 format!("StructureLedger: claim index {idx} out of range ({n} claims)")
733 })?;
734 claim.evidence.absorb_log(log_e_value)
735 }
736
737 pub fn claims(&self) -> &[StructuralClaim] {
738 &self.claims
739 }
740
741 /// The likelihood half of the probe-design loop (work-plan step 4):
742 /// after running a planned probe ([`ProbePlan`] →
743 /// `crate::inference::steering::steer_delta`), evaluate the REALIZED
744 /// outcomes under both hypotheses' predictive densities and absorb the
745 /// log-ratio into the contested claim's e-process.
746 ///
747 /// Validity contract: both predictive densities must be FROZEN before the
748 /// probe outcome is observed — which the design loop satisfies by
749 /// construction, since both hypotheses' dictionaries were fitted before
750 /// the probe was even chosen. For a composite null, the null density must
751 /// be the honest constrained fit (the same rule as
752 /// [`split_likelihood_log_e_value`], which this delegates to); for a
753 /// simple null the predictive density is the sup. Probe outcomes are new
754 /// data by construction (the model was steered to produce them), so they
755 /// compound validly with the claim's prior shard evidence.
756 pub fn absorb_probe_outcome(
757 &mut self,
758 idx: usize,
759 log_lik_alt_on_outcome: f64,
760 log_lik_null_on_outcome: f64,
761 ) -> Result<(), String> {
762 let log_e = split_likelihood_log_e_value(log_lik_alt_on_outcome, log_lik_null_on_outcome)
763 .map_err(|error| error.to_string())?;
764 self.absorb_log(idx, log_e)
765 }
766
767 /// The dictionary certificate: e-BH over the ledger's CURRENT
768 /// e-values at level α. FDR ≤ α over the confirmed set under arbitrary
769 /// dependence — atoms sharing every token is fine — and valid at any
770 /// stopping time because each entry is an e-process. Claims not
771 /// confirmed are CONTESTED, never rejected (demote-never-reject); they
772 /// keep their evidence and are the inputs to the probe-design loop.
773 pub fn certify(&self, alpha: f64) -> Result<StructureCertificate, EBhError> {
774 let log_e: Vec<f64> = self
775 .claims
776 .iter()
777 .map(|c| c.evidence.current_e_value_log())
778 .collect();
779 let confirmed_idx = e_benjamini_hochberg(&log_e, alpha)?;
780 let mut entries: Vec<CertificateEntry> = self
781 .claims
782 .iter()
783 .zip(&log_e)
784 .map(|(c, &le)| CertificateEntry {
785 kind: c.kind.clone(),
786 log_e: le,
787 steps: c.evidence.steps(),
788 confirmed: false,
789 })
790 .collect();
791 for &i in &confirmed_idx {
792 entries[i].confirmed = true;
793 }
794 Ok(StructureCertificate { alpha, entries })
795 }
796}
797
798/// One line of the certificate's e-value ledger: the claim, its
799/// log-evidence at the stop, how many batches produced it, and the e-BH
800/// outcome. The full entry list IS the reproducibility artifact: anyone
801/// holding it can re-run [`e_benjamini_hochberg`] and re-derive the
802/// confirmed set.
803#[derive(Clone, Debug, Serialize, Deserialize)]
804pub struct CertificateEntry {
805 pub kind: ClaimKind,
806 pub log_e: f64,
807 pub steps: usize,
808 pub confirmed: bool,
809}
810
811/// The deliverable: "we found N structures at FDR ≤ α, certificate
812/// attached". Ships next to the identifiability certificate
813/// (`crate::terms::sae::identifiability::residual_gauge`, #981) — that one says
814/// what the GAUGE cannot distinguish, this one says what the DATA can.
815#[derive(Clone, Debug, Serialize, Deserialize)]
816pub struct StructureCertificate {
817 pub alpha: f64,
818 pub entries: Vec<CertificateEntry>,
819}
820
821impl StructureCertificate {
822
823}
824
825/// Calibrate one (super)uniform p-value into a single e-value, in log
826/// space: `e(p) = ½ p^{−1/2}` (the κ = ½ member of the calibrator family
827/// `e_κ(p) = κ p^{κ−1}`; `∫₀¹ e_κ(p) dp = 1`, so `E_{H0}[e(P)] ≤ 1` for
828/// any valid p — superuniformity only, no other conditions).
829///
830/// This is the bridge from p-value-shaped instruments into the ledger —
831/// e.g. the feature-binding Wald test (`terms::structure::anova_atom::carve`'s
832/// `edge_p_value` → a [`ClaimKind::BindingEdge`] entry). It spends
833/// calibration slack (a p of 0.01 becomes e = 5, not 100), which is the
834/// honest price of converting a fixed-sample test into anytime-valid
835/// currency; instruments that can produce e-values natively should.
836/// CONTRACT: one calibrated e-value per INDEPENDENT data batch — feeding
837/// repeated tests of the same accumulating data into one e-process is the
838/// p-hacking this module exists to kill.
839pub fn log_e_from_p_calibrator(p_value: f64) -> Result<f64, String> {
840 if !(p_value > 0.0) || p_value > 1.0 {
841 return Err(format!("p-value must be in (0, 1], got {p_value}"));
842 }
843 Ok(0.5f64.ln() - 0.5 * p_value.ln())
844}
845
846/// A candidate steering probe for resolving one contested structural
847/// claim: the intervention direction (in the steering primitive's
848/// coordinates), and the two hypotheses' PREDICTED output-mean responses
849/// to it.
850pub struct CandidateProbe {
851 /// Steering displacement δ, to be applied via
852 /// `crate::inference::steering` (which enforces its own validity
853 /// radius and reports realized dosimetry).
854 pub delta: Array1<f64>,
855 /// Predicted output-mean response under the null structure, μ₀(δ).
856 pub predicted_mean_null: Array1<f64>,
857 /// Predicted output-mean response under the alternative, μ₁(δ).
858 pub predicted_mean_alt: Array1<f64>,
859}
860
861/// Greedy KL-optimal experimental design under the local Gaussian
862/// output-Fisher model: pick the probe maximizing
863/// `½ (μ₁(δ) − μ₀(δ))ᵀ F (μ₁(δ) − μ₀(δ))` — the expected per-observation
864/// log-growth of the deciding e-process under the alternative.
865///
866/// `fisher` is the output-Fisher metric at the operating point (#980
867/// harvest; the same object steering dosimetry contracts against). Probes
868/// whose hypotheses predict the SAME response score zero no matter how
869/// large their raw effect — the design rule selects for DISCRIMINATION,
870/// not impact, which is the entire point: a maximally-steered output that
871/// both hypotheses predict identically teaches nothing.
872///
873/// Returns the index of the best probe and its expected log-growth (nats
874/// per observation), or None if no probe discriminates.
875pub fn select_probe_by_expected_evidence(
876 probes: &[CandidateProbe],
877 fisher: &Array2<f64>,
878) -> Option<(usize, f64)> {
879 let mut best: Option<(usize, f64)> = None;
880 for (idx, probe) in probes.iter().enumerate() {
881 let diff = &probe.predicted_mean_alt - &probe.predicted_mean_null;
882 if diff.len() != fisher.nrows() {
883 continue;
884 }
885 let f_diff = fisher.dot(&diff);
886 let growth = 0.5 * diff.dot(&f_diff);
887 if growth.is_finite() && growth > 0.0 {
888 match best {
889 Some((_, g)) if g >= growth => {}
890 _ => best = Some((idx, growth)),
891 }
892 }
893 }
894 best
895}
896
897/// Expected number of observations for the chosen probe to push a claim's
898/// e-process across the 1/α Ville threshold, under the alternative: the
899/// design-time budget `log(1/α) / growth_rate`. This is what turns the
900/// abstract guarantee into an experiment plan ("this probe should resolve
901/// the claim in ~N tokens; if it hasn't, the alternative is weaker than
902/// hypothesized — itself evidence").
903pub fn expected_resolution_budget(alpha: f64, growth_nats_per_obs: f64) -> Option<f64> {
904 if alpha <= 0.0 || alpha >= 1.0 || growth_nats_per_obs <= 0.0 {
905 return None;
906 }
907 Some(-(alpha.ln()) / growth_nats_per_obs)
908}
909
910/// The experiment plan for one contested claim: which probe to run, the
911/// expected per-observation evidence growth under the alternative, and the
912/// design-time resolution budget. This is the loop's actionable output —
913/// hand `probes[probe]`'s δ to `crate::inference::steering::steer_delta`
914/// (which enforces the validity radius and reports realized dosimetry),
915/// evaluate both hypotheses' likelihoods on the realized outputs, absorb
916/// the log-ratio into the claim's e-process, re-certify.
917#[derive(Clone, Debug, PartialEq)]
918pub struct ProbePlan {
919 /// Index into the candidate probe list.
920 pub probe: usize,
921 /// Expected log-growth of the deciding e-process, nats/observation,
922 /// under the alternative (the KL of the hypotheses' predicted
923 /// responses in the output-Fisher metric).
924 pub expected_log_growth: f64,
925 /// Expected observations to cross 1/α from ZERO evidence — the
926 /// conservative from-scratch budget.
927 pub budget_from_scratch: f64,
928 /// Expected observations to cross 1/α from the claim's CURRENT
929 /// log-evidence — the remaining budget; 0 when already across.
930 pub budget_remaining: f64,
931}
932
933/// Close the design loop for one contested claim: pick the probe whose
934/// predicted hypothesis-disagreement (not raw effect) buys evidence
935/// fastest, and convert the claim's current evidence into a remaining
936/// budget — "this probe should resolve the claim in ~N more observations
937/// at level α; if it does not, the alternative is weaker than
938/// hypothesized, which is itself evidence."
939///
940/// `current_log_e` is the contested claim's running log-evidence (from its
941/// [`StructuralClaim`] / [`GateVerdict::Contested`]). Returns `None` when
942/// no probe discriminates (all candidates score zero growth: the
943/// hypotheses agree on everything reachable inside the validity radius —
944/// the claim is undecidable by steering and needs a different instrument,
945/// which is a finding, not a failure).
946pub fn plan_probe_for_contested_claim(
947 probes: &[CandidateProbe],
948 fisher: &Array2<f64>,
949 alpha: f64,
950 current_log_e: f64,
951) -> Option<ProbePlan> {
952 let (probe, expected_log_growth) = select_probe_by_expected_evidence(probes, fisher)?;
953 let budget_from_scratch = expected_resolution_budget(alpha, expected_log_growth)?;
954 let nats_remaining = (-(alpha.ln()) - current_log_e).max(0.0);
955 Some(ProbePlan {
956 probe,
957 expected_log_growth,
958 budget_from_scratch,
959 budget_remaining: nats_remaining / expected_log_growth,
960 })
961}
962
963#[cfg(test)]
964mod tests {
965 use super::*;
966 use ndarray::array;
967
968 /// e-BH on a hand-checkable configuration.
969 #[test]
970 fn e_bh_rejects_exactly_the_qualifying_prefix() {
971 // m = 4, α = 0.1 → thresholds m/(αk) = 40, 20, 13.33, 10.
972 let log_e = [45.0f64.ln(), 21.0f64.ln(), 12.0f64.ln(), 1.0f64.ln()];
973 let rejected = e_benjamini_hochberg(&log_e, 0.1).unwrap();
974 // e_(1)=45 ≥ 40 ✓, e_(2)=21 ≥ 20 ✓, e_(3)=12 < 13.33 ✗ → k* = 2.
975 assert_eq!(rejected, vec![0, 1]);
976
977 // A weaker tail cannot drag in a stronger prefix decision.
978 let log_e2 = [45.0f64.ln(), 5.0f64.ln(), 2.0f64.ln(), 1.0f64.ln()];
979 assert_eq!(e_benjamini_hochberg(&log_e2, 0.1).unwrap(), vec![0]);
980 }
981
982 #[test]
983 fn split_likelihood_equal_impossibility_is_neutral_log_evidence() {
984 let log_e = split_likelihood_log_e_value(f64::NEG_INFINITY, f64::NEG_INFINITY).unwrap();
985 assert_eq!(log_e, 0.0);
986 assert!(log_e.is_finite());
987
988 let mut proc = EProcess::new();
989 proc.absorb_log(log_e).unwrap();
990 assert_eq!(proc.log_evidence(), 0.0);
991 assert_eq!(proc.steps(), 1);
992 }
993
994 #[test]
995 fn e_bh_accepts_exact_zero_evidence_but_refuses_positive_infinity() {
996 let log_e = [f64::NEG_INFINITY, 45.0f64.ln(), 1.0f64.ln()];
997 assert_eq!(e_benjamini_hochberg(&log_e, 0.1).unwrap(), vec![1]);
998 assert!(matches!(
999 e_benjamini_hochberg(&[f64::INFINITY], 0.1),
1000 Err(EBhError::InvalidLogEvidence { claim: 0, value }) if value == f64::INFINITY
1001 ));
1002 }
1003
1004 /// A NaN has no e-value meaning. The whole certificate must fail at its
1005 /// deterministic smallest bad claim rather than silently changing that
1006 /// claim to exact zero evidence.
1007 #[test]
1008 fn e_bh_refuses_nan_at_the_certificate_boundary() {
1009 let error = e_benjamini_hochberg(&[45.0f64.ln(), f64::NAN], 0.1).unwrap_err();
1010 assert!(matches!(
1011 error,
1012 EBhError::InvalidLogEvidence { claim: 1, value } if value.is_nan()
1013 ));
1014 }
1015
1016 #[test]
1017 fn e_bh_refuses_invalid_levels_even_for_an_empty_family() {
1018 for alpha in [0.0, 1.0, -0.1, f64::NAN, f64::INFINITY] {
1019 assert!(matches!(
1020 e_benjamini_hochberg(&[], alpha),
1021 Err(EBhError::InvalidAlpha { .. })
1022 ));
1023 }
1024 }
1025
1026 /// The full source→consumer chain: a shard with zero density under both
1027 /// the alternative and the null produces `(−∞) − (−∞)`, which the split-LR
1028 /// resolves to neutral `log E = 0` rather than NaN; banking it and
1029 /// certifying must not panic. A genuinely NaN log-likelihood is refused at
1030 /// the source and never reaches the ledger.
1031 #[test]
1032 fn common_zero_density_is_neutral_but_nan_is_refused() {
1033 // (−∞) − (−∞): zero density under both hypotheses → neutral.
1034 let neutral = split_likelihood_log_e_value(f64::NEG_INFINITY, f64::NEG_INFINITY).unwrap();
1035 assert_eq!(neutral, 0.0);
1036 // A NaN is a failed likelihood evaluation, not evidence that may be
1037 // silently rewritten as neutral.
1038 assert!(split_likelihood_log_e_value(f64::NAN, -3.0).is_err());
1039
1040 let mut ledger = StructureLedger::new();
1041 let degenerate = ledger.register(ClaimKind::AtomExists { atom: 0 });
1042 let strong = ledger.register(ClaimKind::AtomExists { atom: 1 });
1043 // Bank the neutral split-LR on the degenerate claim — no NaN reaches
1044 // the e-process.
1045 ledger.absorb_log(degenerate, neutral).unwrap();
1046 ledger.absorb_log(strong, 45.0f64.ln()).unwrap();
1047 // certify() runs e_benjamini_hochberg internally; must not panic.
1048 let certificate = ledger.certify(0.1).unwrap();
1049 let degenerate_entry = certificate
1050 .entries
1051 .iter()
1052 .find(|e| e.kind == ClaimKind::AtomExists { atom: 0 })
1053 .expect("degenerate claim present");
1054 // Neutral evidence (log_e = 0) never qualifies → contested, not confirmed.
1055 assert!(!degenerate_entry.confirmed);
1056 assert_eq!(degenerate_entry.log_e, 0.0);
1057 }
1058
1059 #[test]
1060 fn e_process_absorb_log_rejects_undefined_log_products() {
1061 let mut proc = EProcess::new();
1062 assert!(proc.absorb_log(f64::NAN).is_err());
1063
1064 proc.absorb_log(f64::INFINITY).unwrap();
1065 assert!(proc.absorb_log(f64::NEG_INFINITY).is_err());
1066 assert_eq!(proc.log_evidence(), f64::INFINITY);
1067 assert_eq!(proc.steps(), 1);
1068 }
1069
1070 /// Ville-style sanity: under H0 (simulated fair e-values from a
1071 /// likelihood ratio of identical Gaussians), the e-process crosses
1072 /// 1/α rarely; under a true alternative it crosses fast and the
1073 /// crossing is PERMANENT (running-sup semantics).
1074 #[test]
1075 fn e_process_crossing_is_permanent_and_directional() {
1076 // Deterministic "stream": per-batch log-LR of N(μ,1) vs N(0,1)
1077 // evaluated at x drawn from the alternative: log e = μ x − μ²/2.
1078 // Use a fixed quasi-random sequence; no RNG state needed.
1079 let mu = 0.6f64;
1080 let mut proc_alt = EProcess::new();
1081 let mut crossed_at: Option<usize> = None;
1082 for t in 0..200 {
1083 // x_t ~ alternative-ish deterministic surrogate around μ
1084 let x = mu + 0.9 * ((t as f64 * 0.7321).sin());
1085 proc_alt.absorb_log(mu * x - 0.5 * mu * mu).unwrap();
1086 if proc_alt.rejects_at(0.05) && crossed_at.is_none() {
1087 crossed_at = Some(t);
1088 }
1089 }
1090 let t_cross = crossed_at.expect("true alternative must cross 1/α");
1091 assert!(t_cross < 100, "evidence should accumulate quickly");
1092 // Permanence: rejection holds at the end even if late evidence dips.
1093 assert!(proc_alt.rejects_at(0.05));
1094
1095 // Null stream: x centered at 0 → expected log e = −μ²/2 < 0.
1096 let mut proc_null = EProcess::new();
1097 for t in 0..200 {
1098 let x = 0.9 * ((t as f64 * 0.7321).sin());
1099 proc_null.absorb_log(mu * x - 0.5 * mu * mu).unwrap();
1100 }
1101 assert!(
1102 !proc_null.rejects_at(0.05),
1103 "null stream must not accumulate evidence (log E = {:.3})",
1104 proc_null.log_evidence()
1105 );
1106 }
1107
1108 /// The design rule selects discrimination, not raw effect.
1109 #[test]
1110 fn probe_selection_prefers_discrimination_over_impact() {
1111 let fisher = array![[2.0, 0.0], [0.0, 0.5]];
1112 let probes = vec![
1113 // Huge effect, but both hypotheses predict it identically.
1114 CandidateProbe {
1115 delta: array![1.0, 0.0],
1116 predicted_mean_null: array![10.0, 10.0],
1117 predicted_mean_alt: array![10.0, 10.0],
1118 },
1119 // Modest effect, hypotheses disagree along the informative axis.
1120 CandidateProbe {
1121 delta: array![0.0, 1.0],
1122 predicted_mean_null: array![0.0, 0.0],
1123 predicted_mean_alt: array![1.0, 0.2],
1124 },
1125 ];
1126 let (idx, growth) =
1127 select_probe_by_expected_evidence(&probes, &fisher).expect("a probe discriminates");
1128 assert_eq!(idx, 1);
1129 // ½·(1,0.2)ᵀ diag(2,0.5) (1,0.2) = ½·(2 + 0.02) = 1.01 nats/obs.
1130 assert!((growth - 1.01).abs() < 1e-12);
1131 // Budget: ~3 observations to certify at α=0.05.
1132 let budget = expected_resolution_budget(0.05, growth).expect("budget");
1133 assert!(budget > 2.0 && budget < 4.0);
1134 }
1135
1136 /// The birth gate certifies under a true alternative, stays contested
1137 /// under the null, and never emits anything but those two verdicts.
1138 #[test]
1139 fn birth_gate_certifies_alternative_and_demotes_never_rejects() {
1140 let mut gate = AtomBirthGate::new(0.05).expect("valid alpha");
1141 // Strong shards: alternative beats the honest null sup by 1 nat each.
1142 for _ in 0..5 {
1143 gate.absorb_shard(-100.0, -101.0);
1144 }
1145 match gate.verdict() {
1146 GateVerdict::Certified { log_e } => assert!((log_e - 5.0).abs() < 1e-12),
1147 v => panic!("5 nats must certify at α=0.05, got {v:?}"),
1148 }
1149 // Permanence: a later evidence retreat cannot un-certify.
1150 gate.absorb_shard(-110.0, -100.0);
1151 assert!(matches!(gate.verdict(), GateVerdict::Certified { .. }));
1152
1153 // Null-ish stream: the prefit alternative loses to the on-shard sup
1154 // (it must, on average — the sup is fit on the eval shard itself).
1155 let mut null_gate = AtomBirthGate::new(0.05).expect("valid alpha");
1156 for _ in 0..50 {
1157 null_gate.absorb_shard(-100.3, -100.0);
1158 }
1159 match null_gate.verdict() {
1160 GateVerdict::Contested { log_e } => assert!(log_e < 0.0),
1161 v => panic!("null stream must stay contested, got {v:?}"),
1162 }
1163 assert!(AtomBirthGate::new(0.0).is_err());
1164 assert!(AtomBirthGate::new(1.0).is_err());
1165 }
1166
1167 /// Resumability: a serialized ledger reloads with its evidence and
1168 /// keeps absorbing — the #973 shard contract.
1169 #[test]
1170 fn ledger_evidence_resumes_across_serialization() {
1171 let mut ledger = StructureLedger::new();
1172 let idx = ledger.register(ClaimKind::GeometryKind {
1173 atom: 3,
1174 kind: "circle".to_string(),
1175 });
1176 ledger.absorb_log(idx, 1.25).unwrap();
1177
1178 let persisted = serde_json::to_string(&ledger).expect("serialize ledger");
1179 let mut resumed: StructureLedger =
1180 serde_json::from_str(&persisted).expect("deserialize ledger");
1181 assert_eq!(resumed.claims()[idx].evidence.steps(), 1);
1182
1183 resumed.absorb_log(idx, 0.75).unwrap();
1184 let log_e = resumed.claims()[idx].evidence.log_evidence();
1185 assert!((log_e - 2.0).abs() < 1e-12);
1186 }
1187
1188 /// The probe plan discounts the remaining budget by evidence already
1189 /// banked, and floors at zero once the claim is across the line.
1190 #[test]
1191 fn probe_plan_discounts_remaining_budget_by_current_evidence() {
1192 let fisher = array![[2.0, 0.0], [0.0, 0.5]];
1193 let probes = vec![CandidateProbe {
1194 delta: array![0.0, 1.0],
1195 predicted_mean_null: array![0.0, 0.0],
1196 predicted_mean_alt: array![1.0, 0.2],
1197 }];
1198 // growth = 1.01 nats/obs (checked above); α=0.05 → need ln(20) ≈ 3.0 nats.
1199 let from_zero = plan_probe_for_contested_claim(&probes, &fisher, 0.05, 0.0).expect("plan");
1200 assert_eq!(from_zero.probe, 0);
1201 assert!((from_zero.budget_remaining - from_zero.budget_from_scratch).abs() < 1e-12);
1202
1203 let halfway = plan_probe_for_contested_claim(&probes, &fisher, 0.05, 1.5).expect("plan");
1204 assert!(halfway.budget_remaining < from_zero.budget_remaining);
1205 assert!((halfway.budget_remaining - (-(0.05f64.ln()) - 1.5) / 1.01).abs() < 1e-12);
1206
1207 let across = plan_probe_for_contested_claim(&probes, &fisher, 0.05, 10.0).expect("plan");
1208 assert_eq!(across.budget_remaining, 0.0);
1209
1210 // No discriminating probe → no plan (undecidable by steering).
1211 let blind = vec![CandidateProbe {
1212 delta: array![1.0, 0.0],
1213 predicted_mean_null: array![5.0, 5.0],
1214 predicted_mean_alt: array![5.0, 5.0],
1215 }];
1216 assert!(plan_probe_for_contested_claim(&blind, &fisher, 0.05, 0.0).is_none());
1217 }
1218
1219 /// The p→e calibrator on hand-checkable values, including its edges.
1220 #[test]
1221 fn p_to_e_calibrator_hand_values() {
1222 // e(p) = ½ p^{−1/2}: p = 1 → e = 0.5; p = 0.04 → e = 2.5; p = 1e-4 → e = 50.
1223 assert!((log_e_from_p_calibrator(1.0).unwrap() - 0.5f64.ln()).abs() < 1e-12);
1224 assert!((log_e_from_p_calibrator(0.04).unwrap() - 2.5f64.ln()).abs() < 1e-12);
1225 assert!((log_e_from_p_calibrator(1e-4).unwrap() - 50.0f64.ln()).abs() < 1e-12);
1226 assert!(log_e_from_p_calibrator(0.0).is_err());
1227 assert!(log_e_from_p_calibrator(1.5).is_err());
1228 assert!(log_e_from_p_calibrator(f64::NAN).is_err());
1229 }
1230
1231 /// The e-value validity condition: under the null `P ~ Uniform(0, 1]`,
1232 /// the calibrated e-value must satisfy `E_{H0}[e(P)] = ∫₀¹ e(p) dp ≤ 1`
1233 /// (Wang–Ramdas e-BH controls FDR ONLY for genuine e-values). The κ = ½
1234 /// member `e(p) = ½ p^{−1/2}` integrates to exactly 1, the boundary of
1235 /// admissibility. We verify this numerically with a midpoint Riemann sum
1236 /// over the unit interval (which UNDER-estimates `∫ p^{−1/2}` slightly
1237 /// because the integrand is convex, so the analytic value 1 sits just
1238 /// above the quadrature estimate — both safely ≤ 1 + tolerance).
1239 ///
1240 /// This is the property `e = 1/p` VIOLATES — `∫₀¹ (1/p) dp = ∞` — which
1241 /// is why the behavioral-head and anova-atom paths route through this
1242 /// calibrator rather than `−ln p`.
1243 #[test]
1244 fn p_to_e_calibrator_null_expectation_at_most_one() {
1245 let n = 2_000_000usize;
1246 let h = 1.0 / n as f64;
1247 let mut mean_e = 0.0_f64;
1248 for i in 0..n {
1249 // Midpoint of cell i: p = (i + 0.5)/n, always in (0, 1).
1250 let p = (i as f64 + 0.5) * h;
1251 let e = log_e_from_p_calibrator(p).unwrap().exp();
1252 mean_e += e * h;
1253 }
1254 // Analytic E_{H0}[e(P)] = 1; allow a small quadrature tolerance, but
1255 // it MUST NOT exceed 1 by more than that (an invalid calibrator like
1256 // 1/p would diverge here, not land near 1).
1257 assert!(
1258 mean_e <= 1.0 + 1e-3,
1259 "calibrated e-value null expectation {mean_e} exceeds 1 — not a valid e-value"
1260 );
1261 assert!(
1262 mean_e > 0.99,
1263 "calibrated e-value null expectation {mean_e} far below the analytic 1.0"
1264 );
1265 }
1266
1267 /// POWER STUDY, null side: the heuristic gate every dictionary paper
1268 /// runs — "accept the K+1-th atom the first time the cumulative
1269 /// likelihood ratio shows improvement" — versus the e-gate, on a
1270 /// family of NULL streams peeked at after every shard. The per-shard
1271 /// log-LR is `μ x_t − μ²/2` with `x_t = A sin(ω t + φ)` (a
1272 /// deterministic null surrogate: mean drift −μ²/2 < 0, bounded
1273 /// fluctuation). The naive gate's false-accept mechanism is exactly
1274 /// optional stopping: any phase whose partial sums wander above zero
1275 /// at ANY peek accepts a nonexistent atom. The e-gate needs
1276 /// log(1/α) ≈ 3.0 nats, and the partial-sum fluctuation is bounded by
1277 /// `μ·A/sin(ω/2) ≈ 1.51` nats (Dirichlet-kernel bound) BEFORE the
1278 /// negative drift — so it can never certify on any phase, which is
1279 /// Ville's inequality made concrete.
1280 #[test]
1281 fn power_study_null_naive_peeking_gate_false_accepts_e_gate_never() {
1282 let mu = 0.6f64;
1283 let amp = 0.9f64;
1284 let omega = 0.7321f64;
1285 let n_phases = 60usize;
1286 let n_shards = 200usize;
1287
1288 let mut naive_false_accepts = 0usize;
1289 let mut e_gate_false_accepts = 0usize;
1290 for k in 0..n_phases {
1291 let phase = 2.0 * std::f64::consts::PI * (k as f64) / (n_phases as f64);
1292 let mut gate = AtomBirthGate::new(0.05).expect("alpha");
1293 let mut cum_log_lr = 0.0f64;
1294 let mut naive_accepted = false;
1295 for t in 0..n_shards {
1296 let x = amp * ((t as f64) * omega + phase).sin();
1297 let log_lr = mu * x - 0.5 * mu * mu;
1298 cum_log_lr += log_lr;
1299 // The broken test: peek, accept on any improvement.
1300 if cum_log_lr > 0.0 {
1301 naive_accepted = true;
1302 }
1303 gate.absorb_shard(log_lr, 0.0);
1304 }
1305 if naive_accepted {
1306 naive_false_accepts += 1;
1307 }
1308 if matches!(gate.verdict(), GateVerdict::Certified { .. }) {
1309 e_gate_false_accepts += 1;
1310 }
1311 }
1312 // The naive gate false-accepts on a large fraction of null phases
1313 // (any phase with early-positive partial sums); the e-gate on none.
1314 assert!(
1315 naive_false_accepts >= n_phases / 3,
1316 "the peeking gate should false-accept often under the null \
1317 (got {naive_false_accepts}/{n_phases})"
1318 );
1319 assert_eq!(
1320 e_gate_false_accepts, 0,
1321 "the e-gate must never certify under the null"
1322 );
1323 }
1324
1325 /// POWER STUDY, alternative side, through the orchestration harness:
1326 /// a planted K+1-th atom worth 0.5 nats/shard CROSSES the single-hypothesis
1327 /// Ville bar in ⌈log(1/α)/0.5⌉ = 6 shards — the realized
1328 /// time-to-certification matching the design-time
1329 /// `expected_resolution_budget`. The gate keeps absorbing past the
1330 /// crossing (the crossing time is latched in `certified_at_step`, so it is
1331 /// not lost), banking the full stream's evidence — which is exactly what
1332 /// the dictionary-level e-BH certificate needs to clear its higher
1333 /// multiplicity bar. The alternative refits on every shard regardless.
1334 #[test]
1335 fn power_study_planted_atom_certifies_at_the_predicted_budget() {
1336 let growth = 0.5f64;
1337 let n_shards = 20usize;
1338 let (gate, alt_state) = run_atom_birth_gate(
1339 0.05,
1340 0usize, // alt state = number of shards folded into the fit
1341 0..n_shards,
1342 |_, _| Ok(-99.5), // prefit alternative log-lik on the shard
1343 |_| Ok(-100.0), // honest null sup on the shard
1344 |folded, _| Ok(folded + 1),
1345 )
1346 .expect("valid alpha");
1347
1348 // Full-stream evidence is banked: 0.5 nats/shard × 20 shards = 10 nats,
1349 // far past the single-hypothesis bar — this surplus is what the e-BH
1350 // dictionary certificate consumes against its `ln(m/(α·k))` bar.
1351 match gate.verdict() {
1352 GateVerdict::Certified { log_e } => {
1353 assert!((log_e - growth * n_shards as f64).abs() < 1e-12)
1354 }
1355 v => panic!("planted atom must certify, got {v:?}"),
1356 }
1357 // Realized time-to-certification (first-passage of the running sup over
1358 // `ln(1/α)`) == the design-time budget, rounded up.
1359 let budget = expected_resolution_budget(0.05, growth).expect("budget");
1360 assert_eq!(gate.certified_at_step(), Some(budget.ceil() as usize));
1361 assert_eq!(gate.certified_at_step(), Some(6));
1362 // Absorption does NOT stop at the crossing: every shard is banked.
1363 assert_eq!(gate.test.process.steps(), n_shards);
1364 // The alternative state saw the whole stream.
1365 assert_eq!(alt_state, n_shards);
1366 }
1367
1368}