"""r11_review.py: analyses completed or added after the independent review of the Phase 4 draft.

Each block answers a review item and is logged in inputs/deviations4.json before its result was computed:
  DV15  planned E1 and E8 outputs the first draft left out: clipping above d9_threshold by machine, sound channel
        and label (labels 0 and 1 only), counted from the same path table as D9;
  DV13  post hoc, exploratory: the AUC of the Sound0 full-scale share itself (no model, no split) on both
        laboratory machines, with family-stratified run-bootstrap intervals;
  DV14  post hoc gauge of the refit variability the registered bootstrap leaves out: delete-one-run jackknife
        standard errors from the refits of check A1 against the standard deviation of the registered bootstrap
        draws, which are recomputed here and asserted to reproduce the registered intervals;
  DV16  post hoc, exploratory: H3 agreement against independence within runs (frozen kappa functions with the
        run as the stratum), with a run-bootstrap interval of the excess and no p-value;
  DV17  post hoc, descriptive: frequencies of the dominant off-harmonic lines of Acc0 and Sound2 (quartiles by
        label, Spearman correlation with the commanded spindle speed among label 1 paths), and the frequencies of
        the two Fig. 2 examples.
It also registers the counts of Table 1 by tool family and workpiece material (review item R12) and the times of
the freeze and of the confirmatory run from the authors' own files (review item R9). Values go to values4.json
with the prefix ph_ for the post hoc analyses; nothing here enters a confirmatory verdict.
"""
import sys, json, os, time, datetime
sys.path.insert(0, "/home/claude/mill/phase4/code")
sys.path.insert(0, "/home/claude/mill/phase3/code")
sys.path.insert(0, "/home/claude/mill/code")
import r_freeze
r_freeze.assert_freeze()

import numpy as np
import pandas as pd
from scipy.stats import spearmanr
from vals import Values
import q_common as QC
import q_confirm as QF
import r_fast as RF

P3, P4 = "/home/claude/mill/phase3", "/home/claude/mill/phase4"
V2 = json.load(open("/home/claude/mill/outputs/values.json"))
V3 = json.load(open(f"{P3}/outputs/values3.json"))
V4 = json.load(open(f"{P4}/outputs/values4.json"))
V = Values(path=f"{P4}/outputs/values4.json", script="r11_review.py")
g = lambda n: (V4[n] if n in V4 else V3[n] if n in V3 else V2[n])["value"]
POP = "labelled cutting paths (label.csv rows)"
B, SEED, LEVEL, FRAC = int(g("boot_b")), int(g("seed_analysis")), float(g("ci_level")), float(g("thr_kappa_frac"))
THR = float(g("d9_threshold"))
METHOD = "percentile interval at ci_level over boot_b run resamples drawn within tool families (q_confirm.boot_weights); post hoc, exploratory, no p-value"
est = lambda f, s: {"population": POP, "filter": f, "statistic": s}
t0 = time.time()


def reg(name, value, unit, e, kind="value"):
    return V.set(name, float(value) if isinstance(value, (float, np.floating)) else (int(value) if isinstance(value, (np.integer,)) else value), unit, e, kind=kind)


def reg_ci(name, point, draws, e):
    lo, hi = QC.percentile_ci(draws, LEVEL)
    V.interval(name, point, lo, hi, LEVEL, METHOD, e)


def log(msg):
    print(f"{msg} ({time.time() - t0:.0f} s)", flush=True)


paths = QC.load_paths()
y = paths.y.values
arr = np.load(f"{P4}/outputs/r01_arrays.npz", allow_pickle=True)
assert (arr["h1_run_id"] == paths.run_id.values.astype(str)).all()
CH_NAME = {"s0": "Sound0", "s1": "Sound1", "s2": "Sound2", "acc0": "Acc0", "acc1": "Acc1"}

