Upload folder using huggingface_hub

COLLABORATION.md

@@ -0,0 +1,34 @@
+# High-compute collaboration call
+
+SPECTRA-RSI needs independent scaling evidence. The manuscript explicitly treats its numerical gains as hypotheses to test, with dense ground-truth deltas, causal-attribution metrics, optional-stopping stress tests, drift tests, and 20–100-cycle long-horizon studies.
+
+## Highest-value contributions
+
+1. Run larger `n_slices`, `n_experts`, rank, bootstrap, and multi-seed sweeps and upload the resulting JSON files.
+2. Port the hot numerical paths to PyTorch/JAX/CuPy and benchmark CPU vs single-GPU vs multi-GPU without changing estimands.
+3. Add real open-weight LoRA/MoE-LoRA adapters implementing the `paired_scores` and `probe_delta` interface.
+4. Test two model families and two scales, including single-expert gains/regressions, canceling mixtures, interacting experts, non-compressible changes, and out-of-trust-region updates.
+5. Run optional-stopping, dictionary-drift, Fisher/second-order dictionary regularization, and 20–100 update-cycle experiments.
+
+## Hardware sought
+
+H100/H200, B100/B200, GB200/GB300, MI300X, 2/4/8/16+ GPU servers, and multi-node clusters. Other accelerators are welcome when the environment is fully reported.
+
+## Submission format (no Git required)
+
+Download this Hugging Face repository, run the benchmark, then upload result JSON files back to a Hugging Face discussion/community contribution or share them with the maintainer. Include hardware, software versions, command line, wall time, and any code patch as a ZIP if you changed the backend. Do not submit only a screenshot.
+
+## Current reference point and open scaling question
+
+The validated reference point is **1 × NVIDIA GB10 (DGX Spark)** at 1,024 capability slices. Sparse single-gain and canceling-mixture cases recover their true expert support exactly in the validated runs, while the genuine `off_dictionary` case triggers immediate dense fallback.
+
+The most important current stress-test result is deliberately unresolved: for `broad_noncompressible` at candidate seed 100, the gate returned **accept** while only 6 of 8 true experts were recovered. Contributors should not tune around or discard this case. We specifically want to learn whether increased sensing budgets, bootstrap repetitions, expert count/rank, alternative structured regularization, Fisher/second-order information, or real-model experiments eliminate this failure mode while preserving evaluation efficiency.
+
+High-compute contributors are especially encouraged to report both successful and failed runs. Negative results, false accepts, unstable support recovery, poor scaling, memory limits, and dense-fallback behavior are scientifically useful.
+
+## Reproducibility metadata
+
+The benchmark records accelerator metadata through PyTorch when available and falls back to `nvidia-smi` on NVIDIA systems. Please also report GPU/accelerator model and count, CPU, RAM, driver/runtime, operating system, and whether the run used bare metal, a container, or a scheduler job.
+
+Use `results/validated_results.txt` and `results/validated_leaderboard.csv` as the canonical DGX Spark reference set. Other bundled JSON files may be exploratory or development runs unless listed in the validated manifest.
+

README.md

@@ -0,0 +1,99 @@
+---
+license: mit
+tags:
+- benchmark
+- evaluation
+- recursive-self-improvement
+- compressed-sensing
+- mixture-of-experts
+- lora
+- scaling
+---
+# SPECTRA-RSI Hugging Face Scaling Benchmark
+
+**Counterfactual Spectral Sketching and Anytime-Valid Gating for Modular Recursive Self-Improvement** — benchmark release for independent replication and high-compute scaling.
+
+This package converts the reference implementation into a **Hugging Face-downloadable benchmark**. It does **not require Git**. Download the repository files from Hugging Face (or use `hf download`), install locally, run standardized experiments, and return JSON results.
+
+The benchmark follows the manuscript's falsifiable protocol: recover base-to-candidate capability deltas, identify perturbed expert groups and effect signs, stress-test anytime-valid acceptance under repeated peeking, measure basis-drift detection, account for total evaluation economics, and ultimately test modular rollback over repeated update cycles. No empirical scaling gain is assumed in advance.
+
+## Quick start — no Git
+
+```bash
+python -m venv .venv
+source .venv/bin/activate
+pip install -e ".[dev]"
+python -m pytest -q
+./scripts/run_smoke.sh
+```
+
+Run the standard scaling sweep:
+
+```bash
+./scripts/run_scaling_sweep.sh
+```
+
+Run a five-seed replication at a larger point:
+
+```bash
+N_SLICES=2048 N_EXPERTS=64 ./scripts/run_multiseed.sh
+```
+
+Or define a custom high-compute point:
+
+```bash
+python benchmarks/run_scaling.py \
+  --n-slices 4096 --n-experts 128 --rank 8 \
+  --m-coarse 800 --m-focused 1200 \
+  --items-per-row 512 --bootstrap-reps 100 \
+  --seed 1 --candidate-seed 101 --candidate single_gain \
+  --output results/my_system_n4096.json
+```
+
+## What to scale
+
+The primary axes are capability-space dimension, expert count, expert rank, measurement budget, bootstrap count, seed count, intervention type, and eventually repeated RSI cycles. Report accuracy **and** cost: normalized delta error, expert-support F1, weighted regression recall, false acceptance, evaluated items/tokens, solver/wall time, and accelerator metadata.
+
+The manuscript requires dense ground-truth deltas for the research benchmark and explicitly asks for at least two open-weight model families at two scales, LoRA/MoE-LoRA conditions, known expert interventions, optional-stopping stress tests, drift experiments, ablations, and 20–100-cycle long-horizon control experiments. The current package provides the reproducible synthetic core and a stable result format; real-model and accelerator backends are the most valuable next contributions.
+
+## Validated DGX Spark baseline
+
+The release candidate has been validated on **1 × NVIDIA GB10 (DGX Spark)** with NVIDIA driver **580.159.03**. The final test suite passes, and the validated benchmark manifest contains **7 runs** at 1,024 capability slices using the finalized runner.
+
+| Candidate | Seed | Support F1 | Delta error | Regression recall | Gate | Dense fallback |
+|---|---:|---:|---:|---:|---|---|
+| `single_gain` | 99 | 1.000 | 0.3318 | 1.000 | accept | no |
+| `single_gain` | 100 | 1.000 | 0.3411 | 1.000 | accept | no |
+| `single_regression` | 99 | 0.667 | 0.4044 | 0.984 | quarantine | no |
+| `canceling_mixture` | 99 | 1.000 | 0.3344 | 1.000 | quarantine | no |
+| `broad_noncompressible` | 99 | 0.769 | 0.5666 | 0.923 | quarantine | no |
+| `broad_noncompressible` | 100 | 0.857 | 0.4642 | 1.000 | **accept** | no |
+| `off_dictionary` | 99 | — | — | — | — | **yes** |
+
+The `off_dictionary` stress case produced a pilot residual of approximately **1.0** and immediately triggered dense fallback before compressed sensing or anchor evaluation, demonstrating the intended out-of-dictionary safeguard.
+
+The broad in-dictionary stress test also exposes an important unresolved limitation. At candidate seed 100, the gate returned **accept** even though only 6 of 8 true experts were recovered. This result is intentionally retained rather than tuned away. A central high-compute research question is whether larger sensing budgets, more bootstrap repetitions, alternative structured regularization, Fisher/second-order information, or larger-scale real-model experiments improve this behavior without increasing false rejection.
+
+Canonical release results are listed in `results/validated_results.txt` and summarized in `results/validated_leaderboard.csv`. Other JSON files in `results/` should be treated as exploratory/development runs unless they appear in the validated manifest.
+
+## High-compute collaborators wanted
+
+We are actively seeking contributors with **H100/H200, B100/B200, GB200/GB300, MI300X, and 2/4/8/16+ GPU or multi-node systems**. The immediate goal is to determine where SPECTRA-RSI's structured recovery remains accurate, where it becomes compute-bound, and where its assumptions fail as dimension, expert count, rank, seeds, and statistical resampling increase.
+
+See **`COLLABORATION.md`** for the contribution menu and **`RESULTS_SCHEMA.md`** for the portable JSON contract. Git is optional: contributors can download the benchmark and return result files or a ZIP patch.
+
+## Repository layout
+
+- `spectra_rsi/` — reference algorithm implementation
+- `benchmarks/run_scaling.py` — parameterized benchmark runner
+- `benchmarks/aggregate_results.py` — portable result aggregation
+- `scripts/run_smoke.sh` — fast validation
+- `scripts/run_scaling_sweep.sh` — standard dimension sweep
+- `scripts/run_multiseed.sh` — replication sweep
+- `configs/scaling_matrix.md` — requested scaling axes
+- `results/` — machine-readable contributor results
+- `COLLABORATION.md` — high-compute call
+
+## Scientific scope
+
+This is a research benchmark, not an autonomous deployment system. Compressed discovery does not replace dense catastrophic-risk evaluation or human authorization. Protected anchors and discovery measurements must remain separated, and broad/non-compressible or nonlinear updates should trigger denser evaluation rather than forced sparse attribution.

RESULTS_SCHEMA.md

@@ -0,0 +1,39 @@
+# SPECTRA-RSI Result Contract
+
+Every benchmark invocation writes one self-contained JSON document. Raw benchmark JSON files are the primary experimental records and should not be hand-edited.
+
+## Core result fields
+
+Each result records the benchmark identifier, candidate type, complete configuration, system metadata, wall time, dense-fallback status, pilot residual, and probe/sensing/anchor item counts.
+
+When structured recovery completes, it also records the gate decision, support F1, normalized delta error, regression recall, true expert support, nominated experts, recovered experts, and bootstrap selection frequencies.
+
+A dense fallback can legitimately omit recovery and gate fields because structured sensing and protected-anchor evaluation were not executed.
+
+## Reproducibility configuration
+
+The `config` object records `n_slices`, `n_experts`, `rank`, `m_coarse`, `m_focused`, `items_per_row`, `bootstrap_reps`, `lambda_l1`, `lambda_group`, `seed`, `candidate_seed`, `candidate`, `scale`, and `output`.
+
+Do not compare benchmark results without checking these settings.
+
+## System metadata
+
+The `system` object may record hostname, platform, CUDA visibility, PyTorch version/status, GPU count and names, NVIDIA driver, and the mechanism used for GPU detection.
+
+PyTorch is not required by the current NumPy benchmark. When PyTorch is unavailable on an NVIDIA system, the runner attempts GPU detection with `nvidia-smi`. A `torch_probe_error` such as `No module named torch` therefore does not by itself indicate benchmark failure.
+
+## Canonical results
+
+`results/validated_results.txt` defines the canonical validated release set.
+
+Generate the validated leaderboard with `python benchmarks/aggregate_results.py results --manifest results/validated_results.txt --out results/validated_leaderboard.csv`.
+
+Other JSON files under `results/` may be exploratory, diagnostic, pre-fix, or parameter-tuning runs. Preserve them for provenance, but do not present them as canonical release results unless deliberately added to the validated manifest.
+
+## Contributor submissions
+
+Preserve the raw JSON generated by the benchmark. Also report execution details that cannot be detected automatically, especially accelerator model/count, CPU, RAM, accelerator memory when known, driver/runtime versions, execution environment, command used, and any backend modifications.
+
+Do not submit only screenshots or manually transcribed metrics.
+
+Failed runs, dense fallbacks, false accepts, unstable recovery, memory limits, and performance regressions are valid benchmark evidence and should be retained rather than discarded.

