From 2a2ef9ce49758eac7c1bee086d4528232d9ed0dc Mon Sep 17 00:00:00 2001 From: Shay Date: Wed, 20 May 2026 21:29:42 -0700 Subject: [PATCH] =?UTF-8?q?perf(salience):=20vectorize=20curvature=20pairw?= =?UTF-8?q?ise=20loop=20=E2=80=94=2057=C3=97=20faster,=2042%=20e2e=20(#96)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 | 75 ++++++++---- tests/test_salience_vectorize_parity.py | 154 ++++++++++++++++++++++++ 2 files changed, 207 insertions(+), 22 deletions(-) create mode 100644 tests/test_salience_vectorize_parity.py diff --git a/core/physics/salience.py b/core/physics/salience.py index 56f77d8c..d0fa14e7 100644 --- a/core/physics/salience.py +++ b/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( diff --git a/tests/test_salience_vectorize_parity.py b/tests/test_salience_vectorize_parity.py new file mode 100644 index 00000000..18744811 --- /dev/null +++ b/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