perf(salience): vectorize curvature pairwise loop — 57× faster, 42% e2e (#96)

cProfile attribution (2026-05-21) identified
``core.physics.salience.SalienceOperator.compute`` as 64% of total
``ChatRuntime.chat()`` time.  Pre-fix it was a nested Python loop
over ``regions × regions`` with one ``np.linalg.norm`` call per
pair.  For N≈500 mounted-vocab regions per turn that meant ~250k
norm calls per turn, dominating end-to-end latency.

Fix: numpy broadcast for pairwise displacement, distance,
pressure-delta, and contribution.  Same math; same contract.
ULP-level reassociation drift is absorbed by the 12-decimal
precision ``_salience_address`` already used for content
addressing, and by the float32 conversion at the downstream
``SalienceMap.scores_arr`` site, so neither the content_address
nor the top-k ordering changes.

Measurements (region set: N=493, dim=5, seeded):

  vectorized:  11.78 ms/call
  old-loop:   672.30 ms/call
  speedup:    57.1×

End-to-end on 8 cognition-shape prompts:

  pre-fix:  ~970 ms/turn
  post-fix:  565 ms/turn   (-42%)

Validation:

  * 15 new tests in ``tests/test_salience_vectorize_parity.py``:
      - parity with a nested-loop reference to 1e-9 absolute on
        curvature_magnitude, gradient_vector, influence_radius
        across N ∈ {1, 2, 8, 32, 128, 493}
      - content_address byte-identical across N ∈ {1, 8, 32, 128}
      - top-16 ordering matches the reference at N ∈ {32, 128, 493}
      - empty regions returns empty map
      - single region has zero curvature
  * ``core eval cognition`` byte-identical: public 100/100/91.7/100.
  * ``core test --suite cognition`` 120/0/1, ``smoke`` 67/0.

The file's pre-existing docstring promised a Rust path
(``core_rs::physics::salience::compute_curvature``) that does not
yet exist — the numpy vectorization realizes the lift now while
keeping the Rust port a future optimization on stable semantics
(CLAUDE.md: "Rust backend parity only after Python semantics are
locked by tests").

core/physics/salience.py

@@ -58,35 +58,66 @@ class SalienceOperator:
     """
 
     def compute(self, regions: Tuple[FieldRegion, ...], cycle_index: int) -> SalienceMap:
-        """Compute local curvature by pairwise pressure-gradient deflection."""
+        """Compute local curvature by pairwise pressure-gradient deflection.
+
+        Vectorized 2026-05-21 — pre-fix this was a nested Python loop
+        over ``regions × regions`` with one ``np.linalg.norm`` call per
+        pair.  For N≈500 mounted-vocab regions per turn that meant
+        ~250k norm calls per turn, dominating ~64% of total turn time
+        (cProfile, 2026-05-21).  The math is unchanged: pairwise
+        pressure-gradient deflection.  The contract — curvature_magnitude,
+        gradient_vector, influence_radius — is preserved exactly, with
+        only ULP-level drift from float-sum reassociation (well below
+        the 12-decimal precision used by ``_salience_address`` and
+        the float32 precision used by downstream score arrays).
+        """
         if not regions:
             return SalienceMap(entries=(), cycle_index=cycle_index, content_address=_salience_address(()))
-        coords = [np.asarray(region.coordinates, dtype=np.float64) for region in regions]
+
+        # (N, D) coordinate matrix and (N,) pressure vector.
+        coords = np.stack(
+            [np.asarray(region.coordinates, dtype=np.float64) for region in regions]
+        )
+        pressures = np.asarray(
+            [region.pressure_magnitude for region in regions], dtype=np.float64
+        )
+
+        # Pairwise displacement: deltas[i, j] = coords[j] - coords[i].
+        deltas = coords[None, :, :] - coords[:, None, :]  # (N, N, D)
+        # Pairwise Euclidean distance, clamped to >= 1e-8 (matches the
+        # historical max(..., 1e-8) per-pair guard).
+        distances = np.linalg.norm(deltas, axis=-1)  # (N, N)
+        np.maximum(distances, 1e-8, out=distances)
+        # Avoid 0/0 on the diagonal; zero its contributions later.
+        np.fill_diagonal(distances, 1.0)
+
+        # Pairwise pressure deltas: |pressures[j] - pressures[i]|.
+        pressure_deltas = np.abs(pressures[None, :] - pressures[:, None])  # (N, N)
+
+        # contribution[i, j] = pressure_delta[i, j] / distance[i, j]^2
+        contributions = pressure_deltas / (distances * distances)
+        np.fill_diagonal(contributions, 0.0)  # exclude i == i
+
+        # direction[i, j] = deltas[i, j] / distance[i, j];  diagonal direction
+        # vectors are zero by construction (deltas[i, i] = 0).
+        directions = deltas / distances[..., None]  # (N, N, D)
+
+        # Per-region aggregates: sum-over-j with diagonal contributions zeroed.
+        # gradient[i] = Σ_j direction[i, j] * contribution[i, j]
+        gradients = np.einsum("ijd,ij->id", directions, contributions)
+        curvatures = contributions.sum(axis=1)  # (N,)
+        radius_num = (distances * contributions).sum(axis=1)
+        radius_den = curvatures  # identical sum
+        radii = np.where(radius_den > 0.0, radius_num / np.where(radius_den > 0, radius_den, 1.0), 0.0)
+
         entries: list[SalienceEntry] = []
         for idx, region in enumerate(regions):
-            gradient = np.zeros_like(coords[idx], dtype=np.float64)
-            curvature = 0.0
-            radius_num = 0.0
-            radius_den = 0.0
-            for jdx, neighbor in enumerate(regions):
-                if idx == jdx:
-                    continue
-                delta = coords[jdx] - coords[idx]
-                distance = max(float(np.linalg.norm(delta)), 1e-8)
-                pressure_delta = abs(float(neighbor.pressure_magnitude) - float(region.pressure_magnitude))
-                contribution = pressure_delta / (distance * distance)
-                direction = delta / distance
-                gradient += direction * contribution
-                curvature += contribution
-                radius_num += distance * contribution
-                radius_den += contribution
-            gradient_tuple = tuple(float(v) for v in gradient)
             entries.append(
                 SalienceEntry(
                     region_id=region.region_id,
-                    curvature_magnitude=float(curvature),
-                    gradient_vector=gradient_tuple,
-                    influence_radius=float(radius_num / radius_den) if radius_den > 0.0 else 0.0,
+                    curvature_magnitude=float(curvatures[idx]),
+                    gradient_vector=tuple(float(v) for v in gradients[idx]),
+                    influence_radius=float(radii[idx]),
                 )
             )
         ordered = tuple(

tests/test_salience_vectorize_parity.py

@@ -0,0 +1,154 @@
+"""Salience curvature vectorization parity (perf 2026-05-21).
+
+The pre-fix ``SalienceOperator.compute`` was a nested Python loop over
+``regions × regions``.  For N≈500 mounted-vocab regions per turn it ran
+~250k ``np.linalg.norm`` calls per turn and dominated ~64% of total
+chat() time (cProfile attribution).
+
+The vectorized version uses numpy broadcast for pairwise distance
+and contribution.  Math is unchanged, but float-sum reassociation
+can shift values at ULP level.  These tests pin:
+
+  1. Parity with a reference implementation matching the pre-fix
+     nested-loop semantics, to better than 1e-9 absolute on
+     curvature_magnitude (well below the 12-decimal precision used
+     by ``_salience_address``).
+  2. ``content_address`` equality on the byte-level — proves the
+     SHA-256 fingerprint stays stable.
+  3. Top-k ordering by curvature_magnitude is preserved.
+"""
+
+from __future__ import annotations
+
+import hashlib
+from typing import Tuple
+
+import numpy as np
+import pytest
+
+from core.physics.salience import (
+    FieldRegion,
+    SalienceEntry,
+    SalienceMap,
+    SalienceOperator,
+)
+
+
+def _reference_compute(
+    regions: Tuple[FieldRegion, ...], cycle_index: int
+) -> SalienceMap:
+    """Pre-fix nested-loop reference for parity comparison."""
+    if not regions:
+        return SalienceMap(entries=(), cycle_index=cycle_index, content_address=_ref_address(()))
+    coords = [np.asarray(region.coordinates, dtype=np.float64) for region in regions]
+    entries: list[SalienceEntry] = []
+    for idx, region in enumerate(regions):
+        gradient = np.zeros_like(coords[idx], dtype=np.float64)
+        curvature = 0.0
+        radius_num = 0.0
+        radius_den = 0.0
+        for jdx, neighbor in enumerate(regions):
+            if idx == jdx:
+                continue
+            delta = coords[jdx] - coords[idx]
+            distance = max(float(np.linalg.norm(delta)), 1e-8)
+            pressure_delta = abs(
+                float(neighbor.pressure_magnitude) - float(region.pressure_magnitude)
+            )
+            contribution = pressure_delta / (distance * distance)
+            direction = delta / distance
+            gradient += direction * contribution
+            curvature += contribution
+            radius_num += distance * contribution
+            radius_den += contribution
+        entries.append(
+            SalienceEntry(
+                region_id=region.region_id,
+                curvature_magnitude=float(curvature),
+                gradient_vector=tuple(float(v) for v in gradient),
+                influence_radius=float(radius_num / radius_den) if radius_den > 0.0 else 0.0,
+            )
+        )
+    ordered = tuple(
+        sorted(entries, key=lambda entry: (-entry.curvature_magnitude, entry.region_id))
+    )
+    return SalienceMap(entries=ordered, cycle_index=cycle_index, content_address=_ref_address(ordered))
+
+
+def _ref_address(entries: Tuple[SalienceEntry, ...]) -> str:
+    h = hashlib.sha256()
+    for entry in entries:
+        h.update(entry.region_id.encode("utf-8"))
+        h.update(f":{entry.curvature_magnitude:.12f}:".encode("ascii"))
+        h.update(",".join(f"{v:.12f}" for v in entry.gradient_vector).encode("ascii"))
+    return h.hexdigest()
+
+
+def _make_regions(n: int, dim: int = 5, seed: int = 0) -> Tuple[FieldRegion, ...]:
+    rng = np.random.default_rng(seed)
+    coords = rng.standard_normal((n, dim))
+    pressures = np.abs(rng.standard_normal(n))
+    return tuple(
+        FieldRegion(
+            region_id=f"r{i:03d}",
+            coordinates=tuple(float(v) for v in coords[i]),
+            pressure_magnitude=float(pressures[i]),
+        )
+        for i in range(n)
+    )
+
+
+@pytest.mark.parametrize("n", [1, 2, 8, 32, 128, 493])
+def test_vectorized_parity_curvature_magnitude(n: int) -> None:
+    """Vectorized curvature matches the nested-loop reference to 1e-9."""
+    regions = _make_regions(n, seed=42)
+    fast = SalienceOperator().compute(regions, cycle_index=0)
+    ref = _reference_compute(regions, cycle_index=0)
+    assert len(fast.entries) == len(ref.entries)
+    # Build region_id → curvature map for both, compare.
+    fast_by_id = {e.region_id: e for e in fast.entries}
+    ref_by_id = {e.region_id: e for e in ref.entries}
+    for rid in fast_by_id:
+        f = fast_by_id[rid]
+        r = ref_by_id[rid]
+        assert f.curvature_magnitude == pytest.approx(r.curvature_magnitude, abs=1e-9, rel=1e-9)
+        assert f.influence_radius == pytest.approx(r.influence_radius, abs=1e-9, rel=1e-9)
+        for fv, rv in zip(f.gradient_vector, r.gradient_vector):
+            assert fv == pytest.approx(rv, abs=1e-9, rel=1e-9)
+
+
+@pytest.mark.parametrize("n", [1, 8, 32, 128])
+def test_vectorized_content_address_byte_stable(n: int) -> None:
+    """SHA-256 content_address is byte-identical (12-decimal precision
+    truncation hides ULP-level reassociation drift)."""
+    regions = _make_regions(n, seed=17)
+    fast = SalienceOperator().compute(regions, cycle_index=0)
+    ref = _reference_compute(regions, cycle_index=0)
+    assert fast.content_address == ref.content_address
+
+
+@pytest.mark.parametrize("n", [32, 128, 493])
+def test_vectorized_top_k_ordering_matches(n: int) -> None:
+    """Top-k by curvature stays identical — load-bearing for the
+    walk's salience candidate set."""
+    regions = _make_regions(n, seed=99)
+    fast = SalienceOperator().compute(regions, cycle_index=0)
+    ref = _reference_compute(regions, cycle_index=0)
+    fast_top = [e.region_id for e in fast.entries[:16]]
+    ref_top = [e.region_id for e in ref.entries[:16]]
+    assert fast_top == ref_top
+
+
+def test_empty_regions_returns_empty_map() -> None:
+    fast = SalienceOperator().compute((), cycle_index=42)
+    assert fast.entries == ()
+    assert fast.cycle_index == 42
+
+
+def test_single_region_has_zero_curvature() -> None:
+    """A region with no neighbors has nothing to curve against."""
+    regions = _make_regions(1, seed=1)
+    fast = SalienceOperator().compute(regions, cycle_index=0)
+    assert len(fast.entries) == 1
+    assert fast.entries[0].curvature_magnitude == 0.0
+    assert fast.entries[0].influence_radius == 0.0