Fix tree averaging issue #712

rmgpy/data/base.py

@@ -451,7 +451,7 @@ def loadOldDictionary(self, path, pattern):
 
     def __loadTree(self, tree):
         """
-        Parse an old-style RMG tree located at `tree`. An RMG tree is an n-ary
+        Parse an group tree located at `tree`. An RMG tree is an n-ary
         tree representing the hierarchy of items in the dictionary.
         """
 
@@ -498,7 +498,10 @@ def __loadTree(self, tree):
                 else:
                     entry.parent = None
                     self.top.append(entry)
-
+
+                # Save the level of the tree into the entry
+                entry.level = level
+
                 # Add node to list of parents for subsequent iteration
                 parents.append(label)
 

rmgpy/data/kinetics/family.py

@@ -1039,15 +1039,15 @@ def getRootTemplate(self):
         else:
             return self.groups.top
 
-    def fillKineticsRulesByAveragingUp(self, rootTemplate=None, alreadyDone=None):
+    def fillKineticsRulesByAveragingUp(self):
         """
-        Fill in gaps in the kinetics rate rules by averaging child nodes.
+        Fill in gaps in the kinetics rate rules by averaging child nodes
+        recursively starting from the top level root template.
         """
-        # If no template is specified, then start at the top-level nodes
-        if rootTemplate is None:
-            rootTemplate = self.getRootTemplate()
-            alreadyDone = {}
-        self.rules.fillRulesByAveragingUp(rootTemplate, alreadyDone)
+
+        self.rules.fillRulesByAveragingUp(self.getRootTemplate(), {})
+
+
 
     def applyRecipe(self, reactantStructures, forward=True, unique=True):
         """

rmgpy/data/kinetics/rules.py

@@ -448,33 +448,26 @@ def fillRulesByAveragingUp(self, rootTemplate, alreadyDone):
 
         if rootLabel in alreadyDone:
             return alreadyDone[rootLabel]
-
-        # See if we already have a rate rule for this exact template 
-        entry = self.getRule(rootTemplate)
-        if entry is not None and entry.rank > 0:
-            # We already have a rate rule for this exact template
-            # If the entry has rank of zero, then we have so little faith
-            # in it that we'd rather use an averaged value if possible
-            # Since this entry does not have a rank of zero, we keep its
-            # value
-            alreadyDone[rootLabel] = entry.data
-            return entry.data
-
-        # Recursively descend to the child nodes
-        childrenList = [[group] for group in rootTemplate]
-        for group in childrenList:
-            parent = group.pop(0)
-            if len(parent.children) > 0:
-                group.extend(parent.children)
-            else:
-                group.append(parent)
-
-        childrenList = getAllCombinations(childrenList)
+
+        # Generate the distance 1 pairings which must be averaged for this root template.
+        # The distance 1 template is created by taking the parent node from one or more trees
+        # and creating the combinations with children from a single remaining tree.  
+        # i.e. for some node (A,B), we want to fetch all combinations for the pairing of (A,B's children) and
+        # (A's children, B).  For node (A,B,C), we would retrieve all combinations of (A,B,C's children) 
+        # (A,B's children,C) etc...  
+        # If a particular node has no children, it is skipped from the children expansion altogether.
+
+        childrenList = []
+        for i, parent in enumerate(rootTemplate):
+            # Start with the root template, and replace the ith member with its children
+            if parent.children:
+                childrenSet = [[group] for group in rootTemplate]
+                childrenSet[i] = parent.children
+                childrenList.extend(getAllCombinations(childrenSet))
+
         kineticsList = []
         for template in childrenList:
             label = ';'.join([g.label for g in template])
-            if template == rootTemplate: 
-                continue
 
             if label in alreadyDone:
                 kinetics = alreadyDone[label]
@@ -483,17 +476,45 @@ def fillRulesByAveragingUp(self, rootTemplate, alreadyDone):
 
             if kinetics is not None:
                 kineticsList.append([kinetics, template])
+
+        # See if we already have a rate rule for this exact template instead
+        # and return it now that we have finished searching its children
+        entry = self.getRule(rootTemplate)
+        if entry is not None and entry.rank > 0:
+            # We already have a rate rule for this exact template
+            # If the entry has rank of zero, then we have so little faith
+            # in it that we'd rather use an averaged value if possible
+            # Since this entry does not have a rank of zero, we keep its
+            # value
+            alreadyDone[rootLabel] = entry.data
+            return entry.data
 
         if len(kineticsList) > 0:
 
-            # We found one or more results! Let's average them together
-            kinetics = self.__getAverageKinetics([k for k, t in kineticsList])
             if len(kineticsList) > 1:
-                kinetics.comment += 'Average of ({0})'.format(
-                    ' + '.join(k.comment if k.comment != '' else ';'.join(g.label for g in t) for k, t in kineticsList))
+                # We found one or more results! Let's average them together
+                kinetics = self.__getAverageKinetics([k for k, t in kineticsList])
+                kinetics.comment = 'Average of ({0})'.format(
+                     ' + '.join(';'.join(g.label for g in t) for k, t in kineticsList))
+
+                # For debug mode: uncomment the following kinetics commenting
+                # lines and use them instead of the lines above. Caution: large memory usage.
+
+                # kinetics.comment += 'Average of ({0})'.format(
+                #     ' + '.join(k.comment if k.comment != '' else ';'.join(g.label for g in t) for k, t in kineticsList))
+
             else:
                 k,t = kineticsList[0]
-                kinetics.comment += k.comment if k.comment != '' else ';'.join(g.label for g in t)
+                kinetics = deepcopy(k)
+                # Even though we are using just a single set of kinetics, it's still considered
+                # an average.  It just happens that the other distance 1 children had no data.
+                kinetics.comment = 'Average of ({0})'.format(';'.join(g.label for g in t))
+
+                # For debug mode: uncomment the following kinetics commenting
+                # lines and use them instead of the lines above. Caution: large memory usage.
+
+                # kinetics.comment += 'Average of ({0}).format(k.comment if k.comment != '' else ';'.join(g.label for g in t))
+
             entry = Entry(
                 index = 0,
                 label = rootLabel,
@@ -544,7 +565,25 @@ def __getAverageKinetics(self, kineticsList):
             E0 = (E0*0.001,"kJ/mol"),
         )
         return averagedKinetics
+
+    def calculateNormDistance(self, template, otherTemplate):
+        """
+        Calculate the norm distance squared between two rate rules with
+        `template` and `otherTemplate`.  The norm distance is 
+        a^2 + b^2 + c^2 .... when a is the distance between the nodes in the
+        first tree, b is the distance between the nodes in the second tree, etc.
+        """
+
+        # Do it the stupid way first and calculate distances from the top 
+        # rather than from each other for now... it's dumb but need to see results first
+        import numpy
+        depth = numpy.array([node.level for node in template])
+        otherDepth = numpy.array([otherNode.level for otherNode in otherTemplate])
 
+        distance = numpy.array(depth-otherDepth)
+        norm = numpy.dot(distance,distance)
+        return norm
+
     def estimateKinetics(self, template, degeneracy=1):
         """
         Determine the appropriate kinetics for a reaction with the given
@@ -567,7 +606,14 @@ def getTemplateLabel(template):
                 kineticsList.append([kinetics, t])
 
             if len(kineticsList) > 0:                 
-
+
+                if len(kineticsList) > 1:
+                    # Filter the kinetics to use templates with the lowest minimum euclidean distance 
+                    # from the specified template
+                    norms = [self.calculateNormDistance(template, t) for kinetics,t in kineticsList]
+                    minNorm = min(norms) 
+                    kineticsList = [pair for pair, norm in zip(kineticsList,norms) if norm == min(norms)]
+
                 if len(kineticsList) == 1:
                     kinetics, t = kineticsList[0]
                     # Check whether the exact rate rule for the original template (most specific