Changed toleranceKeepInEdge, which was initially a flux, into a ratio of edge rate against characteristic rate!

A unitless tolerance makes more sense than flux.

Author: KEHANG

rmgpy/rmg/model.py

@@ -1043,7 +1043,7 @@ def addSpeciesToEdge(self, spec):
         """
         self.edge.species.append(spec)
 
-    def prune(self, reactionSystems, fluxToleranceKeepInEdge, maximumEdgeSpecies):
+    def prune(self, reactionSystems, toleranceKeepInEdge, maximumEdgeSpecies):
         """
         Remove species from the model edge based on the simulation results from
         the list of `reactionSystems`.
@@ -1062,55 +1062,57 @@ def prune(self, reactionSystems, fluxToleranceKeepInEdge, maximumEdgeSpecies):
 
         # Get the maximum species rates (and network leak rates)
         # across all reaction systems
-        maxEdgeSpeciesRates = numpy.zeros((numEdgeSpecies), numpy.float64)
+        maxEdgeSpeciesRateRatios = numpy.zeros((numEdgeSpecies), numpy.float64)
         for reactionSystem in reactionSystems:
             for i in range(numEdgeSpecies):
-                rate = reactionSystem.maxEdgeSpeciesRates[i]
-                if maxEdgeSpeciesRates[i] < rate:
-                    maxEdgeSpeciesRates[i] = rate
+                rateRatio = reactionSystem.maxEdgeSpeciesRateRatios[i]
+                if maxEdgeSpeciesRateRatios[i] < rateRatio:
+                    maxEdgeSpeciesRateRatios[i] = rateRatio
 
             for i in range(len(self.networkList)):
                 network = self.networkList[i]
-                rate = reactionSystem.maxNetworkLeakRates[i]
+                rateRatio = reactionSystem.maxNetworkLeakRateRatios[i]
                 # Add the fraction of the network leak rate contributed by
                 # each unexplored species to that species' rate
                 # This is to ensure we have an overestimate of that species flux
                 ratios = network.getLeakBranchingRatios(reactionSystem.T.value_si,reactionSystem.P.value_si)
                 for spec, frac in ratios.iteritems():
                     index = self.edge.species.index(spec)
-                    maxEdgeSpeciesRates[index] += frac * rate
+                    maxEdgeSpeciesRateRatios[index] += frac * rateRatio
                 # Mark any species that is explored in any partial network as ineligible for pruning
                 for spec in network.explored:
                     if spec not in ineligibleSpecies:
                         ineligibleSpecies.append(spec)
 
         # Sort the edge species rates by index
-        indices = numpy.argsort(maxEdgeSpeciesRates)
-
+        indices = numpy.argsort(maxEdgeSpeciesRateRatios)
         # Determine which species to prune
         speciesToPrune = []
         pruneDueToRateCounter = 0
         for index in indices:
             # Remove the species with rates below the pruning tolerance from the model edge
-            if maxEdgeSpeciesRates[index] < fluxToleranceKeepInEdge and self.edge.species[index] not in ineligibleSpecies:
+            if maxEdgeSpeciesRateRatios[index] < toleranceKeepInEdge and self.edge.species[index] not in ineligibleSpecies:
                 speciesToPrune.append((index, self.edge.species[index]))
                 pruneDueToRateCounter += 1
             # Keep removing species with the lowest rates until we are below the maximum edge species size
             elif numEdgeSpecies - len(speciesToPrune) > maximumEdgeSpecies and self.edge.species[index] not in ineligibleSpecies:
+                logging.info('Pruning species {0} to make numEdgeSpecies smaller than maximumEdgeSpecies'.format(self.edge.species[index]))
                 speciesToPrune.append((index, self.edge.species[index]))
             else:
-                break
+                continue
 
         # Actually do the pruning
         if pruneDueToRateCounter > 0:
-            logging.info('Pruning {0:d} species whose rates did not exceed the minimum threshold of {1:g}'.format(pruneDueToRateCounter, fluxToleranceKeepInEdge))
+            logging.info('Pruning {0:d} species whose rate ratios against characteristic rate did not exceed the minimum threshold of {1:g}'.format(pruneDueToRateCounter, toleranceKeepInEdge))
             for index, spec in speciesToPrune[0:pruneDueToRateCounter]:
-                logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRates[index]))
+                logging.info('Pruning species {0:<56}'.format(spec))
+                logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRateRatios[index]))
                 self.removeSpeciesFromEdge(spec)
         if len(speciesToPrune) - pruneDueToRateCounter > 0:
             logging.info('Pruning {0:d} species to obtain an edge size of {1:d} species'.format(len(speciesToPrune) - pruneDueToRateCounter, maximumEdgeSpecies))
             for index, spec in speciesToPrune[pruneDueToRateCounter:]:
-                logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRates[index]))
+                logging.info('Pruning species {0:<56}'.format(spec))
+                logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRateRatios[index]))
                 self.removeSpeciesFromEdge(spec)
 
         # Delete any networks that became empty as a result of pruning

rmgpy/solver/base.pxd

