uncertainty uses recorded std dev for SIDT nodes

Up until now the uncertainty tool has assumed the kinetics trees are the original hand-built type. This enables the uncertainty tool to grab the calculated standard deviation and number of training reactions from the specific autogen tree node used, instead of incorrectly treating it as an exact match of a single rate rule with a std dev of 0.5

Author: sevyharris

rmgpy/data/kinetics/family.py

@@ -4440,8 +4440,9 @@ def extract_source_from_comments(self, reaction):
             [Family_Label, {'template': originalTemplate,
                             'degeneracy': degeneracy,
                             'exact': boolean_exact?,
-                            'rules': a list of (original rate rule entry, weight in average)
-                            'training': a list of (original rate rule entry associated with training entry, original training entry, weight in average)}]
+                            'rules': a list of (original rate rule entry, weight in average),
+                            'training': a list of (original rate rule entry associated with training entry, original training entry, weight in average),
+                            'autogenerated': boolean for whether kinetics come from autogenerated subgraph isomorphic decision tree}]
 
 
         where Exact is a boolean of whether the rate is an exact match, Template is
@@ -4455,6 +4456,7 @@ def extract_source_from_comments(self, reaction):
         rules = None
         training_entries = None
         degeneracy = 1
+        autogenerated = False
 
         training_reaction_pattern = r'Matched reaction\s*(\d+).*in.*training'
         degeneracy_pattern = r'Multiplied by reaction path degeneracy\s*(\d+)'
@@ -4494,6 +4496,7 @@ def extract_source_from_comments(self, reaction):
         autogen_node_matches = re.search(autogen_node_search_pattern, full_comment_string)
         template_matches = re.search(template_pattern, full_comment_string)
         if autogen_node_matches is not None:  # autogenerated trees
+            autogenerated = True
             template_str = autogen_node_matches.group(1).split('Multiplied by reaction path degeneracy')[0].strip()
             template_str = template_str.split('in family')[0].strip()
             tokens = template_str.split()
@@ -4510,7 +4513,7 @@ def extract_source_from_comments(self, reaction):
             raise ValueError(f'Could not find rate rule in comments for reaction {reaction}.')
         rules, training_entries = self.get_sources_for_template(template)
         source_dict = {'template': template, 'degeneracy': degeneracy, 'exact': exact_rule,
-                       'rules': rules, 'training': training_entries}
+                       'rules': rules, 'training': training_entries, 'autogenerated': autogenerated}
 
         # Source of the kinetics is from rate rules
         return False, [self.label, source_dict]

rmgpy/tools/uncertainty.py

@@ -28,6 +28,7 @@
 ###############################################################################
 
 import os
+import re
 
 import numpy as np
 
@@ -165,19 +166,35 @@ def get_uncertainty_value(self, source):
 
             dlnk += self.dlnk_family ** 2
             N = len(rule_weights) + len(training_weights)
-            if not exact:
-                # nonexactness contribution increases as N increases
+            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']
+                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']
+
+                # Technically every lookup in the autogenerated trees is an "exact" match because
+                # 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
+            else:
+                # Handle hand-made trees
+                if not exact:
+                    # nonexactness contribution increases as N increases
+                    dlnk += np.log10(N + 1) * self.dlnk_nonexact
 
-            # Add the contributions from rules
-            dlnk += np.sum([weight * self.dlnk_rule 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])
+                # Add the contributions from rules
+                dlnk += np.sum([weight * self.dlnk_rule 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])
 
         return dlnk
 
@@ -413,8 +430,20 @@ def extract_sources_from_model(self):
                 # Do nothing here because training source already saves the entry from the training reaction
                 pass
             elif 'Rate Rules' in source:
-                # Do nothing
-                pass
+                # Fetch standard deviation if autogenerated tree
+                if source['Rate Rules'][1]['autogenerated']:
+                    long_desc = source['Rate Rules'][1]['rules'][0][0].long_desc
+                    std_dev_matches = re.search(r'Standard Deviation in ln\(k\): ([0-9]*.[0-9]*)', long_desc)
+                    std_dev = 1.329  # Default value is uniform distribution with upper bound of 10x the nominal value
+                    if std_dev_matches is not None:
+                        std_dev = float(std_dev_matches[1])
+
+                    n_train_matches = re.search('rule fitted to ([0-9]*) training reactions', long_desc)
+                    n_train = -1
+                    if n_train_matches is not None:
+                        n_train = int(n_train_matches[1])
+                    source['Rate Rules'][1]['node_std_dev'] = std_dev
+                    source['Rate Rules'][1]['node_n_train'] = n_train
             else:
                 raise Exception('Source of kinetics must be either Library, PDep, Training, or Rate Rules')
             self.reaction_sources_dict[reaction] = source

test/rmgpy/tools/uncertaintyTest.py

@@ -174,6 +174,7 @@ def test_uncertainty_assignment(self):
         )
         np.testing.assert_allclose(
             kinetic_unc,
-            [0.5, 1.5, 3.169924, 3.169924, 2.553605, 0.5, 2.0, 2.553605, 2.553605, 0.5],
+            # 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],
             rtol=1e-4
         )