From 22c62c9250da9880a7a9e7e5c2046eb6e3eeb661 Mon Sep 17 00:00:00 2001 From: parkyr Date: Mon, 13 Jul 2026 10:39:44 -0400 Subject: [PATCH 1/2] Fix reversible reaction Jacobian for issue #23 --- PharmaPy/Kinetics.py | 55 ++++++- .../test_reversible_jacobian_reverse_term.py | 137 ++++++++++++++++++ 2 files changed, 189 insertions(+), 3 deletions(-) create mode 100644 tests/test_reversible_jacobian_reverse_term.py diff --git a/PharmaPy/Kinetics.py b/PharmaPy/Kinetics.py index 9a59366..ff27103 100644 --- a/PharmaPy/Kinetics.py +++ b/PharmaPy/Kinetics.py @@ -532,12 +532,52 @@ 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): + is_product = self.stoich_matrix > 0 + orders = abs(is_product * self.stoich_matrix) + conc = np.asarray(conc) + conc_correc = np.maximum(eps, conc) + + keq_temp = self.equil_term(temp, self.delta_hrxn) + 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) @@ -567,8 +607,17 @@ def derivatives(self, conc, temp, dstates=True): 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) + + 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 diff --git a/tests/test_reversible_jacobian_reverse_term.py b/tests/test_reversible_jacobian_reverse_term.py new file mode 100644 index 0000000..fc81837 --- /dev/null +++ b/tests/test_reversible_jacobian_reverse_term.py @@ -0,0 +1,137 @@ +# -*- 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 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 _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 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) From 64210286966b6b8c7ed7a02bdede443b92deb855 Mon Sep 17 00:00:00 2001 From: parkyr Date: Mon, 13 Jul 2026 11:44:51 -0400 Subject: [PATCH 2/2] Thread runtime heat of reaction into reversible Jacobian --- PharmaPy/Kinetics.py | 14 +- PharmaPy/Reactors.py | 17 ++- .../test_reversible_jacobian_reverse_term.py | 131 ++++++++++++++++++ 3 files changed, 154 insertions(+), 8 deletions(-) diff --git a/PharmaPy/Kinetics.py b/PharmaPy/Kinetics.py index ff27103..c466def 100644 --- a/PharmaPy/Kinetics.py +++ b/PharmaPy/Kinetics.py @@ -544,13 +544,16 @@ def elem_df_dstates(self, conc): return df_dconc - def _reverse_df_dstates(self, conc, temp): + 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) - keq_temp = self.equil_term(temp, self.delta_hrxn) + 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: @@ -601,14 +604,15 @@ 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) if self.keq_params is not None: - df_dstates = df_dstates - self._reverse_df_dstates(conc, temp) + df_dstates = df_dstates - self._reverse_df_dstates( + conc, temp, delta_hrxn) dr_dstates = df_dstates * temp_terms[..., np.newaxis] @@ -643,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: diff --git a/PharmaPy/Reactors.py b/PharmaPy/Reactors.py index d41bdb5..6d2dce5 100644 --- a/PharmaPy/Reactors.py +++ b/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 diff --git a/tests/test_reversible_jacobian_reverse_term.py b/tests/test_reversible_jacobian_reverse_term.py index fc81837..196c6f7 100644 --- a/tests/test_reversible_jacobian_reverse_term.py +++ b/tests/test_reversible_jacobian_reverse_term.py @@ -11,6 +11,9 @@ 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 @@ -38,6 +41,19 @@ def _reversible_kinetics(data_path): ) +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), @@ -48,6 +64,33 @@ def _irreversible_kinetics(data_path): ) +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) @@ -135,3 +178,91 @@ def net_rates(c): 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)