diff --git a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py b/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py index cf1b8a2de23..e2d172f34f6 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py +++ b/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", ) diff --git a/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py b/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py index 811571e20ed..cfd3891c00e 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py +++ b/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) diff --git a/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py b/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py index 95af53e63d3..6952a7fc733 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py +++ b/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 diff --git a/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py b/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py index 13a40774740..a0fe6f22305 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py +++ b/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,7 +48,7 @@ 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 @@ -56,7 +56,7 @@ def SSE(model, data): # 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" ) diff --git a/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py b/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py index bc564cbdfd3..a92ac626fae 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py +++ b/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 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 334dfa264a4..581d3904c04 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py +++ b/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 diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py index f4cd6c8dbf5..16f65e236eb 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/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) diff --git a/pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv b/pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv new file mode 100644 index 00000000000..79f03e07dcd --- /dev/null +++ b/pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv @@ -0,0 +1,1009 @@ +,k1,k2,E1,E2,obj +0,4,40,29000,38000,667.4023645794207 +1,4,40,29000,38500,665.8312183437167 +2,4,40,29000,39000,672.7539769993407 +3,4,40,29000,39500,684.9503752463216 +4,4,40,29000,40000,699.985589093255 +5,4,40,29000,40500,716.1241770970677 +6,4,40,29000,41000,732.2023201586336 +7,4,40,29000,41500,747.4931745925483 +8,4,40,29500,38000,907.4405527163311 +9,4,40,29500,38500,904.2229271927299 +10,4,40,29500,39000,907.6942345285257 +11,4,40,29500,39500,915.4570013614677 +12,4,40,29500,40000,925.65401444575 +13,4,40,29500,40500,936.9348578520337 +14,4,40,29500,41000,948.3759339765711 +15,4,40,29500,41500,959.386491783636 +16,4,40,30000,38000,1169.8685711377334 +17,4,40,30000,38500,1166.2211505723928 +18,4,40,30000,39000,1167.702295374574 +19,4,40,30000,39500,1172.5517020611685 +20,4,40,30000,40000,1179.3820406408263 +21,4,40,30000,40500,1187.1698633839655 +22,4,40,30000,41000,1195.2047840919602 +23,4,40,30000,41500,1203.0241101248102 +24,4,40,30500,38000,1445.9591944684807 +25,4,40,30500,38500,1442.6632745483 +26,4,40,30500,39000,1443.1982444457385 +27,4,40,30500,39500,1446.2833842279929 +28,4,40,30500,40000,1450.9012120934779 +29,4,40,30500,40500,1456.295140290636 +30,4,40,30500,41000,1461.9350767569827 +31,4,40,30500,41500,1467.4715014446226 +32,4,40,31000,38000,1726.8744994061449 +33,4,40,31000,38500,1724.2679845375048 +34,4,40,31000,39000,1724.4550886870552 +35,4,40,31000,39500,1726.5124587129135 +36,4,40,31000,40000,1729.7061680616455 +37,4,40,31000,40500,1733.48893482641 +38,4,40,31000,41000,1737.4753558920438 +39,4,40,31000,41500,1741.4093763605517 +40,4,40,31500,38000,2004.1978135112938 +41,4,40,31500,38500,2002.2807839860222 +42,4,40,31500,39000,2002.3676405166086 +43,4,40,31500,39500,2003.797808439923 +44,4,40,31500,40000,2006.048051591001 +45,4,40,31500,40500,2008.7281679153625 +46,4,40,31500,41000,2011.5626384878237 +47,4,40,31500,41500,2014.3675286347284 +48,4,80,29000,38000,845.8197358579285 +49,4,80,29000,38500,763.5039795545781 +50,4,80,29000,39000,709.8529964173656 +51,4,80,29000,39500,679.4215539491266 +52,4,80,29000,40000,666.4876088521157 +53,4,80,29000,40500,665.978271760966 +54,4,80,29000,41000,673.7240200504901 +55,4,80,29000,41500,686.4763909417914 +56,4,80,29500,38000,1042.519415429413 +57,4,80,29500,38500,982.8097210678039 +58,4,80,29500,39000,942.2990207573541 +59,4,80,29500,39500,917.9550916645245 +60,4,80,29500,40000,906.3116029967189 +61,4,80,29500,40500,904.0326666308792 +62,4,80,29500,41000,908.1964630052729 +63,4,80,29500,41500,916.4222043837499 +64,4,80,30000,38000,1271.1030403496538 +65,4,80,30000,38500,1227.7527550544085 +66,4,80,30000,39000,1197.433957624904 +67,4,80,30000,39500,1178.447676126182 +68,4,80,30000,40000,1168.645219243497 +69,4,80,30000,40500,1165.7995210546096 +70,4,80,30000,41000,1167.8586496250396 +71,4,80,30000,41500,1173.0949214020527 +72,4,80,30500,38000,1520.8220402652044 +73,4,80,30500,38500,1489.2563260709424 +74,4,80,30500,39000,1466.8099189128857 +75,4,80,30500,39500,1452.4352624958806 +76,4,80,30500,40000,1444.7074679423818 +77,4,80,30500,40500,1442.0820578624343 +78,4,80,30500,41000,1443.099006489627 +79,4,80,30500,41500,1446.5106517200784 +80,4,80,31000,38000,1781.149136032395 +81,4,80,31000,38500,1758.2414369536502 +82,4,80,31000,39000,1741.891639711003 +83,4,80,31000,39500,1731.358661496594 +84,4,80,31000,40000,1725.6231647999593 +85,4,80,31000,40500,1723.5757174297378 +86,4,80,31000,41000,1724.1680229486278 +87,4,80,31000,41500,1726.5050840601884 +88,4,80,31500,38000,2042.8335948845602 +89,4,80,31500,38500,2026.3067503042414 +90,4,80,31500,39000,2014.5720701940838 +91,4,80,31500,39500,2007.0463766643977 +92,4,80,31500,40000,2002.9647983728314 +93,4,80,31500,40500,2001.5163951989875 +94,4,80,31500,41000,2001.9474217001339 +95,4,80,31500,41500,2003.6204088755821 +96,4,120,29000,38000,1176.0713512305115 +97,4,120,29000,38500,1016.8213383282462 +98,4,120,29000,39000,886.0136231565133 +99,4,120,29000,39500,789.0101180066036 +100,4,120,29000,40000,724.5420056133441 +101,4,120,29000,40500,686.6877602625062 +102,4,120,29000,41000,668.8129085873959 +103,4,120,29000,41500,665.1167761036883 +104,4,120,29500,38000,1263.887274509128 +105,4,120,29500,38500,1155.6528408872423 +106,4,120,29500,39000,1066.393539894248 +107,4,120,29500,39500,998.9931006471243 +108,4,120,29500,40000,952.36314487701 +109,4,120,29500,40500,923.4000293372077 +110,4,120,29500,41000,908.407361383214 +111,4,120,29500,41500,903.8136176328255 +112,4,120,30000,38000,1421.1418235449091 +113,4,120,30000,38500,1347.114022652679 +114,4,120,30000,39000,1285.686103704643 +115,4,120,30000,39500,1238.2456448658272 +116,4,120,30000,40000,1204.3526810790904 +117,4,120,30000,40500,1182.4272879027071 +118,4,120,30000,41000,1170.3447810121902 +119,4,120,30000,41500,1165.8422968073423 +120,4,120,30500,38000,1625.5588911535713 +121,4,120,30500,38500,1573.5546642859429 +122,4,120,30500,39000,1530.1592840718379 +123,4,120,30500,39500,1496.2087139473604 +124,4,120,30500,40000,1471.525855239756 +125,4,120,30500,40500,1455.2084749904016 +126,4,120,30500,41000,1445.9160840082027 +127,4,120,30500,41500,1442.1255377330835 +128,4,120,31000,38000,1855.8467211183756 +129,4,120,31000,38500,1818.4368412235558 +130,4,120,31000,39000,1787.25956706785 +131,4,120,31000,39500,1762.8169908546402 +132,4,120,31000,40000,1744.9825741661596 +133,4,120,31000,40500,1733.136625016882 +134,4,120,31000,41000,1726.3352245899828 +135,4,120,31000,41500,1723.492199933745 +136,4,120,31500,38000,2096.6479813687533 +137,4,120,31500,38500,2069.3606691038876 +138,4,120,31500,39000,2046.792043575205 +139,4,120,31500,39500,2029.2128703900223 +140,4,120,31500,40000,2016.4664599897606 +141,4,120,31500,40500,2008.054814885348 +142,4,120,31500,41000,2003.2622557140814 +143,4,120,31500,41500,2001.289784483679 +144,7,40,29000,38000,149.32898706737052 +145,7,40,29000,38500,161.04814413969586 +146,7,40,29000,39000,187.87801343005242 +147,7,40,29000,39500,223.00789161520424 +148,7,40,29000,40000,261.66779887964003 +149,7,40,29000,40500,300.676316191238 +150,7,40,29000,41000,338.04021206995765 +151,7,40,29000,41500,372.6191631389286 +152,7,40,29500,38000,276.6495061185777 +153,7,40,29500,38500,282.1304583501965 +154,7,40,29500,39000,300.91417483065254 +155,7,40,29500,39500,327.24304394350395 +156,7,40,29500,40000,357.0561976596432 +157,7,40,29500,40500,387.61662064170207 +158,7,40,29500,41000,417.1836349752378 +159,7,40,29500,41500,444.73705844573243 +160,7,40,30000,38000,448.0380830353589 +161,7,40,30000,38500,448.8094536459122 +162,7,40,30000,39000,460.77530593327293 +163,7,40,30000,39500,479.342874472736 +164,7,40,30000,40000,501.20694459059405 +165,7,40,30000,40500,524.0971649678811 +166,7,40,30000,41000,546.539334134893 +167,7,40,30000,41500,567.6447156158981 +168,7,40,30500,38000,657.9909416906933 +169,7,40,30500,38500,655.7465129488842 +170,7,40,30500,39000,662.5420970804985 +171,7,40,30500,39500,674.8914651553109 +172,7,40,30500,40000,690.2111920703564 +173,7,40,30500,40500,706.6833639709198 +174,7,40,30500,41000,723.0994507096715 +175,7,40,30500,41500,738.7096013891406 +176,7,40,31000,38000,899.1769906655776 +177,7,40,31000,38500,895.4391505892945 +178,7,40,31000,39000,898.7695629120826 +179,7,40,31000,39500,906.603316771593 +180,7,40,31000,40000,916.9811481373996 +181,7,40,31000,40500,928.4913367709245 +182,7,40,31000,41000,940.1744934710283 +183,7,40,31000,41500,951.4199286075984 +184,7,40,31500,38000,1163.093373675207 +185,7,40,31500,38500,1159.0457727559028 +186,7,40,31500,39000,1160.3831770028223 +187,7,40,31500,39500,1165.2451698296604 +188,7,40,31500,40000,1172.1768190340001 +189,7,40,31500,40500,1180.1105659428963 +190,7,40,31500,41000,1188.3083929833688 +191,7,40,31500,41500,1196.29112579565 +192,7,80,29000,38000,514.0332369183081 +193,7,80,29000,38500,329.3645784712966 +194,7,80,29000,39000,215.73000998706416 +195,7,80,29000,39500,162.37338399591852 +196,7,80,29000,40000,149.8401793263549 +197,7,80,29000,40500,162.96125998112578 +198,7,80,29000,41000,191.173279165834 +199,7,80,29000,41500,227.2781971491003 +200,7,80,29500,38000,623.559246695578 +201,7,80,29500,38500,448.60620511421484 +202,7,80,29500,39000,344.21940687907573 +203,7,80,29500,39500,292.9758707105001 +204,7,80,29500,40000,277.07670134364804 +205,7,80,29500,40500,283.5158840045542 +206,7,80,29500,41000,303.33951582820265 +207,7,80,29500,41500,330.43357046741954 +208,7,80,30000,38000,732.5907387079073 +209,7,80,30000,38500,593.1926567994672 +210,7,80,30000,39000,508.5638538704666 +211,7,80,30000,39500,464.47881763522037 +212,7,80,30000,40000,448.0394620671692 +213,7,80,30000,40500,449.64309860415494 +214,7,80,30000,41000,462.4490598612332 +215,7,80,30000,41500,481.6323506247537 +216,7,80,30500,38000,871.1163930229344 +217,7,80,30500,38500,771.1320563649375 +218,7,80,30500,39000,707.8872660015606 +219,7,80,30500,39500,672.6612145133173 +220,7,80,30500,40000,657.4974157809264 +221,7,80,30500,40500,656.0835852491216 +222,7,80,30500,41000,663.6006958125331 +223,7,80,30500,41500,676.460675405631 +224,7,80,31000,38000,1053.1852617390061 +225,7,80,31000,38500,984.3647109805877 +226,7,80,31000,39000,938.6158531749268 +227,7,80,31000,39500,911.4268280093535 +228,7,80,31000,40000,898.333365348419 +229,7,80,31000,40500,895.3996527486954 +230,7,80,31000,41000,899.3556288533885 +231,7,80,31000,41500,907.6180684887955 +232,7,80,31500,38000,1274.2255948763498 +233,7,80,31500,38500,1226.5236809533717 +234,7,80,31500,39000,1193.4538731398666 +235,7,80,31500,39500,1172.8105398345213 +236,7,80,31500,40000,1162.0692230240734 +237,7,80,31500,40500,1158.7461521476607 +238,7,80,31500,41000,1160.6173577210805 +239,7,80,31500,41500,1165.840315694716 +240,7,120,29000,38000,1325.2409732290193 +241,7,120,29000,38500,900.8063148840154 +242,7,120,29000,39000,629.9300352098937 +243,7,120,29000,39500,413.81648033893424 +244,7,120,29000,40000,257.3116751690404 +245,7,120,29000,40500,177.89217179438947 +246,7,120,29000,41000,151.58366848473491 +247,7,120,29000,41500,157.56967437251706 +248,7,120,29500,38000,1211.2807882170853 +249,7,120,29500,38500,956.936161969002 +250,7,120,29500,39000,753.3050086992201 +251,7,120,29500,39500,528.2452647799327 +252,7,120,29500,40000,382.62610532894917 +253,7,120,29500,40500,308.44199089882375 +254,7,120,29500,41000,280.3893024671524 +255,7,120,29500,41500,280.4028092582749 +256,7,120,30000,38000,1266.5740351143413 +257,7,120,30000,38500,1084.3028700477778 +258,7,120,30000,39000,834.2392498526193 +259,7,120,30000,39500,650.7560171314304 +260,7,120,30000,40000,537.7846910878052 +261,7,120,30000,40500,477.3001078155485 +262,7,120,30000,41000,451.6865380286754 +263,7,120,30000,41500,448.14911508024613 +264,7,120,30500,38000,1319.6603196780936 +265,7,120,30500,38500,1102.3027489012372 +266,7,120,30500,39000,931.2523583659847 +267,7,120,30500,39500,807.0833484596384 +268,7,120,30500,40000,727.4852710400268 +269,7,120,30500,40500,682.1437030344305 +270,7,120,30500,41000,660.7859329989657 +271,7,120,30500,41500,655.6001132492668 +272,7,120,31000,38000,1330.5306924865326 +273,7,120,31000,38500,1195.9190861202942 +274,7,120,31000,39000,1086.0328080422887 +275,7,120,31000,39500,1005.4160637517409 +276,7,120,31000,40000,951.2021706290612 +277,7,120,31000,40500,918.1457644271304 +278,7,120,31000,41000,901.0511005554887 +279,7,120,31000,41500,895.4599964465793 +280,7,120,31500,38000,1447.8365822059013 +281,7,120,31500,38500,1362.3417347939844 +282,7,120,31500,39000,1292.382727215108 +283,7,120,31500,39500,1239.1826828976662 +284,7,120,31500,40000,1201.6474412465277 +285,7,120,31500,40500,1177.5235955796813 +286,7,120,31500,41000,1164.1761722345295 +287,7,120,31500,41500,1158.9997785002718 +288,10,40,29000,38000,33.437068437082054 +289,10,40,29000,38500,58.471249815534996 +290,10,40,29000,39000,101.41937628542912 +291,10,40,29000,39500,153.80690200519626 +292,10,40,29000,40000,209.66451461551316 +293,10,40,29000,40500,265.03070792175197 +294,10,40,29000,41000,317.46079310177566 +295,10,40,29000,41500,365.59950388342645 +296,10,40,29500,38000,70.26818405688635 +297,10,40,29500,38500,87.96463718548947 +298,10,40,29500,39000,122.58188233160993 +299,10,40,29500,39500,166.2478945807132 +300,10,40,29500,40000,213.48669617414316 +301,10,40,29500,40500,260.67953961944477 +302,10,40,29500,41000,305.5877041218316 +303,10,40,29500,41500,346.95612213021155 +304,10,40,30000,38000,153.67588703371362 +305,10,40,30000,38500,164.07504103479005 +306,10,40,30000,39000,190.0800160661499 +307,10,40,30000,39500,224.61382980242837 +308,10,40,30000,40000,262.79232847382445 +309,10,40,30000,40500,301.38687703450415 +310,10,40,30000,41000,338.38536686093164 +311,10,40,30000,41500,372.6399011703545 +312,10,40,30500,38000,284.2936286531718 +313,10,40,30500,38500,288.4690608277705 +314,10,40,30500,39000,306.44667517621144 +315,10,40,30500,39500,332.20122250191986 +316,10,40,30500,40000,361.5566690083291 +317,10,40,30500,40500,391.72755224929614 +318,10,40,30500,41000,420.95317535960476 +319,10,40,30500,41500,448.2049230608669 +320,10,40,31000,38000,459.03140021766137 +321,10,40,31000,38500,458.71477027519967 +322,10,40,31000,39000,469.9910751800656 +323,10,40,31000,39500,488.05850105225426 +324,10,40,31000,40000,509.5204701455629 +325,10,40,31000,40500,532.0674969691778 +326,10,40,31000,41000,554.2088430693509 +327,10,40,31000,41500,575.0485839499048 +328,10,40,31500,38000,672.2476845983564 +329,10,40,31500,38500,669.2240508488649 +330,10,40,31500,39000,675.4956226836405 +331,10,40,31500,39500,687.447764319295 +332,10,40,31500,40000,702.4395430742891 +333,10,40,31500,40500,718.6279487347668 +334,10,40,31500,41000,734.793684592168 +335,10,40,31500,41500,750.1821072409286 +336,10,80,29000,38000,387.7617282731497 +337,10,80,29000,38500,195.33642612593002 +338,10,80,29000,39000,82.7306931465102 +339,10,80,29000,39500,35.13436471793541 +340,10,80,29000,40000,33.521138659248706 +341,10,80,29000,40500,61.47395975053128 +342,10,80,29000,41000,106.71403229340167 +343,10,80,29000,41500,160.56068704487473 +344,10,80,29500,38000,459.63404601804103 +345,10,80,29500,38500,258.7453720995899 +346,10,80,29500,39000,135.96435731320256 +347,10,80,29500,39500,80.2685095017944 +348,10,80,29500,40000,70.86302366453106 +349,10,80,29500,40500,90.43203026480438 +350,10,80,29500,41000,126.7844695901737 +351,10,80,29500,41500,171.63682876805044 +352,10,80,30000,38000,564.1463320344325 +353,10,80,30000,38500,360.75718124523866 +354,10,80,30000,39000,231.70119191254307 +355,10,80,30000,39500,170.74752201483128 +356,10,80,30000,40000,154.7149036950422 +357,10,80,30000,40500,166.10596450541493 +358,10,80,30000,41000,193.3351721194443 +359,10,80,30000,41500,228.78394172417038 +360,10,80,30500,38000,689.6797223218513 +361,10,80,30500,38500,484.8023695265838 +362,10,80,30500,39000,363.5979340028588 +363,10,80,30500,39500,304.67857102688225 +364,10,80,30500,40000,285.29210000833734 +365,10,80,30500,40500,290.0135917456113 +366,10,80,30500,41000,308.8672169492536 +367,10,80,30500,41500,335.3210332569182 +368,10,80,31000,38000,789.946106942773 +369,10,80,31000,38500,625.7722360026959 +370,10,80,31000,39000,528.6063264942235 +371,10,80,31000,39500,478.6863763478618 +372,10,80,31000,40000,459.5026243189753 +373,10,80,31000,40500,459.6982093164963 +374,10,80,31000,41000,471.6790024321937 +375,10,80,31000,41500,490.3034492109124 +376,10,80,31500,38000,912.3540488244158 +377,10,80,31500,38500,798.2135101409633 +378,10,80,31500,39000,727.746684419146 +379,10,80,31500,39500,689.0119464356724 +380,10,80,31500,40000,672.0757202772029 +381,10,80,31500,40500,669.678339553036 +382,10,80,31500,41000,676.5761221409929 +383,10,80,31500,41500,688.9934449650118 +384,10,120,29000,38000,1155.1165164624408 +385,10,120,29000,38500,840.2641727088946 +386,10,120,29000,39000,506.9102636732852 +387,10,120,29000,39500,265.5278912452038 +388,10,120,29000,40000,116.39516513179322 +389,10,120,29000,40500,45.2088092745619 +390,10,120,29000,41000,30.22267557153353 +391,10,120,29000,41500,51.06063746392809 +392,10,120,29500,38000,1343.7868459826054 +393,10,120,29500,38500,977.9852373227346 +394,10,120,29500,39000,594.632756549817 +395,10,120,29500,39500,346.2478773329187 +396,10,120,29500,40000,180.23082247413407 +397,10,120,29500,40500,95.81649989178923 +398,10,120,29500,41000,71.0837801649128 +399,10,120,29500,41500,82.84289818279714 +400,10,120,30000,38000,1532.9333545384934 +401,10,120,30000,38500,1012.2223350568845 +402,10,120,30000,39000,688.4884716222766 +403,10,120,30000,39500,464.6206903113392 +404,10,120,30000,40000,283.5644748300334 +405,10,120,30000,40500,190.27593217865416 +406,10,120,30000,41000,158.0192279691727 +407,10,120,30000,41500,161.3611926772337 +408,10,120,30500,38000,1349.3785399811063 +409,10,120,30500,38500,1014.785480110738 +410,10,120,30500,39000,843.0316833766408 +411,10,120,30500,39500,589.4543896730125 +412,10,120,30500,40000,412.3358512291996 +413,10,120,30500,40500,324.11715620464133 +414,10,120,30500,41000,290.17588242984766 +415,10,120,30500,41500,287.56857384673356 +416,10,120,31000,38000,1328.0973931040146 +417,10,120,31000,38500,1216.5659656437845 +418,10,120,31000,39000,928.4831767181619 +419,10,120,31000,39500,700.3115484040329 +420,10,120,31000,40000,565.0876352458171 +421,10,120,31000,40500,494.44016026435037 +422,10,120,31000,41000,464.38005437182983 +423,10,120,31000,41500,458.7614573733091 +424,10,120,31500,38000,1473.1154650008834 +425,10,120,31500,38500,1195.943614951571 +426,10,120,31500,39000,990.2486604382486 +427,10,120,31500,39500,843.1390407497395 +428,10,120,31500,40000,751.2746391170706 +429,10,120,31500,40500,700.215375503209 +430,10,120,31500,41000,676.1585052687219 +431,10,120,31500,41500,669.5907920932743 +432,13,40,29000,38000,49.96352152045025 +433,13,40,29000,38500,83.75104994958261 +434,13,40,29000,39000,136.8176091795391 +435,13,40,29000,39500,199.91486685466407 +436,13,40,29000,40000,266.4367154860076 +437,13,40,29000,40500,331.97224579940524 +438,13,40,29000,41000,393.8001583706036 +439,13,40,29000,41500,450.42425363084493 +440,13,40,29500,38000,29.775721038786923 +441,13,40,29500,38500,57.37673742631121 +442,13,40,29500,39000,103.49161398239501 +443,13,40,29500,39500,159.3058253852367 +444,13,40,29500,40000,218.60083223764073 +445,13,40,29500,40500,277.2507278183831 +446,13,40,29500,41000,332.7141278886951 +447,13,40,29500,41500,383.58832292300576 +448,13,40,30000,38000,47.72263852005472 +449,13,40,30000,38500,68.07581028940402 +450,13,40,30000,39000,106.13974628945516 +451,13,40,30000,39500,153.58449949683063 +452,13,40,30000,40000,204.62393623358633 +453,13,40,30000,40500,255.44513025602419 +454,13,40,30000,41000,303.69954914051766 +455,13,40,30000,41500,348.0803709720354 +456,13,40,30500,38000,110.9331168284094 +457,13,40,30500,38500,123.63361262704746 +458,13,40,30500,39000,153.02654433825705 +459,13,40,30500,39500,191.40769947472756 +460,13,40,30500,40000,233.503841403055 +461,13,40,30500,40500,275.8557790922913 +462,13,40,30500,41000,316.32529882763697 +463,13,40,30500,41500,353.7060432094809 +464,13,40,31000,38000,221.90608823073939 +465,13,40,31000,38500,227.67026441593657 +466,13,40,31000,39000,248.62107049869064 +467,13,40,31000,39500,277.9507605389158 +468,13,40,31000,40000,311.0267471957685 +469,13,40,31000,40500,344.8024031161673 +470,13,40,31000,41000,377.3761144228052 +471,13,40,31000,41500,407.6529635071056 +472,13,40,31500,38000,378.8738382757093 +473,13,40,31500,38500,379.39748335944216 +474,13,40,31500,39000,393.01223361732553 +475,13,40,31500,39500,414.10238059122855 +476,13,40,31500,40000,438.8024282436204 +477,13,40,31500,40500,464.5348067190265 +478,13,40,31500,41000,489.6621039898805 +479,13,40,31500,41500,513.2163939332803 +480,13,80,29000,38000,364.387588581215 +481,13,80,29000,38500,184.2902007673634 +482,13,80,29000,39000,81.57192155036655 +483,13,80,29000,39500,42.54811210095659 +484,13,80,29000,40000,49.897338772663076 +485,13,80,29000,40500,87.84229516509882 +486,13,80,29000,41000,143.85451969447664 +487,13,80,29000,41500,208.71467984917848 +488,13,80,29500,38000,382.5794635435733 +489,13,80,29500,38500,188.38619353711718 +490,13,80,29500,39000,75.75749359688277 +491,13,80,29500,39500,29.27891251986562 +492,13,80,29500,40000,29.794874961934568 +493,13,80,29500,40500,60.654888662698205 +494,13,80,29500,41000,109.25801388824325 +495,13,80,29500,41500,166.6311093454692 +496,13,80,30000,38000,448.97795526074816 +497,13,80,30000,38500,238.44530107604737 +498,13,80,30000,39000,112.34545890264337 +499,13,80,30000,39500,56.125871791222835 +500,13,80,30000,40000,48.29987461781518 +501,13,80,30000,40500,70.7900626637678 +502,13,80,30000,41000,110.76865376691964 +503,13,80,30000,41500,159.50197316936024 +504,13,80,30500,38000,547.7818730461195 +505,13,80,30500,38500,332.92604070423494 +506,13,80,30500,39000,193.80760050280742 +507,13,80,30500,39500,128.3457644087917 +508,13,80,30500,40000,112.23915895822442 +509,13,80,30500,40500,125.96369396512564 +510,13,80,30500,41000,156.67918617660013 +511,13,80,30500,41500,196.05195109523765 +512,13,80,31000,38000,682.8591931963246 +513,13,80,31000,38500,457.56562267948556 +514,13,80,31000,39000,313.6380169123524 +515,13,80,31000,39500,245.13531819580908 +516,13,80,31000,40000,223.54473391202873 +517,13,80,31000,40500,229.60752111202834 +518,13,80,31000,41000,251.42377424735136 +519,13,80,31000,41500,281.48720903016886 +520,13,80,31500,38000,807.925638050234 +521,13,80,31500,38500,588.686585641994 +522,13,80,31500,39000,464.0488586698228 +523,13,80,31500,39500,402.69214492641095 +524,13,80,31500,40000,380.13626165363934 +525,13,80,31500,40500,380.8064948609387 +526,13,80,31500,41000,395.05186915919086 +527,13,80,31500,41500,416.70193045600774 +528,13,120,29000,38000,1068.8279454397398 +529,13,120,29000,38500,743.0012805963486 +530,13,120,29000,39000,451.2538301167544 +531,13,120,29000,39500,235.4154251166075 +532,13,120,29000,40000,104.73720814447498 +533,13,120,29000,40500,46.91983990671749 +534,13,120,29000,41000,42.81092192562316 +535,13,120,29000,41500,74.33530639171506 +536,13,120,29500,38000,1133.1178848710972 +537,13,120,29500,38500,824.0745323788527 +538,13,120,29500,39000,499.10867111401996 +539,13,120,29500,39500,256.1626809904186 +540,13,120,29500,40000,107.68599585294751 +541,13,120,29500,40500,38.18533662516749 +542,13,120,29500,41000,25.499608203619154 +543,13,120,29500,41500,49.283537699300375 +544,13,120,30000,38000,1292.409871290162 +545,13,120,30000,38500,994.669572829704 +546,13,120,30000,39000,598.9783697712826 +547,13,120,30000,39500,327.47348408537925 +548,13,120,30000,40000,156.82634841081907 +549,13,120,30000,40500,71.30833688875883 +550,13,120,30000,41000,47.72389750130817 +551,13,120,30000,41500,62.1982461882982 +552,13,120,30500,38000,1585.8797221278146 +553,13,120,30500,38500,1144.66688416451 +554,13,120,30500,39000,692.6651441690645 +555,13,120,30500,39500,441.98837639874046 +556,13,120,30500,40000,251.56311435857728 +557,13,120,30500,40500,149.79670413140468 +558,13,120,30500,41000,115.52645596043719 +559,13,120,30500,41500,120.44019473389324 +560,13,120,31000,38000,1702.7625866892163 +561,13,120,31000,38500,1071.7854750250656 +562,13,120,31000,39000,807.8943299034604 +563,13,120,31000,39500,588.672223513561 +564,13,120,31000,40000,376.44658358671404 +565,13,120,31000,40500,269.2159719426485 +566,13,120,31000,41000,229.41660529009877 +567,13,120,31000,41500,226.78274707181976 +568,13,120,31500,38000,1331.3523701291767 +569,13,120,31500,38500,1151.2055268669133 +570,13,120,31500,39000,1006.811285091974 +571,13,120,31500,39500,702.0053094629535 +572,13,120,31500,40000,515.9081891614829 +573,13,120,31500,40500,423.8652275555525 +574,13,120,31500,41000,386.4939696097151 +575,13,120,31500,41500,379.8118453367429 +576,16,40,29000,38000,106.1025746852808 +577,16,40,29000,38500,145.32590128581407 +578,16,40,29000,39000,204.74804378224422 +579,16,40,29000,39500,274.6339266648551 +580,16,40,29000,40000,347.9667393938497 +581,16,40,29000,40500,420.03753452490974 +582,16,40,29000,41000,487.9353932879741 +583,16,40,29000,41500,550.0623063219693 +584,16,40,29500,38000,54.65040870471303 +585,16,40,29500,38500,88.94089091627293 +586,16,40,29500,39000,142.72223808288405 +587,16,40,29500,39500,206.63598763907422 +588,16,40,29500,40000,273.99851593521134 +589,16,40,29500,40500,340.34861536649436 +590,16,40,29500,41000,402.935270882596 +591,16,40,29500,41500,460.2471155081633 +592,16,40,30000,38000,29.788548081995298 +593,16,40,30000,38500,57.96323252610644 +594,16,40,30000,39000,104.92815906834525 +595,16,40,30000,39500,161.71867032726158 +596,16,40,30000,40000,222.01677586338877 +597,16,40,30000,40500,281.6349465235367 +598,16,40,30000,41000,337.99683241119567 +599,16,40,30000,41500,389.68271710858414 +600,16,40,30500,38000,42.06569536892785 +601,16,40,30500,38500,62.95145274276575 +602,16,40,30500,39000,101.93860830594608 +603,16,40,30500,39500,150.47910837525734 +604,16,40,30500,40000,202.65388851823258 +605,16,40,30500,40500,254.5724108541227 +606,16,40,30500,41000,303.84403622726694 +607,16,40,30500,41500,349.1422884543064 +608,16,40,31000,38000,99.21707896667829 +609,16,40,31000,38500,112.24153596941301 +610,16,40,31000,39000,142.5186177618655 +611,16,40,31000,39500,182.02836955332134 +612,16,40,31000,40000,225.3201896575212 +613,16,40,31000,40500,268.83705389232614 +614,16,40,31000,41000,310.3895932135811 +615,16,40,31000,41500,348.7480165565453 +616,16,40,31500,38000,204.30418825821732 +617,16,40,31500,38500,210.0759235359138 +618,16,40,31500,39000,231.7643258544752 +619,16,40,31500,39500,262.1512494310348 +620,16,40,31500,40000,296.3864127264238 +621,16,40,31500,40500,331.30743171999035 +622,16,40,31500,41000,364.95322314895554 +623,16,40,31500,41500,396.20142191205844 +624,16,80,29000,38000,399.5975649320935 +625,16,80,29000,38500,225.6318269911425 +626,16,80,29000,39000,127.97354075513151 +627,16,80,29000,39500,93.73584101549991 +628,16,80,29000,40000,106.43084032022394 +629,16,80,29000,40500,150.51245762256931 +630,16,80,29000,41000,213.24213500046466 +631,16,80,29000,41500,285.0426423013882 +632,16,80,29500,38000,371.37706087096393 +633,16,80,29500,38500,189.77150413822454 +634,16,80,29500,39000,86.22375488959844 +635,16,80,29500,39500,46.98714814001572 +636,16,80,29500,40000,54.596900621760675 +637,16,80,29500,40500,93.12033833747024 +638,16,80,29500,41000,149.89341227947025 +639,16,80,29500,41500,215.5937000584367 +640,16,80,30000,38000,388.43657991253195 +641,16,80,30000,38500,190.77121362008674 +642,16,80,30000,39000,76.28535232335287 +643,16,80,30000,39500,29.152860363695716 +644,16,80,30000,40000,29.820972887404942 +645,16,80,30000,40500,61.320203047752464 +646,16,80,30000,41000,110.82086782062603 +647,16,80,30000,41500,169.197767615573 +648,16,80,30500,38000,458.8964339917103 +649,16,80,30500,38500,239.547928886725 +650,16,80,30500,39000,109.02338779317503 +651,16,80,30500,39500,50.888746196140914 +652,16,80,30500,40000,42.73606982375976 +653,16,80,30500,40500,65.75935122724029 +654,16,80,30500,41000,106.68884313872147 +655,16,80,30500,41500,156.54100549486617 +656,16,80,31000,38000,561.7385153195615 +657,16,80,31000,38500,335.5692026144635 +658,16,80,31000,39000,188.0383015831574 +659,16,80,31000,39500,118.2318539104416 +660,16,80,31000,40000,100.81000168801492 +661,16,80,31000,40500,114.72014539486217 +662,16,80,31000,41000,146.2992492326178 +663,16,80,31000,41500,186.8074429488408 +664,16,80,31500,38000,697.9937997454152 +665,16,80,31500,38500,466.42234442578484 +666,16,80,31500,39000,306.52125608515166 +667,16,80,31500,39500,230.54692639209762 +668,16,80,31500,40000,206.461121102699 +669,16,80,31500,40500,212.23429887269359 +670,16,80,31500,41000,234.70913795495554 +671,16,80,31500,41500,265.8143069252357 +672,16,120,29000,38000,1085.688903883652 +673,16,120,29000,38500,750.2887000017752 +674,16,120,29000,39000,469.92662852990964 +675,16,120,29000,39500,267.1560282754928 +676,16,120,29000,40000,146.06299930062625 +677,16,120,29000,40500,95.28836772053619 +678,16,120,29000,41000,97.41466545178946 +679,16,120,29000,41500,135.3804131941845 +680,16,120,29500,38000,1079.5576154477903 +681,16,120,29500,38500,751.2932384998761 +682,16,120,29500,39000,458.27083477307207 +683,16,120,29500,39500,240.9658024131812 +684,16,120,29500,40000,109.3801465044384 +685,16,120,29500,40500,51.274139057659724 +686,16,120,29500,41000,47.36446629605638 +687,16,120,29500,41500,79.42944320845996 +688,16,120,30000,38000,1139.3792936518537 +689,16,120,30000,38500,833.7979589668842 +690,16,120,30000,39000,507.805443202025 +691,16,120,30000,39500,259.93892964607977 +692,16,120,30000,40000,108.7341499557062 +693,16,120,30000,40500,38.152937143498605 +694,16,120,30000,41000,25.403985123518716 +695,16,120,30000,41500,49.72822589160786 +696,16,120,30500,38000,1285.0396277304772 +697,16,120,30500,38500,1025.254169031627 +698,16,120,30500,39000,622.5890550779666 +699,16,120,30500,39500,333.3353043756717 +700,16,120,30500,40000,155.70268128051293 +701,16,120,30500,40500,66.84125446522368 +702,16,120,30500,41000,42.25187049753978 +703,16,120,30500,41500,56.98314898830595 +704,16,120,31000,38000,1595.7993459811262 +705,16,120,31000,38500,1252.8886556470425 +706,16,120,31000,39000,731.4408383874198 +707,16,120,31000,39500,451.0090473423308 +708,16,120,31000,40000,251.5086563526081 +709,16,120,31000,40500,141.8915050063955 +710,16,120,31000,41000,104.67474675582574 +711,16,120,31000,41500,109.1609567535697 +712,16,120,31500,38000,1942.3896021770768 +713,16,120,31500,38500,1197.207050908449 +714,16,120,31500,39000,812.6818768064074 +715,16,120,31500,39500,611.45532452889 +716,16,120,31500,40000,380.63642711770643 +717,16,120,31500,40500,258.5514125337487 +718,16,120,31500,41000,213.48518421250665 +719,16,120,31500,41500,209.58134396574906 +720,19,40,29000,38000,169.3907733115706 +721,19,40,29000,38500,212.23331960093145 +722,19,40,29000,39000,275.9376503672959 +723,19,40,29000,39500,350.4301397081139 +724,19,40,29000,40000,428.40863665493924 +725,19,40,29000,40500,504.955113902399 +726,19,40,29000,41000,577.023450987656 +727,19,40,29000,41500,642.9410032211753 +728,19,40,29500,38000,102.40889356493292 +729,19,40,29500,38500,141.19036226103668 +730,19,40,29500,39000,200.19333708701748 +731,19,40,29500,39500,269.6750686488757 +732,19,40,29500,40000,342.6217886299377 +733,19,40,29500,40500,414.33044375626207 +734,19,40,29500,41000,481.89521316730713 +735,19,40,29500,41500,543.7211700546151 +736,19,40,30000,38000,51.95330426445395 +737,19,40,30000,38500,85.69656829127965 +738,19,40,30000,39000,138.98376466247876 +739,19,40,30000,39500,202.43251598105033 +740,19,40,30000,40000,269.3557903452929 +741,19,40,30000,40500,335.2960133312316 +742,19,40,30000,41000,397.50658847538665 +743,19,40,30000,41500,454.47903112410967 +744,19,40,30500,38000,28.864802790801026 +745,19,40,30500,38500,56.32899754732796 +746,19,40,30500,39000,102.69825523352162 +747,19,40,30500,39500,158.95118263535466 +748,19,40,30500,40000,218.75241957992617 +749,19,40,30500,40500,277.9122290233915 +750,19,40,30500,41000,333.8561815041273 +751,19,40,30500,41500,385.1662652901447 +752,19,40,31000,38000,43.72359701781447 +753,19,40,31000,38500,63.683967347844224 +754,19,40,31000,39000,101.95579433282329 +755,19,40,31000,39500,149.8826019475827 +756,19,40,31000,40000,201.50605279789198 +757,19,40,31000,40500,252.92391570754876 +758,19,40,31000,41000,301.7431453727685 +759,19,40,31000,41500,346.6368192781496 +760,19,40,31500,38000,104.05710998615942 +761,19,40,31500,38500,115.95783594434451 +762,19,40,31500,39000,145.42181873662554 +763,19,40,31500,39500,184.26373455825217 +764,19,40,31500,40000,226.97066340897095 +765,19,40,31500,40500,269.96403356902357 +766,19,40,31500,41000,311.04753558871505 +767,19,40,31500,41500,348.98866332680115 +768,19,80,29000,38000,453.1314944429312 +769,19,80,29000,38500,281.24067760117225 +770,19,80,29000,39000,185.83730378881882 +771,19,80,29000,39500,154.25726305915472 +772,19,80,29000,40000,170.2912737797755 +773,19,80,29000,40500,218.38979299191152 +774,19,80,29000,41000,285.604024444273 +775,19,80,29000,41500,362.0858325427657 +776,19,80,29500,38000,400.06299682217264 +777,19,80,29500,38500,224.41725666435008 +778,19,80,29500,39000,125.58476107530382 +779,19,80,29500,39500,90.55733834394478 +780,19,80,29500,40000,102.67519971027264 +781,19,80,29500,40500,146.27807815967392 +782,19,80,29500,41000,208.57372904155937 +783,19,80,29500,41500,279.9669583078214 +784,19,80,30000,38000,376.1594584816549 +785,19,80,30000,38500,191.30452808298463 +786,19,80,30000,39000,85.63116084217559 +787,19,80,30000,39500,45.10487847849711 +788,19,80,30000,40000,51.88389644342952 +789,19,80,30000,40500,89.78942817703852 +790,19,80,30000,41000,146.0393555385696 +791,19,80,30000,41500,211.26567367707352 +792,19,80,30500,38000,401.874315275947 +793,19,80,30500,38500,197.55305366608133 +794,19,80,30500,39000,79.00348967857379 +795,19,80,30500,39500,29.602719961568614 +796,19,80,30500,40000,28.980451378502487 +797,19,80,30500,40500,59.63541802023186 +798,19,80,30500,41000,108.48607655362268 +799,19,80,30500,41500,166.30589286399507 +800,19,80,31000,38000,484.930958445979 +801,19,80,31000,38500,254.27552635537404 +802,19,80,31000,39000,116.75543721560439 +803,19,80,31000,39500,54.77547840250418 +804,19,80,31000,40000,44.637472658824976 +805,19,80,31000,40500,66.50466903927668 +806,19,80,31000,41000,106.62737262508298 +807,19,80,31000,41500,155.8310688191254 +808,19,80,31500,38000,595.6094306603337 +809,19,80,31500,38500,359.60040819463063 +810,19,80,31500,39000,201.85328967228585 +811,19,80,31500,39500,126.24442464793601 +812,19,80,31500,40000,106.07388975142673 +813,19,80,31500,40500,118.52358345403363 +814,19,80,31500,41000,149.1597537162607 +815,19,80,31500,41500,188.94964975523197 +816,19,120,29000,38000,1133.9213841599772 +817,19,120,29000,38500,793.9759807804692 +818,19,120,29000,39000,516.5580425563733 +819,19,120,29000,39500,318.60172051726147 +820,19,120,29000,40000,201.662212274693 +821,19,120,29000,40500,154.47522945829064 +822,19,120,29000,41000,160.28049502033574 +823,19,120,29000,41500,202.35345983501588 +824,19,120,29500,38000,1091.6343400395158 +825,19,120,29500,38500,754.9332443184217 +826,19,120,29500,39000,472.1777992591152 +827,19,120,29500,39500,267.03951846894995 +828,19,120,29500,40000,144.25558152688114 +829,19,120,29500,40500,92.40384156679512 +830,19,120,29500,41000,93.81833253459942 +831,19,120,29500,41500,131.24753560710644 +832,19,120,30000,38000,1092.719296892266 +833,19,120,30000,38500,764.7065490850255 +834,19,120,30000,39000,467.2268758064373 +835,19,120,30000,39500,244.9367732985332 +836,19,120,30000,40000,110.00996333393202 +837,19,120,30000,40500,49.96381544207811 +838,19,120,30000,41000,44.9298739569088 +839,19,120,30000,41500,76.25447129089613 +840,19,120,30500,38000,1160.6160120981158 +841,19,120,30500,38500,865.5953188304933 +842,19,120,30500,39000,531.1657093741892 +843,19,120,30500,39500,271.98520008106277 +844,19,120,30500,40000,114.03616090967407 +845,19,120,30500,40500,39.74252227099571 +846,19,120,30500,41000,25.07176465285551 +847,19,120,30500,41500,48.298794094852724 +848,19,120,31000,38000,1304.8870694342509 +849,19,120,31000,38500,1089.6854636757826 +850,19,120,31000,39000,668.6632735260521 +851,19,120,31000,39500,356.7751012890747 +852,19,120,31000,40000,168.32491564142487 +853,19,120,31000,40500,72.82648063377391 +854,19,120,31000,41000,45.02326687759286 +855,19,120,31000,41500,58.13111530831655 +856,19,120,31500,38000,1645.2697164013964 +857,19,120,31500,38500,1373.859712069864 +858,19,120,31500,39000,787.3948673670299 +859,19,120,31500,39500,483.60546305948367 +860,19,120,31500,40000,273.4285373433001 +861,19,120,31500,40500,153.21079535396908 +862,19,120,31500,41000,111.21299419905313 +863,19,120,31500,41500,113.52006337929113 +864,22,40,29000,38000,229.2032513971666 +865,22,40,29000,38500,274.65023153674116 +866,22,40,29000,39000,341.4424739822062 +867,22,40,29000,39500,419.2624324130753 +868,22,40,29000,40000,500.6022690006133 +869,22,40,29000,40500,580.3923016374031 +870,22,40,29000,41000,655.4874207991389 +871,22,40,29000,41500,724.1595537770351 +872,22,40,29500,38000,155.45206306046595 +873,22,40,29500,38500,197.41588482427002 +874,22,40,29500,39000,260.1641484982308 +875,22,40,29500,39500,333.666918810689 +876,22,40,29500,40000,410.66541588422854 +877,22,40,29500,40500,486.276072112155 +878,22,40,29500,41000,557.4760464927683 +879,22,40,29500,41500,622.6057687448293 +880,22,40,30000,38000,90.70026588811803 +881,22,40,30000,38500,128.41239603755494 +882,22,40,30000,39000,186.27261386900233 +883,22,40,30000,39500,254.5802373859711 +884,22,40,30000,40000,326.3686182341553 +885,22,40,30000,40500,396.9735001502319 +886,22,40,30000,41000,463.5155278718613 +887,22,40,30000,41500,524.414569320113 +888,22,40,30500,38000,44.551475763397946 +889,22,40,30500,38500,76.95264448905411 +890,22,40,30500,39000,128.85898727872572 +891,22,40,30500,39500,190.91422001003792 +892,22,40,30500,40000,256.4755613806196 +893,22,40,30500,40500,321.125224208803 +894,22,40,30500,41000,382.14434919800453 +895,22,40,30500,41500,438.03974322333033 +896,22,40,31000,38000,28.101321546315717 +897,22,40,31000,38500,53.867829756398805 +898,22,40,31000,39000,98.57619184859544 +899,22,40,31000,39500,153.19473192134507 +900,22,40,31000,40000,211.4202434313414 +901,22,40,31000,40500,269.09905982026265 +902,22,40,31000,41000,323.68306330754416 +903,22,40,31000,41500,373.76836451736045 +904,22,40,31500,38000,51.648288279447364 +905,22,40,31500,38500,69.56074881661863 +906,22,40,31500,39000,105.91402675097291 +907,22,40,31500,39500,151.99456204656389 +908,22,40,31500,40000,201.85995274525234 +909,22,40,31500,40500,251.63807959916412 +910,22,40,31500,41000,298.9593498669657 +911,22,40,31500,41500,342.50888994628025 +912,22,80,29000,38000,507.5440336860194 +913,22,80,29000,38500,336.42019672232965 +914,22,80,29000,39000,242.21016116765423 +915,22,80,29000,39500,212.33396533224905 +916,22,80,29000,40000,230.67632355958136 +917,22,80,29000,40500,281.6224662955561 +918,22,80,29000,41000,352.0457411487133 +919,22,80,29000,41500,431.89288175778637 +920,22,80,29500,38000,443.2889283037078 +921,22,80,29500,38500,270.0648237630224 +922,22,80,29500,39000,173.57666711629645 +923,22,80,29500,39500,141.06258420240613 +924,22,80,29500,40000,156.18412870159142 +925,22,80,29500,40500,203.33105261575707 +926,22,80,29500,41000,269.5552387411201 +927,22,80,29500,41500,345.03801326123767 +928,22,80,30000,38000,395.34177505602497 +929,22,80,30000,38500,217.11094192826982 +930,22,80,30000,39000,116.38535634181476 +931,22,80,30000,39500,79.94742924888467 +932,22,80,30000,40000,90.84706550421288 +933,22,80,30000,40500,133.26308067939766 +934,22,80,30000,41000,194.36064414396228 +935,22,80,30000,41500,264.56059537656466 +936,22,80,30500,38000,382.0341866812038 +937,22,80,30500,38500,191.65621311671836 +938,22,80,30500,39000,82.3318677587146 +939,22,80,30500,39500,39.44606931321677 +940,22,80,30500,40000,44.476166488763134 +941,22,80,30500,40500,80.84561981845566 +942,22,80,30500,41000,135.62459431793735 +943,22,80,30500,41500,199.42208168600175 +944,22,80,31000,38000,425.5181957619983 +945,22,80,31000,38500,210.2667219741389 +946,22,80,31000,39000,84.97041062888985 +947,22,80,31000,39500,31.593073529038755 +948,22,80,31000,40000,28.407154164211214 +949,22,80,31000,40500,57.05446633976857 +950,22,80,31000,41000,104.10423883907688 +951,22,80,31000,41500,160.23135976433713 +952,22,80,31500,38000,527.5015417150911 +953,22,80,31500,38500,282.29650611769665 +954,22,80,31500,39000,134.62881845323489 +955,22,80,31500,39500,66.62736532046851 +956,22,80,31500,40000,52.9918858786988 +957,22,80,31500,40500,72.36913743145999 +958,22,80,31500,41000,110.38003828747726 +959,22,80,31500,41500,157.65470091455973 +960,22,120,29000,38000,1186.823326813257 +961,22,120,29000,38500,844.3317816964005 +962,22,120,29000,39000,567.7367986440256 +963,22,120,29000,39500,371.79782508970567 +964,22,120,29000,40000,256.9261857702517 +965,22,120,29000,40500,211.85466060592006 +966,22,120,29000,41000,220.09534855737033 +967,22,120,29000,41500,265.02731793490034 +968,22,120,29500,38000,1128.4568915685559 +969,22,120,29500,38500,787.7709648712951 +970,22,120,29500,39000,508.4832626962424 +971,22,120,29500,39500,308.52654841064975 +972,22,120,29500,40000,190.01030358402707 +973,22,120,29500,40500,141.62663282114926 +974,22,120,29500,41000,146.40704203984612 +975,22,120,29500,41500,187.48734389188584 +976,22,120,30000,38000,1094.7007205604846 +977,22,120,30000,38500,757.7313528729464 +978,22,120,30000,39000,471.282561364766 +979,22,120,30000,39500,262.0412520036699 +980,22,120,30000,40000,136.26956239282435 +981,22,120,30000,40500,82.4268827471484 +982,22,120,30000,41000,82.3695177584498 +983,22,120,30000,41500,118.51210034475737 +984,22,120,30500,38000,1111.0872182758205 +985,22,120,30500,38500,787.2204655558988 +986,22,120,30500,39000,481.85960605002055 +987,22,120,30500,39500,250.28740868446397 +988,22,120,30500,40000,109.21968920710272 +989,22,120,30500,40500,45.51600269221681 +990,22,120,30500,41000,38.172157811051115 +991,22,120,30500,41500,67.73748641348168 +992,22,120,31000,38000,1193.3958874354898 +993,22,120,31000,38500,923.0731791194576 +994,22,120,31000,39000,573.4457650536078 +995,22,120,31000,39500,294.2980811757103 +996,22,120,31000,40000,124.86249624679849 +997,22,120,31000,40500,43.948524347749846 +998,22,120,31000,41000,25.582084045731808 +999,22,120,31000,41500,46.36268252714472 +1000,22,120,31500,38000,1336.0993444856913 +1001,22,120,31500,38500,1194.893001664831 +1002,22,120,31500,39000,740.6584250286721 +1003,22,120,31500,39500,397.18127104230757 +1004,22,120,31500,40000,194.20390582893873 +1005,22,120,31500,40500,88.22588964369922 +1006,22,120,31500,41000,54.97797247760634 +1007,22,120,31500,41500,64.88195101638016 diff --git a/pyomo/contrib/parmest/examples/semibatch/semibatch.py b/pyomo/contrib/parmest/examples/semibatch/semibatch.py index 8cda262c019..6762531a338 100644 --- a/pyomo/contrib/parmest/examples/semibatch/semibatch.py +++ b/pyomo/contrib/parmest/examples/semibatch/semibatch.py @@ -34,11 +34,20 @@ def generate_model(data): + # if data is a file name, then load file first + if isinstance(data, str): + file_name = data + try: + with open(file_name, "r") as infile: + data = json.load(infile) + except: + raise RuntimeError(f"Could not read {file_name} as json") + # unpack and fix the data - cameastemp = data['Ca_meas'] - cbmeastemp = data['Cb_meas'] - ccmeastemp = data['Cc_meas'] - trmeastemp = data['Tr_meas'] + cameastemp = data["Ca_meas"] + cbmeastemp = data["Cb_meas"] + ccmeastemp = data["Cc_meas"] + trmeastemp = data["Tr_meas"] cameas = {} cbmeas = {} @@ -79,9 +88,9 @@ def generate_model(data): m.Vc = Param(initialize=0.07) # m^3 m.rhow = Param(initialize=700.0) # kg/m^3 m.cpw = Param(initialize=3.1) # kJ/kg/K - m.Ca0 = Param(initialize=data['Ca0']) # kmol/m^3) - m.Cb0 = Param(initialize=data['Cb0']) # kmol/m^3) - m.Cc0 = Param(initialize=data['Cc0']) # kmol/m^3) + m.Ca0 = Param(initialize=data["Ca0"]) # kmol/m^3) + m.Cb0 = Param(initialize=data["Cb0"]) # kmol/m^3) + m.Cc0 = Param(initialize=data["Cc0"]) # kmol/m^3) m.Tr0 = Param(initialize=300.0) # K m.Vr0 = Param(initialize=1.0) # m^3 @@ -92,9 +101,9 @@ def generate_model(data): # def _initTc(m, t): if t < 10800: - return data['Tc1'] + return data["Tc1"] else: - return data['Tc2'] + return data["Tc2"] m.Tc = Param( m.time, initialize=_initTc, default=_initTc @@ -102,9 +111,9 @@ def _initTc(m, t): def _initFa(m, t): if t < 10800: - return data['Fa1'] + return data["Fa1"] else: - return data['Fa2'] + return data["Fa2"] m.Fa = Param( m.time, initialize=_initFa, default=_initFa @@ -230,7 +239,7 @@ def AllMeasurements(m): ) def MissingMeasurements(m): - if data['experiment'] == 1: + if data["experiment"] == 1: return sum( (m.Ca[t] - m.Ca_meas[t]) ** 2 + (m.Cb[t] - m.Cb_meas[t]) ** 2 @@ -238,7 +247,7 @@ def MissingMeasurements(m): + (m.Tr[t] - m.Tr_meas[t]) ** 2 for t in m.measT ) - elif data['experiment'] == 2: + elif data["experiment"] == 2: return sum((m.Tr[t] - m.Tr_meas[t]) ** 2 for t in m.measT) else: return sum( @@ -254,7 +263,7 @@ def total_cost_rule(model): m.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize) # Discretize model - disc = TransformationFactory('dae.collocation') + disc = TransformationFactory("dae.collocation") disc.apply_to(m, nfe=20, ncp=4) return m @@ -262,17 +271,17 @@ def total_cost_rule(model): def main(): # Data loaded from files file_dirname = dirname(abspath(str(__file__))) - file_name = abspath(join(file_dirname, 'exp2.out')) - with open(file_name, 'r') as infile: + file_name = abspath(join(file_dirname, "exp2.out")) + with open(file_name, "r") as infile: data = json.load(infile) - data['experiment'] = 2 + data["experiment"] = 2 model = generate_model(data) - solver = SolverFactory('ipopt') + solver = SolverFactory("ipopt") solver.solve(model) - print('k1 = ', model.k1()) - print('E1 = ', model.E1()) + print("k1 = ", model.k1()) + print("E1 = ", model.E1()) -if __name__ == '__main__': +if __name__ == "__main__": main() diff --git a/pyomo/contrib/parmest/graphics.py b/pyomo/contrib/parmest/graphics.py index f01622d2d17..b8dfa243b9a 100644 --- a/pyomo/contrib/parmest/graphics.py +++ b/pyomo/contrib/parmest/graphics.py @@ -29,7 +29,7 @@ # (e.g. python 3.5) get released that are either broken not # compatible, resulting in a SyntaxError sns, seaborn_available = attempt_import( - 'seaborn', catch_exceptions=(ImportError, SyntaxError) + "seaborn", catch_exceptions=(ImportError, SyntaxError) ) imports_available = ( @@ -93,17 +93,17 @@ def _get_data_slice(xvar, yvar, columns, data, theta_star): temp[col] = temp[col] + data[col].std() data = pd.concat([data, temp], ignore_index=True) - data_slice['obj'] = scipy.interpolate.griddata( + data_slice["obj"] = scipy.interpolate.griddata( np.array(data[columns]), - np.array(data[['obj']]), + np.array(data[["obj"]]), np.array(data_slice[columns]), - method='linear', + method="linear", rescale=True, ) X = data_slice[xvar] Y = data_slice[yvar] - Z = data_slice['obj'] + Z = data_slice["obj"] return X, Y, Z @@ -178,11 +178,11 @@ def _add_obj_contour(x, y, color, columns, data, theta_star, label=None): X, Y, Z = _get_data_slice(xvar, yvar, columns, data, theta_star) triang = matplotlib.tri.Triangulation(X, Y) - cmap = plt.cm.get_cmap('Greys') + cmap = matplotlib.colormaps["Greys"] plt.tricontourf(triang, Z, cmap=cmap) except: - print('Objective contour plot for', xvar, yvar, 'slice failed') + print("Objective contour plot for", xvar, yvar, "slice failed") def _set_axis_limits(g, axis_limits, theta_vals, theta_star): @@ -277,7 +277,7 @@ def pairwise_plot( assert isinstance(theta_star, (type(None), dict, pd.Series, pd.DataFrame)) assert isinstance(alpha, (type(None), int, float)) assert isinstance(distributions, list) - assert set(distributions).issubset(set(['MVN', 'KDE', 'Rect'])) + assert set(distributions).issubset(set(["MVN", "KDE", "Rect"])) assert isinstance(axis_limits, (type(None), dict)) assert isinstance(title, (type(None), str)) assert isinstance(add_obj_contour, bool) @@ -307,7 +307,7 @@ def pairwise_plot( theta_names = [ col for col in theta_values.columns - if (col not in ['obj']) + if (col not in ["obj"]) and (not isinstance(col, float)) and (not isinstance(col, int)) ] @@ -335,7 +335,7 @@ def pairwise_plot( g.map_diag(sns.distplot, kde=False, hist=True, norm_hist=False) # Plot filled contours using all theta values based on obj - if 'obj' in theta_values.columns and add_obj_contour: + if "obj" in theta_values.columns and add_obj_contour: g.map_offdiag( _add_obj_contour, columns=theta_names, @@ -349,10 +349,10 @@ def pairwise_plot( matplotlib.lines.Line2D( [0], [0], - marker='o', - color='w', - label='thetas', - markerfacecolor='cadetblue', + marker="o", + color="w", + label="thetas", + markerfacecolor="cadetblue", markersize=5, ) ) @@ -360,23 +360,23 @@ def pairwise_plot( # Plot theta* if theta_star is not None: g.map_offdiag( - _add_scatter, color='k', columns=theta_names, theta_star=theta_star + _add_scatter, color="k", columns=theta_names, theta_star=theta_star ) legend_elements.append( matplotlib.lines.Line2D( [0], [0], - marker='o', - color='w', - label='theta*', - markerfacecolor='k', + marker="o", + color="w", + label="theta*", + markerfacecolor="k", markersize=6, ) ) # Plot confidence regions - colors = ['r', 'mediumblue', 'darkgray'] + colors = ["r", "mediumblue", "darkgray"] if (alpha is not None) and (len(distributions) > 0): if theta_star is None: print( @@ -388,7 +388,7 @@ def pairwise_plot( mvn_dist = None kde_dist = None for i, dist in enumerate(distributions): - if dist == 'Rect': + if dist == "Rect": lb, ub = fit_rect_dist(thetas, alpha) g.map_offdiag( _add_rectangle_CI, @@ -401,7 +401,7 @@ def pairwise_plot( matplotlib.lines.Line2D([0], [0], color=colors[i], lw=1, label=dist) ) - elif dist == 'MVN': + elif dist == "MVN": mvn_dist = fit_mvn_dist(thetas) Z = mvn_dist.pdf(thetas) score = stats.scoreatpercentile(Z, (1 - alpha) * 100) @@ -418,7 +418,7 @@ def pairwise_plot( matplotlib.lines.Line2D([0], [0], color=colors[i], lw=1, label=dist) ) - elif dist == 'KDE': + elif dist == "KDE": kde_dist = fit_kde_dist(thetas) Z = kde_dist.pdf(thetas.transpose()) score = stats.scoreatpercentile(Z, (1 - alpha) * 100) @@ -438,12 +438,12 @@ def pairwise_plot( _set_axis_limits(g, axis_limits, thetas, theta_star) for ax in g.axes.flatten(): - ax.ticklabel_format(style='sci', scilimits=(-2, 2), axis='both') + ax.ticklabel_format(style="sci", scilimits=(-2, 2), axis="both") if add_legend: xvar, yvar, loc = _get_variables(ax, theta_names) if loc == (len(theta_names) - 1, 0): - ax.legend(handles=legend_elements, loc='best', prop={'size': 8}) + ax.legend(handles=legend_elements, loc="best", prop={"size": 8}) if title: g.fig.subplots_adjust(top=0.9) g.fig.suptitle(title) @@ -474,7 +474,7 @@ def pairwise_plot( ax.tick_params(reset=True) if add_legend: - ax.legend(handles=legend_elements, loc='best', prop={'size': 8}) + ax.legend(handles=legend_elements, loc="best", prop={"size": 8}) plt.close(g.fig) @@ -563,15 +563,15 @@ def _get_grouped_data(data1, data2, normalize, group_names): # Combine data1 and data2 to create a grouped histogram data = pd.concat({group_names[0]: data1, group_names[1]: data2}) data.reset_index(level=0, inplace=True) - data.rename(columns={'level_0': 'set'}, inplace=True) + data.rename(columns={"level_0": "set"}, inplace=True) - data = data.melt(id_vars='set', value_vars=data1.columns, var_name='columns') + data = data.melt(id_vars="set", value_vars=data1.columns, var_name="columns") return data def grouped_boxplot( - data1, data2, normalize=False, group_names=['data1', 'data2'], filename=None + data1, data2, normalize=False, group_names=["data1", "data2"], filename=None ): """ Plot a grouped boxplot to compare two datasets @@ -600,11 +600,11 @@ def grouped_boxplot( data = _get_grouped_data(data1, data2, normalize, group_names) plt.figure() - sns.boxplot(data=data, hue='set', y='value', x='columns', order=data1.columns) + sns.boxplot(data=data, hue="set", y="value", x="columns", order=data1.columns) - plt.gca().legend().set_title('') - plt.gca().set_xlabel('') - plt.gca().set_ylabel('') + plt.gca().legend().set_title("") + plt.gca().set_xlabel("") + plt.gca().set_ylabel("") if filename is None: plt.show() @@ -614,7 +614,7 @@ def grouped_boxplot( def grouped_violinplot( - data1, data2, normalize=False, group_names=['data1', 'data2'], filename=None + data1, data2, normalize=False, group_names=["data1", "data2"], filename=None ): """ Plot a grouped violinplot to compare two datasets @@ -644,12 +644,12 @@ def grouped_violinplot( plt.figure() sns.violinplot( - data=data, hue='set', y='value', x='columns', order=data1.columns, split=True + data=data, hue="set", y="value", x="columns", order=data1.columns, split=True ) - plt.gca().legend().set_title('') - plt.gca().set_xlabel('') - plt.gca().set_ylabel('') + plt.gca().legend().set_title("") + plt.gca().set_xlabel("") + plt.gca().set_ylabel("") if filename is None: plt.show() diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index f26ecec2fce..2cc8ad36b0a 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -22,7 +22,7 @@ import platform -is_osx = platform.mac_ver()[0] != '' +is_osx = platform.mac_ver()[0] != "" import pyomo.common.unittest as unittest import sys @@ -38,11 +38,11 @@ from pyomo.opt import SolverFactory -ipopt_available = SolverFactory('ipopt').available() +ipopt_available = SolverFactory("ipopt").available() from pyomo.common.fileutils import find_library -pynumero_ASL_available = False if find_library('pynumero_ASL') is None else True +pynumero_ASL_available = False if find_library("pynumero_ASL") is None else True testdir = os.path.dirname(os.path.abspath(__file__)) @@ -61,10 +61,10 @@ def setUp(self): # 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], [7, 19.8]], - columns=['hour', 'y'], + columns=["hour", "y"], ) - theta_names = ['asymptote', 'rate_constant'] + theta_names = ["asymptote", "rate_constant"] def SSE(model, data): expr = sum( @@ -73,7 +73,7 @@ def SSE(model, data): ) return expr - solver_options = {'tol': 1e-8} + solver_options = {"tol": 1e-8} self.data = data self.pest = parmest.Estimator( @@ -90,10 +90,10 @@ def test_theta_est(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper @unittest.skipIf( @@ -105,14 +105,14 @@ def test_bootstrap(self): num_bootstraps = 10 theta_est = self.pest.theta_est_bootstrap(num_bootstraps, return_samples=True) - num_samples = theta_est['samples'].apply(len) + num_samples = theta_est["samples"].apply(len) self.assertTrue(len(theta_est.index), 10) self.assertTrue(num_samples.equals(pd.Series([6] * 10))) - del theta_est['samples'] + del theta_est["samples"] # apply confidence region test - CR = self.pest.confidence_region_test(theta_est, 'MVN', [0.5, 0.75, 1.0]) + CR = self.pest.confidence_region_test(theta_est, "MVN", [0.5, 0.75, 1.0]) self.assertTrue(set(CR.columns) >= set([0.5, 0.75, 1.0])) self.assertTrue(CR[0.5].sum() == 5) @@ -121,7 +121,7 @@ def test_bootstrap(self): graphics.pairwise_plot(theta_est) graphics.pairwise_plot(theta_est, thetavals) - graphics.pairwise_plot(theta_est, thetavals, 0.8, ['MVN', 'KDE', 'Rect']) + graphics.pairwise_plot(theta_est, thetavals, 0.8, ["MVN", "KDE", "Rect"]) @unittest.skipIf( not graphics.imports_available, "parmest.graphics imports are unavailable" @@ -151,7 +151,7 @@ def test_leaveNout(self): self.assertTrue(lNo_theta.shape == (6, 2)) results = self.pest.leaveNout_bootstrap_test( - 1, None, 3, 'Rect', [0.5, 1.0], seed=5436 + 1, None, 3, "Rect", [0.5, 1.0], seed=5436 ) self.assertTrue(len(results) == 6) # 6 lNo samples i = 1 @@ -221,10 +221,10 @@ def test_theta_est_cov(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper # Covariance matrix @@ -239,22 +239,22 @@ def test_theta_est_cov(self): ) # -0.4322 from paper self.assertAlmostEqual(cov.iloc[1, 1], 0.04124, places=2) # 0.04124 from paper - ''' Why does the covariance matrix from parmest not match the paper? Parmest is + """ 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 the paper used an alternative finite difference approximation for the Jacobian. - ''' + """ def model(theta, t): - ''' + """ Model to be fitted y = model(theta, t) Arguments: theta: vector of fitted parameters @@ -262,32 +262,32 @@ def model(theta, t): Returns: y: model predictions [need to check paper for units] - ''' + """ asymptote = theta[0] rate_constant = theta[1] return asymptote * (1 - np.exp(-rate_constant * t)) def residual(theta, t, y): - ''' + """ Calculate residuals Arguments: theta: vector of fitted parameters t: independent variable [hours] y: dependent variable [?] - ''' + """ return y - model(theta, t) # define data - t = self.data['hour'].to_numpy() - y = self.data['y'].to_numpy() + t = self.data["hour"].to_numpy() + y = self.data["y"].to_numpy() # define initial guess theta_guess = np.array([15, 0.5]) ## solve with optimize.least_squares sol = scipy.optimize.least_squares( - residual, theta_guess, method='trf', args=(t, y), verbose=2 + residual, theta_guess, method="trf", args=(t, y), verbose=2 ) theta_hat = sol.x @@ -313,18 +313,18 @@ def residual(theta, t, y): self.assertAlmostEqual(cov[1, 1], 0.04124, places=2) # 0.04124 from paper def test_cov_scipy_curve_fit_comparison(self): - ''' + """ Scipy results differ in the 3rd decimal place from the paper. It is possible the paper used an alternative finite difference approximation for the Jacobian. - ''' + """ ## solve with optimize.curve_fit def model(t, asymptote, rate_constant): return asymptote * (1 - np.exp(-rate_constant * t)) # define data - t = self.data['hour'].to_numpy() - y = self.data['y'].to_numpy() + t = self.data["hour"].to_numpy() + y = self.data["y"].to_numpy() # define initial guess theta_guess = np.array([15, 0.5]) @@ -351,7 +351,7 @@ class TestModelVariants(unittest.TestCase): def setUp(self): self.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'], + columns=["hour", "y"], ) def rooney_biegler_params(data): @@ -371,16 +371,16 @@ def response_rule(m, h): def rooney_biegler_indexed_params(data): model = pyo.ConcreteModel() - model.param_names = pyo.Set(initialize=['asymptote', 'rate_constant']) + model.param_names = pyo.Set(initialize=["asymptote", "rate_constant"]) model.theta = pyo.Param( model.param_names, - initialize={'asymptote': 15, 'rate_constant': 0.5}, + initialize={"asymptote": 15, "rate_constant": 0.5}, mutable=True, ) def response_rule(m, h): - expr = m.theta['asymptote'] * ( - 1 - pyo.exp(-m.theta['rate_constant'] * h) + expr = m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * h) ) return expr @@ -407,20 +407,20 @@ def response_rule(m, h): def rooney_biegler_indexed_vars(data): model = pyo.ConcreteModel() - model.var_names = pyo.Set(initialize=['asymptote', 'rate_constant']) + model.var_names = pyo.Set(initialize=["asymptote", "rate_constant"]) model.theta = pyo.Var( - model.var_names, initialize={'asymptote': 15, 'rate_constant': 0.5} + model.var_names, initialize={"asymptote": 15, "rate_constant": 0.5} ) model.theta[ - 'asymptote' + "asymptote" ].fixed = ( True # parmest will unfix theta variables, even when they are indexed ) - model.theta['rate_constant'].fixed = True + model.theta["rate_constant"].fixed = True def response_rule(m, h): - expr = m.theta['asymptote'] * ( - 1 - pyo.exp(-m.theta['rate_constant'] * h) + expr = m.theta["asymptote"] * ( + 1 - pyo.exp(-m.theta["rate_constant"] * h) ) return expr @@ -437,41 +437,41 @@ def SSE(model, data): self.objective_function = SSE - theta_vals = pd.DataFrame([20, 1], index=['asymptote', 'rate_constant']).T + theta_vals = pd.DataFrame([20, 1], index=["asymptote", "rate_constant"]).T theta_vals_index = pd.DataFrame( [20, 1], index=["theta['asymptote']", "theta['rate_constant']"] ).T self.input = { - 'param': { - 'model': rooney_biegler_params, - 'theta_names': ['asymptote', 'rate_constant'], - 'theta_vals': theta_vals, + "param": { + "model": rooney_biegler_params, + "theta_names": ["asymptote", "rate_constant"], + "theta_vals": theta_vals, }, - 'param_index': { - 'model': rooney_biegler_indexed_params, - 'theta_names': ['theta'], - 'theta_vals': theta_vals_index, + "param_index": { + "model": rooney_biegler_indexed_params, + "theta_names": ["theta"], + "theta_vals": theta_vals_index, }, - 'vars': { - 'model': rooney_biegler_vars, - 'theta_names': ['asymptote', 'rate_constant'], - 'theta_vals': theta_vals, + "vars": { + "model": rooney_biegler_vars, + "theta_names": ["asymptote", "rate_constant"], + "theta_vals": theta_vals, }, - 'vars_index': { - 'model': rooney_biegler_indexed_vars, - 'theta_names': ['theta'], - 'theta_vals': theta_vals_index, + "vars_index": { + "model": rooney_biegler_indexed_vars, + "theta_names": ["theta"], + "theta_vals": theta_vals_index, }, - 'vars_quoted_index': { - 'model': rooney_biegler_indexed_vars, - 'theta_names': ["theta['asymptote']", "theta['rate_constant']"], - 'theta_vals': theta_vals_index, + "vars_quoted_index": { + "model": rooney_biegler_indexed_vars, + "theta_names": ["theta['asymptote']", "theta['rate_constant']"], + "theta_vals": theta_vals_index, }, - 'vars_str_index': { - 'model': rooney_biegler_indexed_vars, - 'theta_names': ["theta[asymptote]", "theta[rate_constant]"], - 'theta_vals': theta_vals_index, + "vars_str_index": { + "model": rooney_biegler_indexed_vars, + "theta_names": ["theta[asymptote]", "theta[rate_constant]"], + "theta_vals": theta_vals_index, }, } @@ -483,9 +483,9 @@ def SSE(model, data): def test_parmest_basics(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) @@ -505,15 +505,15 @@ def test_parmest_basics(self): cov.iloc[1, 1], 0.04193591, places=2 ) # 0.04124 from paper - obj_at_theta = pest.objective_at_theta(parmest_input['theta_vals']) - self.assertAlmostEqual(obj_at_theta['obj'][0], 16.531953, places=2) + obj_at_theta = pest.objective_at_theta(parmest_input["theta_vals"]) + self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) def test_parmest_basics_with_initialize_parmest_model_option(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) @@ -534,22 +534,22 @@ def test_parmest_basics_with_initialize_parmest_model_option(self): ) # 0.04124 from paper obj_at_theta = pest.objective_at_theta( - parmest_input['theta_vals'], initialize_parmest_model=True + parmest_input["theta_vals"], initialize_parmest_model=True ) - self.assertAlmostEqual(obj_at_theta['obj'][0], 16.531953, places=2) + self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) def test_parmest_basics_with_square_problem_solve(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) obj_at_theta = pest.objective_at_theta( - parmest_input['theta_vals'], initialize_parmest_model=True + parmest_input["theta_vals"], initialize_parmest_model=True ) objval, thetavals, cov = pest.theta_est(calc_cov=True, cov_n=6) @@ -568,14 +568,14 @@ def test_parmest_basics_with_square_problem_solve(self): cov.iloc[1, 1], 0.04193591, places=2 ) # 0.04124 from paper - self.assertAlmostEqual(obj_at_theta['obj'][0], 16.531953, places=2) + self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2) def test_parmest_basics_with_square_problem_solve_no_theta_vals(self): for model_type, parmest_input in self.input.items(): pest = parmest.Estimator( - parmest_input['model'], + parmest_input["model"], self.data, - parmest_input['theta_names'], + parmest_input["theta_names"], self.objective_function, ) @@ -632,17 +632,17 @@ def setUp(self): [1.90, 10000, 4491.3, 1049.4, 920.5, 1769.4], [1.95, 10000, 4538.8, 1045.8, 893.9, 1760.8], ], - columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'], + columns=["sv", "caf", "ca", "cb", "cc", "cd"], ) - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] 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 @@ -656,13 +656,13 @@ def test_theta_est(self): # used in data reconciliation objval, thetavals = self.pest.theta_est() - self.assertAlmostEqual(thetavals['k1'], 5.0 / 6.0, places=4) - self.assertAlmostEqual(thetavals['k2'], 5.0 / 3.0, places=4) - self.assertAlmostEqual(thetavals['k3'], 1.0 / 6000.0, places=7) + self.assertAlmostEqual(thetavals["k1"], 5.0 / 6.0, places=4) + self.assertAlmostEqual(thetavals["k2"], 5.0 / 3.0, places=4) + self.assertAlmostEqual(thetavals["k3"], 1.0 / 6000.0, places=7) def test_return_values(self): objval, thetavals, data_rec = self.pest.theta_est( - return_values=['ca', 'cb', 'cc', 'cd', 'caf'] + return_values=["ca", "cb", "cc", "cd", "caf"] ) self.assertAlmostEqual(data_rec["cc"].loc[18], 893.84924, places=3) @@ -679,9 +679,9 @@ class TestReactorDesign_DAE(unittest.TestCase): def setUp(self): def ABC_model(data): - ca_meas = data['ca'] - cb_meas = data['cb'] - cc_meas = data['cc'] + ca_meas = data["ca"] + cb_meas = data["cb"] + cc_meas = data["cc"] if isinstance(data, pd.DataFrame): meas_t = data.index # time index @@ -754,7 +754,7 @@ def total_cost_rule(model): rule=total_cost_rule, sense=pyo.minimize ) - disc = pyo.TransformationFactory('dae.collocation') + disc = pyo.TransformationFactory("dae.collocation") disc.apply_to(m, nfe=20, ncp=2) return m @@ -785,15 +785,15 @@ def total_cost_rule(model): [4.737, 0.004, 0.036, 0.971], [5.000, -0.024, 0.028, 0.985], ] - data = pd.DataFrame(data, columns=['t', 'ca', 'cb', 'cc']) - data_df = data.set_index('t') + data = pd.DataFrame(data, columns=["t", "ca", "cb", "cc"]) + data_df = data.set_index("t") data_dict = { - 'ca': {k: v for (k, v) in zip(data.t, data.ca)}, - 'cb': {k: v for (k, v) in zip(data.t, data.cb)}, - 'cc': {k: v for (k, v) in zip(data.t, data.cc)}, + "ca": {k: v for (k, v) in zip(data.t, data.ca)}, + "cb": {k: v for (k, v) in zip(data.t, data.cb)}, + "cc": {k: v for (k, v) in zip(data.t, data.cc)}, } - theta_names = ['k1', 'k2'] + theta_names = ["k1", "k2"] self.pest_df = parmest.Estimator(ABC_model, [data_df], theta_names) self.pest_dict = parmest.Estimator(ABC_model, [data_dict], theta_names) @@ -815,30 +815,30 @@ def test_dataformats(self): obj2, theta2 = self.pest_dict.theta_est() self.assertAlmostEqual(obj1, obj2, places=6) - self.assertAlmostEqual(theta1['k1'], theta2['k1'], places=6) - self.assertAlmostEqual(theta1['k2'], theta2['k2'], places=6) + self.assertAlmostEqual(theta1["k1"], theta2["k1"], places=6) + self.assertAlmostEqual(theta1["k2"], theta2["k2"], places=6) def test_return_continuous_set(self): - ''' + """ test if ContinuousSet elements are returned correctly from theta_est() - ''' - obj1, theta1, return_vals1 = self.pest_df.theta_est(return_values=['time']) - obj2, theta2, return_vals2 = self.pest_dict.theta_est(return_values=['time']) - self.assertAlmostEqual(return_vals1['time'].loc[0][18], 2.368, places=3) - self.assertAlmostEqual(return_vals2['time'].loc[0][18], 2.368, places=3) + """ + obj1, theta1, return_vals1 = self.pest_df.theta_est(return_values=["time"]) + obj2, theta2, return_vals2 = self.pest_dict.theta_est(return_values=["time"]) + self.assertAlmostEqual(return_vals1["time"].loc[0][18], 2.368, places=3) + self.assertAlmostEqual(return_vals2["time"].loc[0][18], 2.368, places=3) def test_return_continuous_set_multiple_datasets(self): - ''' + """ test if ContinuousSet elements are returned correctly from theta_est() - ''' + """ obj1, theta1, return_vals1 = self.pest_df_multiple.theta_est( - return_values=['time'] + return_values=["time"] ) obj2, theta2, return_vals2 = self.pest_dict_multiple.theta_est( - return_values=['time'] + return_values=["time"] ) - self.assertAlmostEqual(return_vals1['time'].loc[1][18], 2.368, places=3) - self.assertAlmostEqual(return_vals2['time'].loc[1][18], 2.368, places=3) + self.assertAlmostEqual(return_vals1["time"].loc[1][18], 2.368, places=3) + self.assertAlmostEqual(return_vals2["time"].loc[1][18], 2.368, places=3) def test_covariance(self): from pyomo.contrib.interior_point.inverse_reduced_hessian import ( @@ -862,13 +862,13 @@ def test_covariance(self): l = len(vars_list) cov_interior_point = 2 * obj / (n - l) * inv_red_hes cov_interior_point = pd.DataFrame( - cov_interior_point, ['k1', 'k2'], ['k1', 'k2'] + cov_interior_point, ["k1", "k2"], ["k1", "k2"] ) cov_diff = (cov - cov_interior_point).abs().sum().sum() - self.assertTrue(cov.loc['k1', 'k1'] > 0) - self.assertTrue(cov.loc['k2', 'k2'] > 0) + self.assertTrue(cov.loc["k1", "k1"] > 0) + self.assertTrue(cov.loc["k2", "k2"] > 0) self.assertAlmostEqual(cov_diff, 0, places=6) @@ -886,10 +886,10 @@ def setUp(self): # 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], [7, 19.8]], - columns=['hour', 'y'], + columns=["hour", "y"], ) - theta_names = ['asymptote', 'rate_constant'] + theta_names = ["asymptote", "rate_constant"] def SSE(model, data): expr = sum( @@ -898,7 +898,7 @@ def SSE(model, data): ) return expr - solver_options = {'tol': 1e-8} + solver_options = {"tol": 1e-8} self.data = data self.pest = parmest.Estimator( @@ -916,15 +916,15 @@ def test_theta_est_with_square_initialization(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper def test_theta_est_with_square_initialization_and_custom_init_theta(self): theta_vals_init = pd.DataFrame( - data=[[19.0, 0.5]], columns=['asymptote', 'rate_constant'] + data=[[19.0, 0.5]], columns=["asymptote", "rate_constant"] ) obj_init = self.pest.objective_at_theta( theta_values=theta_vals_init, initialize_parmest_model=True @@ -932,10 +932,10 @@ def test_theta_est_with_square_initialization_and_custom_init_theta(self): objval, thetavals = self.pest.theta_est() self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper def test_theta_est_with_square_initialization_diagnostic_mode_true(self): @@ -945,14 +945,14 @@ def test_theta_est_with_square_initialization_diagnostic_mode_true(self): self.assertAlmostEqual(objval, 4.3317112, places=2) self.assertAlmostEqual( - thetavals['asymptote'], 19.1426, places=2 + thetavals["asymptote"], 19.1426, places=2 ) # 19.1426 from the paper self.assertAlmostEqual( - thetavals['rate_constant'], 0.5311, places=2 + thetavals["rate_constant"], 0.5311, places=2 ) # 0.5311 from the paper self.pest.diagnostic_mode = False -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/parmest/tests/test_scenariocreator.py b/pyomo/contrib/parmest/tests/test_scenariocreator.py index a2dcf4c2739..22a851ae32e 100644 --- a/pyomo/contrib/parmest/tests/test_scenariocreator.py +++ b/pyomo/contrib/parmest/tests/test_scenariocreator.py @@ -24,7 +24,7 @@ import pyomo.environ as pyo from pyomo.environ import SolverFactory -ipopt_available = SolverFactory('ipopt').available() +ipopt_available = SolverFactory("ipopt").available() testdir = os.path.dirname(os.path.abspath(__file__)) @@ -63,17 +63,17 @@ def setUp(self): [1.90, 10000, 4491.3, 1049.4, 920.5, 1769.4], [1.95, 10000, 4538.8, 1045.8, 893.9, 1760.8], ], - columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'], + columns=["sv", "caf", "ca", "cb", "cc", "cd"], ) - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] 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 @@ -116,7 +116,7 @@ def setUp(self): import json # Vars to estimate in parmest - theta_names = ['k1', 'k2', 'E1', 'E2'] + theta_names = ["k1", "k2", "E1", "E2"] self.fbase = os.path.join(testdir, "..", "examples", "semibatch") # Data, list of dictionaries @@ -124,7 +124,7 @@ def setUp(self): for exp_num in range(10): fname = "exp" + str(exp_num + 1) + ".out" fullname = os.path.join(self.fbase, fname) - with open(fullname, 'r') as infile: + with open(fullname, "r") as infile: d = json.load(infile) data.append(d) @@ -142,5 +142,5 @@ def test_semibatch_bootstrap(self): self.assertAlmostEqual(tval, 20.64, places=1) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/pyomo/contrib/parmest/tests/test_utils.py b/pyomo/contrib/parmest/tests/test_utils.py index bd0706ac38d..514c14b1e82 100644 --- a/pyomo/contrib/parmest/tests/test_utils.py +++ b/pyomo/contrib/parmest/tests/test_utils.py @@ -16,7 +16,7 @@ import pyomo.contrib.parmest.parmest as parmest from pyomo.opt import SolverFactory -ipopt_available = SolverFactory('ipopt').available() +ipopt_available = SolverFactory("ipopt").available() @unittest.skipIf( @@ -45,13 +45,13 @@ def test_convert_param_to_var(self): [1.10, 10000, 3535.1, 1064.8, 1613.3, 1893.4], [1.15, 10000, 3609.1, 1067.8, 1547.5, 1887.8], ], - columns=['sv', 'caf', 'ca', 'cb', 'cc', 'cd'], + columns=["sv", "caf", "ca", "cb", "cc", "cd"], ) - theta_names = ['k1', 'k2', 'k3'] + theta_names = ["k1", "k2", "k3"] instance = reactor_design_model(data.loc[0]) - solver = pyo.SolverFactory('ipopt') + solver = pyo.SolverFactory("ipopt") solver.solve(instance) instance_vars = parmest.utils.convert_params_to_vars(