diff --git a/pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py b/pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py new file mode 100644 index 00000000000..cd6b0b75748 --- /dev/null +++ b/pyomo/contrib/parmest/examples/reaction_kinetics/__init__.py @@ -0,0 +1,9 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright 2017 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. +# ___________________________________________________________________________ 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 new file mode 100644 index 00000000000..c17cd7fc22e --- /dev/null +++ b/pyomo/contrib/parmest/examples/reaction_kinetics/simple_reaction_parmest_example.py @@ -0,0 +1,108 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright 2017 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. +# ___________________________________________________________________________ +''' +Example from Y. Bard, "Nonlinear Parameter Estimation", (pg. 124) + +This example shows: +1. How to define the unknown (to be regressed parameters) with an index +2. How to call parmest to only estimate some of the parameters (and fix the rest) + +Code provided by Paul Akula. +''' +import pandas as pd +from pandas import DataFrame +from os import path + +from pyomo.environ import (ConcreteModel, Param, Var, PositiveReals, Objective, + Constraint, RangeSet, Expression, minimize, exp, value) + +#from idaes.core.util import get_default_solver +import pyomo.contrib.parmest.parmest as parmest + +# ======================================================================= +''' Data from Table 5.2 in Y. Bard, "Nonlinear Parameter Estimation", (pg. 124) +''' + +data = [{'experiment': 1, 'x1': 0.1, 'x2': 100, 'y': 0.98}, + {'experiment': 2, 'x1': 0.2, 'x2': 100, 'y': 0.983}, + {'experiment': 3, 'x1': 0.3, 'x2': 100, 'y': 0.955}, + {'experiment': 4, 'x1': 0.4, 'x2': 100, 'y': 0.979}, + {'experiment': 5, 'x1': 0.5, 'x2': 100, 'y': 0.993}, + {'experiment': 6, 'x1': 0.05, 'x2': 200, 'y': 0.626}, + {'experiment': 7, 'x1': 0.1, 'x2': 200, 'y': 0.544}, + {'experiment': 8, 'x1': 0.15, 'x2': 200, 'y': 0.455}, + {'experiment': 9, 'x1': 0.2, 'x2': 200, 'y': 0.225}, + {'experiment': 10, 'x1': 0.25, 'x2': 200, 'y': 0.167}, + {'experiment': 11, 'x1': 0.02, 'x2': 300, 'y': 0.566}, + {'experiment': 12, 'x1': 0.04, 'x2': 300, 'y': 0.317}, + {'experiment': 13, 'x1': 0.06, 'x2': 300, 'y': 0.034}, + {'experiment': 14, 'x1': 0.08, 'x2': 300, 'y': 0.016}, + {'experiment': 15, 'x1': 0.1, 'x2': 300, 'y': 0.006}] + +# ======================================================================= + +def simple_reaction_model(data): + + # Create the concrete model + model = ConcreteModel() + + model.x1 = Param(initialize=float(data['x1'])) + model.x2 = Param(initialize=float(data['x2'])) + + # Rate constants + model.rxn = RangeSet(2) + initial_guess = {1: 750, 2: 1200} + model.k = Var(model.rxn, initialize=initial_guess, within=PositiveReals) + + # reaction product + model.y = Expression(expr=exp(-model.k[1] * + model.x1 * exp(-model.k[2] / model.x2))) + + # 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 + +if __name__ == "__main__": + + # ======================================================================= + # Parameter estimation without covariance estimate + # Only estimate the parameter k[1]. The parameter k[2] will remain fixed + # at its initial value + theta_names = ['k[1]'] + pest = parmest.Estimator(simple_reaction_model, data, theta_names) + obj, theta = pest.theta_est() + print(obj) + print(theta) + print() + #======================================================================= + # Estimate both k1 and k2 and compute the covariance matrix + theta_names = ['k'] + pest = parmest.Estimator(simple_reaction_model, data, theta_names) + obj, theta, cov = pest.theta_est(calc_cov=True) + print(obj) + print(theta) + print(cov) diff --git a/pyomo/contrib/parmest/examples/reactor_design/__init__.py b/pyomo/contrib/parmest/examples/reactor_design/__init__.py index 8d1c8b69c3f..6b39dd18d6a 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/__init__.py +++ b/pyomo/contrib/parmest/examples/reactor_design/__init__.py @@ -1 +1,10 @@ - +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright 2017 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. +# ___________________________________________________________________________ + diff --git a/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py b/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py index d9f70706c29..cd6b0b75748 100644 --- a/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py +++ b/pyomo/contrib/parmest/examples/rooney_biegler/__init__.py @@ -1,4 +1,4 @@ - # ___________________________________________________________________________ +# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC diff --git a/pyomo/contrib/parmest/examples/semibatch/__init__.py b/pyomo/contrib/parmest/examples/semibatch/__init__.py index 5296dafcc78..6b39dd18d6a 100644 --- a/pyomo/contrib/parmest/examples/semibatch/__init__.py +++ b/pyomo/contrib/parmest/examples/semibatch/__init__.py @@ -1,4 +1,4 @@ - # ___________________________________________________________________________ +# ___________________________________________________________________________ # # Pyomo: Python Optimization Modeling Objects # Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC diff --git a/pyomo/contrib/parmest/parmest.py b/pyomo/contrib/parmest/parmest.py index 899d04eb2ab..127b7faa39d 100644 --- a/pyomo/contrib/parmest/parmest.py +++ b/pyomo/contrib/parmest/parmest.py @@ -273,7 +273,7 @@ class Estimator(object): model_function: function Function that generates an instance of the Pyomo model using 'data' as the input argument - data: pandas DataFrame, list of dictionaries, or list of json file names + data: pd.DataFrame, list of dictionaries, or list of json file names Data that is used to build an instance of the Pyomo model and build the objective function theta_names: list of strings @@ -337,7 +337,7 @@ def _create_parmest_model(self, data): # If the component that was found is not a variable, # this will generate an exception (and the warning # in the 'except') - var_validate.fixed = False + var_validate.unfix() # We want to standardize on the CUID string # representation self.theta_names[i] = repr(var_cuid) @@ -502,6 +502,8 @@ def _Q_opt(self, ThetaVals=None, solver="ef_ipopt", cov = 2 * sse / (n - l) * inv_red_hes cov = pd.DataFrame(cov, index=thetavals.keys(), columns=thetavals.keys()) + thetavals = pd.Series(thetavals) + if len(return_values) > 0: var_values = [] for exp_i in self.ef_instance.component_objects(Block, descend_into=False): @@ -646,7 +648,7 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov Parameters ---------- solver: string, optional - "ef_ipopt" or "k_aug". Default is "ef_ipopt". + Currently only "ef_ipopt" is supported. Default is "ef_ipopt". return_values: list, optional List of Variable names used to return values from the model bootlist: list, optional @@ -658,12 +660,10 @@ def theta_est(self, solver="ef_ipopt", return_values=[], bootlist=None, calc_cov ------- objectiveval: float The objective function value - thetavals: dict - A dictionary of all values for theta + thetavals: pd.Series + Estimated values for theta variable values: pd.DataFrame Variable values for each variable name in return_values (only for solver='ef_ipopt') - Hessian: dict - A dictionary of dictionaries for the Hessian (only for solver='k_aug') cov: pd.DataFrame Covariance matrix of the fitted parameters (only for solver='ef_ipopt') """ @@ -696,7 +696,7 @@ def theta_est_bootstrap(self, bootstrap_samples, samplesize=None, Returns ------- - bootstrap_theta: DataFrame + bootstrap_theta: pd.DataFrame Theta values for each sample and (if return_samples = True) the sample numbers used in each estimation """ @@ -758,7 +758,7 @@ def theta_est_leaveNout(self, lNo, lNo_samples=None, seed=None, Returns ------- - lNo_theta: DataFrame + lNo_theta: pd.DataFrame Theta values for each sample and (if return_samples = True) the sample numbers left out of each estimation """ @@ -887,12 +887,12 @@ def objective_at_theta(self, theta_values): Parameters ---------- - theta_values: DataFrame, columns=theta_names + theta_values: pd.DataFrame, columns=theta_names Values of theta used to compute the objective Returns ------- - obj_at_theta: DataFrame + obj_at_theta: pd.DataFrame Objective value for each theta (infeasible solutions are omitted). """ @@ -927,7 +927,7 @@ def likelihood_ratio_test(self, obj_at_theta, obj_value, alphas, Parameters ---------- - obj_at_theta: DataFrame, columns = theta_names + 'obj' + obj_at_theta: pd.DataFrame, columns = theta_names + 'obj' Objective values for each theta value (returned by objective_at_theta) obj_value: int or float @@ -939,10 +939,10 @@ def likelihood_ratio_test(self, obj_at_theta, obj_value, alphas, Returns ------- - LR: DataFrame + LR: pd.DataFrame Objective values for each theta value along with True or False for each alpha - thresholds: dictionary + thresholds: pd.Series If return_threshold = True, the thresholds are also returned. """ assert isinstance(obj_at_theta, pd.DataFrame) @@ -958,6 +958,8 @@ def likelihood_ratio_test(self, obj_at_theta, obj_value, alphas, thresholds[a] = obj_value * ((chi2_val / (S - 2)) + 1) LR[a] = LR['obj'] < thresholds[a] + thresholds = pd.Series(thresholds) + if return_thresholds: return LR, thresholds else: @@ -972,7 +974,7 @@ def confidence_region_test(self, theta_values, distribution, alphas, Parameters ---------- - theta_values: DataFrame, columns = theta_names + theta_values: pd.DataFrame, columns = theta_names Theta values used to generate a confidence region (generally returned by theta_est_bootstrap) distribution: string @@ -982,25 +984,24 @@ def confidence_region_test(self, theta_values, distribution, alphas, alphas: list List of alpha values used to determine if theta values are inside or outside the region. - test_theta_values: dictionary or DataFrame, keys/columns = theta_names, optional + test_theta_values: pd.Series or pd.DataFrame, keys/columns = theta_names, optional Additional theta values that are compared to the confidence region to determine if they are inside or outside. Returns - ------- - training_results: DataFrame + training_results: pd.DataFrame Theta value used to generate the confidence region along with True (inside) or False (outside) for each alpha - test_results: DataFrame + test_results: pd.DataFrame If test_theta_values is not None, returns test theta value along with True (inside) or False (outside) for each alpha """ assert isinstance(theta_values, pd.DataFrame) assert distribution in ['Rect', 'MVN', 'KDE'] assert isinstance(alphas, list) - assert isinstance(test_theta_values, (type(None), dict, pd.DataFrame)) + assert isinstance(test_theta_values, (type(None), dict, pd.Series, pd.DataFrame)) - if isinstance(test_theta_values, dict): + if isinstance(test_theta_values, (dict, pd.Series)): test_theta_values = pd.Series(test_theta_values).to_frame().transpose() training_results = theta_values.copy() diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index bac451ee419..c424c9303d2 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -41,13 +41,14 @@ @unittest.skipIf(not parmest.parmest_available, "Cannot test parmest: required dependencies are missing") @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") -class parmest_object_Tester_RB(unittest.TestCase): +class TestRooneyBiegler(unittest.TestCase): def setUp(self): from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import rooney_biegler_model + # Note, the data used in this test has been corrected to use data.loc[5,'hour'] = 7 (instead of 6) data = pd.DataFrame(data=[[1,8.3],[2,10.3],[3,19.0], - [4,16.0],[5,15.6],[6,19.8]], columns=['hour', 'y']) + [4,16.0],[5,15.6],[7,19.8]], columns=['hour', 'y']) theta_names = ['asymptote', 'rate_constant'] @@ -65,9 +66,9 @@ def SSE(model, data): def test_theta_est(self): objval, thetavals = self.pest.theta_est() - self.assertAlmostEqual(objval, 4.4675, places=2) - self.assertAlmostEqual(thetavals['asymptote'], 19.2189, places=2) # 19.1426 from the paper - self.assertAlmostEqual(thetavals['rate_constant'], 0.5312, places=2) # 0.5311 from the paper + 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 @unittest.skipIf(not graphics.imports_available, "parmest.graphics imports are unavailable") @@ -109,8 +110,8 @@ def test_likelihood_ratio(self): LR = self.pest.likelihood_ratio_test(obj_at_theta, objval, [0.8, 0.9, 1.0]) self.assertTrue(set(LR.columns) >= set([0.8, 0.9, 1.0])) - self.assertTrue(LR[0.8].sum() == 7) - self.assertTrue(LR[0.9].sum() == 11) + self.assertTrue(LR[0.8].sum() == 6) + self.assertTrue(LR[0.9].sum() == 10) self.assertTrue(LR[1.0].sum() == 60) # all true graphics.pairwise_plot(LR, thetavals, 0.8) @@ -185,37 +186,6 @@ def test_theta_k_aug_for_Hessian(self): objval, thetavals, Hessian = self.pest.theta_est(solver="k_aug") self.assertAlmostEqual(objval, 4.4675, places=2) - -''' -The test cases above were developed with a transcription mistake in the dataset. -This test works with the correct dataset. -''' -@unittest.skipIf(not parmest.parmest_available, - "Cannot test parmest: required dependencies are missing") -@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") -class parmest_object_Tester_RB_match_paper(unittest.TestCase): - - def setUp(self): - from pyomo.contrib.parmest.examples.rooney_biegler.rooney_biegler import rooney_biegler_model - - 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']) - - theta_names = ['asymptote', 'rate_constant'] - - def SSE(model, data): - expr = sum((data.y[i] - model.response_function[data.hour[i]])**2 for i in data.index) - return expr - - self.pest = parmest.Estimator(rooney_biegler_model, data, theta_names, SSE) - - def test_theta_est(self): - objval, thetavals = self.pest.theta_est(calc_cov=False) - - 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 - @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available") @unittest.skipIf(not parmest.inverse_reduced_hessian_available, "Cannot test covariance matrix: required ASL dependency is missing") @@ -244,7 +214,7 @@ def test_theta_est_cov(self): @unittest.skipIf(not parmest.parmest_available, "Cannot test parmest: required dependencies are missing") @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") -class Test_parmest_indexed_variables(unittest.TestCase): +class TestIndexedVariables(unittest.TestCase): def make_model(self, theta_names): @@ -262,7 +232,9 @@ def rooney_biegler_model_alternate(data): model.var_names = pyo.Set(initialize=['asymptote','rate_constant']) model.theta = pyo.Var(model.var_names, initialize={'asymptote':15, 'rate_constant':0.5}) - + model.theta['asymptote'].fixed = True # parmest will unfix theta variables, even when they are indexed + model.theta['rate_constant'].fixed = True + def response_rule(m, h): expr = m.theta['asymptote'] * (1 - pyo.exp(-m.theta['rate_constant'] * h)) return expr @@ -279,7 +251,18 @@ def SSE(model, data): return expr return parmest.Estimator(rooney_biegler_model_alternate, data, theta_names, SSE) + + def test_theta_est(self): + theta_names = ["theta"] + + pest = self.make_model(theta_names) + objval, thetavals = pest.theta_est(calc_cov=False) + + self.assertAlmostEqual(objval, 4.3317112, places=2) + self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) + self.assertAlmostEqual(thetavals["theta[rate_constant]"], 0.5311, places=2) + def test_theta_est_quotedIndex(self): theta_names = ["theta['asymptote']", "theta['rate_constant']"] @@ -288,8 +271,8 @@ def test_theta_est_quotedIndex(self): objval, thetavals = pest.theta_est(calc_cov=False) self.assertAlmostEqual(objval, 4.3317112, places=2) - self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) # 19.1426 from the paper - self.assertAlmostEqual(thetavals['theta[rate_constant]'], 0.5311, places=2) # 0.5311 from the paper + self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) + self.assertAlmostEqual(thetavals["theta[rate_constant]"], 0.5311, places=2) def test_theta_est_impliedStrIndex(self): @@ -300,7 +283,7 @@ def test_theta_est_impliedStrIndex(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual(thetavals["theta[asymptote]"], 19.1426, places=2) # 19.1426 from the paper - self.assertAlmostEqual(thetavals['theta[rate_constant]'], 0.5311, places=2) # 0.5311 from the paper + self.assertAlmostEqual(thetavals["theta[rate_constant]"], 0.5311, places=2) # 0.5311 from the paper @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available") @@ -308,7 +291,7 @@ def test_theta_est_impliedStrIndex(self): not parmest.inverse_reduced_hessian_available, "Cannot test covariance matrix: required ASL dependency is missing") def test_theta_est_cov(self): - theta_names = ["theta[asymptote]", "theta[rate_constant]"] + theta_names = ["theta"] pest = self.make_model(theta_names) objval, thetavals, cov = pest.theta_est(calc_cov=True) @@ -329,7 +312,7 @@ def test_theta_est_cov(self): "Cannot test parmest: required dependencies are missing") @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available") -class parmest_object_Tester_reactor_design(unittest.TestCase): +class TestReactorDesign(unittest.TestCase): def setUp(self): from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model @@ -389,7 +372,7 @@ def test_return_values(self): @unittest.skipIf(not graphics.imports_available, "parmest.graphics imports are unavailable") @unittest.skipIf(is_osx, "Disabling graphics tests on OSX due to issue in Matplotlib, see Pyomo PR #1337") -class parmest_graphics(unittest.TestCase): +class TestGraphics(unittest.TestCase): def setUp(self): self.A = pd.DataFrame(np.random.randint(0,100,size=(100,4)), columns=list('ABCD')) diff --git a/pyomo/contrib/parmest/tests/test_scenariocreator.py b/pyomo/contrib/parmest/tests/test_scenariocreator.py index ee58ebb8f4b..92794d334f4 100644 --- a/pyomo/contrib/parmest/tests/test_scenariocreator.py +++ b/pyomo/contrib/parmest/tests/test_scenariocreator.py @@ -31,7 +31,7 @@ @unittest.skipIf(not parmest.parmest_available, "Cannot test parmest: required dependencies are missing") @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") -class pamest_Scenario_creator_reactor_design(unittest.TestCase): +class TestScenarioReactorDesign(unittest.TestCase): def setUp(self): from pyomo.contrib.parmest.examples.reactor_design.reactor_design import reactor_design_model @@ -99,7 +99,7 @@ def test_no_csv_if_empty(self): @unittest.skipIf(not parmest.parmest_available, "Cannot test parmest: required dependencies are missing") @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available") -class pamest_Scenario_creator_semibatch(unittest.TestCase): +class TestScenarioSemibatch(unittest.TestCase): def setUp(self): import pyomo.contrib.parmest.examples.semibatch.semibatch as sb