Added FIM test metrics for doe.py

Author: smondal13

pyomo/contrib/doe/doe.py

@@ -1523,7 +1523,7 @@ def compute_FIM_full_factorial(
             FIM = self._computed_FIM
 
             '''
-            # Alex said to makek this a function
+            # Alex said to make this a function
             # This should be a static (?) function
             # Making it a function allows us to perform error tests
             # on the FIM
@@ -1535,9 +1535,11 @@ def compute_metrics(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
             IMG_THERESHOLD = 1e-6  # Threshold for imaginary component
@@ -1564,9 +1566,9 @@ def compute_metrics(FIM):
             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[-1])
-            fim_factorial_results["det_FIM"].append(np.linalg.det(FIM))
-            fim_factorial_results["trace_FIM"].append(np.trace(FIM))
+            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,90 @@
+from pyomo.common.dependencies import (
+    numpy as np,
+    pandas as pd
+)
+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  # Threshold for imaginary component
+    # Warn the user if there is a ``large`` imaginary component (should not be)
+    if abs(E_vals.imag[E_ind]) > IMG_THERESHOLD:
+        print(
+            f"Eigenvalue has imaginary component greater than {IMG_THERESHOLD}, 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-06, contact developers if this issue persists."
+
+    # Ensure expected message is in the output
+    assert expected_message in out