Merge branch 'main' of github.com:Pyomo/pyomo into finalize-release

Author: mrmundt

doc/OnlineDocs/_templates/recursive-module.rst

@@ -94,7 +94,7 @@ Library Reference
    :toctree:
    :template: recursive-module.rst
    :recursive:
-{% for item in modules %}
+{% for item in all_modules %}
 {# Need item != tests for Sphinx >= 8.0; !endswith(.tests) for < 8.0 #}
 {% if item != 'tests' and not item.endswith('.tests')
    and item != 'examples' and not item.endswith('.examples') %}

doc/OnlineDocs/explanation/analysis/parmest/covariance.rst

@@ -1,16 +1,121 @@
 Covariance Matrix Estimation
-=================================
+============================
 
-If the optional argument ``calc_cov=True`` is specified for :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`, 
-parmest will calculate the covariance matrix :math:`V_{\theta}` as follows:
+The uncertainty in the estimated parameters is quantified using the covariance matrix.
+The diagonal of the covariance matrix contains the variance of the estimated parameters.
+Assuming Gaussian independent and identically distributed measurement errors, the
+covariance matrix of the estimated parameters can be computed using the following
+methods which have been implemented in parmest.
 
-.. math::
-   V_{\theta} = 2 \sigma^2 H^{-1} 
+1. Reduced Hessian Method
 
-This formula assumes all measurement errors are independent and identically distributed with 
-variance :math:`\sigma^2`. :math:`H^{-1}` is the inverse of the Hessian matrix for an unweighted 
-sum of least squares problem. Currently, the covariance approximation is only valid if the 
-objective given to parmest is the sum of squared error. Moreover, parmest approximates the 
-variance of the measurement errors as :math:`\sigma^2 = \frac{1}{n-l} \sum e_i^2` where :math:`n` is 
-the number of data points, :math:`l` is the number of fitted parameters, and :math:`e_i` is the 
-residual for experiment :math:`i`.
+    When the objective function is the sum of squared errors (SSE):
+    :math:`\text{SSE} = \sum_{i = 1}^n \left(y_{i} - \hat{y}_{i}\right)^2`,
+    the covariance matrix is:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = 2 \sigma^2 \left(\frac{\partial^2 \text{SSE}}
+        {\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}\right)^{-1}_{\boldsymbol{\theta}
+        = \boldsymbol{\theta}^*}
+
+    When the objective function is the weighted SSE (WSSE):
+    :math:`\text{WSSE} = \frac{1}{2} \left(\mathbf{y} - f(\mathbf{x};\boldsymbol{\theta})\right)^\text{T}
+    \mathbf{W} \left(\mathbf{y} - f(\mathbf{x};\boldsymbol{\theta})\right)`,
+    the covariance matrix is:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left(\frac{\partial^2 \text{WSSE}}
+        {\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}\right)^{-1}_{\boldsymbol{\theta}
+        = \boldsymbol{\theta}^*}
+
+    Where :math:`V_{\boldsymbol{\theta}}` is the covariance matrix of the estimated
+    parameters, :math:`y` are the observed measured variables, :math:`\hat{y}` are the
+    predicted measured variables, :math:`n` is the number of data points,
+    :math:`\boldsymbol{\theta}` are the unknown parameters, :math:`\boldsymbol{\theta^*}`
+    are the estimates of the unknown parameters, :math:`\mathbf{x}` are the decision
+    variables, and :math:`\mathbf{W}` is a diagonal matrix containing the inverse of the
+    variance of the measurement error, :math:`\sigma^2`. When the standard
+    deviation of the measurement error is not supplied by the user, parmest
+    approximates the variance of the measurement error as
+    :math:`\sigma^2 = \frac{1}{n-l} \sum e_i^2` where :math:`l` is the number of
+    fitted parameters, and :math:`e_i` is the residual for experiment :math:`i`.
+
+    In parmest, this method computes the inverse of the Hessian by scaling the
+    objective function (SSE or WSSE) with a constant probability factor.
+
+2. Finite Difference Method
+
+    In this method, the covariance matrix, :math:`V_{\boldsymbol{\theta}}`, is
+    calculated by applying the Gauss-Newton approximation to the Hessian,
+    :math:`\frac{\partial^2 \text{SSE}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}`
+    or
+    :math:`\frac{\partial^2 \text{WSSE}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}}`,
+    leading to:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left(\sum_{i = 1}^n \mathbf{G}_{i}^{\mathrm{T}} \mathbf{W}
+        \mathbf{G}_{i} \right)^{-1}
+
+    This method uses central finite difference to compute the Jacobian matrix,
+    :math:`\mathbf{G}_{i}`, for experiment :math:`i`, which is the sensitivity of
+    the measured variables with respect to the parameters, :math:`\boldsymbol{\theta}`.
+    :math:`\mathbf{W}` is a diagonal matrix containing the inverse of the variance
+    of the measurement errors, :math:`\sigma^2`.
+
+3. Automatic Differentiation Method
+
+    Similar to the finite difference method, the covariance matrix is calculated as:
+
+    .. math::
+       V_{\boldsymbol{\theta}} = \left( \sum_{i = 1}^n \mathbf{G}_{\text{kaug},\, i}^{\mathrm{T}}
+        \mathbf{W} \mathbf{G}_{\text{kaug},\, i} \right)^{-1}
+
+    However, this method uses the model optimality (KKT) condition to compute the
+    Jacobian matrix, :math:`\mathbf{G}_{\text{kaug},\, i}`, for experiment :math:`i`.
+
+The covariance matrix calculation is only supported with the built-in objective
+functions "SSE" or "SSE_weighted".
+
+In parmest, the covariance matrix can be calculated after defining the
+:class:`~pyomo.contrib.parmest.parmest.Estimator` object and estimating the unknown
+parameters using :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`. To
+estimate the covariance matrix, with the default method being "finite_difference", call
+the :class:`~pyomo.contrib.parmest.parmest.Estimator.cov_est` function, e.g.,
+
+.. testsetup:: *
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    # Data
+    import pandas as pd
+    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 the Experiment class
+    from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import RooneyBieglerExperiment
+
+    exp_list = []
+    for i in range(data.shape[0]):
+        exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))
+
+.. doctest::
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    >>> import pyomo.contrib.parmest.parmest as parmest
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
+    >>> obj_val, theta_val = pest.theta_est()
+    >>> cov = pest.cov_est()
+
+Optionally, one of the three methods; "reduced_hessian", "finite_difference",
+and "automatic_differentiation_kaug" can be supplied for the covariance calculation,
+e.g.,
+
+.. doctest::
+    :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
+
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
+    >>> obj_val, theta_val = pest.theta_est()
+    >>> cov_method = "reduced_hessian"
+    >>> cov = pest.cov_est(method=cov_method)

