diff --git a/configs/reproduce_ellphi_main.yaml b/configs/reproduce_ellphi_main.yaml new file mode 100644 index 0000000..b33dbe3 --- /dev/null +++ b/configs/reproduce_ellphi_main.yaml @@ -0,0 +1,48 @@ +# Post-fix ellphi multiseed main run (same hyperparams as reproduce.yaml; lighter viz I/O). +# Pilot calibration uses reproduce.yaml; this config is for the 5-seed main phase. + +meta: + config_id: "reproduce_ellphi_main" + +model: + point_dim: 2 + feature_dim: 128 + ellipse_param_dim: 5 + threshold: 0.5 + topology_loss: + metric_type: "anisotropic" + distance_backend: "mahalanobis" + ellphi_differentiable: true + +loss: + w_class: 1.0 + w_topo: 0.1 + w_aniso: 0.01099204345474479 + w_size: 0.0055785823086202556 + pos_weight: 1.0 + aniso_mode: "linear" + +training: + lr: 0.0004897466143769238 + epochs: 50 + grad_clip_value: 1.0 + visualize_every: 10 + warmup_epochs: 0 + use_amp: true + +data: + max_points: 200 + num_outliers: 20 + seed: 42 + test_size: 1000 + num_workers: 4 + batch_size: 64 + pin_memory: true + +performance: + cudnn_benchmark: true + enable_tf32: true + matmul_precision: high + +outputs: + base_dir: "outputs" diff --git a/experiments/aggregate_paper_results.py b/experiments/aggregate_paper_results.py new file mode 100644 index 0000000..00a2cca --- /dev/null +++ b/experiments/aggregate_paper_results.py @@ -0,0 +1,441 @@ +#!/usr/bin/env python3 +""" +Aggregate paper Table~\\ref{tab:comparison} metrics into ``summary_for_paper.csv``. + +Reads test-split DBSCAN metrics (``logs/paper_metrics_test.json`` from +``evaluate_paper_protocol.py``) for proposed / ablation runs, and +``summary_baselines.csv`` from Phase 3. + +Usage:: + + uv run python experiments/aggregate_paper_results.py \\ + --proposed-dir outputs/paper_reproduce_1week_tuned \\ + --wtopo0-dir outputs/paper_wtopo0 \\ + --baselines-csv outputs/paper_baselines/summary_baselines.csv + +Writes ``summary_for_paper.csv`` and ``MANIFEST.md`` under ``--out-dir`` +(default: ``--proposed-dir``). +""" + +from __future__ import annotations + +import argparse +import csv +import json +from dataclasses import dataclass +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Sequence + +import numpy as np + +from tda_ml.supervised_diagnostics import git_revision + +REPO_ROOT = Path(__file__).resolve().parents[1] +PAPER_SEEDS = [42, 123, 456, 789, 1024] + +SUMMARY_COLUMNS = [ + "method", + "mcc_mean", + "mcc_std", + "gmean_mean", + "gmean_std", + "wdist_mean", + "wdist_std", + "notes", +] + +METHOD_ORDER = [ + "proposed_w_topo_pos", + "proposed_w_topo_0", + "euclidean_dbscan", + "isolation_forest", + "lof", + "adbscan", +] + +# ``evaluate_paper_protocol`` writes SplitMetrics fields (mcc, gmean, wdist). +# Mac-side stubs may use alternate names — try all aliases. +MCC_KEYS = ("mcc", "test_mcc", "mcc_mean", "mean_mcc") +GMEAN_KEYS = ("gmean", "g_mean", "test_gmean", "gmean_mean", "mean_gmean") +WDIST_KEYS = ("wdist", "w_dist", "test_wdist", "wdist_mean", "mean_wdist", "w_dist_mean") + + +@dataclass +class SeedPaperMetrics: + seed: int | None + run_dir: Path + metrics_path: Path + mcc: float + gmean: float + wdist: float + n_clouds: int | None + dbscan_eps: float | None + dbscan_min_samples: int | None + backend: str | None + source_revision: str | None + + +def sample_std(values: Sequence[float]) -> float: + arr = np.asarray(values, dtype=np.float64) + if arr.size < 2: + return 0.0 + return float(np.std(arr, ddof=1)) + + +def _first_key(payload: dict[str, Any], keys: Sequence[str], *, label: str) -> float: + for key in keys: + if key in payload and payload[key] is not None: + return float(payload[key]) + raise KeyError(f"Missing {label} in JSON (tried {list(keys)}); keys={sorted(payload)}") + + +def load_paper_metrics_json(path: Path) -> dict[str, Any]: + payload = json.loads(path.read_text(encoding="utf-8")) + if payload.get("split") not in (None, "test"): + raise ValueError(f"Expected test-split metrics in {path}; got split={payload.get('split')!r}") + return payload + + +def parse_seed_paper_metrics(run_dir: Path, metrics_path: Path) -> SeedPaperMetrics: + payload = load_paper_metrics_json(metrics_path) + manifest_path = run_dir / "logs" / "run_manifest.json" + seed: int | None = None + if manifest_path.is_file(): + manifest = json.loads(manifest_path.read_text(encoding="utf-8")) + if manifest.get("seed") is not None: + seed = int(manifest["seed"]) + + return SeedPaperMetrics( + seed=seed, + run_dir=run_dir.resolve(), + metrics_path=metrics_path.resolve(), + mcc=_first_key(payload, MCC_KEYS, label="mcc"), + gmean=_first_key(payload, GMEAN_KEYS, label="gmean"), + wdist=_first_key(payload, WDIST_KEYS, label="wdist"), + n_clouds=int(payload["n_clouds"]) if payload.get("n_clouds") is not None else None, + dbscan_eps=float(payload["dbscan_eps"]) if payload.get("dbscan_eps") is not None else None, + dbscan_min_samples=( + int(payload["dbscan_min_samples"]) if payload.get("dbscan_min_samples") is not None else None + ), + backend=str(payload["backend"]) if payload.get("backend") is not None else None, + source_revision=str(payload["source_revision"]) if payload.get("source_revision") else None, + ) + + +def discover_run_dirs(out_base: Path) -> list[Path]: + """Return run directories with ``logs/paper_metrics_test.json``.""" + found: list[Path] = [] + progress = out_base / "progress_summary.csv" + if progress.is_file(): + with progress.open(newline="", encoding="utf-8") as f: + for row in csv.DictReader(f): + if row.get("backend") not in (None, "", "ellphi"): + continue + run_dir = Path(row["run_dir"]) + if (run_dir / "logs" / "paper_metrics_test.json").is_file(): + found.append(run_dir.resolve()) + + if not found: + for metrics_path in sorted(out_base.glob("backend_ellphi_seed*/logs/paper_metrics_test.json")): + found.append(metrics_path.parent.parent.resolve()) + for metrics_path in sorted(out_base.glob("reproduce_*/logs/paper_metrics_test.json")): + found.append(metrics_path.parent.parent.resolve()) + + # Deduplicate while preserving order + seen: set[Path] = set() + unique: list[Path] = [] + for rd in found: + if rd not in seen: + seen.add(rd) + unique.append(rd) + return unique + + +def aggregate_seed_metrics( + seeds: list[SeedPaperMetrics], + *, + method: str, + notes: str, +) -> dict[str, Any]: + if not seeds: + raise ValueError(f"No seed metrics for method={method!r}") + mccs = [s.mcc for s in seeds] + gmeans = [s.gmean for s in seeds] + wdist = [s.wdist for s in seeds] + n_clouds = seeds[0].n_clouds + cloud_note = f"test n_clouds={n_clouds} per seed" if n_clouds is not None else "test split" + return { + "method": method, + "mcc_mean": float(np.mean(mccs)), + "mcc_std": sample_std(mccs), + "gmean_mean": float(np.mean(gmeans)), + "gmean_std": sample_std(gmeans), + "wdist_mean": float(np.mean(wdist)), + "wdist_std": sample_std(wdist), + "notes": notes or f"{len(seeds)} seeds; DBSCAN test metrics; {cloud_note}", + } + + +def load_proposed_row( + out_base: Path, + *, + method: str, + default_notes: str, + expected_seeds: Sequence[int], +) -> tuple[dict[str, Any] | None, list[SeedPaperMetrics], list[str]]: + warnings: list[str] = [] + run_dirs = discover_run_dirs(out_base) + if not run_dirs: + warnings.append(f"{method}: no run dirs with paper_metrics_test.json under {out_base}") + return None, [], warnings + + per_seed: list[SeedPaperMetrics] = [] + for run_dir in run_dirs: + metrics_path = run_dir / "logs" / "paper_metrics_test.json" + try: + per_seed.append(parse_seed_paper_metrics(run_dir, metrics_path)) + except (KeyError, ValueError, json.JSONDecodeError) as exc: + warnings.append(f"{method}: skip {run_dir}: {exc}") + + if not per_seed: + warnings.append(f"{method}: no parseable paper_metrics_test.json under {out_base}") + return None, [], warnings + + found_seeds = {s.seed for s in per_seed if s.seed is not None} + missing = [s for s in expected_seeds if s not in found_seeds] + if missing: + warnings.append(f"{method}: missing seeds {missing} (found {sorted(found_seeds)})") + + row = aggregate_seed_metrics(per_seed, method=method, notes=default_notes) + return row, per_seed, warnings + + +def load_baselines_rows(path: Path) -> list[dict[str, Any]]: + if not path.is_file(): + raise FileNotFoundError(f"Baselines CSV not found: {path}") + with path.open(newline="", encoding="utf-8") as f: + rows = list(csv.DictReader(f)) + for row in rows: + for col in SUMMARY_COLUMNS: + if col not in row: + row[col] = "" + return rows + + +def write_summary_csv(path: Path, rows: Sequence[dict[str, Any]]) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + with path.open("w", newline="", encoding="utf-8") as f: + writer = csv.DictWriter(f, fieldnames=SUMMARY_COLUMNS) + writer.writeheader() + for row in rows: + writer.writerow({col: row.get(col, "") for col in SUMMARY_COLUMNS}) + + +def format_pm(mean: float, std: float, digits: int = 3) -> str: + return f"{mean:.{digits}f} ± {std:.{digits}f}" + + +def write_manifest( + path: Path, + *, + summary_csv: Path, + proposed_dir: Path, + wtopo0_dir: Path | None, + baselines_csv: Path, + rows: Sequence[dict[str, Any]], + proposed_seeds: list[SeedPaperMetrics], + wtopo0_seeds: list[SeedPaperMetrics], + warnings: Sequence[str], +) -> None: + now = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M:%S UTC") + lines = [ + "# Paper experiment manifest", + "", + f"- Generated: {now}", + f"- Git HEAD: `{git_revision(REPO_ROOT)}`", + f"- Summary CSV: `{summary_csv}`", + f"- Proposed runs: `{proposed_dir}`", + f"- w_topo=0 ablation: `{wtopo0_dir}`" if wtopo0_dir else "- w_topo=0 ablation: (not provided)", + f"- Baselines CSV: `{baselines_csv}`", + "", + "## Metrics source", + "", + "- Proposed / ablation: `evaluate_paper_protocol.py` → `logs/paper_metrics_test.json`", + " (keys: `mcc`, `gmean`, `wdist` from DBSCAN test evaluation)", + "- Baselines: `evaluate_paper_baselines.py` → `summary_baselines.csv`", + "", + "## summary_for_paper.csv", + "", + "| method | MCC | G-Mean | W-Dist | notes |", + "|--------|-----|--------|--------|-------|", + ] + + for row in rows: + method = row.get("method", "") + if not row.get("mcc_mean"): + lines.append(f"| {method} | — | — | — | pending |") + continue + lines.append( + f"| {method} | {format_pm(float(row['mcc_mean']), float(row['mcc_std']))} " + f"| {format_pm(float(row['gmean_mean']), float(row['gmean_std']))} " + f"| {format_pm(float(row['wdist_mean']), float(row['wdist_std']))} " + f"| {row.get('notes', '')} |" + ) + + def _seed_section(title: str, seeds: Sequence[SeedPaperMetrics]) -> list[str]: + if not seeds: + return [f"## {title}", "", "(not available)", ""] + out = [f"## {title}", ""] + for s in sorted(seeds, key=lambda x: (x.seed is None, x.seed or 0)): + hparam = "" + if s.dbscan_eps is not None and s.dbscan_min_samples is not None: + hparam = f" eps={s.dbscan_eps}, min_samples={s.dbscan_min_samples}" + seed_label = s.seed if s.seed is not None else "?" + out.append( + f"- seed {seed_label}: `{s.run_dir}` — " + f"MCC={s.mcc:.4f}, G-Mean={s.gmean:.4f}, W-Dist={s.wdist:.4f}{hparam}" + ) + out.append("") + return out + + lines.extend(_seed_section("Proposed per-seed (w_topo > 0)", proposed_seeds)) + lines.extend(_seed_section("Ablation per-seed (w_topo = 0)", wtopo0_seeds)) + + if warnings: + lines.append("## Warnings") + lines.append("") + for w in warnings: + lines.append(f"- {w}") + lines.append("") + + path.write_text("\n".join(lines) + "\n", encoding="utf-8") + + +def parse_args() -> argparse.Namespace: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument( + "--proposed-dir", + type=Path, + default=REPO_ROOT / "outputs" / "paper_reproduce_1week_tuned", + help="Main proposed runs (12ep × 5 seed, ellphi).", + ) + p.add_argument( + "--wtopo0-dir", + type=Path, + default=REPO_ROOT / "outputs" / "paper_wtopo0", + help="w_topo=0 ablation output tree (optional until Phase 2 completes).", + ) + p.add_argument( + "--baselines-csv", + type=Path, + default=REPO_ROOT / "outputs" / "paper_baselines" / "summary_baselines.csv", + ) + p.add_argument( + "--out-dir", + type=Path, + default=None, + help="Write summary_for_paper.csv and MANIFEST.md here (default: --proposed-dir).", + ) + p.add_argument("--seeds", type=int, nargs="+", default=PAPER_SEEDS) + p.add_argument( + "--skip-wtopo0", + action="store_true", + help="Do not require w_topo=0 ablation row.", + ) + p.add_argument( + "--require-proposed", + action="store_true", + help="Exit non-zero if proposed runs are missing.", + ) + return p.parse_args() + + +def main() -> int: + args = parse_args() + proposed_dir = args.proposed_dir.resolve() + wtopo0_dir = args.wtopo0_dir.resolve() if args.wtopo0_dir else None + baselines_csv = args.baselines_csv.resolve() + out_dir = (args.out_dir or proposed_dir).resolve() + warnings: list[str] = [] + + proposed_row, proposed_seeds, w1 = load_proposed_row( + proposed_dir, + method="proposed_w_topo_pos", + default_notes="proposed; ellphi DBSCAN test metrics; val-tuned eps/min_samples per seed", + expected_seeds=args.seeds, + ) + warnings.extend(w1) + + wtopo0_row: dict[str, Any] | None = None + wtopo0_seeds: list[SeedPaperMetrics] = [] + if not args.skip_wtopo0 and wtopo0_dir is not None: + wtopo0_row, wtopo0_seeds, w2 = load_proposed_row( + wtopo0_dir, + method="proposed_w_topo_0", + default_notes="ablation w_topo=0; ellphi DBSCAN test metrics; val-tuned eps/min_samples per seed", + expected_seeds=args.seeds, + ) + warnings.extend(w2) + elif not args.skip_wtopo0: + warnings.append("proposed_w_topo_0: --wtopo0-dir not set") + + baseline_rows = load_baselines_rows(baselines_csv) + baseline_by_method = {r["method"]: r for r in baseline_rows} + + rows_by_method: dict[str, dict[str, Any]] = {} + if proposed_row: + rows_by_method["proposed_w_topo_pos"] = proposed_row + if wtopo0_row: + rows_by_method["proposed_w_topo_0"] = wtopo0_row + for method in ("euclidean_dbscan", "isolation_forest", "lof", "adbscan"): + if method in baseline_by_method: + rows_by_method[method] = baseline_by_method[method] + + ordered_rows: list[dict[str, Any]] = [] + for method in METHOD_ORDER: + if method in rows_by_method: + ordered_rows.append(rows_by_method[method]) + else: + warnings.append(f"metrics pending: {method}") + ordered_rows.append({"method": method, "notes": "pending"}) + + summary_path = out_dir / "summary_for_paper.csv" + manifest_path = out_dir / "MANIFEST.md" + write_summary_csv(summary_path, ordered_rows) + write_manifest( + manifest_path, + summary_csv=summary_path, + proposed_dir=proposed_dir, + wtopo0_dir=wtopo0_dir if not args.skip_wtopo0 else None, + baselines_csv=baselines_csv, + rows=ordered_rows, + proposed_seeds=proposed_seeds, + wtopo0_seeds=wtopo0_seeds, + warnings=warnings, + ) + + print(f"Wrote {summary_path}") + print(f"Wrote {manifest_path}") + for row in ordered_rows: + if row.get("mcc_mean"): + print( + f" {row['method']}: MCC={float(row['mcc_mean']):.4f}±{float(row['mcc_std']):.4f} " + f"G-Mean={float(row['gmean_mean']):.4f}±{float(row['gmean_std']):.4f} " + f"W-Dist={float(row['wdist_mean']):.4f}±{float(row['wdist_std']):.4f}" + ) + else: + print(f" {row['method']}: pending") + + if warnings: + print("\nWarnings:") + for w in warnings: + print(f" - {w}") + + if args.require_proposed and proposed_row is None: + return 1 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/experiments/ellphi_postfix_autopipeline.sh b/experiments/ellphi_postfix_autopipeline.sh new file mode 100755 index 0000000..07b9db5 --- /dev/null +++ b/experiments/ellphi_postfix_autopipeline.sh @@ -0,0 +1,303 @@ +#!/usr/bin/env bash +# Wait for 12ep ellphi calibration (seed 42), then launch 5-seed paper main run (fixed 12ep). +# +# Usage (background): +# nohup ./experiments/ellphi_postfix_autopipeline.sh >> outputs/ellphi_postfix_calib/autopipeline.log 2>&1 & +# +set -euo pipefail + +REPO_ROOT="$(cd "$(dirname "$0")/.." && pwd)" +cd "$REPO_ROOT" + +CALIB_OUT="${CALIB_OUT:-outputs/ellphi_postfix_calib}" +PILOT_LOG="${CALIB_OUT}/driver_12ep_seed42.log" +STATE_JSON="${CALIB_OUT}/autopipeline_state.json" +DECISION_JSON="${CALIB_OUT}/epoch_decision.json" +POLL_SEC="${POLL_SEC:-300}" +SEEDS=(42 123 456 789 1024) +MAIN_EPOCHS="${MAIN_EPOCHS:-12}" +MAIN_OUT="${MAIN_OUT:-outputs/paper_reproduce_1week_tuned}" +MAIN_CONFIG="${MAIN_CONFIG:-reproduce}" + +log() { echo "[$(date -Iseconds)] $*"; } + +write_state() { + local phase="$1" + local detail="${2:-}" + python3 - "$phase" "$detail" "$STATE_JSON" <<'PY' +import json, sys +from datetime import datetime, timezone +phase, detail, path = sys.argv[1:4] +state = {} +try: + with open(path) as f: + state = json.load(f) +except FileNotFoundError: + pass +state.update({ + "updated_utc": datetime.now(timezone.utc).isoformat(), + "phase": phase, + "detail": detail, +}) +with open(path, "w") as f: + json.dump(state, f, indent=2, ensure_ascii=False) + f.write("\n") +PY +} + +pilot_running() { + pgrep -f "run_backend_multiseed.py.*out-base ${CALIB_OUT}" >/dev/null 2>&1 \ + || pgrep -f "run_backend_multiseed.py.*${CALIB_OUT}" >/dev/null 2>&1 +} + +pilot_done() { + [[ -f "$PILOT_LOG" ]] && grep -q '^\[DONE\]' "$PILOT_LOG" +} + +pilot_failed() { + if pilot_running; then + return 1 + fi + if pilot_done; then + return 1 + fi + if [[ -f "$PILOT_LOG" ]] && grep -q 'Traceback' "$PILOT_LOG"; then + return 0 + fi + # No process, no DONE — treat as failure if log exists and has START + if [[ -f "$PILOT_LOG" ]] && grep -q '^\[START\]' "$PILOT_LOG"; then + return 0 + fi + return 1 +} + +wait_for_pilot() { + write_state "waiting_pilot" "polling every ${POLL_SEC}s" + log "Waiting for 12ep calibration in ${CALIB_OUT} ..." + while true; do + if pilot_done; then + log "Pilot completed ([DONE] in log)." + write_state "pilot_completed" + return 0 + fi + if pilot_failed; then + log "ERROR: Pilot exited without [DONE]. See ${PILOT_LOG}" + write_state "pilot_failed" "check driver log" + exit 1 + fi + if pilot_running; then + progress="$(python3 - "$PILOT_LOG" <<'PY' || true +import re, sys +t = open(sys.argv[1]).read() +m = list(re.finditer(r'Epoch (\d+):.*?(\d+)/71', t)) +if m: + print(f"epoch {m[-1].group(1)} batch {m[-1].group(2)}/71") +else: + n = len(re.findall(r'Epoch \d+: Val MCC=', t)) + if n: + print(f"completed_epochs={n}") + else: + print("epoch 1 starting") +PY +)" + log "Pilot still running (${progress})" + evaluate_pilot_if_done + else + log "Pilot process not found yet; waiting ..." + fi + sleep "$POLL_SEC" + done +} + +decide_epochs() { + python3 - "$CALIB_OUT" "$DECISION_JSON" <<'PY' +import csv +import json +import sys +from pathlib import Path + +calib = Path(sys.argv[1]) +out_path = Path(sys.argv[2]) + +# Find metrics from progress_summary or latest run dir +metrics_path = None +progress = calib / "progress_summary.csv" +if progress.exists(): + rows = list(csv.DictReader(progress.open())) + if rows: + run_dir = Path(rows[-1]["run_dir"]) + candidate = run_dir / "logs" / "metrics.csv" + if candidate.exists(): + metrics_path = candidate + +if metrics_path is None: + candidates = sorted(calib.glob("backend_ellphi_seed42_*/logs/metrics.csv")) + if candidates: + metrics_path = candidates[-1] + +if metrics_path is None or not metrics_path.exists(): + raise SystemExit("metrics.csv not found for pilot") + +rows = list(csv.DictReader(metrics_path.open())) +if len(rows) < 3: + raise SystemExit(f"Too few epochs in {metrics_path} ({len(rows)} rows)") + +def mcc_at(epoch: int) -> float | None: + for r in rows: + if int(r["epoch"]) == epoch: + return float(r["val_mcc"]) + return None + +best_row = max(rows, key=lambda r: float(r["val_mcc"])) +best_ep = int(best_row["epoch"]) +best_mcc = float(best_row["val_mcc"]) + +m10 = mcc_at(10) +m12 = mcc_at(12) if mcc_at(12) is not None else float(rows[-1]["val_mcc"]) +last_ep = int(rows[-1]["epoch"]) + +gain_10_12 = None +if m10 is not None and m12 is not None: + gain_10_12 = m12 - m10 + +use_20 = False +reason = [] +if gain_10_12 is not None and gain_10_12 < 0.002 and best_ep <= 11: + if m10 is not None and m10 >= 0.99 * best_mcc: + use_20 = True + reason.append("gain_10_12<0.002, best<=11, ep10>=99%best") +if last_ep < 10: + use_20 = False + reason = [f"only {last_ep} epochs logged; default 30"] + +chosen = 20 if use_20 else 30 +decision = { + "chosen_epochs": chosen, + "criteria": { + "gain_10_12": gain_10_12, + "best_epoch": best_ep, + "best_val_mcc": best_mcc, + "val_mcc_ep10": m10, + "val_mcc_ep12_or_last": m12, + "last_logged_epoch": last_ep, + }, + "reason": reason or ["default 30 (still improving or criteria not met)"], + "metrics_path": str(metrics_path), +} +out_path.write_text(json.dumps(decision, indent=2, ensure_ascii=False) + "\n") +print(chosen) +PY +} + +launch_main() { + local epochs="$1" + local out_base="${MAIN_OUT}" + mkdir -p "$out_base" + + if [[ -f "${out_base}/progress_summary.csv" ]]; then + completed="$(python3 - "$out_base/progress_summary.csv" <<'PY' +import csv, sys +rows = list(csv.DictReader(open(sys.argv[1]))) +print(sum(1 for r in rows if r.get("backend") == "ellphi")) +PY +)" + if [[ "$completed" -ge 5 ]]; then + log "Main run already has ${completed} ellphi rows in ${out_base}; skipping." + write_state "main_already_complete" "$out_base" + evaluate_all_runs "$out_base" + return 0 + fi + fi + + if pgrep -f "run_backend_multiseed.py.*out-base ${out_base}" >/dev/null 2>&1; then + log "Main run already in progress under ${out_base}; skipping duplicate launch." + write_state "main_running" "$out_base" + return 0 + fi + + log "Launching main: ${epochs}ep × 5 seeds → ${out_base}" + write_state "main_starting" "${epochs}ep ${out_base}" + + uv run python experiments/run_backend_multiseed.py \ + --base-config "$MAIN_CONFIG" \ + --epochs "$epochs" \ + --seeds "${SEEDS[@]}" \ + --backends ellphi \ + --out-base "$out_base" \ + 2>&1 | tee -a "${out_base}/driver_main.log" + + write_state "main_completed" "$out_base" + log "Main run finished. See ${out_base}/backend_stats.csv" + evaluate_all_runs "$out_base" +} + +evaluate_run_paper_protocol() { + local run_dir="$1" + local base_cfg="${2:-reproduce}" + log "Paper DBSCAN eval: ${run_dir} (config=${base_cfg})" + uv run python experiments/evaluate_paper_protocol.py \ + --run-dir "$run_dir" \ + --base-config "$base_cfg" \ + --split val + uv run python experiments/evaluate_paper_protocol.py \ + --run-dir "$run_dir" \ + --base-config "$base_cfg" \ + --split test +} + +evaluate_all_runs() { + local out_base="$1" + local base_cfg="${2:-$MAIN_CONFIG}" + local progress="${out_base}/progress_summary.csv" + [[ -f "$progress" ]] || return 0 + python3 - "$progress" <<'PY' | while read -r rd; do +import csv, sys +for row in csv.DictReader(open(sys.argv[1])): + if row.get("backend") == "ellphi": + print(row["run_dir"]) +PY + evaluate_run_paper_protocol "$rd" "$base_cfg" + done +} + +evaluate_pilot_if_done() { + if ! pilot_done; then + return 0 + fi + evaluate_all_runs "$CALIB_OUT" "reproduce" +} + +main() { + log "=== ellphi postfix autopipeline start (repo=${REPO_ROOT}) ===" + if ! pilot_done; then + wait_for_pilot + else + log "Pilot already complete; skipping wait." + write_state "pilot_completed" "pre-existing" + fi + + evaluate_pilot_if_done + + epochs="$MAIN_EPOCHS" + log "Paper main: fixed ${epochs}ep → ${MAIN_OUT} (config=${MAIN_CONFIG})" + python3 - "$epochs" "$DECISION_JSON" <<'PY' +import json, sys +from datetime import datetime, timezone +epochs, path = sys.argv[1:3] +decision = { + "chosen_epochs": int(epochs), + "reason": ["fixed 12ep for paper reproduce_1week_tuned (override MAIN_EPOCHS to change)"], + "decided_utc": datetime.now(timezone.utc).isoformat(), +} +with open(path, "w") as f: + json.dump(decision, f, indent=2, ensure_ascii=False) + f.write("\n") +PY + write_state "epoch_decided" "${epochs}ep ${MAIN_OUT}" + + launch_main "$epochs" + write_state "pipeline_finished" "${epochs}ep ${MAIN_OUT}" + log "=== autopipeline finished ===" +} + +main "$@" diff --git a/experiments/evaluate_paper_baselines.py b/experiments/evaluate_paper_baselines.py new file mode 100644 index 0000000..3deb7aa --- /dev/null +++ b/experiments/evaluate_paper_baselines.py @@ -0,0 +1,514 @@ +#!/usr/bin/env python3 +""" +Paper-aligned baseline evaluation (Phase 3). + +Same data split as ``evaluate_paper_protocol.py`` (``configs/reproduce.yaml``, +seeds 42, 123, 456, 789, 1024): tune hyperparameters on validation clouds, +report MCC / G-Mean / W-Dist on test via ``compute_recall_specificity_gmean_mcc_wdist``. + +Methods: + 1. Euclidean DBSCAN (sklearn on raw coordinates) + 2. Isolation Forest + 3. Local Outlier Factor (LOF) + 4. ADBSCAN (local PCA ellipses + Mahalanobis ``apply_anisotropic_dbscan``) + +Usage:: + + uv run python experiments/evaluate_paper_baselines.py \\ + --out-dir outputs/paper_baselines +""" + +from __future__ import annotations + +import argparse +import csv +import json +from collections.abc import Callable, Sequence +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Any + +import numpy as np +import torch +from sklearn.cluster import DBSCAN +from sklearn.ensemble import IsolationForest +from sklearn.neighbors import LocalOutlierFactor +from tqdm import tqdm + +from experiments.evaluate_paper_protocol import ( + CloudMetrics, + _aggregate_cloud_metrics, + _valid_clean_inliers, + build_split_loader, + dbscan_labels_to_outlier_pred, + evaluate_cloud_dbscan, +) +from tda_ml.config import deep_update, load_config +from tda_ml.metrics import compute_recall_specificity_gmean_mcc_wdist +from tda_ml.numerical_eps import EIGENVALUE_FLOOR, PCA_RIDGE_EPS +from tda_ml.supervised_diagnostics import git_revision + +REPO_ROOT = Path(__file__).resolve().parents[1] +PAPER_SEEDS = [42, 123, 456, 789, 1024] +LOCAL_PCA_K = 10 + +DEFAULT_EPS_VALUES = list(np.linspace(0.15, 1.5, 15)) +DEFAULT_MIN_SAMPLES_VALUES = [3, 5, 7, 10, 15] +DEFAULT_CONTAMINATION_VALUES = [0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.20] +DEFAULT_LOF_N_NEIGHBORS = [5, 10, 15, 20, 30] + + +@dataclass +class CloudSample: + points: np.ndarray + labels_gt: np.ndarray + clean_pc: np.ndarray + adbscan_params: np.ndarray | None = None + + +@dataclass +class SeedResult: + seed: int + method: str + hparams: dict[str, Any] + val_mcc: float + test_recall: float + test_specificity: float + test_gmean: float + test_mcc: float + test_wdist: float + n_test_clouds: int + + +def local_pca_ellipse_params(points: np.ndarray, k: int = LOCAL_PCA_K) -> np.ndarray: + """ + Baseline ellipse parameters from local PCA only (no learned deltas). + + Matches ``DecoupledGeometricEncoder`` + zero MLP corrections: + ``[a, b, theta]`` per point, shape ``(N, 3)``. + """ + if points.ndim != 2 or points.shape[1] != 2: + raise ValueError(f"points must be (N, 2); got {points.shape}") + n = points.shape[0] + if n < k: + raise ValueError(f"Need at least k={k} points; got n={n}") + + x = torch.from_numpy(points.astype(np.float32)).unsqueeze(0) # 1, N, 2 + dist_sq = torch.cdist(x, x, p=2) ** 2 + _, idx = torch.topk(-dist_sq, k=k, dim=-1) + + batch_idx = torch.arange(1).view(1, 1, 1).expand(1, n, k) + flat_x = x.view(n, 2) + flat_neighbors = flat_x[idx.view(1, -1) + (batch_idx.view(1, -1) * n), :] + neighbors = flat_neighbors.view(1, n, k, 2) + + relative_coords = neighbors - x.unsqueeze(2) + mean_neighbor = relative_coords.mean(dim=2, keepdim=True) + centered = relative_coords - mean_neighbor + cov = torch.matmul(centered.transpose(-1, -2), centered) / (k - 1) + + eye2 = torch.eye(2, dtype=torch.float32) + e, v = torch.linalg.eigh(cov.float() + eye2 * PCA_RIDGE_EPS) + v1 = v[:, :, :, 1] + base_angle = torch.atan2(v1[:, :, 1], v1[:, :, 0]) + + base_axes = torch.sqrt(torch.clamp(e, min=EIGENVALUE_FLOOR)) + base_axes = torch.flip(base_axes, dims=[-1]) + base_axes = base_axes / (base_axes.max(dim=-1, keepdim=True)[0] + EIGENVALUE_FLOOR) + + params = torch.cat([base_axes, base_angle.unsqueeze(-1)], dim=-1) + return params.squeeze(0).numpy() + + +def load_clouds(config: dict[str, Any], split: str, device: torch.device) -> list[CloudSample]: + loader = build_split_loader(config, split, device) + clouds: list[CloudSample] = [] + for data, labels, clean_pc in loader: + data_np = data.numpy() + labels_np = labels.numpy() + clean_np = clean_pc.numpy() + for b in range(data_np.shape[0]): + points = data_np[b] + params = local_pca_ellipse_params(points) + clouds.append( + CloudSample( + points=points, + labels_gt=labels_np[b], + clean_pc=clean_np[b], + adbscan_params=params, + ) + ) + return clouds + + +def cloud_metrics_from_pred( + labels_gt: np.ndarray, + pred: np.ndarray, + points: np.ndarray, + clean_pc: np.ndarray, +) -> CloudMetrics: + gt_inliers = _valid_clean_inliers(clean_pc) + recall, specificity, gmean, mcc, wdist = compute_recall_specificity_gmean_mcc_wdist( + labels_gt, + pred, + points=points, + gt_inliers=gt_inliers, + ) + return CloudMetrics(recall, specificity, gmean, mcc, wdist) + + +def evaluate_euclidean_dbscan( + cloud: CloudSample, + *, + eps: float, + min_samples: int, +) -> CloudMetrics: + db_labels = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(cloud.points) + pred = dbscan_labels_to_outlier_pred(db_labels) + return cloud_metrics_from_pred(cloud.labels_gt, pred, cloud.points, cloud.clean_pc) + + +def evaluate_isolation_forest( + cloud: CloudSample, + *, + contamination: float, + random_state: int, +) -> CloudMetrics: + iso = IsolationForest(contamination=contamination, random_state=random_state) + sk_pred = iso.fit_predict(cloud.points) + pred = (sk_pred == -1).astype(np.int64) + return cloud_metrics_from_pred(cloud.labels_gt, pred, cloud.points, cloud.clean_pc) + + +def evaluate_lof( + cloud: CloudSample, + *, + n_neighbors: int, + contamination: float, +) -> CloudMetrics: + lof = LocalOutlierFactor( + n_neighbors=n_neighbors, + contamination=contamination, + novelty=False, + ) + sk_pred = lof.fit_predict(cloud.points) + pred = (sk_pred == -1).astype(np.int64) + return cloud_metrics_from_pred(cloud.labels_gt, pred, cloud.points, cloud.clean_pc) + + +def evaluate_adbscan( + cloud: CloudSample, + *, + eps: float, + min_samples: int, +) -> CloudMetrics: + if cloud.adbscan_params is None: + raise ValueError("adbscan_params missing") + return evaluate_cloud_dbscan( + cloud.points, + cloud.adbscan_params, + cloud.labels_gt, + cloud.clean_pc, + eps=eps, + min_samples=min_samples, + backend="mahalanobis", + ) + + +def grid_search_clouds( + clouds: Sequence[CloudSample], + evaluate_fn: Callable[..., CloudMetrics], + param_combos: Sequence[dict[str, Any]], + *, + desc: str, +) -> tuple[dict[str, Any], float, list[dict[str, Any]]]: + best_mcc = -1.0 + best_params = dict(param_combos[0]) + grid_log: list[dict[str, Any]] = [] + + for params in tqdm(param_combos, desc=desc, leave=False): + per_cloud: list[CloudMetrics] = [] + error: str | None = None + for cloud in clouds: + try: + per_cloud.append(evaluate_fn(cloud, **params)) + except Exception as exc: # noqa: BLE001 — skip invalid hparam combos + error = str(exc) + per_cloud = [] + break + if error is not None: + grid_log.append({**params, "error": error}) + continue + _, _, _, mcc, _ = _aggregate_cloud_metrics(per_cloud) + grid_log.append({**params, "mean_mcc": mcc, "n_clouds": len(per_cloud)}) + if mcc > best_mcc: + best_mcc = mcc + best_params = dict(params) + + if best_mcc < 0: + raise RuntimeError(f"Grid search failed for all hparams ({desc}); log={grid_log[:5]}") + return best_params, best_mcc, grid_log + + +def dbscan_param_combos( + eps_values: Sequence[float], + min_samples_values: Sequence[int], +) -> list[dict[str, Any]]: + return [ + {"eps": float(eps), "min_samples": int(ms)} + for eps in eps_values + for ms in min_samples_values + ] + + +def contamination_param_combos(contamination_values: Sequence[float]) -> list[dict[str, Any]]: + return [{"contamination": float(c)} for c in contamination_values] + + +def lof_param_combos( + n_neighbors_values: Sequence[int], + contamination_values: Sequence[float], +) -> list[dict[str, Any]]: + return [ + {"n_neighbors": int(n), "contamination": float(c)} + for n in n_neighbors_values + for c in contamination_values + ] + + +def evaluate_method_on_seed( + method: str, + config: dict[str, Any], + seed: int, + device: torch.device, + *, + eps_values: Sequence[float], + min_samples_values: Sequence[int], + contamination_values: Sequence[float], + lof_n_neighbors_values: Sequence[int], +) -> tuple[SeedResult, list[dict[str, Any]]]: + val_clouds = load_clouds(config, "val", device) + test_clouds = load_clouds(config, "test", device) + + if method == "euclidean_dbscan": + combos = dbscan_param_combos(eps_values, min_samples_values) + best_params, val_mcc, grid_log = grid_search_clouds( + val_clouds, + evaluate_euclidean_dbscan, + combos, + desc=f"seed{seed} euclidean_dbscan val", + ) + test_fn: Callable[..., CloudMetrics] = evaluate_euclidean_dbscan + elif method == "isolation_forest": + combos = [ + {"contamination": float(c), "random_state": seed} + for c in contamination_values + ] + best_params, val_mcc, grid_log = grid_search_clouds( + val_clouds, + evaluate_isolation_forest, + combos, + desc=f"seed{seed} isolation_forest val", + ) + test_fn = evaluate_isolation_forest + elif method == "lof": + combos = lof_param_combos(lof_n_neighbors_values, contamination_values) + best_params, val_mcc, grid_log = grid_search_clouds( + val_clouds, + evaluate_lof, + combos, + desc=f"seed{seed} lof val", + ) + test_fn = evaluate_lof + elif method == "adbscan": + combos = dbscan_param_combos(eps_values, min_samples_values) + best_params, val_mcc, grid_log = grid_search_clouds( + val_clouds, + evaluate_adbscan, + combos, + desc=f"seed{seed} adbscan val", + ) + test_fn = evaluate_adbscan + else: + raise ValueError(f"Unknown method: {method!r}") + + per_test: list[CloudMetrics] = [] + for cloud in test_clouds: + per_test.append(test_fn(cloud, **best_params)) + recall, specificity, gmean, mcc, wdist = _aggregate_cloud_metrics(per_test) + + return SeedResult( + seed=seed, + method=method, + hparams=best_params, + val_mcc=val_mcc, + test_recall=recall, + test_specificity=specificity, + test_gmean=gmean, + test_mcc=mcc, + test_wdist=wdist, + n_test_clouds=len(per_test), + ), grid_log + + +def config_for_seed(base_config: str, seed: int) -> dict[str, Any]: + cfg = load_config(base_config, project_root=REPO_ROOT) + return deep_update(cfg, {"data": {"seed": int(seed)}}) + + +def sample_std(values: Sequence[float]) -> float: + arr = np.asarray(values, dtype=np.float64) + if arr.size < 2: + return 0.0 + return float(np.std(arr, ddof=1)) + + +def write_summary_csv( + out_path: Path, + rows: list[dict[str, Any]], +) -> None: + fieldnames = [ + "method", + "mcc_mean", + "mcc_std", + "gmean_mean", + "gmean_std", + "wdist_mean", + "wdist_std", + "notes", + ] + out_path.parent.mkdir(parents=True, exist_ok=True) + with out_path.open("w", newline="", encoding="utf-8") as f: + writer = csv.DictWriter(f, fieldnames=fieldnames) + writer.writeheader() + for row in rows: + writer.writerow({k: row.get(k, "") for k in fieldnames}) + + +def aggregate_method_results(seed_results: Sequence[SeedResult]) -> dict[str, Any]: + mccs = [r.test_mcc for r in seed_results] + gmeans = [r.test_gmean for r in seed_results] + wdist = [r.test_wdist for r in seed_results] + method = seed_results[0].method + return { + "method": method, + "mcc_mean": float(np.mean(mccs)), + "mcc_std": sample_std(mccs), + "gmean_mean": float(np.mean(gmeans)), + "gmean_std": sample_std(gmeans), + "wdist_mean": float(np.mean(wdist)), + "wdist_std": sample_std(wdist), + "notes": ( + f"5 seeds; val hparam selection by max mean cloud MCC; " + f"test n_clouds={seed_results[0].n_test_clouds} per seed" + ), + } + + +METHOD_ORDER = [ + "euclidean_dbscan", + "isolation_forest", + "lof", + "adbscan", +] + + +def parse_args() -> argparse.Namespace: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--base-config", type=str, default="reproduce") + p.add_argument("--out-dir", type=Path, default=REPO_ROOT / "outputs" / "paper_baselines") + p.add_argument("--seeds", type=int, nargs="+", default=PAPER_SEEDS) + p.add_argument( + "--methods", + nargs="+", + choices=METHOD_ORDER, + default=METHOD_ORDER, + ) + p.add_argument("--eps-values", type=float, nargs="+", default=DEFAULT_EPS_VALUES) + p.add_argument("--min-samples-values", type=int, nargs="+", default=DEFAULT_MIN_SAMPLES_VALUES) + p.add_argument("--contamination-values", type=float, nargs="+", default=DEFAULT_CONTAMINATION_VALUES) + p.add_argument("--lof-n-neighbors", type=int, nargs="+", default=DEFAULT_LOF_N_NEIGHBORS) + return p.parse_args() + + +def main() -> int: + args = parse_args() + out_dir = args.out_dir.resolve() + out_dir.mkdir(parents=True, exist_ok=True) + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + + manifest = { + "source_revision": git_revision(REPO_ROOT), + "base_config": args.base_config, + "seeds": list(args.seeds), + "methods": list(args.methods), + "selection": "max mean cloud MCC on validation split", + "metrics": "compute_recall_specificity_gmean_mcc_wdist", + "local_pca_k": LOCAL_PCA_K, + "grids": { + "eps_values": list(args.eps_values), + "min_samples_values": list(args.min_samples_values), + "contamination_values": list(args.contamination_values), + "lof_n_neighbors": list(args.lof_n_neighbors), + }, + } + + all_seed_results: dict[str, list[SeedResult]] = {m: [] for m in args.methods} + + for seed in args.seeds: + config = config_for_seed(args.base_config, seed) + seed_dir = out_dir / f"seed{seed}" + seed_dir.mkdir(parents=True, exist_ok=True) + + for method in args.methods: + print(f"\n=== seed={seed} method={method} ===") + result, grid_log = evaluate_method_on_seed( + method, + config, + seed, + device, + eps_values=args.eps_values, + min_samples_values=args.min_samples_values, + contamination_values=args.contamination_values, + lof_n_neighbors_values=args.lof_n_neighbors, + ) + all_seed_results[method].append(result) + + payload = { + **asdict(result), + "grid_log": grid_log, + } + out_json = seed_dir / f"{method}.json" + out_json.write_text(json.dumps(payload, indent=2) + "\n") + print( + f" val_mcc={result.val_mcc:.4f} test_mcc={result.test_mcc:.4f} " + f"test_gmean={result.test_gmean:.4f} test_wdist={result.test_wdist:.4f} " + f"hparams={result.hparams}" + ) + + summary_rows = [ + aggregate_method_results(all_seed_results[m]) + for m in args.methods + if all_seed_results[m] + ] + summary_path = out_dir / "summary_baselines.csv" + write_summary_csv(summary_path, summary_rows) + + manifest["summary_csv"] = str(summary_path) + manifest["per_seed_json"] = str(out_dir / "seed{seed}/{method}.json") + manifest_path = out_dir / "MANIFEST_baselines.json" + manifest_path.write_text(json.dumps(manifest, indent=2) + "\n") + + print(f"\nWrote {summary_path}") + print(f"Wrote {manifest_path}") + for row in summary_rows: + print( + f" {row['method']}: MCC={row['mcc_mean']:.4f}±{row['mcc_std']:.4f} " + f"G-Mean={row['gmean_mean']:.4f}±{row['gmean_std']:.4f} " + f"W-Dist={row['wdist_mean']:.4f}±{row['wdist_std']:.4f}" + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/experiments/evaluate_paper_protocol.py b/experiments/evaluate_paper_protocol.py new file mode 100644 index 0000000..cbb67e4 --- /dev/null +++ b/experiments/evaluate_paper_protocol.py @@ -0,0 +1,407 @@ +#!/usr/bin/env python3 +""" +Paper-aligned evaluation: ellphi DBSCAN inference + MCC / G-Mean / W-Dist. + +Trainer ``val_mcc`` uses sigmoid(logit) > threshold; the paper reports metrics +after DBSCAN on the learned ellphi precomputed distance matrix. + +Usage:: + + uv run python experiments/evaluate_paper_protocol.py \\ + --run-dir outputs/paper_reproduce_1week_tuned/backend_ellphi_seed42_* \\ + --base-config reproduce \\ + --split val + + uv run python experiments/evaluate_paper_protocol.py \\ + --run-dir ... --split test \\ + --dbscan-hparams outputs/.../logs/dbscan_hparams.json +""" + +from __future__ import annotations + +import argparse +import json +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Any, Iterator + +import numpy as np +import torch + +from tda_ml.checkpoint_io import extract_model_state_dict, load_torch_checkpoint +from tda_ml.config import deep_update, load_config, model_kwargs_from_config +from tda_ml.data_loader import NoisyMNISTDataset, create_data_loader +from tda_ml.dbscan import apply_anisotropic_dbscan +from tda_ml.metrics import compute_recall_specificity_gmean_mcc_wdist +from tda_ml.models import AnisotropicOutlierClassifier +from tda_ml.seed_utils import set_global_seed +from tda_ml.supervised_diagnostics import git_revision + +REPO_ROOT = Path(__file__).resolve().parents[1] + + +@dataclass +class CloudMetrics: + recall: float + specificity: float + gmean: float + mcc: float + wdist: float + + +@dataclass +class SplitMetrics: + split: str + n_clouds: int + recall: float + specificity: float + gmean: float + mcc: float + wdist: float + dbscan_eps: float | None = None + dbscan_min_samples: int | None = None + backend: str = "ellphi" + + +def _valid_clean_inliers(clean_pc: np.ndarray) -> np.ndarray: + mask = np.abs(clean_pc).sum(axis=1) > 1e-6 + return clean_pc[mask] + + +def dbscan_labels_to_outlier_pred(labels: np.ndarray) -> np.ndarray: + """DBSCAN noise (-1) -> outlier (1); clustered points -> inlier (0).""" + return (labels == -1).astype(np.int64) + + +def evaluate_cloud_dbscan( + points: np.ndarray, + params: np.ndarray, + labels_gt: np.ndarray, + clean_pc: np.ndarray, + *, + eps: float, + min_samples: int, + backend: str = "ellphi", + metric: str = "max", +) -> CloudMetrics: + db_labels = apply_anisotropic_dbscan( + points, + params, + eps=eps, + min_samples=min_samples, + metric=metric, + backend=backend, + ) + pred = dbscan_labels_to_outlier_pred(db_labels) + gt_inliers = _valid_clean_inliers(clean_pc) + recall, specificity, gmean, mcc, wdist = compute_recall_specificity_gmean_mcc_wdist( + labels_gt, + pred, + points=points, + gt_inliers=gt_inliers, + ) + return CloudMetrics(recall, specificity, gmean, mcc, wdist) + + +def _aggregate_cloud_metrics(rows: list[CloudMetrics]) -> tuple[float, float, float, float, float]: + if not rows: + raise ValueError("No clouds to aggregate") + return ( + float(np.mean([r.recall for r in rows])), + float(np.mean([r.specificity for r in rows])), + float(np.mean([r.gmean for r in rows])), + float(np.mean([r.mcc for r in rows])), + float(np.mean([r.wdist for r in rows])), + ) + + +def build_split_loader(config: dict[str, Any], split: str, device: torch.device): + data_cfg = config["data"] + seed = int(data_cfg.get("seed", 42)) + set_global_seed(seed, deterministic_algorithms=False) + + train_size = int(data_cfg.get("train_size", 4500)) + val_size = int(data_cfg.get("val_size", 500)) + test_size = int(data_cfg.get("test_size", 1000)) + generator = torch.Generator().manual_seed(seed) + full_train_indices = torch.randperm(60000, generator=generator)[: train_size + val_size] + train_indices = full_train_indices[:train_size] + val_indices = full_train_indices[train_size:] + test_indices = torch.randperm(10000, generator=generator)[:test_size] + + num_workers = int(data_cfg.get("num_workers", 0)) + pin_memory = bool(data_cfg.get("pin_memory", device.type == "cuda")) + batch_size = int(data_cfg.get("batch_size", 64)) + + if split == "val": + indices = val_indices + train_flag = True + elif split == "test": + indices = test_indices + train_flag = False + else: + raise ValueError(f"split must be 'val' or 'test'; got {split!r}") + + dataset = NoisyMNISTDataset( + root=str(REPO_ROOT / "data"), + train=train_flag, + max_points=data_cfg["max_points"], + num_outliers=data_cfg["num_outliers"], + indices=indices, + deterministic=True, + noise_seed=seed, + preload=True, + ) + loader = create_data_loader( + dataset, + batch_size=batch_size, + shuffle=False, + num_workers=num_workers, + pin_memory=pin_memory, + persistent_workers=False, + prefetch_factor=2 if num_workers > 0 else None, + ) + return loader + + +def load_model_from_run(run_dir: Path, config: dict[str, Any], device: torch.device) -> AnisotropicOutlierClassifier: + ckpt_path = run_dir / "best_model.pth" + if not ckpt_path.is_file(): + raise FileNotFoundError(f"Missing checkpoint: {ckpt_path}") + model = AnisotropicOutlierClassifier(**model_kwargs_from_config(config)) + ckpt = load_torch_checkpoint(str(ckpt_path), map_location="cpu") + model.load_state_dict(extract_model_state_dict(ckpt), strict=False) + model.to(device) + model.eval() + return model + + +def iter_cloud_predictions( + model: AnisotropicOutlierClassifier, + loader, + device: torch.device, +) -> Iterator[tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]]: + with torch.no_grad(): + for data, labels, clean_pc in loader: + data = data.to(device, non_blocking=True) + _, params = model(data) + data_np = data.cpu().numpy() + params_np = params.cpu().numpy() + labels_np = labels.cpu().numpy() + clean_np = clean_pc.cpu().numpy() + batch_size = data_np.shape[0] + for b in range(batch_size): + yield ( + data_np[b], + params_np[b], + labels_np[b], + clean_np[b], + ) + + +def grid_search_dbscan( + clouds: list[tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]], + *, + eps_values: list[float], + min_samples_values: list[int], + backend: str = "ellphi", +) -> tuple[float, int, float]: + best_mcc = -1.0 + best_eps = eps_values[0] + best_min_samples = min_samples_values[0] + grid_log: list[dict[str, Any]] = [] + + for eps in eps_values: + for min_samples in min_samples_values: + per_cloud: list[CloudMetrics] = [] + for points, params, labels_gt, clean_pc in clouds: + try: + m = evaluate_cloud_dbscan( + points, + params, + labels_gt, + clean_pc, + eps=eps, + min_samples=min_samples, + backend=backend, + ) + per_cloud.append(m) + except Exception as exc: # noqa: BLE001 — log and skip bad hparams + grid_log.append( + { + "eps": eps, + "min_samples": min_samples, + "error": str(exc), + } + ) + per_cloud = [] + break + if not per_cloud: + continue + _, _, _, mcc, _ = _aggregate_cloud_metrics(per_cloud) + grid_log.append( + { + "eps": eps, + "min_samples": min_samples, + "mean_mcc": mcc, + "n_clouds": len(per_cloud), + } + ) + if mcc > best_mcc: + best_mcc = mcc + best_eps = eps + best_min_samples = min_samples + + if best_mcc < 0: + raise RuntimeError(f"DBSCAN grid search failed for all hparams; log={grid_log[:5]}") + return best_eps, best_min_samples, best_mcc + + +def evaluate_split( + run_dir: Path, + config: dict[str, Any], + split: str, + device: torch.device, + *, + eps: float | None = None, + min_samples: int | None = None, + eps_values: list[float] | None = None, + min_samples_values: list[int] | None = None, + backend: str = "ellphi", +) -> SplitMetrics: + loader = build_split_loader(config, split, device) + model = load_model_from_run(run_dir, config, device) + clouds = list(iter_cloud_predictions(model, loader, device)) + + if split == "val": + eps_values = eps_values or list(np.linspace(0.15, 1.5, 15)) + min_samples_values = min_samples_values or [3, 5, 7, 10, 15] + eps, min_samples, _ = grid_search_dbscan( + clouds, + eps_values=eps_values, + min_samples_values=min_samples_values, + backend=backend, + ) + else: + if eps is None or min_samples is None: + raise ValueError("test split requires eps and min_samples (from val tuning)") + + per_cloud: list[CloudMetrics] = [] + for points, params, labels_gt, clean_pc in clouds: + per_cloud.append( + evaluate_cloud_dbscan( + points, + params, + labels_gt, + clean_pc, + eps=float(eps), + min_samples=int(min_samples), + backend=backend, + ) + ) + recall, specificity, gmean, mcc, wdist = _aggregate_cloud_metrics(per_cloud) + return SplitMetrics( + split=split, + n_clouds=len(per_cloud), + recall=recall, + specificity=specificity, + gmean=gmean, + mcc=mcc, + wdist=wdist, + dbscan_eps=float(eps), + dbscan_min_samples=int(min_samples), + backend=backend, + ) + + +def load_run_config(run_dir: Path, base_config: str, seed: int | None) -> dict[str, Any]: + manifest_path = run_dir / "logs" / "run_manifest.json" + cfg = load_config(base_config, project_root=REPO_ROOT) + if manifest_path.is_file(): + manifest = json.loads(manifest_path.read_text()) + if seed is None: + seed = manifest.get("seed") + if seed is not None: + cfg = deep_update(cfg, {"data": {"seed": int(seed)}}) + return cfg + + +def parse_args() -> argparse.Namespace: + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--run-dir", type=Path, required=True) + p.add_argument("--base-config", type=str, default="reproduce") + p.add_argument("--split", choices=["val", "test"], required=True) + p.add_argument("--seed", type=int, default=None, help="Override data.seed (else run_manifest)") + p.add_argument("--backend", type=str, default="ellphi", choices=["ellphi", "mahalanobis"]) + p.add_argument( + "--dbscan-hparams", + type=Path, + default=None, + help="JSON with eps/min_samples for test split", + ) + p.add_argument("--out-json", type=Path, default=None, help="Write metrics JSON (default: run_dir/logs/)") + return p.parse_args() + + +def main() -> int: + args = parse_args() + run_dir = args.run_dir.resolve() + if not run_dir.is_dir(): + raise SystemExit(f"run-dir not found: {run_dir}") + + device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + config = load_run_config(run_dir, args.base_config, args.seed) + + eps = min_samples = None + if args.split == "test": + if args.dbscan_hparams is None: + args.dbscan_hparams = run_dir / "logs" / "dbscan_hparams.json" + hparams = json.loads(args.dbscan_hparams.read_text()) + eps = float(hparams["eps"]) + min_samples = int(hparams["min_samples"]) + + metrics = evaluate_split( + run_dir, + config, + args.split, + device, + eps=eps, + min_samples=min_samples, + backend=args.backend, + ) + + log_dir = run_dir / "logs" + log_dir.mkdir(parents=True, exist_ok=True) + + if args.split == "val": + hparams_path = log_dir / "dbscan_hparams.json" + hparams_path.write_text( + json.dumps( + { + "eps": metrics.dbscan_eps, + "min_samples": metrics.dbscan_min_samples, + "backend": metrics.backend, + "mean_val_mcc_dbscan": metrics.mcc, + "selection": "max mean cloud MCC on validation split", + }, + indent=2, + ) + + "\n" + ) + print(f"Saved DBSCAN hparams: {hparams_path}") + + out_path = args.out_json or log_dir / f"paper_metrics_{args.split}.json" + payload = { + "source_revision": git_revision(REPO_ROOT), + "run_dir": str(run_dir), + "split": args.split, + **asdict(metrics), + } + out_path.write_text(json.dumps(payload, indent=2) + "\n") + print(json.dumps(payload, indent=2)) + print(f"Wrote {out_path}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/experiments/run_backend_multiseed.py b/experiments/run_backend_multiseed.py index 1835699..4d9161a 100644 --- a/experiments/run_backend_multiseed.py +++ b/experiments/run_backend_multiseed.py @@ -267,13 +267,14 @@ def run_one( epochs: int, out_base: str, progress_csv: str, + w_topo: float | None = None, ) -> None: cfg = load_config(base_config_name, project_root=REPO_ROOT) config_id = f"backend_{backend}_seed{seed}" topo_cfg = cfg.get("model", {}).get("topology_loss", {}) ellphi_diff = bool(topo_cfg.get("ellphi_differentiable", True)) - overrides = { + overrides: dict[str, Any] = { "meta": {"config_id": config_id}, "training": {"epochs": epochs}, "data": {"seed": seed}, @@ -285,6 +286,8 @@ def run_one( }, "outputs": {"base_dir": out_base}, } + if w_topo is not None: + overrides["loss"] = {"w_topo": float(w_topo)} cfg = deep_update(cfg, overrides) before = set(glob.glob(os.path.join(out_base, f"{config_id}_*"))) @@ -346,6 +349,12 @@ def parse_args() -> argparse.Namespace: "progress_summary.csv." ), ) + p.add_argument( + "--w-topo", + type=float, + default=None, + help="Override loss.w_topo (e.g. 0 for topology ablation).", + ) return p.parse_args() @@ -387,6 +396,7 @@ def main() -> None: epochs=args.epochs, out_base=args.out_base, progress_csv=progress_csv, + w_topo=args.w_topo, ) write_backend_stats(progress_csv, stats_csv)