#!/usr/bin/env python3
"""Instrument check for blur.py. NOT a test of the hypothesis, and not in PREREG-blur.md.

Added after the primary run returned INCONCLUSIVE, because "no effect detected" and "an
instrument that could not detect one" look identical from the outside, and only one of them is
a finding. Two questions:

  1. Calibration -- with 14 clusters, does the percentile cluster bootstrap actually cover 95%?
     A too-narrow interval would make the pre-registered decision table wrong in my favour.
  2. Power -- at the pre-registered effect size (|rho| = 0.20), how often would this design's
     interval exclude zero? If that is low, "inconclusive" says nothing about the world.

Writes blur_power.json. Separate file so the pre-registered analysis stays untouched.

    python3 blur_power.py
"""
import json, os, sys
import numpy as np

# FOOLDET_DIR lets a long run be launched from a copy elsewhere without the copy's location
# silently becoming the data directory -- which is how the first launch of this file died.
sys.path.insert(0, os.environ.get("FOOLDET_DIR", os.path.dirname(os.path.abspath(__file__))))
import blur                                   # reuse the exact estimator under test

HERE = blur.HERE                              # wherever the module really loaded from, not __file__
NSIM, NBOOT, SEED = 150, 800, 20260830   # ~10 min on this one-core box [[box-is-one-core]]
TARGETS = [0.0, 0.10, 0.20, 0.30]


def ci(rows, sources, by_src, rng, nboot):
    b = []
    for _ in range(nboot):
        draw = rng.integers(0, len(sources), len(sources))
        rs = [{**r, "source": f"{r['source']}#{k}"}
              for k, di in enumerate(draw) for r in by_src[sources[di]]]
        b.append(blur.estimate(rs)[0])
    return (float(v) for v in np.percentile(b, [2.5, 97.5]))


def main():
    base = json.load(open(f"{HERE}/blur.json"))
    rows, provenance = blur.load_corpus()   # same resolver as the analysis, so a clone gets the
    print(f"  input: {provenance}")         # shipped measures rather than a machine-local path
    sources = sorted({r["source"] for r in rows})
    by_src = {s: [r for r in rows if r["source"] == s] for s in sources}

    # within-source standardised rank of the real blur measure -- the signal we inject against
    x = blur.within_source(rows, "logvol")
    x = x / x.std()

    rng = np.random.default_rng(SEED)
    out = {"n_sim": NSIM, "n_boot": NBOOT, "seed": SEED, "targets": {}}
    for t in TARGETS:
        excl, rhos = 0, []
        for _ in range(NSIM):
            noise = rng.standard_normal(len(rows))
            y = -t * x + np.sqrt(max(1 - t * t, 0.0)) * noise   # negative = H1's direction
            rs = [{**r, "score": float(v)} for r, v in zip(rows, y)]
            bysrc_s = {s: [r for r in rs if r["source"] == s] for s in sources}
            lo, hi = ci(rs, sources, bysrc_s, rng, NBOOT)
            rhos.append(blur.estimate(rs)[0])
            if not (lo <= 0 <= hi):
                excl += 1
        out["targets"][f"{t:.2f}"] = {
            "injected_rho": -t, "mean_recovered_rho": round(float(np.mean(rhos)), 4),
            "frac_ci_excludes_zero": excl / NSIM,
        }
        lbl = "false-positive rate" if t == 0 else "power"
        print(f"  injected rho {-t:+.2f}  recovered {np.mean(rhos):+.3f}  "
              f"{lbl} {excl / NSIM:.2f}")

    fp = out["targets"]["0.00"]["frac_ci_excludes_zero"]
    pw = out["targets"]["0.20"]["frac_ci_excludes_zero"]
    out["observed_rho"] = base["study_a"]["rho_within_source"]
    out["calibration_ok"] = fp <= 0.10
    out["powered_at_prereg_effect"] = pw >= 0.80
    out["reading"] = (
        "The interval is calibrated and the design would detect the pre-registered effect, so "
        "INCONCLUSIVE is evidence of a small-or-absent effect."
        if out["calibration_ok"] and out["powered_at_prereg_effect"] else
        "The design cannot reliably detect the pre-registered effect, so INCONCLUSIVE means the "
        "study is uninformative about the direction -- it is not evidence against it.")
    json.dump(out, open(f"{HERE}/blur_power.json", "w"), indent=1)
    print(f"\n  calibration_ok={out['calibration_ok']} (false-positive {fp:.2f}, want <=0.10)")
    print(f"  powered_at_0.20={out['powered_at_prereg_effect']} (power {pw:.2f}, want >=0.80)")
    print(f"  {out['reading']}")


if __name__ == "__main__":
    main()