Merge branch 'main' into add-multiexperiment

Author: smondal13

pyomo/contrib/doe/doe.py

@@ -48,9 +48,10 @@
 
 import pyomo.environ as pyo
 from pyomo.contrib.doe.utils import (
+    _SMALL_TOLERANCE_DEFINITENESS,
     check_FIM,
     compute_FIM_metrics,
-    _SMALL_TOLERANCE_DEFINITENESS,
+    regularize_fim_for_cholesky,
 )
 from pyomo.contrib.parmest.utils.model_utils import update_model_from_suffix
 
@@ -380,7 +381,6 @@ def run_doe(self, model=None, results_file=None):
         results_file: string name of the file path to save the results
                       to in the form of a .json file
                       default: None --> don't save
-
         """
         # Check results file name
         if results_file is not None:
@@ -461,11 +461,31 @@ def run_doe(self, model=None, results_file=None):
         model.obj_cons.activate()
 
         if self.use_grey_box:
-            self._initialize_grey_box_block(
-                model.obj_cons.egb_fim_block,
-                np.asarray(self.get_FIM(model=model), dtype=np.float64),
-                model.parameter_names,
-            )
+            # Initialize grey box inputs to be fim values currently
+            for i in model.parameter_names:
+                for j in model.parameter_names:
+                    if list(model.parameter_names).index(i) >= list(
+                        model.parameter_names
+                    ).index(j):
+                        model.obj_cons.egb_fim_block.inputs[(j, i)].set_value(
+                            pyo.value(model.fim[(i, j)])
+                        )
+            # Set objective value
+            if self.objective_option == ObjectiveLib.trace:
+                trace_val = np.trace(np.linalg.pinv(self.get_FIM()))
+                model.obj_cons.egb_fim_block.outputs["A-opt"].set_value(trace_val)
+            elif self.objective_option == ObjectiveLib.determinant:
+                det_val = np.linalg.det(np.array(self.get_FIM()))
+                model.obj_cons.egb_fim_block.outputs["log-D-opt"].set_value(
+                    np.log(det_val)
+                )
+            elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
+                eig, _ = np.linalg.eig(np.array(self.get_FIM()))
+                model.obj_cons.egb_fim_block.outputs["E-opt"].set_value(np.min(eig))
+            elif self.objective_option == ObjectiveLib.condition_number:
+                eig, _ = np.linalg.eig(np.array(self.get_FIM()))
+                cond_number = np.log(np.abs(np.max(eig) / np.min(eig)))
+                model.obj_cons.egb_fim_block.outputs["ME-opt"].set_value(cond_number)
 
         # If the model has L, initialize it with the solved FIM
         if hasattr(model, "L"):
@@ -487,7 +507,7 @@ def run_doe(self, model=None, results_file=None):
             # Check if the FIM is positive definite
             # If not, add jitter to the diagonal
             # to ensure positive definiteness
-            min_eig = np.min(np.real(np.linalg.eigvals(fim_np)))
+            min_eig = np.min(np.linalg.eigvals(fim_np))
 
             if min_eig < _SMALL_TOLERANCE_DEFINITENESS:
                 # Raise the minimum eigenvalue to at
@@ -619,6 +639,92 @@ def run_doe(self, model=None, results_file=None):
             with open(results_file, "w") as file:
                 json.dump(self.results, file)
 
+    def _get_fim_numpy(self, model):
+        """
+        Assemble the current FIM variable values into a NumPy array.
+
+        Parameters
+        ----------
+        model: ConcreteModel
+            DoE model containing variable ``fim``.
+
+        Returns
+        -------
+        ndarray
+            Dense FIM array. If ``only_compute_fim_lower`` is True, the
+            returned array is symmetrized from the lower triangle.
+        """
+        fim_vals = [
+            pyo.value(model.fim[i, j])
+            for i in model.parameter_names
+            for j in model.parameter_names
+        ]
+        fim_np = np.array(fim_vals, dtype=float).reshape(
+            (len(model.parameter_names), len(model.parameter_names))
+        )
+        if self.only_compute_fim_lower:
+            fim_np = fim_np + fim_np.T - np.diag(np.diag(fim_np))
+        return fim_np
+
+    def _initialize_cholesky_from_fim(self, model=None):
+        """
+        Synchronize Cholesky-related variables using the current FIM.
+
+        Parameters
+        ----------
+        model: ConcreteModel, optional
+            DoE model to update. Defaults to ``self.model``.
+
+        Returns
+        -------
+        None
+            Updates model values in place for available variables:
+            ``L``, ``L_inv``, ``fim_inv``, and ``cov_trace``.
+        """
+        if model is None:
+            model = self.model
+        if not hasattr(model, "L"):
+            # The model doesn't have the Cholesky variables, so we can't initialize them.
+            # This happens if the function is called with a model using GreyBox.
+            return
+
+        fim_np = self._get_fim_numpy(model)
+        fim_pd, _ = regularize_fim_for_cholesky(fim_np)
+
+        L_vals = np.linalg.cholesky(fim_pd)
+        for i, c in enumerate(model.parameter_names):
+            for j, d in enumerate(model.parameter_names):
+                if i >= j:
+                    model.L[c, d].value = L_vals[i, j]
+                else:
+                    model.L[c, d].value = 0.0
+
+        if hasattr(model, "L_inv"):
+            L_inv_vals = np.linalg.inv(L_vals)
+            for i, c in enumerate(model.parameter_names):
+                for j, d in enumerate(model.parameter_names):
+                    if i >= j:
+                        model.L_inv[c, d].value = L_inv_vals[i, j]
+                    else:
+                        model.L_inv[c, d].value = 0.0
+
+        if hasattr(model, "fim_inv"):
+            # Use the pseudo-inverse here rather than the strict inverse.
+            # The jittered matrix should be positive definite, but ``pinv``
+            # is safer for borderline ill-conditioned cases and matches the
+            # defensive approach already used when initializing ``fim_inv``
+            # from user-provided starting values.
+            fim_inv_vals = np.linalg.pinv(fim_pd)
+            for i, c in enumerate(model.parameter_names):
+                for j, d in enumerate(model.parameter_names):
+                    if self.only_compute_fim_lower and i < j:
+                        model.fim_inv[c, d].value = 0.0
+                    else:
+                        model.fim_inv[c, d].value = fim_inv_vals[i, j]
+
+        if hasattr(model, "cov_trace"):
+            model.cov_trace.value = np.trace(fim_inv_vals)
+
     def optimize_experiments(
         self,
         results_file=None,
@@ -2297,6 +2403,7 @@ def create_doe_model(
             self.only_compute_fim_lower
             and self.objective_option == ObjectiveLib.determinant
             and not self.Cholesky_option
+            and not self.use_grey_box
         ):
             raise ValueError(
                 "Cannot compute determinant with explicit formula "
@@ -2811,30 +2918,14 @@ def create_objective_function(self, model=None):
         model.obj_cons = pyo.Block()
 
         # Assemble the FIM matrix. This is helpful for initialization!
-        # collect current FIM values in row-major order
-        fim_vals = [
-            model.fim[bu, un].value
-            for bu in model.parameter_names
-            for un in model.parameter_names
-        ]
-        fim = np.array(fim_vals).reshape(
-            len(model.parameter_names), len(model.parameter_names)
-        )
+        fim = self._get_fim_numpy(model)
 
         ### Initialize the Cholesky decomposition matrix
         if self.Cholesky_option and self.objective_option in (
             ObjectiveLib.determinant,
             ObjectiveLib.trace,
         ):
-            # Calculate the eigenvalues of the FIM matrix
-            eig = np.linalg.eigvals(fim)
-
-            # If the smallest eigenvalue is (practically) negative,
-            # add a diagonal matrix to make it positive definite
-            if min(eig) < small_number:
-                fim = fim + np.eye(len(model.parameter_names)) * (
-                    small_number - min(eig)
-                )
+            fim, _ = regularize_fim_for_cholesky(fim)
 
             # Compute the Cholesky decomposition of the FIM matrix
             L = np.linalg.cholesky(fim)
@@ -3466,6 +3557,20 @@ def FIM_egb_cons(m, p1, p2):
         )
 
         if build_objective:
+            if self.objective_option == ObjectiveLib.trace:
+                output_name = "A-opt"
+            elif self.objective_option == ObjectiveLib.pseudo_trace:
+                output_name = "pseudo-A-opt"
+            elif self.objective_option == ObjectiveLib.determinant:
+                output_name = "log-D-opt"
+            elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
+                output_name = "E-opt"
+            elif self.objective_option == ObjectiveLib.condition_number:
+                output_name = "ME-opt"
+            else:
+                # Error path is intentionally deferred to the external model.
+                output_name = "A-opt"
+
             model.objective = pyo.Objective(
                 expr=model.obj_cons.egb_fim_block.outputs[output_name],
                 sense=(

pyomo/contrib/doe/grey_box_utilities.py

@@ -139,9 +139,11 @@ def __init__(
     def _get_FIM(self):
         # Grabs the current FIM subject
         # to the input values.
-        # This function currently assumes
-        # that we use a lower triangular
-        # FIM.
+        # Inputs store one triangular half
+        # of a symmetric FIM. Reconstruct
+        # the full symmetric matrix here,
+        # consistent with manuscript equation S5.
+        # https://arxiv.org/abs/2604.03354v1
         upt_FIM = self._input_values
 
         # Create FIM in the correct way
@@ -195,6 +197,8 @@ def output_names(self):
 
         if self.objective_option == ObjectiveLib.trace:
             obj_name = "A-opt"
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            obj_name = "pseudo-A-opt"
         elif self.objective_option == ObjectiveLib.determinant:
             obj_name = "log-D-opt"
         elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
@@ -227,6 +231,8 @@ def evaluate_outputs(self):
 
         if self.objective_option == ObjectiveLib.trace:
             obj_value = np.trace(np.linalg.pinv(M))
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            obj_value = np.trace(M)
         elif self.objective_option == ObjectiveLib.determinant:
             sign, logdet = np.linalg.slogdet(M)
             obj_value = logdet
@@ -262,6 +268,8 @@ def finalize_block_construction(self, pyomo_block):
         # objective function.
         if self.objective_option == ObjectiveLib.trace:
             pyomo_block.outputs["A-opt"] = output_value
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            pyomo_block.outputs["pseudo-A-opt"] = output_value
         elif self.objective_option == ObjectiveLib.determinant:
             pyomo_block.outputs["log-D-opt"] = output_value
         elif self.objective_option == ObjectiveLib.minimum_eigenvalue:
@@ -299,6 +307,8 @@ def evaluate_jacobian_outputs(self):
             # is -inv(FIM) @ inv(FIM). Add reference to
             # pyomo.DoE 2.0 manuscript S.I.
             jac_M = -Minv @ Minv
+        elif self.objective_option == ObjectiveLib.pseudo_trace:
+            jac_M = np.eye(self._n_params, dtype=np.float64)
         elif self.objective_option == ObjectiveLib.determinant:
             Minv = np.linalg.pinv(M)
             # Derivative formula derived using tensor

pyomo/contrib/doe/tests/test_doe_errors.py

@@ -46,9 +46,6 @@
         RooneyBieglerMultiExperiment,
         RooneyBieglerMultiInputExperimentFlag,
     )
-    from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import (
-        RooneyBieglerExperiment,
-    )
 
 from pyomo.contrib.doe.examples.rooney_biegler_doe_example import run_rooney_biegler_doe
 import pyomo.environ as pyo
@@ -122,7 +119,7 @@ def test_experiment_none_error(self):
             # Experiment provided as None
             DoE_args = get_standard_args(None, fd_method, obj_used, flag_val)
 
-            doe_obj = DesignOfExperiments(**DoE_args)
+            DesignOfExperiments(**DoE_args)
 
     def test_experiment_empty_list_error(self):
         with self.assertRaisesRegex(
@@ -147,7 +144,7 @@ def test_reactor_check_no_get_labeled_model(self):
         ):
             DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val)
 
-            doe_obj = DesignOfExperiments(**DoE_args)
+            DesignOfExperiments(**DoE_args)
 
     def test_reactor_check_no_experiment_outputs(self):
         fd_method = "central"
@@ -1270,6 +1267,25 @@ def test_lhs_missing_bounds_error_message(self):
                 init_method="latin_hypercube_sampling", init_n_samples=2
             )
 
+    def test_invalid_determinant_without_cholesky(self):
+        fd_method = "central"
+        obj_used = "determinant"
+
+        experiment = get_rooney_biegler_experiment_flag()
+
+        DoE_args = get_standard_args(experiment, fd_method, obj_used, flag=None)
+        DoE_args["_Cholesky_option"] = False
+
+        doe_obj = DesignOfExperiments(**DoE_args)
+
+        # The explicit determinant formulation needs the full FIM, so we keep
+        # this regression test focused on the early validation guard.
+        with self.assertRaisesRegex(
+            ValueError,
+            "Cannot compute determinant with explicit formula if only_compute_fim_lower is True.",
+        ):
+            doe_obj.create_doe_model()
+
 
 if __name__ == "__main__":
     unittest.main()