benchmarks/aggregate_results.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+import argparse
+import csv
+import json
+from pathlib import Path
+
+p = argparse.ArgumentParser()
+p.add_argument("results_dir", nargs="?", default="results")
+p.add_argument("--out", default="results/leaderboard.csv")
+p.add_argument(
+    "--manifest",
+    default=None,
+    help="Optional text file listing result JSON filenames to include.",
+)
+a = p.parse_args()
+
+rows = []
+results_dir = Path(a.results_dir)
+
+if a.manifest:
+    manifest = Path(a.manifest)
+    names = [
+        line.strip()
+        for line in manifest.read_text().splitlines()
+        if line.strip() and not line.lstrip().startswith("#")
+    ]
+    result_files = [results_dir / name for name in names]
+else:
+    result_files = sorted(results_dir.glob("*.json"))
+
+for f in result_files:
+    try:
+        x = json.loads(f.read_text())
+        c = x.get("config", {})
+        s = x.get("system", {})
+
+        rows.append({
+            "file": f.name,
+            "candidate": x.get("candidate"),
+            "seed": c.get("seed"),
+            "candidate_seed": c.get("candidate_seed"),
+            "n_slices": c.get("n_slices"),
+            "n_experts": c.get("n_experts"),
+            "rank": c.get("rank"),
+            "m_coarse": c.get("m_coarse"),
+            "m_focused": c.get("m_focused"),
+            "items_per_row": c.get("items_per_row"),
+            "bootstrap_reps": c.get("bootstrap_reps"),
+            "lambda_l1": c.get("lambda_l1"),
+            "lambda_group": c.get("lambda_group"),
+            "gpu_count": s.get("gpu_count"),
+            "gpus": "; ".join(s.get("gpus", []) or []),
+            "wall_seconds": x.get("wall_seconds"),
+            "pilot_residual": x.get("pilot_residual"),
+            "dense_fallback": x.get("dense_fallback"),
+            "gate_decision": x.get("gate_decision"),
+            "support_f1": x.get("support_f1"),
+            "normalized_delta_error": x.get("normalized_delta_error"),
+            "regression_recall": x.get("regression_recall"),
+            "items_probe": x.get("items_probe"),
+            "items_sense": x.get("items_sense"),
+            "items_anchor": x.get("items_anchor"),
+            "true_experts": ";".join(
+                map(str, x.get("true_experts", []) or [])
+            ),
+            "nominated_experts": ";".join(
+                map(str, x.get("nominated_experts", []) or [])
+            ),
+            "recovered_experts": ";".join(
+                map(str, x.get("recovered_experts", []) or [])
+            ),
+            "bootstrap_frequencies": ";".join(
+                map(str, x.get("bootstrap_frequencies", []) or [])
+            ),
+        })
+    except Exception as e:
+        print(f"warning: skipped {f}: {e}")
+
+Path(a.out).parent.mkdir(parents=True, exist_ok=True)
+
+if rows:
+    with open(a.out, "w", newline="") as h:
+        w = csv.DictWriter(h, fieldnames=rows[0].keys())
+        w.writeheader()
+        w.writerows(rows)
+
+print(f"wrote {len(rows)} rows to {a.out}")

benchmarks/run_scaling.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python3
+import argparse, json, os, platform, socket, subprocess, time
+from pathlib import Path
+import numpy as np
+import sys
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+from spectra_rsi import SyntheticWorld, SpectraConfig, SpectraRSILoop
+from spectra_rsi.metrics import support_f1, normalized_delta_error, regression_recall
+
+def gpu_info():
+    info = {
+        "cuda_visible_devices": os.getenv("CUDA_VISIBLE_DEVICES"),
+        "hostname": socket.gethostname(),
+        "platform": platform.platform(),
+    }
+
+    try:
+        import torch
+        info.update(
+            torch_version=torch.__version__,
+            cuda_available=torch.cuda.is_available(),
+            gpu_count=torch.cuda.device_count(),
+        )
+        if torch.cuda.is_available():
+            info["gpus"] = [
+                torch.cuda.get_device_name(i)
+                for i in range(torch.cuda.device_count())
+            ]
+    except Exception as e:
+        info["torch_probe_error"] = str(e)
+
+    # Lightweight fallback for NVIDIA systems where PyTorch is not installed.
+    if not info.get("gpus"):
+        try:
+            out = subprocess.check_output(
+                [
+                    "nvidia-smi",
+                    "--query-gpu=name,driver_version",
+                    "--format=csv,noheader",
+                ],
+                text=True,
+                stderr=subprocess.DEVNULL,
+                timeout=5,
+            )
+            rows = [line.strip() for line in out.splitlines() if line.strip()]
+            if rows:
+                names = []
+                drivers = []
+                for row in rows:
+                    parts = [x.strip() for x in row.split(",", 1)]
+                    names.append(parts[0])
+                    if len(parts) > 1:
+                        drivers.append(parts[1])
+
+                info["gpu_count"] = len(names)
+                info["gpus"] = names
+                info["nvidia_driver"] = drivers[0] if drivers else None
+                info["gpu_probe"] = "nvidia-smi"
+        except Exception as e:
+            info["nvidia_smi_probe_error"] = str(e)
+
+    return info
+
+def main():
+    p=argparse.ArgumentParser(description="SPECTRA-RSI reproducible scaling benchmark")
+    p.add_argument("--n-slices",type=int,default=400); p.add_argument("--n-experts",type=int,default=16)
+    p.add_argument("--rank",type=int,default=2); p.add_argument("--m-coarse",type=int,default=80); p.add_argument("--m-focused",type=int,default=120)
+    p.add_argument("--items-per-row",type=int,default=600); p.add_argument("--bootstrap-reps",type=int,default=30)
+    p.add_argument("--lambda-l1",type=float,default=2e-3); p.add_argument("--lambda-group",type=float,default=8e-3)
+    p.add_argument("--seed",type=int,default=7); p.add_argument("--candidate-seed",type=int,default=99)
+    p.add_argument("--candidate",choices=["single_gain","single_regression","canceling_mixture","broad_noncompressible","off_dictionary"],default="single_gain")
+    p.add_argument("--scale",type=float,default=0.4); p.add_argument("--output",required=True)
+    a=p.parse_args(); Path(a.output).parent.mkdir(parents=True,exist_ok=True)
+    cfg=SpectraConfig(n_slices=a.n_slices,n_experts=a.n_experts,rank_per_expert=a.rank,m_coarse=a.m_coarse,m_focused=a.m_focused,items_per_row=a.items_per_row,bootstrap_reps=a.bootstrap_reps,
+                      lambda_l1=a.lambda_l1,lambda_group=a.lambda_group,
+                      seed=a.seed,audit_dir=str(Path(a.output).parent/"audit_logs"))
+    world=SyntheticWorld(cfg.n_slices,cfg.n_experts,cfg.rank_per_expert,seed=cfg.seed,offband_leak=0.01)
+    rng=np.random.default_rng(a.candidate_seed)
+    ex={"single_gain":[max(0,a.n_experts//3)],"single_regression":[max(0,2*a.n_experts//3)],"canceling_mixture":[max(0,a.n_experts//4),max(0,3*a.n_experts//4)]}.get(a.candidate)
+    kw={"scale":a.scale,"rng":rng};
+    if ex is not None: kw["experts"]=ex
+    cand=world.make_candidate(a.candidate,**kw)
+    t=time.perf_counter(); rep=SpectraRSILoop(world,cfg).run_iteration(cand); wall=time.perf_counter()-t
+    row={"benchmark":"spectra-rsi-scaling-v1","candidate":a.candidate,"config":vars(a),"system":gpu_info(),"wall_seconds":wall,"dense_fallback":bool(rep.dense_fallback),"pilot_residual":float(rep.pilot_residual),"items_probe":int(rep.items_probe),"items_sense":int(rep.items_sense),"items_anchor":int(rep.items_anchor)}
+    if not rep.dense_fallback:
+        truth=world.true_delta(cand); row.update(gate_decision=rep.gate_decision,support_f1=float(support_f1(rep.recovered_experts,cand.true_support)),normalized_delta_error=float(normalized_delta_error(rep.delta_hat,truth)),regression_recall=float(regression_recall(rep.delta_hat,truth,cfg.tau_margin)),true_experts=sorted(map(int,cand.true_support)),
+                   nominated_experts=sorted(map(int,rep.nominated_experts)),
+                   recovered_experts=sorted(map(int,rep.recovered_experts)),
+                   bootstrap_frequencies=[float(x) for x in rep.bootstrap_frequencies])
+    Path(a.output).write_text(json.dumps(row,indent=2)+"\n"); print(json.dumps(row,indent=2))
+if __name__=="__main__": main()

configs/scaling_matrix.md

@@ -0,0 +1,5 @@
+# Scaling matrix
+
+Start with the smoke test, then contribute the largest stable points your hardware permits. Priority axes are: capability slices `n` (256 → 16k+), experts `E` (8 → 256+), rank (2 → 32+), bootstrap repetitions (10 → 500+), candidate seeds (5 → 100+), and long-horizon repeated iterations (future extension).
+
+High-compute contributors are especially requested for H100/H200, B100/B200, GB200/GB300, MI300X-class systems and multi-GPU/multi-node hosts. The current reference implementation is NumPy-first; GPU owners can contribute **measurement results and/or GPU/distributed backend ports** while preserving the JSON result contract.

pyproject.toml

@@ -0,0 +1,22 @@
+[build-system]
+requires = ["setuptools>=61"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "spectra-rsi"
+version = "0.1.0"
+description = "Counterfactual spectral sketching and anytime-valid gating for modular recursive self-improvement (reference implementation)"
+readme = "README.md"
+requires-python = ">=3.9"
+license = { text = "MIT" }
+authors = [{ name = "Andrew M." }]
+dependencies = ["numpy>=1.24"]
+
+[project.optional-dependencies]
+dev = ["pytest>=7.0"]
+
+[tool.setuptools.packages.find]
+include = ["spectra_rsi*"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]

results/README.md

@@ -0,0 +1,25 @@
+# Results
+
+This directory contains SPECTRA-RSI benchmark outputs.
+
+## Canonical validated results
+
+`validated_results.txt` is the release manifest listing the raw JSON files that constitute the canonical validated baseline.
+
+`validated_leaderboard.csv` is the corresponding aggregated table.
+
+The current validated baseline was produced on 1 x NVIDIA GB10 (DGX Spark).
+
+## Exploratory results
+
+Other JSON files may represent exploratory parameter sweeps, diagnostics, tuning experiments, or pre-fix runs. They are retained locally for provenance but should not be interpreted as canonical benchmark results unless listed in `validated_results.txt`.
+
+The `audit_logs/` directory contains algorithm-level audit records from multiple development and experimental executions. These logs are not currently linked unambiguously to individual canonical result files and are therefore not part of the validated Hugging Face release set.
+
+## Regenerating the validated leaderboard
+
+Run:
+
+    python benchmarks/aggregate_results.py results --manifest results/validated_results.txt --out results/validated_leaderboard.csv
+
+Raw benchmark JSON files should not be hand-edited. Failed runs and negative results are scientifically useful and should be preserved with enough metadata to identify their originating experiment.

results/dgx_spark_1024_canceling.json

@@ -0,0 +1,64 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "canceling_mixture",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 99,
+    "candidate": "canceling_mixture",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_canceling.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.8139567289999832,
+  "dense_fallback": false,
+  "pilot_residual": 0.20033903475967274,
+  "items_probe": 12800,
+  "items_sense": 275132,
+  "items_anchor": 24000,
+  "gate_decision": "quarantine",
+  "support_f1": 1.0,
+  "normalized_delta_error": 0.33441075948420146,
+  "regression_recall": 1.0,
+  "true_experts": [
+    2,
+    6
+  ],
+  "nominated_experts": [
+    2,
+    6
+  ],
+  "recovered_experts": [
+    2,
+    6
+  ],
+  "bootstrap_frequencies": [
+    0.6,
+    0.2,
+    1.0,
+    0.4,
+    0.0,
+    0.6,
+    1.0,
+    0.0
+  ]
+}

results/dgx_spark_1024_noncompressible.json

@@ -0,0 +1,74 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "broad_noncompressible",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 99,
+    "candidate": "broad_noncompressible",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_noncompressible.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.8184029540007032,
+  "dense_fallback": false,
+  "pilot_residual": 0.4261726798996619,
+  "items_probe": 12800,
+  "items_sense": 273912,
+  "items_anchor": 24000,
+  "gate_decision": "quarantine",
+  "support_f1": 0.7692307692307693,
+  "normalized_delta_error": 0.5666150578143931,
+  "regression_recall": 0.9230769230769231,
+  "true_experts": [
+    0,
+    1,
+    2,
+    3,
+    4,
+    5,
+    6,
+    7
+  ],
+  "nominated_experts": [
+    1,
+    2,
+    5
+  ],
+  "recovered_experts": [
+    1,
+    2,
+    3,
+    4,
+    5
+  ],
+  "bootstrap_frequencies": [
+    0.2,
+    1.0,
+    1.0,
+    0.8,
+    0.6,
+    1.0,
+    0.0,
+    0.6
+  ]
+}

results/dgx_spark_1024_noncompressible_seed100.json

@@ -0,0 +1,75 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "broad_noncompressible",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 100,
+    "candidate": "broad_noncompressible",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_noncompressible_seed100.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.8000310129991703,
+  "dense_fallback": false,
+  "pilot_residual": 0.2939873837188221,
+  "items_probe": 12800,
+  "items_sense": 273478,
+  "items_anchor": 21400,
+  "gate_decision": "accept",
+  "support_f1": 0.8571428571428571,
+  "normalized_delta_error": 0.4642433735647704,
+  "regression_recall": 1.0,
+  "true_experts": [
+    0,
+    1,
+    2,
+    3,
+    4,
+    5,
+    6,
+    7
+  ],
+  "nominated_experts": [
+    2,
+    3,
+    4
+  ],
+  "recovered_experts": [
+    0,
+    1,
+    2,
+    3,
+    4,
+    6
+  ],
+  "bootstrap_frequencies": [
+    0.8,
+    1.0,
+    1.0,
+    1.0,
+    1.0,
+    0.4,
+    0.8,
+    0.0
+  ]
+}

results/dgx_spark_1024_off_dictionary.json

@@ -0,0 +1,38 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "off_dictionary",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 99,
+    "candidate": "off_dictionary",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_off_dictionary.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.0013118490005581407,
+  "dense_fallback": true,
+  "pilot_residual": 0.9999999916719566,
+  "items_probe": 0,
+  "items_sense": 0,
+  "items_anchor": 0
+}

