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);