fix(mm): stable tQSSA free-substrate root and subtraction-free partials (#89)

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

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