diff --git a/pyomo/contrib/doe/documentation.md b/pyomo/contrib/doe/documentation.md new file mode 100644 index 00000000000..245c53f44c3 --- /dev/null +++ b/pyomo/contrib/doe/documentation.md @@ -0,0 +1,185 @@ +# DesignOfExperiments `optimize_experiments()` Proposed Interface + +This is a brief documentation explaining the changes. This document will not get merged. + +## Why This Change + +The `optimize_experiments()` path has grown from a single-experiment workflow into a +multi-experiment optimization interface with optional initialization strategies and +richer result payloads. This document captures the interface contract, operating modes, +and optimization model at a high level. + +## API Summary + +### Constructor + +```python +DesignOfExperiments( + experiment=..., # single experiment object OR list of experiment objects + ... +) +``` + +- `experiment` now accepts either: + - one experiment object (template mode input), or + - a list of experiment objects (user-initialized mode input). +- Internally, the implementation normalizes this to `self.experiment_list`. + +### Multi-Experiment Solve Entry Point + +```python +optimize_experiments( + results_file=None, + n_exp: int = None, + init_method: InitializationMethod = None, + init_n_samples: int = 5, + init_seed: int = None, + init_parallel: bool = False, + init_combo_parallel: bool = False, + init_n_workers: int = None, + init_combo_chunk_size: int = 5000, + init_combo_parallel_threshold: int = 20000, + init_max_wall_clock_time: float = None, + init_solver=None, +) +``` + +### Short description of arguments + +- `results_file`: Optional path for writing JSON results. +- `n_exp`: Number of experiments to optimize in template mode. +- `init_method`: Initialization strategy (`None` or `"lhs"`). +- `init_n_samples`: Number of LHS samples per experiment-input dimension. +- `init_seed`: Random seed used by LHS initialization. +- `init_parallel`: Enables parallel candidate-point FIM evaluation. +- `init_combo_parallel`: Enables parallel candidate-combination scoring. +- `init_n_workers`: Worker count for LHS parallel evaluation/scoring paths. +- `init_combo_chunk_size`: Number of combinations handled per worker task. +- `init_combo_parallel_threshold`: Minimum combinations required before using combo parallelism. +- `init_max_wall_clock_time`: Optional LHS time budget (seconds); returns best-so-far if exceeded. +- `init_solver`: Optional solver used only during initialization phases (LHS and square init solve). + +## Operating Modes + +### 1) Template Mode + +- Condition: `len(experiment) == 1` +- `n_exp` may be provided to choose how many experiments to optimize simultaneously. +- If `n_exp` is omitted, default is `1`. + +### 2) User-Initialized Mode + +- Condition: `len(experiment) > 1` +- Number of experiments is fixed by the list length. +- `n_exp` must not be provided. + +## Initialization Behavior + +### No Special Initialization (`init_method=None`) + +- Uses current experiment design values from the model labels directly. + +### LHS Initialization (`init_method="lhs"`) + +- Currently supported in template mode. +- Requires explicit lower and upper bounds for all experiment inputs. +- Generates candidate points using per-dimension 1-D LHS, then Cartesian product. +- Scores combinations of candidate points using objective-specific FIM metrics. +- Selects best-scoring set of initial points for the nonlinear solve. + +### `init_solver` (new) + +- If provided, `init_solver` is used for: + - initialization-phase solves (including multi-experiment block-construction + solves and the LHS candidate FIM evaluation path), + - the square initialization solve before final optimization. +- Final optimization solve still uses the main DoE solver (`self.solver`). +- If `init_solver` is `None`, initialization also uses `self.solver`. + +## Optimization Formulation + +The proposed multi-experiment interface follows a simultaneous design formulation. + +### General Form + +Let: + +- $E = \{1, 2, ..., N_{exp}\}$ be the experiment index set, +- $\phi_k$ be the design variables for experiment `k`, +- $\theta$ be the model parameters +- $M_0$ be the prior Fisher Information Matrix (FIM), +- $M_k(\hat{\theta}, \phi_k)$ be the FIM contribution from experiment `k`, +- $\Psi(M)$ be the chosen FIM metric (D-, A-, pseudo-A-, etc.). + +Then: + +```math +\max_{\phi_1,\ldots,\phi_{N_{exp}}} \Psi(\mathbf{M}) +``` + +subject to: + +```math +\mathbf{M} = \sum_{k=1}^{N_{exp}} \mathbf{M}_k(\hat{\theta}, \phi_k) + \mathbf{M}_0 +``` + +```math +\mathbf{M}_k = \mathbf{Q}_k^\top \Sigma_{\bar{y},k}^{-1} \mathbf{Q}_k, \quad \forall k \in E +``` + +```math +\mathbf{m}(\bar{x}_k, \hat{\bar{y}}_k, \phi_k, \hat{\theta}, t) = 0, \quad \forall k \in E +``` + +```math +\mathbf{g}(\bar{x}_k, \hat{\bar{y}}_k, \phi_k, \hat{\theta}, t) \le 0, \quad \forall k \in E +``` + +where $\mathbf{Q}_k$ is the sensitivity matrix for experiment `k`. + +### Symmetry-Breaking Constraint + +To avoid permutation-equivalent solutions in simultaneous design: + +```math +\varphi_{\text{primary},1} \le \varphi_{\text{primary},2} \le \cdots \le \varphi_{\text{primary},N_{exp}} +``` +Here, $\varphi_{\text{primary},\cdot}$ is the primary variable which is considered for symmetry breaking. + +This is implemented in `optimize_experiments()` by using a user-marked primary design +variable passed in a `Pyomo.Suffix` (or a default selection (first variable from +experiment_inputs Suffix) with warning if not marked). + +### Current Implementation Specialization + +The current implementation corresponds to the single-scenario case (`N_s = 1`) with: + +```math +\mathbf{M}_{\text{total}} = \mathbf{M}_0 + \sum_{k=1}^{N_{exp}} \mathbf{M}_k +``` +Future implementation will handle parametric uncertainty with $N_s >1$ + +Objective handling in code: + +- Determinant and pseudo-trace: solved in maximization form (with monotonic transforms used in the NLP objective expressions). +- Trace: solved in minimization form via covariance/FIM-inverse representation. +- Zero objective: feasibility/debug mode. + +Cholesky-based constraints and variables are used in supported objective paths to stabilize determinant/trace formulations. + +## Result Payload Highlights + +The solver output includes both legacy fields and structured fields: + +- `run_info` (API name, solver status/termination info) +- `settings` (objective, finite-difference, initialization config, modeling mode) +- `timing` (build/initialization/solve/total timing) +- `names` (design/output/parameter/error labels) +- `diagnostics` (symmetry and LHS diagnostics) +- `scenarios` (scenario-level objective metrics, total FIM, per-experiment details) + +Notably, initialization settings now include: + +- `settings["initialization"]["solver_name"]` + +to make it explicit which solver was used during initialization. diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py index a13e6785bf4..1f0acf3789d 100644 --- a/pyomo/contrib/doe/doe.py +++ b/pyomo/contrib/doe/doe.py @@ -25,11 +25,20 @@ # ____________________________________________________________________________________ from enum import Enum -from itertools import permutations, product - +from itertools import ( + permutations, + product, + combinations as _combinations, + islice as _islice, +) +import concurrent.futures as _cf import json import logging import math +import os +import threading +import time +import warnings from pyomo.common.dependencies import ( numpy as np, @@ -49,6 +58,7 @@ from pyomo.contrib.doe.grey_box_utilities import FIMExternalGreyBox from pyomo.contrib.pynumero.interfaces.external_grey_box import ExternalGreyBoxBlock + from scipy.stats.qmc import LatinHypercube import pyomo.environ as pyo from pyomo.contrib.doe.utils import ( @@ -56,6 +66,7 @@ compute_FIM_metrics, _SMALL_TOLERANCE_DEFINITENESS, ) +from pyomo.contrib.parmest.utils.model_utils import update_model_from_suffix from pyomo.opt import SolverStatus @@ -75,10 +86,37 @@ class FiniteDifferenceStep(Enum): backward = "backward" +class InitializationMethod(Enum): + lhs = "lhs" + + +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", @@ -109,13 +147,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 @@ -129,7 +171,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`` @@ -186,17 +229,33 @@ def __init__( Specify the level of the logger. Change to logging.DEBUG for all messages. """ + # Validate experiment(s) are provided 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"): raise ValueError( - "The experiment object must have a ``get_labeled_model`` function" + "The 'experiment' parameter is required. " + "Pass a single Experiment object or a 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 + + # Validate list is not empty + if len(experiment_list) == 0: + raise ValueError("The 'experiment' list cannot be empty.") + + # 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) @@ -269,6 +328,11 @@ 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] + # Perform doe def run_doe(self, model=None, results_file=None): """ @@ -285,11 +349,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 +394,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 +420,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() @@ -409,7 +472,7 @@ def run_doe(self, model=None, results_file=None): # 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)) + min_eig = np.min(np.real(np.linalg.eigvals(fim_np))) if min_eig < _SMALL_TOLERANCE_DEFINITENESS: # Raise the minimum eigenvalue to at @@ -478,10 +541,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 @@ -489,29 +556,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)) @@ -525,7 +592,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 @@ -537,6 +604,1583 @@ def run_doe(self, model=None, results_file=None): with open(results_file, "w") as file: json.dump(self.results, file) + 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_parallel: bool = False, + init_combo_parallel: bool = False, + init_n_workers: int = None, + init_combo_chunk_size: int = 5000, + init_combo_parallel_threshold: int = 20000, + init_max_wall_clock_time: float = 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``. + - ``"lhs"`` (or ``InitializationMethod.lhs``): 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="lhs"``. 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="lhs"``. + Default: ``None`` (non-deterministic). + init_parallel: + If True, evaluate candidate-point FIMs in parallel during LHS + initialization. Default: False. + init_combo_parallel: + If True, the scoring of Latin hypercube candidate combinations + (``C(n_candidates, n_exp)`` during ``init_method="lhs"``) + is split across a thread pool. Each worker computes the scalar + objective derived from the FIM for its chunk of combinations. The + flag has no effect unless ``init_method="lhs"`` and the + total number of combinations exceeds ``init_combo_parallel_threshold``. + Default: False. + init_n_workers: + Number of worker threads for combination FIM metric when + ``init_combo_parallel=True``. Default: ``None`` (auto-select). + init_combo_chunk_size: + Number of combinations scored per worker task. Default: 5000. + init_combo_parallel_threshold: + Parallel combo scoring is used only when number of combinations is + at least this value. Default: 20000. + init_max_wall_clock_time: + Optional time budget (seconds) for LHS initialization. If exceeded + during combination scoring, best-so-far is returned. + 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="lhs"): + 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. + + Solver options in LHS worker evaluations: + When ``init_parallel=True``, worker threads construct solver instances + using the same solver name and options as ``self.solver`` (when + available). Therefore, per-solve limits (e.g., iteration/time limits) + configured on the main DoE solver are propagated to candidate FIM + evaluations. + """ + # Check results file name + if results_file is not None: + if not isinstance(results_file, (pathlib.Path, str)): + raise ValueError( + "``results_file`` must be either a Path object or a string." + ) + + # --- Resolve n_exp and determine operating mode --- + n_list = len(self.experiment_list) + if n_list > 1: + # User-initialized mode: experiment list already contains all + # pre-initialized experiment objects. + 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``." + ) + n_exp = n_list + _template_mode = False + else: + # Template mode: single experiment object cloned n_exp times. + if n_exp is None: + n_exp = 1 # default: single-experiment optimization + elif not isinstance(n_exp, int) or n_exp < 1: + raise ValueError( + f"``n_exp`` must be a positive integer, got {n_exp!r}." + ) + _template_mode = True + # --------------------------------------------------- + + # --- Validate initialization arguments --- + 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.lhs: + if not _template_mode: + raise ValueError( + "``init_method='lhs'`` 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='lhs'." + ) + 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 not isinstance(init_parallel, bool): + raise ValueError( + f"``init_parallel`` must be a bool, got {init_parallel!r}." + ) + if not isinstance(init_combo_parallel, bool): + raise ValueError( + "``init_combo_parallel`` must be a bool, " + f"got {init_combo_parallel!r}." + ) + if init_n_workers is not None and ( + not isinstance(init_n_workers, int) or init_n_workers < 1 + ): + raise ValueError( + "``init_n_workers`` must be None or a positive integer, " + f"got {init_n_workers!r}." + ) + if not isinstance(init_combo_chunk_size, int) or init_combo_chunk_size < 1: + raise ValueError( + "``init_combo_chunk_size`` must be a positive integer, " + f"got {init_combo_chunk_size!r}." + ) + if ( + not isinstance(init_combo_parallel_threshold, int) + or init_combo_parallel_threshold < 1 + ): + raise ValueError( + "``init_combo_parallel_threshold`` must be a positive integer, " + f"got {init_combo_parallel_threshold!r}." + ) + if init_max_wall_clock_time is not None and ( + not isinstance(init_max_wall_clock_time, (int, float)) + or not math.isfinite(init_max_wall_clock_time) + or init_max_wall_clock_time <= 0 + ): + raise ValueError( + "``init_max_wall_clock_time`` must be None or a positive number, " + f"got {init_max_wall_clock_time!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 + init_solver_name = getattr( + resolved_init_solver, "name", str(resolved_init_solver) + ) + lhs_init_diagnostics = None + lhs_initialization_time = 0.0 + + # Start timer + sp_timer = TicTocTimer() + sp_timer.tic(msg=None) + self.logger.info( + f"Beginning multi-experiment optimization with {n_exp} experiments." + ) + # Rebuild the multi-experiment model from a clean base each call. + self.model = pyo.ConcreteModel() + + n_param_scenarios = 1 # currently single-scenario optimization + # Use an immutable tuple since weights are not intended to be modified + self.scenario_weights = (1.0,) # Single scenario, weight = 1 + # Add parameter scenario blocks to the model + self.model.param_scenario_blocks = pyo.Block(range(n_param_scenarios)) + symmetry_breaking_info = { + "enabled": n_exp > 1, + "variable": None, + "source": None, + } + diagnostics_warnings = [] + # Route all pre-final solves through the initialization solver, then + # restore the primary solver for the final optimization solve. + self.solver = resolved_init_solver + try: + # Add experiment(s) for each scenario + # TODO: Add s_prev = 0 to handle parameter scenarios + for s in range(n_param_scenarios): + self.model.param_scenario_blocks[s].exp_blocks = pyo.Block(range(n_exp)) + for k in range(n_exp): + # Generate FIM and Sensitivity expressions for each experiment. + # In template mode all experiments share the single template + # (experiment_index=0); in user-initialized mode each experiment + # maps to its own entry in experiment_list (experiment_index=k). + self.create_doe_model( + model=self.model.param_scenario_blocks[s].exp_blocks[k], + experiment_index=0 if _template_mode else k, + _for_multi_experiment=True, # Skip creating L matrix per experiment + ) + # TODO: Update the parameter scenarios for each experiment block + # when using parametric uncertainty + + # Add symmetry breaking constraints to prevent equivalent permutations for + # multiple experiments + if n_exp > 1: + # Check if user provided a symmetry breaking variable via Suffix + # Use first scenario since variable names are the same across all scenarios + first_exp_block = ( + self.model.param_scenario_blocks[0] + .exp_blocks[0] + .fd_scenario_blocks[0] + ) + + # Determine symmetry breaking variable + if ( + hasattr(first_exp_block, 'sym_break_cons') + and len(first_exp_block.sym_break_cons) > 0 + ): + # User provided symmetry breaking variable(s) + sym_break_var_list = list(first_exp_block.sym_break_cons.keys()) + + if len(sym_break_var_list) > 1: + warning_msg = ( + f"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: + # Use first experiment input as default symmetry breaking variable + 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( + f"No symmetry breaking variable specified. Automatically using " + f"'{sym_break_var.local_name}'." + ) + + # Add constraints for each scenario + for s in range(n_param_scenarios): + for k in range(1, n_exp): + # Get the variable from experiment k-1 + 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] + ) + + # Get the variable from experiment k + 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." + ) + + # Add symmetry breaking constraint + 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}" + ) + + # Create aggregated objective for multi-experiment optimization + self.create_multi_experiment_objective_function(self.model) + + # Track time required to build the DoE 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 + ) + + # --- Apply experiment initialization (if requested) --- + # This must be done AFTER the model is built but BEFORE the square solve + # so that the solver uses the correct starting design. + if resolved_init_method == InitializationMethod.lhs: + 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, + lhs_parallel=init_parallel, + lhs_combo_parallel=init_combo_parallel, + lhs_n_workers=init_n_workers, + lhs_combo_chunk_size=init_combo_chunk_size, + lhs_combo_parallel_threshold=init_combo_parallel_threshold, + lhs_max_wall_clock_time=init_max_wall_clock_time, + ) + ) + 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) + + # ------------------------------------------------------ + # Reset delta timing so initialization_time measures only square solve. + sp_timer.tic(msg=None) + + # Solve the square problem first to initialize the FIM and sensitivity constraints + # Deactivate objective expression and objective constraints + self.model.objective.deactivate() + # Deactivate obj_cons for each scenario (holds Cholesky/determinant/trace constraints) + 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() + + # Fix the design variables across all scenarios and experiments + self._set_experiment_inputs_fixed( + n_param_scenarios=n_param_scenarios, n_exp=n_exp, fixed=True + ) + + # Create and solve dummy objective for initialization + self.model.dummy_obj = pyo.Objective(expr=0, sense=pyo.minimize) + self.solver.solve(self.model, tee=self.tee) + + # Track time to initialize the DoE model + initialization_time = sp_timer.toc(msg=None) + self.logger.info( + ( + "Successfully initialized the multi-experiment DoE model." + "\nInitialization time: %0.1f seconds" % initialization_time + ) + ) + + # Deactivate dummy objective + self.model.dummy_obj.deactivate() + self.model.del_component("dummy_obj") + + # Reactivate objective, obj_cons, and unfix experimental design decisions + 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() + + # Initialize scenario-level variables (L, determinant, pseudo_trace) based on + # the solved FIM values from the square solve + parameter_names = ( + self.model.param_scenario_blocks[0].exp_blocks[0].parameter_names + ) + + for s in range(n_param_scenarios): + scenario = self.model.param_scenario_blocks[s] + # Update total_fim values from solved individual experiment FIMs + # Individual experiment FIMs don't include prior_FIM in multi-experiment mode, + # so we add it once here to the aggregated total + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + # When only_compute_fim_lower=True, only the lower triangle is computed + # Upper triangle elements are fixed to 0, so only aggregate lower triangle + 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] + ) + + # Build scenario total FIM once and reuse for all objective initializations. + 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) + + # Initialize scenario-level variables based on total_fim + if hasattr(scenario.obj_cons, "L"): + # Compute Cholesky factorization + # Check positive definiteness and add jitter if needed + 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)) + + # Compute Cholesky decomposition + L_vals = np.linalg.cholesky(fim_jittered) + + # Initialize L values + for i, p in enumerate(parameter_names): + for j, q in enumerate(parameter_names): + scenario.obj_cons.L[p, q].set_value(L_vals[i, j]) + + # If trace objective, also initialize L_inv, fim_inv, and cov_trace + if hasattr(scenario.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: + scenario.obj_cons.L_inv[p, q].set_value( + L_inv_vals[i, j] + ) + else: + scenario.obj_cons.L_inv[p, q].set_value(0.0) + + # Initialize fim_inv + if hasattr(scenario.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): + scenario.obj_cons.fim_inv[p, q].set_value( + fim_inv_np[i, j] + ) + + # Initialize cov_trace + if hasattr(scenario.obj_cons, "cov_trace"): + initial_cov_trace = np.sum(L_inv_vals**2) + scenario.obj_cons.cov_trace.set_value(initial_cov_trace) + + if hasattr(scenario.obj_cons, "determinant"): + # Initialize determinant + scenario.obj_cons.determinant.set_value(np.linalg.det(total_fim_np)) + + if hasattr(scenario.obj_cons, "pseudo_trace"): + # Initialize pseudo_trace + pseudo_trace_val = float(np.trace(total_fim_np)) + scenario.obj_cons.pseudo_trace.set_value(pseudo_trace_val) + finally: + self.solver = primary_solver + + # Solve the full model + res = self.solver.solve(self.model, tee=self.tee) + + # Track time used to solve the DoE model + 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) + ) + + # Collect results + self.results = {} + self.results["Solver Status"] = res.solver.status + self.results["Termination Condition"] = 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 "" + ) + self.results["Termination Message"] = results_message + + 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") + + # Store results for each scenario + self.results["Scenarios"] = [] + scenarios_structured = [] + for s in range(n_param_scenarios): + scenario = self.model.param_scenario_blocks[s] + scenario_results = {} + + # Get aggregated FIM for this scenario + 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)) + ) + + # Complete FIM if only computing lower triangle + if self.only_compute_fim_lower: + total_fim_np = self._symmetrize_lower_tri(total_fim_np) + + # Store the completed (symmetric) FIM + scenario_results["Total FIM"] = total_fim_np.tolist() + + # Statistics on the aggregated FIM (consistent with run_doe), guarded + # against singular/indefinite matrices and numerical failures. + scenario_results["log10 A-opt"] = _safe_metric( + "log10 A-opt", + lambda fim=total_fim_np: np.log10(np.trace(np.linalg.inv(fim))), + s, + ) + scenario_results["log10 pseudo A-opt"] = _safe_metric( + "log10 pseudo A-opt", + lambda fim=total_fim_np: np.log10(np.trace(fim)), + s, + ) + scenario_results["log10 D-opt"] = _safe_metric( + "log10 D-opt", lambda fim=total_fim_np: np.log10(np.linalg.det(fim)), s + ) + scenario_results["log10 E-opt"] = _safe_metric( + "log10 E-opt", + lambda fim=total_fim_np: np.log10(min(np.linalg.eigvalsh(fim))), + s, + ) + scenario_results["log10 ME-opt"] = _safe_metric( + "log10 ME-opt", + lambda fim=total_fim_np: np.log10(np.linalg.cond(fim)), + s, + ) + + # Store unknown parameter values at scenario level (same for all experiments) + # Use first experiment to get the values + scenario_results["Unknown Parameters"] = self.get_unknown_parameter_values( + model=scenario.exp_blocks[0] + ) + + # Store results for each experiment in this scenario + scenario_results["Experiments"] = [] + scenario_structured = { + "id": s, + "total_fim": total_fim_np.tolist(), + "metrics": { + "log10_a_opt": scenario_results["log10 A-opt"], + "log10_pseudo_a_opt": scenario_results["log10 pseudo A-opt"], + "log10_d_opt": scenario_results["log10 D-opt"], + "log10_e_opt": scenario_results["log10 E-opt"], + "log10_me_opt": scenario_results["log10 ME-opt"], + }, + "unknown_parameters": None, + "experiments": [], + } + for k in range(n_exp): + exp_block = scenario.exp_blocks[k] + exp_results = {} + + # Use helper functions for consistent extraction (same as run_doe) + # Store only the VALUES for each experiment (names are at top level) + exp_results["Experiment Design"] = self.get_experiment_input_values( + model=exp_block + ) + exp_results["Experiment Outputs"] = self.get_experiment_output_values( + model=exp_block + ) + exp_results["Measurement Error"] = self.get_measurement_error_values( + model=exp_block + ) + + # Individual experiment FIM (get_FIM handles symmetry completion) + exp_results["FIM"] = self.get_FIM(model=exp_block) + + # Sensitivity matrix for this experiment + if hasattr(exp_block, "sensitivity_jacobian"): + exp_results["Sensitivity Matrix"] = self.get_sensitivity_matrix( + model=exp_block + ) + + scenario_results["Experiments"].append(exp_results) + experiment_structured = { + "id": k, + "design": exp_results["Experiment Design"], + "outputs": exp_results["Experiment Outputs"], + "measurement_error": exp_results["Measurement Error"], + "fim": exp_results["FIM"], + } + if "Sensitivity Matrix" in exp_results: + experiment_structured["sensitivity"] = exp_results[ + "Sensitivity Matrix" + ] + scenario_structured["experiments"].append(experiment_structured) + + self.results["Scenarios"].append(scenario_results) + scenario_structured["unknown_parameters"] = scenario_results[ + "Unknown Parameters" + ] + scenarios_structured.append(scenario_structured) + + # Store variable names once (structural properties, same across all scenarios/experiments) + # Use first scenario's first experiment to get the structure + first_exp_block_fd = ( + self.model.param_scenario_blocks[0].exp_blocks[0].fd_scenario_blocks[0] + ) + self.results["Experiment Design Names"] = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.experiment_inputs + ] + self.results["Experiment Output Names"] = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.experiment_outputs + ] + self.results["Unknown Parameter Names"] = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.unknown_parameters + ] + self.results["Measurement Error Names"] = [ + str(pyo.ComponentUID(comp, context=first_exp_block_fd)) + for comp in first_exp_block_fd.measurement_error + ] + + # Store general settings and info + self.results["Number of Scenarios"] = n_param_scenarios + self.results["Number of Experiments per Scenario"] = n_exp + self.results["Prior FIM"] = self.prior_FIM.tolist() + self.results["Objective expression"] = self._enum_label(self.objective_option) + 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 + + # Initialization info + self.results["Initialization Method"] = ( + resolved_init_method.value if resolved_init_method is not None else "none" + ) + if resolved_init_method == InitializationMethod.lhs: + self.results["LHS Samples Per Dimension"] = init_n_samples + self.results["LHS Seed"] = init_seed + self.results["LHS Best Initial Points"] = best_initial_points + + # Timing statistics + self.results["Build Time"] = build_time + self.results["Initialization Time"] = initialization_time + self.results["LHS Initialization Time"] = lhs_initialization_time + self.results["Solve Time"] = solve_time + self.results["Wall-clock Time"] = ( + build_time + lhs_initialization_time + initialization_time + solve_time + ) + + # Structured result payload + objective_sense = ( + "maximize" + if self.objective_option in self._MAXIMIZE_OBJECTIVES + else "minimize" + ) + + self.results["run_info"] = { + "api": "DesignOfExperiments.optimize_experiments", + "solver": { + "name": getattr(self.solver, "name", str(self.solver)), + "status": self.results["Solver Status"], + "termination_condition": self.results["Termination Condition"], + "message": self.results["Termination Message"], + }, + } + self.results["settings"] = { + "objective": { + "name": self.results["Objective expression"], + "sense": objective_sense, + }, + "finite_difference": { + "scheme": self.results["Finite Difference Scheme"], + "step": self.results["Finite Difference Step"], + }, + "scaling": { + "nominal_parameter_scaling": self.results["Nominal Parameter Scaling"] + }, + "initialization": { + "method": self.results["Initialization Method"], + "solver_name": init_solver_name, + "lhs_n_samples": self.results.get("LHS Samples Per Dimension"), + "lhs_seed": self.results.get("LHS Seed"), + "best_points": self.results.get("LHS Best Initial Points"), + "lhs_parallel": ( + init_parallel + if resolved_init_method == InitializationMethod.lhs + else None + ), + "lhs_combo_parallel": ( + init_combo_parallel + if resolved_init_method == InitializationMethod.lhs + else None + ), + "lhs_n_workers": ( + init_n_workers + if resolved_init_method == InitializationMethod.lhs + else None + ), + "lhs_combo_chunk_size": ( + init_combo_chunk_size + if resolved_init_method == InitializationMethod.lhs + else None + ), + "lhs_combo_parallel_threshold": ( + init_combo_parallel_threshold + if resolved_init_method == InitializationMethod.lhs + else None + ), + "lhs_max_wall_clock_time": ( + init_max_wall_clock_time + if resolved_init_method == InitializationMethod.lhs + else None + ), + }, + "modeling": { + "n_scenarios": self.results["Number of Scenarios"], + "n_experiments_per_scenario": self.results[ + "Number of Experiments per Scenario" + ], + "template_mode": _template_mode, + }, + "prior_fim": self.results["Prior FIM"], + } + self.results["timing"] = { + "build_s": self.results["Build Time"], + "lhs_initialization_s": self.results["LHS Initialization Time"], + "initialization_s": self.results["Initialization Time"], + "solve_s": self.results["Solve Time"], + "total_s": self.results["Wall-clock Time"], + } + self.results["names"] = { + "experiment_design": self.results["Experiment Design Names"], + "experiment_output": self.results["Experiment Output Names"], + "unknown_parameter": self.results["Unknown Parameter Names"], + "measurement_error": self.results["Measurement Error Names"], + } + self.results["diagnostics"] = { + "symmetry_breaking": symmetry_breaking_info, + "warnings": diagnostics_warnings, + "lhs_initialization": lhs_init_diagnostics, + } + self.results["scenarios"] = scenarios_structured + + # Save results to file if requested + if results_file is not None: + with open(results_file, "w") as file: + json.dump(self.results, file, indent=2, cls=_DoEResultsJSONEncoder) + + # 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))) + else: # minimum_eigenvalue, condition_number, 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``) a larger value is better. For minimisation + objectives (``trace`` / A-optimality) a smaller value is better. + LHS initialization is not supported for ``minimum_eigenvalue`` or + ``condition_number``; those return 0.0. + """ + 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, + lhs_parallel=False, + lhs_combo_parallel=False, + lhs_n_workers=None, + lhs_combo_chunk_size=5000, + lhs_combo_parallel_threshold=20000, + lhs_max_wall_clock_time=None, + ): + """ + Use per-dimension Latin Hypercube Sampling to identify a good initial + experiment design for ``optimize_experiments``. + """ + # Use a monotonic wall-clock for progress estimates and deadline checks + # in both serial and threaded LHS evaluation loops. + t_start = time.perf_counter() + + # 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='lhs'``." + ) + + 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 = 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)) + t_after_sampling = time.perf_counter() + 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." + ) + + resolved_workers = ( + lhs_n_workers + if lhs_n_workers is not None + else max(1, min(os.cpu_count() or 1, 8)) + ) + use_parallel_fim = lhs_parallel and resolved_workers > 1 + fim_eval_mode = "parallel" if use_parallel_fim else "serial" + self.logger.info( + f"LHS: evaluating FIM at {n_candidates} candidate designs " + f"({n_candidates * n_scenarios_per_candidate} solver calls expected; " + f"{fim_eval_mode} candidate FIM mode)." + ) + + # Track worker DoE construction failures to avoid repeated logging of the same issue + _solver_fallback_lock = threading.Lock() + _solver_fallback_logged = False + + def _make_worker_solver(): + nonlocal _solver_fallback_logged + solver_name = getattr(self.solver, "name", None) + if solver_name is None: + with _solver_fallback_lock: + if not _solver_fallback_logged: + self.logger.debug( + "LHS parallel: solver has no 'name' attribute; worker DoE " + "will use default solver settings." + ) + _solver_fallback_logged = True + return None + worker_solver = pyo.SolverFactory(solver_name) + if worker_solver is None: + with _solver_fallback_lock: + if not _solver_fallback_logged: + self.logger.debug( + f"LHS parallel: could not construct solver '{solver_name}'; " + "worker DoE will use default solver settings." + ) + _solver_fallback_logged = True + return None + try: + if hasattr(self.solver, "options") and hasattr( + worker_solver, "options" + ): + worker_solver.options.update(self.solver.options) + except Exception: + pass + return worker_solver + + thread_state = threading.local() + n_params = len(first_exp_block.fd_scenario_blocks[0].unknown_parameters) + + def _compute_candidate_fim(idx_pt): + idx, pt = idx_pt + try: + worker_doe = getattr(thread_state, "doe", None) + if worker_doe is None: + # Retry worker DoE construction on subsequent points even if a + # previous point failed; transient failures should not disable + # the rest of this thread's candidate evaluations. + worker_solver = _make_worker_solver() + worker_doe = DesignOfExperiments( + experiment=self.experiment_list, + fd_formula=self.fd_formula.value, + step=self.step, + objective_option=self.objective_option.value, + use_grey_box_objective=self.use_grey_box, + scale_constant_value=self.scale_constant_value, + scale_nominal_param_value=self.scale_nominal_param_value, + prior_FIM=self.prior_FIM, + jac_initial=self.jac_initial, + fim_initial=self.fim_initial, + L_diagonal_lower_bound=self.L_diagonal_lower_bound, + solver=worker_solver, + grey_box_solver=self.grey_box_solver, + tee=False, + grey_box_tee=False, + get_labeled_model_args=self.get_labeled_model_args, + logger_level=logging.ERROR, + improve_cholesky_roundoff_error=self.improve_cholesky_roundoff_error, + _Cholesky_option=self.Cholesky_option, + _only_compute_fim_lower=self.only_compute_fim_lower, + ) + thread_state.doe = worker_doe + # LHS initialization evaluates candidate points against the + # canonical experiment template (experiment_list[0]). + fim = worker_doe._compute_fim_at_point_no_prior( + experiment_index=0, input_values=list(pt) + ) + return idx, fim + except Exception as exc: + if not getattr(thread_state, "construction_failed", False): + thread_state.construction_failed = True + self.logger.error( + f"LHS: worker DoE construction/evaluation failed on thread " + f"{threading.current_thread().name}: {exc}. " + "Using zero FIM for this candidate and continuing." + ) + return idx, np.zeros((n_params, n_params)) + + # 3. Evaluate FIM at every candidate design + candidate_fims = [None] * n_candidates + if lhs_parallel and not use_parallel_fim: + self.logger.warning( + "LHS candidate FIM parallel evaluation requested with " + "``lhs_parallel=True``, but running serially: " + f"resolved_workers={resolved_workers} <= 1." + ) + timed_out = False + deadline = ( + None + if lhs_max_wall_clock_time is None + else (t_start + lhs_max_wall_clock_time) + ) + if use_parallel_fim: + self.logger.info( + f"LHS: using parallel candidate FIM evaluation with {resolved_workers} workers." + ) + idx_iter = iter(enumerate(candidate_points)) + max_pending = max(1, 2 * resolved_workers) + with _cf.ThreadPoolExecutor(max_workers=resolved_workers) as ex: + pending = set() + n_done = 0 + while True: + while len(pending) < max_pending: + if deadline is not None and time.perf_counter() > deadline: + timed_out = True + break + try: + idx_pt = next(idx_iter) + except StopIteration: + break + pending.add(ex.submit(_compute_candidate_fim, idx_pt)) + + if not pending: + break + + done_now, pending = _cf.wait( + pending, timeout=0.1, return_when=_cf.FIRST_COMPLETED + ) + for fut in done_now: + pt_idx, fim = fut.result() + candidate_fims[pt_idx] = fim + n_done += 1 + if n_done % max(1, n_candidates // 10) == 0: + elapsed = time.perf_counter() - t_start + frac_done = n_done / n_candidates + est_total = elapsed / frac_done if frac_done > 0 else 0 + self.logger.info( + f" LHS FIM eval: {n_done}/{n_candidates} " + f"({elapsed:.1f}s elapsed, ~{est_total:.1f}s total)" + ) + + if timed_out: + for fut in pending: + fut.cancel() + break + else: + for pt_idx, pt in enumerate(candidate_points): + if deadline is not None and time.perf_counter() > deadline: + timed_out = True + break + 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: + elapsed = time.perf_counter() - t_start + frac_done = (pt_idx + 1) / n_candidates + est_total = elapsed / frac_done if frac_done > 0 else 0 + self.logger.info( + f" LHS FIM eval: {pt_idx + 1}/{n_candidates} " + f"({elapsed:.1f}s elapsed, ~{est_total:.1f}s total)" + ) + + 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. + if not computed_pairs: + timed_out = True + computed_pairs = [(candidate_points[0], np.zeros((n_params, n_params)))] + + # If timeout stops FIM evaluation early, retain only scored candidates so + # downstream combination scoring does not use missing entries. + if len(computed_pairs) < n_candidates: + self.logger.warning( + "LHS candidate FIM evaluation reached time budget " + f"({lhs_max_wall_clock_time}s). Scored {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) + + t_after_fim = time.perf_counter() + fim_eval_time = t_after_fim - t_after_sampling + self.logger.info(f"LHS: completed FIM evaluations in {fim_eval_time:.1f}s.") + + # 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``." + ) + + 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 + + def _score_chunk(combo_chunk, deadline_ts): + local_best_obj = -np.inf if is_maximize else np.inf + local_best_combo = None + local_timed_out = False + for combo in combo_chunk: + if deadline_ts is not None and time.perf_counter() > deadline_ts: + local_timed_out = True + break + 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 > local_best_obj: + local_best_obj = obj_val + local_best_combo = combo + else: + if obj_val < local_best_obj: + local_best_obj = obj_val + local_best_combo = combo + return local_best_obj, local_best_combo, local_timed_out + + use_parallel_combo = ( + lhs_combo_parallel + and n_combinations >= lhs_combo_parallel_threshold + and resolved_workers > 1 + ) + if lhs_combo_parallel and not use_parallel_combo: + reasons = [] + if n_combinations < lhs_combo_parallel_threshold: + reasons.append( + f"n_combinations={n_combinations} < " + f"lhs_combo_parallel_threshold={lhs_combo_parallel_threshold}" + ) + if resolved_workers <= 1: + reasons.append(f"resolved_workers={resolved_workers} <= 1") + reason_txt = ( + "; ".join(reasons) if reasons else "parallel preconditions not met" + ) + self.logger.warning( + "LHS combination scoring requested with " + "``lhs_combo_parallel=True``, but running serially: " + f"{reason_txt}." + ) + + if use_parallel_combo: + self.logger.info( + f"LHS: using parallel combination scoring with {resolved_workers} workers " + f"(chunk_size={lhs_combo_chunk_size})." + ) + combo_iter = _combinations(range(n_candidates), n_exp) + max_pending = max(1, 2 * resolved_workers) + with _cf.ThreadPoolExecutor(max_workers=resolved_workers) as ex: + pending = set() + while True: + while len(pending) < max_pending: + if deadline is not None and time.perf_counter() > deadline: + timed_out = True + break + chunk = list(_islice(combo_iter, lhs_combo_chunk_size)) + if not chunk: + break + pending.add(ex.submit(_score_chunk, chunk, deadline)) + if not pending: + break + + done, pending = _cf.wait(pending, return_when=_cf.FIRST_COMPLETED) + for fut in done: + local_obj, local_combo, local_timed_out = fut.result() + timed_out = timed_out or local_timed_out + if local_combo is None: + continue + if is_maximize: + if local_obj > best_obj: + best_obj = local_obj + best_combo = local_combo + else: + if local_obj < best_obj: + best_obj = local_obj + best_combo = local_combo + + if deadline is not None and time.perf_counter() > deadline: + timed_out = True + for fut in pending: + fut.cancel() + break + else: + for combo in _combinations(range(n_candidates), n_exp): + if deadline is not None and time.perf_counter() > deadline: + timed_out = True + break + 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 timed_out: + self.logger.warning( + f"LHS combination scoring reached time budget " + f"({lhs_max_wall_clock_time}s). Returning best-so-far." + ) + + t_after_combo = time.perf_counter() + combo_time = t_after_combo - t_after_fim + + 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": "thread" if use_parallel_fim else "serial", + "combo_mode": "thread" if use_parallel_combo else "serial", + "n_workers": resolved_workers, + "n_candidates": n_candidates, + "n_combinations": n_combinations, + "elapsed_sampling_s": t_after_sampling - t_start, + "elapsed_fim_eval_s": fim_eval_time, + "elapsed_combo_scoring_s": combo_time, + "elapsed_total_s": t_after_combo - t_start, + "timed_out": timed_out, + "time_budget_s": lhs_max_wall_clock_time, + "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): raise NotImplementedError("Multiple experiment optimization not yet supported.") @@ -558,14 +2202,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. @@ -595,20 +2253,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 @@ -622,9 +2350,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 @@ -783,9 +2513,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 @@ -864,7 +2596,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 @@ -878,6 +2612,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: @@ -905,25 +2642,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() ] ) @@ -1005,9 +2744,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 @@ -1059,12 +2799,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: @@ -1117,20 +2859,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 @@ -1151,7 +2911,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 @@ -1165,15 +2925,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) @@ -1259,8 +3022,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 @@ -1271,7 +3035,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: @@ -1288,9 +3052,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: @@ -1302,11 +3066,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): @@ -1317,8 +3085,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( @@ -1372,10 +3140,11 @@ def create_objective_function(self, model=None): model.obj_cons = pyo.Block() # Assemble the FIM matrix. This is helpful for initialization! + # collect current FIM values in row-major order fim_vals = [ model.fim[bu, un].value - for i, bu in enumerate(model.parameter_names) - for j, un in enumerate(model.parameter_names) + for bu in model.parameter_names + for un in model.parameter_names ] fim = np.array(fim_vals).reshape( len(model.parameter_names), len(model.parameter_names) @@ -1411,23 +3180,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: @@ -1522,6 +3278,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 = [] @@ -1546,7 +3303,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), @@ -1565,17 +3326,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 @@ -1612,6 +3373,283 @@ def determinant_general(b): # add dummy objective function model.objective = pyo.Objective(expr=0) + 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.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 self.objective_option == ObjectiveLib.trace and not self.Cholesky_option: + raise ValueError( + "objective_option='trace' currently only implemented with " + "``_Cholesky_option=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. Create obj_cons block to hold objective-related constraints + # (similar to single-experiment case in create_objective_function) + scenario.obj_cons = pyo.Block() + # 5. Create variables and constraints (initialization will happen after square solve) + if ( + self.Cholesky_option + and self.objective_option == ObjectiveLib.determinant + ): + # 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: + # 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: + # 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: + # 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.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): # Add external grey box block to a block named ``obj_cons`` to # reuse material for initializing the objective-free square model @@ -1699,30 +3737,71 @@ def check_model_labels(self, model=None): "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".' ) + # 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".' ) + # 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".' ) + # 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.") # Check the FIM shape against what is expected from the model. @@ -1868,9 +3947,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 @@ -2643,7 +4724,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 " @@ -2652,7 +4733,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()] @@ -2680,7 +4761,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 " @@ -2689,7 +4770,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()] @@ -2716,7 +4797,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 " @@ -2725,7 +4806,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()] @@ -2754,7 +4835,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 " @@ -2763,7 +4844,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/tests/experiment_class_example_flags.py b/pyomo/contrib/doe/tests/experiment_class_example_flags.py index c6b3cb9e0ec..33e4cf65010 100644 --- a/pyomo/contrib/doe/tests/experiment_class_example_flags.py +++ b/pyomo/contrib/doe/tests/experiment_class_example_flags.py @@ -70,3 +70,77 @@ 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) + + 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 7ba299dc4d0..5100c3d7401 100644 --- a/pyomo/contrib/doe/tests/test_doe_build.py +++ b/pyomo/contrib/doe/tests/test_doe_build.py @@ -7,7 +7,11 @@ # 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 unittest.mock import patch from pyomo.common.dependencies import ( numpy as np, @@ -19,6 +23,7 @@ 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") @@ -28,6 +33,9 @@ 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, ) @@ -97,9 +105,11 @@ 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.scenario_blocks[0].unknown_parameters.items()]] + [[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)) @@ -123,7 +133,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 @@ -173,16 +183,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 @@ -212,17 +222,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 @@ -250,17 +264,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 @@ -282,7 +300,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_" @@ -315,7 +335,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_" @@ -348,7 +370,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_" @@ -502,11 +526,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) + "]") ) @@ -569,7 +593,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), @@ -583,7 +607,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), @@ -640,5 +664,400 @@ def test_trace_initialization_consistency(self): self.assertAlmostEqual(val, expected, places=4) +class TestOptimizeExperimentsBuildStructure(unittest.TestCase): + """Coverage for optimize_experiments() build, output, and diagnostics behavior.""" + + def _make_solver(self): + # Make solver object with options for DoE runs. + solver = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + return solver + + def test_constructor_accepts_single_experiment_or_list(self): + # Tests that the public ``experiment`` argument accepts either one + # experiment object or a list of experiment objects. + 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) + + @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 = self._make_solver() + init_solver = self._make_solver() + # Use distinct option values so each solver path can be identified. + main_solver.options["max_iter"] = 321 + init_solver.options["max_iter"] = 123 + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=main_solver, + ) + + # Track both how many times each solver is called and the chronological + # order of those calls. The optimize_experiments() contract is that all + # setup solves run first on init_solver and the final NLP solve runs last + # on the primary solver. + main_calls = 0 + init_calls = 0 + call_order = [] + # Record an option value on each solve so the test can verify that the + # call really went through the expected solver object, not just the + # expected phase label. + option_markers = [] + original_main_solve = main_solver.solve + original_init_solve = init_solver.solve + + def _main_solve(*args, **kwargs): + nonlocal main_calls + main_calls += 1 + call_order.append("main") + option_markers.append(main_solver.options.get("max_iter")) + return original_main_solve(*args, **kwargs) + + def _init_solve(*args, **kwargs): + nonlocal init_calls + init_calls += 1 + call_order.append("init") + option_markers.append(init_solver.options.get("max_iter")) + return original_init_solve(*args, **kwargs) + + # Patch both solver objects in place so the real solves still run while + # we collect lightweight diagnostics about solver routing. + with ( + patch.object(main_solver, "solve", side_effect=_main_solve), + patch.object(init_solver, "solve", side_effect=_init_solve), + ): + 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_calls, 1) # At least one initialization solve + self.assertEqual(main_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. + 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 patched solve() calls above. + self.assertEqual( + doe_obj.results["settings"]["initialization"]["solver_name"], + getattr(init_solver, "name", str(init_solver)), + ) + self.assertEqual( + doe_obj.results["run_info"]["solver"]["name"], + 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) + + scenario_result = doe_obj.results["Scenarios"][0] + self.assertEqual(doe_obj.results["Initialization Method"], "none") + self.assertEqual(doe_obj.results["Number of Experiments per Scenario"], 2) + self.assertEqual(len(scenario_result["Experiments"]), 2) + self.assertEqual(len(doe_obj.results["Experiment Design Names"]), 1) + self.assertEqual(len(doe_obj.results["Unknown Parameter Names"]), 2) + + # New structured payload checks + self.assertIn("run_info", doe_obj.results) + self.assertIn("settings", doe_obj.results) + self.assertIn("timing", doe_obj.results) + self.assertIn("names", doe_obj.results) + self.assertIn("scenarios", doe_obj.results) + self.assertEqual( + doe_obj.results["run_info"]["solver"]["status"], + doe_obj.results["Solver Status"], + ) + self.assertEqual( + doe_obj.results["settings"]["modeling"]["n_experiments_per_scenario"], 2 + ) + self.assertFalse(doe_obj.results["settings"]["modeling"]["template_mode"]) + self.assertEqual(len(doe_obj.results["scenarios"]), 1) + self.assertEqual(doe_obj.results["scenarios"][0]["id"], 0) + self.assertEqual(len(doe_obj.results["scenarios"][0]["experiments"]), 2) + self.assertEqual(doe_obj.results["scenarios"][0]["experiments"][0]["id"], 0) + + # hour of exp[0] should be <= hour of exp[1] due to symmetry breaking + h0 = scenario_result["Experiments"][0]["Experiment Design"][0] + h1 = scenario_result["Experiments"][1]["Experiment 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["settings"]["modeling"]["template_mode"]) + self.assertIn("Scenarios", payload) + self.assertIn("run_info", payload) + self.assertIn("settings", payload) + self.assertIn("timing", payload) + self.assertIn("names", payload) + self.assertIn("scenarios", 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["settings"]["modeling"]["template_mode"]) + + @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["Number of Experiments per Scenario"], 1) + self.assertTrue(doe_obj.results["settings"]["modeling"]["template_mode"]) + + @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="lhs", init_n_samples=2, init_seed=11 + ) + + timing = doe_obj.results["timing"] + self.assertIn("lhs_initialization_s", timing) + self.assertGreaterEqual(timing["lhs_initialization_s"], 0.0) + self.assertAlmostEqual( + timing["total_s"], + timing["build_s"] + + timing["lhs_initialization_s"] + + timing["initialization_s"] + + timing["solve_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 threaded 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="lhs", + init_n_samples=2, + init_seed=11, + init_parallel=True, + init_combo_parallel=True, + init_n_workers=2, + init_combo_chunk_size=2, + init_combo_parallel_threshold=1, + init_max_wall_clock_time=60.0, + ) + lhs_diag = doe_obj.results["diagnostics"]["lhs_initialization"] + self.assertIsNotNone(lhs_diag) + self.assertEqual(lhs_diag["candidate_fim_mode"], "thread") + self.assertEqual(lhs_diag["combo_mode"], "thread") + self.assertEqual(lhs_diag["n_workers"], 2) + self.assertFalse(lhs_diag["timed_out"]) + self.assertGreater(lhs_diag["elapsed_total_s"], 0.0) + self.assertIn("best_obj", lhs_diag) + self.assertIsInstance(lhs_diag["best_obj"], float) + self.assertTrue(np.isfinite(lhs_diag["best_obj"])) + self.assertGreater(lhs_diag["best_obj"], 0.0) + self.assertIn("best_obj_log10", lhs_diag) + self.assertIsInstance(lhs_diag["best_obj_log10"], float) + self.assertAlmostEqual( + lhs_diag["best_obj_log10"], np.log10(lhs_diag["best_obj"]), places=12 + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_termination_message_bytes(self): + # Tests that solver termination messages returned as bytes are decoded + # and persisted as plain strings in results. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + original_solve = doe_obj.solver.solve + + def _solve_with_bytes_message(*args, **kwargs): + res = original_solve(*args, **kwargs) + res.solver.message = b"byte-message" + return res + + with patch.object( + doe_obj.solver, "solve", side_effect=_solve_with_bytes_message + ): + doe_obj.optimize_experiments(n_exp=1) + + self.assertEqual(doe_obj.results["Termination Message"], "byte-message") + + @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") + def test_optimize_experiments_termination_message_fallback_to_str(self): + # Tests that non-string termination messages fall back to str(message) so + # results always store a serializable user-facing value. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + original_solve = doe_obj.solver.solve + + class _CustomMessage: + def __str__(self): + return "custom-message" + + def _solve_with_custom_message(*args, **kwargs): + res = original_solve(*args, **kwargs) + res.solver.message = _CustomMessage() + return res + + with patch.object( + doe_obj.solver, "solve", side_effect=_solve_with_custom_message + ): + doe_obj.optimize_experiments(n_exp=1) + + self.assertEqual(doe_obj.results["Termination Message"], "custom-message") + + def test_maximize_objective_set_contents(self): + 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_doe_results_json_encoder_handles_numpy_and_enum(self): + 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)) + + def test_symmetrize_lower_tri_helper(self): + 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 10165aa6267..1252bf8f255 100644 --- a/pyomo/contrib/doe/tests/test_doe_errors.py +++ b/pyomo/contrib/doe/tests/test_doe_errors.py @@ -6,7 +6,8 @@ # 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, @@ -17,26 +18,48 @@ from pyomo.common.errors import DeveloperError import pyomo.common.unittest as unittest +from unittest.mock import patch 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 +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.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.parmest.examples.rooney_biegler.rooney_biegler import ( RooneyBieglerExperiment, ) 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() +class _DummyExperiment: + def get_labeled_model(self, **kwargs): + raise RuntimeError("Should not be called in argument-validation tests") + + def get_rooney_biegler_experiment_flag(): """Get a fresh RooneyBieglerExperimentFlag instance for testing. @@ -53,9 +76,17 @@ def get_rooney_biegler_experiment_flag(): ) +def _make_ipopt_solver(): + solver = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + return solver + + 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 @@ -90,13 +121,23 @@ def test_experiment_none_error(self): 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' parameter is required" ): # Experiment provided as None DoE_args = get_standard_args(None, fd_method, obj_used, flag_val) doe_obj = DesignOfExperiments(**DoE_args) + def test_experiment_empty_list_error(self): + with self.assertRaisesRegex( + ValueError, "The 'experiment' list cannot be empty" + ): + 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" @@ -105,8 +146,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) @@ -677,7 +717,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): @@ -762,6 +802,36 @@ 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) + + def _fake_sequential(*args, **kwargs): + # This is only used if execution reaches the FIM solve call. + doe_obj.seq_FIM = np.eye(2) + + with patch.object(doe_obj, "_sequential_FIM", side_effect=_fake_sequential): + # The mismatch is detected before the second experiment solve, + # when compute_FIM validates unknown parameter values. + 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" @@ -781,6 +851,407 @@ def test_invalid_trace_without_cholesky(self): ): doe_obj.create_objective_function() + def test_optimize_experiments_invalid_init_method(self): + # Tests that an unsupported init_method value is rejected. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_method`` must be one of \[None, 'lhs'\], got 'bad'." + ): + doe_obj.optimize_experiments(init_method="bad") + + def test_optimize_experiments_init_method_enum_accepted(self): + # Tests that enum init_method inputs are accepted and proceed to later validation. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_n_samples`` must be a positive integer, got 0." + ): + doe_obj.optimize_experiments( + init_method=InitializationMethod.lhs, init_n_samples=0 + ) + + def test_optimize_experiments_init_method_enum_invalid_init_n_samples(self): + # Tests that enum init_method path still enforces init_n_samples validation. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_n_samples`` must be a positive integer, got 0." + ): + doe_obj.optimize_experiments( + init_method=InitializationMethod.lhs, init_n_samples=0 + ) + + def test_optimize_experiments_invalid_init_n_samples(self): + # Tests that non-positive init_n_samples is rejected for LHS initialization. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_n_samples`` must be a positive integer, got 0." + ): + doe_obj.optimize_experiments(init_method="lhs", init_n_samples=0) + + def test_optimize_experiments_invalid_init_n_samples_float(self): + # Tests that non-integer init_n_samples values are rejected. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_n_samples`` must be a positive integer, got 2.5." + ): + doe_obj.optimize_experiments(init_method="lhs", init_n_samples=2.5) + + 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='lhs'`` is currently supported only in template mode", + ): + doe_obj.optimize_experiments(init_method="lhs") + + 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" + ) + with patch("pyomo.contrib.doe.doe.scipy_available", False): + with self.assertRaisesRegex( + ImportError, r"LHS initialization requires scipy" + ): + doe_obj.optimize_experiments(init_method="lhs") + + def test_optimize_experiments_init_parallel_requires_bool(self): + # Tests that init_parallel must be a boolean flag. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_parallel`` must be a bool, got 1." + ): + doe_obj.optimize_experiments(init_method="lhs", init_parallel=1) + + def test_optimize_experiments_init_combo_parallel_requires_bool(self): + # Tests that init_combo_parallel must be a boolean flag. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_combo_parallel`` must be a bool" + ): + doe_obj.optimize_experiments(init_method="lhs", init_combo_parallel="yes") + + def test_optimize_experiments_init_n_workers_must_be_positive_integer(self): + # Tests that init_n_workers must be None or a positive integer. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_n_workers`` must be None or a positive integer" + ): + doe_obj.optimize_experiments(init_method="lhs", init_n_workers=0) + + def test_optimize_experiments_init_combo_chunk_size_must_be_positive_integer(self): + # Tests that init_combo_chunk_size must be a positive integer. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_combo_chunk_size`` must be a positive integer" + ): + doe_obj.optimize_experiments(init_method="lhs", init_combo_chunk_size=0) + + def test_optimize_experiments_init_combo_parallel_threshold_positive_integer(self): + # Tests that init_combo_parallel_threshold must be a positive integer. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_combo_parallel_threshold`` must be a positive integer" + ): + doe_obj.optimize_experiments( + init_method="lhs", init_combo_parallel_threshold=0 + ) + + def test_optimize_experiments_init_max_wall_clock_time_must_be_positive(self): + # Tests that init_max_wall_clock_time must be positive when provided. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, + r"``init_max_wall_clock_time`` must be None or a positive number", + ): + doe_obj.optimize_experiments(init_method="lhs", init_max_wall_clock_time=0) + + def test_optimize_experiments_init_max_wall_clock_time_rejects_nan(self): + # Tests that NaN is rejected for init_max_wall_clock_time. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, + r"``init_max_wall_clock_time`` must be None or a positive number", + ): + doe_obj.optimize_experiments( + init_method="lhs", init_max_wall_clock_time=float("nan") + ) + + def test_optimize_experiments_init_max_wall_clock_time_rejects_inf(self): + # Tests that infinity is rejected for init_max_wall_clock_time. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, + r"``init_max_wall_clock_time`` must be None or a positive number", + ): + doe_obj.optimize_experiments( + init_method="lhs", init_max_wall_clock_time=float("inf") + ) + + def test_optimize_experiments_n_exp_with_multi_list(self): + # Tests that n_exp cannot be set when multiple experiments are already provided. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment(), _DummyExperiment()], + objective_option="pseudo_trace", + ) + with self.assertRaisesRegex( + ValueError, + r"``n_exp`` must not be set when the experiment list contains more than one experiment", + ): + doe_obj.optimize_experiments(n_exp=2) + + def test_optimize_experiments_n_exp_not_positive(self): + # Tests that n_exp must be a positive integer. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``n_exp`` must be a positive integer, got 0." + ): + doe_obj.optimize_experiments(n_exp=0) + + def test_optimize_experiments_results_file_bad_type(self): + # Tests that results_file accepts only str or pathlib.Path values. + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``results_file`` must be either a Path object or a string." + ): + doe_obj.optimize_experiments(results_file=5) + + def test_optimize_experiments_init_solver_bad_type(self): + # Tests that init_solver must be None or a solver object with solve(). + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, + r"``init_solver`` must be None or a solver object with a 'solve' method.", + ): + doe_obj.optimize_experiments(init_solver=object()) + + def test_optimize_experiments_init_seed_requires_integer(self): + # Tests that init_seed must be an integer (or None). + doe_obj = DesignOfExperiments( + experiment=[_DummyExperiment()], objective_option="pseudo_trace" + ) + with self.assertRaisesRegex( + ValueError, r"``init_seed`` must be None or an integer" + ): + doe_obj.optimize_experiments( + n_exp=2, init_method="lhs", init_n_samples=2, init_seed=1.5 + ) + + +@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() + + 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) + + with patch( + "pyomo.contrib.doe.doe.pyo.ComponentUID.find_component_on", + autospec=True, + side_effect=_fail_only_symmetry_mapping, + ): + with self.assertRaisesRegex( + RuntimeError, "Failed to map symmetry breaking variable" + ): + doe_obj.optimize_experiments(n_exp=2) + + 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") + with patch( + "pyomo.contrib.doe.doe._combinations", return_value=iter([(0, 1)]) + ): + with patch.object( + doe_obj, + "_compute_fim_at_point_no_prior", + return_value=np.eye(2), + ): + doe_obj.optimize_experiments( + n_exp=2, + init_method="lhs", + init_n_samples=10001, + init_seed=11, + ) + + 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_combo_parallel_requested_but_not_used_warns(self): + # Tests that combo-parallel requests warn when thresholds force serial scoring. + doe_obj = DesignOfExperiments( + experiment=[RooneyBieglerMultiExperiment(hour=2.0, y=10.0)], + objective_option="pseudo_trace", + step=1e-2, + solver=self._make_solver(), + ) + with patch.object( + doe_obj, "_compute_fim_at_point_no_prior", return_value=np.eye(2) + ): + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as cm: + doe_obj.optimize_experiments( + n_exp=2, + init_method="lhs", + init_n_samples=2, + init_seed=11, + init_combo_parallel=True, + init_n_workers=2, + init_combo_parallel_threshold=10_000, + ) + + self.assertTrue( + any( + "lhs_combo_parallel=True" in msg and "running serially" in msg + for msg in 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="lhs", init_n_samples=2) + if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/doe/tests/test_doe_solve.py b/pyomo/contrib/doe/tests/test_doe_solve.py index 80a77b5f012..1435ec9bc9c 100644 --- a/pyomo/contrib/doe/tests/test_doe_solve.py +++ b/pyomo/contrib/doe/tests/test_doe_solve.py @@ -9,7 +9,11 @@ import json import logging import os, os.path +import subprocess +import time +from itertools import combinations, product from glob import glob +from unittest.mock import patch from pyomo.common.dependencies import ( numpy as np, @@ -40,6 +44,7 @@ ) from pyomo.contrib.doe.tests.experiment_class_example_flags import ( RooneyBieglerExperimentBad, + RooneyBieglerMultiExperiment, ) from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( RooneyBieglerExperiment, @@ -54,6 +59,38 @@ ipopt_available = SolverFactory("ipopt").available() k_aug_available = SolverFactory("k_aug", solver_io="nl", validate=False) + +def k_aug_runtime_available(): + """ + Check that k_aug is not only discoverable but also runnable in this + environment (e.g., no missing dynamic libraries). + """ + if not k_aug_available.available(False): + return False + exe = k_aug_available.executable() + if not exe: + return False + + try: + # Trigger dynamic loader checks; return code may be nonzero for usage, + # so we inspect output for runtime-linker failures. + proc = subprocess.run( + [exe, "--help"], + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True, + check=False, + ) + except OSError: + return False + + output = (proc.stdout or "") + (proc.stderr or "") + bad_runtime_markers = ("Library not loaded", "dyld:", "libgfortran") + if any(marker in output for marker in bad_runtime_markers): + return False + return True + + currdir = this_file_dir() file_path = os.path.join(currdir, "..", "examples", "result.json") @@ -113,10 +150,10 @@ 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() ] param_vals = np.array( - [[v for k, v in model.scenario_blocks[0].unknown_parameters.items()]] + [[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)) @@ -139,7 +176,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 @@ -383,12 +420,58 @@ 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" + + fim_exp1 = np.array([[12.0, 3.0], [3.0, 8.0]]) + fim_exp2 = np.array([[5.0, 2.0], [2.0, 4.0]]) + 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) + + def _fake_sequential_fim(*args, **kwargs): + # Return deterministic per-experiment FIMs keyed by the fixed + # experiment input so the aggregation can be asserted exactly. + model = kwargs.get("model") + hour = float(pyo.value(model.hour)) + if np.isclose(hour, 1.5): + doe_multi.seq_FIM = fim_exp1.copy() + elif np.isclose(hour, 3.5): + doe_multi.seq_FIM = fim_exp2.copy() + else: + raise RuntimeError(f"Unexpected hour value in mocked test: {hour}") + + with patch.object( + doe_multi, "_sequential_FIM", side_effect=_fake_sequential_fim + ): + fim_total = doe_multi.compute_FIM(method="sequential") + + fim_expected = fim_exp1 + fim_exp2 + prior_fim + self.assertTrue(np.allclose(fim_total, fim_expected, atol=1e-12)) + self.assertEqual(len(doe_multi._computed_FIM_by_experiment), 2) + self.assertTrue( + np.allclose(doe_multi._computed_FIM_by_experiment[0], fim_exp1, atol=1e-12) + ) + self.assertTrue( + np.allclose(doe_multi._computed_FIM_by_experiment[1], fim_exp2, atol=1e-12) + ) + # 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_runtime_available(), + "The 'k_aug' command is not available or not runnable in this environment", ) def test_compute_FIM_kaug(self): fd_method = "forward" @@ -486,7 +569,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() ] @@ -919,5 +1002,1250 @@ 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 = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + 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, + ) + + 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 + + with patch.object(doe, "_sequential_FIM", side_effect=RuntimeError("boom")): + 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] + ) + + 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): + # Tests that optimization proceeds through the Cholesky jitter branch when needed. + doe = self._make_template_doe("determinant") + + original_solve = doe.solver.solve + solve_count = {"n": 0} + + class _MockSolverInfo: + status = "ok" + termination_condition = "optimal" + message = "mock-solve" + + class _MockResults: + solver = _MockSolverInfo() + + def _solve_first_real_then_mock(*args, **kwargs): + solve_count["n"] += 1 + if solve_count["n"] == 1: + return original_solve(*args, **kwargs) + return _MockResults() + + with patch( + "pyomo.contrib.doe.doe.np.linalg.eigvals", + return_value=np.array([-1.0, 1.0]), + ) as eigvals_mock: + with patch.object( + doe.solver, "solve", side_effect=_solve_first_real_then_mock + ): + doe.optimize_experiments(n_exp=1) + + self.assertEqual(doe.results["Solver Status"], "ok") + self.assertGreaterEqual(eigvals_mock.call_count, 1) + + def test_optimize_experiments_lhs_matches_bruteforce_combo(self): + # Tests that LHS-selected initial points match an independent brute-force oracle. + n_exp = 2 + lhs_n_samples = 3 + 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] + + # Run full optimization with LHS initialization and compare chosen points. + doe = self._make_template_doe("pseudo_trace") + doe.optimize_experiments( + n_exp=n_exp, + init_method="lhs", + init_n_samples=lhs_n_samples, + init_seed=lhs_seed, + ) + + actual_points = doe.results["LHS Best Initial Points"] + 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.assertEqual(doe.results["Initialization Method"], "lhs") + + # Numerical consistency of aggregated FIM in result payload. + scenario = doe.results["Scenarios"][0] + total_fim = np.array(scenario["Total FIM"]) + prior = np.array(doe.results["Prior FIM"]) + exp_fim_sum = sum( + (np.array(exp_data["FIM"]) for exp_data in scenario["Experiments"]), + np.zeros_like(total_fim), + ) + self.assertTrue(np.allclose(total_fim, exp_fim_sum + prior, atol=1e-6)) + + 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 = doe.results["Scenarios"][0]["Experiments"][0][ + "Experiment Design" + ] + first_build_time = doe.results["timing"]["build_s"] + + doe.optimize_experiments(n_exp=1) + second_design = doe.results["Scenarios"][0]["Experiments"][0][ + "Experiment Design" + ] + + self.assertEqual(len(first_design), len(second_design)) + self.assertIn("timing", doe.results) + self.assertGreater(doe.results["timing"]["build_s"], 0.0) + self.assertGreaterEqual(first_build_time, 0.0) + self.assertEqual(len(list(doe.model.param_scenario_blocks.keys())), 1) + + def test_lhs_combo_parallel_matches_serial(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + # Use deterministic synthetic FIMs to isolate combo scorer behavior. + 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]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_serial, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, lhs_seed=123, n_exp=2, lhs_combo_parallel=False + ) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_parallel, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=123, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=2, + lhs_combo_parallel_threshold=1, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_norm) + + def test_lhs_parallel_fim_eval_matches_serial(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + # Patch at class level so both serial path (self) and parallel worker + # DOE objects use the same deterministic synthetic FIM mapping. + def _fake_fim(self_obj, experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 0.5, 0.0], [0.0, 3.0 * x + 0.5]]) + + with patch.object( + DesignOfExperiments, + "_compute_fim_at_point_no_prior", + autospec=True, + side_effect=_fake_fim, + ): + points_serial, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, lhs_seed=321, n_exp=2, lhs_parallel=False + ) + + points_parallel, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=321, + n_exp=2, + lhs_parallel=True, + lhs_n_workers=2, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_norm) + + def test_lhs_parallel_fim_eval_real_path_smoke(self): + doe_serial = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe_serial, n_exp=2) + points_serial, _ = doe_serial._lhs_initialize_experiments( + lhs_n_samples=3, + lhs_seed=9, + n_exp=2, + lhs_parallel=False, + lhs_combo_parallel=False, + ) + + doe_parallel = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe_parallel, n_exp=2) + points_parallel, _ = doe_parallel._lhs_initialize_experiments( + lhs_n_samples=3, + lhs_seed=9, + n_exp=2, + lhs_parallel=True, + lhs_n_workers=2, + lhs_combo_parallel=False, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_norm) + + def test_lhs_parallel_solver_without_name_uses_default_worker_solver(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + doe.solver = object() + + def _fake_fim(self_obj, experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 0.5, 0.0], [0.0, x + 1.0]]) + + with patch.object( + DesignOfExperiments, + "_compute_fim_at_point_no_prior", + autospec=True, + side_effect=_fake_fim, + ): + with self.assertLogs("pyomo.contrib.doe.doe", level="DEBUG") as log_cm: + points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=14, + n_exp=2, + lhs_parallel=True, + lhs_n_workers=2, + lhs_combo_parallel=False, + ) + + self.assertEqual(len(points), 2) + self.assertTrue( + any("solver has no 'name' attribute" in msg for msg in log_cm.output) + ) + + def test_lhs_parallel_solver_factory_failure_uses_default_worker_solver(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + class _NamedSolver: + name = "definitely_missing_solver" + options = {} + + doe.solver = _NamedSolver() + + def _fake_fim(self_obj, experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 0.5, 0.0], [0.0, x + 1.0]]) + + with patch("pyomo.contrib.doe.doe.pyo.SolverFactory", return_value=None): + with patch.object( + DesignOfExperiments, + "_compute_fim_at_point_no_prior", + autospec=True, + side_effect=_fake_fim, + ): + with self.assertLogs("pyomo.contrib.doe.doe", level="DEBUG") as log_cm: + points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=18, + n_exp=2, + lhs_parallel=True, + lhs_n_workers=2, + lhs_combo_parallel=False, + ) + + self.assertEqual(len(points), 2) + self.assertTrue( + any( + "could not construct solver 'definitely_missing_solver'" in msg + for msg in log_cm.output + ) + ) + + def test_lhs_parallel_worker_exception_uses_zero_fim_fallback(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + with patch.object( + DesignOfExperiments, + "_compute_fim_at_point_no_prior", + side_effect=RuntimeError("worker boom"), + ): + with self.assertLogs("pyomo.contrib.doe.doe", level="ERROR") as log_cm: + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=22, + n_exp=2, + lhs_parallel=True, + lhs_n_workers=2, + lhs_combo_parallel=False, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(diag["timed_out"] is False) + self.assertTrue( + any( + "Using zero FIM for this candidate and continuing." in msg + for msg in log_cm.output + ) + ) + + def test_lhs_parallel_candidate_timeout_cancels_pending(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _slow_fim(self_obj, experiment_index, input_values): + time.sleep(0.05) + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 2.0]]) + + with patch.object( + DesignOfExperiments, + "_compute_fim_at_point_no_prior", + autospec=True, + side_effect=_slow_fim, + ): + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=8, + lhs_seed=5, + n_exp=2, + lhs_parallel=True, + lhs_n_workers=2, + lhs_combo_parallel=False, + lhs_max_wall_clock_time=0.001, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(diag["timed_out"]) + + 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) + + with patch.object(doe, "_compute_fim_at_point_no_prior", return_value=None): + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=3, + lhs_seed=8, + n_exp=1, + lhs_parallel=False, + lhs_combo_parallel=False, + ) + + self.assertEqual(len(points), 1) + self.assertTrue(diag["timed_out"]) + + def test_lhs_candidate_subset_padding_to_n_exp(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=3) + call_state = {"n": 0} + + def _fim_with_one_slow_call(experiment_index, input_values): + call_state["n"] += 1 + x = float(input_values[0]) + if call_state["n"] == 2: + time.sleep(0.1) + return np.array([[x + 1.0, 0.0], [0.0, x + 2.0]]) + + with patch.object( + doe, "_compute_fim_at_point_no_prior", side_effect=_fim_with_one_slow_call + ): + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=5, + lhs_seed=16, + n_exp=3, + lhs_parallel=False, + lhs_combo_parallel=False, + lhs_max_wall_clock_time=0.05, + ) + + self.assertEqual(len(points), 3) + self.assertTrue(diag["timed_out"]) + self.assertLess(diag["n_candidates"], 5) + + def test_lhs_combo_parallel_chunk_boundary_matches_serial(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 2.0, 0.0], [0.0, x + 1.0]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_serial, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, lhs_seed=99, n_exp=2, lhs_combo_parallel=False + ) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_parallel, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, + lhs_seed=99, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=3, # C(5,2)=10, non-divisible chunking + lhs_combo_parallel_threshold=1, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_norm) + + def test_lhs_combo_parallel_warning_reports_single_worker_reason(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as log_cm: + points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=27, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=1, + lhs_combo_parallel_threshold=1, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(any("resolved_workers=1 <= 1" in msg for msg in log_cm.output)) + + def test_lhs_combo_parallel_skips_empty_local_combo(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with patch.object( + DesignOfExperiments, + "_evaluate_objective_for_option", + return_value=float("nan"), + ): + points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=4, + lhs_seed=33, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=2, + lhs_combo_parallel_threshold=1, + ) + + self.assertEqual(len(points), 2) + + def test_lhs_combo_scoring_timeout_returns_best_so_far(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + def _slow_obj(fim, objective_option): + time.sleep(0.003) + return float(np.trace(fim)) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with patch.object( + DesignOfExperiments, + "_evaluate_objective_for_option", + side_effect=_slow_obj, + ): + with self.assertLogs( + "pyomo.contrib.doe.doe", level="WARNING" + ) as log_cm: + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=6, + lhs_seed=7, + n_exp=2, + lhs_combo_parallel=False, + lhs_max_wall_clock_time=0.001, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(any("time budget" in msg for msg in log_cm.output)) + self.assertTrue(diag["timed_out"]) + + def test_lhs_combo_scoring_parallel_timeout_returns_best_so_far(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + def _slow_obj(fim, objective_option): + time.sleep(0.003) + return float(np.trace(fim)) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with patch.object( + DesignOfExperiments, + "_evaluate_objective_for_option", + side_effect=_slow_obj, + ): + with self.assertLogs( + "pyomo.contrib.doe.doe", level="WARNING" + ) as log_cm: + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=6, + lhs_seed=12, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=5, + lhs_combo_parallel_threshold=1, + lhs_max_wall_clock_time=0.001, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(any("time budget" in msg for msg in log_cm.output)) + self.assertTrue(diag["timed_out"]) + + def test_lhs_combo_parallel_timeout_cancels_pending(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + def _slow_obj(fim, objective_option): + time.sleep(0.01) + return float(np.trace(fim)) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with patch.object( + DesignOfExperiments, + "_evaluate_objective_for_option", + side_effect=_slow_obj, + ): + with self.assertLogs( + "pyomo.contrib.doe.doe", level="WARNING" + ) as log_cm: + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=9, + lhs_seed=4, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=1, + lhs_combo_parallel_threshold=1, + lhs_max_wall_clock_time=0.02, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(diag["timed_out"]) + self.assertTrue(any("time budget" in msg for msg in log_cm.output)) + + def test_lhs_combo_parallel_submit_loop_timeout_sets_timed_out(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=1) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + times = iter([0.0, 0.0, 0.0, 0.0, 0.01, 0.01, 0.01]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with patch( + "pyomo.contrib.doe.doe.time.perf_counter", + side_effect=lambda: next(times), + ): + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=1, + lhs_seed=19, + n_exp=1, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=1, + lhs_combo_parallel_threshold=1, + lhs_max_wall_clock_time=0.001, + ) + + self.assertEqual(len(points), 1) + self.assertTrue(diag["timed_out"]) + + def test_lhs_combo_parallel_deadline_cancels_pending_futures(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + stage = {"timeout": False} + + class _FakeFuture: + def __init__(self, payload): + self._payload = payload + self.cancelled = False + + def result(self): + return self._payload + + def cancel(self): + self.cancelled = True + return True + + class _FakeExecutor: + def __init__(self, max_workers=None): + self.created = [] + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc, tb): + return False + + def submit(self, fn, *args, **kwargs): + # Return deterministic chunk scores; first future will be "done". + idx = len(self.created) + payload = ( + (10.0 - idx, (0, 1), False) + if idx == 0 + else (9.0 - idx, (0, 2), False) + ) + fut = _FakeFuture(payload) + self.created.append(fut) + return fut + + def _fake_wait(pending, return_when=None): + pending_list = list(pending) + done = {pending_list[0]} + still_pending = set(pending_list[1:]) + stage["timeout"] = True + return done, still_pending + + def _fake_perf_counter(): + return 1.0 if stage["timeout"] else 0.0 + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + with patch("pyomo.contrib.doe.doe._cf.ThreadPoolExecutor", _FakeExecutor): + with patch("pyomo.contrib.doe.doe._cf.wait", side_effect=_fake_wait): + with patch( + "pyomo.contrib.doe.doe.time.perf_counter", + side_effect=_fake_perf_counter, + ): + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=3, + lhs_seed=23, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=1, + lhs_combo_parallel_threshold=1, + lhs_max_wall_clock_time=0.5, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(diag["timed_out"]) + + def test_lhs_combo_parallel_minimize_objective_update(self): + doe = self._make_template_doe("trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, 3.0 * x + 1.0]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_serial, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, lhs_seed=52, n_exp=2, lhs_combo_parallel=False + ) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_parallel, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, + lhs_seed=52, + n_exp=2, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=2, + lhs_combo_parallel_threshold=1, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_norm) + + def test_lhs_combo_serial_minimize_objective_update(self): + doe = self._make_template_doe("trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + lhs_n_samples = 4 + lhs_seed = 61 + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 0.25, 0.0], [0.0, 2.5 * x + 0.25]]) + + 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)) + candidate_fims = [_fake_fim(0, pt) for pt in candidate_points] + + best_combo = None + best_obj = np.inf + for combo in combinations(range(len(candidate_points)), 2): + fim_total = sum((candidate_fims[idx] for idx in combo), np.zeros((2, 2))) + obj_val = doe._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] + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + got_points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=lhs_n_samples, + lhs_seed=lhs_seed, + n_exp=2, + lhs_combo_parallel=False, + ) + + 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_score_chunk_minimize_branch(self): + doe = self._make_template_doe("trace") + self._build_template_model_for_multi_experiment(doe, n_exp=3) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 0.75, 0.0], [0.0, 2.0 * x + 0.75]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_serial, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, lhs_seed=77, n_exp=3, lhs_combo_parallel=False + ) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_parallel, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, + lhs_seed=77, + n_exp=3, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=2, + lhs_combo_parallel_threshold=1, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_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 + unit_samples = np.array([0.2, 0.5, 0.8]) + + class _FakeLHS: + def __init__(self, d, seed=None): + self.d = d + + def random(self, n): + assert self.d == 1 + assert n == lhs_n_samples + return unit_samples.reshape((n, 1)) + + first_exp_block = doe.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe._get_experiment_input_vars(first_exp_block) + lb = float(exp_input_vars[0].lb) + ub = float(exp_input_vars[0].ub) + expected_points = [[float(lb + s * (ub - lb))] for s in unit_samples] + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 1.0]]) + + with patch("pyomo.contrib.doe.doe.LatinHypercube", _FakeLHS): + with patch("pyomo.contrib.doe.doe._combinations", return_value=iter(())): + with patch.object( + doe, "_compute_fim_at_point_no_prior", side_effect=_fake_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=13, + n_exp=3, + lhs_combo_parallel=False, + ) + + 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) + self.assertTrue( + any( + "Falling back to the first n_exp candidate points." in msg + for msg in log_cm.output + ) + ) + + def test_lhs_fim_evaluation_timeout_stops_early(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=2) + + def _slow_fim(experiment_index, input_values): + time.sleep(0.01) + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, x + 2.0]]) + + with patch.object( + doe, "_compute_fim_at_point_no_prior", side_effect=_slow_fim + ) as p: + points, diag = doe._lhs_initialize_experiments( + lhs_n_samples=10, + lhs_seed=3, + n_exp=2, + lhs_parallel=False, + lhs_combo_parallel=False, + lhs_max_wall_clock_time=0.03, + ) + + self.assertEqual(len(points), 2) + self.assertTrue(diag["timed_out"]) + self.assertLess(p.call_count, 10) + self.assertLess(diag["n_candidates"], 10) + + def test_lhs_combo_scoring_n_exp_3_parallel_matches_serial(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=3) + + def _fake_fim(experiment_index, input_values): + x = float(input_values[0]) + return np.array([[x + 1.0, 0.0], [0.0, 0.5 * x + 2.0]]) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_serial, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, lhs_seed=1234, n_exp=3, lhs_combo_parallel=False + ) + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + points_parallel, _ = doe._lhs_initialize_experiments( + lhs_n_samples=5, + lhs_seed=1234, + n_exp=3, + lhs_combo_parallel=True, + lhs_n_workers=2, + lhs_combo_chunk_size=3, + lhs_combo_parallel_threshold=1, + ) + + serial_norm = sorted(tuple(np.round(p, 8)) for p in points_serial) + parallel_norm = sorted(tuple(np.round(p, 8)) for p in points_parallel) + self.assertEqual(serial_norm, parallel_norm) + self.assertEqual(len(points_serial), 3) + + def test_lhs_combo_scoring_n_exp_3_matches_oracle(self): + doe = self._make_template_doe("pseudo_trace") + self._build_template_model_for_multi_experiment(doe, n_exp=3) + lhs_n_samples = 5 + lhs_seed = 2 + + 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 + 0.5]]) + + # Recreate the exact candidate points from LHS generation (independent + # oracle for combination scoring logic). + 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)) + + # Oracle over all combinations of size 3. + fims = [_fake_fim(0, pt) for pt in candidate_points] + best_obj = -np.inf + best_combo = None + for combo in combinations(range(len(candidate_points)), 3): + fim_total = sum((fims[i] for i in combo), np.zeros((2, 2))) + obj = float(np.trace(fim_total)) + if obj > best_obj: + best_obj = obj + best_combo = combo + expected_points = [list(candidate_points[i]) for i in best_combo] + + with patch.object(doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim): + got_points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=lhs_n_samples, + lhs_seed=lhs_seed, + n_exp=3, + lhs_combo_parallel=False, + ) + + 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_independent_oracle_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 + + first_exp_block = doe.model.param_scenario_blocks[0].exp_blocks[0] + exp_input_vars = doe._get_experiment_input_vars(first_exp_block) + self.assertEqual(len(exp_input_vars), 1) + lb = float(exp_input_vars[0].lb) + ub = float(exp_input_vars[0].ub) + + unit_samples = np.array([0.1, 0.4, 0.9]) + scaled_points = [float(lb + s * (ub - lb)) for s in unit_samples] + + class _FakeLHS: + def __init__(self, d, seed=None): + self.d = d + + def random(self, n): + assert self.d == 1 + assert n == lhs_n_samples + return unit_samples.reshape((n, 1)) + + 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]]) + + with patch("pyomo.contrib.doe.doe.LatinHypercube", _FakeLHS): + with patch.object( + doe, "_compute_fim_at_point_no_prior", side_effect=_fake_fim + ): + got_points, _ = doe._lhs_initialize_experiments( + lhs_n_samples=lhs_n_samples, + lhs_seed=123, + n_exp=2, + lhs_combo_parallel=False, + ) + + # Independent brute-force oracle over explicit candidate points + best_obj = -np.inf + best_combo = None + for combo in combinations(range(len(scaled_points)), 2): + f1 = _fake_fim(0, [scaled_points[combo[0]]]) + f2 = _fake_fim(0, [scaled_points[combo[1]]]) + obj_val = float(np.trace(f1 + f2)) + if obj_val > best_obj: + best_obj = obj_val + best_combo = combo + expected_points = [[scaled_points[i]] for i in best_combo] + + 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 = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + + doe = DesignOfExperiments( + experiment=exp_list, + objective_option="determinant", + step=1e-2, + solver=solver, + ) + doe.optimize_experiments() + + scenario = doe.results["Scenarios"][0] + got_hours = sorted( + exp["Experiment 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["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 = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + # 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 = doe.results["Scenarios"][0] + got_hours = sorted( + exp["Experiment 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["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 = doe.results["Scenarios"][0] + 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["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 = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + 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 = doe.results["Scenarios"][0] + 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["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 with a warning. + doe = self._make_template_doe("pseudo_trace") + with patch( + "pyomo.contrib.doe.doe.np.linalg.inv", side_effect=RuntimeError("boom") + ): + with self.assertLogs("pyomo.contrib.doe.doe", level="WARNING") as log_cm: + doe.optimize_experiments(n_exp=1) + + scenario = doe.results["Scenarios"][0] + self.assertTrue(np.isnan(scenario["log10 A-opt"])) + self.assertTrue( + any("failed to compute log10 A-opt" in msg for msg in log_cm.output) + ) + + 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 = SolverFactory("ipopt") + solver.options["linear_solver"] = "ma57" + solver.options["halt_on_ampl_error"] = "yes" + solver.options["max_iter"] = 3000 + doe = DesignOfExperiments( + experiment=[exp], + objective_option="determinant", + step=1e-2, + solver=solver, + _Cholesky_option=False, + _only_compute_fim_lower=False, + ) + original_solve = doe.solver.solve + + class _MockSolverInfo: + status = "ok" + termination_condition = "optimal" + message = "mock-solve" + + class _MockResults: + solver = _MockSolverInfo() + + def _solve_real_for_square_then_mock(*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 original_solve(*args, **kwargs) + return _MockResults() + + with patch.object( + doe.solver, "solve", side_effect=_solve_real_for_square_then_mock + ): + 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(doe.results["Scenarios"][0]["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()