Merge pull request #2 from kaklise/parmest-redesign

Parmest updates, minor code updates and docs for new API

Author: smartin71

Diff for: doc/OnlineDocs/contributed_packages/parmest/datarec.rst

@@ -3,35 +3,27 @@
 Data Reconciliation 
 ====================
 
-The method :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`
-can optionally return model values.  This feature can be used to return
-reconciled data using a user specified objective. In this case, the list
-of variable names the user wants to estimate (theta_names) is set to an
-empty list and the objective function is defined to minimize
+The optional argument ``return_values`` in :class:`~pyomo.contrib.parmest.parmest.Estimator.theta_est`
+can be used for data reconciliation or to return model values based on the specified objective. 
+
+For data reconciliation, the ``m.unknown_parameters`` is empty
+and the objective function is defined to minimize
 measurement to model error. Note that the model used for data
 reconciliation may differ from the model used for parameter estimation.
 
-The following example illustrates the use of parmest for data
-reconciliation.  The functions
+The functions
 :class:`~pyomo.contrib.parmest.graphics.grouped_boxplot` or
 :class:`~pyomo.contrib.parmest.graphics.grouped_violinplot` can be used
 to visually compare the original and reconciled data.
 
-Here's a stylized code snippet showing how box plots might be created:
-
-.. doctest::
-    :skipif: True
-
-    >>> import pyomo.contrib.parmest.parmest as parmest
-    >>> pest = parmest.Estimator(model_function, data, [], objective_function)
-    >>> obj, theta, data_rec = pest.theta_est(return_values=['A', 'B']) 
-    >>> parmest.graphics.grouped_boxplot(data, data_rec)
+The following example from the reactor design subdirectory returns reconciled values for experiment outputs 
+(`ca`, `cb`, `cc`, and `cd`) and then uses those values in 
+parameter estimation (`k1`, `k2`, and `k3`).
 
-Returned Values
-^^^^^^^^^^^^^^^
-
-Here's a full program that can be run to see returned values (in this case it
-is the response function that is defined in the model file):
+.. literalinclude:: ../../../../pyomo/contrib/parmest/examples/reactor_design/datarec_example.py
+   :language: python
+   
+The following example returns model values from a Pyomo Expression. 
 
 .. doctest::
    :skipif: not ipopt_available or not parmest_available
@@ -60,4 +52,3 @@ is the response function that is defined in the model file):
    >>> pest = parmest.Estimator(exp_list, obj_function=SSE, solver_options=None)
    >>> obj, theta, var_values = pest.theta_est(return_values=['response_function'])
    >>> #print(var_values)
-   

Diff for: doc/OnlineDocs/contributed_packages/parmest/driver.rst

@@ -4,7 +4,7 @@ Parameter Estimation
 ==================================
 
 Parameter Estimation using parmest requires a Pyomo model, experimental
-data which defines multiple scenarios, and a list of parameter names
+data which defines multiple scenarios, and parameters
 (thetas) to estimate.  parmest uses Pyomo [PyomoBookII]_ and (optionally) 
 mpi-sppy [mpisppy]_ to solve a
 two-stage stochastic programming problem, where the experimental data is
@@ -36,8 +36,8 @@ which includes the following methods:
    ~pyomo.contrib.parmest.parmest.Estimator.likelihood_ratio_test
    ~pyomo.contrib.parmest.parmest.Estimator.leaveNout_bootstrap_test
 
-Additional functions are available in parmest to group data, plot
-results, and fit distributions to theta values.
+Additional functions are available in parmest to plot
+results and fit distributions to theta values.
 
 .. autosummary::
    :nosignatures:
@@ -92,65 +92,43 @@ Optionally, solver options can be supplied, e.g.,
 
     >>> solver_options = {"max_iter": 6000}
     >>> pest = parmest.Estimator(exp_list, obj_function=SSE, solver_options=solver_options)
-        
-
-
-Model function
---------------
-
-The first argument is a function which uses data for a single scenario
-to return a populated and initialized Pyomo model for that scenario.
-
-Parameters that the user would like to estimate can be defined as
-**mutable parameters (Pyomo `Param`) or variables (Pyomo `Var`)**.  
-Within parmest, any parameters that are to be estimated are converted to unfixed variables. 
-Variables that are to be estimated are also unfixed.
-
-The model does not have to be specifically written as a 
-two-stage stochastic programming problem for parmest. 
-That is, parmest can modify the
-objective, see :ref:`ObjFunction` below.
-
-Data
-----
-
-The second argument is the data which will be used to populate the Pyomo
-model.  Supported data formats include:
-
-* **Pandas Dataframe** where each row is a separate scenario and column
-  names refer to observed quantities.  Pandas DataFrames are easily
-  stored and read in from csv, excel, or databases, or created directly
-  in Python.
-* **List of Pandas Dataframe** where each entry in the list is a separate scenario. 
-  Dataframes store observed quantities, referenced by index and column.
-* **List of dictionaries** where each entry in the list is a separate
-  scenario and the keys (or nested keys) refer to observed quantities.
-  Dictionaries are often preferred over DataFrames when using static and
-  time series data.  Dictionaries are easily stored and read in from
-  json or yaml files, or created directly in Python.
-* **List of json file names** where each entry in the list contains a
-  json file name for a separate scenario.  This format is recommended
-  when using large datasets in parallel computing.
-
-The data must be compatible with the model function that returns a
-populated and initialized Pyomo model for a single scenario.  Data can
-include multiple entries per variable (time series and/or duplicate
-sensors).  This information can be included in custom objective
-functions, see :ref:`ObjFunction` below.
-
-Theta names
------------
-
-The third argument is a list of parameters or variable names that the user wants to
-estimate.  The list contains strings with `Param` and/or `Var` names from the Pyomo
-model.
+
+
+List of experiment objects
+--------------------------
+
+The first argument is a list of experiment objects which is used to
+create one labeled model for each expeirment. 
+The template :class:`~pyomo.contrib.parmest.experiment.Experiment` 
+can be used to generate a list of experiment objects.
+
+A labeled Pyomo model ``m`` has the following additional suffixes (Pyomo `Suffix`):
+
+* ``m.experiment_outputs`` which defines experiment output (Pyomo `Param`, `Var`, or `Expression`)
+  and their associated data values (float, int).
+* ``m.unknown_parameters`` which defines the mutable parameters or variables (Pyomo `Param` or `Var`)
+  to estimate along with their component unique identifier (Pyomo `ComponentUID`). 
+  Within parmest, any parameters that are to be estimated are converted to unfixed variables. 
+  Variables that are to be estimated are also unfixed.
+
+The experiment class has one required method:
+
+* :class:`~pyomo.contrib.parmest.experiment.Experiment.get_labeled_model` which returns the labeled Pyomo model.
+  Note that the model does not have to be specifically written as a 
+  two-stage stochastic programming problem for parmest. 
+  That is, parmest can modify the
+  objective, see :ref:`ObjFunction` below.
+ 
+Parmest comes with several :ref:`examplesection` that illustrates how to set up the list of experiment objects.
+The examples commonly include additional :class:`~pyomo.contrib.parmest.experiment.Experiment` class methods to
+create the model, finalize the model, and label the model.  The user can customize methods to suit their needs.
 
 .. _ObjFunction:
 
 Objective function
 ------------------
 
-The fourth argument is an optional argument which defines the
+The second argument is an optional argument which defines the
 optimization objective function to use in parameter estimation.  
 
 If no objective function is specified, the Pyomo model is used "as is" and
@@ -161,20 +139,27 @@ 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 arguments
-include `model` and `data` and the objective function returns a Pyomo
+defined for parameter estimation.  The objective function 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.
+
 Suggested initialization procedure for parameter estimation problems
 --------------------------------------------------------------------
 
 To check the quality of initial guess values provided for the fitted parameters, we suggest solving a 
 square instance of the problem prior to solving the parameter estimation problem using the following steps:
 
-1. Create :class:`~pyomo.contrib.parmest.parmest.Estimator` object. To initialize the parameter estimation solve from the square problem solution, set optional argument ``solver_options = {bound_push: 1e-8}``.
+1. Create :class:`~pyomo.contrib.parmest.parmest.Estimator` object. To initialize the parameter 
+estimation solve from the square problem solution, set optional argument ``solver_options = {bound_push: 1e-8}``.
 
-2. Call :class:`~pyomo.contrib.parmest.parmest.Estimator.objective_at_theta` with optional argument ``(initialize_parmest_model=True)``. Different initial guess values for the fitted parameters can be provided using optional argument `theta_values` (**Pandas Dataframe**)
+2. Call :class:`~pyomo.contrib.parmest.parmest.Estimator.objective_at_theta` with optional 
+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`