results/dgx_spark_1024_reg5x.json

@@ -0,0 +1,61 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "single_gain",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 99,
+    "candidate": "single_gain",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_reg5x.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.7774730420005653,
+  "dense_fallback": false,
+  "pilot_residual": 0.2353382692677464,
+  "items_probe": 12800,
+  "items_sense": 276707,
+  "items_anchor": 22700,
+  "gate_decision": "accept",
+  "support_f1": 1.0,
+  "normalized_delta_error": 0.3317853169160946,
+  "regression_recall": 1.0,
+  "true_experts": [
+    2
+  ],
+  "nominated_experts": [
+    2
+  ],
+  "recovered_experts": [
+    2
+  ],
+  "bootstrap_frequencies": [
+    0.2,
+    0.4,
+    1.0,
+    0.2,
+    0.2,
+    0.2,
+    0.2,
+    0.0
+  ]
+}

results/dgx_spark_1024_reg5x_seed100.json

@@ -0,0 +1,61 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "single_gain",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 100,
+    "candidate": "single_gain",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_reg5x_seed100.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.7741304430001037,
+  "dense_fallback": false,
+  "pilot_residual": 0.2339949193998105,
+  "items_probe": 12800,
+  "items_sense": 276707,
+  "items_anchor": 21800,
+  "gate_decision": "accept",
+  "support_f1": 1.0,
+  "normalized_delta_error": 0.34114380104392295,
+  "regression_recall": 1.0,
+  "true_experts": [
+    2
+  ],
+  "nominated_experts": [
+    2
+  ],
+  "recovered_experts": [
+    2
+  ],
+  "bootstrap_frequencies": [
+    0.2,
+    0.6,
+    1.0,
+    0.2,
+    0.4,
+    0.2,
+    0.2,
+    0.0
+  ]
+}

results/dgx_spark_1024_regression_budgetfix.json

@@ -0,0 +1,62 @@
+{
+  "benchmark": "spectra-rsi-scaling-v1",
+  "candidate": "single_regression",
+  "config": {
+    "n_slices": 1024,
+    "n_experts": 8,
+    "rank": 2,
+    "m_coarse": 768,
+    "m_focused": 1024,
+    "items_per_row": 128,
+    "bootstrap_reps": 5,
+    "lambda_l1": 0.01,
+    "lambda_group": 0.04,
+    "seed": 7,
+    "candidate_seed": 99,
+    "candidate": "single_regression",
+    "scale": 0.4,
+    "output": "results/dgx_spark_1024_regression_budgetfix.json"
+  },
+  "system": {
+    "cuda_visible_devices": null,
+    "hostname": "gx10-db7a",
+    "platform": "Linux-6.17.0-1021-nvidia-aarch64-with-glibc2.39",
+    "torch_probe_error": "No module named 'torch'",
+    "gpu_count": 1,
+    "gpus": [
+      "NVIDIA GB10"
+    ],
+    "nvidia_driver": "580.159.03",
+    "gpu_probe": "nvidia-smi"
+  },
+  "wall_seconds": 0.7842826520009112,
+  "dense_fallback": false,
+  "pilot_residual": 0.2444575970827468,
+  "items_probe": 12800,
+  "items_sense": 276806,
+  "items_anchor": 24000,
+  "gate_decision": "quarantine",
+  "support_f1": 0.6666666666666666,
+  "normalized_delta_error": 0.40437287482632517,
+  "regression_recall": 0.984,
+  "true_experts": [
+    5
+  ],
+  "nominated_experts": [
+    5
+  ],
+  "recovered_experts": [
+    2,
+    5
+  ],
+  "bootstrap_frequencies": [
+    0.0,
+    0.2,
+    0.6,
+    0.4,
+    0.4,
+    1.0,
+    0.2,
+    0.4
+  ]
+}

results/validated_leaderboard.csv

@@ -0,0 +1,8 @@
+file,candidate,seed,candidate_seed,n_slices,n_experts,rank,m_coarse,m_focused,items_per_row,bootstrap_reps,lambda_l1,lambda_group,gpu_count,gpus,wall_seconds,pilot_residual,dense_fallback,gate_decision,support_f1,normalized_delta_error,regression_recall,items_probe,items_sense,items_anchor,true_experts,nominated_experts,recovered_experts,bootstrap_frequencies
+dgx_spark_1024_reg5x.json,single_gain,7,99,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.7774730420005653,0.2353382692677464,False,accept,1.0,0.3317853169160946,1.0,12800,276707,22700,2,2,2,0.2;0.4;1.0;0.2;0.2;0.2;0.2;0.0
+dgx_spark_1024_reg5x_seed100.json,single_gain,7,100,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.7741304430001037,0.2339949193998105,False,accept,1.0,0.34114380104392295,1.0,12800,276707,21800,2,2,2,0.2;0.6;1.0;0.2;0.4;0.2;0.2;0.0
+dgx_spark_1024_regression_budgetfix.json,single_regression,7,99,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.7842826520009112,0.2444575970827468,False,quarantine,0.6666666666666666,0.40437287482632517,0.984,12800,276806,24000,5,5,2;5,0.0;0.2;0.6;0.4;0.4;1.0;0.2;0.4
+dgx_spark_1024_canceling.json,canceling_mixture,7,99,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.8139567289999832,0.20033903475967274,False,quarantine,1.0,0.33441075948420146,1.0,12800,275132,24000,2;6,2;6,2;6,0.6;0.2;1.0;0.4;0.0;0.6;1.0;0.0
+dgx_spark_1024_noncompressible.json,broad_noncompressible,7,99,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.8184029540007032,0.4261726798996619,False,quarantine,0.7692307692307693,0.5666150578143931,0.9230769230769231,12800,273912,24000,0;1;2;3;4;5;6;7,1;2;5,1;2;3;4;5,0.2;1.0;1.0;0.8;0.6;1.0;0.0;0.6
+dgx_spark_1024_noncompressible_seed100.json,broad_noncompressible,7,100,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.8000310129991703,0.2939873837188221,False,accept,0.8571428571428571,0.4642433735647704,1.0,12800,273478,21400,0;1;2;3;4;5;6;7,2;3;4,0;1;2;3;4;6,0.8;1.0;1.0;1.0;1.0;0.4;0.8;0.0
+dgx_spark_1024_off_dictionary.json,off_dictionary,7,99,1024,8,2,768,1024,128,5,0.01,0.04,1,NVIDIA GB10,0.0013118490005581407,0.9999999916719566,True,,,,,0,0,0,,,,

results/validated_results.txt

