"""p02_profile.py: data quality profile (reported before any finding).

Labels and their structure (D1), segment timing, parameter plausibility against the article's tool
table, MTConnect consistency with label.csv (D6), full-scale sound (D9), stream alignment (D8),
encoding, and the clustering of labels within runs. Writes outputs/profile.json and values.
"""
import json, glob, os, re, sys
import numpy as np, pandas as pd
sys.path.insert(0, "/home/claude/mill/code")
from vals import Values

RAW = "/home/claude/mill/raw"
DER = f"{RAW}/Phase2/derived"
OUT = "/home/claude/mill/outputs"
MACH = ("imi_vm20i", "imi_vmx30ui", "tmf_vf10")
LAB = ("imi_vm20i", "imi_vmx30ui")
V = Values(script="p02_profile.py")
T = pd.read_parquet(f"{OUT}/segments.parquet")
S = json.load(open("/home/claude/mill/inputs/sources.json"))
src = {s["id"]: s for s in S["sources"]}
prof = {}


def E(pop, filt, stat):
    return {"population": pop, "filter": filt, "statistic": stat}


SEGPOP = "labelled cutting paths (label.csv rows)"

# ---------------------------------------------------------------- D1 labels ------------------------------
V.set("n_paths", len(T), "paths", E(SEGPOP, "all machines", "count"))
tab = {}
for mach in MACH:
    sub = T[T.machine == mach]
    V.set(f"n_paths_{mach}", len(sub), "paths", E(SEGPOP, mach, "count"))
    tot = sub.dur.sum()
    for lab in (0, 1, 2):
        s2 = sub[sub.label == lab]
        V.set(f"n_paths_{mach}_l{lab}", len(s2), "paths", E(SEGPOP, f"{mach}, label {lab}", "count"))
        V.set(f"time_share_{mach}_l{lab}", float(s2.dur.sum() / tot), "fraction", E(SEGPOP, f"{mach}", f"share of cutting time with label {lab}"))
        V.set(f"path_share_{mach}_l{lab}", float(len(s2) / len(sub)), "fraction", E(SEGPOP, f"{mach}", f"share of paths with label {lab}"))
        tab[(mach, lab)] = float(s2.dur.sum() / tot)
dev = max(abs(round(100 * tab[(m, l)], 2) - 100 * V.get(f"doc_t4_{m}_l{l}")) for m in MACH for l in (0, 1, 2))
V.set("t4_max_abs_dev_pp", float(dev), "percentage points", E(SEGPOP, "machine by label cells", "max absolute difference between time shares rounded to 2 decimals and Table 4"))
V.set("n_paths_l1", int((T.label == 1).sum()), "paths", E(SEGPOP, "all machines, label 1", "count"))
V.set("n_paths_l2_nondefective_tool", int(((T.label == 2) & ~T.defective).sum()), "paths", E(SEGPOP, "label 2 on a tool without asterisk", "count"))
V.set("n_paths_defective_l1", int(((T.label == 1) & T.defective).sum()), "paths", E(SEGPOP, "asterisk tools labelled 1", "count"))
V.set("n_paths_defective", int(T.defective.sum()), "paths", E(SEGPOP, "asterisk tools", "count"))

# labels by material and tool (lab machines)
L = T[T.machine.isin(LAB)]
bymat = L.groupby("Workpiece").apply(lambda g: pd.Series({"paths": len(g), "l1_time_share": g.dur[g.label == 1].sum() / g.dur.sum(),
                                                         "l1_path_share": (g.label == 1).mean()}), include_groups=False)
prof["lab_by_material"] = bymat.reset_index().to_dict("records")
for wp, r in bymat.iterrows():
    key = re.sub(r"[^a-z0-9]", "", wp.lower())
    V.set(f"lab_l1_path_share_{key}", float(r.l1_path_share), "fraction", E(SEGPOP, f"laboratory machines, workpiece {wp}", "share of paths with label 1"))
    V.set(f"lab_n_paths_{key}", int(r.paths), "paths", E(SEGPOP, f"laboratory machines, workpiece {wp}", "count"))