# ------------------------------------------------------------------ DV15: E1 clipping above the D9 threshold by machine, channel and label
for mach in ("imi_vm20i", "imi_vmx30ui"):
    d = QC.load_paths(mach)
    for ch in ("s0", "s1", "s2"):
        col = f"{ch}_fullscale"
        if col not in d.columns or d[col].notna().sum() == 0:
            continue
        ok = d[col].notna().values
        for lab in (0, 1):
            m = ok & (d.y.values == lab)
            share = float((d.loc[m, col] > THR).mean())
            f = f"{mach}, label {lab}, {CH_NAME[ch]}"
            if mach == "imi_vm20i" and ch == "s0":      # seen in Phase 2 (D9); asserted equal and not registered twice
                assert abs(share - V2[f"d9_paths_gt1pct_imi_vm20i_s0_l{lab}"]["value"]) < 1e-12
            else:
                reg(f"e1_{mach}_{ch}_l{lab}_gt1pct", share, "fraction", est(f, "share of paths with more than d9_threshold of their samples at full scale (E1)"))
            reg(f"e1_{mach}_{ch}_l{lab}_n", int(m.sum()), "paths", est(f, "paths with the channel recorded (E1)"))
    ok = d.s0_fullscale.notna().values
    reg(f"e1_{mach}_s0_l01_gt1pct", float((d.s0_fullscale[ok] > THR).mean()), "fraction",
        est(f"{mach}, labels 0 and 1, Sound0", "share of paths with more than d9_threshold of their samples at full scale (E1; label 2 left out)"))
log("E1 shares")

# ------------------------------------------------------------------ DV13: AUC of the Sound0 full-scale share itself
for mach in ("imi_vm20i", "imi_vmx30ui"):
    d = QC.load_paths(mach).dropna(subset=["s0_fullscale"]).reset_index(drop=True)
    yy, s = d.y.values, d.s0_fullscale.values.astype(float)
    Wm = QF.boot_weights(d, B, SEED)
    RF.check_against_frozen(yy, s, Wm)
    e = est(f"{mach}, labels 0 and 1, Sound0", "AUC of the share of samples at full scale itself for label 1 against label 0 (no model, no split; ties count half)")
    reg(f"ph_clip_auc_raw_{mach}", QC.auc_w(yy, s), "AUC", e)
    reg_ci(f"ph_clip_auc_raw_{mach}_ci", f"ph_clip_auc_raw_{mach}", RF.auc_draws(yy, s, Wm), e)
    reg(f"ph_clip_auc_raw_{mach}_n", int(len(d)), "paths", est(f"{mach}, labels 0 and 1, Sound0", "paths"))
    reg(f"ph_clip_auc_raw_{mach}_runs", int(d.run_id.nunique()), "runs", est(f"{mach}, labels 0 and 1, Sound0", "runs"))
log("DV13 raw clipping AUC")

# ------------------------------------------------------------------ DV14: jackknife over the A1 refits against the registered bootstrap
W1 = QF.boot_weights(paths, B, SEED)
L = pd.read_csv(f"{P4}/outputs/a1_leave_one_run_out.csv")
assert len(L) == g("a1_n_runs") == paths.run_id.nunique()
lab1 = paths[paths.y == 1].reset_index(drop=True)
assert (arr["h3_run_id"] == lab1.run_id.values.astype(str)).all() and np.allclose(arr["h3_start"], lab1.start.values)
ia, ib = arr["h3_acc0"].astype(float), arr["h3_s2"].astype(float)
assert np.isfinite(ia).all() and np.isfinite(ib).all()
ia, ib = ia.astype(int), ib.astype(int)
W3 = QF.boot_weights(lab1, B, SEED)
draws = {}
for H in ("h1", "h2"):
    d_run = RF.auc_rep_draws(y, arr[f"{H}_scores_run"], W1)
    d_fam = RF.auc_draws(y, arr[f"{H}_scores_family"], W1)
    draws[H] = d_run - d_fam
dk = np.array([QC.kappa_w(ia, ib, w) for w in W3])
draws["h3"] = dk
for H, name in (("h1", "h1_drop"), ("h2", "h2_drop"), ("h3", "h3_kappa")):
    lo, hi = QC.percentile_ci(draws[H], LEVEL)
    assert abs(lo - V4[f"{name}_ci"]["lo"]) < 1e-12 and abs(hi - V4[f"{name}_ci"]["hi"]) < 1e-12, (H, lo, hi)
