High: [ParamEstim] Fixed-parameter gradients and IPOPT residual weighting are broken

[sergey-gusev94]

Maintainer triage (June 29, 2026)

• Status: Verified bug against the current master source in this fork.
• Severity: High.
• Area labels: area:param-estimation.
• 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 · Part of #67 (round-2 audit batch)

Both defects are in PharmaPy/ParamEstim.py within the single get_gradient method (and its helper func_aux). They share the gradient-assembly subsystem and the reconstruct_params / residual-ordering / weighting machinery, and would be fixed and tested together in one parameter-estimation-correctness PR. They are distinct mechanisms (finite-difference indexing with optimize_flags vs IPOPT residual ordering+weighting) but tightly co-located and both must be corrected for fixed-parameter / multi-state / weighted fits to work.

This issue bundles 2 defects that share one fix/PR.

1. (High) Numerical-Jacobian branch crashes / feeds wrong params when optimize_flags fixes any parameter

Location: PharmaPy/ParamEstim.py:544-547 (with func_aux at 455-458; helper jac_module.py:105-111)

In the default finite-difference branch, the optimizer-supplied variable-only vector params (length num_params = map_variable.sum()) is passed to numerical_jac_data with pick_x=pick_p=np.where(self.map_variable)[0], which holds GLOBAL indices into the full parameter vector. Inside numerical_jac_data, delx/dx are sized to the variable-only vector, and delx[i]=dx[i] uses i as a global index, so when any fixed parameter shifts a variable parameter's global index >= num_params, delx[i]/dx[i] go out of bounds -> IndexError. Separately, func_aux (455-458) calls self.function(params,...) WITHOUT reconstruct_params (unlike get_objective at line 466), so even without an index crash the model is evaluated on a truncated/misaligned parameter vector (e.g. a 2-arg model receives 1 value -> ValueError 'not enough values to unpack'), making objective and gradient inconsistent.

• Trigger: ParameterEstimation(func, param_seed, x, y, optimize_flags=[True, False]) (or any flags fixing >=1 parameter) with default jac_fun=None and default method='LM'; LM calls get_gradient on the first iteration. Reproduced: fix-last-param -> ValueError (model gets size-1 vector); fix-first-param -> IndexError 'index 2 is out of bounds for axis 0 with size 2'.
• Wrong result: Parameter estimation with any fixed parameter aborts immediately (IndexError/ValueError) and is unrunnable; in the non-crashing edge case the Jacobian is computed on a wrong parameter vector, giving wrong gradient/fitted parameters/covariance.
• Suggested fix: Reconstruct the full parameter set before differentiating and have func_aux operate on the full vector: pass self.reconstruct_params(params) (length num_params_total) as x with pick_x=pick_p, and make func_aux evaluate self.function on the same reconstructed full-length vector used in get_objective so global-index perturbations land correctly.
• Net-new (not a duplicate): Distinct from High: [Reactors] Sensitivity Jacobian writes species block into oversized non-isothermal state block #55 (Reactors.py) and High: [ParamEstim/Calibration] IPOPT result assembly crashes and PCR predict recenters on new data #78 (IPOPT update_self keyword / Calibration.predict). This is the LM/default numerical-Jacobian path interacting with optimize_flags/pick_p and func_aux's missing reconstruct_params.

2. (High) IPOPT gradient pairs a state-major weighted Jacobian with a time-major, unweighted residual vector

Location: PharmaPy/ParamEstim.py:587-589

In the non-out_array branch, gradient = jacobian.T.dot(res) with res = np.concatenate([a.T.ravel() for a in self.residuals]). Two errors: (1) ORDERING -- self.residuals (line 517) is a 2D array (sum_times, n_meas); iterating its rows flattens res TIME-major [s0t0,s1t0,s0t1,...], but jacobian (concat_sens, built via two reorder_sens calls at lines 562,565) is STATE-major [s0t0,s0t1,s0t2,s1t0,...], so jacobian.T.dot(res) pairs each sensitivity row with the wrong residual for any model with >1 measured state. (2) WEIGHTING -- jacobian is sigma_inv-weighted (line 564) but res uses the RAW residuals (line 517); the correct gradient of 1/2||W^{-1/2} r||^2 needs the metric inv(W) (two sigma_inv factors), but the code applies only one. Verified: 2-state identity-weight test gives wrong gradient; single-state W=[[3]] gives code/true ratio = sqrt(3).

• Trigger: paramest.optimize_fn(method='IPOPT') on a model with >=2 measured states (or any single-state model fit with a non-identity weight_matrix); cyipopt calls jac=self.get_gradient so out_array defaults to False and this branch runs every IPOPT iteration.
• Wrong result: The optimizer receives a gradient whose components are mispaired with residual entries (multi-state) and/or omit residual weighting (non-identity weight_matrix). IPOPT converges to wrong parameter estimates (and wrong covariance/CIs) or fails to converge; silently incorrect fit. The LM path is unaffected (it uses out_array=True).
• Suggested fix: Build res to match the jacobian's state-major order and apply the same weighting: reuse the weighted, state-major residual produced in get_objective (residual_out at line 526) by storing it on self, or recompute res = np.concatenate([(r @ self.sigma_inv).T.ravel() for r in self.resid_runs]) so both ordering and weighting match concat_sens.
• Net-new (not a duplicate): Distinct from High: [ParamEstim/Calibration] IPOPT result assembly crashes and PCR predict recenters on new data #78 (post-convergence get_objective keyword crash at 641-642, reached AFTER minimize_ipopt) and High: [Reactors] Sensitivity Jacobian writes species block into oversized non-isothermal state block #55 (reactor sensitivity, different file). This corrupts the gradient DURING the IPOPT solve in the out_array=False branch.

Related existing issues (referenced, not modified): #55 (Reactors.py sensitivity block, different file), #78 (IPOPT update_self keyword crash at 641-642 and PCR_calibration.predict) -- reference only; neither covers these two mechanisms.

---

Found by a multi-agent source audit (round 2); confirmed by re-reading the source and cross-checked against all 80 existing open issues. Each defect was verified by an independent code-truth skeptic and a novelty skeptic.