From 1066d6729cdf24cc673c4e02a7da1e85a22539cf Mon Sep 17 00:00:00 2001 From: parkyr Date: Fri, 10 Jul 2026 12:33:22 -0400 Subject: [PATCH] Fix kinetics custom params and equilibrium indexing --- PharmaPy/Kinetics.py | 45 ++++-- tests/test_kinetics_custom_and_equilibrium.py | 151 ++++++++++++++++++ 2 files changed, 185 insertions(+), 11 deletions(-) create mode 100644 tests/test_kinetics_custom_and_equilibrium.py diff --git a/PharmaPy/Kinetics.py b/PharmaPy/Kinetics.py index 9a59366..74770fd 100644 --- a/PharmaPy/Kinetics.py +++ b/PharmaPy/Kinetics.py @@ -356,11 +356,15 @@ def set_params(self, params): # From 1D params to matrix-shaped params orders = orders[np.newaxis, ...] params_f = orders - elif params_f is None: - raise RuntimeError("For user-defined kinetic function, " - "argument 'params_f' is mandatory.") + else: + params_f = params.get('params_f', None) + if params_f is None: + raise RuntimeError("For user-defined kinetic function, " + "argument 'params_f' is mandatory.") + params_f = np.asarray(params_f) + self.params_f_shape = params_f.shape - self.params_f = orders + self.params_f = params_f self.order_map = self.stoich_matrix < 0 else: @@ -373,6 +377,9 @@ def set_params(self, params): # From 1D params to matrix-shaped params dtype=np.float64) self.params_f[self.order_map] = params[self.num_paramsk:] + else: + params_f = np.asarray(params[self.num_paramsk:]) + self.params_f = params_f.reshape(self.params_f_shape) def nomenclature(self, stoich_matrix, kvals): @@ -387,9 +394,14 @@ def nomenclature(self, stoich_matrix, kvals): name_e = ['E_{a, %i}' % ind for ind in range(1, num_kpar + 1)] if self.fit_paramsf: - num_orders = (stoich_matrix < 0).sum() - name_orders = [r'\alpha_{}'.format(ind) - for ind in range(1, num_orders + 1)] + if self.elem_flag: + num_orders = (stoich_matrix < 0).sum() + name_orders = [r'\alpha_{}'.format(ind) + for ind in range(1, num_orders + 1)] + else: + num_orders = np.asarray(self.params_f).size + name_orders = [r'\theta_{f,%i}' % ind + for ind in range(1, num_orders + 1)] else: name_orders = [] @@ -424,7 +436,8 @@ def concat_params(self): params_concat = params_k_conc else: - params_concat = np.concatenate((params_k_conc, self.params_f)) + params_concat = np.concatenate( + (params_k_conc, np.asarray(self.params_f).ravel())) return params_concat @@ -452,13 +465,18 @@ def temp_term(self, temp): def equil_term(self, temp, deltah_temp): temp = np.asarray(temp) + deltah_temp = np.asarray(deltah_temp) inv_temp = (1/temp - 1/self.tref_hrxn) if temp.ndim == 0: k_eq = self.keq_params * np.exp(-deltah_temp/gas_ct * inv_temp) else: - k_eq = self.keq_params * \ - np.exp(-np.outer(inv_temp, deltah_temp/gas_ct)) + if deltah_temp.ndim <= 1: + exponent = np.outer(inv_temp, deltah_temp/gas_ct) + else: + exponent = inv_temp[:, np.newaxis] * deltah_temp/gas_ct + + k_eq = self.keq_params * np.exp(-exponent) return k_eq @@ -523,8 +541,13 @@ def equilibrium_model(self, conc, temp, deltah_rxn): else: 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**(orders[ind]), axis=1) / keq_temp[ind] + conc**(orders[ind]), axis=1) / keq_rxn overall_rate = f_term - r_term diff --git a/tests/test_kinetics_custom_and_equilibrium.py b/tests/test_kinetics_custom_and_equilibrium.py new file mode 100644 index 0000000..9849ad5 --- /dev/null +++ b/tests/test_kinetics_custom_and_equilibrium.py @@ -0,0 +1,151 @@ +# -*- coding: utf-8 -*- +"""Regression tests for issue #71. + +Two independent RxnKinetics defects, grouped as one kinetics-correctness fix: + +1. Constructing ``RxnKinetics`` with a user-defined ``kinetic_model`` + (``elem_flag=False``) raises ``UnboundLocalError`` in ``set_params``. +2. The vectorized (array-``temp``) equilibrium reverse term divides by ``Keq`` + indexed on the time axis instead of the reaction axis, so every time point + reuses the first row's ``Keq``. + +Both live in ``PharmaPy.Kinetics``, which is pure NumPy (no solver backend), +so these tests need no assimulo marker. +""" + +import numpy as np +import pytest + +from PharmaPy.Kinetics import RxnKinetics + + +STOICH = [[-1, -1, 1, 0]] # A + B -> C +PARTIC = ["A", "B", "C", "solv"] + + +def _db(data_path): + return str(data_path["integration"] / "pfr_test_pure_comp.json") + + +def _custom_rate(conc, params_f): + conc = np.maximum(1e-15, conc) + return np.exp(np.dot(np.log(conc), np.atleast_2d(params_f).T)) + + +def test_custom_kinetic_model_constructs(data_path): + """Defect 1: the documented custom-rate-law workflow must construct.""" + + # Must not raise UnboundLocalError in set_params. + kin = RxnKinetics( + path=_db(data_path), + k_params=[1.0], + ea_params=[0.0], + stoich_matrix=STOICH, + partic_species=PARTIC, + kinetic_model=_custom_rate, + params_f=np.array([[1.0, 1.0]]), + ) + + assert kin.elem_flag is False + assert kin.kinetic_model is _custom_rate + # The user-supplied f-parameters must be stored, not dropped. + assert kin.params_f is not None + np.testing.assert_allclose(np.asarray(kin.params_f).ravel(), [1.0, 1.0]) + + +def test_custom_kinetic_model_requires_params_f(data_path): + """Custom rate laws still require explicit f-parameters.""" + with pytest.raises(RuntimeError, match="params_f"): + RxnKinetics( + path=_db(data_path), + k_params=[1.0], + ea_params=[0.0], + stoich_matrix=STOICH, + partic_species=PARTIC, + kinetic_model=_custom_rate, + params_f=None, + ) + + +def test_custom_kinetic_model_flat_set_params_updates_params_f(data_path): + """Parameter-estimation vector updates custom model f-parameters.""" + kin = RxnKinetics( + path=_db(data_path), + k_params=[1.0], + ea_params=[0.0], + stoich_matrix=STOICH, + partic_species=PARTIC, + kinetic_model=_custom_rate, + params_f=np.array([[1.0, 1.0]]), + ) + + params = kin.concat_params() + params[-2:] = [2.0, 3.0] + + kin.set_params(params) + + np.testing.assert_allclose(np.asarray(kin.params_f).ravel(), [2.0, 3.0]) + + +def test_vectorized_equilibrium_matches_scalar_path(data_path): + """Defect 2: array-temp reverse term must use per-reaction Keq. + + Oracle is the trusted scalar-temp / 1-D branch: evaluating each time point + separately and stacking must equal the single vectorized 2-D evaluation. + """ + kin = 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=[-5e3], + tref_hrxn=298.15, + ) + + temp = np.array([290.0, 310.0, 330.0]) + conc = np.array([[1.0, 1.0, 1.0, 1.0], + [2.0, 2.0, 2.0, 2.0], + [0.5, 0.5, 0.5, 0.5]]) + deltah = np.atleast_1d(-5e3) + + vectorized = np.asarray(kin.equilibrium_model(conc, temp, deltah)) + per_time = np.array([ + np.asarray( + kin.equilibrium_model(conc[i], float(temp[i]), deltah) + ).ravel() + for i in range(temp.size) + ]) + + np.testing.assert_allclose(vectorized.reshape(per_time.shape), per_time) + + +def test_vectorized_equilibrium_accepts_2d_deltah(data_path): + """Array-temp equilibrium must accept per-time heat-of-reaction values.""" + kin = 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=[-5e3], + tref_hrxn=298.15, + ) + + temp = np.array([290.0, 310.0, 330.0]) + conc = np.array([[1.0, 1.0, 1.0, 1.0], + [2.0, 2.0, 2.0, 2.0], + [0.5, 0.5, 0.5, 0.5]]) + deltah = np.array([[-5.0e3], [-4.5e3], [-4.0e3]]) + + vectorized = np.asarray(kin.equilibrium_model(conc, temp, deltah)) + per_time = np.array([ + np.asarray( + kin.equilibrium_model(conc[i], float(temp[i]), deltah[i]) + ).ravel() + for i in range(temp.size) + ]) + + np.testing.assert_allclose(vectorized.reshape(per_time.shape), per_time)