"""Sensor correlation network for SECOM, and whether its hubs make good features. **Honest label:** this is a *data-driven correlation network*, not a process flow graph. SECOM's sensors are anonymised, so nothing here says which process step a sensor belongs to. The genuine process-flow story lives in 03_bosch. What this network can do is find sensors that sit at the centre of the measurement system. Method transfer: correlate every pair of channels, keep the strongest fixed fraction of edges (proportional thresholding), then rank nodes by degree centrality and call the top ones hubs. That is the standard fMRI functional connectivity pipeline applied to fab sensors instead of brain regions. The question this script actually answers is a sceptical one: **does selecting features by network centrality beat selecting them by plain mutual information with the label?** Reported as a curve over k, not a single cherry-picked number, and against a random-k control. Run: python 01_secom/sensor_network.py """ from __future__ import annotations import sys from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import networkx as nx import numpy as np import pandas as pd from sklearn.feature_selection import mutual_info_classif from xgboost import XGBClassifier HERE = Path(__file__).parent sys.path.insert(0, str(HERE.parent)) sys.path.insert(0, str(HERE)) from common import metrics, splits # noqa: E402 from run_pipeline import (COST_FN, COST_FP, load, preprocess_apply, # noqa: E402 preprocess_fit) FIGS = HERE / "figures" OUT = HERE / "outputs" EDGE_DENSITY = 0.05 # keep the strongest 5% of sensor pairs K_VALUES = [10, 20, 30, 50, 80, 120] N_SEEDS = 10 # random-control repeats; their spread IS the noise floor def build_network(X_train: np.ndarray, names: list[str], density: float = EDGE_DENSITY) -> tuple[nx.Graph, pd.DataFrame]: """Proportional-threshold correlation graph + centrality table. Correlations are computed on training rows only — the network is part of the model, so it must not see the test period. """ corr = np.corrcoef(X_train, rowvar=False) corr = np.nan_to_num(corr, nan=0.0) np.fill_diagonal(corr, 0.0) absr = np.abs(corr) iu = np.triu_indices_from(absr, k=1) weights = absr[iu] cutoff = np.quantile(weights, 1 - density) keep = weights >= cutoff G = nx.Graph() G.add_nodes_from(names) for a, b, w in zip(iu[0][keep], iu[1][keep], weights[keep]): G.add_edge(names[a], names[b], weight=float(w)) deg = dict(G.degree()) strength = {n: float(sum(d["weight"] for d in G[n].values())) for n in G} btw = nx.betweenness_centrality(G, k=min(150, G.number_of_nodes()), seed=0) table = pd.DataFrame({ "sensor": names, "degree": [deg.get(n, 0) for n in names], "strength": [strength.get(n, 0.0) for n in names], "betweenness": [btw.get(n, 0.0) for n in names], }) table["degree_centrality"] = table["degree"] / max(1, len(names) - 1) return G, table.sort_values("degree", ascending=False).reset_index(drop=True) def evaluate_subset(cols: list[int], X_fit, y_fit, X_val, y_val, X_te, y_te) -> dict: """Same model, same split, only the feature subset changes.""" spw = (y_fit == 0).sum() / max(1, (y_fit == 1).sum()) mdl = XGBClassifier(n_estimators=300, max_depth=3, learning_rate=0.05, subsample=0.8, colsample_bytree=0.6, scale_pos_weight=spw, eval_metric="aucpr", reg_lambda=2.0, n_jobs=-1, random_state=42) mdl.fit(X_fit[:, cols], y_fit) val_s = mdl.predict_proba(X_val[:, cols])[:, 1] thr, _ = metrics.best_threshold(y_val, val_s, COST_FN, COST_FP) mdl.fit(np.vstack([X_fit, X_val])[:, cols], np.concatenate([y_fit, y_val])) te_s = mdl.predict_proba(X_te[:, cols])[:, 1] return { "pr_auc": metrics.pr_auc(y_te, te_s), "capture@20%": metrics.decile_capture(y_te, te_s, 0.2), "expected_cost": metrics.expected_cost( y_te, (te_s >= thr).astype(int), COST_FN, COST_FP), } def draw_network(G: nx.Graph, table: pd.DataFrame, n_hubs: int = 12) -> None: """Plot the largest connected component; size and colour by degree.""" comp = max(nx.connected_components(G), key=len) H = G.subgraph(comp) deg = dict(H.degree()) hubs = set(table.head(n_hubs)["sensor"]) & set(H.nodes) pos = nx.spring_layout(H, seed=1, k=0.32, iterations=120) fig, ax = plt.subplots(figsize=(12, 9)) nx.draw_networkx_edges(H, pos, ax=ax, alpha=0.10, edge_color="gray", width=0.5) d = np.array([deg[n] for n in H.nodes], dtype=float) nx.draw_networkx_nodes(H, pos, ax=ax, nodelist=list(H.nodes), node_size=12 + 340 * (d / d.max()) ** 1.4, node_color=d, cmap="viridis", edgecolors="black", linewidths=0.3) nx.draw_networkx_nodes(H, pos, ax=ax, nodelist=list(hubs), node_size=[12 + 340 * (deg[n] / d.max()) ** 1.4 + 150 for n in hubs], node_color="none", edgecolors="crimson", linewidths=1.8) nx.draw_networkx_labels(H, pos, ax=ax, font_size=7.5, labels={n: n for n in hubs}) ax.set_title(f"SECOM sensor correlation network (training period only)\n" f"{H.number_of_nodes()} sensors, {H.number_of_edges():,} edges " f"at {EDGE_DENSITY:.0%} density — red rings = top-{n_hubs} degree hubs\n" f"data-driven correlation, NOT a process flow graph", fontsize=11) ax.axis("off") fig.tight_layout() fig.savefig(FIGS / "sensor_correlation_network.png", dpi=150) plt.close(fig) def main(): FIGS.mkdir(exist_ok=True) OUT.mkdir(exist_ok=True) df = load() sensor_cols = [c for c in df.columns if c.startswith("s")] tr_idx, te_idx = splits.temporal_split(df, "timestamp", train_frac=0.7) train, test = df.loc[tr_idx], df.loc[te_idx] v_idx, _ = splits.temporal_split(train, "timestamp", train_frac=0.7) fit, val = train.loc[v_idx], train.drop(v_idx) prep = preprocess_fit(fit, sensor_cols) names = prep["cols"] X_fit, X_val, X_te = (preprocess_apply(d, prep) for d in (fit, val, test)) y_fit, y_val, y_te = (d["y"].to_numpy() for d in (fit, val, test)) print(f"{len(names)} sensors after preprocessing; " f"fit {len(y_fit)} / val {len(y_val)} / test {len(y_te)}") G, table = build_network(X_fit, names) print(f"network: {G.number_of_nodes()} nodes, {G.number_of_edges():,} edges " f"at {EDGE_DENSITY:.0%} density") print("\n=== top network hubs by degree ===") print(table.head(10).to_string(index=False)) table.to_csv(OUT / "sensor_network_centrality.csv", index=False) draw_network(G, table) # --- selection methods, all fit on the training window only mi = mutual_info_classif(X_fit, y_fit, random_state=0) mi_order = np.argsort(mi)[::-1] hub_order = [names.index(s) for s in table["sensor"]] rng = np.random.default_rng(7) rows = [] for k in K_VALUES: for label, cols in [("network hubs (degree)", hub_order[:k]), ("mutual information top-k", list(mi_order[:k]))]: r = evaluate_subset(cols, X_fit, y_fit, X_val, y_val, X_te, y_te) rows.append({"method": label, "k": k, **r}) rand = [evaluate_subset(list(rng.choice(len(names), k, replace=False)), X_fit, y_fit, X_val, y_val, X_te, y_te) for _ in range(N_SEEDS)] row = {"method": f"random control (mean of {N_SEEDS})", "k": k, **{m: float(np.mean([x[m] for x in rand])) for m in rand[0]}} # the spread of equally-sized random subsets is the noise floor any # claimed improvement has to clear row.update({f"{m}_sd": float(np.std([x[m] for x in rand])) for m in rand[0]}) rows.append(row) print(f" k={k}: random PR-AUC {row['pr_auc']:.4f} +/- {row['pr_auc_sd']:.4f}") allf = evaluate_subset(list(range(len(names))), X_fit, y_fit, X_val, y_val, X_te, y_te) res = pd.DataFrame(rows) res.to_csv(OUT / "feature_selection_comparison.csv", index=False) print("\n=== feature selection comparison (test period) ===") print(res.pivot(index="k", columns="method", values="pr_auc").round(4).to_string()) print(f"\nall {len(names)} features: PR-AUC {allf['pr_auc']:.4f}, " f"capture@20% {allf['capture@20%']:.4f}") fig, axes = plt.subplots(1, 2, figsize=(13, 4.8)) colors = {"network hubs (degree)": "crimson", "mutual information top-k": "steelblue", f"random control (mean of {N_SEEDS})": "gray"} for ax, metric in zip(axes, ["pr_auc", "capture@20%"]): for method, g in res.groupby("method"): ax.plot(g["k"], g[metric], marker="o", ms=4, label=method, color=colors.get(method), ls="--" if "random" in method else "-") if "random" in method: ax.fill_between(g["k"], g[metric] - g[f"{metric}_sd"], g[metric] + g[f"{metric}_sd"], color="gray", alpha=0.22, label="random ±1 sd (noise floor)") ax.axhline(allf[metric], color="black", ls=":", lw=1, label=f"all {len(names)} features") if metric == "pr_auc": ax.axhline(y_te.mean(), color="darkgreen", ls=":", lw=1, label="random ranker (prevalence)") ax.set(xlabel="k features selected", ylabel=metric, title=f"{metric} vs number of features (test period)") ax.legend(fontsize=7.5) fig.tight_layout() fig.savefig(FIGS / "feature_selection_comparison.png", dpi=150) plt.close(fig) if __name__ == "__main__": main()