"""seg_table.py: build one row per labelled cutting path (segment) with design variables, high-rate channel features (from the device-side d03 pass) and MTConnect / CT summaries. Used by p02_profile.py and p03_stats.py. Writes outputs/segments.parquet and outputs/spectra.npz. No statistic is computed here; this is the shared data layer. """ import json, glob, os, re import numpy as np, pandas as pd RAW = "/home/claude/mill/raw" SEG = f"{RAW}/Phase2/derived/seg_features" OUT = "/home/claude/mill/outputs" MACH = ("imi_vm20i", "imi_vmx30ui", "tmf_vf10") CH = ("acc0", "acc1", "s0", "s1", "s2") FEAT = ["rms", "peak", "crest", "kurt", "fullscale", "n"] BANDS = [(50, 500), (500, 2000), (2000, 5000), (5000, 12000), (50, 12000)] MIN_SPINDLE_HZ_FOR_COMB = 50.0 # below this the harmonic comb leaves no clean gaps at 6 Hz bins COMB_HALF_BINS = 2 # bins excluded on each side of every spindle harmonic NH_BAND = (200.0, 10000.0) # band searched for the dominant line (Hz) MTC_LAB = ["Spindle_Speed", "Feed_Rate", "S_Axis_Load", "X_Axis_Load", "Y_Axis_Load", "Z_Axis_Load", "watt", "smarms", "smapeak", "smvrms", "smcrest", "Feed_Override", "Spindle_Override"] MTC_VF = ["s_speed", "s_load"] def spectral_features(P, f, fsp): """Band powers, spectral centroid and the largest peak away from spindle harmonics.""" out = {} for lo, hi in BANDS: m = (f >= lo) & (f < hi) out[f"bp_{lo}_{hi}"] = float(np.trapezoid(P[m], f[m])) if m.any() else np.nan m = (f >= 50) & (f < 12000) out["centroid"] = float((f[m] * P[m]).sum() / P[m].sum()) if P[m].sum() > 0 else np.nan out["nh_peak_f"] = np.nan; out["nh_peak_ratio"] = np.nan; out["nh_share"] = np.nan; out["h_peak_ratio"] = np.nan df = f[1] - f[0] if fsp >= MIN_SPINDLE_HZ_FOR_COMB and np.isfinite(P[m]).all(): band = (f >= NH_BAND[0]) & (f < NH_BAND[1]) harm = np.zeros_like(f, dtype=bool) k = np.arange(1, int(10000 / fsp) + 2) idx = np.round(k * fsp / df).astype(int) for j in range(-COMB_HALF_BINS, COMB_HALF_BINS + 1): ii = idx + j ii = ii[(ii >= 0) & (ii < len(f))] harm[ii] = True nh = band & ~harm hh = band & harm med = np.median(P[band]) if nh.any() and med > 0: i = np.argmax(np.where(nh, P, -np.inf)) out["nh_peak_f"] = float(f[i]); out["nh_peak_ratio"] = float(P[i] / med) out["nh_share"] = float(P[nh].sum() / P[band].sum()) ih = np.argmax(np.where(hh, P, -np.inf)) out["h_peak_ratio"] = float(P[ih] / med) return out def mtc_summaries(mach, run, starts, ends): f = glob.glob(f"{RAW}/{mach}/dataset{run}/*_mtc.csv")[0] df = pd.read_csv(f, encoding="utf-8-sig", low_memory=False) cols = MTC_LAB if mach != "tmf_vf10" else MTC_VF t = df.time.values res = [] for s, e in zip(starts, ends): # MTConnect rows hold the latest observation, so a value applies from its time stamp until the next row: # use rows strictly inside the path plus the last row before it (step-function average). i0 = max(np.searchsorted(t, s, side="right") - 1, 0) i1 = np.searchsorted(t, e, side="right") sub = df.iloc[i0:i1] tt = np.clip(np.r_[sub.time.values, e], s, e) w = np.diff(tt) r = {"mtc_rows_in_path": int(((t >= s) & (t <= e)).sum())} for c in cols: v = pd.to_numeric(sub[c], errors="coerce").values ok = np.isfinite(v) & (w > 0) r[f"mtc_{c}"] = float((v[ok] * w[ok]).sum() / w[ok].sum()) if ok.any() else np.nan if c in ("smarms", "watt", "s_load", "S_Axis_Load"): # number of distinct values inside the path: the sensors update about once per second r[f"mtc_{c}_nuniq"] = int(pd.Series(v[ok]).nunique()) # pre-cut level for power and spindle load: same window as the high-rate pre-cut window p0, p1 = s - 1.2, s - 0.2 if p0 >= 0: j0 = max(np.searchsorted(t, p0, side="right") - 1, 0); j1 = np.searchsorted(t, p1, side="right") pre = df.iloc[j0:j1] for c in (("watt", "S_Axis_Load", "Spindle_Speed") if mach != "tmf_vf10" else ("s_load", "s_speed")): v = pd.to_numeric(pre[c], errors="coerce").values r[f"mtc_pre_{c}"] = float(np.nanmean(v)) if len(v) else np.nan res.append(r) return pd.DataFrame(res) def ct_summaries(run, starts, ends): f = glob.glob(f"{RAW}/tmf_vf10/dataset{run}/*_ct.csv")[0] ct = pd.read_csv(f) out = [] for s, e in zip(starts, ends): m = (ct.time >= s) & (ct.time <= e) out.append({"ct_mean": float(ct.value[m].mean()) if m.any() else np.nan, "ct_n": int(m.sum())}) return pd.DataFrame(out) def build(save=True): rows = []; spectra = {} for mach in MACH: for d in sorted(glob.glob(f"{RAW}/{mach}/dataset*"), key=lambda p: int(p.split("dataset")[-1])): run = int(d.split("dataset")[-1]) lab = pd.read_csv(f"{d}/label.csv") z = np.load(f"{SEG}/{mach}_dataset{run}.npz") meta = json.loads(str(z["meta"])) assert np.allclose(z["start"], lab.start.values) and np.allclose(z["end"], lab.end.values) base = lab.copy() base.insert(0, "run", run); base.insert(0, "machine", mach) base["seg"] = np.arange(len(lab)) base["dur"] = base.end - base.start base["tool_base"] = base.Tool.astype(str).str.replace("*", "", regex=False) base["defective"] = base.Tool.astype(str).str.contains(r"\*") base["fz"] = base.F * 60 / (base.N * base.Z) base["vc"] = np.pi * base.D * base.N / 1000 base["mrr"] = base.ap * base.ae * base.F base["f_spindle"] = base.N / 60 base["tpf"] = base.N * base.Z / 60 for ch in CH: if f"{ch}_feat" not in z: continue F = z[f"{ch}_feat"]; FP = z[f"{ch}_pre_feat"] for j, c in enumerate(FEAT): base[f"{ch}_{c}"] = F[:, j] base[f"{ch}_pre_rms"] = FP[:, 0] base[f"{ch}_pre_fullscale"] = FP[:, 4] P = z[f"{ch}_psd"].astype(np.float64) fs = meta["channels"][ch]["fs"]; f = np.arange(P.shape[1]) * fs / meta["nfft"] sf = [spectral_features(P[i], f, base.f_spindle.iloc[i]) for i in range(len(base))] sf = pd.DataFrame(sf).add_prefix(f"{ch}_") base = pd.concat([base.reset_index(drop=True), sf], axis=1) spectra[f"{mach}|{run}|{ch}"] = z[f"{ch}_psd"] spectra[f"{mach}|{run}|{ch}|pre"] = z[f"{ch}_pre_psd"] m = mtc_summaries(mach, run, lab.start.values, lab.end.values) base = pd.concat([base.reset_index(drop=True), m], axis=1) if mach == "tmf_vf10": base = pd.concat([base.reset_index(drop=True), ct_summaries(run, lab.start.values, lab.end.values)], axis=1) rows.append(base) T = pd.concat(rows, ignore_index=True) T["run_id"] = T.machine + ":" + T.run.astype(str) if save: os.makedirs(OUT, exist_ok=True) T.to_parquet(f"{OUT}/segments.parquet", index=False) return T if __name__ == "__main__": T = build() print(T.shape) print(T.groupby("machine").size())