"""r03_explore.py: the exploratory analyses E1 to E11 of the plan (inputs/plan.json, "exploratory").

Estimates and family-stratified run-bootstrap intervals only; no p-value is computed or registered. The
bootstrap weights are the ones of the confirmatory analysis (q_confirm.boot_weights on the same population,
seed and number of resamples), and every statistic uses the frozen functions of q_common.py or the
vectorised helpers of r_fast.py, which are checked against the frozen functions before use.
"""
import sys, json, time
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 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="r03_explore.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 = int(g("boot_b")), int(g("seed_analysis")), float(g("ci_level"))
REP, K, C = int(g("s1_repeats")), int(g("s1_folds")), float(g("lr_c"))
EDGE = float(g("speed_class_edge_rpm"))
METHOD = "percentile interval at ci_level over boot_b run resamples drawn within tool families (q_confirm.boot_weights), predictions held fixed; 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 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
fam = paths.family.values
W = QF.boot_weights(paths, B, SEED)
arr = np.load(f"{P4}/outputs/r01_arrays.npz", allow_pickle=True)
assert (arr["h1_run_id"] == paths.run_id.values.astype(str)).all() and (arr["h2_run_id"] == paths.run_id.values.astype(str)).all()
RF.check_against_frozen(y, arr["h1_scores_family"], W, fam)
BASEF = "imi_vm20i, labels 0 and 1"


# ------------------------------------------------------------------ evaluation schemes beyond the frozen ones
def oof_random_matrix(df, fset, repeats=REP, k=K, seed=SEED):
    """Random path-level folds stratified by label (the Phase 2 screen's random scheme), balanced model."""
    X = QC.design(df, fset); yy = df.y.values
    rng = np.random.default_rng(seed)
    S = np.full((repeats, len(df)), np.nan)
    idx = np.arange(len(df))
    for rep in range(repeats):
        folds = np.zeros(len(df), int)
        for lab in (0, 1):
            ii = idx[yy == lab].copy(); rng.shuffle(ii)
            folds[ii] = np.arange(len(ii)) % k
        for f in range(k):
            tr, te = folds != f, folds == f
            S[rep, te] = QC.fit_predict(X[tr], yy[tr], X[te], "lr", C, False, True)
    return S


def oof_group(df, fset, groups):
    """A fold per group (E7: a fold per workpiece material), balanced model."""
    X = QC.design(df, fset); yy = df.y.values; groups = np.asarray(groups)
    s = np.full(len(df), np.nan)
    for grp in sorted(set(groups)):
        te = groups == grp; tr = ~te
        if len(np.unique(yy[tr])) < 2:
            raise RuntimeError(f"training set without {grp} lacks a label")
        s[te] = QC.fit_predict(X[tr], yy[tr], X[te], "lr", C, False, True)
    return s


FSET_NAME = {"highrate": "log RMS of Acc0, Acc1, Sound0, Sound1 and Sound2", "mtc": "the six MTConnect spindle features",
             "highrate_no_s0": "log RMS of Acc0, Acc1, Sound1 and Sound2 (Sound0 left out)", "clip": "share of Sound0 samples at full scale",
             "meta": "design variables (log spindle speed, depths of cut and feed per tooth, cutting direction, workpiece) with tool family",
             "meta_no_family": "design variables (log spindle speed, depths of cut and feed per tooth, cutting direction, workpiece) without tool family"}
TAG = {"highrate": "highrate", "mtc": "mtc", "highrate_no_s0": "nos0", "clip": "clip", "meta": "meta", "meta_no_family": "metanf"}


def reg_auc(name, auc, draws, fset, scheme):
    e = est(f"{BASEF}, feature set {TAG[fset]} ({FSET_NAME[fset]})",
            f"pooled out-of-fold AUC, logistic model with balanced class weights, {scheme}")
    reg(name, auc, "AUC", e)
    reg_ci(f"{name}_ci", name, draws, e)


RUN_S = "run-held-out scheme (mean over repeats of each repeat's pooled AUC, s1_repeats repeats of s1_folds folds of whole runs)"
FAM_S = "tool-family-held-out scheme (a fold per family)"
RAND_S = "random path-level scheme (s1_repeats repeats of s1_folds folds stratified by label; mean over repeats of each repeat's pooled AUC)"
MAT_S = "workpiece-held-out scheme (a fold per workpiece material)"

