From a4a1adea92e6a745e7997687ab8aa9f5b42ed3d7 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Mon, 16 Jun 2025 16:29:55 -0400 Subject: [PATCH 01/25] added zigzag in utils and compute fim in doe.py --- pyomo/contrib/doe/doe.py | 144 +++++++++++++++++++++++++++++++++++++ pyomo/contrib/doe/utils.py | 57 +++++++++++++++ 2 files changed, 201 insertions(+) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 76c154b07ee..5844920808e 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -45,6 +45,7 @@ from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp import pyomo.environ as pyo +from pyomo.contrib.doe.utils import generate_snake_zigzag_pattern from pyomo.opt import SolverStatus @@ -1617,6 +1618,149 @@ def compute_FIM_full_factorial( return self.fim_factorial_results + def compute_FIM_factorial( + self, model=None, design_ranges=None, method="sequential" + ): + """ + Will run a simulation-based factorial exploration of + the experimental input space (i.e., a ``grid search`` or + ``parameter sweep``) to understand how the FIM metrics + change as a function of the experimental design space. + + Parameters + ---------- + model: model to perform the full factorial exploration on + design_ranges: dict of lists, of the form {: []} + method: string to specify which method should be used + options are ``kaug`` and ``sequential`` + + """ + # Start timer + sp_timer = TicTocTimer() + sp_timer.tic(msg=None) + self.logger.info("Beginning Factorial Design.") + + # Make new model for factorial design + self.factorial_model = self.experiment.get_labeled_model( + **self.get_labeled_model_args + ).clone() + model = self.factorial_model + + # Permute the inputs to be aligned with the experiment input indices + design_ranges_enum = design_ranges # {k: v for k, v in design_ranges.items()} + design_map = { + ind: (k[0].name, k[0]) + for ind, k in enumerate(model.experiment_inputs.items()) + } + + # Check whether the design_ranges keys are in the experiment_inputs + design_keys = set(design_ranges.keys()) + map_keys = set([k.name for k in model.experiment_inputs.keys]) + + if not design_keys.issubset(map_keys): + raise ValueError( + f"design_ranges keys {design_keys} must be a subset of experimental\ + design names: {map_keys}." + ) + + des_ranges = [design_ranges[k] for k in design_ranges.keys()] + + factorial_points = generate_snake_zigzag_pattern(*des_ranges) + factorial_results = {k.name: [] for k in model.experiment_inputs.keys()} + factorial_results.update( + { + "log10 D-opt": [], + "log10 A-opt": [], + "log10 E-opt": [], + "log10 ME-opt": [], + "solve_time": [], + } + ) + + success_count = 0 + failure_count = 0 + total_points = len(factorial_points) + + time_set = [] + curr_point = 1 # Initial current point + for design_point in factorial_points: + # Fix design variables at fixed experimental design point + for i in range(len(design_point)): + design_map[i][1].fix(design_point[i]) + + self.logger.info(f"=======Iteration Number: {curr_point} ======") + iter_timer = TicTocTimer() + iter_timer.tic(msg=None) + + try: + curr_point = success_count + failure_count + 1 + self.logger.info(f"This is run {curr_point} out of {total_points}.") + self.compute_FIM(model=model, method=method) + success_count += 1 + # iteration time + iter_t = iter_timer.toc(msg=None) + time_set.append(iter_t) + + # More logging + self.logger.info( + f"The code has run for {round(sum(time_set), 2)} seconds." + ) + self.logger.info( + "Estimated remaining time: %s seconds", + round( + sum(time_set) / (curr_point) * (total_points - curr_point + 1), + 2, + ), + ) + except: + self.logger.warning( + ":::::::::::Warning: Cannot converge this run.::::::::::::" + ) + failures += 1 + self.logger.warning("failed count:", failures) + + self._computed_FIM = np.zeros(self.prior_FIM.shape) + + iter_t = iter_timer.toc(msg=None) + time_set.append(iter_t) + + 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)) + + for k, v in model.experiment_inputs.items(): + factorial_results[k.name].append(pyo.value(k)) + + factorial_results["log10 D-opt"].append(D_opt) + factorial_results["log10 A-opt"].append(A_opt) + factorial_results["log10 E-opt"].append(E_opt) + factorial_results["log10 ME-opt"].append(ME_opt) + factorial_results["solve_time"].append(time_set[-1]) + + self.factorial_results = factorial_results + return self.factorial_results + + # for k in design_ranges.keys(): + # if k not in [k2 for k2 in model.experiment_inputs.keys()]: + # raise ValueError( + # TODO: Overhaul plotting functions to not use strings # TODO: Make the plotting functionalities work for >2 design features def draw_factorial_figure( diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 24be6fd696a..e2a018dab14 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -100,3 +100,60 @@ def rescale_FIM(FIM, param_vals): # pass # ToDo: Write error for suffix keys that aren't ParamData or VarData # # return param_list + + +def generate_snake_zigzag_pattern(*lists): + """ + Generates a multi-dimensional zigzag pattern for an arbitrary number of lists. + This pattern is useful for generating patterns for sensitivity analysis when we want + to change one variable at a time. + + This function uses recursion and acts as a generator. The direction of + iteration for each list is determined by the sum of the indices from the + preceding lists. If the sum is even, the list is iterated forwards; + if it's odd, the list is iterated in reverse. + + Parameters + ---------- + lists: A variable number of 1D list arguments. + + Yields + ------ + A tuple representing points in the snake-like zigzag pattern. + """ + + # The main logic is in a nested recursive helper function. + def _generate_recursive(depth, index_sum): + """ + Parameters + ---------- + depth : int + corresponds to the index of the current list in the `lists` argument. + Represents which list we are currently processing. + + index_sum : int + It's a running total of the indices chosen from all the previous lists. + The direction of the depth-th list depends entirely on whether index_sum + even or odd. + """ + # Base case: If we've processed all lists, we're at the end of a path. + if depth == len(lists): + yield () + return + + current_list = lists[depth] + + # Determine the iteration direction based on the sum of parent indices. + # This is the mathematical rule for the zigzag. + is_forward = index_sum % 2 == 0 + iterable = current_list if is_forward else reversed(current_list) + + # Enumerate to get the index `i` for the *next* recursive call's sum. + for i, value in enumerate(iterable): + # Recur for the next list, updating the index_sum. + for sub_pattern in _generate_recursive(depth + 1, index_sum + i): + # Prepend the current value to the results from deeper levels. + yield (value,) + sub_pattern + + # Start the recursion at the first list (depth 0) with an initial sum of 0. + yield from _generate_recursive(depth=0, index_sum=0) From 26f8cd34a4e1895a6ea0650501b1dd6e330c12d0 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 17 Jun 2025 00:57:15 -0400 Subject: [PATCH 02/25] The compute_FIM_factorial() using zigzag pattern is now working as expected. Now, I need to add cases when some of the design variables are not changed. --- pyomo/contrib/doe/doe.py | 49 +++++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 16 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 5844920808e..4238f41374b 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1619,7 +1619,11 @@ def compute_FIM_full_factorial( return self.fim_factorial_results def compute_FIM_factorial( - self, model=None, design_ranges=None, method="sequential" + self, + model=None, + design_values: dict = None, + method="sequential", + file_name: str = None, ): """ Will run a simulation-based factorial exploration of @@ -1630,7 +1634,7 @@ def compute_FIM_factorial( Parameters ---------- model: model to perform the full factorial exploration on - design_ranges: dict of lists, of the form {: []} + design_values: dict of lists, of the form {: []} method: string to specify which method should be used options are ``kaug`` and ``sequential`` @@ -1647,25 +1651,31 @@ def compute_FIM_factorial( model = self.factorial_model # Permute the inputs to be aligned with the experiment input indices - design_ranges_enum = design_ranges # {k: v for k, v in design_ranges.items()} design_map = { ind: (k[0].name, k[0]) for ind, k in enumerate(model.experiment_inputs.items()) } # Check whether the design_ranges keys are in the experiment_inputs - design_keys = set(design_ranges.keys()) - map_keys = set([k.name for k in model.experiment_inputs.keys]) + design_keys = set(design_values.keys()) + map_keys = set([k.name for k in model.experiment_inputs.keys()]) + print(f"design_keys: {design_keys}, \nmap_keys: {map_keys}") - if not design_keys.issubset(map_keys): + if design_keys != map_keys: + incorrect_given_keys = design_keys - map_keys + incorrect_map_keys = map_keys - design_keys raise ValueError( - f"design_ranges keys {design_keys} must be a subset of experimental\ - design names: {map_keys}." + f"design_values key(s): {incorrect_given_keys} is incorrect." + f"Please provide values for the following key(s): {incorrect_map_keys}." ) - des_ranges = [design_ranges[k] for k in design_ranges.keys()] + des_ranges = [design_values[k] for k in design_values.keys()] + print(f"des_ranges: {des_ranges}") factorial_points = generate_snake_zigzag_pattern(*des_ranges) + factorial_points_list = list(factorial_points) + print("factorial_points:", factorial_points_list) + factorial_results = {k.name: [] for k in model.experiment_inputs.keys()} factorial_results.update( { @@ -1679,11 +1689,13 @@ def compute_FIM_factorial( success_count = 0 failure_count = 0 - total_points = len(factorial_points) + total_points = len(factorial_points_list) + print(f"Total points: {total_points}") time_set = [] curr_point = 1 # Initial current point - for design_point in factorial_points: + for design_point in factorial_points_list: + print(f"design_point: {design_point}") # Fix design variables at fixed experimental design point for i in range(len(design_point)): design_map[i][1].fix(design_point[i]) @@ -1696,6 +1708,9 @@ def compute_FIM_factorial( curr_point = success_count + failure_count + 1 self.logger.info(f"This is run {curr_point} out of {total_points}.") self.compute_FIM(model=model, method=method) + print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") + print("y_hat", model.y.value) + print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") success_count += 1 # iteration time iter_t = iter_timer.toc(msg=None) @@ -1754,12 +1769,14 @@ def compute_FIM_factorial( factorial_results["log10 ME-opt"].append(ME_opt) factorial_results["solve_time"].append(time_set[-1]) - self.factorial_results = factorial_results - return self.factorial_results + self.factorial_results = factorial_results + + if file_name is not None: + with open(file_name + ".json", "w") as f: + json.dump(self.factorial_results, f, indent=4) # Save results to file + self.logger.info(f"Results saved to {file_name}.json.") - # for k in design_ranges.keys(): - # if k not in [k2 for k2 in model.experiment_inputs.keys()]: - # raise ValueError( + return self.factorial_results # TODO: Overhaul plotting functions to not use strings # TODO: Make the plotting functionalities work for >2 design features From 55bfeae82b395206f9f50b265be7e42024d03a4f Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 17 Jun 2025 17:16:24 -0400 Subject: [PATCH 03/25] Added utils.py as utils_updated.py from adding_eigen_values branch to use compute_FIM_metrics() method. Added both the nested for loop and change_one_design_at_a_time argument to compute_FIM_factorial() in doe.py. Everything is working fine. --- pyomo/contrib/doe/doe.py | 139 +++++++++------ pyomo/contrib/doe/utils_updated.py | 267 +++++++++++++++++++++++++++++ 2 files changed, 356 insertions(+), 50 deletions(-) create mode 100644 pyomo/contrib/doe/utils_updated.py diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 4238f41374b..0acadaf230c 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -46,6 +46,7 @@ import pyomo.environ as pyo from pyomo.contrib.doe.utils import generate_snake_zigzag_pattern +from pyomo.contrib.doe.utils_updated import compute_FIM_metrics from pyomo.opt import SolverStatus @@ -1623,22 +1624,56 @@ def compute_FIM_factorial( model=None, design_values: dict = None, method="sequential", + change_one_design_at_a_time=True, file_name: str = None, ): - """ - Will run a simulation-based factorial exploration of - the experimental input space (i.e., a ``grid search`` or - ``parameter sweep``) to understand how the FIM metrics - change as a function of the experimental design space. + """Will run a simulation-based factorial exploration of the experimental input + space (i.e., a ``grid search`` or ``parameter sweep``) to understand how the + FIM metrics change as a function of the experimental design space. This method + can be used for both full factorial and fractional factorial designs. Parameters ---------- - model: model to perform the full factorial exploration on - design_values: dict of lists, of the form {: []} - method: string to specify which method should be used - options are ``kaug`` and ``sequential`` + model : DoE model, optional + The model to perform the full factorial exploration on. Default: None + design_values : dict, optional + dict of lists, of the form {: []}. Default: None. The + decision variables should be passed in the same order as the model's + `experiment_inputs`. If a particular decision variable is not to be changed, + it must be set to `None`. For example, design_values= {"x1": [1, 2, 3], + "x2": [None], "x3": [7, 8], "x4": [-10]} + method : str, optional + string to specify which method should be used. options are ``kaug`` and + ``sequential`. Default: "sequential" + change_one_design_at_a_time : bool, optional + If True, will generate a snake-zigzag pattern of the design values that + changes only one of the decision variables at a time. If False, will + generate a regular nested for loop that can change one to all the design + values at a time. Default: True + file_name : str, optional + if provided, will save the results to a json file. Default: None + Returns + ------- + factorial_results: dict + A dictionary containing the results of the factorial design with the + following keys: + - keys of model's experiment_inputs + - "log10(D-opt)": list of D-optimality values + - "log10(A-opt)": list of A-optimality values + - "log10(E-opt)": list of E-optimality values + - "log10(ME-opt)": list of ME-optimality values + - "solve_time": list of solve times + - "total_points": total number of points in the factorial design + - "success_count": number of successful runs + - "failure_count": number of failed runs + + Raises + ------ + ValueError + If the design_values' keys do not match the model's experiment_inputs' keys. """ + # Start timer sp_timer = TicTocTimer() sp_timer.tic(msg=None) @@ -1659,30 +1694,34 @@ def compute_FIM_factorial( # Check whether the design_ranges keys are in the experiment_inputs design_keys = set(design_values.keys()) map_keys = set([k.name for k in model.experiment_inputs.keys()]) - print(f"design_keys: {design_keys}, \nmap_keys: {map_keys}") if design_keys != map_keys: incorrect_given_keys = design_keys - map_keys - incorrect_map_keys = map_keys - design_keys + unmatched_map_keys = map_keys - design_keys raise ValueError( - f"design_values key(s): {incorrect_given_keys} is incorrect." - f"Please provide values for the following key(s): {incorrect_map_keys}." + f"design_values key(s): {incorrect_given_keys} are incorrect." + f"Please provide values for the following key(s): {unmatched_map_keys}." ) des_ranges = [design_values[k] for k in design_values.keys()] - print(f"des_ranges: {des_ranges}") + if change_one_design_at_a_time: + factorial_points = generate_snake_zigzag_pattern(*des_ranges) + else: + factorial_points = product(*des_ranges) - factorial_points = generate_snake_zigzag_pattern(*des_ranges) factorial_points_list = list(factorial_points) - print("factorial_points:", factorial_points_list) factorial_results = {k.name: [] for k in model.experiment_inputs.keys()} factorial_results.update( { - "log10 D-opt": [], - "log10 A-opt": [], - "log10 E-opt": [], - "log10 ME-opt": [], + "log10(D-opt)": [], + "log10(A-opt)": [], + "log10(E-opt)": [], + "log10(ME-opt)": [], + "eigval_min": [], + "eigval_max": [], + "det_FIM": [], + "trace_FIM": [], "solve_time": [], } ) @@ -1690,7 +1729,6 @@ def compute_FIM_factorial( success_count = 0 failure_count = 0 total_points = len(factorial_points_list) - print(f"Total points: {total_points}") time_set = [] curr_point = 1 # Initial current point @@ -1698,18 +1736,22 @@ def compute_FIM_factorial( print(f"design_point: {design_point}") # Fix design variables at fixed experimental design point for i in range(len(design_point)): - design_map[i][1].fix(design_point[i]) + if design_point[i] is not None: + design_map[i][1].fix(design_point[i]) + else: + pass - self.logger.info(f"=======Iteration Number: {curr_point} ======") - iter_timer = TicTocTimer() - iter_timer.tic(msg=None) + # Timing and logging objects + self.logger.info(f"=======Iteration Number: {curr_point} =======") + iter_timer = TicTocTimer() + iter_timer.tic(msg=None) try: curr_point = success_count + failure_count + 1 self.logger.info(f"This is run {curr_point} out of {total_points}.") self.compute_FIM(model=model, method=method) print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") - print("y_hat", model.y.value) + print(model.pprint()) print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") success_count += 1 # iteration time @@ -1731,8 +1773,8 @@ def compute_FIM_factorial( self.logger.warning( ":::::::::::Warning: Cannot converge this run.::::::::::::" ) - failures += 1 - self.logger.warning("failed count:", failures) + failure_count += 1 + self.logger.warning("failed count:", failure_count) self._computed_FIM = np.zeros(self.prior_FIM.shape) @@ -1742,33 +1784,30 @@ def compute_FIM_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) + ) - for k, v in model.experiment_inputs.items(): + for k in model.experiment_inputs.keys(): factorial_results[k.name].append(pyo.value(k)) - factorial_results["log10 D-opt"].append(D_opt) - factorial_results["log10 A-opt"].append(A_opt) - factorial_results["log10 E-opt"].append(E_opt) - factorial_results["log10 ME-opt"].append(ME_opt) + factorial_results["log10(D-opt)"].append(D_opt) + factorial_results["log10(A-opt)"].append(A_opt) + factorial_results["log10(E-opt)"].append(E_opt) + factorial_results["log10(ME-opt)"].append(ME_opt) + factorial_results["eigval_min"].append(np.min(E_vals)) + factorial_results["eigval_max"].append(np.max(E_vals)) + factorial_results["det_FIM"].append(det_FIM) + factorial_results["trace_FIM"].append(trace_FIM) factorial_results["solve_time"].append(time_set[-1]) + factorial_results.update( + { + "total_points": total_points, + "success_counts": success_count, + "failure_counts": failure_count, + } + ) self.factorial_results = factorial_results if file_name is not None: diff --git a/pyomo/contrib/doe/utils_updated.py b/pyomo/contrib/doe/utils_updated.py new file mode 100644 index 00000000000..2232fd88c52 --- /dev/null +++ b/pyomo/contrib/doe/utils_updated.py @@ -0,0 +1,267 @@ +# ___________________________________________________________________________ +# +# 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. +# +# Pyomo.DoE was produced under the Department of Energy Carbon Capture Simulation +# Initiative (CCSI), and is copyright (c) 2022 by the software owners: +# TRIAD National Security, LLC., Lawrence Livermore National Security, LLC., +# Lawrence Berkeley National Laboratory, Pacific Northwest National Laboratory, +# Battelle Memorial Institute, University of Notre Dame, +# The University of Pittsburgh, The University of Texas at Austin, +# University of Toledo, West Virginia University, et al. All rights reserved. +# +# NOTICE. This Software was developed under funding from the +# U.S. Department of Energy and the U.S. Government consequently retains +# certain rights. As such, the U.S. Government has been granted for itself +# and others acting on its behalf a paid-up, nonexclusive, irrevocable, +# worldwide license in the Software to reproduce, distribute copies to the +# public, prepare derivative works, and perform publicly and display +# publicly, and to permit other to do so. +# ___________________________________________________________________________ + +import pyomo.environ as pyo + +from pyomo.common.dependencies import numpy as np, numpy_available + +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): + """ + Rescales the FIM based on the input and parameter vals. + It is assumed the parameter vals align with the FIM + dimensions such that (1, i) corresponds to the i-th + column or row of the FIM. + + Parameters + ---------- + FIM: 2D numpy array to be scaled + param_vals: scaling factors for the parameters + + """ + if isinstance(param_vals, list): + param_vals = np.array([param_vals]) + elif isinstance(param_vals, np.ndarray): + if len(param_vals.shape) > 2 or ( + (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) + ) + if len(param_vals.shape) == 1: + param_vals = np.array([param_vals]) + else: + raise ValueError( + "param_vals should be a list or numpy array of dimensions: 1 by `n_params`" + ) + scaling_mat = (1 / param_vals).transpose().dot((1 / param_vals)) + scaled_FIM = np.multiply(FIM, scaling_mat) + return scaled_FIM + + # TODO: Add swapping parameters for variables helper function + # def get_parameters_from_suffix(suffix, fix_vars=False): + # """ + # Finds the Params within the suffix provided. It will also check to see + # if there are Vars in the suffix provided. ``fix_vars`` will indicate + # if we should fix all the Vars in the set or not. + # + # Parameters + # ---------- + # suffix: pyomo Suffix object, contains the components to be checked + # as keys + # fix_vars: boolean, whether or not to fix the Vars, default = False + # + # Returns + # ------- + # param_list: list of Param + # """ + # param_list = [] + # + # # FIX THE MODEL TREE ISSUE WHERE I GET base_model. INSTEAD OF + # # Check keys if they are Param or Var. Fix the vars if ``fix_vars`` is True + # for k, v in suffix.items(): + # if isinstance(k, ParamData): + # param_list.append(k.name) + # elif isinstance(k, VarData): + # if fix_vars: + # k.fix() + # else: + # 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 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, + } From f8105ea733af52d3584c0625ca32cbb42fb5548a Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Thu, 19 Jun 2025 01:02:52 -0400 Subject: [PATCH 04/25] this version of compute_FIM_factorial() runs fine. The argument `design_values` does not need to be in the same order as the `experiment_inputs` as long as the key matches, we are good. Also, if we don't want to change a design variable, just not passing it as a key in the `design_values` dictionary will do the trick. Tested with the `reactor_example.py` --- pyomo/contrib/doe/doe.py | 73 ++++++++++++++++++++++---------------- pyomo/contrib/doe/utils.py | 2 +- 2 files changed, 44 insertions(+), 31 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 0acadaf230c..95a0d75acd6 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1636,20 +1636,24 @@ def compute_FIM_factorial( ---------- model : DoE model, optional The model to perform the full factorial exploration on. Default: None - design_values : dict, optional - dict of lists, of the form {: []}. Default: None. The - decision variables should be passed in the same order as the model's - `experiment_inputs`. If a particular decision variable is not to be changed, - it must be set to `None`. For example, design_values= {"x1": [1, 2, 3], - "x2": [None], "x3": [7, 8], "x4": [-10]} + design_values : dict, + dict of lists, of the form {<"var_name">: []}. Default: None. + The `design_values` should have the same key(s) as the `experiment_inputs`. + If one or more design variables are not to be changed, they do not need to + be passed in the `design_values` dictionary, but if they are passed in the + dictionary, then they must be a list of floats. For example, if our + experiment has 4 design variables (i.e., `experiment_inputs`): model.x1, + model.x2, model.x3, and model.x4, their values may be passed as, + design_values= {"x1": [1, 2, 3], "x3": [7], "x4": [-10, 20]}. In this case, + x2 is not be changed and will be fixed at the default value in the model. method : str, optional string to specify which method should be used. options are ``kaug`` and ``sequential`. Default: "sequential" change_one_design_at_a_time : bool, optional - If True, will generate a snake-zigzag pattern of the design values that - changes only one of the decision variables at a time. If False, will - generate a regular nested for loop that can change one to all the design - values at a time. Default: True + If True, will generate a snake-like zigzag pattern of the design values that + changes only one of the design variables at a time. If False, will + generate a regular nested for loop that can change from one to all the design + variables at a time. Default: True file_name : str, optional if provided, will save the results to a json file. Default: None @@ -1667,11 +1671,14 @@ def compute_FIM_factorial( - "total_points": total number of points in the factorial design - "success_count": number of successful runs - "failure_count": number of failed runs + - "FIM_all": list of all FIMs computed for each point in the factorial + design. Raises ------ ValueError - If the design_values' keys do not match the model's experiment_inputs' keys. + If the design_values' keys is not a subset of the model's experiment_inputs' + keys or if the design_values are not provided as a dictionary of lists. """ # Start timer @@ -1685,25 +1692,29 @@ def compute_FIM_factorial( ).clone() model = self.factorial_model - # Permute the inputs to be aligned with the experiment input indices - design_map = { - ind: (k[0].name, k[0]) - for ind, k in enumerate(model.experiment_inputs.items()) - } + if not design_values: + raise ValueError( + "design_values must be provided as a dictionary of lists " + "in the form {<'var_name'>: []}." + ) # Check whether the design_ranges keys are in the experiment_inputs design_keys = set(design_values.keys()) map_keys = set([k.name for k in model.experiment_inputs.keys()]) - - if design_keys != map_keys: + if not design_keys.issubset(map_keys): incorrect_given_keys = design_keys - map_keys - unmatched_map_keys = map_keys - design_keys raise ValueError( - f"design_values key(s): {incorrect_given_keys} are incorrect." - f"Please provide values for the following key(s): {unmatched_map_keys}." + f"design_values keys: {incorrect_given_keys} are incorrect." + f"The key should be from the following: {map_keys}." ) - des_ranges = [design_values[k] for k in design_values.keys()] + # Get the design map keys that match the design_values keys + design_map_keys = [ + k for k in model.experiment_inputs.keys() if k.name in design_values.keys() + ] + # This ensures that the order of the design_values keys matches the order of the + # design_map_keys so that design_point can be constructed correctly in the loop. + des_ranges = [design_values[k.name] for k in design_map_keys] if change_one_design_at_a_time: factorial_points = generate_snake_zigzag_pattern(*des_ranges) else: @@ -1730,16 +1741,16 @@ def compute_FIM_factorial( failure_count = 0 total_points = len(factorial_points_list) + # save all the FIMs for each point in the factorial design + self.n_parameters = len(model.unknown_parameters) + FIM_all = np.zeros((total_points, self.n_parameters, self.n_parameters)) + time_set = [] curr_point = 1 # Initial current point for design_point in factorial_points_list: - print(f"design_point: {design_point}") # Fix design variables at fixed experimental design point for i in range(len(design_point)): - if design_point[i] is not None: - design_map[i][1].fix(design_point[i]) - else: - pass + design_map_keys[i].fix(design_point[i]) # Timing and logging objects self.logger.info(f"=======Iteration Number: {curr_point} =======") @@ -1750,9 +1761,6 @@ def compute_FIM_factorial( curr_point = success_count + failure_count + 1 self.logger.info(f"This is run {curr_point} out of {total_points}.") self.compute_FIM(model=model, method=method) - print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") - print(model.pprint()) - print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n") success_count += 1 # iteration time iter_t = iter_timer.toc(msg=None) @@ -1783,6 +1791,9 @@ def compute_FIM_factorial( FIM = self._computed_FIM + # Save FIM for the current design point + FIM_all[curr_point - 1, :, :] = FIM + # Compute and record metrics on FIM det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = ( compute_FIM_metrics(FIM) @@ -1806,8 +1817,10 @@ def compute_FIM_factorial( "total_points": total_points, "success_counts": success_count, "failure_counts": failure_count, + "FIM_all": FIM_all.tolist(), # Save all FIMs } ) + self.factorial_results = factorial_results if file_name is not None: diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index e2a018dab14..dff264220f0 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -115,7 +115,7 @@ def generate_snake_zigzag_pattern(*lists): Parameters ---------- - lists: A variable number of 1D list arguments. + *lists: A variable number of 1D list arguments. Yields ------ From 0eb9791ba573139135e01f1f2fa44a1cb3065bef Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Mon, 23 Jun 2025 07:30:46 -0400 Subject: [PATCH 05/25] Changed documentation and improved the error message --- pyomo/contrib/doe/doe.py | 46 +++++++++++++++++++++++--------------- pyomo/contrib/doe/utils.py | 9 ++------ 2 files changed, 30 insertions(+), 25 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 95a0d75acd6..eb57ac0f577 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -45,7 +45,13 @@ from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp import pyomo.environ as pyo -from pyomo.contrib.doe.utils import generate_snake_zigzag_pattern +from pyomo.contrib.doe.utils import ( + generate_snake_zigzag_pattern, +) # , compute_FIM_metrics + +# utils_updated.py is the utils.py from my open PR # 3525. +# When that PR is merged, compute_FIM_metrics will be imported from utils.py and this +# import will be removed. from pyomo.contrib.doe.utils_updated import compute_FIM_metrics from pyomo.opt import SolverStatus @@ -1637,23 +1643,25 @@ def compute_FIM_factorial( model : DoE model, optional The model to perform the full factorial exploration on. Default: None design_values : dict, - dict of lists, of the form {<"var_name">: []}. Default: None. - The `design_values` should have the same key(s) as the `experiment_inputs`. - If one or more design variables are not to be changed, they do not need to - be passed in the `design_values` dictionary, but if they are passed in the - dictionary, then they must be a list of floats. For example, if our - experiment has 4 design variables (i.e., `experiment_inputs`): model.x1, - model.x2, model.x3, and model.x4, their values may be passed as, - design_values= {"x1": [1, 2, 3], "x3": [7], "x4": [-10, 20]}. In this case, - x2 is not be changed and will be fixed at the default value in the model. + dict of lists, of the form {"var_name": [var_values]}. Default: None. + The `design_values` should have the key(s) passed as strings that is a + subset of the `experiment_inputs`. If one or more design variables are not + to be changed, then they should not be passed in the `design_values` + dictionary, but if they are passed in the dictionary, then they must be a + list of floats. For example, if our experiment has 4 design variables + (i.e., `experiment_inputs`): model.x1, model.x2, model.x3, and model.x4, + their values may be passed as, design_values= {"x1": [1, 2, 3], "x3": [7], + "x4": [-10, 20]}. In this case, x2 will not be changed and will be fixed at + the value in the model. method : str, optional string to specify which method should be used. options are ``kaug`` and ``sequential`. Default: "sequential" change_one_design_at_a_time : bool, optional - If True, will generate a snake-like zigzag pattern of the design values that - changes only one of the design variables at a time. If False, will - generate a regular nested for loop that can change from one to all the design - variables at a time. Default: True + If True, will generate a snake-like zigzag combination of the design values + thatchanges only one of the design variables at a time. This combination + may help with the convergence in some scenarios. If False, will + generate a regular nested for loop that can change from one to all the + design variables at a time. Default: True file_name : str, optional if provided, will save the results to a json file. Default: None @@ -1677,8 +1685,8 @@ def compute_FIM_factorial( Raises ------ ValueError - If the design_values' keys is not a subset of the model's experiment_inputs' - keys or if the design_values are not provided as a dictionary of lists. + If the design_values' keys are not a subset of the model's + `experiment_inputs` keys or if the design_values are not provided. """ # Start timer @@ -1703,9 +1711,10 @@ def compute_FIM_factorial( map_keys = set([k.name for k in model.experiment_inputs.keys()]) if not design_keys.issubset(map_keys): incorrect_given_keys = design_keys - map_keys + suggested_keys = map_keys - design_keys raise ValueError( f"design_values keys: {incorrect_given_keys} are incorrect." - f"The key should be from the following: {map_keys}." + f"The keys should be from the following keys: {suggested_keys}." ) # Get the design map keys that match the design_values keys @@ -1823,9 +1832,10 @@ def compute_FIM_factorial( self.factorial_results = factorial_results + # Save the results to a json file based on the file_name provided if file_name is not None: with open(file_name + ".json", "w") as f: - json.dump(self.factorial_results, f, indent=4) # Save results to file + json.dump(self.factorial_results, f, indent=4) self.logger.info(f"Results saved to {file_name}.json.") return self.factorial_results diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index dff264220f0..2a82ad2014f 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -106,16 +106,11 @@ def generate_snake_zigzag_pattern(*lists): """ Generates a multi-dimensional zigzag pattern for an arbitrary number of lists. This pattern is useful for generating patterns for sensitivity analysis when we want - to change one variable at a time. - - This function uses recursion and acts as a generator. The direction of - iteration for each list is determined by the sum of the indices from the - preceding lists. If the sum is even, the list is iterated forwards; - if it's odd, the list is iterated in reverse. + to change one variable at a time. This function uses recursion and acts as a generator. Parameters ---------- - *lists: A variable number of 1D list arguments. + *lists: A variable number of 1D arraylike arguments. Yields ------ From d75e66733ba4d256d50457d9fbe05641a453141a Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Mon, 23 Jun 2025 10:43:58 -0400 Subject: [PATCH 06/25] updated comments --- pyomo/contrib/doe/doe.py | 9 ++++++--- pyomo/contrib/doe/utils.py | 12 ------------ 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index eb57ac0f577..4e1058f7f38 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -49,11 +49,14 @@ generate_snake_zigzag_pattern, ) # , compute_FIM_metrics -# utils_updated.py is the utils.py from my open PR # 3525. -# When that PR is merged, compute_FIM_metrics will be imported from utils.py and this -# import will be removed. from pyomo.contrib.doe.utils_updated import compute_FIM_metrics +""" +utils_updated.py is the utils.py from my open PR # 3525. +When that PR is merged, compute_FIM_metrics will be imported from utils.py and this +import will be removed. +""" + from pyomo.opt import SolverStatus # This small and positive tolerance is used when checking diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 2a82ad2014f..1948346e8d3 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -119,18 +119,6 @@ def generate_snake_zigzag_pattern(*lists): # The main logic is in a nested recursive helper function. def _generate_recursive(depth, index_sum): - """ - Parameters - ---------- - depth : int - corresponds to the index of the current list in the `lists` argument. - Represents which list we are currently processing. - - index_sum : int - It's a running total of the indices chosen from all the previous lists. - The direction of the depth-th list depends entirely on whether index_sum - even or odd. - """ # Base case: If we've processed all lists, we're at the end of a path. if depth == len(lists): yield () From e5b2e73d0701f71ecf6825418ae0826f12fd76c9 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 24 Jun 2025 14:40:56 -0400 Subject: [PATCH 07/25] Added argument checking for generate zigzag pattern. Also added test for patter generator --- pyomo/contrib/doe/doe.py | 9 +- pyomo/contrib/doe/tests/test_utils.py | 140 ++++++++++++++++++++++++++ pyomo/contrib/doe/utils.py | 21 +++- 3 files changed, 161 insertions(+), 9 deletions(-) create mode 100644 pyomo/contrib/doe/tests/test_utils.py diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 4e1058f7f38..cc6a7e4d9e3 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1646,12 +1646,12 @@ def compute_FIM_factorial( model : DoE model, optional The model to perform the full factorial exploration on. Default: None design_values : dict, - dict of lists, of the form {"var_name": [var_values]}. Default: None. + dict of lists or other array-like objects, of the form {"var_name": }. Default: None. The `design_values` should have the key(s) passed as strings that is a subset of the `experiment_inputs`. If one or more design variables are not to be changed, then they should not be passed in the `design_values` dictionary, but if they are passed in the dictionary, then they must be a - list of floats. For example, if our experiment has 4 design variables + array-like object of floats. For example, if our experiment has 4 design variables (i.e., `experiment_inputs`): model.x1, model.x2, model.x3, and model.x4, their values may be passed as, design_values= {"x1": [1, 2, 3], "x3": [7], "x4": [-10, 20]}. In this case, x2 will not be changed and will be fixed at @@ -1705,8 +1705,8 @@ def compute_FIM_factorial( if not design_values: raise ValueError( - "design_values must be provided as a dictionary of lists " - "in the form {<'var_name'>: []}." + "design_values must be provided as a dictionary of array-like objects " + "in the form {<'var_name'>: }." ) # Check whether the design_ranges keys are in the experiment_inputs @@ -1726,6 +1726,7 @@ def compute_FIM_factorial( ] # This ensures that the order of the design_values keys matches the order of the # design_map_keys so that design_point can be constructed correctly in the loop. + # TODO: define an Enum to add different sensitivity analysis sequences des_ranges = [design_values[k.name] for k in design_map_keys] if change_one_design_at_a_time: factorial_points = generate_snake_zigzag_pattern(*des_ranges) diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py new file mode 100644 index 00000000000..3ab34c231ef --- /dev/null +++ b/pyomo/contrib/doe/tests/test_utils.py @@ -0,0 +1,140 @@ +# ___________________________________________________________________________ +# +# 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 generate_snake_zigzag_pattern + + +@unittest.skipIf(not numpy_available, "Numpy is not available") +class TestUtilsFIM(unittest.TestCase): + def test_generate_snake_zigzag_pattern_errors(self): + # Test the error handling with lists + list_2d = [[1, 2, 3], [4, 5, 6]] + with self.assertRaises(ValueError) as cm: + list(generate_snake_zigzag_pattern(list_2d)) + self.assertEqual( + str(cm.exception), + "Argument at position 0 is not 1D. Got shape (2, 3).", + ) + + list_2d_wrong_shape = [[1, 2, 3], [4, 5, 6, 7]] + with self.assertRaises(ValueError) as cm: + list(generate_snake_zigzag_pattern(list_2d_wrong_shape)) + self.assertEqual( + str(cm.exception), + "Argument at position 0 is not 1D array-like.", + ) + + # Test the error handling with tuples + tuple_2d = ((1, 2, 3), (4, 5, 6)) + with self.assertRaises(ValueError) as cm: + list(generate_snake_zigzag_pattern(tuple_2d)) + self.assertEqual( + str(cm.exception), + "Argument at position 0 is not 1D. Got shape (2, 3).", + ) + + tuple_2d_wrong_shape = ((1, 2, 3), (4, 5, 6, 7)) + with self.assertRaises(ValueError) as cm: + list(generate_snake_zigzag_pattern(tuple_2d_wrong_shape)) + self.assertEqual( + str(cm.exception), + "Argument at position 0 is not 1D array-like.", + ) + + # Test the error handling with numpy arrays + array_2d = np.array([[1, 2, 3], [4, 5, 6]]) + with self.assertRaises(ValueError) as cm: + list(generate_snake_zigzag_pattern(array_2d)) + self.assertEqual( + str(cm.exception), + "Argument at position 0 is not 1D. Got shape (2, 3).", + ) + + def test_generate_snake_zigzag_pattern_values(self): + # Test with lists + # Test with a single list + list1 = [1, 2, 3] + result_list1 = list(generate_snake_zigzag_pattern(list1)) + expected_list1 = [(1,), (2,), (3,)] + self.assertEqual(result_list1, expected_list1) + + # Test with two lists + list2 = [4, 5, 6] + result_list2 = list(generate_snake_zigzag_pattern(list1, list2)) + expected_list2 = [ + (1, 4), + (1, 5), + (1, 6), + (2, 6), + (2, 5), + (2, 4), + (3, 4), + (3, 5), + (3, 6), + ] + self.assertEqual(result_list2, expected_list2) + + # Test with three lists + list3 = [7, 8] + result_list3 = list(generate_snake_zigzag_pattern(list1, list2, list3)) + expected_list3 = [ + (1, 4, 7), + (1, 4, 8), + (1, 5, 8), + (1, 5, 7), + (1, 6, 7), + (1, 6, 8), + (2, 6, 8), + (2, 6, 7), + (2, 5, 7), + (2, 5, 8), + (2, 4, 8), + (2, 4, 7), + (3, 4, 7), + (3, 4, 8), + (3, 5, 8), + (3, 5, 7), + (3, 6, 7), + (3, 6, 8), + ] + self.assertEqual(result_list3, expected_list3) + + # Test with tuples + tuple1 = (1, 2, 3) + result_tuple1 = list(generate_snake_zigzag_pattern(tuple1)) + tuple2 = (4, 5, 6) + result_tuple2 = list(generate_snake_zigzag_pattern(tuple1, tuple2)) + tuple3 = (7, 8) + result_tuple3 = list(generate_snake_zigzag_pattern(tuple1, tuple2, tuple3)) + self.assertEqual(result_tuple1, expected_list1) + self.assertEqual(result_tuple2, expected_list2) + self.assertEqual(result_tuple3, expected_list3) + + # Test with numpy arrays + array1 = np.array([1, 2, 3]) + array2 = np.array([4, 5, 6]) + array3 = np.array([7, 8]) + result_array1 = list(generate_snake_zigzag_pattern(array1)) + result_array2 = list(generate_snake_zigzag_pattern(array1, array2)) + result_array3 = list(generate_snake_zigzag_pattern(array1, array2, array3)) + self.assertEqual(result_array1, expected_list1) + self.assertEqual(result_array2, expected_list2) + self.assertEqual(result_array3, expected_list3) + + # Test with mixed types(List, Tuple, numpy array) + result_mixed = list(generate_snake_zigzag_pattern(list1, tuple2, array3)) + self.assertEqual(result_mixed, expected_list3) + + +if __name__ == "__main__": + unittest.main() diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 1948346e8d3..c893c458f98 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -102,29 +102,40 @@ def rescale_FIM(FIM, param_vals): # return param_list -def generate_snake_zigzag_pattern(*lists): +def generate_snake_zigzag_pattern(*array_like_args): """ - Generates a multi-dimensional zigzag pattern for an arbitrary number of lists. + Generates a multi-dimensional zigzag pattern for an arbitrary number of 1D array-like arguments. This pattern is useful for generating patterns for sensitivity analysis when we want to change one variable at a time. This function uses recursion and acts as a generator. Parameters ---------- - *lists: A variable number of 1D arraylike arguments. + *array_like_args: A variable number of 1D array-like arguments. Yields ------ A tuple representing points in the snake-like zigzag pattern. """ + # throw an error if the array_like_args are not array-like or all 1D + for i, arg in enumerate(array_like_args): + try: + arr_np = np.asarray(arg) + except Exception: + raise ValueError(f"Argument at position {i} is not 1D array-like.") + + if arr_np.ndim != 1: + raise ValueError( + f"Argument at position {i} is not 1D. Got shape {arr_np.shape}." + ) # The main logic is in a nested recursive helper function. def _generate_recursive(depth, index_sum): # Base case: If we've processed all lists, we're at the end of a path. - if depth == len(lists): + if depth == len(array_like_args): yield () return - current_list = lists[depth] + current_list = array_like_args[depth] # Determine the iteration direction based on the sum of parent indices. # This is the mathematical rule for the zigzag. From 13848fa1188536a811d1887996b4d0da598b417e Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Wed, 25 Jun 2025 12:47:15 -0400 Subject: [PATCH 08/25] changed the name of the pattern generator --- pyomo/contrib/doe/tests/test_utils.py | 38 +++++++++++++-------------- pyomo/contrib/doe/utils.py | 4 +-- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py index 3ab34c231ef..446eeeb4ddc 100644 --- a/pyomo/contrib/doe/tests/test_utils.py +++ b/pyomo/contrib/doe/tests/test_utils.py @@ -11,16 +11,16 @@ from pyomo.common.dependencies import numpy as np, numpy_available import pyomo.common.unittest as unittest -from pyomo.contrib.doe.utils import generate_snake_zigzag_pattern +from pyomo.contrib.doe.utils import serpentine_traversal_sampling @unittest.skipIf(not numpy_available, "Numpy is not available") class TestUtilsFIM(unittest.TestCase): - def test_generate_snake_zigzag_pattern_errors(self): + def test_serpentine_traversal_sampling_errors(self): # Test the error handling with lists - list_2d = [[1, 2, 3], [4, 5, 6]] + list_2d_bad = [[1, 2, 3], [4, 5, 6]] with self.assertRaises(ValueError) as cm: - list(generate_snake_zigzag_pattern(list_2d)) + list(serpentine_traversal_sampling(list_2d_bad)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3).", @@ -28,7 +28,7 @@ def test_generate_snake_zigzag_pattern_errors(self): list_2d_wrong_shape = [[1, 2, 3], [4, 5, 6, 7]] with self.assertRaises(ValueError) as cm: - list(generate_snake_zigzag_pattern(list_2d_wrong_shape)) + list(serpentine_traversal_sampling(list_2d_wrong_shape)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D array-like.", @@ -37,7 +37,7 @@ def test_generate_snake_zigzag_pattern_errors(self): # Test the error handling with tuples tuple_2d = ((1, 2, 3), (4, 5, 6)) with self.assertRaises(ValueError) as cm: - list(generate_snake_zigzag_pattern(tuple_2d)) + list(serpentine_traversal_sampling(tuple_2d)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3).", @@ -45,7 +45,7 @@ def test_generate_snake_zigzag_pattern_errors(self): tuple_2d_wrong_shape = ((1, 2, 3), (4, 5, 6, 7)) with self.assertRaises(ValueError) as cm: - list(generate_snake_zigzag_pattern(tuple_2d_wrong_shape)) + list(serpentine_traversal_sampling(tuple_2d_wrong_shape)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D array-like.", @@ -54,23 +54,23 @@ def test_generate_snake_zigzag_pattern_errors(self): # Test the error handling with numpy arrays array_2d = np.array([[1, 2, 3], [4, 5, 6]]) with self.assertRaises(ValueError) as cm: - list(generate_snake_zigzag_pattern(array_2d)) + list(serpentine_traversal_sampling(array_2d)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3).", ) - def test_generate_snake_zigzag_pattern_values(self): + def test_serpentine_traversal_sampling_values(self): # Test with lists # Test with a single list list1 = [1, 2, 3] - result_list1 = list(generate_snake_zigzag_pattern(list1)) + result_list1 = list(serpentine_traversal_sampling(list1)) expected_list1 = [(1,), (2,), (3,)] self.assertEqual(result_list1, expected_list1) # Test with two lists list2 = [4, 5, 6] - result_list2 = list(generate_snake_zigzag_pattern(list1, list2)) + result_list2 = list(serpentine_traversal_sampling(list1, list2)) expected_list2 = [ (1, 4), (1, 5), @@ -86,7 +86,7 @@ def test_generate_snake_zigzag_pattern_values(self): # Test with three lists list3 = [7, 8] - result_list3 = list(generate_snake_zigzag_pattern(list1, list2, list3)) + result_list3 = list(serpentine_traversal_sampling(list1, list2, list3)) expected_list3 = [ (1, 4, 7), (1, 4, 8), @@ -111,11 +111,11 @@ def test_generate_snake_zigzag_pattern_values(self): # Test with tuples tuple1 = (1, 2, 3) - result_tuple1 = list(generate_snake_zigzag_pattern(tuple1)) + result_tuple1 = list(serpentine_traversal_sampling(tuple1)) tuple2 = (4, 5, 6) - result_tuple2 = list(generate_snake_zigzag_pattern(tuple1, tuple2)) + result_tuple2 = list(serpentine_traversal_sampling(tuple1, tuple2)) tuple3 = (7, 8) - result_tuple3 = list(generate_snake_zigzag_pattern(tuple1, tuple2, tuple3)) + result_tuple3 = list(serpentine_traversal_sampling(tuple1, tuple2, tuple3)) self.assertEqual(result_tuple1, expected_list1) self.assertEqual(result_tuple2, expected_list2) self.assertEqual(result_tuple3, expected_list3) @@ -124,15 +124,15 @@ def test_generate_snake_zigzag_pattern_values(self): array1 = np.array([1, 2, 3]) array2 = np.array([4, 5, 6]) array3 = np.array([7, 8]) - result_array1 = list(generate_snake_zigzag_pattern(array1)) - result_array2 = list(generate_snake_zigzag_pattern(array1, array2)) - result_array3 = list(generate_snake_zigzag_pattern(array1, array2, array3)) + result_array1 = list(serpentine_traversal_sampling(array1)) + result_array2 = list(serpentine_traversal_sampling(array1, array2)) + result_array3 = list(serpentine_traversal_sampling(array1, array2, array3)) self.assertEqual(result_array1, expected_list1) self.assertEqual(result_array2, expected_list2) self.assertEqual(result_array3, expected_list3) # Test with mixed types(List, Tuple, numpy array) - result_mixed = list(generate_snake_zigzag_pattern(list1, tuple2, array3)) + result_mixed = list(serpentine_traversal_sampling(list1, tuple2, array3)) self.assertEqual(result_mixed, expected_list3) diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index c893c458f98..b42cfcd4b00 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -102,9 +102,9 @@ def rescale_FIM(FIM, param_vals): # return param_list -def generate_snake_zigzag_pattern(*array_like_args): +def serpentine_traversal_sampling(*array_like_args): """ - Generates a multi-dimensional zigzag pattern for an arbitrary number of 1D array-like arguments. + Generates a multi-dimensional pattern for an arbitrary number of 1D array-like arguments. This pattern is useful for generating patterns for sensitivity analysis when we want to change one variable at a time. This function uses recursion and acts as a generator. From 0c87484adbad9419901aaf802e2dd306656c405d Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Fri, 27 Jun 2025 10:38:03 -0400 Subject: [PATCH 09/25] Changed the name of the grid sampling function and deleted the utils_updated.py file --- pyomo/contrib/doe/doe.py | 13 +- pyomo/contrib/doe/examples/reactor_example.py | 2 - pyomo/contrib/doe/tests/test_utils.py | 64 ++--- pyomo/contrib/doe/utils.py | 7 +- pyomo/contrib/doe/utils_updated.py | 267 ------------------ 5 files changed, 28 insertions(+), 325 deletions(-) delete mode 100644 pyomo/contrib/doe/utils_updated.py diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 23c7da4c9d1..79bad93dfdf 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -49,18 +49,9 @@ check_FIM, compute_FIM_metrics, _SMALL_TOLERANCE_DEFINITENESS, + snake_traversal_grid_sampling, ) -from pyomo.contrib.doe.utils import ( - generate_snake_zigzag_pattern, -) # , compute_FIM_metrics - -from pyomo.contrib.doe.utils_updated import compute_FIM_metrics -""" -utils_updated.py is the utils.py from my open PR # 3525. -When that PR is merged, compute_FIM_metrics will be imported from utils.py and this -import will be removed. -""" from pyomo.opt import SolverStatus @@ -1722,7 +1713,7 @@ def compute_FIM_factorial( # TODO: define an Enum to add different sensitivity analysis sequences des_ranges = [design_values[k.name] for k in design_map_keys] if change_one_design_at_a_time: - factorial_points = generate_snake_zigzag_pattern(*des_ranges) + factorial_points = snake_traversal_grid_sampling(*des_ranges) else: factorial_points = product(*des_ranges) diff --git a/pyomo/contrib/doe/examples/reactor_example.py b/pyomo/contrib/doe/examples/reactor_example.py index 88635b5c80e..0038ed114de 100644 --- a/pyomo/contrib/doe/examples/reactor_example.py +++ b/pyomo/contrib/doe/examples/reactor_example.py @@ -14,8 +14,6 @@ from pyomo.contrib.doe import DesignOfExperiments import pyomo.environ as pyo -import pyomo.environ as pyo -import idaes.core.solvers.get_solver import json diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py index b0981609d46..eb3956176c8 100644 --- a/pyomo/contrib/doe/tests/test_utils.py +++ b/pyomo/contrib/doe/tests/test_utils.py @@ -15,6 +15,7 @@ check_FIM, compute_FIM_metrics, get_FIM_metrics, + snake_traversal_grid_sampling, _SMALL_TOLERANCE_DEFINITENESS, _SMALL_TOLERANCE_SYMMETRY, _SMALL_TOLERANCE_IMG, @@ -141,82 +142,61 @@ def test_get_FIM_metrics(self): self.assertEqual(fim_metrics["log10(E-Optimality)"], E_opt_expected) self.assertEqual(fim_metrics["log10(Modified E-Optimality)"], ME_opt_expected) - -if __name__ == "__main__": - unittest.main() -# ___________________________________________________________________________ -# -# 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 serpentine_traversal_sampling - - -@unittest.skipIf(not numpy_available, "Numpy is not available") -class TestUtilsFIM(unittest.TestCase): - def test_serpentine_traversal_sampling_errors(self): + def test_snake_traversal_grid_sampling_errors(self): # Test the error handling with lists list_2d_bad = [[1, 2, 3], [4, 5, 6]] with self.assertRaises(ValueError) as cm: - list(serpentine_traversal_sampling(list_2d_bad)) + list(snake_traversal_grid_sampling(list_2d_bad)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3).", ) - list_2d_wrong_shape = [[1, 2, 3], [4, 5, 6, 7]] + list_2d_wrong_shape_bad = [[1, 2, 3], [4, 5, 6, 7]] with self.assertRaises(ValueError) as cm: - list(serpentine_traversal_sampling(list_2d_wrong_shape)) + list(snake_traversal_grid_sampling(list_2d_wrong_shape_bad)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D array-like.", ) # Test the error handling with tuples - tuple_2d = ((1, 2, 3), (4, 5, 6)) + tuple_2d_bad = ((1, 2, 3), (4, 5, 6)) with self.assertRaises(ValueError) as cm: - list(serpentine_traversal_sampling(tuple_2d)) + list(snake_traversal_grid_sampling(tuple_2d_bad)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3).", ) - tuple_2d_wrong_shape = ((1, 2, 3), (4, 5, 6, 7)) + tuple_2d_wrong_shape_bad = ((1, 2, 3), (4, 5, 6, 7)) with self.assertRaises(ValueError) as cm: - list(serpentine_traversal_sampling(tuple_2d_wrong_shape)) + list(snake_traversal_grid_sampling(tuple_2d_wrong_shape_bad)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D array-like.", ) # Test the error handling with numpy arrays - array_2d = np.array([[1, 2, 3], [4, 5, 6]]) + array_2d_bad = np.array([[1, 2, 3], [4, 5, 6]]) with self.assertRaises(ValueError) as cm: - list(serpentine_traversal_sampling(array_2d)) + list(snake_traversal_grid_sampling(array_2d_bad)) self.assertEqual( str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3).", ) - def test_serpentine_traversal_sampling_values(self): + def test_snake_traversal_grid_sampling_values(self): # Test with lists # Test with a single list list1 = [1, 2, 3] - result_list1 = list(serpentine_traversal_sampling(list1)) + result_list1 = list(snake_traversal_grid_sampling(list1)) expected_list1 = [(1,), (2,), (3,)] self.assertEqual(result_list1, expected_list1) # Test with two lists list2 = [4, 5, 6] - result_list2 = list(serpentine_traversal_sampling(list1, list2)) + result_list2 = list(snake_traversal_grid_sampling(list1, list2)) expected_list2 = [ (1, 4), (1, 5), @@ -232,7 +212,7 @@ def test_serpentine_traversal_sampling_values(self): # Test with three lists list3 = [7, 8] - result_list3 = list(serpentine_traversal_sampling(list1, list2, list3)) + result_list3 = list(snake_traversal_grid_sampling(list1, list2, list3)) expected_list3 = [ (1, 4, 7), (1, 4, 8), @@ -257,11 +237,11 @@ def test_serpentine_traversal_sampling_values(self): # Test with tuples tuple1 = (1, 2, 3) - result_tuple1 = list(serpentine_traversal_sampling(tuple1)) + result_tuple1 = list(snake_traversal_grid_sampling(tuple1)) tuple2 = (4, 5, 6) - result_tuple2 = list(serpentine_traversal_sampling(tuple1, tuple2)) + result_tuple2 = list(snake_traversal_grid_sampling(tuple1, tuple2)) tuple3 = (7, 8) - result_tuple3 = list(serpentine_traversal_sampling(tuple1, tuple2, tuple3)) + result_tuple3 = list(snake_traversal_grid_sampling(tuple1, tuple2, tuple3)) self.assertEqual(result_tuple1, expected_list1) self.assertEqual(result_tuple2, expected_list2) self.assertEqual(result_tuple3, expected_list3) @@ -270,15 +250,15 @@ def test_serpentine_traversal_sampling_values(self): array1 = np.array([1, 2, 3]) array2 = np.array([4, 5, 6]) array3 = np.array([7, 8]) - result_array1 = list(serpentine_traversal_sampling(array1)) - result_array2 = list(serpentine_traversal_sampling(array1, array2)) - result_array3 = list(serpentine_traversal_sampling(array1, array2, array3)) + result_array1 = list(snake_traversal_grid_sampling(array1)) + result_array2 = list(snake_traversal_grid_sampling(array1, array2)) + result_array3 = list(snake_traversal_grid_sampling(array1, array2, array3)) self.assertEqual(result_array1, expected_list1) self.assertEqual(result_array2, expected_list2) self.assertEqual(result_array3, expected_list3) # Test with mixed types(List, Tuple, numpy array) - result_mixed = list(serpentine_traversal_sampling(list1, tuple2, array3)) + result_mixed = list(snake_traversal_grid_sampling(list1, tuple2, array3)) self.assertEqual(result_mixed, expected_list3) diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 332c7fe6a48..2a3465a04ef 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -267,7 +267,7 @@ def get_FIM_metrics(FIM): } -def serpentine_traversal_sampling(*array_like_args): +def snake_traversal_grid_sampling(*array_like_args): """ Generates a multi-dimensional pattern for an arbitrary number of 1D array-like arguments. This pattern is useful for generating patterns for sensitivity analysis when we want @@ -275,11 +275,12 @@ def serpentine_traversal_sampling(*array_like_args): Parameters ---------- - *array_like_args: A variable number of 1D array-like arguments. + *array_like_args: A variable number of 1D array-like objects as arguments. Yields ------ - A tuple representing points in the snake-like zigzag pattern. + A tuple representing points in the snake pattern. + Naming source of the pattern: https://dev.to/thesanjeevsharma/dsa-101-matrix-30lf """ # throw an error if the array_like_args are not array-like or all 1D for i, arg in enumerate(array_like_args): diff --git a/pyomo/contrib/doe/utils_updated.py b/pyomo/contrib/doe/utils_updated.py deleted file mode 100644 index 2232fd88c52..00000000000 --- a/pyomo/contrib/doe/utils_updated.py +++ /dev/null @@ -1,267 +0,0 @@ -# ___________________________________________________________________________ -# -# 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. -# -# Pyomo.DoE was produced under the Department of Energy Carbon Capture Simulation -# Initiative (CCSI), and is copyright (c) 2022 by the software owners: -# TRIAD National Security, LLC., Lawrence Livermore National Security, LLC., -# Lawrence Berkeley National Laboratory, Pacific Northwest National Laboratory, -# Battelle Memorial Institute, University of Notre Dame, -# The University of Pittsburgh, The University of Texas at Austin, -# University of Toledo, West Virginia University, et al. All rights reserved. -# -# NOTICE. This Software was developed under funding from the -# U.S. Department of Energy and the U.S. Government consequently retains -# certain rights. As such, the U.S. Government has been granted for itself -# and others acting on its behalf a paid-up, nonexclusive, irrevocable, -# worldwide license in the Software to reproduce, distribute copies to the -# public, prepare derivative works, and perform publicly and display -# publicly, and to permit other to do so. -# ___________________________________________________________________________ - -import pyomo.environ as pyo - -from pyomo.common.dependencies import numpy as np, numpy_available - -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): - """ - Rescales the FIM based on the input and parameter vals. - It is assumed the parameter vals align with the FIM - dimensions such that (1, i) corresponds to the i-th - column or row of the FIM. - - Parameters - ---------- - FIM: 2D numpy array to be scaled - param_vals: scaling factors for the parameters - - """ - if isinstance(param_vals, list): - param_vals = np.array([param_vals]) - elif isinstance(param_vals, np.ndarray): - if len(param_vals.shape) > 2 or ( - (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) - ) - if len(param_vals.shape) == 1: - param_vals = np.array([param_vals]) - else: - raise ValueError( - "param_vals should be a list or numpy array of dimensions: 1 by `n_params`" - ) - scaling_mat = (1 / param_vals).transpose().dot((1 / param_vals)) - scaled_FIM = np.multiply(FIM, scaling_mat) - return scaled_FIM - - # TODO: Add swapping parameters for variables helper function - # def get_parameters_from_suffix(suffix, fix_vars=False): - # """ - # Finds the Params within the suffix provided. It will also check to see - # if there are Vars in the suffix provided. ``fix_vars`` will indicate - # if we should fix all the Vars in the set or not. - # - # Parameters - # ---------- - # suffix: pyomo Suffix object, contains the components to be checked - # as keys - # fix_vars: boolean, whether or not to fix the Vars, default = False - # - # Returns - # ------- - # param_list: list of Param - # """ - # param_list = [] - # - # # FIX THE MODEL TREE ISSUE WHERE I GET base_model. INSTEAD OF - # # Check keys if they are Param or Var. Fix the vars if ``fix_vars`` is True - # for k, v in suffix.items(): - # if isinstance(k, ParamData): - # param_list.append(k.name) - # elif isinstance(k, VarData): - # if fix_vars: - # k.fix() - # else: - # 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 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, - } From 6db9e6666e8c32504a3439113f6fd74be6866581 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Mon, 30 Jun 2025 14:27:51 -0400 Subject: [PATCH 10/25] added rel_change and abs_change for change in design --- pyomo/contrib/doe/doe.py | 89 +++++++++++++++++++++++++++----------- pyomo/contrib/doe/utils.py | 63 ++++++++++++++------------- 2 files changed, 96 insertions(+), 56 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 79bad93dfdf..0ce09a6722c 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -50,6 +50,7 @@ compute_FIM_metrics, _SMALL_TOLERANCE_DEFINITENESS, snake_traversal_grid_sampling, + update_model_from_suffix, ) @@ -1615,7 +1616,8 @@ def compute_FIM_full_factorial( def compute_FIM_factorial( self, model=None, - design_values: dict = None, + abs_change: list = None, + rel_change: list = None, method="sequential", change_one_design_at_a_time=True, file_name: str = None, @@ -1629,6 +1631,14 @@ def compute_FIM_factorial( ---------- model : DoE model, optional The model to perform the full factorial exploration on. Default: None + # TODO: Update doc string for absolute and relative change + abs_change : list, optional + Absolute change in the design variable values. Default: None. + If provided, will use this value to generate the design values. + If not provided, will use the `design_values` parameter. + rel_change : list, optional + Relative change in the design variable values. Default: None. + If provided, will use this value to generate the design values. design_values : dict, dict of lists or other array-like objects, of the form {"var_name": }. Default: None. The `design_values` should have the key(s) passed as strings that is a @@ -1687,35 +1697,47 @@ def compute_FIM_factorial( ).clone() model = self.factorial_model - if not design_values: - raise ValueError( - "design_values must be provided as a dictionary of array-like objects " - "in the form {<'var_name'>: }." - ) - - # Check whether the design_ranges keys are in the experiment_inputs - design_keys = set(design_values.keys()) - map_keys = set([k.name for k in model.experiment_inputs.keys()]) - if not design_keys.issubset(map_keys): - incorrect_given_keys = design_keys - map_keys - suggested_keys = map_keys - design_keys - raise ValueError( - f"design_values keys: {incorrect_given_keys} are incorrect." - f"The keys should be from the following keys: {suggested_keys}." - ) - # Get the design map keys that match the design_values keys - design_map_keys = [ - k for k in model.experiment_inputs.keys() if k.name in design_values.keys() - ] + # design_map_keys = [ + # k for k in model.experiment_inputs.keys() if k.name in design_values.keys() + # ] # This ensures that the order of the design_values keys matches the order of the # design_map_keys so that design_point can be constructed correctly in the loop. # TODO: define an Enum to add different sensitivity analysis sequences - des_ranges = [design_values[k.name] for k in design_map_keys] + # des_ranges = [design_values[k.name] for k in design_map_keys] + + design_keys = [k for k in model.experiment_inputs.keys()] + + design_values = [] + for i, comp in enumerate(design_keys): + lb = comp.lb + ub = comp.ub + if lb is None or ub is None: + raise ValueError(f"{comp.name} does not have a lower or upper bound.") + + if abs_change[i] is None and rel_change[i] is None: + n_des = 5 # Default number of points in the design value + des_val = np.linspace(lb, ub, n_des) + + elif abs_change[i] is not None and rel_change[i] is not None: + des_val = [] + del_val = comp.lb * rel_change[i] + abs_change[i] + if del_val == 0: + raise ValueError( + f"Design variable {comp.name} has no change in value - check " + "abs_change and rel_change values." + ) + val = lb + while val <= ub: + des_val.append(val) + val += del_val + + design_values.append(des_val) + if change_one_design_at_a_time: - factorial_points = snake_traversal_grid_sampling(*des_ranges) + factorial_points = snake_traversal_grid_sampling(*design_values) else: - factorial_points = product(*des_ranges) + factorial_points = product(*design_values) factorial_points_list = list(factorial_points) @@ -1747,7 +1769,7 @@ def compute_FIM_factorial( for design_point in factorial_points_list: # Fix design variables at fixed experimental design point for i in range(len(design_point)): - design_map_keys[i].fix(design_point[i]) + design_keys[i].fix(design_point[i]) # Timing and logging objects self.logger.info(f"=======Iteration Number: {curr_point} =======") @@ -1817,6 +1839,23 @@ def compute_FIM_factorial( "FIM_all": FIM_all.tolist(), # Save all FIMs } ) + if self.tee: + exclude_keys = { + "total_points", + "success_counts", + "failure_counts", + "FIM_all", + } + dict_for_df = { + k: v for k, v in factorial_results.items() if k not in exclude_keys + } + res_df = pd.DataFrame(dict_for_df) + print("\n\n=========Factorial results DataFrame===========") + print(res_df) + print("\n\n") + print("Total points:", total_points) + print("Success counts:", success_count) + print("Failure counts:", failure_count) self.factorial_results = factorial_results diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 356444eca80..bb1fcf881f7 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -91,37 +91,38 @@ def rescale_FIM(FIM, param_vals): scaled_FIM = np.multiply(FIM, scaling_mat) return scaled_FIM - # TODO: Add swapping parameters for variables helper function - # def get_parameters_from_suffix(suffix, fix_vars=False): - # """ - # Finds the Params within the suffix provided. It will also check to see - # if there are Vars in the suffix provided. ``fix_vars`` will indicate - # if we should fix all the Vars in the set or not. - # - # Parameters - # ---------- - # suffix: pyomo Suffix object, contains the components to be checked - # as keys - # fix_vars: boolean, whether or not to fix the Vars, default = False - # - # Returns - # ------- - # param_list: list of Param - # """ - # param_list = [] - # - # # FIX THE MODEL TREE ISSUE WHERE I GET base_model. INSTEAD OF - # # Check keys if they are Param or Var. Fix the vars if ``fix_vars`` is True - # for k, v in suffix.items(): - # if isinstance(k, ParamData): - # param_list.append(k.name) - # elif isinstance(k, VarData): - # if fix_vars: - # k.fix() - # else: - # pass # ToDo: Write error for suffix keys that aren't ParamData or VarData - # - # return param_list + +# TODO: Add swapping parameters for variables helper function +# def get_parameters_from_suffix(suffix, fix_vars=False): +# """ +# Finds the Params within the suffix provided. It will also check to see +# if there are Vars in the suffix provided. ``fix_vars`` will indicate +# if we should fix all the Vars in the set or not. +# +# Parameters +# ---------- +# suffix: pyomo Suffix object, contains the components to be checked +# as keys +# fix_vars: boolean, whether or not to fix the Vars, default = False +# +# Returns +# ------- +# param_list: list of Param +# """ +# param_list = [] +# +# # FIX THE MODEL TREE ISSUE WHERE I GET base_model. INSTEAD OF +# # Check keys if they are Param or Var. Fix the vars if ``fix_vars`` is True +# for k, v in suffix.items(): +# if isinstance(k, ParamData): +# param_list.append(k.name) +# elif isinstance(k, VarData): +# if fix_vars: +# k.fix() +# else: +# pass # ToDo: Write error for suffix keys that aren't ParamData or VarData +# +# return param_list # Adding utility to update parameter values in a model based on the suffix From f12d9083391f93180ade153d2bd1b470100eb33d Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 1 Jul 2025 09:03:05 -0400 Subject: [PATCH 11/25] added update_suffix function in compute the sensitivity --- pyomo/contrib/doe/doe.py | 92 ++++++++++++++++++++++++++-------------- 1 file changed, 60 insertions(+), 32 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 0ce09a6722c..e3f8b9b8956 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -70,6 +70,11 @@ class FiniteDifferenceStep(Enum): backward = "backward" +class SensitivityInitialization(Enum): + snake_traversal = "snake_traversal" + nested_for_loop = "nested_for_loop" + + class DesignOfExperiments: def __init__( self, @@ -1618,8 +1623,9 @@ def compute_FIM_factorial( model=None, abs_change: list = None, rel_change: list = None, + n_designs: int = 5, method="sequential", - change_one_design_at_a_time=True, + initialization_scheme=SensitivityInitialization.snake_traversal, file_name: str = None, ): """Will run a simulation-based factorial exploration of the experimental input @@ -1635,21 +1641,11 @@ def compute_FIM_factorial( abs_change : list, optional Absolute change in the design variable values. Default: None. If provided, will use this value to generate the design values. - If not provided, will use the `design_values` parameter. + If not provided, will use the n_designs to generate the design values. + If `abs_change` is provided, `rel_change` must also be provided. rel_change : list, optional Relative change in the design variable values. Default: None. If provided, will use this value to generate the design values. - design_values : dict, - dict of lists or other array-like objects, of the form {"var_name": }. Default: None. - The `design_values` should have the key(s) passed as strings that is a - subset of the `experiment_inputs`. If one or more design variables are not - to be changed, then they should not be passed in the `design_values` - dictionary, but if they are passed in the dictionary, then they must be a - array-like object of floats. For example, if our experiment has 4 design variables - (i.e., `experiment_inputs`): model.x1, model.x2, model.x3, and model.x4, - their values may be passed as, design_values= {"x1": [1, 2, 3], "x3": [7], - "x4": [-10, 20]}. In this case, x2 will not be changed and will be fixed at - the value in the model. method : str, optional string to specify which method should be used. options are ``kaug`` and ``sequential`. Default: "sequential" @@ -1697,29 +1693,44 @@ def compute_FIM_factorial( ).clone() model = self.factorial_model - # Get the design map keys that match the design_values keys - # design_map_keys = [ - # k for k in model.experiment_inputs.keys() if k.name in design_values.keys() - # ] - # This ensures that the order of the design_values keys matches the order of the - # design_map_keys so that design_point can be constructed correctly in the loop. - # TODO: define an Enum to add different sensitivity analysis sequences - # des_ranges = [design_values[k.name] for k in design_map_keys] - design_keys = [k for k in model.experiment_inputs.keys()] + # check if abs_change and rel_change are provided and have the same length as + # design_keys + if abs_change: + if len(abs_change) != len(design_keys): + raise ValueError( + "`abs_change` must have the same length of " + f"`{len(design_keys)}` as `design_keys`." + ) + if rel_change: + if len(rel_change) != len(design_keys): + raise ValueError( + "`rel_change` must have the same length of " + f"`{len(design_keys)}` as `design_keys`." + ) + + # if either abs_change or rel_change is not provided, set it to list of + # zeros + if abs_change or rel_change: + if not abs_change: + abs_change = [0] * len(design_keys) + elif not rel_change: + rel_change = [0] * len(design_keys) + design_values = [] + # loop over design keys and generate design values for i, comp in enumerate(design_keys): lb = comp.lb ub = comp.ub + # Check if the component has finite lower and upper bounds if lb is None or ub is None: raise ValueError(f"{comp.name} does not have a lower or upper bound.") - if abs_change[i] is None and rel_change[i] is None: - n_des = 5 # Default number of points in the design value - des_val = np.linspace(lb, ub, n_des) + if abs_change is None and rel_change is None: + des_val = np.linspace(lb, ub, n_designs) - elif abs_change[i] is not None and rel_change[i] is not None: + elif abs_change or rel_change: des_val = [] del_val = comp.lb * rel_change[i] + abs_change[i] if del_val == 0: @@ -1732,12 +1743,25 @@ def compute_FIM_factorial( des_val.append(val) val += del_val + else: + raise ValueError( + "Unexpected error in generating design values. Please check the " + "input parameters." + ) + design_values.append(des_val) - if change_one_design_at_a_time: + # generate the factorial points based on the initialization scheme + if initialization_scheme == SensitivityInitialization.snake_traversal: factorial_points = snake_traversal_grid_sampling(*design_values) - else: + elif initialization_scheme == SensitivityInitialization.nested_for_loop: factorial_points = product(*design_values) + else: + self.logger.warning( + "initialization_scheme not recognized. Using `snake_traversal` as " + "default." + ) + factorial_points = snake_traversal_grid_sampling(*design_values) factorial_points_list = list(factorial_points) @@ -1760,7 +1784,7 @@ def compute_FIM_factorial( failure_count = 0 total_points = len(factorial_points_list) - # save all the FIMs for each point in the factorial design + # save the FIM for each point in the factorial design self.n_parameters = len(model.unknown_parameters) FIM_all = np.zeros((total_points, self.n_parameters, self.n_parameters)) @@ -1771,6 +1795,10 @@ def compute_FIM_factorial( for i in range(len(design_point)): design_keys[i].fix(design_point[i]) + # TODO: check the following lines of code instead of the for loop above + # update_model_from_suffix( + # model, 'experiment_inputs', design_point) + # Timing and logging objects self.logger.info(f"=======Iteration Number: {curr_point} =======") iter_timer = TicTocTimer() @@ -1850,12 +1878,12 @@ def compute_FIM_factorial( k: v for k, v in factorial_results.items() if k not in exclude_keys } res_df = pd.DataFrame(dict_for_df) - print("\n\n=========Factorial results DataFrame===========") - print(res_df) - print("\n\n") + print("\n\n=========Factorial results===========") print("Total points:", total_points) print("Success counts:", success_count) print("Failure counts:", failure_count) + print("\n\n") + print(res_df) self.factorial_results = factorial_results From fda812e5ac5b017da41ae97bd188410d5c0f7829 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 1 Jul 2025 14:02:10 -0400 Subject: [PATCH 12/25] working code. added design vals and also changed the settings of pd df. The code is run by using a simple model. --- pyomo/contrib/doe/doe.py | 236 +++++++++++++++++++++++---------------- 1 file changed, 140 insertions(+), 96 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index e3f8b9b8956..386d12b6813 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1621,11 +1621,13 @@ def compute_FIM_full_factorial( def compute_FIM_factorial( self, model=None, + design_vals: list = None, abs_change: list = None, rel_change: list = None, n_designs: int = 5, method="sequential", - initialization_scheme=SensitivityInitialization.snake_traversal, + df_settings=(True, None, None, 500), + initialization_scheme="snake_traversal", file_name: str = None, ): """Will run a simulation-based factorial exploration of the experimental input @@ -1637,24 +1639,41 @@ def compute_FIM_factorial( ---------- model : DoE model, optional The model to perform the full factorial exploration on. Default: None - # TODO: Update doc string for absolute and relative change + design_vals : list, optional + A list of design values to use for the full factorial exploration. + Default: None abs_change : list, optional Absolute change in the design variable values. Default: None. If provided, will use this value to generate the design values. - If not provided, will use the n_designs to generate the design values. - If `abs_change` is provided, `rel_change` must also be provided. + If `abs_change` is provided, but `rel_change` is not provided, `rel_change` + will be set to zero. + Formula to calculate the design values: + change_in_value = lower_bound * rel_change + abs_change` + design_value += design_value + change_in_value rel_change : list, optional Relative change in the design variable values. Default: None. If provided, will use this value to generate the design values. + If `rel_change` is provided, but `abs_change` is not provided, `abs_change` + will be set to zero. + n_designs : int, optional + Number of designs to generate for each design variable. Default: 5. + If `abs_change` and/or `rel_change` are provided, this value will be ignored. method : str, optional string to specify which method should be used. options are ``kaug`` and ``sequential`. Default: "sequential" - change_one_design_at_a_time : bool, optional - If True, will generate a snake-like zigzag combination of the design values - thatchanges only one of the design variables at a time. This combination - may help with the convergence in some scenarios. If False, will - generate a regular nested for loop that can change from one to all the - design variables at a time. Default: True + df_settings : tuple, optional + A tuple containing the settings for set_option() method of the pandas + DataFrame. Default: (True, None, None, 500) + - first element: whether to return a pandas DataFrame (True/False) + - second element: number of max_columns for the DataFrame. Default: None, + i.e., no limit on the number of columns. + - third element: number of max_rows for the DataFrame. Default: None, + i.e., no limit on the number of rows. + - fourth element: display width for the DataFrame. Default: 500. + initialization_scheme : str, optional + Which scheme to use for initializing the design variables. + Options are ``"snake_traversal"`` and ``"nested_for_loop"``. + Default: "snake_traversal" file_name : str, optional if provided, will save the results to a json file. Default: None @@ -1668,18 +1687,15 @@ def compute_FIM_factorial( - "log10(A-opt)": list of A-optimality values - "log10(E-opt)": list of E-optimality values - "log10(ME-opt)": list of ME-optimality values + - "eigval_min": list of minimum eigenvalues + - "eigval_max": list of maximum eigenvalues + - "det_FIM": list of determinants of the FIM + - "trace_FIM": list of traces of the FIM - "solve_time": list of solve times - "total_points": total number of points in the factorial design - "success_count": number of successful runs - "failure_count": number of failed runs - "FIM_all": list of all FIMs computed for each point in the factorial - design. - - Raises - ------ - ValueError - If the design_values' keys are not a subset of the model's - `experiment_inputs` keys or if the design_values are not provided. """ # Start timer @@ -1695,73 +1711,91 @@ def compute_FIM_factorial( design_keys = [k for k in model.experiment_inputs.keys()] - # check if abs_change and rel_change are provided and have the same length as - # design_keys - if abs_change: - if len(abs_change) != len(design_keys): - raise ValueError( - "`abs_change` must have the same length of " - f"`{len(design_keys)}` as `design_keys`." - ) - if rel_change: - if len(rel_change) != len(design_keys): + # check if design_values, abs_change and rel_change are provided and have the + # same length as design_keys + # Design values are of higher priority than abs_change and rel_change. + if design_vals is not None: + if len(design_vals) != len(design_keys): raise ValueError( - "`rel_change` must have the same length of " + "`design_values` must have the same length of " f"`{len(design_keys)}` as `design_keys`." ) + design_values = design_vals - # if either abs_change or rel_change is not provided, set it to list of - # zeros - if abs_change or rel_change: - if not abs_change: - abs_change = [0] * len(design_keys) - elif not rel_change: - rel_change = [0] * len(design_keys) - - design_values = [] - # loop over design keys and generate design values - for i, comp in enumerate(design_keys): - lb = comp.lb - ub = comp.ub - # Check if the component has finite lower and upper bounds - if lb is None or ub is None: - raise ValueError(f"{comp.name} does not have a lower or upper bound.") - - if abs_change is None and rel_change is None: - des_val = np.linspace(lb, ub, n_designs) - - elif abs_change or rel_change: - des_val = [] - del_val = comp.lb * rel_change[i] + abs_change[i] - if del_val == 0: + else: + if abs_change: + if len(abs_change) != len(design_keys): raise ValueError( - f"Design variable {comp.name} has no change in value - check " - "abs_change and rel_change values." + "`abs_change` must have the same length of " + f"`{len(design_keys)}` as `design_keys`." ) - val = lb - while val <= ub: - des_val.append(val) - val += del_val - else: - raise ValueError( - "Unexpected error in generating design values. Please check the " - "input parameters." - ) + if rel_change: + if len(rel_change) != len(design_keys): + raise ValueError( + "`rel_change` must have the same length of " + f"`{len(design_keys)}` as `design_keys`." + ) + + # if either abs_change or rel_change is not provided, set it to list of + # zeros + if abs_change or rel_change: + if not abs_change: + abs_change = [0] * len(design_keys) + elif not rel_change: + rel_change = [0] * len(design_keys) + + design_values = [] + # loop over design keys and generate design values + for i, comp in enumerate(design_keys): + lb = comp.lb + ub = comp.ub + # Check if the component has finite lower and upper bounds + if lb is None or ub is None: + raise ValueError( + f"{comp.name} does not have a lower or upper bound." + ) - design_values.append(des_val) + if abs_change is None and rel_change is None: + des_val = np.linspace(lb, ub, n_designs) + + elif abs_change or rel_change: + des_val = [] + del_val = comp.lb * rel_change[i] + abs_change[i] + if del_val == 0: + raise ValueError( + f"Design variable {comp.name} has no change in value - check " + "abs_change and rel_change values." + ) + val = lb + while val <= ub: + des_val.append(val) + val += del_val + + else: + raise ValueError( + "Unexpected error in generating design values. Please check the " + "input parameters." + ) + + design_values.append(des_val) # generate the factorial points based on the initialization scheme - if initialization_scheme == SensitivityInitialization.snake_traversal: - factorial_points = snake_traversal_grid_sampling(*design_values) - elif initialization_scheme == SensitivityInitialization.nested_for_loop: - factorial_points = product(*design_values) - else: + try: + scheme_enum = SensitivityInitialization(initialization_scheme) + except ValueError: self.logger.warning( - "initialization_scheme not recognized. Using `snake_traversal` as " - "default." + f"Initialization scheme '{initialization_scheme}' is not recognized. " + "Using `snake_traversal` as default." ) + scheme_enum = SensitivityInitialization.snake_traversal + + if scheme_enum == SensitivityInitialization.snake_traversal: factorial_points = snake_traversal_grid_sampling(*design_values) + elif scheme_enum == SensitivityInitialization.nested_for_loop: + factorial_points = product(*design_values) + + # TODO: Add more initialization schemes factorial_points_list = list(factorial_points) @@ -1791,13 +1825,14 @@ def compute_FIM_factorial( time_set = [] curr_point = 1 # Initial current point for design_point in factorial_points_list: + # kept the following code to check later whether it is faster. + # In a simple model, this code took 15.9s to compute 125 points in JN + # for the same condition, `update_model_from_suffix` took 16.5s in JN # Fix design variables at fixed experimental design point - for i in range(len(design_point)): - design_keys[i].fix(design_point[i]) + # for i in range(len(design_point)): + # design_keys[i].fix(design_point[i]) - # TODO: check the following lines of code instead of the for loop above - # update_model_from_suffix( - # model, 'experiment_inputs', design_point) + update_model_from_suffix(model, "experiment_inputs", design_point) # Timing and logging objects self.logger.info(f"=======Iteration Number: {curr_point} =======") @@ -1862,31 +1897,40 @@ def compute_FIM_factorial( factorial_results.update( { "total_points": total_points, - "success_counts": success_count, - "failure_counts": failure_count, + "success_count": success_count, + "failure_count": failure_count, "FIM_all": FIM_all.tolist(), # Save all FIMs } ) - if self.tee: - exclude_keys = { - "total_points", - "success_counts", - "failure_counts", - "FIM_all", - } - dict_for_df = { - k: v for k, v in factorial_results.items() if k not in exclude_keys - } - res_df = pd.DataFrame(dict_for_df) - print("\n\n=========Factorial results===========") - print("Total points:", total_points) - print("Success counts:", success_count) - print("Failure counts:", failure_count) - print("\n\n") - print(res_df) - self.factorial_results = factorial_results + # Set pandas DataFrame options + if df_settings[0]: + with pd.option_context( + "display.max_columns", + df_settings[1], + "display.max_rows", + df_settings[2], + "display.width", + df_settings[3], + ): + exclude_keys = { + "total_points", + "success_counts", + "failure_counts", + "FIM_all", + } + dict_for_df = { + k: v for k, v in factorial_results.items() if k not in exclude_keys + } + res_df = pd.DataFrame(dict_for_df) + print("\n\n=========Factorial results===========") + print("Total points:", total_points) + print("Success counts:", success_count) + print("Failure counts:", failure_count) + print("\n") + print(res_df) + # Save the results to a json file based on the file_name provided if file_name is not None: with open(file_name + ".json", "w") as f: From 043e7de97d004f1a887282b6904642c752aba7ab Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 1 Jul 2025 16:00:08 -0400 Subject: [PATCH 13/25] changed some minor comments and docstrings --- pyomo/contrib/doe/doe.py | 32 ++++++++++++++++----------- pyomo/contrib/doe/tests/test_utils.py | 15 +++++-------- 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 386d12b6813..3ccf25ef78c 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1640,7 +1640,9 @@ def compute_FIM_factorial( model : DoE model, optional The model to perform the full factorial exploration on. Default: None design_vals : list, optional - A list of design values to use for the full factorial exploration. + A list of design values to use for the full factorial exploration. If + provided, will use this value to generate the design values and ignore + `abs_change`, `rel_change`, and `n_designs`. Default: None Default: None abs_change : list, optional Absolute change in the design variable values. Default: None. @@ -1759,13 +1761,16 @@ def compute_FIM_factorial( if abs_change is None and rel_change is None: des_val = np.linspace(lb, ub, n_designs) + # if abs_change and/or rel_change is provided, generate design values + # using the formula: + # change_in_value = lower_bound * rel_change + abs_change elif abs_change or rel_change: des_val = [] del_val = comp.lb * rel_change[i] + abs_change[i] if del_val == 0: raise ValueError( - f"Design variable {comp.name} has no change in value - check " - "abs_change and rel_change values." + f"Design variable {comp.name} has no change in value - " + "check abs_change and rel_change values." ) val = lb while val <= ub: @@ -1774,8 +1779,8 @@ def compute_FIM_factorial( else: raise ValueError( - "Unexpected error in generating design values. Please check the " - "input parameters." + "Unexpected error in generating design values. Please check the" + " input parameters." ) design_values.append(des_val) @@ -1786,7 +1791,7 @@ def compute_FIM_factorial( except ValueError: self.logger.warning( f"Initialization scheme '{initialization_scheme}' is not recognized. " - "Using `snake_traversal` as default." + "Using `snake_traversal` as the default initialization scheme." ) scheme_enum = SensitivityInitialization.snake_traversal @@ -1825,9 +1830,10 @@ def compute_FIM_factorial( time_set = [] curr_point = 1 # Initial current point for design_point in factorial_points_list: - # kept the following code to check later whether it is faster. - # In a simple model, this code took 15.9s to compute 125 points in JN - # for the same condition, `update_model_from_suffix` took 16.5s in JN + # kept (commented out) the following code to check later whether it will + # be considerably faster for complex models. + # In a simple model, this code took 15.9s to compute 125 points in JN. + # For the same condition, `update_model_from_suffix` took 16.5s in JN # Fix design variables at fixed experimental design point # for i in range(len(design_point)): # design_keys[i].fix(design_point[i]) @@ -1916,8 +1922,8 @@ def compute_FIM_factorial( ): exclude_keys = { "total_points", - "success_counts", - "failure_counts", + "success_count", + "failure_count", "FIM_all", } dict_for_df = { @@ -1926,8 +1932,8 @@ def compute_FIM_factorial( res_df = pd.DataFrame(dict_for_df) print("\n\n=========Factorial results===========") print("Total points:", total_points) - print("Success counts:", success_count) - print("Failure counts:", failure_count) + print("Success count:", success_count) + print("Failure count:", failure_count) print("\n") print(res_df) diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py index eb3956176c8..5503f8c0c44 100644 --- a/pyomo/contrib/doe/tests/test_utils.py +++ b/pyomo/contrib/doe/tests/test_utils.py @@ -148,16 +148,14 @@ def test_snake_traversal_grid_sampling_errors(self): with self.assertRaises(ValueError) as cm: list(snake_traversal_grid_sampling(list_2d_bad)) self.assertEqual( - str(cm.exception), - "Argument at position 0 is not 1D. Got shape (2, 3).", + str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3)." ) list_2d_wrong_shape_bad = [[1, 2, 3], [4, 5, 6, 7]] with self.assertRaises(ValueError) as cm: list(snake_traversal_grid_sampling(list_2d_wrong_shape_bad)) self.assertEqual( - str(cm.exception), - "Argument at position 0 is not 1D array-like.", + str(cm.exception), "Argument at position 0 is not 1D array-like." ) # Test the error handling with tuples @@ -165,16 +163,14 @@ def test_snake_traversal_grid_sampling_errors(self): with self.assertRaises(ValueError) as cm: list(snake_traversal_grid_sampling(tuple_2d_bad)) self.assertEqual( - str(cm.exception), - "Argument at position 0 is not 1D. Got shape (2, 3).", + str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3)." ) tuple_2d_wrong_shape_bad = ((1, 2, 3), (4, 5, 6, 7)) with self.assertRaises(ValueError) as cm: list(snake_traversal_grid_sampling(tuple_2d_wrong_shape_bad)) self.assertEqual( - str(cm.exception), - "Argument at position 0 is not 1D array-like.", + str(cm.exception), "Argument at position 0 is not 1D array-like." ) # Test the error handling with numpy arrays @@ -182,8 +178,7 @@ def test_snake_traversal_grid_sampling_errors(self): with self.assertRaises(ValueError) as cm: list(snake_traversal_grid_sampling(array_2d_bad)) self.assertEqual( - str(cm.exception), - "Argument at position 0 is not 1D. Got shape (2, 3).", + str(cm.exception), "Argument at position 0 is not 1D. Got shape (2, 3)." ) def test_snake_traversal_grid_sampling_values(self): From 33e1bd0717f9d7800c5c71994366c5c05a5b447e Mon Sep 17 00:00:00 2001 From: Stephen Cini <114932899+sscini@users.noreply.github.com> Date: Tue, 1 Jul 2025 16:31:27 -0400 Subject: [PATCH 14/25] updated with sscini's `utils.py`. --- pyomo/contrib/doe/utils.py | 46 +++++++++++++++----------------------- 1 file changed, 18 insertions(+), 28 deletions(-) diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index bb1fcf881f7..d5e428dee14 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -126,47 +126,37 @@ def rescale_FIM(FIM, param_vals): # Adding utility to update parameter values in a model based on the suffix -def update_model_from_suffix(model, suffix_name, values): +def update_model_from_suffix(suffix_obj: pyo.Suffix, values): """ - Iterate over the components (variables or parameters) referenced by the - given suffix in the model, and assign each a new value from the provided iterable. + Overwrite each variable/parameter referenced by ``suffix_obj`` with the + corresponding value in ``values``. Parameters ---------- - model : pyomo.environ.ConcreteModel - The Pyomo model containing the suffix and components to update. - suffix_name : str - The name of the Suffix attribute on the model whose items will be updated. - Must be one of: 'experiment_outputs', 'experiment_inputs', 'unknown_parameters', or 'measurement_error'. + suffix_obj : pyomo.core.base.suffix.Suffix + The suffix whose *keys* are the components you want to update. + Call like ``update_from_suffix(model.unknown_parameters, vals)``. values : iterable of numbers - The new values to assign to each component referenced by the suffix. The length of this - iterable must match the number of items in the suffix. + New numerical values for the components referenced by the suffix. + Must be the same length as ``suffix_obj``. """ - # Allowed suffix names - allowed = { - "experiment_outputs", - "experiment_inputs", - "unknown_parameters", - "measurement_error", - } - # Validate input is an allowed suffix name - if suffix_name not in allowed: - raise ValueError(f"suffix_name must be one of {sorted(allowed)}") - # Check if the model has the specified suffix - suffix_obj = getattr(model, suffix_name, None) - if suffix_obj is None: - raise AttributeError(f"Model has no attribute '{suffix_name}'") - # Check if the suffix is a Suffix object + # Check that the length of values matches the suffix length items = list(suffix_obj.items()) if len(items) != len(values): raise ValueError("values length does not match suffix length") - # Set the new values for the suffix items + + # Iterate through the items in the suffix and update their values + # Note: items are tuples of (component, suffix_value) for (comp, _), new_val in zip(items, values): - # Update the variable/parameter itself if it is VarData or ParamData + + # update the component value + # Check if the component is a VarData or ParamData if isinstance(comp, (VarData, ParamData)): comp.set_value(new_val) else: - raise TypeError(f"Unsupported component type: {type(comp)}") + raise TypeError( + f"Unsupported component type {type(comp)}; expected VarData or ParamData." + ) def check_FIM(FIM): From 553ab82d0658a27d570ceba422ad1047f07d072c Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Wed, 9 Jul 2025 07:32:41 -0400 Subject: [PATCH 15/25] added new func and changed test files added a new correlation matrix computation in doe/utils.py, added TODOs. --- pyomo/contrib/doe/doe.py | 7 +- pyomo/contrib/doe/tests/test_doe_errors.py | 4 +- pyomo/contrib/doe/tests/test_utils.py | 60 ++++++++++---- pyomo/contrib/doe/utils.py | 95 +++++++++++++++++----- 4 files changed, 128 insertions(+), 38 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 88d09720b9d..50c00508a6b 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -46,7 +46,7 @@ import pyomo.environ as pyo from pyomo.contrib.doe.utils import ( - check_FIM, + check_matrix, compute_FIM_metrics, _SMALL_TOLERANCE_DEFINITENESS, snake_traversal_grid_sampling, @@ -629,6 +629,9 @@ def _sequential_FIM(self, model=None): if self.scale_nominal_param_value: scale_factor *= v + # TODO: scale by response values (i.e., measurement values) + # if self.scale_response_values: + # scale_factor /= measurement_vals_np[:, col_1].mean() # Calculate the column of the sensitivity matrix self.seq_jac[:, i] = ( measurement_vals_np[:, col_1] - measurement_vals_np[:, col_2] @@ -1384,7 +1387,7 @@ def check_model_FIM(self, model=None, FIM=None): ) # Check FIM is positive definite and symmetric - check_FIM(FIM) + check_matrix(FIM) self.logger.info( "FIM provided matches expected dimensions from model and is approximately positive (semi) definite." diff --git a/pyomo/contrib/doe/tests/test_doe_errors.py b/pyomo/contrib/doe/tests/test_doe_errors.py index ce77fb4f553..578c3c253c5 100644 --- a/pyomo/contrib/doe/tests/test_doe_errors.py +++ b/pyomo/contrib/doe/tests/test_doe_errors.py @@ -183,7 +183,7 @@ def test_reactor_check_bad_prior_negative_eigenvalue(self): with self.assertRaisesRegex( ValueError, - r"FIM provided is not positive definite. It has one or more negative eigenvalue\(s\) less than -{:.1e}".format( + r"Matrix provided is not positive definite. It has one or more negative eigenvalue\(s\) less than -{:.1e}".format( _SMALL_TOLERANCE_DEFINITENESS ), ): @@ -208,7 +208,7 @@ def test_reactor_check_bad_prior_not_symmetric(self): with self.assertRaisesRegex( ValueError, - "FIM provided is not symmetric using absolute tolerance {}".format( + "Matrix provided is not symmetric using absolute tolerance {}".format( _SMALL_TOLERANCE_SYMMETRY ), ): diff --git a/pyomo/contrib/doe/tests/test_utils.py b/pyomo/contrib/doe/tests/test_utils.py index cba102bce4b..20c63b1060d 100644 --- a/pyomo/contrib/doe/tests/test_utils.py +++ b/pyomo/contrib/doe/tests/test_utils.py @@ -8,14 +8,20 @@ # rights in this software. # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -from pyomo.common.dependencies import numpy as np, numpy_available +from pyomo.common.dependencies import ( + numpy as np, + numpy_available, + pandas as pd, + pandas_available, +) import pyomo.common.unittest as unittest from pyomo.contrib.doe.utils import ( - check_FIM, + check_matrix, compute_FIM_metrics, get_FIM_metrics, snake_traversal_grid_sampling, + compute_correlation_matrix as compcorr, _SMALL_TOLERANCE_DEFINITENESS, _SMALL_TOLERANCE_SYMMETRY, _SMALL_TOLERANCE_IMG, @@ -23,45 +29,49 @@ @unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not pandas_available, "Pandas is not available") class TestUtilsFIM(unittest.TestCase): - """Test the check_FIM() from utils.py.""" + """Test the check_matrix() from utils.py.""" - def test_check_FIM_valid(self): + # TODO: add tests when `check_pos_def = False` is used in check_matrix() + def test_check_matrix_valid(self): """Test case where the FIM is valid (square, positive definite, symmetric).""" FIM = np.array([[4, 1], [1, 3]]) try: - check_FIM(FIM) + check_matrix(FIM) except ValueError as e: self.fail(f"Unexpected error: {e}") - def test_check_FIM_non_square(self): + def test_check_matrix_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) + with self.assertRaisesRegex( + ValueError, "argument mat must be a 2D square matrix" + ): + check_matrix(FIM) - def test_check_FIM_non_positive_definite(self): + def test_check_matrix_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 " + "Matrix provided is not positive definite. It has one or more negative " + r"eigenvalue\(s\) less than -{:.1e}".format( _SMALL_TOLERANCE_DEFINITENESS ), ): - check_FIM(FIM) + check_matrix(FIM) - def test_check_FIM_non_symmetric(self): + def test_check_matrix_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( + "Matrix provided is not symmetric using absolute tolerance {}".format( _SMALL_TOLERANCE_SYMMETRY ), ): - check_FIM(FIM) + check_matrix(FIM) """Test the compute_FIM_metrics() from utils.py.""" @@ -259,6 +269,28 @@ def test_snake_traversal_grid_sampling_values(self): result_mixed = list(snake_traversal_grid_sampling(list1, tuple2, array3)) self.assertEqual(result_mixed, expected_list3) + # TODO: Add more tests as needed + def test_compute_correlation_matrix(self): + # Create a sample covariance matrix + covariance_matrix = np.array([[4, 2], [2, 3]]) + var_name = ["X1", "X2"] + + # Compute the correlation matrix + correlation_matrix = compcorr(covariance_matrix, var_name) + + # Expected correlation matrix + expected_correlation_matrix = pd.DataFrame( + [[1.0, 0.577], [0.577, 1.0]], index=var_name, columns=var_name + ) + + # Check if the computed correlation matrix matches the expected one + pd.testing.assert_frame_equal( + correlation_matrix, + expected_correlation_matrix, + check_exact=False, + atol=1e-6, + ) + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 222abf93d9a..309ad4320b3 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -27,7 +27,7 @@ import pyomo.environ as pyo -from pyomo.common.dependencies import numpy as np, numpy_available +from pyomo.common.dependencies import numpy as np, numpy_available, pandas as pd from pyomo.core.base.param import ParamData from pyomo.core.base.var import VarData @@ -159,41 +159,46 @@ def update_model_from_suffix(suffix_obj: pyo.Suffix, values): ) -def check_FIM(FIM): +def check_matrix(mat, check_pos_def=True): """ - Checks that the FIM is square, positive definite, and symmetric. + Checks that the Matrix is square, positive definite, and symmetric. Parameters ---------- - FIM: 2D numpy array representing the FIM + mat: 2D numpy array representing the matrix 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") + # Check that the matrix is a square matrix + if mat.ndim != 2 or mat.shape[0] != mat.shape[1]: + raise ValueError("argument mat must be a 2D square matrix") - # Compute the eigenvalues of the FIM - evals = np.linalg.eigvals(FIM) + if check_pos_def: + # Compute the eigenvalues of the matrix + evals = np.linalg.eigvals(mat) - # 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 matrix is positive definite + if np.min(evals) < -_SMALL_TOLERANCE_DEFINITENESS: + raise ValueError( + "Matrix 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): + # Check if the matrix is symmetric + if not np.allclose(mat, mat.T, atol=_SMALL_TOLERANCE_SYMMETRY): raise ValueError( - "FIM provided is not symmetric using absolute tolerance {}".format( + "Matrix provided is not symmetric using absolute tolerance {}".format( _SMALL_TOLERANCE_SYMMETRY ) ) +# TODO: probably can merge compute_FIM_metrics() and get_FIM_metrics() in a single + +# TODO: function with an argument (e.g., return_dict = True) to compute the metrics # Functions to compute FIM metrics def compute_FIM_metrics(FIM): """ @@ -225,7 +230,7 @@ def compute_FIM_metrics(FIM): """ # Check whether the FIM is square, positive definite, and symmetric - check_FIM(FIM) + check_matrix(FIM) # Compute FIM metrics det_FIM = np.linalg.det(FIM) @@ -315,7 +320,6 @@ def snake_traversal_grid_sampling(*array_like_args): Yields ------ A tuple representing points in the snake pattern. - Naming source of the pattern: https://dev.to/thesanjeevsharma/dsa-101-matrix-30lf """ # throw an error if the array_like_args are not array-like or all 1D for i, arg in enumerate(array_like_args): @@ -352,3 +356,54 @@ def _generate_recursive(depth, index_sum): # Start the recursion at the first list (depth 0) with an initial sum of 0. yield from _generate_recursive(depth=0, index_sum=0) + + +def compute_correlation_matrix( + covariance_matrix, var_name: list = None, significant_digits=3 +): + """ + Computes the correlation matrix from a covariance matrix. + + Parameters + ---------- + covariance_matrix : numpy.ndarray + 2D array representing the covariance matrix. + var_name : list, optional + List of variable names corresponding to the rows/columns of the covariance matrix. + significant_digits : int, optional + Number of significant digits to round the correlation matrix to. Default: 3. + + Returns + ------- + pandas.DataFrame/numpy.ndarray + If `var_name` is provided, returns a pandas DataFrame with the correlation matrix + and the specified variable names as both index and columns. If `var_name` is not + provided, returns a numpy array representing the correlation matrix. + """ + # Check if covariance matrix is symmetric and square + check_matrix(covariance_matrix, check_pos_def=False) + + if var_name: + assert len(var_name) == covariance_matrix.shape[0], ( + "Length of var_name must match the number of rows/columns in the " + "covariance matrix." + ) + + if not np.all(np.isfinite(covariance_matrix)): + raise ValueError("Covariance matrix contains non-finite values.") + + std_dev = np.sqrt(np.diag(covariance_matrix)) + + std_dev_matrix = np.outer(std_dev, std_dev) + + correlation_matrix = covariance_matrix / std_dev_matrix + + # Set the index to the variable names if provided, + if var_name is not None: + corr_df = pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) + else: + corr_df = correlation_matrix + + return ( + corr_df.round(significant_digits) if significant_digits else correlation_matrix + ) From a063c3384cfe44da9ba642bdf01e20314f7ae548 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Wed, 9 Jul 2025 07:56:37 -0400 Subject: [PATCH 16/25] changed if statement to ternary operator and minor changes in return statement --- pyomo/contrib/doe/tests/test_solve.ipynb | 258 +++++++++++++++++++++++ pyomo/contrib/doe/utils.py | 18 +- 2 files changed, 268 insertions(+), 8 deletions(-) create mode 100644 pyomo/contrib/doe/tests/test_solve.ipynb diff --git a/pyomo/contrib/doe/tests/test_solve.ipynb b/pyomo/contrib/doe/tests/test_solve.ipynb new file mode 100644 index 00000000000..55a5977a5fc --- /dev/null +++ b/pyomo/contrib/doe/tests/test_solve.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "e1ef738a", + "metadata": {}, + "outputs": [], + "source": [ + "from pyomo.contrib.doe.examples.reactor_example import run_reactor_doe" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "1e4828b4", + "metadata": {}, + "outputs": [], + "source": [ + "ff = run_reactor_doe(\n", + " n_points_for_design=2,\n", + " compute_FIM_full_factorial=False,\n", + " plot_factorial_results=False,\n", + " save_plots=False,\n", + " run_optimal_doe=False,\n", + ")\n", + "ff_results = ff.compute_FIM_full_factorial(\n", + " design_ranges={\n", + " \"CA[0]\": [1, 1.5, 2],\n", + " \"T[0]\": [350, 400, 2],\n", + " }\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "58c50f5a", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "473279cc", + "metadata": {}, + "outputs": [], + "source": [ + "log10_D_opt_expected = [float(x) for x in ff_results[\"log10 D-opt\"]]\n", + "log10_A_opt_expected = [float(x) for x in ff_results[\"log10 A-opt\"]]\n", + "log10_E_opt_expected = [float(x) for x in ff_results[\"log10 E-opt\"]]\n", + "log10_ME_opt_expected = [float(x) for x in ff_results[\"log10 ME-opt\"]]\n", + "eigval_min_expected = [float(x) for x in ff_results[\"eigval_min\"]]\n", + "eigval_max_expected = [float(x) for x in ff_results[\"eigval_max\"]]\n", + "det_FIM_expected = [float(x) for x in ff_results[\"det_FIM\"]]\n", + "trace_FIM_expected = [float(x) for x in ff_results[\"trace_FIM\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "3df9f882", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[3.7734377852467524, 5.137792359070963, 5.182167857710023, 6.546522431509408]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "log10_D_opt_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "d73b39b2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[3.5935726800929695, 3.6133186151486933, 3.945755198204365, 3.9655011332598367]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "log10_A_opt_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b84bba5f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[-1.7201873126109162,\n", + " -0.691340497355524,\n", + " -1.3680047944877138,\n", + " -0.3391579792516522]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "log10_E_opt_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "75a3aff9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[5.221185311075697, 4.244741560076784, 5.221185311062606, 4.244741560083524]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "log10_ME_opt_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "aa441b5c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[0.019046390638130666,\n", + " 0.20354456134677426,\n", + " 0.04285437893696232,\n", + " 0.45797526302234304]" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigval_min_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "ec0bdc18", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[3169.552855492114, 3576.0292523637977, 7131.493924857995, 8046.0658178139165]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "eigval_max_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "148b371d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[5935.233170586055, 137338.51875774842, 152113.5345070818, 3519836.021699428]" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "det_FIM_expected" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "461b488d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[3922.5878617108597, 4105.051549241871, 8825.822688850109, 9236.36598578955]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "trace_FIM_expected" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pyodev", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 309ad4320b3..20c93fc73d5 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -161,11 +161,14 @@ def update_model_from_suffix(suffix_obj: pyo.Suffix, values): def check_matrix(mat, check_pos_def=True): """ - Checks that the Matrix is square, positive definite, and symmetric. + Checks that the matrix is square, positive definite, and symmetric. Parameters ---------- mat: 2D numpy array representing the matrix + check_pos_def: bool, optional + If True, checks if the matrix is positive definite. + Default: True. Returns ------- @@ -399,11 +402,10 @@ def compute_correlation_matrix( correlation_matrix = covariance_matrix / std_dev_matrix # Set the index to the variable names if provided, - if var_name is not None: - corr_df = pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) - else: - corr_df = correlation_matrix - - return ( - corr_df.round(significant_digits) if significant_digits else correlation_matrix + corr_df = ( + pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) + if var_name + else correlation_matrix ) + + return corr_df.round(significant_digits) if significant_digits else corr_df From 9af160d87615e3ff50e4f9bcc11150af75ce7a5c Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Wed, 9 Jul 2025 07:59:41 -0400 Subject: [PATCH 17/25] Revert "changed if statement to ternary operator and minor changes in return statement" This reverts commit a063c3384cfe44da9ba642bdf01e20314f7ae548. --- pyomo/contrib/doe/utils.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 20c93fc73d5..309ad4320b3 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -161,14 +161,11 @@ def update_model_from_suffix(suffix_obj: pyo.Suffix, values): def check_matrix(mat, check_pos_def=True): """ - Checks that the matrix is square, positive definite, and symmetric. + Checks that the Matrix is square, positive definite, and symmetric. Parameters ---------- mat: 2D numpy array representing the matrix - check_pos_def: bool, optional - If True, checks if the matrix is positive definite. - Default: True. Returns ------- @@ -402,10 +399,11 @@ def compute_correlation_matrix( correlation_matrix = covariance_matrix / std_dev_matrix # Set the index to the variable names if provided, - corr_df = ( - pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) - if var_name - else correlation_matrix - ) + if var_name is not None: + corr_df = pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) + else: + corr_df = correlation_matrix - return corr_df.round(significant_digits) if significant_digits else corr_df + return ( + corr_df.round(significant_digits) if significant_digits else correlation_matrix + ) From 1f849a05a4410f2b63dcd199abd00773d5e8f0b7 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Wed, 9 Jul 2025 08:01:14 -0400 Subject: [PATCH 18/25] Reapply "changed if statement to ternary operator and minor changes in return statement" This reverts commit 9af160d87615e3ff50e4f9bcc11150af75ce7a5c. --- pyomo/contrib/doe/utils.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 309ad4320b3..20c93fc73d5 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -161,11 +161,14 @@ def update_model_from_suffix(suffix_obj: pyo.Suffix, values): def check_matrix(mat, check_pos_def=True): """ - Checks that the Matrix is square, positive definite, and symmetric. + Checks that the matrix is square, positive definite, and symmetric. Parameters ---------- mat: 2D numpy array representing the matrix + check_pos_def: bool, optional + If True, checks if the matrix is positive definite. + Default: True. Returns ------- @@ -399,11 +402,10 @@ def compute_correlation_matrix( correlation_matrix = covariance_matrix / std_dev_matrix # Set the index to the variable names if provided, - if var_name is not None: - corr_df = pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) - else: - corr_df = correlation_matrix - - return ( - corr_df.round(significant_digits) if significant_digits else correlation_matrix + corr_df = ( + pd.DataFrame(correlation_matrix, index=var_name, columns=var_name) + if var_name + else correlation_matrix ) + + return corr_df.round(significant_digits) if significant_digits else corr_df From da91f2f4318acc7447297431f8b4322062b20b3b Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Wed, 9 Jul 2025 08:01:59 -0400 Subject: [PATCH 19/25] deleted test_solve --- pyomo/contrib/doe/tests/test_solve.ipynb | 258 ----------------------- 1 file changed, 258 deletions(-) delete mode 100644 pyomo/contrib/doe/tests/test_solve.ipynb diff --git a/pyomo/contrib/doe/tests/test_solve.ipynb b/pyomo/contrib/doe/tests/test_solve.ipynb deleted file mode 100644 index 55a5977a5fc..00000000000 --- a/pyomo/contrib/doe/tests/test_solve.ipynb +++ /dev/null @@ -1,258 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "e1ef738a", - "metadata": {}, - "outputs": [], - "source": [ - "from pyomo.contrib.doe.examples.reactor_example import run_reactor_doe" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "1e4828b4", - "metadata": {}, - "outputs": [], - "source": [ - "ff = run_reactor_doe(\n", - " n_points_for_design=2,\n", - " compute_FIM_full_factorial=False,\n", - " plot_factorial_results=False,\n", - " save_plots=False,\n", - " run_optimal_doe=False,\n", - ")\n", - "ff_results = ff.compute_FIM_full_factorial(\n", - " design_ranges={\n", - " \"CA[0]\": [1, 1.5, 2],\n", - " \"T[0]\": [350, 400, 2],\n", - " }\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "58c50f5a", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "473279cc", - "metadata": {}, - "outputs": [], - "source": [ - "log10_D_opt_expected = [float(x) for x in ff_results[\"log10 D-opt\"]]\n", - "log10_A_opt_expected = [float(x) for x in ff_results[\"log10 A-opt\"]]\n", - "log10_E_opt_expected = [float(x) for x in ff_results[\"log10 E-opt\"]]\n", - "log10_ME_opt_expected = [float(x) for x in ff_results[\"log10 ME-opt\"]]\n", - "eigval_min_expected = [float(x) for x in ff_results[\"eigval_min\"]]\n", - "eigval_max_expected = [float(x) for x in ff_results[\"eigval_max\"]]\n", - "det_FIM_expected = [float(x) for x in ff_results[\"det_FIM\"]]\n", - "trace_FIM_expected = [float(x) for x in ff_results[\"trace_FIM\"]]" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "3df9f882", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[3.7734377852467524, 5.137792359070963, 5.182167857710023, 6.546522431509408]" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "log10_D_opt_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "d73b39b2", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[3.5935726800929695, 3.6133186151486933, 3.945755198204365, 3.9655011332598367]" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "log10_A_opt_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "b84bba5f", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[-1.7201873126109162,\n", - " -0.691340497355524,\n", - " -1.3680047944877138,\n", - " -0.3391579792516522]" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "log10_E_opt_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "75a3aff9", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[5.221185311075697, 4.244741560076784, 5.221185311062606, 4.244741560083524]" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "log10_ME_opt_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "aa441b5c", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[0.019046390638130666,\n", - " 0.20354456134677426,\n", - " 0.04285437893696232,\n", - " 0.45797526302234304]" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "eigval_min_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "id": "ec0bdc18", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[3169.552855492114, 3576.0292523637977, 7131.493924857995, 8046.0658178139165]" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "eigval_max_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "148b371d", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[5935.233170586055, 137338.51875774842, 152113.5345070818, 3519836.021699428]" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "det_FIM_expected" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "461b488d", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[3922.5878617108597, 4105.051549241871, 8825.822688850109, 9236.36598578955]" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "trace_FIM_expected" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "pyodev", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} From fd8d3eaea87dd2fa7bb8ccaf1f71d51c23bdb4be Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Mon, 28 Jul 2025 14:24:53 -0400 Subject: [PATCH 20/25] added a strategy to compute the fractional factorial as well. --- pyomo/contrib/doe/doe.py | 63 +++++++++++++++++++++++++--------------- 1 file changed, 39 insertions(+), 24 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 50c00508a6b..7e6f9f5fae8 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1625,7 +1625,7 @@ def compute_FIM_full_factorial( def compute_FIM_factorial( self, model=None, - design_vals: list = None, + design_vals: dict = None, abs_change: list = None, rel_change: list = None, n_designs: int = 5, @@ -1643,11 +1643,19 @@ def compute_FIM_factorial( ---------- model : DoE model, optional The model to perform the full factorial exploration on. Default: None - design_vals : list, optional - A list of design values to use for the full factorial exploration. If - provided, will use this value to generate the design values and ignore - `abs_change`, `rel_change`, and `n_designs`. Default: None - Default: None + design_values : dict, + dict of lists or other array-like objects, of the form + {"var_name": }. Default: None. + The `design_values` should have the key(s) passed as strings that is a + subset of the `experiment_inputs`. If one or more design variables are not + to be changed, then they should not be passed in the `design_values` + dictionary, but if they are passed in the dictionary, then they must be a + array-like object of floats. For example, if our experiment has 4 design variables + (i.e., `experiment_inputs`): model.x1, model.x2, model.x3, and model.x4, + their values may be passed as, design_values= {"x1": [1, 2, 3], "x3": [7], + "x4": [-10, 20]}. In this case, x2 will not be changed and will be fixed at + the value in the model. + If design_values is provided, then `abs_change` and `rel_change` are ignored. abs_change : list, optional Absolute change in the design variable values. Default: None. If provided, will use this value to generate the design values. @@ -1715,20 +1723,30 @@ def compute_FIM_factorial( ).clone() model = self.factorial_model - design_keys = [k for k in model.experiment_inputs.keys()] - - # check if design_values, abs_change and rel_change are provided and have the - # same length as design_keys - # Design values are of higher priority than abs_change and rel_change. if design_vals is not None: - if len(design_vals) != len(design_keys): + # Check whether the design_ranges keys are in the experiment_inputs + design_keys = set(design_vals.keys()) + map_keys = set([k.name for k in model.experiment_inputs.keys()]) + if not design_keys.issubset(map_keys): + incorrect_given_keys = design_keys - map_keys + suggested_keys = map_keys - design_keys raise ValueError( - "`design_values` must have the same length of " - f"`{len(design_keys)}` as `design_keys`." + f"design_values keys: {incorrect_given_keys} are incorrect." + f"The keys should be from the following keys: {suggested_keys}." ) - design_values = design_vals + + # Get the design map keys that match the design_values keys + design_map_keys = [ + k + for k in model.experiment_inputs.keys() + if k.name in design_vals.keys() + ] + # This ensures that the order of the design_values keys matches the order of the + # design_map_keys so that design_point can be constructed correctly in the loop. + design_values = [design_vals[k.name] for k in design_map_keys] else: + design_keys = [k for k in model.experiment_inputs.keys()] if abs_change: if len(abs_change) != len(design_keys): raise ValueError( @@ -1834,15 +1852,12 @@ def compute_FIM_factorial( time_set = [] curr_point = 1 # Initial current point for design_point in factorial_points_list: - # kept (commented out) the following code to check later whether it will - # be considerably faster for complex models. - # In a simple model, this code took 15.9s to compute 125 points in JN. - # For the same condition, `update_model_from_suffix` took 16.5s in JN - # Fix design variables at fixed experimental design point - # for i in range(len(design_point)): - # design_keys[i].fix(design_point[i]) - - update_model_from_suffix(model, "experiment_inputs", design_point) + if design_vals: + for i in range(len(design_point)): + # Set the design variable value from the design_vals dictionary + design_map_keys[i].set_value(design_point[i]) + else: + update_model_from_suffix(model.experiment_inputs, design_point) # Timing and logging objects self.logger.info(f"=======Iteration Number: {curr_point} =======") From 0ee68567f505abc57cfb0b72c6d220362330006a Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 29 Jul 2025 12:35:17 -0400 Subject: [PATCH 21/25] Created a temporary suffix `design_suff` to contain the keys of the user provided design variables --- pyomo/contrib/doe/doe.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 7e6f9f5fae8..f5ed077a990 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1643,7 +1643,7 @@ def compute_FIM_factorial( ---------- model : DoE model, optional The model to perform the full factorial exploration on. Default: None - design_values : dict, + design_vals : dict, dict of lists or other array-like objects, of the form {"var_name": }. Default: None. The `design_values` should have the key(s) passed as strings that is a @@ -1745,6 +1745,10 @@ def compute_FIM_factorial( # design_map_keys so that design_point can be constructed correctly in the loop. design_values = [design_vals[k.name] for k in design_map_keys] + # Create a temporary suffix to pass in `update_model_from_suffix` + design_suff = pyo.Suffix(direction=pyo.Suffix.LOCAL) + design_suff.update((k, None) for k in design_map_keys) + else: design_keys = [k for k in model.experiment_inputs.keys()] if abs_change: @@ -1853,9 +1857,7 @@ def compute_FIM_factorial( curr_point = 1 # Initial current point for design_point in factorial_points_list: if design_vals: - for i in range(len(design_point)): - # Set the design variable value from the design_vals dictionary - design_map_keys[i].set_value(design_point[i]) + update_model_from_suffix(design_suff, design_point) else: update_model_from_suffix(model.experiment_inputs, design_point) From c6e81ea9770182c77f214b924abeba5fa8413def Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 16 Sep 2025 21:34:21 -0400 Subject: [PATCH 22/25] changed the import location for `update_model_from_suffix` and resolved conflict of doe/utils.py --- pyomo/contrib/doe/doe.py | 3 +- pyomo/contrib/doe/utils.py | 74 ++++++++++++++++++++++++++++++++++---- 2 files changed, 69 insertions(+), 8 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 81af87aae84..0fb379b88d8 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -57,9 +57,10 @@ compute_FIM_metrics, _SMALL_TOLERANCE_DEFINITENESS, snake_traversal_grid_sampling, - update_model_from_suffix, ) +from pyomo.contrib.parmest.utils.model_utils import update_model_from_suffix + from pyomo.opt import SolverStatus diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 67519284624..20c93fc73d5 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -92,11 +92,74 @@ def rescale_FIM(FIM, param_vals): return scaled_FIM -<<<<<<< HEAD +# TODO: Add swapping parameters for variables helper function +# def get_parameters_from_suffix(suffix, fix_vars=False): +# """ +# Finds the Params within the suffix provided. It will also check to see +# if there are Vars in the suffix provided. ``fix_vars`` will indicate +# if we should fix all the Vars in the set or not. +# +# Parameters +# ---------- +# suffix: pyomo Suffix object, contains the components to be checked +# as keys +# fix_vars: boolean, whether or not to fix the Vars, default = False +# +# Returns +# ------- +# param_list: list of Param +# """ +# param_list = [] +# +# # FIX THE MODEL TREE ISSUE WHERE I GET base_model. INSTEAD OF +# # Check keys if they are Param or Var. Fix the vars if ``fix_vars`` is True +# for k, v in suffix.items(): +# if isinstance(k, ParamData): +# param_list.append(k.name) +# elif isinstance(k, VarData): +# if fix_vars: +# k.fix() +# else: +# pass # ToDo: Write error for suffix keys that aren't ParamData or VarData +# +# return param_list + + +# Adding utility to update parameter values in a model based on the suffix +def update_model_from_suffix(suffix_obj: pyo.Suffix, values): + """ + Overwrite each variable/parameter referenced by ``suffix_obj`` with the + corresponding value in ``values``. + + Parameters + ---------- + suffix_obj : pyomo.core.base.suffix.Suffix + The suffix whose *keys* are the components you want to update. + Call like ``update_from_suffix(model.unknown_parameters, vals)``. + values : iterable of numbers + New numerical values for the components referenced by the suffix. + Must be the same length as ``suffix_obj``. + """ + # Check that the length of values matches the suffix length + items = list(suffix_obj.items()) + if len(items) != len(values): + raise ValueError("values length does not match suffix length") + + # Iterate through the items in the suffix and update their values + # Note: items are tuples of (component, suffix_value) + for (comp, _), new_val in zip(items, values): + + # update the component value + # Check if the component is a VarData or ParamData + if isinstance(comp, (VarData, ParamData)): + comp.set_value(new_val) + else: + raise TypeError( + f"Unsupported component type {type(comp)}; expected VarData or ParamData." + ) + + def check_matrix(mat, check_pos_def=True): -======= -def check_FIM(FIM): ->>>>>>> upstream/main """ Checks that the matrix is square, positive definite, and symmetric. @@ -245,7 +308,6 @@ def get_FIM_metrics(FIM): "log10(E-Optimality)": E_opt, "log10(Modified E-Optimality)": ME_opt, } -<<<<<<< HEAD def snake_traversal_grid_sampling(*array_like_args): @@ -347,5 +409,3 @@ def compute_correlation_matrix( ) return corr_df.round(significant_digits) if significant_digits else corr_df -======= ->>>>>>> upstream/main From aa3961f757d4ed9212d1c28cd941ce378d2b1421 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 16 Sep 2025 21:41:27 -0400 Subject: [PATCH 23/25] cleaned doe/utils.py --- pyomo/contrib/doe/utils.py | 67 -------------------------------------- 1 file changed, 67 deletions(-) diff --git a/pyomo/contrib/doe/utils.py b/pyomo/contrib/doe/utils.py index 20c93fc73d5..cc396a00288 100644 --- a/pyomo/contrib/doe/utils.py +++ b/pyomo/contrib/doe/utils.py @@ -92,73 +92,6 @@ def rescale_FIM(FIM, param_vals): return scaled_FIM -# TODO: Add swapping parameters for variables helper function -# def get_parameters_from_suffix(suffix, fix_vars=False): -# """ -# Finds the Params within the suffix provided. It will also check to see -# if there are Vars in the suffix provided. ``fix_vars`` will indicate -# if we should fix all the Vars in the set or not. -# -# Parameters -# ---------- -# suffix: pyomo Suffix object, contains the components to be checked -# as keys -# fix_vars: boolean, whether or not to fix the Vars, default = False -# -# Returns -# ------- -# param_list: list of Param -# """ -# param_list = [] -# -# # FIX THE MODEL TREE ISSUE WHERE I GET base_model. INSTEAD OF -# # Check keys if they are Param or Var. Fix the vars if ``fix_vars`` is True -# for k, v in suffix.items(): -# if isinstance(k, ParamData): -# param_list.append(k.name) -# elif isinstance(k, VarData): -# if fix_vars: -# k.fix() -# else: -# pass # ToDo: Write error for suffix keys that aren't ParamData or VarData -# -# return param_list - - -# Adding utility to update parameter values in a model based on the suffix -def update_model_from_suffix(suffix_obj: pyo.Suffix, values): - """ - Overwrite each variable/parameter referenced by ``suffix_obj`` with the - corresponding value in ``values``. - - Parameters - ---------- - suffix_obj : pyomo.core.base.suffix.Suffix - The suffix whose *keys* are the components you want to update. - Call like ``update_from_suffix(model.unknown_parameters, vals)``. - values : iterable of numbers - New numerical values for the components referenced by the suffix. - Must be the same length as ``suffix_obj``. - """ - # Check that the length of values matches the suffix length - items = list(suffix_obj.items()) - if len(items) != len(values): - raise ValueError("values length does not match suffix length") - - # Iterate through the items in the suffix and update their values - # Note: items are tuples of (component, suffix_value) - for (comp, _), new_val in zip(items, values): - - # update the component value - # Check if the component is a VarData or ParamData - if isinstance(comp, (VarData, ParamData)): - comp.set_value(new_val) - else: - raise TypeError( - f"Unsupported component type {type(comp)}; expected VarData or ParamData." - ) - - def check_matrix(mat, check_pos_def=True): """ Checks that the matrix is square, positive definite, and symmetric. From 555a3078a326c858644d813967e04a42328cea7f Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 23 Sep 2025 10:59:29 -0400 Subject: [PATCH 24/25] changed the docstring --- pyomo/contrib/doe/doe.py | 47 +++++++++++++++++++++++++++------------- 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 0fb379b88d8..9bfb49166bc 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -1828,7 +1828,7 @@ def compute_FIM_factorial( file_name: str = None, ): """Will run a simulation-based factorial exploration of the experimental input - space (i.e., a ``grid search`` or ``parameter sweep``) to understand how the + space (a.k.a., a ``grid search`` or ``parameter sweep``) to understand how the FIM metrics change as a function of the experimental design space. This method can be used for both full factorial and fractional factorial designs. @@ -1840,15 +1840,18 @@ def compute_FIM_factorial( dict of lists or other array-like objects, of the form {"var_name": }. Default: None. The `design_values` should have the key(s) passed as strings that is a - subset of the `experiment_inputs`. If one or more design variables are not - to be changed, then they should not be passed in the `design_values` - dictionary, but if they are passed in the dictionary, then they must be a - array-like object of floats. For example, if our experiment has 4 design variables + subset of the `experiment_inputs` Suffix of the get_labeled_model(). If one + or more design variables are not to be changed, then they should not be + passed in the `design_values` dictionary, but if they are passed in the + dictionary, then they must be an array-like object of floats. For example, + if our experiment has 4 design variables (i.e., `experiment_inputs`): model.x1, model.x2, model.x3, and model.x4, their values may be passed as, design_values= {"x1": [1, 2, 3], "x3": [7], "x4": [-10, 20]}. In this case, x2 will not be changed and will be fixed at - the value in the model. - If design_values is provided, then `abs_change` and `rel_change` are ignored. + the value in the model. Of course, we can pass x2 here as well if we want to + change it. + If design_values argument is provided, then `abs_change`, `rel_change`, + `n_designs` arguments will be ignored. abs_change : list, optional Absolute change in the design variable values. Default: None. If provided, will use this value to generate the design values. @@ -1857,19 +1860,27 @@ def compute_FIM_factorial( Formula to calculate the design values: change_in_value = lower_bound * rel_change + abs_change` design_value += design_value + change_in_value + Also, the design variables need to be bounded for this argument to be used. rel_change : list, optional Relative change in the design variable values. Default: None. If provided, will use this value to generate the design values. If `rel_change` is provided, but `abs_change` is not provided, `abs_change` - will be set to zero. + will be set to zero. Also, the design variables need to be bounded for this + argument to be used. n_designs : int, optional - Number of designs to generate for each design variable. Default: 5. - If `abs_change` and/or `rel_change` are provided, this value will be ignored. + Number of design points to generate for each design variable. Default: 5. + if n_designs = 3 and we have design variables x1 ( 1 x1 = [1, 3, 5] + x2 = np.linspace(10, 20, 3) -> x2 = [10, 15, 20] + If `abs_change` and/or `rel_change` are provided, this argument will be + ignored. Also, the design variables need to be bounded for this argument + to be used. method : str, optional - string to specify which method should be used. options are ``kaug`` and - ``sequential`. Default: "sequential" + string to specify which method should be used to compute Fisher Information + Matrix (FIM). Options are ``kaug`` and ``sequential`. Default: "sequential" df_settings : tuple, optional - A tuple containing the settings for set_option() method of the pandas + A tuple containing the arguments for option_context() method of the pandas DataFrame. Default: (True, None, None, 500) - first element: whether to return a pandas DataFrame (True/False) - second element: number of max_columns for the DataFrame. Default: None, @@ -1877,12 +1888,18 @@ def compute_FIM_factorial( - third element: number of max_rows for the DataFrame. Default: None, i.e., no limit on the number of rows. - fourth element: display width for the DataFrame. Default: 500. + If the first element is True, the results (FIM metrics) will be printed as + a pandas DataFrame as well. initialization_scheme : str, optional Which scheme to use for initializing the design variables. Options are ``"snake_traversal"`` and ``"nested_for_loop"``. - Default: "snake_traversal" + If ``"snake_traversal"``, the design variables will be initialized in a + snake-like pattern where only one design variables change at a time. + If ``"nested_for_loop"``, the design variables will be initialized using a + nested for loop. Default: "snake_traversal" file_name : str, optional - if provided, will save the results to a json file. Default: None + if provided, will save the results to a json file. Default: None (will not + save the results). Returns ------- From 40ac36c16fd8c4ad4450e424a76acc2f42a301f9 Mon Sep 17 00:00:00 2001 From: Shuvashish Mondal Date: Tue, 23 Sep 2025 17:19:31 -0400 Subject: [PATCH 25/25] implemented suggestions from 9/23 pyomo dev meeting. --- pyomo/contrib/doe/doe.py | 135 +++++++++++++++++++++------------------ 1 file changed, 72 insertions(+), 63 deletions(-) diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 9bfb49166bc..9ee20854371 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -79,7 +79,7 @@ class FiniteDifferenceStep(Enum): backward = "backward" -class SensitivityInitialization(Enum): +class DesignSpaceTraversal(Enum): snake_traversal = "snake_traversal" nested_for_loop = "nested_for_loop" @@ -1821,10 +1821,10 @@ def compute_FIM_factorial( design_vals: dict = None, abs_change: list = None, rel_change: list = None, - n_designs: int = 5, + n_design_points: int = None, method="sequential", - df_settings=(True, None, None, 500), - initialization_scheme="snake_traversal", + return_df=True, + traversal_scheme="snake_traversal", file_name: str = None, ): """Will run a simulation-based factorial exploration of the experimental input @@ -1838,68 +1838,65 @@ def compute_FIM_factorial( The model to perform the full factorial exploration on. Default: None design_vals : dict, dict of lists or other array-like objects, of the form - {"var_name": }. Default: None. + {"var_name": []}. Default: None. The `design_values` should have the key(s) passed as strings that is a subset of the `experiment_inputs` Suffix of the get_labeled_model(). If one or more design variables are not to be changed, then they should not be passed in the `design_values` dictionary, but if they are passed in the dictionary, then they must be an array-like object of floats. For example, - if our experiment has 4 design variables - (i.e., `experiment_inputs`): model.x1, model.x2, model.x3, and model.x4, + if our experiment has 4 design variables (i.e., `experiment_inputs`): + model.x1, model.x2, model.x3, and model.x4, their values may be passed as, design_values= {"x1": [1, 2, 3], "x3": [7], "x4": [-10, 20]}. In this case, x2 will not be changed and will be fixed at the value in the model. Of course, we can pass x2 here as well if we want to change it. - If design_values argument is provided, then `abs_change`, `rel_change`, - `n_designs` arguments will be ignored. abs_change : list, optional Absolute change in the design variable values. Default: None. + Length of the list should be equal to the number of design variables in the + experiment_inputs Suffix of the get_labeled_model(), and the order of the values + in the list should be the same as the order of the design variables in the + experiment_inputs Suffix of the get_labeled_model(). If provided, will use this value to generate the design values. - If `abs_change` is provided, but `rel_change` is not provided, `rel_change` - will be set to zero. + If `abs_change` is provided, but `rel_change` is not provided, then + `rel_change` will be set to zero. Formula to calculate the design values: change_in_value = lower_bound * rel_change + abs_change` design_value += design_value + change_in_value Also, the design variables need to be bounded for this argument to be used. rel_change : list, optional Relative change in the design variable values. Default: None. + Length of the list should be equal to the number of design variables in the + experiment_inputs Suffix of the get_labeled_model(), and the order of the values + in the list should be the same as the order of the design variables in the + experiment_inputs Suffix of the get_labeled_model(). If provided, will use this value to generate the design values. If `rel_change` is provided, but `abs_change` is not provided, `abs_change` will be set to zero. Also, the design variables need to be bounded for this argument to be used. - n_designs : int, optional - Number of design points to generate for each design variable. Default: 5. - if n_designs = 3 and we have design variables x1 ( 1 x1 = [1, 3, 5] x2 = np.linspace(10, 20, 3) -> x2 = [10, 15, 20] - If `abs_change` and/or `rel_change` are provided, this argument will be - ignored. Also, the design variables need to be bounded for this argument - to be used. + Also, the design variables need to be bounded for this argument to be used. method : str, optional string to specify which method should be used to compute Fisher Information Matrix (FIM). Options are ``kaug`` and ``sequential`. Default: "sequential" - df_settings : tuple, optional - A tuple containing the arguments for option_context() method of the pandas - DataFrame. Default: (True, None, None, 500) - - first element: whether to return a pandas DataFrame (True/False) - - second element: number of max_columns for the DataFrame. Default: None, - i.e., no limit on the number of columns. - - third element: number of max_rows for the DataFrame. Default: None, - i.e., no limit on the number of rows. - - fourth element: display width for the DataFrame. Default: 500. - If the first element is True, the results (FIM metrics) will be printed as + return_df : bool, optional + If True, the results (FIM metrics) will be printed as a pandas DataFrame as well. - initialization_scheme : str, optional + traversal_scheme : str, optional Which scheme to use for initializing the design variables. Options are ``"snake_traversal"`` and ``"nested_for_loop"``. - If ``"snake_traversal"``, the design variables will be initialized in a + If ``"snake_traversal"`` is used, the design variables will be initialized in a snake-like pattern where only one design variables change at a time. - If ``"nested_for_loop"``, the design variables will be initialized using a + If ``"nested_for_loop"`` is used, the design variables will be initialized using a nested for loop. Default: "snake_traversal" file_name : str, optional - if provided, will save the results to a json file. Default: None (will not - save the results). + If provided, will save the results to a json file. + Example: file_name = "myresult". + Default: None (will not save the results). Returns ------- @@ -1933,6 +1930,29 @@ def compute_FIM_factorial( ).clone() model = self.factorial_model + # Group the mutually exclusive arguments for passing design values + # abs_change, and rel_change can be provided together or separately + num_des_provided = ( + (design_vals is not None) + + (n_design_points is not None) + + (abs_change is not None or rel_change is not None) + ) + + # Check the count and raise an error if more than one was provided + if num_des_provided > 1: + raise ValueError( + "Please provide only one of the following arguments: " + "`design_vals`, `n_design_points`, or `abs_change` and/or `rel_change`." + ) + + # Optional: Check if at least one was provided, if that's a requirement + if num_des_provided == 0: + raise ValueError( + "At least one of the following arguments: " + "`design_vals`, `n_design_points`, or `abs_change` and/or `rel_change`" + " must be provided." + ) + if design_vals is not None: # Check whether the design_ranges keys are in the experiment_inputs design_keys = set(design_vals.keys()) @@ -2002,6 +2022,7 @@ def compute_FIM_factorial( # change_in_value = lower_bound * rel_change + abs_change elif abs_change or rel_change: des_val = [] + # Calculate the change in value, delta value del_val = comp.lb * rel_change[i] + abs_change[i] if del_val == 0: raise ValueError( @@ -2023,20 +2044,20 @@ def compute_FIM_factorial( # generate the factorial points based on the initialization scheme try: - scheme_enum = SensitivityInitialization(initialization_scheme) + scheme_enum = DesignSpaceTraversal(traversal_scheme) except ValueError: self.logger.warning( - f"Initialization scheme '{initialization_scheme}' is not recognized. " + f"Initialization scheme '{traversal_scheme}' is not recognized. " "Using `snake_traversal` as the default initialization scheme." ) - scheme_enum = SensitivityInitialization.snake_traversal + scheme_enum = DesignSpaceTraversal.snake_traversal - if scheme_enum == SensitivityInitialization.snake_traversal: + if scheme_enum == DesignSpaceTraversal.snake_traversal: factorial_points = snake_traversal_grid_sampling(*design_values) - elif scheme_enum == SensitivityInitialization.nested_for_loop: + elif scheme_enum == DesignSpaceTraversal.nested_for_loop: factorial_points = product(*design_values) - # TODO: Add more initialization schemes + # TODO: Add more DesignSpaceTraversal schemes factorial_points_list = list(factorial_points) @@ -2142,31 +2163,19 @@ def compute_FIM_factorial( self.factorial_results = factorial_results # Set pandas DataFrame options - if df_settings[0]: - with pd.option_context( - "display.max_columns", - df_settings[1], - "display.max_rows", - df_settings[2], - "display.width", - df_settings[3], - ): - exclude_keys = { - "total_points", - "success_count", - "failure_count", - "FIM_all", - } - dict_for_df = { - k: v for k, v in factorial_results.items() if k not in exclude_keys - } - res_df = pd.DataFrame(dict_for_df) - print("\n\n=========Factorial results===========") - print("Total points:", total_points) - print("Success count:", success_count) - print("Failure count:", failure_count) - print("\n") - print(res_df) + if return_df: + # exclude matrix-like entries, and unchanging results from the DataFrame + exclude_keys = {"total_points", "success_count", "failure_count", "FIM_all"} + dict_for_df = { + k: v for k, v in factorial_results.items() if k not in exclude_keys + } + res_df = pd.DataFrame(dict_for_df) + print("\n\n=========Factorial results===========") + print("Total points:", total_points) + print("Success count:", success_count) + print("Failure count:", failure_count) + print("\n") + print(res_df) # Save the results to a json file based on the file_name provided if file_name is not None: