<commit_before>#!/usr/bin/env python
# encoding: utf-8

################################################################################
#
#   RMG - Reaction Mechanism Generator
#
#   Copyright (c) 2009-2011 by the RMG Team (<EMAIL>)
#
#   Permission is hereby granted, free of charge, to any person obtaining a
#   copy of this software and associated documentation files (the 'Software'),
#   to deal in the Software without restriction, including without limitation
#   the rights to use, copy, modify, merge, publish, distribute, sublicense,
#   and/or sell copies of the Software, and to permit persons to whom the
#   Software is furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#   DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
This module contains classes and methods for working with kinetics models. All 
such models derive from the :class:`KineticsModel` base class:

* :class:`Arrhenius` - A kinetics model based on the modified Arrhenius
  equation

* :class:`ArrheniusEP` - A kinetics model based on the modified Arrhenius
  equation with Evans-Polanyi correction to the activation energy

* :class:`MultiKinetics` - A kinetics model based on a summation of
  several other kinetics expressions

* :class:`ThirdBody` - A pressure-dependent kinetics model based on the  
  modified Arrhenius equation, but with an additional factor for the third body 
  concentration

* :class:`Lindemann` - A pressure-dependent kinetics model based on the
  Lindemann equation

* :class:`Troe` - A pressure-dependent kinetics model based on the
  Lindemann equation with improved Troe falloff factor

* :class:`PDepArrhenius` - A pressure-dependent kinetics model based on a
  set of modified Arrhenius equations at various pressures, which are then
  interpolated between on a logarithmic pressure scale

* :class:`Chebyshev` - A pressure-dependent kinetics model using an array
  of Chebyshev polynomials in inverse temperature and logarithmic pressure