# ------------------------------------------------------------------ E2: without Sound0, and the clipping-only score
for fset in ("highrate_no_s0", "clip"):
    res, sc = QF.auc_drop(paths, fset, B, SEED, LEVEL, REP, K, "lr", C, False, W=W)
    reg_auc(f"e2_{TAG[fset]}_auc_run", res["auc_run"], RF.auc_rep_draws(y, sc["scores_run"], W), fset, RUN_S)
    reg_auc(f"e2_{TAG[fset]}_auc_family", res["auc_family"], RF.auc_draws(y, sc["scores_family"], W), fset, FAM_S)
    e = est(f"{BASEF}, feature set {TAG[fset]} ({FSET_NAME[fset]})", "run-held-out AUC minus tool-family-held-out AUC, exploratory")
    reg(f"e2_{TAG[fset]}_drop", res["drop"], "AUC", e)
    V.interval(f"e2_{TAG[fset]}_drop_ci", f"e2_{TAG[fset]}_drop", res["drop_lo"], res["drop_hi"], LEVEL, METHOD, e)
    assert abs(res["auc_run"] - QC.auc_over_repeats(y, sc["scores_run"])) < 1e-12
    log(f"E2 {fset}: run {res['auc_run']:.3f} family {res['auc_family']:.3f}")

# ------------------------------------------------------------------ E3: design variables alone (run-held-out), and E4 family scheme for the set without family
for fset in ("meta", "meta_no_family"):
    S = QC.oof_run_matrix(paths, fset, REP, K, "lr", C, False, True, SEED)
    reg_auc(f"e3_{TAG[fset]}_auc_run", QC.auc_over_repeats(y, S), RF.auc_rep_draws(y, S, W), fset, RUN_S)
    log(f"E3 {fset}: run {QC.auc_over_repeats(y, S):.3f}")
s = QC.oof_family(paths, "meta_no_family", "lr", C, False, True)
reg_auc("e4_metanf_auc_family", QC.auc_w(y, s), RF.auc_draws(y, s, W), "meta_no_family", FAM_S)

# ------------------------------------------------------------------ E4: random path-level scheme for every feature set
for fset in ("highrate", "mtc", "highrate_no_s0", "clip", "meta", "meta_no_family"):
    S = oof_random_matrix(paths, fset)
    reg_auc(f"e4_{TAG[fset]}_auc_random", QC.auc_over_repeats(y, S), RF.auc_rep_draws(y, S, W), fset, RAND_S)
    log(f"E4 {fset}: random {QC.auc_over_repeats(y, S):.3f}")

# ------------------------------------------------------------------ E7: a fold per workpiece material
reg("e7_n_materials", int(paths.Workpiece.nunique()), "materials", est(BASEF, "workpiece materials (folds of the workpiece-held-out scheme)"))
for fset in ("highrate", "mtc", "highrate_no_s0", "clip"):
    s = oof_group(paths, fset, paths.Workpiece.values)
    reg_auc(f"e7_{TAG[fset]}_auc_material", QC.auc_w(y, s), RF.auc_draws(y, s, W), fset, MAT_S)
    log(f"E7 {fset}: material {QC.auc_w(y, s):.3f}")