doc/OnlineDocs/explanation/analysis/parmest/datarec.rst

@@ -44,8 +44,8 @@ The following example returns model values from a Pyomo Expression.
 
    >>> # Define objective
    >>> def SSE(model):
-   ...     expr = (model.experiment_outputs[model.y]
-   ...             - model.response_function[model.experiment_outputs[model.hour]]
+   ...     expr = (model.experiment_outputs[model.y[model.hour]]
+   ...             - model.y[model.hour]
    ...            ) ** 2
    ...     return expr
 

doc/OnlineDocs/explanation/analysis/parmest/driver.rst

@@ -16,19 +16,21 @@ the model and the observations (typically defined as the sum of squared
 deviation between model values and observed values).
 
 If the Pyomo model is not formatted as a two-stage stochastic
-programming problem in this format, the user can supply a custom
-function to use as the second stage cost and the Pyomo model will be
+programming problem in this format, the user can choose either the
+built-in "SSE" or "SSE_weighted" objective functions, or supply a custom
+objective function to use as the second stage cost. The Pyomo model will then be
 modified within parmest to match the required specifications.
-The stochastic programming callback function is also defined within parmest.  The callback
-function returns a populated and initialized model for each scenario.
+The stochastic programming callback function is also defined within parmest.
+The callback function returns a populated and initialized model for each scenario.
 
-To use parmest, the user creates a :class:`~pyomo.contrib.parmest.parmest.Estimator` object 
-which includes the following methods:
+To use parmest, the user creates a :class:`~pyomo.contrib.parmest.parmest.Estimator`
+object which includes the following methods:
 
 .. autosummary::
    :nosignatures:
 
    ~pyomo.contrib.parmest.parmest.Estimator.theta_est
+   ~pyomo.contrib.parmest.parmest.Estimator.cov_est
    ~pyomo.contrib.parmest.parmest.Estimator.theta_est_bootstrap
    ~pyomo.contrib.parmest.parmest.Estimator.theta_est_leaveNout
    ~pyomo.contrib.parmest.parmest.Estimator.objective_at_theta
@@ -65,16 +67,9 @@ Section.
         columns=['hour', 'y'],
     )
 
