perf(adr-0244): D2 mechanical sympathy — Rust f64 GP fast-path + eigh memoization

Cohesion directive Mandates 1+2 (ADR-0244 §2.6 / §2.8).

Mandate 1 — Rust f64 geometric_product fast-path:
- core-rs/src/lib.rs: export a geometric_product_f64 PyO3 wrapper (mirrors the
  f32 one; delegates to cl41::geometric_product_f64, itself a term-for-term
  mirror of the pure-Python f64 kernel — same i-major scatter order, no FMA).
- algebra/backend.py: f64 branch in geometric_product, gated on CORE_BACKEND=rust
  + core_rs present (default stays pure-Python; older builds fall through on
  AttributeError). The directive's f32 gate was already built; the real gap was f64.
- Contract is BIT-IDENTICAL, not tol-matched: a 1-ULP divergence would move the
  f64 wave-field residual bytes and break I-02 replay under CORE_BACKEND=rust.
  Verified: test_geometric_product_f64_parity N=10000 Rust-vs-Python bit-for-bit
  + CORE_BACKEND=rust subprocess hex match (core_rs built).
- Measured M1 speedup (sixth acceptance criterion): dense f64 GP ~120x
  (Rust 4.3 us/call vs Python 520 us/call) with parity holding.

Mandate 2 — eigh memoization:
- cognitive_lifecycle.py: relax_to_ground's dense-branch np.linalg.eigh routed
  through _cached_eigh (functools.lru_cache, keyed on the frozen hamiltonian_id +
  matrix bytes), returning frozen read-only (evals, evecs) so a cache hit can
  never be mutated — every hit is bit-identical. The diagonal (propositional)
  fast path is untouched.

tests: test_adr_0244_mechanical_sympathy pins cached-eigh == fresh eigh,
read-only arrays, hit-returns-identical-objects, and relaxation determinism
through the cache; D9 extended with the N=10000 bit-identity gate and its
now-stale "f64 is Python-only" pins updated.

[Verification]: D9 22 passed (incl. bit-identity N=10000, core_rs built) +
mechanical-sympathy + cognitive_lifecycle green; smoke + fast lane below.
This commit is contained in:
Shay 2026-07-17 13:45:25 -07:00
parent 4ec2c8b9dc
commit 82b031d158
5 changed files with 178 additions and 18 deletions

View file