bytool = L.groupby("Tool").apply(lambda g: pd.Series({"paths": len(g), "runs": g.run_id.nunique(), "l0": (g.label == 0).mean(), "l1": (g.label == 1).mean(), "l2": (g.label == 2).mean()}), include_groups=False)
prof["lab_by_tool"] = bytool.reset_index().to_dict("records")
for tool, r in bytool.iterrows():
    key = tool.lower().replace("*", "star")
    et = E(SEGPOP, f"laboratory machines, tool {tool}", "count")
    V.set(f"d1_paths_tool_{key}", int(r.paths), "paths", et)
    V.set(f"d1_runs_tool_{key}", int(r.runs), "runs", {**et, "statistic": "runs"})
    V.set(f"d1_l1_share_tool_{key}", float(r.l1), "fraction", {**et, "statistic": "share of paths with label 1"})
    V.set(f"d1_l2_share_tool_{key}", float(r.l2), "fraction", {**et, "statistic": "share of paths with label 2"})
    g = L[L.Tool == tool]
    V.set(f"d1_l1_time_share_tool_{key}", float(g.dur[g.label == 1].sum() / g.dur.sum()), "fraction", {**et, "statistic": "share of cutting time with label 1"})
isA = L.Tool.str.replace("*", "", regex=False) == "A"
l1 = L.label == 1
alu = L.Workpiece == "AA 6061"
V.set("d1_toola_share_of_lab_l1", float((isA & l1).sum() / l1.sum()), "fraction", E(SEGPOP, "laboratory machines, label 1", "share of label 1 paths cut with tool A or A*"))
V.set("d1_toola_share_of_lab_l1_time", float(L.dur[isA & l1].sum() / L.dur[l1].sum()), "fraction", E(SEGPOP, "laboratory machines, label 1", "share of label 1 cutting time with tool A or A*"))
V.set("d1_n_lab_l1", int(l1.sum()), "paths", E(SEGPOP, "laboratory machines, label 1", "count"))
V.set("d1_l1_paths_aa6061_toola", int((alu & isA & l1).sum()), "paths", E(SEGPOP, "laboratory machines, AA 6061, tool A or A*", "label 1 paths"))
V.set("d1_paths_aa6061_toola", int((alu & isA).sum()), "paths", E(SEGPOP, "laboratory machines, AA 6061, tool A or A*", "paths"))
V.set("d1_l1_paths_aa6061_other", int((alu & ~isA & l1).sum()), "paths", E(SEGPOP, "laboratory machines, AA 6061, tools other than A and A*", "label 1 paths"))
V.set("d1_paths_aa6061_other", int((alu & ~isA).sum()), "paths", E(SEGPOP, "laboratory machines, AA 6061, tools other than A and A*", "paths"))
V.set("d1_l1_paths_aa6061", int((alu & l1).sum()), "paths", E(SEGPOP, "laboratory machines, AA 6061", "label 1 paths"))
V.set("d1_toola_paths_on_aa6061_share", float((alu & isA).sum() / isA.sum()), "fraction", E(SEGPOP, "laboratory machines, tool A or A*", "share of paths cut in AA 6061"))
V.set("lab_tools_per_run_max", int(L.groupby("run_id").Tool.nunique().max()), "tools", E("laboratory runs", "label.csv", "most distinct tools in one run"))
VF = T[T.machine == "tmf_vf10"]
V.set("vf10_tools_per_run_min", int(VF.groupby("run").Tool.nunique().min()), "tools", E("tmf_vf10 runs", "label.csv", "fewest distinct tools in one run"))
V.set("vf10_tools_per_run_max", int(VF.groupby("run").Tool.nunique().max()), "tools", E("tmf_vf10 runs", "label.csv", "most distinct tools in one run"))
V.set("vf10_max_f_per_tool_run", int(VF.groupby(["run", "Tool"]).F.nunique().max()), "feeds", E("tmf_vf10 runs", "label.csv", "most distinct F values for one tool within one run"))
INVR = pd.read_csv(f"{OUT}/inventory_runs.csv")
rem = INVR.set_index(["machine", "run"]).summary_remark.fillna("")
V.set("vf10_runs_feed_turned_down", int(rem.loc["tmf_vf10"].str.contains("Feed rate turned down", case=False).sum()), "runs",
      E("dataset_summary.xlsx", "tmf_vf10", "runs whose remark records a feed reduction"))

# clustering of label 1 within runs
runs = T.groupby(["machine", "run"]).apply(lambda g: pd.Series({"paths": len(g), "l1_paths": int((g.label == 1).sum()),
                                                               "l1_time": g.dur[g.label == 1].sum(), "time": g.dur.sum()}), include_groups=False).reset_index()
