Merge pull request #3028 from smartin71/cleanup-parmest

Cleanup parmest

pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py

@@ -19,20 +19,20 @@
 
 def main():
     # Vars to estimate
-    theta_names = ['k1', 'k2', 'k3']
+    theta_names = ["k1", "k2", "k3"]
 
     # Data
     file_dirname = dirname(abspath(str(__file__)))
-    file_name = abspath(join(file_dirname, 'reactor_data.csv'))
+    file_name = abspath(join(file_dirname, "reactor_data.csv"))
     data = pd.read_csv(file_name)
 
     # Sum of squared error function
     def SSE(model, data):
         expr = (
-            (float(data['ca']) - model.ca) ** 2
-            + (float(data['cb']) - model.cb) ** 2
-            + (float(data['cc']) - model.cc) ** 2
-            + (float(data['cd']) - model.cd) ** 2
+            (float(data.iloc[0]["ca"]) - model.ca) ** 2
+            + (float(data.iloc[0]["cb"]) - model.cb) ** 2
+            + (float(data.iloc[0]["cc"]) - model.cc) ** 2
+            + (float(data.iloc[0]["cd"]) - model.cd) ** 2
         )
         return expr
 
@@ -46,13 +46,13 @@ def SSE(model, data):
     bootstrap_theta = pest.theta_est_bootstrap(50)
 
     # Plot results
-    parmest.graphics.pairwise_plot(bootstrap_theta, title='Bootstrap theta')
+    parmest.graphics.pairwise_plot(bootstrap_theta, title="Bootstrap theta")
     parmest.graphics.pairwise_plot(
         bootstrap_theta,
         theta,
         0.8,
-        ['MVN', 'KDE', 'Rect'],
-        title='Bootstrap theta with confidence regions',
+        ["MVN", "KDE", "Rect"],
+        title="Bootstrap theta with confidence regions",
     )
 
 

pyomo/contrib/parmest/examples/reactor_design/datarec_example.py

@@ -39,16 +39,16 @@ def generate_data():
     data = pd.DataFrame()
     ndata = 200
     # Normal distribution, mean = 3400, std = 500
-    data['ca'] = 500 * np.random.randn(ndata) + 3400
+    data["ca"] = 500 * np.random.randn(ndata) + 3400
     # Random distribution between 500 and 1500
-    data['cb'] = np.random.rand(ndata) * 1000 + 500
+    data["cb"] = np.random.rand(ndata) * 1000 + 500
     # Lognormal distribution
-    data['cc'] = np.random.lognormal(np.log(1600), 0.25, ndata)
+    data["cc"] = np.random.lognormal(np.log(1600), 0.25, ndata)
     # Triangular distribution between 1000 and 2000
-    data['cd'] = np.random.triangular(1000, 1800, 3000, size=ndata)
+    data["cd"] = np.random.triangular(1000, 1800, 3000, size=ndata)
 
-    data['sv'] = sv_real
-    data['caf'] = caf_real
+    data["sv"] = sv_real
+    data["caf"] = caf_real
 
     return data
 
@@ -61,10 +61,10 @@ def main():
     # Define sum of squared error objective function for data rec
     def SSE(model, data):
         expr = (
-            ((float(data['ca']) - model.ca) / float(data_std['ca'])) ** 2
-            + ((float(data['cb']) - model.cb) / float(data_std['cb'])) ** 2
-            + ((float(data['cc']) - model.cc) / float(data_std['cc'])) ** 2
-            + ((float(data['cd']) - model.cd) / float(data_std['cd'])) ** 2
+            ((float(data.iloc[0]["ca"]) - model.ca) / float(data_std["ca"])) ** 2
+            + ((float(data.iloc[0]["cb"]) - model.cb) / float(data_std["cb"])) ** 2
+            + ((float(data.iloc[0]["cc"]) - model.cc) / float(data_std["cc"])) ** 2
+            + ((float(data.iloc[0]["cd"]) - model.cd) / float(data_std["cd"])) ** 2
         )
         return expr
 
