diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py index 72eeb606fc2..301c2bebb30 100644 --- a/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py +++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py @@ -52,23 +52,6 @@ def simple_reaction_model(data): # fix all of the regressed parameters model.k.fix() - # =================================================================== - # Stage-specific cost computations - def ComputeFirstStageCost_rule(model): - return 0 - - model.FirstStageCost = Expression(rule=ComputeFirstStageCost_rule) - - def AllMeasurements(m): - return (float(data['y']) - m.y) ** 2 - - model.SecondStageCost = Expression(rule=AllMeasurements) - - def total_cost_rule(m): - return m.FirstStageCost + m.SecondStageCost - - model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize) - return model @@ -90,6 +73,8 @@ def label_model(self): m.experiment_outputs.update( [(m.x1, self.data['x1']), (m.x2, self.data['x2']), (m.y, self.data['y'])] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None), (m.x1, None), (m.x2, None)]) return m @@ -161,7 +146,7 @@ def main(): # Only estimate the parameter k[1]. The parameter k[2] will remain fixed # at its initial value - pest = parmest.Estimator(exp_list) + pest = parmest.Estimator(exp_list, obj_function="SSE") obj, theta = pest.theta_est() print(obj) print(theta) @@ -174,9 +159,11 @@ def main(): # ======================================================================= # Estimate both k1 and k2 and compute the covariance matrix - pest = parmest.Estimator(exp_list) - n = 15 # total number of data points used in the objective (y in 15 scenarios) - obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=n) + pest = parmest.Estimator(exp_list, obj_function="SSE") + # Calculate the objective value and estimated parameters + obj, theta = pest.theta_est() + # Compute the covariance matrix using the reduced Hessian method + cov = pest.cov_est(method="reduced_hessian") print(obj) print(theta) print(cov) diff --git a/pyomo/contrib/parmest/examples/reactor_design/multistart_example.py b/pyomo/contrib/parmest/examples/reactor_design/multistart_example.py new file mode 100644 index 00000000000..cfb9eea9798 --- /dev/null +++ b/pyomo/contrib/parmest/examples/reactor_design/multistart_example.py @@ -0,0 +1,56 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from pyomo.common.dependencies import numpy as np, pandas as pd +from itertools import product +from os.path import join, abspath, dirname +import pyomo.contrib.parmest.parmest as parmest +from pyomo.contrib.parmest.examples.reactor_design.reactor_design import ( + ReactorDesignExperiment, +) + + +def main(): + + # Read in data + file_dirname = dirname(abspath(str(__file__))) + file_name = abspath(join(file_dirname, "reactor_data.csv")) + data = pd.read_csv(file_name) + + # Create an experiment list + exp_list = [ReactorDesignExperiment(data, i) for i in range(data.shape[0])] + + # Solver options belong here (Ipopt options shown as example) + solver_options = {"max_iter": 1000, "tol": 1e-6} + + pest = parmest.Estimator( + exp_list, obj_function="SSE", solver_options=solver_options + ) + + # Single-start estimation + obj, theta = pest.theta_est() + print("Single-start objective:", obj) + print("Single-start theta:\n", theta) + + # Multistart estimation + results_df, best_theta, best_obj = pest.theta_est_multistart( + n_restarts=10, + multistart_sampling_method="uniform_random", + seed=42, + save_results=False, # True if you want CSV via file_name= + ) + + print("\nMultistart best objective:", best_obj) + print("Multistart best theta:", best_theta) + print("\nAll multistart results:") + print(results_df) + + +if __name__ == "__main__": + main() diff --git a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py index 630802ddba9..a5a644a4c7b 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py @@ -33,11 +33,17 @@ def main(): pest = parmest.Estimator(exp_list, obj_function='SSE') - # Parameter estimation with covariance - obj, theta, cov = pest.theta_est(calc_cov=True, cov_n=17) - print(obj) + # Parameter estimation + obj, theta = pest.theta_est() + print("Least squares objective value:", obj) + print("Estimated parameters (theta):\n") print(theta) + # Compute the covariance matrix at the estimated parameter + cov = pest.cov_est() + print("Covariance matrix:\n") + print(cov) + if __name__ == "__main__": main() diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py index 1b3360a2d93..6a065d20a63 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py @@ -22,15 +22,15 @@ def reactor_design_model(): # Create the concrete model model = pyo.ConcreteModel() - # Rate constants - model.k1 = pyo.Param( - initialize=5.0 / 6.0, within=pyo.PositiveReals, mutable=True + # Rate constants, make unknown parameters variables + model.k1 = pyo.Var( + initialize=5.0 / 6.0, within=pyo.PositiveReals, bounds=(0.1, 10.0) ) # min^-1 - model.k2 = pyo.Param( - initialize=5.0 / 3.0, within=pyo.PositiveReals, mutable=True + model.k2 = pyo.Var( + initialize=5.0 / 3.0, within=pyo.PositiveReals, bounds=(0.1, 10.0) ) # min^-1 - model.k3 = pyo.Param( - initialize=1.0 / 6000.0, within=pyo.PositiveReals, mutable=True + model.k3 = pyo.Var( + initialize=1.0 / 6000.0, within=pyo.PositiveReals, bounds=(1e-5, 1e-3) ) # m^3/(gmol min) # Inlet concentration of A, gmol/m^3 @@ -119,6 +119,11 @@ def label_model(self): (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2, m.k3] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update( + [(m.ca, None), (m.cb, None), (m.cc, None), (m.cd, None)] + ) + return m def get_labeled_model(self): diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/multistart_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/multistart_example.py new file mode 100644 index 00000000000..7ff8f55ad4b --- /dev/null +++ b/pyomo/contrib/parmest/examples/rooney_biegler/multistart_example.py @@ -0,0 +1,62 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from pyomo.common.dependencies import numpy as np, pandas as pd +from itertools import product +import pyomo.contrib.parmest.parmest as parmest +from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( + RooneyBieglerExperiment, +) + + +def main(): + + # Data + data = pd.DataFrame( + data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], + columns=['hour', 'y'], + ) + + # Create an experiment list + exp_list = [] + for i in range(data.shape[0]): + exp_list.append(RooneyBieglerExperiment(data.loc[i, :])) + + # View one model + # exp0_model = exp_list[0].get_labeled_model() + # exp0_model.pprint() + + # Solver options belong here (Ipopt options shown as example) + solver_options = {"max_iter": 1000, "tol": 1e-6} + + pest = parmest.Estimator( + exp_list, obj_function="SSE", solver_options=solver_options + ) + + # Single-start estimation + obj, theta = pest.theta_est() + print("Single-start objective:", obj) + print("Single-start theta:\n", theta) + + # Multistart estimation + results_df, best_theta, best_obj = pest.theta_est_multistart( + n_restarts=10, + multistart_sampling_method="uniform_random", + seed=42, + save_results=False, # True if you want CSV via file_name= + ) + + print("\nMultistart best objective:", best_obj) + print("Multistart best theta:", best_theta) + print("\nAll multistart results:") + print(results_df) + + +if __name__ == "__main__": + main() diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/profile_likelihood_example.py b/pyomo/contrib/parmest/examples/rooney_biegler/profile_likelihood_example.py new file mode 100644 index 00000000000..8d3bb2b22ef --- /dev/null +++ b/pyomo/contrib/parmest/examples/rooney_biegler/profile_likelihood_example.py @@ -0,0 +1,74 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from pyomo.common.dependencies import pandas as pd +import pyomo.contrib.parmest.parmest as parmest +from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import ( + RooneyBieglerExperiment, +) + + +def main(): + # Data + data = pd.DataFrame( + data=[[1, 8.3], [2, 10.3], [3, 19.0], [4, 16.0], [5, 15.6], [7, 19.8]], + columns=["hour", "y"], + ) + + # Build experiment list + exp_list = [RooneyBieglerExperiment(data.loc[i, :]) for i in range(data.shape[0])] + + # Create estimator + pest = parmest.Estimator(exp_list, obj_function="SSE") + + # Compute profile likelihood for both unknown parameters. + # Use a small grid for quick terminal runs. + profile_results = pest.profile_likelihood( + profiled_theta=["asymptote", "rate_constant"], + n_grid=9, + solver="ef_ipopt", + warmstart="neighbor", + # Demonstrate baseline from multistart integration: + use_multistart_for_baseline=True, + baseline_multistart_kwargs={ + "n_restarts": 5, + "multistart_sampling_method": "uniform_random", + "seed": 7, + }, + ) + + # Display a compact summary table + profiles = profile_results["profiles"] + print("\nBaseline:") + print(profile_results["baseline"]) + print("\nProfile results (first 12 rows):") + print( + profiles[ + [ + "profiled_theta", + "theta_value", + "obj", + "delta_obj", + "lr_stat", + "status", + "success", + ] + ].head(12) + ) + + # Plot profile curves to file for terminal/non-GUI usage + out_file = "rooney_biegler_profile_likelihood.png" + parmest.graphics.profile_likelihood_plot( + profile_results, alpha=0.95, filename=out_file + ) + print(f"\nSaved plot to: {out_file}") + + +if __name__ == "__main__": + main() diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py index bdb494bca03..4ed4e8fc947 100644 --- a/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py +++ b/pyomo/contrib/parmest/examples/rooney_biegler/rooney_biegler.py @@ -43,8 +43,8 @@ def rooney_biegler_model(data, theta=None): if theta is None: theta = {'asymptote': 15, 'rate_constant': 0.5} - model.asymptote = pyo.Var(initialize=theta['asymptote']) - model.rate_constant = pyo.Var(initialize=theta['rate_constant']) + model.asymptote = pyo.Var(initialize=theta['asymptote'], bounds=(0.1, 100)) + model.rate_constant = pyo.Var(initialize=theta['rate_constant'], bounds=(0, 10)) # Fix the unknown parameters model.asymptote.fix() diff --git a/pyomo/contrib/parmest/graphics.py b/pyomo/contrib/parmest/graphics.py index bed9f9eb824..9361c737b61 100644 --- a/pyomo/contrib/parmest/graphics.py +++ b/pyomo/contrib/parmest/graphics.py @@ -488,6 +488,99 @@ def pairwise_plot( plt.close() +def profile_likelihood_plot( + profile_results, alpha=None, filename=None, by="profiled_theta", y="lr_stat" +): + """ + Plot profile likelihood curves from Estimator.profile_likelihood results. + + Parameters + ---------- + profile_results: dict or DataFrame + Result dict returned by ``Estimator.profile_likelihood()`` or the + ``profiles`` DataFrame directly. + alpha: float, optional + If provided, draw the chi-square threshold line using dof=1. + filename: string, optional + Filename used to save the figure. + by: string, optional + Grouping column for separate panels. Default is ``profiled_theta``. + y: string, optional + Y-axis column to plot. Default is ``lr_stat``. + """ + assert isinstance(profile_results, (dict, pd.DataFrame)) + assert isinstance(alpha, (type(None), int, float)) + assert isinstance(filename, (type(None), str)) + assert isinstance(by, str) + assert isinstance(y, str) + + if isinstance(profile_results, dict): + if "profiles" not in profile_results: + raise KeyError("profile_results dict must contain 'profiles'.") + profiles = profile_results["profiles"] + else: + profiles = profile_results + + if not isinstance(profiles, pd.DataFrame): + raise TypeError("profiles must be a pandas DataFrame.") + if profiles.empty: + raise ValueError("profiles DataFrame is empty.") + required_cols = {"theta_value", y, by} + missing = required_cols.difference(profiles.columns) + if missing: + raise KeyError( + f"profiles DataFrame is missing required columns: {sorted(missing)}" + ) + + groups = list(profiles[by].dropna().unique()) + if len(groups) == 0: + raise ValueError(f"No non-null group values found in '{by}'.") + + fig, axes = plt.subplots( + nrows=len(groups), + ncols=1, + figsize=(6.5, max(3.0, 2.8 * len(groups))), + squeeze=False, + ) + + threshold = None + if alpha is not None: + threshold = stats.chi2.ppf(float(alpha), df=1) + + for i, grp in enumerate(groups): + ax = axes[i, 0] + subset = profiles[profiles[by] == grp].sort_values("theta_value") + ax.plot(subset["theta_value"], subset[y], marker="o", linewidth=1.5) + if "success" in subset.columns: + failed = subset[~subset["success"].astype(bool)] + if len(failed) > 0: + ax.scatter( + failed["theta_value"], + failed[y], + marker="x", + color="tab:red", + s=24, + label="failed", + ) + if threshold is not None: + ax.axhline(threshold, color="tab:orange", linestyle="--", linewidth=1.0) + ax.set_xlabel(str(grp)) + ax.set_ylabel(y) + ax.set_title(f"{by}: {grp}") + if ( + "success" in subset.columns + and len(subset[~subset["success"].astype(bool)]) > 0 + ): + ax.legend(loc="best", prop={"size": 8}) + + fig.tight_layout() + if filename is None: + plt.show() + else: + plt.savefig(filename) + plt.close() + + def fit_rect_dist(theta_values, alpha): """ Fit an alpha-level rectangular distribution to theta values diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 21138d472da..3885985f00b 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -11,6 +11,7 @@ #### Wrapping mpi-sppy functionality and local option Jan 2021, Feb 2021 #### Redesign with Experiment class Dec 2023 +# Options for using mpi-sppy or local EF only used in the deprecatedEstimator class # TODO: move use_mpisppy to a Pyomo configuration option # False implies always use the EF that is local to parmest use_mpisppy = True # Use it if we can but use local if not. @@ -41,6 +42,7 @@ import logging import types import json +import time from collections.abc import Callable from itertools import combinations from functools import singledispatchmethod @@ -80,6 +82,7 @@ logger = logging.getLogger(__name__) +# Only used in the deprecatedEstimator class def ef_nonants(ef): # Wrapper to call someone's ef_nonants # (the function being called is very short, but it might be changed) @@ -89,6 +92,7 @@ def ef_nonants(ef): return local_ef.ef_nonants(ef) +# Only used in the deprecatedEstimator class def _experiment_instance_creation_callback( scenario_name, node_names=None, cb_data=None ): @@ -354,14 +358,24 @@ def _count_total_experiments(experiment_list): Returns ------- - total_number_data : int + total_data_points : int The total number of data points in the list of experiments """ - total_number_data = 0 + total_data_points = 0 for experiment in experiment_list: - total_number_data += len(experiment.get_labeled_model().experiment_outputs) + # 1. Identify the first parent component of the experiment outputs + output_vars = experiment.get_labeled_model().experiment_outputs + + # 1. Identify the first parent component + first_var_key = list(output_vars.keys())[0] + first_parent = first_var_key.parent_component() + # 2. Count only the keys that belong to this specific parent + first_parent_indices = [ + v for v in output_vars.keys() if v.parent_component() is first_parent + ] + total_data_points += len(first_parent_indices) - return total_number_data + return total_data_points class CovarianceMethod(Enum): @@ -911,7 +925,32 @@ def _expand_indexed_unknowns(self, model_temp): def _create_parmest_model(self, experiment_number): """ - Modify the Pyomo model for parameter estimation + Build a parmest-ready model for a single experiment. + + This helper retrieves the labeled experiment model, prepares objective + components needed by parmest, and converts unknown parameters to + decision variables. The returned model is the one used to populate EF + scenario blocks. + + Parameters + ---------- + experiment_number : int + Index into ``self.exp_list`` selecting which experiment model to + load. + + Returns + ------- + ConcreteModel + A model configured for parmest optimization, including: + 1. a ``Total_Cost_Objective`` (if ``self.obj_function`` is set) + 2. converted unknown-parameter variables (unfixed) + + Notes + ----- + - Existing user objectives are deactivated before parmest objective + components are attached. + - Reserved component names are checked to avoid overriding user model + components. """ model = _get_labeled_model(self.exp_list[experiment_number]) @@ -969,197 +1008,500 @@ def _instance_creation_callback(self, experiment_number=None, cb_data=None): model = self._create_parmest_model(experiment_number) return model - def _Q_opt( + def _create_scenario_blocks( self, - ThetaVals=None, - solver="ef_ipopt", - return_values=[], bootlist=None, - calc_cov=NOTSET, - cov_n=NOTSET, + theta_vals=None, + fix_theta=False, + multistart=False, + fixed_theta_values=None, ): """ - Set up all thetas as first stage Vars, return resulting theta - values as well as the objective function value. + Build the block-based extensive form (EF) model for estimation. - """ - if solver == "k_aug": - raise RuntimeError("k_aug no longer supported.") + The EF includes: + 1. a master theta variable container (``model.parmest_theta``), + 2. one child block per selected experiment (``model.exp_scenarios``), + 3. optional theta-linking constraints between master and child blocks, + 4. a single aggregate objective over all child blocks. - # (Bootstrap scenarios will use indirection through the bootlist) - if bootlist is None: - scenario_numbers = list(range(len(self.exp_list))) - scen_names = ["Scenario{}".format(i) for i in scenario_numbers] - else: - scen_names = ["Scenario{}".format(i) for i in range(len(bootlist))] + In multistart mode, this method also refreshes experiment-level cached + model state before rebuilding each scenario so per-start initializations + are applied to the model that is actually solved. - # get the probability constant that is applied to the objective function - # parmest solves the estimation problem by applying equal probabilities to - # the objective function of all the scenarios from the experiment list - self.obj_probability_constant = len(scen_names) + Parameters + ---------- + bootlist : list, optional + Experiment indices to include. If ``None``, all experiments in + ``self.exp_list`` are used. + theta_vals : dict, optional + Theta values to apply as initial values to parent and child theta + variables. When ``multistart=True``, these values are also pushed to + experiment ``theta_initial`` (if present) before rebuilding. + fix_theta : bool, optional + If ``True``, theta variables are fixed in each scenario and no + linking constraints are created. + multistart : bool, optional + If ``True``, force experiment model refresh between starts to avoid + stale cached model reuse. - # tree_model.CallbackModule = None - outer_cb_data = dict() - outer_cb_data["callback"] = self._instance_creation_callback - if ThetaVals is not None: - outer_cb_data["ThetaVals"] = ThetaVals - if bootlist is not None: - outer_cb_data["BootList"] = bootlist - outer_cb_data["cb_data"] = None # None is OK - outer_cb_data["theta_names"] = self.estimator_theta_names + Returns + ------- + ConcreteModel + EF model used by parmest solve routines. - options = {"solver": "ipopt"} - scenario_creator_options = {"cb_data": outer_cb_data} - if use_mpisppy: - ef = sputils.create_EF( - scen_names, - _experiment_instance_creation_callback, - EF_name="_Q_opt", - suppress_warnings=True, - scenario_creator_kwargs=scenario_creator_options, - ) - else: - ef = local_ef.create_EF( - scen_names, - _experiment_instance_creation_callback, - EF_name="_Q_opt", - suppress_warnings=True, - scenario_creator_kwargs=scenario_creator_options, + Raises + ------ + ValueError + If the selected scenario set is empty. + """ + # Build a clean parent EF container and attach one scenario model per block. + model = pyo.ConcreteModel() + if multistart: + template_experiment = self.exp_list[0] + if theta_vals is not None and hasattr(template_experiment, "theta_initial"): + template_experiment.theta_initial = dict(theta_vals) + if hasattr(template_experiment, "model"): + template_experiment.model = None + template_model = self._create_parmest_model(0) + expanded_theta_names = self._expand_indexed_unknowns(template_model) + model._parmest_theta_names = tuple(expanded_theta_names) + model.parmest_theta = pyo.Var(model._parmest_theta_names) + + fixed_theta_values = fixed_theta_values or {} + invalid_fixed_theta = set(fixed_theta_values).difference(expanded_theta_names) + if invalid_fixed_theta: + raise ValueError( + f"Unknown theta name(s) in fixed_theta_values: {sorted(invalid_fixed_theta)}" ) - self.ef_instance = ef + fixed_theta_names = set(fixed_theta_values.keys()) + + for name in expanded_theta_names: + template_theta_var = template_model.find_component(name) + parent_theta_var = model.parmest_theta[name] + parent_theta_var.set_value(pyo.value(template_theta_var)) + if theta_vals is not None and name in theta_vals: + parent_theta_var.set_value(theta_vals[name]) + if name in fixed_theta_values: + parent_theta_var.set_value(fixed_theta_values[name]) + if fix_theta: + parent_theta_var.fix() + elif name in fixed_theta_names: + parent_theta_var.fix() + else: + parent_theta_var.unfix() - # Solve the extensive form with ipopt - if solver == "ef_ipopt": - if calc_cov is NOTSET or not calc_cov: - # Do not calculate the reduced hessian + # Set the number of experiments to use, either from bootlist or all experiments + scenario_numbers = ( + list(bootlist) if bootlist is not None else list(range(len(self.exp_list))) + ) + self.obj_probability_constant = len(scenario_numbers) + if self.obj_probability_constant <= 0: + raise ValueError("At least one scenario is required to build the EF model.") + + # Create indexed block for holding scenario models + model.exp_scenarios = pyo.Block(range(self.obj_probability_constant)) + for i, experiment_number in enumerate(scenario_numbers): + if multistart: + experiment = self.exp_list[experiment_number] + if theta_vals is not None and hasattr(experiment, "theta_initial"): + experiment.theta_initial = dict(theta_vals) + if hasattr(experiment, "model"): + experiment.model = None + parmest_model = self._create_parmest_model(experiment_number) + for name in expanded_theta_names: + child_theta_var = parmest_model.find_component(name) + parent_theta_var = model.parmest_theta[name] + if theta_vals is not None and name in theta_vals: + child_theta_var.set_value(theta_vals[name]) + else: + child_theta_var.set_value(pyo.value(parent_theta_var)) + if name in fixed_theta_values: + child_theta_var.set_value(fixed_theta_values[name]) + if fix_theta: + child_theta_var.fix() + elif name in fixed_theta_names: + child_theta_var.fix() + else: + child_theta_var.unfix() + model.exp_scenarios[i].transfer_attributes_from(parmest_model) + + model.theta_link_constraints = pyo.ConstraintList() + if not fix_theta: + for name in expanded_theta_names: + if name in fixed_theta_names: + continue + parent_theta_var = model.parmest_theta[name] + for i in range(self.obj_probability_constant): + child_theta_var = model.exp_scenarios[i].find_component(name) + model.theta_link_constraints.add( + child_theta_var == parent_theta_var + ) - solver = SolverFactory('ipopt') - if self.solver_options is not None: - for key in self.solver_options: - solver.options[key] = self.solver_options[key] + for block in model.exp_scenarios.values(): + for obj in block.component_objects(pyo.Objective): + obj.deactivate() - solve_result = solver.solve(self.ef_instance, tee=self.tee) - assert_optimal_termination(solve_result) - elif calc_cov is not NOTSET and calc_cov: - # parmest makes the fitted parameters stage 1 variables - ind_vars = [] - for nd_name, Var, sol_val in ef_nonants(ef): - ind_vars.append(Var) - # calculate the reduced hessian - solve_result, inv_red_hes = ( - inverse_reduced_hessian.inv_reduced_hessian_barrier( - self.ef_instance, - independent_variables=ind_vars, - solver_options=self.solver_options, - tee=self.tee, - ) + # Make an objective that sums over all scenario blocks and divides by number of experiments + def total_obj(m): + return ( + sum( + block.Total_Cost_Objective.expr + for block in m.exp_scenarios.values() ) + / self.obj_probability_constant + ) - if self.diagnostic_mode: - print( - ' Solver termination condition = ', - str(solve_result.solver.termination_condition), - ) + model.Obj = pyo.Objective(rule=total_obj, sense=pyo.minimize) + self.ef_instance = model + return model - # assume all first stage are thetas... - theta_vals = {} - for nd_name, Var, sol_val in ef_nonants(ef): - # process the name - # the scenarios are blocks, so strip the scenario name - var_name = Var.name[Var.name.find(".") + 1 :] - theta_vals[var_name] = sol_val + def _generate_initial_theta( + self, + seed=None, + n_restarts=None, + multistart_sampling_method=None, + user_provided_df=None, + experiment_number=0, + ): + """ + Create the canonical multistart initialization/results DataFrame. - obj_val = pyo.value(ef.EF_Obj) - self.obj_value = obj_val - self.estimated_theta = theta_vals + Output schema is: + 1. theta columns (canonical order, quote-normalized names), + 2. ``converged_`` columns, + 3. ``final objective``, ``solver termination``, ``solve_time``. - if calc_cov is not NOTSET and calc_cov: - # Calculate the covariance matrix + Initial theta rows are either sampled from bounds or taken from a + user-provided DataFrame. - if not isinstance(cov_n, int): - raise TypeError( - f"Expected an integer for the 'cov_n' argument. " - f"Got {type(cov_n)}." - ) - num_unknowns = max( - [ - len(experiment.get_labeled_model().unknown_parameters) - for experiment in self.exp_list - ] - ) - assert cov_n > num_unknowns, ( - "The number of datapoints must be greater than the " - "number of parameters to estimate." - ) + Parameters + ---------- + seed : int, optional + Random seed used by stochastic samplers. + n_restarts : int, optional + Number of starts to generate for sampled methods. Ignored when + ``user_provided_df`` is provided. + multistart_sampling_method : str, optional + Sampling method. Supported values: + ``uniform_random``, ``latin_hypercube``, ``sobol_sampling``. + user_provided_df : DataFrame, optional + Explicit initialization table. Must contain exactly the theta + columns (order may vary). Values must be finite and within bounds. + experiment_number : int, optional + Experiment index used to discover canonical theta names and bounds. - # Number of data points considered - n = cov_n + Returns + ------- + DataFrame + Canonical initialization/results table ready for multistart solve + bookkeeping. + + Raises + ------ + ValueError + For missing/invalid bounds, invalid sampling method, malformed + user-provided starts, non-finite values, or out-of-bound starts. + TypeError + For invalid input types (for example, non-DataFrame + ``user_provided_df`` or non-integer ``n_restarts`` when required). + RuntimeError + If expected theta components cannot be located on the model. + """ + parmest_model = self._create_parmest_model(experiment_number) - # Extract number of fitted parameters - l = len(theta_vals) + raw_theta_names = self._expand_indexed_unknowns(parmest_model) + theta_names = [n.replace("'", "") for n in raw_theta_names] + if len(theta_names) != len(set(theta_names)): + raise ValueError(f"Duplicate theta names are not allowed: {theta_names}") - # Assumption: Objective value is sum of squared errors - sse = obj_val + theta_vars = [parmest_model.find_component(name) for name in raw_theta_names] + if any(v is None for v in theta_vars): + raise RuntimeError( + "Failed to locate one or more theta components on model." + ) - '''Calculate covariance assuming experimental observation errors - are independent and follow a Gaussian distribution - with constant variance. + lower_bound = np.array([v.lb for v in theta_vars], dtype=float) + upper_bound = np.array([v.ub for v in theta_vars], dtype=float) - The formula used in parmest was verified against equations - (7-5-15) and (7-5-16) in "Nonlinear Parameter Estimation", - Y. Bard, 1974. + if np.any(np.isnan(lower_bound)) or np.any(np.isnan(upper_bound)): + raise ValueError( + "The lower and upper bounds for the theta values must be defined." + ) + if np.any(lower_bound > upper_bound): + raise ValueError( + "Each lower bound must be less than or equal to the corresponding upper bound." + ) - This formula is also applicable if the objective is scaled by a - constant; the constant cancels out. - (was scaled by 1/n because it computes an expected value.) - ''' - cov = 2 * sse / (n - l) * inv_red_hes - cov = pd.DataFrame( - cov, index=theta_vals.keys(), columns=theta_vals.keys() + if user_provided_df is not None: + if not isinstance(user_provided_df, pd.DataFrame): + raise TypeError("user_provided_df must be a pandas DataFrame.") + if user_provided_df.shape[1] != len(theta_names): + raise ValueError( + "user_provided_df must have exactly one column per theta name." + ) + clean_cols = [str(c).replace("'", "") for c in user_provided_df.columns] + if len(clean_cols) != len(set(clean_cols)): + raise ValueError("Duplicate theta columns are not allowed.") + if set(clean_cols) != set(theta_names): + raise ValueError( + f"Provided columns {clean_cols} do not match expected theta names {theta_names}." + ) + df_multistart = user_provided_df.copy() + df_multistart.columns = clean_cols + df_multistart = df_multistart.reindex(columns=theta_names) + if df_multistart.shape[0] == 0: + raise ValueError("user_provided_df must contain at least one row.") + if n_restarts is not None and n_restarts != df_multistart.shape[0]: + raise ValueError( + "n_restarts must match the number of rows in user_provided_df." ) + theta_vals_multistart = df_multistart.to_numpy(dtype=float) + n_restarts = df_multistart.shape[0] + elif multistart_sampling_method == "uniform_random": + if not isinstance(n_restarts, int): + raise TypeError("n_restarts must be an integer.") + if n_restarts <= 0: + raise ValueError("n_restarts must be greater than zero.") + # Use a local RNG to avoid mutating global random state. + rng = np.random.default_rng(seed) + theta_vals_multistart = rng.uniform( + low=lower_bound, high=upper_bound, size=(n_restarts, len(theta_names)) + ) - theta_vals = pd.Series(theta_vals) + elif multistart_sampling_method == "latin_hypercube": + if not isinstance(n_restarts, int): + raise TypeError("n_restarts must be an integer.") + if n_restarts <= 0: + raise ValueError("n_restarts must be greater than zero.") + # Generate theta values using Latin hypercube sampling or Sobol sampling + # Generate theta values using Latin hypercube sampling + # Create a Latin Hypercube sampler that uses the dimensions of the theta names + sampler = scipy.stats.qmc.LatinHypercube(d=len(theta_names), seed=seed) + # Generate random samples in the range of [0, 1] for number of restarts + samples = sampler.random(n=n_restarts) + # Resulting samples should be size (n_restarts, len(theta_names)) + + elif multistart_sampling_method == "sobol_sampling": + if not isinstance(n_restarts, int): + raise TypeError("n_restarts must be an integer.") + if n_restarts <= 0: + raise ValueError("n_restarts must be greater than zero.") + sampler = scipy.stats.qmc.Sobol(d=len(theta_names), seed=seed) + # Generate theta values using Sobol sampling + # The first value of the Sobol sequence is 0, so we skip it + samples = sampler.random(n=n_restarts + 1)[1:] + + elif multistart_sampling_method == "user_provided_values": + raise ValueError( + "multistart_sampling_method='user_provided_values' requires user_provided_df." + ) - if len(return_values) > 0: - var_values = [] - if len(scen_names) > 1: # multiple scenarios - block_objects = self.ef_instance.component_objects( - Block, descend_into=False - ) - else: # single scenario - block_objects = [self.ef_instance] - for exp_i in block_objects: - vals = {} - for var in return_values: - exp_i_var = exp_i.find_component(str(var)) - if ( - exp_i_var is None - ): # we might have a block such as _mpisppy_data - continue - # if value to return is ContinuousSet - if type(exp_i_var) == ContinuousSet: - temp = list(exp_i_var) - else: - temp = [pyo.value(_) for _ in exp_i_var.values()] - if len(temp) == 1: - vals[var] = temp[0] - else: - vals[var] = temp - if len(vals) > 0: - var_values.append(vals) - var_values = pd.DataFrame(var_values) - if calc_cov is not NOTSET and calc_cov: - return obj_val, theta_vals, var_values, cov - elif calc_cov is NOTSET or not calc_cov: - return obj_val, theta_vals, var_values + else: + raise ValueError( + "Invalid sampling method. Choose 'uniform_random', 'latin_hypercube', 'sobol_sampling' or 'user_provided_values'." + ) + + if ( + multistart_sampling_method == "sobol_sampling" + or multistart_sampling_method == "latin_hypercube" + ): + # Scale the samples to the range of the lower and upper bounds for each theta in theta_names + # The samples are in the range [0, 1], so we scale them to the range of the lower and upper bounds + theta_vals_multistart = np.array( + [lower_bound + (upper_bound - lower_bound) * theta for theta in samples] + ) + + # Create a DataFrame where each row is an initial theta vector for a restart + if user_provided_df is None: + # Ensure theta_vals_multistart is 2D (n_restarts, len(theta_names)) + arr = np.atleast_2d(theta_vals_multistart) + if arr.shape[0] == 1 and n_restarts > 1: + arr = np.tile(arr, (n_restarts, 1)) + df_multistart = pd.DataFrame(arr, columns=theta_names) + + theta_arr = df_multistart[theta_names].to_numpy(dtype=float) + if not np.isfinite(theta_arr).all(): + raise ValueError("Initial theta values must be finite.") + if np.any(theta_arr < lower_bound) or np.any(theta_arr > upper_bound): + raise ValueError("Initial theta values must be within model bounds.") + + # Add columns for output info, initialized as nan + for name in theta_names: + df_multistart[f'converged_{name}'] = np.nan + df_multistart["final objective"] = np.nan + df_multistart["solver termination"] = "" + df_multistart["solve_time"] = np.nan + + # Debugging output + # print(df_multistart) + + return df_multistart + + # Redesigned _Q_opt method using scenario blocks, and combined with + # _Q_at_theta structure. + def _Q_opt( + self, + return_values=None, + bootlist=None, + solver="ef_ipopt", + theta_vals=None, + fix_theta=False, + multistart=False, + fixed_theta_values=None, + ): + """ + Solve the EF parameter-estimation problem and return objective/theta data. + + This routine creates the EF model via ``_create_scenario_blocks``, + solves it with the requested solver, and returns objective value plus + theta estimates. Depending on mode, it can also return variable values + and covariance estimates. + + Parameters + ---------- + return_values : list, optional + List of variable names to return values for. Default is None. + bootlist : list, optional + List of bootstrap experiment numbers to use. If None, use all experiments in exp_list. + Default is None. + theta_vals : dict, optional + Dictionary of theta values to set in the model. If None, use default values from experiment class. + Default is None. + solver : str, optional + Solver to use for optimization. Default is "ef_ipopt". + calc_cov : bool, optional + If True, calculate covariance matrix of estimated parameters. Default is NOTSET. + cov_n : int, optional + Number of data points to use for covariance calculation. Required if calc_cov is True. Default is NOTSET. + fix_theta : bool, optional + If True, fix the theta values in the model. If False, leave them free. + Default is False. + multistart : bool, optional + If True, run in multistart mode. Non-optimal termination is + returned instead of raising assertion failure. + + Returns + ------- + tuple + Return shape depends on mode: + 1. Standard solve: ``(obj_value, theta_series)`` + 2. Standard + return values: ``(obj_value, theta_series, var_values)`` + 3. Standard + covariance: ``(obj_value, theta_series, cov)`` or + ``(obj_value, theta_series, var_values, cov)`` + 4. Fixed-theta or multistart: ``(obj_value, theta_dict, worst_status)`` + """ + # Create extended form model with scenario blocks + model = self._create_scenario_blocks( + bootlist=bootlist, + theta_vals=theta_vals, + fix_theta=fix_theta, + fixed_theta_values=fixed_theta_values, + multistart=multistart, + ) + expanded_theta_names = list(model._parmest_theta_names) - if calc_cov is not NOTSET and calc_cov: - return obj_val, theta_vals, cov - elif calc_cov is NOTSET or not calc_cov: - return obj_val, theta_vals + # Print model if in diagnostic mode + if self.diagnostic_mode: + print("Parmest _Q_opt model with scenario blocks:") + model.pprint() + # Check solver and set options + if solver == "k_aug": + raise RuntimeError("k_aug no longer supported.") + if solver == "ef_ipopt": + sol = SolverFactory('ipopt') else: raise RuntimeError("Unknown solver in Q_Opt=" + solver) + # Currently, parmest is only tested with ipopt via ef_ipopt + # No other pyomo solvers have been verified to work with parmest from current release + # to my knowledge. + + if self.solver_options is not None: + for key in self.solver_options: + sol.options[key] = self.solver_options[key] + + # Solve model + solve_result = sol.solve(model, tee=self.tee) + partial_fix_mode = bool(fixed_theta_values) + + # Separate handling of termination conditions for _Q_at_theta vs _Q_opt + # If not fixing theta, ensure optimal termination of the solve to return result + if not fix_theta and not multistart and not partial_fix_mode: + # Ensure optimal termination + assert_optimal_termination(solve_result) + # If fixing theta, capture termination condition if not optimal unless infeasible + else: + # Initialize worst_status to optimal, update if not optimal + worst_status = pyo.TerminationCondition.optimal + # Get termination condition from solve result + status = solve_result.solver.termination_condition + + # In case of fixing theta, just log a warning if not optimal + if status != pyo.TerminationCondition.optimal: + # logger.warning( + # "Solver did not terminate optimally when thetas were fixed. " + # "Termination condition: %s", + # str(status), + # ) + # Unless infeasible, update worst_status + if worst_status != pyo.TerminationCondition.infeasible: + worst_status = status + + # Extract objective value + obj_value = pyo.value(model.Obj) + theta_estimates = {} + # Extract theta estimates from parent model + for name in expanded_theta_names: + theta_estimates[name] = pyo.value(model.parmest_theta[name]) + + self.obj_value = obj_value + self.estimated_theta = theta_estimates + + # If fixing theta, return objective value, theta estimates, and worst status + if fix_theta or multistart or partial_fix_mode: + return obj_value, theta_estimates, worst_status + + # Return theta estimates as a pandas Series + theta_estimates = pd.Series(theta_estimates) + + # Extract return values if requested + # Assumes the model components are named the same in each block, and are pyo.Vars. + if return_values is not None and len(return_values) > 0: + var_values = [] + # In the scenario blocks structure, exp_scenarios is an IndexedBlock + exp_blocks = self.ef_instance.exp_scenarios.values() + # Loop over each experiment block and extract requested variable values + for exp_i in exp_blocks: + # In each block, extract requested variables + vals = {} + for var in return_values: + # Find the variable in the block + exp_i_var = exp_i.find_component(str(var)) + # Check if variable exists in the block + if exp_i_var is None: + continue + # Extract value(s) from variable + if type(exp_i_var) == ContinuousSet: + temp = list(exp_i_var) + else: + temp = [pyo.value(_) for _ in exp_i_var.values()] + if len(temp) == 1: + vals[var] = temp[0] + else: + vals[var] = temp + # Only append if vals is not empty + if len(vals) > 0: + var_values.append(vals) + # Convert to DataFrame + var_values = pd.DataFrame(var_values) + + if return_values is not None and len(return_values) > 0: + return obj_value, theta_estimates, var_values + else: + return obj_value, theta_estimates + + # Removed old _Q_opt function def _cov_at_theta(self, method, solver, step): """ @@ -1181,14 +1523,20 @@ def _cov_at_theta(self, method, solver, step): cov : pd.DataFrame Covariance matrix of the estimated parameters """ + if hasattr(self.ef_instance, "exp_scenarios"): + ref_model = self.ef_instance.exp_scenarios[0] + else: + ref_model = self.ef_instance + if method == CovarianceMethod.reduced_hessian.value: # compute the inverse reduced hessian to be used # in the "reduced_hessian" method - # parmest makes the fitted parameters stage 1 variables + # retrieve the independent variables (i.e., estimated parameters) ind_vars = [] - for nd_name, Var, sol_val in ef_nonants(self.ef_instance): - ind_vars.append(Var) - # calculate the reduced hessian + for name in self.ef_instance._parmest_theta_names: + var = self.ef_instance.parmest_theta[name] + ind_vars.append(var) + solve_result, inv_red_hes = ( inverse_reduced_hessian.inv_reduced_hessian_barrier( self.ef_instance, @@ -1199,6 +1547,43 @@ def _cov_at_theta(self, method, solver, step): ) self.inv_red_hes = inv_red_hes + else: + # if not using the 'reduced_hessian' method, calculate the sum of squared errors + # using 'finite_difference' method or 'automatic_differentiation_kaug' + sse_vals = [] + for experiment in self.exp_list: + model = _get_labeled_model(experiment) + + # fix the value of the unknown parameters to the estimated values + for param in model.unknown_parameters: + param.fix(self.estimated_theta[param.name]) + + # re-solve the model with the estimated parameters + results = pyo.SolverFactory(solver).solve(model, tee=self.tee) + assert_optimal_termination(results) + + # choose and evaluate the sum of squared errors expression + if self.obj_function == ObjectiveType.SSE: + sse_expr = SSE(model) + elif self.obj_function == ObjectiveType.SSE_weighted: + sse_expr = SSE_weighted(model) + else: + raise ValueError( + f"Invalid objective function for covariance calculation. " + f"The covariance matrix can only be calculated using the built-in " + f"objective functions: {[e.value for e in ObjectiveType]}. Supply " + f"the Estimator object one of these built-in objectives and " + f"re-run the code." + ) + + # evaluate the numerical SSE and store it + sse_val = pyo.value(sse_expr) + sse_vals.append(sse_val) + + sse = sum(sse_vals) + logger.info( + f"The sum of squared errors at the estimated parameter(s) is: {sse}" + ) # Number of data points considered n = self.number_exp @@ -1206,42 +1591,6 @@ def _cov_at_theta(self, method, solver, step): # Extract the number of fitted parameters l = len(self.estimated_theta) - # calculate the sum of squared errors at the estimated parameter values - sse_vals = [] - for experiment in self.exp_list: - model = _get_labeled_model(experiment) - - # fix the value of the unknown parameters to the estimated values - for param in model.unknown_parameters: - param.fix(self.estimated_theta[param.name]) - - # re-solve the model with the estimated parameters - results = pyo.SolverFactory(solver).solve(model, tee=self.tee) - assert_optimal_termination(results) - - # choose and evaluate the sum of squared errors expression - if self.obj_function == ObjectiveType.SSE: - sse_expr = SSE(model) - elif self.obj_function == ObjectiveType.SSE_weighted: - sse_expr = SSE_weighted(model) - else: - raise ValueError( - f"Invalid objective function for covariance calculation. " - f"The covariance matrix can only be calculated using the built-in " - f"objective functions: {[e.value for e in ObjectiveType]}. Supply " - f"the Estimator object one of these built-in objectives and " - f"re-run the code." - ) - - # evaluate the numerical SSE and store it - sse_val = pyo.value(sse_expr) - sse_vals.append(sse_val) - - sse = sum(sse_vals) - logger.info( - f"The sum of squared errors at the estimated parameter(s) is: {sse}" - ) - """Calculate covariance assuming experimental observation errors are independent and follow a Gaussian distribution with constant variance. @@ -1264,11 +1613,11 @@ def _cov_at_theta(self, method, solver, step): # check if the user specified 'SSE' or 'SSE_weighted' as the objective function if self.obj_function == ObjectiveType.SSE: # check if the user defined the 'measurement_error' attribute - if hasattr(model, "measurement_error"): + if hasattr(ref_model, "measurement_error"): # get the measurement errors meas_error = [ - model.measurement_error[y_hat] - for y_hat, y in model.experiment_outputs.items() + ref_model.measurement_error[y_hat] + for y_hat, y in ref_model.experiment_outputs.items() ] # check if the user supplied the values of the measurement errors @@ -1340,10 +1689,10 @@ def _cov_at_theta(self, method, solver, step): ) elif self.obj_function == ObjectiveType.SSE_weighted: # check if the user defined the 'measurement_error' attribute - if hasattr(model, "measurement_error"): + if hasattr(ref_model, "measurement_error"): meas_error = [ - model.measurement_error[y_hat] - for y_hat, y in model.experiment_outputs.items() + ref_model.measurement_error[y_hat] + for y_hat, y in ref_model.experiment_outputs.items() ] # check if the user supplied the values for the measurement errors @@ -1380,241 +1729,57 @@ def _cov_at_theta(self, method, solver, step): raise AttributeError( 'Experiment model does not have suffix "measurement_error".' ) + else: + raise ValueError( + f"Invalid objective function for covariance calculation. " + f"The covariance matrix can only be calculated using the built-in " + f"objective functions: {[e.value for e in ObjectiveType]}. Supply " + f"the Estimator object one of these built-in objectives and " + f"re-run the code." + ) return cov - def _Q_at_theta(self, thetavals, initialize_parmest_model=False): + def _get_sample_list(self, samplesize, num_samples, replacement=True): + samplelist = list() + + scenario_numbers = list(range(len(self.exp_list))) + + if num_samples is None: + # This could get very large + for i, l in enumerate(combinations(scenario_numbers, samplesize)): + samplelist.append((i, np.sort(l))) + else: + for i in range(num_samples): + attempts = 0 + unique_samples = 0 # check for duplicates in each sample + duplicate = False # check for duplicates between samples + while (unique_samples <= len(self._return_theta_names())) and ( + not duplicate + ): + sample = np.random.choice( + scenario_numbers, samplesize, replace=replacement + ) + sample = np.sort(sample).tolist() + unique_samples = len(np.unique(sample)) + if sample in samplelist: + duplicate = True + + attempts += 1 + if attempts > num_samples: # arbitrary timeout limit + raise RuntimeError("""Internal error: timeout constructing + a sample, the dim of theta may be too + close to the samplesize""") + + samplelist.append((i, sample)) + + return samplelist + + def theta_est( + self, solver="ef_ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET + ): """ - Return the objective function value with fixed theta values. - - Parameters - ---------- - thetavals: dict - A dictionary of theta values. - - initialize_parmest_model: boolean - If True: Solve square problem instance, build extensive form of the model for - parameter estimation, and set flag model_initialized to True. Default is False. - - Returns - ------- - objectiveval: float - The objective function value. - thetavals: dict - A dictionary of all values for theta that were input. - solvertermination: Pyomo TerminationCondition - Tries to return the "worst" solver status across the scenarios. - pyo.TerminationCondition.optimal is the best and - pyo.TerminationCondition.infeasible is the worst. - """ - - optimizer = pyo.SolverFactory('ipopt') - - if len(thetavals) > 0: - dummy_cb = { - "callback": self._instance_creation_callback, - "ThetaVals": thetavals, - "theta_names": self._return_theta_names(), - "cb_data": None, - } - else: - dummy_cb = { - "callback": self._instance_creation_callback, - "theta_names": self._return_theta_names(), - "cb_data": None, - } - - if self.diagnostic_mode: - if len(thetavals) > 0: - print(' Compute objective at theta = ', str(thetavals)) - else: - print(' Compute objective at initial theta') - - # start block of code to deal with models with no constraints - # (ipopt will crash or complain on such problems without special care) - instance = _experiment_instance_creation_callback("FOO0", None, dummy_cb) - try: # deal with special problems so Ipopt will not crash - first = next(instance.component_objects(pyo.Constraint, active=True)) - active_constraints = True - except: - active_constraints = False - # end block of code to deal with models with no constraints - - WorstStatus = pyo.TerminationCondition.optimal - totobj = 0 - scenario_numbers = list(range(len(self.exp_list))) - if initialize_parmest_model: - # create dictionary to store pyomo model instances (scenarios) - scen_dict = dict() - - for snum in scenario_numbers: - sname = "scenario_NODE" + str(snum) - instance = _experiment_instance_creation_callback(sname, None, dummy_cb) - model_theta_names = self._expand_indexed_unknowns(instance) - - if initialize_parmest_model: - # list to store fitted parameter names that will be unfixed - # after initialization - theta_init_vals = [] - # use appropriate theta_names member - theta_ref = model_theta_names - - for i, theta in enumerate(theta_ref): - # Use parser in ComponentUID to locate the component - var_cuid = ComponentUID(theta) - var_validate = var_cuid.find_component_on(instance) - if var_validate is None: - logger.warning( - "theta_name %s was not found on the model", (theta) - ) - else: - try: - if len(thetavals) == 0: - var_validate.fix() - else: - var_validate.fix(thetavals[theta]) - theta_init_vals.append(var_validate) - except: - logger.warning( - 'Unable to fix model parameter value for %s (not a Pyomo model Var)', - (theta), - ) - - if active_constraints: - if self.diagnostic_mode: - print(' Experiment = ', snum) - print(' First solve with special diagnostics wrapper') - status_obj, solved, iters, time, regu = ( - utils.ipopt_solve_with_stats( - instance, optimizer, max_iter=500, max_cpu_time=120 - ) - ) - print( - " status_obj, solved, iters, time, regularization_stat = ", - str(status_obj), - str(solved), - str(iters), - str(time), - str(regu), - ) - - results = optimizer.solve(instance) - if self.diagnostic_mode: - print( - 'standard solve solver termination condition=', - str(results.solver.termination_condition), - ) - - if ( - results.solver.termination_condition - != pyo.TerminationCondition.optimal - ): - # DLW: Aug2018: not distinguishing "middlish" conditions - if WorstStatus != pyo.TerminationCondition.infeasible: - WorstStatus = results.solver.termination_condition - if initialize_parmest_model: - if self.diagnostic_mode: - print( - "Scenario {:d} infeasible with initialized parameter values".format( - snum - ) - ) - else: - if initialize_parmest_model: - if self.diagnostic_mode: - print( - "Scenario {:d} initialization successful with initial parameter values".format( - snum - ) - ) - if initialize_parmest_model: - # unfix parameters after initialization - for theta in theta_init_vals: - theta.unfix() - scen_dict[sname] = instance - else: - if initialize_parmest_model: - # unfix parameters after initialization - for theta in theta_init_vals: - theta.unfix() - scen_dict[sname] = instance - - objobject = getattr(instance, self._second_stage_cost_exp) - objval = pyo.value(objobject) - totobj += objval - - retval = totobj / len(scenario_numbers) # -1?? - if initialize_parmest_model and not hasattr(self, 'ef_instance'): - # create extensive form of the model using scenario dictionary - if len(scen_dict) > 0: - for scen in scen_dict.values(): - scen._mpisppy_probability = 1 / len(scen_dict) - - if use_mpisppy: - EF_instance = sputils._create_EF_from_scen_dict( - scen_dict, - EF_name="_Q_at_theta", - # suppress_warnings=True - ) - else: - EF_instance = local_ef._create_EF_from_scen_dict( - scen_dict, EF_name="_Q_at_theta", nonant_for_fixed_vars=True - ) - - self.ef_instance = EF_instance - # set self.model_initialized flag to True to skip extensive form model - # creation using theta_est() - self.model_initialized = True - - # return initialized theta values - if len(thetavals) == 0: - # use appropriate theta_names member - theta_ref = self._return_theta_names() - for i, theta in enumerate(theta_ref): - thetavals[theta] = theta_init_vals[i]() - - return retval, thetavals, WorstStatus - - def _get_sample_list(self, samplesize, num_samples, replacement=True): - samplelist = list() - - scenario_numbers = list(range(len(self.exp_list))) - - if num_samples is None: - # This could get very large - for i, l in enumerate(combinations(scenario_numbers, samplesize)): - samplelist.append((i, np.sort(l))) - else: - for i in range(num_samples): - attempts = 0 - unique_samples = 0 # check for duplicates in each sample - duplicate = False # check for duplicates between samples - while (unique_samples <= len(self._return_theta_names())) and ( - not duplicate - ): - sample = np.random.choice( - scenario_numbers, samplesize, replace=replacement - ) - sample = np.sort(sample).tolist() - unique_samples = len(np.unique(sample)) - if sample in samplelist: - duplicate = True - - attempts += 1 - if attempts > num_samples: # arbitrary timeout limit - raise RuntimeError("""Internal error: timeout constructing - a sample, the dim of theta may be too - close to the samplesize""") - - samplelist.append((i, sample)) - - return samplelist - - def theta_est( - self, solver="ef_ipopt", return_values=[], calc_cov=NOTSET, cov_n=NOTSET - ): - """ - Parameter estimation using all scenarios in the data + Parameter estimation using all scenarios in the data Parameters ---------- @@ -1671,13 +1836,7 @@ def theta_est( solver=solver, return_values=return_values ) - return self._Q_opt( - solver=solver, - return_values=return_values, - bootlist=None, - calc_cov=calc_cov, - cov_n=cov_n, - ) + return self._Q_opt(solver=solver, return_values=return_values, bootlist=None) def cov_est(self, method="finite_difference", solver="ipopt", step=1e-3): """ @@ -1944,9 +2103,398 @@ def leaveNout_bootstrap_test( return results + def theta_est_multistart( + self, + n_restarts=20, + multistart_sampling_method="uniform_random", + user_provided_df=None, + seed=None, + save_results=False, + file_name="multistart_results.csv", + ): + """ + Run multistart parameter estimation and aggregate per-start results. + + A canonical starts/results table is created first, then each start is + solved (potentially in parallel with ``ParallelTaskManager``), and the + output table is populated with objective values, solver terminations, + solve times, and converged theta values. + + Parameters + ---------- + n_restarts : int, optional + Number of starts for sampled methods. Ignored when + ``user_provided_df`` is provided. + multistart_sampling_method : str, optional + Sampling method for generated starts. + user_provided_df : DataFrame, optional + User-provided starts. If provided, these rows define the restart + set directly. + seed : int, optional + Seed used by sampling methods. + save_results : bool, optional + If True, write the full results DataFrame to ``file_name``. + file_name : str, optional + Output CSV path used when ``save_results`` is True. + + Returns + ------- + tuple + ``(results_df, best_theta, best_obj)``, where: + - ``results_df`` contains one row per start plus converged metadata + - ``best_theta`` is the selected best feasible theta dictionary or + ``None`` if no finite objective exists + - ``best_obj`` is the selected objective value or ``np.nan`` + + Notes + ----- + Best-run selection prioritizes acceptable solver terminations + (optimal/locallyOptimal/globallyOptimal) and then minimizes objective. + If no acceptable statuses exist, finite-objective rows are considered. + """ + if self.pest_deprecated is not None: + raise RuntimeError( + "Multistart is not supported in the deprecated parmest interface." + ) + if user_provided_df is None: + if not isinstance(n_restarts, int): + raise TypeError("n_restarts must be an integer.") + if n_restarts <= 0: + raise ValueError("n_restarts must be greater than zero.") + + n_restarts_for_generation = None if user_provided_df is not None else n_restarts + results_df = self._generate_initial_theta( + seed=seed, + n_restarts=n_restarts_for_generation, + multistart_sampling_method=multistart_sampling_method, + user_provided_df=user_provided_df, + experiment_number=0, + ) + theta_names = [ + c + for c in results_df.columns + if not c.startswith("converged_") + and c not in {"final objective", "solver termination", "solve_time"} + ] + n_restarts = results_df.shape[0] + + # Convert each row to (row_index, theta_dict) + tasks = [] + for i in range(results_df.shape[0]): + Theta = {name: float(results_df.iloc[i][name]) for name in theta_names} + tasks.append((i, Theta)) + + task_mgr = utils.ParallelTaskManager(len(tasks)) + local_tasks = task_mgr.global_to_local_data(tasks) + + # Solve in parallel + local_results = [] + for i, Theta in local_tasks: + import time + + t0 = time.time() + try: + final_obj, theta_hat, worst = self._Q_opt( + theta_vals=Theta, multistart=True + ) + solve_time = time.time() - t0 + local_results.append((i, final_obj, str(worst), solve_time, theta_hat)) + except Exception as exc: + solve_time = time.time() - t0 + local_results.append( + ( + i, + np.nan, + f"exception(start={i}, sampler={multistart_sampling_method}): {exc}", + solve_time, + None, + ) + ) + + global_results = task_mgr.allgather_global_data(local_results) + + # Fill results_df + for i, final_obj, term, solve_time, theta_hat in global_results: + results_df.at[i, "final objective"] = final_obj + results_df.at[i, "solver termination"] = term + results_df.at[i, "solve_time"] = solve_time + + if theta_hat is not None: + for name in theta_names: + if name in theta_hat: + results_df.at[i, f"converged_{name}"] = float(theta_hat[name]) + + # Best solution: + # prioritize starts with acceptable solver terminations, then minimum objective. + acceptable_terms = { + str(pyo.TerminationCondition.optimal), + str(pyo.TerminationCondition.locallyOptimal), + str(pyo.TerminationCondition.globallyOptimal), + } + finite_obj_mask = np.isfinite( + results_df["final objective"].to_numpy(dtype=float) + ) + acceptable_mask = results_df["solver termination"].isin(acceptable_terms) + ranked = results_df[finite_obj_mask & acceptable_mask] + if ranked.empty: + ranked = results_df[finite_obj_mask] + + if ranked.empty: + best_theta = None + best_obj = np.nan + else: + best_idx = ranked["final objective"].astype(float).idxmin() + best_obj = float(results_df.loc[best_idx, "final objective"]) + best_theta = { + name: float(results_df.loc[best_idx, f"converged_{name}"]) + for name in theta_names + } + + if save_results: + results_df.to_csv(file_name, index=False) + + return results_df, best_theta, best_obj + + def _build_profile_grid(self, profiled_theta, grid, n_grid, theta_hat): + """Build sorted profile grid for a single profiled parameter.""" + if not isinstance(profiled_theta, str): + raise TypeError("profiled_theta must be a string.") + if grid is not None and not isinstance( + grid, (list, tuple, np.ndarray, pd.Series) + ): + raise TypeError("grid must be a sequence of numeric values or None.") + if not isinstance(n_grid, int): + raise TypeError("n_grid must be an integer.") + if n_grid < 2: + raise ValueError("n_grid must be at least 2.") + if not isinstance(theta_hat, dict): + raise TypeError("theta_hat must be a dictionary.") + if profiled_theta not in theta_hat: + raise ValueError(f"theta_hat does not include '{profiled_theta}'.") + + template_model = self._create_parmest_model(0) + theta_var = template_model.find_component(profiled_theta) + if theta_var is None: + raise ValueError(f"Unknown theta name '{profiled_theta}'.") + + theta_hat_val = float(theta_hat[profiled_theta]) + if grid is None: + if theta_var.lb is None or theta_var.ub is None: + raise ValueError( + f"Cannot auto-build grid for '{profiled_theta}' without lower and upper bounds." + ) + values = np.linspace(float(theta_var.lb), float(theta_var.ub), n_grid) + else: + values = np.asarray(list(grid), dtype=float) + if values.size == 0: + raise ValueError("Provided grid must contain at least one value.") + + values = np.append(values, theta_hat_val) + values = np.unique(np.round(values.astype(float), decimals=14)) + values.sort() + return values + + def _order_grid_for_continuation(self, grid_values, center): + """Return center-out ordering for continuation warm starts.""" + ordered = sorted( + [float(v) for v in grid_values], key=lambda v: (abs(v - center), v) + ) + return ordered + + def profile_likelihood( + self, + profiled_theta, + grid=None, + n_grid=21, + theta_hat=None, + obj_hat=None, + use_multistart_for_baseline=False, + baseline_multistart_kwargs=None, + solver="ef_ipopt", + warmstart="neighbor", + max_consecutive_failures=None, + return_theta_paths=True, + seed=None, + ): + """ + Compute one-dimensional profile likelihood curves by fixing one parameter + at grid values and re-optimizing all other parameters. + """ + if self.pest_deprecated is not None: + raise RuntimeError( + "Profile likelihood is not supported in the deprecated parmest interface." + ) + if isinstance(profiled_theta, str): + profiled_theta_list = [profiled_theta] + elif isinstance(profiled_theta, (list, tuple)): + profiled_theta_list = list(profiled_theta) + else: + raise TypeError( + "profiled_theta must be a string or a list/tuple of strings." + ) + if len(profiled_theta_list) == 0: + raise ValueError("At least one profiled theta must be provided.") + if not all(isinstance(p, str) for p in profiled_theta_list): + raise TypeError("Each profiled theta name must be a string.") + if warmstart not in ("neighbor", "none"): + raise ValueError("warmstart must be either 'neighbor' or 'none'.") + if max_consecutive_failures is not None: + if not isinstance(max_consecutive_failures, int): + raise TypeError("max_consecutive_failures must be an integer or None.") + if max_consecutive_failures <= 0: + raise ValueError("max_consecutive_failures must be greater than zero.") + if not isinstance(return_theta_paths, bool): + raise TypeError("return_theta_paths must be a bool.") + if seed is not None and not isinstance(seed, int): + raise TypeError("seed must be an integer or None.") + + theta_names = self._expand_indexed_unknowns(self._create_parmest_model(0)) + unknown_requested = set(profiled_theta_list).difference(theta_names) + if unknown_requested: + raise ValueError( + f"Unknown profile theta name(s): {sorted(unknown_requested)}. " + f"Known names: {theta_names}." + ) + + if seed is not None: + np.random.seed(seed) + + if theta_hat is None or obj_hat is None: + if use_multistart_for_baseline: + ms_kwargs = dict(baseline_multistart_kwargs or {}) + _, best_theta, best_obj = self.theta_est_multistart(**ms_kwargs) + if best_theta is None or not np.isfinite(best_obj): + raise RuntimeError( + "Failed to compute baseline from multistart: no feasible solution found." + ) + theta_hat = best_theta + obj_hat = best_obj + else: + obj_hat, theta_hat_series = self.theta_est(solver=solver) + theta_hat = theta_hat_series.to_dict() + + theta_hat = {k: float(v) for k, v in theta_hat.items()} + obj_hat = float(obj_hat) + + rows = [] + acceptable_terms = { + str(pyo.TerminationCondition.optimal), + str(pyo.TerminationCondition.locallyOptimal), + str(pyo.TerminationCondition.globallyOptimal), + } + started_at = time.time() + + for pname in profiled_theta_list: + if isinstance(grid, dict): + grid_spec = grid.get(pname, None) + else: + grid_spec = grid + grid_values = self._build_profile_grid( + profiled_theta=pname, grid=grid_spec, n_grid=n_grid, theta_hat=theta_hat + ) + + if warmstart == "neighbor": + ordered_grid = self._order_grid_for_continuation( + grid_values, center=theta_hat[pname] + ) + else: + ordered_grid = [float(v) for v in grid_values] + + previous_theta = dict(theta_hat) + consecutive_failures = 0 + + for g in ordered_grid: + init_theta = dict(theta_hat) + if warmstart == "neighbor": + init_theta.update(previous_theta) + init_theta[pname] = float(g) + + t0 = time.time() + status = "" + success = False + obj_val = np.nan + theta_conv = None + try: + obj_val, theta_conv, term = self._Q_opt( + theta_vals=init_theta, + fixed_theta_values={pname: float(g)}, + solver=solver, + ) + status = str(term) + success = status in acceptable_terms and np.isfinite(obj_val) + if success: + consecutive_failures = 0 + if warmstart == "neighbor": + previous_theta = dict(theta_conv) + else: + consecutive_failures += 1 + except Exception as exc: + status = f"exception: {exc}" + consecutive_failures += 1 + + row = { + "profiled_theta": pname, + "theta_value": float(g), + "obj": float(obj_val) if np.isfinite(obj_val) else np.nan, + "status": status, + "success": bool(success), + "solve_time": float(time.time() - t0), + } + if return_theta_paths and isinstance(theta_conv, dict): + for tname, tval in theta_conv.items(): + row[f"theta__{tname}"] = float(tval) + rows.append(row) + + if ( + max_consecutive_failures is not None + and consecutive_failures >= max_consecutive_failures + ): + break + + profiles = pd.DataFrame(rows) + if profiles.empty: + profiles = pd.DataFrame( + columns=[ + "profiled_theta", + "theta_value", + "obj", + "status", + "success", + "solve_time", + "delta_obj", + "lr_stat", + ] + ) + else: + profiles = profiles.sort_values( + by=["profiled_theta", "theta_value"], kind="stable" + ).reset_index(drop=True) + profiles["delta_obj"] = profiles["obj"] - obj_hat + profiles["lr_stat"] = 2.0 * profiles["delta_obj"] + + return { + "profiles": profiles, + "baseline": { + "theta_hat": theta_hat, + "obj_hat": obj_hat, + "solver": solver, + "used_multistart": bool(use_multistart_for_baseline), + }, + "metadata": { + "grid_strategy": "user_provided" if grid is not None else "auto_bounds", + "warmstart": warmstart, + "seed": seed, + "profiled_theta": list(profiled_theta_list), + "started_at_epoch": float(started_at), + "finished_at_epoch": float(time.time()), + }, + } + + # Updated version that uses _Q_opt def objective_at_theta(self, theta_values=None, initialize_parmest_model=False): """ - Objective value for each theta + Objective value for each theta, solving extensive form problem with + fixed theta values. Parameters ---------- @@ -1958,7 +2506,6 @@ def objective_at_theta(self, theta_values=None, initialize_parmest_model=False): of the model for parameter estimation, and set flag model_initialized to True. Default is False. - Returns ------- obj_at_theta: pd.DataFrame @@ -1973,65 +2520,73 @@ def objective_at_theta(self, theta_values=None, initialize_parmest_model=False): initialize_parmest_model=initialize_parmest_model, ) - if len(self.estimator_theta_names) == 0: - pass # skip assertion if model has no fitted parameters - else: - # create a local instance of the pyomo model to access model variables and parameters - model_temp = self._create_parmest_model(0) - model_theta_list = self._expand_indexed_unknowns(model_temp) - - # if self.estimator_theta_names is not the same as temp model_theta_list, - # create self.theta_names_updated - if set(self.estimator_theta_names) == set(model_theta_list) and len( - self.estimator_theta_names - ) == len(set(model_theta_list)): - pass - else: - self.theta_names_updated = model_theta_list + if initialize_parmest_model: + # Print deprecation warning, that this option will be removed in + # future releases. + deprecation_warning( + "The `initialize_parmest_model` option in `objective_at_theta()` is " + "deprecated and will be removed in future releases. Please ensure the" + "model is initialized within the Experiment class definition.", + version="6.9.5", + ) if theta_values is None: all_thetas = {} # dictionary to store fitted variables # use appropriate theta names member - theta_names = model_theta_list + # Get theta names from fresh parmest model, assuming this can be called + # directly after creating Estimator. + theta_names = self._expand_indexed_unknowns(self._create_parmest_model(0)) else: assert isinstance(theta_values, pd.DataFrame) # for parallel code we need to use lists and dicts in the loop theta_names = theta_values.columns # # check if theta_names are in model - for theta in list(theta_names): - theta_temp = theta.replace("'", "") # cleaning quotes from theta_names - assert theta_temp in [ - t.replace("'", "") for t in model_theta_list - ], "Theta name {} in 'theta_values' not in 'theta_names' {}".format( - theta_temp, model_theta_list + # Clean names, ignore quotes, and compare sets + clean_provided = [t.replace("'", "") for t in theta_names] + if len(clean_provided) != len(set(clean_provided)): + raise ValueError( + f"Duplicate theta names are not allowed: {clean_provided}" ) + clean_expected = [ + t.replace("'", "") + for t in self._expand_indexed_unknowns(self._create_parmest_model(0)) + ] + # If they do not match, raise error + if (len(clean_provided) != len(clean_expected)) or ( + set(clean_provided) != set(clean_expected) + ): + raise ValueError( + f"Provided theta values {clean_provided} do not match expected theta names {clean_expected}." + ) + # Rename columns using cleaned names + if list(clean_provided) != list(theta_names): + theta_values = theta_values.copy() + theta_values.columns = clean_provided - assert len(list(theta_names)) == len(model_theta_list) - + # Convert to list of dicts for parallel processing all_thetas = theta_values.to_dict('records') + # Initialize task manager + num_tasks = len(all_thetas) if all_thetas else 1 + task_mgr = utils.ParallelTaskManager(num_tasks) + + # Use local theta values for each task if all_thetas is provided, else empty list if all_thetas: - task_mgr = utils.ParallelTaskManager(len(all_thetas)) local_thetas = task_mgr.global_to_local_data(all_thetas) - else: - if initialize_parmest_model: - task_mgr = utils.ParallelTaskManager( - 1 - ) # initialization performed using just 1 set of theta values + elif initialize_parmest_model: + local_thetas = [] + # walk over the mesh, return objective function all_obj = list() if len(all_thetas) > 0: for Theta in local_thetas: - obj, thetvals, worststatus = self._Q_at_theta( - Theta, initialize_parmest_model=initialize_parmest_model + obj, thetvals, worststatus = self._Q_opt( + theta_vals=Theta, fix_theta=True ) if worststatus != pyo.TerminationCondition.infeasible: all_obj.append(list(Theta.values()) + [obj]) - # DLW, Aug2018: should we also store the worst solver status? else: - obj, thetvals, worststatus = self._Q_at_theta( - thetavals={}, initialize_parmest_model=initialize_parmest_model - ) + obj, thetvals, worststatus = self._Q_opt(theta_vals=None, fix_theta=True) if worststatus != pyo.TerminationCondition.infeasible: all_obj.append(list(thetvals.values()) + [obj]) diff --git a/pyomo/contrib/parmest/tests/test_graphics.py b/pyomo/contrib/parmest/tests/test_graphics.py index 42c0b2a9f1f..e9a7b1c75fd 100644 --- a/pyomo/contrib/parmest/tests/test_graphics.py +++ b/pyomo/contrib/parmest/tests/test_graphics.py @@ -61,6 +61,28 @@ def test_grouped_boxplot(self): def test_grouped_violinplot(self): graphics.grouped_violinplot(self.A, self.B) + def test_profile_plot_smoke_single_parameter(self): + prof = pd.DataFrame( + { + "profiled_theta": ["theta"] * 3, + "theta_value": [0.8, 1.0, 1.2], + "lr_stat": [3.0, 0.0, 4.0], + "success": [True, True, False], + } + ) + graphics.profile_likelihood_plot({"profiles": prof}, alpha=0.95) + + def test_profile_plot_smoke_multi_parameter(self): + prof = pd.DataFrame( + { + "profiled_theta": ["theta_a", "theta_a", "theta_b", "theta_b"], + "theta_value": [1.0, 2.0, -1.0, 0.0], + "lr_stat": [2.0, 0.0, 1.0, 0.0], + "success": [True, True, True, True], + } + ) + graphics.profile_likelihood_plot({"profiles": prof}, alpha=0.9) + if __name__ == '__main__': unittest.main() diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index dfeab769e71..dfd60b8704f 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -11,7 +11,6 @@ import os import subprocess from itertools import product - from pyomo.common.unittest import pytest from parameterized import parameterized, parameterized_class import pyomo.common.unittest as unittest @@ -30,6 +29,7 @@ pynumero_ASL_available = AmplInterface.available() testdir = this_file_dir() +# TESTS HERE WILL BE MODIFIED FOR _Q_OPT_BLOCKS LATER # Set the global seed for random number generation in tests _RANDOM_SEED_FOR_TESTING = 524 @@ -508,40 +508,6 @@ def test_parallel_parmest(self): retcode = subprocess.call(rlist) self.assertEqual(retcode, 0) - @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available") - def test_theta_est_cov(self): - objval, thetavals, cov = self.pest.theta_est(calc_cov=True, cov_n=6) - - self.assertAlmostEqual(objval, 4.3317112, places=2) - self.assertAlmostEqual( - thetavals["asymptote"], 19.1426, places=2 - ) # 19.1426 from the paper - self.assertAlmostEqual( - thetavals["rate_constant"], 0.5311, places=2 - ) # 0.5311 from the paper - - # Covariance matrix - self.assertAlmostEqual( - cov["asymptote"]["asymptote"], 6.155892, places=2 - ) # 6.22864 from paper - self.assertAlmostEqual( - cov["asymptote"]["rate_constant"], -0.425232, places=2 - ) # -0.4322 from paper - self.assertAlmostEqual( - cov["rate_constant"]["asymptote"], -0.425232, places=2 - ) # -0.4322 from paper - self.assertAlmostEqual( - cov["rate_constant"]["rate_constant"], 0.040571, places=2 - ) # 0.04124 from paper - - """ Why does the covariance matrix from parmest not match the paper? Parmest is - calculating the exact reduced Hessian. The paper (Rooney and Bielger, 2001) likely - employed the first order approximation common for nonlinear regression. The paper - values were verified with Scipy, which uses the same first order approximation. - The formula used in parmest was verified against equations (7-5-15) and (7-5-16) in - "Nonlinear Parameter Estimation", Y. Bard, 1974. - """ - def test_cov_scipy_least_squares_comparison(self): """ Scipy results differ in the 3rd decimal place from the paper. It is possible @@ -655,20 +621,27 @@ def setUp(self): columns=["hour", "y"], ) + # Updated models to use Vars for experiment output, and Constraints def rooney_biegler_params(data): model = pyo.ConcreteModel() model.asymptote = pyo.Param(initialize=15, mutable=True) model.rate_constant = pyo.Param(initialize=0.5, mutable=True) - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - def response_rule(m, h): - expr = m.asymptote * (1 - pyo.exp(-m.rate_constant * h)) - return expr + # Fix the experiment inputs + model.h.fix() - model.response_function = pyo.Expression(data.hour, rule=response_rule) + # Add experiment outputs + model.y = pyo.Var(initialize=data['y'].iloc[0], within=pyo.PositiveReals) + + # Define the model equations + def response_rule(m): + return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.h)) + + model.response_con = pyo.Constraint(rule=response_rule) return model @@ -683,14 +656,14 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + m.experiment_outputs.update([(m.y, self.data["y"])]) m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) m.unknown_parameters.update( (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) rooney_biegler_params_exp_list = [] for i in range(self.data.shape[0]): @@ -701,24 +674,30 @@ def label_model(self): def rooney_biegler_indexed_params(data): model = pyo.ConcreteModel() + # Define the indexed parameters model.param_names = pyo.Set(initialize=["asymptote", "rate_constant"]) model.theta = pyo.Param( model.param_names, initialize={"asymptote": 15, "rate_constant": 0.5}, mutable=True, ) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Fix the experiment inputs + model.h.fix() - def response_rule(m, h): - expr = m.theta["asymptote"] * ( - 1 - pyo.exp(-m.theta["rate_constant"] * h) - ) - return expr + # Add experiment outputs + model.y = pyo.Var(initialize=data['y'].iloc[0], within=pyo.PositiveReals) - model.response_function = pyo.Expression(data.hour, rule=response_rule) + # Define the model equations + def response_rule(m): + return m.y == m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * m.h) + ) + # Add the model equations to the model + model.response_con = pyo.Constraint(rule=response_rule) return model class RooneyBieglerExperimentIndexedParams(RooneyBieglerExperiment): @@ -732,13 +711,14 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + m.experiment_outputs.update([(m.y, self.data["y"])]) m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) m.unknown_parameters.update((k, pyo.ComponentUID(k)) for k in [m.theta]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + rooney_biegler_indexed_params_exp_list = [] for i in range(self.data.shape[0]): rooney_biegler_indexed_params_exp_list.append( @@ -753,14 +733,20 @@ def rooney_biegler_vars(data): model.asymptote.fixed = True # parmest will unfix theta variables model.rate_constant.fixed = True - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - def response_rule(m, h): - expr = m.asymptote * (1 - pyo.exp(-m.rate_constant * h)) - return expr + # Fix the experiment inputs + model.h.fix() - model.response_function = pyo.Expression(data.hour, rule=response_rule) + # Add experiment outputs + model.y = pyo.Var(initialize=data['y'].iloc[0], within=pyo.PositiveReals) + + # Define the model equations + def response_rule(m): + return m.y == m.asymptote * (1 - pyo.exp(-m.rate_constant * m.h)) + + model.response_con = pyo.Constraint(rule=response_rule) return model @@ -775,14 +761,14 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + m.experiment_outputs.update([(m.y, self.data["y"])]) m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) m.unknown_parameters.update( (k, pyo.ComponentUID(k)) for k in [m.asymptote, m.rate_constant] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) rooney_biegler_vars_exp_list = [] for i in range(self.data.shape[0]): @@ -802,16 +788,22 @@ def rooney_biegler_indexed_vars(data): ) model.theta["rate_constant"].fixed = True - model.hour = pyo.Param(within=pyo.PositiveReals, mutable=True) - model.y = pyo.Param(within=pyo.PositiveReals, mutable=True) + # Add the experiment inputs + model.h = pyo.Var(initialize=data["hour"].iloc[0], bounds=(0, 10)) - def response_rule(m, h): - expr = m.theta["asymptote"] * ( - 1 - pyo.exp(-m.theta["rate_constant"] * h) + # Fix the experiment inputs + model.h.fix() + + # Add experiment outputs + model.y = pyo.Var(initialize=data['y'].iloc[0], within=pyo.PositiveReals) + + # Define the model equations + def response_rule(m): + return m.y == m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * m.h) ) - return expr - model.response_function = pyo.Expression(data.hour, rule=response_rule) + model.response_con = pyo.Constraint(rule=response_rule) return model @@ -826,28 +818,21 @@ def label_model(self): m = self.model m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) - m.experiment_outputs.update( - [(m.hour, self.data["hour"]), (m.y, self.data["y"])] - ) + m.experiment_outputs.update([(m.y, self.data["y"])]) m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) m.unknown_parameters.update((k, pyo.ComponentUID(k)) for k in [m.theta]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + rooney_biegler_indexed_vars_exp_list = [] for i in range(self.data.shape[0]): rooney_biegler_indexed_vars_exp_list.append( RooneyBieglerExperimentIndexedVars(self.data.loc[i, :]) ) - # Sum of squared error function - def SSE(model): - expr = ( - model.experiment_outputs[model.y] - - model.response_function[model.experiment_outputs[model.hour]] - ) ** 2 - return expr - - self.objective_function = SSE + self.objective_function = 'SSE' theta_vals = pd.DataFrame([20, 1], index=["asymptote", "rate_constant"]).T theta_vals_index = pd.DataFrame( @@ -899,16 +884,16 @@ def check_rooney_biegler_results(self, objval, cov): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - cov.iloc[asymptote_index, asymptote_index], 6.30579403, places=2 + cov.iloc[asymptote_index, asymptote_index], 6.155892, places=2 ) # 6.22864 from paper self.assertAlmostEqual( - cov.iloc[asymptote_index, rate_constant_index], -0.4395341, places=2 + cov.iloc[asymptote_index, rate_constant_index], -0.425232, places=2 ) # -0.4322 from paper self.assertAlmostEqual( - cov.iloc[rate_constant_index, asymptote_index], -0.4395341, places=2 + cov.iloc[rate_constant_index, asymptote_index], -0.425232, places=2 ) # -0.4322 from paper self.assertAlmostEqual( - cov.iloc[rate_constant_index, rate_constant_index], 0.04193591, places=2 + cov.iloc[rate_constant_index, rate_constant_index], 0.040571, places=2 ) # 0.04124 from paper @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available') @@ -918,8 +903,13 @@ def test_parmest_basics(self): pest = parmest.Estimator( parmest_input["exp_list"], obj_function=self.objective_function ) - - objval, thetavals, cov = pest.theta_est(calc_cov=True, cov_n=6) + # estimate the parameters and covariance matrix + objval, thetavals = pest.theta_est() + # For covariance, using reduced_hessian method since finite difference + # and automatic differentiation may differ from paper results in the + # 3rd decimal place, likely due to differences in finite difference + # approximation of the Jacobian + cov = pest.cov_est(method="reduced_hessian") self.check_rooney_biegler_results(objval, cov) obj_at_theta = pest.objective_at_theta(parmest_input["theta_vals"]) @@ -933,7 +923,8 @@ def test_parmest_basics_with_initialize_parmest_model_option(self): parmest_input["exp_list"], obj_function=self.objective_function ) - objval, thetavals, cov = pest.theta_est(calc_cov=True, cov_n=6) + objval, thetavals = pest.theta_est() + cov = pest.cov_est(method="reduced_hessian") self.check_rooney_biegler_results(objval, cov) obj_at_theta = pest.objective_at_theta( @@ -954,7 +945,8 @@ def test_parmest_basics_with_square_problem_solve(self): parmest_input["theta_vals"], initialize_parmest_model=True ) - objval, thetavals, cov = pest.theta_est(calc_cov=True, cov_n=6) + objval, thetavals = pest.theta_est() + cov = pest.cov_est(method="reduced_hessian") self.check_rooney_biegler_results(objval, cov) self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) @@ -970,7 +962,8 @@ def test_parmest_basics_with_square_problem_solve_no_theta_vals(self): obj_at_theta = pest.objective_at_theta(initialize_parmest_model=True) - objval, thetavals, cov = pest.theta_est(calc_cov=True, cov_n=6) + objval, thetavals = pest.theta_est() + cov = pest.cov_est(method="reduced_hessian") self.check_rooney_biegler_results(objval, cov) @@ -1107,7 +1100,8 @@ def _dccrate(m, t): def ComputeFirstStageCost_rule(m): return 0 - m.FirstStageCost = pyo.Expression(rule=ComputeFirstStageCost_rule) + # Model objective component names adjusted to prevent reserved name error. + m.FirstStage = pyo.Expression(rule=ComputeFirstStageCost_rule) def ComputeSecondStageCost_rule(m): return sum( @@ -1117,14 +1111,12 @@ def ComputeSecondStageCost_rule(m): for t in meas_t ) - m.SecondStageCost = pyo.Expression(rule=ComputeSecondStageCost_rule) + m.SecondStage = pyo.Expression(rule=ComputeSecondStageCost_rule) def total_cost_rule(model): - return model.FirstStageCost + model.SecondStageCost + return model.FirstStage + model.SecondStage - m.Total_Cost_Objective = pyo.Objective( - rule=total_cost_rule, sense=pyo.minimize - ) + m.Total_Cost = pyo.Objective(rule=total_cost_rule, sense=pyo.minimize) disc = pyo.TransformationFactory("dae.collocation") disc.apply_to(m, nfe=20, ncp=2) @@ -1165,6 +1157,10 @@ def label_model(self): m.unknown_parameters.update( (k, pyo.ComponentUID(k)) for k in [m.k1, m.k2] ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update((m.ca[t], None) for t in meas_time_points) + m.measurement_error.update((m.cb[t], None) for t in meas_time_points) + m.measurement_error.update((m.cc[t], None) for t in meas_time_points) def get_labeled_model(self): self.create_model() @@ -1210,8 +1206,8 @@ def get_labeled_model(self): exp_list_df = [ReactorDesignExperimentDAE(data_df)] exp_list_dict = [ReactorDesignExperimentDAE(data_dict)] - self.pest_df = parmest.Estimator(exp_list_df) - self.pest_dict = parmest.Estimator(exp_list_dict) + self.pest_df = parmest.Estimator(exp_list_df, obj_function="SSE") + self.pest_dict = parmest.Estimator(exp_list_dict, obj_function="SSE") # Estimator object with multiple scenarios exp_list_df_multiple = [ @@ -1223,8 +1219,12 @@ def get_labeled_model(self): ReactorDesignExperimentDAE(data_dict), ] - self.pest_df_multiple = parmest.Estimator(exp_list_df_multiple) - self.pest_dict_multiple = parmest.Estimator(exp_list_dict_multiple) + self.pest_df_multiple = parmest.Estimator( + exp_list_df_multiple, obj_function="SSE" + ) + self.pest_dict_multiple = parmest.Estimator( + exp_list_dict_multiple, obj_function="SSE" + ) # Create an instance of the model self.m_df = ABC_model(data_df) @@ -1315,10 +1315,11 @@ def test_covariance(self): # 3 data components (ca, cb, cc), 20 timesteps, 1 scenario = 60 # In this example, this is the number of data points in data_df, but that's # only because the data is indexed by time and contains no additional information. - n = 60 + n = 20 # Compute covariance using parmest - obj, theta, cov = self.pest_df.theta_est(calc_cov=True, cov_n=n) + obj, theta = self.pest_df.theta_est() + cov = self.pest_df.cov_est(method="reduced_hessian") # Compute covariance using interior_point vars_list = [self.m_df.k1, self.m_df.k2] @@ -1416,6 +1417,131 @@ def test_theta_est_with_square_initialization_diagnostic_mode_true(self): self.pest.diagnostic_mode = False +class LinearThetaExperiment(Experiment): + def __init__(self, x, y, include_second_output=False): + self.x_data = x + self.y_data = y + self.include_second_output = include_second_output + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.theta = pyo.Var(initialize=0.0, bounds=(-10.0, 10.0)) + m.x = pyo.Param(initialize=float(self.x_data), mutable=False) + m.y = pyo.Var(initialize=float(self.y_data)) + m.y_link = pyo.Constraint(expr=m.y == m.theta + m.x) + if self.include_second_output: + m.z = pyo.Var(initialize=2.0 * self.y_data) + m.z_link = pyo.Constraint(expr=m.z == 2.0 * m.theta + m.x) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, float(self.y_data))]) + if self.include_second_output: + m.experiment_outputs.update([(m.z, float(2.0 * self.y_data))]) + + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + if self.include_second_output: + m.measurement_error.update([(m.z, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + +def _build_estimator(data, include_second_output=False): + exp_list = [ + LinearThetaExperiment(x=x, y=y, include_second_output=include_second_output) + for x, y in data + ] + return parmest.Estimator(exp_list, obj_function="SSE") + + +@unittest.skipIf( + not parmest.parmest_available, + "Cannot test parmest: required dependencies are missing", +) +class TestParmestBlockEF(unittest.TestCase): + def test_block_ef_structure_counts(self): + pest = _build_estimator([(1.0, 2.0), (2.0, 4.0)]) + model = pest._create_scenario_blocks() + + theta_names = model._parmest_theta_names + self.assertEqual(len(list(model.exp_scenarios.keys())), 2) + self.assertEqual(len(model.theta_link_constraints), 2 * len(theta_names)) + self.assertTrue(hasattr(model, "Obj")) + for block in model.exp_scenarios.values(): + self.assertFalse(block.Total_Cost_Objective.active) + + def test_block_isolation_no_component_leakage(self): + pest = _build_estimator([(1.0, 2.0), (5.0, 6.0)]) + model = pest._create_scenario_blocks() + + block0 = model.exp_scenarios[0] + block1 = model.exp_scenarios[1] + self.assertIsNot(block0.y, block1.y) + block0.y.set_value(123.0) + self.assertNotEqual(pyo.value(block1.y), 123.0) + self.assertNotEqual(pyo.value(block0.x), pyo.value(block1.x)) + + def test_fix_theta_sets_all_scenario_theta_values(self): + pest = _build_estimator([(1.0, 2.0), (2.0, 4.0)]) + model = pest._create_scenario_blocks(theta_vals={"theta": 1.0}, fix_theta=True) + + self.assertTrue(model.parmest_theta["theta"].fixed) + self.assertAlmostEqual(pyo.value(model.parmest_theta["theta"]), 1.0, places=10) + for block in model.exp_scenarios.values(): + self.assertTrue(block.theta.fixed) + self.assertAlmostEqual(pyo.value(block.theta), 1.0, places=10) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_objective_at_theta_fixed_value(self): + pest = _build_estimator([(1.0, 2.0), (2.0, 4.0)]) + theta_values = pd.DataFrame([[1.0]], columns=["theta"]) + obj_at_theta = pest.objective_at_theta(theta_values=theta_values) + # residuals at theta=1 are [0, 1], objective is averaged over two scenarios + self.assertAlmostEqual(obj_at_theta.loc[0, "obj"], 0.5, places=8) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_objective_at_theta_none_uses_initial_theta(self): + pest = _build_estimator([(1.0, 2.0), (2.0, 3.0)]) + obj_at_theta = pest.objective_at_theta() + # with theta initialized to 0, predictions are [1,2], residuals [1,1], avg objective 1 + self.assertAlmostEqual(obj_at_theta.loc[0, "obj"], 1.0, places=8) + self.assertAlmostEqual(obj_at_theta.loc[0, "theta"], 0.0, places=8) + + def test_invalid_solver_name_raises_runtimeerror(self): + pest = _build_estimator([(1.0, 2.0), (2.0, 4.0)]) + with self.assertRaisesRegex( + RuntimeError, "Unknown solver in Q_Opt=not_a_solver" + ): + pest.theta_est(solver="not_a_solver") + + def test_theta_values_duplicate_columns_rejected(self): + pest = _build_estimator([(1.0, 2.0), (2.0, 4.0)]) + duplicate_cols = pd.DataFrame([[1.0, 2.0]], columns=["theta", "theta"]) + with self.assertRaisesRegex( + ValueError, "Duplicate theta names are not allowed" + ): + pest.objective_at_theta(theta_values=duplicate_cols) + + def test_count_total_experiments_multi_output(self): + exp_list = [ + LinearThetaExperiment(1.0, 2.0, include_second_output=True), + LinearThetaExperiment(2.0, 4.0, include_second_output=True), + ] + total_points = parmest._count_total_experiments(exp_list) + # The current parmest convention counts datapoints for one output family. + self.assertEqual(total_points, 2) + + ########################### # tests for deprecated UI # ########################### diff --git a/pyomo/contrib/parmest/tests/test_parmest_multistart.py b/pyomo/contrib/parmest/tests/test_parmest_multistart.py new file mode 100644 index 00000000000..e8770e07221 --- /dev/null +++ b/pyomo/contrib/parmest/tests/test_parmest_multistart.py @@ -0,0 +1,433 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of +# Sandia, LLC Under the terms of Contract DE-NA0003525 with National +# Technology and Engineering Solutions of Sandia, LLC, the U.S. Government +# retains certain rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import math + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo +from pyomo.common.dependencies import numpy as np, pandas as pd +from unittest.mock import patch + +import pyomo.contrib.parmest.parmest as parmest +from pyomo.contrib.parmest.experiment import Experiment + +ipopt_available = pyo.SolverFactory("ipopt").available() + + +class LinearThetaExperiment(Experiment): + def __init__(self, x, y): + self.x_data = x + self.y_data = y + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.theta = pyo.Var(initialize=0.0, bounds=(-10.0, 10.0)) + m.x = pyo.Param(initialize=float(self.x_data), mutable=False) + m.y = pyo.Var(initialize=float(self.y_data)) + m.eq = pyo.Constraint(expr=m.y == m.theta + m.x) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, float(self.y_data))]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + +class IndexedThetaExperiment(Experiment): + def __init__(self): + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.I = pyo.Set(initialize=["a", "b"]) + m.theta = pyo.Var(m.I, initialize={"a": 1.0, "b": 2.0}) + m.theta["a"].setlb(0.0) + m.theta["a"].setub(5.0) + m.theta["b"].setlb(0.0) + m.theta["b"].setub(5.0) + m.theta["a"].fix() + m.theta["b"].fix() + m.y = pyo.Var(initialize=3.0) + m.eq = pyo.Constraint(expr=m.y == m.theta["a"] + m.theta["b"]) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, 3.0)]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + +class NoBoundsExperiment(Experiment): + def __init__(self): + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.theta = pyo.Var(initialize=1.0) + m.y = pyo.Var(initialize=2.0) + m.eq = pyo.Constraint(expr=m.y == m.theta + 1.0) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, 2.0)]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + +class StartCoupledExperiment(Experiment): + """ + Model intentionally couples a fixed term ("bias") to theta_initial at + build time. This exposes stale-model bugs in multistart paths. + """ + + def __init__(self, theta_initial=None): + self.theta_initial = ( + theta_initial if theta_initial is not None else {"theta": 0.0} + ) + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.theta = pyo.Var( + initialize=float(self.theta_initial["theta"]), bounds=(-10.0, 10.0) + ) + m.bias = pyo.Param(initialize=float(self.theta_initial["theta"]), mutable=False) + m.y = pyo.Var(initialize=0.0) + m.eq = pyo.Constraint(expr=m.y == m.theta + m.bias) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, 0.0)]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + + def get_labeled_model(self): + if self.model is None: + self.create_model() + self.label_model() + return self.model + + +def _build_linear_estimator(): + exp_list = [LinearThetaExperiment(1.0, 2.0), LinearThetaExperiment(2.0, 3.0)] + return parmest.Estimator(exp_list, obj_function="SSE") + + +@unittest.skipIf( + not parmest.parmest_available, + "Cannot test parmest: required dependencies are missing", +) +class TestParmestMultistart(unittest.TestCase): + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_multistart_baseline_equivalence_n1(self): + pest = _build_linear_estimator() + obj1, theta1 = pest.theta_est() + _, best_theta, best_obj = pest.theta_est_multistart( + n_restarts=1, multistart_sampling_method="uniform_random", seed=7 + ) + self.assertAlmostEqual(obj1, best_obj, places=7) + self.assertAlmostEqual(theta1["theta"], best_theta["theta"], places=7) + + def test_uniform_sampling_is_deterministic_with_seed(self): + pest = _build_linear_estimator() + df1 = pest._generate_initial_theta( + seed=4, n_restarts=5, multistart_sampling_method="uniform_random" + ) + df2 = pest._generate_initial_theta( + seed=4, n_restarts=5, multistart_sampling_method="uniform_random" + ) + self.assertTrue(df1[["theta"]].equals(df2[["theta"]])) + + def test_uniform_sampling_changes_with_different_seed(self): + pest = _build_linear_estimator() + df1 = pest._generate_initial_theta( + seed=4, n_restarts=5, multistart_sampling_method="uniform_random" + ) + df2 = pest._generate_initial_theta( + seed=5, n_restarts=5, multistart_sampling_method="uniform_random" + ) + self.assertFalse(df1[["theta"]].equals(df2[["theta"]])) + + def test_latin_hypercube_sampling_is_deterministic(self): + pest = _build_linear_estimator() + df1 = pest._generate_initial_theta( + seed=11, n_restarts=4, multistart_sampling_method="latin_hypercube" + ) + df2 = pest._generate_initial_theta( + seed=11, n_restarts=4, multistart_sampling_method="latin_hypercube" + ) + self.assertTrue(df1[["theta"]].equals(df2[["theta"]])) + + def test_sobol_sampling_is_deterministic(self): + pest = _build_linear_estimator() + df1 = pest._generate_initial_theta( + seed=12, n_restarts=4, multistart_sampling_method="sobol_sampling" + ) + df2 = pest._generate_initial_theta( + seed=12, n_restarts=4, multistart_sampling_method="sobol_sampling" + ) + self.assertTrue(df1[["theta"]].equals(df2[["theta"]])) + + def test_generated_starts_are_within_bounds(self): + pest = _build_linear_estimator() + for method in ("uniform_random", "latin_hypercube", "sobol_sampling"): + df = pest._generate_initial_theta( + seed=1, n_restarts=8, multistart_sampling_method=method + ) + self.assertTrue(((df["theta"] >= -10.0) & (df["theta"] <= 10.0)).all()) + + def test_missing_bounds_raise_error(self): + pest = parmest.Estimator([NoBoundsExperiment()], obj_function="SSE") + with self.assertRaisesRegex( + ValueError, "lower and upper bounds for the theta values must be defined" + ): + pest._generate_initial_theta( + seed=1, n_restarts=2, multistart_sampling_method="uniform_random" + ) + + def test_invalid_bounds_raise_error(self): + class InvalidBoundsExperiment(Experiment): + def __init__(self): + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.theta = pyo.Var(initialize=1.0) + m.theta.setlb(2.0) + m.theta.setub(1.0) + m.y = pyo.Var(initialize=2.0) + m.eq = pyo.Constraint(expr=m.y == m.theta + 1.0) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, 2.0)]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + pest = parmest.Estimator([InvalidBoundsExperiment()], obj_function="SSE") + with self.assertRaisesRegex(ValueError, "lower bound must be less than"): + pest._generate_initial_theta( + seed=1, n_restarts=2, multistart_sampling_method="uniform_random" + ) + + def test_user_provided_values_dimension_mismatch_raises(self): + pest = _build_linear_estimator() + user_df = pd.DataFrame([[1.0, 2.0]], columns=["theta", "extra"]) + with self.assertRaisesRegex(ValueError, "exactly one column per theta name"): + pest.theta_est_multistart( + n_restarts=1, + multistart_sampling_method="user_provided_values", + user_provided_df=user_df, + ) + + def test_user_provided_values_column_order_maps_by_name(self): + pest = parmest.Estimator([IndexedThetaExperiment()], obj_function="SSE") + user_df = pd.DataFrame( + [[0.3, 4.2], [0.4, 4.1]], columns=["theta[b]", "theta[a]"] + ) + results_df, _, _ = pest.theta_est_multistart( + n_restarts=2, + multistart_sampling_method="user_provided_values", + user_provided_df=user_df, + ) + self.assertAlmostEqual(results_df.loc[0, "theta[a]"], 4.2, places=12) + self.assertAlmostEqual(results_df.loc[0, "theta[b]"], 0.3, places=12) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_state_isolation_between_starts(self): + pest = _build_linear_estimator() + init = pd.DataFrame([[-9.0], [9.0]], columns=["theta"]) + results_df, _, _ = pest.theta_est_multistart( + user_provided_df=init, save_results=False + ) + # Initial starts should remain exactly as supplied. + self.assertAlmostEqual(results_df.loc[0, "theta"], -9.0, places=12) + self.assertAlmostEqual(results_df.loc[1, "theta"], 9.0, places=12) + # Both runs converge to the same optimum, showing no cross-start contamination. + self.assertAlmostEqual( + results_df.loc[0, "converged_theta"], + results_df.loc[1, "converged_theta"], + places=8, + ) + + def test_one_start_failure_returns_best_feasible(self): + pest = _build_linear_estimator() + theta_values = pd.DataFrame([[-1.0], [2.0]], columns=["theta"]) + + def fake_q_opt(*args, **kwargs): + theta = kwargs["theta_vals"]["theta"] + if theta < 0: + raise RuntimeError("boom") + return 1.25, {"theta": 1.0}, pyo.TerminationCondition.optimal + + with patch.object(pest, "_Q_opt", side_effect=fake_q_opt): + results_df, best_theta, best_obj = pest.theta_est_multistart( + user_provided_df=theta_values, save_results=False + ) + + self.assertTrue( + str(results_df.loc[0, "solver termination"]).startswith("exception(start=0") + ) + self.assertAlmostEqual(best_obj, 1.25, places=12) + self.assertAlmostEqual(best_theta["theta"], 1.0, places=12) + + def test_all_starts_fail_returns_diagnostics(self): + pest = _build_linear_estimator() + theta_values = pd.DataFrame([[1.0], [2.0]], columns=["theta"]) + + def fake_q_opt(*args, **kwargs): + raise RuntimeError("all failed") + + with patch.object(pest, "_Q_opt", side_effect=fake_q_opt): + results_df, best_theta, best_obj = pest.theta_est_multistart( + user_provided_df=theta_values, save_results=False + ) + + self.assertIsNone(best_theta) + self.assertTrue(math.isnan(best_obj)) + self.assertTrue( + results_df["solver termination"] + .astype(str) + .str.contains("exception\\(start=", regex=True) + .all() + ) + + def test_best_selection_filters_nonoptimal_status(self): + pest = _build_linear_estimator() + theta_values = pd.DataFrame([[1.0], [2.0]], columns=["theta"]) + + def fake_q_opt(*args, **kwargs): + theta = kwargs["theta_vals"]["theta"] + if theta < 1.5: + return 0.1, {"theta": 0.1}, pyo.TerminationCondition.maxIterations + return 0.2, {"theta": 0.2}, pyo.TerminationCondition.optimal + + with patch.object(pest, "_Q_opt", side_effect=fake_q_opt): + _, best_theta, best_obj = pest.theta_est_multistart( + user_provided_df=theta_values, save_results=False + ) + + self.assertAlmostEqual(best_obj, 0.2, places=12) + self.assertAlmostEqual(best_theta["theta"], 0.2, places=12) + + def test_tie_breaking_is_deterministic_first_index(self): + pest = _build_linear_estimator() + theta_values = pd.DataFrame([[5.0], [6.0], [7.0]], columns=["theta"]) + + def fake_q_opt(*args, **kwargs): + theta = kwargs["theta_vals"]["theta"] + return 1.0, {"theta": theta}, pyo.TerminationCondition.optimal + + with patch.object(pest, "_Q_opt", side_effect=fake_q_opt): + _, best_theta, best_obj = pest.theta_est_multistart( + user_provided_df=theta_values, save_results=False + ) + + self.assertAlmostEqual(best_obj, 1.0, places=12) + self.assertAlmostEqual(best_theta["theta"], 5.0, places=12) + + def test_indexed_unknown_parameters_supported_in_sampling(self): + pest = parmest.Estimator([IndexedThetaExperiment()], obj_function="SSE") + df = pest._generate_initial_theta( + seed=10, n_restarts=3, multistart_sampling_method="uniform_random" + ) + self.assertTrue({"theta[a]", "theta[b]"}.issubset(set(df.columns))) + + def test_count_total_experiments_uses_one_output_family(self): + class MultiOutputExperiment(Experiment): + def create_model(self): + m = pyo.ConcreteModel() + m.theta = pyo.Var(initialize=0.0, bounds=(-10, 10)) + m.y = pyo.Var(initialize=1.0) + m.z = pyo.Var(initialize=2.0) + m.c1 = pyo.Constraint(expr=m.y == m.theta + 1.0) + m.c2 = pyo.Constraint(expr=m.z == 2.0 * m.theta + 2.0) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, 1.0), (m.z, 2.0)]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update([(m.theta, pyo.ComponentUID(m.theta))]) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None), (m.z, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + total_points = parmest._count_total_experiments( + [MultiOutputExperiment(), MultiOutputExperiment()] + ) + self.assertEqual(total_points, 2) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_multistart_results_reproducible_when_rerun_from_recorded_init(self): + pest = parmest.Estimator([StartCoupledExperiment()], obj_function="SSE") + init_df = pd.DataFrame([[-2.0], [1.5], [3.0]], columns=["theta"]) + results_df, _, _ = pest.theta_est_multistart( + user_provided_df=init_df, save_results=False + ) + + for _, row in results_df.iterrows(): + theta_init = {"theta": float(row["theta"])} + exp = StartCoupledExperiment(theta_initial=theta_init) + rerun = parmest.Estimator([exp], obj_function="SSE") + obj, theta = rerun.theta_est() + + self.assertTrue( + np.isclose(obj, row["final objective"], rtol=1e-6, atol=1e-8) + ) + self.assertTrue( + np.isclose(theta["theta"], row["converged_theta"], rtol=1e-6, atol=1e-8) + ) diff --git a/pyomo/contrib/parmest/tests/test_parmest_profile.py b/pyomo/contrib/parmest/tests/test_parmest_profile.py new file mode 100644 index 00000000000..6e63af13e73 --- /dev/null +++ b/pyomo/contrib/parmest/tests/test_parmest_profile.py @@ -0,0 +1,273 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of +# Sandia, LLC Under the terms of Contract DE-NA0003525 with National +# Technology and Engineering Solutions of Sandia, LLC, the U.S. Government +# retains certain rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo +from pyomo.common.dependencies import numpy as np, pandas as pd +from unittest.mock import patch + +import pyomo.contrib.parmest.parmest as parmest +from pyomo.contrib.parmest.experiment import Experiment + +ipopt_available = pyo.SolverFactory("ipopt").available() + + +class AffineTwoThetaExperiment(Experiment): + def __init__(self, x, y): + self.x_data = x + self.y_data = y + self.model = None + + def create_model(self): + m = pyo.ConcreteModel() + m.theta_a = pyo.Var(initialize=1.0, bounds=(0.0, 4.0)) + m.theta_b = pyo.Var(initialize=0.0, bounds=(-3.0, 3.0)) + m.x = pyo.Param(initialize=float(self.x_data), mutable=False) + m.y = pyo.Var(initialize=float(self.y_data)) + m.eq = pyo.Constraint(expr=m.y == m.theta_a * m.x + m.theta_b) + self.model = m + + def label_model(self): + m = self.model + m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.experiment_outputs.update([(m.y, float(self.y_data))]) + m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.unknown_parameters.update( + [ + (m.theta_a, pyo.ComponentUID(m.theta_a)), + (m.theta_b, pyo.ComponentUID(m.theta_b)), + ] + ) + m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL) + m.measurement_error.update([(m.y, None)]) + + def get_labeled_model(self): + self.create_model() + self.label_model() + return self.model + + +def _build_two_theta_estimator(): + exp_list = [ + AffineTwoThetaExperiment(1.0, 3.0), + AffineTwoThetaExperiment(2.0, 5.0), + AffineTwoThetaExperiment(3.0, 7.0), + ] + return parmest.Estimator(exp_list, obj_function="SSE") + + +@unittest.skipIf( + not parmest.parmest_available, + "Cannot test parmest: required dependencies are missing", +) +class TestParmestProfileLikelihood(unittest.TestCase): + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_profile_contains_baseline_minimum_point(self): + pest = _build_two_theta_estimator() + res = pest.profile_likelihood( + "theta_a", grid=[1.5, 2.0, 2.5], solver="ef_ipopt" + ) + prof = res["profiles"] + self.assertAlmostEqual(res["baseline"]["obj_hat"], prof["obj"].min(), places=8) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_profile_partial_fix_enforced(self): + pest = _build_two_theta_estimator() + res = pest.profile_likelihood( + "theta_a", grid=[1.5, 2.0, 2.5], solver="ef_ipopt" + ) + prof = res["profiles"] + self.assertTrue(np.allclose(prof["theta_value"], prof["theta__theta_a"])) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_profile_other_thetas_unfixed(self): + pest = _build_two_theta_estimator() + res = pest.profile_likelihood( + "theta_a", grid=[1.5, 2.0, 2.5], solver="ef_ipopt" + ) + prof = res["profiles"].sort_values("theta_value") + self.assertGreater( + prof["theta__theta_b"].max() - prof["theta__theta_b"].min(), 1e-6 + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_profile_repeatability_same_seed(self): + pest = _build_two_theta_estimator() + res1 = pest.profile_likelihood( + "theta_a", grid=[1.5, 2.0, 2.5], solver="ef_ipopt", seed=9 + ) + res2 = pest.profile_likelihood( + "theta_a", grid=[1.5, 2.0, 2.5], solver="ef_ipopt", seed=9 + ) + cols = [ + "theta_value", + "obj", + "status", + "success", + "theta__theta_a", + "theta__theta_b", + ] + self.assertTrue(res1["profiles"][cols].equals(res2["profiles"][cols])) + + def test_profile_warmstart_neighbor_values_used(self): + pest = _build_two_theta_estimator() + call_inits = [] + + def fake_q_opt(*args, **kwargs): + init = dict(kwargs["theta_vals"]) + fixed = kwargs["fixed_theta_values"] + call_inits.append(init) + return ( + float((fixed["theta_a"] - 2.0) ** 2), + {"theta_a": fixed["theta_a"], "theta_b": fixed["theta_a"] + 10.0}, + pyo.TerminationCondition.optimal, + ) + + with patch.object(pest, "_Q_opt", side_effect=fake_q_opt): + pest.profile_likelihood( + "theta_a", + grid=[2.0, 2.1, 2.2], + theta_hat={"theta_a": 2.0, "theta_b": 1.0}, + obj_hat=0.0, + warmstart="neighbor", + ) + + self.assertGreaterEqual(len(call_inits), 2) + self.assertAlmostEqual(call_inits[1]["theta_b"], 12.0, places=12) + + def test_profile_failure_recorded_continue(self): + pest = _build_two_theta_estimator() + + def fake_q_opt(*args, **kwargs): + a = kwargs["fixed_theta_values"]["theta_a"] + if abs(a - 2.0) < 1e-12: + raise RuntimeError("boom") + return ( + 1.0 + abs(a), + {"theta_a": a, "theta_b": 0.0}, + pyo.TerminationCondition.optimal, + ) + + with patch.object(pest, "_Q_opt", side_effect=fake_q_opt): + res = pest.profile_likelihood( + "theta_a", + grid=[1.9, 2.0, 2.1], + theta_hat={"theta_a": 2.0, "theta_b": 0.0}, + obj_hat=1.0, + ) + + prof = res["profiles"].sort_values("theta_value").reset_index(drop=True) + self.assertEqual(len(prof), 3) + self.assertIn("exception", str(prof.loc[1, "status"])) + self.assertFalse(bool(prof.loc[1, "success"])) + + def test_profile_all_failures_returns_structure(self): + pest = _build_two_theta_estimator() + + with patch.object(pest, "_Q_opt", side_effect=RuntimeError("all failed")): + res = pest.profile_likelihood( + "theta_a", + grid=[1.9, 2.0, 2.1], + theta_hat={"theta_a": 2.0, "theta_b": 0.0}, + obj_hat=1.0, + ) + + prof = res["profiles"] + self.assertEqual(len(prof), 3) + self.assertFalse(prof["success"].any()) + self.assertTrue(prof["status"].astype(str).str.contains("exception").all()) + + def test_profile_user_grid_preserved(self): + pest = _build_two_theta_estimator() + + with patch.object( + pest, + "_Q_opt", + side_effect=lambda *args, **kwargs: ( + 1.0, + {"theta_a": kwargs["fixed_theta_values"]["theta_a"], "theta_b": 0.0}, + pyo.TerminationCondition.optimal, + ), + ): + res = pest.profile_likelihood( + "theta_a", + grid=[1.2, 2.4, 2.0], + theta_hat={"theta_a": 2.0, "theta_b": 0.0}, + obj_hat=1.0, + ) + + attempted = sorted(res["profiles"]["theta_value"].tolist()) + self.assertEqual(attempted, [1.2, 2.0, 2.4]) + + def test_profile_auto_grid_includes_theta_hat(self): + pest = _build_two_theta_estimator() + grid = pest._build_profile_grid( + profiled_theta="theta_a", + grid=[1.0, 1.5, 2.5], + n_grid=5, + theta_hat={"theta_a": 2.0, "theta_b": 0.0}, + ) + self.assertIn(2.0, set(np.round(grid, 12))) + + def test_profile_result_columns_schema(self): + pest = _build_two_theta_estimator() + + with patch.object( + pest, + "_Q_opt", + side_effect=lambda *args, **kwargs: ( + 1.0, + {"theta_a": kwargs["fixed_theta_values"]["theta_a"], "theta_b": 0.0}, + pyo.TerminationCondition.optimal, + ), + ): + res = pest.profile_likelihood( + "theta_a", + grid=[1.9, 2.0], + theta_hat={"theta_a": 2.0, "theta_b": 0.0}, + obj_hat=1.0, + ) + prof = res["profiles"] + self.assertTrue( + set( + [ + "profiled_theta", + "theta_value", + "obj", + "delta_obj", + "lr_stat", + "status", + "success", + "solve_time", + ] + ).issubset(prof.columns) + ) + + @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") + def test_profile_baseline_from_multistart(self): + pest = _build_two_theta_estimator() + res = pest.profile_likelihood( + profiled_theta="theta_a", + n_grid=5, + use_multistart_for_baseline=True, + baseline_multistart_kwargs={ + "n_restarts": 3, + "multistart_sampling_method": "uniform_random", + "seed": 7, + }, + ) + self.assertIn("baseline", res) + self.assertIn("profiles", res) + self.assertTrue(np.isfinite(res["baseline"]["obj_hat"])) + self.assertGreaterEqual(res["profiles"].shape[0], 1) + + +if __name__ == "__main__": + unittest.main() diff --git a/pyomo/contrib/parmest/utils/create_ef.py b/pyomo/contrib/parmest/utils/create_ef.py index 56048fced10..17336114668 100644 --- a/pyomo/contrib/parmest/utils/create_ef.py +++ b/pyomo/contrib/parmest/utils/create_ef.py @@ -21,6 +21,7 @@ from pyomo.core import Objective +# File no longer used in parmest; retained for possible future use. def get_objs(scenario_instance): """return the list of objective functions for scenario_instance""" scenario_objs = scenario_instance.component_data_objects( diff --git a/pyomo/contrib/parmest/utils/mpi_utils.py b/pyomo/contrib/parmest/utils/mpi_utils.py index 82376062c45..30560a6a77d 100644 --- a/pyomo/contrib/parmest/utils/mpi_utils.py +++ b/pyomo/contrib/parmest/utils/mpi_utils.py @@ -10,6 +10,8 @@ from collections import OrderedDict import importlib +# ParallelTaskManager is used, MPI Interface is not. + """ This module is a collection of classes that provide a friendlier interface to MPI (through mpi4py). They help