scripts/bootstrap.py

#!/usr/bin/env python
"""Paired-difference bootstrap CI (Claim 3 statistical test).

Reproduces the paper's bootstrap-CI machinery (Table 6): given two per-instance
relative-to-best gap vectors over the 70 subset instances, compute the paired
difference Delta = mean(A - B) and a 95% CI by resampling the 70 paired
differences with replacement (percentile method), plus a two-sided bootstrap
p-value.

CPU-reproducible target (strong-vs-weak, mirrors the paper's Qwen-vs-SPT+SPT
Delta ~ -50pp, p<0.001):
    LIFO:LIT  -  SPT:SPT   ->   Delta ~ -50pp, CI excludes 0, p < 0.001

Reported-not-run target (needs Qwen3-8B rollouts):
    Qwen3-8B  -  LIFO:LIT  ->   Delta = -0.099pp, 95% CI [-0.531, 0.389]

If a Qwen gap vector is available (outputs/qwen_gaps.json), the Qwen-vs-LIFO+LIT
paired test is run directly; otherwise it is labeled reported-not-run.
"""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np

ROOT = Path("/home/ubuntu/samuel/dynasched-repro")
GAPS = ROOT / "outputs" / "pdr_gaps.json"
QWEN = ROOT / "outputs" / "qwen_gaps.json"
OUT = ROOT / "outputs" / "bootstrap_ci.json"

N_BOOT = 100_000
SEED = 12345


def paired_bootstrap(a, b, n_boot=N_BOOT, seed=SEED):
    """Paired bootstrap of Delta = mean(a - b). Percentile 95% CI + 2-sided p."""
    a = np.asarray(a, float)
    b = np.asarray(b, float)
    assert a.shape == b.shape
    d = a - b
    n = len(d)
    delta = float(d.mean())
    rng = np.random.default_rng(seed)
    idx = rng.integers(0, n, size=(n_boot, n))
    boot = d[idx].mean(axis=1)
    lo, hi = np.percentile(boot, [2.5, 97.5])
    # two-sided bootstrap p-value: reflect boot dist about 0, fraction as
    # extreme as observed. Use the fraction of resamples on the opposite side
    # of 0 from delta, doubled (clamped to >= 1/n_boot).
    if delta < 0:
        p = 2.0 * float((boot >= 0).mean())
    else:
        p = 2.0 * float((boot <= 0).mean())
    p = max(p, 1.0 / n_boot)
    p = min(p, 1.0)
    return {
        "delta_pp": delta,
        "ci95_lo": float(lo),
        "ci95_hi": float(hi),
        "boot_p_two_sided": p,
        "n": n,
        "n_boot": n_boot,
    }


def main():
    g = json.loads(GAPS.read_text())
    lifo = g["lifo_lit_gap_vector"]
    spt = g["spt_spt_gap_vector"]
    gbi = g["gap_by_instance"]
    instances = g["instances"]
    ranked = g["ranked_rules"]

    results = {}

    # --- CPU-reproducible: strong (LIFO:LIT) vs weak (SPT:SPT) ---
    r = paired_bootstrap(lifo, spt)
    results["LIFO+LIT_minus_SPT+SPT"] = r
    print("=== Paired bootstrap: LIFO+LIT - SPT+SPT (CPU-reproducible) ===")
    print(
        f"  Delta = {r['delta_pp']:.3f} pp   95% CI [{r['ci95_lo']:.3f}, {r['ci95_hi']:.3f}]"
        f"   p = {r['boot_p_two_sided']:.2e}"
    )
    print("  (paper reports strong-vs-SPT+SPT Delta ~ -50pp, p<0.001)")

    # --- Methodology demo mirroring the Qwen-vs-LIFO+LIT near-zero test:
    #     LIFO+LIT (best rule) vs the 2nd-best rule, a small-Delta paired
    #     comparison expected to be statistically indistinguishable (CI spans 0).
    second = ranked[1] if ranked[0] == "LIFO:LIT" else ranked[0]
    second_vec = [gbi[k][second] for k in instances]
    rr = paired_bootstrap(lifo, second_vec)
    results["LIFO+LIT_minus_2nd_best"] = dict(rr)
    results["LIFO+LIT_minus_2nd_best"]["second_best_rule"] = second
    print(
        f"\n=== Paired bootstrap: LIFO+LIT - {second.replace(':', '+')} "
        "(2nd-best; methodology demo, CPU-reproducible) ==="
    )
    print(
        f"  Delta = {rr['delta_pp']:.3f} pp   95% CI [{rr['ci95_lo']:.3f}, "
        f"{rr['ci95_hi']:.3f}]   p = {rr['boot_p_two_sided']:.3f}"
    )
    print(
        "  Same paired-percentile-bootstrap machinery the paper applies to "
        "Qwen3-8B vs LIFO+LIT; a CI spanning 0 => indistinguishable."
    )

    # --- Qwen3-8B vs LIFO+LIT ---
    if QWEN.exists():
        q = json.loads(QWEN.read_text())
        qvec = q["qwen_gap_vector"]
        rq = paired_bootstrap(qvec, lifo)
        results["Qwen3-8B_minus_LIFO+LIT"] = rq
        results["Qwen3-8B_minus_LIFO+LIT"]["source"] = "RUN (local Qwen3-8B slice)"
        print("\n=== Paired bootstrap: Qwen3-8B - LIFO+LIT (RUN) ===")
        print(
            f"  Delta = {rq['delta_pp']:.3f} pp   95% CI [{rq['ci95_lo']:.3f}, {rq['ci95_hi']:.3f}]"
            f"   p = {rq['boot_p_two_sided']:.3f}"
        )
        print("  (paper target: Delta = -0.099pp, 95% CI [-0.531, 0.389])")
    else:
        results["Qwen3-8B_minus_LIFO+LIT"] = {
            "source": "REPORTED-NOT-RUN",
            "paper_delta_pp": -0.099,
            "paper_ci95": [-0.531, 0.389],
            "paper_fdr_p": 0.397,
            "note": (
                "Requires Qwen3-8B L1-Direct rollouts on the 70 subset "
                "instances. No gap vector available; the bootstrap "
                "machinery is reproduced above on the CPU pair."
            ),
        }
        print("\n=== Qwen3-8B - LIFO+LIT: REPORTED-NOT-RUN ===")
        print("  paper: Delta = -0.099pp, 95% CI [-0.531, 0.389], FDR p = 0.397")

    OUT.write_text(json.dumps(results, indent=2))
    print(f"\nWrote {OUT}")


if __name__ == "__main__":
    main()