pyomo.doe adding more verbose output for sensitivity analysis

pyomo/contrib/doe/doe.py

@@ -61,6 +61,13 @@
 # below and the tests. The user should not need to adjust it.
 _SMALL_TOLERANCE_SYMMETRY = 1e-6
 
+# This small and positive tolerance is used to check
+# if the imaginary part of the eigenvalues of the FIM is
+# greater than a small tolerance. It is defined as a
+# tolerance here to ensure consistency between the code
+# below and the tests. The user should not need to adjust it.
+_SMALL_TOLERANCE_IMG = 1e-6
+
 
 class ObjectiveLib(Enum):
     determinant = "determinant"
@@ -1446,7 +1453,7 @@ def update_FIM_prior(self, model=None, FIM=None):
 
         self.logger.info("FIM prior has been updated.")
 
-    # ToDo: Add an update function for the parameter values? --> closed loop parameter estimation?
+    # TODO: Add an update function for the parameter values? --> closed loop parameter estimation?
     # Or leave this to the user?????
     def update_unknown_parameter_values(self, model=None, param_vals=None):
         raise NotImplementedError(
@@ -1505,15 +1512,19 @@ def compute_FIM_full_factorial(
                 "Design ranges keys must be a subset of experimental design names."
             )
 
-        # ToDo: Add more objective types? i.e., modified-E; G-opt; V-opt; etc?
-        # ToDo: Also, make this a result object, or more user friendly.
+        # TODO: Add more objective types? i.e., modified-E; G-opt; V-opt; etc?
+        # TODO: Also, make this a result object, or more user friendly.
         fim_factorial_results = {k.name: [] for k, v in model.experiment_inputs.items()}
         fim_factorial_results.update(
             {
                 "log10 D-opt": [],
                 "log10 A-opt": [],
                 "log10 E-opt": [],
                 "log10 ME-opt": [],
+                "eigval_min": [],
+                "eigval_max": [],
+                "det_FIM": [],
+                "trace_FIM": [],
                 "solve_time": [],
             }
         )
@@ -1576,14 +1587,18 @@ def compute_FIM_full_factorial(
             FIM = self._computed_FIM
 
             # Compute and record metrics on FIM
-            D_opt = np.log10(np.linalg.det(FIM))
-            A_opt = np.log10(np.trace(FIM))
-            E_vals, E_vecs = np.linalg.eig(FIM)  # Grab eigenvalues
+            det_FIM = np.linalg.det(FIM)  # Determinant of FIM
+            D_opt = np.log10(det_FIM)
+            trace_FIM = np.trace(FIM)  # Trace of FIM
+            A_opt = np.log10(trace_FIM)
+            E_vals, E_vecs = np.linalg.eig(FIM)  # Grab eigenvalues and eigenvectors
+
             E_ind = np.argmin(E_vals.real)  # Grab index of minima to check imaginary
+
             # Warn the user if there is a ``large`` imaginary component (should not be)
-            if abs(E_vals.imag[E_ind]) > 1e-8:
+            if abs(E_vals.imag[E_ind]) > _SMALL_TOLERANCE_IMG:
                 self.logger.warning(
-                    "Eigenvalue has imaginary component greater than 1e-6, contact developers if this issue persists."
+                    f"Eigenvalue has imaginary component greater than {_SMALL_TOLERANCE_IMG}, contact developers if this issue persists."
                 )
 
             # If the real value is less than or equal to zero, set the E_opt value to nan
@@ -1602,6 +1617,10 @@ def compute_FIM_full_factorial(
             fim_factorial_results["log10 A-opt"].append(A_opt)
             fim_factorial_results["log10 E-opt"].append(E_opt)
             fim_factorial_results["log10 ME-opt"].append(ME_opt)
+            fim_factorial_results["eigval_min"].append(E_vals[0])
+            fim_factorial_results["eigval_max"].append(E_vals.max())
+            fim_factorial_results["det_FIM"].append(det_FIM)
+            fim_factorial_results["trace_FIM"].append(trace_FIM)
             fim_factorial_results["solve_time"].append(time_set[-1])
 
         self.fim_factorial_results = fim_factorial_results

pyomo/contrib/doe/tests/test_doe_FIM_metrics.py

@@ -0,0 +1,89 @@
+import numpy as np
+
+# from pyomo.common.dependencies import numpy as np
+from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_IMG
+import pytest
+
+
+def compute_FIM_metrics(FIM):
+    # Compute and record metrics on FIM
+    det_FIM = np.linalg.det(FIM)  # Determinant of FIM
+    D_opt = np.log10(det_FIM)
+    trace_FIM = np.trace(FIM)  # Trace of FIM
+    A_opt = np.log10(trace_FIM)
+    E_vals, E_vecs = np.linalg.eig(FIM)  # Grab eigenvalues and eigenvectors
+
+    E_ind = np.argmin(E_vals.real)  # Grab index of minima to check imaginary
+    IMG_THERESHOLD = 1e-6  # Instead of creating a new constant, use `SMALL_DIFF` by importiing it form `doe.py`
+    # Warn the user if there is a ``large`` imaginary component (should not be)
+    if abs(E_vals.imag[E_ind]) > _SMALL_TOLERANCE_IMG:
+        print(
+            "Eigenvalue has imaginary component greater than 1e-6, contact developers if this issue persists."
+        )
+
+    # If the real value is less than or equal to zero, set the E_opt value to nan
+    if E_vals.real[E_ind] <= 0:
+        E_opt = np.nan
+    else:
+        E_opt = np.log10(E_vals.real[E_ind])
+
+    ME_opt = np.log10(np.linalg.cond(FIM))
+
+    return {
+        "det_FIM": det_FIM,
+        "trace_FIM": trace_FIM,
+        "E_vals": E_vals,
+        "D_opt": D_opt,
+        "A_opt": A_opt,
+        "E_opt": E_opt,
+        "ME_opt": ME_opt,
+    }
+
+
+def test_FIM_metrics():
+    # Create a sample Fisher Information Matrix (FIM)
+    FIM = np.array([[4, 2], [2, 3]])
+
+    # Call the function to compute metrics
+    results = compute_FIM_metrics(FIM)
+
+    # Use known values for assertions
+    det_expected = np.linalg.det(FIM)
+    D_opt_expected = np.log10(det_expected)
+
+    trace_expected = np.trace(FIM)
+    A_opt_expected = np.log10(trace_expected)
+
+    E_vals_expected, _ = np.linalg.eig(FIM)
+    min_eigval = np.min(E_vals_expected.real)
+
+    cond_expected = np.linalg.cond(FIM)
+
+    assert np.isclose(results["det_FIM"], det_expected)
+    assert np.isclose(results["trace_FIM"], trace_expected)
+    assert np.allclose(results["E_vals"], E_vals_expected)
+    assert np.isclose(results["D_opt"], D_opt_expected)
+    assert np.isclose(results["A_opt"], A_opt_expected)
+    if min_eigval.real > 0:
+        assert np.isclose(results["E_opt"], np.log10(min_eigval))
+    else:
+        assert np.isnan(results["E_opt"])
+
+    assert np.isclose(results["ME_opt"], np.log10(cond_expected))
+
+
+def test_FIM_metrics_warning_printed(capfd):
+    # Create a matrix with an imaginary component large enough to trigger the warning
+    FIM = np.array([[9, -2], [9, 3]])
+
+    # Call the function
+    compute_FIM_metrics(FIM)
+
+    # Capture stdout and stderr
+    out, err = capfd.readouterr()
+
+    # Correct expected message
+    expected_message = "Eigenvalue has imaginary component greater than 1e-6, contact developers if this issue persists."
+
+    # Ensure expected message is in the output
+    assert expected_message in out