@@ -42,6 +42,8 @@ cdef class ReactionSystem(DASSL):
     cdef public numpy.ndarray maxCoreSpeciesRates
     cdef public numpy.ndarray maxEdgeSpeciesRates
     cdef public numpy.ndarray maxNetworkLeakRates
+    cdef public numpy.ndarray maxEdgeSpeciesRateRatios
+    cdef public numpy.ndarray maxNetworkLeakRateRatios
     cdef public numpy.ndarray sensitivityCoefficients
 
     cdef public list termination

rmgpy/solver/base.pyx

@@ -64,6 +64,8 @@ cdef class ReactionSystem(DASSL):
         self.maxCoreSpeciesRates = None
         self.maxEdgeSpeciesRates = None
         self.maxNetworkLeakRates = None
+        self.maxEdgeSpeciesRateRatios = None
+        self.maxNetworkLeakRateRatios = None
         self.sensitivityCoefficients = None
         self.termination = termination or []
 
@@ -93,6 +95,8 @@ cdef class ReactionSystem(DASSL):
         self.maxCoreSpeciesRates = numpy.zeros((numCoreSpecies), numpy.float64)
         self.maxEdgeSpeciesRates = numpy.zeros((numEdgeSpecies), numpy.float64)
         self.maxNetworkLeakRates = numpy.zeros((numPdepNetworks), numpy.float64)
+        self.maxEdgeSpeciesRateRatios = numpy.zeros((numEdgeSpecies), numpy.float64)
+        self.maxNetworkLeakRateRatios = numpy.zeros((numPdepNetworks), numpy.float64)
         self.sensitivityCoefficients = numpy.zeros((numCoreSpecies, numCoreReactions), numpy.float64)
 
 
@@ -126,7 +130,7 @@ cdef class ReactionSystem(DASSL):
         cdef double stepTime, charRate, maxSpeciesRate, maxNetworkRate, prevTime, volume
         cdef numpy.ndarray[numpy.float64_t, ndim=1] y0, dydk  #: Vector containing the number of moles of each species
         cdef numpy.ndarray[numpy.float64_t, ndim=1] coreSpeciesRates, edgeSpeciesRates, networkLeakRates
-        cdef numpy.ndarray[numpy.float64_t, ndim=1] maxCoreSpeciesRates, maxEdgeSpeciesRates, maxNetworkLeakRates
+        cdef numpy.ndarray[numpy.float64_t, ndim=1] maxCoreSpeciesRates, maxEdgeSpeciesRates, maxNetworkLeakRates,maxEdgeSpeciesRateRatios, maxNetworkLeakRateRatios
         cdef bint terminated
         cdef object maxSpecies, maxNetwork
         cdef int iteration, i
@@ -158,6 +162,8 @@ cdef class ReactionSystem(DASSL):
         maxCoreSpeciesRates = self.maxCoreSpeciesRates
         maxEdgeSpeciesRates = self.maxEdgeSpeciesRates
         maxNetworkLeakRates = self.maxNetworkLeakRates
+        maxEdgeSpeciesRateRatios = self.maxEdgeSpeciesRateRatios
+        maxNetworkLeakRateRatios = self.maxNetworkLeakRateRatios
 
         # Copy the initial conditions to use in evaluating conversions
         y0 = self.y.copy()
@@ -220,6 +226,8 @@ cdef class ReactionSystem(DASSL):
             coreSpeciesRates = numpy.abs(self.coreSpeciesRates)
             edgeSpeciesRates = numpy.abs(self.edgeSpeciesRates)
             networkLeakRates = numpy.abs(self.networkLeakRates)
+            edgeSpeciesRateRatios = numpy.abs(self.edgeSpeciesRates/charRate)
+            networkLeakRateRatios = numpy.abs(self.networkLeakRates/charRate)
 
             # Update the maximum species rate and maximum network leak rate arrays
             for index in range(numCoreSpecies):
@@ -231,6 +239,12 @@ cdef class ReactionSystem(DASSL):
             for index in range(numPdepNetworks):
                 if maxNetworkLeakRates[index] < networkLeakRates[index]:
                     maxNetworkLeakRates[index] = networkLeakRates[index]
+            for index in range(numEdgeSpecies):
+                if maxEdgeSpeciesRateRatios[index] < edgeSpeciesRateRatios[index]:
+                    maxEdgeSpeciesRateRatios[index] = edgeSpeciesRateRatios[index]
+            for index in range(numPdepNetworks):
+                if maxNetworkLeakRateRatios[index] < networkLeakRateRatios[index]:
+                    maxNetworkLeakRateRatios[index] = networkLeakRateRatios[index]
 
             # Get the edge species with the highest flux
             if numEdgeSpecies > 0:
@@ -313,6 +327,8 @@ cdef class ReactionSystem(DASSL):
         self.maxCoreSpeciesRates = maxCoreSpeciesRates
         self.maxEdgeSpeciesRates = maxEdgeSpeciesRates
         self.maxNetworkLeakRates = maxNetworkLeakRates
+        self.maxEdgeSpeciesRateRatios = maxEdgeSpeciesRateRatios
+        self.maxNetworkLeakRateRatios = maxNetworkLeakRateRatios
 
         # Return the invalid object (if the simulation was invalid) or None
         # (if the simulation was valid)