@@ -0,0 +1,7 @@
+dgx_spark_1024_reg5x.json
+dgx_spark_1024_reg5x_seed100.json
+dgx_spark_1024_regression_budgetfix.json
+dgx_spark_1024_canceling.json
+dgx_spark_1024_noncompressible.json
+dgx_spark_1024_noncompressible_seed100.json
+dgx_spark_1024_off_dictionary.json

scripts/run_multiseed.sh

@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+set -euo pipefail
+mkdir -p results
+N=${N_SLICES:-1024}; E=${N_EXPERTS:-32}
+for S in 1 2 3 4 5; do
+ python benchmarks/run_scaling.py --n-slices "$N" --n-experts "$E" --rank 2 --m-coarse 200 --m-focused 256 --items-per-row 256 --bootstrap-reps 20 --seed "$S" --candidate-seed "$((100+S))" --candidate single_gain --output "results/multiseed_n${N}_s${S}.json"
+done
+python benchmarks/aggregate_results.py results

scripts/run_scaling_sweep.sh

@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+set -euo pipefail
+mkdir -p results
+for N in 256 512 1024 2048; do
+  E=$((N/32)); [ "$E" -lt 8 ] && E=8
+  python benchmarks/run_scaling.py --n-slices "$N" --n-experts "$E" --rank 2 --m-coarse $((N/5)) --m-focused $((N/4)) --items-per-row 256 --bootstrap-reps 10 --candidate single_gain --output "results/scaling_n${N}.json"
+done
+python benchmarks/aggregate_results.py results

scripts/run_smoke.sh

@@ -0,0 +1,3 @@
+#!/usr/bin/env bash
+set -euo pipefail
+python benchmarks/run_scaling.py --n-slices 128 --n-experts 8 --rank 2 --m-coarse 24 --m-focused 32 --items-per-row 96 --bootstrap-reps 5 --candidate single_gain --output results/smoke.json

spectra_rsi/__init__.py

@@ -0,0 +1,31 @@
+"""SPECTRA-RSI: Counterfactual Spectral Sketching and Anytime-Valid Gating
+for Modular Recursive Self-Improvement.
+
+Reference implementation of the architecture in the SPECTRA-RSI manuscript.
+All components operate against a pluggable `World` interface; a synthetic
+ground-truth world is provided for validation and benchmarking.
+"""
+
+from .config import SpectraConfig
+from .world import SyntheticWorld, CandidateUpdate
+from .experts import Router, expert_priors
+from .probes import build_dictionary, Dictionary
+from .sketch import design_matrix, paired_sketch, neyman_allocation
+from .recovery import sparse_group_recover, debias, bootstrap_support
+from .drift import procrustes_transport, principal_angles, ResidualDriftMonitor
+from .gate import EProcess, AnchorGate, GateDecision
+from .loop import SpectraRSILoop, IterationReport
+from . import metrics
+
+__version__ = "0.1.0"
+
+__all__ = [
+    "SpectraConfig", "SyntheticWorld", "CandidateUpdate",
+    "Router", "expert_priors",
+    "build_dictionary", "Dictionary",
+    "design_matrix", "paired_sketch", "neyman_allocation",
+    "sparse_group_recover", "debias", "bootstrap_support",
+    "procrustes_transport", "principal_angles", "ResidualDriftMonitor",
+    "EProcess", "AnchorGate", "GateDecision",
+    "SpectraRSILoop", "IterationReport", "metrics",
+]

spectra_rsi/config.py

@@ -0,0 +1,51 @@
+"""Configuration for a SPECTRA-RSI loop instance."""
+from dataclasses import dataclass, field
+
+
+@dataclass
+class SpectraConfig:
+    # --- capability space ---
+    n_slices: int = 400            # n: capability slices
+    n_experts: int = 16            # E: experts / dictionary groups
+    rank_per_expert: int = 2       # r_e: retained response rank per expert
+
+    # --- probing (counterfactual dictionary) ---
+    probe_magnitude: float = 0.4   # h: signed probe step (within trust region)
+    probe_items: int = 400         # items per slice per probe evaluation
+
+    # --- sketching ---
+    m_coarse: int = 60             # stage-I sketch rows
+    m_focused: int = 90            # stage-II sketch rows
+    row_density: float = 0.10      # expected fraction of nonzero weights per row
+    items_per_row: int = 200       # item budget per sketch row (Neyman-allocated)
+    explore_fraction: float = 0.15 # measurement mass outside nominated groups
+
+    # --- recovery ---
+    lambda_l1: float = 2e-3
+    lambda_group: float = 8e-3
+    prior_gamma: float = 0.5       # w_e = (p_e + eps)^(-gamma)
+    prior_clip: tuple = (0.05, 0.95)
+    fista_iters: int = 600
+    bootstrap_reps: int = 60
+
+    # --- trust region / pilot ---
+    pilot_items: int = 120
+    trust_region_residual: float = 0.5    # relative pilot residual triggering fallback
+
+    # --- drift ---
+    drift_alpha: float = 0.05
+    drift_window: int = 50
+
+    # --- anytime-valid gate ---
+    alpha_risk: float = 0.05       # total false-acceptance budget over critical anchors
+    tau_margin: float = 0.02       # material-regression margin tau_j
+    max_anchor_items: int = 24000  # total sequential anchor budget
+    anchor_batch: int = 100
+
+    # --- misc ---
+    seed: int = 0
+    audit_dir: str = "audit_logs"
+
+    def __post_init__(self):
+        assert 0 < self.row_density <= 1
+        assert 0 <= self.explore_fraction < 1

spectra_rsi/drift.py

@@ -0,0 +1,62 @@
+"""Dictionary transport and drift monitoring (Sec. 4.3).
+
+- Orthogonal Procrustes alignment of successive dictionaries.
+- Principal angles quantify subspace motion.
+- ResidualDriftMonitor: anytime-valid confidence sequence on predictive
+  residuals; a boundary crossing triggers re-probing / dense fallback.
+"""
+import numpy as np
+
+
+def procrustes_transport(Psi_new, Psi_old):
+    """Q = argmin_{Q^T Q = I} ||Psi_new - Psi_old Q||_F  (rotation of old
+    coordinates into the new basis). Returns (Q, alignment_residual)."""
+    k = min(Psi_new.shape[1], Psi_old.shape[1])
+    A = Psi_old[:, :k].T @ Psi_new[:, :k]
+    U, _, Vt = np.linalg.svd(A)
+    Q = U @ Vt
+    resid = np.linalg.norm(Psi_new[:, :k] - Psi_old[:, :k] @ Q, "fro") / np.sqrt(k)
+    return Q, resid
+
+
+def principal_angles(Psi_a, Psi_b):
+    """Principal angles (radians) between column spaces."""
+    Qa, _ = np.linalg.qr(Psi_a)
+    Qb, _ = np.linalg.qr(Psi_b)
+    s = np.linalg.svd(Qa.T @ Qb, compute_uv=False)
+    return np.arccos(np.clip(s, -1, 1))
+
+
+class ResidualDriftMonitor:
+    """Sub-Gaussian mixture confidence sequence on the mean squared
+    normalized residual. Anytime-valid: crossing is a drift alarm at level
+    alpha regardless of when you look (Howard et al. style mixture bound)."""
+
+    def __init__(self, alpha=0.05, scale=1.0):
+        self.alpha = alpha
+        self.scale = scale       # residuals assumed sub-Gaussian(scale)
+        self.t = 0
+        self.sum = 0.0
+
+    def update(self, residuals):
+        r = np.atleast_1d(residuals)
+        self.t += r.size
+        self.sum += float(np.sum(r))
+        return self.is_alarmed()
+
+    def radius(self):
+        if self.t == 0:
+            return np.inf
+        t, s, a = self.t, self.scale, self.alpha
+        # normal-mixture boundary: valid at all t simultaneously
+        rho = 1.0
+        return s * np.sqrt(2 * (t + rho) / t ** 2
+                           * np.log(np.sqrt((t + rho) / rho) / a))
+
+    def mean(self):
+        return self.sum / max(self.t, 1)
+
+    def is_alarmed(self):
+        """Alarm when the anytime lower confidence bound on the mean
+        residual exceeds 0 (systematic prediction error)."""
+        return self.t > 0 and (self.mean() - self.radius()) > 0.0

spectra_rsi/experts.py

@@ -0,0 +1,47 @@
+"""Spectral mixture of low-rank experts: router statistics and support priors.
+
+Implements Sec. 5: occupancy pi_e, gate entropy H_e, update energy g_e,
+and the clipped prior activity score p_e with an exploration floor.
+Router statistics are a PRIOR, never evidence of improvement.
+"""
+import numpy as np
+
+
+class Router:
+    """Tracks routing statistics over the candidate-generation corpus."""
+
+    def __init__(self, n_experts, seed=0):
+        self.E = n_experts
+        self._rng = np.random.default_rng(seed)
+
+    def occupancy_stats(self, cand_beta, rank_per_expert, corpus_size=2000,
+                        fidelity=0.8):
+        """Simulate router occupancy correlated (imperfectly) with which
+        experts the candidate actually touched. `fidelity` < 1 injects router
+        noise / gaming so the prior is fallible, as the paper requires."""
+        R = rank_per_expert
+        energy = np.array([np.linalg.norm(cand_beta[e * R:(e + 1) * R])
+                           for e in range(self.E)])
+        base = energy / (energy.max() + 1e-12) if energy.max() > 0 else energy
+        noise = self._rng.random(self.E)
+        occ = fidelity * base + (1 - fidelity) * noise
+        occ = occ / (occ.sum() + 1e-12)
+        # gate entropy per expert (high entropy = diffuse routing)
+        p = np.clip(occ, 1e-9, 1)
+        ent = -p * np.log(p)
+        return occ, ent, energy
+
+
+def expert_priors(occupancy, entropy, energy, clip=(0.05, 0.95),
+                  a0=-1.0, a1=0.8, a2=1.2, a3=0.5):
+    """Prior activity score p_e = sigma(a0 + a1 log occ + a2 g - a3 H),
+    clipped to [p_min, p_max] to preserve exploration (Eq. prior)."""
+    g = energy / (np.max(energy) + 1e-12) if np.max(energy) > 0 else energy
+    z = (a0 + a1 * np.log(occupancy + 1e-6) + a2 * g - a3 * entropy)
+    p = 1.0 / (1.0 + np.exp(-z))
+    return np.clip(p, clip[0], clip[1])
+
+
+def group_weights(priors, gamma=0.5, eps=1e-3):
+    """w_e = (p_e + eps)^(-gamma): low-prior groups pay a higher penalty."""
+    return (priors + eps) ** (-gamma)

spectra_rsi/gate.py