COVER = 0.5  # share of all label 1 cutting time that the fewest runs must hold
V.set("l1_time_cover_share", COVER, "fraction", E("runs", "definition", "share of label 1 cutting time used for n_runs_half_l1_time_*"))
for mach in MACH:
    r = runs[runs.machine == mach]
    V.set(f"n_runs_any_l1_{mach}", int((r.l1_paths > 0).sum()), "runs", E("runs", mach, "runs with at least one label 1 path"))
    tot = r.l1_time.sum()
    srt = np.sort(r.l1_time.values)[::-1]
    k50 = int(np.searchsorted(np.cumsum(srt), COVER * tot) + 1) if tot > 0 else 0
    V.set(f"n_runs_half_l1_time_{mach}", k50, "runs", E("runs", mach, "fewest runs holding at least the cover share (l1_time_cover_share) of all label 1 cutting time"))
    if tot > 0:
        V.set(f"l1_time_share_top_runs_{mach}", float(np.cumsum(srt)[k50 - 1] / tot), "fraction", E("runs", mach, "share of all label 1 cutting time held by those fewest runs"))
    V.set(f"max_run_l1_path_share_{mach}", float((r.l1_paths / r.paths).max()), "fraction", E("runs", mach, "largest within-run share of label 1 paths"))
V.set("n_runs_all_l1", int((runs.l1_paths == runs.paths).sum()), "runs", E("runs", "all machines", "runs in which every path is label 1"))
_inv = pd.read_csv(f"{OUT}/inventory_runs.csv").set_index(["machine", "run"]).summary_remark.fillna("")
_all1 = runs[runs.l1_paths == runs.paths]
V.set("n_runs_all_l1_early_stopped", int(sum("early stop" in _inv.loc[(m_, r_)].lower() for m_, r_ in zip(_all1.machine, _all1.run))), "runs",
      E("runs", "all machines, runs in which every path is label 1", "runs whose summary remark records an early stop"))
V.set("n_runs_all_l1_lab", int((_all1.machine != "tmf_vf10").sum()), "runs", E("runs", "laboratory machines", "runs in which every path is label 1"))
V.set("n_runs_any_l1", int((runs.l1_paths > 0).sum()), "runs", E("runs", "all machines", "runs with at least one label 1 path"))
prof["runs"] = runs.to_dict("records")

# ---------------------------------------------------------------- segment timing ------------------------------
ov = 0; dups = 0; unsorted = 0; gaps = []; gap_mach = []; gap_labchg = []
for (m, r), g in T.groupby(["machine", "run"]):
    st, en = g.start.values, g.end.values
    if not (np.diff(st) >= 0).all():
        unsorted += 1
    o = np.argsort(st, kind="stable")
    st, en, lb = st[o], en[o], g.label.values[o]
    ov += int((st[1:] < en[:-1] - 1e-9).sum())
    dups += int(pd.DataFrame({"s": st, "e": en}).duplicated().sum())
    gaps.extend((st[1:] - en[:-1]).tolist())
    gap_mach.extend([m] * (len(st) - 1))
    gap_labchg.extend((lb[1:] != lb[:-1]).tolist())
gaps = np.array(gaps); gap_mach = np.array(gap_mach); gap_labchg = np.array(gap_labchg)
zg = np.abs(gaps) < 1e-9
V.set("n_zero_gaps_tmf_vf10", int((zg & (gap_mach == "tmf_vf10")).sum()), "pairs", E(SEGPOP, "consecutive paths within a tmf_vf10 run", "pairs with no gap"))
V.set("n_zero_gaps_label_change", int((zg & gap_labchg).sum()), "pairs", E(SEGPOP, "consecutive paths within a run", "pairs with no gap whose two labels differ"))
V.set("n_path_overlaps", ov, "pairs", E(SEGPOP, "consecutive paths within a run", "pairs where a path starts before the previous one ends"))
V.set("n_path_duplicates", dups, "paths", E(SEGPOP, "within a run", "rows repeating an earlier start and end"))
V.set("n_runs_unsorted", unsorted, "runs", E("runs", "label.csv", "runs whose rows are not in start order"))
V.set("n_zero_gaps", int((np.abs(gaps) < 1e-9).sum()), "pairs", E(SEGPOP, "consecutive paths within a run", "pairs with no gap (end equals next start)"))
V.set("n_nonpositive_dur", int((T.dur <= 0).sum()), "paths", E(SEGPOP, "all", "paths with end not after start"))
for mach in MACH:
    d = T[T.machine == mach].dur
    V.set(f"dur_median_{mach}", float(d.median()), "s", E(SEGPOP, mach, "median path duration"))
    V.set(f"dur_min_{mach}", float(d.min()), "s", E(SEGPOP, mach, "min path duration"))
    V.set(f"dur_max_{mach}", float(d.max()), "s", E(SEGPOP, mach, "max path duration"))

