Fix reversible reaction Jacobian

PharmaPy/Kinetics.py

@@ -532,12 +532,55 @@ def equilibrium_model(self, conc, temp, deltah_rxn):
 
     def elem_df_dstates(self, conc):
 
+        conc = np.asarray(conc)
         f_term = self.elem_f_model(conc, self.params_f)
 
-        df_dconc = f_term[..., np.newaxis] * self.params_f / (conc + eps)
+        if conc.ndim == 1:
+            conc_term = conc + eps
+        else:
+            conc_term = conc[:, np.newaxis, :] + eps
+
+        df_dconc = f_term[..., np.newaxis] * self.params_f / conc_term
 
         return df_dconc
 
+    def _reverse_df_dstates(self, conc, temp, deltah_rxn=None):
+        is_product = self.stoich_matrix > 0
+        orders = abs(is_product * self.stoich_matrix)
+        conc = np.asarray(conc)
+        conc_correc = np.maximum(eps, conc)
+
+        if deltah_rxn is None:
+            deltah_rxn = self.delta_hrxn
+
+        keq_temp = self.equil_term(temp, deltah_rxn)
+        keq_temp[keq_temp == 0] = 1e20
+
+        if conc.ndim == 1:
+            r_term = np.zeros(self.num_rxns)
+            for ind in range(self.num_rxns):
+                r_term[ind] = np.prod(
+                    conc_correc**orders[ind]) / keq_temp[ind]
+
+            dr_dconc = r_term[:, np.newaxis] * orders / (conc + eps)
+        else:
+            n_conc = len(conc)
+            r_term = np.zeros((n_conc, self.num_rxns))
+            for ind in range(self.num_rxns):
+                if keq_temp.ndim == 1:
+                    keq_rxn = keq_temp[ind]
+                else:
+                    keq_rxn = keq_temp[:, ind]
+
+                r_term[:, ind] = np.prod(
+                    conc_correc**orders[ind], axis=1) / keq_rxn
+
+            dr_dconc = (
+                r_term[..., np.newaxis] * orders /
+                (conc[:, np.newaxis, :] + eps))
+
+        return dr_dconc
+
     def elem_df_dtheta(self, conc):
 
         f_term = self.elem_f_model(conc, self.params_f)
@@ -561,14 +604,24 @@ def elem_df_dtheta(self, conc):
 
         return drate_dorder
 
-    def derivatives(self, conc, temp, dstates=True):
+    def derivatives(self, conc, temp, dstates=True, delta_hrxn=None):
         temp_terms = self.temp_term(temp)
         f_terms = self.kinetic_model(conc, self.params_f, *self.args_kin)
 
         if dstates:  # --------------- wrt states
             df_dstates = self.df_dstates(conc, *self.args_kin)
-            dr_dstates = (temp_terms * df_dstates.T).T
-            jac_states = np.dot(self.normalized_stoich, dr_dstates)
+            if self.keq_params is not None:
+                df_dstates = df_dstates - self._reverse_df_dstates(
+                    conc, temp, delta_hrxn)
+
+            dr_dstates = df_dstates * temp_terms[..., np.newaxis]
+
+            if dr_dstates.ndim == 2:
+                jac_states = np.dot(self.normalized_stoich, dr_dstates)
+            else:
+                jac_states = np.einsum(
+                    'sr,trc->tsc', self.normalized_stoich, dr_dstates)
+
             return jac_states
         else:  # --------------- wrt parameters
             dk_dphi = self.dk_dkparams(temp).T
@@ -594,7 +647,7 @@ def get_rxn_rates(self, conc, temp=298.15, overall_rates=True, jac=False,
                       delta_hrxn=None):
 
         if jac:
-            jac_states = self.derivatives(conc, temp)
+            jac_states = self.derivatives(conc, temp, delta_hrxn=delta_hrxn)
             return jac_states
 
         else:

PharmaPy/Reactors.py

@@ -429,19 +429,30 @@ def get_jacobians(self, time, states, sw, sens, params, wrt_states=True):
         conc = states[:num_species]
         if self.isothermal:
             temp = self.Liquid_1.temp
-            jac_states = self.Kinetics.derivatives(conc, temp)
         else:
             temp = states[num_species]
             jac_states = np.zeros((num_states, num_states))
 
-            jac_r_kin = self.Kinetics.derivatives(conc, temp)
+        if self.Kinetics.keq_params is None:
+            deltah_rxn = None
+        else:
+            deltah_rxn = self.Liquid_1.getHeatOfRxn(
+                self.Kinetics.stoich_matrix, temp, self.mask_species,
+                self.Kinetics.delta_hrxn, self.Kinetics.tref_hrxn)
+
+        jac_r_kin = self.Kinetics.derivatives(
+            conc, temp, delta_hrxn=deltah_rxn)
+
+        if self.isothermal:
+            jac_states = jac_r_kin
+        else:
             jac_states[:num_states - 1, :num_states - 1] = jac_r_kin
 
         if wrt_states:
             return jac_states
         else:  # ---------- w.r.t params
             jac_theta_kin = self.Kinetics.derivatives(
-                conc, temp, dstates=False)
+                conc, temp, dstates=False, delta_hrxn=deltah_rxn)
 
             if self.isothermal:
                 jac_params = jac_theta_kin