# ------------------------------------------------------------------ E10: pairs within and across families, and the within-family AUC of A and E
qual = [f.strip() for f in g("pop_within_families_primary_list").split(",")]
for H in ("h1", "h2"):
    Srun, sfam = arr[f"{H}_scores_run"], arr[f"{H}_scores_family"]
    parts_run = [RF.within_across(y, srow, fam, W) for srow in Srun]
    pooled_r, within_r, across_r = (np.mean([p[i] for p in parts_run], axis=0) for i in range(3))
    pooled_f, within_f, across_f, share_f = RF.within_across(y, sfam, fam, W)
    pt_run = [RF.within_across(y, srow, fam, np.ones((1, len(y)))) for srow in Srun]
    pt_fam = RF.within_across(y, sfam, fam, np.ones((1, len(y))))
    for scheme, pt, dr in (("run", [np.mean([p[i][0] for p in pt_run]) for i in range(4)], (within_r, across_r)),
                           ("family", [x[0] for x in pt_fam], (within_f, across_f))):
        for j, part in enumerate(("within", "across")):
            e = est(f"{BASEF}, {'highrate' if H == 'h1' else 'mtc'} features, primary model of {H.upper()}",
                    f"AUC over the label 1 and label 0 pairs {'inside one tool family' if part == 'within' else 'from different tool families'}, {RUN_S if scheme == 'run' else FAM_S}")
            reg(f"e10_{H}_{scheme}_{part}", pt[1 + j], "AUC", e)
            reg_ci(f"e10_{H}_{scheme}_{part}_ci", f"e10_{H}_{scheme}_{part}", dr[j], e)
        assert abs(pt[0] - (V4[f"{H}_auc_run"]["value"] if scheme == "run" else V4[f"{H}_auc_family"]["value"])) < 1e-9
    if H == "h1":
        reg("e10_share_across", pt_fam[3][0], "fraction", est(BASEF, "share of the label 1 and label 0 pairs whose paths come from different tool families"))
    for f in qual:
        m = fam == f
        for scheme in ("run", "family"):
            if scheme == "run":
                val = float(np.mean([QC.auc_w(y[m], srow[m]) for srow in Srun]))
                dr = RF.auc_rep_draws(y[m], Srun[:, m], W[:, m])
            else:
                val = QC.auc_w(y[m], sfam[m])
                dr = RF.auc_draws(y[m], sfam[m], W[:, m])
            e = est(f"imi_vm20i, labels 0 and 1, family {f}, {'highrate' if H == 'h1' else 'mtc'} features, primary model of {H.upper()}",
                    f"AUC inside family {f}, {RUN_S if scheme == 'run' else FAM_S + ' (the family is held out)'}")
            reg(f"e10_{H}_{scheme}_{f.lower()}", val, "AUC", e)
            reg_ci(f"e10_{H}_{scheme}_{f.lower()}_ci", f"e10_{H}_{scheme}_{f.lower()}", dr, e)
    log(f"E10 {H} done")

# ------------------------------------------------------------------ E5: off-harmonic shares by label and channel across comb and band settings
SETTINGS = {"base": {}, "comb1": {"comb_half_bins": int(g("comb_bins_v6a"))}, "comb3": {"comb_half_bins": int(g("comb_bins_v6b"))},
            "rel": {"rel_speed_tol": float(g("rel_speed_tol"))}, "band5k": {"band": (float(g("d5_band_lo_hz")), float(g("band_hi_v7_hz")))}}
SET_DESC = {"base": "D5 settings", "comb1": "comb of comb_bins_v6a bins (V6a)", "comb3": "comb of comb_bins_v6b bins (V6b)",
            "rel": "relative comb, rel_speed_tol (V6c)", "band5k": "band up to band_hi_v7_hz (V7)"}
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, rel_speed_tol=0.0)
IND = {}
for ch in ("acc0", "acc1", "s1", "s2"):
    for sname, chg in SETTINGS.items():
        IND[(ch, sname)] = QC.nh_indicator(paths, ch, **{**BASE_IND, **chg})
    log(f"E5 indicators {ch}")
for ch in ("acc0", "acc1", "s1", "s2"):
    base = IND[(ch, "base")]
    ok = np.isfinite(base)
    assert ok.sum() == V2[f"d5_n_imi_vm20i_{ch}"]["value"], ch
    for lab in (0, 1):
        m = ok & (y == lab)
        assert abs(base[m].mean() - V2[f"d5_nhdom_imi_vm20i_{ch}_l{lab}"]["value"]) < 1e-12, (ch, lab)
    for sname in SETTINGS:
        ind = IND[(ch, sname)]
        for lab in (0, 1):
            m = np.isfinite(ind) & (y == lab)
            e = est(f"imi_vm20i, label {lab}, spindle speed at least d5_min_rpm, {ch} spectrum", f"share of paths whose dominant line lies off the spindle harmonics, {SET_DESC[sname]}")
            name = f"e5_{ch}_{sname}_l{lab}"
            reg(name, float(ind[m].mean()), "fraction", e)
            reg_ci(f"{name}_ci", name, RF.share_draws(np.where(m, ind, 0.0), W * m), e)
        if sname != "base":
            both = np.isfinite(ind) & ok
            reg(f"e5_{ch}_{sname}_changed", float((ind[both] != base[both]).mean()), "fraction",
                est(f"imi_vm20i, labels 0 and 1, spindle speed at least d5_min_rpm, {ch} spectrum", f"share of paths whose indicator under the {SET_DESC[sname]} differs from the D5 settings"))
