From f6e211e67f9c8d0e4b2fdeddb0357da750da8f9e Mon Sep 17 00:00:00 2001 From: Shawn Martin Date: Tue, 31 Oct 2023 09:18:32 -0600 Subject: [PATCH 1/4] Fixed Pandas warning in parmest examples plus bug in semibatch parallel example. --- .../reactor_design/bootstrap_example.py | 8 +- .../reactor_design/datarec_example.py | 8 +- .../reactor_design/leaveNout_example.py | 8 +- .../likelihood_ratio_example.py | 8 +- .../multisensor_data_example.py | 14 +- .../parameter_estimation_example.py | 8 +- .../examples/reactor_design/reactor_design.py | 12 +- .../examples/semibatch/obj_at_theta.csv | 1009 +++++++++++++++++ .../parmest/examples/semibatch/semibatch.py | 10 + 9 files changed, 1055 insertions(+), 30 deletions(-) create mode 100644 pyomo/contrib/parmest/examples/semibatch/obj_at_theta.csv diff --git a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py b/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py index cf1b8a2de23..67724644ef5 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py @@ -29,10 +29,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 diff --git a/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py b/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py index 811571e20ed..b50ee46d9b9 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/datarec_example.py @@ -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 diff --git a/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py b/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py index 95af53e63d3..1e14e1fb329 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/leaveNout_example.py @@ -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 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..5224097c13f 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/likelihood_ratio_example.py @@ -31,10 +31,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 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..af7620b47b3 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/multisensor_data_example.py @@ -31,13 +31,13 @@ def main(): # 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..070c5934be5 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py +++ b/pyomo/contrib/parmest/examples/reactor_design/parameter_estimation_example.py @@ -29,10 +29,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 diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py index f4cd6c8dbf5..3284a174e93 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py @@ -37,10 +37,16 @@ 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): + model.caf = Param(initialize=float(data['caf']), within=PositiveReals) + else: + model.caf = Param(initialize=float(data.iloc[0]['caf']), within=PositiveReals) # Space velocity (flowrate/volume) - model.sv = Param(initialize=float(data['sv']), within=PositiveReals) + if isinstance(data, dict): + model.sv = Param(initialize=float(data['sv']), within=PositiveReals) + else: + model.sv = Param(initialize=float(data.iloc[0]['sv']), within=PositiveReals) # Outlet concentration of each component model.ca = Var(initialize=5000.0, within=PositiveReals) @@ -81,7 +87,7 @@ 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}) + 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()]) 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..b3da21ed993 100644 --- a/pyomo/contrib/parmest/examples/semibatch/semibatch.py +++ b/pyomo/contrib/parmest/examples/semibatch/semibatch.py @@ -34,6 +34,16 @@ 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'] From 42d825ef86c74633702840e2159c49bf7531b2ab Mon Sep 17 00:00:00 2001 From: Shawn Martin Date: Tue, 31 Oct 2023 11:26:53 -0600 Subject: [PATCH 2/4] Fixed deprecations in parmest tests. --- .../examples/reactor_design/reactor_design.py | 14 +++++++++----- pyomo/contrib/parmest/graphics.py | 3 ++- pyomo/contrib/parmest/tests/test_parmest.py | 8 ++++---- .../contrib/parmest/tests/test_scenariocreator.py | 8 ++++---- pyomo/contrib/parmest/tests/test_utils.py | 1 + 5 files changed, 20 insertions(+), 14 deletions(-) diff --git a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py index 3284a174e93..80df6fb3c12 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py @@ -37,17 +37,21 @@ def reactor_design_model(data): ) # m^3/(gmol min) # Inlet concentration of A, gmol/m^3 - if isinstance(data, dict): + if isinstance(data, dict) or isinstance(data, pd.Series): model.caf = Param(initialize=float(data['caf']), within=PositiveReals) - else: + 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) - if isinstance(data, dict): + if isinstance(data, dict) or isinstance(data, pd.Series): model.sv = Param(initialize=float(data['sv']), within=PositiveReals) - else: + 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) model.cb = Var(initialize=2000.0, within=PositiveReals) diff --git a/pyomo/contrib/parmest/graphics.py b/pyomo/contrib/parmest/graphics.py index f01622d2d17..991395a556d 100644 --- a/pyomo/contrib/parmest/graphics.py +++ b/pyomo/contrib/parmest/graphics.py @@ -178,7 +178,8 @@ 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 = plt.cm.get_cmap('Greys') + cmap = matplotlib.colormaps['Greys'] plt.tricontourf(triang, Z, cmap=cmap) except: diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py index f26ecec2fce..b09c48e9709 100644 --- a/pyomo/contrib/parmest/tests/test_parmest.py +++ b/pyomo/contrib/parmest/tests/test_parmest.py @@ -639,10 +639,10 @@ def setUp(self): 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 diff --git a/pyomo/contrib/parmest/tests/test_scenariocreator.py b/pyomo/contrib/parmest/tests/test_scenariocreator.py index a2dcf4c2739..fe4528120f6 100644 --- a/pyomo/contrib/parmest/tests/test_scenariocreator.py +++ b/pyomo/contrib/parmest/tests/test_scenariocreator.py @@ -70,10 +70,10 @@ def setUp(self): 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 diff --git a/pyomo/contrib/parmest/tests/test_utils.py b/pyomo/contrib/parmest/tests/test_utils.py index bd0706ac38d..03b48c0326c 100644 --- a/pyomo/contrib/parmest/tests/test_utils.py +++ b/pyomo/contrib/parmest/tests/test_utils.py @@ -50,6 +50,7 @@ def test_convert_param_to_var(self): theta_names = ['k1', 'k2', 'k3'] + print(data.loc[0]) instance = reactor_design_model(data.loc[0]) solver = pyo.SolverFactory('ipopt') solver.solve(instance) From fe31c4049cd20af7a97ae3a81a431a3b47de4a23 Mon Sep 17 00:00:00 2001 From: Shawn Martin Date: Tue, 31 Oct 2023 15:37:44 -0600 Subject: [PATCH 3/4] Removed extraneous code. --- pyomo/contrib/parmest/graphics.py | 1 - pyomo/contrib/parmest/tests/test_utils.py | 1 - 2 files changed, 2 deletions(-) diff --git a/pyomo/contrib/parmest/graphics.py b/pyomo/contrib/parmest/graphics.py index 991395a556d..99eda8aad7a 100644 --- a/pyomo/contrib/parmest/graphics.py +++ b/pyomo/contrib/parmest/graphics.py @@ -178,7 +178,6 @@ 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) diff --git a/pyomo/contrib/parmest/tests/test_utils.py b/pyomo/contrib/parmest/tests/test_utils.py index 03b48c0326c..bd0706ac38d 100644 --- a/pyomo/contrib/parmest/tests/test_utils.py +++ b/pyomo/contrib/parmest/tests/test_utils.py @@ -50,7 +50,6 @@ def test_convert_param_to_var(self): theta_names = ['k1', 'k2', 'k3'] - print(data.loc[0]) instance = reactor_design_model(data.loc[0]) solver = pyo.SolverFactory('ipopt') solver.solve(instance) From a79912c0a0c4fecc59d1c40185d5d36619088f8b Mon Sep 17 00:00:00 2001 From: Shawn Martin Date: Thu, 2 Nov 2023 12:55:30 -0600 Subject: [PATCH 4/4] Ran black on modified files. --- .../reactor_design/bootstrap_example.py | 18 +- .../reactor_design/datarec_example.py | 34 +-- .../reactor_design/leaveNout_example.py | 18 +- .../likelihood_ratio_example.py | 16 +- .../multisensor_data_example.py | 18 +- .../parameter_estimation_example.py | 18 +- .../examples/reactor_design/reactor_design.py | 20 +- .../parmest/examples/semibatch/semibatch.py | 47 ++-- pyomo/contrib/parmest/graphics.py | 76 +++--- pyomo/contrib/parmest/tests/test_parmest.py | 256 +++++++++--------- .../parmest/tests/test_scenariocreator.py | 20 +- pyomo/contrib/parmest/tests/test_utils.py | 8 +- 12 files changed, 274 insertions(+), 275 deletions(-) diff --git a/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py b/pyomo/contrib/parmest/examples/reactor_design/bootstrap_example.py index 67724644ef5..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.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 + (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 b50ee46d9b9..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.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 + ((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 1e14e1fb329..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.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 + (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 5224097c13f..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.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 + (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 af7620b47b3..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.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 + ((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 070c5934be5..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.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 + (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 80df6fb3c12..16f65e236eb 100644 --- a/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py +++ b/pyomo/contrib/parmest/examples/reactor_design/reactor_design.py @@ -38,20 +38,20 @@ def reactor_design_model(data): # Inlet concentration of A, gmol/m^3 if isinstance(data, dict) or isinstance(data, pd.Series): - model.caf = Param(initialize=float(data['caf']), within=PositiveReals) + 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) + model.caf = Param(initialize=float(data.iloc[0]["caf"]), within=PositiveReals) else: - raise ValueError('Unrecognized data type.') + raise ValueError("Unrecognized data type.") # Space velocity (flowrate/volume) if isinstance(data, dict) or isinstance(data, pd.Series): - model.sv = Param(initialize=float(data['sv']), within=PositiveReals) + 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) + model.sv = Param(initialize=float(data.iloc[0]["sv"]), within=PositiveReals) else: - raise ValueError('Unrecognized data type.') - + raise ValueError("Unrecognized data type.") + # Outlet concentration of each component model.ca = Var(initialize=5000.0, within=PositiveReals) model.cb = Var(initialize=2000.0, within=PositiveReals) @@ -91,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(pd.DataFrame(data={'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/semibatch.py b/pyomo/contrib/parmest/examples/semibatch/semibatch.py index b3da21ed993..6762531a338 100644 --- a/pyomo/contrib/parmest/examples/semibatch/semibatch.py +++ b/pyomo/contrib/parmest/examples/semibatch/semibatch.py @@ -34,21 +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: + with open(file_name, "r") as infile: data = json.load(infile) except: - raise RuntimeError(f'Could not read {file_name} as json') + 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 = {} @@ -89,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 @@ -102,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 @@ -112,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 @@ -240,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 @@ -248,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( @@ -264,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 @@ -272,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 99eda8aad7a..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 = matplotlib.colormaps['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 b09c48e9709..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.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 + (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 fe4528120f6..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.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 + (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(