-    # 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
-
     # Create an experiment list
     from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import RooneyBieglerExperiment
+
     exp_list = []
     for i in range(data.shape[0]):
         exp_list.append(RooneyBieglerExperiment(data.loc[i, :]))
@@ -83,15 +78,15 @@ Section.
     :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
 
     >>> import pyomo.contrib.parmest.parmest as parmest
-    >>> pest = parmest.Estimator(exp_list, obj_function=SSE)
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
 
 Optionally, solver options can be supplied, e.g.,
 
 .. doctest::
     :skipif: not __import__('pyomo.contrib.parmest.parmest').contrib.parmest.parmest.parmest_available
 
     >>> solver_options = {"max_iter": 6000}
-    >>> pest = parmest.Estimator(exp_list, obj_function=SSE, solver_options=solver_options)
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE", solver_options=solver_options)
 
 
 List of experiment objects
@@ -137,17 +132,20 @@ expressions that are used to build an objective for the two-stage
 stochastic programming problem.  
 
 If the Pyomo model is not written as a two-stage stochastic programming problem in
-this format, and/or if the user wants to use an objective that is
-different than the original model, a custom objective function can be
-defined for parameter estimation.  The objective function has a single argument, 
-which is the model from a single experiment.
+this format, the user can select the "SSE" or "SSE_weighted" built-in objective
+functions. If the user wants to use an objective that is different from the built-in
+options, a custom objective function can be defined for parameter estimation. However,
+covariance matrix estimation will not support this custom objective function. The objective
+function (built-in or custom) has a single argument, which is the model from a single
+experiment.
 The objective function returns a Pyomo
 expression which is used to define "SecondStageCost".  The objective
 function can be used to customize data points and weights that are used
 in parameter estimation.
 
-Parmest includes one built in objective function to compute the sum of squared errors ("SSE") between the 
-``m.experiment_outputs`` model values and data values.
+Parmest includes two built-in objective functions ("SSE" and "SSE_weighted") to compute
+the sum of squared errors between the ``m.experiment_outputs`` model values and
+data values.
 
 Suggested initialization procedure for parameter estimation problems
 --------------------------------------------------------------------
@@ -162,4 +160,11 @@ estimation solve from the square problem solution, set optional argument ``solve
 argument ``(initialize_parmest_model=True)``. Different initial guess values for the fitted 
 parameters can be provided using optional argument `theta_values` (**Pandas Dataframe**)
 
-3. Solve parameter estimation problem by calling :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`
+3. Solve parameter estimation problem by calling
+:class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`, e.g.,
+
+.. doctest::
+    :skipif: not parmest_available
+
+    >>> pest = parmest.Estimator(exp_list, obj_function="SSE")
+    >>> obj_val, theta_val = pest.theta_est()

doc/OnlineDocs/explanation/analysis/parmest/overview.rst

@@ -41,23 +41,32 @@ that for most experiments, only small parts of :math:`x` will change
 from one experiment to the next.
 
 The following least squares objective can be used to estimate parameter
-values, where data points are indexed by :math:`s=1,\ldots,S`
+values assuming Gaussian independent and identically distributed measurement
+errors, where data points are indexed by :math:`s=1,\ldots,S`
 
 .. math::
 
    \min_{{\theta}} Q({\theta};{\tilde{x}}, {\tilde{y}}) \equiv \sum_{s=1}^{S}q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) \;\;
 
-where
+where :math:`q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s})` can be:
 