reg("e5_n_paths", int(np.isfinite(IND[("acc0", "base")]).sum()), "paths", est("imi_vm20i, labels 0 and 1, spindle speed at least d5_min_rpm", "paths with an indicator"))

# ------------------------------------------------------------------ E6: label 0 paths with both lines, label 1 paths with neither
a, b = IND[("acc0", "base")], IND[("s2", "base")]
okab = np.isfinite(a) & np.isfinite(b)
for lab, combo, cond in ((0, "both", (a == 1) & (b == 1)), (1, "neither", (a == 0) & (b == 0))):
    m = okab & (y == lab)
    hit = m & cond
    f = f"imi_vm20i, label {lab}, spindle speed at least d5_min_rpm, Acc0 and Sound2 spectra"
    stat = "Acc0 and Sound2 indicators both on" if combo == "both" else "Acc0 and Sound2 indicators both off"
    reg(f"e6_l{lab}_{combo}_n", int(hit.sum()), "paths", est(f, f"paths with the {stat} (D5 settings)"))
    reg(f"e6_l{lab}_n", int(m.sum()), "paths", est(f, "paths with both indicators defined"))
    e = est(f, f"share of paths with the {stat} (D5 settings)")
    reg(f"e6_l{lab}_{combo}_share", float(hit.sum() / m.sum()), "fraction", e)
    reg_ci(f"e6_l{lab}_{combo}_share_ci", f"e6_l{lab}_{combo}_share", RF.share_draws(np.where(m, cond, 0).astype(float), W * m), e)
    reg(f"e6_l{lab}_{combo}_runs", int(paths[hit].run.nunique()), "runs", est(f, f"runs holding paths with the {stat}"))
    for fm in sorted(set(fam)):
        reg(f"e6_l{lab}_{combo}_fam_{fm.lower()}", int((hit & (fam == fm)).sum()), "paths", est(f"{f}, family {fm}", f"paths with the {stat}"))
    top = paths[hit].groupby("run_id").size().sort_values(ascending=False)
    reg(f"e6_l{lab}_{combo}_top_run_n", int(top.iloc[0]) if len(top) else 0, "paths", est(f, f"paths with the {stat} in the run that holds most of them"))
    reg(f"e6_l{lab}_{combo}_top_run", str(top.index[0]).replace("imi_vm20i:", "run ") if len(top) else "none", "text", est(f, "that run"), kind="text")
log("E6 done")


# ------------------------------------------------------------------ E9: do the two lines share a frequency?
def dominant_nh(df, ch):
    """Frequency of the dominant line where it lies off the harmonics (D5 settings), else NaN; same bins as
    q_common.nh_indicator, and the on and off decision is asserted equal to it below."""
    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


lab1 = paths[paths.y == 1].reset_index(drop=True)
m1 = (paths.y == 1).values
fa, ia = dominant_nh(lab1, "acc0")
fb, ib = dominant_nh(lab1, "s2")
assert np.array_equal(ia, a[m1], equal_nan=True) and np.array_equal(ib, b[m1], equal_nan=True)
both = (ia == 1) & (ib == 1)
TOL = float(g("freq_match_tol_hz"))
match = both & (np.abs(fa - fb) <= TOL)
W1 = W[:, m1]
f9 = "imi_vm20i, label 1, Acc0 and Sound2 indicators both on (D5 settings)"
reg("e9_n_both", int(both.sum()), "paths", est(f9, "paths"))
e = est(f9, "share of paths whose Acc0 and Sound2 dominant lines lie within freq_match_tol_hz of each other")
reg("e9_share_match", float(match.sum() / both.sum()), "fraction", e)
reg_ci("e9_share_match_ci", "e9_share_match", RF.share_draws(match[both].astype(float), W1[:, both]), e)
runs1 = lab1.run_id.values
num_r, den_r = {}, {}
for r in np.unique(runs1[both]):
    ii = np.where(both & (runs1 == r))[0]
    if len(ii) < 2:
        continue
    D = np.abs(fa[ii][:, None] - fb[ii][None, :]) <= TOL
    np.fill_diagonal(D, False)
    num_r[r], den_r[r] = D.sum(), len(ii) * (len(ii) - 1)