# ---------------------------------------------------------------- parameters against the tool table ---------------
t2 = {r[0].replace("*", ""): r[1] for r in src["S_table2"]["table"]["rows"]}
mism_Z = mism_D = 0
tooltab = {}
for tb, desc in t2.items():
    D = float(re.search(r"([\d.]+) mm diameter", desc).group(1))
    Z = int(re.search(r"(\d+)-(?:flute|insert)", desc).group(1))
    tooltab[tb] = (D, Z)
for tb, g in T.groupby("tool_base"):
    D, Z = tooltab[tb]
    mism_Z += int((g.Z != Z).sum()); mism_D += int((np.abs(g.D - D) > 1e-6).sum())
V.set("n_paths_z_mismatch_table2", mism_Z, "paths", E(SEGPOP, "all", "paths whose flute count differs from the article's Table 2"))
V.set("n_paths_d_mismatch_table2", mism_D, "paths", E(SEGPOP, "all", "paths whose tool diameter differs from the article's Table 2"))
V.set("n_paths_ae_gt_d", int((T.ae > T.D + 1e-9).sum()), "paths", E(SEGPOP, "all", "paths with radial width of cut above the tool diameter"))
prof["tool_table"] = {k: {"D": v[0], "Z": v[1]} for k, v in tooltab.items()}
nF = T[T.machine == "tmf_vf10"].groupby("run").F.nunique()
V.set("vf10_runs_multiple_f", int((nF > 1).sum()), "runs", E("tmf_vf10 runs", "label.csv", "runs with more than one feed value"))
for mach in LAB:
    sub = T[T.machine == mach]
    V.set(f"fz_min_{mach}", float(sub.fz.min()), "mm per tooth", E(SEGPOP, mach, "min feed per tooth from F, N, Z"))
    V.set(f"fz_max_{mach}", float(sub.fz.max()), "mm per tooth", E(SEGPOP, mach, "max feed per tooth from F, N, Z"))
    V.set(f"n_min_{mach}", float(sub.N.min()), "rpm", E(SEGPOP, mach, "min spindle speed"))
    V.set(f"n_max_{mach}", float(sub.N.max()), "rpm", E(SEGPOP, mach, "max spindle speed"))

# ---------------------------------------------------------------- D6 MTConnect vs label.csv ----------------------
D6_SPEED_TOL, D6_FEED_TOL, D6_SHORT, D6_LONG, GAP_THR = 0.02, 0.05, 1.0, 2.0, 1.0
V.set("d6_speed_tol", D6_SPEED_TOL, "fraction", E(SEGPOP, "D6 definition", "relative tolerance for spindle speed agreement"))
V.set("d6_feed_tol", D6_FEED_TOL, "fraction", E(SEGPOP, "D6 definition", "relative tolerance for feed agreement"))
V.set("d6_short_s", D6_SHORT, "s", E(SEGPOP, "D6 breakdown", "upper duration of a short path"))
V.set("d6_long_s", D6_LONG, "s", E(SEGPOP, "D6 breakdown", "lower duration of a long path"))
V.set("mtc_gap_thr_s", GAP_THR, "s", E("mtc.csv rows", "profile definition", "time step above which a step counts as a gap"))
for mach in LAB:
    sub = T[T.machine == mach]
    fr = (sub.mtc_Feed_Rate - sub.F) / sub.F
    V.set(f"d6_feed_relerr_median_short_{mach}", float(fr[sub.dur <= D6_SHORT].median()), "fraction", E(SEGPOP, f"{mach}, paths up to the short-path duration", "median of (recorded minus F) over F"))
    V.set(f"d6_feed_relerr_median_long_{mach}", float(fr[sub.dur > D6_LONG].median()), "fraction", E(SEGPOP, f"{mach}, paths above the long-path duration", "median of (recorded minus F) over F"))
    V.set(f"d6_n_short_{mach}", int((sub.dur <= D6_SHORT).sum()), "paths", E(SEGPOP, f"{mach}, paths up to the short-path duration", "count"))
    lg = fr[sub.dur > D6_LONG].dropna()
    V.set(f"d6_feed_within5pct_long_{mach}", float((lg.abs() <= D6_FEED_TOL).mean()), "fraction", E(SEGPOP, f"{mach}, paths above the long-path duration with an MTConnect feed", "share with recorded feed within the feed tolerance of F"))