log("registered bootstrap intervals reproduced")
top_runs = set()
for H, name, col in (("h1", "h1_drop", "h1"), ("h2", "h2_drop", "h2"), ("h3", "h3_kappa", "h3")):
    th = L[col].values.astype(float); n = len(th)
    dev2 = (th - th.mean()) ** 2
    se = float(np.sqrt((n - 1) / n * dev2.sum()))
    sd = float(np.std(draws[H], ddof=1))
    f = "imi_vm20i, labels 0 and 1" if H != "h3" else "imi_vm20i, label 1, Acc0 and Sound2 spectra"
    reg(f"ph_jack_se_{name}", se, "AUC" if H != "h3" else "kappa",
        est(f + ", each of the a1_n_runs runs left out in turn (refits of check A1)", f"delete-one-run jackknife standard error of {name}: square root of (n minus 1)/n times the sum of squared deviations of the leave-one-run-out values from their mean"))
    reg(f"ph_boot_sd_{name}", sd, "AUC" if H != "h3" else "kappa",
        est(f, f"standard deviation of the boot_b registered bootstrap draws of {name} (predictions held fixed; recomputed from the saved scores and asserted to reproduce the registered interval)"))
    reg(f"ph_jack_ratio_{name}", se / sd, "ratio", est(f, f"ph_jack_se_{name} divided by ph_boot_sd_{name}"))
    if H != "h3":   # where the estimate lies on the jackknife scale, against zero and against the threshold
        pt, thr_ = g(name), g("thr_auc_drop")
        reg(f"ph_jack_z0_{name}", pt / se, "standard errors", est(f, f"{name} divided by ph_jack_se_{name}"))
        reg(f"ph_jack_zthr_{name}", (thr_ - pt) / se, "standard errors", est(f, f"thr_auc_drop minus {name}, divided by ph_jack_se_{name}"))
    order = np.argsort(-dev2)
    if H != "h3":
        top_runs |= {int(L.run.values[i].split(":")[1]) for i in order[:3]}
        reg(f"ph_jack_top3_share_{name}", float(dev2[order[:3]].sum() / dev2.sum()), "fraction",
            est(f + ", refits of check A1", f"share of the jackknife sum of squares of {name} that comes from the three runs with the largest deviations"))
tr = sorted(top_runs)
reg("ph_jack_top_runs", ", ".join(str(r) for r in tr[:-1]) + " and " + str(tr[-1]), "text",
    est("imi_vm20i, labels 0 and 1, refits of check A1", "runs among the three largest jackknife deviations of h1_drop or of h2_drop"), kind="text")
reg("ph_jack_top_runs_n", len(tr), "runs", est("imi_vm20i, labels 0 and 1, refits of check A1", "number of runs in ph_jack_top_runs"))
log("DV14 jackknife")

# ------------------------------------------------------------------ DV16: H3 against independence within runs
strata = lab1.run_id.values
k = QC.kappa_w(ia, ib)
assert abs(k - g("h3_kappa")) < 1e-12
kn, km = QC.kappa_null_w(ia, ib, strata), QC.kappa_max_w(ia, ib, strata)
thr = QC.kappa_threshold(kn, km, FRAC)
dn = np.array([QC.kappa_null_w(ia, ib, strata, w) for w in W3])
f16 = "imi_vm20i, label 1, Acc0 and Sound2 spectra (D5 settings), strata: runs"
reg("ph_h3_run_kappa_null", kn, "kappa", est(f16, "kappa expected if the Acc0 and Sound2 indicators were independent within each run (q_common.kappa_null_w)"))
reg("ph_h3_run_kappa_max", km, "kappa", est(f16, "largest kappa the per-run shares allow (q_common.kappa_max_w)"))
reg("ph_h3_run_thr", thr, "kappa", est(f16, "ph_h3_run_kappa_null plus thr_kappa_frac of the distance to ph_h3_run_kappa_max (the H3 rule with runs as strata)"))
e = est(f16, "h3_kappa minus ph_h3_run_kappa_null, the null recomputed in every resample")
reg("ph_h3_run_excess", k - kn, "kappa", e)
reg_ci("ph_h3_run_excess_ci", "ph_h3_run_excess", dk - dn, e)
reg("ph_h3_run_margin", k - thr, "kappa", est(f16, "h3_kappa minus ph_h3_run_thr"))
reg("ph_h3_run_n_runs", int(len(np.unique(strata))), "runs", est(f16, "runs (strata)"))
log("DV16 within-run null")

