add tests for uncertainty covariance matrices

Author: sevyharris

rmgpy/tools/uncertainty.py

@@ -944,7 +944,8 @@ def assign_intermediate_uncertainties(self, g_param_engine=None, k_param_engine=
 
                     for ruleEntry, trainingEntry, weight in training:
                         # TODO - test that training reactions in a tree are correlated with the exact match kind of training reaction
-                        label = '{}Training {} {}'.format(surface_prefix, family, ruleEntry)  # ruleEntry should probably be the reaction equation itself
+                        # for now, we follow the old convention of treating these as rate rules
+                        label = '{}Rate Rule {} {}'.format(surface_prefix, family, ruleEntry)  # ruleEntry should probably be the reaction equation itself
                         dlnkdq[label] = weight  # dlnk/dlnk_training = weight, because the training entry is scaled by the weight when it is used in the kinetics estimation
 
                     # There is also estimation error if rate rules are used

test/rmgpy/tools/uncertaintyTest.py

@@ -28,7 +28,7 @@
 ###############################################################################
 
 import os
-
+import copy
 
 import numpy as np
 
@@ -302,43 +302,108 @@ def test_local_analysis(self):
 
         # -------------------- repeat the exact same test for new formulation --------------------------
         # uncorrelated analysis first
-        self.uncertainty.assign_intermediate_uncertainties(correlated=False)
-        output = self.uncertainty.local_analysis_intermediate(sensitive_species=sensitive_species, correlated=False)
+        self.uncertainty.assign_intermediate_uncertainties()
+        output = self.uncertainty.local_analysis_intermediate(sensitive_species=sensitive_species)
         total_variance, kinetic_uncertainty, thermo_uncertainty = output[sensitive_species[0]]
         assert np.isclose(total_variance, expected_uncorrelated_total_variance)
 
-        # order of kinetic or thermo uncertainty is not guaranteed, this sorts by name
+        # order of kinetic or thermo uncertainty is not guaranteed, this sorts by contribution
         kinetic_variances = [r[2] for r in kinetic_uncertainty]
         kinetics_names = [r[0] for r in kinetic_uncertainty]
-        sorted_kinetic_variances = [x for _, x in sorted(zip(kinetics_names, kinetic_variances))]
+        sorted_kinetics_names = [x for _, x in sorted(zip(kinetic_variances, kinetics_names))][::-1]
+        sorted_kinetic_variances = sorted(kinetic_variances, reverse=True)
         assert np.isclose(sorted_kinetic_variances, expected_uncorrelated_kinetics_variances).all()
-        assert set(kinetics_names) == expected_uncorrelated_kinetics_labels
+        assert sorted_kinetics_names == expected_uncorrelated_kinetics_labels
 
         thermo_variances = [s[2] for s in thermo_uncertainty]
         thermo_names = [s[0] for s in thermo_uncertainty]
-        sorted_thermo_variances = [x for _, x in sorted(zip(thermo_names, thermo_variances))]
+        sorted_thermo_names = [x for _, x in sorted(zip(thermo_variances, thermo_names))][::-1]
+        sorted_thermo_variances = sorted(thermo_variances, reverse=True)
         assert np.isclose(sorted_thermo_variances, expected_uncorrelated_thermo_variances).all()
-        assert set(thermo_names) == expected_uncorrelated_thermo_labels
+        assert sorted_thermo_names == expected_uncorrelated_thermo_labels
 
         # now repeat for correlated analysis
         self.uncertainty.assign_intermediate_uncertainties(correlated=True)
         output = self.uncertainty.local_analysis_intermediate(sensitive_species=sensitive_species, correlated=True)
         total_variance, kinetic_uncertainty, thermo_uncertainty = output[sensitive_species[0]]
         assert np.isclose(total_variance, expected_correlated_total_variance)
 
-        # order of kinetic or thermo uncertainty is not guaranteed, this sorts by name
+        # order of kinetic or thermo uncertainty is not guaranteed, this sorts by contribution
         kinetic_variances = [r[2] for r in kinetic_uncertainty]
         kinetics_names = [r[0] for r in kinetic_uncertainty]
-        sorted_kinetic_variances = [x for _, x in sorted(zip(kinetics_names, kinetic_variances))]
+        sorted_kinetic_variances = sorted(kinetic_variances, reverse=True)
+        sorted_kinetics_names = [x for _, x in sorted(zip(kinetic_variances, kinetics_names))][::-1]
         assert np.isclose(sorted_kinetic_variances, expected_correlated_kinetics_variances).all()
-        assert set(kinetics_names) == expected_correlated_kinetics_labels
+        assert sorted_kinetics_names == expected_correlated_kinetics_labels
 
         thermo_variances = [s[2] for s in thermo_uncertainty]
         thermo_names = [s[0] for s in thermo_uncertainty]
-        sorted_thermo_variances = [x for _, x in sorted(zip(thermo_names, thermo_variances))]
+        sorted_thermo_variances = sorted(thermo_variances, reverse=True)
+        sorted_thermo_names = [x for _, x in sorted(zip(thermo_variances, thermo_names))][::-1]
         assert np.isclose(sorted_thermo_variances, expected_correlated_thermo_variances).all()
-        assert set(thermo_names) == expected_correlated_thermo_labels
+        assert sorted_thermo_names == expected_correlated_thermo_labels
 
+    def test_covariance_matrices(self):
+        """
+        Test that the covariance matrices are being constructed correctly, and that the correlated uncertainties are different from the uncorrelated ones
+        """
+
+        # have to add an extra reaction to see any kinetic correlations
+        # copy reaction 4 and change the index so it is a new reaction, but with the same source (rate rule) as the original reaction
+        extra_reaction = copy.deepcopy(self.uncertainty.reaction_list[4])
+        extra_reaction.index = len(self.uncertainty.reaction_list)
+        self.uncertainty.reaction_list.append(extra_reaction)
+        self.uncertainty.extract_sources_from_model()  # this will assign the same source to the new reaction as the original reaction
+
+        self.uncertainty.assign_parameter_uncertainties(correlated=False)
+        uncorrelated_thermo_inputs = np.array(self.uncertainty.thermo_input_uncertainties)
+        uncorrelated_kinetic_inputs = np.array(self.uncertainty.kinetic_input_uncertainties)
+
+        self.uncertainty.assign_intermediate_uncertainties(correlated=False)
+        uncorrelated_thermo_covariance = self.uncertainty.get_thermo_covariance_matrix()
+        uncorrelated_kinetic_covariance = self.uncertainty.get_kinetic_covariance_matrix()
+        
+        self.uncertainty.assign_intermediate_uncertainties(correlated=True)
+        correlated_thermo_covariance = self.uncertainty.get_thermo_covariance_matrix()
+        correlated_kinetic_covariance = self.uncertainty.get_kinetic_covariance_matrix()
+        Sigma_ww_thermo = self.uncertainty._get_intermediate_thermo_covariance_matrix()
+        Sigma_ww_kinetics = self.uncertainty._get_intermediate_kinetics_covariance_matrix()
+        
+        # check that the diagonal elements of the correlated and uncorrelated covariance matrices are the same and equal to the squares of the input uncertainties
+        np.testing.assert_allclose(np.diag(uncorrelated_thermo_covariance), np.float_power(uncorrelated_thermo_inputs, 2.0), rtol=1e-4)
+        np.testing.assert_allclose(np.diag(correlated_thermo_covariance), np.float_power(uncorrelated_thermo_inputs, 2.0), rtol=1e-4)
+        np.testing.assert_allclose(np.diag(uncorrelated_kinetic_covariance), np.float_power(uncorrelated_kinetic_inputs, 2.0), rtol=1e-4)
+        np.testing.assert_allclose(np.diag(correlated_kinetic_covariance), np.float_power(uncorrelated_kinetic_inputs, 2.0), rtol=1e-4)
+
+        # check that the off-diagonal elements of the uncorrelated covariance matrix are zero
+        off_diagonal_kinetic_uncorrelated = uncorrelated_kinetic_covariance - np.diag(np.diag(uncorrelated_kinetic_covariance))
+        assert np.allclose(off_diagonal_kinetic_uncorrelated, 0, atol=1e-8)
+        off_diagonal_thermo_uncorrelated = uncorrelated_thermo_covariance - np.diag(np.diag(uncorrelated_thermo_covariance))
+        assert np.allclose(off_diagonal_thermo_uncorrelated, 0, atol=1e-8)
+
+        # check that the off-diagonal elements of the correlated covariance matrix are not all zero
+        off_diagonal_kinetic_correlated = correlated_kinetic_covariance - np.diag(np.diag(correlated_kinetic_covariance))
+        assert not np.allclose(off_diagonal_kinetic_correlated, 0, atol=1e-8)
+        off_diagonal_thermo_correlated = correlated_thermo_covariance - np.diag(np.diag(correlated_thermo_covariance))
+        assert not np.allclose(off_diagonal_thermo_correlated, 0, atol=1e-8)
+
+        # check that the correlated covariance matrices are symmetric
+        assert np.allclose(correlated_kinetic_covariance, correlated_kinetic_covariance.T, atol=1e-8)
+        assert np.allclose(correlated_thermo_covariance, correlated_thermo_covariance.T, atol=1e-8)
+        assert np.allclose(Sigma_ww_kinetics, Sigma_ww_kinetics.T, atol=1e-8)
+        assert np.allclose(Sigma_ww_thermo, Sigma_ww_thermo.T, atol=1e-8)
+
+        # check that the matrix is positive semi-definite by confirming that all eigenvalues are non-negative
+        kinetic_eigenvalues = np.linalg.eigvals(correlated_kinetic_covariance)
+        assert np.all(kinetic_eigenvalues >= -1e-8)  # allow for small numerical errors
+        thermo_eigenvalues = np.linalg.eigvals(correlated_thermo_covariance)
+        assert np.all(thermo_eigenvalues >= -1e-8)  # allow for small numerical errors
+        intermediate_kinetic_eigenvalues = np.linalg.eigvals(Sigma_ww_kinetics)
+        assert np.all(intermediate_kinetic_eigenvalues >= -1e-8)
+        intermediate_thermo_eigenvalues = np.linalg.eigvals(Sigma_ww_thermo)
+        assert np.all(intermediate_thermo_eigenvalues >= -1e-8)
+
+        self.uncertainty.reaction_list.pop()  # remove the extra reaction so it doesn't affect other tests
 
     def test_specific_species_uncertainties(self):
         """