fix(convert): carry #89's stable tQSSA root into the SBML export

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

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})))"

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 ──────────────────────────────────