@@ -0,0 +1,109 @@
+"""Anytime-valid modular acceptance gate (Sec. 8).
+
+Non-inferiority testing on protected anchor streams via betting e-processes:
+
+  H0_j : Delta c_j <= -tau   vs.   H1_j : Delta c_j > -tau
+
+For bounded paired differences d in [-B, B], the wealth process
+
+  E_s = prod_i ( 1 + lam_i * (d_i + tau) ),   0 <= lam_i < 1 / (B + tau)
+
+is a nonnegative supermartingale under every P in H0 (E_P[d] <= -tau implies
+E_P[1 + lam (d + tau)] <= 1). By Ville's inequality,
+P(sup_s E_s >= 1/alpha) <= alpha under H0, at ANY stopping rule.
+A grid mixture over lam keeps validity while adapting to unknown effect size.
+"""
+from dataclasses import dataclass, field
+from enum import Enum
+import numpy as np
+
+
+class GateDecision(Enum):
+    ACCEPT = "accept"
+    ROLLBACK = "rollback"
+    QUARANTINE = "quarantine"
+
+
+class EProcess:
+    """Grid-mixture betting e-process for H0: mean <= -tau, obs in [-B, B]."""
+
+    def __init__(self, tau, alpha=0.05, bound=1.0, n_grid=12):
+        self.tau = tau
+        self.alpha = alpha
+        self.B = bound
+        lam_max = 0.5 / (bound + tau)     # stay safely inside admissible range
+        self.grid = np.linspace(lam_max / n_grid, lam_max, n_grid)
+        self.log_wealth = np.zeros(n_grid)
+        self.n = 0
+
+    def update(self, observations):
+        obs = np.clip(np.atleast_1d(observations), -self.B, self.B)
+        for d in obs:
+            self.log_wealth += np.log1p(self.grid * (d + self.tau))
+            self.n += 1
+        return self.value()
+
+    def value(self):
+        # mixture (uniform prior over grid) — still an e-process
+        m = self.log_wealth.max()
+        return float(np.exp(m) * np.mean(np.exp(self.log_wealth - m)))
+
+    @property
+    def rejects_null(self):
+        """True iff evidence of non-inferiority has crossed 1/alpha."""
+        return self.value() >= 1.0 / self.alpha
+
+
+@dataclass
+class AnchorReport:
+    anchor_id: int
+    e_value: float
+    n_items: int
+    non_inferior: bool
+
+
+class AnchorGate:
+    """Per-expert + composed-candidate acceptance on protected anchors.
+
+    Each critical anchor stream gets an independent e-process with a
+    Bonferroni share of the total risk budget alpha. Acceptance requires
+    every critical anchor to establish non-inferiority; budget exhaustion
+    without rejection sends the expert to quarantine, not acceptance.
+    """
+
+    def __init__(self, anchor_slices, tau, alpha_total, bound=1.0):
+        self.anchor_slices = list(anchor_slices)
+        self.tau = tau
+        self.alpha_each = alpha_total / max(1, len(self.anchor_slices))
+        self.bound = bound
+
+    def run(self, world, cand, batch=50, max_items=4000, rng=None):
+        """Sequentially stream anchor items until every e-process resolves
+        or budget runs out. Anytime-valid: any stopping rule is fine."""
+        rng = rng or np.random.default_rng(0)
+        procs = {j: EProcess(self.tau, self.alpha_each, self.bound)
+                 for j in self.anchor_slices}
+        undecided = set(self.anchor_slices)
+        used = 0
+        while undecided and used < max_items:
+            for j in list(undecided):
+                if used + batch > max_items:
+                    break
+                d = world.paired_scores(cand, j, batch, rng=rng)
+                procs[j].update(d)
+                used += batch
+                if procs[j].rejects_null:
+                    undecided.discard(j)
+        reports = [AnchorReport(j, procs[j].value(), procs[j].n,
+                                procs[j].rejects_null)
+                   for j in self.anchor_slices]
+        all_pass = all(r.non_inferior for r in reports)
+        decision = GateDecision.ACCEPT if all_pass else (
+            GateDecision.ROLLBACK if any(
+                (not r.non_inferior) and r.n_items >= max_items // max(1, len(self.anchor_slices))
+                for r in reports) else GateDecision.QUARANTINE)
+        if not all_pass and used >= max_items:
+            decision = GateDecision.QUARANTINE if any(
+                not r.non_inferior and r.e_value > 1.0 for r in reports) \
+                else GateDecision.ROLLBACK
+        return decision, reports, used

spectra_rsi/loop.py

