High: [Kinetics] Custom rate-law setup and vectorized equilibrium reverse term crash

[sergey-gusev94]

Maintainer triage (June 29, 2026)

• Status: Verified bug against the current master source in this fork.
• Severity: High.
• Area labels: area:kinetics.
• Tackle batch: P3 - documented-mode crashes and tooling breakages.

The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.

---

Original audit details

Severity: High · Category: crash + indexing · Part of #67

Two Kinetics.py defects, grouped as one kinetics-correctness fix.

---

1. (High) RxnKinetics with a user-defined kinetic_model crashes at construction

Location: PharmaPy/Kinetics.py:337-363

In the dict branch of set_params, the locals params_f/orders are assigned only inside if self.elem_flag:. For a custom kinetic_model, elem_flag=False, so the block is skipped and execution hits elif params_f is None: (359) referencing an unbound local → UnboundLocalError. Even guarded, line 363 self.params_f = orders references undefined orders. The constructor sets elem_flag=False whenever kinetic_model is not None (262-272) and unconditionally calls set_params(params_dict) (288).

• Trigger: the documented custom-rate-law workflow, RxnKinetics(..., kinetic_model=my_func, params_f=array, df_dstates=..., df_dtheta=...).
• Wrong result: constructor raises UnboundLocalError immediately; the entire user-defined rate-law feature is unusable.
• Secondary: the array branch (366-375) also never updates self.params_f when elem_flag=False, freezing custom f-parameters during parameter estimation.
• Fix: on the elem_flag=False path set params_f = params.get('params_f', None) before the elif, then self.params_f = np.asarray(params_f) (raise the intended RuntimeError when None); handle elem_flag=False in the array branch too.

2. (High) Vectorized equilibrium reverse term divides by Keq indexed on the time axis, not the reaction axis

Location: PharmaPy/Kinetics.py:519-531 (with equil_term 453-463)

The 2-D branch builds r_term[:, ind] = np.prod(conc**orders[ind], axis=1) / keq_temp[ind] for ind in range(num_rxns) (524-527). For array temp, equil_term returns shape (n_temp, n_rxns) (via np.outer(inv_temp, deltah_temp/gas_ct), 460), so keq_temp[ind] indexes the first (time) axis — the per-reaction Keq vector of row ind — not the scalar Keq of reaction ind. Worse, when deltah_rxn is 2-D (as getHeatOfRxn returns for array temp), np.outer flattens it to (n_temp, n_temp*n_rxns), so keq_temp is malformed and the divide can raise a broadcast error. The 1-D ODE-integration branch (519-522) is fine.

• Trigger: a reversible reaction (keq_params not None) in a CSTR/PFR/Batch reactor, then retrieve results / build heat profiles (2-D mole_conc + 1-D temp into the 2-D branch, e.g. CSTR.energy_balances with heat_prof=True).
• Wrong result: reverse (and net) rates for reversible reactions use the wrong Keq per time point → incorrect rate profiles and heat-of-reaction source terms; with 2-D deltah_rxn it can raise ValueError: operands could not be broadcast together.
• Fix: index the reaction axis (keq_temp[:, ind]) in the 2-D branch; and fix equil_term to avoid np.outer flattening when deltah_rxn is already 2-D (broadcast inv_temp[:,None]*deltah/gas_ct).

---

Related (not duplicates): #36 (reverse-term exponents vs custom orders) and #23 (analytical Jacobian ignores the reverse term) are different defects in the same function; #44 targets integer-dtype truncation on the elem_flag=True path.

[parkyr]

Triage: reproduced ✅ (both defects)

Both defects in the issue audit reproduce against current master. The predicted symptoms are exact — Defect 1 raises UnboundLocalError at Kinetics.py:359 (custom kinetic_model skips the if self.elem_flag: block, so the elif params_f is None: reads an unbound local), and Defect 2's array-temp branch indexes keq_temp on the time axis at Kinetics.py:527.

Environment

• Base commit: 1a3ab3155216f07e9a18c60687a7b6b1c5959202 (origin/master)
• Test command: .venv/bin/python -m pytest tests/test_kinetics_custom_and_equilibrium.py -v
• PharmaPy.Kinetics is pure NumPy — no assimulo/solver backend needed.

Observed vs. expected

#	Trigger	Observed on base	Expected
1	RxnKinetics(..., kinetic_model=fn, params_f=arr)	UnboundLocalError: ... 'params_f' ... at Kinetics.py:359	constructs; self.params_f holds the passed f-parameters
2	equilibrium_model(2-D conc, array temp), reversible rxn	net rate reuses row-0 Keq: [0.5276, 3.0551, 0.0138]	per-time Keq: [0.5276, 2.9198, -0.0537]

Failing regression test

Save as tests/test_kinetics_custom_and_equilibrium.py:

# -*- 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

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 test_custom_kinetic_model_constructs(data_path):
    """Defect 1: the documented custom-rate-law workflow must construct."""

    def my_rate(conc, params_f):
        conc = np.maximum(1e-15, conc)
        return np.exp(np.dot(np.log(conc), np.atleast_2d(params_f).T))

    # 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=my_rate,
        params_f=np.array([[1.0, 1.0]]),
    )

    assert kin.elem_flag is False
    assert kin.kinetic_model is my_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_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)

Red output on unchanged base:

tests/test_kinetics_custom_and_equilibrium.py::test_custom_kinetic_model_constructs FAILED
tests/test_kinetics_custom_and_equilibrium.py::test_vectorized_equilibrium_matches_scalar_path FAILED

E   UnboundLocalError: cannot access local variable 'params_f' where it is not
    associated with a value                             (Kinetics.py:359)

E   AssertionError: Not equal to tolerance rtol=1e-07, atol=0
E   Mismatched elements: 2 / 3 (66.7%)
E    x: array([[0.527555], [3.05511 ], [0.013778]])     # vectorized
E    y: array([[0.527555], [2.919845], [-0.05373 ]])    # trusted scalar path

The second test uses the trusted scalar/1-D branch as its oracle, so the expected value is independent of the buggy vectorized expression. Confirmed non-vacuous: applying the minimal fix (Defect 1 — assign params_f/orders on the elem_flag=False path; Defect 2 — keq_temp[:, ind]) turns both green.

Acceptance criteria

1. RxnKinetics(..., kinetic_model=fn, params_f=arr) constructs without error and stores the passed f-parameters in self.params_f (and the array branch of set_params updates them during parameter estimation when elem_flag=False).
2. A custom-model construction with params_f=None still raises the intended RuntimeError.
3. The vectorized (array-temp, 2-D conc) equilibrium_model result equals the scalar/1-D path evaluated per time point (per-reaction Keq indexing).
4. equil_term handles a 2-D deltah_rxn (as getHeatOfRxn returns for array temp) without np.outer flattening or a broadcast error.
5. Existing tests continue to pass.

Next: gh-issue-to-pr 71.