@@ -73,26 +73,26 @@ def SSE(model, data):
 
     pest = parmest.Estimator(reactor_design_model_for_datarec, data, theta_names, SSE)
 
-    obj, theta, data_rec = pest.theta_est(return_values=['ca', 'cb', 'cc', 'cd', 'caf'])
+    obj, theta, data_rec = pest.theta_est(return_values=["ca", "cb", "cc", "cd", "caf"])
     print(obj)
     print(theta)
 
     parmest.graphics.grouped_boxplot(
-        data[['ca', 'cb', 'cc', 'cd']],
-        data_rec[['ca', 'cb', 'cc', 'cd']],
-        group_names=['Data', 'Data Rec'],
+        data[["ca", "cb", "cc", "cd"]],
+        data_rec[["ca", "cb", "cc", "cd"]],
+        group_names=["Data", "Data Rec"],
     )
 
     ### Parameter estimation using reconciled data
-    theta_names = ['k1', 'k2', 'k3']
-    data_rec['sv'] = data['sv']
+    theta_names = ["k1", "k2", "k3"]
+    data_rec["sv"] = data["sv"]
 
     pest = parmest.Estimator(reactor_design_model, data_rec, theta_names, SSE)
     obj, theta = pest.theta_est()
     print(obj)
     print(theta)
 
-    theta_real = {'k1': 5.0 / 6.0, 'k2': 5.0 / 3.0, 'k3': 1.0 / 6000.0}
+    theta_real = {"k1": 5.0 / 6.0, "k2": 5.0 / 3.0, "k3": 1.0 / 6000.0}
     print(theta_real)
 
 

pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py

@@ -20,11 +20,11 @@
 
 def main():
     # Vars to estimate
-    theta_names = ['k1', 'k2', 'k3']
+    theta_names = ["k1", "k2", "k3"]
 
     # Data
     file_dirname = dirname(abspath(str(__file__)))
-    file_name = abspath(join(file_dirname, 'reactor_data.csv'))
+    file_name = abspath(join(file_dirname, "reactor_data.csv"))
     data = pd.read_csv(file_name)
 
     # Create more data for the example
@@ -37,10 +37,10 @@ def main():
     # Sum of squared error function
     def SSE(model, data):
         expr = (
-            (float(data['ca']) - model.ca) ** 2
-            + (float(data['cb']) - model.cb) ** 2
-            + (float(data['cc']) - model.cc) ** 2
-            + (float(data['cd']) - model.cd) ** 2
+            (float(data.iloc[0]["ca"]) - model.ca) ** 2
+            + (float(data.iloc[0]["cb"]) - model.cb) ** 2
+            + (float(data.iloc[0]["cc"]) - model.cc) ** 2
+            + (float(data.iloc[0]["cd"]) - model.cd) ** 2
         )
         return expr
 
@@ -68,7 +68,7 @@ def SSE(model, data):
     lNo = 25
     lNo_samples = 5
     bootstrap_samples = 20
-    dist = 'MVN'
+    dist = "MVN"
     alphas = [0.7, 0.8, 0.9]
 
     results = pest.leaveNout_bootstrap_test(
@@ -84,8 +84,8 @@ def SSE(model, data):
             bootstrap_results,
             theta_est_N,
             alpha,
-            ['MVN'],
-            title='Alpha: ' + str(alpha) + ', ' + str(theta_est_N.loc[0, alpha]),
+            ["MVN"],
+            title="Alpha: " + str(alpha) + ", " + str(theta_est_N.loc[0, alpha]),
         )
 
     # Extract the percent of points that are within the alpha region

pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py

@@ -21,20 +21,20 @@
 
 def main():
     # Vars to estimate
