Skip to content

pyomo.doe adding more verbose output for sensitivity analysis #3525

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 30 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
2f35d58
added eigenvalues and determinant of the FIM
smondal13 Mar 18, 2025
0327879
added trace as a output in compute_FIM_full_factorial()
smondal13 Mar 18, 2025
eb78589
added eigvalsh for calculating
smondal13 Mar 18, 2025
5b2e90e
changed the method to calculate the minimum eigenvalues
smondal13 Mar 19, 2025
f3b0f7c
Merge branch 'Pyomo:main' into adding_eigen_values
smondal13 Mar 20, 2025
701947b
Added img_thereshold instead of the embedded constant for imaginary v…
smondal13 Mar 26, 2025
114a86e
Alex suggested some objectives and tests
smondal13 Mar 26, 2025
cb694be
Added `IMG_THRESHOLD` instead of embedded variable
smondal13 Mar 31, 2025
0a938ac
Added FIM test metrics for doe.py
smondal13 Apr 4, 2025
acf06f0
Changed the ``isclose() `` method to ``==``
smondal13 Apr 16, 2025
41a6272
Changed back to ``isclose()``
smondal13 Apr 16, 2025
b55786f
Added todo list
smondal13 Apr 17, 2025
46f75c1
Merge branch 'Pyomo:main' into adding_eigen_values
smondal13 Apr 17, 2025
3cc77ef
imported `_SMALL_TOLERANCE_IMG` from `doe.py` and used it in the `com…
smondal13 Apr 17, 2025
3ca1c4e
Deleted comment "Alex said..."
smondal13 Apr 24, 2025
7c2535f
added 1e-6 as the _SMALL_TOLERANCE_IMG
smondal13 Apr 24, 2025
5216672
Deleted the figures that was generated by runnig the example
smondal13 Apr 24, 2025
708a21d
adding unittest
smondal13 Apr 28, 2025
39f26e6
Merge branch 'Pyomo:main' into adding_eigen_values
smondal13 Apr 29, 2025
e9e83b6
to test doe , added the "doe_test_example.py" script, which returns a…
smondal13 Apr 30, 2025
7064966
The test_example file is not required. It is deleted. in `test_doe_FI…
smondal13 Apr 30, 2025
a475ada
Merge branch 'Pyomo:main' into adding_eigen_values
smondal13 Apr 30, 2025
5011800
Coding with Alex about `compute_FIM_metrics` funciton.
smondal13 Apr 30, 2025
b9ea35b
Merge branch 'adding_eigen_values' of github.com:smondal13/pyomo into…
smondal13 Apr 30, 2025
6d26fdf
Changed the test file to test the new function and added logging for …
smondal13 Apr 30, 2025
32079fe
added working test code for _compute_FIM_metric function and _check_F…
smondal13 May 1, 2025
e1842cf
Splitted the test files into 2 separate files. Commented out the code…
smondal13 May 1, 2025
898e87f
Deleted TODO: "make static method", and the chunk of code in `compute…
smondal13 May 2, 2025
ec2a361
Added `Returns` in the `compute_FIM_full_factorial()` docstring. Thou…
smondal13 May 2, 2025
f0da48e
Merge branch 'main' into adding_eigen_values
smondal13 May 2, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -34,3 +34,9 @@ cplex.log

# Mac tracking files
*.DS_Store*

# DOE example .png files
example_reactor_compute_FIM_E_opt.png
example_reactor_compute_FIM_A_opt.png
example_reactor_compute_FIM_D_opt.png
example_reactor_compute_FIM_ME_opt.png
177 changes: 152 additions & 25 deletions pyomo/contrib/doe/doe.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -1374,6 +1381,30 @@ def check_model_FIM(self, model=None, FIM=None):
)
)

# Check FIM is positive definite and symmetric
self._check_FIM(FIM)

self.logger.info(
"FIM provided matches expected dimensions from model and is approximately positive (semi) definite."
)

@staticmethod
def _check_FIM(FIM):
"""Private method for basic diagonists on FIM to ensure that the FIM is square, positive definite and symmetric.

Parameters
----------
FIM: 2D numpy array representing the FIM

Returns
-------
None, but will raise error messages as needed

"""
# Ensure that the FIM is a square matrix
if FIM.shape[0] != FIM.shape[1]:
raise ValueError("FIM must be a square matrix")

# Compute the eigenvalues of the FIM
evals = np.linalg.eigvals(FIM)

Expand All @@ -1393,10 +1424,6 @@ def check_model_FIM(self, model=None, FIM=None):
)
)

self.logger.info(
"FIM provided matches expected dimensions from model and is approximately positive (semi) definite."
)

# Check the jacobian shape against what is expected from the model.
def check_model_jac(self, jac=None):
if jac.shape != (self.n_experiment_outputs, self.n_parameters):
Expand Down Expand Up @@ -1446,7 +1473,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(
Expand All @@ -1470,7 +1497,22 @@ def compute_FIM_full_factorial(
method: string to specify which method should be used
options are ``kaug`` and ``sequential``

Returns
-------
fim_factorial_results: a dictionary of the results with the following keys and
their corresponding values as a list. Each element in the list corresponds to
a different design point in the full factorial space.
"log10 D-opt": list of log10(D-optimality)
"log10 A-opt": list of log10(A-optimality)
"log10 E-opt": list of log10(E-optimality)
"log10 ME-opt": list of log10(ME-optimality)
"eigval_min": list of minimum eigenvalues
"eigval_max": list of maximum eigenvalues
"det_FIM": list of determinants
"trace_FIM": list of traces
"solve_time": list of solve times
"""

# Start timer
sp_timer = TicTocTimer()
sp_timer.tic(msg=None)
Expand Down Expand Up @@ -1505,15 +1547,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": [],
}
)
Expand Down Expand Up @@ -1575,24 +1621,9 @@ 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
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:
self.logger.warning(
"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))
det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
_compute_FIM_metrics(FIM)
)

# Append the values for each of the experiment inputs
for k, v in model.experiment_inputs.items():
Expand All @@ -1602,6 +1633,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.min())
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
Expand Down Expand Up @@ -2314,3 +2349,95 @@ def _sgn(self, p):
return 1
else:
return -1


# loggers for the functions
function_logger = logging.getLogger(__name__)
function_logger.setLevel(logging.WARNING)


# Functions to compute FIM metrics
def _compute_FIM_metrics(FIM):
"""
This private function calculates the FIM metrics and returns them as a tuple.

Parameters
----------
FIM: 2D numpy array of the FIM

Returns
-------
det_FIM: determinant of the FIM
trace_FIM: trace of the FIM
E_vals: eigenvalues of the FIM
E_vecs: eigenvectors of the FIM
D_opt: log10(D-optimality) metric
A_opt: log10(A-optimality) metric
E_opt: log10(E-optimality) metric
ME_opt: log10(Modified E-optimality) metric
"""

# Check whether the FIM is square, positive definite, and symmetric
DesignOfExperiments._check_FIM(FIM)

# Compute FIM metrics
det_FIM = np.linalg.det(FIM)
D_opt = np.log10(det_FIM)

trace_FIM = np.trace(FIM)
A_opt = np.log10(trace_FIM)

E_vals, E_vecs = np.linalg.eig(FIM)
E_ind = np.argmin(E_vals.real) # index of smallest eigenvalue

# Warn the user if there is a ``large`` imaginary component (should not be)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will you ever get to this warning? Or will _check_FIM throw an error first?

Copy link
Author

@smondal13 smondal13 May 2, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@adowling2 _check_FIM() will not trigger this warning. _check_FIM() only checks whether the FIM is square, symmetric, and positive definite. Although this method checks for the symmetry of the FIM with a tolerance, and it is unlikely that the FIM will have large imaginary eigenvalues, we have kept the warning there just in case there is a large imaginary eigenvalue.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@smondal13 Thank you.

@blnicho @mrmundt Is there a best practice for testing this warning message is properly thrown?

if abs(E_vals.imag[E_ind]) > _SMALL_TOLERANCE_IMG:
function_logger.warning(
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
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, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt


# Standalone Function for user to calculate FIM metrics directly without using the class
def get_FIM_metrics(FIM):
"""
This function calculates the FIM metrics and returns them as a dictionary.

Parameters
----------
FIM: 2D numpy array of the FIM

Returns
-------
det_FIM: determinant of the FIM
trace_FIM: trace of the FIM
E_vals: eigenvalues of the FIM
E_vecs: eigenvectors of the FIM
D_opt: D-optimality metric
A_opt: A-optimality metric
E_opt: E-optimality metric
ME_opt: Modified E-optimality metric
"""
det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
_compute_FIM_metrics(FIM)
)

return {
"Determinanat of FIM": det_FIM,
"Trace of FIM": trace_FIM,
"Eigenvalues": E_vals,
"Eigen vectors": E_vecs,
"log10(D-Optimality)": D_opt,
"log10(A-Optimality)": A_opt,
"log10(E-Optimality)": E_opt,
"log10(Modified E-Optimality)": ME_opt,
}
73 changes: 73 additions & 0 deletions pyomo/contrib/doe/tests/test_doe_FIM_metrics.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2025
# National Technology and Engineering Solutions of Sandia, LLC
# Under the terms of Contract DE-NA0003525 with National Technology and
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
# rights in this software.
# This software is distributed under the 3-clause BSD License.
# ___________________________________________________________________________
from pyomo.common.dependencies import (
numpy as np,
numpy_available,
)

import pyomo.common.unittest as unittest
from pyomo.contrib.doe.doe import (
_SMALL_TOLERANCE_IMG,
_compute_FIM_metrics,
)


@unittest.skipIf(not numpy_available, "Numpy is not available")
class TestComputeFIMMetrics(unittest.TestCase):
def test_compute_FIM_metrics(self):
# Create a sample Fisher Information Matrix (FIM)
FIM = np.array([[10, 2], [2, 3]])

det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
_compute_FIM_metrics(FIM)
)

# expected results
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, E_vecs_expected = np.linalg.eig(FIM)
min_eigval = np.min(E_vals_expected.real)
if min_eigval <= 0:
E_opt_expected = np.nan
else:
E_opt_expected = np.log10(min_eigval)

cond_expected = np.linalg.cond(FIM)

ME_opt_expected = np.log10(cond_expected)

# Test results
self.assertEqual(det_FIM, det_expected)
self.assertEqual(trace_FIM, trace_expected)
self.assertTrue(np.allclose(E_vals, E_vals_expected))
self.assertTrue(np.allclose(E_vecs, E_vecs_expected))
self.assertEqual(D_opt, D_opt_expected)
self.assertEqual(A_opt, A_opt_expected)
self.assertEqual(E_opt, E_opt_expected)
self.assertEqual(ME_opt, ME_opt_expected)


class TestFIMWarning(unittest.TestCase):
def test_FIM_eigenvalue_warning(self):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@blnicho @mrmundt I have not tested warning messages before, so your input here is appreciated.

# Create a matrix with an imaginary component large enough to trigger the warning
FIM = np.array([[6, 5j], [5j, 7]])
with self.assertLogs("pyomo.contrib.doe", level="WARNING") as cm:
_compute_FIM_metrics(FIM)
expected_warning = f"Eigenvalue has imaginary component greater than {_SMALL_TOLERANCE_IMG}, contact developers if this issue persists."
self.assertIn(expected_warning, cm.output[0])


if __name__ == "__main__":
unittest.main()
69 changes: 69 additions & 0 deletions pyomo/contrib/doe/tests/test_doe_check_FIM.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2025
# National Technology and Engineering Solutions of Sandia, LLC
# Under the terms of Contract DE-NA0003525 with National Technology and
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
# rights in this software.
# This software is distributed under the 3-clause BSD License.
# ___________________________________________________________________________
from pyomo.common.dependencies import (
numpy as np,
numpy_available,
)

import pyomo.common.unittest as unittest
from pyomo.contrib.doe import DesignOfExperiments

# Not need? from pyomo.contrib.doe.tests import doe_test_example
from pyomo.contrib.doe.doe import (
_SMALL_TOLERANCE_DEFINITENESS,
_SMALL_TOLERANCE_SYMMETRY,
)


@unittest.skipIf(not numpy_available, "Numpy is not available")
class TestDesignOfExperimentsCheckFIM(unittest.TestCase):
"""Test the check_FIM method of the DesignOfExperiments class."""

def test_check_FIM_valid(self):
"""Test case where the FIM is valid (square, positive definite, symmetric)."""
FIM = np.array([[4, 1], [1, 3]])
try:
# Call the static method directly
DesignOfExperiments._check_FIM(FIM)
except ValueError as e:
self.fail(f"Unexpected error: {e}")

def test_check_FIM_non_square(self):
"""Test case where the FIM is not square."""
FIM = np.array([[4, 1], [1, 3], [2, 1]])
with self.assertRaisesRegex(ValueError, "FIM must be a square matrix"):
DesignOfExperiments._check_FIM(FIM)

def test_check_FIM_non_positive_definite(self):
"""Test case where the FIM is not positive definite."""
FIM = np.array([[1, 0], [0, -2]])
with self.assertRaisesRegex(
ValueError,
r"FIM provided is not positive definite. It has one or more negative eigenvalue\(s\) less than -{:.1e}".format(
_SMALL_TOLERANCE_DEFINITENESS
),
):
DesignOfExperiments._check_FIM(FIM)

def test_check_FIM_non_symmetric(self):
"""Test case where the FIM is not symmetric."""
FIM = np.array([[4, 1], [0, 3]])
with self.assertRaisesRegex(
ValueError,
"FIM provided is not symmetric using absolute tolerance {}".format(
_SMALL_TOLERANCE_SYMMETRY
),
):
DesignOfExperiments._check_FIM(FIM)


if __name__ == "__main__":
unittest.main()