Diff for: doc/OnlineDocs/contributed_packages/parmest/examples.rst

@@ -20,7 +20,7 @@ Additional use cases include:
 * Parameter estimation using mpi4py, the example saves results to a file
   for later analysis/graphics (semibatch example)
 
-The description below uses the reactor design example.  The file
+The example below uses the reactor design example.  The file
 **reactor_design.py** includes a function which returns an populated
 instance of the Pyomo model.  Note that the model is defined to maximize
 `cb` and that `k1`, `k2`, and `k3` are fixed.  The _main_ program is

Diff for: doc/OnlineDocs/contributed_packages/parmest/scencreate.rst

@@ -18,5 +18,5 @@ scenarios to the screen, accessing them via the ``ScensItator`` a ``print``
    :language: python
 
 .. note::
-   This example may produce an error message your version of Ipopt is not based
+   This example may produce an error message if your version of Ipopt is not based
    on a good linear solver.

Diff for: pyomo/contrib/parmest/examples/reactor_design/datarec_example.py

@@ -20,24 +20,18 @@
 np.random.seed(1234)
 
 
-def reactor_design_model_for_datarec():
-
-    # Unfix inlet concentration for data rec
-    model = reactor_design_model()
-    model.caf.fixed = False
-
-    return model
-
-
-class ReactorDesignExperimentPreDataRec(ReactorDesignExperiment):
+class ReactorDesignExperimentDataRec(ReactorDesignExperiment):
 
     def __init__(self, data, data_std, experiment_number):
 
         super().__init__(data, experiment_number)
         self.data_std = data_std
 
     def create_model(self):