-    theta_names = ['k1', 'k2', 'k3']
+    theta_names = ["k1", "k2", "k3"]
 
     # Data
     file_dirname = dirname(abspath(str(__file__)))
-    file_name = abspath(join(file_dirname, 'reactor_data.csv'))
+    file_name = abspath(join(file_dirname, "reactor_data.csv"))
     data = pd.read_csv(file_name)
 
     # Sum of squared error function
     def SSE(model, data):
         expr = (
-            (float(data['ca']) - model.ca) ** 2
-            + (float(data['cb']) - model.cb) ** 2
-            + (float(data['cc']) - model.cc) ** 2
-            + (float(data['cd']) - model.cd) ** 2
+            (float(data.iloc[0]["ca"]) - model.ca) ** 2
+            + (float(data.iloc[0]["cb"]) - model.cb) ** 2
+            + (float(data.iloc[0]["cc"]) - model.cc) ** 2
+            + (float(data.iloc[0]["cd"]) - model.cd) ** 2
         )
         return expr
 
@@ -48,15 +48,15 @@ def SSE(model, data):
     k1 = [0.8, 0.85, 0.9]
     k2 = [1.6, 1.65, 1.7]
     k3 = [0.00016, 0.000165, 0.00017]
-    theta_vals = pd.DataFrame(list(product(k1, k2, k3)), columns=['k1', 'k2', 'k3'])
+    theta_vals = pd.DataFrame(list(product(k1, k2, k3)), columns=["k1", "k2", "k3"])
     obj_at_theta = pest.objective_at_theta(theta_vals)
 
     # Run the likelihood ratio test
     LR = pest.likelihood_ratio_test(obj_at_theta, obj, [0.8, 0.85, 0.9, 0.95])
 
     # Plot results
     parmest.graphics.pairwise_plot(
-        LR, theta, 0.9, title='LR results within 90% confidence region'
+        LR, theta, 0.9, title="LR results within 90% confidence region"
     )
 
 

pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py

@@ -21,23 +21,23 @@ def main():
     # Parameter estimation using multisensor data
 
     # Vars to estimate
-    theta_names = ['k1', 'k2', 'k3']
+    theta_names = ["k1", "k2", "k3"]
 
     # Data, includes multiple sensors for ca and cc
     file_dirname = dirname(abspath(str(__file__)))
-    file_name = abspath(join(file_dirname, 'reactor_data_multisensor.csv'))
+    file_name = abspath(join(file_dirname, "reactor_data_multisensor.csv"))
     data = pd.read_csv(file_name)
 
     # Sum of squared error function
     def SSE_multisensor(model, data):
         expr = (
-            ((float(data['ca1']) - model.ca) ** 2) * (1 / 3)
-            + ((float(data['ca2']) - model.ca) ** 2) * (1 / 3)
-            + ((float(data['ca3']) - model.ca) ** 2) * (1 / 3)
-            + (float(data['cb']) - model.cb) ** 2
-            + ((float(data['cc1']) - model.cc) ** 2) * (1 / 2)
-            + ((float(data['cc2']) - model.cc) ** 2) * (1 / 2)
-            + (float(data['cd']) - model.cd) ** 2
+            ((float(data.iloc[0]["ca1"]) - model.ca) ** 2) * (1 / 3)
+            + ((float(data.iloc[0]["ca2"]) - model.ca) ** 2) * (1 / 3)
+            + ((float(data.iloc[0]["ca3"]) - model.ca) ** 2) * (1 / 3)
+            + (float(data.iloc[0]["cb"]) - model.cb) ** 2
+            + ((float(data.iloc[0]["cc1"]) - model.cc) ** 2) * (1 / 2)
+            + ((float(data.iloc[0]["cc2"]) - model.cc) ** 2) * (1 / 2)
+            + (float(data.iloc[0]["cd"]) - model.cd) ** 2
         )
         return expr
 

pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py

@@ -19,20 +19,20 @@
 
 def main():
     # Vars to estimate