rw = {r: W1[:, np.where(runs1 == r)[0][0]] for r in num_r}
bnum = sum(rw[r] * num_r[r] for r in num_r); bden = sum(rw[r] * den_r[r] for r in num_r)
e = est(f"{f9}, pairs of different paths from the same run", "share of pairs whose Acc0 line of one path lies within freq_match_tol_hz of the Sound2 line of the other")
reg("e9_share_baseline", float(sum(num_r.values()) / sum(den_r.values())), "fraction", e)
reg_ci("e9_share_baseline_ci", "e9_share_baseline", RF.ratio(bnum, bden), e)
reg("e9_n_pairs_baseline", int(sum(den_r.values())), "pairs", est(f"{f9}, pairs of different paths from the same run", "ordered pairs"))
log("E9 done")

# ------------------------------------------------------------------ E11: Acc0 line share in family A against family E (label 1)
ae = (lab1.family.isin(["A", "E"])).values & np.isfinite(ia)
d11 = lab1[ae].reset_index(drop=True)
i11 = ia[ae]
W11 = QF.boot_weights(d11, B, SEED)
fa11 = (d11.family == "A").values
f11 = "imi_vm20i, label 1, families A and E, Acc0 spectrum"
sa = RF.share_draws(np.where(fa11, i11, 0.0), W11 * fa11)
se = RF.share_draws(np.where(~fa11, i11, 0.0), W11 * ~fa11)
for tag, val, dr, what in (("a", i11[fa11].mean(), sa, "family A"), ("e", i11[~fa11].mean(), se, "family E")):
    e = est(f"{f11}, {what}", "share of paths whose dominant Acc0 line lies off the spindle harmonics (D5 settings)")
    reg(f"e11_share_{tag}", float(val), "fraction", e)
    reg_ci(f"e11_share_{tag}_ci", f"e11_share_{tag}", dr, e)
    sub = d11[fa11 if tag == "a" else ~fa11]
    reg(f"e11_n_{tag}", int(len(sub)), "paths", est(f"{f11}, {what}", "paths"))
    reg(f"e11_runs_{tag}", int(sub.run.nunique()), "runs", est(f"{f11}, {what}", "runs"))
    reg(f"e11_rpm_median_{tag}", float(sub.N.median()), "r/min", est(f"{f11}, {what}", "median commanded spindle speed"))
    reg(f"e11_rpm_hi_share_{tag}", float((sub.N >= EDGE).mean()), "fraction", est(f"{f11}, {what}", "share of paths at or above speed_class_edge_rpm"))
e = est(f11, "share in family A minus share in family E (runs resampled within each family)")
reg("e11_diff", float(i11[fa11].mean() - i11[~fa11].mean()), "fraction", e)
reg_ci("e11_diff_ci", "e11_diff", sa - se, e)
log("E11 done")

