diff --git a/CHANGELOG.md b/CHANGELOG.md index 08646f9a..ed814d3c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,6 +35,24 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). eigensolver paths use, so one quantity is not measured two different ways. ### Changed +- **Relaxation propagation now computes the action `exp(tL) rho0` with + `scipy.sparse.linalg.expm_multiply` instead of materialising the full dense + matrix exponential (issue #162).** Exact NumPy-`linspace` grids use SciPy's + interval action API; arbitrary grids use independent actions from the same + initial state, avoiding stepwise error accumulation. The low-level propagator + accepts dense and SciPy sparse matrices and preserves the exact `t=0` state. + The #156 fail-closed contract remains in force: non-finite `L*t`, an + unrepresentable 1-norm, action failures, or non-finite output states raise + `UnrepresentableTrajectoryError` instead of flowing into fitting. The + relaxation result now records the backend plus trajectory-wide trace error, + Hermiticity defect and minimum Hermitian eigenvalue; positivity drift is + measured, never clipped. Dense-reference, sparse-vs-dense, non-uniform-grid, + and `L -> cL, t -> t/c` metamorphic tests bound the change. This is a + NUMERICAL METHODOLOGY change: floating-point trajectories may differ at + round-off level from the old full-`expm` route, although the mathematical + evolution is unchanged. The top-level orchestrator remains dense-only; sparse + support here is scoped to the propagation primitive. No manifest/schema + field changes. - **The Gaussian likelihood behind AICc is evaluated in log-RSS space, and an exact-zero RSS is now an explicit abstention (issue #135).** This is a METHODOLOGY change with user-visible consequences: it can reorder AICc, diff --git a/CITATION.cff b/CITATION.cff index ddfe2715..c2a90cae 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -153,6 +153,18 @@ abstract: >- # than reduction_rtol * (eps * ||L||_F + d**3 * 2**-1074), reduction_rtol # defaulting to the existing ZERO_MODE_EPS_FACTOR. Like the restriction itself # this moves no reported number: no diagnostic consumes it yet. +# Also pending for the next cut: relaxation trajectories use the +# Al-Mohy--Higham exponential-action route exposed by SciPy +# `sparse.linalg.expm_multiply` instead of materialising `expm(tL)` and then +# multiplying by the initial state (issue #162). Exact linspace grids use the +# interval API; non-uniform grids use independent actions from the same initial +# state. The mathematical semigroup is unchanged, but floating-point values may +# differ at round-off level, so this is a NUMERICAL METHODOLOGY change rather +# than a pure refactor. Backend identity plus trace/Hermiticity/positivity-drift +# measurements are recorded on RelaxationResult. The low-level action path is +# dense/sparse; this must NOT be described as making diagnose() sparse-native. +# The #156 representability guard remains a prerequisite and fails closed before +# non-finite or unrepresentable action inputs can enter downstream fits. keywords: - open quantum systems - Lindblad diff --git a/README.md b/README.md index 9190dcdd..deafcf43 100644 --- a/README.md +++ b/README.md @@ -237,6 +237,17 @@ LiouScope is built for paper-grade reproducibility: - **Anchor tests.** `tests/test_anchors.py` locks the numerical anchors that paper figures depend on; changes to physics code that move these values are caught in CI. - **Paper-figure pipeline.** `figures/generate_all.py` regenerates Fig 1-3 deterministically. +- **Relaxation trajectories use exponential action, not a materialised propagator.** + `diagnostics.relaxation` evolves `exp(tL) rho0` with + `scipy.sparse.linalg.expm_multiply`: exact `linspace` grids use its interval + API, while arbitrary/non-uniform grids evaluate each requested time directly + from the same initial state. `RelaxationResult.trajectory_backend` records + which path was used, together with maximum trace error, maximum Hermiticity + defect, and the minimum eigenvalue of the Hermitian part seen along the + trajectory. These are audit measurements, not silent positivity repairs. + The low-level propagation primitive accepts dense and SciPy sparse matrices; + the top-level `diagnose()` orchestrator remains dense-only until its separate + sparse solver path is wired. Unrepresentable propagation fails closed. --- diff --git a/benchmarks/issue162_expm_action.py b/benchmarks/issue162_expm_action.py new file mode 100644 index 00000000..e570da4d --- /dev/null +++ b/benchmarks/issue162_expm_action.py @@ -0,0 +1,81 @@ +"""Issue #162 benchmark: full dense expm vs exponential action. + +Run manually: + python benchmarks/issue162_expm_action.py + +The script is evidence generation, not a performance gate. It prints elapsed +time, peak RSS (where ``resource`` is available), and dense/action agreement. +Each backend runs in a fresh child process so peak RSS is not contaminated by +the backend that ran first. +""" + +from __future__ import annotations + +import multiprocessing as mp +import time + +import numpy as np +import scipy.linalg as sla +import scipy.sparse.linalg as spla + + +def _matrix(n: int, seed: int = 162) -> tuple[np.ndarray, np.ndarray]: + rng = np.random.default_rng(seed) + A = rng.standard_normal((n, n)) + 1j * rng.standard_normal((n, n)) + # Shift left so exp(A) is not dominated by explosive growth. + A -= (float(n) + 1.0) * np.eye(n) + b = rng.standard_normal(n) + 1j * rng.standard_normal(n) + return np.asarray(A, dtype=complex), np.asarray(b, dtype=complex) + + +def _peak_rss_mb() -> float: + try: + import resource + + rss = float(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) + # Linux reports KiB; macOS reports bytes. This repo CI/benchmark target + # is Linux, but keep the output approximately useful on macOS. + return rss / (1024.0 if rss < 1.0e8 else 1024.0**2) + except ImportError: + return float("nan") + + +def _worker(kind: str, n: int, q: mp.Queue) -> None: + A, b = _matrix(n) + started = time.perf_counter() + if kind == "dense": + out = sla.expm(A) @ b + elif kind == "action": + out = spla.expm_multiply(A, b) + else: + raise ValueError(kind) + q.put((kind, time.perf_counter() - started, _peak_rss_mb(), out)) + + +def _run(kind: str, n: int) -> tuple[float, float, np.ndarray]: + q: mp.Queue = mp.Queue() + proc = mp.Process(target=_worker, args=(kind, n, q)) + proc.start() + proc.join() + if proc.exitcode != 0: + raise RuntimeError(f"{kind} child failed with exit code {proc.exitcode}") + _kind, elapsed, rss, out = q.get() + return float(elapsed), float(rss), np.asarray(out) + + +def main() -> None: + for n in (64, 128, 256): + dense_t, dense_rss, dense = _run("dense", n) + action_t, action_rss, action = _run("action", n) + rel = float( + np.linalg.norm(action - dense) + / max(np.linalg.norm(dense), np.finfo(float).tiny) + ) + print( + f"n={n:4d} dense={dense_t:9.4f}s/{dense_rss:9.1f}MB " + f"action={action_t:9.4f}s/{action_rss:9.1f}MB relerr={rel:.3e}" + ) + + +if __name__ == "__main__": + main() diff --git a/src/liouscope/_types.py b/src/liouscope/_types.py index 67ddbf22..d71d461f 100644 --- a/src/liouscope/_types.py +++ b/src/liouscope/_types.py @@ -223,6 +223,13 @@ class RelaxationResult: # uses. Additive + defaulted so older callers / serialised reports stay valid. beta_D_linear: float = float("nan") linear_fit_model: str = "none" + # Issue #162: trajectory backend and numerical-state audit. Additive + + # defaulted so older synthetic callers and serialised report objects remain + # valid. These fields are measurements, not physicality gates. + trajectory_backend: str = "legacy_dense_expm" + trajectory_max_trace_error: float = float("nan") + trajectory_max_hermiticity_defect: float = float("nan") + trajectory_min_eigenvalue: float = float("nan") @dataclass(frozen=True, slots=True, kw_only=True) diff --git a/src/liouscope/diagnostics/relaxation.py b/src/liouscope/diagnostics/relaxation.py index cdf0870e..71040609 100644 --- a/src/liouscope/diagnostics/relaxation.py +++ b/src/liouscope/diagnostics/relaxation.py @@ -5,9 +5,12 @@ from __future__ import annotations import warnings +from typing import Any import numpy as np import scipy.linalg as sla +import scipy.sparse as sp +import scipy.sparse.linalg as spla from .._consts import EPS_SUPP from .._types import FitResult, RelaxationResult @@ -36,45 +39,131 @@ class UnrepresentableTrajectoryError(RuntimeError): """The requested relaxation propagation is not representable reliably. This is a numerical-domain failure, not a statement about the underlying - GKSL dynamics. Returning a non-finite propagator or state would launder an + GKSL dynamics. Returning a non-finite propagated state would launder an arithmetic failure into entropy, fitting and uncertainty calculations, so the relaxation layer fails closed. """ -def _evolve(L_super: np.ndarray, rho0: np.ndarray, t_grid: np.ndarray) -> np.ndarray: - """Propagate ``rho(t) = expm(L t) rho_0``, failing closed on non-finite work. +def _operator_is_finite(operator: Any) -> bool: + """Exact stored-entry finiteness check without densifying sparse matrices.""" + if sp.issparse(operator): + data = np.asarray(operator.data) + return bool(np.all(np.isfinite(data))) + return bool(np.all(np.isfinite(np.asarray(operator)))) - A finite generator and finite time do not imply that the dimensionless - product ``L*t``, the dense matrix exponential, or its action on the state - is representable in float64. Each boundary is checked separately so later - diagnostics cannot consume NaN/inf as if it were physical data. + +def _scaled_action_operand( + L_super: Any, t: float +) -> Any: + """Return t*L after representability checks needed by expm_multiply. + + The action algorithm needs norm information internally. A matrix whose + represented entries are finite can still have an unrepresentable 1-norm; + passing that case onward can turn a numerical-domain failure into an + effectively unbounded scaling loop. Refuse only when the matrix or its + mathematical 1-norm is not representable -- no physics threshold is used. + """ + with np.errstate(over="ignore", invalid="ignore", under="ignore"): + scaled = L_super * t + if not _operator_is_finite(scaled): + raise UnrepresentableTrajectoryError( + "relaxation trajectory: L*t contains non-finite entries at " + f"t={float(t):.6g}" + ) + + try: + with warnings.catch_warnings(): + warnings.filterwarnings("error", category=RuntimeWarning) + with np.errstate(over="ignore", invalid="ignore", under="ignore"): + norm_1 = ( + float(spla.norm(scaled, ord=1)) + if sp.issparse(scaled) + else float(np.linalg.norm(np.asarray(scaled), ord=1)) + ) + except (RuntimeWarning, OverflowError, ValueError) as exc: + raise UnrepresentableTrajectoryError( + "relaxation trajectory: ||L*t||_1 is not representable in float64 " + f"at t={float(t):.6g}" + ) from exc + if not np.isfinite(norm_1): + raise UnrepresentableTrajectoryError( + "relaxation trajectory: ||L*t||_1 is not representable in float64 " + f"at t={float(t):.6g}" + ) + return scaled + + +def _expm_action(operator: Any, state: np.ndarray, *, t: float) -> np.ndarray: + """Compute one exponential action and normalise numerical failures.""" + scaled = _scaled_action_operand(operator, t) + try: + with warnings.catch_warnings(): + warnings.filterwarnings("error", category=RuntimeWarning) + with np.errstate(over="ignore", invalid="ignore", under="ignore"): + out = spla.expm_multiply(scaled, state) + except ( + RuntimeWarning, + OverflowError, + ValueError, + np.linalg.LinAlgError, + sla.LinAlgError, + ) as exc: + raise UnrepresentableTrajectoryError( + "relaxation trajectory: scipy.sparse.linalg.expm_multiply could not " + f"represent the exponential action at t={float(t):.6g}" + ) from exc + out = np.asarray(out, dtype=complex) + if not np.all(np.isfinite(out)): + raise UnrepresentableTrajectoryError( + "relaxation trajectory: scipy.sparse.linalg.expm_multiply returned " + f"a non-finite state at t={float(t):.6g}" + ) + return out + + +def _is_exact_linspace(t_grid: np.ndarray) -> bool: + """Whether t_grid is exactly reproducible by NumPy linspace.""" + if t_grid.ndim != 1 or t_grid.size < 2: + return False + expected = np.linspace( + float(t_grid[0]), float(t_grid[-1]), int(t_grid.size), endpoint=True + ) + return bool(np.array_equal(np.asarray(t_grid, dtype=float), expected)) + + +def _evolve_with_backend( + L_super: Any, rho0: np.ndarray, t_grid: np.ndarray +) -> tuple[np.ndarray, str]: + """Propagate by exponential action without materialising exp(tL). + + Exactly linspace-generated grids use SciPy interval mode so setup work can + be reused. Arbitrary grids use independent actions from the same initial + state; this avoids accumulated stepwise error on non-uniform grids. """ rho_vec0 = vec(rho0) d = rho0.shape[0] - traj = np.empty((t_grid.size, d, d), dtype=complex) - for k, t in enumerate(t_grid): - if t == 0.0: - traj[k] = rho0 - continue + times = np.asarray(t_grid, dtype=float) + traj = np.empty((times.size, d, d), dtype=complex) - with np.errstate(over="ignore", invalid="ignore", under="ignore"): - scaled = np.asarray(L_super * t) - if not np.all(np.isfinite(scaled)): - raise UnrepresentableTrajectoryError( - "relaxation trajectory: L*t contains non-finite entries at " - f"t={float(t):.6g}; the requested dimensionless propagation " - "is outside the current float64 dense-expm domain" - ) + if times.size == 0: + return traj, "expm_multiply_pointwise" + if _is_exact_linspace(times) and times.size >= 2: + max_t = float(np.max(np.abs(times))) + _scaled_action_operand(L_super, max_t) try: - # SciPy's scaling-and-squaring path can emit RuntimeWarning outside - # NumPy's errstate. Convert that into the same explicit domain - # failure so warnings-as-errors and ordinary callers agree. with warnings.catch_warnings(): warnings.filterwarnings("error", category=RuntimeWarning) with np.errstate(over="ignore", invalid="ignore", under="ignore"): - propagator = sla.expm(scaled) + states = spla.expm_multiply( + L_super, + rho_vec0, + start=float(times[0]), + stop=float(times[-1]), + num=int(times.size), + endpoint=True, + ) except ( RuntimeWarning, OverflowError, @@ -83,26 +172,58 @@ def _evolve(L_super: np.ndarray, rho0: np.ndarray, t_grid: np.ndarray) -> np.nda sla.LinAlgError, ) as exc: raise UnrepresentableTrajectoryError( - "relaxation trajectory: scipy.linalg.expm could not represent " - f"a finite propagator at t={float(t):.6g}; use a different " - "numerical propagation method or time window" + "relaxation trajectory: scipy.sparse.linalg.expm_multiply " + "interval propagation failed" ) from exc - - if not np.all(np.isfinite(propagator)): + states = np.asarray(states, dtype=complex) + if states.shape != (times.size, rho_vec0.size) or not np.all( + np.isfinite(states) + ): raise UnrepresentableTrajectoryError( - "relaxation trajectory: scipy.linalg.expm returned a non-finite " - f"propagator at t={float(t):.6g} although L*t is finite" + "relaxation trajectory: expm_multiply interval propagation " + "returned a non-finite or malformed state sequence" ) + for k in range(times.size): + traj[k] = unvec(states[k], d=d) + if times[0] == 0.0: + traj[0] = rho0 + return traj, "expm_multiply_interval" - with np.errstate(over="ignore", invalid="ignore", under="ignore"): - rho_vec_t = propagator @ rho_vec0 - if not np.all(np.isfinite(rho_vec_t)): - raise UnrepresentableTrajectoryError( - "relaxation trajectory: the propagated state became non-finite " - f"at t={float(t):.6g}" - ) - traj[k] = unvec(rho_vec_t, d=d) - return traj + for k, t in enumerate(times): + if t == 0.0: + traj[k] = rho0 + continue + traj[k] = unvec(_expm_action(L_super, rho_vec0, t=float(t)), d=d) + return traj, "expm_multiply_pointwise" + + +def _evolve(L_super: Any, rho0: np.ndarray, t_grid: np.ndarray) -> np.ndarray: + """Compatibility wrapper returning only the propagated trajectory.""" + return _evolve_with_backend(L_super, rho0, t_grid)[0] + + +def _trajectory_audit(traj: np.ndarray) -> tuple[float, float, float]: + """Return trace error, Hermiticity defect, and minimum Hermitian eigenvalue. + + These are measurements, not positivity repairs or gates. In particular a + small negative eigenvalue is recorded as drift rather than silently clipped. + """ + max_trace_error = 0.0 + max_hermiticity_defect = 0.0 + min_eigenvalue = float("inf") + for rho in np.asarray(traj): + max_trace_error = max( + max_trace_error, float(abs(np.trace(rho) - 1.0)) + ) + anti = rho - rho.conj().T + max_hermiticity_defect = max( + max_hermiticity_defect, float(np.linalg.norm(anti, ord="fro")) + ) + herm = 0.5 * (rho + rho.conj().T) + min_eigenvalue = min( + min_eigenvalue, float(np.min(np.linalg.eigvalsh(herm)).real) + ) + return max_trace_error, max_hermiticity_defect, min_eigenvalue def von_neumann_entropy(rho: np.ndarray) -> float: @@ -387,7 +508,12 @@ def compute_relaxation_layer( if t_grid is None: t_grid = np.linspace(0.0, 10.0, 80) - traj = _evolve(L_super, rho_initial, t_grid) + traj, trajectory_backend = _evolve_with_backend(L_super, rho_initial, t_grid) + ( + trajectory_max_trace_error, + trajectory_max_hermiticity_defect, + trajectory_min_eigenvalue, + ) = _trajectory_audit(traj) final_rho = traj[-1] rel_entropy = np.array( [relative_entropy(traj[k], rho_steady_state) for k in range(traj.shape[0])] @@ -488,4 +614,8 @@ def compute_relaxation_layer( bca_ci_beta=(bca_lo, bca_hi), beta_D_linear=float(beta_D_linear), linear_fit_model=linear_fit_model, + trajectory_backend=trajectory_backend, + trajectory_max_trace_error=trajectory_max_trace_error, + trajectory_max_hermiticity_defect=trajectory_max_hermiticity_defect, + trajectory_min_eigenvalue=trajectory_min_eigenvalue, ) diff --git a/tests/test_issue162_expm_multiply.py b/tests/test_issue162_expm_multiply.py new file mode 100644 index 00000000..884b6b28 --- /dev/null +++ b/tests/test_issue162_expm_multiply.py @@ -0,0 +1,115 @@ +"""Issue #162: exponential-action trajectory backend.""" + +from __future__ import annotations + +import numpy as np +import pytest +import scipy.linalg as sla +import scipy.sparse as sp + +from liouscope.core.lindblad import build_liouvillian +from liouscope.diagnostics.relaxation import ( + _evolve, + _evolve_with_backend, + _trajectory_audit, + compute_relaxation_layer, +) +from liouscope.numerics.kronecker import unvec, vec + +_SM = np.array([[0.0, 1.0], [0.0, 0.0]], dtype=complex) +_RHO_PLUS = 0.5 * np.array([[1.0, 1.0], [1.0, 1.0]], dtype=complex) + + +def _generator(scale: float = 1.0) -> np.ndarray: + return build_liouvillian( + np.zeros((2, 2), dtype=complex), + [np.sqrt(0.7 * scale) * _SM], + ) + + +def _dense_reference(L: np.ndarray, grid: np.ndarray) -> np.ndarray: + state = vec(_RHO_PLUS) + out = np.empty((grid.size, 2, 2), dtype=complex) + for k, t in enumerate(grid): + if t == 0.0: + out[k] = _RHO_PLUS + else: + out[k] = unvec(sla.expm(L * t) @ state, d=2) + return out + + +@pytest.mark.parametrize( + "grid", + [ + np.linspace(0.0, 4.0, 17), + np.array([0.0, 0.01, 0.1, 0.7, 2.0, 4.0]), + ], +) +def test_action_backend_agrees_with_dense_reference(grid: np.ndarray) -> None: + """Action propagation matches the former full-exponential formula.""" + L = _generator() + actual = _evolve(L, _RHO_PLUS, grid) + expected = _dense_reference(L, grid) + np.testing.assert_allclose(actual, expected, rtol=2.0e-13, atol=2.0e-14) + + +def test_exact_linspace_uses_interval_backend_and_preserves_t0_bitwise() -> None: + grid = np.linspace(0.0, 3.0, 31) + traj, backend = _evolve_with_backend(_generator(), _RHO_PLUS, grid) + assert backend == "expm_multiply_interval" + np.testing.assert_array_equal(traj[0], _RHO_PLUS) + + +def test_nonuniform_grid_uses_pointwise_backend() -> None: + grid = np.array([0.0, 0.01, 0.13, 0.7, 3.0]) + _traj, backend = _evolve_with_backend(_generator(), _RHO_PLUS, grid) + assert backend == "expm_multiply_pointwise" + + +def test_sparse_and_dense_action_paths_agree() -> None: + L = _generator() + grid = np.linspace(0.0, 2.0, 21) + dense = _evolve(L, _RHO_PLUS, grid) + sparse = _evolve(sp.csr_matrix(L), _RHO_PLUS, grid) + np.testing.assert_allclose(sparse, dense, rtol=3.0e-13, atol=3.0e-14) + + +def test_rate_unit_metamorphism_L_to_cL_t_to_t_over_c() -> None: + """Pure rate-unit changes must leave the propagated states invariant.""" + L = _generator() + grid = np.array([0.0, 0.02, 0.11, 0.8, 3.0]) + reference = _evolve(L, _RHO_PLUS, grid) + for c in (1.0e-6, 1.0e6): + changed = _evolve(c * L, _RHO_PLUS, grid / c) + np.testing.assert_allclose( + changed, reference, rtol=5.0e-12, atol=5.0e-13 + ) + + +def test_trajectory_audit_measures_physicality_drift_without_repairing() -> None: + traj = _evolve(_generator(), _RHO_PLUS, np.linspace(0.0, 3.0, 17)) + trace_error, hermiticity, min_eval = _trajectory_audit(traj) + assert trace_error < 1.0e-12 + assert hermiticity < 1.0e-12 + assert min_eval > -1.0e-12 + + perturbed = traj.copy() + perturbed[-1] = np.diag([-1.0e-5, 1.0 + 1.0e-5]).astype(complex) + _te, _hd, drift = _trajectory_audit(perturbed) + assert drift == pytest.approx(-1.0e-5) + # The audit records the input; it does not clip it back to positivity. + assert perturbed[-1, 0, 0] == pytest.approx(-1.0e-5) + + +def test_relaxation_result_records_backend_and_audit_fields() -> None: + report = compute_relaxation_layer( + _generator(), + rho_initial=_RHO_PLUS, + t_grid=np.linspace(0.0, 2.0, 16), + bootstrap_B=5, + seed=1, + ) + assert report.trajectory_backend == "expm_multiply_interval" + assert report.trajectory_max_trace_error < 1.0e-12 + assert report.trajectory_max_hermiticity_defect < 1.0e-12 + assert report.trajectory_min_eigenvalue > -1.0e-12 diff --git a/tests/test_trajectory_representability.py b/tests/test_trajectory_representability.py index 8f57b409..564320c8 100644 --- a/tests/test_trajectory_representability.py +++ b/tests/test_trajectory_representability.py @@ -41,8 +41,8 @@ def _extreme_thermalising_qubit() -> np.ndarray: return L -def test_ordinary_dense_trajectory_is_bitwise_the_existing_formula() -> None: - """The guard must not change representable dense propagation.""" +def test_ordinary_action_trajectory_agrees_with_dense_reference() -> None: + """The backend change is bounded against the former dense formula.""" L = _ordinary_generator() grid = np.array([0.0, 0.2, 1.0, 2.0]) actual = _evolve(L, _RHO_PLUS, grid) @@ -54,7 +54,7 @@ def test_ordinary_dense_trajectory_is_bitwise_the_existing_formula() -> None: expected[k] = _RHO_PLUS else: expected[k] = unvec(sla.expm(L * t) @ initial, d=2) - np.testing.assert_array_equal(actual, expected) + np.testing.assert_allclose(actual, expected, rtol=2.0e-13, atol=2.0e-14) def test_nonfinite_scaled_generator_fails_before_matrix_exponential() -> None: @@ -68,15 +68,15 @@ def test_nonfinite_scaled_generator_fails_before_matrix_exponential() -> None: _evolve(L, _RHO_PLUS, grid) -def test_extreme_finite_Lt_fails_closed_when_dense_expm_is_unusable() -> None: - """Finite L*t is necessary, not sufficient, for a usable dense propagator.""" +def test_extreme_finite_Lt_stays_fail_closed_under_action_backend() -> None: + """The action backend must not turn the #156 extreme into a hang or NaN.""" L = _extreme_thermalising_qubit() grid = np.array([0.0, 1.1]) assert np.all(np.isfinite(L * grid[-1])) with pytest.raises( UnrepresentableTrajectoryError, - match=r"scipy\.linalg\.expm", + match=r"\|\|L\*t\|\|_1|expm_multiply", ): _evolve(L, _RHO_PLUS, grid) @@ -86,24 +86,29 @@ def test_scipy_runtime_warning_is_normalised_to_the_domain_error( ) -> None: """Warnings-as-errors and ordinary callers must see the same contract.""" - def _warn(_scaled: np.ndarray) -> np.ndarray: - raise RuntimeWarning("synthetic scaling-and-squaring failure") + def _warn(*args: object, **kwargs: object) -> np.ndarray: + raise RuntimeWarning("synthetic exponential-action failure") - monkeypatch.setattr(relaxation.sla, "expm", _warn) - with pytest.raises(UnrepresentableTrajectoryError, match=r"scipy\.linalg\.expm"): + monkeypatch.setattr(relaxation.spla, "expm_multiply", _warn) + with pytest.raises(UnrepresentableTrajectoryError, match=r"expm_multiply"): _evolve(_ordinary_generator(), _RHO_PLUS, np.array([0.0, 1.0])) -def test_nonfinite_action_is_rejected_even_with_finite_propagator( +def test_nonfinite_action_result_is_rejected( monkeypatch: pytest.MonkeyPatch, ) -> None: - """A finite propagator can still overflow when applied to the state.""" - huge = np.full((4, 4), 1.0e308, dtype=complex) - assert np.all(np.isfinite(huge)) - monkeypatch.setattr(relaxation.sla, "expm", lambda _scaled: huge) + """The action backend must not pass NaN/inf states downstream.""" + def _nonfinite(*args: object, **kwargs: object) -> np.ndarray: + return np.full(4, np.inf, dtype=complex) + + monkeypatch.setattr(relaxation.spla, "expm_multiply", _nonfinite) with pytest.raises( UnrepresentableTrajectoryError, - match=r"propagated state became non-finite", + match=r"non-finite", ): - _evolve(_ordinary_generator(), _RHO_PLUS, np.array([0.0, 1.0])) + _evolve( + _ordinary_generator(), + _RHO_PLUS, + np.array([0.0, 0.3, 1.0]), + )