"""Stage 2: EB-variant sanity + H3 placebo + H4 CIs + H5-H8 -> side_results.json (merge)."""
import json, math, random, statistics, sys
sys.path.insert(0, ".")
from collections import defaultdict

exec(open("run3/side_locked.py").read().split("# ---------- A1a")[0])  # frame builders

rated = [f for f in frame if f["rated"]]
def pearson(xs, ys):
    n = len(xs); mx = sum(xs)/n; my = sum(ys)/n
    cov = sum((x-mx)*(y-my) for x, y in zip(xs, ys))
    vx = sum((x-mx)**2 for x in xs); vy = sum((y-my)**2 for y in ys)
    return cov/(vx*vy)**.5
def zscore(xs):
    m = sum(xs)/len(xs); v = sum((x-m)**2 for x in xs)/len(xs)
    s = v**.5 if v > 0 else 1.0
    return [(x-m)/s for x in xs]

# reuse ols_fe from stage1 file
_src = open("run3/side_locked.py").read()
_ols = _src.split("def ols_fe(rows_):")[1].split("# bootstrap CI over reviewers")[0]
exec("def ols_fe(rows_):" + _ols)
out = json.load(open("run3/side_results.json"))

# ---- EB sanity ladder: w in {0 (=base), 0.5, 1.0 (=no FE)} ----
groups0 = defaultdict(list)
for i, f in enumerate(rated):
    groups0[f["reviewer"]].append(i)
ys0 = [f["q"] for f in rated]
gm = sum(ys0)/len(ys0)
eb_ladder = {}
for w in (0.0, 0.5, 1.0):
    rows2 = [dict(f) for f in rated]
    for g, idxs in groups0.items():
        mg = sum(ys0[i] for i in idxs)/len(idxs)
        for i in idxs:
            rows2[i]["q"] = gm + (1-w)*(ys0[i]-mg)*0 + (ys0[i]-mg)*0 + ((ys0[i]-mg)) # placeholder
    # correct: transformed y = gm + w_norm*(...) -- implement directly:
    for g, idxs in groups0.items():
        mg = sum(ys0[i] for i in idxs)/len(idxs)
        for i in idxs:
            # w=0 -> pure demeaning removed later anyway; define transformed = gm*(1-c)+ (mg + c*(y-mg))
            c = 1.0 - w
            rows2[i]["q"] = (1-c)*mg + c*ys0[i]
    bb, _, _ = ols_fe(rows2)
    eb_ladder[str(w)] = round(bb, 4)
print("EB ladder (w=0 full-demean ... w=1 no-demean):", eb_ladder)
out["H2_eb_ladder_sanity"] = eb_ladder

# ---- H3 placebo: permute q within reviewer blocks ----
rng = random.Random(20260823)
r_obs = out["A1a_r_absdev_quality"]
B = 5000
cnt = 0
by_rev = defaultdict(list)
for f in rated:
    by_rev[f["reviewer"]].append(f)
qs_by_rev = {g: [f["q"] for f in lst] for g, lst in by_rev.items()}
xs_all = [f["absdev"] for f in rated]
for _ in range(B):
    perm_q = []
    for g, lst in by_rev.items():
        qs = qs_by_rev[g][:]
        rng.shuffle(qs)
        perm_q.extend(qs)
    # keep alignment: rebuild in same order as rated list grouped iteration
    rr = pearson(xs_all, perm_q)
    if abs(rr) >= abs(r_obs) - 1e-12:
        cnt += 1
p_h3 = (cnt+1)/(B+1)
out["H3_placebo_p_within_reviewer"] = round(p_h3, 4)
print("H3 within-reviewer permutation p =", round(p_h3, 4))

# ---- H4 CIs ----
rng2 = random.Random(20260823)
B2 = 10000
rs_iid, rs_clu = [], []
for _ in range(B2):
    samp = [rated[rng2.randrange(len(rated))] for _ in rated]
    rs_iid.append(pearson([f["absdev"] for f in samp], [f["q"] for f in samp]))
revs_list = sorted(by_rev)
for _ in range(B2):
    samp = []
    for _ in revs_list:
        samp.extend(by_rev[rng2.choice(revs_list)])
    rs_clu.append(pearson([f["absdev"] for f in samp], [f["q"] for f in samp]))
rs_iid.sort(); rs_clu.sort()
out["H4_r_iid_ci95"] = [round(rs_iid[int(.025*B2)], 4), round(rs_iid[int(.975*B2)-1], 4)]
out["H4_r_cluster_ci95"] = [round(rs_clu[int(.025*B2)], 4), round(rs_clu[int(.975*B2)-1], 4)]
print("H4 r CI iid:", out["H4_r_iid_ci95"], "cluster:", out["H4_r_cluster_ci95"])

# ---- H5 nonlinearity ----
srt = sorted(rated, key=lambda f: f["absdev"])
k4 = len(srt)//4
quart = {}
for qi in range(4):
    chunk = srt[qi*k4:(qi+1)*k4] if qi < 3 else srt[3*k4:]
    mq = statistics.mean([f["q"] for f in chunk])
    mad = statistics.mean([f["absdev"] for f in chunk])
    boot = []
    r3 = random.Random(77+qi)
    for _ in range(2000):
        smp = [chunk[r3.randrange(len(chunk))] for _ in chunk]
        boot.append(statistics.mean([f["q"] for f in smp]))
    boot.sort()
    quart[f"Q{qi+1}"] = {"mean_absdev": round(mad, 3), "mean_q": round(mq, 3),
                          "ci95": [round(boot[50], 2), round(boot[1949], 2)]}
