Merge pull request #3723 from djlaky/parmest-doe-bugfixes

Pyomo.DoE - Measurement Error Convention Correction

Author: mrmundt

pyomo/contrib/doe/doe.py

@@ -731,7 +731,7 @@ def _sequential_FIM(self, model=None):
         cov_y = np.zeros((len(model.measurement_error), len(model.measurement_error)))
         count = 0
         for k, v in model.measurement_error.items():
-            cov_y[count, count] = 1 / v
+            cov_y[count, count] = 1 / v**2
             count += 1
 
         # Compute and record FIM
@@ -822,7 +822,7 @@ def _kaug_FIM(self, model=None):
         cov_y = np.zeros((len(model.measurement_error), len(model.measurement_error)))
         count = 0
         for k, v in model.measurement_error.items():
-            cov_y[count, count] = 1 / v
+            cov_y[count, count] = 1 / v**2
             count += 1
 
         # TODO: need to add a covariance matrix for measurements (sigma inverse)
@@ -1055,6 +1055,7 @@ def fim_rule(m, p, q):
                         / m.scenario_blocks[0].measurement_error[
                             pyo.ComponentUID(n).find_component_on(m.scenario_blocks[0])
                         ]
+                        ** 2
                         * m.sensitivity_jacobian[n, p]
                         * m.sensitivity_jacobian[n, q]
                         for n in m.output_names
@@ -1779,7 +1780,7 @@ def compute_FIM_full_factorial(
 
             FIM = self._computed_FIM
 
-            det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
+            (det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt) = (
                 compute_FIM_metrics(FIM)
             )
 

pyomo/contrib/doe/tests/test_doe_solve.py

@@ -83,7 +83,9 @@ def get_FIM_Q_L(doe_obj=None):
         for i in model.output_names
         for j in model.parameter_names
     ]
-    sigma_inv = [1 / v for k, v in model.scenario_blocks[0].measurement_error.items()]
+    sigma_inv = [
+        1 / v**2 for k, v in model.scenario_blocks[0].measurement_error.items()
+    ]
     param_vals = np.array(
         [[v for k, v in model.scenario_blocks[0].unknown_parameters.items()]]
     )
@@ -469,24 +471,24 @@ def test_reactor_grid_search_bad_model(self):
 class TestDoe(unittest.TestCase):
     def test_doe_full_factorial(self):
         log10_D_opt_expected = [
-            3.7734377852467524,
-            5.137792359070963,
-            5.182167857710023,
-            6.546522431509408,
+            11.77343778527225,
+            13.137792359064383,
+            13.182167857699808,
+            14.54652243150573,
         ]
 
         log10_A_opt_expected = [
-            3.5935726800929695,
-            3.6133186151486933,
-            3.945755198204365,
-            3.9655011332598367,
+            5.59357268009304,
+            5.613318615148643,
+            5.945755198204368,
+            5.965501133259909,
         ]
 
         log10_E_opt_expected = [
-            -1.7201873126109162,
-            -0.691340497355524,
-            -1.3680047944877138,
-            -0.3391579792516522,
+            0.27981268741620413,
+            1.3086595026369012,
+            0.6319952055040333,
+            1.6608420207466377,
         ]
 
         log10_ME_opt_expected = [
@@ -497,31 +499,31 @@ def test_doe_full_factorial(self):
         ]
 
         eigval_min_expected = [
-            0.019046390638130666,
-            0.20354456134677426,
-            0.04285437893696232,
-            0.45797526302234304,
+            1.9046390638130666,
+            20.354456134677426,
+            4.285437893696232,
+            45.797526302234304,
         ]
 
         eigval_max_expected = [
-            3169.552855492114,
-            3576.0292523637977,
-            7131.493924857995,
-            8046.0658178139165,
+            316955.2855492114,
+            357602.92523637977,
+            713149.3924857995,
+            804606.58178139165,
         ]
 
         det_FIM_expected = [
-            5935.233170586055,
-            137338.51875774842,
-            152113.5345070818,
-            3519836.021699428,
+            593523317093.4525,
+            13733851875566.766,
+            15211353450350.424,
+            351983602166961.56,
         ]
 
         trace_FIM_expected = [
-            3922.5878617108597,
-            4105.051549241871,
-            8825.822688850109,
-            9236.36598578955,
+            392258.78617108597,
+            410505.1549241871,
+            882582.2688850109,
+            923636.598578955,
         ]
         ff = run_reactor_doe(
             n_points_for_design=2,
@@ -551,15 +553,11 @@ def test_doe_full_factorial(self):
         self.assertStructuredAlmostEqual(
             ff_results["eigval_min"], eigval_min_expected, abstol=1e-4
         )
-        self.assertStructuredAlmostEqual(
-            ff_results["eigval_max"], eigval_max_expected, abstol=1e-4
-        )
-        self.assertStructuredAlmostEqual(
-            ff_results["det_FIM"], det_FIM_expected, abstol=1e-4
-        )
-        self.assertStructuredAlmostEqual(
-            ff_results["trace_FIM"], trace_FIM_expected, abstol=1e-4
-        )
+        # abstol of 1e-4 removed for the following values as
+        # their non-log values are large (e.g., >1e10)
+        self.assertStructuredAlmostEqual(ff_results["eigval_max"], eigval_max_expected)
+        self.assertStructuredAlmostEqual(ff_results["det_FIM"], det_FIM_expected)
+        self.assertStructuredAlmostEqual(ff_results["trace_FIM"], trace_FIM_expected)
 
 
 if __name__ == "__main__":

pyomo/contrib/doe/tests/test_utils.py

@@ -68,7 +68,7 @@ def test_compute_FIM_metrics(self):
         # Create a sample Fisher Information Matrix (FIM)
         FIM = np.array([[10, 2], [2, 3]])
 
-        det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
+        (det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt) = (
             compute_FIM_metrics(FIM)
         )
 

pyomo/contrib/doe/utils.py

@@ -219,7 +219,7 @@ def get_FIM_metrics(FIM):
         log10(Modified E-optimality) metric
     """
 
-    det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt = (
+    (det_FIM, trace_FIM, E_vals, E_vecs, D_opt, A_opt, E_opt, ME_opt) = (
         compute_FIM_metrics(FIM)
     )