Add tight-pair bond ultrametricity scores

For each tight pair (a,b) — mutual nearest neighbors or the k smallest-D pairs —
average the Rammal index u over triples where (a,b) is the shortest edge; the
bond score 1-<u> separates the ultrametric backbone (high) from specific
contacts (low). Compared against a null computed the same way on a distance
matrix with shuffled off-diagonals; low-score pairs are expected to match the
most-negative D-C residual pairs, and bimodality marks the backbone-vs-contact
split.

tight_pair_scores returns the pair indices, per-bond scores, and null so they
can be cross-referenced against the R eigenvectors. Persisted to the .npz/summary
and surfaced in the dashboard with live method (MNN / smallest-D) and k controls,
a score-vs-null histogram, and a flagged low-score table tagged with residual
overlap.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

apps/ultrametricity_view.py

@@ -23,7 +23,7 @@ from scipy.cluster.hierarchy import dendrogram, linkage, cophenet
 from scipy.spatial.distance import squareform
 
 from transformer_analysis.ultrametricity import (
-    to_distance, synthetic_ultrametric, null_distance_samples)
+    to_distance, synthetic_ultrametric, null_distance_samples, tight_pair_scores)
 
 METRICS = {"Frobenius cosine": "frob_cosine",
            "Jensen-Shannon": "hist_jensen_shannon"}
@@ -59,6 +59,13 @@ def _null(d_head: int, d_model: int, metric: str):
     return null_distance_samples(d_head, d_model, metric)
 
 
+@st.cache_data
+def _tight(exp_dir: str, run_key: str, metric: str, method: str, k: int):
+    _, U, D, _ = _load(exp_dir, run_key, metric)
+    keys = [tuple(x) for x in U["head_index"]]
+    return tight_pair_scores(D, keys, method=method, k_smallest=k)
+
+
 def _grid(values, keys, n_layers, n_heads):
     g = np.full((n_layers, n_heads), np.nan)
     for v, (l, h) in zip(values, keys):
@@ -403,6 +410,51 @@ def render(exp_dir: Path):
                    "energy in a few eigenvalues ⇒ a handful of global factors on top of "
                    "the tree; high Gini/kurtosis ⇒ sparse specific pairs; otherwise noise.")
 
+    # ── Tight-pair bond scores (backbone vs contact) ──
+    st.subheader("Tight-pair bond scores")
+    tcol = st.columns(2)
+    method_label = tcol[0].selectbox("Tight pairs", ["smallest-D", "mutual NN"],
+                                     key="tight_method")
+    method = "smallest" if method_label == "smallest-D" else "mnn"
+    kk = tcol[1].slider("k (smallest-D)", 20, 500, 200, 20, key="tight_k",
+                        disabled=(method != "smallest"))
+    tg = _tight(str(exp_dir), run_key, metric, method, kk)
+    sc = tg["scores"][np.isfinite(tg["scores"])]
+    nl = tg["null_scores"][np.isfinite(tg["null_scores"])]
+    if len(sc):
+        bins = np.linspace(0, 1, 41)
+        cen = 0.5 * (bins[:-1] + bins[1:])
+        figt = go.Figure()
+        for vals, name, color, fc in (
+            (sc, f"tight pairs ({tg['method']})", "#19D3F3", "rgba(25,211,243,0.40)"),
+            (nl, "null (shuffled D)", "#FF6692", "rgba(255,102,146,0.35)")):
+            if len(vals):
+                h, _ = np.histogram(vals, bins=bins, density=True)
+                figt.add_scatter(x=cen, y=h, name=name, line_shape="hv",
+                                 fill="tozeroy", line=dict(color=color), fillcolor=fc)
+        figt.update_layout(title="bond score 1−⟨u⟩ over tight pairs", height=340,
+                           xaxis_title="1−⟨u⟩  (high = backbone, low = contact)",
+                           yaxis_title="density", legend=dict(x=0.01, y=0.99))
+        st.plotly_chart(figt, use_container_width=True)
+        st.caption("Per tight pair (a,b): 1−⟨u⟩ averaged over triples where (a,b) is the "
+                   "shortest edge. High ⇒ clean ultrametric sibling (backbone); low ⇒ a "
+                   "close pair whose triples break ultrametricity (a specific contact). "
+                   "Bimodality is the backbone-vs-contact split. Null = same procedure on "
+                   "a distance matrix with shuffled off-diagonals.")
+
+        low = pd.DataFrame(tg.get("low_flagged", []))
+        if len(low):
+            respairs = {frozenset((rp_["i"], rp_["j"]))
+                        for rp_ in m.get("residual_pairs", []) if "i" in rp_}
+            low["in_R_residual_top"] = [frozenset((r.i, r.j)) in respairs
+                                        for r in low.itertuples()]
+            st.caption(f"Lowest-score tight pairs (median {tg['score_median']:.3f} vs "
+                       f"null {tg['null_median']:.3f}, {tg['n_pairs']} pairs) — these "
+                       "should match the most-negative R residual pairs.")
+            st.dataframe(low[["layer_1", "head_1", "layer_2", "head_2", "delta_layer",
+                              "distance", "score", "in_R_residual_top"]],
+                         use_container_width=True, height=320)
+
     # ── Flat geometry alternative: classical MDS ──
     st.subheader("Euclidean alternative (classical MDS)")
     em = m.get("embedding")