-        self.model = m = reactor_design_model_for_datarec()
+
+        self.model = m = reactor_design_model()
+        m.caf.fixed = False
+
         return m
 
     def label_model(self):
@@ -121,23 +115,18 @@ def main():
     # Create an experiment list
     exp_list = []
     for i in range(data.shape[0]):
-        exp_list.append(ReactorDesignExperimentPreDataRec(data, data_std, i))
+        exp_list.append(ReactorDesignExperimentDataRec(data, data_std, i))
 
     # Define sum of squared error objective function for data rec
-    def SSE(model):
+    def SSE_with_std(model):
         expr = sum(
             ((y - yhat) / model.experiment_outputs_std[y]) ** 2
             for y, yhat in model.experiment_outputs.items()
         )
         return expr
 
-    # View one model & SSE
-    # exp0_model = exp_list[0].get_labeled_model()
-    # print(exp0_model.pprint())
-    # print(SSE(exp0_model))
-
     ### Data reconciliation
-    pest = parmest.Estimator(exp_list, obj_function=SSE)
+    pest = parmest.Estimator(exp_list, obj_function=SSE_with_std)
 
     obj, theta, data_rec = pest.theta_est(return_values=["ca", "cb", "cc", "cd", "caf"])
     print(obj)
@@ -157,7 +146,7 @@ def SSE(model):
     for i in range(data_rec.shape[0]):
         exp_list.append(ReactorDesignExperimentPostDataRec(data_rec, data_std, i))
 
-    pest = parmest.Estimator(exp_list, obj_function=SSE)
+    pest = parmest.Estimator(exp_list, obj_function=SSE_with_std)
     obj, theta = pest.theta_est()
     print(obj)
     print(theta)

Diff for: pyomo/contrib/parmest/experiment.py

@@ -1,14 +1,29 @@
-# The experiment class is a template for making experiment lists
-# to pass to parmest.  An experiment is a pyomo model "m" which has
-# additional suffixes:
-#   m.experiment_outputs -- which variables are experiment outputs
-#   m.unknown_parameters -- which variables are parameters to estimate
-# The experiment class has only one required method:
-#   get_labeled_model()
-# which returns the labeled pyomo model.
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2024
+#  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.
+#  ___________________________________________________________________________
 
 
 class Experiment:
+    """
+    The experiment class is a template for making experiment lists
+    to pass to parmest.
+
+    An experiment is a Pyomo model "m" which is labeled
+    with additional suffixes:
+    * m.experiment_outputs which defines experiment outputs
+    * m.unknown_parameters which defines parameters to estimate
+
+    The experiment class has one required method:
+    * get_labeled_model() which returns the labeled Pyomo model
+    """
+
     def __init__(self, model=None):
         self.model = model