# ------------------------------------------------------------------ E1: clipping and how the Sound0 level flattens
FS = 2 ** (int(V2["wav_bits_imi_vm20i_s0"]["value"]) - 1) if "wav_bits_imi_vm20i_s0" in V2 else 2 ** 15
reg("e1_fullscale_counts", FS, "counts", {"population": "definition", "filter": "16-bit sound files", "statistic": "full-scale magnitude of a 16-bit sample"})
MID = V.set("e1_bin_edge_mid", 0.1, "fraction", {"population": "analyst judgment", "filter": "E1 presentation", "statistic": "middle edge of the full-scale share bins"}, kind="judgment")
HIGH = V.set("e1_bin_edge_high", 0.3, "fraction", {"population": "analyst judgment", "filter": "E1 presentation", "statistic": "upper edge of the full-scale share bins"}, kind="judgment")
EDGES = [0.0, float(g("d9_threshold")), MID, HIGH, 1.0]
for mach in ("imi_vm20i", "imi_vmx30ui"):
    Pm = QC.load_paths(mach)
    for lab in (0, 1):
        s = Pm[Pm.y == lab]
        f1 = f"{mach}, label {lab}, Sound0"
        reg(f"e1_{mach}_l{lab}_fs_median", float(s.s0_fullscale.median()), "fraction", est(f1, "median share of samples at full scale"))
        reg(f"e1_{mach}_l{lab}_any_clip", float((s.s0_fullscale > 0).mean()), "fraction", est(f1, "share of paths with at least one sample at full scale"))
        reg(f"e1_{mach}_l{lab}_rmsfrac_median", float((s.s0_rms / FS).median()), "fraction", est(f1, "median RMS as a fraction of e1_fullscale_counts"))
    for k in range(len(EDGES) - 1):
        lo, hi = EDGES[k], EDGES[k + 1]
        m = (Pm.s0_fullscale > lo) & (Pm.s0_fullscale <= hi) if k else (Pm.s0_fullscale <= hi)
        f1 = f"{mach}, labels 0 and 1, Sound0 share of samples at full scale in bin {k + 1} of the E1 bins"
        reg(f"e1_{mach}_bin{k + 1}_n", int(m.sum()), "paths", est(f1, "paths"))
        if m.sum():
            reg(f"e1_{mach}_bin{k + 1}_rmsfrac_median", float((Pm.s0_rms[m] / FS).median()), "fraction", est(f1, "median RMS as a fraction of e1_fullscale_counts"))
            reg(f"e1_{mach}_bin{k + 1}_l1_share", float(Pm.y[m].mean()), "fraction", est(f1, "share of label 1 paths"))
log("E1 done")

# ------------------------------------------------------------------ E8: imi_vmx30ui, the MTConnect model of imi_vm20i transferred
vmx = QC.load_paths("imi_vmx30ui").dropna(subset=QC.FEATS["mtc"]).reset_index(drop=True)
s8 = QC.fit_predict(QC.design(paths, "mtc"), y, QC.design(vmx, "mtc"), "lr", C, False, True)
W8 = QF.boot_weights(vmx, B, SEED)
f8 = "imi_vmx30ui, labels 0 and 1, MTConnect features; model fitted on every imi_vm20i path (labels 0 and 1)"
e = est(f8, "pooled AUC of the transferred balanced logistic model on the six MTConnect spindle features")
reg("e8_mtc_transfer_auc", QC.auc_w(vmx.y.values, s8), "AUC", e)
reg_ci("e8_mtc_transfer_auc_ci", "e8_mtc_transfer_auc", RF.auc_draws(vmx.y.values, s8, W8), e)
reg("e8_n_paths", int(len(vmx)), "paths", est(f8, "paths"))
reg("e8_n_l1", int(vmx.y.sum()), "paths", est(f8, "label 1 paths"))
reg("e8_n_runs", int(vmx.run.nunique()), "runs", est(f8, "runs"))
reg("e8_n_runs_l1", int(vmx[vmx.y == 1].run.nunique()), "runs", est(f8, "runs holding label 1"))
log("E8 done")

# ------------------------------------------------------------------ per-path indicator table (Online Resource ESM_3) and arrays for the figures
tab = paths[["machine", "run", "start", "end", "label", "Tool", "family", "Workpiece", "N", "dur"]].copy()
for ch in ("acc0", "acc1", "s1", "s2"):
    tab[f"{ch}_offharmonic"] = IND[(ch, "base")]
tab.to_csv(f"{P4}/outputs/per_path_indicators.csv", index=False)
np.savez_compressed(f"{P4}/outputs/r03_arrays.npz", **{f"{ch}_{s_}": v for (ch, s_), v in IND.items()},
                    e9_fa=fa, e9_fb=fb, e9_both=both)
V.save()
log(f"r03_explore: registered {sum(1 for r in V.v.values() if r['source'] == 'r03_explore.py')} values")