diff --git a/CHANGELOG.md b/CHANGELOG.md index ebc4c6f..d01756f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -379,6 +379,57 @@ in `CMakeLists.txt`) is derived from it. the rate and returned the sign-flipped `∂rate/∂E` above. The C++ suite asserts that sign directly, so the sharpest symptom is the first thing to fail. +- **The SBML exporter kept the cancelling tQSSA root after #89 fixed the engine, + so `net2sbml` emitted a model its own validator rejects (follow-up to issue + #89).** `_mm_formula` in `python/bngsim/convert/_sbml_writer.py` was the one + site deliberately left out of #89, on the reading that an interchange artifact + should stay a faithful transcription of the canonical literature formula. That + reading does not survive contact with the measurement: SBML MathML is evaluated + by the *consumer*, in float64, so the exported law performs the same + `0.5*(delta + D)` subtraction the engine had just stopped performing. + + With #89 in place the two sides disagree, and the gate says so. On + `tests/data/mm_tqssa_stiff.net` (ratio ~1e7), before → after: + + | | textbook export | stable export | + |---|---|---| + | RHS vs source model | 2.17e-02 | **0.0** (bit-for-bit) | + | conversion gate | L2 **fail**, L3 **fail**, `ok=False` | L0–L3 pass, `ok=True` | + | libRoadRunner on the artifact | `CV_ERR_FAILURE`, no trajectory | agrees with bngsim to 2.9e-11 | + | rate at t=0 (vs mpmath, 50 dps) | 2.09e-02 | 1.90e-16 | + + So the faithful-transcription reading was shipping a model that bngsim itself + marks not-ok and that RoadRunner cannot integrate. The export now carries the + same branch as `_mm_sfree_c_lines`, spelled as the L3v2 `piecewise` the writer + already targets: + + ``` + sFree = piecewise(0.5*(delta + D), delta >= 0, 2*Km*S/(D - delta)) + ``` + + It omits only that helper's `(D - delta) > 0` degenerate guard, unreachable + once `delta < 0` and backstopped by the enclosing `max(..., 0)`. + + The `piecewise` costs less than it looks. It is not a new construct for this + writer — the target is L3v2 chosen for exactly this MathML, `piecewise` is + already in `_EXPRTK_SUPPORTED_CALLS`, and every BNGL `if()` already exports as + one. It does not move L4 either: that level is *already* `inconclusive` on MM + (pinned by `test_full_gate_l4_is_non_gating`) because the law already contains + `max`, and `_validate.py` lists `Max` and `Piecewise` in the same `_NONSMOOTH` + tuple. And unlike that `max` — a genuine kink — this branch is *removable*: + both arms agree to the last ulp at `delta == 0` (measured across + `delta = ±1e-3 … ±1e-12`, worst 1.11e-16), so no consumer's integrator gains a + discontinuity to resolve. The cost that is real is size: the emitted document + grows 6519 → 12121 bytes on the stiff fixture, because `sFree` appears twice in + the law and SBML has no `let` binding. + + Where the corpus lives nothing moves at all: at ratio ~1 and on the `delta >= 0` + branch the emitted law is **bit-identical** to the textbook one, so both corpus + MM models (`test_MM`, `mCaMKII_Ca_Spike`) export unchanged. Pinned by + `test_mm_tqssa_stiff_rhs_exact` and `test_mm_tqssa_stiff_full_gate_passes`; the + reasoning is recorded in the `_mm_formula` docstring so the scope question does + not get re-opened from first principles. + - **The conservation-law reduction chose dependent species it could not solve for, so the reduced system was singular on models that are perfectly well posed (follow-up to issue #63).** `detect_conservation_laws` found every law — the diff --git a/python/bngsim/convert/_sbml_writer.py b/python/bngsim/convert/_sbml_writer.py index 22be0c1..2e4e5a7 100644 --- a/python/bngsim/convert/_sbml_writer.py +++ b/python/bngsim/convert/_sbml_writer.py @@ -1158,13 +1158,53 @@ def _mm_formula( """Closed-form Michaelis–Menten (tQSSA) kinetic law, matching ``model.cpp``. With E = first reactant, S = second reactant, and (kcat, Km) the two rate - parameters:: + parameters, ``delta = S - Km - E`` and ``D = sqrt(delta^2 + 4*Km*S)``:: - sFree = 0.5 * ((S - Km - E) + sqrt((S - Km - E)^2 + 4*Km*S)) + sFree = piecewise(0.5 * (delta + D), delta >= 0, 2*Km*S / (D - delta)) rate = kcat * stat_factor * max(sFree, 0) * E / (Km + max(sFree, 0)) (The engine floors a negative ``sFree`` to 0; ``max`` reproduces that. SBML expresses the whole law explicitly — net2sbml is the faithful half.) + + **Why the export carries the branch and not the textbook one-liner (GH #89).** + ``sFree`` is the positive root of ``x^2 - delta*x - Km*S = 0``. Spelled + ``0.5*(delta + D)`` it subtracts two nearly-equal positive numbers once + ``delta < 0``, losing ~2 significant digits per decade of + ``|delta|/sqrt(4*Km*S)`` — no correct digit left by 1e8. This is not a + bngsim-internal concern that stops at the engine boundary: the *consumer* + evaluating the exported MathML in float64 performs the same subtraction. So + "stay a faithful transcription of the canonical literature formula" is not a + neutral choice — measured on ``tests/data/mm_tqssa_stiff.net`` (ratio ~1e7), + once the engine took the stable root in #89 the textbook export + + * disagreed with its own source model by 2.17e-02 relative RHS (the branch + above reproduces the engine **bit-for-bit**, delta 0.0), + * therefore **failed its own conversion gate** at L2 (reverse round-trip) and + L3 (numerical), so ``net_to_sbml`` emitted an artifact ``report.ok`` marked + False, and + * made libRoadRunner's CVODE abort with ``CV_ERR_FAILURE`` instead of + integrating — where the branch above agrees with bngsim to 2.9e-11. + + Both regressions are pinned by ``test_mm_tqssa_stiff_*`` in + ``python/tests/test_net2sbml.py``. + + The branch mirrors ``_mm_sfree_c_lines`` in ``python/bngsim/_codegen.py`` + operation-for-operation, so engine and export agree bit-for-bit; it omits only + that helper's ``(D - delta) > 0`` degenerate guard, which is unreachable here + (``delta < 0`` and ``D >= 0`` force ``D - delta > 0``) and which the enclosing + ``max(..., 0)`` would backstop anyway. + + Cost of the ``piecewise``, weighed before adopting it: it is *not* a new + construct for this writer — the target is L3v2 chosen for exactly this MathML + (see ``_SBML_VERSION``), ``piecewise`` is already in + ``_EXPRTK_SUPPORTED_CALLS``, and every BNGL ``if()`` already exports as one. + It costs nothing at L4 either: that level is *already* ``inconclusive`` for MM + (pinned by ``test_full_gate_l4_is_non_gating``) because the law already + contains ``max``, and ``_validate.py`` lists ``Max`` and ``Piecewise`` in the + same ``_NONSMOOTH`` tuple. Note the asymmetry between the two: ``max(sFree, 0)`` + is a genuine kink, whereas this branch is *removable* — both arms agree to the + last ulp at ``delta == 0`` and the function is smooth across it — so it adds no + discontinuity a consumer's integrator has to resolve. """ reactants = list(rxn["reactants"]) rate_idx = list(rxn.get("rate_param_indices", [])) @@ -1182,7 +1222,11 @@ def _mm_formula( sf = float(rxn["stat_factor"]) delta = f"(({S}) - ({Km}) - ({E}))" - s_free = f"(0.5 * ({delta} + sqrt({delta}^2 + 4 * ({Km}) * ({S}))))" + d_root = f"sqrt({delta}^2 + 4 * ({Km}) * ({S}))" + s_free = ( + f"piecewise(0.5 * ({delta} + {d_root}), {delta} >= 0, " + f"2 * ({Km}) * ({S}) / ({d_root} - {delta}))" + ) s_free = f"max({s_free}, 0)" pre = f"({_fmt(sf)} * " if sf != 1.0 else "(" return f"{pre}({kcat}) * ({s_free}) * ({E}) / (({Km}) + ({s_free})))" diff --git a/python/tests/test_net2sbml.py b/python/tests/test_net2sbml.py index d0c4a47..ac30d5d 100644 --- a/python/tests/test_net2sbml.py +++ b/python/tests/test_net2sbml.py @@ -126,6 +126,47 @@ def test_mm_tqssa_rhs_exact(data_dir: Path, tmp_path: Path) -> None: assert _rhs_max_delta(src_model, sbml_model) <= 1e-9 +def test_mm_tqssa_stiff_rhs_exact(data_dir: Path, tmp_path: Path) -> None: + """The exported law must carry GH #89's *stable* free-substrate root. + + ``mm_tqssa_stiff.net`` sits at ``|delta|/sqrt(4*Km*S) ~ 1e7``, where spelling + ``sFree`` as the textbook ``0.5*(delta + D)`` subtracts two nearly-equal + positive numbers. The engine takes the conjugate root there (#89); an export + that keeps the textbook form disagrees with it by ~2e-2 relative — this + assertion measured 2.17e-02 against the textbook form, and exactly 0.0 + (bit-for-bit) against the piecewise the writer now emits. + """ + src = data_dir / "mm_tqssa_stiff.net" + if not src.is_file(): + pytest.skip("mm_tqssa_stiff.net not present") + with warnings.catch_warnings(): + warnings.simplefilter("ignore", bngsim.ConversionWarning) + src_model = bngsim.Model.from_net(src) + write_sbml(src_model, tmp_path / "mm_stiff.xml", strict=True) + sbml_model = bngsim.Model.from_sbml(tmp_path / "mm_stiff.xml") + assert _rhs_max_delta(src_model, sbml_model) <= 1e-9 + + +def test_mm_tqssa_stiff_full_gate_passes(data_dir: Path, tmp_path: Path) -> None: + """The stiff fixture must clear the gating levels, not just L4's punt. + + This is the end-to-end tripwire for the same regression: with the textbook + root the emitted artifact diverges from its source badly enough that **L2 + (reverse round-trip) and L3 (numerical) both fail** and ``report.ok`` is + False — i.e. net2sbml would emit a model its own validator rejects. + """ + src = data_dir / "mm_tqssa_stiff.net" + if not src.is_file(): + pytest.skip("mm_tqssa_stiff.net not present") + with warnings.catch_warnings(): + warnings.simplefilter("ignore", bngsim.ConversionWarning) + report = net_to_sbml(src, tmp_path / "mm_stiff.xml", validate="full", strict=True) + for lv in ("L0", "L1", "L2", "L3"): + level = report.validation.level(lv) + assert level is not None and level.status == "pass", report.summary() + assert report.ok, report.summary() + + # ─── ExprTk→MathML translator unit checks ──────────────────────────────────