diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 76c154b07ee..ab68ae19814 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -45,22 +45,14 @@ from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp import pyomo.environ as pyo +from pyomo.contrib.doe.utils import ( + check_FIM, + compute_FIM_metrics, + _SMALL_TOLERANCE_DEFINITENESS, +) from pyomo.opt import SolverStatus -# This small and positive tolerance is used when checking -# if the prior is negative definite or approximately -# indefinite. 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_DEFINITENESS = 1e-6 - -# This small and positive tolerance is used to check -# the FIM is approximately symmetric. 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_SYMMETRY = 1e-6 - class ObjectiveLib(Enum): determinant = "determinant" @@ -1383,24 +1375,8 @@ def check_model_FIM(self, model=None, FIM=None): ) ) - # Compute the eigenvalues of the FIM - evals = np.linalg.eigvals(FIM) - - # Check if the FIM is positive definite - if np.min(evals) < -_SMALL_TOLERANCE_DEFINITENESS: - raise ValueError( - "FIM provided is not positive definite. It has one or more negative eigenvalue(s) less than -{:.1e}".format( - _SMALL_TOLERANCE_DEFINITENESS - ) - ) - - # Check if the FIM is symmetric - if not np.allclose(FIM, FIM.T, atol=_SMALL_TOLERANCE_SYMMETRY): - raise ValueError( - "FIM provided is not symmetric using absolute tolerance {}".format( - _SMALL_TOLERANCE_SYMMETRY - ) - ) + # Check FIM is positive definite and symmetric + check_FIM(FIM) self.logger.info( "FIM provided matches expected dimensions from model and is approximately positive (semi) definite." @@ -1455,7 +1431,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( @@ -1474,12 +1450,37 @@ def compute_FIM_full_factorial( Parameters ---------- - model: model to perform the full factorial exploration on - design_ranges: dict of lists, of the form {: [start, stop, numsteps]} - method: string to specify which method should be used - options are ``kaug`` and ``sequential`` + model: DoE model, optional + model to perform the full factorial exploration on + design_ranges: dict + dictionary of lists, of the form {: [start, stop, numsteps]} + method: str, optional + to specify which method should be used. + Options are ``kaug`` and ``sequential`` + Returns + ------- + fim_factorial_results: dict + A dictionary of the results with the following keys and their corresponding + values as a list: + + - keys of model's experiment_inputs + - "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 + + Raises + ------ + ValueError + If the design_ranges' keys do not match the model's experiment_inputs' keys. """ + # Start timer sp_timer = TicTocTimer() sp_timer.tic(msg=None) @@ -1514,8 +1515,8 @@ 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( { @@ -1523,6 +1524,10 @@ def compute_FIM_full_factorial( "log10 A-opt": [], "log10 E-opt": [], "log10 ME-opt": [], + "eigval_min": [], + "eigval_max": [], + "det_FIM": [], + "trace_FIM": [], "solve_time": [], } ) @@ -1584,24 +1589,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(): @@ -1611,6 +1601,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 diff --git a/pyomo/contrib/doe/examples/reactor_example.py b/pyomo/contrib/doe/examples/reactor_example.py index 0fe81c723dd..0038ed114de 100644 --- a/pyomo/contrib/doe/examples/reactor_example.py +++ b/pyomo/contrib/doe/examples/reactor_example.py @@ -20,7 +20,29 @@ # Example for sensitivity analysis on the reactor experiment # After sensitivity analysis is done, we perform optimal DoE -def run_reactor_doe(): +def run_reactor_doe( + n_points_for_design=9, + compute_FIM_full_factorial=True, + plot_factorial_results=True, + save_plots=True, + run_optimal_doe=True, +): + """ + This function demonstrates how to perform sensitivity analysis on the reactor + + Parameters + ---------- + n_points_for_design : int, optional + number of points to use for the design ranges, by default 9 + compute_FIM_full_factorial : bool, optional + whether to compute the full factorial design, by default True + plot_factorial_results : bool, optional + whether to plot the results of the full factorial design, by default True + save_plots : bool, optional + whether to save draw_factorial_figure plots, by default True + run_optimal_doe : bool, optional + whether to run the optimal DoE, by default True + """ # Read in file DATA_DIR = pathlib.Path(__file__).parent file_path = DATA_DIR / "result.json" @@ -66,63 +88,76 @@ def run_reactor_doe(): _Cholesky_option=True, _only_compute_fim_lower=True, ) - - # Make design ranges to compute the full factorial design - design_ranges = {"CA[0]": [1, 5, 9], "T[0]": [300, 700, 9]} - - # Compute the full factorial design with the sequential FIM calculation - doe_obj.compute_FIM_full_factorial(design_ranges=design_ranges, method="sequential") - - # Plot the results - doe_obj.draw_factorial_figure( - sensitivity_design_variables=["CA[0]", "T[0]"], - fixed_design_variables={ - "T[0.125]": 300, - "T[0.25]": 300, - "T[0.375]": 300, - "T[0.5]": 300, - "T[0.625]": 300, - "T[0.75]": 300, - "T[0.875]": 300, - "T[1]": 300, - }, - title_text="Reactor Example", - xlabel_text="Concentration of A (M)", - ylabel_text="Initial Temperature (K)", - figure_file_name="example_reactor_compute_FIM", - log_scale=False, - ) + if compute_FIM_full_factorial: + # Make design ranges to compute the full factorial design + design_ranges = { + "CA[0]": [1, 5, n_points_for_design], + "T[0]": [300, 700, n_points_for_design], + } + + # Compute the full factorial design with the sequential FIM calculation + doe_obj.compute_FIM_full_factorial( + design_ranges=design_ranges, method="sequential" + ) + if plot_factorial_results: + if save_plots: + figure_file_name = "example_reactor_compute_FIM" + else: + figure_file_name = None + + # Plot the results + doe_obj.draw_factorial_figure( + sensitivity_design_variables=["CA[0]", "T[0]"], + fixed_design_variables={ + "T[0.125]": 300, + "T[0.25]": 300, + "T[0.375]": 300, + "T[0.5]": 300, + "T[0.625]": 300, + "T[0.75]": 300, + "T[0.875]": 300, + "T[1]": 300, + }, + title_text="Reactor Example", + xlabel_text="Concentration of A (M)", + ylabel_text="Initial Temperature (K)", + figure_file_name=figure_file_name, + log_scale=False, + ) ########################### # End sensitivity analysis # Begin optimal DoE #################### - doe_obj.run_doe() - - # Print out a results summary - print("Optimal experiment values: ") - print( - "\tInitial concentration: {:.2f}".format( - doe_obj.results["Experiment Design"][0] + if run_optimal_doe: + doe_obj.run_doe() + + # Print out a results summary + print("Optimal experiment values: ") + print( + "\tInitial concentration: {:.2f}".format( + doe_obj.results["Experiment Design"][0] + ) ) - ) - print( - ("\tTemperature values: [" + "{:.2f}, " * 8 + "{:.2f}]").format( - *doe_obj.results["Experiment Design"][1:] + print( + ("\tTemperature values: [" + "{:.2f}, " * 8 + "{:.2f}]").format( + *doe_obj.results["Experiment Design"][1:] + ) ) - ) - print("FIM at optimal design:\n {}".format(np.array(doe_obj.results["FIM"]))) - print( - "Objective value at optimal design: {:.2f}".format( - pyo.value(doe_obj.model.objective) + print("FIM at optimal design:\n {}".format(np.array(doe_obj.results["FIM"]))) + print( + "Objective value at optimal design: {:.2f}".format( + pyo.value(doe_obj.model.objective) + ) ) - ) - print(doe_obj.results["Experiment Design Names"]) + print(doe_obj.results["Experiment Design Names"]) + + ################### + # End optimal DoE - ################### - # End optimal DoE + return doe_obj if __name__ == "__main__": diff --git a/pyomo/contrib/doe/tests/test_doe_errors.py b/pyomo/contrib/doe/tests/test_doe_errors.py index 4c9823a251d..ce77fb4f553 100644 --- a/pyomo/contrib/doe/tests/test_doe_errors.py +++ b/pyomo/contrib/doe/tests/test_doe_errors.py @@ -190,7 +190,7 @@ def test_reactor_check_bad_prior_negative_eigenvalue(self): doe_obj.create_doe_model() def test_reactor_check_bad_prior_not_symmetric(self): - from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_SYMMETRY + from pyomo.contrib.doe.utils import _SMALL_TOLERANCE_SYMMETRY fd_method = "central" obj_used = "trace" diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py index 83884be0b2c..ea66093253b 100644 --- a/pyomo/contrib/doe/tests/test_doe_solve.py +++ b/pyomo/contrib/doe/tests/test_doe_solve.py @@ -23,8 +23,10 @@ import pyomo.common.unittest as unittest from pyomo.contrib.doe import DesignOfExperiments +from pyomo.contrib.doe.examples.reactor_experiment import ReactorExperiment from pyomo.contrib.doe.examples.reactor_example import ( ReactorExperiment as FullReactorExperiment, + run_reactor_doe, ) from pyomo.contrib.doe.tests.experiment_class_example_flags import ( FullReactorExperimentBad, @@ -37,7 +39,7 @@ ipopt_available = SolverFactory("ipopt").available() -k_aug_available = SolverFactory('k_aug', solver_io='nl', validate=False) +k_aug_available = SolverFactory("k_aug", solver_io="nl", validate=False) currdir = this_file_dir() file_path = os.path.join(currdir, "..", "examples", "result.json") @@ -100,21 +102,21 @@ def get_FIM_Q_L(doe_obj=None): def get_standard_args(experiment, fd_method, obj_used): args = {} - args['experiment'] = experiment - args['fd_formula'] = fd_method - args['step'] = 1e-3 - args['objective_option'] = obj_used - args['scale_constant_value'] = 1 - args['scale_nominal_param_value'] = True - args['prior_FIM'] = None - args['jac_initial'] = None - args['fim_initial'] = None - args['L_diagonal_lower_bound'] = 1e-7 - args['solver'] = None - args['tee'] = False - args['get_labeled_model_args'] = None - args['_Cholesky_option'] = True - args['_only_compute_fim_lower'] = True + args["experiment"] = experiment + args["fd_formula"] = fd_method + args["step"] = 1e-3 + args["objective_option"] = obj_used + args["scale_constant_value"] = 1 + args["scale_nominal_param_value"] = True + args["prior_FIM"] = None + args["jac_initial"] = None + args["fim_initial"] = None + args["L_diagonal_lower_bound"] = 1e-7 + args["solver"] = None + args["tee"] = False + args["get_labeled_model_args"] = None + args["_Cholesky_option"] = True + args["_only_compute_fim_lower"] = True return args @@ -195,17 +197,17 @@ def test_reactor_obj_det_solve(self): experiment = FullReactorExperiment(data_ex, 10, 3) DoE_args = get_standard_args(experiment, fd_method, obj_used) - DoE_args['scale_nominal_param_value'] = ( + DoE_args["scale_nominal_param_value"] = ( False # Vanilla determinant solve needs this ) - DoE_args['_Cholesky_option'] = False - DoE_args['_only_compute_fim_lower'] = False + DoE_args["_Cholesky_option"] = False + DoE_args["_only_compute_fim_lower"] = False doe_obj = DesignOfExperiments(**DoE_args) doe_obj.run_doe() - self.assertEqual(doe_obj.results['Solver Status'], "ok") + self.assertEqual(doe_obj.results["Solver Status"], "ok") def test_reactor_obj_cholesky_solve(self): fd_method = "central" @@ -244,7 +246,7 @@ def test_reactor_obj_cholesky_solve_bad_prior(self): # Specify a prior that is slightly negative definite # Because it is less than the tolerance, it should be adjusted to be positive definite # No error should be thrown - DoE_args['prior_FIM'] = -(_SMALL_TOLERANCE_DEFINITENESS / 100) * np.eye(4) + DoE_args["prior_FIM"] = -(_SMALL_TOLERANCE_DEFINITENESS / 100) * np.eye(4) doe_obj = DesignOfExperiments(**DoE_args) @@ -366,7 +368,7 @@ def test_rescale_FIM(self): # Without parameter scaling DoE_args2 = get_standard_args(experiment, fd_method, obj_used) - DoE_args2['scale_nominal_param_value'] = False + DoE_args2["scale_nominal_param_value"] = False doe_obj2 = DesignOfExperiments(**DoE_args2) # Run both problems @@ -418,7 +420,7 @@ def test_reactor_grid_search_bad_model(self): experiment = FullReactorExperimentBad(data_ex, 10, 3) DoE_args = get_standard_args(experiment, fd_method, obj_used) - DoE_args['logger_level'] = logging.ERROR + DoE_args["logger_level"] = logging.ERROR doe_obj = DesignOfExperiments(**DoE_args) @@ -443,5 +445,103 @@ def test_reactor_grid_search_bad_model(self): ) +@unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +class TestDoe(unittest.TestCase): + def test_doe_full_factorial(self): + log10_D_opt_expected = [ + 3.7734377852467524, + 5.137792359070963, + 5.182167857710023, + 6.546522431509408, + ] + + log10_A_opt_expected = [ + 3.5935726800929695, + 3.6133186151486933, + 3.945755198204365, + 3.9655011332598367, + ] + + log10_E_opt_expected = [ + -1.7201873126109162, + -0.691340497355524, + -1.3680047944877138, + -0.3391579792516522, + ] + + log10_ME_opt_expected = [ + 5.221185311075697, + 4.244741560076784, + 5.221185311062606, + 4.244741560083524, + ] + + eigval_min_expected = [ + 0.019046390638130666, + 0.20354456134677426, + 0.04285437893696232, + 0.45797526302234304, + ] + + eigval_max_expected = [ + 3169.552855492114, + 3576.0292523637977, + 7131.493924857995, + 8046.0658178139165, + ] + + det_FIM_expected = [ + 5935.233170586055, + 137338.51875774842, + 152113.5345070818, + 3519836.021699428, + ] + + trace_FIM_expected = [ + 3922.5878617108597, + 4105.051549241871, + 8825.822688850109, + 9236.36598578955, + ] + ff = run_reactor_doe( + n_points_for_design=2, + compute_FIM_full_factorial=False, + plot_factorial_results=False, + save_plots=False, + run_optimal_doe=False, + ) + ff.compute_FIM_full_factorial( + design_ranges={"CA[0]": [1, 1.5, 2], "T[0]": [350, 400, 2]} + ) + + ff_results = ff.fim_factorial_results + + self.assertStructuredAlmostEqual( + ff_results["log10 D-opt"], log10_D_opt_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["log10 A-opt"], log10_A_opt_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["log10 E-opt"], log10_E_opt_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["log10 ME-opt"], log10_ME_opt_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["eigval_min"], eigval_min_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["eigval_max"], eigval_max_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["det_FIM"], det_FIM_expected, abstol=1e-4 + ) + self.assertStructuredAlmostEqual( + ff_results["trace_FIM"], trace_FIM_expected, abstol=1e-4 + ) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py new file mode 100644 index 00000000000..47df75c401c --- /dev/null +++ b/pyomo/contrib/doe/tests/test_utils.py @@ -0,0 +1,149 @@ +# ___________________________________________________________________________ +# +# 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.utils import ( + check_FIM, + compute_FIM_metrics, + get_FIM_metrics, + _SMALL_TOLERANCE_DEFINITENESS, + _SMALL_TOLERANCE_SYMMETRY, + _SMALL_TOLERANCE_IMG, +) + + +@unittest.skipIf(not numpy_available, "Numpy is not available") +class TestUtilsFIM(unittest.TestCase): + """Test the check_FIM() from utils.py.""" + + 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: + 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"): + 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, + "FIM provided is not positive definite. It has one or more negative " + + r"eigenvalue\(s\) less than -{:.1e}".format( + _SMALL_TOLERANCE_DEFINITENESS + ), + ): + 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 + ), + ): + check_FIM(FIM) + + """Test the compute_FIM_metrics() from utils.py.""" + + 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 = 26.000000000000004 + D_opt_expected = 1.414973347970818 + + trace_expected = 13 + A_opt_expected = 1.1139433523068367 + + E_vals_expected = np.array([10.53112887, 2.46887113]) + E_vecs_expected = np.array( + [[0.96649965, -0.25666794], [0.25666794, 0.96649965]] + ) + E_opt_expected = 0.3924984205140895 + + ME_opt_expected = 0.6299765069426388 + + # Test results + self.assertAlmostEqual(det_FIM, det_expected) + self.assertAlmostEqual(trace_FIM, trace_expected) + self.assertTrue(np.allclose(E_vals, E_vals_expected)) + self.assertTrue(np.allclose(E_vecs, E_vecs_expected)) + self.assertAlmostEqual(D_opt, D_opt_expected) + self.assertAlmostEqual(A_opt, A_opt_expected) + self.assertAlmostEqual(E_opt, E_opt_expected) + self.assertAlmostEqual(ME_opt, ME_opt_expected) + + def test_FIM_eigenvalue_warning(self): + # 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.utils", level="WARNING") as cm: + compute_FIM_metrics(FIM) + expected_warning = ( + "Eigenvalue has imaginary component greater than " + + f"{_SMALL_TOLERANCE_IMG}, contact the developers if this issue " + + "persists." + ) + self.assertIn(expected_warning, cm.output[0]) + + """Test the get_FIM_metrics() from utils.py.""" + + def test_get_FIM_metrics(self): + # Create a sample Fisher Information Matrix (FIM) + FIM = np.array([[10, 2], [2, 3]]) + fim_metrics = get_FIM_metrics(FIM) + + # expected results + det_expected = 26.000000000000004 + D_opt_expected = 1.414973347970818 + + trace_expected = 13 + A_opt_expected = 1.1139433523068367 + + E_vals_expected = np.array([10.53112887, 2.46887113]) + E_vecs_expected = np.array( + [[0.96649965, -0.25666794], [0.25666794, 0.96649965]] + ) + E_opt_expected = 0.3924984205140895 + + ME_opt_expected = 0.6299765069426388 + + # Test results + self.assertAlmostEqual(fim_metrics["Determinant of FIM"], det_expected) + self.assertAlmostEqual(fim_metrics["Trace of FIM"], trace_expected) + self.assertTrue(np.allclose(fim_metrics["Eigenvalues"], E_vals_expected)) + self.assertTrue(np.allclose(fim_metrics["Eigenvectors"], E_vecs_expected)) + self.assertAlmostEqual(fim_metrics["log10(D-Optimality)"], D_opt_expected) + self.assertAlmostEqual(fim_metrics["log10(A-Optimality)"], A_opt_expected) + self.assertAlmostEqual(fim_metrics["log10(E-Optimality)"], E_opt_expected) + self.assertAlmostEqual( + fim_metrics["log10(Modified E-Optimality)"], ME_opt_expected + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 24be6fd696a..26b991a753a 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -32,6 +32,30 @@ from pyomo.core.base.param import ParamData from pyomo.core.base.var import VarData +import logging + +logger = logging.getLogger(__name__) + +# This small and positive tolerance is used when checking +# if the prior is negative definite or approximately +# indefinite. 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_DEFINITENESS = 1e-6 + +# This small and positive tolerance is used to check +# the FIM is approximately symmetric. 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_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 + # Rescale FIM (a scaling function to help rescale FIM from parameter values) def rescale_FIM(FIM, param_vals): @@ -54,9 +78,8 @@ def rescale_FIM(FIM, param_vals): (len(param_vals.shape) == 2) and (param_vals.shape[0] != 1) ): raise ValueError( - "param_vals should be a vector of dimensions: 1 by `n_params`. The shape you provided is {}.".format( - param_vals.shape - ) + "param_vals should be a vector of dimensions: 1 by `n_params`. " + + "The shape you provided is {}.".format(param_vals.shape) ) if len(param_vals.shape) == 1: param_vals = np.array([param_vals]) @@ -100,3 +123,146 @@ def rescale_FIM(FIM, param_vals): # pass # ToDo: Write error for suffix keys that aren't ParamData or VarData # # return param_list + + +def check_FIM(FIM): + """ + Checks 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 + """ + # Check 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) + + # Check if the FIM is positive definite + if np.min(evals) < -_SMALL_TOLERANCE_DEFINITENESS: + raise ValueError( + "FIM provided is not positive definite. It has one or more negative " + + "eigenvalue(s) less than -{:.1e}".format(_SMALL_TOLERANCE_DEFINITENESS) + ) + + # Check if the FIM is symmetric + if not np.allclose(FIM, FIM.T, atol=_SMALL_TOLERANCE_SYMMETRY): + raise ValueError( + "FIM provided is not symmetric using absolute tolerance {}".format( + _SMALL_TOLERANCE_SYMMETRY + ) + ) + + +# Functions to compute FIM metrics +def compute_FIM_metrics(FIM): + """ + Parameters + ---------- + FIM : numpy.ndarray + 2D array representing the Fisher Information Matrix (FIM). + + Returns + ------- + Returns the following metrics as a tuple in the order shown below: + + det_FIM : float + Determinant of the FIM. + trace_FIM : float + Trace of the FIM. + E_vals : numpy.ndarray + 1D array of eigenvalues of the FIM. + E_vecs : numpy.ndarray + 2D array of eigenvectors of the FIM. + D_opt : float + log10(D-optimality) metric. + A_opt : float + log10(A-optimality) metric. + E_opt : float + log10(E-optimality) metric. + ME_opt : float + log10(Modified E-optimality) metric. + """ + + # Check whether the FIM is square, positive definite, and symmetric + 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) + if abs(E_vals.imag[E_ind]) > _SMALL_TOLERANCE_IMG: + logger.warning( + "Eigenvalue has imaginary component greater than " + + f"{_SMALL_TOLERANCE_IMG}, contact the 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 : numpy.ndarray + 2D numpy array of the FIM + + Returns + ------- + A dictionary containing the following keys: + + "Determinant of FIM" : float + determinant of the FIM + "Trace of FIM" : float + trace of the FIM + "Eigenvalues" : numpy.ndarray + eigenvalues of the FIM + "Eigenvectors" : numpy.ndarray + eigenvectors of the FIM + "log10(D-Optimality)" : float + log10(D-optimality) metric + "log10(A-Optimality)" : float + log10(A-optimality) metric + "log10(E-Optimality)" : float + log10(E-optimality) metric + "log10(Modified E-Optimality)" : float + log10(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 { + "Determinant of FIM": det_FIM, + "Trace of FIM": trace_FIM, + "Eigenvalues": E_vals, + "Eigenvectors": E_vecs, + "log10(D-Optimality)": D_opt, + "log10(A-Optimality)": A_opt, + "log10(E-Optimality)": E_opt, + "log10(Modified E-Optimality)": ME_opt, + }