Merge pull request #1705 from ReactionMechanismGenerator/catfix

Fix a few bugs and add tests: KineticsModel and uncertainty handling.

rmgpy/kinetics/chebyshev.pyx

@@ -249,7 +249,7 @@ cdef class Chebyshev(PDepKineticsModel):
             return False
         if self.degreeT != other_kinetics.degreeT or self.degreeP != other_kinetics.degreeP:
             return False
-        if self.kunits != other_kinetics.kunits or (self.coeffs != other_kinetics.coeffs).any():
+        if self.kunits != other_kinetics.kunits or not np.array_equal(self.coeffs.value, other_kinetics.coeffs.value):
             return False
 
         return True

rmgpy/kinetics/chebyshevTest.py

@@ -239,6 +239,96 @@ def test_change_rate(self):
             kact = self.chebyshev.get_rate_coefficient(T, 1e5)
             self.assertAlmostEqual(2 * kexp, kact, delta=1e-6 * kexp)
 
+    def test_is_identical_to(self):
+        """
+        Test the Chebyshev.is_identical_to() method.
+        """
+        # Trivial case, compare to a KineticsModel
+        from rmgpy.kinetics.model import KineticsModel
+        self.assertFalse(self.chebyshev.is_identical_to(KineticsModel()))
+
+        # Compare to identical Chebyshev
+        new_chebyshev = Chebyshev(
+            coeffs=self.coeffs,
+            kunits="cm^3/(mol*s)",
+            Tmin=(self.Tmin, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertTrue(self.chebyshev.is_identical_to(new_chebyshev))
+
+        # Compare to Chebyshev with different Tmin/Tmax
+        new_chebyshev = Chebyshev(
+            coeffs=self.coeffs,
+            kunits="cm^3/(mol*s)",
+            Tmin=(200, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertFalse(self.chebyshev.is_identical_to(new_chebyshev))
+
+        new_chebyshev = Chebyshev(
+            coeffs=self.coeffs,
+            kunits="cm^3/(mol*s)",
+            Tmin=(self.Tmin, "K"),
+            Tmax=(2500, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertFalse(self.chebyshev.is_identical_to(new_chebyshev))
+
+        # Compare to Chebyshev with different degreeT/degreeP
+        new_chebyshev = Chebyshev(
+            coeffs=self.coeffs[0:-1, :],  # Remove one T dimension
+            kunits="cm^3/(mol*s)",
+            Tmin=(self.Tmin, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertFalse(self.chebyshev.is_identical_to(new_chebyshev))
+
+        new_chebyshev = Chebyshev(
+            coeffs=self.coeffs[:, 0:-1],  # Remove one P dimension
+            kunits="cm^3/(mol*s)",
+            Tmin=(self.Tmin, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertFalse(self.chebyshev.is_identical_to(new_chebyshev))
+
+        # Compare to Chebyshev with different units
+        new_chebyshev = Chebyshev(
+            coeffs=self.coeffs,
+            kunits="m^3/(mol*s)",
+            Tmin=(self.Tmin, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertFalse(self.chebyshev.is_identical_to(new_chebyshev))
+
+        # Compare to Chebyshev with slightly different coefficients
+        new_chebyshev = Chebyshev(
+            coeffs=np.copy(self.coeffs) * 0.01,
+            kunits="cm^3/(mol*s)",
+            Tmin=(self.Tmin, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            comment=self.comment,
+        )
+        self.assertFalse(self.chebyshev.is_identical_to(new_chebyshev))
+
 
 ################################################################################
 

rmgpy/kinetics/model.pyx

@@ -121,7 +121,7 @@ cdef class KineticsModel:
         Return a string representation that can be used to reconstruct the
         KineticsModel object.
         """
-        return 'KineticsModel(Tmin={0!r}, Tmax={1!r}, Pmin={0!r}, Pmax={1!r}, comment="""{2}""")'.format(
+        return 'KineticsModel(Tmin={0!r}, Tmax={1!r}, Pmin={2!r}, Pmax={3!r}, uncertainty={4!r}, comment="""{5}""")'.format(
             self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.uncertainty, self.comment)
 
     def __reduce__(self):
@@ -230,15 +230,15 @@ cdef class KineticsModel:
         Returns ``True`` if Tmin, Tmax for both objects match.
         Otherwise returns ``False``
         """
-        if self.Tmin is not None and other_kinetics.Tmin is not None and not self.Tmin.equals(other_kinetics.Tmin):
-            return False
+        if self.Tmin is not None and other_kinetics.Tmin is not None and self.Tmin.equals(other_kinetics.Tmin):
+            pass
         elif self.Tmin is None and other_kinetics.Tmin is None:
             pass
         else:
             return False
 
-        if self.Tmax is not None and other_kinetics.Tmax is not None and not self.Tmax.equals(other_kinetics.Tmax):
-            return False
+        if self.Tmax is not None and other_kinetics.Tmax is not None and self.Tmax.equals(other_kinetics.Tmax):
+            pass
         elif self.Tmax is None and other_kinetics.Tmax is None:
             pass
         else:

rmgpy/kinetics/modelTest.py

@@ -35,7 +35,76 @@
 import unittest
 
 from rmgpy.kinetics.model import get_reaction_order_from_rate_coefficient_units, \
-                                 get_rate_coefficient_units_from_reaction_order
+                                 get_rate_coefficient_units_from_reaction_order, \
+                                 KineticsModel
+from rmgpy.kinetics.uncertainties import RateUncertainty
+
+
+class TestKineticsModel(unittest.TestCase):
+    """
+    Contains unit tests of the KineticsModel class
+    """
+
+    def setUp(self):
+        """
+        A function run before each unit test in this class.
+        """
+        self.Tmin = 300.
+        self.Tmax = 3000.
+        self.Pmin = 0.1
+        self.Pmax = 100.
+        self.comment = 'foo bar'
+        self.uncertainty = RateUncertainty(mu=0.3, var=0.6, Tref=1000.0, N=1, correlation="ab")
+        self.km = KineticsModel(
+            Tmin=(self.Tmin, "K"),
+            Tmax=(self.Tmax, "K"),
+            Pmin=(self.Pmin, "bar"),
+            Pmax=(self.Pmax, "bar"),
+            uncertainty=self.uncertainty,
+            comment=self.comment,
+        )
+
+    def test_is_identical_to(self):
+        """
+        Test that the KineticsModel.is_identical_to method works on itself.
+
+        This just checks the Temperature range
+        """
+        self.assertTrue(self.km.is_identical_to(self.km))
+
+        import copy
+        km = copy.deepcopy(self.km)
+        self.assertTrue(self.km.is_identical_to(self.km))
+
+        km.Tmax = (self.Tmax - 50, 'K')  # discrepancy must be more than 1%!
+        self.assertFalse(self.km.is_identical_to(km))
+
+    def test_repr(self):
+        """
+        Test that an KineticsModel object can be reconstructed from its repr()
+        output with no loss of information.
+        """
+        namespace = {}
+        exec('km = {0!r}'.format(self.km), globals(), namespace)
+        self.assertIn('km', namespace)
+        km = namespace['km']
+        self.assertTrue(self.km.is_identical_to(km))
+        self.assertEqual(dir(self.km), dir(km))
+        for att in 'Tmax Tmin Pmax Pmin comment uncertainty'.split():
+            self.assertEqual(repr(getattr(self.km, att)), repr(getattr(km, att)))
+
+    def test_pickle(self):
+        """
+        Test that an KineticsModel object can be pickled and unpickled
+        with no loss of information.
+        """
+        import pickle
+        km = pickle.loads(pickle.dumps(self.km, -1))
+        self.assertTrue(self.km.is_identical_to(km))
+        self.assertEqual(dir(self.km), dir(km))
+        for att in 'Tmax Tmin Pmax Pmin comment uncertainty'.split():
+            self.assertEqual(repr(getattr(self.km, att)), repr(getattr(km, att)))
+
 
 ################################################################################
 
@@ -44,7 +113,7 @@ class TestOrder(unittest.TestCase):
     Contains unit tests of the functions for converting rate coefficient units
     to/from reaction orders.
     """
-    
+
     def test_to_order_zeroth(self):
         """
         Test the conversion of zeroth-order rate coefficient units to an integer
@@ -54,14 +123,14 @@ def test_to_order_zeroth(self):
         self.assertEqual(0, get_reaction_order_from_rate_coefficient_units('mol/(cm^3*s)'))
         self.assertEqual(0, get_reaction_order_from_rate_coefficient_units('molecule/(m^3*s)'))
         self.assertEqual(0, get_reaction_order_from_rate_coefficient_units('molecule/(cm^3*s)'))
-        
+
     def test_to_order_first(self):
         """
         Test the conversion of first-order rate coefficient units to an integer
         reaction order.
         """
         self.assertEqual(1, get_reaction_order_from_rate_coefficient_units('s^-1'))
-        
+
     def test_to_order_second(self):
         """
         Test the conversion of second-order rate coefficient units to an integer
@@ -71,7 +140,7 @@ def test_to_order_second(self):
         self.assertEqual(2, get_reaction_order_from_rate_coefficient_units('cm^3/(mol*s)'))
         self.assertEqual(2, get_reaction_order_from_rate_coefficient_units('m^3/(molecule*s)'))
         self.assertEqual(2, get_reaction_order_from_rate_coefficient_units('cm^3/(molecule*s)'))
-        
+
     def test_to_order_third(self):
         """
         Test the conversion of third-order rate coefficient units to an integer
@@ -81,28 +150,28 @@ def test_to_order_third(self):
         self.assertEqual(3, get_reaction_order_from_rate_coefficient_units('cm^6/(mol^2*s)'))
         self.assertEqual(3, get_reaction_order_from_rate_coefficient_units('m^6/(molecule^2*s)'))
         self.assertEqual(3, get_reaction_order_from_rate_coefficient_units('cm^6/(molecule^2*s)'))
-        
+
     def test_to_units_zeroth(self):
         """
         Test the conversion of a reaction order of zero to rate coefficient
         units.
         """
         self.assertEqual('mol/(m^3*s)', get_rate_coefficient_units_from_reaction_order(0))
-        
+
     def test_to_units_first(self):
         """
         Test the conversion of a reaction order of one to rate coefficient
         units.
         """
         self.assertEqual('s^-1', get_rate_coefficient_units_from_reaction_order(1))
-        
+
     def test_to_units_second(self):
         """
         Test the conversion of a reaction order of two to rate coefficient
         units.
         """
         self.assertEqual('m^3/(mol*s)', get_rate_coefficient_units_from_reaction_order(2))
-        
+
     def test_to_units_third(self):
         """
         Test the conversion of a reaction order of three to rate coefficient

rmgpy/kinetics/surface.pyx

@@ -393,6 +393,7 @@ cdef class SurfaceArrhenius(Arrhenius):
     `Tmax`          The maximum temperature at which the model is valid, or zero if unknown or undefined
     `Pmin`          The minimum pressure at which the model is valid, or zero if unknown or undefined
     `Pmax`          The maximum pressure at which the model is valid, or zero if unknown or undefined
+    `uncertainty`   Uncertainty information
     `comment`       Information about the model (e.g. its source)
     =============== =============================================================
     """
@@ -415,6 +416,7 @@ cdef class SurfaceArrhenius(Arrhenius):
         if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax)
         if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin)
         if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax)
+        if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty)
         if self.comment != '': string += ', comment="""{0}"""'.format(self.comment)
         string += ')'
         return string
@@ -424,7 +426,7 @@ cdef class SurfaceArrhenius(Arrhenius):
         A helper function used when pickling a SurfaceArrhenius object.
         """
         return (SurfaceArrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax,
-                                   self.comment))
+                                   self.uncertainty, self.comment))
 
 
 ################################################################################
@@ -451,6 +453,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP):
     `Tmax`          The maximum temperature at which the model is valid, or zero if unknown or undefined
     `Pmin`          The minimum pressure at which the model is valid, or zero if unknown or undefined
     `Pmax`          The maximum pressure at which the model is valid, or zero if unknown or undefined
+    `uncertainty`   Uncertainty information
     `comment`       Information about the model (e.g. its source)
     =============== =============================================================
     
@@ -475,6 +478,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP):
         if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax)
         if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin)
         if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax)
+        if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty)
         if self.comment != '': string += ', comment="""{0}"""'.format(self.comment)
         string += ')'
         return string
@@ -484,7 +488,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP):
         A helper function used when pickling an SurfaceArrheniusBEP object.
         """
         return (SurfaceArrheniusBEP, (self.A, self.n, self.alpha, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax,
-                                      self.comment))
+                                      self.uncertainty, self.comment))
 
     cpdef SurfaceArrhenius to_arrhenius(self, double dHrxn):
         """
@@ -502,5 +506,6 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP):
             T0=(1, "K"),
             Tmin=self.Tmin,
             Tmax=self.Tmax,
+            uncertainty = self.uncertainty,
             comment=self.comment,
         )