Fix kinetics custom params and equilibrium indexing

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
 

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)