scripts/run_ultrametricity.py

@@ -28,6 +28,7 @@ from transformer_analysis.ultrametricity import (
     residual_decomposition,
     embedding_analysis,
     linkage_comparison,
+    tight_pair_scores,
 )
 
 
@@ -107,6 +108,8 @@ def main():
             quad = quadruple_statistics(D, sample_cap=args.quad_sample_cap,
                                         seed=args.seed, tol=args.tol)
             lk = linkage_comparison(D)
+            tight = tight_pair_scores(D, keys, method="smallest", k_smallest=300,
+                                      seed=args.seed)
 
             dense = D.shape[0] <= args.max_dense
             resid = residual_decomposition(D, clust["linkage"]) if dense else None
@@ -131,6 +134,9 @@ def main():
                 mds_eigs=emb["eigenvalues"] if emb else np.array([]),
                 mds_top_eigvecs=emb["top_eigvecs"] if emb else np.zeros((0, 0)),
                 mds_top_eigs=emb["top_eigs"] if emb else np.array([]),
+                tight_scores=tight["scores"],
+                tight_null_scores=tight["null_scores"],
+                tight_pair_idx=tight["pair_idx"],
             )
 
             summary["metrics"][metric] = {
@@ -155,6 +161,8 @@ def main():
                      "gini", "kurtosis", "verdict")},
                 "embedding": emb and {k: emb[k] for k in
                     ("neg_energy_frac", "var_frac_2d", "var_frac_3d", "eff_dim")},
+                "tight": {k: tight[k] for k in
+                    ("method", "n_pairs", "score_median", "null_median", "low_flagged")},
             }
             log = {
                 f"{metric}.mean_u": stats["mean_u"],
@@ -163,6 +171,7 @@ def main():
                 f"{metric}.triangle_violation_frac": stats["triangle_violation_frac"],
                 f"{metric}.mean_u4": quad["mean_u4"],
                 f"{metric}.frac_additive": quad["frac_additive"],
+                f"{metric}.tight_score_median": tight["score_median"],
             }
             if resid:
                 log[f"{metric}.res_energy_frac"] = resid["res_energy_frac"]

src/transformer_analysis/ultrametricity.py

@@ -229,6 +229,91 @@ def cluster_analysis(D, keys, top_k=50):
     }
 
 
