Fix IPOPT result assembly and PCR prediction centering

PharmaPy/Calibration.py

@@ -47,6 +47,25 @@ def __init__(self, data, num_comp=None, standardize=True, snv=False,
         self.y_suffixes = y_suffixes
 
     def __center_data(self, data=None, mean=None, std=None):
+        """Center and optionally standardize predictor data.
+
+        Parameters
+        ----------
+        data : numpy.ndarray
+            Predictor rows in the same units as the calibration data.
+        mean : numpy.ndarray, optional
+            Column means in predictor units. When omitted with ``std``, both
+            statistics are computed from ``data``.
+        std : numpy.ndarray, optional
+            Column standard deviations in predictor units.
+
+        Returns
+        -------
+        data_centered : numpy.ndarray
+            Centered predictors. Values are dimensionless [-] when
+            ``standardize`` or ``snv`` scaling is active.
+
+        """
         if mean is None and std is None:
             mean = data.mean(axis=0)
             std = data.std(axis=0)
@@ -189,12 +208,42 @@ def get_regression(self, y_data, num_comp=None, update_instance=True):
 
     def predict(self, inputs, num_comp=None, regression_coeff=None,
                 full_output=False):
+        """Predict responses for new predictor rows.
+
+        Parameters
+        ----------
+        inputs : array-like
+            New predictor data in the same units as the calibration data.
+        num_comp : int, optional
+            Number of principal components used for prediction.
+        regression_coeff : numpy.ndarray, optional
+            Regression coefficients with respect to principal-component
+            scores.
+        full_output : bool, optional
+            If True, return projections, predictions, and MSE diagnostics.
+
+        Returns
+        -------
+        response : numpy.ndarray or dict
+            Predicted responses when ``full_output`` is False. If
+            ``full_output`` is True, return a dictionary containing
+            projections, predictions, and MSE. Prediction units match the
+            responses used to fit ``regression_coeff``.
+
+        Notes
+        -----
+        For non-SNV models, new predictors are centered with the training
+        mean and standard deviation so the scores remain on the fitted
+        principal-component basis.
+
+        """
         inputs = np.atleast_2d(inputs)
 
         if self.snv:
             inputs_centered = self.__center_data(inputs)
         else:
-            inputs_centered = self.__center_data(inputs)
+            inputs_centered = self.__center_data(inputs, self.data_mean,
+                                                 self.data_std)
 
         if regression_coeff is None:
             coeff = self.regression_coeff

PharmaPy/ParamEstim.py

@@ -598,6 +598,36 @@ def get_cond_number(self, sens_matrix):
 
     def optimize_fn(self, optim_options=None, simulate=False, verbose=True,
                     store_iter=True, method='LM', bounds=None):
+        """Optimize variable parameters and assemble fit statistics.
+
+        Parameters
+        ----------
+        optim_options : dict, optional
+            Solver options passed to the selected optimization method.
+        simulate : bool, optional
+            If True, evaluate residuals in simulation mode by disabling the
+            optimization residual contribution.
+        verbose : bool, optional
+            If True, request verbose solver output when supported.
+        store_iter : bool, optional
+            If True, store unique parameter/objective iterates.
+        method : {'LM', 'IPOPT'}, optional
+            Optimization method used for fitting.
+        bounds : sequence, optional
+            Parameter bounds passed to IPOPT.
+
+        Returns
+        -------
+        opt_par : numpy.ndarray
+            Optimized variable parameters.
+        covar_params : numpy.ndarray
+            Estimated covariance matrix for the variable parameters.
+        info : dict
+            Solver Jacobian and residual information. For IPOPT,
+            ``info['fun']`` stores weighted residuals, dimensionless [-] after
+            applying ``sigma_inv``.
+
+        """
 
         self.optimize_flag = not simulate
         self.opt_method = method
@@ -638,8 +668,10 @@ def optimize_fn(self, optim_options=None, simulate=False, verbose=True,
             # final_sens = np.vstack(self.sens_runs)
             # final_fun = np.concatenate(self.resid_runs)
 
+            # Assemble final weighted residuals without overwriting the
+            # solved-state residuals stored during the IPOPT callbacks.
             resid_multidim = self.get_objective(opt_par, out_array=True,
-                                                update_self=False)
+                                                set_self=False)
             jac_multidim = self.get_gradient(opt_par, out_array=True)
             info = {'jac': jac_multidim, 'fun': resid_multidim}
 

tests/test_paramestim_calibration_fit_predict.py

@@ -0,0 +1,120 @@
+"""Regression tests for issue #78 fit and prediction paths.
+
+Dimensional fixture names carry units, and comments call out dimensionless
+intermediates as [-] where normalization or projection removes units.
+"""
+
+import numpy as np
+import pytest
+
+from PharmaPy import ParamEstim
+from PharmaPy.Calibration import PCR_calibration
+
+
+pytestmark = pytest.mark.unit
+
+
+def test_parameter_estimation_ipopt_result_assembly_uses_base_keyword(
+        monkeypatch):
+    """Exercise IPOPT post-solve assembly while stubbing only the solver.
+
+    The optional cyipopt/IPOPT dependency is replaced with a deterministic
+    boundary fake, but the ParameterEstimation objective, gradient, residual,
+    y_model, and covariance assembly paths remain real. Time is [s],
+    concentration is [mol/L], and the fitted rate is [mol/L/s].
+    """
+    time_s = np.array([0.0, 1.0, 2.0])
+    rate_seed_mol_l_s = 1.0
+    rate_mol_l_s = 2.0
+    y_model_mol_l = rate_mol_l_s * time_s
+    residual_offset_mol_l = np.array([0.10, -0.05, 0.20])
+    y_obs_mol_l = y_model_mol_l + residual_offset_mol_l
+
+    def linear_model(params, x_data_s):
+        """Return concentration [mol/L] from rate [mol/L/s] and time [s]."""
+        return params[0] * x_data_s
+
+    def linear_jacobian(params, x_data_s):
+        """Return d(concentration)/d(rate) sensitivities with units [s]."""
+        return x_data_s[np.newaxis, :]
+
+    estimator = ParamEstim.ParameterEstimation(
+        linear_model,
+        param_seed=[rate_seed_mol_l_s],
+        x_data=time_s,
+        y_data=y_obs_mol_l,
+        name_params=["rate_mol_l_s"],
+        jac_fun=linear_jacobian,
+    )
+
+    def fake_minimize_ipopt(objective, params_var, jac=None, bounds=None,
+                            options=None, kwargs=None):
+        """Mimic IPOPT returning the solved rate [mol/L/s]."""
+        optimum_mol_l_s = np.array([rate_mol_l_s])
+        # Match IPOPT's solved-state callback: residuals have units of the
+        # measured response before weighting, here [mol/L].
+        objective(optimum_mol_l_s, **(kwargs or {}))
+        return {"x": optimum_mol_l_s}
+
+    # cyipopt/IPOPT is an optional external solver stack absent from the core
+    # test lane. Patch only that boundary; objective, gradient, and covariance
+    # assembly stay on the real ParameterEstimation methods.
+    monkeypatch.setattr(ParamEstim, "have_cyipopt", True)
+    monkeypatch.setattr(ParamEstim, "minimize_ipopt", fake_minimize_ipopt,
+                        raising=False)
+
+    opt_par_mol_l_s, covar_rate, info = estimator.optimize_fn(
+        method="IPOPT", verbose=False)
+
+    # The default identity weight matrix leaves the [mol/L] residual and [s]
+    # sensitivity values numerically unchanged after sigma_inv weighting.
+    expected_weighted_residuals = -residual_offset_mol_l
+    expected_weighted_jacobian_s = time_s[np.newaxis, :]
+
+    np.testing.assert_allclose(opt_par_mol_l_s, [rate_mol_l_s])
+    np.testing.assert_allclose(info["fun"], expected_weighted_residuals)
+    np.testing.assert_allclose(info["jac"], expected_weighted_jacobian_s)
+    np.testing.assert_allclose(estimator.y_model[0].ravel(), y_model_mol_l)
+    # Covariance entries correspond to rate variance units [(mol/L/s)^2].
+    assert covar_rate.shape == (1, 1)
+
+
+def test_pcr_predict_uses_training_centering_for_single_new_spectrum():
+    """Predict a single absorbance spectrum [AU] with training statistics."""
+    # Calibration predictor rows are spectra [AU] at three wavelengths.
+    spectra_au = np.array([
+        [0.20, 1.10, 2.40],
+        [0.45, 1.35, 2.95],
+        [0.80, 1.85, 3.45],
+        [1.10, 2.10, 4.05],
+    ])
+    # Response concentrations are [g/L].
+    concentration_g_l = np.array([1.2, 1.8, 2.6, 3.1])
+
+    num_comp = 2  # number of retained principal components [-]
+    calibration = PCR_calibration(spectra_au, num_comp=num_comp,
+                                  standardize=True)
+    calibration.get_regression(concentration_g_l, num_comp=num_comp)
+
+    # Single prediction spectrum is in the same absorbance units [AU].
+    new_spectrum_au = np.array([[0.70, 1.70, 3.20]])
+    prediction_g_l = calibration.predict(new_spectrum_au)
+
+    # Predictor centering uses the training absorbance statistics [AU], and
+    # division by the training standard deviation makes the predictors [-].
+    training_mean_au = spectra_au.mean(axis=0)
+    training_std_au = spectra_au.std(axis=0)
+    centered_new = (
+        (new_spectrum_au - training_mean_au) / training_std_au
+    )
+    # SVD loadings and principal-component scores are dimensionless [-].
+    loadings = calibration.svd_dict["V"][:, :num_comp]
+    scores = centered_new @ loadings
+    # Regression coefficients convert dimensionless scores [-] to [g/L], and
+    # y_means is the response offset [g/L].
+    regression_coeff_g_l = calibration.regression_coeff
+    response_offset_g_l = calibration.y_means
+    expected_g_l = scores @ regression_coeff_g_l + response_offset_g_l
+
+    assert np.all(np.isfinite(prediction_g_l))
+    np.testing.assert_allclose(prediction_g_l, expected_g_l)