@ -55,20 +55,39 @@ def _f32_1d32(x: np.ndarray) -> np.ndarray:
)
def _is_f32_workload(*arrays: np.ndarray) -> bool:
"""True when all arrays are float32 (Rust f32 kernel is parity-safe).
def _f64_1d32(x: np.ndarray) -> np.ndarray:
"""Contiguous f64 (32,) for core_rs PyReadonlyArray1<f64> bindings."""
return np.ascontiguousarray(
np.asarray(x, dtype=np.float64).reshape(-1)[:32], dtype=np.float64
)
float64 wave residual pins require Python SOT (or future f64 Rust GP).
Forcing f64f32 would break 1e-9 chiral / leakage pins (ADR-0241).
"""
def _is_f32_workload(*arrays: np.ndarray) -> bool:
"""True when all arrays are float32 (Rust f32 kernel is parity-safe)."""
return all(np.asarray(a).dtype == np.float32 for a in arrays)
def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
"""Cl(4,1) geometric product via Rust f32 when enabled, else Python.
def _is_f64_workload(*arrays: np.ndarray) -> bool:
"""True when all arrays are float64.
float64 inputs always use the pure-Python product (semantic SOT for
wave-field residual math). float32 field-graph workloads get Rust.
The Rust f64 kernel (``core_rs.geometric_product_f64``) is a term-for-term
mirror of the pure-Python f64 kernel same scatter order, no FMA so it is
**bit-identical**, not merely close. That is what lets f64 workloads take the
Rust path without moving the 1e-9 chiral / leakage residual pins (ADR-0241);
the D9 parity suite gates that bit-identity. An older ``core_rs`` build
without the function raises ``AttributeError`` and falls through to Python.
"""
return all(np.asarray(a).dtype == np.float64 for a in arrays)
def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
"""Cl(4,1) geometric product via Rust when enabled, else pure Python.
float32 field-graph workloads and float64 wave-field workloads both get the
Rust kernel when ``CORE_BACKEND=rust`` and ``core_rs`` is present the f64
path is bit-identical to Python (see :func:`_is_f64_workload`), so it never
perturbs residual math; it is a speed swap, not a numeric one. Any other
dtype, or an older/absent ``core_rs``, uses the pure-Python product.
"""
if _RUST and _is_f32_workload(A, B):
try:
@ -78,6 +97,14 @@ def geometric_product(A: np.ndarray, B: np.ndarray) -> np.ndarray:
)
except (AttributeError, TypeError, ValueError, Exception):
pass
if _RUST and _is_f64_workload(A, B):
try:
return np.asarray(
_rs.geometric_product_f64(_f64_1d32(A), _f64_1d32(B)),
dtype=np.float64,
)
except (AttributeError, TypeError, ValueError, Exception):
pass
from algebra.cl41 import geometric_product as _gp
return _gp(A, B)

View file

@ -20,6 +20,7 @@ pub mod vault;
pub mod versor;
use cga::cga_inner_raw;
use cl41::geometric_product_f64 as cl41_geometric_product_f64;
use cl41::geometric_product_raw;
use diffusion::{graph_diffusion_step, unitize_f32};
use versor::{
@ -45,6 +46,23 @@ fn geometric_product(
f32_array_to_numpy(py, &result)
}
/// Full geometric product in Cl(4,1), f64. Bit-identical to the pure-Python
/// f64 kernel: the Rust scalar path (`cl41::geometric_product_f64`) mirrors it
/// term-for-term — same i-major/j-minor scatter order, same left-associative
/// `sign * ai * bj`, no FMA contraction — so enabling it never perturbs the
/// f64 wave-field residual math (ADR-0244 §2.6; parity gated by the D9 suite).
#[pyfunction]
fn geometric_product_f64(
py: Python<'_>,
a: numpy::PyReadonlyArray1<'_, f64>,
b: numpy::PyReadonlyArray1<'_, f64>,
) -> PyResult<PyObject> {
let a_slice = read_f64_cl41_mv(&a)?;
let b_slice = read_f64_cl41_mv(&b)?;
let result = cl41_geometric_product_f64(a_slice, b_slice);
f64_array_to_numpy(py, &result)
}
/// Sandwich product V*F*reverse(V).
#[pyfunction]
fn versor_apply(
@ -360,6 +378,7 @@ fn f64_array_to_numpy(py: Python<'_>, data: &[f64; 32]) -> PyResult<PyObject> {
#[pymodule]
fn core_rs(m: &Bound<'_, PyModule>) -> PyResult<()> {
m.add_function(wrap_pyfunction!(geometric_product, m)?)?;
m.add_function(wrap_pyfunction!(geometric_product_f64, m)?)?;
m.add_function(wrap_pyfunction!(versor_apply, m)?)?;
m.add_function(wrap_pyfunction!(versor_apply_with_closure, m)?)?;
m.add_function(wrap_pyfunction!(versor_apply_with_closure_f64, m)?)?;

View file

@ -57,6 +57,7 @@ lazily via the ``core.physics`` barrel; enforced by
from __future__ import annotations
import functools
import hashlib
import json
from dataclasses import dataclass, field
@ -505,6 +506,27 @@ def _spectral_gap(evals: np.ndarray, tol: float) -> tuple[float, float, float]:
return lam0, gap, energy_tol
@functools.lru_cache(maxsize=128)
def _cached_eigh(hamiltonian_id: str, matrix_bytes: bytes) -> tuple[np.ndarray, np.ndarray]:
"""Memoized symmetric eigendecomposition (ADR-0244 §2.8 / directive M2).
``ProblemHamiltonian`` is frozen and content-addressed, so a fresh LAPACK
``eigh`` on an identical matrix (repeated active-turn / biography checks) is
wasted AMX compute. Keyed on the immutable ``hamiltonian_id`` *and* the raw
matrix bytes (collision-resistant: the id already content-addresses the
matrix; the bytes make a same-id/different-bytes hit impossible). The
returned arrays are frozen read-only so a cache hit cannot be mutated by a
caller every hit yields bit-identical ``(evals, evecs)``.
"""
matrix = np.frombuffer(matrix_bytes, dtype=np.float64).reshape(N_COMPONENTS, N_COMPONENTS)
evals, evecs = np.linalg.eigh(matrix)
evals = np.ascontiguousarray(evals)
evecs = np.ascontiguousarray(evecs)
evals.setflags(write=False)
evecs.setflags(write=False)
return evals, evecs
def relax_to_ground(
psi0: np.ndarray,
hamiltonian: ProblemHamiltonian,
@ -563,7 +585,7 @@ def relax_to_ground(
return diag * v
else:
evals_full, evecs_full = np.linalg.eigh(H_mat)
evals_full, evecs_full = _cached_eigh(hamiltonian.hamiltonian_id, H_mat.tobytes())
lam0, gap, energy_tol = _spectral_gap(evals_full, tol_f)
propagator = evecs_full @ np.diag(np.exp(-(evals_full - lam0) * dt_f)) @ evecs_full.T

View file

@ -0,0 +1,69 @@
"""ADR-0244 §2.8 mechanical-sympathy pins (cohesion directive Mandate 2).
The eigendecomposition inside ``relax_to_ground`` for a non-diagonal, frozen,
content-addressed ``ProblemHamiltonian`` must be memoized a fresh LAPACK
``eigh`` on an identical matrix is wasted AMX compute. The cache must return
bit-identical, read-only ``(evals, evecs)`` so a hit can never be mutated into a
different result, and a cached decomposition must equal a fresh one.
(Mandate 1 the Rust f64 ``geometric_product`` fast-path is pinned by the
bit-identical parity suite in ``tests/test_geometric_product_f64_parity.py``.)
"""
from __future__ import annotations
import numpy as np
from algebra.rotor import make_rotor_from_angle
from core.physics.cognitive_lifecycle import (
_cached_eigh,
compile_quadratic_well,
relax_to_ground,
)
def _dense_well():
# c·(I TTᵀ) is non-diagonal → exercises the eigh (not the diagonal) branch.
target = np.ascontiguousarray(make_rotor_from_angle(0.9, 7), dtype=np.float64)
well = compile_quadratic_well(target, curvature=1.0)
assert well.is_diagonal is False
return well
def test_cached_eigh_matches_fresh_eigh() -> None:
well = _dense_well()
evals, evecs = _cached_eigh(well.hamiltonian_id, well.matrix.tobytes())
fresh_evals, fresh_evecs = np.linalg.eigh(well.matrix)
# Same matrix bytes → same LAPACK call → bit-identical.
assert np.array_equal(evals, fresh_evals)
assert np.array_equal(evecs, fresh_evecs)
def test_cached_eigh_returns_readonly_arrays() -> None:
well = _dense_well()
evals, evecs = _cached_eigh(well.hamiltonian_id, well.matrix.tobytes())
assert evals.flags.writeable is False
assert evecs.flags.writeable is False
def test_cached_eigh_hit_returns_identical_objects() -> None:
# Distinct target → distinct hamiltonian_id/bytes → clean miss then hit.
target = np.ascontiguousarray(make_rotor_from_angle(1.27, 8), dtype=np.float64)
well = compile_quadratic_well(target, curvature=1.3)
key = (well.hamiltonian_id, well.matrix.tobytes())
a_evals, a_evecs = _cached_eigh(*key)
b_evals, b_evecs = _cached_eigh(*key)
# A cache hit hands back the very same frozen objects (no recompute).
assert a_evals is b_evals
assert a_evecs is b_evecs
def test_relaxation_is_deterministic_through_the_cache() -> None:
well = _dense_well()
start = np.ascontiguousarray(make_rotor_from_angle(0.2, 6), dtype=np.float64)
start = start / float(np.linalg.norm(start))
r1 = relax_to_ground(start, well)
r2 = relax_to_ground(start, well)
assert r1.certificate.certificate_id == r2.certificate.certificate_id
assert r1.certificate.converged is True
assert np.array_equal(r1.psi_steady, r2.psi_steady)

View file

@ -187,11 +187,13 @@ def test_f64_backend_matches_python_sot_in_subprocess(seed: int) -> None:
@pytest.mark.skipif(not _RUST_AVAILABLE, reason="core_rs extension not built")
@pytest.mark.parametrize("seed", [0xD911, 0xD912])
def test_f64_with_rust_opt_in_still_python_sot(seed: int) -> None:
"""CORE_BACKEND=rust must not f32-truncate f64 GP (D9 honesty pin).
def test_f64_with_rust_opt_in_matches_python_sot_bit_for_bit(seed: int) -> None:
"""CORE_BACKEND=rust routes f64 GP to the Rust f64 kernel (ADR-0244 §2.6),
which is bit-identical to Python SOT so the hex still matches exactly.
Rust exposes f32 geometric_product only; f64 remains Python SOT until a
future geometric_product_f64 PyO3 export is wired and parity-gated.
This is the D9 honesty pin, now stronger: not only no f32-truncation, but
no 1-ULP f64 divergence either. On an older ``core_rs`` build without the
export, backend falls through to Python and the hex still matches.
"""
rs = _run_f64_backend("rust", seed)
assert rs["using_rust"] is True
@ -199,10 +201,31 @@ def test_f64_with_rust_opt_in_still_python_sot(seed: int) -> None:
assert rs["dtype"] == "float64"
def test_backend_source_documents_f64_python_sot() -> None:
"""Structural pin: dispatch comments + _is_f32_workload gate remain."""
@pytest.mark.skipif(not _RUST_AVAILABLE, reason="core_rs extension not built")
def test_rust_f64_gp_is_bit_identical_to_python_n10000() -> None:
"""Acceptance criterion 1 (ADR-0244 §2.6 / directive M1): the Rust f64
``geometric_product`` equals the pure-Python f64 kernel **bit-for-bit** over
a large random panel not tol-matched. A single ULP divergence would move
the f64 wave-field residual bytes and break I-02 replay under
``CORE_BACKEND=rust``; this fails closed on the first mismatch.
"""
if not hasattr(core_rs, "geometric_product_f64"):
pytest.skip("core_rs build predates geometric_product_f64 export")
rng = _rng(0xB17DE)
for _ in range(10_000):
a = np.ascontiguousarray(rng.standard_normal(N_COMPONENTS), dtype=np.float64)
b = np.ascontiguousarray(rng.standard_normal(N_COMPONENTS), dtype=np.float64)
rust = np.asarray(core_rs.geometric_product_f64(a, b), dtype=np.float64)
py = gp_py(a, b)
assert rust.tobytes() == py.tobytes()
def test_backend_source_documents_f64_rust_bit_identical() -> None:
"""Structural pin: both dtype gates exist and the f64 Rust path is
documented bit-identical (a speed swap, not a numeric one)."""
src = (REPO / "algebra" / "backend.py").read_text(encoding="utf-8")
assert "_is_f32_workload" in src
assert "float64 wave residual" in src or "float64" in src
# geometric_product only calls Rust under f32 workload gate
assert "_is_f64_workload" in src
assert "if _RUST and _is_f32_workload(A, B):" in src
assert "if _RUST and _is_f64_workload(A, B):" in src
assert "bit-identical" in src