# missing values in MTConnect other than the asset-event columns
unav_other = 0
for f in glob.glob(f"{RAW}/*/dataset*/*_mtc.csv"):
    df = pd.read_csv(f, encoding="utf-8-sig", dtype=str, keep_default_na=False)
    for c in df.columns:
        if "asset" in c:
            continue
        unav_other += int(((df[c] == "UNAVAILABLE") | (df[c] == "")).sum())
V.set("n_mtc_unavailable_other", unav_other, "cells", E("mtc.csv files", "all columns except the asset-change and asset-removal events", "cells that are UNAVAILABLE or empty"))
for mach in MACH:
    sub = T[T.machine == mach]
    sp = sub["mtc_Spindle_Speed"] if mach != "tmf_vf10" else sub["mtc_s_speed"]
    rel = (sp - sub.N) / sub.N
    ok = rel.dropna()
    V.set(f"d6_speed_within2pct_{mach}", float((ok.abs() <= D6_SPEED_TOL).mean()), "fraction", E(SEGPOP, f"{mach}, paths with an MTConnect speed", "share with recorded spindle speed within 2 percent of N"))
    V.set(f"d6_speed_relerr_median_{mach}", float(ok.median()), "fraction", E(SEGPOP, f"{mach}, paths with an MTConnect speed", "median of (recorded minus N) over N"))
    V.set(f"d6_speed_n_{mach}", int(len(ok)), "paths", E(SEGPOP, f"{mach}", "paths with an MTConnect speed"))
    if mach != "tmf_vf10":
        fr = (sub.mtc_Feed_Rate - sub.F) / sub.F
        okf = fr.dropna()
        V.set(f"d6_feed_within5pct_{mach}", float((okf.abs() <= D6_FEED_TOL).mean()), "fraction", E(SEGPOP, f"{mach}, paths with an MTConnect feed", "share with recorded feed within 5 percent of F"))
        V.set(f"d6_feed_relerr_median_{mach}", float(okf.median()), "fraction", E(SEGPOP, f"{mach}, paths with an MTConnect feed", "median of (recorded minus F) over F"))
        V.set(f"d6_feed_override_not100_{mach}", float((sub.mtc_Feed_Override.round(3) != 100).mean()), "fraction", E(SEGPOP, mach, "share of paths whose time-weighted feed override differs from 100 percent"))
        V.set(f"d6_spindle_override_not100_{mach}", float((sub.mtc_Spindle_Override.round(3) != 100).mean()), "fraction", E(SEGPOP, mach, "share of paths whose time-weighted spindle override differs from 100 percent"))
        V.set(f"ifm_single_value_share_{mach}", float((sub.mtc_smarms_nuniq <= 1).mean()), "fraction", E(SEGPOP, mach, "share of paths in which the spindle-motor vibration RMS takes one value (no update inside the path)"))
        V.set(f"watt_single_value_share_{mach}", float((sub.mtc_watt_nuniq <= 1).mean()), "fraction", E(SEGPOP, mach, "share of paths in which true power takes one value"))