# ------------------------------------------------------------------ DV17: frequencies of the dominant off-harmonic lines
BASE_IND = dict(comb_half_bins=int(g("d5_comb_bins")), band=(float(g("d5_band_lo_hz")), float(g("d5_band_hi_hz"))),
                min_spindle_hz=float(g("d5_min_rpm")) / 60.0)


def dominant_nh(df, ch):
    """As r03_explore.dominant_nh: frequency of the dominant line where it lies off the harmonics (D5 settings)."""
    out_f = np.full(len(df), np.nan); out_i = np.full(len(df), np.nan); cache = {}
    for i, r in enumerate(df.itertuples()):
        key = (r.machine, int(r.run))
        if key not in cache:
            z = np.load(f"{QC.SEG}/{r.machine}_dataset{int(r.run)}.npz", allow_pickle=True)
            cache[key] = (z, json.loads(str(z["meta"])), z[f"{ch}_psd"] if f"{ch}_psd" in z.files else None)
        z, meta, PSD = cache[key]
        if PSD is None:
            continue
        j = int(np.argmin(np.abs(z["start"] - r.start) + np.abs(z["end"] - r.end)))
        P = PSD[j].astype(float); dfh = meta["channels"][ch]["df"]
        fr = np.arange(len(P)) * dfh; fsp = r.N / 60.0
        chk = (fr >= 50.0) & (fr < 12000.0)
        if fsp < BASE_IND["min_spindle_hz"] or not np.isfinite(P[chk]).all():
            continue
        inband = (fr >= BASE_IND["band"][0]) & (fr < BASE_IND["band"][1])
        harm = np.zeros(len(fr), bool)
        idx = np.round(np.arange(1, int(10000.0 / fsp) + 2) * fsp / dfh).astype(int)
        for dd in range(-BASE_IND["comb_half_bins"], BASE_IND["comb_half_bins"] + 1):
            ii = idx + dd; harm[ii[(ii >= 0) & (ii < len(fr))]] = True
        nh, hh = inband & ~harm, inband & harm
        if not nh.any() or not hh.any():
            continue
        on = P[nh].max() > P[hh].max()
        out_i[i] = float(on)
        if on:
            out_f[i] = fr[np.where(nh)[0][np.argmax(P[nh])]]
    return out_f, out_i


A3 = np.load(f"{P4}/outputs/r03_arrays.npz")
FQ = {}
for ch in ("acc0", "s2"):
    fq, ind = dominant_nh(paths, ch)
    assert np.array_equal(ind, A3[f"{ch}_base"], equal_nan=True), ch
    m1 = (paths.y == 1).values
    assert np.array_equal(fq[m1], A3["e9_fa" if ch == "acc0" else "e9_fb"], equal_nan=True), ch
    FQ[ch] = fq
    for lab in (0, 1):
        m = (ind == 1) & (y == lab)
        f = f"imi_vm20i, label {lab}, {CH_NAME[ch]} indicator on (D5 settings)"
        q25, q50, q75 = np.quantile(fq[m], [0.25, 0.5, 0.75])
        for nm, val in (("q25", q25), ("median", q50), ("q75", q75)):
            reg(f"ph_freq_{ch}_l{lab}_{nm}", float(val), "Hz", est(f, {"q25": "lower quartile", "median": "median", "q75": "upper quartile"}[nm] + " of the frequency of the dominant off-harmonic line"))
        reg(f"ph_freq_{ch}_l{lab}_n", int(m.sum()), "paths", est(f, "paths"))
    m = (ind == 1) & (y == 1)
    rho = spearmanr(fq[m], paths.N.values[m]).statistic
    reg(f"ph_freq_{ch}_l1_rho_rpm", float(rho), "rank correlation", est(f"imi_vm20i, label 1, {CH_NAME[ch]} indicator on (D5 settings)", "Spearman rank correlation between the frequency of the dominant off-harmonic line and the commanded spindle speed"))
    reg(f"ph_freq_{ch}_l1_rpm_levels", int(len(np.unique(paths.N.values[m]))), "speeds", est(f"imi_vm20i, label 1, {CH_NAME[ch]} indicator on (D5 settings)", "distinct commanded spindle speeds among these paths"))
