diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index 8c5f5f86ada..103bdc96225 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -25,11 +25,11 @@ # ____________________________________________________________________________________ from enum import Enum -from itertools import permutations, product - +from itertools import permutations, product, combinations import json import logging import math +import warnings from pyomo.common.dependencies import ( numpy as np, @@ -37,6 +37,7 @@ pandas as pd, pathlib, matplotlib as plt, + scipy, scipy_available, ) @@ -45,11 +46,6 @@ from pyomo.contrib.sensitivity_toolbox.sens import get_dsdp -if numpy_available and scipy_available: - from pyomo.contrib.doe.grey_box_utilities import FIMExternalGreyBox - - from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock - import pyomo.environ as pyo from pyomo.contrib.doe.utils import ( _SMALL_TOLERANCE_DEFINITENESS, @@ -57,6 +53,7 @@ compute_FIM_metrics, regularize_fim_for_cholesky, ) +from pyomo.contrib.parmest.utils.model_utils import update_model_from_suffix from pyomo.opt import SolverStatus @@ -76,10 +73,37 @@ class FiniteDifferenceStep(Enum): backward = "backward" +class InitializationMethod(Enum): + latin_hypercube_sampling = "latin_hypercube_sampling" + + +class _DoEResultsJSONEncoder(json.JSONEncoder): + """JSON encoder for DoE result payloads with numpy/Pyomo objects.""" + + def default(self, obj): + if isinstance(obj, np.generic): + return obj.item() + if isinstance(obj, np.ndarray): + return obj.tolist() + if isinstance(obj, Enum): + return str(obj) + return super().default(obj) + + class DesignOfExperiments: + # Objective options whose scalar score is compared with "larger is better" + # in initialization and diagnostics paths. + _MAXIMIZE_OBJECTIVES = frozenset( + { + ObjectiveLib.determinant, + ObjectiveLib.pseudo_trace, + ObjectiveLib.minimum_eigenvalue, + } + ) + def __init__( self, - experiment=None, + experiment: list = None, fd_formula="central", step=1e-3, objective_option="determinant", @@ -110,13 +134,17 @@ def __init__( Parameters ---------- experiment: - Experiment object that holds the model and labels all the components. The - object should have a ``get_labeled_model`` where a model is returned with - the following labeled sets: + Experiment object(s) that hold the model and labels all the components. + Can be a single Experiment object or a list of Experiment objects. + For single experiments, you can pass the object directly: + ``experiment=exp`` or as a list: ``experiment=[exp]``. + Each object should have a ``get_labeled_model`` method that returns a model with + the following labeled Pyomo Suffixes: - ``unknown_parameters``, - ``experimental_inputs``, - ``experimental_outputs`` + - ``measurement_error``. fd_formula: Finite difference formula for computing the sensitivity matrix. Must be @@ -130,7 +158,8 @@ def __init__( - ``determinant`` (for determinant, or D-optimality), - ``trace`` (for trace of covariance matrix, or A-optimality), - ``pseudo_trace`` (for trace of Fisher Information Matrix(FIM), or pseudo A-optimality), - - ``minimum_eigenvalue``, (for E-optimality), or ``condition_number`` (for ME-optimality) + - ``minimum_eigenvalue``, (for E-optimality), or + - ``condition_number`` (for ME-optimality) Note: E-optimality and ME-optimality are only supported when using the grey box objective (i.e., ``grey_box_solver`` is True) default: ``determinant`` @@ -187,17 +216,30 @@ def __init__( Specify the level of the logger. Change to logging.DEBUG for all messages. """ - if experiment is None: - raise ValueError("Experiment object must be provided to perform DoE.") - - # Check if the Experiment object has callable ``get_labeled_model`` function - if not hasattr(experiment, "get_labeled_model"): + # Validate experiment(s) are provided + if experiment is None or (isinstance(experiment, list) and not experiment): raise ValueError( - "The experiment object must have a ``get_labeled_model`` function" + "The 'experiment' argument is required and cannot be an empty list. " + "Pass a single Experiment object or a non-empty list of Experiment " + "objects." ) - # Set the experiment object from the user - self.experiment = experiment + # Auto-convert single experiment to list + if not isinstance(experiment, list): + experiment_list = [experiment] + else: + experiment_list = experiment + + # Check each experiment has get_labeled_model method + for idx, exp in enumerate(experiment_list): + if not hasattr(exp, "get_labeled_model"): + raise ValueError( + f"Experiment at index {idx} in 'experiment' must have a " + f"'get_labeled_model' method" + ) + + # Store experiment_list + self.experiment_list = experiment_list # Set the finite difference and subsequent step size self.fd_formula = FiniteDifferenceStep(fd_formula) @@ -270,6 +312,67 @@ def __init__( # (i.e., no model rebuilding for large models with sequential) self._built_scenarios = False + @staticmethod + def _enum_label(value): + """Return a stable short label for enum-like values.""" + return str(value).split(".")[-1] + + def _grey_box_output_name(self): + """Return the output label exposed by the FIM grey box model.""" + if self.objective_option == ObjectiveLib.trace: + return "A-opt" + if self.objective_option == ObjectiveLib.pseudo_trace: + return "pseudo-A-opt" + if self.objective_option == ObjectiveLib.determinant: + return "log-D-opt" + if self.objective_option == ObjectiveLib.minimum_eigenvalue: + return "E-opt" + if self.objective_option == ObjectiveLib.condition_number: + return "ME-opt" + raise ValueError( + "Grey-box objective support is only available for " + "objective_option in ['determinant', 'trace', 'pseudo_trace', " + "'minimum_eigenvalue', 'condition_number']." + ) + + def _initialize_grey_box_block(self, egb_block, fim_np, parameter_names): + """ + Seed a grey box block from the FIM computed by the square solve. + + The square solve initializes the finite-difference scenarios and the + aggregated FIM, but the external grey box block stores that FIM through + its own input/output variables. We therefore push the current FIM values + into the grey box explicitly before the final NLP solve so the external + model starts from a consistent state. + """ + param_list = list(parameter_names) + + # The external model only takes the symmetric triangular FIM entries as + # inputs. That keeps the grey box interface compact while still + # reconstructing the full symmetric FIM internally. + for i, p1 in enumerate(param_list): + for j, p2 in enumerate(param_list): + if i >= j: + egb_block.inputs[(p2, p1)].set_value(float(fim_np[i, j])) + + # Initialize the single grey box output so the final solve begins with + # an objective value consistent with the current square-solve FIM. + if self.objective_option == ObjectiveLib.trace: + output_value = np.trace(np.linalg.pinv(fim_np)) + elif self.objective_option == ObjectiveLib.pseudo_trace: + output_value = np.trace(fim_np) + elif self.objective_option == ObjectiveLib.determinant: + output_value = np.log(np.linalg.det(fim_np)) + elif self.objective_option == ObjectiveLib.minimum_eigenvalue: + output_value = np.min(np.linalg.eigvalsh(fim_np)) + elif self.objective_option == ObjectiveLib.condition_number: + eig = np.linalg.eigvalsh(fim_np) + output_value = np.log(np.abs(np.max(eig) / np.min(eig))) + else: + output_value = 0.0 + + egb_block.outputs[self._grey_box_output_name()].set_value(float(output_value)) + # Perform doe def run_doe(self, model=None, results_file=None): """ @@ -285,11 +388,10 @@ def run_doe(self, model=None, results_file=None): """ # Check results file name if results_file is not None: - if type(results_file) not in [pathlib.Path, str]: + if not isinstance(results_file, (pathlib.Path, str)): raise ValueError( "``results_file`` must be either a Path object or a string." ) - # Start timer sp_timer = TicTocTimer() sp_timer.tic(msg=None) @@ -331,7 +433,7 @@ def run_doe(self, model=None, results_file=None): # and fix the design variables. model.objective.deactivate() model.obj_cons.deactivate() - for comp in model.scenario_blocks[0].experiment_inputs: + for comp in model.fd_scenario_blocks[0].experiment_inputs: comp.fix() # TODO: safeguard solver call to see if solver terminated successfully @@ -357,7 +459,7 @@ def run_doe(self, model=None, results_file=None): model.dummy_obj.deactivate() # Reactivate objective and unfix experimental design decisions - for comp in model.scenario_blocks[0].experiment_inputs: + for comp in model.fd_scenario_blocks[0].experiment_inputs: comp.unfix() model.objective.activate() model.obj_cons.activate() @@ -376,11 +478,6 @@ def run_doe(self, model=None, results_file=None): if self.objective_option == ObjectiveLib.trace: trace_val = np.trace(np.linalg.pinv(self.get_FIM())) model.obj_cons.egb_fim_block.outputs["A-opt"].set_value(trace_val) - elif self.objective_option == ObjectiveLib.pseudo_trace: - pseudo_trace_val = np.trace(np.array(self.get_FIM())) - model.obj_cons.egb_fim_block.outputs["pseudo-A-opt"].set_value( - pseudo_trace_val - ) elif self.objective_option == ObjectiveLib.determinant: det_val = np.linalg.det(np.array(self.get_FIM())) model.obj_cons.egb_fim_block.outputs["log-D-opt"].set_value( @@ -394,8 +491,63 @@ def run_doe(self, model=None, results_file=None): cond_number = np.log(np.abs(np.max(eig) / np.min(eig))) model.obj_cons.egb_fim_block.outputs["ME-opt"].set_value(cond_number) - # Keep Cholesky-related variables synchronized with current FIM values - self._initialize_cholesky_from_fim(model=model) + # If the model has L, initialize it with the solved FIM + if hasattr(model, "L"): + # Get the FIM values + fim_vals = [ + pyo.value(model.fim[i, j]) + for i in model.parameter_names + for j in model.parameter_names + ] + fim_np = np.array(fim_vals).reshape( + (len(model.parameter_names), len(model.parameter_names)) + ) + + # Need to compute the full FIM before + # initializing the Cholesky factorization + if self.only_compute_fim_lower: + fim_np = fim_np + fim_np.T - np.diag(np.diag(fim_np)) + + # Check if the FIM is positive definite + # If not, add jitter to the diagonal + # to ensure positive definiteness + min_eig = np.min(np.linalg.eigvals(fim_np)) + + if min_eig < _SMALL_TOLERANCE_DEFINITENESS: + # Raise the minimum eigenvalue to at + # least _SMALL_TOLERANCE_DEFINITENESS + jitter = np.min( + [ + -min_eig + _SMALL_TOLERANCE_DEFINITENESS, + _SMALL_TOLERANCE_DEFINITENESS, + ] + ) + else: + # No jitter needed + jitter = 0 + + # Add jitter to the diagonal to ensure positive definiteness + L_vals_sq = np.linalg.cholesky( + fim_np + jitter * np.eye(len(model.parameter_names)) + ) + for i, c in enumerate(model.parameter_names): + for j, d in enumerate(model.parameter_names): + model.L[c, d].value = L_vals_sq[i, j] + + # Initialize the inverse of L if it exists + if hasattr(model, "L_inv"): + L_inv_vals = np.linalg.inv(L_vals_sq) + + for i, c in enumerate(model.parameter_names): + for j, d in enumerate(model.parameter_names): + if i >= j: + model.L_inv[c, d].value = L_inv_vals[i, j] + else: + model.L_inv[c, d].value = 0.0 + # Initialize the cov_trace if it exists + if hasattr(model, "cov_trace"): + initial_cov_trace = np.sum(L_inv_vals**2) + model.cov_trace.value = initial_cov_trace if hasattr(model, "determinant"): model.determinant.value = np.linalg.det(np.array(self.get_FIM())) @@ -428,10 +580,14 @@ def run_doe(self, model=None, results_file=None): self.results["Solver Status"] = res.solver.status self.results["Termination Condition"] = res.solver.termination_condition - if type(res.solver.message) is str: + if isinstance(res.solver.message, str): results_message = res.solver.message - elif type(res.solver.message) is bytes: + elif isinstance(res.solver.message, bytes): results_message = res.solver.message.decode("utf-8") + else: + results_message = ( + str(res.solver.message) if res.solver.message is not None else "" + ) self.results["Termination Message"] = results_message # Important quantities for optimal design @@ -439,29 +595,29 @@ def run_doe(self, model=None, results_file=None): self.results["Sensitivity Matrix"] = self.get_sensitivity_matrix() self.results["Experiment Design"] = self.get_experiment_input_values() self.results["Experiment Design Names"] = [ - str(pyo.ComponentUID(k, context=model.scenario_blocks[0])) - for k in model.scenario_blocks[0].experiment_inputs + str(pyo.ComponentUID(k, context=model.fd_scenario_blocks[0])) + for k in model.fd_scenario_blocks[0].experiment_inputs ] self.results["Experiment Outputs"] = self.get_experiment_output_values() self.results["Experiment Output Names"] = [ - str(pyo.ComponentUID(k, context=model.scenario_blocks[0])) - for k in model.scenario_blocks[0].experiment_outputs + str(pyo.ComponentUID(k, context=model.fd_scenario_blocks[0])) + for k in model.fd_scenario_blocks[0].experiment_outputs ] self.results["Unknown Parameters"] = self.get_unknown_parameter_values() self.results["Unknown Parameter Names"] = [ - str(pyo.ComponentUID(k, context=model.scenario_blocks[0])) - for k in model.scenario_blocks[0].unknown_parameters + str(pyo.ComponentUID(k, context=model.fd_scenario_blocks[0])) + for k in model.fd_scenario_blocks[0].unknown_parameters ] self.results["Measurement Error"] = self.get_measurement_error_values() self.results["Measurement Error Names"] = [ - str(pyo.ComponentUID(k, context=model.scenario_blocks[0])) - for k in model.scenario_blocks[0].measurement_error + str(pyo.ComponentUID(k, context=model.fd_scenario_blocks[0])) + for k in model.fd_scenario_blocks[0].measurement_error ] - self.results["Prior FIM"] = [list(row) for row in list(self.prior_FIM)] + self.results["Prior FIM"] = self.prior_FIM.tolist() # Saving some stats on the FIM for convenience - self.results["Objective expression"] = str(self.objective_option).split(".")[-1] + self.results["Objective expression"] = self._enum_label(self.objective_option) self.results["log10 A-opt"] = np.log10(np.trace(np.linalg.inv(fim_local))) self.results["log10 pseudo A-opt"] = np.log10(np.trace(fim_local)) self.results["log10 D-opt"] = np.log10(np.linalg.det(fim_local)) @@ -475,7 +631,7 @@ def run_doe(self, model=None, results_file=None): self.results["Wall-clock Time"] = build_time + initialization_time + solve_time # Settings used to generate the optimal DoE - self.results["Finite Difference Scheme"] = str(self.fd_formula).split(".")[-1] + self.results["Finite Difference Scheme"] = self._enum_label(self.fd_formula) self.results["Finite Difference Step"] = self.step self.results["Nominal Parameter Scaling"] = self.scale_nominal_param_value @@ -536,42 +692,1266 @@ def _initialize_cholesky_from_fim(self, model=None): # This happens if the function is called with a model using GreyBox. return - fim_np = self._get_fim_numpy(model) - fim_pd, _ = regularize_fim_for_cholesky(fim_np) + fim_np = self._get_fim_numpy(model) + fim_pd, _ = regularize_fim_for_cholesky(fim_np) + + L_vals = np.linalg.cholesky(fim_pd) + for i, c in enumerate(model.parameter_names): + for j, d in enumerate(model.parameter_names): + if i >= j: + model.L[c, d].value = L_vals[i, j] + else: + model.L[c, d].value = 0.0 + + if hasattr(model, "L_inv"): + L_inv_vals = np.linalg.inv(L_vals) + for i, c in enumerate(model.parameter_names): + for j, d in enumerate(model.parameter_names): + if i >= j: + model.L_inv[c, d].value = L_inv_vals[i, j] + else: + model.L_inv[c, d].value = 0.0 + + if hasattr(model, "fim_inv"): + # Use the pseudo-inverse here rather than the strict inverse. + # The jittered matrix should be positive definite, but ``pinv`` + # is safer for borderline ill-conditioned cases and matches the + # defensive approach already used when initializing ``fim_inv`` + # from user-provided starting values. + fim_inv_vals = np.linalg.pinv(fim_pd) + for i, c in enumerate(model.parameter_names): + for j, d in enumerate(model.parameter_names): + if self.only_compute_fim_lower and i < j: + model.fim_inv[c, d].value = 0.0 + else: + model.fim_inv[c, d].value = fim_inv_vals[i, j] + + if hasattr(model, "cov_trace"): + model.cov_trace.value = np.trace(fim_inv_vals) + + def optimize_experiments( + self, + results_file=None, + n_exp: int = None, + init_method: InitializationMethod = None, + init_n_samples: int = 5, + init_seed: int = None, + init_solver=None, + ): + """ + Optimize single experiment or multiple experiments simultaneously for + Design of Experiments. + + The number of experiments is determined by the length of the + experiment list provided through the ``experiment`` argument. + + Parameters + ---------- + results_file: + string name of the file path to save the results to in the form + of a .json file + + init_method: + Method used to initialize the experiment design variables before + optimization. Options are: + + - ``None`` (default): No special initialization; use the initial + values from ``get_labeled_model()``. To provide a custom starting + point, initialize the ``Experiment`` objects with the desired + design values before passing them in ``experiment``. + - ``"latin_hypercube_sampling"`` (or ``InitializationMethod.latin_hypercube_sampling``): Use Latin Hypercube Sampling (LHS) to find a good + initial design. For each experiment-input dimension, ``init_n_samples`` + points are sampled independently using 1-D LHS, and their Cartesian + product forms the set of candidate experiment designs. The FIM is + evaluated at every candidate, and the combination of ``n_exp`` + candidates (without replacement) that best satisfies the chosen + objective is selected as the starting point for the optimization. + + init_n_samples: + Number of LHS samples per experiment-input dimension when + ``init_method="latin_hypercube_sampling"``. The total number of candidate + designs is ``init_n_samples ** n_exp_inputs``. A warning is issued + when this exceeds 10,000. Default: 5. + + init_seed: + Integer seed for the LHS random-number generator (for + reproducibility). Used only when ``init_method="latin_hypercube_sampling"``. + Default: ``None`` (non-deterministic). + init_solver: + Optional solver object used only for initialization phases + (including multi-experiment block-construction solves, LHS + candidate-FIM evaluations, and the square initialization solve). + The final optimization solve always uses the primary DoE solver + (``self.solver``). If ``init_solver = None``, initialization also uses + ``self.solver``. + + Notes + ----- + Number of Experiments: + When ``len(experiment) == 1`` (template mode), pass ``n_exp`` + to specify how many experiments to optimize. When + ``len(experiment) > 1`` (user-initialized mode), the list + length determines the number of experiments and ``n_exp`` must + not be set. + + Symmetry Breaking (for multiple experiments): + To prevent equivalent permutations of identical experiments, you must + mark a "primary" design variable using a Pyomo Suffix in your experiment's + `label_experiment()` method: + + Example:: + + m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.sym_break_cons[m.CA[0]] = None # Mark CA[0] as primary variable + + This will add constraints: exp[k-1].primary_var <= exp[k].primary_var + for k = 1, ..., n_exp-1, which breaks permutation symmetry and can + significantly reduce solve times. + + LHS Initialization (init_method="latin_hypercube_sampling"): + Each dimension of the experiment inputs is sampled independently + using a 1-D Latin Hypercube, giving ``init_n_samples`` evenly-spaced + stratified samples across the variable bounds. The joint candidate + set is the Cartesian product of these per-dimension samples (i.e., + a ``init_n_samples^n_inputs`` grid with good marginal coverage). The + FIM is evaluated sequentially at each candidate, then all + ``C(n_candidates, n_exp)`` combinations are scored and the best one + is used as the initial point for the NLP solver. This can + significantly improve solution quality when the problem has multiple + local optima. + + """ + if results_file is not None and not isinstance( + results_file, (pathlib.Path, str) + ): + raise ValueError( + "``results_file`` must be either a Path object or a string." + ) + + n_exp, _template_mode = self._resolve_n_exp_and_mode(n_exp=n_exp) + resolved_init_method, resolved_init_solver = ( + self._resolve_optimize_experiments_initialization( + init_method=init_method, + init_n_samples=init_n_samples, + init_seed=init_seed, + init_solver=init_solver, + template_mode=_template_mode, + ) + ) + init_solver_name = getattr( + resolved_init_solver, "name", str(resolved_init_solver) + ) + + sp_timer = TicTocTimer() + sp_timer.tic(msg=None) + self.logger.info( + f"Beginning multi-experiment optimization with {n_exp} experiments." + ) + + self.model = pyo.ConcreteModel() + n_param_scenarios = 1 + self.scenario_weights = (1.0,) + self.model.param_scenario_blocks = pyo.Block(range(n_param_scenarios)) + + symmetry_breaking_info = { + "enabled": n_exp > 1, + "variable": None, + "source": None, + } + diagnostics_warnings = [] + lhs_init_diagnostics = None + lhs_initialization_time = 0.0 + best_initial_points = None + + primary_solver = self.solver + self.solver = resolved_init_solver + try: + self._build_multi_experiment_blocks( + n_param_scenarios=n_param_scenarios, + n_exp=n_exp, + template_mode=_template_mode, + ) + self._add_symmetry_breaking_constraints( + n_param_scenarios=n_param_scenarios, + n_exp=n_exp, + symmetry_breaking_info=symmetry_breaking_info, + diagnostics_warnings=diagnostics_warnings, + ) + + self.create_multi_experiment_objective_function(self.model) + build_time = sp_timer.toc(msg=None) + self.logger.info( + "Successfully built the multi-experiment DoE model.\nBuild time: %0.1f seconds" + % build_time + ) + + best_initial_points, lhs_init_diagnostics, lhs_initialization_time = ( + self._apply_optimize_experiment_initialization( + resolved_init_method=resolved_init_method, + init_n_samples=init_n_samples, + init_seed=init_seed, + n_param_scenarios=n_param_scenarios, + n_exp=n_exp, + ) + ) + + sp_timer.tic(msg=None) + initialization_time = self._run_square_initialization_solve( + n_param_scenarios=n_param_scenarios, n_exp=n_exp, timer=sp_timer + ) + + parameter_names = ( + self.model.param_scenario_blocks[0].exp_blocks[0].parameter_names + ) + self._initialize_scenario_quantities_from_square_solve( + n_param_scenarios=n_param_scenarios, + n_exp=n_exp, + parameter_names=parameter_names, + ) + finally: + self.solver = primary_solver + + final_solver = self.grey_box_solver if self.use_grey_box else self.solver + final_solver_name = getattr(final_solver, "name", str(final_solver)) + res = final_solver.solve( + self.model, tee=self.grey_box_tee if self.use_grey_box else self.tee + ) + + solve_time = sp_timer.toc(msg=None) + self.logger.info( + ( + "Successfully optimized multi-experiment design." + "\nSolve time: %0.1f seconds" % solve_time + ) + ) + self.logger.info( + "Total time for build, initialization, and solve: %0.1f seconds" + % (build_time + initialization_time + solve_time) + ) + + self.results = self._build_optimize_experiments_results( + res=res, + n_param_scenarios=n_param_scenarios, + n_exp=n_exp, + parameter_names=parameter_names, + template_mode=_template_mode, + resolved_init_method=resolved_init_method, + init_solver_name=init_solver_name, + init_n_samples=init_n_samples, + init_seed=init_seed, + best_initial_points=best_initial_points, + lhs_init_diagnostics=lhs_init_diagnostics, + lhs_initialization_time=lhs_initialization_time, + symmetry_breaking_info=symmetry_breaking_info, + diagnostics_warnings=diagnostics_warnings, + final_solver_name=final_solver_name, + build_time=build_time, + initialization_time=initialization_time, + solve_time=solve_time, + ) + + if results_file is not None: + with open(results_file, "w") as file: + json.dump(self.results, file, indent=2, cls=_DoEResultsJSONEncoder) + + def _resolve_n_exp_and_mode(self, n_exp): + """Resolve experiment count and whether optimize_experiments() is in template mode.""" + n_list = len(self.experiment_list) + if n_list > 1: + if n_exp is not None: + raise ValueError( + "``n_exp`` must not be set when the experiment list contains " + f"more than one experiment (got {n_list} experiments in the " + "list). Either pass a single template experiment and set " + "``n_exp``, or pass a fully-initialized list and omit ``n_exp``." + ) + return n_list, False + + if n_exp is None: + return 1, True + if not isinstance(n_exp, int) or n_exp < 1: + raise ValueError(f"``n_exp`` must be a positive integer, got {n_exp!r}.") + return n_exp, True + + def _resolve_optimize_experiments_initialization( + self, init_method, init_n_samples, init_seed, init_solver, template_mode + ): + """Validate and normalize optimize_experiments() initialization options.""" + if init_method is None: + resolved_init_method = None + else: + try: + resolved_init_method = InitializationMethod(init_method) + except ValueError: + valid = ", ".join(f"'{m.value}'" for m in InitializationMethod) + raise ValueError( + "``init_method`` must be one of [None, " + + valid + + f"], got {init_method!r}." + ) + + if resolved_init_method == InitializationMethod.latin_hypercube_sampling: + if not template_mode: + raise ValueError( + "``init_method='latin_hypercube_sampling'`` is currently supported only in " + "template mode (``len(experiment) == 1``)." + ) + if not scipy_available: + raise ImportError( + "LHS initialization requires scipy. " + "Please install scipy to use init_method='latin_hypercube_sampling'." + ) + if not isinstance(init_n_samples, int) or init_n_samples < 1: + raise ValueError( + "``init_n_samples`` must be a positive integer, " + f"got {init_n_samples!r}." + ) + if init_seed is not None and not isinstance(init_seed, int): + raise ValueError( + "``init_seed`` must be None or an integer, " f"got {init_seed!r}." + ) + + if init_solver is not None and not hasattr(init_solver, "solve"): + raise ValueError( + "``init_solver`` must be None or a solver object with a 'solve' method." + ) + + primary_solver = self.solver + resolved_init_solver = primary_solver if init_solver is None else init_solver + return resolved_init_method, resolved_init_solver + + def _build_multi_experiment_blocks(self, n_param_scenarios, n_exp, template_mode): + """Build scenario/experiment blocks and validate parameter alignment.""" + for s in range(n_param_scenarios): + scenario_block = self.model.param_scenario_blocks[s] + scenario_block.exp_blocks = pyo.Block(range(n_exp)) + reference_param_names = None + reference_param_values = None + + for k in range(n_exp): + self.create_doe_model( + model=scenario_block.exp_blocks[k], + experiment_index=0 if template_mode else k, + _for_multi_experiment=True, + ) + + if reference_param_names is None: + reference_fd_block = scenario_block.exp_blocks[ + k + ].fd_scenario_blocks[0] + reference_param_names = [ + str(pyo.ComponentUID(param, context=reference_fd_block)) + for param in reference_fd_block.unknown_parameters + ] + reference_param_values = [ + float(value) + for value in reference_fd_block.unknown_parameters.values() + ] + continue + + current_fd_block = scenario_block.exp_blocks[k].fd_scenario_blocks[0] + current_param_names = [ + str(pyo.ComponentUID(param, context=current_fd_block)) + for param in current_fd_block.unknown_parameters + ] + if current_param_names != reference_param_names: + raise ValueError( + "All experiments passed to optimize_experiments() " + "must define the same unknown_parameters in the same " + "order. " + f"Experiment 0 uses {reference_param_names}, while " + f"experiment {k} uses {current_param_names}." + ) + + current_param_values = [ + float(value) + for value in current_fd_block.unknown_parameters.values() + ] + for name, ref_val, cur_val in zip( + reference_param_names, reference_param_values, current_param_values + ): + if not math.isclose(ref_val, cur_val, rel_tol=1e-12, abs_tol=1e-12): + raise ValueError( + "All experiments passed to optimize_experiments() " + "must use the same nominal values for " + "unknown_parameters. " + f"Parameter '{name}' has value {ref_val} in " + f"experiment 0 and {cur_val} in experiment {k}." + ) + + def _add_symmetry_breaking_constraints( + self, n_param_scenarios, n_exp, symmetry_breaking_info, diagnostics_warnings + ): + """Add optional experiment-ordering constraints for permutation symmetry.""" + if n_exp <= 1: + return + + first_exp_block = ( + self.model.param_scenario_blocks[0].exp_blocks[0].fd_scenario_blocks[0] + ) + + if ( + hasattr(first_exp_block, 'sym_break_cons') + and len(first_exp_block.sym_break_cons) > 0 + ): + sym_break_var_list = list(first_exp_block.sym_break_cons.keys()) + if len(sym_break_var_list) > 1: + warning_msg = ( + "Multiple variables marked in sym_break_cons. " + f"Using {sym_break_var_list[0].local_name} for symmetry breaking." + ) + self.logger.warning(warning_msg) + diagnostics_warnings.append(warning_msg) + + sym_break_var = sym_break_var_list[0] + if not any( + sym_break_var is inp for inp in first_exp_block.experiment_inputs.keys() + ): + raise ValueError( + "Variable selected in ``sym_break_cons`` must also be an " + "experiment input variable. " + f"Got non-input variable '{sym_break_var.local_name}'." + ) + symmetry_breaking_info["variable"] = sym_break_var.local_name + symmetry_breaking_info["source"] = "user" + self.logger.info( + f"Using user-specified variable '{sym_break_var.local_name}' for symmetry breaking." + ) + else: + sym_break_var = next(iter(first_exp_block.experiment_inputs)) + symmetry_breaking_info["variable"] = sym_break_var.local_name + symmetry_breaking_info["source"] = "auto" + self.logger.warning( + "No symmetry breaking variable specified. Automatically using the first " + f"experiment input '{sym_break_var.local_name}' for ordering constraints. " + "To specify a different variable, add: " + "m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL); " + "m.sym_break_cons[m.your_variable] = None" + ) + diagnostics_warnings.append( + "No symmetry breaking variable specified. Automatically using " + f"'{sym_break_var.local_name}'." + ) + + for s in range(n_param_scenarios): + for k in range(1, n_exp): + var_prev = pyo.ComponentUID( + sym_break_var, context=first_exp_block + ).find_component_on( + self.model.param_scenario_blocks[s] + .exp_blocks[k - 1] + .fd_scenario_blocks[0] + ) + var_curr = pyo.ComponentUID( + sym_break_var, context=first_exp_block + ).find_component_on( + self.model.param_scenario_blocks[s] + .exp_blocks[k] + .fd_scenario_blocks[0] + ) + if var_prev is None or var_curr is None: + raise RuntimeError( + "Failed to map symmetry breaking variable " + f"'{sym_break_var.local_name}' onto scenario {s}, " + f"experiment pair ({k - 1}, {k}). Ensure the variable " + "exists on all experiment blocks with compatible labels." + ) + + con_name = f"symmetry_breaking_s{s}_exp{k}" + self.model.param_scenario_blocks[s].add_component( + con_name, pyo.Constraint(expr=var_prev <= var_curr) + ) + + self.logger.info( + f"Added {n_exp - 1} symmetry breaking constraints for scenario {s} " + f"using variable: {sym_break_var.local_name}" + ) + + def _apply_optimize_experiment_initialization( + self, resolved_init_method, init_n_samples, init_seed, n_param_scenarios, n_exp + ): + """Apply optional initialization strategy and return diagnostics.""" + if resolved_init_method != InitializationMethod.latin_hypercube_sampling: + return None, None, 0.0 + + lhs_timer = TicTocTimer() + lhs_timer.tic(msg=None) + self.logger.info( + f"Applying LHS initialization with {init_n_samples} samples per " + f"experiment-input dimension..." + ) + best_initial_points, lhs_init_diagnostics = self._lhs_initialize_experiments( + lhs_n_samples=init_n_samples, lhs_seed=init_seed, n_exp=n_exp + ) + self.logger.info( + "Setting LHS best-found initial design in the optimization model..." + ) + for s in range(n_param_scenarios): + for k in range(n_exp): + exp_input_vars = self._get_experiment_input_vars( + self.model.param_scenario_blocks[s].exp_blocks[k] + ) + for var, val in zip(exp_input_vars, best_initial_points[k]): + var.set_value(val) + lhs_initialization_time = lhs_timer.toc(msg=None) + return best_initial_points, lhs_init_diagnostics, lhs_initialization_time + + def _run_square_initialization_solve(self, n_param_scenarios, n_exp, timer): + """Run objective-free square solve used to initialize scenario quantities.""" + self.model.objective.deactivate() + for s in range(n_param_scenarios): + if hasattr(self.model.param_scenario_blocks[s], "obj_cons"): + self.model.param_scenario_blocks[s].obj_cons.deactivate() + + self._set_experiment_inputs_fixed( + n_param_scenarios=n_param_scenarios, n_exp=n_exp, fixed=True + ) + + self.model.dummy_obj = pyo.Objective(expr=0, sense=pyo.minimize) + self.solver.solve(self.model, tee=self.tee) + initialization_time = timer.toc(msg=None) + + self.logger.info( + "Successfully initialized the multi-experiment DoE model." + "\nInitialization time: %0.1f seconds" % initialization_time + ) + + self.model.dummy_obj.deactivate() + self.model.del_component("dummy_obj") + self._set_experiment_inputs_fixed( + n_param_scenarios=n_param_scenarios, n_exp=n_exp, fixed=False + ) + self.model.objective.activate() + for s in range(n_param_scenarios): + if hasattr(self.model.param_scenario_blocks[s], "obj_cons"): + self.model.param_scenario_blocks[s].obj_cons.activate() + + return initialization_time + + def _initialize_scenario_quantities_from_square_solve( + self, n_param_scenarios, n_exp, parameter_names + ): + """Initialize objective-specific scenario variables from square-solve FIMs.""" + for s in range(n_param_scenarios): + scenario = self.model.param_scenario_blocks[s] + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if self.only_compute_fim_lower and i < j: + continue + fim_sum = sum( + pyo.value(scenario.exp_blocks[k].fim[p, q]) or 0 + for k in range(n_exp) + ) + scenario.total_fim[p, q].set_value(fim_sum + self.prior_FIM[i, j]) + + total_fim_vals = [ + pyo.value(scenario.total_fim[p, q]) + for p in parameter_names + for q in parameter_names + ] + total_fim_np = np.array(total_fim_vals).reshape( + (len(parameter_names), len(parameter_names)) + ) + if self.only_compute_fim_lower: + total_fim_np = self._symmetrize_lower_tri(total_fim_np) + obj_cons = getattr(scenario, "obj_cons", None) + + if self.use_grey_box: + self._initialize_grey_box_block( + obj_cons.egb_fim_block, total_fim_np, parameter_names + ) + elif obj_cons is not None and hasattr(obj_cons, "L"): + min_eig = np.min(np.real(np.linalg.eigvals(total_fim_np))) + if min_eig < _SMALL_TOLERANCE_DEFINITENESS: + jitter = np.min( + [ + -min_eig + _SMALL_TOLERANCE_DEFINITENESS, + _SMALL_TOLERANCE_DEFINITENESS, + ] + ) + else: + jitter = 0 + + fim_jittered = total_fim_np + jitter * np.eye(len(parameter_names)) + L_vals = np.linalg.cholesky(fim_jittered) + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + obj_cons.L[p, q].set_value(L_vals[i, j]) + + if hasattr(obj_cons, "L_inv"): + L_inv_vals = np.linalg.inv(L_vals) + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if i >= j: + obj_cons.L_inv[p, q].set_value(L_inv_vals[i, j]) + else: + obj_cons.L_inv[p, q].set_value(0.0) + + if hasattr(obj_cons, "fim_inv"): + fim_inv_np = np.linalg.inv(fim_jittered) + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + obj_cons.fim_inv[p, q].set_value(fim_inv_np[i, j]) + + if hasattr(obj_cons, "cov_trace"): + obj_cons.cov_trace.set_value(np.sum(L_inv_vals**2)) + + if obj_cons is not None and hasattr(obj_cons, "determinant"): + obj_cons.determinant.set_value(np.linalg.det(total_fim_np)) + + if obj_cons is not None and hasattr(obj_cons, "pseudo_trace"): + obj_cons.pseudo_trace.set_value(float(np.trace(total_fim_np))) + + def _build_optimize_experiments_results( + self, + res, + n_param_scenarios, + n_exp, + parameter_names, + template_mode, + resolved_init_method, + init_solver_name, + init_n_samples, + init_seed, + best_initial_points, + lhs_init_diagnostics, + lhs_initialization_time, + symmetry_breaking_info, + diagnostics_warnings, + final_solver_name, + build_time, + initialization_time, + solve_time, + ): + """Create the optimize_experiments() structured results payload.""" + solver_status = str(res.solver.status) + termination_condition = str(res.solver.termination_condition) + if isinstance(res.solver.message, str): + results_message = res.solver.message + elif isinstance(res.solver.message, bytes): + results_message = res.solver.message.decode("utf-8") + else: + results_message = ( + str(res.solver.message) if res.solver.message is not None else "" + ) + + def _safe_metric(metric_name, compute_fn, scenario_index): + try: + val = float(compute_fn()) + return val if np.isfinite(val) else float("nan") + except Exception as exc: + self.logger.warning( + f"Scenario {scenario_index}: failed to compute {metric_name}: {exc}. " + "Setting metric to NaN." + ) + return float("nan") + + param_scenarios = [] + for s in range(n_param_scenarios): + scenario = self.model.param_scenario_blocks[s] + total_fim_vals = [ + pyo.value(scenario.total_fim[p, q]) + for p in parameter_names + for q in parameter_names + ] + total_fim_np = np.array(total_fim_vals).reshape( + (len(parameter_names), len(parameter_names)) + ) + if self.only_compute_fim_lower: + total_fim_np = self._symmetrize_lower_tri(total_fim_np) + + log10_a_opt = _safe_metric( + "log10 A-opt", + lambda fim=total_fim_np: np.log10(np.trace(np.linalg.inv(fim))), + s, + ) + log10_pseudo_a_opt = _safe_metric( + "log10 pseudo A-opt", + lambda fim=total_fim_np: np.log10(np.trace(fim)), + s, + ) + log10_d_opt = _safe_metric( + "log10 D-opt", lambda fim=total_fim_np: np.log10(np.linalg.det(fim)), s + ) + log10_e_opt = _safe_metric( + "log10 E-opt", + lambda fim=total_fim_np: np.log10(min(np.linalg.eigvalsh(fim))), + s, + ) + log10_me_opt = _safe_metric( + "log10 ME-opt", + lambda fim=total_fim_np: np.log10(np.linalg.cond(fim)), + s, + ) + + scenario_structured = { + "param_scenario_id": s, + "param_scenario_weight": float(self.scenario_weights[s]), + "total_fim": total_fim_np.tolist(), + "parameter_values": self.get_unknown_parameter_values( + model=scenario.exp_blocks[0] + ), + "quality_metrics": { + "log10_a_opt": log10_a_opt, + "log10_pseudo_a_opt": log10_pseudo_a_opt, + "log10_d_opt": log10_d_opt, + "log10_e_opt": log10_e_opt, + "log10_me_opt": log10_me_opt, + }, + "experiments": [], + } + for k in range(n_exp): + exp_block = scenario.exp_blocks[k] + experiment_structured = { + "exp_id": k, + "design": self.get_experiment_input_values(model=exp_block), + "outputs": self.get_experiment_output_values(model=exp_block), + "fim": self.get_FIM(model=exp_block), + "sensitivity": None, + } + if hasattr(exp_block, "sensitivity_jacobian"): + experiment_structured["sensitivity"] = self.get_sensitivity_matrix( + model=exp_block + ) + scenario_structured["experiments"].append(experiment_structured) + + param_scenarios.append(scenario_structured) + + first_exp_block_fd = ( + self.model.param_scenario_blocks[0].exp_blocks[0].fd_scenario_blocks[0] + ) + design_variable_labels = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.experiment_inputs + ] + output_labels = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.experiment_outputs + ] + parameter_labels = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.unknown_parameters + ] + measurement_error_labels = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.measurement_error + ] + measurement_error_values = self.get_measurement_error_values( + model=self.model.param_scenario_blocks[0].exp_blocks[0] + ) + total_time = ( + build_time + lhs_initialization_time + initialization_time + solve_time + ) + initialization_method = ( + resolved_init_method.value if resolved_init_method is not None else "none" + ) + lhs_details = lhs_init_diagnostics or {} + + return { + "problem": { + "number_of_param_scenarios": n_param_scenarios, + "number_of_experiments_per_scenario": n_exp, + "used_template_experiment": template_mode, + "finite_difference_scheme": self._enum_label(self.fd_formula), + "finite_difference_step": self.step, + "scaled_nominal_parameters": self.scale_nominal_param_value, + "prior_fim": self.prior_FIM.tolist(), + "measurement_error_values": measurement_error_values, + "design_variables": design_variable_labels, + "outputs": output_labels, + "parameters": parameter_labels, + "measurement_errors": measurement_error_labels, + }, + "initialization": { + "method": initialization_method, + "solver": init_solver_name, + "samples_per_design_variable": ( + init_n_samples + if resolved_init_method + == InitializationMethod.latin_hypercube_sampling + else None + ), + "random_seed": ( + init_seed + if resolved_init_method + == InitializationMethod.latin_hypercube_sampling + else None + ), + "candidate_fim_evaluation_mode": lhs_details.get("candidate_fim_mode"), + "combination_scoring_mode": lhs_details.get("combo_mode"), + "number_of_candidate_designs": lhs_details.get("n_candidates"), + "number_of_design_combinations": lhs_details.get("n_combinations"), + "sampling_time_s": lhs_details.get("elapsed_sampling_s"), + "candidate_fim_evaluation_time_s": lhs_details.get( + "elapsed_fim_eval_s" + ), + "combination_scoring_time_s": lhs_details.get( + "elapsed_combo_scoring_s" + ), + "time_s": lhs_details.get("elapsed_total_s"), + "timed_out": lhs_details.get("timed_out"), + "selected_initial_designs": ( + best_initial_points + if resolved_init_method + == InitializationMethod.latin_hypercube_sampling + else None + ), + "best_initial_objective_value": lhs_details.get("best_obj"), + "best_initial_objective_value_log10": lhs_details.get("best_obj_log10"), + }, + "experiment_ordering": { + "design_variable": symmetry_breaking_info["variable"], + "chosen_by": symmetry_breaking_info["source"], + }, + "optimization_solve": { + "solver": final_solver_name, + "status": solver_status, + "termination_condition": termination_condition, + "message": results_message, + }, + "timing": { + "build_time_s": build_time, + "initialization_time_s": initialization_time, + "lhs_initialization_time_s": lhs_initialization_time, + "optimization_solve_time_s": solve_time, + "total_time_s": total_time, + }, + "warnings": diagnostics_warnings, + "solution": { + "objective": self._enum_label(self.objective_option), + "param_scenarios": param_scenarios, + }, + } + + # LHS-initialization helpers + def _get_experiment_input_vars(self, exp_block): + """ + Return the experiment-input Pyomo variable objects for an experiment + block, abstracting over the specific sensitivity-computation structure + (FD, AD, etc.). + + When the block exposes ``experiment_inputs`` directly (e.g. in a future + automatic-differentiation path), those are used. Otherwise the method + falls back to the FD structure (``exp_block.fd_scenario_blocks[0]``). + + Parameters + ---------- + exp_block : Pyomo Block + An ``exp_blocks[k]`` sub-block of the multi-experiment model. + + Returns + ------- + list + Ordered list of Pyomo :class:`Var` objects corresponding to the + experiment inputs. + """ + if hasattr(exp_block, "experiment_inputs"): + return list(exp_block.experiment_inputs.keys()) + # FD structure: inputs live on the base finite-difference scenario block + return list(exp_block.fd_scenario_blocks[0].experiment_inputs.keys()) + + def _set_experiment_inputs_fixed(self, n_param_scenarios, n_exp, fixed): + """Fix or unfix experiment inputs across all scenarios and experiments.""" + for s in range(n_param_scenarios): + for k in range(n_exp): + for comp in ( + self.model.param_scenario_blocks[s] + .exp_blocks[k] + .fd_scenario_blocks[0] + .experiment_inputs + ): + if fixed: + comp.fix() + else: + comp.unfix() + + @staticmethod + def _evaluate_objective_for_option(fim_matrix, objective_option): + _bad = ( + -np.inf + if objective_option in DesignOfExperiments._MAXIMIZE_OBJECTIVES + else np.inf + ) + + try: + if objective_option == ObjectiveLib.determinant: + return float(np.linalg.det(fim_matrix)) + elif objective_option == ObjectiveLib.pseudo_trace: + return float(np.trace(fim_matrix)) + elif objective_option == ObjectiveLib.trace: + return float(np.trace(np.linalg.inv(fim_matrix))) + elif objective_option == ObjectiveLib.minimum_eigenvalue: + return float(np.min(np.linalg.eigvalsh(fim_matrix))) + elif objective_option == ObjectiveLib.condition_number: + return float(np.linalg.cond(fim_matrix)) + else: # zero or unknown + return 0.0 + except (np.linalg.LinAlgError, ValueError): + return _bad + + def _evaluate_objective_from_fim(self, fim_matrix): + """ + Compute the scalar DoE objective from a numpy FIM matrix. + + Parameters + ---------- + fim_matrix : np.ndarray + Square FIM to score. + + Returns + ------- + float + Objective value. For maximisation objectives (``determinant``, + ``pseudo_trace``, ``minimum_eigenvalue``) a larger value is better. + For minimisation objectives (``trace`` / A-optimality and + ``condition_number`` / ME-optimality) a smaller value is better. + """ + return self._evaluate_objective_for_option(fim_matrix, self.objective_option) + + @staticmethod + def _symmetrize_lower_tri(mat): + """Mirror lower-triangle FIM entries to the upper triangle.""" + return mat + mat.T - np.diag(np.diag(mat)) + + @staticmethod + def _make_cholesky_rule(fim_expr, L_expr, parameter_names): + """ + Return a constraint rule that enforces ``fim_expr = L_expr * L_expr^T`` + on the lower-triangular portion defined by ``parameter_names``. + + The produced rule follows the Pyomo signature ``rule(block, p, q)`` + but does **not** actually use `block` in its body; the two matrix + expressions are captured from the enclosing scope. + + Parameters + ---------- + fim_expr : Var-like + Indexed by ``(p, q)``; usually ``model.fim`` or + ``scenario.total_fim``. + L_expr : Var-like + Indexed by ``(p, q)``; the corresponding lower-triangular + Cholesky factors. + parameter_names : Set + Pyomo Set listing the parameter indices in order. + """ + + def rule(_b, p, q): + p_idx = list(parameter_names).index(p) + q_idx = list(parameter_names).index(q) + if p_idx >= q_idx: + return fim_expr[p, q] == sum( + L_expr[p, parameter_names.at(k + 1)] + * L_expr[q, parameter_names.at(k + 1)] + for k in range(q_idx + 1) + ) + else: + return pyo.Constraint.Skip + + return rule + + @staticmethod + def _make_cholesky_inv_rule(fim_inv_expr, L_inv_expr, parameter_names): + """ + Return a rule that enforces ``fim_inv_expr = L_inv_expr^T * L_inv_expr`` + over the lower-triangular index region. + """ + + def rule(_b, p, q): + p_idx = list(parameter_names).index(p) + q_idx = list(parameter_names).index(q) + if p_idx >= q_idx: + return fim_inv_expr[p, q] == sum( + L_inv_expr[parameter_names.at(k + 1), p] + * L_inv_expr[parameter_names.at(k + 1), q] + for k in range(p_idx, len(parameter_names)) + ) + return pyo.Constraint.Skip + + return rule + + @staticmethod + def _make_cholesky_LLinv_rule(L_expr, L_inv_expr, parameter_names): + """ + Return a rule that enforces ``L_expr * L_inv_expr = I`` over the + lower-triangular index region. + """ + + def rule(_b, p, q): + p_idx = list(parameter_names).index(p) + q_idx = list(parameter_names).index(q) + if p_idx < q_idx: + return pyo.Constraint.Skip + target = 1 if p_idx == q_idx else 0 + return ( + sum( + L_expr[p, parameter_names.at(k + 1)] + * L_inv_expr[parameter_names.at(k + 1), q] + for k in range(len(parameter_names)) + ) + == target + ) + + return rule + + def _compute_fim_at_point_no_prior(self, experiment_index, input_values): + """ + Compute the FIM (without the prior FIM contribution) for the given + experiment at the specified experiment-input values using the + sequential finite-difference method. + + Parameters + ---------- + experiment_index : int + Index of the experiment in ``self.experiment_list`` to evaluate. + input_values : list + Numeric values for each experiment input variable (same order as + ``model.experiment_inputs``). + + Returns + ------- + np.ndarray + ``(n_params, n_params)`` FIM matrix, **excluding** the prior. + A zero matrix is returned on solver failure (with a warning). + """ + # Get a fresh labeled model for this experiment + model = ( + self.experiment_list[experiment_index] + .get_labeled_model(**self.get_labeled_model_args) + .clone() + ) + self.check_model_labels(model=model) + n_params = len(model.unknown_parameters) + + # Override experiment input values + update_model_from_suffix( + suffix_obj=model.experiment_inputs, values=input_values + ) + + # Temporarily zero the prior so that seq_FIM = Q^T * 𝚺^{-1} * Q only + saved_prior = self.prior_FIM + self.prior_FIM = np.zeros((n_params, n_params)) + + try: + self._sequential_FIM(model=model) + fim = self.seq_FIM.copy() + except Exception as exc: + self.logger.warning( + f"FIM evaluation failed at point {input_values}: {exc}. " + "Using zero matrix as fallback." + ) + fim = np.zeros((n_params, n_params)) + finally: + self.prior_FIM = saved_prior + + return fim + + def _lhs_initialize_experiments(self, lhs_n_samples, lhs_seed, n_exp): + """ + Use per-dimension Latin Hypercube Sampling to identify a good initial + experiment design for ``optimize_experiments``. + """ + sampling_timer = TicTocTimer() + + # 1. Get experiment-input bounds from the already-built model + first_exp_block = self.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = self._get_experiment_input_vars(first_exp_block) + n_inputs = len(exp_input_vars) + + missing = [v.name for v in exp_input_vars if v.lb is None or v.ub is None] + if missing: + raise ValueError( + "LHS initialization requires explicit lower and upper bounds on " + "all experiment input variables. The following variables are " + f"missing bounds: {missing}. " + "Set bounds in your experiment input variables before " + "calling ``optimize_experiments`` with " + "``init_method='latin_hypercube_sampling'``." + ) + + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + + # 2. Generate per-dimension 1-D LHS samples and take Cartesian product + # Split the user seed into per-dimension seeds so each 1-D LHS draw + # is independent while remaining reproducible for a fixed lhs_seed. + rng = np.random.default_rng(lhs_seed) + per_dim_samples = [] + for i in range(n_inputs): + dim_seed = int(rng.integers(0, 2**31)) + sampler = scipy.stats.qmc.LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + + candidate_points = tuple(product(*per_dim_samples)) + sampling_time = sampling_timer.toc(None) + n_candidates = len(candidate_points) + + if n_candidates > 10_000: + warnings.warn( + f"LHS initialization generated {n_candidates:,} candidate " + f"experiment designs (lhs_n_samples={lhs_n_samples}, " + f"n_inputs={n_inputs}). Evaluating FIM at all candidates may " + "take a long time. Consider reducing ``lhs_n_samples``.", + UserWarning, + stacklevel=2, + ) + + if hasattr(first_exp_block, "fd_scenario_blocks"): + n_scenarios_per_candidate = len(list(first_exp_block.fd_scenario_blocks)) + else: + n_scenarios_per_candidate = 1 + # Change the following block if we add support for LHS initialization with + # non-FD structures (e.g. AD) + if ( + not hasattr(first_exp_block, "fd_scenario_blocks") + or len(first_exp_block.fd_scenario_blocks) == 0 + ): + raise RuntimeError( + "_lhs_initialize_experiments requires finite-difference scenario " + "blocks on the experiment model. Ensure optimize_experiments is " + "using the sequential FIM path." + ) + + self.logger.info( + f"LHS: evaluating FIM at {n_candidates} candidate designs " + f"({n_candidates * n_scenarios_per_candidate} solver calls expected; " + "serial candidate FIM mode)." + ) + + n_params = len(first_exp_block.fd_scenario_blocks[0].unknown_parameters) + + # 3. Evaluate FIM at every candidate design + candidate_fims = [None] * n_candidates + fim_timer = TicTocTimer() + for pt_idx, pt in enumerate(candidate_points): + fim = self._compute_fim_at_point_no_prior( + experiment_index=0, input_values=list(pt) + ) + candidate_fims[pt_idx] = fim + if (pt_idx + 1) % max(1, n_candidates // 10) == 0: + self.logger.info(f" LHS FIM eval: {pt_idx + 1}/{n_candidates}") + + computed_pairs = [ + (pt, fim) + for pt, fim in zip(candidate_points, candidate_fims) + if fim is not None + ] + + # If no candidates were successfully evaluated, use the first candidate with + # a zero FIM to allow downstream code to run without missing data. + # This is an extreme fallback that should only occur if every candidate + # evaluation failed or if the time budget was too small to evaluate + # any candidates; in either case we want to avoid crashing and allow + # the user to get something back (with warnings) rather than nothing. + timed_out = False + if not computed_pairs: + computed_pairs = [(candidate_points[0], np.zeros((n_params, n_params)))] + + if len(computed_pairs) < n_candidates: + self.logger.warning( + "LHS candidate FIM evaluation scored " + f"{len(computed_pairs)}/{n_candidates} candidates; continuing " + "with best available subset." + ) + if len(computed_pairs) < n_exp: + first_pt, first_fim = computed_pairs[0] + computed_pairs.extend( + (first_pt, first_fim.copy()) + for _ in range(n_exp - len(computed_pairs)) + ) + _pts, _fims = zip(*computed_pairs) + candidate_points = tuple(_pts) + candidate_fims = tuple(_fims) + n_candidates = len(candidate_points) - L_vals = np.linalg.cholesky(fim_pd) - for i, c in enumerate(model.parameter_names): - for j, d in enumerate(model.parameter_names): - if i >= j: - model.L[c, d].value = L_vals[i, j] - else: - model.L[c, d].value = 0.0 + fim_eval_time = fim_timer.toc(None) + self.logger.info(f"LHS: completed FIM evaluations in {fim_eval_time:.1f}s.") - if hasattr(model, "L_inv"): - L_inv_vals = np.linalg.inv(L_vals) - for i, c in enumerate(model.parameter_names): - for j, d in enumerate(model.parameter_names): - if i >= j: - model.L_inv[c, d].value = L_inv_vals[i, j] - else: - model.L_inv[c, d].value = 0.0 + # 4. Enumerate combinations and score + n_combinations = math.comb(n_candidates, n_exp) + self.logger.info( + f"LHS: scoring {n_combinations:,} combinations of {n_exp} experiments..." + ) + if n_combinations > 100_000: + self.logger.warning( + f"LHS: {n_combinations:,} combinations to evaluate. " + "This may be slow. Consider reducing ``lhs_n_samples``." + ) - if hasattr(model, "fim_inv"): - # Use the pseudo-inverse here rather than the strict inverse. - # The jittered matrix should be positive definite, but ``pinv`` - # is safer for borderline ill-conditioned cases and matches the - # defensive approach already used when initializing ``fim_inv`` - # from user-provided starting values. - fim_inv_vals = np.linalg.pinv(fim_pd) - for i, c in enumerate(model.parameter_names): - for j, d in enumerate(model.parameter_names): - if self.only_compute_fim_lower and i < j: - model.fim_inv[c, d].value = 0.0 - else: - model.fim_inv[c, d].value = fim_inv_vals[i, j] + combo_timer = TicTocTimer() + prior = self.prior_FIM.copy() + _obj_option = self.objective_option + is_maximize = _obj_option in self._MAXIMIZE_OBJECTIVES + best_obj = -np.inf if is_maximize else np.inf + best_combo = None + _score_obj = DesignOfExperiments._evaluate_objective_for_option + + for combo in combinations(range(n_candidates), n_exp): + if n_exp == 2: + i, j = combo + fim_total = prior + candidate_fims[i] + candidate_fims[j] + else: + fim_total = prior.copy() + for idx in combo: + fim_total = fim_total + candidate_fims[idx] + obj_val = _score_obj(fim_total, _obj_option) + if is_maximize: + if obj_val > best_obj: + best_obj = obj_val + best_combo = combo + else: + if obj_val < best_obj: + best_obj = obj_val + best_combo = combo - if hasattr(model, "cov_trace"): - model.cov_trace.value = np.trace(fim_inv_vals) + combo_time = combo_timer.toc(None) + + if best_combo is None: + self.logger.warning( + "LHS combination scoring ended before any combination was scored. " + "Falling back to the first n_exp candidate points." + ) + best_combo = tuple(range(n_exp)) + if n_exp == 2: + i, j = best_combo + fim_total = prior + candidate_fims[i] + candidate_fims[j] + else: + fim_total = prior.copy() + for idx in best_combo: + fim_total = fim_total + candidate_fims[idx] + best_obj = float(_score_obj(fim_total, _obj_option)) + + best_obj_log10 = ( + float(np.log10(best_obj)) + if np.isfinite(best_obj) and best_obj > 0 + else None + ) + self.logger.info( + f"LHS: best {self.objective_option.value} objective = {best_obj:.6g} " + f"(combo scoring took {combo_time:.1f}s)." + ) + + best_initial_points = [list(candidate_points[i]) for i in best_combo] + self.logger.info( + f"LHS initial design: " + + ", ".join(f"exp[{k}]={best_initial_points[k]}" for k in range(n_exp)) + ) + + lhs_diagnostics = { + "candidate_fim_mode": "serial", + "combo_mode": "serial", + "n_candidates": n_candidates, + "n_combinations": n_combinations, + "elapsed_sampling_s": sampling_time, + "elapsed_fim_eval_s": fim_eval_time, + "elapsed_combo_scoring_s": combo_time, + "elapsed_total_s": sampling_time + fim_eval_time + combo_time, + "timed_out": timed_out, + "best_obj": best_obj, + "best_obj_log10": best_obj_log10, + } + return best_initial_points, lhs_diagnostics # Perform multi-experiment doe (sequential, or ``greedy`` approach) def run_multi_doe_sequential(self, N_exp=1): @@ -594,14 +1974,28 @@ def compute_FIM(self, model=None, method="sequential"): method: string to specify which method should be used options are ``kaug`` and ``sequential`` + Notes + ----- + When ``model is None`` and ``experiment_list`` contains multiple + experiments, this method computes each experiment FIM and returns + the aggregate: + + ``total_fim = sum(fim_i for each experiment i) + prior_FIM`` + + where each ``fim_i`` excludes the prior contribution. + Returns ------- computed FIM: 2D numpy array of the FIM """ + aggregate_all_experiments = model is None and len(self.experiment_list) > 1 + if model is None: - self.compute_FIM_model = self.experiment.get_labeled_model( - **self.get_labeled_model_args - ).clone() + self.compute_FIM_model = ( + self.experiment_list[0] + .get_labeled_model(**self.get_labeled_model_args) + .clone() + ) model = self.compute_FIM_model else: # TODO: Add safe naming when a model is passed by the user. @@ -631,20 +2025,90 @@ def compute_FIM(self, model=None, method="sequential"): # TODO: Add a check to see if the model has an objective and deactivate it. # This solve should only be a square solve without any obj function. - if method == "sequential": - self._sequential_FIM(model=model) - self._computed_FIM = self.seq_FIM - elif method == "kaug": - self._kaug_FIM(model=model) - self._computed_FIM = self.kaug_FIM - else: - raise ValueError( - ( - "The method provided, {}, must be either `sequential` " - "or `kaug`".format(method) + def _compute_fim_for_model(eval_model): + if method == "sequential": + self._sequential_FIM(model=eval_model) + return np.array(self.seq_FIM, copy=True) + elif method == "kaug": + self._kaug_FIM(model=eval_model) + return np.array(self.kaug_FIM, copy=True) + else: + raise ValueError( + ( + "The method provided, {}, must be either `sequential` " + "or `kaug`".format(method) + ) ) + + def _unknown_parameter_signature(eval_model): + # Use stable model-local component identifiers and values so we can + # verify all experiments are consistent before aggregating FIMs. + names = [ + str(pyo.ComponentUID(param, context=eval_model)) + for param in eval_model.unknown_parameters + ] + values = np.array( + [ + float(pyo.value(val)) + for val in eval_model.unknown_parameters.values() + ] + ) + return names, values + + if aggregate_all_experiments: + saved_prior = self.prior_FIM + self._computed_FIM_by_experiment = [] + # Capture the baseline parameter labels/values from experiment 0. + reference_param_names, reference_param_values = ( + _unknown_parameter_signature(model) ) + try: + # Compute each experiment FIM without prior so the aggregate adds + # prior exactly once at the end. + self.prior_FIM = np.zeros(saved_prior.shape) + for idx, exp in enumerate(self.experiment_list): + if idx == 0: + # Reuse the already-built model for experiment 0. + exp_model = model + else: + exp_model = exp.get_labeled_model( + **self.get_labeled_model_args + ).clone() + self.check_model_labels(model=exp_model) + + param_names, param_values = _unknown_parameter_signature(exp_model) + # Reject heterogeneous parameter spaces before solving so + # summed FIMs remain well-defined. + if param_names != reference_param_names: + raise ValueError( + "All experiments in 'experiment_list' must share the same " + "unknown parameter labels and order for compute_FIM " + f"aggregation. Mismatch detected at experiment index {idx}." + ) + if not np.allclose( + param_values, reference_param_values, atol=1e-12, rtol=1e-12 + ): + raise ValueError( + "All experiments in 'experiment_list' must share the same " + "unknown parameter values for compute_FIM aggregation. " + f"Mismatch detected at experiment index {idx}." + ) + + fim_i = _compute_fim_for_model(exp_model) + self._computed_FIM_by_experiment.append(fim_i) + finally: + self.prior_FIM = saved_prior + + # Aggregate all experiment FIMs and add the saved prior once. + total_fim = np.zeros(saved_prior.shape) + for fim_i in self._computed_FIM_by_experiment: + total_fim = total_fim + fim_i + self._computed_FIM = total_fim + saved_prior + else: + self._computed_FIM = _compute_fim_for_model(model) + self._computed_FIM_by_experiment = [np.array(self._computed_FIM, copy=True)] + return self._computed_FIM # Use a sequential method to get the FIM @@ -658,9 +2122,11 @@ def _sequential_FIM(self, model=None): """ # Build a single model instance if model is None: - self.compute_FIM_model = self.experiment.get_labeled_model( - **self.get_labeled_model_args - ).clone() + self.compute_FIM_model = ( + self.experiment_list[0] + .get_labeled_model(**self.get_labeled_model_args) + .clone() + ) model = self.compute_FIM_model # Create suffix to keep track of parameter scenarios @@ -819,9 +2285,11 @@ def _kaug_FIM(self, model=None): # Remake compute_FIM_model if model is None. # compute_FIM_model needs to be the right version for function to work. if model is None: - self.compute_FIM_model = self.experiment.get_labeled_model( - **self.get_labeled_model_args - ).clone() + self.compute_FIM_model = ( + self.experiment_list[0] + .get_labeled_model(**self.get_labeled_model_args) + .clone() + ) model = self.compute_FIM_model # add zero (dummy/placeholder) objective function @@ -900,7 +2368,9 @@ def _kaug_FIM(self, model=None): self.kaug_FIM = self.kaug_jac.T @ cov_y @ self.kaug_jac + self.prior_FIM # Create the DoE model (with ``scenarios`` from finite differencing scheme) - def create_doe_model(self, model=None): + def create_doe_model( + self, model=None, experiment_index=0, _for_multi_experiment=False + ): """ Add equations to compute sensitivities, FIM, and objective. Builds the DoE model. Adds the scenarios, the sensitivity matrix @@ -914,6 +2384,9 @@ def create_doe_model(self, model=None): Parameters ---------- model: model to add finite difference scenarios + experiment_index: index of experiment in experiment_list to use for this model (default: 0) + _for_multi_experiment: if True, skip creating L matrix and other objective-related + variables that will be created at the aggregated level (default: False) """ if model is None: @@ -942,25 +2415,27 @@ def create_doe_model(self, model=None): ) # Generate scenarios for finite difference formulae - self._generate_scenario_blocks(model=model) + self._generate_fd_scenario_blocks( + model=model, experiment_index=experiment_index + ) # Set names for indexing sensitivity matrix (jacobian) and FIM scen_block_ind = min( [ - k.name.split(".").index("scenario_blocks[0]") - for k in model.scenario_blocks[0].unknown_parameters.keys() + k.name.split(".").index("fd_scenario_blocks[0]") + for k in model.fd_scenario_blocks[0].unknown_parameters.keys() ] ) model.parameter_names = pyo.Set( initialize=[ ".".join(k.name.split(".")[(scen_block_ind + 1) :]) - for k in model.scenario_blocks[0].unknown_parameters.keys() + for k in model.fd_scenario_blocks[0].unknown_parameters.keys() ] ) model.output_names = pyo.Set( initialize=[ ".".join(k.name.split(".")[(scen_block_ind + 1) :]) - for k in model.scenario_blocks[0].experiment_outputs.keys() + for k in model.fd_scenario_blocks[0].experiment_outputs.keys() ] ) @@ -1042,9 +2517,10 @@ def initialize_fim_inv(m, j, d): # To-Do: Look into this functionality..... # if cholesky, define L elements as variables - if self.Cholesky_option and self.objective_option in ( - ObjectiveLib.determinant, - ObjectiveLib.trace, + if ( + not _for_multi_experiment + and self.Cholesky_option + and self.objective_option in (ObjectiveLib.determinant, ObjectiveLib.trace) ): model.L = pyo.Var( model.parameter_names, model.parameter_names, initialize=identity_matrix @@ -1096,12 +2572,14 @@ def jacobian_rule(m, n, p): s1 = 0 s2 = param_ind + 1 - var_up = cuid.find_component_on(m.scenario_blocks[s1]) - var_lo = cuid.find_component_on(m.scenario_blocks[s2]) + var_up = cuid.find_component_on(m.fd_scenario_blocks[s1]) + var_lo = cuid.find_component_on(m.fd_scenario_blocks[s2]) param = m.parameter_scenarios[max(s1, s2)] - param_loc = pyo.ComponentUID(param).find_component_on(m.scenario_blocks[0]) - param_val = m.scenario_blocks[0].unknown_parameters[param_loc] + param_loc = pyo.ComponentUID(param).find_component_on( + m.fd_scenario_blocks[0] + ) + param_val = m.fd_scenario_blocks[0].unknown_parameters[param_loc] param_diff = param_val * fd_step_mult * self.step if self.scale_nominal_param_value: @@ -1154,20 +2632,38 @@ def fim_rule(m, p, q): else: return m.fim[p, q] == m.fim[q, p] else: - return ( - m.fim[p, q] - == sum( + # For multi-experiment optimization, prior_FIM is added to the + # aggregated total_fim, not to each individual experiment's FIM + if _for_multi_experiment: + return m.fim[p, q] == sum( 1 - / m.scenario_blocks[0].measurement_error[ - pyo.ComponentUID(n).find_component_on(m.scenario_blocks[0]) + / m.fd_scenario_blocks[0].measurement_error[ + pyo.ComponentUID(n).find_component_on( + m.fd_scenario_blocks[0] + ) ] ** 2 * m.sensitivity_jacobian[n, p] * m.sensitivity_jacobian[n, q] for n in m.output_names ) - + m.prior_FIM[p, q] - ) + else: + return ( + m.fim[p, q] + == sum( + 1 + / m.fd_scenario_blocks[0].measurement_error[ + pyo.ComponentUID(n).find_component_on( + m.fd_scenario_blocks[0] + ) + ] + ** 2 + * m.sensitivity_jacobian[n, p] + * m.sensitivity_jacobian[n, q] + for n in m.output_names + ) + + m.prior_FIM[p, q] + ) model.jacobian_constraint = pyo.Constraint( model.output_names, model.parameter_names, rule=jacobian_rule @@ -1188,7 +2684,7 @@ def fim_rule(m, p, q): model.fim_inv[p, q].fix(0.0) # Create scenario block structure - def _generate_scenario_blocks(self, model=None): + def _generate_fd_scenario_blocks(self, model=None, experiment_index=0): """ Generates the modeling blocks corresponding to the scenarios for the finite differencing scheme to compute the sensitivity jacobian @@ -1202,15 +2698,18 @@ def _generate_scenario_blocks(self, model=None): Parameters ---------- model: model to add finite difference scenarios + experiment_index: index of experiment in experiment_list to use for this model (default: 0) """ # If model is none, assume it is self.model if model is None: model = self.model # Generate initial scenario to populate unknown parameter values - model.base_model = self.experiment.get_labeled_model( - **self.get_labeled_model_args - ).clone() + model.base_model = ( + self.experiment_list[experiment_index] + .get_labeled_model(**self.get_labeled_model_args) + .clone() + ) # Check the model that labels are correct self.check_model_labels(model=model.base_model) @@ -1296,8 +2795,9 @@ def _generate_scenario_blocks(self, model=None): # Generate blocks for finite difference scenarios def build_block_scenarios(b, s): # Generate model for the finite difference scenario - m = b.model() - b.transfer_attributes_from(m.base_model.clone()) + # Get the parent block that contains base_model + parent_block = b.parent_block() + b.transfer_attributes_from(parent_block.base_model.clone()) # Forward/Backward difference have a stationary # case (s == 0), no parameter to perturb @@ -1308,7 +2808,7 @@ def build_block_scenarios(b, s): if s == 0: return - param = m.parameter_scenarios[s] + param = parent_block.parameter_scenarios[s] # Perturbation to be (1 + diff) * param_value if self.fd_formula == FiniteDifferenceStep.central: @@ -1325,9 +2825,9 @@ def build_block_scenarios(b, s): pass # Update parameter values for the given finite difference scenario - pyo.ComponentUID(param, context=m.base_model).find_component_on( + pyo.ComponentUID(param, context=parent_block.base_model).find_component_on( b - ).set_value(m.base_model.unknown_parameters[param] * (1 + diff)) + ).set_value(parent_block.base_model.unknown_parameters[param] * (1 + diff)) # Fix experiment inputs before solve (enforce square solve) for comp in b.experiment_inputs: @@ -1339,11 +2839,15 @@ def build_block_scenarios(b, s): for comp in b.experiment_inputs: comp.unfix() - model.scenario_blocks = pyo.Block(model.scenarios, rule=build_block_scenarios) + model.fd_scenario_blocks = pyo.Block( + model.scenarios, rule=build_block_scenarios + ) # TODO: this might have to change if experiment inputs have # a different value in the Suffix (currently it is the CUID) - design_vars = [k for k, v in model.scenario_blocks[0].experiment_inputs.items()] + design_vars = [ + k for k, v in model.fd_scenario_blocks[0].experiment_inputs.items() + ] # Add constraints to equate block design with global design: for ind, d in enumerate(design_vars): @@ -1354,8 +2858,8 @@ def global_design_fixing(m, s): if s == 0: return pyo.Constraint.Skip block_design_var = pyo.ComponentUID( - d, context=m.scenario_blocks[0] - ).find_component_on(m.scenario_blocks[s]) + d, context=m.fd_scenario_blocks[0] + ).find_component_on(m.fd_scenario_blocks[s]) return d == block_design_var model.add_component( @@ -1386,6 +2890,15 @@ def create_objective_function(self, model=None): if model is None: model = self.model + if self.objective_option in [ + ObjectiveLib.minimum_eigenvalue, + ObjectiveLib.condition_number, + ]: + raise ValueError( + f"objective_option='{self.objective_option.value}' requires " + "use_grey_box_objective=True." + ) + if self.objective_option not in [ ObjectiveLib.determinant, ObjectiveLib.trace, @@ -1433,23 +2946,10 @@ def create_objective_function(self, model=None): for j, d in enumerate(model.parameter_names): model.L_inv[c, d].value = L_inv[i, j] - def cholesky_imp(b, c, d): - """ - Calculate Cholesky L matrix using algebraic constraints - """ - # If the row is greater than or equal to the column, we are in the - # lower triangle region of the L and FIM matrices. - # This region is where our equations are well-defined. - m = b.model() - if list(m.parameter_names).index(c) >= list(m.parameter_names).index(d): - return m.fim[c, d] == sum( - m.L[c, m.parameter_names.at(k + 1)] - * m.L[d, m.parameter_names.at(k + 1)] - for k in range(list(m.parameter_names).index(d) + 1) - ) - else: - # This is the empty half of L above the diagonal - return pyo.Constraint.Skip + if self.objective_option == ObjectiveLib.trace and not self.Cholesky_option: + raise ValueError( + "objective_option='trace' currently only implemented with ``_Cholesky option=True``." + ) # If trace objective, need L inverse constraints if self.Cholesky_option and self.objective_option == ObjectiveLib.trace: @@ -1544,6 +3044,7 @@ def determinant_general(b): list_p = list(object_p) # generate a name_order to iterate \sigma_i + # NOTE: Not used in calculation. Delete? det_perm = 0 for i in range(len(list_p)): name_order = [] @@ -1568,7 +3069,11 @@ def determinant_general(b): if self.Cholesky_option and self.objective_option == ObjectiveLib.determinant: model.obj_cons.cholesky_cons = pyo.Constraint( - model.parameter_names, model.parameter_names, rule=cholesky_imp + model.parameter_names, + model.parameter_names, + rule=self._make_cholesky_rule( + model.fim, model.L, model.parameter_names + ), ) model.objective = pyo.Objective( expr=2 * sum(pyo.log10(model.L[j, j]) for j in model.parameter_names), @@ -1587,17 +3092,17 @@ def determinant_general(b): ) elif self.objective_option == ObjectiveLib.trace: - if not self.Cholesky_option: - raise ValueError( - "objective_option='trace' currently only implemented with ``_Cholesky option=True``." - ) # if Cholesky and trace, calculating # the OBJ with trace model.cov_trace = pyo.Var( initialize=np.trace(np.linalg.inv(fim)), bounds=(small_number, None) ) model.obj_cons.cholesky_cons = pyo.Constraint( - model.parameter_names, model.parameter_names, rule=cholesky_imp + model.parameter_names, + model.parameter_names, + rule=self._make_cholesky_rule( + model.fim, model.L, model.parameter_names + ), ) model.obj_cons.cholesky_inv_cons = pyo.Constraint( model.parameter_names, model.parameter_names, rule=cholesky_inv_imp @@ -1634,11 +3139,386 @@ def determinant_general(b): # add dummy objective function model.objective = pyo.Objective(expr=0) - def create_grey_box_objective_function(self, model=None): + def create_multi_experiment_objective_function(self, model): + """ + Create objective for multi-experiment optimization. + + For each scenario s: + 1. Creates total_fim[s] = sum of exp_blocks[k].fim + prior_FIM + 2. Creates Cholesky/determinant/trace variables and constraints per scenario + 3. Creates single top-level objective summing across scenarios + + Parameters + ---------- + model: model with param_scenario_blocks structure + """ + # Validate objective option for multi-experiment. + if self.use_grey_box: + if self.objective_option == ObjectiveLib.zero: + raise ValueError( + "Grey-box objective support in optimize_experiments() is only " + "available for objective_option in ['determinant', 'trace', " + "'pseudo_trace', 'minimum_eigenvalue', 'condition_number']." + ) + self._grey_box_output_name() + else: + if self.objective_option in [ + ObjectiveLib.minimum_eigenvalue, + ObjectiveLib.condition_number, + ]: + raise ValueError( + f"objective_option='{self.objective_option.value}' requires " + "use_grey_box_objective=True." + ) + if self.objective_option not in [ + ObjectiveLib.determinant, + ObjectiveLib.trace, + ObjectiveLib.pseudo_trace, + ObjectiveLib.zero, + ]: + raise DeveloperError( + "Objective option not recognized for multi-experiment optimization. " + "Please contact the developers as you should not see this error." + ) + + # Validate trace objective requires Cholesky option + if ( + not self.use_grey_box + and self.objective_option == ObjectiveLib.trace + and not self.Cholesky_option + ): + raise ValueError( + "objective_option='trace' currently only implemented with " + "``_Cholesky_option=True`` or ``use_grey_box_objective=True``." + ) + + small_number = 1e-10 + n_scenarios = len(model.param_scenario_blocks) + + # Get weights from instance attribute (set in optimize_experiments) and + # default to uniform weights if not provided + # retrieve weights; default to uniform tuple of appropriate length + default_weights = tuple([1.0 / n_scenarios] * n_scenarios) + scenario_weights = getattr(self, 'scenario_weights', default_weights) + + # Get parameter names from first experiment (same across all) + parameter_names = model.param_scenario_blocks[0].exp_blocks[0].parameter_names + n_exp = len(model.param_scenario_blocks[0].exp_blocks) + + # For each scenario: create aggregated FIM and constraints + for s in range(n_scenarios): + scenario = model.param_scenario_blocks[s] + + # 1. Create aggregated FIM variable for each scenario: + # total_fim = sum of all exp FIMs + prior_FIM + scenario.total_fim = pyo.Var(parameter_names, parameter_names) + + # 2. Constraint: total_fim[p,q] = sum_k (exp_blocks[k].fim[p,q]) + # + prior_FIM[p,q] + def total_fim_rule(b, p, q): + p_idx = list(parameter_names).index(p) + q_idx = list(parameter_names).index(q) + + # When only_compute_fim_lower=True, only compute lower triangle + # Upper triangle elements will be handled through symmetry + if self.only_compute_fim_lower and p_idx < q_idx: + return pyo.Constraint.Skip + + return b.total_fim[p, q] == ( + sum(b.exp_blocks[k].fim[p, q] for k in range(n_exp)) + + self.prior_FIM[p_idx, q_idx] + ) + + scenario.total_fim_cons = pyo.Constraint( + parameter_names, parameter_names, rule=total_fim_rule + ) + + # 3. Fix upper triangle elements to 0 if only_compute_fim_lower=True + # Initialize lower triangle from sum of individual FIMs + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if self.only_compute_fim_lower and i < j: + # Fix upper triangle to 0 + scenario.total_fim[p, q].fix(0.0) + else: + # Initialize lower triangle and diagonal + fim_sum = sum( + pyo.value(scenario.exp_blocks[k].fim[p, q]) or 0 + for k in range(n_exp) + ) + scenario.total_fim[p, q].value = fim_sum + self.prior_FIM[i, j] + + # 4. Build the objective representation for this scenario. + if self.use_grey_box: + # In multi-experiment mode, the grey box should see the + # aggregated scenario FIM rather than any one experiment block. + # That keeps the objective mathematically aligned with the + # public API while reusing the same grey box implementation as + # the single-experiment path. + self.create_grey_box_objective_function( + model=scenario, + fim_expr=scenario.total_fim, + parameter_names=parameter_names, + build_objective=False, + ) + # 5. Create variables and constraints (initialization will happen after square solve) + elif ( + self.Cholesky_option + and self.objective_option == ObjectiveLib.determinant + ): + # (similar to single-experiment case in create_objective_function) + scenario.obj_cons = pyo.Block() + # Add lower triangular Cholesky variables per scenario + scenario.obj_cons.L = pyo.Var(parameter_names, parameter_names) + + # Fix upper triangle to 0 and set lower bound on diagonal + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + # Fix upper triangle to 0 + if i < j: + scenario.obj_cons.L[p, q].fix(0.0) + # Lower bound on diagonal + elif i == j and self.L_diagonal_lower_bound: + scenario.obj_cons.L[p, q].setlb(self.L_diagonal_lower_bound) + + # reuse shared helper to create the constraint rule + cholesky_rule = self._make_cholesky_rule( + scenario.total_fim, scenario.obj_cons.L, parameter_names + ) + scenario.obj_cons.cholesky_cons = pyo.Constraint( + parameter_names, parameter_names, rule=cholesky_rule + ) + + elif self.Cholesky_option and self.objective_option == ObjectiveLib.trace: + scenario.obj_cons = pyo.Block() + # Add lower triangular Cholesky variables per scenario + scenario.obj_cons.L = pyo.Var(parameter_names, parameter_names) + scenario.obj_cons.L_inv = pyo.Var(parameter_names, parameter_names) + scenario.obj_cons.fim_inv = pyo.Var(parameter_names, parameter_names) + scenario.obj_cons.cov_trace = pyo.Var(bounds=(small_number, None)) + + # Fix upper triangle of L and L_inv to 0 and set lower bound on diagonal + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + # Fix upper triangle to 0 + if i < j: + scenario.obj_cons.L[p, q].fix(0.0) + scenario.obj_cons.L_inv[p, q].fix(0.0) + # Lower bound on diagonal + elif i == j and self.L_diagonal_lower_bound: + scenario.obj_cons.L[p, q].setlb(self.L_diagonal_lower_bound) + + # reuse shared helper to create the constraint rule + cholesky_rule = self._make_cholesky_rule( + scenario.total_fim, scenario.obj_cons.L, parameter_names + ) + scenario.obj_cons.cholesky_cons = pyo.Constraint( + parameter_names, parameter_names, rule=cholesky_rule + ) + + # reuse helpers for the inverse and identity rules instead of + # re-implementing the logic in-place + cholesky_inv_rule = self._make_cholesky_inv_rule( + scenario.obj_cons.fim_inv, scenario.obj_cons.L_inv, parameter_names + ) + + cholesky_LLinv_rule = self._make_cholesky_LLinv_rule( + scenario.obj_cons.L, scenario.obj_cons.L_inv, parameter_names + ) + + # Covariance trace calculation + def cov_trace_rule(b): + return b.cov_trace == sum(b.fim_inv[j, j] for j in parameter_names) + + # Add all constraints + scenario.obj_cons.cholesky_inv_cons = pyo.Constraint( + parameter_names, parameter_names, rule=cholesky_inv_rule + ) + scenario.obj_cons.cholesky_LLinv_cons = pyo.Constraint( + parameter_names, parameter_names, rule=cholesky_LLinv_rule + ) + scenario.obj_cons.cov_trace_cons = pyo.Constraint(rule=cov_trace_rule) + + # Optional: improve Cholesky roundoff error + if self.improve_cholesky_roundoff_error: + + def cholesky_fim_diag(b, p, q): + return scenario.total_fim[p, p] >= b.L[p, q] ** 2 + + def cholesky_fim_inv_diag(b, p, q): + return b.fim_inv[p, p] >= b.L_inv[p, q] ** 2 + + scenario.obj_cons.cholesky_fim_diag_cons = pyo.Constraint( + parameter_names, parameter_names, rule=cholesky_fim_diag + ) + scenario.obj_cons.cholesky_fim_inv_diag_cons = pyo.Constraint( + parameter_names, parameter_names, rule=cholesky_fim_inv_diag + ) + + elif self.objective_option == ObjectiveLib.determinant: + scenario.obj_cons = pyo.Block() + # Non-Cholesky determinant: create determinant var per scenario + scenario.obj_cons.determinant = pyo.Var(bounds=(small_number, None)) + + # Determinant constraint (explicit formula) + def determinant_general(b): + r_list = list(range(len(parameter_names))) + object_p = permutations(r_list) + list_p = list(object_p) + + det_perm = sum( + self._sgn(list_p[d]) + * math.prod( + scenario.total_fim[ + parameter_names.at(val + 1), parameter_names.at(ind + 1) + ] + for ind, val in enumerate(list_p[d]) + ) + for d in range(len(list_p)) + ) + return b.determinant == det_perm + + scenario.obj_cons.determinant_cons = pyo.Constraint( + rule=determinant_general + ) + + elif self.objective_option == ObjectiveLib.pseudo_trace: + scenario.obj_cons = pyo.Block() + # Pseudo trace objective (Trace of FIM) + scenario.obj_cons.pseudo_trace = pyo.Var(bounds=(small_number, None)) + + # Pseudo trace constraint + def pseudo_trace_rule(b): + return b.pseudo_trace == sum( + scenario.total_fim[j, j] for j in parameter_names + ) + + scenario.obj_cons.pseudo_trace_cons = pyo.Constraint( + rule=pseudo_trace_rule + ) + + # 5. Create single top-level objective summing across scenarios + if self.use_grey_box: + model.objective = pyo.Objective( + expr=sum( + scenario_weights[s] + * model.param_scenario_blocks[s].obj_cons.egb_fim_block.outputs[ + self._grey_box_output_name() + ] + for s in range(n_scenarios) + ), + sense=( + pyo.maximize + if self.objective_option in self._MAXIMIZE_OBJECTIVES + else pyo.minimize + ), + ) + + elif self.Cholesky_option and self.objective_option == ObjectiveLib.determinant: + model.objective = pyo.Objective( + expr=sum( + scenario_weights[s] + * 2 + * sum( + pyo.log10(model.param_scenario_blocks[s].obj_cons.L[j, j]) + for j in parameter_names + ) + for s in range(n_scenarios) + ), + sense=pyo.maximize, + ) + + elif self.objective_option == ObjectiveLib.determinant: + model.objective = pyo.Objective( + expr=sum( + scenario_weights[s] + * pyo.log10( + model.param_scenario_blocks[s].obj_cons.determinant + + _SMALL_TOLERANCE_DEFINITENESS # to avoid log(0) + ) + for s in range(n_scenarios) + ), + sense=pyo.maximize, + ) + + elif self.Cholesky_option and self.objective_option == ObjectiveLib.trace: + model.objective = pyo.Objective( + expr=sum( + scenario_weights[s] + * model.param_scenario_blocks[s].obj_cons.cov_trace + for s in range(n_scenarios) + ), + sense=pyo.minimize, + ) + + elif self.objective_option == ObjectiveLib.pseudo_trace: + model.objective = pyo.Objective( + expr=sum( + scenario_weights[s] + * pyo.log10( + model.param_scenario_blocks[s].obj_cons.pseudo_trace + + _SMALL_TOLERANCE_DEFINITENESS # to avoid log(0) + ) + for s in range(n_scenarios) + ), + sense=pyo.maximize, + ) + + elif self.objective_option == ObjectiveLib.zero: + # Dummy objective + model.objective = pyo.Objective(expr=0) + + def create_grey_box_objective_function( + self, model=None, fim_expr=None, parameter_names=None, build_objective=True + ): + """ + Attach the FIM grey box block to a model or scenario block. + + The same helper is used for both public APIs. In the single-experiment + path it connects to ``model.fim`` and also creates the top-level + objective. In multi-experiment mode it connects to ``scenario.total_fim`` + so each scenario contributes its aggregated FIM metric to the outer + objective sum. + """ + from pyomo.contrib.doe.grey_box_utilities import FIMExternalGreyBox + from pyomo.contrib.pynumero.interfaces.external_grey_box import ( + ExternalGreyBoxBlock, + ) + # Add external grey box block to a block named ``obj_cons`` to # reuse material for initializing the objective-free square model if model is None: - model = model = self.model + model = self.model + if fim_expr is None: + if not hasattr(model, "fim"): + raise RuntimeError( + "Model provided does not have variable `fim`. Please make " + "sure the model is built properly before creating the grey " + "box objective." + ) + fim_expr = model.fim + if parameter_names is None: + if not hasattr(model, "parameter_names"): + raise RuntimeError( + "Model provided does not define `parameter_names`. Please " + "make sure the model is built properly before creating the " + "grey box objective." + ) + parameter_names = model.parameter_names + + output_name = self._grey_box_output_name() + + # The external model expects a dense symmetric FIM as its initializer. + # When the algebraic model stores only the lower triangle, mirror those + # values here so the grey box starts from the same physical matrix. + fim_vals = [ + pyo.value(fim_expr[p, q]) for p in parameter_names for q in parameter_names + ] + fim_initial = np.array(fim_vals, dtype=np.float64).reshape( + (len(parameter_names), len(parameter_names)) + ) + if self.only_compute_fim_lower: + fim_initial = self._symmetrize_lower_tri(fim_initial) # TODO: Make this naming convention robust model.obj_cons = pyo.Block() @@ -1646,6 +3526,8 @@ def create_grey_box_objective_function(self, model=None): # Create FIM External Grey Box object grey_box_FIM = FIMExternalGreyBox( doe_object=self, + parameter_names=parameter_names, + fim_initial=fim_initial, objective_option=self.objective_option, logger_level=self.logger.getEffectiveLevel(), ) @@ -1664,48 +3546,43 @@ def FIM_egb_cons(m, p1, p2): p2: parameter 2 """ - # Using upper triangular FIM to - # make numerics better. - if list(model.parameter_names).index(p1) >= list( - model.parameter_names - ).index(p2): - return model.fim[(p1, p2)] == m.egb_fim_block.inputs[(p2, p1)] + # The grey box receives only the triangular portion of the FIM. We + # therefore map the lower-triangular algebraic entries into the + # upper-triangular input names expected by ExternalGreyBoxBlock. + if list(parameter_names).index(p1) >= list(parameter_names).index(p2): + return fim_expr[(p1, p2)] == m.egb_fim_block.inputs[(p2, p1)] else: return pyo.Constraint.Skip # Add the FIM and External Grey # Box inputs constraints model.obj_cons.FIM_equalities = pyo.Constraint( - model.parameter_names, model.parameter_names, rule=FIM_egb_cons + parameter_names, parameter_names, rule=FIM_egb_cons ) - # Add objective based on user provided - # type within ObjectiveLib - if self.objective_option == ObjectiveLib.trace: - model.objective = pyo.Objective( - expr=model.obj_cons.egb_fim_block.outputs["A-opt"], sense=pyo.minimize - ) - elif self.objective_option == ObjectiveLib.pseudo_trace: - model.objective = pyo.Objective( - expr=model.obj_cons.egb_fim_block.outputs["pseudo-A-opt"], - sense=pyo.maximize, - ) - elif self.objective_option == ObjectiveLib.determinant: - model.objective = pyo.Objective( - expr=model.obj_cons.egb_fim_block.outputs["log-D-opt"], - sense=pyo.maximize, - ) - elif self.objective_option == ObjectiveLib.minimum_eigenvalue: - model.objective = pyo.Objective( - expr=model.obj_cons.egb_fim_block.outputs["E-opt"], sense=pyo.maximize - ) - elif self.objective_option == ObjectiveLib.condition_number: + if build_objective: + if self.objective_option == ObjectiveLib.trace: + output_name = "A-opt" + elif self.objective_option == ObjectiveLib.pseudo_trace: + output_name = "pseudo-A-opt" + elif self.objective_option == ObjectiveLib.determinant: + output_name = "log-D-opt" + elif self.objective_option == ObjectiveLib.minimum_eigenvalue: + output_name = "E-opt" + elif self.objective_option == ObjectiveLib.condition_number: + output_name = "ME-opt" + else: + # Error path is intentionally deferred to the external model. + output_name = "A-opt" + model.objective = pyo.Objective( - expr=model.obj_cons.egb_fim_block.outputs["ME-opt"], sense=pyo.minimize + expr=model.obj_cons.egb_fim_block.outputs[output_name], + sense=( + pyo.maximize + if self.objective_option in self._MAXIMIZE_OBJECTIVES + else pyo.minimize + ), ) - # Else error not needed for spurious objective - # options as the error will always appear - # when creating the FIMExternalGreyBox block # Check to see if the model has all the required suffixes def check_model_labels(self, model=None): @@ -1723,31 +3600,72 @@ def check_model_labels(self, model=None): outputs = [k.name for k, v in model.experiment_outputs.items()] except: raise RuntimeError( - "Experiment model does not have suffix " + '"experiment_outputs".' + "Experiment model does not have suffix 'experiment_outputs'." + ) + + # Check that experiment_outputs is not empty + if len(outputs) == 0: + raise ValueError( + "No experiment outputs found. Design of Experiments requires at least " + "one experiment output (measurement) to optimize. Please add an " + "'experiment_outputs' Suffix to your model with at least one variable." ) # Check that experimental inputs exist try: - outputs = [k.name for k, v in model.experiment_inputs.items()] + inputs = [k.name for k, v in model.experiment_inputs.items()] except: raise RuntimeError( - "Experiment model does not have suffix " + '"experiment_inputs".' + "Experiment model does not have suffix 'experiment_inputs'." + ) + + # Check that experiment_inputs is not empty + if len(inputs) == 0: + raise ValueError( + "No experiment inputs found. Design of Experiments requires at least " + "one experiment input (design variable) to optimize. Please add an " + "'experiment_inputs' Suffix to your model with at least one variable." ) # Check that unknown parameters exist try: - outputs = [k.name for k, v in model.unknown_parameters.items()] + parameters = [k.name for k, v in model.unknown_parameters.items()] except: raise RuntimeError( - "Experiment model does not have suffix " + '"unknown_parameters".' + "Experiment model does not have suffix 'unknown_parameters'." + ) + + # Check that unknown_parameters is not empty + if len(parameters) == 0: + raise ValueError( + "No unknown parameters found. Design of Experiments requires at least " + "one unknown parameter to estimate. Please add an " + "'unknown_parameters' Suffix to your model with at least one variable." ) # Check that measurement errors exist try: - outputs = [k.name for k, v in model.measurement_error.items()] + errors = [k.name for k, v in model.measurement_error.items()] except: raise RuntimeError( - "Experiment model does not have suffix " + '"measurement_error".' + "Experiment model does not have suffix 'measurement_error'." + ) + + # Check that measurement_error is not empty + if len(errors) == 0: + raise ValueError( + "No measurement errors found. Design of Experiments requires at least " + "one measurement error specification. Please add a " + "'measurement_error' Suffix to your model with at least one variable." + ) + + # Check that measurement_error and experiment_outputs have the same length + if len(errors) != len(outputs): + raise ValueError( + "Number of experiment outputs, {}, and length of measurement error, " + "{}, do not match. Please check model labeling.".format( + len(outputs), len(errors) + ) ) self.logger.info("Model has expected labels.") @@ -1895,9 +3813,11 @@ def compute_FIM_full_factorial( self.logger.info("Beginning Full Factorial Design.") # Make new model for factorial design - self.factorial_model = self.experiment.get_labeled_model( - **self.get_labeled_model_args - ).clone() + self.factorial_model = ( + self.experiment_list[0] + .get_labeled_model(**self.get_labeled_model_args) + .clone() + ) model = self.factorial_model # Permute the inputs to be aligned with the experiment input indices @@ -2670,7 +4590,7 @@ def get_experiment_input_values(self, model=None): model = self.model if not hasattr(model, "experiment_inputs"): - if not hasattr(model, "scenario_blocks"): + if not hasattr(model, "fd_scenario_blocks"): raise RuntimeError( "Model provided does not have expected structure. " "Please make sure model is built properly before " @@ -2679,7 +4599,7 @@ def get_experiment_input_values(self, model=None): d_vals = [ pyo.value(k) - for k, v in model.scenario_blocks[0].experiment_inputs.items() + for k, v in model.fd_scenario_blocks[0].experiment_inputs.items() ] else: d_vals = [pyo.value(k) for k, v in model.experiment_inputs.items()] @@ -2707,7 +4627,7 @@ def get_unknown_parameter_values(self, model=None): model = self.model if not hasattr(model, "unknown_parameters"): - if not hasattr(model, "scenario_blocks"): + if not hasattr(model, "fd_scenario_blocks"): raise RuntimeError( "Model provided does not have expected structure. Please make " "sure model is built properly before calling " @@ -2716,7 +4636,7 @@ def get_unknown_parameter_values(self, model=None): theta_vals = [ pyo.value(k) - for k, v in model.scenario_blocks[0].unknown_parameters.items() + for k, v in model.fd_scenario_blocks[0].unknown_parameters.items() ] else: theta_vals = [pyo.value(k) for k, v in model.unknown_parameters.items()] @@ -2743,7 +4663,7 @@ def get_experiment_output_values(self, model=None): model = self.model if not hasattr(model, "experiment_outputs"): - if not hasattr(model, "scenario_blocks"): + if not hasattr(model, "fd_scenario_blocks"): raise RuntimeError( "Model provided does not have expected structure. Please make " "sure model is built properly before calling " @@ -2752,7 +4672,7 @@ def get_experiment_output_values(self, model=None): y_hat_vals = [ pyo.value(k) - for k, v in model.scenario_blocks[0].experiment_outputs.items() + for k, v in model.fd_scenario_blocks[0].experiment_outputs.items() ] else: y_hat_vals = [pyo.value(k) for k, v in model.experiment_outputs.items()] @@ -2781,7 +4701,7 @@ def get_measurement_error_values(self, model=None): model = self.model if not hasattr(model, "measurement_error"): - if not hasattr(model, "scenario_blocks"): + if not hasattr(model, "fd_scenario_blocks"): raise RuntimeError( "Model provided does not have expected structure. Please make " "sure model is built properly before calling " @@ -2790,7 +4710,7 @@ def get_measurement_error_values(self, model=None): sigma_vals = [ pyo.value(k) - for k, v in model.scenario_blocks[0].measurement_error.items() + for k, v in model.fd_scenario_blocks[0].measurement_error.items() ] else: sigma_vals = [pyo.value(k) for k, v in model.measurement_error.items()] diff --git a/pyomo/contrib/doe/grey_box_utilities.py b/pyomo/contrib/doe/grey_box_utilities.py index 84eb6187a71..9f09ccedc4c 100644 --- a/pyomo/contrib/doe/grey_box_utilities.py +++ b/pyomo/contrib/doe/grey_box_utilities.py @@ -44,7 +44,14 @@ class FIMExternalGreyBox( ExternalGreyBoxModel if (scipy_available and numpy_available) else object ): - def __init__(self, doe_object, objective_option="determinant", logger_level=None): + def __init__( + self, + doe_object=None, + objective_option="determinant", + logger_level=None, + parameter_names=None, + fim_initial=None, + ): """ Grey box model for metrics on the FIM. This methodology reduces numerical complexity for the computation of FIM metrics related @@ -68,17 +75,35 @@ def __init__(self, doe_object, objective_option="determinant", logger_level=None default: None, or equivalently, use the logging level of doe_object. NOTE: Use logging.DEBUG for all messages. + parameter_names: + Optional ordered iterable of parameter labels. When provided, this + lets the grey box object operate on any FIM source with the same + ordering instead of assuming the data must come from + ``doe_object.model.parameter_names``. This is needed for the + multi-experiment grey box path because the linked FIM lives on a + scenario block (``scenario.total_fim``), while ``doe_object.model`` + is the top-level container and does not own ``parameter_names`` + directly. + fim_initial: + Optional dense, symmetric FIM used to seed the grey box inputs. This + is required when ``doe_object`` is not provided. """ - if doe_object is None: + if doe_object is None and (parameter_names is None or fim_initial is None): raise ValueError( - "DoE Object must be provided to build external grey box of the FIM." + "Either ``doe_object`` or both ``parameter_names`` and " + "``fim_initial`` must be provided to build the FIM grey box." ) self.doe_object = doe_object - # Grab parameter list from the doe_object model - self._param_names = [i for i in self.doe_object.model.parameter_names] + # Grab parameter ordering from the explicit arguments when available. + # Multi-experiment optimization passes the aggregated scenario FIM + # directly, so we should not assume the linked FIM always shares the + # same location as self.doe_object.model. + if parameter_names is None: + parameter_names = self.doe_object.model.parameter_names + self._param_names = [i for i in parameter_names] self._n_params = len(self._param_names) # Check if the doe_object has model components that are required @@ -93,15 +118,21 @@ def __init__(self, doe_object, objective_option="determinant", logger_level=None # If logger level is None, use doe_object's logger level if logger_level is None: - logger_level = doe_object.logger.level - + if doe_object is not None: + logger_level = doe_object.logger.level + else: + logger_level = logging.WARNING self.logger.setLevel(level=logger_level) # Set initial values for inputs # Need a mask structure - self._masking_matrix = np.triu(np.ones_like(self.doe_object.fim_initial)) + if fim_initial is None: + fim_initial = self.doe_object.fim_initial + fim_initial = np.asarray(fim_initial, dtype=np.float64) + + self._masking_matrix = np.triu(np.ones_like(fim_initial)) self._input_values = np.asarray( - self.doe_object.fim_initial[self._masking_matrix > 0], dtype=np.float64 + fim_initial[self._masking_matrix > 0], dtype=np.float64 ) self._n_inputs = len(self._input_values) diff --git a/pyomo/contrib/doe/tests/experiment_class_example_flags.py b/pyomo/contrib/doe/tests/experiment_class_example_flags.py index c6b3cb9e0ec..38cea05ee4a 100644 --- a/pyomo/contrib/doe/tests/experiment_class_example_flags.py +++ b/pyomo/contrib/doe/tests/experiment_class_example_flags.py @@ -70,3 +70,80 @@ def get_labeled_model(self): m.bad_con_2 = pyo.Constraint(expr=m.hour <= 0.0) return m + + +class RooneyBieglerMultiExperiment(RooneyBieglerExperiment): + """ + Experiment class based on the multi-experiment Rooney-Biegler prototype. + + This mirrors the implementation in + ``examples/multiexperiment-prototype/rooney_biegler_multiexperiment.py`` + while allowing test-time control over initial hour and bounds. + """ + + def __init__( + self, hour=2.0, y=10.0, theta=None, measure_error=0.1, hour_bounds=(1.0, 10.0) + ): + data = {'hour': hour, 'y': y} + super().__init__(data=data, measure_error=measure_error, theta=theta) + self.hour_bounds = hour_bounds + + def get_labeled_model(self): + m = super().get_labeled_model() + hour_lb, hour_ub = self.hour_bounds + m.hour.setlb(hour_lb) + m.hour.setub(hour_ub) + + if hasattr(m, "sym_break_cons"): + m.sym_break_cons.clear() + else: + m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.sym_break_cons[m.hour] = None + return m + + +class RooneyBieglerMultiInputExperimentFlag(RooneyBieglerExperiment): + """ + Two-input Rooney-Biegler style experiment for symmetry-breaking tests. + + Parameters + ---------- + sym_break_flag : int + 0 -> do not add ``sym_break_cons`` suffix + 1 -> add single marker (hour) + 2 -> add multiple markers (hour and temp) + """ + + def __init__(self, hour=2.0, temp=300.0, y=10.0, sym_break_flag=1): + data = {'hour': hour, 'y': y} + super().__init__( + data=data, measure_error=0.1, theta={'asymptote': 15, 'rate_constant': 0.5} + ) + self.hour = hour + self.temp = temp + self.sym_break_flag = sym_break_flag + + def get_labeled_model(self): + m = super().get_labeled_model() + + m.hour.setlb(1.0) + m.hour.setub(10.0) + m.temp = pyo.Var(initialize=self.temp, bounds=(280.0, 340.0)) + m.temp.fix() + + # Replace base Rooney-Biegler response with two-input synthetic variant + # used only for symmetry-breaking tests. + m.del_component(m.response_function) + m.response_function = pyo.Constraint( + expr=m.y + == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.hour)) + 0.01 * m.temp + ) + m.experiment_inputs[m.temp] = self.temp + + if self.sym_break_flag in (1, 2): + m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.sym_break_cons[m.hour] = None + if self.sym_break_flag == 2: + m.sym_break_cons[m.temp] = None + + return m diff --git a/pyomo/contrib/doe/tests/test_doe_build.py b/pyomo/contrib/doe/tests/test_doe_build.py index 333ac1669ef..7e821e91408 100644 --- a/pyomo/contrib/doe/tests/test_doe_build.py +++ b/pyomo/contrib/doe/tests/test_doe_build.py @@ -7,7 +7,10 @@ # software. This software is distributed under the 3-clause BSD License. # ____________________________________________________________________________________ import json +import os import os.path +import tempfile +from pathlib import Path from pyomo.common.dependencies import ( numpy as np, @@ -19,15 +22,23 @@ from pyomo.common.fileutils import this_file_dir import pyomo.common.unittest as unittest +from pyomo.contrib.doe.doe import ObjectiveLib, _DoEResultsJSONEncoder if not (numpy_available and scipy_available): raise unittest.SkipTest("Pyomo.DoE needs scipy and numpy to run tests") if scipy_available: from pyomo.contrib.doe import DesignOfExperiments + from pyomo.contrib.doe.examples.reactor_example import ( + ReactorExperiment as FullReactorExperiment, + ) + from pyomo.contrib.doe.tests.experiment_class_example_flags import ( + RooneyBieglerMultiExperiment, + ) from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( RooneyBieglerExperiment, ) +from pyomo.contrib.doe.tests.utils_for_doe_tests import make_ipopt_solver from pyomo.contrib.doe.examples.rooney_biegler_doe_example import run_rooney_biegler_doe import pyomo.environ as pyo @@ -45,6 +56,26 @@ data_ex["control_points"] = {float(k): v for k, v in data_ex["control_points"].items()} +class _TrackingSolver: + """Proxy solver that records solve-call metadata while running real solves.""" + + def __init__(self, solver, phase_tag, call_order): + self._solver = solver + self._phase_tag = phase_tag + self._call_order = call_order + self.calls = 0 + self.option_markers = [] + + def solve(self, *args, **kwargs): + self.calls += 1 + self._call_order.append(self._phase_tag) + self.option_markers.append(self._solver.options.get("max_iter")) + return self._solver.solve(*args, **kwargs) + + def __getattr__(self, name): + return getattr(self._solver, name) + + def get_rooney_biegler_experiment(): """Get a fresh RooneyBieglerExperiment instance for testing. @@ -94,7 +125,13 @@ def get_FIM_FIMPrior_Q_L(doe_obj=None): for i in model.output_names for j in model.parameter_names ] - sigma_inv = [1 / v for k, v in model.scenario_blocks[0].measurement_error.items()] + sigma_inv = [ + 1 / v for k, v in model.fd_scenario_blocks[0].measurement_error.items() + ] + param_vals = np.array( + [[v for k, v in model.fd_scenario_blocks[0].unknown_parameters.items()]] + ) + FIM_vals_np = np.array(FIM_vals).reshape((n_param, n_param)) FIM_prior_vals_np = np.array(FIM_prior_vals).reshape((n_param, n_param)) @@ -116,7 +153,7 @@ def get_FIM_FIMPrior_Q_L(doe_obj=None): def get_standard_args(experiment, fd_method, obj_used): args = {} - args['experiment'] = experiment + args['experiment'] = None if experiment is None else [experiment] args['fd_formula'] = fd_method args['step'] = 1e-3 args['objective_option'] = obj_used @@ -126,13 +163,8 @@ def get_standard_args(experiment, fd_method, obj_used): args['jac_initial'] = None args['fim_initial'] = None args['L_diagonal_lower_bound'] = 1e-7 - # Make solver object with - # good linear subroutines - solver = SolverFactory("ipopt") - solver.options["linear_solver"] = "ma57" - solver.options["halt_on_ampl_error"] = "yes" - solver.options["max_iter"] = 3000 - args['solver'] = solver + # Make solver object with good linear subroutines. + args['solver'] = make_ipopt_solver() args['tee'] = False args['get_labeled_model_args'] = None args['_Cholesky_option'] = True @@ -145,6 +177,20 @@ def get_standard_args(experiment, fd_method, obj_used): @unittest.skipIf(not scipy_available, "scipy is not available") @unittest.skipIf(not pandas_available, "pandas is not available") class TestDoeBuild(unittest.TestCase): + def test_constructor_accepts_single_experiment_or_list(self): + # The public constructor should normalize either form into the same + # internal experiment_list representation. + single = DesignOfExperiments( + experiment=RooneyBieglerMultiExperiment(hour=2.0), + objective_option="pseudo_trace", + ) + as_list = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + ) + self.assertEqual(len(single.experiment_list), 1) + self.assertEqual(len(as_list.experiment_list), 1) + def test_rooney_biegler_fd_central_check_fd_eqns(self): fd_method = "central" obj_used = "pseudo_trace" @@ -166,16 +212,16 @@ def test_rooney_biegler_fd_central_check_fd_eqns(self): diff = (-1) ** s * doe_obj.step param_val = pyo.value( - pyo.ComponentUID(param).find_component_on(model.scenario_blocks[s]) + pyo.ComponentUID(param).find_component_on(model.fd_scenario_blocks[s]) ) - param_val_from_step = model.scenario_blocks[0].unknown_parameters[ - pyo.ComponentUID(param).find_component_on(model.scenario_blocks[0]) + param_val_from_step = model.fd_scenario_blocks[0].unknown_parameters[ + pyo.ComponentUID(param).find_component_on(model.fd_scenario_blocks[0]) ] * (1 + diff) - for k, v in model.scenario_blocks[s].unknown_parameters.items(): + for k, v in model.fd_scenario_blocks[s].unknown_parameters.items(): if pyo.ComponentUID( - k, context=model.scenario_blocks[s] + k, context=model.fd_scenario_blocks[s] ) == pyo.ComponentUID(param): continue @@ -205,17 +251,21 @@ def test_rooney_biegler_fd_backward_check_fd_eqns(self): param = model.parameter_scenarios[s] param_val = pyo.value( - pyo.ComponentUID(param).find_component_on(model.scenario_blocks[s]) + pyo.ComponentUID(param).find_component_on( + model.fd_scenario_blocks[s] + ) ) - param_val_from_step = model.scenario_blocks[0].unknown_parameters[ - pyo.ComponentUID(param).find_component_on(model.scenario_blocks[0]) + param_val_from_step = model.fd_scenario_blocks[0].unknown_parameters[ + pyo.ComponentUID(param).find_component_on( + model.fd_scenario_blocks[0] + ) ] * (1 + diff) self.assertAlmostEqual(param_val, param_val_from_step) - for k, v in model.scenario_blocks[s].unknown_parameters.items(): + for k, v in model.fd_scenario_blocks[s].unknown_parameters.items(): if (s != 0) and pyo.ComponentUID( - k, context=model.scenario_blocks[s] + k, context=model.fd_scenario_blocks[s] ) == pyo.ComponentUID(param): continue @@ -243,17 +293,21 @@ def test_rooney_biegler_fd_forward_check_fd_eqns(self): param = model.parameter_scenarios[s] param_val = pyo.value( - pyo.ComponentUID(param).find_component_on(model.scenario_blocks[s]) + pyo.ComponentUID(param).find_component_on( + model.fd_scenario_blocks[s] + ) ) - param_val_from_step = model.scenario_blocks[0].unknown_parameters[ - pyo.ComponentUID(param).find_component_on(model.scenario_blocks[0]) + param_val_from_step = model.fd_scenario_blocks[0].unknown_parameters[ + pyo.ComponentUID(param).find_component_on( + model.fd_scenario_blocks[0] + ) ] * (1 + diff) self.assertAlmostEqual(param_val, param_val_from_step) - for k, v in model.scenario_blocks[s].unknown_parameters.items(): + for k, v in model.fd_scenario_blocks[s].unknown_parameters.items(): if (s != 0) and pyo.ComponentUID( - k, context=model.scenario_blocks[s] + k, context=model.fd_scenario_blocks[s] ) == pyo.ComponentUID(param): continue @@ -275,7 +329,9 @@ def test_rooney_biegler_fd_central_design_fixing(self): model = doe_obj.model # Check that the design fixing constraints are generated - design_vars = [k for k, v in model.scenario_blocks[0].experiment_inputs.items()] + design_vars = [ + k for k, v in model.fd_scenario_blocks[0].experiment_inputs.items() + ] con_name_base = "global_design_eq_con_" @@ -308,7 +364,9 @@ def test_rooney_biegler_fd_backward_design_fixing(self): model = doe_obj.model # Check that the design fixing constraints are generated - design_vars = [k for k, v in model.scenario_blocks[0].experiment_inputs.items()] + design_vars = [ + k for k, v in model.fd_scenario_blocks[0].experiment_inputs.items() + ] con_name_base = "global_design_eq_con_" @@ -341,7 +399,9 @@ def test_rooney_biegler_fd_forward_design_fixing(self): model = doe_obj.model # Check that the design fixing constraints are generated - design_vars = [k for k, v in model.scenario_blocks[0].experiment_inputs.items()] + design_vars = [ + k for k, v in model.fd_scenario_blocks[0].experiment_inputs.items() + ] con_name_base = "global_design_eq_con_" @@ -495,11 +555,11 @@ def test_generate_blocks_without_model(self): doe_obj = DesignOfExperiments(**DoE_args) - doe_obj._generate_scenario_blocks() + doe_obj._generate_fd_scenario_blocks() for i in doe_obj.model.parameter_scenarios: self.assertTrue( - doe_obj.model.find_component("scenario_blocks[" + str(i) + "]") + doe_obj.model.find_component("fd_scenario_blocks[" + str(i) + "]") ) @@ -520,6 +580,18 @@ def test_reactor_update_suffix_items(self): @unittest.skipIf(not numpy_available, "Numpy is not available") @unittest.skipIf(not pandas_available, "pandas is not available") class TestDoEObjectiveOptions(unittest.TestCase): + def test_maximize_objective_set_contents(self): + # The objective-sense partition drives maximize/minimize scoring logic + # in initialization and solve helpers, so keep a direct regression on + # the public enum membership of the maximize set. + maximize = DesignOfExperiments._MAXIMIZE_OBJECTIVES + self.assertIn(ObjectiveLib.determinant, maximize) + self.assertIn(ObjectiveLib.pseudo_trace, maximize) + self.assertIn(ObjectiveLib.minimum_eigenvalue, maximize) + self.assertNotIn(ObjectiveLib.trace, maximize) + self.assertNotIn(ObjectiveLib.condition_number, maximize) + self.assertNotIn(ObjectiveLib.zero, maximize) + def test_trace_constraints(self): fd_method = "central" obj_used = "trace" @@ -562,7 +634,7 @@ def test_trace_constraints(self): for i, c in enumerate(params): for j, d in enumerate(params): - # cholesky_inv_imp: only defined for i >= j + # inverse constraint only exists for lower triangle (i >= j) if i >= j: self.assertIn( (c, d), @@ -576,7 +648,7 @@ def test_trace_constraints(self): msg=f"Unexpected cholesky_inv_cons[{c},{d}]", ) - # cholesky_LLinv_imp: only defined for i >= j + # identity constraint only defined for lower triangle (i >= j) if i >= j: self.assertIn( (c, d), @@ -633,5 +705,404 @@ def test_trace_initialization_consistency(self): self.assertAlmostEqual(val, expected, places=4) +@unittest.skipIf(not pandas_available, "pandas is not available") +class TestOptimizeExperimentsBuildStructure(unittest.TestCase): + """Coverage for optimize_experiments() build, output, and diagnostics behavior.""" + + def _make_solver(self): + return make_ipopt_solver() + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_init_solver_used_for_initialization_only(self): + # Tests that all pre-final solves use init_solver while the final + # optimization solve still uses the primary solver. + main_solver_inner = self._make_solver() + init_solver_inner = self._make_solver() + # Use distinct option values so each solver path can be identified. + main_solver_inner.options["max_iter"] = 321 + init_solver_inner.options["max_iter"] = 123 + call_order = [] + main_solver = _TrackingSolver( + solver=main_solver_inner, phase_tag="main", call_order=call_order + ) + init_solver = _TrackingSolver( + solver=init_solver_inner, phase_tag="init", call_order=call_order + ) + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=main_solver, + ) + + doe_obj.optimize_experiments(n_exp=2, init_solver=init_solver) + + # The exact number of initialization solves is implementation-dependent, + # but they must all occur before the one final main-solver call. + self.assertGreaterEqual( + init_solver.calls, 1 + ) # At least one initialization solve + self.assertEqual(main_solver.calls, 1) # Exactly one main optimization solve + self.assertEqual(call_order[-1], "main") + self.assertTrue(all(tag == "init" for tag in call_order[:-1])) + # Distinct option markers provide a second check that solver routing + # matches the expected init-versus-final phase split. + option_markers = init_solver.option_markers + main_solver.option_markers + self.assertTrue(all(marker == 123 for marker in option_markers[:-1])) + self.assertEqual(option_markers[-1], 321) + # Result payloads should report the same phase-specific solver names that + # were observed through the tracked solve() calls above. + self.assertEqual( + doe_obj.results["initialization"]["solver"], + getattr(init_solver, "name", str(init_solver)), + ) + self.assertEqual( + doe_obj.results["optimization_solve"]["solver"], + getattr(main_solver, "name", str(main_solver)), + ) + + def test_get_experiment_input_vars_direct_and_fd_fallback(self): + # Test the helper used by optimize_experiments() finds input vars for both + # direct models and finite-difference scenario-block models. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + ) + + model_direct = pyo.ConcreteModel() + model_direct.x = pyo.Var(initialize=2.0, bounds=(1.0, 10.0)) + model_direct.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + model_direct.experiment_inputs[model_direct.x] = 2.0 + vars_direct = doe_obj._get_experiment_input_vars(model_direct) + self.assertEqual([v.name for v in vars_direct], [model_direct.x.name]) + + model_fd = pyo.ConcreteModel() + model_fd.fd_scenario_blocks = pyo.Block([0]) + model_fd.fd_scenario_blocks[0].x = pyo.Var(initialize=3.0, bounds=(1.0, 10.0)) + model_fd.fd_scenario_blocks[0].experiment_inputs = pyo.Suffix( + direction=pyo.Suffix.LOCAL + ) + model_fd.fd_scenario_blocks[0].experiment_inputs[ + model_fd.fd_scenario_blocks[0].x + ] = 3.0 + vars_fd = doe_obj._get_experiment_input_vars(model_fd) + self.assertEqual( + [v.name for v in vars_fd], [model_fd.fd_scenario_blocks[0].x.name] + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_multi_experiment_structure_and_results(self): + # Test that the multi-experiment optimize_experiments() run builds the expected + # scenario/experiment structure and structured results. + solver = self._make_solver() + + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiExperiment(hour=2.0), + RooneyBieglerMultiExperiment(hour=4.0), + ], + objective_option="pseudo_trace", + step=1e-2, + solver=solver, + ) + doe_obj.optimize_experiments() + + scenario_block = doe_obj.model.param_scenario_blocks[0] + self.assertTrue(hasattr(scenario_block, "symmetry_breaking_s0_exp1")) + self.assertEqual(len(list(scenario_block.exp_blocks.keys())), 2) + + param_scenario = doe_obj.results["solution"]["param_scenarios"][0] + self.assertEqual(doe_obj.results["initialization"]["method"], "none") + self.assertEqual( + doe_obj.results["problem"]["number_of_experiments_per_scenario"], 2 + ) + self.assertEqual(len(param_scenario["experiments"]), 2) + self.assertEqual(len(doe_obj.results["problem"]["design_variables"]), 1) + self.assertEqual(len(doe_obj.results["problem"]["parameters"]), 2) + + # Results should expose a single structured parameter-scenario payload. + self.assertIn("problem", doe_obj.results) + self.assertIn("initialization", doe_obj.results) + self.assertIn("optimization_solve", doe_obj.results) + self.assertIn("timing", doe_obj.results) + self.assertIn("solution", doe_obj.results) + self.assertEqual(doe_obj.results["solution"]["objective"], "pseudo_trace") + self.assertEqual(doe_obj.results["optimization_solve"]["status"], "ok") + self.assertFalse(doe_obj.results["problem"]["used_template_experiment"]) + self.assertEqual(len(doe_obj.results["solution"]["param_scenarios"]), 1) + self.assertEqual(param_scenario["param_scenario_id"], 0) + self.assertEqual(param_scenario["param_scenario_weight"], 1.0) + self.assertEqual(len(param_scenario["experiments"]), 2) + self.assertEqual(param_scenario["experiments"][0]["exp_id"], 0) + self.assertEqual( + doe_obj.results["problem"]["measurement_error_values"], + doe_obj.get_measurement_error_values(model=scenario_block.exp_blocks[0]), + ) + self.assertNotIn("measurement_errors", param_scenario["experiments"][0]) + + # hour of exp[0] should be <= hour of exp[1] due to symmetry breaking + h0 = param_scenario["experiments"][0]["design"][0] + h1 = param_scenario["experiments"][1]["design"][0] + self.assertLessEqual(h0, h1) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_writes_results_file(self): + # Tests that optimize_experiments() writes JSON results when given either + # a string path or a pathlib.Path for results_file. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + fd, results_path = tempfile.mkstemp(suffix=".json") + os.close(fd) + self.addCleanup( + lambda: os.path.exists(results_path) and os.remove(results_path) + ) + + doe_obj.optimize_experiments(n_exp=1, results_file=results_path) + + with open(results_path) as f: + payload = json.load(f) + self.assertEqual(payload["initialization"]["method"], "none") + self.assertTrue(payload["problem"]["used_template_experiment"]) + self.assertIn("solution", payload) + self.assertEqual(payload["solution"]["objective"], "pseudo_trace") + self.assertIn("optimization_solve", payload) + self.assertIn("problem", payload) + self.assertIn("timing", payload) + path_payload = Path(results_path) + doe_obj.optimize_experiments(n_exp=1, results_file=path_payload) + with open(results_path) as f: + payload_path = json.load(f) + self.assertEqual(payload_path["initialization"]["method"], "none") + self.assertTrue(payload_path["problem"]["used_template_experiment"]) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_single_experiment_defaults_to_template_mode(self): + # Tests that optimize_experiments() uses template mode by default when + # n_exp=1. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + doe_obj.optimize_experiments() + self.assertEqual( + doe_obj.results["problem"]["number_of_experiments_per_scenario"], 1 + ) + self.assertTrue(doe_obj.results["problem"]["used_template_experiment"]) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_zero_objective_works_without_obj_cons(self): + # The zero objective intentionally builds no scenario.obj_cons block, so + # optimize_experiments() must skip the deactivate/reactivate cycle and + # still produce the usual results payload from the square solve. + init_solver = self._make_solver() + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="zero", + step=1e-2, + solver=self._make_solver(), + ) + doe_obj.optimize_experiments(n_exp=2, init_solver=init_solver) + + scenario = doe_obj.model.param_scenario_blocks[0] + self.assertFalse(hasattr(scenario, "obj_cons")) + self.assertIn( + doe_obj.results["optimization_solve"]["status"], {"ok", "warning"} + ) + self.assertTrue( + isinstance(doe_obj.results["optimization_solve"]["message"], str) + ) + self.assertEqual( + len(doe_obj.results["solution"]["param_scenarios"][0]["experiments"]), 2 + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_trace_roundoff_flag_builds_extra_constraints(self): + # The multi-experiment trace path optionally adds extra Cholesky/FIM + # diagonal constraints to reduce roundoff drift. + init_solver = self._make_solver() + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="trace", + step=1e-2, + solver=self._make_solver(), + improve_cholesky_roundoff_error=True, + ) + doe_obj.optimize_experiments(n_exp=2, init_solver=init_solver) + + scenario = doe_obj.model.param_scenario_blocks[0] + parameter_names = list(scenario.exp_blocks[0].parameter_names) + + self.assertIn( + doe_obj.results["optimization_solve"]["status"], {"ok", "warning"} + ) + self.assertTrue(hasattr(scenario.obj_cons, "cholesky_fim_diag_cons")) + self.assertTrue(hasattr(scenario.obj_cons, "cholesky_fim_inv_diag_cons")) + self.assertEqual( + len(scenario.obj_cons.cholesky_fim_diag_cons), len(parameter_names) ** 2 + ) + self.assertEqual( + len(scenario.obj_cons.cholesky_fim_inv_diag_cons), len(parameter_names) ** 2 + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_timing_includes_lhs_phase_separately(self): + # Tests that LHS initialization timing is tracked separately and + # contributes additively to total runtime accounting. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + doe_obj.optimize_experiments( + n_exp=2, + init_method="latin_hypercube_sampling", + init_n_samples=2, + # Set random seed to keep LHS initialization deterministic. + init_seed=11, + ) + + timing = doe_obj.results["timing"] + self.assertIn("lhs_initialization_time_s", timing) + self.assertGreaterEqual(timing["lhs_initialization_time_s"], 0.0) + self.assertAlmostEqual( + timing["total_time_s"], + timing["build_time_s"] + + timing["lhs_initialization_time_s"] + + timing["initialization_time_s"] + + timing["optimization_solve_time_s"], + places=8, + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_symmetry_log_once_per_scenario(self): + # Tests that symmetry-breaking constraint logging is emitted once per + # scenario (not once per generated constraint). + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + with self.assertLogs("pyomo.contrib.doe.doe", level="INFO") as log_cm: + doe_obj.optimize_experiments(n_exp=3) + + matching = [ + m + for m in log_cm.output + if "Added 2 symmetry breaking constraints for scenario 0" in m + ] + self.assertEqual(len(matching), 1) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_lhs_diagnostics_populated(self): + # Tests that LHS initialization records diagnostics needed for debugging + # and performance visibility. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + doe_obj.optimize_experiments( + n_exp=2, + init_method="latin_hypercube_sampling", + init_n_samples=2, + # Set random seed to keep LHS initialization deterministic. + init_seed=11, + ) + lhs_init = doe_obj.results["initialization"] + self.assertEqual(lhs_init["candidate_fim_evaluation_mode"], "serial") + self.assertEqual(lhs_init["combination_scoring_mode"], "serial") + self.assertFalse(lhs_init["timed_out"]) + self.assertGreater(lhs_init["time_s"], 0.0) + self.assertIn("best_initial_objective_value", lhs_init) + self.assertIsInstance(lhs_init["best_initial_objective_value"], float) + self.assertTrue(np.isfinite(lhs_init["best_initial_objective_value"])) + self.assertGreater(lhs_init["best_initial_objective_value"], 0.0) + self.assertIn("best_initial_objective_value_log10", lhs_init) + self.assertIsInstance(lhs_init["best_initial_objective_value_log10"], float) + self.assertAlmostEqual( + lhs_init["best_initial_objective_value_log10"], + np.log10(lhs_init["best_initial_objective_value"]), + places=4, + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_lhs_rooney_biegler_integration_records_consistent_counts( + self, + ): + # End-to-end Rooney-Biegler integration check for LHS initialization: + # selected initial designs should be recorded with consistent candidate/ + # combination counts for n_exp=2 and a single design variable. + lhs_n_samples = 3 + n_exp = 2 + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + doe_obj.optimize_experiments( + n_exp=n_exp, + init_method="latin_hypercube_sampling", + init_n_samples=lhs_n_samples, + # Set random seed to keep LHS initialization deterministic. + init_seed=17, + ) + lhs_init = doe_obj.results["initialization"] + + self.assertEqual(lhs_init["method"], "latin_hypercube_sampling") + self.assertEqual(len(lhs_init["selected_initial_designs"]), n_exp) + for point in lhs_init["selected_initial_designs"]: + self.assertEqual(len(point), 1) + self.assertGreaterEqual(point[0], 1.0) + self.assertLessEqual(point[0], 10.0) + + # For Rooney-Biegler there is one design variable (hour), so: + # n_candidates = lhs_n_samples^1 and n_combinations = C(n_candidates, n_exp). + self.assertEqual(lhs_init["number_of_candidate_designs"], lhs_n_samples) + self.assertEqual(lhs_init["number_of_design_combinations"], 3) + + +@unittest.skipIf(not numpy_available, "Numpy is not available") +class TestDoeResultsSerialization(unittest.TestCase): + """Coverage for DoE results payload serialization helpers.""" + + def test_doe_results_json_encoder_handles_numpy_and_enum(self): + # Results payloads are written to JSON from optimize_experiments(), so + # the encoder must normalize numpy scalars/arrays and ObjectiveLib enums. + payload = { + "scalar": np.int64(7), + "array": np.array([1.0, 2.0]), + "objective": ObjectiveLib.trace, + } + encoded = json.dumps(payload, cls=_DoEResultsJSONEncoder) + decoded = json.loads(encoded) + + self.assertEqual(decoded["scalar"], 7) + self.assertEqual(decoded["array"], [1.0, 2.0]) + self.assertEqual(decoded["objective"], str(ObjectiveLib.trace)) + + +@unittest.skipIf(not numpy_available, "Numpy is not available") +class TestDoeBuildHelpers(unittest.TestCase): + """Coverage for small matrix helpers used by build/solve paths.""" + + def test_symmetrize_lower_tri_helper(self): + # Only the lower-triangular FIM is stored in several code paths; this + # helper must mirror it to a full symmetric matrix without doubling the diagonal. + m = np.array([[1.0, 0.0, 0.0], [2.0, 3.0, 0.0], [4.0, 5.0, 6.0]]) + got = DesignOfExperiments._symmetrize_lower_tri(m) + expected = np.array([[1.0, 2.0, 4.0], [2.0, 3.0, 5.0], [4.0, 5.0, 6.0]]) + self.assertTrue(np.allclose(got, expected, atol=1e-12)) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/doe/tests/test_doe_errors.py b/pyomo/contrib/doe/tests/test_doe_errors.py index b5587760148..a4d40d6e0bc 100644 --- a/pyomo/contrib/doe/tests/test_doe_errors.py +++ b/pyomo/contrib/doe/tests/test_doe_errors.py @@ -6,13 +6,15 @@ # Solutions of Sandia, LLC, the U.S. Government retains certain rights in this # software. This software is distributed under the 3-clause BSD License. # ____________________________________________________________________________________ - +import json +import warnings from pyomo.common.dependencies import ( numpy as np, numpy_available, pandas as pd, pandas_available, scipy_available, + attempt_import, ) from pyomo.common.errors import DeveloperError @@ -21,17 +23,41 @@ if not (numpy_available and scipy_available): raise unittest.SkipTest("Pyomo.DoE needs scipy and numpy to run tests") +from pyomo.contrib.doe import DesignOfExperiments +import pyomo.contrib.doe.doe as doe_module +from pyomo.contrib.doe.doe import InitializationMethod, _DoEResultsJSONEncoder +from pyomo.contrib.doe.tests.experiment_class_example_flags import ( + BadExperiment, + RooneyBieglerExperimentFlag, + RooneyBieglerMultiExperiment, + RooneyBieglerMultiInputExperimentFlag, +) +from pyomo.contrib.doe.tests.utils_for_doe_tests import make_ipopt_solver +from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( + RooneyBieglerExperiment, +) + if scipy_available: from pyomo.contrib.doe import DesignOfExperiments + from pyomo.contrib.doe.doe import InitializationMethod, _DoEResultsJSONEncoder from pyomo.contrib.doe.tests.experiment_class_example_flags import ( BadExperiment, RooneyBieglerExperimentFlag, + RooneyBieglerMultiExperiment, + RooneyBieglerMultiInputExperimentFlag, ) from pyomo.contrib.doe.examples.rooney_biegler_doe_example import run_rooney_biegler_doe +import pyomo.environ as pyo from pyomo.opt import SolverFactory ipopt_available = SolverFactory("ipopt").available() +parameterized, parameterized_available = attempt_import("parameterized") + + +class _DummyExperiment: + def get_labeled_model(self, **kwargs): + raise RuntimeError("Should not be called in argument-validation tests") def get_rooney_biegler_experiment_flag(): @@ -52,7 +78,7 @@ def get_rooney_biegler_experiment_flag(): def get_standard_args(experiment, fd_method, obj_used, flag): args = {} - args['experiment'] = experiment + args['experiment'] = None if experiment is None else [experiment] args['fd_formula'] = fd_method args['step'] = 1e-3 args['objective_option'] = obj_used @@ -62,13 +88,8 @@ def get_standard_args(experiment, fd_method, obj_used, flag): args['jac_initial'] = None args['fim_initial'] = None args['L_diagonal_lower_bound'] = 1e-7 - # Make solver object with - # good linear subroutines - solver = SolverFactory("ipopt") - solver.options["linear_solver"] = "ma57" - solver.options["halt_on_ampl_error"] = "yes" - solver.options["max_iter"] = 3000 - args['solver'] = solver + # Make solver object with good linear subroutines. + args['solver'] = make_ipopt_solver() args['tee'] = False if flag is not None: args['get_labeled_model_args'] = {"flag": flag} @@ -81,19 +102,36 @@ def get_standard_args(experiment, fd_method, obj_used, flag): @unittest.skipIf(not scipy_available, "scipy is not available") @unittest.skipIf(not pandas_available, "pandas is not available") class TestDoEErrors(unittest.TestCase): + def _make_dummy_optimize_experiments_doe(self, n_experiments=1): + return DesignOfExperiments( + experiment=[_DummyExperiment() for _ in range(n_experiments)], + objective_option="pseudo_trace", + ) + def test_experiment_none_error(self): fd_method = "central" obj_used = "pseudo_trace" flag_val = 1 # Value for faulty model build mode - 1: No exp outputs with self.assertRaisesRegex( - ValueError, "Experiment object must be provided to perform DoE." + ValueError, "The 'experiment' argument is required" ): # Experiment provided as None DoE_args = get_standard_args(None, fd_method, obj_used, flag_val) DesignOfExperiments(**DoE_args) + def test_experiment_empty_list_error(self): + with self.assertRaisesRegex( + ValueError, + "The 'experiment' argument is required and cannot be an empty list", + ): + DesignOfExperiments(experiment=[], objective_option="pseudo_trace") + + def test_doe_results_json_encoder_unsupported_object_raises(self): + with self.assertRaises(TypeError): + json.dumps({"x": object()}, cls=_DoEResultsJSONEncoder) + def test_reactor_check_no_get_labeled_model(self): fd_method = "central" obj_used = "pseudo_trace" @@ -102,8 +140,7 @@ def test_reactor_check_no_get_labeled_model(self): experiment = BadExperiment() with self.assertRaisesRegex( - ValueError, - "The experiment object must have a ``get_labeled_model`` function", + ValueError, "Experiment at index .* must have a.*get_labeled_model" ): DoE_args = get_standard_args(experiment, fd_method, obj_used, flag_val) @@ -121,8 +158,7 @@ def test_reactor_check_no_experiment_outputs(self): doe_obj = DesignOfExperiments(**DoE_args) with self.assertRaisesRegex( - RuntimeError, - "Experiment model does not have suffix " + '"experiment_outputs".', + RuntimeError, "Experiment model does not have suffix 'experiment_outputs'." ): doe_obj.create_doe_model() @@ -138,8 +174,7 @@ def test_reactor_check_no_measurement_error(self): doe_obj = DesignOfExperiments(**DoE_args) with self.assertRaisesRegex( - RuntimeError, - "Experiment model does not have suffix " + '"measurement_error".', + RuntimeError, "Experiment model does not have suffix 'measurement_error'." ): doe_obj.create_doe_model() @@ -155,8 +190,7 @@ def test_reactor_check_no_experiment_inputs(self): doe_obj = DesignOfExperiments(**DoE_args) with self.assertRaisesRegex( - RuntimeError, - "Experiment model does not have suffix " + '"experiment_inputs".', + RuntimeError, "Experiment model does not have suffix 'experiment_inputs'." ): doe_obj.create_doe_model() @@ -172,8 +206,7 @@ def test_reactor_check_no_unknown_parameters(self): doe_obj = DesignOfExperiments(**DoE_args) with self.assertRaisesRegex( - RuntimeError, - "Experiment model does not have suffix " + '"unknown_parameters".', + RuntimeError, "Experiment model does not have suffix 'unknown_parameters'." ): doe_obj.create_doe_model() @@ -674,7 +707,7 @@ def test_bad_FD_generate_scens(self): "Please report this to the Pyomo Developers.", ): doe_obj.fd_formula = "bad things" - doe_obj._generate_scenario_blocks() + doe_obj._generate_fd_scenario_blocks() @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") def test_bad_FD_seq_compute_FIM(self): @@ -759,6 +792,31 @@ def test_bad_compute_FIM_option(self): ): doe_obj.compute_FIM(method="Bad Method") + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_compute_FIM_multi_experiment_parameter_value_mismatch(self): + fd_method = "central" + obj_used = "pseudo_trace" + + DoE_args = get_standard_args( + RooneyBieglerMultiExperiment(hour=1.5, y=9.0), fd_method, obj_used, None + ) + DoE_args["experiment"] = [ + RooneyBieglerMultiExperiment( + hour=1.5, y=9.0, theta={'asymptote': 15, 'rate_constant': 0.5} + ), + RooneyBieglerMultiExperiment( + hour=3.5, y=12.0, theta={'asymptote': 16, 'rate_constant': 0.5} + ), + ] + doe_obj = DesignOfExperiments(**DoE_args) + + # The mismatch is detected before the second experiment solve, when + # compute_FIM validates unknown parameter values across experiments. + with self.assertRaisesRegex( + ValueError, "must share the same unknown parameter values" + ): + doe_obj.compute_FIM(method="sequential") + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") def test_invalid_trace_without_cholesky(self): fd_method = "central" @@ -774,10 +832,442 @@ def test_invalid_trace_without_cholesky(self): with self.assertRaisesRegex( ValueError, - "objective_option='trace' currently only implemented with ``_Cholesky option=True``.", + "objective_option='trace' currently only implemented with " + "``_Cholesky option=True``.", ): doe_obj.create_objective_function() + @unittest.skipIf( + not parameterized_available, "The 'parameterized' package is not available" + ) + @parameterized.parameterized.expand( + [ + ( + "unsupported_init_method", + {"init_method": "bad"}, + ValueError, + r"``init_method`` must be one of \[None, 'latin_hypercube_sampling'\], got 'bad'.", + ), + ( + "enum_init_method_still_validates_init_n_samples", + { + "init_method": InitializationMethod.latin_hypercube_sampling, + "init_n_samples": 0, + }, + ValueError, + r"``init_n_samples`` must be a positive integer, got 0.", + ), + ( + "non_positive_init_n_samples", + {"init_method": "latin_hypercube_sampling", "init_n_samples": 0}, + ValueError, + r"``init_n_samples`` must be a positive integer, got 0.", + ), + ( + "non_integer_init_n_samples", + {"init_method": "latin_hypercube_sampling", "init_n_samples": 2.5}, + ValueError, + r"``init_n_samples`` must be a positive integer, got 2.5.", + ), + ( + "init_seed_must_be_integer", + { + "n_exp": 2, + "init_method": "latin_hypercube_sampling", + "init_n_samples": 2, + "init_seed": 1.5, + }, + ValueError, + r"``init_seed`` must be None or an integer", + ), + ] + ) + def test_optimize_experiments_init_argument_validation_cases( + self, _label, kwargs, exc_type, regex + ): + # These argument checks all fail before any model build, so this + # parameterized test keeps contracts aligned without duplicated setup. + doe_obj = self._make_dummy_optimize_experiments_doe() + with self.assertRaisesRegex(exc_type, regex): + doe_obj.optimize_experiments(**kwargs) + + def test_optimize_experiments_lhs_requires_template_mode(self): + # Tests that LHS initialization is disallowed in user-initialized + # multi-experiment mode. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment(), _DummyExperiment()], + objective_option="pseudo_trace", + ) + with self.assertRaisesRegex( + ValueError, + r"``init_method='latin_hypercube_sampling'`` is currently supported only in template mode", + ): + doe_obj.optimize_experiments(init_method="latin_hypercube_sampling") + + def test_optimize_experiments_lhs_requires_scipy(self): + # Tests that LHS initialization requires scipy to be available. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + old_scipy_available = doe_module.scipy_available + doe_module.scipy_available = False + try: + with self.assertRaisesRegex( + ImportError, r"LHS initialization requires scipy" + ): + doe_obj.optimize_experiments(init_method="latin_hypercube_sampling") + finally: + doe_module.scipy_available = old_scipy_available + + @unittest.skipIf( + not parameterized_available, "The 'parameterized' package is not available" + ) + @parameterized.parameterized.expand( + [ + ( + "n_exp_disallowed_with_multi_experiment_list", + 2, + {"n_exp": 2}, + ValueError, + r"``n_exp`` must not be set when the experiment list contains more " + "than one experiment", + ), + ( + "n_exp_must_be_positive", + 1, + {"n_exp": 0}, + ValueError, + r"``n_exp`` must be a positive integer, got 0.", + ), + ( + "results_file_must_be_str_or_path", + 1, + {"results_file": 5}, + ValueError, + r"``results_file`` must be either a Path object or a string.", + ), + ( + "init_solver_must_have_solve", + 1, + {"init_solver": object()}, + ValueError, + r"``init_solver`` must be None or a solver object with a 'solve' " + "method.", + ), + ] + ) + def test_optimize_experiments_general_argument_validation_cases( + self, _label, n_experiments, kwargs, exc_type, regex + ): + # These are remaining lightweight API validations whose only contract + # is the raised error for a bad user-facing kwarg/value pair. + doe_obj = self._make_dummy_optimize_experiments_doe(n_experiments=n_experiments) + with self.assertRaisesRegex(exc_type, regex): + doe_obj.optimize_experiments(**kwargs) + + +@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not scipy_available, "scipy is not available") +@unittest.skipIf(not pandas_available, "pandas is not available") +class TestDoEErrorsRequiringSolver(unittest.TestCase): + def _make_solver(self): + return make_ipopt_solver() + + @unittest.skipIf( + not parameterized_available, "The 'parameterized' package is not available" + ) + @parameterized.parameterized.expand( + [("minimum_eigenvalue",), ("condition_number",)] + ) + def test_optimize_experiments_non_greybox_rejects_e_and_me_objectives( + self, objective_option + ): + # E-opt and ME-opt require the greybox objective path in + # optimize_experiments(); the standard algebraic multi-experiment build + # should fail fast with a user-facing validation error instead. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option=objective_option, + step=1e-2, + solver=self._make_solver(), + ) + + with self.assertRaisesRegex( + ValueError, + rf"objective_option='{objective_option}' requires " + r"use_grey_box_objective=True\.", + ): + doe_obj.optimize_experiments(n_exp=2) + + def test_optimize_experiments_trace_requires_cholesky_or_greybox(self): + # Multi-experiment trace uses the Cholesky-based build unless the + # greybox objective path is enabled, so optimize_experiments() should + # reject ``_Cholesky_option=False`` before any solve phase begins. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option="trace", + step=1e-2, + solver=self._make_solver(), + _Cholesky_option=False, + ) + + with self.assertRaisesRegex( + ValueError, + r"objective_option='trace' currently only implemented with " + r"``_Cholesky_option=True`` or ``use_grey_box_objective=True``\.", + ): + doe_obj.optimize_experiments(n_exp=2) + + def test_optimize_experiments_requires_matching_unknown_parameter_values(self): + # Tests that user-initialized multi-experiment mode rejects experiments + # that linearize around different nominal theta values. + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiExperiment( + hour=2.0, y=10.0, theta={'asymptote': 15.0, 'rate_constant': 0.5} + ), + RooneyBieglerMultiExperiment( + hour=3.0, y=11.0, theta={'asymptote': 15.0, 'rate_constant': 0.6} + ), + ], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + + with self.assertRaisesRegex( + ValueError, "must use the same nominal values for unknown_parameters" + ): + doe_obj.optimize_experiments() + + def test_optimize_experiments_requires_matching_unknown_parameter_labels(self): + # Tests that user-initialized multi-experiment mode rejects experiments + # whose unknown-parameter sets differ. + class _DifferentUnknownParameterExperiment: + def __init__(self, base_exp): + self._base_exp = base_exp + + def get_labeled_model(self, **kwargs): + m = self._base_exp.get_labeled_model(**kwargs) + m.fake_theta = pyo.Var(initialize=1.0) + m.fake_theta.fix() + m.del_component(m.unknown_parameters) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update( + [ + (m.asymptote, pyo.value(m.asymptote)), + (m.fake_theta, pyo.value(m.fake_theta)), + ] + ) + return m + + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiExperiment(hour=2.0, y=10.0), + _DifferentUnknownParameterExperiment( + RooneyBieglerMultiExperiment(hour=3.0, y=11.0) + ), + ], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + + with self.assertRaisesRegex( + ValueError, "must define the same unknown_parameters in the same order" + ): + doe_obj.optimize_experiments() + + def test_optimize_experiments_requires_matching_unknown_parameter_order(self): + # Tests that user-initialized multi-experiment mode rejects experiments + # whose unknown parameters appear in a different order. + class _ReorderedUnknownParameterExperiment: + def __init__(self, base_exp): + self._base_exp = base_exp + + def get_labeled_model(self, **kwargs): + m = self._base_exp.get_labeled_model(**kwargs) + m.del_component(m.unknown_parameters) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update( + [ + (m.rate_constant, pyo.value(m.rate_constant)), + (m.asymptote, pyo.value(m.asymptote)), + ] + ) + return m + + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiExperiment(hour=2.0, y=10.0), + _ReorderedUnknownParameterExperiment( + RooneyBieglerMultiExperiment(hour=3.0, y=11.0) + ), + ], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + + with self.assertRaisesRegex( + ValueError, "must define the same unknown_parameters in the same order" + ): + doe_obj.optimize_experiments() + + def test_optimize_experiments_sym_break_var_must_be_input(self): + # Tests that symmetry-breaking marker variables must also be experiment inputs. + class _BadSymBreakExperiment: + def __init__(self, base_exp): + self._base_exp = base_exp + + def get_labeled_model(self, **kwargs): + m = self._base_exp.get_labeled_model(**kwargs) + m.sym_break_cons = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.sym_break_cons[next(iter(m.unknown_parameters.keys()))] = None + return m + + exp = _BadSymBreakExperiment(RooneyBieglerMultiExperiment(hour=2.0, y=10.0)) + doe_obj = DesignOfExperiments( + experiment=[exp], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + + with self.assertRaisesRegex( + ValueError, "sym_break_cons.*must also be an experiment input variable" + ): + doe_obj.optimize_experiments(n_exp=2) + + def test_optimize_experiments_symmetry_mapping_failure_raises(self): + # Tests that failure to map symmetry variable across experiment blocks raises. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + ) + probe_model = doe_obj.experiment_list[0].get_labeled_model( + **doe_obj.get_labeled_model_args + ) + sym_var_name = next(iter(probe_model.experiment_inputs.keys())).local_name + original_find = pyo.ComponentUID.find_component_on + + def _fail_only_symmetry_mapping(cuid, block): + if ( + sym_var_name in str(cuid) + and hasattr(block, "experiment_inputs") + and block.index() == 0 + ): + return None + return original_find(cuid, block) + + pyo.ComponentUID.find_component_on = _fail_only_symmetry_mapping + try: + with self.assertRaisesRegex( + RuntimeError, "Failed to map symmetry breaking variable" + ): + doe_obj.optimize_experiments(n_exp=2) + finally: + pyo.ComponentUID.find_component_on = original_find + + def test_optimize_experiments_symmetry_breaking_default_variable_warning(self): + # Tests that missing explicit symmetry marker triggers warning and default + # choice. + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiInputExperimentFlag(hour=2.0, sym_break_flag=0), + RooneyBieglerMultiInputExperimentFlag(hour=4.0, sym_break_flag=0), + ], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as cm: + doe_obj.optimize_experiments() + self.assertTrue( + any("No symmetry breaking variable specified" in msg for msg in cm.output) + ) + self.assertTrue( + hasattr(doe_obj.model.param_scenario_blocks[0], "symmetry_breaking_s0_exp1") + ) + + def test_optimize_experiments_symmetry_breaking_multiple_markers_warning(self): + # Tests that multiple symmetry markers trigger an ambiguity warning. + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiInputExperimentFlag(hour=2.0, sym_break_flag=2), + RooneyBieglerMultiInputExperimentFlag(hour=4.0, sym_break_flag=2), + ], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as cm: + doe_obj.optimize_experiments() + self.assertTrue( + any( + "Multiple variables marked in sym_break_cons" in msg + for msg in cm.output + ) + ) + + def test_lhs_initialization_large_space_emits_warnings(self): + # Tests that very large LHS candidate/combo spaces emit user-facing warnings. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as log_cm: + with warnings.catch_warnings(record=True) as warn_cm: + warnings.simplefilter("always") + # Keep this warning-path test fast and deterministic: with + # init_n_samples=10001 and n_exp=2, the real path would do + # 10,001 candidate FIM evaluations and score + # C(10,001, 2)=50,005,000 combinations, which is too expensive + # for a unit test. We only need to verify warning contracts. + original_combinations = doe_module.combinations + original_compute_fim = doe_obj._compute_fim_at_point_no_prior + doe_module.combinations = lambda *_args, **_kwargs: iter([(0, 1)]) + doe_obj._compute_fim_at_point_no_prior = lambda *args, **kwargs: np.eye( + 2 + ) + try: + doe_obj.optimize_experiments( + n_exp=2, + init_method="latin_hypercube_sampling", + init_n_samples=10001, + # Set random seed to keep LHS initialization deterministic. + init_seed=11, + ) + finally: + doe_module.combinations = original_combinations + doe_obj._compute_fim_at_point_no_prior = original_compute_fim + + self.assertTrue( + any("candidate experiment designs" in str(w.message) for w in warn_cm) + ) + self.assertTrue(any("combinations to evaluate" in msg for msg in log_cm.output)) + + def test_lhs_missing_bounds_error_message(self): + # Tests that LHS initialization fails fast when experiment inputs lack bounds. + doe_obj = DesignOfExperiments( + experiment=[ + RooneyBieglerMultiExperiment(hour=2.0, hour_bounds=(None, 10.0)) + ], + objective_option="pseudo_trace", + ) + with self.assertRaisesRegex( + ValueError, + r"LHS initialization requires explicit lower and upper bounds on all " + "experiment input variables", + ): + doe_obj.optimize_experiments( + init_method="latin_hypercube_sampling", init_n_samples=2 + ) + def test_invalid_determinant_without_cholesky(self): fd_method = "central" obj_used = "determinant" diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py index 4329ddf37ab..e155a7e52ac 100644 --- a/pyomo/contrib/doe/tests/test_doe_solve.py +++ b/pyomo/contrib/doe/tests/test_doe_solve.py @@ -9,6 +9,7 @@ import json import logging import os, os.path +from itertools import combinations, product from glob import glob from pyomo.common.dependencies import ( @@ -39,10 +40,12 @@ ) from pyomo.contrib.doe.tests.experiment_class_example_flags import ( RooneyBieglerExperimentBad, + RooneyBieglerMultiExperiment, ) from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( RooneyBieglerExperiment, ) +from pyomo.contrib.doe.tests.utils_for_doe_tests import make_ipopt_solver from pyomo.contrib.doe.utils import rescale_FIM from pyomo.contrib.doe.examples.rooney_biegler_doe_example import run_rooney_biegler_doe @@ -83,6 +86,11 @@ def get_rooney_biegler_experiment(): ) +def _optimize_experiments_param_scenario(results, index=0): + """Return one parameter-scenario entry from optimize_experiments() results.""" + return results["solution"]["param_scenarios"][index] + + def get_FIM_Q_L(doe_obj=None): """ Helper function to retrieve results to compare. @@ -112,7 +120,7 @@ def get_FIM_Q_L(doe_obj=None): for j in model.parameter_names ] sigma_inv = [ - 1 / v**2 for k, v in model.scenario_blocks[0].measurement_error.items() + 1 / v**2 for k, v in model.fd_scenario_blocks[0].measurement_error.items() ] FIM_vals_np = np.array(FIM_vals).reshape((n_param, n_param)) @@ -134,7 +142,7 @@ def get_FIM_Q_L(doe_obj=None): def get_standard_args(experiment, fd_method, obj_used): args = {} - args['experiment'] = experiment + args['experiment'] = None if experiment is None else [experiment] args['fd_formula'] = fd_method args['step'] = 1e-3 args['objective_option'] = obj_used @@ -145,10 +153,7 @@ def get_standard_args(experiment, fd_method, obj_used): args['fim_initial'] = None args['L_diagonal_lower_bound'] = 1e-7 # Make solver object with good linear subroutines - solver = SolverFactory("ipopt") - solver.options["linear_solver"] = "ma57" - solver.options["halt_on_ampl_error"] = "yes" - solver.options["max_iter"] = 3000 + solver = make_ipopt_solver() args['solver'] = solver args['tee'] = False args['get_labeled_model_args'] = None @@ -378,12 +383,43 @@ def test_compute_FIM_seq_forward(self): doe_obj.compute_FIM(method="sequential") + def test_compute_FIM_multi_experiment_is_sum_of_experiment_fims(self): + fd_method = "central" + obj_used = "pseudo_trace" + + prior_fim = np.array([[1.5, 0.1], [0.1, 0.5]]) + + multi_args = get_standard_args( + RooneyBieglerMultiExperiment(hour=1.5, y=9.0), fd_method, obj_used + ) + multi_args["experiment"] = [ + RooneyBieglerMultiExperiment(hour=1.5, y=9.0), + RooneyBieglerMultiExperiment(hour=3.5, y=12.0), + ] + multi_args["prior_FIM"] = prior_fim + doe_multi = DesignOfExperiments(**multi_args) + + fim_total = doe_multi.compute_FIM(method="sequential") + + self.assertEqual(len(doe_multi._computed_FIM_by_experiment), 2) + fim_expected = prior_fim.copy() + for fim_exp in doe_multi._computed_FIM_by_experiment: + fim_expected = fim_expected + np.array(fim_exp) + + self.assertTrue(np.allclose(fim_total, fim_expected, atol=1e-8)) + self.assertTrue( + np.all(np.isfinite(np.array(doe_multi._computed_FIM_by_experiment[0]))) + ) + self.assertTrue( + np.all(np.isfinite(np.array(doe_multi._computed_FIM_by_experiment[1]))) + ) + # This test ensure that compute FIM runs without error using the # `kaug` option. kaug computes the FIM directly so no finite difference # scheme is needed. @unittest.skipIf(not scipy_available, "Scipy is not available") @unittest.skipIf( - not k_aug_available.available(False), "The 'k_aug' command is not available" + not k_aug_available, "The 'k_aug' is not available in this environment" ) def test_compute_FIM_kaug(self): fd_method = "forward" @@ -481,7 +517,7 @@ def test_rescale_FIM(self): [ [ v - for k, v in doe_obj.model.scenario_blocks[ + for k, v in doe_obj.model.fd_scenario_blocks[ 0 ].unknown_parameters.items() ] @@ -917,5 +953,621 @@ def cleanup_files(): ) +@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not scipy_available, "scipy is not available") +class TestOptimizeExperimentsAlgorithm(unittest.TestCase): + def _make_template_doe(self, objective_option="pseudo_trace"): + exp = RooneyBieglerMultiExperiment(hour=2.0, y=10.0) + solver = make_ipopt_solver() + return DesignOfExperiments( + experiment=[exp], + objective_option=objective_option, + step=1e-2, + solver=solver, + ) + + def _build_template_model_for_multi_experiment(self, doe_obj, n_exp): + doe_obj.model.param_scenario_blocks = pyo.Block(range(1)) + doe_obj.model.param_scenario_blocks[0].exp_blocks = pyo.Block(range(n_exp)) + for k in range(n_exp): + doe_obj.create_doe_model( + model=doe_obj.model.param_scenario_blocks[0].exp_blocks[k], + experiment_index=0, + _for_multi_experiment=True, + ) + + @staticmethod + def _make_square_solve_then_stub_solver(real_solver): + def _stub_results(): + class _SolverInfo: + status = "ok" + termination_condition = "optimal" + message = "stub-solve" + + class _SolverResults: + solver = _SolverInfo() + + return _SolverResults() + + class _SquareSolveThenStubSolver: + def __init__(self, wrapped_solver): + self._wrapped_solver = wrapped_solver + self.name = getattr(wrapped_solver, "name", "stub-solver") + self.options = getattr(wrapped_solver, "options", {}) + + def solve(self, *args, **kwargs): + model = args[0] if args else kwargs.get("model", None) + if model is not None and hasattr(model, "dummy_obj"): + # Keep square-solve path real so model state initializes correctly. + return self._wrapped_solver.solve(*args, **kwargs) + return _stub_results() + + return _SquareSolveThenStubSolver(real_solver) + + def test_evaluate_objective_from_fim_numerical_values(self): + fim = np.array([[4.0, 1.0], [1.0, 3.0]]) + expected_det = 11.0 + expected_pseudo_trace = 7.0 + expected_trace = 7.0 / 11.0 + + doe_det = self._make_template_doe("determinant") + self.assertAlmostEqual( + doe_det._evaluate_objective_from_fim(fim), expected_det, places=10 + ) + + doe_ptr = self._make_template_doe("pseudo_trace") + self.assertAlmostEqual( + doe_ptr._evaluate_objective_from_fim(fim), expected_pseudo_trace, places=10 + ) + + doe_tr = self._make_template_doe("trace") + self.assertAlmostEqual( + doe_tr._evaluate_objective_from_fim(fim), expected_trace, places=10 + ) + + def test_evaluate_objective_from_fim_fallback_paths(self): + singular_fim = np.zeros((2, 2)) + + doe_trace = self._make_template_doe("trace") + self.assertEqual(doe_trace._evaluate_objective_from_fim(singular_fim), np.inf) + + doe_zero = self._make_template_doe("zero") + self.assertEqual(doe_zero._evaluate_objective_from_fim(singular_fim), 0.0) + + def test_compute_fim_at_point_no_prior_restores_prior(self): + doe_no_prior = self._make_template_doe("pseudo_trace") + doe_no_prior.prior_FIM = np.zeros((2, 2)) + expected = doe_no_prior.compute_FIM(method="sequential") + + doe = self._make_template_doe("pseudo_trace") + saved_prior = np.array([[7.0, 0.0], [0.0, 5.0]]) + doe.prior_FIM = saved_prior + got = doe._compute_fim_at_point_no_prior(experiment_index=0, input_values=[2.0]) + + self.assertTrue(np.allclose(got, expected, atol=1e-8)) + self.assertIs(doe.prior_FIM, saved_prior) + self.assertTrue(np.allclose(doe.prior_FIM, saved_prior, atol=1e-12)) + + def test_compute_fim_at_point_no_prior_exception_fallback_zero(self): + doe = self._make_template_doe("pseudo_trace") + saved_prior = np.array([[3.0, 0.0], [0.0, 4.0]]) + doe.prior_FIM = saved_prior + + original_sequential_fim = doe._sequential_FIM + + def _raise_boom(*args, **kwargs): + # Intentionally inject a sequential FIM failure so we can verify the + # fallback behavior (zero FIM + warning + prior restoration). + raise RuntimeError("boom") + + doe._sequential_FIM = _raise_boom + try: + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as log_cm: + got = doe._compute_fim_at_point_no_prior( + experiment_index=0, input_values=[2.0] + ) + finally: + doe._sequential_FIM = original_sequential_fim + + self.assertTrue(np.allclose(got, np.zeros((2, 2)))) + self.assertIs(doe.prior_FIM, saved_prior) + self.assertTrue( + any("Using zero matrix as fallback" in msg for msg in log_cm.output) + ) + + def test_optimize_experiments_cholesky_jitter_branch(self): + # Force the positive-definiteness check down the jitter path and verify + # the matrix passed into Cholesky includes the expected diagonal shift. + doe = self._make_template_doe("determinant") + # Keep the square initialization solve real (path under test), but + # stub the final solve. Without this, this test setup can fail + # before assertions with solver status=error on the final solve because + # rooney-biegler's model will not generate a negative eigval. + doe.solver = self._make_square_solve_then_stub_solver(doe.solver) + + from pyomo.contrib.doe.doe import _SMALL_TOLERANCE_DEFINITENESS + + doe.optimize_experiments(n_exp=1) + + scenario = doe.model.param_scenario_blocks[0] + total_fim = np.array( + _optimize_experiments_param_scenario(doe.results)["total_fim"] + ) + min_eig = np.min(np.real(np.linalg.eigvals(total_fim))) + if min_eig < _SMALL_TOLERANCE_DEFINITENESS: + jitter = min( + -min_eig + _SMALL_TOLERANCE_DEFINITENESS, _SMALL_TOLERANCE_DEFINITENESS + ) + else: + jitter = 0.0 + expected_cholesky_input = total_fim + jitter * np.eye(total_fim.shape[0]) + param_names = list(scenario.exp_blocks[0].parameter_names) + L_vals = np.array( + [ + [pyo.value(scenario.obj_cons.L[p, q]) for q in param_names] + for p in param_names + ] + ) + + self.assertEqual(doe.results["optimization_solve"]["status"], "ok") + self.assertGreater(jitter, 0.0) + self.assertTrue( + np.allclose(L_vals @ L_vals.T, expected_cholesky_input, atol=1e-8) + ) + + def test_lhs_initialize_experiments_matches_bruteforce_combo(self): + # Compare the helper's chosen initial design directly against an + # independent brute-force scorer so this test isolates LHS selection + # instead of the later NLP solve and result-payload bookkeeping. + n_exp = 2 + lhs_n_samples = 3 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 17 + + # Build expected best combo using the same helper path and objective scoring. + expected_obj = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(expected_obj, n_exp=n_exp) + + first_exp_block = expected_obj.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = expected_obj._get_experiment_input_vars(first_exp_block) + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + + rng = np.random.default_rng(lhs_seed) + from scipy.stats.qmc import LatinHypercube + + per_dim_samples = [] + for i in range(len(exp_input_vars)): + dim_seed = int(rng.integers(0, 2**31)) + sampler = LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + + candidate_points = list(product(*per_dim_samples)) + candidate_fims = [ + expected_obj._compute_fim_at_point_no_prior(0, list(pt)) + for pt in candidate_points + ] + + best_combo = None + best_obj = -np.inf + for combo in combinations(range(len(candidate_points)), n_exp): + fim_total = sum((candidate_fims[idx] for idx in combo), np.zeros((2, 2))) + obj_val = expected_obj._evaluate_objective_from_fim(fim_total) + if obj_val > best_obj: + best_obj = obj_val + best_combo = combo + + expected_points = [list(candidate_points[i]) for i in best_combo] + + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=n_exp) + actual_points, lhs_diag = doe._lhs_initialize_experiments( + lhs_n_samples=lhs_n_samples, lhs_seed=lhs_seed, n_exp=n_exp + ) + + actual_points_norm = sorted(tuple(np.round(p, 8)) for p in actual_points) + expected_points_norm = sorted(tuple(np.round(p, 8)) for p in expected_points) + self.assertEqual(actual_points_norm, expected_points_norm) + self.assertAlmostEqual(lhs_diag["best_obj"], best_obj, places=12) + + def test_optimize_experiments_is_reentrant_on_same_object(self): + # Tests that optimize_experiments() can be run repeatedly on one DoE object. + doe = self._make_template_doe("pseudo_trace") + doe.optimize_experiments(n_exp=1) + first_design = _optimize_experiments_param_scenario(doe.results)["experiments"][ + 0 + ]["design"] + first_build_time = doe.results["timing"]["build_time_s"] + + doe.optimize_experiments(n_exp=1) + second_design = _optimize_experiments_param_scenario(doe.results)[ + "experiments" + ][0]["design"] + + self.assertEqual(len(first_design), len(second_design)) + self.assertIn("timing", doe.results) + self.assertGreater(doe.results["timing"]["build_time_s"], 0.0) + self.assertGreaterEqual(first_build_time, 0.0) + self.assertEqual(len(list(doe.model.param_scenario_blocks.keys())), 1) + + def test_lhs_no_candidate_fim_scored_uses_zero_fallback(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=1) + + # Force candidate FIM evaluations to return None so we exercise the + # LHS fallback path when no candidate FIM can be scored. + doe._compute_fim_at_point_no_prior = lambda experiment_index, input_values: None + # Set random seed to keep LHS initialization deterministic. + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=3, lhs_seed=8, n_exp=1 + ) + + self.assertEqual(len(points), 1) + self.assertFalse(diag["timed_out"]) + + def test_lhs_combo_serial_minimize_objective_update(self): + # Rooney-Biegler integration check for the LHS combination scorer on a + # minimize objective (trace/A-opt). Verify that optimize_experiments() + # picks the same initial-design pair as an exhaustive reference over + # the generated candidates. + n_exp = 2 + lhs_n_samples = 4 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 61 + + doe_ref = self._make_template_doe("trace") + self._build_template_model_for_multi_experiment(doe_ref, n_exp=n_exp) + + first_exp_block = doe_ref.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe_ref._get_experiment_input_vars(first_exp_block) + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + + rng = np.random.default_rng(lhs_seed) + from scipy.stats.qmc import LatinHypercube + + per_dim_samples = [] + for i in range(len(exp_input_vars)): + dim_seed = int(rng.integers(0, 2**31)) + sampler = LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + candidate_points = list(product(*per_dim_samples)) + candidate_fims = [ + doe_ref._compute_fim_at_point_no_prior(0, list(pt)) + for pt in candidate_points + ] + + best_combo = None + best_obj = np.inf + for combo in combinations(range(len(candidate_points)), n_exp): + fim_total = sum((candidate_fims[idx] for idx in combo), np.zeros((2, 2))) + obj_val = doe_ref._evaluate_objective_from_fim(fim_total) + if obj_val < best_obj: + best_obj = obj_val + best_combo = combo + expected_points = [list(candidate_points[i]) for i in best_combo] + + doe = self._make_template_doe("trace") + doe.optimize_experiments( + n_exp=n_exp, + init_method="latin_hypercube_sampling", + init_n_samples=lhs_n_samples, + init_seed=lhs_seed, + ) + got_points = doe.results["initialization"]["selected_initial_designs"] + + # Normalize (round + sort) so we compare selected points robustly + # despite floating-point noise and ordering differences. + got_norm = sorted(tuple(np.round(p, 8)) for p in got_points) + exp_norm = sorted(tuple(np.round(p, 8)) for p in expected_points) + self.assertEqual(got_norm, exp_norm) + + def test_lhs_combo_no_scored_combo_falls_back_to_first_n_exp(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=3) + lhs_n_samples = 3 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 13 + first_exp_block = doe.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe._get_experiment_input_vars(first_exp_block) + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + rng = np.random.default_rng(lhs_seed) + from scipy.stats.qmc import LatinHypercube + + per_dim_samples = [] + for i in range(len(exp_input_vars)): + dim_seed = int(rng.integers(0, 2**31)) + sampler = LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + candidate_points = list(product(*per_dim_samples)) + expected_points = [list(candidate_points[i]) for i in range(3)] + + def _nan_fim(experiment_index, input_values): + # Any comparison with NaN is False, so best_combo remains None and + # we exercise the fallback branch that selects the first n_exp points. + return np.array([[np.nan, 0.0], [0.0, np.nan]]) + + doe._compute_fim_at_point_no_prior = _nan_fim + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as log_cm: + got_points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=lhs_n_samples, lhs_seed=lhs_seed, n_exp=3 + ) + + # Normalize (round + sort) so we compare selected points robustly + # despite floating-point noise and ordering differences. + got_norm = sorted(tuple(np.round(p, 8)) for p in got_points) + exp_norm = sorted(tuple(np.round(p, 8)) for p in expected_points) + self.assertEqual(got_norm, exp_norm) + # The fallback branch must be user-visible: no scored combinations + # should emit the warning that first n_exp candidates were used. + self.assertTrue( + any( + "Falling back to the first n_exp candidate points." in msg + for msg in log_cm.output + ) + ) + + def test_lhs_combo_scoring_n_exp_3_matches_exhaustive_reference(self): + n_exp = 3 + doe_ref = self._make_template_doe("pseudo_trace") + # Keep n_exp=3 to exercise the general combination-scoring path + # (different from the optimized n_exp==2 branch). + self._build_template_model_for_multi_experiment(doe_ref, n_exp=n_exp) + lhs_n_samples = 5 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 2 + + # Recreate the exact candidate points from LHS generation and compute + # an independent exhaustive reference for combination scoring. + first_exp_block = doe_ref.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe_ref._get_experiment_input_vars(first_exp_block) + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + rng = np.random.default_rng(lhs_seed) + from scipy.stats.qmc import LatinHypercube + + per_dim_samples = [] + for i in range(len(exp_input_vars)): + dim_seed = int(rng.integers(0, 2**31)) + sampler = LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + candidate_points = list(product(*per_dim_samples)) + + # Exhaustive reference over all combinations of size 3. + fims = [ + doe_ref._compute_fim_at_point_no_prior(0, list(pt)) + for pt in candidate_points + ] + best_obj = -np.inf + best_combo = None + for combo in combinations(range(len(candidate_points)), n_exp): + fim_total = sum((fims[i] for i in combo), np.zeros((2, 2))) + obj = doe_ref._evaluate_objective_from_fim(fim_total) + if obj > best_obj: + best_obj = obj + best_combo = combo + expected_points = [list(candidate_points[i]) for i in best_combo] + + doe = self._make_template_doe("pseudo_trace") + doe.optimize_experiments( + n_exp=n_exp, + init_method="latin_hypercube_sampling", + init_n_samples=lhs_n_samples, + init_seed=lhs_seed, + ) + got_points = doe.results["initialization"]["selected_initial_designs"] + + # Normalize (round + sort) so we compare selected points robustly + # despite floating-point noise and ordering differences. + got_norm = sorted(tuple(np.round(p, 8)) for p in got_points) + exp_norm = sorted(tuple(np.round(p, 8)) for p in expected_points) + self.assertEqual(got_norm, exp_norm) + + def test_lhs_matches_exhaustive_reference_with_fixed_samples(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + lhs_n_samples = 3 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 123 + first_exp_block = doe.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe._get_experiment_input_vars(first_exp_block) + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + rng = np.random.default_rng(lhs_seed) + from scipy.stats.qmc import LatinHypercube + + per_dim_samples = [] + for i in range(len(exp_input_vars)): + dim_seed = int(rng.integers(0, 2**31)) + sampler = LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + candidate_points = list(product(*per_dim_samples)) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, 2.0 * x + 1.0]]) + + doe._compute_fim_at_point_no_prior = _fake_fim + got_points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=lhs_n_samples, lhs_seed=lhs_seed, n_exp=2 + ) + + # Independent exhaustive reference over explicit candidate points. + best_obj = -np.inf + best_combo = None + for combo in combinations(range(len(candidate_points)), 2): + f1 = _fake_fim(0, list(candidate_points[combo[0]])) + f2 = _fake_fim(0, list(candidate_points[combo[1]])) + obj_val = float(np.trace(f1 + f2)) + if obj_val > best_obj: + best_obj = obj_val + best_combo = combo + expected_points = [list(candidate_points[i]) for i in best_combo] + + # Normalize (round + sort) so we compare selected points robustly + # despite floating-point noise and ordering differences. + got_norm = sorted(tuple(np.round(p, 8)) for p in got_points) + exp_norm = sorted(tuple(np.round(p, 8)) for p in expected_points) + self.assertEqual(got_norm, exp_norm) + + def test_optimize_experiments_determinant_expected_values(self): + # Tests determinant-objective optimization against known expected design/metric values. + # Match the multi-experiment example style (explicit experiment list) + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.0, y=8.3), + RooneyBieglerMultiExperiment(hour=2.0, y=10.3), + ] + solver = make_ipopt_solver() + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="determinant", + step=1e-2, + solver=solver, + ) + doe.optimize_experiments() + + scenario = _optimize_experiments_param_scenario(doe.results) + got_hours = sorted(exp["design"][0] for exp in scenario["experiments"]) + expected_hours = [1.9321985035514362, 9.999999685577139] + + self.assertStructuredAlmostEqual(got_hours, expected_hours, abstol=1e-3) + self.assertAlmostEqual( + scenario["quality_metrics"]["log10_d_opt"], 6.028152580313302, places=3 + ) + + def test_optimize_experiments_trace_expected_values(self): + # Tests trace-objective optimization against known expected design/metric values. + # Match the multi-experiment example style (explicit experiment list) + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.0, y=8.3), + RooneyBieglerMultiExperiment(hour=2.0, y=10.3), + ] + solver = make_ipopt_solver() + # prior_FIM from data `hour = 1, y = 8.3` with default value of parameters, which + # is theta = {'asymptote': 15, 'rate_constant': 0.5} + prior_FIM = np.array( + [[15.48181217, 357.97684273], [357.97684273, 8277.28811613]] + ) + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="trace", + step=1e-2, + solver=solver, + prior_FIM=prior_FIM, + ) + doe.optimize_experiments() + + scenario = _optimize_experiments_param_scenario(doe.results) + got_hours = sorted(exp["design"][0] for exp in scenario["experiments"]) + expected_hours = [10.0, 10.0] + + self.assertStructuredAlmostEqual(got_hours, expected_hours, abstol=1e-3) + self.assertAlmostEqual( + scenario["quality_metrics"]["log10_a_opt"], -2.2347, places=3 + ) + + def test_optimize_experiments_prior_fim_aggregation_non_lhs_template_mode(self): + # Tests that total FIM equals sum(experiment FIMs) + prior in template mode. + prior_fim = np.array([[2.0, 0.1], [0.1, 1.5]]) + doe = self._make_template_doe("pseudo_trace") + doe.prior_FIM = prior_fim.copy() + + doe.optimize_experiments(n_exp=2, init_method=None) + + scenario = _optimize_experiments_param_scenario(doe.results) + total_fim = np.array(scenario["total_fim"]) + exp_fim_sum = sum( + (np.array(exp_data["fim"]) for exp_data in scenario["experiments"]), + np.zeros_like(total_fim), + ) + stored_prior = np.array(doe.results["problem"]["prior_fim"]) + + self.assertTrue(np.allclose(total_fim, exp_fim_sum + prior_fim, atol=1e-6)) + self.assertTrue(np.allclose(total_fim, exp_fim_sum + stored_prior, atol=1e-6)) + + def test_optimize_experiments_prior_fim_aggregation_non_lhs_user_initialized_mode( + self, + ): + # Tests that total FIM aggregation with prior is correct in user-initialized mode. + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.5, y=9.0), + RooneyBieglerMultiExperiment(hour=3.5, y=12.0), + ] + solver = make_ipopt_solver() + prior_fim = np.array([[1.25, 0.05], [0.05, 0.9]]) + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="pseudo_trace", + step=1e-2, + solver=solver, + prior_FIM=prior_fim, + ) + doe.optimize_experiments(init_method=None) + + scenario = _optimize_experiments_param_scenario(doe.results) + total_fim = np.array(scenario["total_fim"]) + exp_fim_sum = sum( + (np.array(exp_data["fim"]) for exp_data in scenario["experiments"]), + np.zeros_like(total_fim), + ) + stored_prior = np.array(doe.results["problem"]["prior_fim"]) + + self.assertTrue(np.allclose(total_fim, exp_fim_sum + prior_fim, atol=1e-6)) + self.assertTrue(np.allclose(total_fim, exp_fim_sum + stored_prior, atol=1e-6)) + + def test_optimize_experiments_safe_metric_failure_sets_nan(self): + # Tests that metric-computation failures are captured as NaN. + # For this template setup, the A-opt metric naturally becomes non-finite + # and the API should surface it as NaN. + doe = self._make_template_doe("pseudo_trace") + doe.optimize_experiments(n_exp=1) + + scenario = _optimize_experiments_param_scenario(doe.results) + self.assertTrue(np.isnan(scenario["quality_metrics"]["log10_a_opt"])) + + def test_optimize_experiments_non_cholesky_determinant_initialization(self): + # Tests determinant initialization correctness when Cholesky formulation is disabled. + exp = RooneyBieglerMultiExperiment(hour=2.0, y=10.0) + solver = make_ipopt_solver() + doe = DesignOfExperiments( + experiment=[exp], + objective_option="determinant", + step=1e-2, + solver=solver, + _Cholesky_option=False, + _only_compute_fim_lower=False, + ) + # Keep the square initialization solve real and stub the final solve. + # Without this, this case can terminate at max iterations and the + # determinant assertion becomes solver-fragile/non-deterministic. + doe.solver = self._make_square_solve_then_stub_solver(doe.solver) + doe.optimize_experiments(n_exp=1) + + scenario_block = doe.model.param_scenario_blocks[0] + self.assertTrue(hasattr(scenario_block.obj_cons, "determinant")) + total_fim = np.array( + _optimize_experiments_param_scenario(doe.results)["total_fim"] + ) + expected_det = np.linalg.det(total_fim) + self.assertAlmostEqual( + pyo.value(scenario_block.obj_cons.determinant), expected_det, places=6 + ) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/doe/tests/test_greybox.py b/pyomo/contrib/doe/tests/test_greybox.py index f1bcbbe9f64..7948ed30ac0 100644 --- a/pyomo/contrib/doe/tests/test_greybox.py +++ b/pyomo/contrib/doe/tests/test_greybox.py @@ -17,8 +17,11 @@ pandas as pd, pandas_available, scipy_available, + attempt_import, ) +parameterized, parameterized_available = attempt_import('parameterized') + from pyomo.common.fileutils import this_file_dir import pyomo.common.unittest as unittest @@ -30,9 +33,13 @@ from pyomo.contrib.doe.examples.reactor_example import ( ReactorExperiment as FullReactorExperiment, ) + from pyomo.contrib.doe.tests.experiment_class_example_flags import ( + RooneyBieglerMultiExperiment, + ) from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( RooneyBieglerExperiment, ) +from pyomo.contrib.doe.tests.utils_for_doe_tests import make_ipopt_solver import pyomo.environ as pyo @@ -288,15 +295,8 @@ def get_standard_args(experiment, fd_method, obj_used): args['jac_initial'] = None args['fim_initial'] = None args['L_diagonal_lower_bound'] = 1e-7 - solver = SolverFactory("ipopt") - args['solver'] = solver - # Make solver object with - # good linear subroutines - solver = SolverFactory("ipopt") - solver.options["linear_solver"] = "ma57" - solver.options["halt_on_ampl_error"] = "yes" - solver.options["max_iter"] = 3000 - args['solver'] = solver + # Make solver object with good linear subroutines. + args['solver'] = make_ipopt_solver() # Make greybox solver object with # good linear subroutines grey_box_solver = SolverFactory("cyipopt") @@ -396,6 +396,131 @@ def make_greybox_and_doe_objects_rooney_biegler(objective_option): return doe_obj, grey_box_object +class _MockGreyBoxSolver: + def __init__(self, name="mock-greybox"): + self.name = name + self.calls = [] + + def solve(self, model, tee=False): + self.calls.append({"model": model, "tee": tee}) + + class _MockSolverInfo: + status = "ok" + termination_condition = "optimal" + message = "mock-greybox-solve" + + class _MockResults: + solver = _MockSolverInfo() + + return _MockResults() + + +class _TrackingSolverWrapper: + """Proxy solver that records solve routing while delegating real solves.""" + + def __init__(self, solver, phase_tag, call_order, forbid=False): + self._solver = solver + self._phase_tag = phase_tag + self._call_order = call_order + self._forbid = forbid + self.solve_calls_count = 0 + + def solve(self, *args, **kwargs): + self.solve_calls_count += 1 + self._call_order.append(self._phase_tag) + if self._forbid: + raise AssertionError( + "Primary solver should not be used in greybox optimize_experiments()." + ) + return self._solver.solve(*args, **kwargs) + + def __getattr__(self, name): + return getattr(self._solver, name) + + +def _make_cyipopt_solver(tol=1e-4): + grey_box_solver = SolverFactory("cyipopt") + grey_box_solver.config.options["linear_solver"] = "ma57" + grey_box_solver.config.options['tol'] = tol + grey_box_solver.config.options['mu_strategy'] = "monotone" + return grey_box_solver + + +def _make_multiexperiment_greybox_doe( + objective_option, prior_FIM=None, grey_box_solver=None +): + if prior_FIM is None: + prior_FIM = np.zeros((2, 2)) + return DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option=objective_option, + use_grey_box_objective=True, + step=1e-2, + solver=make_ipopt_solver(), + grey_box_solver=( + grey_box_solver if grey_box_solver is not None else _MockGreyBoxSolver() + ), + prior_FIM=prior_FIM, + ) + + +def _get_multiexperiment_scenario_data(doe_obj): + scenario = doe_obj.model.param_scenario_blocks[0] + total_fim = np.array(doe_obj.results["solution"]["param_scenarios"][0]["total_fim"]) + parameter_names = list(scenario.exp_blocks[0].parameter_names) + return scenario, total_fim, parameter_names + + +def _generate_lhs_candidate_points(doe_obj, lhs_n_samples, lhs_seed): + from scipy.stats.qmc import LatinHypercube + + first_exp_block = doe_obj.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe_obj._get_experiment_input_vars(first_exp_block) + lb_vals = np.array([v.lb for v in exp_input_vars]) + ub_vals = np.array([v.ub for v in exp_input_vars]) + + # Set random seed to keep LHS candidate generation deterministic. + rng = np.random.default_rng(lhs_seed) + per_dim_samples = [] + for i in range(len(exp_input_vars)): + dim_seed = int(rng.integers(0, 2**31)) + sampler = LatinHypercube(d=1, seed=dim_seed) + s_unit = sampler.random(n=lhs_n_samples).flatten() + s_scaled = lb_vals[i] + s_unit * (ub_vals[i] - lb_vals[i]) + per_dim_samples.append(s_scaled.tolist()) + + return list(itertools.product(*per_dim_samples)) + + +def _expected_multiexperiment_greybox_output(objective_option, fim_np): + if objective_option == "trace": + return float(np.trace(np.linalg.pinv(fim_np))) + if objective_option == "determinant": + return float(np.log(np.linalg.det(fim_np))) + if objective_option == "minimum_eigenvalue": + return float(np.min(np.linalg.eigvalsh(fim_np))) + if objective_option == "condition_number": + eig = np.linalg.eigvalsh(fim_np) + return float(np.log(np.abs(np.max(eig) / np.min(eig)))) + raise AssertionError(f"Unexpected greybox objective: {objective_option!r}") + + +def _reconstruct_fim_from_egb_inputs(egb_block, parameter_names): + dim = len(parameter_names) + fim = np.zeros((dim, dim)) + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if i >= j: + fim[i, j] = pyo.value(egb_block.inputs[(q, p)]) + fim[j, i] = fim[i, j] + return fim + + +def _diagonal_hour_fim_reference(experiment_index, input_values): + hour = float(input_values[0]) + return np.eye(2) * (hour + 1.0) + + # Test whether or not cyipopt # is appropriately calling the # linear solvers. @@ -1066,7 +1191,8 @@ def test_constructor_doe_object_error(self): """Ensure constructor rejects missing DoE object dependencies.""" with self.assertRaisesRegex( ValueError, - "DoE Object must be provided to build external grey box of the FIM.", + "Either ``doe_object`` or both ``parameter_names`` and ``fim_initial`` " + "must be provided to build the FIM grey box.", ): FIMExternalGreyBox(doe_object=None) @@ -1374,5 +1500,601 @@ def test_solve_ME_optimality_condition_number(self): ) +@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not scipy_available, "scipy is not available") +@unittest.skipIf(not pandas_available, "pandas is not available") +@unittest.skipIf(not cyipopt_available, "'cyipopt' is not available") +class TestMultiexperimentBuild(unittest.TestCase): + def test_optimize_experiments_greybox_uses_aggregated_fim(self): + # Check that the multi-experiment greybox block is seeded from the + # aggregated scenario FIM and that the final solve is routed through + # the greybox solver interface. + doe_obj = _make_multiexperiment_greybox_doe( + objective_option="minimum_eigenvalue", + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + ) + + doe_obj.optimize_experiments(n_exp=2) + + scenario, total_fim, parameter_names = _get_multiexperiment_scenario_data( + doe_obj + ) + self.assertTrue(hasattr(scenario.obj_cons, "egb_fim_block")) + self.assertFalse(hasattr(scenario.obj_cons, "L")) + + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if i >= j: + self.assertAlmostEqual( + pyo.value(scenario.obj_cons.egb_fim_block.inputs[(q, p)]), + total_fim[i, j], + places=7, + ) + + self.assertAlmostEqual( + pyo.value(scenario.obj_cons.egb_fim_block.outputs["E-opt"]), + np.min(np.linalg.eigvalsh(total_fim)), + places=7, + ) + + @unittest.skipIf( + not parameterized_available, "The 'parameterized' package is not available" + ) + @parameterized.parameterized.expand( + [("determinant",), ("trace",), ("minimum_eigenvalue",), ("condition_number",)] + ) + def test_optimize_experiments_greybox_outputs_match_numpy_for_supported_objective( + self, objective_option + ): + # Validate the deterministic wiring using a real greybox solve: + # for each supported greybox metric, the external block output + # should match direct NumPy on scenario.total_fim. + prior_fim = np.array([[6.0, 0.75], [0.75, 4.0]]) + doe_obj = _make_multiexperiment_greybox_doe( + objective_option=objective_option, + prior_FIM=prior_fim.copy(), + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + ) + + doe_obj.optimize_experiments(n_exp=2) + + scenario, total_fim, parameter_names = _get_multiexperiment_scenario_data( + doe_obj + ) + egb_block = scenario.obj_cons.egb_fim_block + + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if i >= j: + self.assertAlmostEqual( + pyo.value(egb_block.inputs[(q, p)]), total_fim[i, j], places=7 + ) + + self.assertAlmostEqual( + pyo.value(egb_block.outputs[doe_obj._grey_box_output_name()]), + _expected_multiexperiment_greybox_output(objective_option, total_fim), + places=7, + ) + + def test_optimize_experiments_greybox_prior_fim_is_included_in_inputs_and_output( + self, + ): + # Use a large prior so the external block must see + # total_fim = sum(experiment_fim) + prior_FIM. + prior_fim = np.array([[100.0, 12.0], [12.0, 80.0]]) + doe_obj = _make_multiexperiment_greybox_doe( + objective_option="determinant", + prior_FIM=prior_fim.copy(), + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + ) + + doe_obj.optimize_experiments(n_exp=2) + + scenario, total_fim, parameter_names = _get_multiexperiment_scenario_data( + doe_obj + ) + exp_fim_sum = sum( + ( + np.array(exp_data["fim"]) + for exp_data in doe_obj.results["solution"]["param_scenarios"][0][ + "experiments" + ] + ), + np.zeros_like(total_fim), + ) + egb_block = scenario.obj_cons.egb_fim_block + + self.assertTrue(np.allclose(total_fim, exp_fim_sum + prior_fim, atol=1e-6)) + self.assertFalse(np.allclose(total_fim, exp_fim_sum, atol=1e-6)) + + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + if i >= j: + self.assertAlmostEqual( + pyo.value(egb_block.inputs[(q, p)]), total_fim[i, j], places=7 + ) + + output_with_prior = pyo.value(egb_block.outputs["log-D-opt"]) + self.assertAlmostEqual( + output_with_prior, + _expected_multiexperiment_greybox_output("determinant", total_fim), + places=7, + ) + self.assertFalse( + np.isclose( + output_with_prior, + _expected_multiexperiment_greybox_output("determinant", exp_fim_sum), + rtol=1e-6, + atol=1e-6, + ) + ) + + def test_optimize_experiments_greybox_uses_init_solver_for_square_solve_and_grey_box_solver_for_final_solve( + self, + ): + # Initialization should use init_solver; the greybox solver should be + # reserved for the final optimize_experiments NLP solve. + # Capture phase labels in solve-call order so we can assert that all + # initialization solves happen before the single final greybox solve. + solve_phase_order = [] + main_solver_inner = make_ipopt_solver() + init_solver_inner = make_ipopt_solver() + init_solver_inner.options["max_iter"] = 123 + grey_box_solver_inner = _make_cyipopt_solver(tol=1e-6) + main_solver = _TrackingSolverWrapper( + solver=main_solver_inner, + phase_tag="main", + call_order=solve_phase_order, + forbid=True, + ) + init_solver = _TrackingSolverWrapper( + solver=init_solver_inner, phase_tag="init", call_order=solve_phase_order + ) + grey_box_solver = _TrackingSolverWrapper( + solver=grey_box_solver_inner, + phase_tag="greybox", + call_order=solve_phase_order, + ) + + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="minimum_eigenvalue", + use_grey_box_objective=True, + step=1e-2, + solver=main_solver, + grey_box_solver=grey_box_solver, + ) + + doe_obj.optimize_experiments(n_exp=2, init_solver=init_solver) + + self.assertGreaterEqual(init_solver.solve_calls_count, 1) + self.assertEqual(grey_box_solver.solve_calls_count, 1) + self.assertEqual(solve_phase_order[-1], "greybox") + self.assertTrue(all(tag == "init" for tag in solve_phase_order[:-1])) + self.assertEqual( + doe_obj.results["initialization"]["solver"], + getattr(init_solver, "name", str(init_solver)), + ) + self.assertEqual( + doe_obj.results["optimization_solve"]["solver"], + getattr(grey_box_solver, "name", str(grey_box_solver)), + ) + + def test_optimize_experiments_greybox_is_reentrant_on_same_object(self): + # Re-running the same greybox DoE object should rebuild a fresh + # external block and reseed it from the current aggregated total FIM. + grey_box_solver_inner = _make_cyipopt_solver(tol=1e-6) + grey_box_solver = _TrackingSolverWrapper( + solver=grey_box_solver_inner, phase_tag="greybox", call_order=[] + ) + doe_obj = _make_multiexperiment_greybox_doe( + objective_option="minimum_eigenvalue", + prior_FIM=np.zeros((2, 2)), + grey_box_solver=grey_box_solver, + ) + + doe_obj.optimize_experiments(n_exp=2) + first_scenario, first_total_fim, _ = _get_multiexperiment_scenario_data(doe_obj) + first_egb_block = first_scenario.obj_cons.egb_fim_block + + self.assertAlmostEqual( + pyo.value(first_egb_block.outputs["E-opt"]), + _expected_multiexperiment_greybox_output( + "minimum_eigenvalue", first_total_fim + ), + places=7, + ) + + doe_obj.prior_FIM = np.array([[20.0, 2.0], [2.0, 15.0]]) + doe_obj.optimize_experiments(n_exp=2) + + second_scenario, second_total_fim, second_parameter_names = ( + _get_multiexperiment_scenario_data(doe_obj) + ) + second_egb_block = second_scenario.obj_cons.egb_fim_block + + self.assertIsNot(first_egb_block, second_egb_block) + self.assertEqual(len(list(doe_obj.model.param_scenario_blocks.keys())), 1) + self.assertEqual(grey_box_solver.solve_calls_count, 2) + self.assertFalse(np.allclose(first_total_fim, second_total_fim, atol=1e-6)) + + for i, p in enumerate(second_parameter_names): + for j, q in enumerate(second_parameter_names): + if i >= j: + self.assertAlmostEqual( + pyo.value(second_egb_block.inputs[(q, p)]), + second_total_fim[i, j], + places=7, + ) + + self.assertAlmostEqual( + pyo.value(second_egb_block.outputs["E-opt"]), + _expected_multiexperiment_greybox_output( + "minimum_eigenvalue", second_total_fim + ), + places=7, + ) + + @unittest.skipIf( + not parameterized_available, "The 'parameterized' package is not available" + ) + @parameterized.parameterized.expand( + [("minimum_eigenvalue",), ("condition_number",)] + ) + def test_optimize_experiments_greybox_lhs_initialization_scores_e_opt_and_me_opt( + self, objective_option + ): + # Rooney-Biegler integration check for the LHS candidate-combination + # scorer on greybox-only objectives (E-opt and ME-opt). Compare the + # selected initial designs against an independent exhaustive reference + # built from real candidate FIM evaluations. + lhs_n_samples = 4 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 19 + + doe_obj = _make_multiexperiment_greybox_doe( + objective_option=objective_option, + prior_FIM=np.zeros((2, 2)), + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + ) + + doe_obj.optimize_experiments( + n_exp=2, + init_method="latin_hypercube_sampling", + init_n_samples=lhs_n_samples, + init_seed=lhs_seed, + ) + + lhs_init = doe_obj.results["initialization"] + lhs_diag = { + "best_obj": lhs_init["best_initial_objective_value"], + "timed_out": lhs_init["timed_out"], + } + actual_points = lhs_init["selected_initial_designs"] + candidate_points = _generate_lhs_candidate_points( + doe_obj, lhs_n_samples=lhs_n_samples, lhs_seed=lhs_seed + ) + candidate_norm = {tuple(np.round(point, 8)) for point in candidate_points} + + self.assertEqual( + doe_obj.results["initialization"]["method"], "latin_hypercube_sampling" + ) + self.assertTrue(np.isfinite(lhs_diag["best_obj"])) + self.assertGreater(lhs_diag["best_obj"], 0.0) + + for point in actual_points: + self.assertIn(tuple(np.round(point, 8)), candidate_norm) + + if doe_obj.objective_option in DesignOfExperiments._MAXIMIZE_OBJECTIVES: + best_obj = -np.inf + is_better = lambda new, best: new > best + else: + best_obj = np.inf + is_better = lambda new, best: new < best + + # Build exhaustive reference over all candidate pairs. + candidate_fims = [ + doe_obj._compute_fim_at_point_no_prior(0, list(point)) + for point in candidate_points + ] + for combo in itertools.combinations(range(len(candidate_points)), 2): + fim_total = sum((candidate_fims[idx] for idx in combo), np.zeros((2, 2))) + obj_val = doe_obj._evaluate_objective_from_fim(fim_total) + if is_better(obj_val, best_obj): + best_obj = obj_val + + actual_fim_total = sum( + ( + doe_obj._compute_fim_at_point_no_prior(0, list(point)) + for point in actual_points + ), + np.zeros((2, 2)), + ) + actual_obj = doe_obj._evaluate_objective_from_fim(actual_fim_total) + + self.assertAlmostEqual(actual_obj, best_obj, places=12) + self.assertAlmostEqual(lhs_diag["best_obj"], best_obj, places=12) + + def test_optimize_experiments_greybox_initialization_refreshes_inputs_after_square_solve( + self, + ): + # Build-time greybox inputs reflect the template design, but after LHS + # changes the starting point the square-solve refresh should reseed the + # block from the new aggregated FIM before the final solve. + lhs_n_samples = 4 + # Set random seed to keep LHS initialization deterministic. + lhs_seed = 29 + captured = {} + doe_obj = _make_multiexperiment_greybox_doe( + objective_option="minimum_eigenvalue", + prior_FIM=np.zeros((2, 2)), + grey_box_solver=_MockGreyBoxSolver(), + ) + original_initialize_impl = doe_obj._initialize_grey_box_block + + def _capture_initialize(egb_block, fim_np, parameter_names): + captured["before"] = _reconstruct_fim_from_egb_inputs( + egb_block, parameter_names + ) + captured["refreshed"] = np.array(fim_np, copy=True) + return original_initialize_impl(egb_block, fim_np, parameter_names) + + original_compute_fim = doe_obj._compute_fim_at_point_no_prior + original_initialize_for_restore = doe_obj._initialize_grey_box_block + doe_obj._compute_fim_at_point_no_prior = _diagonal_hour_fim_reference + doe_obj._initialize_grey_box_block = _capture_initialize + try: + doe_obj.optimize_experiments( + n_exp=1, + init_method="latin_hypercube_sampling", + init_n_samples=lhs_n_samples, + init_seed=lhs_seed, + ) + finally: + doe_obj._compute_fim_at_point_no_prior = original_compute_fim + doe_obj._initialize_grey_box_block = original_initialize_for_restore + + scenario, total_fim, parameter_names = _get_multiexperiment_scenario_data( + doe_obj + ) + final_fim = _reconstruct_fim_from_egb_inputs( + scenario.obj_cons.egb_fim_block, parameter_names + ) + + self.assertIn("before", captured) + self.assertIn("refreshed", captured) + self.assertFalse(np.allclose(captured["before"], captured["refreshed"])) + self.assertTrue(np.allclose(final_fim, captured["refreshed"], atol=1e-7)) + self.assertTrue(np.allclose(final_fim, total_fim, atol=1e-7)) + # The Rooney-Biegler template is built with hour=2.0, so this confirms + # LHS moved away from the build-time design and made the refresh check meaningful. + self.assertNotAlmostEqual( + doe_obj.results["initialization"]["selected_initial_designs"][0][0], + 2.0, + places=6, + ) + + def test_optimize_experiments_greybox_tee_flag_reaches_solver(self): + # grey_box_tee is only meaningful if it propagates to the external + # solver interface, so capture the mocked solve() call and verify tee. + grey_box_solver = _MockGreyBoxSolver() + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="minimum_eigenvalue", + use_grey_box_objective=True, + step=1e-2, + solver=make_ipopt_solver(), + grey_box_solver=grey_box_solver, + grey_box_tee=True, + ) + + doe_obj.optimize_experiments(n_exp=2) + + self.assertEqual(len(grey_box_solver.calls), 1) + self.assertTrue(grey_box_solver.calls[0]["tee"]) + + +@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not scipy_available, "scipy is not available") +@unittest.skipIf(not cyipopt_available, "'cyipopt' is not available") +@unittest.skipIf( + not cyipopt_call_working, "cyipopt is not properly accessing linear solvers" +) +@unittest.skipIf(not pandas_available, "pandas is not available") +class TestMultiexperimentSolve(unittest.TestCase): + def test_optimize_experiments_determinant_expected_values_greybox(self): + # The scenario objective is order-invariant, so compare the chosen + # experiment hours in sorted order. + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.0, y=8.3), + RooneyBieglerMultiExperiment(hour=2.0, y=10.3), + ] + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="determinant", + step=1e-2, + use_grey_box_objective=True, + grey_box_solver=_make_cyipopt_solver(tol=1e-4), + grey_box_tee=False, + ) + doe.optimize_experiments() + + scenario = doe.results["solution"]["param_scenarios"][0] + got_hours = sorted(exp["design"][0] for exp in scenario["experiments"]) + self.assertStructuredAlmostEqual( + got_hours, sorted([1.9321985035514362, 9.999999685577139]), abstol=1e-3 + ) + self.assertAlmostEqual( + scenario["quality_metrics"]["log10_d_opt"], 6.028152580313302, places=3 + ) + + def test_optimize_experiments_trace_expected_values_greybox(self): + # This is a regression test for the cyipopt-backed multi-experiment + # greybox solve on A-opt: the chosen design and reported objective + # should stay near the known Rooney-Biegler reference solution. + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.0, y=8.3), + RooneyBieglerMultiExperiment(hour=2.0, y=10.3), + ] + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="trace", + step=1e-2, + use_grey_box_objective=True, + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + grey_box_tee=False, + ) + doe.optimize_experiments() + + scenario = doe.results["solution"]["param_scenarios"][0] + got_hours = sorted(exp["design"][0] for exp in scenario["experiments"]) + self.assertStructuredAlmostEqual(got_hours, sorted([1.01, 10.0]), abstol=1e-3) + self.assertAlmostEqual( + scenario["quality_metrics"]["log10_a_opt"], -1.9438, places=3 + ) + + def test_optimize_experiments_min_eig_expected_values_greybox(self): + # This checks the end-to-end greybox E-opt solve against a stable + # reference solution so future greybox wiring changes do not silently + # alter the chosen experiment pair or final metric. + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.0, y=8.3), + RooneyBieglerMultiExperiment(hour=2.0, y=10.3), + ] + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="minimum_eigenvalue", + step=1e-2, + solver=make_ipopt_solver(), + use_grey_box_objective=True, + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + grey_box_tee=False, + ) + doe.optimize_experiments() + + scenario = doe.results["solution"]["param_scenarios"][0] + got_hours = sorted(exp["design"][0] for exp in scenario["experiments"]) + self.assertStructuredAlmostEqual(got_hours, sorted([1.0, 10.0]), abstol=1e-2) + self.assertAlmostEqual( + scenario["quality_metrics"]["log10_e_opt"], 1.9532, places=2 + ) + + def test_optimize_experiments_me_opt_expected_values_greybox(self): + # ME-opt is greybox-only in optimize_experiments(), so keep a dedicated + # solve regression here to guard the condition-number objective path. + exp_list = [ + RooneyBieglerMultiExperiment(hour=1.0, y=8.3), + RooneyBieglerMultiExperiment(hour=2.0, y=10.3), + ] + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="condition_number", + step=1e-2, + solver=make_ipopt_solver(), + use_grey_box_objective=True, + grey_box_solver=_make_cyipopt_solver(tol=1e-6), + grey_box_tee=False, + ) + doe.optimize_experiments() + + scenario = doe.results["solution"]["param_scenarios"][0] + got_hours = sorted(exp["design"][0] for exp in scenario["experiments"]) + self.assertStructuredAlmostEqual(got_hours, sorted([7.13, 10.0]), abstol=1e-2) + self.assertAlmostEqual( + scenario["quality_metrics"]["log10_me_opt"], 1.5289, places=2 + ) + + +@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not scipy_available, "scipy is not available") +@unittest.skipIf(not cyipopt_available, "'cyipopt' is not available") +@unittest.skipIf( + not cyipopt_call_working, "cyipopt is not properly accessing linear solvers" +) +@unittest.skipIf(not pandas_available, "pandas is not available") +class TestSingleExperimentSolve(unittest.TestCase): + def test_optimize_experiments_single_experiment_greybox_path_works(self): + # Even with n_exp=1, optimize_experiments() should take the template-mode + # greybox path, solve with the grey_box_solver, and keep the greybox + # block synchronized with the final scenario FIM. + grey_box_solver = _make_cyipopt_solver(tol=1e-6) + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option="minimum_eigenvalue", + step=1e-2, + solver=make_ipopt_solver(), + use_grey_box_objective=True, + grey_box_solver=grey_box_solver, + grey_box_tee=False, + ) + + doe_obj.optimize_experiments(n_exp=1) + + scenario, total_fim, _ = _get_multiexperiment_scenario_data(doe_obj) + scenario_results = doe_obj.results["solution"]["param_scenarios"][0] + design_hour = scenario_results["experiments"][0]["design"][0] + + self.assertEqual(doe_obj.results["optimization_solve"]["status"], "ok") + self.assertEqual( + doe_obj.results["problem"]["number_of_experiments_per_scenario"], 1 + ) + self.assertTrue(doe_obj.results["problem"]["used_template_experiment"]) + self.assertEqual(len(scenario_results["experiments"]), 1) + self.assertEqual( + doe_obj.results["optimization_solve"]["solver"], + getattr(grey_box_solver, "name", str(grey_box_solver)), + ) + self.assertGreaterEqual(design_hour, 1.0) + self.assertLessEqual(design_hour, 10.0) + self.assertTrue(np.isfinite(scenario_results["quality_metrics"]["log10_e_opt"])) + self.assertAlmostEqual( + pyo.value(scenario.obj_cons.egb_fim_block.outputs["E-opt"]), + _expected_multiexperiment_greybox_output("minimum_eigenvalue", total_fim), + places=7, + ) + + +@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") +@unittest.skipIf(not numpy_available, "Numpy is not available") +@unittest.skipIf(not scipy_available, "scipy is not available") +@unittest.skipIf(not pandas_available, "pandas is not available") +@unittest.skipIf( + not parameterized_available, "The 'parameterized' package is not available" +) +class TestMultiexperimentError(unittest.TestCase): + def test_optimize_experiments_greybox_unsupported_objectives_are_rejected(self): + # These unsupported objectives share the same early-validation path, so + # keep them in one table-driven test and verify none reaches the + # external grey_box_solver interface. + objective_option = "zero" + + class _UnusedGreyBoxSolver: + def solve(self, model, tee=False): + raise AssertionError("grey_box_solver.solve should not be reached") + + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option=objective_option, + use_grey_box_objective=True, + step=1e-2, + solver=make_ipopt_solver(), + grey_box_solver=_UnusedGreyBoxSolver(), + ) + + with self.assertRaisesRegex( + ValueError, + "Grey-box objective support in optimize_experiments\\(\\) is only " + "available", + ): + doe_obj.optimize_experiments(n_exp=2) + + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/doe/tests/utils_for_doe_tests.py b/pyomo/contrib/doe/tests/utils_for_doe_tests.py new file mode 100644 index 00000000000..cf0dc938000 --- /dev/null +++ b/pyomo/contrib/doe/tests/utils_for_doe_tests.py @@ -0,0 +1,18 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 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.opt import SolverFactory + + +def make_ipopt_solver(): + """Create an IPOPT solver configured for DOE test runs.""" + solver = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + return solver