# feed override in the laboratory: label.csv F against the effective path length (duration times F)
LO = T[T.machine.isin(LAB)].copy()
LO["len_eff"] = LO.dur * LO.F
run_ovr = LO.groupby(["machine", "run"]).mtc_Feed_Override.median()
ovr_runs = run_ovr[run_ovr.round(3) != 100]
V.set("n_lab_runs_override_not100", int(len(ovr_runs)), "runs", E("laboratory runs", "median time-weighted feed override over the run's paths", "runs whose median differs from 100 percent"))
ovr_detail = []
for (mach, run), ovr in ovr_runs.items():
    g = LO[(LO.machine == mach) & (LO.run == run)]
    ref = LO[(LO.machine == mach) & (LO.run != run) & (LO.tool_base.isin(g.tool_base.unique())) & (LO.Workpiece.isin(g.Workpiece.unique()))]
    ref = ref[ref.groupby(["machine", "run"]).mtc_Feed_Override.transform("median").round(3) == 100]
    key = ["F", "ae"]
    m = g.merge(ref.groupby(key).len_eff.median().rename("len_ref").reset_index(), on=key, how="inner")
    ratio = float(np.median(m.len_eff / m.len_ref)) if len(m) else float("nan")
    ovr_detail.append({"machine": mach, "run": int(run), "override_pct": float(ovr), "paths": int(len(g)), "matched_paths": int(len(m)), "len_ratio": ratio})
    V.set(f"ovr_{mach}_run{int(run)}_pct", float(ovr), "percent", E("laboratory run", f"{mach} dataset{int(run)}", "median time-weighted feed override over the run's paths"))
    V.set(f"ovr_{mach}_run{int(run)}_paths", int(len(g)), "paths", E("laboratory run", f"{mach} dataset{int(run)}", "paths"))
    V.set(f"ovr_{mach}_run{int(run)}_len_ratio", ratio, "ratio", E("laboratory run", f"{mach} dataset{int(run)}, paths matched on tool, workpiece, F and ae to runs at 100 percent override",
                                                                  "median ratio of duration times F to the matched median"))
    V.set(f"ovr_{mach}_run{int(run)}_matched", int(len(m)), "paths", E("laboratory run", f"{mach} dataset{int(run)}", "paths with a matched reference"))
prof["override_runs"] = ovr_detail

prof["d6_worst_speed"] =T.assign(rel=(T.mtc_Spindle_Speed.fillna(T.mtc_s_speed) - T.N) / T.N).sort_values("rel", key=np.abs, ascending=False)[["machine", "run", "seg", "N", "mtc_Spindle_Speed", "mtc_s_speed", "rel", "label"]].head(15).to_dict("records")

# ---------------------------------------------------------------- D9 full-scale sound ---------------------------------
D9_THR = 0.01
V.set("d9_threshold", D9_THR, "fraction", E(SEGPOP, "D9 definition", "share of samples at full scale above which a path counts as clipped"))
d9 = []
for mach in MACH:
    for ch in ("s0", "s1", "s2"):
        col = f"{ch}_fullscale"
        sub = T[(T.machine == mach) & T[col].notna()]
        if not len(sub):
            continue
        w = sub[col] * sub[f"{ch}_n"]
        V.set(f"d9_fs_share_{mach}_{ch}", float(w.sum() / sub[f"{ch}_n"].sum()), "fraction", E(SEGPOP, f"{mach}, {ch}", "share of samples at full scale during cutting"))
        V.set(f"d9_paths_gt1pct_{mach}_{ch}", float((sub[col] > D9_THR).mean()), "fraction", E(SEGPOP, f"{mach}, {ch}", "share of paths with more than 1 percent of samples at full scale"))
        pre = sub[f"{ch}_pre_fullscale"].dropna()
        V.set(f"d9_pre_paths_gt1pct_{mach}_{ch}", float((pre > D9_THR).mean()), "fraction", E("pre-cut windows", f"{mach}, {ch}, windows available", "share of windows with more than 1 percent of samples at full scale"))
        for lab in (0, 1, 2):
            s2 = sub[sub.label == lab]
            if len(s2):
                d9.append({"machine": mach, "ch": ch, "label": lab, "paths": len(s2), "share_gt1pct": float((s2[col] > D9_THR).mean()),
                           "sample_share": float((s2[col] * s2[f"{ch}_n"]).sum() / s2[f"{ch}_n"].sum())})
prof["d9"] = d9
for r in d9:
    if r["ch"] == "s0" and r["machine"] == "imi_vm20i" and r["label"] in (0, 1):
        V.set(f"d9_paths_gt1pct_imi_vm20i_s0_l{r['label']}", r["share_gt1pct"], "fraction", E(SEGPOP, f"imi_vm20i, s0, label {r['label']}", "share of paths with more than 1 percent of samples at full scale"))

# ---------------------------------------------------------------- D8 alignment ---------------------------------------
def xcorr_lag(a, b, maxlag):
    a = (a - a.mean()) / (a.std() + 1e-12); b = (b - b.mean()) / (b.std() + 1e-12)
    n = min(len(a), len(b)); a, b = a[:n], b[:n]
    best, bl = -np.inf, 0
    for L in range(-maxlag, maxlag + 1):
        if L >= 0:
            c = np.mean(a[L:] * b[:n - L])
        else:
            c = np.mean(a[:n + L] * b[-L:])
        if c > best:
            best, bl = c, L
    return bl, best


lags = []
for r in sorted(glob.glob(f"{DER}/acc_stats/*_env10hz.npy")):
    tag = os.path.basename(r).replace("_env10hz.npy", "")
    acc = np.load(r)[:, 0].astype(float)
    for ch in ("s0", "s1", "s2"):
        f = f"{DER}/audio_stats/{tag}_{ch}_env10hz.npy"
        if os.path.exists(f):
            snd = np.load(f).astype(float)
            # acc time zero is the first sample time (about 2e-5 s), negligible at 10 Hz
            L, c = xcorr_lag(np.log(acc + 1e-9), np.log(snd + 1e-9), 50)
            lags.append({"run": tag, "ch": ch, "lag_s": L / 10.0, "corr": c})
LG = pd.DataFrame(lags)
prof["d8_acc_sound_lags"] = LG.to_dict("records")
V.set("d8_n_pairs", len(LG), "pairs", E("imi_vm20i runs", "Acc0 envelope against each sound channel envelope, 10 Hz", "count"))
V.set("d8_lag_abs_max_s", float(LG.lag_s.abs().max()), "s", E("imi_vm20i runs", "Acc0 envelope against each sound channel envelope, 10 Hz, lags within 5 s", "max absolute lag at peak correlation"))
V.set("d8_share_lag0", float((LG.lag_s.abs() < 0.05).mean()), "fraction", E("imi_vm20i runs", "Acc0 envelope against each sound channel envelope, 10 Hz, lags within 5 s", "share of pairs with peak correlation at zero lag"))
V.set("d8_share_lag_le01", float((LG.lag_s.abs() <= 0.1 + 1e-9).mean()), "fraction", E("imi_vm20i runs", "Acc0 envelope against each sound channel envelope, 10 Hz, lags within 5 s", "share of pairs with peak correlation within one 0.1 s step of zero"))
V.set("d8_corr_median", float(LG["corr"].median()), "correlation", E("imi_vm20i runs", "Acc0 envelope against each sound channel envelope, 10 Hz", "median peak correlation of log envelopes"))

# MTConnect spindle load onset against sound onset at path starts (lab machines): lag of the load rise
ons = []
for (mach, run), g in T[T.machine.isin(LAB)].groupby(["machine", "run"]):
    f = glob.glob(f"{RAW}/{mach}/dataset{run}/*_mtc.csv")[0]
    df = pd.read_csv(f, encoding="utf-8-sig", usecols=["time", "S_Axis_Load"])
    env = np.load(f"{DER}/audio_stats/{mach}_dataset{run}_s2_env10hz.npy").astype(float)
    t10 = np.arange(len(env)) / 10.0
    load10 = np.interp(t10, df.time.values, df.S_Axis_Load.values.astype(float))
    L, c = xcorr_lag(load10, np.log(env + 1e-9), 50)
    ons.append({"machine": mach, "run": run, "lag_s": L / 10.0, "corr": c})