-.. math::
+1. Sum of squared errors
+
+    .. math::
+
+       q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) =
+        \sum_{i=1}^{m}\left({\tilde{y}}_{s,i} - g_{i}({\tilde{x}}_{s};{\theta})\right)^{2}
+
+2. Weighted sum of squared errors
+
+    .. math::
 
-   q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) = \sum_{i=1}^{m}w_{i}\left[{\tilde{y}}_{si} - g_{i}({\tilde{x}}_{s};{\theta})\right]^{2}, 
+       q_{s}({\theta};{\tilde{x}}_{s}, {\tilde{y}}_{s}) =
+        \sum_{i=1}^{m}\left(\frac{{\tilde{y}}_{s,i} - g_{i}({\tilde{x}}_{s};{\theta})}{w_i}\right)^{2}
 
 i.e., the contribution of sample :math:`s` to :math:`Q`, where :math:`w
-\in \Re^{m}` is a vector of weights for the responses. For
-multi-dimensional :math:`y`, this is the squared weighted :math:`L_{2}`
-norm and for univariate :math:`y` the weighted squared deviation.
-Custom objectives can also be defined for parameter estimation.
+\in \Re^{m}` is a vector containing the standard deviation of the measurement
+errors of :math:`y`. Custom objectives can also be defined for parameter estimation.
 
 In the applications of interest to us, the function :math:`g(\cdot)` is
 usually defined as an optimization problem with a large number of

pyomo/_archive/__init__.py

@@ -13,6 +13,6 @@
 official Pyomo API.
 
 These modules are still importable through their old names via
-:func:`pyomo.common.moved_module()`
+:func:`pyomo.common.deprecation.moved_module()`
 
 """

pyomo/common/collections/__init__.py

@@ -10,10 +10,10 @@
 #  ___________________________________________________________________________
 
 
-from collections.abc import MutableMapping, MutableSet, Mapping, Set, Sequence
-from collections import UserDict
+from collections import OrderedDict, UserDict
+from collections.abc import Mapping, MutableMapping, MutableSet, Sequence, Set
 
-from .orderedset import OrderedDict, OrderedSet
+from .bunch import Bunch
 from .component_map import ComponentMap, DefaultComponentMap
 from .component_set import ComponentSet
-from .bunch import Bunch
+from .orderedset import OrderedSet

pyomo/common/collections/_hasher.py

@@ -0,0 +1,72 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2025
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+from collections import defaultdict
+
+
+class HashDispatcher(defaultdict):
+    """Dispatch table for generating "universal" hashing of all Python objects.
+
+    This class manages a dispatch table for providing hash functions for all Python
+    types.  When an object is passed to the Hasher, it determines the appropriate
+    hashing strategy based on the object's type:
+
+      - If a custom hashing function is registered for the type, it is used.
+      - If the object is natively hashable, the default hash is used.
+      - If the object is unhashable, the object's :func:`id()` is used as a fallback.
+
+    The Hasher also includes special handling for tuples by recursively applying the
+    appropriate hashing strategy to each element within the tuple.
+    """
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(lambda: self._missing_impl, *args, **kwargs)
+        self[tuple] = self._tuple
+
+    def _missing_impl(self, val):
+        try:
+            hash(val)
+            self[val.__class__] = self._hashable
+        except:
+            self[val.__class__] = self._unhashable
+        return self[val.__class__](val)
+
+    @staticmethod
+    def _hashable(val):
+        return val
+
+    @staticmethod
+    def _unhashable(val):
+        return id(val)
+
+    def _tuple(self, val):
+        return tuple(self[i.__class__](i) for i in val)
+
+    def hashable(self, obj, hashable=None):
+        if isinstance(obj, type):
+            cls = obj
+        else:
+            cls = type(obj)
+        if hashable is None:
+            fcn = self.get(cls, None)
+            if fcn is None:
+                raise KeyError(obj)
+            return fcn is self._hashable
+        self[cls] = self._hashable if hashable else self._unhashable
+
+
+#: The global 'hasher' instance for managing "universal" hashing.
+#:
+#: This instance of the :class:`HashDispatcher` is used by
+#: :class:`~pyomo.common.collections.component_map.ComponentMap` and
+#: :class:`~pyomo.common.collections.component_set.ComponentSet` for
+#: generating hashes for all Python and Pyomo types.
+hasher = HashDispatcher()