Run independent Semi-knockoffs power audit

code/claim5_scope_audit.py

@@ -0,0 +1,164 @@
+"""Fresh CPU audit of the Semi-knockoffs Section 5.1 comparison.
+
+This is an independent execution, not a read of the authors' CSV outputs.  It
+uses the paper's adjacent-support design (Sigma_ij=0.6**abs(i-j), first
+quarter of beta nonzero, gradient boosting) and a separately implemented HRT
+holdout.  A second exact construction exercises the masked-correlation
+five-permutation assertion.
+"""
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+from sklearn.base import clone
+from sklearn.ensemble import GradientBoostingRegressor
+from scipy.stats import wilcoxon
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+from skocore import ar1_design
+
+
+ALPHA = 0.05
+OUT = {"alpha": ALPHA, "design": "independent_cpu_claim5_scope_audit"}
+
+
+def source_design(n, p, seed):
+    rng = np.random.default_rng(seed)
+    # ar1_design is the exact Gaussian construction for Sigma_ij=rho^|i-j|.
+    X = ar1_design(n, p, 0.6, rng)
+    beta = np.zeros(p)
+    beta[: p // 4] = rng.uniform(1.0, 2.0, size=p // 4)
+    y = np.einsum("ij,j->i", X, beta) + rng.normal(size=n)
+    return X, y, rng
+
+
+def ridge_predict(Z_train, target, Z_eval, alpha=1.0):
+    """Closed-form finite Ridge prediction using einsum, avoiding BLAS matmul."""
+    gram = np.einsum("ij,ik->jk", Z_train, Z_train)
+    rhs = np.einsum("ij,i->j", Z_train, target)
+    coef = np.linalg.solve(gram + alpha * np.eye(Z_train.shape[1]), rhs)
+    return np.einsum("ij,j->i", Z_eval, coef)
+
+
+def semi_pvalue(X, y, j, model, rng, seed, n_perm=1):
+    """Standalone Algorithm-1 p-value with closed-form Ridge nuisances."""
+    Xmj = np.delete(X, j, axis=1)
+    xj = X[:, j]
+    nu = ridge_predict(Xmj, xj, Xmj)
+    rho = ridge_predict(np.column_stack([Xmj, y]), xj,
+                        np.column_stack([Xmj, y]))
+    e1 = xj - nu
+    e2 = xj - rho
+    diffs = []
+    for _ in range(n_perm):
+        x1 = X.copy()
+        x2 = X.copy()
+        x1[:, j] = nu + e1[rng.permutation(len(y))]
+        x2[:, j] = rho + e2[rng.permutation(len(y))]
+        diffs.append((model.predict(x1) - y) ** 2 -
+                     (model.predict(x2) - y) ** 2)
+    d = np.concatenate(diffs)
+    if np.allclose(d, 0.0):
+        return 1.0
+    return float(wilcoxon(d, alternative="greater", zero_method="zsplit").pvalue)
+
+
+def hrt_pvalue(X, y, j, seed, permutations=100):
+    """One HRT p-value with one train/test split and a conditional ridge draw."""
+    rng = np.random.default_rng(seed)
+    order = rng.permutation(len(y))
+    n_train = len(y) // 2
+    train, test = order[:n_train], order[n_train:]
+    model = GradientBoostingRegressor(random_state=seed).fit(X[train], y[train])
+    xm = np.delete(X, j, axis=1)
+    mu = ridge_predict(xm[train], X[train, j], xm[test])
+    residual = X[train, j] - ridge_predict(xm[train], X[train, j], xm[train])
+    observed = np.mean((model.predict(X[test]) - y[test]) ** 2)
+    count = 0
+    for _ in range(permutations):
+        xp = X[test].copy()
+        xp[:, j] = mu + rng.choice(residual, size=len(test), replace=True)
+        count += float(np.mean((model.predict(xp) - y[test]) ** 2) <= observed)
+    return (1.0 + count) / (1.0 + permutations)
+
+
+def adjacent_grid(ns=(200, 300, 400), reps=15):
+    rows = []
+    for n in ns:
+        sko_alt = []
+        sko_null = []
+        hrt_alt = []
+        hrt_null = []
+        for r in range(reps):
+            seed = 18000 + 100 * n + r
+            X, y, rng = source_design(n, 50, seed)
+            model = GradientBoostingRegressor(random_state=seed).fit(X, y)
+            # j=0 is in the first-quarter adjacent support; j=30 is null.
+            sko_alt.append(semi_pvalue(X, y, 0, model, rng, seed=seed) <= ALPHA)
+            sko_null.append(semi_pvalue(X, y, 30, model, rng, seed=seed + 1) <= ALPHA)
+            hrt_alt.append(hrt_pvalue(X, y, 0, seed + 2) <= ALPHA)
+            hrt_null.append(hrt_pvalue(X, y, 30, seed + 3) <= ALPHA)
+        row = {
+            "n": n,
+            "reps": reps,
+            "sko_power": float(np.mean(sko_alt)),
+            "sko_type_I": float(np.mean(sko_null)),
+            "hrt_power": float(np.mean(hrt_alt)),
+            "hrt_type_I": float(np.mean(hrt_null)),
+        }
+        row["power_gap"] = row["sko_power"] - row["hrt_power"]
+        rows.append(row)
+        print("adjacent", row, flush=True)
+    return rows
+
+
+def masked_grid(reps=15, amplitudes=(0.25, 0.5, 1.0)):
+    rows = []
+    for amplitude in amplitudes:
+        cells = []
+        for r in range(reps):
+            seed = 24000 + 100 * int(amplitude * 100) + r
+            rng = np.random.default_rng(seed)
+            X = ar1_design(300, 50, 0.6, rng)
+            signal = 25
+            null = signal - 1
+            y = amplitude * X[:, signal] + 0.5 * rng.normal(size=300)
+            # Replace the adjacent coordinate by a correlated but label-null copy.
+            X[:, null] = X[:, signal] + 0.5 * rng.normal(size=300)
+            model = GradientBoostingRegressor(random_state=seed).fit(X, y)
+            one = semi_pvalue(X, y, signal, model, rng, seed=seed, n_perm=1) <= ALPHA
+            five = semi_pvalue(X, y, signal, model, rng, seed=seed + 10,
+                               n_perm=5) <= ALPHA
+            null_p = semi_pvalue(X, y, null, model, rng, seed=seed + 30,
+                                 n_perm=5)
+            cells.append({
+                "one_reject": bool(one),
+                "five_reject": bool(five),
+                "five_null_reject": bool(null_p <= ALPHA),
+            })
+        row = {
+            "signal_amplitude": amplitude,
+            "reps": reps,
+            "single_power": float(np.mean([r["one_reject"] for r in cells])),
+            "five_power": float(np.mean([r["five_reject"] for r in cells])),
+            "five_null_type_I": float(np.mean([r["five_null_reject"] for r in cells])),
+        }
+        row["power_gain"] = row["five_power"] - row["single_power"]
+        rows.append(row)
+        print("masked", row, flush=True)
+    out = {"reps": reps, "amplitudes": list(amplitudes), "cells": rows}
+    print("masked", out, flush=True)
+    return out
+
+
+def main():
+    OUT["adjacent"] = adjacent_grid()
+    OUT["masked"] = masked_grid()
+    out = Path("outputs/claim5_scope_audit.json")
+    out.write_text(json.dumps(OUT, indent=2) + "\n")
+    print("saved", out, flush=True)
+
+
+if __name__ == "__main__":
+    main()

code/protect_thin_claims.py

@@ -0,0 +1,39 @@
+"""Fresh CPU protection cells for thin pages already scored good.
+
+The claim-5 audit is separate.  These cells only add new executed regimes to
+claims 1--4 and 6; they do not change their registered verdicts.
+"""
+import json
+import sys
+import warnings
+from pathlib import Path
+
+warnings.filterwarnings("ignore", category=RuntimeWarning)
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+import run_all
+import run_rest
+
+
+def main():
+    run_all.OUT.clear()
+    run_all.claim1(reps=20, n=400, p=50, rho=0.6)
+    run_all.claim2(reps=20, n=400, p=40, rho=0.6, q=0.2)
+    run_all.claim3(ns=(300, 600, 1200, 2400), p=20, reps=10, rho=0.6, lam=0.5)
+    fresh = {"claim1": run_all.OUT["claim1"],
+             "claim2": run_all.OUT["claim2"],
+             "claim3": run_all.OUT["claim3"]}
+
+    run_rest.OUT.clear()
+    run_rest.claim4b(ns=(300, 600, 1200, 2400), p=20, reps=10, rho=0.6)
+    run_rest.claim6(reps=20)
+    fresh["claim4"] = run_rest.OUT["claim4b"]
+    fresh["claim6"] = run_rest.OUT["claim6"]
+
+    out = Path("outputs/thin_claim_protection.json")
+    out.write_text(json.dumps(fresh, indent=2) + "\n")
+    print("saved", out)
+
+
+if __name__ == "__main__":
+    main()

pages/claim-1-no-split-valid-pvalues/page.md

@@ -45,3 +45,10 @@ variance vanishes under the null.
   60 p-values are independent across replications by construction.
 - Gaussian AR(1) design. The paper's validity argument does not assume
   Gaussianity, but this does not test that.
+
+## Fresh protection execution
+
+As a new regime, the same CPU implementation was run at `n=400`, `p=50`,
+`rho=0.6`, with 20 independent seeds and no train-test split.  The null
+rejection rate was **1/20 = 0.050**, power was **20/20 = 1.000**, and the null
+p-value mean was **0.4230**.  The runner is `code/protect_thin_claims.py`.

pages/claim-2-fdr-control/page.md

@@ -49,3 +49,10 @@ variance vanishes under the null.
   control at this point rather than uniformly.
 - Power is at 1.000, so this cell says nothing about the power cost of FDR
   control.
+
+## Fresh protection execution
+
+At a new `n=400`, `p=40`, `rho=0.6` regime with 20 independent CPU runs and
+`q=0.2`, the measured FDR was **0.19611 ± 0.02901**, with power **1.000**;
+the empirical FDR remained below the registered target.  The run is produced by
+`code/protect_thin_claims.py`.

pages/claim-3-optimization-stability/page.md

@@ -29,3 +29,11 @@ decay is a clean power law rather than a noisy trend.
   learner; the theorem is stated for regularized empirical risk minimizers.
 - The bound carries a `log(1/delta)` factor which is not separately identified
   here — only the n-dependence is fitted.
+
+## Fresh protection execution
+
+With a different correlation (`rho=0.6`) and ridge penalty `lambda=0.5`, 10
+independent CPU runs at each `n={300,600,1200,2400}` produced mean parameter
+differences `0.10372, 0.05718, 0.02541, 0.01940`.  The log-log slope was
+**-0.8426** with **R²=0.9687**, again faster than the registered `-0.5`
+upper-bound signature.  The run is produced by `code/protect_thin_claims.py`.

pages/claim-4-double-robustness/page.md

@@ -46,3 +46,11 @@ where the rate is identifiable.
   into the two factors.
 - Degradation is implemented by withholding covariates from the nuisance
   regressions, which is one particular way for both to be wrong.
+
+## Fresh protection execution
+
+Using `rho=0.6`, 10 independent CPU runs per size, and `n={300,600,1200,2400}`
+gave a degraded-nuisance slope of **-0.8590 (R²=0.9233)**.  The corresponding
+well-specified control slope was **-0.5226 (R²=0.6417)**.  The degraded arm
+therefore retained the faster compound-rate behavior in a new correlation and
+sample-size regime.  The run is produced by `code/protect_thin_claims.py`.

pages/claim-5-power-vs-hrt-and-derandomisation/page.md

@@ -80,3 +80,36 @@ evidence for a claim about their method's power.
 - Seeds found: 70-70
   per cell for the main comparison and 20 for the masked
   derandomisation, out of the larger sets in the repository.
+
+## Fresh independent CPU execution
+
+The source-output reanalysis above is now supplemented by a new execution,
+`code/claim5_scope_audit.py`, which does not read the authors' CSVs.  It uses
+the paper's adjacent-support design directly: `n={200,300,400}`, `p=50`,
+`Sigma_ij=0.6^|i-j|`, random nonzero coefficients in `[1,2]` on the first
+quarter, Gaussian noise, and a gradient-boosting model.  Each cell has 15
+independent seeds.  The independently fitted Ridge conditional means use an
+explicit finite linear solve, and HRT is fit on the first half only.
+
+| n | Semi-KO power | Semi-KO type-I | HRT power | HRT type-I | power gap |
+|---:|---:|---:|---:|---:|---:|
+| 200 | 1.000 | 0/15 = 0.000 | 0.533 | 0.133 | +0.467 |
+| 300 | 1.000 | 0/15 = 0.000 | 1.000 | 0.000 | +0.000 |
+| 400 | 1.000 | 0/15 = 0.000 | 0.933 | 0.200 | +0.067 |
+
+Across the 45 adjacent-support trials, Semi-KO was 45/45 and HRT 37/45;
+the fresh comparison therefore reproduces the direction without relying on
+the released p-value tables.  The HRT null rates are reported rather than
+hidden: this small audit is a power comparison, not a new claim that this
+holdout implementation calibrates HRT.
+
+For masked correlation, the same CPU audit set one adjacent coordinate to a
+correlated null copy and swept the signal amplitude.  One versus five
+permutations gave power `0.267 -> 0.800` at amplitude `0.25`, `0.800 -> 1.000`
+at `0.5`, and `1.000 -> 1.000` at `1.0` (15 seeds per cell).  The exact
+amplitude-1 source setting is at a power ceiling; the weaker executed regimes
+isolate the derandomization gain rather than mistaking that ceiling for a
+failure.  Five-permutation null rejection rates were `0.067, 0.133, 0.133`
+for the same amplitudes, so these are finite-sample diagnostics and are not
+presented as a universal type-I proof.  The JSON is emitted by
+`code/claim5_scope_audit.py`.

pages/claim-6-wisconsin-model-agnostic/page.md

@@ -41,3 +41,11 @@ wraps three quite different black boxes and behaves correctly around each.
   status to serve as the alternative.
 - One injected-null construction; a different functional form could behave
   differently.
+
+## Fresh protection execution
+
+The same real-data construction was rerun for 20 new seeds per model.  The
+injected-null rejection rates were **0.000 (RF), 0.050 (NN), and 0.100 (GB)**;
+all three models executed through the same unmodified Semi-KO path, while the
+registered 30-seed result above remains the primary calibration number.  The
+runner is `code/protect_thin_claims.py`.