log("DV17 frequencies")
# the two Fig. 2 examples (selected by the rule of the plan in r04_figures.py)
for L_ in (1, 0):
    j = np.where((paths.run.values == g("fig3_run")) & np.isclose(paths.start.values, g(f"fig3_l{L_}_start_s")))[0]
    assert len(j) == 1
    for ch in ("acc0", "s2"):
        assert A3[f"{ch}_base"][j[0]] == 1
        reg(f"ph_fig2_l{L_}_{ch}_hz", float(FQ[ch][j[0]]), "Hz", est(f"imi_vm20i, run fig3_run, the label {L_} example of Fig. 2, {CH_NAME[ch]}", "frequency of the dominant off-harmonic line"))
pd.DataFrame({"machine": paths.machine, "run": paths.run.astype(int), "start": paths.start, "end": paths.end,
              "acc0_line_hz": FQ["acc0"], "s2_line_hz": FQ["s2"]}).to_csv(f"{P4}/outputs/line_frequencies.csv", index=False)

# ------------------------------------------------------------------ Table 1 by tool family and workpiece material (review item R12)
MAT = {"AA 6061": "aa6061", "AISI 1018": "aisi1018", "AISI 4140": "aisi4140"}
for mach, mtag in (("imi_vm20i", "vm20i"), ("imi_vmx30ui", "vmx30ui")):
    d = QC.load_paths(mach)
    for (fam, mat), s in d.groupby(["family", "Workpiece"]):
        f = f"{mach}, labels 0 and 1, family {fam}, workpiece {mat}"
        reg(f"tab1_{mtag}_{fam.lower()}_{MAT[mat]}_runs", int(s.run.nunique()), "runs", est(f, "runs"))
        for lab in (0, 1):
            reg(f"tab1_{mtag}_{fam.lower()}_{MAT[mat]}_l{lab}", int((s.y == lab).sum()), "paths", est(f, f"label {lab} paths"))
    assert (d.groupby("run").Workpiece.nunique() == 1).all()
log("Table 1 counts")

# ------------------------------------------------------------------ times of the freeze and of the confirmatory run (the authors' own files; review item R9)
ft = datetime.datetime.strptime(V3["freeze_time"]["value"], "%Y-%m-%dT%H:%M:%SZ")
ct = datetime.datetime.fromtimestamp(os.path.getmtime(f"{P4}/outputs/r01_arrays.npz"), datetime.UTC).replace(tzinfo=None)
assert ct > ft
fmt = lambda t: f"{t:%H:%M} UTC on {t.day} {t:%B %Y}"
reg("facts_freeze_time_text", fmt(ft), "text", {"population": "documented facts", "filter": "values3.json freeze_time (q06_freeze.py)", "statistic": "time the SHA-256 digests of the frozen files were taken"}, kind="text")
reg("facts_confirm_time_text", fmt(ct), "text", {"population": "documented facts", "filter": "modification time of outputs/r01_arrays.npz", "statistic": "time the confirmatory script r01_confirm.py wrote its results"}, kind="text")
reg("facts_freeze_to_confirm_min", int(round((ct - ft).total_seconds() / 60)), "minutes", {"population": "documented facts", "filter": "the two times above", "statistic": "minutes from the freeze to the confirmatory results"})
V.save()
log("r11_review: registered " + str(sum(1 for r in V.v.values() if r["source"] == "r11_review.py")) + " values")