make uncertainties add in quadrature

Author: sevyharris

rmgpy/tools/uncertainty.py

@@ -43,7 +43,7 @@ class ThermoParameterUncertainty(object):
     This class is an engine that generates the species uncertainty based on its thermo sources.
     """
 
-    def __init__(self, dG_library=1.5, dG_QM=3.0, dG_GAV=1.5, dG_group=0.1, dG_ADS_correction=6.918, dG_surf_lib=6.918):
+    def __init__(self, dG_library=1.5, dG_QM=3.0, dG_GAV=1.5, dG_group=0.316, dG_ADS_correction=6.918, dG_surf_lib=6.918):
         """
         Initialize the different uncertainties dG_library, dG_QM, dG_GAV, and dG_other with set values
         in units of kcal/mol.
@@ -62,22 +62,22 @@ def get_uncertainty_value(self, source):
         """
         Retrieve the uncertainty value in kcal/mol when the source of the thermo of a species is given.
         """
-        dG = 0.0
+        varG = 0.0
         if 'Library' in source:
-            dG += self.dG_library
+            varG += self.dG_library ** 2
         if 'Surface_Library' in source:
-            dG += self.dG_surf_lib
+            varG += self.dG_surf_lib ** 2
         if 'QM' in source:
-            dG += self.dG_QM
+            varG += self.dG_QM ** 2
         if 'GAV' in source:
-            dG += self.dG_GAV  # Add a fixed uncertainty for the GAV method
+            varG += self.dG_GAV ** 2  # Add a fixed uncertainty for the GAV method
             for group_type, group_entries in source['GAV'].items():
                 group_weights = [groupTuple[-1] for groupTuple in group_entries]
-                dG += np.sum([weight * self.dG_group for weight in group_weights])
+                varG += np.sum([weight ** 2 * self.dG_group ** 2 for weight in group_weights])
         if 'ADS' in source:
-            dG += self.dG_ADS_correction  # Add adsorption correction uncertainty
+            varG += self.dG_ADS_correction ** 2  # Add adsorption correction uncertainty
 
-        return dG
+        return np.sqrt(varG)
 
     def get_partial_uncertainty_value(self, source, corr_source_type, corr_param=None, corr_group_type=None):
         """
@@ -169,21 +169,21 @@ def get_uncertainty_value(self, source):
         """
         Retrieve the dlnk uncertainty when the source of the reaction kinetics are given
         """
-        dlnk = 0.0
+        varlnk = 0.0
         if 'Library' in source:
             # Should be a single library reaction source
-            dlnk += self.dlnk_library
+            varlnk += self.dlnk_library ** 2
         elif 'Surface_Library' in source:
             # Should be a single library reaction source
-            dlnk += self.dlnk_Surface_Library
+            varlnk += self.dlnk_Surface_Library ** 2
         elif 'PDep' in source:
             # Should be a single pdep reaction source
-            dlnk += self.dlnk_pdep
+            varlnk += self.dlnk_pdep ** 2
         elif 'Training' in source:
             # Should be a single training reaction
             # Although some training entries may be used in reverse,
             # We still consider the kinetics to be directly dependent
-            dlnk += self.dlnk_training
+            varlnk += self.dlnk_training ** 2
         elif 'Rate Rules' in source:
             family_label = source['Rate Rules'][0]
             source_dict = source['Rate Rules'][1]
@@ -192,16 +192,16 @@ def get_uncertainty_value(self, source):
             training_weights = [trainingTuple[-1] for trainingTuple in source_dict['training']]
 
             if 'surface' in family_label.lower():
-                dlnk += self.dlnk_family_surface
+                varlnk += self.dlnk_family_surface ** 2
             else:
-                dlnk += self.dlnk_family
+                varlnk += self.dlnk_family ** 2
 
             N = len(rule_weights) + len(training_weights)
             if 'node_std_dev' in source_dict:
                 # Handle autogen BM trees
                 if source_dict['node_std_dev'] < 0:
                     raise ValueError('Invalid value for std dev of kinetics family rule node')
-                dlnk += source_dict['node_std_dev']
+                varlnk += np.float_power(source_dict['node_std_dev'], 2.0)
                 if source_dict['node_n_train'] < 0:
                     raise ValueError('Invalid number of training reactions for kinetics family rule node')
                 N = source_dict['node_n_train']
@@ -210,24 +210,24 @@ def get_uncertainty_value(self, source):
                 # every node template has its own fitted rate rule by definition, but here we use the
                 # number of training reactions as an approximation of the node's specificity/generality
                 # and add a penalty for being too general (large # of training reactions)
-                dlnk += np.log10(N + 1) * self.dlnk_nonexact
+                varlnk += (np.log10(N + 1) * self.dlnk_nonexact) ** 2
             else:
                 # Handle hand-made trees
                 if not exact:
                     # nonexactness contribution increases as N increases
-                    dlnk += np.log10(N + 1) * self.dlnk_nonexact
+                    varlnk += (np.log10(N + 1) * self.dlnk_nonexact) ** 2
 
                 # Add the contributions from rules
-                dlnk += np.sum([weight * self.dlnk_rule for weight in rule_weights])
+                varlnk += np.sum([weight ** 2 * self.dlnk_rule ** 2 for weight in rule_weights])
                 # Add the contributions from training
                 # Even though these source from training reactions, we actually
                 # use the uncertainty for rate rules, since these are now approximations
                 # of the original reaction.  We consider these to be independent of original the training
                 # parameters because the rate rules may be reversing the training reactions,
                 # which leads to more complicated dependence
-                dlnk += np.sum([weight * self.dlnk_rule for weight in training_weights])
+                varlnk += np.sum([weight ** 2 * self.dlnk_rule ** 2 for weight in training_weights])
 
-        return dlnk
+        return np.sqrt(varlnk)
 
     def get_partial_uncertainty_value(self, source, corr_source_type, corr_param=None, corr_family=None):
         """

test/rmgpy/tools/uncertaintyTest.py

@@ -169,12 +169,11 @@ def test_uncertainty_assignment(self):
 
         np.testing.assert_allclose(
             thermo_unc,
-            [1.5, 1.5, 2.0, 1.9, 3.1, 1.5, 1.9, 2.0, 2.0, 1.9, 2.2, 1.9, 2.0, 1.5, 3.1, 1.9, 1.5, 2.0, 1.7, 1.8, 1.8, 1.9, 1.8, 1.9, 1.9],
+            [1.5, 1.5, 1.7745, 1.7461, 2.1447, 1.5, 1.6879, 1.7173, 1.7745, 1.7461, 1.7745, 1.7461, 1.7745, 1.5, 2.1447, 1.6879, 1.5, 1.7745, 1.6277, 1.6581, 1.6581, 1.6879, 1.5967, 1.6277, 1.6277],
             rtol=1e-4,
         )
         np.testing.assert_allclose(
             kinetic_unc,
-            # the 13.598 value comes from the SIDT tree node uncertainty: 11.54 + non-exact penalty for N=1: log10(1+1) * 3.5 + family uncertainty: 1.0
-            [0.5, 1.5, 3.169924, 3.169924, 2.553605, 0.5, 2.0, 13.5938, 13.5938, 0.5],
+            [0.5, 1.118, 1.9783, 1.9783, 1.5363, 0.5, 2.0, 11.6312, 11.6312, 0.5],
             rtol=1e-4
         )