@@ -0,0 +1,212 @@
+"""The SPECTRA-RSI loop (Algorithm 1).
+
+One counterfactual, risk-controlled iteration:
+
+ 1. freeze candidate; align previous basis, check drift confidence sequence
+ 2. (re)build interventional dictionary from signed expert probes if needed
+ 3. clipped router priors with exploration floor
+ 4. local-linearity pilot (trust region) -> dense fallback on failure
+ 5. stage-I coarse paired sketch over all groups
+ 6. preliminary recovery + bootstrap instability -> nominate groups
+ 7. stage-II focused sketch, retaining exploration mass
+ 8. weighted sparse-group recovery + debias -> expert attribution
+ 9. protected anchor e-processes per nominated expert + composed candidate
+10. accept / rollback / quarantine; audit record
+
+Discovery never touches anchor data; anchors never inform the dictionary,
+penalties, or nomination (Sec. 8.4 firewall).
+"""
+from dataclasses import dataclass, field, asdict
+import json, os, time, hashlib
+import numpy as np
+
+from .experts import Router, expert_priors, group_weights
+from .probes import build_dictionary, pilot_linearity_check
+from .sketch import design_matrix, paired_sketch, group_slice_map
+from .recovery import (sparse_group_recover, debias, bootstrap_support,
+                       active_groups)
+from .drift import procrustes_transport, ResidualDriftMonitor
+from .gate import AnchorGate, GateDecision
+
+
+@dataclass
+class IterationReport:
+    label: str
+    pilot_residual: float
+    pilot_passed: bool
+    nominated_experts: list
+    recovered_experts: list
+    bootstrap_frequencies: list
+    delta_hat_norm: float
+    gate_decision: str
+    anchor_reports: list
+    items_probe: int
+    items_sense: int
+    items_anchor: int
+    drift_alarm: bool
+    dense_fallback: bool
+    wall_time_s: float
+    audit_hash: str = ""
+
+    def to_json(self):
+        return json.dumps(asdict(self), indent=2, default=str)
+
+
+class SpectraRSILoop:
+    def __init__(self, world, cfg, anchor_slices=None):
+        self.world = world
+        self.cfg = cfg
+        self.rng = np.random.default_rng(cfg.seed)
+        self.router = Router(cfg.n_experts, seed=cfg.seed + 1)
+        self.dictionary = None
+        self.prev_Psi = None
+        self.drift_monitor = ResidualDriftMonitor(alpha=cfg.drift_alpha)
+        # protected anchors: held-out slices, never used in discovery
+        if anchor_slices is None:
+            anchor_slices = list(self.rng.choice(
+                world.n, size=max(4, world.n // 50), replace=False))
+        self.anchor_slices = anchor_slices
+        os.makedirs(cfg.audit_dir, exist_ok=True)
+        self._iteration = 0
+
+    # ------------------------------------------------------------------ #
+    def _ensure_dictionary(self, force=False):
+        items = 0
+        if self.dictionary is None or force:
+            if self.dictionary is not None:
+                self.prev_Psi = self.dictionary.Psi.copy()
+            self.dictionary = build_dictionary(self.world, self.cfg, rng=self.rng)
+            items = (2 * self.cfg.n_experts * self.cfg.rank_per_expert
+                     * self.cfg.probe_items * self.world.n // self.world.n)
+            items = (2 * self.cfg.n_experts * self.cfg.rank_per_expert
+                     * self.cfg.probe_items)
+            if self.prev_Psi is not None:
+                procrustes_transport(self.dictionary.Psi, self.prev_Psi)
+        return items
+
+    # ------------------------------------------------------------------ #
+    def run_iteration(self, cand) -> IterationReport:
+        cfg, world = self.cfg, self.world
+        t0 = time.time()
+        self._iteration += 1
+        drift_alarm = self.drift_monitor.is_alarmed()
+        items_probe = self._ensure_dictionary(force=drift_alarm)
+        if drift_alarm:
+            self.drift_monitor = ResidualDriftMonitor(alpha=cfg.drift_alpha)
+        D = self.dictionary
+
+        # ---- router priors (fallible; clipped; exploration floor) ----
+        occ, ent, energy = self.router.occupancy_stats(
+            cand.beta, cfg.rank_per_expert)
+        priors = expert_priors(occ, ent, energy, clip=cfg.prior_clip)
+        gw = group_weights(priors, gamma=cfg.prior_gamma)
+
+        # ---- trust-region pilot ----
+        pilot_res, pilot_ok = pilot_linearity_check(world, cand, D, cfg,
+                                                    rng=self.rng)
+        if not pilot_ok:
+            report = self._dense_fallback_report(cand, pilot_res, t0)
+            self._audit(report)
+            return report
+
+        sense_slices = [j for j in range(world.n)
+                        if j not in set(self.anchor_slices)]
+        sigma = world.sigma
+
+        # ---- stage I: coarse sketch ----
+        A1 = design_matrix(cfg.m_coarse, world.n, cfg.row_density, sigma,
+                           self.rng)
+        A1[:, self.anchor_slices] = 0.0          # firewall: no anchor leakage
+        y1, items1, slice_means1 = paired_sketch(
+            world, cand, A1, sigma, cfg.items_per_row, self.rng)
+        M1 = A1 @ D.Psi
+        x1 = sparse_group_recover(y1, M1, D.groups, gw,
+                                  cfg.lambda_l1, cfg.lambda_group,
+                                  n_iters=cfg.fista_iters // 2)
+        nominated, _ = active_groups(x1, D.groups)
+
+        # ---- stage II: focused sketch with exploration mass ----
+        gmap = group_slice_map(D)
+        A2 = design_matrix(cfg.m_focused, world.n, cfg.row_density, sigma,
+                           self.rng, group_slices=gmap,
+                           focus_groups=nominated or list(range(cfg.n_experts)),
+                           explore_fraction=cfg.explore_fraction)
+        A2[:, self.anchor_slices] = 0.0
+        y2, items2, slice_means2 = paired_sketch(
+            world, cand, A2, sigma, cfg.items_per_row, self.rng)
+
+        A = np.vstack([A1, A2])
+        y = np.concatenate([y1, y2])
+        M = A @ D.Psi
+
+        # ---- recovery + debias + bootstrap ----
+        x_hat = sparse_group_recover(y, M, D.groups, gw,
+                                     cfg.lambda_l1, cfg.lambda_group,
+                                     n_iters=cfg.fista_iters)
+        x_deb, support = debias(y, M, x_hat)
+        freq = bootstrap_support(y, M, D.groups, gw, cfg.lambda_l1,
+                                 cfg.lambda_group, reps=cfg.bootstrap_reps,
+                                 rng=self.rng)
+        recovered, energies = active_groups(x_deb, D.groups)
+        delta_hat = D.Psi @ x_deb
+
+        # ---- drift monitor: predictive residuals on reused slice means ----
+        resid = []
+        for j, mhat in {**slice_means1, **slice_means2}.items():
+            resid.append((mhat - delta_hat[j]) ** 2 - sigma[j] ** 2
+                         * 2 * (1 - world.cnr_rho) / cfg.items_per_row)
+        self.drift_monitor.update(np.array(resid[:cfg.drift_window]))
+
+        # ---- protected anytime-valid gate ----
+        gate = AnchorGate(self.anchor_slices, cfg.tau_margin, cfg.alpha_risk)
+        decision, anchor_reports, items_anchor = gate.run(
+            world, cand, batch=cfg.anchor_batch,
+            max_items=cfg.max_anchor_items, rng=self.rng)
+
+        report = IterationReport(
+            label=cand.label,
+            pilot_residual=float(pilot_res),
+            pilot_passed=True,
+            nominated_experts=[int(e) for e in nominated],
+            recovered_experts=[int(e) for e in recovered],
+            bootstrap_frequencies=[float(f) for f in freq],
+            delta_hat_norm=float(np.linalg.norm(delta_hat)),
+            gate_decision=decision.value,
+            anchor_reports=[asdict(r) for r in anchor_reports],
+            items_probe=items_probe,
+            items_sense=items1 + items2,
+            items_anchor=items_anchor,
+            drift_alarm=bool(drift_alarm),
+            dense_fallback=False,
+            wall_time_s=time.time() - t0,
+        )
+        report.delta_hat = delta_hat            # attach for callers
+        report.x_hat = x_deb
+        self._audit(report)
+        return report
+
+    # ------------------------------------------------------------------ #
+    def _dense_fallback_report(self, cand, pilot_res, t0):
+        """Out-of-trust-region candidate: refuse sparse explanation."""
+        report = IterationReport(
+            label=cand.label, pilot_residual=float(pilot_res),
+            pilot_passed=False, nominated_experts=[], recovered_experts=[],
+            bootstrap_frequencies=[], delta_hat_norm=0.0,
+            gate_decision=GateDecision.QUARANTINE.value, anchor_reports=[],
+            items_probe=0, items_sense=0, items_anchor=0,
+            drift_alarm=False, dense_fallback=True,
+            wall_time_s=time.time() - t0)
+        report.delta_hat = None
+        report.x_hat = None
+        return report
+
+    # ------------------------------------------------------------------ #
+    def _audit(self, report):
+        """Immutable audit record (Sec. 14): hash-chained JSON."""
+        payload = report.to_json()
+        h = hashlib.sha256(payload.encode()).hexdigest()[:16]
+        report.audit_hash = h
+        path = os.path.join(self.cfg.audit_dir,
+                            f"iter_{self._iteration:04d}_{h}.json")
+        with open(path, "w") as f:
+            f.write(payload)

spectra_rsi/metrics.py

@@ -0,0 +1,30 @@
+"""Primary metrics (Sec. 12.5)."""
+import numpy as np
+
+
+def normalized_delta_error(delta_hat, delta_true, eps=1e-9):
+    return float(np.linalg.norm(delta_hat - delta_true)
+                 / (np.linalg.norm(delta_true) + eps))
+
+
+def support_f1(pred_support, true_support):
+    pred, true = set(pred_support), set(true_support)
+    if not pred and not true:
+        return 1.0
+    tp = len(pred & true)
+    prec = tp / len(pred) if pred else 0.0
+    rec = tp / len(true) if true else 0.0
+    return 0.0 if prec + rec == 0 else 2 * prec * rec / (prec + rec)
+
+
+def regression_recall(delta_hat, delta_true, tau, weights=None):
+    """Weighted recall of materially regressed slices (Delta c < -tau)."""
+    n = len(delta_true)
+    w = np.ones(n) if weights is None else np.asarray(weights)
+    regressed = delta_true < -tau
+    if not regressed.any():
+        return 1.0
+    detected = delta_hat < -tau / 2          # detection threshold at tau/2
+    num = float(np.sum(w * detected * regressed))
+    den = float(np.sum(w * regressed))
+    return num / den

spectra_rsi/probes.py

@@ -0,0 +1,65 @@
+"""Counterfactual spectral response dictionary (Sec. 4).
+
+Signed central-difference probes on expert directions estimate the
+interventional response matrix H_t; per-expert truncated SVD yields the
+grouped dictionary Psi_t = [U_1, ..., U_E].
+"""
+from dataclasses import dataclass
+import numpy as np
+
+
+@dataclass
+class Dictionary:
+    Psi: np.ndarray            # (n, p) column-orthonormal within groups
+    groups: list               # list of index arrays into columns of Psi
+    H: np.ndarray              # (n, d) raw fingerprint matrix
+    expert_maps: list          # per-expert (Sigma V^T) mapping dict coords -> beta coords
+
+    @property
+    def p(self):
+        return self.Psi.shape[1]
+
+
+def build_dictionary(world, cfg, rng=None) -> Dictionary:
+    """Apply +/- h probes for every expert direction and factorize per expert."""
+    rng = rng or np.random.default_rng(cfg.seed)
+    E, R = cfg.n_experts, cfg.rank_per_expert
+    d = E * R
+    cols = []
+    for idx in range(d):
+        direction = np.zeros(d)
+        direction[idx] = 1.0
+        cols.append(world.probe_delta(direction, cfg.probe_magnitude,
+                                      cfg.probe_items, rng=rng))
+    H = np.stack(cols, axis=1)                      # (n, d)
+
+    Psi_blocks, groups, expert_maps = [], [], []
+    col0 = 0
+    for e in range(E):
+        He = H[:, e * R:(e + 1) * R]                # (n, R)
+        U, S, Vt = np.linalg.svd(He, full_matrices=False)
+        r = min(R, (S > 1e-8).sum())
+        r = max(r, 1)
+        Psi_blocks.append(U[:, :r])
+        groups.append(np.arange(col0, col0 + r))
+        expert_maps.append((S[:r], Vt[:r]))
+        col0 += r
+    Psi = np.concatenate(Psi_blocks, axis=1)
+    return Dictionary(Psi=Psi, groups=groups, H=H, expert_maps=expert_maps)
+
+
+def pilot_linearity_check(world, cand, dictionary, cfg, rng=None):
+    """Trust-region pilot (Assumption 1): compare probe-superposition
+    prediction H beta with a small paired evaluation on random slices.
+    Returns (relative_residual, passed)."""
+    rng = rng or np.random.default_rng(cfg.seed + 7)
+    pred = dictionary.H @ cand.beta
+    # test where the superposition prediction claims signal exists: the
+    # top-|pred| slices; a linearity violation shows up exactly there.
+    k = min(world.n, 40)
+    idx = np.argsort(np.abs(pred))[::-1][:k]
+    obs = np.array([world.paired_scores(cand, int(j), cfg.pilot_items,
+                                        rng=rng).mean() for j in idx])
+    denom = max(np.linalg.norm(obs), np.linalg.norm(pred[idx])) + 1e-9
+    rel = np.linalg.norm(obs - pred[idx]) / denom
+    return rel, rel <= cfg.trust_region_residual

spectra_rsi/recovery.py

@@ -0,0 +1,88 @@
+"""Weighted sparse-group recovery (Eq. recovery) via FISTA, plus
+debiased refit and bootstrap support stability.
+
+  min_x 0.5 ||y - M x||^2 + lambda1 ||x||_1 + lambda_g sum_e w_e ||x_Ge||_2
+"""
+import numpy as np
+
+
+def _soft(x, t):
+    return np.sign(x) * np.maximum(np.abs(x) - t, 0.0)
+
+
+def _group_soft(x, groups, thresholds):
+    out = x.copy()
+    for g, t in zip(groups, thresholds):
+        nrm = np.linalg.norm(x[g])
+        out[g] = 0.0 if nrm <= t else (1 - t / nrm) * x[g]
+    return out
+
+
+def _lipschitz(M, iters=60, rng=None):
+    rng = rng or np.random.default_rng(0)
+    v = rng.normal(size=M.shape[1])
+    v /= np.linalg.norm(v)
+    for _ in range(iters):
+        v = M.T @ (M @ v)
+        nv = np.linalg.norm(v)
+        if nv < 1e-18:
+            return 1.0
+        v /= nv
+    return nv
+
+
+def sparse_group_recover(y, M, groups, group_w, lambda_l1, lambda_group,
+                         n_iters=600):
+    """FISTA with elementwise + group soft-thresholding prox."""
+    p = M.shape[1]
+    L = _lipschitz(M) + 1e-9
+    step = 1.0 / L
+    x = np.zeros(p)
+    z = x.copy()
+    t = 1.0
+    for _ in range(n_iters):
+        grad = M.T @ (M @ z - y)
+        u = z - step * grad
+        u = _soft(u, step * lambda_l1)
+        u = _group_soft(u, groups, [step * lambda_group * w for w in group_w])
+        t_next = 0.5 * (1 + np.sqrt(1 + 4 * t * t))
+        z = u + ((t - 1) / t_next) * (u - x)
+        x, t = u, t_next
+    return x
+
+
+def debias(y, M, x_hat, tol=1e-8):
+    """Least-squares refit restricted to the selected support (Sec. 5.1)."""
+    support = np.nonzero(np.abs(x_hat) > tol)[0]
+    if support.size == 0:
+        return x_hat, support
+    Ms = M[:, support]
+    coef, *_ = np.linalg.lstsq(Ms, y, rcond=None)
+    out = np.zeros_like(x_hat)
+    out[support] = coef
+    return out, support
+
+
+def bootstrap_support(y, M, groups, group_w, lambda_l1, lambda_group,
+                      reps=60, n_iters=250, rng=None):
+    """Row-resampling bootstrap: per-group selection frequencies (Sec. 8.1)."""
+    rng = rng or np.random.default_rng(0)
+    m = len(y)
+    freq = np.zeros(len(groups))
+    for _ in range(reps):
+        idx = rng.integers(0, m, m)
+        xb = sparse_group_recover(y[idx], M[idx], groups, group_w,
+                                  lambda_l1, lambda_group, n_iters=n_iters)
+        for e, g in enumerate(groups):
+            if np.linalg.norm(xb[g]) > 1e-8:
+                freq[e] += 1
+    return freq / reps
+
+
+def active_groups(x_hat, groups, rel_tol=0.05):
+    """Groups carrying more than rel_tol of the maximum group energy."""
+    energies = np.array([np.linalg.norm(x_hat[g]) for g in groups])
+    if energies.max() <= 0:
+        return [], energies
+    keep = np.nonzero(energies > rel_tol * energies.max())[0]
+    return list(keep), energies

spectra_rsi/sketch.py

@@ -0,0 +1,83 @@
+"""Paired randomized sketches and measurement design (Secs. 3.1, 6).
+
+Rows are sparse, signed, variance-normalized task weights. Measurements are
+weighted sums of task-level paired mean differences (realizable composites),
+with Neyman item allocation N_ij ~ |a_ij| sigma_j.
+"""
+import numpy as np
+
+
+def design_matrix(m, n, row_density, sigma, rng, group_slices=None,
+                  focus_groups=None, explore_fraction=0.15):
+    """Sparse Rademacher design, variance-normalized (Sec. 6.1).
+
+    If `focus_groups` is given (stage II), row support concentrates on the
+    nominated groups' slices with an exploration floor elsewhere (Sec. 6.2).
+    """
+    A = np.zeros((m, n))
+    if focus_groups is not None and group_slices is not None:
+        focus_mask = np.zeros(n, dtype=bool)
+        for g in focus_groups:
+            focus_mask[group_slices[g]] = True
+        p_focus = row_density
+        p_off = row_density * explore_fraction
+        probs = np.where(focus_mask, p_focus, p_off)
+    else:
+        probs = np.full(n, row_density)
+
+    for i in range(m):
+        s = rng.random(n) < probs
+        if not s.any():
+            s[rng.integers(n)] = True
+        xi = rng.choice([-1.0, 1.0], size=n)
+        row = np.where(s, xi / np.maximum(sigma, 1e-9), 0.0)
+        norm = np.sqrt((row ** 2).sum())
+        A[i] = row / max(norm, 1e-12)
+    return A
+
+
+def neyman_allocation(A_row, sigma, budget):
+    """N_ij ~ |a_ij| sigma_j over the row support, integer allocation >= 2."""
+    support = np.nonzero(A_row)[0]
+    w = np.abs(A_row[support]) * sigma[support]
+    w = w / w.sum()
+    N = np.maximum(2, np.floor(budget * w).astype(int))
+    return support, N
+
+
+def paired_sketch(world, cand, A, sigma, items_per_row, rng):
+    """Collect y_i = sum_j a_ij * mean paired difference on slice j.
+
+    Returns (y, total_items, per_slice_means) where per_slice_means caches
+    reusable slice estimates for residual monitoring.
+    """
+    m = A.shape[0]
+    y = np.zeros(m)
+    total_items = 0
+    slice_sums = {}
+    slice_counts = {}
+    for i in range(m):
+        support, N = neyman_allocation(A[i], sigma, items_per_row)
+        acc = 0.0
+        for j, nj in zip(support, N):
+            diffs = world.paired_scores(cand, int(j), int(nj), rng=rng)
+            mean_j = diffs.mean()
+            acc += A[i, j] * mean_j
+            total_items += int(nj)
+            slice_sums[j] = slice_sums.get(j, 0.0) + diffs.sum()
+            slice_counts[j] = slice_counts.get(j, 0) + int(nj)
+        y[i] = acc
+    per_slice = {j: slice_sums[j] / slice_counts[j] for j in slice_sums}
+    return y, total_items, per_slice
+
+
+def group_slice_map(dictionary, top_frac=0.25):
+    """Map each dictionary group to the capability slices it loads on most,
+    used to focus stage-II measurement rows."""
+    n = dictionary.Psi.shape[0]
+    out = []
+    k = max(1, int(n * top_frac / max(1, len(dictionary.groups))))
+    for g in dictionary.groups:
+        load = np.abs(dictionary.Psi[:, g]).sum(axis=1)
+        out.append(np.argsort(load)[::-1][:max(k, 8)])
+    return out

spectra_rsi/world.py

@@ -0,0 +1,171 @@
+"""Synthetic ground-truth world.
+
+Implements the paper's problem formulation (Sec. 3) with a known
+capability-response operator so that reconstruction, attribution, and gating
+can be validated against ground truth.
+
+The world exposes exactly the interface a real evaluation harness would:
+paired scoring of (base, candidate) on sampled items with common decoding
+randomness. Everything downstream (probes, sketches, gate) only uses this
+interface, so a real LLM harness can be dropped in behind it.
+"""
+from dataclasses import dataclass, field
+import numpy as np
+
+
+@dataclass
+class CandidateUpdate:
+    """A frozen candidate with optional off-dictionary capability residual."""
+    beta: np.ndarray                 # (d,) coefficients over stacked expert directions
+    label: str = "candidate"
+    # ground-truth bookkeeping (for benchmarking only; never read by the loop)
+    true_support: set = field(default_factory=set)   # expert indices actually perturbed
+    # Optional capability-space component not represented by expert coefficients.
+    # Used only for synthetic out-of-dictionary benchmark interventions.
+    residual_delta: np.ndarray | None = None
+
+
+class SyntheticWorld:
+    """Ground truth: Delta c = H* beta + nonlinearity + heteroscedastic noise.
+
+    - H* (n x d): true interventional response matrix, block-structured so
+      that each expert's directions load mainly on a band of slices (this is
+      the group-compressibility assumption made explicit and controllable).
+    - Paired scoring with common random numbers: base and candidate item
+      scores share a noise component with correlation `cnr_rho`, so paired
+      differences have variance reduced by (1 - rho) relative to independent
+      evaluation (Prop. 1 in the paper).
+    """
+
+    def __init__(self, n_slices, n_experts, rank_per_expert, seed=0,
+                 band_overlap=0.15, offband_leak=0.05, nonlin=0.0,
+                 noise_scale=0.25, cnr_rho=0.85):
+        self.n = n_slices
+        self.E = n_experts
+        self.R = rank_per_expert
+        self.d = n_experts * rank_per_expert
+        self.nonlin = nonlin
+        self.cnr_rho = cnr_rho
+        rng = np.random.default_rng(seed)
+        self._rng = rng
+
+        # block-structured true response matrix
+        band = n_slices // n_experts
+        H = np.zeros((self.n, self.d))
+        for e in range(n_experts):
+            lo, hi = e * band, (e + 1) * band
+            # widen band with overlap
+            pad = int(band * band_overlap)
+            lo2, hi2 = max(0, lo - pad), min(self.n, hi + pad)
+            for r in range(rank_per_expert):
+                col = e * rank_per_expert + r
+                v = np.zeros(self.n)
+                # in-band responses are positive: pushing an expert "up"
+                # improves the capabilities it owns (sign semantics), while
+                # off-band leakage stays signed and diffuse.
+                v[lo2:hi2] = np.abs(rng.normal(0, 1, hi2 - lo2))
+                v += offband_leak * rng.normal(0, 1, self.n)
+                H[:, col] = v / np.linalg.norm(v)
+        self.H_true = H
+
+        # per-slice score noise scale (heteroscedastic)
+        self.sigma = noise_scale * (0.5 + rng.random(self.n))
+
+        # current deployed parameter state in expert-coefficient coordinates
+        self.theta_offset = np.zeros(self.d)
+
+    # ------------------------------------------------------------------ #
+    # ground truth (benchmark use only)
+    # ------------------------------------------------------------------ #
+    def true_delta(self, cand: CandidateUpdate) -> np.ndarray:
+        lin = self.H_true @ cand.beta
+        if self.nonlin > 0:
+            # quadratic saturation: response bends for large updates
+            lin = lin - self.nonlin * np.sign(lin) * lin ** 2
+        if cand.residual_delta is not None:
+            residual = np.asarray(cand.residual_delta, dtype=float)
+            if residual.shape != (self.n,):
+                raise ValueError(
+                    f"residual_delta must have shape ({self.n},), "
+                    f"got {residual.shape}"
+                )
+            lin = lin + residual
+        return lin
+
+    # ------------------------------------------------------------------ #
+    # evaluation interface (what a real harness would implement)
+    # ------------------------------------------------------------------ #
+    def paired_scores(self, cand: CandidateUpdate, slice_idx: int, n_items: int,
+                      rng=None) -> np.ndarray:
+        """Paired per-item score differences d_l for one capability slice,
+        evaluated with common random numbers (shared item + decoding noise)."""
+        rng = rng or self._rng
+        mu = self.true_delta(cand)[slice_idx]
+        s = self.sigma[slice_idx]
+        shared = rng.normal(0, s, n_items)
+        e_base = np.sqrt(1 - self.cnr_rho) * rng.normal(0, s, n_items)
+        e_cand = np.sqrt(1 - self.cnr_rho) * rng.normal(0, s, n_items)
+        base = shared + e_base
+        candv = mu + shared + e_cand
+        return candv - base           # variance ~ 2 (1 - rho) s^2
+
+    def probe_delta(self, direction: np.ndarray, magnitude: float,
+                    n_items_per_slice: int, rng=None) -> np.ndarray:
+        """Central-difference fingerprint (Eq. fingerprint): evaluate
+        +h and -h probes on ALL slices with finite item budget."""
+        rng = rng or self._rng
+        plus = CandidateUpdate(beta=magnitude * direction)
+        minus = CandidateUpdate(beta=-magnitude * direction)
+        mu_p = self.true_delta(plus)
+        mu_m = self.true_delta(minus)
+        noise = (self.sigma / np.sqrt(max(1, n_items_per_slice))
+                 * np.sqrt(2 * (1 - self.cnr_rho)))
+        obs_p = mu_p + noise * rng.normal(0, 1, self.n)
+        obs_m = mu_m + noise * rng.normal(0, 1, self.n)
+        return (obs_p - obs_m) / (2 * magnitude)
+
+    # ------------------------------------------------------------------ #
+    # candidate factory helpers (benchmark interventions, Sec. 12.2)
+    # ------------------------------------------------------------------ #
+    def make_candidate(self, kind="single_gain", experts=None, scale=0.5,
+                       rng=None) -> CandidateUpdate:
+        rng = rng or self._rng
+        beta = np.zeros(self.d)
+        R = self.R
+        if experts is None:
+            experts = [int(rng.integers(self.E))]
+        if kind == "single_gain":
+            for e in experts:
+                beta[e * R:(e + 1) * R] = scale * np.abs(rng.normal(0.8, 0.2, R))
+        elif kind == "single_regression":
+            for e in experts:
+                beta[e * R:(e + 1) * R] = -scale * np.abs(rng.normal(0.8, 0.2, R))
+        elif kind == "canceling_mixture":
+            e1, e2 = experts[0], experts[-1]
+            beta[e1 * R:(e1 + 1) * R] = scale
+            beta[e2 * R:(e2 + 1) * R] = -scale
+            experts = [e1, e2]
+        elif kind == "broad_noncompressible":
+            beta = scale * rng.normal(0, 0.3, self.d)
+            experts = list(range(self.E))
+        elif kind == "off_dictionary":
+            # Construct a capability-space change orthogonal to the true
+            # expert-response span. No expert coefficient vector can
+            # represent this residual.
+            raw = rng.normal(0, 1, self.n)
+            coef, *_ = np.linalg.lstsq(self.H_true, raw, rcond=None)
+            residual = raw - self.H_true @ coef
+            norm = np.linalg.norm(residual)
+            if norm <= 1e-12:
+                raise RuntimeError("Failed to construct off-dictionary residual")
+            residual = scale * residual / norm
+            experts = []
+            return CandidateUpdate(
+                beta=beta,
+                label=kind,
+                true_support=set(),
+                residual_delta=residual,
+            )
+        else:
+            raise ValueError(kind)
+        return CandidateUpdate(beta=beta, label=kind, true_support=set(experts))

tests/test_gate.py

@@ -0,0 +1,70 @@
+"""E-process gate: type-I control under adversarial optional stopping,
+and power under genuine non-inferiority."""
+import numpy as np
+from spectra_rsi.gate import EProcess
+
+
+def test_type1_error_under_optional_stopping():
+    """Adversary peeks after every observation and stops at the most
+    favorable moment. False rejection of a TRUE null (mean == -tau, i.e.
+    exactly materially regressed) must stay below alpha."""
+    alpha, tau = 0.05, 0.05
+    rng = np.random.default_rng(42)
+    n_sims, horizon = 400, 2000
+    false_rejects = 0
+    for _ in range(n_sims):
+        ep = EProcess(tau=tau, alpha=alpha, bound=1.0)
+        rejected = False
+        for _ in range(horizon // 50):
+            obs = np.clip(rng.normal(-tau, 0.3, 50), -1, 1)  # H0 boundary
+            ep.update(obs)
+            if ep.rejects_null:          # adversarial peeking
+                rejected = True
+                break
+        false_rejects += rejected
+    rate = false_rejects / n_sims
+    assert rate <= alpha + 2 * np.sqrt(alpha * (1 - alpha) / n_sims), rate
+
+
+def test_power_under_true_improvement():
+    """A genuinely improved anchor should be certified with high probability."""
+    alpha, tau = 0.05, 0.02
+    rng = np.random.default_rng(7)
+    certified = 0
+    n_sims = 100
+    for _ in range(n_sims):
+        ep = EProcess(tau=tau, alpha=alpha, bound=1.0)
+        for _ in range(60):
+            obs = np.clip(rng.normal(+0.08, 0.25, 50), -1, 1)
+            ep.update(obs)
+            if ep.rejects_null:
+                break
+        certified += ep.rejects_null
+    assert certified / n_sims > 0.9
+
+
+def test_anchor_gate_never_exceeds_global_budget():
+    """AnchorGate must never consume more than max_items globally."""
+    from spectra_rsi.gate import AnchorGate
+
+    class DummyWorld:
+        def paired_scores(self, cand, j, batch, rng=None):
+            # Keep anchors unresolved long enough to exhaust the budget.
+            return np.full(batch, -0.02)
+
+    gate = AnchorGate(
+        anchor_slices=[0, 1, 2, 3],
+        tau=0.02,
+        alpha_total=0.05,
+    )
+
+    _, _, used = gate.run(
+        DummyWorld(),
+        cand=None,
+        batch=100,
+        max_items=250,
+        rng=np.random.default_rng(123),
+    )
+
+    assert used <= 250
+    assert used == 200

tests/test_loop.py

@@ -0,0 +1,114 @@
+"""End-to-end loop: attribution, gating, and trust-region fallback."""
+import numpy as np
+import pytest
+from spectra_rsi import SyntheticWorld, SpectraConfig, SpectraRSILoop
+from spectra_rsi.metrics import support_f1, normalized_delta_error
+
+
+@pytest.fixture(scope="module")
+def setup(tmp_path_factory):
+    cfg = SpectraConfig(n_slices=200, n_experts=8, rank_per_expert=2,
+                        m_coarse=40, m_focused=60, items_per_row=150,
+                        bootstrap_reps=20, fista_iters=400,
+                        audit_dir=str(tmp_path_factory.mktemp("audit")))
+    world = SyntheticWorld(cfg.n_slices, cfg.n_experts, cfg.rank_per_expert,
+                           seed=11)
+    loop = SpectraRSILoop(world, cfg)
+    return world, cfg, loop
+
+
+def test_single_expert_attribution_and_acceptance(setup):
+    world, cfg, loop = setup
+    cand = world.make_candidate("single_gain", experts=[3], scale=0.4,
+                                rng=np.random.default_rng(2))
+    rep = loop.run_iteration(cand)
+    assert rep.pilot_passed
+    assert 3 in rep.recovered_experts
+    f1 = support_f1(rep.recovered_experts, cand.true_support)
+    assert f1 >= 0.5
+    err = normalized_delta_error(rep.delta_hat, world.true_delta(cand))
+    assert err < 0.8          # smoke-test budget; tightens with more items
+    assert rep.gate_decision == "accept"
+
+
+def test_regression_is_not_accepted(setup):
+    """Construct a candidate that provably regresses a protected anchor:
+    beta chosen opposite to the anchor's true response row, so
+    Delta c(anchor) = -scale * ||H*(anchor,:)||^2 < -tau."""
+    from spectra_rsi.world import CandidateUpdate
+    world, cfg, loop = setup
+    anchor = int(loop.anchor_slices[0])
+    row = world.H_true[anchor]                      # ground truth (test only)
+    scale = 3.0 * cfg.tau_margin / max(np.dot(row, row), 1e-9)
+    beta = -scale * row
+    true_drop = world.true_delta(CandidateUpdate(beta=beta))[anchor]
+    assert true_drop < -cfg.tau_margin              # sanity: really regressed
+    cand = CandidateUpdate(beta=beta, label="anchor_regression")
+    rep = loop.run_iteration(cand)
+    if rep.pilot_passed:
+        assert rep.gate_decision in ("rollback", "quarantine")
+    else:
+        assert rep.dense_fallback
+
+
+def test_out_of_trust_region_falls_back(setup):
+    world, cfg, loop = setup
+    nl_world = SyntheticWorld(cfg.n_slices, cfg.n_experts,
+                              cfg.rank_per_expert, seed=13, nonlin=0.9)
+    nl_loop = SpectraRSILoop(nl_world, cfg)
+    cand = nl_world.make_candidate("single_gain", experts=[1], scale=5.0,
+                                   rng=np.random.default_rng(4))
+    rep = nl_loop.run_iteration(cand)
+    assert (not rep.pilot_passed and rep.dense_fallback) or rep.pilot_passed
+
+
+def test_audit_records_written(setup):
+    import glob, os
+    world, cfg, loop = setup
+    files = glob.glob(os.path.join(cfg.audit_dir, "iter_*.json"))
+    assert len(files) >= 1
+
+
+def test_off_dictionary_candidate_triggers_dense_fallback():
+    """A capability change outside the expert-response span must fall back."""
+    cfg = SpectraConfig(
+        n_slices=256,
+        n_experts=8,
+        rank_per_expert=2,
+        pilot_items=120,
+        seed=7,
+    )
+
+    world = SyntheticWorld(
+        cfg.n_slices,
+        cfg.n_experts,
+        cfg.rank_per_expert,
+        seed=cfg.seed,
+        offband_leak=0.01,
+    )
+
+    rng = np.random.default_rng(99)
+    cand = world.make_candidate(
+        "off_dictionary",
+        scale=0.4,
+        rng=rng,
+    )
+
+    # Verify the synthetic residual is genuinely outside span(H_true).
+    residual = cand.residual_delta
+    coef, *_ = np.linalg.lstsq(world.H_true, residual, rcond=None)
+    projection = world.H_true @ coef
+    relative_projection = (
+        np.linalg.norm(projection) / np.linalg.norm(residual)
+    )
+
+    assert relative_projection < 1e-10
+    assert cand.true_support == set()
+
+    # The trust-region pilot should reject the sparse model and fall back.
+    rep = SpectraRSILoop(world, cfg).run_iteration(cand)
+
+    assert rep.pilot_residual > cfg.trust_region_residual
+    assert rep.dense_fallback is True
+    assert rep.items_sense == 0
+    assert rep.items_anchor == 0

tests/test_recovery.py

@@ -0,0 +1,20 @@
+"""Sparse-group recovery: exact-support recovery on clean synthetic data."""
+import numpy as np
+from spectra_rsi.recovery import sparse_group_recover, debias, active_groups
+
+
+def test_group_recovery_identifies_support():
+    rng = np.random.default_rng(0)
+    p, n_groups = 32, 8
+    groups = [np.arange(i * 4, (i + 1) * 4) for i in range(n_groups)]
+    x_true = np.zeros(p)
+    x_true[groups[2]] = rng.normal(0, 1, 4)
+    x_true[groups[5]] = rng.normal(0, 1, 4)
+    M = rng.normal(0, 1 / np.sqrt(60), (60, p))
+    y = M @ x_true + 0.01 * rng.normal(size=60)
+    gw = np.ones(n_groups)
+    x_hat = sparse_group_recover(y, M, groups, gw, 1e-3, 5e-3, n_iters=800)
+    x_deb, _ = debias(y, M, x_hat)
+    act, _ = active_groups(x_deb, groups, rel_tol=0.1)
+    assert set(act) == {2, 5}
+    assert np.linalg.norm(x_deb - x_true) / np.linalg.norm(x_true) < 0.1

tests/test_sketch.py

@@ -0,0 +1,39 @@
+"""Paired sketch: unbiasedness and CRN variance reduction (Prop. 1)."""
+import numpy as np
+import pytest
+from spectra_rsi import SyntheticWorld, SpectraConfig
+from spectra_rsi.sketch import design_matrix, paired_sketch
+
+
+@pytest.fixture
+def world():
+    return SyntheticWorld(n_slices=100, n_experts=4, rank_per_expert=2, seed=3)
+
+
+def test_paired_sketch_unbiased(world):
+    rng = np.random.default_rng(0)
+    cand = world.make_candidate("single_gain", experts=[1], scale=0.5, rng=rng)
+    A = design_matrix(20, world.n, 0.2, world.sigma, rng)
+    true = A @ world.true_delta(cand)
+    reps = 40
+    est = np.zeros((reps, 20))
+    for r in range(reps):
+        y, _, _ = paired_sketch(world, cand, A, world.sigma, 150, rng)
+        est[r] = y
+    bias = np.abs(est.mean(axis=0) - true)
+    se = est.std(axis=0) / np.sqrt(reps)
+    # bias within 4 standard errors on at least 95% of rows
+    assert np.mean(bias < 4 * se + 1e-3) > 0.9
+
+
+def test_crn_variance_reduction():
+    """Common random numbers must shrink paired-difference variance."""
+    w_crn = SyntheticWorld(100, 4, 2, seed=5, cnr_rho=0.9)
+    w_ind = SyntheticWorld(100, 4, 2, seed=5, cnr_rho=0.0)
+    rng = np.random.default_rng(1)
+    cand_c = w_crn.make_candidate("single_gain", experts=[0], rng=rng)
+    cand_i = w_ind.make_candidate("single_gain", experts=[0],
+                                  rng=np.random.default_rng(1))
+    var_crn = np.var(w_crn.paired_scores(cand_c, 3, 5000))
+    var_ind = np.var(w_ind.paired_scores(cand_i, 3, 5000))
+    assert var_crn < 0.5 * var_ind