tests/test_reversible_jacobian_reverse_term.py

@@ -0,0 +1,268 @@
+# -*- coding: utf-8 -*-
+"""Regression test for issue #23.
+
+For a reversible reaction the analytical state Jacobian returned by
+``RxnKinetics.get_rxn_rates(..., jac=True)`` (``derivatives``) must equal the
+finite-difference Jacobian of the net rate ``get_rxn_rates(..., jac=False)``.
+
+``derivatives`` differentiates only the forward term (``self.df_dstates`` ->
+``elem_df_dstates``), so for any reversible (``keq_params``) reaction the
+reverse-term contribution ``-d/dC[prod(C_prod**beta)/Keq]`` is missing from the
+Jacobian. This is pure NumPy (no solver backend), so no assimulo marker.
+"""
+
+import sys
+from types import ModuleType, SimpleNamespace
+
+import numpy as np
+
+from PharmaPy.Kinetics import RxnKinetics
+
+
+# A <-> B (reversible); C and solv are inert padding species in the database.
+STOICH = [[-1, 1, 0, 0]]
+PARTIC = ["A", "B", "C", "solv"]
+
+
+def _db(data_path):
+    return str(data_path["integration"] / "pfr_test_pure_comp.json")
+
+
+def _reversible_kinetics(data_path):
+    return RxnKinetics(
+        path=_db(data_path),
+        k_params=[1.0],        # kf = 1 at temp_ref = inf
+        ea_params=[0.0],
+        stoich_matrix=STOICH,
+        partic_species=PARTIC,
+        keq_params=[2.0],
+        delta_hrxn=[0.0],      # Keq constant = 2, independent of temperature
+        tref_hrxn=298.15,
+    )
+
+
+def _temperature_sensitive_reversible_kinetics(data_path):
+    return RxnKinetics(
+        path=_db(data_path),
+        k_params=[1.0],
+        ea_params=[0.0],
+        stoich_matrix=STOICH,
+        partic_species=PARTIC,
+        keq_params=[2.0],
+        delta_hrxn=[-5.0e4],
+        tref_hrxn=298.15,
+    )
+
+
+def _irreversible_kinetics(data_path):
+    return RxnKinetics(
+        path=_db(data_path),
+        k_params=[1.0],
+        ea_params=[0.0],
+        stoich_matrix=STOICH,
+        partic_species=PARTIC,
+    )
+
+
+def _stub_assimulo_modules(monkeypatch):
+    assimulo = ModuleType("assimulo")
+
+    solvers = ModuleType("assimulo.solvers")
+    solvers.CVode = object
+    solvers.LSODAR = object
+
+    problem = ModuleType("assimulo.problem")
+    problem.Explicit_Problem = object
+
+    monkeypatch.setitem(sys.modules, "assimulo", assimulo)
+    monkeypatch.setitem(sys.modules, "assimulo.solvers", solvers)
+    monkeypatch.setitem(sys.modules, "assimulo.problem", problem)
+
+
+def _import_base_reactor(monkeypatch):
+    try:
+        from PharmaPy.Reactors import _BaseReactor
+    except ModuleNotFoundError as exc:
+        if exc.name != "assimulo":
+            raise
+        _stub_assimulo_modules(monkeypatch)
+        from PharmaPy.Reactors import _BaseReactor
+
+    return _BaseReactor
+
+
+def test_reversible_jacobian_includes_reverse_term(data_path):
+    kin = _reversible_kinetics(data_path)
+
+    temp = 298.15
+    conc = np.array([1.0, 1.0, 1.0, 1.0])   # reverse rate nonzero at C_B = 1
+    deltah = np.array([0.0])
+
+    analytical = np.asarray(
+        kin.get_rxn_rates(conc, temp, overall_rates=True, jac=True)
+    )
+
+    # Central finite-difference Jacobian of the net species rates.
+    def net_rates(c):
+        return np.asarray(
+            kin.get_rxn_rates(c, temp, overall_rates=True, delta_hrxn=deltah)
+        )
+
+    n = conc.size
+    fd = np.zeros((n, n))
+    h = 1e-6
+    for j in range(n):
+        cp = conc.copy(); cp[j] += h
+        cm = conc.copy(); cm[j] -= h
+        fd[:, j] = (net_rates(cp) - net_rates(cm)) / (2 * h)
+
+    # Hand-computed truth for the participating (A, B) block: Keq = 2, kf = 1.
+    # net = [-(C_A - C_B/2), +(C_A - C_B/2), 0, 0]
+    #   d/dC_A = [-1, 1];  d/dC_B = [+0.5, -0.5]
+    np.testing.assert_allclose(fd[0, 1], 0.5, atol=1e-5)
+    np.testing.assert_allclose(fd[1, 1], -0.5, atol=1e-5)
+
+    # The analytical Jacobian must reproduce the reverse-term B column.
+    np.testing.assert_allclose(analytical, fd, atol=1e-5)
+
+
+def test_irreversible_forward_jacobian_is_unchanged(data_path):
+    kin = _irreversible_kinetics(data_path)
+
+    analytical = np.asarray(
+        kin.get_rxn_rates(
+            np.array([1.0, 1.0, 1.0, 1.0]),
+            298.15,
+            overall_rates=True,
+            jac=True,
+        )
+    )
+
+    expected = np.array([
+        [-1.0, 0.0, 0.0, 0.0],
+        [1.0, 0.0, 0.0, 0.0],
+        [0.0, 0.0, 0.0, 0.0],
+        [0.0, 0.0, 0.0, 0.0],
+    ])
+
+    np.testing.assert_allclose(analytical, expected, atol=1e-12)
+
+
+def test_reversible_jacobian_supports_batched_concentrations(data_path):
+    kin = _reversible_kinetics(data_path)
+
+    temp = np.array([298.15, 310.0])
+    conc = np.array([
+        [1.0, 1.0, 1.0, 1.0],
+        [2.0, 0.5, 1.0, 1.0],
+    ])
+    deltah = np.array([0.0])
+
+    analytical = np.asarray(
+        kin.get_rxn_rates(conc, temp, overall_rates=True, jac=True)
+    )
+
+    def net_rates(c):
+        return np.asarray(
+            kin.get_rxn_rates(c, temp, overall_rates=True, delta_hrxn=deltah)
+        )
+
+    n_times, n_species = conc.shape
+    fd = np.zeros((n_times, n_species, n_species))
+    h = 1e-6
+    for i in range(n_times):
+        for j in range(n_species):
+            cp = conc.copy(); cp[i, j] += h
+            cm = conc.copy(); cm[i, j] -= h
+            fd[i, :, j] = (net_rates(cp) - net_rates(cm))[i] / (2 * h)
+
+    assert analytical.shape == fd.shape
+    np.testing.assert_allclose(analytical, fd, atol=1e-5)
+
+
+def test_reversible_jacobian_uses_runtime_delta_hrxn(data_path):
+    kin = _temperature_sensitive_reversible_kinetics(data_path)
+
+    temp = 350.0
+    conc = np.array([1.0, 1.0, 1.0, 1.0])
+    runtime_deltah = np.array([-4.0e4])
+
+    analytical = np.asarray(
+        kin.get_rxn_rates(
+            conc,
+            temp,
+            overall_rates=True,
+            jac=True,
+            delta_hrxn=runtime_deltah,
+        )
+    )
+
+    def net_rates(c):
+        return np.asarray(
+            kin.get_rxn_rates(
+                c,
+                temp,
+                overall_rates=True,
+                delta_hrxn=runtime_deltah,
+            )
+        )
+
+    n = conc.size
+    fd = np.zeros((n, n))
+    h = 1e-6
+    for j in range(n):
+        cp = conc.copy(); cp[j] += h
+        cm = conc.copy(); cm[j] -= h
+        fd[:, j] = (net_rates(cp) - net_rates(cm)) / (2 * h)
+
+    np.testing.assert_allclose(analytical, fd, rtol=1e-5, atol=1e-5)
+
+
+def test_reactor_jacobian_passes_runtime_delta_hrxn(monkeypatch):
+    base_reactor = _import_base_reactor(monkeypatch)
+    runtime_deltah = np.array([-4.0e4])
+
+    class CaptureKinetics:
+        num_species = 2
+        keq_params = np.array([2.0])
+        stoich_matrix = np.array([[-1.0, 1.0]])
+        delta_hrxn = np.array([-5.0e4])
+        tref_hrxn = 298.15
+
+        def __init__(self):
+            self.calls = []
+
+        def derivatives(self, conc, temp, dstates=True, delta_hrxn=None):
+            self.calls.append((conc.copy(), temp, dstates, delta_hrxn.copy()))
+            return np.eye(2)
+
+    class CaptureLiquid:
+        temp = 350.0
+
+        def __init__(self):
+            self.calls = []
+
+        def getHeatOfRxn(self, stoich_matrix, temp, mask, heat_rxn_ref,
+                         tref_hrxn):
+            self.calls.append(
+                (stoich_matrix.copy(), temp, mask.copy(),
+                 heat_rxn_ref.copy(), tref_hrxn))
+            return runtime_deltah
+
+    kinetics = CaptureKinetics()
+    liquid = CaptureLiquid()
+    reactor = SimpleNamespace(
+        Kinetics=kinetics,
+        Liquid_1=liquid,
+        isothermal=True,
+        mask_species=np.array([True, True]),
+    )
+
+    jac = base_reactor.get_jacobians(
+        reactor, time=0.0, states=np.array([1.0, 1.0]), sw=None,
+        sens=None, params=None)
+
+    np.testing.assert_allclose(jac, np.eye(2))
+    assert len(liquid.calls) == 1
+    assert len(kinetics.calls) == 1
+    np.testing.assert_allclose(kinetics.calls[0][3], runtime_deltah)