From 1373893c74a3340d000cbdcbf3e29387eba9aded Mon Sep 17 00:00:00 2001 From: Bill Hlavacek Date: Mon, 27 Jul 2026 20:41:51 -0600 Subject: [PATCH] fix(mm): stable tQSSA free-substrate root and subtraction-free partials (#89) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Michaelis-Menten free substrate was computed as delta = S - Km - E D = sqrt(delta*delta + 4*Km*S) sFree = 0.5*(delta + D) When delta < 0 that last line subtracts two nearly-equal positive numbers and loses about two significant digits per decade of |delta|/sqrt(4*Km*S): 8 correct digits at ratio 1e4, none at 1e8. This is in the RHS, so a model in deep enzyme excess had a wrong trajectory before it had a wrong gradient, and every derivative inherited it. sFree is the positive root of x^2 - delta*x - Km*S = 0, so the delta < 0 branch now uses the conjugate form 2*Km*S/(D - delta), which multiplies out to the same value with no subtraction. delta >= 0 keeps the textbook expression bit-for-bit. Fixing the root was necessary and not sufficient. The derivatives were the chain rule through sFree, dsF/dE = (-1 - delta/D)/2 dsF/dS = ( 1 + (delta + 2*Km)/D)/2 dsF/dKm = (-1 + (2*S - delta)/D)/2 and those cancel in exactly the same regime, so d(rate)/dE stayed at 1e+10 relative error with a correct sFree in hand. On the stiff fixture the cancellation went past magnitude and flipped the sign: -3.83e-05 where the true value is +4.96e-15. A negative d(rate)/dE says adding enzyme slows the reaction, and that entry drives CVODE's Newton step, so it steps confidently the wrong way rather than merely slowly. Differentiating the symmetric form of the rate collapses each partial to one subtraction-free quotient. The tQSSA complex is c = (A - D)/2 with A = E + S + Km, and that same D is sqrt(A^2 - 4*E*S), so dc/dE = (S - c)/D with S - c = sFree: d(rate)/dE = kcat*stat*sFree/D d(rate)/dS = kcat*stat*E*Km/((Km + sFree)*D) d(rate)/dKm = -rate/D Each is identical to sympy.diff of the shipped rate. Measured against mpmath at 60 digits over the four regime sweeps from the issue, worst relative error in float64, before -> after: sweep sFree rate d/dE uniform O(1) 4.5e-14 -> 4.1e-16 4.3e-14 -> 4.9e-16 5.6e-13 -> 6.1e-16 log-spread 6.2e-04 -> 3.5e-15 6.2e-04 -> 5.8e-16 3.0e+03 -> 1.8e-15 Km >> S,E 2.6e-09 -> 5.4e-16 2.6e-09 -> 6.5e-16 2.6e-09 -> 7.6e-16 Km << S,E 2.1e-02 -> 2.4e-13 2.1e-02 -> 5.1e-16 2.3e+09 -> 4.3e-14 sweep d/dS d/dKm uniform O(1) 2.6e-15 -> 7.3e-16 1.4e-12 -> 7.6e-16 log-spread 2.5e-10 -> 6.1e-15 2.6e+03 -> 2.4e-15 Km >> S,E 7.6e-16 -> 4.7e-16 9.4e-09 -> 7.9e-16 Km << S,E 5.3e-05 -> 4.7e-13 8.7e+09 -> 2.3e-13 The three lines were open-coded at every site rather than shared, which is why the issue asks for a grep before starting. They are now emitted from one helper per language: mm_tqssa() in include/bngsim/mm_jacobian.hpp, which compute_rxn_rate and mm_tqssa_derivatives both call, and _mm_sfree_c_lines() in python/bngsim/_codegen.py, which every emitter calls. Covered: the interpreted RHS and the SSA propensity emitter (src/model.cpp), the closed-form analytical Jacobian shared by the dense and sparse CVODE paths, both compiled-RHS emitters, the analytical Jacobian and its fused J*v, #55's analytic df/dp, and the JAX RHS (where the denominator is masked before the divide so the unselected branch cannot inject a NaN into a tangent). _CODEGEN_VERSION -> 26: MM models' emitted C changes, and a cached v25 .so would keep serving pre-fix numbers. The corpus cannot gate this. Both MM models (test_MM, mCaMKII_Ca_Spike) do reach delta < 0, so they take the new branch, but at |delta|/sqrt(4*Km*S) <= 4.95 and <= 1.32, where the two forms differ by 1-3 ulp — largest sFree shift measured over their trajectories is 1.9e-16 and 6.0e-16. generate_rhs_c is byte-identical on 120 sampled non-MM corpus models and differs only on those two. So the regression fixture is a deliberately stiff one, tests/data/mm_tqssa_stiff.net at ratio ~1e7, shared by the C++ and Python suites. Both halves are independently gated. Reverting only the root fails test_mm_tqssa_stiff_root at rel=0.0199; reverting only the derivatives fails it at rel=1.36e10, with the sign assertion on its own line so the sharpest symptom is the first thing to fail. TestNumericalStability pins each rejected form against mpmath — the textbook root, the chain rule through sFree, and #55's already-rejected "simplified" d/dKm — so a future tidy-up fails there instead of silently in a trajectory or a gradient. TestTheEmittedCAtAStiffRatio compiles the emitted RHS and df/dp and checks them against mpmath rather than against a finite difference of themselves, which is blind to a cancellation both sides share. Not touched: _mm_formula in python/bngsim/convert/_sbml_writer.py still exports the cancelling form. Fixing it means a piecewise() in the exported MathML, which changes the shape of every exported MM kinetic law and touches both the symbolic round-trip check and downstream consumers; that is its own change, not a rider on this one. --- CHANGELOG.md | 66 ++++++ include/bngsim/mm_jacobian.hpp | 79 +++++-- python/bngsim/_codegen.py | 212 +++++++++++-------- python/bngsim/_jax_rhs.py | 10 +- python/tests/test_codegen_mm_sens.py | 302 +++++++++++++++++++++++---- src/model.cpp | 15 +- tests/data/mm_tqssa_stiff.net | 36 ++++ tests/test_bngsim.cpp | 67 ++++++ 8 files changed, 633 insertions(+), 154 deletions(-) create mode 100644 tests/data/mm_tqssa_stiff.net diff --git a/CHANGELOG.md b/CHANGELOG.md index 7f9a4ca..ebc4c6f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -313,6 +313,72 @@ in `CMakeLists.txt`) is derived from it. fields. ### Fixed +- **The Michaelis–Menten (tQSSA) free substrate lost about two significant digits + per decade of `|δ|/√(4·Km·S)`, in the rate itself and not only in a derivative + (issue #89).** `sFree` was computed as + + ``` + delta = S - Km - E; D = sqrt(delta*delta + 4*Km*S); sFree = 0.5*(delta + D); + ``` + + When `delta < 0` — a large enzyme excess, or a small `Km·S` — that last line + subtracts two nearly-equal positive numbers. At `|δ|/√(4·Km·S) = 1e4` the result + has 8 correct digits, at 1e8 it has **none**. `sFree` is the positive root of + `x² − δ·x − Km·S = 0`, so the `delta < 0` branch now uses the conjugate form + `2·Km·S/(D − delta)`, which multiplies out to the same value with no subtraction + and is exact to 0–1 ulp at every ratio measured. `delta ≥ 0` keeps the textbook + expression bit-for-bit. + + Fixing the root was necessary but not sufficient. The derivatives were written + as the chain rule through `sFree` — `∂sFree/∂E = ½(−1 − δ/D)`, + `∂sFree/∂S = ½(1 + (δ + 2·Km)/D)`, `∂sFree/∂Km = ½(−1 + (2S − δ)/D)` — and those + cancel in exactly the same regime, so `∂rate/∂E` stayed at a relative error of + **1e+10** with a correct `sFree` in hand. On the regression fixture the + cancellation went past magnitude and **flipped the sign**: `∂rate/∂E` came out + as `-3.83e-05` where the true value is `+4.96e-15`. A negative `∂rate/∂E` says + adding enzyme slows the reaction, and this entry feeds CVODE's Newton + iteration — a sign error there steps confidently the wrong way, which is worse + than converging slowly. Differentiating the symmetric form of + the rate instead (the tQSSA complex is `c = ½(A − D)` with `A = E + S + Km`, and + that same `D` is `√(A² − 4·E·S)`) collapses each partial to one subtraction-free + quotient: + + ``` + ∂rate/∂E = kcat·stat·sFree/D + ∂rate/∂S = kcat·stat·E·Km/((Km + sFree)·D) + ∂rate/∂Km = −rate/D + ``` + + Each is identical to `sympy.diff` of the shipped rate. Measured against mpmath + at 60 digits over the four regime sweeps from the issue — worst relative error + in float64, before → after: + + | sweep | `sFree` | `rate` | `∂/∂E` | `∂/∂S` | `∂/∂Km` | + |---|---|---|---|---|---| + | uniform, O(1) | 4.5e-14 → 4.1e-16 | 4.3e-14 → 4.9e-16 | 5.6e-13 → 6.1e-16 | 2.6e-15 → 7.3e-16 | 1.4e-12 → 7.6e-16 | + | log-spread 1e-4..1e3 | 6.2e-04 → 3.5e-15 | 6.2e-04 → 5.8e-16 | 3.0e+03 → 1.8e-15 | 2.5e-10 → 6.1e-15 | 2.6e+03 → 2.4e-15 | + | Km ≫ S,E | 2.6e-09 → 5.4e-16 | 2.6e-09 → 6.5e-16 | 2.6e-09 → 7.6e-16 | 7.6e-16 → 4.7e-16 | 9.4e-09 → 7.9e-16 | + | Km ≪ S,E (deep saturation) | 2.1e-02 → 2.4e-13 | 2.1e-02 → 5.1e-16 | 2.3e+09 → 4.3e-14 | 5.3e-05 → 4.7e-13 | 8.7e+09 → 2.3e-13 | + + This is a **trajectory** change, not only a gradient one, and it lands in every + place the three lines were open-coded: `compute_rxn_rate`'s MM branch and the + SSA propensity emitter (`src/model.cpp`), the closed-form analytical Jacobian + shared by the dense and sparse CVODE paths (`include/bngsim/mm_jacobian.hpp`), + the two compiled-RHS emitters, the analytical Jacobian / fused `J·v`, and the + analytic `∂f/∂p` from #55 (`python/bngsim/_codegen.py`), plus the JAX RHS + (`python/bngsim/_jax_rhs.py`). The expression was duplicated at each site; it is + now emitted from one helper per language. `_CODEGEN_VERSION` → 26, since a + cached v25 `.so` would keep serving pre-fix numbers. + + Both corpus MM models (`test_MM`, `mCaMKII_Ca_Spike`) run at + `|δ|/√(4·Km·S) ≤ 5`, where the two forms differ by 1–3 ulp — measured largest + `sFree` shift 1.9e-16 and 6.0e-16 — so the corpus shows no trajectory change and + the emitted C is byte-identical on the 120 non-MM models sampled. The + regression fixture is therefore a deliberately stiff one + (`tests/data/mm_tqssa_stiff.net`, ratio ~1e7), where the old code was 2% out on + 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 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/include/bngsim/mm_jacobian.hpp b/include/bngsim/mm_jacobian.hpp index 13beb4a..d8f0b58 100644 --- a/include/bngsim/mm_jacobian.hpp +++ b/include/bngsim/mm_jacobian.hpp @@ -1,21 +1,40 @@ -// bngsim/include/bngsim/mm_jacobian.hpp — Michaelis–Menten (tQSSA) closed-form -// Jacobian derivative (GH #76 task 3) +// bngsim/include/bngsim/mm_jacobian.hpp — Michaelis–Menten (tQSSA) closed form +// (GH #76 task 3, made numerically stable in GH #89) // -// Single source of truth for ∂rate/∂E and ∂rate/∂S of the tQSSA MM rate law, -// shared by the dense path (NetworkModel::fill_dense_analytical_jacobian) and the -// sparse CVODE callback (cvode_analytical_jac). Derived by hand to match -// compute_rxn_rate's MM formula exactly: +// Single source of truth for the tQSSA free substrate and for ∂rate/∂E and +// ∂rate/∂S: the free substrate is shared with compute_rxn_rate's MM branch, the +// derivatives with the dense path (NetworkModel::fill_dense_analytical_jacobian) +// and the sparse CVODE callback (cvode_analytical_jac). // // delta = S - Km - E // D = sqrt(delta² + 4·Km·S) // sFree = ½·(delta + D) (clamped to 0 in the RHS) // rate = kcat·stat·sFree·E/(Km + sFree) // -// Chain rule through sFree (∂delta/∂E = -1, ∂delta/∂S = +1): -// ∂sFree/∂E = ½·(-1 - delta/D) -// ∂sFree/∂S = ½·(+1 + (delta + 2·Km)/D) -// ∂rate/∂E = kcat·stat·[ sFree/(Km+sFree) + E·Km/(Km+sFree)²·∂sFree/∂E ] -// ∂rate/∂S = kcat·stat·E·Km/(Km+sFree)²·∂sFree/∂S +// **Stability (GH #89).** Two cancellations used to live in these lines, both +// fatal once delta < 0 with |delta| ≫ √(4·Km·S) — the deep enzyme excess regime. +// +// 1. ½·(delta + D) subtracts two nearly-equal positive numbers, losing about two +// significant digits per decade of that ratio (no correct digit left at 1e8). +// sFree is the positive root of x² − delta·x − Km·S = 0, so for delta < 0 the +// conjugate form 2·Km·S/(D − delta) multiplies out to the same value with no +// subtraction at all. For delta ≥ 0 the textbook form is already +// cancellation-free and is kept bit-for-bit. +// +// 2. The chain rule through sFree used to be written ∂sFree/∂E = ½·(−1 − delta/D) +// and ∂sFree/∂S = ½·(1 + (delta + 2·Km)/D), which cancel in exactly the same +// regime — so fixing sFree alone still left ∂rate/∂E with a relative error of +// 1e+10 on a deep-saturation sweep. Differentiating the *symmetric* form of the +// rate collapses each partial to a single subtraction-free quotient (the tQSSA +// complex is c = ½·(A − D) with A = E + S + Km, and that same D is √(A² − 4·E·S), +// so ∂c/∂E = (S − c)/D = sFree/D and ∂c/∂Km = −c/D): +// +// ∂rate/∂E = kcat·stat·sFree/D +// ∂rate/∂S = kcat·stat·E·Km/((Km + sFree)·D) +// ∂rate/∂Km = −rate/D (emitted by python/bngsim/_codegen.py) +// +// Verified identical to sympy.diff of the shipped rate, and measured against +// mpmath at 60 digits: machine precision in every column of all four sweeps. // // Where the RHS clamps sFree to 0 (delta + D ≤ 0), the rate is identically 0 in a // neighbourhood, so both derivatives are 0 — matching the clamped flat region. @@ -26,24 +45,42 @@ namespace bngsim { +struct MMTqssa { + double delta; // S - Km - E + double D; // sqrt(delta² + 4·Km·S) + double sFree; // positive root of x² - delta·x - Km·S = 0 (unclamped) +}; + +// Free substrate as the *stable* quadratic root — see note 1 in the header +// comment. Callers that feed a rate still apply the RHS's own `sFree < 0` clamp; +// this returns the root as computed so the derivative guard can see it. +inline MMTqssa mm_tqssa(double Km, double E, double S) { + MMTqssa q; + q.delta = S - Km - E; + q.D = std::sqrt(q.delta * q.delta + 4.0 * Km * S); + if (q.delta >= 0.0) { + q.sFree = 0.5 * (q.delta + q.D); + } else { + // D - delta > 0 whenever delta < 0 and D ≥ 0, so this cannot divide by + // zero on real inputs; the guard covers a NaN/negative-S degeneracy. + double denom = q.D - q.delta; + q.sFree = denom > 0.0 ? 2.0 * Km * S / denom : 0.0; + } + return q; +} + inline void mm_tqssa_derivatives(double kcat, double Km, double E, double S, double stat, double &dE, double &dS) { - double delta = S - Km - E; - double D = std::sqrt(delta * delta + 4.0 * Km * S); - double sFree = 0.5 * (delta + D); - if (sFree <= 0.0 || D <= 0.0) { + MMTqssa q = mm_tqssa(Km, E, S); + if (q.sFree <= 0.0 || q.D <= 0.0) { // Clamped (rate ≡ 0 locally) or degenerate — flat, zero derivative. dE = 0.0; dS = 0.0; return; } - double dsF_dE = 0.5 * (-1.0 - delta / D); - double dsF_dS = 0.5 * (1.0 + (delta + 2.0 * Km) / D); double C = kcat * stat; - double KpsF = Km + sFree; - double common = C * E * Km / (KpsF * KpsF); // C·E·Km/(Km+sFree)² - dE = C * sFree / KpsF + common * dsF_dE; - dS = common * dsF_dS; + dE = C * q.sFree / q.D; + dS = C * E * Km / ((Km + q.sFree) * q.D); } } // namespace bngsim diff --git a/python/bngsim/_codegen.py b/python/bngsim/_codegen.py index 8c78b56..ada0951 100644 --- a/python/bngsim/_codegen.py +++ b/python/bngsim/_codegen.py @@ -204,7 +204,13 @@ def _codegen_jit_backend() -> str: # builder the analytical Jacobian uses), so an MM model emits a # bngsim_codegen_sens_rhs where v24 declined to CVODES' difference quotient. # Invalidate v24. -_CODEGEN_VERSION = "25" +# v26: lanl/bngsim #89 — every emitted Michaelis–Menten site changes. The free +# substrate is now the stable quadratic root (the conjugate form for δ < 0), and +# the Jacobian / ∂f/∂p partials are the subtraction-free quotients rather than +# the chain rule through ∂sFree/∂E,S,Km. This moves MM *trajectories*, not only +# gradients, so a cached v25 .so would keep serving pre-fix numbers. +# Invalidate v25. +_CODEGEN_VERSION = "26" # Modules whose *source* determines the emitted C. ``_codegen`` holds the @@ -1805,6 +1811,7 @@ def generate_rhs_c(net_path: str) -> str: grp: list[str] = [] g = grp.append kind = _classify_rate_law(rate_law, func_names) + rate_expr: str | None = None if kind[0] == "elementary": _, pname, sf = kind rate_expr = _rate_elementary(pname, sf, reactants, param_idx, func_idx) @@ -1812,11 +1819,25 @@ def generate_rhs_c(net_path: str) -> str: _, fname, sf = kind rate_expr = _rate_functional(fname, sf, reactants, func_idx, use_array=chunk) elif kind[0] == "mm": + # A braced block, not an expression: the stable free-substrate root + # is a branch on delta's sign (GH #89), and inlining it would repeat + # the sqrt four times. _, kcat, km, sf = kind - rate_expr = _rate_mm(kcat, km, sf, reactants, param_idx) + if len(reactants) >= 2: + for ln in _mm_rate_lines( + f"p[{param_idx[kcat]}]" if kcat in param_idx else "0.0", + f"p[{param_idx[km]}]" if km in param_idx else "0.0", + sf, + reactants[0] - 1, + reactants[1] - 1, + ): + g(ln) + else: + rate_expr = "0.0" else: rate_expr = "0.0" - g(f" rate = {rate_expr};") + if rate_expr is not None: + g(f" rate = {rate_expr};") # Subtract from reactants (index 0 = null reactant, skip) for ri in reactants: if ri > 0: @@ -2201,35 +2222,49 @@ def _rate_functional( return " * ".join(parts) -def _rate_mm( - kcat: str, - km: str, - sf: float, - reactants: list, - param_idx: dict, -) -> str: - """Generate C expression for Michaelis-Menten tQSSA rate.""" - kcat_c = f"p[{param_idx[kcat]}]" if kcat in param_idx else "0.0" - km_c = f"p[{param_idx[km]}]" if km in param_idx else "0.0" +def _mm_sfree_c_lines(km_c: str, e_idx: int, s_idx: int, indent: str) -> list[str]: + """C lines declaring ``Km``/``E``/``S``/``delta``/``Dmm``/``sFree`` for one + tQSSA Michaelis–Menten reaction — the single source of the free-substrate + root for every emitter on this path (GH #89). - if len(reactants) >= 2: - e_idx = reactants[0] - 1 - s_idx = reactants[1] - 1 - else: - return "0.0" + ``sFree`` is the positive root of ``x² − delta·x − Km·S = 0``. Writing it as + ``½(delta + D)`` subtracts two nearly-equal positive numbers once + ``delta < 0``, losing about two significant digits per decade of + ``|delta|/√(4·Km·S)`` — no correct digit left at 1e8, in the *rate* and so in + everything derived from it. The conjugate form ``2·Km·S/(D − delta)`` + multiplies out to the same value with no subtraction. ``delta ≥ 0`` keeps the + textbook form, which is already cancellation-free, bit-for-bit. + """ + return [ + f"{indent}double Km = {km_c}, E = y[{e_idx}], S = y[{s_idx}];", + f"{indent}double delta = S - Km - E;", + f"{indent}double Dmm = sqrt(delta*delta + 4.0*Km*S);", + f"{indent}double sFree = (delta >= 0.0) ? 0.5*(delta + Dmm)", + f"{indent} : ((Dmm - delta) > 0.0 ? 2.0*Km*S/(Dmm - delta) : 0.0);", + ] - sf_c = f"{sf} * " if sf != 1.0 else "" - return ( - f"({sf_c}{kcat_c} * (0.5 * ((y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"sqrt((y[{s_idx}] - {km_c} - y[{e_idx}]) * " - f"(y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"4.0 * {km_c} * y[{s_idx}]))) * y[{e_idx}] / " - f"({km_c} + 0.5 * ((y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"sqrt((y[{s_idx}] - {km_c} - y[{e_idx}]) * " - f"(y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"4.0 * {km_c} * y[{s_idx}]))))" - ) +def _mm_rate_lines( + kcat_c: str, + km_c: str, + sf: float, + e_idx: int, + s_idx: int, + indent: str = " ", +) -> list[str]: + """C lines assigning the Michaelis–Menten tQSSA ``rate`` (a braced block, so + the locals cannot collide with the caller's). + + ``e_idx``/``s_idx`` are 0-based species indices — the enzyme is the *first* + reactant and the substrate the second, matching ``src/model.cpp``'s MM branch. + """ + sf_c = f"{sf} * " if sf != 1.0 else "" + return [ + f"{indent}{{", + *_mm_sfree_c_lines(km_c, e_idx, s_idx, indent + " "), + f"{indent} rate = {sf_c}{kcat_c} * sFree * E / (Km + sFree);", + f"{indent}}}", + ] # ─── Sensitivity RHS code generation ──────────────────────────────────── @@ -4395,23 +4430,17 @@ def generate_rhs_from_model(model) -> str: g(f" rate = {' * '.join(parts)};") elif rtype == "mm": - # tQSSA Michaelis-Menten + # tQSSA Michaelis-Menten — the free substrate through the stable + # quadratic root (GH #89), shared with generate_rhs_c. if len(rate_params) >= 2 and len(reactants) >= 2: - kcat_c = f"p[{rate_params[0]}]" - km_c = f"p[{rate_params[1]}]" - e_idx = reactants[0] - s_idx = reactants[1] - sf_c = f"{sf} * " if sf != 1.0 else "" - g( - f" rate = ({sf_c}{kcat_c} * (0.5 * ((y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"sqrt((y[{s_idx}] - {km_c} - y[{e_idx}]) * " - f"(y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"4.0 * {km_c} * y[{s_idx}]))) * y[{e_idx}] / " - f"({km_c} + 0.5 * ((y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"sqrt((y[{s_idx}] - {km_c} - y[{e_idx}]) * " - f"(y[{s_idx}] - {km_c} - y[{e_idx}]) + " - f"4.0 * {km_c} * y[{s_idx}]))));" - ) + for ln in _mm_rate_lines( + f"p[{rate_params[0]}]", + f"p[{rate_params[1]}]", + sf, + reactants[0], + reactants[1], + ): + g(ln) else: g(" rate = 0.0; /* malformed MM */") else: @@ -4800,11 +4829,18 @@ def _mm_v_lines( Written as a braced block with locals rather than one inline expression, and grouped exactly as ``_mm_jacobian_groups`` groups ∂rate/∂E and ∂rate/∂S, - because the grouping is load-bearing: the algebraically-identical "simplified" - form of ∂/∂Km (with the ½ factors cancelled) loses every significant digit on - log-spread parameters, while this one tracks ``sFree`` itself. No derivative - can beat the value it differentiates — see the note in - ``_mm_dfdp_terms`` about ``sFree``'s own cancellation. + because the grouping is load-bearing: + + ∂rate/∂kcat = stat·E·sFree/(Km + sFree) ( = rate/kcat) + ∂rate/∂Km = −kcat·stat·E·sFree/((Km + sFree)·D) ( = −rate/D) + + Every algebraically-identical alternative tried so far is worse in float64. + The "simplified" ∂/∂Km with the ½ factors cancelled loses every significant + digit on log-spread parameters; the chain rule through + ``∂sFree/∂Km = ½(−1 + (2S − δ)/D)`` cancels catastrophically in deep + saturation (1e+10 relative error there, vs machine precision for the form + above — GH #89). Both are pinned by + ``test_the_emitted_grouping_is_the_stable_one``. ``chain_c`` is the #15/#41 factor ∂(kcat or Km)/∂primary when the rate constant is a derived parameter; ``None`` for the direct column. @@ -4813,14 +4849,9 @@ def _mm_v_lines( stat_c = "" if stat == 1.0 else f"{_jac_c_float(stat)} * " out = [" {"] if wrt == "Km": - out.append(f" double kcat = {kcat_c}, Km = {km_c};") - else: - out.append(f" double Km = {km_c};") + out.append(f" double kcat = {kcat_c};") + out += _mm_sfree_c_lines(km_c, e_idx, s_idx, " ") out += [ - f" double E = y[{e_idx}], S = y[{s_idx}];", - " double delta = S - Km - E;", - " double Dmm = sqrt(delta*delta + 4.0*Km*S);", - " double sFree = 0.5*(delta + Dmm);", " v = 0.0;", # At sFree == 0 the rate is 0 for every kcat/Km (S = 0 forces sFree = 0 # whatever Km is), so both partials are genuinely 0 — the same guard the @@ -4830,11 +4861,7 @@ def _mm_v_lines( if wrt == "kcat": out.append(f" v = {stat_c}E*sFree/(Km + sFree){tail};") else: - out += [ - " double dsF_dKm = 0.5*(-1.0 + (2.0*S - delta)/Dmm);", - " double KpsF = Km + sFree;", - f" v = kcat * {stat_c}E*(dsF_dKm*Km - sFree)/(KpsF*KpsF){tail};", - ] + out.append(f" v = -kcat * {stat_c}E*sFree/((Km + sFree)*Dmm){tail};") out += [" }", " }"] return out @@ -4847,21 +4874,22 @@ def _mm_dfdp_terms(data, plan_mm, param_idx_by_name, primary_param_names, derive else on this path. The tQSSA rate is closed form, so unlike the Functional path (GH #66) there is - no sympy here at all: ∂/∂kcat = rate/kcat and ∂/∂Km follows from - ``∂sFree/∂Km = ½(−1 + (2S − δ)/D)``. Both were checked against ``sympy.diff`` - symbolically and over random parameter points before being written down. + no sympy here at all: ∂/∂kcat = rate/kcat and ∂/∂Km = −rate/D. Both were + checked against ``sympy.diff`` symbolically and over random parameter points + before being written down. Cross-checked against ``plan_mm`` entry by entry rather than trusted: the plan is what ``_mm_jacobian_groups`` builds ``J·v`` from, so if the two disagree about which species is the enzyme, this declines instead of emitting a ∂f/∂p that belongs to a different reaction than the ``J·v`` beside it. - **Accuracy floor.** ``sFree = ½(δ + D)`` cancels catastrophically once - ``|δ| ≫ √(4·Km·S)`` — roughly two digits per decade — and that is shipped - behaviour in the RHS and the analytical Jacobian, not something introduced - here. These partials are exactly as accurate as the rate they differentiate - and as the Jacobian's own ∂rate/∂E; a model in that regime has a wrong - trajectory before it has a wrong gradient. + **Accuracy.** ``sFree`` used to be written ``½(δ + D)``, which cancels + catastrophically once ``|δ| ≫ √(4·Km·S)`` — roughly two digits per decade, in + the RHS and the analytical Jacobian as much as here, so a model in that regime + had a wrong trajectory before it had a wrong gradient. GH #89 replaced the + root and these partials' grouping both; see :func:`_mm_sfree_c_lines` and + :func:`_mm_v_lines`. These partials remain exactly as accurate as the rate + they differentiate and as the Jacobian's own ∂rate/∂E. """ mm_rxns = [(i, r) for i, r in enumerate(data["reactions"]) if r["type"] == "mm"] if not mm_rxns: @@ -4957,10 +4985,24 @@ def _mm_jacobian_groups(plan_mm, add) -> list[list[str]] | None: """Reconstruct every Michaelis–Menten reaction's Jacobian contribution as a balanced ``{ … }`` C line-group, scattered through the caller's ``add``. - The tQSSA rate is ``kcat·stat·E·sFree/(Km + sFree)`` with - ``sFree = ½(δ + D)``, ``δ = S − Km − E``, ``D = √(δ² + 4·Km·S)`` — closed - form, so unlike the Functional path there is no symbolic differentiation - here, just the two partials written out. + The tQSSA rate is ``kcat·stat·E·sFree/(Km + sFree)`` with ``sFree`` the + positive root of ``x² − δ·x − Km·S = 0``, ``δ = S − Km − E``, + ``D = √(δ² + 4·Km·S)`` — closed form, so unlike the Functional path there is + no symbolic differentiation here, just the two partials written out: + + ∂rate/∂E = kcat·stat·sFree/D + ∂rate/∂S = kcat·stat·E·Km/((Km + sFree)·D) + + Both are single subtraction-free quotients, and that is load-bearing (GH #89). + The obvious chain rule through ``sFree`` — ``∂sFree/∂E = ½(−1 − δ/D)``, + ``∂sFree/∂S = ½(1 + (δ + 2·Km)/D)`` — is algebraically identical and cancels + catastrophically wherever ``δ < 0`` with ``|δ| ≫ √(4·Km·S)``, which is the + same regime that used to break ``sFree`` itself: measured against mpmath at + 60 digits, that grouping reached a relative error of 1e+10 in ∂rate/∂E on a + deep-saturation sweep *after* the root was fixed, while these forms stay at + machine precision. They fall out of differentiating the symmetric form of the + rate — the tQSSA complex is ``c = ½(A − D)`` with ``A = E + S + Km`` and the + same ``D = √(A² − 4·E·S)``, so ``∂c/∂E = (S − c)/D`` and ``S − c = sFree``. ``add(col, row, value_c, prefix)`` writes one accumulation of ``value_c`` into entry ``(row, col)`` of ∂f/∂x and returns the C line, or ``None`` if the @@ -4985,23 +5027,15 @@ def _mm_jacobian_groups(plan_mm, add) -> list[list[str]] | None: grp: list[str] = [] g = grp.append g(" {") - g( - f" double kcat = p[{int(mt['kcat_param_idx0'])}], " - f"Km = p[{int(mt['km_param_idx0'])}];" - ) - g(f" double E = y[{e}], S = y[{s}];") - g(" double delta = S - Km - E;") - g(" double Dmm = sqrt(delta*delta + 4.0*Km*S);") - g(" double sFree = 0.5*(delta + Dmm);") + # Km/E/S are declared by _mm_sfree_c_lines; only kcat is extra here. + g(f" double kcat = p[{int(mt['kcat_param_idx0'])}];") + for ln in _mm_sfree_c_lines(f"p[{int(mt['km_param_idx0'])}]", e, s, " "): + g(ln) g(" double dE = 0.0, dS = 0.0;") g(" if (sFree > 0.0 && Dmm > 0.0) {") - g(" double dsF_dE = 0.5*(-1.0 - delta/Dmm);") - g(" double dsF_dS = 0.5*(1.0 + (delta + 2.0*Km)/Dmm);") g(f" double Cmm = kcat * {_jac_c_float(mt['stat_factor'])};") - g(" double KpsF = Km + sFree;") - g(" double common = Cmm*E*Km/(KpsF*KpsF);") - g(" dE = Cmm*sFree/KpsF + common*dsF_dE;") - g(" dS = common*dsF_dS;") + g(" dE = Cmm*sFree/Dmm;") + g(" dS = Cmm*E*Km/((Km + sFree)*Dmm);") g(" }") for col, affected, deriv in ((e, mt["e_affected"], "dE"), (s, mt["s_affected"], "dS")): if not affected: diff --git a/python/bngsim/_jax_rhs.py b/python/bngsim/_jax_rhs.py index 652f746..1bf59f8 100644 --- a/python/bngsim/_jax_rhs.py +++ b/python/bngsim/_jax_rhs.py @@ -298,7 +298,15 @@ def rhs(y, t, params): if len(reactants) >= 2: e = y[reactants[0] - 1] s = y[reactants[1] - 1] - s_free = 0.5 * ((s - km - e) + jnp.sqrt((s - km - e) ** 2 + 4.0 * km * s)) + # Stable positive root of x² − delta·x − km·s = 0 (GH #89): + # ½(delta + D) cancels for delta < 0, so that branch uses the + # conjugate form. The denominator is masked before the divide + # so the unselected branch cannot inject a NaN into a tangent. + delta = s - km - e + d_mm = jnp.sqrt(delta * delta + 4.0 * km * s) + neg = delta < 0.0 + denom = jnp.where(neg, d_mm - delta, 1.0) + s_free = jnp.where(neg, 2.0 * km * s / denom, 0.5 * (delta + d_mm)) rate = sf * kcat * s_free * e / (km + s_free) else: rate = 0.0 diff --git a/python/tests/test_codegen_mm_sens.py b/python/tests/test_codegen_mm_sens.py index 2bd9f9f..efb96eb 100644 --- a/python/tests/test_codegen_mm_sens.py +++ b/python/tests/test_codegen_mm_sens.py @@ -6,34 +6,34 @@ modeller who writes one should not silently lose the analytic gradient. Unlike the Functional path there is no symbolic differentiation here. The tQSSA -rate is closed form, +rate is closed form — ``sFree`` is the positive root of ``x² − δ·x − Km·S = 0``, - δ = S − Km − E, D = √(δ² + 4·Km·S), sFree = ½(δ + D) + δ = S − Km − E, D = √(δ² + 4·Km·S) rate = kcat·stat·E·sFree/(Km + sFree) so both partials are written out: ∂rate/∂kcat = stat·E·sFree/(Km + sFree) ( = rate/kcat ) - ∂rate/∂Km = kcat·stat·E·(∂sFree/∂Km·Km − sFree)/(Km + sFree)² - with ∂sFree/∂Km = ½(−1 + (2S − δ)/D) + ∂rate/∂Km = −kcat·stat·E·sFree/((Km + sFree)·D) ( = −rate/D ) Both were checked against ``sympy.diff`` — symbolically and over random parameter points — before being written down, and ``test_the_partials_match_sympy`` keeps that check in the suite rather than in a notebook someone has to trust. -Two things are worth knowing about this file's tolerances. +**The grouping is load-bearing**, and that is the one thing to keep in mind when +touching this file. Every algebraically identical rewrite tried so far is worse +in float64, some catastrophically: -**The grouping is load-bearing.** An algebraically identical form of ∂/∂Km with -the ½ factors cancelled is 14 digits worse on log-spread parameters. What is -emitted mirrors the shipped Jacobian's own ∂rate/∂E grouping, and -``test_the_emitted_grouping_is_the_stable_one`` pins that. +* ``sFree = ½(δ + D)`` — the textbook root, and what shipped until GH #89 — + subtracts two nearly-equal positive numbers once ``δ < 0``, losing ~2 digits per + decade of ``|δ|/√(4·Km·S)`` and *every* digit by 1e8. It is now the conjugate + form ``2·Km·S/(D − δ)`` on that branch. +* The chain rule through ``∂sFree/∂Km = ½(−1 + (2S − δ)/D)`` cancels in the same + regime, so it stayed at 1e+10 relative error even after the root was fixed. +* ∂/∂Km with the ½ factors cancelled is 14 digits worse on log-spread parameters. -**There is an accuracy floor, and it is not this stage's.** ``sFree = ½(δ + D)`` -cancels catastrophically once ``|δ| ≫ √(4·Km·S)``, roughly two digits per decade. -That is shipped behaviour in the RHS and the analytical Jacobian; a model in that -regime has a wrong trajectory before it has a wrong gradient. So the oracles here -run in the well-conditioned regime, and the degenerate corner is asserted only to -the extent that ∂f/∂p is no worse than the rate it differentiates. +``TestNumericalStability`` pins all three against mpmath, so a future tidy-up +fails there instead of silently in a trajectory or a gradient. """ from __future__ import annotations @@ -113,6 +113,32 @@ def _has_cc() -> bool: """ +# Deliberately stiff root ratio: a large enzyme excess over a small Km·S puts +# |δ|/√(4·Km·S) at ~1e7, where the pre-#89 root and partials fall apart. Both +# corpus MM models sit at ratio ~1, so this is the only way to show the fix in +# the emitted C. Same numbers as tests/data/mm_tqssa_stiff.net, which the C++ +# suite uses for the interpreted RHS and the native Jacobian. +MM_STIFF = """\ +begin parameters + 1 kcat 2.0 # Constant + 2 Km 1.0e-8 # Constant +end parameters +begin species + 1 E() 9000 + 2 S() 20 + 3 P() 0 +end species +begin reactions + 1 1,2 1,3 MM kcat Km #_R1 +end reactions +begin groups + 1 Et 1 + 2 St 2 + 3 Pt 3 +end groups +""" + + def _model(tmp_path, text, name="m.net"): net = tmp_path / name net.write_text(text) @@ -122,17 +148,23 @@ def _model(tmp_path, text, name="m.net"): # ─── the derivative, against sympy ───────────────────────────────────────── -def _tqssa_partials(kcat, Km, E, S, stat=1.0): - """Exactly the arithmetic the emitted C performs.""" +def _tqssa_sfree(Km, E, S): + """The free substrate exactly as ``_mm_sfree_c_lines`` emits it.""" delta = S - Km - E D = math.sqrt(delta * delta + 4.0 * Km * S) - sF = 0.5 * (delta + D) + if delta >= 0.0: + return delta, D, 0.5 * (delta + D) + return delta, D, (2.0 * Km * S / (D - delta) if D - delta > 0.0 else 0.0) + + +def _tqssa_partials(kcat, Km, E, S, stat=1.0): + """Exactly the arithmetic the emitted C performs.""" + _delta, D, sF = _tqssa_sfree(Km, E, S) if not (sF > 0.0 and D > 0.0): return 0.0, 0.0 - dkcat = stat * E * sF / (Km + sF) - dsF_dKm = 0.5 * (-1.0 + (2.0 * S - delta) / D) KpsF = Km + sF - dKm = kcat * stat * E * (dsF_dKm * Km - sF) / (KpsF * KpsF) + dkcat = stat * E * sF / KpsF + dKm = -kcat * stat * E * sF / (KpsF * D) return dkcat, dKm @@ -150,18 +182,127 @@ def test_the_partials_match_sympy(self): rate = kcat * stat * E * sF / (Km + sF) mine_kcat = stat * E * sF / (Km + sF) - dsF_dKm = (-1 + (2 * S - delta) / D) / 2 - mine_Km = kcat * stat * E * (dsF_dKm * Km - sF) / (Km + sF) ** 2 + mine_Km = -kcat * stat * E * sF / ((Km + sF) * D) assert sp.simplify(sp.diff(rate, kcat) - mine_kcat) == 0 assert sp.simplify(sp.diff(rate, Km) - mine_Km) == 0 - def test_the_emitted_grouping_is_the_stable_one(self): - """Algebraic identity is not enough. The 'simplified' ∂/∂Km — the ½ - factors cancelled, ``2·kcat·E·(E−S−D+Km·B/D)/(A+D)²`` — is identical on - paper and worthless in float64 on log-spread parameters. This pins that - the *emitted* grouping is the one that survives, so a future tidy-up - that 'simplifies' it fails here instead of silently in a gradient.""" + # ...and the conjugate root the emitted C actually evaluates is the same + # sFree, so differentiating ½(δ + D) above describes what ships. + assert sp.simplify(2 * Km * S / (D - delta) - sF) == 0 + + def test_both_partials_vanish_where_the_rate_does(self): + """At S = 0 the rate is 0 for *every* kcat and Km (S = 0 forces + sFree = 0 whatever Km is), so the guard returning 0 is the correct + derivative, not a fallback.""" + for Km in (0.1, 1.0, 100.0): + assert _tqssa_partials(2.0, Km, 25.0, 0.0) == (0.0, 0.0) + + +# ─── the groupings, against mpmath (GH #89) ──────────────────────────────── +# +# Reference is mpmath at 60 decimal digits; every candidate below is evaluated in +# plain float64, which is exactly the arithmetic the emitted C performs. Symbolic +# identity is checked above and is not enough on its own — each rejected form +# here simplifies to zero against the one that ships. + + +def _ref(kcat, Km, E, S, stat=1.0): + """The tQSSA rate and its three partials at 60 digits.""" + import mpmath as mp + + with mp.workdps(60): + kcat, Km, E, S, stat = (mp.mpf(repr(v)) for v in (kcat, Km, E, S, stat)) + delta = S - Km - E + D = mp.sqrt(delta * delta + 4 * Km * S) + u = (delta + D) / 2 + rate = kcat * stat * E * u / (Km + u) + return { + "sFree": u, + "rate": rate, + "dE": kcat * stat * u / D, + "dS": kcat * stat * E * Km / ((Km + u) * D), + "dKm": -rate / D, + } + + +def _rel(got, want): + import mpmath as mp + + with mp.workdps(60): + return float(abs((mp.mpf(repr(got)) - want) / want)) if want != 0 else float(got != 0) + + +class TestNumericalStability: + """Each test states a rejected form, shows it is the same expression, and + measures both against mpmath. The point of the pairing is that no reviewer + can 'simplify' one of these back without the measurement disagreeing.""" + + def test_the_root_survives_a_huge_negative_delta(self): + """``½(δ + D)`` subtracts two nearly-equal positive numbers when δ < 0. + At ``|δ|/√(4·Km·S) = 1e8`` it has no correct digit left — and this is the + *rate*, so a trajectory in deep enzyme excess was wrong before any + derivative was.""" + Km = S = 1.0 + for exponent, floor in ((4, 1e-9), (8, 0.5)): + ratio = 10.0**exponent + E = ratio * math.sqrt(4.0 * Km * S) + S - Km # δ = −ratio·√(4·Km·S) + want = _ref(1.0, Km, E, S)["sFree"] + + delta = S - Km - E + D = math.sqrt(delta * delta + 4.0 * Km * S) + assert _rel(0.5 * (delta + D), want) > floor # the old root + assert _rel(_tqssa_sfree(Km, E, S)[2], want) < 1e-15 # what ships + + def test_the_root_is_bit_identical_for_a_non_negative_delta(self): + """The conjugate form is a branch, not a replacement: where ``½(δ + D)`` + is already cancellation-free it is still what evaluates, bit for bit.""" + for Km, E, S in ((50.0, 10.0, 100.0), (35.0, 25.0, 120.0), (1e-3, 1.0, 4.0)): + delta, D, sFree = _tqssa_sfree(Km, E, S) + assert delta >= 0.0 + assert sFree == 0.5 * (delta + D) + + def test_a_corpus_scale_ratio_moves_by_ulps_not_digits(self): + """Both MM models in the corpus *do* reach δ < 0, so they take the new + branch — the claim is not that nothing changed, it is that at their + ``|δ|/√(4·Km·S) ≤ 5`` the two forms differ at round-off scale. ~2 digits + per decade puts a ratio under 10 at ~1e-15 relative, against 2e-2 at the + 1e7 of the stiff fixture. Measured over the corpus trajectories + themselves, the largest ``sFree`` shift was 1.9e-16 and 6.0e-16.""" + for Km, E, S in ((3.0, 5.0, 1.0), (10.0, 40.0, 2.0), (0.5, 8.0, 1.5)): + delta, D, sFree = _tqssa_sfree(Km, E, S) + assert delta < 0.0 and abs(delta) / math.sqrt(4.0 * Km * S) < 10.0 + assert abs(sFree - 0.5 * (delta + D)) <= 1e-13 * sFree + + def test_the_partials_beat_the_chain_rule_through_sfree(self): + """Fixing the root alone is not enough. ``∂sFree/∂E = ½(−1 − δ/D)`` and + ``∂sFree/∂Km = ½(−1 + (2S − δ)/D)`` — the obvious chain rule, and what + shipped — cancel in exactly the regime the root did, so they stay wrong + after the root is right. Deep saturation, ``Km ≪ S,E``.""" + kcat, Km, E, S = 2.0, 1e-8, 9000.0, 20.0 + want = _ref(kcat, Km, E, S) + delta, D, sFree = _tqssa_sfree(Km, E, S) + KpsF = Km + sFree + common = kcat * E * Km / (KpsF * KpsF) + + # ∂rate/∂E: the chain rule through sFree, vs the emitted single quotient. + chained = kcat * sFree / KpsF + common * (0.5 * (-1.0 - delta / D)) + assert _rel(chained, want["dE"]) > 1e6 # ten orders of magnitude out + assert _rel(kcat * sFree / D, want["dE"]) < 1e-14 + + # ∂rate/∂Km: same story, same regime — 2.6% rather than 1e+10, because + # the surviving term is smaller, but still nothing a gradient can use. + chained_km = ( + kcat * E * ((0.5 * (-1.0 + (2.0 * S - delta) / D)) * Km - sFree) / (KpsF * KpsF) + ) + assert _rel(chained_km, want["dKm"]) > 1e-2 + assert _rel(_tqssa_partials(kcat, Km, E, S)[1], want["dKm"]) < 1e-14 + + def test_the_simplified_dkm_is_still_rejected(self): + """The other identical-on-paper ∂/∂Km — the ½ factors cancelled, + ``2·kcat·E·(E−S−D+Km·B/D)/(A+D)²`` — remains worthless in float64 on + log-spread parameters. Kept from #55: it is a different rewrite from the + chain rule above and fails for a different reason.""" kcat, Km, E, S = 29.9, 2.88e-4, 510.0, 2.4e-4 emitted = _tqssa_partials(kcat, Km, E, S)[1] @@ -172,13 +313,50 @@ def test_the_emitted_grouping_is_the_stable_one(self): # They are the same expression; they are not the same number. assert abs(simplified - emitted) > 0.1 * abs(emitted) - - def test_both_partials_vanish_where_the_rate_does(self): - """At S = 0 the rate is 0 for *every* kcat and Km (S = 0 forces - sFree = 0 whatever Km is), so the guard returning 0 is the correct - derivative, not a fallback.""" - for Km in (0.1, 1.0, 100.0): - assert _tqssa_partials(2.0, Km, 25.0, 0.0) == (0.0, 0.0) + assert _rel(emitted, _ref(kcat, Km, E, S)["dKm"]) < 1e-12 + + def test_every_emitted_quantity_holds_up_across_the_four_sweeps(self): + """The regime sweep from the issue, as an assertion rather than a table. + Worst relative error over each sweep, for every quantity the MM path + emits — rate, both Jacobian partials, both ∂f/∂p columns.""" + import random + + sweeps = { + "uniform": lambda r: [r.uniform(0.1, 10.0) for _ in range(4)], + "log-spread": lambda r: [10.0 ** r.uniform(-4, 3) for _ in range(4)], + "Km >> S,E": lambda r: [ + 10.0 ** r.uniform(-2, 2), + 10.0 ** r.uniform(3, 6), + 10.0 ** r.uniform(-2, 1), + 10.0 ** r.uniform(-2, 1), + ], + "Km << S,E": lambda r: [ + 10.0 ** r.uniform(-2, 2), + 10.0 ** r.uniform(-8, -4), + 10.0 ** r.uniform(1, 4), + 10.0 ** r.uniform(1, 4), + ], + } + for name, gen in sweeps.items(): + rng = random.Random(12345) # fixed: a flaky numerics test is useless + worst = dict.fromkeys(("rate", "dE", "dS", "dKm", "dkcat"), 0.0) + for _ in range(300): + kcat, Km, E, S = gen(rng) + want = _ref(kcat, Km, E, S) + _delta, D, sFree = _tqssa_sfree(Km, E, S) + KpsF = Km + sFree + dkcat, dKm = _tqssa_partials(kcat, Km, E, S) + got = { + "rate": kcat * E * sFree / KpsF, + "dE": kcat * sFree / D, + "dS": kcat * E * Km / (KpsF * D), + "dKm": dKm, + "dkcat": dkcat, + } + for k, v in got.items(): + ref_k = want[k] if k != "dkcat" else want["rate"] / kcat + worst[k] = max(worst[k], _rel(v, ref_k)) + assert max(worst.values()) < 1e-11, (name, worst) # ─── the gate ────────────────────────────────────────────────────────────── @@ -281,7 +459,10 @@ def test_both_rate_constants_get_a_column(self, tmp_path): names = [p["name"] for p in model._core.codegen_data()["parameters"]] for name in ("kcat", "Km"): assert f" case {names.index(name)}:" in src - assert "double sFree = 0.5*(delta + Dmm);" in src + # The free substrate is the stable root, not ½(δ + D) — see + # TestNumericalStability (GH #89). + assert "double sFree = (delta >= 0.0) ? 0.5*(delta + Dmm)" in src + assert "2.0*Km*S/(Dmm - delta)" in src def test_the_enzyme_row_is_not_scattered(self, tmp_path): """An MM enzyme is on both sides, so its net stoichiometry is 0. Emitting @@ -401,7 +582,9 @@ class TestAgainstFiniteDifference: """Central FD of the *emitted* RHS against the *emitted* ∂f/∂p, both called through ctypes on the same ``p[]``. No integrator, no tolerance tuning.""" - @pytest.mark.parametrize("text", [MM, MM_DERIVED]) + # Named, or pytest builds the id from the .net text and every -rA line in the + # CI log carries the whole fixture. + @pytest.mark.parametrize("text", [MM, MM_DERIVED], ids=["direct", "derived"]) def test_dfdp_matches_a_finite_difference_of_the_rhs(self, tmp_path, monkeypatch, text): model = _model(tmp_path, text) comp = _Compiled(model, tmp_path, monkeypatch) @@ -444,6 +627,49 @@ def test_dfdp_matches_a_finite_difference_of_the_rhs(self, tmp_path, monkeypatch assert worst < 1e-6, f"worst {worst:.3e} at {where}" +@requires_cc +class TestTheEmittedCAtAStiffRatio: + """The other tests in this file compare the emitted C to a finite difference + *of itself*, which is blind to a cancellation both sides share. This one + compiles the emitted RHS and ∂f/∂p and checks them against mpmath at a + ``|δ|/√(4·Km·S)`` of 1e7 — where the pre-#89 emission was 2% out on the rate + and 1e+10 out on ∂rate/∂E (GH #89).""" + + def test_the_emitted_rhs_and_dfdp_hold_at_ratio_1e7(self, tmp_path, monkeypatch): + model = _model(tmp_path, MM_STIFF) + comp = _Compiled(model, tmp_path, monkeypatch) + params = comp.data["parameters"] + names = [p["name"] for p in params] + base = [float(p["value"]) for p in params] + kcat, Km, E, S = 2.0, 1.0e-8, 9000.0, 20.0 + + # The regime really is the pathological one, not a mild one. + delta, D, _sFree = _tqssa_sfree(Km, E, S) + assert delta < 0.0 and abs(delta) / math.sqrt(4.0 * Km * S) > 1e6 + + want = _ref(kcat, Km, E, S) + # E + S -> E + P, so the enzyme row is 0 and P gains what S loses. + ydot = comp.f(0.0, [E, S, 0.0], base) + assert ydot[0] == 0.0 + assert _rel(-ydot[1], want["rate"]) < 1e-14 + assert _rel(ydot[2], want["rate"]) < 1e-14 + + for name, ref in (("kcat", want["rate"] / kcat), ("Km", want["dKm"])): + col = comp.dfdp(0.0, [E, S, 0.0], base, names.index(name)) + assert col[0] == 0.0 + assert _rel(col[2], ref) < 1e-13, name + assert _rel(-col[1], ref) < 1e-13, name + + # What the emission this replaced would have produced at this same point, + # so the numbers in the issue are reproducible from the suite. + old_sfree = 0.5 * (delta + D) + assert _rel(old_sfree, want["sFree"]) > 1e-3 + old_dE = kcat * old_sfree / (Km + old_sfree) + (kcat * E * Km / (Km + old_sfree) ** 2) * ( + 0.5 * (-1.0 - delta / D) + ) + assert _rel(old_dE, want["dE"]) > 1e9 + + @requires_cc class TestEndToEnd: def test_mm_sensitivities_match_a_resolved_trajectory(self, tmp_path): diff --git a/src/model.cpp b/src/model.cpp index 9196ebe..0db9ab1 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1603,9 +1603,10 @@ compute_rxn_rate(const Reaction &rxn, const std::vector ¶ms, cons double E = conc[e_si]; double S = conc[s_si]; - // tQSSA: compute free substrate - double delta = S - Km - E; - double sFree = 0.5 * (delta + std::sqrt(delta * delta + 4.0 * Km * S)); + // tQSSA: compute free substrate. mm_tqssa takes the conjugate + // root when delta < 0 — writing it as ½(delta + D) there loses + // ~2 digits per decade of |delta|/√(4·Km·S) (GH #89). + double sFree = mm_tqssa(Km, E, S).sFree; // Guard against numerical issues (sFree should be non-negative) if (sFree < 0.0) @@ -1841,7 +1842,8 @@ std::pair NetworkModel::emit_ssa_propensity_source_structure() // MichaelisMenten tQSSA — mirrors compute_rxn_rate's SSA MM branch exactly // (raw conc for E/S, no ssa_volume_factor): kcat=p[k0], Km=p[km0], - // E=x[e0], S=x[s0]; sFree = max(0, 0.5·((S−Km−E) + √((S−Km−E)²+4·Km·S))). + // E=x[e0], S=x[s0]; sFree = max(0, root of x²−d·x−Km·S with d = S−Km−E), + // taken through mm_tqssa's stable branch (GH #89). const bool mm_shape = rxn.rate_law_type == RateLawType::MichaelisMenten && rxn.rate_law_param_indices.size() >= 2 && rxn.reactant_indices.size() >= 2; @@ -1857,7 +1859,10 @@ std::pair NetworkModel::emit_ssa_propensity_source_structure() body += " { double Km = p[" + std::to_string(km0) + "]; double E = x[" + std::to_string(e0) + "]; double S = x[" + std::to_string(s0) + "];\n"; body += " double d = S - Km - E;\n"; - body += " double sf = 0.5 * (d + sqrt(d * d + 4.0 * Km * S));\n"; + body += " double D = sqrt(d * d + 4.0 * Km * S);\n"; + body += " double sf = (d >= 0.0) ? 0.5 * (d + D)\n"; + body += + " : ((D - d) > 0.0 ? 2.0 * Km * S / (D - d) : 0.0);\n"; body += " if (sf < 0.0) sf = 0.0;\n"; body += " a[" + std::to_string(r) + "] = p[" + std::to_string(kcat0) + "]"; if (rxn.stat_factor != 1.0) diff --git a/tests/data/mm_tqssa_stiff.net b/tests/data/mm_tqssa_stiff.net new file mode 100644 index 0000000..d3be347 --- /dev/null +++ b/tests/data/mm_tqssa_stiff.net @@ -0,0 +1,36 @@ +# Michaelis-Menten tQSSA at a deliberately stiff root ratio — GH #89 regression. +# +# Same rate law as mm_tqssa.net, but a large enzyme excess over a small Km*S: +# delta = S - Km - E = 20 - 1e-8 - 9000 = -8980.00000001 +# |delta| / sqrt(4*Km*S) ~= 1.0e7 +# +# In that regime sFree = 0.5*(delta + sqrt(delta^2 + 4*Km*S)) subtracts two +# nearly-equal positive numbers. Reference values at 50 decimal digits (mpmath), +# against what plain float64 produced before the conjugate root landed: +# +# sFree = 2.2271714922024141071e-11 old: 2.1e-02 relative error +# rate = 39.99999999995545657 old: 2.1e-02 +# dE = 4.9602928556791156198e-15 old: 7.7e+09 (and the wrong SIGN) +# dS = 1.9999999999977678682 old: 5.6e-05 +# +# The corpus's two MM models both sit at ratio ~1, so nothing there moves; this +# fixture is the honest way to show the fix does something. +begin parameters + 1 kcat 2.0 + 2 Km 1.0e-8 + 3 E_0 9000.0 + 4 S_0 20.0 +end parameters +begin species + 1 E() E_0 + 2 S() S_0 + 3 P() 0 +end species +begin reactions + 1 1,2 1,3 MM kcat Km #Rule1 +end reactions +begin groups + 1 Etot 1 + 2 Stot 2 + 3 Ptot 3 +end groups diff --git a/tests/test_bngsim.cpp b/tests/test_bngsim.cpp index 20a34bd..8faef42 100644 --- a/tests/test_bngsim.cpp +++ b/tests/test_bngsim.cpp @@ -5,6 +5,7 @@ #include #include +#include // test_mm_tqssa_stiff_root calls mm_tqssa directly #include #include @@ -947,6 +948,71 @@ int test_mm_tqssa() { return 0; } +// ═══════════════════════════════════════════════════════════════════════════════ +// Test: MM tQSSA at a stiff root ratio (GH #89) +// +// mm_tqssa.net sits at |delta|/sqrt(4*Km*S) ~ 1, where every way of writing the +// root agrees; so does every MM model in the corpus. mm_tqssa_stiff.net sits at +// ~1e7, where the textbook 0.5*(delta + D) has lost all but two digits of the +// RATE and the chain rule through sFree has the wrong SIGN for d(rate)/dE. +// +// Reference values are mpmath at 50 decimal digits (see the fixture header). They +// are hard-coded rather than recomputed in long double, which is 64-bit on arm64 +// and would silently become a self-comparison. +// ═══════════════════════════════════════════════════════════════════════════════ + +int test_mm_tqssa_stiff_root() { + // 50-digit references, rounded to double. See tests/data/mm_tqssa_stiff.net. + const double ref_sfree = 2.2271714922024141071e-11; + const double ref_rate = 39.99999999995545657; + const double ref_dE = 4.9602928556791156198e-15; + const double ref_dS = 1.9999999999977678682; + + auto rel = [](double got, double want) { return std::abs(got - want) / std::abs(want); }; + + // 1. The shared root helper, called directly. + { + auto q = bngsim::mm_tqssa(1.0e-8, 9000.0, 20.0); + CHECK(rel(q.sFree, ref_sfree) < 1e-14, + "mm_tqssa root at ratio 1e7 (rel=" + std::to_string(rel(q.sFree, ref_sfree)) + ")"); + // The form it replaced, evaluated here so the regression is visible in + // the failure message rather than only in the issue. + double delta = 20.0 - 1.0e-8 - 9000.0; + double textbook = 0.5 * (delta + std::sqrt(delta * delta + 4.0 * 1.0e-8 * 20.0)); + CHECK(rel(textbook, ref_sfree) > 1e-3, "the textbook root really is wrong here"); + } + + // 2. The interpreted RHS. Reaction is E + S -> E + P, so ydot = (0, -rate, +rate). + auto model = bngsim::NetworkModel::from_net(data_path("mm_tqssa_stiff.net")); + CHECK(model.n_species() == 3, "Expected 3 species"); + + std::vector y = {9000.0, 20.0, 0.0}, f(3); + model.compute_derivs(0.0, y.data(), f.data()); + CHECK(rel(-f[1], ref_rate) < 1e-14, + "interpreted MM rate at ratio 1e7 (rel=" + std::to_string(rel(-f[1], ref_rate)) + ")"); + CHECK_CLOSE(f[2], -f[1], 1e-12, "P gains exactly what S loses"); + CHECK_CLOSE(f[0], 0.0, 1e-30, "the enzyme is conserved"); + + // 3. The closed-form analytical Jacobian at the same point. d(rate)/dE is + // ~1e-15 against a d(rate)/dS of ~2, so this is where the old grouping + // returned -3.8e-05 — nine orders out, and negative. + const auto& ajd = model.analytical_jacobian(); + CHECK(ajd.available && ajd.mm_reactions.size() == 1, "MM analytical Jac built"); + std::vector J(9); + model.fill_dense_analytical_jacobian(0.0, y.data(), J.data()); + double dE = J[0 * 3 + 2]; // d(ydot[P])/d(E) = +d(rate)/dE + double dS = J[1 * 3 + 2]; // d(ydot[P])/d(S) = +d(rate)/dS + CHECK(dE > 0.0, "d(rate)/dE is positive — more enzyme cannot slow the reaction"); + CHECK(rel(dE, ref_dE) < 1e-12, + "MM d(rate)/dE at ratio 1e7 (rel=" + std::to_string(rel(dE, ref_dE)) + ")"); + CHECK(rel(dS, ref_dS) < 1e-14, + "MM d(rate)/dS at ratio 1e7 (rel=" + std::to_string(rel(dS, ref_dS)) + ")"); + CHECK_CLOSE(J[0 * 3 + 1], -dE, 1e-30, "S row mirrors the P row"); + CHECK_CLOSE(J[1 * 3 + 1], -dS, 1e-30, "S row mirrors the P row"); + + return 0; +} + // ═══════════════════════════════════════════════════════════════════════════════ // Test: TFUN time-indexed — load .net with tfun(), simulate, verify interpolation // ═══════════════════════════════════════════════════════════════════════════════ @@ -2142,6 +2208,7 @@ int main() { RUN_TEST(test_sat_loads_with_rewrite_warning); RUN_TEST(test_hill_loads_with_rewrite_warning); RUN_TEST(test_mm_tqssa); + RUN_TEST(test_mm_tqssa_stiff_root); // Table functions (ADR-001 §3.12) RUN_TEST(test_tfun_time_indexed);