"""N1, P1 and D1, printed before the real statistic exists. Step 4 of five. python research/red-sky/controls.py [--distance 300] `rule.py` orders these before `analyse.py` and this file obeys that by never computing the real delta at all: it prints the null, the floor and the downwind dummy, and stops. The real number is `analyse.py`'s and nothing here reads it. N1 the label permuted within (site, calendar month) blocks, 200 draws. The stratified delta must come back within +-1.0 point of zero and "the CI excludes zero" must fire in at most 10 % of draws. Otherwise the pipeline is measuring itself and no verdict may be published. P1 an effect of exactly X points injected by flipping a random subset of exposed days' outcomes to fine, X in {0, 2, 5, 10}, 200 injections each. The floor is the smallest X recovered with a CI excluding zero in >= 90 % of them. X = 0 is the failing baseline and MUST NOT fire above about 5 %; a positive control without one is a formality. THE INJECTION IS BUILT ON A PERMUTED LABEL, NOT THE REAL ONE. `rule.py` says "starting from the real labels", and taken literally that makes the X = 0 case the real statistic with no randomness in it — 200 identical draws, and the real answer printed in the controls file, which is the one thing this file exists to avoid. So each draw permutes the label within blocks first (destroying whatever the real association is, and keeping every site and seasonal rate) and injects X on top of that. The floor this measures is the detector's sensitivity at this sample size, which is what the floor is for, and it cannot see the answer. D1 the same statistic with the DOWNWIND point substituted for the upwind one. Not a null — a clear eastern horizon at sunset correlates with a clear western one — so it is printed as a comparison and never as a test. The bootstrap inside a control uses 1,000 resamples rather than 10,000: these are rates over 200 draws, not a published interval, and the saving is what makes 200 draws affordable. `analyse.py` uses the full 10,000. """ import os import sys import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import stats as S # noqa: E402 (the pack is run as scripts, not a package) OUT = os.path.expanduser("~/data/red-sky") CONTROL_RESAMPLES = 1000 DRAWS = 200 BAR = 5.0 def ci_excludes_zero(packed, seed): point, lo, hi = S.interval(packed, CONTROL_RESAMPLES, seed) return point, (lo > 0) or (hi < 0) def permute_within_blocks(rows, label_key, rng): """A copy of the labels, shuffled inside each (site, calendar month).""" by_block = {} for i, r in enumerate(rows): by_block.setdefault((r["site"], r["date"][:7]), []).append(i) labels = [r[label_key] for r in rows] for idx in by_block.values(): values = [labels[i] for i in idx] rng.shuffle(values) for i, v in zip(idx, values): labels[i] = v return labels def with_labels(rows, labels, label_key): return [dict(r, **{label_key: l}) for r, l in zip(rows, labels)] def inject(rows, label_key, outcome_key, points, rng): """Flip exposed days' outcomes to fine until the effect is +points.""" exposed = [i for i, r in enumerate(rows) if r[label_key] == 1 and r[outcome_key] == 0] n_exposed = sum(1 for r in rows if r[label_key] == 1 and r[outcome_key] in (0, 1)) k = int(round(n_exposed * points / 100.0)) out = [dict(r) for r in rows] if k > 0 and exposed: for i in rng.choice(exposed, size=min(k, len(exposed)), replace=False): out[i][outcome_key] = 1 return out def main(): km = 300 if "--distance" in sys.argv: km = int(sys.argv[sys.argv.index("--distance") + 1]) rows = S.read_rows(os.path.join(OUT, f"site-days-{km}.csv")) print(f"{len(rows)} usable site-days at {km} km\n") rng = np.random.default_rng(S.SEED) print("N1 the null: label permuted within (site, calendar month)") points, fired = [], 0 for d in range(DRAWS): labels = permute_within_blocks(rows, "red_night", rng) packed = S.blocks(with_labels(rows, labels, "red_night"), "red_night", "fine_next") p, excl = ci_excludes_zero(packed, S.SEED + d) points.append(p) fired += excl mean = float(np.mean(points)) rate = 100.0 * fired / DRAWS print(f" mean delta over {DRAWS} draws {mean:+.3f} points " f"(rule: within +-1.0)") print(f" CI excluded zero in {rate:.1f} % (rule: at most 10 %)") ok = abs(mean) <= 1.0 and rate <= 10.0 print(f" N1 {'PASSES' if ok else 'FAILS — everything is void'}\n") if not ok: return 1 print("P1 injected effect, and the floor") floor = None for x in (0, 2, 5, 10): fired = 0 for d in range(DRAWS): labels = permute_within_blocks(rows, "red_night", rng) base = with_labels(rows, labels, "red_night") injected = inject(base, "red_night", "fine_next", x, rng) packed = S.blocks(injected, "red_night", "fine_next") _, excl = ci_excludes_zero(packed, S.SEED + 7919 * d) fired += excl rate = 100.0 * fired / DRAWS mark = "" if x == 0: mark = " <- the failing baseline; above ~5 % means the gate is hollow" elif floor is None and rate >= 90.0: floor = x mark = " <- floor" print(f" +{x:2d} points recovered in {rate:5.1f} % of {DRAWS}{mark}") print(f" floor {floor} points against a bar of {BAR:.0f}") if floor is None or floor > BAR: print(" THE SURVIVED BRANCH CANNOT FIRE: this arm is " "VOID: UNFAILABLE and must be published as such.") print() print("D1 the downwind dummy: the EAST point at sunset, in place of west") dummy = [dict(r, red_night=int(r["clear_e_at_set"] == 1 and r["lit"] == 1)) for r in rows] point, lo, hi = S.interval(S.blocks(dummy, "red_night", "fine_next"), CONTROL_RESAMPLES, S.SEED) print(f" delta {point:+.3f} points, 95 % [{lo:+.3f}, {hi:+.3f}]") print(" A comparison and not a test: a clear eastern horizon at sunset " "is correlated with a clear western one.") print("\nThe real statistic is not computed in this file. Next: analyse.py.") return 0 if __name__ == "__main__": raise SystemExit(main())