uncertainty: separate family and non-exact error

The previous implementation combined family and non-exact error into a total estimation error. But this requires assuming that family and non-exact error are perfectly correlated. I don't think that was the paper's original intention, nor do I think it makes sense, but more importantly the log10(N+1) doesn't play nicely with the new uncertainty formulation if you write it this way. If you look at the local analysis example in uncertaintyTest.py, you'll see that the correlated variance goes from about 3 to 1.8, which makes sense because we are no longer stacking non-exact on top of family, but considering them as unrelated.

Author: sevyharris

rmgpy/tools/uncertainty.py

@@ -223,7 +223,7 @@ def get_partial_uncertainty_value(self, source, corr_source_type, corr_param=Non
         Obtain the partial uncertainty dlnk/dlnk_corr*dlnk_corr, where dlnk_corr is the correlated parameter
 
         `corr_param` is the parameter identifier itself, which is the string identifier of the rate rule
-        `corr_source_type` is a string, being either 'Rate Rules', 'Library', 'PDep', 'Training' or 'Estimation'
+        `corr_source_type` is a string, being either 'Rate Rules', 'Library', 'PDep', 'Training', 'Estimation Nonexact', or 'Estimation Family'
         `corr_family` is a string used only when the source type is 'Rate Rules' and indicates the family
         """
 
@@ -275,24 +275,22 @@ def get_partial_uncertainty_value(self, source, corr_source_type, corr_param=Non
                     else:
                         return self.dlnk_training
 
-        elif corr_source_type == 'Estimation':
-            # Return all the uncorrelated uncertainty associated with using an estimation scheme
 
+        elif corr_source_type == 'Estimation Nonexact':
+            # Return the uncorrelated uncertainty associated with using a non-exact rate rule
             if 'Rate Rules' in source:
                 source_dict = source['Rate Rules'][1]
                 exact = source_dict['exact']
-
-                family_label = source['Rate Rules'][0]
-                dlnk = self.dlnk_family  # Base uncorrelated uncertainty just from using rate rule estimation
-
-                # Additional uncertainty from using non-exact rate rule
                 N = len(source_dict['rules']) + len(source_dict['training'])
                 if not exact:
                     # nonexactness contribution increases as N increases
-                    dlnk += np.log10(N + 1) * self.dlnk_nonexact
-                return dlnk
+                    return np.log10(N + 1) * self.dlnk_nonexact
+        elif corr_source_type == 'Estimation Family':
+            # Return the uncertainty associated with using a family decision tree
+            if 'Rate Rules' in source:
+                return self.dlnk_family
         else:
-            raise Exception('Kinetics correlated source must be Rate Rules, Library, PDep, Training, or Estimation')
+            raise Exception('Kinetics correlated source must be Rate Rules, Library, PDep, Training, Estimation Nonexact, or Estimation Family')
 
         # If we get here, it means that we did not find the correlated parameter in the source
         return None
@@ -733,11 +731,16 @@ def assign_parameter_uncertainties(self, g_param_engine=None, k_param_engine=Non
                         label = '{}Rate Rule {} {}'.format(surface_prefix, family, ruleEntry)
                         dlnk[label] = dplnk
 
-                    # There is also estimation error if rate rules are used
-                    est_dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Estimation')
-                    if est_dplnk:
-                        label = 'Estimation {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
-                        dlnk[label] = est_dplnk
+                    # There is also estimation error if rate rules are used (nonexact and family contribute to this)
+                    nonexact_dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Estimation Nonexact', corr_family=family)
+                    if nonexact_dplnk:
+                        label = 'Estimation Nonexact {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
+                        dlnk[label] = nonexact_dplnk
+
+                    family_dplnk = k_param_engine.get_partial_uncertainty_value(source, 'Estimation Family', corr_family=family)
+                    if family_dplnk:
+                        label = 'Estimation Family {}'.format(reaction.to_chemkin(self.species_list, kinetics=False))
+                        dlnk[label] = family_dplnk
 
                 elif 'PDep' in source:
                     dplnk = k_param_engine.get_partial_uncertainty_value(source, 'PDep', source['PDep'])

test/rmgpy/tools/uncertaintyTest.py

@@ -189,9 +189,9 @@ def test_local_analysis(self):
         expected_uncorrelated_total_variance = 1.8329056941266446
         expected_uncorrelated_thermo_variances = np.array([0.17092419, 0.09781627, 0.06186124, 0.04856985, 0.00391013, 0.00306632, 0.00041446, 9.953e-05])
         expected_uncorrelated_kinetics_variances = np.array([1.1311145, 0.15888459, 0.085189, 0.02022449, 0.01687337, 0.01605351, 0.01588366, 0.0010829, 0.00080811, 0.00012957])
-        expected_correlated_total_variance = 3.006637059881831
+        expected_correlated_total_variance = 1.7732795017083922
         expected_correlated_thermo_variances = np.array([0.09145902, 0.07672388, 0.04856985, 0.04573167, 0.03236887, 0.01747643, 0.01098087, 0.00143231, 0.0013764, 0.00031352, 0.00028968, 0.00014685, 0.0001338, 3.263e-05])
-        expected_correlated_kinetics_variances = np.array([2.0212174, 0.28939964, 0.15516713, 0.11981721, 0.02838292, 0.01687337, 0.0161796, 0.01605351, 0.0040449, 0.00253735, 0.00253735, 0.00193506, 0.00168253, 0.00136045, 0.00136045, 0.00080811, 0.00012957, 0.00011471])
+        expected_correlated_kinetics_variances = np.array([0.53202843, 0.47926886, 0.11981721, 0.11321232, 0.06070094, 0.04059757, 0.02176716, 0.01687337, 0.0161796, 0.01605351, 0.007471, 0.00673013, 0.0040449, 0.00253735, 0.00253735, 0.00168253, 0.00136045, 0.00136045, 0.00080811, 0.00050935, 0.00045884, 0.00012957, 0.00011471])
         expected_uncorrelated_thermo_labels = [
             'dln[C2H6(18)]/dG[CH(4)]',
             'dln[C2H6(18)]/dG[C2H3(20)]',
@@ -231,22 +231,27 @@ def test_local_analysis(self):
             'Estimation C2H5(12)',
         ]
         expected_correlated_kinetics_labels = [
-            'Estimation C2H5(12)+CH3CHCH3(21)<=>C2H6(18)+C3H6(22)',
-            'Estimation C2H6(18)+PC3H7(15)<=>C2H5(12)+C3H8(19)',
-            'Estimation C2H3(20)+C3H8(19)<=>C2H4(11)+PC3H7(15)',
+            'Estimation Nonexact C2H5(12)+CH3CHCH3(21)<=>C2H6(18)+C3H6(22)',
+            'Estimation Family C2H5(12)+CH3CHCH3(21)<=>C2H6(18)+C3H6(22)',
             'Rate Rule Disproportionation Root_Ext-2R!H-R_2R!H->C_4R->C',
-            'Estimation CH3(14)+C3H8(19)<=>CH4(16)+PC3H7(15)',
+            'Estimation Nonexact C2H6(18)+PC3H7(15)<=>C2H5(12)+C3H8(19)',
+            'Estimation Nonexact C2H3(20)+C3H8(19)<=>C2H4(11)+PC3H7(15)',
+            'Estimation Family C2H6(18)+PC3H7(15)<=>C2H5(12)+C3H8(19)',
+            'Estimation Family C2H3(20)+C3H8(19)<=>C2H4(11)+PC3H7(15)',
             'Library O(0)+H2O2(3)<=>OH(1)+HO2(2)',
-            'Estimation CH3(14)+PC3H7(15)<=>CH4(16)+CH2CH2CH2(17)',
+            'Estimation Family CH3(14)+PC3H7(15)<=>CH4(16)+CH2CH2CH2(17)',
             'PDep HCCO(10)(+M)<=>O(0)+C2H(8)(+M)',
+            'Estimation Nonexact CH3(14)+C3H8(19)<=>CH4(16)+PC3H7(15)',
+            'Estimation Family CH3(14)+C3H8(19)<=>CH4(16)+PC3H7(15)',
             'Rate Rule H_Abstraction C/H3/Cs;C_methyl',
             'Rate Rule H_Abstraction C/H3/Cs\\H3;C_rad/H2/Cs\\H\\Cs\\Cs|O',
             'Rate Rule H_Abstraction C/H3/Cs\\H3;C_rad/H2/Cs\\H3',
-            'Estimation C3H5(24)+CH2CH2CH2(17)<=>C3H5(23)+C3H6(22)',
             'Rate Rule H_Abstraction C/H3/Cs\\H2\\O;C_methyl',
             'Rate Rule H_Abstraction C/H3/Cs\\H3;Cd_Cd\\H2_pri_rad',
             'Rate Rule H_Abstraction C/H3/Cs\\H2\\Cs|O;Cd_Cd\\H2_rad/Cs',
             'Training Disproportionation CH3(14)+C2H5(12)<=>CH4(16)+C2H4(11)',
+            'Estimation Nonexact C3H5(24)+CH2CH2CH2(17)<=>C3H5(23)+C3H6(22)',
+            'Estimation Family C3H5(24)+CH2CH2CH2(17)<=>C3H5(23)+C3H6(22)',
             'Training H_Abstraction CH3(14)+C2H6(18)<=>CH4(16)+C2H5(12)',
             'Rate Rule Disproportionation Root_Ext-1R!H-R_N-4R->O_N-Sp-5R!H=1R!H_Ext-4CHNS-R_N-6R!H->S_4CHNS->C_N-Sp-6BrBrBrCCCClClClFFFIIINNNOOOPPPSiSiSi#4C_6BrCClFINOPSi->C_N-1R!H-inRing_N-Sp-6C-4C',
         ]