out["H5_quartiles"] = quart
sp_x = [f["absdev"] for f in rated]
def rankify(v):
    order = sorted(range(len(v)), key=lambda i: v[i])
    rk = [0.0]*len(v)
    i = 0
    while i < len(order):
        j = i
        while j+1 < len(order) and v[order[j+1]] == v[order[i]]:
            j += 1
        avg = (i+j)/2 + 1
        for t2 in range(i, j+1):
            rk[order[t2]] = avg
        i = j+1
    return rk
out["H5_spearman"] = round(pearson(rankify(sp_x), rankify([f["q"] for f in rated])), 4)

# ---- H6 secondary dims ----
sec = {}
for d in ("nov", "cla", "sig"):
    sec[d] = round(pearson([f[f"absdev_{d}"] for f in rated], [f["q"] for f in rated]), 4)
out["H6_secondary_dims_r"] = sec

# ---- H7 sparsity ----
allr = list(rev.values())
n0 = sum(1 for r in allr if r.get("ratings_count", 0) == 0)
out["H7_share_unrated"] = round(n0/len(allr), 4)
unr_q = [r["quality_score"] for r in allr if r.get("ratings_count", 0) == 0 and r.get("quality_score") is not None]
rat_q = [r["quality_score"] for r in allr if r.get("ratings_count", 0) >= 1 and r.get("quality_score") is not None]
vr_un = statistics.pvariance(unr_q); vr_ra = statistics.pvariance(rat_q)
out["H7_var_ratio_un_over_rated"] = round(vr_un/vr_ra, 3)
out["H7_unrated_mean_sd"] = [round(statistics.mean(unr_q), 2), round(statistics.pstdev(unr_q), 2)]
# Mann-Whitney rated vs unrated on text_len / position / n_other / reviewer rank
def mannwhitney_u(a, b):
    allv = [(v, 0) for v in a] + [(v, 1) for v in b]
    allv.sort(key=lambda t: t[0])
    ranks = [0.0]*len(allv)
    i = 0
    while i < len(allv):
        j = i
        while j+1 < len(allv) and allv[j+1][0] == allv[i][0]:
            j += 1
        avg = (i+j)/2 + 1
        for t2 in range(i, j+1):
            ranks[t2] = avg
        i = j+1
    Ra = sum(ranks[t2] for t2 in range(len(allv)) if allv[t2][1] == 0)
    na, nb = len(a), len(b)
    U = Ra - na*(na+1)/2
    mu = na*nb/2
    sd = (na*nb*(na+nb+1)/12)**.5
    z = (U-mu)/sd if sd > 0 else 0
    from math import erf
    p2 = 2*(1-0.5*(1+erf(abs(z)/math.sqrt(2))))
    return U, z, p2

rated_set = set(id(f) for f in rated)
unrated_f = [f for f in frame if not f["rated"]]
mw = {}
U, z, p2 = mannwhitney_u([f["text_len"] for f in rated], [f["text_len"] for f in unrated_f])
mw["text_len"] = {"z": round(z, 2), "p": round(p2, 4)}
U, z, p2 = mannwhitney_u([f["position"] for f in rated], [f["position"] for f in unrated_f])
mw["position"] = {"z": round(z, 2), "p": round(p2, 4)}
U, z, p2 = mannwhitney_u([f["n_other"] for f in rated], [f["n_other"] for f in unrated_f])
mw["n_other"] = {"z": round(z, 2), "p": round(p2, 4)}
out["H7_mw_rated_vs_unrated"] = mw
first_share_rated = statistics.mean([1.0 if f["position"] == 0 else 0.0 for f in rated])
first_share_unrated = statistics.mean([1.0 if f["position"] == 0 else 0.0 for f in unrated_f])
out["H7_first_review_share"] = {"rated": round(first_share_rated, 3), "unrated": round(first_share_unrated, 3)}

# ---- H8 fidelity ----
agg = {x["rid"]: x for x in json.load(open("run3/agg_rows.json"))}
rr = [f for f in rated if f["rid"] in agg]
cm = [agg[f["rid"]]["comp_mean"] for f in rr]
qq = [f["q"] for f in rr]
r_cq = pearson(cm, qq)
b0 = statistics.mean(qq) - r_cq*(statistics.pstdev(qq)/statistics.pstdev(cm))*statistics.mean(cm)
sl = r_cq*statistics.pstdev(qq)/statistics.pstdev(cm)
resid = [qq[i] - (b0 + sl*cm[i]) for i in range(len(rr))]
nr = [agg[f["rid"]]["n_raters"] for f in rr]
ad = [f["absdev"] for f in rr]
t_res_mean = statistics.mean(resid)/(statistics.pstdev(resid)/math.sqrt(len(resid)))
slope_nr = pearson(nr, resid)*statistics.pstdev(resid)/statistics.pstdev(nr) if statistics.pstdev(nr) > 0 else float('nan')
out["H8"] = {"n": len(rr), "r_comp_quality": round(r_cq, 4),
             "resid_mean_t": round(t_res_mean, 2),
             "slope_resid_on_n_raters": round(slope_nr, 4),
             "r_resid_absdev": round(pearson(resid, ad), 4),
             "share_abs_gt2": round(sum(1 for e in resid if abs(e) > 2)/len(resid), 3)}
print("H8:", json.dumps(out["H8"]))

json.dump(out, open("run3/side_results.json", "w"), indent=1)
print("\n=== SIDE RESULTS COMPLETE ===")
print(json.dumps({k: out[k] for k in ["A1a_r_absdev_quality", "A1b_b1", "A1b_boot_ci95",
                                       "H3_placebo_p_within_reviewer", "H6_secondary_dims_r",
                                       "H7_share_unrated"]}, indent=1))