-    theta_names = ['k1', 'k2', 'k3']
+    theta_names = ["k1", "k2", "k3"]
 
     # Data
     file_dirname = dirname(abspath(str(__file__)))
-    file_name = abspath(join(file_dirname, 'reactor_data.csv'))
+    file_name = abspath(join(file_dirname, "reactor_data.csv"))
     data = pd.read_csv(file_name)
 
     # Sum of squared error function
     def SSE(model, data):
         expr = (
-            (float(data['ca']) - model.ca) ** 2
-            + (float(data['cb']) - model.cb) ** 2
-            + (float(data['cc']) - model.cc) ** 2
-            + (float(data['cd']) - model.cd) ** 2
+            (float(data.iloc[0]["ca"]) - model.ca) ** 2
+            + (float(data.iloc[0]["cb"]) - model.cb) ** 2
+            + (float(data.iloc[0]["cc"]) - model.cc) ** 2
+            + (float(data.iloc[0]["cd"]) - model.cd) ** 2
         )
         return expr
 
@@ -46,11 +46,11 @@ def SSE(model, data):
     k1_expected = 5.0 / 6.0
     k2_expected = 5.0 / 3.0
     k3_expected = 1.0 / 6000.0
-    relative_error = abs(theta['k1'] - k1_expected) / k1_expected
+    relative_error = abs(theta["k1"] - k1_expected) / k1_expected
     assert relative_error < 0.05
-    relative_error = abs(theta['k2'] - k2_expected) / k2_expected
+    relative_error = abs(theta["k2"] - k2_expected) / k2_expected
     assert relative_error < 0.05
-    relative_error = abs(theta['k3'] - k3_expected) / k3_expected
+    relative_error = abs(theta["k3"] - k3_expected) / k3_expected
     assert relative_error < 0.05
 
 

pyomo/contrib/parmest/examples/reactor_design/reactor_design.py

@@ -37,10 +37,20 @@ def reactor_design_model(data):
     )  # m^3/(gmol min)
 
     # Inlet concentration of A, gmol/m^3
-    model.caf = Param(initialize=float(data['caf']), within=PositiveReals)
+    if isinstance(data, dict) or isinstance(data, pd.Series):
+        model.caf = Param(initialize=float(data["caf"]), within=PositiveReals)
+    elif isinstance(data, pd.DataFrame):
+        model.caf = Param(initialize=float(data.iloc[0]["caf"]), within=PositiveReals)
+    else:
+        raise ValueError("Unrecognized data type.")
 
     # Space velocity (flowrate/volume)
-    model.sv = Param(initialize=float(data['sv']), within=PositiveReals)
+    if isinstance(data, dict) or isinstance(data, pd.Series):
+        model.sv = Param(initialize=float(data["sv"]), within=PositiveReals)
+    elif isinstance(data, pd.DataFrame):
+        model.sv = Param(initialize=float(data.iloc[0]["sv"]), within=PositiveReals)
+    else:
+        raise ValueError("Unrecognized data type.")
 
     # Outlet concentration of each component
     model.ca = Var(initialize=5000.0, within=PositiveReals)
@@ -81,12 +91,12 @@ def main():
     sv_values = [1.0 + v * 0.05 for v in range(1, 20)]
     caf = 10000
     for sv in sv_values:
-        model = reactor_design_model({'caf': caf, 'sv': sv})
-        solver = SolverFactory('ipopt')
+        model = reactor_design_model(pd.DataFrame(data={"caf": [caf], "sv": [sv]}))
+        solver = SolverFactory("ipopt")
         solver.solve(model)
         results.append([sv, caf, model.ca(), model.cb(), model.cc(), model.cd()])
 
-    results = pd.DataFrame(results, columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'])
+    results = pd.DataFrame(results, columns=["sv", "caf", "ca", "cb", "cc", "cd"])
     print(results)