diff --git a/4-Infrastructure/shim/label_canary_theorems.py b/4-Infrastructure/shim/label_canary_theorems.py
new file mode 100644
index 00000000..e5062ba7
--- /dev/null
+++ b/4-Infrastructure/shim/label_canary_theorems.py
@@ -0,0 +1,259 @@
+#!/usr/bin/env python3
+"""Generate ground-truth multi-labels for the 42 canary theorems.
+
+Reads canary receipts + results, infers labels from Lean code and proof output.
+Outputs shared-data/pist_canary_ground_truth_labels.jsonl
+"""
+
+import json
+import os
+import re
+import sys
+from collections import Counter, defaultdict
+
+RESULTS_PATH = os.path.join(os.path.dirname(__file__), "../..",
+                            "shared-data/pist_canary_results.jsonl")
+RECEIPTS_PATH = os.path.join(os.path.dirname(__file__), "../..",
+                             "shared-data/pist_canary_receipts.jsonl")
+LABELS_PATH = os.path.join(os.path.dirname(__file__), "../..",
+                           "shared-data/pist_canary_ground_truth_labels.jsonl")
+
+
+def parse_domain(code: str) -> str:
+    """Infer theorem domain from Lean code."""
+    code_lower = code.lower()
+    if "list" in code_lower and ("map" in code_lower or "reverse" in code_lower or "length" in code_lower):
+        return "list"
+    if "set" in code_lower and ("union" in code_lower or "inter" in code_lower):
+        return "set"
+    if "prop" in code_lower and ("∧" in code or "∨" in code or "→" in code or "¬" in code):
+        return "logic"
+    if "≤" in code or "≥" in code or "order" in code_lower or "le_refl" in code or "le_trans" in code:
+        return "order"
+    if "ring" in code_lower or "mul_comm" in code or "*" in code:
+        return "algebra"
+    if "gcd" in code_lower or "nat." in code_lower or "+" in code or "omega" in code_lower:
+        return "arithmetic"
+    if "def " in code:
+        return "fixed_point"
+    if "induction" in code_lower:
+        return "induction"
+    return "other"
+
+
+def parse_proof_method(code: str) -> str:
+    """Infer primary proof method from proof script.
+    
+    Handles both full theorem statements and bare by-blocks.
+    """
+    # Extract the by-block if present
+    m = re.search(r'by\s+(.+?)$', code)
+    if m:
+        proof = m.group(1).strip()
+    else:
+        proof = code.strip()
+
+    if not proof:
+        return "unknown"
+
+    if "omega" in proof:
+        return "omega"
+    if "nlinarith" in proof:
+        return "nlinarith"
+    if "ring" in proof:
+        return "ring"
+    if "simp" in proof and "induction" not in proof:
+        return "simp"
+    if "induction" in proof:
+        return "induction"
+    if "rw " in proof or "rewrite" in proof:
+        return "rewrite"
+    if "rfl" in proof:
+        return "rfl"
+    if "cases" in proof:
+        return "cases"
+    if "apply" in proof or "exact" in proof:
+        return "apply"
+    if "intro" in proof:
+        return "intro"
+    if "calc" in proof:
+        return "calc"
+    if "native_decide" in proof:
+        return "native_decide"
+    if "sorry" in proof:
+        return "incomplete"
+    return "unknown"
+
+
+def parse_obstruction(stdout: str, stderr: str, status: str, code: str) -> str | None:
+    """Infer obstruction type from proof output."""
+    if status == "verified" or "verified" in status:
+        return None
+    combined = (stdout + " " + stderr).lower()
+    if "timeout" in status or "timeout" in combined:
+        return "timeout"
+    if "type mismatch" in combined or "typeMismatch" in combined:
+        return "type_mismatch"
+    if "failed to synthesize" in combined or "synthInstance" in combined:
+        return "missing_instance"
+    if "unknown identifier" in combined or "unknownConstant" in combined:
+        return "missing_lemma"
+    if "unsat" in combined or "contradiction" in combined:
+        return "unsatisfiable_goal"
+    if "omega" in combined and "could not prove" in combined:
+        return "nonlinear_arithmetic"
+    if "coercion" in combined:
+        return "coercion_gap"
+    if "simp" in combined and "did not simplify" in combined:
+        return "simplifier_gap"
+    if "stack overflow" in combined or "recursion" in combined:
+        return "recursion_limit"
+    if "no applicable" in combined:
+        return "no_applicable_tactic"
+    return "other"
+
+
+def infer_joint_label(domain: str, proof_method: str, obstruction: str | None, code: str) -> str:
+    """Infer the type of proof joint from domain + method."""
+    if obstruction:
+        return f"{domain}_{obstruction}"
+    if proof_method == "rfl":
+        return "definitional_equality"
+    if proof_method == "simp":
+        return "rewrite_normalization"
+    if proof_method == "omega":
+        return "linear_constraint_closure"
+    if proof_method == "nlinarith":
+        return "nonlinear_constraint_closure"
+    if proof_method == "induction":
+        return "inductive_step"
+    if proof_method == "rewrite":
+        return "equality_chaining"
+    if proof_method == "apply":
+        return "rule_application"
+    if proof_method == "cases":
+        return "case_analysis"
+    if proof_method == "intro":
+        return "implication_introduction"
+    if proof_method == "calc":
+        return "calculational_proof"
+    if proof_method == "incomplete":
+        return "incomplete_proof"
+    return f"{domain}_{proof_method}"
+
+
+def infer_rrc_shape(domain: str, proof_method: str, obstruction: str | None, code: str) -> str:
+    """Map domain + proof method to RRCShape."""
+    if obstruction:
+        return "HoldForUnlawfulOrUnderspecifiedShape"
+    if proof_method in ("rfl", "simp"):
+        return "CognitiveLoadField"
+    if proof_method in ("omega", "nlinarith", "ring"):
+        return "SignalShapedRouteCompiler"
+    if proof_method == "induction":
+        return "ProjectableGeometryTopology"
+    if proof_method in ("rewrite", "calc"):
+        return "CadForceProbeReceipt"
+    if proof_method in ("apply", "cases", "intro"):
+        return "LogogramProjection"
+    if domain in ("logic", "order"):
+        return "LogogramProjection"
+    return "HoldForUnlawfulOrUnderspecifiedShape"
+
+
+def find_in_receipt(receipts: list[dict], name: str) -> dict | None:
+    for r in receipts:
+        if r.get("theorem_name") == name:
+            return r
+    return None
+
+
+def main():
+    # Load results
+    results = []
+    with open(RESULTS_PATH) as f:
+        for line in f:
+            results.append(json.loads(line))
+
+    # Load receipts
+    receipts = []
+    with open(RECEIPTS_PATH) as f:
+        for line in f:
+            receipts.append(json.loads(line))
+
+    print(f"Loaded {len(results)} results, {len(receipts)} receipts", flush=True)
+
+    labels = []
+    for r in results:
+        name = r["name"]
+        receipt = find_in_receipt(receipts, name)
+        code = receipt.get("theorem_statement", "") if receipt else ""
+        stdout = receipt.get("theorem_statement", "") if receipt else ""
+
+        status = r["status"]
+        proof_method = parse_proof_method(code)
+        domain = parse_domain(code)
+        obstruction = parse_obstruction(
+            r.get("stdout", ""),
+            r.get("stderr", ""),
+            status,
+            code,
+        )
+        joint = infer_joint_label(domain, proof_method, obstruction, code)
+        rrc = infer_rrc_shape(domain, proof_method, obstruction, code)
+
+        labels.append({
+            "theorem_name": name,
+            "canonical_hash": r.get("canonical_hash", ""),
+            "proof_status": "verified" if "verified" in status else status,
+            "domain_label": domain,
+            "proof_method_label": proof_method,
+            "obstruction_label": obstruction,
+            "joint_label": joint,
+            "manual_rrc_shape": rrc,
+            "label_confidence": "auto",
+            "notes": code[:80] if code else "",
+        })
+
+    # Print distribution
+    print(f"\nGround-truth labels generated: {len(labels)}", flush=True)
+
+    print("\nDomain distribution:")
+    for label, count in sorted(Counter(l["domain_label"] for l in labels).items()):
+        print(f"  {label:20s}: {count:3d}")
+
+    print("\nProof method distribution:")
+    for label, count in sorted(Counter(l["proof_method_label"] for l in labels).items()):
+        print(f"  {label:20s}: {count:3d}")
+    print("\nObstruction distribution:")
+    for label, count in sorted(Counter(l["obstruction_label"] or "none" for l in labels).items()):
+        print(f"  {label:30s}: {count:3d}")
+
+    print("\nManual RRCShape distribution:")
+    for label, count in sorted(Counter(l["manual_rrc_shape"] for l in labels).items()):
+        print(f"  {label:35s}: {count:3d}")
+
+    # Write labels
+    with open(LABELS_PATH, "w") as f:
+        for lbl in labels:
+            f.write(json.dumps(lbl) + "\n")
+    print(f"\nLabels: {LABELS_PATH}", flush=True)
+
+    # Cross-tab: domain × proof method
+    print("\nCross-tab: domain × proof method")
+    domains = sorted(set(l["domain_label"] for l in labels))
+    methods = sorted(set(l["proof_method_label"] for l in labels))
+    header = "  " + "".join(f"{m:18s}" for m in methods)
+    print(header)
+    for dom in domains:
+        row = f"  {dom:12s}"
+        for meth in methods:
+            c = sum(1 for l in labels if l["domain_label"] == dom and l["proof_method_label"] == meth)
+            row += f"{c:>18d}"
+        print(row)
+
+    return 0
+
+
+if __name__ == "__main__":
+    main()
diff --git a/4-Infrastructure/shim/pist_enrich_and_train.py b/4-Infrastructure/shim/pist_enrich_and_train.py
new file mode 100644
index 00000000..5b354946
--- /dev/null
+++ b/4-Infrastructure/shim/pist_enrich_and_train.py
@@ -0,0 +1,301 @@
+#!/usr/bin/env python3
+"""Enrich canary receipts with full spectral data, then train classifiers.
+
+Usage:
+    python3 pist_enrich_and_train.py
+"""
+
+import json
+import os
+import subprocess
+import sys
+import tempfile
+from collections import defaultdict
+from math import sqrt
+
+PIST_DECOMPOSE = os.environ.get(
+    "PIST_DECOMPOSE_BIN",
+    "/home/allaun/.local/share/opencode/worktree/"
+    "0b42981cf7f7d5e172b1e93f8d4bb64a3dd63962/Turn-and-Burn/infra/rust/"
+    "ene-rds/target/release/pist-decompose",
+)
+
+FEATURE_NAMES = [
+    "zero_mode_proxy_count",
+    "rank_estimate",
+    "laplacian_zero_count",
+    "spectral_gap",
+    "crossing_density",
+    "strand_entropy",
+]
+EIGEN_LEN = 8
+SINGULAR_LEN = 8
+
+
+def extract(pist_out: dict) -> dict:
+    """Extract flattened feature vector from pist-decompose output."""
+    spectral = pist_out.get("spectral", {})
+    braid = pist_out.get("braid", {})
+    gamma = pist_out.get("gamma_packet", {})
+    zmp = spectral.get("zero_mode_proxy_count", 0)
+    rank = spectral.get("rank_estimate", 0)
+    lap0 = spectral.get("laplacian_zero_count", 0)
+    gap = spectral.get("symmetric_spectral_gap", 0)
+    cd = braid.get("crossing_density", 0)
+    sent = braid.get("strand_entropy", 0)
+    ev = spectral.get("symmetric_eigenvalues")
+    if not ev:
+        ev = [0.0] * 8
+    sv = spectral.get("singular_values")
+    if not sv:
+        sv = [0.0] * 8
+    slack = braid.get("sidon_slack", 0)
+    yb = braid.get("yang_baxter_valid", True)
+    steps = braid.get("step_count", 0)
+    gamma_v = gamma.get("gamma", {}).get("value", 0)
+    chi = gamma.get("chi", 0)
+    kappa = gamma.get("kappa", 0)
+    tau = gamma.get("tau", 0)
+    theta = gamma.get("theta", 0)
+    eps = gamma.get("epsilon", 0)
+    mhash = braid.get("matrix_hash", "?")[:16]
+    chash = pist_out.get("canonical_hash", "?")[:16]
+
+    vec = [zmp, rank, lap0, gap, cd, sent]
+    for v in ev[:EIGEN_LEN]:
+        vec.append(float(v))
+    for v in sv[:SINGULAR_LEN]:
+        vec.append(float(v))
+    vec.extend([float(slack), 1.0 if yb else 0.0, float(steps),
+                float(gamma_v), float(chi), float(kappa), float(tau), float(theta), float(eps)])
+
+    return {
+        "vector": vec,
+        "features": dict(zip(FEATURE_NAMES, [zmp, rank, lap0, gap, cd, sent])),
+        "eigenvalues": [round(float(v), 6) for v in ev[:EIGEN_LEN]],
+        "singular_values": [round(float(v), 6) for v in sv[:SINGULAR_LEN]],
+        "matrix_hash": mhash,
+        "canonical_hash": chash,
+    }
+
+
+def run_pist(receipt: dict) -> dict | None:
+    """Run pist-decompose on a receipt dict."""
+    with tempfile.NamedTemporaryFile(mode="w", suffix=".json", delete=False) as f:
+        json.dump(receipt, f)
+        fpath = f.name
+    try:
+        r = subprocess.run([PIST_DECOMPOSE, fpath, "--num-leaves", "8"],
+                           capture_output=True, text=True, timeout=30)
+        if r.returncode != 0:
+            return None
+        return json.loads(r.stdout)
+    finally:
+        os.unlink(fpath)
+
+
+def normalize(vectors):
+    n = len(vectors)
+    if n == 0:
+        return vectors, [], []
+    dim = len(vectors[0])
+    means = [sum(v[i] for v in vectors) / n for i in range(dim)]
+    stds = [sqrt(sum((v[i] - means[i]) ** 2 for v in vectors) / max(n - 1, 1)) for i in range(dim)]
+    stds = [s if s > 1e-9 else 1.0 for s in stds]
+    return [[(v[i] - means[i]) / stds[i] for i in range(dim)] for v in vectors], means, stds
+
+
+def centroid(vecs):
+    if not vecs:
+        return []
+    return [sum(v[i] for v in vecs) / len(vecs) for i in range(len(vecs[0]))]
+
+
+def euclidean(a, b):
+    return sqrt(sum((a[i] - b[i]) ** 2 for i in range(len(a))))
+
+
+def knn(train_v, train_l, test_v, k):
+    dists = [(euclidean(test_v, tv), tl) for tv, tl in zip(train_v, train_l)]
+    dists.sort(key=lambda x: x[0])
+    nearest = dists[:k]
+    votes = defaultdict(int)
+    for _, lbl in nearest:
+        votes[lbl] += 1
+    return max(votes, key=votes.get)
+
+
+def eval_loocv(vectors, labels, method="centroid", k=3):
+    n = len(vectors)
+    correct = 0
+    top2_correct = 0
+    confusion = defaultdict(lambda: defaultdict(int))
+    for i in range(n):
+        train_v = vectors[:i] + vectors[i + 1:]
+        train_l = labels[:i] + labels[i + 1:]
+        test_v = vectors[i]
+        test_l = labels[i]
+
+        if method == "centroid":
+            cls_vecs = defaultdict(list)
+            for v, lbl in zip(train_v, train_l):
+                cls_vecs[lbl].append(v)
+            centroids_dict = {lbl: centroid(vecs) for lbl, vecs in cls_vecs.items()}
+            pred = min(centroids_dict, key=lambda lbl: euclidean(test_v, centroids_dict[lbl]))
+        else:
+            pred = knn(train_v, train_l, test_v, k)
+
+        confusion[test_l][pred] += 1
+        if pred == test_l:
+            correct += 1
+        # top-2 check
+        dists = sorted([(euclidean(test_v, cent), lbl) for lbl, cent in centroids_dict.items()]) if method == "centroid" else sorted([(euclidean(test_v, train_v[j]), train_l[j]) for j in range(len(train_v))])[:2]
+        top2_labels = [d[1] for d in dists[:2]]
+        if test_l in top2_labels:
+            top2_correct += 1
+
+    return correct / n, top2_correct / n, dict(confusion)
+
+
+def main():
+    receipts_path = os.path.join(os.path.dirname(__file__), "../..",
+                                 "shared-data/pist_canary_receipts.jsonl")
+    results_path = os.path.join(os.path.dirname(__file__), "../..",
+                                "shared-data/pist_canary_results.jsonl")
+
+    # Load results to get status/label info
+    results = []
+    with open(results_path) as f:
+        for line in f:
+            row = json.loads(line)
+            results.append(row)
+
+    # Load receipts and run PIST
+    enriched = []
+    with open(receipts_path) as f:
+        for line in f:
+            receipt = json.loads(line)
+            enriched.append(receipt)
+
+    print(f"Loaded {len(results)} results, {len(enriched)} receipts", flush=True)
+
+    # Re-run PIST on each receipt and collect rich features
+    records = []
+    for i, (result, receipt) in enumerate(zip(results, enriched)):
+        print(f"  [{i+1}/{len(results)}] {result['name']:30s} ... ", end="", flush=True)
+        pist = run_pist(receipt)
+        if pist is None:
+            print("PIST FAILED", flush=True)
+            continue
+        ext = extract(pist)
+        records.append({
+            "name": result["name"],
+            "status": result["status"],
+            "ok": result["ok"],
+            "rrc_shape": result["exact_shape"],
+            "vector": ext["vector"],
+            "features": ext["features"],
+            "eigenvalues": ext["eigenvalues"],
+            "singular_values": ext["singular_values"],
+            "matrix_hash": ext["matrix_hash"],
+            "canonical_hash": ext["canonical_hash"],
+        })
+        print(f"ok ZMP={ext['features']['zero_mode_proxy_count']} rank={ext['features']['rank_estimate']}", flush=True)
+
+    n = len(records)
+    print(f"\nEnriched: {n} records", flush=True)
+
+    # Prepare feature matrix
+    vectors = [r["vector"] for r in records]
+    dim = len(vectors[0])
+    normed, means, stds = normalize(vectors)
+
+    # Feature variance
+    print(f"\nFeature dimensions: {dim}", flush=True)
+    for i, name in enumerate(FEATURE_NAMES):
+        vals = [r["features"][name] for r in records]
+        uniq = len(set(vals))
+        print(f"  {name:25s} uniq={uniq:2d} vals=[{min(vals)},{max(vals)}]", flush=True)
+
+    # ── Train on rrc_shape ──
+    print("\n" + "=" * 60, flush=True)
+    print("TARGET: RRCShape (exact classifier prediction)", flush=True)
+    print("=" * 60, flush=True)
+
+    labels_rrc = [r["rrc_shape"] for r in records]
+    unique_labels = sorted(set(labels_rrc))
+    print(f"Labels: {unique_labels}", flush=True)
+    print(f"Distribution: {dict(Counter(labels_rrc))}", flush=True)
+
+    for method, k in [("centroid", None), ("knn_1", 1), ("knn_3", 3), ("knn_5", 5)]:
+        acc, top2, conf = eval_loocv(normed, labels_rrc, "centroid" if method == "centroid" else "knn", k or 3)
+        print(f"\n  {method:15s} LOOCV accuracy: {acc:.1%} ({int(acc*n)}/{n}) top-2: {top2:.1%}", flush=True)
+
+    # ── Train on proof status ──
+    print("\n" + "=" * 60, flush=True)
+    print("TARGET: Proof Status (verified vs failed)", flush=True)
+    print("=" * 60, flush=True)
+
+    # Map status to binary
+    def status_binary(s):
+        return "verified" if "verified" in s else "failed"
+
+    labels_status = [status_binary(r["status"]) for r in records]
+    print(f"Distribution: {dict(Counter(labels_status))}", flush=True)
+
+    for method, k in [("centroid", None), ("knn_3", 3), ("knn_5", 5)]:
+        acc, top2, conf = eval_loocv(normed, labels_status, "centroid" if method == "centroid" else "knn", k or 3)
+        print(f"  {method:15s} LOOCV accuracy: {acc:.1%} ({int(acc*n)}/{n})", flush=True)
+
+    # ── Save enriched features ──
+    vec_path = os.path.join(os.path.dirname(__file__), "../..",
+                            "shared-data/pist_canary_feature_vectors.jsonl")
+    with open(vec_path, "w") as f:
+        for r in records:
+            f.write(json.dumps({
+                "name": r["name"],
+                "status": r["status"],
+                "rrc_shape": r["rrc_shape"],
+                "vector": [round(x, 6) for x in r["vector"]],
+                "features": r["features"],
+                "eigenvalues": r["eigenvalues"],
+                "singular_values": r["singular_values"],
+            }) + "\n")
+    print(f"\nFeature vectors: {vec_path}", flush=True)
+
+    # Summary report
+    rrc_acc_centroid, rrc_top2, _ = eval_loocv(normed, labels_rrc, "centroid")
+    status_acc_centroid, status_top2, _ = eval_loocv(normed, labels_status, "centroid", 3)
+
+    report = {
+        "n_samples": n,
+        "dimension": dim,
+        "targets": ["rrc_shape", "proof_status"],
+        "rrc_shape": {
+            "unique_labels": unique_labels,
+            "distribution": dict(Counter(labels_rrc)),
+            "centroid_loocv_accuracy": round(rrc_acc_centroid, 4),
+            "centroid_top2_accuracy": round(rrc_top2, 4),
+        },
+        "proof_status": {
+            "distribution": dict(Counter(labels_status)),
+            "centroid_loocv_accuracy": round(status_acc_centroid, 4),
+            "centroid_top2_accuracy": round(status_top2, 4),
+        },
+        "warnings": [
+            f"Only {n} samples; all results are calibration-only.",
+            "Do not promote to production until 278+ labeled artifacts.",
+        ],
+    }
+
+    report_path = os.path.join(os.path.dirname(__file__), "../..",
+                               "shared-data/pist_canary_training_report.json")
+    with open(report_path, "w") as f:
+        json.dump(report, f, indent=2)
+    print(f"Report: {report_path}", flush=True)
+    return 0
+
+
+if __name__ == "__main__":
+    from collections import Counter
+    sys.exit(main())
diff --git a/4-Infrastructure/shim/pist_train_ground_truth.py b/4-Infrastructure/shim/pist_train_ground_truth.py
new file mode 100644
index 00000000..1e8277c2
--- /dev/null
+++ b/4-Infrastructure/shim/pist_train_ground_truth.py
@@ -0,0 +1,224 @@
+#!/usr/bin/env python3
+"""Train on ground-truth labels (proof_method, domain, RRCShape).
+
+Usage:
+    python3 pist_train_ground_truth.py
+"""
+
+import json
+import os
+import sys
+from collections import Counter, defaultdict
+from math import sqrt
+
+FEATURE_VECTORS_PATH = os.path.join(os.path.dirname(__file__), "../..",
+                                    "shared-data/pist_canary_feature_vectors.jsonl")
+GROUND_TRUTH_PATH = os.path.join(os.path.dirname(__file__), "../..",
+                                 "shared-data/pist_canary_ground_truth_labels.jsonl")
+REPORT_PATH = os.path.join(os.path.dirname(__file__), "../..",
+                           "shared-data/pist_canary_ground_truth_report.json")
+
+
+def load_vectors(path: str) -> list[dict]:
+    rows = []
+    with open(path) as f:
+        for line in f:
+            rows.append(json.loads(line))
+    return rows
+
+
+def load_labels(path: str) -> dict[str, dict]:
+    labels = {}
+    with open(path) as f:
+        for line in f:
+            row = json.loads(line)
+            labels[row["theorem_name"]] = row
+    return labels
+
+
+def normalize(vectors):
+    n = len(vectors)
+    if n == 0:
+        return vectors, [], []
+    dim = len(vectors[0])
+    means = [sum(v[i] for v in vectors) / n for i in range(dim)]
+    stds = [sqrt(sum((v[i] - means[i]) ** 2 for v in vectors) / max(n - 1, 1)) for i in range(dim)]
+    stds = [s if s > 1e-9 else 1.0 for s in stds]
+    return [[(v[i] - means[i]) / stds[i] for i in range(dim)] for v in vectors], means, stds
+
+
+def centroid(vecs):
+    if not vecs:
+        return []
+    return [sum(v[i] for v in vecs) / len(vecs) for i in range(len(vecs[0]))]
+
+
+def euclidean(a, b):
+    return sqrt(sum((a[i] - b[i]) ** 2 for i in range(len(a))))
+
+
+def knn(train_v, train_l, test_v, k):
+    dists = sorted(
+        [(euclidean(test_v, tv), tl) for tv, tl in zip(train_v, train_l)],
+        key=lambda x: x[0],
+    )
+    votes = Counter(tl for _, tl in dists[:k])
+    return votes.most_common(1)[0][0]
+
+
+def eval_loocv(vectors, labels, method="centroid", k=3):
+    n = len(vectors)
+    correct = 0
+    top2_correct = 0
+    n_classes = len(set(labels))
+    confusion = defaultdict(lambda: defaultdict(int))
+
+    for i in range(n):
+        tv = vectors[:i] + vectors[i + 1 :]
+        tl = labels[:i] + labels[i + 1 :]
+        test_v = vectors[i]
+        test_l = labels[i]
+
+        if method == "centroid":
+            cls_vecs = defaultdict(list)
+            for v, l in zip(tv, tl):
+                cls_vecs[l].append(v)
+            cents = {l: centroid(vecs) for l, vecs in cls_vecs.items()}
+            preds = sorted(
+                cents.keys(), key=lambda l: euclidean(test_v, cents[l])
+            )
+            pred = preds[0]
+            top2 = preds[:2]
+        else:
+            dists = sorted(
+                [(euclidean(test_v, tv[j]), tl[j]) for j in range(len(tv))],
+                key=lambda x: x[0],
+            )
+            votes = Counter(tl for _, tl in dists[:k])
+            pred = votes.most_common(1)[0][0]
+            top2 = [t[1] for t in dists[:2]]
+
+        confusion[test_l][pred] += 1
+        if pred == test_l:
+            correct += 1
+        if test_l in top2:
+            top2_correct += 1
+
+    baseline = max(Counter(labels).values()) / n
+    accuracy = correct / n
+    top2_acc = top2_correct / n
+
+    return accuracy, top2_acc, dict(confusion), baseline
+
+
+def main():
+    vectors = load_vectors(FEATURE_VECTORS_PATH)
+    gt_labels = load_labels(GROUND_TRUTH_PATH)
+
+    # Align by theorem name
+    aligned = []
+    for v in vectors:
+        name = v.get("name", "")
+        if name in gt_labels:
+            aligned.append((v["vector"], gt_labels[name]))
+
+    n = len(aligned)
+    print(f"Aligned vectors: {n}", flush=True)
+
+    raw_vecs = [a[0] for a in aligned]
+    normed, _, _ = normalize(raw_vecs)
+
+    TARGETS = [
+        ("proof_method_label", "Proof method"),
+        ("domain_label", "Domain"),
+        ("manual_rrc_shape", "Manual RRCShape"),
+    ]
+
+    results = {}
+
+    for key, label in TARGETS:
+        print(f"\n{'=' * 60}", flush=True)
+        print(f"TARGET: {label} ({key})", flush=True)
+        print(f"{'=' * 60}", flush=True)
+
+        targets = []
+        for _, lbl in aligned:
+            val = lbl.get(key, "unknown")
+            if val is None:
+                val = "none"
+            targets.append(val)
+
+        classes = sorted(set(targets))
+        dist = Counter(targets)
+        print(f"Classes ({len(classes)}): {dict(dist)}", flush=True)
+        print(f"Majority baseline: {max(dist.values()) / n:.1%}", flush=True)
+
+        for method, k in [
+            ("centroid", None),
+            ("knn_1", 1),
+            ("knn_3", 3),
+            ("knn_5", 5),
+        ]:
+            acc, top2, conf, baseline = eval_loocv(
+                normed, targets, "centroid" if method == "centroid" else "knn", k or 3
+            )
+            print(f"  {method:15s} LOOCV: {acc:.1%} ({int(acc * n)}/{n})  "
+                  f"top-2: {top2:.1%}  baseline: {baseline:.1%}",
+                  flush=True)
+
+            results[f"{key}_{method}"] = {
+                "accuracy": round(acc, 4),
+                "top2_accuracy": round(top2, 4),
+                "baseline": round(baseline, 4),
+                "n_classes": len(classes),
+                "improvement_over_baseline": round(
+                    (acc - baseline) / max(baseline, 0.01), 3
+                ),
+            }
+
+        # Per-class accuracy for centroid
+        acc, top2, conf, baseline = eval_loocv(normed, targets, "centroid", 3)
+        print(f"\n  Per-class (centroid):")
+        for cls in sorted(classes):
+            tp = conf.get(cls, {}).get(cls, 0)
+            total = dist[cls]
+            print(f"    {cls:30s} {tp:3d}/{total:3d} = {tp / max(total, 1):.0%}")
+
+    # Summary
+    print(f"\n{'=' * 60}", flush=True)
+    print("SUMMARY", flush=True)
+    print(f"{'=' * 60}", flush=True)
+    print(f"{'Target':35s} {'Method':12s} {'Acc':>6} {'Top-2':>6} {'Base':>6} {'Impr':>6}",
+          flush=True)
+    print(f"{'-' * 35} {'-' * 12} {'-' * 6} {'-' * 6} {'-' * 6} {'-' * 6}", flush=True)
+    for key, label in TARGETS:
+        for method in ["centroid", "knn_3"]:
+            r = results.get(f"{key}_{method}", {})
+            impr = r.get("improvement_over_baseline", 0)
+            print(
+                f"{label:35s} {method:12s} "
+                f"{r.get('accuracy', 0):6.1%} "
+                f"{r.get('top2_accuracy', 0):6.1%} "
+                f"{r.get('baseline', 0):6.1%} "
+                f"{impr:>+5.1f}x",
+                flush=True,
+            )
+
+    # Save report
+    report = {
+        "n_samples": n,
+        "targets": {key: label for key, label in TARGETS},
+        "results": results,
+        "class_distributions": {
+            key: dict(Counter([lbl.get(key, "unknown") for _, lbl in aligned]))
+            for key, _ in TARGETS
+        },
+    }
+    with open(REPORT_PATH, "w") as f:
+        json.dump(report, f, indent=2)
+    print(f"\nReport: {REPORT_PATH}", flush=True)
+    return 0
+
+
+if __name__ == "__main__":
+    main()
diff --git a/shared-data/pist_flexure_library_report.json b/shared-data/pist_flexure_library_report.json
index e953f0e0..83ae5be4 100644
--- a/shared-data/pist_flexure_library_report.json
+++ b/shared-data/pist_flexure_library_report.json
@@ -1,5 +1,5 @@
 {
-  "session_id": "d94c6353-5ed9-42a4-b2b7-d0fee8b36a8e",
+  "session_id": "ae31d595-0535-4a0c-9d41-af9c0357dba1",
   "total_flexures_inserted": 38,
   "total_skipped": 0,
   "motifs_discovered": 10,