"""

################################################################################

import math
import numpy
import numpy.linalg
import cython
import re
import logging

from quantity import Quantity, constants
from molecule import Molecule

################################################################################

class KineticsError(Exception):
    """
    An exception to be raised when an error occurs while working with 
    kinetics models and data. Pass a string describing the circumstances of the
    exceptional behavior.
    """
    pass

################################################################################

class KineticsModel:
    """
    A base class for kinetics models, containing several attributes common to 
    all models:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `Tmin`          :class:`Quantity`   The minimum absolute temperature in K at which the model is valid
    `Tmax`          :class:`Quantity`   The maximum absolute temperature in K at which the model is valid
    `Pmin`          :class:`Quantity`   The minimum absolute pressure in Pa at which the model is valid
    `Pmax`          :class:`Quantity`   The maximum absolute pressure in Pa at which the model is valid
    `comment`       :class:`str`        A string containing information about the model (e.g. its source)
    =============== =================== ========================================

    """

    def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''):
        if Tmin is not None:
            self.Tmin = Quantity(Tmin)
        else:
            self.Tmin = None
        if Tmax is not None:
            self.Tmax = Quantity(Tmax)
        else:
            self.Tmax = None
        if Pmin is not None:
            self.Pmin = Quantity(Pmin)
        else:
            self.Pmin = None
        if Pmax is not None:
            self.Pmax = Quantity(Pmax)
        else:
            self.Pmax = None
        self.comment = comment

    def __repr__(self):
        """
        Return a string representation that can be used to reconstruct the
        KineticsModel object.
        """
        string = self.toPrettyRepr()
        string = re.sub(r'\(\n    ', '(', string)
        string = re.sub(r',\n    ', ', ', string)
        string = re.sub(r',\n\)', ')', string)
        string = re.sub(r' = ', '=', string)
        return string

    def toPrettyRepr(self):
        """
        Return a string representation that can be used to reconstruct the
        KineticsModel object.
        """
        raise NotImplementedError('You must implement this method in your derived class.')
    
    def __reduce__(self):
        """
        A helper function used when pickling a KineticsModel object.
        """
        return (KineticsModel, (self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.comment))

    def isTemperatureValid(self, T):
        """
        Return ``True`` if the temperature `T` in K is within the valid
        temperature range of the kinetic data, or ``False`` if not. If
        the minimum and maximum temperature are not defined, ``True`` is 
        returned.
        """
        return self.Tmin is None or self.Tmax is None or (self.Tmin.value <= T and T <= self.Tmax.value)

    def isPressureValid(self, P):
        """
        Return ``True`` if the pressure `P` in Pa is within the valid
        pressure range of the kinetic data, or ``False`` if not. If
        the minimum and maximum pressure are not defined, ``True`` is 
        returned.
        """
        return self.Pmin is None or self.Pmax is None or (self.Pmin.value <= P and P <= self.Pmax.value)

    def isPressureDependent(self):
        """
        Return ``True`` if the kinetics are pressure-dependent or ``False`` if
        they are pressure-independent. This method must be overloaded in the 
        derived class.
        """
        raise KineticsError('Unexpected call to KineticsModel.isPressureDependent(); you should be using a class derived from KineticsModel.')

    def getRateCoefficients(self, Tlist):
        """
        Return the rate coefficients k(T) in SI units at temperatures
        `Tlist` in K.
        """
        return numpy.array([self.getRateCoefficient(T) for T in Tlist], numpy.float64)
    
    def toHTML(self):
        """
        Return an HTML rendering.
        """    
        cython.declare(T=cython.double, k=cython.double)
        cython.declare(string=str)
        cython.declare(kdata1=list, Tdata=list, kdata10=list)
        string = '<table class="KineticsData"><tr class="KineticsData_Tdata"><th>T/[K]</th>\n   '
        kdata1 = []
        kdata10 = []
        Tdata = [500,1000,1500,2000]
        for T in Tdata:
            string += '<td>{0:.0f}</td>'.format(T)
            kdata1.append(self.getRateCoefficient(T,P=1e5))
            if self.isPressureDependent():
                kdata10.append(self.getRateCoefficient(T,P=1e6))
                
        if self.isPressureDependent():
            string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k(1 bar)/[mole,m,s])\n    '
            for k in kdata1:
                string += '<td>{0:+.1f}'.format(math.log10(k))
            string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k(10 bar)/[mole,m,s])\n    '
            for k in kdata10:
                string += '<td>{0:+.1f}'.format(math.log10(k))
        else:
            string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k/[mole,m,s])\n    '
            for k in kdata1:
                string += '<td>{0:+.1f}'.format(math.log10(k))
     
        string += '\n</tr></table>'
        string += "<span class='KineticsData_repr'>{0!r}</span>".format(self)
        return string

    def isSimilarTo(self,otherKinetics):
        """
        Returns ``True`` if rates of reaction at temperatures 500,1000,1500,2000 K
        and 1 and 10 bar are within +/ .3 for log(k), in other words, within a factor of 2.
        """
        cython.declare(T=cython.double)
        cython.declare(Tdata=list)
        Tdata = [500,1000,1500,2000]

        if self.isPressureDependent() and otherKinetics.isPressureDependent():
            for T in Tdata:
                if (abs(math.log10(self.getRateCoefficient(T,P=1e5))-math.log10(otherKinetics.getRateCoefficient(T,P=1e5))) <= .5
                and abs(math.log10(self.getRateCoefficient(T,P=1e6))-math.log10(otherKinetics.getRateCoefficient(T,P=1e6))) <= .5):
                    pass
                else:
                    return False
            return True

        elif not self.isPressureDependent() and not otherKinetics.isPressureDependent():
            for T in Tdata:
                if abs(math.log10(self.getRateCoefficient(T,P=1e5))-math.log10(otherKinetics.getRateCoefficient(T,P=1e5))) <= .5:
                    pass
                else:
                    return False
            return True

        else:
            return False


    def isIdenticalTo(self,otherKinetics):
        """
        Returns ``True`` if Tmin, Tmax, Pmin, Pmax for both objects match.
        Otherwise returns ``False``
        """
        if self.Tmin is not None and otherKinetics.Tmin is not None:
            if self.Tmin.equals(otherKinetics.Tmin):
                pass
            else:
                return False
        elif self.Tmin is None and otherKinetics.Tmin is None:
            pass
        else:
            return False

        if self.Tmax is not None and otherKinetics.Tmax is not None:
            if self.Tmax.equals(otherKinetics.Tmax):
                pass
            else:
                return False
        elif self.Tmax is None and otherKinetics.Tmax is None:
            pass
        else:
            return False

        if self.Pmin is not None and otherKinetics.Pmin is not None:
            if self.Pmin.equals(otherKinetics.Pmin):
                pass
            else:
                return False
        elif self.Pmin is None and otherKinetics.Pmin is None:
            pass
        else:
            return False

        if self.Pmax is not None and otherKinetics.Pmax is not None:
            if self.Pmax.equals(otherKinetics.Pmax):
                pass
            else:
                return False
        elif self.Pmax is None and otherKinetics.Pmax is None:
            pass
        else:
            return False

        return True

################################################################################

class KineticsData(KineticsModel):
    """
    A kinetics model based around a set of discrete (high-pressure limit)
    rate coefficients at various temperatures. The attributes are:

    =========== =================== ============================================
    Attribute   Type                Description
    =========== =================== ============================================
    `Tdata`     :class:`Quantity`   The temperatures at which the heat capacity data is provided
    `kdata`     :class:`Quantity`   The rate coefficients in SI units at each temperature in `Tdata`
    =========== =================== ============================================
    
    """

    def __init__(self, Tdata=None, kdata=None, Tmin=None, Tmax=None, comment=''):
        KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, comment=comment)
        self.Tdata = Quantity(Tdata)
        self.kdata = Quantity(kdata)

    def toPrettyRepr(self):
        """
        Return a string representation of the reference that can be used to
        reconstruct the object.
        """
        string = u'KineticsData(\n'
        string += u'    Tdata = {0!r},\n'.format(self.Tdata)
        string += u'    kdata = {0!r},\n'.format(self.kdata)
        if self.Tmin is not None: string += '    Tmin = {0!r},\n'.format(self.Tmin)
        if self.Tmax is not None: string += '    Tmax = {0!r},\n'.format(self.Tmax)
        if self.comment != '': string += '    comment = """{0}""",\n'.format(self.comment)
        return string + u')'
    
    def __reduce__(self):
        """
        A helper function used when pickling a KineticsData object.
        """
        return (KineticsData, (self.Tdata, self.kdata, self.Tmin, self.Tmax, self.comment))
    
    def toHTML(self):
        """
        Return an HTML rendering.
        """
        cython.declare(T=cython.double, k=cython.double)
        cython.declare(string=str)
        string = '<table class="KineticsData"><tr class="KineticsData_Tdata"><th>T/[{0!s}]</th>\n   '.format(self.Tdata.units)
        for T in self.Tdata.values:
            string += '<td>{0:.0f}</td>'.format(T)
        string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k/[{0!s}])\n    '.format(self.kdata.units)
        for k in self.kdata.values:
            string += '<td>{0:+.1f}'.format(math.log10(k))
        string += '\n</tr></table>'
        return string

    def isPressureDependent(self):
        """
        Returns ``False`` since KineticsData kinetics are not 
        pressure-dependent.
        """
        return False

    def getRateCoefficient(self, T, P=1e5):
        """
        Return the rate coefficient k(T) in SI units at temperature
        `T` in K.
        """
        cython.declare(Tmin=cython.double, Tmax=cython.double, kmin=cython.double, kmax=cython.double)
        cython.declare(k=cython.double, slope=cython.double)
        k = 0.0
        if not self.isTemperatureValid(T):
            raise KineticsError('Invalid temperature "{0:g} K" for heat capacity estimation.'.format(T))
        if T < numpy.min(self.Tdata.values):
            k = self.kdata.values[0]
        elif T >= numpy.max(self.Tdata.values):
            k = self.kdata.values[-1]
        else:
            for Tmin, Tmax, kmin, kmax in zip(self.Tdata.values[:-1], self.Tdata.values[1:], self.kdata.values[:-1], self.kdata.values[1:]):
                if Tmin <= T and T < Tmax:
                    slope = (1.0/T - 1.0/Tmin) / (1.0/Tmax - 1.0/Tmin)
                    k = kmin * (kmax / kmin)**slope
        return k
    
    def toArrhenius(self):
        """
        Return an :class:`Arrhenius` expression fitted to this data
        """
        return Arrhenius(comment=self.comment, Tmin=self.Tmin, Tmax=self.Tmax).fitToData(Tlist=self.Tdata.values, klist=self.kdata.values, kunits=self.kdata.units)

    def isIdenticalTo(self,otherKinetics):
        """
        Returns ``True`` if the kdata and Tdata match. Returns ``False`` otherwise.
        """
        if isinstance(otherKinetics,KineticsData):
            if KineticsModel.isIdenticalTo(self,otherKinetics):
                if self.Tdata.equals(otherKinetics.Tdata) and self.kdata.equals(otherKinetics.kdata):
                    return True

        return False
################################################################################

class Arrhenius(KineticsModel):
    """
    Represent a set of modified Arrhenius kinetics. The kinetic expression has
    the form

    .. math:: k(T) = A \\left( \\frac{T}{T_0} \\right)^n \\exp \\left( - \\frac{E_\\mathrm{a}}{RT} \\right)

    where :math:`A`, :math:`n`, :math:`E_\\mathrm{a}`, and :math:`T_0` are the
    parameters to be set, :math:`T` is absolute temperature, and :math:`R` is
    the gas law constant. The attributes are:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `A`             :class:`Quantity`   The preexponential factor in s^-1, m^3/mol*s, etc.
    `T0`            :class:`Quantity`   The reference temperature in K
    `n`             :class:`Quantity`   The temperature exponent
    `Ea`            :class:`Quantity`   The activation energy in J/mol
    =============== =================== ========================================
    
    """
    
    def __init__(self, A=0.0, n=0.0, Ea=0.0, T0=1.0, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''):
        KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, comment=comment)
        self.A = Quantity(A)
        self.T0 = Quantity(T0)
        self.n = Quantity(n)
        self.Ea = Quantity(Ea)
    
    def toPrettyRepr(self):
        """
        Return a string representation of the reference that can be used to
        reconstruct the object.
        """
        string = u'Arrhenius(\n'
        string += u'    A = {0!r},\n'.format(self.A)
        string += u'    n = {0!r},\n'.format(self.n)
        string += u'    Ea = {0!r},\n'.format(self.Ea)
        string += u'    T0 = {0!r},\n'.format(self.T0)
        if self.Tmin is not None: string += '    Tmin = {0!r},\n'.format(self.Tmin)
        if self.Tmax is not None: string += '    Tmax = {0!r},\n'.format(self.Tmax)
        if self.Pmin is not None: string += '    Pmin = {0!r},\n'.format(self.Pmin)
        if self.Pmax is not None: string += '    Pmax = {0!r},\n'.format(self.Pmax)
        if self.comment != '': string += '    comment = """{0}""",\n'.format(self.comment)
        return string + u')'
    
    def __str__(self):
        """
        Return a string representation that is a bit shorter and prettier than __repr__.
        """
        string = 'Arrhenius(A={0!r}, n={1!r}, Ea={2!r}, T0={3!r})'.format(self.A, self.n, self.Ea, self.T0)
        return string

    def __reduce__(self):
        """
        A helper function used when pickling an Arrhenius object.
        """
        return (Arrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.comment))

    def isPressureDependent(self):
        """
        Returns ``False`` since Arrhenius kinetics are not pressure-dependent.
        """
        return False

    def getRateCoefficient(self, T, P=1e5):
        """
        Return the rate coefficient k(T) in SI units at temperature 
        `T` in K.
        """
        return self.A.value * (T / self.T0.value)** self.n.value * math.exp(-self.Ea.value / constants.R / T)

    def changeT0(self, T0):
        """
        Changes the reference temperature used in the exponent to `T0` in K, and
        adjusts the preexponential accordingly.
        """
        self.A.value /= (self.T0.value / T0)**self.n.value
        self.T0.value = T0

    def fitToData(self, Tlist, klist, kunits, T0=1, weights=None):
        """
        Fit the Arrhenius parameters to a set of rate coefficient data `klist`
        in units of `kunits` corresponding to a set of temperatures `Tlist` in 
        K. A linear least-squares fit is used, which guarantees that the 
        resulting parameters provide the best possible approximation to the 
        data.
        """
        import numpy.linalg
        import scipy.stats
        A = numpy.zeros((len(Tlist),3), numpy.float64)
        A[:,0] = numpy.ones_like(Tlist)
        A[:,1] = numpy.log(Tlist / T0)
        A[:,2] = -1.0 / constants.R / Tlist
        b = numpy.log(klist)
        if weights is not None:
            for n in range(b.size):
                A[n,:] *= weights[n]
                b[n] *= weights[n]
        x, residues, rank, s = numpy.linalg.lstsq(A,b)
        
        # Determine covarianace matrix to obtain parameter uncertainties
        count = klist.size
        cov = residues[0] / (count - 3) * numpy.linalg.inv(numpy.dot(A.T, A))
        t = scipy.stats.t.ppf(0.975, count - 3)
            
        self.A = Quantity(math.exp(x[0]), kunits, '*|/', t * math.exp(math.sqrt(cov[0,0])))
        self.n = Quantity(x[1], '', '+|-', t * math.sqrt(cov[1,1]))
        self.Ea = Quantity(x[2], "J/mol", '+|-', t * math.sqrt(cov[2,2]))
        self.T0 = Quantity(T0, "K")
        return self

    def isIdenticalTo(self, otherKinetics):
        """
        Returns ``True`` if kinetics matches that of another kinetics model.  Must match temperature
        and pressure range of kinetics model, as well as parameters: A, n, Ea, T0. (Shouldn't have pressure
        range if it's Arrhenius.) Otherwise returns ``False``.
        """
        if isinstance(otherKinetics,Arrhenius):
            if KineticsModel.isIdenticalTo(self, otherKinetics):
                if (self.A.equals(otherKinetics.A) and self.n.equals(otherKinetics.n)
                and self.Ea.equals(otherKinetics.Ea) and self.T0.equals(otherKinetics.T0)):
                    return True
                
        return False

    
################################################################################

class ArrheniusEP(KineticsModel):
    """
    Represent a set of modified Arrhenius kinetics with Evans-Polanyi data. The
    kinetic expression has the form

    .. math:: k(T) = A T^n \\exp \\left( - \\frac{E_0 + \\alpha \\Delta H_\\mathrm{rxn}}{RT} \\right)
<commit_msg>Fixed residual handling bug in Arrhenius.fitToData<commit_after>#!/usr/bin/env python
# encoding: utf-8

################################################################################
#
#   RMG - Reaction Mechanism Generator
#
#   Copyright (c) 2009-2011 by the RMG Team (<EMAIL>)
#
#   Permission is hereby granted, free of charge, to any person obtaining a
#   copy of this software and associated documentation files (the 'Software'),
#   to deal in the Software without restriction, including without limitation
#   the rights to use, copy, modify, merge, publish, distribute, sublicense,
#   and/or sell copies of the Software, and to permit persons to whom the
#   Software is furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#   DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
This module contains classes and methods for working with kinetics models. All 
such models derive from the :class:`KineticsModel` base class:

* :class:`Arrhenius` - A kinetics model based on the modified Arrhenius
  equation

* :class:`ArrheniusEP` - A kinetics model based on the modified Arrhenius
  equation with Evans-Polanyi correction to the activation energy

* :class:`MultiKinetics` - A kinetics model based on a summation of
  several other kinetics expressions

* :class:`ThirdBody` - A pressure-dependent kinetics model based on the  
  modified Arrhenius equation, but with an additional factor for the third body 
  concentration

* :class:`Lindemann` - A pressure-dependent kinetics model based on the
  Lindemann equation

* :class:`Troe` - A pressure-dependent kinetics model based on the
  Lindemann equation with improved Troe falloff factor

* :class:`PDepArrhenius` - A pressure-dependent kinetics model based on a
  set of modified Arrhenius equations at various pressures, which are then
  interpolated between on a logarithmic pressure scale

* :class:`Chebyshev` - A pressure-dependent kinetics model using an array
  of Chebyshev polynomials in inverse temperature and logarithmic pressure

"""

################################################################################

import math
import numpy
import numpy.linalg
import cython
import re
import logging

from quantity import Quantity, constants
from molecule import Molecule

################################################################################

class KineticsError(Exception):
    """
    An exception to be raised when an error occurs while working with 
    kinetics models and data. Pass a string describing the circumstances of the
    exceptional behavior.
    """
    pass

################################################################################

class KineticsModel:
    """
    A base class for kinetics models, containing several attributes common to 
    all models:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `Tmin`          :class:`Quantity`   The minimum absolute temperature in K at which the model is valid
    `Tmax`          :class:`Quantity`   The maximum absolute temperature in K at which the model is valid
    `Pmin`          :class:`Quantity`   The minimum absolute pressure in Pa at which the model is valid
    `Pmax`          :class:`Quantity`   The maximum absolute pressure in Pa at which the model is valid
    `comment`       :class:`str`        A string containing information about the model (e.g. its source)
    =============== =================== ========================================

    """

    def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''):
        if Tmin is not None:
            self.Tmin = Quantity(Tmin)
        else:
            self.Tmin = None
        if Tmax is not None:
            self.Tmax = Quantity(Tmax)
        else:
            self.Tmax = None
        if Pmin is not None:
            self.Pmin = Quantity(Pmin)
        else:
            self.Pmin = None
        if Pmax is not None:
            self.Pmax = Quantity(Pmax)
        else:
            self.Pmax = None
        self.comment = comment

    def __repr__(self):
        """
        Return a string representation that can be used to reconstruct the
        KineticsModel object.
        """
        string = self.toPrettyRepr()
        string = re.sub(r'\(\n    ', '(', string)
        string = re.sub(r',\n    ', ', ', string)
        string = re.sub(r',\n\)', ')', string)
        string = re.sub(r' = ', '=', string)
        return string

    def toPrettyRepr(self):
        """
        Return a string representation that can be used to reconstruct the
        KineticsModel object.
        """
        raise NotImplementedError('You must implement this method in your derived class.')
    
    def __reduce__(self):
        """
        A helper function used when pickling a KineticsModel object.
        """
        return (KineticsModel, (self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.comment))

    def isTemperatureValid(self, T):
        """
        Return ``True`` if the temperature `T` in K is within the valid
        temperature range of the kinetic data, or ``False`` if not. If
        the minimum and maximum temperature are not defined, ``True`` is 
        returned.
        """
        return self.Tmin is None or self.Tmax is None or (self.Tmin.value <= T and T <= self.Tmax.value)

    def isPressureValid(self, P):
        """
        Return ``True`` if the pressure `P` in Pa is within the valid
        pressure range of the kinetic data, or ``False`` if not. If
        the minimum and maximum pressure are not defined, ``True`` is 
        returned.
        """
        return self.Pmin is None or self.Pmax is None or (self.Pmin.value <= P and P <= self.Pmax.value)

    def isPressureDependent(self):
        """
        Return ``True`` if the kinetics are pressure-dependent or ``False`` if
        they are pressure-independent. This method must be overloaded in the 
        derived class.
        """
        raise KineticsError('Unexpected call to KineticsModel.isPressureDependent(); you should be using a class derived from KineticsModel.')

    def getRateCoefficients(self, Tlist):
        """
        Return the rate coefficients k(T) in SI units at temperatures
        `Tlist` in K.
        """
        return numpy.array([self.getRateCoefficient(T) for T in Tlist], numpy.float64)
    
    def toHTML(self):
        """
        Return an HTML rendering.
        """    
        cython.declare(T=cython.double, k=cython.double)
        cython.declare(string=str)
        cython.declare(kdata1=list, Tdata=list, kdata10=list)
        string = '<table class="KineticsData"><tr class="KineticsData_Tdata"><th>T/[K]</th>\n   '
        kdata1 = []
        kdata10 = []
        Tdata = [500,1000,1500,2000]
        for T in Tdata:
            string += '<td>{0:.0f}</td>'.format(T)
            kdata1.append(self.getRateCoefficient(T,P=1e5))
            if self.isPressureDependent():
                kdata10.append(self.getRateCoefficient(T,P=1e6))
                
        if self.isPressureDependent():
            string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k(1 bar)/[mole,m,s])\n    '
            for k in kdata1:
                string += '<td>{0:+.1f}'.format(math.log10(k))
            string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k(10 bar)/[mole,m,s])\n    '
            for k in kdata10:
                string += '<td>{0:+.1f}'.format(math.log10(k))
        else:
            string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k/[mole,m,s])\n    '
            for k in kdata1:
                string += '<td>{0:+.1f}'.format(math.log10(k))
     
        string += '\n</tr></table>'
        string += "<span class='KineticsData_repr'>{0!r}</span>".format(self)
        return string

    def isSimilarTo(self,otherKinetics):
        """
        Returns ``True`` if rates of reaction at temperatures 500,1000,1500,2000 K
        and 1 and 10 bar are within +/ .3 for log(k), in other words, within a factor of 2.
        """
        cython.declare(T=cython.double)
        cython.declare(Tdata=list)
        Tdata = [500,1000,1500,2000]

        if self.isPressureDependent() and otherKinetics.isPressureDependent():
            for T in Tdata:
                if (abs(math.log10(self.getRateCoefficient(T,P=1e5))-math.log10(otherKinetics.getRateCoefficient(T,P=1e5))) <= .5
                and abs(math.log10(self.getRateCoefficient(T,P=1e6))-math.log10(otherKinetics.getRateCoefficient(T,P=1e6))) <= .5):
                    pass
                else:
                    return False
            return True

        elif not self.isPressureDependent() and not otherKinetics.isPressureDependent():
            for T in Tdata:
                if abs(math.log10(self.getRateCoefficient(T,P=1e5))-math.log10(otherKinetics.getRateCoefficient(T,P=1e5))) <= .5:
                    pass
                else:
                    return False
            return True

        else:
            return False


    def isIdenticalTo(self,otherKinetics):
        """
        Returns ``True`` if Tmin, Tmax, Pmin, Pmax for both objects match.
        Otherwise returns ``False``
        """
        if self.Tmin is not None and otherKinetics.Tmin is not None:
            if self.Tmin.equals(otherKinetics.Tmin):
                pass
            else:
                return False
        elif self.Tmin is None and otherKinetics.Tmin is None:
            pass
        else:
            return False

        if self.Tmax is not None and otherKinetics.Tmax is not None:
            if self.Tmax.equals(otherKinetics.Tmax):
                pass
            else:
                return False
        elif self.Tmax is None and otherKinetics.Tmax is None:
            pass
        else:
            return False

        if self.Pmin is not None and otherKinetics.Pmin is not None:
            if self.Pmin.equals(otherKinetics.Pmin):
                pass
            else:
                return False
        elif self.Pmin is None and otherKinetics.Pmin is None:
            pass
        else:
            return False

        if self.Pmax is not None and otherKinetics.Pmax is not None:
            if self.Pmax.equals(otherKinetics.Pmax):
                pass
            else:
                return False
        elif self.Pmax is None and otherKinetics.Pmax is None:
            pass
        else:
            return False

        return True

################################################################################

class KineticsData(KineticsModel):
    """
    A kinetics model based around a set of discrete (high-pressure limit)
    rate coefficients at various temperatures. The attributes are:

    =========== =================== ============================================
    Attribute   Type                Description
    =========== =================== ============================================
    `Tdata`     :class:`Quantity`   The temperatures at which the heat capacity data is provided
    `kdata`     :class:`Quantity`   The rate coefficients in SI units at each temperature in `Tdata`
    =========== =================== ============================================
    
    """

    def __init__(self, Tdata=None, kdata=None, Tmin=None, Tmax=None, comment=''):
        KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, comment=comment)
        self.Tdata = Quantity(Tdata)
        self.kdata = Quantity(kdata)

    def toPrettyRepr(self):
        """
        Return a string representation of the reference that can be used to
        reconstruct the object.
        """
        string = u'KineticsData(\n'
        string += u'    Tdata = {0!r},\n'.format(self.Tdata)
        string += u'    kdata = {0!r},\n'.format(self.kdata)
        if self.Tmin is not None: string += '    Tmin = {0!r},\n'.format(self.Tmin)
        if self.Tmax is not None: string += '    Tmax = {0!r},\n'.format(self.Tmax)
        if self.comment != '': string += '    comment = """{0}""",\n'.format(self.comment)
        return string + u')'
    
    def __reduce__(self):
        """
        A helper function used when pickling a KineticsData object.
        """
        return (KineticsData, (self.Tdata, self.kdata, self.Tmin, self.Tmax, self.comment))
    
    def toHTML(self):
        """
        Return an HTML rendering.
        """
        cython.declare(T=cython.double, k=cython.double)
        cython.declare(string=str)
        string = '<table class="KineticsData"><tr class="KineticsData_Tdata"><th>T/[{0!s}]</th>\n   '.format(self.Tdata.units)
        for T in self.Tdata.values:
            string += '<td>{0:.0f}</td>'.format(T)
        string += '\n</tr><tr class="KineticsData_kdata"><th>log<sub>10</sub>(k/[{0!s}])\n    '.format(self.kdata.units)
        for k in self.kdata.values:
            string += '<td>{0:+.1f}'.format(math.log10(k))
        string += '\n</tr></table>'
        return string

    def isPressureDependent(self):
        """
        Returns ``False`` since KineticsData kinetics are not 
        pressure-dependent.
        """
        return False

    def getRateCoefficient(self, T, P=1e5):
        """
        Return the rate coefficient k(T) in SI units at temperature
        `T` in K.
        """
        cython.declare(Tmin=cython.double, Tmax=cython.double, kmin=cython.double, kmax=cython.double)
        cython.declare(k=cython.double, slope=cython.double)
        k = 0.0
        if not self.isTemperatureValid(T):
            raise KineticsError('Invalid temperature "{0:g} K" for heat capacity estimation.'.format(T))
        if T < numpy.min(self.Tdata.values):
            k = self.kdata.values[0]
        elif T >= numpy.max(self.Tdata.values):
            k = self.kdata.values[-1]
        else:
            for Tmin, Tmax, kmin, kmax in zip(self.Tdata.values[:-1], self.Tdata.values[1:], self.kdata.values[:-1], self.kdata.values[1:]):
                if Tmin <= T and T < Tmax:
                    slope = (1.0/T - 1.0/Tmin) / (1.0/Tmax - 1.0/Tmin)
                    k = kmin * (kmax / kmin)**slope
        return k
    
    def toArrhenius(self):
        """
        Return an :class:`Arrhenius` expression fitted to this data
        """
        return Arrhenius(comment=self.comment, Tmin=self.Tmin, Tmax=self.Tmax).fitToData(Tlist=self.Tdata.values, klist=self.kdata.values, kunits=self.kdata.units)

    def isIdenticalTo(self,otherKinetics):
        """
        Returns ``True`` if the kdata and Tdata match. Returns ``False`` otherwise.
        """
        if isinstance(otherKinetics,KineticsData):
            if KineticsModel.isIdenticalTo(self,otherKinetics):
                if self.Tdata.equals(otherKinetics.Tdata) and self.kdata.equals(otherKinetics.kdata):
                    return True

        return False
################################################################################

class Arrhenius(KineticsModel):
    """
    Represent a set of modified Arrhenius kinetics. The kinetic expression has
    the form

    .. math:: k(T) = A \\left( \\frac{T}{T_0} \\right)^n \\exp \\left( - \\frac{E_\\mathrm{a}}{RT} \\right)

    where :math:`A`, :math:`n`, :math:`E_\\mathrm{a}`, and :math:`T_0` are the
    parameters to be set, :math:`T` is absolute temperature, and :math:`R` is
    the gas law constant. The attributes are:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `A`             :class:`Quantity`   The preexponential factor in s^-1, m^3/mol*s, etc.
    `T0`            :class:`Quantity`   The reference temperature in K
    `n`             :class:`Quantity`   The temperature exponent
    `Ea`            :class:`Quantity`   The activation energy in J/mol
    =============== =================== ========================================
    
    """
    
    def __init__(self, A=0.0, n=0.0, Ea=0.0, T0=1.0, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''):
        KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, comment=comment)
        self.A = Quantity(A)
        self.T0 = Quantity(T0)
        self.n = Quantity(n)
        self.Ea = Quantity(Ea)
    
    def toPrettyRepr(self):
        """
        Return a string representation of the reference that can be used to
        reconstruct the object.
        """
        string = u'Arrhenius(\n'
        string += u'    A = {0!r},\n'.format(self.A)
        string += u'    n = {0!r},\n'.format(self.n)
        string += u'    Ea = {0!r},\n'.format(self.Ea)
        string += u'    T0 = {0!r},\n'.format(self.T0)
        if self.Tmin is not None: string += '    Tmin = {0!r},\n'.format(self.Tmin)
        if self.Tmax is not None: string += '    Tmax = {0!r},\n'.format(self.Tmax)
        if self.Pmin is not None: string += '    Pmin = {0!r},\n'.format(self.Pmin)
        if self.Pmax is not None: string += '    Pmax = {0!r},\n'.format(self.Pmax)
        if self.comment != '': string += '    comment = """{0}""",\n'.format(self.comment)
        return string + u')'
    
    def __str__(self):
        """
        Return a string representation that is a bit shorter and prettier than __repr__.
        """
        string = 'Arrhenius(A={0!r}, n={1!r}, Ea={2!r}, T0={3!r})'.format(self.A, self.n, self.Ea, self.T0)
        return string

    def __reduce__(self):
        """
        A helper function used when pickling an Arrhenius object.
        """
        return (Arrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, self.comment))

    def isPressureDependent(self):
        """
        Returns ``False`` since Arrhenius kinetics are not pressure-dependent.
        """
        return False

    def getRateCoefficient(self, T, P=1e5):
        """
        Return the rate coefficient k(T) in SI units at temperature 
        `T` in K.
        """
        return self.A.value * (T / self.T0.value)** self.n.value * math.exp(-self.Ea.value / constants.R / T)

    def changeT0(self, T0):
        """
        Changes the reference temperature used in the exponent to `T0` in K, and
        adjusts the preexponential accordingly.
        """
        self.A.value /= (self.T0.value / T0)**self.n.value
        self.T0.value = T0

    def fitToData(self, Tlist, klist, kunits, T0=1, weights=None):
        """
        Fit the Arrhenius parameters to a set of rate coefficient data `klist`
        in units of `kunits` corresponding to a set of temperatures `Tlist` in 
        K. A linear least-squares fit is used, which guarantees that the 
        resulting parameters provide the best possible approximation to the 
        data.
        """
        import numpy.linalg
        import scipy.stats
        A = numpy.zeros((len(Tlist),3), numpy.float64)
        A[:,0] = numpy.ones_like(Tlist)
        A[:,1] = numpy.log(Tlist / T0)
        A[:,2] = -1.0 / constants.R / Tlist
        b = numpy.log(klist)
        if weights is not None:
            for n in range(b.size):
                A[n,:] *= weights[n]
                b[n] *= weights[n]
        x, residues, rank, s = numpy.linalg.lstsq(A,b)
        
        if residues:
            # Determine covarianace matrix to obtain parameter uncertainties
            count = klist.size
            t = scipy.stats.t.ppf(0.975, count - 3)
            cov = residues[0] / (count - 3) * numpy.linalg.inv(numpy.dot(A.T, A))
            
            self.A = Quantity(math.exp(x[0]), kunits, '*|/', t * math.exp(math.sqrt(cov[0,0])))
            self.n = Quantity(x[1], '', '+|-', t * math.sqrt(cov[1,1]))
            self.Ea = Quantity(x[2], "J/mol", '+|-', t * math.sqrt(cov[2,2]))
        else:
            # Ignore uncertainty if fitting to three or less points
            self.A = Quantity(math.exp(x[0]), kunits)
            self.n = Quantity(x[1], '')
            self.Ea = Quantity(x[2], "J/mol")
        self.T0 = Quantity(T0, "K")
        return self

    def isIdenticalTo(self, otherKinetics):
        """
        Returns ``True`` if kinetics matches that of another kinetics model.  Must match temperature
        and pressure range of kinetics model, as well as parameters: A, n, Ea, T0. (Shouldn't have pressure
        range if it's Arrhenius.) Otherwise returns ``False``.
        """
        if isinstance(otherKinetics,Arrhenius):
            if KineticsModel.isIdenticalTo(self, otherKinetics):
                if (self.A.equals(otherKinetics.A) and self.n.equals(otherKinetics.n)
                and self.Ea.equals(otherKinetics.Ea) and self.T0.equals(otherKinetics.T0)):
                    return True
                
        return False

    
################################################################################

class ArrheniusEP(KineticsModel):
    """
    Represent a set of modified Arrhenius kinetics with Evans-Polanyi data. The
    kinetic expression has the form

    .. math:: k(T) = A T^n \\exp \\left( - \\frac{E_0 + \\alpha \\Delta H_\\mathrm{rxn}}{RT} \\right)

{'fraction_non_alphanumeric': 0.1685375166766343, 'fraction_numerical': 0.011110088788701292, 'mean_word_length': 2.3898635477582846, 'pattern_counts': {'":': 0, '<': 109, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5350020', 'n_tokens_mistral': 13242, 'n_tokens_neox': 12415, 'n_words': 3669}