diff --git a/cmd/semantic-router-experiment/experiment.py b/cmd/semantic-router-experiment/experiment.py
new file mode 100644
index 0000000..266315e
--- /dev/null
+++ b/cmd/semantic-router-experiment/experiment.py
@@ -0,0 +1,856 @@
+#!/usr/bin/env python3
+"""
+Semantic Router Linear Head Experiment
+======================================
+
+Evaluates whether the six-way residual routing boundary is linearly learnable
+from Maven's existing e5-small representation.
+
+Architecture under test:
+    NormalizedInput.MatchText
+    → existing multilingual-e5-small embedder
+    → tiny 6-class linear head (softmax logistic regression)
+    → SemanticRouteDecision
+
+No new encoder. No fine-tuning. Linear separability floor only.
+"""
+
+import json
+import hashlib
+import sys
+import warnings
+from collections import defaultdict
+from pathlib import Path
+
+import numpy as np
+from sklearn.exceptions import ConvergenceWarning
+from sklearn.linear_model import LogisticRegression
+from sklearn.metrics import (
+    accuracy_score,
+    f1_score,
+    precision_recall_fscore_support,
+    confusion_matrix,
+    brier_score_loss,
+    log_loss,
+)
+from sklearn.preprocessing import LabelEncoder
+
+warnings.filterwarnings("ignore", category=ConvergenceWarning)
+
+# ─── Constants ──────────────────────────────────────────────────────────────
+
+ROUTES = ["action", "conversation", "knowledge", "memory_write", "system", "uncertain"]
+ROUTE_IDX = {r: i for i, r in enumerate(ROUTES)}
+
+# Regularization grid
+C_VALUES = [0.01, 0.1, 1.0, 10.0, 100.0]
+
+# Abstention thresholds
+THRESHOLDS = [0.40, 0.50, 0.60, 0.70, 0.80, 0.90]
+
+# Contrast families
+CONTRAST_FAMILIES = [
+    "negation", "question", "reported_speech", "quotation",
+    "hypothetical", "capability_question",
+]
+
+# ─── Data Loading ───────────────────────────────────────────────────────────
+
+def load_embeddings(path):
+    with open(path) as f:
+        data = json.load(f)
+    meta = data["meta"]
+    examples = data["examples"]
+    return meta, examples
+
+
+def filter_dev_pool(examples):
+    """Return only development pool examples (not frozen holdout)."""
+    return [e for e in examples if e["dev_pool"]]
+
+
+def filter_residual(examples):
+    """Return only router-residual examples (fast_path_resolved == false)."""
+    return [e for e in examples if not e["fast_path_resolved"]]
+
+
+def extract_Xy(examples):
+    """Extract feature matrix and label vector from cached examples."""
+    X = np.array([e["embedding"] for e in examples])
+    y = np.array([e["route"] for e in examples])
+    return X, y
+
+
+def get_fold_groups(examples):
+    """Return fold assignment array matching the Go-generated CV folds."""
+    return np.array([e["cv_fold"] for e in examples])
+
+
+# ─── Grouped CV ─────────────────────────────────────────────────────────────
+
+def grouped_cv_experiment(X, y, fold_ids, C_values, examples_meta):
+    """
+    Run grouped cross-validation with logistic regression.
+    Returns best C, per-fold results, and out-of-fold predictions.
+    """
+    unique_folds = sorted(set(fold_ids))
+    n_classes = len(ROUTES)
+
+    results_by_C = {}
+    for C in C_values:
+        fold_metrics = []
+        oof_rows = []
+
+        for test_fold in unique_folds:
+            train_mask = fold_ids != test_fold
+            test_mask = fold_ids == test_fold
+
+            X_train, y_train = X[train_mask], y[train_mask]
+            X_test, y_test = X[test_mask], y[test_mask]
+
+            model = LogisticRegression(
+                C=C, max_iter=2000, solver="lbfgs",
+                random_state=42,
+            )
+            model.fit(X_train, y_train)
+
+            y_pred = model.predict(X_test)
+            y_proba = model.predict_proba(X_test)
+            classes = model.classes_
+
+            acc = accuracy_score(y_test, y_pred)
+            macro_f1 = f1_score(y_test, y_pred, average="macro", zero_division=0)
+
+            prec, rec, f1, sup = precision_recall_fscore_support(
+                y_test, y_pred, labels=ROUTES, zero_division=0
+            )
+
+            false_action = 0
+            for true, pred in zip(y_test, y_pred):
+                if true != "action" and pred == "action":
+                    false_action += 1
+
+            fold_metrics.append({
+                "fold": test_fold,
+                "train_size": len(X_train),
+                "test_size": len(X_test),
+                "accuracy": acc,
+                "macro_f1": macro_f1,
+                "per_class": {
+                    route: {"precision": prec[i], "recall": rec[i], "f1": f1[i]}
+                    for i, route in enumerate(ROUTES)
+                },
+                "false_action": false_action,
+            })
+
+            for i, (true, pred) in enumerate(zip(y_test, y_pred)):
+                test_indices = np.where(test_mask)[0]
+                meta = examples_meta[test_indices[i]]
+                proba_dict = {cls: float(y_proba[i][j]) for j, cls in enumerate(classes)}
+                max_proba = max(proba_dict.values()) if proba_dict else 0.0
+                oof_rows.append({
+                    "source_id": meta["source_id"],
+                    "fold": test_fold,
+                    "true": true,
+                    "predicted": pred,
+                    "correct": true == pred,
+                    "max_proba": max_proba,
+                    "proba": proba_dict,
+                    "fast_path_resolved": meta.get("fast_path_resolved", False),
+                    "tags": meta.get("tags", []),
+                    "text": meta["text"],
+                })
+
+        # Aggregate across folds
+        mean_acc = np.mean([m["accuracy"] for m in fold_metrics])
+        mean_f1 = np.mean([m["macro_f1"] for m in fold_metrics])
+        std_acc = np.std([m["accuracy"] for m in fold_metrics])
+        std_f1 = np.std([m["macro_f1"] for m in fold_metrics])
+        total_fa = sum(m["false_action"] for m in fold_metrics)
+
+        results_by_C[C] = {
+            "mean_accuracy": mean_acc,
+            "std_accuracy": std_acc,
+            "mean_macro_f1": mean_f1,
+            "std_macro_f1": std_f1,
+            "total_false_action": total_fa,
+            "fold_metrics": fold_metrics,
+            "oof_predictions": oof_rows,
+        }
+
+    # Select best C by mean macro F1
+    best_C = max(results_by_C, key=lambda c: results_by_C[c]["mean_macro_f1"])
+    return best_C, results_by_C
+
+
+# ─── Metrics Computation ────────────────────────────────────────────────────
+
+def compute_full_metrics(y_true, y_pred, y_proba=None):
+    """Compute all required metrics from out-of-fold predictions."""
+    acc = accuracy_score(y_true, y_pred)
+    macro_f1 = f1_score(y_true, y_pred, average="macro", zero_division=0)
+
+    # Per-class P/R/F1
+    prec, rec, f1, sup = precision_recall_fscore_support(
+        y_true, y_pred, labels=ROUTES, zero_division=0
+    )
+
+    # Confusion matrix
+    cm = confusion_matrix(y_true, y_pred, labels=ROUTES)
+
+    # False action
+    false_action = 0
+    false_action_predicted = []
+    for true, pred in zip(y_true, y_pred):
+        if true != "action" and pred == "action":
+            false_action += 1
+
+    # Action precision/recall
+    action_tp = sum(1 for t, p in zip(y_true, y_pred) if t == "action" and p == "action")
+    action_fp = sum(1 for t, p in zip(y_true, y_pred) if t != "action" and p == "action")
+    action_fn = sum(1 for t, p in zip(y_true, y_pred) if t == "action" and p != "action")
+    action_precision = action_tp / max(action_tp + action_fp, 1)
+    action_recall = action_tp / max(action_tp + action_fn, 1)
+
+    # Uncertain precision/recall
+    unc_tp = sum(1 for t, p in zip(y_true, y_pred) if t == "uncertain" and p == "uncertain")
+    unc_fp = sum(1 for t, p in zip(y_true, y_pred) if t != "uncertain" and p == "uncertain")
+    unc_fn = sum(1 for t, p in zip(y_true, y_pred) if t == "uncertain" and p != "uncertain")
+    unc_precision = unc_tp / max(unc_tp + unc_fp, 1)
+    unc_recall = unc_tp / max(unc_tp + unc_fn, 1)
+
+    metrics = {
+        "accuracy": acc,
+        "macro_f1": macro_f1,
+        "false_action_count": false_action,
+        "false_action_rate": false_action / max(len(y_true), 1),
+        "action_precision": action_precision,
+        "action_recall": action_recall,
+        "uncertain_precision": unc_precision,
+        "uncertain_recall": unc_recall,
+        "per_class": {},
+        "confusion_matrix": cm.tolist(),
+    }
+
+    for i, route in enumerate(ROUTES):
+        metrics["per_class"][route] = {
+            "precision": float(prec[i]),
+            "recall": float(rec[i]),
+            "f1": float(f1[i]),
+            "support": int(sup[i]),
+        }
+
+    return metrics
+
+
+def compute_calibration(y_true, y_pred, oof_rows):
+    """Compute ECE, Brier score, and per-threshold abstention curves."""
+    # ECE (Expected Calibration Error) with 10 bins
+    n_bins = 10
+    bin_boundaries = np.linspace(0, 1, n_bins + 1)
+    ece = 0.0
+    total = len(y_true)
+
+    confidences = np.array([r["max_proba"] for r in oof_rows])
+    correct = np.array([r["correct"] for r in oof_rows])
+
+    for i in range(n_bins):
+        lo, hi = bin_boundaries[i], bin_boundaries[i + 1]
+        mask = (confidences > lo) & (confidences <= hi)
+        if mask.sum() == 0:
+            continue
+        bin_acc = correct[mask].mean()
+        bin_conf = confidences[mask].mean()
+        ece += mask.sum() / total * abs(bin_acc - bin_conf)
+
+    # Brier score (multiclass one-hot encoding)
+    y_true_idx = np.array([ROUTES.index(t) for t in y_true])
+    n_samples = len(y_true)
+    n_classes = len(ROUTES)
+    y_true_oh = np.zeros((n_samples, n_classes))
+    y_true_oh[np.arange(n_samples), y_true_idx] = 1.0
+
+    # Build probability matrix
+    y_proba_matrix = np.zeros((n_samples, n_classes))
+    for i, r in enumerate(oof_rows):
+        for j, route in enumerate(ROUTES):
+            y_proba_matrix[i, j] = r["proba"].get(route, 0.0)
+
+    brier = np.mean(np.sum((y_proba_matrix - y_true_oh) ** 2, axis=1))
+    ll = log_loss(y_true, y_proba_matrix, labels=ROUTES)
+
+    # Abstention curves
+    abstention_curves = []
+    for thr in THRESHOLDS:
+        accepted_mask = confidences >= thr
+        n_accepted = accepted_mask.sum()
+        coverage = n_accepted / max(total, 1)
+
+        if n_accepted > 0:
+            acc_accepted = accuracy_score(y_true[accepted_mask], y_pred[accepted_mask])
+            f1_accepted = f1_score(
+                y_true[accepted_mask], y_pred[accepted_mask],
+                average="macro", zero_division=0,
+            )
+            fa_count = sum(
+                1 for t, p in zip(y_true[accepted_mask], y_pred[accepted_mask])
+                if t != "action" and p == "action"
+            )
+        else:
+            acc_accepted = 0.0
+            f1_accepted = 0.0
+            fa_count = 0
+
+        abstention_curves.append({
+            "threshold": thr,
+            "n_accepted": int(n_accepted),
+            "coverage": coverage,
+            "accuracy": acc_accepted,
+            "macro_f1": f1_accepted,
+            "false_action_count": fa_count,
+        })
+
+    return {
+        "ece": float(ece),
+        "brier": float(brier),
+        "log_loss": float(ll),
+        "abstention_curves": abstention_curves,
+    }
+
+
+def compute_action_threshold(oof_rows):
+    """Evaluate action-specific safety gate thresholds."""
+    action_curves = []
+    for thr in THRESHOLDS:
+        action_pred = []
+        for r in oof_rows:
+            p = r["predicted"]
+            proba = r["proba"].get("action", 0.0)
+            if p == "action" and proba < thr:
+                # Demote action prediction
+                # Find next best route that isn't action
+                sorted_routes = sorted(r["proba"].items(), key=lambda x: -x[1])
+                for route, _ in sorted_routes:
+                    if route != "action":
+                        p = route
+                        break
+            action_pred.append(p)
+
+        y_true = np.array([r["true"] for r in oof_rows])
+        y_pred = np.array(action_pred)
+
+        action_tp = sum(1 for t, p in zip(y_true, y_pred) if t == "action" and p == "action")
+        action_fp = sum(1 for t, p in zip(y_true, y_pred) if t != "action" and p == "action")
+        action_fn = sum(1 for t, p in zip(y_true, y_pred) if t == "action" and p != "action")
+        false_action = sum(1 for t, p in zip(y_true, y_pred) if t != "action" and p == "action")
+
+        action_curves.append({
+            "threshold": thr,
+            "action_precision": action_tp / max(action_tp + action_fp, 1),
+            "action_recall": action_tp / max(action_tp + action_fn, 1),
+            "false_action_count": false_action,
+        })
+
+    return action_curves
+
+
+def compute_disagreement(y_true, y_pred, legacy_pred, oof_rows):
+    """Analyze disagreements between legacy and learned router."""
+    results = {
+        "legacy_wrong_learned_right": [],
+        "legacy_right_learned_wrong": [],
+        "both_wrong_differently": [],
+    }
+
+    legacy_false_actions = []
+    learned_false_actions = []
+    shared_false_actions = []
+
+    for i, r in enumerate(oof_rows):
+        sid = r["source_id"]
+        text = r["text"]
+        true = y_true[i]
+        learned = y_pred[i]
+        legacy = legacy_pred[i]
+
+        if legacy != true and learned == true:
+            results["legacy_wrong_learned_right"].append({
+                "source_id": sid, "text": text,
+                "true": true, "legacy": legacy, "learned": learned,
+            })
+        elif legacy == true and learned != true:
+            results["legacy_right_learned_wrong"].append({
+                "source_id": sid, "text": text,
+                "true": true, "legacy": legacy, "learned": learned,
+            })
+        elif legacy != true and learned != true and legacy != learned:
+            results["both_wrong_differently"].append({
+                "source_id": sid, "text": text,
+                "true": true, "legacy": legacy, "learned": learned,
+            })
+
+        # False action tracking
+        if true != "action" and legacy == "action":
+            legacy_false_actions.append(sid)
+        if true != "action" and learned == "action":
+            learned_false_actions.append(sid)
+        if true != "action" and legacy == "action" and learned == "action":
+            shared_false_actions.append(sid)
+
+    # Repaired false actions
+    repaired = [sid for sid in legacy_false_actions if sid not in learned_false_actions]
+    new_errors = [sid for sid in learned_false_actions if sid not in legacy_false_actions]
+    shared = shared_false_actions
+
+    return {
+        "details": results,
+        "legacy_false_actions": legacy_false_actions,
+        "learned_false_actions": learned_false_actions,
+        "repaired": repaired,
+        "new_errors": new_errors,
+        "shared": shared,
+    }
+
+
+def compute_contrast_family(oof_rows):
+    """Analyze performance per contrast family."""
+    family_results = {}
+    for family in CONTRAST_FAMILIES:
+        members = [r for r in oof_rows if family in r.get("tags", [])]
+        if not members:
+            continue
+        y_true = [r["true"] for r in members]
+        y_pred = [r["predicted"] for r in members]
+        correct = sum(1 for t, p in zip(y_true, y_pred) if t == p)
+        false_act = sum(1 for t, p in zip(y_true, y_pred) if t != "action" and p == "action")
+        family_results[family] = {
+            "count": len(members),
+            "correct": correct,
+            "accuracy": correct / len(members),
+            "false_action": false_act,
+        }
+    return family_results
+
+
+def compute_legacy_baseline(examples):
+    """
+    Compute legacy baseline by mapping each example through the known
+    fast-path and classifier behavior. Since we don't have the actual
+    router running, we use the corpus metadata:
+    - fast_path_resolved examples are correct (grammar handles them)
+    - We simulate the legacy baseline from the test output numbers
+    """
+    # The actual legacy baseline was measured in TestLegacyBaseline:
+    # overall accuracy: 52.2%, residual: 40.8%, false-action rate: 19.9%
+    # We need per-example predictions. We'll approximate from the corpus structure.
+    #
+    # For fast_path_resolved examples, the legacy router is correct (stage-0 grammar).
+    # For residual examples, we need to simulate the hash-embedder classifier.
+    # Since we don't have the hash embedder running, we use the known aggregate.
+    #
+    # This is a limitation: the legacy baseline numbers come from the Go test,
+    # and we only have the aggregate. We'll use the aggregate for comparison.
+    pass
+
+
+# ─── Report Generation ──────────────────────────────────────────────────────
+
+def fmt_pct(v, decimals=1):
+    return f"{100 * v:.{decimals}f}%"
+
+
+def fmt_float(v, decimals=3):
+    return f"{v:.{decimals}f}"
+
+
+def generate_report(meta, all_results, residual_results, oof_all, oof_residual):
+    """Generate the full experiment report."""
+    lines = []
+    lines.append("# Semantic Router Linear Head Experiment — Report")
+    lines.append("")
+    lines.append("## 1. Exact e5 representation used")
+    lines.append("")
+    lines.append(f"- **Model**: {meta['embedder_id']}")
+    lines.append(f"- **Checkpoint**: {meta['model_path']}")
+    lines.append(f"- **Tokenizer**: {meta['tokenizer_path']}")
+    lines.append(f"- **Dimension**: {meta['dimension']}")
+    lines.append(f"- **Pooling**: {meta['pooling']}")
+    lines.append(f"- **Normalization**: {meta['normalization']}")
+    lines.append(f"- **Input template**: {meta['input_template']}")
+    lines.append("")
+
+    lines.append("## 2. Development/residual row counts")
+    lines.append("")
+    lines.append(f"- Total corpus: {meta['total_examples']}")
+    lines.append(f"- Frozen holdout: {meta['frozen_count']}")
+    lines.append(f"- Development pool: {meta['dev_count']}")
+    lines.append(f"- Fast-path resolved: {meta['fast_path_count']}")
+    lines.append(f"- Router-residual: {meta['residual_count']}")
+    lines.append("")
+    lines.append("Route distribution (full corpus):")
+    for route, count in sorted(meta["route_counts"].items()):
+        lines.append(f"  - {route}: {count}")
+    lines.append("")
+
+    lines.append("## 3. Grouped fold composition")
+    lines.append("")
+    lines.append(f"Folds: {meta['cv_folds']}")
+    for fold_id, stats in sorted(meta["fold_composition"].items()):
+        lines.append(f"  - Fold {fold_id}: eval={stats['eval_count']} train={stats['train_count']} routes={stats['eval_routes']}")
+    lines.append("")
+
+    # Regularization selection
+    lines.append("## 4. Selected regularization")
+    lines.append("")
+
+    lines.append("### Experiment A: All development examples")
+    best_C_all = all_results["best_C"]
+    lines.append(f"- Best C: {best_C_all}")
+    lines.append(f"- Mean accuracy: {fmt_pct(all_results['results_by_C'][best_C_all]['mean_accuracy'])} ± {fmt_pct(all_results['results_by_C'][best_C_all]['std_accuracy'])}")
+    lines.append(f"- Mean macro F1: {fmt_float(all_results['results_by_C'][best_C_all]['mean_macro_f1'])} ± {fmt_float(all_results['results_by_C'][best_C_all]['std_macro_f1'])}")
+    lines.append(f"- Total false actions (CV): {all_results['results_by_C'][best_C_all]['total_false_action']}")
+    lines.append("")
+
+    lines.append("### Experiment B: Router-residual only")
+    best_C_res = residual_results["best_C"]
+    lines.append(f"- Best C: {best_C_res}")
+    lines.append(f"- Mean accuracy: {fmt_pct(residual_results['results_by_C'][best_C_res]['mean_accuracy'])} ± {fmt_pct(residual_results['results_by_C'][best_C_res]['std_accuracy'])}")
+    lines.append(f"- Mean macro F1: {fmt_float(residual_results['results_by_C'][best_C_res]['mean_macro_f1'])} ± {fmt_float(residual_results['results_by_C'][best_C_res]['std_macro_f1'])}")
+    lines.append(f"- Total false actions (CV): {residual_results['results_by_C'][best_C_res]['total_false_action']}")
+    lines.append("")
+
+    # Stability across folds
+    lines.append("### Stability across folds")
+    lines.append("")
+    for C_val in C_VALUES:
+        r = all_results["results_by_C"][C_val]
+        fold_accs = [m["accuracy"] for m in r["fold_metrics"]]
+        fold_f1s = [m["macro_f1"] for m in r["fold_metrics"]]
+        lines.append(f"  C={C_val:<6} acc={fmt_pct(r['mean_accuracy'])}±{fmt_pct(r['std_accuracy'])}  f1={fmt_float(r['mean_macro_f1'])}±{fmt_float(r['std_macro_f1'])}  folds_acc={[fmt_pct(a) for a in fold_accs]}")
+    lines.append("")
+
+    # Experiment A metrics
+    lines.append("## 5. All-example CV metrics")
+    lines.append("")
+    metrics_all = all_results["full_metrics"]
+    lines.append(f"- Accuracy: {fmt_pct(metrics_all['accuracy'])}")
+    lines.append(f"- Macro F1: {fmt_float(metrics_all['macro_f1'])}")
+    lines.append(f"- False-action count: {metrics_all['false_action_count']}")
+    lines.append(f"- False-action rate: {fmt_pct(metrics_all['false_action_rate'])}")
+    lines.append(f"- Action precision: {fmt_float(metrics_all['action_precision'])}")
+    lines.append(f"- Action recall: {fmt_float(metrics_all['action_recall'])}")
+    lines.append(f"- Uncertain precision: {fmt_float(metrics_all['uncertain_precision'])}")
+    lines.append(f"- Uncertain recall: {fmt_float(metrics_all['uncertain_recall'])}")
+    lines.append("")
+    lines.append("Per-class metrics:")
+    for route in ROUTES:
+        pc = metrics_all["per_class"][route]
+        lines.append(f"  {route:<15} P={fmt_float(pc['precision'])} R={fmt_float(pc['recall'])} F1={fmt_float(pc['f1'])} (n={pc['support']})")
+    lines.append("")
+    lines.append("Confusion matrix (rows=expected, cols=predicted):")
+    header = f"{'':>15}" + "".join(f"{r:>15}" for r in ROUTES)
+    lines.append(header)
+    for i, route in enumerate(ROUTES):
+        row = f"{route:>15}" + "".join(f"{metrics_all['confusion_matrix'][i][j]:>15}" for j in range(len(ROUTES)))
+        lines.append(row)
+    lines.append("")
+
+    # Experiment B metrics
+    lines.append("## 6. Residual-only CV metrics")
+    lines.append("")
+    metrics_res = residual_results["full_metrics"]
+    lines.append(f"- Accuracy: {fmt_pct(metrics_res['accuracy'])}")
+    lines.append(f"- Macro F1: {fmt_float(metrics_res['macro_f1'])}")
+    lines.append(f"- False-action count: {metrics_res['false_action_count']}")
+    lines.append(f"- False-action rate: {fmt_pct(metrics_res['false_action_rate'])}")
+    lines.append(f"- Action precision: {fmt_float(metrics_res['action_precision'])}")
+    lines.append(f"- Action recall: {fmt_float(metrics_res['action_recall'])}")
+    lines.append(f"- Uncertain precision: {fmt_float(metrics_res['uncertain_precision'])}")
+    lines.append(f"- Uncertain recall: {fmt_float(metrics_res['uncertain_recall'])}")
+    lines.append("")
+    lines.append("Per-class metrics:")
+    for route in ROUTES:
+        pc = metrics_res["per_class"][route]
+        lines.append(f"  {route:<15} P={fmt_float(pc['precision'])} R={fmt_float(pc['recall'])} F1={fmt_float(pc['f1'])} (n={pc['support']})")
+    lines.append("")
+    lines.append("Confusion matrix (rows=expected, cols=predicted):")
+    header = f"{'':>15}" + "".join(f"{r:>15}" for r in ROUTES)
+    lines.append(header)
+    for i, route in enumerate(ROUTES):
+        row = f"{route:>15}" + "".join(f"{metrics_res['confusion_matrix'][i][j]:>15}" for j in range(len(ROUTES)))
+        lines.append(row)
+    lines.append("")
+
+    # Legacy comparison
+    lines.append("## 7. Legacy-vs-linear comparison")
+    lines.append("")
+    lines.append("### All examples")
+    lines.append(f"{'metric':<25} {'legacy':>10} {'linear e5':>10} {'delta':>10}")
+    lines.append("-" * 55)
+    # Legacy baseline from test: 52.2% overall, 40.8% residual, 19.9% false-action
+    # These are approximate since we don't have per-example legacy predictions
+    legacy_acc = 0.522
+    legacy_fa_rate = 0.199
+    legacy_macro_f1 = 0.0  # unknown precisely
+    lines.append(f"{'accuracy':<25} {fmt_pct(legacy_acc):>10} {fmt_pct(metrics_all['accuracy']):>10} {fmt_pct(metrics_all['accuracy'] - legacy_acc):>10}")
+    lines.append(f"{'macro F1':<25} {'—':>10} {fmt_float(metrics_all['macro_f1']):>10} {'—':>10}")
+    lines.append(f"{'action precision':<25} {'—':>10} {fmt_float(metrics_all['action_precision']):>10} {'—':>10}")
+    lines.append(f"{'false-action rate':<25} {fmt_pct(legacy_fa_rate):>10} {fmt_pct(metrics_all['false_action_rate']):>10} {fmt_pct(metrics_all['false_action_rate'] - legacy_fa_rate):>10}")
+    lines.append(f"{'uncertain F1':<25} {fmt_float(0.0):>10} {fmt_float(metrics_all['per_class']['uncertain']['f1']):>10} {fmt_float(metrics_all['per_class']['uncertain']['f1']):>10}")
+    lines.append("")
+
+    lines.append("### Router-residual only")
+    lines.append(f"{'metric':<25} {'legacy':>10} {'linear e5':>10} {'delta':>10}")
+    lines.append("-" * 55)
+    legacy_res_acc = 0.408
+    lines.append(f"{'accuracy':<25} {fmt_pct(legacy_res_acc):>10} {fmt_pct(metrics_res['accuracy']):>10} {fmt_pct(metrics_res['accuracy'] - legacy_res_acc):>10}")
+    lines.append(f"{'macro F1':<25} {'—':>10} {fmt_float(metrics_res['macro_f1']):>10} {'—':>10}")
+    lines.append(f"{'false-action rate':<25} {'—':>10} {fmt_pct(metrics_res['false_action_rate']):>10} {'—':>10}")
+    lines.append("")
+
+    # Fold variance
+    lines.append("## 8. Fold variance")
+    lines.append("")
+    lines.append("All-example CV:")
+    for m in all_results["results_by_C"][best_C_all]["fold_metrics"]:
+        lines.append(f"  Fold {m['fold']}: acc={fmt_pct(m['accuracy'])} f1={fmt_float(m['macro_f1'])} false_action={m['false_action']}")
+    lines.append("")
+    lines.append("Residual-only CV:")
+    for m in residual_results["results_by_C"][best_C_res]["fold_metrics"]:
+        lines.append(f"  Fold {m['fold']}: acc={fmt_pct(m['accuracy'])} f1={fmt_float(m['macro_f1'])} false_action={m['false_action']}")
+    lines.append("")
+
+    # Disagreement analysis
+    lines.append("## 9. False-action repair/new-error analysis")
+    lines.append("")
+    # The disagreement analysis requires legacy per-example predictions.
+    # Since we don't have those, we report what we can from the out-of-fold data.
+    lines.append("Note: Legacy per-example predictions were not available for this experiment.")
+    lines.append("The legacy baseline was measured in aggregate in the Go test suite.")
+    lines.append("")
+    lines.append("Learned router false-action cases (out-of-fold):")
+    for r in oof_all:
+        if r["true"] != "action" and r["predicted"] == "action":
+            lines.append(f"  {r['source_id']}: '{r['text']}' (true={r['true']}, proba(action)={r['proba'].get('action', 0):.3f})")
+    lines.append("")
+
+    # Contrast family
+    lines.append("## 10. Contrast-family results")
+    lines.append("")
+    lines.append("### Experiment A (all dev)")
+    contrast_all = all_results["contrast_family"]
+    lines.append(f"{'family':<25} {'count':>6} {'correct':>8} {'accuracy':>10} {'false_act':>10}")
+    lines.append("-" * 60)
+    for family in CONTRAST_FAMILIES:
+        if family in contrast_all:
+            c = contrast_all[family]
+            lines.append(f"{family:<25} {c['count']:>6} {c['correct']:>8} {fmt_pct(c['accuracy']):>10} {c['false_action']:>10}")
+    lines.append("")
+
+    lines.append("### Experiment B (residual only)")
+    contrast_res = residual_results["contrast_family"]
+    lines.append(f"{'family':<25} {'count':>6} {'correct':>8} {'accuracy':>10} {'false_act':>10}")
+    lines.append("-" * 60)
+    for family in CONTRAST_FAMILIES:
+        if family in contrast_res:
+            c = contrast_res[family]
+            lines.append(f"{family:<25} {c['count']:>6} {c['correct']:>8} {fmt_pct(c['accuracy']):>10} {c['false_action']:>10}")
+    lines.append("")
+
+    # Calibration
+    lines.append("## 11. Calibration metrics")
+    lines.append("")
+    lines.append("### Experiment A")
+    cal_all = all_results["calibration"]
+    lines.append(f"- ECE: {fmt_float(cal_all['ece'])}")
+    lines.append(f"- Brier score: {fmt_float(cal_all['brier'])}")
+    lines.append(f"- Log loss: {fmt_float(cal_all['log_loss'])}")
+    lines.append("")
+    lines.append("### Experiment B")
+    cal_res = residual_results["calibration"]
+    lines.append(f"- ECE: {fmt_float(cal_res['ece'])}")
+    lines.append(f"- Brier score: {fmt_float(cal_res['brier'])}")
+    lines.append(f"- Log loss: {fmt_float(cal_res['log_loss'])}")
+    lines.append("")
+
+    # Abstention curves
+    lines.append("## 12. Abstention curves")
+    lines.append("")
+    lines.append("### Experiment A (all dev)")
+    lines.append(f"{'threshold':>10} {'n_accepted':>11} {'coverage':>10} {'accuracy':>10} {'macro_f1':>10} {'false_act':>10}")
+    lines.append("-" * 62)
+    for curve in cal_all["abstention_curves"]:
+        lines.append(f"{curve['threshold']:>10.2f} {curve['n_accepted']:>11} {fmt_pct(curve['coverage']):>10} {fmt_pct(curve['accuracy']):>10} {fmt_float(curve['macro_f1']):>10} {curve['false_action_count']:>10}")
+    lines.append("")
+
+    lines.append("### Experiment B (residual only)")
+    lines.append(f"{'threshold':>10} {'n_accepted':>11} {'coverage':>10} {'accuracy':>10} {'macro_f1':>10} {'false_act':>10}")
+    lines.append("-" * 62)
+    for curve in cal_res["abstention_curves"]:
+        lines.append(f"{curve['threshold']:>10.2f} {curve['n_accepted']:>11} {fmt_pct(curve['coverage']):>10} {fmt_pct(curve['accuracy']):>10} {fmt_float(curve['macro_f1']):>10} {curve['false_action_count']:>10}")
+    lines.append("")
+
+    # Action threshold
+    lines.append("## 13. Action-threshold curve")
+    lines.append("")
+    lines.append("### Experiment A")
+    lines.append(f"{'threshold':>10} {'action_P':>10} {'action_R':>10} {'false_act':>10}")
+    lines.append("-" * 40)
+    for curve in all_results["action_threshold"]:
+        lines.append(f"{curve['threshold']:>10.2f} {fmt_float(curve['action_precision']):>10} {fmt_float(curve['action_recall']):>10} {curve['false_action_count']:>10}")
+    lines.append("")
+
+    lines.append("### Experiment B")
+    lines.append(f"{'threshold':>10} {'action_P':>10} {'action_R':>10} {'false_act':>10}")
+    lines.append("-" * 40)
+    for curve in residual_results["action_threshold"]:
+        lines.append(f"{curve['threshold']:>10.2f} {fmt_float(curve['action_precision']):>10} {fmt_float(curve['action_recall']):>10} {curve['false_action_count']:>10}")
+    lines.append("")
+
+    # Model artifact size
+    lines.append("## 14. Model artifact size and runtime cost")
+    lines.append("")
+    # Logistic regression: 6 classes × 384 features + 6 biases = 2310 parameters
+    n_params = len(ROUTES) * meta["dimension"] + len(ROUTES)
+    serialized_bytes = n_params * 4  # float32
+    lines.append(f"- Trainable parameters: {n_params}")
+    lines.append(f"  - {len(ROUTES)} classes × {meta['dimension']} features = {len(ROUTES) * meta['dimension']} weights")
+    lines.append(f"  - {len(ROUTES)} bias terms")
+    lines.append(f"- Serialized head size: {serialized_bytes} bytes ({serialized_bytes / 1024:.1f} KB)")
+    lines.append(f"- Additional inference FLOPs: {len(ROUTES) * meta['dimension']} multiply-accumulates")
+    lines.append(f"- Incremental cost (e5 already computed): ~{len(ROUTES) * meta['dimension']} FLOPs, <1µs")
+    lines.append(f"- Cost if semantic router must trigger its own e5: full ONNX inference (~{meta['dimension']} × 128 × 12 = ~590K FLOPs)")
+    lines.append("")
+
+    # Recommendation
+    lines.append("## 16. Recommendation")
+    lines.append("")
+    # Decision logic
+    all_f1 = metrics_all["macro_f1"]
+    res_f1 = metrics_res["macro_f1"]
+    res_acc = metrics_res["accuracy"]
+
+    if res_f1 > 0.5 and res_acc > 0.55:
+        verdict = "linear head sufficient"
+        detail = (f"Residual macro F1 of {fmt_float(res_f1)} and accuracy of {fmt_pct(res_acc)} "
+                  f"exceed the legacy baseline (40.8% residual accuracy) by a meaningful margin. "
+                  f"A linear head over frozen e5-small embeddings is a viable first production candidate.")
+    elif all_f1 > 0.5 and res_f1 < 0.5:
+        verdict = "need more data"
+        detail = (f"All-example F1 ({fmt_float(all_f1)}) is acceptable but residual-only F1 "
+                  f"({fmt_float(res_f1)}) drops, suggesting the contrast-family examples are "
+                  f"hard for a linear classifier. More contrastive training data may help.")
+    else:
+        verdict = "representation inadequate"
+        detail = (f"Neither all-example ({fmt_float(all_f1)}) nor residual-only ({fmt_float(res_f1)}) "
+                  f"F1 reaches the minimum viable threshold. The e5-small linear separability floor "
+                  f"is insufficient for this 6-way task. Consider a non-linear head or a different encoder.")
+
+    lines.append(f"**{verdict}**")
+    lines.append("")
+    lines.append(detail)
+    lines.append("")
+
+    # Commit hash
+    lines.append("## 17. Commit hash for experiment tooling")
+    lines.append("")
+    import subprocess
+    try:
+        result = subprocess.run(
+            ["git", "rev-parse", "HEAD"],
+            capture_output=True, text=True, cwd="/home/kami/apps/Maven"
+        )
+        lines.append(f"`{result.stdout.strip()}`")
+    except Exception:
+        lines.append("(unable to determine)")
+    lines.append("")
+
+    return "\n".join(lines)
+
+
+# ─── Main ───────────────────────────────────────────────────────────────────
+
+def run_experiment(X, y, fold_ids, examples_meta, label):
+    """Run the full experiment pipeline for one population."""
+    print(f"\n{'='*60}")
+    print(f"  Experiment: {label}")
+    print(f"  Samples: {len(y)}  Features: {X.shape[1]}  Folds: {len(set(fold_ids))}")
+    print(f"{'='*60}")
+
+    # 1. Grouped CV with regularization search
+    best_C, results_by_C = grouped_cv_experiment(X, y, fold_ids, C_VALUES, examples_meta)
+
+    # 2. Use best C to get full metrics from OOF predictions
+    oof_rows = results_by_C[best_C]["oof_predictions"]
+    oof_true = np.array([r["true"] for r in oof_rows])
+    oof_pred = np.array([r["predicted"] for r in oof_rows])
+
+    full_metrics = compute_full_metrics(oof_true, oof_pred)
+
+    # 3. Calibration
+    calibration = compute_calibration(oof_true, oof_pred, oof_rows)
+
+    # 4. Action threshold
+    action_threshold = compute_action_threshold(oof_rows)
+
+    # 5. Contrast family
+    contrast_family = compute_contrast_family(oof_rows)
+
+    return {
+        "best_C": best_C,
+        "results_by_C": results_by_C,
+        "full_metrics": full_metrics,
+        "calibration": calibration,
+        "action_threshold": action_threshold,
+        "contrast_family": contrast_family,
+        "oof_rows": oof_rows,
+    }
+
+
+if __name__ == "__main__":
+    EMBEDDING_PATH = "/tmp/mvn-experiment/embeddings.json"
+    REPORT_PATH = "/tmp/mvn-experiment/report.md"
+
+    # Load data
+    meta, examples = load_embeddings(EMBEDDING_PATH)
+    print(f"Loaded {len(examples)} examples, embedder={meta['embedder_id']}, dim={meta['dimension']}")
+
+    # Development pool only
+    dev_examples = filter_dev_pool(examples)
+    print(f"Development pool: {len(dev_examples)} examples")
+
+    # Extract features
+    X_all, y_all = extract_Xy(dev_examples)
+    fold_ids_all = get_fold_groups(dev_examples)
+
+    # Experiment A: all dev examples
+    all_results = run_experiment(X_all, y_all, fold_ids_all, dev_examples, "All development examples")
+
+    # Experiment B: residual only
+    dev_residual = filter_residual(dev_examples)
+    X_res, y_res = extract_Xy(dev_residual)
+    fold_ids_res = get_fold_groups(dev_residual)
+    residual_results = run_experiment(X_res, y_res, fold_ids_res, dev_residual, "Router-residual only")
+
+    # Generate report
+    report = generate_report(
+        meta, all_results, residual_results,
+        all_results["oof_rows"], residual_results["oof_rows"],
+    )
+
+    with open(REPORT_PATH, "w") as f:
+        f.write(report)
+    print(f"\nReport written to {REPORT_PATH}")
+
+    # Also print summary
+    print("\n" + "="*60)
+    print("  SUMMARY")
+    print("="*60)
+    print(f"  All-example:  acc={fmt_pct(all_results['full_metrics']['accuracy'])}  macro_f1={fmt_float(all_results['full_metrics']['macro_f1'])}  false_action={all_results['full_metrics']['false_action_count']}")
+    print(f"  Residual:     acc={fmt_pct(residual_results['full_metrics']['accuracy'])}  macro_f1={fmt_float(residual_results['full_metrics']['macro_f1'])}  false_action={residual_results['full_metrics']['false_action_count']}")
+    print(f"  Best C (all): {all_results['best_C']}")
+    print(f"  Best C (res): {residual_results['best_C']}")
+    print(f"  ECE (all):    {fmt_float(all_results['calibration']['ece'])}")
+    print(f"  ECE (res):    {fmt_float(residual_results['calibration']['ece'])}")
+    print(f"  Brier (all):  {fmt_float(all_results['calibration']['brier'])}")
+    print(f"  Brier (res):  {fmt_float(residual_results['calibration']['brier'])}")