ON = pd.DataFrame(ons)
prof["d8_load_sound_lags"] = ON.to_dict("records")
D8_MAXLAG_S, D8_WIN = 5.0, (0.0, 1.0)  # search range of the cross-correlation; descriptive window for a plausible delay
V.set("d8_load_n_runs", len(ON), "runs", E("laboratory runs", "MTConnect spindle load against Sound2 envelope", "count"))
V.set("d8_maxlag_s", D8_MAXLAG_S, "s", E("D8 definition", "cross-correlation search range", "largest lag searched in either direction"))
V.set("d8_win_lo_s", D8_WIN[0], "s", E("D8 description (after seeing the lags)", "load against Sound2", "lower edge, exclusive, of the window counted as a plausible delay"))
V.set("d8_win_hi_s", D8_WIN[1], "s", E("D8 description (after seeing the lags)", "load against Sound2", "upper edge, inclusive, of the window counted as a plausible delay"))
for mach in LAB:
    o = ON[ON.machine == mach]
    w_ = o[(o.lag_s > D8_WIN[0]) & (o.lag_s <= D8_WIN[1])]
    ed = {"population": "laboratory runs", "filter": f"{mach}, MTConnect spindle load (step-held, 10 Hz) against Sound2 envelope"}
    V.set(f"d8_load_n_runs_{mach}", len(o), "runs", {**ed, "statistic": "count"})
    V.set(f"d8_load_n_inwin_{mach}", len(w_), "runs", {**ed, "statistic": "runs whose peak lag lies in the plausible-delay window"})
    V.set(f"d8_load_n_edge_{mach}", int((o.lag_s.abs() >= D8_MAXLAG_S - 1e-9).sum()), "runs", {**ed, "statistic": "runs whose peak lag sits at the edge of the search range"})
    V.set(f"d8_load_n_neg_{mach}", int((o.lag_s < 0).sum()), "runs", {**ed, "statistic": "runs whose peak lag is negative (load earlier)"})
    V.set(f"d8_load_corr_median_{mach}", float(o["corr"].median()), "correlation", {**ed, "statistic": "median peak correlation"})
    if len(w_):
        V.set(f"d8_load_lag_median_inwin_{mach}", float(w_.lag_s.median()), "s", {**ed, "statistic": "median peak lag among runs in the plausible-delay window (positive: load later)"})
        V.set(f"d8_load_corr_median_inwin_{mach}", float(w_["corr"].median()), "correlation", {**ed, "statistic": "median peak correlation among runs in the plausible-delay window"})
if len(ON[(ON.machine == "imi_vm20i") & (ON.lag_s > D8_WIN[0]) & (ON.lag_s <= D8_WIN[1])]):
    V.set("d8_lag_over_median_path_imi_vm20i", float(V.get("d8_load_lag_median_inwin_imi_vm20i") / V.get("dur_median_imi_vm20i")), "fraction",
          E(SEGPOP, "imi_vm20i", "median in-window load lag over the median path duration"))

# ---------------------------------------------------------------- encoding -----------------------------------------
bom = crlf = 0; nfiles = 0
for f in glob.glob(f"{RAW}/*/dataset*/*.csv"):
    raw = open(f, "rb").read(4096); nfiles += 1
    bom += raw.startswith(b"\xef\xbb\xbf"); crlf += b"\r\n" in raw
V.set("n_csv_small", nfiles, "files", E("CSV files other than acc.csv", "all", "count"))
V.set("n_csv_bom", bom, "files", E("CSV files other than acc.csv", "all", "files starting with a UTF-8 byte order mark"))
V.set("n_csv_crlf", crlf, "files", E("CSV files other than acc.csv", "all", "files with CRLF line endings"))
mtc_bom = sum(open(f, "rb").read(3) == b"\xef\xbb\xbf" for f in glob.glob(f"{RAW}/*/dataset*/*_mtc.csv"))
V.set("n_mtc_bom", int(mtc_bom), "files", E("mtc.csv files", "all", "files starting with a UTF-8 byte order mark"))

# ---------------------------------------------------------------- accelerometer extremes ---------------------------
A = pd.DataFrame([json.load(open(f)) for f in glob.glob(f"{DER}/acc_stats/*.json")])
a0 = pd.DataFrame(list(A.acc0)); a1 = pd.DataFrame(list(A.acc1))
V.set("acc0_abs_max_g", float(max(a0["max"].max(), -a0["min"].min()) / 9806.65), "g", E("acc.csv files", "Acc0, all samples", "largest absolute value, mm/s2 divided by 9806.65"))
V.set("acc1_abs_max_g", float(max(a1["max"].max(), -a1["min"].min()) / 9806.65), "g", E("acc.csv files", "Acc1, all samples", "largest absolute value, mm/s2 divided by 9806.65"))
V.set("acc0_mean_max_abs", float(a0["mean"].abs().max()), "mm/s2", E("acc.csv files", "Acc0", "largest absolute file mean (DC offset)"))
prof["acc_extreme_runs"] = A.assign(a0max=a0["max"], a0min=a0["min"], a1max=a1["max"]).sort_values("a0max", ascending=False)[["run", "a0max", "a0min", "a1max"]].head(5).to_dict("records")

json.dump(prof, open(f"{OUT}/profile.json", "w"), indent=1, default=float)
V.save()
print("p02_profile: registered", len(V._touched), "values")