+# ── Tight-pair bond scores (backbone vs contact) ──────────────────────
+
+def _candidate_pairs(M, method, k_smallest):
+    n = M.shape[0]
+    if method == "smallest":
+        iu = np.triu_indices(n, k=1)
+        order = np.argsort(M[iu])[:k_smallest]
+        return [(int(iu[0][o]), int(iu[1][o])) for o in order]
+    # mutual nearest neighbors
+    nn = (M + np.eye(n) * (M.max() + 1.0)).argmin(axis=1)
+    return [(i, int(nn[i])) for i in range(n) if nn[nn[i]] == i and i < nn[i]]
+
+
+def _bond_scores(M, pairs):
+    scores = np.full(len(pairs), np.nan)
+    counts = np.zeros(len(pairs), dtype=int)
+    for idx, (a, b) in enumerate(pairs):
+        dab = M[a, b]
+        dac, dbc = M[a], M[b]
+        mask = (dac >= dab) & (dbc >= dab)     # triples where (a,b) is shortest edge
+        mask[a] = mask[b] = False
+        if mask.any():
+            d2 = np.minimum(dac, dbc)[mask]
+            d3 = np.maximum(dac, dbc)[mask]
+            denom = d3 - dab
+            u = np.where(denom > 0, (d3 - d2) / denom, 0.0)
+            scores[idx] = 1.0 - float(u.mean())
+            counts[idx] = int(mask.sum())
+    return scores, counts
+
+
+def tight_pair_scores(D, keys, method="mnn", k_smallest=300, seed=0, top_flag=25):
+    """Per-bond ultrametricity score for tight (close) head pairs.
+
+    For each tight pair (a,b) — mutual nearest neighbors, or the k smallest-D
+    pairs — average the Rammal index u over all triples (a,b,c) in which (a,b)
+    is the shortest edge; the bond score is 1 - <u>.  High = a clean ultrametric
+    sibling (backbone); low = a close pair whose triples break ultrametricity
+    (a specific contact, matching the most-negative D-C residuals).  A null is
+    computed the same way on a distance matrix with shuffled off-diagonals.
+    Bimodality in the score distribution is the backbone-vs-contact split.
+    """
+    n = D.shape[0]
+    layers = np.array([k[0] for k in keys])
+    heads = np.array([k[1] for k in keys])
+
+    pairs = _candidate_pairs(D, method, k_smallest)
+    scores, counts = _bond_scores(D, pairs)
+
+    rng = np.random.default_rng(seed)
+    iu = np.triu_indices(n, k=1)
+    vals = D[iu].copy()
+    rng.shuffle(vals)
+    Dn = np.zeros_like(D)
+    Dn[iu] = vals
+    Dn[(iu[1], iu[0])] = vals
+    null_scores, _ = _bond_scores(Dn, _candidate_pairs(Dn, method, k_smallest))
+
+    recs = []
+    for (a, b), s, c in zip(pairs, scores, counts):
+        recs.append({
+            "i": int(a), "j": int(b),
+            "layer_1": int(layers[a]), "head_1": int(heads[a]),
+            "layer_2": int(layers[b]), "head_2": int(heads[b]),
+            "delta_layer": int(abs(layers[a] - layers[b])),
+            "distance": float(D[a, b]),
+            "score": float(s) if np.isfinite(s) else None,
+            "n_triples": int(c),
+        })
+    valid = [r for r in recs if r["score"] is not None]
+    low_flagged = sorted(valid, key=lambda r: r["score"])[:top_flag]
+
+    return {
+        "method": method,
+        "n_pairs": len(pairs),
+        "pair_idx": np.array([[r["i"], r["j"]] for r in recs], dtype=int)
+                    if recs else np.zeros((0, 2), dtype=int),
+        "scores": scores,
+        "null_scores": null_scores,
+        "score_median": float(np.nanmedian(scores)) if len(scores) else float("nan"),
+        "null_median": float(np.nanmedian(null_scores)) if len(null_scores) else float("nan"),
+        "low_flagged": low_flagged,
+    }
+
+
 # ── Is it the wrong tree, or extra structure on top? ──────────────────
 
 def _gini(x):