From e51a327b46a4bafcd137b35086813431064b939c Mon Sep 17 00:00:00 2001 From: tjkessler Date: Thu, 15 Oct 2026 15:57:46 -0400 Subject: [PATCH 1/5] Reject non-stationary Newton stalls in reaction equilibrium equilibrium() clamps Newton iterates to c >= 0. On the converted Cantera h2o2 system at 1200 K it stopped at a near-boundary point (H, O ~ 0) where every flux is tiny and the absolute residual (~5e-11) looked converged, but net production equalled gross turnover: flux ran one way, so it was not a steady state. Cantera's thermochemical equilibrium has net/gross ~ 4e-15 there. After convergence, require max_s |dc_s/dt| / gross turnover_s <= 1e-6 and raise RuntimeError otherwise (suggesting a better guess). Exact boundary equilibria (e.g. Michaelis-Menten, zero flux) still pass. Co-Authored-By: Claude Opus 5.5 --- .../datasets/reaction/oracles.py | 34 ++++++++++++++++++- tests/datasets/test_reaction_oracles.py | 21 ++++++++++++ 2 files changed, 54 insertions(+), 1 deletion(-) diff --git a/src/koopman_graph/datasets/reaction/oracles.py b/src/koopman_graph/datasets/reaction/oracles.py index 58c8044..b6ab905 100644 --- a/src/koopman_graph/datasets/reaction/oracles.py +++ b/src/koopman_graph/datasets/reaction/oracles.py @@ -194,7 +194,12 @@ def equilibrium( accumulating into a species with no outflow), or arguments are invalid. RuntimeError - If Newton does not converge within ``max_iter`` iterations. + If Newton does not converge within ``max_iter`` iterations, or + stalls at a non-negative point that is not stationary. Clamping can + stop Newton on the boundary where every flux is tiny but flows one + way; pass a ``guess`` closer to the interior steady state, or use a + thermochemical equilibrium (e.g. Cantera's ``equilibrate``) as the + reference. """ c0 = _state(system, initial, "initial") if not math.isfinite(float(tol)) or float(tol) <= 0.0: @@ -242,6 +247,7 @@ def residual(c: Tensor) -> Tensor: moved = float((candidate - c).abs().max()) c, f = candidate, f_candidate if moved <= tol * (1.0 + float(c.abs().max())): + _require_stationary(system, c) return c msg = ( f"Newton did not converge for {system.name!r} in {max_iter} iterations " @@ -307,6 +313,32 @@ def oracle_timescales(system: ReactionSystem, state: Tensor) -> Tensor: return oracle_spectrum(system, state).timescales +_STATIONARY_RTOL = 1e-6 + + +def _require_stationary(system: ReactionSystem, c: Tensor) -> None: + """Raise unless net rates are negligible against gross species turnover. + + Parameters + ---------- + system + See signature. + c + See signature.""" + rates = system.reaction_rates(c) + net = system.time_derivative(c).abs() + gross = rates @ (system.reactant_stoichiometry + system.product_stoichiometry) + tiny = torch.finfo(torch.float64).tiny + ratio = float((net / gross.clamp(min=tiny)).max()) + if ratio > _STATIONARY_RTOL: + msg = ( + f"Newton stalled at a non-stationary point for {system.name!r} " + f"(max net/gross species turnover {ratio:.3e}); pass a guess " + "closer to the steady state" + ) + raise RuntimeError(msg) + + def _linear_equilibrium(system: ReactionSystem, c0: Tensor, w: Tensor) -> Tensor: """Solve ``A c + b = 0`` with ``W c = W c0`` for a linear system. diff --git a/tests/datasets/test_reaction_oracles.py b/tests/datasets/test_reaction_oracles.py index 16a89ad..9e638d0 100644 --- a/tests/datasets/test_reaction_oracles.py +++ b/tests/datasets/test_reaction_oracles.py @@ -185,3 +185,24 @@ def test_equilibrium_validation() -> None: for dt in (0.0, -1.0, math.inf): with pytest.raises(ValueError, match="time_step must be"): spectrum.discrete_eigenvalues(dt) + + +def test_equilibrium_rejects_one_way_boundary_point() -> None: + """A clamped Newton stall where flux runs one way is not a steady state.""" + # A -> B with a slow B -> A: from (1, 0) the only steady state is interior. + system = ReactionSystem( + species=("A", "B", "C"), + reactant_stoichiometry=torch.tensor( + [[1.0, 1.0, 0.0], [0.0, 0.0, 1.0]], dtype=_F64 + ), + product_stoichiometry=torch.tensor( + [[0.0, 0.0, 1.0], [1.0, 1.0, 0.0]], dtype=_F64 + ), + rate_constants=torch.tensor([1.0, 1e-3], dtype=_F64), + ) + c_eq = equilibrium(system, torch.tensor([1.0, 1.0, 0.0], dtype=_F64)) + assert float(system.time_derivative(c_eq).abs().max()) < 1e-12 + from koopman_graph.datasets.reaction.oracles import _require_stationary + + with pytest.raises(RuntimeError, match="non-stationary point"): + _require_stationary(system, torch.tensor([1.0, 0.0, 1e-9], dtype=_F64)) From e1008d5deccfce8ac84c9f8f93f25d53d85a07fe Mon Sep 17 00:00:00 2001 From: tjkessler Date: Thu, 15 Oct 2026 16:09:02 -0400 Subject: [PATCH 2/5] Add ReactionNetworkBenchmark with packaged FAIR cards MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Refs #65. Fills in koopman_graph.datasets.reaction.benchmark: - Six FAIR JSON cards under reaction/data/ define the frozen datasets (chain, Michaelis-Menten, Brusselator, random network, Van de Vusse CSTR with controls, Cantera h2o2). Cards record generator, simulator, grid, seed, initial ranges, log floor, graph encoding, oracle references, units, license, limitations, and Cantera provenance. - ReactionNetworkBenchmark.generate -> ReactionDataset (float64 concentrations, controls, per-trajectory GraphSnapshotSequence with species nodes, timestamps, control_inputs), optional SSA simulator, deterministic trajectory split, versioned save/load (weights_only). - compute_card_oracle + scripts/build_reaction_cards.py (--check) write the oracle sections; tests recompute them. - eigenvalue_timescales maps numerically zero Re(λ) to inf; the stationarity check gains a flux-scaled absolute tolerance so boundary steady states with vanishing leftovers are accepted. - Exported from koopman_graph.datasets (root __all__ unchanged). Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 46 +- docs/source/api.rst | 1 + scripts/build_reaction_cards.py | 290 ++++++ src/koopman_graph/datasets/__init__.py | 8 + .../datasets/reaction/__init__.py | 20 + .../datasets/reaction/benchmark.py | 974 ++++++++++++++++++ .../datasets/reaction/data/.gitkeep | 0 .../data/reaction_brusselator_v1.json | 67 ++ .../data/reaction_cantera_h2o2_v1.json | 167 +++ .../reaction/data/reaction_chain_v1.json | 95 ++ .../data/reaction_michaelis_menten_v1.json | 98 ++ .../data/reaction_random_network_v1.json | 131 +++ .../data/reaction_van_de_vusse_v1.json | 108 ++ .../datasets/reaction/oracles.py | 82 +- tests/datasets/test_reaction_benchmark.py | 426 ++++++++ .../test_reaction_import_isolation.py | 9 + tests/datasets/test_reaction_oracles.py | 14 + 17 files changed, 2513 insertions(+), 23 deletions(-) create mode 100644 scripts/build_reaction_cards.py delete mode 100644 src/koopman_graph/datasets/reaction/data/.gitkeep create mode 100644 src/koopman_graph/datasets/reaction/data/reaction_brusselator_v1.json create mode 100644 src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json create mode 100644 src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json create mode 100644 src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json create mode 100644 src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json create mode 100644 src/koopman_graph/datasets/reaction/data/reaction_van_de_vusse_v1.json create mode 100644 tests/datasets/test_reaction_benchmark.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 5e0d63d..5d1be36 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -80,14 +80,28 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ideal-gas reactors, adiabatic or isothermal. - Also: `require_cantera`, `resolve_mechanism_path`, and `mechanism_sha256`. -- `[reaction]` extra (`scipy>=1.11.0`) for stiff reaction-network solvers - and `[cantera]` extra (`cantera>=3.2,<3.3`) for Cantera mechanisms. -- `cantera` pytest marker; marked tests are skipped automatically when the - `[cantera]` extra is not installed. -- CI test jobs and the release verify job install `[reaction]` and - `[cantera]`, so `cantera`-marked tests run there. -- `tests/test_import_hygiene.py`: `import koopman_graph` must not load - SciPy, Cantera, or torchdiffeq. +- `ReactionNetworkBenchmark` (exported from `koopman_graph.datasets` and + `koopman_graph.datasets.reaction`): card-defined, seeded reaction-network + datasets. + - Six packaged FAIR JSON cards under `reaction/data/`: + `reaction_chain_v1`, `reaction_michaelis_menten_v1`, + `reaction_brusselator_v1`, `reaction_random_network_v1`, + `reaction_van_de_vusse_v1` (piecewise-constant temperature / + dilution-rate controls), and `reaction_cantera_h2o2_v1` (`[cantera]`). + Each card records the generator, simulator, time grid, seed, + initial-condition ranges, log floor, graph encoding, oracle references + (equilibrium, dynamic eigenvalues, timescales, conservation laws), + units, license, limitations, and Cantera provenance. + - `generate` returns a `ReactionDataset` (raw `float64` concentrations, + optional controls, per-trajectory `GraphSnapshotSequence` with species + nodes, `timestamps`, and global `control_inputs`), with an optional + Gillespie `simulator="ssa"`. Also: deterministic trajectory `split`, and + versioned `save` / `load` payloads (`weights_only` loading). + - `compute_card_oracle` and `scripts/build_reaction_cards.py` (with + `--check`) regenerate the oracle sections; tests recompute them to catch + drift. +- `eigenvalue_timescales`: `-1/Re(λ)` with numerically zero real parts + (relative to `max|λ|`) mapped to `inf`; `OracleSpectrum.timescales` uses it. ### Changed @@ -98,6 +112,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `uv.lock` relocked; it was missing extras already declared in `pyproject.toml` (`export`, `control-lmi`, ...). +### Fixed + +- `equilibrium` no longer returns a clamped Newton stall that is not + stationary. It now requires `|dc_s/dt| <= 1e-6 · gross turnover_s + + 1e-12 · (peak initial turnover)` and otherwise raises `RuntimeError`. + Previously the converted Cantera `h2o2` system could report a + near-boundary point with one-way flux as its equilibrium. +- `[reaction]` extra (`scipy>=1.11.0`) for stiff reaction-network solvers + and `[cantera]` extra (`cantera>=3.2,<3.3`) for Cantera mechanisms. +- `cantera` pytest marker; marked tests are skipped automatically when the + `[cantera]` extra is not installed. +- CI test jobs and the release verify job install `[reaction]` and + `[cantera]`, so `cantera`-marked tests run there. +- `tests/test_import_hygiene.py`: `import koopman_graph` must not load + SciPy, Cantera, or torchdiffeq. + ## [0.15.0] - 2026-08-19 Integration increment on the 0.14 stack: opt-in closed-form diff --git a/docs/source/api.rst b/docs/source/api.rst index b378b0c..360013b 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -1046,6 +1046,7 @@ Datasets .. automodule:: koopman_graph.datasets.reaction :members: :show-inheritance: + :exclude-members: ReactionNetworkBenchmark .. automodule:: koopman_graph.datasets.era5 :members: diff --git a/scripts/build_reaction_cards.py b/scripts/build_reaction_cards.py new file mode 100644 index 0000000..567e816 --- /dev/null +++ b/scripts/build_reaction_cards.py @@ -0,0 +1,290 @@ +#!/usr/bin/env python3 +"""Write the FAIR JSON cards for ``koopman_graph.datasets.reaction``. + +Maintainer tool. The specs below define every frozen ``reaction_*_v1`` +dataset. Running this script computes each card's ``oracle`` section with +:func:`koopman_graph.datasets.reaction.benchmark.compute_card_oracle` and +writes ``src/koopman_graph/datasets/reaction/data/.json``. +The committed JSON is the source of truth at runtime; +``tests/datasets/test_reaction_benchmark.py`` recomputes the oracles to +catch drift. The Cantera card needs the ``[cantera]`` extra (Cantera 3.2). + +Usage:: + + python scripts/build_reaction_cards.py # write all cards + python scripts/build_reaction_cards.py --check # exit 1 if any differ +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path +from typing import Any + +from koopman_graph.datasets.reaction.benchmark import ( + compute_card_oracle, + system_from_card, +) + +DATA_DIR = ( + Path(__file__).resolve().parents[1] + / "src" + / "koopman_graph" + / "datasets" + / "reaction" + / "data" +) + +_LICENSE = "Apache-2.0" +# SHA-256 of h2o2.yaml in the Cantera 3.2.0 wheels (issue #65). +_H2O2_SHA256 = "0efc6c52862741a29e0c29b65d979c7d8cb409db5282bca83b9c5437b3d8c8d4" +_UNITS = {"time": "s", "concentration": "mol/L"} +_GENERATED = "In-repo mass-action generator (koopman_graph.datasets.reaction)." + +SPECS: list[dict[str, Any]] = [ + { + "dataset_id": "reaction_chain_v1", + "scope": "Linear first-order chain S0 -> S1 -> S2 -> S3 (rates 1, " + "0.5, 0.25); exact matrix-exponential trajectories with closed-form " + "dynamic spectrum {-0.25, -0.5, -1} plus one conservation zero mode.", + "generator": { + "kind": "catalog", + "name": "first_order_chain", + "params": {"num_species": 4}, + }, + "simulator": {"kind": "expm"}, + "time_step": 0.25, + "default_num_timesteps": 64, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "S0": [0.5, 1.5], + "S1": [0.0, 0.2], + "S2": [0.0, 0.2], + "S3": [0.0, 0.2], + }, + "transform": {"kind": "log_concentration", "floor": 1e-8}, + "graph_encoding": "species_influence", + "limitations": "Linear and closed; every mode is a pure decay. Use as " + "the exact-spectrum sanity check, not as a nonlinear benchmark.", + }, + { + "dataset_id": "reaction_michaelis_menten_v1", + "scope": "Mass-action Michaelis-Menten E + S <-> ES -> E + P with two " + "conservation laws (enzyme, substrate backbone) and a boundary " + "equilibrium.", + "generator": {"kind": "catalog", "name": "michaelis_menten", "params": {}}, + "simulator": {"kind": "rk4", "substeps": 50}, + "time_step": 0.5, + "default_num_timesteps": 64, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "E": [0.05, 0.15], + "S": [0.5, 1.5], + "ES": [0.0, 0.05], + "P": [0.0, 0.1], + }, + "transform": {"kind": "log_concentration", "floor": 1e-8}, + "graph_encoding": "reaction_hyperedge", + "limitations": "Irreversible catalysis: trajectories end on the " + "boundary S = ES = 0, where the log transform saturates at the floor.", + }, + { + "dataset_id": "reaction_brusselator_v1", + "scope": "Open Brusselator (a = 1, b = 3) with an unstable focus and a " + "limit cycle; tests oscillatory, non-decaying dynamics.", + "generator": {"kind": "catalog", "name": "brusselator", "params": {}}, + "simulator": {"kind": "rk4", "substeps": 10}, + "time_step": 0.1, + "default_num_timesteps": 128, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": {"X": [0.5, 1.5], "Y": [2.0, 4.0]}, + "transform": None, + "graph_encoding": "species_influence", + "limitations": "Open system without conservation laws; the oracle " + "fixed point is unstable, so its timescales are negative (growth).", + }, + { + "dataset_id": "reaction_random_network_v1", + "scope": "Seeded random element-balanced mass-action network (6 " + "species, 2 elements, 8 reactions, order <= 2).", + "generator": { + "kind": "random_network", + "params": { + "num_species": 6, + "num_elements": 2, + "num_reactions": 8, + "max_atoms": 4, + "max_order": 2, + "rate_log10_range": [-1.0, 1.0], + "seed": 0, + }, + }, + "simulator": {"kind": "rk4", "substeps": 20}, + "time_step": 0.1, + "default_num_timesteps": 64, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": {f"X{i}": [0.1, 1.0] for i in range(6)}, + "oracle_relaxation_time": 200.0, + "transform": {"kind": "log_concentration", "floor": 1e-8}, + "graph_encoding": "species_influence", + "limitations": "Synthetic chemistry with no physical species; rates " + "are log-uniform draws, not measured kinetics.", + }, + { + "dataset_id": "reaction_van_de_vusse_v1", + "scope": "Van de Vusse CSTR (A -> B -> C, 2 A -> D) with seeded " + "piecewise-constant reactor temperature and dilution-rate inputs " + "stored as global control_inputs.", + "generator": {"kind": "cstr", "params": {}}, + "simulator": {"kind": "cstr", "substeps": 20}, + "time_step": 0.01, + "default_num_timesteps": 100, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "A": [0.5, 3.0], + "B": [0.0, 1.5], + "C": [0.0, 1.5], + "D": [0.0, 0.5], + }, + "controls": { + "names": ["temperature_k", "dilution_rate"], + "ranges": {"temperature_k": [380.0, 400.0], "dilution_rate": [5.0, 30.0]}, + "hold_steps": 10, + "units": {"temperature_k": "K", "dilution_rate": "1/h"}, + }, + "transform": None, + "graph_encoding": "species_influence", + "units": {"time": "h", "concentration": "mol/L"}, + "references": [ + "K.-U. Klatt and S. Engell, Comput. Chem. Eng. 22(4-5), 491-502, " + "1998 (kinetics)", + "H. Chen, A. Kremling and F. Allgower, Proc. 3rd ECC, 1995 " + "(operating point)", + ], + "limitations": "Reactor temperature is a direct input (no energy " + "balance); the graph uses the reaction network only, not the feed.", + }, + { + "dataset_id": "reaction_cantera_h2o2_v1", + "scope": "Isothermal constant-volume H2/O2/Ar chemistry at 1200 K from " + "Cantera's bundled h2o2.yaml (10 species, 29 reactions); " + "trajectories come from Cantera's ReactorNet.", + "generator": {"kind": "cantera", "params": {}}, + "simulator": {"kind": "cantera"}, + "time_step": 2e-5, + "default_num_timesteps": 50, + "default_num_trajectories": 16, + "default_seed": 0, + "cantera": { + "mechanism": "h2o2.yaml", + "sha256": _H2O2_SHA256, + "cantera_version": "3.2.0", + "temperature_k": 1200.0, + "pressure_pa": 101325.0, + "reactor": "const_volume", + "isothermal": True, + "state_dependent": "freeze", + "reference_mole_fractions": "H2:2, O2:1, AR:4", + "initial_mole_fraction_ranges": { + "H2": [1.0, 3.0], + "O2": [0.5, 1.5], + "AR": [3.0, 5.0], + }, + "note": "Cantera's bundled input files are provided for " + "illustration only; they are not validated kinetic models for " + "quantitative use.", + "citation_doi": "10.5281/zenodo.17620923", + "license": "Cantera is BSD-3-Clause", + }, + "transform": {"kind": "log_concentration", "floor": 1e-12}, + "graph_encoding": "species_influence", + "limitations": "Illustration-only mechanism; the graph and oracle use " + "the mass-action conversion with the falloff reaction frozen at the " + "reference state, while trajectories are exact Cantera integrations.", + }, +] + + +def build_card(spec: dict[str, Any]) -> dict[str, Any]: + """Return the full card for ``spec`` with species, provenance, and oracle. + + Parameters + ---------- + spec : dict + One entry of :data:`SPECS`. + + Returns + ------- + dict + Card ready to serialize. + """ + card = { + "dataset_id": spec["dataset_id"], + "version": spec["dataset_id"].rsplit("_", 1)[-1], + **{k: v for k, v in spec.items() if k != "dataset_id"}, + } + card.setdefault("units", dict(_UNITS)) + card.setdefault("controls", None) + card.setdefault("cantera", None) + card["format"] = ( + "Per-trajectory GraphSnapshotSequence (one node per species, float32 " + "features of shape (num_species, 1), timestamps, optional global " + "control_inputs); raw float64 concentrations in ReactionDataset." + ) + card["provenance"] = ( + _GENERATED + if card["cantera"] is None + else "Cantera ReactorNet integration of the bundled mechanism; see the " + "cantera section for version and SHA-256." + ) + card["license"] = _LICENSE + system = system_from_card(card) + card["species"] = list(system.species) + card["elements"] = None if system.elements is None else list(system.elements) + card["num_reactions"] = system.num_reactions + card["oracle"] = compute_card_oracle(card) + return card + + +def main(argv: list[str] | None = None) -> int: + """Write (or check) every card. + + Parameters + ---------- + argv : list of str or None, optional + Command-line arguments. + + Returns + ------- + int + Exit status (``1`` when ``--check`` finds a difference). + """ + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--check", action="store_true", help="compare, do not write") + args = parser.parse_args(argv) + DATA_DIR.mkdir(parents=True, exist_ok=True) + stale = [] + for spec in SPECS: + text = json.dumps(build_card(spec), indent=2) + "\n" + path = DATA_DIR / f"{spec['dataset_id']}.json" + if args.check: + if not path.is_file() or path.read_text(encoding="utf-8") != text: + stale.append(path.name) + else: + path.write_text(text, encoding="utf-8") + print(f"wrote {path.relative_to(DATA_DIR.parents[3])}") + if stale: + print("stale cards: " + ", ".join(stale), file=sys.stderr) + return 1 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/koopman_graph/datasets/__init__.py b/src/koopman_graph/datasets/__init__.py index 639748e..fbcd74c 100644 --- a/src/koopman_graph/datasets/__init__.py +++ b/src/koopman_graph/datasets/__init__.py @@ -69,6 +69,12 @@ ERA5SliceCard / generate_tiny_era5_mesh Documented ERA5 teaching-slice card and a tiny generated mesh (not production GraphCast / ERA5 training). +ReactionNetworkBenchmark + Card-defined mass-action reaction networks (catalog systems, random + networks, Van de Vusse CSTR, Cantera ``h2o2``) with oracle equilibria, + spectra, and timescales. The rest of the reaction toolkit lives in + :mod:`koopman_graph.datasets.reaction`; SciPy / Cantera stay behind the + ``[reaction]`` / ``[cantera]`` extras. """ from koopman_graph.datasets.contact_epidemic import ContactEpidemicBenchmark @@ -91,6 +97,7 @@ from koopman_graph.datasets.lorenz96 import Lorenz96GraphBenchmark from koopman_graph.datasets.metr_la import MetrLaTrafficBenchmark from koopman_graph.datasets.pems import PemsBayTrafficBenchmark, PemsTrafficBenchmark +from koopman_graph.datasets.reaction.benchmark import ReactionNetworkBenchmark from koopman_graph.datasets.synthetic import SyntheticDynamicGraphBenchmark from koopman_graph.datasets.topology import TopologyPayload @@ -107,6 +114,7 @@ "MetrLaTrafficBenchmark", "PemsBayTrafficBenchmark", "PemsTrafficBenchmark", + "ReactionNetworkBenchmark", "SyntheticDynamicGraphBenchmark", "TopologyPayload", "TypedIEEE118Topology", diff --git a/src/koopman_graph/datasets/reaction/__init__.py b/src/koopman_graph/datasets/reaction/__init__.py index 8ff840a..caa0a66 100644 --- a/src/koopman_graph/datasets/reaction/__init__.py +++ b/src/koopman_graph/datasets/reaction/__init__.py @@ -38,6 +38,16 @@ replacement for a production kinetics solver. """ +from koopman_graph.datasets.reaction.benchmark import ( + PAYLOAD_FORMAT, + PAYLOAD_FORMAT_VERSION, + GraphEncoding, + ReactionDataset, + ReactionDatasetSplit, + ReactionNetworkBenchmark, + compute_card_oracle, + system_from_card, +) from koopman_graph.datasets.reaction.cantera_io import ( CANTERA_CITATION_DOI, CANTERA_ILLUSTRATION_NOTE, @@ -102,6 +112,7 @@ from koopman_graph.datasets.reaction.oracles import ( OracleSpectrum, conservation_laws, + eigenvalue_timescales, equilibrium, oracle_spectrum, oracle_timescales, @@ -137,6 +148,8 @@ "DEFAULT_MECHANISM", "DEFAULT_MICHAELIS_MENTEN_RATES", "GAS_CONSTANT_J_PER_MOL_K", + "PAYLOAD_FORMAT", + "PAYLOAD_FORMAT_VERSION", "PRODUCT_RELATION", "REACTANT_RELATION", "REACTION_NODE", @@ -149,9 +162,13 @@ "CanteraProvenance", "CanteraReactorTrajectory", "ElementBalanceReport", + "GraphEncoding", "IntegrationMethod", "LogConcentrationTransform", "OracleSpectrum", + "ReactionDataset", + "ReactionDatasetSplit", + "ReactionNetworkBenchmark", "ReactionSystem", "SSATrajectory", "SampledReactionNetwork", @@ -160,7 +177,9 @@ "bin_to_grid", "brusselator", "check_element_balance", + "compute_card_oracle", "conservation_laws", + "eigenvalue_timescales", "element_balance_violation", "equilibrium", "first_order_chain", @@ -189,5 +208,6 @@ "species_influence_graph", "species_reaction_bipartite", "stoichiometric_rank", + "system_from_card", "van_de_vusse_cstr", ] diff --git a/src/koopman_graph/datasets/reaction/benchmark.py b/src/koopman_graph/datasets/reaction/benchmark.py index e69de29..9bee9f3 100644 --- a/src/koopman_graph/datasets/reaction/benchmark.py +++ b/src/koopman_graph/datasets/reaction/benchmark.py @@ -0,0 +1,974 @@ +"""``ReactionNetworkBenchmark``: card-defined reaction-network datasets. + +Each frozen dataset is defined by a FAIR JSON card under +``koopman_graph/datasets/reaction/data/`` (``.json``), read with +:mod:`importlib.resources`. The card is the source of truth. It records: + +* the generator (catalog constructor, random-network sampler, Van de Vusse + CSTR, or Cantera mechanism) and its parameters; +* the simulator, time grid, default trajectory count, and seed; +* per-species initial-concentration ranges (and control ranges for the + CSTR); +* the optional log-concentration floor and the graph encoding; +* oracle references (equilibrium, dynamic Jacobian eigenvalues, + timescales, conservation laws) that a learned operator can be checked + against; +* provenance, units, license, limitations, and Cantera provenance where + relevant. + +Oracle sections are computed by :func:`compute_card_oracle`; +``scripts/build_reaction_cards.py`` writes them, and the tests recompute +them so cards and code cannot drift apart. + +:meth:`ReactionNetworkBenchmark.generate` draws seeded initial conditions +(and control schedules), integrates every trajectory, and returns a +:class:`ReactionDataset`. Raw concentrations stay ``float64``; +:meth:`ReactionDataset.sequence` exposes each trajectory as a +:class:`~koopman_graph.data.GraphSnapshotSequence` with one node per +species, with ``timestamps`` and (for the CSTR) global ``control_inputs``. + +Importing this module does not import scipy or cantera; the Cantera card +needs the ``[cantera]`` extra at generation time. +""" + +from __future__ import annotations + +import json +import math +from dataclasses import dataclass +from importlib import resources +from pathlib import Path +from typing import Any, Literal + +import torch +from torch import Tensor +from torch_geometric.data import Data + +from koopman_graph.data import GraphSnapshotSequence +from koopman_graph.datasets.reaction.catalog import ( + brusselator, + first_order_chain, + lotka_volterra, + michaelis_menten, + reversible_dimerization, +) +from koopman_graph.datasets.reaction.cstr import van_de_vusse_cstr +from koopman_graph.datasets.reaction.gillespie import gillespie_ssa +from koopman_graph.datasets.reaction.graphs import ( + reaction_hyperedge_index, + species_influence_graph, +) +from koopman_graph.datasets.reaction.integrate import ( + integrate_reaction_system, + mass_action_jacobian, +) +from koopman_graph.datasets.reaction.oracles import ( + conservation_laws, + eigenvalue_timescales, + equilibrium, + oracle_spectrum, +) +from koopman_graph.datasets.reaction.random_network import sample_random_network +from koopman_graph.datasets.reaction.system import ReactionSystem +from koopman_graph.datasets.reaction.transform import LogConcentrationTransform + +PAYLOAD_FORMAT = "koopman_graph.reaction_dataset" +"""``format`` tag written into saved dataset payloads.""" + +PAYLOAD_FORMAT_VERSION = 1 +"""Current payload schema version accepted by :meth:`ReactionDataset.load`.""" + +GraphEncoding = Literal["species_influence", "reaction_hyperedge"] +"""Topology attached to each snapshot (see :mod:`.graphs`).""" + +_PACKAGE = "koopman_graph.datasets.reaction" +_CATALOG_FACTORIES = { + "first_order_chain": first_order_chain, + "reversible_dimerization": reversible_dimerization, + "michaelis_menten": michaelis_menten, + "brusselator": brusselator, + "lotka_volterra": lotka_volterra, +} + + +@dataclass(frozen=True, eq=False) +class ReactionDataset: + """Generated reaction-network trajectories plus the card that defines them. + + Attributes + ---------- + dataset_id : str + Card identifier (e.g. ``"reaction_chain_v1"``). + card : dict + Parsed FAIR card. + system : ReactionSystem + Mass-action system the data were generated from (for the CSTR, rate + constants are at the nominal temperature; for Cantera, the + converted mechanism at the card's reference state). + times : Tensor + ``float64`` sample times ``(num_timesteps,)`` in the card's time unit. + concentrations : Tensor + ``float64`` concentrations ``(num_trajectories, num_timesteps, + num_species)`` in the card's concentration unit. + controls : Tensor or None + ``float64`` inputs ``(num_trajectories, num_timesteps, num_controls)``; + row ``t`` drives ``t → t+1`` (the last row repeats the last + interval). ``None`` for uncontrolled datasets. + transform : LogConcentrationTransform or None + Feature transform applied by :meth:`sequence`. + seed : int + Seed of the generator that drew initial conditions and controls. + simulator : str + Simulator actually used (``"expm"``, ``"rk4"``, ``"ssa"``, + ``"cstr"``, or ``"cantera"``). + volume : float or None + SSA volume when ``simulator == "ssa"``. + """ + + dataset_id: str + card: dict[str, Any] + system: ReactionSystem + times: Tensor + concentrations: Tensor + controls: Tensor | None + transform: LogConcentrationTransform | None + seed: int + simulator: str + volume: float | None = None + + @property + def num_trajectories(self) -> int: + """Return the number of trajectories. + + Returns + ------- + int + Leading dimension of :attr:`concentrations`. + """ + return int(self.concentrations.shape[0]) + + @property + def graph_encoding(self) -> str: + """Return the card's graph encoding name. + + Returns + ------- + str + ``"species_influence"`` or ``"reaction_hyperedge"``. + """ + return str(self.card["graph_encoding"]) + + def features(self, index: int, *, dtype: torch.dtype = torch.float32) -> Tensor: + """Return node features ``(num_timesteps, num_species, 1)`` of one trajectory. + + Parameters + ---------- + index : int + Trajectory index. + dtype : torch.dtype, optional + Output dtype. Default is ``torch.float32``. + + Returns + ------- + Tensor + Transformed concentrations when :attr:`transform` is set, raw + concentrations otherwise. + """ + c = self.concentrations[index] + values = c if self.transform is None else self.transform.forward(c) + return values.unsqueeze(-1).to(dtype) + + def sequence( + self, index: int, *, dtype: torch.dtype = torch.float32 + ) -> GraphSnapshotSequence: + """Return one trajectory as a species-node snapshot sequence. + + Parameters + ---------- + index : int + Trajectory index in ``[0, num_trajectories)``. + dtype : torch.dtype, optional + Feature dtype. Default is ``torch.float32``. + + Returns + ------- + GraphSnapshotSequence + ``num_timesteps`` snapshots with ``x`` of shape + ``(num_species, 1)``, shared ``edge_index`` / ``edge_weight`` + from :func:`species_influence_graph` (plus ``hyperedge_index`` + for ``"reaction_hyperedge"``), ``timestamps``, and global + ``control_inputs`` when present. + + Raises + ------ + IndexError + If ``index`` is out of range. + """ + if not 0 <= index < self.num_trajectories: + msg = f"index {index} out of range for {self.num_trajectories} trajectories" + raise IndexError(msg) + edge_index, edge_weight = species_influence_graph(self.system) + hyperedge_index = ( + reaction_hyperedge_index(self.system) + if self.graph_encoding == "reaction_hyperedge" + else None + ) + features = self.features(index, dtype=dtype) + snapshots = [] + for x in features: + data = Data(x=x, edge_index=edge_index, edge_weight=edge_weight.to(dtype)) + if hyperedge_index is not None: + data.hyperedge_index = hyperedge_index + snapshots.append(data) + controls = None if self.controls is None else self.controls[index].to(dtype) + return GraphSnapshotSequence( + snapshots, + timestamps=self.times.to(dtype=torch.float64), + control_inputs=controls, + ) + + def subset(self, indices: list[int] | Tensor) -> ReactionDataset: + """Return a dataset restricted to the given trajectory indices. + + Parameters + ---------- + indices : list of int or Tensor + Trajectory indices, in the order to keep. + + Returns + ------- + ReactionDataset + Same card, system, and times; selected trajectories only. + """ + index = torch.as_tensor(indices, dtype=torch.long) + return ReactionDataset( + dataset_id=self.dataset_id, + card=self.card, + system=self.system, + times=self.times, + concentrations=self.concentrations[index], + controls=None if self.controls is None else self.controls[index], + transform=self.transform, + seed=self.seed, + simulator=self.simulator, + volume=self.volume, + ) + + def save(self, path: str | Path) -> Path: + """Write the dataset to a versioned ``torch.save`` payload. + + Parameters + ---------- + path : str or Path + Output file (parent directories are created). + + Returns + ------- + Path + The written path. + """ + target = Path(path) + target.parent.mkdir(parents=True, exist_ok=True) + payload = { + "format": PAYLOAD_FORMAT, + "format_version": PAYLOAD_FORMAT_VERSION, + "dataset_id": self.dataset_id, + "card_json": json.dumps(self.card), + "system_json": json.dumps(self.system.to_dict()), + "times": self.times.clone(), + "concentrations": self.concentrations.clone(), + "controls": None if self.controls is None else self.controls.clone(), + "transform_floor": None if self.transform is None else self.transform.floor, + "seed": self.seed, + "simulator": self.simulator, + "volume": self.volume, + } + torch.save(payload, target) + return target + + @classmethod + def load(cls, path: str | Path) -> ReactionDataset: + """Read a payload written by :meth:`save` (``weights_only`` loading). + + Parameters + ---------- + path : str or Path + Payload file. + + Returns + ------- + ReactionDataset + The saved dataset. + + Raises + ------ + ValueError + If the payload format or version is not recognised. + """ + payload = torch.load(Path(path), weights_only=True) + if not isinstance(payload, dict) or payload.get("format") != PAYLOAD_FORMAT: + msg = f"{path} is not a {PAYLOAD_FORMAT} payload" + raise ValueError(msg) + version = payload.get("format_version") + if version != PAYLOAD_FORMAT_VERSION: + msg = ( + f"unsupported reaction dataset format_version {version!r}; " + f"expected {PAYLOAD_FORMAT_VERSION}" + ) + raise ValueError(msg) + floor = payload["transform_floor"] + return cls( + dataset_id=str(payload["dataset_id"]), + card=json.loads(payload["card_json"]), + system=ReactionSystem.from_dict(json.loads(payload["system_json"])), + times=payload["times"], + concentrations=payload["concentrations"], + controls=payload["controls"], + transform=None if floor is None else LogConcentrationTransform(floor), + seed=int(payload["seed"]), + simulator=str(payload["simulator"]), + volume=payload["volume"], + ) + + +@dataclass(frozen=True, eq=False) +class ReactionDatasetSplit: + """Trajectory-level train / validation / test split. + + Attributes + ---------- + train : ReactionDataset + First trajectories by index. + val : ReactionDataset + Following trajectories. + test : ReactionDataset + Last trajectories. + """ + + train: ReactionDataset + val: ReactionDataset + test: ReactionDataset + + +class ReactionNetworkBenchmark: + """Card-defined, seeded mass-action reaction-network datasets. + + Notes + ----- + Dataset card (per ``dataset_id``; see the JSON for exact values): + + * **Scope:** diagnostic / teaching data for checking whether a learned + latent operator recovers slow kinetic timescales and respects + conservation; not a production kinetics benchmark + * **Format:** per-trajectory ``GraphSnapshotSequence`` (one node per + species, ``float32`` features), plus raw ``float64`` concentrations + and oracle references in the card + * **Source:** in-repo generators; the Cantera card reads ``h2o2.yaml`` + from the installed Cantera package (illustration-only input) + * **License:** Apache-2.0 for generated data; Cantera is BSD-3-Clause + * **Version:** encoded in ``dataset_id`` (``*_v1``) and + :data:`PAYLOAD_FORMAT_VERSION` + """ + + @staticmethod + def dataset_ids() -> tuple[str, ...]: + """List the packaged dataset cards. + + Returns + ------- + tuple of str + Sorted card identifiers. + """ + root = resources.files(_PACKAGE).joinpath("data") + return tuple( + sorted( + entry.name.removesuffix(".json") + for entry in root.iterdir() + if entry.name.startswith("reaction_") and entry.name.endswith(".json") + ) + ) + + @staticmethod + def load_card(dataset_id: str) -> dict[str, Any]: + """Load a packaged FAIR card. + + Parameters + ---------- + dataset_id : str + Card identifier from :meth:`dataset_ids`. + + Returns + ------- + dict + Parsed card. + + Raises + ------ + ValueError + If no card has this identifier. + """ + if dataset_id not in ReactionNetworkBenchmark.dataset_ids(): + known = ", ".join(ReactionNetworkBenchmark.dataset_ids()) + msg = f"unknown reaction dataset {dataset_id!r}; known: {known}" + raise ValueError(msg) + resource = resources.files(_PACKAGE).joinpath(f"data/{dataset_id}.json") + return json.loads(resource.read_text(encoding="utf-8")) + + @classmethod + def generate( + cls, + dataset_id: str = "reaction_chain_v1", + *, + num_trajectories: int | None = None, + num_timesteps: int | None = None, + seed: int | None = None, + simulator: Literal["card", "ssa"] = "card", + volume: float = 1000.0, + ) -> ReactionDataset: + """Generate the trajectories a card defines. + + Parameters + ---------- + dataset_id : str, optional + Card identifier. Default is ``"reaction_chain_v1"``. + num_trajectories : int or None, optional + Trajectory count. Default is the card's + ``default_num_trajectories``. + num_timesteps : int or None, optional + Samples per trajectory (>= 2). Default is the card's + ``default_num_timesteps``. + seed : int or None, optional + Seed for initial conditions, controls, and SSA paths. Default is + the card's ``default_seed``. + simulator : {"card", "ssa"}, optional + ``"card"`` (default) uses the card's deterministic simulator; + ``"ssa"`` replaces it with Gillespie paths at ``volume`` + (uncontrolled, non-Cantera cards only). + volume : float, optional + SSA volume (counts = ``round(c0 * volume)``). Default is + ``1000.0``. + + Returns + ------- + ReactionDataset + Concentrations, controls, and the card. + + Raises + ------ + ValueError + If arguments are invalid or ``simulator="ssa"`` does not apply. + ImportError + For the Cantera card when the ``[cantera]`` extra is missing. + """ + card = cls.load_card(dataset_id) + n_traj = _count( + "num_trajectories", + card["default_num_trajectories"] + if num_trajectories is None + else num_trajectories, + low=1, + ) + n_steps = _count( + "num_timesteps", + card["default_num_timesteps"] if num_timesteps is None else num_timesteps, + low=2, + ) + resolved_seed = _count( + "seed", card["default_seed"] if seed is None else seed, low=0 + ) + kind = card["simulator"]["kind"] + if simulator not in ("card", "ssa"): + msg = f"simulator must be 'card' or 'ssa', got {simulator!r}" + raise ValueError(msg) + if simulator == "ssa" and kind not in ("expm", "rk4"): + msg = f"simulator='ssa' is not available for {kind!r} cards" + raise ValueError(msg) + + generator = torch.Generator().manual_seed(resolved_seed) + times = float(card["time_step"]) * torch.arange(n_steps, dtype=torch.float64) + system = system_from_card(card) + controls: Tensor | None = None + used_volume: float | None = None + if kind == "cantera": + concentrations = _simulate_cantera(card, times, n_traj, generator) + else: + initial = _sample_ranges( + card["initial_ranges"], system.species, n_traj, generator + ) + if simulator == "ssa": + used_volume = float(volume) + concentrations = _simulate_ssa( + system, initial, times, used_volume, resolved_seed + ) + elif kind == "cstr": + concentrations, controls = _simulate_cstr( + card, initial, times, generator + ) + else: + method = "expm" if kind == "expm" else "rk4" + options = ( + {} + if method == "expm" + else {"substeps": int(card["simulator"]["substeps"])} + ) + states = integrate_reaction_system( + system, initial, times, method=method, **options + ) + concentrations = states.permute(1, 0, 2).contiguous() + floor = (card.get("transform") or {}).get("floor") + return ReactionDataset( + dataset_id=dataset_id, + card=card, + system=system, + times=times, + concentrations=concentrations, + controls=controls, + transform=None if floor is None else LogConcentrationTransform(floor), + seed=resolved_seed, + simulator="ssa" if simulator == "ssa" else kind, + volume=used_volume, + ) + + @staticmethod + def split( + dataset: ReactionDataset, + *, + train_ratio: float = 0.7, + val_ratio: float = 0.1, + test_ratio: float = 0.2, + ) -> ReactionDatasetSplit: + """Split trajectories deterministically by index. + + Trajectories are i.i.d. draws, so contiguous index blocks are an + unbiased split; the same seed always yields the same split. + + Parameters + ---------- + dataset : ReactionDataset + Dataset with at least three trajectories. + train_ratio, val_ratio, test_ratio : float, optional + Fractions summing to ``1``. Defaults are ``0.7 / 0.1 / 0.2``. + + Returns + ------- + ReactionDatasetSplit + Each part holds at least one trajectory. + + Raises + ------ + ValueError + If ratios are invalid or there are too few trajectories. + """ + ratios = (train_ratio, val_ratio, test_ratio) + if ( + any(r < 0 or not math.isfinite(r) for r in ratios) + or abs(sum(ratios) - 1.0) > 1e-6 + ): + msg = f"ratios must be non-negative and sum to 1, got {ratios}" + raise ValueError(msg) + n = dataset.num_trajectories + if n < 3: + msg = f"need at least 3 trajectories to split, got {n}" + raise ValueError(msg) + n_val = max(1, round(val_ratio * n)) + n_test = max(1, round(test_ratio * n)) + n_train = n - n_val - n_test + if n_train < 1: + msg = f"{n} trajectories leave no training data for ratios {ratios}" + raise ValueError(msg) + return ReactionDatasetSplit( + train=dataset.subset(list(range(n_train))), + val=dataset.subset(list(range(n_train, n_train + n_val))), + test=dataset.subset(list(range(n_train + n_val, n))), + ) + + +def system_from_card(card: dict[str, Any]) -> ReactionSystem: + """Build the mass-action system a card's generator section describes. + + Parameters + ---------- + card : dict + Parsed card. + + Returns + ------- + ReactionSystem + The system (Cantera cards need the ``[cantera]`` extra). + + Raises + ------ + ValueError + If the generator kind or catalog name is unknown. + """ + generator = card["generator"] + kind = generator["kind"] + params = dict(generator.get("params", {})) + if kind == "catalog": + name = generator["name"] + if name not in _CATALOG_FACTORIES: + msg = f"unknown catalog system {name!r}" + raise ValueError(msg) + return _CATALOG_FACTORIES[name](**params) + if kind == "random_network": + return sample_random_network(**params).system + if kind == "cstr": + return van_de_vusse_cstr(**params).system + if kind == "cantera": + from koopman_graph.datasets.reaction.cantera_io import load_cantera_mechanism + + cantera = card["cantera"] + return load_cantera_mechanism( + cantera["mechanism"], + temperature_k=cantera["temperature_k"], + pressure_pa=cantera["pressure_pa"], + mole_fractions=cantera["reference_mole_fractions"], + state_dependent=cantera["state_dependent"], + ).system + msg = f"unknown generator kind {kind!r}" + raise ValueError(msg) + + +def compute_card_oracle(card: dict[str, Any]) -> dict[str, Any]: + """Compute the oracle section for a card from its generator spec. + + Parameters + ---------- + card : dict + Card with ``generator``, ``initial_ranges`` (or ``cantera``) and, + for the CSTR, ``controls``. + + Returns + ------- + dict + ``state`` (how the reference point was found), ``reference_state``, + ``num_conservation_laws``, ``conservation_laws``, + ``dynamic_eigenvalues_real`` / ``_imag`` (slowest first), and + ``timescales`` (``None`` for infinite). For the CSTR the spectrum + is that of the controlled field ``J - D I`` at the nominal-input + steady state. + """ + kind = card["generator"]["kind"] + system = system_from_card(card) + if kind == "cstr": + controls = card["controls"]["ranges"] + temperature = sum(controls["temperature_k"]) / 2.0 + dilution = sum(controls["dilution_rate"]) / 2.0 + reactor = van_de_vusse_cstr(**card["generator"].get("params", {})) + state = reactor.steady_state(temperature_k=temperature, dilution_rate=dilution) + jac = mass_action_jacobian(reactor.system.at_temperature(temperature), state) + jac = jac - dilution * torch.eye(system.num_species, dtype=torch.float64) + eigenvalues = torch.linalg.eigvals(jac) + eigenvalues = eigenvalues[torch.argsort(eigenvalues.real, descending=True)] + return _oracle_dict( + "cstr_steady_state_at_control_midpoint", + state, + eigenvalues, + torch.zeros(0, system.num_species, dtype=torch.float64), + ) + if kind == "cantera": + from koopman_graph.datasets.reaction.cantera_io import require_cantera + + ct = require_cantera() + cantera = card["cantera"] + gas = ct.Solution(cantera["mechanism"]) + gas.TPX = ( + cantera["temperature_k"], + cantera["pressure_pa"], + cantera["reference_mole_fractions"], + ) + gas.equilibrate("TV") + state = torch.tensor(gas.concentrations, dtype=torch.float64) + label = "cantera_thermochemical_equilibrium_TV" + else: + reference = _range_midpoints(card["initial_ranges"], system.species) + relax = card.get("oracle_relaxation_time") + guess = None + label = "equilibrium_from_range_midpoints" + if relax is not None: + # Warm-start Newton from an RK4 relaxation (networks that drain + # species onto a boundary steady state converge slowly otherwise). + steps_per_unit = int(card["simulator"]["substeps"]) / float( + card["time_step"] + ) + horizon = torch.tensor([0.0, float(relax)], dtype=torch.float64) + guess = integrate_reaction_system( + system, + reference, + horizon, + method="rk4", + substeps=max(1, math.ceil(steps_per_unit * float(relax))), + )[-1] + label = "equilibrium_from_range_midpoints_after_relaxation" + state = equilibrium(system, reference, guess=guess) + spectrum = oracle_spectrum(system, state) + return _oracle_dict( + label, state, spectrum.dynamic_eigenvalues, conservation_laws(system) + ) + + +def _oracle_dict( + label: str, state: Tensor, eigenvalues: Tensor, laws: Tensor +) -> dict[str, Any]: + """Pack oracle quantities into JSON-compatible primitives. + + Parameters + ---------- + label + See signature. + state + See signature. + eigenvalues + See signature. + laws + See signature. + + Returns + ------- + See signature.""" + real = eigenvalues.real + timescales = [ + None if math.isinf(t) else t + for t in eigenvalue_timescales(eigenvalues).tolist() + ] + return { + "state": label, + "reference_state": state.tolist(), + "num_conservation_laws": int(laws.shape[0]), + "conservation_laws": laws.tolist(), + "dynamic_eigenvalues_real": real.tolist(), + "dynamic_eigenvalues_imag": eigenvalues.imag.tolist(), + "timescales": timescales, + } + + +def _range_midpoints( + ranges: dict[str, list[float]], species: tuple[str, ...] +) -> Tensor: + """Return per-species range midpoints in species order. + + Parameters + ---------- + ranges + See signature. + species + See signature. + + Returns + ------- + See signature.""" + _check_ranges(ranges, species) + return torch.tensor( + [(ranges[s][0] + ranges[s][1]) / 2.0 for s in species], dtype=torch.float64 + ) + + +def _sample_ranges( + ranges: dict[str, list[float]], + species: tuple[str, ...], + count: int, + generator: torch.Generator, +) -> Tensor: + """Draw ``count`` uniform samples inside per-species ranges. + + Parameters + ---------- + ranges + See signature. + species + See signature. + count + See signature. + generator + See signature. + + Returns + ------- + See signature.""" + _check_ranges(ranges, species) + low = torch.tensor([ranges[s][0] for s in species], dtype=torch.float64) + high = torch.tensor([ranges[s][1] for s in species], dtype=torch.float64) + u = torch.rand(count, len(species), generator=generator, dtype=torch.float64) + return low + (high - low) * u + + +def _check_ranges(ranges: dict[str, list[float]], species: tuple[str, ...]) -> None: + """Validate that ``ranges`` covers exactly ``species`` with ``0 <= lo <= hi``. + + Parameters + ---------- + ranges + See signature. + species + See signature.""" + if set(ranges) != set(species): + msg = ( + f"initial_ranges keys {sorted(ranges)} do not match species {list(species)}" + ) + raise ValueError(msg) + for name, (lo, hi) in ranges.items(): + if not (0.0 <= lo <= hi) or not math.isfinite(hi): + msg = ( + f"initial range for {name!r} must satisfy 0 <= low <= high, " + f"got {[lo, hi]}" + ) + raise ValueError(msg) + + +def _simulate_ssa( + system: ReactionSystem, + initial: Tensor, + times: Tensor, + volume: float, + seed: int, +) -> Tensor: + """Run one Gillespie path per initial state and sample it on ``times``. + + Parameters + ---------- + system + See signature. + initial + See signature. + times + See signature. + volume + See signature. + seed + See signature. + + Returns + ------- + See signature.""" + if not math.isfinite(volume) or volume <= 0.0: + msg = f"volume must be finite and > 0, got {volume!r}" + raise ValueError(msg) + t_end = float(times[-1]) + paths = [] + for i, c0 in enumerate(initial): + counts = torch.round(c0 * volume).to(torch.long) + path = gillespie_ssa( + system, counts, t_end, volume=volume, seed=seed * 100_003 + i + ) + paths.append(path.concentrations_on_grid(times)) + return torch.stack(paths) + + +def _simulate_cstr( + card: dict[str, Any], + initial: Tensor, + times: Tensor, + generator: torch.Generator, +) -> tuple[Tensor, Tensor]: + """Simulate the CSTR under seeded piecewise-constant control schedules. + + Parameters + ---------- + card + See signature. + initial + See signature. + times + See signature. + generator + See signature. + + Returns + ------- + See signature.""" + reactor = van_de_vusse_cstr(**card["generator"].get("params", {})) + spec = card["controls"] + names = spec["names"] + hold = int(spec["hold_steps"]) + substeps = int(card["simulator"]["substeps"]) + num_intervals = int(times.shape[0]) - 1 + num_blocks = math.ceil(num_intervals / hold) + states, schedules = [], [] + for c0 in initial: + columns = [] + for name in names: + lo, hi = spec["ranges"][name] + block = lo + (hi - lo) * torch.rand( + num_blocks, generator=generator, dtype=torch.float64 + ) + columns.append(block.repeat_interleave(hold)[:num_intervals]) + schedule = torch.stack(columns, dim=-1) # (T - 1, C) + states.append( + reactor.simulate( + c0, + times, + temperature_k=schedule[:, names.index("temperature_k")], + dilution_rate=schedule[:, names.index("dilution_rate")], + substeps=substeps, + ) + ) + schedules.append(torch.cat([schedule, schedule[-1:]], dim=0)) + return torch.stack(states), torch.stack(schedules) + + +def _simulate_cantera( + card: dict[str, Any], + times: Tensor, + count: int, + generator: torch.Generator, +) -> Tensor: + """Run seeded Cantera reactors from mole-fraction ranges. + + Parameters + ---------- + card + See signature. + times + See signature. + count + See signature. + generator + See signature. + + Returns + ------- + See signature.""" + from koopman_graph.datasets.reaction.cantera_io import simulate_cantera_reactor + + cantera = card["cantera"] + ranges = cantera["initial_mole_fraction_ranges"] + names = sorted(ranges) + runs = [] + for _ in range(count): + u = torch.rand(len(names), generator=generator, dtype=torch.float64) + mix = { + name: ranges[name][0] + (ranges[name][1] - ranges[name][0]) * float(u[k]) + for k, name in enumerate(names) + } + trajectory = simulate_cantera_reactor( + cantera["mechanism"], + times=times, + temperature_k=cantera["temperature_k"], + pressure_pa=cantera["pressure_pa"], + mole_fractions=mix, + reactor=cantera["reactor"], + isothermal=cantera["isothermal"], + ) + runs.append(trajectory.concentrations) + return torch.stack(runs) + + +def _count(name: str, value: Any, *, low: int) -> int: + """Validate an integer argument ``>= low``. + + Parameters + ---------- + name + See signature. + value + See signature. + low + See signature. + + Returns + ------- + See signature.""" + if isinstance(value, bool) or not isinstance(value, int) or value < low: + msg = f"{name} must be an int >= {low}, got {value!r}" + raise ValueError(msg) + return value diff --git a/src/koopman_graph/datasets/reaction/data/.gitkeep b/src/koopman_graph/datasets/reaction/data/.gitkeep deleted file mode 100644 index e69de29..0000000 diff --git a/src/koopman_graph/datasets/reaction/data/reaction_brusselator_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_brusselator_v1.json new file mode 100644 index 0000000..f744660 --- /dev/null +++ b/src/koopman_graph/datasets/reaction/data/reaction_brusselator_v1.json @@ -0,0 +1,67 @@ +{ + "dataset_id": "reaction_brusselator_v1", + "version": "v1", + "scope": "Open Brusselator (a = 1, b = 3) with an unstable focus and a limit cycle; tests oscillatory, non-decaying dynamics.", + "generator": { + "kind": "catalog", + "name": "brusselator", + "params": {} + }, + "simulator": { + "kind": "rk4", + "substeps": 10 + }, + "time_step": 0.1, + "default_num_timesteps": 128, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "X": [ + 0.5, + 1.5 + ], + "Y": [ + 2.0, + 4.0 + ] + }, + "transform": null, + "graph_encoding": "species_influence", + "limitations": "Open system without conservation laws; the oracle fixed point is unstable, so its timescales are negative (growth).", + "units": { + "time": "s", + "concentration": "mol/L" + }, + "controls": null, + "cantera": null, + "format": "Per-trajectory GraphSnapshotSequence (one node per species, float32 features of shape (num_species, 1), timestamps, optional global control_inputs); raw float64 concentrations in ReactionDataset.", + "provenance": "In-repo mass-action generator (koopman_graph.datasets.reaction).", + "license": "Apache-2.0", + "species": [ + "X", + "Y" + ], + "elements": null, + "num_reactions": 4, + "oracle": { + "state": "equilibrium_from_range_midpoints", + "reference_state": [ + 1.0, + 3.0 + ], + "num_conservation_laws": 0, + "conservation_laws": [], + "dynamic_eigenvalues_real": [ + 0.5, + 0.5 + ], + "dynamic_eigenvalues_imag": [ + 0.8660254037844388, + -0.8660254037844388 + ], + "timescales": [ + -2.0, + -2.0 + ] + } +} diff --git a/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json new file mode 100644 index 0000000..ce0fe23 --- /dev/null +++ b/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json @@ -0,0 +1,167 @@ +{ + "dataset_id": "reaction_cantera_h2o2_v1", + "version": "v1", + "scope": "Isothermal constant-volume H2/O2/Ar chemistry at 1200 K from Cantera's bundled h2o2.yaml (10 species, 29 reactions); trajectories come from Cantera's ReactorNet.", + "generator": { + "kind": "cantera", + "params": {} + }, + "simulator": { + "kind": "cantera" + }, + "time_step": 2e-05, + "default_num_timesteps": 50, + "default_num_trajectories": 16, + "default_seed": 0, + "cantera": { + "mechanism": "h2o2.yaml", + "sha256": "0efc6c52862741a29e0c29b65d979c7d8cb409db5282bca83b9c5437b3d8c8d4", + "cantera_version": "3.2.0", + "temperature_k": 1200.0, + "pressure_pa": 101325.0, + "reactor": "const_volume", + "isothermal": true, + "state_dependent": "freeze", + "reference_mole_fractions": "H2:2, O2:1, AR:4", + "initial_mole_fraction_ranges": { + "H2": [ + 1.0, + 3.0 + ], + "O2": [ + 0.5, + 1.5 + ], + "AR": [ + 3.0, + 5.0 + ] + }, + "note": "Cantera's bundled input files are provided for illustration only; they are not validated kinetic models for quantitative use.", + "citation_doi": "10.5281/zenodo.17620923", + "license": "Cantera is BSD-3-Clause" + }, + "transform": { + "kind": "log_concentration", + "floor": 1e-12 + }, + "graph_encoding": "species_influence", + "limitations": "Illustration-only mechanism; the graph and oracle use the mass-action conversion with the falloff reaction frozen at the reference state, while trajectories are exact Cantera integrations.", + "units": { + "time": "s", + "concentration": "mol/L" + }, + "controls": null, + "format": "Per-trajectory GraphSnapshotSequence (one node per species, float32 features of shape (num_species, 1), timestamps, optional global control_inputs); raw float64 concentrations in ReactionDataset.", + "provenance": "Cantera ReactorNet integration of the bundled mechanism; see the cantera section for version and SHA-256.", + "license": "Apache-2.0", + "species": [ + "H2", + "H", + "O", + "O2", + "OH", + "H2O", + "HO2", + "H2O2", + "AR", + "N2" + ], + "elements": [ + "O", + "H", + "Ar", + "N" + ], + "num_reactions": 136, + "oracle": { + "state": "cantera_thermochemical_equilibrium_TV", + "reference_state": [ + 3.0669032201594336e-08, + 3.482844896613705e-12, + 3.059422020259053e-13, + 1.4677299229464722e-08, + 2.6314766744630864e-09, + 0.0029015385787945157, + 3.1032512185622856e-14, + 8.412261154894816e-14, + 0.00580314113081214, + 0.0 + ], + "num_conservation_laws": 4, + "conservation_laws": [ + [ + 0.28064053942470113, + 0.1403202697123506, + -0.3167114270242802, + -0.6334228540485629, + -0.17639115731193028, + -0.036070887599579333, + -0.493102584336212, + -0.3527823146238605, + 0.0, + 0.0 + ], + [ + 0.5812981629921767, + 0.2906490814960885, + -0.062133233127772275, + -0.12426646625554473, + 0.2285158483683161, + 0.5191649298644041, + 0.16638261524054332, + 0.4570316967366318, + 0.0, + 0.0 + ], + [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0 + ], + [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 0.0 + ] + ], + "dynamic_eigenvalues_real": [ + -0.000469766753797218, + -0.023066363423131318, + -390.50523712580406, + -4998.956241475473, + -48124.93946401666, + -53087.23619072366 + ], + "dynamic_eigenvalues_imag": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "timescales": [ + 2128.7159891942147, + 43.35317109402621, + 0.0025607851186841904, + 0.0002000417590582557, + 2.0779246917239383e-05, + 1.8836919601678898e-05 + ] + } +} diff --git a/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json new file mode 100644 index 0000000..b4826d1 --- /dev/null +++ b/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json @@ -0,0 +1,95 @@ +{ + "dataset_id": "reaction_chain_v1", + "version": "v1", + "scope": "Linear first-order chain S0 -> S1 -> S2 -> S3 (rates 1, 0.5, 0.25); exact matrix-exponential trajectories with closed-form dynamic spectrum {-0.25, -0.5, -1} plus one conservation zero mode.", + "generator": { + "kind": "catalog", + "name": "first_order_chain", + "params": { + "num_species": 4 + } + }, + "simulator": { + "kind": "expm" + }, + "time_step": 0.25, + "default_num_timesteps": 64, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "S0": [ + 0.5, + 1.5 + ], + "S1": [ + 0.0, + 0.2 + ], + "S2": [ + 0.0, + 0.2 + ], + "S3": [ + 0.0, + 0.2 + ] + }, + "transform": { + "kind": "log_concentration", + "floor": 1e-08 + }, + "graph_encoding": "species_influence", + "limitations": "Linear and closed; every mode is a pure decay. Use as the exact-spectrum sanity check, not as a nonlinear benchmark.", + "units": { + "time": "s", + "concentration": "mol/L" + }, + "controls": null, + "cantera": null, + "format": "Per-trajectory GraphSnapshotSequence (one node per species, float32 features of shape (num_species, 1), timestamps, optional global control_inputs); raw float64 concentrations in ReactionDataset.", + "provenance": "In-repo mass-action generator (koopman_graph.datasets.reaction).", + "license": "Apache-2.0", + "species": [ + "S0", + "S1", + "S2", + "S3" + ], + "elements": [ + "M" + ], + "num_reactions": 3, + "oracle": { + "state": "equilibrium_from_range_midpoints", + "reference_state": [ + -3.514983438626458e-17, + -8.784035482230695e-17, + -1.378654434722274e-16, + 1.300000000000001 + ], + "num_conservation_laws": 1, + "conservation_laws": [ + [ + 0.5000000000000001, + 0.5000000000000001, + 0.5, + 0.49999999999999983 + ] + ], + "dynamic_eigenvalues_real": [ + -0.25, + -0.5, + -1.0 + ], + "dynamic_eigenvalues_imag": [ + 0.0, + 0.0, + 0.0 + ], + "timescales": [ + 4.0, + 2.0, + 1.0 + ] + } +} diff --git a/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json new file mode 100644 index 0000000..e728cbf --- /dev/null +++ b/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json @@ -0,0 +1,98 @@ +{ + "dataset_id": "reaction_michaelis_menten_v1", + "version": "v1", + "scope": "Mass-action Michaelis-Menten E + S <-> ES -> E + P with two conservation laws (enzyme, substrate backbone) and a boundary equilibrium.", + "generator": { + "kind": "catalog", + "name": "michaelis_menten", + "params": {} + }, + "simulator": { + "kind": "rk4", + "substeps": 50 + }, + "time_step": 0.5, + "default_num_timesteps": 64, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "E": [ + 0.05, + 0.15 + ], + "S": [ + 0.5, + 1.5 + ], + "ES": [ + 0.0, + 0.05 + ], + "P": [ + 0.0, + 0.1 + ] + }, + "transform": { + "kind": "log_concentration", + "floor": 1e-08 + }, + "graph_encoding": "reaction_hyperedge", + "limitations": "Irreversible catalysis: trajectories end on the boundary S = ES = 0, where the log transform saturates at the floor.", + "units": { + "time": "s", + "concentration": "mol/L" + }, + "controls": null, + "cantera": null, + "format": "Per-trajectory GraphSnapshotSequence (one node per species, float32 features of shape (num_species, 1), timestamps, optional global control_inputs); raw float64 concentrations in ReactionDataset.", + "provenance": "In-repo mass-action generator (koopman_graph.datasets.reaction).", + "license": "Apache-2.0", + "species": [ + "E", + "S", + "ES", + "P" + ], + "elements": [ + "Enz", + "Sub" + ], + "num_reactions": 3, + "oracle": { + "state": "equilibrium_from_range_midpoints", + "reference_state": [ + 0.12499999999999997, + 0.0, + 0.0, + 1.0750000000000002 + ], + "num_conservation_laws": 2, + "conservation_laws": [ + [ + -0.7123282872413622, + 0.010643277802640164, + -0.701685009438722, + 0.01064327780264014 + ], + [ + -0.3042834389114655, + 0.6323659704930491, + 0.32808253158158357, + 0.632365970493049 + ] + ], + "dynamic_eigenvalues_real": [ + -0.25, + -2.5 + ], + "dynamic_eigenvalues_imag": [ + 0.0, + 0.0 + ], + "timescales": [ + 4.0, + 0.4 + ] + } +} diff --git a/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json new file mode 100644 index 0000000..422a223 --- /dev/null +++ b/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json @@ -0,0 +1,131 @@ +{ + "dataset_id": "reaction_random_network_v1", + "version": "v1", + "scope": "Seeded random element-balanced mass-action network (6 species, 2 elements, 8 reactions, order <= 2).", + "generator": { + "kind": "random_network", + "params": { + "num_species": 6, + "num_elements": 2, + "num_reactions": 8, + "max_atoms": 4, + "max_order": 2, + "rate_log10_range": [ + -1.0, + 1.0 + ], + "seed": 0 + } + }, + "simulator": { + "kind": "rk4", + "substeps": 20 + }, + "time_step": 0.1, + "default_num_timesteps": 64, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "X0": [ + 0.1, + 1.0 + ], + "X1": [ + 0.1, + 1.0 + ], + "X2": [ + 0.1, + 1.0 + ], + "X3": [ + 0.1, + 1.0 + ], + "X4": [ + 0.1, + 1.0 + ], + "X5": [ + 0.1, + 1.0 + ] + }, + "oracle_relaxation_time": 200.0, + "transform": { + "kind": "log_concentration", + "floor": 1e-08 + }, + "graph_encoding": "species_influence", + "limitations": "Synthetic chemistry with no physical species; rates are log-uniform draws, not measured kinetics.", + "units": { + "time": "s", + "concentration": "mol/L" + }, + "controls": null, + "cantera": null, + "format": "Per-trajectory GraphSnapshotSequence (one node per species, float32 features of shape (num_species, 1), timestamps, optional global control_inputs); raw float64 concentrations in ReactionDataset.", + "provenance": "In-repo mass-action generator (koopman_graph.datasets.reaction).", + "license": "Apache-2.0", + "species": [ + "X0", + "X1", + "X2", + "X3", + "X4", + "X5" + ], + "elements": [ + "e0", + "e1" + ], + "num_reactions": 8, + "oracle": { + "state": "equilibrium_from_range_midpoints_after_relaxation", + "reference_state": [ + 0.9371973830410871, + 2.668598691520542, + 7.203240409913488e-32, + 1.181401308479456, + 2.8981443113177526e-30, + 1.6307572951339986e-30 + ], + "num_conservation_laws": 2, + "conservation_laws": [ + [ + 0.2889695340165011, + -0.4048248518165172, + -0.11585531780001589, + 0.1731142162164851, + -0.23171063560003166, + -0.8096497036330341 + ], + [ + 0.31116631768677233, + 0.00715999869017006, + 0.3183263163769422, + 0.6294926340637144, + 0.6366526327538841, + 0.014319997380339834 + ] + ], + "dynamic_eigenvalues_real": [ + -1.1843321716463685e-17, + -0.8165944002656043, + -0.8249568943108102, + -30.345958028796638 + ], + "dynamic_eigenvalues_imag": [ + 0.0, + 0.0, + 0.0, + 0.0 + ], + "timescales": [ + null, + 1.2245981599613487, + 1.2121845479398354, + 0.032953317837290066 + ] + } +} diff --git a/src/koopman_graph/datasets/reaction/data/reaction_van_de_vusse_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_van_de_vusse_v1.json new file mode 100644 index 0000000..2d80565 --- /dev/null +++ b/src/koopman_graph/datasets/reaction/data/reaction_van_de_vusse_v1.json @@ -0,0 +1,108 @@ +{ + "dataset_id": "reaction_van_de_vusse_v1", + "version": "v1", + "scope": "Van de Vusse CSTR (A -> B -> C, 2 A -> D) with seeded piecewise-constant reactor temperature and dilution-rate inputs stored as global control_inputs.", + "generator": { + "kind": "cstr", + "params": {} + }, + "simulator": { + "kind": "cstr", + "substeps": 20 + }, + "time_step": 0.01, + "default_num_timesteps": 100, + "default_num_trajectories": 32, + "default_seed": 0, + "initial_ranges": { + "A": [ + 0.5, + 3.0 + ], + "B": [ + 0.0, + 1.5 + ], + "C": [ + 0.0, + 1.5 + ], + "D": [ + 0.0, + 0.5 + ] + }, + "controls": { + "names": [ + "temperature_k", + "dilution_rate" + ], + "ranges": { + "temperature_k": [ + 380.0, + 400.0 + ], + "dilution_rate": [ + 5.0, + 30.0 + ] + }, + "hold_steps": 10, + "units": { + "temperature_k": "K", + "dilution_rate": "1/h" + } + }, + "transform": null, + "graph_encoding": "species_influence", + "units": { + "time": "h", + "concentration": "mol/L" + }, + "references": [ + "K.-U. Klatt and S. Engell, Comput. Chem. Eng. 22(4-5), 491-502, 1998 (kinetics)", + "H. Chen, A. Kremling and F. Allgower, Proc. 3rd ECC, 1995 (operating point)" + ], + "limitations": "Reactor temperature is a direct input (no energy balance); the graph uses the reaction network only, not the feed.", + "cantera": null, + "format": "Per-trajectory GraphSnapshotSequence (one node per species, float32 features of shape (num_species, 1), timestamps, optional global control_inputs); raw float64 concentrations in ReactionDataset.", + "provenance": "In-repo mass-action generator (koopman_graph.datasets.reaction).", + "license": "Apache-2.0", + "species": [ + "A", + "B", + "C", + "D" + ], + "elements": null, + "num_reactions": 3, + "oracle": { + "state": "cstr_steady_state_at_control_midpoint", + "reference_state": [ + 2.187229135690542, + 1.0935337017908153, + 1.0933719935983206, + 0.3629325844601609 + ], + "num_conservation_laws": 0, + "conservation_laws": [], + "dynamic_eigenvalues_real": [ + -17.5, + -17.5, + -34.99741215715252, + -46.6126931067593 + ], + "dynamic_eigenvalues_imag": [ + 0.0, + 0.0, + 0.0, + 0.0 + ], + "timescales": [ + 0.05714285714285714, + 0.05714285714285714, + 0.028573541252410208, + 0.021453383903600073 + ] + } +} diff --git a/src/koopman_graph/datasets/reaction/oracles.py b/src/koopman_graph/datasets/reaction/oracles.py index b6ab905..00bcda5 100644 --- a/src/koopman_graph/datasets/reaction/oracles.py +++ b/src/koopman_graph/datasets/reaction/oracles.py @@ -67,15 +67,12 @@ def timescales(self) -> Tensor: ------- Tensor ``float64`` tensor aligned with :attr:`dynamic_eigenvalues`: - positive for decaying modes, ``inf`` when ``Re(λ) = 0`` (e.g. a - neutral centre), and negative for growing modes. + positive for decaying modes, ``inf`` when ``Re(λ)`` is zero to + within ``1e-12 · max|λ|`` (a neutral centre, or a species + drained onto a zero-flux boundary), and negative for growing + modes. """ - real = self.dynamic_eigenvalues.real - return torch.where( - real == 0, - torch.full_like(real, math.inf), - -1.0 / real, - ) + return eigenvalue_timescales(self.dynamic_eigenvalues) @property def is_stable(self) -> bool: @@ -115,6 +112,30 @@ def discrete_eigenvalues(self, time_step: float) -> Tensor: return torch.exp(self.dynamic_eigenvalues * dt) +_ZERO_EIGENVALUE_RTOL = 1e-12 + + +def eigenvalue_timescales(eigenvalues: Tensor) -> Tensor: + """Return relaxation timescales ``-1 / Re(λ)`` with numerical zeros as ``inf``. + + Parameters + ---------- + eigenvalues : Tensor + Complex (or real) eigenvalues ``(k,)``. + + Returns + ------- + Tensor + ``float64`` timescales ``(k,)``. ``|Re(λ)| <= 1e-12 · max|λ|`` maps + to ``inf``. + """ + real = eigenvalues.real.to(torch.float64) + scale = float(eigenvalues.abs().max()) if eigenvalues.numel() else 0.0 + zero = real.abs() <= _ZERO_EIGENVALUE_RTOL * scale + safe = torch.where(zero, torch.ones_like(real), real) + return torch.where(zero, torch.full_like(real, math.inf), -1.0 / safe) + + def stoichiometric_rank(system: ReactionSystem) -> int: """Return the rank of the stoichiometric matrix ``N``. @@ -232,6 +253,8 @@ def residual(c: Tensor) -> Tensor: start = c0 if guess is None else _state(system, guess, "guess") c = start.clamp(min=0) + # Turnover at the initial state sets the absolute stationarity scale. + flux_scale = float(_gross_turnover(system, c0.clamp(min=0)).max()) f = residual(c) for _ in range(max_iter): jac = torch.cat([q.T @ mass_action_jacobian(system, c), w]) @@ -247,7 +270,7 @@ def residual(c: Tensor) -> Tensor: moved = float((candidate - c).abs().max()) c, f = candidate, f_candidate if moved <= tol * (1.0 + float(c.abs().max())): - _require_stationary(system, c) + _require_stationary(system, c, flux_scale=flux_scale) return c msg = ( f"Newton did not converge for {system.name!r} in {max_iter} iterations " @@ -316,21 +339,50 @@ def oracle_timescales(system: ReactionSystem, state: Tensor) -> Tensor: _STATIONARY_RTOL = 1e-6 -def _require_stationary(system: ReactionSystem, c: Tensor) -> None: - """Raise unless net rates are negligible against gross species turnover. +_STATIONARY_ATOL = 1e-12 + + +def _gross_turnover(system: ReactionSystem, c: Tensor) -> Tensor: + """Return per-species gross turnover ``Σ_j r_j (α_js + β_js)``. Parameters ---------- system See signature. c + See signature. + + Returns + ------- See signature.""" rates = system.reaction_rates(c) + return rates @ (system.reactant_stoichiometry + system.product_stoichiometry) + + +def _require_stationary( + system: ReactionSystem, c: Tensor, *, flux_scale: float = 0.0 +) -> None: + """Raise unless net rates are negligible against species turnover. + + A species passes when ``|dc_s/dt| <= rtol * gross_s + atol * flux_scale``: + the relative term rejects one-way flux at a clamped boundary point, and + the absolute term (scaled by the peak gross turnover at ``initial``) + accepts vanishing leftovers on the way to a boundary steady state. + + Parameters + ---------- + system + See signature. + c + See signature. + flux_scale + See signature.""" net = system.time_derivative(c).abs() - gross = rates @ (system.reactant_stoichiometry + system.product_stoichiometry) - tiny = torch.finfo(torch.float64).tiny - ratio = float((net / gross.clamp(min=tiny)).max()) - if ratio > _STATIONARY_RTOL: + gross = _gross_turnover(system, c) + allowed = _STATIONARY_RTOL * gross + _STATIONARY_ATOL * flux_scale + if bool((net > allowed).any()): + tiny = torch.finfo(torch.float64).tiny + ratio = float((net / gross.clamp(min=tiny)).max()) msg = ( f"Newton stalled at a non-stationary point for {system.name!r} " f"(max net/gross species turnover {ratio:.3e}); pass a guess " diff --git a/tests/datasets/test_reaction_benchmark.py b/tests/datasets/test_reaction_benchmark.py new file mode 100644 index 0000000..52120a5 --- /dev/null +++ b/tests/datasets/test_reaction_benchmark.py @@ -0,0 +1,426 @@ +"""Tests for :mod:`koopman_graph.datasets.reaction.benchmark`.""" + +from __future__ import annotations + +import importlib.util +import math +import subprocess +import sys +from pathlib import Path + +import pytest +import torch +from tests.helpers import REPO_ROOT + +from koopman_graph.data import GraphSnapshotSequence +from koopman_graph.datasets.reaction import ( + PAYLOAD_FORMAT_VERSION, + ReactionDataset, + ReactionNetworkBenchmark, + check_element_balance, + compute_card_oracle, + integrate_reaction_system, + oracle_spectrum, +) + +_NON_CANTERA = ( + "reaction_brusselator_v1", + "reaction_chain_v1", + "reaction_michaelis_menten_v1", + "reaction_random_network_v1", + "reaction_van_de_vusse_v1", +) +_CLOSED = ( + "reaction_chain_v1", + "reaction_michaelis_menten_v1", + "reaction_random_network_v1", +) +_REQUIRED_CARD_KEYS = { + "dataset_id", + "version", + "scope", + "generator", + "simulator", + "time_step", + "default_num_timesteps", + "default_num_trajectories", + "default_seed", + "transform", + "graph_encoding", + "species", + "num_reactions", + "oracle", + "units", + "format", + "provenance", + "license", + "limitations", +} + + +def _close(a: object, b: object, *, rtol: float = 1e-8, atol: float = 1e-12) -> bool: + """Compare nested JSON values with float tolerance. + + Parameters + ---------- + a + See signature. + b + See signature. + rtol + See signature. + atol + See signature. + + Returns + ------- + See signature.""" + if isinstance(a, list) and isinstance(b, list): + return len(a) == len(b) and all( + _close(x, y, rtol=rtol, atol=atol) for x, y in zip(a, b, strict=True) + ) + if isinstance(a, float) and isinstance(b, float): + return math.isclose(a, b, rel_tol=rtol, abs_tol=atol) + return a == b + + +def test_packaged_cards_are_listed_and_complete() -> None: + ids = ReactionNetworkBenchmark.dataset_ids() + assert set(ids) == {*_NON_CANTERA, "reaction_cantera_h2o2_v1"} + for dataset_id in ids: + card = ReactionNetworkBenchmark.load_card(dataset_id) + assert set(card) >= _REQUIRED_CARD_KEYS, dataset_id + assert card["dataset_id"] == dataset_id + assert card["version"] == "v1" + assert card["license"] == "Apache-2.0" + assert card["graph_encoding"] in {"species_influence", "reaction_hyperedge"} + with pytest.raises(ValueError, match="unknown reaction dataset"): + ReactionNetworkBenchmark.load_card("reaction_nope_v1") + + +@pytest.mark.parametrize("dataset_id", _NON_CANTERA) +def test_card_oracle_matches_recomputation(dataset_id: str) -> None: + card = ReactionNetworkBenchmark.load_card(dataset_id) + fresh = compute_card_oracle(card) + for key, value in card["oracle"].items(): + assert _close(fresh[key], value), f"{dataset_id}: oracle[{key!r}] drifted" + + +@pytest.mark.cantera +def test_cantera_card_oracle_and_provenance() -> None: + from koopman_graph.datasets.reaction import ( + load_cantera_mechanism, + mechanism_sha256, + resolve_mechanism_path, + ) + + card = ReactionNetworkBenchmark.load_card("reaction_cantera_h2o2_v1") + cantera = card["cantera"] + assert cantera["sha256"] == mechanism_sha256(resolve_mechanism_path("h2o2.yaml")) + assert "illustration only" in cantera["note"] + fresh = compute_card_oracle(card) + for key, value in card["oracle"].items(): + assert _close(fresh[key], value, rtol=1e-6, atol=1e-18), key + mechanism = load_cantera_mechanism( + cantera["mechanism"], + temperature_k=cantera["temperature_k"], + mole_fractions=cantera["reference_mole_fractions"], + ) + assert mechanism.system.species == tuple(card["species"]) + dataset = ReactionNetworkBenchmark.generate( + "reaction_cantera_h2o2_v1", num_trajectories=2, num_timesteps=6 + ) + assert dataset.concentrations.shape == (2, 6, 10) + assert dataset.simulator == "cantera" + assert dataset.transform is not None + assert dataset.transform.floor == 1e-12 + + +def test_chain_oracle_is_closed_form() -> None: + oracle = ReactionNetworkBenchmark.load_card("reaction_chain_v1")["oracle"] + assert oracle["num_conservation_laws"] == 1 + assert _close(oracle["dynamic_eigenvalues_real"], [-0.25, -0.5, -1.0]) + assert _close(oracle["timescales"], [4.0, 2.0, 1.0]) + + +def test_random_network_card_marks_drained_mode_infinite() -> None: + oracle = ReactionNetworkBenchmark.load_card("reaction_random_network_v1")["oracle"] + assert oracle["state"].endswith("after_relaxation") + assert oracle["timescales"][0] is None # numerically zero drained mode + assert all(t is None or t > 0 for t in oracle["timescales"]) + + +@pytest.mark.parametrize("dataset_id", _NON_CANTERA) +def test_generation_is_deterministic_by_seed(dataset_id: str) -> None: + a = ReactionNetworkBenchmark.generate( + dataset_id, num_trajectories=3, num_timesteps=8 + ) + b = ReactionNetworkBenchmark.generate( + dataset_id, num_trajectories=3, num_timesteps=8 + ) + c = ReactionNetworkBenchmark.generate( + dataset_id, num_trajectories=3, num_timesteps=8, seed=1 + ) + assert torch.equal(a.concentrations, b.concentrations) + assert not torch.equal(a.concentrations, c.concentrations) + if a.controls is not None: + assert b.controls is not None + assert torch.equal(a.controls, b.controls) + assert a.concentrations.dtype == torch.float64 + assert bool(torch.isfinite(a.concentrations).all()) + assert a.times.tolist() == pytest.approx( + [a.card["time_step"] * i for i in range(8)] + ) + + +def test_defaults_come_from_card() -> None: + dataset = ReactionNetworkBenchmark.generate("reaction_chain_v1") + card = dataset.card + assert dataset.concentrations.shape == ( + card["default_num_trajectories"], + card["default_num_timesteps"], + 4, + ) + assert dataset.seed == card["default_seed"] + assert dataset.simulator == "expm" + + +@pytest.mark.parametrize("dataset_id", _CLOSED) +def test_closed_datasets_conserve_elements(dataset_id: str) -> None: + dataset = ReactionNetworkBenchmark.generate(dataset_id, num_trajectories=4) + system = dataset.system + assert system.element_matrix is not None + assert check_element_balance(system).balanced + totals = dataset.concentrations @ system.element_matrix.T + drift = (totals - totals[:, :1]).abs().max() + assert float(drift) < 1e-9 + + +def test_initial_conditions_respect_card_ranges() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_michaelis_menten_v1", num_trajectories=50 + ) + ranges = dataset.card["initial_ranges"] + initial = dataset.concentrations[:, 0] + for s, name in enumerate(dataset.system.species): + lo, hi = ranges[name] + assert float(initial[:, s].min()) >= lo + assert float(initial[:, s].max()) <= hi + + +def test_trajectories_follow_the_card_simulator() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_brusselator_v1", num_trajectories=2, num_timesteps=12 + ) + rerun = integrate_reaction_system( + dataset.system, + dataset.concentrations[:, 0], + dataset.times, + method="rk4", + substeps=dataset.card["simulator"]["substeps"], + ).permute(1, 0, 2) + assert torch.allclose(rerun, dataset.concentrations) + + +def test_long_runs_relax_to_card_equilibrium() -> None: + card = ReactionNetworkBenchmark.load_card("reaction_michaelis_menten_v1") + dataset = ReactionNetworkBenchmark.generate( + "reaction_michaelis_menten_v1", num_trajectories=1, num_timesteps=400 + ) + final = dataset.concentrations[0, -1] + # Boundary equilibrium: all substrate converted, enzyme freed. + assert float(final[1]) < 1e-6 and float(final[2]) < 1e-6 + spectrum = oracle_spectrum(dataset.system, final) + assert spectrum.num_conservation_laws == card["oracle"]["num_conservation_laws"] + + +def test_sequence_layout_and_encodings() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_michaelis_menten_v1", num_trajectories=2, num_timesteps=5 + ) + sequence = dataset.sequence(1) + assert isinstance(sequence, GraphSnapshotSequence) + assert sequence.num_timesteps == 5 + snapshot = sequence[0] + assert snapshot.x.shape == (4, 1) + assert snapshot.x.dtype == torch.float32 + assert snapshot.hyperedge_index.shape[0] == 2 # reaction_hyperedge card + assert dataset.transform is not None + expected = torch.log(dataset.concentrations[1, 0].clamp(min=1e-8)).to(torch.float32) + assert torch.allclose(snapshot.x.squeeze(-1), expected) + assert torch.allclose(sequence.timestamps, dataset.times) + raw = ReactionNetworkBenchmark.generate( + "reaction_brusselator_v1", num_trajectories=1, num_timesteps=3 + ) + assert raw.transform is None + assert torch.allclose( + raw.sequence(0)[2].x.squeeze(-1), + raw.concentrations[0, 2].to(torch.float32), + ) + assert not hasattr(raw.sequence(0)[0], "hyperedge_index") + with pytest.raises(IndexError, match="out of range"): + raw.sequence(5) + + +def test_cstr_controls_are_piecewise_constant_and_in_range() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_van_de_vusse_v1", num_trajectories=2, num_timesteps=25 + ) + assert dataset.controls is not None + assert dataset.controls.shape == (2, 25, 2) + spec = dataset.card["controls"] + for k, name in enumerate(spec["names"]): + lo, hi = spec["ranges"][name] + values = dataset.controls[..., k] + assert float(values.min()) >= lo and float(values.max()) <= hi + hold = spec["hold_steps"] + first = dataset.controls[0, :hold] + assert torch.allclose(first, first[:1].expand_as(first)) + assert not torch.allclose(dataset.controls[0, 0], dataset.controls[0, hold]) + sequence = dataset.sequence(0) + assert sequence.has_controls + assert sequence.control_inputs.shape == (25, 2) + + +def test_ssa_simulator_option() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_chain_v1", + num_trajectories=2, + num_timesteps=6, + simulator="ssa", + volume=200.0, + ) + assert dataset.simulator == "ssa" + assert dataset.volume == 200.0 + counts = dataset.concentrations * 200.0 + assert torch.allclose(counts, counts.round()) + repeat = ReactionNetworkBenchmark.generate( + "reaction_chain_v1", + num_trajectories=2, + num_timesteps=6, + simulator="ssa", + volume=200.0, + ) + assert torch.equal(dataset.concentrations, repeat.concentrations) + with pytest.raises(ValueError, match="not available for 'cstr'"): + ReactionNetworkBenchmark.generate("reaction_van_de_vusse_v1", simulator="ssa") + with pytest.raises(ValueError, match="volume must be"): + ReactionNetworkBenchmark.generate( + "reaction_chain_v1", num_trajectories=1, simulator="ssa", volume=0.0 + ) + + +def test_split_is_deterministic_and_covers_all() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_chain_v1", num_trajectories=10 + ) + split = ReactionNetworkBenchmark.split(dataset) + sizes = ( + split.train.num_trajectories, + split.val.num_trajectories, + split.test.num_trajectories, + ) + assert sizes == (7, 1, 2) + rebuilt = torch.cat( + [ + split.train.concentrations, + split.val.concentrations, + split.test.concentrations, + ] + ) + assert torch.equal(rebuilt, dataset.concentrations) + tiny = ReactionNetworkBenchmark.split(dataset.subset([0, 1, 2])) + assert tiny.train.num_trajectories == 1 + with pytest.raises(ValueError, match="at least 3 trajectories"): + ReactionNetworkBenchmark.split(dataset.subset([0, 1])) + with pytest.raises(ValueError, match="sum to 1"): + ReactionNetworkBenchmark.split(dataset, train_ratio=0.5) + with pytest.raises(ValueError, match="leave no training data"): + ReactionNetworkBenchmark.split( + dataset.subset([0, 1, 2]), train_ratio=0.0, val_ratio=0.5, test_ratio=0.5 + ) + + +def test_payload_round_trip(tmp_path: Path) -> None: + for dataset_id in ("reaction_van_de_vusse_v1", "reaction_chain_v1"): + dataset = ReactionNetworkBenchmark.generate( + dataset_id, num_trajectories=2, num_timesteps=6 + ) + path = dataset.save(tmp_path / "nested" / f"{dataset_id}.pt") + loaded = ReactionDataset.load(path) + assert loaded.dataset_id == dataset_id + assert loaded.card == dataset.card + assert loaded.system.to_dict() == dataset.system.to_dict() + assert torch.equal(loaded.concentrations, dataset.concentrations) + assert torch.equal(loaded.times, dataset.times) + if dataset.controls is None: + assert loaded.controls is None + else: + assert loaded.controls is not None + assert torch.equal(loaded.controls, dataset.controls) + assert (loaded.transform is None) == (dataset.transform is None) + assert loaded.seed == dataset.seed + assert torch.equal(loaded.sequence(0)[1].x, dataset.sequence(0)[1].x) + + +def test_payload_rejects_wrong_format_and_version(tmp_path: Path) -> None: + bad = tmp_path / "bad.pt" + torch.save({"format": "other"}, bad) + with pytest.raises(ValueError, match="is not a"): + ReactionDataset.load(bad) + dataset = ReactionNetworkBenchmark.generate( + "reaction_chain_v1", num_trajectories=1, num_timesteps=2 + ) + path = dataset.save(tmp_path / "ok.pt") + payload = torch.load(path, weights_only=True) + payload["format_version"] = PAYLOAD_FORMAT_VERSION + 1 + torch.save(payload, path) + with pytest.raises(ValueError, match="unsupported reaction dataset format_version"): + ReactionDataset.load(path) + + +@pytest.mark.parametrize( + ("kwargs", "match"), + [ + ({"num_trajectories": 0}, "num_trajectories must be"), + ({"num_timesteps": 1}, "num_timesteps must be"), + ({"seed": -1}, "seed must be"), + ({"simulator": "euler"}, "simulator must be"), + ], +) +def test_generate_validation(kwargs: dict[str, object], match: str) -> None: + with pytest.raises(ValueError, match=match): + ReactionNetworkBenchmark.generate("reaction_chain_v1", **kwargs) # type: ignore[arg-type] + + +def test_card_validation_helpers() -> None: + from koopman_graph.datasets.reaction.benchmark import ( + _check_ranges, + system_from_card, + ) + + with pytest.raises(ValueError, match="do not match species"): + _check_ranges({"A": [0.0, 1.0]}, ("A", "B")) + with pytest.raises(ValueError, match="0 <= low <= high"): + _check_ranges({"A": [1.0, 0.0]}, ("A",)) + with pytest.raises(ValueError, match="unknown generator kind"): + system_from_card({"generator": {"kind": "magic"}}) + with pytest.raises(ValueError, match="unknown catalog system"): + system_from_card({"generator": {"kind": "catalog", "name": "nope"}}) + + +def test_card_builder_script_check_is_clean() -> None: + """Committed cards equal what ``scripts/build_reaction_cards.py`` writes.""" + if importlib.util.find_spec("cantera") is None: + pytest.skip("card builder needs the [cantera] extra for the h2o2 card") + result = subprocess.run( + [ + sys.executable, + str(REPO_ROOT / "scripts" / "build_reaction_cards.py"), + "--check", + ], + capture_output=True, + text=True, + check=False, + ) + assert result.returncode == 0, result.stderr diff --git a/tests/datasets/test_reaction_import_isolation.py b/tests/datasets/test_reaction_import_isolation.py index d1ba6cf..3ddfbfc 100644 --- a/tests/datasets/test_reaction_import_isolation.py +++ b/tests/datasets/test_reaction_import_isolation.py @@ -40,3 +40,12 @@ def test_root_all_does_not_export_reaction_names() -> None: def test_cantera_marker_runs_only_with_extra() -> None: """``tests/conftest.py`` skips this unless ``[cantera]`` is installed.""" assert importlib.util.find_spec("cantera") is not None + + +def test_datasets_exports_only_the_benchmark() -> None: + import koopman_graph.datasets as datasets + from koopman_graph.datasets.reaction import ReactionNetworkBenchmark + + assert datasets.ReactionNetworkBenchmark is ReactionNetworkBenchmark + assert "ReactionNetworkBenchmark" in datasets.__all__ + assert "ReactionSystem" not in datasets.__all__ diff --git a/tests/datasets/test_reaction_oracles.py b/tests/datasets/test_reaction_oracles.py index 9e638d0..6b715b0 100644 --- a/tests/datasets/test_reaction_oracles.py +++ b/tests/datasets/test_reaction_oracles.py @@ -206,3 +206,17 @@ def test_equilibrium_rejects_one_way_boundary_point() -> None: with pytest.raises(RuntimeError, match="non-stationary point"): _require_stationary(system, torch.tensor([1.0, 0.0, 1e-9], dtype=_F64)) + + +def test_equilibrium_accepts_vanishing_leftovers_on_boundary() -> None: + """Draining onto a zero-flux boundary is stationary despite tiny leftovers.""" + from koopman_graph.datasets.reaction import sample_random_network + + system = sample_random_network(seed=0).system + c0 = torch.full((6,), 0.55, dtype=_F64) + relaxed = integrate_rk4(system, c0, torch.tensor([0.0, 200.0]), substeps=4000)[-1] + c_eq = equilibrium(system, c0, guess=relaxed) + assert float(system.time_derivative(c_eq).abs().max()) < 1e-12 + assert system.element_matrix is not None + drift = system.element_matrix @ (c_eq - c0) + assert float(drift.abs().max()) < 1e-10 From 6e68ae0c5f219591e8ba06f057293346b1479be6 Mon Sep 17 00:00:00 2001 From: tjkessler Date: Thu, 15 Oct 2026 16:25:37 -0400 Subject: [PATCH 3/5] Address review findings on ReactionNetworkBenchmark - build_reaction_cards.py --check compares card numbers with the test tolerances instead of byte-exact text, so platform/BLAS last-digit differences no longer fail CI, and flags reaction_*.json files that no spec produces. - system_from_card (and so generate / compute_card_oracle) refuses a Cantera card whose recorded mechanism SHA-256 differs from the installed file. - eigenvalue_timescales treats only round-off-level real parts (100 * eps * max|lambda|) as zero, so slow modes in stiff spectra stay finite. - split rejects zero ratios instead of silently assigning one trajectory. - SSA path seeds are drawn from the dataset generator (no aliasing across dataset seeds). - ReactionDataset.features bounds-checks its index like sequence. - compute_card_oracle builds the CSTR reactor once. - CHANGELOG: move the extras / marker / CI / import-hygiene entries back under Added. - Tests: equilibrium raises on a forced Newton stall, stiff-spectrum timescales, tolerant card comparison, orphan cards, SHA mismatch, index checks, split zero ratios, SSA seed independence. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 23 +++-- scripts/build_reaction_cards.py | 68 +++++++++++++- .../datasets/reaction/benchmark.py | 74 +++++++++++---- .../datasets/reaction/oracles.py | 16 ++-- tests/datasets/test_reaction_benchmark.py | 93 ++++++++++++++++++- tests/datasets/test_reaction_oracles.py | 32 ++++++- 6 files changed, 265 insertions(+), 41 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d1be36..d4c63a4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -80,6 +80,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ideal-gas reactors, adiabatic or isothermal. - Also: `require_cantera`, `resolve_mechanism_path`, and `mechanism_sha256`. +- `[reaction]` extra (`scipy>=1.11.0`) for stiff reaction-network solvers + and `[cantera]` extra (`cantera>=3.2,<3.3`) for Cantera mechanisms. +- `cantera` pytest marker; marked tests are skipped automatically when the + `[cantera]` extra is not installed. +- CI test jobs and the release verify job install `[reaction]` and + `[cantera]`, so `cantera`-marked tests run there. +- `tests/test_import_hygiene.py`: `import koopman_graph` must not load + SciPy, Cantera, or torchdiffeq. - `ReactionNetworkBenchmark` (exported from `koopman_graph.datasets` and `koopman_graph.datasets.reaction`): card-defined, seeded reaction-network datasets. @@ -100,8 +108,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `compute_card_oracle` and `scripts/build_reaction_cards.py` (with `--check`) regenerate the oracle sections; tests recompute them to catch drift. -- `eigenvalue_timescales`: `-1/Re(λ)` with numerically zero real parts - (relative to `max|λ|`) mapped to `inf`; `OracleSpectrum.timescales` uses it. +- `eigenvalue_timescales`: `-1/Re(λ)` with real parts below eigensolver + round-off (`100 · eps · max|λ|`) mapped to `inf`; + `OracleSpectrum.timescales` uses it. `ReactionNetworkBenchmark.generate` + refuses a Cantera card whose recorded mechanism SHA-256 differs from the + installed file. ### Changed @@ -119,14 +130,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 1e-12 · (peak initial turnover)` and otherwise raises `RuntimeError`. Previously the converted Cantera `h2o2` system could report a near-boundary point with one-way flux as its equilibrium. -- `[reaction]` extra (`scipy>=1.11.0`) for stiff reaction-network solvers - and `[cantera]` extra (`cantera>=3.2,<3.3`) for Cantera mechanisms. -- `cantera` pytest marker; marked tests are skipped automatically when the - `[cantera]` extra is not installed. -- CI test jobs and the release verify job install `[reaction]` and - `[cantera]`, so `cantera`-marked tests run there. -- `tests/test_import_hygiene.py`: `import koopman_graph` must not load - SciPy, Cantera, or torchdiffeq. ## [0.15.0] - 2026-08-19 diff --git a/scripts/build_reaction_cards.py b/scripts/build_reaction_cards.py index 567e816..4a596c5 100644 --- a/scripts/build_reaction_cards.py +++ b/scripts/build_reaction_cards.py @@ -13,12 +13,18 @@ python scripts/build_reaction_cards.py # write all cards python scripts/build_reaction_cards.py --check # exit 1 if any differ + +``--check`` compares numbers with a tolerance (oracle values come from +eigensolvers, Newton, RK4, and Cantera, whose last digits vary across +platforms and BLAS builds) and also flags ``reaction_*.json`` files in the +data directory that no spec produces. """ from __future__ import annotations import argparse import json +import math import sys from pathlib import Path from typing import Any @@ -253,6 +259,45 @@ def build_card(spec: dict[str, Any]) -> dict[str, Any]: return card +# Tolerances for ``--check``; they match tests/datasets/test_reaction_benchmark.py. +_RTOL, _ATOL = 1e-8, 1e-12 +_CANTERA_RTOL, _CANTERA_ATOL = 1e-6, 1e-18 + + +def cards_match(a: Any, b: Any, *, rtol: float, atol: float) -> bool: + """Compare two parsed cards, with a tolerance on floating-point leaves. + + Parameters + ---------- + a, b : Any + Parsed JSON values. + rtol, atol : float + Relative and absolute tolerance for floats. + + Returns + ------- + bool + ``True`` when the structures are equal and every float pair is close. + """ + if isinstance(a, dict) and isinstance(b, dict): + return a.keys() == b.keys() and all( + cards_match(a[k], b[k], rtol=rtol, atol=atol) for k in a + ) + if isinstance(a, list) and isinstance(b, list): + return len(a) == len(b) and all( + cards_match(x, y, rtol=rtol, atol=atol) for x, y in zip(a, b, strict=True) + ) + numeric = (int, float) + if ( + isinstance(a, numeric) + and isinstance(b, numeric) + and not isinstance(a, bool) + and not isinstance(b, bool) + ): + return math.isclose(a, b, rel_tol=rtol, abs_tol=atol) + return type(a) is type(b) and a == b + + def main(argv: list[str] | None = None) -> int: """Write (or check) every card. @@ -272,18 +317,31 @@ def main(argv: list[str] | None = None) -> int: DATA_DIR.mkdir(parents=True, exist_ok=True) stale = [] for spec in SPECS: - text = json.dumps(build_card(spec), indent=2) + "\n" + card = build_card(spec) path = DATA_DIR / f"{spec['dataset_id']}.json" if args.check: - if not path.is_file() or path.read_text(encoding="utf-8") != text: + cantera = card["cantera"] is not None + if not path.is_file() or not cards_match( + json.loads(path.read_text(encoding="utf-8")), + card, + rtol=_CANTERA_RTOL if cantera else _RTOL, + atol=_CANTERA_ATOL if cantera else _ATOL, + ): stale.append(path.name) else: - path.write_text(text, encoding="utf-8") + path.write_text(json.dumps(card, indent=2) + "\n", encoding="utf-8") print(f"wrote {path.relative_to(DATA_DIR.parents[3])}") + expected = {f"{spec['dataset_id']}.json" for spec in SPECS} + orphans = sorted( + path.name + for path in DATA_DIR.glob("reaction_*.json") + if path.name not in expected + ) if stale: print("stale cards: " + ", ".join(stale), file=sys.stderr) - return 1 - return 0 + if orphans: + print("cards without a spec: " + ", ".join(orphans), file=sys.stderr) + return 1 if stale or orphans else 0 if __name__ == "__main__": diff --git a/src/koopman_graph/datasets/reaction/benchmark.py b/src/koopman_graph/datasets/reaction/benchmark.py index 9bee9f3..d8b3fc9 100644 --- a/src/koopman_graph/datasets/reaction/benchmark.py +++ b/src/koopman_graph/datasets/reaction/benchmark.py @@ -173,7 +173,13 @@ def features(self, index: int, *, dtype: torch.dtype = torch.float32) -> Tensor: Tensor Transformed concentrations when :attr:`transform` is set, raw concentrations otherwise. + + Raises + ------ + IndexError + If ``index`` is out of range. """ + self._check_index(index) c = self.concentrations[index] values = c if self.transform is None else self.transform.forward(c) return values.unsqueeze(-1).to(dtype) @@ -204,9 +210,7 @@ def sequence( IndexError If ``index`` is out of range. """ - if not 0 <= index < self.num_trajectories: - msg = f"index {index} out of range for {self.num_trajectories} trajectories" - raise IndexError(msg) + self._check_index(index) edge_index, edge_weight = species_influence_graph(self.system) hyperedge_index = ( reaction_hyperedge_index(self.system) @@ -227,6 +231,17 @@ def sequence( control_inputs=controls, ) + def _check_index(self, index: int) -> None: + """Raise ``IndexError`` unless ``0 <= index < num_trajectories``. + + Parameters + ---------- + index + See signature.""" + if not 0 <= index < self.num_trajectories: + msg = f"index {index} out of range for {self.num_trajectories} trajectories" + raise IndexError(msg) + def subset(self, indices: list[int] | Tensor) -> ReactionDataset: """Return a dataset restricted to the given trajectory indices. @@ -459,6 +474,9 @@ def generate( If arguments are invalid or ``simulator="ssa"`` does not apply. ImportError For the Cantera card when the ``[cantera]`` extra is missing. + RuntimeError + For the Cantera card when the installed mechanism's SHA-256 does + not match the card. """ card = cls.load_card(dataset_id) n_traj = _count( @@ -498,7 +516,7 @@ def generate( if simulator == "ssa": used_volume = float(volume) concentrations = _simulate_ssa( - system, initial, times, used_volume, resolved_seed + system, initial, times, used_volume, generator ) elif kind == "cstr": concentrations, controls = _simulate_cstr( @@ -547,7 +565,8 @@ def split( dataset : ReactionDataset Dataset with at least three trajectories. train_ratio, val_ratio, test_ratio : float, optional - Fractions summing to ``1``. Defaults are ``0.7 / 0.1 / 0.2``. + Positive fractions summing to ``1``. Defaults are + ``0.7 / 0.1 / 0.2``. Returns ------- @@ -561,10 +580,10 @@ def split( """ ratios = (train_ratio, val_ratio, test_ratio) if ( - any(r < 0 or not math.isfinite(r) for r in ratios) + any(r <= 0 or not math.isfinite(r) for r in ratios) or abs(sum(ratios) - 1.0) > 1e-6 ): - msg = f"ratios must be non-negative and sum to 1, got {ratios}" + msg = f"ratios must be positive and sum to 1, got {ratios}" raise ValueError(msg) n = dataset.num_trajectories if n < 3: @@ -600,6 +619,9 @@ def system_from_card(card: dict[str, Any]) -> ReactionSystem: ------ ValueError If the generator kind or catalog name is unknown. + RuntimeError + For a Cantera card whose recorded mechanism SHA-256 differs from the + installed mechanism file. """ generator = card["generator"] kind = generator["kind"] @@ -615,9 +637,22 @@ def system_from_card(card: dict[str, Any]) -> ReactionSystem: if kind == "cstr": return van_de_vusse_cstr(**params).system if kind == "cantera": - from koopman_graph.datasets.reaction.cantera_io import load_cantera_mechanism + from koopman_graph.datasets.reaction.cantera_io import ( + load_cantera_mechanism, + mechanism_sha256, + resolve_mechanism_path, + ) cantera = card["cantera"] + installed = mechanism_sha256(resolve_mechanism_path(cantera["mechanism"])) + if installed != cantera["sha256"]: + msg = ( + f"installed {cantera['mechanism']} has SHA-256 {installed}, but " + f"card {card.get('dataset_id', '?')!r} was built from " + f"{cantera['sha256']} (Cantera {cantera['cantera_version']}); " + "the generated data would not match the card" + ) + raise RuntimeError(msg) return load_cantera_mechanism( cantera["mechanism"], temperature_k=cantera["temperature_k"], @@ -649,23 +684,24 @@ def compute_card_oracle(card: dict[str, Any]) -> dict[str, Any]: steady state. """ kind = card["generator"]["kind"] - system = system_from_card(card) if kind == "cstr": controls = card["controls"]["ranges"] temperature = sum(controls["temperature_k"]) / 2.0 dilution = sum(controls["dilution_rate"]) / 2.0 reactor = van_de_vusse_cstr(**card["generator"].get("params", {})) + num_species = reactor.system.num_species state = reactor.steady_state(temperature_k=temperature, dilution_rate=dilution) jac = mass_action_jacobian(reactor.system.at_temperature(temperature), state) - jac = jac - dilution * torch.eye(system.num_species, dtype=torch.float64) + jac = jac - dilution * torch.eye(num_species, dtype=torch.float64) eigenvalues = torch.linalg.eigvals(jac) eigenvalues = eigenvalues[torch.argsort(eigenvalues.real, descending=True)] return _oracle_dict( "cstr_steady_state_at_control_midpoint", state, eigenvalues, - torch.zeros(0, system.num_species, dtype=torch.float64), + torch.zeros(0, num_species, dtype=torch.float64), ) + system = system_from_card(card) if kind == "cantera": from koopman_graph.datasets.reaction.cantera_io import require_cantera @@ -820,10 +856,13 @@ def _simulate_ssa( initial: Tensor, times: Tensor, volume: float, - seed: int, + generator: torch.Generator, ) -> Tensor: """Run one Gillespie path per initial state and sample it on ``times``. + Each path gets its own seed drawn from ``generator``, so paths are + independent within and across dataset seeds. + Parameters ---------- system @@ -834,7 +873,7 @@ def _simulate_ssa( See signature. volume See signature. - seed + generator See signature. Returns @@ -844,12 +883,13 @@ def _simulate_ssa( msg = f"volume must be finite and > 0, got {volume!r}" raise ValueError(msg) t_end = float(times[-1]) + path_seeds = torch.randint( + 0, 2**62, (initial.shape[0],), generator=generator, dtype=torch.long + ).tolist() paths = [] - for i, c0 in enumerate(initial): + for c0, path_seed in zip(initial, path_seeds, strict=True): counts = torch.round(c0 * volume).to(torch.long) - path = gillespie_ssa( - system, counts, t_end, volume=volume, seed=seed * 100_003 + i - ) + path = gillespie_ssa(system, counts, t_end, volume=volume, seed=path_seed) paths.append(path.concentrations_on_grid(times)) return torch.stack(paths) diff --git a/src/koopman_graph/datasets/reaction/oracles.py b/src/koopman_graph/datasets/reaction/oracles.py index 00bcda5..a558dd4 100644 --- a/src/koopman_graph/datasets/reaction/oracles.py +++ b/src/koopman_graph/datasets/reaction/oracles.py @@ -68,9 +68,9 @@ def timescales(self) -> Tensor: Tensor ``float64`` tensor aligned with :attr:`dynamic_eigenvalues`: positive for decaying modes, ``inf`` when ``Re(λ)`` is zero to - within ``1e-12 · max|λ|`` (a neutral centre, or a species - drained onto a zero-flux boundary), and negative for growing - modes. + within eigensolver round-off, ``100 · eps · max|λ|`` (a neutral + centre, or a species drained onto a zero-flux boundary), and + negative for growing modes. """ return eigenvalue_timescales(self.dynamic_eigenvalues) @@ -112,7 +112,11 @@ def discrete_eigenvalues(self, time_step: float) -> Tensor: return torch.exp(self.dynamic_eigenvalues * dt) -_ZERO_EIGENVALUE_RTOL = 1e-12 +# Eigenvalues are only resolved to about ``eps · ‖J‖``; real parts below a +# small multiple of that are numerical zeros. Keeping the cutoff at round-off +# level preserves genuinely slow modes in stiff mechanisms, whose eigenvalue +# magnitudes can span ten or more decades. +_ZERO_EIGENVALUE_RTOL = 100.0 * torch.finfo(torch.float64).eps def eigenvalue_timescales(eigenvalues: Tensor) -> Tensor: @@ -126,8 +130,8 @@ def eigenvalue_timescales(eigenvalues: Tensor) -> Tensor: Returns ------- Tensor - ``float64`` timescales ``(k,)``. ``|Re(λ)| <= 1e-12 · max|λ|`` maps - to ``inf``. + ``float64`` timescales ``(k,)``. ``|Re(λ)| <= 100 · eps · max|λ|`` + (about ``2.2e-14 · max|λ|``) maps to ``inf``. """ real = eigenvalues.real.to(torch.float64) scale = float(eigenvalues.abs().max()) if eigenvalues.numel() else 0.0 diff --git a/tests/datasets/test_reaction_benchmark.py b/tests/datasets/test_reaction_benchmark.py index 52120a5..ddb1fe1 100644 --- a/tests/datasets/test_reaction_benchmark.py +++ b/tests/datasets/test_reaction_benchmark.py @@ -3,6 +3,7 @@ from __future__ import annotations import importlib.util +import json import math import subprocess import sys @@ -337,7 +338,57 @@ def test_split_is_deterministic_and_covers_all() -> None: ReactionNetworkBenchmark.split(dataset, train_ratio=0.5) with pytest.raises(ValueError, match="leave no training data"): ReactionNetworkBenchmark.split( - dataset.subset([0, 1, 2]), train_ratio=0.0, val_ratio=0.5, test_ratio=0.5 + dataset.subset([0, 1, 2, 3]), + train_ratio=0.1, + val_ratio=0.45, + test_ratio=0.45, + ) + # A zero ratio would be silently raised to one trajectory; reject it. + for ratios in ((0.8, 0.0, 0.2), (0.8, 0.2, 0.0), (0.0, 0.5, 0.5)): + with pytest.raises(ValueError, match="positive and sum to 1"): + ReactionNetworkBenchmark.split( + dataset, + train_ratio=ratios[0], + val_ratio=ratios[1], + test_ratio=ratios[2], + ) + + +def test_features_and_sequence_reject_out_of_range_index() -> None: + dataset = ReactionNetworkBenchmark.generate( + "reaction_chain_v1", num_trajectories=2, num_timesteps=3 + ) + for index in (-1, 2): + with pytest.raises(IndexError, match="out of range"): + dataset.features(index) + with pytest.raises(IndexError, match="out of range"): + dataset.sequence(index) + + +def test_ssa_paths_differ_across_seeds_and_trajectories() -> None: + def ssa(seed: int) -> torch.Tensor: + return ReactionNetworkBenchmark.generate( + "reaction_chain_v1", + num_trajectories=2, + num_timesteps=8, + seed=seed, + simulator="ssa", + volume=50.0, + ).concentrations + + assert not torch.equal(ssa(0), ssa(1)) + + +def test_cantera_card_rejects_mechanism_sha_mismatch( + monkeypatch: pytest.MonkeyPatch, +) -> None: + from koopman_graph.datasets.reaction import cantera_io + + monkeypatch.setattr(cantera_io, "resolve_mechanism_path", lambda name: Path(name)) + monkeypatch.setattr(cantera_io, "mechanism_sha256", lambda path: "0" * 64) + with pytest.raises(RuntimeError, match="would not match the card"): + ReactionNetworkBenchmark.generate( + "reaction_cantera_h2o2_v1", num_trajectories=1, num_timesteps=2 ) @@ -409,6 +460,46 @@ def test_card_validation_helpers() -> None: system_from_card({"generator": {"kind": "catalog", "name": "nope"}}) +def _load_card_builder() -> object: + """Import ``scripts/build_reaction_cards.py`` as a module.""" + path = REPO_ROOT / "scripts" / "build_reaction_cards.py" + spec = importlib.util.spec_from_file_location("build_reaction_cards", path) + assert spec is not None and spec.loader is not None + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def test_card_builder_compares_floats_with_tolerance() -> None: + builder = _load_card_builder() + card = {"oracle": {"x": [1.0, -0.5], "label": "a"}, "n": 3, "flag": True} + near = {"oracle": {"x": [1.0 + 1e-12, -0.5], "label": "a"}, "n": 3, "flag": True} + assert builder.cards_match(card, near, rtol=1e-8, atol=1e-12) + far = {"oracle": {"x": [1.001, -0.5], "label": "a"}, "n": 3, "flag": True} + assert not builder.cards_match(card, far, rtol=1e-8, atol=1e-12) + relabeled = {"oracle": {"x": [1.0, -0.5], "label": "b"}, "n": 3, "flag": True} + assert not builder.cards_match(card, relabeled, rtol=1e-8, atol=1e-12) + assert not builder.cards_match({"a": True}, {"a": 1}, rtol=1e-8, atol=1e-12) + assert not builder.cards_match([1.0], [1.0, 2.0], rtol=1e-8, atol=1e-12) + + +def test_card_builder_check_flags_orphan_cards( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch, capsys: pytest.CaptureFixture +) -> None: + builder = _load_card_builder() + spec = next(s for s in builder.SPECS if s["dataset_id"] == "reaction_chain_v1") + monkeypatch.setattr(builder, "SPECS", [spec]) + monkeypatch.setattr(builder, "DATA_DIR", tmp_path) + assert builder.main(["--check"]) == 1 # card missing + (tmp_path / "reaction_chain_v1.json").write_text( + json.dumps(builder.build_card(spec), indent=2), encoding="utf-8" + ) + assert builder.main(["--check"]) == 0 + (tmp_path / "reaction_old_v1.json").write_text("{}", encoding="utf-8") + assert builder.main(["--check"]) == 1 + assert "reaction_old_v1.json" in capsys.readouterr().err + + def test_card_builder_script_check_is_clean() -> None: """Committed cards equal what ``scripts/build_reaction_cards.py`` writes.""" if importlib.util.find_spec("cantera") is None: diff --git a/tests/datasets/test_reaction_oracles.py b/tests/datasets/test_reaction_oracles.py index 6b715b0..47f186a 100644 --- a/tests/datasets/test_reaction_oracles.py +++ b/tests/datasets/test_reaction_oracles.py @@ -187,9 +187,11 @@ def test_equilibrium_validation() -> None: spectrum.discrete_eigenvalues(dt) -def test_equilibrium_rejects_one_way_boundary_point() -> None: +def test_equilibrium_rejects_one_way_boundary_point( + monkeypatch: pytest.MonkeyPatch, +) -> None: """A clamped Newton stall where flux runs one way is not a steady state.""" - # A -> B with a slow B -> A: from (1, 0) the only steady state is interior. + # A + B -> C with a slow C -> A + B: the only steady state is interior. system = ReactionSystem( species=("A", "B", "C"), reactant_stoichiometry=torch.tensor( @@ -207,6 +209,32 @@ def test_equilibrium_rejects_one_way_boundary_point() -> None: with pytest.raises(RuntimeError, match="non-stationary point"): _require_stationary(system, torch.tensor([1.0, 0.0, 1e-9], dtype=_F64)) + # Force Newton to stall where it starts (zero Jacobian -> zero step): the + # step test passes immediately, and equilibrium must still refuse the point. + from koopman_graph.datasets.reaction import oracles + + monkeypatch.setattr( + oracles, + "mass_action_jacobian", + lambda system, c: torch.zeros(3, 3, dtype=_F64), + ) + with pytest.raises(RuntimeError, match="non-stationary point"): + equilibrium(system, torch.tensor([1.0, 1.0, 0.0], dtype=_F64)) + + +def test_eigenvalue_timescales_keep_slow_modes_in_stiff_spectra() -> None: + """Only round-off-level real parts become infinite, not slow stiff modes.""" + from koopman_graph.datasets.reaction import eigenvalue_timescales + + # max|λ| = 1e10, so the round-off cutoff is about 2.2e-4. + eigenvalues = torch.tensor([-1e10, -1e-2, -1e-3, 1e-6], dtype=_F64) + timescales = eigenvalue_timescales(eigenvalues) + assert torch.allclose( + timescales[:3], torch.tensor([1e-10, 1e2, 1e3], dtype=_F64), rtol=1e-12 + ) + assert math.isinf(float(timescales[3])) + assert eigenvalue_timescales(torch.zeros(0, dtype=_F64)).shape == (0,) + def test_equilibrium_accepts_vanishing_leftovers_on_boundary() -> None: """Draining onto a zero-flux boundary is stationary despite tiny leftovers.""" From f20605baade8a6e2c6f71b1aa9f9f714d2375400 Mon Sep 17 00:00:00 2001 From: tjkessler Date: Thu, 15 Oct 2026 17:09:40 -0400 Subject: [PATCH 4/5] Make conservation-law basis canonical so cards reproduce on Linux Ubuntu CI failed test_card_oracle_matches_recomputation (Michaelis-Menten, random network), test_cantera_card_oracle_and_provenance, and the card builder --check: oracle['conservation_laws'] drifted. With two or more conservation laws any rotation of an orthonormal basis is valid, and torch.linalg.svd returns a different one under Linux LAPACK than under macOS Accelerate, so the committed cards held a platform-specific basis. The single-law chain card was unaffected. _stoichiometric_basis now canonicalizes the conservation laws: Gram-Schmidt in species order on the columns of the projector P = W^T W, which is basis independent, so every platform gets the same rows (to rounding) with a fixed sign. Cards regenerated; apart from the rotated laws, oracle values move only at rounding level (<= 3e-14). Adds a test that a rotated, sign-flipped basis canonicalizes back to the same rows. Also merges main (#70-#73) into the branch, which GitHub reported as behind. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 7 +++ .../data/reaction_cantera_h2o2_v1.json | 40 +++++++-------- .../reaction/data/reaction_chain_v1.json | 8 +-- .../data/reaction_michaelis_menten_v1.json | 28 +++++------ .../data/reaction_random_network_v1.json | 46 ++++++++--------- .../datasets/reaction/oracles.py | 49 ++++++++++++++++++- tests/datasets/test_reaction_oracles.py | 21 ++++++++ 7 files changed, 136 insertions(+), 63 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 62221fc..27195fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -132,6 +132,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- `conservation_laws` returns a canonical basis (Gram-Schmidt in species + order on the basis-independent projector) instead of whatever an SVD + returns. With two or more conservation laws the SVD basis differed + between LAPACK builds, so the Michaelis-Menten, random-network, and + Cantera cards did not reproduce on Linux CI. Card oracle values change + only at rounding level (<= 3e-14) apart from the rotated laws. + - `equilibrium` no longer returns a clamped Newton stall that is not stationary. It now requires `|dc_s/dt| <= 1e-6 · gross turnover_s + 1e-12 · (peak initial turnover)` and otherwise raises `RuntimeError`. diff --git a/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json index ce0fe23..4461e1a 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json @@ -91,26 +91,26 @@ "num_conservation_laws": 4, "conservation_laws": [ [ - 0.28064053942470113, - 0.1403202697123506, - -0.3167114270242802, - -0.6334228540485629, - -0.17639115731193028, - -0.036070887599579333, - -0.493102584336212, - -0.3527823146238605, + 0.6454972243679027, + 0.32274861218395146, + -0.19364916731037016, + -0.3872983346207416, + 0.12909944487358088, + 0.45184805705753206, + -0.06454972243679041, + 0.2581988897471615, 0.0, 0.0 ], [ - 0.5812981629921767, - 0.2906490814960885, - -0.062133233127772275, - -0.12426646625554473, - 0.2285158483683161, - 0.5191649298644041, - 0.16638261524054332, - 0.4570316967366318, + 1.0749688212373174e-16, + 5.374844106186587e-17, + 0.25819888974716043, + 0.5163977794943231, + 0.25819888974716104, + 0.25819888974716054, + 0.5163977794943226, + 0.5163977794943219, 0.0, 0.0 ], @@ -123,8 +123,8 @@ 0.0, 0.0, 0.0, - 0.0, - 1.0 + 1.0, + 0.0 ], [ 0.0, @@ -135,8 +135,8 @@ 0.0, 0.0, 0.0, - 1.0, - 0.0 + 0.0, + 1.0 ] ], "dynamic_eigenvalues_real": [ diff --git a/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json index b4826d1..61eacf6 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_chain_v1.json @@ -62,10 +62,10 @@ "oracle": { "state": "equilibrium_from_range_midpoints", "reference_state": [ - -3.514983438626458e-17, - -8.784035482230695e-17, - -1.378654434722274e-16, - 1.300000000000001 + 0.0, + 0.0, + 0.0, + 1.3000000000000007 ], "num_conservation_laws": 1, "conservation_laws": [ diff --git a/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json index e728cbf..295fa86 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json @@ -62,37 +62,37 @@ "oracle": { "state": "equilibrium_from_range_midpoints", "reference_state": [ - 0.12499999999999997, + 0.12499999999999993, 0.0, 0.0, - 1.0750000000000002 + 1.075 ], "num_conservation_laws": 2, "conservation_laws": [ [ - -0.7123282872413622, - 0.010643277802640164, - -0.701685009438722, - 0.01064327780264014 + 0.7745966692414834, + -0.258198889747161, + 0.5163977794943224, + -0.25819888974716093 ], [ - -0.3042834389114655, - 0.6323659704930491, - 0.32808253158158357, - 0.632365970493049 + -1.9229626863835636e-16, + 0.5773502691896258, + 0.5773502691896255, + 0.5773502691896257 ] ], "dynamic_eigenvalues_real": [ - -0.25, - -2.5 + -0.2499999999999999, + -2.499999999999999 ], "dynamic_eigenvalues_imag": [ 0.0, 0.0 ], "timescales": [ - 4.0, - 0.4 + 4.000000000000002, + 0.40000000000000013 ] } } diff --git a/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json index 422a223..8709ec9 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json @@ -83,37 +83,37 @@ "oracle": { "state": "equilibrium_from_range_midpoints_after_relaxation", "reference_state": [ - 0.9371973830410871, - 2.668598691520542, - 7.203240409913488e-32, + 0.937197383041087, + 2.6685986915205424, + 0.0, 1.181401308479456, - 2.8981443113177526e-30, - 1.6307572951339986e-30 + 1.0829309696332013e-30, + 0.0 ], "num_conservation_laws": 2, "conservation_laws": [ [ - 0.2889695340165011, - -0.4048248518165172, - -0.11585531780001589, - 0.1731142162164851, - -0.23171063560003166, - -0.8096497036330341 + 0.42465029006520066, + -0.2702320027687639, + 0.1544182872964367, + 0.5790685773616372, + 0.3088365745928733, + -0.5404640055375279 ], [ - 0.31116631768677233, - 0.00715999869017006, - 0.3183263163769422, - 0.6294926340637144, - 0.6366526327538841, - 0.014319997380339834 + 4.60274150786893e-17, + 0.30151134457776385, + 0.3015113445777636, + 0.30151134457776363, + 0.6030226891555269, + 0.6030226891555275 ] ], "dynamic_eigenvalues_real": [ - -1.1843321716463685e-17, - -0.8165944002656043, - -0.8249568943108102, - -30.345958028796638 + 0.0, + -0.8165944002656039, + -0.8249568943108101, + -30.34595802879661 ], "dynamic_eigenvalues_imag": [ 0.0, @@ -123,9 +123,9 @@ ], "timescales": [ null, - 1.2245981599613487, + 1.2245981599613491, 1.2121845479398354, - 0.032953317837290066 + 0.0329533178372901 ] } } diff --git a/src/koopman_graph/datasets/reaction/oracles.py b/src/koopman_graph/datasets/reaction/oracles.py index 9a956a4..7bf9c55 100644 --- a/src/koopman_graph/datasets/reaction/oracles.py +++ b/src/koopman_graph/datasets/reaction/oracles.py @@ -193,6 +193,18 @@ def conservation_laws(system: ReactionSystem) -> Tensor: num_species)`` and orthonormal rows, ``W N = 0``. It has zero rows when nothing is conserved (e.g. an open system with independent sources). + + Notes + ----- + With two or more laws any rotation of an orthonormal basis is equally + valid, and the one an SVD returns depends on the LAPACK build. The rows + here are canonical instead: Gram-Schmidt in species order on the + columns of the projector ``P = WᵀW``, which does not depend on the + starting basis. Row ``k`` is the normalized component of the ``k``-th + independent species unit vector inside the conserved subspace, with a + positive entry for that species. The same network therefore gives the + same ``W`` (to rounding) on every platform, which keeps dataset cards + reproducible. """ return _stoichiometric_basis(system)[2] @@ -448,10 +460,43 @@ def _stoichiometric_basis(system: ReactionSystem) -> tuple[int, Tensor, Tensor]: ------- ``(rank, q, w)`` with ``q`` of shape ``(num_species, rank)`` and the conservation laws ``w`` of shape ``(num_species - rank, num_species)``, - from one SVD.""" + from one SVD; ``w`` is canonicalized (see :func:`conservation_laws`).""" _require_system(system) u, _vh, rank = _svd_rank(system.stoichiometric_matrix) - return rank, u[:, :rank], u[:, rank:].T.contiguous() + return rank, u[:, :rank], _canonical_row_basis(u[:, rank:].T.contiguous()) + + +_CANONICAL_PIVOT_TOL = 1e-6 + + +def _canonical_row_basis(w: Tensor) -> Tensor: + """Return a basis-independent orthonormal basis of ``rowspace(w)``. + + Parameters + ---------- + w + See signature. + + Returns + ------- + Orthonormal rows spanning the same subspace, chosen by Gram-Schmidt in + coordinate order on the projector ``wᵀw`` (independent of which + orthonormal ``w`` was passed).""" + num_laws, num_species = int(w.shape[0]), int(w.shape[1]) + if num_laws == 0: + return w + projector = w.T @ w # basis independent + basis: list[Tensor] = [] + for j in range(num_species): + v = projector[:, j].clone() + for b in basis: + v = v - (b @ v) * b + norm = float(v.norm()) + if norm > _CANONICAL_PIVOT_TOL: + basis.append(v / norm) + if len(basis) == num_laws: + break + return torch.stack(basis) def _svd_rank(matrix: Tensor) -> tuple[Tensor, Tensor, int]: diff --git a/tests/datasets/test_reaction_oracles.py b/tests/datasets/test_reaction_oracles.py index 6e5e9a9..bfb3aca 100644 --- a/tests/datasets/test_reaction_oracles.py +++ b/tests/datasets/test_reaction_oracles.py @@ -310,3 +310,24 @@ def test_equilibrium_accepts_vanishing_leftovers_on_boundary() -> None: assert system.element_matrix is not None drift = system.element_matrix @ (c_eq - c0) assert float(drift.abs().max()) < 1e-10 + + +def test_conservation_laws_are_canonical_across_bases() -> None: + """The basis does not depend on which orthonormal basis an SVD returns.""" + from koopman_graph.datasets.reaction.oracles import _canonical_row_basis + + system = michaelis_menten() + w = conservation_laws(system) + theta = torch.tensor(0.7, dtype=_F64) + rotation = torch.stack( + [ + torch.stack([torch.cos(theta), -torch.sin(theta)]), + torch.stack([torch.sin(theta), torch.cos(theta)]), + ] + ) + rotated = -(rotation @ w) # another valid orthonormal basis, sign flipped + assert torch.allclose(_canonical_row_basis(rotated), w, atol=1e-12) + # Substrate backbone S + ES + P is the second canonical row. + expected = torch.tensor([0.0, 1.0, 1.0, 1.0], dtype=_F64) / math.sqrt(3.0) + assert torch.allclose(w[1], expected, atol=1e-12) + assert _canonical_row_basis(torch.zeros(0, 4, dtype=_F64)).shape == (0, 4) From bc06eefc6b5f93de6b026c40b5229a5c666e7d89 Mon Sep 17 00:00:00 2001 From: tjkessler Date: Thu, 15 Oct 2026 17:24:33 -0400 Subject: [PATCH 5/5] Snap round-off zeros in canonical conservation laws After f20605b, Linux CI still failed test_cantera_card_oracle_and_provenance and the card builder --check on reaction_cantera_h2o2_v1, again on oracle['conservation_laws']. The basis itself now matches across platforms; what differed were entries that are exactly zero in theory but come out as round-off (1.07e-16, 5.37e-17 on macOS). Their values depend on the LAPACK build, and the Cantera card is compared with atol=1e-18 (appropriate for its ~1e-14 concentrations, not for O(1) basis entries). _canonical_row_basis now sets entries below 1e-12 to exactly zero, so the stored basis is identical across platforms rather than merely close. Cards regenerated: changes are those zeros plus rounding-level shifts (<= 1.8e-15). Test asserts that injected 1e-16 noise canonicalizes to an exact zero. Co-Authored-By: Claude Opus 5.5 --- .../data/reaction_cantera_h2o2_v1.json | 4 ++-- .../data/reaction_michaelis_menten_v1.json | 18 +++++++++--------- .../data/reaction_random_network_v1.json | 4 ++-- src/koopman_graph/datasets/reaction/oracles.py | 11 +++++++++-- tests/datasets/test_reaction_oracles.py | 4 ++++ 5 files changed, 26 insertions(+), 15 deletions(-) diff --git a/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json index 4461e1a..42cdfef 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_cantera_h2o2_v1.json @@ -103,8 +103,8 @@ 0.0 ], [ - 1.0749688212373174e-16, - 5.374844106186587e-17, + 0.0, + 0.0, 0.25819888974716043, 0.5163977794943231, 0.25819888974716104, diff --git a/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json index 295fa86..2f68827 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_michaelis_menten_v1.json @@ -62,10 +62,10 @@ "oracle": { "state": "equilibrium_from_range_midpoints", "reference_state": [ - 0.12499999999999993, - 0.0, - 0.0, - 1.075 + 0.125, + 9.919563547554191e-18, + 8.207559790687939e-18, + 1.0750000000000002 ], "num_conservation_laws": 2, "conservation_laws": [ @@ -76,23 +76,23 @@ -0.25819888974716093 ], [ - -1.9229626863835636e-16, + 0.0, 0.5773502691896258, 0.5773502691896255, 0.5773502691896257 ] ], "dynamic_eigenvalues_real": [ - -0.2499999999999999, - -2.499999999999999 + -0.25, + -2.5 ], "dynamic_eigenvalues_imag": [ 0.0, 0.0 ], "timescales": [ - 4.000000000000002, - 0.40000000000000013 + 4.0, + 0.4 ] } } diff --git a/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json b/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json index 8709ec9..1e7f792 100644 --- a/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json +++ b/src/koopman_graph/datasets/reaction/data/reaction_random_network_v1.json @@ -87,7 +87,7 @@ 2.6685986915205424, 0.0, 1.181401308479456, - 1.0829309696332013e-30, + 5.6833502756612666e-31, 0.0 ], "num_conservation_laws": 2, @@ -101,7 +101,7 @@ -0.5404640055375279 ], [ - 4.60274150786893e-17, + 0.0, 0.30151134457776385, 0.3015113445777636, 0.30151134457776363, diff --git a/src/koopman_graph/datasets/reaction/oracles.py b/src/koopman_graph/datasets/reaction/oracles.py index 7bf9c55..6d2702e 100644 --- a/src/koopman_graph/datasets/reaction/oracles.py +++ b/src/koopman_graph/datasets/reaction/oracles.py @@ -468,6 +468,11 @@ def _stoichiometric_basis(system: ReactionSystem) -> tuple[int, Tensor, Tensor]: _CANONICAL_PIVOT_TOL = 1e-6 +# Entries of the canonical rows that should be exactly zero come out as +# round-off (~1e-16) whose value depends on the LAPACK build; snap them so +# the stored basis is identical across platforms, not just close. +_CANONICAL_ZERO_TOL = 1e-12 + def _canonical_row_basis(w: Tensor) -> Tensor: """Return a basis-independent orthonormal basis of ``rowspace(w)``. @@ -481,7 +486,8 @@ def _canonical_row_basis(w: Tensor) -> Tensor: ------- Orthonormal rows spanning the same subspace, chosen by Gram-Schmidt in coordinate order on the projector ``wᵀw`` (independent of which - orthonormal ``w`` was passed).""" + orthonormal ``w`` was passed), with round-off entries below + ``1e-12`` set to exactly zero.""" num_laws, num_species = int(w.shape[0]), int(w.shape[1]) if num_laws == 0: return w @@ -496,7 +502,8 @@ def _canonical_row_basis(w: Tensor) -> Tensor: basis.append(v / norm) if len(basis) == num_laws: break - return torch.stack(basis) + rows = torch.stack(basis) + return torch.where(rows.abs() < _CANONICAL_ZERO_TOL, torch.zeros_like(rows), rows) def _svd_rank(matrix: Tensor) -> tuple[Tensor, Tensor, int]: diff --git a/tests/datasets/test_reaction_oracles.py b/tests/datasets/test_reaction_oracles.py index bfb3aca..e765428 100644 --- a/tests/datasets/test_reaction_oracles.py +++ b/tests/datasets/test_reaction_oracles.py @@ -331,3 +331,7 @@ def test_conservation_laws_are_canonical_across_bases() -> None: expected = torch.tensor([0.0, 1.0, 1.0, 1.0], dtype=_F64) / math.sqrt(3.0) assert torch.allclose(w[1], expected, atol=1e-12) assert _canonical_row_basis(torch.zeros(0, 4, dtype=_F64)).shape == (0, 4) + # Round-off zeros are exact zeros, so stored bases match across platforms. + assert float(w[1, 0]) == 0.0 + noisy = w + 1e-16 * torch.tensor([[0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0]]) + assert float(_canonical_row_basis(noisy.to(_F64))[1, 0]) == 0.0