#!/usr/bin/env python3 """Best possible fine-grid IoU when each fine node inherits one coarse-node prediction. This is a representation ceiling for the current nearest-node evaluation map, not a score achieved by a trained model or a bound on other decoders. """ from __future__ import annotations import csv from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np from scipy.spatial import cKDTree HERE = Path(__file__).resolve().parent ROOT = HERE.parents[1] RUNS = (("AC28_16h", ROOT, "16h"), ("AC28_18h01", ROOT, "18h01"), ("ME1a_14h06", ROOT / "experiments/ME1a_DT", "14h06")) BUDGETS = (3000, 5000, 6000, 12000) def best_projected_binary(target: np.ndarray, group: np.ndarray, groups: int) -> dict: """Maximize TP/(P+FP) over all binary coarse-node assignments. At candidate IoU lambda, independently turn on groups for which TP_group - lambda*FP_group > 0; Dinkelbach iteration finds the optimum. """ p = int(target.sum()) tp_by_group = np.bincount(group, weights=target.astype(np.int64), minlength=groups).astype(np.int64) total_by_group = np.bincount(group, minlength=groups) fp_by_group = total_by_group - tp_by_group positive_groups = np.flatnonzero(tp_by_group > 0) # AP is the weighted precision at each newly retrieved positive group. # All-negative groups go last. Dynamic programming finds the best ranking # of the positive-bearing groups, respecting tied fine nodes per group. ranked_tp = tp_by_group[positive_groups] ranked_total = total_by_group[positive_groups] if np.array_equal(ranked_tp, ranked_total): oracle_ap = 1.0 else: count = len(positive_groups) if count > 16: raise ValueError("Exact AP ranking needs a smaller positive-group count") subset_tp = np.zeros(1 << count, dtype=np.int64) subset_total = np.zeros(1 << count, dtype=np.int64) best_ap = np.full(1 << count, -np.inf) best_ap[0] = 0.0 for subset in range(1, 1 << count): bit = subset & -subset index = bit.bit_length() - 1 previous = subset ^ bit subset_tp[subset] = subset_tp[previous] + ranked_tp[index] subset_total[subset] = subset_total[previous] + ranked_total[index] for index in range(count): if subset & (1 << index): previous = subset ^ (1 << index) gain = (ranked_tp[index] / p) * (subset_tp[subset] / subset_total[subset]) best_ap[subset] = max(best_ap[subset], best_ap[previous] + gain) oracle_ap = float(best_ap[-1]) score = 0.0 for _ in range(100): chosen = tp_by_group > score * fp_by_group tp = int(tp_by_group[chosen].sum()) fp = int(fp_by_group[chosen].sum()) updated = tp / (p + fp) if abs(updated - score) < 1e-14: score = updated break score = updated else: raise RuntimeError("IoU optimization did not converge") fn = p - tp return {"weak_positive_fine": p, "positive_bearing_groups": int((tp_by_group > 0).sum()), "oracle_tp": tp, "oracle_fp": fp, "oracle_fn": fn, "oracle_iou": score, "oracle_dice": 2 * tp / max(1, 2 * tp + fp + fn), "oracle_ap": oracle_ap} def main() -> None: with (HERE / "gnn_5000_vs_12000.csv").open(newline="") as f: observed = {(int(r["requested_nodes"]), r["run"]): r for r in csv.DictReader(f) if r["fit"] == "individual"} rows = [] for run, root, run_id in RUNS: with np.load(root / f"data_processed/paired/{run_id}_cad_ply_targets.npz", allow_pickle=False) as d: fine_pos = np.asarray(d["cad_pos_mm"], np.float32) target = np.asarray(d["weak_defect_mask"][:, 1], bool) for budget in BUDGETS: with np.load(root / f"data_processed/graphs/graph_{budget}.npz", allow_pickle=False) as d: coarse_pos = np.asarray(d["pos_mm"], np.float32) group = cKDTree(coarse_pos).query(fine_pos, workers=-1)[1] bound = best_projected_binary(target, group, len(coarse_pos)) training = observed.get((budget, run)) rows.append({"run": run, "requested_nodes": budget, "actual_nodes": len(coarse_pos), **bound, "observed_hard_individual_ap": training["same_or_held_out_weak_ap"] if training else "", "observed_hard_individual_iou": training["same_or_held_out_weak_iou"] if training else "", "interpretation": "IoU ceiling for current nearest-coarse-to-fine binary map only"}) output = HERE / "projection_oracle.csv" with output.open("w", newline="") as f: writer = csv.DictWriter(f, fieldnames=list(rows[0])) writer.writeheader(); writer.writerows(rows) fig, axs = plt.subplots(2, 3, figsize=(14, 7.5), sharey="row", sharex="col", constrained_layout=True) for column, (run, _, _) in enumerate(RUNS): subset = [r for r in rows if r["run"] == run] budgets = [r["requested_nodes"] for r in subset] trained = [r for r in subset if r["observed_hard_individual_iou"] != ""] for row, metric in enumerate(("iou", "ap")): ax = axs[row, column] ax.plot(budgets, [r[f"oracle_{metric}"] for r in subset], "o-", color="#087f8c", linewidth=2.2, label="Best possible with nearest-node transfer") ax.scatter([r["requested_nodes"] for r in trained], [float(r[f"observed_hard_individual_{metric}"]) for r in trained], color="#d76827", s=75, marker="D", zorder=5, label="Observed 180-epoch hard GNN") ax.set_xticks(budgets) ax.set_xlabel("Requested graph nodes") ax.set_ylim(0, 1.06) ax.grid(alpha=.2) axs[0, column].set_title(run.replace("_", " ")) if run == "ME1a_14h06": axs[0, column].annotate("5k max 0.458 < 12k fit 0.476", xy=(5000, subset[1]["oracle_iou"]), xytext=(3100, .68), arrowprops={"arrowstyle": "->", "color": "#37485a"}, fontsize=9) axs[1, column].annotate("5k max 0.626 < 12k fit 0.629", xy=(5000, subset[1]["oracle_ap"]), xytext=(3100, .82), arrowprops={"arrowstyle": "->", "color": "#37485a"}, fontsize=9) axs[0, 0].set_ylabel("Fine-grid IoU") axs[1, 0].set_ylabel("Fine-grid average precision") axs[1, -1].legend(loc="lower right", fontsize=8) fig.suptitle("Projection ceilings under the existing nearest-node evaluation map (not trained scores)", fontsize=13, weight="bold") fig.savefig(HERE / "projection_oracle.png", dpi=165, facecolor="white") plt.close(fig) print(output) for r in rows: if r["run"] == "ME1a_14h06": print(r["requested_nodes"], r["actual_nodes"], f"oracle IoU {r['oracle_iou']:.4f}", f"oracle AP {r['oracle_ap']:.4f}", f"TP/FP/FN {r['oracle_tp']}/{r['oracle_fp']}/{r['oracle_fn']}") if __name__ == "__main__": main()