# Verification checks for P4 (reads out/results.json, out/table1.csv, data/design.csv)
import json, numpy as np, pandas as pd
r = json.load(open('out/results.json')); t = pd.read_csv('out/table1.csv'); d = pd.read_csv('data/design.csv')
ok = lambda c, m: print(("PASS " if c else "FAIL ") + m)
ok(len(t) == 27 and set(t.estimator) == {'cr0','cr1','wcb'}, "table1 has 27 rows, 3 estimators")
ok(len(r['coverage']) == 27 and len(r['mc_se_coverage']) == 27 and len(r['rel_bias_var']) == 18, "key counts 27/27/18")
ok(t[t.estimator=='wcb'].rel_bias_var.isna().all(), "rel_bias_var empty for wcb")
ok(all(0 <= v <= 1 for v in r['coverage'].values()), "coverages in [0,1]")
ok(all(v > 0 for v in r['mc_se_coverage'].values()), "MC SEs positive")
ok(all(abs(r['mc_se_coverage'][k] - np.sqrt(p*(1-p)/r['R'])) < 1e-12 for k, p in r['coverage'].items()), "MC SE = sqrt(p(1-p)/R)")
ok(all(abs(v*r['R'] - round(v*r['R'])) < 1e-9 for v in r['coverage'].values()), "coverage*R integer")
f = sorted(k for k, v in r['coverage'].items() if v < 0.9)
ok(f == sorted(r['failure_cells']) and r['n_failure_cells'] == len(f), "failure_cells consistent")
nm = sorted(k for k, v in r['coverage'].items() if abs(v-0.95) <= 2*r['mc_se_coverage'][k])
ok(nm == sorted(r['nominal_cells']), "nominal_cells consistent")
ok(r['n_analysed'] == r['n_cells']*r['R'] == 18000 and r['R']==2000 and r['B']==399, "n_analysed=9*2000, R,B unchanged")
ok(r['random_seed'] == int(d.seed[d.cell_id=='C1'].iloc[0]), "random_seed = C1 seed")
ex = pd.read_csv('out/exclusions.csv'); ok(18000 - ex.n_removed.sum() == r['n_analysed'], "exclusion ledger sums")
# CR1/CR0 ratio must be G/(G-1)
ok(all(abs(r['rel_bias_var']['cr1'+k[3:]]/r['rel_bias_var'][k] - (int(k.split('_')[1][1:])/(int(k.split('_')[1][1:])-1))) < 1e-9
       for k in r['rel_bias_var'] if k.startswith('cr0')), "rel_bias cr1/cr0 = G/(G-1)")
# analytic E[V_CR0]/Var(ybar) vs MC estimate
print("Analytic vs MC relative bias of CR0 (z = (MC-exact)/MCSE):")
for _, row in d.iterrows():
    n = np.array([int(x) for x in row.cluster_sizes.split('|')], float); N = n.sum(); icc = row.icc
    V = (icc*np.sum(n**2) + (1-icc)*N)/N**2
    VT = icc*n**2 + (1-icc)*n
    EV0 = np.sum(VT*(1-2*n/N) + n**2*V)/N**2
    k = f"cr0_G{row.G}_icc{icc:.1f}"
    z = (r['rel_bias_var'][k] - EV0/V)/r['mc_se_rel_bias_var'][k]
    ok(abs(z) < 3, f"  {k}: exact {EV0/V:.4f}  MC {r['rel_bias_var'][k]:.4f}  z={z:.2f}")
    dg = r['cell_diagnostics'][row.cell_id]
    ok(abs(dg['var_exact'] - V) < 1e-15, f"  {row.cell_id} var_exact matches formula")
# monotonicity sanity: cr1 halfwidth > cr0 halfwidth
w = t.pivot(index='cell_id', columns='estimator', values='mean_ci_halfwidth')
ok((w.cr1 > w.cr0).all(), "CR1 intervals wider than CR0 in every cell")